Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/cake.py: 85%
2023 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-23 12:04 +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) == 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) == 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) == 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) == 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 filled(v, len(z))
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) == 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 filled(self.interface.z, len(p))
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):
2609 '''
2610 Get geometrical representation of ray path.
2611 '''
2613 if self._is_headwave:
2614 assert p.size == 1
2615 x, t = self.xt(p, endgaps)
2616 xstretch = x_for_headwave-x
2617 nout = xstretch.size
2618 else:
2619 nout = p.size
2621 dxl, dtl = self.xt_endgaps(p, endgaps, which='left')
2622 dxr, dtr = self.xt_endgaps(p, endgaps, which='right')
2624 # first create full path including the endgaps
2625 sx = num.zeros(nout) - dxl
2626 st = num.zeros(nout) - dtl
2627 zxt = []
2628 for s in self.straights():
2629 n = points_per_straight
2631 back = None
2632 zin, zout = s.z_in(), s.z_out()
2633 if type(s) is HeadwaveStraight:
2634 z = zin
2635 for i in range(n):
2636 xs = float(i)/(n-1) * xstretch
2637 ts = s.x2t_headwave(xs)
2638 zxt.append((filled(z, xstretch.size), sx+xs, st+ts))
2639 else:
2640 if zin != zout: # normal traversal
2641 zs = num.linspace(zin, zout, n).tolist()
2642 for z in zs:
2643 x, t = s.xt(p, zpart=sorted([zin, z]))
2644 zxt.append((filled(z, nout), sx + x, st + t))
2646 else: # ray turns in layer
2647 zturn = s.zturn(p)
2648 back = []
2649 for i in range(n):
2650 z = zin + (zturn - zin) * num.sin(
2651 float(i)/(n-1)*math.pi/2.0) * 0.999
2653 if zturn[0] >= zin:
2654 x, t = s.xt(p, zpart=[zin, z])
2655 else:
2656 x, t = s.xt(p, zpart=[z, zin])
2657 zxt.append((z, sx + x, st + t))
2658 back.append((z, x, t))
2660 if type(s) is HeadwaveStraight:
2661 x = xstretch
2662 t = s.x2t_headwave(xstretch)
2663 else:
2664 x, t = s.xt(p)
2666 sx += x
2667 st += t
2668 if back:
2669 for z, x, t in reversed(back):
2670 zxt.append((z, sx - x, st - t))
2672 # gather results as arrays with such that x[ip, ipoint]
2673 fanz, fanx, fant = [], [], []
2674 for z, x, t in zxt:
2675 fanz.append(z)
2676 fanx.append(x)
2677 fant.append(t)
2679 z = num.array(fanz).T
2680 x = num.array(fanx).T
2681 t = num.array(fant).T
2683 # cut off the endgaps, add exact endpoints
2684 xmax = x[:, -1] - dxr
2685 tmax = t[:, -1] - dtr
2686 zstart, zstop = endgaps[:2]
2687 zs, xs, ts = [], [], []
2688 for i in range(nout):
2689 t_ = t[i]
2690 indices = num.where(num.logical_and(0. <= t_, t_ <= tmax[i]))[0]
2691 n = indices.size + 2
2692 zs_, xs_, ts_ = [num.empty(n, dtype=float) for j in range(3)]
2693 zs_[1:-1] = z[i, indices]
2694 xs_[1:-1] = x[i, indices]
2695 ts_[1:-1] = t[i, indices]
2696 zs_[0], zs_[-1] = zstart, zstop
2697 xs_[0], xs_[-1] = 0., xmax[i]
2698 ts_[0], ts_[-1] = 0., tmax[i]
2699 zs.append(zs_)
2700 xs.append(xs_)
2701 ts.append(ts_)
2703 return zs, xs, ts
2705 def _analyse(self):
2706 if self._p is not None:
2707 return
2709 p = self.make_p(nmin=20)
2710 xmin, xmax, tmin, tmax = self.xt_limits(p)
2712 self._x, self._t, self._p = xmax, tmax, p
2713 self._xmin, self._xmax = xmin.min(), xmax.max()
2714 self._tmin, self._tmax = tmin.min(), tmax.max()
2716 def draft_pxt(self, endgaps):
2717 self._analyse()
2719 if not self._is_headwave:
2720 cp, cx, ct = self._p, self._x, self._t
2721 pcrit = min(
2722 self.critical_pstart(endgaps),
2723 self.critical_pstop(endgaps))
2725 if pcrit < self._pmin:
2726 empty = num.array([], dtype=float)
2727 return empty, empty, empty
2729 elif pcrit >= self._pmax:
2730 dx, dt = self.xt_endgaps(cp, endgaps)
2731 return cp, cx-dx, ct-dt
2733 else:
2734 n = num.searchsorted(cp, pcrit) + 1
2735 rp, rx, rt = num.empty((3, n), dtype=float)
2736 rp[:-1] = cp[:n-1]
2737 rx[:-1] = cx[:n-1]
2738 rt[:-1] = ct[:n-1]
2739 rp[-1] = pcrit
2740 rx[-1], rt[-1] = self.xt(pcrit, endgaps)
2741 dx, dt = self.xt_endgaps(rp, endgaps)
2742 rx[:-1] -= dx[:-1]
2743 rt[:-1] -= dt[:-1]
2744 return rp, rx, rt
2746 else:
2747 dx, dt = self.xt_endgaps(self._p, endgaps)
2748 p, x, t = self._p, self._x - dx, self._t - dt
2749 p, x, t = p[0], x[0], t[0]
2750 xh = num.linspace(0., x*10-x, 10)
2751 th = self.headwave_straight().x2t_headwave(xh)
2752 return filled(p, xh.size), x+xh, t+th
2754 def interpolate_x2pt_linear(self, x, endgaps):
2755 '''
2756 Get approximate ray parameter and traveltime for distance.
2757 '''
2759 self._analyse()
2761 if self._is_headwave:
2762 dx, dt = self.xt_endgaps(self._p, endgaps)
2763 xmin = self._x[0] - dx[0]
2764 tmin = self._t[0] - dt[0]
2765 el = self.headwave_straight()
2766 xok = x[x >= xmin]
2767 th = el.x2t_headwave(xstretch=(xok-xmin)) + tmin
2768 return [
2769 (x_, self._p[0], t, None) for (x_, t) in zip(xok, th)]
2771 else:
2772 if num.all(x < self._xmin) or num.all(self._xmax < x):
2773 return []
2775 rp, rx, rt = self.draft_pxt(endgaps)
2777 xp = interp(x, rx, rp, 0)
2778 xt = interp(x, rx, rt, 0)
2780 if (rp.size and
2781 len(xp) == 0 and
2782 rx[0] == 0.0 and
2783 any(x == 0.0) and
2784 rp[0] == 0.0):
2786 xp = [(0.0, rp[0])]
2787 xt = [(0.0, rt[0])]
2789 return [
2790 (x_, p, t, (rp, rx, rt)) for ((x_, p), (_, t)) in zip(xp, xt)]
2792 def __eq__(self, other):
2793 if len(self.elements) != len(other.elements):
2794 return False
2796 return all(a == b for a, b in zip(self.elements, other.elements))
2798 def __hash__(self):
2799 return hash(
2800 tuple(hash(x) for x in self.elements) +
2801 (self.phase.definition(), ))
2803 def __str__(self, p=None, eps=1.):
2804 x = []
2805 start_i = None
2806 end_i = None
2807 turn_i = None
2809 def append_layers(si, ei, ti):
2810 if si == ei and (ti is None or ti == si):
2811 x.append('%i' % si)
2812 else:
2813 if ti is not None:
2814 x.append('(%i-%i-%i)' % (si, ti, ei))
2815 else:
2816 x.append('(%i-%i)' % (si, ei))
2818 for el in self.elements:
2819 if type(el) is Straight:
2820 if start_i is None:
2821 start_i = el.layer.ilayer
2822 if el._direction_in != el._direction_out:
2823 turn_i = el.layer.ilayer
2824 end_i = el.layer.ilayer
2826 elif isinstance(el, Kink):
2827 if start_i is not None:
2828 append_layers(start_i, end_i, turn_i)
2829 start_i = None
2830 turn_i = None
2832 x.append(str(el))
2834 if start_i is not None:
2835 append_layers(start_i, end_i, turn_i)
2837 su = '(%s)' % self.used_phase(p=p, eps=eps).used_repr()
2839 return '%-15s %-17s %s' % (self.phase.definition(), su, ''.join(x))
2841 def critical_pstart(self, endgaps):
2842 '''
2843 Get critical ray parameter for source depth choice.
2844 '''
2846 return self.first_straight().critical_p_in(endgaps)
2848 def critical_pstop(self, endgaps):
2849 '''
2850 Get critical ray parameter for receiver depth choice.
2851 '''
2853 return self.last_straight().critical_p_out(endgaps)
2855 def ranges(self, endgaps):
2856 '''
2857 Get valid ranges of ray parameter, distance, and traveltime.
2858 '''
2859 p, x, t = self.draft_pxt(endgaps)
2860 return p.min(), p.max(), x.min(), x.max(), t.min(), t.max()
2862 def describe(self, endgaps=None, as_degrees=False):
2863 '''
2864 Get textual representation.
2865 '''
2867 self._analyse()
2869 if as_degrees:
2870 xunit = 'deg'
2871 xfact = 1.
2872 else:
2873 xunit = 'km'
2874 xfact = d2m/km
2876 sg = ''' Ranges for all depths in source and receiver layers:
2877 - x [%g, %g] %s
2878 - t [%g, %g] s
2879 - p [%g, %g] s/deg
2880''' % (
2881 self._xmin*xfact,
2882 self._xmax*xfact,
2883 xunit,
2884 self._tmin,
2885 self._tmax,
2886 self._pmin/r2d,
2887 self._pmax/r2d)
2889 if endgaps is not None:
2890 pmin, pmax, xmin, xmax, tmin, tmax = self.ranges(endgaps)
2891 ss = ''' Ranges for given source and receiver depths:
2892\n - x [%g, %g] %s
2893\n - t [%g, %g] s
2894\n - p [%g, %g] s/deg
2895\n''' % (xmin*xfact, xmax*xfact, xunit, tmin, tmax, pmin/r2d, pmax/r2d)
2897 else:
2898 ss = ''
2900 return '%s\n' % self + ss + sg
2903class RefineFailed(CakeError):
2904 pass
2907class Ray(object):
2908 '''
2909 Representation of a ray with a specific (path, ray parameter, distance,
2910 arrival time) choice.
2912 **Attributes:**
2914 .. py:attribute:: path
2916 :py:class:`RayPath` object containing complete propagation history.
2918 .. py:attribute:: p
2920 Ray parameter (spherical) [s/rad]
2922 .. py:attribute:: x
2924 Radial distance [deg]
2926 .. py:attribute:: t
2928 Traveltime [s]
2930 .. py:attribute:: endgaps
2932 Needed for source/receiver depth adjustments in many
2933 :py:class:`RayPath` methods.
2934 '''
2936 def __init__(self, path, p, x, t, endgaps, draft_pxt):
2937 self.path = path
2938 self.p = p
2939 self.x = x
2940 self.t = t
2941 self.endgaps = endgaps
2942 self.draft_pxt = draft_pxt
2944 def given_phase(self):
2945 '''
2946 Get phase definition which was used to create the ray.
2948 :returns: :py:class:`PhaseDef` object
2949 '''
2951 return self.path.phase
2953 def used_phase(self):
2954 '''
2955 Compute phase definition from propagation path.
2957 :returns: :py:class:`PhaseDef` object
2958 '''
2960 return self.path.used_phase(self.p)
2962 def refine(self):
2963 if self.path._is_headwave:
2964 return
2966 if self.t == 0.0 and self.p == 0.0 and self.x == 0.0:
2967 return
2969 cp, cx, ct = self.draft_pxt
2970 ip = num.searchsorted(cp, self.p)
2971 if not (0 < ip < cp.size):
2972 raise RefineFailed()
2974 pl, ph = cp[ip-1], cp[ip]
2975 p_to_t = {}
2976 i = [0]
2978 def f(p):
2979 i[0] += 1
2980 x, t = self.path.xt(p, self.endgaps)
2981 p_to_t[p] = t
2982 return self.x - x
2984 try:
2985 self.p = brentq(f, pl, ph)
2986 self.t = p_to_t[self.p]
2988 except ValueError:
2989 raise RefineFailed()
2991 def t_and_attributes(self, attributes):
2992 d = {
2993 'takeoff_angle': lambda: self.takeoff_angle(),
2994 'incidence_angle': lambda: self.incidence_angle(),
2995 't': lambda: self.t,
2996 'p': lambda: self.p}
2998 return tuple(d[k]() for k in ['t'] + attributes)
3000 def takeoff_angle(self):
3001 '''
3002 Get takeoff angle of ray.
3004 The angle is returned in [degrees].
3005 '''
3007 return self.path.first_straight().angle_in(self.p, self.endgaps)
3009 def incidence_angle(self):
3010 '''
3011 Get incidence angle of ray.
3013 The angle is returned in [degrees].
3014 '''
3016 return self.path.last_straight().angle_out(self.p, self.endgaps)
3018 def efficiency(self):
3019 '''
3020 Get conversion/reflection efficiency of the ray.
3022 A value between 0 and 1 is returned, reflecting the relative amount of
3023 energy which is transmitted along the ray and not lost by reflections
3024 or conversions.
3025 '''
3027 return self.path.efficiency(self.p)
3029 def spreading(self):
3030 '''
3031 Get geometrical spreading factor.
3032 '''
3034 return self.path.spreading(self.p, self.endgaps)
3036 def surface_sphere(self):
3037 x1, y1 = 0., earthradius - self.endgaps[0]
3038 r2 = earthradius - self.endgaps[1]
3039 x2, y2 = r2*math.sin(self.x*d2r), r2*math.cos(self.x*d2r)
3040 return ((x2-x1)**2 + (y2-y1)**2)*4.0*math.pi
3042 def zxt_path_subdivided(self, points_per_straight=20):
3043 '''
3044 Get geometrical representation of ray path.
3046 Three arrays (depth, distance, time) with points on the ray's path of
3047 propagation are returned. The number of points which are used in each
3048 ray segment (passage through one layer) may be controlled by the
3049 ``points_per_straight`` parameter.
3050 '''
3051 return self.path.zxt_path_subdivided(
3052 num.atleast_1d(self.p), self.endgaps,
3053 points_per_straight=points_per_straight,
3054 x_for_headwave=num.atleast_1d(self.x))
3056 def __str__(self, as_degrees=False):
3057 if as_degrees:
3058 sd = '%6.3g deg' % self.x
3059 else:
3060 sd = '%7.5g km' % (self.x*(d2r*earthradius/km))
3062 return '%7.5g s/deg %s %6.4g s %5.1f %5.1f %3.0f%% %3.0f%% %s' % (
3063 self.p/r2d,
3064 sd,
3065 self.t,
3066 self.takeoff_angle(),
3067 self.incidence_angle(),
3068 100*self.efficiency(),
3069 100*self.spreading()*self.surface_sphere(),
3070 self.path.__str__(p=self.p))
3073def anything_to_crust2_profile(crust2_profile):
3074 from pyrocko.dataset import crust2x2
3075 if isinstance(crust2_profile, tuple):
3076 lat, lon = [float(x) for x in crust2_profile]
3077 return crust2x2.get_profile(lat, lon)
3078 elif isinstance(crust2_profile, str):
3079 return crust2x2.get_profile(crust2_profile)
3080 elif isinstance(crust2_profile, crust2x2.Crust2Profile):
3081 return crust2_profile
3082 else:
3083 assert False, 'crust2_profile must be (lat, lon) a profile ' \
3084 'key or a crust2x2 Profile object)'
3087class DiscontinuityNotFound(CakeError):
3088 def __init__(self, depth_or_name):
3089 CakeError.__init__(self)
3090 self.depth_or_name = depth_or_name
3092 def __str__(self):
3093 return 'Cannot find discontinuity from given depth or name: %s' % \
3094 self.depth_or_name
3097class LayeredModelError(CakeError):
3098 pass
3101class LayeredModel(object):
3102 '''
3103 Representation of a layer cake model.
3105 There are several ways to initialize an instance of this class.
3107 1. Use the module function :py:func:`load_model` to read a model from a
3108 file.
3109 2. Create an empty model with the default constructor and append layers and
3110 discontinuities with the :py:meth:`append` method (from top to bottom).
3111 3. Use the constructor :py:meth:`LayeredModel.from_scanlines`, to
3112 automatically create the :py:class:`Layer` and :py:class:`Discontinuity`
3113 objects from a given velocity profile.
3115 An earth model is represented by as stack of :py:class:`Layer` and
3116 :py:class:`Discontinuity` objects. The method :py:meth:`arrivals` returns
3117 :py:class:`Ray` objects which may be e.g. queried for arrival times of
3118 specific phases. Each ray is associated with a :py:class:`RayPath` object.
3119 Ray objects share common ray paths if they have the same
3120 conversion/reflection/propagation history. Creating the ray path objects is
3121 relatively expensive (this is done in :py:meth:`gather_paths`), but they
3122 are cached for reuse in successive invocations.
3123 '''
3125 def __init__(self):
3126 self._surface_material = None
3127 self._elements = []
3128 self.nlayers = 0
3129 self._np = 10000
3130 self._pdepth = 5
3131 self._pathcache = {}
3133 def copy_with_elevation(self, elevation):
3134 '''
3135 Get copy of the model with surface layer stretched to given elevation.
3137 :param elevation: new surface elevation in [m]
3139 Elevation is positiv upward, contrary to the layered models downward
3140 `z` axis.
3141 '''
3143 c = copy.deepcopy(self)
3144 c._pathcache = {}
3145 surface = c._elements[0]
3146 toplayer = c._elements[1]
3148 assert toplayer.zbot > -elevation
3150 surface.z = -elevation
3151 c._elements[1] = toplayer.copy(ztop=-elevation)
3152 c._elements[1].ilayer = 0
3153 return c
3155 def zeq(self, z1, z2):
3156 return abs(z1-z2) < ZEPS
3158 def append(self, element):
3159 '''
3160 Add a layer or discontinuity at bottom of model.
3162 :param element: object of subclass of :py:class:`Layer` or
3163 :py:class:`Discontinuity`.
3164 '''
3166 if isinstance(element, Layer):
3167 if element.zbot >= earthradius:
3168 element.zbot = earthradius - 1.
3170 if element.ztop >= earthradius:
3171 raise CakeError('Layer deeper than earthradius')
3173 element.ilayer = self.nlayers
3174 self.nlayers += 1
3176 self._elements.append(element)
3178 def elements(self, direction=DOWN):
3179 '''
3180 Iterate over all elements of the model.
3182 :param direction: direction of traversal :py:const:`DOWN` or
3183 :py:const:`UP`.
3185 Objects derived from the :py:class:`Discontinuity` and
3186 :py:class:`Layer` classes are yielded.
3187 '''
3189 if direction == DOWN:
3190 return iter(self._elements)
3191 else:
3192 return reversed(self._elements)
3194 def layers(self, direction=DOWN):
3195 '''
3196 Iterate over all layers of model.
3198 :param direction: direction of traversal :py:const:`DOWN` or
3199 :py:const:`UP`.
3201 Objects derived from the :py:class:`Layer` class are yielded.
3202 '''
3204 if direction == DOWN:
3205 return (el for el in self._elements if isinstance(el, Layer))
3206 else:
3207 return (
3208 el for el in reversed(self._elements) if isinstance(el, Layer))
3210 def layer(self, z, direction=DOWN):
3211 '''
3212 Get layer for given depth.
3214 :param z: depth [m]
3215 :param direction: direction of traversal :py:const:`DOWN` or
3216 :py:const:`UP`.
3218 Returns first layer which touches depth ``z`` (tolerant at boundaries).
3219 '''
3221 for layer in self.layers(direction):
3222 if layer.contains(z):
3223 return layer
3224 else:
3225 raise CakeError('Failed extracting layer at depth z=%s' % z)
3227 def walker(self):
3228 return Walker(self._elements)
3230 def material(self, z, direction=DOWN):
3231 '''
3232 Get material at given depth.
3234 :param z: depth [m]
3235 :param direction: direction of traversal :py:const:`DOWN` or
3236 :py:const:`UP`
3237 :returns: object of type :py:class:`Material`
3239 If given depth ``z`` happens to be at an interface, the material of the
3240 first layer with respect to the the traversal ordering is returned.
3241 '''
3243 lyr = self.layer(z, direction)
3244 return lyr.material(z)
3246 def discontinuities(self):
3247 '''
3248 Iterate over all discontinuities of the model.'''
3250 return (el for el in self._elements if isinstance(el, Discontinuity))
3252 def discontinuity(self, name_or_z):
3253 '''
3254 Get discontinuity by name or depth.
3256 :param name_or_z: name of discontinuity or depth [m] as float value
3257 '''
3259 if isinstance(name_or_z, float):
3260 candi = sorted(
3261 self.discontinuities(), key=lambda i: abs(i.z-name_or_z))
3262 else:
3263 candi = [i for i in self.discontinuities() if i.name == name_or_z]
3265 if not candi:
3266 raise DiscontinuityNotFound(name_or_z)
3268 return candi[0]
3270 def adapt_phase(self, phase):
3271 '''
3272 Adapt a phase definition for use with this model.
3274 This returns a copy of the phase definition, where named
3275 discontinuities are replaced with the actual depth of these, as defined
3276 in the model.
3277 '''
3279 phase = phase.copy()
3280 for knee in phase.knees():
3281 if knee.depth != 'surface':
3282 knee.depth = self.discontinuity(knee.depth).z
3283 for leg in phase.legs():
3284 if leg.depthmax is not None and isinstance(leg.depthmax, str):
3285 leg.depthmax = self.discontinuity(leg.depthmax).z
3287 return phase
3289 def path(self, p, phase, layer_start, layer_stop):
3290 '''
3291 Get ray path for given combination of ray parameter, phase definition,
3292 source and receiver layers.
3294 :param p: ray parameter (spherical) [s/rad]
3295 :param phase: phase definition (:py:class:`PhaseDef` object)
3296 :param layer_start: layer with source
3297 :param layer_stop: layer with receiver
3298 :returns: :py:class:`RayPath` object
3300 If it is not possible to find a solution, an exception of type
3301 :py:exc:`NotPhaseConform`, :py:exc:`MinDepthReached`,
3302 :py:exc:`MaxDepthReached`, :py:exc:`CannotPropagate`,
3303 :py:exc:`BottomReached` or :py:exc:`SurfaceReached` is raised.
3304 '''
3306 phase = self.adapt_phase(phase)
3307 knees = phase.knees()
3308 legs = phase.legs()
3309 next_knee = next_or_none(knees)
3310 leg = next_or_none(legs)
3311 assert leg is not None
3313 direction = leg.departure
3314 direction_stop = phase.direction_stop()
3315 mode = leg.mode
3316 mode_stop = phase.last_leg().mode
3318 walker = self.walker()
3319 walker.goto_layer(layer_start)
3320 current = walker.current()
3322 ttop, tbot = current.tests(p, mode)
3323 if not ttop and not tbot:
3324 raise CannotPropagate(direction, current.ilayer)
3326 if (direction == DOWN and not ttop) or (direction == UP and not tbot):
3327 direction = -direction
3329 path = RayPath(phase)
3330 trapdetect = set()
3331 while True:
3332 at_layer = isinstance(current, Layer)
3333 at_discontinuity = isinstance(current, Discontinuity)
3335 # detect trapped wave
3336 k = (id(next_knee), id(current), direction, mode)
3337 if k in trapdetect:
3338 raise Trapped()
3340 trapdetect.add(k)
3342 if at_discontinuity:
3343 oldmode, olddirection = mode, direction
3344 headwave = False
3345 if next_knee is not None and next_knee.matches(
3346 current, mode, direction):
3348 headwave = next_knee.headwave
3349 direction = next_knee.out_direction()
3350 mode = next_knee.out_mode
3351 next_knee = next_or_none(knees)
3352 leg = next(legs)
3354 else: # implicit reflection/transmission
3355 direction = current.propagate(p, mode, direction)
3357 if headwave:
3358 path.set_is_headwave(True)
3360 path.append(Kink(
3361 olddirection, olddirection, oldmode, oldmode, current))
3363 path.append(HeadwaveStraight(
3364 olddirection, direction, oldmode, current))
3366 path.append(Kink(
3367 olddirection, direction, oldmode, mode, current))
3369 else:
3370 path.append(Kink(
3371 olddirection, direction, oldmode, mode, current))
3373 if at_layer:
3374 direction_in = direction
3375 direction = current.propagate(p, mode, direction_in)
3377 zturn = None
3378 if direction_in != direction:
3379 zturn = current.zturn(p, mode)
3381 zmin, zmax = leg.depthmin, leg.depthmax
3382 if zmin is not None or zmax is not None:
3383 if direction_in != direction:
3384 if zmin is not None and zturn <= zmin:
3385 raise MinDepthReached()
3386 if zmax is not None and zturn >= zmax:
3387 raise MaxDepthReached()
3388 else:
3389 if zmin is not None and current.ztop <= zmin:
3390 raise MinDepthReached()
3391 if zmax is not None and current.zbot >= zmax:
3392 raise MaxDepthReached()
3394 path.append(Straight(direction_in, direction, mode, current))
3396 if next_knee is None and mode == mode_stop and \
3397 current is layer_stop:
3399 if zturn is None:
3400 if direction == direction_stop:
3401 break
3402 else:
3403 break
3405 walker.go(direction)
3406 current = walker.current()
3408 return path
3410 def gather_paths(self, phases=PhaseDef('P'), zstart=0.0, zstop=0.0):
3411 '''
3412 Get all possible ray paths for given source and receiver depths for one
3413 or more phase definitions.
3415 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3416 Comma-separated strings and lists of such strings are also accepted
3417 and are converted to :py:class:`PhaseDef` objects for convenience.
3418 :param zstart: source depth [m]
3419 :param zstop: receiver depth [m]
3420 :returns: a list of :py:class:`RayPath` objects
3422 Results of this method are cached internally. Cached results are
3423 returned, when a given combination of source layer, receiver layer and
3424 phase definition has been used before.
3425 '''
3427 eps = 1e-7 # num.finfo(float).eps * 1000.
3429 phases = to_phase_defs(phases)
3431 paths = []
3432 for phase in phases:
3434 layer_start = self.layer(zstart, -phase.direction_start())
3435 layer_stop = self.layer(zstop, phase.direction_stop())
3437 pathcachekey = (phase.definition(), layer_start, layer_stop)
3439 if pathcachekey in self._pathcache:
3440 phase_paths = self._pathcache[pathcachekey]
3441 else:
3442 hwknee = phase.headwave_knee()
3443 if hwknee:
3444 name_or_z = hwknee.depth
3445 interface = self.discontinuity(name_or_z)
3446 mode = hwknee.in_mode
3447 in_direction = hwknee.direction
3449 pabove, pbelow = interface.critical_ps(mode)
3451 p = min_not_none(pabove, pbelow)
3453 # diffracted wave:
3454 if in_direction == DOWN and (
3455 pbelow is None or pbelow >= pabove):
3457 p *= (1.0 - eps)
3459 path = self.path(p, phase, layer_start, layer_stop)
3460 path.set_prange(p, p, 1.)
3462 phase_paths = [path]
3464 else:
3465 try:
3466 pmax_start = max([
3467 radius(z)/layer_start.v(phase.first_leg().mode, z)
3468 for z in (layer_start.ztop, layer_start.zbot)])
3470 pmax_stop = max([
3471 radius(z)/layer_stop.v(phase.last_leg().mode, z)
3472 for z in (layer_stop.ztop, layer_stop.zbot)])
3474 pmax = min(pmax_start, pmax_stop)
3476 pedges = [0.]
3477 for layer in self.layers():
3478 for z in (layer.ztop, layer.zbot):
3479 for mode in (P, S):
3480 for eps2 in [eps]:
3481 v = layer.v(mode, z)
3482 if v != 0.0:
3483 p = radius(z)/v
3484 if p <= pmax:
3485 pedges.append(p*(1.0-eps2))
3486 pedges.append(p)
3487 pedges.append(p*(1.0+eps2))
3489 pedges = num.unique(sorted(pedges))
3491 phase_paths = {}
3492 cached = {}
3493 counter = [0]
3495 def p_to_path(p):
3496 if p in cached:
3497 return cached[p]
3499 try:
3500 counter[0] += 1
3501 path = self.path(
3502 p, phase, layer_start, layer_stop)
3504 if path not in phase_paths:
3505 phase_paths[path] = []
3507 phase_paths[path].append(p)
3509 except PathFailed:
3510 path = None
3512 cached[p] = path
3513 return path
3515 def recurse(pmin, pmax, i=0):
3516 if i > self._pdepth:
3517 return
3518 path1 = p_to_path(pmin)
3519 path2 = p_to_path(pmax)
3520 if path1 is None and path2 is None and i > 0:
3521 return
3522 if path1 is None or path2 is None or \
3523 hash(path1) != hash(path2):
3525 recurse(pmin, (pmin+pmax)/2., i+1)
3526 recurse((pmin+pmax)/2., pmax, i+1)
3528 for (pl, ph) in zip(pedges[:-1], pedges[1:]):
3529 recurse(pl, ph)
3531 for path, ps in phase_paths.items():
3532 path.set_prange(
3533 min(ps), max(ps), pmax/(self._np-1))
3535 phase_paths = list(phase_paths.keys())
3537 except ZeroDivisionError:
3538 phase_paths = []
3540 self._pathcache[pathcachekey] = phase_paths
3542 paths.extend(phase_paths)
3544 paths.sort(key=lambda x: x.pmin())
3545 return paths
3547 def arrivals(
3548 self,
3549 distances=[],
3550 phases=PhaseDef('P'),
3551 zstart=0.0,
3552 zstop=0.0,
3553 refine=True):
3555 '''
3556 Compute rays and traveltimes for given distances.
3558 :param distances: list or array of distances [deg]
3559 :param phases: a :py:class:`PhaseDef` object or a list of such objects.
3560 Comma-separated strings and lists of such strings are also accepted
3561 and are converted to :py:class:`PhaseDef` objects for convenience.
3562 :param zstart: source depth [m]
3563 :param zstop: receiver depth [m]
3564 :param refine: bool flag, whether to use bisectioning to improve
3565 (p, x, t) estimated from interpolation
3566 :returns: a list of :py:class:`Ray` objects, sorted by
3567 (distance, arrival time)
3568 '''
3570 distances = num.asarray(distances, dtype=float)
3572 arrivals = []
3573 for path in self.gather_paths(phases, zstart=zstart, zstop=zstop):
3575 endgaps = path.endgaps(zstart, zstop)
3576 for x, p, t, draft_pxt in path.interpolate_x2pt_linear(
3577 distances, endgaps):
3579 arrivals.append(Ray(path, p, x, t, endgaps, draft_pxt))
3581 if refine:
3582 refined = []
3583 for ray in arrivals:
3585 if ray.path._is_headwave:
3586 refined.append(ray)
3588 try:
3589 ray.refine()
3590 refined.append(ray)
3592 except RefineFailed:
3593 pass
3595 arrivals = refined
3597 arrivals.sort(key=lambda x: (x.x, x.t))
3598 return arrivals
3600 @classmethod
3601 def from_scanlines(cls, producer):
3602 '''
3603 Create layer cake model from sequence of materials at depths.
3605 :param producer: iterable yielding (depth, material, name) tuples
3607 Creates a new :py:class:`LayeredModel` object and uses its
3608 :py:meth:`append` method to add layers and discontinuities as needed.
3609 '''
3611 self = cls()
3612 for z, material, name in producer:
3614 if not self._elements:
3615 self.append(Surface(z, material))
3616 else:
3617 element = self._elements[-1]
3618 if self.zeq(element.zbot, z):
3619 assert isinstance(element, Layer)
3620 self.append(
3621 Interface(z, element.mbot, material, name=name))
3623 else:
3624 if isinstance(element, Discontinuity):
3625 ztop = element.z
3626 mtop = element.mbelow
3627 elif isinstance(element, Layer):
3628 ztop = element.zbot
3629 mtop = element.mbot
3631 if mtop == material:
3632 layer = HomogeneousLayer(
3633 ztop, z, material, name=name)
3634 else:
3635 layer = GradientLayer(
3636 ztop, z, mtop, material, name=name)
3638 self.append(layer)
3640 return self
3642 def to_scanlines(self, get_burgers=False):
3643 def fmt(z, m):
3644 if not m._has_default_burgers() or get_burgers:
3645 return (z, m.vp, m.vs, m.rho, m.qp, m.qs,
3646 m.burger_eta1, m.burger_eta2, m.burger_valpha)
3647 return (z, m.vp, m.vs, m.rho, m.qp, m.qs)
3649 last = None
3650 lines = []
3651 for element in self.elements():
3652 if isinstance(element, Layer):
3653 if not isinstance(last, Layer):
3654 lines.append(fmt(element.ztop, element.mtop))
3656 lines.append(fmt(element.zbot, element.mbot))
3658 last = element
3660 if not isinstance(last, Layer):
3661 lines.append(fmt(last.z, last.mbelow))
3663 return lines
3665 def iter_material_parameter(self, get):
3666 assert get in ('vp', 'vs', 'rho', 'qp', 'qs', 'z')
3667 if get == 'z':
3668 for layer in self.layers():
3669 yield layer.ztop
3670 yield layer.zbot
3671 else:
3672 getter = operator.attrgetter(get)
3673 for layer in self.layers():
3674 yield getter(layer.mtop)
3675 yield getter(layer.mbot)
3677 def profile(self, get):
3678 '''
3679 Get parameter profile along depth of the earthmodel.
3681 :param get: property to be queried (
3682 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3683 :type get: str
3684 '''
3686 return num.array(list(self.iter_material_parameter(get)))
3688 def min(self, get='vp'):
3689 '''
3690 Find minimum value of a material property or depth.
3692 :param get: property to be queried (
3693 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``)
3694 '''
3696 return min(self.iter_material_parameter(get))
3698 def max(self, get='vp'):
3699 '''
3700 Find maximum value of a material property or depth.
3702 :param get: property to be queried (
3703 ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, ``'qs'``, or ``'z'``)
3704 '''
3706 return max(self.iter_material_parameter(get))
3708 def simplify_layers(self, layers, max_rel_error=0.001):
3709 if len(layers) <= 1:
3710 return layers
3712 ztop = layers[0].ztop
3713 zbot = layers[-1].zbot
3714 zorigs = [layer.ztop for layer in layers]
3715 zorigs.append(zbot)
3716 zs = num.linspace(ztop, zbot, 100)
3717 data = []
3718 for z in zs:
3719 if z == ztop:
3720 direction = UP
3721 else:
3722 direction = DOWN
3724 mat = self.material(z, direction)
3725 data.append(mat.astuple())
3727 data = num.array(data, dtype=float)
3728 data_means = num.mean(data, axis=0)
3729 nmax = len(layers) // 2
3730 accept = False
3732 zcut_best = []
3733 for n in range(1, nmax+1):
3734 ncutintervals = 20
3735 zdelta = (zbot-ztop)/ncutintervals
3736 if n == 2:
3737 zcuts = [
3738 [ztop, ztop + i*zdelta, zbot]
3739 for i in range(1, ncutintervals)]
3740 elif n == 3:
3741 zcuts = []
3742 for j in range(1, ncutintervals):
3743 for i in range(j+1, ncutintervals):
3744 zcuts.append(
3745 [ztop, ztop + j*zdelta, ztop + i*zdelta, zbot])
3746 else:
3747 zcuts = []
3748 zcuts.append(num.linspace(ztop, zbot, n+1))
3749 if zcut_best:
3750 zcuts.append(sorted(num.linspace(
3751 ztop, zbot, n).tolist() + zcut_best[1]))
3752 zcuts.append(sorted(num.linspace(
3753 ztop, zbot, n-1).tolist() + zcut_best[2]))
3755 best = None
3756 for icut, zcut in enumerate(zcuts):
3757 rel_par_errors = num.zeros(5)
3758 mpar_nodes = num.zeros((n+1, 5))
3760 for ipar in range(5):
3761 znodes, vnodes, error_rms = util.polylinefit(
3762 zs, data[:, ipar], zcut)
3764 mpar_nodes[:, ipar] = vnodes
3765 if data_means[ipar] == 0.0:
3766 rel_par_errors[ipar] = -1
3767 else:
3768 rel_par_errors[ipar] = error_rms/data_means[ipar]
3770 rel_error = rel_par_errors.max()
3771 if best is None or rel_error < best[0]:
3772 best = (rel_error, zcut, mpar_nodes)
3774 rel_error, zcut, mpar_nodes = best
3776 zcut_best.append(list(zcut))
3777 zcut_best[-1].pop(0)
3778 zcut_best[-1].pop()
3780 if rel_error <= max_rel_error:
3781 accept = True
3782 break
3784 if not accept:
3785 return layers
3787 rel_error, zcut, mpar_nodes = best
3789 material_nodes = []
3790 for i in range(n+1):
3791 material_nodes.append(Material(*mpar_nodes[i, :]))
3793 out_layers = []
3794 for i in range(n):
3795 mtop = material_nodes[i]
3796 mbot = material_nodes[i+1]
3797 ztop = zcut[i]
3798 zbot = zcut[i+1]
3799 if mtop == mbot:
3800 lyr = HomogeneousLayer(ztop, zbot, mtop)
3801 else:
3802 lyr = GradientLayer(ztop, zbot, mtop, mbot)
3804 out_layers.append(lyr)
3805 return out_layers
3807 def simplify(self, max_rel_error=0.001):
3808 '''
3809 Get representation of model with lower resolution.
3811 Returns an approximation of the model. All discontinuities are kept,
3812 but layer stacks with continuous model parameters are represented, if
3813 possible, by a lower number of layers. Piecewise linear functions are
3814 fitted against the original model parameter's piecewise linear
3815 functions. Successively larger numbers of layers are tried, until the
3816 difference to the original model is below ``max_rel_error``. The
3817 difference is measured as the RMS error of the fit normalized by the
3818 mean of the input (i.e. the fitted curves should deviate, on average,
3819 less than 0.1% from the input curves if ``max_rel_error`` = 0.001).
3820 '''
3822 mod_simple = LayeredModel()
3824 glayers = []
3825 for element in self.elements():
3827 if isinstance(element, Discontinuity):
3828 for layer in self.simplify_layers(
3829 glayers, max_rel_error=max_rel_error):
3831 mod_simple.append(layer)
3833 glayers = []
3834 mod_simple.append(element)
3835 else:
3836 glayers.append(element)
3838 for layer in self.simplify_layers(
3839 glayers, max_rel_error=max_rel_error):
3840 mod_simple.append(layer)
3842 return mod_simple
3844 def extract(self, depth_min=None, depth_max=None):
3845 '''
3846 Extract :py:class:`LayeredModel` from :py:class:`LayeredModel`.
3848 :param depth_min: depth of upper cut or name of :py:class:`Interface`
3849 :param depth_max: depth of lower cut or name of :py:class:`Interface`
3851 Interpolates a :py:class:`GradientLayer` at ``depth_min`` and/or
3852 ``depth_max``.
3853 '''
3855 if isinstance(depth_min, str):
3856 depth_min = self.discontinuity(depth_min).z
3858 if isinstance(depth_max, str):
3859 depth_max = self.discontinuity(depth_max).z
3861 mod_extracted = LayeredModel()
3863 for element in self.elements():
3864 element = element.copy()
3865 do_append = False
3866 if (depth_min is None or depth_min <= element.ztop) \
3867 and (depth_max is None or depth_max >= element.zbot):
3868 mod_extracted.append(element)
3869 continue
3871 if depth_min is not None:
3872 if element.ztop < depth_min and depth_min < element.zbot:
3873 _, element = element.split(depth_min)
3874 do_append = True
3876 if depth_max is not None:
3877 if element.zbot > depth_max and depth_max > element.ztop:
3878 element, _ = element.split(depth_max)
3879 do_append = True
3881 if do_append:
3882 mod_extracted.append(element)
3884 return mod_extracted
3886 def replaced_crust(self, crust2_profile=None, crustmod=None):
3887 if crust2_profile is not None:
3888 profile = anything_to_crust2_profile(crust2_profile)
3889 crustmod = LayeredModel.from_scanlines(
3890 from_crust2x2_profile(profile))
3892 newmod = LayeredModel()
3893 for element in crustmod.extract(depth_max='moho').elements():
3894 if element.name != 'moho':
3895 newmod.append(element)
3896 else:
3897 moho1 = element
3899 mod = self.extract(depth_min='moho')
3900 first = True
3901 for element in mod.elements():
3902 if element.name == 'moho':
3903 if element.z <= moho1.z:
3904 mbelow = mod.material(moho1.z, direction=UP)
3905 else:
3906 mbelow = element.mbelow
3908 moho = Interface(moho1.z, moho1.mabove, mbelow, name='moho')
3909 newmod.append(moho)
3910 else:
3911 if first:
3912 if isinstance(element, Layer) and element.zbot > moho.z:
3913 newmod.append(GradientLayer(
3914 moho.z,
3915 element.zbot,
3916 moho.mbelow,
3917 element.mbot,
3918 name=element.name))
3920 first = False
3921 else:
3922 newmod.append(element)
3923 return newmod
3925 def perturb(self, rstate=None, keep_vp_vs=False, zmax=None, **kwargs):
3926 '''
3927 Create a perturbed variant of the earth model.
3929 Randomly change the thickness and material parameters of the earth
3930 model from a uniform distribution.
3932 :param kwargs: Maximum change fraction (e.g. 0.1) of the parameters.
3933 Name the parameter, prefixed by ``p``. Supported parameters are
3934 ``ph, pvp, pvs, prho, pqs, pqp``.
3935 :type kwargs: dict
3936 :param rstate: Random state to draw from, defaults to ``None``
3937 :type rstate: :class:`numpy.random.RandomState`
3938 :param keep_vp_vs: Keep the Vp/Vs ratio, defaults to False
3939 :type keep_vp_vs: bool
3941 :returns: A new, perturbed earth model
3942 :rtype: :class:`~pyrocko.cake.LayeredModel`
3944 .. code-block :: python
3946 perturbed_model = model.perturb(ph=.1, pvp=.05, prho=.1)
3947 '''
3949 if rstate is None:
3950 rstate = num.random.RandomState()
3952 def factor(k):
3953 try:
3954 pval = kwargs[k]
3955 return rstate.uniform(1.0-pval, 1.0+pval)
3956 except KeyError:
3957 return 1.0
3959 def perturb_material(mat):
3960 mat_new = copy.deepcopy(mat)
3961 for param in kwargs:
3962 if param != 'ph':
3963 value = getattr(mat, param[1:])
3964 setattr(mat_new, param[1:], value*factor(param))
3966 if keep_vp_vs:
3967 mat_new.vs = mat_new.vp * (mat.vs / mat.vp)
3969 return mat_new
3971 model_new = LayeredModel()
3972 elements = list(self.elements())
3973 assert isinstance(elements[0], Surface)
3974 z = elements[0].z
3975 mat = None
3976 for element in elements:
3977 element_new = copy.deepcopy(element)
3978 if isinstance(element_new, Discontinuity):
3979 element_new.z = z
3981 elif isinstance(element_new, Layer):
3982 element_new.ztop = z
3983 if zmax is None or element.zbot < zmax:
3984 z += (element_new.zbot-element_new.ztop) * factor('ph')
3985 else:
3986 z += (element_new.zbot-element_new.ztop)
3987 element_new.zbot = z
3989 if mat is None or element.mtop != mat:
3990 mat = element.mtop
3991 if zmax is None or element.ztop < zmax:
3992 mat_new = perturb_material(mat)
3993 else:
3994 mat_new = copy.deepcopy(mat)
3996 element_new.mtop = mat_new
3998 if element.mbot != mat:
3999 mat = element.mbot
4000 if zmax is None or element.zbot < zmax:
4001 mat_new = perturb_material(mat)
4002 else:
4003 mat_new = copy.deepcopy(mat)
4005 element_new.mbot = mat_new
4007 model_new.append(element_new)
4009 return model_new
4011 def require_homogeneous(self):
4012 elements = list(self.elements())
4014 if len(elements) != 2:
4015 raise LayeredModelError(
4016 'Homogeneous model required but found more than one layer in '
4017 'earthmodel.')
4019 if not isinstance(elements[1], HomogeneousLayer):
4020 raise LayeredModelError(
4021 'Homogeneous model required but element #1 is not of type '
4022 'HomogeneousLayer.')
4024 return elements[1].m
4026 def __str__(self):
4027 return '\n'.join(str(element) for element in self._elements)
4030def read_hyposat_model(fn):
4031 '''
4032 Reader for HYPOSAT earth model files.
4034 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
4036 Interface names are translated as follows: ``'MOHO'`` -> ``'moho'``,
4037 ``'CONR'`` -> ``'conrad'``
4038 '''
4040 with open(fn, 'r') as f:
4041 translate = {'MOHO': 'moho', 'CONR': 'conrad'}
4042 lname = None
4043 for iline, line in enumerate(f):
4044 if iline == 0:
4045 continue
4047 z, vp, vs, name = util.unpack_fixed('f10, f10, f10, a4', line)
4048 if not name:
4049 name = None
4050 material = Material(vp*1000., vs*1000.)
4052 tname = translate.get(lname, lname)
4053 yield z*1000., material, tname
4055 lname = name
4058def read_nd_model(fn):
4059 '''
4060 Reader for TauP style '.nd' (named discontinuity) files.
4062 To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
4064 Interface names are translated as follows: ``'mantle'`` -> ``'moho'``,
4065 ``'outer-core'`` -> ``'cmb'``, ``'inner-core'`` -> ``'icb'``.
4067 The format has been modified to include Burgers materials parameters in
4068 columns 7 (burger_eta1), 8 (burger_eta2) and 9. eta(3).
4069 '''
4070 with open(fn, 'r') as f:
4071 for x in read_nd_model_fh(f):
4072 yield x
4075def read_nd_model_str(s):
4076 f = StringIO(s)
4077 for x in read_nd_model_fh(f):
4078 yield x
4079 f.close()
4082def read_nd_model_fh(f):
4083 translate = {'mantle': 'moho', 'outer-core': 'cmb', 'inner-core': 'icb'}
4084 name = None
4085 for line in f:
4086 toks = line.split()
4087 if len(toks) == 9 or len(toks) == 6 or len(toks) == 4:
4088 z, vp, vs, rho = [float(x) for x in toks[:4]]
4089 qp, qs = None, None
4090 burgers = None
4091 if len(toks) == 6 or len(toks) == 9:
4092 qp, qs = [float(x) for x in toks[4:6]]
4093 if len(toks) == 9:
4094 burgers = \
4095 [float(x) for x in toks[6:]]
4097 material = Material(
4098 vp*1000., vs*1000., rho*1000., qp, qs,
4099 burgers=burgers)
4101 yield z*1000., material, name
4102 name = None
4103 elif len(toks) == 1:
4104 name = translate.get(toks[0], toks[0])
4106 f.close()
4109def from_crust2x2_profile(profile, depthmantle=50000):
4110 from pyrocko.dataset import crust2x2
4112 default_qp_qs = {
4113 'soft sed.': (50., 50.),
4114 'hard sed.': (200., 200.),
4115 'upper crust': (600., 400.),
4116 }
4118 z = 0.
4119 for i in range(8):
4120 dz, vp, vs, rho = profile.get_layer(i)
4121 name = crust2x2.Crust2Profile.layer_names[i]
4122 if name in default_qp_qs:
4123 qp, qs = default_qp_qs[name]
4124 else:
4125 qp, qs = None, None
4127 material = Material(vp, vs, rho, qp, qs)
4128 iname = None
4129 if i == 7:
4130 iname = 'moho'
4131 if dz != 0.0:
4132 yield z, material, iname
4133 if i != 7:
4134 yield z+dz, material, name
4135 else:
4136 yield z+depthmantle, material, name
4138 z += dz
4141def write_nd_model_fh(mod, fh):
4142 def fmt(z, mat):
4143 rstr = ' '.join(
4144 util.gform(x, 4)
4145 for x in (
4146 z/1000.,
4147 mat.vp/1000.,
4148 mat.vs/1000.,
4149 mat.rho/1000.,
4150 mat.qp, mat.qs))
4151 if not mat._has_default_burgers():
4152 rstr += ' '.join(
4153 util.gform(x, 4)
4154 for x in (
4155 mat.burger_eta1,
4156 mat.burger_eta2,
4157 mat.burger_valpha))
4158 return rstr.rstrip() + '\n'
4160 translate = {
4161 'moho': 'mantle',
4162 'cmb': 'outer-core',
4163 'icb': 'inner-core'}
4165 last = None
4166 for element in mod.elements():
4167 if isinstance(element, Interface):
4168 if element.name is not None:
4169 n = translate.get(element.name, element.name)
4170 fh.write('%s\n' % n)
4172 elif isinstance(element, Layer):
4173 if not isinstance(last, Layer):
4174 fh.write(fmt(element.ztop, element.mtop))
4176 fh.write(fmt(element.zbot, element.mbot))
4178 last = element
4180 if not isinstance(last, Layer):
4181 fh.write(fmt(last.z, last.mbelow))
4184def write_nd_model_str(mod):
4185 f = StringIO()
4186 write_nd_model_fh(mod, f)
4187 return f.getvalue()
4190def write_nd_model(mod, fn):
4191 with open(fn, 'w') as f:
4192 write_nd_model_fh(mod, f)
4195def builtin_models():
4196 return sorted([
4197 os.path.splitext(os.path.basename(x))[0]
4198 for x in glob.glob(builtin_model_filename('*'))])
4201def builtin_model_filename(modelname):
4202 return util.data_file(os.path.join('earthmodels', modelname+'.nd'))
4205def load_model(fn='ak135-f-continental.m', format='nd', crust2_profile=None):
4206 '''
4207 Load layered earth model from file.
4209 :param fn: filename
4210 :param format: format
4211 :param crust2_profile: ``(lat, lon)`` or
4212 :py:class:`pyrocko.dataset.crust2x2.Crust2Profile` object, merge model
4213 with crustal profile. If ``fn`` is forced to be ``None`` only the
4214 converted CRUST2.0 profile is returned.
4215 :returns: object of type :py:class:`LayeredModel`
4217 The following formats are currently supported:
4219 ============== ===========================================================
4220 format description
4221 ============== ===========================================================
4222 ``'nd'`` 'named discontinuity' format used by the TauP programs
4223 ``'hyposat'`` format used by the HYPOSAT location program
4224 ============== ===========================================================
4226 The naming of interfaces is translated from the file format's native naming
4227 to Cake's own convention (See :py:func:`read_nd_model` and
4228 :py:func:`read_hyposat_model` for details). Cake likes the following
4229 internal names: ``'conrad'``, ``'moho'``, ``'cmb'`` (core-mantle boundary),
4230 ``'icb'`` (inner core boundary).
4231 '''
4233 if fn is not None:
4234 if format == 'nd':
4235 if not os.path.exists(fn) and fn in builtin_models():
4236 fn = builtin_model_filename(fn)
4237 reader = read_nd_model(fn)
4238 elif format == 'hyposat':
4239 reader = read_hyposat_model(fn)
4240 else:
4241 assert False, 'unsupported model format'
4243 mod = LayeredModel.from_scanlines(reader)
4244 if crust2_profile is not None:
4245 return mod.replaced_crust(crust2_profile)
4247 return mod
4249 else:
4250 assert crust2_profile is not None
4251 profile = anything_to_crust2_profile(crust2_profile)
4252 return LayeredModel.from_scanlines(
4253 from_crust2x2_profile(profile))
4256def castagna_vs_to_vp(vs):
4257 '''
4258 Calculate vp from vs using Castagna's relation.
4260 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4261 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4262 1985]
4264 vp = 1.16 * vs + 1360 [m/s]
4266 :param vs: S-wave velocity [m/s]
4267 :returns: P-wave velocity [m/s]
4268 '''
4270 return vs*1.16 + 1360.0
4273def castagna_vp_to_vs(vp):
4274 '''
4275 Calculate vp from vs using Castagna's relation.
4277 Castagna's relation (the mudrock line) is an empirical relation for vp/vs
4278 for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al.,
4279 1985]
4281 vp = 1.16 * vs + 1360 [m/s]
4283 :param vp: P-wave velocity [m/s]
4284 :returns: S-wave velocity [m/s]
4285 '''
4287 return (vp - 1360.0) / 1.16
4290def evenize(x, y, minsize=10):
4291 if x.size < minsize:
4292 return x
4293 ry = (y.max()-y.min())
4294 if ry == 0:
4295 return x
4296 dx = (x[1:] - x[:-1])/(x.max()-x.min())
4297 dy = (y[1:] + y[:-1])/ry
4299 s = num.zeros(x.size)
4300 s[1:] = num.cumsum(num.sqrt(dy**2 + dx**2))
4301 s2 = num.linspace(0, s[-1], x.size)
4302 x2 = num.interp(s2, s, x)
4303 x2[0] = x[0]
4304 x2[-1] = x[-1]
4305 return x2
4308def filled(v, *args, **kwargs):
4309 '''
4310 Create NumPy array filled with given value.
4312 This works like :py:func:`numpy.ones` but initializes the array with ``v``
4313 instead of ones.
4314 '''
4315 x = num.empty(*args, **kwargs)
4316 x.fill(v)
4317 return x
4320def next_or_none(i):
4321 try:
4322 return next(i)
4323 except StopIteration:
4324 return None
4327def reci_or_none(x):
4328 try:
4329 return 1./x
4330 except ZeroDivisionError:
4331 return None
4334def mult_or_none(a, b):
4335 if a is None or b is None:
4336 return None
4337 return a*b
4340def min_not_none(a, b):
4341 if a is None:
4342 return b
4343 if b is None:
4344 return a
4345 return min(a, b)
4348def xytups(xx, yy):
4349 d = []
4350 for x, y in zip(xx, yy):
4351 if num.isfinite(y):
4352 d.append((x, y))
4353 return d
4356def interp(x, xp, fp, monoton):
4357 if monoton == 1:
4358 return xytups(
4359 x, num.interp(x, xp, fp, left=num.nan, right=num.nan))
4360 elif monoton == -1:
4361 return xytups(
4362 x, num.interp(x, xp[::-1], fp[::-1], left=num.nan, right=num.nan))
4363 else:
4364 fs = []
4365 for xv in x:
4366 indices = num.where(num.logical_or(
4367 num.logical_and(xp[:-1] >= xv, xv > xp[1:]),
4368 num.logical_and(xp[:-1] <= xv, xv < xp[1:])))[0]
4370 for i in indices:
4371 xr = (xv - xp[i])/(xp[i+1]-xp[i])
4372 fv = xr*fp[i] + (1.-xr)*fp[i+1]
4373 fs.append((xv, fv))
4375 return fs
4378def float_or_none(x):
4379 if x is not None:
4380 return float(x)
4383def parstore_float(thelocals, obj, *args):
4384 for k, v in thelocals.items():
4385 if k != 'self' and (not args or k in args):
4386 setattr(obj, k, float_or_none(v))