1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Classical seismic ray theory for layered earth models (*layer cake* models).
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.
14The main classes defined in this module are:
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`.
24 * :py:class:`HomogeneousLayer` - A homogeneous :py:class:`Layer`.
25 * :py:class:`GradientLayer` - A gradient :py:class:`Layer`.
27* :py:class:`Discontinuity` - A discontinuity in a :py:class:`LayeredModel`.
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`.
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`.
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'''
46from __future__ import absolute_import
47from functools import reduce
49import os
50import logging
51import copy
52import math
53import cmath
54import operator
55try:
56 from StringIO import StringIO
57except ImportError:
58 from io import StringIO
60import glob
61import numpy as num
62from scipy.optimize import bisect, brentq
64from . import util, config
66try:
67 newstr = unicode
68except NameError:
69 newstr = str
71logger = logging.getLogger('cake')
73ZEPS = 0.01
74P = 1
75S = 2
76DOWN = 4
77UP = -4
79DEFAULT_BURGERS = (0., 0., 1.)
81earthradius = config.config().earthradius
83r2d = 180./math.pi
84d2r = 1./r2d
85km = 1000.
86d2m = d2r*earthradius
87m2d = 1./d2m
88sprad2spm = 1.0/(r2d*d2m)
89sprad2spkm = 1.0/(r2d*d2m/km)
90spm2sprad = 1.0/sprad2spm
91spkm2sprad = 1.0/sprad2spkm
94class CakeError(Exception):
95 pass
98class InvalidArguments(CakeError):
99 pass
102class Material(object):
103 '''
104 Isotropic elastic material.
106 :param vp: P-wave velocity [m/s]
107 :param vs: S-wave velocity [m/s]
108 :param rho: density [kg/m^3]
109 :param qp: P-wave attenuation Qp
110 :param qs: S-wave attenuation Qs
111 :param poisson: Poisson ratio
112 :param lame: tuple with Lame parameter `lambda` and `shear modulus` [Pa]
113 :param qk: bulk attenuation Qk
114 :param qmu: shear attenuation Qmu
116 :param burgers: Burgers rheology paramerters as `tuple`.
117 `transient viscosity` [Pa], <= 0 means infinite value,
118 `steady-state viscosity` [Pa] and `alpha`, the ratio between the
119 effective and unreleaxed shear modulus, mu1/(mu1 + mu2).
120 :type burgers: tuple
122 If no velocities and no lame parameters are given, standard crustal values
123 of vp = 5800 m/s and vs = 3200 m/s are used. If no Q values are given,
124 standard crustal values of qp = 1456 and qs = 600 are used. If no Burgers
125 material parameters are given, transient and steady-state viscosities are
126 0 and alpha=1.
128 Everything is in SI units (m/s, Pa, kg/m^3) unless explicitly stated.
130 The main material properties are considered independant and are accessible
131 as attributes (it is allowed to assign to these):
133 .. py:attribute:: vp, vs, rho, qp, qs
135 Other material properties are considered dependant and can be queried by
136 instance methods.
137 '''
139 def __init__(
140 self, vp=None, vs=None, rho=2600., qp=None, qs=None, poisson=None,
141 lame=None, qk=None, qmu=None, burgers=None):
143 parstore_float(locals(), self, 'vp', 'vs', 'rho', 'qp', 'qs')
145 if vp is not None and vs is not None:
146 if poisson is not None or lame is not None:
147 raise InvalidArguments(
148 'If vp and vs are given, poisson ratio and lame paramters '
149 'should not be given.')
151 elif vp is None and vs is None and lame is None:
152 self.vp = 5800.
153 if poisson is None:
154 poisson = 0.25
155 self.vs = self.vp / math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson))
157 elif vp is None and vs is None and lame is not None:
158 if poisson is not None:
159 raise InvalidArguments(
160 'Poisson ratio should not be given, when lame parameters '
161 'are given.')
163 lam, mu = float(lame[0]), float(lame[1])
164 self.vp = math.sqrt((lam + 2.0*mu)/rho)
165 self.vs = math.sqrt(mu/rho)
167 elif vp is not None and vs is None:
168 if poisson is None:
169 poisson = 0.25
171 if lame is not None:
172 raise InvalidArguments(
173 'If vp is given, Lame parameters should not be given.')
175 poisson = float(poisson)
176 self.vs = vp / math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson))
178 elif vp is None and vs is not None:
179 if poisson is None:
180 poisson = 0.25
181 if lame is not None:
182 raise InvalidArguments(
183 'If vs is given, Lame parameters should not be given.')
185 poisson = float(poisson)
186 self.vp = vs * math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson))
188 else:
189 raise InvalidArguments(
190 'Invalid combination of input parameters in material '
191 'definition.')
193 if qp is not None or qs is not None:
194 if not (qk is None and qmu is None):
195 raise InvalidArguments(
196 'If qp or qs are given, qk and qmu should not be given.')
198 if qp is None:
199 if self.vs != 0.0:
200 s = (4.0/3.0)*(self.vs/self.vp)**2
201 self.qp = self.qs / s
202 else:
203 self.qp = 1456.
205 if qs is None:
206 if self.vs != 0.0:
207 s = (4.0/3.0)*(self.vs/self.vp)**2
208 self.qs = self.qp * s
209 else:
210 self.vs = 600.
212 elif qp is None and qs is None and qk is None and qmu is None:
213 if self.vs == 0.:
214 self.qs = 0.
215 self.qp = 5782e4
216 else:
217 self.qs = 600.
218 s = (4.0/3.0)*(self.vs/self.vp)**2
219 self.qp = self.qs/s
221 elif qp is None and qs is None and qk is not None and qmu is not None:
222 s = (4.0/3.0)*(self.vs/self.vp)**2
223 if qmu == 0. and self.vs == 0.:
224 self.qp = qk
225 else:
226 if num.isinf(qk):
227 self.qp = qmu/s
228 else:
229 self.qp = 1.0 / (s/qmu + (1.0-s)/qk)
230 self.qs = qmu
231 else:
232 raise InvalidArguments(
233 'Invalid combination of input parameters in material '
234 'definition.')
236 if burgers is None:
237 burgers = DEFAULT_BURGERS
239 self.burger_eta1 = burgers[0]
240 self.burger_eta2 = burgers[1]
241 self.burger_valpha = burgers[2]
243 def astuple(self):
244 '''
245 Get independant material properties as a tuple.
247 Returns a tuple with ``(vp, vs, rho, qp, qs)``.
248 '''
249 return self.vp, self.vs, self.rho, self.qp, self.qs
251 def __eq__(self, other):
252 return self.astuple() == other.astuple()
254 def lame(self):
255 '''
256 Get Lame's parameter lambda and shear modulus.
257 '''
258 mu = self.vs**2 * self.rho
259 lam = self.vp**2 * self.rho - 2.0*mu
260 return lam, mu
262 def lame_lambda(self):
263 '''
264 Get Lame's parameter lambda.
266 Returned units are [Pa].
267 '''
268 lam, _ = self.lame()
269 return lam
271 def shear_modulus(self):
272 '''
273 Get shear modulus.
275 Returned units are [Pa].
276 '''
277 return self.vs**2 * self.rho
279 def poisson(self):
280 '''
281 Get Poisson's ratio.
282 '''
283 lam, mu = self.lame()
284 return lam / (2.0*(lam+mu))
286 def bulk(self):
287 '''
288 Get bulk modulus.
289 '''
290 lam, mu = self.lame()
291 return lam + 2.0*mu/3.0
293 def youngs(self):
294 '''
295 Get Young's modulus.
296 '''
297 lam, mu = self.lame()
298 return mu * (3.0*lam + 2.0*mu) / (lam+mu)
300 def vp_vs_ratio(self):
301 '''
302 Get vp/vs ratio.
303 '''
304 return self.vp/self.vs
306 def qmu(self):
307 '''
308 Get shear attenuation coefficient Qmu.
309 '''
310 return self.qs
312 def qk(self):
313 '''
314 Get bulk attenuation coefficient Qk.
315 '''
316 if self.vs == 0. and self.qs == 0.:
317 return self.qp
318 else:
319 s = (4.0/3.0)*(self.vs/self.vp)**2
320 denom = (1/self.qp - s/self.qs)
321 if denom <= 0.0:
322 return num.inf
323 else:
324 return (1.-s)/(1.0/self.qp - s/self.qs)
326 def burgers(self):
327 '''
328 Get Burger parameters.
329 '''
330 return self.burger_eta1, self.burger_eta2, self.burger_valpha
332 def _rayleigh_equation(self, cr):
333 cr_a = (cr/self.vp)**2
334 cr_b = (cr/self.vs)**2
335 if cr_a > 1.0 or cr_b > 1.0:
336 return None
338 return (2.0-cr_b)**2 - 4.0 * math.sqrt(1.0-cr_a) * math.sqrt(1.0-cr_b)
340 def rayleigh(self):
341 '''
342 Get Rayleigh velocity assuming a homogenous halfspace.
344 Returned units are [m/s].
345 '''
346 return bisect(self._rayleigh_equation, 0.001*self.vs, self.vs)
348 def _has_default_burgers(self):
349 if self.burger_eta1 == DEFAULT_BURGERS[0] and \
350 self.burger_eta2 == DEFAULT_BURGERS[1] and \
351 self.burger_valpha == DEFAULT_BURGERS[2]:
352 return True
353 return False
355 def describe(self):
356 '''
357 Get a readable listing of the material properties.
358 '''
359 template = '''
360P wave velocity [km/s] : %12g
361S wave velocity [km/s] : %12g
362P/S wave vel. ratio : %12g
363Lame lambda [GPa] : %12g
364Lame shear modulus [GPa] : %12g
365Poisson ratio : %12g
366Bulk modulus [GPa] : %12g
367Young's modulus [GPa] : %12g
368Rayleigh wave vel. [km/s] : %12g
369Density [g/cm**3] : %12g
370Qp P-wave attenuation : %12g
371Qs S-wave attenuation (Qmu) : %12g
372Qk bulk attenuation : %12g
373transient viscos., eta1 [GPa] : %12g
374st.-state viscos., eta2 [GPa] : %12g
375relaxation: valpha : %12g
376'''.strip()
378 return template % (
379 self.vp/km,
380 self.vs/km,
381 self.vp/self.vs,
382 self.lame_lambda()*1e-9,
383 self.shear_modulus()*1e-9,
384 self.poisson(),
385 self.bulk()*1e-9,
386 self.youngs()*1e-9,
387 self.rayleigh()/km,
388 self.rho/km,
389 self.qp,
390 self.qs,
391 self.qk(),
392 self.burger_eta1*1e-9,
393 self.burger_eta2*1e-9,
394 self.burger_valpha)
396 def __str__(self):
397 vp, vs, rho, qp, qs = self.astuple()
398 return '%10g km/s %10g km/s %10g g/cm^3 %10g %10g' % (
399 vp/km, vs/km, rho/km, qp, qs)
401 def __repr__(self):
402 return 'Material(vp=%s, vs=%s, rho=%s, qp=%s, qs=%s)' % \
403 tuple(repr(x) for x in (
404 self.vp, self.vs, self.rho, self.qp, self.qs))
407class Leg(object):
408 '''
409 Represents a continuous piece of wave propagation in a phase definition.
411 **Attributes:**
413 To be considered as read-only.
415 .. py:attribute:: departure
417 One of the constants :py:const:`UP` or :py:const:`DOWN` indicating
418 upward or downward departure.
420 .. py:attribute:: mode
422 One of the constants :py:const:`P` or :py:const:`S`, indicating the
423 propagation mode.
425 .. py:attribute:: depthmin
427 ``None``, a number (a depth in [m]) or a string (an interface name),
428 minimum depth.
430 .. py:attribute:: depthmax
432 ``None``, a number (a depth in [m]) or a string (an interface name),
433 maximum depth.
435 '''
437 def __init__(self, departure=None, mode=None):
438 self.departure = departure
439 self.mode = mode
440 self.depthmin = None
441 self.depthmax = None
443 def set_depthmin(self, depthmin):
444 self.depthmin = depthmin
446 def set_depthmax(self, depthmax):
447 self.depthmax = depthmax
449 def __str__(self):
450 def sd(d):
451 if isinstance(d, float):
452 return '%g km' % (d/km)
453 else:
454 return 'interface %s' % d
456 s = '%s mode propagation, departing %s' % (
457 smode(self.mode).upper(), {
458 UP: 'upward', DOWN: 'downward'}[self.departure])
460 sc = []
461 if self.depthmax is not None:
462 sc.append('deeper than %s' % sd(self.depthmax))
463 if self.depthmin is not None:
464 sc.append('shallower than %s' % sd(self.depthmin))
466 if sc:
467 s = s + ' (may not propagate %s)' % ' or '.join(sc)
469 return s
472class InvalidKneeDef(CakeError):
473 pass
476class Knee(object):
477 '''
478 Represents a change in wave propagation within a :py:class:`PhaseDef`.
480 **Attributes:**
482 To be considered as read-only.
484 .. py:attribute:: depth
486 Depth at which the conversion/reflection should happen. this can be
487 a string or a number.
489 .. py:attribute:: direction
491 One of the constants :py:const:`UP` or :py:const:`DOWN` to indicate
492 the incoming direction.
494 .. py:attribute:: in_mode
496 One of the constants :py:const:`P` or :py:const:`S` to indicate the
497 type of mode of the incoming wave.
499 .. py:attribute:: out_mode
501 One of the constants :py:const:`P` or :py:const:`S` to indicate the
502 type of mode of the outgoing wave.
504 .. py:attribute:: conversion
506 Boolean, whether there is a mode conversion involved.
508 .. py:attribute:: reflection
510 Boolean, whether there is a reflection involved.
512 .. py:attribute:: headwave
514 Boolean, whether there is headwave propagation involved.
516 '''
518 defaults = dict(
519 depth='surface',
520 direction=UP,
521 conversion=True,
522 reflection=False,
523 headwave=False,
524 in_setup_state=True)
526 defaults_surface = dict(
527 depth='surface',
528 direction=UP,
529 conversion=False,
530 reflection=True,
531 headwave=False,
532 in_setup_state=True)
534 def __init__(self, *args):
535 if args:
536 (self.depth, self.direction, self.reflection, self.in_mode,
537 self.out_mode) = args
539 self.conversion = self.in_mode != self.out_mode
540 self.in_setup_state = False
542 def default(self, k):
543 depth = self.__dict__.get('depth', 'surface')
544 if depth == 'surface':
545 return Knee.defaults_surface[k]
546 else:
547 return Knee.defaults[k]
549 def __setattr__(self, k, v):
550 if self.in_setup_state and k in self.__dict__:
551 raise InvalidKneeDef('%s has already been set' % k)
552 else:
553 self.__dict__[k] = v
555 def __getattr__(self, k):
556 if k.startswith('__'):
557 raise AttributeError(k)
559 if k not in self.__dict__:
560 return self.default(k)
562 def set_modes(self, in_leg, out_leg):
564 if out_leg.departure == UP and (
565 (self.direction == UP) == self.reflection):
567 raise InvalidKneeDef(
568 'cannot enter %s from %s and emit ray upwards' % (
569 ['conversion', 'reflection'][self.reflection],
570 {UP: 'below', DOWN: 'above'}[self.direction]))
572 if out_leg.departure == DOWN and (
573 (self.direction == DOWN) == self.reflection):
575 raise InvalidKneeDef(
576 'cannot enter %s from %s and emit ray downwards' % (
577 ['conversion', 'reflection'][self.reflection],
578 {UP: 'below', DOWN: 'above'}[self.direction]))
580 self.in_mode = in_leg.mode
581 self.out_mode = out_leg.mode
583 def at_surface(self):
584 return self.depth == 'surface'
586 def matches(self, discontinuity, mode, direction):
587 '''
588 Check whether it is relevant to a given combination of interface,
589 propagation mode, and direction.
590 '''
592 if isinstance(self.depth, float):
593 if abs(self.depth - discontinuity.z) > ZEPS:
594 return False
595 else:
596 if discontinuity.name != self.depth:
597 return False
599 return self.direction == direction and self.in_mode == mode
601 def out_direction(self):
602 '''
603 Get outgoing direction.
605 Returns one of the constants :py:const:`UP` or :py:const:`DOWN`.
606 '''
608 if self.reflection:
609 return - self.direction
610 else:
611 return self.direction
613 def __str__(self):
614 x = []
615 if self.reflection:
616 if self.at_surface():
617 x.append('surface')
618 else:
619 if not self.headwave:
620 if self.direction == UP:
621 x.append('underside')
622 else:
623 x.append('upperside')
625 if self.headwave:
626 x.append('headwave propagation along')
627 elif self.reflection and self.conversion:
628 x.append('reflection with conversion from %s to %s' % (
629 smode(self.in_mode).upper(), smode(self.out_mode).upper()))
630 if not self.at_surface():
631 x.append('at')
632 elif self.reflection:
633 x.append('reflection')
634 if not self.at_surface():
635 x.append('at')
636 elif self.conversion:
637 x.append('conversion from %s to %s at' % (
638 smode(self.in_mode).upper(), smode(self.out_mode).upper()))
639 else:
640 x.append('passing through')
642 if isinstance(self.depth, float):
643 x.append('interface in %g km depth' % (self.depth/1000.))
644 else:
645 if not self.at_surface():
646 x.append('%s' % self.depth)
648 if not self.reflection:
649 if self.direction == UP:
650 x.append('on upgoing path')
651 else:
652 x.append('on downgoing path')
654 return ' '.join(x)
657class Head(Knee):
658 def __init__(self, *args):
659 if args:
660 z, in_direction, mode = args
661 Knee.__init__(self, z, in_direction, True, mode, mode)
662 else:
663 Knee.__init__(self)
665 def __str__(self):
666 x = ['propagation as headwave']
667 if isinstance(self.depth, float):
668 x.append('at interface in %g km depth' % (self.depth/1000.))
669 else:
670 x.append('at %s' % self.depth)
672 return ' '.join(x)
675class UnknownClassicPhase(CakeError):
676 def __init__(self, phasename):
677 self.phasename = phasename
679 def __str__(self):
680 return 'Unknown classic phase name: %s' % self.phasename
683class PhaseDefParseError(CakeError):
684 '''
685 Exception raised when an error occures during parsing of a phase
686 definition string.
687 '''
689 def __init__(self, definition, position, exception):
690 self.definition = definition
691 self.position = position
692 self.exception = exception
694 def __str__(self):
695 return 'Invalid phase definition: "%s" (at character %i: %s)' % (
696 self.definition, self.position+1, str(self.exception))
699class PhaseDef(object):
701 '''
702 Definition of a seismic phase arrival, based on ray propagation path.
704 :param definition: string representation of the phase in Cake's phase
705 syntax
707 Seismic phases are conventionally named e.g. P, Pn, PP, PcP, etc. In Cake,
708 a slightly different terminology is adapted, which allows to specify
709 arbitrary conversion/reflection histories for seismic ray paths. The
710 conventions used here are inspired by those used in the TauP toolkit, but
711 are not completely compatible with those.
713 The definition of a seismic ray propagation path in Cake's phase syntax is
714 a string consisting of an alternating sequence of *legs* and *knees*.
716 A *leg* represents seismic wave propagation without any conversions,
717 encountering only super-critical reflections. Legs are denoted by ``P``,
718 ``p``, ``S``, or ``s``. The capital letters are used when the take-off of
719 the *leg* is in downward direction, while the lower case letters indicate a
720 take-off in upward direction.
722 A *knee* is an interaction with an interface. It can be a mode conversion,
723 a reflection, or propagation as a headwave or diffracted wave.
725 * conversion is simply denoted as: ``(INTERFACE)`` or ``DEPTH``
726 * upperside reflection: ``v(INTERFACE)`` or ``vDEPTH``
727 * underside reflection: ``^(INTERFACE)`` or ``^DEPTH``
728 * normal kind headwave or diffracted wave: ``v_(INTERFACE)`` or
729 ``v_DEPTH``
731 The interface may be given by name or by depth: INTERFACE is the name of an
732 interface defined in the model, DEPTH is the depth of an interface in
733 [km] (the interface closest to that depth is chosen). If two legs appear
734 consecutively without an explicit *knee*, surface interaction is assumed.
736 The phase definition may end with a backslash ``\\``, to indicate that the
737 ray should arrive at the receiver from above instead of from below. It is
738 possible to restrict the maximum and minimum depth of a *leg* by appending
739 ``<(INTERFACE)`` or ``<DEPTH`` or ``>(INTERFACE)`` or ``>DEPTH`` after the
740 leg character, respectively.
742 **Examples:**
744 * ``P`` - like the classical P, but includes PKP, PKIKP, Pg
745 * ``P<(moho)`` - like classical Pg, but must leave source downwards
746 * ``pP`` - leaves source upward, reflects at surface, then travels as P
747 * ``P(moho)s`` - conversion from P to S at the Moho on upgoing path
748 * ``P(moho)S`` - conversion from P to S at the Moho on downgoing path
749 * ``Pv12p`` - P with reflection at 12 km deep interface (or the
750 interface closest to that)
751 * ``Pv_(moho)p`` - classical Pn
752 * ``Pv_(cmb)p`` - classical Pdiff
753 * ``P^(conrad)P`` - underside reflection of P at the Conrad
754 discontinuity
756 **Usage:**
758 >>> from pyrocko.cake import PhaseDef
759 # must escape the backslash
760 >>> my_crazy_phase = PhaseDef('pPv(moho)sP\\\\')
761 >>> print my_crazy_phase
762 Phase definition "pPv(moho)sP\":
763 - P mode propagation, departing upward
764 - surface reflection
765 - P mode propagation, departing downward
766 - upperside reflection with conversion from P to S at moho
767 - S mode propagation, departing upward
768 - surface reflection with conversion from S to P
769 - P mode propagation, departing downward
770 - arriving at target from above
772 .. note::
774 (1) These conventions might be extended in a way to allow to fix wave
775 propagation to SH mode, possibly by specifying SH, or a single
776 character (e.g. H) instead of S. This would be benificial for the
777 selection of conversion and reflection coefficients, which
778 currently only deal with the P-SV case.
779 '''
781 allowed_characters_pattern = r'[0-9a-zA-Z_()<>^v\\.]+'
782 allowed_characters_pattern_classic = r'[a-zA-Z0-9]+'
784 @staticmethod
785 def classic_definitions():
786 defs = {}
787 # PmP, PmS, PcP, PcS, SmP, ...
788 for r in 'mc':
789 for a, b in 'PP PS SS SP'.split():
790 defs[a+r+b] = [
791 '%sv(%s)%s' % (a, {'m': 'moho', 'c': 'cmb'}[r], b.lower())]
793 # Pg, P, S, Sg
794 for a in 'PS':
795 defs[a+'g'] = ['%s<(moho)' % x for x in (a, a.lower())]
796 defs[a] = ['%s<(cmb)(moho)%s' % (x, x.lower()) for x in (
797 a, a.lower())]
799 defs[a.lower()] = [a.lower()]
801 for a, b in 'PP PS SS SP'.split():
802 defs[a+'K'+b] = ['%s(cmb)P<(icb)(cmb)%s' % (a, b.lower())]
803 defs[a+'KIK'+b] = ['%s(cmb)P(icb)P(icb)p(cmb)%s' % (a, b.lower())]
804 defs[a+'KJK'+b] = ['%s(cmb)P(icb)S(icb)p(cmb)%s' % (a, b.lower())]
805 defs[a+'KiK'+b] = ['%s(cmb)Pv(icb)p(cmb)%s' % (a, b.lower())]
807 # PP, SS, PS, SP, PPP, ...
808 for a in 'PS':
809 for b in 'PS':
810 for c in 'PS':
811 defs[a+b+c] = [''.join(defs[x][0] for x in a+b+c)]
813 defs[a+b] = [''.join(defs[x][0] for x in a+b)]
815 # Pc, Pdiff, Sc, ...
816 for x in 'PS':
817 defs[x+'c'] = defs[x+'diff'] = [x+'v_(cmb)'+x.lower()]
818 defs[x+'n'] = [x+'v_(moho)'+x.lower()]
820 # depth phases
821 for k in list(defs.keys()):
822 if k not in 'ps':
823 for x in 'ps':
824 defs[x+k] = [x + defs[k][0]]
826 return defs
828 @staticmethod
829 def classic(phasename):
830 '''
831 Get phase definitions based on classic phase name.
833 :param phasename: classic name of a phase
834 :returns: list of PhaseDef objects
836 This returns a list of PhaseDef objects, because some classic phases
837 (like e.g. Pg) can only be represented by two Cake style PhaseDef
838 objects (one with downgoing and one with upgoing first leg).
839 '''
841 defs = PhaseDef.classic_definitions()
842 if phasename not in defs:
843 raise UnknownClassicPhase(phasename)
845 return [PhaseDef(d, classicname=phasename) for d in defs[phasename]]
847 def __init__(self, definition=None, classicname=None):
849 state = 0
850 sdepth = ''
851 sinterface = ''
852 depthmax = depthmin = None
853 depthlim = None
854 depthlimtype = None
855 sdepthlim = ''
856 events = []
857 direction_stop = UP
858 need_leg = True
859 ic = 0
860 if definition is not None:
861 knee = Knee()
862 try:
863 for ic, c in enumerate(definition):
865 if state in (0, 1):
867 if c in '0123456789.':
868 need_leg = True
869 state = 1
870 sdepth += c
871 continue
873 elif state == 1:
874 knee.depth = float(sdepth)*1000.
875 state = 0
877 if state == 2:
878 if c == ')':
879 knee.depth = sinterface
880 state = 0
881 else:
882 sinterface += c
884 continue
886 if state in (3, 4):
888 if state == 3:
889 if c in '0123456789.':
890 sdepthlim += c
891 continue
892 elif c == '(':
893 state = 4
894 continue
895 else:
896 depthlim = float(sdepthlim)*1000.
897 if depthlimtype == '<':
898 depthmax = depthlim
899 else:
900 depthmin = depthlim
901 state = 0
903 elif state == 4:
904 if c == ')':
905 depthlim = sdepthlim
906 if depthlimtype == '<':
907 depthmax = depthlim
908 else:
909 depthmin = depthlim
910 state = 0
911 continue
912 else:
913 sdepthlim += c
914 continue
916 if state == 0:
918 if c == '(':
919 need_leg = True
920 state = 2
921 continue
923 elif c in '<>':
924 state = 3
925 depthlim = None
926 sdepthlim = ''
927 depthlimtype = c
928 continue
930 elif c in 'psPS':
931 leg = Leg()
932 if c in 'ps':
933 leg.departure = UP
934 else:
935 leg.departure = DOWN
936 leg.mode = imode(c)
938 if events:
939 in_leg = events[-1]
940 if depthmin is not None:
941 in_leg.set_depthmin(depthmin)
942 depthmin = None
943 if depthmax is not None:
944 in_leg.set_depthmax(depthmax)
945 depthmax = None
947 if in_leg.mode != leg.mode:
948 knee.conversion = True
949 else:
950 knee.conversion = False
952 if not knee.reflection:
953 if c in 'ps':
954 knee.direction = UP
955 else:
956 knee.direction = DOWN
958 knee.set_modes(in_leg, leg)
959 knee.in_setup_state = False
960 events.append(knee)
961 knee = Knee()
962 sdepth = ''
963 sinterface = ''
965 events.append(leg)
966 need_leg = False
967 continue
969 elif c == '^':
970 need_leg = True
971 knee.direction = UP
972 knee.reflection = True
973 continue
975 elif c == 'v':
976 need_leg = True
977 knee.direction = DOWN
978 knee.reflection = True
979 continue
981 elif c == '_':
982 need_leg = True
983 knee.headwave = True
984 continue
986 elif c == '\\':
987 direction_stop = DOWN
988 continue
990 else:
991 raise PhaseDefParseError(
992 definition, ic, 'invalid character: "%s"' % c)
994 if state == 3:
995 depthlim = float(sdepthlim)*1000.
996 if depthlimtype == '<':
997 depthmax = depthlim
998 else:
999 depthmin = depthlim
1000 state = 0
1002 except (ValueError, InvalidKneeDef) as e:
1003 raise PhaseDefParseError(definition, ic, e)
1005 if state != 0 or need_leg:
1006 raise PhaseDefParseError(
1007 definition, ic, 'unfinished expression')
1009 if events and depthmin is not None:
1010 events[-1].set_depthmin(depthmin)
1011 if events and depthmax is not None:
1012 events[-1].set_depthmax(depthmax)
1014 self._definition = definition
1015 self._classicname = classicname
1016 self._events = events
1017 self._direction_stop = direction_stop
1019 def __iter__(self):
1020 for ev in self._events:
1021 yield ev
1023 def append(self, ev):
1024 self._events.append(ev)
1026 def first_leg(self):
1027 '''
1028 Get the first leg in phase definition.
1029 '''
1030 return self._events[0]
1032 def last_leg(self):
1033 '''
1034 Get the last leg in phase definition.
1035 '''
1036 return self._events[-1]
1038 def legs(self):
1039 '''
1040 Iterate over the continuous pieces of wave propagation (legs) defined
1041 within this phase definition.
1042 '''
1044 return (leg for leg in self if isinstance(leg, Leg))
1046 def knees(self):
1047 '''
1048 Iterate over conversions and reflections (knees) defined within this
1049 phase definition.
1050 '''
1051 return (knee for knee in self if isinstance(knee, Knee))
1053 def definition(self):
1054 '''
1055 Get original definition of the phase.
1056 '''
1057 return self._definition
1059 def given_name(self):
1060 '''
1061 Get entered classic name if any, or original definition of the phase.
1062 '''
1064 if self._classicname:
1065 return self._classicname
1066 else:
1067 return self._definition
1069 def direction_start(self):
1070 return self.first_leg().departure
1072 def direction_stop(self):
1073 return self._direction_stop
1075 def headwave_knee(self):
1076 for el in self:
1077 if type(el) == Knee and el.headwave:
1078 return el
1079 return None
1081 def used_repr(self):
1082 '''
1083 Translate into textual representation (cake phase syntax).
1084 '''
1085 def strdepth(x):
1086 if isinstance(x, float):
1087 return '%g' % (x/1000.)
1088 else:
1089 return '(%s)' % x
1091 x = []
1092 for el in self:
1093 if type(el) == Leg:
1094 if el.departure == UP:
1095 x.append(smode(el.mode).lower())
1096 else:
1097 x.append(smode(el.mode).upper())
1099 if el.depthmax is not None:
1100 x.append('<'+strdepth(el.depthmax))
1102 if el.depthmin is not None:
1103 x.append('>'+strdepth(el.depthmin))
1105 elif type(el) == Knee:
1106 if el.reflection and not el.at_surface():
1107 if el.direction == DOWN:
1108 x.append('v')
1109 else:
1110 x.append('^')
1111 if el.headwave:
1112 x.append('_')
1113 if not el.at_surface():
1114 x.append(strdepth(el.depth))
1116 elif type(el) == Head:
1117 x.append('_')
1118 x.append(strdepth(el.depth))
1120 if self._direction_stop == DOWN:
1121 x.append('\\')
1123 return ''.join(x)
1125 def __repr__(self):
1126 if self._definition is not None:
1127 return "PhaseDef('%s')" % self._definition
1128 else:
1129 return "PhaseDef('%s')" % self.used_repr()
1131 def __str__(self):
1132 orig = ''
1133 used = self.used_repr()
1134 if self._definition != used:
1135 orig = ' (entered as "%s")' % self._definition
1137 sarrive = '\n - arriving at target from %s' % ('below', 'above')[
1138 self._direction_stop == DOWN]
1140 return 'Phase definition "%s"%s:\n - ' % (used, orig) + \
1141 '\n - '.join(str(ev) for ev in self) + sarrive
1143 def copy(self):
1144 '''
1145 Get a deep copy of it.
1146 '''
1147 return copy.deepcopy(self)
1150def to_phase_defs(phases):
1151 if isinstance(phases, (str, newstr, PhaseDef)):
1152 phases = [phases]
1154 phases_out = []
1155 for phase in phases:
1156 if isinstance(phase, (str, newstr)):
1157 phases_out.extend(PhaseDef(x.strip()) for x in phase.split(','))
1158 elif isinstance(phase, PhaseDef):
1159 phases_out.append(phase)
1160 else:
1161 raise PhaseDefParseError('invalid phase definition')
1163 return phases_out
1166def csswap(x):
1167 return cmath.sqrt(1.-x**2)
1170def psv_surface_ind(in_mode, out_mode):
1171 '''
1172 Get indices to select the appropriate element from scatter matrix for free
1173 surface.
1174 '''
1176 return (int(in_mode == S), int(out_mode == S))
1179def psv_surface(material, p, energy=False):
1180 '''
1181 Scatter matrix for free surface reflection/conversions.
1183 :param material: material, object of type :py:class:`Material`
1184 :param p: flat ray parameter [s/m]
1185 :param energy: bool, when ``True`` energy normalized coefficients are
1186 returned
1187 :returns: Scatter matrix
1189 The scatter matrix is ordered as follows::
1191 [[PP, PS],
1192 [SP, SS]]
1194 The formulas given in Aki & Richards are used.
1195 '''
1197 vp, vs, rho = material.vp, material.vs, material.rho
1198 sinphi = p * vp
1199 sinlam = p * vs
1200 cosphi = csswap(sinphi)
1201 coslam = csswap(sinlam)
1203 if vs == 0.0:
1204 scatter = num.array([[-1.0, 0.0], [0.0, 1.0]])
1206 else:
1207 vsp_term = (1.0/vs**2 - 2.0*p**2)
1208 pcc_term = 4.0 * p**2 * cosphi/vp * coslam/vs
1209 denom = vsp_term**2 + pcc_term
1211 scatter = num.array([
1212 [- vsp_term**2 + pcc_term, 4.0*p*coslam/vp*vsp_term],
1213 [4.0*p*cosphi/vs*vsp_term, vsp_term**2 - pcc_term]],
1214 dtype=complex) / denom
1216 if not energy:
1217 return scatter
1218 else:
1219 eps = 1e-16
1220 normvec = num.array([vp*rho*cosphi+eps, vs*rho*coslam+eps])
1221 escatter = scatter*num.conj(scatter) * num.real(
1222 (normvec[:, num.newaxis]) / (normvec[num.newaxis, :]))
1223 return num.real(escatter)
1226def psv_solid_ind(in_direction, out_direction, in_mode, out_mode):
1227 '''
1228 Get indices to select the appropriate element from scatter matrix for
1229 solid-solid interface.
1230 '''
1232 return (
1233 (out_direction == DOWN)*2 + (out_mode == S),
1234 (in_direction == UP)*2 + (in_mode == S))
1237def psv_solid(material1, material2, p, energy=False):
1238 '''
1239 Scatter matrix for solid-solid interface.
1241 :param material1: material above, object of type :py:class:`Material`
1242 :param material2: material below, object of type :py:class:`Material`
1243 :param p: flat ray parameter [s/m]
1244 :param energy: bool, when ``True`` energy normalized coefficients are
1245 returned
1246 :returns: Scatter matrix
1248 The scatter matrix is ordered as follows::
1250 [[P1P1, S1P1, P2P1, S2P1],
1251 [P1S1, S1S1, P2S1, S2S1],
1252 [P1P2, S1P2, P2P2, S2P2],
1253 [P1S2, S1S2, P2S2, S2S2]]
1255 The formulas given in Aki & Richards are used.
1256 '''
1258 vp1, vs1, rho1 = material1.vp, material1.vs, material1.rho
1259 vp2, vs2, rho2 = material2.vp, material2.vs, material2.rho
1261 sinphi1 = p * vp1
1262 cosphi1 = csswap(sinphi1)
1263 sinlam1 = p * vs1
1264 coslam1 = csswap(sinlam1)
1265 sinphi2 = p * vp2
1266 cosphi2 = csswap(sinphi2)
1267 sinlam2 = p * vs2
1268 coslam2 = csswap(sinlam2)
1270 # from aki and richards
1271 M = num.array([
1272 [-vp1*p, -coslam1, vp2*p, coslam2],
1273 [cosphi1, -vs1*p, cosphi2, -vs2*p],
1274 [2.0*rho1*vs1**2*p*cosphi1, rho1*vs1*(1.0-2.0*vs1**2*p**2),
1275 2.0*rho2*vs2**2*p*cosphi2, rho2*vs2*(1.0-2.0*vs2**2*p**2)],
1276 [-rho1*vp1*(1.0-2.0*vs1**2*p**2), 2.0*rho1*vs1**2*p*coslam1,
1277 rho2*vp2*(1.0-2.0*vs2**2*p**2), -2.0*rho2*vs2**2*p*coslam2]],
1278 dtype=complex)
1279 N = M.copy()
1280 N[0] *= -1.0
1281 N[3] *= -1.0
1283 scatter = num.dot(num.linalg.inv(M), N)
1285 if not energy:
1286 return scatter
1287 else:
1288 eps = 1e-16
1289 if vs1 == 0.:
1290 vs1 = vp1*1e-16
1291 if vs2 == 0.:
1292 vs2 = vp2*1e-16
1293 normvec = num.array([
1294 vp1*rho1*(cosphi1+eps), vs1*rho1*(coslam1+eps),
1295 vp2*rho2*(cosphi2+eps), vs2*rho2*(coslam2+eps)], dtype=complex)
1296 escatter = scatter*num.conj(scatter) * num.real(
1297 normvec[:, num.newaxis] / normvec[num.newaxis, :])
1299 return num.real(escatter)
1302class BadPotIntCoefs(CakeError):
1303 pass
1306def potint_coefs(c1, c2, r1, r2): # r2 > r1
1307 eps = r2*1e-9
1308 if c1 == 0. and c2 == 0.:
1309 c1c2 = 1.
1310 else:
1311 c1c2 = c1/c2
1312 b = math.log(c1c2)/math.log((r1+eps)/r2)
1313 if abs(b) > 10.:
1314 raise BadPotIntCoefs()
1315 a = c1/(r1+eps)**b
1316 return a, b
1319def imode(s):
1320 if s.lower() == 'p':
1321 return P
1322 elif s.lower() == 's':
1323 return S
1326def smode(i):
1327 if i == P:
1328 return 'p'
1329 elif i == S:
1330 return 's'
1333class PathFailed(CakeError):
1334 pass
1337class SurfaceReached(PathFailed):
1338 pass
1341class BottomReached(PathFailed):
1342 pass
1345class MaxDepthReached(PathFailed):
1346 pass
1349class MinDepthReached(PathFailed):
1350 pass
1353class Trapped(PathFailed):
1354 pass
1357class NotPhaseConform(PathFailed):
1358 pass
1361class CannotPropagate(PathFailed):
1362 def __init__(self, direction, ilayer):
1363 PathFailed.__init__(self)
1364 self._direction = direction
1365 self._ilayer = ilayer
1367 def __str__(self):
1368 return 'Cannot enter layer %i from %s' % (
1369 self._ilayer, {
1370 UP: 'below',
1371 DOWN: 'above'}[self._direction])
1374class Layer(object):
1375 '''
1376 Representation of a layer in a layered earth model.
1378 :param ztop: depth of top of layer
1379 :param zbot: depth of bottom of layer
1380 :param name: name of layer (optional)
1382 Subclasses are: :py:class:`HomogeneousLayer` and :py:class:`GradientLayer`.
1383 '''
1385 def __init__(self, ztop, zbot, name=None):
1386 self.ztop = ztop
1387 self.zbot = zbot
1388 self.zmid = (self.ztop + self.zbot) * 0.5
1389 self.name = name
1390 self.ilayer = None
1392 def _update_potint_coefs(self):
1393 potint_p = potint_s = False
1394 try:
1395 self._ppic = potint_coefs(
1396 self.mbot.vp, self.mtop.vp,
1397 radius(self.zbot), radius(self.ztop))
1398 potint_p = True
1399 except BadPotIntCoefs:
1400 pass
1402 potint_s = False
1403 try:
1404 self._spic = potint_coefs(
1405 self.mbot.vs, self.mtop.vs,
1406 radius(self.zbot), radius(self.ztop))
1407 potint_s = True
1408 except BadPotIntCoefs:
1409 pass
1411 assert P == 1 and S == 2
1412 self._use_potential_interpolation = (None, potint_p, potint_s)
1414 def potint_coefs(self, mode):
1415 '''
1416 Get coefficients for potential interpolation.
1418 :param mode: mode of wave propagation, :py:const:`P` or :py:const:`S`
1419 :returns: coefficients ``(a, b)``
1420 '''
1422 if mode == P:
1423 return self._ppic
1424 else:
1425 return self._spic
1427 def contains(self, z):
1428 '''
1429 Tolerantly check if a given depth is within the layer
1430 (including boundaries).
1431 '''
1433 return self.ztop <= z <= self.zbot or \
1434 self.at_bottom(z) or self.at_top(z)
1436 def inner(self, z):
1437 '''
1438 Tolerantly check if a given depth is within the layer
1439 (not including boundaries).
1440 '''
1442 return self.ztop <= z <= self.zbot and not \
1443 self.at_bottom(z) and not \
1444 self.at_top(z)
1446 def at_bottom(self, z):
1447 '''
1448 Tolerantly check if given depth is at the bottom of the layer.
1449 '''
1451 return abs(self.zbot - z) < ZEPS
1453 def at_top(self, z):
1454 '''
1455 Tolerantly check if given depth is at the top of the layer.
1456 '''
1457 return abs(self.ztop - z) < ZEPS
1459 def pflat_top(self, p):
1460 '''
1461 Convert spherical ray parameter to local flat ray parameter for top of
1462 layer.
1463 '''
1464 return p / (earthradius-self.ztop)
1466 def pflat_bottom(self, p):
1467 '''
1468 Convert spherical ray parameter to local flat ray parameter for bottom
1469 of layer.
1470 '''
1471 return p / (earthradius-self.zbot)
1473 def pflat(self, p, z):
1474 '''
1475 Convert spherical ray parameter to local flat ray parameter for
1476 given depth.
1477 '''
1478 return p / (earthradius-z)
1480 def v_potint(self, mode, z):
1481 a, b = self.potint_coefs(mode)
1482 return a*(earthradius-z)**b
1484 def u_potint(self, mode, z):
1485 a, b = self.potint_coefs(mode)
1486 return 1./(a*(earthradius-z)**b)
1488 def xt_potint(self, p, mode, zpart=None):
1489 '''
1490 Get travel time and distance for for traversal with given mode and ray
1491 parameter.
1493 :param p: ray parameter (spherical)
1494 :param mode: mode of propagation (:py:const:`P` or :py:const:`S`)
1495 :param zpart: if given, tuple with two depths to restrict computation
1496 to a part of the layer
1498 This implementation uses analytic formulas valid for a spherical earth
1499 in the case where the velocity c within the layer is given by potential
1500 interpolation of the form
1502 c(z) = a*z^b
1503 '''
1504 utop, ubot = self.u_top_bottom(mode)
1505 a, b = self.potint_coefs(mode)
1506 ztop = self.ztop
1507 zbot = self.zbot
1508 if zpart is not None:
1509 utop = self.u(mode, zpart[0])
1510 ubot = self.u(mode, zpart[1])
1511 ztop, zbot = zpart
1512 utop = 1./(a*(earthradius-ztop)**b)
1513 ubot = 1./(a*(earthradius-zbot)**b)
1515 r1 = radius(zbot)
1516 r2 = radius(ztop)
1517 burger_eta1 = r1 * ubot
1518 burger_eta2 = r2 * utop
1519 if b != 1:
1520 def cpe(eta):
1521 return num.arccos(num.minimum(p/num.maximum(eta, p/2), 1.0))
1523 def sep(eta):
1524 return num.sqrt(num.maximum(eta**2 - p**2, 0.0))
1526 x = (cpe(burger_eta2)-cpe(burger_eta1))/(1.0-b)
1527 t = (sep(burger_eta2)-sep(burger_eta1))/(1.0-b)
1528 else:
1529 lr = math.log(r2/r1)
1530 sap = num.sqrt(1.0/a**2 - p**2)
1531 x = p/sap * lr
1532 t = 1./(a**2 * sap)
1534 x *= r2d
1536 return x, t
1538 def test(self, p, mode, z):
1539 '''
1540 Check if wave mode can exist for given ray parameter at given depth
1541 within the layer.
1542 '''
1543 return (self.u(mode, z)*radius(z) - p) > 0.
1545 def tests(self, p, mode):
1546 utop, ubot = self.u_top_bottom(mode)
1547 return (
1548 (utop * radius(self.ztop) - p) > 0.,
1549 (ubot * radius(self.zbot) - p) > 0.)
1551 def zturn_potint(self, p, mode):
1552 '''
1553 Get turning depth for given ray parameter and propagation mode.
1554 '''
1556 a, b = self.potint_coefs(mode)
1557 r = num.exp(num.log(a*p)/(1.0-b))
1558 return earthradius-r
1560 def propagate(self, p, mode, direction):
1561 '''
1562 Propagate ray through layer.
1564 :param p: ray parameter
1565 :param mode: propagation mode
1566 :param direction: in direction (:py:const:`UP` or :py:const:`DOWN`
1567 '''
1568 if direction == DOWN:
1569 zin, zout = self.ztop, self.zbot
1570 else:
1571 zin, zout = self.zbot, self.ztop
1573 if self.v(mode, zin) == 0.0 or not self.test(p, mode, zin):
1574 raise CannotPropagate(direction, self.ilayer)
1576 if not self.test(p, mode, zout):
1577 return -direction
1578 else:
1579 return direction
1581 def resize(self, depth_min=None, depth_max=None):
1582 '''
1583 Change layer thinkness and interpolate material if required.
1584 '''
1585 if depth_min:
1586 mtop = self.material(depth_min)
1588 if depth_max:
1589 mbot = self.material(depth_max)
1591 self.mtop = mtop if depth_min else self.mtop
1592 self.mbot = mbot if depth_max else self.mbot
1593 self.ztop = depth_min if depth_min else self.ztop
1594 self.zbot = depth_max if depth_max else self.zbot
1595 self.zmid = self.ztop + (self.zbot - self.ztop)/2.
1598class DoesNotTurn(CakeError):
1599 pass
1602def radius(z):
1603 return earthradius - z
1606class HomogeneousLayer(Layer):
1607 '''
1608 Representation of a homogeneous layer in a layered earth model.
1610 Base class: :py:class:`Layer`.
1611 '''
1613 def __init__(self, ztop, zbot, m, name=None):
1614 Layer.__init__(self, ztop, zbot, name=name)
1615 self.m = m
1616 self.mtop = m
1617 self.mbot = m
1618 self._update_potint_coefs()
1620 def copy(self, ztop=None, zbot=None):
1621 if ztop is None:
1622 ztop = self.ztop
1624 if zbot is None:
1625 zbot = self.zbot
1627 return HomogeneousLayer(ztop, zbot, self.m, name=self.name)
1629 def material(self, z):
1630 return self.m
1632 def u(self, mode, z=None):
1633 if self._use_potential_interpolation[mode] and z is not None:
1634 return self.u_potint(mode, z)
1636 if mode == P:
1637 return 1./self.m.vp
1638 if mode == S:
1639 return 1./self.m.vs
1641 def u_top_bottom(self, mode):
1642 u = self.u(mode)
1643 return u, u
1645 def v(self, mode, z=None):
1646 if self._use_potential_interpolation[mode] and z is not None:
1647 return self.v_potint(mode, z)
1649 if mode == P:
1650 v = self.m.vp
1651 if mode == S:
1652 v = self.m.vs
1654 if num.isscalar(z):
1655 return v
1656 else:
1657 return filled(v, len(z))
1659 def v_top_bottom(self, mode):
1660 v = self.v(mode)
1661 return v, v
1663 def xt(self, p, mode, zpart=None):
1664 if self._use_potential_interpolation[mode]:
1665 return self.xt_potint(p, mode, zpart)
1667 u = self.u(mode)
1668 pflat = self.pflat_bottom(p)
1669 if zpart is None:
1670 dz = (self.zbot - self.ztop)
1671 else:
1672 dz = abs(zpart[1]-zpart[0])
1674 u = self.u(mode)
1675 eps = u*0.001
1676 denom = num.sqrt(u**2 - pflat**2) + eps
1678 x = r2d*pflat/(earthradius-self.zmid) * dz / denom
1679 t = u**2 * dz / denom
1680 return x, t
1682 def zturn(self, p, mode):
1683 if self._use_potential_interpolation[mode]:
1684 return self.zturn_potint(p, mode)
1686 raise DoesNotTurn()
1688 def split(self, z):
1689 upper = HomogeneousLayer(self.ztop, z, self.m, name=self.name)
1690 lower = HomogeneousLayer(z, self.zbot, self.m, name=self.name)
1691 upper.ilayer = self.ilayer
1692 lower.ilayer = self.ilayer
1693 return upper, lower
1695 def __str__(self):
1696 if self.name:
1697 name = self.name + ' '
1698 else:
1699 name = ''
1701 calcmode = ''.join('HP'[self._use_potential_interpolation[mode]]
1702 for mode in (P, S))
1704 return ' (%i) homogeneous layer %s(%g km - %g km) [%s]\n %s' % (
1705 self.ilayer, name, self.ztop/km, self.zbot/km, calcmode, self.m)
1708class GradientLayer(Layer):
1709 '''
1710 Representation of a gradient layer in a layered earth model.
1712 Base class: :py:class:`Layer`.
1713 '''
1715 def __init__(self, ztop, zbot, mtop, mbot, name=None):
1716 Layer.__init__(self, ztop, zbot, name=name)
1717 self.mtop = mtop
1718 self.mbot = mbot
1719 self._update_potint_coefs()
1721 def copy(self, ztop=None, zbot=None):
1722 if ztop is None:
1723 ztop = self.ztop
1725 if zbot is None:
1726 zbot = self.zbot
1728 return GradientLayer(ztop, zbot, self.mtop, self.mbot, name=self.name)
1730 def interpolate(self, z, ptop, pbot):
1731 return ptop + (z - self.ztop)*(pbot - ptop)/(self.zbot-self.ztop)
1733 def material(self, z):
1734 dtop = self.mtop.astuple()
1735 dbot = self.mbot.astuple()
1736 d = [
1737 self.interpolate(z, ptop, pbot)
1738 for (ptop, pbot) in zip(dtop, dbot)]
1740 return Material(*d)
1742 def u_top_bottom(self, mode):
1743 if mode == P:
1744 return 1./self.mtop.vp, 1./self.mbot.vp
1745 if mode == S:
1746 return 1./self.mtop.vs, 1./self.mbot.vs
1748 def u(self, mode, z):
1749 if self._use_potential_interpolation[mode]:
1750 return self.u_potint(mode, z)
1752 if mode == P:
1753 return 1./self.interpolate(z, self.mtop.vp, self.mbot.vp)
1754 if mode == S:
1755 return 1./self.interpolate(z, self.mtop.vs, self.mbot.vs)
1757 def v_top_bottom(self, mode):
1758 if mode == P:
1759 return self.mtop.vp, self.mbot.vp
1760 if mode == S:
1761 return self.mtop.vs, self.mbot.vs
1763 def v(self, mode, z):
1764 if self._use_potential_interpolation[mode]:
1765 return self.v_potint(mode, z)
1767 if mode == P:
1768 return self.interpolate(z, self.mtop.vp, self.mbot.vp)
1769 if mode == S:
1770 return self.interpolate(z, self.mtop.vs, self.mbot.vs)
1772 def xt(self, p, mode, zpart=None):
1773 if self._use_potential_interpolation[mode]:
1774 return self.xt_potint(p, mode, zpart)
1776 utop, ubot = self.u_top_bottom(mode)
1777 b = (1./ubot - 1./utop)/(self.zbot - self.ztop)
1779 pflat = self.pflat_bottom(p)
1780 if zpart is not None:
1781 utop = self.u(mode, zpart[0])
1782 ubot = self.u(mode, zpart[1])
1784 peps = 1e-16
1785 pdp = pflat + peps
1787 def func(u):
1788 eta = num.sqrt(num.maximum(u**2 - pflat**2, 0.0))
1789 xx = eta/u
1790 tt = num.where(
1791 pflat <= u,
1792 num.log(u+eta) - num.log(pdp) - eta/u,
1793 0.0)
1795 return xx, tt
1797 xxtop, tttop = func(utop)
1798 xxbot, ttbot = func(ubot)
1800 x = (xxtop - xxbot) / (b*pdp)
1801 t = (tttop - ttbot) / b + pflat*x
1803 x *= r2d/(earthradius - self.zmid)
1804 return x, t
1806 def zturn(self, p, mode):
1807 if self._use_potential_interpolation[mode]:
1808 return self.zturn_potint(p, mode)
1809 pflat = self.pflat_bottom(p)
1810 vtop, vbot = self.v_top_bottom(mode)
1811 return (1./pflat - vtop) * (self.zbot - self.ztop) / \
1812 (vbot-vtop) + self.ztop
1814 def split(self, z):
1815 mmid = self.material(z)
1816 upper = GradientLayer(self.ztop, z, self.mtop, mmid, name=self.name)
1817 lower = GradientLayer(z, self.zbot, mmid, self.mbot, name=self.name)
1818 upper.ilayer = self.ilayer
1819 lower.ilayer = self.ilayer
1820 return upper, lower
1822 def __str__(self):
1823 if self.name:
1824 name = self.name + ' '
1825 else:
1826 name = ''
1828 calcmode = ''.join('HP'[self._use_potential_interpolation[mode]]
1829 for mode in (P, S))
1831 return ''' (%i) gradient layer %s(%g km - %g km) [%s]
1832 %s
1833 %s''' % (
1834 self.ilayer,
1835 name,
1836 self.ztop/km,
1837 self.zbot/km,
1838 calcmode,
1839 self.mtop,
1840 self.mbot)
1843class Discontinuity(object):
1844 '''
1845 Base class for discontinuities in layered earth model.
1847 Subclasses are: :py:class:`Interface` and :py:class:`Surface`.
1848 '''
1850 def __init__(self, z, name=None):
1851 self.z = z
1852 self.zbot = z
1853 self.ztop = z
1854 self.name = name
1856 def change_depth(self, z):
1857 self.z = z
1858 self.zbot = z
1859 self.ztop = z
1861 def copy(self):
1862 return copy.deepcopy(self)
1865class Interface(Discontinuity):
1866 '''
1867 Representation of an interface in a layered earth model.
1869 Base class: :py:class:`Discontinuity`.
1870 '''
1872 def __init__(self, z, mabove, mbelow, name=None):
1873 Discontinuity.__init__(self, z, name)
1874 self.mabove = mabove
1875 self.mbelow = mbelow
1877 def __str__(self):
1878 if self.name is None:
1879 return 'interface'
1880 else:
1881 return 'interface "%s"' % self.name
1883 def u_top_bottom(self, mode):
1884 if mode == P:
1885 return reci_or_none(self.mabove.vp), reci_or_none(self.mbelow.vp)
1886 if mode == S:
1887 return reci_or_none(self.mabove.vs), reci_or_none(self.mbelow.vs)
1889 def critical_ps(self, mode):
1890 uabove, ubelow = self.u_top_bottom(mode)
1891 return (
1892 mult_or_none(uabove, radius(self.z)),
1893 mult_or_none(ubelow, radius(self.z)))
1895 def propagate(self, p, mode, direction):
1896 uabove, ubelow = self.u_top_bottom(mode)
1897 if direction == DOWN:
1898 if ubelow is not None and ubelow*radius(self.z) - p >= 0:
1899 return direction
1900 else:
1901 return -direction
1902 if direction == UP:
1903 if uabove is not None and uabove*radius(self.z) - p >= 0:
1904 return direction
1905 else:
1906 return -direction
1908 def pflat(self, p):
1909 return p / (earthradius-self.z)
1911 def efficiency(self, in_direction, out_direction, in_mode, out_mode, p):
1912 scatter = psv_solid(
1913 self.mabove, self.mbelow, self.pflat(p), energy=True)
1914 return scatter[
1915 psv_solid_ind(in_direction, out_direction, in_mode, out_mode)]
1918class Surface(Discontinuity):
1919 '''
1920 Representation of the surface discontinuity in a layered earth model.
1922 Base class: :py:class:`Discontinuity`.
1923 '''
1925 def __init__(self, z, mbelow):
1926 Discontinuity.__init__(self, z, 'surface')
1927 self.z = z
1928 self.mbelow = mbelow
1930 def propagate(self, p, mode, direction):
1931 return direction # no implicit reflection at surface
1933 def u_top_bottom(self, mode):
1934 if mode == P:
1935 return None, reci_or_none(self.mbelow.vp)
1936 if mode == S:
1937 return None, reci_or_none(self.mbelow.vs)
1939 def critical_ps(self, mode):
1940 _, ubelow = self.u_top_bottom(mode)
1941 return None, mult_or_none(ubelow, radius(self.z))
1943 def pflat(self, p):
1944 return p / (earthradius-self.z)
1946 def efficiency(self, in_direction, out_direction, in_mode, out_mode, p):
1947 if in_direction == DOWN or out_direction == UP:
1948 return 0.0
1949 else:
1950 return psv_surface(
1951 self.mbelow, self.pflat(p), energy=True)[
1952 psv_surface_ind(in_mode, out_mode)]
1954 def __str__(self):
1955 return 'surface'
1958class Walker(object):
1959 def __init__(self, elements):
1960 self._elements = elements
1961 self._i = 0
1963 def current(self):
1964 return self._elements[self._i]
1966 def go(self, direction):
1967 if direction == UP:
1968 self.up()
1969 else:
1970 self.down()
1972 def down(self):
1973 if self._i < len(self._elements)-1:
1974 self._i += 1
1975 else:
1976 raise BottomReached()
1978 def up(self):
1979 if self._i > 0:
1980 self._i -= 1
1981 else:
1982 raise SurfaceReached()
1984 def goto_layer(self, layer):
1985 self._i = self._elements.index(layer)
1988class RayElement(object):
1989 '''
1990 An element of a :py:class:`RayPath`.
1991 '''
1993 def __eq__(self, other):
1994 return type(self) == type(other) and self.__dict__ == other.__dict__
1996 def is_straight(self):
1997 return isinstance(self, Straight)
1999 def is_kink(self):
2000 return isinstance(self, Kink)
2003class Straight(RayElement):
2004 '''
2005 A ray segment representing wave propagation through one :py:class:`Layer`.
2006 '''
2008 def __init__(self, direction_in, direction_out, mode, layer):
2009 self.mode = mode
2010 self._direction_in = direction_in
2011 self._direction_out = direction_out
2012 self.layer = layer
2014 def angle_in(self, p, endgaps=None):
2015 z = self.z_in(endgaps)
2016 dir = self.eff_direction_in(endgaps)
2017 v = self.layer.v(self.mode, z)
2018 pf = self.layer.pflat(p, z)
2020 if dir == DOWN:
2021 return num.arcsin(v*pf)*r2d
2022 else:
2023 return 180.-num.arcsin(v*pf)*r2d
2025 def angle_out(self, p, endgaps=None):
2026 z = self.z_out(endgaps)
2027 dir = self.eff_direction_out(endgaps)
2028 v = self.layer.v(self.mode, z)
2029 pf = self.layer.pflat(p, z)
2031 if dir == DOWN:
2032 return 180.-num.arcsin(v*pf)*r2d
2033 else:
2034 return num.arcsin(v*pf)*r2d
2036 def pflat_in(self, p, endgaps=None):
2037 return p / (earthradius-self.z_in(endgaps))
2039 def pflat_out(self, p, endgaps=None):
2040 return p / (earthradius-self.z_out(endgaps))
2042 def test(self, p, z):
2043 return self.layer.test(p, self.mode, z)
2045 def z_in(self, endgaps=None):
2046 if endgaps is not None:
2047 return endgaps[0]
2048 else:
2049 lyr = self.layer
2050 return (lyr.ztop, lyr.zbot)[self._direction_in == UP]
2052 def z_out(self, endgaps=None):
2053 if endgaps is not None:
2054 return endgaps[1]
2055 else:
2056 lyr = self.layer
2057 return (lyr.ztop, lyr.zbot)[self._direction_out == DOWN]
2059 def turns(self):
2060 return self._direction_in != self._direction_out
2062 def eff_direction_in(self, endgaps=None):
2063 if endgaps is None:
2064 return self._direction_in
2065 else:
2066 return endgaps[2]
2068 def eff_direction_out(self, endgaps=None):
2069 if endgaps is None:
2070 return self._direction_out
2071 else:
2072 return endgaps[3]
2074 def zturn(self, p):
2075 lyr = self.layer
2076 return lyr.zturn(p, self.mode)
2078 def u_in(self, endgaps=None):
2079 return self.layer.u(self.mode, self.z_in(endgaps))
2081 def u_out(self, endgaps=None):
2082 return self.layer.u(self.mode, self.z_out(endgaps))
2084 def critical_p_in(self, endgaps=None):
2085 z = self.z_in(endgaps)
2086 return self.layer.u(self.mode, z)*radius(z)
2088 def critical_p_out(self, endgaps=None):
2089 z = self.z_out(endgaps)
2090 return self.layer.u(self.mode, z)*radius(z)
2092 def xt(self, p, zpart=None):
2093 x, t = self.layer.xt(p, self.mode, zpart=zpart)
2094 if self._direction_in != self._direction_out and zpart is None:
2095 x *= 2.
2096 t *= 2.
2097 return x, t
2099 def xt_gap(self, p, zstart, zstop, samedir):
2100 z1, z2 = zstart, zstop
2101 if z1 > z2:
2102 z1, z2 = z2, z1
2104 x, t = self.layer.xt(p, self.mode, zpart=(z1, z2))
2106 if samedir:
2107 return x, t
2108 else:
2109 xfull, tfull = self.xt(p)
2110 return xfull-x, tfull-t
2112 def __hash__(self):
2113 return hash((
2114 self._direction_in,
2115 self._direction_out,
2116 self.mode,
2117 id(self.layer)))
2120class HeadwaveStraight(Straight):
2121 def __init__(self, direction_in, direction_out, mode, interface):
2122 Straight.__init__(self, direction_in, direction_out, mode, None)
2124 self.interface = interface
2126 def z_in(self, zpart=None):
2127 return self.interface.z
2129 def z_out(self, zpart=None):
2130 return self.interface.z
2132 def zturn(self, p):
2133 return filled(self.interface.z, len(p))
2135 def xt(self, p, zpart=None):
2136 return 0., 0.
2138 def x2t_headwave(self, xstretch):
2139 xstretch_m = xstretch*d2r*radius(self.interface.z)
2140 return min_not_none(*self.interface.u_top_bottom(self.mode))*xstretch_m
2143class Kink(RayElement):
2144 '''
2145 An interaction of a ray with a :py:class:`Discontinuity`.
2146 '''
2148 def __init__(
2149 self,
2150 in_direction,
2151 out_direction,
2152 in_mode,
2153 out_mode,
2154 discontinuity):
2156 self.in_direction = in_direction
2157 self.out_direction = out_direction
2158 self.in_mode = in_mode
2159 self.out_mode = out_mode
2160 self.discontinuity = discontinuity
2162 def reflection(self):
2163 return self.in_direction != self.out_direction
2165 def conversion(self):
2166 return self.in_mode != self.out_mode
2168 def efficiency(self, p, out_direction=None, out_mode=None):
2170 if out_direction is None:
2171 out_direction = self.out_direction
2173 if out_mode is None:
2174 out_mode = self.out_mode
2176 return self.discontinuity.efficiency(
2177 self.in_direction, out_direction, self.in_mode, out_mode, p)
2179 def __str__(self):
2180 r, c = self.reflection(), self.conversion()
2181 if r and c:
2182 return '|~'
2183 if r:
2184 return '|'
2185 if c:
2186 return '~'
2187 return '_'
2189 def __hash__(self):
2190 return hash((
2191 self.in_direction,
2192 self.out_direction,
2193 self.in_mode,
2194 self.out_mode,
2195 id(self.discontinuity)))
2198class PRangeNotSet(CakeError):
2199 pass
2202class RayPath(object):
2203 '''
2204 Representation of a fan of rays running through a common sequence of
2205 layers / interfaces.
2206 '''
2208 def __init__(self, phase):
2209 self.elements = []
2210 self.phase = phase
2211 self._pmax = None
2212 self._pmin = None
2213 self._p = None
2214 self._is_headwave = False
2216 def set_is_headwave(self, is_headwave):
2217 self._is_headwave = is_headwave
2219 def copy(self):
2220 '''
2221 Get a copy of it.
2222 '''
2224 c = copy.copy(self)
2225 c.elements = list(self.elements)
2226 return c
2228 def endgaps(self, zstart, zstop):
2229 '''
2230 Get information needed for end point adjustments.
2231 '''
2233 return (
2234 zstart,
2235 zstop,
2236 self.phase.direction_start(),
2237 self.phase.direction_stop())
2239 def append(self, element):
2240 self.elements.append(element)
2242 def _check_have_prange(self):
2243 if self._pmax is None:
2244 raise PRangeNotSet()
2246 def set_prange(self, pmin, pmax, dp):
2247 self._pmin, self._pmax = pmin, pmax
2248 self._prange_dp = dp
2250 def used_phase(self, p=None, eps=1.):
2251 '''
2252 Calculate phase definition from ray path.
2253 '''
2255 used = PhaseDef()
2256 fleg = self.phase.first_leg()
2257 used.append(Leg(fleg.departure, fleg.mode))
2258 n_elements_n = [None] + self.elements + [None]
2259 for before, element, after in zip(
2260 n_elements_n[:-2],
2261 n_elements_n[1:-1],
2262 n_elements_n[2:]):
2264 if element.is_kink() and HeadwaveStraight not in (
2265 type(before),
2266 type(after)):
2268 if element.reflection() or element.conversion():
2269 z = element.discontinuity.z
2270 used.append(Knee(
2271 z,
2272 element.in_direction,
2273 element.out_direction != element.in_direction,
2274 element.in_mode,
2275 element.out_mode))
2277 used.append(Leg(element.out_direction, element.out_mode))
2279 elif type(element) is HeadwaveStraight:
2280 z = element.interface.z
2281 k = Knee(
2282 z,
2283 before.in_direction,
2284 after.out_direction != before.in_direction,
2285 before.in_mode,
2286 after.out_mode)
2288 k.headwave = True
2289 used.append(k)
2290 used.append(Leg(after.out_direction, after.out_mode))
2292 if (p is not None and before and after
2293 and element.is_straight()
2294 and before.is_kink()
2295 and after.is_kink()
2296 and element.turns()
2297 and not before.reflection() and not before.conversion()
2298 and not after.reflection() and not after.conversion()):
2300 ai = element.angle_in(p)
2301 if 90.0-eps < ai and ai < 90+eps:
2302 used.append(
2303 Head(
2304 before.discontinuity.z,
2305 before.out_direction,
2306 element.mode))
2307 used.append(
2308 Leg(-before.out_direction, element.mode))
2310 used._direction_stop = self.phase.direction_stop()
2311 used._definition = self.phase.definition()
2313 return used
2315 def pmax(self):
2316 '''
2317 Get maximum valid ray parameter.
2318 '''
2319 self._check_have_prange()
2320 return self._pmax
2322 def pmin(self):
2323 '''
2324 Get minimum valid ray parameter.
2325 '''
2326 self._check_have_prange()
2327 return self._pmin
2329 def xmin(self):
2330 '''
2331 Get minimal distance.
2332 '''
2333 self._analyse()
2334 return self._xmin
2336 def xmax(self):
2337 '''
2338 Get maximal distance.
2339 '''
2340 self._analyse()
2341 return self._xmax
2343 def kinks(self):
2344 '''
2345 Iterate over propagation mode changes (reflections/transmissions).
2346 '''
2347 return (k for k in self.elements if isinstance(k, Kink))
2349 def straights(self):
2350 '''
2351 Iterate over ray segments.
2352 '''
2353 return (s for s in self.elements if isinstance(s, Straight))
2355 def headwave_straight(self):
2356 for s in self.elements:
2357 if type(s) is HeadwaveStraight:
2358 return s
2360 def first_straight(self):
2361 '''
2362 Get first ray segment.
2363 '''
2364 for s in self.elements:
2365 if isinstance(s, Straight):
2366 return s
2368 def last_straight(self):
2369 '''
2370 Get last ray segment.
2371 '''
2372 for s in reversed(self.elements):
2373 if isinstance(s, Straight):
2374 return s
2376 def efficiency(self, p):
2377 '''
2378 Get product of all conversion/reflection coefficients encountered on
2379 path.
2380 '''
2381 return reduce(
2382 operator.mul, (k.efficiency(p) for k in self.kinks()), 1.)
2384 def spreading(self, p, endgaps):
2385 '''
2386 Get geometrical spreading factor.
2387 '''
2388 if self._is_headwave:
2389 return 0.0
2391 self._check_have_prange()
2392 dp = self._prange_dp * 0.01
2393 assert self._pmax - self._pmin > dp
2395 if p + dp > self._pmax:
2396 p = p-dp
2398 x0, t = self.xt(p, endgaps)
2399 x1, t = self.xt(p+dp, endgaps)
2400 x0 *= d2r
2401 x1 *= d2r
2402 if x1 == x0:
2403 return num.nan
2405 dp_dx = dp/(x1-x0)
2407 x = x0
2408 if x == 0.:
2409 x = x1
2410 p = dp
2412 first = self.first_straight()
2413 last = self.last_straight()
2414 return num.abs(dp_dx) * first.pflat_in(p, endgaps) / (
2415 4.0 * math.pi * num.sin(x) *
2416 (earthradius-first.z_in(endgaps)) *
2417 (earthradius-last.z_out(endgaps))**2 *
2418 first.u_in(endgaps)**2 *
2419 num.abs(num.cos(first.angle_in(p, endgaps)*d2r)) *
2420 num.abs(num.cos(last.angle_out(p, endgaps)*d2r)))
2422 def make_p(self, dp=None, n=None, nmin=None):
2423 assert dp is None or n is None
2425 if self._pmin == self._pmax:
2426 return num.array([self._pmin])
2428 if dp is None:
2429 dp = self._prange_dp
2431 if n is None:
2432 n = int(round((self._pmax-self._pmin)/dp)) + 1
2434 if nmin is not None:
2435 n = max(n, nmin)
2437 ppp = num.linspace(self._pmin, self._pmax, n)
2438 return ppp
2440 def xt_endgaps(self, p, endgaps, which='both'):
2441 '''
2442 Get amount of distance/traveltime to be subtracted at the generic ray
2443 path's ends.
2444 '''
2446 zstart, zstop, dirstart, dirstop = endgaps
2447 firsts = self.first_straight()
2448 lasts = self.last_straight()
2449 xs, ts = firsts.xt_gap(
2450 p, zstart, firsts.z_in(), dirstart == firsts._direction_in)
2451 xe, te = lasts.xt_gap(
2452 p, zstop, lasts.z_out(), dirstop == lasts._direction_out)
2454 if which == 'both':
2455 return xs + xe, ts + te
2456 elif which == 'left':
2457 return xs, ts
2458 elif which == 'right':
2459 return xe, te
2461 def xt_endgaps_ptest(self, p, endgaps):
2462 '''
2463 Check if ray parameter is valid at source and receiver.
2464 '''
2466 zstart, zstop, dirstart, dirstop = endgaps
2467 firsts = self.first_straight()
2468 lasts = self.last_straight()
2469 return num.logical_and(firsts.test(p, zstart), lasts.test(p, zstop))
2471 def xt(self, p, endgaps):
2472 '''
2473 Calculate distance and traveltime for given ray parameter.
2474 '''
2476 if isinstance(p, num.ndarray):
2477 sx = num.zeros(p.size)
2478 st = num.zeros(p.size)
2479 else:
2480 sx = 0.0
2481 st = 0.0
2483 for s in self.straights():
2484 x, t = s.xt(p)
2485 sx += x
2486 st += t
2488 if endgaps:
2489 dx, dt = self.xt_endgaps(p, endgaps)
2490 sx -= dx
2491 st -= dt
2493 return sx, st
2495 def xt_limits(self, p):
2496 '''
2497 Calculate limits of distance and traveltime for given ray parameter.
2498 '''
2500 if isinstance(p, num.ndarray):
2501 sx = num.zeros(p.size)
2502 st = num.zeros(p.size)
2503 sxe = num.zeros(p.size)
2504 ste = num.zeros(p.size)
2505 else:
2506 sx = 0.0
2507 st = 0.0
2508 sxe = 0.0
2509 ste = 0.0
2511 sfirst = self.first_straight()
2512 slast = self.last_straight()
2514 for s in self.straights():
2515 if s is not sfirst and s is not slast:
2516 x, t = s.xt(p)
2517 sx += x
2518 st += t
2520 sends = [sfirst]
2521 if sfirst is not slast:
2522 sends.append(slast)
2524 for s in sends:
2525 x, t = s.xt(p)
2526 sxe += x
2527 ste += t
2529 return sx, (sx + sxe), st, (st + ste)
2531 def iter_zxt(self, p):
2532 '''
2533 Iterate over (depth, distance, traveltime) at each layer interface on
2534 ray path.
2535 '''
2537 sx = num.zeros(p.size)
2538 st = num.zeros(p.size)
2539 ok = False
2540 for s in self.straights():
2541 yield s.z_in(), sx.copy(), st.copy()
2543 x, t = s.xt(p)
2544 sx += x
2545 st += t
2546 ok = True
2548 if ok:
2549 yield s.z_out(), sx.copy(), st.copy()
2551 def zxt_path_subdivided(
2552 self, p, endgaps,
2553 points_per_straight=20,
2554 x_for_headwave=None):
2556 '''
2557 Get geometrical representation of ray path.
2558 '''
2560 if self._is_headwave:
2561 assert p.size == 1
2562 x, t = self.xt(p, endgaps)
2563 xstretch = x_for_headwave-x
2564 nout = xstretch.size
2565 else:
2566 nout = p.size
2568 dxl, dtl = self.xt_endgaps(p, endgaps, which='left')
2569 dxr, dtr = self.xt_endgaps(p, endgaps, which='right')
2571 # first create full path including the endgaps
2572 sx = num.zeros(nout) - dxl
2573 st = num.zeros(nout) - dtl
2574 zxt = []
2575 for s in self.straights():
2576 n = points_per_straight
2578 back = None
2579 zin, zout = s.z_in(), s.z_out()
2580 if type(s) is HeadwaveStraight:
2581 z = zin
2582 for i in range(n):
2583 xs = float(i)/(n-1) * xstretch
2584 ts = s.x2t_headwave(xs)
2585 zxt.append((filled(z, xstretch.size), sx+xs, st+ts))
2586 else:
2587 if zin != zout: # normal traversal
2588 zs = num.linspace(zin, zout, n).tolist()
2589 for z in zs:
2590 x, t = s.xt(p, zpart=sorted([zin, z]))
2591 zxt.append((filled(z, nout), sx + x, st + t))
2593 else: # ray turns in layer
2594 zturn = s.zturn(p)
2595 back = []
2596 for i in range(n):
2597 z = zin + (zturn - zin) * num.sin(
2598 float(i)/(n-1)*math.pi/2.0) * 0.999
2600 if zturn[0] >= zin:
2601 x, t = s.xt(p, zpart=[zin, z])
2602 else:
2603 x, t = s.xt(p, zpart=[z, zin])
2604 zxt.append((z, sx + x, st + t))
2605 back.append((z, x, t))
2607 if type(s) is HeadwaveStraight:
2608 x = xstretch
2609 t = s.x2t_headwave(xstretch)
2610 else:
2611 x, t = s.xt(p)
2613 sx += x
2614 st += t
2615 if back:
2616 for z, x, t in reversed(back):
2617 zxt.append((z, sx - x, st - t))
2619 # gather results as arrays with such that x[ip, ipoint]
2620 fanz, fanx, fant = [], [], []
2621 for z, x, t in zxt:
2622 fanz.append(z)
2623 fanx.append(x)
2624 fant.append(t)
2626 z = num.array(fanz).T
2627 x = num.array(fanx).T
2628 t = num.array(fant).T
2630 # cut off the endgaps, add exact endpoints
2631 xmax = x[:, -1] - dxr
2632 tmax = t[:, -1] - dtr
2633 zstart, zstop = endgaps[:2]
2634 zs, xs, ts = [], [], []
2635 for i in range(nout):
2636 t_ = t[i]
2637 indices = num.where(num.logical_and(0. <= t_, t_ <= tmax[i]))[0]
2638 n = indices.size + 2
2639 zs_, xs_, ts_ = [num.empty(n, dtype=float) for j in range(3)]
2640 zs_[1:-1] = z[i, indices]
2641 xs_[1:-1] = x[i, indices]
2642 ts_[1:-1] = t[i, indices]
2643 zs_[0], zs_[-1] = zstart, zstop
2644 xs_[0], xs_[-1] = 0., xmax[i]
2645 ts_[0], ts_[-1] = 0., tmax[i]
2646 zs.append(zs_)
2647 xs.append(xs_)
2648 ts.append(ts_)
2650 return zs, xs, ts
2652 def _analyse(self):
2653 if self._p is not None:
2654 return
2656 p = self.make_p(nmin=20)
2657 xmin, xmax, tmin, tmax = self.xt_limits(p)
2659 self._x, self._t, self._p = xmax, tmax, p
2660 self._xmin, self._xmax = xmin.min(), xmax.max()
2661 self._tmin, self._tmax = tmin.min(), tmax.max()
2663 def draft_pxt(self, endgaps):
2664 self._analyse()
2666 if not self._is_headwave:
2667 cp, cx, ct = self._p, self._x, self._t
2668 pcrit = min(
2669 self.critical_pstart(endgaps),
2670 self.critical_pstop(endgaps))
2672 if pcrit < self._pmin:
2673 empty = num.array([], dtype=float)
2674 return empty, empty, empty
2676 elif pcrit >= self._pmax:
2677 dx, dt = self.xt_endgaps(cp, endgaps)
2678 return cp, cx-dx, ct-dt
2680 else:
2681 n = num.searchsorted(cp, pcrit) + 1
2682 rp, rx, rt = num.empty((3, n), dtype=float)
2683 rp[:-1] = cp[:n-1]
2684 rx[:-1] = cx[:n-1]
2685 rt[:-1] = ct[:n-1]
2686 rp[-1] = pcrit
2687 rx[-1], rt[-1] = self.xt(pcrit, endgaps)
2688 dx, dt = self.xt_endgaps(rp, endgaps)
2689 rx[:-1] -= dx[:-1]
2690 rt[:-1] -= dt[:-1]
2691 return rp, rx, rt
2693 else:
2694 dx, dt = self.xt_endgaps(self._p, endgaps)
2695 p, x, t = self._p, self._x - dx, self._t - dt
2696 p, x, t = p[0], x[0], t[0]
2697 xh = num.linspace(0., x*10-x, 10)
2698 th = self.headwave_straight().x2t_headwave(xh)
2699 return filled(p, xh.size), x+xh, t+th
2701 def interpolate_x2pt_linear(self, x, endgaps):
2702 '''
2703 Get approximate ray parameter and traveltime for distance.
2704 '''
2706 self._analyse()
2708 if self._is_headwave:
2709 dx, dt = self.xt_endgaps(self._p, endgaps)
2710 xmin = self._x[0] - dx[0]
2711 tmin = self._t[0] - dt[0]
2712 el = self.headwave_straight()
2713 xok = x[x >= xmin]
2714 th = el.x2t_headwave(xstretch=(xok-xmin)) + tmin
2715 return [
2716 (x_, self._p[0], t, None) for (x_, t) in zip(xok, th)]
2718 else:
2719 if num.all(x < self._xmin) or num.all(self._xmax < x):
2720 return []
2722 rp, rx, rt = self.draft_pxt(endgaps)
2724 xp = interp(x, rx, rp, 0)
2725 xt = interp(x, rx, rt, 0)
2727 if (rp.size and
2728 len(xp) == 0 and
2729 rx[0] == 0.0 and
2730 any(x == 0.0) and
2731 rp[0] == 0.0):
2733 xp = [(0.0, rp[0])]
2734 xt = [(0.0, rt[0])]
2736 return [
2737 (x_, p, t, (rp, rx, rt)) for ((x_, p), (_, t)) in zip(xp, xt)]
2739 def __eq__(self, other):
2740 if len(self.elements) != len(other.elements):
2741 return False
2743 return all(a == b for a, b in zip(self.elements, other.elements))
2745 def __hash__(self):
2746 return hash(
2747 tuple(hash(x) for x in self.elements) +
2748 (self.phase.definition(), ))
2750 def __str__(self, p=None, eps=1.):
2751 x = []
2752 start_i = None
2753 end_i = None
2754 turn_i = None
2756 def append_layers(si, ei, ti):
2757 if si == ei and (ti is None or ti == si):
2758 x.append('%i' % si)
2759 else:
2760 if ti is not None:
2761 x.append('(%i-%i-%i)' % (si, ti, ei))
2762 else:
2763 x.append('(%i-%i)' % (si, ei))
2765 for el in self.elements:
2766 if type(el) is Straight:
2767 if start_i is None:
2768 start_i = el.layer.ilayer
2769 if el._direction_in != el._direction_out:
2770 turn_i = el.layer.ilayer
2771 end_i = el.layer.ilayer
2773 elif isinstance(el, Kink):
2774 if start_i is not None:
2775 append_layers(start_i, end_i, turn_i)
2776 start_i = None
2777 turn_i = None
2779 x.append(str(el))
2781 if start_i is not None:
2782 append_layers(start_i, end_i, turn_i)
2784 su = '(%s)' % self.used_phase(p=p, eps=eps).used_repr()
2786 return '%-15s %-17s %s' % (self.phase.definition(), su, ''.join(x))
2788 def critical_pstart(self, endgaps):
2789 '''
2790 Get critical ray parameter for source depth choice.
2791 '''
2793 return self.first_straight().critical_p_in(endgaps)
2795 def critical_pstop(self, endgaps):
2796 '''
2797 Get critical ray parameter for receiver depth choice.
2798 '''
2800 return self.last_straight().critical_p_out(endgaps)
2802 def ranges(self, endgaps):
2803 '''
2804 Get valid ranges of ray parameter, distance, and traveltime.
2805 '''
2806 p, x, t = self.draft_pxt(endgaps)
2807 return p.min(), p.max(), x.min(), x.max(), t.min(), t.max()
2809 def describe(self, endgaps=None, as_degrees=False):
2810 '''
2811 Get textual representation.
2812 '''
2814 self._analyse()
2816 if as_degrees:
2817 xunit = 'deg'
2818 xfact = 1.
2819 else:
2820 xunit = 'km'
2821 xfact = d2m/km
2823 sg = ''' Ranges for all depths in source and receiver layers:
2824 - x [%g, %g] %s
2825 - t [%g, %g] s
2826 - p [%g, %g] s/deg
2827''' % (
2828 self._xmin*xfact,
2829 self._xmax*xfact,
2830 xunit,
2831 self._tmin,
2832 self._tmax,
2833 self._pmin/r2d,
2834 self._pmax/r2d)
2836 if endgaps is not None:
2837 pmin, pmax, xmin, xmax, tmin, tmax = self.ranges(endgaps)
2838 ss = ''' Ranges for given source and receiver depths:
2839\n - x [%g, %g] %s
2840\n - t [%g, %g] s
2841\n - p [%g, %g] s/deg
2842\n''' % (xmin*xfact, xmax*xfact, xunit, tmin, tmax, pmin/r2d, pmax/r2d)
2844 else:
2845 ss = ''
2847 return '%s\n' % self + ss + sg
2850class RefineFailed(CakeError):
2851 pass
2854class Ray(object):
2855 '''
2856 Representation of a ray with a specific (path, ray parameter, distance,
2857 arrival time) choice.
2859 **Attributes:**
2861 .. py:attribute:: path
2863 :py:class:`RayPath` object containing complete propagation history.
2865 .. py:attribute:: p
2867 Ray parameter (spherical) [s/rad]
2869 .. py:attribute:: x
2871 Radial distance [deg]
2873 .. py:attribute:: t
2875 Traveltime [s]
2877 .. py:attribute:: endgaps
2879 Needed for source/receiver depth adjustments in many
2880 :py:class:`RayPath` methods.
2881 '''
2883 def __init__(self, path, p, x, t, endgaps, draft_pxt):
2884 self.path = path
2885 self.p = p
2886 self.x = x
2887 self.t = t
2888 self.endgaps = endgaps
2889 self.draft_pxt = draft_pxt
2891 def given_phase(self):
2892 '''
2893 Get phase definition which was used to create the ray.
2895 :returns: :py:class:`PhaseDef` object
2896 '''
2898 return self.path.phase
2900 def used_phase(self):
2901 '''
2902 Compute phase definition from propagation path.
2904 :returns: :py:class:`PhaseDef` object
2905 '''
2907 return self.path.used_phase(self.p)
2909 def refine(self):
2910 if self.path._is_headwave:
2911 return
2913 if self.t == 0.0 and self.p == 0.0 and self.x == 0.0:
2914 return
2916 cp, cx, ct = self.draft_pxt
2917 ip = num.searchsorted(cp, self.p)
2918 if not (0 < ip < cp.size):
2919 raise RefineFailed()
2921 pl, ph = cp[ip-1], cp[ip]
2922 p_to_t = {}
2923 i = [0]
2925 def f(p):
2926 i[0] += 1
2927 x, t = self.path.xt(p, self.endgaps)
2928 p_to_t[p] = t
2929 return self.x - x
2931 try:
2932 self.p = brentq(f, pl, ph)
2933 self.t = p_to_t[self.p]
2935 except ValueError:
2936 raise RefineFailed()
2938 def takeoff_angle(self):
2939 '''
2940 Get takeoff angle of ray.
2942 The angle is returned in [degrees].
2943 '''
2945 return self.path.first_straight().angle_in(self.p, self.endgaps)
2947 def incidence_angle(self):
2948 '''
2949 Get incidence angle of ray.
2951 The angle is returned in [degrees].
2952 '''
2954 return self.path.last_straight().angle_out(self.p, self.endgaps)
2956 def efficiency(self):
2957 '''
2958 Get conversion/reflection efficiency of the ray.
2960 A value between 0 and 1 is returned, reflecting the relative amount of
2961 energy which is transmitted along the ray and not lost by reflections
2962 or conversions.
2963 '''
2965 return self.path.efficiency(self.p)
2967 def spreading(self):
2968 '''
2969 Get geometrical spreading factor.
2970 '''
2972 return self.path.spreading(self.p, self.endgaps)
2974 def surface_sphere(self):
2975 x1, y1 = 0., earthradius - self.endgaps[0]
2976 r2 = earthradius - self.endgaps[1]
2977 x2, y2 = r2*math.sin(self.x*d2r), r2*math.cos(self.x*d2r)
2978 return ((x2-x1)**2 + (y2-y1)**2)*4.0*math.pi
2980 def zxt_path_subdivided(self, points_per_straight=20):
2981 '''
2982 Get geometrical representation of ray path.
2984 Three arrays (depth, distance, time) with points on the ray's path of
2985 propagation are returned. The number of points which are used in each
2986 ray segment (passage through one layer) may be controlled by the
2987 ``points_per_straight`` parameter.
2988 '''
2989 return self.path.zxt_path_subdivided(
2990 num.atleast_1d(self.p), self.endgaps,
2991 points_per_straight=points_per_straight,
2992 x_for_headwave=num.atleast_1d(self.x))
2994 def __str__(self, as_degrees=False):
2995 if as_degrees:
2996 sd = '%6.3g deg' % self.x
2997 else:
2998 sd = '%7.5g km' % (self.x*(d2r*earthradius/km))
3000 return '%7.5g s/deg %s %6.4g s %5.1f %5.1f %3.0f%% %3.0f%% %s' % (
3001 self.p/r2d,
3002 sd,
3003 self.t,
3004 self.takeoff_angle(),
3005 self.incidence_angle(),
3006 100*self.efficiency(),
3007 100*self.spreading()*self.surface_sphere(),
3008 self.path.__str__(p=self.p))
3011def anything_to_crust2_profile(crust2_profile):
3012 from pyrocko.dataset import crust2x2
3013 if isinstance(crust2_profile, tuple):
3014 lat, lon = [float(x) for x in crust2_profile]
3015 return crust2x2.get_profile(lat, lon)
3016 elif isinstance(crust2_profile, (str, newstr)):
3017 return crust2x2.get_profile(crust2_profile)
3018 elif isinstance(crust2_profile, crust2x2.Crust2Profile):
3019 return crust2_profile
3020 else:
3021 assert False, 'crust2_profile must be (lat, lon) a profile ' \
3022 'key or a crust2x2 Profile object)'
3025class DiscontinuityNotFound(CakeError):
3026 def __init__(self, depth_or_name):
3027 CakeError.__init__(self)
3028 self.depth_or_name = depth_or_name
3030 def __str__(self):
3031 return 'Cannot find discontinuity from given depth or name: %s' % \
3032 self.depth_or_name
3035class LayeredModelError(CakeError):
3036 pass
3039class LayeredModel(object):
3040 '''
3041 Representation of a layer cake model.
3043 There are several ways to initialize an instance of this class.
3045 1. Use the module function :py:func:`load_model` to read a model from a
3046 file.
3047 2. Create an empty model with the default constructor and append layers and
3048 discontinuities with the :py:meth:`append` method (from top to bottom).
3049 3. Use the constructor :py:meth:`LayeredModel.from_scanlines`, to
3050 automatically create the :py:class:`Layer` and :py:class:`Discontinuity`
3051 objects from a given velocity profile.
3053 An earth model is represented by as stack of :py:class:`Layer` and
3054 :py:class:`Discontinuity` objects. The method :py:meth:`arrivals` returns
3055 :py:class:`Ray` objects which may be e.g. queried for arrival times of
3056 specific phases. Each ray is associated with a :py:class:`RayPath` object.
3057 Ray objects share common ray paths if they have the same
3058 conversion/reflection/propagation history. Creating the ray path objects is
3059 relatively expensive (this is done in :py:meth:`gather_paths`), but they
3060 are cached for reuse in successive invocations.
3061 '''
3063 def __init__(self):
3064 self._surface_material = None
3065 self._elements = []
3066 self.nlayers = 0
3067 self._np = 10000
3068 self._pdepth = 5
3069 self._pathcache = {}
3071 def copy_with_elevation(self, elevation):
3072 '''
3073 Get copy of the model with surface layer stretched to given elevation.
3075 :param elevation: new surface elevation in [m]
3077 Elevation is positiv upward, contrary to the layered models downward
3078 `z` axis.
3079 '''
3081 c = copy.deepcopy(self)
3082 c._pathcache = {}
3083 surface = c._elements[0]
3084 toplayer = c._elements[1]
3086 assert toplayer.zbot > -elevation
3088 surface.z = -elevation
3089 c._elements[1] = toplayer.copy(ztop=-elevation)
3090 c._elements[1].ilayer = 0
3091 return c
3093 def zeq(self, z1, z2):
3094 return abs(z1-z2) < ZEPS
3096 def append(self, element):
3097 '''
3098 Add a layer or discontinuity at bottom of model.
3100 :param element: object of subclass of :py:class:`Layer` or
3101 :py:class:`Discontinuity`.
3102 '''
3104 if isinstance(element, Layer):
3105 if element.zbot >= earthradius:
3106 element.zbot = earthradius - 1.
3108 if element.ztop >= earthradius:
3109 raise CakeError('Layer deeper than earthradius')
3111 element.ilayer = self.nlayers
3112 self.nlayers += 1
3114 self._elements.append(element)
3116 def elements(self, direction=DOWN):
3117 '''
3118 Iterate over all elements of the model.
3120 :param direction: direction of traversal :py:const:`DOWN` or
3121 :py:const:`UP`.
3123 Objects derived from the :py:class:`Discontinuity` and
3124 :py:class:`Layer` classes are yielded.
3125 '''
3127 if direction == DOWN:
3128 return iter(self._elements)
3129 else:
3130 return reversed(self._elements)
3132 def layers(self, direction=DOWN):
3133 '''
3134 Iterate over all layers of model.
3136 :param direction: direction of traversal :py:const:`DOWN` or
3137 :py:const:`UP`.
3139 Objects derived from the :py:class:`Layer` class are yielded.
3140 '''
3142 if direction == DOWN:
3143 return (el for el in self._elements if isinstance(el, Layer))
3144 else:
3145 return (
3146 el for el in reversed(self._elements) if isinstance(el, Layer))
3148 def layer(self, z, direction=DOWN):
3149 '''
3150 Get layer for given depth.
3152 :param z: depth [m]
3153 :param direction: direction of traversal :py:const:`DOWN` or
3154 :py:const:`UP`.
3156 Returns first layer which touches depth ``z`` (tolerant at boundaries).
3157 '''
3159 for layer in self.layers(direction):
3160 if layer.contains(z):
3161 return layer
3162 else:
3163 raise CakeError('Failed extracting layer at depth z=%s' % z)
3165 def walker(self):
3166 return Walker(self._elements)
3168 def material(self, z, direction=DOWN):
3169 '''
3170 Get material at given depth.
3172 :param z: depth [m]
3173 :param direction: direction of traversal :py:const:`DOWN` or
3174 :py:const:`UP`
3175 :returns: object of type :py:class:`Material`
3177 If given depth ``z`` happens to be at an interface, the material of the
3178 first layer with respect to the the traversal ordering is returned.
3179 '''
3181 lyr = self.layer(z, direction)
3182 return lyr.material(z)
3184 def discontinuities(self):
3185 '''
3186 Iterate over all discontinuities of the model.'''
3188 return (el for el in self._elements if isinstance(el, Discontinuity))
3190 def discontinuity(self, name_or_z):
3191 '''
3192 Get discontinuity by name or depth.
3194 :param name_or_z: name of discontinuity or depth [m] as float value
3195 '''
3197 if isinstance(name_or_z, float):
3198 candi = sorted(
3199 self.discontinuities(), key=lambda i: abs(i.z-name_or_z))
3200 else:
3201 candi = [i for i in self.discontinuities() if i.name == name_or_z]
3203 if not candi:
3204 raise DiscontinuityNotFound(name_or_z)
3206 return candi[0]
3208 def adapt_phase(self, phase):
3209 '''
3210 Adapt a phase definition for use with this model.
3212 This returns a copy of the phase definition, where named
3213 discontinuities are replaced with the actual depth of these, as defined
3214 in the model.
3215 '''
3217 phase = phase.copy()
3218 for knee in phase.knees():
3219 if knee.depth != 'surface':
3220 knee.depth = self.discontinuity(knee.depth).z
3221 for leg in phase.legs():
3222 if leg.depthmax is not None and isinstance(leg.depthmax, str):
3223 leg.depthmax = self.discontinuity(leg.depthmax).z
3225 return phase
3227 def path(self, p, phase, layer_start, layer_stop):
3228 '''
3229 Get ray path for given combination of ray parameter, phase definition,
3230 source and receiver layers.
3232 :param p: ray parameter (spherical) [s/rad]
3233 :param phase: phase definition (:py:class:`PhaseDef` object)
3234 :param layer_start: layer with source
3235 :param layer_stop: layer with receiver
3236 :returns: :py:class:`RayPath` object
3238 If it is not possible to find a solution, an exception of type
3239 :py:exc:`NotPhaseConform`, :py:exc:`MinDepthReached`,
3240 :py:exc:`MaxDepthReached`, :py:exc:`CannotPropagate`,
3241 :py:exc:`BottomReached` or :py:exc:`SurfaceReached` is raised.
3242 '''
3244 phase = self.adapt_phase(phase)
3245 knees = phase.knees()
3246 legs = phase.legs()
3247 next_knee = next_or_none(knees)
3248 leg = next_or_none(legs)
3249 assert leg is not None
3251 direction = leg.departure
3252 direction_stop = phase.direction_stop()
3253 mode = leg.mode
3254 mode_stop = phase.last_leg().mode
3256 walker = self.walker()
3257 walker.goto_layer(layer_start)
3258 current = walker.current()
3260 ttop, tbot = current.tests(p, mode)
3261 if not ttop and not tbot:
3262 raise CannotPropagate(direction, current.ilayer)
3264 if (direction == DOWN and not ttop) or (direction == UP and not tbot):
3265 direction = -direction
3267 path = RayPath(phase)
3268 trapdetect = set()
3269 while True:
3270 at_layer = isinstance(current, Layer)
3271 at_discontinuity = isinstance(current, Discontinuity)
3273 # detect trapped wave
3274 k = (id(next_knee), id(current), direction, mode)
3275 if k in trapdetect:
3276 raise Trapped()
3278 trapdetect.add(k)
3280 if at_discontinuity:
3281 oldmode, olddirection = mode, direction
3282 headwave = False
3283 if next_knee is not None and next_knee.matches(
3284 current, mode, direction):
3286 headwave = next_knee.headwave
3287 direction = next_knee.out_direction()
3288 mode = next_knee.out_mode
3289 next_knee = next_or_none(knees)
3290 leg = next(legs)
3292 else: # implicit reflection/transmission
3293 direction = current.propagate(p, mode, direction)
3295 if headwave:
3296 path.set_is_headwave(True)
3298 path.append(Kink(
3299 olddirection, olddirection, oldmode, oldmode, current))
3301 path.append(HeadwaveStraight(
3302 olddirection, direction, oldmode, current))
3304 path.append(Kink(
3305 olddirection, direction, oldmode, mode, current))
3307 else:
3308 path.append(Kink(
3309 olddirection, direction, oldmode, mode, current))
3311 if at_layer:
3312 direction_in = direction
3313 direction = current.propagate(p, mode, direction_in)
3315 zturn = None
3316 if direction_in != direction:
3317 zturn = current.zturn(p, mode)
3319 zmin, zmax = leg.depthmin, leg.depthmax
3320 if zmin is not None or zmax is not None:
3321 if direction_in != direction:
3322 if zmin is not None and zturn <= zmin:
3323 raise MinDepthReached()
3324 if zmax is not None and zturn >= zmax:
3325 raise MaxDepthReached()
3326 else:
3327 if zmin is not None and current.ztop <= zmin:
3328 raise MinDepthReached()
3329 if zmax is not None and current.zbot >= zmax:
3330 raise MaxDepthReached()
3332 path.append(Straight(direction_in, direction, mode, current))
3334 if next_knee is None and mode == mode_stop and \
3335 current is layer_stop:
3337 if zturn is None:
3338 if direction == direction_stop:
3339 break
3340 else:
3341 break
3343 walker.go(direction)
3344 current = walker.current()
3346 return path
3348 def gather_paths(self, phases=PhaseDef('P'), zstart=0.0, zstop=0.0):
3349 '''
3350 Get all possible ray paths for given source and receiver depths for one
3351 or more phase definitions.
3353 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3354 Comma-separated strings and lists of such strings are also accepted
3355 and are converted to :py:class:`PhaseDef` objects for convenience.
3356 :param zstart: source depth [m]
3357 :param zstop: receiver depth [m]
3358 :returns: a list of :py:class:`RayPath` objects
3360 Results of this method are cached internally. Cached results are
3361 returned, when a given combination of source layer, receiver layer and
3362 phase definition has been used before.
3363 '''
3365 eps = 1e-7 # num.finfo(float).eps * 1000.
3367 phases = to_phase_defs(phases)
3369 paths = []
3370 for phase in phases:
3372 layer_start = self.layer(zstart, -phase.direction_start())
3373 layer_stop = self.layer(zstop, phase.direction_stop())
3375 pathcachekey = (phase.definition(), layer_start, layer_stop)
3377 if pathcachekey in self._pathcache:
3378 phase_paths = self._pathcache[pathcachekey]
3379 else:
3380 hwknee = phase.headwave_knee()
3381 if hwknee:
3382 name_or_z = hwknee.depth
3383 interface = self.discontinuity(name_or_z)
3384 mode = hwknee.in_mode
3385 in_direction = hwknee.direction
3387 pabove, pbelow = interface.critical_ps(mode)
3389 p = min_not_none(pabove, pbelow)
3391 # diffracted wave:
3392 if in_direction == DOWN and (
3393 pbelow is None or pbelow >= pabove):
3395 p *= (1.0 - eps)
3397 path = self.path(p, phase, layer_start, layer_stop)
3398 path.set_prange(p, p, 1.)
3400 phase_paths = [path]
3402 else:
3403 try:
3404 pmax_start = max([
3405 radius(z)/layer_start.v(phase.first_leg().mode, z)
3406 for z in (layer_start.ztop, layer_start.zbot)])
3408 pmax_stop = max([
3409 radius(z)/layer_stop.v(phase.last_leg().mode, z)
3410 for z in (layer_stop.ztop, layer_stop.zbot)])
3412 pmax = min(pmax_start, pmax_stop)
3414 pedges = [0.]
3415 for layer in self.layers():
3416 for z in (layer.ztop, layer.zbot):
3417 for mode in (P, S):
3418 for eps2 in [eps]:
3419 v = layer.v(mode, z)
3420 if v != 0.0:
3421 p = radius(z)/v
3422 if p <= pmax:
3423 pedges.append(p*(1.0-eps2))
3424 pedges.append(p)
3425 pedges.append(p*(1.0+eps2))
3427 pedges = num.unique(sorted(pedges))
3429 phase_paths = {}
3430 cached = {}
3431 counter = [0]
3433 def p_to_path(p):
3434 if p in cached:
3435 return cached[p]
3437 try:
3438 counter[0] += 1
3439 path = self.path(
3440 p, phase, layer_start, layer_stop)
3442 if path not in phase_paths:
3443 phase_paths[path] = []
3445 phase_paths[path].append(p)
3447 except PathFailed:
3448 path = None
3450 cached[p] = path
3451 return path
3453 def recurse(pmin, pmax, i=0):
3454 if i > self._pdepth:
3455 return
3456 path1 = p_to_path(pmin)
3457 path2 = p_to_path(pmax)
3458 if path1 is None and path2 is None and i > 0:
3459 return
3460 if path1 is None or path2 is None or \
3461 hash(path1) != hash(path2):
3463 recurse(pmin, (pmin+pmax)/2., i+1)
3464 recurse((pmin+pmax)/2., pmax, i+1)
3466 for (pl, ph) in zip(pedges[:-1], pedges[1:]):
3467 recurse(pl, ph)
3469 for path, ps in phase_paths.items():
3470 path.set_prange(
3471 min(ps), max(ps), pmax/(self._np-1))
3473 phase_paths = list(phase_paths.keys())
3475 except ZeroDivisionError:
3476 phase_paths = []
3478 self._pathcache[pathcachekey] = phase_paths
3480 paths.extend(phase_paths)
3482 paths.sort(key=lambda x: x.pmin())
3483 return paths
3485 def arrivals(
3486 self,
3487 distances=[],
3488 phases=PhaseDef('P'),
3489 zstart=0.0,
3490 zstop=0.0,
3491 refine=True):
3493 '''
3494 Compute rays and traveltimes for given distances.
3496 :param distances: list or array of distances [deg]
3497 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3498 Comma-separated strings and lists of such strings are also accepted
3499 and are converted to :py:class:`PhaseDef` objects for convenience.
3500 :param zstart: source depth [m]
3501 :param zstop: receiver depth [m]
3502 :param refine: bool flag, whether to use bisectioning to improve
3503 (p, x, t) estimated from interpolation
3504 :returns: a list of :py:class:`Ray` objects, sorted by
3505 (distance, arrival time)
3506 '''
3508 distances = num.asarray(distances, dtype=float)
3510 arrivals = []
3511 for path in self.gather_paths(phases, zstart=zstart, zstop=zstop):
3513 endgaps = path.endgaps(zstart, zstop)
3514 for x, p, t, draft_pxt in path.interpolate_x2pt_linear(
3515 distances, endgaps):
3517 arrivals.append(Ray(path, p, x, t, endgaps, draft_pxt))
3519 if refine:
3520 refined = []
3521 for ray in arrivals:
3523 if ray.path._is_headwave:
3524 refined.append(ray)
3526 try:
3527 ray.refine()
3528 refined.append(ray)
3530 except RefineFailed:
3531 pass
3533 arrivals = refined
3535 arrivals.sort(key=lambda x: (x.x, x.t))
3536 return arrivals
3538 @classmethod
3539 def from_scanlines(cls, producer):
3540 '''
3541 Create layer cake model from sequence of materials at depths.
3543 :param producer: iterable yielding (depth, material, name) tuples
3545 Creates a new :py:class:`LayeredModel` object and uses its
3546 :py:meth:`append` method to add layers and discontinuities as needed.
3547 '''
3549 self = cls()
3550 for z, material, name in producer:
3552 if not self._elements:
3553 self.append(Surface(z, material))
3554 else:
3555 element = self._elements[-1]
3556 if self.zeq(element.zbot, z):
3557 assert isinstance(element, Layer)
3558 self.append(
3559 Interface(z, element.mbot, material, name=name))
3561 else:
3562 if isinstance(element, Discontinuity):
3563 ztop = element.z
3564 mtop = element.mbelow
3565 elif isinstance(element, Layer):
3566 ztop = element.zbot
3567 mtop = element.mbot
3569 if mtop == material:
3570 layer = HomogeneousLayer(
3571 ztop, z, material, name=name)
3572 else:
3573 layer = GradientLayer(
3574 ztop, z, mtop, material, name=name)
3576 self.append(layer)
3578 return self
3580 def to_scanlines(self, get_burgers=False):
3581 def fmt(z, m):
3582 if not m._has_default_burgers() or get_burgers:
3583 return (z, m.vp, m.vs, m.rho, m.qp, m.qs,
3584 m.burger_eta1, m.burger_eta2, m.burger_valpha)
3585 return (z, m.vp, m.vs, m.rho, m.qp, m.qs)
3587 last = None
3588 lines = []
3589 for element in self.elements():
3590 if isinstance(element, Layer):
3591 if not isinstance(last, Layer):
3592 lines.append(fmt(element.ztop, element.mtop))
3594 lines.append(fmt(element.zbot, element.mbot))
3596 last = element
3598 if not isinstance(last, Layer):
3599 lines.append(fmt(last.z, last.mbelow))
3601 return lines
3603 def iter_material_parameter(self, get):
3604 assert get in ('vp', 'vs', 'rho', 'qp', 'qs', 'z')
3605 if get == 'z':
3606 for layer in self.layers():
3607 yield layer.ztop
3608 yield layer.zbot
3609 else:
3610 getter = operator.attrgetter(get)
3611 for layer in self.layers():
3612 yield getter(layer.mtop)
3613 yield getter(layer.mbot)
3615 def profile(self, get):
3616 '''
3617 Get parameter profile along depth of the earthmodel.
3619 :param get: property to be queried (
3620 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3621 :type get: string
3622 '''
3624 return num.array(list(self.iter_material_parameter(get)))
3626 def min(self, get='vp'):
3627 '''
3628 Find minimum value of a material property or depth.
3630 :param get: property to be queried (
3631 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3632 '''
3634 return min(self.iter_material_parameter(get))
3636 def max(self, get='vp'):
3637 '''
3638 Find maximum value of a material property or depth.
3640 :param get: property to be queried (
3641 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, ``'qs'``, or ``'z'``)
3642 '''
3644 return max(self.iter_material_parameter(get))
3646 def simplify_layers(self, layers, max_rel_error=0.001):
3647 if len(layers) <= 1:
3648 return layers
3650 ztop = layers[0].ztop
3651 zbot = layers[-1].zbot
3652 zorigs = [layer.ztop for layer in layers]
3653 zorigs.append(zbot)
3654 zs = num.linspace(ztop, zbot, 100)
3655 data = []
3656 for z in zs:
3657 if z == ztop:
3658 direction = UP
3659 else:
3660 direction = DOWN
3662 mat = self.material(z, direction)
3663 data.append(mat.astuple())
3665 data = num.array(data, dtype=float)
3666 data_means = num.mean(data, axis=0)
3667 nmax = len(layers) // 2
3668 accept = False
3670 zcut_best = []
3671 for n in range(1, nmax+1):
3672 ncutintervals = 20
3673 zdelta = (zbot-ztop)/ncutintervals
3674 if n == 2:
3675 zcuts = [
3676 [ztop, ztop + i*zdelta, zbot]
3677 for i in range(1, ncutintervals)]
3678 elif n == 3:
3679 zcuts = []
3680 for j in range(1, ncutintervals):
3681 for i in range(j+1, ncutintervals):
3682 zcuts.append(
3683 [ztop, ztop + j*zdelta, ztop + i*zdelta, zbot])
3684 else:
3685 zcuts = []
3686 zcuts.append(num.linspace(ztop, zbot, n+1))
3687 if zcut_best:
3688 zcuts.append(sorted(num.linspace(
3689 ztop, zbot, n).tolist() + zcut_best[1]))
3690 zcuts.append(sorted(num.linspace(
3691 ztop, zbot, n-1).tolist() + zcut_best[2]))
3693 best = None
3694 for icut, zcut in enumerate(zcuts):
3695 rel_par_errors = num.zeros(5)
3696 mpar_nodes = num.zeros((n+1, 5))
3698 for ipar in range(5):
3699 znodes, vnodes, error_rms = util.polylinefit(
3700 zs, data[:, ipar], zcut)
3702 mpar_nodes[:, ipar] = vnodes
3703 if data_means[ipar] == 0.0:
3704 rel_par_errors[ipar] = -1
3705 else:
3706 rel_par_errors[ipar] = error_rms/data_means[ipar]
3708 rel_error = rel_par_errors.max()
3709 if best is None or rel_error < best[0]:
3710 best = (rel_error, zcut, mpar_nodes)
3712 rel_error, zcut, mpar_nodes = best
3714 zcut_best.append(list(zcut))
3715 zcut_best[-1].pop(0)
3716 zcut_best[-1].pop()
3718 if rel_error <= max_rel_error:
3719 accept = True
3720 break
3722 if not accept:
3723 return layers
3725 rel_error, zcut, mpar_nodes = best
3727 material_nodes = []
3728 for i in range(n+1):
3729 material_nodes.append(Material(*mpar_nodes[i, :]))
3731 out_layers = []
3732 for i in range(n):
3733 mtop = material_nodes[i]
3734 mbot = material_nodes[i+1]
3735 ztop = zcut[i]
3736 zbot = zcut[i+1]
3737 if mtop == mbot:
3738 lyr = HomogeneousLayer(ztop, zbot, mtop)
3739 else:
3740 lyr = GradientLayer(ztop, zbot, mtop, mbot)
3742 out_layers.append(lyr)
3743 return out_layers
3745 def simplify(self, max_rel_error=0.001):
3746 '''
3747 Get representation of model with lower resolution.
3749 Returns an approximation of the model. All discontinuities are kept,
3750 but layer stacks with continuous model parameters are represented, if
3751 possible, by a lower number of layers. Piecewise linear functions are
3752 fitted against the original model parameter's piecewise linear
3753 functions. Successively larger numbers of layers are tried, until the
3754 difference to the original model is below ``max_rel_error``. The
3755 difference is measured as the RMS error of the fit normalized by the
3756 mean of the input (i.e. the fitted curves should deviate, on average,
3757 less than 0.1% from the input curves if ``max_rel_error`` = 0.001).
3758 '''
3760 mod_simple = LayeredModel()
3762 glayers = []
3763 for element in self.elements():
3765 if isinstance(element, Discontinuity):
3766 for layer in self.simplify_layers(
3767 glayers, max_rel_error=max_rel_error):
3769 mod_simple.append(layer)
3771 glayers = []
3772 mod_simple.append(element)
3773 else:
3774 glayers.append(element)
3776 for layer in self.simplify_layers(
3777 glayers, max_rel_error=max_rel_error):
3778 mod_simple.append(layer)
3780 return mod_simple
3782 def extract(self, depth_min=None, depth_max=None):
3783 '''
3784 Extract :py:class:`LayeredModel` from :py:class:`LayeredModel`.
3786 :param depth_min: depth of upper cut or name of :py:class:`Interface`
3787 :param depth_max: depth of lower cut or name of :py:class:`Interface`
3789 Interpolates a :py:class:`GradientLayer` at ``depth_min`` and/or
3790 ``depth_max``.
3791 '''
3793 if isinstance(depth_min, (str, newstr)):
3794 depth_min = self.discontinuity(depth_min).z
3796 if isinstance(depth_max, (str, newstr)):
3797 depth_max = self.discontinuity(depth_max).z
3799 mod_extracted = LayeredModel()
3801 for element in self.elements():
3802 element = element.copy()
3803 do_append = False
3804 if (depth_min is None or depth_min <= element.ztop) \
3805 and (depth_max is None or depth_max >= element.zbot):
3806 mod_extracted.append(element)
3807 continue
3809 if depth_min is not None:
3810 if element.ztop < depth_min and depth_min < element.zbot:
3811 _, element = element.split(depth_min)
3812 do_append = True
3814 if depth_max is not None:
3815 if element.zbot > depth_max and depth_max > element.ztop:
3816 element, _ = element.split(depth_max)
3817 do_append = True
3819 if do_append:
3820 mod_extracted.append(element)
3822 return mod_extracted
3824 def replaced_crust(self, crust2_profile=None, crustmod=None):
3825 if crust2_profile is not None:
3826 profile = anything_to_crust2_profile(crust2_profile)
3827 crustmod = LayeredModel.from_scanlines(
3828 from_crust2x2_profile(profile))
3830 newmod = LayeredModel()
3831 for element in crustmod.extract(depth_max='moho').elements():
3832 if element.name != 'moho':
3833 newmod.append(element)
3834 else:
3835 moho1 = element
3837 mod = self.extract(depth_min='moho')
3838 first = True
3839 for element in mod.elements():
3840 if element.name == 'moho':
3841 if element.z <= moho1.z:
3842 mbelow = mod.material(moho1.z, direction=UP)
3843 else:
3844 mbelow = element.mbelow
3846 moho = Interface(moho1.z, moho1.mabove, mbelow, name='moho')
3847 newmod.append(moho)
3848 else:
3849 if first:
3850 if isinstance(element, Layer) and element.zbot > moho.z:
3851 newmod.append(GradientLayer(
3852 moho.z,
3853 element.zbot,
3854 moho.mbelow,
3855 element.mbot,
3856 name=element.name))
3858 first = False
3859 else:
3860 newmod.append(element)
3861 return newmod
3863 def perturb(self, rstate=None, keep_vp_vs=False, **kwargs):
3864 '''
3865 Create a perturbed variant of the earth model.
3867 Randomly change the thickness and material parameters of the earth
3868 model from a uniform distribution.
3870 :param kwargs: Maximum change fraction (e.g. 0.1) of the parameters.
3871 Name the parameter, prefixed by ``p``. Supported parameters are
3872 ``ph, pvp, pvs, prho, pqs, pqp``.
3873 :type kwargs: dict
3874 :param rstate: Random state to draw from, defaults to ``None``
3875 :type rstate: :class:`numpy.random.RandomState`, optional
3876 :param keep_vp_vs: Keep the Vp/Vs ratio, defaults to False
3877 :type keep_vp_vs: bool, optional
3879 :returns: A new, perturbed earth model
3880 :rtype: :class:`~pyrocko.cake.LayeredModel`
3882 .. code-block :: python
3884 perturbed_model = model.perturb(ph=.1, pvp=.05, prho=.1)
3885 '''
3886 _pargs = set(['ph', 'pvp', 'pvs', 'prho', 'pqs', 'pqp'])
3887 earthmod = copy.deepcopy(self)
3889 if rstate is None:
3890 rstate = num.random.RandomState()
3892 layers = earthmod.layers()
3893 discont = earthmod.discontinuities()
3894 prev_layer = None
3896 def get_change_ratios():
3897 values = dict.fromkeys([p[1:] for p in _pargs], 0.)
3899 for param, pval in kwargs.items():
3900 if param not in _pargs:
3901 continue
3902 values[param[1:]] = float(rstate.uniform(-pval, pval, size=1))
3903 return values
3905 # skip Surface
3906 while True:
3907 disc = next(discont)
3908 if isinstance(disc, Surface):
3909 break
3911 while True:
3912 try:
3913 layer = next(layers)
3914 m = layer.material(None)
3915 h = layer.zbot - layer.ztop
3916 except StopIteration:
3917 break
3919 if not isinstance(layer, HomogeneousLayer):
3920 raise NotImplementedError(
3921 'Can only perturbate homogeneous layers!')
3923 changes = get_change_ratios()
3925 # Changing thickness
3926 dh = h * changes['h']
3927 changes['h'] = dh
3929 layer.resize(depth_max=layer.zbot + dh,
3930 depth_min=prev_layer.zbot if prev_layer else None)
3932 try:
3933 disc = next(discont)
3934 disc.change_depth(disc.z + dh)
3935 except StopIteration:
3936 pass
3938 # Setting material parameters
3939 for param, change_ratio in changes.items():
3940 if param == 'h':
3941 continue
3943 value = m.__getattribute__(param)
3944 changes[param] = value * change_ratio
3946 if keep_vp_vs and changes['vp'] != 0.:
3947 changes['vs'] = (m.vp + changes['vp']) / m.vp_vs_ratio() - m.vs
3949 for param, change in changes.items():
3950 if param == 'h':
3951 continue
3952 value = m.__getattribute__(param)
3953 m.__setattr__(param, value + change)
3955 logger.info(
3956 'perturbating earthmodel: {}'.format(
3957 ' '.join(['{param}: {change:{len}.2f}'.format(
3958 param=p, change=c, len=8)
3959 for p, c in changes.items()])))
3961 prev_layer = layer
3963 return earthmod
3965 def require_homogeneous(self):
3966 elements = list(self.elements())
3968 if len(elements) != 2:
3969 raise LayeredModelError(
3970 'Homogeneous model required but found more than one layer in '
3971 'earthmodel.')
3973 if not isinstance(elements[1], HomogeneousLayer):
3974 raise LayeredModelError(
3975 'Homogeneous model required but element #1 is not of type '
3976 'HomogeneousLayer.')
3978 return elements[1].m
3980 def __str__(self):
3981 return '\n'.join(str(element) for element in self._elements)
3984def read_hyposat_model(fn):
3985 '''
3986 Reader for HYPOSAT earth model files.
3988 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
3990 Interface names are translated as follows: ``'MOHO'`` -> ``'moho'``,
3991 ``'CONR'`` -> ``'conrad'``
3992 '''
3994 with open(fn, 'r') as f:
3995 translate = {'MOHO': 'moho', 'CONR': 'conrad'}
3996 lname = None
3997 for iline, line in enumerate(f):
3998 if iline == 0:
3999 continue
4001 z, vp, vs, name = util.unpack_fixed('f10, f10, f10, a4', line)
4002 if not name:
4003 name = None
4004 material = Material(vp*1000., vs*1000.)
4006 tname = translate.get(lname, lname)
4007 yield z*1000., material, tname
4009 lname = name
4012def read_nd_model(fn):
4013 '''
4014 Reader for TauP style '.nd' (named discontinuity) files.
4016 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
4018 Interface names are translated as follows: ``'mantle'`` -> ``'moho'``,
4019 ``'outer-core'`` -> ``'cmb'``, ``'inner-core'`` -> ``'icb'``.
4021 The format has been modified to include Burgers materials parameters in
4022 columns 7 (burger_eta1), 8 (burger_eta2) and 9. eta(3).
4023 '''
4024 with open(fn, 'r') as f:
4025 for x in read_nd_model_fh(f):
4026 yield x
4029def read_nd_model_str(s):
4030 f = StringIO(s)
4031 for x in read_nd_model_fh(f):
4032 yield x
4033 f.close()
4036def read_nd_model_fh(f):
4037 translate = {'mantle': 'moho', 'outer-core': 'cmb', 'inner-core': 'icb'}
4038 name = None
4039 for line in f:
4040 toks = line.split()
4041 if len(toks) == 9 or len(toks) == 6 or len(toks) == 4:
4042 z, vp, vs, rho = [float(x) for x in toks[:4]]
4043 qp, qs = None, None
4044 burgers = None
4045 if len(toks) == 6 or len(toks) == 9:
4046 qp, qs = [float(x) for x in toks[4:6]]
4047 if len(toks) == 9:
4048 burgers = \
4049 [float(x) for x in toks[6:]]
4051 material = Material(
4052 vp*1000., vs*1000., rho*1000., qp, qs,
4053 burgers=burgers)
4055 yield z*1000., material, name
4056 name = None
4057 elif len(toks) == 1:
4058 name = translate.get(toks[0], toks[0])
4060 f.close()
4063def from_crust2x2_profile(profile, depthmantle=50000):
4064 from pyrocko.dataset import crust2x2
4066 default_qp_qs = {
4067 'soft sed.': (50., 50.),
4068 'hard sed.': (200., 200.),
4069 'upper crust': (600., 400.),
4070 }
4072 z = 0.
4073 for i in range(8):
4074 dz, vp, vs, rho = profile.get_layer(i)
4075 name = crust2x2.Crust2Profile.layer_names[i]
4076 if name in default_qp_qs:
4077 qp, qs = default_qp_qs[name]
4078 else:
4079 qp, qs = None, None
4081 material = Material(vp, vs, rho, qp, qs)
4082 iname = None
4083 if i == 7:
4084 iname = 'moho'
4085 if dz != 0.0:
4086 yield z, material, iname
4087 if i != 7:
4088 yield z+dz, material, name
4089 else:
4090 yield z+depthmantle, material, name
4092 z += dz
4095def write_nd_model_fh(mod, fh):
4096 def fmt(z, mat):
4097 rstr = ' '.join(
4098 util.gform(x, 4)
4099 for x in (
4100 z/1000.,
4101 mat.vp/1000.,
4102 mat.vs/1000.,
4103 mat.rho/1000.,
4104 mat.qp, mat.qs))
4105 if not mat._has_default_burgers():
4106 rstr += ' '.join(
4107 util.gform(x, 4)
4108 for x in (
4109 mat.burger_eta1,
4110 mat.burger_eta2,
4111 mat.burger_valpha))
4112 return rstr.rstrip() + '\n'
4114 translate = {
4115 'moho': 'mantle',
4116 'cmb': 'outer-core',
4117 'icb': 'inner-core'}
4119 last = None
4120 for element in mod.elements():
4121 if isinstance(element, Interface):
4122 if element.name is not None:
4123 n = translate.get(element.name, element.name)
4124 fh.write('%s\n' % n)
4126 elif isinstance(element, Layer):
4127 if not isinstance(last, Layer):
4128 fh.write(fmt(element.ztop, element.mtop))
4130 fh.write(fmt(element.zbot, element.mbot))
4132 last = element
4134 if not isinstance(last, Layer):
4135 fh.write(fmt(last.z, last.mbelow))
4138def write_nd_model_str(mod):
4139 f = StringIO()
4140 write_nd_model_fh(mod, f)
4141 return f.getvalue()
4144def write_nd_model(mod, fn):
4145 with open(fn, 'w') as f:
4146 write_nd_model_fh(mod, f)
4149def builtin_models():
4150 return sorted([
4151 os.path.splitext(os.path.basename(x))[0]
4152 for x in glob.glob(builtin_model_filename('*'))])
4155def builtin_model_filename(modelname):
4156 return util.data_file(os.path.join('earthmodels', modelname+'.nd'))
4159def load_model(fn='ak135-f-continental.m', format='nd', crust2_profile=None):
4160 '''
4161 Load layered earth model from file.
4163 :param fn: filename
4164 :param format: format
4165 :param crust2_profile: ``(lat, lon)`` or
4166 :py:class:`pyrocko.crust2x2.Crust2Profile` object, merge model with
4167 crustal profile. If ``fn`` is forced to be ``None`` only the converted
4168 CRUST2.0 profile is returned.
4169 :returns: object of type :py:class:`LayeredModel`
4171 The following formats are currently supported:
4173 ============== ===========================================================
4174 format description
4175 ============== ===========================================================
4176 ``'nd'`` 'named discontinuity' format used by the TauP programs
4177 ``'hyposat'`` format used by the HYPOSAT location program
4178 ============== ===========================================================
4180 The naming of interfaces is translated from the file format's native naming
4181 to Cake's own convention (See :py:func:`read_nd_model` and
4182 :py:func:`read_hyposat_model` for details). Cake likes the following
4183 internal names: ``'conrad'``, ``'moho'``, ``'cmb'`` (core-mantle boundary),
4184 ``'icb'`` (inner core boundary).
4185 '''
4187 if fn is not None:
4188 if format == 'nd':
4189 if not os.path.exists(fn) and fn in builtin_models():
4190 fn = builtin_model_filename(fn)
4191 reader = read_nd_model(fn)
4192 elif format == 'hyposat':
4193 reader = read_hyposat_model(fn)
4194 else:
4195 assert False, 'unsupported model format'
4197 mod = LayeredModel.from_scanlines(reader)
4198 if crust2_profile is not None:
4199 return mod.replaced_crust(crust2_profile)
4201 return mod
4203 else:
4204 assert crust2_profile is not None
4205 profile = anything_to_crust2_profile(crust2_profile)
4206 return LayeredModel.from_scanlines(
4207 from_crust2x2_profile(profile))
4210def castagna_vs_to_vp(vs):
4211 '''
4212 Calculate vp from vs using Castagna's relation.
4214 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4215 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4216 1985]
4218 vp = 1.16 * vs + 1360 [m/s]
4220 :param vs: S-wave velocity [m/s]
4221 :returns: P-wave velocity [m/s]
4222 '''
4224 return vs*1.16 + 1360.0
4227def castagna_vp_to_vs(vp):
4228 '''
4229 Calculate vp from vs using Castagna's relation.
4231 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4232 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4233 1985]
4235 vp = 1.16 * vs + 1360 [m/s]
4237 :param vp: P-wave velocity [m/s]
4238 :returns: S-wave velocity [m/s]
4239 '''
4241 return (vp - 1360.0) / 1.16
4244def evenize(x, y, minsize=10):
4245 if x.size < minsize:
4246 return x
4247 ry = (y.max()-y.min())
4248 if ry == 0:
4249 return x
4250 dx = (x[1:] - x[:-1])/(x.max()-x.min())
4251 dy = (y[1:] + y[:-1])/ry
4253 s = num.zeros(x.size)
4254 s[1:] = num.cumsum(num.sqrt(dy**2 + dx**2))
4255 s2 = num.linspace(0, s[-1], x.size)
4256 x2 = num.interp(s2, s, x)
4257 x2[0] = x[0]
4258 x2[-1] = x[-1]
4259 return x2
4262def filled(v, *args, **kwargs):
4263 '''
4264 Create NumPy array filled with given value.
4266 This works like :py:func:`numpy.ones` but initializes the array with ``v``
4267 instead of ones.
4268 '''
4269 x = num.empty(*args, **kwargs)
4270 x.fill(v)
4271 return x
4274def next_or_none(i):
4275 try:
4276 return next(i)
4277 except StopIteration:
4278 return None
4281def reci_or_none(x):
4282 try:
4283 return 1./x
4284 except ZeroDivisionError:
4285 return None
4288def mult_or_none(a, b):
4289 if a is None or b is None:
4290 return None
4291 return a*b
4294def min_not_none(a, b):
4295 if a is None:
4296 return b
4297 if b is None:
4298 return a
4299 return min(a, b)
4302def xytups(xx, yy):
4303 d = []
4304 for x, y in zip(xx, yy):
4305 if num.isfinite(y):
4306 d.append((x, y))
4307 return d
4310def interp(x, xp, fp, monoton):
4311 if monoton == 1:
4312 return xytups(
4313 x, num.interp(x, xp, fp, left=num.nan, right=num.nan))
4314 elif monoton == -1:
4315 return xytups(
4316 x, num.interp(x, xp[::-1], fp[::-1], left=num.nan, right=num.nan))
4317 else:
4318 fs = []
4319 for xv in x:
4320 indices = num.where(num.logical_or(
4321 num.logical_and(xp[:-1] >= xv, xv > xp[1:]),
4322 num.logical_and(xp[:-1] <= xv, xv < xp[1:])))[0]
4324 for i in indices:
4325 xr = (xv - xp[i])/(xp[i+1]-xp[i])
4326 fv = xr*fp[i] + (1.-xr)*fp[i+1]
4327 fs.append((xv, fv))
4329 return fs
4332def float_or_none(x):
4333 if x is not None:
4334 return float(x)
4337def parstore_float(thelocals, obj, *args):
4338 for k, v in thelocals.items():
4339 if k != 'self' and (not args or k in args):
4340 setattr(obj, k, float_or_none(v))