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 t_and_attributes(self, attributes):
2934 d = {
2935 'takeoff_angle': lambda: self.takeoff_angle(),
2936 'incidence_angle': lambda: self.incidence_angle(),
2937 't': lambda: self.t,
2938 'p': lambda: self.p}
2940 return tuple(d[k]() for k in ['t'] + attributes)
2942 def takeoff_angle(self):
2943 '''
2944 Get takeoff angle of ray.
2946 The angle is returned in [degrees].
2947 '''
2949 return self.path.first_straight().angle_in(self.p, self.endgaps)
2951 def incidence_angle(self):
2952 '''
2953 Get incidence angle of ray.
2955 The angle is returned in [degrees].
2956 '''
2958 return self.path.last_straight().angle_out(self.p, self.endgaps)
2960 def efficiency(self):
2961 '''
2962 Get conversion/reflection efficiency of the ray.
2964 A value between 0 and 1 is returned, reflecting the relative amount of
2965 energy which is transmitted along the ray and not lost by reflections
2966 or conversions.
2967 '''
2969 return self.path.efficiency(self.p)
2971 def spreading(self):
2972 '''
2973 Get geometrical spreading factor.
2974 '''
2976 return self.path.spreading(self.p, self.endgaps)
2978 def surface_sphere(self):
2979 x1, y1 = 0., earthradius - self.endgaps[0]
2980 r2 = earthradius - self.endgaps[1]
2981 x2, y2 = r2*math.sin(self.x*d2r), r2*math.cos(self.x*d2r)
2982 return ((x2-x1)**2 + (y2-y1)**2)*4.0*math.pi
2984 def zxt_path_subdivided(self, points_per_straight=20):
2985 '''
2986 Get geometrical representation of ray path.
2988 Three arrays (depth, distance, time) with points on the ray's path of
2989 propagation are returned. The number of points which are used in each
2990 ray segment (passage through one layer) may be controlled by the
2991 ``points_per_straight`` parameter.
2992 '''
2993 return self.path.zxt_path_subdivided(
2994 num.atleast_1d(self.p), self.endgaps,
2995 points_per_straight=points_per_straight,
2996 x_for_headwave=num.atleast_1d(self.x))
2998 def __str__(self, as_degrees=False):
2999 if as_degrees:
3000 sd = '%6.3g deg' % self.x
3001 else:
3002 sd = '%7.5g km' % (self.x*(d2r*earthradius/km))
3004 return '%7.5g s/deg %s %6.4g s %5.1f %5.1f %3.0f%% %3.0f%% %s' % (
3005 self.p/r2d,
3006 sd,
3007 self.t,
3008 self.takeoff_angle(),
3009 self.incidence_angle(),
3010 100*self.efficiency(),
3011 100*self.spreading()*self.surface_sphere(),
3012 self.path.__str__(p=self.p))
3015def anything_to_crust2_profile(crust2_profile):
3016 from pyrocko.dataset import crust2x2
3017 if isinstance(crust2_profile, tuple):
3018 lat, lon = [float(x) for x in crust2_profile]
3019 return crust2x2.get_profile(lat, lon)
3020 elif isinstance(crust2_profile, str):
3021 return crust2x2.get_profile(crust2_profile)
3022 elif isinstance(crust2_profile, crust2x2.Crust2Profile):
3023 return crust2_profile
3024 else:
3025 assert False, 'crust2_profile must be (lat, lon) a profile ' \
3026 'key or a crust2x2 Profile object)'
3029class DiscontinuityNotFound(CakeError):
3030 def __init__(self, depth_or_name):
3031 CakeError.__init__(self)
3032 self.depth_or_name = depth_or_name
3034 def __str__(self):
3035 return 'Cannot find discontinuity from given depth or name: %s' % \
3036 self.depth_or_name
3039class LayeredModelError(CakeError):
3040 pass
3043class LayeredModel(object):
3044 '''
3045 Representation of a layer cake model.
3047 There are several ways to initialize an instance of this class.
3049 1. Use the module function :py:func:`load_model` to read a model from a
3050 file.
3051 2. Create an empty model with the default constructor and append layers and
3052 discontinuities with the :py:meth:`append` method (from top to bottom).
3053 3. Use the constructor :py:meth:`LayeredModel.from_scanlines`, to
3054 automatically create the :py:class:`Layer` and :py:class:`Discontinuity`
3055 objects from a given velocity profile.
3057 An earth model is represented by as stack of :py:class:`Layer` and
3058 :py:class:`Discontinuity` objects. The method :py:meth:`arrivals` returns
3059 :py:class:`Ray` objects which may be e.g. queried for arrival times of
3060 specific phases. Each ray is associated with a :py:class:`RayPath` object.
3061 Ray objects share common ray paths if they have the same
3062 conversion/reflection/propagation history. Creating the ray path objects is
3063 relatively expensive (this is done in :py:meth:`gather_paths`), but they
3064 are cached for reuse in successive invocations.
3065 '''
3067 def __init__(self):
3068 self._surface_material = None
3069 self._elements = []
3070 self.nlayers = 0
3071 self._np = 10000
3072 self._pdepth = 5
3073 self._pathcache = {}
3075 def copy_with_elevation(self, elevation):
3076 '''
3077 Get copy of the model with surface layer stretched to given elevation.
3079 :param elevation: new surface elevation in [m]
3081 Elevation is positiv upward, contrary to the layered models downward
3082 `z` axis.
3083 '''
3085 c = copy.deepcopy(self)
3086 c._pathcache = {}
3087 surface = c._elements[0]
3088 toplayer = c._elements[1]
3090 assert toplayer.zbot > -elevation
3092 surface.z = -elevation
3093 c._elements[1] = toplayer.copy(ztop=-elevation)
3094 c._elements[1].ilayer = 0
3095 return c
3097 def zeq(self, z1, z2):
3098 return abs(z1-z2) < ZEPS
3100 def append(self, element):
3101 '''
3102 Add a layer or discontinuity at bottom of model.
3104 :param element: object of subclass of :py:class:`Layer` or
3105 :py:class:`Discontinuity`.
3106 '''
3108 if isinstance(element, Layer):
3109 if element.zbot >= earthradius:
3110 element.zbot = earthradius - 1.
3112 if element.ztop >= earthradius:
3113 raise CakeError('Layer deeper than earthradius')
3115 element.ilayer = self.nlayers
3116 self.nlayers += 1
3118 self._elements.append(element)
3120 def elements(self, direction=DOWN):
3121 '''
3122 Iterate over all elements of the model.
3124 :param direction: direction of traversal :py:const:`DOWN` or
3125 :py:const:`UP`.
3127 Objects derived from the :py:class:`Discontinuity` and
3128 :py:class:`Layer` classes are yielded.
3129 '''
3131 if direction == DOWN:
3132 return iter(self._elements)
3133 else:
3134 return reversed(self._elements)
3136 def layers(self, direction=DOWN):
3137 '''
3138 Iterate over all layers of model.
3140 :param direction: direction of traversal :py:const:`DOWN` or
3141 :py:const:`UP`.
3143 Objects derived from the :py:class:`Layer` class are yielded.
3144 '''
3146 if direction == DOWN:
3147 return (el for el in self._elements if isinstance(el, Layer))
3148 else:
3149 return (
3150 el for el in reversed(self._elements) if isinstance(el, Layer))
3152 def layer(self, z, direction=DOWN):
3153 '''
3154 Get layer for given depth.
3156 :param z: depth [m]
3157 :param direction: direction of traversal :py:const:`DOWN` or
3158 :py:const:`UP`.
3160 Returns first layer which touches depth ``z`` (tolerant at boundaries).
3161 '''
3163 for layer in self.layers(direction):
3164 if layer.contains(z):
3165 return layer
3166 else:
3167 raise CakeError('Failed extracting layer at depth z=%s' % z)
3169 def walker(self):
3170 return Walker(self._elements)
3172 def material(self, z, direction=DOWN):
3173 '''
3174 Get material at given depth.
3176 :param z: depth [m]
3177 :param direction: direction of traversal :py:const:`DOWN` or
3178 :py:const:`UP`
3179 :returns: object of type :py:class:`Material`
3181 If given depth ``z`` happens to be at an interface, the material of the
3182 first layer with respect to the the traversal ordering is returned.
3183 '''
3185 lyr = self.layer(z, direction)
3186 return lyr.material(z)
3188 def discontinuities(self):
3189 '''
3190 Iterate over all discontinuities of the model.'''
3192 return (el for el in self._elements if isinstance(el, Discontinuity))
3194 def discontinuity(self, name_or_z):
3195 '''
3196 Get discontinuity by name or depth.
3198 :param name_or_z: name of discontinuity or depth [m] as float value
3199 '''
3201 if isinstance(name_or_z, float):
3202 candi = sorted(
3203 self.discontinuities(), key=lambda i: abs(i.z-name_or_z))
3204 else:
3205 candi = [i for i in self.discontinuities() if i.name == name_or_z]
3207 if not candi:
3208 raise DiscontinuityNotFound(name_or_z)
3210 return candi[0]
3212 def adapt_phase(self, phase):
3213 '''
3214 Adapt a phase definition for use with this model.
3216 This returns a copy of the phase definition, where named
3217 discontinuities are replaced with the actual depth of these, as defined
3218 in the model.
3219 '''
3221 phase = phase.copy()
3222 for knee in phase.knees():
3223 if knee.depth != 'surface':
3224 knee.depth = self.discontinuity(knee.depth).z
3225 for leg in phase.legs():
3226 if leg.depthmax is not None and isinstance(leg.depthmax, str):
3227 leg.depthmax = self.discontinuity(leg.depthmax).z
3229 return phase
3231 def path(self, p, phase, layer_start, layer_stop):
3232 '''
3233 Get ray path for given combination of ray parameter, phase definition,
3234 source and receiver layers.
3236 :param p: ray parameter (spherical) [s/rad]
3237 :param phase: phase definition (:py:class:`PhaseDef` object)
3238 :param layer_start: layer with source
3239 :param layer_stop: layer with receiver
3240 :returns: :py:class:`RayPath` object
3242 If it is not possible to find a solution, an exception of type
3243 :py:exc:`NotPhaseConform`, :py:exc:`MinDepthReached`,
3244 :py:exc:`MaxDepthReached`, :py:exc:`CannotPropagate`,
3245 :py:exc:`BottomReached` or :py:exc:`SurfaceReached` is raised.
3246 '''
3248 phase = self.adapt_phase(phase)
3249 knees = phase.knees()
3250 legs = phase.legs()
3251 next_knee = next_or_none(knees)
3252 leg = next_or_none(legs)
3253 assert leg is not None
3255 direction = leg.departure
3256 direction_stop = phase.direction_stop()
3257 mode = leg.mode
3258 mode_stop = phase.last_leg().mode
3260 walker = self.walker()
3261 walker.goto_layer(layer_start)
3262 current = walker.current()
3264 ttop, tbot = current.tests(p, mode)
3265 if not ttop and not tbot:
3266 raise CannotPropagate(direction, current.ilayer)
3268 if (direction == DOWN and not ttop) or (direction == UP and not tbot):
3269 direction = -direction
3271 path = RayPath(phase)
3272 trapdetect = set()
3273 while True:
3274 at_layer = isinstance(current, Layer)
3275 at_discontinuity = isinstance(current, Discontinuity)
3277 # detect trapped wave
3278 k = (id(next_knee), id(current), direction, mode)
3279 if k in trapdetect:
3280 raise Trapped()
3282 trapdetect.add(k)
3284 if at_discontinuity:
3285 oldmode, olddirection = mode, direction
3286 headwave = False
3287 if next_knee is not None and next_knee.matches(
3288 current, mode, direction):
3290 headwave = next_knee.headwave
3291 direction = next_knee.out_direction()
3292 mode = next_knee.out_mode
3293 next_knee = next_or_none(knees)
3294 leg = next(legs)
3296 else: # implicit reflection/transmission
3297 direction = current.propagate(p, mode, direction)
3299 if headwave:
3300 path.set_is_headwave(True)
3302 path.append(Kink(
3303 olddirection, olddirection, oldmode, oldmode, current))
3305 path.append(HeadwaveStraight(
3306 olddirection, direction, oldmode, current))
3308 path.append(Kink(
3309 olddirection, direction, oldmode, mode, current))
3311 else:
3312 path.append(Kink(
3313 olddirection, direction, oldmode, mode, current))
3315 if at_layer:
3316 direction_in = direction
3317 direction = current.propagate(p, mode, direction_in)
3319 zturn = None
3320 if direction_in != direction:
3321 zturn = current.zturn(p, mode)
3323 zmin, zmax = leg.depthmin, leg.depthmax
3324 if zmin is not None or zmax is not None:
3325 if direction_in != direction:
3326 if zmin is not None and zturn <= zmin:
3327 raise MinDepthReached()
3328 if zmax is not None and zturn >= zmax:
3329 raise MaxDepthReached()
3330 else:
3331 if zmin is not None and current.ztop <= zmin:
3332 raise MinDepthReached()
3333 if zmax is not None and current.zbot >= zmax:
3334 raise MaxDepthReached()
3336 path.append(Straight(direction_in, direction, mode, current))
3338 if next_knee is None and mode == mode_stop and \
3339 current is layer_stop:
3341 if zturn is None:
3342 if direction == direction_stop:
3343 break
3344 else:
3345 break
3347 walker.go(direction)
3348 current = walker.current()
3350 return path
3352 def gather_paths(self, phases=PhaseDef('P'), zstart=0.0, zstop=0.0):
3353 '''
3354 Get all possible ray paths for given source and receiver depths for one
3355 or more phase definitions.
3357 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3358 Comma-separated strings and lists of such strings are also accepted
3359 and are converted to :py:class:`PhaseDef` objects for convenience.
3360 :param zstart: source depth [m]
3361 :param zstop: receiver depth [m]
3362 :returns: a list of :py:class:`RayPath` objects
3364 Results of this method are cached internally. Cached results are
3365 returned, when a given combination of source layer, receiver layer and
3366 phase definition has been used before.
3367 '''
3369 eps = 1e-7 # num.finfo(float).eps * 1000.
3371 phases = to_phase_defs(phases)
3373 paths = []
3374 for phase in phases:
3376 layer_start = self.layer(zstart, -phase.direction_start())
3377 layer_stop = self.layer(zstop, phase.direction_stop())
3379 pathcachekey = (phase.definition(), layer_start, layer_stop)
3381 if pathcachekey in self._pathcache:
3382 phase_paths = self._pathcache[pathcachekey]
3383 else:
3384 hwknee = phase.headwave_knee()
3385 if hwknee:
3386 name_or_z = hwknee.depth
3387 interface = self.discontinuity(name_or_z)
3388 mode = hwknee.in_mode
3389 in_direction = hwknee.direction
3391 pabove, pbelow = interface.critical_ps(mode)
3393 p = min_not_none(pabove, pbelow)
3395 # diffracted wave:
3396 if in_direction == DOWN and (
3397 pbelow is None or pbelow >= pabove):
3399 p *= (1.0 - eps)
3401 path = self.path(p, phase, layer_start, layer_stop)
3402 path.set_prange(p, p, 1.)
3404 phase_paths = [path]
3406 else:
3407 try:
3408 pmax_start = max([
3409 radius(z)/layer_start.v(phase.first_leg().mode, z)
3410 for z in (layer_start.ztop, layer_start.zbot)])
3412 pmax_stop = max([
3413 radius(z)/layer_stop.v(phase.last_leg().mode, z)
3414 for z in (layer_stop.ztop, layer_stop.zbot)])
3416 pmax = min(pmax_start, pmax_stop)
3418 pedges = [0.]
3419 for layer in self.layers():
3420 for z in (layer.ztop, layer.zbot):
3421 for mode in (P, S):
3422 for eps2 in [eps]:
3423 v = layer.v(mode, z)
3424 if v != 0.0:
3425 p = radius(z)/v
3426 if p <= pmax:
3427 pedges.append(p*(1.0-eps2))
3428 pedges.append(p)
3429 pedges.append(p*(1.0+eps2))
3431 pedges = num.unique(sorted(pedges))
3433 phase_paths = {}
3434 cached = {}
3435 counter = [0]
3437 def p_to_path(p):
3438 if p in cached:
3439 return cached[p]
3441 try:
3442 counter[0] += 1
3443 path = self.path(
3444 p, phase, layer_start, layer_stop)
3446 if path not in phase_paths:
3447 phase_paths[path] = []
3449 phase_paths[path].append(p)
3451 except PathFailed:
3452 path = None
3454 cached[p] = path
3455 return path
3457 def recurse(pmin, pmax, i=0):
3458 if i > self._pdepth:
3459 return
3460 path1 = p_to_path(pmin)
3461 path2 = p_to_path(pmax)
3462 if path1 is None and path2 is None and i > 0:
3463 return
3464 if path1 is None or path2 is None or \
3465 hash(path1) != hash(path2):
3467 recurse(pmin, (pmin+pmax)/2., i+1)
3468 recurse((pmin+pmax)/2., pmax, i+1)
3470 for (pl, ph) in zip(pedges[:-1], pedges[1:]):
3471 recurse(pl, ph)
3473 for path, ps in phase_paths.items():
3474 path.set_prange(
3475 min(ps), max(ps), pmax/(self._np-1))
3477 phase_paths = list(phase_paths.keys())
3479 except ZeroDivisionError:
3480 phase_paths = []
3482 self._pathcache[pathcachekey] = phase_paths
3484 paths.extend(phase_paths)
3486 paths.sort(key=lambda x: x.pmin())
3487 return paths
3489 def arrivals(
3490 self,
3491 distances=[],
3492 phases=PhaseDef('P'),
3493 zstart=0.0,
3494 zstop=0.0,
3495 refine=True):
3497 '''
3498 Compute rays and traveltimes for given distances.
3500 :param distances: list or array of distances [deg]
3501 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3502 Comma-separated strings and lists of such strings are also accepted
3503 and are converted to :py:class:`PhaseDef` objects for convenience.
3504 :param zstart: source depth [m]
3505 :param zstop: receiver depth [m]
3506 :param refine: bool flag, whether to use bisectioning to improve
3507 (p, x, t) estimated from interpolation
3508 :returns: a list of :py:class:`Ray` objects, sorted by
3509 (distance, arrival time)
3510 '''
3512 distances = num.asarray(distances, dtype=float)
3514 arrivals = []
3515 for path in self.gather_paths(phases, zstart=zstart, zstop=zstop):
3517 endgaps = path.endgaps(zstart, zstop)
3518 for x, p, t, draft_pxt in path.interpolate_x2pt_linear(
3519 distances, endgaps):
3521 arrivals.append(Ray(path, p, x, t, endgaps, draft_pxt))
3523 if refine:
3524 refined = []
3525 for ray in arrivals:
3527 if ray.path._is_headwave:
3528 refined.append(ray)
3530 try:
3531 ray.refine()
3532 refined.append(ray)
3534 except RefineFailed:
3535 pass
3537 arrivals = refined
3539 arrivals.sort(key=lambda x: (x.x, x.t))
3540 return arrivals
3542 @classmethod
3543 def from_scanlines(cls, producer):
3544 '''
3545 Create layer cake model from sequence of materials at depths.
3547 :param producer: iterable yielding (depth, material, name) tuples
3549 Creates a new :py:class:`LayeredModel` object and uses its
3550 :py:meth:`append` method to add layers and discontinuities as needed.
3551 '''
3553 self = cls()
3554 for z, material, name in producer:
3556 if not self._elements:
3557 self.append(Surface(z, material))
3558 else:
3559 element = self._elements[-1]
3560 if self.zeq(element.zbot, z):
3561 assert isinstance(element, Layer)
3562 self.append(
3563 Interface(z, element.mbot, material, name=name))
3565 else:
3566 if isinstance(element, Discontinuity):
3567 ztop = element.z
3568 mtop = element.mbelow
3569 elif isinstance(element, Layer):
3570 ztop = element.zbot
3571 mtop = element.mbot
3573 if mtop == material:
3574 layer = HomogeneousLayer(
3575 ztop, z, material, name=name)
3576 else:
3577 layer = GradientLayer(
3578 ztop, z, mtop, material, name=name)
3580 self.append(layer)
3582 return self
3584 def to_scanlines(self, get_burgers=False):
3585 def fmt(z, m):
3586 if not m._has_default_burgers() or get_burgers:
3587 return (z, m.vp, m.vs, m.rho, m.qp, m.qs,
3588 m.burger_eta1, m.burger_eta2, m.burger_valpha)
3589 return (z, m.vp, m.vs, m.rho, m.qp, m.qs)
3591 last = None
3592 lines = []
3593 for element in self.elements():
3594 if isinstance(element, Layer):
3595 if not isinstance(last, Layer):
3596 lines.append(fmt(element.ztop, element.mtop))
3598 lines.append(fmt(element.zbot, element.mbot))
3600 last = element
3602 if not isinstance(last, Layer):
3603 lines.append(fmt(last.z, last.mbelow))
3605 return lines
3607 def iter_material_parameter(self, get):
3608 assert get in ('vp', 'vs', 'rho', 'qp', 'qs', 'z')
3609 if get == 'z':
3610 for layer in self.layers():
3611 yield layer.ztop
3612 yield layer.zbot
3613 else:
3614 getter = operator.attrgetter(get)
3615 for layer in self.layers():
3616 yield getter(layer.mtop)
3617 yield getter(layer.mbot)
3619 def profile(self, get):
3620 '''
3621 Get parameter profile along depth of the earthmodel.
3623 :param get: property to be queried (
3624 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3625 :type get: string
3626 '''
3628 return num.array(list(self.iter_material_parameter(get)))
3630 def min(self, get='vp'):
3631 '''
3632 Find minimum value of a material property or depth.
3634 :param get: property to be queried (
3635 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3636 '''
3638 return min(self.iter_material_parameter(get))
3640 def max(self, get='vp'):
3641 '''
3642 Find maximum value of a material property or depth.
3644 :param get: property to be queried (
3645 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, ``'qs'``, or ``'z'``)
3646 '''
3648 return max(self.iter_material_parameter(get))
3650 def simplify_layers(self, layers, max_rel_error=0.001):
3651 if len(layers) <= 1:
3652 return layers
3654 ztop = layers[0].ztop
3655 zbot = layers[-1].zbot
3656 zorigs = [layer.ztop for layer in layers]
3657 zorigs.append(zbot)
3658 zs = num.linspace(ztop, zbot, 100)
3659 data = []
3660 for z in zs:
3661 if z == ztop:
3662 direction = UP
3663 else:
3664 direction = DOWN
3666 mat = self.material(z, direction)
3667 data.append(mat.astuple())
3669 data = num.array(data, dtype=float)
3670 data_means = num.mean(data, axis=0)
3671 nmax = len(layers) // 2
3672 accept = False
3674 zcut_best = []
3675 for n in range(1, nmax+1):
3676 ncutintervals = 20
3677 zdelta = (zbot-ztop)/ncutintervals
3678 if n == 2:
3679 zcuts = [
3680 [ztop, ztop + i*zdelta, zbot]
3681 for i in range(1, ncutintervals)]
3682 elif n == 3:
3683 zcuts = []
3684 for j in range(1, ncutintervals):
3685 for i in range(j+1, ncutintervals):
3686 zcuts.append(
3687 [ztop, ztop + j*zdelta, ztop + i*zdelta, zbot])
3688 else:
3689 zcuts = []
3690 zcuts.append(num.linspace(ztop, zbot, n+1))
3691 if zcut_best:
3692 zcuts.append(sorted(num.linspace(
3693 ztop, zbot, n).tolist() + zcut_best[1]))
3694 zcuts.append(sorted(num.linspace(
3695 ztop, zbot, n-1).tolist() + zcut_best[2]))
3697 best = None
3698 for icut, zcut in enumerate(zcuts):
3699 rel_par_errors = num.zeros(5)
3700 mpar_nodes = num.zeros((n+1, 5))
3702 for ipar in range(5):
3703 znodes, vnodes, error_rms = util.polylinefit(
3704 zs, data[:, ipar], zcut)
3706 mpar_nodes[:, ipar] = vnodes
3707 if data_means[ipar] == 0.0:
3708 rel_par_errors[ipar] = -1
3709 else:
3710 rel_par_errors[ipar] = error_rms/data_means[ipar]
3712 rel_error = rel_par_errors.max()
3713 if best is None or rel_error < best[0]:
3714 best = (rel_error, zcut, mpar_nodes)
3716 rel_error, zcut, mpar_nodes = best
3718 zcut_best.append(list(zcut))
3719 zcut_best[-1].pop(0)
3720 zcut_best[-1].pop()
3722 if rel_error <= max_rel_error:
3723 accept = True
3724 break
3726 if not accept:
3727 return layers
3729 rel_error, zcut, mpar_nodes = best
3731 material_nodes = []
3732 for i in range(n+1):
3733 material_nodes.append(Material(*mpar_nodes[i, :]))
3735 out_layers = []
3736 for i in range(n):
3737 mtop = material_nodes[i]
3738 mbot = material_nodes[i+1]
3739 ztop = zcut[i]
3740 zbot = zcut[i+1]
3741 if mtop == mbot:
3742 lyr = HomogeneousLayer(ztop, zbot, mtop)
3743 else:
3744 lyr = GradientLayer(ztop, zbot, mtop, mbot)
3746 out_layers.append(lyr)
3747 return out_layers
3749 def simplify(self, max_rel_error=0.001):
3750 '''
3751 Get representation of model with lower resolution.
3753 Returns an approximation of the model. All discontinuities are kept,
3754 but layer stacks with continuous model parameters are represented, if
3755 possible, by a lower number of layers. Piecewise linear functions are
3756 fitted against the original model parameter's piecewise linear
3757 functions. Successively larger numbers of layers are tried, until the
3758 difference to the original model is below ``max_rel_error``. The
3759 difference is measured as the RMS error of the fit normalized by the
3760 mean of the input (i.e. the fitted curves should deviate, on average,
3761 less than 0.1% from the input curves if ``max_rel_error`` = 0.001).
3762 '''
3764 mod_simple = LayeredModel()
3766 glayers = []
3767 for element in self.elements():
3769 if isinstance(element, Discontinuity):
3770 for layer in self.simplify_layers(
3771 glayers, max_rel_error=max_rel_error):
3773 mod_simple.append(layer)
3775 glayers = []
3776 mod_simple.append(element)
3777 else:
3778 glayers.append(element)
3780 for layer in self.simplify_layers(
3781 glayers, max_rel_error=max_rel_error):
3782 mod_simple.append(layer)
3784 return mod_simple
3786 def extract(self, depth_min=None, depth_max=None):
3787 '''
3788 Extract :py:class:`LayeredModel` from :py:class:`LayeredModel`.
3790 :param depth_min: depth of upper cut or name of :py:class:`Interface`
3791 :param depth_max: depth of lower cut or name of :py:class:`Interface`
3793 Interpolates a :py:class:`GradientLayer` at ``depth_min`` and/or
3794 ``depth_max``.
3795 '''
3797 if isinstance(depth_min, str):
3798 depth_min = self.discontinuity(depth_min).z
3800 if isinstance(depth_max, str):
3801 depth_max = self.discontinuity(depth_max).z
3803 mod_extracted = LayeredModel()
3805 for element in self.elements():
3806 element = element.copy()
3807 do_append = False
3808 if (depth_min is None or depth_min <= element.ztop) \
3809 and (depth_max is None or depth_max >= element.zbot):
3810 mod_extracted.append(element)
3811 continue
3813 if depth_min is not None:
3814 if element.ztop < depth_min and depth_min < element.zbot:
3815 _, element = element.split(depth_min)
3816 do_append = True
3818 if depth_max is not None:
3819 if element.zbot > depth_max and depth_max > element.ztop:
3820 element, _ = element.split(depth_max)
3821 do_append = True
3823 if do_append:
3824 mod_extracted.append(element)
3826 return mod_extracted
3828 def replaced_crust(self, crust2_profile=None, crustmod=None):
3829 if crust2_profile is not None:
3830 profile = anything_to_crust2_profile(crust2_profile)
3831 crustmod = LayeredModel.from_scanlines(
3832 from_crust2x2_profile(profile))
3834 newmod = LayeredModel()
3835 for element in crustmod.extract(depth_max='moho').elements():
3836 if element.name != 'moho':
3837 newmod.append(element)
3838 else:
3839 moho1 = element
3841 mod = self.extract(depth_min='moho')
3842 first = True
3843 for element in mod.elements():
3844 if element.name == 'moho':
3845 if element.z <= moho1.z:
3846 mbelow = mod.material(moho1.z, direction=UP)
3847 else:
3848 mbelow = element.mbelow
3850 moho = Interface(moho1.z, moho1.mabove, mbelow, name='moho')
3851 newmod.append(moho)
3852 else:
3853 if first:
3854 if isinstance(element, Layer) and element.zbot > moho.z:
3855 newmod.append(GradientLayer(
3856 moho.z,
3857 element.zbot,
3858 moho.mbelow,
3859 element.mbot,
3860 name=element.name))
3862 first = False
3863 else:
3864 newmod.append(element)
3865 return newmod
3867 def perturb(self, rstate=None, keep_vp_vs=False, **kwargs):
3868 '''
3869 Create a perturbed variant of the earth model.
3871 Randomly change the thickness and material parameters of the earth
3872 model from a uniform distribution.
3874 :param kwargs: Maximum change fraction (e.g. 0.1) of the parameters.
3875 Name the parameter, prefixed by ``p``. Supported parameters are
3876 ``ph, pvp, pvs, prho, pqs, pqp``.
3877 :type kwargs: dict
3878 :param rstate: Random state to draw from, defaults to ``None``
3879 :type rstate: :class:`numpy.random.RandomState`, optional
3880 :param keep_vp_vs: Keep the Vp/Vs ratio, defaults to False
3881 :type keep_vp_vs: bool, optional
3883 :returns: A new, perturbed earth model
3884 :rtype: :class:`~pyrocko.cake.LayeredModel`
3886 .. code-block :: python
3888 perturbed_model = model.perturb(ph=.1, pvp=.05, prho=.1)
3889 '''
3890 _pargs = set(['ph', 'pvp', 'pvs', 'prho', 'pqs', 'pqp'])
3891 earthmod = copy.deepcopy(self)
3893 if rstate is None:
3894 rstate = num.random.RandomState()
3896 layers = earthmod.layers()
3897 discont = earthmod.discontinuities()
3898 prev_layer = None
3900 def get_change_ratios():
3901 values = dict.fromkeys([p[1:] for p in _pargs], 0.)
3903 for param, pval in kwargs.items():
3904 if param not in _pargs:
3905 continue
3906 values[param[1:]] = float(rstate.uniform(-pval, pval, size=1))
3907 return values
3909 # skip Surface
3910 while True:
3911 disc = next(discont)
3912 if isinstance(disc, Surface):
3913 break
3915 while True:
3916 try:
3917 layer = next(layers)
3918 m = layer.material(None)
3919 h = layer.zbot - layer.ztop
3920 except StopIteration:
3921 break
3923 if not isinstance(layer, HomogeneousLayer):
3924 raise NotImplementedError(
3925 'Can only perturbate homogeneous layers!')
3927 changes = get_change_ratios()
3929 # Changing thickness
3930 dh = h * changes['h']
3931 changes['h'] = dh
3933 layer.resize(depth_max=layer.zbot + dh,
3934 depth_min=prev_layer.zbot if prev_layer else None)
3936 try:
3937 disc = next(discont)
3938 disc.change_depth(disc.z + dh)
3939 except StopIteration:
3940 pass
3942 # Setting material parameters
3943 for param, change_ratio in changes.items():
3944 if param == 'h':
3945 continue
3947 value = m.__getattribute__(param)
3948 changes[param] = value * change_ratio
3950 if keep_vp_vs and changes['vp'] != 0.:
3951 changes['vs'] = (m.vp + changes['vp']) / m.vp_vs_ratio() - m.vs
3953 for param, change in changes.items():
3954 if param == 'h':
3955 continue
3956 value = m.__getattribute__(param)
3957 m.__setattr__(param, value + change)
3959 logger.info(
3960 'perturbating earthmodel: {}'.format(
3961 ' '.join(['{param}: {change:{len}.2f}'.format(
3962 param=p, change=c, len=8)
3963 for p, c in changes.items()])))
3965 prev_layer = layer
3967 return earthmod
3969 def require_homogeneous(self):
3970 elements = list(self.elements())
3972 if len(elements) != 2:
3973 raise LayeredModelError(
3974 'Homogeneous model required but found more than one layer in '
3975 'earthmodel.')
3977 if not isinstance(elements[1], HomogeneousLayer):
3978 raise LayeredModelError(
3979 'Homogeneous model required but element #1 is not of type '
3980 'HomogeneousLayer.')
3982 return elements[1].m
3984 def __str__(self):
3985 return '\n'.join(str(element) for element in self._elements)
3988def read_hyposat_model(fn):
3989 '''
3990 Reader for HYPOSAT earth model files.
3992 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
3994 Interface names are translated as follows: ``'MOHO'`` -> ``'moho'``,
3995 ``'CONR'`` -> ``'conrad'``
3996 '''
3998 with open(fn, 'r') as f:
3999 translate = {'MOHO': 'moho', 'CONR': 'conrad'}
4000 lname = None
4001 for iline, line in enumerate(f):
4002 if iline == 0:
4003 continue
4005 z, vp, vs, name = util.unpack_fixed('f10, f10, f10, a4', line)
4006 if not name:
4007 name = None
4008 material = Material(vp*1000., vs*1000.)
4010 tname = translate.get(lname, lname)
4011 yield z*1000., material, tname
4013 lname = name
4016def read_nd_model(fn):
4017 '''
4018 Reader for TauP style '.nd' (named discontinuity) files.
4020 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
4022 Interface names are translated as follows: ``'mantle'`` -> ``'moho'``,
4023 ``'outer-core'`` -> ``'cmb'``, ``'inner-core'`` -> ``'icb'``.
4025 The format has been modified to include Burgers materials parameters in
4026 columns 7 (burger_eta1), 8 (burger_eta2) and 9. eta(3).
4027 '''
4028 with open(fn, 'r') as f:
4029 for x in read_nd_model_fh(f):
4030 yield x
4033def read_nd_model_str(s):
4034 f = StringIO(s)
4035 for x in read_nd_model_fh(f):
4036 yield x
4037 f.close()
4040def read_nd_model_fh(f):
4041 translate = {'mantle': 'moho', 'outer-core': 'cmb', 'inner-core': 'icb'}
4042 name = None
4043 for line in f:
4044 toks = line.split()
4045 if len(toks) == 9 or len(toks) == 6 or len(toks) == 4:
4046 z, vp, vs, rho = [float(x) for x in toks[:4]]
4047 qp, qs = None, None
4048 burgers = None
4049 if len(toks) == 6 or len(toks) == 9:
4050 qp, qs = [float(x) for x in toks[4:6]]
4051 if len(toks) == 9:
4052 burgers = \
4053 [float(x) for x in toks[6:]]
4055 material = Material(
4056 vp*1000., vs*1000., rho*1000., qp, qs,
4057 burgers=burgers)
4059 yield z*1000., material, name
4060 name = None
4061 elif len(toks) == 1:
4062 name = translate.get(toks[0], toks[0])
4064 f.close()
4067def from_crust2x2_profile(profile, depthmantle=50000):
4068 from pyrocko.dataset import crust2x2
4070 default_qp_qs = {
4071 'soft sed.': (50., 50.),
4072 'hard sed.': (200., 200.),
4073 'upper crust': (600., 400.),
4074 }
4076 z = 0.
4077 for i in range(8):
4078 dz, vp, vs, rho = profile.get_layer(i)
4079 name = crust2x2.Crust2Profile.layer_names[i]
4080 if name in default_qp_qs:
4081 qp, qs = default_qp_qs[name]
4082 else:
4083 qp, qs = None, None
4085 material = Material(vp, vs, rho, qp, qs)
4086 iname = None
4087 if i == 7:
4088 iname = 'moho'
4089 if dz != 0.0:
4090 yield z, material, iname
4091 if i != 7:
4092 yield z+dz, material, name
4093 else:
4094 yield z+depthmantle, material, name
4096 z += dz
4099def write_nd_model_fh(mod, fh):
4100 def fmt(z, mat):
4101 rstr = ' '.join(
4102 util.gform(x, 4)
4103 for x in (
4104 z/1000.,
4105 mat.vp/1000.,
4106 mat.vs/1000.,
4107 mat.rho/1000.,
4108 mat.qp, mat.qs))
4109 if not mat._has_default_burgers():
4110 rstr += ' '.join(
4111 util.gform(x, 4)
4112 for x in (
4113 mat.burger_eta1,
4114 mat.burger_eta2,
4115 mat.burger_valpha))
4116 return rstr.rstrip() + '\n'
4118 translate = {
4119 'moho': 'mantle',
4120 'cmb': 'outer-core',
4121 'icb': 'inner-core'}
4123 last = None
4124 for element in mod.elements():
4125 if isinstance(element, Interface):
4126 if element.name is not None:
4127 n = translate.get(element.name, element.name)
4128 fh.write('%s\n' % n)
4130 elif isinstance(element, Layer):
4131 if not isinstance(last, Layer):
4132 fh.write(fmt(element.ztop, element.mtop))
4134 fh.write(fmt(element.zbot, element.mbot))
4136 last = element
4138 if not isinstance(last, Layer):
4139 fh.write(fmt(last.z, last.mbelow))
4142def write_nd_model_str(mod):
4143 f = StringIO()
4144 write_nd_model_fh(mod, f)
4145 return f.getvalue()
4148def write_nd_model(mod, fn):
4149 with open(fn, 'w') as f:
4150 write_nd_model_fh(mod, f)
4153def builtin_models():
4154 return sorted([
4155 os.path.splitext(os.path.basename(x))[0]
4156 for x in glob.glob(builtin_model_filename('*'))])
4159def builtin_model_filename(modelname):
4160 return util.data_file(os.path.join('earthmodels', modelname+'.nd'))
4163def load_model(fn='ak135-f-continental.m', format='nd', crust2_profile=None):
4164 '''
4165 Load layered earth model from file.
4167 :param fn: filename
4168 :param format: format
4169 :param crust2_profile: ``(lat, lon)`` or
4170 :py:class:`pyrocko.crust2x2.Crust2Profile` object, merge model with
4171 crustal profile. If ``fn`` is forced to be ``None`` only the converted
4172 CRUST2.0 profile is returned.
4173 :returns: object of type :py:class:`LayeredModel`
4175 The following formats are currently supported:
4177 ============== ===========================================================
4178 format description
4179 ============== ===========================================================
4180 ``'nd'`` 'named discontinuity' format used by the TauP programs
4181 ``'hyposat'`` format used by the HYPOSAT location program
4182 ============== ===========================================================
4184 The naming of interfaces is translated from the file format's native naming
4185 to Cake's own convention (See :py:func:`read_nd_model` and
4186 :py:func:`read_hyposat_model` for details). Cake likes the following
4187 internal names: ``'conrad'``, ``'moho'``, ``'cmb'`` (core-mantle boundary),
4188 ``'icb'`` (inner core boundary).
4189 '''
4191 if fn is not None:
4192 if format == 'nd':
4193 if not os.path.exists(fn) and fn in builtin_models():
4194 fn = builtin_model_filename(fn)
4195 reader = read_nd_model(fn)
4196 elif format == 'hyposat':
4197 reader = read_hyposat_model(fn)
4198 else:
4199 assert False, 'unsupported model format'
4201 mod = LayeredModel.from_scanlines(reader)
4202 if crust2_profile is not None:
4203 return mod.replaced_crust(crust2_profile)
4205 return mod
4207 else:
4208 assert crust2_profile is not None
4209 profile = anything_to_crust2_profile(crust2_profile)
4210 return LayeredModel.from_scanlines(
4211 from_crust2x2_profile(profile))
4214def castagna_vs_to_vp(vs):
4215 '''
4216 Calculate vp from vs using Castagna's relation.
4218 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4219 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4220 1985]
4222 vp = 1.16 * vs + 1360 [m/s]
4224 :param vs: S-wave velocity [m/s]
4225 :returns: P-wave velocity [m/s]
4226 '''
4228 return vs*1.16 + 1360.0
4231def castagna_vp_to_vs(vp):
4232 '''
4233 Calculate vp from vs using Castagna's relation.
4235 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4236 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4237 1985]
4239 vp = 1.16 * vs + 1360 [m/s]
4241 :param vp: P-wave velocity [m/s]
4242 :returns: S-wave velocity [m/s]
4243 '''
4245 return (vp - 1360.0) / 1.16
4248def evenize(x, y, minsize=10):
4249 if x.size < minsize:
4250 return x
4251 ry = (y.max()-y.min())
4252 if ry == 0:
4253 return x
4254 dx = (x[1:] - x[:-1])/(x.max()-x.min())
4255 dy = (y[1:] + y[:-1])/ry
4257 s = num.zeros(x.size)
4258 s[1:] = num.cumsum(num.sqrt(dy**2 + dx**2))
4259 s2 = num.linspace(0, s[-1], x.size)
4260 x2 = num.interp(s2, s, x)
4261 x2[0] = x[0]
4262 x2[-1] = x[-1]
4263 return x2
4266def filled(v, *args, **kwargs):
4267 '''
4268 Create NumPy array filled with given value.
4270 This works like :py:func:`numpy.ones` but initializes the array with ``v``
4271 instead of ones.
4272 '''
4273 x = num.empty(*args, **kwargs)
4274 x.fill(v)
4275 return x
4278def next_or_none(i):
4279 try:
4280 return next(i)
4281 except StopIteration:
4282 return None
4285def reci_or_none(x):
4286 try:
4287 return 1./x
4288 except ZeroDivisionError:
4289 return None
4292def mult_or_none(a, b):
4293 if a is None or b is None:
4294 return None
4295 return a*b
4298def min_not_none(a, b):
4299 if a is None:
4300 return b
4301 if b is None:
4302 return a
4303 return min(a, b)
4306def xytups(xx, yy):
4307 d = []
4308 for x, y in zip(xx, yy):
4309 if num.isfinite(y):
4310 d.append((x, y))
4311 return d
4314def interp(x, xp, fp, monoton):
4315 if monoton == 1:
4316 return xytups(
4317 x, num.interp(x, xp, fp, left=num.nan, right=num.nan))
4318 elif monoton == -1:
4319 return xytups(
4320 x, num.interp(x, xp[::-1], fp[::-1], left=num.nan, right=num.nan))
4321 else:
4322 fs = []
4323 for xv in x:
4324 indices = num.where(num.logical_or(
4325 num.logical_and(xp[:-1] >= xv, xv > xp[1:]),
4326 num.logical_and(xp[:-1] <= xv, xv < xp[1:])))[0]
4328 for i in indices:
4329 xr = (xv - xp[i])/(xp[i+1]-xp[i])
4330 fv = xr*fp[i] + (1.-xr)*fp[i+1]
4331 fs.append((xv, fv))
4333 return fs
4336def float_or_none(x):
4337 if x is not None:
4338 return float(x)
4341def parstore_float(thelocals, obj, *args):
4342 for k, v in thelocals.items():
4343 if k != 'self' and (not args or k in args):
4344 setattr(obj, k, float_or_none(v))