%ANISOPRO calcule les vitesses de propagation dans un milieu isotrope transverse
% et les coefficients de reflexion-transmission du deplacement des particules.
% voir : http://jmmeost.free.fr/Sismique/osis.html#IT
% fig1  = http://jmmeost.free.fr/Sismique/anisopro.gif
% fig2  = http://jmmeost.free.fr/Sismique/anisorefl.gif

clear all

c11 = 45e9; c12 = 10e9; c13 = 8e9; c33 = 28e9; c44 = 11e9; c66 = (c11-c12)/2; rho = 2380;
aa = c11-c44; bb = c33-c44; cc = c13+c44;

thd = 1:1:89; th = thd*pi/180; UN = ones(1,length(th));

sth = sin(th); cth = cos(th); s2th = sth.^2; c2th = cth.^2; s22th = (sin(2*th)).^2;

% vitesses de phase
D = sqrt((aa*s2th-bb*c2th).^2+cc^2*s22th);
VSH = sqrt((c66*s2th+c44*c2th)/rho);
VQP = sqrt((c44+c11*s2th+c33*c2th+D)/2/rho);
VQS = sqrt((c44+c11*s2th+c33*c2th-D)/2/rho);

%vitesses de groupe
VGH = sqrt(c66^2*s2th+c44^2*c2th)/rho./VSH;
dDdth = ((aa*s2th-bb*c2th).*((c11+c33-2*c44)*sin(2*th))+cc^2*sin(4*th))./D;
dVPdth = ((c11-c33)*sin(2*th)+dDdth)./(4*rho*VQP); VGP = sqrt(VQP.^2+dVPdth.^2);
dVSdth = ((c11-c33)*sin(2*th)-dDdth)./(4*rho*VQS); VGS = sqrt(VQS.^2+dVSdth.^2);

% angle polarisation
al = atan(cc*sth.*cth./(rho*VQP.^2-c11*s2th-c44*c2th)); ald = 180*al/pi;
be = atan(cc*sth.*cth./(rho*VQS.^2-c11*s2th-c44*c2th)); bed = 180*be/pi;
alm = pi+atan(-cc*sth.*cth./(rho*VQP.^2-c11*s2th-c44*c2th)); almd = 180*alm/pi;
bem = atan(-cc*sth.*cth./(rho*VQS.^2-c11*s2th-c44*c2th)); bemd = 180*bem/pi;

% angle groupe
eg = atan((c66/c44)*tan(th)); egd = 180*eg/pi;
thPg = atan((VQP.*sth+cth.*dVPdth)./(VQP.*cth-sth.*dVPdth)); thPgd = 180*thPg/pi;
thSg = atan((VQS.*sth+cth.*dVSdth)./(VQS.*cth-sth.*dVSdth)); thSgd = 180*thSg/pi;

% parametre Thomsen et vitesses phases pour faible anisotropie
alt = sqrt(c33/rho); bet = sqrt(c44/rho); gat = (c66/c44-1)/2; ept = (c11/c33-1)/2; 
det = (cc^2-bb^2)/(2*c33*bb);
VSHT = bet*(1+gat*s2th);
VQPT = alt*(1+det*s2th+(ept-det)*s2th.^2);
VQST = bet*(1+(alt^2/bet^2)*(ept-det)*s2th.*c2th);

figure(1), clf

subplot(2,2,1)
plot(thd,VQP,'r',thd,VQS,'b',thd,VSH,'m',thd,VQPT,'-.r',thd,VQST,'-.b',thd,VSHT,'-.m'), grid on, axis([0 90 2000 5000])
legend('VQP','VQS','VSH','VQPTH','VQSTH','VSHTH'), xlabel('ANGLE INCIDENCE'), ylabel('VITESSE PHASE ET THOMSEN (M/S)')
title('ISOTROPIE TRANSVERSE HORIZONTALE')

subplot(2,2,2)
plot(thd,ald-thd,'r',thd,bed-(thd-90),'b-.'), grid on 
legend('ALPHAD - THETA','BETAD - (THETA-90)') 
xlabel('ANGLE INCIDENCE'), ylabel('ANGLE POLARISATION ANISO - ISO')
title(['C11 = ' int2str(c11*1e-9) , ' , C33 = ' int2str(c33*1e-9) ' , C12 = ' int2str(c12*1e-9) ' , C13 = ' int2str(c13*1e-9) ...
    ' , C44 = ' int2str(c44*1e-9) ' , C66 = ' int2str(c66*1e-9)])

subplot(2,2,3)
plot(thd,100*(VGP-VQP)./VQP,'r',thd,100*(VGS-VQS)./VQS,'b',thd,100*(VGH-VSH)./VSH,'m'), grid on, axis([0 90 0 5])
legend('(VGP/VQP-1)%','(VGS/VQS-1)%','(VGH/VSH-1)%'), xlabel('ANGLE INCIDENCE'), ylabel('ECART VITESSES GROUPE ET PHASE (%)')

subplot(2,2,4)
plot(thd,thPgd-thd,'r',thd,thSgd-thd,'b',thd,egd-thd,'m'), grid on, axis([0 90 -20 20])
legend('QP','QS','SH'), xlabel('ANGLE INCIDENCE'), ylabel('ANGLE GROUPE - ANGLE PHASE')

% coef. reflexion transmission
d11 = c11; d33 = c33; d12 = c12; d13 = c13; d44 = c44; d66 = c66; % milieu 2 anisotrope

VP1 = 3000; VS1 = 2000; rho1 = 2500; % milieu 1 isotrope
c11 = rho1*VP1^2; c33 = c11; c12 = c11-2*rho1*VS1^2; c13 = c12; c44 = (c11-c12)/2; c66 = c44; VQP1 = VP1; VQS1 = VS1;
VQP2H = sqrt(d11/rho); th1c = asin(VP1/VQP2H); th1cd = 180*th1c/pi;

th1d = 2:2:90; th1 = th1d*pi/180; p = sin(th1)/VP1; p2 = p.^2; np = length(p); % p pour onde P incidence theta1
reth2 = find(th1 < th1c); imth2 = find(th1 >= th1c); %indices correspondant a QP2 transmise ou evanescente
qP1 = sqrt(1/VP1^2-p2); qS1 = sqrt(1/VS1^2-p2);
et1 = asin(p*VS1); % eta1 pour S reflechie
al1 = th1; be1 = et1; % polarisation P et S isotrope calculee en fonction angle incidence

AA = d33*d44 ; BB = -d33*(rho-d11*p2)-d44*(rho-d44*p2)-(d13+d44)^2*p2 ; CC = (rho-d11*p2).*(rho-d44*p2); DD = sqrt(BB.^2-4*AA*CC);
qP22 = (-BB-DD)/2/AA; qP2 = sqrt(qP22); %resolution equation dispersion pour qP(p)
qS22 = (-BB+DD)/2/AA; qS2 = sqrt(qS22); %resolution equation dispersion pour qS(p)
th2 = atan(p./qP2); th2d = 180*th2/pi; %incidence QP2
et2 = atan(p./qS2); et2d = 180*et2/pi; %incidence QS2
al2 = atan(cc*p.*qP2./(rho-d11*p2-d44*qP22)); al2d = 180*al2/pi; %polarisation QP
be2 = atan(cc*p.*qS2./(rho-d11*p2-d44*qS22))+pi/2; be2d = 180*be2/pi; %polarisation QS

sal1 = sin(al1); cal1 = cos(al1); 
sbe1 = sin(be1); cbe1 = cos(be1); 
sal2 = sin(al2); cal2 = cos(al2); 
sbe2 = sin(be2); cbe2 = cos(be2); 

A1 = [ sal1 ; -cal1 ;  (c13*sal1.*p+c33*cal1.*qP1) ; c44*(sal1.*qP1+cal1.*p)];
A2 = [ cbe1 ;  sbe1 ;  (c13*cbe1.*p-c33*sbe1.*qS1) ; c44*(cbe1.*qS1-sbe1.*p)];
A3 = [-sal2 ; -cal2 ; -(d13*sal2.*p+d33*cal2.*qP2) ; d44*(sal2.*qP2+cal2.*p)];
A4 = [ cbe2 ; -sbe2 ;  (d13*cbe2.*p-d33*sbe2.*qS2) ; d44*(-cbe2.*qS2+sbe2.*p)];

A = [reshape(A1,4,1,np) , reshape(A2,4,1,np) , reshape(A3,4,1,np) , reshape(A4,4,1,np)];

B = [-sal1 ; -cal1 ; -(c13*sal1.*p+c33*cal1.*qP1) ; c44*(sal1.*qP1+cal1.*p)]; 

for ip = 1:np
    R(:,ip) = A(:,:,ip)\B(:,ip);
end

figure(2), clf

subplot(2,2,1), hold on, grid on, box on
plot(-p,qS1,'bo',-p,qS2,'go',-p,qP1,'ro',-p(reth2),qP2(reth2),'mo',-p(imth2),imag(qP2(imth2)),'m+')
polar(th+pi/2,UN/VS1,'b'), polar(th+pi/2,1./VQS,'g'), polar(th+pi/2,UN/VP1,'r'), polar(th+pi/2,1./VQP,'m')
axis equal, axis([-5e-4 0 0 5e-4]), xlabel('P = SIN (THETA) / V (S/M)') , ylabel('Q = COS (THETA) / V (S/M)')
legend('QS1','QS2','QP1','QP2','IM(QP2)'), title('VITESSES ANISOTROPES')

subplot(2,2,3), hold on, grid on, box on
plot(th1d(reth2),al2d(reth2),'m.',th1d(reth2),th2d(reth2),'m',th1d,be2d,'g.',th1d,et2d,'g')
axis([0 90 0 90]),xlabel('THETA1'), ylabel('ANGLE (DEG)')
legend('ALPHA2','THETA2','BETA2','ETA2'), title('ANGLES INCIDENCE ET POLARISATION QP2 et QS2')

subplot(2,2,2), hold on, grid on, box on
plot(th1d,abs(R(1,:)),'r',th1d,abs(R(2,:)),'b',th1d,abs(R(3,:)),'m',th1d,abs(R(4,:)),'g')
plot(th1cd,0,'mo')
axis([0 90 0 2]), xlabel('THETA1'), ylabel('MODULE')
legend('RPP','RPS','TPQP','TPQS'), title('COEFFICIENTS DE REFLEXION-TRANSMISSION')

subplot(2,2,4), hold on, grid on, box on
plot(th1d,angle(R(1,:))*180/pi,'r',th1d,angle(R(2,:))*180/pi,'b',th1d,angle(R(3,:))*180/pi,'m',th1d,angle(R(4,:))*180/pi,'g')
axis([0 90 -180 180]),xlabel('THETA1'), ylabel('PHASE (DEG)')
legend('RPP','RPS','TPQP','TPQS');