% SLB Project 
% Program to calculate cross correlation between Rz and At  

% Pham D. N., H. Igel, and M. Kaeser, 
% Geophysics, LMU Munich.
% Jan. 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';
            Tmax=1/15;
            Res=input('Enter the source index (1 to 5):');
        else
            f_name1='Layer';
            Tmax=1/30;
            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 theoretical back-azimuth
baz_5sour=load('Baz_z0_theo.txt');
b_azo=baz_5sour(N_Rec,Res);
baz=b_azo*pi/180;




numvars=6; % Number of variables;

% 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=SEIS(1,:);
Vn=SEIS(2,:);
Ve=SEIS(3,:);
Vz=SEIS(4,:);
%Rn=SEIS(5,:);
%Re=SEIS(6,:);
Rz=SEIS(7,:);

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;
    At=An*sin(baz)-Ae*cos(baz);


% vertical acceleration
Az=diff(Vz)/dts;
Az(length(Az)+1)=0;


%%%%%%%%%%%%%%%%%%%%%%
noise_amp=0.03;
% Generate a series of random noise perturb from -1 to 1 around a mean value of 0
pert1=(rand(1,length(Ttime))-0.5).*2;
% Max. Amp. of random noise is equal to 3% of max. Amp. of Rz
pert_Rz=noise_amp*max(Rz).*pert1;

pert2=(rand(1,length(Ttime))-0.5).*2;
% Max. Amp. of random noise is equal to 3% of max. Amp. of At
pert_At=noise_amp*max(At).*pert2;

At_n=At+pert_At;
Rz_n=Rz+pert_Rz;

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


% Filtering   

addpath .\MatLabCodes;
disp(' Filtering ... ')
fAt=highpass(At_n,dts,1/Tmax,1);
fRz=highpass(Rz_n,dts,1/Tmax,1);



% Quantify semblance between data
% d: length of time window (s) for correlation
d=2*Tmax;
id=round(d/dts);
id2=round(id/2);
id_m=round(id/10); % Sliding step


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

% Extracting window
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;
disp(' Calculating semblance R-T ');

while i < i1,
    ia=i-id2+1;iz=i+id2; 
    if iz<nts,
    trt(ii)=Ttime(i);
    xrt(ii)=max(xcorr(fAt(ia:iz),-fRz(ia:iz),'coeff'));
    ii=ii+1;
    end
    i=i+id_m;
end


figure(1)
subplot(411),plot(Ttime(i0:i1),Az(i0:i1),'k-'), ylabel(' A_Z (m/s^2)'),axis([ta, te, -max(abs(Az(i0:i1))), max(abs(Az(i0:i1)))])
subplot(412),plot(Ttime(i0:i1),At(i0:i1),'b-'), ylabel(' A_T (m/s^2) '),axis([ta, te, -max(abs(At(i0:i1))), max(abs(At(i0:i1)))])
subplot(413),plot(Ttime(i0:i1),Rz(i0:i1),'r-'), ylabel('R_Z (rad/s)'),axis([ta, te, -max(abs(Rz(i0:i1))), max(abs(Rz(i0:i1)))])
subplot(414),plot(trt,abs(xrt)), axis([ta, te, 0, 1]), ylabel('Corr. Coef.');
xlabel(' Time (s) ')
drawnow


%disp(' Plotting to file ... ')
%pfile=sprintf('./Figs/Fig5_Xco_AtRz_%i_%i_%i.eps',mtype,Res,N_Rec); eval([' print -depsc2 ' pfile ]);
%pfile=sprintf('./Figs/Fig5_Xco_AtRz_%i_%i_%i.jpg',mtype,Res,N_Rec); eval([' print -djpeg100 ' pfile ]);

cd MatLabCodes; % Back to program directory

