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 [
770 num.dot(rotmat.T, cvec(1., 0., 0.)) for rotmat in self._rotmats]
772 def m(self):
773 '''
774 Get plain moment tensor as 3x3 matrix.
775 '''
776 return self._m.copy()
778 def m6(self):
779 '''
780 Get the moment tensor as a six-element array.
782 :returns: ``(mnn, mee, mdd, mne, mnd, med)``
783 '''
784 return to6(self._m)
786 def m_up_south_east(self):
787 '''
788 Get moment tensor in up-south-east convention as 3x3 matrix.
790 .. math::
791 :nowrap:
793 \\begin{align*}
794 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
795 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
796 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
797 \\end{align*}
798 '''
800 return num.dot(
801 self._to_up_south_east.T,
802 num.dot(self._m, self._to_up_south_east))
804 def m_east_north_up(self):
805 '''
806 Get moment tensor in east-north-up convention as 3x3 matrix.
807 '''
809 return num.dot(
810 self._to_east_north_up.T, num.dot(self._m, self._to_east_north_up))
812 def m6_up_south_east(self):
813 '''
814 Get moment tensor in up-south-east convention as a six-element array.
816 :returns: ``(muu, mss, mee, mus, mue, mse)``
817 '''
818 return to6(self.m_up_south_east())
820 def m6_east_north_up(self):
821 '''
822 Get moment tensor in east-north-up convention as a six-element array.
824 :returns: ``(mee, mnn, muu, men, meu, mnu)``
825 '''
826 return to6(self.m_east_north_up())
828 def m_plain_double_couple(self):
829 '''
830 Get plain double couple with same scalar moment as moment tensor.
831 '''
832 rotmat1 = self._rotmats[0]
833 m = self.scalar_moment() * num.dot(
834 rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
835 return m
837 def moment_magnitude(self):
838 '''
839 Get moment magnitude of moment tensor.
840 '''
841 return moment_to_magnitude(self.scalar_moment())
843 def scalar_moment(self):
844 '''
845 Get the scalar moment of the moment tensor (Frobenius norm based)
847 .. math::
849 M_0 = \\frac{1}{\\sqrt{2}}\\sqrt{\\sum_{i,j} |M_{ij}|^2}
851 The scalar moment is calculated based on the Euclidean (Frobenius) norm
852 (Silver and Jordan, 1982). The scalar moment returned by this function
853 differs from the standard decomposition based definition of the scalar
854 moment for non-double-couple moment tensors.
855 '''
856 return num.linalg.norm(self._m_eigenvals) / math.sqrt(2.)
858 @property
859 def moment(self):
860 return float(self.scalar_moment())
862 @moment.setter
863 def moment(self, value):
864 self._m *= value / self.moment
865 self._update()
867 @property
868 def magnitude(self):
869 return float(self.moment_magnitude())
871 @magnitude.setter
872 def magnitude(self, value):
873 self._m *= magnitude_to_moment(value) / self.moment
874 self._update()
876 def __str__(self):
877 mexp = pow(10, math.ceil(num.log10(num.max(num.abs(self._m)))))
878 m = self._m / mexp
879 s = '''Scalar Moment [Nm]: M0 = %g (Mw = %3.1f)
880Moment Tensor [Nm]: Mnn = %6.3f, Mee = %6.3f, Mdd = %6.3f,
881 Mne = %6.3f, Mnd = %6.3f, Med = %6.3f [ x %g ]
882''' % (
883 self.scalar_moment(),
884 self.moment_magnitude(),
885 m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2],
886 mexp)
888 s += self.str_fault_planes()
889 return s
891 def str_fault_planes(self):
892 s = ''
893 for i, sdr in enumerate(self.both_strike_dip_rake()):
894 s += 'Fault plane %i [deg]: ' \
895 'strike = %3.0f, dip = %3.0f, slip-rake = %4.0f\n' \
896 % (i+1, sdr[0], sdr[1], sdr[2])
898 return s
900 def deviatoric(self):
901 '''
902 Get deviatoric part of moment tensor.
904 Returns a new moment tensor object with zero trace.
905 '''
907 m = self.m()
909 trace_m = num.trace(m)
910 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
911 m_devi = m - m_iso
912 mt = MomentTensor(m=m_devi)
913 return mt
915 def standard_decomposition(self):
916 '''
917 Decompose moment tensor into isotropic, DC and CLVD components.
919 Standard decomposition according to e.g. Jost and Herrmann 1989 is
920 returned as::
922 [
923 (moment_iso, ratio_iso, m_iso),
924 (moment_dc, ratio_dc, m_dc),
925 (moment_clvd, ratio_clvd, m_clvd),
926 (moment_devi, ratio_devi, m_devi),
927 (moment, 1.0, m)
928 ]
929 '''
931 epsilon = 1e-6
933 m = self.m()
935 trace_m = num.trace(m)
936 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
937 moment_iso = abs(trace_m / 3.)
939 m_devi = m - m_iso
941 evals, evecs = eigh_check(m_devi)
943 moment_devi = num.max(num.abs(evals))
944 moment = moment_iso + moment_devi
946 iorder = num.argsort(num.abs(evals))
947 evals_sorted = evals[iorder]
948 evecs_sorted = (evecs.T[iorder]).T
950 if moment_devi < epsilon * moment_iso:
951 signed_moment_dc = 0.
952 else:
953 assert -epsilon <= -evals_sorted[0] / evals_sorted[2] <= 0.5
954 signed_moment_dc = evals_sorted[2] * (1.0 + 2.0 * (
955 min(0.0, evals_sorted[0] / evals_sorted[2])))
957 moment_dc = abs(signed_moment_dc)
958 m_dc_es = num.dot(signed_moment_dc, num.diag([0., -1.0, 1.0]))
959 m_dc = num.dot(evecs_sorted, num.dot(m_dc_es, evecs_sorted.T))
961 m_clvd = m_devi - m_dc
963 moment_clvd = moment_devi - moment_dc
965 ratio_dc = moment_dc / moment
966 ratio_clvd = moment_clvd / moment
967 ratio_iso = moment_iso / moment
968 ratio_devi = moment_devi / moment
970 return [
971 (moment_iso, ratio_iso, m_iso),
972 (moment_dc, ratio_dc, m_dc),
973 (moment_clvd, ratio_clvd, m_clvd),
974 (moment_devi, ratio_devi, m_devi),
975 (moment, 1.0, m)]
977 def rotated(self, rot):
978 '''
979 Get rotated moment tensor.
981 :param rot: ratation matrix, coordinate system is NED
982 :returns: new :py:class:`MomentTensor` object
983 '''
985 rotmat = num.array(rot)
986 return MomentTensor(m=num.dot(rotmat, num.dot(self.m(), rotmat.T)))
988 def random_rotated(self, angle_std=None, angle=None, rstate=None):
989 '''
990 Get distorted MT by rotation around random axis and angle.
992 :param angle_std: angles are drawn from a normal distribution with
993 zero mean and given standard deviation [degrees]
994 :param angle: set angle [degrees], only axis will be random
995 :param rstate: :py:class:`numpy.random.RandomState` object, can be
996 used to create reproducible pseudo-random sequences
997 :returns: new :py:class:`MomentTensor` object
998 '''
1000 assert (angle_std is None) != (angle is None), \
1001 'either angle or angle_std must be given'
1003 if angle_std is not None:
1004 rstate = rstate or num.random
1005 angle = rstate.normal(scale=angle_std)
1007 axis = random_axis(rstate=rstate)
1008 rot = rotation_from_angle_and_axis(angle, axis)
1009 return self.rotated(rot)
1012def other_plane(strike, dip, rake):
1013 '''
1014 Get the respectively other plane in the double-couple ambiguity.
1015 '''
1017 mt = MomentTensor(strike=strike, dip=dip, rake=rake)
1018 both_sdr = mt.both_strike_dip_rake()
1019 w = [sum([abs(x-y) for x, y in zip(both_sdr[i], (strike, dip, rake))])
1020 for i in (0, 1)]
1022 if w[0] < w[1]:
1023 return both_sdr[1]
1024 else:
1025 return both_sdr[0]
1028def dsdr(sdr1, sdr2):
1029 s1, d1, r1 = sdr1
1030 s2, d2, r2 = sdr2
1032 s1 = s1 % 360.
1033 s2 = s2 % 360.
1034 r1 = r1 % 360.
1035 r2 = r2 % 360.
1037 ds = abs(s1 - s2)
1038 ds = ds if ds <= 180. else 360. - ds
1040 dr = abs(r1 - r2)
1041 dr = dr if dr <= 180. else 360. - dr
1043 dd = abs(d1 - d2)
1045 return math.sqrt(ds**2 + dr**2 + dd**2)
1048def order_like(sdrs, sdrs_ref):
1049 '''
1050 Order strike-dip-rake pair post closely to a given reference pair.
1052 :param sdrs: tuple, ``((strike1, dip1, rake1), (strike2, dip2, rake2))``
1053 :param sdrs_ref: as above but with reference pair
1054 '''
1056 d1 = min(dsdr(sdrs[0], sdrs_ref[0]), dsdr(sdrs[1], sdrs_ref[1]))
1057 d2 = min(dsdr(sdrs[0], sdrs_ref[1]), dsdr(sdrs[1], sdrs_ref[0]))
1058 if d1 < d2:
1059 return sdrs
1060 else:
1061 return sdrs[::-1]
1064def _tpb2q(t, p, b):
1065 eps = 0.001
1066 tqw = 1. + t[0] + p[1] + b[2]
1067 tqx = 1. + t[0] - p[1] - b[2]
1068 tqy = 1. - t[0] + p[1] - b[2]
1069 tqz = 1. - t[0] - p[1] + b[2]
1071 q = num.zeros(4)
1072 if tqw > eps:
1073 q[0] = 0.5 * math.sqrt(tqw)
1074 q[1:] = p[2] - b[1], b[0] - t[2], t[1] - p[0]
1075 elif tqx > eps:
1076 q[0] = 0.5 * math.sqrt(tqx)
1077 q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1078 elif tqy > eps:
1079 q[0] = 0.5 * math.sqrt(tqy)
1080 q[1:] = b[0] - t[2], p[0] + t[1], b[1] + p[2]
1081 elif tqz > eps:
1082 q[0] = 0.5 * math.sqrt(tqz)
1083 q[1:] = t[1] - p[0], b[0] + t[2], b[1] + p[2]
1084 else:
1085 raise Exception('should not reach this line, check theory!')
1086 # q[0] = max(0.5 * math.sqrt(tqx), eps)
1087 # q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1089 q[1:] /= 4.0 * q[0]
1091 q /= math.sqrt(num.sum(q**2))
1093 return q
1096_pbt2tpb = num.array(((0., 0., 1.), (1., 0., 0.), (0., 1., 0.)))
1099def kagan_angle(mt1, mt2):
1100 '''
1101 Given two moment tensors, return the Kagan angle in degrees.
1103 After Kagan (1991) and Tape & Tape (2012).
1104 '''
1106 ai = num.dot(_pbt2tpb, mt1._m_eigenvecs.T)
1107 aj = num.dot(_pbt2tpb, mt2._m_eigenvecs.T)
1109 u = num.dot(ai, aj.T)
1111 tk, pk, bk = u
1113 qk = _tpb2q(tk, pk, bk)
1115 return 2. * r2d * math.acos(num.max(num.abs(qk)))
1118def rand_to_gutenberg_richter(rand, b_value, magnitude_min):
1119 '''
1120 Draw magnitude from Gutenberg Richter distribution.
1121 '''
1122 return magnitude_min + num.log10(1.-rand) / -b_value
1125def magnitude_to_duration_gcmt(magnitudes):
1126 '''
1127 Scaling relation used by Global CMT catalog for most of its events.
1128 '''
1130 mom = magnitude_to_moment(magnitudes)
1131 return (mom / 1.1e16)**(1./3.)
1134if __name__ == '__main__':
1136 import sys
1137 v = list(map(float, sys.argv[1:]))
1138 mt = MomentTensor.from_values(v)
1139 print(mt)