%REFL3D trace des rayons PP et PS sur reflecteur plan 3D
% La source et le geophone sont au choix dans un puits dans le socle ou a
% l'interface socle-sediment. Les rayons sont prolonges dans les sediments
% si l'incidence est inferieure à l'incidence critique.
% voir : http://jmmeost.free.fr/Sismique/osis.html#VX
% fig1 = http://jmmeost.free.fr/Sismique/refl3da.gif
% fig2 = http://jmmeost.free.fr/Sismique/refl3db.gif

clear all, clf
rad = pi/180;
xx = [-500:100:500]; yy = xx; [XX YY] = meshgrid(xx,yy);

% vitesses sédiments et h surface
VP1 = 3000; VPS1 = 1.7; VS1 = VP1/VPS1; h1 = 500;
% vitesses socle
VP0 = 5000; VPS0 = 1.7; VS0 = VP0/VPS0;

% pendage, azimut du pendage par rapport à l'axe x en °, 
% profondeur à laquelle l'interface séd-socle intersecte x=y=0
ga1 = 10; az1 = -30; z1 = 0;  ga1 = ga1*rad; az1 = az1*rad;
% réflecteur dans socle (z0<0)
ga0 = 40; az0 = 30; z0 = -1000;  ga0 = ga0*rad; az0 = az0*rad;

% normales E et profondeur interfaces
e1x = cos(az1)*sin(ga1); e1y = sin(az1)*sin(ga1); e1z = cos(ga1); E1 = [e1x ; e1y; e1z]; ZZ1 = (z1*e1z-e1x*XX-e1y*YY)/e1z;
e0x = cos(az0)*sin(ga0); e0y = sin(az0)*sin(ga0); e0z = cos(ga0); E0 = [e0x ; e0y; e0z]; ZZ0 = (z0*e0z-e0x*XX-e0y*YY)/e0z; 

% source dans puits (xs=ys=0, zs<0, puits=1) ou en z1 (puits=0)
xs = -200; ys = 0;  zg = -100; puits = 0;
if puits == 0, zs = (z1*e1z-e1x*xs-e1y*ys)/e1z; end
XS = [xs ; ys; zs];

% geophones dans puits (xg=yg=0, zg<0, puitg=1) ou en z1 (puitg=0)
xg = 400; yg = 0; zg = -100; puitg = 0;
if puitg == 0, zg = (z1*e1z-e1x*xg-e1y*yg)/e1z;  end
XG = [xg ; yg; zg];

% projection de la source sur le reflecteur et source virtuelle
A = [e0x, e0y, e0z ; 0, e0z, -e0y ; e0z, 0, -e0x]; B = [z0*e0z ; e0z*ys-e0y*zs ; e0z*xs-e0x*zs];
XSR = A\B; xsr = XSR(1); ysr = XSR(2); zsr = XSR(3);
XSV = 2*XSR-XS; xsv = XSV(1); ysv = XSV(2); zsv = XSV(3);

% projection du geophone sur le reflecteur
A = [e0x, e0y, e0z ; 0, e0z, -e0y ; e0z, 0, -e0x]; B = [z0*e0z ; e0z*yg-e0y*zg ; e0z*xg-e0x*zg];
XGR = A\B; xgr = XGR(1); ygr = XGR(2); zgr = XGR(3);

%  point refl PP = intersection droite source virtuelle - geophone avec reflecteur
EGV = (XG-XSV)/norm(XG-XSV);
A = [e0x, e0y, e0z ; 0, EGV(3), -EGV(2) ; EGV(3), 0, -EGV(1)]; B = [z0*e0z ; EGV(3)*ysv-EGV(2)*zsv ; EGV(3)*xsv-EGV(1)*zsv];
XPP = A\B; xpp = XPP(1); ypp = XPP(2); zpp = XPP(3);
% ESPP = (XS-XPP)/norm(XS-XPP); snellPP = sin(acos(dot(ESPP,E0)))-sin(acos(dot(EGPP,E0)))
tpp = (norm(XPP-XS)+norm(XPP-XG))/VP0;

% point refl PS = temps minimum pour points sur segment point refl PP - proj geophone sur reflecteur
ESGR = (XGR-XPP)/norm(XGR-XPP);
bin = 1; nrs = fix(norm(XGR-XPP)/bin);
for ir = 1:nrs
XD = XPP+(ir-1)*bin*ESGR;
TDS(ir) = sqrt((XD(1)-xs).^2+(XD(2)-ys).^2+(XD(3)-zs).^2)/VP0+sqrt((XD(1)-xg).^2+(XD(2)-yg).^2+(XD(3)-zg).^2)/VS0;
end
tps = min(min(TDS)); ir = find(TDS == tps);
XPS = XPP+(ir-1)*bin*ESGR; xps = XPS(1); yps = XPS(2); zps = XPS(3);
% ESPS = (XS-XPS)/norm(XS-XPS); snellPS = sin(acos(dot(ESPS,E0)))-VPS0*sin(acos(dot(EGPS,E0)))

% prolongation P dans couche sédiment
EGPP = (XG-XPP)/norm(XG-XPP); cth0 = dot(EGPP,E1); 
aa = (VP0/VP1)^2-(1-cth0^2); 
if (aa >= 0 && puitg ~= 1)
EGP1 = (VP1/VP0)*(EGPP+E1*(sqrt(aa)-cth0)); cth1 = dot(EGP1,E1);
A = [0, 0, 1 ; 0, EGP1(3), -EGP1(2) ; EGP1(3), 0, -EGP1(1)]; B = [h1 ; EGP1(3)*yg-EGP1(2)*zg ; EGP1(3)*xg-EGP1(1)*zg];
XGP1 = A\B; xgp1 = XGP1(1); ygp1 = XGP1(2); zgp1 = XGP1(3);
else
xgp1 = xg; ygp1 = yg; zgp1 = zg; 
end
%snellGP1 = sin(acos(cth0))/VP0-sin(acos(cth1))/VP1

ESPP = (XS-XPP)/norm(XS-XPP); cth0 = dot(ESPP,E1);
aa = (VP0/VP1)^2-(1-cth0^2); 
if (aa >= 0 && puits ~= 1)
ESP1 = (VP1/VP0)*(ESPP+E1*(sqrt(aa)-cth0)); cth1 = dot(ESP1,E1);
A = [0, 0, 1 ; 0, ESP1(3), -ESP1(2) ; ESP1(3), 0, -ESP1(1)]; B = [h1 ; ESP1(3)*ys-ESP1(2)*zs ; ESP1(3)*xs-ESP1(1)*zs];
XSP1 = A\B; xsp1 = XSP1(1); ysp1 = XSP1(2); zsp1 = XSP1(3);
else
xsp1 = xs; ysp1 = ys; zsp1 = zs; 
end
%snellSP1 = sin(acos(cth0))/VP0-sin(acos(cth1))/VP1

ESPS = (XS-XPS)/norm(XS-XPS); cth0 = dot(ESPS,E1);
aa = (VP0/VP1)^2-(1-cth0^2); 
if (aa >= 0 && puits ~= 1)
ESPS1 = (VP1/VP0)*(ESPS+E1*(sqrt(aa)-cth0)); cth1 = dot(ESPS1,E1);
A = [0, 0, 1 ; 0, ESPS1(3), -ESPS1(2) ; ESPS1(3), 0, -ESPS1(1)]; B = [h1 ; ESPS1(3)*ys-ESPS1(2)*zs ; ESPS1(3)*xs-ESPS1(1)*zs];
XSPS1 = A\B ; xsps1 = XSPS1(1); ysps1 = XSPS1(2); zsps1 = XSPS1(3);
else
xsps1 = xs; ysps1 = ys; zsps1 = zs; 
end
%snellSPS1 = sin(acos(cth0))/VP0-sin(acos(cth1))/VP1

% prolongation S dans couche sédiment
EGPS = (XG-XPS)/norm(XG-XPS); cet0 = dot(EGPS,E1);
aa = (VS0/VS1)^2-(1-cet0^2);
if (aa >= 0 && puitg ~= 1)
EGS1 = (VS1/VS0)*(EGPS+E1*(sqrt(aa)-cet0)); cet1 = dot(EGS1,E1);
A = [0, 0, 1 ; 0, EGS1(3), -EGS1(2) ; EGS1(3), 0, -EGS1(1)]; B = [h1 ; EGS1(3)*yg-EGS1(2)*zg ; EGS1(3)*xg-EGS1(1)*zg];
XGS1 = A\B ; xgs1 = XGS1(1); ygs1 = XGS1(2); zgs1 = XGS1(3);
else
xgs1 = xg; ygs1 = yg; zgs1 = zg; 
end
%snellGS1 = sin(acos(cet0))/VS0-sin(acos(cet1))/VS1

figure(1), box on, grid on, hold on
xmax = 1000; ymax = 500;
% puits
plot3([0 0],[0 0],[h1 z1],'-y',[0 0],[0 0],[z1 z0],'-c','Linewidth',3)
% reflecteurs plan
surf(XX,YY,ZZ1), surf(XX,YY,ZZ0), 
% puits, z=0
% plot3([0 0],[0 0],[zg z0],'-k','Linewidth',3)
% source, geophone,source virtuelle,
plot3(xs,ys,zs,'ro',xg,yg,zg,'rs',xsv,ysv,zsv,'bo','Linewidth',2)
plot3([xsv xs],[ysv ys],[zsv zs],'b')
% rayons PP, PS, PS pour legende
plot3([xpp xs],[ypp ys],[zpp zs],'-r',[xps xs],[yps ys],[zps zs],'-m',[xg xgs1],[yg ygs1],[zg zgs1],'-g','Linewidth',2)
% rayons réfléchis PP
plot3([xpp xs],[ypp ys],[zpp zs],'-r',[xpp xg],[ypp yg],[zpp zg],'-r','Linewidth',2)
plot3([xg xgp1],[yg ygp1],[zg zgp1],'-r','Linewidth',2)
plot3([xs xsp1],[ys ysp1],[zs zsp1],'-r','Linewidth',2)
% rayons réfléchis PS
plot3([xps xs],[yps ys],[zps zs],'-m',[xps xg],[yps yg],[zps zg],'-g','Linewidth',2)
plot3([xg xgs1],[yg ygs1],[zg zgs1],'-g','Linewidth',2)
plot3([xs xsps1],[ys ysps1],[zs zsps1],'-m','Linewidth',2)
% plan propagation dans socle
plot3([xsv xs],[ysv ys],[zsv zs],'b',[xsv xpp],[ysv ypp],[zsv zpp],'b')
plot3([xsr xgr],[ysr ygr],[zsr zgr],'b',[xgr xg],[ygr yg],[zgr zg],'b',[xs xg],[ys yg],[zs zg],'b')

set(gca,'XLim',[-xmax xmax],'YLim',[-ymax ymax],'ZLim',[-1500 h1],'DataAspectRatio',[1 1 1])
xlabel('X (M)'),ylabel('Y (M)'),zlabel('Z (M)'), colormap('bone'), colorbar
legend(['VP1 = ' num2str(VP1/1000) ' KM/S , VS1 = ' num2str(VS1/1000,2) ' KM/S'],...
       ['VP0 = ' num2str(VP0/1000) ' KM/S , VS0 = ' num2str(VS0/1000,2) ' KM/S'],...
       ['PEND = ' num2str(ga1/rad) '° AZ = ' num2str(az1/rad) '°'],...
       ['PEND = ' num2str(ga0/rad) '° AZ = ' num2str(az0/rad) '°'],...
       ['XS = ' num2str(xs,3) ' YS = ' num2str(ys,3) ' ZS = '  num2str(zs,3) ' M'],...
       ['XG = ' num2str(xg,3) ' YG = ' num2str(yg,3) ' ZG = '  num2str(zg,3) ' M'],...
       'SOURCE VIRTUELLE','PLAN REFLEXION','RAYON PP','RAYON P DE PS','RAYON S DE PS')
title (['REFLEXION-TRANSMISSION 3D PP ET PS'])
view(30,20)


% temps trajet PP pour points diffractants autour point reflexion PP
bin = 2; lx = 10*bin; ly = 10*bin;
[XDP YDP] = meshgrid([xpp-lx:bin:xpp+lx],[ypp-ly:bin:ypp+ly]); ZDP = (z0*e0z-e0x*XDP-e0y*YDP)/e0z;
TDP = (sqrt((XDP-xs).^2+(YDP-ys).^2+(ZDP-zs).^2)+sqrt((XDP-xg).^2+(YDP-yg).^2+(ZDP-zg).^2))/VP0;

% temps trajet PS pour points diffractants autour point reflexion PS
[XDS YDS] = meshgrid([xps-lx:bin:xps+lx],[yps-ly:bin:yps+ly]); ZDS = (z0*e0z-e0x*XDS-e0y*YDS)/e0z;
TDS = sqrt((XDS-xs).^2+(YDS-ys).^2+(ZDS-zs).^2)/VP0+sqrt((XDS-xg).^2+(YDS-yg).^2+(ZDS-zg).^2)/VS0;

figure(2), clf
subplot(2,1,1),box on, grid on, hold on
surf(XDP-xpp,YDP-ypp,TDP-tpp)
set(gca,'XLim',[-lx lx],'YLim',[-ly ly])
xlabel('XDIFPP - XREFLPP (M)'), ylabel('YDIFPP - YREFLPP (M)'),  zlabel('S'), colorbar, colormap('hsv')
title ('TEMPS DIFFRACTION - REFLEXION PP')
view(60,80)

subplot(2,1,2),box on, grid on, hold on
surf(XDS-xps,YDS-yps,TDS-tps)
set(gca,'XLim',[-lx lx],'YLim',[-ly ly])
xlabel('XDIFPS - XREFLPS (M)'), ylabel('YDIFPS - YREFLPS (M)'), zlabel('S'), colorbar
title ('TEMPS DIFFRACTION - REFLEXION PS')
view(60,80)