% SLB Project 

% Program to figure out the geometry set up for SLB Project 

% then estimate the theoretical back-azimuth 
% from locations of sources and receivers

% Pham D. N., H. Igel, M. Kaeser
% Jan. 2009

clear all;
close all;

cd ..; % Shift to data directory

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);

Sta1=load('Stage1.txt');Sta1=Sta1*0.3048;
Sta2=load('Stage2.txt');Sta2=Sta2*0.3048;

plot3(Sta1(:,2),Sta1(:,1),Sta1(:,3),'k+','MarkerSize',5);
hold on;
plot3(Sta2(:,2),Sta2(:,1),Sta2(:,3),'b+','MarkerSize',5);
hold on;
plot3(Rec(:,2),Rec(:,3),Rec(:,4),'gv');
hold on;
plot3(Sor([13 1 6 17 9],2),Sor([13 1 6 17 9],3),Sor([13 1 6 17 9],4),'ro','MarkerSize',12);

grid on;
xlabel('Northing (m)');ylabel('Easting (m)');zlabel('Depth (m)');


for aaa=[1 8 9 15],  % Receiver index
    if aaa<9
        h=text(Rec(aaa,2),Rec(aaa,3)-35,Rec(aaa,4),sprintf('%i',aaa));set(h,'FontSize', 14);
    else
        h=text(Rec(aaa,2),Rec(aaa,3),Rec(aaa,4)-35,sprintf('%i',aaa));set(h,'FontSize', 14);
    end
end


ccc=[13 1 6 17 9]; % Source index according to the SLB's order
for bbb=1:5
    bbbb=ccc(bbb);
    h=text(Sor(bbbb,2),Sor(bbbb,3),Sor(bbbb,4)-35,sprintf('%i',bbb));set(h,'Color','red','FontSize', 14);
end



% Estimating the back-azimuth 
% baz=angle between Ox(1,0,0) and the projection of (Sor, Rec) on xOy

for sss=1:5, % Loop for different sources
    Sor_b=Sor(ccc(sss),2:4);
    
    for rrr=1:15, % Loop for different receivers

        Rec_b=Rec(rrr,2:4);
        Evt=-1.*[Rec_b(1)-Sor_b(1) Rec_b(2)-Sor_b(2) 0]; % Backward direction

        Northing=[1 0 0]; % Unit vector of Northing axis

        C_ang=dot(Evt,Northing)/(norm(Evt)*norm(Northing)); % Using the definition of dot product 

        ax=cross(Evt,Northing);
        if ax(3)<0
            ang(sss,rrr)=acos(C_ang)/pi*180;
        else
            ang(sss,rrr)=360-acos(C_ang)/pi*180;
        end
        
        
    end  % End the loop for different receivers
    
end  % End the loop for different sources


% Save the calculated result to file
%ang=ang'; save Baz_z0_theo.txt -ASCII -TABS ang;


cd MatLabCodes; % Back to the program directory