Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/cake.py: 85%
2041 statements
« prev ^ index » next coverage.py v7.6.0, created at 2025-12-04 10:41 +0000
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
70P = 1
71'''
72Constant indicating P wave propagation.
73'''
75S = 2
76'''
77Constant indicating S wave propagation
78'''
80DOWN = 4
81'''
82Constant indicating downward direction.
83'''
85UP = -4
86'''
87Constant indicating upward direction.
88'''
90DEFAULT_BURGERS = (0., 0., 1.)
92earthradius = config.config().earthradius
94r2d = 180./math.pi
95d2r = 1./r2d
96km = 1000.
97d2m = d2r*earthradius
98m2d = 1./d2m
99sprad2spm = 1.0/(r2d*d2m)
100sprad2spkm = 1.0/(r2d*d2m/km)
101spm2sprad = 1.0/sprad2spm
102spkm2sprad = 1.0/sprad2spkm
105class CakeError(Exception):
106 '''
107 Base class for exceptions raised in :py:mod:`pyrocko.cake`.
108 '''
109 pass
112class InvalidArguments(CakeError):
113 '''
114 Invalid arguments.
115 '''
116 pass
119class Material(object):
120 '''
121 Isotropic elastic material.
123 :param vp: P-wave velocity [m/s]
124 :param vs: S-wave velocity [m/s]
125 :param rho: density [kg/m^3]
126 :param qp: P-wave attenuation Qp
127 :param qs: S-wave attenuation Qs
128 :param poisson: Poisson ratio
129 :param lame: tuple with Lame parameter `lambda` and `shear modulus` [Pa]
130 :param qk: bulk attenuation Qk
131 :param qmu: shear attenuation Qmu
133 :param burgers: Burgers rheology paramerters as `tuple`.
134 `transient viscosity` [Pa], <= 0 means infinite value,
135 `steady-state viscosity` [Pa] and `alpha`, the ratio between the
136 effective and unreleaxed shear modulus, mu1/(mu1 + mu2).
137 :type burgers: tuple
139 If no velocities and no lame parameters are given, standard crustal values
140 of vp = 5800 m/s and vs = 3200 m/s are used. If no Q values are given,
141 standard crustal values of qp = 1456 and qs = 600 are used. If no Burgers
142 material parameters are given, transient and steady-state viscosities are
143 0 and alpha=1.
145 Everything is in SI units (m/s, Pa, kg/m^3) unless explicitly stated.
147 The main material properties are considered independant and are accessible
148 as attributes (it is allowed to assign to these):
150 .. py:attribute:: vp, vs, rho, qp, qs
152 Other material properties are considered dependant and can be queried by
153 instance methods.
154 '''
156 def __init__(
157 self, vp=None, vs=None, rho=2600., qp=None, qs=None, poisson=None,
158 lame=None, qk=None, qmu=None, burgers=None):
160 parstore_float(locals(), self, 'vp', 'vs', 'rho', 'qp', 'qs')
162 if vp is not None and vs is not None:
163 if poisson is not None or lame is not None:
164 raise InvalidArguments(
165 'If vp and vs are given, poisson ratio and lame paramters '
166 'should not be given.')
168 elif vp is None and vs is None and lame is None:
169 self.vp = 5800.
170 if poisson is None:
171 poisson = 0.25
172 self.vs = self.vp / math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson))
174 elif vp is None and vs is None and lame is not None:
175 if poisson is not None:
176 raise InvalidArguments(
177 'Poisson ratio should not be given, when lame parameters '
178 'are given.')
180 lam, mu = float(lame[0]), float(lame[1])
181 self.vp = math.sqrt((lam + 2.0*mu)/rho)
182 self.vs = math.sqrt(mu/rho)
184 elif vp is not None and vs is None:
185 if poisson is None:
186 poisson = 0.25
188 if lame is not None:
189 raise InvalidArguments(
190 'If vp is given, Lame parameters should not be given.')
192 poisson = float(poisson)
193 self.vs = vp / math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson))
195 elif vp is None and vs is not None:
196 if poisson is None:
197 poisson = 0.25
198 if lame is not None:
199 raise InvalidArguments(
200 'If vs is given, Lame parameters should not be given.')
202 poisson = float(poisson)
203 self.vp = vs * math.sqrt(2.0*(1.0-poisson)/(1.0-2.0*poisson))
205 else:
206 raise InvalidArguments(
207 'Invalid combination of input parameters in material '
208 'definition.')
210 if qp is not None or qs is not None:
211 if not (qk is None and qmu is None):
212 raise InvalidArguments(
213 'If qp or qs are given, qk and qmu should not be given.')
215 if qp is None:
216 if self.vs != 0.0:
217 s = (4.0/3.0)*(self.vs/self.vp)**2
218 self.qp = self.qs / s
219 else:
220 self.qp = 1456.
222 if qs is None:
223 if self.vs != 0.0:
224 s = (4.0/3.0)*(self.vs/self.vp)**2
225 self.qs = self.qp * s
226 else:
227 self.vs = 600.
229 elif qp is None and qs is None and qk is None and qmu is None:
230 if self.vs == 0.:
231 self.qs = 0.
232 self.qp = 5782e4
233 else:
234 self.qs = 600.
235 s = (4.0/3.0)*(self.vs/self.vp)**2
236 self.qp = self.qs/s
238 elif qp is None and qs is None and qk is not None and qmu is not None:
239 s = (4.0/3.0)*(self.vs/self.vp)**2
240 if qmu == 0. and self.vs == 0.:
241 self.qp = qk
242 else:
243 if num.isinf(qk):
244 self.qp = qmu/s
245 else:
246 self.qp = 1.0 / (s/qmu + (1.0-s)/qk)
247 self.qs = qmu
248 else:
249 raise InvalidArguments(
250 'Invalid combination of input parameters in material '
251 'definition.')
253 if burgers is None:
254 burgers = DEFAULT_BURGERS
256 self.burger_eta1 = burgers[0]
257 self.burger_eta2 = burgers[1]
258 self.burger_valpha = burgers[2]
260 def astuple(self):
261 '''
262 Get independant material properties as a tuple.
264 Returns a tuple with ``(vp, vs, rho, qp, qs)``.
265 '''
266 return self.vp, self.vs, self.rho, self.qp, self.qs
268 def __eq__(self, other):
269 return self.astuple() == other.astuple()
271 def lame(self):
272 '''
273 Get Lame's parameter lambda and shear modulus.
274 '''
275 mu = self.vs**2 * self.rho
276 lam = self.vp**2 * self.rho - 2.0*mu
277 return lam, mu
279 def lame_lambda(self):
280 '''
281 Get Lame's parameter lambda.
283 Returned units are [Pa].
284 '''
285 lam, _ = self.lame()
286 return lam
288 def shear_modulus(self):
289 '''
290 Get shear modulus.
292 Returned units are [Pa].
293 '''
294 return self.vs**2 * self.rho
296 def poisson(self):
297 '''
298 Get Poisson's ratio.
299 '''
300 lam, mu = self.lame()
301 return lam / (2.0*(lam+mu))
303 def bulk(self):
304 '''
305 Get bulk modulus.
306 '''
307 lam, mu = self.lame()
308 return lam + 2.0*mu/3.0
310 def youngs(self):
311 '''
312 Get Young's modulus.
313 '''
314 lam, mu = self.lame()
315 return mu * (3.0*lam + 2.0*mu) / (lam+mu)
317 def vp_vs_ratio(self):
318 '''
319 Get vp/vs ratio.
320 '''
321 return self.vp/self.vs
323 def qmu(self):
324 '''
325 Get shear attenuation coefficient Qmu.
326 '''
327 return self.qs
329 def qk(self):
330 '''
331 Get bulk attenuation coefficient Qk.
332 '''
333 if self.vs == 0. and self.qs == 0.:
334 return self.qp
335 else:
336 s = (4.0/3.0)*(self.vs/self.vp)**2
337 denom = (1/self.qp - s/self.qs)
338 if denom <= 0.0:
339 return num.inf
340 else:
341 return (1.-s)/(1.0/self.qp - s/self.qs)
343 def burgers(self):
344 '''
345 Get Burger parameters.
346 '''
347 return self.burger_eta1, self.burger_eta2, self.burger_valpha
349 def _rayleigh_equation(self, cr):
350 cr_a = (cr/self.vp)**2
351 cr_b = (cr/self.vs)**2
352 if cr_a > 1.0 or cr_b > 1.0:
353 return None
355 return (2.0-cr_b)**2 - 4.0 * math.sqrt(1.0-cr_a) * math.sqrt(1.0-cr_b)
357 def rayleigh(self):
358 '''
359 Get Rayleigh velocity assuming a homogenous halfspace.
361 Returned units are [m/s].
362 '''
363 return bisect(self._rayleigh_equation, 0.001*self.vs, self.vs)
365 def _has_default_burgers(self):
366 if self.burger_eta1 == DEFAULT_BURGERS[0] and \
367 self.burger_eta2 == DEFAULT_BURGERS[1] and \
368 self.burger_valpha == DEFAULT_BURGERS[2]:
369 return True
370 return False
372 def describe(self):
373 '''
374 Get a readable listing of the material properties.
375 '''
376 template = '''
377P wave velocity [km/s] : %12g
378S wave velocity [km/s] : %12g
379P/S wave vel. ratio : %12g
380Lame lambda [GPa] : %12g
381Lame shear modulus [GPa] : %12g
382Poisson ratio : %12g
383Bulk modulus [GPa] : %12g
384Young's modulus [GPa] : %12g
385Rayleigh wave vel. [km/s] : %12g
386Density [g/cm**3] : %12g
387Qp P-wave attenuation : %12g
388Qs S-wave attenuation (Qmu) : %12g
389Qk bulk attenuation : %12g
390transient viscos., eta1 [GPa] : %12g
391st.-state viscos., eta2 [GPa] : %12g
392relaxation: valpha : %12g
393'''.strip()
395 return template % (
396 self.vp/km,
397 self.vs/km,
398 self.vp/self.vs,
399 self.lame_lambda()*1e-9,
400 self.shear_modulus()*1e-9,
401 self.poisson(),
402 self.bulk()*1e-9,
403 self.youngs()*1e-9,
404 self.rayleigh()/km,
405 self.rho/km,
406 self.qp,
407 self.qs,
408 self.qk(),
409 self.burger_eta1*1e-9,
410 self.burger_eta2*1e-9,
411 self.burger_valpha)
413 def __str__(self):
414 vp, vs, rho, qp, qs = self.astuple()
415 return '%10g km/s %10g km/s %10g g/cm^3 %10g %10g' % (
416 vp/km, vs/km, rho/km, qp, qs)
418 def __repr__(self):
419 return 'Material(vp=%s, vs=%s, rho=%s, qp=%s, qs=%s)' % \
420 tuple(repr(x) for x in (
421 self.vp, self.vs, self.rho, self.qp, self.qs))
424class Leg(object):
425 '''
426 Represents a continuous piece of wave propagation in a phase definition.
428 **Attributes:**
430 To be considered as read-only.
432 .. py:attribute:: departure
434 One of the constants :py:const:`UP` or :py:const:`DOWN` indicating
435 upward or downward departure.
437 .. py:attribute:: mode
439 One of the constants :py:const:`P` or :py:const:`S`, indicating the
440 propagation mode.
442 .. py:attribute:: depthmin
444 ``None``, a number (a depth in [m]) or a string (an interface name),
445 minimum depth.
447 .. py:attribute:: depthmax
449 ``None``, a number (a depth in [m]) or a string (an interface name),
450 maximum depth.
452 '''
454 def __init__(self, departure=None, mode=None):
455 self.departure = departure
456 self.mode = mode
457 self.depthmin = None
458 self.depthmax = None
460 def set_depthmin(self, depthmin):
461 self.depthmin = depthmin
463 def set_depthmax(self, depthmax):
464 self.depthmax = depthmax
466 def __str__(self):
467 def sd(d):
468 if isinstance(d, float):
469 return '%g km' % (d/km)
470 else:
471 return 'interface %s' % d
473 s = '%s mode propagation, departing %s' % (
474 smode(self.mode).upper(), {
475 UP: 'upward', DOWN: 'downward'}[self.departure])
477 sc = []
478 if self.depthmax is not None:
479 sc.append('deeper than %s' % sd(self.depthmax))
480 if self.depthmin is not None:
481 sc.append('shallower than %s' % sd(self.depthmin))
483 if sc:
484 s = s + ' (may not propagate %s)' % ' or '.join(sc)
486 return s
489class InvalidKneeDef(CakeError):
490 '''
491 Invalid definition of a cake :py:class:`Knee`.
492 '''
493 pass
496class Knee(object):
497 '''
498 Represents a change in wave propagation within a :py:class:`PhaseDef`.
500 **Attributes:**
502 To be considered as read-only.
504 .. py:attribute:: depth
506 Depth at which the conversion/reflection should happen. this can be
507 a string or a number.
509 .. py:attribute:: direction
511 One of the constants :py:const:`UP` or :py:const:`DOWN` to indicate
512 the incoming direction.
514 .. py:attribute:: in_mode
516 One of the constants :py:const:`P` or :py:const:`S` to indicate the
517 type of mode of the incoming wave.
519 .. py:attribute:: out_mode
521 One of the constants :py:const:`P` or :py:const:`S` to indicate the
522 type of mode of the outgoing wave.
524 .. py:attribute:: conversion
526 Boolean, whether there is a mode conversion involved.
528 .. py:attribute:: reflection
530 Boolean, whether there is a reflection involved.
532 .. py:attribute:: headwave
534 Boolean, whether there is headwave propagation involved.
536 '''
538 defaults = dict(
539 depth='surface',
540 direction=UP,
541 conversion=True,
542 reflection=False,
543 headwave=False,
544 in_setup_state=True)
546 defaults_surface = dict(
547 depth='surface',
548 direction=UP,
549 conversion=False,
550 reflection=True,
551 headwave=False,
552 in_setup_state=True)
554 def __init__(self, *args):
555 if args:
556 (self.depth, self.direction, self.reflection, self.in_mode,
557 self.out_mode) = args
559 self.conversion = self.in_mode != self.out_mode
560 self.in_setup_state = False
562 def default(self, k):
563 depth = self.__dict__.get('depth', 'surface')
564 if depth == 'surface':
565 return Knee.defaults_surface[k]
566 else:
567 return Knee.defaults[k]
569 def __setattr__(self, k, v):
570 if self.in_setup_state and k in self.__dict__:
571 raise InvalidKneeDef('%s has already been set' % k)
572 else:
573 self.__dict__[k] = v
575 def __getattr__(self, k):
576 if k.startswith('__'):
577 raise AttributeError(k)
579 if k not in self.__dict__:
580 return self.default(k)
582 def set_modes(self, in_leg, out_leg):
584 if out_leg.departure == UP and (
585 (self.direction == UP) == self.reflection):
587 raise InvalidKneeDef(
588 'cannot enter %s from %s and emit ray upwards' % (
589 ['conversion', 'reflection'][self.reflection],
590 {UP: 'below', DOWN: 'above'}[self.direction]))
592 if out_leg.departure == DOWN and (
593 (self.direction == DOWN) == self.reflection):
595 raise InvalidKneeDef(
596 'cannot enter %s from %s and emit ray downwards' % (
597 ['conversion', 'reflection'][self.reflection],
598 {UP: 'below', DOWN: 'above'}[self.direction]))
600 self.in_mode = in_leg.mode
601 self.out_mode = out_leg.mode
603 def at_surface(self):
604 return self.depth == 'surface'
606 def matches(self, discontinuity, mode, direction):
607 '''
608 Check whether it is relevant to a given combination of interface,
609 propagation mode, and direction.
610 '''
612 if isinstance(self.depth, float):
613 if abs(self.depth - discontinuity.z) > ZEPS:
614 return False
615 else:
616 if discontinuity.name != self.depth:
617 return False
619 return self.direction == direction and self.in_mode == mode
621 def out_direction(self):
622 '''
623 Get outgoing direction.
625 Returns one of the constants :py:const:`UP` or :py:const:`DOWN`.
626 '''
628 if self.reflection:
629 return - self.direction
630 else:
631 return self.direction
633 def __str__(self):
634 x = []
635 if self.reflection:
636 if self.at_surface():
637 x.append('surface')
638 else:
639 if not self.headwave:
640 if self.direction == UP:
641 x.append('underside')
642 else:
643 x.append('upperside')
645 if self.headwave:
646 x.append('headwave propagation along')
647 elif self.reflection and self.conversion:
648 x.append('reflection with conversion from %s to %s' % (
649 smode(self.in_mode).upper(), smode(self.out_mode).upper()))
650 if not self.at_surface():
651 x.append('at')
652 elif self.reflection:
653 x.append('reflection')
654 if not self.at_surface():
655 x.append('at')
656 elif self.conversion:
657 x.append('conversion from %s to %s at' % (
658 smode(self.in_mode).upper(), smode(self.out_mode).upper()))
659 else:
660 x.append('passing through')
662 if isinstance(self.depth, float):
663 x.append('interface in %g km depth' % (self.depth/1000.))
664 else:
665 if not self.at_surface():
666 x.append('%s' % self.depth)
668 if not self.reflection:
669 if self.direction == UP:
670 x.append('on upgoing path')
671 else:
672 x.append('on downgoing path')
674 return ' '.join(x)
677class Head(Knee):
678 def __init__(self, *args):
679 if args:
680 z, in_direction, mode = args
681 Knee.__init__(self, z, in_direction, True, mode, mode)
682 else:
683 Knee.__init__(self)
685 def __str__(self):
686 x = ['propagation as headwave']
687 if isinstance(self.depth, float):
688 x.append('at interface in %g km depth' % (self.depth/1000.))
689 else:
690 x.append('at %s' % self.depth)
692 return ' '.join(x)
695class UnknownClassicPhase(CakeError):
696 '''
697 Raised when an invalid classic phase name has been specified.'
698 '''
699 def __init__(self, phasename):
700 self.phasename = phasename
702 def __str__(self):
703 return 'Unknown classic phase name: %s' % self.phasename
706class PhaseDefParseError(CakeError):
707 '''
708 Exception raised when an error occures during parsing of a phase
709 definition string.
710 '''
712 def __init__(self, definition, position, exception):
713 self.definition = definition
714 self.position = position
715 self.exception = exception
717 def __str__(self):
718 return 'Invalid phase definition: "%s" (at character %i: %s)' % (
719 self.definition, self.position+1, str(self.exception))
722class PhaseDef(object):
724 '''
725 Definition of a seismic phase arrival, based on ray propagation path.
727 :param definition: string representation of the phase in Cake's phase
728 syntax
730 Seismic phases are conventionally named e.g. P, Pn, PP, PcP, etc. In Cake,
731 a slightly different terminology is adapted, which allows to specify
732 arbitrary conversion/reflection histories for seismic ray paths. The
733 conventions used here are inspired by those used in the TauP toolkit, but
734 are not completely compatible with those.
736 The definition of a seismic ray propagation path in Cake's phase syntax is
737 a string consisting of an alternating sequence of *legs* and *knees*.
739 A *leg* represents seismic wave propagation without any conversions,
740 encountering only super-critical reflections. Legs are denoted by ``P``,
741 ``p``, ``S``, or ``s``. The capital letters are used when the take-off of
742 the *leg* is in downward direction, while the lower case letters indicate a
743 take-off in upward direction.
745 A *knee* is an interaction with an interface. It can be a mode conversion,
746 a reflection, or propagation as a headwave or diffracted wave.
748 * conversion is simply denoted as: ``(INTERFACE)`` or ``DEPTH``
749 * upperside reflection: ``v(INTERFACE)`` or ``vDEPTH``
750 * underside reflection: ``^(INTERFACE)`` or ``^DEPTH``
751 * normal kind headwave or diffracted wave: ``v_(INTERFACE)`` or
752 ``v_DEPTH``
754 The interface may be given by name or by depth: INTERFACE is the name of an
755 interface defined in the model, DEPTH is the depth of an interface in
756 [km] (the interface closest to that depth is chosen). If two legs appear
757 consecutively without an explicit *knee*, surface interaction is assumed.
759 The phase definition may end with a backslash ``\\``, to indicate that the
760 ray should arrive at the receiver from above instead of from below. It is
761 possible to restrict the maximum and minimum depth of a *leg* by appending
762 ``<(INTERFACE)`` or ``<DEPTH`` or ``>(INTERFACE)`` or ``>DEPTH`` after the
763 leg character, respectively.
765 **Examples:**
767 * ``P`` - like the classical P, but includes PKP, PKIKP, Pg
768 * ``P<(moho)`` - like classical Pg, but must leave source downwards
769 * ``pP`` - leaves source upward, reflects at surface, then travels as P
770 * ``P(moho)s`` - conversion from P to S at the Moho on upgoing path
771 * ``P(moho)S`` - conversion from P to S at the Moho on downgoing path
772 * ``Pv12p`` - P with reflection at 12 km deep interface (or the
773 interface closest to that)
774 * ``Pv_(moho)p`` - classical Pn
775 * ``Pv_(cmb)p`` - classical Pdiff
776 * ``P^(conrad)P`` - underside reflection of P at the Conrad
777 discontinuity
779 **Usage:**
781 >>> from pyrocko.cake import PhaseDef
782 # must escape the backslash
783 >>> my_crazy_phase = PhaseDef('pPv(moho)sP\\\\')
784 >>> print my_crazy_phase
785 Phase definition "pPv(moho)sP\":
786 - P mode propagation, departing upward
787 - surface reflection
788 - P mode propagation, departing downward
789 - upperside reflection with conversion from P to S at moho
790 - S mode propagation, departing upward
791 - surface reflection with conversion from S to P
792 - P mode propagation, departing downward
793 - arriving at target from above
795 .. note::
797 (1) These conventions might be extended in a way to allow to fix wave
798 propagation to SH mode, possibly by specifying SH, or a single
799 character (e.g. H) instead of S. This would be benificial for the
800 selection of conversion and reflection coefficients, which
801 currently only deal with the P-SV case.
802 '''
804 allowed_characters_pattern = r'[0-9a-zA-Z_()<>^v\\.]+'
805 allowed_characters_pattern_classic = r'[a-zA-Z0-9]+'
807 @staticmethod
808 def classic_definitions():
809 defs = {}
810 # PmP, PmS, PcP, PcS, SmP, ...
811 for r in 'mc':
812 for a, b in 'PP PS SS SP'.split():
813 defs[a+r+b] = [
814 '%sv(%s)%s' % (a, {'m': 'moho', 'c': 'cmb'}[r], b.lower())]
816 # Pg, P, S, Sg
817 for a in 'PS':
818 defs[a+'g'] = ['%s<(moho)' % x for x in (a, a.lower())]
819 defs[a] = ['%s<(cmb)(moho)%s' % (x, x.lower()) for x in (
820 a, a.lower())]
822 defs[a.lower()] = [a.lower()]
824 for a, b in 'PP PS SS SP'.split():
825 defs[a+'K'+b] = ['%s(cmb)P<(icb)(cmb)%s' % (a, b.lower())]
826 defs[a+'KIK'+b] = ['%s(cmb)P(icb)P(icb)p(cmb)%s' % (a, b.lower())]
827 defs[a+'KJK'+b] = ['%s(cmb)P(icb)S(icb)p(cmb)%s' % (a, b.lower())]
828 defs[a+'KiK'+b] = ['%s(cmb)Pv(icb)p(cmb)%s' % (a, b.lower())]
830 # PP, SS, PS, SP, PPP, ...
831 for a in 'PS':
832 for b in 'PS':
833 for c in 'PS':
834 defs[a+b+c] = [''.join(defs[x][0] for x in a+b+c)]
836 defs[a+b] = [''.join(defs[x][0] for x in a+b)]
838 # Pc, Pdiff, Sc, ...
839 for x in 'PS':
840 defs[x+'c'] = defs[x+'diff'] = [x+'v_(cmb)'+x.lower()]
841 defs[x+'n'] = [x+'v_(moho)'+x.lower()]
843 # depth phases
844 for k in list(defs.keys()):
845 if k not in 'ps':
846 for x in 'ps':
847 defs[x+k] = [x + defs[k][0]]
849 return defs
851 @staticmethod
852 def classic(phasename):
853 '''
854 Get phase definitions based on classic phase name.
856 :param phasename: classic name of a phase
857 :returns: list of PhaseDef objects
859 This returns a list of PhaseDef objects, because some classic phases
860 (like e.g. Pg) can only be represented by two Cake style PhaseDef
861 objects (one with downgoing and one with upgoing first leg).
862 '''
864 defs = PhaseDef.classic_definitions()
865 if phasename not in defs:
866 raise UnknownClassicPhase(phasename)
868 return [PhaseDef(d, classicname=phasename) for d in defs[phasename]]
870 def __init__(self, definition=None, classicname=None):
872 state = 0
873 sdepth = ''
874 sinterface = ''
875 depthmax = depthmin = None
876 depthlim = None
877 depthlimtype = None
878 sdepthlim = ''
879 events = []
880 direction_stop = UP
881 need_leg = True
882 ic = 0
883 if definition is not None:
884 knee = Knee()
885 try:
886 for ic, c in enumerate(definition):
888 if state in (0, 1):
890 if c in '0123456789.':
891 need_leg = True
892 state = 1
893 sdepth += c
894 continue
896 elif state == 1:
897 knee.depth = float(sdepth)*1000.
898 state = 0
900 if state == 2:
901 if c == ')':
902 knee.depth = sinterface
903 state = 0
904 else:
905 sinterface += c
907 continue
909 if state in (3, 4):
911 if state == 3:
912 if c in '0123456789.':
913 sdepthlim += c
914 continue
915 elif c == '(':
916 state = 4
917 continue
918 else:
919 depthlim = float(sdepthlim)*1000.
920 if depthlimtype == '<':
921 depthmax = depthlim
922 else:
923 depthmin = depthlim
924 state = 0
926 elif state == 4:
927 if c == ')':
928 depthlim = sdepthlim
929 if depthlimtype == '<':
930 depthmax = depthlim
931 else:
932 depthmin = depthlim
933 state = 0
934 continue
935 else:
936 sdepthlim += c
937 continue
939 if state == 0:
941 if c == '(':
942 need_leg = True
943 state = 2
944 continue
946 elif c in '<>':
947 state = 3
948 depthlim = None
949 sdepthlim = ''
950 depthlimtype = c
951 continue
953 elif c in 'psPS':
954 leg = Leg()
955 if c in 'ps':
956 leg.departure = UP
957 else:
958 leg.departure = DOWN
959 leg.mode = imode(c)
961 if events:
962 in_leg = events[-1]
963 if depthmin is not None:
964 in_leg.set_depthmin(depthmin)
965 depthmin = None
966 if depthmax is not None:
967 in_leg.set_depthmax(depthmax)
968 depthmax = None
970 if in_leg.mode != leg.mode:
971 knee.conversion = True
972 else:
973 knee.conversion = False
975 if not knee.reflection:
976 if c in 'ps':
977 knee.direction = UP
978 else:
979 knee.direction = DOWN
981 knee.set_modes(in_leg, leg)
982 knee.in_setup_state = False
983 events.append(knee)
984 knee = Knee()
985 sdepth = ''
986 sinterface = ''
988 events.append(leg)
989 need_leg = False
990 continue
992 elif c == '^':
993 need_leg = True
994 knee.direction = UP
995 knee.reflection = True
996 continue
998 elif c == 'v':
999 need_leg = True
1000 knee.direction = DOWN
1001 knee.reflection = True
1002 continue
1004 elif c == '_':
1005 need_leg = True
1006 knee.headwave = True
1007 continue
1009 elif c == '\\':
1010 direction_stop = DOWN
1011 continue
1013 else:
1014 raise PhaseDefParseError(
1015 definition, ic, 'invalid character: "%s"' % c)
1017 if state == 3:
1018 depthlim = float(sdepthlim)*1000.
1019 if depthlimtype == '<':
1020 depthmax = depthlim
1021 else:
1022 depthmin = depthlim
1023 state = 0
1025 except (ValueError, InvalidKneeDef) as e:
1026 raise PhaseDefParseError(definition, ic, e)
1028 if state != 0 or need_leg:
1029 raise PhaseDefParseError(
1030 definition, ic, 'unfinished expression')
1032 if events and depthmin is not None:
1033 events[-1].set_depthmin(depthmin)
1034 if events and depthmax is not None:
1035 events[-1].set_depthmax(depthmax)
1037 self._definition = definition
1038 self._classicname = classicname
1039 self._events = events
1040 self._direction_stop = direction_stop
1042 def __iter__(self):
1043 for ev in self._events:
1044 yield ev
1046 def append(self, ev):
1047 self._events.append(ev)
1049 def first_leg(self):
1050 '''
1051 Get the first leg in phase definition.
1052 '''
1053 return self._events[0]
1055 def last_leg(self):
1056 '''
1057 Get the last leg in phase definition.
1058 '''
1059 return self._events[-1]
1061 def legs(self):
1062 '''
1063 Iterate over the continuous pieces of wave propagation (legs) defined
1064 within this phase definition.
1065 '''
1067 return (leg for leg in self if isinstance(leg, Leg))
1069 def knees(self):
1070 '''
1071 Iterate over conversions and reflections (knees) defined within this
1072 phase definition.
1073 '''
1074 return (knee for knee in self if isinstance(knee, Knee))
1076 def definition(self):
1077 '''
1078 Get original definition of the phase.
1079 '''
1080 return self._definition
1082 def given_name(self):
1083 '''
1084 Get entered classic name if any, or original definition of the phase.
1085 '''
1087 if self._classicname:
1088 return self._classicname
1089 else:
1090 return self._definition
1092 def direction_start(self):
1093 return self.first_leg().departure
1095 def direction_stop(self):
1096 return self._direction_stop
1098 def headwave_knee(self):
1099 for el in self:
1100 if type(el) is Knee and el.headwave:
1101 return el
1102 return None
1104 def used_repr(self):
1105 '''
1106 Translate into textual representation (cake phase syntax).
1107 '''
1108 def strdepth(x):
1109 if isinstance(x, float):
1110 return '%g' % (x/1000.)
1111 else:
1112 return '(%s)' % x
1114 x = []
1115 for el in self:
1116 if type(el) is Leg:
1117 if el.departure == UP:
1118 x.append(smode(el.mode).lower())
1119 else:
1120 x.append(smode(el.mode).upper())
1122 if el.depthmax is not None:
1123 x.append('<'+strdepth(el.depthmax))
1125 if el.depthmin is not None:
1126 x.append('>'+strdepth(el.depthmin))
1128 elif type(el) is Knee:
1129 if el.reflection and not el.at_surface():
1130 if el.direction == DOWN:
1131 x.append('v')
1132 else:
1133 x.append('^')
1134 if el.headwave:
1135 x.append('_')
1136 if not el.at_surface():
1137 x.append(strdepth(el.depth))
1139 elif type(el) is Head:
1140 x.append('_')
1141 x.append(strdepth(el.depth))
1143 if self._direction_stop == DOWN:
1144 x.append('\\')
1146 return ''.join(x)
1148 def __repr__(self):
1149 if self._definition is not None:
1150 return "PhaseDef('%s')" % self._definition
1151 else:
1152 return "PhaseDef('%s')" % self.used_repr()
1154 def __str__(self):
1155 orig = ''
1156 used = self.used_repr()
1157 if self._definition != used:
1158 orig = ' (entered as "%s")' % self._definition
1160 sarrive = '\n - arriving at target from %s' % ('below', 'above')[
1161 self._direction_stop == DOWN]
1163 return 'Phase definition "%s"%s:\n - ' % (used, orig) + \
1164 '\n - '.join(str(ev) for ev in self) + sarrive
1166 def copy(self):
1167 '''
1168 Get a deep copy of it.
1169 '''
1170 return copy.deepcopy(self)
1173def to_phase_defs(phases):
1174 if isinstance(phases, (str, PhaseDef)):
1175 phases = [phases]
1177 phases_out = []
1178 for phase in phases:
1179 if isinstance(phase, str):
1180 phases_out.extend(PhaseDef(x.strip()) for x in phase.split(','))
1181 elif isinstance(phase, PhaseDef):
1182 phases_out.append(phase)
1183 else:
1184 raise PhaseDefParseError('invalid phase definition')
1186 return phases_out
1189def csswap(x):
1190 return cmath.sqrt(1.-x**2)
1193def psv_surface_ind(in_mode, out_mode):
1194 '''
1195 Get indices to select the appropriate element from scatter matrix for free
1196 surface.
1197 '''
1199 return (int(in_mode == S), int(out_mode == S))
1202def psv_surface(material, p, energy=False):
1203 '''
1204 Scatter matrix for free surface reflection/conversions.
1206 :param material: material, object of type :py:class:`Material`
1207 :param p: flat ray parameter [s/m]
1208 :param energy: bool, when ``True`` energy normalized coefficients are
1209 returned
1210 :returns: Scatter matrix
1212 The scatter matrix is ordered as follows::
1214 [[PP, PS],
1215 [SP, SS]]
1217 The formulas given in Aki & Richards are used.
1218 '''
1220 vp, vs, rho = material.vp, material.vs, material.rho
1221 sinphi = p * vp
1222 sinlam = p * vs
1223 cosphi = csswap(sinphi)
1224 coslam = csswap(sinlam)
1226 if vs == 0.0:
1227 scatter = num.array([[-1.0, 0.0], [0.0, 1.0]])
1229 else:
1230 vsp_term = (1.0/vs**2 - 2.0*p**2)
1231 pcc_term = 4.0 * p**2 * cosphi/vp * coslam/vs
1232 denom = vsp_term**2 + pcc_term
1234 scatter = num.array([
1235 [- vsp_term**2 + pcc_term, 4.0*p*coslam/vp*vsp_term],
1236 [4.0*p*cosphi/vs*vsp_term, vsp_term**2 - pcc_term]],
1237 dtype=complex) / denom
1239 if not energy:
1240 return scatter
1241 else:
1242 eps = 1e-16
1243 normvec = num.array([vp*rho*cosphi+eps, vs*rho*coslam+eps])
1244 escatter = scatter*num.conj(scatter) * num.real(
1245 (normvec[:, num.newaxis]) / (normvec[num.newaxis, :]))
1246 return num.real(escatter)
1249def psv_solid_ind(in_direction, out_direction, in_mode, out_mode):
1250 '''
1251 Get indices to select the appropriate element from scatter matrix for
1252 solid-solid interface.
1253 '''
1255 return (
1256 (out_direction == DOWN)*2 + (out_mode == S),
1257 (in_direction == UP)*2 + (in_mode == S))
1260def psv_solid(material1, material2, p, energy=False):
1261 '''
1262 Scatter matrix for solid-solid interface.
1264 :param material1: material above, object of type :py:class:`Material`
1265 :param material2: material below, object of type :py:class:`Material`
1266 :param p: flat ray parameter [s/m]
1267 :param energy: bool, when ``True`` energy normalized coefficients are
1268 returned
1269 :returns: Scatter matrix
1271 The scatter matrix is ordered as follows::
1273 [[P1P1, S1P1, P2P1, S2P1],
1274 [P1S1, S1S1, P2S1, S2S1],
1275 [P1P2, S1P2, P2P2, S2P2],
1276 [P1S2, S1S2, P2S2, S2S2]]
1278 The formulas given in Aki & Richards are used.
1279 '''
1281 vp1, vs1, rho1 = material1.vp, material1.vs, material1.rho
1282 vp2, vs2, rho2 = material2.vp, material2.vs, material2.rho
1284 sinphi1 = p * vp1
1285 cosphi1 = csswap(sinphi1)
1286 sinlam1 = p * vs1
1287 coslam1 = csswap(sinlam1)
1288 sinphi2 = p * vp2
1289 cosphi2 = csswap(sinphi2)
1290 sinlam2 = p * vs2
1291 coslam2 = csswap(sinlam2)
1293 # from aki and richards
1294 M = num.array([
1295 [-vp1*p, -coslam1, vp2*p, coslam2],
1296 [cosphi1, -vs1*p, cosphi2, -vs2*p],
1297 [2.0*rho1*vs1**2*p*cosphi1, rho1*vs1*(1.0-2.0*vs1**2*p**2),
1298 2.0*rho2*vs2**2*p*cosphi2, rho2*vs2*(1.0-2.0*vs2**2*p**2)],
1299 [-rho1*vp1*(1.0-2.0*vs1**2*p**2), 2.0*rho1*vs1**2*p*coslam1,
1300 rho2*vp2*(1.0-2.0*vs2**2*p**2), -2.0*rho2*vs2**2*p*coslam2]],
1301 dtype=complex)
1302 N = M.copy()
1303 N[0] *= -1.0
1304 N[3] *= -1.0
1306 scatter = num.dot(num.linalg.inv(M), N)
1308 if not energy:
1309 return scatter
1310 else:
1311 eps = 1e-16
1312 if vs1 == 0.:
1313 vs1 = vp1*1e-16
1314 if vs2 == 0.:
1315 vs2 = vp2*1e-16
1316 normvec = num.array([
1317 vp1*rho1*(cosphi1+eps), vs1*rho1*(coslam1+eps),
1318 vp2*rho2*(cosphi2+eps), vs2*rho2*(coslam2+eps)], dtype=complex)
1319 escatter = scatter*num.conj(scatter) * num.real(
1320 normvec[:, num.newaxis] / normvec[num.newaxis, :])
1322 return num.real(escatter)
1325class BadPotIntCoefs(CakeError):
1326 '''
1327 Raised when the potential interpolation for gradient layers has failed.
1328 '''
1329 pass
1332def potint_coefs(c1, c2, r1, r2): # r2 > r1
1333 eps = r2*1e-9
1334 if c1 == 0. and c2 == 0.:
1335 c1c2 = 1.
1336 else:
1337 c1c2 = c1/c2
1338 b = math.log(c1c2)/math.log((r1+eps)/r2)
1339 if abs(b) > 10.:
1340 raise BadPotIntCoefs()
1341 a = c1/(r1+eps)**b
1342 return a, b
1345def imode(s):
1346 if s.lower() == 'p':
1347 return P
1348 elif s.lower() == 's':
1349 return S
1352def smode(i):
1353 if i == P:
1354 return 'p'
1355 elif i == S:
1356 return 's'
1359class PathFailed(CakeError):
1360 '''
1361 Raised when the ray path computation failed.
1362 '''
1363 pass
1366class SurfaceReached(PathFailed):
1367 '''
1368 Raised when the ray hits the surface before completing the phase
1369 requirements.
1370 '''
1371 pass
1374class BottomReached(PathFailed):
1375 '''
1376 Raised when the ray exits the bottom of the model before completing the
1377 phase requirements.
1378 '''
1379 pass
1382class MaxDepthReached(PathFailed):
1383 '''
1384 Raised when the phase's maximum depth has been exceeded.
1385 '''
1386 pass
1389class MinDepthReached(PathFailed):
1390 '''
1391 Raised when the phase's minimum depth has been underrun.
1392 '''
1393 pass
1396class Trapped(PathFailed):
1397 '''
1398 Raised when the ray has been trapped and therefore cannot satisfy the phase
1399 requirements.
1400 '''
1401 pass
1404class NotPhaseConform(PathFailed):
1405 '''
1406 Raised when the phase requirements cannot be met.
1407 '''
1408 pass
1411class CannotPropagate(PathFailed):
1412 '''
1413 Raised when the phase requirements cannot be met.
1414 '''
1415 def __init__(self, direction, ilayer):
1416 PathFailed.__init__(self)
1417 self._direction = direction
1418 self._ilayer = ilayer
1420 def __str__(self):
1421 return 'Cannot enter layer %i from %s' % (
1422 self._ilayer, {
1423 UP: 'below',
1424 DOWN: 'above'}[self._direction])
1427class Layer(object):
1428 '''
1429 Representation of a layer in a layered earth model.
1431 :param ztop: depth of top of layer
1432 :param zbot: depth of bottom of layer
1433 :param name: name of layer (optional)
1435 Subclasses are: :py:class:`HomogeneousLayer` and :py:class:`GradientLayer`.
1436 '''
1438 def __init__(self, ztop, zbot, name=None):
1439 self.ztop = ztop
1440 self.zbot = zbot
1441 self.zmid = (self.ztop + self.zbot) * 0.5
1442 self.name = name
1443 self.ilayer = None
1445 def _update_potint_coefs(self):
1446 potint_p = potint_s = False
1447 try:
1448 self._ppic = potint_coefs(
1449 self.mbot.vp, self.mtop.vp,
1450 radius(self.zbot), radius(self.ztop))
1451 potint_p = True
1452 except BadPotIntCoefs:
1453 pass
1455 potint_s = False
1456 try:
1457 self._spic = potint_coefs(
1458 self.mbot.vs, self.mtop.vs,
1459 radius(self.zbot), radius(self.ztop))
1460 potint_s = True
1461 except BadPotIntCoefs:
1462 pass
1464 assert P == 1 and S == 2
1465 self._use_potential_interpolation = (None, potint_p, potint_s)
1467 def potint_coefs(self, mode):
1468 '''
1469 Get coefficients for potential interpolation.
1471 :param mode: mode of wave propagation, :py:const:`P` or :py:const:`S`
1472 :returns: coefficients ``(a, b)``
1473 '''
1475 if mode == P:
1476 return self._ppic
1477 else:
1478 return self._spic
1480 def contains(self, z):
1481 '''
1482 Tolerantly check if a given depth is within the layer
1483 (including boundaries).
1484 '''
1486 return self.ztop <= z <= self.zbot or \
1487 self.at_bottom(z) or self.at_top(z)
1489 def inner(self, z):
1490 '''
1491 Tolerantly check if a given depth is within the layer
1492 (not including boundaries).
1493 '''
1495 return self.ztop <= z <= self.zbot and not \
1496 self.at_bottom(z) and not \
1497 self.at_top(z)
1499 def at_bottom(self, z):
1500 '''
1501 Tolerantly check if given depth is at the bottom of the layer.
1502 '''
1504 return abs(self.zbot - z) < ZEPS
1506 def at_top(self, z):
1507 '''
1508 Tolerantly check if given depth is at the top of the layer.
1509 '''
1510 return abs(self.ztop - z) < ZEPS
1512 def pflat_top(self, p):
1513 '''
1514 Convert spherical ray parameter to local flat ray parameter for top of
1515 layer.
1516 '''
1517 return p / (earthradius-self.ztop)
1519 def pflat_bottom(self, p):
1520 '''
1521 Convert spherical ray parameter to local flat ray parameter for bottom
1522 of layer.
1523 '''
1524 return p / (earthradius-self.zbot)
1526 def pflat(self, p, z):
1527 '''
1528 Convert spherical ray parameter to local flat ray parameter for
1529 given depth.
1530 '''
1531 return p / (earthradius-z)
1533 def v_potint(self, mode, z):
1534 a, b = self.potint_coefs(mode)
1535 return a*(earthradius-z)**b
1537 def u_potint(self, mode, z):
1538 a, b = self.potint_coefs(mode)
1539 return 1./(a*(earthradius-z)**b)
1541 def xt_potint(self, p, mode, zpart=None):
1542 '''
1543 Get travel time and distance for for traversal with given mode and ray
1544 parameter.
1546 :param p: ray parameter (spherical)
1547 :param mode: mode of propagation (:py:const:`P` or :py:const:`S`)
1548 :param zpart: if given, tuple with two depths to restrict computation
1549 to a part of the layer
1551 This implementation uses analytic formulas valid for a spherical earth
1552 in the case where the velocity c within the layer is given by potential
1553 interpolation of the form
1555 c(z) = a*z^b
1556 '''
1557 utop, ubot = self.u_top_bottom(mode)
1558 a, b = self.potint_coefs(mode)
1559 ztop = self.ztop
1560 zbot = self.zbot
1561 if zpart is not None:
1562 utop = self.u(mode, zpart[0])
1563 ubot = self.u(mode, zpart[1])
1564 ztop, zbot = zpart
1565 utop = 1./(a*(earthradius-ztop)**b)
1566 ubot = 1./(a*(earthradius-zbot)**b)
1568 r1 = radius(zbot)
1569 r2 = radius(ztop)
1570 burger_eta1 = r1 * ubot
1571 burger_eta2 = r2 * utop
1572 if b != 1:
1573 def cpe(eta):
1574 return num.arccos(num.minimum(p/num.maximum(eta, p/2), 1.0))
1576 def sep(eta):
1577 return num.sqrt(num.maximum(eta**2 - p**2, 0.0))
1579 x = (cpe(burger_eta2)-cpe(burger_eta1))/(1.0-b)
1580 t = (sep(burger_eta2)-sep(burger_eta1))/(1.0-b)
1581 else:
1582 lr = math.log(r2/r1)
1583 sap = num.sqrt(1.0/a**2 - p**2)
1584 x = p/sap * lr
1585 t = 1./(a**2 * sap)
1587 x *= r2d
1589 return x, t
1591 def test(self, p, mode, z):
1592 '''
1593 Check if wave mode can exist for given ray parameter at given depth
1594 within the layer.
1595 '''
1596 return (self.u(mode, z)*radius(z) - p) > 0.
1598 def tests(self, p, mode):
1599 utop, ubot = self.u_top_bottom(mode)
1600 return (
1601 (utop * radius(self.ztop) - p) > 0.,
1602 (ubot * radius(self.zbot) - p) > 0.)
1604 def zturn_potint(self, p, mode):
1605 '''
1606 Get turning depth for given ray parameter and propagation mode.
1607 '''
1609 a, b = self.potint_coefs(mode)
1610 r = num.exp(num.log(a*p)/(1.0-b))
1611 return earthradius-r
1613 def propagate(self, p, mode, direction):
1614 '''
1615 Propagate ray through layer.
1617 :param p: ray parameter
1618 :param mode: propagation mode
1619 :param direction: in direction (:py:const:`UP` or :py:const:`DOWN`
1620 '''
1621 if direction == DOWN:
1622 zin, zout = self.ztop, self.zbot
1623 else:
1624 zin, zout = self.zbot, self.ztop
1626 if self.v(mode, zin) == 0.0 or not self.test(p, mode, zin):
1627 raise CannotPropagate(direction, self.ilayer)
1629 if not self.test(p, mode, zout):
1630 return -direction
1631 else:
1632 return direction
1634 def resize(self, depth_min=None, depth_max=None):
1635 '''
1636 Change layer thinkness and interpolate material if required.
1637 '''
1638 if depth_min:
1639 mtop = self.material(depth_min)
1641 if depth_max:
1642 mbot = self.material(depth_max)
1644 self.mtop = mtop if depth_min else self.mtop
1645 self.mbot = mbot if depth_max else self.mbot
1646 self.ztop = depth_min if depth_min else self.ztop
1647 self.zbot = depth_max if depth_max else self.zbot
1648 self.zmid = self.ztop + (self.zbot - self.ztop)/2.
1651class DoesNotTurn(CakeError):
1652 pass
1655def radius(z):
1656 return earthradius - z
1659class HomogeneousLayer(Layer):
1660 '''
1661 Representation of a homogeneous layer in a layered earth model.
1663 Base class: :py:class:`Layer`.
1664 '''
1666 def __init__(self, ztop, zbot, m, name=None):
1667 Layer.__init__(self, ztop, zbot, name=name)
1668 self.m = m
1669 self.mtop = m
1670 self.mbot = m
1671 self._update_potint_coefs()
1673 def copy(self, ztop=None, zbot=None):
1674 if ztop is None:
1675 ztop = self.ztop
1677 if zbot is None:
1678 zbot = self.zbot
1680 return HomogeneousLayer(ztop, zbot, self.m, name=self.name)
1682 def material(self, z):
1683 return self.m
1685 def u(self, mode, z=None):
1686 if self._use_potential_interpolation[mode] and z is not None:
1687 return self.u_potint(mode, z)
1689 if mode == P:
1690 return 1./self.m.vp
1691 if mode == S:
1692 return 1./self.m.vs
1694 def u_top_bottom(self, mode):
1695 u = self.u(mode)
1696 return u, u
1698 def v(self, mode, z=None):
1699 if self._use_potential_interpolation[mode] and z is not None:
1700 return self.v_potint(mode, z)
1702 if mode == P:
1703 v = self.m.vp
1704 if mode == S:
1705 v = self.m.vs
1707 if num.isscalar(z):
1708 return v
1709 else:
1710 return num.full(len(z), v)
1712 def v_top_bottom(self, mode):
1713 v = self.v(mode)
1714 return v, v
1716 def xt(self, p, mode, zpart=None):
1717 if self._use_potential_interpolation[mode]:
1718 return self.xt_potint(p, mode, zpart)
1720 u = self.u(mode)
1721 pflat = self.pflat_bottom(p)
1722 if zpart is None:
1723 dz = (self.zbot - self.ztop)
1724 else:
1725 dz = abs(zpart[1]-zpart[0])
1727 u = self.u(mode)
1728 eps = u*0.001
1729 denom = num.sqrt(u**2 - pflat**2) + eps
1731 x = r2d*pflat/(earthradius-self.zmid) * dz / denom
1732 t = u**2 * dz / denom
1733 return x, t
1735 def zturn(self, p, mode):
1736 if self._use_potential_interpolation[mode]:
1737 return self.zturn_potint(p, mode)
1739 raise DoesNotTurn()
1741 def split(self, z):
1742 upper = HomogeneousLayer(self.ztop, z, self.m, name=self.name)
1743 lower = HomogeneousLayer(z, self.zbot, self.m, name=self.name)
1744 upper.ilayer = self.ilayer
1745 lower.ilayer = self.ilayer
1746 return upper, lower
1748 def __str__(self):
1749 if self.name:
1750 name = self.name + ' '
1751 else:
1752 name = ''
1754 calcmode = ''.join('HP'[self._use_potential_interpolation[mode]]
1755 for mode in (P, S))
1757 return ' (%i) homogeneous layer %s(%g km - %g km) [%s]\n %s' % (
1758 self.ilayer, name, self.ztop/km, self.zbot/km, calcmode, self.m)
1761class GradientLayer(Layer):
1762 '''
1763 Representation of a gradient layer in a layered earth model.
1765 Base class: :py:class:`Layer`.
1766 '''
1768 def __init__(self, ztop, zbot, mtop, mbot, name=None):
1769 Layer.__init__(self, ztop, zbot, name=name)
1770 self.mtop = mtop
1771 self.mbot = mbot
1772 self._update_potint_coefs()
1774 def copy(self, ztop=None, zbot=None):
1775 if ztop is None:
1776 ztop = self.ztop
1778 if zbot is None:
1779 zbot = self.zbot
1781 return GradientLayer(ztop, zbot, self.mtop, self.mbot, name=self.name)
1783 def interpolate(self, z, ptop, pbot):
1784 return ptop + (z - self.ztop)*(pbot - ptop)/(self.zbot-self.ztop)
1786 def material(self, z):
1787 dtop = self.mtop.astuple()
1788 dbot = self.mbot.astuple()
1789 d = [
1790 self.interpolate(z, ptop, pbot)
1791 for (ptop, pbot) in zip(dtop, dbot)]
1793 return Material(*d)
1795 def u_top_bottom(self, mode):
1796 if mode == P:
1797 return 1./self.mtop.vp, 1./self.mbot.vp
1798 if mode == S:
1799 return 1./self.mtop.vs, 1./self.mbot.vs
1801 def u(self, mode, z):
1802 if self._use_potential_interpolation[mode]:
1803 return self.u_potint(mode, z)
1805 if mode == P:
1806 return 1./self.interpolate(z, self.mtop.vp, self.mbot.vp)
1807 if mode == S:
1808 return 1./self.interpolate(z, self.mtop.vs, self.mbot.vs)
1810 def v_top_bottom(self, mode):
1811 if mode == P:
1812 return self.mtop.vp, self.mbot.vp
1813 if mode == S:
1814 return self.mtop.vs, self.mbot.vs
1816 def v(self, mode, z):
1817 if self._use_potential_interpolation[mode]:
1818 return self.v_potint(mode, z)
1820 if mode == P:
1821 return self.interpolate(z, self.mtop.vp, self.mbot.vp)
1822 if mode == S:
1823 return self.interpolate(z, self.mtop.vs, self.mbot.vs)
1825 def xt(self, p, mode, zpart=None):
1826 if self._use_potential_interpolation[mode]:
1827 return self.xt_potint(p, mode, zpart)
1829 utop, ubot = self.u_top_bottom(mode)
1830 b = (1./ubot - 1./utop)/(self.zbot - self.ztop)
1832 pflat = self.pflat_bottom(p)
1833 if zpart is not None:
1834 utop = self.u(mode, zpart[0])
1835 ubot = self.u(mode, zpart[1])
1837 peps = 1e-16
1838 pdp = pflat + peps
1840 def func(u):
1841 eta = num.sqrt(num.maximum(u**2 - pflat**2, 0.0))
1842 xx = eta/u
1843 tt = num.where(
1844 pflat <= u,
1845 num.log(u+eta) - num.log(pdp) - eta/u,
1846 0.0)
1848 return xx, tt
1850 xxtop, tttop = func(utop)
1851 xxbot, ttbot = func(ubot)
1853 x = (xxtop - xxbot) / (b*pdp)
1854 t = (tttop - ttbot) / b + pflat*x
1856 x *= r2d/(earthradius - self.zmid)
1857 return x, t
1859 def zturn(self, p, mode):
1860 if self._use_potential_interpolation[mode]:
1861 return self.zturn_potint(p, mode)
1862 pflat = self.pflat_bottom(p)
1863 vtop, vbot = self.v_top_bottom(mode)
1864 return (1./pflat - vtop) * (self.zbot - self.ztop) / \
1865 (vbot-vtop) + self.ztop
1867 def split(self, z):
1868 mmid = self.material(z)
1869 upper = GradientLayer(self.ztop, z, self.mtop, mmid, name=self.name)
1870 lower = GradientLayer(z, self.zbot, mmid, self.mbot, name=self.name)
1871 upper.ilayer = self.ilayer
1872 lower.ilayer = self.ilayer
1873 return upper, lower
1875 def __str__(self):
1876 if self.name:
1877 name = self.name + ' '
1878 else:
1879 name = ''
1881 calcmode = ''.join('HP'[self._use_potential_interpolation[mode]]
1882 for mode in (P, S))
1884 return ''' (%i) gradient layer %s(%g km - %g km) [%s]
1885 %s
1886 %s''' % (
1887 self.ilayer,
1888 name,
1889 self.ztop/km,
1890 self.zbot/km,
1891 calcmode,
1892 self.mtop,
1893 self.mbot)
1896class Discontinuity(object):
1897 '''
1898 Base class for discontinuities in layered earth model.
1900 Subclasses are: :py:class:`Interface` and :py:class:`Surface`.
1901 '''
1903 def __init__(self, z, name=None):
1904 self.z = z
1905 self.zbot = z
1906 self.ztop = z
1907 self.name = name
1909 def change_depth(self, z):
1910 self.z = z
1911 self.zbot = z
1912 self.ztop = z
1914 def copy(self):
1915 return copy.deepcopy(self)
1918class Interface(Discontinuity):
1919 '''
1920 Representation of an interface in a layered earth model.
1922 Base class: :py:class:`Discontinuity`.
1923 '''
1925 def __init__(self, z, mabove, mbelow, name=None):
1926 Discontinuity.__init__(self, z, name)
1927 self.mabove = mabove
1928 self.mbelow = mbelow
1930 def __str__(self):
1931 if self.name is None:
1932 return 'interface'
1933 else:
1934 return 'interface "%s"' % self.name
1936 def u_top_bottom(self, mode):
1937 if mode == P:
1938 return reci_or_none(self.mabove.vp), reci_or_none(self.mbelow.vp)
1939 if mode == S:
1940 return reci_or_none(self.mabove.vs), reci_or_none(self.mbelow.vs)
1942 def critical_ps(self, mode):
1943 uabove, ubelow = self.u_top_bottom(mode)
1944 return (
1945 mult_or_none(uabove, radius(self.z)),
1946 mult_or_none(ubelow, radius(self.z)))
1948 def propagate(self, p, mode, direction):
1949 uabove, ubelow = self.u_top_bottom(mode)
1950 if direction == DOWN:
1951 if ubelow is not None and ubelow*radius(self.z) - p >= 0:
1952 return direction
1953 else:
1954 return -direction
1955 if direction == UP:
1956 if uabove is not None and uabove*radius(self.z) - p >= 0:
1957 return direction
1958 else:
1959 return -direction
1961 def pflat(self, p):
1962 return p / (earthradius-self.z)
1964 def efficiency(self, in_direction, out_direction, in_mode, out_mode, p):
1965 scatter = psv_solid(
1966 self.mabove, self.mbelow, self.pflat(p), energy=True)
1967 return scatter[
1968 psv_solid_ind(in_direction, out_direction, in_mode, out_mode)]
1971class Surface(Discontinuity):
1972 '''
1973 Representation of the surface discontinuity in a layered earth model.
1975 Base class: :py:class:`Discontinuity`.
1976 '''
1978 def __init__(self, z, mbelow):
1979 Discontinuity.__init__(self, z, 'surface')
1980 self.z = z
1981 self.mbelow = mbelow
1983 def propagate(self, p, mode, direction):
1984 return direction # no implicit reflection at surface
1986 def u_top_bottom(self, mode):
1987 if mode == P:
1988 return None, reci_or_none(self.mbelow.vp)
1989 if mode == S:
1990 return None, reci_or_none(self.mbelow.vs)
1992 def critical_ps(self, mode):
1993 _, ubelow = self.u_top_bottom(mode)
1994 return None, mult_or_none(ubelow, radius(self.z))
1996 def pflat(self, p):
1997 return p / (earthradius-self.z)
1999 def efficiency(self, in_direction, out_direction, in_mode, out_mode, p):
2000 if in_direction == DOWN or out_direction == UP:
2001 return 0.0
2002 else:
2003 return psv_surface(
2004 self.mbelow, self.pflat(p), energy=True)[
2005 psv_surface_ind(in_mode, out_mode)]
2007 def __str__(self):
2008 return 'surface'
2011class Walker(object):
2012 def __init__(self, elements):
2013 self._elements = elements
2014 self._i = 0
2016 def current(self):
2017 return self._elements[self._i]
2019 def go(self, direction):
2020 if direction == UP:
2021 self.up()
2022 else:
2023 self.down()
2025 def down(self):
2026 if self._i < len(self._elements)-1:
2027 self._i += 1
2028 else:
2029 raise BottomReached()
2031 def up(self):
2032 if self._i > 0:
2033 self._i -= 1
2034 else:
2035 raise SurfaceReached()
2037 def goto_layer(self, layer):
2038 self._i = self._elements.index(layer)
2041class RayElement(object):
2042 '''
2043 An element of a :py:class:`RayPath`.
2044 '''
2046 def __eq__(self, other):
2047 return type(self) is type(other) and self.__dict__ == other.__dict__
2049 def is_straight(self):
2050 return isinstance(self, Straight)
2052 def is_kink(self):
2053 return isinstance(self, Kink)
2056class Straight(RayElement):
2057 '''
2058 A ray segment representing wave propagation through one :py:class:`Layer`.
2059 '''
2061 def __init__(self, direction_in, direction_out, mode, layer):
2062 self.mode = mode
2063 self._direction_in = direction_in
2064 self._direction_out = direction_out
2065 self.layer = layer
2067 def angle_in(self, p, endgaps=None):
2068 z = self.z_in(endgaps)
2069 dir = self.eff_direction_in(endgaps)
2070 v = self.layer.v(self.mode, z)
2071 pf = self.layer.pflat(p, z)
2073 if dir == DOWN:
2074 return num.arcsin(v*pf)*r2d
2075 else:
2076 return 180.-num.arcsin(v*pf)*r2d
2078 def angle_out(self, p, endgaps=None):
2079 z = self.z_out(endgaps)
2080 dir = self.eff_direction_out(endgaps)
2081 v = self.layer.v(self.mode, z)
2082 pf = self.layer.pflat(p, z)
2084 if dir == DOWN:
2085 return 180.-num.arcsin(v*pf)*r2d
2086 else:
2087 return num.arcsin(v*pf)*r2d
2089 def pflat_in(self, p, endgaps=None):
2090 return p / (earthradius-self.z_in(endgaps))
2092 def pflat_out(self, p, endgaps=None):
2093 return p / (earthradius-self.z_out(endgaps))
2095 def test(self, p, z):
2096 return self.layer.test(p, self.mode, z)
2098 def z_in(self, endgaps=None):
2099 if endgaps is not None:
2100 return endgaps[0]
2101 else:
2102 lyr = self.layer
2103 return (lyr.ztop, lyr.zbot)[self._direction_in == UP]
2105 def z_out(self, endgaps=None):
2106 if endgaps is not None:
2107 return endgaps[1]
2108 else:
2109 lyr = self.layer
2110 return (lyr.ztop, lyr.zbot)[self._direction_out == DOWN]
2112 def turns(self):
2113 return self._direction_in != self._direction_out
2115 def eff_direction_in(self, endgaps=None):
2116 if endgaps is None:
2117 return self._direction_in
2118 else:
2119 return endgaps[2]
2121 def eff_direction_out(self, endgaps=None):
2122 if endgaps is None:
2123 return self._direction_out
2124 else:
2125 return endgaps[3]
2127 def zturn(self, p):
2128 lyr = self.layer
2129 return lyr.zturn(p, self.mode)
2131 def u_in(self, endgaps=None):
2132 return self.layer.u(self.mode, self.z_in(endgaps))
2134 def u_out(self, endgaps=None):
2135 return self.layer.u(self.mode, self.z_out(endgaps))
2137 def critical_p_in(self, endgaps=None):
2138 z = self.z_in(endgaps)
2139 return self.layer.u(self.mode, z)*radius(z)
2141 def critical_p_out(self, endgaps=None):
2142 z = self.z_out(endgaps)
2143 return self.layer.u(self.mode, z)*radius(z)
2145 def xt(self, p, zpart=None):
2146 x, t = self.layer.xt(p, self.mode, zpart=zpart)
2147 if self._direction_in != self._direction_out and zpart is None:
2148 x *= 2.
2149 t *= 2.
2150 return x, t
2152 def xt_gap(self, p, zstart, zstop, samedir):
2153 z1, z2 = zstart, zstop
2154 if z1 > z2:
2155 z1, z2 = z2, z1
2157 x, t = self.layer.xt(p, self.mode, zpart=(z1, z2))
2159 if samedir:
2160 return x, t
2161 else:
2162 xfull, tfull = self.xt(p)
2163 return xfull-x, tfull-t
2165 def __hash__(self):
2166 return hash((
2167 self._direction_in,
2168 self._direction_out,
2169 self.mode,
2170 id(self.layer)))
2173class HeadwaveStraight(Straight):
2174 def __init__(self, direction_in, direction_out, mode, interface):
2175 Straight.__init__(self, direction_in, direction_out, mode, None)
2177 self.interface = interface
2179 def z_in(self, zpart=None):
2180 return self.interface.z
2182 def z_out(self, zpart=None):
2183 return self.interface.z
2185 def zturn(self, p):
2186 return num.full(len(p), self.interface.z)
2188 def xt(self, p, zpart=None):
2189 return 0., 0.
2191 def x2t_headwave(self, xstretch):
2192 xstretch_m = xstretch*d2r*radius(self.interface.z)
2193 return min_not_none(*self.interface.u_top_bottom(self.mode))*xstretch_m
2196class Kink(RayElement):
2197 '''
2198 An interaction of a ray with a :py:class:`Discontinuity`.
2199 '''
2201 def __init__(
2202 self,
2203 in_direction,
2204 out_direction,
2205 in_mode,
2206 out_mode,
2207 discontinuity):
2209 self.in_direction = in_direction
2210 self.out_direction = out_direction
2211 self.in_mode = in_mode
2212 self.out_mode = out_mode
2213 self.discontinuity = discontinuity
2215 def reflection(self):
2216 return self.in_direction != self.out_direction
2218 def conversion(self):
2219 return self.in_mode != self.out_mode
2221 def efficiency(self, p, out_direction=None, out_mode=None):
2223 if out_direction is None:
2224 out_direction = self.out_direction
2226 if out_mode is None:
2227 out_mode = self.out_mode
2229 return self.discontinuity.efficiency(
2230 self.in_direction, out_direction, self.in_mode, out_mode, p)
2232 def __str__(self):
2233 r, c = self.reflection(), self.conversion()
2234 if r and c:
2235 return '|~'
2236 if r:
2237 return '|'
2238 if c:
2239 return '~'
2240 return '_'
2242 def __hash__(self):
2243 return hash((
2244 self.in_direction,
2245 self.out_direction,
2246 self.in_mode,
2247 self.out_mode,
2248 id(self.discontinuity)))
2251class PRangeNotSet(CakeError):
2252 pass
2255class RayPath(object):
2256 '''
2257 Representation of a fan of rays running through a common sequence of
2258 layers / interfaces.
2259 '''
2261 def __init__(self, phase):
2262 self.elements = []
2263 self.phase = phase
2264 self._pmax = None
2265 self._pmin = None
2266 self._p = None
2267 self._is_headwave = False
2269 def set_is_headwave(self, is_headwave):
2270 self._is_headwave = is_headwave
2272 def copy(self):
2273 '''
2274 Get a copy of it.
2275 '''
2277 c = copy.copy(self)
2278 c.elements = list(self.elements)
2279 return c
2281 def endgaps(self, zstart, zstop):
2282 '''
2283 Get information needed for end point adjustments.
2284 '''
2286 return (
2287 zstart,
2288 zstop,
2289 self.phase.direction_start(),
2290 self.phase.direction_stop())
2292 def append(self, element):
2293 self.elements.append(element)
2295 def _check_have_prange(self):
2296 if self._pmax is None:
2297 raise PRangeNotSet()
2299 def set_prange(self, pmin, pmax, dp):
2300 self._pmin, self._pmax = pmin, pmax
2301 self._prange_dp = dp
2303 def used_phase(self, p=None, eps=1.):
2304 '''
2305 Calculate phase definition from ray path.
2306 '''
2308 used = PhaseDef()
2309 fleg = self.phase.first_leg()
2310 used.append(Leg(fleg.departure, fleg.mode))
2311 n_elements_n = [None] + self.elements + [None]
2312 for before, element, after in zip(
2313 n_elements_n[:-2],
2314 n_elements_n[1:-1],
2315 n_elements_n[2:]):
2317 if element.is_kink() and HeadwaveStraight not in (
2318 type(before),
2319 type(after)):
2321 if element.reflection() or element.conversion():
2322 z = element.discontinuity.z
2323 used.append(Knee(
2324 z,
2325 element.in_direction,
2326 element.out_direction != element.in_direction,
2327 element.in_mode,
2328 element.out_mode))
2330 used.append(Leg(element.out_direction, element.out_mode))
2332 elif type(element) is HeadwaveStraight:
2333 z = element.interface.z
2334 k = Knee(
2335 z,
2336 before.in_direction,
2337 after.out_direction != before.in_direction,
2338 before.in_mode,
2339 after.out_mode)
2341 k.headwave = True
2342 used.append(k)
2343 used.append(Leg(after.out_direction, after.out_mode))
2345 if (p is not None and before and after
2346 and element.is_straight()
2347 and before.is_kink()
2348 and after.is_kink()
2349 and element.turns()
2350 and not before.reflection() and not before.conversion()
2351 and not after.reflection() and not after.conversion()):
2353 ai = element.angle_in(p)
2354 if 90.0-eps < ai and ai < 90+eps:
2355 used.append(
2356 Head(
2357 before.discontinuity.z,
2358 before.out_direction,
2359 element.mode))
2360 used.append(
2361 Leg(-before.out_direction, element.mode))
2363 used._direction_stop = self.phase.direction_stop()
2364 used._definition = self.phase.definition()
2366 return used
2368 def pmax(self):
2369 '''
2370 Get maximum valid ray parameter.
2371 '''
2372 self._check_have_prange()
2373 return self._pmax
2375 def pmin(self):
2376 '''
2377 Get minimum valid ray parameter.
2378 '''
2379 self._check_have_prange()
2380 return self._pmin
2382 def xmin(self):
2383 '''
2384 Get minimal distance.
2385 '''
2386 self._analyse()
2387 return self._xmin
2389 def xmax(self):
2390 '''
2391 Get maximal distance.
2392 '''
2393 self._analyse()
2394 return self._xmax
2396 def kinks(self):
2397 '''
2398 Iterate over propagation mode changes (reflections/transmissions).
2399 '''
2400 return (k for k in self.elements if isinstance(k, Kink))
2402 def straights(self):
2403 '''
2404 Iterate over ray segments.
2405 '''
2406 return (s for s in self.elements if isinstance(s, Straight))
2408 def headwave_straight(self):
2409 for s in self.elements:
2410 if type(s) is HeadwaveStraight:
2411 return s
2413 def first_straight(self):
2414 '''
2415 Get first ray segment.
2416 '''
2417 for s in self.elements:
2418 if isinstance(s, Straight):
2419 return s
2421 def last_straight(self):
2422 '''
2423 Get last ray segment.
2424 '''
2425 for s in reversed(self.elements):
2426 if isinstance(s, Straight):
2427 return s
2429 def efficiency(self, p):
2430 '''
2431 Get product of all conversion/reflection coefficients encountered on
2432 path.
2433 '''
2434 return reduce(
2435 operator.mul, (k.efficiency(p) for k in self.kinks()), 1.)
2437 def spreading(self, p, endgaps):
2438 '''
2439 Get geometrical spreading factor.
2440 '''
2441 if self._is_headwave:
2442 return 0.0
2444 self._check_have_prange()
2445 dp = self._prange_dp * 0.01
2446 assert self._pmax - self._pmin > dp
2448 if p + dp > self._pmax:
2449 p = p-dp
2451 x0, t = self.xt(p, endgaps)
2452 x1, t = self.xt(p+dp, endgaps)
2453 x0 *= d2r
2454 x1 *= d2r
2455 if x1 == x0:
2456 return num.nan
2458 dp_dx = dp/(x1-x0)
2460 x = x0
2461 if x == 0.:
2462 x = x1
2463 p = dp
2465 first = self.first_straight()
2466 last = self.last_straight()
2467 return num.abs(dp_dx) * first.pflat_in(p, endgaps) / (
2468 4.0 * math.pi * num.sin(x) *
2469 (earthradius-first.z_in(endgaps)) *
2470 (earthradius-last.z_out(endgaps))**2 *
2471 first.u_in(endgaps)**2 *
2472 num.abs(num.cos(first.angle_in(p, endgaps)*d2r)) *
2473 num.abs(num.cos(last.angle_out(p, endgaps)*d2r)))
2475 def make_p(self, dp=None, n=None, nmin=None):
2476 assert dp is None or n is None
2478 if self._pmin == self._pmax:
2479 return num.array([self._pmin])
2481 if dp is None:
2482 dp = self._prange_dp
2484 if n is None:
2485 n = int(round((self._pmax-self._pmin)/dp)) + 1
2487 if nmin is not None:
2488 n = max(n, nmin)
2490 ppp = num.linspace(self._pmin, self._pmax, n)
2491 return ppp
2493 def xt_endgaps(self, p, endgaps, which='both'):
2494 '''
2495 Get amount of distance/traveltime to be subtracted at the generic ray
2496 path's ends.
2497 '''
2499 zstart, zstop, dirstart, dirstop = endgaps
2500 firsts = self.first_straight()
2501 lasts = self.last_straight()
2502 xs, ts = firsts.xt_gap(
2503 p, zstart, firsts.z_in(), dirstart == firsts._direction_in)
2504 xe, te = lasts.xt_gap(
2505 p, zstop, lasts.z_out(), dirstop == lasts._direction_out)
2507 if which == 'both':
2508 return xs + xe, ts + te
2509 elif which == 'left':
2510 return xs, ts
2511 elif which == 'right':
2512 return xe, te
2514 def xt_endgaps_ptest(self, p, endgaps):
2515 '''
2516 Check if ray parameter is valid at source and receiver.
2517 '''
2519 zstart, zstop, dirstart, dirstop = endgaps
2520 firsts = self.first_straight()
2521 lasts = self.last_straight()
2522 return num.logical_and(firsts.test(p, zstart), lasts.test(p, zstop))
2524 def xt(self, p, endgaps):
2525 '''
2526 Calculate distance and traveltime for given ray parameter.
2527 '''
2529 if isinstance(p, num.ndarray):
2530 sx = num.zeros(p.size)
2531 st = num.zeros(p.size)
2532 else:
2533 sx = 0.0
2534 st = 0.0
2536 for s in self.straights():
2537 x, t = s.xt(p)
2538 sx += x
2539 st += t
2541 if endgaps:
2542 dx, dt = self.xt_endgaps(p, endgaps)
2543 sx -= dx
2544 st -= dt
2546 return sx, st
2548 def xt_limits(self, p):
2549 '''
2550 Calculate limits of distance and traveltime for given ray parameter.
2551 '''
2553 if isinstance(p, num.ndarray):
2554 sx = num.zeros(p.size)
2555 st = num.zeros(p.size)
2556 sxe = num.zeros(p.size)
2557 ste = num.zeros(p.size)
2558 else:
2559 sx = 0.0
2560 st = 0.0
2561 sxe = 0.0
2562 ste = 0.0
2564 sfirst = self.first_straight()
2565 slast = self.last_straight()
2567 for s in self.straights():
2568 if s is not sfirst and s is not slast:
2569 x, t = s.xt(p)
2570 sx += x
2571 st += t
2573 sends = [sfirst]
2574 if sfirst is not slast:
2575 sends.append(slast)
2577 for s in sends:
2578 x, t = s.xt(p)
2579 sxe += x
2580 ste += t
2582 return sx, (sx + sxe), st, (st + ste)
2584 def iter_zxt(self, p):
2585 '''
2586 Iterate over (depth, distance, traveltime) at each layer interface on
2587 ray path.
2588 '''
2590 sx = num.zeros(p.size)
2591 st = num.zeros(p.size)
2592 ok = False
2593 for s in self.straights():
2594 yield s.z_in(), sx.copy(), st.copy()
2596 x, t = s.xt(p)
2597 sx += x
2598 st += t
2599 ok = True
2601 if ok:
2602 yield s.z_out(), sx.copy(), st.copy()
2604 def zxt_path_subdivided(
2605 self, p, endgaps,
2606 points_per_straight=20,
2607 x_for_headwave=None,
2608 annotated=False):
2610 '''
2611 Get geometrical representation of ray path.
2612 '''
2614 if self._is_headwave:
2615 assert p.size == 1
2616 x, t = self.xt(p, endgaps)
2617 xstretch = x_for_headwave-x
2618 nout = xstretch.size
2619 else:
2620 nout = p.size
2622 dxl, dtl = self.xt_endgaps(p, endgaps, which='left')
2623 dxr, dtr = self.xt_endgaps(p, endgaps, which='right')
2625 # first create full path including the endgaps
2626 sx = num.zeros(nout) - dxl
2627 st = num.zeros(nout) - dtl
2628 zxt = []
2629 for s in self.straights():
2630 n = points_per_straight
2632 back = None
2633 zin, zout = s.z_in(), s.z_out()
2634 if type(s) is HeadwaveStraight:
2635 z = zin
2636 for i in range(n):
2637 xs = float(i)/(n-1) * xstretch
2638 ts = s.x2t_headwave(xs)
2639 zs = num.full(xstretch.size, z)
2640 zxt.append((
2641 zs,
2642 sx+xs,
2643 st+ts,
2644 num.full(xstretch.size, s.mode),
2645 num.full(xstretch.size, 2),
2646 s.layer.v(s.mode, zs)))
2647 else:
2648 if zin != zout: # normal traversal
2649 zs = num.linspace(zin, zout, n)
2650 vs = s.layer.v(s.mode, zs)
2651 for z, v in zip(zs, vs):
2652 x, t = s.xt(p, zpart=sorted([zin, z]))
2653 zxt.append((
2654 num.full(nout, z),
2655 sx + x,
2656 st + t,
2657 num.full(nout, s.mode),
2658 num.full(nout, 0),
2659 num.full(nout, v)))
2661 else: # ray turns in layer
2662 zturns = s.zturn(p)
2663 back = []
2664 for i in range(n):
2665 zs = zin + (zturns - zin) * num.sin(
2666 float(i)/(n-1)*math.pi/2.0) * 0.999
2668 vs = s.layer.v(s.mode, zs)
2670 if zturns[0] >= zin:
2671 x, t = s.xt(p, zpart=[zin, zs])
2672 else:
2673 x, t = s.xt(p, zpart=[zs, zin])
2675 zxt.append((
2676 zs, sx + x, st + t,
2677 num.full(zs.size, s.mode),
2678 num.full(zs.size, 1),
2679 vs))
2681 back.append((
2682 zs, x, t,
2683 num.full(zs.size, s.mode),
2684 num.full(zs.size, 1),
2685 vs))
2687 if type(s) is HeadwaveStraight:
2688 x = xstretch
2689 t = s.x2t_headwave(xstretch)
2690 else:
2691 x, t = s.xt(p)
2693 sx += x
2694 st += t
2695 if back:
2696 for z, x, t, mode, kind, v in reversed(back):
2697 zxt.append((z, sx - x, st - t, mode, kind, v))
2699 # gather results as arrays with such that x[ip, ipoint]
2700 fan_z, fan_x, fan_t, fan_mode, fan_kind, fan_v = [], [], [], [], [], []
2701 for z, x, t, mode, kind, v in zxt:
2702 fan_z.append(z)
2703 fan_x.append(x)
2704 fan_t.append(t)
2705 fan_mode.append(mode)
2706 fan_kind.append(kind)
2707 fan_v.append(v)
2709 z = num.array(fan_z).T
2710 x = num.array(fan_x).T
2711 t = num.array(fan_t).T
2712 mode = num.array(fan_mode).T
2713 kind = num.array(fan_kind).T
2714 v = num.array(fan_v).T
2716 # cut off the endgaps, add exact endpoints
2717 xmax = x[:, -1] - dxr
2718 tmax = t[:, -1] - dtr
2719 zstart, zstop = endgaps[:2]
2720 zs, xs, ts, modes, kinds, vs = [], [], [], [], [], []
2721 for i in range(nout):
2722 t_ = t[i]
2723 indices = num.where(num.logical_and(0. <= t_, t_ <= tmax[i]))[0]
2724 n = indices.size + 2
2725 zs_, xs_, ts_, modes_, kinds_, vs_ = [
2726 num.empty(n, dtype=float) for j in range(6)]
2727 zs_[1:-1] = z[i, indices]
2728 xs_[1:-1] = x[i, indices]
2729 ts_[1:-1] = t[i, indices]
2730 modes_[1:-1] = mode[i, indices]
2731 kinds_[1:-1] = kind[i, indices]
2732 vs_[1:-1] = v[i, indices]
2733 zs_[0], zs_[-1] = zstart, zstop
2734 xs_[0], xs_[-1] = 0., xmax[i]
2735 ts_[0], ts_[-1] = 0., tmax[i]
2736 modes_[0] = modes_[1]
2737 modes_[-1] = modes_[-2]
2738 kinds_[0] = kinds_[1]
2739 kinds_[-1] = kinds_[-2]
2740 vs_[0] = vs_[1]
2741 vs_[-1] = vs_[-2]
2742 zs.append(zs_)
2743 xs.append(xs_)
2744 ts.append(ts_)
2745 modes.append(modes_)
2746 kinds.append(kinds_)
2747 vs.append(vs_)
2749 if annotated:
2750 return zs, xs, ts, modes, kinds, vs
2751 else:
2752 return zs, xs, ts
2754 def _analyse(self):
2755 if self._p is not None:
2756 return
2758 p = self.make_p(nmin=20)
2759 xmin, xmax, tmin, tmax = self.xt_limits(p)
2761 self._x, self._t, self._p = xmax, tmax, p
2762 self._xmin, self._xmax = xmin.min(), xmax.max()
2763 self._tmin, self._tmax = tmin.min(), tmax.max()
2765 def draft_pxt(self, endgaps):
2766 self._analyse()
2768 if not self._is_headwave:
2769 cp, cx, ct = self._p, self._x, self._t
2770 pcrit = min(
2771 self.critical_pstart(endgaps),
2772 self.critical_pstop(endgaps))
2774 if pcrit < self._pmin:
2775 empty = num.array([], dtype=float)
2776 return empty, empty, empty
2778 elif pcrit >= self._pmax:
2779 dx, dt = self.xt_endgaps(cp, endgaps)
2780 return cp, cx-dx, ct-dt
2782 else:
2783 n = num.searchsorted(cp, pcrit) + 1
2784 rp, rx, rt = num.empty((3, n), dtype=float)
2785 rp[:-1] = cp[:n-1]
2786 rx[:-1] = cx[:n-1]
2787 rt[:-1] = ct[:n-1]
2788 rp[-1] = pcrit
2789 rx[-1], rt[-1] = self.xt(pcrit, endgaps)
2790 dx, dt = self.xt_endgaps(rp, endgaps)
2791 rx[:-1] -= dx[:-1]
2792 rt[:-1] -= dt[:-1]
2793 return rp, rx, rt
2795 else:
2796 dx, dt = self.xt_endgaps(self._p, endgaps)
2797 p, x, t = self._p, self._x - dx, self._t - dt
2798 p, x, t = p[0], x[0], t[0]
2799 xh = num.linspace(0., x*10-x, 10)
2800 th = self.headwave_straight().x2t_headwave(xh)
2801 return num.full(xh.size, p), x+xh, t+th
2803 def interpolate_x2pt_linear(self, x, endgaps):
2804 '''
2805 Get approximate ray parameter and traveltime for distance.
2806 '''
2808 self._analyse()
2810 if self._is_headwave:
2811 dx, dt = self.xt_endgaps(self._p, endgaps)
2812 xmin = self._x[0] - dx[0]
2813 tmin = self._t[0] - dt[0]
2814 el = self.headwave_straight()
2815 xok = x[x >= xmin]
2816 th = el.x2t_headwave(xstretch=(xok-xmin)) + tmin
2817 return [
2818 (x_, self._p[0], t, None) for (x_, t) in zip(xok, th)]
2820 else:
2821 if num.all(x < self._xmin) or num.all(self._xmax < x):
2822 return []
2824 rp, rx, rt = self.draft_pxt(endgaps)
2826 xp = interp(x, rx, rp, 0)
2827 xt = interp(x, rx, rt, 0)
2829 if (rp.size and
2830 len(xp) == 0 and
2831 rx[0] == 0.0 and
2832 any(x == 0.0) and
2833 rp[0] == 0.0):
2835 xp = [(0.0, rp[0])]
2836 xt = [(0.0, rt[0])]
2838 return [
2839 (x_, p, t, (rp, rx, rt)) for ((x_, p), (_, t)) in zip(xp, xt)]
2841 def __eq__(self, other):
2842 if len(self.elements) != len(other.elements):
2843 return False
2845 return all(a == b for a, b in zip(self.elements, other.elements))
2847 def __hash__(self):
2848 return hash(
2849 tuple(hash(x) for x in self.elements) +
2850 (self.phase.definition(), ))
2852 def __str__(self, p=None, eps=1.):
2853 x = []
2854 start_i = None
2855 end_i = None
2856 turn_i = None
2858 def append_layers(si, ei, ti):
2859 if si == ei and (ti is None or ti == si):
2860 x.append('%i' % si)
2861 else:
2862 if ti is not None:
2863 x.append('(%i-%i-%i)' % (si, ti, ei))
2864 else:
2865 x.append('(%i-%i)' % (si, ei))
2867 for el in self.elements:
2868 if type(el) is Straight:
2869 if start_i is None:
2870 start_i = el.layer.ilayer
2871 if el._direction_in != el._direction_out:
2872 turn_i = el.layer.ilayer
2873 end_i = el.layer.ilayer
2875 elif isinstance(el, Kink):
2876 if start_i is not None:
2877 append_layers(start_i, end_i, turn_i)
2878 start_i = None
2879 turn_i = None
2881 x.append(str(el))
2883 if start_i is not None:
2884 append_layers(start_i, end_i, turn_i)
2886 su = '(%s)' % self.used_phase(p=p, eps=eps).used_repr()
2888 return '%-15s %-17s %s' % (self.phase.definition(), su, ''.join(x))
2890 def critical_pstart(self, endgaps):
2891 '''
2892 Get critical ray parameter for source depth choice.
2893 '''
2895 return self.first_straight().critical_p_in(endgaps)
2897 def critical_pstop(self, endgaps):
2898 '''
2899 Get critical ray parameter for receiver depth choice.
2900 '''
2902 return self.last_straight().critical_p_out(endgaps)
2904 def ranges(self, endgaps):
2905 '''
2906 Get valid ranges of ray parameter, distance, and traveltime.
2907 '''
2908 p, x, t = self.draft_pxt(endgaps)
2909 return p.min(), p.max(), x.min(), x.max(), t.min(), t.max()
2911 def describe(self, endgaps=None, as_degrees=False):
2912 '''
2913 Get textual representation.
2914 '''
2916 self._analyse()
2918 if as_degrees:
2919 xunit = 'deg'
2920 xfact = 1.
2921 else:
2922 xunit = 'km'
2923 xfact = d2m/km
2925 sg = ''' Ranges for all depths in source and receiver layers:
2926 - x [%g, %g] %s
2927 - t [%g, %g] s
2928 - p [%g, %g] s/deg
2929''' % (
2930 self._xmin*xfact,
2931 self._xmax*xfact,
2932 xunit,
2933 self._tmin,
2934 self._tmax,
2935 self._pmin/r2d,
2936 self._pmax/r2d)
2938 if endgaps is not None:
2939 pmin, pmax, xmin, xmax, tmin, tmax = self.ranges(endgaps)
2940 ss = ''' Ranges for given source and receiver depths:
2941\n - x [%g, %g] %s
2942\n - t [%g, %g] s
2943\n - p [%g, %g] s/deg
2944\n''' % (xmin*xfact, xmax*xfact, xunit, tmin, tmax, pmin/r2d, pmax/r2d)
2946 else:
2947 ss = ''
2949 return '%s\n' % self + ss + sg
2952class RefineFailed(CakeError):
2953 pass
2956class Ray(object):
2957 '''
2958 Representation of a ray with a specific (path, ray parameter, distance,
2959 arrival time) choice.
2961 **Attributes:**
2963 .. py:attribute:: path
2965 :py:class:`RayPath` object containing complete propagation history.
2967 .. py:attribute:: p
2969 Ray parameter (spherical) [s/rad]
2971 .. py:attribute:: x
2973 Radial distance [deg]
2975 .. py:attribute:: t
2977 Traveltime [s]
2979 .. py:attribute:: endgaps
2981 Needed for source/receiver depth adjustments in many
2982 :py:class:`RayPath` methods.
2983 '''
2985 def __init__(self, path, p, x, t, endgaps, draft_pxt):
2986 self.path = path
2987 self.p = p
2988 self.x = x
2989 self.t = t
2990 self.endgaps = endgaps
2991 self.draft_pxt = draft_pxt
2993 def given_phase(self):
2994 '''
2995 Get phase definition which was used to create the ray.
2997 :returns: :py:class:`PhaseDef` object
2998 '''
3000 return self.path.phase
3002 def used_phase(self):
3003 '''
3004 Compute phase definition from propagation path.
3006 :returns: :py:class:`PhaseDef` object
3007 '''
3009 return self.path.used_phase(self.p)
3011 def refine(self):
3012 if self.path._is_headwave:
3013 return
3015 if self.t == 0.0 and self.p == 0.0 and self.x == 0.0:
3016 return
3018 cp, cx, ct = self.draft_pxt
3019 ip = num.searchsorted(cp, self.p)
3020 if not (0 < ip < cp.size):
3021 raise RefineFailed()
3023 pl, ph = cp[ip-1], cp[ip]
3024 p_to_t = {}
3025 i = [0]
3027 def f(p):
3028 i[0] += 1
3029 x, t = self.path.xt(p, self.endgaps)
3030 p_to_t[p] = t
3031 return self.x - x
3033 try:
3034 self.p = brentq(f, pl, ph)
3035 self.t = p_to_t[self.p]
3037 except ValueError:
3038 raise RefineFailed()
3040 def t_and_attributes(self, attributes):
3041 d = {
3042 'takeoff_angle': lambda: self.takeoff_angle(),
3043 'incidence_angle': lambda: self.incidence_angle(),
3044 't': lambda: self.t,
3045 'p': lambda: self.p}
3047 return tuple(d[k]() for k in ['t'] + attributes)
3049 def takeoff_angle(self):
3050 '''
3051 Get takeoff angle of ray.
3053 The angle is returned in [degrees].
3054 '''
3056 return self.path.first_straight().angle_in(self.p, self.endgaps)
3058 def incidence_angle(self):
3059 '''
3060 Get incidence angle of ray.
3062 The angle is returned in [degrees].
3063 '''
3065 return self.path.last_straight().angle_out(self.p, self.endgaps)
3067 def efficiency(self):
3068 '''
3069 Get conversion/reflection efficiency of the ray.
3071 A value between 0 and 1 is returned, reflecting the relative amount of
3072 energy which is transmitted along the ray and not lost by reflections
3073 or conversions.
3074 '''
3076 return self.path.efficiency(self.p)
3078 def spreading(self):
3079 '''
3080 Get geometrical spreading factor.
3081 '''
3083 return self.path.spreading(self.p, self.endgaps)
3085 def surface_sphere(self):
3086 x1, y1 = 0., earthradius - self.endgaps[0]
3087 r2 = earthradius - self.endgaps[1]
3088 x2, y2 = r2*math.sin(self.x*d2r), r2*math.cos(self.x*d2r)
3089 return ((x2-x1)**2 + (y2-y1)**2)*4.0*math.pi
3091 def zxt_path_subdivided(self, points_per_straight=20, annotated=False):
3092 '''
3093 Get geometrical representation of ray path.
3095 Three arrays (depth, distance, time) with points on the ray's path of
3096 propagation are returned. The number of points which are used in each
3097 ray segment (passage through one layer) may be controlled by the
3098 ``points_per_straight`` parameter.
3099 '''
3100 return self.path.zxt_path_subdivided(
3101 num.atleast_1d(self.p), self.endgaps,
3102 points_per_straight=points_per_straight,
3103 x_for_headwave=num.atleast_1d(self.x),
3104 annotated=annotated)
3106 def __str__(self, as_degrees=False):
3107 if as_degrees:
3108 sd = '%6.3g deg' % self.x
3109 else:
3110 sd = '%7.5g km' % (self.x*(d2r*earthradius/km))
3112 return '%7.5g s/deg %s %6.4g s %5.1f %5.1f %3.0f%% %3.0f%% %s' % (
3113 self.p/r2d,
3114 sd,
3115 self.t,
3116 self.takeoff_angle(),
3117 self.incidence_angle(),
3118 100*self.efficiency(),
3119 100*self.spreading()*self.surface_sphere(),
3120 self.path.__str__(p=self.p))
3123def anything_to_crust2_profile(crust2_profile):
3124 from pyrocko.dataset import crust2x2
3125 if isinstance(crust2_profile, tuple):
3126 lat, lon = [float(x) for x in crust2_profile]
3127 return crust2x2.get_profile(lat, lon)
3128 elif isinstance(crust2_profile, str):
3129 return crust2x2.get_profile(crust2_profile)
3130 elif isinstance(crust2_profile, crust2x2.Crust2Profile):
3131 return crust2_profile
3132 else:
3133 assert False, 'crust2_profile must be (lat, lon) a profile ' \
3134 'key or a crust2x2 Profile object)'
3137class DiscontinuityNotFound(CakeError):
3138 def __init__(self, depth_or_name):
3139 CakeError.__init__(self)
3140 self.depth_or_name = depth_or_name
3142 def __str__(self):
3143 return 'Cannot find discontinuity from given depth or name: %s' % \
3144 self.depth_or_name
3147class LayeredModelError(CakeError):
3148 pass
3151class LayeredModel(object):
3152 '''
3153 Representation of a layer cake model.
3155 There are several ways to initialize an instance of this class.
3157 1. Use the module function :py:func:`load_model` to read a model from a
3158 file.
3159 2. Create an empty model with the default constructor and append layers and
3160 discontinuities with the :py:meth:`append` method (from top to bottom).
3161 3. Use the constructor :py:meth:`LayeredModel.from_scanlines`, to
3162 automatically create the :py:class:`Layer` and :py:class:`Discontinuity`
3163 objects from a given velocity profile.
3165 An earth model is represented by as stack of :py:class:`Layer` and
3166 :py:class:`Discontinuity` objects. The method :py:meth:`arrivals` returns
3167 :py:class:`Ray` objects which may be e.g. queried for arrival times of
3168 specific phases. Each ray is associated with a :py:class:`RayPath` object.
3169 Ray objects share common ray paths if they have the same
3170 conversion/reflection/propagation history. Creating the ray path objects is
3171 relatively expensive (this is done in :py:meth:`gather_paths`), but they
3172 are cached for reuse in successive invocations.
3173 '''
3175 def __init__(self):
3176 self._surface_material = None
3177 self._elements = []
3178 self.nlayers = 0
3179 self._np = 10000
3180 self._pdepth = 5
3181 self._pathcache = {}
3183 def copy_with_elevation(self, elevation):
3184 '''
3185 Get copy of the model with surface layer stretched to given elevation.
3187 :param elevation: new surface elevation in [m]
3189 Elevation is positiv upward, contrary to the layered models downward
3190 `z` axis.
3191 '''
3193 c = copy.deepcopy(self)
3194 c._pathcache = {}
3195 surface = c._elements[0]
3196 toplayer = c._elements[1]
3198 assert toplayer.zbot > -elevation
3200 surface.z = -elevation
3201 c._elements[1] = toplayer.copy(ztop=-elevation)
3202 c._elements[1].ilayer = 0
3203 return c
3205 def zeq(self, z1, z2):
3206 return abs(z1-z2) < ZEPS
3208 def append(self, element):
3209 '''
3210 Add a layer or discontinuity at bottom of model.
3212 :param element: object of subclass of :py:class:`Layer` or
3213 :py:class:`Discontinuity`.
3214 '''
3216 if isinstance(element, Layer):
3217 if element.zbot >= earthradius:
3218 element.zbot = earthradius - 1.
3220 if element.ztop >= earthradius:
3221 raise CakeError('Layer deeper than earthradius')
3223 element.ilayer = self.nlayers
3224 self.nlayers += 1
3226 self._elements.append(element)
3228 def elements(self, direction=DOWN):
3229 '''
3230 Iterate over all elements of the model.
3232 :param direction: direction of traversal :py:const:`DOWN` or
3233 :py:const:`UP`.
3235 Objects derived from the :py:class:`Discontinuity` and
3236 :py:class:`Layer` classes are yielded.
3237 '''
3239 if direction == DOWN:
3240 return iter(self._elements)
3241 else:
3242 return reversed(self._elements)
3244 def layers(self, direction=DOWN):
3245 '''
3246 Iterate over all layers of model.
3248 :param direction: direction of traversal :py:const:`DOWN` or
3249 :py:const:`UP`.
3251 Objects derived from the :py:class:`Layer` class are yielded.
3252 '''
3254 if direction == DOWN:
3255 return (el for el in self._elements if isinstance(el, Layer))
3256 else:
3257 return (
3258 el for el in reversed(self._elements) if isinstance(el, Layer))
3260 def layer(self, z, direction=DOWN):
3261 '''
3262 Get layer for given depth.
3264 :param z: depth [m]
3265 :param direction: direction of traversal :py:const:`DOWN` or
3266 :py:const:`UP`.
3268 Returns first layer which touches depth ``z`` (tolerant at boundaries).
3269 '''
3271 for layer in self.layers(direction):
3272 if layer.contains(z):
3273 return layer
3274 else:
3275 raise CakeError('Failed extracting layer at depth z=%s' % z)
3277 def walker(self):
3278 return Walker(self._elements)
3280 def material(self, z, direction=DOWN):
3281 '''
3282 Get material at given depth.
3284 :param z: depth [m]
3285 :param direction: direction of traversal :py:const:`DOWN` or
3286 :py:const:`UP`
3287 :returns: object of type :py:class:`Material`
3289 If given depth ``z`` happens to be at an interface, the material of the
3290 first layer with respect to the the traversal ordering is returned.
3291 '''
3293 lyr = self.layer(z, direction)
3294 return lyr.material(z)
3296 def discontinuities(self):
3297 '''
3298 Iterate over all discontinuities of the model.'''
3300 return (el for el in self._elements if isinstance(el, Discontinuity))
3302 def discontinuity(self, name_or_z):
3303 '''
3304 Get discontinuity by name or depth.
3306 :param name_or_z: name of discontinuity or depth [m] as float value
3307 '''
3309 if isinstance(name_or_z, float):
3310 candi = sorted(
3311 self.discontinuities(), key=lambda i: abs(i.z-name_or_z))
3312 else:
3313 candi = [i for i in self.discontinuities() if i.name == name_or_z]
3315 if not candi:
3316 raise DiscontinuityNotFound(name_or_z)
3318 return candi[0]
3320 def adapt_phase(self, phase):
3321 '''
3322 Adapt a phase definition for use with this model.
3324 This returns a copy of the phase definition, where named
3325 discontinuities are replaced with the actual depth of these, as defined
3326 in the model.
3327 '''
3329 phase = phase.copy()
3330 for knee in phase.knees():
3331 if knee.depth != 'surface':
3332 knee.depth = self.discontinuity(knee.depth).z
3333 for leg in phase.legs():
3334 if leg.depthmax is not None and isinstance(leg.depthmax, str):
3335 leg.depthmax = self.discontinuity(leg.depthmax).z
3337 return phase
3339 def path(self, p, phase, layer_start, layer_stop):
3340 '''
3341 Get ray path for given combination of ray parameter, phase definition,
3342 source and receiver layers.
3344 :param p: ray parameter (spherical) [s/rad]
3345 :param phase: phase definition (:py:class:`PhaseDef` object)
3346 :param layer_start: layer with source
3347 :param layer_stop: layer with receiver
3348 :returns: :py:class:`RayPath` object
3350 If it is not possible to find a solution, an exception of type
3351 :py:exc:`NotPhaseConform`, :py:exc:`MinDepthReached`,
3352 :py:exc:`MaxDepthReached`, :py:exc:`CannotPropagate`,
3353 :py:exc:`BottomReached` or :py:exc:`SurfaceReached` is raised.
3354 '''
3356 phase = self.adapt_phase(phase)
3357 knees = phase.knees()
3358 legs = phase.legs()
3359 next_knee = next_or_none(knees)
3360 leg = next_or_none(legs)
3361 assert leg is not None
3363 direction = leg.departure
3364 direction_stop = phase.direction_stop()
3365 mode = leg.mode
3366 mode_stop = phase.last_leg().mode
3368 walker = self.walker()
3369 walker.goto_layer(layer_start)
3370 current = walker.current()
3372 ttop, tbot = current.tests(p, mode)
3373 if not ttop and not tbot:
3374 raise CannotPropagate(direction, current.ilayer)
3376 if (direction == DOWN and not ttop) or (direction == UP and not tbot):
3377 direction = -direction
3379 path = RayPath(phase)
3380 trapdetect = set()
3381 while True:
3382 at_layer = isinstance(current, Layer)
3383 at_discontinuity = isinstance(current, Discontinuity)
3385 # detect trapped wave
3386 k = (id(next_knee), id(current), direction, mode)
3387 if k in trapdetect:
3388 raise Trapped()
3390 trapdetect.add(k)
3392 if at_discontinuity:
3393 oldmode, olddirection = mode, direction
3394 headwave = False
3395 if next_knee is not None and next_knee.matches(
3396 current, mode, direction):
3398 headwave = next_knee.headwave
3399 direction = next_knee.out_direction()
3400 mode = next_knee.out_mode
3401 next_knee = next_or_none(knees)
3402 leg = next(legs)
3404 else: # implicit reflection/transmission
3405 direction = current.propagate(p, mode, direction)
3407 if headwave:
3408 path.set_is_headwave(True)
3410 path.append(Kink(
3411 olddirection, olddirection, oldmode, oldmode, current))
3413 path.append(HeadwaveStraight(
3414 olddirection, direction, oldmode, current))
3416 path.append(Kink(
3417 olddirection, direction, oldmode, mode, current))
3419 else:
3420 path.append(Kink(
3421 olddirection, direction, oldmode, mode, current))
3423 if at_layer:
3424 direction_in = direction
3425 direction = current.propagate(p, mode, direction_in)
3427 zturn = None
3428 if direction_in != direction:
3429 zturn = current.zturn(p, mode)
3431 zmin, zmax = leg.depthmin, leg.depthmax
3432 if zmin is not None or zmax is not None:
3433 if direction_in != direction:
3434 if zmin is not None and zturn <= zmin:
3435 raise MinDepthReached()
3436 if zmax is not None and zturn >= zmax:
3437 raise MaxDepthReached()
3438 else:
3439 if zmin is not None and current.ztop <= zmin:
3440 raise MinDepthReached()
3441 if zmax is not None and current.zbot >= zmax:
3442 raise MaxDepthReached()
3444 path.append(Straight(direction_in, direction, mode, current))
3446 if next_knee is None and mode == mode_stop and \
3447 current is layer_stop:
3449 if zturn is None:
3450 if direction == direction_stop:
3451 break
3452 else:
3453 break
3455 walker.go(direction)
3456 current = walker.current()
3458 return path
3460 def gather_paths(self, phases=PhaseDef('P'), zstart=0.0, zstop=0.0):
3461 '''
3462 Get all possible ray paths for given source and receiver depths for one
3463 or more phase definitions.
3465 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3466 Comma-separated strings and lists of such strings are also accepted
3467 and are converted to :py:class:`PhaseDef` objects for convenience.
3468 :param zstart: source depth [m]
3469 :param zstop: receiver depth [m]
3470 :returns: a list of :py:class:`RayPath` objects
3472 Results of this method are cached internally. Cached results are
3473 returned, when a given combination of source layer, receiver layer and
3474 phase definition has been used before.
3475 '''
3477 eps = 1e-7 # num.finfo(float).eps * 1000.
3479 phases = to_phase_defs(phases)
3481 paths = []
3482 for phase in phases:
3484 layer_start = self.layer(zstart, -phase.direction_start())
3485 layer_stop = self.layer(zstop, phase.direction_stop())
3487 pathcachekey = (phase.definition(), layer_start, layer_stop)
3489 if pathcachekey in self._pathcache:
3490 phase_paths = self._pathcache[pathcachekey]
3491 else:
3492 hwknee = phase.headwave_knee()
3493 if hwknee:
3494 name_or_z = hwknee.depth
3495 interface = self.discontinuity(name_or_z)
3496 mode = hwknee.in_mode
3497 in_direction = hwknee.direction
3499 pabove, pbelow = interface.critical_ps(mode)
3501 p = min_not_none(pabove, pbelow)
3503 # diffracted wave:
3504 if in_direction == DOWN and (
3505 pbelow is None or pbelow >= pabove):
3507 p *= (1.0 - eps)
3509 path = self.path(p, phase, layer_start, layer_stop)
3510 path.set_prange(p, p, 1.)
3512 phase_paths = [path]
3514 else:
3515 try:
3516 pmax_start = max([
3517 radius(z)/layer_start.v(phase.first_leg().mode, z)
3518 for z in (layer_start.ztop, layer_start.zbot)])
3520 pmax_stop = max([
3521 radius(z)/layer_stop.v(phase.last_leg().mode, z)
3522 for z in (layer_stop.ztop, layer_stop.zbot)])
3524 pmax = min(pmax_start, pmax_stop)
3526 pedges = [0.]
3527 for layer in self.layers():
3528 for z in (layer.ztop, layer.zbot):
3529 for mode in (P, S):
3530 for eps2 in [eps]:
3531 v = layer.v(mode, z)
3532 if v != 0.0:
3533 p = radius(z)/v
3534 if p <= pmax:
3535 pedges.append(p*(1.0-eps2))
3536 pedges.append(p)
3537 pedges.append(p*(1.0+eps2))
3539 pedges = num.unique(sorted(pedges))
3541 phase_paths = {}
3542 cached = {}
3543 counter = [0]
3545 def p_to_path(p):
3546 if p in cached:
3547 return cached[p]
3549 try:
3550 counter[0] += 1
3551 path = self.path(
3552 p, phase, layer_start, layer_stop)
3554 if path not in phase_paths:
3555 phase_paths[path] = []
3557 phase_paths[path].append(p)
3559 except PathFailed:
3560 path = None
3562 cached[p] = path
3563 return path
3565 def recurse(pmin, pmax, i=0):
3566 if i > self._pdepth:
3567 return
3568 path1 = p_to_path(pmin)
3569 path2 = p_to_path(pmax)
3570 if path1 is None and path2 is None and i > 0:
3571 return
3572 if path1 is None or path2 is None or \
3573 hash(path1) != hash(path2):
3575 recurse(pmin, (pmin+pmax)/2., i+1)
3576 recurse((pmin+pmax)/2., pmax, i+1)
3578 for (pl, ph) in zip(pedges[:-1], pedges[1:]):
3579 recurse(pl, ph)
3581 for path, ps in phase_paths.items():
3582 path.set_prange(
3583 min(ps), max(ps), pmax/(self._np-1))
3585 phase_paths = list(phase_paths.keys())
3587 except ZeroDivisionError:
3588 phase_paths = []
3590 self._pathcache[pathcachekey] = phase_paths
3592 paths.extend(phase_paths)
3594 paths.sort(key=lambda x: x.pmin())
3595 return paths
3597 def arrivals(
3598 self,
3599 distances=[],
3600 phases=PhaseDef('P'),
3601 zstart=0.0,
3602 zstop=0.0,
3603 refine=True):
3605 '''
3606 Compute rays and traveltimes for given distances.
3608 :param distances: list or array of distances [deg]
3609 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3610 Comma-separated strings and lists of such strings are also accepted
3611 and are converted to :py:class:`PhaseDef` objects for convenience.
3612 :param zstart: source depth [m]
3613 :param zstop: receiver depth [m]
3614 :param refine: bool flag, whether to use bisectioning to improve
3615 (p, x, t) estimated from interpolation
3616 :returns: a list of :py:class:`Ray` objects, sorted by
3617 (distance, arrival time)
3618 '''
3620 distances = num.asarray(distances, dtype=float)
3622 arrivals = []
3623 for path in self.gather_paths(phases, zstart=zstart, zstop=zstop):
3625 endgaps = path.endgaps(zstart, zstop)
3626 for x, p, t, draft_pxt in path.interpolate_x2pt_linear(
3627 distances, endgaps):
3629 arrivals.append(Ray(path, p, x, t, endgaps, draft_pxt))
3631 if refine:
3632 refined = []
3633 for ray in arrivals:
3635 if ray.path._is_headwave:
3636 refined.append(ray)
3638 try:
3639 ray.refine()
3640 refined.append(ray)
3642 except RefineFailed:
3643 pass
3645 arrivals = refined
3647 arrivals.sort(key=lambda x: (x.x, x.t))
3648 return arrivals
3650 @classmethod
3651 def from_scanlines(cls, producer):
3652 '''
3653 Create layer cake model from sequence of materials at depths.
3655 :param producer: iterable yielding (depth, material, name) tuples
3657 Creates a new :py:class:`LayeredModel` object and uses its
3658 :py:meth:`append` method to add layers and discontinuities as needed.
3659 '''
3661 self = cls()
3662 for z, material, name in producer:
3664 if not self._elements:
3665 self.append(Surface(z, material))
3666 else:
3667 element = self._elements[-1]
3668 if self.zeq(element.zbot, z):
3669 assert isinstance(element, Layer)
3670 self.append(
3671 Interface(z, element.mbot, material, name=name))
3673 else:
3674 if isinstance(element, Discontinuity):
3675 ztop = element.z
3676 mtop = element.mbelow
3677 elif isinstance(element, Layer):
3678 ztop = element.zbot
3679 mtop = element.mbot
3681 if mtop == material:
3682 layer = HomogeneousLayer(
3683 ztop, z, material, name=name)
3684 else:
3685 layer = GradientLayer(
3686 ztop, z, mtop, material, name=name)
3688 self.append(layer)
3690 return self
3692 def to_scanlines(self, get_burgers=False):
3693 def fmt(z, m):
3694 if not m._has_default_burgers() or get_burgers:
3695 return (z, m.vp, m.vs, m.rho, m.qp, m.qs,
3696 m.burger_eta1, m.burger_eta2, m.burger_valpha)
3697 return (z, m.vp, m.vs, m.rho, m.qp, m.qs)
3699 last = None
3700 lines = []
3701 for element in self.elements():
3702 if isinstance(element, Layer):
3703 if not isinstance(last, Layer):
3704 lines.append(fmt(element.ztop, element.mtop))
3706 lines.append(fmt(element.zbot, element.mbot))
3708 last = element
3710 if not isinstance(last, Layer):
3711 lines.append(fmt(last.z, last.mbelow))
3713 return lines
3715 def iter_material_parameter(self, get):
3716 assert get in ('vp', 'vs', 'rho', 'qp', 'qs', 'z')
3717 if get == 'z':
3718 for layer in self.layers():
3719 yield layer.ztop
3720 yield layer.zbot
3721 else:
3722 getter = operator.attrgetter(get)
3723 for layer in self.layers():
3724 yield getter(layer.mtop)
3725 yield getter(layer.mbot)
3727 def profile(self, get):
3728 '''
3729 Get parameter profile along depth of the earthmodel.
3731 :param get: property to be queried (
3732 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3733 :type get: str
3734 '''
3736 return num.array(list(self.iter_material_parameter(get)))
3738 def min(self, get='vp'):
3739 '''
3740 Find minimum value of a material property or depth.
3742 :param get: property to be queried (
3743 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3744 '''
3746 return min(self.iter_material_parameter(get))
3748 def max(self, get='vp'):
3749 '''
3750 Find maximum value of a material property or depth.
3752 :param get: property to be queried (
3753 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, ``'qs'``, or ``'z'``)
3754 '''
3756 return max(self.iter_material_parameter(get))
3758 def simplify_layers(self, layers, max_rel_error=0.001):
3759 if len(layers) <= 1:
3760 return layers
3762 ztop = layers[0].ztop
3763 zbot = layers[-1].zbot
3764 zorigs = [layer.ztop for layer in layers]
3765 zorigs.append(zbot)
3766 zs = num.linspace(ztop, zbot, 100)
3767 data = []
3768 for z in zs:
3769 if z == ztop:
3770 direction = UP
3771 else:
3772 direction = DOWN
3774 mat = self.material(z, direction)
3775 data.append(mat.astuple())
3777 data = num.array(data, dtype=float)
3778 data_means = num.mean(data, axis=0)
3779 nmax = len(layers) // 2
3780 accept = False
3782 zcut_best = []
3783 for n in range(1, nmax+1):
3784 ncutintervals = 20
3785 zdelta = (zbot-ztop)/ncutintervals
3786 if n == 2:
3787 zcuts = [
3788 [ztop, ztop + i*zdelta, zbot]
3789 for i in range(1, ncutintervals)]
3790 elif n == 3:
3791 zcuts = []
3792 for j in range(1, ncutintervals):
3793 for i in range(j+1, ncutintervals):
3794 zcuts.append(
3795 [ztop, ztop + j*zdelta, ztop + i*zdelta, zbot])
3796 else:
3797 zcuts = []
3798 zcuts.append(num.linspace(ztop, zbot, n+1))
3799 if zcut_best:
3800 zcuts.append(sorted(num.linspace(
3801 ztop, zbot, n).tolist() + zcut_best[1]))
3802 zcuts.append(sorted(num.linspace(
3803 ztop, zbot, n-1).tolist() + zcut_best[2]))
3805 best = None
3806 for icut, zcut in enumerate(zcuts):
3807 rel_par_errors = num.zeros(5)
3808 mpar_nodes = num.zeros((n+1, 5))
3810 for ipar in range(5):
3811 znodes, vnodes, error_rms = util.polylinefit(
3812 zs, data[:, ipar], zcut)
3814 mpar_nodes[:, ipar] = vnodes
3815 if data_means[ipar] == 0.0:
3816 rel_par_errors[ipar] = -1
3817 else:
3818 rel_par_errors[ipar] = error_rms/data_means[ipar]
3820 rel_error = rel_par_errors.max()
3821 if best is None or rel_error < best[0]:
3822 best = (rel_error, zcut, mpar_nodes)
3824 rel_error, zcut, mpar_nodes = best
3826 zcut_best.append(list(zcut))
3827 zcut_best[-1].pop(0)
3828 zcut_best[-1].pop()
3830 if rel_error <= max_rel_error:
3831 accept = True
3832 break
3834 if not accept:
3835 return layers
3837 rel_error, zcut, mpar_nodes = best
3839 material_nodes = []
3840 for i in range(n+1):
3841 material_nodes.append(Material(*mpar_nodes[i, :]))
3843 out_layers = []
3844 for i in range(n):
3845 mtop = material_nodes[i]
3846 mbot = material_nodes[i+1]
3847 ztop = zcut[i]
3848 zbot = zcut[i+1]
3849 if mtop == mbot:
3850 lyr = HomogeneousLayer(ztop, zbot, mtop)
3851 else:
3852 lyr = GradientLayer(ztop, zbot, mtop, mbot)
3854 out_layers.append(lyr)
3855 return out_layers
3857 def simplify(self, max_rel_error=0.001):
3858 '''
3859 Get representation of model with lower resolution.
3861 Returns an approximation of the model. All discontinuities are kept,
3862 but layer stacks with continuous model parameters are represented, if
3863 possible, by a lower number of layers. Piecewise linear functions are
3864 fitted against the original model parameter's piecewise linear
3865 functions. Successively larger numbers of layers are tried, until the
3866 difference to the original model is below ``max_rel_error``. The
3867 difference is measured as the RMS error of the fit normalized by the
3868 mean of the input (i.e. the fitted curves should deviate, on average,
3869 less than 0.1% from the input curves if ``max_rel_error`` = 0.001).
3870 '''
3872 mod_simple = LayeredModel()
3874 glayers = []
3875 for element in self.elements():
3877 if isinstance(element, Discontinuity):
3878 for layer in self.simplify_layers(
3879 glayers, max_rel_error=max_rel_error):
3881 mod_simple.append(layer)
3883 glayers = []
3884 mod_simple.append(element)
3885 else:
3886 glayers.append(element)
3888 for layer in self.simplify_layers(
3889 glayers, max_rel_error=max_rel_error):
3890 mod_simple.append(layer)
3892 return mod_simple
3894 def extract(self, depth_min=None, depth_max=None):
3895 '''
3896 Extract :py:class:`LayeredModel` from :py:class:`LayeredModel`.
3898 :param depth_min: depth of upper cut or name of :py:class:`Interface`
3899 :param depth_max: depth of lower cut or name of :py:class:`Interface`
3901 Interpolates a :py:class:`GradientLayer` at ``depth_min`` and/or
3902 ``depth_max``.
3903 '''
3905 if isinstance(depth_min, str):
3906 depth_min = self.discontinuity(depth_min).z
3908 if isinstance(depth_max, str):
3909 depth_max = self.discontinuity(depth_max).z
3911 mod_extracted = LayeredModel()
3913 for element in self.elements():
3914 element = element.copy()
3915 do_append = False
3916 if (depth_min is None or depth_min <= element.ztop) \
3917 and (depth_max is None or depth_max >= element.zbot):
3918 mod_extracted.append(element)
3919 continue
3921 if depth_min is not None:
3922 if element.ztop < depth_min and depth_min < element.zbot:
3923 _, element = element.split(depth_min)
3924 do_append = True
3926 if depth_max is not None:
3927 if element.zbot > depth_max and depth_max > element.ztop:
3928 element, _ = element.split(depth_max)
3929 do_append = True
3931 if do_append:
3932 mod_extracted.append(element)
3934 return mod_extracted
3936 def replaced_crust(self, crust2_profile=None, crustmod=None):
3937 if crust2_profile is not None:
3938 profile = anything_to_crust2_profile(crust2_profile)
3939 crustmod = LayeredModel.from_scanlines(
3940 from_crust2x2_profile(profile))
3942 newmod = LayeredModel()
3943 for element in crustmod.extract(depth_max='moho').elements():
3944 if element.name != 'moho':
3945 newmod.append(element)
3946 else:
3947 moho1 = element
3949 mod = self.extract(depth_min='moho')
3950 first = True
3951 for element in mod.elements():
3952 if element.name == 'moho':
3953 if element.z <= moho1.z:
3954 mbelow = mod.material(moho1.z, direction=UP)
3955 else:
3956 mbelow = element.mbelow
3958 moho = Interface(moho1.z, moho1.mabove, mbelow, name='moho')
3959 newmod.append(moho)
3960 else:
3961 if first:
3962 if isinstance(element, Layer) and element.zbot > moho.z:
3963 newmod.append(GradientLayer(
3964 moho.z,
3965 element.zbot,
3966 moho.mbelow,
3967 element.mbot,
3968 name=element.name))
3970 first = False
3971 else:
3972 newmod.append(element)
3973 return newmod
3975 def perturb(self, rstate=None, keep_vp_vs=False, zmax=None, **kwargs):
3976 '''
3977 Create a perturbed variant of the earth model.
3979 Randomly change the thickness and material parameters of the earth
3980 model from a uniform distribution.
3982 :param kwargs: Maximum change fraction (e.g. 0.1) of the parameters.
3983 Name the parameter, prefixed by ``p``. Supported parameters are
3984 ``ph, pvp, pvs, prho, pqs, pqp``.
3985 :type kwargs: dict
3986 :param rstate: Random state to draw from, defaults to ``None``
3987 :type rstate: :class:`numpy.random.RandomState`
3988 :param keep_vp_vs: Keep the Vp/Vs ratio, defaults to False
3989 :type keep_vp_vs: bool
3991 :returns: A new, perturbed earth model
3992 :rtype: :class:`~pyrocko.cake.LayeredModel`
3994 .. code-block :: python
3996 perturbed_model = model.perturb(ph=.1, pvp=.05, prho=.1)
3997 '''
3999 if rstate is None:
4000 rstate = num.random.RandomState()
4002 def factor(k):
4003 try:
4004 pval = kwargs[k]
4005 return rstate.uniform(1.0-pval, 1.0+pval)
4006 except KeyError:
4007 return 1.0
4009 def perturb_material(mat):
4010 mat_new = copy.deepcopy(mat)
4011 for param in kwargs:
4012 if param != 'ph':
4013 value = getattr(mat, param[1:])
4014 setattr(mat_new, param[1:], value*factor(param))
4016 if keep_vp_vs:
4017 mat_new.vs = mat_new.vp * (mat.vs / mat.vp)
4019 return mat_new
4021 model_new = LayeredModel()
4022 elements = list(self.elements())
4023 assert isinstance(elements[0], Surface)
4024 z = elements[0].z
4025 mat = None
4026 for element in elements:
4027 element_new = copy.deepcopy(element)
4028 if isinstance(element_new, Discontinuity):
4029 element_new.z = z
4031 elif isinstance(element_new, Layer):
4032 element_new.ztop = z
4033 if zmax is None or element.zbot < zmax:
4034 z += (element_new.zbot-element_new.ztop) * factor('ph')
4035 else:
4036 z += (element_new.zbot-element_new.ztop)
4037 element_new.zbot = z
4039 if mat is None or element.mtop != mat:
4040 mat = element.mtop
4041 if zmax is None or element.ztop < zmax:
4042 mat_new = perturb_material(mat)
4043 else:
4044 mat_new = copy.deepcopy(mat)
4046 element_new.mtop = mat_new
4048 if element.mbot != mat:
4049 mat = element.mbot
4050 if zmax is None or element.zbot < zmax:
4051 mat_new = perturb_material(mat)
4052 else:
4053 mat_new = copy.deepcopy(mat)
4055 element_new.mbot = mat_new
4057 model_new.append(element_new)
4059 return model_new
4061 def require_homogeneous(self):
4062 elements = list(self.elements())
4064 if len(elements) != 2:
4065 raise LayeredModelError(
4066 'Homogeneous model required but found more than one layer in '
4067 'earthmodel.')
4069 if not isinstance(elements[1], HomogeneousLayer):
4070 raise LayeredModelError(
4071 'Homogeneous model required but element #1 is not of type '
4072 'HomogeneousLayer.')
4074 return elements[1].m
4076 def __str__(self):
4077 return '\n'.join(str(element) for element in self._elements)
4080def read_hyposat_model(fn):
4081 '''
4082 Reader for HYPOSAT earth model files.
4084 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
4086 Interface names are translated as follows: ``'MOHO'`` -> ``'moho'``,
4087 ``'CONR'`` -> ``'conrad'``
4088 '''
4090 with open(fn, 'r') as f:
4091 translate = {'MOHO': 'moho', 'CONR': 'conrad'}
4092 lname = None
4093 for iline, line in enumerate(f):
4094 if iline == 0:
4095 continue
4097 z, vp, vs, name = util.unpack_fixed('f10, f10, f10, a4', line)
4098 if not name:
4099 name = None
4100 material = Material(vp*1000., vs*1000.)
4102 tname = translate.get(lname, lname)
4103 yield z*1000., material, tname
4105 lname = name
4108def read_nd_model(fn):
4109 '''
4110 Reader for TauP style '.nd' (named discontinuity) files.
4112 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
4114 Interface names are translated as follows: ``'mantle'`` -> ``'moho'``,
4115 ``'outer-core'`` -> ``'cmb'``, ``'inner-core'`` -> ``'icb'``.
4117 The format has been modified to include Burgers materials parameters in
4118 columns 7 (burger_eta1), 8 (burger_eta2) and 9. eta(3).
4119 '''
4120 with open(fn, 'r') as f:
4121 for x in read_nd_model_fh(f):
4122 yield x
4125def read_nd_model_str(s):
4126 f = StringIO(s)
4127 for x in read_nd_model_fh(f):
4128 yield x
4129 f.close()
4132def read_nd_model_fh(f):
4133 translate = {'mantle': 'moho', 'outer-core': 'cmb', 'inner-core': 'icb'}
4134 name = None
4135 for line in f:
4136 toks = line.split()
4137 if len(toks) == 9 or len(toks) == 6 or len(toks) == 4:
4138 z, vp, vs, rho = [float(x) for x in toks[:4]]
4139 qp, qs = None, None
4140 burgers = None
4141 if len(toks) == 6 or len(toks) == 9:
4142 qp, qs = [float(x) for x in toks[4:6]]
4143 if len(toks) == 9:
4144 burgers = \
4145 [float(x) for x in toks[6:]]
4147 material = Material(
4148 vp*1000., vs*1000., rho*1000., qp, qs,
4149 burgers=burgers)
4151 yield z*1000., material, name
4152 name = None
4153 elif len(toks) == 1:
4154 name = translate.get(toks[0], toks[0])
4156 f.close()
4159def from_crust2x2_profile(profile, depthmantle=50000):
4160 from pyrocko.dataset import crust2x2
4162 default_qp_qs = {
4163 'soft sed.': (50., 50.),
4164 'hard sed.': (200., 200.),
4165 'upper crust': (600., 400.),
4166 }
4168 z = 0.
4169 for i in range(8):
4170 dz, vp, vs, rho = profile.get_layer(i)
4171 name = crust2x2.Crust2Profile.layer_names[i]
4172 if name in default_qp_qs:
4173 qp, qs = default_qp_qs[name]
4174 else:
4175 qp, qs = None, None
4177 material = Material(vp, vs, rho, qp, qs)
4178 iname = None
4179 if i == 7:
4180 iname = 'moho'
4181 if dz != 0.0:
4182 yield z, material, iname
4183 if i != 7:
4184 yield z+dz, material, name
4185 else:
4186 yield z+depthmantle, material, name
4188 z += dz
4191def write_nd_model_fh(mod, fh):
4192 def fmt(z, mat):
4193 rstr = ' '.join(
4194 util.gform(x, 4)
4195 for x in (
4196 z/1000.,
4197 mat.vp/1000.,
4198 mat.vs/1000.,
4199 mat.rho/1000.,
4200 mat.qp, mat.qs))
4201 if not mat._has_default_burgers():
4202 rstr += ' '.join(
4203 util.gform(x, 4)
4204 for x in (
4205 mat.burger_eta1,
4206 mat.burger_eta2,
4207 mat.burger_valpha))
4208 return rstr.rstrip() + '\n'
4210 translate = {
4211 'moho': 'mantle',
4212 'cmb': 'outer-core',
4213 'icb': 'inner-core'}
4215 last = None
4216 for element in mod.elements():
4217 if isinstance(element, Interface):
4218 if element.name is not None:
4219 n = translate.get(element.name, element.name)
4220 fh.write('%s\n' % n)
4222 elif isinstance(element, Layer):
4223 if not isinstance(last, Layer):
4224 fh.write(fmt(element.ztop, element.mtop))
4226 fh.write(fmt(element.zbot, element.mbot))
4228 last = element
4230 if not isinstance(last, Layer):
4231 fh.write(fmt(last.z, last.mbelow))
4234def write_nd_model_str(mod):
4235 f = StringIO()
4236 write_nd_model_fh(mod, f)
4237 return f.getvalue()
4240def write_nd_model(mod, fn):
4241 with open(fn, 'w') as f:
4242 write_nd_model_fh(mod, f)
4245def builtin_models():
4246 return sorted([
4247 os.path.splitext(os.path.basename(x))[0]
4248 for x in glob.glob(builtin_model_filename('*'))])
4251def builtin_model_filename(modelname):
4252 return util.data_file(os.path.join('earthmodels', modelname+'.nd'))
4255def load_model(fn='ak135-f-continental.m', format='nd', crust2_profile=None):
4256 '''
4257 Load layered earth model from file.
4259 :param fn: filename
4260 :param format: format
4261 :param crust2_profile: ``(lat, lon)`` or
4262 :py:class:`pyrocko.dataset.crust2x2.Crust2Profile` object, merge model
4263 with crustal profile. If ``fn`` is forced to be ``None`` only the
4264 converted CRUST2.0 profile is returned.
4265 :returns: object of type :py:class:`LayeredModel`
4267 The following formats are currently supported:
4269 ============== ===========================================================
4270 format description
4271 ============== ===========================================================
4272 ``'nd'`` 'named discontinuity' format used by the TauP programs
4273 ``'hyposat'`` format used by the HYPOSAT location program
4274 ============== ===========================================================
4276 The naming of interfaces is translated from the file format's native naming
4277 to Cake's own convention (See :py:func:`read_nd_model` and
4278 :py:func:`read_hyposat_model` for details). Cake likes the following
4279 internal names: ``'conrad'``, ``'moho'``, ``'cmb'`` (core-mantle boundary),
4280 ``'icb'`` (inner core boundary).
4281 '''
4283 if fn is not None:
4284 if format == 'nd':
4285 if not os.path.exists(fn) and fn in builtin_models():
4286 fn = builtin_model_filename(fn)
4287 reader = read_nd_model(fn)
4288 elif format == 'hyposat':
4289 reader = read_hyposat_model(fn)
4290 else:
4291 assert False, 'unsupported model format'
4293 mod = LayeredModel.from_scanlines(reader)
4294 if crust2_profile is not None:
4295 return mod.replaced_crust(crust2_profile)
4297 return mod
4299 else:
4300 assert crust2_profile is not None
4301 profile = anything_to_crust2_profile(crust2_profile)
4302 return LayeredModel.from_scanlines(
4303 from_crust2x2_profile(profile))
4306def castagna_vs_to_vp(vs):
4307 '''
4308 Calculate vp from vs using Castagna's relation.
4310 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4311 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4312 1985]
4314 vp = 1.16 * vs + 1360 [m/s]
4316 :param vs: S-wave velocity [m/s]
4317 :returns: P-wave velocity [m/s]
4318 '''
4320 return vs*1.16 + 1360.0
4323def castagna_vp_to_vs(vp):
4324 '''
4325 Calculate vp from vs using Castagna's relation.
4327 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4328 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4329 1985]
4331 vp = 1.16 * vs + 1360 [m/s]
4333 :param vp: P-wave velocity [m/s]
4334 :returns: S-wave velocity [m/s]
4335 '''
4337 return (vp - 1360.0) / 1.16
4340def evenize(x, y, minsize=10):
4341 if x.size < minsize:
4342 return x
4343 ry = (y.max()-y.min())
4344 if ry == 0:
4345 return x
4346 dx = (x[1:] - x[:-1])/(x.max()-x.min())
4347 dy = (y[1:] + y[:-1])/ry
4349 s = num.zeros(x.size)
4350 s[1:] = num.cumsum(num.sqrt(dy**2 + dx**2))
4351 s2 = num.linspace(0, s[-1], x.size)
4352 x2 = num.interp(s2, s, x)
4353 x2[0] = x[0]
4354 x2[-1] = x[-1]
4355 return x2
4358def next_or_none(i):
4359 try:
4360 return next(i)
4361 except StopIteration:
4362 return None
4365def reci_or_none(x):
4366 try:
4367 return 1./x
4368 except ZeroDivisionError:
4369 return None
4372def mult_or_none(a, b):
4373 if a is None or b is None:
4374 return None
4375 return a*b
4378def min_not_none(a, b):
4379 if a is None:
4380 return b
4381 if b is None:
4382 return a
4383 return min(a, b)
4386def xytups(xx, yy):
4387 d = []
4388 for x, y in zip(xx, yy):
4389 if num.isfinite(y):
4390 d.append((x, y))
4391 return d
4394def interp(x, xp, fp, monoton):
4395 if monoton == 1:
4396 return xytups(
4397 x, num.interp(x, xp, fp, left=num.nan, right=num.nan))
4398 elif monoton == -1:
4399 return xytups(
4400 x, num.interp(x, xp[::-1], fp[::-1], left=num.nan, right=num.nan))
4401 else:
4402 fs = []
4403 for xv in x:
4404 indices = num.where(num.logical_or(
4405 num.logical_and(xp[:-1] >= xv, xv > xp[1:]),
4406 num.logical_and(xp[:-1] <= xv, xv < xp[1:])))[0]
4408 for i in indices:
4409 xr = (xv - xp[i])/(xp[i+1]-xp[i])
4410 fv = xr*fp[i] + (1.-xr)*fp[i+1]
4411 fs.append((xv, fv))
4413 return fs
4416def float_or_none(x):
4417 if x is not None:
4418 return float(x)
4421def parstore_float(thelocals, obj, *args):
4422 for k, v in thelocals.items():
4423 if k != 'self' and (not args or k in args):
4424 setattr(obj, k, float_or_none(v))