%ondePS dessine les 3 composantes du déplacement des particules le long des
% axes x,y,z pour des ondes planes harmoniques P, SV et SH de fréquence ff, 
% ayant une amplitude maximum du déplacement des particules U,
% se propageant dans une direction e définie par les angles thz et azx
%
% Documenter l'effet des changements de la direction e sur :
% - les longueurs d'onde apparentes le long des axes x,y,z
% - les amplitudes relatives des composantes du déplacement des particules
% Dans chaque cas, vérifier que les valeurs mesurées sur les figures 
% sont égales aux valeurs théoriques

clear all

% PARAMETRES **************************************************************
VP = 5000; % vitesse de propagation des ondes P (m/s)
nu = 0.25; % coefficient Poisson pour calcul de VS
ff = 10; om = 2*pi*ff; TT = 1/ff; % fréquence (Hz), pulsation, période (s)

U = .01; % amplitude déplacement particule (m)

thz = 30*pi/180; % angle (axe z - direction propagation)
azx = 0*pi/180; % azimut (axe x - plan de propagation) 

dx = 10; xl = 1000; % echantillonnage et distance max pour plot (m)

tt = TT/4; % temps auquel l'onde est dessinée en fraction de période

% dessine onde P si OP = 'O', onde SV si OSV = 'O', onde SH si OSH = 'O'
OP = 'O'; OSV = 'O'; OSH = 'O';
%**************************************************************************


VS = sqrt((1-2*nu)/2/(1-nu)*VP^2); % vitesse ondes S
Phi = U*VP/om; Psi = U*VS/om; % amplitudes potentiels phi onde P,  psi onde SV ( m^2)

ex = sin(thz)*cos(azx); ey = sin(thz)*sin(azx); ez = cos(thz); % composantes direction propagation

xx = 0:dx:xl; yy = xx; zz = xx; % distances le long des axes

kp = om/VP; kpx = kp*ex; kpy = kp*ey; kpz = kp*ez; %nombre d'onde pour onde P
ks = om/VS; ksx = ks*ex; ksy = ks*ey; ksz = ks*ez; %nombre d'onde pour onde S
lap = 2*pi/kp; lapx = 2*pi/kpx; lapy = 2*pi/kpy; lapz = 2*pi/kpz; %longueur onde P
las = 2*pi/ks; lasx = 2*pi/ksx; lasy = 2*pi/ksy; lasz = 2*pi/ksz; %longueur onde P

if OP == 'O' % ONDE PLANE P ************************************************************
% composantes du déplacement des particules (ux,uy,uz)  le long de x
phi = Phi*exp(1i*(kpx*xx-om*tt)); % phi(x, y = 0, z = 0, t = 0)
uxx = real(1i*kpx*phi); uyx = real(1i*kpy*phi); uzx = real(1i*kpz*phi); % u = grad(phi)
% composantes du déplacement des particules (ux,uy,uz) le long de y
phi = Phi*exp(1i*(kpy*yy-om*tt)); % phi(x = 0, y, z = 0, t = 0)
uxy = real(1i*kpx*phi); uyy = real(1i*kpy*phi); uzy = real(1i*kpz*phi); % u = grad(phi)
% composantes du déplacement des particules (ux,uy,uz) le long de z
phi = Phi*exp(1i*(kpz*zz-om*tt)); % phi(x = 0, y = 0, z , t = 0)
uxz = real(1i*kpx*phi); uyz = real(1i*kpy*phi); uzz = real(1i*kpz*phi); % u = grad(phi)

figure(1), clf, 

subplot(2,2,1), hold on, box on, grid on
plot3([0 ex],[0 ey],[0 ez],'r',[0 ex],[0 ey],[0 0],'r--',[0 ex],[0 0],[0 0],'r--',[0 0],[0 ey],[0 0],'r--',[0 0],[0 0],[0 ez],'r--','Linewidth',3)
plot3([0 cos(azx)],[0 sin(azx)],[0 0],'k')
y1 = sqrt(1-(0:.1:1).^2);  plot3((0:.1:1) , y1, 0*(0:.1:1), 'k') % cercle unité dans plan z = 0
th1 = (0:90)*pi/180; plot3(sin(th1)*cos(azx),sin(th1)*sin(azx),cos(th1),'k') % cercle unité dans plan propagation
legend(['THZ = ' num2str(thz/pi*180) ' °'],['AZX = ' num2str(azx/pi*180) ' °'])
set(gca,'XLim',[0 1],'YLim',[0 1],'ZLim',[0 1],'DataAspectRatio',[1,1,1]), view([120,30])
xlabel('X'), ylabel('Y'), zlabel('Z')
title(['ONDE PLANE P, VP = ' num2str(VP) ' M/S , F = ' num2str(ff), ' Hz'])

subplot(2,2,2), hold on, box on, grid on
plot(zz,uxz,'b',zz,uyz,'g',zz,uzz,'r','Linewidth',2)
legend('UX', 'UY', 'UZ')
set(gca,'XLim',[0 xl],'YLim',Phi*kp*[-1 1]) % échelle déplacement ajustée à valeur max
xlabel('Z (M)'), ylabel('DEPLACEMENT (M)'), 
title(['U( X = 0, Y = 0, Z , T = ' num2str(tt) ' S)'])

subplot(2,2,3), hold on, box on, grid on
plot(xx,uxx,'b',xx,uyx,'g',xx,uzx,'r','Linewidth',2)
legend('UX', 'UY', 'UZ')
set(gca,'XLim',[0 xl],'YLim',Phi*kp*[-1 1])
xlabel('X (M)'), ylabel('DEPLACEMENT (M)'), 
title(['U( X , Y = 0, Z = 0, T = ' num2str(tt) ' S)'])

subplot(2,2,4), hold on, box on, grid on
plot(yy,uxy,'b',yy,uyy,'g',yy,uzy,'r','Linewidth',2)
legend('UX', 'UY', 'UZ')
set(gca,'XLim',[0 xl],'YLim',Phi*kp*[-1 1])
xlabel('Y (M)'), ylabel('DEPLACEMENT (M)'), 
title(['U( X = 0, Y , Z = 0, T = ' num2str(tt) ' S)'])
end %**************************************************************************


if OSV == 'O' % ONDE PLANE SV ***********************************************************
% composantes du déplacement des particules (ux,uy,uz)  le long de x
psi = Psi*exp(1i*(ksx*xx-om*tt)); % psi(x, y = 0, z = 0, t = 0)
psix = sin(azx)*psi; psiy = -cos(azx)*psi; % psi orthogonal au plan de propagation
uxx = real(-1i*ksz*psiy); uyx = real(1i*ksz*psix); uzx = real(1i*(ksx*psiy-ksy*psix)); % u = rot(psi)
% composantes du déplacement des particules (ux,uy,uz) le long de y
psi = Psi*exp(1i*(ksy*yy-om*tt)); % psi(x = 0, y, z = 0, t = 0)
psix = sin(azx)*psi; psiy = -cos(azx)*psi; 
uxy = real(-1i*ksz*psiy); uyy = real(1i*ksz*psix); uzy = real(1i*(ksx*psiy-ksy*psix));
% composantes du déplacement des particules (ux,uy,uz) le long de z
psi = Psi*exp(1i*(ksz*zz-om*tt)); % psi(x = 0, y = 0, z, t = 0)
psix = sin(azx)*psi; psiy = -cos(azx)*psi; 
uxz = real(-1i*ksz*psiy); uyz = real(1i*ksz*psix); uzz = real(1i*(ksx*psiy-ksy*psix));

figure(2), clf, 

subplot(2,2,1), hold on, box on, grid on
plot3([0 cos(azx)],[0 sin(azx)],[0 0],'k')
y1 = sqrt(1-(0:.1:1).^2);  plot3((0:.1:1) , y1, 0*(0:.1:1), 'k') % cercle unité dans plan z = 0
th1 = (0:10:90)*pi/180; plot3(sin(th1)*cos(azx),sin(th1)*sin(azx),cos(th1),'k') % cercle unité dans plan propagation
plot3([0 ex],[0 ey],[0 ez],'r',[0 ex],[0 ey],[0 0],'r--',[0 ex],[0 0],[0 0],'r--',[0 0],[0 ey],[0 0],'r--',[0 0],[0 0],[0 ez],'r--','Linewidth',3)
legend(['THZ = ' num2str(thz/pi*180) ' °'],['AZX = ' num2str(azx/pi*180) ' °'])
set(gca,'XLim',[0 1],'YLim',[0 1],'ZLim',[0 1],'DataAspectRatio',[1,1,1]), view([120,30])
xlabel('X'), ylabel('Y'), zlabel('Z')
title(['ONDE PLANE SV , VS = ' num2str(VS,4) ' M/S , F = ' num2str(ff), ' Hz'])

subplot(2,2,2), hold on, box on, grid on
plot(zz,uxz,'b',zz,uyz,'g',zz,uzz,'r','Linewidth',2)
legend('UX', 'UY', 'UZ')
set(gca,'XLim',[0 xl],'YLim',Psi*ks*[-1 1])  % échelle déplacement ajustée à valeur max
xlabel('Z (M)'), ylabel('DEPLACEMENT (M)'), 
title(['U( X = 0, Y = 0, Z , T = ' num2str(tt) ' S)'])

subplot(2,2,3), hold on, box on, grid on
plot(xx,uxx,'b',xx,uyx,'g',xx,uzx,'r','Linewidth',2)
legend('UX', 'UY', 'UZ')
set(gca,'XLim',[0 xl],'YLim',Psi*ks*[-1 1])
xlabel('X (M)'), ylabel(' DEPLACEMENT (M)'), 
title(['U( X , Y = 0, Z = 0, T = ' num2str(tt) ' S)'])

subplot(2,2,4), hold on, box on, grid on
plot(yy,uxy,'b',yy,uyy,'g',yy,uzy,'r','Linewidth',2)
legend('UX', 'UY', 'UZ')
set(gca,'XLim',[0 xl],'YLim',Psi*ks*[-1 1])
xlabel('Y (M)'), ylabel('DEPLACEMENT (M)'), 
title(['U( X = 0, Y , Z = 0, T = ' num2str(tt) ' S)'])
end %**************************************************************************



if OSH == 'O' %ONDE PLANE SH ***********************************************************
% composantes du déplacement des particules (ux,uy,uz)  le long de x
ush = U*exp(1i*(ksx*xx-om*tt)); % ush(x, y = 0, z = 0, t = 0)
ushx = sin(azx)*ush; ushy = -cos(azx)*ush; % ush orthogonal au plan de propagation
uxx = real(ushx); uyx = real(ushy); uzx = 0;
% composantes du déplacement des particules (ux,uy,uz) le long de y
ush = U*exp(1i*(ksy*yy-om*tt)); % ush(x = 0, y, z = 0, t = 0)
ushx = sin(azx)*ush; ushy = -cos(azx)*ush; % ush orthogonal au plan de propagation
uxy = real(ushx); uyy = real(ushy); uzy = 0;
% composantes du déplacement des particules (ux,uy,uz) le long de z
ush = U*exp(1i*(ksz*zz-om*tt)); % ush(x = 0, y = 0, z, t = 0)
ushx = sin(azx)*ush; ushy = -cos(azx)*ush; % ush orthogonal au plan de propagation
uxz = real(ushx); uyz = real(ushy); uzz = 0;

figure(3), clf, 

subplot(2,2,1), hold on, box on, grid on
plot3([0 ex],[0 ey],[0 ez],'r',[0 ex],[0 ey],[0 0],'r--',[0 ex],[0 0],[0 0],'r--',[0 0],[0 ey],[0 0],'r--',[0 0],[0 0],[0 ez],'r--','Linewidth',3)
plot3([0 cos(azx)],[0 sin(azx)],[0 0],'k')
y1 = sqrt(1-(0:.1:1).^2);  plot3((0:.1:1) , y1, 0*(0:.1:1), 'k') % cercle unité dans plan z = 0
th1 = (0:10:90)*pi/180; plot3(sin(th1)*cos(azx),sin(th1)*sin(azx),cos(th1),'k') % cercle unité dans plan propagation
legend(['THZ = ' num2str(thz/pi*180) ' °'],['AZX = ' num2str(azx/pi*180) ' °'])
set(gca,'XLim',[0 1],'YLim',[0 1],'ZLim',[0 1],'DataAspectRatio',[1,1,1]), view([120,30])
xlabel('X'), ylabel('Y'), zlabel('Z')
title(['ONDE PLANE SH , VS = ' num2str(VS,4) ' M/S , F = ' num2str(ff), ' Hz'])

subplot(2,2,2), hold on, box on, grid on
plot(zz,uxz,'b',zz,uyz,'g',zz,uzz,'r','Linewidth',2)
legend('UX', 'UY', 'UZ')
set(gca,'XLim',[0 xl],'YLim',U*[-1 1])  % échelle déplacement ajustée à valeur max
xlabel('Z (M)'), ylabel('DEPLACEMENT (M)'), 
title(['U( X = 0, Y = 0, Z , T = ' num2str(tt) ' S)'])

subplot(2,2,3), hold on, box on, grid on
plot(xx,uxx,'b',xx,uyx,'g',xx,uzx,'r','Linewidth',2)
legend('UX', 'UY', 'UZ')
set(gca,'XLim',[0 xl],'YLim',U*[-1 1])
xlabel('X (M)'), ylabel(' DEPLACEMENT (M)'), 
title(['U( X , Y = 0, Z = 0, T = ' num2str(tt) ' S)'])

subplot(2,2,4), hold on, box on, grid on
plot(yy,uxy,'b',yy,uyy,'g',yy,uzy,'r','Linewidth',2)
legend('UX', 'UY', 'UZ')
set(gca,'XLim',[0 xl],'YLim',U*[-1 1])
xlabel('Y (M)'), ylabel('DEPLACEMENT (M)'), 
title(['U( X = 0, Y , Z = 0, T = ' num2str(tt) ' S)'])
end %**************************************************************************