% This code is modified and fixed by Dr. Sherif Hanafy Mahmoud
%
% This Matlab script deals with
% 1) reciprocity error check ;
% 2) creat a 'coord.temp' file for traveltime tomography;
% -------------------------------------------------------------------
% Input:
% TTimes.dat --- picked traveltimes, which has 3 columns: 
%                 |src# |rec#| |time in Sec|
%
% Topo.dat --- Elevation file, contains two column |X in meters| and |Z in
% meters|. Note that, all receivers and sources need to be within the
% elevation data range and do not need to be the same, the current code can
% interpolate for sources and receivers locations, but can not extrapolate.
% e.x.
% o X   ^   X   ^ o X   ^ o X o ^   X   ^   X   ^   X o ^   X o
% here, o -> elevation point, X -> source point, and ^ -> receiver point;
% elevation point does not need to be regular, but I need to start the
% profile and end it with elevation point, sor and/or receiver can share
% the first and/or last elevation point.
%
% The X locations of the sources and the receivers mearsured in the 
%  field are apparent values. The current code will convert these apparent values
%  to true values using the following relations:
%
%
%                                                Measurement:(dis,elev)
%                                         ____/| zs=elev;
%                                    ____/     | xs=sqrt(dx^2-dz^2);
%                               ____/          | zg=elev;
%                              /_______________| xg=sqrt(dx^2-dz^2);
% 
% cutoff --- cutoff value for the maximum reciprocity error (default=0.004)
% tshift --- time shift between the picked and the first break traveltimes
%            If you pick the first zero-crossing, you don't need this.
%
%----------------------------------------------------------------------------
% OUTPUT:
% coord.temp

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

clear all;
close all;
clc;

% Load traveltime table, source coordinate, reciever coordinate.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Change the following values according to your data
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Travel time table is a 3 column ASCII file with Sor#, Rec#, and Time in
% Sec for 1st, 2nd, and 3rd column, respectively.
pick=load ('p_2.pik');

% Topography file is a 2-column ASCII file with X and Z values as measured
% in the field. Note, minimum X in topography file should be <= 1st
% shot/receiver X location (whichever is smaller), and maximum X in
% topography file should be >= last source/receiver (whichever is bigger).
Topo_D=load('topo.txt');

%   Input Values, change according to your data
windowSize=1;     % for topography smoothing, = 1 for no smoothing
ds=1;             % Source interval; if constant, if not ignore this value
s1_x=0;           % X position of first shot
ns=61;           % No. of sources
dg=1;             % Receiver interval; if constant, if not ignore this value
g1_x=0;           % X location of first receiver
ng=61;           % No of receivers
zmax=10;          % Tomogram max. required depth
cutoff=0.002;     % Accepted reciprocity error (in Sec)
tshift=0.000;     % Const. value to subtrace from all travel times (in Sec)
Rm_Rec=1;         % =1 to remove reciprocity error, = 0 to keep all data

Tr_unpicked=0;
for i=1:ns*ng
    sx=(pick(i,1)-1)*ds+s1_x;
    gx=(pick(i,2)-1)*dg+g1_x;
    if sx==gx
        pick(i,3)=0;
    end
    if pick(i,3)<0
        pick(i,3)=0;
    end
    if pick(i,3)==0
        Tr_unpicked=Tr_unpicked+1;
    end
end

%   If sor/rec intervals are constant, then
Rec_loc=(g1_x:dg:ng*dg+g1_x)';
Rec_loc=Rec_loc(1:ng);
Sor_loc=(s1_x:ds:ns*ds+s1_x)';
Sor_loc=Sor_loc(1:ns);
%   If sor/receiver are at variable distance, read them from an ASCII file
%   It is one column ASCII files shows the X location of receivers
%   (Rec_loc.dat) and the X location of the source (Sor_loc.dat).
% load Rec_loc.dat
% load Sor_loc.dat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       You should NOT do any changes in the following part
%       Unless you know what you are doing :)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   To Smooth topography data using a moving avrage equation. 
s(:,2)=filter(ones(1,windowSize)/windowSize,1,Topo_D(:,2));
s(1:windowSize,2)=Topo_D(1:windowSize,2);
s(:,1)=Topo_D(:,1);

%   Correct the Sor and Rec x and z from apparant values to true values
Rec_loc(:,2)=interp1(Topo_D(:,1),Topo_D(:,2),Rec_loc(:,1));
cr=Rec_loc*0;
cr(1,1)=Rec_loc(1,1);
for i=2:ng
    dx=Rec_loc(i,1)-Rec_loc(i-1,1);
    dz=Rec_loc(i,2)-Rec_loc(i-1,2);
    cr(i,1)=sqrt(dx^2-dz^2)+cr(i-1,1);
end
cr(:,2)=Rec_loc(:,2);
xr=cr;

Sor_loc(:,2)=interp1(Topo_D(:,1),Topo_D(:,2),Sor_loc(:,1));
cs=Sor_loc*0;
cs(1,1)=Sor_loc(1,1);
for i=2:ns
    dx=Sor_loc(i,1)-Sor_loc(i-1,1);
    dz=Sor_loc(i,2)-Sor_loc(i-1,2);
    cs(i,1)=sqrt(dx^2-dz^2)+cs(i-1,1);
end
cs(:,2)=Sor_loc(:,2);
xs=cs;

%   Find minimum and maximum values
xmin_src=min(cs(:,1));
xmin_rec=min(cr(:,1));
xmax_src=max(cs(:,1));
xmax_rec=max(cr(:,1));
zmin_src=min(cs(:,2));
zmin_rec=min(cr(:,2));

xmin=min(xmin_src,xmin_rec);
xmax=max(xmax_src,xmax_rec);
zmin=min(zmin_src,zmin_rec);

%   Make sure that the minimum in both X and Z are = 0
if xmin>0
    cr(:,1)=cr(:,1)-xmin;
    cs(:,1)=cs(:,1)-xmin;
    xmax=xmax-xmin;
    xmin=0;
end

if zmin>0
    cr(:,2)=cr(:,2)-zmin;
    cs(:,2)=cs(:,2)-zmin;
    zmin=0;
end
    
%%%% Before write into coord.temp, make a reciprocity check, and delete the
%%%% bad picks (> cutoff).
[nrow,ncol]=size(pick);
tt=zeros(ns,ng);

for i=1:nrow
    tt(pick(i,1),pick(i,2))=pick(i,3);
end

imagesc(tt);
title(['Original Travel Times. ',num2str(Tr_unpicked),' traces have not been picked']);
xlabel('Receiver #');ylabel('Source #');

tt_recip=tt;
Tr_recip=0;
% reciprocity check
for is=1:ns
    g1x=(is-1)*ds+s1_x;
    g1=round(g1x/dg)+1;
    for ig=g1+1:ng
        s2x=(ig-1)*dg+g1_x;
        s2=s2x/ds+1;
        g2x=(is-1)*ds+s1_x;
        g2=g2x/dg+1;
        if s2==round(s2) && g2==round(g2)
            if tt_recip(is,ig)>0 && tt_recip(s2,g2)>0
                t_diff=abs(tt_recip(is,ig)-tt_recip(s2,g2));
                if t_diff>cutoff
                    tt_recip(is,ig)=0;
                    tt_recip(s2,g2)=0;
                    Tr_recip=Tr_recip+1;
                end
            end
        end
    end
end

pick2=zeros(ns*ng,3);
k=1;
for i=1:ns
    for j=1:ng
        pick2(k,1)=i;
        pick2(k,2)=j;
        pick2(k,3)=tt_recip(i,j);
        k=k+1;
    end
end

if Rm_Rec==0
    pick2=pick;
end

b=size(pick2);
rm=0;
i=0;
while i+rm<=b(1)-1
    i=i+1;
    if pick2(i,3)<=0
        pick2(i:end-1,:)=pick2(i+1:end,:);
        rm=rm+1;
        i=i-1;
    end
end
pick2=pick2(1:end-rm,:);
pick2(:,3)=pick2(:,3)+tshift;
figure;
imagesc(tt_recip);
title(['Travel Times after reciprocity test. ', num2str(Tr_recip),' traces did not pass reciprocity test']);
xlabel('Receiver #');ylabel('Source #');

%   Find number of receivers/each shot
nr=zeros(ns,1);
[nrow,ncol]=size(pick2);

k=1;
for i=2:nrow
    if pick2(i,1)-pick2(i-1,1)==0
        k=k+1;
    else
        nr(pick2(i-1,1))=k;
        k=1;
    end
end
nr(ns)=k;
nrmax=max(nr);

for is=1:ns
    coord_src(is,1)=is;
    coord_src(is,2:3)=cs(is,:);
    coord_src(is,4)=nr(is);
end

for i=1:nrow
    coord_rec(i,1:2)=cr(pick2(i,2),:);
    coord_rec(i,3)=pick2(i,3);
    coord_rec(i,4)=1;
end

%%% write into coord.temp
file='coord.temp';
fid=fopen(file,'w','b'); 
fprintf(fid,'%4i %12.4f %12.4f %12.4f %12.4f\n',ns,xmin,xmax,xmax-xmin,zmin);
fprintf(fid,'%12.4f %12.4f %4i\n',zmax,max(coord_src(:,3))-min(coord_src(:,3)),nrmax);

a=size(coord_src);
for i=1:a(1)
    fprintf(fid,'%4i %12.4f %12.4f %4i\n',coord_src(i,1:4));
end

b=size(coord_rec);
for i=1:b(1)
    fprintf(fid,'%12.4f %12.4f %12.4f %4i\n',coord_rec(i,1:4));
end
fclose(fid);

tot_un_pick=Tr_unpicked+Tr_recip;
tot_pick=ng*ns-(Tr_unpicked+Tr_recip);
per1=tot_un_pick/(ng*ns);
per1=(round(per1*1000)/10);
per2=tot_pick/(ng*ns);
per2=(round(per2*1000)/10);

display('Statistics');
display(['Unpicked traces: ', num2str(Tr_unpicked)]);
display(['Traces did''t pass reciprocity : ', num2str(Tr_recip)]);
display(['Total unpicked traces: ', num2str(tot_un_pick)]);
display(['Number of picked traces: ',num2str(tot_pick)]);
display([num2str(per1),' % rejected, and ',num2str(per2),' % accepted']);