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

2023 statements  

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

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Classical seismic ray theory for layered earth models (*layer cake* models). 

8 

9This module can be used to e.g. calculate arrival times, ray paths, reflection 

10and transmission coefficients, take-off and incidence angles and geometrical 

11spreading factors for arbitrary seismic phases. Computations are done for a 

12spherical earth, even though the module name may suggests something flat. 

13 

14The main classes defined in this module are: 

15 

16* :py:class:`Material` - Defines an isotropic elastic material. 

17* :py:class:`PhaseDef` - Defines a seismic phase arrival / wave propagation 

18 history. 

19* :py:class:`Leg` - Continuous propagation in a :py:class:`PhaseDef`. 

20* :py:class:`Knee` - Conversion/reflection in a :py:class:`PhaseDef`. 

21* :py:class:`LayeredModel` - Representation of a layer cake model. 

22* :py:class:`Layer` - A layer in a :py:class:`LayeredModel`. 

23 

24 * :py:class:`HomogeneousLayer` - A homogeneous :py:class:`Layer`. 

25 * :py:class:`GradientLayer` - A gradient :py:class:`Layer`. 

26 

27* :py:class:`Discontinuity` - A discontinuity in a :py:class:`LayeredModel`. 

28 

29 * :py:class:`Interface` - A :py:class:`Discontinuity` between two 

30 :py:class:`Layer` instances. 

31 * :py:class:`Surface` - The surface :py:class:`Discontinuity` on top of 

32 a :py:class:`LayeredModel`. 

33 

34* :py:class:`RayPath` - A fan of rays running through a common sequence of 

35 layers / interfaces. 

36* :py:class:`Ray` - A specific ray with a specific (ray parameter, distance, 

37 arrival time) choice. 

38* :py:class:`RayElement` - An element of a :py:class:`RayPath`. 

39 

40 * :py:class:`Straight` - A ray segment representing propagation through 

41 one :py:class:`Layer`. 

42 * :py:class:`Kink` - An interaction of a ray with a 

43 :py:class:`Discontinuity`. 

44''' 

45 

46from functools import reduce 

47 

48import os 

49import logging 

50import copy 

51import math 

52import cmath 

53import operator 

54try: 

55 from StringIO import StringIO 

56except ImportError: 

57 from io import StringIO 

58 

59import glob 

60import numpy as num 

61from scipy.optimize import bisect, brentq 

62 

63from . import util, config 

64 

65 

66logger = logging.getLogger('cake') 

67 

68ZEPS = 0.01 

69 

70P = 1 

71''' 

72Constant indicating P wave propagation. 

73''' 

74 

75S = 2 

76''' 

77Constant indicating S wave propagation 

78''' 

79 

80DOWN = 4 

81''' 

82Constant indicating downward direction. 

83''' 

84 

85UP = -4 

86''' 

87Constant indicating upward direction. 

88''' 

89 

90DEFAULT_BURGERS = (0., 0., 1.) 

91 

92earthradius = config.config().earthradius 

93 

94r2d = 180./math.pi 

95d2r = 1./r2d 

96km = 1000. 

97d2m = d2r*earthradius 

98m2d = 1./d2m 

99sprad2spm = 1.0/(r2d*d2m) 

100sprad2spkm = 1.0/(r2d*d2m/km) 

101spm2sprad = 1.0/sprad2spm 

102spkm2sprad = 1.0/sprad2spkm 

103 

104 

105class CakeError(Exception): 

106 ''' 

107 Base class for exceptions raised in :py:mod:`pyrocko.cake`. 

108 ''' 

109 pass 

110 

111 

112class InvalidArguments(CakeError): 

113 ''' 

114 Invalid arguments. 

115 ''' 

116 pass 

117 

118 

119class Material(object): 

120 ''' 

121 Isotropic elastic material. 

122 

123 :param vp: P-wave velocity [m/s] 

124 :param vs: S-wave velocity [m/s] 

125 :param rho: density [kg/m^3] 

126 :param qp: P-wave attenuation Qp 

127 :param qs: S-wave attenuation Qs 

128 :param poisson: Poisson ratio 

129 :param lame: tuple with Lame parameter `lambda` and `shear modulus` [Pa] 

130 :param qk: bulk attenuation Qk 

131 :param qmu: shear attenuation Qmu 

132 

133 :param burgers: Burgers rheology paramerters as `tuple`. 

134 `transient viscosity` [Pa], <= 0 means infinite value, 

135 `steady-state viscosity` [Pa] and `alpha`, the ratio between the 

136 effective and unreleaxed shear modulus, mu1/(mu1 + mu2). 

137 :type burgers: tuple 

138 

139 If no velocities and no lame parameters are given, standard crustal values 

140 of vp = 5800 m/s and vs = 3200 m/s are used. If no Q values are given, 

141 standard crustal values of qp = 1456 and qs = 600 are used. If no Burgers 

142 material parameters are given, transient and steady-state viscosities are 

143 0 and alpha=1. 

144 

145 Everything is in SI units (m/s, Pa, kg/m^3) unless explicitly stated. 

146 

147 The main material properties are considered independant and are accessible 

148 as attributes (it is allowed to assign to these): 

149 

150 .. py:attribute:: vp, vs, rho, qp, qs 

151 

152 Other material properties are considered dependant and can be queried by 

153 instance methods. 

154 ''' 

155 

156 def __init__( 

157 self, vp=None, vs=None, rho=2600., qp=None, qs=None, poisson=None, 

158 lame=None, qk=None, qmu=None, burgers=None): 

159 

160 parstore_float(locals(), self, 'vp', 'vs', 'rho', 'qp', 'qs') 

161 

162 if vp is not None and vs is not None: 

163 if poisson is not None or lame is not None: 

164 raise InvalidArguments( 

165 'If vp and vs are given, poisson ratio and lame paramters ' 

166 'should not be given.') 

167 

168 elif vp is None and vs is None and lame is None: 

169 self.vp = 5800. 

170 if poisson is None: 

171 poisson = 0.25 

172 self.vs = self.vp / math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson)) 

173 

174 elif vp is None and vs is None and lame is not None: 

175 if poisson is not None: 

176 raise InvalidArguments( 

177 'Poisson ratio should not be given, when lame parameters ' 

178 'are given.') 

179 

180 lam, mu = float(lame[0]), float(lame[1]) 

181 self.vp = math.sqrt((lam + 2.0*mu)/rho) 

182 self.vs = math.sqrt(mu/rho) 

183 

184 elif vp is not None and vs is None: 

185 if poisson is None: 

186 poisson = 0.25 

187 

188 if lame is not None: 

189 raise InvalidArguments( 

190 'If vp is given, Lame parameters should not be given.') 

191 

192 poisson = float(poisson) 

193 self.vs = vp / math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson)) 

194 

195 elif vp is None and vs is not None: 

196 if poisson is None: 

197 poisson = 0.25 

198 if lame is not None: 

199 raise InvalidArguments( 

200 'If vs is given, Lame parameters should not be given.') 

201 

202 poisson = float(poisson) 

203 self.vp = vs * math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson)) 

204 

205 else: 

206 raise InvalidArguments( 

207 'Invalid combination of input parameters in material ' 

208 'definition.') 

209 

210 if qp is not None or qs is not None: 

211 if not (qk is None and qmu is None): 

212 raise InvalidArguments( 

213 'If qp or qs are given, qk and qmu should not be given.') 

214 

215 if qp is None: 

216 if self.vs != 0.0: 

217 s = (4.0/3.0)*(self.vs/self.vp)**2 

218 self.qp = self.qs / s 

219 else: 

220 self.qp = 1456. 

221 

222 if qs is None: 

223 if self.vs != 0.0: 

224 s = (4.0/3.0)*(self.vs/self.vp)**2 

225 self.qs = self.qp * s 

226 else: 

227 self.vs = 600. 

228 

229 elif qp is None and qs is None and qk is None and qmu is None: 

230 if self.vs == 0.: 

231 self.qs = 0. 

232 self.qp = 5782e4 

233 else: 

234 self.qs = 600. 

235 s = (4.0/3.0)*(self.vs/self.vp)**2 

236 self.qp = self.qs/s 

237 

238 elif qp is None and qs is None and qk is not None and qmu is not None: 

239 s = (4.0/3.0)*(self.vs/self.vp)**2 

240 if qmu == 0. and self.vs == 0.: 

241 self.qp = qk 

242 else: 

243 if num.isinf(qk): 

244 self.qp = qmu/s 

245 else: 

246 self.qp = 1.0 / (s/qmu + (1.0-s)/qk) 

247 self.qs = qmu 

248 else: 

249 raise InvalidArguments( 

250 'Invalid combination of input parameters in material ' 

251 'definition.') 

252 

253 if burgers is None: 

254 burgers = DEFAULT_BURGERS 

255 

256 self.burger_eta1 = burgers[0] 

257 self.burger_eta2 = burgers[1] 

258 self.burger_valpha = burgers[2] 

259 

260 def astuple(self): 

261 ''' 

262 Get independant material properties as a tuple. 

263 

264 Returns a tuple with ``(vp, vs, rho, qp, qs)``. 

265 ''' 

266 return self.vp, self.vs, self.rho, self.qp, self.qs 

267 

268 def __eq__(self, other): 

269 return self.astuple() == other.astuple() 

270 

271 def lame(self): 

272 ''' 

273 Get Lame's parameter lambda and shear modulus. 

274 ''' 

275 mu = self.vs**2 * self.rho 

276 lam = self.vp**2 * self.rho - 2.0*mu 

277 return lam, mu 

278 

279 def lame_lambda(self): 

280 ''' 

281 Get Lame's parameter lambda. 

282 

283 Returned units are [Pa]. 

284 ''' 

285 lam, _ = self.lame() 

286 return lam 

287 

288 def shear_modulus(self): 

289 ''' 

290 Get shear modulus. 

291 

292 Returned units are [Pa]. 

293 ''' 

294 return self.vs**2 * self.rho 

295 

296 def poisson(self): 

297 ''' 

298 Get Poisson's ratio. 

299 ''' 

300 lam, mu = self.lame() 

301 return lam / (2.0*(lam+mu)) 

302 

303 def bulk(self): 

304 ''' 

305 Get bulk modulus. 

306 ''' 

307 lam, mu = self.lame() 

308 return lam + 2.0*mu/3.0 

309 

310 def youngs(self): 

311 ''' 

312 Get Young's modulus. 

313 ''' 

314 lam, mu = self.lame() 

315 return mu * (3.0*lam + 2.0*mu) / (lam+mu) 

316 

317 def vp_vs_ratio(self): 

318 ''' 

319 Get vp/vs ratio. 

320 ''' 

321 return self.vp/self.vs 

322 

323 def qmu(self): 

324 ''' 

325 Get shear attenuation coefficient Qmu. 

326 ''' 

327 return self.qs 

328 

329 def qk(self): 

330 ''' 

331 Get bulk attenuation coefficient Qk. 

332 ''' 

333 if self.vs == 0. and self.qs == 0.: 

334 return self.qp 

335 else: 

336 s = (4.0/3.0)*(self.vs/self.vp)**2 

337 denom = (1/self.qp - s/self.qs) 

338 if denom <= 0.0: 

339 return num.inf 

340 else: 

341 return (1.-s)/(1.0/self.qp - s/self.qs) 

342 

343 def burgers(self): 

344 ''' 

345 Get Burger parameters. 

346 ''' 

347 return self.burger_eta1, self.burger_eta2, self.burger_valpha 

348 

349 def _rayleigh_equation(self, cr): 

350 cr_a = (cr/self.vp)**2 

351 cr_b = (cr/self.vs)**2 

352 if cr_a > 1.0 or cr_b > 1.0: 

353 return None 

354 

355 return (2.0-cr_b)**2 - 4.0 * math.sqrt(1.0-cr_a) * math.sqrt(1.0-cr_b) 

356 

357 def rayleigh(self): 

358 ''' 

359 Get Rayleigh velocity assuming a homogenous halfspace. 

360 

361 Returned units are [m/s]. 

362 ''' 

363 return bisect(self._rayleigh_equation, 0.001*self.vs, self.vs) 

364 

365 def _has_default_burgers(self): 

366 if self.burger_eta1 == DEFAULT_BURGERS[0] and \ 

367 self.burger_eta2 == DEFAULT_BURGERS[1] and \ 

368 self.burger_valpha == DEFAULT_BURGERS[2]: 

369 return True 

370 return False 

371 

372 def describe(self): 

373 ''' 

374 Get a readable listing of the material properties. 

375 ''' 

376 template = ''' 

377P wave velocity [km/s] : %12g 

378S wave velocity [km/s] : %12g 

379P/S wave vel. ratio : %12g 

380Lame lambda [GPa] : %12g 

381Lame shear modulus [GPa] : %12g 

382Poisson ratio : %12g 

383Bulk modulus [GPa] : %12g 

384Young's modulus [GPa] : %12g 

385Rayleigh wave vel. [km/s] : %12g 

386Density [g/cm**3] : %12g 

387Qp P-wave attenuation : %12g 

388Qs S-wave attenuation (Qmu) : %12g 

389Qk bulk attenuation : %12g 

390transient viscos., eta1 [GPa] : %12g 

391st.-state viscos., eta2 [GPa] : %12g 

392relaxation: valpha : %12g 

393'''.strip() 

394 

395 return template % ( 

396 self.vp/km, 

397 self.vs/km, 

398 self.vp/self.vs, 

399 self.lame_lambda()*1e-9, 

400 self.shear_modulus()*1e-9, 

401 self.poisson(), 

402 self.bulk()*1e-9, 

403 self.youngs()*1e-9, 

404 self.rayleigh()/km, 

405 self.rho/km, 

406 self.qp, 

407 self.qs, 

408 self.qk(), 

409 self.burger_eta1*1e-9, 

410 self.burger_eta2*1e-9, 

411 self.burger_valpha) 

412 

413 def __str__(self): 

414 vp, vs, rho, qp, qs = self.astuple() 

415 return '%10g km/s %10g km/s %10g g/cm^3 %10g %10g' % ( 

416 vp/km, vs/km, rho/km, qp, qs) 

417 

418 def __repr__(self): 

419 return 'Material(vp=%s, vs=%s, rho=%s, qp=%s, qs=%s)' % \ 

420 tuple(repr(x) for x in ( 

421 self.vp, self.vs, self.rho, self.qp, self.qs)) 

422 

423 

424class Leg(object): 

425 ''' 

426 Represents a continuous piece of wave propagation in a phase definition. 

427 

428 **Attributes:** 

429 

430 To be considered as read-only. 

431 

432 .. py:attribute:: departure 

433 

434 One of the constants :py:const:`UP` or :py:const:`DOWN` indicating 

435 upward or downward departure. 

436 

437 .. py:attribute:: mode 

438 

439 One of the constants :py:const:`P` or :py:const:`S`, indicating the 

440 propagation mode. 

441 

442 .. py:attribute:: depthmin 

443 

444 ``None``, a number (a depth in [m]) or a string (an interface name), 

445 minimum depth. 

446 

447 .. py:attribute:: depthmax 

448 

449 ``None``, a number (a depth in [m]) or a string (an interface name), 

450 maximum depth. 

451 

452 ''' 

453 

454 def __init__(self, departure=None, mode=None): 

455 self.departure = departure 

456 self.mode = mode 

457 self.depthmin = None 

458 self.depthmax = None 

459 

460 def set_depthmin(self, depthmin): 

461 self.depthmin = depthmin 

462 

463 def set_depthmax(self, depthmax): 

464 self.depthmax = depthmax 

465 

466 def __str__(self): 

467 def sd(d): 

468 if isinstance(d, float): 

469 return '%g km' % (d/km) 

470 else: 

471 return 'interface %s' % d 

472 

473 s = '%s mode propagation, departing %s' % ( 

474 smode(self.mode).upper(), { 

475 UP: 'upward', DOWN: 'downward'}[self.departure]) 

476 

477 sc = [] 

478 if self.depthmax is not None: 

479 sc.append('deeper than %s' % sd(self.depthmax)) 

480 if self.depthmin is not None: 

481 sc.append('shallower than %s' % sd(self.depthmin)) 

482 

483 if sc: 

484 s = s + ' (may not propagate %s)' % ' or '.join(sc) 

485 

486 return s 

487 

488 

489class InvalidKneeDef(CakeError): 

490 ''' 

491 Invalid definition of a cake :py:class:`Knee`. 

492 ''' 

493 pass 

494 

495 

496class Knee(object): 

497 ''' 

498 Represents a change in wave propagation within a :py:class:`PhaseDef`. 

499 

500 **Attributes:** 

501 

502 To be considered as read-only. 

503 

504 .. py:attribute:: depth 

505 

506 Depth at which the conversion/reflection should happen. this can be 

507 a string or a number. 

508 

509 .. py:attribute:: direction 

510 

511 One of the constants :py:const:`UP` or :py:const:`DOWN` to indicate 

512 the incoming direction. 

513 

514 .. py:attribute:: in_mode 

515 

516 One of the constants :py:const:`P` or :py:const:`S` to indicate the 

517 type of mode of the incoming wave. 

518 

519 .. py:attribute:: out_mode 

520 

521 One of the constants :py:const:`P` or :py:const:`S` to indicate the 

522 type of mode of the outgoing wave. 

523 

524 .. py:attribute:: conversion 

525 

526 Boolean, whether there is a mode conversion involved. 

527 

528 .. py:attribute:: reflection 

529 

530 Boolean, whether there is a reflection involved. 

531 

532 .. py:attribute:: headwave 

533 

534 Boolean, whether there is headwave propagation involved. 

535 

536 ''' 

537 

538 defaults = dict( 

539 depth='surface', 

540 direction=UP, 

541 conversion=True, 

542 reflection=False, 

543 headwave=False, 

544 in_setup_state=True) 

545 

546 defaults_surface = dict( 

547 depth='surface', 

548 direction=UP, 

549 conversion=False, 

550 reflection=True, 

551 headwave=False, 

552 in_setup_state=True) 

553 

554 def __init__(self, *args): 

555 if args: 

556 (self.depth, self.direction, self.reflection, self.in_mode, 

557 self.out_mode) = args 

558 

559 self.conversion = self.in_mode != self.out_mode 

560 self.in_setup_state = False 

561 

562 def default(self, k): 

563 depth = self.__dict__.get('depth', 'surface') 

564 if depth == 'surface': 

565 return Knee.defaults_surface[k] 

566 else: 

567 return Knee.defaults[k] 

568 

569 def __setattr__(self, k, v): 

570 if self.in_setup_state and k in self.__dict__: 

571 raise InvalidKneeDef('%s has already been set' % k) 

572 else: 

573 self.__dict__[k] = v 

574 

575 def __getattr__(self, k): 

576 if k.startswith('__'): 

577 raise AttributeError(k) 

578 

579 if k not in self.__dict__: 

580 return self.default(k) 

581 

582 def set_modes(self, in_leg, out_leg): 

583 

584 if out_leg.departure == UP and ( 

585 (self.direction == UP) == self.reflection): 

586 

587 raise InvalidKneeDef( 

588 'cannot enter %s from %s and emit ray upwards' % ( 

589 ['conversion', 'reflection'][self.reflection], 

590 {UP: 'below', DOWN: 'above'}[self.direction])) 

591 

592 if out_leg.departure == DOWN and ( 

593 (self.direction == DOWN) == self.reflection): 

594 

595 raise InvalidKneeDef( 

596 'cannot enter %s from %s and emit ray downwards' % ( 

597 ['conversion', 'reflection'][self.reflection], 

598 {UP: 'below', DOWN: 'above'}[self.direction])) 

599 

600 self.in_mode = in_leg.mode 

601 self.out_mode = out_leg.mode 

602 

603 def at_surface(self): 

604 return self.depth == 'surface' 

605 

606 def matches(self, discontinuity, mode, direction): 

607 ''' 

608 Check whether it is relevant to a given combination of interface, 

609 propagation mode, and direction. 

610 ''' 

611 

612 if isinstance(self.depth, float): 

613 if abs(self.depth - discontinuity.z) > ZEPS: 

614 return False 

615 else: 

616 if discontinuity.name != self.depth: 

617 return False 

618 

619 return self.direction == direction and self.in_mode == mode 

620 

621 def out_direction(self): 

622 ''' 

623 Get outgoing direction. 

624 

625 Returns one of the constants :py:const:`UP` or :py:const:`DOWN`. 

626 ''' 

627 

628 if self.reflection: 

629 return - self.direction 

630 else: 

631 return self.direction 

632 

633 def __str__(self): 

634 x = [] 

635 if self.reflection: 

636 if self.at_surface(): 

637 x.append('surface') 

638 else: 

639 if not self.headwave: 

640 if self.direction == UP: 

641 x.append('underside') 

642 else: 

643 x.append('upperside') 

644 

645 if self.headwave: 

646 x.append('headwave propagation along') 

647 elif self.reflection and self.conversion: 

648 x.append('reflection with conversion from %s to %s' % ( 

649 smode(self.in_mode).upper(), smode(self.out_mode).upper())) 

650 if not self.at_surface(): 

651 x.append('at') 

652 elif self.reflection: 

653 x.append('reflection') 

654 if not self.at_surface(): 

655 x.append('at') 

656 elif self.conversion: 

657 x.append('conversion from %s to %s at' % ( 

658 smode(self.in_mode).upper(), smode(self.out_mode).upper())) 

659 else: 

660 x.append('passing through') 

661 

662 if isinstance(self.depth, float): 

663 x.append('interface in %g km depth' % (self.depth/1000.)) 

664 else: 

665 if not self.at_surface(): 

666 x.append('%s' % self.depth) 

667 

668 if not self.reflection: 

669 if self.direction == UP: 

670 x.append('on upgoing path') 

671 else: 

672 x.append('on downgoing path') 

673 

674 return ' '.join(x) 

675 

676 

677class Head(Knee): 

678 def __init__(self, *args): 

679 if args: 

680 z, in_direction, mode = args 

681 Knee.__init__(self, z, in_direction, True, mode, mode) 

682 else: 

683 Knee.__init__(self) 

684 

685 def __str__(self): 

686 x = ['propagation as headwave'] 

687 if isinstance(self.depth, float): 

688 x.append('at interface in %g km depth' % (self.depth/1000.)) 

689 else: 

690 x.append('at %s' % self.depth) 

691 

692 return ' '.join(x) 

693 

694 

695class UnknownClassicPhase(CakeError): 

696 ''' 

697 Raised when an invalid classic phase name has been specified.' 

698 ''' 

699 def __init__(self, phasename): 

700 self.phasename = phasename 

701 

702 def __str__(self): 

703 return 'Unknown classic phase name: %s' % self.phasename 

704 

705 

706class PhaseDefParseError(CakeError): 

707 ''' 

708 Exception raised when an error occures during parsing of a phase 

709 definition string. 

710 ''' 

711 

712 def __init__(self, definition, position, exception): 

713 self.definition = definition 

714 self.position = position 

715 self.exception = exception 

716 

717 def __str__(self): 

718 return 'Invalid phase definition: "%s" (at character %i: %s)' % ( 

719 self.definition, self.position+1, str(self.exception)) 

720 

721 

722class PhaseDef(object): 

723 

724 ''' 

725 Definition of a seismic phase arrival, based on ray propagation path. 

726 

727 :param definition: string representation of the phase in Cake's phase 

728 syntax 

729 

730 Seismic phases are conventionally named e.g. P, Pn, PP, PcP, etc. In Cake, 

731 a slightly different terminology is adapted, which allows to specify 

732 arbitrary conversion/reflection histories for seismic ray paths. The 

733 conventions used here are inspired by those used in the TauP toolkit, but 

734 are not completely compatible with those. 

735 

736 The definition of a seismic ray propagation path in Cake's phase syntax is 

737 a string consisting of an alternating sequence of *legs* and *knees*. 

738 

739 A *leg* represents seismic wave propagation without any conversions, 

740 encountering only super-critical reflections. Legs are denoted by ``P``, 

741 ``p``, ``S``, or ``s``. The capital letters are used when the take-off of 

742 the *leg* is in downward direction, while the lower case letters indicate a 

743 take-off in upward direction. 

744 

745 A *knee* is an interaction with an interface. It can be a mode conversion, 

746 a reflection, or propagation as a headwave or diffracted wave. 

747 

748 * conversion is simply denoted as: ``(INTERFACE)`` or ``DEPTH`` 

749 * upperside reflection: ``v(INTERFACE)`` or ``vDEPTH`` 

750 * underside reflection: ``^(INTERFACE)`` or ``^DEPTH`` 

751 * normal kind headwave or diffracted wave: ``v_(INTERFACE)`` or 

752 ``v_DEPTH`` 

753 

754 The interface may be given by name or by depth: INTERFACE is the name of an 

755 interface defined in the model, DEPTH is the depth of an interface in 

756 [km] (the interface closest to that depth is chosen). If two legs appear 

757 consecutively without an explicit *knee*, surface interaction is assumed. 

758 

759 The phase definition may end with a backslash ``\\``, to indicate that the 

760 ray should arrive at the receiver from above instead of from below. It is 

761 possible to restrict the maximum and minimum depth of a *leg* by appending 

762 ``<(INTERFACE)`` or ``<DEPTH`` or ``>(INTERFACE)`` or ``>DEPTH`` after the 

763 leg character, respectively. 

764 

765 **Examples:** 

766 

767 * ``P`` - like the classical P, but includes PKP, PKIKP, Pg 

768 * ``P<(moho)`` - like classical Pg, but must leave source downwards 

769 * ``pP`` - leaves source upward, reflects at surface, then travels as P 

770 * ``P(moho)s`` - conversion from P to S at the Moho on upgoing path 

771 * ``P(moho)S`` - conversion from P to S at the Moho on downgoing path 

772 * ``Pv12p`` - P with reflection at 12 km deep interface (or the 

773 interface closest to that) 

774 * ``Pv_(moho)p`` - classical Pn 

775 * ``Pv_(cmb)p`` - classical Pdiff 

776 * ``P^(conrad)P`` - underside reflection of P at the Conrad 

777 discontinuity 

778 

779 **Usage:** 

780 

781 >>> from pyrocko.cake import PhaseDef 

782 # must escape the backslash 

783 >>> my_crazy_phase = PhaseDef('pPv(moho)sP\\\\') 

784 >>> print my_crazy_phase 

785 Phase definition "pPv(moho)sP\": 

786 - P mode propagation, departing upward 

787 - surface reflection 

788 - P mode propagation, departing downward 

789 - upperside reflection with conversion from P to S at moho 

790 - S mode propagation, departing upward 

791 - surface reflection with conversion from S to P 

792 - P mode propagation, departing downward 

793 - arriving at target from above 

794 

795 .. note:: 

796 

797 (1) These conventions might be extended in a way to allow to fix wave 

798 propagation to SH mode, possibly by specifying SH, or a single 

799 character (e.g. H) instead of S. This would be benificial for the 

800 selection of conversion and reflection coefficients, which 

801 currently only deal with the P-SV case. 

802 ''' 

803 

804 allowed_characters_pattern = r'[0-9a-zA-Z_()<>^v\\.]+' 

805 allowed_characters_pattern_classic = r'[a-zA-Z0-9]+' 

806 

807 @staticmethod 

808 def classic_definitions(): 

809 defs = {} 

810 # PmP, PmS, PcP, PcS, SmP, ... 

811 for r in 'mc': 

812 for a, b in 'PP PS SS SP'.split(): 

813 defs[a+r+b] = [ 

814 '%sv(%s)%s' % (a, {'m': 'moho', 'c': 'cmb'}[r], b.lower())] 

815 

816 # Pg, P, S, Sg 

817 for a in 'PS': 

818 defs[a+'g'] = ['%s<(moho)' % x for x in (a, a.lower())] 

819 defs[a] = ['%s<(cmb)(moho)%s' % (x, x.lower()) for x in ( 

820 a, a.lower())] 

821 

822 defs[a.lower()] = [a.lower()] 

823 

824 for a, b in 'PP PS SS SP'.split(): 

825 defs[a+'K'+b] = ['%s(cmb)P<(icb)(cmb)%s' % (a, b.lower())] 

826 defs[a+'KIK'+b] = ['%s(cmb)P(icb)P(icb)p(cmb)%s' % (a, b.lower())] 

827 defs[a+'KJK'+b] = ['%s(cmb)P(icb)S(icb)p(cmb)%s' % (a, b.lower())] 

828 defs[a+'KiK'+b] = ['%s(cmb)Pv(icb)p(cmb)%s' % (a, b.lower())] 

829 

830 # PP, SS, PS, SP, PPP, ... 

831 for a in 'PS': 

832 for b in 'PS': 

833 for c in 'PS': 

834 defs[a+b+c] = [''.join(defs[x][0] for x in a+b+c)] 

835 

836 defs[a+b] = [''.join(defs[x][0] for x in a+b)] 

837 

838 # Pc, Pdiff, Sc, ... 

839 for x in 'PS': 

840 defs[x+'c'] = defs[x+'diff'] = [x+'v_(cmb)'+x.lower()] 

841 defs[x+'n'] = [x+'v_(moho)'+x.lower()] 

842 

843 # depth phases 

844 for k in list(defs.keys()): 

845 if k not in 'ps': 

846 for x in 'ps': 

847 defs[x+k] = [x + defs[k][0]] 

848 

849 return defs 

850 

851 @staticmethod 

852 def classic(phasename): 

853 ''' 

854 Get phase definitions based on classic phase name. 

855 

856 :param phasename: classic name of a phase 

857 :returns: list of PhaseDef objects 

858 

859 This returns a list of PhaseDef objects, because some classic phases 

860 (like e.g. Pg) can only be represented by two Cake style PhaseDef 

861 objects (one with downgoing and one with upgoing first leg). 

862 ''' 

863 

864 defs = PhaseDef.classic_definitions() 

865 if phasename not in defs: 

866 raise UnknownClassicPhase(phasename) 

867 

868 return [PhaseDef(d, classicname=phasename) for d in defs[phasename]] 

869 

870 def __init__(self, definition=None, classicname=None): 

871 

872 state = 0 

873 sdepth = '' 

874 sinterface = '' 

875 depthmax = depthmin = None 

876 depthlim = None 

877 depthlimtype = None 

878 sdepthlim = '' 

879 events = [] 

880 direction_stop = UP 

881 need_leg = True 

882 ic = 0 

883 if definition is not None: 

884 knee = Knee() 

885 try: 

886 for ic, c in enumerate(definition): 

887 

888 if state in (0, 1): 

889 

890 if c in '0123456789.': 

891 need_leg = True 

892 state = 1 

893 sdepth += c 

894 continue 

895 

896 elif state == 1: 

897 knee.depth = float(sdepth)*1000. 

898 state = 0 

899 

900 if state == 2: 

901 if c == ')': 

902 knee.depth = sinterface 

903 state = 0 

904 else: 

905 sinterface += c 

906 

907 continue 

908 

909 if state in (3, 4): 

910 

911 if state == 3: 

912 if c in '0123456789.': 

913 sdepthlim += c 

914 continue 

915 elif c == '(': 

916 state = 4 

917 continue 

918 else: 

919 depthlim = float(sdepthlim)*1000. 

920 if depthlimtype == '<': 

921 depthmax = depthlim 

922 else: 

923 depthmin = depthlim 

924 state = 0 

925 

926 elif state == 4: 

927 if c == ')': 

928 depthlim = sdepthlim 

929 if depthlimtype == '<': 

930 depthmax = depthlim 

931 else: 

932 depthmin = depthlim 

933 state = 0 

934 continue 

935 else: 

936 sdepthlim += c 

937 continue 

938 

939 if state == 0: 

940 

941 if c == '(': 

942 need_leg = True 

943 state = 2 

944 continue 

945 

946 elif c in '<>': 

947 state = 3 

948 depthlim = None 

949 sdepthlim = '' 

950 depthlimtype = c 

951 continue 

952 

953 elif c in 'psPS': 

954 leg = Leg() 

955 if c in 'ps': 

956 leg.departure = UP 

957 else: 

958 leg.departure = DOWN 

959 leg.mode = imode(c) 

960 

961 if events: 

962 in_leg = events[-1] 

963 if depthmin is not None: 

964 in_leg.set_depthmin(depthmin) 

965 depthmin = None 

966 if depthmax is not None: 

967 in_leg.set_depthmax(depthmax) 

968 depthmax = None 

969 

970 if in_leg.mode != leg.mode: 

971 knee.conversion = True 

972 else: 

973 knee.conversion = False 

974 

975 if not knee.reflection: 

976 if c in 'ps': 

977 knee.direction = UP 

978 else: 

979 knee.direction = DOWN 

980 

981 knee.set_modes(in_leg, leg) 

982 knee.in_setup_state = False 

983 events.append(knee) 

984 knee = Knee() 

985 sdepth = '' 

986 sinterface = '' 

987 

988 events.append(leg) 

989 need_leg = False 

990 continue 

991 

992 elif c == '^': 

993 need_leg = True 

994 knee.direction = UP 

995 knee.reflection = True 

996 continue 

997 

998 elif c == 'v': 

999 need_leg = True 

1000 knee.direction = DOWN 

1001 knee.reflection = True 

1002 continue 

1003 

1004 elif c == '_': 

1005 need_leg = True 

1006 knee.headwave = True 

1007 continue 

1008 

1009 elif c == '\\': 

1010 direction_stop = DOWN 

1011 continue 

1012 

1013 else: 

1014 raise PhaseDefParseError( 

1015 definition, ic, 'invalid character: "%s"' % c) 

1016 

1017 if state == 3: 

1018 depthlim = float(sdepthlim)*1000. 

1019 if depthlimtype == '<': 

1020 depthmax = depthlim 

1021 else: 

1022 depthmin = depthlim 

1023 state = 0 

1024 

1025 except (ValueError, InvalidKneeDef) as e: 

1026 raise PhaseDefParseError(definition, ic, e) 

1027 

1028 if state != 0 or need_leg: 

1029 raise PhaseDefParseError( 

1030 definition, ic, 'unfinished expression') 

1031 

1032 if events and depthmin is not None: 

1033 events[-1].set_depthmin(depthmin) 

1034 if events and depthmax is not None: 

1035 events[-1].set_depthmax(depthmax) 

1036 

1037 self._definition = definition 

1038 self._classicname = classicname 

1039 self._events = events 

1040 self._direction_stop = direction_stop 

1041 

1042 def __iter__(self): 

1043 for ev in self._events: 

1044 yield ev 

1045 

1046 def append(self, ev): 

1047 self._events.append(ev) 

1048 

1049 def first_leg(self): 

1050 ''' 

1051 Get the first leg in phase definition. 

1052 ''' 

1053 return self._events[0] 

1054 

1055 def last_leg(self): 

1056 ''' 

1057 Get the last leg in phase definition. 

1058 ''' 

1059 return self._events[-1] 

1060 

1061 def legs(self): 

1062 ''' 

1063 Iterate over the continuous pieces of wave propagation (legs) defined 

1064 within this phase definition. 

1065 ''' 

1066 

1067 return (leg for leg in self if isinstance(leg, Leg)) 

1068 

1069 def knees(self): 

1070 ''' 

1071 Iterate over conversions and reflections (knees) defined within this 

1072 phase definition. 

1073 ''' 

1074 return (knee for knee in self if isinstance(knee, Knee)) 

1075 

1076 def definition(self): 

1077 ''' 

1078 Get original definition of the phase. 

1079 ''' 

1080 return self._definition 

1081 

1082 def given_name(self): 

1083 ''' 

1084 Get entered classic name if any, or original definition of the phase. 

1085 ''' 

1086 

1087 if self._classicname: 

1088 return self._classicname 

1089 else: 

1090 return self._definition 

1091 

1092 def direction_start(self): 

1093 return self.first_leg().departure 

1094 

1095 def direction_stop(self): 

1096 return self._direction_stop 

1097 

1098 def headwave_knee(self): 

1099 for el in self: 

1100 if type(el) == Knee and el.headwave: 

1101 return el 

1102 return None 

1103 

1104 def used_repr(self): 

1105 ''' 

1106 Translate into textual representation (cake phase syntax). 

1107 ''' 

1108 def strdepth(x): 

1109 if isinstance(x, float): 

1110 return '%g' % (x/1000.) 

1111 else: 

1112 return '(%s)' % x 

1113 

1114 x = [] 

1115 for el in self: 

1116 if type(el) == Leg: 

1117 if el.departure == UP: 

1118 x.append(smode(el.mode).lower()) 

1119 else: 

1120 x.append(smode(el.mode).upper()) 

1121 

1122 if el.depthmax is not None: 

1123 x.append('<'+strdepth(el.depthmax)) 

1124 

1125 if el.depthmin is not None: 

1126 x.append('>'+strdepth(el.depthmin)) 

1127 

1128 elif type(el) == Knee: 

1129 if el.reflection and not el.at_surface(): 

1130 if el.direction == DOWN: 

1131 x.append('v') 

1132 else: 

1133 x.append('^') 

1134 if el.headwave: 

1135 x.append('_') 

1136 if not el.at_surface(): 

1137 x.append(strdepth(el.depth)) 

1138 

1139 elif type(el) == Head: 

1140 x.append('_') 

1141 x.append(strdepth(el.depth)) 

1142 

1143 if self._direction_stop == DOWN: 

1144 x.append('\\') 

1145 

1146 return ''.join(x) 

1147 

1148 def __repr__(self): 

1149 if self._definition is not None: 

1150 return "PhaseDef('%s')" % self._definition 

1151 else: 

1152 return "PhaseDef('%s')" % self.used_repr() 

1153 

1154 def __str__(self): 

1155 orig = '' 

1156 used = self.used_repr() 

1157 if self._definition != used: 

1158 orig = ' (entered as "%s")' % self._definition 

1159 

1160 sarrive = '\n - arriving at target from %s' % ('below', 'above')[ 

1161 self._direction_stop == DOWN] 

1162 

1163 return 'Phase definition "%s"%s:\n - ' % (used, orig) + \ 

1164 '\n - '.join(str(ev) for ev in self) + sarrive 

1165 

1166 def copy(self): 

1167 ''' 

1168 Get a deep copy of it. 

1169 ''' 

1170 return copy.deepcopy(self) 

1171 

1172 

1173def to_phase_defs(phases): 

1174 if isinstance(phases, (str, PhaseDef)): 

1175 phases = [phases] 

1176 

1177 phases_out = [] 

1178 for phase in phases: 

1179 if isinstance(phase, str): 

1180 phases_out.extend(PhaseDef(x.strip()) for x in phase.split(',')) 

1181 elif isinstance(phase, PhaseDef): 

1182 phases_out.append(phase) 

1183 else: 

1184 raise PhaseDefParseError('invalid phase definition') 

1185 

1186 return phases_out 

1187 

1188 

1189def csswap(x): 

1190 return cmath.sqrt(1.-x**2) 

1191 

1192 

1193def psv_surface_ind(in_mode, out_mode): 

1194 ''' 

1195 Get indices to select the appropriate element from scatter matrix for free 

1196 surface. 

1197 ''' 

1198 

1199 return (int(in_mode == S), int(out_mode == S)) 

1200 

1201 

1202def psv_surface(material, p, energy=False): 

1203 ''' 

1204 Scatter matrix for free surface reflection/conversions. 

1205 

1206 :param material: material, object of type :py:class:`Material` 

1207 :param p: flat ray parameter [s/m] 

1208 :param energy: bool, when ``True`` energy normalized coefficients are 

1209 returned 

1210 :returns: Scatter matrix 

1211 

1212 The scatter matrix is ordered as follows:: 

1213 

1214 [[PP, PS], 

1215 [SP, SS]] 

1216 

1217 The formulas given in Aki & Richards are used. 

1218 ''' 

1219 

1220 vp, vs, rho = material.vp, material.vs, material.rho 

1221 sinphi = p * vp 

1222 sinlam = p * vs 

1223 cosphi = csswap(sinphi) 

1224 coslam = csswap(sinlam) 

1225 

1226 if vs == 0.0: 

1227 scatter = num.array([[-1.0, 0.0], [0.0, 1.0]]) 

1228 

1229 else: 

1230 vsp_term = (1.0/vs**2 - 2.0*p**2) 

1231 pcc_term = 4.0 * p**2 * cosphi/vp * coslam/vs 

1232 denom = vsp_term**2 + pcc_term 

1233 

1234 scatter = num.array([ 

1235 [- vsp_term**2 + pcc_term, 4.0*p*coslam/vp*vsp_term], 

1236 [4.0*p*cosphi/vs*vsp_term, vsp_term**2 - pcc_term]], 

1237 dtype=complex) / denom 

1238 

1239 if not energy: 

1240 return scatter 

1241 else: 

1242 eps = 1e-16 

1243 normvec = num.array([vp*rho*cosphi+eps, vs*rho*coslam+eps]) 

1244 escatter = scatter*num.conj(scatter) * num.real( 

1245 (normvec[:, num.newaxis]) / (normvec[num.newaxis, :])) 

1246 return num.real(escatter) 

1247 

1248 

1249def psv_solid_ind(in_direction, out_direction, in_mode, out_mode): 

1250 ''' 

1251 Get indices to select the appropriate element from scatter matrix for 

1252 solid-solid interface. 

1253 ''' 

1254 

1255 return ( 

1256 (out_direction == DOWN)*2 + (out_mode == S), 

1257 (in_direction == UP)*2 + (in_mode == S)) 

1258 

1259 

1260def psv_solid(material1, material2, p, energy=False): 

1261 ''' 

1262 Scatter matrix for solid-solid interface. 

1263 

1264 :param material1: material above, object of type :py:class:`Material` 

1265 :param material2: material below, object of type :py:class:`Material` 

1266 :param p: flat ray parameter [s/m] 

1267 :param energy: bool, when ``True`` energy normalized coefficients are 

1268 returned 

1269 :returns: Scatter matrix 

1270 

1271 The scatter matrix is ordered as follows:: 

1272 

1273 [[P1P1, S1P1, P2P1, S2P1], 

1274 [P1S1, S1S1, P2S1, S2S1], 

1275 [P1P2, S1P2, P2P2, S2P2], 

1276 [P1S2, S1S2, P2S2, S2S2]] 

1277 

1278 The formulas given in Aki & Richards are used. 

1279 ''' 

1280 

1281 vp1, vs1, rho1 = material1.vp, material1.vs, material1.rho 

1282 vp2, vs2, rho2 = material2.vp, material2.vs, material2.rho 

1283 

1284 sinphi1 = p * vp1 

1285 cosphi1 = csswap(sinphi1) 

1286 sinlam1 = p * vs1 

1287 coslam1 = csswap(sinlam1) 

1288 sinphi2 = p * vp2 

1289 cosphi2 = csswap(sinphi2) 

1290 sinlam2 = p * vs2 

1291 coslam2 = csswap(sinlam2) 

1292 

1293 # from aki and richards 

1294 M = num.array([ 

1295 [-vp1*p, -coslam1, vp2*p, coslam2], 

1296 [cosphi1, -vs1*p, cosphi2, -vs2*p], 

1297 [2.0*rho1*vs1**2*p*cosphi1, rho1*vs1*(1.0-2.0*vs1**2*p**2), 

1298 2.0*rho2*vs2**2*p*cosphi2, rho2*vs2*(1.0-2.0*vs2**2*p**2)], 

1299 [-rho1*vp1*(1.0-2.0*vs1**2*p**2), 2.0*rho1*vs1**2*p*coslam1, 

1300 rho2*vp2*(1.0-2.0*vs2**2*p**2), -2.0*rho2*vs2**2*p*coslam2]], 

1301 dtype=complex) 

1302 N = M.copy() 

1303 N[0] *= -1.0 

1304 N[3] *= -1.0 

1305 

1306 scatter = num.dot(num.linalg.inv(M), N) 

1307 

1308 if not energy: 

1309 return scatter 

1310 else: 

1311 eps = 1e-16 

1312 if vs1 == 0.: 

1313 vs1 = vp1*1e-16 

1314 if vs2 == 0.: 

1315 vs2 = vp2*1e-16 

1316 normvec = num.array([ 

1317 vp1*rho1*(cosphi1+eps), vs1*rho1*(coslam1+eps), 

1318 vp2*rho2*(cosphi2+eps), vs2*rho2*(coslam2+eps)], dtype=complex) 

1319 escatter = scatter*num.conj(scatter) * num.real( 

1320 normvec[:, num.newaxis] / normvec[num.newaxis, :]) 

1321 

1322 return num.real(escatter) 

1323 

1324 

1325class BadPotIntCoefs(CakeError): 

1326 ''' 

1327 Raised when the potential interpolation for gradient layers has failed. 

1328 ''' 

1329 pass 

1330 

1331 

1332def potint_coefs(c1, c2, r1, r2): # r2 > r1 

1333 eps = r2*1e-9 

1334 if c1 == 0. and c2 == 0.: 

1335 c1c2 = 1. 

1336 else: 

1337 c1c2 = c1/c2 

1338 b = math.log(c1c2)/math.log((r1+eps)/r2) 

1339 if abs(b) > 10.: 

1340 raise BadPotIntCoefs() 

1341 a = c1/(r1+eps)**b 

1342 return a, b 

1343 

1344 

1345def imode(s): 

1346 if s.lower() == 'p': 

1347 return P 

1348 elif s.lower() == 's': 

1349 return S 

1350 

1351 

1352def smode(i): 

1353 if i == P: 

1354 return 'p' 

1355 elif i == S: 

1356 return 's' 

1357 

1358 

1359class PathFailed(CakeError): 

1360 ''' 

1361 Raised when the ray path computation failed. 

1362 ''' 

1363 pass 

1364 

1365 

1366class SurfaceReached(PathFailed): 

1367 ''' 

1368 Raised when the ray hits the surface before completing the phase 

1369 requirements. 

1370 ''' 

1371 pass 

1372 

1373 

1374class BottomReached(PathFailed): 

1375 ''' 

1376 Raised when the ray exits the bottom of the model before completing the 

1377 phase requirements. 

1378 ''' 

1379 pass 

1380 

1381 

1382class MaxDepthReached(PathFailed): 

1383 ''' 

1384 Raised when the phase's maximum depth has been exceeded. 

1385 ''' 

1386 pass 

1387 

1388 

1389class MinDepthReached(PathFailed): 

1390 ''' 

1391 Raised when the phase's minimum depth has been underrun. 

1392 ''' 

1393 pass 

1394 

1395 

1396class Trapped(PathFailed): 

1397 ''' 

1398 Raised when the ray has been trapped and therefore cannot satisfy the phase 

1399 requirements. 

1400 ''' 

1401 pass 

1402 

1403 

1404class NotPhaseConform(PathFailed): 

1405 ''' 

1406 Raised when the phase requirements cannot be met. 

1407 ''' 

1408 pass 

1409 

1410 

1411class CannotPropagate(PathFailed): 

1412 ''' 

1413 Raised when the phase requirements cannot be met. 

1414 ''' 

1415 def __init__(self, direction, ilayer): 

1416 PathFailed.__init__(self) 

1417 self._direction = direction 

1418 self._ilayer = ilayer 

1419 

1420 def __str__(self): 

1421 return 'Cannot enter layer %i from %s' % ( 

1422 self._ilayer, { 

1423 UP: 'below', 

1424 DOWN: 'above'}[self._direction]) 

1425 

1426 

1427class Layer(object): 

1428 ''' 

1429 Representation of a layer in a layered earth model. 

1430 

1431 :param ztop: depth of top of layer 

1432 :param zbot: depth of bottom of layer 

1433 :param name: name of layer (optional) 

1434 

1435 Subclasses are: :py:class:`HomogeneousLayer` and :py:class:`GradientLayer`. 

1436 ''' 

1437 

1438 def __init__(self, ztop, zbot, name=None): 

1439 self.ztop = ztop 

1440 self.zbot = zbot 

1441 self.zmid = (self.ztop + self.zbot) * 0.5 

1442 self.name = name 

1443 self.ilayer = None 

1444 

1445 def _update_potint_coefs(self): 

1446 potint_p = potint_s = False 

1447 try: 

1448 self._ppic = potint_coefs( 

1449 self.mbot.vp, self.mtop.vp, 

1450 radius(self.zbot), radius(self.ztop)) 

1451 potint_p = True 

1452 except BadPotIntCoefs: 

1453 pass 

1454 

1455 potint_s = False 

1456 try: 

1457 self._spic = potint_coefs( 

1458 self.mbot.vs, self.mtop.vs, 

1459 radius(self.zbot), radius(self.ztop)) 

1460 potint_s = True 

1461 except BadPotIntCoefs: 

1462 pass 

1463 

1464 assert P == 1 and S == 2 

1465 self._use_potential_interpolation = (None, potint_p, potint_s) 

1466 

1467 def potint_coefs(self, mode): 

1468 ''' 

1469 Get coefficients for potential interpolation. 

1470 

1471 :param mode: mode of wave propagation, :py:const:`P` or :py:const:`S` 

1472 :returns: coefficients ``(a, b)`` 

1473 ''' 

1474 

1475 if mode == P: 

1476 return self._ppic 

1477 else: 

1478 return self._spic 

1479 

1480 def contains(self, z): 

1481 ''' 

1482 Tolerantly check if a given depth is within the layer 

1483 (including boundaries). 

1484 ''' 

1485 

1486 return self.ztop <= z <= self.zbot or \ 

1487 self.at_bottom(z) or self.at_top(z) 

1488 

1489 def inner(self, z): 

1490 ''' 

1491 Tolerantly check if a given depth is within the layer 

1492 (not including boundaries). 

1493 ''' 

1494 

1495 return self.ztop <= z <= self.zbot and not \ 

1496 self.at_bottom(z) and not \ 

1497 self.at_top(z) 

1498 

1499 def at_bottom(self, z): 

1500 ''' 

1501 Tolerantly check if given depth is at the bottom of the layer. 

1502 ''' 

1503 

1504 return abs(self.zbot - z) < ZEPS 

1505 

1506 def at_top(self, z): 

1507 ''' 

1508 Tolerantly check if given depth is at the top of the layer. 

1509 ''' 

1510 return abs(self.ztop - z) < ZEPS 

1511 

1512 def pflat_top(self, p): 

1513 ''' 

1514 Convert spherical ray parameter to local flat ray parameter for top of 

1515 layer. 

1516 ''' 

1517 return p / (earthradius-self.ztop) 

1518 

1519 def pflat_bottom(self, p): 

1520 ''' 

1521 Convert spherical ray parameter to local flat ray parameter for bottom 

1522 of layer. 

1523 ''' 

1524 return p / (earthradius-self.zbot) 

1525 

1526 def pflat(self, p, z): 

1527 ''' 

1528 Convert spherical ray parameter to local flat ray parameter for 

1529 given depth. 

1530 ''' 

1531 return p / (earthradius-z) 

1532 

1533 def v_potint(self, mode, z): 

1534 a, b = self.potint_coefs(mode) 

1535 return a*(earthradius-z)**b 

1536 

1537 def u_potint(self, mode, z): 

1538 a, b = self.potint_coefs(mode) 

1539 return 1./(a*(earthradius-z)**b) 

1540 

1541 def xt_potint(self, p, mode, zpart=None): 

1542 ''' 

1543 Get travel time and distance for for traversal with given mode and ray 

1544 parameter. 

1545 

1546 :param p: ray parameter (spherical) 

1547 :param mode: mode of propagation (:py:const:`P` or :py:const:`S`) 

1548 :param zpart: if given, tuple with two depths to restrict computation 

1549 to a part of the layer 

1550 

1551 This implementation uses analytic formulas valid for a spherical earth 

1552 in the case where the velocity c within the layer is given by potential 

1553 interpolation of the form 

1554 

1555 c(z) = a*z^b 

1556 ''' 

1557 utop, ubot = self.u_top_bottom(mode) 

1558 a, b = self.potint_coefs(mode) 

1559 ztop = self.ztop 

1560 zbot = self.zbot 

1561 if zpart is not None: 

1562 utop = self.u(mode, zpart[0]) 

1563 ubot = self.u(mode, zpart[1]) 

1564 ztop, zbot = zpart 

1565 utop = 1./(a*(earthradius-ztop)**b) 

1566 ubot = 1./(a*(earthradius-zbot)**b) 

1567 

1568 r1 = radius(zbot) 

1569 r2 = radius(ztop) 

1570 burger_eta1 = r1 * ubot 

1571 burger_eta2 = r2 * utop 

1572 if b != 1: 

1573 def cpe(eta): 

1574 return num.arccos(num.minimum(p/num.maximum(eta, p/2), 1.0)) 

1575 

1576 def sep(eta): 

1577 return num.sqrt(num.maximum(eta**2 - p**2, 0.0)) 

1578 

1579 x = (cpe(burger_eta2)-cpe(burger_eta1))/(1.0-b) 

1580 t = (sep(burger_eta2)-sep(burger_eta1))/(1.0-b) 

1581 else: 

1582 lr = math.log(r2/r1) 

1583 sap = num.sqrt(1.0/a**2 - p**2) 

1584 x = p/sap * lr 

1585 t = 1./(a**2 * sap) 

1586 

1587 x *= r2d 

1588 

1589 return x, t 

1590 

1591 def test(self, p, mode, z): 

1592 ''' 

1593 Check if wave mode can exist for given ray parameter at given depth 

1594 within the layer. 

1595 ''' 

1596 return (self.u(mode, z)*radius(z) - p) > 0. 

1597 

1598 def tests(self, p, mode): 

1599 utop, ubot = self.u_top_bottom(mode) 

1600 return ( 

1601 (utop * radius(self.ztop) - p) > 0., 

1602 (ubot * radius(self.zbot) - p) > 0.) 

1603 

1604 def zturn_potint(self, p, mode): 

1605 ''' 

1606 Get turning depth for given ray parameter and propagation mode. 

1607 ''' 

1608 

1609 a, b = self.potint_coefs(mode) 

1610 r = num.exp(num.log(a*p)/(1.0-b)) 

1611 return earthradius-r 

1612 

1613 def propagate(self, p, mode, direction): 

1614 ''' 

1615 Propagate ray through layer. 

1616 

1617 :param p: ray parameter 

1618 :param mode: propagation mode 

1619 :param direction: in direction (:py:const:`UP` or :py:const:`DOWN` 

1620 ''' 

1621 if direction == DOWN: 

1622 zin, zout = self.ztop, self.zbot 

1623 else: 

1624 zin, zout = self.zbot, self.ztop 

1625 

1626 if self.v(mode, zin) == 0.0 or not self.test(p, mode, zin): 

1627 raise CannotPropagate(direction, self.ilayer) 

1628 

1629 if not self.test(p, mode, zout): 

1630 return -direction 

1631 else: 

1632 return direction 

1633 

1634 def resize(self, depth_min=None, depth_max=None): 

1635 ''' 

1636 Change layer thinkness and interpolate material if required. 

1637 ''' 

1638 if depth_min: 

1639 mtop = self.material(depth_min) 

1640 

1641 if depth_max: 

1642 mbot = self.material(depth_max) 

1643 

1644 self.mtop = mtop if depth_min else self.mtop 

1645 self.mbot = mbot if depth_max else self.mbot 

1646 self.ztop = depth_min if depth_min else self.ztop 

1647 self.zbot = depth_max if depth_max else self.zbot 

1648 self.zmid = self.ztop + (self.zbot - self.ztop)/2. 

1649 

1650 

1651class DoesNotTurn(CakeError): 

1652 pass 

1653 

1654 

1655def radius(z): 

1656 return earthradius - z 

1657 

1658 

1659class HomogeneousLayer(Layer): 

1660 ''' 

1661 Representation of a homogeneous layer in a layered earth model. 

1662 

1663 Base class: :py:class:`Layer`. 

1664 ''' 

1665 

1666 def __init__(self, ztop, zbot, m, name=None): 

1667 Layer.__init__(self, ztop, zbot, name=name) 

1668 self.m = m 

1669 self.mtop = m 

1670 self.mbot = m 

1671 self._update_potint_coefs() 

1672 

1673 def copy(self, ztop=None, zbot=None): 

1674 if ztop is None: 

1675 ztop = self.ztop 

1676 

1677 if zbot is None: 

1678 zbot = self.zbot 

1679 

1680 return HomogeneousLayer(ztop, zbot, self.m, name=self.name) 

1681 

1682 def material(self, z): 

1683 return self.m 

1684 

1685 def u(self, mode, z=None): 

1686 if self._use_potential_interpolation[mode] and z is not None: 

1687 return self.u_potint(mode, z) 

1688 

1689 if mode == P: 

1690 return 1./self.m.vp 

1691 if mode == S: 

1692 return 1./self.m.vs 

1693 

1694 def u_top_bottom(self, mode): 

1695 u = self.u(mode) 

1696 return u, u 

1697 

1698 def v(self, mode, z=None): 

1699 if self._use_potential_interpolation[mode] and z is not None: 

1700 return self.v_potint(mode, z) 

1701 

1702 if mode == P: 

1703 v = self.m.vp 

1704 if mode == S: 

1705 v = self.m.vs 

1706 

1707 if num.isscalar(z): 

1708 return v 

1709 else: 

1710 return filled(v, len(z)) 

1711 

1712 def v_top_bottom(self, mode): 

1713 v = self.v(mode) 

1714 return v, v 

1715 

1716 def xt(self, p, mode, zpart=None): 

1717 if self._use_potential_interpolation[mode]: 

1718 return self.xt_potint(p, mode, zpart) 

1719 

1720 u = self.u(mode) 

1721 pflat = self.pflat_bottom(p) 

1722 if zpart is None: 

1723 dz = (self.zbot - self.ztop) 

1724 else: 

1725 dz = abs(zpart[1]-zpart[0]) 

1726 

1727 u = self.u(mode) 

1728 eps = u*0.001 

1729 denom = num.sqrt(u**2 - pflat**2) + eps 

1730 

1731 x = r2d*pflat/(earthradius-self.zmid) * dz / denom 

1732 t = u**2 * dz / denom 

1733 return x, t 

1734 

1735 def zturn(self, p, mode): 

1736 if self._use_potential_interpolation[mode]: 

1737 return self.zturn_potint(p, mode) 

1738 

1739 raise DoesNotTurn() 

1740 

1741 def split(self, z): 

1742 upper = HomogeneousLayer(self.ztop, z, self.m, name=self.name) 

1743 lower = HomogeneousLayer(z, self.zbot, self.m, name=self.name) 

1744 upper.ilayer = self.ilayer 

1745 lower.ilayer = self.ilayer 

1746 return upper, lower 

1747 

1748 def __str__(self): 

1749 if self.name: 

1750 name = self.name + ' ' 

1751 else: 

1752 name = '' 

1753 

1754 calcmode = ''.join('HP'[self._use_potential_interpolation[mode]] 

1755 for mode in (P, S)) 

1756 

1757 return ' (%i) homogeneous layer %s(%g km - %g km) [%s]\n %s' % ( 

1758 self.ilayer, name, self.ztop/km, self.zbot/km, calcmode, self.m) 

1759 

1760 

1761class GradientLayer(Layer): 

1762 ''' 

1763 Representation of a gradient layer in a layered earth model. 

1764 

1765 Base class: :py:class:`Layer`. 

1766 ''' 

1767 

1768 def __init__(self, ztop, zbot, mtop, mbot, name=None): 

1769 Layer.__init__(self, ztop, zbot, name=name) 

1770 self.mtop = mtop 

1771 self.mbot = mbot 

1772 self._update_potint_coefs() 

1773 

1774 def copy(self, ztop=None, zbot=None): 

1775 if ztop is None: 

1776 ztop = self.ztop 

1777 

1778 if zbot is None: 

1779 zbot = self.zbot 

1780 

1781 return GradientLayer(ztop, zbot, self.mtop, self.mbot, name=self.name) 

1782 

1783 def interpolate(self, z, ptop, pbot): 

1784 return ptop + (z - self.ztop)*(pbot - ptop)/(self.zbot-self.ztop) 

1785 

1786 def material(self, z): 

1787 dtop = self.mtop.astuple() 

1788 dbot = self.mbot.astuple() 

1789 d = [ 

1790 self.interpolate(z, ptop, pbot) 

1791 for (ptop, pbot) in zip(dtop, dbot)] 

1792 

1793 return Material(*d) 

1794 

1795 def u_top_bottom(self, mode): 

1796 if mode == P: 

1797 return 1./self.mtop.vp, 1./self.mbot.vp 

1798 if mode == S: 

1799 return 1./self.mtop.vs, 1./self.mbot.vs 

1800 

1801 def u(self, mode, z): 

1802 if self._use_potential_interpolation[mode]: 

1803 return self.u_potint(mode, z) 

1804 

1805 if mode == P: 

1806 return 1./self.interpolate(z, self.mtop.vp, self.mbot.vp) 

1807 if mode == S: 

1808 return 1./self.interpolate(z, self.mtop.vs, self.mbot.vs) 

1809 

1810 def v_top_bottom(self, mode): 

1811 if mode == P: 

1812 return self.mtop.vp, self.mbot.vp 

1813 if mode == S: 

1814 return self.mtop.vs, self.mbot.vs 

1815 

1816 def v(self, mode, z): 

1817 if self._use_potential_interpolation[mode]: 

1818 return self.v_potint(mode, z) 

1819 

1820 if mode == P: 

1821 return self.interpolate(z, self.mtop.vp, self.mbot.vp) 

1822 if mode == S: 

1823 return self.interpolate(z, self.mtop.vs, self.mbot.vs) 

1824 

1825 def xt(self, p, mode, zpart=None): 

1826 if self._use_potential_interpolation[mode]: 

1827 return self.xt_potint(p, mode, zpart) 

1828 

1829 utop, ubot = self.u_top_bottom(mode) 

1830 b = (1./ubot - 1./utop)/(self.zbot - self.ztop) 

1831 

1832 pflat = self.pflat_bottom(p) 

1833 if zpart is not None: 

1834 utop = self.u(mode, zpart[0]) 

1835 ubot = self.u(mode, zpart[1]) 

1836 

1837 peps = 1e-16 

1838 pdp = pflat + peps 

1839 

1840 def func(u): 

1841 eta = num.sqrt(num.maximum(u**2 - pflat**2, 0.0)) 

1842 xx = eta/u 

1843 tt = num.where( 

1844 pflat <= u, 

1845 num.log(u+eta) - num.log(pdp) - eta/u, 

1846 0.0) 

1847 

1848 return xx, tt 

1849 

1850 xxtop, tttop = func(utop) 

1851 xxbot, ttbot = func(ubot) 

1852 

1853 x = (xxtop - xxbot) / (b*pdp) 

1854 t = (tttop - ttbot) / b + pflat*x 

1855 

1856 x *= r2d/(earthradius - self.zmid) 

1857 return x, t 

1858 

1859 def zturn(self, p, mode): 

1860 if self._use_potential_interpolation[mode]: 

1861 return self.zturn_potint(p, mode) 

1862 pflat = self.pflat_bottom(p) 

1863 vtop, vbot = self.v_top_bottom(mode) 

1864 return (1./pflat - vtop) * (self.zbot - self.ztop) / \ 

1865 (vbot-vtop) + self.ztop 

1866 

1867 def split(self, z): 

1868 mmid = self.material(z) 

1869 upper = GradientLayer(self.ztop, z, self.mtop, mmid, name=self.name) 

1870 lower = GradientLayer(z, self.zbot, mmid, self.mbot, name=self.name) 

1871 upper.ilayer = self.ilayer 

1872 lower.ilayer = self.ilayer 

1873 return upper, lower 

1874 

1875 def __str__(self): 

1876 if self.name: 

1877 name = self.name + ' ' 

1878 else: 

1879 name = '' 

1880 

1881 calcmode = ''.join('HP'[self._use_potential_interpolation[mode]] 

1882 for mode in (P, S)) 

1883 

1884 return ''' (%i) gradient layer %s(%g km - %g km) [%s] 

1885 %s 

1886 %s''' % ( 

1887 self.ilayer, 

1888 name, 

1889 self.ztop/km, 

1890 self.zbot/km, 

1891 calcmode, 

1892 self.mtop, 

1893 self.mbot) 

1894 

1895 

1896class Discontinuity(object): 

1897 ''' 

1898 Base class for discontinuities in layered earth model. 

1899 

1900 Subclasses are: :py:class:`Interface` and :py:class:`Surface`. 

1901 ''' 

1902 

1903 def __init__(self, z, name=None): 

1904 self.z = z 

1905 self.zbot = z 

1906 self.ztop = z 

1907 self.name = name 

1908 

1909 def change_depth(self, z): 

1910 self.z = z 

1911 self.zbot = z 

1912 self.ztop = z 

1913 

1914 def copy(self): 

1915 return copy.deepcopy(self) 

1916 

1917 

1918class Interface(Discontinuity): 

1919 ''' 

1920 Representation of an interface in a layered earth model. 

1921 

1922 Base class: :py:class:`Discontinuity`. 

1923 ''' 

1924 

1925 def __init__(self, z, mabove, mbelow, name=None): 

1926 Discontinuity.__init__(self, z, name) 

1927 self.mabove = mabove 

1928 self.mbelow = mbelow 

1929 

1930 def __str__(self): 

1931 if self.name is None: 

1932 return 'interface' 

1933 else: 

1934 return 'interface "%s"' % self.name 

1935 

1936 def u_top_bottom(self, mode): 

1937 if mode == P: 

1938 return reci_or_none(self.mabove.vp), reci_or_none(self.mbelow.vp) 

1939 if mode == S: 

1940 return reci_or_none(self.mabove.vs), reci_or_none(self.mbelow.vs) 

1941 

1942 def critical_ps(self, mode): 

1943 uabove, ubelow = self.u_top_bottom(mode) 

1944 return ( 

1945 mult_or_none(uabove, radius(self.z)), 

1946 mult_or_none(ubelow, radius(self.z))) 

1947 

1948 def propagate(self, p, mode, direction): 

1949 uabove, ubelow = self.u_top_bottom(mode) 

1950 if direction == DOWN: 

1951 if ubelow is not None and ubelow*radius(self.z) - p >= 0: 

1952 return direction 

1953 else: 

1954 return -direction 

1955 if direction == UP: 

1956 if uabove is not None and uabove*radius(self.z) - p >= 0: 

1957 return direction 

1958 else: 

1959 return -direction 

1960 

1961 def pflat(self, p): 

1962 return p / (earthradius-self.z) 

1963 

1964 def efficiency(self, in_direction, out_direction, in_mode, out_mode, p): 

1965 scatter = psv_solid( 

1966 self.mabove, self.mbelow, self.pflat(p), energy=True) 

1967 return scatter[ 

1968 psv_solid_ind(in_direction, out_direction, in_mode, out_mode)] 

1969 

1970 

1971class Surface(Discontinuity): 

1972 ''' 

1973 Representation of the surface discontinuity in a layered earth model. 

1974 

1975 Base class: :py:class:`Discontinuity`. 

1976 ''' 

1977 

1978 def __init__(self, z, mbelow): 

1979 Discontinuity.__init__(self, z, 'surface') 

1980 self.z = z 

1981 self.mbelow = mbelow 

1982 

1983 def propagate(self, p, mode, direction): 

1984 return direction # no implicit reflection at surface 

1985 

1986 def u_top_bottom(self, mode): 

1987 if mode == P: 

1988 return None, reci_or_none(self.mbelow.vp) 

1989 if mode == S: 

1990 return None, reci_or_none(self.mbelow.vs) 

1991 

1992 def critical_ps(self, mode): 

1993 _, ubelow = self.u_top_bottom(mode) 

1994 return None, mult_or_none(ubelow, radius(self.z)) 

1995 

1996 def pflat(self, p): 

1997 return p / (earthradius-self.z) 

1998 

1999 def efficiency(self, in_direction, out_direction, in_mode, out_mode, p): 

2000 if in_direction == DOWN or out_direction == UP: 

2001 return 0.0 

2002 else: 

2003 return psv_surface( 

2004 self.mbelow, self.pflat(p), energy=True)[ 

2005 psv_surface_ind(in_mode, out_mode)] 

2006 

2007 def __str__(self): 

2008 return 'surface' 

2009 

2010 

2011class Walker(object): 

2012 def __init__(self, elements): 

2013 self._elements = elements 

2014 self._i = 0 

2015 

2016 def current(self): 

2017 return self._elements[self._i] 

2018 

2019 def go(self, direction): 

2020 if direction == UP: 

2021 self.up() 

2022 else: 

2023 self.down() 

2024 

2025 def down(self): 

2026 if self._i < len(self._elements)-1: 

2027 self._i += 1 

2028 else: 

2029 raise BottomReached() 

2030 

2031 def up(self): 

2032 if self._i > 0: 

2033 self._i -= 1 

2034 else: 

2035 raise SurfaceReached() 

2036 

2037 def goto_layer(self, layer): 

2038 self._i = self._elements.index(layer) 

2039 

2040 

2041class RayElement(object): 

2042 ''' 

2043 An element of a :py:class:`RayPath`. 

2044 ''' 

2045 

2046 def __eq__(self, other): 

2047 return type(self) == type(other) and self.__dict__ == other.__dict__ 

2048 

2049 def is_straight(self): 

2050 return isinstance(self, Straight) 

2051 

2052 def is_kink(self): 

2053 return isinstance(self, Kink) 

2054 

2055 

2056class Straight(RayElement): 

2057 ''' 

2058 A ray segment representing wave propagation through one :py:class:`Layer`. 

2059 ''' 

2060 

2061 def __init__(self, direction_in, direction_out, mode, layer): 

2062 self.mode = mode 

2063 self._direction_in = direction_in 

2064 self._direction_out = direction_out 

2065 self.layer = layer 

2066 

2067 def angle_in(self, p, endgaps=None): 

2068 z = self.z_in(endgaps) 

2069 dir = self.eff_direction_in(endgaps) 

2070 v = self.layer.v(self.mode, z) 

2071 pf = self.layer.pflat(p, z) 

2072 

2073 if dir == DOWN: 

2074 return num.arcsin(v*pf)*r2d 

2075 else: 

2076 return 180.-num.arcsin(v*pf)*r2d 

2077 

2078 def angle_out(self, p, endgaps=None): 

2079 z = self.z_out(endgaps) 

2080 dir = self.eff_direction_out(endgaps) 

2081 v = self.layer.v(self.mode, z) 

2082 pf = self.layer.pflat(p, z) 

2083 

2084 if dir == DOWN: 

2085 return 180.-num.arcsin(v*pf)*r2d 

2086 else: 

2087 return num.arcsin(v*pf)*r2d 

2088 

2089 def pflat_in(self, p, endgaps=None): 

2090 return p / (earthradius-self.z_in(endgaps)) 

2091 

2092 def pflat_out(self, p, endgaps=None): 

2093 return p / (earthradius-self.z_out(endgaps)) 

2094 

2095 def test(self, p, z): 

2096 return self.layer.test(p, self.mode, z) 

2097 

2098 def z_in(self, endgaps=None): 

2099 if endgaps is not None: 

2100 return endgaps[0] 

2101 else: 

2102 lyr = self.layer 

2103 return (lyr.ztop, lyr.zbot)[self._direction_in == UP] 

2104 

2105 def z_out(self, endgaps=None): 

2106 if endgaps is not None: 

2107 return endgaps[1] 

2108 else: 

2109 lyr = self.layer 

2110 return (lyr.ztop, lyr.zbot)[self._direction_out == DOWN] 

2111 

2112 def turns(self): 

2113 return self._direction_in != self._direction_out 

2114 

2115 def eff_direction_in(self, endgaps=None): 

2116 if endgaps is None: 

2117 return self._direction_in 

2118 else: 

2119 return endgaps[2] 

2120 

2121 def eff_direction_out(self, endgaps=None): 

2122 if endgaps is None: 

2123 return self._direction_out 

2124 else: 

2125 return endgaps[3] 

2126 

2127 def zturn(self, p): 

2128 lyr = self.layer 

2129 return lyr.zturn(p, self.mode) 

2130 

2131 def u_in(self, endgaps=None): 

2132 return self.layer.u(self.mode, self.z_in(endgaps)) 

2133 

2134 def u_out(self, endgaps=None): 

2135 return self.layer.u(self.mode, self.z_out(endgaps)) 

2136 

2137 def critical_p_in(self, endgaps=None): 

2138 z = self.z_in(endgaps) 

2139 return self.layer.u(self.mode, z)*radius(z) 

2140 

2141 def critical_p_out(self, endgaps=None): 

2142 z = self.z_out(endgaps) 

2143 return self.layer.u(self.mode, z)*radius(z) 

2144 

2145 def xt(self, p, zpart=None): 

2146 x, t = self.layer.xt(p, self.mode, zpart=zpart) 

2147 if self._direction_in != self._direction_out and zpart is None: 

2148 x *= 2. 

2149 t *= 2. 

2150 return x, t 

2151 

2152 def xt_gap(self, p, zstart, zstop, samedir): 

2153 z1, z2 = zstart, zstop 

2154 if z1 > z2: 

2155 z1, z2 = z2, z1 

2156 

2157 x, t = self.layer.xt(p, self.mode, zpart=(z1, z2)) 

2158 

2159 if samedir: 

2160 return x, t 

2161 else: 

2162 xfull, tfull = self.xt(p) 

2163 return xfull-x, tfull-t 

2164 

2165 def __hash__(self): 

2166 return hash(( 

2167 self._direction_in, 

2168 self._direction_out, 

2169 self.mode, 

2170 id(self.layer))) 

2171 

2172 

2173class HeadwaveStraight(Straight): 

2174 def __init__(self, direction_in, direction_out, mode, interface): 

2175 Straight.__init__(self, direction_in, direction_out, mode, None) 

2176 

2177 self.interface = interface 

2178 

2179 def z_in(self, zpart=None): 

2180 return self.interface.z 

2181 

2182 def z_out(self, zpart=None): 

2183 return self.interface.z 

2184 

2185 def zturn(self, p): 

2186 return filled(self.interface.z, len(p)) 

2187 

2188 def xt(self, p, zpart=None): 

2189 return 0., 0. 

2190 

2191 def x2t_headwave(self, xstretch): 

2192 xstretch_m = xstretch*d2r*radius(self.interface.z) 

2193 return min_not_none(*self.interface.u_top_bottom(self.mode))*xstretch_m 

2194 

2195 

2196class Kink(RayElement): 

2197 ''' 

2198 An interaction of a ray with a :py:class:`Discontinuity`. 

2199 ''' 

2200 

2201 def __init__( 

2202 self, 

2203 in_direction, 

2204 out_direction, 

2205 in_mode, 

2206 out_mode, 

2207 discontinuity): 

2208 

2209 self.in_direction = in_direction 

2210 self.out_direction = out_direction 

2211 self.in_mode = in_mode 

2212 self.out_mode = out_mode 

2213 self.discontinuity = discontinuity 

2214 

2215 def reflection(self): 

2216 return self.in_direction != self.out_direction 

2217 

2218 def conversion(self): 

2219 return self.in_mode != self.out_mode 

2220 

2221 def efficiency(self, p, out_direction=None, out_mode=None): 

2222 

2223 if out_direction is None: 

2224 out_direction = self.out_direction 

2225 

2226 if out_mode is None: 

2227 out_mode = self.out_mode 

2228 

2229 return self.discontinuity.efficiency( 

2230 self.in_direction, out_direction, self.in_mode, out_mode, p) 

2231 

2232 def __str__(self): 

2233 r, c = self.reflection(), self.conversion() 

2234 if r and c: 

2235 return '|~' 

2236 if r: 

2237 return '|' 

2238 if c: 

2239 return '~' 

2240 return '_' 

2241 

2242 def __hash__(self): 

2243 return hash(( 

2244 self.in_direction, 

2245 self.out_direction, 

2246 self.in_mode, 

2247 self.out_mode, 

2248 id(self.discontinuity))) 

2249 

2250 

2251class PRangeNotSet(CakeError): 

2252 pass 

2253 

2254 

2255class RayPath(object): 

2256 ''' 

2257 Representation of a fan of rays running through a common sequence of 

2258 layers / interfaces. 

2259 ''' 

2260 

2261 def __init__(self, phase): 

2262 self.elements = [] 

2263 self.phase = phase 

2264 self._pmax = None 

2265 self._pmin = None 

2266 self._p = None 

2267 self._is_headwave = False 

2268 

2269 def set_is_headwave(self, is_headwave): 

2270 self._is_headwave = is_headwave 

2271 

2272 def copy(self): 

2273 ''' 

2274 Get a copy of it. 

2275 ''' 

2276 

2277 c = copy.copy(self) 

2278 c.elements = list(self.elements) 

2279 return c 

2280 

2281 def endgaps(self, zstart, zstop): 

2282 ''' 

2283 Get information needed for end point adjustments. 

2284 ''' 

2285 

2286 return ( 

2287 zstart, 

2288 zstop, 

2289 self.phase.direction_start(), 

2290 self.phase.direction_stop()) 

2291 

2292 def append(self, element): 

2293 self.elements.append(element) 

2294 

2295 def _check_have_prange(self): 

2296 if self._pmax is None: 

2297 raise PRangeNotSet() 

2298 

2299 def set_prange(self, pmin, pmax, dp): 

2300 self._pmin, self._pmax = pmin, pmax 

2301 self._prange_dp = dp 

2302 

2303 def used_phase(self, p=None, eps=1.): 

2304 ''' 

2305 Calculate phase definition from ray path. 

2306 ''' 

2307 

2308 used = PhaseDef() 

2309 fleg = self.phase.first_leg() 

2310 used.append(Leg(fleg.departure, fleg.mode)) 

2311 n_elements_n = [None] + self.elements + [None] 

2312 for before, element, after in zip( 

2313 n_elements_n[:-2], 

2314 n_elements_n[1:-1], 

2315 n_elements_n[2:]): 

2316 

2317 if element.is_kink() and HeadwaveStraight not in ( 

2318 type(before), 

2319 type(after)): 

2320 

2321 if element.reflection() or element.conversion(): 

2322 z = element.discontinuity.z 

2323 used.append(Knee( 

2324 z, 

2325 element.in_direction, 

2326 element.out_direction != element.in_direction, 

2327 element.in_mode, 

2328 element.out_mode)) 

2329 

2330 used.append(Leg(element.out_direction, element.out_mode)) 

2331 

2332 elif type(element) is HeadwaveStraight: 

2333 z = element.interface.z 

2334 k = Knee( 

2335 z, 

2336 before.in_direction, 

2337 after.out_direction != before.in_direction, 

2338 before.in_mode, 

2339 after.out_mode) 

2340 

2341 k.headwave = True 

2342 used.append(k) 

2343 used.append(Leg(after.out_direction, after.out_mode)) 

2344 

2345 if (p is not None and before and after 

2346 and element.is_straight() 

2347 and before.is_kink() 

2348 and after.is_kink() 

2349 and element.turns() 

2350 and not before.reflection() and not before.conversion() 

2351 and not after.reflection() and not after.conversion()): 

2352 

2353 ai = element.angle_in(p) 

2354 if 90.0-eps < ai and ai < 90+eps: 

2355 used.append( 

2356 Head( 

2357 before.discontinuity.z, 

2358 before.out_direction, 

2359 element.mode)) 

2360 used.append( 

2361 Leg(-before.out_direction, element.mode)) 

2362 

2363 used._direction_stop = self.phase.direction_stop() 

2364 used._definition = self.phase.definition() 

2365 

2366 return used 

2367 

2368 def pmax(self): 

2369 ''' 

2370 Get maximum valid ray parameter. 

2371 ''' 

2372 self._check_have_prange() 

2373 return self._pmax 

2374 

2375 def pmin(self): 

2376 ''' 

2377 Get minimum valid ray parameter. 

2378 ''' 

2379 self._check_have_prange() 

2380 return self._pmin 

2381 

2382 def xmin(self): 

2383 ''' 

2384 Get minimal distance. 

2385 ''' 

2386 self._analyse() 

2387 return self._xmin 

2388 

2389 def xmax(self): 

2390 ''' 

2391 Get maximal distance. 

2392 ''' 

2393 self._analyse() 

2394 return self._xmax 

2395 

2396 def kinks(self): 

2397 ''' 

2398 Iterate over propagation mode changes (reflections/transmissions). 

2399 ''' 

2400 return (k for k in self.elements if isinstance(k, Kink)) 

2401 

2402 def straights(self): 

2403 ''' 

2404 Iterate over ray segments. 

2405 ''' 

2406 return (s for s in self.elements if isinstance(s, Straight)) 

2407 

2408 def headwave_straight(self): 

2409 for s in self.elements: 

2410 if type(s) is HeadwaveStraight: 

2411 return s 

2412 

2413 def first_straight(self): 

2414 ''' 

2415 Get first ray segment. 

2416 ''' 

2417 for s in self.elements: 

2418 if isinstance(s, Straight): 

2419 return s 

2420 

2421 def last_straight(self): 

2422 ''' 

2423 Get last ray segment. 

2424 ''' 

2425 for s in reversed(self.elements): 

2426 if isinstance(s, Straight): 

2427 return s 

2428 

2429 def efficiency(self, p): 

2430 ''' 

2431 Get product of all conversion/reflection coefficients encountered on 

2432 path. 

2433 ''' 

2434 return reduce( 

2435 operator.mul, (k.efficiency(p) for k in self.kinks()), 1.) 

2436 

2437 def spreading(self, p, endgaps): 

2438 ''' 

2439 Get geometrical spreading factor. 

2440 ''' 

2441 if self._is_headwave: 

2442 return 0.0 

2443 

2444 self._check_have_prange() 

2445 dp = self._prange_dp * 0.01 

2446 assert self._pmax - self._pmin > dp 

2447 

2448 if p + dp > self._pmax: 

2449 p = p-dp 

2450 

2451 x0, t = self.xt(p, endgaps) 

2452 x1, t = self.xt(p+dp, endgaps) 

2453 x0 *= d2r 

2454 x1 *= d2r 

2455 if x1 == x0: 

2456 return num.nan 

2457 

2458 dp_dx = dp/(x1-x0) 

2459 

2460 x = x0 

2461 if x == 0.: 

2462 x = x1 

2463 p = dp 

2464 

2465 first = self.first_straight() 

2466 last = self.last_straight() 

2467 return num.abs(dp_dx) * first.pflat_in(p, endgaps) / ( 

2468 4.0 * math.pi * num.sin(x) * 

2469 (earthradius-first.z_in(endgaps)) * 

2470 (earthradius-last.z_out(endgaps))**2 * 

2471 first.u_in(endgaps)**2 * 

2472 num.abs(num.cos(first.angle_in(p, endgaps)*d2r)) * 

2473 num.abs(num.cos(last.angle_out(p, endgaps)*d2r))) 

2474 

2475 def make_p(self, dp=None, n=None, nmin=None): 

2476 assert dp is None or n is None 

2477 

2478 if self._pmin == self._pmax: 

2479 return num.array([self._pmin]) 

2480 

2481 if dp is None: 

2482 dp = self._prange_dp 

2483 

2484 if n is None: 

2485 n = int(round((self._pmax-self._pmin)/dp)) + 1 

2486 

2487 if nmin is not None: 

2488 n = max(n, nmin) 

2489 

2490 ppp = num.linspace(self._pmin, self._pmax, n) 

2491 return ppp 

2492 

2493 def xt_endgaps(self, p, endgaps, which='both'): 

2494 ''' 

2495 Get amount of distance/traveltime to be subtracted at the generic ray 

2496 path's ends. 

2497 ''' 

2498 

2499 zstart, zstop, dirstart, dirstop = endgaps 

2500 firsts = self.first_straight() 

2501 lasts = self.last_straight() 

2502 xs, ts = firsts.xt_gap( 

2503 p, zstart, firsts.z_in(), dirstart == firsts._direction_in) 

2504 xe, te = lasts.xt_gap( 

2505 p, zstop, lasts.z_out(), dirstop == lasts._direction_out) 

2506 

2507 if which == 'both': 

2508 return xs + xe, ts + te 

2509 elif which == 'left': 

2510 return xs, ts 

2511 elif which == 'right': 

2512 return xe, te 

2513 

2514 def xt_endgaps_ptest(self, p, endgaps): 

2515 ''' 

2516 Check if ray parameter is valid at source and receiver. 

2517 ''' 

2518 

2519 zstart, zstop, dirstart, dirstop = endgaps 

2520 firsts = self.first_straight() 

2521 lasts = self.last_straight() 

2522 return num.logical_and(firsts.test(p, zstart), lasts.test(p, zstop)) 

2523 

2524 def xt(self, p, endgaps): 

2525 ''' 

2526 Calculate distance and traveltime for given ray parameter. 

2527 ''' 

2528 

2529 if isinstance(p, num.ndarray): 

2530 sx = num.zeros(p.size) 

2531 st = num.zeros(p.size) 

2532 else: 

2533 sx = 0.0 

2534 st = 0.0 

2535 

2536 for s in self.straights(): 

2537 x, t = s.xt(p) 

2538 sx += x 

2539 st += t 

2540 

2541 if endgaps: 

2542 dx, dt = self.xt_endgaps(p, endgaps) 

2543 sx -= dx 

2544 st -= dt 

2545 

2546 return sx, st 

2547 

2548 def xt_limits(self, p): 

2549 ''' 

2550 Calculate limits of distance and traveltime for given ray parameter. 

2551 ''' 

2552 

2553 if isinstance(p, num.ndarray): 

2554 sx = num.zeros(p.size) 

2555 st = num.zeros(p.size) 

2556 sxe = num.zeros(p.size) 

2557 ste = num.zeros(p.size) 

2558 else: 

2559 sx = 0.0 

2560 st = 0.0 

2561 sxe = 0.0 

2562 ste = 0.0 

2563 

2564 sfirst = self.first_straight() 

2565 slast = self.last_straight() 

2566 

2567 for s in self.straights(): 

2568 if s is not sfirst and s is not slast: 

2569 x, t = s.xt(p) 

2570 sx += x 

2571 st += t 

2572 

2573 sends = [sfirst] 

2574 if sfirst is not slast: 

2575 sends.append(slast) 

2576 

2577 for s in sends: 

2578 x, t = s.xt(p) 

2579 sxe += x 

2580 ste += t 

2581 

2582 return sx, (sx + sxe), st, (st + ste) 

2583 

2584 def iter_zxt(self, p): 

2585 ''' 

2586 Iterate over (depth, distance, traveltime) at each layer interface on 

2587 ray path. 

2588 ''' 

2589 

2590 sx = num.zeros(p.size) 

2591 st = num.zeros(p.size) 

2592 ok = False 

2593 for s in self.straights(): 

2594 yield s.z_in(), sx.copy(), st.copy() 

2595 

2596 x, t = s.xt(p) 

2597 sx += x 

2598 st += t 

2599 ok = True 

2600 

2601 if ok: 

2602 yield s.z_out(), sx.copy(), st.copy() 

2603 

2604 def zxt_path_subdivided( 

2605 self, p, endgaps, 

2606 points_per_straight=20, 

2607 x_for_headwave=None): 

2608 

2609 ''' 

2610 Get geometrical representation of ray path. 

2611 ''' 

2612 

2613 if self._is_headwave: 

2614 assert p.size == 1 

2615 x, t = self.xt(p, endgaps) 

2616 xstretch = x_for_headwave-x 

2617 nout = xstretch.size 

2618 else: 

2619 nout = p.size 

2620 

2621 dxl, dtl = self.xt_endgaps(p, endgaps, which='left') 

2622 dxr, dtr = self.xt_endgaps(p, endgaps, which='right') 

2623 

2624 # first create full path including the endgaps 

2625 sx = num.zeros(nout) - dxl 

2626 st = num.zeros(nout) - dtl 

2627 zxt = [] 

2628 for s in self.straights(): 

2629 n = points_per_straight 

2630 

2631 back = None 

2632 zin, zout = s.z_in(), s.z_out() 

2633 if type(s) is HeadwaveStraight: 

2634 z = zin 

2635 for i in range(n): 

2636 xs = float(i)/(n-1) * xstretch 

2637 ts = s.x2t_headwave(xs) 

2638 zxt.append((filled(z, xstretch.size), sx+xs, st+ts)) 

2639 else: 

2640 if zin != zout: # normal traversal 

2641 zs = num.linspace(zin, zout, n).tolist() 

2642 for z in zs: 

2643 x, t = s.xt(p, zpart=sorted([zin, z])) 

2644 zxt.append((filled(z, nout), sx + x, st + t)) 

2645 

2646 else: # ray turns in layer 

2647 zturn = s.zturn(p) 

2648 back = [] 

2649 for i in range(n): 

2650 z = zin + (zturn - zin) * num.sin( 

2651 float(i)/(n-1)*math.pi/2.0) * 0.999 

2652 

2653 if zturn[0] >= zin: 

2654 x, t = s.xt(p, zpart=[zin, z]) 

2655 else: 

2656 x, t = s.xt(p, zpart=[z, zin]) 

2657 zxt.append((z, sx + x, st + t)) 

2658 back.append((z, x, t)) 

2659 

2660 if type(s) is HeadwaveStraight: 

2661 x = xstretch 

2662 t = s.x2t_headwave(xstretch) 

2663 else: 

2664 x, t = s.xt(p) 

2665 

2666 sx += x 

2667 st += t 

2668 if back: 

2669 for z, x, t in reversed(back): 

2670 zxt.append((z, sx - x, st - t)) 

2671 

2672 # gather results as arrays with such that x[ip, ipoint] 

2673 fanz, fanx, fant = [], [], [] 

2674 for z, x, t in zxt: 

2675 fanz.append(z) 

2676 fanx.append(x) 

2677 fant.append(t) 

2678 

2679 z = num.array(fanz).T 

2680 x = num.array(fanx).T 

2681 t = num.array(fant).T 

2682 

2683 # cut off the endgaps, add exact endpoints 

2684 xmax = x[:, -1] - dxr 

2685 tmax = t[:, -1] - dtr 

2686 zstart, zstop = endgaps[:2] 

2687 zs, xs, ts = [], [], [] 

2688 for i in range(nout): 

2689 t_ = t[i] 

2690 indices = num.where(num.logical_and(0. <= t_, t_ <= tmax[i]))[0] 

2691 n = indices.size + 2 

2692 zs_, xs_, ts_ = [num.empty(n, dtype=float) for j in range(3)] 

2693 zs_[1:-1] = z[i, indices] 

2694 xs_[1:-1] = x[i, indices] 

2695 ts_[1:-1] = t[i, indices] 

2696 zs_[0], zs_[-1] = zstart, zstop 

2697 xs_[0], xs_[-1] = 0., xmax[i] 

2698 ts_[0], ts_[-1] = 0., tmax[i] 

2699 zs.append(zs_) 

2700 xs.append(xs_) 

2701 ts.append(ts_) 

2702 

2703 return zs, xs, ts 

2704 

2705 def _analyse(self): 

2706 if self._p is not None: 

2707 return 

2708 

2709 p = self.make_p(nmin=20) 

2710 xmin, xmax, tmin, tmax = self.xt_limits(p) 

2711 

2712 self._x, self._t, self._p = xmax, tmax, p 

2713 self._xmin, self._xmax = xmin.min(), xmax.max() 

2714 self._tmin, self._tmax = tmin.min(), tmax.max() 

2715 

2716 def draft_pxt(self, endgaps): 

2717 self._analyse() 

2718 

2719 if not self._is_headwave: 

2720 cp, cx, ct = self._p, self._x, self._t 

2721 pcrit = min( 

2722 self.critical_pstart(endgaps), 

2723 self.critical_pstop(endgaps)) 

2724 

2725 if pcrit < self._pmin: 

2726 empty = num.array([], dtype=float) 

2727 return empty, empty, empty 

2728 

2729 elif pcrit >= self._pmax: 

2730 dx, dt = self.xt_endgaps(cp, endgaps) 

2731 return cp, cx-dx, ct-dt 

2732 

2733 else: 

2734 n = num.searchsorted(cp, pcrit) + 1 

2735 rp, rx, rt = num.empty((3, n), dtype=float) 

2736 rp[:-1] = cp[:n-1] 

2737 rx[:-1] = cx[:n-1] 

2738 rt[:-1] = ct[:n-1] 

2739 rp[-1] = pcrit 

2740 rx[-1], rt[-1] = self.xt(pcrit, endgaps) 

2741 dx, dt = self.xt_endgaps(rp, endgaps) 

2742 rx[:-1] -= dx[:-1] 

2743 rt[:-1] -= dt[:-1] 

2744 return rp, rx, rt 

2745 

2746 else: 

2747 dx, dt = self.xt_endgaps(self._p, endgaps) 

2748 p, x, t = self._p, self._x - dx, self._t - dt 

2749 p, x, t = p[0], x[0], t[0] 

2750 xh = num.linspace(0., x*10-x, 10) 

2751 th = self.headwave_straight().x2t_headwave(xh) 

2752 return filled(p, xh.size), x+xh, t+th 

2753 

2754 def interpolate_x2pt_linear(self, x, endgaps): 

2755 ''' 

2756 Get approximate ray parameter and traveltime for distance. 

2757 ''' 

2758 

2759 self._analyse() 

2760 

2761 if self._is_headwave: 

2762 dx, dt = self.xt_endgaps(self._p, endgaps) 

2763 xmin = self._x[0] - dx[0] 

2764 tmin = self._t[0] - dt[0] 

2765 el = self.headwave_straight() 

2766 xok = x[x >= xmin] 

2767 th = el.x2t_headwave(xstretch=(xok-xmin)) + tmin 

2768 return [ 

2769 (x_, self._p[0], t, None) for (x_, t) in zip(xok, th)] 

2770 

2771 else: 

2772 if num.all(x < self._xmin) or num.all(self._xmax < x): 

2773 return [] 

2774 

2775 rp, rx, rt = self.draft_pxt(endgaps) 

2776 

2777 xp = interp(x, rx, rp, 0) 

2778 xt = interp(x, rx, rt, 0) 

2779 

2780 if (rp.size and 

2781 len(xp) == 0 and 

2782 rx[0] == 0.0 and 

2783 any(x == 0.0) and 

2784 rp[0] == 0.0): 

2785 

2786 xp = [(0.0, rp[0])] 

2787 xt = [(0.0, rt[0])] 

2788 

2789 return [ 

2790 (x_, p, t, (rp, rx, rt)) for ((x_, p), (_, t)) in zip(xp, xt)] 

2791 

2792 def __eq__(self, other): 

2793 if len(self.elements) != len(other.elements): 

2794 return False 

2795 

2796 return all(a == b for a, b in zip(self.elements, other.elements)) 

2797 

2798 def __hash__(self): 

2799 return hash( 

2800 tuple(hash(x) for x in self.elements) + 

2801 (self.phase.definition(), )) 

2802 

2803 def __str__(self, p=None, eps=1.): 

2804 x = [] 

2805 start_i = None 

2806 end_i = None 

2807 turn_i = None 

2808 

2809 def append_layers(si, ei, ti): 

2810 if si == ei and (ti is None or ti == si): 

2811 x.append('%i' % si) 

2812 else: 

2813 if ti is not None: 

2814 x.append('(%i-%i-%i)' % (si, ti, ei)) 

2815 else: 

2816 x.append('(%i-%i)' % (si, ei)) 

2817 

2818 for el in self.elements: 

2819 if type(el) is Straight: 

2820 if start_i is None: 

2821 start_i = el.layer.ilayer 

2822 if el._direction_in != el._direction_out: 

2823 turn_i = el.layer.ilayer 

2824 end_i = el.layer.ilayer 

2825 

2826 elif isinstance(el, Kink): 

2827 if start_i is not None: 

2828 append_layers(start_i, end_i, turn_i) 

2829 start_i = None 

2830 turn_i = None 

2831 

2832 x.append(str(el)) 

2833 

2834 if start_i is not None: 

2835 append_layers(start_i, end_i, turn_i) 

2836 

2837 su = '(%s)' % self.used_phase(p=p, eps=eps).used_repr() 

2838 

2839 return '%-15s %-17s %s' % (self.phase.definition(), su, ''.join(x)) 

2840 

2841 def critical_pstart(self, endgaps): 

2842 ''' 

2843 Get critical ray parameter for source depth choice. 

2844 ''' 

2845 

2846 return self.first_straight().critical_p_in(endgaps) 

2847 

2848 def critical_pstop(self, endgaps): 

2849 ''' 

2850 Get critical ray parameter for receiver depth choice. 

2851 ''' 

2852 

2853 return self.last_straight().critical_p_out(endgaps) 

2854 

2855 def ranges(self, endgaps): 

2856 ''' 

2857 Get valid ranges of ray parameter, distance, and traveltime. 

2858 ''' 

2859 p, x, t = self.draft_pxt(endgaps) 

2860 return p.min(), p.max(), x.min(), x.max(), t.min(), t.max() 

2861 

2862 def describe(self, endgaps=None, as_degrees=False): 

2863 ''' 

2864 Get textual representation. 

2865 ''' 

2866 

2867 self._analyse() 

2868 

2869 if as_degrees: 

2870 xunit = 'deg' 

2871 xfact = 1. 

2872 else: 

2873 xunit = 'km' 

2874 xfact = d2m/km 

2875 

2876 sg = ''' Ranges for all depths in source and receiver layers: 

2877 - x [%g, %g] %s 

2878 - t [%g, %g] s 

2879 - p [%g, %g] s/deg 

2880''' % ( 

2881 self._xmin*xfact, 

2882 self._xmax*xfact, 

2883 xunit, 

2884 self._tmin, 

2885 self._tmax, 

2886 self._pmin/r2d, 

2887 self._pmax/r2d) 

2888 

2889 if endgaps is not None: 

2890 pmin, pmax, xmin, xmax, tmin, tmax = self.ranges(endgaps) 

2891 ss = ''' Ranges for given source and receiver depths: 

2892\n - x [%g, %g] %s 

2893\n - t [%g, %g] s 

2894\n - p [%g, %g] s/deg 

2895\n''' % (xmin*xfact, xmax*xfact, xunit, tmin, tmax, pmin/r2d, pmax/r2d) 

2896 

2897 else: 

2898 ss = '' 

2899 

2900 return '%s\n' % self + ss + sg 

2901 

2902 

2903class RefineFailed(CakeError): 

2904 pass 

2905 

2906 

2907class Ray(object): 

2908 ''' 

2909 Representation of a ray with a specific (path, ray parameter, distance, 

2910 arrival time) choice. 

2911 

2912 **Attributes:** 

2913 

2914 .. py:attribute:: path 

2915 

2916 :py:class:`RayPath` object containing complete propagation history. 

2917 

2918 .. py:attribute:: p 

2919 

2920 Ray parameter (spherical) [s/rad] 

2921 

2922 .. py:attribute:: x 

2923 

2924 Radial distance [deg] 

2925 

2926 .. py:attribute:: t 

2927 

2928 Traveltime [s] 

2929 

2930 .. py:attribute:: endgaps 

2931 

2932 Needed for source/receiver depth adjustments in many 

2933 :py:class:`RayPath` methods. 

2934 ''' 

2935 

2936 def __init__(self, path, p, x, t, endgaps, draft_pxt): 

2937 self.path = path 

2938 self.p = p 

2939 self.x = x 

2940 self.t = t 

2941 self.endgaps = endgaps 

2942 self.draft_pxt = draft_pxt 

2943 

2944 def given_phase(self): 

2945 ''' 

2946 Get phase definition which was used to create the ray. 

2947 

2948 :returns: :py:class:`PhaseDef` object 

2949 ''' 

2950 

2951 return self.path.phase 

2952 

2953 def used_phase(self): 

2954 ''' 

2955 Compute phase definition from propagation path. 

2956 

2957 :returns: :py:class:`PhaseDef` object 

2958 ''' 

2959 

2960 return self.path.used_phase(self.p) 

2961 

2962 def refine(self): 

2963 if self.path._is_headwave: 

2964 return 

2965 

2966 if self.t == 0.0 and self.p == 0.0 and self.x == 0.0: 

2967 return 

2968 

2969 cp, cx, ct = self.draft_pxt 

2970 ip = num.searchsorted(cp, self.p) 

2971 if not (0 < ip < cp.size): 

2972 raise RefineFailed() 

2973 

2974 pl, ph = cp[ip-1], cp[ip] 

2975 p_to_t = {} 

2976 i = [0] 

2977 

2978 def f(p): 

2979 i[0] += 1 

2980 x, t = self.path.xt(p, self.endgaps) 

2981 p_to_t[p] = t 

2982 return self.x - x 

2983 

2984 try: 

2985 self.p = brentq(f, pl, ph) 

2986 self.t = p_to_t[self.p] 

2987 

2988 except ValueError: 

2989 raise RefineFailed() 

2990 

2991 def t_and_attributes(self, attributes): 

2992 d = { 

2993 'takeoff_angle': lambda: self.takeoff_angle(), 

2994 'incidence_angle': lambda: self.incidence_angle(), 

2995 't': lambda: self.t, 

2996 'p': lambda: self.p} 

2997 

2998 return tuple(d[k]() for k in ['t'] + attributes) 

2999 

3000 def takeoff_angle(self): 

3001 ''' 

3002 Get takeoff angle of ray. 

3003 

3004 The angle is returned in [degrees]. 

3005 ''' 

3006 

3007 return self.path.first_straight().angle_in(self.p, self.endgaps) 

3008 

3009 def incidence_angle(self): 

3010 ''' 

3011 Get incidence angle of ray. 

3012 

3013 The angle is returned in [degrees]. 

3014 ''' 

3015 

3016 return self.path.last_straight().angle_out(self.p, self.endgaps) 

3017 

3018 def efficiency(self): 

3019 ''' 

3020 Get conversion/reflection efficiency of the ray. 

3021 

3022 A value between 0 and 1 is returned, reflecting the relative amount of 

3023 energy which is transmitted along the ray and not lost by reflections 

3024 or conversions. 

3025 ''' 

3026 

3027 return self.path.efficiency(self.p) 

3028 

3029 def spreading(self): 

3030 ''' 

3031 Get geometrical spreading factor. 

3032 ''' 

3033 

3034 return self.path.spreading(self.p, self.endgaps) 

3035 

3036 def surface_sphere(self): 

3037 x1, y1 = 0., earthradius - self.endgaps[0] 

3038 r2 = earthradius - self.endgaps[1] 

3039 x2, y2 = r2*math.sin(self.x*d2r), r2*math.cos(self.x*d2r) 

3040 return ((x2-x1)**2 + (y2-y1)**2)*4.0*math.pi 

3041 

3042 def zxt_path_subdivided(self, points_per_straight=20): 

3043 ''' 

3044 Get geometrical representation of ray path. 

3045 

3046 Three arrays (depth, distance, time) with points on the ray's path of 

3047 propagation are returned. The number of points which are used in each 

3048 ray segment (passage through one layer) may be controlled by the 

3049 ``points_per_straight`` parameter. 

3050 ''' 

3051 return self.path.zxt_path_subdivided( 

3052 num.atleast_1d(self.p), self.endgaps, 

3053 points_per_straight=points_per_straight, 

3054 x_for_headwave=num.atleast_1d(self.x)) 

3055 

3056 def __str__(self, as_degrees=False): 

3057 if as_degrees: 

3058 sd = '%6.3g deg' % self.x 

3059 else: 

3060 sd = '%7.5g km' % (self.x*(d2r*earthradius/km)) 

3061 

3062 return '%7.5g s/deg %s %6.4g s %5.1f %5.1f %3.0f%% %3.0f%% %s' % ( 

3063 self.p/r2d, 

3064 sd, 

3065 self.t, 

3066 self.takeoff_angle(), 

3067 self.incidence_angle(), 

3068 100*self.efficiency(), 

3069 100*self.spreading()*self.surface_sphere(), 

3070 self.path.__str__(p=self.p)) 

3071 

3072 

3073def anything_to_crust2_profile(crust2_profile): 

3074 from pyrocko.dataset import crust2x2 

3075 if isinstance(crust2_profile, tuple): 

3076 lat, lon = [float(x) for x in crust2_profile] 

3077 return crust2x2.get_profile(lat, lon) 

3078 elif isinstance(crust2_profile, str): 

3079 return crust2x2.get_profile(crust2_profile) 

3080 elif isinstance(crust2_profile, crust2x2.Crust2Profile): 

3081 return crust2_profile 

3082 else: 

3083 assert False, 'crust2_profile must be (lat, lon) a profile ' \ 

3084 'key or a crust2x2 Profile object)' 

3085 

3086 

3087class DiscontinuityNotFound(CakeError): 

3088 def __init__(self, depth_or_name): 

3089 CakeError.__init__(self) 

3090 self.depth_or_name = depth_or_name 

3091 

3092 def __str__(self): 

3093 return 'Cannot find discontinuity from given depth or name: %s' % \ 

3094 self.depth_or_name 

3095 

3096 

3097class LayeredModelError(CakeError): 

3098 pass 

3099 

3100 

3101class LayeredModel(object): 

3102 ''' 

3103 Representation of a layer cake model. 

3104 

3105 There are several ways to initialize an instance of this class. 

3106 

3107 1. Use the module function :py:func:`load_model` to read a model from a 

3108 file. 

3109 2. Create an empty model with the default constructor and append layers and 

3110 discontinuities with the :py:meth:`append` method (from top to bottom). 

3111 3. Use the constructor :py:meth:`LayeredModel.from_scanlines`, to 

3112 automatically create the :py:class:`Layer` and :py:class:`Discontinuity` 

3113 objects from a given velocity profile. 

3114 

3115 An earth model is represented by as stack of :py:class:`Layer` and 

3116 :py:class:`Discontinuity` objects. The method :py:meth:`arrivals` returns 

3117 :py:class:`Ray` objects which may be e.g. queried for arrival times of 

3118 specific phases. Each ray is associated with a :py:class:`RayPath` object. 

3119 Ray objects share common ray paths if they have the same 

3120 conversion/reflection/propagation history. Creating the ray path objects is 

3121 relatively expensive (this is done in :py:meth:`gather_paths`), but they 

3122 are cached for reuse in successive invocations. 

3123 ''' 

3124 

3125 def __init__(self): 

3126 self._surface_material = None 

3127 self._elements = [] 

3128 self.nlayers = 0 

3129 self._np = 10000 

3130 self._pdepth = 5 

3131 self._pathcache = {} 

3132 

3133 def copy_with_elevation(self, elevation): 

3134 ''' 

3135 Get copy of the model with surface layer stretched to given elevation. 

3136 

3137 :param elevation: new surface elevation in [m] 

3138 

3139 Elevation is positiv upward, contrary to the layered models downward 

3140 `z` axis. 

3141 ''' 

3142 

3143 c = copy.deepcopy(self) 

3144 c._pathcache = {} 

3145 surface = c._elements[0] 

3146 toplayer = c._elements[1] 

3147 

3148 assert toplayer.zbot > -elevation 

3149 

3150 surface.z = -elevation 

3151 c._elements[1] = toplayer.copy(ztop=-elevation) 

3152 c._elements[1].ilayer = 0 

3153 return c 

3154 

3155 def zeq(self, z1, z2): 

3156 return abs(z1-z2) < ZEPS 

3157 

3158 def append(self, element): 

3159 ''' 

3160 Add a layer or discontinuity at bottom of model. 

3161 

3162 :param element: object of subclass of :py:class:`Layer` or 

3163 :py:class:`Discontinuity`. 

3164 ''' 

3165 

3166 if isinstance(element, Layer): 

3167 if element.zbot >= earthradius: 

3168 element.zbot = earthradius - 1. 

3169 

3170 if element.ztop >= earthradius: 

3171 raise CakeError('Layer deeper than earthradius') 

3172 

3173 element.ilayer = self.nlayers 

3174 self.nlayers += 1 

3175 

3176 self._elements.append(element) 

3177 

3178 def elements(self, direction=DOWN): 

3179 ''' 

3180 Iterate over all elements of the model. 

3181 

3182 :param direction: direction of traversal :py:const:`DOWN` or 

3183 :py:const:`UP`. 

3184 

3185 Objects derived from the :py:class:`Discontinuity` and 

3186 :py:class:`Layer` classes are yielded. 

3187 ''' 

3188 

3189 if direction == DOWN: 

3190 return iter(self._elements) 

3191 else: 

3192 return reversed(self._elements) 

3193 

3194 def layers(self, direction=DOWN): 

3195 ''' 

3196 Iterate over all layers of model. 

3197 

3198 :param direction: direction of traversal :py:const:`DOWN` or 

3199 :py:const:`UP`. 

3200 

3201 Objects derived from the :py:class:`Layer` class are yielded. 

3202 ''' 

3203 

3204 if direction == DOWN: 

3205 return (el for el in self._elements if isinstance(el, Layer)) 

3206 else: 

3207 return ( 

3208 el for el in reversed(self._elements) if isinstance(el, Layer)) 

3209 

3210 def layer(self, z, direction=DOWN): 

3211 ''' 

3212 Get layer for given depth. 

3213 

3214 :param z: depth [m] 

3215 :param direction: direction of traversal :py:const:`DOWN` or 

3216 :py:const:`UP`. 

3217 

3218 Returns first layer which touches depth ``z`` (tolerant at boundaries). 

3219 ''' 

3220 

3221 for layer in self.layers(direction): 

3222 if layer.contains(z): 

3223 return layer 

3224 else: 

3225 raise CakeError('Failed extracting layer at depth z=%s' % z) 

3226 

3227 def walker(self): 

3228 return Walker(self._elements) 

3229 

3230 def material(self, z, direction=DOWN): 

3231 ''' 

3232 Get material at given depth. 

3233 

3234 :param z: depth [m] 

3235 :param direction: direction of traversal :py:const:`DOWN` or 

3236 :py:const:`UP` 

3237 :returns: object of type :py:class:`Material` 

3238 

3239 If given depth ``z`` happens to be at an interface, the material of the 

3240 first layer with respect to the the traversal ordering is returned. 

3241 ''' 

3242 

3243 lyr = self.layer(z, direction) 

3244 return lyr.material(z) 

3245 

3246 def discontinuities(self): 

3247 ''' 

3248 Iterate over all discontinuities of the model.''' 

3249 

3250 return (el for el in self._elements if isinstance(el, Discontinuity)) 

3251 

3252 def discontinuity(self, name_or_z): 

3253 ''' 

3254 Get discontinuity by name or depth. 

3255 

3256 :param name_or_z: name of discontinuity or depth [m] as float value 

3257 ''' 

3258 

3259 if isinstance(name_or_z, float): 

3260 candi = sorted( 

3261 self.discontinuities(), key=lambda i: abs(i.z-name_or_z)) 

3262 else: 

3263 candi = [i for i in self.discontinuities() if i.name == name_or_z] 

3264 

3265 if not candi: 

3266 raise DiscontinuityNotFound(name_or_z) 

3267 

3268 return candi[0] 

3269 

3270 def adapt_phase(self, phase): 

3271 ''' 

3272 Adapt a phase definition for use with this model. 

3273 

3274 This returns a copy of the phase definition, where named 

3275 discontinuities are replaced with the actual depth of these, as defined 

3276 in the model. 

3277 ''' 

3278 

3279 phase = phase.copy() 

3280 for knee in phase.knees(): 

3281 if knee.depth != 'surface': 

3282 knee.depth = self.discontinuity(knee.depth).z 

3283 for leg in phase.legs(): 

3284 if leg.depthmax is not None and isinstance(leg.depthmax, str): 

3285 leg.depthmax = self.discontinuity(leg.depthmax).z 

3286 

3287 return phase 

3288 

3289 def path(self, p, phase, layer_start, layer_stop): 

3290 ''' 

3291 Get ray path for given combination of ray parameter, phase definition, 

3292 source and receiver layers. 

3293 

3294 :param p: ray parameter (spherical) [s/rad] 

3295 :param phase: phase definition (:py:class:`PhaseDef` object) 

3296 :param layer_start: layer with source 

3297 :param layer_stop: layer with receiver 

3298 :returns: :py:class:`RayPath` object 

3299 

3300 If it is not possible to find a solution, an exception of type 

3301 :py:exc:`NotPhaseConform`, :py:exc:`MinDepthReached`, 

3302 :py:exc:`MaxDepthReached`, :py:exc:`CannotPropagate`, 

3303 :py:exc:`BottomReached` or :py:exc:`SurfaceReached` is raised. 

3304 ''' 

3305 

3306 phase = self.adapt_phase(phase) 

3307 knees = phase.knees() 

3308 legs = phase.legs() 

3309 next_knee = next_or_none(knees) 

3310 leg = next_or_none(legs) 

3311 assert leg is not None 

3312 

3313 direction = leg.departure 

3314 direction_stop = phase.direction_stop() 

3315 mode = leg.mode 

3316 mode_stop = phase.last_leg().mode 

3317 

3318 walker = self.walker() 

3319 walker.goto_layer(layer_start) 

3320 current = walker.current() 

3321 

3322 ttop, tbot = current.tests(p, mode) 

3323 if not ttop and not tbot: 

3324 raise CannotPropagate(direction, current.ilayer) 

3325 

3326 if (direction == DOWN and not ttop) or (direction == UP and not tbot): 

3327 direction = -direction 

3328 

3329 path = RayPath(phase) 

3330 trapdetect = set() 

3331 while True: 

3332 at_layer = isinstance(current, Layer) 

3333 at_discontinuity = isinstance(current, Discontinuity) 

3334 

3335 # detect trapped wave 

3336 k = (id(next_knee), id(current), direction, mode) 

3337 if k in trapdetect: 

3338 raise Trapped() 

3339 

3340 trapdetect.add(k) 

3341 

3342 if at_discontinuity: 

3343 oldmode, olddirection = mode, direction 

3344 headwave = False 

3345 if next_knee is not None and next_knee.matches( 

3346 current, mode, direction): 

3347 

3348 headwave = next_knee.headwave 

3349 direction = next_knee.out_direction() 

3350 mode = next_knee.out_mode 

3351 next_knee = next_or_none(knees) 

3352 leg = next(legs) 

3353 

3354 else: # implicit reflection/transmission 

3355 direction = current.propagate(p, mode, direction) 

3356 

3357 if headwave: 

3358 path.set_is_headwave(True) 

3359 

3360 path.append(Kink( 

3361 olddirection, olddirection, oldmode, oldmode, current)) 

3362 

3363 path.append(HeadwaveStraight( 

3364 olddirection, direction, oldmode, current)) 

3365 

3366 path.append(Kink( 

3367 olddirection, direction, oldmode, mode, current)) 

3368 

3369 else: 

3370 path.append(Kink( 

3371 olddirection, direction, oldmode, mode, current)) 

3372 

3373 if at_layer: 

3374 direction_in = direction 

3375 direction = current.propagate(p, mode, direction_in) 

3376 

3377 zturn = None 

3378 if direction_in != direction: 

3379 zturn = current.zturn(p, mode) 

3380 

3381 zmin, zmax = leg.depthmin, leg.depthmax 

3382 if zmin is not None or zmax is not None: 

3383 if direction_in != direction: 

3384 if zmin is not None and zturn <= zmin: 

3385 raise MinDepthReached() 

3386 if zmax is not None and zturn >= zmax: 

3387 raise MaxDepthReached() 

3388 else: 

3389 if zmin is not None and current.ztop <= zmin: 

3390 raise MinDepthReached() 

3391 if zmax is not None and current.zbot >= zmax: 

3392 raise MaxDepthReached() 

3393 

3394 path.append(Straight(direction_in, direction, mode, current)) 

3395 

3396 if next_knee is None and mode == mode_stop and \ 

3397 current is layer_stop: 

3398 

3399 if zturn is None: 

3400 if direction == direction_stop: 

3401 break 

3402 else: 

3403 break 

3404 

3405 walker.go(direction) 

3406 current = walker.current() 

3407 

3408 return path 

3409 

3410 def gather_paths(self, phases=PhaseDef('P'), zstart=0.0, zstop=0.0): 

3411 ''' 

3412 Get all possible ray paths for given source and receiver depths for one 

3413 or more phase definitions. 

3414 

3415 :param phases: a :py:class:`PhaseDef` object or a list of such objects. 

3416 Comma-separated strings and lists of such strings are also accepted 

3417 and are converted to :py:class:`PhaseDef` objects for convenience. 

3418 :param zstart: source depth [m] 

3419 :param zstop: receiver depth [m] 

3420 :returns: a list of :py:class:`RayPath` objects 

3421 

3422 Results of this method are cached internally. Cached results are 

3423 returned, when a given combination of source layer, receiver layer and 

3424 phase definition has been used before. 

3425 ''' 

3426 

3427 eps = 1e-7 # num.finfo(float).eps * 1000. 

3428 

3429 phases = to_phase_defs(phases) 

3430 

3431 paths = [] 

3432 for phase in phases: 

3433 

3434 layer_start = self.layer(zstart, -phase.direction_start()) 

3435 layer_stop = self.layer(zstop, phase.direction_stop()) 

3436 

3437 pathcachekey = (phase.definition(), layer_start, layer_stop) 

3438 

3439 if pathcachekey in self._pathcache: 

3440 phase_paths = self._pathcache[pathcachekey] 

3441 else: 

3442 hwknee = phase.headwave_knee() 

3443 if hwknee: 

3444 name_or_z = hwknee.depth 

3445 interface = self.discontinuity(name_or_z) 

3446 mode = hwknee.in_mode 

3447 in_direction = hwknee.direction 

3448 

3449 pabove, pbelow = interface.critical_ps(mode) 

3450 

3451 p = min_not_none(pabove, pbelow) 

3452 

3453 # diffracted wave: 

3454 if in_direction == DOWN and ( 

3455 pbelow is None or pbelow >= pabove): 

3456 

3457 p *= (1.0 - eps) 

3458 

3459 path = self.path(p, phase, layer_start, layer_stop) 

3460 path.set_prange(p, p, 1.) 

3461 

3462 phase_paths = [path] 

3463 

3464 else: 

3465 try: 

3466 pmax_start = max([ 

3467 radius(z)/layer_start.v(phase.first_leg().mode, z) 

3468 for z in (layer_start.ztop, layer_start.zbot)]) 

3469 

3470 pmax_stop = max([ 

3471 radius(z)/layer_stop.v(phase.last_leg().mode, z) 

3472 for z in (layer_stop.ztop, layer_stop.zbot)]) 

3473 

3474 pmax = min(pmax_start, pmax_stop) 

3475 

3476 pedges = [0.] 

3477 for layer in self.layers(): 

3478 for z in (layer.ztop, layer.zbot): 

3479 for mode in (P, S): 

3480 for eps2 in [eps]: 

3481 v = layer.v(mode, z) 

3482 if v != 0.0: 

3483 p = radius(z)/v 

3484 if p <= pmax: 

3485 pedges.append(p*(1.0-eps2)) 

3486 pedges.append(p) 

3487 pedges.append(p*(1.0+eps2)) 

3488 

3489 pedges = num.unique(sorted(pedges)) 

3490 

3491 phase_paths = {} 

3492 cached = {} 

3493 counter = [0] 

3494 

3495 def p_to_path(p): 

3496 if p in cached: 

3497 return cached[p] 

3498 

3499 try: 

3500 counter[0] += 1 

3501 path = self.path( 

3502 p, phase, layer_start, layer_stop) 

3503 

3504 if path not in phase_paths: 

3505 phase_paths[path] = [] 

3506 

3507 phase_paths[path].append(p) 

3508 

3509 except PathFailed: 

3510 path = None 

3511 

3512 cached[p] = path 

3513 return path 

3514 

3515 def recurse(pmin, pmax, i=0): 

3516 if i > self._pdepth: 

3517 return 

3518 path1 = p_to_path(pmin) 

3519 path2 = p_to_path(pmax) 

3520 if path1 is None and path2 is None and i > 0: 

3521 return 

3522 if path1 is None or path2 is None or \ 

3523 hash(path1) != hash(path2): 

3524 

3525 recurse(pmin, (pmin+pmax)/2., i+1) 

3526 recurse((pmin+pmax)/2., pmax, i+1) 

3527 

3528 for (pl, ph) in zip(pedges[:-1], pedges[1:]): 

3529 recurse(pl, ph) 

3530 

3531 for path, ps in phase_paths.items(): 

3532 path.set_prange( 

3533 min(ps), max(ps), pmax/(self._np-1)) 

3534 

3535 phase_paths = list(phase_paths.keys()) 

3536 

3537 except ZeroDivisionError: 

3538 phase_paths = [] 

3539 

3540 self._pathcache[pathcachekey] = phase_paths 

3541 

3542 paths.extend(phase_paths) 

3543 

3544 paths.sort(key=lambda x: x.pmin()) 

3545 return paths 

3546 

3547 def arrivals( 

3548 self, 

3549 distances=[], 

3550 phases=PhaseDef('P'), 

3551 zstart=0.0, 

3552 zstop=0.0, 

3553 refine=True): 

3554 

3555 ''' 

3556 Compute rays and traveltimes for given distances. 

3557 

3558 :param distances: list or array of distances [deg] 

3559 :param phases: a :py:class:`PhaseDef` object or a list of such objects. 

3560 Comma-separated strings and lists of such strings are also accepted 

3561 and are converted to :py:class:`PhaseDef` objects for convenience. 

3562 :param zstart: source depth [m] 

3563 :param zstop: receiver depth [m] 

3564 :param refine: bool flag, whether to use bisectioning to improve 

3565 (p, x, t) estimated from interpolation 

3566 :returns: a list of :py:class:`Ray` objects, sorted by 

3567 (distance, arrival time) 

3568 ''' 

3569 

3570 distances = num.asarray(distances, dtype=float) 

3571 

3572 arrivals = [] 

3573 for path in self.gather_paths(phases, zstart=zstart, zstop=zstop): 

3574 

3575 endgaps = path.endgaps(zstart, zstop) 

3576 for x, p, t, draft_pxt in path.interpolate_x2pt_linear( 

3577 distances, endgaps): 

3578 

3579 arrivals.append(Ray(path, p, x, t, endgaps, draft_pxt)) 

3580 

3581 if refine: 

3582 refined = [] 

3583 for ray in arrivals: 

3584 

3585 if ray.path._is_headwave: 

3586 refined.append(ray) 

3587 

3588 try: 

3589 ray.refine() 

3590 refined.append(ray) 

3591 

3592 except RefineFailed: 

3593 pass 

3594 

3595 arrivals = refined 

3596 

3597 arrivals.sort(key=lambda x: (x.x, x.t)) 

3598 return arrivals 

3599 

3600 @classmethod 

3601 def from_scanlines(cls, producer): 

3602 ''' 

3603 Create layer cake model from sequence of materials at depths. 

3604 

3605 :param producer: iterable yielding (depth, material, name) tuples 

3606 

3607 Creates a new :py:class:`LayeredModel` object and uses its 

3608 :py:meth:`append` method to add layers and discontinuities as needed. 

3609 ''' 

3610 

3611 self = cls() 

3612 for z, material, name in producer: 

3613 

3614 if not self._elements: 

3615 self.append(Surface(z, material)) 

3616 else: 

3617 element = self._elements[-1] 

3618 if self.zeq(element.zbot, z): 

3619 assert isinstance(element, Layer) 

3620 self.append( 

3621 Interface(z, element.mbot, material, name=name)) 

3622 

3623 else: 

3624 if isinstance(element, Discontinuity): 

3625 ztop = element.z 

3626 mtop = element.mbelow 

3627 elif isinstance(element, Layer): 

3628 ztop = element.zbot 

3629 mtop = element.mbot 

3630 

3631 if mtop == material: 

3632 layer = HomogeneousLayer( 

3633 ztop, z, material, name=name) 

3634 else: 

3635 layer = GradientLayer( 

3636 ztop, z, mtop, material, name=name) 

3637 

3638 self.append(layer) 

3639 

3640 return self 

3641 

3642 def to_scanlines(self, get_burgers=False): 

3643 def fmt(z, m): 

3644 if not m._has_default_burgers() or get_burgers: 

3645 return (z, m.vp, m.vs, m.rho, m.qp, m.qs, 

3646 m.burger_eta1, m.burger_eta2, m.burger_valpha) 

3647 return (z, m.vp, m.vs, m.rho, m.qp, m.qs) 

3648 

3649 last = None 

3650 lines = [] 

3651 for element in self.elements(): 

3652 if isinstance(element, Layer): 

3653 if not isinstance(last, Layer): 

3654 lines.append(fmt(element.ztop, element.mtop)) 

3655 

3656 lines.append(fmt(element.zbot, element.mbot)) 

3657 

3658 last = element 

3659 

3660 if not isinstance(last, Layer): 

3661 lines.append(fmt(last.z, last.mbelow)) 

3662 

3663 return lines 

3664 

3665 def iter_material_parameter(self, get): 

3666 assert get in ('vp', 'vs', 'rho', 'qp', 'qs', 'z') 

3667 if get == 'z': 

3668 for layer in self.layers(): 

3669 yield layer.ztop 

3670 yield layer.zbot 

3671 else: 

3672 getter = operator.attrgetter(get) 

3673 for layer in self.layers(): 

3674 yield getter(layer.mtop) 

3675 yield getter(layer.mbot) 

3676 

3677 def profile(self, get): 

3678 ''' 

3679 Get parameter profile along depth of the earthmodel. 

3680 

3681 :param get: property to be queried ( 

3682 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``) 

3683 :type get: str 

3684 ''' 

3685 

3686 return num.array(list(self.iter_material_parameter(get))) 

3687 

3688 def min(self, get='vp'): 

3689 ''' 

3690 Find minimum value of a material property or depth. 

3691 

3692 :param get: property to be queried ( 

3693 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``) 

3694 ''' 

3695 

3696 return min(self.iter_material_parameter(get)) 

3697 

3698 def max(self, get='vp'): 

3699 ''' 

3700 Find maximum value of a material property or depth. 

3701 

3702 :param get: property to be queried ( 

3703 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, ``'qs'``, or ``'z'``) 

3704 ''' 

3705 

3706 return max(self.iter_material_parameter(get)) 

3707 

3708 def simplify_layers(self, layers, max_rel_error=0.001): 

3709 if len(layers) <= 1: 

3710 return layers 

3711 

3712 ztop = layers[0].ztop 

3713 zbot = layers[-1].zbot 

3714 zorigs = [layer.ztop for layer in layers] 

3715 zorigs.append(zbot) 

3716 zs = num.linspace(ztop, zbot, 100) 

3717 data = [] 

3718 for z in zs: 

3719 if z == ztop: 

3720 direction = UP 

3721 else: 

3722 direction = DOWN 

3723 

3724 mat = self.material(z, direction) 

3725 data.append(mat.astuple()) 

3726 

3727 data = num.array(data, dtype=float) 

3728 data_means = num.mean(data, axis=0) 

3729 nmax = len(layers) // 2 

3730 accept = False 

3731 

3732 zcut_best = [] 

3733 for n in range(1, nmax+1): 

3734 ncutintervals = 20 

3735 zdelta = (zbot-ztop)/ncutintervals 

3736 if n == 2: 

3737 zcuts = [ 

3738 [ztop, ztop + i*zdelta, zbot] 

3739 for i in range(1, ncutintervals)] 

3740 elif n == 3: 

3741 zcuts = [] 

3742 for j in range(1, ncutintervals): 

3743 for i in range(j+1, ncutintervals): 

3744 zcuts.append( 

3745 [ztop, ztop + j*zdelta, ztop + i*zdelta, zbot]) 

3746 else: 

3747 zcuts = [] 

3748 zcuts.append(num.linspace(ztop, zbot, n+1)) 

3749 if zcut_best: 

3750 zcuts.append(sorted(num.linspace( 

3751 ztop, zbot, n).tolist() + zcut_best[1])) 

3752 zcuts.append(sorted(num.linspace( 

3753 ztop, zbot, n-1).tolist() + zcut_best[2])) 

3754 

3755 best = None 

3756 for icut, zcut in enumerate(zcuts): 

3757 rel_par_errors = num.zeros(5) 

3758 mpar_nodes = num.zeros((n+1, 5)) 

3759 

3760 for ipar in range(5): 

3761 znodes, vnodes, error_rms = util.polylinefit( 

3762 zs, data[:, ipar], zcut) 

3763 

3764 mpar_nodes[:, ipar] = vnodes 

3765 if data_means[ipar] == 0.0: 

3766 rel_par_errors[ipar] = -1 

3767 else: 

3768 rel_par_errors[ipar] = error_rms/data_means[ipar] 

3769 

3770 rel_error = rel_par_errors.max() 

3771 if best is None or rel_error < best[0]: 

3772 best = (rel_error, zcut, mpar_nodes) 

3773 

3774 rel_error, zcut, mpar_nodes = best 

3775 

3776 zcut_best.append(list(zcut)) 

3777 zcut_best[-1].pop(0) 

3778 zcut_best[-1].pop() 

3779 

3780 if rel_error <= max_rel_error: 

3781 accept = True 

3782 break 

3783 

3784 if not accept: 

3785 return layers 

3786 

3787 rel_error, zcut, mpar_nodes = best 

3788 

3789 material_nodes = [] 

3790 for i in range(n+1): 

3791 material_nodes.append(Material(*mpar_nodes[i, :])) 

3792 

3793 out_layers = [] 

3794 for i in range(n): 

3795 mtop = material_nodes[i] 

3796 mbot = material_nodes[i+1] 

3797 ztop = zcut[i] 

3798 zbot = zcut[i+1] 

3799 if mtop == mbot: 

3800 lyr = HomogeneousLayer(ztop, zbot, mtop) 

3801 else: 

3802 lyr = GradientLayer(ztop, zbot, mtop, mbot) 

3803 

3804 out_layers.append(lyr) 

3805 return out_layers 

3806 

3807 def simplify(self, max_rel_error=0.001): 

3808 ''' 

3809 Get representation of model with lower resolution. 

3810 

3811 Returns an approximation of the model. All discontinuities are kept, 

3812 but layer stacks with continuous model parameters are represented, if 

3813 possible, by a lower number of layers. Piecewise linear functions are 

3814 fitted against the original model parameter's piecewise linear 

3815 functions. Successively larger numbers of layers are tried, until the 

3816 difference to the original model is below ``max_rel_error``. The 

3817 difference is measured as the RMS error of the fit normalized by the 

3818 mean of the input (i.e. the fitted curves should deviate, on average, 

3819 less than 0.1% from the input curves if ``max_rel_error`` = 0.001). 

3820 ''' 

3821 

3822 mod_simple = LayeredModel() 

3823 

3824 glayers = [] 

3825 for element in self.elements(): 

3826 

3827 if isinstance(element, Discontinuity): 

3828 for layer in self.simplify_layers( 

3829 glayers, max_rel_error=max_rel_error): 

3830 

3831 mod_simple.append(layer) 

3832 

3833 glayers = [] 

3834 mod_simple.append(element) 

3835 else: 

3836 glayers.append(element) 

3837 

3838 for layer in self.simplify_layers( 

3839 glayers, max_rel_error=max_rel_error): 

3840 mod_simple.append(layer) 

3841 

3842 return mod_simple 

3843 

3844 def extract(self, depth_min=None, depth_max=None): 

3845 ''' 

3846 Extract :py:class:`LayeredModel` from :py:class:`LayeredModel`. 

3847 

3848 :param depth_min: depth of upper cut or name of :py:class:`Interface` 

3849 :param depth_max: depth of lower cut or name of :py:class:`Interface` 

3850 

3851 Interpolates a :py:class:`GradientLayer` at ``depth_min`` and/or 

3852 ``depth_max``. 

3853 ''' 

3854 

3855 if isinstance(depth_min, str): 

3856 depth_min = self.discontinuity(depth_min).z 

3857 

3858 if isinstance(depth_max, str): 

3859 depth_max = self.discontinuity(depth_max).z 

3860 

3861 mod_extracted = LayeredModel() 

3862 

3863 for element in self.elements(): 

3864 element = element.copy() 

3865 do_append = False 

3866 if (depth_min is None or depth_min <= element.ztop) \ 

3867 and (depth_max is None or depth_max >= element.zbot): 

3868 mod_extracted.append(element) 

3869 continue 

3870 

3871 if depth_min is not None: 

3872 if element.ztop < depth_min and depth_min < element.zbot: 

3873 _, element = element.split(depth_min) 

3874 do_append = True 

3875 

3876 if depth_max is not None: 

3877 if element.zbot > depth_max and depth_max > element.ztop: 

3878 element, _ = element.split(depth_max) 

3879 do_append = True 

3880 

3881 if do_append: 

3882 mod_extracted.append(element) 

3883 

3884 return mod_extracted 

3885 

3886 def replaced_crust(self, crust2_profile=None, crustmod=None): 

3887 if crust2_profile is not None: 

3888 profile = anything_to_crust2_profile(crust2_profile) 

3889 crustmod = LayeredModel.from_scanlines( 

3890 from_crust2x2_profile(profile)) 

3891 

3892 newmod = LayeredModel() 

3893 for element in crustmod.extract(depth_max='moho').elements(): 

3894 if element.name != 'moho': 

3895 newmod.append(element) 

3896 else: 

3897 moho1 = element 

3898 

3899 mod = self.extract(depth_min='moho') 

3900 first = True 

3901 for element in mod.elements(): 

3902 if element.name == 'moho': 

3903 if element.z <= moho1.z: 

3904 mbelow = mod.material(moho1.z, direction=UP) 

3905 else: 

3906 mbelow = element.mbelow 

3907 

3908 moho = Interface(moho1.z, moho1.mabove, mbelow, name='moho') 

3909 newmod.append(moho) 

3910 else: 

3911 if first: 

3912 if isinstance(element, Layer) and element.zbot > moho.z: 

3913 newmod.append(GradientLayer( 

3914 moho.z, 

3915 element.zbot, 

3916 moho.mbelow, 

3917 element.mbot, 

3918 name=element.name)) 

3919 

3920 first = False 

3921 else: 

3922 newmod.append(element) 

3923 return newmod 

3924 

3925 def perturb(self, rstate=None, keep_vp_vs=False, zmax=None, **kwargs): 

3926 ''' 

3927 Create a perturbed variant of the earth model. 

3928 

3929 Randomly change the thickness and material parameters of the earth 

3930 model from a uniform distribution. 

3931 

3932 :param kwargs: Maximum change fraction (e.g. 0.1) of the parameters. 

3933 Name the parameter, prefixed by ``p``. Supported parameters are 

3934 ``ph, pvp, pvs, prho, pqs, pqp``. 

3935 :type kwargs: dict 

3936 :param rstate: Random state to draw from, defaults to ``None`` 

3937 :type rstate: :class:`numpy.random.RandomState` 

3938 :param keep_vp_vs: Keep the Vp/Vs ratio, defaults to False 

3939 :type keep_vp_vs: bool 

3940 

3941 :returns: A new, perturbed earth model 

3942 :rtype: :class:`~pyrocko.cake.LayeredModel` 

3943 

3944 .. code-block :: python 

3945 

3946 perturbed_model = model.perturb(ph=.1, pvp=.05, prho=.1) 

3947 ''' 

3948 

3949 if rstate is None: 

3950 rstate = num.random.RandomState() 

3951 

3952 def factor(k): 

3953 try: 

3954 pval = kwargs[k] 

3955 return rstate.uniform(1.0-pval, 1.0+pval) 

3956 except KeyError: 

3957 return 1.0 

3958 

3959 def perturb_material(mat): 

3960 mat_new = copy.deepcopy(mat) 

3961 for param in kwargs: 

3962 if param != 'ph': 

3963 value = getattr(mat, param[1:]) 

3964 setattr(mat_new, param[1:], value*factor(param)) 

3965 

3966 if keep_vp_vs: 

3967 mat_new.vs = mat_new.vp * (mat.vs / mat.vp) 

3968 

3969 return mat_new 

3970 

3971 model_new = LayeredModel() 

3972 elements = list(self.elements()) 

3973 assert isinstance(elements[0], Surface) 

3974 z = elements[0].z 

3975 mat = None 

3976 for element in elements: 

3977 element_new = copy.deepcopy(element) 

3978 if isinstance(element_new, Discontinuity): 

3979 element_new.z = z 

3980 

3981 elif isinstance(element_new, Layer): 

3982 element_new.ztop = z 

3983 if zmax is None or element.zbot < zmax: 

3984 z += (element_new.zbot-element_new.ztop) * factor('ph') 

3985 else: 

3986 z += (element_new.zbot-element_new.ztop) 

3987 element_new.zbot = z 

3988 

3989 if mat is None or element.mtop != mat: 

3990 mat = element.mtop 

3991 if zmax is None or element.ztop < zmax: 

3992 mat_new = perturb_material(mat) 

3993 else: 

3994 mat_new = copy.deepcopy(mat) 

3995 

3996 element_new.mtop = mat_new 

3997 

3998 if element.mbot != mat: 

3999 mat = element.mbot 

4000 if zmax is None or element.zbot < zmax: 

4001 mat_new = perturb_material(mat) 

4002 else: 

4003 mat_new = copy.deepcopy(mat) 

4004 

4005 element_new.mbot = mat_new 

4006 

4007 model_new.append(element_new) 

4008 

4009 return model_new 

4010 

4011 def require_homogeneous(self): 

4012 elements = list(self.elements()) 

4013 

4014 if len(elements) != 2: 

4015 raise LayeredModelError( 

4016 'Homogeneous model required but found more than one layer in ' 

4017 'earthmodel.') 

4018 

4019 if not isinstance(elements[1], HomogeneousLayer): 

4020 raise LayeredModelError( 

4021 'Homogeneous model required but element #1 is not of type ' 

4022 'HomogeneousLayer.') 

4023 

4024 return elements[1].m 

4025 

4026 def __str__(self): 

4027 return '\n'.join(str(element) for element in self._elements) 

4028 

4029 

4030def read_hyposat_model(fn): 

4031 ''' 

4032 Reader for HYPOSAT earth model files. 

4033 

4034 To be used as producer in :py:meth:`LayeredModel.from_scanlines`. 

4035 

4036 Interface names are translated as follows: ``'MOHO'`` -> ``'moho'``, 

4037 ``'CONR'`` -> ``'conrad'`` 

4038 ''' 

4039 

4040 with open(fn, 'r') as f: 

4041 translate = {'MOHO': 'moho', 'CONR': 'conrad'} 

4042 lname = None 

4043 for iline, line in enumerate(f): 

4044 if iline == 0: 

4045 continue 

4046 

4047 z, vp, vs, name = util.unpack_fixed('f10, f10, f10, a4', line) 

4048 if not name: 

4049 name = None 

4050 material = Material(vp*1000., vs*1000.) 

4051 

4052 tname = translate.get(lname, lname) 

4053 yield z*1000., material, tname 

4054 

4055 lname = name 

4056 

4057 

4058def read_nd_model(fn): 

4059 ''' 

4060 Reader for TauP style '.nd' (named discontinuity) files. 

4061 

4062 To be used as producer in :py:meth:`LayeredModel.from_scanlines`. 

4063 

4064 Interface names are translated as follows: ``'mantle'`` -> ``'moho'``, 

4065 ``'outer-core'`` -> ``'cmb'``, ``'inner-core'`` -> ``'icb'``. 

4066 

4067 The format has been modified to include Burgers materials parameters in 

4068 columns 7 (burger_eta1), 8 (burger_eta2) and 9. eta(3). 

4069 ''' 

4070 with open(fn, 'r') as f: 

4071 for x in read_nd_model_fh(f): 

4072 yield x 

4073 

4074 

4075def read_nd_model_str(s): 

4076 f = StringIO(s) 

4077 for x in read_nd_model_fh(f): 

4078 yield x 

4079 f.close() 

4080 

4081 

4082def read_nd_model_fh(f): 

4083 translate = {'mantle': 'moho', 'outer-core': 'cmb', 'inner-core': 'icb'} 

4084 name = None 

4085 for line in f: 

4086 toks = line.split() 

4087 if len(toks) == 9 or len(toks) == 6 or len(toks) == 4: 

4088 z, vp, vs, rho = [float(x) for x in toks[:4]] 

4089 qp, qs = None, None 

4090 burgers = None 

4091 if len(toks) == 6 or len(toks) == 9: 

4092 qp, qs = [float(x) for x in toks[4:6]] 

4093 if len(toks) == 9: 

4094 burgers = \ 

4095 [float(x) for x in toks[6:]] 

4096 

4097 material = Material( 

4098 vp*1000., vs*1000., rho*1000., qp, qs, 

4099 burgers=burgers) 

4100 

4101 yield z*1000., material, name 

4102 name = None 

4103 elif len(toks) == 1: 

4104 name = translate.get(toks[0], toks[0]) 

4105 

4106 f.close() 

4107 

4108 

4109def from_crust2x2_profile(profile, depthmantle=50000): 

4110 from pyrocko.dataset import crust2x2 

4111 

4112 default_qp_qs = { 

4113 'soft sed.': (50., 50.), 

4114 'hard sed.': (200., 200.), 

4115 'upper crust': (600., 400.), 

4116 } 

4117 

4118 z = 0. 

4119 for i in range(8): 

4120 dz, vp, vs, rho = profile.get_layer(i) 

4121 name = crust2x2.Crust2Profile.layer_names[i] 

4122 if name in default_qp_qs: 

4123 qp, qs = default_qp_qs[name] 

4124 else: 

4125 qp, qs = None, None 

4126 

4127 material = Material(vp, vs, rho, qp, qs) 

4128 iname = None 

4129 if i == 7: 

4130 iname = 'moho' 

4131 if dz != 0.0: 

4132 yield z, material, iname 

4133 if i != 7: 

4134 yield z+dz, material, name 

4135 else: 

4136 yield z+depthmantle, material, name 

4137 

4138 z += dz 

4139 

4140 

4141def write_nd_model_fh(mod, fh): 

4142 def fmt(z, mat): 

4143 rstr = ' '.join( 

4144 util.gform(x, 4) 

4145 for x in ( 

4146 z/1000., 

4147 mat.vp/1000., 

4148 mat.vs/1000., 

4149 mat.rho/1000., 

4150 mat.qp, mat.qs)) 

4151 if not mat._has_default_burgers(): 

4152 rstr += ' '.join( 

4153 util.gform(x, 4) 

4154 for x in ( 

4155 mat.burger_eta1, 

4156 mat.burger_eta2, 

4157 mat.burger_valpha)) 

4158 return rstr.rstrip() + '\n' 

4159 

4160 translate = { 

4161 'moho': 'mantle', 

4162 'cmb': 'outer-core', 

4163 'icb': 'inner-core'} 

4164 

4165 last = None 

4166 for element in mod.elements(): 

4167 if isinstance(element, Interface): 

4168 if element.name is not None: 

4169 n = translate.get(element.name, element.name) 

4170 fh.write('%s\n' % n) 

4171 

4172 elif isinstance(element, Layer): 

4173 if not isinstance(last, Layer): 

4174 fh.write(fmt(element.ztop, element.mtop)) 

4175 

4176 fh.write(fmt(element.zbot, element.mbot)) 

4177 

4178 last = element 

4179 

4180 if not isinstance(last, Layer): 

4181 fh.write(fmt(last.z, last.mbelow)) 

4182 

4183 

4184def write_nd_model_str(mod): 

4185 f = StringIO() 

4186 write_nd_model_fh(mod, f) 

4187 return f.getvalue() 

4188 

4189 

4190def write_nd_model(mod, fn): 

4191 with open(fn, 'w') as f: 

4192 write_nd_model_fh(mod, f) 

4193 

4194 

4195def builtin_models(): 

4196 return sorted([ 

4197 os.path.splitext(os.path.basename(x))[0] 

4198 for x in glob.glob(builtin_model_filename('*'))]) 

4199 

4200 

4201def builtin_model_filename(modelname): 

4202 return util.data_file(os.path.join('earthmodels', modelname+'.nd')) 

4203 

4204 

4205def load_model(fn='ak135-f-continental.m', format='nd', crust2_profile=None): 

4206 ''' 

4207 Load layered earth model from file. 

4208 

4209 :param fn: filename 

4210 :param format: format 

4211 :param crust2_profile: ``(lat, lon)`` or 

4212 :py:class:`pyrocko.dataset.crust2x2.Crust2Profile` object, merge model 

4213 with crustal profile. If ``fn`` is forced to be ``None`` only the 

4214 converted CRUST2.0 profile is returned. 

4215 :returns: object of type :py:class:`LayeredModel` 

4216 

4217 The following formats are currently supported: 

4218 

4219 ============== =========================================================== 

4220 format description 

4221 ============== =========================================================== 

4222 ``'nd'`` 'named discontinuity' format used by the TauP programs 

4223 ``'hyposat'`` format used by the HYPOSAT location program 

4224 ============== =========================================================== 

4225 

4226 The naming of interfaces is translated from the file format's native naming 

4227 to Cake's own convention (See :py:func:`read_nd_model` and 

4228 :py:func:`read_hyposat_model` for details). Cake likes the following 

4229 internal names: ``'conrad'``, ``'moho'``, ``'cmb'`` (core-mantle boundary), 

4230 ``'icb'`` (inner core boundary). 

4231 ''' 

4232 

4233 if fn is not None: 

4234 if format == 'nd': 

4235 if not os.path.exists(fn) and fn in builtin_models(): 

4236 fn = builtin_model_filename(fn) 

4237 reader = read_nd_model(fn) 

4238 elif format == 'hyposat': 

4239 reader = read_hyposat_model(fn) 

4240 else: 

4241 assert False, 'unsupported model format' 

4242 

4243 mod = LayeredModel.from_scanlines(reader) 

4244 if crust2_profile is not None: 

4245 return mod.replaced_crust(crust2_profile) 

4246 

4247 return mod 

4248 

4249 else: 

4250 assert crust2_profile is not None 

4251 profile = anything_to_crust2_profile(crust2_profile) 

4252 return LayeredModel.from_scanlines( 

4253 from_crust2x2_profile(profile)) 

4254 

4255 

4256def castagna_vs_to_vp(vs): 

4257 ''' 

4258 Calculate vp from vs using Castagna's relation. 

4259 

4260 Castagna's relation (the mudrock line) is an empirical relation for vp/vs 

4261 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al., 

4262 1985] 

4263 

4264 vp = 1.16 * vs + 1360 [m/s] 

4265 

4266 :param vs: S-wave velocity [m/s] 

4267 :returns: P-wave velocity [m/s] 

4268 ''' 

4269 

4270 return vs*1.16 + 1360.0 

4271 

4272 

4273def castagna_vp_to_vs(vp): 

4274 ''' 

4275 Calculate vp from vs using Castagna's relation. 

4276 

4277 Castagna's relation (the mudrock line) is an empirical relation for vp/vs 

4278 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al., 

4279 1985] 

4280 

4281 vp = 1.16 * vs + 1360 [m/s] 

4282 

4283 :param vp: P-wave velocity [m/s] 

4284 :returns: S-wave velocity [m/s] 

4285 ''' 

4286 

4287 return (vp - 1360.0) / 1.16 

4288 

4289 

4290def evenize(x, y, minsize=10): 

4291 if x.size < minsize: 

4292 return x 

4293 ry = (y.max()-y.min()) 

4294 if ry == 0: 

4295 return x 

4296 dx = (x[1:] - x[:-1])/(x.max()-x.min()) 

4297 dy = (y[1:] + y[:-1])/ry 

4298 

4299 s = num.zeros(x.size) 

4300 s[1:] = num.cumsum(num.sqrt(dy**2 + dx**2)) 

4301 s2 = num.linspace(0, s[-1], x.size) 

4302 x2 = num.interp(s2, s, x) 

4303 x2[0] = x[0] 

4304 x2[-1] = x[-1] 

4305 return x2 

4306 

4307 

4308def filled(v, *args, **kwargs): 

4309 ''' 

4310 Create NumPy array filled with given value. 

4311 

4312 This works like :py:func:`numpy.ones` but initializes the array with ``v`` 

4313 instead of ones. 

4314 ''' 

4315 x = num.empty(*args, **kwargs) 

4316 x.fill(v) 

4317 return x 

4318 

4319 

4320def next_or_none(i): 

4321 try: 

4322 return next(i) 

4323 except StopIteration: 

4324 return None 

4325 

4326 

4327def reci_or_none(x): 

4328 try: 

4329 return 1./x 

4330 except ZeroDivisionError: 

4331 return None 

4332 

4333 

4334def mult_or_none(a, b): 

4335 if a is None or b is None: 

4336 return None 

4337 return a*b 

4338 

4339 

4340def min_not_none(a, b): 

4341 if a is None: 

4342 return b 

4343 if b is None: 

4344 return a 

4345 return min(a, b) 

4346 

4347 

4348def xytups(xx, yy): 

4349 d = [] 

4350 for x, y in zip(xx, yy): 

4351 if num.isfinite(y): 

4352 d.append((x, y)) 

4353 return d 

4354 

4355 

4356def interp(x, xp, fp, monoton): 

4357 if monoton == 1: 

4358 return xytups( 

4359 x, num.interp(x, xp, fp, left=num.nan, right=num.nan)) 

4360 elif monoton == -1: 

4361 return xytups( 

4362 x, num.interp(x, xp[::-1], fp[::-1], left=num.nan, right=num.nan)) 

4363 else: 

4364 fs = [] 

4365 for xv in x: 

4366 indices = num.where(num.logical_or( 

4367 num.logical_and(xp[:-1] >= xv, xv > xp[1:]), 

4368 num.logical_and(xp[:-1] <= xv, xv < xp[1:])))[0] 

4369 

4370 for i in indices: 

4371 xr = (xv - xp[i])/(xp[i+1]-xp[i]) 

4372 fv = xr*fp[i] + (1.-xr)*fp[i+1] 

4373 fs.append((xv, fv)) 

4374 

4375 return fs 

4376 

4377 

4378def float_or_none(x): 

4379 if x is not None: 

4380 return float(x) 

4381 

4382 

4383def parstore_float(thelocals, obj, *args): 

4384 for k, v in thelocals.items(): 

4385 if k != 'self' and (not args or k in args): 

4386 setattr(obj, k, float_or_none(v))