%MIGSDC migre à déport nul par superposition de demi-cercles dans le plan xz
% et par point de tir 
% voir : http://jmmeost.free.fr/Sismique/isis.html
% fig1 = http://jmmeost.free.fr/Sismique/migsdc1.gif

clear all

V0 = 1500;

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

%type modélisation 1 = points diffractants, 2 = réflexion
mod = 2;
% position des points diffractants ou reflecteur
%xd = [500]; nd = length(xd); zd = [300]; %fig1
%xd = [200:50:800]; nd = length(xd); zd = 300*ones(1,nd); %fig2
%xd = [0:30:300]; nd = length(xd); al = 30*pi/180; zd = 200+tan(al)*(xd-xd(1)); %fig3
%xd = [0:30:1000]; nd = length(xd); zd = 300-200*cos(2*pi*(xd-xd(1))/1000); %fig4
al = -10*pi/180; xd = [0 xl]; zd = 300+[0 xl*tan(al)];

%type migration 1 = déport nul, 2 = point tir
mig = 2;
%choix trace migrée pour subplot4
ixmig = nx/2+1; 
if mod ==1, tix = (sqrt(xd(nd)^2+zd(nd)^2)+sqrt((xx(ixmig)-xd(nd))^2+zd(nd)^2))/V0; zix = sqrt((V0*tix/2)^2-(xd(nd)/2)^2); end

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

pxt = zeros(nt,nx);
% MODELISATION
if mod ==1 
% DIFFRACTIONS
% onde diffractéee en xt
for id = 1:nd
if mig ==1, tmod = 2*(sqrt((xx-xd(id)).^2+zd(id)^2))/V0; 
else tmod = (sqrt(xd(id)^2+zd(id)^2)+sqrt((xx-xd(id)).^2+zd(id)^2))/V0; end
for ix = 1:nx
omt = om*(tt-tmod(ix)); 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
end

else
%REFLEXION
if mig ==1, tmod = 2*(xx*sin(al)+zd(1)*cos(al))/V0; dv = zd(2)*cos(al);xv = xl-dv*sin(al); zv = dv*cos(al);
else, dv = 2*zd(1)*cos(al); xv = -dv*sin(al); zv = dv*cos(al); tmod = sqrt((xx-xv).^2+zv^2)/V0; end
for ix = 1:nx
omt = om*(tt-tmod(ix)); 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
end

% MIGRATION A V0
VM = V0;
pxz = zeros(nz,nx); pxz1 = zeros(nz,nx); pxtm = zeros(nt,nx); 
[X2 Z2] = meshgrid(xx.^2,zz.^2);
% temps diffractés
txz = sqrt(X2+Z2)/VM;
% demi-cercles pour chaque trace et sommation
for ix = 1:nx
if mig ==1, 
% temps diffractés symétriques par rapport à trace courante
tm = 2*txz;
itm = round(tm/dt); itm(find(itm>nt-1)) = nt-1; itm(find(itm<1)) = 1;
itm = [itm(:,ix:-1:1) itm(:,2:nx-ix+1)]; 
end
if mig == 2
% temps source-diffracteur + temps diffracté
tm = txz(:,[1:ix ix+1:nx])+txz(:,[ix:-1:1 2:nx-ix+1]);
itm = round(tm/dt); itm(find(itm>nt-1)) = nt-1; itm(find(itm<1)) = 1;
end
% copie de la trace à migrer dans matrice pxtm
pxtm = pxt(:,ix)*ones(1,nx);
% migration
pxz = pxz+pxtm(itm);
if ix == ixmig, pxz1 = pxtm(itm); end
end
    
figure(1), clf, colormap('copper')
subplot(2,2,1), hold on
imagesc([0 xl],[0 zl],zeros(nz,nx))
if mod == 1, plot([0 xd(nd)],[0 zd(nd)],'b',[xd(nd)/2 xd(nd)],[0 zd(nd)],'r',[xd(nd) xd(nd)],[0 zd(nd)],'r',[3*xd(nd)/2 xd(nd)],[0 zd(nd)],'r',[xd(nd) xl],[zd(nd) 0],'r','Linewidth',2)
plot(xd,zd,'ro','MarkerFaceColor','r'), legend('RAYON SOURCE-POINT DIFFRACTANT','RAYONS POINT DIFFRACTANT-RECEPTEURS','Location','South')
else plot(xd,zd,'y',[0 xv],[0 zv],'b',[xl xv],[0 zv],'r','Linewidth',2), end
xlabel('X (M)'), ylabel('Z (M)'), 
if mig ==1, title('POINTS DIFFRACTANTS / REFLECTEUR A T=0'), else, 
if mod == 1, title('POINTS DIFFRACTANTS ET RAYONS'), else, title('REFLECTEUR ET TRIANGLE S-SV-G'), end
end
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,2,2), hold on
imagesc([0 xl],[0 zl],pxz)
if mod == 1, plot(xd,zd,'ro'), else plot(xd,zd,'y',[0 xv],[0 zv],'b',[xl xv],[0 zv],'r','Linewidth',2), end
xlabel('X (M)'), ylabel('Z (M)'), 
if mig ==1, title('MIGRATION PAR SUPERPOSITION DE 1/2 CERCLES'), else, title('MIGRATION DU POINT DE TIR'), end
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,2,3), hold on
imagesc([0 xl],[0 tl],pxt)
%plot([0 xl],2*zd*cos(al)/V0,'r','Linewidth',2)
plot(xx,tmod,'r','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), 
if mig ==1, title(['DIFFRACTIONS / REFLEXION A Z=0 (F = ' num2str(f0) ' - ' num2str(nf*f0) ' HZ, V = ' num2str(V0) ' M/S)']), 
else, title(['POINT DE TIR (F = ' num2str(f0) ' - ' num2str(nf*f0) ' HZ, V = ' num2str(V0) ' M/S)']), end
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .002 1])

subplot(2,2,4), hold on
imagesc([0 xl],[0 zl],pxz1)
if mod == 1, plot(xd,zd,'ro','MarkerFaceColor','r'), 
if nd ==1, plot([0 xd(nd)],[0 zd(nd)],'b',[xd(nd) xd(nd)],[0 zd(nd)],'r',[dx dx],[0 zd(nd)],'b',[0 xd(nd)],[zd(nd) 0],'r',[0 xx(ixmig)/2],[0 zix],'b',[xx(ixmig) xx(ixmig)/2],[0 zix],'r','Linewidth',2), end
else plot(xd,zd,'y',[0 xv],[0 zv],'b',[xl xv],[0 zv],'r','Linewidth',2), end
xlabel('X (M)'), ylabel('Z (M)'), title(['MIGRATION DE LA TRACE A X = ' num2str((ixmig-1)*dx) ' M'])
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

% %MIGRATION A VM = 1.1*V0 et 0.9*V0
% 
% figure(2), clf, colormap('jet')
% subplot(2,2,1), hold on
% imagesc([0 xl],[0 tl],pxt), colorbar
% xlabel('X (M)'), ylabel('T (S)'), title(['DIFFRACTIONS A Z=0 (F = ' num2str(f0) ' - ' num2str(nf*f0) ' HZ, V = ' num2str(V0) ' M/S)'])
% set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])
% 
% subplot(2,2,3), hold on
% imagesc([0 xl],[0 zl],pxz), colorbar
% plot(xd,zd,'ro','MarkerFaceColor','r')
% xlabel('X (M)'), ylabel('Z (M)'), title(['MIGRATION A VM = ' num2str(V0) ' M/S'])
% set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])
% 
% VM = 1.1*V0;
% 
% pxz = zeros(nz,nx); pxtm = zeros(nt,nx); 
% % temps diffractés
% txz = 2*sqrt(X2+Z2)/VM;
% itxz = round(txz/dt); itxz(find(itxz>nt-1)) = nt-1; itxz(find(itxz<1)) = 1;
% % demi-cercles pour chaque trace et sommation
% for ix = 1:nx
% % copie de la trace à migrer dans matrice pxtm
% for ixm = 1:nx, pxtm(:,ixm) = pxt(:,ix); end
% % temps diffractés symétriques par rapport à trace courante
% itxzm = [itxz(:,ix:-1:1) itxz(:,2:nx-ix+1)];
% pxz = pxz+pxtm(itxzm);
% end
% pxz = pxz/nx;
% 
% subplot(2,2,2), hold on
% imagesc([0 xl],[0 zl],pxz), colorbar
% plot(xd,zd,'ro','MarkerFaceColor','r')
% xlabel('X (M)'), ylabel('Z (M)'), title(['MIGRATION A VM = ' num2str(VM) ' M/S'])
% set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])
% 
% VM = 0.9*V0;
% 
% pxz = zeros(nz,nx); pxtm = zeros(nt,nx); 
% % temps diffractés
% txz = 2*sqrt(X2+Z2)/VM;
% itxz = round(txz/dt); itxz(find(itxz>nt-1)) = nt-1; itxz(find(itxz<1)) = 1;
% % demi-cercles pour chaque trace et sommation
% for ix = 1:nx
% % copie de la trace à migrer dans matrice pxtm
% for ixm = 1:nx, pxtm(:,ixm) = pxt(:,ix); end
% % temps diffractés symétriques par rapport à trace courante
% itxzm = [itxz(:,ix:-1:1) itxz(:,2:nx-ix+1)];
% pxz = pxz+pxtm(itxzm);
% end
% pxz = pxz/nx;
% 
% subplot(2,2,4), hold on
% imagesc([0 xl],[0 zl],pxz), colorbar
% plot(xd,zd,'ro','MarkerFaceColor','r')
% xlabel('X (M)'), ylabel('Z (M)'), title(['MIGRATION A VM = ' num2str(VM) ' M/S'])
% set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])