1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7This module provides various moment tensor related utility functions. 

8 

9It can be used to convert between strike-dip-rake and moment tensor 

10representations and provides different options to produce random moment 

11tensors. 

12 

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). 

17 

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''' 

23 

24from __future__ import division, print_function, absolute_import 

25 

26import math 

27import numpy as num 

28 

29from .guts import Object, Float 

30 

31guts_prefix = 'pf' 

32 

33dynecm = 1e-7 

34 

35 

36def random_axis(rstate=None): 

37 ''' 

38 Get randomly oriented unit vector. 

39 

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 

49 

50 

51def rotation_from_angle_and_axis(angle, axis): 

52 ''' 

53 Build rotation matrix based on axis and angle. 

54 

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 ''' 

60 

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) 

66 

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 

73 

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 ]) 

81 

82 

83def random_rotation(x=None): 

84 ''' 

85 Get random rotation matrix. 

86 

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". 

90 

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 ''' 

96 

97 if x is not None: 

98 x1, x2, x3 = x 

99 else: 

100 x1, x2, x3 = num.random.random(3) 

101 

102 phi = math.pi*2.0*x1 

103 

104 zrot = num.array([ 

105 [math.cos(phi), math.sin(phi), 0.], 

106 [-math.sin(phi), math.cos(phi), 0.], 

107 [0., 0., 1.]]) 

108 

109 lam = math.pi*2.0*x2 

110 

111 v = num.array([[ 

112 math.cos(lam)*math.sqrt(x3), 

113 math.sin(lam)*math.sqrt(x3), 

114 math.sqrt(1.-x3)]]).T 

115 

116 house = num.identity(3) - 2.0 * num.dot(v, v.T) 

117 return -num.dot(house, zrot) 

118 

119 

120def rand(mi, ma): 

121 return float(num.random.uniform(mi, ma)) 

122 

123 

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 

128 

129 

130def random_strike_dip_rake( 

131 strikemin=0., strikemax=360., 

132 dipmin=0., dipmax=90., 

133 rakemin=-180., rakemax=180.): 

134 

135 ''' 

136 Get random strike, dip, rake triplet. 

137 

138 .. note:: 

139 

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 ''' 

144 

145 strike = rand(strikemin, strikemax) 

146 dip = randdip(dipmin, dipmax) 

147 rake = rand(rakemin, rakemax) 

148 

149 return strike, dip, rake 

150 

151 

152def to6(m): 

153 ''' 

154 Get non-redundant components from symmetric 3x3 matrix. 

155 

156 :returns: 1D NumPy array with entries ordered like 

157 ``(a_xx, a_yy, a_zz, a_xy, a_xz, a_yz)`` 

158 ''' 

159 

160 return num.array([m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2]]) 

161 

162 

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 ''' 

167 

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) 

171 

172 

173def values_to_matrix(values): 

174 ''' 

175 Convert anything to moment tensor represented as a NumPy array. 

176 

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. 

179 

180 The ``values`` argument is interpreted depending on shape and type as 

181 follows: 

182 

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 ''' 

190 

191 if isinstance(values, (tuple, list)): 

192 values = num.asarray(values, dtype=float) 

193 

194 if isinstance(values, MomentTensor): 

195 return values.m() 

196 

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)) 

202 

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)) 

209 

210 elif values.shape == (6,): 

211 return symmat6(*values) 

212 

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 

219 

220 elif values.shape == (3, 3): 

221 return num.asarray(values, dtype=float) 

222 

223 raise Exception('cannot convert object to 3x3 matrix') 

224 

225 

226def moment_to_magnitude(moment): 

227 ''' 

228 Convert scalar moment to moment magnitude Mw. 

229 

230 :param moment: scalar moment [Nm] 

231 :returns: moment magnitude Mw 

232 

233 Moment magnitude is defined as 

234 

235 .. math:: 

236 

237 M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0 \\cdot 10^7) - 10.7 

238 

239 where :math:`M_0` is the scalar moment given in [Nm]. 

240 

241 .. note:: 

242 

243 Global CMT uses 10.7333333 instead of 10.7, based on [Kanamori 1977], 

244 10.7 is from [Hanks and Kanamori 1979]. 

245 ''' 

246 

247 return num.log10(moment*1.0e7) / 1.5 - 10.7 

248 

249 

250def magnitude_to_moment(magnitude): 

251 ''' 

252 Convert moment magnitude Mw to scalar moment. 

253 

254 :param magnitude: moment magnitude 

255 :returns: scalar moment [Nm] 

256 

257 See :py:func:`moment_to_magnitude`. 

258 ''' 

259 

260 return 10.0**(1.5*(magnitude+10.7))*1.0e-7 

261 

262 

263magnitude_1Nm = moment_to_magnitude(1.0) 

264 

265 

266def euler_to_matrix(alpha, beta, gamma): 

267 ''' 

268 Given euler angle triplet, create rotation matrix 

269 

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. 

273 

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] 

278 

279 Usage for moment tensors:: 

280 

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)) 

284 

285 ''' 

286 

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) 

293 

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 

298 

299 

300def matrix_to_euler(rotmat): 

301 ''' 

302 Get eulerian angle triplet from rotation matrix. 

303 ''' 

304 

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. 

316 

317 if cos_alpha < -1.: 

318 cos_alpha = -1. 

319 

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.) 

323 

324 return unique_euler(alpha, beta, gamma) 

325 

326 

327def unique_euler(alpha, beta, gamma): 

328 ''' 

329 Uniquify eulerian angle triplet. 

330 

331 Put eulerian angle triplet into ranges compatible with 

332 ``(dip, strike, -rake)`` conventions in seismology:: 

333 

334 alpha (dip) : [0, pi/2] 

335 beta (strike) : [0, 2*pi) 

336 gamma (-rake) : [-pi, pi) 

337 

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. 

340 

341 If ``alpha`` is near to ``pi/2``, ``beta`` is put into the range 

342 ``[0,pi)``. 

343 ''' 

344 

345 pi = math.pi 

346 

347 alpha = num.mod(alpha, 2.0*pi) 

348 

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 

360 

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 

364 

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) 

368 

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 

372 

373 if abs(beta - pi) < 1e-10: 

374 beta = pi 

375 

376 if abs(beta - 2.*pi) < 1e-10: 

377 beta = 0. 

378 

379 if abs(beta) < 1e-10: 

380 beta = 0. 

381 

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 

388 

389 if alpha < 1e-7: 

390 beta = num.mod(beta + gamma, 2.0*pi) 

391 gamma = 0. 

392 

393 return (alpha, beta, gamma) 

394 

395 

396def cvec(x, y, z): 

397 return num.array([[x, y, z]], dtype=float).T 

398 

399 

400def rvec(x, y, z): 

401 return num.array([[x, y, z]], dtype=float) 

402 

403 

404def eigh_check(a): 

405 evals, evecs = num.linalg.eigh(a) 

406 assert evals[0] <= evals[1] <= evals[2] 

407 return evals, evecs 

408 

409 

410r2d = 180. / math.pi 

411d2r = 1. / r2d 

412 

413 

414def random_mt(x=None, scalar_moment=1.0, magnitude=None): 

415 

416 if magnitude is not None: 

417 scalar_moment = magnitude_to_moment(magnitude) 

418 

419 if x is None: 

420 x = num.random.random(6) 

421 

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)) 

427 

428 

429def random_m6(*args, **kwargs): 

430 return to6(random_mt(*args, **kwargs)) 

431 

432 

433def random_dc(x=None, scalar_moment=1.0, magnitude=None): 

434 if magnitude is not None: 

435 scalar_moment = magnitude_to_moment(magnitude) 

436 

437 rotmat = random_rotation(x) 

438 return scalar_moment * num.dot( 

439 rotmat, num.dot(MomentTensor._m_unrot, rotmat.T)) 

440 

441 

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]) 

446 

447 

448def as_mt(mt): 

449 ''' 

450 Convenience function to convert various objects to moment tensor object. 

451 

452 Like :py:meth:``MomentTensor.from_values``, but does not create a new 

453 :py:class:`MomentTensor` object when ``mt`` already is one. 

454 ''' 

455 

456 if isinstance(mt, MomentTensor): 

457 return mt 

458 else: 

459 return MomentTensor.from_values(mt) 

460 

461 

462class MomentTensor(Object): 

463 

464 ''' 

465 Moment tensor object 

466 

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 

473 

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). 

477 

478 .. math:: 

479 :nowrap: 

480 

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*} 

486 

487 ''' 

488 

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 

503 

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) 

512 

513 _u_evals, _u_evecs = eigh_check(_m_unrot) 

514 

515 @classmethod 

516 def random_dc(cls, x=None, scalar_moment=1.0, magnitude=None): 

517 ''' 

518 Create random oriented double-couple moment tensor 

519 

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)) 

524 

525 @classmethod 

526 def random_mt(cls, x=None, scalar_moment=1.0, magnitude=None): 

527 ''' 

528 Create random moment tensor 

529 

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)) 

536 

537 @classmethod 

538 def from_values(cls, values): 

539 ''' 

540 Alternative constructor for moment tensor objects 

541 

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. 

545 

546 The ``values`` argument is interpreted depending on shape and type as 

547 follows: 

548 

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 ''' 

556 

557 m = values_to_matrix(values) 

558 return MomentTensor(m=m) 

559 

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): 

567 

568 Object.__init__(self, init_props=False) 

569 

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) 

572 

573 strike = d2r*strike 

574 dip = d2r*dip 

575 rake = d2r*rake 

576 

577 if isinstance(m, num.matrix): 

578 m = m.A 

579 

580 if isinstance(m_up_south_east, num.matrix): 

581 m_up_south_east = m_up_south_east.A 

582 

583 if isinstance(m_east_north_up, num.matrix): 

584 m_east_north_up = m_east_north_up.A 

585 

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)) 

590 

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)) 

595 

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') 

602 

603 if moment is not None: 

604 scalar_moment = moment 

605 

606 if magnitude is not None: 

607 scalar_moment = magnitude_to_moment(magnitude) 

608 

609 rotmat1 = euler_to_matrix(dip, strike, -rake) 

610 m = scalar_moment * num.dot( 

611 rotmat1.T, 

612 num.dot(MomentTensor._m_unrot, rotmat1)) 

613 

614 self._m = m 

615 self._update() 

616 

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. 

621 

622 rotmat1 = num.dot(m_evecs, MomentTensor._u_evecs.T).T 

623 if num.linalg.det(rotmat1) < 0.: 

624 rotmat1 *= -1. 

625 

626 self._m_eigenvals = m_evals 

627 self._m_eigenvecs = m_evecs 

628 

629 self._rotmats = sorted( 

630 [rotmat1, num.dot(MomentTensor._flip_dc, rotmat1)], 

631 key=lambda m: num.abs(m.flat).tolist()) 

632 

633 @property 

634 def mnn(self): 

635 return float(self._m[0, 0]) 

636 

637 @mnn.setter 

638 def mnn(self, value): 

639 self._m[0, 0] = value 

640 self._update() 

641 

642 @property 

643 def mee(self): 

644 return float(self._m[1, 1]) 

645 

646 @mee.setter 

647 def mee(self, value): 

648 self._m[1, 1] = value 

649 self._update() 

650 

651 @property 

652 def mdd(self): 

653 return float(self._m[2, 2]) 

654 

655 @mdd.setter 

656 def mdd(self, value): 

657 self._m[2, 2] = value 

658 self._update() 

659 

660 @property 

661 def mne(self): 

662 return float(self._m[0, 1]) 

663 

664 @mne.setter 

665 def mne(self, value): 

666 self._m[0, 1] = value 

667 self._m[1, 0] = value 

668 self._update() 

669 

670 @property 

671 def mnd(self): 

672 return float(self._m[0, 2]) 

673 

674 @mnd.setter 

675 def mnd(self, value): 

676 self._m[0, 2] = value 

677 self._m[2, 0] = value 

678 self._update() 

679 

680 @property 

681 def med(self): 

682 return float(self._m[1, 2]) 

683 

684 @med.setter 

685 def med(self, value): 

686 self._m[1, 2] = value 

687 self._m[2, 1] = value 

688 self._update() 

689 

690 @property 

691 def strike1(self): 

692 return float(self.both_strike_dip_rake()[0][0]) 

693 

694 @property 

695 def dip1(self): 

696 return float(self.both_strike_dip_rake()[0][1]) 

697 

698 @property 

699 def rake1(self): 

700 return float(self.both_strike_dip_rake()[0][2]) 

701 

702 @property 

703 def strike2(self): 

704 return float(self.both_strike_dip_rake()[1][0]) 

705 

706 @property 

707 def dip2(self): 

708 return float(self.both_strike_dip_rake()[1][1]) 

709 

710 @property 

711 def rake2(self): 

712 return float(self.both_strike_dip_rake()[1][2]) 

713 

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)) 

722 

723 return results 

724 

725 def p_axis(self): 

726 ''' 

727 Get direction of p axis. 

728 ''' 

729 return (self._m_eigenvecs.T)[0] 

730 

731 def t_axis(self): 

732 ''' 

733 Get direction of t axis. 

734 ''' 

735 return (self._m_eigenvecs.T)[2] 

736 

737 def null_axis(self): 

738 ''' 

739 Get diretion of the null axis. 

740 ''' 

741 return self._m_eigenvecs.T[1] 

742 

743 def eigenvals(self): 

744 ''' 

745 Get the eigenvalues of the moment tensor in accending order. 

746 

747 :returns: ``(ep, en, et)`` 

748 ''' 

749 

750 return self._m_eigenvals 

751 

752 def eigensystem(self): 

753 ''' 

754 Get the eigenvalues and eigenvectors of the moment tensor. 

755 

756 :returns: ``(ep, en, et, vp, vn, vt)`` 

757 ''' 

758 

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 

764 

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] 

771 

772 def m(self): 

773 ''' 

774 Get plain moment tensor as 3x3 matrix. 

775 ''' 

776 return self._m.copy() 

777 

778 def m6(self): 

779 ''' 

780 Get the moment tensor as a six-element array. 

781 

782 :returns: ``(mnn, mee, mdd, mne, mnd, med)`` 

783 ''' 

784 return to6(self._m) 

785 

786 def m_up_south_east(self): 

787 ''' 

788 Get moment tensor in up-south-east convention as 3x3 matrix. 

789 

790 .. math:: 

791 :nowrap: 

792 

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 ''' 

799 

800 return num.dot( 

801 self._to_up_south_east.T, 

802 num.dot(self._m, self._to_up_south_east)) 

803 

804 def m_east_north_up(self): 

805 ''' 

806 Get moment tensor in east-north-up convention as 3x3 matrix. 

807 ''' 

808 

809 return num.dot( 

810 self._to_east_north_up.T, num.dot(self._m, self._to_east_north_up)) 

811 

812 def m6_up_south_east(self): 

813 ''' 

814 Get moment tensor in up-south-east convention as a six-element array. 

815 

816 :returns: ``(muu, mss, mee, mus, mue, mse)`` 

817 ''' 

818 return to6(self.m_up_south_east()) 

819 

820 def m6_east_north_up(self): 

821 ''' 

822 Get moment tensor in east-north-up convention as a six-element array. 

823 

824 :returns: ``(mee, mnn, muu, men, meu, mnu)`` 

825 ''' 

826 return to6(self.m_east_north_up()) 

827 

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 

836 

837 def moment_magnitude(self): 

838 ''' 

839 Get moment magnitude of moment tensor. 

840 ''' 

841 return moment_to_magnitude(self.scalar_moment()) 

842 

843 def scalar_moment(self): 

844 ''' 

845 Get the scalar moment of the moment tensor (Frobenius norm based) 

846 

847 .. math:: 

848 

849 M_0 = \\frac{1}{\\sqrt{2}}\\sqrt{\\sum_{i,j} |M_{ij}|^2} 

850 

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.) 

857 

858 @property 

859 def moment(self): 

860 return float(self.scalar_moment()) 

861 

862 @moment.setter 

863 def moment(self, value): 

864 self._m *= value / self.moment 

865 self._update() 

866 

867 @property 

868 def magnitude(self): 

869 return float(self.moment_magnitude()) 

870 

871 @magnitude.setter 

872 def magnitude(self, value): 

873 self._m *= magnitude_to_moment(value) / self.moment 

874 self._update() 

875 

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) 

887 

888 s += self.str_fault_planes() 

889 return s 

890 

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]) 

897 

898 return s 

899 

900 def deviatoric(self): 

901 ''' 

902 Get deviatoric part of moment tensor. 

903 

904 Returns a new moment tensor object with zero trace. 

905 ''' 

906 

907 m = self.m() 

908 

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 

914 

915 def standard_decomposition(self): 

916 ''' 

917 Decompose moment tensor into isotropic, DC and CLVD components. 

918 

919 Standard decomposition according to e.g. Jost and Herrmann 1989 is 

920 returned as:: 

921 

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 ''' 

930 

931 epsilon = 1e-6 

932 

933 m = self.m() 

934 

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.) 

938 

939 m_devi = m - m_iso 

940 

941 evals, evecs = eigh_check(m_devi) 

942 

943 moment_devi = num.max(num.abs(evals)) 

944 moment = moment_iso + moment_devi 

945 

946 iorder = num.argsort(num.abs(evals)) 

947 evals_sorted = evals[iorder] 

948 evecs_sorted = (evecs.T[iorder]).T 

949 

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]))) 

956 

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)) 

960 

961 m_clvd = m_devi - m_dc 

962 

963 moment_clvd = moment_devi - moment_dc 

964 

965 ratio_dc = moment_dc / moment 

966 ratio_clvd = moment_clvd / moment 

967 ratio_iso = moment_iso / moment 

968 ratio_devi = moment_devi / moment 

969 

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)] 

976 

977 def rotated(self, rot): 

978 ''' 

979 Get rotated moment tensor. 

980 

981 :param rot: ratation matrix, coordinate system is NED 

982 :returns: new :py:class:`MomentTensor` object 

983 ''' 

984 

985 rotmat = num.array(rot) 

986 return MomentTensor(m=num.dot(rotmat, num.dot(self.m(), rotmat.T))) 

987 

988 def random_rotated(self, angle_std=None, angle=None, rstate=None): 

989 ''' 

990 Get distorted MT by rotation around random axis and angle. 

991 

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 ''' 

999 

1000 assert (angle_std is None) != (angle is None), \ 

1001 'either angle or angle_std must be given' 

1002 

1003 if angle_std is not None: 

1004 rstate = rstate or num.random 

1005 angle = rstate.normal(scale=angle_std) 

1006 

1007 axis = random_axis(rstate=rstate) 

1008 rot = rotation_from_angle_and_axis(angle, axis) 

1009 return self.rotated(rot) 

1010 

1011 

1012def other_plane(strike, dip, rake): 

1013 ''' 

1014 Get the respectively other plane in the double-couple ambiguity. 

1015 ''' 

1016 

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)] 

1021 

1022 if w[0] < w[1]: 

1023 return both_sdr[1] 

1024 else: 

1025 return both_sdr[0] 

1026 

1027 

1028def dsdr(sdr1, sdr2): 

1029 s1, d1, r1 = sdr1 

1030 s2, d2, r2 = sdr2 

1031 

1032 s1 = s1 % 360. 

1033 s2 = s2 % 360. 

1034 r1 = r1 % 360. 

1035 r2 = r2 % 360. 

1036 

1037 ds = abs(s1 - s2) 

1038 ds = ds if ds <= 180. else 360. - ds 

1039 

1040 dr = abs(r1 - r2) 

1041 dr = dr if dr <= 180. else 360. - dr 

1042 

1043 dd = abs(d1 - d2) 

1044 

1045 return math.sqrt(ds**2 + dr**2 + dd**2) 

1046 

1047 

1048def order_like(sdrs, sdrs_ref): 

1049 ''' 

1050 Order strike-dip-rake pair post closely to a given reference pair. 

1051 

1052 :param sdrs: tuple, ``((strike1, dip1, rake1), (strike2, dip2, rake2))`` 

1053 :param sdrs_ref: as above but with reference pair 

1054 ''' 

1055 

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] 

1062 

1063 

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] 

1070 

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] 

1088 

1089 q[1:] /= 4.0 * q[0] 

1090 

1091 q /= math.sqrt(num.sum(q**2)) 

1092 

1093 return q 

1094 

1095 

1096_pbt2tpb = num.array(((0., 0., 1.), (1., 0., 0.), (0., 1., 0.))) 

1097 

1098 

1099def kagan_angle(mt1, mt2): 

1100 ''' 

1101 Given two moment tensors, return the Kagan angle in degrees. 

1102 

1103 After Kagan (1991) and Tape & Tape (2012). 

1104 ''' 

1105 

1106 ai = num.dot(_pbt2tpb, mt1._m_eigenvecs.T) 

1107 aj = num.dot(_pbt2tpb, mt2._m_eigenvecs.T) 

1108 

1109 u = num.dot(ai, aj.T) 

1110 

1111 tk, pk, bk = u 

1112 

1113 qk = _tpb2q(tk, pk, bk) 

1114 

1115 return 2. * r2d * math.acos(num.max(num.abs(qk))) 

1116 

1117 

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 

1123 

1124 

1125def magnitude_to_duration_gcmt(magnitudes): 

1126 ''' 

1127 Scaling relation used by Global CMT catalog for most of its events. 

1128 ''' 

1129 

1130 mom = magnitude_to_moment(magnitudes) 

1131 return (mom / 1.1e16)**(1./3.) 

1132 

1133 

1134if __name__ == '__main__': 

1135 

1136 import sys 

1137 v = list(map(float, sys.argv[1:])) 

1138 mt = MomentTensor.from_values(v) 

1139 print(mt)