%SDR illustre la propagation-migration pour un point diffractant
% voir : http://jmmeost.free.fr/Sismique/isis.html
% fig1 = http://jmmeost.free.fr/Sismique/sdr1.gif  
% fig2 = http://jmmeost.free.fr/Sismique/sdm1.gif

clear all

V0 = 2000;

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

xz = zeros(nz,nx); xt = zeros(nt,nx);

% fréquence fondamentale et nombre d'harmoniques
f0 = 1; om = 2*pi*f0; k0 = om/V0; nf = 20;
% position du point diffractant
xd = 500; zd = 300; for ii = xd/dx+[0 1]; xz(zd/dz+[0 1],ii) = 1; end
% position de la source
xs = 0; zs = 0; ts = sqrt((xd-xs)^2+(zd-zs)^2)/V0;

% PROPAGATION
% T = TS
% distance source-diffracteur
rsd = sqrt((xd-xs)^2+(zd-zs)^2);
% onde directe en xt
for ix = max([1 fix((-rsd+xs)/dx)]):min([nx fix((rsd+xs)/dx)])
tdir = abs(xx(ix)-xs)/V0; omt = om*(tt-tdir);
ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
xt(:,ix) = xt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end
xt = xt/nf;
% onde directe en xz
zdir = sqrt(rsd^2-(xx-xs).^2); zdir(find(imag(zdir)~=0)) = NaN;
for ix = 1:min([nx fix((rsd+xs)/dx)])
kz = k0*zdir(ix)/rsd; kzz = kz*(zz-zdir(ix));
ond = cos(kzz); for ii = 2:nf, ond = ond+cos(ii*kzz); end
xz(:,ix) = xz(:,ix) + (ond.*exp(-(nf*kzz/2/pi).^2))';
end
xz = xz/nf;

figure(1), clf, colormap('copper')
subplot(2,3,1), hold on
imagesc([0 xl],[0 zl],xz)
plot(xs,zs,'bo','MarkerFaceColor','b'), plot(xd,zd,'ro','MarkerFaceColor','r'), plot(xx,zdir,'b','Linewidth',2)
xlabel('X (M)'), ylabel('Z (M)'), title('ONDE DIR. A TS')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,4), hold on
imagesc([0 xl],[0 tl],xt), plot([0 xs],[xs/V0 0],'b',[xs xl],[0 (xl-xs)/V0],'b','Linewidth',2)
plot(xx,(rsd+sqrt((xx-xd).^2+zd^2))/V0,'r','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), title('ONDE DIR. A Z=0')
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

% T = TS+TR
xz = zeros(nz,nx); xt = zeros(nt,nx);
zr = zd;
% onde directe en xt
for ix = max([1 fix((-rsd-zr+xs)/dx)]):min([nx fix((rsd+zr+xs)/dx)])
tdir = abs(xx(ix)-xs)/V0;  omt = om*(tt-tdir);
ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
xt(:,ix) = xt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end
% onde diffractéee en xt dans zone Fresnel
xfre = sqrt(zr*V0/f0/nf/2);
for ix = ceil((xd-xfre)/dx):fix((xd+xfre)/dx)
tdif = (rsd+sqrt((xx(ix)-xd).^2+zr^2))/V0; omt = om*(tt-tdif);
ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
xt(:,ix) = xt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end

xt = xt/nf;
% onde directe en xz
zdir = sqrt((rsd+zr)^2-(xx-xs).^2); zdir(find(imag(zdir)~=0)) = NaN;
for ix = 1:min([nx fix((rsd+zr+xs)/dx)])
kz = k0*zdir(ix)/rsd; kzz = kz*(zz-zdir(ix));
ond = cos(kzz); for ii = 2:nf, ond = ond+cos(ii*kzz); end
xz(:,ix) = xz(:,ix) + (ond.*exp(-(nf*kzz/2/pi).^2))';
end
% onde diffractée en xz
zdif = sqrt(zr^2-(xx-xd).^2); non = find(imag(zdif)~=0); zdif(non) = NaN;
for ix = ceil((xd-zr)/dx)+2:fix((xd+zr)/dx)
kz = k0*zdif(ix)/zr; kzzm = kz*(zz-zd+zdif(ix)); kzzd = kz*(zz-zd-zdif(ix));
ondm = cos(kzzm); for ii = 2:nf, ondm = ondm+cos(ii*kzzm); end
ondd = cos(kzzd); for ii = 2:nf, ondd = ondd+cos(ii*kzzd); end
xz(:,ix) = xz(:,ix) + (ondm.*exp(-(nf*kzzm/2/pi).^2))'+ (ondd.*exp(-(nf*kzzd/2/pi).^2))';
end

xz = xz/nf;

subplot(2,3,2), hold on
imagesc([0 xl],[0 zl],xz)
plot(xx,zdir,'b',xx,zd-zdif,'r',xx,zd+zdif,'r','Linewidth',2)
plot(xs,zs,'bo','MarkerFaceColor','b'), plot(xd,zd,'ro','MarkerFaceColor','r'), plot(xd,0,'go','MarkerFaceColor','g')
xlabel('X (M)'), ylabel('Z (M)'), title('ONDE DIR + DIFR. A TS+TR')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,5), hold on
imagesc([0 xl],[0 tl],xt), plot([0 xs],[xs/V0 0],'b',[xs xl],[0 (xl-xs)/V0],'b','Linewidth',2)
plot(xx,(rsd+sqrt((xx-xd).^2+zd^2))/V0,'r','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), title('ONDE DIR + DIFR.  A Z=0')
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

% T = TS+TR


xz = zeros(nz,nx); xt = zeros(nt,nx);

rdl = sqrt(xd^2+zd^2)-100;
% onde directe en xt
for ix = 1:nx
tdir = abs(xx(ix)-xs)/V0;  omt = om*(tt-tdir);
ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
xt(:,ix) = xt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end
% onde diffractéee en xt 
for ix = 1:nx
tdif = (rsd+sqrt((xx(ix)-xd).^2+zd^2))/V0; omt = om*(tt-tdif);
ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
xt(:,ix) = xt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end

xt = xt/nf;
% onde directe en xz
zdir = sqrt((rsd+rdl)^2-(xx-xs).^2); zdir(find(imag(zdir)~=0)) = NaN;
for ix = 1:nx
kz = k0*zdir(ix)/rsd; kzz = kz*(zz-zdir(ix));
ond = cos(kzz); for ii = 2:nf, ond = ond+cos(ii*kzz); end
xz(:,ix) = xz(:,ix) + (ond.*exp(-(nf*kzz/2/pi).^2))';
end
% onde diffractéee en xz
zdif = sqrt(rdl^2-(xx-xd).^2); non = find(imag(zdif)~=0); zdif(non) = NaN;
for ix = 1:nx
kz = k0*zdif(ix)/zd; kzzm = kz*(zz-zd+zdif(ix)); kzzd = kz*(zz-zd-zdif(ix));
ondm = cos(kzzm); for ii = 2:nf, ondm = ondm+cos(ii*kzzm); end
ondd = cos(kzzd); for ii = 2:nf, ondd = ondd+cos(ii*kzzd); end
xz(:,ix) = xz(:,ix) + (ondm.*exp(-(nf*kzzm/2/pi).^2))'+ (ondd.*exp(-(nf*kzzd/2/pi).^2))';
end

xz = xz/nf;

subplot(2,3,3), hold on
imagesc([0 xl],[0 zl],xz)
plot(xx,zdir,'b',xx,zd-zdif,'r',xx,zd+zdif,'r','Linewidth',2)
plot(xs,zs,'bo','MarkerFaceColor','b'), plot(xd,zd,'ro','MarkerFaceColor','r'), plot(sqrt(rdl^2-zd^2)+xd,0,'go','MarkerFaceColor','g')
xlabel('X (M)'), ylabel('Z (M)'), title('ONDE DIR + DIFR.  A TS+TR')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,6), hold on
imagesc([0 xl],[0 tl],xt), plot([0 xs],[xs/V0 0],'b',[xs xl],[0 (xl-xs)/V0],'b','Linewidth',2)
plot(xx,(rsd+sqrt((xx-xd).^2+zd^2))/V0,'r','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), title('ONDE DIR + DIFR.  A Z=0')
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

% MIGRATION 

figure(2), clf, hold on, colormap('copper')

% T = TS+TR


xz = zeros(nz,nx); xt = zeros(nt,nx);

rdl = sqrt(xd^2+zd^2)-100;
% onde diffractéee en xt
for ix = 1:nx
tdif = (rsd+sqrt((xx(ix)-xd).^2+zd^2))/V0; omt = om*(tt-tdif);
ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
xt(:,ix) = xt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end

xt = xt/nf;
% onde diffractéee en xz
zdif = sqrt(rdl^2-(xx-xd).^2); non = find(imag(zdif)~=0); zdif(non) = NaN;
for ix = 1:nx
kz = k0*zdif(ix)/zd; kzzm = kz*(zz-zd+zdif(ix)); kzzd = kz*(zz-zd-zdif(ix));
ondm = cos(kzzm); for ii = 2:nf, ondm = ondm+cos(ii*kzzm); end
ondd = cos(kzzd); for ii = 2:nf, ondd = ondd+cos(ii*kzzd); end
xz(:,ix) = xz(:,ix) + (ondm.*exp(-(nf*kzzm/2/pi).^2))'+ (ondd.*exp(-(nf*kzzd/2/pi).^2))';
end

xz = xz/nf;

subplot(2,3,1), hold on
imagesc([0 xl],[0 zl],xz)
plot(xs,zs,'bo','MarkerFaceColor','b'), plot(xd,zd,'ro','MarkerFaceColor','r'), plot(sqrt(rdl^2-zd^2)+xd,0,'go','MarkerFaceColor','g')
xlabel('X (M)'), ylabel('Z (M)'), title('ONDE DIFR. A TS+TR')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,4), hold on
imagesc([0 xl],[0 tl],xt), plot([0 xl],[ts ts],'b','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), title('ONDE DIFR.  A Z=0')
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

% T = TS+TR
xz = zeros(nz,nx); xt = zeros(nt,nx);

zr = zd/2;
% onde diffractéee en xt
xfre = sqrt(zr*V0/f0/nf/2);
for ix = ceil((xd-xfre)/dx):fix((xd+xfre)/dx)
tdif = (rsd+sqrt((xx(ix)-xd).^2+zr^2))/V0; omt = om*(tt-tdif);
ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
xt(:,ix) = xt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end

xt = xt/nf;

% onde diffractéee en xz
zdif = sqrt(zr^2-(xx-xd).^2); non = find(imag(zdif)~=0); zdif(non) = NaN;
for ix = ceil((xd-zr)/dx)+2:fix((xd+zr)/dx)
kz = k0*zdif(ix)/zr; kzzm = kz*(zz-zd+zdif(ix)); kzzd = kz*(zz-zd-zdif(ix));
ondm = cos(kzzm); for ii = 2:nf, ondm = ondm+cos(ii*kzzm); end
ondd = cos(kzzd); for ii = 2:nf, ondd = ondd+cos(ii*kzzd); end
xz(:,ix) = xz(:,ix) + (ondm.*exp(-(nf*kzzm/2/pi).^2))'+ (ondd.*exp(-(nf*kzzd/2/pi).^2))';
end

xz = xz/nf;

subplot(2,3,2), hold on
imagesc([0 xl],[0 zl],xz)
plot(xs,zs,'bo','MarkerFaceColor','b'), plot(xd,zd,'ro','MarkerFaceColor','r'), plot(xd,zd-zr,'go','MarkerFaceColor','g')
xlabel('X (M)'), ylabel('Z (M)'), title('ONDE DIFR. A TS+TR')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,5), hold on
imagesc([0 xl],[0 tl],xt), plot([0 xl],[ts ts],'b','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), title('ONDE DIFR.  A Z=ZD/2')
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

% T = TS+TR
xz = zeros(nz,nx); xt = zeros(nt,nx);
zr = zd/10;

% onde diffractéee en xt
xfre = sqrt(zr*V0/f0/nf/2);
for ix = ceil((xd-xfre)/dx):fix((xd+xfre)/dx)
tdif = (rsd+sqrt((xx(ix)-xd).^2+zr^2))/V0; omt = om*(tt-tdif);
ond = cos(omt); for ii = 2:nf, ond = ond+cos(ii*omt); end
xt(:,ix) = xt(:,ix) + (ond.*exp(-(nf*omt/2/pi).^2))';
end

xt = xt/nf;
% onde diffractéee en xz
zdif = sqrt(zr^2-(xx-xd).^2); non = find(imag(zdif)~=0); zdif(non) = NaN;
for ix = ceil((xd-zr)/dx)+2:fix((xd+zr)/dx)
kz = k0*zdif(ix)/zr; kzzm = kz*(zz-zd+zdif(ix)); kzzd = kz*(zz-zd-zdif(ix));
ondm = cos(kzzm); for ii = 2:nf, ondm = ondm+cos(ii*kzzm); end
ondd = cos(kzzd); for ii = 2:nf, ondd = ondd+cos(ii*kzzd); end
xz(:,ix) = xz(:,ix) + (ondm.*exp(-(nf*kzzm/2/pi).^2))'+ (ondd.*exp(-(nf*kzzd/2/pi).^2))';
end

xz = xz/nf;

subplot(2,3,3), hold on
imagesc([0 xl],[0 zl],xz)
plot(xs,zs,'bo','MarkerFaceColor','b'), plot(xd,zd,'ro','MarkerFaceColor','r'), plot(xd,zd-zr,'go','MarkerFaceColor','g')
xlabel('X (M)'), ylabel('Z (M)'), title('ONDE DIFR. A TS+TR')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,6), hold on
imagesc([0 xl],[0 tl],xt), plot([0 xl],[ts ts],'b','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), title('ONDE DIFR.  A Z=ZD')
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])