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)
462def pt_axes_to_strike_dip_rake(p_axis, t_axis):
463 n_axis = num.cross(p_axis, t_axis)
464 m_evecs = num.vstack([p_axis, n_axis, t_axis]).T
465 if num.linalg.det(m_evecs) < 0.:
466 m_evecs *= -1.
467 rotmat = num.dot(m_evecs, MomentTensor._u_evecs.T).T
468 if num.linalg.det(rotmat) < 0.:
469 rotmat *= -1.
470 alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)]
471 return (beta, alpha, -gamma)
474class MomentTensor(Object):
476 '''
477 Moment tensor object
479 :param m: NumPy array in north-east-down convention
480 :param m_up_south_east: NumPy array in up-south-east convention
481 :param m_east_north_up: NumPy array in east-north-up convention
482 :param strike,dip,rake: fault plane angles in [degrees]
483 :param p_axis,t_axis: initialize double-couple from p and t axes
484 :param scalar_moment: scalar moment in [Nm]
485 :param magnitude: moment magnitude Mw
487 Global CMT catalog moment tensors use the up-south-east (USE) coordinate
488 system convention with :math:`r` (up), :math:`\\theta` (south), and
489 :math:`\\phi` (east).
491 .. math::
492 :nowrap:
494 \\begin{align*}
495 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
496 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
497 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
498 \\end{align*}
500 '''
502 mnn__ = Float.T(default=0.0)
503 mee__ = Float.T(default=0.0)
504 mdd__ = Float.T(default=0.0)
505 mne__ = Float.T(default=0.0)
506 mnd__ = Float.T(default=-1.0)
507 med__ = Float.T(default=0.0)
508 strike1__ = Float.T(default=None, optional=True) # read-only
509 dip1__ = Float.T(default=None, optional=True) # read-only
510 rake1__ = Float.T(default=None, optional=True) # read-only
511 strike2__ = Float.T(default=None, optional=True) # read-only
512 dip2__ = Float.T(default=None, optional=True) # read-only
513 rake2__ = Float.T(default=None, optional=True) # read-only
514 moment__ = Float.T(default=None, optional=True) # read-only
515 magnitude__ = Float.T(default=None, optional=True) # read-only
517 _flip_dc = num.array(
518 [[0., 0., -1.], [0., -1., 0.], [-1., 0., 0.]], dtype=float)
519 _to_up_south_east = num.array(
520 [[0., 0., -1.], [-1., 0., 0.], [0., 1., 0.]], dtype=float).T
521 _to_east_north_up = num.array(
522 [[0., 1., 0.], [1., 0., 0.], [0., 0., -1.]], dtype=float)
523 _m_unrot = num.array(
524 [[0., 0., -1.], [0., 0., 0.], [-1., 0., 0.]], dtype=float)
526 _u_evals, _u_evecs = eigh_check(_m_unrot)
528 @classmethod
529 def random_dc(cls, x=None, scalar_moment=1.0, magnitude=None):
530 '''
531 Create random oriented double-couple moment tensor
533 The rotations used are uniformly distributed in the space of rotations.
534 '''
535 return MomentTensor(
536 m=random_dc(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
538 @classmethod
539 def random_mt(cls, x=None, scalar_moment=1.0, magnitude=None):
540 '''
541 Create random moment tensor
543 Moment tensors produced by this function appear uniformly distributed
544 when shown in a Hudson's diagram. The rotations used are unifomly
545 distributed in the space of rotations.
546 '''
547 return MomentTensor(
548 m=random_mt(x=x, scalar_moment=scalar_moment, magnitude=magnitude))
550 @classmethod
551 def from_values(cls, values):
552 '''
553 Alternative constructor for moment tensor objects
555 This constructor takes a :py:class:`MomentTensor` object, a tuple, list
556 or NumPy array with 3x3 or 3, 4, 6, or 7 elements to build a Moment
557 tensor object.
559 The ``values`` argument is interpreted depending on shape and type as
560 follows:
562 * ``(strike, dip, rake)``
563 * ``(strike, dip, rake, magnitude)``
564 * ``(mnn, mee, mdd, mne, mnd, med)``
565 * ``(mnn, mee, mdd, mne, mnd, med, magnitude)``
566 * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))``
567 * :py:class:`MomentTensor` object
568 '''
570 m = values_to_matrix(values)
571 return MomentTensor(m=m)
573 def __init__(
574 self, m=None, m_up_south_east=None, m_east_north_up=None,
575 strike=0., dip=0., rake=0., scalar_moment=1.,
576 mnn=None, mee=None, mdd=None, mne=None, mnd=None, med=None,
577 strike1=None, dip1=None, rake1=None,
578 strike2=None, dip2=None, rake2=None,
579 p_axis=None, t_axis=None,
580 magnitude=None, moment=None):
582 Object.__init__(self, init_props=False)
584 if any(mxx is not None for mxx in (mnn, mee, mdd, mne, mnd, med)):
585 m = symmat6(mnn, mee, mdd, mne, mnd, med)
587 if p_axis is not None and t_axis is not None:
588 strike, dip, rake = pt_axes_to_strike_dip_rake(p_axis, t_axis)
590 strike = d2r*strike
591 dip = d2r*dip
592 rake = d2r*rake
594 if isinstance(m, num.matrix):
595 m = m.A
597 if isinstance(m_up_south_east, num.matrix):
598 m_up_south_east = m_up_south_east.A
600 if isinstance(m_east_north_up, num.matrix):
601 m_east_north_up = m_east_north_up.A
603 if m_up_south_east is not None:
604 m = num.dot(
605 self._to_up_south_east,
606 num.dot(m_up_south_east, self._to_up_south_east.T))
608 if m_east_north_up is not None:
609 m = num.dot(
610 self._to_east_north_up,
611 num.dot(m_east_north_up, self._to_east_north_up.T))
613 if m is None:
614 if any(x is not None for x in (
615 strike1, dip1, rake1, strike2, dip2, rake2)):
616 raise Exception(
617 'strike1, dip1, rake1, strike2, dip2, rake2 are read-only '
618 'properties')
620 if moment is not None:
621 scalar_moment = moment
623 if magnitude is not None:
624 scalar_moment = magnitude_to_moment(magnitude)
626 rotmat1 = euler_to_matrix(dip, strike, -rake)
627 m = scalar_moment * num.dot(
628 rotmat1.T,
629 num.dot(MomentTensor._m_unrot, rotmat1))
631 self._m = m
632 self._update()
634 def _update(self):
635 m_evals, m_evecs = eigh_check(self._m)
636 if num.linalg.det(m_evecs) < 0.:
637 m_evecs *= -1.
639 rotmat1 = num.dot(m_evecs, MomentTensor._u_evecs.T).T
640 if num.linalg.det(rotmat1) < 0.:
641 rotmat1 *= -1.
643 self._m_eigenvals = m_evals
644 self._m_eigenvecs = m_evecs
646 self._rotmats = sorted(
647 [rotmat1, num.dot(MomentTensor._flip_dc, rotmat1)],
648 key=lambda m: num.abs(m.flat).tolist())
650 @property
651 def mnn(self):
652 return float(self._m[0, 0])
654 @mnn.setter
655 def mnn(self, value):
656 self._m[0, 0] = value
657 self._update()
659 @property
660 def mee(self):
661 return float(self._m[1, 1])
663 @mee.setter
664 def mee(self, value):
665 self._m[1, 1] = value
666 self._update()
668 @property
669 def mdd(self):
670 return float(self._m[2, 2])
672 @mdd.setter
673 def mdd(self, value):
674 self._m[2, 2] = value
675 self._update()
677 @property
678 def mne(self):
679 return float(self._m[0, 1])
681 @mne.setter
682 def mne(self, value):
683 self._m[0, 1] = value
684 self._m[1, 0] = value
685 self._update()
687 @property
688 def mnd(self):
689 return float(self._m[0, 2])
691 @mnd.setter
692 def mnd(self, value):
693 self._m[0, 2] = value
694 self._m[2, 0] = value
695 self._update()
697 @property
698 def med(self):
699 return float(self._m[1, 2])
701 @med.setter
702 def med(self, value):
703 self._m[1, 2] = value
704 self._m[2, 1] = value
705 self._update()
707 @property
708 def strike1(self):
709 return float(self.both_strike_dip_rake()[0][0])
711 @property
712 def dip1(self):
713 return float(self.both_strike_dip_rake()[0][1])
715 @property
716 def rake1(self):
717 return float(self.both_strike_dip_rake()[0][2])
719 @property
720 def strike2(self):
721 return float(self.both_strike_dip_rake()[1][0])
723 @property
724 def dip2(self):
725 return float(self.both_strike_dip_rake()[1][1])
727 @property
728 def rake2(self):
729 return float(self.both_strike_dip_rake()[1][2])
731 def both_strike_dip_rake(self):
732 '''
733 Get both possible (strike,dip,rake) triplets.
734 '''
735 results = []
736 for rotmat in self._rotmats:
737 alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)]
738 results.append((beta, alpha, -gamma))
740 return results
742 def p_axis(self):
743 '''
744 Get direction of p axis.
745 '''
746 return (self._m_eigenvecs.T)[0]
748 def t_axis(self):
749 '''
750 Get direction of t axis.
751 '''
752 return (self._m_eigenvecs.T)[2]
754 def null_axis(self):
755 '''
756 Get diretion of the null axis.
757 '''
758 return self._m_eigenvecs.T[1]
760 def eigenvals(self):
761 '''
762 Get the eigenvalues of the moment tensor in accending order.
764 :returns: ``(ep, en, et)``
765 '''
767 return self._m_eigenvals
769 def eigensystem(self):
770 '''
771 Get the eigenvalues and eigenvectors of the moment tensor.
773 :returns: ``(ep, en, et, vp, vn, vt)``
774 '''
776 vp = self.p_axis().flatten()
777 vn = self.null_axis().flatten()
778 vt = self.t_axis().flatten()
779 ep, en, et = self._m_eigenvals
780 return ep, en, et, vp, vn, vt
782 def both_slip_vectors(self):
783 '''
784 Get both possible slip directions.
785 '''
786 return [
787 num.dot(rotmat.T, cvec(1., 0., 0.)) for rotmat in self._rotmats]
789 def m(self):
790 '''
791 Get plain moment tensor as 3x3 matrix.
792 '''
793 return self._m.copy()
795 def m6(self):
796 '''
797 Get the moment tensor as a six-element array.
799 :returns: ``(mnn, mee, mdd, mne, mnd, med)``
800 '''
801 return to6(self._m)
803 def m_up_south_east(self):
804 '''
805 Get moment tensor in up-south-east convention as 3x3 matrix.
807 .. math::
808 :nowrap:
810 \\begin{align*}
811 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\
812 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\
813 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne}
814 \\end{align*}
815 '''
817 return num.dot(
818 self._to_up_south_east.T,
819 num.dot(self._m, self._to_up_south_east))
821 def m_east_north_up(self):
822 '''
823 Get moment tensor in east-north-up convention as 3x3 matrix.
824 '''
826 return num.dot(
827 self._to_east_north_up.T, num.dot(self._m, self._to_east_north_up))
829 def m6_up_south_east(self):
830 '''
831 Get moment tensor in up-south-east convention as a six-element array.
833 :returns: ``(muu, mss, mee, mus, mue, mse)``
834 '''
835 return to6(self.m_up_south_east())
837 def m6_east_north_up(self):
838 '''
839 Get moment tensor in east-north-up convention as a six-element array.
841 :returns: ``(mee, mnn, muu, men, meu, mnu)``
842 '''
843 return to6(self.m_east_north_up())
845 def m_plain_double_couple(self):
846 '''
847 Get plain double couple with same scalar moment as moment tensor.
848 '''
849 rotmat1 = self._rotmats[0]
850 m = self.scalar_moment() * num.dot(
851 rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1))
852 return m
854 def moment_magnitude(self):
855 '''
856 Get moment magnitude of moment tensor.
857 '''
858 return moment_to_magnitude(self.scalar_moment())
860 def scalar_moment(self):
861 '''
862 Get the scalar moment of the moment tensor (Frobenius norm based)
864 .. math::
866 M_0 = \\frac{1}{\\sqrt{2}}\\sqrt{\\sum_{i,j} |M_{ij}|^2}
868 The scalar moment is calculated based on the Euclidean (Frobenius) norm
869 (Silver and Jordan, 1982). The scalar moment returned by this function
870 differs from the standard decomposition based definition of the scalar
871 moment for non-double-couple moment tensors.
872 '''
873 return num.linalg.norm(self._m_eigenvals) / math.sqrt(2.)
875 @property
876 def moment(self):
877 return float(self.scalar_moment())
879 @moment.setter
880 def moment(self, value):
881 self._m *= value / self.moment
882 self._update()
884 @property
885 def magnitude(self):
886 return float(self.moment_magnitude())
888 @magnitude.setter
889 def magnitude(self, value):
890 self._m *= magnitude_to_moment(value) / self.moment
891 self._update()
893 def __str__(self):
894 mexp = pow(10, math.ceil(num.log10(num.max(num.abs(self._m)))))
895 m = self._m / mexp
896 s = '''Scalar Moment [Nm]: M0 = %g (Mw = %3.1f)
897Moment Tensor [Nm]: Mnn = %6.3f, Mee = %6.3f, Mdd = %6.3f,
898 Mne = %6.3f, Mnd = %6.3f, Med = %6.3f [ x %g ]
899''' % (
900 self.scalar_moment(),
901 self.moment_magnitude(),
902 m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2],
903 mexp)
905 s += self.str_fault_planes()
906 return s
908 def str_fault_planes(self):
909 s = ''
910 for i, sdr in enumerate(self.both_strike_dip_rake()):
911 s += 'Fault plane %i [deg]: ' \
912 'strike = %3.0f, dip = %3.0f, slip-rake = %4.0f\n' \
913 % (i+1, sdr[0], sdr[1], sdr[2])
915 return s
917 def deviatoric(self):
918 '''
919 Get deviatoric part of moment tensor.
921 Returns a new moment tensor object with zero trace.
922 '''
924 m = self.m()
926 trace_m = num.trace(m)
927 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
928 m_devi = m - m_iso
929 mt = MomentTensor(m=m_devi)
930 return mt
932 def standard_decomposition(self):
933 '''
934 Decompose moment tensor into isotropic, DC and CLVD components.
936 Standard decomposition according to e.g. Jost and Herrmann 1989 is
937 returned as::
939 [
940 (moment_iso, ratio_iso, m_iso),
941 (moment_dc, ratio_dc, m_dc),
942 (moment_clvd, ratio_clvd, m_clvd),
943 (moment_devi, ratio_devi, m_devi),
944 (moment, 1.0, m)
945 ]
946 '''
948 epsilon = 1e-6
950 m = self.m()
952 trace_m = num.trace(m)
953 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
954 moment_iso = abs(trace_m / 3.)
956 m_devi = m - m_iso
958 evals, evecs = eigh_check(m_devi)
960 moment_devi = num.max(num.abs(evals))
961 moment = moment_iso + moment_devi
963 iorder = num.argsort(num.abs(evals))
964 evals_sorted = evals[iorder]
965 evecs_sorted = (evecs.T[iorder]).T
967 if moment_devi < epsilon * moment_iso:
968 signed_moment_dc = 0.
969 else:
970 assert -epsilon <= -evals_sorted[0] / evals_sorted[2] <= 0.5
971 signed_moment_dc = evals_sorted[2] * (1.0 + 2.0 * (
972 min(0.0, evals_sorted[0] / evals_sorted[2])))
974 moment_dc = abs(signed_moment_dc)
975 m_dc_es = num.dot(signed_moment_dc, num.diag([0., -1.0, 1.0]))
976 m_dc = num.dot(evecs_sorted, num.dot(m_dc_es, evecs_sorted.T))
978 m_clvd = m_devi - m_dc
980 moment_clvd = moment_devi - moment_dc
982 ratio_dc = moment_dc / moment
983 ratio_clvd = moment_clvd / moment
984 ratio_iso = moment_iso / moment
985 ratio_devi = moment_devi / moment
987 return [
988 (moment_iso, ratio_iso, m_iso),
989 (moment_dc, ratio_dc, m_dc),
990 (moment_clvd, ratio_clvd, m_clvd),
991 (moment_devi, ratio_devi, m_devi),
992 (moment, 1.0, m)]
994 def rotated(self, rot):
995 '''
996 Get rotated moment tensor.
998 :param rot: ratation matrix, coordinate system is NED
999 :returns: new :py:class:`MomentTensor` object
1000 '''
1002 rotmat = num.array(rot)
1003 return MomentTensor(m=num.dot(rotmat, num.dot(self.m(), rotmat.T)))
1005 def random_rotated(self, angle_std=None, angle=None, rstate=None):
1006 '''
1007 Get distorted MT by rotation around random axis and angle.
1009 :param angle_std: angles are drawn from a normal distribution with
1010 zero mean and given standard deviation [degrees]
1011 :param angle: set angle [degrees], only axis will be random
1012 :param rstate: :py:class:`numpy.random.RandomState` object, can be
1013 used to create reproducible pseudo-random sequences
1014 :returns: new :py:class:`MomentTensor` object
1015 '''
1017 assert (angle_std is None) != (angle is None), \
1018 'either angle or angle_std must be given'
1020 if angle_std is not None:
1021 rstate = rstate or num.random
1022 angle = rstate.normal(scale=angle_std)
1024 axis = random_axis(rstate=rstate)
1025 rot = rotation_from_angle_and_axis(angle, axis)
1026 return self.rotated(rot)
1029def other_plane(strike, dip, rake):
1030 '''
1031 Get the respectively other plane in the double-couple ambiguity.
1032 '''
1034 mt = MomentTensor(strike=strike, dip=dip, rake=rake)
1035 both_sdr = mt.both_strike_dip_rake()
1036 w = [sum([abs(x-y) for x, y in zip(both_sdr[i], (strike, dip, rake))])
1037 for i in (0, 1)]
1039 if w[0] < w[1]:
1040 return both_sdr[1]
1041 else:
1042 return both_sdr[0]
1045def dsdr(sdr1, sdr2):
1046 s1, d1, r1 = sdr1
1047 s2, d2, r2 = sdr2
1049 s1 = s1 % 360.
1050 s2 = s2 % 360.
1051 r1 = r1 % 360.
1052 r2 = r2 % 360.
1054 ds = abs(s1 - s2)
1055 ds = ds if ds <= 180. else 360. - ds
1057 dr = abs(r1 - r2)
1058 dr = dr if dr <= 180. else 360. - dr
1060 dd = abs(d1 - d2)
1062 return math.sqrt(ds**2 + dr**2 + dd**2)
1065def order_like(sdrs, sdrs_ref):
1066 '''
1067 Order strike-dip-rake pair post closely to a given reference pair.
1069 :param sdrs: tuple, ``((strike1, dip1, rake1), (strike2, dip2, rake2))``
1070 :param sdrs_ref: as above but with reference pair
1071 '''
1073 d1 = min(dsdr(sdrs[0], sdrs_ref[0]), dsdr(sdrs[1], sdrs_ref[1]))
1074 d2 = min(dsdr(sdrs[0], sdrs_ref[1]), dsdr(sdrs[1], sdrs_ref[0]))
1075 if d1 < d2:
1076 return sdrs
1077 else:
1078 return sdrs[::-1]
1081def _tpb2q(t, p, b):
1082 eps = 0.001
1083 tqw = 1. + t[0] + p[1] + b[2]
1084 tqx = 1. + t[0] - p[1] - b[2]
1085 tqy = 1. - t[0] + p[1] - b[2]
1086 tqz = 1. - t[0] - p[1] + b[2]
1088 q = num.zeros(4)
1089 if tqw > eps:
1090 q[0] = 0.5 * math.sqrt(tqw)
1091 q[1:] = p[2] - b[1], b[0] - t[2], t[1] - p[0]
1092 elif tqx > eps:
1093 q[0] = 0.5 * math.sqrt(tqx)
1094 q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1095 elif tqy > eps:
1096 q[0] = 0.5 * math.sqrt(tqy)
1097 q[1:] = b[0] - t[2], p[0] + t[1], b[1] + p[2]
1098 elif tqz > eps:
1099 q[0] = 0.5 * math.sqrt(tqz)
1100 q[1:] = t[1] - p[0], b[0] + t[2], b[1] + p[2]
1101 else:
1102 raise Exception('should not reach this line, check theory!')
1103 # q[0] = max(0.5 * math.sqrt(tqx), eps)
1104 # q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2]
1106 q[1:] /= 4.0 * q[0]
1108 q /= math.sqrt(num.sum(q**2))
1110 return q
1113_pbt2tpb = num.array(((0., 0., 1.), (1., 0., 0.), (0., 1., 0.)))
1116def kagan_angle(mt1, mt2):
1117 '''
1118 Given two moment tensors, return the Kagan angle in degrees.
1120 After Kagan (1991) and Tape & Tape (2012).
1121 '''
1123 ai = num.dot(_pbt2tpb, mt1._m_eigenvecs.T)
1124 aj = num.dot(_pbt2tpb, mt2._m_eigenvecs.T)
1126 u = num.dot(ai, aj.T)
1128 tk, pk, bk = u
1130 qk = _tpb2q(tk, pk, bk)
1132 return 2. * r2d * math.acos(num.max(num.abs(qk)))
1135def rand_to_gutenberg_richter(rand, b_value, magnitude_min):
1136 '''
1137 Draw magnitude from Gutenberg Richter distribution.
1138 '''
1139 return magnitude_min + num.log10(1.-rand) / -b_value
1142def magnitude_to_duration_gcmt(magnitudes):
1143 '''
1144 Scaling relation used by Global CMT catalog for most of its events.
1145 '''
1147 mom = magnitude_to_moment(magnitudes)
1148 return (mom / 1.1e16)**(1./3.)
1151if __name__ == '__main__':
1153 import sys
1154 v = list(map(float, sys.argv[1:]))
1155 mt = MomentTensor.from_values(v)
1156 print(mt)