% SLB Project 
% Program to estimate back azimuth 
% using polarizations of trans. and rot. motions of S waves
% the error of estimated back azimuth is then quantified

% TEST FOR ALL SOURCES AND RECEIVERS 

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


clear all;
close all;


cd ..; % Shift to data directory


% The theoretical back azimuth data
baz_5sour=load('Baz_z0_theo.txt');
    


s_ind=1:7; 
    % s_ind= 1 - 5 : Homogeneous medium
    % s_ind= 6 - 7 : Layered structure
    
for sss=s_ind % Loop for different sources & media

    if sss < 6
        f_name1='Homo';
        Res=sss; % Source index
        f_name2=[num2str(Res)];
    else
        f_name1='Layer';
        if sss==6
            Res=1; % Source index
        else
            Res=5; % Source index
        end
        f_name2=[num2str(Res)];   
    end


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Scale factor for the source so that in all case Mw=0; 
        if (sss==1)
            Amp_coef=8.4074865549767914e+005; 
        else
            Amp_coef=1;
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
    for rrr=1:15  % Loop for different receivers
    
        N_Rec=rrr; % Receiver index
        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


        numvars=6; % Number of components;

        % 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(rrr,:)=SEIS(1,:);
        Vn(rrr,:)=SEIS(2,:);
        Ve(rrr,:)=SEIS(3,:);
        Vz(rrr,:)=SEIS(4,:);

        Rn(rrr,:)=SEIS(5,:);
        Re(rrr,:)=SEIS(6,:);
        Rz(rrr,:)=SEIS(7,:);
    
        Er(rrr,:)=(Rn(rrr,:)).^2+(Re(rrr,:)).^2+(Rz(rrr,:)).^2;

        dts(rrr)=Ttime(rrr,2)-Ttime(rrr,1);

        Vnn=SEIS(2,:);
        Vee=SEIS(3,:);
        Vzz=SEIS(4,:);

        Aee= diff(Vee)/dts(rrr);Aee(length(Aee)+1)=0;Ae(rrr,:)=Aee;
        Ann= diff(Vnn)/dts(rrr);Ann(length(Ann)+1)=0;An(rrr,:)=Ann;
        Azz= diff(Vzz)/dts(rrr);Azz(length(Azz)+1)=0;Az(rrr,:)=Azz;

    

        Er_norm=Er(rrr,:)./max(abs(Er(rrr,:)));

        [Y1,i1]=max(abs(Er_norm));
        aaa=length(Er_norm);
        Er_norm([1:i1-1 i1+1:aaa])=0;%Er_norm(i1+5:aaa)=0;


        % Normalized acceleration
        An_norm=An(rrr,:).*Er_norm;
        Ae_norm=Ae(rrr,:).*Er_norm;
        Az_norm=Az(rrr,:).*Er_norm;
        Amax=max([max(abs(An_norm)) max(abs(Ae_norm)) max(abs(Az_norm))]);
        An_norm=An_norm./Amax;
        Ae_norm=Ae_norm./Amax;
        Az_norm=Az_norm./Amax;

        % Normalized rotation rate
        Rn_norm=Rn(rrr,:).*Er_norm;
        Re_norm=Re(rrr,:).*Er_norm;
        Rz_norm=Rz(rrr,:).*Er_norm;
        Rmax=max([max(abs(Rn_norm)) max(abs(Re_norm)) max(abs(Rz_norm))]);
        Rn_norm=Rn_norm./Rmax;
        Re_norm=Re_norm./Rmax;
        Rz_norm=Rz_norm./Rmax;

        % Cross product of polarizations of rotational motions and translational motion  
        for k=1:length(An_norm)
            Rdir(k,:)=cross([Rn_norm(k) Re_norm(k) Rz_norm(k)],...
                        [An_norm(k) Ae_norm(k) Az_norm(k)]);
        end
                    

        [Y2,i2]=max(abs(Rdir(:,1)));

 
        % Back azimuth estimation
        Est_Rdir=-1.*[Rdir(i2,1) Rdir(i2,2) 0]; % Use two horizontal components only
        Northing=[1 0 0]; % Northing is coincided with x axis
        C_ang=dot(Est_Rdir,Northing)/(norm(Est_Rdir)*norm(Northing)); % Using the definition of dot product 
        ax=cross(Est_Rdir,Northing);
        if ax(3)<0
            baz(sss,rrr)=acos(C_ang)/pi*180;
        else
            baz(sss,rrr)=360-acos(C_ang)/pi*180;
        end

        % Theoretical back azimuth        
        b_azo=baz_5sour(rrr,Res);
        % Back azimuth error
        baz_Err(sss,rrr)=baz(sss,rrr)-b_azo;


    
    end % End of the loop for different receivers


end % End of the loop for different sources & media



cl_pl=['b*';...
    'bs';...
    'bx';...
    'bp';...
    'bv'
    'r*';...
    'rv'];

figure(1)
Rec=1:15;
for pf=1:length(s_ind)
    plot(Rec,baz_Err(pf,:),cl_pl(pf,:));hold on;
end
xlabel('Receiver');
ylabel('Backazimuth error (deg)');
legend('Homo - Source 1','Homo - Source 2','Homo - Source 3','Homo - Source 4','Homo - Source 5',...
    'Layer - Source 1','Layer - Source 5','Location','NorthWest');

plot([0:16],[0:16].*0,'k-');
axis([0 16 -45 45]);
%pfile=sprintf('./Figs/Fig9_Baz_Err_crp.jpg'); eval([' print -djpeg100 ' pfile]);


cd MatLabCodes; % Back to program directory
