1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7This module provides various moment tensor related utility functions.
9It can be used to convert between strike-dip-rake and moment tensor
10representations and provides different options to produce random moment
11tensors.
13Moment tensors are represented by :py:class:`MomentTensor` instances. The
14internal representation uses a north-east-down (NED) coordinate system, but it
15can convert from/to the conventions used by the Global CMT catalog
16(up-south-east, USE).
18If not otherwise noted, scalar moment is interpreted as the Frobenius norm
19based scalar moment (see :py:meth:`MomentTensor.scalar_moment`. The scalar
20moment according to the "standard decomposition" can be found in the
21output of :py:meth:`MomentTensor.standard_decomposition`.
22'''
24from __future__ import division, print_function, absolute_import
26import math
27import numpy as num
29from .guts import Object, Float
31guts_prefix = 'pf'
33dynecm = 1e-7
36def random_axis(rstate=None):
37 '''
38 Get randomly oriented unit vector.
40 :param rstate: :py:class:`numpy.random.RandomState` object, can be used to
41 create reproducible pseudo-random sequences
42 '''
43 rstate = rstate or num.random
44 while True:
45 axis = rstate.uniform(size=3) * 2.0 - 1.0
46 uabs = math.sqrt(num.sum(axis**2))
47 if 0.001 < uabs < 1.0:
48 return axis / uabs
51def rotation_from_angle_and_axis(angle, axis):
52 '''
53 Build rotation matrix based on axis and angle.
55 :param angle: rotation angle [degrees]
56 :param axis: orientation of rotation axis, either in spherical
57 coordinates ``(theta, phi)`` [degrees], or as a unit vector
58 ``(ux, uy, uz)``.
59 '''
61 if len(axis) == 2:
62 theta, phi = axis
63 ux = math.sin(d2r*theta)*math.cos(d2r*phi)
64 uy = math.sin(d2r*theta)*math.sin(d2r*phi)
65 uz = math.cos(d2r*theta)
67 elif len(axis) == 3:
68 axis = num.asarray(axis)
69 uabs = math.sqrt(num.sum(axis**2))
70 ux, uy, uz = axis / uabs
71 else:
72 assert False
74 ct = math.cos(d2r*angle)
75 st = math.sin(d2r*angle)
76 return num.matrix([
77 [ct + ux**2*(1.-ct), ux*uy*(1.-ct)-uz*st, ux*uz*(1.-ct)+uy*st],
78 [uy*ux*(1.-ct)+uz*st, ct+uy**2*(1.-ct), uy*uz*(1.-ct)-ux*st],
79 [uz*ux*(1.-ct)-uy*st, uz*uy*(1.-ct)+ux*st, ct+uz**2*(1.-ct)]
80 ])
83def random_rotation(x=None):
84 '''
85 Get random rotation matrix.
87 A random rotation matrix, drawn from a uniform distrubution in the space
88 of rotations is returned, after Avro 1992 - "Fast random rotation
89 matrices".
91 :param x: three (uniform random) numbers in the range [0, 1[ used as input
92 to the distribution tranformation. If ``None``, random numbers are
93 used. Can be used to create grids of random rotations with uniform
94 density in rotation space.
95 '''
97 if x is not None:
98 x1, x2, x3 = x
99 else:
100 x1, x2, x3 = num.random.random(3)
102 phi = math.pi*2.0*x1
104 zrot = num.matrix([
105 [math.cos(phi), math.sin(phi), 0.],
106 [-math.sin(phi), math.cos(phi), 0.],
107 [0., 0., 1.]])
109 lam = math.pi*2.0*x2
111 v = num.matrix([[
112 math.cos(lam)*math.sqrt(x3),
113 math.sin(lam)*math.sqrt(x3),
114 math.sqrt(1.-x3)]]).T
116 house = num.identity(3) - 2.0 * v * v.T
117 return -house*zrot
120def rand(mi, ma):
121 return float(num.random.uniform(mi, ma))
124def randdip(mi, ma):
125 mi_ = 0.5*(math.cos(mi * math.pi/180.)+1.)
126 ma_ = 0.5*(math.cos(ma * math.pi/180.)+1.)
127 return math.acos(rand(mi_, ma_)*2.-1.)*180./math.pi
130def random_strike_dip_rake(
131 strikemin=0., strikemax=360.,
132 dipmin=0., dipmax=90.,
133 rakemin=-180., rakemax=180.):
135 '''
136 Get random strike, dip, rake triplet.
138 .. note::
140 Might not produce a homogeneous distribution of mechanisms. Better use
141 :py:meth:`MomentTensor.random_dc` which is based on
142 :py:func:`random_rotation`.
143 '''
145 strike = rand(strikemin, strikemax)
146 dip = randdip(dipmin, dipmax)
147 rake = rand(rakemin, rakemax)
149 return strike, dip, rake
152def to6(m):
153 '''
154 Get non-redundant components from symmetric 3x3 matrix.
156 :returns: 1D NumPy array with entries ordered like
157 ``(a_xx, a_yy, a_zz, a_xy, a_xz, a_yz)``
158 '''
160 return num.array([m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2]])
163def symmat6(a_xx, a_yy, a_zz, a_xy, a_xz, a_yz):
164 '''
165 Create symmetric 3x3 matrix from its 6 non-redundant values.
166 '''
168 return num.matrix([[a_xx, a_xy, a_xz],
169 [a_xy, a_yy, a_yz],
170 [a_xz, a_yz, a_zz]], dtype=float)
173def values_to_matrix(values):
174 '''
175 Convert anything to moment tensor represented as a NumPy matrix.
177 Transforms :py:class:`MomentTensor` objects, tuples, lists and NumPy arrays
178 with 3x3 or 3, 4, 6, or 7 elements into NumPy 3x3 matrix objects.
180 The ``values`` argument is interpreted depending on shape and type as
181 follows:
183 * ``(strike, dip, rake)``
184 * ``(strike, dip, rake, magnitude)``
185 * ``(mnn, mee, mdd, mne, mnd, med)``
186 * ``(mnn, mee, mdd, mne, mnd, med, magnitude)``
187 * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))``
188 * :py:class:`MomentTensor`
189 '''
191 if isinstance(values, (tuple, list)):
192 values = num.asarray(values, dtype=float)
194 if isinstance(values, MomentTensor):
195 return values.m()
197 elif isinstance(values, num.ndarray):
198 if values.shape == (3,):
199 strike, dip, rake = values
200 rotmat1 = euler_to_matrix(d2r*dip, d2r*strike, -d2r*rake)
201 return rotmat1.T * MomentTensor._m_unrot * rotmat1
203 elif values.shape == (4,):
204 strike, dip, rake, magnitude = values
205 moment = magnitude_to_moment(magnitude)
206 rotmat1 = euler_to_matrix(d2r*dip, d2r*strike, -d2r*rake)
207 return rotmat1.T * MomentTensor._m_unrot * rotmat1 * moment
209 elif values.shape == (6,):
210 return symmat6(*values)
212 elif values.shape == (7,):
213 magnitude = values[6]
214 moment = magnitude_to_moment(magnitude)
215 mt = symmat6(*values[:6])
216 mt *= moment / (num.linalg.norm(mt) / math.sqrt(2.0))
217 return mt
219 elif values.shape == (3, 3):
220 return num.asmatrix(values, dtype=float)
222 raise Exception('cannot convert object to 3x3 matrix')
225def moment_to_magnitude(moment):
226 '''
227 Convert scalar moment to moment magnitude Mw.
229 :param moment: scalar moment [Nm]
230 :returns: moment magnitude Mw
232 Moment magnitude is defined as
234 .. math::
236 M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0 \\cdot 10^7) - 10.7
238 where :math:`M_0` is the scalar moment given in [Nm].
240 .. note::
242 Global CMT uses 10.7333333 instead of 10.7, based on [Kanamori 1977],
243 10.7 is from [Hanks and Kanamori 1979].
244 '''
246 return num.log10(moment*1.0e7) / 1.5 - 10.7
249def magnitude_to_moment(magnitude):
250 '''
251 Convert moment magnitude Mw to scalar moment.
253 :param magnitude: moment magnitude
254 :returns: scalar moment [Nm]
256 See :py:func:`moment_to_magnitude`.
257 '''
259 return 10.0**(1.5*(magnitude+10.7))*1.0e-7
262magnitude_1Nm = moment_to_magnitude(1.0)
265def euler_to_matrix(alpha, beta, gamma):
266 '''
267 Given euler angle triplet, create rotation matrix
269 Given coordinate system `(x,y,z)` and rotated system `(xs,ys,zs)`
270 the line of nodes is the intersection between the `x,y` and the `xs,ys`
271 planes.
273 :param alpha: is the angle between the `z`-axis and the `zs`-axis [rad]
274 :param beta: is the angle between the `x`-axis and the line of nodes [rad]
275 :param gamma: is the angle between the line of nodes and the `xs`-axis
276 [rad]
278 Usage for moment tensors::
280 m_unrot = numpy.matrix([[0,0,-1],[0,0,0],[-1,0,0]])
281 euler_to_matrix(dip,strike,-rake, rotmat)
282 m = rotmat.T * m_unrot * rotmat
284 '''
286 ca = math.cos(alpha)
287 cb = math.cos(beta)
288 cg = math.cos(gamma)
289 sa = math.sin(alpha)
290 sb = math.sin(beta)
291 sg = math.sin(gamma)
293 mat = num.matrix([[cb*cg-ca*sb*sg, sb*cg+ca*cb*sg, sa*sg],
294 [-cb*sg-ca*sb*cg, -sb*sg+ca*cb*cg, sa*cg],
295 [sa*sb, -sa*cb, ca]], dtype=float)
296 return mat
299def matrix_to_euler(rotmat):
300 '''
301 Get eulerian angle triplet from rotation matrix.
302 '''
304 ex = cvec(1., 0., 0.)
305 ez = cvec(0., 0., 1.)
306 exs = rotmat.T * ex
307 ezs = rotmat.T * ez
308 enodes = num.cross(ez.T, ezs.T).T
309 if num.linalg.norm(enodes) < 1e-10:
310 enodes = exs
311 enodess = rotmat*enodes
312 cos_alpha = float((ez.T*ezs))
313 if cos_alpha > 1.:
314 cos_alpha = 1.
316 if cos_alpha < -1.:
317 cos_alpha = -1.
319 alpha = math.acos(cos_alpha)
320 beta = num.mod(math.atan2(enodes[1, 0], enodes[0, 0]), math.pi*2.)
321 gamma = num.mod(-math.atan2(enodess[1, 0], enodess[0, 0]), math.pi*2.)
323 return unique_euler(alpha, beta, gamma)
326def unique_euler(alpha, beta, gamma):
327 '''
328 Uniquify eulerian angle triplet.
330 Put eulerian angle triplet into ranges compatible with
331 ``(dip, strike, -rake)`` conventions in seismology::
333 alpha (dip) : [0, pi/2]
334 beta (strike) : [0, 2*pi)
335 gamma (-rake) : [-pi, pi)
337 If ``alpha1`` is near to zero, ``beta`` is replaced by ``beta+gamma`` and
338 ``gamma`` is set to zero, to prevent this additional ambiguity.
340 If ``alpha`` is near to ``pi/2``, ``beta`` is put into the range
341 ``[0,pi)``.
342 '''
344 pi = math.pi
346 alpha = num.mod(alpha, 2.0*pi)
348 if 0.5*pi < alpha and alpha <= pi:
349 alpha = pi - alpha
350 beta = beta + pi
351 gamma = 2.0*pi - gamma
352 elif pi < alpha and alpha <= 1.5*pi:
353 alpha = alpha - pi
354 gamma = pi - gamma
355 elif 1.5*pi < alpha and alpha <= 2.0*pi:
356 alpha = 2.0*pi - alpha
357 beta = beta + pi
358 gamma = pi + gamma
360 alpha = num.mod(alpha, 2.0*pi)
361 beta = num.mod(beta, 2.0*pi)
362 gamma = num.mod(gamma+pi, 2.0*pi)-pi
364 # If dip is exactly 90 degrees, one is still
365 # free to choose between looking at the plane from either side.
366 # Choose to look at such that beta is in the range [0,180)
368 # This should prevent some problems, when dip is close to 90 degrees:
369 if abs(alpha - 0.5*pi) < 1e-10:
370 alpha = 0.5*pi
372 if abs(beta - pi) < 1e-10:
373 beta = pi
375 if abs(beta - 2.*pi) < 1e-10:
376 beta = 0.
378 if abs(beta) < 1e-10:
379 beta = 0.
381 if alpha == 0.5*pi and beta >= pi:
382 gamma = - gamma
383 beta = num.mod(beta-pi, 2.0*pi)
384 gamma = num.mod(gamma+pi, 2.0*pi)-pi
385 assert 0. <= beta < pi
386 assert -pi <= gamma < pi
388 if alpha < 1e-7:
389 beta = num.mod(beta + gamma, 2.0*pi)
390 gamma = 0.
392 return (alpha, beta, gamma)
395def cvec(x, y, z):
396 return num.matrix([[x, y, z]], dtype=float).T
399def rvec(x, y, z):
400 return num.matrix([[x, y, z]], dtype=float)
403def eigh_check(a):
404 evals, evecs = num.linalg.eigh(a)
405 assert evals[0] <= evals[1] <= evals[2]
406 return evals, evecs
409r2d = 180. / math.pi
410d2r = 1. / r2d
413def random_mt(x=None, scalar_moment=1.0, magnitude=None):
415 if magnitude is not None:
416 scalar_moment = magnitude_to_moment(magnitude)
418 if x is None:
419 x = num.random.random(6)
421 evals = x[:3] * 2. - 1.0
422 evals /= num.sqrt(num.sum(evals**2)) / math.sqrt(2.0)
423 rotmat = random_rotation(x[3:])
424 return scalar_moment * rotmat * num.matrix(num.diag(evals)) * rotmat.T
427def random_m6(*args, **kwargs):
428 return to6(random_mt(*args, **kwargs))
431def random_dc(x=None, scalar_moment=1.0, magnitude=None):
432 if magnitude is not None:
433 scalar_moment = magnitude_to_moment(magnitude)
435 rotmat = random_rotation(x)
436 return scalar_moment * (rotmat * MomentTensor._m_unrot * rotmat.T)
439def sm(m):
440 return "/ %5.2F %5.2F %5.2F \\\n" % (m[0, 0], m[0, 1], m[0, 2]) + \
441 "| %5.2F %5.2F %5.2F |\n" % (m[1, 0], m[1, 1], m[1, 2]) + \
442 "\\ %5.2F %5.2F %5.2F /\n" % (m[2, 0], m[2, 1], m[2, 2])
445def as_mt(mt):
446 '''
447 Convenience function to convert various objects to moment tensor object.
449 Like :py:meth:``MomentTensor.from_values``, but does not create a new
450 :py:class:`MomentTensor` object when ``mt`` already is one.
451 '''
453 if isinstance(mt, MomentTensor):
454 return mt
455 else:
456 return MomentTensor.from_values(mt)
459class MomentTensor(Object):
461 '''
462 Moment tensor object
464 :param m: NumPy matrix in north-east-down convention
465 :param m_up_south_east: NumPy matrix in up-south-east convention
466 :param m_east_north_up: NumPy matrix in east-north-up convention
467 :param strike,dip,rake: fault plane angles in [degrees]
468 :param scalar_moment: scalar moment in [Nm]
469 :param magnitude: moment magnitude Mw
471 Global CMT catalog moment tensors use the up-south-east (USE) coordinate
472 system convention with :math:`r` (up), :math:`\\theta` (south), and
473 :math:`\\phi` (east).
475 .. math::
476 :nowrap:
478 \\begin{align*}
479 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
480 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
481 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
482 \\end{align*}
484 '''
486 mnn__ = Float.T(default=0.0)
487 mee__ = Float.T(default=0.0)
488 mdd__ = Float.T(default=0.0)
489 mne__ = Float.T(default=0.0)
490 mnd__ = Float.T(default=-1.0)
491 med__ = Float.T(default=0.0)
492 strike1__ = Float.T(default=None, optional=True) # read-only
493 dip1__ = Float.T(default=None, optional=True) # read-only
494 rake1__ = Float.T(default=None, optional=True) # read-only
495 strike2__ = Float.T(default=None, optional=True) # read-only
496 dip2__ = Float.T(default=None, optional=True) # read-only
497 rake2__ = Float.T(default=None, optional=True) # read-only
498 moment__ = Float.T(default=None, optional=True) # read-only
499 magnitude__ = Float.T(default=None, optional=True) # read-only
501 _flip_dc = num.matrix(
502 [[0., 0., -1.], [0., -1., 0.], [-1., 0., 0.]], dtype=float)
503 _to_up_south_east = num.matrix(
504 [[0., 0., -1.], [-1., 0., 0.], [0., 1., 0.]], dtype=float).T
505 _to_east_north_up = num.matrix(
506 [[0., 1., 0.], [1., 0., 0.], [0., 0., -1.]], dtype=float)
507 _m_unrot = num.matrix(
508 [[0., 0., -1.], [0., 0., 0.], [-1., 0., 0.]], dtype=float)
510 _u_evals, _u_evecs = eigh_check(_m_unrot)
512 @classmethod
513 def random_dc(cls, x=None, scalar_moment=1.0, magnitude=None):
514 '''
515 Create random oriented double-couple moment tensor
517 The rotations used are uniformly distributed in the space of rotations.
518 '''
519 return MomentTensor(
520 m=random_dc(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
522 @classmethod
523 def random_mt(cls, x=None, scalar_moment=1.0, magnitude=None):
524 '''
525 Create random moment tensor
527 Moment tensors produced by this function appear uniformly distributed
528 when shown in a Hudson's diagram. The rotations used are unifomly
529 distributed in the space of rotations.
530 '''
531 return MomentTensor(
532 m=random_mt(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
534 @classmethod
535 def from_values(cls, values):
536 '''
537 Alternative constructor for moment tensor objects
539 This constructor takes a :py:class:`MomentTensor` object, a tuple, list
540 or NumPy array with 3x3 or 3, 4, 6, or 7 elements to build a Moment
541 tensor object.
543 The ``values`` argument is interpreted depending on shape and type as
544 follows:
546 * ``(strike, dip, rake)``
547 * ``(strike, dip, rake, magnitude)``
548 * ``(mnn, mee, mdd, mne, mnd, med)``
549 * ``(mnn, mee, mdd, mne, mnd, med, magnitude)``
550 * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))``
551 * :py:class:`MomentTensor` object
552 '''
554 m = values_to_matrix(values)
555 return MomentTensor(m=m)
557 def __init__(
558 self, m=None, m_up_south_east=None, m_east_north_up=None,
559 strike=0., dip=0., rake=0., scalar_moment=1.,
560 mnn=None, mee=None, mdd=None, mne=None, mnd=None, med=None,
561 strike1=None, dip1=None, rake1=None,
562 strike2=None, dip2=None, rake2=None,
563 magnitude=None, moment=None):
565 Object.__init__(self, init_props=False)
567 if any(mxx is not None for mxx in (mnn, mee, mdd, mne, mnd, med)):
568 m = symmat6(mnn, mee, mdd, mne, mnd, med)
570 strike = d2r*strike
571 dip = d2r*dip
572 rake = d2r*rake
574 if m_up_south_east is not None:
575 m = self._to_up_south_east * m_up_south_east * \
576 self._to_up_south_east.T
578 if m_east_north_up is not None:
579 m = self._to_east_north_up * m_east_north_up * \
580 self._to_east_north_up.T
582 if m is None:
583 if any(x is not None for x in (
584 strike1, dip1, rake1, strike2, dip2, rake2)):
585 raise Exception(
586 'strike1, dip1, rake1, strike2, dip2, rake2 are read-only '
587 'properties')
589 if moment is not None:
590 scalar_moment = moment
592 if magnitude is not None:
593 scalar_moment = magnitude_to_moment(magnitude)
595 rotmat1 = euler_to_matrix(dip, strike, -rake)
596 m = rotmat1.T * MomentTensor._m_unrot * rotmat1 * scalar_moment
598 self._m = m
599 self._update()
601 def _update(self):
602 m_evals, m_evecs = eigh_check(self._m)
603 if num.linalg.det(m_evecs) < 0.:
604 m_evecs *= -1.
606 rotmat1 = (m_evecs * MomentTensor._u_evecs.T).T
607 if num.linalg.det(rotmat1) < 0.:
608 rotmat1 *= -1.
610 self._m_eigenvals = m_evals
611 self._m_eigenvecs = m_evecs
613 self._rotmats = sorted(
614 [rotmat1, MomentTensor._flip_dc * rotmat1],
615 key=lambda m: num.abs(m.flat).tolist())
617 @property
618 def mnn(self):
619 return float(self._m[0, 0])
621 @mnn.setter
622 def mnn(self, value):
623 self._m[0, 0] = value
624 self._update()
626 @property
627 def mee(self):
628 return float(self._m[1, 1])
630 @mee.setter
631 def mee(self, value):
632 self._m[1, 1] = value
633 self._update()
635 @property
636 def mdd(self):
637 return float(self._m[2, 2])
639 @mdd.setter
640 def mdd(self, value):
641 self._m[2, 2] = value
642 self._update()
644 @property
645 def mne(self):
646 return float(self._m[0, 1])
648 @mne.setter
649 def mne(self, value):
650 self._m[0, 1] = value
651 self._m[1, 0] = value
652 self._update()
654 @property
655 def mnd(self):
656 return float(self._m[0, 2])
658 @mnd.setter
659 def mnd(self, value):
660 self._m[0, 2] = value
661 self._m[2, 0] = value
662 self._update()
664 @property
665 def med(self):
666 return float(self._m[1, 2])
668 @med.setter
669 def med(self, value):
670 self._m[1, 2] = value
671 self._m[2, 1] = value
672 self._update()
674 @property
675 def strike1(self):
676 return float(self.both_strike_dip_rake()[0][0])
678 @property
679 def dip1(self):
680 return float(self.both_strike_dip_rake()[0][1])
682 @property
683 def rake1(self):
684 return float(self.both_strike_dip_rake()[0][2])
686 @property
687 def strike2(self):
688 return float(self.both_strike_dip_rake()[1][0])
690 @property
691 def dip2(self):
692 return float(self.both_strike_dip_rake()[1][1])
694 @property
695 def rake2(self):
696 return float(self.both_strike_dip_rake()[1][2])
698 def both_strike_dip_rake(self):
699 '''
700 Get both possible (strike,dip,rake) triplets.
701 '''
702 results = []
703 for rotmat in self._rotmats:
704 alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)]
705 results.append((beta, alpha, -gamma))
707 return results
709 def p_axis(self):
710 '''
711 Get direction of p axis.
712 '''
713 return (self._m_eigenvecs.T)[0]
715 def t_axis(self):
716 '''
717 Get direction of t axis.
718 '''
719 return (self._m_eigenvecs.T)[2]
721 def null_axis(self):
722 '''
723 Get diretion of the null axis.
724 '''
725 return self._m_eigenvecs.T[1]
727 def eigenvals(self):
728 '''
729 Get the eigenvalues of the moment tensor in accending order.
731 :returns: ``(ep, en, et)``
732 '''
734 return self._m_eigenvals
736 def eigensystem(self):
737 '''
738 Get the eigenvalues and eigenvectors of the moment tensor.
740 :returns: ``(ep, en, et, vp, vn, vt)``
741 '''
743 vp = self.p_axis().A.flatten()
744 vn = self.null_axis().A.flatten()
745 vt = self.t_axis().A.flatten()
746 ep, en, et = self._m_eigenvals
747 return ep, en, et, vp, vn, vt
749 def both_slip_vectors(self):
750 '''
751 Get both possible slip directions.
752 '''
753 return [rotmat*cvec(1., 0., 0.) for rotmat in self._rotmats]
755 def m(self):
756 '''
757 Get plain moment tensor as 3x3 matrix.
758 '''
759 return self._m.copy()
761 def m6(self):
762 '''
763 Get the moment tensor as a six-element array.
765 :returns: ``(mnn, mee, mdd, mne, mnd, med)``
766 '''
767 return to6(self._m)
769 def m_up_south_east(self):
770 '''
771 Get moment tensor in up-south-east convention as 3x3 matrix.
773 .. math::
774 :nowrap:
776 \\begin{align*}
777 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
778 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
779 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
780 \\end{align*}
781 '''
783 return self._to_up_south_east.T * self._m * self._to_up_south_east
785 def m_east_north_up(self):
786 '''
787 Get moment tensor in east-north-up convention as 3x3 matrix.
788 '''
790 return self._to_east_north_up.T * self._m * self._to_east_north_up
792 def m6_up_south_east(self):
793 '''
794 Get moment tensor in up-south-east convention as a six-element array.
796 :returns: ``(muu, mss, mee, mus, mue, mse)``
797 '''
798 return to6(self.m_up_south_east())
800 def m6_east_north_up(self):
801 '''
802 Get moment tensor in east-north-up convention as a six-element array.
804 :returns: ``(mee, mnn, muu, men, meu, mnu)``
805 '''
806 return to6(self.m_east_north_up())
808 def m_plain_double_couple(self):
809 '''
810 Get plain double couple with same scalar moment as moment tensor.
811 '''
812 rotmat1 = self._rotmats[0]
813 m = rotmat1.T * MomentTensor._m_unrot * rotmat1 * self.scalar_moment()
814 return m
816 def moment_magnitude(self):
817 '''
818 Get moment magnitude of moment tensor.
819 '''
820 return moment_to_magnitude(self.scalar_moment())
822 def scalar_moment(self):
823 '''
824 Get the scalar moment of the moment tensor (Frobenius norm based)
826 .. math::
828 M_0 = \\frac{1}{\\sqrt{2}}\\sqrt{\\sum_{i,j} |M_{ij}|^2}
830 The scalar moment is calculated based on the Euclidean (Frobenius) norm
831 (Silver and Jordan, 1982). The scalar moment returned by this function
832 differs from the standard decomposition based definition of the scalar
833 moment for non-double-couple moment tensors.
834 '''
835 return num.linalg.norm(self._m_eigenvals) / math.sqrt(2.)
837 @property
838 def moment(self):
839 return float(self.scalar_moment())
841 @moment.setter
842 def moment(self, value):
843 self._m *= value / self.moment
844 self._update()
846 @property
847 def magnitude(self):
848 return float(self.moment_magnitude())
850 @magnitude.setter
851 def magnitude(self, value):
852 self._m *= magnitude_to_moment(value) / self.moment
853 self._update()
855 def __str__(self):
856 mexp = pow(10, math.ceil(num.log10(num.max(num.abs(self._m)))))
857 m = self._m / mexp
858 s = '''Scalar Moment [Nm]: M0 = %g (Mw = %3.1f)
859Moment Tensor [Nm]: Mnn = %6.3f, Mee = %6.3f, Mdd = %6.3f,
860 Mne = %6.3f, Mnd = %6.3f, Med = %6.3f [ x %g ]
861''' % (
862 self.scalar_moment(),
863 self.moment_magnitude(),
864 m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2],
865 mexp)
867 s += self.str_fault_planes()
868 return s
870 def str_fault_planes(self):
871 s = ''
872 for i, sdr in enumerate(self.both_strike_dip_rake()):
873 s += 'Fault plane %i [deg]: ' \
874 'strike = %3.0f, dip = %3.0f, slip-rake = %4.0f\n' \
875 % (i+1, sdr[0], sdr[1], sdr[2])
877 return s
879 def deviatoric(self):
880 '''
881 Get deviatoric part of moment tensor.
883 Returns a new moment tensor object with zero trace.
884 '''
886 m = self.m()
888 trace_m = num.trace(m)
889 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
890 m_devi = m - m_iso
891 mt = MomentTensor(m=m_devi)
892 return mt
894 def standard_decomposition(self):
895 '''
896 Decompose moment tensor into isotropic, DC and CLVD components.
898 Standard decomposition according to e.g. Jost and Herrmann 1989 is
899 returned as::
901 [
902 (moment_iso, ratio_iso, m_iso),
903 (moment_dc, ratio_dc, m_dc),
904 (moment_clvd, ratio_clvd, m_clvd),
905 (moment_devi, ratio_devi, m_devi),
906 (moment, 1.0, m)
907 ]
908 '''
910 epsilon = 1e-6
912 m = self.m()
914 trace_m = num.trace(m)
915 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
916 moment_iso = abs(trace_m / 3.)
918 m_devi = m - m_iso
920 evals, evecs = eigh_check(m_devi)
922 moment_devi = num.max(num.abs(evals))
923 moment = moment_iso + moment_devi
925 iorder = num.argsort(num.abs(evals))
926 evals_sorted = evals[iorder]
927 evecs_sorted = (evecs.T[iorder]).T
929 if moment_devi < epsilon * moment_iso:
930 signed_moment_dc = 0.
931 else:
932 assert -epsilon <= -evals_sorted[0] / evals_sorted[2] <= 0.5
933 signed_moment_dc = evals_sorted[2] * (1.0 + 2.0 * (
934 min(0.0, evals_sorted[0] / evals_sorted[2])))
936 moment_dc = abs(signed_moment_dc)
937 m_dc_es = signed_moment_dc * num.diag([0., -1.0, 1.0])
938 m_dc = num.dot(evecs_sorted, num.dot(m_dc_es, evecs_sorted.T))
940 m_clvd = m_devi - m_dc
942 moment_clvd = moment_devi - moment_dc
944 ratio_dc = moment_dc / moment
945 ratio_clvd = moment_clvd / moment
946 ratio_iso = moment_iso / moment
947 ratio_devi = moment_devi / moment
949 return [
950 (moment_iso, ratio_iso, m_iso),
951 (moment_dc, ratio_dc, m_dc),
952 (moment_clvd, ratio_clvd, m_clvd),
953 (moment_devi, ratio_devi, m_devi),
954 (moment, 1.0, m)]
956 def rotated(self, rot):
957 '''
958 Get rotated moment tensor.
960 :param rot: ratation matrix, coordinate system is NED
961 :returns: new :py:class:`MomentTensor` object
962 '''
964 rotmat = num.matrix(rot)
965 return MomentTensor(m=rotmat * self.m() * rotmat.T)
967 def random_rotated(self, angle_std=None, angle=None, rstate=None):
968 '''
969 Get distorted MT by rotation around random axis and angle.
971 :param angle_std: angles are drawn from a normal distribution with
972 zero mean and given standard deviation [degrees]
973 :param angle: set angle [degrees], only axis will be random
974 :param rstate: :py:class:`numpy.random.RandomState` object, can be
975 used to create reproducible pseudo-random sequences
976 :returns: new :py:class:`MomentTensor` object
977 '''
979 assert (angle_std is None) != (angle is None), \
980 'either angle or angle_std must be given'
982 if angle_std is not None:
983 rstate = rstate or num.random
984 angle = rstate.normal(scale=angle_std)
986 axis = random_axis(rstate=rstate)
987 rot = rotation_from_angle_and_axis(angle, axis)
988 return self.rotated(rot)
991def other_plane(strike, dip, rake):
992 '''
993 Get the respectively other plane in the double-couple ambiguity.
994 '''
996 mt = MomentTensor(strike=strike, dip=dip, rake=rake)
997 both_sdr = mt.both_strike_dip_rake()
998 w = [sum([abs(x-y) for x, y in zip(both_sdr[i], (strike, dip, rake))])
999 for i in (0, 1)]
1001 if w[0] < w[1]:
1002 return both_sdr[1]
1003 else:
1004 return both_sdr[0]
1007def dsdr(sdr1, sdr2):
1008 s1, d1, r1 = sdr1
1009 s2, d2, r2 = sdr2
1011 s1 = s1 % 360.
1012 s2 = s2 % 360.
1013 r1 = r1 % 360.
1014 r2 = r2 % 360.
1016 ds = abs(s1 - s2)
1017 ds = ds if ds <= 180. else 360. - ds
1019 dr = abs(r1 - r2)
1020 dr = dr if dr <= 180. else 360. - dr
1022 dd = abs(d1 - d2)
1024 return math.sqrt(ds**2 + dr**2 + dd**2)
1027def order_like(sdrs, sdrs_ref):
1028 '''
1029 Order strike-dip-rake pair post closely to a given reference pair.
1031 :param sdrs: tuple, ``((strike1, dip1, rake1), (strike2, dip2, rake2))``
1032 :param sdrs_ref: as above but with reference pair
1033 '''
1035 d1 = min(dsdr(sdrs[0], sdrs_ref[0]), dsdr(sdrs[1], sdrs_ref[1]))
1036 d2 = min(dsdr(sdrs[0], sdrs_ref[1]), dsdr(sdrs[1], sdrs_ref[0]))
1037 if d1 < d2:
1038 return sdrs
1039 else:
1040 return sdrs[::-1]
1043def _tpb2q(t, p, b):
1044 eps = 0.001
1045 tqw = 1. + t[0] + p[1] + b[2]
1046 tqx = 1. + t[0] - p[1] - b[2]
1047 tqy = 1. - t[0] + p[1] - b[2]
1048 tqz = 1. - t[0] - p[1] + b[2]
1050 q = num.zeros(4)
1051 if tqw > eps:
1052 q[0] = 0.5 * math.sqrt(tqw)
1053 q[1:] = p[2] - b[1], b[0] - t[2], t[1] - p[0]
1054 elif tqx > eps:
1055 q[0] = 0.5 * math.sqrt(tqx)
1056 q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1057 elif tqy > eps:
1058 q[0] = 0.5 * math.sqrt(tqy)
1059 q[1:] = b[0] - t[2], p[0] + t[1], b[1] + p[2]
1060 elif tqz > eps:
1061 q[0] = 0.5 * math.sqrt(tqz)
1062 q[1:] = t[1] - p[0], b[0] + t[2], b[1] + p[2]
1063 else:
1064 raise Exception('should not reach this line, check theory!')
1065 # q[0] = max(0.5 * math.sqrt(tqx), eps)
1066 # q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1068 q[1:] /= 4.0 * q[0]
1070 q /= math.sqrt(num.sum(q**2))
1072 return q
1075_pbt2tpb = num.matrix(((0., 0., 1.), (1., 0., 0.), (0., 1., 0.)))
1078def kagan_angle(mt1, mt2):
1079 '''
1080 Given two moment tensors, return the Kagan angle in degrees.
1082 After Kagan (1991) and Tape & Tape (2012).
1083 '''
1085 ai = _pbt2tpb * mt1._m_eigenvecs.T
1086 aj = _pbt2tpb * mt2._m_eigenvecs.T
1088 u = ai * aj.T
1090 tk, pk, bk = u.A
1092 qk = _tpb2q(tk, pk, bk)
1094 return 2. * r2d * math.acos(num.max(num.abs(qk)))
1097def rand_to_gutenberg_richter(rand, b_value, magnitude_min):
1098 '''
1099 Draw magnitude from Gutenberg Richter distribution.
1100 '''
1101 return magnitude_min + num.log10(1.-rand) / -b_value
1104def magnitude_to_duration_gcmt(magnitudes):
1105 '''
1106 Scaling relation used by Global CMT catalog for most of its events.
1107 '''
1109 mom = magnitude_to_moment(magnitudes)
1110 return (mom / 1.1e16)**(1./3.)
1113if __name__ == '__main__':
1115 import sys
1116 v = list(map(float, sys.argv[1:]))
1117 mt = MomentTensor.from_values(v)
1118 print(mt)