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 functools import reduce
48import os
49import logging
50import copy
51import math
52import cmath
53import operator
54try:
55 from StringIO import StringIO
56except ImportError:
57 from io import StringIO
59import glob
60import numpy as num
61from scipy.optimize import bisect, brentq
63from . import util, config
66logger = logging.getLogger('cake')
68ZEPS = 0.01
69P = 1
70S = 2
71DOWN = 4
72UP = -4
74DEFAULT_BURGERS = (0., 0., 1.)
76earthradius = config.config().earthradius
78r2d = 180./math.pi
79d2r = 1./r2d
80km = 1000.
81d2m = d2r*earthradius
82m2d = 1./d2m
83sprad2spm = 1.0/(r2d*d2m)
84sprad2spkm = 1.0/(r2d*d2m/km)
85spm2sprad = 1.0/sprad2spm
86spkm2sprad = 1.0/sprad2spkm
89class CakeError(Exception):
90 pass
93class InvalidArguments(CakeError):
94 pass
97class Material(object):
98 '''
99 Isotropic elastic material.
101 :param vp: P-wave velocity [m/s]
102 :param vs: S-wave velocity [m/s]
103 :param rho: density [kg/m^3]
104 :param qp: P-wave attenuation Qp
105 :param qs: S-wave attenuation Qs
106 :param poisson: Poisson ratio
107 :param lame: tuple with Lame parameter `lambda` and `shear modulus` [Pa]
108 :param qk: bulk attenuation Qk
109 :param qmu: shear attenuation Qmu
111 :param burgers: Burgers rheology paramerters as `tuple`.
112 `transient viscosity` [Pa], <= 0 means infinite value,
113 `steady-state viscosity` [Pa] and `alpha`, the ratio between the
114 effective and unreleaxed shear modulus, mu1/(mu1 + mu2).
115 :type burgers: tuple
117 If no velocities and no lame parameters are given, standard crustal values
118 of vp = 5800 m/s and vs = 3200 m/s are used. If no Q values are given,
119 standard crustal values of qp = 1456 and qs = 600 are used. If no Burgers
120 material parameters are given, transient and steady-state viscosities are
121 0 and alpha=1.
123 Everything is in SI units (m/s, Pa, kg/m^3) unless explicitly stated.
125 The main material properties are considered independant and are accessible
126 as attributes (it is allowed to assign to these):
128 .. py:attribute:: vp, vs, rho, qp, qs
130 Other material properties are considered dependant and can be queried by
131 instance methods.
132 '''
134 def __init__(
135 self, vp=None, vs=None, rho=2600., qp=None, qs=None, poisson=None,
136 lame=None, qk=None, qmu=None, burgers=None):
138 parstore_float(locals(), self, 'vp', 'vs', 'rho', 'qp', 'qs')
140 if vp is not None and vs is not None:
141 if poisson is not None or lame is not None:
142 raise InvalidArguments(
143 'If vp and vs are given, poisson ratio and lame paramters '
144 'should not be given.')
146 elif vp is None and vs is None and lame is None:
147 self.vp = 5800.
148 if poisson is None:
149 poisson = 0.25
150 self.vs = self.vp / math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson))
152 elif vp is None and vs is None and lame is not None:
153 if poisson is not None:
154 raise InvalidArguments(
155 'Poisson ratio should not be given, when lame parameters '
156 'are given.')
158 lam, mu = float(lame[0]), float(lame[1])
159 self.vp = math.sqrt((lam + 2.0*mu)/rho)
160 self.vs = math.sqrt(mu/rho)
162 elif vp is not None and vs is None:
163 if poisson is None:
164 poisson = 0.25
166 if lame is not None:
167 raise InvalidArguments(
168 'If vp is given, Lame parameters should not be given.')
170 poisson = float(poisson)
171 self.vs = vp / math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson))
173 elif vp is None and vs is not None:
174 if poisson is None:
175 poisson = 0.25
176 if lame is not None:
177 raise InvalidArguments(
178 'If vs is given, Lame parameters should not be given.')
180 poisson = float(poisson)
181 self.vp = vs * math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson))
183 else:
184 raise InvalidArguments(
185 'Invalid combination of input parameters in material '
186 'definition.')
188 if qp is not None or qs is not None:
189 if not (qk is None and qmu is None):
190 raise InvalidArguments(
191 'If qp or qs are given, qk and qmu should not be given.')
193 if qp is None:
194 if self.vs != 0.0:
195 s = (4.0/3.0)*(self.vs/self.vp)**2
196 self.qp = self.qs / s
197 else:
198 self.qp = 1456.
200 if qs is None:
201 if self.vs != 0.0:
202 s = (4.0/3.0)*(self.vs/self.vp)**2
203 self.qs = self.qp * s
204 else:
205 self.vs = 600.
207 elif qp is None and qs is None and qk is None and qmu is None:
208 if self.vs == 0.:
209 self.qs = 0.
210 self.qp = 5782e4
211 else:
212 self.qs = 600.
213 s = (4.0/3.0)*(self.vs/self.vp)**2
214 self.qp = self.qs/s
216 elif qp is None and qs is None and qk is not None and qmu is not None:
217 s = (4.0/3.0)*(self.vs/self.vp)**2
218 if qmu == 0. and self.vs == 0.:
219 self.qp = qk
220 else:
221 if num.isinf(qk):
222 self.qp = qmu/s
223 else:
224 self.qp = 1.0 / (s/qmu + (1.0-s)/qk)
225 self.qs = qmu
226 else:
227 raise InvalidArguments(
228 'Invalid combination of input parameters in material '
229 'definition.')
231 if burgers is None:
232 burgers = DEFAULT_BURGERS
234 self.burger_eta1 = burgers[0]
235 self.burger_eta2 = burgers[1]
236 self.burger_valpha = burgers[2]
238 def astuple(self):
239 '''
240 Get independant material properties as a tuple.
242 Returns a tuple with ``(vp, vs, rho, qp, qs)``.
243 '''
244 return self.vp, self.vs, self.rho, self.qp, self.qs
246 def __eq__(self, other):
247 return self.astuple() == other.astuple()
249 def lame(self):
250 '''
251 Get Lame's parameter lambda and shear modulus.
252 '''
253 mu = self.vs**2 * self.rho
254 lam = self.vp**2 * self.rho - 2.0*mu
255 return lam, mu
257 def lame_lambda(self):
258 '''
259 Get Lame's parameter lambda.
261 Returned units are [Pa].
262 '''
263 lam, _ = self.lame()
264 return lam
266 def shear_modulus(self):
267 '''
268 Get shear modulus.
270 Returned units are [Pa].
271 '''
272 return self.vs**2 * self.rho
274 def poisson(self):
275 '''
276 Get Poisson's ratio.
277 '''
278 lam, mu = self.lame()
279 return lam / (2.0*(lam+mu))
281 def bulk(self):
282 '''
283 Get bulk modulus.
284 '''
285 lam, mu = self.lame()
286 return lam + 2.0*mu/3.0
288 def youngs(self):
289 '''
290 Get Young's modulus.
291 '''
292 lam, mu = self.lame()
293 return mu * (3.0*lam + 2.0*mu) / (lam+mu)
295 def vp_vs_ratio(self):
296 '''
297 Get vp/vs ratio.
298 '''
299 return self.vp/self.vs
301 def qmu(self):
302 '''
303 Get shear attenuation coefficient Qmu.
304 '''
305 return self.qs
307 def qk(self):
308 '''
309 Get bulk attenuation coefficient Qk.
310 '''
311 if self.vs == 0. and self.qs == 0.:
312 return self.qp
313 else:
314 s = (4.0/3.0)*(self.vs/self.vp)**2
315 denom = (1/self.qp - s/self.qs)
316 if denom <= 0.0:
317 return num.inf
318 else:
319 return (1.-s)/(1.0/self.qp - s/self.qs)
321 def burgers(self):
322 '''
323 Get Burger parameters.
324 '''
325 return self.burger_eta1, self.burger_eta2, self.burger_valpha
327 def _rayleigh_equation(self, cr):
328 cr_a = (cr/self.vp)**2
329 cr_b = (cr/self.vs)**2
330 if cr_a > 1.0 or cr_b > 1.0:
331 return None
333 return (2.0-cr_b)**2 - 4.0 * math.sqrt(1.0-cr_a) * math.sqrt(1.0-cr_b)
335 def rayleigh(self):
336 '''
337 Get Rayleigh velocity assuming a homogenous halfspace.
339 Returned units are [m/s].
340 '''
341 return bisect(self._rayleigh_equation, 0.001*self.vs, self.vs)
343 def _has_default_burgers(self):
344 if self.burger_eta1 == DEFAULT_BURGERS[0] and \
345 self.burger_eta2 == DEFAULT_BURGERS[1] and \
346 self.burger_valpha == DEFAULT_BURGERS[2]:
347 return True
348 return False
350 def describe(self):
351 '''
352 Get a readable listing of the material properties.
353 '''
354 template = '''
355P wave velocity [km/s] : %12g
356S wave velocity [km/s] : %12g
357P/S wave vel. ratio : %12g
358Lame lambda [GPa] : %12g
359Lame shear modulus [GPa] : %12g
360Poisson ratio : %12g
361Bulk modulus [GPa] : %12g
362Young's modulus [GPa] : %12g
363Rayleigh wave vel. [km/s] : %12g
364Density [g/cm**3] : %12g
365Qp P-wave attenuation : %12g
366Qs S-wave attenuation (Qmu) : %12g
367Qk bulk attenuation : %12g
368transient viscos., eta1 [GPa] : %12g
369st.-state viscos., eta2 [GPa] : %12g
370relaxation: valpha : %12g
371'''.strip()
373 return template % (
374 self.vp/km,
375 self.vs/km,
376 self.vp/self.vs,
377 self.lame_lambda()*1e-9,
378 self.shear_modulus()*1e-9,
379 self.poisson(),
380 self.bulk()*1e-9,
381 self.youngs()*1e-9,
382 self.rayleigh()/km,
383 self.rho/km,
384 self.qp,
385 self.qs,
386 self.qk(),
387 self.burger_eta1*1e-9,
388 self.burger_eta2*1e-9,
389 self.burger_valpha)
391 def __str__(self):
392 vp, vs, rho, qp, qs = self.astuple()
393 return '%10g km/s %10g km/s %10g g/cm^3 %10g %10g' % (
394 vp/km, vs/km, rho/km, qp, qs)
396 def __repr__(self):
397 return 'Material(vp=%s, vs=%s, rho=%s, qp=%s, qs=%s)' % \
398 tuple(repr(x) for x in (
399 self.vp, self.vs, self.rho, self.qp, self.qs))
402class Leg(object):
403 '''
404 Represents a continuous piece of wave propagation in a phase definition.
406 **Attributes:**
408 To be considered as read-only.
410 .. py:attribute:: departure
412 One of the constants :py:const:`UP` or :py:const:`DOWN` indicating
413 upward or downward departure.
415 .. py:attribute:: mode
417 One of the constants :py:const:`P` or :py:const:`S`, indicating the
418 propagation mode.
420 .. py:attribute:: depthmin
422 ``None``, a number (a depth in [m]) or a string (an interface name),
423 minimum depth.
425 .. py:attribute:: depthmax
427 ``None``, a number (a depth in [m]) or a string (an interface name),
428 maximum depth.
430 '''
432 def __init__(self, departure=None, mode=None):
433 self.departure = departure
434 self.mode = mode
435 self.depthmin = None
436 self.depthmax = None
438 def set_depthmin(self, depthmin):
439 self.depthmin = depthmin
441 def set_depthmax(self, depthmax):
442 self.depthmax = depthmax
444 def __str__(self):
445 def sd(d):
446 if isinstance(d, float):
447 return '%g km' % (d/km)
448 else:
449 return 'interface %s' % d
451 s = '%s mode propagation, departing %s' % (
452 smode(self.mode).upper(), {
453 UP: 'upward', DOWN: 'downward'}[self.departure])
455 sc = []
456 if self.depthmax is not None:
457 sc.append('deeper than %s' % sd(self.depthmax))
458 if self.depthmin is not None:
459 sc.append('shallower than %s' % sd(self.depthmin))
461 if sc:
462 s = s + ' (may not propagate %s)' % ' or '.join(sc)
464 return s
467class InvalidKneeDef(CakeError):
468 pass
471class Knee(object):
472 '''
473 Represents a change in wave propagation within a :py:class:`PhaseDef`.
475 **Attributes:**
477 To be considered as read-only.
479 .. py:attribute:: depth
481 Depth at which the conversion/reflection should happen. this can be
482 a string or a number.
484 .. py:attribute:: direction
486 One of the constants :py:const:`UP` or :py:const:`DOWN` to indicate
487 the incoming direction.
489 .. py:attribute:: in_mode
491 One of the constants :py:const:`P` or :py:const:`S` to indicate the
492 type of mode of the incoming wave.
494 .. py:attribute:: out_mode
496 One of the constants :py:const:`P` or :py:const:`S` to indicate the
497 type of mode of the outgoing wave.
499 .. py:attribute:: conversion
501 Boolean, whether there is a mode conversion involved.
503 .. py:attribute:: reflection
505 Boolean, whether there is a reflection involved.
507 .. py:attribute:: headwave
509 Boolean, whether there is headwave propagation involved.
511 '''
513 defaults = dict(
514 depth='surface',
515 direction=UP,
516 conversion=True,
517 reflection=False,
518 headwave=False,
519 in_setup_state=True)
521 defaults_surface = dict(
522 depth='surface',
523 direction=UP,
524 conversion=False,
525 reflection=True,
526 headwave=False,
527 in_setup_state=True)
529 def __init__(self, *args):
530 if args:
531 (self.depth, self.direction, self.reflection, self.in_mode,
532 self.out_mode) = args
534 self.conversion = self.in_mode != self.out_mode
535 self.in_setup_state = False
537 def default(self, k):
538 depth = self.__dict__.get('depth', 'surface')
539 if depth == 'surface':
540 return Knee.defaults_surface[k]
541 else:
542 return Knee.defaults[k]
544 def __setattr__(self, k, v):
545 if self.in_setup_state and k in self.__dict__:
546 raise InvalidKneeDef('%s has already been set' % k)
547 else:
548 self.__dict__[k] = v
550 def __getattr__(self, k):
551 if k.startswith('__'):
552 raise AttributeError(k)
554 if k not in self.__dict__:
555 return self.default(k)
557 def set_modes(self, in_leg, out_leg):
559 if out_leg.departure == UP and (
560 (self.direction == UP) == self.reflection):
562 raise InvalidKneeDef(
563 'cannot enter %s from %s and emit ray upwards' % (
564 ['conversion', 'reflection'][self.reflection],
565 {UP: 'below', DOWN: 'above'}[self.direction]))
567 if out_leg.departure == DOWN and (
568 (self.direction == DOWN) == self.reflection):
570 raise InvalidKneeDef(
571 'cannot enter %s from %s and emit ray downwards' % (
572 ['conversion', 'reflection'][self.reflection],
573 {UP: 'below', DOWN: 'above'}[self.direction]))
575 self.in_mode = in_leg.mode
576 self.out_mode = out_leg.mode
578 def at_surface(self):
579 return self.depth == 'surface'
581 def matches(self, discontinuity, mode, direction):
582 '''
583 Check whether it is relevant to a given combination of interface,
584 propagation mode, and direction.
585 '''
587 if isinstance(self.depth, float):
588 if abs(self.depth - discontinuity.z) > ZEPS:
589 return False
590 else:
591 if discontinuity.name != self.depth:
592 return False
594 return self.direction == direction and self.in_mode == mode
596 def out_direction(self):
597 '''
598 Get outgoing direction.
600 Returns one of the constants :py:const:`UP` or :py:const:`DOWN`.
601 '''
603 if self.reflection:
604 return - self.direction
605 else:
606 return self.direction
608 def __str__(self):
609 x = []
610 if self.reflection:
611 if self.at_surface():
612 x.append('surface')
613 else:
614 if not self.headwave:
615 if self.direction == UP:
616 x.append('underside')
617 else:
618 x.append('upperside')
620 if self.headwave:
621 x.append('headwave propagation along')
622 elif self.reflection and self.conversion:
623 x.append('reflection with conversion from %s to %s' % (
624 smode(self.in_mode).upper(), smode(self.out_mode).upper()))
625 if not self.at_surface():
626 x.append('at')
627 elif self.reflection:
628 x.append('reflection')
629 if not self.at_surface():
630 x.append('at')
631 elif self.conversion:
632 x.append('conversion from %s to %s at' % (
633 smode(self.in_mode).upper(), smode(self.out_mode).upper()))
634 else:
635 x.append('passing through')
637 if isinstance(self.depth, float):
638 x.append('interface in %g km depth' % (self.depth/1000.))
639 else:
640 if not self.at_surface():
641 x.append('%s' % self.depth)
643 if not self.reflection:
644 if self.direction == UP:
645 x.append('on upgoing path')
646 else:
647 x.append('on downgoing path')
649 return ' '.join(x)
652class Head(Knee):
653 def __init__(self, *args):
654 if args:
655 z, in_direction, mode = args
656 Knee.__init__(self, z, in_direction, True, mode, mode)
657 else:
658 Knee.__init__(self)
660 def __str__(self):
661 x = ['propagation as headwave']
662 if isinstance(self.depth, float):
663 x.append('at interface in %g km depth' % (self.depth/1000.))
664 else:
665 x.append('at %s' % self.depth)
667 return ' '.join(x)
670class UnknownClassicPhase(CakeError):
671 def __init__(self, phasename):
672 self.phasename = phasename
674 def __str__(self):
675 return 'Unknown classic phase name: %s' % self.phasename
678class PhaseDefParseError(CakeError):
679 '''
680 Exception raised when an error occures during parsing of a phase
681 definition string.
682 '''
684 def __init__(self, definition, position, exception):
685 self.definition = definition
686 self.position = position
687 self.exception = exception
689 def __str__(self):
690 return 'Invalid phase definition: "%s" (at character %i: %s)' % (
691 self.definition, self.position+1, str(self.exception))
694class PhaseDef(object):
696 '''
697 Definition of a seismic phase arrival, based on ray propagation path.
699 :param definition: string representation of the phase in Cake's phase
700 syntax
702 Seismic phases are conventionally named e.g. P, Pn, PP, PcP, etc. In Cake,
703 a slightly different terminology is adapted, which allows to specify
704 arbitrary conversion/reflection histories for seismic ray paths. The
705 conventions used here are inspired by those used in the TauP toolkit, but
706 are not completely compatible with those.
708 The definition of a seismic ray propagation path in Cake's phase syntax is
709 a string consisting of an alternating sequence of *legs* and *knees*.
711 A *leg* represents seismic wave propagation without any conversions,
712 encountering only super-critical reflections. Legs are denoted by ``P``,
713 ``p``, ``S``, or ``s``. The capital letters are used when the take-off of
714 the *leg* is in downward direction, while the lower case letters indicate a
715 take-off in upward direction.
717 A *knee* is an interaction with an interface. It can be a mode conversion,
718 a reflection, or propagation as a headwave or diffracted wave.
720 * conversion is simply denoted as: ``(INTERFACE)`` or ``DEPTH``
721 * upperside reflection: ``v(INTERFACE)`` or ``vDEPTH``
722 * underside reflection: ``^(INTERFACE)`` or ``^DEPTH``
723 * normal kind headwave or diffracted wave: ``v_(INTERFACE)`` or
724 ``v_DEPTH``
726 The interface may be given by name or by depth: INTERFACE is the name of an
727 interface defined in the model, DEPTH is the depth of an interface in
728 [km] (the interface closest to that depth is chosen). If two legs appear
729 consecutively without an explicit *knee*, surface interaction is assumed.
731 The phase definition may end with a backslash ``\\``, to indicate that the
732 ray should arrive at the receiver from above instead of from below. It is
733 possible to restrict the maximum and minimum depth of a *leg* by appending
734 ``<(INTERFACE)`` or ``<DEPTH`` or ``>(INTERFACE)`` or ``>DEPTH`` after the
735 leg character, respectively.
737 **Examples:**
739 * ``P`` - like the classical P, but includes PKP, PKIKP, Pg
740 * ``P<(moho)`` - like classical Pg, but must leave source downwards
741 * ``pP`` - leaves source upward, reflects at surface, then travels as P
742 * ``P(moho)s`` - conversion from P to S at the Moho on upgoing path
743 * ``P(moho)S`` - conversion from P to S at the Moho on downgoing path
744 * ``Pv12p`` - P with reflection at 12 km deep interface (or the
745 interface closest to that)
746 * ``Pv_(moho)p`` - classical Pn
747 * ``Pv_(cmb)p`` - classical Pdiff
748 * ``P^(conrad)P`` - underside reflection of P at the Conrad
749 discontinuity
751 **Usage:**
753 >>> from pyrocko.cake import PhaseDef
754 # must escape the backslash
755 >>> my_crazy_phase = PhaseDef('pPv(moho)sP\\\\')
756 >>> print my_crazy_phase
757 Phase definition "pPv(moho)sP\":
758 - P mode propagation, departing upward
759 - surface reflection
760 - P mode propagation, departing downward
761 - upperside reflection with conversion from P to S at moho
762 - S mode propagation, departing upward
763 - surface reflection with conversion from S to P
764 - P mode propagation, departing downward
765 - arriving at target from above
767 .. note::
769 (1) These conventions might be extended in a way to allow to fix wave
770 propagation to SH mode, possibly by specifying SH, or a single
771 character (e.g. H) instead of S. This would be benificial for the
772 selection of conversion and reflection coefficients, which
773 currently only deal with the P-SV case.
774 '''
776 allowed_characters_pattern = r'[0-9a-zA-Z_()<>^v\\.]+'
777 allowed_characters_pattern_classic = r'[a-zA-Z0-9]+'
779 @staticmethod
780 def classic_definitions():
781 defs = {}
782 # PmP, PmS, PcP, PcS, SmP, ...
783 for r in 'mc':
784 for a, b in 'PP PS SS SP'.split():
785 defs[a+r+b] = [
786 '%sv(%s)%s' % (a, {'m': 'moho', 'c': 'cmb'}[r], b.lower())]
788 # Pg, P, S, Sg
789 for a in 'PS':
790 defs[a+'g'] = ['%s<(moho)' % x for x in (a, a.lower())]
791 defs[a] = ['%s<(cmb)(moho)%s' % (x, x.lower()) for x in (
792 a, a.lower())]
794 defs[a.lower()] = [a.lower()]
796 for a, b in 'PP PS SS SP'.split():
797 defs[a+'K'+b] = ['%s(cmb)P<(icb)(cmb)%s' % (a, b.lower())]
798 defs[a+'KIK'+b] = ['%s(cmb)P(icb)P(icb)p(cmb)%s' % (a, b.lower())]
799 defs[a+'KJK'+b] = ['%s(cmb)P(icb)S(icb)p(cmb)%s' % (a, b.lower())]
800 defs[a+'KiK'+b] = ['%s(cmb)Pv(icb)p(cmb)%s' % (a, b.lower())]
802 # PP, SS, PS, SP, PPP, ...
803 for a in 'PS':
804 for b in 'PS':
805 for c in 'PS':
806 defs[a+b+c] = [''.join(defs[x][0] for x in a+b+c)]
808 defs[a+b] = [''.join(defs[x][0] for x in a+b)]
810 # Pc, Pdiff, Sc, ...
811 for x in 'PS':
812 defs[x+'c'] = defs[x+'diff'] = [x+'v_(cmb)'+x.lower()]
813 defs[x+'n'] = [x+'v_(moho)'+x.lower()]
815 # depth phases
816 for k in list(defs.keys()):
817 if k not in 'ps':
818 for x in 'ps':
819 defs[x+k] = [x + defs[k][0]]
821 return defs
823 @staticmethod
824 def classic(phasename):
825 '''
826 Get phase definitions based on classic phase name.
828 :param phasename: classic name of a phase
829 :returns: list of PhaseDef objects
831 This returns a list of PhaseDef objects, because some classic phases
832 (like e.g. Pg) can only be represented by two Cake style PhaseDef
833 objects (one with downgoing and one with upgoing first leg).
834 '''
836 defs = PhaseDef.classic_definitions()
837 if phasename not in defs:
838 raise UnknownClassicPhase(phasename)
840 return [PhaseDef(d, classicname=phasename) for d in defs[phasename]]
842 def __init__(self, definition=None, classicname=None):
844 state = 0
845 sdepth = ''
846 sinterface = ''
847 depthmax = depthmin = None
848 depthlim = None
849 depthlimtype = None
850 sdepthlim = ''
851 events = []
852 direction_stop = UP
853 need_leg = True
854 ic = 0
855 if definition is not None:
856 knee = Knee()
857 try:
858 for ic, c in enumerate(definition):
860 if state in (0, 1):
862 if c in '0123456789.':
863 need_leg = True
864 state = 1
865 sdepth += c
866 continue
868 elif state == 1:
869 knee.depth = float(sdepth)*1000.
870 state = 0
872 if state == 2:
873 if c == ')':
874 knee.depth = sinterface
875 state = 0
876 else:
877 sinterface += c
879 continue
881 if state in (3, 4):
883 if state == 3:
884 if c in '0123456789.':
885 sdepthlim += c
886 continue
887 elif c == '(':
888 state = 4
889 continue
890 else:
891 depthlim = float(sdepthlim)*1000.
892 if depthlimtype == '<':
893 depthmax = depthlim
894 else:
895 depthmin = depthlim
896 state = 0
898 elif state == 4:
899 if c == ')':
900 depthlim = sdepthlim
901 if depthlimtype == '<':
902 depthmax = depthlim
903 else:
904 depthmin = depthlim
905 state = 0
906 continue
907 else:
908 sdepthlim += c
909 continue
911 if state == 0:
913 if c == '(':
914 need_leg = True
915 state = 2
916 continue
918 elif c in '<>':
919 state = 3
920 depthlim = None
921 sdepthlim = ''
922 depthlimtype = c
923 continue
925 elif c in 'psPS':
926 leg = Leg()
927 if c in 'ps':
928 leg.departure = UP
929 else:
930 leg.departure = DOWN
931 leg.mode = imode(c)
933 if events:
934 in_leg = events[-1]
935 if depthmin is not None:
936 in_leg.set_depthmin(depthmin)
937 depthmin = None
938 if depthmax is not None:
939 in_leg.set_depthmax(depthmax)
940 depthmax = None
942 if in_leg.mode != leg.mode:
943 knee.conversion = True
944 else:
945 knee.conversion = False
947 if not knee.reflection:
948 if c in 'ps':
949 knee.direction = UP
950 else:
951 knee.direction = DOWN
953 knee.set_modes(in_leg, leg)
954 knee.in_setup_state = False
955 events.append(knee)
956 knee = Knee()
957 sdepth = ''
958 sinterface = ''
960 events.append(leg)
961 need_leg = False
962 continue
964 elif c == '^':
965 need_leg = True
966 knee.direction = UP
967 knee.reflection = True
968 continue
970 elif c == 'v':
971 need_leg = True
972 knee.direction = DOWN
973 knee.reflection = True
974 continue
976 elif c == '_':
977 need_leg = True
978 knee.headwave = True
979 continue
981 elif c == '\\':
982 direction_stop = DOWN
983 continue
985 else:
986 raise PhaseDefParseError(
987 definition, ic, 'invalid character: "%s"' % c)
989 if state == 3:
990 depthlim = float(sdepthlim)*1000.
991 if depthlimtype == '<':
992 depthmax = depthlim
993 else:
994 depthmin = depthlim
995 state = 0
997 except (ValueError, InvalidKneeDef) as e:
998 raise PhaseDefParseError(definition, ic, e)
1000 if state != 0 or need_leg:
1001 raise PhaseDefParseError(
1002 definition, ic, 'unfinished expression')
1004 if events and depthmin is not None:
1005 events[-1].set_depthmin(depthmin)
1006 if events and depthmax is not None:
1007 events[-1].set_depthmax(depthmax)
1009 self._definition = definition
1010 self._classicname = classicname
1011 self._events = events
1012 self._direction_stop = direction_stop
1014 def __iter__(self):
1015 for ev in self._events:
1016 yield ev
1018 def append(self, ev):
1019 self._events.append(ev)
1021 def first_leg(self):
1022 '''
1023 Get the first leg in phase definition.
1024 '''
1025 return self._events[0]
1027 def last_leg(self):
1028 '''
1029 Get the last leg in phase definition.
1030 '''
1031 return self._events[-1]
1033 def legs(self):
1034 '''
1035 Iterate over the continuous pieces of wave propagation (legs) defined
1036 within this phase definition.
1037 '''
1039 return (leg for leg in self if isinstance(leg, Leg))
1041 def knees(self):
1042 '''
1043 Iterate over conversions and reflections (knees) defined within this
1044 phase definition.
1045 '''
1046 return (knee for knee in self if isinstance(knee, Knee))
1048 def definition(self):
1049 '''
1050 Get original definition of the phase.
1051 '''
1052 return self._definition
1054 def given_name(self):
1055 '''
1056 Get entered classic name if any, or original definition of the phase.
1057 '''
1059 if self._classicname:
1060 return self._classicname
1061 else:
1062 return self._definition
1064 def direction_start(self):
1065 return self.first_leg().departure
1067 def direction_stop(self):
1068 return self._direction_stop
1070 def headwave_knee(self):
1071 for el in self:
1072 if type(el) == Knee and el.headwave:
1073 return el
1074 return None
1076 def used_repr(self):
1077 '''
1078 Translate into textual representation (cake phase syntax).
1079 '''
1080 def strdepth(x):
1081 if isinstance(x, float):
1082 return '%g' % (x/1000.)
1083 else:
1084 return '(%s)' % x
1086 x = []
1087 for el in self:
1088 if type(el) == Leg:
1089 if el.departure == UP:
1090 x.append(smode(el.mode).lower())
1091 else:
1092 x.append(smode(el.mode).upper())
1094 if el.depthmax is not None:
1095 x.append('<'+strdepth(el.depthmax))
1097 if el.depthmin is not None:
1098 x.append('>'+strdepth(el.depthmin))
1100 elif type(el) == Knee:
1101 if el.reflection and not el.at_surface():
1102 if el.direction == DOWN:
1103 x.append('v')
1104 else:
1105 x.append('^')
1106 if el.headwave:
1107 x.append('_')
1108 if not el.at_surface():
1109 x.append(strdepth(el.depth))
1111 elif type(el) == Head:
1112 x.append('_')
1113 x.append(strdepth(el.depth))
1115 if self._direction_stop == DOWN:
1116 x.append('\\')
1118 return ''.join(x)
1120 def __repr__(self):
1121 if self._definition is not None:
1122 return "PhaseDef('%s')" % self._definition
1123 else:
1124 return "PhaseDef('%s')" % self.used_repr()
1126 def __str__(self):
1127 orig = ''
1128 used = self.used_repr()
1129 if self._definition != used:
1130 orig = ' (entered as "%s")' % self._definition
1132 sarrive = '\n - arriving at target from %s' % ('below', 'above')[
1133 self._direction_stop == DOWN]
1135 return 'Phase definition "%s"%s:\n - ' % (used, orig) + \
1136 '\n - '.join(str(ev) for ev in self) + sarrive
1138 def copy(self):
1139 '''
1140 Get a deep copy of it.
1141 '''
1142 return copy.deepcopy(self)
1145def to_phase_defs(phases):
1146 if isinstance(phases, (str, PhaseDef)):
1147 phases = [phases]
1149 phases_out = []
1150 for phase in phases:
1151 if isinstance(phase, str):
1152 phases_out.extend(PhaseDef(x.strip()) for x in phase.split(','))
1153 elif isinstance(phase, PhaseDef):
1154 phases_out.append(phase)
1155 else:
1156 raise PhaseDefParseError('invalid phase definition')
1158 return phases_out
1161def csswap(x):
1162 return cmath.sqrt(1.-x**2)
1165def psv_surface_ind(in_mode, out_mode):
1166 '''
1167 Get indices to select the appropriate element from scatter matrix for free
1168 surface.
1169 '''
1171 return (int(in_mode == S), int(out_mode == S))
1174def psv_surface(material, p, energy=False):
1175 '''
1176 Scatter matrix for free surface reflection/conversions.
1178 :param material: material, object of type :py:class:`Material`
1179 :param p: flat ray parameter [s/m]
1180 :param energy: bool, when ``True`` energy normalized coefficients are
1181 returned
1182 :returns: Scatter matrix
1184 The scatter matrix is ordered as follows::
1186 [[PP, PS],
1187 [SP, SS]]
1189 The formulas given in Aki & Richards are used.
1190 '''
1192 vp, vs, rho = material.vp, material.vs, material.rho
1193 sinphi = p * vp
1194 sinlam = p * vs
1195 cosphi = csswap(sinphi)
1196 coslam = csswap(sinlam)
1198 if vs == 0.0:
1199 scatter = num.array([[-1.0, 0.0], [0.0, 1.0]])
1201 else:
1202 vsp_term = (1.0/vs**2 - 2.0*p**2)
1203 pcc_term = 4.0 * p**2 * cosphi/vp * coslam/vs
1204 denom = vsp_term**2 + pcc_term
1206 scatter = num.array([
1207 [- vsp_term**2 + pcc_term, 4.0*p*coslam/vp*vsp_term],
1208 [4.0*p*cosphi/vs*vsp_term, vsp_term**2 - pcc_term]],
1209 dtype=complex) / denom
1211 if not energy:
1212 return scatter
1213 else:
1214 eps = 1e-16
1215 normvec = num.array([vp*rho*cosphi+eps, vs*rho*coslam+eps])
1216 escatter = scatter*num.conj(scatter) * num.real(
1217 (normvec[:, num.newaxis]) / (normvec[num.newaxis, :]))
1218 return num.real(escatter)
1221def psv_solid_ind(in_direction, out_direction, in_mode, out_mode):
1222 '''
1223 Get indices to select the appropriate element from scatter matrix for
1224 solid-solid interface.
1225 '''
1227 return (
1228 (out_direction == DOWN)*2 + (out_mode == S),
1229 (in_direction == UP)*2 + (in_mode == S))
1232def psv_solid(material1, material2, p, energy=False):
1233 '''
1234 Scatter matrix for solid-solid interface.
1236 :param material1: material above, object of type :py:class:`Material`
1237 :param material2: material below, object of type :py:class:`Material`
1238 :param p: flat ray parameter [s/m]
1239 :param energy: bool, when ``True`` energy normalized coefficients are
1240 returned
1241 :returns: Scatter matrix
1243 The scatter matrix is ordered as follows::
1245 [[P1P1, S1P1, P2P1, S2P1],
1246 [P1S1, S1S1, P2S1, S2S1],
1247 [P1P2, S1P2, P2P2, S2P2],
1248 [P1S2, S1S2, P2S2, S2S2]]
1250 The formulas given in Aki & Richards are used.
1251 '''
1253 vp1, vs1, rho1 = material1.vp, material1.vs, material1.rho
1254 vp2, vs2, rho2 = material2.vp, material2.vs, material2.rho
1256 sinphi1 = p * vp1
1257 cosphi1 = csswap(sinphi1)
1258 sinlam1 = p * vs1
1259 coslam1 = csswap(sinlam1)
1260 sinphi2 = p * vp2
1261 cosphi2 = csswap(sinphi2)
1262 sinlam2 = p * vs2
1263 coslam2 = csswap(sinlam2)
1265 # from aki and richards
1266 M = num.array([
1267 [-vp1*p, -coslam1, vp2*p, coslam2],
1268 [cosphi1, -vs1*p, cosphi2, -vs2*p],
1269 [2.0*rho1*vs1**2*p*cosphi1, rho1*vs1*(1.0-2.0*vs1**2*p**2),
1270 2.0*rho2*vs2**2*p*cosphi2, rho2*vs2*(1.0-2.0*vs2**2*p**2)],
1271 [-rho1*vp1*(1.0-2.0*vs1**2*p**2), 2.0*rho1*vs1**2*p*coslam1,
1272 rho2*vp2*(1.0-2.0*vs2**2*p**2), -2.0*rho2*vs2**2*p*coslam2]],
1273 dtype=complex)
1274 N = M.copy()
1275 N[0] *= -1.0
1276 N[3] *= -1.0
1278 scatter = num.dot(num.linalg.inv(M), N)
1280 if not energy:
1281 return scatter
1282 else:
1283 eps = 1e-16
1284 if vs1 == 0.:
1285 vs1 = vp1*1e-16
1286 if vs2 == 0.:
1287 vs2 = vp2*1e-16
1288 normvec = num.array([
1289 vp1*rho1*(cosphi1+eps), vs1*rho1*(coslam1+eps),
1290 vp2*rho2*(cosphi2+eps), vs2*rho2*(coslam2+eps)], dtype=complex)
1291 escatter = scatter*num.conj(scatter) * num.real(
1292 normvec[:, num.newaxis] / normvec[num.newaxis, :])
1294 return num.real(escatter)
1297class BadPotIntCoefs(CakeError):
1298 pass
1301def potint_coefs(c1, c2, r1, r2): # r2 > r1
1302 eps = r2*1e-9
1303 if c1 == 0. and c2 == 0.:
1304 c1c2 = 1.
1305 else:
1306 c1c2 = c1/c2
1307 b = math.log(c1c2)/math.log((r1+eps)/r2)
1308 if abs(b) > 10.:
1309 raise BadPotIntCoefs()
1310 a = c1/(r1+eps)**b
1311 return a, b
1314def imode(s):
1315 if s.lower() == 'p':
1316 return P
1317 elif s.lower() == 's':
1318 return S
1321def smode(i):
1322 if i == P:
1323 return 'p'
1324 elif i == S:
1325 return 's'
1328class PathFailed(CakeError):
1329 pass
1332class SurfaceReached(PathFailed):
1333 pass
1336class BottomReached(PathFailed):
1337 pass
1340class MaxDepthReached(PathFailed):
1341 pass
1344class MinDepthReached(PathFailed):
1345 pass
1348class Trapped(PathFailed):
1349 pass
1352class NotPhaseConform(PathFailed):
1353 pass
1356class CannotPropagate(PathFailed):
1357 def __init__(self, direction, ilayer):
1358 PathFailed.__init__(self)
1359 self._direction = direction
1360 self._ilayer = ilayer
1362 def __str__(self):
1363 return 'Cannot enter layer %i from %s' % (
1364 self._ilayer, {
1365 UP: 'below',
1366 DOWN: 'above'}[self._direction])
1369class Layer(object):
1370 '''
1371 Representation of a layer in a layered earth model.
1373 :param ztop: depth of top of layer
1374 :param zbot: depth of bottom of layer
1375 :param name: name of layer (optional)
1377 Subclasses are: :py:class:`HomogeneousLayer` and :py:class:`GradientLayer`.
1378 '''
1380 def __init__(self, ztop, zbot, name=None):
1381 self.ztop = ztop
1382 self.zbot = zbot
1383 self.zmid = (self.ztop + self.zbot) * 0.5
1384 self.name = name
1385 self.ilayer = None
1387 def _update_potint_coefs(self):
1388 potint_p = potint_s = False
1389 try:
1390 self._ppic = potint_coefs(
1391 self.mbot.vp, self.mtop.vp,
1392 radius(self.zbot), radius(self.ztop))
1393 potint_p = True
1394 except BadPotIntCoefs:
1395 pass
1397 potint_s = False
1398 try:
1399 self._spic = potint_coefs(
1400 self.mbot.vs, self.mtop.vs,
1401 radius(self.zbot), radius(self.ztop))
1402 potint_s = True
1403 except BadPotIntCoefs:
1404 pass
1406 assert P == 1 and S == 2
1407 self._use_potential_interpolation = (None, potint_p, potint_s)
1409 def potint_coefs(self, mode):
1410 '''
1411 Get coefficients for potential interpolation.
1413 :param mode: mode of wave propagation, :py:const:`P` or :py:const:`S`
1414 :returns: coefficients ``(a, b)``
1415 '''
1417 if mode == P:
1418 return self._ppic
1419 else:
1420 return self._spic
1422 def contains(self, z):
1423 '''
1424 Tolerantly check if a given depth is within the layer
1425 (including boundaries).
1426 '''
1428 return self.ztop <= z <= self.zbot or \
1429 self.at_bottom(z) or self.at_top(z)
1431 def inner(self, z):
1432 '''
1433 Tolerantly check if a given depth is within the layer
1434 (not including boundaries).
1435 '''
1437 return self.ztop <= z <= self.zbot and not \
1438 self.at_bottom(z) and not \
1439 self.at_top(z)
1441 def at_bottom(self, z):
1442 '''
1443 Tolerantly check if given depth is at the bottom of the layer.
1444 '''
1446 return abs(self.zbot - z) < ZEPS
1448 def at_top(self, z):
1449 '''
1450 Tolerantly check if given depth is at the top of the layer.
1451 '''
1452 return abs(self.ztop - z) < ZEPS
1454 def pflat_top(self, p):
1455 '''
1456 Convert spherical ray parameter to local flat ray parameter for top of
1457 layer.
1458 '''
1459 return p / (earthradius-self.ztop)
1461 def pflat_bottom(self, p):
1462 '''
1463 Convert spherical ray parameter to local flat ray parameter for bottom
1464 of layer.
1465 '''
1466 return p / (earthradius-self.zbot)
1468 def pflat(self, p, z):
1469 '''
1470 Convert spherical ray parameter to local flat ray parameter for
1471 given depth.
1472 '''
1473 return p / (earthradius-z)
1475 def v_potint(self, mode, z):
1476 a, b = self.potint_coefs(mode)
1477 return a*(earthradius-z)**b
1479 def u_potint(self, mode, z):
1480 a, b = self.potint_coefs(mode)
1481 return 1./(a*(earthradius-z)**b)
1483 def xt_potint(self, p, mode, zpart=None):
1484 '''
1485 Get travel time and distance for for traversal with given mode and ray
1486 parameter.
1488 :param p: ray parameter (spherical)
1489 :param mode: mode of propagation (:py:const:`P` or :py:const:`S`)
1490 :param zpart: if given, tuple with two depths to restrict computation
1491 to a part of the layer
1493 This implementation uses analytic formulas valid for a spherical earth
1494 in the case where the velocity c within the layer is given by potential
1495 interpolation of the form
1497 c(z) = a*z^b
1498 '''
1499 utop, ubot = self.u_top_bottom(mode)
1500 a, b = self.potint_coefs(mode)
1501 ztop = self.ztop
1502 zbot = self.zbot
1503 if zpart is not None:
1504 utop = self.u(mode, zpart[0])
1505 ubot = self.u(mode, zpart[1])
1506 ztop, zbot = zpart
1507 utop = 1./(a*(earthradius-ztop)**b)
1508 ubot = 1./(a*(earthradius-zbot)**b)
1510 r1 = radius(zbot)
1511 r2 = radius(ztop)
1512 burger_eta1 = r1 * ubot
1513 burger_eta2 = r2 * utop
1514 if b != 1:
1515 def cpe(eta):
1516 return num.arccos(num.minimum(p/num.maximum(eta, p/2), 1.0))
1518 def sep(eta):
1519 return num.sqrt(num.maximum(eta**2 - p**2, 0.0))
1521 x = (cpe(burger_eta2)-cpe(burger_eta1))/(1.0-b)
1522 t = (sep(burger_eta2)-sep(burger_eta1))/(1.0-b)
1523 else:
1524 lr = math.log(r2/r1)
1525 sap = num.sqrt(1.0/a**2 - p**2)
1526 x = p/sap * lr
1527 t = 1./(a**2 * sap)
1529 x *= r2d
1531 return x, t
1533 def test(self, p, mode, z):
1534 '''
1535 Check if wave mode can exist for given ray parameter at given depth
1536 within the layer.
1537 '''
1538 return (self.u(mode, z)*radius(z) - p) > 0.
1540 def tests(self, p, mode):
1541 utop, ubot = self.u_top_bottom(mode)
1542 return (
1543 (utop * radius(self.ztop) - p) > 0.,
1544 (ubot * radius(self.zbot) - p) > 0.)
1546 def zturn_potint(self, p, mode):
1547 '''
1548 Get turning depth for given ray parameter and propagation mode.
1549 '''
1551 a, b = self.potint_coefs(mode)
1552 r = num.exp(num.log(a*p)/(1.0-b))
1553 return earthradius-r
1555 def propagate(self, p, mode, direction):
1556 '''
1557 Propagate ray through layer.
1559 :param p: ray parameter
1560 :param mode: propagation mode
1561 :param direction: in direction (:py:const:`UP` or :py:const:`DOWN`
1562 '''
1563 if direction == DOWN:
1564 zin, zout = self.ztop, self.zbot
1565 else:
1566 zin, zout = self.zbot, self.ztop
1568 if self.v(mode, zin) == 0.0 or not self.test(p, mode, zin):
1569 raise CannotPropagate(direction, self.ilayer)
1571 if not self.test(p, mode, zout):
1572 return -direction
1573 else:
1574 return direction
1576 def resize(self, depth_min=None, depth_max=None):
1577 '''
1578 Change layer thinkness and interpolate material if required.
1579 '''
1580 if depth_min:
1581 mtop = self.material(depth_min)
1583 if depth_max:
1584 mbot = self.material(depth_max)
1586 self.mtop = mtop if depth_min else self.mtop
1587 self.mbot = mbot if depth_max else self.mbot
1588 self.ztop = depth_min if depth_min else self.ztop
1589 self.zbot = depth_max if depth_max else self.zbot
1590 self.zmid = self.ztop + (self.zbot - self.ztop)/2.
1593class DoesNotTurn(CakeError):
1594 pass
1597def radius(z):
1598 return earthradius - z
1601class HomogeneousLayer(Layer):
1602 '''
1603 Representation of a homogeneous layer in a layered earth model.
1605 Base class: :py:class:`Layer`.
1606 '''
1608 def __init__(self, ztop, zbot, m, name=None):
1609 Layer.__init__(self, ztop, zbot, name=name)
1610 self.m = m
1611 self.mtop = m
1612 self.mbot = m
1613 self._update_potint_coefs()
1615 def copy(self, ztop=None, zbot=None):
1616 if ztop is None:
1617 ztop = self.ztop
1619 if zbot is None:
1620 zbot = self.zbot
1622 return HomogeneousLayer(ztop, zbot, self.m, name=self.name)
1624 def material(self, z):
1625 return self.m
1627 def u(self, mode, z=None):
1628 if self._use_potential_interpolation[mode] and z is not None:
1629 return self.u_potint(mode, z)
1631 if mode == P:
1632 return 1./self.m.vp
1633 if mode == S:
1634 return 1./self.m.vs
1636 def u_top_bottom(self, mode):
1637 u = self.u(mode)
1638 return u, u
1640 def v(self, mode, z=None):
1641 if self._use_potential_interpolation[mode] and z is not None:
1642 return self.v_potint(mode, z)
1644 if mode == P:
1645 v = self.m.vp
1646 if mode == S:
1647 v = self.m.vs
1649 if num.isscalar(z):
1650 return v
1651 else:
1652 return filled(v, len(z))
1654 def v_top_bottom(self, mode):
1655 v = self.v(mode)
1656 return v, v
1658 def xt(self, p, mode, zpart=None):
1659 if self._use_potential_interpolation[mode]:
1660 return self.xt_potint(p, mode, zpart)
1662 u = self.u(mode)
1663 pflat = self.pflat_bottom(p)
1664 if zpart is None:
1665 dz = (self.zbot - self.ztop)
1666 else:
1667 dz = abs(zpart[1]-zpart[0])
1669 u = self.u(mode)
1670 eps = u*0.001
1671 denom = num.sqrt(u**2 - pflat**2) + eps
1673 x = r2d*pflat/(earthradius-self.zmid) * dz / denom
1674 t = u**2 * dz / denom
1675 return x, t
1677 def zturn(self, p, mode):
1678 if self._use_potential_interpolation[mode]:
1679 return self.zturn_potint(p, mode)
1681 raise DoesNotTurn()
1683 def split(self, z):
1684 upper = HomogeneousLayer(self.ztop, z, self.m, name=self.name)
1685 lower = HomogeneousLayer(z, self.zbot, self.m, name=self.name)
1686 upper.ilayer = self.ilayer
1687 lower.ilayer = self.ilayer
1688 return upper, lower
1690 def __str__(self):
1691 if self.name:
1692 name = self.name + ' '
1693 else:
1694 name = ''
1696 calcmode = ''.join('HP'[self._use_potential_interpolation[mode]]
1697 for mode in (P, S))
1699 return ' (%i) homogeneous layer %s(%g km - %g km) [%s]\n %s' % (
1700 self.ilayer, name, self.ztop/km, self.zbot/km, calcmode, self.m)
1703class GradientLayer(Layer):
1704 '''
1705 Representation of a gradient layer in a layered earth model.
1707 Base class: :py:class:`Layer`.
1708 '''
1710 def __init__(self, ztop, zbot, mtop, mbot, name=None):
1711 Layer.__init__(self, ztop, zbot, name=name)
1712 self.mtop = mtop
1713 self.mbot = mbot
1714 self._update_potint_coefs()
1716 def copy(self, ztop=None, zbot=None):
1717 if ztop is None:
1718 ztop = self.ztop
1720 if zbot is None:
1721 zbot = self.zbot
1723 return GradientLayer(ztop, zbot, self.mtop, self.mbot, name=self.name)
1725 def interpolate(self, z, ptop, pbot):
1726 return ptop + (z - self.ztop)*(pbot - ptop)/(self.zbot-self.ztop)
1728 def material(self, z):
1729 dtop = self.mtop.astuple()
1730 dbot = self.mbot.astuple()
1731 d = [
1732 self.interpolate(z, ptop, pbot)
1733 for (ptop, pbot) in zip(dtop, dbot)]
1735 return Material(*d)
1737 def u_top_bottom(self, mode):
1738 if mode == P:
1739 return 1./self.mtop.vp, 1./self.mbot.vp
1740 if mode == S:
1741 return 1./self.mtop.vs, 1./self.mbot.vs
1743 def u(self, mode, z):
1744 if self._use_potential_interpolation[mode]:
1745 return self.u_potint(mode, z)
1747 if mode == P:
1748 return 1./self.interpolate(z, self.mtop.vp, self.mbot.vp)
1749 if mode == S:
1750 return 1./self.interpolate(z, self.mtop.vs, self.mbot.vs)
1752 def v_top_bottom(self, mode):
1753 if mode == P:
1754 return self.mtop.vp, self.mbot.vp
1755 if mode == S:
1756 return self.mtop.vs, self.mbot.vs
1758 def v(self, mode, z):
1759 if self._use_potential_interpolation[mode]:
1760 return self.v_potint(mode, z)
1762 if mode == P:
1763 return self.interpolate(z, self.mtop.vp, self.mbot.vp)
1764 if mode == S:
1765 return self.interpolate(z, self.mtop.vs, self.mbot.vs)
1767 def xt(self, p, mode, zpart=None):
1768 if self._use_potential_interpolation[mode]:
1769 return self.xt_potint(p, mode, zpart)
1771 utop, ubot = self.u_top_bottom(mode)
1772 b = (1./ubot - 1./utop)/(self.zbot - self.ztop)
1774 pflat = self.pflat_bottom(p)
1775 if zpart is not None:
1776 utop = self.u(mode, zpart[0])
1777 ubot = self.u(mode, zpart[1])
1779 peps = 1e-16
1780 pdp = pflat + peps
1782 def func(u):
1783 eta = num.sqrt(num.maximum(u**2 - pflat**2, 0.0))
1784 xx = eta/u
1785 tt = num.where(
1786 pflat <= u,
1787 num.log(u+eta) - num.log(pdp) - eta/u,
1788 0.0)
1790 return xx, tt
1792 xxtop, tttop = func(utop)
1793 xxbot, ttbot = func(ubot)
1795 x = (xxtop - xxbot) / (b*pdp)
1796 t = (tttop - ttbot) / b + pflat*x
1798 x *= r2d/(earthradius - self.zmid)
1799 return x, t
1801 def zturn(self, p, mode):
1802 if self._use_potential_interpolation[mode]:
1803 return self.zturn_potint(p, mode)
1804 pflat = self.pflat_bottom(p)
1805 vtop, vbot = self.v_top_bottom(mode)
1806 return (1./pflat - vtop) * (self.zbot - self.ztop) / \
1807 (vbot-vtop) + self.ztop
1809 def split(self, z):
1810 mmid = self.material(z)
1811 upper = GradientLayer(self.ztop, z, self.mtop, mmid, name=self.name)
1812 lower = GradientLayer(z, self.zbot, mmid, self.mbot, name=self.name)
1813 upper.ilayer = self.ilayer
1814 lower.ilayer = self.ilayer
1815 return upper, lower
1817 def __str__(self):
1818 if self.name:
1819 name = self.name + ' '
1820 else:
1821 name = ''
1823 calcmode = ''.join('HP'[self._use_potential_interpolation[mode]]
1824 for mode in (P, S))
1826 return ''' (%i) gradient layer %s(%g km - %g km) [%s]
1827 %s
1828 %s''' % (
1829 self.ilayer,
1830 name,
1831 self.ztop/km,
1832 self.zbot/km,
1833 calcmode,
1834 self.mtop,
1835 self.mbot)
1838class Discontinuity(object):
1839 '''
1840 Base class for discontinuities in layered earth model.
1842 Subclasses are: :py:class:`Interface` and :py:class:`Surface`.
1843 '''
1845 def __init__(self, z, name=None):
1846 self.z = z
1847 self.zbot = z
1848 self.ztop = z
1849 self.name = name
1851 def change_depth(self, z):
1852 self.z = z
1853 self.zbot = z
1854 self.ztop = z
1856 def copy(self):
1857 return copy.deepcopy(self)
1860class Interface(Discontinuity):
1861 '''
1862 Representation of an interface in a layered earth model.
1864 Base class: :py:class:`Discontinuity`.
1865 '''
1867 def __init__(self, z, mabove, mbelow, name=None):
1868 Discontinuity.__init__(self, z, name)
1869 self.mabove = mabove
1870 self.mbelow = mbelow
1872 def __str__(self):
1873 if self.name is None:
1874 return 'interface'
1875 else:
1876 return 'interface "%s"' % self.name
1878 def u_top_bottom(self, mode):
1879 if mode == P:
1880 return reci_or_none(self.mabove.vp), reci_or_none(self.mbelow.vp)
1881 if mode == S:
1882 return reci_or_none(self.mabove.vs), reci_or_none(self.mbelow.vs)
1884 def critical_ps(self, mode):
1885 uabove, ubelow = self.u_top_bottom(mode)
1886 return (
1887 mult_or_none(uabove, radius(self.z)),
1888 mult_or_none(ubelow, radius(self.z)))
1890 def propagate(self, p, mode, direction):
1891 uabove, ubelow = self.u_top_bottom(mode)
1892 if direction == DOWN:
1893 if ubelow is not None and ubelow*radius(self.z) - p >= 0:
1894 return direction
1895 else:
1896 return -direction
1897 if direction == UP:
1898 if uabove is not None and uabove*radius(self.z) - p >= 0:
1899 return direction
1900 else:
1901 return -direction
1903 def pflat(self, p):
1904 return p / (earthradius-self.z)
1906 def efficiency(self, in_direction, out_direction, in_mode, out_mode, p):
1907 scatter = psv_solid(
1908 self.mabove, self.mbelow, self.pflat(p), energy=True)
1909 return scatter[
1910 psv_solid_ind(in_direction, out_direction, in_mode, out_mode)]
1913class Surface(Discontinuity):
1914 '''
1915 Representation of the surface discontinuity in a layered earth model.
1917 Base class: :py:class:`Discontinuity`.
1918 '''
1920 def __init__(self, z, mbelow):
1921 Discontinuity.__init__(self, z, 'surface')
1922 self.z = z
1923 self.mbelow = mbelow
1925 def propagate(self, p, mode, direction):
1926 return direction # no implicit reflection at surface
1928 def u_top_bottom(self, mode):
1929 if mode == P:
1930 return None, reci_or_none(self.mbelow.vp)
1931 if mode == S:
1932 return None, reci_or_none(self.mbelow.vs)
1934 def critical_ps(self, mode):
1935 _, ubelow = self.u_top_bottom(mode)
1936 return None, mult_or_none(ubelow, radius(self.z))
1938 def pflat(self, p):
1939 return p / (earthradius-self.z)
1941 def efficiency(self, in_direction, out_direction, in_mode, out_mode, p):
1942 if in_direction == DOWN or out_direction == UP:
1943 return 0.0
1944 else:
1945 return psv_surface(
1946 self.mbelow, self.pflat(p), energy=True)[
1947 psv_surface_ind(in_mode, out_mode)]
1949 def __str__(self):
1950 return 'surface'
1953class Walker(object):
1954 def __init__(self, elements):
1955 self._elements = elements
1956 self._i = 0
1958 def current(self):
1959 return self._elements[self._i]
1961 def go(self, direction):
1962 if direction == UP:
1963 self.up()
1964 else:
1965 self.down()
1967 def down(self):
1968 if self._i < len(self._elements)-1:
1969 self._i += 1
1970 else:
1971 raise BottomReached()
1973 def up(self):
1974 if self._i > 0:
1975 self._i -= 1
1976 else:
1977 raise SurfaceReached()
1979 def goto_layer(self, layer):
1980 self._i = self._elements.index(layer)
1983class RayElement(object):
1984 '''
1985 An element of a :py:class:`RayPath`.
1986 '''
1988 def __eq__(self, other):
1989 return type(self) == type(other) and self.__dict__ == other.__dict__
1991 def is_straight(self):
1992 return isinstance(self, Straight)
1994 def is_kink(self):
1995 return isinstance(self, Kink)
1998class Straight(RayElement):
1999 '''
2000 A ray segment representing wave propagation through one :py:class:`Layer`.
2001 '''
2003 def __init__(self, direction_in, direction_out, mode, layer):
2004 self.mode = mode
2005 self._direction_in = direction_in
2006 self._direction_out = direction_out
2007 self.layer = layer
2009 def angle_in(self, p, endgaps=None):
2010 z = self.z_in(endgaps)
2011 dir = self.eff_direction_in(endgaps)
2012 v = self.layer.v(self.mode, z)
2013 pf = self.layer.pflat(p, z)
2015 if dir == DOWN:
2016 return num.arcsin(v*pf)*r2d
2017 else:
2018 return 180.-num.arcsin(v*pf)*r2d
2020 def angle_out(self, p, endgaps=None):
2021 z = self.z_out(endgaps)
2022 dir = self.eff_direction_out(endgaps)
2023 v = self.layer.v(self.mode, z)
2024 pf = self.layer.pflat(p, z)
2026 if dir == DOWN:
2027 return 180.-num.arcsin(v*pf)*r2d
2028 else:
2029 return num.arcsin(v*pf)*r2d
2031 def pflat_in(self, p, endgaps=None):
2032 return p / (earthradius-self.z_in(endgaps))
2034 def pflat_out(self, p, endgaps=None):
2035 return p / (earthradius-self.z_out(endgaps))
2037 def test(self, p, z):
2038 return self.layer.test(p, self.mode, z)
2040 def z_in(self, endgaps=None):
2041 if endgaps is not None:
2042 return endgaps[0]
2043 else:
2044 lyr = self.layer
2045 return (lyr.ztop, lyr.zbot)[self._direction_in == UP]
2047 def z_out(self, endgaps=None):
2048 if endgaps is not None:
2049 return endgaps[1]
2050 else:
2051 lyr = self.layer
2052 return (lyr.ztop, lyr.zbot)[self._direction_out == DOWN]
2054 def turns(self):
2055 return self._direction_in != self._direction_out
2057 def eff_direction_in(self, endgaps=None):
2058 if endgaps is None:
2059 return self._direction_in
2060 else:
2061 return endgaps[2]
2063 def eff_direction_out(self, endgaps=None):
2064 if endgaps is None:
2065 return self._direction_out
2066 else:
2067 return endgaps[3]
2069 def zturn(self, p):
2070 lyr = self.layer
2071 return lyr.zturn(p, self.mode)
2073 def u_in(self, endgaps=None):
2074 return self.layer.u(self.mode, self.z_in(endgaps))
2076 def u_out(self, endgaps=None):
2077 return self.layer.u(self.mode, self.z_out(endgaps))
2079 def critical_p_in(self, endgaps=None):
2080 z = self.z_in(endgaps)
2081 return self.layer.u(self.mode, z)*radius(z)
2083 def critical_p_out(self, endgaps=None):
2084 z = self.z_out(endgaps)
2085 return self.layer.u(self.mode, z)*radius(z)
2087 def xt(self, p, zpart=None):
2088 x, t = self.layer.xt(p, self.mode, zpart=zpart)
2089 if self._direction_in != self._direction_out and zpart is None:
2090 x *= 2.
2091 t *= 2.
2092 return x, t
2094 def xt_gap(self, p, zstart, zstop, samedir):
2095 z1, z2 = zstart, zstop
2096 if z1 > z2:
2097 z1, z2 = z2, z1
2099 x, t = self.layer.xt(p, self.mode, zpart=(z1, z2))
2101 if samedir:
2102 return x, t
2103 else:
2104 xfull, tfull = self.xt(p)
2105 return xfull-x, tfull-t
2107 def __hash__(self):
2108 return hash((
2109 self._direction_in,
2110 self._direction_out,
2111 self.mode,
2112 id(self.layer)))
2115class HeadwaveStraight(Straight):
2116 def __init__(self, direction_in, direction_out, mode, interface):
2117 Straight.__init__(self, direction_in, direction_out, mode, None)
2119 self.interface = interface
2121 def z_in(self, zpart=None):
2122 return self.interface.z
2124 def z_out(self, zpart=None):
2125 return self.interface.z
2127 def zturn(self, p):
2128 return filled(self.interface.z, len(p))
2130 def xt(self, p, zpart=None):
2131 return 0., 0.
2133 def x2t_headwave(self, xstretch):
2134 xstretch_m = xstretch*d2r*radius(self.interface.z)
2135 return min_not_none(*self.interface.u_top_bottom(self.mode))*xstretch_m
2138class Kink(RayElement):
2139 '''
2140 An interaction of a ray with a :py:class:`Discontinuity`.
2141 '''
2143 def __init__(
2144 self,
2145 in_direction,
2146 out_direction,
2147 in_mode,
2148 out_mode,
2149 discontinuity):
2151 self.in_direction = in_direction
2152 self.out_direction = out_direction
2153 self.in_mode = in_mode
2154 self.out_mode = out_mode
2155 self.discontinuity = discontinuity
2157 def reflection(self):
2158 return self.in_direction != self.out_direction
2160 def conversion(self):
2161 return self.in_mode != self.out_mode
2163 def efficiency(self, p, out_direction=None, out_mode=None):
2165 if out_direction is None:
2166 out_direction = self.out_direction
2168 if out_mode is None:
2169 out_mode = self.out_mode
2171 return self.discontinuity.efficiency(
2172 self.in_direction, out_direction, self.in_mode, out_mode, p)
2174 def __str__(self):
2175 r, c = self.reflection(), self.conversion()
2176 if r and c:
2177 return '|~'
2178 if r:
2179 return '|'
2180 if c:
2181 return '~'
2182 return '_'
2184 def __hash__(self):
2185 return hash((
2186 self.in_direction,
2187 self.out_direction,
2188 self.in_mode,
2189 self.out_mode,
2190 id(self.discontinuity)))
2193class PRangeNotSet(CakeError):
2194 pass
2197class RayPath(object):
2198 '''
2199 Representation of a fan of rays running through a common sequence of
2200 layers / interfaces.
2201 '''
2203 def __init__(self, phase):
2204 self.elements = []
2205 self.phase = phase
2206 self._pmax = None
2207 self._pmin = None
2208 self._p = None
2209 self._is_headwave = False
2211 def set_is_headwave(self, is_headwave):
2212 self._is_headwave = is_headwave
2214 def copy(self):
2215 '''
2216 Get a copy of it.
2217 '''
2219 c = copy.copy(self)
2220 c.elements = list(self.elements)
2221 return c
2223 def endgaps(self, zstart, zstop):
2224 '''
2225 Get information needed for end point adjustments.
2226 '''
2228 return (
2229 zstart,
2230 zstop,
2231 self.phase.direction_start(),
2232 self.phase.direction_stop())
2234 def append(self, element):
2235 self.elements.append(element)
2237 def _check_have_prange(self):
2238 if self._pmax is None:
2239 raise PRangeNotSet()
2241 def set_prange(self, pmin, pmax, dp):
2242 self._pmin, self._pmax = pmin, pmax
2243 self._prange_dp = dp
2245 def used_phase(self, p=None, eps=1.):
2246 '''
2247 Calculate phase definition from ray path.
2248 '''
2250 used = PhaseDef()
2251 fleg = self.phase.first_leg()
2252 used.append(Leg(fleg.departure, fleg.mode))
2253 n_elements_n = [None] + self.elements + [None]
2254 for before, element, after in zip(
2255 n_elements_n[:-2],
2256 n_elements_n[1:-1],
2257 n_elements_n[2:]):
2259 if element.is_kink() and HeadwaveStraight not in (
2260 type(before),
2261 type(after)):
2263 if element.reflection() or element.conversion():
2264 z = element.discontinuity.z
2265 used.append(Knee(
2266 z,
2267 element.in_direction,
2268 element.out_direction != element.in_direction,
2269 element.in_mode,
2270 element.out_mode))
2272 used.append(Leg(element.out_direction, element.out_mode))
2274 elif type(element) is HeadwaveStraight:
2275 z = element.interface.z
2276 k = Knee(
2277 z,
2278 before.in_direction,
2279 after.out_direction != before.in_direction,
2280 before.in_mode,
2281 after.out_mode)
2283 k.headwave = True
2284 used.append(k)
2285 used.append(Leg(after.out_direction, after.out_mode))
2287 if (p is not None and before and after
2288 and element.is_straight()
2289 and before.is_kink()
2290 and after.is_kink()
2291 and element.turns()
2292 and not before.reflection() and not before.conversion()
2293 and not after.reflection() and not after.conversion()):
2295 ai = element.angle_in(p)
2296 if 90.0-eps < ai and ai < 90+eps:
2297 used.append(
2298 Head(
2299 before.discontinuity.z,
2300 before.out_direction,
2301 element.mode))
2302 used.append(
2303 Leg(-before.out_direction, element.mode))
2305 used._direction_stop = self.phase.direction_stop()
2306 used._definition = self.phase.definition()
2308 return used
2310 def pmax(self):
2311 '''
2312 Get maximum valid ray parameter.
2313 '''
2314 self._check_have_prange()
2315 return self._pmax
2317 def pmin(self):
2318 '''
2319 Get minimum valid ray parameter.
2320 '''
2321 self._check_have_prange()
2322 return self._pmin
2324 def xmin(self):
2325 '''
2326 Get minimal distance.
2327 '''
2328 self._analyse()
2329 return self._xmin
2331 def xmax(self):
2332 '''
2333 Get maximal distance.
2334 '''
2335 self._analyse()
2336 return self._xmax
2338 def kinks(self):
2339 '''
2340 Iterate over propagation mode changes (reflections/transmissions).
2341 '''
2342 return (k for k in self.elements if isinstance(k, Kink))
2344 def straights(self):
2345 '''
2346 Iterate over ray segments.
2347 '''
2348 return (s for s in self.elements if isinstance(s, Straight))
2350 def headwave_straight(self):
2351 for s in self.elements:
2352 if type(s) is HeadwaveStraight:
2353 return s
2355 def first_straight(self):
2356 '''
2357 Get first ray segment.
2358 '''
2359 for s in self.elements:
2360 if isinstance(s, Straight):
2361 return s
2363 def last_straight(self):
2364 '''
2365 Get last ray segment.
2366 '''
2367 for s in reversed(self.elements):
2368 if isinstance(s, Straight):
2369 return s
2371 def efficiency(self, p):
2372 '''
2373 Get product of all conversion/reflection coefficients encountered on
2374 path.
2375 '''
2376 return reduce(
2377 operator.mul, (k.efficiency(p) for k in self.kinks()), 1.)
2379 def spreading(self, p, endgaps):
2380 '''
2381 Get geometrical spreading factor.
2382 '''
2383 if self._is_headwave:
2384 return 0.0
2386 self._check_have_prange()
2387 dp = self._prange_dp * 0.01
2388 assert self._pmax - self._pmin > dp
2390 if p + dp > self._pmax:
2391 p = p-dp
2393 x0, t = self.xt(p, endgaps)
2394 x1, t = self.xt(p+dp, endgaps)
2395 x0 *= d2r
2396 x1 *= d2r
2397 if x1 == x0:
2398 return num.nan
2400 dp_dx = dp/(x1-x0)
2402 x = x0
2403 if x == 0.:
2404 x = x1
2405 p = dp
2407 first = self.first_straight()
2408 last = self.last_straight()
2409 return num.abs(dp_dx) * first.pflat_in(p, endgaps) / (
2410 4.0 * math.pi * num.sin(x) *
2411 (earthradius-first.z_in(endgaps)) *
2412 (earthradius-last.z_out(endgaps))**2 *
2413 first.u_in(endgaps)**2 *
2414 num.abs(num.cos(first.angle_in(p, endgaps)*d2r)) *
2415 num.abs(num.cos(last.angle_out(p, endgaps)*d2r)))
2417 def make_p(self, dp=None, n=None, nmin=None):
2418 assert dp is None or n is None
2420 if self._pmin == self._pmax:
2421 return num.array([self._pmin])
2423 if dp is None:
2424 dp = self._prange_dp
2426 if n is None:
2427 n = int(round((self._pmax-self._pmin)/dp)) + 1
2429 if nmin is not None:
2430 n = max(n, nmin)
2432 ppp = num.linspace(self._pmin, self._pmax, n)
2433 return ppp
2435 def xt_endgaps(self, p, endgaps, which='both'):
2436 '''
2437 Get amount of distance/traveltime to be subtracted at the generic ray
2438 path's ends.
2439 '''
2441 zstart, zstop, dirstart, dirstop = endgaps
2442 firsts = self.first_straight()
2443 lasts = self.last_straight()
2444 xs, ts = firsts.xt_gap(
2445 p, zstart, firsts.z_in(), dirstart == firsts._direction_in)
2446 xe, te = lasts.xt_gap(
2447 p, zstop, lasts.z_out(), dirstop == lasts._direction_out)
2449 if which == 'both':
2450 return xs + xe, ts + te
2451 elif which == 'left':
2452 return xs, ts
2453 elif which == 'right':
2454 return xe, te
2456 def xt_endgaps_ptest(self, p, endgaps):
2457 '''
2458 Check if ray parameter is valid at source and receiver.
2459 '''
2461 zstart, zstop, dirstart, dirstop = endgaps
2462 firsts = self.first_straight()
2463 lasts = self.last_straight()
2464 return num.logical_and(firsts.test(p, zstart), lasts.test(p, zstop))
2466 def xt(self, p, endgaps):
2467 '''
2468 Calculate distance and traveltime for given ray parameter.
2469 '''
2471 if isinstance(p, num.ndarray):
2472 sx = num.zeros(p.size)
2473 st = num.zeros(p.size)
2474 else:
2475 sx = 0.0
2476 st = 0.0
2478 for s in self.straights():
2479 x, t = s.xt(p)
2480 sx += x
2481 st += t
2483 if endgaps:
2484 dx, dt = self.xt_endgaps(p, endgaps)
2485 sx -= dx
2486 st -= dt
2488 return sx, st
2490 def xt_limits(self, p):
2491 '''
2492 Calculate limits of distance and traveltime for given ray parameter.
2493 '''
2495 if isinstance(p, num.ndarray):
2496 sx = num.zeros(p.size)
2497 st = num.zeros(p.size)
2498 sxe = num.zeros(p.size)
2499 ste = num.zeros(p.size)
2500 else:
2501 sx = 0.0
2502 st = 0.0
2503 sxe = 0.0
2504 ste = 0.0
2506 sfirst = self.first_straight()
2507 slast = self.last_straight()
2509 for s in self.straights():
2510 if s is not sfirst and s is not slast:
2511 x, t = s.xt(p)
2512 sx += x
2513 st += t
2515 sends = [sfirst]
2516 if sfirst is not slast:
2517 sends.append(slast)
2519 for s in sends:
2520 x, t = s.xt(p)
2521 sxe += x
2522 ste += t
2524 return sx, (sx + sxe), st, (st + ste)
2526 def iter_zxt(self, p):
2527 '''
2528 Iterate over (depth, distance, traveltime) at each layer interface on
2529 ray path.
2530 '''
2532 sx = num.zeros(p.size)
2533 st = num.zeros(p.size)
2534 ok = False
2535 for s in self.straights():
2536 yield s.z_in(), sx.copy(), st.copy()
2538 x, t = s.xt(p)
2539 sx += x
2540 st += t
2541 ok = True
2543 if ok:
2544 yield s.z_out(), sx.copy(), st.copy()
2546 def zxt_path_subdivided(
2547 self, p, endgaps,
2548 points_per_straight=20,
2549 x_for_headwave=None):
2551 '''
2552 Get geometrical representation of ray path.
2553 '''
2555 if self._is_headwave:
2556 assert p.size == 1
2557 x, t = self.xt(p, endgaps)
2558 xstretch = x_for_headwave-x
2559 nout = xstretch.size
2560 else:
2561 nout = p.size
2563 dxl, dtl = self.xt_endgaps(p, endgaps, which='left')
2564 dxr, dtr = self.xt_endgaps(p, endgaps, which='right')
2566 # first create full path including the endgaps
2567 sx = num.zeros(nout) - dxl
2568 st = num.zeros(nout) - dtl
2569 zxt = []
2570 for s in self.straights():
2571 n = points_per_straight
2573 back = None
2574 zin, zout = s.z_in(), s.z_out()
2575 if type(s) is HeadwaveStraight:
2576 z = zin
2577 for i in range(n):
2578 xs = float(i)/(n-1) * xstretch
2579 ts = s.x2t_headwave(xs)
2580 zxt.append((filled(z, xstretch.size), sx+xs, st+ts))
2581 else:
2582 if zin != zout: # normal traversal
2583 zs = num.linspace(zin, zout, n).tolist()
2584 for z in zs:
2585 x, t = s.xt(p, zpart=sorted([zin, z]))
2586 zxt.append((filled(z, nout), sx + x, st + t))
2588 else: # ray turns in layer
2589 zturn = s.zturn(p)
2590 back = []
2591 for i in range(n):
2592 z = zin + (zturn - zin) * num.sin(
2593 float(i)/(n-1)*math.pi/2.0) * 0.999
2595 if zturn[0] >= zin:
2596 x, t = s.xt(p, zpart=[zin, z])
2597 else:
2598 x, t = s.xt(p, zpart=[z, zin])
2599 zxt.append((z, sx + x, st + t))
2600 back.append((z, x, t))
2602 if type(s) is HeadwaveStraight:
2603 x = xstretch
2604 t = s.x2t_headwave(xstretch)
2605 else:
2606 x, t = s.xt(p)
2608 sx += x
2609 st += t
2610 if back:
2611 for z, x, t in reversed(back):
2612 zxt.append((z, sx - x, st - t))
2614 # gather results as arrays with such that x[ip, ipoint]
2615 fanz, fanx, fant = [], [], []
2616 for z, x, t in zxt:
2617 fanz.append(z)
2618 fanx.append(x)
2619 fant.append(t)
2621 z = num.array(fanz).T
2622 x = num.array(fanx).T
2623 t = num.array(fant).T
2625 # cut off the endgaps, add exact endpoints
2626 xmax = x[:, -1] - dxr
2627 tmax = t[:, -1] - dtr
2628 zstart, zstop = endgaps[:2]
2629 zs, xs, ts = [], [], []
2630 for i in range(nout):
2631 t_ = t[i]
2632 indices = num.where(num.logical_and(0. <= t_, t_ <= tmax[i]))[0]
2633 n = indices.size + 2
2634 zs_, xs_, ts_ = [num.empty(n, dtype=float) for j in range(3)]
2635 zs_[1:-1] = z[i, indices]
2636 xs_[1:-1] = x[i, indices]
2637 ts_[1:-1] = t[i, indices]
2638 zs_[0], zs_[-1] = zstart, zstop
2639 xs_[0], xs_[-1] = 0., xmax[i]
2640 ts_[0], ts_[-1] = 0., tmax[i]
2641 zs.append(zs_)
2642 xs.append(xs_)
2643 ts.append(ts_)
2645 return zs, xs, ts
2647 def _analyse(self):
2648 if self._p is not None:
2649 return
2651 p = self.make_p(nmin=20)
2652 xmin, xmax, tmin, tmax = self.xt_limits(p)
2654 self._x, self._t, self._p = xmax, tmax, p
2655 self._xmin, self._xmax = xmin.min(), xmax.max()
2656 self._tmin, self._tmax = tmin.min(), tmax.max()
2658 def draft_pxt(self, endgaps):
2659 self._analyse()
2661 if not self._is_headwave:
2662 cp, cx, ct = self._p, self._x, self._t
2663 pcrit = min(
2664 self.critical_pstart(endgaps),
2665 self.critical_pstop(endgaps))
2667 if pcrit < self._pmin:
2668 empty = num.array([], dtype=float)
2669 return empty, empty, empty
2671 elif pcrit >= self._pmax:
2672 dx, dt = self.xt_endgaps(cp, endgaps)
2673 return cp, cx-dx, ct-dt
2675 else:
2676 n = num.searchsorted(cp, pcrit) + 1
2677 rp, rx, rt = num.empty((3, n), dtype=float)
2678 rp[:-1] = cp[:n-1]
2679 rx[:-1] = cx[:n-1]
2680 rt[:-1] = ct[:n-1]
2681 rp[-1] = pcrit
2682 rx[-1], rt[-1] = self.xt(pcrit, endgaps)
2683 dx, dt = self.xt_endgaps(rp, endgaps)
2684 rx[:-1] -= dx[:-1]
2685 rt[:-1] -= dt[:-1]
2686 return rp, rx, rt
2688 else:
2689 dx, dt = self.xt_endgaps(self._p, endgaps)
2690 p, x, t = self._p, self._x - dx, self._t - dt
2691 p, x, t = p[0], x[0], t[0]
2692 xh = num.linspace(0., x*10-x, 10)
2693 th = self.headwave_straight().x2t_headwave(xh)
2694 return filled(p, xh.size), x+xh, t+th
2696 def interpolate_x2pt_linear(self, x, endgaps):
2697 '''
2698 Get approximate ray parameter and traveltime for distance.
2699 '''
2701 self._analyse()
2703 if self._is_headwave:
2704 dx, dt = self.xt_endgaps(self._p, endgaps)
2705 xmin = self._x[0] - dx[0]
2706 tmin = self._t[0] - dt[0]
2707 el = self.headwave_straight()
2708 xok = x[x >= xmin]
2709 th = el.x2t_headwave(xstretch=(xok-xmin)) + tmin
2710 return [
2711 (x_, self._p[0], t, None) for (x_, t) in zip(xok, th)]
2713 else:
2714 if num.all(x < self._xmin) or num.all(self._xmax < x):
2715 return []
2717 rp, rx, rt = self.draft_pxt(endgaps)
2719 xp = interp(x, rx, rp, 0)
2720 xt = interp(x, rx, rt, 0)
2722 if (rp.size and
2723 len(xp) == 0 and
2724 rx[0] == 0.0 and
2725 any(x == 0.0) and
2726 rp[0] == 0.0):
2728 xp = [(0.0, rp[0])]
2729 xt = [(0.0, rt[0])]
2731 return [
2732 (x_, p, t, (rp, rx, rt)) for ((x_, p), (_, t)) in zip(xp, xt)]
2734 def __eq__(self, other):
2735 if len(self.elements) != len(other.elements):
2736 return False
2738 return all(a == b for a, b in zip(self.elements, other.elements))
2740 def __hash__(self):
2741 return hash(
2742 tuple(hash(x) for x in self.elements) +
2743 (self.phase.definition(), ))
2745 def __str__(self, p=None, eps=1.):
2746 x = []
2747 start_i = None
2748 end_i = None
2749 turn_i = None
2751 def append_layers(si, ei, ti):
2752 if si == ei and (ti is None or ti == si):
2753 x.append('%i' % si)
2754 else:
2755 if ti is not None:
2756 x.append('(%i-%i-%i)' % (si, ti, ei))
2757 else:
2758 x.append('(%i-%i)' % (si, ei))
2760 for el in self.elements:
2761 if type(el) is Straight:
2762 if start_i is None:
2763 start_i = el.layer.ilayer
2764 if el._direction_in != el._direction_out:
2765 turn_i = el.layer.ilayer
2766 end_i = el.layer.ilayer
2768 elif isinstance(el, Kink):
2769 if start_i is not None:
2770 append_layers(start_i, end_i, turn_i)
2771 start_i = None
2772 turn_i = None
2774 x.append(str(el))
2776 if start_i is not None:
2777 append_layers(start_i, end_i, turn_i)
2779 su = '(%s)' % self.used_phase(p=p, eps=eps).used_repr()
2781 return '%-15s %-17s %s' % (self.phase.definition(), su, ''.join(x))
2783 def critical_pstart(self, endgaps):
2784 '''
2785 Get critical ray parameter for source depth choice.
2786 '''
2788 return self.first_straight().critical_p_in(endgaps)
2790 def critical_pstop(self, endgaps):
2791 '''
2792 Get critical ray parameter for receiver depth choice.
2793 '''
2795 return self.last_straight().critical_p_out(endgaps)
2797 def ranges(self, endgaps):
2798 '''
2799 Get valid ranges of ray parameter, distance, and traveltime.
2800 '''
2801 p, x, t = self.draft_pxt(endgaps)
2802 return p.min(), p.max(), x.min(), x.max(), t.min(), t.max()
2804 def describe(self, endgaps=None, as_degrees=False):
2805 '''
2806 Get textual representation.
2807 '''
2809 self._analyse()
2811 if as_degrees:
2812 xunit = 'deg'
2813 xfact = 1.
2814 else:
2815 xunit = 'km'
2816 xfact = d2m/km
2818 sg = ''' Ranges for all depths in source and receiver layers:
2819 - x [%g, %g] %s
2820 - t [%g, %g] s
2821 - p [%g, %g] s/deg
2822''' % (
2823 self._xmin*xfact,
2824 self._xmax*xfact,
2825 xunit,
2826 self._tmin,
2827 self._tmax,
2828 self._pmin/r2d,
2829 self._pmax/r2d)
2831 if endgaps is not None:
2832 pmin, pmax, xmin, xmax, tmin, tmax = self.ranges(endgaps)
2833 ss = ''' Ranges for given source and receiver depths:
2834\n - x [%g, %g] %s
2835\n - t [%g, %g] s
2836\n - p [%g, %g] s/deg
2837\n''' % (xmin*xfact, xmax*xfact, xunit, tmin, tmax, pmin/r2d, pmax/r2d)
2839 else:
2840 ss = ''
2842 return '%s\n' % self + ss + sg
2845class RefineFailed(CakeError):
2846 pass
2849class Ray(object):
2850 '''
2851 Representation of a ray with a specific (path, ray parameter, distance,
2852 arrival time) choice.
2854 **Attributes:**
2856 .. py:attribute:: path
2858 :py:class:`RayPath` object containing complete propagation history.
2860 .. py:attribute:: p
2862 Ray parameter (spherical) [s/rad]
2864 .. py:attribute:: x
2866 Radial distance [deg]
2868 .. py:attribute:: t
2870 Traveltime [s]
2872 .. py:attribute:: endgaps
2874 Needed for source/receiver depth adjustments in many
2875 :py:class:`RayPath` methods.
2876 '''
2878 def __init__(self, path, p, x, t, endgaps, draft_pxt):
2879 self.path = path
2880 self.p = p
2881 self.x = x
2882 self.t = t
2883 self.endgaps = endgaps
2884 self.draft_pxt = draft_pxt
2886 def given_phase(self):
2887 '''
2888 Get phase definition which was used to create the ray.
2890 :returns: :py:class:`PhaseDef` object
2891 '''
2893 return self.path.phase
2895 def used_phase(self):
2896 '''
2897 Compute phase definition from propagation path.
2899 :returns: :py:class:`PhaseDef` object
2900 '''
2902 return self.path.used_phase(self.p)
2904 def refine(self):
2905 if self.path._is_headwave:
2906 return
2908 if self.t == 0.0 and self.p == 0.0 and self.x == 0.0:
2909 return
2911 cp, cx, ct = self.draft_pxt
2912 ip = num.searchsorted(cp, self.p)
2913 if not (0 < ip < cp.size):
2914 raise RefineFailed()
2916 pl, ph = cp[ip-1], cp[ip]
2917 p_to_t = {}
2918 i = [0]
2920 def f(p):
2921 i[0] += 1
2922 x, t = self.path.xt(p, self.endgaps)
2923 p_to_t[p] = t
2924 return self.x - x
2926 try:
2927 self.p = brentq(f, pl, ph)
2928 self.t = p_to_t[self.p]
2930 except ValueError:
2931 raise RefineFailed()
2933 def takeoff_angle(self):
2934 '''
2935 Get takeoff angle of ray.
2937 The angle is returned in [degrees].
2938 '''
2940 return self.path.first_straight().angle_in(self.p, self.endgaps)
2942 def incidence_angle(self):
2943 '''
2944 Get incidence angle of ray.
2946 The angle is returned in [degrees].
2947 '''
2949 return self.path.last_straight().angle_out(self.p, self.endgaps)
2951 def efficiency(self):
2952 '''
2953 Get conversion/reflection efficiency of the ray.
2955 A value between 0 and 1 is returned, reflecting the relative amount of
2956 energy which is transmitted along the ray and not lost by reflections
2957 or conversions.
2958 '''
2960 return self.path.efficiency(self.p)
2962 def spreading(self):
2963 '''
2964 Get geometrical spreading factor.
2965 '''
2967 return self.path.spreading(self.p, self.endgaps)
2969 def surface_sphere(self):
2970 x1, y1 = 0., earthradius - self.endgaps[0]
2971 r2 = earthradius - self.endgaps[1]
2972 x2, y2 = r2*math.sin(self.x*d2r), r2*math.cos(self.x*d2r)
2973 return ((x2-x1)**2 + (y2-y1)**2)*4.0*math.pi
2975 def zxt_path_subdivided(self, points_per_straight=20):
2976 '''
2977 Get geometrical representation of ray path.
2979 Three arrays (depth, distance, time) with points on the ray's path of
2980 propagation are returned. The number of points which are used in each
2981 ray segment (passage through one layer) may be controlled by the
2982 ``points_per_straight`` parameter.
2983 '''
2984 return self.path.zxt_path_subdivided(
2985 num.atleast_1d(self.p), self.endgaps,
2986 points_per_straight=points_per_straight,
2987 x_for_headwave=num.atleast_1d(self.x))
2989 def __str__(self, as_degrees=False):
2990 if as_degrees:
2991 sd = '%6.3g deg' % self.x
2992 else:
2993 sd = '%7.5g km' % (self.x*(d2r*earthradius/km))
2995 return '%7.5g s/deg %s %6.4g s %5.1f %5.1f %3.0f%% %3.0f%% %s' % (
2996 self.p/r2d,
2997 sd,
2998 self.t,
2999 self.takeoff_angle(),
3000 self.incidence_angle(),
3001 100*self.efficiency(),
3002 100*self.spreading()*self.surface_sphere(),
3003 self.path.__str__(p=self.p))
3006def anything_to_crust2_profile(crust2_profile):
3007 from pyrocko.dataset import crust2x2
3008 if isinstance(crust2_profile, tuple):
3009 lat, lon = [float(x) for x in crust2_profile]
3010 return crust2x2.get_profile(lat, lon)
3011 elif isinstance(crust2_profile, str):
3012 return crust2x2.get_profile(crust2_profile)
3013 elif isinstance(crust2_profile, crust2x2.Crust2Profile):
3014 return crust2_profile
3015 else:
3016 assert False, 'crust2_profile must be (lat, lon) a profile ' \
3017 'key or a crust2x2 Profile object)'
3020class DiscontinuityNotFound(CakeError):
3021 def __init__(self, depth_or_name):
3022 CakeError.__init__(self)
3023 self.depth_or_name = depth_or_name
3025 def __str__(self):
3026 return 'Cannot find discontinuity from given depth or name: %s' % \
3027 self.depth_or_name
3030class LayeredModelError(CakeError):
3031 pass
3034class LayeredModel(object):
3035 '''
3036 Representation of a layer cake model.
3038 There are several ways to initialize an instance of this class.
3040 1. Use the module function :py:func:`load_model` to read a model from a
3041 file.
3042 2. Create an empty model with the default constructor and append layers and
3043 discontinuities with the :py:meth:`append` method (from top to bottom).
3044 3. Use the constructor :py:meth:`LayeredModel.from_scanlines`, to
3045 automatically create the :py:class:`Layer` and :py:class:`Discontinuity`
3046 objects from a given velocity profile.
3048 An earth model is represented by as stack of :py:class:`Layer` and
3049 :py:class:`Discontinuity` objects. The method :py:meth:`arrivals` returns
3050 :py:class:`Ray` objects which may be e.g. queried for arrival times of
3051 specific phases. Each ray is associated with a :py:class:`RayPath` object.
3052 Ray objects share common ray paths if they have the same
3053 conversion/reflection/propagation history. Creating the ray path objects is
3054 relatively expensive (this is done in :py:meth:`gather_paths`), but they
3055 are cached for reuse in successive invocations.
3056 '''
3058 def __init__(self):
3059 self._surface_material = None
3060 self._elements = []
3061 self.nlayers = 0
3062 self._np = 10000
3063 self._pdepth = 5
3064 self._pathcache = {}
3066 def copy_with_elevation(self, elevation):
3067 '''
3068 Get copy of the model with surface layer stretched to given elevation.
3070 :param elevation: new surface elevation in [m]
3072 Elevation is positiv upward, contrary to the layered models downward
3073 `z` axis.
3074 '''
3076 c = copy.deepcopy(self)
3077 c._pathcache = {}
3078 surface = c._elements[0]
3079 toplayer = c._elements[1]
3081 assert toplayer.zbot > -elevation
3083 surface.z = -elevation
3084 c._elements[1] = toplayer.copy(ztop=-elevation)
3085 c._elements[1].ilayer = 0
3086 return c
3088 def zeq(self, z1, z2):
3089 return abs(z1-z2) < ZEPS
3091 def append(self, element):
3092 '''
3093 Add a layer or discontinuity at bottom of model.
3095 :param element: object of subclass of :py:class:`Layer` or
3096 :py:class:`Discontinuity`.
3097 '''
3099 if isinstance(element, Layer):
3100 if element.zbot >= earthradius:
3101 element.zbot = earthradius - 1.
3103 if element.ztop >= earthradius:
3104 raise CakeError('Layer deeper than earthradius')
3106 element.ilayer = self.nlayers
3107 self.nlayers += 1
3109 self._elements.append(element)
3111 def elements(self, direction=DOWN):
3112 '''
3113 Iterate over all elements of the model.
3115 :param direction: direction of traversal :py:const:`DOWN` or
3116 :py:const:`UP`.
3118 Objects derived from the :py:class:`Discontinuity` and
3119 :py:class:`Layer` classes are yielded.
3120 '''
3122 if direction == DOWN:
3123 return iter(self._elements)
3124 else:
3125 return reversed(self._elements)
3127 def layers(self, direction=DOWN):
3128 '''
3129 Iterate over all layers of model.
3131 :param direction: direction of traversal :py:const:`DOWN` or
3132 :py:const:`UP`.
3134 Objects derived from the :py:class:`Layer` class are yielded.
3135 '''
3137 if direction == DOWN:
3138 return (el for el in self._elements if isinstance(el, Layer))
3139 else:
3140 return (
3141 el for el in reversed(self._elements) if isinstance(el, Layer))
3143 def layer(self, z, direction=DOWN):
3144 '''
3145 Get layer for given depth.
3147 :param z: depth [m]
3148 :param direction: direction of traversal :py:const:`DOWN` or
3149 :py:const:`UP`.
3151 Returns first layer which touches depth ``z`` (tolerant at boundaries).
3152 '''
3154 for layer in self.layers(direction):
3155 if layer.contains(z):
3156 return layer
3157 else:
3158 raise CakeError('Failed extracting layer at depth z=%s' % z)
3160 def walker(self):
3161 return Walker(self._elements)
3163 def material(self, z, direction=DOWN):
3164 '''
3165 Get material at given depth.
3167 :param z: depth [m]
3168 :param direction: direction of traversal :py:const:`DOWN` or
3169 :py:const:`UP`
3170 :returns: object of type :py:class:`Material`
3172 If given depth ``z`` happens to be at an interface, the material of the
3173 first layer with respect to the the traversal ordering is returned.
3174 '''
3176 lyr = self.layer(z, direction)
3177 return lyr.material(z)
3179 def discontinuities(self):
3180 '''
3181 Iterate over all discontinuities of the model.'''
3183 return (el for el in self._elements if isinstance(el, Discontinuity))
3185 def discontinuity(self, name_or_z):
3186 '''
3187 Get discontinuity by name or depth.
3189 :param name_or_z: name of discontinuity or depth [m] as float value
3190 '''
3192 if isinstance(name_or_z, float):
3193 candi = sorted(
3194 self.discontinuities(), key=lambda i: abs(i.z-name_or_z))
3195 else:
3196 candi = [i for i in self.discontinuities() if i.name == name_or_z]
3198 if not candi:
3199 raise DiscontinuityNotFound(name_or_z)
3201 return candi[0]
3203 def adapt_phase(self, phase):
3204 '''
3205 Adapt a phase definition for use with this model.
3207 This returns a copy of the phase definition, where named
3208 discontinuities are replaced with the actual depth of these, as defined
3209 in the model.
3210 '''
3212 phase = phase.copy()
3213 for knee in phase.knees():
3214 if knee.depth != 'surface':
3215 knee.depth = self.discontinuity(knee.depth).z
3216 for leg in phase.legs():
3217 if leg.depthmax is not None and isinstance(leg.depthmax, str):
3218 leg.depthmax = self.discontinuity(leg.depthmax).z
3220 return phase
3222 def path(self, p, phase, layer_start, layer_stop):
3223 '''
3224 Get ray path for given combination of ray parameter, phase definition,
3225 source and receiver layers.
3227 :param p: ray parameter (spherical) [s/rad]
3228 :param phase: phase definition (:py:class:`PhaseDef` object)
3229 :param layer_start: layer with source
3230 :param layer_stop: layer with receiver
3231 :returns: :py:class:`RayPath` object
3233 If it is not possible to find a solution, an exception of type
3234 :py:exc:`NotPhaseConform`, :py:exc:`MinDepthReached`,
3235 :py:exc:`MaxDepthReached`, :py:exc:`CannotPropagate`,
3236 :py:exc:`BottomReached` or :py:exc:`SurfaceReached` is raised.
3237 '''
3239 phase = self.adapt_phase(phase)
3240 knees = phase.knees()
3241 legs = phase.legs()
3242 next_knee = next_or_none(knees)
3243 leg = next_or_none(legs)
3244 assert leg is not None
3246 direction = leg.departure
3247 direction_stop = phase.direction_stop()
3248 mode = leg.mode
3249 mode_stop = phase.last_leg().mode
3251 walker = self.walker()
3252 walker.goto_layer(layer_start)
3253 current = walker.current()
3255 ttop, tbot = current.tests(p, mode)
3256 if not ttop and not tbot:
3257 raise CannotPropagate(direction, current.ilayer)
3259 if (direction == DOWN and not ttop) or (direction == UP and not tbot):
3260 direction = -direction
3262 path = RayPath(phase)
3263 trapdetect = set()
3264 while True:
3265 at_layer = isinstance(current, Layer)
3266 at_discontinuity = isinstance(current, Discontinuity)
3268 # detect trapped wave
3269 k = (id(next_knee), id(current), direction, mode)
3270 if k in trapdetect:
3271 raise Trapped()
3273 trapdetect.add(k)
3275 if at_discontinuity:
3276 oldmode, olddirection = mode, direction
3277 headwave = False
3278 if next_knee is not None and next_knee.matches(
3279 current, mode, direction):
3281 headwave = next_knee.headwave
3282 direction = next_knee.out_direction()
3283 mode = next_knee.out_mode
3284 next_knee = next_or_none(knees)
3285 leg = next(legs)
3287 else: # implicit reflection/transmission
3288 direction = current.propagate(p, mode, direction)
3290 if headwave:
3291 path.set_is_headwave(True)
3293 path.append(Kink(
3294 olddirection, olddirection, oldmode, oldmode, current))
3296 path.append(HeadwaveStraight(
3297 olddirection, direction, oldmode, current))
3299 path.append(Kink(
3300 olddirection, direction, oldmode, mode, current))
3302 else:
3303 path.append(Kink(
3304 olddirection, direction, oldmode, mode, current))
3306 if at_layer:
3307 direction_in = direction
3308 direction = current.propagate(p, mode, direction_in)
3310 zturn = None
3311 if direction_in != direction:
3312 zturn = current.zturn(p, mode)
3314 zmin, zmax = leg.depthmin, leg.depthmax
3315 if zmin is not None or zmax is not None:
3316 if direction_in != direction:
3317 if zmin is not None and zturn <= zmin:
3318 raise MinDepthReached()
3319 if zmax is not None and zturn >= zmax:
3320 raise MaxDepthReached()
3321 else:
3322 if zmin is not None and current.ztop <= zmin:
3323 raise MinDepthReached()
3324 if zmax is not None and current.zbot >= zmax:
3325 raise MaxDepthReached()
3327 path.append(Straight(direction_in, direction, mode, current))
3329 if next_knee is None and mode == mode_stop and \
3330 current is layer_stop:
3332 if zturn is None:
3333 if direction == direction_stop:
3334 break
3335 else:
3336 break
3338 walker.go(direction)
3339 current = walker.current()
3341 return path
3343 def gather_paths(self, phases=PhaseDef('P'), zstart=0.0, zstop=0.0):
3344 '''
3345 Get all possible ray paths for given source and receiver depths for one
3346 or more phase definitions.
3348 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3349 Comma-separated strings and lists of such strings are also accepted
3350 and are converted to :py:class:`PhaseDef` objects for convenience.
3351 :param zstart: source depth [m]
3352 :param zstop: receiver depth [m]
3353 :returns: a list of :py:class:`RayPath` objects
3355 Results of this method are cached internally. Cached results are
3356 returned, when a given combination of source layer, receiver layer and
3357 phase definition has been used before.
3358 '''
3360 eps = 1e-7 # num.finfo(float).eps * 1000.
3362 phases = to_phase_defs(phases)
3364 paths = []
3365 for phase in phases:
3367 layer_start = self.layer(zstart, -phase.direction_start())
3368 layer_stop = self.layer(zstop, phase.direction_stop())
3370 pathcachekey = (phase.definition(), layer_start, layer_stop)
3372 if pathcachekey in self._pathcache:
3373 phase_paths = self._pathcache[pathcachekey]
3374 else:
3375 hwknee = phase.headwave_knee()
3376 if hwknee:
3377 name_or_z = hwknee.depth
3378 interface = self.discontinuity(name_or_z)
3379 mode = hwknee.in_mode
3380 in_direction = hwknee.direction
3382 pabove, pbelow = interface.critical_ps(mode)
3384 p = min_not_none(pabove, pbelow)
3386 # diffracted wave:
3387 if in_direction == DOWN and (
3388 pbelow is None or pbelow >= pabove):
3390 p *= (1.0 - eps)
3392 path = self.path(p, phase, layer_start, layer_stop)
3393 path.set_prange(p, p, 1.)
3395 phase_paths = [path]
3397 else:
3398 try:
3399 pmax_start = max([
3400 radius(z)/layer_start.v(phase.first_leg().mode, z)
3401 for z in (layer_start.ztop, layer_start.zbot)])
3403 pmax_stop = max([
3404 radius(z)/layer_stop.v(phase.last_leg().mode, z)
3405 for z in (layer_stop.ztop, layer_stop.zbot)])
3407 pmax = min(pmax_start, pmax_stop)
3409 pedges = [0.]
3410 for layer in self.layers():
3411 for z in (layer.ztop, layer.zbot):
3412 for mode in (P, S):
3413 for eps2 in [eps]:
3414 v = layer.v(mode, z)
3415 if v != 0.0:
3416 p = radius(z)/v
3417 if p <= pmax:
3418 pedges.append(p*(1.0-eps2))
3419 pedges.append(p)
3420 pedges.append(p*(1.0+eps2))
3422 pedges = num.unique(sorted(pedges))
3424 phase_paths = {}
3425 cached = {}
3426 counter = [0]
3428 def p_to_path(p):
3429 if p in cached:
3430 return cached[p]
3432 try:
3433 counter[0] += 1
3434 path = self.path(
3435 p, phase, layer_start, layer_stop)
3437 if path not in phase_paths:
3438 phase_paths[path] = []
3440 phase_paths[path].append(p)
3442 except PathFailed:
3443 path = None
3445 cached[p] = path
3446 return path
3448 def recurse(pmin, pmax, i=0):
3449 if i > self._pdepth:
3450 return
3451 path1 = p_to_path(pmin)
3452 path2 = p_to_path(pmax)
3453 if path1 is None and path2 is None and i > 0:
3454 return
3455 if path1 is None or path2 is None or \
3456 hash(path1) != hash(path2):
3458 recurse(pmin, (pmin+pmax)/2., i+1)
3459 recurse((pmin+pmax)/2., pmax, i+1)
3461 for (pl, ph) in zip(pedges[:-1], pedges[1:]):
3462 recurse(pl, ph)
3464 for path, ps in phase_paths.items():
3465 path.set_prange(
3466 min(ps), max(ps), pmax/(self._np-1))
3468 phase_paths = list(phase_paths.keys())
3470 except ZeroDivisionError:
3471 phase_paths = []
3473 self._pathcache[pathcachekey] = phase_paths
3475 paths.extend(phase_paths)
3477 paths.sort(key=lambda x: x.pmin())
3478 return paths
3480 def arrivals(
3481 self,
3482 distances=[],
3483 phases=PhaseDef('P'),
3484 zstart=0.0,
3485 zstop=0.0,
3486 refine=True):
3488 '''
3489 Compute rays and traveltimes for given distances.
3491 :param distances: list or array of distances [deg]
3492 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3493 Comma-separated strings and lists of such strings are also accepted
3494 and are converted to :py:class:`PhaseDef` objects for convenience.
3495 :param zstart: source depth [m]
3496 :param zstop: receiver depth [m]
3497 :param refine: bool flag, whether to use bisectioning to improve
3498 (p, x, t) estimated from interpolation
3499 :returns: a list of :py:class:`Ray` objects, sorted by
3500 (distance, arrival time)
3501 '''
3503 distances = num.asarray(distances, dtype=float)
3505 arrivals = []
3506 for path in self.gather_paths(phases, zstart=zstart, zstop=zstop):
3508 endgaps = path.endgaps(zstart, zstop)
3509 for x, p, t, draft_pxt in path.interpolate_x2pt_linear(
3510 distances, endgaps):
3512 arrivals.append(Ray(path, p, x, t, endgaps, draft_pxt))
3514 if refine:
3515 refined = []
3516 for ray in arrivals:
3518 if ray.path._is_headwave:
3519 refined.append(ray)
3521 try:
3522 ray.refine()
3523 refined.append(ray)
3525 except RefineFailed:
3526 pass
3528 arrivals = refined
3530 arrivals.sort(key=lambda x: (x.x, x.t))
3531 return arrivals
3533 @classmethod
3534 def from_scanlines(cls, producer):
3535 '''
3536 Create layer cake model from sequence of materials at depths.
3538 :param producer: iterable yielding (depth, material, name) tuples
3540 Creates a new :py:class:`LayeredModel` object and uses its
3541 :py:meth:`append` method to add layers and discontinuities as needed.
3542 '''
3544 self = cls()
3545 for z, material, name in producer:
3547 if not self._elements:
3548 self.append(Surface(z, material))
3549 else:
3550 element = self._elements[-1]
3551 if self.zeq(element.zbot, z):
3552 assert isinstance(element, Layer)
3553 self.append(
3554 Interface(z, element.mbot, material, name=name))
3556 else:
3557 if isinstance(element, Discontinuity):
3558 ztop = element.z
3559 mtop = element.mbelow
3560 elif isinstance(element, Layer):
3561 ztop = element.zbot
3562 mtop = element.mbot
3564 if mtop == material:
3565 layer = HomogeneousLayer(
3566 ztop, z, material, name=name)
3567 else:
3568 layer = GradientLayer(
3569 ztop, z, mtop, material, name=name)
3571 self.append(layer)
3573 return self
3575 def to_scanlines(self, get_burgers=False):
3576 def fmt(z, m):
3577 if not m._has_default_burgers() or get_burgers:
3578 return (z, m.vp, m.vs, m.rho, m.qp, m.qs,
3579 m.burger_eta1, m.burger_eta2, m.burger_valpha)
3580 return (z, m.vp, m.vs, m.rho, m.qp, m.qs)
3582 last = None
3583 lines = []
3584 for element in self.elements():
3585 if isinstance(element, Layer):
3586 if not isinstance(last, Layer):
3587 lines.append(fmt(element.ztop, element.mtop))
3589 lines.append(fmt(element.zbot, element.mbot))
3591 last = element
3593 if not isinstance(last, Layer):
3594 lines.append(fmt(last.z, last.mbelow))
3596 return lines
3598 def iter_material_parameter(self, get):
3599 assert get in ('vp', 'vs', 'rho', 'qp', 'qs', 'z')
3600 if get == 'z':
3601 for layer in self.layers():
3602 yield layer.ztop
3603 yield layer.zbot
3604 else:
3605 getter = operator.attrgetter(get)
3606 for layer in self.layers():
3607 yield getter(layer.mtop)
3608 yield getter(layer.mbot)
3610 def profile(self, get):
3611 '''
3612 Get parameter profile along depth of the earthmodel.
3614 :param get: property to be queried (
3615 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3616 :type get: string
3617 '''
3619 return num.array(list(self.iter_material_parameter(get)))
3621 def min(self, get='vp'):
3622 '''
3623 Find minimum value of a material property or depth.
3625 :param get: property to be queried (
3626 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3627 '''
3629 return min(self.iter_material_parameter(get))
3631 def max(self, get='vp'):
3632 '''
3633 Find maximum value of a material property or depth.
3635 :param get: property to be queried (
3636 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, ``'qs'``, or ``'z'``)
3637 '''
3639 return max(self.iter_material_parameter(get))
3641 def simplify_layers(self, layers, max_rel_error=0.001):
3642 if len(layers) <= 1:
3643 return layers
3645 ztop = layers[0].ztop
3646 zbot = layers[-1].zbot
3647 zorigs = [layer.ztop for layer in layers]
3648 zorigs.append(zbot)
3649 zs = num.linspace(ztop, zbot, 100)
3650 data = []
3651 for z in zs:
3652 if z == ztop:
3653 direction = UP
3654 else:
3655 direction = DOWN
3657 mat = self.material(z, direction)
3658 data.append(mat.astuple())
3660 data = num.array(data, dtype=float)
3661 data_means = num.mean(data, axis=0)
3662 nmax = len(layers) // 2
3663 accept = False
3665 zcut_best = []
3666 for n in range(1, nmax+1):
3667 ncutintervals = 20
3668 zdelta = (zbot-ztop)/ncutintervals
3669 if n == 2:
3670 zcuts = [
3671 [ztop, ztop + i*zdelta, zbot]
3672 for i in range(1, ncutintervals)]
3673 elif n == 3:
3674 zcuts = []
3675 for j in range(1, ncutintervals):
3676 for i in range(j+1, ncutintervals):
3677 zcuts.append(
3678 [ztop, ztop + j*zdelta, ztop + i*zdelta, zbot])
3679 else:
3680 zcuts = []
3681 zcuts.append(num.linspace(ztop, zbot, n+1))
3682 if zcut_best:
3683 zcuts.append(sorted(num.linspace(
3684 ztop, zbot, n).tolist() + zcut_best[1]))
3685 zcuts.append(sorted(num.linspace(
3686 ztop, zbot, n-1).tolist() + zcut_best[2]))
3688 best = None
3689 for icut, zcut in enumerate(zcuts):
3690 rel_par_errors = num.zeros(5)
3691 mpar_nodes = num.zeros((n+1, 5))
3693 for ipar in range(5):
3694 znodes, vnodes, error_rms = util.polylinefit(
3695 zs, data[:, ipar], zcut)
3697 mpar_nodes[:, ipar] = vnodes
3698 if data_means[ipar] == 0.0:
3699 rel_par_errors[ipar] = -1
3700 else:
3701 rel_par_errors[ipar] = error_rms/data_means[ipar]
3703 rel_error = rel_par_errors.max()
3704 if best is None or rel_error < best[0]:
3705 best = (rel_error, zcut, mpar_nodes)
3707 rel_error, zcut, mpar_nodes = best
3709 zcut_best.append(list(zcut))
3710 zcut_best[-1].pop(0)
3711 zcut_best[-1].pop()
3713 if rel_error <= max_rel_error:
3714 accept = True
3715 break
3717 if not accept:
3718 return layers
3720 rel_error, zcut, mpar_nodes = best
3722 material_nodes = []
3723 for i in range(n+1):
3724 material_nodes.append(Material(*mpar_nodes[i, :]))
3726 out_layers = []
3727 for i in range(n):
3728 mtop = material_nodes[i]
3729 mbot = material_nodes[i+1]
3730 ztop = zcut[i]
3731 zbot = zcut[i+1]
3732 if mtop == mbot:
3733 lyr = HomogeneousLayer(ztop, zbot, mtop)
3734 else:
3735 lyr = GradientLayer(ztop, zbot, mtop, mbot)
3737 out_layers.append(lyr)
3738 return out_layers
3740 def simplify(self, max_rel_error=0.001):
3741 '''
3742 Get representation of model with lower resolution.
3744 Returns an approximation of the model. All discontinuities are kept,
3745 but layer stacks with continuous model parameters are represented, if
3746 possible, by a lower number of layers. Piecewise linear functions are
3747 fitted against the original model parameter's piecewise linear
3748 functions. Successively larger numbers of layers are tried, until the
3749 difference to the original model is below ``max_rel_error``. The
3750 difference is measured as the RMS error of the fit normalized by the
3751 mean of the input (i.e. the fitted curves should deviate, on average,
3752 less than 0.1% from the input curves if ``max_rel_error`` = 0.001).
3753 '''
3755 mod_simple = LayeredModel()
3757 glayers = []
3758 for element in self.elements():
3760 if isinstance(element, Discontinuity):
3761 for layer in self.simplify_layers(
3762 glayers, max_rel_error=max_rel_error):
3764 mod_simple.append(layer)
3766 glayers = []
3767 mod_simple.append(element)
3768 else:
3769 glayers.append(element)
3771 for layer in self.simplify_layers(
3772 glayers, max_rel_error=max_rel_error):
3773 mod_simple.append(layer)
3775 return mod_simple
3777 def extract(self, depth_min=None, depth_max=None):
3778 '''
3779 Extract :py:class:`LayeredModel` from :py:class:`LayeredModel`.
3781 :param depth_min: depth of upper cut or name of :py:class:`Interface`
3782 :param depth_max: depth of lower cut or name of :py:class:`Interface`
3784 Interpolates a :py:class:`GradientLayer` at ``depth_min`` and/or
3785 ``depth_max``.
3786 '''
3788 if isinstance(depth_min, str):
3789 depth_min = self.discontinuity(depth_min).z
3791 if isinstance(depth_max, str):
3792 depth_max = self.discontinuity(depth_max).z
3794 mod_extracted = LayeredModel()
3796 for element in self.elements():
3797 element = element.copy()
3798 do_append = False
3799 if (depth_min is None or depth_min <= element.ztop) \
3800 and (depth_max is None or depth_max >= element.zbot):
3801 mod_extracted.append(element)
3802 continue
3804 if depth_min is not None:
3805 if element.ztop < depth_min and depth_min < element.zbot:
3806 _, element = element.split(depth_min)
3807 do_append = True
3809 if depth_max is not None:
3810 if element.zbot > depth_max and depth_max > element.ztop:
3811 element, _ = element.split(depth_max)
3812 do_append = True
3814 if do_append:
3815 mod_extracted.append(element)
3817 return mod_extracted
3819 def replaced_crust(self, crust2_profile=None, crustmod=None):
3820 if crust2_profile is not None:
3821 profile = anything_to_crust2_profile(crust2_profile)
3822 crustmod = LayeredModel.from_scanlines(
3823 from_crust2x2_profile(profile))
3825 newmod = LayeredModel()
3826 for element in crustmod.extract(depth_max='moho').elements():
3827 if element.name != 'moho':
3828 newmod.append(element)
3829 else:
3830 moho1 = element
3832 mod = self.extract(depth_min='moho')
3833 first = True
3834 for element in mod.elements():
3835 if element.name == 'moho':
3836 if element.z <= moho1.z:
3837 mbelow = mod.material(moho1.z, direction=UP)
3838 else:
3839 mbelow = element.mbelow
3841 moho = Interface(moho1.z, moho1.mabove, mbelow, name='moho')
3842 newmod.append(moho)
3843 else:
3844 if first:
3845 if isinstance(element, Layer) and element.zbot > moho.z:
3846 newmod.append(GradientLayer(
3847 moho.z,
3848 element.zbot,
3849 moho.mbelow,
3850 element.mbot,
3851 name=element.name))
3853 first = False
3854 else:
3855 newmod.append(element)
3856 return newmod
3858 def perturb(self, rstate=None, keep_vp_vs=False, **kwargs):
3859 '''
3860 Create a perturbed variant of the earth model.
3862 Randomly change the thickness and material parameters of the earth
3863 model from a uniform distribution.
3865 :param kwargs: Maximum change fraction (e.g. 0.1) of the parameters.
3866 Name the parameter, prefixed by ``p``. Supported parameters are
3867 ``ph, pvp, pvs, prho, pqs, pqp``.
3868 :type kwargs: dict
3869 :param rstate: Random state to draw from, defaults to ``None``
3870 :type rstate: :class:`numpy.random.RandomState`, optional
3871 :param keep_vp_vs: Keep the Vp/Vs ratio, defaults to False
3872 :type keep_vp_vs: bool, optional
3874 :returns: A new, perturbed earth model
3875 :rtype: :class:`~pyrocko.cake.LayeredModel`
3877 .. code-block :: python
3879 perturbed_model = model.perturb(ph=.1, pvp=.05, prho=.1)
3880 '''
3881 _pargs = set(['ph', 'pvp', 'pvs', 'prho', 'pqs', 'pqp'])
3882 earthmod = copy.deepcopy(self)
3884 if rstate is None:
3885 rstate = num.random.RandomState()
3887 layers = earthmod.layers()
3888 discont = earthmod.discontinuities()
3889 prev_layer = None
3891 def get_change_ratios():
3892 values = dict.fromkeys([p[1:] for p in _pargs], 0.)
3894 for param, pval in kwargs.items():
3895 if param not in _pargs:
3896 continue
3897 values[param[1:]] = float(rstate.uniform(-pval, pval, size=1))
3898 return values
3900 # skip Surface
3901 while True:
3902 disc = next(discont)
3903 if isinstance(disc, Surface):
3904 break
3906 while True:
3907 try:
3908 layer = next(layers)
3909 m = layer.material(None)
3910 h = layer.zbot - layer.ztop
3911 except StopIteration:
3912 break
3914 if not isinstance(layer, HomogeneousLayer):
3915 raise NotImplementedError(
3916 'Can only perturbate homogeneous layers!')
3918 changes = get_change_ratios()
3920 # Changing thickness
3921 dh = h * changes['h']
3922 changes['h'] = dh
3924 layer.resize(depth_max=layer.zbot + dh,
3925 depth_min=prev_layer.zbot if prev_layer else None)
3927 try:
3928 disc = next(discont)
3929 disc.change_depth(disc.z + dh)
3930 except StopIteration:
3931 pass
3933 # Setting material parameters
3934 for param, change_ratio in changes.items():
3935 if param == 'h':
3936 continue
3938 value = m.__getattribute__(param)
3939 changes[param] = value * change_ratio
3941 if keep_vp_vs and changes['vp'] != 0.:
3942 changes['vs'] = (m.vp + changes['vp']) / m.vp_vs_ratio() - m.vs
3944 for param, change in changes.items():
3945 if param == 'h':
3946 continue
3947 value = m.__getattribute__(param)
3948 m.__setattr__(param, value + change)
3950 logger.info(
3951 'perturbating earthmodel: {}'.format(
3952 ' '.join(['{param}: {change:{len}.2f}'.format(
3953 param=p, change=c, len=8)
3954 for p, c in changes.items()])))
3956 prev_layer = layer
3958 return earthmod
3960 def require_homogeneous(self):
3961 elements = list(self.elements())
3963 if len(elements) != 2:
3964 raise LayeredModelError(
3965 'Homogeneous model required but found more than one layer in '
3966 'earthmodel.')
3968 if not isinstance(elements[1], HomogeneousLayer):
3969 raise LayeredModelError(
3970 'Homogeneous model required but element #1 is not of type '
3971 'HomogeneousLayer.')
3973 return elements[1].m
3975 def __str__(self):
3976 return '\n'.join(str(element) for element in self._elements)
3979def read_hyposat_model(fn):
3980 '''
3981 Reader for HYPOSAT earth model files.
3983 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
3985 Interface names are translated as follows: ``'MOHO'`` -> ``'moho'``,
3986 ``'CONR'`` -> ``'conrad'``
3987 '''
3989 with open(fn, 'r') as f:
3990 translate = {'MOHO': 'moho', 'CONR': 'conrad'}
3991 lname = None
3992 for iline, line in enumerate(f):
3993 if iline == 0:
3994 continue
3996 z, vp, vs, name = util.unpack_fixed('f10, f10, f10, a4', line)
3997 if not name:
3998 name = None
3999 material = Material(vp*1000., vs*1000.)
4001 tname = translate.get(lname, lname)
4002 yield z*1000., material, tname
4004 lname = name
4007def read_nd_model(fn):
4008 '''
4009 Reader for TauP style '.nd' (named discontinuity) files.
4011 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
4013 Interface names are translated as follows: ``'mantle'`` -> ``'moho'``,
4014 ``'outer-core'`` -> ``'cmb'``, ``'inner-core'`` -> ``'icb'``.
4016 The format has been modified to include Burgers materials parameters in
4017 columns 7 (burger_eta1), 8 (burger_eta2) and 9. eta(3).
4018 '''
4019 with open(fn, 'r') as f:
4020 for x in read_nd_model_fh(f):
4021 yield x
4024def read_nd_model_str(s):
4025 f = StringIO(s)
4026 for x in read_nd_model_fh(f):
4027 yield x
4028 f.close()
4031def read_nd_model_fh(f):
4032 translate = {'mantle': 'moho', 'outer-core': 'cmb', 'inner-core': 'icb'}
4033 name = None
4034 for line in f:
4035 toks = line.split()
4036 if len(toks) == 9 or len(toks) == 6 or len(toks) == 4:
4037 z, vp, vs, rho = [float(x) for x in toks[:4]]
4038 qp, qs = None, None
4039 burgers = None
4040 if len(toks) == 6 or len(toks) == 9:
4041 qp, qs = [float(x) for x in toks[4:6]]
4042 if len(toks) == 9:
4043 burgers = \
4044 [float(x) for x in toks[6:]]
4046 material = Material(
4047 vp*1000., vs*1000., rho*1000., qp, qs,
4048 burgers=burgers)
4050 yield z*1000., material, name
4051 name = None
4052 elif len(toks) == 1:
4053 name = translate.get(toks[0], toks[0])
4055 f.close()
4058def from_crust2x2_profile(profile, depthmantle=50000):
4059 from pyrocko.dataset import crust2x2
4061 default_qp_qs = {
4062 'soft sed.': (50., 50.),
4063 'hard sed.': (200., 200.),
4064 'upper crust': (600., 400.),
4065 }
4067 z = 0.
4068 for i in range(8):
4069 dz, vp, vs, rho = profile.get_layer(i)
4070 name = crust2x2.Crust2Profile.layer_names[i]
4071 if name in default_qp_qs:
4072 qp, qs = default_qp_qs[name]
4073 else:
4074 qp, qs = None, None
4076 material = Material(vp, vs, rho, qp, qs)
4077 iname = None
4078 if i == 7:
4079 iname = 'moho'
4080 if dz != 0.0:
4081 yield z, material, iname
4082 if i != 7:
4083 yield z+dz, material, name
4084 else:
4085 yield z+depthmantle, material, name
4087 z += dz
4090def write_nd_model_fh(mod, fh):
4091 def fmt(z, mat):
4092 rstr = ' '.join(
4093 util.gform(x, 4)
4094 for x in (
4095 z/1000.,
4096 mat.vp/1000.,
4097 mat.vs/1000.,
4098 mat.rho/1000.,
4099 mat.qp, mat.qs))
4100 if not mat._has_default_burgers():
4101 rstr += ' '.join(
4102 util.gform(x, 4)
4103 for x in (
4104 mat.burger_eta1,
4105 mat.burger_eta2,
4106 mat.burger_valpha))
4107 return rstr.rstrip() + '\n'
4109 translate = {
4110 'moho': 'mantle',
4111 'cmb': 'outer-core',
4112 'icb': 'inner-core'}
4114 last = None
4115 for element in mod.elements():
4116 if isinstance(element, Interface):
4117 if element.name is not None:
4118 n = translate.get(element.name, element.name)
4119 fh.write('%s\n' % n)
4121 elif isinstance(element, Layer):
4122 if not isinstance(last, Layer):
4123 fh.write(fmt(element.ztop, element.mtop))
4125 fh.write(fmt(element.zbot, element.mbot))
4127 last = element
4129 if not isinstance(last, Layer):
4130 fh.write(fmt(last.z, last.mbelow))
4133def write_nd_model_str(mod):
4134 f = StringIO()
4135 write_nd_model_fh(mod, f)
4136 return f.getvalue()
4139def write_nd_model(mod, fn):
4140 with open(fn, 'w') as f:
4141 write_nd_model_fh(mod, f)
4144def builtin_models():
4145 return sorted([
4146 os.path.splitext(os.path.basename(x))[0]
4147 for x in glob.glob(builtin_model_filename('*'))])
4150def builtin_model_filename(modelname):
4151 return util.data_file(os.path.join('earthmodels', modelname+'.nd'))
4154def load_model(fn='ak135-f-continental.m', format='nd', crust2_profile=None):
4155 '''
4156 Load layered earth model from file.
4158 :param fn: filename
4159 :param format: format
4160 :param crust2_profile: ``(lat, lon)`` or
4161 :py:class:`pyrocko.crust2x2.Crust2Profile` object, merge model with
4162 crustal profile. If ``fn`` is forced to be ``None`` only the converted
4163 CRUST2.0 profile is returned.
4164 :returns: object of type :py:class:`LayeredModel`
4166 The following formats are currently supported:
4168 ============== ===========================================================
4169 format description
4170 ============== ===========================================================
4171 ``'nd'`` 'named discontinuity' format used by the TauP programs
4172 ``'hyposat'`` format used by the HYPOSAT location program
4173 ============== ===========================================================
4175 The naming of interfaces is translated from the file format's native naming
4176 to Cake's own convention (See :py:func:`read_nd_model` and
4177 :py:func:`read_hyposat_model` for details). Cake likes the following
4178 internal names: ``'conrad'``, ``'moho'``, ``'cmb'`` (core-mantle boundary),
4179 ``'icb'`` (inner core boundary).
4180 '''
4182 if fn is not None:
4183 if format == 'nd':
4184 if not os.path.exists(fn) and fn in builtin_models():
4185 fn = builtin_model_filename(fn)
4186 reader = read_nd_model(fn)
4187 elif format == 'hyposat':
4188 reader = read_hyposat_model(fn)
4189 else:
4190 assert False, 'unsupported model format'
4192 mod = LayeredModel.from_scanlines(reader)
4193 if crust2_profile is not None:
4194 return mod.replaced_crust(crust2_profile)
4196 return mod
4198 else:
4199 assert crust2_profile is not None
4200 profile = anything_to_crust2_profile(crust2_profile)
4201 return LayeredModel.from_scanlines(
4202 from_crust2x2_profile(profile))
4205def castagna_vs_to_vp(vs):
4206 '''
4207 Calculate vp from vs using Castagna's relation.
4209 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4210 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4211 1985]
4213 vp = 1.16 * vs + 1360 [m/s]
4215 :param vs: S-wave velocity [m/s]
4216 :returns: P-wave velocity [m/s]
4217 '''
4219 return vs*1.16 + 1360.0
4222def castagna_vp_to_vs(vp):
4223 '''
4224 Calculate vp from vs using Castagna's relation.
4226 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4227 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4228 1985]
4230 vp = 1.16 * vs + 1360 [m/s]
4232 :param vp: P-wave velocity [m/s]
4233 :returns: S-wave velocity [m/s]
4234 '''
4236 return (vp - 1360.0) / 1.16
4239def evenize(x, y, minsize=10):
4240 if x.size < minsize:
4241 return x
4242 ry = (y.max()-y.min())
4243 if ry == 0:
4244 return x
4245 dx = (x[1:] - x[:-1])/(x.max()-x.min())
4246 dy = (y[1:] + y[:-1])/ry
4248 s = num.zeros(x.size)
4249 s[1:] = num.cumsum(num.sqrt(dy**2 + dx**2))
4250 s2 = num.linspace(0, s[-1], x.size)
4251 x2 = num.interp(s2, s, x)
4252 x2[0] = x[0]
4253 x2[-1] = x[-1]
4254 return x2
4257def filled(v, *args, **kwargs):
4258 '''
4259 Create NumPy array filled with given value.
4261 This works like :py:func:`numpy.ones` but initializes the array with ``v``
4262 instead of ones.
4263 '''
4264 x = num.empty(*args, **kwargs)
4265 x.fill(v)
4266 return x
4269def next_or_none(i):
4270 try:
4271 return next(i)
4272 except StopIteration:
4273 return None
4276def reci_or_none(x):
4277 try:
4278 return 1./x
4279 except ZeroDivisionError:
4280 return None
4283def mult_or_none(a, b):
4284 if a is None or b is None:
4285 return None
4286 return a*b
4289def min_not_none(a, b):
4290 if a is None:
4291 return b
4292 if b is None:
4293 return a
4294 return min(a, b)
4297def xytups(xx, yy):
4298 d = []
4299 for x, y in zip(xx, yy):
4300 if num.isfinite(y):
4301 d.append((x, y))
4302 return d
4305def interp(x, xp, fp, monoton):
4306 if monoton == 1:
4307 return xytups(
4308 x, num.interp(x, xp, fp, left=num.nan, right=num.nan))
4309 elif monoton == -1:
4310 return xytups(
4311 x, num.interp(x, xp[::-1], fp[::-1], left=num.nan, right=num.nan))
4312 else:
4313 fs = []
4314 for xv in x:
4315 indices = num.where(num.logical_or(
4316 num.logical_and(xp[:-1] >= xv, xv > xp[1:]),
4317 num.logical_and(xp[:-1] <= xv, xv < xp[1:])))[0]
4319 for i in indices:
4320 xr = (xv - xp[i])/(xp[i+1]-xp[i])
4321 fv = xr*fp[i] + (1.-xr)*fp[i+1]
4322 fs.append((xv, fv))
4324 return fs
4327def float_or_none(x):
4328 if x is not None:
4329 return float(x)
4332def parstore_float(thelocals, obj, *args):
4333 for k, v in thelocals.items():
4334 if k != 'self' and (not args or k in args):
4335 setattr(obj, k, float_or_none(v))