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'''
24import math
25import numpy as num
27from .guts import Object, Float
29guts_prefix = 'pf'
31dynecm = 1e-7
34def random_axis(rstate=None):
35 '''
36 Get randomly oriented unit vector.
38 :param rstate: :py:class:`numpy.random.RandomState` object, can be used to
39 create reproducible pseudo-random sequences
40 '''
41 rstate = rstate or num.random
42 while True:
43 axis = rstate.uniform(size=3) * 2.0 - 1.0
44 uabs = math.sqrt(num.sum(axis**2))
45 if 0.001 < uabs < 1.0:
46 return axis / uabs
49def rotation_from_angle_and_axis(angle, axis):
50 '''
51 Build rotation matrix based on axis and angle.
53 :param angle: rotation angle [degrees]
54 :param axis: orientation of rotation axis, either in spherical
55 coordinates ``(theta, phi)`` [degrees], or as a unit vector
56 ``(ux, uy, uz)``.
57 '''
59 if len(axis) == 2:
60 theta, phi = axis
61 ux = math.sin(d2r*theta)*math.cos(d2r*phi)
62 uy = math.sin(d2r*theta)*math.sin(d2r*phi)
63 uz = math.cos(d2r*theta)
65 elif len(axis) == 3:
66 axis = num.asarray(axis)
67 uabs = math.sqrt(num.sum(axis**2))
68 ux, uy, uz = axis / uabs
69 else:
70 assert False
72 ct = math.cos(d2r*angle)
73 st = math.sin(d2r*angle)
74 return num.array([
75 [ct + ux**2*(1.-ct), ux*uy*(1.-ct)-uz*st, ux*uz*(1.-ct)+uy*st],
76 [uy*ux*(1.-ct)+uz*st, ct+uy**2*(1.-ct), uy*uz*(1.-ct)-ux*st],
77 [uz*ux*(1.-ct)-uy*st, uz*uy*(1.-ct)+ux*st, ct+uz**2*(1.-ct)]
78 ])
81def random_rotation(x=None):
82 '''
83 Get random rotation matrix.
85 A random rotation matrix, drawn from a uniform distrubution in the space
86 of rotations is returned, after Avro 1992 - "Fast random rotation
87 matrices".
89 :param x: three (uniform random) numbers in the range [0, 1[ used as input
90 to the distribution tranformation. If ``None``, random numbers are
91 used. Can be used to create grids of random rotations with uniform
92 density in rotation space.
93 '''
95 if x is not None:
96 x1, x2, x3 = x
97 else:
98 x1, x2, x3 = num.random.random(3)
100 phi = math.pi*2.0*x1
102 zrot = num.array([
103 [math.cos(phi), math.sin(phi), 0.],
104 [-math.sin(phi), math.cos(phi), 0.],
105 [0., 0., 1.]])
107 lam = math.pi*2.0*x2
109 v = num.array([[
110 math.cos(lam)*math.sqrt(x3),
111 math.sin(lam)*math.sqrt(x3),
112 math.sqrt(1.-x3)]]).T
114 house = num.identity(3) - 2.0 * num.dot(v, v.T)
115 return -num.dot(house, zrot)
118def rand(mi, ma):
119 return float(num.random.uniform(mi, ma))
122def randdip(mi, ma):
123 mi_ = 0.5*(math.cos(mi * math.pi/180.)+1.)
124 ma_ = 0.5*(math.cos(ma * math.pi/180.)+1.)
125 return math.acos(rand(mi_, ma_)*2.-1.)*180./math.pi
128def random_strike_dip_rake(
129 strikemin=0., strikemax=360.,
130 dipmin=0., dipmax=90.,
131 rakemin=-180., rakemax=180.):
133 '''
134 Get random strike, dip, rake triplet.
136 .. note::
138 Might not produce a homogeneous distribution of mechanisms. Better use
139 :py:meth:`MomentTensor.random_dc` which is based on
140 :py:func:`random_rotation`.
141 '''
143 strike = rand(strikemin, strikemax)
144 dip = randdip(dipmin, dipmax)
145 rake = rand(rakemin, rakemax)
147 return strike, dip, rake
150def to6(m):
151 '''
152 Get non-redundant components from symmetric 3x3 matrix.
154 :returns: 1D NumPy array with entries ordered like
155 ``(a_xx, a_yy, a_zz, a_xy, a_xz, a_yz)``
156 '''
158 return num.array([m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2]])
161def symmat6(a_xx, a_yy, a_zz, a_xy, a_xz, a_yz):
162 '''
163 Create symmetric 3x3 matrix from its 6 non-redundant values.
164 '''
166 return num.array([[a_xx, a_xy, a_xz],
167 [a_xy, a_yy, a_yz],
168 [a_xz, a_yz, a_zz]], dtype=float)
171def values_to_matrix(values):
172 '''
173 Convert anything to moment tensor represented as a NumPy array.
175 Transforms :py:class:`MomentTensor` objects, tuples, lists and NumPy arrays
176 with 3x3 or 3, 4, 6, or 7 elements into NumPy 3x3 array objects.
178 The ``values`` argument is interpreted depending on shape and type as
179 follows:
181 * ``(strike, dip, rake)``
182 * ``(strike, dip, rake, magnitude)``
183 * ``(mnn, mee, mdd, mne, mnd, med)``
184 * ``(mnn, mee, mdd, mne, mnd, med, magnitude)``
185 * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))``
186 * :py:class:`MomentTensor`
187 '''
189 if isinstance(values, (tuple, list)):
190 values = num.asarray(values, dtype=float)
192 if isinstance(values, MomentTensor):
193 return values.m()
195 elif isinstance(values, num.ndarray):
196 if values.shape == (3,):
197 strike, dip, rake = values
198 rotmat1 = euler_to_matrix(d2r*dip, d2r*strike, -d2r*rake)
199 return num.dot(rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
201 elif values.shape == (4,):
202 strike, dip, rake, magnitude = values
203 moment = magnitude_to_moment(magnitude)
204 rotmat1 = euler_to_matrix(d2r*dip, d2r*strike, -d2r*rake)
205 return moment * num.dot(
206 rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
208 elif values.shape == (6,):
209 return symmat6(*values)
211 elif values.shape == (7,):
212 magnitude = values[6]
213 moment = magnitude_to_moment(magnitude)
214 mt = symmat6(*values[:6])
215 mt *= moment / (num.linalg.norm(mt) / math.sqrt(2.0))
216 return mt
218 elif values.shape == (3, 3):
219 return num.asarray(values, dtype=float)
221 raise Exception('cannot convert object to 3x3 matrix')
224def moment_to_magnitude(moment):
225 '''
226 Convert scalar moment to moment magnitude Mw.
228 :param moment: scalar moment [Nm]
229 :returns: moment magnitude Mw
231 Moment magnitude is defined as
233 .. math::
235 M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0 \\cdot 10^7) - 10.7
237 where :math:`M_0` is the scalar moment given in [Nm].
239 .. note::
241 Global CMT uses 10.7333333 instead of 10.7, based on [Kanamori 1977],
242 10.7 is from [Hanks and Kanamori 1979].
243 '''
245 return num.log10(moment*1.0e7) / 1.5 - 10.7
248def magnitude_to_moment(magnitude):
249 '''
250 Convert moment magnitude Mw to scalar moment.
252 :param magnitude: moment magnitude
253 :returns: scalar moment [Nm]
255 See :py:func:`moment_to_magnitude`.
256 '''
258 return 10.0**(1.5*(magnitude+10.7))*1.0e-7
261magnitude_1Nm = moment_to_magnitude(1.0)
264def euler_to_matrix(alpha, beta, gamma):
265 '''
266 Given euler angle triplet, create rotation matrix
268 Given coordinate system `(x,y,z)` and rotated system `(xs,ys,zs)`
269 the line of nodes is the intersection between the `x,y` and the `xs,ys`
270 planes.
272 :param alpha: is the angle between the `z`-axis and the `zs`-axis [rad]
273 :param beta: is the angle between the `x`-axis and the line of nodes [rad]
274 :param gamma: is the angle between the line of nodes and the `xs`-axis
275 [rad]
277 Usage for moment tensors::
279 m_unrot = numpy.array([[0,0,-1],[0,0,0],[-1,0,0]])
280 euler_to_matrix(dip,strike,-rake, rotmat)
281 m = num.dot(rotmat.T, num.dot(m_unrot, rotmat))
283 '''
285 ca = math.cos(alpha)
286 cb = math.cos(beta)
287 cg = math.cos(gamma)
288 sa = math.sin(alpha)
289 sb = math.sin(beta)
290 sg = math.sin(gamma)
292 mat = num.array([[cb*cg-ca*sb*sg, sb*cg+ca*cb*sg, sa*sg],
293 [-cb*sg-ca*sb*cg, -sb*sg+ca*cb*cg, sa*cg],
294 [sa*sb, -sa*cb, ca]], dtype=float)
295 return mat
298def matrix_to_euler(rotmat):
299 '''
300 Get eulerian angle triplet from rotation matrix.
301 '''
303 ex = cvec(1., 0., 0.)
304 ez = cvec(0., 0., 1.)
305 exs = num.dot(rotmat.T, ex)
306 ezs = num.dot(rotmat.T, ez)
307 enodes = num.cross(ez.T, ezs.T).T
308 if num.linalg.norm(enodes) < 1e-10:
309 enodes = exs
310 enodess = num.dot(rotmat, enodes)
311 cos_alpha = float((num.dot(ez.T, ezs)))
312 if cos_alpha > 1.:
313 cos_alpha = 1.
315 if cos_alpha < -1.:
316 cos_alpha = -1.
318 alpha = math.acos(cos_alpha)
319 beta = num.mod(math.atan2(enodes[1, 0], enodes[0, 0]), math.pi*2.)
320 gamma = num.mod(-math.atan2(enodess[1, 0], enodess[0, 0]), math.pi*2.)
322 return unique_euler(alpha, beta, gamma)
325def unique_euler(alpha, beta, gamma):
326 '''
327 Uniquify eulerian angle triplet.
329 Put eulerian angle triplet into ranges compatible with
330 ``(dip, strike, -rake)`` conventions in seismology::
332 alpha (dip) : [0, pi/2]
333 beta (strike) : [0, 2*pi)
334 gamma (-rake) : [-pi, pi)
336 If ``alpha1`` is near to zero, ``beta`` is replaced by ``beta+gamma`` and
337 ``gamma`` is set to zero, to prevent this additional ambiguity.
339 If ``alpha`` is near to ``pi/2``, ``beta`` is put into the range
340 ``[0,pi)``.
341 '''
343 pi = math.pi
345 alpha = num.mod(alpha, 2.0*pi)
347 if 0.5*pi < alpha and alpha <= pi:
348 alpha = pi - alpha
349 beta = beta + pi
350 gamma = 2.0*pi - gamma
351 elif pi < alpha and alpha <= 1.5*pi:
352 alpha = alpha - pi
353 gamma = pi - gamma
354 elif 1.5*pi < alpha and alpha <= 2.0*pi:
355 alpha = 2.0*pi - alpha
356 beta = beta + pi
357 gamma = pi + gamma
359 alpha = num.mod(alpha, 2.0*pi)
360 beta = num.mod(beta, 2.0*pi)
361 gamma = num.mod(gamma+pi, 2.0*pi)-pi
363 # If dip is exactly 90 degrees, one is still
364 # free to choose between looking at the plane from either side.
365 # Choose to look at such that beta is in the range [0,180)
367 # This should prevent some problems, when dip is close to 90 degrees:
368 if abs(alpha - 0.5*pi) < 1e-10:
369 alpha = 0.5*pi
371 if abs(beta - pi) < 1e-10:
372 beta = pi
374 if abs(beta - 2.*pi) < 1e-10:
375 beta = 0.
377 if abs(beta) < 1e-10:
378 beta = 0.
380 if alpha == 0.5*pi and beta >= pi:
381 gamma = - gamma
382 beta = num.mod(beta-pi, 2.0*pi)
383 gamma = num.mod(gamma+pi, 2.0*pi)-pi
384 assert 0. <= beta < pi
385 assert -pi <= gamma < pi
387 if alpha < 1e-7:
388 beta = num.mod(beta + gamma, 2.0*pi)
389 gamma = 0.
391 return (alpha, beta, gamma)
394def cvec(x, y, z):
395 return num.array([[x, y, z]], dtype=float).T
398def rvec(x, y, z):
399 return num.array([[x, y, z]], dtype=float)
402def eigh_check(a):
403 evals, evecs = num.linalg.eigh(a)
404 assert evals[0] <= evals[1] <= evals[2]
405 return evals, evecs
408r2d = 180. / math.pi
409d2r = 1. / r2d
412def random_mt(x=None, scalar_moment=1.0, magnitude=None):
414 if magnitude is not None:
415 scalar_moment = magnitude_to_moment(magnitude)
417 if x is None:
418 x = num.random.random(6)
420 evals = x[:3] * 2. - 1.0
421 evals /= num.sqrt(num.sum(evals**2)) / math.sqrt(2.0)
422 rotmat = random_rotation(x[3:])
423 return scalar_moment * num.dot(
424 rotmat, num.dot(num.array(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 * num.dot(
437 rotmat, num.dot(MomentTensor._m_unrot, rotmat.T))
440def sm(m):
441 return '/ %5.2F %5.2F %5.2F \\\n' % (m[0, 0], m[0, 1], m[0, 2]) + \
442 '| %5.2F %5.2F %5.2F |\n' % (m[1, 0], m[1, 1], m[1, 2]) + \
443 '\\ %5.2F %5.2F %5.2F /\n' % (m[2, 0], m[2, 1], m[2, 2])
446def as_mt(mt):
447 '''
448 Convenience function to convert various objects to moment tensor object.
450 Like :py:meth:`MomentTensor.from_values`, but does not create a new
451 :py:class:`MomentTensor` object when ``mt`` already is one.
452 '''
454 if isinstance(mt, MomentTensor):
455 return mt
456 else:
457 return MomentTensor.from_values(mt)
460def pt_axes_to_strike_dip_rake(p_axis, t_axis):
461 n_axis = num.cross(p_axis, t_axis)
462 m_evecs = num.vstack([p_axis, n_axis, t_axis]).T
463 if num.linalg.det(m_evecs) < 0.:
464 m_evecs *= -1.
465 rotmat = num.dot(m_evecs, MomentTensor._u_evecs.T).T
466 if num.linalg.det(rotmat) < 0.:
467 rotmat *= -1.
468 alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)]
469 return (beta, alpha, -gamma)
472class MomentTensor(Object):
474 '''
475 Moment tensor object
477 :param m: NumPy array in north-east-down convention
478 :param m_up_south_east: NumPy array in up-south-east convention
479 :param m_east_north_up: NumPy array in east-north-up convention
480 :param strike,dip,rake: fault plane angles in [degrees]
481 :param p_axis,t_axis: initialize double-couple from p and t axes
482 :param scalar_moment: scalar moment in [Nm]
483 :param magnitude: moment magnitude Mw
485 Global CMT catalog moment tensors use the up-south-east (USE) coordinate
486 system convention with :math:`r` (up), :math:`\\theta` (south), and
487 :math:`\\phi` (east).
489 .. math::
490 :nowrap:
492 \\begin{align*}
493 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
494 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
495 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
496 \\end{align*}
498 '''
500 mnn__ = Float.T(default=0.0)
501 mee__ = Float.T(default=0.0)
502 mdd__ = Float.T(default=0.0)
503 mne__ = Float.T(default=0.0)
504 mnd__ = Float.T(default=-1.0)
505 med__ = Float.T(default=0.0)
506 strike1__ = Float.T(default=None, optional=True) # read-only
507 dip1__ = Float.T(default=None, optional=True) # read-only
508 rake1__ = Float.T(default=None, optional=True) # read-only
509 strike2__ = Float.T(default=None, optional=True) # read-only
510 dip2__ = Float.T(default=None, optional=True) # read-only
511 rake2__ = Float.T(default=None, optional=True) # read-only
512 moment__ = Float.T(default=None, optional=True) # read-only
513 magnitude__ = Float.T(default=None, optional=True) # read-only
515 _flip_dc = num.array(
516 [[0., 0., -1.], [0., -1., 0.], [-1., 0., 0.]], dtype=float)
517 _to_up_south_east = num.array(
518 [[0., 0., -1.], [-1., 0., 0.], [0., 1., 0.]], dtype=float).T
519 _to_east_north_up = num.array(
520 [[0., 1., 0.], [1., 0., 0.], [0., 0., -1.]], dtype=float)
521 _m_unrot = num.array(
522 [[0., 0., -1.], [0., 0., 0.], [-1., 0., 0.]], dtype=float)
524 _u_evals, _u_evecs = eigh_check(_m_unrot)
526 @classmethod
527 def random_dc(cls, x=None, scalar_moment=1.0, magnitude=None):
528 '''
529 Create random oriented double-couple moment tensor
531 The rotations used are uniformly distributed in the space of rotations.
532 '''
533 return MomentTensor(
534 m=random_dc(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
536 @classmethod
537 def random_mt(cls, x=None, scalar_moment=1.0, magnitude=None):
538 '''
539 Create random moment tensor
541 Moment tensors produced by this function appear uniformly distributed
542 when shown in a Hudson's diagram. The rotations used are unifomly
543 distributed in the space of rotations.
544 '''
545 return MomentTensor(
546 m=random_mt(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
548 @classmethod
549 def from_values(cls, values):
550 '''
551 Alternative constructor for moment tensor objects
553 This constructor takes a :py:class:`MomentTensor` object, a tuple, list
554 or NumPy array with 3x3 or 3, 4, 6, or 7 elements to build a Moment
555 tensor object.
557 The ``values`` argument is interpreted depending on shape and type as
558 follows:
560 * ``(strike, dip, rake)``
561 * ``(strike, dip, rake, magnitude)``
562 * ``(mnn, mee, mdd, mne, mnd, med)``
563 * ``(mnn, mee, mdd, mne, mnd, med, magnitude)``
564 * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))``
565 * :py:class:`MomentTensor` object
566 '''
568 m = values_to_matrix(values)
569 return MomentTensor(m=m)
571 def __init__(
572 self, m=None, m_up_south_east=None, m_east_north_up=None,
573 strike=0., dip=0., rake=0., scalar_moment=1.,
574 mnn=None, mee=None, mdd=None, mne=None, mnd=None, med=None,
575 strike1=None, dip1=None, rake1=None,
576 strike2=None, dip2=None, rake2=None,
577 p_axis=None, t_axis=None,
578 magnitude=None, moment=None):
580 Object.__init__(self, init_props=False)
582 if any(mxx is not None for mxx in (mnn, mee, mdd, mne, mnd, med)):
583 m = symmat6(mnn, mee, mdd, mne, mnd, med)
585 if p_axis is not None and t_axis is not None:
586 strike, dip, rake = pt_axes_to_strike_dip_rake(p_axis, t_axis)
588 strike = d2r*strike
589 dip = d2r*dip
590 rake = d2r*rake
592 if isinstance(m, num.matrix):
593 m = m.A
595 if isinstance(m_up_south_east, num.matrix):
596 m_up_south_east = m_up_south_east.A
598 if isinstance(m_east_north_up, num.matrix):
599 m_east_north_up = m_east_north_up.A
601 if m_up_south_east is not None:
602 m = num.dot(
603 self._to_up_south_east,
604 num.dot(m_up_south_east, self._to_up_south_east.T))
606 if m_east_north_up is not None:
607 m = num.dot(
608 self._to_east_north_up,
609 num.dot(m_east_north_up, self._to_east_north_up.T))
611 if m is None:
612 if any(x is not None for x in (
613 strike1, dip1, rake1, strike2, dip2, rake2)):
614 raise Exception(
615 'strike1, dip1, rake1, strike2, dip2, rake2 are read-only '
616 'properties')
618 if moment is not None:
619 scalar_moment = moment
621 if magnitude is not None:
622 scalar_moment = magnitude_to_moment(magnitude)
624 rotmat1 = euler_to_matrix(dip, strike, -rake)
625 m = scalar_moment * num.dot(
626 rotmat1.T,
627 num.dot(MomentTensor._m_unrot, rotmat1))
629 self._m = m
630 self._update()
632 def _update(self):
633 m_evals, m_evecs = eigh_check(self._m)
634 if num.linalg.det(m_evecs) < 0.:
635 m_evecs *= -1.
637 rotmat1 = num.dot(m_evecs, MomentTensor._u_evecs.T).T
638 if num.linalg.det(rotmat1) < 0.:
639 rotmat1 *= -1.
641 self._m_eigenvals = m_evals
642 self._m_eigenvecs = m_evecs
644 self._rotmats = sorted(
645 [rotmat1, num.dot(MomentTensor._flip_dc, rotmat1)],
646 key=lambda m: num.abs(m.flat).tolist())
648 @property
649 def mnn(self):
650 return float(self._m[0, 0])
652 @mnn.setter
653 def mnn(self, value):
654 self._m[0, 0] = value
655 self._update()
657 @property
658 def mee(self):
659 return float(self._m[1, 1])
661 @mee.setter
662 def mee(self, value):
663 self._m[1, 1] = value
664 self._update()
666 @property
667 def mdd(self):
668 return float(self._m[2, 2])
670 @mdd.setter
671 def mdd(self, value):
672 self._m[2, 2] = value
673 self._update()
675 @property
676 def mne(self):
677 return float(self._m[0, 1])
679 @mne.setter
680 def mne(self, value):
681 self._m[0, 1] = value
682 self._m[1, 0] = value
683 self._update()
685 @property
686 def mnd(self):
687 return float(self._m[0, 2])
689 @mnd.setter
690 def mnd(self, value):
691 self._m[0, 2] = value
692 self._m[2, 0] = value
693 self._update()
695 @property
696 def med(self):
697 return float(self._m[1, 2])
699 @med.setter
700 def med(self, value):
701 self._m[1, 2] = value
702 self._m[2, 1] = value
703 self._update()
705 @property
706 def strike1(self):
707 return float(self.both_strike_dip_rake()[0][0])
709 @property
710 def dip1(self):
711 return float(self.both_strike_dip_rake()[0][1])
713 @property
714 def rake1(self):
715 return float(self.both_strike_dip_rake()[0][2])
717 @property
718 def strike2(self):
719 return float(self.both_strike_dip_rake()[1][0])
721 @property
722 def dip2(self):
723 return float(self.both_strike_dip_rake()[1][1])
725 @property
726 def rake2(self):
727 return float(self.both_strike_dip_rake()[1][2])
729 def both_strike_dip_rake(self):
730 '''
731 Get both possible (strike,dip,rake) triplets.
732 '''
733 results = []
734 for rotmat in self._rotmats:
735 alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)]
736 results.append((beta, alpha, -gamma))
738 return results
740 def p_axis(self):
741 '''
742 Get direction of p axis.
743 '''
744 return (self._m_eigenvecs.T)[0]
746 def t_axis(self):
747 '''
748 Get direction of t axis.
749 '''
750 return (self._m_eigenvecs.T)[2]
752 def null_axis(self):
753 '''
754 Get diretion of the null axis.
755 '''
756 return self._m_eigenvecs.T[1]
758 def eigenvals(self):
759 '''
760 Get the eigenvalues of the moment tensor in accending order.
762 :returns: ``(ep, en, et)``
763 '''
765 return self._m_eigenvals
767 def eigensystem(self):
768 '''
769 Get the eigenvalues and eigenvectors of the moment tensor.
771 :returns: ``(ep, en, et, vp, vn, vt)``
772 '''
774 vp = self.p_axis().flatten()
775 vn = self.null_axis().flatten()
776 vt = self.t_axis().flatten()
777 ep, en, et = self._m_eigenvals
778 return ep, en, et, vp, vn, vt
780 def both_slip_vectors(self):
781 '''
782 Get both possible slip directions.
783 '''
784 return [
785 num.dot(rotmat.T, cvec(1., 0., 0.)) for rotmat in self._rotmats]
787 def m(self):
788 '''
789 Get plain moment tensor as 3x3 matrix.
790 '''
791 return self._m.copy()
793 def m6(self):
794 '''
795 Get the moment tensor as a six-element array.
797 :returns: ``(mnn, mee, mdd, mne, mnd, med)``
798 '''
799 return to6(self._m)
801 def m_up_south_east(self):
802 '''
803 Get moment tensor in up-south-east convention as 3x3 matrix.
805 .. math::
806 :nowrap:
808 \\begin{align*}
809 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
810 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
811 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
812 \\end{align*}
813 '''
815 return num.dot(
816 self._to_up_south_east.T,
817 num.dot(self._m, self._to_up_south_east))
819 def m_east_north_up(self):
820 '''
821 Get moment tensor in east-north-up convention as 3x3 matrix.
822 '''
824 return num.dot(
825 self._to_east_north_up.T, num.dot(self._m, self._to_east_north_up))
827 def m6_up_south_east(self):
828 '''
829 Get moment tensor in up-south-east convention as a six-element array.
831 :returns: ``(muu, mss, mee, mus, mue, mse)``
832 '''
833 return to6(self.m_up_south_east())
835 def m6_east_north_up(self):
836 '''
837 Get moment tensor in east-north-up convention as a six-element array.
839 :returns: ``(mee, mnn, muu, men, meu, mnu)``
840 '''
841 return to6(self.m_east_north_up())
843 def m_plain_double_couple(self):
844 '''
845 Get plain double couple with same scalar moment as moment tensor.
846 '''
847 rotmat1 = self._rotmats[0]
848 m = self.scalar_moment() * num.dot(
849 rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
850 return m
852 def moment_magnitude(self):
853 '''
854 Get moment magnitude of moment tensor.
855 '''
856 return moment_to_magnitude(self.scalar_moment())
858 def scalar_moment(self):
859 '''
860 Get the scalar moment of the moment tensor (Frobenius norm based)
862 .. math::
864 M_0 = \\frac{1}{\\sqrt{2}}\\sqrt{\\sum_{i,j} |M_{ij}|^2}
866 The scalar moment is calculated based on the Euclidean (Frobenius) norm
867 (Silver and Jordan, 1982). The scalar moment returned by this function
868 differs from the standard decomposition based definition of the scalar
869 moment for non-double-couple moment tensors.
870 '''
871 return num.linalg.norm(self._m_eigenvals) / math.sqrt(2.)
873 @property
874 def moment(self):
875 return float(self.scalar_moment())
877 @moment.setter
878 def moment(self, value):
879 self._m *= value / self.moment
880 self._update()
882 @property
883 def magnitude(self):
884 return float(self.moment_magnitude())
886 @magnitude.setter
887 def magnitude(self, value):
888 self._m *= magnitude_to_moment(value) / self.moment
889 self._update()
891 def __str__(self):
892 mexp = pow(10, math.ceil(num.log10(num.max(num.abs(self._m)))))
893 m = self._m / mexp
894 s = '''Scalar Moment [Nm]: M0 = %g (Mw = %3.1f)
895Moment Tensor [Nm]: Mnn = %6.3f, Mee = %6.3f, Mdd = %6.3f,
896 Mne = %6.3f, Mnd = %6.3f, Med = %6.3f [ x %g ]
897''' % (
898 self.scalar_moment(),
899 self.moment_magnitude(),
900 m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2],
901 mexp)
903 s += self.str_fault_planes()
904 return s
906 def str_fault_planes(self):
907 s = ''
908 for i, sdr in enumerate(self.both_strike_dip_rake()):
909 s += 'Fault plane %i [deg]: ' \
910 'strike = %3.0f, dip = %3.0f, slip-rake = %4.0f\n' \
911 % (i+1, sdr[0], sdr[1], sdr[2])
913 return s
915 def deviatoric(self):
916 '''
917 Get deviatoric part of moment tensor.
919 Returns a new moment tensor object with zero trace.
920 '''
922 m = self.m()
924 trace_m = num.trace(m)
925 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
926 m_devi = m - m_iso
927 mt = MomentTensor(m=m_devi)
928 return mt
930 def standard_decomposition(self):
931 '''
932 Decompose moment tensor into isotropic, DC and CLVD components.
934 Standard decomposition according to e.g. Jost and Herrmann 1989 is
935 returned as::
937 [
938 (moment_iso, ratio_iso, m_iso),
939 (moment_dc, ratio_dc, m_dc),
940 (moment_clvd, ratio_clvd, m_clvd),
941 (moment_devi, ratio_devi, m_devi),
942 (moment, 1.0, m)
943 ]
944 '''
946 epsilon = 1e-6
948 m = self.m()
950 trace_m = num.trace(m)
951 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
952 moment_iso = abs(trace_m / 3.)
954 m_devi = m - m_iso
956 evals, evecs = eigh_check(m_devi)
958 moment_devi = num.max(num.abs(evals))
959 moment = moment_iso + moment_devi
961 iorder = num.argsort(num.abs(evals))
962 evals_sorted = evals[iorder]
963 evecs_sorted = (evecs.T[iorder]).T
965 if moment_devi < epsilon * moment_iso:
966 signed_moment_dc = 0.
967 else:
968 assert -epsilon <= -evals_sorted[0] / evals_sorted[2] <= 0.5
969 signed_moment_dc = evals_sorted[2] * (1.0 + 2.0 * (
970 min(0.0, evals_sorted[0] / evals_sorted[2])))
972 moment_dc = abs(signed_moment_dc)
973 m_dc_es = num.dot(signed_moment_dc, num.diag([0., -1.0, 1.0]))
974 m_dc = num.dot(evecs_sorted, num.dot(m_dc_es, evecs_sorted.T))
976 m_clvd = m_devi - m_dc
978 moment_clvd = moment_devi - moment_dc
980 ratio_dc = moment_dc / moment
981 ratio_clvd = moment_clvd / moment
982 ratio_iso = moment_iso / moment
983 ratio_devi = moment_devi / moment
985 return [
986 (moment_iso, ratio_iso, m_iso),
987 (moment_dc, ratio_dc, m_dc),
988 (moment_clvd, ratio_clvd, m_clvd),
989 (moment_devi, ratio_devi, m_devi),
990 (moment, 1.0, m)]
992 def rotated(self, rot):
993 '''
994 Get rotated moment tensor.
996 :param rot: ratation matrix, coordinate system is NED
997 :returns: new :py:class:`MomentTensor` object
998 '''
1000 rotmat = num.array(rot)
1001 return MomentTensor(m=num.dot(rotmat, num.dot(self.m(), rotmat.T)))
1003 def random_rotated(self, angle_std=None, angle=None, rstate=None):
1004 '''
1005 Get distorted MT by rotation around random axis and angle.
1007 :param angle_std: angles are drawn from a normal distribution with
1008 zero mean and given standard deviation [degrees]
1009 :param angle: set angle [degrees], only axis will be random
1010 :param rstate: :py:class:`numpy.random.RandomState` object, can be
1011 used to create reproducible pseudo-random sequences
1012 :returns: new :py:class:`MomentTensor` object
1013 '''
1015 assert (angle_std is None) != (angle is None), \
1016 'either angle or angle_std must be given'
1018 if angle_std is not None:
1019 rstate = rstate or num.random
1020 angle = rstate.normal(scale=angle_std)
1022 axis = random_axis(rstate=rstate)
1023 rot = rotation_from_angle_and_axis(angle, axis)
1024 return self.rotated(rot)
1027def other_plane(strike, dip, rake):
1028 '''
1029 Get the respectively other plane in the double-couple ambiguity.
1030 '''
1032 mt = MomentTensor(strike=strike, dip=dip, rake=rake)
1033 both_sdr = mt.both_strike_dip_rake()
1034 w = [sum([abs(x-y) for x, y in zip(both_sdr[i], (strike, dip, rake))])
1035 for i in (0, 1)]
1037 if w[0] < w[1]:
1038 return both_sdr[1]
1039 else:
1040 return both_sdr[0]
1043def dsdr(sdr1, sdr2):
1044 s1, d1, r1 = sdr1
1045 s2, d2, r2 = sdr2
1047 s1 = s1 % 360.
1048 s2 = s2 % 360.
1049 r1 = r1 % 360.
1050 r2 = r2 % 360.
1052 ds = abs(s1 - s2)
1053 ds = ds if ds <= 180. else 360. - ds
1055 dr = abs(r1 - r2)
1056 dr = dr if dr <= 180. else 360. - dr
1058 dd = abs(d1 - d2)
1060 return math.sqrt(ds**2 + dr**2 + dd**2)
1063def order_like(sdrs, sdrs_ref):
1064 '''
1065 Order strike-dip-rake pair post closely to a given reference pair.
1067 :param sdrs: tuple, ``((strike1, dip1, rake1), (strike2, dip2, rake2))``
1068 :param sdrs_ref: as above but with reference pair
1069 '''
1071 d1 = min(dsdr(sdrs[0], sdrs_ref[0]), dsdr(sdrs[1], sdrs_ref[1]))
1072 d2 = min(dsdr(sdrs[0], sdrs_ref[1]), dsdr(sdrs[1], sdrs_ref[0]))
1073 if d1 < d2:
1074 return sdrs
1075 else:
1076 return sdrs[::-1]
1079def _tpb2q(t, p, b):
1080 eps = 0.001
1081 tqw = 1. + t[0] + p[1] + b[2]
1082 tqx = 1. + t[0] - p[1] - b[2]
1083 tqy = 1. - t[0] + p[1] - b[2]
1084 tqz = 1. - t[0] - p[1] + b[2]
1086 q = num.zeros(4)
1087 if tqw > eps:
1088 q[0] = 0.5 * math.sqrt(tqw)
1089 q[1:] = p[2] - b[1], b[0] - t[2], t[1] - p[0]
1090 elif tqx > eps:
1091 q[0] = 0.5 * math.sqrt(tqx)
1092 q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1093 elif tqy > eps:
1094 q[0] = 0.5 * math.sqrt(tqy)
1095 q[1:] = b[0] - t[2], p[0] + t[1], b[1] + p[2]
1096 elif tqz > eps:
1097 q[0] = 0.5 * math.sqrt(tqz)
1098 q[1:] = t[1] - p[0], b[0] + t[2], b[1] + p[2]
1099 else:
1100 raise Exception('should not reach this line, check theory!')
1101 # q[0] = max(0.5 * math.sqrt(tqx), eps)
1102 # q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1104 q[1:] /= 4.0 * q[0]
1106 q /= math.sqrt(num.sum(q**2))
1108 return q
1111_pbt2tpb = num.array(((0., 0., 1.), (1., 0., 0.), (0., 1., 0.)))
1114def kagan_angle(mt1, mt2):
1115 '''
1116 Given two moment tensors, return the Kagan angle in degrees.
1118 After Kagan (1991) and Tape & Tape (2012).
1119 '''
1121 ai = num.dot(_pbt2tpb, mt1._m_eigenvecs.T)
1122 aj = num.dot(_pbt2tpb, mt2._m_eigenvecs.T)
1124 u = num.dot(ai, aj.T)
1126 tk, pk, bk = u
1128 qk = _tpb2q(tk, pk, bk)
1130 return 2. * r2d * math.acos(num.max(num.abs(qk)))
1133def rand_to_gutenberg_richter(rand, b_value, magnitude_min):
1134 '''
1135 Draw magnitude from Gutenberg Richter distribution.
1136 '''
1137 return magnitude_min + num.log10(1.-rand) / -b_value
1140def magnitude_to_duration_gcmt(magnitudes):
1141 '''
1142 Scaling relation used by Global CMT catalog for most of its events.
1143 '''
1145 mom = magnitude_to_moment(magnitudes)
1146 return (mom / 1.1e16)**(1./3.)
1149if __name__ == '__main__':
1151 import sys
1152 v = list(map(float, sys.argv[1:]))
1153 mt = MomentTensor.from_values(v)
1154 print(mt)