% SLB Project 
% Program to estimate propagation direction 
% (back azimuth & incidence angle)
% using polarizations of trans. and rot. motions of S waves

% Apply for one receiver, one source only 
% SLIDING TIME WINDOWS ARE USED (Solution 2: for each time window all samples are considered)

% Pham D. N., H. Igel, and M. Kaeser, 
% Geophysics, LMU Munich.
% Arp. 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';
            d=(1/15); %Sliding window length for the Homo. case
            Res=input('Enter the source index (1 to 5):');
        else
            f_name1='Layer';
            d=(1/30); %Sliding window length for the layered structure
            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 data of Locations of sources and receivers
Rec=load('ReceiversPlot.txt'); % Receiver
Rec(:,2:4)=Rec(:,2:4)*0.3048; % Convert from foot to meter
Rec(:,4)=-Rec(:,4);

Sor=load('SourcesPlot.txt'); % Source
Sor(:,2:4)=Sor(:,2:4)*0.3048; % Convert from foot to meter
Sor(:,4)=-Sor(:,4);


Northing=[1 0 0]; % Northing axis (used to estimate back azimuth)
Vert=[0 0 1]; % Vertical axis (used to estimate incidence angle)

    
    numvars=6; % Number of components;


    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,:);
    
    Er=Rn.^2+Re.^2+Rz.^2;
    Er_m=max(abs(Er));

    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;
    

id=round(d/dts);
id2=round(id/2);
id_m=round(id/2);


event_window(1)=min(Ttime)+(id2+1)*dts;
event_window(2)=max(Ttime)-(id2+1)*dts;

i0=round(event_window(1)/dts);
ishift=round(min(Ttime)/dts);
i0=i0-ishift;

i1=round(event_window(2)/dts);
i1=i1-ishift;


ta=Ttime(i0);
te=Ttime(i1);

i=i0;
ii=1;


while i < i1,
    ia=i-id2+1;iz=i+id2; 
    if iz<nts,
        trt(ii)=Ttime(i);
        Er_m_i=max(abs(Er(ia:iz)));
        if Er_m_i < 0.01*Er_m
            baz(ii)=NaN;
            inc(ii)=NaN;
        else

            Er_w=Er(ia:iz)./max(abs(Er(ia:iz)));
        
            % acceleration
            An_w=An(ia:iz).*Er_w;
            Ae_w=Ae(ia:iz).*Er_w;
            Az_w=Az(ia:iz).*Er_w;

            % rotation rate
            Rn_w=Rn(ia:iz).*Er_w;
            Re_w=Re(ia:iz).*Er_w;
            Rz_w=Rz(ia:iz).*Er_w;



            % Cross product of acceleration and rotation rate
            for k=1:length(An_w)
                Rdir=cross([Rn_w(k) Re_w(k) Rz_w(k)],...
                        [An_w(k) Ae_w(k) Az_w(k)]);
                Rdir_n(k)=norm(Rdir);
                teta(k)=dot(Rdir, Vert)/norm(Rdir);
                teta(k)=acos(teta(k))/pi*180;
                phi(k)=dot(-1.*[Rdir(1) Rdir(2) 0], Northing)/norm([Rdir(1) Rdir(2) 0]);
                phi(k)=acos(phi(k))/pi*180;
            end
            
            
            % Weighting factor            
            r_min=min(Rdir_n);
            for k=1:length(An_w)
                ft_r(k)=round(((Rdir_n(k))^2)/(r_min^2));
            end

           
            % Estimated azimuth and incidence angle for each sliding time window
            teta_m=sum(teta.*ft_r)/sum(ft_r);
            phi_m=sum(phi.*ft_r)/sum(ft_r);;

            % Azimuth (from 0-180)
            baz(ii)=phi_m;

            % Incidence 
            inc(ii)=teta_m; 

            
        end
        
        ii=ii+1;
    end
    i=i+id_m;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Rec_ind=N_Rec;
    
% Locations of sources and receivers (geometry input)
    Rec_i=[Rec(Rec_ind,2),Rec(Rec_ind,3),Rec(Rec_ind,4)]; % Receiver
    
    S_I=[13 1 6 17 9];  % Source index according to the SLB's order
                        % coressponding with sources [1 2 3 4 5] used in the report

    Sor_ind=S_I(Res);
    Sor_i=[Sor(Sor_ind,2),Sor(Sor_ind,3),Sor(Sor_ind,4)]; % Source
    

% Theoretical azimuth estimated from locations of sources and receivers
        Theo_Dir=[Rec_i(1)-Sor_i(1) Rec_i(2)-Sor_i(2) Rec_i(3)-Sor_i(3)];
        

        Theo_Bang=-1.*[Theo_Dir(1) Theo_Dir(2) 0]; % Use two horizontal components (of backward direction) only

        Theo_B_ang=dot(Theo_Bang,Northing)/(norm(Theo_Bang)*norm(Northing)); % Using the definition of dot product 
            b_azo=acos(Theo_B_ang)/pi*180;

% Theoretical incidence estimated from locations of sources and receivers
        Theo_I_ang=dot(Theo_Dir,Vert)/(norm(Theo_Dir)*norm(Vert)); % Using the definition of dot product 
        inc_o=acos(Theo_I_ang)/pi*180;


% Plot figure
        
figure(1)
subplot(511),
plot(Ttime(i0:i1),Vn(i0:i1)./max(abs(Vn(i0:i1))),'k-'), text(ta+0.02, 0.4,'V_N');
set(gca,'YTick',[]);
axis([ta, te, -1, 3]);box on;
hold on;
plot(Ttime(i0:i1),1+Ve(i0:i1)./max(abs(Ve(i0:i1))),'b-'), text(ta+0.02, 1.4,'V_E');
hold on;
plot(Ttime(i0:i1),2+Vz(i0:i1)./max(abs(Vz(i0:i1))),'r-'), text(ta+0.02, 2.4,'V_Z');

subplot(512),
plot(Ttime(i0:i1),Rn(i0:i1)./max(abs(Rn(i0:i1))),'k-'), text(ta+0.02, 0.4,'R_N');
set(gca,'YTick',[]);
axis([ta, te, -1, 3]);box on;
hold on;
plot(Ttime(i0:i1),1+Re(i0:i1)./max(abs(Re(i0:i1))),'b-'), text(ta+0.02, 1.4,'R_E');
hold on;
plot(Ttime(i0:i1),2+Rz(i0:i1)./max(abs(Rz(i0:i1))),'r-'), text(ta+0.02, 2.4,'R_Z');

subplot(513),
plot(Ttime(i0:i1),Er(i0:i1)./max(abs(Er(i0:i1))),'r-'),
set(gca,'YTick',[]);
axis([ta te 0 1]); box on;
text(ta+0.02, 0.1,'Er');

subplot(514),
plot(trt,baz,'b*'), ylabel('az. (deg)'),axis([ta, te, 0,180]);
hold on;
plot(Ttime(i0:i1),ones(size(Ttime(i0:i1))).*b_azo,'r-');


subplot(515),
plot(trt,inc,'r+'), ylabel('incidence (deg)'),axis([ta, te, 0,180]);
hold on;
plot(Ttime(i0:i1),ones(size(Ttime(i0:i1))).*inc_o,'b-');
xlabel(' Time (s) ');

%pfile=sprintf('./Figs/Fig8_Baz_Inc_%i_%i_%i.jpg',mtype,Res,N_Rec);eval([' print -djpeg100 ' pfile ]);    

cd MatLabCodes; % Back to program directory

