% This function is used to generate the input elastic stiffness matrix for the main function ASCI3D.m.
% The inputs includes:
%     i) The type of anisotropic media: MediaType =1 for isotropic cases, =2 for Transversely isotropic cases, or =3 for Orthorhombic cases;
%     ii) The type of inputs: InputType =1 for inputing Moduli, =2 for inputing velocities, or =3 for inputing elastic stiffnesses;
%     iii) the input parameters;
% The outputs are the stress in cylindrical co-ordinates including:
%     elastic stiffness matrix- Cm.
clear all;
%% Input the type of the media and inputs
MediaType=1;%Media are =1 isotropic; =2 Transversely isotropic; =3 Orthorhombic
InputType=1;%Inputs are =1 Modulus; =2 velocities; =3 elastic stiffnesses
%% Inputs
% We offer examples for each possible inputs here.
% 1) Input the moduli of isotropic cases.
Y=10;% Young's modulus
P=0.2;% Poission's ratio

%2) Input the velocities of isotropic cases.
% density=2.65;%density of the rock.
% vpp=4.5;%Phase velocity of P-wave
% vps=2.3;%Phase velocity of S-wave

%3) Input the elastic stiffnesses of isotropic cases.
% c11=5;%The elastic stiffness C(1,1), GPa
% c12=3;%The elastic stiffness C(1,2), GPa

%4) Input the moduli of transversely isotropic cases.
% Yx=15;% Young's modulus in the direction parallel to the foliation
% Yz=8;% Young's modulus in the direction perpendicular to the foliation
% Pxy=0.2;% Poisson's ratio of strain in y direction over that in x direction.
% Pzx=0.1;% Poisson's ratio of strain in x direction over that in z direction.
% G=5;% Shear modulus

%5) Input the velocities of transversely isotropic cases.
% density=2.48;%density of the media, g/cm3
% vp0=3.7513;%0 degree P-wave phase velocty(vertical to the foliation), km/s
% vp90=5.2521;%90 degree P-wave phase velocty(parallel to the foliation), km/s
% vsh0=2.4013;%0 degree SH-wave phase velocty(vertical to the foliation), km/s
% vsh90=3.3014;%90 degree SH-wave phase velocty(parallel to the foliation), km/s
% vp45=4.4252;%45 degree P-wave phase velocty (45 degree oblique to the foliation), km/s

%6) Input the elastic stiffnesses of transversely isotropic cases.
% c11=68.41;%The elastic stiffness C(1,1), GPa
% c33=34.9;%The elastic stiffness C(3,3), GPa
% c44=14.3;%The elastic stiffness C(4,4), GPa
% c66=27.03;%The elastic stiffness C(6,6), GPa
% c13=12.0;%The elastic stiffness C(1,3), GPa
%% Convert the inputs into elastic stiffness matrix
% The equations come from Mah and Schmitt(2003) and Chan and Schmitt(2014).
if (abs(MediaType-1)<10^(-10) && abs(InputType-1)<10^(-10))
c11=Y*(1-P)/(1+P)/(1-2*P);
c12=Y*P/(1+P)/(1-2*P);
Cm=[c11 c12 c12 0 0 0;
    c12 c11 c12 0 0 0;
    c12 c12 c11 0 0 0;
    0 0 0 (c11-c12)/2 0 0;
    0 0 0 0 (c11-c12)/2 0;
    0 0 0 0 0 (c11-c12)/2];
end
if (abs(MediaType-1)<10^(-10) && abs(InputType-2)<10^(-10))
c11=density*vpp*vpp;
c44=density*vps*vps;
Cm=[c11 c11-2*c44 c11-2*c44 0 0 0;
    c11-2*c44 c11 c11-2*c44 0 0 0;
    c11-2*c44 c11-2*c44 c11 0 0 0;
    0 0 0 c44 0 0;
    0 0 0 0 c44 0;
    0 0 0 0 0 c44];
end
if (abs(MediaType-1)<10^(-10) && abs(InputType-3)<10^(-10))
Cm=[c11 c12 c12 0 0 0;
    c12 c11 c12 0 0 0;
    c12 c12 c11 0 0 0;
    0 0 0 (c11-c12)/2 0 0;
    0 0 0 0 (c11-c12)/2 0;
    0 0 0 0 0 (c11-c12)/2];
end
if (abs(MediaType-2)<10^(-10) && abs(InputType-1)<10^(-10))
Sm=[1/Yx -Pxy/Yx -Pzx/Yz 0 0 0;
    -Pxy/Yx 1/Yx -Pzx/Yz 0 0 0;
    -Pzx/Yz -Pzx/Yz 1/Yz 0 0 0;
    0 0 0 1/G 0 0;
    0 0 0 0 1/G 0;
    0 0 0 0 0 2*(1+Pxy)/Yx;];
Cm=Sm^(-1);
end
if (abs(MediaType-2)<10^(-10) && abs(InputType-2)<10^(-10))
c11=density*vp90^2.0;
c33=density*vp0^2.0;
c44=density*vsh0^2.0;
c66=density*vsh90^2.0;
c13= -c44+ sqrt( (c11+c44-2*density*vp45^2)*(c33+c44-2*density*vp45^2) );
c12=c11-2*c66;
Cm=[c11 c12 c13 0 0 0;
    c12 c11 c13 0 0 0;
    c13 c13 c33 0 0 0;
    0 0 0 c44 0 0;
    0 0 0 0 c44 0;
    0 0 0 0 0 c66];
end
if (abs(MediaType-2)<10^(-10) && abs(InputType-3)<10^(-10))
c12=c11-2*c66;
Cm=[c11 c12 c13 0 0 0;
    c12 c11 c13 0 0 0;
    c13 c13 c33 0 0 0;
    0 0 0 c44 0 0;
    0 0 0 0 c44 0;
    0 0 0 0 0 c66];
end