|
# http://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg----------
This module can be used to e.g. calculate arrival times, ray paths, reflection and transmission coefficients, take-off and incidence angles and geometrical spreading factors for arbitrary seismic phases. Computations are done for a spherical earth, even though the module name may suggests something flat.
The main classes defined in this module are:
* :py:class:`Material` - Defines an isotropic elastic material. * :py:class:`PhaseDef` - Defines a seismic phase arrival / wave propagation history. * :py:class:`Leg` - Continuous propagation in a :py:class:`PhaseDef`. * :py:class:`Knee` - Conversion/reflection in a :py:class:`PhaseDef`. * :py:class:`LayeredModel` - Representation of a layer cake model. * :py:class:`Layer` - A layer in a :py:class:`LayeredModel`.
* :py:class:`HomogeneousLayer` - A homogeneous :py:class:`Layer`. * :py:class:`GradientLayer` - A gradient :py:class:`Layer`.
* :py:class:`Discontinuity` - A discontinuity in a :py:class:`LayeredModel`.
* :py:class:`Interface` - A :py:class:`Discontinuity` between two :py:class:`Layer` instances. * :py:class:`Surface` - The surface :py:class:`Discontinuity` on top of a :py:class:`LayeredModel`.
* :py:class:`RayPath` - A fan of rays running through a common sequence of layers / interfaces. * :py:class:`Ray` - A specific ray with a specific (ray parameter, distance, arrival time) choice. * :py:class:`RayElement` - An element of a :py:class:`RayPath`.
* :py:class:`Straight` - A ray segment representing propagation through one :py:class:`Layer`. * :py:class:`Kink` - An interaction of a ray with a :py:class:`Discontinuity`. '''
'''Isotropic elastic material.
:param vp: P-wave velocity [m/s] :param vs: S-wave velocity [m/s] :param rho: density [kg/m^3] :param qp: P-wave attenuation Qp :param qs: S-wave attenuation Qs :param poisson: Poisson ratio :param lame: tuple with Lame parameter `lambda` and `shear modulus` [Pa] :param qk: bulk attenuation Qk :param qmu: shear attenuation Qmu
:param burgers: Burgers rheology paramerters as `tuple`. `transient viscosity` [Pa], <= 0 means infinite value, `steady-state viscosity` [Pa] and `alpha`, the ratio between the effective and unreleaxed shear modulus, mu1/(mu1 + mu2). :type burgers: tuple
If no velocities and no lame parameters are given, standard crustal values of vp = 5800 m/s and vs = 3200 m/s are used. If no Q values are given, standard crustal values of qp = 1456 and qs = 600 are used. If no Burgers material parameters are given, transient and steady-state viscosities are 0 and alpha=1.
Everything is in SI units (m/s, Pa, kg/m^3) unless explicitly stated.
The main material properties are considered independant and are accessible as attributes (it is allowed to assign to these):
.. py:attribute:: vp, vs, rho, qp, qs
Other material properties are considered dependant and can be queried by instance methods. '''
self, vp=None, vs=None, rho=2600., qp=None, qs=None, poisson=None, lame=None, qk=None, qmu=None, burgers=None):
raise InvalidArguments( 'If vp and vs are given, poisson ratio and lame paramters ' 'should not be given.')
poisson = 0.25
raise InvalidArguments( 'Poisson ratio should not be given, when lame parameters ' 'are given.')
poisson = 0.25
raise InvalidArguments( 'If vp is given, Lame parameters should not be given.')
poisson = 0.25 raise InvalidArguments( 'If vs is given, Lame parameters should not be given.')
else: raise InvalidArguments( 'Invalid combination of input parameters in material ' 'definition.')
raise InvalidArguments( 'if qp or qs are given, qk and qmu should not be given.')
else: self.qp = 1456.
else: self.vs = 600.
else:
self.qp = qk else: else: self.qp = 1.0 / (s/qmu + (1.0-s)/qk) else: raise InvalidArguments( 'Invalid combination of input parameters in material ' 'definition.')
'''Get independant material properties as a tuple.
Returns a tuple with ``(vp, vs, rho, qp, qs)``. '''
'''Get Lame's parameter lambda and shear modulus.'''
'''Get Lame's parameter lambda.
Returned units are [Pa]. '''
'''Get shear modulus.
Returned units are [Pa]. '''
'''Get Poisson's ratio.'''
'''Get bulk modulus.'''
'''Get Young's modulus.'''
'''Get vp/vs ratio.''' return self.vp/self.vs
'''Get shear attenuation coefficient Qmu.'''
'''Get bulk attenuation coefficient Qk.''' return self.qp else: else: return (1.-s)/(1.0/self.qp - s/self.qs)
'''Get Burger parameters.''' return self.burger_eta1, self.burger_eta2, self.burger_valpha
return None
'''Get rayleigh velocity assuming a homogenous halfspace.
Returned units are [m/s].'''
self.burger_eta2 == DEFAULT_BURGERS[1] and \ self.burger_valpha == DEFAULT_BURGERS[2]: return False
'''Get a readable listing of the material properties.''' P wave velocity [km/s] : %12g S wave velocity [km/s] : %12g P/S wave vel. ratio : %12g Lame lambda [GPa] : %12g Lame shear modulus [GPa] : %12g Poisson ratio : %12g Bulk modulus [GPa] : %12g Young's modulus [GPa] : %12g Rayleigh wave vel. [km/s] : %12g Density [g/cm**3] : %12g Qp P-wave attenuation : %12g Qs S-wave attenuation (Qmu) : %12g Qk bulk attenuation : %12g transient viscos., eta1 [GPa] : %12g st.-state viscos., eta2 [GPa] : %12g relaxation: valpha : %12g '''.strip()
self.vp/km, self.vs/km, self.vp/self.vs, self.lame_lambda()*1e-9, self.shear_modulus()*1e-9, self.poisson(), self.bulk()*1e-9, self.youngs()*1e-9, self.rayleigh()/km, self.rho/km, self.qp, self.qs, self.qk(), self.burger_eta1*1e-9, self.burger_eta2*1e-9, self.burger_valpha)
def __str__(self): vp, vs, rho, qp, qs = self.astuple() return '%10g km/s %10g km/s %10g g/cm^3 %10g %10g' % ( vp/km, vs/km, rho/km, qp, qs)
def __repr__(self): return 'Material(vp=%s, vs=%s, rho=%s, qp=%s, qs=%s)' % \ tuple(repr(x) for x in ( self.vp, self.vs, self.rho, self.qp, self.qs))
''' Represents a continuous piece of wave propagation in a phase definition
**Attributes:**
To be considered as read-only.
.. py:attribute:: departure
One of the constants :py:const:`UP` or :py:const:`DOWN` indicating upward or downward departure.
.. py:attribute:: mode
One of the constants :py:const:`P` or :py:const:`S`, indicating the propagation mode.
.. py:attribute:: depthmin
``None``, a number (a depth in [m]) or a string (an interface name), minimum depth.
.. py:attribute:: depthmax
``None``, a number (a depth in [m]) or a string (an interface name), maximum depth.
'''
self.depthmin = depthmin
def __str__(self): def sd(d): if isinstance(d, float): return '%g km' % (d/km) else: return 'interface %s' % d
s = '%s mode propagation, departing %s' % ( smode(self.mode).upper(), { UP: 'upward', DOWN: 'downward'}[self.departure])
sc = [] if self.depthmax is not None: sc.append('deeper than %s' % sd(self.depthmax)) if self.depthmin is not None: sc.append('shallower than %s' % sd(self.depthmin))
if sc: s = s + ' (may not propagate %s)' % ' or '.join(sc)
return s
'''Represents a change in wave propagation within a :py:class:`PhaseDef`.
**Attributes:**
To be considered as read-only.
.. py:attribute:: depth
Depth at which the conversion/reflection should happen. this can be a string or a number.
.. py:attribute:: direction
One of the constants :py:const:`UP` or :py:const:`DOWN` to indicate the incoming direction.
.. py:attribute:: in_mode
One of the constants :py:const:`P` or :py:const:`S` to indicate the type of mode of the incoming wave.
.. py:attribute:: out_mode
One of the constants :py:const:`P` or :py:const:`S` to indicate the type of mode of the outgoing wave.
.. py:attribute:: conversion
Boolean, whether there is a mode conversion involved.
.. py:attribute:: reflection
Boolean, whether there is a reflection involved.
.. py:attribute:: headwave
Boolean, whether there is headwave propagation involved.
'''
depth='surface', direction=UP, conversion=True, reflection=False, headwave=False, in_setup_state=True)
depth='surface', direction=UP, conversion=False, reflection=True, headwave=False, in_setup_state=True)
self.out_mode) = args
else:
raise InvalidKneeDef('%s has already been set' % k) else:
raise AttributeError(k)
(self.direction == UP) == self.reflection):
raise InvalidKneeDef( 'cannot enter %s from %s and emit ray upwards' % ( ['conversion', 'reflection'][self.reflection], {UP: 'below', DOWN: 'above'}[self.direction]))
(self.direction == DOWN) == self.reflection):
raise InvalidKneeDef( 'cannot enter %s from %s and emit ray downwards' % ( ['conversion', 'reflection'][self.reflection], {UP: 'below', DOWN: 'above'}[self.direction]))
''' Check whether it is relevant to a given combination of interface, propagation mode, and direction. '''
else:
'''Get outgoing direction.
Returns one of the constants :py:const:`UP` or :py:const:`DOWN`. '''
else:
def __str__(self): x = [] if self.reflection: if self.at_surface(): x.append('surface') else: if not self.headwave: if self.direction == UP: x.append('underside') else: x.append('upperside')
if self.headwave: x.append('headwave propagation along') elif self.reflection and self.conversion: x.append('reflection with conversion from %s to %s' % ( smode(self.in_mode).upper(), smode(self.out_mode).upper())) if not self.at_surface(): x.append('at') elif self.reflection: x.append('reflection') if not self.at_surface(): x.append('at') elif self.conversion: x.append('conversion from %s to %s at' % ( smode(self.in_mode).upper(), smode(self.out_mode).upper())) else: x.append('passing through')
if isinstance(self.depth, float): x.append('interface in %g km depth' % (self.depth/1000.)) else: if not self.at_surface(): x.append('%s' % self.depth)
if not self.reflection: if self.direction == UP: x.append('on upgoing path') else: x.append('on downgoing path')
return ' '.join(x)
else: Knee.__init__(self)
def __str__(self): x = ['propagation as headwave'] if isinstance(self.depth, float): x.append('at interface in %g km depth' % (self.depth/1000.)) else: x.append('at %s' % self.depth)
return ' '.join(x)
self.phasename = phasename
def __str__(self): return 'Unknown classic phase name: %s' % self.phasename
''' Exception raised when an error occures during parsing of a phase definition string. '''
self.definition = definition self.position = position self.exception = exception
def __str__(self): return 'Invalid phase definition: "%s" (at character %i: %s)' % ( self.definition, self.position+1, str(self.exception))
'''Definition of a seismic phase arrival, based on ray propagation path.
:param definition: string representation of the phase in Cake's phase syntax
Seismic phases are conventionally named e.g. P, Pn, PP, PcP, etc. In Cake, a slightly different terminology is adapted, which allows to specify arbitrary conversion/reflection histories for seismic ray paths. The conventions used here are inspired by those used in the TauP toolkit, but are not completely compatible with those.
The definition of a seismic ray propagation path in Cake's phase syntax is a string consisting of an alternating sequence of *legs* and *knees*.
A *leg* represents seismic wave propagation without any conversions, encountering only super-critical reflections. Legs are denoted by ``P``, ``p``, ``S``, or ``s``. The capital letters are used when the take-off of the *leg* is in downward direction, while the lower case letters indicate a take-off in upward direction.
A *knee* is an interaction with an interface. It can be a mode conversion, a reflection, or propagation as a headwave or diffracted wave.
* conversion is simply denoted as: ``(INTERFACE)`` or ``DEPTH`` * upperside reflection: ``v(INTERFACE)`` or ``vDEPTH`` * underside reflection: ``^(INTERFACE)`` or ``^DEPTH`` * normal kind headwave or diffracted wave: ``v_(INTERFACE)`` or ``v_DEPTH``
The interface may be given by name or by depth: INTERFACE is the name of an interface defined in the model, DEPTH is the depth of an interface in [km] (the interface closest to that depth is chosen). If two legs appear consecutively without an explicit *knee*, surface interaction is assumed.
The phase definition may end with a backslash ``\\``, to indicate that the ray should arrive at the receiver from above instead of from below. It is possible to restrict the maximum and minimum depth of a *leg* by appending ``<(INTERFACE)`` or ``<DEPTH`` or ``>(INTERFACE)`` or ``>DEPTH`` after the leg character, respectively.
**Examples:**
* ``P`` - like the classical P, but includes PKP, PKIKP, Pg * ``P<(moho)`` - like classical Pg, but must leave source downwards * ``pP`` - leaves source upward, reflects at surface, then travels as P * ``P(moho)s`` - conversion from P to S at the Moho on upgoing path * ``P(moho)S`` - conversion from P to S at the Moho on downgoing path * ``Pv12p`` - P with reflection at 12 km deep interface (or the interface closest to that) * ``Pv_(moho)p`` - classical Pn * ``Pv_(cmb)p`` - classical Pdiff * ``P^(conrad)P`` - underside reflection of P at the Conrad discontinuity
**Usage:**
>>> from pyrocko.cake import PhaseDef # must escape the backslash >>> my_crazy_phase = PhaseDef('pPv(moho)sP\\\\') >>> print my_crazy_phase Phase definition "pPv(moho)sP\": - P mode propagation, departing upward - surface reflection - P mode propagation, departing downward - upperside reflection with conversion from P to S at moho - S mode propagation, departing upward - surface reflection with conversion from S to P - P mode propagation, departing downward - arriving at target from above
.. note::
(1) These conventions might be extended in a way to allow to fix wave propagation to SH mode, possibly by specifying SH, or a single character (e.g. H) instead of S. This would be benificial for the selection of conversion and reflection coefficients, which currently only deal with the P-SV case. '''
def classic_definitions(): # PmP, PmS, PcP, PcS, SmP, ... '%sv(%s)%s' % (a, {'m': 'moho', 'c': 'cmb'}[r], b.lower())]
# Pg, P, S, Sg a, a.lower())]
# PP, SS, PS, SP, PPP, ...
# Pc, Pdiff, Sc, ...
# depth phases
def classic(phasename): '''Get phase definitions based on classic phase name.
:param phasename: classic name of a phase :returns: list of PhaseDef objects
This returns a list of PhaseDef objects, because some classic phases (like e.g. Pg) can only be represented by two Cake style PhaseDef objects (one with downgoing and one with upgoing first leg). '''
raise UnknownClassicPhase(phasename)
need_leg = True state = 1 sdepth += c continue
knee.depth = float(sdepth)*1000. state = 0
else:
sdepthlim += c continue else: depthlim = float(sdepthlim)*1000. if depthlimtype == '<': depthmax = depthlim else: depthmin = depthlim state = 0
else: depthmin = depthlim else:
else:
in_leg.set_depthmin(depthmin) depthmin = None
else:
else: knee.direction = DOWN
need_leg = True knee.direction = UP knee.reflection = True continue
else: raise PhaseDefParseError( definition, ic, 'invalid character: "%s"' % c)
depthlim = float(sdepthlim)*1000. if depthlimtype == '<': depthmax = depthlim else: depthmin = depthlim state = 0
except (ValueError, InvalidKneeDef) as e: raise PhaseDefParseError(definition, ic, e)
raise PhaseDefParseError( definition, ic, 'unfinished expression')
events[-1].set_depthmin(depthmin)
'''Get the first leg in phase definition.'''
'''Get the last leg in phase definition.'''
''' Iterate over the continuous pieces of wave propagation (legs) defined within this phase definition. '''
''' Iterate over conversions and reflections (knees) defined within this phase definition. '''
'''Get original definition of the phase.'''
''' Get entered classic name if any, or original definition of the phase. '''
return self._classicname else:
'''Translate into textual representation (cake phase syntax).''' else:
else:
x.append('>'+strdepth(el.depthmin))
else: x.append('_')
x.append('\\')
def __repr__(self): if self._definition is not None: return "PhaseDef('%s')" % self._definition else: return "PhaseDef('%s')" % self.used_repr()
def __str__(self): orig = '' used = self.used_repr() if self._definition != used: orig = ' (entered as "%s")' % self._definition
sarrive = '\n - arriving at target from %s' % ('below', 'above')[ self._direction_stop == DOWN]
return 'Phase definition "%s"%s:\n - ' % (used, orig) + \ '\n - '.join(str(ev) for ev in self) + sarrive
'''Get a deep copy of it.'''
else: raise PhaseDefParseError('invalid phase definition')
''' Get indices to select the appropriate element from scatter matrix for free surface. '''
'''Scatter matrix for free surface reflection/conversions.
:param material: material, object of type :py:class:`Material` :param p: flat ray parameter [s/m] :param energy: bool, when ``True`` energy normalized coefficients are returned :returns: Scatter matrix
The scatter matrix is ordered as follows::
[[PP, PS], [SP, SS]]
The formulas given in Aki & Richards are used. '''
scatter = num.array([[-1.0, 0.0], [0.0, 1.0]])
else:
[- vsp_term**2 + pcc_term, 4.0*p*coslam/vp*vsp_term], [4.0*p*cosphi/vs*vsp_term, vsp_term**2 - pcc_term]], dtype=num.complex) / denom
return scatter else: (normvec[:, num.newaxis]) / (normvec[num.newaxis, :]))
''' Get indices to select the appropriate element from scatter matrix for solid-solid interface. '''
(out_direction == DOWN)*2 + (out_mode == S), (in_direction == UP)*2 + (in_mode == S))
'''Scatter matrix for solid-solid interface.
:param material1: material above, object of type :py:class:`Material` :param material2: material below, object of type :py:class:`Material` :param p: flat ray parameter [s/m] :param energy: bool, when ``True`` energy normalized coefficients are returned :returns: Scatter matrix
The scatter matrix is ordered as follows::
[[P1P1, S1P1, P2P1, S2P1], [P1S1, S1S1, P2S1, S2S1], [P1P2, S1P2, P2P2, S2P2], [P1S2, S1S2, P2S2, S2S2]]
The formulas given in Aki & Richards are used. '''
# from aki and richards [-vp1*p, -coslam1, vp2*p, coslam2], [cosphi1, -vs1*p, cosphi2, -vs2*p], [2.0*rho1*vs1**2*p*cosphi1, rho1*vs1*(1.0-2.0*vs1**2*p**2), 2.0*rho2*vs2**2*p*cosphi2, rho2*vs2*(1.0-2.0*vs2**2*p**2)], [-rho1*vp1*(1.0-2.0*vs1**2*p**2), 2.0*rho1*vs1**2*p*coslam1, rho2*vp2*(1.0-2.0*vs2**2*p**2), -2.0*rho2*vs2**2*p*coslam2]], dtype=num.complex)
return scatter else: vp1*rho1*(cosphi1+eps), vs1*rho1*(coslam1+eps), vp2*rho2*(cosphi2+eps), vs2*rho2*(coslam2+eps)], dtype=num.complex) normvec[:, num.newaxis] / normvec[num.newaxis, :])
else: raise BadPotIntCoefs()
elif i == S: return 's'
def __str__(self): return 'Cannot enter layer %i from %s' % ( self._ilayer, { UP: 'below', DOWN: 'above'}[self._direction])
'''Representation of a layer in a layered earth model.
:param ztop: depth of top of layer :param zbot: depth of bottom of layer :param name: name of layer (optional)
Subclasses are: :py:class:`HomogeneousLayer` and :py:class:`GradientLayer`. '''
self.mbot.vp, self.mtop.vp, radius(self.zbot), radius(self.ztop))
self.mbot.vs, self.mtop.vs, radius(self.zbot), radius(self.ztop))
'''Get coefficients for potential interpolation.
:param mode: mode of wave propagation, :py:const:`P` or :py:const:`S` :returns: coefficients ``(a, b)`` '''
else:
''' Tolerantly check if a given depth is within the layer (including boundaries). '''
self.at_bottom(z) or self.at_top(z)
''' Tolerantly check if a given depth is within the layer (not including boundaries). '''
return self.ztop <= z <= self.zbot and not \ self.at_bottom(z) and not \ self.at_top(z)
'''Tolerantly check if given depth is at the bottom of the layer.'''
'''Tolerantly check if given depth is at the top of the layer.'''
''' Convert spherical ray parameter to local flat ray parameter for top of layer. ''' return p / (earthradius-self.ztop)
''' Convert spherical ray parameter to local flat ray parameter for bottom of layer. ''' return p / (earthradius-self.zbot)
''' Convert spherical ray parameter to local flat ray parameter for given depth. '''
''' Get travel time and distance for for traversal with given mode and ray parameter.
:param p: ray parameter (spherical) :param mode: mode of propagation (:py:const:`P` or :py:const:`S`) :param zpart: if given, tuple with two depths to restrict computation to a part of the layer
This implementation uses analytic formulas valid for a spherical earth in the case where the velocity c within the layer is given by potential interpolation of the form
c(z) = a*z^b '''
else: lr = math.log(r2/r1) sap = num.sqrt(1.0/a**2 - p**2) x = p/sap * lr t = 1./(a**2 * sap)
''' Check if wave mode can exist for given ray parameter at given depth within the layer. '''
(utop * radius(self.ztop) - p) > 0., (ubot * radius(self.zbot) - p) > 0.)
'''Get turning depth for given ray parameter and propagation mode.'''
'''Propagate ray through layer.
:param p: ray parameter :param mode: propagation mode :param direction: in direction (:py:const:`UP` or :py:const:`DOWN`''' else:
raise CannotPropagate(direction, self.ilayer)
else:
'''Change layer thinkness and interpolate :py:class:`Material` if required.'''
'''Representation of a homogeneous layer in a layered earth model.
Base class: :py:class:`Layer`. '''
if mode == P: v = self.m.vp if mode == S: v = self.m.vs
if num.isscalar(z): return v else: return filled(v, len(z))
v = self.v(mode) return v, v
u = self.u(mode) pflat = self.pflat_bottom(p) if zpart is None: dz = (self.zbot - self.ztop) else: dz = abs(zpart[1]-zpart[0])
u = self.u(mode) eps = u*0.001 denom = num.sqrt(u**2 - pflat**2) + eps
x = r2d*pflat/(earthradius-self.zmid) * dz / denom t = u**2 * dz / denom return x, t
raise DoesNotTurn()
def __str__(self): if self.name: name = self.name + ' ' else: name = ''
calcmode = ''.join('HP'[self._use_potential_interpolation[mode]] for mode in (P, S))
return ' (%i) homogeneous layer %s(%g km - %g km) [%s]\n %s' % ( self.ilayer, name, self.ztop/km, self.zbot/km, calcmode, self.m)
'''Representation of a gradient layer in a layered earth model.
Base class: :py:class:`Layer`. '''
self.interpolate(z, ptop, pbot) for (ptop, pbot) in zip(dtop, dbot)]
if mode == P: return 1./self.interpolate(z, self.mtop.vp, self.mbot.vp) if mode == S: return 1./self.interpolate(z, self.mtop.vs, self.mbot.vs)
if mode == P: return self.mtop.vp, self.mbot.vp if mode == S: return self.mtop.vs, self.mbot.vs
if mode == P: return self.interpolate(z, self.mtop.vp, self.mbot.vp) if mode == S: return self.interpolate(z, self.mtop.vs, self.mbot.vs)
utop, ubot = self.u_top_bottom(mode) b = (1./ubot - 1./utop)/(self.zbot - self.ztop)
pflat = self.pflat_bottom(p) if zpart is not None: utop = self.u(mode, zpart[0]) ubot = self.u(mode, zpart[1])
peps = 1e-16 pdp = pflat + peps
def func(u): eta = num.sqrt(num.maximum(u**2 - pflat**2, 0.0)) xx = eta/u tt = num.where( pflat <= u, num.log(u+eta) - num.log(pdp) - eta/u, 0.0)
return xx, tt
xxtop, tttop = func(utop) xxbot, ttbot = func(ubot)
x = (xxtop - xxbot) / (b*pdp) t = (tttop - ttbot) / b + pflat*x
x *= r2d/(earthradius - self.zmid) return x, t
pflat = self.pflat_bottom(p) vtop, vbot = self.v_top_bottom(mode) return (1./pflat - vtop) * (self.zbot - self.ztop) / \ (vbot-vtop) + self.ztop
def __str__(self): if self.name: name = self.name + ' ' else: name = ''
calcmode = ''.join('HP'[self._use_potential_interpolation[mode]] for mode in (P, S))
return ''' (%i) gradient layer %s(%g km - %g km) [%s] %s %s''' % ( self.ilayer, name, self.ztop/km, self.zbot/km, calcmode, self.mtop, self.mbot)
'''Base class for discontinuities in layered earth model.
Subclasses are: :py:class:`Interface` and :py:class:`Surface`. '''
self.z = z self.zbot = z self.ztop = z
'''Representation of an interface in a layered earth model.
Base class: :py:class:`Discontinuity`. '''
def __str__(self): if self.name is None: return 'interface' else: return 'interface "%s"' % self.name
mult_or_none(uabove, radius(self.z)), mult_or_none(ubelow, radius(self.z)))
else: else: return -direction
self.mabove, self.mbelow, self.pflat(p), energy=True) psv_solid_ind(in_direction, out_direction, in_mode, out_mode)]
'''Representation of the surface discontinuity in a layered earth model.
Base class: :py:class:`Discontinuity`. '''
else: self.mbelow, self.pflat(p), energy=True)[ psv_surface_ind(in_mode, out_mode)]
def __str__(self): return 'surface'
else:
else: raise BottomReached()
else: raise SurfaceReached()
'''An element of a :py:class:`RayPath`.'''
''' A ray segment representing wave propagation through one :py:class:`Layer`. '''
else:
else:
return p / (earthradius-self.z_out(endgaps))
return self.layer.test(p, self.mode, z)
else:
else:
else:
return self._direction_out else:
return self.layer.u(self.mode, self.z_out(endgaps))
else:
self._direction_in, self._direction_out, self.mode, id(self.layer)))
return self.interface.z
return self.interface.z
return filled(self.interface.z, len(p))
'''An interaction of a ray with a :py:class:`Discontinuity`.'''
self, in_direction, out_direction, in_mode, out_mode, discontinuity):
self.in_direction, out_direction, self.in_mode, out_mode, p)
def __str__(self): r, c = self.reflection(), self.conversion() if r and c: return '|~' if r: return '|' if c: return '~' return '_'
self.in_direction, self.out_direction, self.in_mode, self.out_mode, id(self.discontinuity)))
''' Representation of a fan of rays running through a common sequence of layers / interfaces. '''
'''Get a copy of it.'''
c = copy.copy(self) c.elements = list(self.elements) return c
'''Get information needed for end point adjustments.'''
zstart, zstop, self.phase.direction_start(), self.phase.direction_stop())
raise PRangeNotSet()
'''Calculate phase definition from ray path.'''
n_elements_n[:-2], n_elements_n[1:-1], n_elements_n[2:]):
type(before), type(after)):
z, element.in_direction, element.out_direction != element.in_direction, element.in_mode, element.out_mode))
z = element.interface.z k = Knee( z, before.in_direction, after.out_direction != before.in_direction, before.in_mode, after.out_mode)
k.headwave = True used.append(k) used.append(Leg(after.out_direction, after.out_mode))
and element.is_straight() and before.is_kink() and after.is_kink() and element.turns() and not before.reflection() and not before.conversion() and not after.reflection() and not after.conversion()):
Head( before.discontinuity.z, before.out_direction, element.mode)) Leg(-before.out_direction, element.mode))
'''Get maximum valid ray parameter.''' self._check_have_prange() return self._pmax
'''Get minimum valid ray parameter.'''
'''Get minimal distance.''' self._analyse() return self._xmin
'''Get maximal distance.''' self._analyse() return self._xmax
''' Iterate over propagation mode changes (reflections/transmissions). '''
'''Iterate over ray segments.'''
'''Get first ray segment.'''
'''Get last ray segment.'''
''' Get product of all conversion/reflection coefficients encountered on path. ''' operator.mul, (k.efficiency(p) for k in self.kinks()), 1.)
'''Get geometrical spreading factor.''' return 0.0
return num.nan
x = x1 p = dp
4.0 * math.pi * num.sin(x) * (earthradius-first.z_in(endgaps)) * (earthradius-last.z_out(endgaps))**2 * first.u_in(endgaps)**2 * num.abs(num.cos(first.angle_in(p, endgaps)*d2r)) * num.abs(num.cos(last.angle_out(p, endgaps)*d2r)))
''' Get amount of distance/traveltime to be subtracted at the generic ray path's ends. '''
p, zstart, firsts.z_in(), dirstart == firsts._direction_in) p, zstop, lasts.z_out(), dirstop == lasts._direction_out)
'''Check if ray parameter is valid at source and receiver.'''
zstart, zstop, dirstart, dirstop = endgaps firsts = self.first_straight() lasts = self.last_straight() return num.logical_and(firsts.test(p, zstart), lasts.test(p, zstop))
'''Calculate distance and traveltime for given ray parameter.'''
sx = num.zeros(p.size) st = num.zeros(p.size) else:
''' Calculate limits of distance and traveltime for given ray parameter. '''
else: sx = 0.0 st = 0.0 sxe = 0.0 ste = 0.0
''' Iterate over (depth, distance, traveltime) at each layer interface on ray path. '''
sx = num.zeros(p.size) st = num.zeros(p.size) ok = False for s in self.straights(): yield s.z_in(), sx.copy(), st.copy()
x, t = s.xt(p) sx += x st += t ok = True
if ok: yield s.z_out(), sx.copy(), st.copy()
self, p, endgaps, points_per_straight=20, x_for_headwave=None):
'''Get geometrical representation of ray path.'''
assert p.size == 1 x, t = self.xt(p, endgaps) xstretch = x_for_headwave-x nout = xstretch.size else:
# first create full path including the endgaps
z = zin for i in range(n): xs = float(i)/(n-1) * xstretch ts = s.x2t_headwave(xs) zxt.append((filled(z, xstretch.size), sx+xs, st+ts)) else:
else: # ray turns in layer float(i)/(n-1)*math.pi/2.0) * 0.999
else: x, t = s.xt(p, zpart=[z, zin])
x = xstretch t = s.x2t_headwave(xstretch) else:
# gather results as arrays with such that x[ip, ipoint]
# cut off the endgaps, add exact endpoints
self.critical_pstart(endgaps), self.critical_pstop(endgaps))
empty = num.array([], dtype=num.float) return empty, empty, empty
else:
else: dx, dt = self.xt_endgaps(self._p, endgaps) p, x, t = self._p, self._x - dx, self._t - dt p, x, t = p[0], x[0], t[0] xh = num.linspace(0., x*10-x, 10) th = self.headwave_straight().x2t_headwave(xh) return filled(p, xh.size), x+xh, t+th
'''Get approximate ray parameter and traveltime for distance.'''
(x_, self._p[0], t, None) for (x_, t) in zip(xok, th)]
else:
len(xp) == 0 and rx[0] == 0.0 and any(x == 0.0) and rp[0] == 0.0):
(x_, p, t, (rp, rx, rt)) for ((x_, p), (_, t)) in zip(xp, xt)]
return False
tuple(hash(x) for x in self.elements) + (self.phase.definition(), ))
def __str__(self, p=None, eps=1.): x = [] start_i = None end_i = None turn_i = None
def append_layers(si, ei, ti): if si == ei and (ti is None or ti == si): x.append('%i' % si) else: if ti is not None: x.append('(%i-%i-%i)' % (si, ti, ei)) else: x.append('(%i-%i)' % (si, ei))
for el in self.elements: if type(el) is Straight: if start_i is None: start_i = el.layer.ilayer if el._direction_in != el._direction_out: turn_i = el.layer.ilayer end_i = el.layer.ilayer
elif isinstance(el, Kink): if start_i is not None: append_layers(start_i, end_i, turn_i) start_i = None turn_i = None
x.append(str(el))
if start_i is not None: append_layers(start_i, end_i, turn_i)
su = '(%s)' % self.used_phase(p=p, eps=eps).used_repr()
return '%-15s %-17s %s' % (self.phase.definition(), su, ''.join(x))
'''Get critical ray parameter for source depth choice.'''
'''Get critical ray parameter for receiver depth choice.'''
'''Get valid ranges of ray parameter, distance, and traveltime.'''
'''Get textual representation.'''
xunit = 'deg' xfact = 1. else:
- x [%g, %g] %s - t [%g, %g] s - p [%g, %g] s/deg ''' % ( self._xmin*xfact, self._xmax*xfact, xunit, self._tmin, self._tmax, self._pmin/r2d, self._pmax/r2d)
\n - x [%g, %g] %s \n - t [%g, %g] s \n - p [%g, %g] s/deg \n''' % (xmin*xfact, xmax*xfact, xunit, tmin, tmax, pmin/r2d, pmax/r2d)
else: ss = ''
''' Representation of a ray with a specific (path, ray parameter, distance, arrival time) choice.
**Attributes:**
.. py:attribute:: path
:py:class:`RayPath` object containing complete propagation history.
.. py:attribute:: p
Ray parameter (spherical) [s/rad]
.. py:attribute:: x
Radial distance [deg]
.. py:attribute:: t
Traveltime [s]
.. py:attribute:: endgaps
Needed for source/receiver depth adjustments in many :py:class:`RayPath` methods. '''
'''Get phase definition which was used to create the ray.
:returns: :py:class:`PhaseDef` object '''
return self.path.phase
'''Compute phase definition from propagation path.
:returns: :py:class:`PhaseDef` object '''
raise RefineFailed()
except ValueError: raise RefineFailed()
'''Get takeoff angle of ray.
The angle is returned in [degrees]. '''
'''Get incidence angle of ray.
The angle is returned in [degrees]. '''
'''Get conversion/reflection efficiency of the ray.
A value between 0 and 1 is returned, reflecting the relative amount of energy which is transmitted along the ray and not lost by reflections or conversions. '''
'''Get geometrical spreading factor.'''
'''Get geometrical representation of ray path.
Three arrays (depth, distance, time) with points on the ray's path of propagation are returned. The number of points which are used in each ray segment (passage through one layer) may be controlled by the ``points_per_straight`` parameter. ''' num.atleast_1d(self.p), self.endgaps, points_per_straight=points_per_straight, x_for_headwave=num.atleast_1d(self.x))
def __str__(self, as_degrees=False): if as_degrees: sd = '%6.3g deg' % self.x else: sd = '%7.5g km' % (self.x*(d2r*earthradius/km))
return '%7.5g s/deg %s %6.4g s %5.1f %5.1f %3.0f%% %3.0f%% %s' % ( self.p/r2d, sd, self.t, self.takeoff_angle(), self.incidence_angle(), 100*self.efficiency(), 100*self.spreading()*self.surface_sphere(), self.path.__str__(p=self.p))
elif isinstance(crust2_profile, (str, newstr)): return crust2x2.get_profile(crust2_profile) elif isinstance(crust2_profile, crust2x2.Crust2Profile): return crust2_profile else: assert False, 'crust2_profile must be (lat, lon) a profile ' \ 'key or a crust2x2 Profile object)'
CakeError.__init__(self) self.depth_or_name = depth_or_name
def __str__(self): return 'Cannot find discontinuity from given depth or name: %s' % \ self.depth_or_name
'''Representation of a layer cake model.
There are several ways to initialize an instance of this class.
1. Use the module function :py:func:`load_model` to read a model from a file. 2. Create an empty model with the default constructor and append layers and discontinuities with the :py:meth:`append` method (from top to bottom). 3. Use the constructor :py:meth:`LayeredModel.from_scanlines`, to automatically create the :py:class:`Layer` and :py:class:`Discontinuity` objects from a given velocity profile.
An earth model is represented by as stack of :py:class:`Layer` and :py:class:`Discontinuity` objects. The method :py:meth:`arrivals` returns :py:class:`Ray` objects which may be e.g. queried for arrival times of specific phases. Each ray is associated with a :py:class:`RayPath` object. Ray objects share common ray paths if they have the same conversion/reflection/propagation history. Creating the ray path objects is relatively expensive (this is done in :py:meth:`gather_paths`), but they are cached for reuse in successive invocations. '''
'''Get a copy of the model with surface layer stretched to given elevation.
:param elevation: new surface elevation in [m]
Elevation is positiv upward, contrary to the layered models downward `z` axis. '''
c = copy.deepcopy(self) c._pathcache = {} surface = c._elements[0] toplayer = c._elements[1]
assert toplayer.zbot > -elevation
surface.z = -elevation c._elements[1] = toplayer.copy(ztop=-elevation) c._elements[1].ilayer = 0 return c
'''Add a layer or discontinuity at bottom of model.
:param element: object of subclass of :py:class:`Layer` or :py:class:`Discontinuity`. '''
raise CakeError('Layer deeper than earthradius')
'''Iterate over all elements of the model.
:param direction: direction of traversal :py:const:`DOWN` or :py:const:`UP`.
Objects derived from the :py:class:`Discontinuity` and :py:class:`Layer` classes are yielded. '''
else: return reversed(self._elements)
'''Iterate over all layers of model.
:param direction: direction of traversal :py:const:`DOWN` or :py:const:`UP`.
Objects derived from the :py:class:`Layer` class are yielded. '''
else: el for el in reversed(self._elements) if isinstance(el, Layer))
'''Get layer for given depth.
:param z: depth [m] :param direction: direction of traversal :py:const:`DOWN` or :py:const:`UP`.
Returns first layer which touches depth ``z`` (tolerant at boundaries). '''
else: raise CakeError('Failed extracting layer at depth z=%s' % z)
'''Get material at given depth.
:param z: depth [m] :param direction: direction of traversal :py:const:`DOWN` or :py:const:`UP` :returns: object of type :py:class:`Material`
If given depth ``z`` happens to be at an interface, the material of the first layer with respect to the the traversal ordering is returned. '''
'''Iterate over all discontinuities of the model.'''
'''Get discontinuity by name or depth.
:param name_or_z: name of discontinuity or depth [m] as float value '''
candi = sorted( self.discontinuities(), key=lambda i: abs(i.z-name_or_z)) else:
raise DiscontinuityNotFound(name_or_z)
'''Adapt a phase definition for use with this model.
This returns a copy of the phase definition, where named discontinuities are replaced with the actual depth of these, as defined in the model. '''
''' Get ray path for given combination of ray parameter, phase definition, source and receiver layers.
:param p: ray parameter (spherical) [s/rad] :param phase: phase definition (:py:class:`PhaseDef` object) :param layer_start: layer with source :param layer_stop: layer with receiver :returns: :py:class:`RayPath` object
If it is not possible to find a solution, an exception of type :py:exc:`NotPhaseConform`, :py:exc:`MinDepthReached`, :py:exc:`MaxDepthReached`, :py:exc:`CannotPropagate`, :py:exc:`BottomReached` or :py:exc:`SurfaceReached` is raised. '''
raise CannotPropagate(direction, current.ilayer)
# detect trapped wave raise Trapped()
current, mode, direction):
else: # implicit reflection/transmission
olddirection, olddirection, oldmode, oldmode, current))
olddirection, direction, oldmode, current))
olddirection, direction, oldmode, mode, current))
else: olddirection, direction, oldmode, mode, current))
raise MinDepthReached() raise MaxDepthReached() else: raise MinDepthReached() raise MaxDepthReached()
current is layer_stop:
else:
''' Get all possible ray paths for given source and receiver depths for one or more phase definitions.
:param phases: a :py:class:`PhaseDef` object or a list of such objects. Comma-separated strings and lists of such strings are also accepted and are converted to :py:class:`PhaseDef` objects for convenience. :param zstart: source depth [m] :param zstop: receiver depth [m] :returns: a list of :py:class:`RayPath` objects
Results of this method are cached internally. Cached results are returned, when a given combination of source layer, receiver layer and phase definition has been used before. '''
else:
# diffracted wave: pbelow is None or pbelow >= pabove):
else: radius(z)/layer_start.v(phase.first_leg().mode, z) for z in (layer_start.ztop, layer_start.zbot)])
radius(z)/layer_stop.v(phase.last_leg().mode, z) for z in (layer_stop.ztop, layer_stop.zbot)])
p, phase, layer_start, layer_stop)
hash(path1) != hash(path2):
min(ps), max(ps), pmax/(self._np-1))
except ZeroDivisionError: phase_paths = []
self, distances=[], phases=PhaseDef('P'), zstart=0.0, zstop=0.0, refine=True):
'''Compute rays and traveltimes for given distances.
:param distances: list or array of distances [deg] :param phases: a :py:class:`PhaseDef` object or a list of such objects. Comma-separated strings and lists of such strings are also accepted and are converted to :py:class:`PhaseDef` objects for convenience. :param zstart: source depth [m] :param zstop: receiver depth [m] :param refine: bool flag, whether to use bisectioning to improve (p, x, t) estimated from interpolation :returns: a list of :py:class:`Ray` objects, sorted by (distance, arrival time) '''
distances, endgaps):
except RefineFailed: pass
def from_scanlines(cls, producer): '''Create layer cake model from sequence of materials at depths.
:param producer: iterable yielding (depth, material, name) tuples
Creates a new :py:class:`LayeredModel` object and uses its :py:meth:`append` method to add layers and discontinuities as needed. '''
else: Interface(z, element.mbot, material, name=name))
else:
ztop, z, material, name=name) else: ztop, z, mtop, material, name=name)
m.burger_eta1, m.burger_eta2, m.burger_valpha)
else:
''' Get parameter profile along depth of the earthmodel.
:param get: property to be queried ( ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``) :type get: string '''
''' Find minimum value of a material property or depth.
:param get: property to be queried ( ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, or ``'qs'``, or ``'z'``) '''
return min(self.iter_material_parameter(get))
''' Find maximum value of a material property or depth.
:param get: property to be queried ( ``'vp'``, ``'vs'``, ``'rho'``, ``'qp'``, ``'qs'``, or ``'z'``) '''
return max(self.iter_material_parameter(get))
else:
[ztop, ztop + i*zdelta, zbot] for i in range(1, ncutintervals)] [ztop, ztop + j*zdelta, ztop + i*zdelta, zbot]) else: ztop, zbot, n).tolist() + zcut_best[1])) ztop, zbot, n-1).tolist() + zcut_best[2]))
zs, data[:, ipar], zcut)
else:
lyr = HomogeneousLayer(ztop, zbot, mtop) else:
'''Get representation of model with lower resolution.
Returns an approximation of the model. All discontinuities are kept, but layer stacks with continuous model parameters are represented, if possible, by a lower number of layers. Piecewise linear functions are fitted against the original model parameter's piecewise linear functions. Successively larger numbers of layers are tried, until the difference to the original model is below ``max_rel_error``. The difference is measured as the RMS error of the fit normalized by the mean of the input (i.e. the fitted curves should deviate, on average, less than 0.1% from the input curves if ``max_rel_error`` = 0.001).'''
glayers, max_rel_error=max_rel_error):
else:
glayers, max_rel_error=max_rel_error):
'''Extract :py:class:`LayeredModel` from :py:class:`LayeredModel`.
:param depth_min: depth of upper cut or name of :py:class:`Interface` :param depth_max: depth of lower cut or name of :py:class:`Interface`
Interpolates a :py:class:`GradientLayer` at ``depth_min`` and/or ``depth_max``.'''
and (depth_max is None or depth_max >= element.zbot):
from_crust2x2_profile(profile))
else:
mbelow = mod.material(moho1.z, direction=UP) else:
else: moho.z, element.zbot, moho.mbelow, element.mbot, name=element.name))
else:
''' Create a perturbed variant of the earth model.
Randomly change the thickness and material parameters of the earth model from a uniform distribution.
:param kwargs: Maximum change fraction (e.g. 0.1) of the parameters. Name the parameter, prefixed by ``p``. Supported parameters are ``ph, pvp, pvs, prho, pqs, pqp``. :type kwargs: dict :param rstate: Random state to draw from, defaults to ``None`` :type rstate: :class:`numpy.random.RandomState`, optional :param keep_vp_vs: Keep the Vp/Vs ratio, defaults to False :type keep_vp_vs: bool, optional
:returns: A new, perturbed earth model :rtype: :class:`~pyrocko.cake.LayeredModel`
.. code-block :: python
perturbed_model = model.perturb(ph=.1, pvp=.05, prho=.1) ''' _pargs = set(['ph', 'pvp', 'pvs', 'prho', 'pqs', 'pqp']) earthmod = copy.deepcopy(self)
if rstate is None: rstate = num.random.RandomState()
layers = earthmod.layers() discont = earthmod.discontinuities() prev_layer = None
def get_change_ratios(): values = dict.fromkeys([p[1:] for p in _pargs], 0.)
for param, pval in kwargs.items(): if param not in _pargs: continue values[param[1:]] = float(rstate.uniform(-pval, pval, size=1)) return values
# skip Surface while True: disc = next(discont) if isinstance(disc, Surface): break
while True: try: layer = next(layers) m = layer.material(None) h = layer.zbot - layer.ztop except StopIteration: break
if not isinstance(layer, HomogeneousLayer): raise NotImplementedError( 'Can only perturbate homogeneous layers!')
changes = get_change_ratios()
# Changing thickness dh = h * changes['h'] changes['h'] = dh
layer.resize(depth_max=layer.zbot + dh, depth_min=prev_layer.zbot if prev_layer else None)
try: disc = next(discont) disc.change_depth(disc.z + dh) except StopIteration: pass
# Setting material parameters for param, change_ratio in changes.items(): if param == 'h': continue
value = m.__getattribute__(param) changes[param] = value * change_ratio
if keep_vp_vs and changes['vp'] != 0.: changes['vs'] = (m.vp + changes['vp']) / m.vp_vs_ratio() - m.vs
for param, change in changes.items(): if param == 'h': continue value = m.__getattribute__(param) m.__setattr__(param, value + change)
logger.info( 'perturbating earthmodel: {}'.format( ' '.join(['{param}: {change:{len}.2f}'.format( param=p, change=c, len=8) for p, c in changes.items()])))
prev_layer = layer
return earthmod
raise LayeredModelError('More than one layer in earthmodel') raise LayeredModelError('Layer has to be a HomogeneousLayer')
def __str__(self): return '\n'.join(str(element) for element in self._elements)
'''Reader for HYPOSAT earth model files.
To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
Interface names are translated as follows: ``'MOHO'`` -> ``'moho'``, ``'CONR'`` -> ``'conrad'`` '''
with open(fn, 'r') as f: translate = {'MOHO': 'moho', 'CONR': 'conrad'} lname = None for iline, line in enumerate(f): if iline == 0: continue
z, vp, vs, name = util.unpack_fixed('f10, f10, f10, a4', line) if not name: name = None material = Material(vp*1000., vs*1000.)
tname = translate.get(lname, lname) yield z*1000., material, tname
lname = name
'''Reader for TauP style '.nd' (named discontinuity) files.
To be used as producer in :py:meth:`LayeredModel.from_scanlines`.
Interface names are translated as follows: ``'mantle'`` -> ``'moho'``, ``'outer-core'`` -> ``'cmb'``, ``'inner-core'`` -> ``'icb'``.
The format has been modified to include Burgers materials parameters in columns 7 (burger_eta1), 8 (burger_eta2) and 9. eta(3). '''
burgers = \ [float(x) for x in toks[6:]]
vp*1000., vs*1000., rho*1000., qp, qs, burgers=burgers)
'soft sed.': (50., 50.), 'hard sed.': (200., 200.), 'upper crust': (600., 400.), }
else:
else:
util.gform(x, 4) for x in ( z/1000., mat.vp/1000., mat.vs/1000., mat.rho/1000., mat.qp, mat.qs)) rstr += ' '.join( util.gform(x, 4) for x in ( mat.burger_eta1, mat.burger_eta2, mat.burger_valpha))
'moho': 'mantle', 'cmb': 'outer-core', 'icb': 'inner-core'}
with open(fn, 'w') as f: write_nd_model_fh(mod, f)
os.path.splitext(os.path.basename(x))[0] for x in glob.glob(builtin_model_filename('*'))])
'''Load layered earth model from file.
:param fn: filename :param format: format :param crust2_profile: ``(lat, lon)`` or :py:class:`pyrocko.crust2x2.Crust2Profile` object, merge model with crustal profile. If ``fn`` is forced to be ``None`` only the converted CRUST2.0 profile is returned. :returns: object of type :py:class:`LayeredModel`
The following formats are currently supported:
============== =========================================================== format description ============== =========================================================== ``'nd'`` 'named discontinuity' format used by the TauP programs ``'hyposat'`` format used by the HYPOSAT location program ============== ===========================================================
The naming of interfaces is translated from the file format's native naming to Cake's own convention (See :py:func:`read_nd_model` and :py:func:`read_hyposat_model` for details). Cake likes the following internal names: ``'conrad'``, ``'moho'``, ``'cmb'`` (core-mantle boundary), ``'icb'`` (inner core boundary). '''
elif format == 'hyposat': reader = read_hyposat_model(fn) else: assert False, 'unsupported model format'
else: assert crust2_profile is not None profile = anything_to_crust2_profile(crust2_profile) return LayeredModel.from_scanlines( from_crust2x2_profile(profile))
'''Calculate vp from vs using castagna's relation.
Castagna's relation (the mudrock line) is an empirical relation for vp/vs for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al., 1985]
vp = 1.16 * vs + 1360 [m/s]
:param vs: S-wave velocity [m/s] :returns: P-wave velocity [m/s] '''
return vs*1.16 + 1360.0
'''Calculate vp from vs using castagna's relation.
Castagna's relation (the mudrock line) is an empirical relation for vp/vs for siliciclastic rocks (i.e. sandstones and shales). [Castagna et al., 1985]
vp = 1.16 * vs + 1360 [m/s]
:param vp: P-wave velocity [m/s] :returns: S-wave velocity [m/s] '''
return (vp - 1360.0) / 1.16
if x.size < minsize: return x ry = (y.max()-y.min()) if ry == 0: return x dx = (x[1:] - x[:-1])/(x.max()-x.min()) dy = (y[1:] + y[:-1])/ry
s = num.zeros(x.size) s[1:] = num.cumsum(num.sqrt(dy**2 + dx**2)) s2 = num.linspace(0, s[-1], x.size) x2 = num.interp(s2, s, x) x2[0] = x[0] x2[-1] = x[-1] return x2
''' Create NumPy array filled with given value.
This works like :py:func:`numpy.ones` but initializes the array with ``v`` instead of ones. '''
return b return a
d = [] for x, y in zip(xx, yy): if num.isfinite(y): d.append((x, y)) return d
return xytups( x, num.interp(x, xp, fp, left=num.nan, right=num.nan)) return xytups( x, num.interp(x, xp[::-1], fp[::-1], left=num.nan, right=num.nan)) else: num.logical_and(xp[:-1] >= xv, xv > xp[1:]), num.logical_and(xp[:-1] <= xv, xv < xp[1:])))[0]
|