% generate a synthetic radially polarized (e.g. SKS) wave, then split it 
% using prescribed dt and phi, 
% then unsplit it following Silver and Chan, 1991 to estimate dt and phi.
% J.A. Olive, MIT/WHOI Joint Program in Oceanography, Dec. 2010

clear all
close all

% time 
N=1500; % number of points 
Tmax=30; % (s) so that sampling frequency is 50 Hz
t=linspace(0,Tmax,N);
nus=Tmax/N; % sampling rate (s)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODIFY HERE
az=40; % azimuth (deg) of incoming sks wave polarization
% SHEAR WAVE SPLITTING PARAMETERS
dt=1.3; % delay time through anisotropic layer (s)
phi=80; % azimuth of fast direction (deg)

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


p=[sind(az); cosd(az)]; % incoming sks wave polarization vector
s=[-cosd(phi); sind(phi)]; % slow direction
f=[sind(phi); cosd(phi)]; % fast direction


% ISOTROPIC INCOMING SIGNAL
iso=exp(-(t-Tmax/2).^2/(Tmax/20)^2).*cos(2*pi*t/(Tmax/7)-2);

isoE=iso*sind(az);
isoN=iso*cosd(az);

isoR=zeros(1,length(isoE));
isoT=isoR;
Rgw=[sind(az) cosd(az); -cosd(az) sind(az)];
% convert corrected east-north into radial-transverse
for vv=1:length(isoE)
    geo=[isoE(vv); isoN(vv)];
    wav=Rgw*geo;
    isoR(vv)=wav(1);
    isoT(vv)=wav(2);
end


% SIMULATING SPLITTING FROM AN ANISOTROPIC LAYER, BY:
% 1/ APPLYING THE SHIFTING OPERATOR
fast=dot(p,f)*ttshift(iso,nus,-dt/2);
slow=dot(p,s)*ttshift(iso,nus,dt/2);
% 2/ PROJECTING THE SHIFTED TRACES ON FAST AND SLOW DIRECTIONS
splitE=fast*sind(phi)-slow*cosd(phi);
splitN=fast*cosd(phi)+slow*sind(phi);



splitR=zeros(1,length(splitE));
splitT=splitR;
Rgw=[sind(az) cosd(az); -cosd(az) sind(az)];
% convert corrected east-north into radial-transverse
for vv=1:length(splitE)
    geo=[splitE(vv); splitN(vv)];
    wav=Rgw*geo;
    splitR(vv)=wav(1);
    splitT(vv)=wav(2);
end


% AT THIS POINT, WE HAVE THE ORIGINAL SIGNAL 
% split in geographic coordinates: splitE, splitN
% split in wave coordinates: splitR, splitT


% HERE WE PERFORM THE UNSPLITTING METHOD OF Silver and Chan, 1991
% TO ESTIMATE dt AND phi FROM THE SPLITTED TRACES

[D,tdt,tphi]=split(splitE,splitN,nus);
minval=min(min(D'));
[ii jj]=find(D'==minval);
dtbest=tdt(jj);
phibest=tphi(ii);
hold on
plot(dtbest,phibest,'wx','linewidth',2)

disp('---- FAST AZIMUTH (deg)')
disp('real')
disp(phi)
disp('estimated')
disp(phibest)
disp('---- DELAY TIME (s)')
disp('real')
disp(dt)
disp('estimated')
disp(dtbest)




% SOME DIAGNOSTIC PLOTS

smp=nus;
time=t;
se=splitE;
sn=splitN;
sf=zeros(1,length(se));
ss=sf;
Rga=[sind(phibest) cosd(phibest); -cosd(phibest) sind(phibest)];
% go to fast-slow coordinates
for vv=1:length(se)
    geo=[se(vv); sn(vv)];
    ani=Rga*geo;
    sf(vv)=ani(1);
    ss(vv)=ani(2);
end
% time (un-)shift the fast-slow components
sf_sh=ttshift(sf,smp,+dtbest/2);
ss_sh=ttshift(ss,smp,-dtbest/2);
corE=zeros(1,length(sf_sh));
corN=corE;
% convert corrected fast and slow into east-north
for vv=1:length(sf_sh)
    ani=[sf_sh(vv); ss_sh(vv)];
    geo=Rga'*ani;
    corE(vv)=geo(1);
    corN(vv)=geo(2);
end
corR=zeros(1,length(sf_sh));
corT=corR;
% convert corrected east-north into radial-transverse
for vv=1:length(sf_sh)
    geo=[corE(vv); corN(vv)];
    wav=Rgw*geo;
    corR(vv)=wav(1);
    corT(vv)=wav(2);
end

% NOW WE HAVE THE UNSPLIT SIGNAL
% in geographic coordinates: corE, corN
% in wave coordinates: corR, corT
% TO BE COMPARED WITH THE ORIGINAL ISOTROPIC WAVE
% in geographic coordinates: isoE, isoN
% in wave coordinates: isoR, isoT


figure
subplot(211)
plot(t,isoR,'g--','linewidth',1.2)
hold on
plot(t,splitR,'r')
plot(t,corR,'k')
legend('isotropic','split','unsplit')
hold on
set(gca,'PlotBoxAspectRatio',[3 1 1])
grid on
ylabel('radial','fontsize',20)
xlabel('time (s)','fontsize',20)
subplot(212)
plot(t,isoT,'g--','linewidth',1.2)
hold on
plot(t,splitT,'r')
plot(t,corT,'k')
legend('isotropic','split','unsplit')
hold on
set(gca,'PlotBoxAspectRatio',[3 1 1])
grid on
ylabel('transverse','fontsize',20)
xlabel('time (s)','fontsize',20)