%MIGSVZ migre à déport nul pour V(Z)= V0+az
% voir : http://jmmeost.free.fr/Sismique/isis.html
% fig1 = http://jmmeost.free.fr/Sismique/migsvz1.gif

clear all

V0 = 1500; aa = 2;

dx = 10; xl = 3500; xx = 0:dx:xl-dx; nx = length(xx);
dt = .008; tl = 3; tt = 0:dt:tl-dt; nt = length(tt); 
dz = V0*dt/2; zl = 1000; zz = 0:dz:zl-dz; nz = length(zz);

Vz = V0+aa*zz;

% position des points diffractants et reflecteur
%xd = [100]; nd = length(xd); zd = [100];
al = 45*pi/180; zd = [0 500]; xd = 2500+[0 zd(2)/tan(al)]; Vd = V0+aa*zd(2); ixd = fix(xd(1)/dx); nd = length(xd);
x0 = xd(1)-V0/aa/tan(al); % intersection du réflecteur avec z = -V0/a = centre des rayons réfléchis
xxm = x0+V0/aa/tan(asin(V0/Vd*sin(al))); ixm = fix(xxm/dx); % abscisse max des rayons réfléchis vers haut
xxmb = 2*x0-xxm; ixmb = fix(xxmb/dx); xxmc = 2*x0-xd(1); ixmc = fix(xxmc/dx); % abscisse max et min des rayons réfléchis vers bas


% fréquence fondamentale et nombre d'harmoniques
f0 = 1; om = 2*pi*f0; k0 = om/V0; nf = 10;

pxt = zeros(nt,nx);
% MODELISATION

% DIFFRACTIONS
% onde diffractée en xt
for id = nd
% calcul des paramètres de rayons p(xx-xd)
Vd = V0+aa*zd(id); Va2 = (V0-Vd)^2; Vb2 = (V0+Vd)^2; Xp = sqrt(Vd^2-V0^2)/aa; Tp = -2/aa*log(V0/Vd/(1+sqrt(1-(V0/Vd)^2))); iXp = fix(Xp/dx);
ax2 = (aa*(xx-xd(id))).^2; p2 = 4*ax2./(ax2+Va2)./(ax2+Vb2); pp = sqrt(p2);
% calcul des temps diffractés vers haut
tg0 = pp*V0./(1+sqrt(1-p2*V0^2)); tgz = pp*Vd./(1+sqrt(1-p2*Vd^2)); 
tmod = -2/aa*log(tg0./tgz); tmod((xx == xd(id))) = 2/aa*log(Vd/V0);
% calcul des temps diffractés vers bas
ixp = find(abs(xx-xd(id)) >= Xp); tmod(ixp) = tmod(ixp) - 4/aa*log(tgz(ixp));
% somme d'harmoniques + fenêtre gaussienne
amp = .2;
for ix = 1:nx
omt = om*(tt-tmod(ix)); ond = amp*cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
pxt(:,ix) = pxt(:,ix) + amp*(ond.*exp(-(nf*omt/2/pi).^2))';
end
end

%REFLEXION
% le paramètre des rayons réfléchis est sin(al)/Vz
tam = -log(tan(al/2))/aa;
% angle d'incidence des rayons réfléchis en z = 0
th0 = atan(V0/aa./(xx(ixd:ixm)-x0)); ta0 = -log(tan(th0/2))/aa;
% temps réfléchis par dessus
tmod = 2*(ta0-tam); 
for ix = ixd:ixm
omt = om*(tt-tmod(ix-ixd+1)); ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
pxt(:,ix) = pxt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end
% temps réfléchis par dessous
tmod = 2*(ta0+tam);
for ix = ixmb+1:ixmc
omt = om*(tt-tmod(ixmc-ix+1)); ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
pxt(:,ix) = pxt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end


% MIGRATION 
% temps diffracté
txz = zeros(nz,nx); pxz = zeros(nz,nx); pxz1 = zeros(nz,nx); pxtm = zeros(nt,nx); 
for iz = 1:nz
Va2 = (V0-Vz(iz))^2; Vb2 = (V0+Vz(iz))^2; Xz = sqrt(Vz(iz)^2-V0^2)/aa; Tz = -2/aa*log(V0/Vz(iz)/(1+sqrt(1-(V0/Vz(iz))^2)));
ax2 = (aa*xx).^2; p2 = 4*ax2./(ax2+Va2)./(ax2+Vb2); pp = sqrt(p2);
tg0 = pp*V0./(1+sqrt(1-p2*V0^2)); tgz = pp*Vz(iz)./(1+sqrt(1-p2*Vz(iz)^2)); 
txz(iz,2:nx) = -2/aa*log(tg0(2:nx)./tgz(2:nx)); txz(iz,1) = 2/aa*log(Vz(iz)/V0);
ixp = find(xx >= Xz); txz(iz,ixp) = txz(iz,ixp) - 4/aa*log(tgz(ixp));
end
% limite rayons vers haut et bas pour plot séparé des migrations par dessus et dessous
xhb = xd(2)-Xp; ixhb = fix(xhb/dx);
for ix = 1:nx
% temps diffractés symétriques par rapport à trace courante
tm = txz;
itm = round(tm/dt); itm((itm>nt-1)) = nt-1; itm((itm<1)) = 1; itm(isnan(itm) ==1) = 1;
itm = [itm(:,ix:-1:1) itm(:,2:nx-ix+1)]; 
% copie de la trace à migrer dans matrice pxtm
pxtm = pxt(:,ix)*ones(1,nx);
% migration
%if ix == iXp+fix(xd(id)/dx), pxz1 = pxtm(itm); end
if (ix <= ixhb), pxz = pxz+pxtm(itm); else,  pxz1 = pxz1+pxtm(itm); end
end
    
figure(1), clf, colormap('copper')
subplot(3,2,1), hold on
imagesc([0 xl],[0 zl],zeros(nz,nx))
% trace rayon diffracte issu de zd 
plot([xd(1) xd(1)],[0 zd(1)],'m', xd(2)-Xp,0,'go','Linewidth',2)
for pp = 1/Vd*[sin(30*pi/180) sin(45*pi/180) sin(60*pi/180) 1]
R = 1/aa/pp;  tg0 = pp*V0/(1+sqrt(1-(pp*V0)^2)); tgz = pp*Vd./(1+sqrt(1-(pp*Vd)^2)); 
th0 = 2*atan(tg0); thz = 2*atan(tgz); th = [th0:.001:thz thz:.001:pi-th0];
plot(R*(cos(th)-cos(thz))+xd(id),R*sin(th)-V0/aa,'m',R*(cos(th)+cos(thz))+xd(id),R*sin(th)-V0/aa,'m') %,'Linewidth',2);
end
% rayons réfléchis
za = 0:100:zd(2); xa = xd(1)+za/tan(al); na = length(za); plot(xa,za,'bo')
Vza = V0+aa*za; th0 = asin(V0./Vza*sin(al)); R = Vza/sin(al)/aa;
xam = R.*cos(al); xa0 = R.*cos(th0); xah = xa0-xam; xab = xa0+xam; % distance hor de za a z=0 pour rai vers haut et bas
for ia = 1:na
xrh = 0:10:xah(ia); zrh = sqrt(R(ia)^2-(xrh+xam(ia)).^2)+za(ia)-Vza(ia)/aa; % point courant du rai vers haut
xrb = 0:10:xab(ia); zrb = sqrt(R(ia)^2-(xrb-xam(ia)).^2)+za(ia)-Vza(ia)/aa; % point courant du rai vers bas
plot([xa(ia)+xrh xa(ia)+xah(ia)],[zrh 0],'b','Linewidth',2)
plot([xa(ia)-xrb xa(ia)-xab(ia)],[zrb 0],'r','Linewidth',2)
end
plot(xd,zd,'y','Linewidth',2)
xlabel('X (M)'), ylabel('Z (M)'), 
%title('POINTS DIFFRACTANTS ET RAYONS POUR Theta = 0, 30, 45, 60, 90°'), 
title(['REFLECTEUR {\alpha} = ' num2str(al*180/pi) ' ° , V = ' num2str(V0) ' + ' num2str(aa) ' Z M/S'])
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(3,2,3), hold on
imagesc([0 xl],[0 tl],[zeros(nt, ixhb)  pxt(:,ixhb+1:nx)])
plot([xhb xhb],[0 zl],'g','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), 
title(['REFLEXION/DIFFRACTIONS  PAR DESSUS (F = ' num2str(f0) ' - ' num2str(nf*f0) ' HZ)']), 
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .002 1])

subplot(3,2,5), hold on
imagesc([0 xl],[0 tl],[pxt(:,1:ixhb) zeros(nt, nx-ixhb)])
plot([xhb xhb],[0 zl],'g','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), 
title(['REFLEXION/DIFFRACTIONS  PAR DESSOUS  (F = ' num2str(f0) ' - ' num2str(nf*f0) ' HZ)']), 
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .002 1])

subplot(3,2,2), hold on
imagesc([0 xl],[0 zl],pxz+pxz1)
plot(xd,zd,'y','Linewidth',2)
xlabel('X (M)'), ylabel('Z (M)'), 
title('MIGRATION PAR DESSUS + PAR DESSOUS')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(3,2,4), hold on
imagesc([0 xl],[0 zl],pxz1)
plot(xd,zd,'y','Linewidth',2)
xlabel('X (M)'), ylabel('Z (M)'), 
title('MIGRATION PAR DESSUS')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(3,2,6), hold on
imagesc([0 xl],[0 zl],pxz)
plot(xd,zd,'y','Linewidth',2)
xlabel('X (M)'), ylabel('Z (M)'), 
title('MIGRATION PAR DESSOUS')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])