%DIF15 effectue migration par resolution par differences finies implicites
%de l'équation paraxiale 15° pour les ondes montantes.
%voir http://jmmeost.free.fr/Sismique/isis.html#DF

clear all

V0= 2000; z0 = 400;
dd = 10; zmax = 1000; xmax = 1000; xx = 00:dd:xmax; nx = xmax/dd+1; zz = 00:dd:zmax; nz = zmax/dd+1;
f0 = V0/2/zmax; nf = 20;
theta = 30; % angle pour cond. limite

figure(1), clf, colormap('copper'), ifig = -1;

for theta = 15:10:45 % reflexions planes
% for z0 = 400:200:600 % impulsions
t0 = z0/V0;

pf = zeros(nz,nx);

for fq = f0:f0:nf*f0
om = 2*pi*fq; edd = exp(-i*om/V0*dd);
th = theta*pi/180; kx = om*sin(th)/V0; edx = exp(i*kx*dd);

%differences finies 15° implicite
aa = V0/4/i/om/dd; 
bb = 1/aa+2; cc = 1/aa-2; % coef. tridiagonaux divises par aa
ee(2) = 1/edx; % cond. limite pour pendage theta
for ix = 2:nx
ee(ix+1) = -1/(ee(ix)-bb); % p(ix-1) = ee(ix)p(ix)+ff(ix)
end
pp = zeros(nz,nx); 
%pp(1,:) = exp(i*(kx*xx+om*t0)); % p(x,t) = reflexions planes
ix0 = fix(nx/2); pp(:,ix0) = exp(i*om*(t0-zz/V0)); % p(x,t) = impulsions
for itr = 1:5, pp(:,ix0+itr) = pp(:,ix0)/itr/2;  pp(:,ix0-itr) = pp(:,ix0)/itr/2; end % impulsion elargie
qq(1) = 0; qq(nx) = 0;
for iz = 1:nz-1
ff(2) = 0; 
for ix = 2:nx-1
qq(ix) = pp(iz,ix-1)+cc*pp(iz,ix)+pp(iz,ix+1); 
ff(ix+1) = (qq(ix)+ff(ix))*ee(ix+1);
end
pp(iz+1,nx) = ff(nx)/(1-ee(nx)/edx); %cond. limite
for ix = nx:-1:3
pp(iz+1,ix-1) = ee(ix)*pp(iz+1,ix)+ff(ix);
end
pp(iz+1,1) = pp(iz+1,2)/edx; % cond. limite
pp(iz+1,:) = edd*pp(iz+1,:); %prolong. vertical
end
pf = pf + pp;
end

ifig = ifig+1;
%reflecteurs plans
tgth = tan(th); zm = (xmax*tgth+z0)/(1+tgth^2); xm = xmax-zm*tgth; %sol. geometrique
subplot(2,2,1+ifig), hold on
imagesc([0 xmax],[0 zmax],real(pf/nf),[-0.2 1]), colorbar 
plot([0 xm],[z0 zm],'r','Linewidth',2)
xlim([0 xmax]), ylim([0 zmax]), set(gca, 'YDir', 'reverse'), axis('square')
title(['DIF. FINIES 15° , THETA = ' num2str(theta) '°']), xlabel('X (M)'), ylabel('Z (M)')

% reponses impulsionnelles
% zz0 = sqrt(z0^2-(xx-xmax/2).^2); zz0(find(real(zz0) == 0)) = NaN; %sol. geometrique
% subplot(2,2,1+ifig), hold on 
% imagesc([0 xmax],[0 zmax],real(pf)), colorbar 
% plot(xx,zz0,'r','Linewidth',2)
% xlim([0 xmax]), ylim([0 zmax]), set(gca, 'YDir', 'reverse'), axis('square')
% title(['DIF. FINIES 15° , V = ' num2str(V0) ' M/S']), xlabel('X (M)'), ylabel('Z (M)')
% subplot(2,2,3+ifig), hold on 
% imagesc([0 xmax],[0 zmax],real(pf),[0.2 2]), colorbar 
% plot(xx,zz0,'r','Linewidth',2)
% xlim([0 xmax]), ylim([0 zmax]), set(gca, 'YDir', 'reverse'), axis('square')
% title(['DIF. FINIES 15° , V = ' num2str(V0) ' M/S']), xlabel('X (M)'), ylabel('Z (M)')

end