% SLB Project 

% Program to show the direct comparison between 
% the transverse acceleration and the vertical rotation rate

% 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


cd ..; % Shift to data directory


% Load the theoretical back-azimuth
baz_5sour=load('Baz_z0_theo.txt');
b_azo=baz_5sour(N_Rec,Res);
baz=b_azo*pi/180;


numvars=6; % Number of variables;


% Load the data of the simulated seismograms
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;
    At=An*sin(baz)-Ae*cos(baz);

Rz_norm=Rz/max(abs(Rz));
At_norm=At/max(abs(At));
    
% vertical acceleration
Az=diff(Vz)/dts;Az(length(Az)+1)=0;

% Apparent S Velocity 
app_Vel=max(abs(At))/max(abs(Rz))/2;app_Vel=round(app_Vel);


figure(1)
subplot(411),plot(Ttime,Az,'k-'), ylabel(' A_Z (m/s^2)'),axis([-Inf, Inf, -max(abs(Az)), max(abs(Az))]);
title(sprintf('Apparent S velocity %g (m/s)', app_Vel));
subplot(412),plot(Ttime,At,'b-'), ylabel(' A_T (m/s^2) '),axis([-Inf, Inf, -max(abs(At)), max(abs(At))]);
subplot(413),plot(Ttime,Rz,'r-'), ylabel('R_Z (rad/s)'),axis([-Inf, Inf, -max(abs(Rz)), max(abs(Rz))]);
subplot(414),plot(Ttime,-Rz_norm,'r-',Ttime,At_norm,'b-'), axis([-Inf, Inf, -1, 1]); 
legend('Normalized R_Z','Normalized A_T','Location','SouthWest');
%legend('Normalized R_Z','Normalized A_T');
xlabel(' Time (s) ')
drawnow


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


cd MatLabCodes; % Back to program directory

