% SLB Project 
% Program to estimate S wave slowness
% from joint measurements of rotational and translational motions
% at linear receivers in a vertical borehole

% Then, back-azimuths are extracted and quantified 
% DEMONSTRATION THE BACK AZIMUTH ERROR 


% CASE STUDY WITH ALL SOURCES AND RECEIVERS IN THE VERTICAL BOREHOLE
% A SLIDING TIME WINDOW IS USED


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


clear all;
close all;

cd ..; % Shift to data directory


% Load the Locations of the receivers
Rec=load('ReceiversPlot.txt'); % Receiver
Rec(:,2:4)=Rec(:,2:4)*0.3048; % Convert from foot to meter
Rec(:,4)=-Rec(:,4);

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

% Distances between receivers in the vertical borehole
for ddd=1:7
    dz(ddd)=Rec(ddd,4)-Rec(ddd+1,4);
end

% Load the information of the theoretical back azimuth 
baz_5sour=load('Baz_z0_theo.txt');


Northing=[1 0 0]; % Northing is coincided with x axis    
%Vert=[0 0 1]; % Vertical axis        



% 1) Load the simulated seismograms

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';
        d=0.2*(1/15); %Sliding window length for the Homo. case
        Res=sss; % Source index
        f_name2=[num2str(Res)];
    else
        f_name1='Layer';
        d=0.2*(1/30); %Sliding window length for the layered case
        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:8 % Loop for 8 receivers in the vertical borehole
    
        N_Rec=rrr;
        f_name3=[num2str(N_Rec)];

        % 1) Load the data of the simulated seismograms

        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;


        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);
        nts(rrr)=length(Ttime(rrr,:));

        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;
    

    end



    % Index of Receiver location that we choose to investigate
    for Rec_ind=1:8 % Loop for 8 receivers in the vertical borehole

        % 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


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

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

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

        i1=round(event_window(2)/dts(Rec_ind));
        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(Rec_ind),
                trt(ii)=Ttime(Rec_ind,i);
                Er_m_i(ii)=max(abs(Er(Rec_ind,ia:iz)));
                if Er_m_i(ii) < 0.01*max(abs(Er(Rec_ind,:)));
                    baz(ii)=NaN;
                    %inc(ii)=NaN;
                else

                    if Rec_ind==1
                        dVx(ii)=-(mean(Vn(Rec_ind,ia:iz))-mean(Vn(Rec_ind+1,ia:iz)));
                        Sz1(ii)=(dVx(ii)/dz(Rec_ind))/mean(An(Rec_ind,ia:iz));
    
                        dVy(ii)=-(mean(Ve(Rec_ind,ia:iz))-mean(Ve(Rec_ind+1,ia:iz)));
                        Sz2(ii)=(dVy(ii)/dz(Rec_ind))/mean(Ae(Rec_ind,ia:iz));
                        Sz(ii)=(Sz1(ii)+Sz2(ii))/2;
                    elseif Rec_ind==8
                        dVx(ii)=-(mean(Vn(Rec_ind-1,ia:iz))-mean(Vn(Rec_ind,ia:iz)));
                        Sz1(ii)=(dVx(ii)/dz(Rec_ind-1))/mean(An(Rec_ind,ia:iz));
    
                        dVy(ii)=-(mean(Ve(Rec_ind-1,ia:iz))-mean(Ve(Rec_ind,ia:iz)));
                        Sz2(ii)=(dVy(ii)/dz(Rec_ind-1))/mean(Ae(Rec_ind,ia:iz));
                        Sz(ii)=(Sz1(ii)+Sz2(ii))/2;
                    else
                        dVx1(ii)=-(mean(Vn(Rec_ind,ia:iz))-mean(Vn(Rec_ind+1,ia:iz)));
                        Sz1(ii)=(dVx1(ii)/dz(Rec_ind))/mean(An(Rec_ind,ia:iz));
                        dVx2(ii)=-(mean(Vn(Rec_ind-1,ia:iz))-mean(Vn(Rec_ind,ia:iz)));
                        Sz2(ii)=(dVx2(ii)/dz(Rec_ind-1))/mean(An(Rec_ind,ia:iz));
                        dVy1(ii)=-(mean(Ve(Rec_ind,ia:iz))-mean(Ve(Rec_ind+1,ia:iz)));
                        Sz3(ii)=(dVy1(ii)/dz(Rec_ind))/mean(Ae(Rec_ind,ia:iz));
                        dVy2(ii)=-(mean(Ve(Rec_ind-1,ia:iz))-mean(Ve(Rec_ind,ia:iz)));
                        Sz4(ii)=(dVy2(ii)/dz(Rec_ind-1))/mean(Ae(Rec_ind,ia:iz));
                        Sz(ii)=(Sz1(ii)+Sz2(ii)+Sz3(ii)+Sz4(ii))/4;
            
                    end
        
                    sAx(ii)=mean(An(Rec_ind,ia:iz));
                    sAy(ii)=mean(Ae(Rec_ind,ia:iz));
                    sAz(ii)=mean(Az(Rec_ind,ia:iz));
                    sRx(ii)=mean(Rn(Rec_ind,ia:iz));
                    sRy(ii)=mean(Re(Rec_ind,ia:iz));
                    sRz(ii)=mean(Rz(Rec_ind,ia:iz));

                    Sy(ii)=(sAy(ii)*Sz(ii)-2*sRx(ii))/sAz(ii);
                    Sx(ii)=(sAx(ii)*Sz(ii)+2*sRy(ii))/sAz(ii);

        

                    % Back azimuth estimation
                    Est_Rdir=-1.*[Sx(ii) Sy(ii) 0]; % Use two horizontal components only

                    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(ii)=acos(C_ang)/pi*180;
                    else
                        baz(ii)=360-acos(C_ang)/pi*180;
                    end

                    % Theoretical back azimuth estimated from locations of sources and receivers
                    b_azo=baz_5sour(Rec_ind,Res);
                    baz_Err(ii)=baz(ii)-b_azo;    
    
                end
        
                ii=ii+1;
            end
            i=i+id_m;
        end

        % Back-azimuth error at the window that Er -> max
        [dvS idS]=max(Er_m_i);
        baz_Err_f(sss,Rec_ind)=baz_Err(idS);


    end

end         % End of the loop for different sources & media


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

figure(1)
Rec=1:8;
for pf=1:length(s_ind)
    plot(Rec,baz_Err_f(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');

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


cd MatLabCodes; % Back to program directory

