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 

69P = 1 

70S = 2 

71DOWN = 4 

72UP = -4 

73 

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

75 

76earthradius = config.config().earthradius 

77 

78r2d = 180./math.pi 

79d2r = 1./r2d 

80km = 1000. 

81d2m = d2r*earthradius 

82m2d = 1./d2m 

83sprad2spm = 1.0/(r2d*d2m) 

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

85spm2sprad = 1.0/sprad2spm 

86spkm2sprad = 1.0/sprad2spkm 

87 

88 

89class CakeError(Exception): 

90 pass 

91 

92 

93class InvalidArguments(CakeError): 

94 pass 

95 

96 

97class Material(object): 

98 ''' 

99 Isotropic elastic material. 

100 

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

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

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

104 :param qp: P-wave attenuation Qp 

105 :param qs: S-wave attenuation Qs 

106 :param poisson: Poisson ratio 

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

108 :param qk: bulk attenuation Qk 

109 :param qmu: shear attenuation Qmu 

110 

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

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

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

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

115 :type burgers: tuple 

116 

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

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

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

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

121 0 and alpha=1. 

122 

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

124 

125 The main material properties are considered independant and are accessible 

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

127 

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

129 

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

131 instance methods. 

132 ''' 

133 

134 def __init__( 

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

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

137 

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

139 

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

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

142 raise InvalidArguments( 

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

144 'should not be given.') 

145 

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

147 self.vp = 5800. 

148 if poisson is None: 

149 poisson = 0.25 

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

151 

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

153 if poisson is not None: 

154 raise InvalidArguments( 

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

156 'are given.') 

157 

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

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

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

161 

162 elif vp is not None and vs is None: 

163 if poisson is None: 

164 poisson = 0.25 

165 

166 if lame is not None: 

167 raise InvalidArguments( 

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

169 

170 poisson = float(poisson) 

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

172 

173 elif vp is None and vs is not None: 

174 if poisson is None: 

175 poisson = 0.25 

176 if lame is not None: 

177 raise InvalidArguments( 

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

179 

180 poisson = float(poisson) 

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

182 

183 else: 

184 raise InvalidArguments( 

185 'Invalid combination of input parameters in material ' 

186 'definition.') 

187 

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

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

190 raise InvalidArguments( 

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

192 

193 if qp is None: 

194 if self.vs != 0.0: 

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

196 self.qp = self.qs / s 

197 else: 

198 self.qp = 1456. 

199 

200 if qs is None: 

201 if self.vs != 0.0: 

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

203 self.qs = self.qp * s 

204 else: 

205 self.vs = 600. 

206 

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

208 if self.vs == 0.: 

209 self.qs = 0. 

210 self.qp = 5782e4 

211 else: 

212 self.qs = 600. 

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

214 self.qp = self.qs/s 

215 

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

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

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

219 self.qp = qk 

220 else: 

221 if num.isinf(qk): 

222 self.qp = qmu/s 

223 else: 

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

225 self.qs = qmu 

226 else: 

227 raise InvalidArguments( 

228 'Invalid combination of input parameters in material ' 

229 'definition.') 

230 

231 if burgers is None: 

232 burgers = DEFAULT_BURGERS 

233 

234 self.burger_eta1 = burgers[0] 

235 self.burger_eta2 = burgers[1] 

236 self.burger_valpha = burgers[2] 

237 

238 def astuple(self): 

239 ''' 

240 Get independant material properties as a tuple. 

241 

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

243 ''' 

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

245 

246 def __eq__(self, other): 

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

248 

249 def lame(self): 

250 ''' 

251 Get Lame's parameter lambda and shear modulus. 

252 ''' 

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

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

255 return lam, mu 

256 

257 def lame_lambda(self): 

258 ''' 

259 Get Lame's parameter lambda. 

260 

261 Returned units are [Pa]. 

262 ''' 

263 lam, _ = self.lame() 

264 return lam 

265 

266 def shear_modulus(self): 

267 ''' 

268 Get shear modulus. 

269 

270 Returned units are [Pa]. 

271 ''' 

272 return self.vs**2 * self.rho 

273 

274 def poisson(self): 

275 ''' 

276 Get Poisson's ratio. 

277 ''' 

278 lam, mu = self.lame() 

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

280 

281 def bulk(self): 

282 ''' 

283 Get bulk modulus. 

284 ''' 

285 lam, mu = self.lame() 

286 return lam + 2.0*mu/3.0 

287 

288 def youngs(self): 

289 ''' 

290 Get Young's modulus. 

291 ''' 

292 lam, mu = self.lame() 

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

294 

295 def vp_vs_ratio(self): 

296 ''' 

297 Get vp/vs ratio. 

298 ''' 

299 return self.vp/self.vs 

300 

301 def qmu(self): 

302 ''' 

303 Get shear attenuation coefficient Qmu. 

304 ''' 

305 return self.qs 

306 

307 def qk(self): 

308 ''' 

309 Get bulk attenuation coefficient Qk. 

310 ''' 

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

312 return self.qp 

313 else: 

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

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

316 if denom <= 0.0: 

317 return num.inf 

318 else: 

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

320 

321 def burgers(self): 

322 ''' 

323 Get Burger parameters. 

324 ''' 

325 return self.burger_eta1, self.burger_eta2, self.burger_valpha 

326 

327 def _rayleigh_equation(self, cr): 

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

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

330 if cr_a > 1.0 or cr_b > 1.0: 

331 return None 

332 

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

334 

335 def rayleigh(self): 

336 ''' 

337 Get Rayleigh velocity assuming a homogenous halfspace. 

338 

339 Returned units are [m/s]. 

340 ''' 

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

342 

343 def _has_default_burgers(self): 

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

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

346 self.burger_valpha == DEFAULT_BURGERS[2]: 

347 return True 

348 return False 

349 

350 def describe(self): 

351 ''' 

352 Get a readable listing of the material properties. 

353 ''' 

354 template = ''' 

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

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

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

358Lame lambda [GPa] : %12g 

359Lame shear modulus [GPa] : %12g 

360Poisson ratio : %12g 

361Bulk modulus [GPa] : %12g 

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

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

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

365Qp P-wave attenuation : %12g 

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

367Qk bulk attenuation : %12g 

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

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

370relaxation: valpha : %12g 

371'''.strip() 

372 

373 return template % ( 

374 self.vp/km, 

375 self.vs/km, 

376 self.vp/self.vs, 

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

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

379 self.poisson(), 

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

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

382 self.rayleigh()/km, 

383 self.rho/km, 

384 self.qp, 

385 self.qs, 

386 self.qk(), 

387 self.burger_eta1*1e-9, 

388 self.burger_eta2*1e-9, 

389 self.burger_valpha) 

390 

391 def __str__(self): 

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

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

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

395 

396 def __repr__(self): 

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

398 tuple(repr(x) for x in ( 

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

400 

401 

402class Leg(object): 

403 ''' 

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

405 

406 **Attributes:** 

407 

408 To be considered as read-only. 

409 

410 .. py:attribute:: departure 

411 

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

413 upward or downward departure. 

414 

415 .. py:attribute:: mode 

416 

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

418 propagation mode. 

419 

420 .. py:attribute:: depthmin 

421 

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

423 minimum depth. 

424 

425 .. py:attribute:: depthmax 

426 

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

428 maximum depth. 

429 

430 ''' 

431 

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

433 self.departure = departure 

434 self.mode = mode 

435 self.depthmin = None 

436 self.depthmax = None 

437 

438 def set_depthmin(self, depthmin): 

439 self.depthmin = depthmin 

440 

441 def set_depthmax(self, depthmax): 

442 self.depthmax = depthmax 

443 

444 def __str__(self): 

445 def sd(d): 

446 if isinstance(d, float): 

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

448 else: 

449 return 'interface %s' % d 

450 

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

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

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

454 

455 sc = [] 

456 if self.depthmax is not None: 

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

458 if self.depthmin is not None: 

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

460 

461 if sc: 

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

463 

464 return s 

465 

466 

467class InvalidKneeDef(CakeError): 

468 pass 

469 

470 

471class Knee(object): 

472 ''' 

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

474 

475 **Attributes:** 

476 

477 To be considered as read-only. 

478 

479 .. py:attribute:: depth 

480 

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

482 a string or a number. 

483 

484 .. py:attribute:: direction 

485 

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

487 the incoming direction. 

488 

489 .. py:attribute:: in_mode 

490 

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

492 type of mode of the incoming wave. 

493 

494 .. py:attribute:: out_mode 

495 

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

497 type of mode of the outgoing wave. 

498 

499 .. py:attribute:: conversion 

500 

501 Boolean, whether there is a mode conversion involved. 

502 

503 .. py:attribute:: reflection 

504 

505 Boolean, whether there is a reflection involved. 

506 

507 .. py:attribute:: headwave 

508 

509 Boolean, whether there is headwave propagation involved. 

510 

511 ''' 

512 

513 defaults = dict( 

514 depth='surface', 

515 direction=UP, 

516 conversion=True, 

517 reflection=False, 

518 headwave=False, 

519 in_setup_state=True) 

520 

521 defaults_surface = dict( 

522 depth='surface', 

523 direction=UP, 

524 conversion=False, 

525 reflection=True, 

526 headwave=False, 

527 in_setup_state=True) 

528 

529 def __init__(self, *args): 

530 if args: 

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

532 self.out_mode) = args 

533 

534 self.conversion = self.in_mode != self.out_mode 

535 self.in_setup_state = False 

536 

537 def default(self, k): 

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

539 if depth == 'surface': 

540 return Knee.defaults_surface[k] 

541 else: 

542 return Knee.defaults[k] 

543 

544 def __setattr__(self, k, v): 

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

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

547 else: 

548 self.__dict__[k] = v 

549 

550 def __getattr__(self, k): 

551 if k.startswith('__'): 

552 raise AttributeError(k) 

553 

554 if k not in self.__dict__: 

555 return self.default(k) 

556 

557 def set_modes(self, in_leg, out_leg): 

558 

559 if out_leg.departure == UP and ( 

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

561 

562 raise InvalidKneeDef( 

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

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

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

566 

567 if out_leg.departure == DOWN and ( 

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

569 

570 raise InvalidKneeDef( 

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

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

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

574 

575 self.in_mode = in_leg.mode 

576 self.out_mode = out_leg.mode 

577 

578 def at_surface(self): 

579 return self.depth == 'surface' 

580 

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

582 ''' 

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

584 propagation mode, and direction. 

585 ''' 

586 

587 if isinstance(self.depth, float): 

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

589 return False 

590 else: 

591 if discontinuity.name != self.depth: 

592 return False 

593 

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

595 

596 def out_direction(self): 

597 ''' 

598 Get outgoing direction. 

599 

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

601 ''' 

602 

603 if self.reflection: 

604 return - self.direction 

605 else: 

606 return self.direction 

607 

608 def __str__(self): 

609 x = [] 

610 if self.reflection: 

611 if self.at_surface(): 

612 x.append('surface') 

613 else: 

614 if not self.headwave: 

615 if self.direction == UP: 

616 x.append('underside') 

617 else: 

618 x.append('upperside') 

619 

620 if self.headwave: 

621 x.append('headwave propagation along') 

622 elif self.reflection and self.conversion: 

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

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

625 if not self.at_surface(): 

626 x.append('at') 

627 elif self.reflection: 

628 x.append('reflection') 

629 if not self.at_surface(): 

630 x.append('at') 

631 elif self.conversion: 

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

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

634 else: 

635 x.append('passing through') 

636 

637 if isinstance(self.depth, float): 

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

639 else: 

640 if not self.at_surface(): 

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

642 

643 if not self.reflection: 

644 if self.direction == UP: 

645 x.append('on upgoing path') 

646 else: 

647 x.append('on downgoing path') 

648 

649 return ' '.join(x) 

650 

651 

652class Head(Knee): 

653 def __init__(self, *args): 

654 if args: 

655 z, in_direction, mode = args 

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

657 else: 

658 Knee.__init__(self) 

659 

660 def __str__(self): 

661 x = ['propagation as headwave'] 

662 if isinstance(self.depth, float): 

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

664 else: 

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

666 

667 return ' '.join(x) 

668 

669 

670class UnknownClassicPhase(CakeError): 

671 def __init__(self, phasename): 

672 self.phasename = phasename 

673 

674 def __str__(self): 

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

676 

677 

678class PhaseDefParseError(CakeError): 

679 ''' 

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

681 definition string. 

682 ''' 

683 

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

685 self.definition = definition 

686 self.position = position 

687 self.exception = exception 

688 

689 def __str__(self): 

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

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

692 

693 

694class PhaseDef(object): 

695 

696 ''' 

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

698 

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

700 syntax 

701 

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

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

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

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

706 are not completely compatible with those. 

707 

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

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

710 

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

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

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

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

715 take-off in upward direction. 

716 

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

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

719 

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

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

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

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

724 ``v_DEPTH`` 

725 

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

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

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

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

730 

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

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

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

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

735 leg character, respectively. 

736 

737 **Examples:** 

738 

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

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

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

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

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

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

745 interface closest to that) 

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

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

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

749 discontinuity 

750 

751 **Usage:** 

752 

753 >>> from pyrocko.cake import PhaseDef 

754 # must escape the backslash 

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

756 >>> print my_crazy_phase 

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

758 - P mode propagation, departing upward 

759 - surface reflection 

760 - P mode propagation, departing downward 

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

762 - S mode propagation, departing upward 

763 - surface reflection with conversion from S to P 

764 - P mode propagation, departing downward 

765 - arriving at target from above 

766 

767 .. note:: 

768 

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

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

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

772 selection of conversion and reflection coefficients, which 

773 currently only deal with the P-SV case. 

774 ''' 

775 

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

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

778 

779 @staticmethod 

780 def classic_definitions(): 

781 defs = {} 

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

783 for r in 'mc': 

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

785 defs[a+r+b] = [ 

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

787 

788 # Pg, P, S, Sg 

789 for a in 'PS': 

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

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

792 a, a.lower())] 

793 

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

795 

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

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

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

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

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

801 

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

803 for a in 'PS': 

804 for b in 'PS': 

805 for c in 'PS': 

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

807 

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

809 

810 # Pc, Pdiff, Sc, ... 

811 for x in 'PS': 

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

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

814 

815 # depth phases 

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

817 if k not in 'ps': 

818 for x in 'ps': 

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

820 

821 return defs 

822 

823 @staticmethod 

824 def classic(phasename): 

825 ''' 

826 Get phase definitions based on classic phase name. 

827 

828 :param phasename: classic name of a phase 

829 :returns: list of PhaseDef objects 

830 

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

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

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

834 ''' 

835 

836 defs = PhaseDef.classic_definitions() 

837 if phasename not in defs: 

838 raise UnknownClassicPhase(phasename) 

839 

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

841 

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

843 

844 state = 0 

845 sdepth = '' 

846 sinterface = '' 

847 depthmax = depthmin = None 

848 depthlim = None 

849 depthlimtype = None 

850 sdepthlim = '' 

851 events = [] 

852 direction_stop = UP 

853 need_leg = True 

854 ic = 0 

855 if definition is not None: 

856 knee = Knee() 

857 try: 

858 for ic, c in enumerate(definition): 

859 

860 if state in (0, 1): 

861 

862 if c in '0123456789.': 

863 need_leg = True 

864 state = 1 

865 sdepth += c 

866 continue 

867 

868 elif state == 1: 

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

870 state = 0 

871 

872 if state == 2: 

873 if c == ')': 

874 knee.depth = sinterface 

875 state = 0 

876 else: 

877 sinterface += c 

878 

879 continue 

880 

881 if state in (3, 4): 

882 

883 if state == 3: 

884 if c in '0123456789.': 

885 sdepthlim += c 

886 continue 

887 elif c == '(': 

888 state = 4 

889 continue 

890 else: 

891 depthlim = float(sdepthlim)*1000. 

892 if depthlimtype == '<': 

893 depthmax = depthlim 

894 else: 

895 depthmin = depthlim 

896 state = 0 

897 

898 elif state == 4: 

899 if c == ')': 

900 depthlim = sdepthlim 

901 if depthlimtype == '<': 

902 depthmax = depthlim 

903 else: 

904 depthmin = depthlim 

905 state = 0 

906 continue 

907 else: 

908 sdepthlim += c 

909 continue 

910 

911 if state == 0: 

912 

913 if c == '(': 

914 need_leg = True 

915 state = 2 

916 continue 

917 

918 elif c in '<>': 

919 state = 3 

920 depthlim = None 

921 sdepthlim = '' 

922 depthlimtype = c 

923 continue 

924 

925 elif c in 'psPS': 

926 leg = Leg() 

927 if c in 'ps': 

928 leg.departure = UP 

929 else: 

930 leg.departure = DOWN 

931 leg.mode = imode(c) 

932 

933 if events: 

934 in_leg = events[-1] 

935 if depthmin is not None: 

936 in_leg.set_depthmin(depthmin) 

937 depthmin = None 

938 if depthmax is not None: 

939 in_leg.set_depthmax(depthmax) 

940 depthmax = None 

941 

942 if in_leg.mode != leg.mode: 

943 knee.conversion = True 

944 else: 

945 knee.conversion = False 

946 

947 if not knee.reflection: 

948 if c in 'ps': 

949 knee.direction = UP 

950 else: 

951 knee.direction = DOWN 

952 

953 knee.set_modes(in_leg, leg) 

954 knee.in_setup_state = False 

955 events.append(knee) 

956 knee = Knee() 

957 sdepth = '' 

958 sinterface = '' 

959 

960 events.append(leg) 

961 need_leg = False 

962 continue 

963 

964 elif c == '^': 

965 need_leg = True 

966 knee.direction = UP 

967 knee.reflection = True 

968 continue 

969 

970 elif c == 'v': 

971 need_leg = True 

972 knee.direction = DOWN 

973 knee.reflection = True 

974 continue 

975 

976 elif c == '_': 

977 need_leg = True 

978 knee.headwave = True 

979 continue 

980 

981 elif c == '\\': 

982 direction_stop = DOWN 

983 continue 

984 

985 else: 

986 raise PhaseDefParseError( 

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

988 

989 if state == 3: 

990 depthlim = float(sdepthlim)*1000. 

991 if depthlimtype == '<': 

992 depthmax = depthlim 

993 else: 

994 depthmin = depthlim 

995 state = 0 

996 

997 except (ValueError, InvalidKneeDef) as e: 

998 raise PhaseDefParseError(definition, ic, e) 

999 

1000 if state != 0 or need_leg: 

1001 raise PhaseDefParseError( 

1002 definition, ic, 'unfinished expression') 

1003 

1004 if events and depthmin is not None: 

1005 events[-1].set_depthmin(depthmin) 

1006 if events and depthmax is not None: 

1007 events[-1].set_depthmax(depthmax) 

1008 

1009 self._definition = definition 

1010 self._classicname = classicname 

1011 self._events = events 

1012 self._direction_stop = direction_stop 

1013 

1014 def __iter__(self): 

1015 for ev in self._events: 

1016 yield ev 

1017 

1018 def append(self, ev): 

1019 self._events.append(ev) 

1020 

1021 def first_leg(self): 

1022 ''' 

1023 Get the first leg in phase definition. 

1024 ''' 

1025 return self._events[0] 

1026 

1027 def last_leg(self): 

1028 ''' 

1029 Get the last leg in phase definition. 

1030 ''' 

1031 return self._events[-1] 

1032 

1033 def legs(self): 

1034 ''' 

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

1036 within this phase definition. 

1037 ''' 

1038 

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

1040 

1041 def knees(self): 

1042 ''' 

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

1044 phase definition. 

1045 ''' 

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

1047 

1048 def definition(self): 

1049 ''' 

1050 Get original definition of the phase. 

1051 ''' 

1052 return self._definition 

1053 

1054 def given_name(self): 

1055 ''' 

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

1057 ''' 

1058 

1059 if self._classicname: 

1060 return self._classicname 

1061 else: 

1062 return self._definition 

1063 

1064 def direction_start(self): 

1065 return self.first_leg().departure 

1066 

1067 def direction_stop(self): 

1068 return self._direction_stop 

1069 

1070 def headwave_knee(self): 

1071 for el in self: 

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

1073 return el 

1074 return None 

1075 

1076 def used_repr(self): 

1077 ''' 

1078 Translate into textual representation (cake phase syntax). 

1079 ''' 

1080 def strdepth(x): 

1081 if isinstance(x, float): 

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

1083 else: 

1084 return '(%s)' % x 

1085 

1086 x = [] 

1087 for el in self: 

1088 if type(el) == Leg: 

1089 if el.departure == UP: 

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

1091 else: 

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

1093 

1094 if el.depthmax is not None: 

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

1096 

1097 if el.depthmin is not None: 

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

1099 

1100 elif type(el) == Knee: 

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

1102 if el.direction == DOWN: 

1103 x.append('v') 

1104 else: 

1105 x.append('^') 

1106 if el.headwave: 

1107 x.append('_') 

1108 if not el.at_surface(): 

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

1110 

1111 elif type(el) == Head: 

1112 x.append('_') 

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

1114 

1115 if self._direction_stop == DOWN: 

1116 x.append('\\') 

1117 

1118 return ''.join(x) 

1119 

1120 def __repr__(self): 

1121 if self._definition is not None: 

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

1123 else: 

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

1125 

1126 def __str__(self): 

1127 orig = '' 

1128 used = self.used_repr() 

1129 if self._definition != used: 

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

1131 

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

1133 self._direction_stop == DOWN] 

1134 

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

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

1137 

1138 def copy(self): 

1139 ''' 

1140 Get a deep copy of it. 

1141 ''' 

1142 return copy.deepcopy(self) 

1143 

1144 

1145def to_phase_defs(phases): 

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

1147 phases = [phases] 

1148 

1149 phases_out = [] 

1150 for phase in phases: 

1151 if isinstance(phase, str): 

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

1153 elif isinstance(phase, PhaseDef): 

1154 phases_out.append(phase) 

1155 else: 

1156 raise PhaseDefParseError('invalid phase definition') 

1157 

1158 return phases_out 

1159 

1160 

1161def csswap(x): 

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

1163 

1164 

1165def psv_surface_ind(in_mode, out_mode): 

1166 ''' 

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

1168 surface. 

1169 ''' 

1170 

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

1172 

1173 

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

1175 ''' 

1176 Scatter matrix for free surface reflection/conversions. 

1177 

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

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

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

1181 returned 

1182 :returns: Scatter matrix 

1183 

1184 The scatter matrix is ordered as follows:: 

1185 

1186 [[PP, PS], 

1187 [SP, SS]] 

1188 

1189 The formulas given in Aki & Richards are used. 

1190 ''' 

1191 

1192 vp, vs, rho = material.vp, material.vs, material.rho 

1193 sinphi = p * vp 

1194 sinlam = p * vs 

1195 cosphi = csswap(sinphi) 

1196 coslam = csswap(sinlam) 

1197 

1198 if vs == 0.0: 

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

1200 

1201 else: 

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

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

1204 denom = vsp_term**2 + pcc_term 

1205 

1206 scatter = num.array([ 

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

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

1209 dtype=complex) / denom 

1210 

1211 if not energy: 

1212 return scatter 

1213 else: 

1214 eps = 1e-16 

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

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

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

1218 return num.real(escatter) 

1219 

1220 

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

1222 ''' 

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

1224 solid-solid interface. 

1225 ''' 

1226 

1227 return ( 

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

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

1230 

1231 

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

1233 ''' 

1234 Scatter matrix for solid-solid interface. 

1235 

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

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

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

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

1240 returned 

1241 :returns: Scatter matrix 

1242 

1243 The scatter matrix is ordered as follows:: 

1244 

1245 [[P1P1, S1P1, P2P1, S2P1], 

1246 [P1S1, S1S1, P2S1, S2S1], 

1247 [P1P2, S1P2, P2P2, S2P2], 

1248 [P1S2, S1S2, P2S2, S2S2]] 

1249 

1250 The formulas given in Aki & Richards are used. 

1251 ''' 

1252 

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

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

1255 

1256 sinphi1 = p * vp1 

1257 cosphi1 = csswap(sinphi1) 

1258 sinlam1 = p * vs1 

1259 coslam1 = csswap(sinlam1) 

1260 sinphi2 = p * vp2 

1261 cosphi2 = csswap(sinphi2) 

1262 sinlam2 = p * vs2 

1263 coslam2 = csswap(sinlam2) 

1264 

1265 # from aki and richards 

1266 M = num.array([ 

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

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

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

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

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

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

1273 dtype=complex) 

1274 N = M.copy() 

1275 N[0] *= -1.0 

1276 N[3] *= -1.0 

1277 

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

1279 

1280 if not energy: 

1281 return scatter 

1282 else: 

1283 eps = 1e-16 

1284 if vs1 == 0.: 

1285 vs1 = vp1*1e-16 

1286 if vs2 == 0.: 

1287 vs2 = vp2*1e-16 

1288 normvec = num.array([ 

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

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

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

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

1293 

1294 return num.real(escatter) 

1295 

1296 

1297class BadPotIntCoefs(CakeError): 

1298 pass 

1299 

1300 

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

1302 eps = r2*1e-9 

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

1304 c1c2 = 1. 

1305 else: 

1306 c1c2 = c1/c2 

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

1308 if abs(b) > 10.: 

1309 raise BadPotIntCoefs() 

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

1311 return a, b 

1312 

1313 

1314def imode(s): 

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

1316 return P 

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

1318 return S 

1319 

1320 

1321def smode(i): 

1322 if i == P: 

1323 return 'p' 

1324 elif i == S: 

1325 return 's' 

1326 

1327 

1328class PathFailed(CakeError): 

1329 pass 

1330 

1331 

1332class SurfaceReached(PathFailed): 

1333 pass 

1334 

1335 

1336class BottomReached(PathFailed): 

1337 pass 

1338 

1339 

1340class MaxDepthReached(PathFailed): 

1341 pass 

1342 

1343 

1344class MinDepthReached(PathFailed): 

1345 pass 

1346 

1347 

1348class Trapped(PathFailed): 

1349 pass 

1350 

1351 

1352class NotPhaseConform(PathFailed): 

1353 pass 

1354 

1355 

1356class CannotPropagate(PathFailed): 

1357 def __init__(self, direction, ilayer): 

1358 PathFailed.__init__(self) 

1359 self._direction = direction 

1360 self._ilayer = ilayer 

1361 

1362 def __str__(self): 

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

1364 self._ilayer, { 

1365 UP: 'below', 

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

1367 

1368 

1369class Layer(object): 

1370 ''' 

1371 Representation of a layer in a layered earth model. 

1372 

1373 :param ztop: depth of top of layer 

1374 :param zbot: depth of bottom of layer 

1375 :param name: name of layer (optional) 

1376 

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

1378 ''' 

1379 

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

1381 self.ztop = ztop 

1382 self.zbot = zbot 

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

1384 self.name = name 

1385 self.ilayer = None 

1386 

1387 def _update_potint_coefs(self): 

1388 potint_p = potint_s = False 

1389 try: 

1390 self._ppic = potint_coefs( 

1391 self.mbot.vp, self.mtop.vp, 

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

1393 potint_p = True 

1394 except BadPotIntCoefs: 

1395 pass 

1396 

1397 potint_s = False 

1398 try: 

1399 self._spic = potint_coefs( 

1400 self.mbot.vs, self.mtop.vs, 

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

1402 potint_s = True 

1403 except BadPotIntCoefs: 

1404 pass 

1405 

1406 assert P == 1 and S == 2 

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

1408 

1409 def potint_coefs(self, mode): 

1410 ''' 

1411 Get coefficients for potential interpolation. 

1412 

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

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

1415 ''' 

1416 

1417 if mode == P: 

1418 return self._ppic 

1419 else: 

1420 return self._spic 

1421 

1422 def contains(self, z): 

1423 ''' 

1424 Tolerantly check if a given depth is within the layer 

1425 (including boundaries). 

1426 ''' 

1427 

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

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

1430 

1431 def inner(self, z): 

1432 ''' 

1433 Tolerantly check if a given depth is within the layer 

1434 (not including boundaries). 

1435 ''' 

1436 

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

1438 self.at_bottom(z) and not \ 

1439 self.at_top(z) 

1440 

1441 def at_bottom(self, z): 

1442 ''' 

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

1444 ''' 

1445 

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

1447 

1448 def at_top(self, z): 

1449 ''' 

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

1451 ''' 

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

1453 

1454 def pflat_top(self, p): 

1455 ''' 

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

1457 layer. 

1458 ''' 

1459 return p / (earthradius-self.ztop) 

1460 

1461 def pflat_bottom(self, p): 

1462 ''' 

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

1464 of layer. 

1465 ''' 

1466 return p / (earthradius-self.zbot) 

1467 

1468 def pflat(self, p, z): 

1469 ''' 

1470 Convert spherical ray parameter to local flat ray parameter for 

1471 given depth. 

1472 ''' 

1473 return p / (earthradius-z) 

1474 

1475 def v_potint(self, mode, z): 

1476 a, b = self.potint_coefs(mode) 

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

1478 

1479 def u_potint(self, mode, z): 

1480 a, b = self.potint_coefs(mode) 

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

1482 

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

1484 ''' 

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

1486 parameter. 

1487 

1488 :param p: ray parameter (spherical) 

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

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

1491 to a part of the layer 

1492 

1493 This implementation uses analytic formulas valid for a spherical earth 

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

1495 interpolation of the form 

1496 

1497 c(z) = a*z^b 

1498 ''' 

1499 utop, ubot = self.u_top_bottom(mode) 

1500 a, b = self.potint_coefs(mode) 

1501 ztop = self.ztop 

1502 zbot = self.zbot 

1503 if zpart is not None: 

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

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

1506 ztop, zbot = zpart 

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

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

1509 

1510 r1 = radius(zbot) 

1511 r2 = radius(ztop) 

1512 burger_eta1 = r1 * ubot 

1513 burger_eta2 = r2 * utop 

1514 if b != 1: 

1515 def cpe(eta): 

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

1517 

1518 def sep(eta): 

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

1520 

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

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

1523 else: 

1524 lr = math.log(r2/r1) 

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

1526 x = p/sap * lr 

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

1528 

1529 x *= r2d 

1530 

1531 return x, t 

1532 

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

1534 ''' 

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

1536 within the layer. 

1537 ''' 

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

1539 

1540 def tests(self, p, mode): 

1541 utop, ubot = self.u_top_bottom(mode) 

1542 return ( 

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

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

1545 

1546 def zturn_potint(self, p, mode): 

1547 ''' 

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

1549 ''' 

1550 

1551 a, b = self.potint_coefs(mode) 

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

1553 return earthradius-r 

1554 

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

1556 ''' 

1557 Propagate ray through layer. 

1558 

1559 :param p: ray parameter 

1560 :param mode: propagation mode 

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

1562 ''' 

1563 if direction == DOWN: 

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

1565 else: 

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

1567 

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

1569 raise CannotPropagate(direction, self.ilayer) 

1570 

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

1572 return -direction 

1573 else: 

1574 return direction 

1575 

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

1577 ''' 

1578 Change layer thinkness and interpolate material if required. 

1579 ''' 

1580 if depth_min: 

1581 mtop = self.material(depth_min) 

1582 

1583 if depth_max: 

1584 mbot = self.material(depth_max) 

1585 

1586 self.mtop = mtop if depth_min else self.mtop 

1587 self.mbot = mbot if depth_max else self.mbot 

1588 self.ztop = depth_min if depth_min else self.ztop 

1589 self.zbot = depth_max if depth_max else self.zbot 

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

1591 

1592 

1593class DoesNotTurn(CakeError): 

1594 pass 

1595 

1596 

1597def radius(z): 

1598 return earthradius - z 

1599 

1600 

1601class HomogeneousLayer(Layer): 

1602 ''' 

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

1604 

1605 Base class: :py:class:`Layer`. 

1606 ''' 

1607 

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

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

1610 self.m = m 

1611 self.mtop = m 

1612 self.mbot = m 

1613 self._update_potint_coefs() 

1614 

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

1616 if ztop is None: 

1617 ztop = self.ztop 

1618 

1619 if zbot is None: 

1620 zbot = self.zbot 

1621 

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

1623 

1624 def material(self, z): 

1625 return self.m 

1626 

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

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

1629 return self.u_potint(mode, z) 

1630 

1631 if mode == P: 

1632 return 1./self.m.vp 

1633 if mode == S: 

1634 return 1./self.m.vs 

1635 

1636 def u_top_bottom(self, mode): 

1637 u = self.u(mode) 

1638 return u, u 

1639 

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

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

1642 return self.v_potint(mode, z) 

1643 

1644 if mode == P: 

1645 v = self.m.vp 

1646 if mode == S: 

1647 v = self.m.vs 

1648 

1649 if num.isscalar(z): 

1650 return v 

1651 else: 

1652 return filled(v, len(z)) 

1653 

1654 def v_top_bottom(self, mode): 

1655 v = self.v(mode) 

1656 return v, v 

1657 

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

1659 if self._use_potential_interpolation[mode]: 

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

1661 

1662 u = self.u(mode) 

1663 pflat = self.pflat_bottom(p) 

1664 if zpart is None: 

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

1666 else: 

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

1668 

1669 u = self.u(mode) 

1670 eps = u*0.001 

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

1672 

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

1674 t = u**2 * dz / denom 

1675 return x, t 

1676 

1677 def zturn(self, p, mode): 

1678 if self._use_potential_interpolation[mode]: 

1679 return self.zturn_potint(p, mode) 

1680 

1681 raise DoesNotTurn() 

1682 

1683 def split(self, z): 

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

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

1686 upper.ilayer = self.ilayer 

1687 lower.ilayer = self.ilayer 

1688 return upper, lower 

1689 

1690 def __str__(self): 

1691 if self.name: 

1692 name = self.name + ' ' 

1693 else: 

1694 name = '' 

1695 

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

1697 for mode in (P, S)) 

1698 

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

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

1701 

1702 

1703class GradientLayer(Layer): 

1704 ''' 

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

1706 

1707 Base class: :py:class:`Layer`. 

1708 ''' 

1709 

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

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

1712 self.mtop = mtop 

1713 self.mbot = mbot 

1714 self._update_potint_coefs() 

1715 

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

1717 if ztop is None: 

1718 ztop = self.ztop 

1719 

1720 if zbot is None: 

1721 zbot = self.zbot 

1722 

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

1724 

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

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

1727 

1728 def material(self, z): 

1729 dtop = self.mtop.astuple() 

1730 dbot = self.mbot.astuple() 

1731 d = [ 

1732 self.interpolate(z, ptop, pbot) 

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

1734 

1735 return Material(*d) 

1736 

1737 def u_top_bottom(self, mode): 

1738 if mode == P: 

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

1740 if mode == S: 

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

1742 

1743 def u(self, mode, z): 

1744 if self._use_potential_interpolation[mode]: 

1745 return self.u_potint(mode, z) 

1746 

1747 if mode == P: 

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

1749 if mode == S: 

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

1751 

1752 def v_top_bottom(self, mode): 

1753 if mode == P: 

1754 return self.mtop.vp, self.mbot.vp 

1755 if mode == S: 

1756 return self.mtop.vs, self.mbot.vs 

1757 

1758 def v(self, mode, z): 

1759 if self._use_potential_interpolation[mode]: 

1760 return self.v_potint(mode, z) 

1761 

1762 if mode == P: 

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

1764 if mode == S: 

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

1766 

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

1768 if self._use_potential_interpolation[mode]: 

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

1770 

1771 utop, ubot = self.u_top_bottom(mode) 

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

1773 

1774 pflat = self.pflat_bottom(p) 

1775 if zpart is not None: 

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

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

1778 

1779 peps = 1e-16 

1780 pdp = pflat + peps 

1781 

1782 def func(u): 

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

1784 xx = eta/u 

1785 tt = num.where( 

1786 pflat <= u, 

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

1788 0.0) 

1789 

1790 return xx, tt 

1791 

1792 xxtop, tttop = func(utop) 

1793 xxbot, ttbot = func(ubot) 

1794 

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

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

1797 

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

1799 return x, t 

1800 

1801 def zturn(self, p, mode): 

1802 if self._use_potential_interpolation[mode]: 

1803 return self.zturn_potint(p, mode) 

1804 pflat = self.pflat_bottom(p) 

1805 vtop, vbot = self.v_top_bottom(mode) 

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

1807 (vbot-vtop) + self.ztop 

1808 

1809 def split(self, z): 

1810 mmid = self.material(z) 

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

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

1813 upper.ilayer = self.ilayer 

1814 lower.ilayer = self.ilayer 

1815 return upper, lower 

1816 

1817 def __str__(self): 

1818 if self.name: 

1819 name = self.name + ' ' 

1820 else: 

1821 name = '' 

1822 

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

1824 for mode in (P, S)) 

1825 

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

1827 %s 

1828 %s''' % ( 

1829 self.ilayer, 

1830 name, 

1831 self.ztop/km, 

1832 self.zbot/km, 

1833 calcmode, 

1834 self.mtop, 

1835 self.mbot) 

1836 

1837 

1838class Discontinuity(object): 

1839 ''' 

1840 Base class for discontinuities in layered earth model. 

1841 

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

1843 ''' 

1844 

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

1846 self.z = z 

1847 self.zbot = z 

1848 self.ztop = z 

1849 self.name = name 

1850 

1851 def change_depth(self, z): 

1852 self.z = z 

1853 self.zbot = z 

1854 self.ztop = z 

1855 

1856 def copy(self): 

1857 return copy.deepcopy(self) 

1858 

1859 

1860class Interface(Discontinuity): 

1861 ''' 

1862 Representation of an interface in a layered earth model. 

1863 

1864 Base class: :py:class:`Discontinuity`. 

1865 ''' 

1866 

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

1868 Discontinuity.__init__(self, z, name) 

1869 self.mabove = mabove 

1870 self.mbelow = mbelow 

1871 

1872 def __str__(self): 

1873 if self.name is None: 

1874 return 'interface' 

1875 else: 

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

1877 

1878 def u_top_bottom(self, mode): 

1879 if mode == P: 

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

1881 if mode == S: 

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

1883 

1884 def critical_ps(self, mode): 

1885 uabove, ubelow = self.u_top_bottom(mode) 

1886 return ( 

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

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

1889 

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

1891 uabove, ubelow = self.u_top_bottom(mode) 

1892 if direction == DOWN: 

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

1894 return direction 

1895 else: 

1896 return -direction 

1897 if direction == UP: 

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

1899 return direction 

1900 else: 

1901 return -direction 

1902 

1903 def pflat(self, p): 

1904 return p / (earthradius-self.z) 

1905 

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

1907 scatter = psv_solid( 

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

1909 return scatter[ 

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

1911 

1912 

1913class Surface(Discontinuity): 

1914 ''' 

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

1916 

1917 Base class: :py:class:`Discontinuity`. 

1918 ''' 

1919 

1920 def __init__(self, z, mbelow): 

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

1922 self.z = z 

1923 self.mbelow = mbelow 

1924 

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

1926 return direction # no implicit reflection at surface 

1927 

1928 def u_top_bottom(self, mode): 

1929 if mode == P: 

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

1931 if mode == S: 

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

1933 

1934 def critical_ps(self, mode): 

1935 _, ubelow = self.u_top_bottom(mode) 

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

1937 

1938 def pflat(self, p): 

1939 return p / (earthradius-self.z) 

1940 

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

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

1943 return 0.0 

1944 else: 

1945 return psv_surface( 

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

1947 psv_surface_ind(in_mode, out_mode)] 

1948 

1949 def __str__(self): 

1950 return 'surface' 

1951 

1952 

1953class Walker(object): 

1954 def __init__(self, elements): 

1955 self._elements = elements 

1956 self._i = 0 

1957 

1958 def current(self): 

1959 return self._elements[self._i] 

1960 

1961 def go(self, direction): 

1962 if direction == UP: 

1963 self.up() 

1964 else: 

1965 self.down() 

1966 

1967 def down(self): 

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

1969 self._i += 1 

1970 else: 

1971 raise BottomReached() 

1972 

1973 def up(self): 

1974 if self._i > 0: 

1975 self._i -= 1 

1976 else: 

1977 raise SurfaceReached() 

1978 

1979 def goto_layer(self, layer): 

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

1981 

1982 

1983class RayElement(object): 

1984 ''' 

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

1986 ''' 

1987 

1988 def __eq__(self, other): 

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

1990 

1991 def is_straight(self): 

1992 return isinstance(self, Straight) 

1993 

1994 def is_kink(self): 

1995 return isinstance(self, Kink) 

1996 

1997 

1998class Straight(RayElement): 

1999 ''' 

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

2001 ''' 

2002 

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

2004 self.mode = mode 

2005 self._direction_in = direction_in 

2006 self._direction_out = direction_out 

2007 self.layer = layer 

2008 

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

2010 z = self.z_in(endgaps) 

2011 dir = self.eff_direction_in(endgaps) 

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

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

2014 

2015 if dir == DOWN: 

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

2017 else: 

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

2019 

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

2021 z = self.z_out(endgaps) 

2022 dir = self.eff_direction_out(endgaps) 

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

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

2025 

2026 if dir == DOWN: 

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

2028 else: 

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

2030 

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

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

2033 

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

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

2036 

2037 def test(self, p, z): 

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

2039 

2040 def z_in(self, endgaps=None): 

2041 if endgaps is not None: 

2042 return endgaps[0] 

2043 else: 

2044 lyr = self.layer 

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

2046 

2047 def z_out(self, endgaps=None): 

2048 if endgaps is not None: 

2049 return endgaps[1] 

2050 else: 

2051 lyr = self.layer 

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

2053 

2054 def turns(self): 

2055 return self._direction_in != self._direction_out 

2056 

2057 def eff_direction_in(self, endgaps=None): 

2058 if endgaps is None: 

2059 return self._direction_in 

2060 else: 

2061 return endgaps[2] 

2062 

2063 def eff_direction_out(self, endgaps=None): 

2064 if endgaps is None: 

2065 return self._direction_out 

2066 else: 

2067 return endgaps[3] 

2068 

2069 def zturn(self, p): 

2070 lyr = self.layer 

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

2072 

2073 def u_in(self, endgaps=None): 

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

2075 

2076 def u_out(self, endgaps=None): 

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

2078 

2079 def critical_p_in(self, endgaps=None): 

2080 z = self.z_in(endgaps) 

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

2082 

2083 def critical_p_out(self, endgaps=None): 

2084 z = self.z_out(endgaps) 

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

2086 

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

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

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

2090 x *= 2. 

2091 t *= 2. 

2092 return x, t 

2093 

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

2095 z1, z2 = zstart, zstop 

2096 if z1 > z2: 

2097 z1, z2 = z2, z1 

2098 

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

2100 

2101 if samedir: 

2102 return x, t 

2103 else: 

2104 xfull, tfull = self.xt(p) 

2105 return xfull-x, tfull-t 

2106 

2107 def __hash__(self): 

2108 return hash(( 

2109 self._direction_in, 

2110 self._direction_out, 

2111 self.mode, 

2112 id(self.layer))) 

2113 

2114 

2115class HeadwaveStraight(Straight): 

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

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

2118 

2119 self.interface = interface 

2120 

2121 def z_in(self, zpart=None): 

2122 return self.interface.z 

2123 

2124 def z_out(self, zpart=None): 

2125 return self.interface.z 

2126 

2127 def zturn(self, p): 

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

2129 

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

2131 return 0., 0. 

2132 

2133 def x2t_headwave(self, xstretch): 

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

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

2136 

2137 

2138class Kink(RayElement): 

2139 ''' 

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

2141 ''' 

2142 

2143 def __init__( 

2144 self, 

2145 in_direction, 

2146 out_direction, 

2147 in_mode, 

2148 out_mode, 

2149 discontinuity): 

2150 

2151 self.in_direction = in_direction 

2152 self.out_direction = out_direction 

2153 self.in_mode = in_mode 

2154 self.out_mode = out_mode 

2155 self.discontinuity = discontinuity 

2156 

2157 def reflection(self): 

2158 return self.in_direction != self.out_direction 

2159 

2160 def conversion(self): 

2161 return self.in_mode != self.out_mode 

2162 

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

2164 

2165 if out_direction is None: 

2166 out_direction = self.out_direction 

2167 

2168 if out_mode is None: 

2169 out_mode = self.out_mode 

2170 

2171 return self.discontinuity.efficiency( 

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

2173 

2174 def __str__(self): 

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

2176 if r and c: 

2177 return '|~' 

2178 if r: 

2179 return '|' 

2180 if c: 

2181 return '~' 

2182 return '_' 

2183 

2184 def __hash__(self): 

2185 return hash(( 

2186 self.in_direction, 

2187 self.out_direction, 

2188 self.in_mode, 

2189 self.out_mode, 

2190 id(self.discontinuity))) 

2191 

2192 

2193class PRangeNotSet(CakeError): 

2194 pass 

2195 

2196 

2197class RayPath(object): 

2198 ''' 

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

2200 layers / interfaces. 

2201 ''' 

2202 

2203 def __init__(self, phase): 

2204 self.elements = [] 

2205 self.phase = phase 

2206 self._pmax = None 

2207 self._pmin = None 

2208 self._p = None 

2209 self._is_headwave = False 

2210 

2211 def set_is_headwave(self, is_headwave): 

2212 self._is_headwave = is_headwave 

2213 

2214 def copy(self): 

2215 ''' 

2216 Get a copy of it. 

2217 ''' 

2218 

2219 c = copy.copy(self) 

2220 c.elements = list(self.elements) 

2221 return c 

2222 

2223 def endgaps(self, zstart, zstop): 

2224 ''' 

2225 Get information needed for end point adjustments. 

2226 ''' 

2227 

2228 return ( 

2229 zstart, 

2230 zstop, 

2231 self.phase.direction_start(), 

2232 self.phase.direction_stop()) 

2233 

2234 def append(self, element): 

2235 self.elements.append(element) 

2236 

2237 def _check_have_prange(self): 

2238 if self._pmax is None: 

2239 raise PRangeNotSet() 

2240 

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

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

2243 self._prange_dp = dp 

2244 

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

2246 ''' 

2247 Calculate phase definition from ray path. 

2248 ''' 

2249 

2250 used = PhaseDef() 

2251 fleg = self.phase.first_leg() 

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

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

2254 for before, element, after in zip( 

2255 n_elements_n[:-2], 

2256 n_elements_n[1:-1], 

2257 n_elements_n[2:]): 

2258 

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

2260 type(before), 

2261 type(after)): 

2262 

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

2264 z = element.discontinuity.z 

2265 used.append(Knee( 

2266 z, 

2267 element.in_direction, 

2268 element.out_direction != element.in_direction, 

2269 element.in_mode, 

2270 element.out_mode)) 

2271 

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

2273 

2274 elif type(element) is HeadwaveStraight: 

2275 z = element.interface.z 

2276 k = Knee( 

2277 z, 

2278 before.in_direction, 

2279 after.out_direction != before.in_direction, 

2280 before.in_mode, 

2281 after.out_mode) 

2282 

2283 k.headwave = True 

2284 used.append(k) 

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

2286 

2287 if (p is not None and before and after 

2288 and element.is_straight() 

2289 and before.is_kink() 

2290 and after.is_kink() 

2291 and element.turns() 

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

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

2294 

2295 ai = element.angle_in(p) 

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

2297 used.append( 

2298 Head( 

2299 before.discontinuity.z, 

2300 before.out_direction, 

2301 element.mode)) 

2302 used.append( 

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

2304 

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

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

2307 

2308 return used 

2309 

2310 def pmax(self): 

2311 ''' 

2312 Get maximum valid ray parameter. 

2313 ''' 

2314 self._check_have_prange() 

2315 return self._pmax 

2316 

2317 def pmin(self): 

2318 ''' 

2319 Get minimum valid ray parameter. 

2320 ''' 

2321 self._check_have_prange() 

2322 return self._pmin 

2323 

2324 def xmin(self): 

2325 ''' 

2326 Get minimal distance. 

2327 ''' 

2328 self._analyse() 

2329 return self._xmin 

2330 

2331 def xmax(self): 

2332 ''' 

2333 Get maximal distance. 

2334 ''' 

2335 self._analyse() 

2336 return self._xmax 

2337 

2338 def kinks(self): 

2339 ''' 

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

2341 ''' 

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

2343 

2344 def straights(self): 

2345 ''' 

2346 Iterate over ray segments. 

2347 ''' 

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

2349 

2350 def headwave_straight(self): 

2351 for s in self.elements: 

2352 if type(s) is HeadwaveStraight: 

2353 return s 

2354 

2355 def first_straight(self): 

2356 ''' 

2357 Get first ray segment. 

2358 ''' 

2359 for s in self.elements: 

2360 if isinstance(s, Straight): 

2361 return s 

2362 

2363 def last_straight(self): 

2364 ''' 

2365 Get last ray segment. 

2366 ''' 

2367 for s in reversed(self.elements): 

2368 if isinstance(s, Straight): 

2369 return s 

2370 

2371 def efficiency(self, p): 

2372 ''' 

2373 Get product of all conversion/reflection coefficients encountered on 

2374 path. 

2375 ''' 

2376 return reduce( 

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

2378 

2379 def spreading(self, p, endgaps): 

2380 ''' 

2381 Get geometrical spreading factor. 

2382 ''' 

2383 if self._is_headwave: 

2384 return 0.0 

2385 

2386 self._check_have_prange() 

2387 dp = self._prange_dp * 0.01 

2388 assert self._pmax - self._pmin > dp 

2389 

2390 if p + dp > self._pmax: 

2391 p = p-dp 

2392 

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

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

2395 x0 *= d2r 

2396 x1 *= d2r 

2397 if x1 == x0: 

2398 return num.nan 

2399 

2400 dp_dx = dp/(x1-x0) 

2401 

2402 x = x0 

2403 if x == 0.: 

2404 x = x1 

2405 p = dp 

2406 

2407 first = self.first_straight() 

2408 last = self.last_straight() 

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

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

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

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

2413 first.u_in(endgaps)**2 * 

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

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

2416 

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

2418 assert dp is None or n is None 

2419 

2420 if self._pmin == self._pmax: 

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

2422 

2423 if dp is None: 

2424 dp = self._prange_dp 

2425 

2426 if n is None: 

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

2428 

2429 if nmin is not None: 

2430 n = max(n, nmin) 

2431 

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

2433 return ppp 

2434 

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

2436 ''' 

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

2438 path's ends. 

2439 ''' 

2440 

2441 zstart, zstop, dirstart, dirstop = endgaps 

2442 firsts = self.first_straight() 

2443 lasts = self.last_straight() 

2444 xs, ts = firsts.xt_gap( 

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

2446 xe, te = lasts.xt_gap( 

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

2448 

2449 if which == 'both': 

2450 return xs + xe, ts + te 

2451 elif which == 'left': 

2452 return xs, ts 

2453 elif which == 'right': 

2454 return xe, te 

2455 

2456 def xt_endgaps_ptest(self, p, endgaps): 

2457 ''' 

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

2459 ''' 

2460 

2461 zstart, zstop, dirstart, dirstop = endgaps 

2462 firsts = self.first_straight() 

2463 lasts = self.last_straight() 

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

2465 

2466 def xt(self, p, endgaps): 

2467 ''' 

2468 Calculate distance and traveltime for given ray parameter. 

2469 ''' 

2470 

2471 if isinstance(p, num.ndarray): 

2472 sx = num.zeros(p.size) 

2473 st = num.zeros(p.size) 

2474 else: 

2475 sx = 0.0 

2476 st = 0.0 

2477 

2478 for s in self.straights(): 

2479 x, t = s.xt(p) 

2480 sx += x 

2481 st += t 

2482 

2483 if endgaps: 

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

2485 sx -= dx 

2486 st -= dt 

2487 

2488 return sx, st 

2489 

2490 def xt_limits(self, p): 

2491 ''' 

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

2493 ''' 

2494 

2495 if isinstance(p, num.ndarray): 

2496 sx = num.zeros(p.size) 

2497 st = num.zeros(p.size) 

2498 sxe = num.zeros(p.size) 

2499 ste = num.zeros(p.size) 

2500 else: 

2501 sx = 0.0 

2502 st = 0.0 

2503 sxe = 0.0 

2504 ste = 0.0 

2505 

2506 sfirst = self.first_straight() 

2507 slast = self.last_straight() 

2508 

2509 for s in self.straights(): 

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

2511 x, t = s.xt(p) 

2512 sx += x 

2513 st += t 

2514 

2515 sends = [sfirst] 

2516 if sfirst is not slast: 

2517 sends.append(slast) 

2518 

2519 for s in sends: 

2520 x, t = s.xt(p) 

2521 sxe += x 

2522 ste += t 

2523 

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

2525 

2526 def iter_zxt(self, p): 

2527 ''' 

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

2529 ray path. 

2530 ''' 

2531 

2532 sx = num.zeros(p.size) 

2533 st = num.zeros(p.size) 

2534 ok = False 

2535 for s in self.straights(): 

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

2537 

2538 x, t = s.xt(p) 

2539 sx += x 

2540 st += t 

2541 ok = True 

2542 

2543 if ok: 

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

2545 

2546 def zxt_path_subdivided( 

2547 self, p, endgaps, 

2548 points_per_straight=20, 

2549 x_for_headwave=None): 

2550 

2551 ''' 

2552 Get geometrical representation of ray path. 

2553 ''' 

2554 

2555 if self._is_headwave: 

2556 assert p.size == 1 

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

2558 xstretch = x_for_headwave-x 

2559 nout = xstretch.size 

2560 else: 

2561 nout = p.size 

2562 

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

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

2565 

2566 # first create full path including the endgaps 

2567 sx = num.zeros(nout) - dxl 

2568 st = num.zeros(nout) - dtl 

2569 zxt = [] 

2570 for s in self.straights(): 

2571 n = points_per_straight 

2572 

2573 back = None 

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

2575 if type(s) is HeadwaveStraight: 

2576 z = zin 

2577 for i in range(n): 

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

2579 ts = s.x2t_headwave(xs) 

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

2581 else: 

2582 if zin != zout: # normal traversal 

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

2584 for z in zs: 

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

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

2587 

2588 else: # ray turns in layer 

2589 zturn = s.zturn(p) 

2590 back = [] 

2591 for i in range(n): 

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

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

2594 

2595 if zturn[0] >= zin: 

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

2597 else: 

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

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

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

2601 

2602 if type(s) is HeadwaveStraight: 

2603 x = xstretch 

2604 t = s.x2t_headwave(xstretch) 

2605 else: 

2606 x, t = s.xt(p) 

2607 

2608 sx += x 

2609 st += t 

2610 if back: 

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

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

2613 

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

2615 fanz, fanx, fant = [], [], [] 

2616 for z, x, t in zxt: 

2617 fanz.append(z) 

2618 fanx.append(x) 

2619 fant.append(t) 

2620 

2621 z = num.array(fanz).T 

2622 x = num.array(fanx).T 

2623 t = num.array(fant).T 

2624 

2625 # cut off the endgaps, add exact endpoints 

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

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

2628 zstart, zstop = endgaps[:2] 

2629 zs, xs, ts = [], [], [] 

2630 for i in range(nout): 

2631 t_ = t[i] 

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

2633 n = indices.size + 2 

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

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

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

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

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

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

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

2641 zs.append(zs_) 

2642 xs.append(xs_) 

2643 ts.append(ts_) 

2644 

2645 return zs, xs, ts 

2646 

2647 def _analyse(self): 

2648 if self._p is not None: 

2649 return 

2650 

2651 p = self.make_p(nmin=20) 

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

2653 

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

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

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

2657 

2658 def draft_pxt(self, endgaps): 

2659 self._analyse() 

2660 

2661 if not self._is_headwave: 

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

2663 pcrit = min( 

2664 self.critical_pstart(endgaps), 

2665 self.critical_pstop(endgaps)) 

2666 

2667 if pcrit < self._pmin: 

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

2669 return empty, empty, empty 

2670 

2671 elif pcrit >= self._pmax: 

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

2673 return cp, cx-dx, ct-dt 

2674 

2675 else: 

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

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

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

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

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

2681 rp[-1] = pcrit 

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

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

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

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

2686 return rp, rx, rt 

2687 

2688 else: 

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

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

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

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

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

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

2695 

2696 def interpolate_x2pt_linear(self, x, endgaps): 

2697 ''' 

2698 Get approximate ray parameter and traveltime for distance. 

2699 ''' 

2700 

2701 self._analyse() 

2702 

2703 if self._is_headwave: 

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

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

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

2707 el = self.headwave_straight() 

2708 xok = x[x >= xmin] 

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

2710 return [ 

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

2712 

2713 else: 

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

2715 return [] 

2716 

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

2718 

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

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

2721 

2722 if (rp.size and 

2723 len(xp) == 0 and 

2724 rx[0] == 0.0 and 

2725 any(x == 0.0) and 

2726 rp[0] == 0.0): 

2727 

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

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

2730 

2731 return [ 

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

2733 

2734 def __eq__(self, other): 

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

2736 return False 

2737 

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

2739 

2740 def __hash__(self): 

2741 return hash( 

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

2743 (self.phase.definition(), )) 

2744 

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

2746 x = [] 

2747 start_i = None 

2748 end_i = None 

2749 turn_i = None 

2750 

2751 def append_layers(si, ei, ti): 

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

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

2754 else: 

2755 if ti is not None: 

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

2757 else: 

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

2759 

2760 for el in self.elements: 

2761 if type(el) is Straight: 

2762 if start_i is None: 

2763 start_i = el.layer.ilayer 

2764 if el._direction_in != el._direction_out: 

2765 turn_i = el.layer.ilayer 

2766 end_i = el.layer.ilayer 

2767 

2768 elif isinstance(el, Kink): 

2769 if start_i is not None: 

2770 append_layers(start_i, end_i, turn_i) 

2771 start_i = None 

2772 turn_i = None 

2773 

2774 x.append(str(el)) 

2775 

2776 if start_i is not None: 

2777 append_layers(start_i, end_i, turn_i) 

2778 

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

2780 

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

2782 

2783 def critical_pstart(self, endgaps): 

2784 ''' 

2785 Get critical ray parameter for source depth choice. 

2786 ''' 

2787 

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

2789 

2790 def critical_pstop(self, endgaps): 

2791 ''' 

2792 Get critical ray parameter for receiver depth choice. 

2793 ''' 

2794 

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

2796 

2797 def ranges(self, endgaps): 

2798 ''' 

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

2800 ''' 

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

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

2803 

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

2805 ''' 

2806 Get textual representation. 

2807 ''' 

2808 

2809 self._analyse() 

2810 

2811 if as_degrees: 

2812 xunit = 'deg' 

2813 xfact = 1. 

2814 else: 

2815 xunit = 'km' 

2816 xfact = d2m/km 

2817 

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

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

2820 - t [%g, %g] s 

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

2822''' % ( 

2823 self._xmin*xfact, 

2824 self._xmax*xfact, 

2825 xunit, 

2826 self._tmin, 

2827 self._tmax, 

2828 self._pmin/r2d, 

2829 self._pmax/r2d) 

2830 

2831 if endgaps is not None: 

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

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

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

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

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

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

2838 

2839 else: 

2840 ss = '' 

2841 

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

2843 

2844 

2845class RefineFailed(CakeError): 

2846 pass 

2847 

2848 

2849class Ray(object): 

2850 ''' 

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

2852 arrival time) choice. 

2853 

2854 **Attributes:** 

2855 

2856 .. py:attribute:: path 

2857 

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

2859 

2860 .. py:attribute:: p 

2861 

2862 Ray parameter (spherical) [s/rad] 

2863 

2864 .. py:attribute:: x 

2865 

2866 Radial distance [deg] 

2867 

2868 .. py:attribute:: t 

2869 

2870 Traveltime [s] 

2871 

2872 .. py:attribute:: endgaps 

2873 

2874 Needed for source/receiver depth adjustments in many 

2875 :py:class:`RayPath` methods. 

2876 ''' 

2877 

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

2879 self.path = path 

2880 self.p = p 

2881 self.x = x 

2882 self.t = t 

2883 self.endgaps = endgaps 

2884 self.draft_pxt = draft_pxt 

2885 

2886 def given_phase(self): 

2887 ''' 

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

2889 

2890 :returns: :py:class:`PhaseDef` object 

2891 ''' 

2892 

2893 return self.path.phase 

2894 

2895 def used_phase(self): 

2896 ''' 

2897 Compute phase definition from propagation path. 

2898 

2899 :returns: :py:class:`PhaseDef` object 

2900 ''' 

2901 

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

2903 

2904 def refine(self): 

2905 if self.path._is_headwave: 

2906 return 

2907 

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

2909 return 

2910 

2911 cp, cx, ct = self.draft_pxt 

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

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

2914 raise RefineFailed() 

2915 

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

2917 p_to_t = {} 

2918 i = [0] 

2919 

2920 def f(p): 

2921 i[0] += 1 

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

2923 p_to_t[p] = t 

2924 return self.x - x 

2925 

2926 try: 

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

2928 self.t = p_to_t[self.p] 

2929 

2930 except ValueError: 

2931 raise RefineFailed() 

2932 

2933 def takeoff_angle(self): 

2934 ''' 

2935 Get takeoff angle of ray. 

2936 

2937 The angle is returned in [degrees]. 

2938 ''' 

2939 

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

2941 

2942 def incidence_angle(self): 

2943 ''' 

2944 Get incidence angle of ray. 

2945 

2946 The angle is returned in [degrees]. 

2947 ''' 

2948 

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

2950 

2951 def efficiency(self): 

2952 ''' 

2953 Get conversion/reflection efficiency of the ray. 

2954 

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

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

2957 or conversions. 

2958 ''' 

2959 

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

2961 

2962 def spreading(self): 

2963 ''' 

2964 Get geometrical spreading factor. 

2965 ''' 

2966 

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

2968 

2969 def surface_sphere(self): 

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

2971 r2 = earthradius - self.endgaps[1] 

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

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

2974 

2975 def zxt_path_subdivided(self, points_per_straight=20): 

2976 ''' 

2977 Get geometrical representation of ray path. 

2978 

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

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

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

2982 ``points_per_straight`` parameter. 

2983 ''' 

2984 return self.path.zxt_path_subdivided( 

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

2986 points_per_straight=points_per_straight, 

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

2988 

2989 def __str__(self, as_degrees=False): 

2990 if as_degrees: 

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

2992 else: 

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

2994 

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

2996 self.p/r2d, 

2997 sd, 

2998 self.t, 

2999 self.takeoff_angle(), 

3000 self.incidence_angle(), 

3001 100*self.efficiency(), 

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

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

3004 

3005 

3006def anything_to_crust2_profile(crust2_profile): 

3007 from pyrocko.dataset import crust2x2 

3008 if isinstance(crust2_profile, tuple): 

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

3010 return crust2x2.get_profile(lat, lon) 

3011 elif isinstance(crust2_profile, str): 

3012 return crust2x2.get_profile(crust2_profile) 

3013 elif isinstance(crust2_profile, crust2x2.Crust2Profile): 

3014 return crust2_profile 

3015 else: 

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

3017 'key or a crust2x2 Profile object)' 

3018 

3019 

3020class DiscontinuityNotFound(CakeError): 

3021 def __init__(self, depth_or_name): 

3022 CakeError.__init__(self) 

3023 self.depth_or_name = depth_or_name 

3024 

3025 def __str__(self): 

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

3027 self.depth_or_name 

3028 

3029 

3030class LayeredModelError(CakeError): 

3031 pass 

3032 

3033 

3034class LayeredModel(object): 

3035 ''' 

3036 Representation of a layer cake model. 

3037 

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

3039 

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

3041 file. 

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

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

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

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

3046 objects from a given velocity profile. 

3047 

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

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

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

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

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

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

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

3055 are cached for reuse in successive invocations. 

3056 ''' 

3057 

3058 def __init__(self): 

3059 self._surface_material = None 

3060 self._elements = [] 

3061 self.nlayers = 0 

3062 self._np = 10000 

3063 self._pdepth = 5 

3064 self._pathcache = {} 

3065 

3066 def copy_with_elevation(self, elevation): 

3067 ''' 

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

3069 

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

3071 

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

3073 `z` axis. 

3074 ''' 

3075 

3076 c = copy.deepcopy(self) 

3077 c._pathcache = {} 

3078 surface = c._elements[0] 

3079 toplayer = c._elements[1] 

3080 

3081 assert toplayer.zbot > -elevation 

3082 

3083 surface.z = -elevation 

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

3085 c._elements[1].ilayer = 0 

3086 return c 

3087 

3088 def zeq(self, z1, z2): 

3089 return abs(z1-z2) < ZEPS 

3090 

3091 def append(self, element): 

3092 ''' 

3093 Add a layer or discontinuity at bottom of model. 

3094 

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

3096 :py:class:`Discontinuity`. 

3097 ''' 

3098 

3099 if isinstance(element, Layer): 

3100 if element.zbot >= earthradius: 

3101 element.zbot = earthradius - 1. 

3102 

3103 if element.ztop >= earthradius: 

3104 raise CakeError('Layer deeper than earthradius') 

3105 

3106 element.ilayer = self.nlayers 

3107 self.nlayers += 1 

3108 

3109 self._elements.append(element) 

3110 

3111 def elements(self, direction=DOWN): 

3112 ''' 

3113 Iterate over all elements of the model. 

3114 

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

3116 :py:const:`UP`. 

3117 

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

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

3120 ''' 

3121 

3122 if direction == DOWN: 

3123 return iter(self._elements) 

3124 else: 

3125 return reversed(self._elements) 

3126 

3127 def layers(self, direction=DOWN): 

3128 ''' 

3129 Iterate over all layers of model. 

3130 

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

3132 :py:const:`UP`. 

3133 

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

3135 ''' 

3136 

3137 if direction == DOWN: 

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

3139 else: 

3140 return ( 

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

3142 

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

3144 ''' 

3145 Get layer for given depth. 

3146 

3147 :param z: depth [m] 

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

3149 :py:const:`UP`. 

3150 

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

3152 ''' 

3153 

3154 for layer in self.layers(direction): 

3155 if layer.contains(z): 

3156 return layer 

3157 else: 

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

3159 

3160 def walker(self): 

3161 return Walker(self._elements) 

3162 

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

3164 ''' 

3165 Get material at given depth. 

3166 

3167 :param z: depth [m] 

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

3169 :py:const:`UP` 

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

3171 

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

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

3174 ''' 

3175 

3176 lyr = self.layer(z, direction) 

3177 return lyr.material(z) 

3178 

3179 def discontinuities(self): 

3180 ''' 

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

3182 

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

3184 

3185 def discontinuity(self, name_or_z): 

3186 ''' 

3187 Get discontinuity by name or depth. 

3188 

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

3190 ''' 

3191 

3192 if isinstance(name_or_z, float): 

3193 candi = sorted( 

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

3195 else: 

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

3197 

3198 if not candi: 

3199 raise DiscontinuityNotFound(name_or_z) 

3200 

3201 return candi[0] 

3202 

3203 def adapt_phase(self, phase): 

3204 ''' 

3205 Adapt a phase definition for use with this model. 

3206 

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

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

3209 in the model. 

3210 ''' 

3211 

3212 phase = phase.copy() 

3213 for knee in phase.knees(): 

3214 if knee.depth != 'surface': 

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

3216 for leg in phase.legs(): 

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

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

3219 

3220 return phase 

3221 

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

3223 ''' 

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

3225 source and receiver layers. 

3226 

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

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

3229 :param layer_start: layer with source 

3230 :param layer_stop: layer with receiver 

3231 :returns: :py:class:`RayPath` object 

3232 

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

3234 :py:exc:`NotPhaseConform`, :py:exc:`MinDepthReached`, 

3235 :py:exc:`MaxDepthReached`, :py:exc:`CannotPropagate`, 

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

3237 ''' 

3238 

3239 phase = self.adapt_phase(phase) 

3240 knees = phase.knees() 

3241 legs = phase.legs() 

3242 next_knee = next_or_none(knees) 

3243 leg = next_or_none(legs) 

3244 assert leg is not None 

3245 

3246 direction = leg.departure 

3247 direction_stop = phase.direction_stop() 

3248 mode = leg.mode 

3249 mode_stop = phase.last_leg().mode 

3250 

3251 walker = self.walker() 

3252 walker.goto_layer(layer_start) 

3253 current = walker.current() 

3254 

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

3256 if not ttop and not tbot: 

3257 raise CannotPropagate(direction, current.ilayer) 

3258 

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

3260 direction = -direction 

3261 

3262 path = RayPath(phase) 

3263 trapdetect = set() 

3264 while True: 

3265 at_layer = isinstance(current, Layer) 

3266 at_discontinuity = isinstance(current, Discontinuity) 

3267 

3268 # detect trapped wave 

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

3270 if k in trapdetect: 

3271 raise Trapped() 

3272 

3273 trapdetect.add(k) 

3274 

3275 if at_discontinuity: 

3276 oldmode, olddirection = mode, direction 

3277 headwave = False 

3278 if next_knee is not None and next_knee.matches( 

3279 current, mode, direction): 

3280 

3281 headwave = next_knee.headwave 

3282 direction = next_knee.out_direction() 

3283 mode = next_knee.out_mode 

3284 next_knee = next_or_none(knees) 

3285 leg = next(legs) 

3286 

3287 else: # implicit reflection/transmission 

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

3289 

3290 if headwave: 

3291 path.set_is_headwave(True) 

3292 

3293 path.append(Kink( 

3294 olddirection, olddirection, oldmode, oldmode, current)) 

3295 

3296 path.append(HeadwaveStraight( 

3297 olddirection, direction, oldmode, current)) 

3298 

3299 path.append(Kink( 

3300 olddirection, direction, oldmode, mode, current)) 

3301 

3302 else: 

3303 path.append(Kink( 

3304 olddirection, direction, oldmode, mode, current)) 

3305 

3306 if at_layer: 

3307 direction_in = direction 

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

3309 

3310 zturn = None 

3311 if direction_in != direction: 

3312 zturn = current.zturn(p, mode) 

3313 

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

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

3316 if direction_in != direction: 

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

3318 raise MinDepthReached() 

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

3320 raise MaxDepthReached() 

3321 else: 

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

3323 raise MinDepthReached() 

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

3325 raise MaxDepthReached() 

3326 

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

3328 

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

3330 current is layer_stop: 

3331 

3332 if zturn is None: 

3333 if direction == direction_stop: 

3334 break 

3335 else: 

3336 break 

3337 

3338 walker.go(direction) 

3339 current = walker.current() 

3340 

3341 return path 

3342 

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

3344 ''' 

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

3346 or more phase definitions. 

3347 

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

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

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

3351 :param zstart: source depth [m] 

3352 :param zstop: receiver depth [m] 

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

3354 

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

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

3357 phase definition has been used before. 

3358 ''' 

3359 

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

3361 

3362 phases = to_phase_defs(phases) 

3363 

3364 paths = [] 

3365 for phase in phases: 

3366 

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

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

3369 

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

3371 

3372 if pathcachekey in self._pathcache: 

3373 phase_paths = self._pathcache[pathcachekey] 

3374 else: 

3375 hwknee = phase.headwave_knee() 

3376 if hwknee: 

3377 name_or_z = hwknee.depth 

3378 interface = self.discontinuity(name_or_z) 

3379 mode = hwknee.in_mode 

3380 in_direction = hwknee.direction 

3381 

3382 pabove, pbelow = interface.critical_ps(mode) 

3383 

3384 p = min_not_none(pabove, pbelow) 

3385 

3386 # diffracted wave: 

3387 if in_direction == DOWN and ( 

3388 pbelow is None or pbelow >= pabove): 

3389 

3390 p *= (1.0 - eps) 

3391 

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

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

3394 

3395 phase_paths = [path] 

3396 

3397 else: 

3398 try: 

3399 pmax_start = max([ 

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

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

3402 

3403 pmax_stop = max([ 

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

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

3406 

3407 pmax = min(pmax_start, pmax_stop) 

3408 

3409 pedges = [0.] 

3410 for layer in self.layers(): 

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

3412 for mode in (P, S): 

3413 for eps2 in [eps]: 

3414 v = layer.v(mode, z) 

3415 if v != 0.0: 

3416 p = radius(z)/v 

3417 if p <= pmax: 

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

3419 pedges.append(p) 

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

3421 

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

3423 

3424 phase_paths = {} 

3425 cached = {} 

3426 counter = [0] 

3427 

3428 def p_to_path(p): 

3429 if p in cached: 

3430 return cached[p] 

3431 

3432 try: 

3433 counter[0] += 1 

3434 path = self.path( 

3435 p, phase, layer_start, layer_stop) 

3436 

3437 if path not in phase_paths: 

3438 phase_paths[path] = [] 

3439 

3440 phase_paths[path].append(p) 

3441 

3442 except PathFailed: 

3443 path = None 

3444 

3445 cached[p] = path 

3446 return path 

3447 

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

3449 if i > self._pdepth: 

3450 return 

3451 path1 = p_to_path(pmin) 

3452 path2 = p_to_path(pmax) 

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

3454 return 

3455 if path1 is None or path2 is None or \ 

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

3457 

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

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

3460 

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

3462 recurse(pl, ph) 

3463 

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

3465 path.set_prange( 

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

3467 

3468 phase_paths = list(phase_paths.keys()) 

3469 

3470 except ZeroDivisionError: 

3471 phase_paths = [] 

3472 

3473 self._pathcache[pathcachekey] = phase_paths 

3474 

3475 paths.extend(phase_paths) 

3476 

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

3478 return paths 

3479 

3480 def arrivals( 

3481 self, 

3482 distances=[], 

3483 phases=PhaseDef('P'), 

3484 zstart=0.0, 

3485 zstop=0.0, 

3486 refine=True): 

3487 

3488 ''' 

3489 Compute rays and traveltimes for given distances. 

3490 

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

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

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

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

3495 :param zstart: source depth [m] 

3496 :param zstop: receiver depth [m] 

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

3498 (p, x, t) estimated from interpolation 

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

3500 (distance, arrival time) 

3501 ''' 

3502 

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

3504 

3505 arrivals = [] 

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

3507 

3508 endgaps = path.endgaps(zstart, zstop) 

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

3510 distances, endgaps): 

3511 

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

3513 

3514 if refine: 

3515 refined = [] 

3516 for ray in arrivals: 

3517 

3518 if ray.path._is_headwave: 

3519 refined.append(ray) 

3520 

3521 try: 

3522 ray.refine() 

3523 refined.append(ray) 

3524 

3525 except RefineFailed: 

3526 pass 

3527 

3528 arrivals = refined 

3529 

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

3531 return arrivals 

3532 

3533 @classmethod 

3534 def from_scanlines(cls, producer): 

3535 ''' 

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

3537 

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

3539 

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

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

3542 ''' 

3543 

3544 self = cls() 

3545 for z, material, name in producer: 

3546 

3547 if not self._elements: 

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

3549 else: 

3550 element = self._elements[-1] 

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

3552 assert isinstance(element, Layer) 

3553 self.append( 

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

3555 

3556 else: 

3557 if isinstance(element, Discontinuity): 

3558 ztop = element.z 

3559 mtop = element.mbelow 

3560 elif isinstance(element, Layer): 

3561 ztop = element.zbot 

3562 mtop = element.mbot 

3563 

3564 if mtop == material: 

3565 layer = HomogeneousLayer( 

3566 ztop, z, material, name=name) 

3567 else: 

3568 layer = GradientLayer( 

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

3570 

3571 self.append(layer) 

3572 

3573 return self 

3574 

3575 def to_scanlines(self, get_burgers=False): 

3576 def fmt(z, m): 

3577 if not m._has_default_burgers() or get_burgers: 

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

3579 m.burger_eta1, m.burger_eta2, m.burger_valpha) 

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

3581 

3582 last = None 

3583 lines = [] 

3584 for element in self.elements(): 

3585 if isinstance(element, Layer): 

3586 if not isinstance(last, Layer): 

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

3588 

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

3590 

3591 last = element 

3592 

3593 if not isinstance(last, Layer): 

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

3595 

3596 return lines 

3597 

3598 def iter_material_parameter(self, get): 

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

3600 if get == 'z': 

3601 for layer in self.layers(): 

3602 yield layer.ztop 

3603 yield layer.zbot 

3604 else: 

3605 getter = operator.attrgetter(get) 

3606 for layer in self.layers(): 

3607 yield getter(layer.mtop) 

3608 yield getter(layer.mbot) 

3609 

3610 def profile(self, get): 

3611 ''' 

3612 Get parameter profile along depth of the earthmodel. 

3613 

3614 :param get: property to be queried ( 

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

3616 :type get: string 

3617 ''' 

3618 

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

3620 

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

3622 ''' 

3623 Find minimum value of a material property or depth. 

3624 

3625 :param get: property to be queried ( 

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

3627 ''' 

3628 

3629 return min(self.iter_material_parameter(get)) 

3630 

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

3632 ''' 

3633 Find maximum value of a material property or depth. 

3634 

3635 :param get: property to be queried ( 

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

3637 ''' 

3638 

3639 return max(self.iter_material_parameter(get)) 

3640 

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

3642 if len(layers) <= 1: 

3643 return layers 

3644 

3645 ztop = layers[0].ztop 

3646 zbot = layers[-1].zbot 

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

3648 zorigs.append(zbot) 

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

3650 data = [] 

3651 for z in zs: 

3652 if z == ztop: 

3653 direction = UP 

3654 else: 

3655 direction = DOWN 

3656 

3657 mat = self.material(z, direction) 

3658 data.append(mat.astuple()) 

3659 

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

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

3662 nmax = len(layers) // 2 

3663 accept = False 

3664 

3665 zcut_best = [] 

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

3667 ncutintervals = 20 

3668 zdelta = (zbot-ztop)/ncutintervals 

3669 if n == 2: 

3670 zcuts = [ 

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

3672 for i in range(1, ncutintervals)] 

3673 elif n == 3: 

3674 zcuts = [] 

3675 for j in range(1, ncutintervals): 

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

3677 zcuts.append( 

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

3679 else: 

3680 zcuts = [] 

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

3682 if zcut_best: 

3683 zcuts.append(sorted(num.linspace( 

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

3685 zcuts.append(sorted(num.linspace( 

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

3687 

3688 best = None 

3689 for icut, zcut in enumerate(zcuts): 

3690 rel_par_errors = num.zeros(5) 

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

3692 

3693 for ipar in range(5): 

3694 znodes, vnodes, error_rms = util.polylinefit( 

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

3696 

3697 mpar_nodes[:, ipar] = vnodes 

3698 if data_means[ipar] == 0.0: 

3699 rel_par_errors[ipar] = -1 

3700 else: 

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

3702 

3703 rel_error = rel_par_errors.max() 

3704 if best is None or rel_error < best[0]: 

3705 best = (rel_error, zcut, mpar_nodes) 

3706 

3707 rel_error, zcut, mpar_nodes = best 

3708 

3709 zcut_best.append(list(zcut)) 

3710 zcut_best[-1].pop(0) 

3711 zcut_best[-1].pop() 

3712 

3713 if rel_error <= max_rel_error: 

3714 accept = True 

3715 break 

3716 

3717 if not accept: 

3718 return layers 

3719 

3720 rel_error, zcut, mpar_nodes = best 

3721 

3722 material_nodes = [] 

3723 for i in range(n+1): 

3724 material_nodes.append(Material(*mpar_nodes[i, :])) 

3725 

3726 out_layers = [] 

3727 for i in range(n): 

3728 mtop = material_nodes[i] 

3729 mbot = material_nodes[i+1] 

3730 ztop = zcut[i] 

3731 zbot = zcut[i+1] 

3732 if mtop == mbot: 

3733 lyr = HomogeneousLayer(ztop, zbot, mtop) 

3734 else: 

3735 lyr = GradientLayer(ztop, zbot, mtop, mbot) 

3736 

3737 out_layers.append(lyr) 

3738 return out_layers 

3739 

3740 def simplify(self, max_rel_error=0.001): 

3741 ''' 

3742 Get representation of model with lower resolution. 

3743 

3744 Returns an approximation of the model. All discontinuities are kept, 

3745 but layer stacks with continuous model parameters are represented, if 

3746 possible, by a lower number of layers. Piecewise linear functions are 

3747 fitted against the original model parameter's piecewise linear 

3748 functions. Successively larger numbers of layers are tried, until the 

3749 difference to the original model is below ``max_rel_error``. The 

3750 difference is measured as the RMS error of the fit normalized by the 

3751 mean of the input (i.e. the fitted curves should deviate, on average, 

3752 less than 0.1% from the input curves if ``max_rel_error`` = 0.001). 

3753 ''' 

3754 

3755 mod_simple = LayeredModel() 

3756 

3757 glayers = [] 

3758 for element in self.elements(): 

3759 

3760 if isinstance(element, Discontinuity): 

3761 for layer in self.simplify_layers( 

3762 glayers, max_rel_error=max_rel_error): 

3763 

3764 mod_simple.append(layer) 

3765 

3766 glayers = [] 

3767 mod_simple.append(element) 

3768 else: 

3769 glayers.append(element) 

3770 

3771 for layer in self.simplify_layers( 

3772 glayers, max_rel_error=max_rel_error): 

3773 mod_simple.append(layer) 

3774 

3775 return mod_simple 

3776 

3777 def extract(self, depth_min=None, depth_max=None): 

3778 ''' 

3779 Extract :py:class:`LayeredModel` from :py:class:`LayeredModel`. 

3780 

3781 :param depth_min: depth of upper cut or name of :py:class:`Interface` 

3782 :param depth_max: depth of lower cut or name of :py:class:`Interface` 

3783 

3784 Interpolates a :py:class:`GradientLayer` at ``depth_min`` and/or 

3785 ``depth_max``. 

3786 ''' 

3787 

3788 if isinstance(depth_min, str): 

3789 depth_min = self.discontinuity(depth_min).z 

3790 

3791 if isinstance(depth_max, str): 

3792 depth_max = self.discontinuity(depth_max).z 

3793 

3794 mod_extracted = LayeredModel() 

3795 

3796 for element in self.elements(): 

3797 element = element.copy() 

3798 do_append = False 

3799 if (depth_min is None or depth_min <= element.ztop) \ 

3800 and (depth_max is None or depth_max >= element.zbot): 

3801 mod_extracted.append(element) 

3802 continue 

3803 

3804 if depth_min is not None: 

3805 if element.ztop < depth_min and depth_min < element.zbot: 

3806 _, element = element.split(depth_min) 

3807 do_append = True 

3808 

3809 if depth_max is not None: 

3810 if element.zbot > depth_max and depth_max > element.ztop: 

3811 element, _ = element.split(depth_max) 

3812 do_append = True 

3813 

3814 if do_append: 

3815 mod_extracted.append(element) 

3816 

3817 return mod_extracted 

3818 

3819 def replaced_crust(self, crust2_profile=None, crustmod=None): 

3820 if crust2_profile is not None: 

3821 profile = anything_to_crust2_profile(crust2_profile) 

3822 crustmod = LayeredModel.from_scanlines( 

3823 from_crust2x2_profile(profile)) 

3824 

3825 newmod = LayeredModel() 

3826 for element in crustmod.extract(depth_max='moho').elements(): 

3827 if element.name != 'moho': 

3828 newmod.append(element) 

3829 else: 

3830 moho1 = element 

3831 

3832 mod = self.extract(depth_min='moho') 

3833 first = True 

3834 for element in mod.elements(): 

3835 if element.name == 'moho': 

3836 if element.z <= moho1.z: 

3837 mbelow = mod.material(moho1.z, direction=UP) 

3838 else: 

3839 mbelow = element.mbelow 

3840 

3841 moho = Interface(moho1.z, moho1.mabove, mbelow, name='moho') 

3842 newmod.append(moho) 

3843 else: 

3844 if first: 

3845 if isinstance(element, Layer) and element.zbot > moho.z: 

3846 newmod.append(GradientLayer( 

3847 moho.z, 

3848 element.zbot, 

3849 moho.mbelow, 

3850 element.mbot, 

3851 name=element.name)) 

3852 

3853 first = False 

3854 else: 

3855 newmod.append(element) 

3856 return newmod 

3857 

3858 def perturb(self, rstate=None, keep_vp_vs=False, **kwargs): 

3859 ''' 

3860 Create a perturbed variant of the earth model. 

3861 

3862 Randomly change the thickness and material parameters of the earth 

3863 model from a uniform distribution. 

3864 

3865 :param kwargs: Maximum change fraction (e.g. 0.1) of the parameters. 

3866 Name the parameter, prefixed by ``p``. Supported parameters are 

3867 ``ph, pvp, pvs, prho, pqs, pqp``. 

3868 :type kwargs: dict 

3869 :param rstate: Random state to draw from, defaults to ``None`` 

3870 :type rstate: :class:`numpy.random.RandomState`, optional 

3871 :param keep_vp_vs: Keep the Vp/Vs ratio, defaults to False 

3872 :type keep_vp_vs: bool, optional 

3873 

3874 :returns: A new, perturbed earth model 

3875 :rtype: :class:`~pyrocko.cake.LayeredModel` 

3876 

3877 .. code-block :: python 

3878 

3879 perturbed_model = model.perturb(ph=.1, pvp=.05, prho=.1) 

3880 ''' 

3881 _pargs = set(['ph', 'pvp', 'pvs', 'prho', 'pqs', 'pqp']) 

3882 earthmod = copy.deepcopy(self) 

3883 

3884 if rstate is None: 

3885 rstate = num.random.RandomState() 

3886 

3887 layers = earthmod.layers() 

3888 discont = earthmod.discontinuities() 

3889 prev_layer = None 

3890 

3891 def get_change_ratios(): 

3892 values = dict.fromkeys([p[1:] for p in _pargs], 0.) 

3893 

3894 for param, pval in kwargs.items(): 

3895 if param not in _pargs: 

3896 continue 

3897 values[param[1:]] = float(rstate.uniform(-pval, pval, size=1)) 

3898 return values 

3899 

3900 # skip Surface 

3901 while True: 

3902 disc = next(discont) 

3903 if isinstance(disc, Surface): 

3904 break 

3905 

3906 while True: 

3907 try: 

3908 layer = next(layers) 

3909 m = layer.material(None) 

3910 h = layer.zbot - layer.ztop 

3911 except StopIteration: 

3912 break 

3913 

3914 if not isinstance(layer, HomogeneousLayer): 

3915 raise NotImplementedError( 

3916 'Can only perturbate homogeneous layers!') 

3917 

3918 changes = get_change_ratios() 

3919 

3920 # Changing thickness 

3921 dh = h * changes['h'] 

3922 changes['h'] = dh 

3923 

3924 layer.resize(depth_max=layer.zbot + dh, 

3925 depth_min=prev_layer.zbot if prev_layer else None) 

3926 

3927 try: 

3928 disc = next(discont) 

3929 disc.change_depth(disc.z + dh) 

3930 except StopIteration: 

3931 pass 

3932 

3933 # Setting material parameters 

3934 for param, change_ratio in changes.items(): 

3935 if param == 'h': 

3936 continue 

3937 

3938 value = m.__getattribute__(param) 

3939 changes[param] = value * change_ratio 

3940 

3941 if keep_vp_vs and changes['vp'] != 0.: 

3942 changes['vs'] = (m.vp + changes['vp']) / m.vp_vs_ratio() - m.vs 

3943 

3944 for param, change in changes.items(): 

3945 if param == 'h': 

3946 continue 

3947 value = m.__getattribute__(param) 

3948 m.__setattr__(param, value + change) 

3949 

3950 logger.info( 

3951 'perturbating earthmodel: {}'.format( 

3952 ' '.join(['{param}: {change:{len}.2f}'.format( 

3953 param=p, change=c, len=8) 

3954 for p, c in changes.items()]))) 

3955 

3956 prev_layer = layer 

3957 

3958 return earthmod 

3959 

3960 def require_homogeneous(self): 

3961 elements = list(self.elements()) 

3962 

3963 if len(elements) != 2: 

3964 raise LayeredModelError( 

3965 'Homogeneous model required but found more than one layer in ' 

3966 'earthmodel.') 

3967 

3968 if not isinstance(elements[1], HomogeneousLayer): 

3969 raise LayeredModelError( 

3970 'Homogeneous model required but element #1 is not of type ' 

3971 'HomogeneousLayer.') 

3972 

3973 return elements[1].m 

3974 

3975 def __str__(self): 

3976 return '\n'.join(str(element) for element in self._elements) 

3977 

3978 

3979def read_hyposat_model(fn): 

3980 ''' 

3981 Reader for HYPOSAT earth model files. 

3982 

3983 To be used as producer in :py:meth:`LayeredModel.from_scanlines`. 

3984 

3985 Interface names are translated as follows: ``'MOHO'`` -> ``'moho'``, 

3986 ``'CONR'`` -> ``'conrad'`` 

3987 ''' 

3988 

3989 with open(fn, 'r') as f: 

3990 translate = {'MOHO': 'moho', 'CONR': 'conrad'} 

3991 lname = None 

3992 for iline, line in enumerate(f): 

3993 if iline == 0: 

3994 continue 

3995 

3996 z, vp, vs, name = util.unpack_fixed('f10, f10, f10, a4', line) 

3997 if not name: 

3998 name = None 

3999 material = Material(vp*1000., vs*1000.) 

4000 

4001 tname = translate.get(lname, lname) 

4002 yield z*1000., material, tname 

4003 

4004 lname = name 

4005 

4006 

4007def read_nd_model(fn): 

4008 ''' 

4009 Reader for TauP style '.nd' (named discontinuity) files. 

4010 

4011 To be used as producer in :py:meth:`LayeredModel.from_scanlines`. 

4012 

4013 Interface names are translated as follows: ``'mantle'`` -> ``'moho'``, 

4014 ``'outer-core'`` -> ``'cmb'``, ``'inner-core'`` -> ``'icb'``. 

4015 

4016 The format has been modified to include Burgers materials parameters in 

4017 columns 7 (burger_eta1), 8 (burger_eta2) and 9. eta(3). 

4018 ''' 

4019 with open(fn, 'r') as f: 

4020 for x in read_nd_model_fh(f): 

4021 yield x 

4022 

4023 

4024def read_nd_model_str(s): 

4025 f = StringIO(s) 

4026 for x in read_nd_model_fh(f): 

4027 yield x 

4028 f.close() 

4029 

4030 

4031def read_nd_model_fh(f): 

4032 translate = {'mantle': 'moho', 'outer-core': 'cmb', 'inner-core': 'icb'} 

4033 name = None 

4034 for line in f: 

4035 toks = line.split() 

4036 if len(toks) == 9 or len(toks) == 6 or len(toks) == 4: 

4037 z, vp, vs, rho = [float(x) for x in toks[:4]] 

4038 qp, qs = None, None 

4039 burgers = None 

4040 if len(toks) == 6 or len(toks) == 9: 

4041 qp, qs = [float(x) for x in toks[4:6]] 

4042 if len(toks) == 9: 

4043 burgers = \ 

4044 [float(x) for x in toks[6:]] 

4045 

4046 material = Material( 

4047 vp*1000., vs*1000., rho*1000., qp, qs, 

4048 burgers=burgers) 

4049 

4050 yield z*1000., material, name 

4051 name = None 

4052 elif len(toks) == 1: 

4053 name = translate.get(toks[0], toks[0]) 

4054 

4055 f.close() 

4056 

4057 

4058def from_crust2x2_profile(profile, depthmantle=50000): 

4059 from pyrocko.dataset import crust2x2 

4060 

4061 default_qp_qs = { 

4062 'soft sed.': (50., 50.), 

4063 'hard sed.': (200., 200.), 

4064 'upper crust': (600., 400.), 

4065 } 

4066 

4067 z = 0. 

4068 for i in range(8): 

4069 dz, vp, vs, rho = profile.get_layer(i) 

4070 name = crust2x2.Crust2Profile.layer_names[i] 

4071 if name in default_qp_qs: 

4072 qp, qs = default_qp_qs[name] 

4073 else: 

4074 qp, qs = None, None 

4075 

4076 material = Material(vp, vs, rho, qp, qs) 

4077 iname = None 

4078 if i == 7: 

4079 iname = 'moho' 

4080 if dz != 0.0: 

4081 yield z, material, iname 

4082 if i != 7: 

4083 yield z+dz, material, name 

4084 else: 

4085 yield z+depthmantle, material, name 

4086 

4087 z += dz 

4088 

4089 

4090def write_nd_model_fh(mod, fh): 

4091 def fmt(z, mat): 

4092 rstr = ' '.join( 

4093 util.gform(x, 4) 

4094 for x in ( 

4095 z/1000., 

4096 mat.vp/1000., 

4097 mat.vs/1000., 

4098 mat.rho/1000., 

4099 mat.qp, mat.qs)) 

4100 if not mat._has_default_burgers(): 

4101 rstr += ' '.join( 

4102 util.gform(x, 4) 

4103 for x in ( 

4104 mat.burger_eta1, 

4105 mat.burger_eta2, 

4106 mat.burger_valpha)) 

4107 return rstr.rstrip() + '\n' 

4108 

4109 translate = { 

4110 'moho': 'mantle', 

4111 'cmb': 'outer-core', 

4112 'icb': 'inner-core'} 

4113 

4114 last = None 

4115 for element in mod.elements(): 

4116 if isinstance(element, Interface): 

4117 if element.name is not None: 

4118 n = translate.get(element.name, element.name) 

4119 fh.write('%s\n' % n) 

4120 

4121 elif isinstance(element, Layer): 

4122 if not isinstance(last, Layer): 

4123 fh.write(fmt(element.ztop, element.mtop)) 

4124 

4125 fh.write(fmt(element.zbot, element.mbot)) 

4126 

4127 last = element 

4128 

4129 if not isinstance(last, Layer): 

4130 fh.write(fmt(last.z, last.mbelow)) 

4131 

4132 

4133def write_nd_model_str(mod): 

4134 f = StringIO() 

4135 write_nd_model_fh(mod, f) 

4136 return f.getvalue() 

4137 

4138 

4139def write_nd_model(mod, fn): 

4140 with open(fn, 'w') as f: 

4141 write_nd_model_fh(mod, f) 

4142 

4143 

4144def builtin_models(): 

4145 return sorted([ 

4146 os.path.splitext(os.path.basename(x))[0] 

4147 for x in glob.glob(builtin_model_filename('*'))]) 

4148 

4149 

4150def builtin_model_filename(modelname): 

4151 return util.data_file(os.path.join('earthmodels', modelname+'.nd')) 

4152 

4153 

4154def load_model(fn='ak135-f-continental.m', format='nd', crust2_profile=None): 

4155 ''' 

4156 Load layered earth model from file. 

4157 

4158 :param fn: filename 

4159 :param format: format 

4160 :param crust2_profile: ``(lat, lon)`` or 

4161 :py:class:`pyrocko.crust2x2.Crust2Profile` object, merge model with 

4162 crustal profile. If ``fn`` is forced to be ``None`` only the converted 

4163 CRUST2.0 profile is returned. 

4164 :returns: object of type :py:class:`LayeredModel` 

4165 

4166 The following formats are currently supported: 

4167 

4168 ============== =========================================================== 

4169 format description 

4170 ============== =========================================================== 

4171 ``'nd'`` 'named discontinuity' format used by the TauP programs 

4172 ``'hyposat'`` format used by the HYPOSAT location program 

4173 ============== =========================================================== 

4174 

4175 The naming of interfaces is translated from the file format's native naming 

4176 to Cake's own convention (See :py:func:`read_nd_model` and 

4177 :py:func:`read_hyposat_model` for details). Cake likes the following 

4178 internal names: ``'conrad'``, ``'moho'``, ``'cmb'`` (core-mantle boundary), 

4179 ``'icb'`` (inner core boundary). 

4180 ''' 

4181 

4182 if fn is not None: 

4183 if format == 'nd': 

4184 if not os.path.exists(fn) and fn in builtin_models(): 

4185 fn = builtin_model_filename(fn) 

4186 reader = read_nd_model(fn) 

4187 elif format == 'hyposat': 

4188 reader = read_hyposat_model(fn) 

4189 else: 

4190 assert False, 'unsupported model format' 

4191 

4192 mod = LayeredModel.from_scanlines(reader) 

4193 if crust2_profile is not None: 

4194 return mod.replaced_crust(crust2_profile) 

4195 

4196 return mod 

4197 

4198 else: 

4199 assert crust2_profile is not None 

4200 profile = anything_to_crust2_profile(crust2_profile) 

4201 return LayeredModel.from_scanlines( 

4202 from_crust2x2_profile(profile)) 

4203 

4204 

4205def castagna_vs_to_vp(vs): 

4206 ''' 

4207 Calculate vp from vs using Castagna's relation. 

4208 

4209 Castagna's relation (the mudrock line) is an empirical relation for vp/vs 

4210 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al., 

4211 1985] 

4212 

4213 vp = 1.16 * vs + 1360 [m/s] 

4214 

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

4216 :returns: P-wave velocity [m/s] 

4217 ''' 

4218 

4219 return vs*1.16 + 1360.0 

4220 

4221 

4222def castagna_vp_to_vs(vp): 

4223 ''' 

4224 Calculate vp from vs using Castagna's relation. 

4225 

4226 Castagna's relation (the mudrock line) is an empirical relation for vp/vs 

4227 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al., 

4228 1985] 

4229 

4230 vp = 1.16 * vs + 1360 [m/s] 

4231 

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

4233 :returns: S-wave velocity [m/s] 

4234 ''' 

4235 

4236 return (vp - 1360.0) / 1.16 

4237 

4238 

4239def evenize(x, y, minsize=10): 

4240 if x.size < minsize: 

4241 return x 

4242 ry = (y.max()-y.min()) 

4243 if ry == 0: 

4244 return x 

4245 dx = (x[1:] - x[:-1])/(x.max()-x.min()) 

4246 dy = (y[1:] + y[:-1])/ry 

4247 

4248 s = num.zeros(x.size) 

4249 s[1:] = num.cumsum(num.sqrt(dy**2 + dx**2)) 

4250 s2 = num.linspace(0, s[-1], x.size) 

4251 x2 = num.interp(s2, s, x) 

4252 x2[0] = x[0] 

4253 x2[-1] = x[-1] 

4254 return x2 

4255 

4256 

4257def filled(v, *args, **kwargs): 

4258 ''' 

4259 Create NumPy array filled with given value. 

4260 

4261 This works like :py:func:`numpy.ones` but initializes the array with ``v`` 

4262 instead of ones. 

4263 ''' 

4264 x = num.empty(*args, **kwargs) 

4265 x.fill(v) 

4266 return x 

4267 

4268 

4269def next_or_none(i): 

4270 try: 

4271 return next(i) 

4272 except StopIteration: 

4273 return None 

4274 

4275 

4276def reci_or_none(x): 

4277 try: 

4278 return 1./x 

4279 except ZeroDivisionError: 

4280 return None 

4281 

4282 

4283def mult_or_none(a, b): 

4284 if a is None or b is None: 

4285 return None 

4286 return a*b 

4287 

4288 

4289def min_not_none(a, b): 

4290 if a is None: 

4291 return b 

4292 if b is None: 

4293 return a 

4294 return min(a, b) 

4295 

4296 

4297def xytups(xx, yy): 

4298 d = [] 

4299 for x, y in zip(xx, yy): 

4300 if num.isfinite(y): 

4301 d.append((x, y)) 

4302 return d 

4303 

4304 

4305def interp(x, xp, fp, monoton): 

4306 if monoton == 1: 

4307 return xytups( 

4308 x, num.interp(x, xp, fp, left=num.nan, right=num.nan)) 

4309 elif monoton == -1: 

4310 return xytups( 

4311 x, num.interp(x, xp[::-1], fp[::-1], left=num.nan, right=num.nan)) 

4312 else: 

4313 fs = [] 

4314 for xv in x: 

4315 indices = num.where(num.logical_or( 

4316 num.logical_and(xp[:-1] >= xv, xv > xp[1:]), 

4317 num.logical_and(xp[:-1] <= xv, xv < xp[1:])))[0] 

4318 

4319 for i in indices: 

4320 xr = (xv - xp[i])/(xp[i+1]-xp[i]) 

4321 fv = xr*fp[i] + (1.-xr)*fp[i+1] 

4322 fs.append((xv, fv)) 

4323 

4324 return fs 

4325 

4326 

4327def float_or_none(x): 

4328 if x is not None: 

4329 return float(x) 

4330 

4331 

4332def parstore_float(thelocals, obj, *args): 

4333 for k, v in thelocals.items(): 

4334 if k != 'self' and (not args or k in args): 

4335 setattr(obj, k, float_or_none(v))