%FRESNEL superpose ondes plane et sphérique harmonique dans le plan (x,z,t=0) et (x,z=0,t)
%Les interférences font apparaitre la zone de Fresnel
% voir : http://jmmeost.free.fr/Sismique/isis.html
% fig1 = http://jmmeost.free.fr/Sismique/fresnel1.gif

clear all

V0 = 1000;

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

[XT TT] = meshgrid(xx,tt);
[XZ ZZ] = meshgrid(xx,zz);

% fréquence fondamentale
f0 = 4; om = 2*pi*f0; k0 = om/V0;
% point diffractant
xd = xl/2; zd = 1000;

% zone Fresnel
xfz = sqrt(V0*zz/f0/2);
xft = sqrt(V0*zd/f0/2);


figure(1), clf, colormap('jet')

subplot(2,3,1), hold on
imagesc([0 xl],[0 zl],sin(-k0*(ZZ-zd)))
xlabel('X (M)'), ylabel('Z (M)'), title('REFLECTEUR A T=0')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,2), hold on
imagesc([0 xl],[0 zl],sin(k0*sqrt((XZ-xd).^2+(ZZ-zd).^2)))
xlabel('X (M)'), ylabel('Z (M)'), title('DIFFRACTION A T=0')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,3), hold on
imagesc([0 xl],[0 zl],sin(-k0*(ZZ-zd))+sin(k0*sqrt((XZ-xd).^2+(ZZ-zd).^2)))
plot(xfz+xd,-zz+zd,'k',-xfz+xd,-zz+zd,'k','Linewidth',2)
xlabel('X (M)'), ylabel('Z (M)'), title(['ZONE FRESNEL, \Lambda = ' num2str(V0/f0) ' M'])
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],sin(-om*TT+k0*zd))
xlabel('X (M)'), ylabel('T (S)'), title('REFLEXION A Z=0')
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

subplot(2,3,5), hold on
imagesc([0 xl],[0 tl],sin(-om*TT+k0*sqrt((XT-xd).^2+zd.^2)))
xlabel('X (M)'), ylabel('T (S)'), title(['DIFFRACTION A Z=0'])
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

subplot(2,3,6), hold on
imagesc([0 xl],[0 tl],sin(-om*TT+k0*zd)+sin(-om*TT+k0*sqrt((XT-xd).^2+zd.^2)))
plot([xft xft]+xd,[0 tl],'k',-[xft xft]+xd,[0 tl],'k','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'), title(['ZONE FRESNEL, F = ' num2str(f0) ' HZ'])
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

%CAS REFLECTEUR PENTE
% point diffractant
xd = 0; zd = 1000;
al = 30*pi/180;
% calcul largeur zone Fresnel en z=0
zza = [0:dz:2*zl]; xfza = sqrt(V0*zza/f0/2);
ixft = find((-zza*cos(al)+xfza*sin(al)+zd)<=0); xftd = zza(ixft(1))*sin(al)+xfza(ixft(1))*cos(al)+xd;
ixft = find((-zza*cos(al)-xfza*sin(al)+zd)<=0); xftg = zza(ixft(1))*sin(al)-xfza(ixft(1))*cos(al)+xd;


figure(2), clf, colormap('jet')

subplot(2,3,1), hold on
imagesc([0 xl],[0 zl],sin(k0*(sin(al)*(XZ-xd)-cos(al)*(ZZ-zd))))
plot([0 zl*tan(al)],[zl 0],'w','Linewidth',2)
xlabel('X (M)'), ylabel('Z (M)'), title(['REFLECTEUR A T=0, \alpha = ' num2str(al*180/pi) ' °'])
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,2), hold on
imagesc([0 xl],[0 zl],sin(k0*sqrt((XZ-xd).^2+(ZZ-zd).^2)))
plot([0 zl*tan(al)],[zl 0],'w','Linewidth',2)
xlabel('X (M)'), ylabel('Z (M)'), title('DIFFRACTION A T=0')
set(gca,'Xlim',[0 xl],'YLim',[0 zl],'YDir','reverse','DataAspectRatio',[1 1 1])

subplot(2,3,3), hold on
imagesc([0 xl],[0 zl],sin(k0*(sin(al)*(XZ-xd)-cos(al)*(ZZ-zd)))+sin(k0*sqrt((XZ-xd).^2+(ZZ-zd).^2)))
plot([0 zl*tan(al)],[zl 0],'w','Linewidth',2)
plot(zza*sin(al)+xfza*cos(al)+xd,-zza*cos(al)+xfza*sin(al)+zd,'k',zza*sin(al)-xfza*cos(al)+xd,-zza*cos(al)-xfza*sin(al)+zd,'k','Linewidth',2)
xlabel('X (M)'), ylabel('Z (M)'), title(['ZONE FRESNEL, \Lambda = ' num2str(V0/f0) ' M'])
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],sin(-om*TT+k0*(sin(al)*(XZ-xd)+cos(al)*zd)))
xlabel('X (M)'), ylabel('T (S)'), title('REFLEXION = ONDE PLANE A Z=0')
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

subplot(2,3,5), hold on
imagesc([0 xl],[0 tl],sin(-om*TT+k0*sqrt((XT-xd).^2+zd.^2)))
xlabel('X (M)'), ylabel('T (S)'), title(['DIFFRACTION A Z=0'])
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])

subplot(2,3,6), hold on
imagesc([0 xl],[0 tl],sin(-om*TT+k0*(sin(al)*(XZ-xd)+cos(al)*zd))+sin(-om*TT+k0*sqrt((XT-xd).^2+zd.^2)))
plot([xftg xftg],[0 tl],'k',[xftd xftd],[0 tl],'k','Linewidth',2)
xlabel('X (M)'), ylabel('T (S)'),  title(['ZONE FRESNEL, F = ' num2str(f0) ' HZ'])
set(gca,'Xlim',[0 xl],'YLim',[0 tl],'YDir','reverse','DataAspectRatio',[1 .001 1])