%FORCE calcule les déplacements et polarisation P et S en champ proche et lointain 
%pour une force verticale ponctuelle harmonique 
% voir : http://jmmeost.free.fr/Sismique/osis.html#FP
% fig1  = http://jmmeost.free.fr/Sismique/force.gif
% fig2  = http://jmmeost.free.fr/Sismique/forcep.gif
% fig3  = http://jmmeost.free.fr/Sismique/forces.gif

clear all


VP = 3000; VS = VP/1.7; VP2 = VP^2; VS2 = VS^2;
F = 1; rho = 2600; UP = F/4/pi/rho/VP2; US = UP*VP2/VS2;
ff = 100; om = 2*pi*ff; phi =0;
kp = om/VP; lp = 2*pi/kp;
ks = om/VS; ls = 2*pi/ks; 

% ondes P et S pour x < lambdaP/2
xmax = lp/2; ymax = lp/2;
xx = -xmax:.2:xmax; yy = -ymax:.2:ymax; [X,Y] = meshgrid(xx,yy);

rr = sqrt(X.^2+Y.^2); 
theta = atan(X./Y);  ith = find(Y < 0 & X >= 0); theta(ith) = theta(ith)+pi; ith = find(Y < 0 & X < 0); theta(ith) = theta(ith)-pi; 
ir0 = find(rr < .2);  rr(ir0) = NaN; theta(ir0) = NaN;

kpr = kp*rr; eikpr = cos(kpr-phi)./rr; ieikpr = -sin(kpr-phi)./rr; pcp = ((ieikpr-eikpr./kpr)./kpr+1./rr.^3);
ksr = ks*rr; eiksr = cos(ksr-phi)./rr; ieiksr = -sin(ksr-phi)./rr; scp = ((ieiksr-eiksr./ksr)./ksr+1./rr.^3);

% déplacements P
upcl = UP*cos(theta).*eikpr; upcpr = UP*2*cos(theta).*pcp; upcpt = UP*sin(theta).*pcp; 
upn = sqrt((upcl+upcpr).^2+upcpt.^2); mp = max(max(upn));
% déplacements S
uscl = -US*sin(theta).*eiksr; uscpr = -US*2*cos(theta).*scp; uscpt = -US*sin(theta).*scp;
usn = sqrt(uscpr.^2+(uscl+uscpt).^2); ms = max(max(usn));
% polarisation P
rr = lp*[.08 .16 .25 .35 .49]; 
kpr = kp*rr; eikpr = cos(kpr-phi)./rr; ieikpr = -sin(kpr-phi)./rr; pcp = ((ieikpr-eikpr./kpr)./kpr+1./rr.^3);

ith = 0;
for theta = [0:pi/8:2*pi]
ith = ith+1;
xqp(ith,:) = rr*sin(theta); yqp(ith,:) = rr*cos(theta);
upcl = UP*cos(theta).*eikpr; upcpr = UP*2*cos(theta).*pcp; upcpt = UP*sin(theta).*pcp; upq = sqrt((upcl+upcpr).^2+upcpt.^2);
px(ith,:) = ((upcpr+upcl).*sin(theta)+upcpt.*cos(theta))./upq;
pz(ith,:) = ((upcpr+upcl).*cos(theta)-upcpt.*sin(theta))./upq;
end
% polarisation S
rr = ls*[.2 .4 .6 .9]; 
ksr = ks*rr; eiksr = cos(ksr-phi)./rr; ieiksr = -sin(ksr-phi)./rr; scp = ((ieiksr-eiksr./ksr)./ksr+1./rr.^3);

ith = 0;
for theta = [0:pi/8:2*pi]
ith = ith+1;
xqs(ith,:) = rr*sin(theta); yqs(ith,:) = rr*cos(theta);
uscl = -US*sin(theta).*eiksr; uscpr = -US*2*cos(theta).*scp; uscpt =-US*sin(theta).*scp; usq = sqrt(uscpr.^2+(uscl+uscpt).^2);
sx(ith,:) = (uscpr.*sin(theta)+(uscl+uscpt).*cos(theta))./usq;
sz(ith,:) = (uscpr.*cos(theta)-(uscl+uscpt).*sin(theta))./usq;
end

figure(1), clf
subplot(2,2,1)
pcolor(xx/lp,yy/lp,upn), caxis([0,mp/10000]), shading('interp'), colorbar, hold on
quiver(xqp/lp,yqp/lp,px,pz,.25,'w','Linewidth',1.2)
quiver(0,0,0,cos(phi),.15,'m','Linewidth',2)
set(gca,'XLim',xmax*[-1 1]/lp,'YLim',ymax*[-1 1]/lp,'DataAspectRatio',[1,1,1])
xlabel('X / {\Lambda}'), ylabel('Z / {\Lambda}'), 


subplot(2,2,2)
pcolor(xx/ls,yy/ls,usn), caxis([0,ms/9400]),  shading('interp'), colorbar, hold on
quiver(xqs/ls,yqs/ls,sx,sz,.25,'w','Linewidth',1.2)
quiver(0,0,0,cos(phi),.25,'m','Linewidth',2)
set(gca,'XLim',xmax*[-1 1]/ls,'YLim',ymax*[-1 1]/ls,'DataAspectRatio',[1,1,1])
xlabel('X / {\Lambda}'), ylabel('Z / {\Lambda}')

clear xqp yqp xqs yqs px pz sx sz

% ondes P et S pour x < 2*lambdaP
xmax = 2*lp; ymax = 2*lp;
xx = -xmax:.2:xmax; yy = -ymax:.2:ymax; [X,Y] = meshgrid(xx,yy);

rr = sqrt(X.^2+Y.^2); theta = atan(X./Y); ir0 = find(rr < 1);  rr(ir0) = NaN; theta(ir0) = NaN;

kpr = kp*rr; eikpr = cos(kpr-phi)./rr; ieikpr = -sin(kpr-phi)./rr; pcp = ((ieikpr-eikpr./kpr)./kpr+1./rr.^3);
ksr = ks*rr; eiksr = cos(ksr-phi)./rr; ieiksr = -sin(ksr-phi)./rr; scp = ((ieiksr-eiksr./ksr)./ksr+1./rr.^3);
% déplacements P
upcl = UP*cos(theta).*eikpr; upcpr = UP*2*cos(theta).*pcp; upcpt = UP*sin(theta).*pcp; 
upn = sqrt((upcl+upcpr).^2+upcpt.^2); mp = max(max(upn));
% déplacements S
uscl = -US*sin(theta).*eiksr; uscpr = -US*2*cos(theta).*scp;uscpt =-US*sin(theta).*scp;
usn = sqrt(uscpr.^2+(uscl+uscpt).^2); ms = max(max(usn));
% polarisation P
rr = lp*[.4 .95 1.45 1.97 2.5]; 
kpr = kp*rr; eikpr = cos(kpr-phi)./rr; ieikpr = -sin(kpr-phi)./rr; pcp = ((ieikpr-eikpr./kpr)./kpr+1./rr.^3);

ith = 0;
for theta = [-pi/2+pi/8:pi/8:pi/2-pi/8 pi/2+pi/8:pi/8:3*pi/2-pi/8]
ith = ith+1;
xqp(ith,:) = rr*sin(theta); yqp(ith,:) = rr*cos(theta);

upcl = UP*cos(theta).*eikpr; upcpr = UP*2*cos(theta).*pcp; upcpt = UP*sin(theta).*pcp; upq = sqrt((upcl+upcpr).^2+upcpt.^2);
px(ith,:) = ((upcpr+upcl).*sin(theta)+upcpt.*cos(theta))./upq;
pz(ith,:) = ((upcpr+upcl).*cos(theta)-upcpt.*sin(theta))./upq;
end
% polarisation S
rr = ls*[.5 1 1.5 2 2.5 3 3.5 4]; 
ksr = ks*rr; eiksr = cos(ksr-phi)./rr; ieiksr = -sin(ksr-phi)./rr; scp = ((ieiksr-eiksr./ksr)./ksr+1./rr.^3);

ith = 0;
for theta = [pi/8:pi/8:pi-pi/8  pi+pi/8:pi/8:2*pi-pi/8]
ith = ith+1;
xqs(ith,:) = rr*sin(theta); yqs(ith,:) = rr*cos(theta);

uscl = -US*sin(theta).*eiksr; uscpr = -US*2*cos(theta).*scp; uscpt =-US*sin(theta).*scp;usq = sqrt(uscpr.^2+(uscl+uscpt).^2);
sx(ith,:) = (uscpr.*sin(theta)+(uscl+uscpt).*cos(theta))./usq;
sz(ith,:) = (uscpr.*cos(theta)-(uscl+uscpt).*sin(theta))./usq;
end

figure(1)
subplot(2,2,3)
pcolor(xx/lp,yy/lp,upn), caxis([0,mp/700]), shading('interp'), colorbar, hold on
quiver(xqp/lp,yqp/lp,px,pz,.25,'w','Linewidth',1.2)
quiver(0,0,0,cos(phi),.7,'m','Linewidth',2)
set(gca,'XLim',xmax*[-1 1]/lp,'YLim',ymax*[-1 1]/lp,'DataAspectRatio',[1,1,1])
xlabel('X / {\Lambda}'), ylabel('Z / {\Lambda}'), title(['ONDE P : VP = ' num2str(VP) ' M/S '])

subplot(2,2,4)
pcolor(xx/ls,yy/ls,usn), caxis([0,ms/700]), shading('interp'), colorbar, hold on
quiver(xqs/ls,yqs/ls,sx,sz,.5,'w','Linewidth',1.2)
quiver(0,0,0,cos(phi),1.1,'m','Linewidth',2)
set(gca,'XLim',xmax*[-1 1]/ls,'YLim',ymax*[-1 1]/ls,'DataAspectRatio',[1,1,1])
xlabel('X / {\Lambda}'), ylabel('Z / {\Lambda}'), title(['ONDE S : VS = ' num2str(VS,4) ' M/S'])

%gtext(['ONDES P ET S : DEPLACEMENT (M/N) ET POLARISATION POUR FORCE PONCTUELLE VERTICALE HARMONIQUE (F = 100 HZ)'])


% onde P pour r < lambdaP/2
rmin = lp/20; rmax = lp/2; rr = rmin:.1:rmax; 
kpr = kp*rr; eikpr = cos(kpr-phi)./rr; ieikpr = -sin(kpr-phi)./rr; pcp = ((ieikpr-eikpr./kpr)./kpr+1./rr.^3); kprr = kpr.*rr;

figure(2), clf
ith = -1;
for theta = [0 -30 -60]*pi/180;
ith = ith+2;
upcl = cos(theta).*eikpr; upcpr = 2*cos(theta).*pcp; upcpt = sin(theta).*pcp; 

subplot(3,2,ith), hold on,  box on, grid on
plot(rr/lp,kprr.*upcpr,'r',rr/lp,kprr.*upcpt,'b',rr/lp,kprr.*upcl,'r--','Linewidth',2)
plot(rr/lp,kprr.*(upcl+upcpr),'r','Linewidth',3)
plot(rr/lp,1./kpr,'k',rr/lp,-1./kpr,'k')
xlim([rmin rmax]/lp), ylim([-4 4])
if ith ==1, legend('UCPRadial','UCPTheta','UCLRadial','UC(P+L)Radial','1 / (KR)'), end
title(['R^2 * DEPLACEMENT(R/{\Lambda}, Theta = ' num2str(theta*180/pi) ' °)']), xlabel('R/{\Lambda}')
end

% onde P pour r < 3*lambdaP
rmin = lp/10; rmax = 3*lp; rr = rmin:rmax; 
kpr = kp*rr; eikpr = cos(kpr-phi)./rr; ieikpr = -sin(kpr-phi)./rr; pcp = ((ieikpr-eikpr./kpr)./kpr+1./rr.^3); 

ith = 0;
for theta = [0 -30 -60]*pi/180;
ith = ith+2;
upcl = cos(theta).*eikpr; upcpr = 2*cos(theta).*pcp; upcpt = sin(theta).*pcp; 

subplot(3,2,ith), hold on, box on, grid on
plot(rr/lp,rr.*upcpr,'r',rr/lp,rr.*upcpt,'b',rr/lp,rr.*upcl,'r--','Linewidth',2)
plot(rr/lp,rr.*(upcl+upcpr),'r','Linewidth',3)
plot(rr/lp,1./kpr,'k',rr/lp,-1./kpr,'k')
xlim([rmin rmax]/lp),ylim([-2 2])
title(['R * DEPLACEMENT(R/{\Lambda}, Theta = ' num2str(theta*180/pi) ' °)']), xlabel('R/{\Lambda}')

end

%gtext('ONDE P  : TERMES DU DEPLACEMENT EN CHAMP PROCHE ET LOINTAIN, RADIAL ET TANGENTIEL POUR UNE FORCE PONCTUELLE VERTICALE')

% onde S pour r < lambdaS/2
rmin = ls/20; rmax = ls/2; rr = rmin:.1:rmax; 
ksr = ks*rr; eiksr = cos(ksr-phi)./rr; ieiksr = -sin(ksr-phi)./rr; scp = ((ieiksr-eiksr./ksr)./ksr+1./rr.^3); ksrr = ksr.*rr;

figure(3), clf
ith = -1;
for theta = -[90 60 30]*pi/180;
ith = ith+2;
uscl = -sin(theta).*eiksr; uscpr = -2*cos(theta).*scp; uscpt =-sin(theta).*scp;

subplot(3,2,ith), hold on,  box on, grid on
plot(rr/ls,ksrr.*uscpr,'r',rr/ls,ksrr.*uscpt,'b',rr/ls,ksrr.*uscl,'b--','Linewidth',2)
plot(rr/ls,ksrr.*(uscl+uscpt),'b','Linewidth',3)
plot(rr/ls,1./ksr,'k',rr/ls,-1./ksr,'k')
xlim([rmin rmax]/ls), ylim([-4 4])
if ith ==1, legend('UCPRadial','UCPTheta','UCLTheta','UC(P+L)Theta','1 / (KR)'), end
title(['R^2 * DEPLACEMENT (R/{\Lambda}, Theta = ' num2str(theta*180/pi) ' °)']), xlabel('R/{\Lambda}')

end

% onde S pour r < 3*lambdaS
rmin = ls/10; rmax = 3*ls; rr = rmin:rmax; 
ksr = ks*rr; eiksr = cos(ksr-phi)./rr; ieiksr = -sin(ksr-phi)./rr; scp = ((ieiksr-eiksr./ksr)./ksr+1./rr.^3); 

ith = 0;
for theta = -[90 60 30]*pi/180;
ith = ith+2;

uscl = -sin(theta).*eiksr; uscpr = -2*cos(theta).*scp; uscpt =-sin(theta).*scp;

subplot(3,2,ith), hold on,  box on, grid on
plot(rr/ls,rr.*uscpr,'r',rr/ls,rr.*uscpt,'b',rr/ls,rr.*uscl,'b--','Linewidth',2)
plot(rr/ls,rr.*(uscl+uscpt),'b','Linewidth',3)
plot(rr/ls,1./ksr,'k',rr/ls,-1./ksr,'k')
xlim([rmin rmax]/ls), ylim([-2 2])
title(['R * DEPLACEMENT (R/{\Lambda}, Theta = ' num2str(theta*180/pi) ' °)']), xlabel('R/{\Lambda}')

end

gtext('ONDE S  : TERMES DU DEPLACEMENT EN CHAMP PROCHE ET LOINTAIN, RADIAL ET TANGENTIEL, POUR UNE FORCE PONCTUELLE VERTICALE')