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.array([
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.array([
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.array([[
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 * num.dot(v, v.T)
117 return -num.dot(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.array([[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 array.
177 Transforms :py:class:`MomentTensor` objects, tuples, lists and NumPy arrays
178 with 3x3 or 3, 4, 6, or 7 elements into NumPy 3x3 array 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 num.dot(rotmat1.T, num.dot(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 moment * num.dot(
208 rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
210 elif values.shape == (6,):
211 return symmat6(*values)
213 elif values.shape == (7,):
214 magnitude = values[6]
215 moment = magnitude_to_moment(magnitude)
216 mt = symmat6(*values[:6])
217 mt *= moment / (num.linalg.norm(mt) / math.sqrt(2.0))
218 return mt
220 elif values.shape == (3, 3):
221 return num.asarray(values, dtype=float)
223 raise Exception('cannot convert object to 3x3 matrix')
226def moment_to_magnitude(moment):
227 '''
228 Convert scalar moment to moment magnitude Mw.
230 :param moment: scalar moment [Nm]
231 :returns: moment magnitude Mw
233 Moment magnitude is defined as
235 .. math::
237 M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0 \\cdot 10^7) - 10.7
239 where :math:`M_0` is the scalar moment given in [Nm].
241 .. note::
243 Global CMT uses 10.7333333 instead of 10.7, based on [Kanamori 1977],
244 10.7 is from [Hanks and Kanamori 1979].
245 '''
247 return num.log10(moment*1.0e7) / 1.5 - 10.7
250def magnitude_to_moment(magnitude):
251 '''
252 Convert moment magnitude Mw to scalar moment.
254 :param magnitude: moment magnitude
255 :returns: scalar moment [Nm]
257 See :py:func:`moment_to_magnitude`.
258 '''
260 return 10.0**(1.5*(magnitude+10.7))*1.0e-7
263magnitude_1Nm = moment_to_magnitude(1.0)
266def euler_to_matrix(alpha, beta, gamma):
267 '''
268 Given euler angle triplet, create rotation matrix
270 Given coordinate system `(x,y,z)` and rotated system `(xs,ys,zs)`
271 the line of nodes is the intersection between the `x,y` and the `xs,ys`
272 planes.
274 :param alpha: is the angle between the `z`-axis and the `zs`-axis [rad]
275 :param beta: is the angle between the `x`-axis and the line of nodes [rad]
276 :param gamma: is the angle between the line of nodes and the `xs`-axis
277 [rad]
279 Usage for moment tensors::
281 m_unrot = numpy.array([[0,0,-1],[0,0,0],[-1,0,0]])
282 euler_to_matrix(dip,strike,-rake, rotmat)
283 m = num.dot(rotmat.T, num.dot(m_unrot, rotmat))
285 '''
287 ca = math.cos(alpha)
288 cb = math.cos(beta)
289 cg = math.cos(gamma)
290 sa = math.sin(alpha)
291 sb = math.sin(beta)
292 sg = math.sin(gamma)
294 mat = num.array([[cb*cg-ca*sb*sg, sb*cg+ca*cb*sg, sa*sg],
295 [-cb*sg-ca*sb*cg, -sb*sg+ca*cb*cg, sa*cg],
296 [sa*sb, -sa*cb, ca]], dtype=float)
297 return mat
300def matrix_to_euler(rotmat):
301 '''
302 Get eulerian angle triplet from rotation matrix.
303 '''
305 ex = cvec(1., 0., 0.)
306 ez = cvec(0., 0., 1.)
307 exs = num.dot(rotmat.T, ex)
308 ezs = num.dot(rotmat.T, ez)
309 enodes = num.cross(ez.T, ezs.T).T
310 if num.linalg.norm(enodes) < 1e-10:
311 enodes = exs
312 enodess = num.dot(rotmat, enodes)
313 cos_alpha = float((num.dot(ez.T, ezs)))
314 if cos_alpha > 1.:
315 cos_alpha = 1.
317 if cos_alpha < -1.:
318 cos_alpha = -1.
320 alpha = math.acos(cos_alpha)
321 beta = num.mod(math.atan2(enodes[1, 0], enodes[0, 0]), math.pi*2.)
322 gamma = num.mod(-math.atan2(enodess[1, 0], enodess[0, 0]), math.pi*2.)
324 return unique_euler(alpha, beta, gamma)
327def unique_euler(alpha, beta, gamma):
328 '''
329 Uniquify eulerian angle triplet.
331 Put eulerian angle triplet into ranges compatible with
332 ``(dip, strike, -rake)`` conventions in seismology::
334 alpha (dip) : [0, pi/2]
335 beta (strike) : [0, 2*pi)
336 gamma (-rake) : [-pi, pi)
338 If ``alpha1`` is near to zero, ``beta`` is replaced by ``beta+gamma`` and
339 ``gamma`` is set to zero, to prevent this additional ambiguity.
341 If ``alpha`` is near to ``pi/2``, ``beta`` is put into the range
342 ``[0,pi)``.
343 '''
345 pi = math.pi
347 alpha = num.mod(alpha, 2.0*pi)
349 if 0.5*pi < alpha and alpha <= pi:
350 alpha = pi - alpha
351 beta = beta + pi
352 gamma = 2.0*pi - gamma
353 elif pi < alpha and alpha <= 1.5*pi:
354 alpha = alpha - pi
355 gamma = pi - gamma
356 elif 1.5*pi < alpha and alpha <= 2.0*pi:
357 alpha = 2.0*pi - alpha
358 beta = beta + pi
359 gamma = pi + gamma
361 alpha = num.mod(alpha, 2.0*pi)
362 beta = num.mod(beta, 2.0*pi)
363 gamma = num.mod(gamma+pi, 2.0*pi)-pi
365 # If dip is exactly 90 degrees, one is still
366 # free to choose between looking at the plane from either side.
367 # Choose to look at such that beta is in the range [0,180)
369 # This should prevent some problems, when dip is close to 90 degrees:
370 if abs(alpha - 0.5*pi) < 1e-10:
371 alpha = 0.5*pi
373 if abs(beta - pi) < 1e-10:
374 beta = pi
376 if abs(beta - 2.*pi) < 1e-10:
377 beta = 0.
379 if abs(beta) < 1e-10:
380 beta = 0.
382 if alpha == 0.5*pi and beta >= pi:
383 gamma = - gamma
384 beta = num.mod(beta-pi, 2.0*pi)
385 gamma = num.mod(gamma+pi, 2.0*pi)-pi
386 assert 0. <= beta < pi
387 assert -pi <= gamma < pi
389 if alpha < 1e-7:
390 beta = num.mod(beta + gamma, 2.0*pi)
391 gamma = 0.
393 return (alpha, beta, gamma)
396def cvec(x, y, z):
397 return num.array([[x, y, z]], dtype=float).T
400def rvec(x, y, z):
401 return num.array([[x, y, z]], dtype=float)
404def eigh_check(a):
405 evals, evecs = num.linalg.eigh(a)
406 assert evals[0] <= evals[1] <= evals[2]
407 return evals, evecs
410r2d = 180. / math.pi
411d2r = 1. / r2d
414def random_mt(x=None, scalar_moment=1.0, magnitude=None):
416 if magnitude is not None:
417 scalar_moment = magnitude_to_moment(magnitude)
419 if x is None:
420 x = num.random.random(6)
422 evals = x[:3] * 2. - 1.0
423 evals /= num.sqrt(num.sum(evals**2)) / math.sqrt(2.0)
424 rotmat = random_rotation(x[3:])
425 return scalar_moment * num.dot(
426 rotmat, num.dot(num.array(num.diag(evals)), rotmat.T))
429def random_m6(*args, **kwargs):
430 return to6(random_mt(*args, **kwargs))
433def random_dc(x=None, scalar_moment=1.0, magnitude=None):
434 if magnitude is not None:
435 scalar_moment = magnitude_to_moment(magnitude)
437 rotmat = random_rotation(x)
438 return scalar_moment * num.dot(
439 rotmat, num.dot(MomentTensor._m_unrot, rotmat.T))
442def sm(m):
443 return "/ %5.2F %5.2F %5.2F \\\n" % (m[0, 0], m[0, 1], m[0, 2]) + \
444 "| %5.2F %5.2F %5.2F |\n" % (m[1, 0], m[1, 1], m[1, 2]) + \
445 "\\ %5.2F %5.2F %5.2F /\n" % (m[2, 0], m[2, 1], m[2, 2])
448def as_mt(mt):
449 '''
450 Convenience function to convert various objects to moment tensor object.
452 Like :py:meth:``MomentTensor.from_values``, but does not create a new
453 :py:class:`MomentTensor` object when ``mt`` already is one.
454 '''
456 if isinstance(mt, MomentTensor):
457 return mt
458 else:
459 return MomentTensor.from_values(mt)
462class MomentTensor(Object):
464 '''
465 Moment tensor object
467 :param m: NumPy array in north-east-down convention
468 :param m_up_south_east: NumPy array in up-south-east convention
469 :param m_east_north_up: NumPy array in east-north-up convention
470 :param strike,dip,rake: fault plane angles in [degrees]
471 :param scalar_moment: scalar moment in [Nm]
472 :param magnitude: moment magnitude Mw
474 Global CMT catalog moment tensors use the up-south-east (USE) coordinate
475 system convention with :math:`r` (up), :math:`\\theta` (south), and
476 :math:`\\phi` (east).
478 .. math::
479 :nowrap:
481 \\begin{align*}
482 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
483 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
484 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
485 \\end{align*}
487 '''
489 mnn__ = Float.T(default=0.0)
490 mee__ = Float.T(default=0.0)
491 mdd__ = Float.T(default=0.0)
492 mne__ = Float.T(default=0.0)
493 mnd__ = Float.T(default=-1.0)
494 med__ = Float.T(default=0.0)
495 strike1__ = Float.T(default=None, optional=True) # read-only
496 dip1__ = Float.T(default=None, optional=True) # read-only
497 rake1__ = Float.T(default=None, optional=True) # read-only
498 strike2__ = Float.T(default=None, optional=True) # read-only
499 dip2__ = Float.T(default=None, optional=True) # read-only
500 rake2__ = Float.T(default=None, optional=True) # read-only
501 moment__ = Float.T(default=None, optional=True) # read-only
502 magnitude__ = Float.T(default=None, optional=True) # read-only
504 _flip_dc = num.array(
505 [[0., 0., -1.], [0., -1., 0.], [-1., 0., 0.]], dtype=float)
506 _to_up_south_east = num.array(
507 [[0., 0., -1.], [-1., 0., 0.], [0., 1., 0.]], dtype=float).T
508 _to_east_north_up = num.array(
509 [[0., 1., 0.], [1., 0., 0.], [0., 0., -1.]], dtype=float)
510 _m_unrot = num.array(
511 [[0., 0., -1.], [0., 0., 0.], [-1., 0., 0.]], dtype=float)
513 _u_evals, _u_evecs = eigh_check(_m_unrot)
515 @classmethod
516 def random_dc(cls, x=None, scalar_moment=1.0, magnitude=None):
517 '''
518 Create random oriented double-couple moment tensor
520 The rotations used are uniformly distributed in the space of rotations.
521 '''
522 return MomentTensor(
523 m=random_dc(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
525 @classmethod
526 def random_mt(cls, x=None, scalar_moment=1.0, magnitude=None):
527 '''
528 Create random moment tensor
530 Moment tensors produced by this function appear uniformly distributed
531 when shown in a Hudson's diagram. The rotations used are unifomly
532 distributed in the space of rotations.
533 '''
534 return MomentTensor(
535 m=random_mt(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
537 @classmethod
538 def from_values(cls, values):
539 '''
540 Alternative constructor for moment tensor objects
542 This constructor takes a :py:class:`MomentTensor` object, a tuple, list
543 or NumPy array with 3x3 or 3, 4, 6, or 7 elements to build a Moment
544 tensor object.
546 The ``values`` argument is interpreted depending on shape and type as
547 follows:
549 * ``(strike, dip, rake)``
550 * ``(strike, dip, rake, magnitude)``
551 * ``(mnn, mee, mdd, mne, mnd, med)``
552 * ``(mnn, mee, mdd, mne, mnd, med, magnitude)``
553 * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))``
554 * :py:class:`MomentTensor` object
555 '''
557 m = values_to_matrix(values)
558 return MomentTensor(m=m)
560 def __init__(
561 self, m=None, m_up_south_east=None, m_east_north_up=None,
562 strike=0., dip=0., rake=0., scalar_moment=1.,
563 mnn=None, mee=None, mdd=None, mne=None, mnd=None, med=None,
564 strike1=None, dip1=None, rake1=None,
565 strike2=None, dip2=None, rake2=None,
566 magnitude=None, moment=None):
568 Object.__init__(self, init_props=False)
570 if any(mxx is not None for mxx in (mnn, mee, mdd, mne, mnd, med)):
571 m = symmat6(mnn, mee, mdd, mne, mnd, med)
573 strike = d2r*strike
574 dip = d2r*dip
575 rake = d2r*rake
577 if isinstance(m, num.matrix):
578 m = m.A
580 if isinstance(m_up_south_east, num.matrix):
581 m_up_south_east = m_up_south_east.A
583 if isinstance(m_east_north_up, num.matrix):
584 m_east_north_up = m_east_north_up.A
586 if m_up_south_east is not None:
587 m = num.dot(
588 self._to_up_south_east,
589 num.dot(m_up_south_east, self._to_up_south_east.T))
591 if m_east_north_up is not None:
592 m = num.dot(
593 self._to_east_north_up,
594 num.dot(m_east_north_up, self._to_east_north_up.T))
596 if m is None:
597 if any(x is not None for x in (
598 strike1, dip1, rake1, strike2, dip2, rake2)):
599 raise Exception(
600 'strike1, dip1, rake1, strike2, dip2, rake2 are read-only '
601 'properties')
603 if moment is not None:
604 scalar_moment = moment
606 if magnitude is not None:
607 scalar_moment = magnitude_to_moment(magnitude)
609 rotmat1 = euler_to_matrix(dip, strike, -rake)
610 m = scalar_moment * num.dot(
611 rotmat1.T,
612 num.dot(MomentTensor._m_unrot, rotmat1))
614 self._m = m
615 self._update()
617 def _update(self):
618 m_evals, m_evecs = eigh_check(self._m)
619 if num.linalg.det(m_evecs) < 0.:
620 m_evecs *= -1.
622 rotmat1 = num.dot(m_evecs, MomentTensor._u_evecs.T).T
623 if num.linalg.det(rotmat1) < 0.:
624 rotmat1 *= -1.
626 self._m_eigenvals = m_evals
627 self._m_eigenvecs = m_evecs
629 self._rotmats = sorted(
630 [rotmat1, num.dot(MomentTensor._flip_dc, rotmat1)],
631 key=lambda m: num.abs(m.flat).tolist())
633 @property
634 def mnn(self):
635 return float(self._m[0, 0])
637 @mnn.setter
638 def mnn(self, value):
639 self._m[0, 0] = value
640 self._update()
642 @property
643 def mee(self):
644 return float(self._m[1, 1])
646 @mee.setter
647 def mee(self, value):
648 self._m[1, 1] = value
649 self._update()
651 @property
652 def mdd(self):
653 return float(self._m[2, 2])
655 @mdd.setter
656 def mdd(self, value):
657 self._m[2, 2] = value
658 self._update()
660 @property
661 def mne(self):
662 return float(self._m[0, 1])
664 @mne.setter
665 def mne(self, value):
666 self._m[0, 1] = value
667 self._m[1, 0] = value
668 self._update()
670 @property
671 def mnd(self):
672 return float(self._m[0, 2])
674 @mnd.setter
675 def mnd(self, value):
676 self._m[0, 2] = value
677 self._m[2, 0] = value
678 self._update()
680 @property
681 def med(self):
682 return float(self._m[1, 2])
684 @med.setter
685 def med(self, value):
686 self._m[1, 2] = value
687 self._m[2, 1] = value
688 self._update()
690 @property
691 def strike1(self):
692 return float(self.both_strike_dip_rake()[0][0])
694 @property
695 def dip1(self):
696 return float(self.both_strike_dip_rake()[0][1])
698 @property
699 def rake1(self):
700 return float(self.both_strike_dip_rake()[0][2])
702 @property
703 def strike2(self):
704 return float(self.both_strike_dip_rake()[1][0])
706 @property
707 def dip2(self):
708 return float(self.both_strike_dip_rake()[1][1])
710 @property
711 def rake2(self):
712 return float(self.both_strike_dip_rake()[1][2])
714 def both_strike_dip_rake(self):
715 '''
716 Get both possible (strike,dip,rake) triplets.
717 '''
718 results = []
719 for rotmat in self._rotmats:
720 alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)]
721 results.append((beta, alpha, -gamma))
723 return results
725 def p_axis(self):
726 '''
727 Get direction of p axis.
728 '''
729 return (self._m_eigenvecs.T)[0]
731 def t_axis(self):
732 '''
733 Get direction of t axis.
734 '''
735 return (self._m_eigenvecs.T)[2]
737 def null_axis(self):
738 '''
739 Get diretion of the null axis.
740 '''
741 return self._m_eigenvecs.T[1]
743 def eigenvals(self):
744 '''
745 Get the eigenvalues of the moment tensor in accending order.
747 :returns: ``(ep, en, et)``
748 '''
750 return self._m_eigenvals
752 def eigensystem(self):
753 '''
754 Get the eigenvalues and eigenvectors of the moment tensor.
756 :returns: ``(ep, en, et, vp, vn, vt)``
757 '''
759 vp = self.p_axis().flatten()
760 vn = self.null_axis().flatten()
761 vt = self.t_axis().flatten()
762 ep, en, et = self._m_eigenvals
763 return ep, en, et, vp, vn, vt
765 def both_slip_vectors(self):
766 '''
767 Get both possible slip directions.
768 '''
769 return [num.dot(rotmat, cvec(1., 0., 0.)) for rotmat in self._rotmats]
771 def m(self):
772 '''
773 Get plain moment tensor as 3x3 matrix.
774 '''
775 return self._m.copy()
777 def m6(self):
778 '''
779 Get the moment tensor as a six-element array.
781 :returns: ``(mnn, mee, mdd, mne, mnd, med)``
782 '''
783 return to6(self._m)
785 def m_up_south_east(self):
786 '''
787 Get moment tensor in up-south-east convention as 3x3 matrix.
789 .. math::
790 :nowrap:
792 \\begin{align*}
793 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
794 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
795 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
796 \\end{align*}
797 '''
799 return num.dot(
800 self._to_up_south_east.T,
801 num.dot(self._m, self._to_up_south_east))
803 def m_east_north_up(self):
804 '''
805 Get moment tensor in east-north-up convention as 3x3 matrix.
806 '''
808 return num.dot(
809 self._to_east_north_up.T, num.dot(self._m, self._to_east_north_up))
811 def m6_up_south_east(self):
812 '''
813 Get moment tensor in up-south-east convention as a six-element array.
815 :returns: ``(muu, mss, mee, mus, mue, mse)``
816 '''
817 return to6(self.m_up_south_east())
819 def m6_east_north_up(self):
820 '''
821 Get moment tensor in east-north-up convention as a six-element array.
823 :returns: ``(mee, mnn, muu, men, meu, mnu)``
824 '''
825 return to6(self.m_east_north_up())
827 def m_plain_double_couple(self):
828 '''
829 Get plain double couple with same scalar moment as moment tensor.
830 '''
831 rotmat1 = self._rotmats[0]
832 m = self.scalar_moment() * num.dot(
833 rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
834 return m
836 def moment_magnitude(self):
837 '''
838 Get moment magnitude of moment tensor.
839 '''
840 return moment_to_magnitude(self.scalar_moment())
842 def scalar_moment(self):
843 '''
844 Get the scalar moment of the moment tensor (Frobenius norm based)
846 .. math::
848 M_0 = \\frac{1}{\\sqrt{2}}\\sqrt{\\sum_{i,j} |M_{ij}|^2}
850 The scalar moment is calculated based on the Euclidean (Frobenius) norm
851 (Silver and Jordan, 1982). The scalar moment returned by this function
852 differs from the standard decomposition based definition of the scalar
853 moment for non-double-couple moment tensors.
854 '''
855 return num.linalg.norm(self._m_eigenvals) / math.sqrt(2.)
857 @property
858 def moment(self):
859 return float(self.scalar_moment())
861 @moment.setter
862 def moment(self, value):
863 self._m *= value / self.moment
864 self._update()
866 @property
867 def magnitude(self):
868 return float(self.moment_magnitude())
870 @magnitude.setter
871 def magnitude(self, value):
872 self._m *= magnitude_to_moment(value) / self.moment
873 self._update()
875 def __str__(self):
876 mexp = pow(10, math.ceil(num.log10(num.max(num.abs(self._m)))))
877 m = self._m / mexp
878 s = '''Scalar Moment [Nm]: M0 = %g (Mw = %3.1f)
879Moment Tensor [Nm]: Mnn = %6.3f, Mee = %6.3f, Mdd = %6.3f,
880 Mne = %6.3f, Mnd = %6.3f, Med = %6.3f [ x %g ]
881''' % (
882 self.scalar_moment(),
883 self.moment_magnitude(),
884 m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2],
885 mexp)
887 s += self.str_fault_planes()
888 return s
890 def str_fault_planes(self):
891 s = ''
892 for i, sdr in enumerate(self.both_strike_dip_rake()):
893 s += 'Fault plane %i [deg]: ' \
894 'strike = %3.0f, dip = %3.0f, slip-rake = %4.0f\n' \
895 % (i+1, sdr[0], sdr[1], sdr[2])
897 return s
899 def deviatoric(self):
900 '''
901 Get deviatoric part of moment tensor.
903 Returns a new moment tensor object with zero trace.
904 '''
906 m = self.m()
908 trace_m = num.trace(m)
909 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
910 m_devi = m - m_iso
911 mt = MomentTensor(m=m_devi)
912 return mt
914 def standard_decomposition(self):
915 '''
916 Decompose moment tensor into isotropic, DC and CLVD components.
918 Standard decomposition according to e.g. Jost and Herrmann 1989 is
919 returned as::
921 [
922 (moment_iso, ratio_iso, m_iso),
923 (moment_dc, ratio_dc, m_dc),
924 (moment_clvd, ratio_clvd, m_clvd),
925 (moment_devi, ratio_devi, m_devi),
926 (moment, 1.0, m)
927 ]
928 '''
930 epsilon = 1e-6
932 m = self.m()
934 trace_m = num.trace(m)
935 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
936 moment_iso = abs(trace_m / 3.)
938 m_devi = m - m_iso
940 evals, evecs = eigh_check(m_devi)
942 moment_devi = num.max(num.abs(evals))
943 moment = moment_iso + moment_devi
945 iorder = num.argsort(num.abs(evals))
946 evals_sorted = evals[iorder]
947 evecs_sorted = (evecs.T[iorder]).T
949 if moment_devi < epsilon * moment_iso:
950 signed_moment_dc = 0.
951 else:
952 assert -epsilon <= -evals_sorted[0] / evals_sorted[2] <= 0.5
953 signed_moment_dc = evals_sorted[2] * (1.0 + 2.0 * (
954 min(0.0, evals_sorted[0] / evals_sorted[2])))
956 moment_dc = abs(signed_moment_dc)
957 m_dc_es = num.dot(signed_moment_dc, num.diag([0., -1.0, 1.0]))
958 m_dc = num.dot(evecs_sorted, num.dot(m_dc_es, evecs_sorted.T))
960 m_clvd = m_devi - m_dc
962 moment_clvd = moment_devi - moment_dc
964 ratio_dc = moment_dc / moment
965 ratio_clvd = moment_clvd / moment
966 ratio_iso = moment_iso / moment
967 ratio_devi = moment_devi / moment
969 return [
970 (moment_iso, ratio_iso, m_iso),
971 (moment_dc, ratio_dc, m_dc),
972 (moment_clvd, ratio_clvd, m_clvd),
973 (moment_devi, ratio_devi, m_devi),
974 (moment, 1.0, m)]
976 def rotated(self, rot):
977 '''
978 Get rotated moment tensor.
980 :param rot: ratation matrix, coordinate system is NED
981 :returns: new :py:class:`MomentTensor` object
982 '''
984 rotmat = num.array(rot)
985 return MomentTensor(m=num.dot(rotmat, num.dot(self.m(), rotmat.T)))
987 def random_rotated(self, angle_std=None, angle=None, rstate=None):
988 '''
989 Get distorted MT by rotation around random axis and angle.
991 :param angle_std: angles are drawn from a normal distribution with
992 zero mean and given standard deviation [degrees]
993 :param angle: set angle [degrees], only axis will be random
994 :param rstate: :py:class:`numpy.random.RandomState` object, can be
995 used to create reproducible pseudo-random sequences
996 :returns: new :py:class:`MomentTensor` object
997 '''
999 assert (angle_std is None) != (angle is None), \
1000 'either angle or angle_std must be given'
1002 if angle_std is not None:
1003 rstate = rstate or num.random
1004 angle = rstate.normal(scale=angle_std)
1006 axis = random_axis(rstate=rstate)
1007 rot = rotation_from_angle_and_axis(angle, axis)
1008 return self.rotated(rot)
1011def other_plane(strike, dip, rake):
1012 '''
1013 Get the respectively other plane in the double-couple ambiguity.
1014 '''
1016 mt = MomentTensor(strike=strike, dip=dip, rake=rake)
1017 both_sdr = mt.both_strike_dip_rake()
1018 w = [sum([abs(x-y) for x, y in zip(both_sdr[i], (strike, dip, rake))])
1019 for i in (0, 1)]
1021 if w[0] < w[1]:
1022 return both_sdr[1]
1023 else:
1024 return both_sdr[0]
1027def dsdr(sdr1, sdr2):
1028 s1, d1, r1 = sdr1
1029 s2, d2, r2 = sdr2
1031 s1 = s1 % 360.
1032 s2 = s2 % 360.
1033 r1 = r1 % 360.
1034 r2 = r2 % 360.
1036 ds = abs(s1 - s2)
1037 ds = ds if ds <= 180. else 360. - ds
1039 dr = abs(r1 - r2)
1040 dr = dr if dr <= 180. else 360. - dr
1042 dd = abs(d1 - d2)
1044 return math.sqrt(ds**2 + dr**2 + dd**2)
1047def order_like(sdrs, sdrs_ref):
1048 '''
1049 Order strike-dip-rake pair post closely to a given reference pair.
1051 :param sdrs: tuple, ``((strike1, dip1, rake1), (strike2, dip2, rake2))``
1052 :param sdrs_ref: as above but with reference pair
1053 '''
1055 d1 = min(dsdr(sdrs[0], sdrs_ref[0]), dsdr(sdrs[1], sdrs_ref[1]))
1056 d2 = min(dsdr(sdrs[0], sdrs_ref[1]), dsdr(sdrs[1], sdrs_ref[0]))
1057 if d1 < d2:
1058 return sdrs
1059 else:
1060 return sdrs[::-1]
1063def _tpb2q(t, p, b):
1064 eps = 0.001
1065 tqw = 1. + t[0] + p[1] + b[2]
1066 tqx = 1. + t[0] - p[1] - b[2]
1067 tqy = 1. - t[0] + p[1] - b[2]
1068 tqz = 1. - t[0] - p[1] + b[2]
1070 q = num.zeros(4)
1071 if tqw > eps:
1072 q[0] = 0.5 * math.sqrt(tqw)
1073 q[1:] = p[2] - b[1], b[0] - t[2], t[1] - p[0]
1074 elif tqx > eps:
1075 q[0] = 0.5 * math.sqrt(tqx)
1076 q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1077 elif tqy > eps:
1078 q[0] = 0.5 * math.sqrt(tqy)
1079 q[1:] = b[0] - t[2], p[0] + t[1], b[1] + p[2]
1080 elif tqz > eps:
1081 q[0] = 0.5 * math.sqrt(tqz)
1082 q[1:] = t[1] - p[0], b[0] + t[2], b[1] + p[2]
1083 else:
1084 raise Exception('should not reach this line, check theory!')
1085 # q[0] = max(0.5 * math.sqrt(tqx), eps)
1086 # q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1088 q[1:] /= 4.0 * q[0]
1090 q /= math.sqrt(num.sum(q**2))
1092 return q
1095_pbt2tpb = num.array(((0., 0., 1.), (1., 0., 0.), (0., 1., 0.)))
1098def kagan_angle(mt1, mt2):
1099 '''
1100 Given two moment tensors, return the Kagan angle in degrees.
1102 After Kagan (1991) and Tape & Tape (2012).
1103 '''
1105 ai = num.dot(_pbt2tpb, mt1._m_eigenvecs.T)
1106 aj = num.dot(_pbt2tpb, mt2._m_eigenvecs.T)
1108 u = num.dot(ai, aj.T)
1110 tk, pk, bk = u
1112 qk = _tpb2q(tk, pk, bk)
1114 return 2. * r2d * math.acos(num.max(num.abs(qk)))
1117def rand_to_gutenberg_richter(rand, b_value, magnitude_min):
1118 '''
1119 Draw magnitude from Gutenberg Richter distribution.
1120 '''
1121 return magnitude_min + num.log10(1.-rand) / -b_value
1124def magnitude_to_duration_gcmt(magnitudes):
1125 '''
1126 Scaling relation used by Global CMT catalog for most of its events.
1127 '''
1129 mom = magnitude_to_moment(magnitudes)
1130 return (mom / 1.1e16)**(1./3.)
1133if __name__ == '__main__':
1135 import sys
1136 v = list(map(float, sys.argv[1:]))
1137 mt = MomentTensor.from_values(v)
1138 print(mt)