In [11]:
import numpy as np
np.set_printoptions(precision=2)
def Cij(e, d, y, p, Vp, Vs):
"""
Returns the elastic stiffness elements C11, C33, C13, C55, and C66 that
characterize transversely isotropic materials, using the Thomsen parameters
and elastic parameters.
"""
C = np.zeros(shape=(6, 6))
f = 1 - (Vs/Vp)**2
dtil = f*(np.sqrt(1 + 2*d/f) - 1)
C[0][0] = p*Vp**2*(1 + 2*e)
C[1][1] = C[0][0]
C[2][2] = p*Vp**2
C[3][3] = p*Vp**2*(1 - f)
C[4][4] = C[3][3]
C[5][5] = p*Vp**2*(1 - f)*(1 + 2*y)
C[0][2] = p*Vp**2*(2*f + dtil - 1)
C[2][0] = C[0][2]
C[0][1] = C[0][0] - 2*C[5][5]
C[1][0] = C[0][1]
return(C)
def inverse3(M):
M11=M[0][0]
M22=M[1][1]
M33=M[2][2]
M12=M[0][1]
M13=M[0][2]
M23=M[1][2]
M21=M[1][0]
M31=M[2][0]
M32=M[2][1]
temp = M11*M22*M33-M11*M23*M32-M21*M12*M33+M21*M13*M32+M31*M12*M23-M31*M13*M22
temp = 1.0 / temp
iM = np.zeros(shape=(3, 3))
iM[0][0] = (M22*M33-M23*M32) * temp
iM[0][1] = (M13*M32-M12*M33) * temp
iM[0][2] = (M12*M23-M13*M22) * temp
iM[1][0] = (M23*M31-M21*M33) * temp
iM[1][1] = (M11*M33-M13*M31) * temp
iM[1][2] = (M21*M13-M23*M11) * temp
iM[2][0] = (M32*M21-M31*M22) * temp
iM[2][1] = (M31*M12-M32*M11) * temp
iM[2][2] = (M11*M22-M12*M21) * temp
return iM
def order3incr(temp3):
A = temp3[0]
B = temp3[1]
C = temp3[2]
P = np.min(A,B)
Q = np.max(A,B)
A = P
B = Q
P = np.min(B,C)
Q = np.max(B,C)
B = P
temp3[2] = Q
temp3[0] = np.min(A,B)
temp3[1] = np.max(A,B)
return temp3;
def eigvec3symm(M, ev):
M11=M[0][0]
M22=M[1][1]
M33=M[2][2]
M12=M[0][1]
M13=M[0][2]
M23=M[1][2]
M12sq = M12*M12
M13sq = M13*M13
M23sq = M23*M23
x1 = (M22-ev)*M13 - M12*M23
x2 = (M11-ev)*M23 - M12*M13
x3 = - (ev*ev -(M11+M22)*ev + M11*M22-M12sq)
temp = 1 / np.sqrt( x1*x1 + x2*x2 + x3*x3 )
eigvec = np.zeros(shape=(3, 1))
eigvec[0] = x1*temp
eigvec[1] = x2*temp
eigvec[2] = x3*temp
return(eigvec)
def cubic_num_real_roots(A, B, C, D):
a1 = B/A
a2 = C/A
a3 = D/A
Q = (a1*a1-3*a2)/9.0
R = (2*a1*a1*a1-9*a1*a2+27*a3)/54.0
a1 = Q*Q*Q - R*R
num = 0
if (a1>=0):
num = 3
if (a1<0):
num = 1
return(num)
def cubic_three_real_roots(A, B, C, D):
a1 = B/A
a2 = C/A
a3 = D/A
Q = (a1*a1-3*a2)/9.0
R = (2*a1*a1*a1-9*a1*a2+27*a3)/54.0
theta = np.arccos( R / np.sqrt(Q*Q*Q) )
roots = np.zeros(shape=(3, 1))
roots[0] = -2.0 * np.sqrt(Q) * np.cos( theta/3.0 ) - a1/3.0
roots[1] = -2.0 * np.sqrt(Q) * np.cos( (theta+2.0*3.14159)/3.0 ) - a1/3.0
roots[2] = -2.0 * np.sqrt(Q) * np.cos( (theta+4.0*3.14159)/3.0 ) - a1/3.0
roots = np.roots(np.array([A, B, C, D]))
return roots
def eigval3symm(M):
M11=M[0][0]
M22=M[1][1]
M33=M[2][2]
M12=M[0][1]
M13=M[0][2]
M23=M[1][2]
M12sq = M12*M12
M13sq = M13*M13
M23sq = M23*M23
A = -1.
B = M11 + M22 + M33
C = M12sq + M13sq + M23sq - M11*M22 - M11*M33 - M22*M33
D = M11*M22*M33 - M11*M23sq - M22*M13sq - M33*M12sq + 2.0*M12*M13*M23
num_roots = cubic_num_real_roots(A,B,C,D)
mv=0
if (num_roots == 1):
eigvals[0] = cubic_one_real_root(A,B,C,D)
temp22 = cubic_two_imag_roots(A,B,C,D,eigvals[0]);
eigvals[1] = temp22[0][0]
eigvals[2] = temp22[1][0]
if (num_roots == 3):
eigvals = cubic_three_real_roots(A,B,C,D)
return(eigvals)
def bicubic_coeffs_monoclinic(s1, s2, rho,c11, c12, c13, c16, c22, c23, c26, c33, c36, c44, c45, c55, c66):
t1 = c45*c45
t2 = c55 * c44
a3 = c33 * (-t1 + t2)
t4 = s1*s1
t5 = c13 * c13
t6 = t4 * t5
t8 = c45 * s1
t9 = s2 * c12
t10 = t9 * c33
t13 = c55 * t4
t14 = c36 * c36
t16 = rho * c44
t18 = s2 * s2
t19 = t18 * t14
t21 = t1 * t4
t24 = c55 * t18
t25 = c23 * c23
t29 = c55 * rho
t31 = t1 * t18
t34 = c11 * t4
t35 = c44 * c33
t37 = c66 * t18
t39 = c23 * c44
t42 = c16 * t4
t43 = c45 * c33
t46 = -t6 * c44 - 2 * t8 * t10 - t13 * t14 - t16 * c33 - t19 * c44 + 2 * t21 * c13 - t24 * t25 - t2 * rho + t1 * rho - t29 * c33 + 2 * t31 * c23 + t34 * t35 + t37 * t35 - 2 * t24 * t39 - 2 * t42 * t43
t47 = c26 * t18
t51 = c13 * c36
t54 = c55 * c66
t55 = t4 * c33
t57 = c55 * c22
t58 = t18 * c33
t61 = c36 * c23
t64 = c45 * s2
t68 = t1 * s2
t72 = t4 * c13
t75 = c55 * c26
t76 = s1 * s2
t77 = t76 * c33
t80 = c55 * s2
t81 = c23 * s1
t82 = t81 * c36
t86 = c44 * s1 * c36
t89 = c16 * s1
t90 = s2 * c44
t94 = s1 * c13
t95 = s2 * c36
t96 = t95 * c44
t99 = s2 * c66
t103 = c13 * s2
t104 = t103 * c23
t107 = -2 * t47 * t43 + 2 * c45 * t4 * t51 + t54 * t55 + t57 * t58 + 2 * c45 * t18 * t61 + 2 * t64 * t14 * s1 + 4 * t68 * c36 * s1 - 2 * t72 * t2 + 2 * t75 * t77 - 2 * t80 * t82 - 4 * t80 * t86 + 2 * t89 * t90 * c33 - 2 * t94 * t96 - 2 * t8 * t99 * c33 + 2 * t8 * t104
a2 = t46 + t107
t108 = s1 * c55
t109 = t18 * s2
t111 = t109 * c36 * c22
t114 = t4 * t4
t117 = rho * t18
t119 = t18 * t18
t122 = c26 * c26
t123 = t122 * t119
t125 = rho * t4
t131 = c16 * c16
t132 = t131 * t114
t134 = c11 * t114
t137 = rho * rho
t140 = -2 * t108 * t111 - t114 * t5 * c66 + t117 * t25 - t119 * t14 * c22 - t123 * c33 + t125 * t1 - t119 * t1 * c22 + t125 * t14 + t6 * rho - t132 * c33 - t134 * t1 + t31 * rho + c55 * t137 + t19 * rho
t142 = c26 * t109
t143 = s1 * c12
t148 = c66 * t119
t150 = t4 * s1
t151 = c11 * t150
t152 = c26 * s2
t156 = s2 * rho
t167 = c22 * t18
t170 = t4 * t18
t171 = c12 * c13
t177 = c36 * c45
t180 = rho * c33
t182 = c44 * t137 - 2 * t142 * t143 * c33 - t134 * t14 - t148 * t25 + 2 * t151 * t152 * c33 - 2 * t89 * t156 * c33 + 2 * t8 * t99 * rho + 2 * t142 * t94 * c23 + t137 * c33 + t34 * t167 * c33 + 2 * t170 * t171 * c44 + 2 * t117 * t39 + 2 * t125 * t177 - t34 * t180
t202 = c44 * t18
t204 = c16 * t150
t207 = -2 * t134 * t177 - t34 * t18 * t25 - t34 * t16 - t37 * t180 - 2 * t148 * t39 - 4 * t37 * t21 - t37 * t16 - t54 * t125 + t134 * c66 * c33 + t134 * t2 + t148 * c22 * c33 - t57 * t117 - t29 * t202 - 4 * t204 * t68
t218 = c16 * t114
t227 = c45 * rho
t230 = c26 * t119
t246 = -2 * t89 * t109 * t25 - rho * c66 * t55 - rho * c22 * t58 - t16 * t13 + 2 * t72 * t29 + 2 * t218 * t51 + 2 * t218 * c13 * c45 + 2 * t218 * c55 * c36 + 2 * t42 * t227 + 2 * t230 * t61 + 2 * t230 * c36 * c44 + 2 * t230 * c45 * c23 - 4 * t142 * t1 * s1 + 2 * t47 * t227 + t57 * t119 * c44
t255 = c12 * c12
t272 = c12 * c36
t279 = t18 * c23
t283 = s2 * c23
t292 = -2 * t21 * t18 * c12 + 2 * t170 * c12 * t14 - t170 * t255 * c33 - 2 * t114 * c13 * t54 - t6 * t167 - 2 * t150 * t5 * t152 - 2 * t119 * c36 * c45 * c22 + 2 * t18 * c36 * t227 + 4 * t170 * t272 * c45 + 2 * t170 * t54 * c23 - 2 * t34 * t279 * c44 - 2 * t151 * t283 * c36 - 2 * t151 * t283 * c45 - 2 * t151 * t96
t296 = s1 * t109
t307 = t109 * c22
t339 = 4 * t37 * t2 * t4 + 4 * t75 * t296 * c44 - 2 * t75 * t76 * rho - 4 * t89 * t109 * c23 * c44 + 2 * t89 * t307 * c33 + 4 * t204 * t90 * c55 - 2 * t89 * t90 * rho + 2 * t42 * t47 * c33 - 2 * t42 * t279 * c36 - 2 * t42 * t279 * c45 - 2 * t42 * t202 * c36 - 2 * rho * c26 * t77 + 2 * t156 * t82 + 2 * t156 * t81 * c45 + 2 * t156 * t86
t344 = t95 * rho
t366 = t150 * s2
t369 = -t94 * t109 * c45 * c22 + t94 * t344 + t94 * t64 * rho - t204 * t10 + t204 * t104 + t204 * t103 * c44 + t204 * t80 * c23 - t47 * t72 * c36 - t47 * t72 * c45 - t47 * t13 * c36 + t142 * t94 * c44 + t142 * t108 * c23 + t8 * t9 * rho + t366 * t171 * c36
t374 = c12 * c55
t388 = c12 * c66
t401 = c66 * c13
t419 = 2 * t366 * t171 * c45 + 2 * t366 * t374 * c36 + 2 * t296 * t272 * c23 + 2 * t296 * t272 * c44 + 2 * t296 * c12 * c45 * c23 - 2 * t170 * t388 * c33 + 2 * t170 * t171 * c23 + 2 * t170 * t374 * c23 + 2 * t170 * t374 * c44 + 2 * t170 * t401 * c23 + 2 * t170 * t401 * c44 - 2 * t94 * t111 - 2 * t72 * t57 * t18 - 4 * t150 * c13 * t75 * s2 + 2 * t108 * t344
a1 = t140 + t182 + t207 + t246 + t292 + t339 + 2 * t369 + t419
t451 = -t132 + t134 * c66 + 2 * t151 * t152 - 2 * t204 * t9 - t34 * rho - t170 * t255 - 2 * t170 * t388 + t34 * t167 - c66 * t4 * rho + 2 * t42 * t47 + 2 * t89 * t307 - 2 * t89 * t156 - 2 * t142 * t143 - 2 * c26 * s1 * t156 - t123 + t137 - t167 * rho + t148 * c22 - t37 * rho
a0 = (t202 + 2 * t8 * s2 + t13 - rho) * t451
temp4 = np.array([a0, a1, a2, a3])
return(temp4)
def axial_rotation(chi, c):
# Calculates a rotated elastic constant matrix for a rotation of angle chi
# about the z-axis
# see www.engin.brown.edu/courses/en175/Notes/Elastic_material/Elastic_material.htm
# see also Winterstein and Hanten, 1990, (Geophysics) p. 1075
cchi = np.cos(chi)
schi = np.sin(chi)
K=np.array([[cchi*cchi, schi*schi, 0, 0, 0, 2.*cchi*schi],
[schi*schi, cchi*cchi, 0, 0, 0, -2.*cchi*schi],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, cchi, -schi, 0],
[0, 0, 0, schi, cchi, 0],
[-cchi*schi, cchi*schi, 0, 0, 0, cchi*cchi-schi*schi]])
Kp = K.T
c_rot = K.dot(c).dot(Kp)
return(c_rot)
def ortho_Rpp(cU, cL,rhoU, rhoL, chiU, chiL, phi, tht):
phi = np.radians(phi)
tht = np.radians(tht)
cphi = np.cos(phi)
sphi = np.sin(phi)
stht = np.sin(tht)
ctht = np.cos(tht)
n1 = cphi*stht
n2 = sphi*stht
n3 = ctht
# Rotate by chi
cU = axial_rotation(chiU,cU)
c11Ur=cU[0][0]
c12Ur=cU[0][1]
c13Ur=cU[0][2]
c16Ur=cU[0][5]
c22Ur=cU[1][1]
c23Ur=cU[1][2]
c26Ur=cU[1][5]
c33Ur=cU[2][2]
c36Ur=cU[2][5]
c44Ur=cU[3][3]
c45Ur=cU[3][4]
c55Ur=cU[4][4]
c66Ur=cU[5][5]
# System.out.println("cU: "+Arrays.deepToString(cU));
cL = axial_rotation(chiL,cL)
c11Lr=cL[0][0]
c12Lr=cL[0][1]
c13Lr=cL[0][2]
c16Lr=cL[0][5]
c22Lr=cL[1][1]
c23Lr=cL[1][2]
c26Lr=cL[1][5]
c33Lr=cL[2][2]
c36Lr=cL[2][5]
c44Lr=cL[3][3]
c45Lr=cL[3][4]
c55Lr=cL[4][4]
c66Lr=cL[5][5]
# Form the velocity form of the Christoffel matrix for the upper layer
CM=np.array([[c11Ur*n1*n1 + c66Ur*n2*n2 + c55Ur*n3*n3 + 2*c16Ur*n1*n2,
c16Ur*n1*n1 + c26Ur*n2*n2 + c45Ur*n3*n3 + (c12Ur+c66Ur)*n1*n2,
(c13Ur+c55Ur)*n1*n3 + (c36Ur+c45Ur)*n2*n3],
[c16Ur*n1*n1 + c26Ur*n2*n2 + c45Ur*n3*n3 + (c12Ur+c66Ur)*n1*n2,
c66Ur*n1*n1 + c22Ur*n2*n2 + c44Ur*n3*n3 + 2*c26Ur*n1*n2,
(c36Ur+c45Ur)*n1*n3 + (c23Ur+c44Ur)*n2*n3],
[(c13Ur+c55Ur)*n1*n3 + (c36Ur+c45Ur)*n2*n3,
(c36Ur+c45Ur)*n1*n3 + (c23Ur+c44Ur)*n2*n3,
c55Ur*n1*n1+c44Ur*n2*n2+c33Ur*n3*n3]])
CM = CM/rhoU
#System.out.println("CM: "+Arrays.deepToString(CM))
# Find eigenvalues (phase velocities) and eigenvectors (phase polarization
# vectors) from the Christoffel matrix just obtained, and use to create the
# slowness vector for the incident wave
temp3 = eigval3symm(CM)
#System.out.println("temp3: "+Arrays.toString(temp3));
# double Vinc = Math.sqrt(temp3[1-1]); // this assumes an incident T-wave
# double Vinc = Math.sqrt(temp3[2-1]); // this assumes an incident S-wave
Vinc = np.sqrt(max(temp3)) # this assumes an incident P-wave
s1 = n1/Vinc
s2 = n2/Vinc
s3inc = n3/Vinc
print('Rotated Stiffness', cU)
print('n', n1, n2, n3)
print('L', CM)
print('vp1', Vinc)
print('s', s1, s2, s3inc)
# Generate coefficients for bicubic for upper layer
# [a3 a2 a1 a0]=bicubic_coeffs_orthotropic(s1,s2,rhoL,c11L,c12L,c13L,c22L,c23L,c33L,c44L,c55L,c66L);
temp4=bicubic_coeffs_monoclinic(s1,s2,rhoU,c11Ur,c12Ur,c13Ur,c16Ur,c22Ur,c23Ur,c26Ur,c33Ur,c36Ur,c44Ur,c45Ur,c55Ur,c66Ur)
# Now solve the bicubic
num_roots = cubic_num_real_roots(temp4[3],temp4[2],temp4[1],temp4[0])
if num_roots == 1:
print("Guess what, you only have one real root in your bicubic (upper)")
temp3 = cubic_three_realified_roots(temp4[3],temp4[2],temp4[1],temp4[0])
if num_roots == 3:
temp3 = cubic_three_real_roots(temp4[3],temp4[2],temp4[1],temp4[0])
print('Coeffs')
print(temp4)
print('Roots')
print(temp3)
# Assign roots
#temp3 = order3incr(temp3)
temp = np.sort(temp3)
s3U2P = temp3[2]
s3U2S = temp3[1]
s3U2T = temp3[0]
# Form vertical slownesses
# HERE - could have imaginary vertical slownesses
s3UP = np.sqrt(s3U2P)
s3US = np.sqrt(s3U2S)
s3UT = np.sqrt(s3U2T)
print('s1P', s1, s2, s3UP)
print('s1S', s1, s2, s3US)
print('s1T', s1, s2, s3UT)
# Generate coefficients for bicubic for lower layer
#[a3 a2 a1 a0]=bicubic_coeffs_orthotropic(s1,s2,rhoL,c11L,c12L,c13L,c22L,c23L,c33L,c44L,c55L,c66L);
temp4=bicubic_coeffs_monoclinic(s1,s2,rhoL,c11Lr,c12Lr,c13Lr,c16Lr,c22Lr,c23Lr,c26Lr,c33Lr,c36Lr,c44Lr,c45Lr,c55Lr,c66Lr)
print('lower layer coef')
print(temp4)
# Now solve the bicubic
#num_roots = cubic_num_real_roots(temp4[0],temp4[1],temp4[2],temp4[3]);
num_roots = cubic_num_real_roots(temp4[3],temp4[2],temp4[1],temp4[0])
if num_roots == 1:
print("Guess what, you only have one real root in your bicubic (upper)")
temp3 = cubic_three_realified_roots(temp4[3],temp4[2],temp4[1],temp4[0])
if (num_roots == 3):
temp3 = cubic_three_real_roots(temp4[3],temp4[2],temp4[1],temp4[0])
# Assign roots
#temp3 = order3incr(temp3);
temp3 = np.sort(temp3)
s3L2P = temp3[0]
s3L2S = temp3[1]
s3L2T = temp3[2]
# Form vertical slownesses
# HERE - could have imaginary vertical slownesses
s3LP = np.sqrt(s3L2P)
s3LS = np.sqrt(s3L2S)
s3LT = np.sqrt(s3L2T)
print('LL roots')
print(temp3)
print('Slowness of Trans.')
print('s3LP', s3LP)
print('s3LS', s3LS)
print('s3LT', s3LT)
# Form the velocity form of the Christoffel matrix for the reflected P-wave
CMUP = np.array([[c11Ur*s1*s1 + c66Ur*s2*s2 + c55Ur*s3U2P + 2*c16Ur*s1*s2,
c16Ur*s1*s1 + c26Ur*s2*s2 + c45Ur*s3U2P + (c12Ur+c66Ur)*s1*s2,
(c13Ur+c55Ur)*s1*s3UP + (c36Ur+c45Ur)*s2*s3UP],
[c16Ur*s1*s1 + c26Ur*s2*s2 + c45Ur*s3U2P + (c12Ur+c66Ur)*s1*s2,
c66Ur*s1*s1 + c22Ur*s2*s2 + c44Ur*s3U2P + 2*c26Ur*s1*s2,
(c36Ur+c45Ur)*s1*s3UP + (c23Ur+c44Ur)*s2*s3UP],
[(c13Ur+c55Ur)*s1*s3UP + (c36Ur+c45Ur)*s2*s3UP,
(c36Ur+c45Ur)*s1*s3UP + (c23Ur+c44Ur)*s2*s3UP,
c55Ur*s1*s1+c44Ur*s2*s2+c33Ur*s3U2P]])
# HERE - could have imaginary vertical slowness above - affects next step?
eigUP = eigvec3symm(CMUP,rhoU)
eigprint = eigUP
eigvp = eigval3symm(CMUP)
temp1 = np.abs(np.sqrt(eigUP[0]*eigUP[0] + eigUP[1]*eigUP[1] + eigUP[2]*eigUP[2]))
temp1 = (eigUP[0]*s1+eigUP[1]*s2)/np.abs(eigUP[0]*s1+eigUP[1]*s2) / temp1
eigUP = eigUP*temp1
print('CMUP', CMUP)
print()
print('eigUP')
print(eigUP)
# Form the velocity form of the Christoffel matrix for the reflected S-wave
CMUS = np.array([[c11Ur*s1*s1 + c66Ur*s2*s2 + c55Ur*s3U2S + 2*c16Ur*s1*s2,
c16Ur*s1*s1 + c26Ur*s2*s2 + c45Ur*s3U2S + (c12Ur+c66Ur)*s1*s2,
(c13Ur+c55Ur)*s1*s3US + (c36Ur+c45Ur)*s2*s3US],
[c16Ur*s1*s1 + c26Ur*s2*s2 + c45Ur*s3U2S + (c12Ur+c66Ur)*s1*s2,
c66Ur*s1*s1 + c22Ur*s2*s2 + c44Ur*s3U2S + 2*c26Ur*s1*s2,
(c36Ur+c45Ur)*s1*s3US + (c23Ur+c44Ur)*s2*s3US],
[(c13Ur+c55Ur)*s1*s3US + (c36Ur+c45Ur)*s2*s3US,
(c36Ur+c45Ur)*s1*s3US + (c23Ur+c44Ur)*s2*s3US,
c55Ur*s1*s1+c44Ur*s2*s2+c33Ur*s3U2S]])
eigUS = eigvec3symm(CMUS,rhoU)
temp1 = np.abs(np.sqrt( eigUS[0]*eigUS[0] + eigUS[1]*eigUS[1] + eigUS[2]*eigUS[2]))
temp1 = (eigUS[0]*s1+eigUS[1]*s2)/np.abs(eigUS[0]*s1+eigUS[1]*s2) / temp1;
eigUS = eigUS*temp1;
# Form the velocity form of the Christoffel matrix for the reflected T-wave
CMUT = np.array([[c11Ur*s1*s1 + c66Ur*s2*s2 + c55Ur*s3U2T + 2*c16Ur*s1*s2,
c16Ur*s1*s1 + c26Ur*s2*s2 + c45Ur*s3U2T + (c12Ur+c66Ur)*s1*s2,
(c13Ur+c55Ur)*s1*s3UT + (c36Ur+c45Ur)*s2*s3UT],
[c16Ur*s1*s1 + c26Ur*s2*s2 + c45Ur*s3U2T + (c12Ur+c66Ur)*s1*s2,
c66Ur*s1*s1 + c22Ur*s2*s2 + c44Ur*s3U2T + 2*c26Ur*s1*s2,
(c36Ur+c45Ur)*s1*s3UT + (c23Ur+c44Ur)*s2*s3UT],
[(c13Ur+c55Ur)*s1*s3UT + (c36Ur+c45Ur)*s2*s3UT,
(c36Ur+c45Ur)*s1*s3UT + (c23Ur+c44Ur)*s2*s3UT,
c55Ur*s1*s1+c44Ur*s2*s2+c33Ur*s3U2T]])
eigUT = eigvec3symm(CMUT,rhoU)
temp1 = np.abs(np.sqrt( eigUT[0]*eigUT[0] + eigUT[1]*eigUT[1] + eigUT[2]*eigUT[2] ))
temp1 = (eigUT[0]*s1+eigUT[1]*s2)/np.abs(eigUT[0]*s1+eigUT[1]*s2) / temp1
eigUT = eigUT*temp1
# try and match up qSV and qSH with the most appropriate eigenvalues/eigenvectors
temp1 = eigUS[0]*n3 + eigUS[1]*n2 - eigUS[2]*n1
temp = eigUT[0]*n3 + eigUT[1]*n2 - eigUT[2]*n1
if (temp > temp1):
temp1 = s3US
s3US = s3UT
s3UT = temp1
temp3 = eigUS
eigUS = eigUT
eigUT = temp3
# Form the velocity form of the Christoffel matrix for the transmitted P-wave
CMLP = np.array([[c11Lr*s1*s1 + c66Lr*s2*s2 + c55Lr*s3L2P + 2*c16Lr*s1*s2,
c16Lr*s1*s1 + c26Lr*s2*s2 + c45Lr*s3L2P + (c12Lr+c66Lr)*s1*s2,
(c13Lr+c55Lr)*s1*s3LP + (c36Lr+c45Lr)*s2*s3LP],
[c16Lr*s1*s1 + c26Lr*s2*s2 + c45Lr*s3L2P + (c12Lr+c66Lr)*s1*s2,
c66Lr*s1*s1 + c22Lr*s2*s2 + c44Lr*s3L2P + 2*c26Lr*s1*s2,
(c36Lr+c45Lr)*s1*s3LP + (c23Lr+c44Lr)*s2*s3LP],
[(c13Lr+c55Lr)*s1*s3LP + (c36Lr+c45Lr)*s2*s3LP,
(c36Lr+c45Lr)*s1*s3LP + (c23Lr+c44Lr)*s2*s3LP,
c55Lr*s1*s1+c44Lr*s2*s2+c33Lr*s3L2P]])
eigLP = eigvec3symm(CMLP,rhoL)
temp1 = np.abs(np.sqrt( eigLP[0]*eigLP[0] + eigLP[1]*eigLP[1] + eigLP[2]*eigLP[2] ))
temp1 = (eigLP[0]*s1+eigLP[1]*s2)/np.abs(eigLP[0]*s1+eigLP[1]*s2) / temp1
eigLP = eigLP*temp1
# Form the velocity form of the Christoffel matrix for the transmitted S-wave
CMLS = np.array([[c11Lr*s1*s1 + c66Lr*s2*s2 + c55Lr*s3L2S + 2*c16Lr*s1*s2,
c16Lr*s1*s1 + c26Lr*s2*s2 + c45Lr*s3L2S + (c12Lr+c66Lr)*s1*s2,
(c13Lr+c55Lr)*s1*s3LS + (c36Lr+c45Lr)*s2*s3LS],
[c16Lr*s1*s1 + c26Lr*s2*s2 + c45Lr*s3L2S + (c12Lr+c66Lr)*s1*s2,
c66Lr*s1*s1 + c22Lr*s2*s2 + c44Lr*s3L2S + 2*c26Lr*s1*s2,
(c36Lr+c45Lr)*s1*s3LS + (c23Lr+c44Lr)*s2*s3LS],
[(c13Lr+c55Lr)*s1*s3LS + (c36Lr+c45Lr)*s2*s3LS,
(c36Lr+c45Lr)*s1*s3LS + (c23Lr+c44Lr)*s2*s3LS,
c55Lr*s1*s1+c44Lr*s2*s2+c33Lr*s3L2S]])
eigLS = eigvec3symm(CMLS,rhoL)
temp1 = np.abs(np.sqrt( eigLS[0]*eigLS[0] + eigLS[1]*eigLS[1] + eigLS[2]*eigLS[2] ))
temp1 = (eigLS[0]*s1+eigLS[1]*s2)/np.abs(eigLS[0]*s1+eigLS[1]*s2) / temp1
eigLS = eigLS*temp1
# Form the velocity form of the Christoffel matrix for the transmitted T-wave
CMLT = np.array([[c11Lr*s1*s1 + c66Lr*s2*s2 + c55Lr*s3L2T + 2*c16Lr*s1*s2,
c16Lr*s1*s1 + c26Lr*s2*s2 + c45Lr*s3L2T + (c12Lr+c66Lr)*s1*s2,
(c13Lr+c55Lr)*s1*s3LT + (c36Lr+c45Lr)*s2*s3LT],
[c16Lr*s1*s1 + c26Lr*s2*s2 + c45Lr*s3L2T + (c12Lr+c66Lr)*s1*s2,
c66Lr*s1*s1 + c22Lr*s2*s2 + c44Lr*s3L2T + 2*c26Lr*s1*s2,
(c36Lr+c45Lr)*s1*s3LT + (c23Lr+c44Lr)*s2*s3LT],
[(c13Lr+c55Lr)*s1*s3LT + (c36Lr+c45Lr)*s2*s3LT,
(c36Lr+c45Lr)*s1*s3LT + (c23Lr+c44Lr)*s2*s3LT,
c55Lr*s1*s1+c44Lr*s2*s2+c33Lr*s3L2T]])
eigLT = eigvec3symm(CMLT,rhoL)
temp1 = np.abs(np.sqrt( eigLT[0]*eigLT[0] + eigLT[1]*eigLT[1] + eigLT[2]*eigLT[2] ))
temp1 = (eigLT[0]*s1+eigLT[1]*s2)/np.abs(eigLT[0]*s1+eigLT[1]*s2) / temp1
eigLT = eigLT*temp1
# try and match up qSV and qSH with the most appropriate eigenvalues/eigenvectors
temp1 = eigLS[0]*n3 + eigLS[1]*n2 - eigLS[2]*n1
temp2 = eigLT[0]*n3 + eigLT[1]*n2 - eigLT[2]*n1
if (temp2 > temp1):
temp = s3LS
s3LS = s3LT
s3LT = temp
temp3 = eigLS
eigLS = eigLT
eigLT = temp3
print('eigUS')
print(eigUS)
print('eigUT')
print(eigUT)
print('eigLP')
print(eigLP)
print('eigLS')
print(eigLS)
print('eigLT')
print(eigLT)
# Construct X and Y matrices
XU = np.zeros(shape=(3, 3))
XU[0][0] = eigUP[0]
XU[0][1] = eigUS[0]
XU[0][2] = eigUT[0]
XU[1][0] = eigUP[1]
XU[1][1] = eigUS[1]
XU[1][2] = eigUT[1]
XU[2][0] = -(c13Ur*eigUP[0]+c36Ur*eigUP[1])*s1-(c23Ur*eigUP[1]+c36Ur*eigUP[0])*s2-c33Ur*eigUP[2]*s3UP
XU[2][1] = -(c13Ur*eigUS[0]+c36Ur*eigUS[1])*s1-(c23Ur*eigUS[1]+c36Ur*eigUS[0])*s2-c33Ur*eigUS[2]*s3US
XU[2][2] = -(c13Ur*eigUT[0]+c36Ur*eigUT[1])*s1-(c23Ur*eigUT[1]+c36Ur*eigUT[0])*s2-c33Ur*eigUT[2]*s3UT
YU = np.zeros(shape=(3, 3))
YU[0][0] = -(c55Ur*s1+c45Ur*s2)*eigUP[2]-(c55Ur*eigUP[0]+c45Ur*eigUP[1])*s3UP
YU[0][1] = -(c55Ur*s1+c45Ur*s2)*eigUS[2]-(c55Ur*eigUS[0]+c45Ur*eigUS[1])*s3US
YU[0][2] = -(c55Ur*s1+c45Ur*s2)*eigUT[2]-(c55Ur*eigUT[0]+c45Ur*eigUT[1])*s3UT
YU[1][0] = -(c45Ur*s1+c44Ur*s2)*eigUP[2]-(c45Ur*eigUP[0]+c44Ur*eigUP[1])*s3UP
YU[1][1] = -(c45Ur*s1+c44Ur*s2)*eigUS[2]-(c45Ur*eigUS[0]+c44Ur*eigUS[1])*s3US
YU[1][2] = -(c45Ur*s1+c44Ur*s2)*eigUT[2]-(c45Ur*eigUT[0]+c44Ur*eigUT[1])*s3UT
YU[2][0] = eigUP[2]
YU[2][1] = eigUS[2]
YU[2][2] = eigUT[2]
XL = np.zeros(shape=(3, 3))
XL[0][0] = eigLP[0]
XL[0][1] = eigLS[0]
XL[0][2] = eigLT[0]
XL[1][0] = eigLP[1]
XL[1][1] = eigLS[1]
XL[1][2] = eigLT[1]
XL[2][0] = -(c13Lr*eigLP[0]+c36Lr*eigLP[1])*s1-(c23Lr*eigLP[1]+c36Lr*eigLP[0])*s2-c33Lr*eigLP[2]*s3LP
XL[2][1] = -(c13Lr*eigLS[0]+c36Lr*eigLS[1])*s1-(c23Lr*eigLS[1]+c36Lr*eigLS[0])*s2-c33Lr*eigLS[2]*s3LS
XL[2][2] = -(c13Lr*eigLT[0]+c36Lr*eigLT[1])*s1-(c23Lr*eigLT[1]+c36Lr*eigLT[0])*s2-c33Lr*eigLT[2]*s3LT
YL = np.zeros(shape=(3, 3))
YL[0][0] = -(c55Lr*s1+c45Lr*s2)*eigLP[2]-(c55Lr*eigLP[0]+c45Lr*eigLP[1])*s3LP
YL[0][1] = -(c55Lr*s1+c45Lr*s2)*eigLS[2]-(c55Lr*eigLS[0]+c45Lr*eigLS[1])*s3LS
YL[0][2] = -(c55Lr*s1+c45Lr*s2)*eigLT[2]-(c55Lr*eigLT[0]+c45Lr*eigLT[1])*s3LT
YL[1][0] = -(c45Lr*s1+c44Lr*s2)*eigLP[2]-(c45Lr*eigLP[0]+c44Lr*eigLP[1])*s3LP
YL[1][1] = -(c45Lr*s1+c44Lr*s2)*eigLS[2]-(c45Lr*eigLS[0]+c44Lr*eigLS[1])*s3LS
YL[1][2] = -(c45Lr*s1+c44Lr*s2)*eigLT[2]-(c45Lr*eigLT[0]+c44Lr*eigLT[1])*s3LT
YL[2][0] = eigLP[2]
YL[2][1] = eigLS[2]
YL[2][2] = eigLT[2]
iXU = np.linalg.inv(XU)
iYU = np.linalg.inv(YU)
XUL = iXU.dot(XL)
YUL = iYU.dot(YL)
XpY = XUL + YUL
XmY = XUL - YUL
Tran = np.linalg.inv(XpY)
Refl = XmY.dot(Tran)
Tran = 2.0*Tran
print('#### Impedance matrices')
print('XU')
print(XU)
print()
print('YU')
print(YU)
print()
print('XL')
print(XL)
print()
print('YL')
print(YL)
print()
return(Refl[0][0])
thetas = 1
phi = 60
p1 = 2600
chi1 = 0
C1 = Cij(0, 0, 0, p1, 2260, 1428)
p2 = 2700
chi2 = 0
C2 = Cij(0.05, 0.02, 0.1, p1, 2370, 1360)
Rpp = ortho_Rpp(C1, C2, p1, p2, chi1, chi2, phi, theta)
print('***Rpp***')
print(Rpp)
In [ ]: