3.1.8. MTfit.convertΒΆ
"""
moment_tensor_conversion.py
***************************
Module containing moment tensor conversion functions. Acts on the parameters of the moment tensor 3x3 form or the modified 6-vector form, dependent on the name
The function naming is OriginalVariables_NewVariables
The coordinate system is North (X), East (Y), Down (Z)
"""
# **Restricted: For Non-Commercial Use Only**
# This code is protected intellectual property and is available solely for teaching
# and non-commercially funded academic research purposes.
#
# Applications for commercial use should be made to Schlumberger or the University of Cambridge.
import logging
import numpy as np
from scipy.optimize import fsolve
from ..utilities import C_EXTENSION_FALLBACK_LOG_MSG
logger = logging.getLogger('MTfit.convert')
try:
from . import cmoment_tensor_conversion
except ImportError:
cmoment_tensor_conversion = None
except Exception:
logger.exception('Error importing c extension')
cmoment_tensor_conversion = None
def MT33_MT6(MT33):
"""
Convert a 3x3 matrix to six vector maintaining normalisation. 6-vector has the form::
Mxx
Myy
Mzz
sqrt(2)*Mxy
sqrt(2)*Mxz
sqrt(2)*Myz
Args
M33: 3x3 numpy matrix
Returns
numpy.matrix: MT 6-vector
"""
MT6 = np.matrix([[MT33[0, 0]], [MT33[1, 1]], [MT33[2, 2]], [np.sqrt(2)*MT33[0, 1]],
[np.sqrt(2)*MT33[0, 2]], [np.sqrt(2)*MT33[1, 2]]])
MT6 = np.matrix(MT6/np.sqrt(np.sum(np.multiply(MT6, MT6), axis=0)))
return MT6
def MT6_MT33(MT6):
"""
Convert a six vector to a 3x3 MT maintaining normalisation. 6-vector has the form::
Mxx
Myy
Mzz
sqrt(2)*Mxy
sqrt(2)*Mxz
sqrt(2)*Myz
Args
MT6: numpy matrix Moment tensor 6-vector
Returns
numpy.matrix: 3x3 Moment Tensor
"""
if np.prod(MT6.shape) != 6:
raise ValueError("Input MT must be 6 vector not {}".format(MT6.shape))
if len(MT6.shape) == 1:
MT6 = np.matrix(MT6)
if len(MT6.shape) == 2 and MT6.shape[1] == 6:
MT6 = MT6.T
return np.matrix([[MT6[0, 0], (1/np.sqrt(2))*MT6[3, 0], (1/np.sqrt(2))*MT6[4, 0]],
[(1/np.sqrt(2))*MT6[3, 0], MT6[1, 0],
(1/np.sqrt(2))*MT6[5, 0]],
[(1/np.sqrt(2))*MT6[4, 0], (1/np.sqrt(2))*MT6[5, 0], MT6[2, 0]]])
def MT6_TNPE(MT6):
"""
Convert the 6xn Moment Tensor to the T,N,P vectors and the eigenvalues.
Args
MT6: 6xn numpy matrix
Returns
(numpy.matrix, numpy.matrix, numpy.matrix, numpy.array): tuple of T, N, P
vectors and Eigenvalue array
"""
if cmoment_tensor_conversion:
try:
if not isinstance(MT6, np.ndarray):
MT6 = np.array(MT6)
if MT6.ndim < 2:
MT6 = np.array([MT6])
if MT6.shape[1] == 6 and MT6.shape[0] != 6:
MT6 = MT6.T
return cmoment_tensor_conversion.MT6_TNPE(MT6.astype(np.float64))
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
try:
n = MT6.shape[1]
except Exception:
MT6 = np.array([MT6]).T
n = MT6.shape[1]
T = np.matrix(np.empty((3, n)))
N = np.matrix(np.empty((3, n)))
P = np.matrix(np.empty((3, n)))
E = np.empty((3, n))
for i in range(n):
T[:, i], N[:, i], P[:, i], E[:, i] = MT33_TNPE(MT6_MT33(MT6[:, i]))
return T, N, P, E
def MT6_Tape(MT6):
"""
Convert the moment tensor 6-vector to the Tape parameters.
6-vector has the form::
Mxx
Myy
Mzz
sqrt(2)*Mxy
sqrt(2)*Mxz
sqrt(2)*Myz
Args
MT6: numpy matrix six-vector
Returns
(numpy.array, numpy.array, numpy.array, numpy.array, numpy.array): tuple of
gamma, delta, strike, cos(dip) and slip (angles in radians)
"""
if len(MT6.shape) > 1 and MT6.shape[1] > 1:
gamma = np.array(np.empty((MT6.shape[1],)))
delta = np.array(np.empty((MT6.shape[1],)))
kappa = np.array(np.empty((MT6.shape[1],)))
h = np.array(np.empty((MT6.shape[1],)))
sigma = np.array(np.empty((MT6.shape[1],)))
for i in range(MT6.shape[1]):
gamma[i], delta[i], kappa[i], h[i], sigma[i] = MT6_Tape(MT6[:, i])
return gamma, delta, kappa, h, sigma
MT33 = MT6_MT33(MT6)
T, N, P, E = MT33_TNPE(MT33)
gamma, delta = E_GD(E)
kappa, dip, sigma = TNP_SDR(T, N, P)
if np.abs(sigma) > np.pi/2:
kappa, dip, sigma = SDR_SDR(kappa, dip, sigma)
h = np.cos(dip)
return (np.array(gamma).flatten(), np.array(delta).flatten(),
np.array(kappa).flatten(), np.array(h).flatten(),
np.array(sigma).flatten())
def MT33_TNPE(MT33):
"""
Convert the 3x3 Moment Tensor to the T,N,P vectors and the eigenvalues.
Args
MT33: 3x3 numpy matrix
Returns
(numpy.matrix, numpy.matrix, numpy.matrix, numpy.array): tuple of T, N, P
vectors and Eigenvalue array
"""
E, L = np.linalg.eig(MT33)
idx = E.argsort()[::-1]
E = E[idx]
L = L[:, idx]
T = L[:, 0]
P = L[:, 2]
N = L[:, 1]
return (T, N, P, E)
def MT33_SDR(MT33):
"""
Convert the 3x3 Moment Tensor to the strike, dip and rake.
Args
MT33: 3x3 numpy matrix
Returns
(float, float, float): tuple of strike, dip, rake angles in radians
"""
T, N, P, E = MT33_TNPE(MT33)
N1, N2 = TP_FP(T, P)
return FP_SDR(N1, N2)
def MT33_GD(MT33):
"""
Convert the 3x3 Moment Tensor to theTape parameterisation gamma and delta.
Args
MT33: 3x3 numpy matrix
Returns
(numpy.array, numpy.array): tuple of gamma, delta
"""
E, L = np.linalg.eig(MT33)
return E_GD(E)
def E_tk(E):
"""
Convert the moment tensor eigenvalues to the Hudson tau, k parameters
Args
E: indexable list/array (e.g numpy.array) of moment tensor eigenvalues
Returns
(float, float): tau, k tuple
"""
if isinstance(E, np.ndarray):
E = np.squeeze(np.array(E))
if len(E.shape) > 1 and E.shape[1] > 1:
tau = np.empty(E.shape[1])
k = np.empty(E.shape[1])
for i in range(E.shape[1]):
tau[i], k[i] = E_tk(E[:, i])
else:
idx = [0, 2, 1] # Odd sorting
E = E[idx]
iso = (E[0]+E[1]+E[2])/3
dev0 = E[0]-iso
dev1 = E[1]-iso
dev2 = E[2]-iso
if dev2 > 0:
k = iso/(abs(iso)-dev1)
T = -2*(dev2)/(dev1)
elif dev2 < 0:
k = iso/(abs(iso)+dev0)
T = 2*dev2/dev0
else:
k = iso/(abs(iso)+dev0)
T = 0
tau = T*(1-abs(k))
return tau, k
def tk_uv(tau, k):
"""
Convert the Hudson tau, k parameters to the Hudson u, v parameters
Args
tau: float, Hudson tau parameter
k: float, Hudson k parameter
Returns
(float, float): u, v tuple
"""
try:
u = tau.copy()
v = k.copy()
idx = (tau > 0)*(k > 0)*(tau < 4*k)
u[idx] = tau[idx]/(1-(tau[idx]/2))
v[idx] = k[idx]/(1-(tau[idx]/2))
idx = (tau > 0)*(k > 0)*~(tau < 4*k)
u[idx] = tau[idx]/(1-2*k[idx])
v[idx] = k[idx]/(1-2*k[idx])
idx = (tau < 0)*(k < 0)*(tau > 4*k)
u[idx] = tau[idx]/(1+(tau[idx]/2))
v[idx] = k[idx]/(1+(tau[idx]/2))
idx = (tau < 0)*(k < 0)*~(tau > 4*k)
u[idx] = tau[idx]/(1+2*k[idx])
v[idx] = k[idx]/(1+2*k[idx])
except Exception:
if tau > 0 and k > 0:
if tau < 4*k:
u = tau/(1-(tau/2))
v = k/(1-(tau/2))
else:
u = tau/(1-2*k)
v = k/(1-2*k)
elif tau < 0 and k < 0:
if tau > 4*k:
u = tau/(1+(tau/2))
v = k/(1+(tau/2))
else:
u = tau/(1+2*k)
v = k/(1+2*k)
else:
u = tau
v = k
return u, v
def E_uv(E):
"""
Convert the eigenvalues to the Hudson i, v parameters
Args
E: indexable list/array (e.g numpy.array) of moment tensor eigenvalues
Returns
(float, float): u, v tuple
"""
return tk_uv(*E_tk(E))
def E_GD(E):
"""
Convert the eigenvalues to the Tape parameterisation gamma and delta.
Args
E: array of eigenvalues
Returns
(numpy.array, numpy.array): tuple of gamma, delta
"""
if cmoment_tensor_conversion:
try:
# Expect E to be an array of arrays
if E.ndim < 2:
E = np.array([E])
if E.shape[1] == 3 and E.shape != (3, 3):
E = E.T
return cmoment_tensor_conversion.E_GD(E.astype(np.float64))
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
if isinstance(E, np.ndarray):
E = np.squeeze(np.array(E))
if len(E.shape) > 1 and E.shape[1] > 1:
gamma = np.empty(E.shape[1])
delta = np.empty(E.shape[1])
for i in range(E.shape[1]):
gamma[i], delta[i] = E_GD(E[:, i])
else:
if len(np.unique(E)) == 1:
return 0, np.sign(np.unique(E))*np.pi/2
idx = E.argsort()[::-1]
E = E[idx]
gamma = np.arctan2(-E[0]+2*E[1]-E[2], np.sqrt(3)*(E[0]-E[2]))
beta = np.arccos(np.sum(E)/(np.sqrt(3)*np.sqrt(np.sum(E*E))))
delta = np.real(np.pi/2-beta)
return gamma, delta
def GD_basic_cdc(gamma, delta):
"""
Convert gamma, delta to basic crack+double-couple parameters
Gamma and delta are the source type parameters from the Tape parameterisation.
Args
gamma: numpy array of gamma values
delta: numpy array of delta values
Returns:
(numpy.array, numpy.array): tuple of alpha, poisson
"""
alpha = np.arccos(-np.sqrt(3)*np.tan(gamma))
poisson = (1+np.sqrt(2) * (np.tan((np.pi/2)-delta) * np.sin(gamma)))
poisson /= (2-(np.sqrt(2) * (np.tan((np.pi/2)-delta) * np.sin(gamma))))
return alpha, poisson
def TNP_SDR(T, N, P):
"""
Convert the T,N,P vectors to the strike, dip and rake in radians
Args
T: numpy matrix of T vectors.
N: numpy matrix of N vectors.
P: numpy matrix of P vectors.
Returns
(float, float, float): tuple of strike, dip and rake angles of fault plane in radians
"""
if cmoment_tensor_conversion:
try:
return cmoment_tensor_conversion.TP_SDR(T.astype(np.float64), P.astype(np.float64))
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
N1, N2 = TP_FP(T, P)
return FP_SDR(N1, N2)
def TP_FP(T, P):
"""
Convert the 3x3 Moment Tensor to the fault normal and slip vectors.
Args
T: numpy matrix of T vectors.
P: numpy matrix of P vectors.
Returns
(numpy.matrix, numpy.matrix): tuple of Normal and slip vectors
"""
if T.ndim == 1:
T = np.matrix(T)
if P.ndim == 1:
P = np.matrix(P)
if T.shape[0] != 3:
T = T.T
if P.shape[0] != 3:
P = P.T
TP1 = T+P
TP2 = T-P
N1 = (TP1)/np.sqrt(np.einsum('ij,ij->j', TP1, TP1))
N2 = (TP2)/np.sqrt(np.einsum('ij,ij->j', TP2, TP2))
return (N1, N2)
def FP_SDR(normal, slip):
"""
Convert fault normal and slip to strike, dip and rake
Coordinate system is North East Down.
Args
normal: numpy matrix - Normal vector
slip: numpy matrix - Slip vector
Returns
(float, float, float): tuple of strike, dip and rake angles in radians
"""
if not isinstance(slip, np.matrixlib.defmatrix.matrix):
slip = slip/np.sqrt(np.sum(slip*slip, axis=0))
else:
# Do we need to replace this with einsum
slip = slip/np.sqrt(np.einsum('ij,ij->j', slip, slip))
if not isinstance(normal, np.matrixlib.defmatrix.matrix):
normal = normal/np.sqrt(np.sum(normal*normal, axis=0))
else:
normal = normal/np.sqrt(np.einsum('ij,ij->j', normal, normal))
slip[:, np.array(normal[2, :] > 0).flatten()] *= -1
normal[:, np.array(normal[2, :] > 0).flatten()] *= -1
normal = np.array(normal)
slip = np.array(slip)
strike, dip = normal_SD(normal)
rake = np.arctan2(-slip[2], slip[0]*normal[1]-slip[1]*normal[0])
strike[dip > np.pi/2] += np.pi
rake[dip > np.pi/2] = 2*np.pi-rake[dip > np.pi/2]
dip[dip > np.pi/2] = np.pi-dip[dip > np.pi/2]
strike = np.mod(strike, 2*np.pi)
rake[rake > np.pi] -= 2*np.pi
rake[rake < -np.pi] += 2*np.pi
return (np.array(strike).flatten(), np.array(dip).flatten(), np.array(rake).flatten())
def basic_cdc_GD(alpha, poisson=0.25):
"""
Convert alpha and poisson ratio to gamma and delta
alpha is opening angle, poisson : ratio lambda/(2(lambda+mu)) Defaults to 0.25.
Uses basic crack+double-couple model of Minson et al (Seismically and geodetically
determined nondouble-couple source mechanisms from the 2000 Miyakejima volcanic
earthquake swarm, Minson et al, 2007, JGR 112) and Tape and Tape (The classical
model for moment tensors, Tape and Tape, 2013, GJI)
Args
alpha: Opening angle in radians (between 0 and pi/2)
poisson:[0.25] Poisson ratio on the fault surface.
Returns:
(numpy.array, numpy.array): tuple of gamma, delta
"""
gamma = np.arctan((-1/np.sqrt(3.))*np.cos(alpha))
beta = np.arccos(np.sqrt(2./3)*np.cos(alpha)*(1.+poisson) /
np.sqrt(((1.-2.*poisson)*(1.-2.*poisson)) +
np.cos(alpha)*np.cos(alpha)*(1.+2.*poisson*poisson)))
try:
# fix for odd cos division problems
beta[alpha == np.pi/2] = np.pi/2
except Exception:
if alpha == np.pi/2:
beta = np.pi/2
delta = np.pi/2-beta
return (gamma, delta)
def GD_E(gamma, delta):
"""
Convert the Tape parameterisation gamma and delta to the eigenvalues.
Args
gamma: numpy array of gamma values
delta: numpy array of delta values
Returns
numpy.array: array of eigenvalues
"""
U = (1/np.sqrt(6))*np.matrix([[np.sqrt(3), 0, -np.sqrt(3)],
[-1, 2, -1],
[np.sqrt(2), np.sqrt(2), np.sqrt(2)]])
gamma = np.array(gamma).flatten()
delta = np.array(delta).flatten()
X = np.matrix([np.multiply(np.cos(gamma), np.sin((np.pi/2)-delta)),
np.multiply(np.sin(gamma), np.sin((np.pi/2)-delta)),
np.cos((np.pi/2)-delta)]).T
if X.shape[1] == 3:
X = X.T
return np.array((U.T*X))
def SDR_TNP(strike, dip, rake):
"""
Convert strike, dip rake to TNP vectors
Coordinate system is North East Down.
Args
strike: float radians
dip: float radians
rake: float radians
Returns
(numpy.matrix, numpy.matrix, numpy.matrix): tuple of T,N,P vectors.
"""
strike = np.array(strike).flatten()
dip = np.array(dip).flatten()
rake = np.array(rake).flatten()
N1 = np.matrix([(np.cos(strike)*np.cos(rake))+(np.sin(strike)*np.cos(dip)*np.sin(rake)),
(np.sin(strike)*np.cos(rake)) -
np.cos(strike)*np.cos(dip)*np.sin(rake),
-np.sin(dip)*np.sin(rake)])
N2 = np.matrix([-np.sin(strike)*np.sin(dip), np.cos(strike)*np.sin(dip), -np.cos(dip)])
return FP_TNP(N1, N2)
def SDR_SDR(strike, dip, rake):
"""
Convert strike, dip rake to strike, dip rake for other fault plane
Coordinate system is North East Down.
Args
strike: float radians
dip: float radians
rake: float radians
Returns
(float, float, float): tuple of strike, dip and rake angles of alternate fault
plane in radians
"""
if cmoment_tensor_conversion:
try:
return cmoment_tensor_conversion.SDR_SDR(strike.astype(np.float64), dip.astype(np.float64), rake.astype(np.float64))
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
# Handle multiple inputs
N1, N2 = SDR_FP(strike, dip, rake)
s1, d1, r1 = FP_SDR(N1, N2)
s2, d2, r2 = FP_SDR(N2, N1)
# This should be ok to return s2,d2,r2 but doesn't seem to work
try:
r2[np.abs(strike-s2) < 1] = r1[np.abs(strike-s2) < 1]
d2[np.abs(strike-s2) < 1] = d1[np.abs(strike-s2) < 1]
s2[np.abs(strike-s2) < 1] = s1[np.abs(strike-s2) < 1]
return s2, d2, r2
except Exception:
if np.abs(strike-s1) < 1:
return (s2, d2, r2)
else:
return (s1, d1, r1)
def FP_TNP(normal, slip):
"""
Convert fault normal and slip to TNP axes
Coordinate system is North East Down.
Args
normal: numpy matrix - normal vector
slip: numpy matrix - slip vector
Returns
(numpy.matrix, numpy.matrix, numpy.matrix): tuple of T, N, P vectors
"""
T = (normal+slip)
T = T/np.sqrt(np.einsum('ij,ij->j', T, T))
P = (normal-slip)
P = P/np.sqrt(np.einsum('ij,ij->j', P, P))
N = np.matrix(-np.cross(T.T, P.T)).T
return (T, N, P)
def SDSD_FP(strike1, dip1, strike2, dip2):
"""
Convert strike and dip pairs to fault normal and slip
Converts the strike and dip pairs in radians to the fault normal and slip.
Args
strike1: float strike angle of fault plane 1 in radians
dip1: float dip angle of fault plane 1 in radians
strike2: float strike angle of fault plane 2 in radians
dip2: float dip of fault plane 2 in radians
Returns
(numpy.matrix, numpy.matrix): tuple of Normal and slip vectors
"""
strike1 = np.array(strike1).flatten()
dip1 = np.array(dip1).flatten()
strike2 = np.array(strike2).flatten()
dip2 = np.array(dip2).flatten()
N1 = np.matrix([-np.sin(strike2)*np.sin(dip2),
np.cos(strike2)*np.sin(dip2),
-np.cos(dip2)])
N2 = np.matrix([-np.sin(strike1)*np.sin(dip1),
np.cos(strike1)*np.sin(dip1),
-np.cos(dip1)])
return (N1, N2)
def SDR_FP(strike, dip, rake):
"""
Convert the strike, dip and rake in radians to the fault normal and slip.
Args
strike: float strike angle of fault plane in radians
dip: float dip angle of fault plane in radians
rake: float rake angle of fault plane in radians
Returns
(numpy.matrix, numpy.matrix): tuple of Normal and slip vectors
"""
T, N, P = SDR_TNP(strike, dip, rake)
return TP_FP(T, P)
def SDR_SDSD(strike, dip, rake):
"""
Convert the strike, dip and rake to the strike and dip pairs (all angles in radians).
Args
strike: float strike angle of fault plane in radians
dip: float dip angle of fault plane in radians
rake: float rake angle of fault plane in radians
Returns
(float, float, float, float): tuple of strike1, dip1, strike2, dip2 angles in radians
"""
N1, N2 = SDR_FP(strike, dip, rake)
return FP_SDSD(N1, N2)
def FP_SDSD(N1, N2):
"""
Convert the the fault normal and slip vectors to the strike and dip pairs
(all angles in radians).
Args
Normal: numpy matrix - Normal vector
Slip: numpy matrix - Slip vector
Returns
(float, float, float, float): tuple of strike1, dip1, strike2, dip2 angles
in radians
"""
s1, d1 = normal_SD(N1)
s2, d2 = normal_SD(N2)
return (s1, d1, s2, d2)
def Tape_MT33(gamma, delta, kappa, h, sigma, **kwargs):
"""
Convert Tape parameters to a 3x3 moment tensor
Args
gamma: float, gamma parameter (longitude on funamental lune takes values
between -pi/6 and pi/6).
delta: float, delta parameter (latitude on funamental lune takes values
between -pi/2 and pi/2)
kappa: float, strike (takes values between 0 and 2*pi)
h: float, cos(dip) (takes values between 0 and 1)
sigma: float, slip angle (takes values between -pi/2 and pi/2)
Returns
numpy.matrix: 3x3 moment tensor
"""
E = GD_E(gamma, delta)
D = np.diag(E.flatten())
T, N, P = SDR_TNP(kappa, np.arccos(h), sigma)
L = np.matrix(
[np.array(T).flatten(), np.array(N).flatten(), np.array(P).flatten()]).T
MT33 = L*D*L.T
return MT33
def Tape_MT6(gamma, delta, kappa, h, sigma):
"""
Convert the Tape parameterisation to the moment tensor six-vectors.
Args
gamma: Gamma parameter (longitude on funamental lune takes values
between -pi/6 and pi/6).
delta: Delta parameter (latitude on funamental lune takes values
between -pi/2 and pi/2)
kappa: Strike (takes values between 0 and 2*pi)
h: Cos(dip) (takes values between 0 and 1)
sigma: Slip angle (takes values between -pi/2 and pi/2)
Returns
np.array: Array of MT 6-vectors
"""
if cmoment_tensor_conversion:
try:
return cmoment_tensor_conversion.Tape_MT6(gamma.astype(np.float64),
delta.astype(np.float64),
kappa.astype(np.float64),
h.astype(np.float64),
sigma.astype(np.float64))
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
if isinstance(gamma, (list, np.ndarray)):
MT6 = np.empty((6, len(gamma)))
for i in range(len(gamma)):
MT6[:, i] = np.squeeze(MT33_MT6(Tape_MT33(gamma[i],
delta[i],
kappa[i],
h[i],
sigma[i])))
return MT6
MT33 = Tape_MT33(gamma, delta, kappa, h, sigma)
return MT33_MT6(MT33)
def Tape_TNPE(gamma, delta, kappa, h, sigma):
"""
Convert the Tape parameterisation to the T,N,P vectors and the eigenvalues.
Args
gamma: Gamma parameter (longitude on funamental lune takes values
between -pi/6 and pi/6).
delta: Delta parameter (latitude on funamental lune takes values
between -pi/2 and pi/2)
kappa: Strike (takes values between 0 and 2*pi)
h: Cos(dip) (takes values between 0 and 1)
sigma: Slip angle (takes values between -pi/2 and pi/2)
Returns
(numpy.matrix, numpy.matrix, numpy.matrix, numpy.array): T,N,P vectors
and Eigenvalues tuple
"""
E = GD_E(gamma, delta)
T, N, P = SDR_TNP(kappa, np.arccos(h), sigma)
return (T, N, P, E)
def normal_SD(normal):
"""
Convert a plane normal to strike and dip
Coordinate system is North East Down.
Args
normal: numpy matrix - Normal vector
Returns
(float, float): tuple of strike and dip angles in radians
"""
if not isinstance(normal, np.matrixlib.defmatrix.matrix):
normal = np.array(normal)/np.sqrt(np.sum(normal*normal, axis=0))
else:
normal = normal/np.sqrt(np.diag(normal.T*normal))
normal[:, np.array(normal[2, :] > 0).flatten()] *= -1
normal = np.array(normal)
strike = np.arctan2(-normal[0], normal[1])
dip = np.arctan2((normal[1]**2+normal[0]**2),
np.sqrt((normal[0]*normal[2])**2+(normal[1]*normal[2])**2))
strike = np.mod(strike, 2*np.pi)
return strike, dip
def toa_vec(azimuth, plunge, radians=False):
"""
Convert the azimuth and plunge of a vector to a cartesian description of the vector
Args
azimuth: float, vector azimuth
plunge: float, vector plunge
Keyword Arguments
radians: boolean, flag to use radians [default = False]
Returns
np.matrix: vector
"""
if not radians:
azimuth = np.pi*np.array(azimuth)/180.
plunge = np.pi*np.array(plunge)/180.
if not isinstance(plunge, np.ndarray):
plunge = np.array([plunge])
try:
return np.matrix([np.cos(azimuth)*np.sin(plunge),
np.sin(azimuth)*np.sin(plunge),
np.cos(plunge)])
except Exception:
return np.array([np.cos(azimuth)*np.sin(plunge),
np.sin(azimuth)*np.sin(plunge),
np.cos(plunge)])
def output_convert(mts):
"""
Convert the moment tensors into several different parameterisations
The moment tensor six-vectors are converted into the Tape gamma,delta,kappa,h,sigma
parameterisation; the Hudson u,v parameterisation; and the strike, dip and rakes
of the two fault planes are calculated.
Args
mts: numpy array of moment tensor six-vectors
Returns
dict: dictionary of numpy arrays for each parameter
"""
if cmoment_tensor_conversion:
try:
return cmoment_tensor_conversion.MT_output_convert(mts.astype(np.float64))
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
g = np.empty((mts.shape[1],))
d = np.empty((mts.shape[1],))
k = np.empty((mts.shape[1],))
h = np.empty((mts.shape[1],))
s = np.empty((mts.shape[1],))
u = np.empty((mts.shape[1],))
v = np.empty((mts.shape[1],))
s1 = np.empty((mts.shape[1],))
d1 = np.empty((mts.shape[1],))
r1 = np.empty((mts.shape[1],))
s2 = np.empty((mts.shape[1],))
d2 = np.empty((mts.shape[1],))
r2 = np.empty((mts.shape[1],))
rad_cor = 180/np.pi
for i in range(mts.shape[1]):
MT33 = MT6_MT33(mts[:, i])
T, N, P, E = MT33_TNPE(MT33)
u[i], v[i] = tk_uv(*E_tk(E))
g[i], d[i] = E_GD(E)
k[i], dip, s[i] = TNP_SDR(T, N, P)
s1[i] = k[i]
d1[i] = dip
r1[i] = s[i]
if np.abs(s[i]) > np.pi/2:
k[i], dip, s[i] = SDR_SDR(k[i], dip, s[i])
h[i] = np.cos(dip)
s2[i], d2[i], r2[i] = SDR_SDR(s1[i], d1[i], r1[i])
s1[i] *= rad_cor
d1[i] *= rad_cor
r1[i] *= rad_cor
s2[i] *= rad_cor
d2[i] *= rad_cor
r2[i] *= rad_cor
g = np.array(g).flatten()
d = np.array(d).flatten()
k = np.array(k).flatten()
h = np.array(h).flatten()
s = np.array(s).flatten()
u = np.array(u).flatten()
v = np.array(v).flatten()
s1 = np.array(s1).flatten()
d1 = np.array(d1).flatten()
r1 = np.array(r1).flatten()
s2 = np.array(s2).flatten()
d2 = np.array(d2).flatten()
r2 = np.array(r2).flatten()
return {'g': g, 'd': d, 'k': k, 'h': h, 's': s, 'u': u, 'v': v,
'S1': s1, 'D1': d1, 'R1': r1, 'S2': s2, 'D2': d2, 'R2': r2}
# Bi-axes
def isotropic_c(lambda_=1, mu=1, c=False):
"""
Calculate the isotropic stiffness tensor
Calculate isotropic stiffness parameters. The input parameters are either
the two lame parameters lambda and mu, or is a full 21 element stiffness tensor::
C(21) = ( C_11, C_12, C_13, C_14, C_15, C_16,
C_22, C_23, C_24, C_25, C_26,
C_33, C_34, C_35, C_36,
C_44, C_45, C_46,
C_55, C_56,
C_66 )
(upper triangular part of Voigt stiffness matrix)
If the full stiffness tensor is used, the "average" isotropic approximation is
calculated using Eqns 81a and 81b from Chapman, C and Leaney,S, 2011. A new
moment-tensor decomposition for seismic events in anisotropic media, GJI, 188(1),
343-370.
Args
lambda_: lambda value
mu: mu value
c: list or numpy array of the 21 element input stiffness tensor (overrides
lambda and mu arguments)
Returns
list: list of 21 elements of the stiffness tensor
"""
if cmoment_tensor_conversion:
try:
if isinstance(c, bool):
c = []
return cmoment_tensor_conversion.isotropic_c(lambda_, mu, c)
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
# Calculate isotropic approximation
if not isinstance(c, bool) and len(c) == 21:
# Eqns 81a and 81b from chapman and leaney 2011
mu = ((c[0]+c[6]+c[11])+3*(c[15]+c[18]+c[20])-(c[1]+c[2]+c[7]))/15
lambda_ = ((c[0]+c[6]+c[11])-2*(c[15]+c[18]+c[20])+4*(c[1]+c[2]+c[7]))/15
n = lambda_+2*mu
c = [n, lambda_, lambda_, 0, 0, 0, n, lambda_, 0, 0, 0, n, 0, 0, 0, mu, 0, 0, mu, 0, mu]
return c
def MT6_biaxes(MT6, c=isotropic_c(lambda_=1, mu=1)):
"""
Convert six vector to bi-axes
Convert the moment tensor 6-vector to the bi-axes decomposition from Chapman, C and
Leaney,S, 2011. A new moment-tensor decomposition for seismic events in anisotropic
media, GJI, 188(1), 343-370.
The 6-vector has the form::
Mxx
Myy
Mzz
sqrt(2)*Mxy
sqrt(2)*Mxz
sqrt(2)*Myz
The stiffness tensor can be provided as an input, as a list of the 21 elements of the
upper triangular part of the Voigt stiffness matrix::
C(21) = ( C_11, C_12, C_13, C_14, C_15, C_16,
C_22, C_23, C_24, C_25, C_26,
C_33, C_34, C_35, C_36,
C_44, C_45, C_46,
C_55, C_56,
C_66 )
(upper triangular part of Voigt stiffness matrix)
Alternatively, the default isotropic parameters can be used or a possible isotropic
stiffness tensor can be genereated using (isotropic_c)
Args
MT6: numpy matrix Moment tensor 6-vector
c: list or numpy array of the 21 element input stiffness tensor
Returns
(numpy.array, numpy.array, numpy.array): tuple of phi (bi-axes) vectors,
explosion value and area_displacement value.
"""
if isinstance(c, bool):
c = isotropic_c(1, 1)
if len(MT6.shape) > 1 and MT6.shape[1] > 1:
phi = np.empty((MT6.shape[1], 3, 2))
explosion = np.empty((MT6.shape[1],))
area_displacement = np.empty((MT6.shape[1],))
for i in range(MT6.shape[1]):
try:
phi[i, :, :], explosion[i], area_displacement[i] = MT6_biaxes(MT6[:, i], c[i])
except Exception:
phi[i, :, :], explosion[i], area_displacement[i] = MT6_biaxes(MT6[:, i], c)
else:
if np.prod(MT6.shape) != 6:
MT6 = np.asarray([MT6]).T
if cmoment_tensor_conversion:
try:
return cmoment_tensor_conversion.MT6_biaxes(MT6.flatten().astype(np.float64), c)
except Exception:
pass
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
lambda2mu = (
3*(c[0]+c[6]+c[11])+4*(c[15]+c[18]+c[20])+2*(c[1]+c[2]+c[7]))/15
mu = ((c[0]+c[6]+c[11])+3*(c[15]+c[18]+c[20])-(c[1]+c[2]+c[7]))/15
lambda_ = lambda2mu-2*mu
T, N, P, E = MT6_TNPE(MT6)
isotropic = (lambda_+mu)*E[1]/mu-lambda_*(E[0]+E[2])/(2*mu)
if is_isotropic_c(c):
explosion = isotropic
else:
def isotropic_solve(iso):
iso6 = np.squeeze(iso)*np.array([[1], [1], [1], [0], [0], [0]])
if iso6.shape != MT6.shape:
iso6 = iso6.T
T, N, P, E = MT6_TNPE(MT6c_D6(np.squeeze(MT6-iso6), c))
return E[1]
explosion = fsolve(isotropic_solve, isotropic)
explosion6 = np.squeeze(explosion)*np.array([[1], [1], [1], [0], [0], [0]])
if explosion6.shape != MT6.shape:
explosion6 = explosion6.T
T, N, P, E = MT6_TNPE(MT6c_D6(np.squeeze(MT6-explosion6), c))
area_displacement = E[0]-E[2]
phi = np.zeros((3, 2))
if area_displacement != 0: # to avoid undefined
cphi = np.squeeze(np.sqrt(E[0]/area_displacement))
sphi = np.squeeze(np.sqrt(-E[2]/area_displacement))
phi[:, 0] = np.array(cphi*T+sphi*P).flatten()
phi[:, 1] = np.array(cphi*T-sphi*P).flatten()
return phi, explosion, area_displacement
def MT6c_D6(MT6, c=isotropic_c(lambda_=1, mu=1)):
"""
Convert the moment tensor 6-vector to the potency tensor.
The 6-vector has the form::
Mxx
Myy
Mzz
sqrt(2)*Mxy
sqrt(2)*Mxz
sqrt(2)*Myz
The stiffness tensor can be provided as an input, as a list of the 21 elements of the
upper triangular part of the Voigt stiffness matrix::
C(21) = ( C_11, C_12, C_13, C_14, C_15, C_16,
C_22, C_23, C_24, C_25, C_26,
C_33, C_34, C_35, C_36,
C_44, C_45, C_46,
C_55, C_56,
C_66 )
(upper triangular part of Voigt stiffness matrix)
Alternatively, the default isotropic parameters can be used or a possible isotropic
stiffness tensor can be genereated using (isotropic_c).
Args
MT6: numpy matrix Moment tensor 6-vector
c: list or numpy array of the 21 element input stiffness tensor
Returns
numpy.array: numpy array of the potency 6 vector (in the same ordering as the
moment tensor six vector)
"""
if cmoment_tensor_conversion:
try:
return np.asarray(cmoment_tensor_conversion.MT6c_D6(MT6.astype(np.float64), c))
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
mtvoigt = MT6[np.array([0, 1, 2, 5, 4, 3])]
mtvoigt = np.matrix(mtvoigt)
if mtvoigt.shape[1] == 6:
mtvoigt = mtvoigt.T
# Convert to voigt
dvoigt = np.linalg.solve(np.matrix(c21_cvoigt(c)), mtvoigt)
return dvoigt[np.array([0, 1, 2, 5, 4, 3])]
def is_isotropic_c(c):
"""
Evaluate if an input stiffness tensor is isotropic
The stiffness tensor needs to be provided as a list of the 21 elements of the
upper triangular part of the Voigt stiffness matrix::
C(21) = ( C_11, C_12, C_13, C_14, C_15, C_16,
C_22, C_23, C_24, C_25, C_26,
C_33, C_34, C_35, C_36,
C_44, C_45, C_46,
C_55, C_56,
C_66 )
(upper triangular part of Voigt stiffness matrix)
Args
c: input list of stiffness parameters (21 required)
Returns
bool: result of is_isotropic check
"""
if cmoment_tensor_conversion:
try:
return cmoment_tensor_conversion.is_isotropic_c(c)
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
tol = 1.e-6*c_norm(c)
return ((abs(c[3]) < tol) and (abs(c[4]) < tol) and (abs(c[5]) < tol) and
(abs(c[8]) < tol) and (abs(c[9]) < tol) and (abs(c[10]) < tol) and
(abs(c[12]) < tol) and (abs(c[13]) < tol) and (abs(c[14]) < tol) and
(abs(c[16]) < tol) and (abs(c[17]) < tol) and (abs(c[19]) < tol) and
(abs(c[0]-c[6]) < tol) and (abs(c[6]-c[11]) < tol) and
(abs(c[15]-c[18]) < tol) and (abs(c[18]-c[20]) < tol) and
(abs(c[1]-c[2]) < tol) and (abs(c[2]-c[7]) < tol) and
(abs(c[0]-c[1]-2*c[15]) < tol))
def c21_cvoigt(c):
"""
Convert an input stiffness tensor to voigt form
The stiffness tensor needs to be provided as a list of the 21 elements of the
upper triangular part of the Voigt stiffness matrix::
C(21) = ( C_11, C_12, C_13, C_14, C_15, C_16,
C_22, C_23, C_24, C_25, C_26,
C_33, C_34, C_35, C_36,
C_44, C_45, C_46,
C_55, C_56,
C_66 )
(upper triangular part of Voigt stiffness matrix)
Args
c: input list of stiffness parameters (21 required)
Returns
numpy.array: voigt form of the stiffness tensor
"""
if cmoment_tensor_conversion:
try:
return np.asarray(cmoment_tensor_conversion.c21_cvoigt(c))
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
return np.array([[c[0], c[1], c[2], np.sqrt(2)*c[3], np.sqrt(2)*c[4], np.sqrt(2)*c[5]],
[c[1], c[6], c[7], np.sqrt(
2)*c[8], np.sqrt(2)*c[9], np.sqrt(2)*c[10]],
[c[2], c[7], c[11], np.sqrt(
2)*c[12], np.sqrt(2)*c[13], np.sqrt(2)*c[14]],
[np.sqrt(2)*c[3], np.sqrt(2)*c[8], np.sqrt(2)*c[12],
2*c[15], 2*c[16], 2*c[17]],
[np.sqrt(2)*c[4], np.sqrt(2)*c[9], np.sqrt(2)*c[13],
2*c[16], 2*c[18], 2*c[19]],
[np.sqrt(2)*c[5], np.sqrt(2)*c[10], np.sqrt(2)*c[14],
2*c[17], 2*c[19], 2*c[20]]])
def c_norm(c):
"""
Calculate norm of the stiffness tensor.
The stiffness tensor needs to be provided as a list of the 21 elements of the
upper triangular part of the Voigt stiffness matrix::
C(21) = ( C_11, C_12, C_13, C_14, C_15, C_16,
C_22, C_23, C_24, C_25, C_26,
C_33, C_34, C_35, C_36,
C_44, C_45, C_46,
C_55, C_56,
C_66 )
(upper triangular part of Voigt stiffness matrix)
Args
c: input list of stiffness parameters (21 required)
Returns
float: Euclidean norm of the tensor
"""
if cmoment_tensor_conversion:
try:
return cmoment_tensor_conversion.c_norm(c)
except Exception:
logger.exception('Error with C extension')
else:
logger.info(C_EXTENSION_FALLBACK_LOG_MSG)
return np.sqrt(c[0]**2 + c[6]**2 + c[11]**2 + 2*(c[1]**2+c[2]**2+c[7]**2) +
4*(c[3]**2+c[4]**2+c[5]**2+c[8]**2+c[9]**2+c[10]**2 +
c[12]**2+c[13]**2+c[14]**2+c[15]**2+c[18]**2+c[20]**2) +
8*(c[16]**2+c[17]**2+c[19]**2))
There are also Cython functions:
#!python
# cython: infer_types=True
# **Restricted: For Non-Commercial Use Only**
# This code is protected intellectual property and is available solely for teaching
# and non-commercially funded academic research purposes.
#
# Applications for commercial use should be made to Schlumberger or the University of Cambridge.
from libc.math cimport fabs
from libc.math cimport fmod
from libc.math cimport cos
from libc.math cimport acos
from libc.math cimport atan2
from libc.math cimport sin
from libc.math cimport sqrt
from libc.math cimport M_PI as pi
from libc.math cimport M_SQRT2 as sqrt2
from libc.stdlib cimport free
import unittest
cimport cython
cimport numpy as np
from scipy.linalg import eigh
from scipy.optimize import fsolve
from scipy import __version__ as __scipy_version__
import numpy as np
from cpython cimport bool
from ..utilities.unittest_utils import TestCase
cdef DTYPE_t PI2=2*pi
cdef DTYPE_t sqrt3=sqrt(3)
cdef DTYPE_t rad_cor=180/pi
cdef int check_finite=1
if int(__scipy_version__.split('.')[0])==0 and int(__scipy_version__.split('.')[1])<13:
check_finite=0
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef inline DTYPE_t[::1] cE_tk(DTYPE_t[::1] E,DTYPE_t[::1] results) nogil:
cdef DTYPE_t iso=(E[0]+E[1]+E[2])/3
cdef DTYPE_t dev0=E[0]-iso
cdef DTYPE_t dev1=E[2]-iso#Odd sorting from hudson paper E[0]>=E[2]>=E[1]
cdef DTYPE_t dev2=E[1]-iso
if dev2>0:
results[5]=iso/(fabs(iso)-dev1)
results[6]=-2*(dev2)/(dev1)
elif dev2<0:
results[5]=iso/(fabs(iso)+dev0)
results[6]=2*dev2/dev0
else:
results[5]=iso/(fabs(iso)+dev0)
results[6]=0
results[6]=results[6]*(1-fabs(results[5]));
return results
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef inline DTYPE_t[::1] ctk_uv(DTYPE_t[::1] results) nogil:
cdef DTYPE_t k=results[5]
cdef DTYPE_t tau=results[6]
if tau>0 and k>0:
if tau<4*k:
results[5]=tau/(1-(tau/2))
results[6]=k/(1-(tau/2))
else:
results[5]=tau/(1-2*k)
results[6]=k/(1-2*k)
elif tau<0 and k<0:
if tau>4*k:
results[5]=tau/(1+(tau/2))
results[6]=k/(1+(tau/2))
else:
results[5]=tau/(1+2*k)
results[6]=k/(1+2*k)
else:
results[5]=tau
results[6]=k
return results
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef inline void cE_gd(DTYPE_t[:] E,DTYPE_t*g,DTYPE_t*d) nogil:
if E[0]==E[2]:
if E[0]>0:
g[0]=0
d[0]=pi/2
elif E[0]<0:
g[0]=0
d[0]=-pi/2
else:
g[0]=atan2(-E[0]+2*E[1]-E[2],sqrt3*(E[0]-E[2]))
d[0]=pi/2-acos((E[0]+E[1]+E[2])/(sqrt3*sqrt((E[0]*E[0]+E[1]*E[1]+E[2]*E[2]))))
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef inline void cN_SDR(DTYPE_t N0,DTYPE_t N1,DTYPE_t N2,DTYPE_t S0,DTYPE_t S1,DTYPE_t S2,DTYPE_t* strike,DTYPE_t*dip,DTYPE_t*rake) nogil:
if N2 > 0:
N0 = -N0
N1 = -N1
N2 = -N2
S0 = -S0
S1 = -S1
S2 = -S2
strike[0] = atan2(-N0, N1)
dip[0] = atan2((N1**2+N0**2), sqrt((N0*N2)**2+(N1*N2)**2))
rake[0] = atan2(-S2, S0*N1 - S1*N0)
if dip[0] > pi/2:
dip[0] = pi-dip[0]
strike[0] = strike[0]+pi
rake[0] = PI2-rake[0]
if fabs(strike[0]) > PI2:
strike[0] = fmod(strike[0],PI2)
if strike[0] < 0.0:
strike[0] += PI2
if rake[0] > pi:
rake[0] -= PI2
elif rake[0] < -pi:
rake[0] += PI2
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef DTYPE_t[::1] cTP_SDR(DTYPE_t[::1] T,DTYPE_t[::1] P,DTYPE_t[::1] results) nogil:
cdef DTYPE_t Nt=sqrt((T[0]+P[0])*(T[0]+P[0])+(T[1]+P[1])*(T[1]+P[1])+(T[2]+P[2])*(T[2]+P[2]))
cdef DTYPE_t St=sqrt((T[0]-P[0])*(T[0]-P[0])+(T[1]-P[1])*(T[1]-P[1])+(T[2]-P[2])*(T[2]-P[2]))
cdef DTYPE_t N2=(T[2]+P[2])/Nt
cdef DTYPE_t N0=(T[0]+P[0])/Nt
cdef DTYPE_t N1=(T[1]+P[1])/Nt
cdef DTYPE_t S2=(T[2]-P[2])/St
cdef DTYPE_t S0=(T[0]-P[0])/St
cdef DTYPE_t S1=(T[1]-P[1])/St
cdef DTYPE_t strike1=0.
cdef DTYPE_t dip1=0.
cdef DTYPE_t rake1=0.
cdef DTYPE_t strike2=0.
cdef DTYPE_t dip2=0.
cdef DTYPE_t rake2=0.
cN_SDR(N0,N1,N2,S0,S1,S2,&strike1,&dip1,&rake1)
cN_SDR(S0,S1,S2,N0,N1,N2,&strike2,&dip2,&rake2)
#Switch for sdr2
if fabs(rake1)>pi/2:
results[2]=strike2
results[3]=cos(dip2)
results[4]=rake2
else:
results[2]=strike1
results[3]=cos(dip1)
results[4]=rake1
results[7]=strike1*rad_cor
results[8]=dip1*rad_cor
results[9]=rake1*rad_cor
results[10]=strike2*rad_cor
results[11]=dip2*rad_cor
results[12]=rake2*rad_cor
return results
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef DTYPE_t[::1] cTPE_convert(DTYPE_t[::1] T,DTYPE_t[::1] P,DTYPE_t[::1]E,DTYPE_t[::1]results) nogil:
cE_gd(E,&results[0],&results[1])
results=cE_tk(E,results)
results=ctk_uv(results)
results=cTP_SDR(T,P,results)
return results
#Tape to Moment tensor conversion
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef void cTape_MT6(DTYPE_t*M, DTYPE_t gamma,DTYPE_t delta,DTYPE_t kappa,DTYPE_t h,DTYPE_t sigma) nogil:
#Get T N P axes
cdef DTYPE_t ck=cos(kappa)
cdef DTYPE_t cs=cos(sigma)
cdef DTYPE_t sk=sin(kappa)
cdef DTYPE_t ss=sin(sigma)
cdef DTYPE_t cg=cos(gamma)
cdef DTYPE_t sg=sin(gamma)
cdef DTYPE_t cd=cos(delta)
cdef DTYPE_t sh=sqrt(1-h*h)
cdef DTYPE_t NT=sqrt((ck*cs+sk*h*ss-sk*sh)*(ck*cs+sk*h*ss-sk*sh)+(sk*cs-ck*h*ss+ck*sh)*(sk*cs-ck*h*ss+ck*sh)+(-sh*ss-h)*(-sh*ss-h))
cdef DTYPE_t NP=sqrt((ck*cs+sk*h*ss+sk*sh)*(ck*cs+sk*h*ss+sk*sh)+(sk*cs-ck*h*ss-ck*sh)*(sk*cs-ck*h*ss-ck*sh)+(-sh*ss+h)*(-sh*ss+h))
cdef DTYPE_t T1=(ck*cs+sk*h*ss-sk*sh)/NT
cdef DTYPE_t T2=(sk*cs-ck*h*ss+ck*sh)/NT
cdef DTYPE_t T3=(-sh*ss-h)/NT
cdef DTYPE_t P1=(ck*cs+sk*h*ss+sk*sh)/NP
cdef DTYPE_t P2=(sk*cs-ck*h*ss-ck*sh)/NP
cdef DTYPE_t P3=(-sh*ss+h)/NP
cdef DTYPE_t N1=T2*P3-P2*T3
cdef DTYPE_t N2=-T1*P3+P1*T3
cdef DTYPE_t N3=T1*P2-T2*P1
#Normalised TNP Values
cdef DTYPE_t sd=sin(delta)
cdef DTYPE_t L1=(sqrt(3)*cg*cd-sg*cd+sqrt(2)*sd)/sqrt(6)
cdef DTYPE_t L2=(2*sg*cd+sqrt(2)*sd)/sqrt(6)
cdef DTYPE_t L3=(-sqrt(3)*cg*cd-sg*cd+sqrt(2)*sd)/sqrt(6)
M[0]=L1*T1*T1+L2*N1*N1+L3*P1*P1
M[1]=L1*T2*T2+L2*N2*N2+L3*P2*P2
M[2]=L1*T3*T3+L2*N3*N3+L3*P3*P3
M[3]=sqrt2*(L1*T1*T2+L2*N1*N2+L3*P1*P2)
M[4]=sqrt2*(L1*T1*T3+L2*N1*N3+L3*P1*P3)
M[5]=sqrt2*(L1*T2*T3+L2*N2*N3+L3*P2*P3)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef void cMultipleTape_MT6(DTYPE_t*M,DTYPE_t* gamma,DTYPE_t* delta,DTYPE_t* kappa,DTYPE_t* h,DTYPE_t* sigma,Py_ssize_t n) nogil:
for i from 0<=i<n:
cTape_MT6(&M[i*6],gamma[i],delta[i],kappa[i],h[i],sigma[i])
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
def Tape_MT6(DTYPE_t[::1] gamma,DTYPE_t[::1] delta,DTYPE_t[::1] kappa,DTYPE_t[::1] h,DTYPE_t[::1] sigma):
cdef Py_ssize_t n=gamma.shape[0]
if not (gamma.shape[0]==delta.shape[0] and gamma.shape[0]==kappa.shape[0] and gamma.shape[0]==h.shape[0] and gamma.shape[0]==sigma.shape[0]):
raise ValueError('Arguments different size')
cdef DTYPE_t[:,::1] MTs=np.empty((n,6))
cMultipleTape_MT6(&MTs[0,0],&gamma[0],&delta[0],&kappa[0],&h[0],&sigma[0],n)
return np.ascontiguousarray(np.asarray(MTs).T)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cpdef MT_convert(DTYPE_t[:] MT,DTYPE_t [::1] results):
T,N,P,E=cMT6_TNPE(MT)
return cTPE_convert(T,P,E,results)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
def MT_output_convert(DTYPE_t[:,:] MTs):
cdef Py_ssize_t imax=MTs.shape[1]
cdef Py_ssize_t i
cdef DTYPE_t [::1] g=np.empty((imax,))
cdef DTYPE_t [::1] d=np.empty((imax,))
cdef DTYPE_t [::1] k=np.empty((imax,))
cdef DTYPE_t [::1] h=np.empty((imax,))
cdef DTYPE_t [::1] s=np.empty((imax,))
cdef DTYPE_t [::1] u=np.empty((imax,))
cdef DTYPE_t [::1] v=np.empty((imax,))
cdef DTYPE_t [::1] s1=np.empty((imax,))
cdef DTYPE_t [::1] d1=np.empty((imax,))
cdef DTYPE_t [::1] r1=np.empty((imax,))
cdef DTYPE_t [::1] s2=np.empty((imax,))
cdef DTYPE_t [::1] d2=np.empty((imax,))
cdef DTYPE_t [::1] r2=np.empty((imax,))
cdef DTYPE_t [::1] results=np.empty((13,))
for i in range(imax):
results=MT_convert(MTs[:,i],np.empty((13,)))
g[i]=results[0]
d[i]=results[1]
k[i]=results[2]
h[i]=results[3]
s[i]=results[4]
u[i]=results[5]
v[i]=results[6]
s1[i]=results[7]
d1[i]=results[8]
r1[i]=results[9]
s2[i]=results[10]
d2[i]=results[11]
r2[i]=results[12]
return {'g':np.asarray(g),'d':np.asarray(d),'k':np.asarray(k),'h':np.asarray(h),'s':np.asarray(s),'u':np.asarray(u),'v':np.asarray(v),'S1':np.asarray(s1),'D1':np.asarray(d1),'R1':np.asarray(r1),'S2':np.asarray(s2),'D2':np.asarray(d2),'R2':np.asarray(r2)}
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef cMT6_TNPE(DTYPE_t[:] MT):
cdef DTYPE_t [:,::1] MT33=np.empty((3,3))
cdef DTYPE_t [::1] E=np.empty((3,))
cdef DTYPE_t [::1] e=np.empty((3,))
cdef DTYPE_t [::1] T=np.empty((3,))
cdef DTYPE_t [::1] P=np.empty((3,))
cdef DTYPE_t [::1] N=np.empty((3,))
cdef DTYPE_t [:,:] L=np.empty((3,3))
MT33[0,0]=MT[0]
MT33[1,1]=MT[1]
MT33[2,2]=MT[2]
MT33[0,1]=MT[3]/sqrt2
MT33[0,2]=MT[4]/sqrt2
MT33[1,2]=MT[5]/sqrt2
MT33[1,0]=MT[3]/sqrt2
MT33[2,0]=MT[4]/sqrt2
MT33[2,1]=MT[5]/sqrt2
if check_finite>0:
e,L=eigh(MT33,overwrite_a=True,check_finite=False)
else:
e,L=eigh(MT33,overwrite_a=True)
maxi=0
mini=0
cdef Py_ssize_t i
for i in [0,1,2]:
if e[i]>e[maxi]:
maxi=i
if e[i]<e[mini]:
mini=i
E[0]=e[maxi]
E[1]=e[3-maxi-mini]
E[2]=e[mini]
T[0]=L[0,maxi]
T[1]=L[1,maxi]
T[2]=L[2,maxi]
N[0]=L[0,3-maxi-mini]
N[1]=L[1,3-maxi-mini]
N[2]=L[2,3-maxi-mini]
P[0]=L[0,mini]
P[1]=L[1,mini]
P[2]=L[2,mini]
return T,N,P,E
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cpdef MT6_TNPE(DTYPE_t[:,:] MTs):
cdef Py_ssize_t imax=MTs.shape[1]
cdef Py_ssize_t i
cdef DTYPE_t [:,::1] T=np.empty((3,imax))
cdef DTYPE_t [:,::1] N=np.empty((3,imax))
cdef DTYPE_t [:,::1] P=np.empty((3,imax))
cdef DTYPE_t [:,::1] E=np.empty((3,imax))
for i in range(imax):
t,n,p,e=cMT6_TNPE(MTs[:,i])
T[0,i]=t[0]
N[0,i]=n[0]
P[0,i]=p[0]
E[0,i]=e[0]
T[1,i]=t[1]
N[1,i]=n[1]
P[1,i]=p[1]
E[1,i]=e[1]
T[2,i]=t[2]
N[2,i]=n[2]
P[2,i]=p[2]
E[2,i]=e[2]
return np.asarray(T),np.asarray(N),np.asarray(P),np.asarray(E)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cpdef E_GD(DTYPE_t[:,:]E):
cdef Py_ssize_t imax=E.shape[1]
cdef Py_ssize_t i
cdef DTYPE_t [::1] g=np.empty((imax,))
cdef DTYPE_t [::1] d=np.empty((imax,))
for i in range(imax):
cE_gd(E[:,i],&g[i],&d[i])
return np.asarray(g),np.asarray(d)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef cSingleTP_SDR(DTYPE_t [:] T,DTYPE_t[:]P):
cdef DTYPE_t t0 = T[0]
cdef DTYPE_t t1 = T[1]
cdef DTYPE_t t2 = T[2]
cdef DTYPE_t p0 = P[0]
cdef DTYPE_t p1 = P[1]
cdef DTYPE_t p2 = P[2]
cdef DTYPE_t Nt = sqrt((t0+p0)*(t0+p0)+(t1+p1)*(t1+p1)+(t2+p2)*(t2+p2))
cdef DTYPE_t St = sqrt((t0-p0)*(t0-p0)+(t1-p1)*(t1-p1)+(t2-p2)*(t2-p2))
cdef DTYPE_t N0 = (t0+p0)/Nt
cdef DTYPE_t N1 = (t1+p1)/Nt
cdef DTYPE_t N2 = (t2+p2)/Nt
cdef DTYPE_t S0 = (t0-p0)/St
cdef DTYPE_t S1 = (t1-p1)/St
cdef DTYPE_t S2 = (t2-p2)/St
cdef DTYPE_t strike1
cdef DTYPE_t dip1
cdef DTYPE_t rake1
cN_SDR(N0, N1, N2, S0, S1, S2, &strike1, &dip1, &rake1)
return strike1, dip1, rake1
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cpdef TP_SDR(DTYPE_t[:,:]T,DTYPE_t[:,:]P):
cdef Py_ssize_t imax=T.shape[1]
cdef Py_ssize_t i
cdef DTYPE_t [::1] s=np.empty((imax,))
cdef DTYPE_t [::1] d=np.empty((imax,))
cdef DTYPE_t [::1] r=np.empty((imax,))
for i in range(imax):
s[i],d[i],r[i]=cSingleTP_SDR(T[:,i],P[:,i])
return np.asarray(s),np.asarray(d),np.asarray(r)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef void csingleSDR_SDR(DTYPE_t s,DTYPE_t d,DTYPE_t r,DTYPE_t*s2,DTYPE_t*d2,DTYPE_t*r2):
cdef DTYPE_t ck=cos(s)
cdef DTYPE_t cs=cos(r)
cdef DTYPE_t sk=sin(s)
cdef DTYPE_t ss=sin(r)
cdef DTYPE_t h=cos(d)
cdef DTYPE_t sh=sqrt(1-h*h)
cdef DTYPE_t NT=sqrt((ck*cs+sk*h*ss-sk*sh)*(ck*cs+sk*h*ss-sk*sh)+(sk*cs-ck*h*ss+ck*sh)*(sk*cs-ck*h*ss+ck*sh)+(-sh*ss-h)*(-sh*ss-h))
cdef DTYPE_t NP=sqrt((ck*cs+sk*h*ss+sk*sh)*(ck*cs+sk*h*ss+sk*sh)+(sk*cs-ck*h*ss-ck*sh)*(sk*cs-ck*h*ss-ck*sh)+(-sh*ss+h)*(-sh*ss+h))
cdef DTYPE_t T1=(ck*cs+sk*h*ss-sk*sh)/NT
cdef DTYPE_t T2=(sk*cs-ck*h*ss+ck*sh)/NT
cdef DTYPE_t T3=(-sh*ss-h)/NT
cdef DTYPE_t P1=(ck*cs+sk*h*ss+sk*sh)/NP
cdef DTYPE_t P2=(sk*cs-ck*h*ss-ck*sh)/NP
cdef DTYPE_t P3=(-sh*ss+h)/NP
cdef DTYPE_t N0=(T1-P1)
cdef DTYPE_t N1=(T2-P2)
cdef DTYPE_t N2=(T3-P3)
cdef DTYPE_t S0=(T1+P1)
cdef DTYPE_t S1=(T2+P2)
cdef DTYPE_t S2=(T3+P3)
cdef DTYPE_t St=sqrt(S0*S0+S1*S1+S2*S2)
cdef DTYPE_t Nt=sqrt(N0*N0+N1*N1+N2*N2)
cN_SDR(S0/St,S1/St,S2/St,N0/Nt,N1/Nt,N2/Nt,s2,d2,r2)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cpdef SDR_SDR(DTYPE_t[::1] s,DTYPE_t[::1] d,DTYPE_t[::1] r):
cdef Py_ssize_t imax=s.shape[0]
cdef Py_ssize_t i
cdef DTYPE_t [::1] s2=np.empty((imax,))
cdef DTYPE_t [::1] d2=np.empty((imax,))
cdef DTYPE_t [::1] r2=np.empty((imax,))
for i in range(imax):
csingleSDR_SDR(s[i],d[i],r[i],&s2[i],&d2[i],&r2[i])
return np.asarray(s2),np.asarray(d2),np.asarray(r2)
#Bi-axes
cpdef list isotropic_c(DTYPE_t l=1.,DTYPE_t mu=1.,list c=[]):
if len(c)==21: #Calculate isotropic approacximation
#Eqns 81a and 81b from chapman and leaney 2011
mu=((c[0]+c[6]+c[11])+3*(c[15]+c[18]+c[20])-(c[1]+c[2]+c[7]))/15
l=((c[0]+c[6]+c[11])-2*(c[15]+c[18]+c[20])+4*(c[1]+c[2]+c[7]))/15
n=l+2*mu
return [n,l,l,0,0,0,n,l,0,0,0,n,0,0,0,mu,0,0,mu,0,mu]
def MT6_biaxes(DTYPE_t[:]MT6,list c):
lambda2mu=(3*(c[0]+c[6]+c[11])+4*(c[15]+c[18]+c[20])+2*(c[1]+c[2]+c[7]))/15
mu=((c[0]+c[6]+c[11])+3*(c[15]+c[18]+c[20])-(c[1]+c[2]+c[7]))/15
l=lambda2mu-2*mu
T,N,P,E=cMT6_TNPE(MT6)
isotropic=(l+mu)*E[1]/mu-l*(E[0]+E[2])/(2*mu)
if is_isotropic_c(c):
explosion=isotropic
else:
def isotropic_solve(iso):
iso6=np.squeeze(iso)*np.array([[1],[1],[1],[0],[0],[0]])
if iso6.shape!=MT6.shape:
iso6=iso6.T
T,N,P,E=cMT6_TNPE(MT6c_D6(np.squeeze(MT6-iso6),c).flatten())
return E[1]
explosion=fsolve(isotropic_solve,isotropic)
explosion6=np.squeeze(explosion)*np.array([[1],[1],[1],[0],[0],[0]])
if explosion6.shape!=MT6.shape:
explosion6=explosion6.T
T,N,P,E=cMT6_TNPE(MT6c_D6(np.squeeze(MT6-explosion6),c).flatten())
area_displacement = E[0]-E[2]
phi=np.zeros((3,2))
if area_displacement!=0: # to avoid undefined
cphi=np.squeeze(np.sqrt(E[0]/area_displacement))
sphi=np.squeeze(np.sqrt(-E[2]/area_displacement))
phi[:,0]=np.array(cphi*T+sphi*P).flatten()
phi[:,1]=np.array(cphi*T-sphi*P).flatten()
return phi,explosion,area_displacement
cpdef MT6c_D6(mt6,list c=isotropic_c(l=1,mu=1)):
mtvoigt=mt6[np.array([0,1,2,5,4,3])]
mtvoigt=np.matrix(mtvoigt)
if mtvoigt.shape[1]==6:
mtvoigt=mtvoigt.T
#Convert to voigt
dvoigt=np.linalg.solve(np.matrix(c21_cvoigt(c)),mtvoigt)
return np.asarray(dvoigt[np.array([0,1,2,5,4,3])])
cpdef bool is_isotropic_c(list c):
cdef DTYPE_t tol = 1.e-6*c_norm(c);
return ((fabs(c[ 3])<tol)and(fabs(c[ 4])<tol)and(fabs(c[ 5])<tol)and \
(fabs(c[ 8])<tol)and(fabs(c[ 9])<tol)and(fabs(c[10])<tol)and \
(fabs(c[12])<tol)and(fabs(c[13])<tol)and(fabs(c[14])<tol)and \
(fabs(c[16])<tol)and(fabs(c[17])<tol)and(fabs(c[19])<tol)and \
(fabs(c[ 0]-c[ 6])<tol)and(fabs(c[ 6]-c[11])<tol)and \
(fabs(c[15]-c[18])<tol)and(fabs(c[18]-c[20])<tol)and \
(fabs(c[ 1]-c[ 2])<tol)and(fabs(c[ 2]-c[ 7])<tol)and \
(fabs(c[0]-c[1]-2*c[15])<tol))
cpdef DTYPE_t [:,::1] c21_cvoigt(list c):
return np.array([[c[0],c[1],c[2],sqrt2*c[3],sqrt2*c[4],sqrt2*c[5]],
[c[1],c[6],c[7],sqrt2*c[8],sqrt2*c[9],sqrt2*c[10]],
[c[2],c[7],c[11],sqrt2*c[12],sqrt2*c[13],sqrt2*c[14]],
[sqrt2*c[3],sqrt2*c[8],sqrt2*c[12],2*c[15],2*c[16],2*c[17]],
[sqrt2*c[4],sqrt2*c[9],sqrt2*c[13],2*c[16],2*c[18],2*c[19]],
[sqrt2*c[5],sqrt2*c[10],sqrt2*c[14],2*c[17],2*c[19],2*c[20]]])
cpdef DTYPE_t c_norm(list c):
return sqrt(c[0]**2+c[6]**2+c[11]**2+2*(c[1]**2+c[2]**2+c[7]**2)+
4*(c[3]**2+c[4]**2+c[5]**2+c[8]**2+c[9]**2+c[10]**2+
c[12]**2+c[13]**2+c[14]**2+c[15]**2+c[18]**2+c[20]**2)+
8*(c[16]**2+c[17]**2+c[19]**2));