# http://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg----------
This module provides various moment tensor related utility functions.
It can be used to convert between strike-dip-rake and moment tensor representations and provides different options to produce random moment tensors.
Moment tensors are represented by :py:class:`MomentTensor` instances. The internal representation uses a north-east-down (NED) coordinate system, but it can convert from/to the conventions used by the Global CMT catalog (up-south-east, USE).
If not otherwise noted, scalar moment is interpreted as the Frobenius norm based scalar moment (see :py:meth:`MomentTensor.scalar_moment`. The scalar moment according to the "standard decomposition" can be found in the output of :py:meth:`MomentTensor.standard_decomposition`. '''
''' Get randomly oriented unit vector.
:param rstate: :py:class:`numpy.random.RandomState` object, can be used to create reproducible pseudo-random sequences ''' rstate = rstate or num.random while True: axis = rstate.uniform(size=3) * 2.0 - 1.0 uabs = math.sqrt(num.sum(axis**2)) if 0.001 < uabs < 1.0: return axis / uabs
''' Build rotation matrix based on axis and angle.
:param angle: rotation angle [degrees] :param axis: orientation of rotation axis, either in spherical coordinates ``(theta, phi)`` [degrees], or as a unit vector ``(ux, uy, uz)``. '''
if len(axis) == 2: theta, phi = axis ux = math.sin(d2r*theta)*math.cos(d2r*phi) uy = math.sin(d2r*theta)*math.sin(d2r*phi) uz = math.cos(d2r*theta)
elif len(axis) == 3: axis = num.asarray(axis) uabs = math.sqrt(num.sum(axis**2)) ux, uy, uz = axis / uabs else: assert False
ct = math.cos(d2r*angle) st = math.sin(d2r*angle) return num.matrix([ [ct + ux**2*(1.-ct), ux*uy*(1.-ct)-uz*st, ux*uz*(1.-ct)+uy*st], [uy*ux*(1.-ct)+uz*st, ct+uy**2*(1.-ct), uy*uz*(1.-ct)-ux*st], [uz*ux*(1.-ct)-uy*st, uz*uy*(1.-ct)+ux*st, ct+uz**2*(1.-ct)] ])
'''Get random rotation matrix.
A random rotation matrix, drawn from a uniform distrubution in the space of rotations is returned, after Avro 1992 - "Fast random rotation matrices".
:param x: three (uniform random) numbers in the range [0, 1[ used as input to the distribution tranformation. If ``None``, random numbers are used. Can be used to create grids of random rotations with uniform density in rotation space. '''
if x is not None: x1, x2, x3 = x else: x1, x2, x3 = num.random.random(3)
phi = math.pi*2.0*x1
zrot = num.matrix([ [math.cos(phi), math.sin(phi), 0.], [-math.sin(phi), math.cos(phi), 0.], [0., 0., 1.]])
lam = math.pi*2.0*x2
v = num.matrix([[ math.cos(lam)*math.sqrt(x3), math.sin(lam)*math.sqrt(x3), math.sqrt(1.-x3)]]).T
house = num.identity(3) - 2.0 * v * v.T return -house*zrot
return float(num.random.uniform(mi, ma))
mi_ = 0.5*(math.cos(mi * math.pi/180.)+1.) ma_ = 0.5*(math.cos(ma * math.pi/180.)+1.) return math.acos(rand(mi_, ma_)*2.-1.)*180./math.pi
strikemin=0., strikemax=360., dipmin=0., dipmax=90., rakemin=-180., rakemax=180.):
''' Get random strike, dip, rake triplet.
.. note::
Might not produce a homogeneous distribution of mechanisms. Better use :py:meth:`MomentTensor.random_dc` which is based on :py:func:`random_rotation`. '''
strike = rand(strikemin, strikemax) dip = randdip(dipmin, dipmax) rake = rand(rakemin, rakemax)
return strike, dip, rake
'''Get non-redundant components from symmetric 3x3 matrix
:returns: 1D NumPy array with entries ordered like ``(a_xx, a_yy, a_zz, a_xy, a_xz, a_yz)`` '''
return num.array([m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2]])
''' Create symmetric 3x3 matrix from its 6 non-redundant values. '''
return num.matrix([[a_xx, a_xy, a_xz], [a_xy, a_yy, a_yz], [a_xz, a_yz, a_zz]], dtype=num.float)
'''Convert anything to moment tensor represented as a NumPy matrix.
Transforms :py:class:`MomentTensor` objects, tuples, lists and NumPy arrays with 3x3 or 3, 4, 6, or 7 elements into NumPy 3x3 matrix objects.
The ``values`` argument is interpreted depending on shape and type as follows:
* ``(strike, dip, rake)`` * ``(strike, dip, rake, magnitude)`` * ``(mnn, mee, mdd, mne, mnd, med)`` * ``(mnn, mee, mdd, mne, mnd, med, magnitude)`` * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))`` * :py:class:`MomentTensor` '''
if isinstance(values, (tuple, list)): values = num.asarray(values, dtype=num.float)
if isinstance(values, MomentTensor): return values.m()
elif isinstance(values, num.ndarray): if values.shape == (3,): strike, dip, rake = values rotmat1 = euler_to_matrix(d2r*dip, d2r*strike, -d2r*rake) return rotmat1.T * MomentTensor._m_unrot * rotmat1
elif values.shape == (4,): strike, dip, rake, magnitude = values moment = magnitude_to_moment(magnitude) rotmat1 = euler_to_matrix(d2r*dip, d2r*strike, -d2r*rake) return rotmat1.T * MomentTensor._m_unrot * rotmat1 * moment
elif values.shape == (6,): return symmat6(*values)
elif values.shape == (7,): magnitude = values[6] moment = magnitude_to_moment(magnitude) mt = symmat6(*values[:6]) mt *= moment / (num.linalg.norm(mt) / math.sqrt(2.0)) return mt
elif values.shape == (3, 3): return num.asmatrix(values, dtype=num.float)
raise Exception('cannot convert object to 3x3 matrix')
''' Convert scalar moment to moment magnitude Mw.
:param moment: scalar moment [Nm] :returns: moment magnitude Mw
Moment magnitude is defined as
.. math::
M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0) - 10.7
where :math:`M_0` is the scalar moment given in [Nm].
.. note::
Global CMT uses 10.7333333 instead of 10.7, based on [Kanamori 1977], 10.7 is from [Hanks and Kanamori 1979]. '''
''' Convert moment magnitude Mw to scalar moment.
:param magnitude: moment magnitude :returns: scalar moment [Nm]
See :py:func:`moment_to_magnitude`. '''
return 10.0**(1.5*(magnitude+10.7))*1.0e-7
'''Given euler angle triplet, create rotation matrix
Given coordinate system `(x,y,z)` and rotated system `(xs,ys,zs)` the line of nodes is the intersection between the `x,y` and the `xs,ys` planes.
:param alpha: is the angle between the `z`-axis and the `zs`-axis [rad] :param beta: is the angle between the `x`-axis and the line of nodes [rad] :param gamma: is the angle between the line of nodes and the `xs`-axis [rad]
Usage for moment tensors::
m_unrot = numpy.matrix([[0,0,-1],[0,0,0],[-1,0,0]]) euler_to_matrix(dip,strike,-rake, rotmat) m = rotmat.T * m_unrot * rotmat
'''
ca = math.cos(alpha) cb = math.cos(beta) cg = math.cos(gamma) sa = math.sin(alpha) sb = math.sin(beta) sg = math.sin(gamma)
mat = num.matrix([[cb*cg-ca*sb*sg, sb*cg+ca*cb*sg, sa*sg], [-cb*sg-ca*sb*cg, -sb*sg+ca*cb*cg, sa*cg], [sa*sb, -sa*cb, ca]], dtype=num.float) return mat
'''Get eulerian angle triplet from rotation matrix.'''
ex = cvec(1., 0., 0.) ez = cvec(0., 0., 1.) exs = rotmat.T * ex ezs = rotmat.T * ez enodes = num.cross(ez.T, ezs.T).T if num.linalg.norm(enodes) < 1e-10: enodes = exs enodess = rotmat*enodes cos_alpha = float((ez.T*ezs)) if cos_alpha > 1.: cos_alpha = 1.
if cos_alpha < -1.: cos_alpha = -1.
alpha = math.acos(cos_alpha) beta = num.mod(math.atan2(enodes[1, 0], enodes[0, 0]), math.pi*2.) gamma = num.mod(-math.atan2(enodess[1, 0], enodess[0, 0]), math.pi*2.)
return unique_euler(alpha, beta, gamma)
'''Uniquify eulerian angle triplet.
Put eulerian angle triplet into ranges compatible with ``(dip, strike, -rake)`` conventions in seismology::
alpha (dip) : [0, pi/2] beta (strike) : [0, 2*pi) gamma (-rake) : [-pi, pi)
If ``alpha1`` is near to zero, ``beta`` is replaced by ``beta+gamma`` and ``gamma`` is set to zero, to prevent this additional ambiguity.
If ``alpha`` is near to ``pi/2``, ``beta`` is put into the range ``[0,pi)``. '''
pi = math.pi
alpha = num.mod(alpha, 2.0*pi)
if 0.5*pi < alpha and alpha <= pi: alpha = pi - alpha beta = beta + pi gamma = 2.0*pi - gamma elif pi < alpha and alpha <= 1.5*pi: alpha = alpha - pi gamma = pi - gamma elif 1.5*pi < alpha and alpha <= 2.0*pi: alpha = 2.0*pi - alpha beta = beta + pi gamma = pi + gamma
alpha = num.mod(alpha, 2.0*pi) beta = num.mod(beta, 2.0*pi) gamma = num.mod(gamma+pi, 2.0*pi)-pi
# If dip is exactly 90 degrees, one is still # free to choose between looking at the plane from either side. # Choose to look at such that beta is in the range [0,180)
# This should prevent some problems, when dip is close to 90 degrees: if abs(alpha - 0.5*pi) < 1e-10: alpha = 0.5*pi
if abs(beta - pi) < 1e-10: beta = pi
if abs(beta - 2.*pi) < 1e-10: beta = 0.
if abs(beta) < 1e-10: beta = 0.
if alpha == 0.5*pi and beta >= pi: gamma = - gamma beta = num.mod(beta-pi, 2.0*pi) gamma = num.mod(gamma+pi, 2.0*pi)-pi assert 0. <= beta < pi assert -pi <= gamma < pi
if alpha < 1e-7: beta = num.mod(beta + gamma, 2.0*pi) gamma = 0.
return (alpha, beta, gamma)
return num.matrix([[x, y, z]], dtype=num.float).T
return num.matrix([[x, y, z]], dtype=num.float)
if magnitude is not None: scalar_moment = magnitude_to_moment(magnitude)
if x is None: x = num.random.random(6)
evals = x[:3] * 2. - 1.0 evals /= num.sqrt(num.sum(evals**2)) / math.sqrt(2.0) rotmat = random_rotation(x[3:]) return scalar_moment * rotmat * num.matrix(num.diag(evals)) * rotmat.T
return to6(random_mt(*args, **kwargs))
if magnitude is not None: scalar_moment = magnitude_to_moment(magnitude)
rotmat = random_rotation(x) return scalar_moment * (rotmat * MomentTensor._m_unrot * rotmat.T)
return "/ %5.2F %5.2F %5.2F \\\n" % (m[0, 0], m[0, 1], m[0, 2]) + \ "| %5.2F %5.2F %5.2F |\n" % (m[1, 0], m[1, 1], m[1, 2]) + \ "\\ %5.2F %5.2F %5.2F /\n" % (m[2, 0], m[2, 1], m[2, 2])
''' Convenience function to convert various objects to moment tensor object.
Like :py:meth:``MomentTensor.from_values``, but does not create a new :py:class:`MomentTensor` object when ``mt`` already is one. '''
if isinstance(mt, MomentTensor): return mt else: return MomentTensor.from_values(mt)
''' Moment tensor object
:param m: NumPy matrix in north-east-down convention :param m_up_south_east: NumPy matrix in up-south-east convention :param strike,dip,rake: fault plane angles in [degrees] :param scalar_moment: scalar moment in [Nm] :param magnitude: moment magnitude Mw
Global CMT catalog moment tensors use the up-south-east (USE) coordinate system convention with :math:`r` (up), :math:`\\theta` (south), and :math:`\\phi` (east).
.. math:: :nowrap:
\\begin{align*} M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\ M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\ M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne} \\end{align*}
'''
[[0., 0., -1.], [0., -1., 0.], [-1., 0., 0.]], dtype=num.float) [[0., 0., -1.], [-1., 0., 0.], [0., 1., 0.]], dtype=num.float).T [[0., 0., -1.], [0., 0., 0.], [-1., 0., 0.]], dtype=num.float)
''' Create random oriented double-couple moment tensor
The rotations used are uniformly distributed in the space of rotations. ''' return MomentTensor( m=random_dc(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
''' Create random moment tensor
Moment tensors produced by this function appear uniformly distributed when shown in a Hudson's diagram. The rotations used are unifomly distributed in the space of rotations. ''' return MomentTensor( m=random_mt(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
def from_values(cls, values): ''' Alternative constructor for moment tensor objects
This constructor takes a :py:class:`MomentTensor` object, a tuple, list or NumPy array with 3x3 or 3, 4, 6, or 7 elements to build a Moment tensor object.
The ``values`` argument is interpreted depending on shape and type as follows:
* ``(strike, dip, rake)`` * ``(strike, dip, rake, magnitude)`` * ``(mnn, mee, mdd, mne, mnd, med)`` * ``(mnn, mee, mdd, mne, mnd, med, magnitude)`` * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))`` * :py:class:`MomentTensor` object '''
m = values_to_matrix(values) return MomentTensor(m=m)
self, m=None, m_up_south_east=None, strike=0., dip=0., rake=0., scalar_moment=1., mnn=None, mee=None, mdd=None, mne=None, mnd=None, med=None, strike1=None, dip1=None, rake1=None, strike2=None, dip2=None, rake2=None, magnitude=None, moment=None):
Object.__init__(self, init_props=False)
if any(mxx is not None for mxx in (mnn, mee, mdd, mne, mnd, med)): m = symmat6(mnn, mee, mdd, mne, mnd, med)
strike = d2r*strike dip = d2r*dip rake = d2r*rake
if m_up_south_east is not None: m = self._to_up_south_east * m_up_south_east * \ self._to_up_south_east.T
if m is None: if any(x is not None for x in ( strike1, dip1, rake1, strike2, dip2, rake2)): raise Exception( 'strike1, dip1, rake1, strike2, dip2, rake2 are read-only ' 'properties')
if moment is not None: scalar_moment = moment
if magnitude is not None: scalar_moment = magnitude_to_moment(magnitude)
rotmat1 = euler_to_matrix(dip, strike, -rake) m = rotmat1.T * MomentTensor._m_unrot * rotmat1 * scalar_moment
self._m = m self._update()
m_evals, m_evecs = eigh_check(self._m) if num.linalg.det(m_evecs) < 0.: m_evecs *= -1.
rotmat1 = (m_evecs * MomentTensor._u_evecs.T).T if num.linalg.det(rotmat1) < 0.: rotmat1 *= -1.
self._m_eigenvals = m_evals self._m_eigenvecs = m_evecs
self._rotmats = sorted( [rotmat1, MomentTensor._flip_dc * rotmat1], key=lambda m: num.abs(m.flat).tolist())
def mnn(self): return float(self._m[0, 0])
def mnn(self, value): self._m[0, 0] = value self._update()
def mee(self): return float(self._m[1, 1])
def mee(self, value): self._m[1, 1] = value self._update()
def mdd(self): return float(self._m[2, 2])
def mdd(self, value): self._m[2, 2] = value self._update()
def mne(self): return float(self._m[0, 1])
def mne(self, value): self._m[0, 1] = value self._m[1, 0] = value self._update()
def mnd(self): return float(self._m[0, 2])
def mnd(self, value): self._m[0, 2] = value self._m[2, 0] = value self._update()
def med(self): return float(self._m[1, 2])
def med(self, value): self._m[1, 2] = value self._m[2, 1] = value self._update()
def strike1(self): return float(self.both_strike_dip_rake()[0][0])
def dip1(self): return float(self.both_strike_dip_rake()[0][1])
def rake1(self): return float(self.both_strike_dip_rake()[0][2])
def strike2(self): return float(self.both_strike_dip_rake()[1][0])
def dip2(self): return float(self.both_strike_dip_rake()[1][1])
def rake2(self): return float(self.both_strike_dip_rake()[1][2])
'''Get both possible (strike,dip,rake) triplets.''' results = [] for rotmat in self._rotmats: alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)] results.append((beta, alpha, -gamma))
return results
'''Get direction of p axis.''' return (self._m_eigenvecs.T)[0]
'''Get direction of t axis.''' return (self._m_eigenvecs.T)[2]
'''Get diretion of the null axis.''' return self._m_eigenvecs.T[1]
''' Get the eigenvalues of the moment tensor in accending order.
:returns: ``(ep, en, et)`` '''
return self._m_eigenvals
''' Get the eigenvalues and eigenvectors of the moment tensor.
:returns: ``(ep, en, et, vp, vn, vt)``'''
vp = self.p_axis().A.flatten() vn = self.null_axis().A.flatten() vt = self.t_axis().A.flatten() ep, en, et = self._m_eigenvals return ep, en, et, vp, vn, vt
'''Get both possible slip directions.''' return [rotmat*cvec(1., 0., 0.) for rotmat in self._rotmats]
'''Get plain moment tensor as 3x3 matrix.''' return self._m.copy()
''' Get the moment tensor as a six-element array.
:returns: ``(mnn, mee, mdd, mne, mnd, med)`` ''' return to6(self._m)
'''Get moment tensor in up-south-east convention as 3x3 matrix.
.. math:: :nowrap:
\\begin{align*} M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\ M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\ M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne} \\end{align*} '''
return self._to_up_south_east.T * self._m * self._to_up_south_east
'''Get moment tensor in up-south-east convention as a six-element array.
:returns: ``(muu, mss, mee, mus, mue, mse)`` ''' return to6(self.m_up_south_east())
'''Get plain double couple with same scalar moment as moment tensor.''' rotmat1 = self._rotmats[0] m = rotmat1.T * MomentTensor._m_unrot * rotmat1 * self.scalar_moment() return m
'''Get moment magnitude of moment tensor.''' return moment_to_magnitude(self.scalar_moment())
''' Get the scalar moment of the moment tensor (Frobenius norm based)
.. math::
M0 = \\frac{1}{\\sqrt{2}}\\sqrt{\\sum_{i,j} |M_{ij}|^2}
The scalar moment is calculated based on the Euclidean (Frobenius) norm (Silver and Jordan, 1982). The scalar moment returned by this function differs from the standard decomposition based definition of the scalar moment for non-double-couple moment tensors. ''' return num.linalg.norm(self._m_eigenvals) / math.sqrt(2.)
def moment(self): return float(self.scalar_moment())
def moment(self, value): self._m *= value / self.moment self._update()
def magnitude(self): return float(self.moment_magnitude())
def magnitude(self, value): self._m *= magnitude_to_moment(value) / self.moment self._update()
def __str__(self): mexp = pow(10, math.ceil(num.log10(num.max(num.abs(self._m))))) m = self._m / mexp s = '''Scalar Moment [Nm]: M0 = %g (Mw = %3.1f) Moment Tensor [Nm]: Mnn = %6.3f, Mee = %6.3f, Mdd = %6.3f, Mne = %6.3f, Mnd = %6.3f, Med = %6.3f [ x %g ] ''' % ( self.scalar_moment(), self.moment_magnitude(), m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2], mexp)
s += self.str_fault_planes() return s
s = '' for i, sdr in enumerate(self.both_strike_dip_rake()): s += 'Fault plane %i [deg]: ' \ 'strike = %3.0f, dip = %3.0f, slip-rake = %4.0f\n' \ % (i+1, sdr[0], sdr[1], sdr[2])
return s
''' Get deviatoric part of moment tensor.
Returns a new moment tensor object with zero trace. '''
m = self.m()
trace_m = num.trace(m) m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.]) m_devi = m - m_iso mt = MomentTensor(m=m_devi) return mt
'''Decompose moment tensor into isotropic, DC and CLVD components.
Standard decomposition according to e.g. Jost and Herrmann 1989 is returned as::
[ (moment_iso, ratio_iso, m_iso), (moment_dc, ratio_dc, m_dc), (moment_clvd, ratio_clvd, m_clvd), (moment_devi, ratio_devi, m_devi), (moment, 1.0, m) ] '''
epsilon = 1e-6
m = self.m()
trace_m = num.trace(m) m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.]) moment_iso = abs(trace_m / 3.)
m_devi = m - m_iso
evals, evecs = eigh_check(m_devi)
moment_devi = num.max(num.abs(evals)) moment = moment_iso + moment_devi
iorder = num.argsort(num.abs(evals)) evals_sorted = evals[iorder] evecs_sorted = (evecs.T[iorder]).T
if moment_devi < epsilon * moment_iso: signed_moment_dc = 0. else: assert -epsilon <= -evals_sorted[0] / evals_sorted[2] <= 0.5 signed_moment_dc = evals_sorted[2] * (1.0 + 2.0 * ( min(0.0, evals_sorted[0] / evals_sorted[2])))
moment_dc = abs(signed_moment_dc) m_dc_es = signed_moment_dc * num.diag([0., -1.0, 1.0]) m_dc = num.dot(evecs_sorted, num.dot(m_dc_es, evecs_sorted.T))
m_clvd = m_devi - m_dc
moment_clvd = moment_devi - moment_dc
ratio_dc = moment_dc / moment ratio_clvd = moment_clvd / moment ratio_iso = moment_iso / moment ratio_devi = moment_devi / moment
return [ (moment_iso, ratio_iso, m_iso), (moment_dc, ratio_dc, m_dc), (moment_clvd, ratio_clvd, m_clvd), (moment_devi, ratio_devi, m_devi), (moment, 1.0, m)]
''' Get rotated moment tensor.
:param rot: ratation matrix, coordinate system is NED :returns: new :py:class:`MomentTensor` object '''
rotmat = num.matrix(rot) return MomentTensor(m=rotmat * self.m() * rotmat.T)
''' Get distorted MT by rotation around random axis and angle.
:param angle_std: angles are drawn from a normal distribution with zero mean and given standard deviation [degrees] :param angle: set angle [degrees], only axis will be random :param rstate: :py:class:`numpy.random.RandomState` object, can be used to create reproducible pseudo-random sequences :returns: new :py:class:`MomentTensor` object '''
assert (angle_std is None) != (angle is None), \ 'either angle or angle_std must be given'
if angle_std is not None: rstate = rstate or num.random angle = rstate.normal(scale=angle_std)
axis = random_axis(rstate=rstate) rot = rotation_from_angle_and_axis(angle, axis) return self.rotated(rot)
''' Get the respectively other plane in the double-couple ambiguity. '''
mt = MomentTensor(strike=strike, dip=dip, rake=rake) both_sdr = mt.both_strike_dip_rake() w = [sum([abs(x-y) for x, y in zip(both_sdr[i], (strike, dip, rake))]) for i in (0, 1)]
if w[0] < w[1]: return both_sdr[1] else: return both_sdr[0]
s1, d1, r1 = sdr1 s2, d2, r2 = sdr2
s1 = s1 % 360. s2 = s2 % 360. r1 = r1 % 360. r2 = r2 % 360.
ds = abs(s1 - s2) ds = ds if ds <= 180. else 360. - ds
dr = abs(r1 - r2) dr = dr if dr <= 180. else 360. - dr
dd = abs(d1 - d2)
return math.sqrt(ds**2 + dr**2 + dd**2)
''' Order strike-dip-rake pair post closely to a given reference pair.
:param sdrs: tuple, ``((strike1, dip1, rake1), (strike2, dip2, rake2))`` :param sdrs_ref: as above but with reference pair '''
d1 = min(dsdr(sdrs[0], sdrs_ref[0]), dsdr(sdrs[1], sdrs_ref[1])) d2 = min(dsdr(sdrs[0], sdrs_ref[1]), dsdr(sdrs[1], sdrs_ref[0])) if d1 < d2: return sdrs else: return sdrs[::-1]
eps = 0.001 tqw = 1. + t[0] + p[1] + b[2] tqx = 1. + t[0] - p[1] - b[2] tqy = 1. - t[0] + p[1] - b[2] tqz = 1. - t[0] - p[1] + b[2]
q = num.zeros(4) if tqw > eps: q[0] = 0.5 * math.sqrt(tqw) q[1:] = p[2] - b[1], b[0] - t[2], t[1] - p[0] elif tqx > eps: q[0] = 0.5 * math.sqrt(tqx) q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2] elif tqy > eps: q[0] = 0.5 * math.sqrt(tqy) q[1:] = b[0] - t[2], p[0] + t[1], b[1] + p[2] elif tqz > eps: q[0] = 0.5 * math.sqrt(tqz) q[1:] = t[1] - p[0], b[0] + t[2], b[1] + p[2] else: raise Exception('should not reach this line, check theory!') # q[0] = max(0.5 * math.sqrt(tqx), eps) # q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
q[1:] /= 4.0 * q[0]
q /= math.sqrt(num.sum(q**2))
return q
''' Given two moment tensors, return the Kagan angle in degrees.
After Kagan (1991) and Tape & Tape (2012). '''
ai = _pbt2tpb * mt1._m_eigenvecs.T aj = _pbt2tpb * mt2._m_eigenvecs.T
u = ai * aj.T
tk, pk, bk = u.A
qk = _tpb2q(tk, pk, bk)
return 2. * r2d * math.acos(num.max(num.abs(qk)))
''' Draw magnitude from Gutenberg Richter distribution. ''' return magnitude_min + num.log10(1.-rand) / -b_value
''' Scaling relation used by Global CMT catalog for most of its events. '''
mom = magnitude_to_moment(magnitudes) return (mom / 1.1e16)**(1./3.)
if __name__ == '__main__':
import sys v = list(map(float, sys.argv[1:])) mt = MomentTensor.from_values(v) print(mt) |