%ANISOMOD calcule les modules elastiques anisotropes equivalents pour stratification
% voir : http://jmmeost.free.fr/Sismique/osis.html#IT
% fig1  = http://jmmeost.free.fr/Sismique/anisomod.gif

clear all, clf

VP1 = 7000; VP2 = 3000; VPVS1 = 1.7; VPVS2 = 5; rho1 = 2700; rho2 = 1800; 
p12 = 1/2; i12 = 1/p12; iz1 = 1:i12:10;
iz2 = iz1(1)+1:iz1(1)+i12-1; for ii = 2:length(iz1), iz2 = [iz2 , iz1(ii)+1:iz1(ii)+i12-1]; end, iz2 =iz2(find(iz2<=10));
SAP = 100; % contrainte appliquee

VP = [VP1,VP2]; VS = [VP1/VPVS1,VP2/VPVS2]; RHO = [rho1,rho2];
MU = RHO.*VS.^2; L2M = RHO.*VP.^2; LA = L2M-2*MU; 
PI = [p12,1-p12];

z = 0:100:1000;
x1 = 0:1000:1000; [X1,Z1] = meshgrid(x1,z); UN = ones(size(X1));  
x2 = 1000:1000:2000; [X2,Z2] = meshgrid(x2,z);

figure(1), clf
emax = SAP/min([LA,MU])/2;

% C33 : szz = SAP , exx = eyy = 0

szz = SAP;  C33 = 1/sum(PI./L2M);
EZZ = szz./L2M; SXX = LA.*EZZ;

val = EZZ(1)*UN; val(iz2) = EZZ(2)*UN(iz2); 
subplot(5,3,1), pcolor([X1,X2],[Z1,Z2],[0*UN,val]), axis off, set(gca,'CLim',[-emax,emax]), colorbar
title('EXX = EYY = 0 | EZZ')

val = SXX(1)*UN; val(iz2) = SXX(2)*UN(iz2);
subplot(5,3,2), pcolor([X1,X2],[Z1,Z2],[val,szz*UN]), axis off, set(gca,'CLim',[-SAP,SAP]), colorbar
title(['SXX , SYY | SZZ = ' int2str(SAP)])

val = L2M(1)*UN; val(iz2) = L2M(2)*UN(iz2);
subplot(5,3,3), pcolor([X1,X2],[Z1,Z2],[val,C33*UN]), axis off, set(gca,'CLim',[min([LA,MU]),max(L2M)]),colorbar
title(['LA+2MU  = ' num2str(L2M(1)*1e-9,'%3.3g') ' , ' num2str(L2M(2)*1e-9,'%3.3g') '  |  C33 = ' num2str(C33*1e-9,'%3.3g')])

% C11 : sxx = SAP , eyy = ezz = 0

sxx = SAP;  C11 = sum(PI.*LA./L2M)^2*C33+4*sum(PI.*MU.*(LA+MU)./L2M);
exx = sxx/C11; szz = exx*sum(PI.*LA./L2M)*C33;
EZZ = (szz-LA*exx)./L2M; SXX = L2M*exx+LA.*EZZ;

val = EZZ(1)*UN; val(iz2) = EZZ(2)*UN(iz2); 
subplot(5,3,4), pcolor([X1,X2],[Z1,Z2],[exx*UN,val]), axis off, set(gca,'CLim',[-emax,emax]), colorbar
title('EXX (EYY = 0) | EZZ = 0')

val = SXX(1)*UN; val(iz2) = SXX(2)*UN(iz2);
subplot(5,3,5), pcolor([X1,X2],[Z1,Z2],[val,szz*UN]), axis off, set(gca,'CLim',[-SAP,SAP]), colorbar
title(['SXX = ' int2str(SAP) '  |  SZZ'])

val = L2M(1)*UN; val(iz2) = L2M(2)*UN(iz2);
subplot(5,3,6), pcolor([X1,X2],[Z1,Z2],[val,C11*UN]), axis off, set(gca,'CLim',[min([LA,MU]),max(L2M)]),colorbar
title(['LA+2MU  = ' num2str(L2M(1)*1e-9,'%3.3g') ' , ' num2str(L2M(2)*1e-9,'%3.3g') '   |  C11 = ' num2str(C11*1e-9,'%3.3g')])


% C44 : syz = SAP , eyy = ezz = 0

syz = SAP; C44 = 1/sum(PI./MU); EYZ = syz./MU/2;

val = EYZ(1)*UN; val(iz2) = EYZ(2)*UN(iz2); 
subplot(5,3,7), pcolor([X1,X2],[Z1,Z2],[NaN*UN,val]), axis off, set(gca,'CLim',[-emax,emax]), colorbar
title('EXY  |  EYZ')

subplot(5,3,8), pcolor([X1,X2],[Z1,Z2],[NaN*UN,syz*UN,]), axis off, set(gca,'CLim',[-SAP,SAP]), colorbar
title(['SXY  |  SYZ = ' int2str(SAP)])

val = MU(1)*UN; val(iz2) = MU(2)*UN(iz2);
subplot(5,3,9), pcolor([X1,X2],[Z1,Z2],[val,C44*UN]), axis off, set(gca,'CLim',[min([LA,MU]),max(L2M)]),colorbar
title(['MU  = ' num2str(MU(1)*1e-9,'%3.3g') ' , ' num2str(MU(2)*1e-9,'%3.3g') '  |  C44 = ' num2str(C44*1e-9,'%3.3g')])

% C66 , C12 = C11-2C66 : sxx = SAP , eyy = ezz = 0

sxy = SAP; C66 = sum(PI.*MU); C12 = C11-2*C66; 
exy = sxy/C66/2; SXY = 2*MU*exy;

subplot(5,3,10), pcolor([X1,X2],[Z1,Z2],[exy*UN,NaN*UN]), axis off, set(gca,'CLim',[-emax,emax]), colorbar
title('EXY  |  EYZ')

val = SXY(1)*UN; val(iz2) = SXY(2)*UN(iz2);
subplot(5,3,11), pcolor([X1,X2],[Z1,Z2],[val,NaN*UN]), axis off, set(gca,'CLim',[-SAP,SAP]), colorbar
title(['SXY = ' int2str(SAP) ' |  SYZ'])

subplot(5,3,12), pcolor([X1,X2],[Z1,Z2],[C66*UN,C12*UN]), axis off, set(gca,'CLim',[min([LA,MU]),max(L2M)]),colorbar
title(['C66 = ' num2str(C66*1e-9,'%3.3g') '  |  C12 = ' num2str(C12*1e-9,'%3.3g')])

% C13 : sxx = SAP , szz = 0, eyy = 0

sxx = SAP; szz = 0; C13 = sum(PI.*LA./L2M)*C33;
exx = sxx/(C11-C13^2/C33); EZZ = -exx*LA./L2M; SXX = L2M*exx+LA.*EZZ;

val = EZZ(1)*UN; val(iz2) = EZZ(2)*UN(iz2); 
subplot(5,3,13), pcolor([X1,X2],[Z1,Z2],[exx*UN,val]), set(gca,'CLim',[-emax,emax]), colorbar
title('EXX (EYY = 0) |  EZZ') , xlabel('DEFORMATIONS');

val = SXX(1)*UN; val(iz2) = SXX(2)*UN(iz2);
subplot(5,3,14), pcolor([X1,X2],[Z1,Z2],[val,0*UN]), set(gca,'CLim',[-SAP,SAP]), colorbar
title(['SXX = ' int2str(SAP) ' |  SZZ = 0']), xlabel('CONTRAINTES')

val = LA(1)*UN; val(iz2) = LA(2)*UN(iz2);
subplot(5,3,15), pcolor([X1,X2],[Z1,Z2],[val,C13*UN]),  set(gca,'CLim',[min([LA,MU]),max(L2M)]),colorbar
title(['LA = ' num2str(LA(1)*1e-9,'%3.3g') ' , ' num2str(LA(2)*1e-9,'%3.3g') '  |  C13 = ' num2str(C13*1e-9,'%3.3g')]), xlabel('MODULES ELASTIQUES')

gtext(['VP1 = ' int2str(VP1) ' M/S , VP2 = ' int2str(VP2) ' M/S , VS1 = ' int2str(VS(1)) ' ,  VS2 = ' int2str(VS(2)) '  , rho1 = ' int2str(rho1) ' KG/M3 , rho2 = ' int2str(rho2) ' KG/M3'])
colormap('jet'),