Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/moment_tensor.py: 86%

493 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 15:01 +0000

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 

24import math 

25import numpy as num 

26 

27from .guts import Object, Float 

28 

29guts_prefix = 'pf' 

30 

31dynecm = 1e-7 

32 

33 

34def random_axis(rstate=None): 

35 ''' 

36 Get randomly oriented unit vector. 

37 

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 

47 

48 

49def rotation_from_angle_and_axis(angle, axis): 

50 ''' 

51 Build rotation matrix based on axis and angle. 

52 

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

58 

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) 

64 

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 

71 

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

79 

80 

81def random_rotation(x=None): 

82 ''' 

83 Get random rotation matrix. 

84 

85 A random rotation matrix, drawn from a uniform distribution in the space 

86 of rotations is returned, after Avro 1992 - "Fast random rotation 

87 matrices". 

88 

89 :param x: three (uniform random) numbers in the range [0, 1[ used as input 

90 to the distribution transformation. If ``None``, random numbers are 

91 used. Can be used to create grids of random rotations with uniform 

92 density in rotation space. 

93 ''' 

94 

95 if x is not None: 

96 x1, x2, x3 = x 

97 else: 

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

99 

100 phi = math.pi*2.0*x1 

101 

102 zrot = num.array([ 

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

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

105 [0., 0., 1.]]) 

106 

107 lam = math.pi*2.0*x2 

108 

109 v = num.array([[ 

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

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

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

113 

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

115 return -num.dot(house, zrot) 

116 

117 

118def rand(mi, ma): 

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

120 

121 

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 

126 

127 

128def random_strike_dip_rake( 

129 strikemin=0., strikemax=360., 

130 dipmin=0., dipmax=90., 

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

132 

133 ''' 

134 Get random strike, dip, rake triplet. 

135 

136 .. note:: 

137 

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

142 

143 strike = rand(strikemin, strikemax) 

144 dip = randdip(dipmin, dipmax) 

145 rake = rand(rakemin, rakemax) 

146 

147 return strike, dip, rake 

148 

149 

150def to6(m): 

151 ''' 

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

153 

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

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

156 ''' 

157 

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

159 

160 

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

165 

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) 

169 

170 

171def values_to_matrix(values): 

172 ''' 

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

174 

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. 

177 

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

179 follows: 

180 

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

188 

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

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

191 

192 if isinstance(values, MomentTensor): 

193 return values.m() 

194 

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

200 

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

207 

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

209 return symmat6(*values) 

210 

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 

217 

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

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

220 

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

222 

223 

224def moment_to_magnitude(moment): 

225 ''' 

226 Convert scalar moment to moment magnitude Mw. 

227 

228 :param moment: scalar moment [Nm] 

229 :returns: moment magnitude Mw 

230 

231 Moment magnitude is defined as 

232 

233 .. math:: 

234 

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

236 

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

238 

239 .. note:: 

240 

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

242 10.7 is from [Hanks and Kanamori 1979]. 

243 ''' 

244 

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

246 

247 

248def magnitude_to_moment(magnitude): 

249 ''' 

250 Convert moment magnitude Mw to scalar moment. 

251 

252 :param magnitude: moment magnitude 

253 :returns: scalar moment [Nm] 

254 

255 See :py:func:`moment_to_magnitude`. 

256 ''' 

257 

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

259 

260 

261magnitude_1Nm = moment_to_magnitude(1.0) 

262 

263 

264def euler_to_matrix(alpha, beta, gamma): 

265 ''' 

266 Given euler angle triplet, create rotation matrix 

267 

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. 

271 

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] 

276 

277 Usage for moment tensors:: 

278 

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

282 

283 ''' 

284 

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) 

291 

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 

296 

297 

298def matrix_to_euler(rotmat): 

299 ''' 

300 Get Eulerian angle triplet from rotation matrix. 

301 ''' 

302 

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)[0, 0]) 

312 if cos_alpha > 1.: 

313 cos_alpha = 1. 

314 

315 if cos_alpha < -1.: 

316 cos_alpha = -1. 

317 

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

321 

322 return unique_euler(alpha, beta, gamma) 

323 

324 

325def unique_euler(alpha, beta, gamma): 

326 ''' 

327 Uniquify Eulerian angle triplet. 

328 

329 Put Eulerian angle triplet into ranges compatible with 

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

331 

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

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

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

335 

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. 

338 

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

340 ``[0,pi)``. 

341 ''' 

342 

343 pi = math.pi 

344 

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

346 

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 

358 

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 

362 

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) 

366 

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 

370 

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

372 beta = pi 

373 

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

375 beta = 0. 

376 

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

378 beta = 0. 

379 

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 

386 

387 if alpha < 1e-7: 

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

389 gamma = 0. 

390 

391 return (alpha, beta, gamma) 

392 

393 

394def cvec(x, y, z): 

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

396 

397 

398def rvec(x, y, z): 

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

400 

401 

402def eigh_check(a): 

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

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

405 return evals, evecs 

406 

407 

408r2d = 180. / math.pi 

409d2r = 1. / r2d 

410 

411 

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

413 

414 if magnitude is not None: 

415 scalar_moment = magnitude_to_moment(magnitude) 

416 

417 if x is None: 

418 x = num.random.random(6) 

419 

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

425 

426 

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

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

429 

430 

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

432 if magnitude is not None: 

433 scalar_moment = magnitude_to_moment(magnitude) 

434 

435 rotmat = random_rotation(x) 

436 return scalar_moment * num.dot( 

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

438 

439 

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

444 

445 

446def as_mt(mt): 

447 ''' 

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

449 

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

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

452 ''' 

453 

454 if isinstance(mt, MomentTensor): 

455 return mt 

456 else: 

457 return MomentTensor.from_values(mt) 

458 

459 

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) 

470 

471 

472def ned_to_rta(ned): 

473 ''' 

474 Convert north-east-down coordinate vectors to radial-takeoff-azimuth. 

475 

476 Output coordinate system has coordinates radial, takeoff angle [deg] 

477 (downward is zero), and azimuth [deg] (northward is zero). 

478 

479 :param ned: 

480 Coordinate vector or array of coordinate vectors (north, east, down). 

481 :type ned: 

482 :py:class:`numpy.ndarray` of shape ``(N, 3)`` or ``(3,)`` 

483 

484 :returns: 

485 Coordinate vector or array of coordinate vectors (radial, takeoff, 

486 azimuth) 

487 :rtype: 

488 :py:class:`numpy.ndarray` of shape ``(N, 3)`` or ``(3,)`` 

489 ''' 

490 if ned.ndim == 1: 

491 return ned_to_rta(ned[num.newaxis, :])[0] 

492 

493 rta = num.empty_like(ned) 

494 rta[:, 0] = num.sqrt(num.sum(ned**2, axis=1)) 

495 rta[:, 1] = num.arccos(ned[:, 2]) * r2d 

496 rta[:, 2] = num.arctan2(ned[:, 1], ned[:, 0]) * r2d 

497 return rta 

498 

499 

500class MomentTensor(Object): 

501 

502 ''' 

503 Moment tensor object 

504 

505 :param m: NumPy array in north-east-down convention 

506 :param m_up_south_east: NumPy array in up-south-east convention 

507 :param m_east_north_up: NumPy array in east-north-up convention 

508 :param strike,dip,rake: fault plane angles in [degrees] 

509 :param p_axis,t_axis: initialize double-couple from p and t axes 

510 :param scalar_moment: scalar moment in [Nm] 

511 :param magnitude: moment magnitude Mw 

512 

513 Global CMT catalog moment tensors use the up-south-east (USE) coordinate 

514 system convention with :math:`r` (up), :math:`\\theta` (south), and 

515 :math:`\\phi` (east). 

516 

517 .. math:: 

518 :nowrap: 

519 

520 \\begin{align*} 

521 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\ 

522 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\ 

523 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne} 

524 \\end{align*} 

525 

526 ''' 

527 

528 mnn__ = Float.T(default=0.0) 

529 mee__ = Float.T(default=0.0) 

530 mdd__ = Float.T(default=0.0) 

531 mne__ = Float.T(default=0.0) 

532 mnd__ = Float.T(default=-1.0) 

533 med__ = Float.T(default=0.0) 

534 strike1__ = Float.T(default=None, optional=True) # read-only 

535 dip1__ = Float.T(default=None, optional=True) # read-only 

536 rake1__ = Float.T(default=None, optional=True) # read-only 

537 strike2__ = Float.T(default=None, optional=True) # read-only 

538 dip2__ = Float.T(default=None, optional=True) # read-only 

539 rake2__ = Float.T(default=None, optional=True) # read-only 

540 moment__ = Float.T(default=None, optional=True) # read-only 

541 magnitude__ = Float.T(default=None, optional=True) # read-only 

542 

543 _flip_dc = num.array( 

544 [[0., 0., -1.], [0., -1., 0.], [-1., 0., 0.]], dtype=float) 

545 _to_up_south_east = num.array( 

546 [[0., 0., -1.], [-1., 0., 0.], [0., 1., 0.]], dtype=float).T 

547 _to_east_north_up = num.array( 

548 [[0., 1., 0.], [1., 0., 0.], [0., 0., -1.]], dtype=float) 

549 _m_unrot = num.array( 

550 [[0., 0., -1.], [0., 0., 0.], [-1., 0., 0.]], dtype=float) 

551 

552 _u_evals, _u_evecs = eigh_check(_m_unrot) 

553 

554 @classmethod 

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

556 ''' 

557 Create random oriented double-couple moment tensor 

558 

559 The rotations used are uniformly distributed in the space of rotations. 

560 ''' 

561 return MomentTensor( 

562 m=random_dc(x=x, scalar_moment=scalar_moment, magnitude=magnitude)) 

563 

564 @classmethod 

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

566 ''' 

567 Create random moment tensor 

568 

569 Moment tensors produced by this function appear uniformly distributed 

570 when shown in a Hudson's diagram. The rotations used are uniformly 

571 distributed in the space of rotations. 

572 ''' 

573 return MomentTensor( 

574 m=random_mt(x=x, scalar_moment=scalar_moment, magnitude=magnitude)) 

575 

576 @classmethod 

577 def from_values(cls, values): 

578 ''' 

579 Alternative constructor for moment tensor objects 

580 

581 This constructor takes a :py:class:`MomentTensor` object, a tuple, list 

582 or NumPy array with 3x3 or 3, 4, 6, or 7 elements to build a Moment 

583 tensor object. 

584 

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

586 follows: 

587 

588 * ``(strike, dip, rake)`` 

589 * ``(strike, dip, rake, magnitude)`` 

590 * ``(mnn, mee, mdd, mne, mnd, med)`` 

591 * ``(mnn, mee, mdd, mne, mnd, med, magnitude)`` 

592 * ``((mnn, mne, mnd), (mne, mee, med), (mnd, med, mdd))`` 

593 * :py:class:`MomentTensor` object 

594 ''' 

595 

596 m = values_to_matrix(values) 

597 return MomentTensor(m=m) 

598 

599 def __init__( 

600 self, m=None, m_up_south_east=None, m_east_north_up=None, 

601 strike=0., dip=0., rake=0., scalar_moment=1., 

602 mnn=None, mee=None, mdd=None, mne=None, mnd=None, med=None, 

603 strike1=None, dip1=None, rake1=None, 

604 strike2=None, dip2=None, rake2=None, 

605 p_axis=None, t_axis=None, 

606 magnitude=None, moment=None): 

607 

608 Object.__init__(self, init_props=False) 

609 

610 if any(mxx is not None for mxx in (mnn, mee, mdd, mne, mnd, med)): 

611 m = symmat6(mnn, mee, mdd, mne, mnd, med) 

612 

613 if p_axis is not None and t_axis is not None: 

614 strike, dip, rake = pt_axes_to_strike_dip_rake(p_axis, t_axis) 

615 

616 strike = d2r*strike 

617 dip = d2r*dip 

618 rake = d2r*rake 

619 

620 if isinstance(m, num.matrix): 

621 m = m.A 

622 

623 if isinstance(m_up_south_east, num.matrix): 

624 m_up_south_east = m_up_south_east.A 

625 

626 if isinstance(m_east_north_up, num.matrix): 

627 m_east_north_up = m_east_north_up.A 

628 

629 if m_up_south_east is not None: 

630 m = num.dot( 

631 self._to_up_south_east, 

632 num.dot(m_up_south_east, self._to_up_south_east.T)) 

633 

634 if m_east_north_up is not None: 

635 m = num.dot( 

636 self._to_east_north_up, 

637 num.dot(m_east_north_up, self._to_east_north_up.T)) 

638 

639 if m is None: 

640 if any(x is not None for x in ( 

641 strike1, dip1, rake1, strike2, dip2, rake2)): 

642 raise Exception( 

643 'strike1, dip1, rake1, strike2, dip2, rake2 are read-only ' 

644 'properties') 

645 

646 if moment is not None: 

647 scalar_moment = moment 

648 

649 if magnitude is not None: 

650 scalar_moment = magnitude_to_moment(magnitude) 

651 

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

653 m = scalar_moment * num.dot( 

654 rotmat1.T, 

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

656 

657 self._m = m 

658 self._update() 

659 

660 def _update(self): 

661 m_evals, m_evecs = eigh_check(self._m) 

662 if num.linalg.det(m_evecs) < 0.: 

663 m_evecs *= -1. 

664 

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

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

667 rotmat1 *= -1. 

668 

669 self._m_eigenvals = m_evals 

670 self._m_eigenvecs = m_evecs 

671 

672 self._rotmats = sorted( 

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

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

675 

676 @property 

677 def mnn(self): 

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

679 

680 @mnn.setter 

681 def mnn(self, value): 

682 self._m[0, 0] = value 

683 self._update() 

684 

685 @property 

686 def mee(self): 

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

688 

689 @mee.setter 

690 def mee(self, value): 

691 self._m[1, 1] = value 

692 self._update() 

693 

694 @property 

695 def mdd(self): 

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

697 

698 @mdd.setter 

699 def mdd(self, value): 

700 self._m[2, 2] = value 

701 self._update() 

702 

703 @property 

704 def mne(self): 

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

706 

707 @mne.setter 

708 def mne(self, value): 

709 self._m[0, 1] = value 

710 self._m[1, 0] = value 

711 self._update() 

712 

713 @property 

714 def mnd(self): 

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

716 

717 @mnd.setter 

718 def mnd(self, value): 

719 self._m[0, 2] = value 

720 self._m[2, 0] = value 

721 self._update() 

722 

723 @property 

724 def med(self): 

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

726 

727 @med.setter 

728 def med(self, value): 

729 self._m[1, 2] = value 

730 self._m[2, 1] = value 

731 self._update() 

732 

733 @property 

734 def strike1(self): 

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

736 

737 @property 

738 def dip1(self): 

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

740 

741 @property 

742 def rake1(self): 

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

744 

745 @property 

746 def strike2(self): 

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

748 

749 @property 

750 def dip2(self): 

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

752 

753 @property 

754 def rake2(self): 

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

756 

757 def both_strike_dip_rake(self): 

758 ''' 

759 Get both possible (strike,dip,rake) triplets. 

760 ''' 

761 results = [] 

762 for rotmat in self._rotmats: 

763 alpha, beta, gamma = [r2d*x for x in matrix_to_euler(rotmat)] 

764 results.append((beta, alpha, -gamma)) 

765 

766 return results 

767 

768 def p_axis(self): 

769 ''' 

770 Get direction of p axis. 

771 

772 Use :py:func:`ned_to_rta` to get takeoff angle and azimuth of the 

773 returned vectors. 

774 

775 :returns: 

776 Direction of P (pressure) axis as a (north, east, down) vector. 

777 :rtype: 

778 :py:class:`numpy.ndarray` of shape ``(3,)`` 

779 ''' 

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

781 

782 def t_axis(self): 

783 ''' 

784 Get direction of t axis. 

785 

786 Use :py:func:`ned_to_rta` to get takeoff angle and azimuth of the 

787 returned vectors. 

788 

789 :returns: 

790 Direction of T (tension) axis as a (north, east, down) vector. 

791 :rtype: 

792 :py:class:`numpy.ndarray` of shape ``(3,)`` 

793 ''' 

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

795 

796 def null_axis(self): 

797 ''' 

798 Get direction of the null axis. 

799 

800 Use :py:func:`ned_to_rta` to get takeoff angle and azimuth of the 

801 returned vectors. 

802 

803 :returns: 

804 Direction of null (B) axis as a (north, east, down) vector. 

805 :rtype: 

806 :py:class:`numpy.ndarray` of shape ``(3,)`` 

807 ''' 

808 return self._m_eigenvecs.T[1] 

809 

810 def eigenvals(self): 

811 ''' 

812 Get the eigenvalues of the moment tensor in ascending order. 

813 

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

815 ''' 

816 

817 return self._m_eigenvals 

818 

819 def eigensystem(self): 

820 ''' 

821 Get the eigenvalues and eigenvectors of the moment tensor. 

822 

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

824 ''' 

825 

826 vp = self.p_axis().flatten() 

827 vn = self.null_axis().flatten() 

828 vt = self.t_axis().flatten() 

829 ep, en, et = self._m_eigenvals 

830 return ep, en, et, vp, vn, vt 

831 

832 def both_slip_vectors(self): 

833 ''' 

834 Get both possible slip directions. 

835 ''' 

836 return [ 

837 num.dot(rotmat.T, cvec(1., 0., 0.)) for rotmat in self._rotmats] 

838 

839 def m(self): 

840 ''' 

841 Get plain moment tensor as 3x3 matrix. 

842 ''' 

843 return self._m.copy() 

844 

845 def m6(self): 

846 ''' 

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

848 

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

850 ''' 

851 return to6(self._m) 

852 

853 def m_up_south_east(self): 

854 ''' 

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

856 

857 .. math:: 

858 :nowrap: 

859 

860 \\begin{align*} 

861 M_{rr} &= M_{dd}, & M_{ r\\theta} &= M_{nd},\\\\ 

862 M_{\\theta\\theta} &= M_{ nn}, & M_{r\\phi} &= -M_{ed},\\\\ 

863 M_{\\phi\\phi} &= M_{ee}, & M_{\\theta\\phi} &= -M_{ne} 

864 \\end{align*} 

865 ''' 

866 

867 return num.dot( 

868 self._to_up_south_east.T, 

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

870 

871 def m_east_north_up(self): 

872 ''' 

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

874 ''' 

875 

876 return num.dot( 

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

878 

879 def m6_up_south_east(self): 

880 ''' 

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

882 

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

884 ''' 

885 return to6(self.m_up_south_east()) 

886 

887 def m6_east_north_up(self): 

888 ''' 

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

890 

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

892 ''' 

893 return to6(self.m_east_north_up()) 

894 

895 def m_plain_double_couple(self): 

896 ''' 

897 Get plain double couple with same scalar moment as moment tensor. 

898 ''' 

899 rotmat1 = self._rotmats[0] 

900 m = self.scalar_moment() * num.dot( 

901 rotmat1.T, num.dot(MomentTensor._m_unrot, rotmat1)) 

902 return m 

903 

904 def moment_magnitude(self): 

905 ''' 

906 Get moment magnitude of moment tensor. 

907 ''' 

908 return moment_to_magnitude(self.scalar_moment()) 

909 

910 def scalar_moment(self): 

911 ''' 

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

913 

914 .. math:: 

915 

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

917 

918 The scalar moment is calculated based on the Euclidean (Frobenius) norm 

919 (Silver and Jordan, 1982). The scalar moment returned by this function 

920 differs from the standard decomposition based definition of the scalar 

921 moment for non-double-couple moment tensors. 

922 ''' 

923 return num.linalg.norm(self._m_eigenvals) / math.sqrt(2.) 

924 

925 @property 

926 def moment(self): 

927 return float(self.scalar_moment()) 

928 

929 @moment.setter 

930 def moment(self, value): 

931 self._m *= value / self.moment 

932 self._update() 

933 

934 @property 

935 def magnitude(self): 

936 return float(self.moment_magnitude()) 

937 

938 @magnitude.setter 

939 def magnitude(self, value): 

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

941 self._update() 

942 

943 def __str__(self): 

944 mexp = pow(10, math.ceil(num.log10(num.max(num.abs(self._m))))) 

945 m = self._m / mexp 

946 s = '''Scalar Moment [Nm]: M0 = %g (Mw = %3.1f) 

947Moment Tensor [Nm]: Mnn = %6.3f, Mee = %6.3f, Mdd = %6.3f, 

948 Mne = %6.3f, Mnd = %6.3f, Med = %6.3f [ x %g ] 

949''' % ( 

950 self.scalar_moment(), 

951 self.moment_magnitude(), 

952 m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2], 

953 mexp) 

954 

955 s += self.str_fault_planes() 

956 return s 

957 

958 def str_fault_planes(self): 

959 s = '' 

960 for i, sdr in enumerate(self.both_strike_dip_rake()): 

961 s += 'Fault plane %i [deg]: ' \ 

962 'strike = %3.0f, dip = %3.0f, slip-rake = %4.0f\n' \ 

963 % (i+1, sdr[0], sdr[1], sdr[2]) 

964 

965 return s 

966 

967 def deviatoric(self): 

968 ''' 

969 Get deviatoric part of moment tensor. 

970 

971 Returns a new moment tensor object with zero trace. 

972 ''' 

973 

974 m = self.m() 

975 

976 trace_m = num.trace(m) 

977 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.]) 

978 m_devi = m - m_iso 

979 mt = MomentTensor(m=m_devi) 

980 return mt 

981 

982 def standard_decomposition(self): 

983 ''' 

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

985 

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

987 returned as:: 

988 

989 [ 

990 (moment_iso, ratio_iso, m_iso), 

991 (moment_dc, ratio_dc, m_dc), 

992 (moment_clvd, ratio_clvd, m_clvd), 

993 (moment_devi, ratio_devi, m_devi), 

994 (moment, 1.0, m) 

995 ] 

996 ''' 

997 

998 epsilon = 1e-6 

999 

1000 m = self.m() 

1001 

1002 trace_m = num.trace(m) 

1003 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.]) 

1004 moment_iso = abs(trace_m / 3.) 

1005 

1006 m_devi = m - m_iso 

1007 

1008 evals, evecs = eigh_check(m_devi) 

1009 

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

1011 moment = moment_iso + moment_devi 

1012 

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

1014 evals_sorted = evals[iorder] 

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

1016 

1017 if moment_devi < epsilon * moment_iso: 

1018 signed_moment_dc = 0. 

1019 else: 

1020 assert -epsilon <= -evals_sorted[0] / evals_sorted[2] <= 0.5 

1021 signed_moment_dc = evals_sorted[2] * (1.0 + 2.0 * ( 

1022 min(0.0, evals_sorted[0] / evals_sorted[2]))) 

1023 

1024 moment_dc = abs(signed_moment_dc) 

1025 m_dc_es = num.dot(signed_moment_dc, num.diag([0., -1.0, 1.0])) 

1026 m_dc = num.dot(evecs_sorted, num.dot(m_dc_es, evecs_sorted.T)) 

1027 

1028 m_clvd = m_devi - m_dc 

1029 

1030 moment_clvd = moment_devi - moment_dc 

1031 

1032 ratio_dc = moment_dc / moment 

1033 ratio_clvd = moment_clvd / moment 

1034 ratio_iso = moment_iso / moment 

1035 ratio_devi = moment_devi / moment 

1036 

1037 return [ 

1038 (moment_iso, ratio_iso, m_iso), 

1039 (moment_dc, ratio_dc, m_dc), 

1040 (moment_clvd, ratio_clvd, m_clvd), 

1041 (moment_devi, ratio_devi, m_devi), 

1042 (moment, 1.0, m)] 

1043 

1044 def rotated(self, rot): 

1045 ''' 

1046 Get rotated moment tensor. 

1047 

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

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

1050 ''' 

1051 

1052 rotmat = num.array(rot) 

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

1054 

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

1056 ''' 

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

1058 

1059 :param angle_std: angles are drawn from a normal distribution with 

1060 zero mean and given standard deviation [degrees] 

1061 :param angle: set angle [degrees], only axis will be random 

1062 :param rstate: :py:class:`numpy.random.RandomState` object, can be 

1063 used to create reproducible pseudo-random sequences 

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

1065 ''' 

1066 

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

1068 'either angle or angle_std must be given' 

1069 

1070 if angle_std is not None: 

1071 rstate = rstate or num.random 

1072 angle = rstate.normal(scale=angle_std) 

1073 

1074 axis = random_axis(rstate=rstate) 

1075 rot = rotation_from_angle_and_axis(angle, axis) 

1076 return self.rotated(rot) 

1077 

1078 

1079def other_plane(strike, dip, rake): 

1080 ''' 

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

1082 ''' 

1083 

1084 mt = MomentTensor(strike=strike, dip=dip, rake=rake) 

1085 both_sdr = mt.both_strike_dip_rake() 

1086 w = [sum([abs(x-y) for x, y in zip(both_sdr[i], (strike, dip, rake))]) 

1087 for i in (0, 1)] 

1088 

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

1090 return both_sdr[1] 

1091 else: 

1092 return both_sdr[0] 

1093 

1094 

1095def dsdr(sdr1, sdr2): 

1096 s1, d1, r1 = sdr1 

1097 s2, d2, r2 = sdr2 

1098 

1099 s1 = s1 % 360. 

1100 s2 = s2 % 360. 

1101 r1 = r1 % 360. 

1102 r2 = r2 % 360. 

1103 

1104 ds = abs(s1 - s2) 

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

1106 

1107 dr = abs(r1 - r2) 

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

1109 

1110 dd = abs(d1 - d2) 

1111 

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

1113 

1114 

1115def order_like(sdrs, sdrs_ref): 

1116 ''' 

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

1118 

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

1120 :param sdrs_ref: as above but with reference pair 

1121 ''' 

1122 

1123 d1 = min(dsdr(sdrs[0], sdrs_ref[0]), dsdr(sdrs[1], sdrs_ref[1])) 

1124 d2 = min(dsdr(sdrs[0], sdrs_ref[1]), dsdr(sdrs[1], sdrs_ref[0])) 

1125 if d1 < d2: 

1126 return sdrs 

1127 else: 

1128 return sdrs[::-1] 

1129 

1130 

1131def _tpb2q(t, p, b): 

1132 eps = 0.001 

1133 tqw = 1. + t[0] + p[1] + b[2] 

1134 tqx = 1. + t[0] - p[1] - b[2] 

1135 tqy = 1. - t[0] + p[1] - b[2] 

1136 tqz = 1. - t[0] - p[1] + b[2] 

1137 

1138 q = num.zeros(4) 

1139 if tqw > eps: 

1140 q[0] = 0.5 * math.sqrt(tqw) 

1141 q[1:] = p[2] - b[1], b[0] - t[2], t[1] - p[0] 

1142 elif tqx > eps: 

1143 q[0] = 0.5 * math.sqrt(tqx) 

1144 q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2] 

1145 elif tqy > eps: 

1146 q[0] = 0.5 * math.sqrt(tqy) 

1147 q[1:] = b[0] - t[2], p[0] + t[1], b[1] + p[2] 

1148 elif tqz > eps: 

1149 q[0] = 0.5 * math.sqrt(tqz) 

1150 q[1:] = t[1] - p[0], b[0] + t[2], b[1] + p[2] 

1151 else: 

1152 raise Exception('should not reach this line, check theory!') 

1153 # q[0] = max(0.5 * math.sqrt(tqx), eps) 

1154 # q[1:] = p[2] - b[1], p[0] + t[1], b[0] + t[2] 

1155 

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

1157 

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

1159 

1160 return q 

1161 

1162 

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

1164 

1165 

1166def kagan_angle(mt1, mt2): 

1167 ''' 

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

1169 

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

1171 ''' 

1172 

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

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

1175 

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

1177 

1178 tk, pk, bk = u 

1179 

1180 qk = _tpb2q(tk, pk, bk) 

1181 

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

1183 

1184 

1185def rand_to_gutenberg_richter(rand, b_value, magnitude_min): 

1186 ''' 

1187 Draw magnitude from Gutenberg Richter distribution. 

1188 ''' 

1189 return magnitude_min + num.log10(1.-rand) / -b_value 

1190 

1191 

1192def magnitude_to_duration_gcmt(magnitudes): 

1193 ''' 

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

1195 ''' 

1196 

1197 mom = magnitude_to_moment(magnitudes) 

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

1199 

1200 

1201if __name__ == '__main__': 

1202 

1203 import sys 

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

1205 mt = MomentTensor.from_values(v) 

1206 print(mt)