% SLB Project 

% Program to show the direct comparison between 
% the translational accelerations and the rotation rates

% Pham D. N., H. Igel, and M. Kaeser, 
% Geophysics, LMU Munich.
% Jan. 2009

clear all;
close all;

mtype=input('Enter the type of the medium: 1 or 2 (1- Homo, 2 - Layer)):');
        if mtype==1
            f_name1='Homo';
            Res=input('Enter the source index (1 to 5):');
        else
            f_name1='Layer';
            Res=input('Enter the source index (1 or 5):');
        end

        f_name2=[num2str(Res)];

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Scale factor for the source so that in all case Mw=0; 
        if (Res==1 & mtype==1)
            Amp_coef=8.4074865549767914e+005; 
        else
            Amp_coef=1;
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        
N_Rec=input('Enter the receiver index (1 to 15):');
f_name3=[num2str(N_Rec)];

        if Res<10
            filename=[f_name1 '0' f_name2 ' (' f_name3 ').dat'];
        else
            filename=[f_name1 f_name2 ' (' f_name3 ').dat'];
        end

% 1) Load the data of the simulated seismograms

cd ..; % Shift to data directory



numvars=6; % Number of variables;

fid = fopen(filename,'r');
tmp=fgets(fid);
tmp=fgets(fid);
tmp=fgets(fid);
tmp=fgets(fid);
tmp=fgets(fid);
S=fscanf(fid,'%f');
fclose('all');

for j=1:length(S)/(numvars+1)
    for i=1:numvars+1
         SEIS(i,j)=S(i+(j-1)*(numvars+1));
    end
end

SEIS(2:i,:)=SEIS(2:i,:)*Amp_coef;
   
Ttime=SEIS(1,:);
Vn=SEIS(2,:);
Ve=SEIS(3,:);
Vz=SEIS(4,:);
Rn=SEIS(5,:);
Re=SEIS(6,:);
Rz=SEIS(7,:);

dts=Ttime(2)-Ttime(1);
nts=length(Ttime);

    Ae= diff(Ve)/dts;Ae(length(Ae)+1)=0;
    An= diff(Vn)/dts;An(length(An)+1)=0;
    Az= diff(Vz)/dts;Az(length(Az)+1)=0;


Rn_norm=Rn/max(abs(Rn));    
Re_norm=Re/max(abs(Re));
Rz_norm=Rz/max(abs(Rz));
An_norm=An/max(abs(An));    
Ae_norm=Ae/max(abs(Ae));
Az_norm=Az/max(abs(Az));


figure(1)
subplot(421),plot(Ttime,An,'k-'), ylabel(' A_N (m/s^2)'),axis([-Inf, Inf, -max(abs(An)), max(abs(An))]);
subplot(423),plot(Ttime,Ae,'r-'), ylabel(' A_E (m/s^2)'),axis([-Inf, Inf, -max(abs(Ae)), max(abs(Ae))]);
subplot(425),plot(Ttime,Az,'g-'), ylabel(' A_Z (m/s^2)'),axis([-Inf, Inf, -max(abs(Az)), max(abs(Az))]);
xlabel(' Time (s) ');

subplot(422),plot(Ttime,Rn,'b-'), ylabel('R_N (rad/s)'),axis([-Inf, Inf, -max(abs(Rn)), max(abs(Rn))]);
subplot(424),plot(Ttime,Re,'b-'), ylabel('R_E (rad/s)'),axis([-Inf, Inf, -max(abs(Re)), max(abs(Re))]);
subplot(426),plot(Ttime,Rz,'b-'), ylabel('R_Z (rad/s)'),axis([-Inf, Inf, -max(abs(Rz)), max(abs(Rz))]);
xlabel(' Time (s) ');

subplot(4,2,7:8),plot(Ttime,-Rz_norm,'b-',Ttime,Az_norm,'g-',Ttime,Ae_norm,'r-'), axis([-Inf, Inf, -1, 1]); 
h=legend('Normalized R_Z','Normalized A_Z','Normalized A_E','Location','SouthWest');set(h,'FontSize',8);

xlabel(' Time (s) ')
drawnow


%disp(' Plotting to file ... ')
%pfile=sprintf('./Figs/Fig3_AR_%i_%i_%i.eps',mtype,Res,N_Rec); eval([' print -depsc2 ' pfile ]);
%pfile=sprintf('./Figs/Fig3_AR_%i_%i_%i.jpg',mtype,Res,N_Rec); eval([' print -djpeg100 ' pfile ]);


cd MatLabCodes; % Back to program directory

