%DISPLOV dessine vitesse phase et groupe pour onde de Love dans une couche
% d'epaisseur H et vitesse VS1 sur demi-espace de vitesse VS2
% voir : http://jmmeost.free.fr/Sismique/osis.html#LO
% fig1 = http://jmmeost.free.fr/Sismique/lovdisp.gif  
% fig2 = http://jmmeost.free.fr/Sismique/lovg.gif

clear all

V1 = 2000 ; V2 = 2600; H = 100;

% angles incidences   eta critique < eta1< 90 et parametres des ondes qui
% rebondissent dans la couche
etac = asin(V1/V2); xc = 2*H*tan(etac);
eta1 = asin(V1/V2):.001:pi/2; p = sin(eta1)/V1; 
% la vitesse de phase VL de l'onde de Love est la vitesse apparente
% horizontale des ondes dans la couche
VL = 1./p;
% demi-dephasage a la reflexion totale -chi(p)/2
I1 = V1*sqrt(1-(p*V1).^2); I2 = V2*sqrt((p*V2).^2-1);
C2 = atan(I2./I1);
% relation de dispersion omega(p) et frequences coupure 
omp0 = C2./sqrt(1-(p*V1).^2)*V1/H; % mode fondamental
omp1 = (C2+pi)./sqrt(1-(p*V1).^2)*V1/H; omc1 = pi/H/sqrt(1/V1^2-1/V2^2); % harmoniques 
omp2 = (C2+2*pi)./sqrt(1-(p*V1).^2)*V1/H; omc2 = 2*omc1;
omp3 = (C2+3*pi)./sqrt(1-(p*V1).^2)*V1/H; omc3 = 3*omc1;

figure(1), clf, hold on

plot(omp0/2/pi,VL,'b',omp1/2/pi,VL,'r',omp2/2/pi,VL,'g',omp3/2/pi,VL,'m','Linewidth',2)
xlim([0,100]), ylim([V1,V2])
plot(omc1/2/pi,V1,'r+',omc2/2/pi,V1,'g+',omc3/2/pi,V1,'m+')
title('VITESSE PHASE ONDE LOVE V1 =2000 M/S V2 =2600 M/S H =100 M')
legend('n = 0','n = 1','n = 2','n = 3','fc1','fc2','fc3')
xlabel('FREQUENCE (HZ)'), ylabel('VL (M/S)'), grid on, box on

% vitesse de groupe
J1 = I1/V1; J2 = I2/V2;
DPO = H/V1.*J1.^2; % dp/dom
% mode fondamental
DP0 = DPO./(H*omp0.*p*V1+V2/V1.*(cos(H*omp0.*J1/V1)).^2.*(p*V2^2./J2+p*V1^2.*J2/J1.^2));
VG0 = 1./(p+omp0.*DP0);
% harmoniques
DP1 = DPO./(H*omp1.*p*V1+V2/V1.*(cos(H*omp1.*J1/V1)).^2.*(p*V2^2./J2+p*V1^2.*J2/J1.^2));
VG1 = 1./(p+omp1.*DP1);

DP2 = DPO./(H*omp2.*p*V1+V2/V1.*(cos(H*omp2.*J1/V1)).^2.*(p*V2^2./J2+p*V1^2.*J2/J1.^2));
VG2 = 1./(p+omp2.*DP2);

DP3 = DPO./(H*omp3.*p*V1+V2/V1.*(cos(H*omp3.*J1/V1)).^2.*(p*V2^2./J2+p*V1^2.*J2/J1.^2));
VG3 = 1./(p+omp3.*DP3);

figure(2), clf, hold on
plot(omp0/2/pi,VL,'b',omp0/2/pi,VG0,'b--',omp1/2/pi,VL,'r',omp1/2/pi,VG1,'r--',omp2/2/pi,VL,'g',omp2/2/pi,VG2,'g--',omp3/2/pi,VL,'m',omp3/2/pi,VG3,'m--','Linewidth',2)
xlim([0,100]), ylim([1700,V2])
title('VITESSE PHASE VL ET GROUPE VG ONDE LOVE V1 =2000 M/S V2 =2600 M/S H =100 M')
legend('VL n=0','VG','VL n=1','VG','VL n=2','VG','VL n=3','VG')
xlabel('FREQUENCE (HZ)'), ylabel('VG (M/S) VL (M/S)'), grid on, box on

xmin = 2000; xmax = 10000;
dom = .2;
figure(3), clf

for ifr = 1:3
switch ifr 
    case 1, f1 = 2; f2 = 8; fm = 6; 
    case 2, f1 = 8; f2 = 16; fm = 12; 
    case 3, f1 = 16; f2 = 32; fm = 24;
end
om0 = [f1:dom:f2]*2*pi; nom = length(om0); dx = V2/f2/4; ifm = (fm-f1)/dom+1;
p0 = interp1(omp0,p,om0); g0 = interp1(omp0,VG0,om0); 
xx = xmin:dx:xmax; tt = 0:.002:1; [XX, TT] = meshgrid(xx,tt);
u0 = cos(om0(1)*(p0(1)*XX-TT-XX/V2));
for iom = 2:nom
u0 = u0+cos(om0(iom)*(p0(iom)*XX-TT-XX/V2));
end
u0 = u0/nom;
subplot(2,2,ifr), hold on
imagesc(xx,tt,u0)
plot([xmin xmax],[xmin xmax]/g0(ifm)-[xmin xmax]/V2,'k','LineWidth',2)
for icy = 1:f2
plot([xmin xmax],[xmin xmax]*p0(ifm)-[xmin xmax]/V2+(icy-1)/fm,'k')
end
xlim([xmin xmax-dx]), ylim([0 1]), caxis([-1.2 1.2]), colormap('jet'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('T - X / V2 (S)'),
title(['MODE FONDAMENTAL F = ' num2str(f1) '  -  ' num2str(f2) ' HZ'])
legend(['VG @ ' num2str(fm) ' HZ'],['VL @ ' num2str(fm) ' HZ'])
end


p1 = interp1(omp1,p,om0); g1 = interp1(omp1,VG1,om0); 
u1 = cos(om0(1)*(p1(1)*XX-TT-XX/V2));
for iom = 2:nom
u1 = u1+cos(om0(iom)*(p1(iom)*XX-TT-XX/V2));
end
u1 = u1/nom;
subplot(2,2,4),hold on
imagesc(xx,tt,u1)
plot([xmin xmax],[xmin xmax]/g1(ifm)-[xmin xmax]/V2,'k','LineWidth',2)
for icy = 1:f2
plot([xmin xmax],[xmin xmax]*p1(ifm)-[xmin xmax]/V2+(icy-1)/fm,'k')
end
xlim([xmin xmax-dx]), ylim([0 1]), caxis([-1 1]), colormap('jet'), set(gca,'YDir','reverse')
xlabel('X (M)'), ylabel('T - X / V2 (S)'),
title(['PREMIER HARMONIQUE F = ' num2str(f1) '  -  ' num2str(f2) ' HZ'])
legend(['VG @ ' num2str(fm) ' HZ'],['VL @ ' num2str(fm) ' HZ'])
gtext('ONDE DE LOVE V1 =2000 M/S V2 =2600 M/S H =100 M : VITESSES DE GROUPE VG ET DE PHASE VL')