# -*- coding: utf-8 -*-
"""
Slowness model class.
"""
from copy import deepcopy
import math
import numpy as np
from . import _DEFAULT_VALUES
from .helper_classes import (CriticalDepth, DepthRange, SlownessLayer,
SlownessModelError, SplitLayerInfo, TimeDist)
from .slowness_layer import (bullen_depth_for,
bullen_radial_slowness, create_from_vlayer,
evaluate_at_bullen)
from .velocity_layer import (VelocityLayer, evaluate_velocity_at_bottom,
evaluate_velocity_at_top)
from .utils import _flat
[docs]
def _fix_critical_depths(critical_depths, layer_num, is_p_wave):
name = 'p_layer_num' if is_p_wave else 's_layer_num'
mask = critical_depths[name] > layer_num
critical_depths[name][mask] += 1
[docs]
class SlownessModel(object):
"""
Storage and methods for generating slowness-depth pairs.
"""
[docs]
def __init__(self, v_mod, min_delta_p=0.1, max_delta_p=11,
max_depth_interval=115,
max_range_interval=2.5 * math.pi / 180,
max_interp_error=0.05, allow_inner_core_s=True,
slowness_tolerance=_DEFAULT_VALUES["slowness_tolerance"],
skip_model_creation=False):
self.debug = False
# NB if the following are actually cleared (lists are mutable) every
# time create_sample is called, maybe it would be better to just put
# these initialisations into the relevant methods? They do have to be
# persistent across method calls in create_sample though, so don't.
# Stores the layer number for layers in the velocity model with a
# critical point at their top. These form the "branches" of slowness
# sampling.
self.critical_depths = None # will be list of CriticalDepth objects
# Store depth ranges that contains a high slowness zone for P/S. Stored
# as DepthRange objects, containing the top depth and bottom depth.
self.high_slowness_layer_depths_p = [] # will be list of DepthRanges
self.high_slowness_layer_depths_s = []
# Stores depth ranges that are fluid, ie S velocity is zero. Stored as
# DepthRange objects, containing the top depth and bottom depth.
self.fluid_layer_depths = []
self.p_layers = None
self.s_layers = None
# For methods that have an is_p_wave parameter
self.s_wave = False
self.p_wave = True
self.v_mod = v_mod
self.min_delta_p = min_delta_p
self.max_delta_p = max_delta_p
self.max_depth_interval = max_depth_interval
self.max_range_interval = max_range_interval
self.max_interp_error = max_interp_error
self.allow_inner_core_s = allow_inner_core_s
self.slowness_tolerance = slowness_tolerance
if skip_model_creation:
return
self.create_sample()
[docs]
def __str__(self):
desc = "".join([
"radius_of_planet=", str(self.radius_of_planet), "\n max_delta_p=",
str(self.max_delta_p),
"\n min_delta_p=", str(self.min_delta_p), "\n max_depth_interval=",
str(self.max_depth_interval), "\n max_range_interval=",
str(self.max_range_interval),
"\n allow_inner_core_s=", str(self.allow_inner_core_s),
"\n slownessTolerance=", str(self.slowness_tolerance),
"\n get_num_layers('P')=", str(self.get_num_layers(self.p_wave)),
"\n get_num_layers('S')=", str(self.get_num_layers(self.s_wave)),
"\n fluid_layer_depths.size()=", str(len(self.fluid_layer_depths)),
"\n high_slowness_layer_depths_p.size()=",
str(len(self.high_slowness_layer_depths_p)),
"\n high_slowness_layer_depths_s.size()=",
str(len(self.high_slowness_layer_depths_s)),
"\n critical_depths.size()=",
(str(len(self.critical_depths))
if self.critical_depths is not None else 'N/A'),
"\n"])
desc += "**** Critical Depth Layers ************************\n"
desc += str(self.critical_depths)
desc += "\n"
desc += "\n**** Fluid Layer Depths ************************\n"
for fl in self.fluid_layer_depths:
desc += str(fl.top_depth) + "," + str(fl.bot_depth) + " "
desc += "\n"
desc += "\n**** P High Slowness Layer Depths ****************\n"
for fl in self.high_slowness_layer_depths_p:
desc += str(fl.top_depth) + "," + str(fl.bot_depth) + " "
desc += "\n"
desc += "\n**** S High Slowness Layer Depths ****************\n"
for fl in self.high_slowness_layer_depths_s:
desc += str(fl.top_depth) + "," + str(fl.bot_depth) + " "
desc += "\n"
desc += "\n**** P Layers ****************\n"
for l in self.p_layers: # NOQA
desc += str(l) + "\n"
return desc
[docs]
def create_sample(self):
"""
Create slowness-depth layers from a velocity model.
This method takes a velocity model and creates a vector containing
slowness-depth layers that, hopefully, adequately sample both slowness
and depth so that the travel time as a function of distance can be
reconstructed from the theta function.
"""
# Some checks on the velocity model
self.v_mod.validate()
if len(self.v_mod) == 0:
raise SlownessModelError("velModel.get_num_layers()==0")
if self.v_mod.layers[0]['top_s_velocity'] == 0:
raise SlownessModelError(
"Unable to handle zero S velocity layers at surface. "
"This should be fixed at some point, but is a limitation of "
"TauP at this point.")
if self.debug:
print("start create_sample")
self.radius_of_planet = self.v_mod.radius_of_planet
if self.debug:
print("find_critical_points")
self.find_critical_points()
if self.debug:
print("coarse_sample")
self.coarse_sample()
if self.debug:
self.validate()
if self.debug:
print("ray_paramCheck")
self.ray_param_inc_check()
if self.debug:
print("depth_inc_check")
self.depth_inc_check()
if self.debug:
print("distance_check")
self.distance_check()
if self.debug:
print("fix_critical_points")
self.fix_critical_points()
self.validate()
if self.debug:
print("create_sample seems to be done successfully.")
[docs]
def find_critical_points(self):
"""
Find all critical points within a velocity model.
Critical points are first order discontinuities in velocity/slowness,
local extrema in slowness. A high slowness zone is a low velocity zone,
but it is possible to have a slightly low velocity zone within a
spherical planet that is not a high slowness zone and thus does not
exhibit any of the pathological behavior of a low velocity zone.
"""
high_slowness_zone_p = DepthRange()
high_slowness_zone_s = DepthRange()
fluid_zone = DepthRange()
in_fluid_zone = False
below_outer_core = False
in_high_slowness_zone_p = False
in_high_slowness_zone_s = False
# just some very big values (java had max possible of type,
# but these should do)
min_p_so_far = 1.1e300
min_s_so_far = 1.1e300
# First remove any critical points previously stored
# so these are effectively re-initialised... it's probably silly
self.critical_depths = np.zeros(len(self.v_mod.layers) + 1,
dtype=CriticalDepth)
cd_count = 0
self.high_slowness_layer_depths_p = [] # lists of DepthRange
self.high_slowness_layer_depths_s = []
self.fluid_layer_depths = []
# Initialize the current velocity layer
# to be zero thickness layer with values at the surface
curr_v_layer = self.v_mod.layers[0]
curr_v_layer = np.array([(
curr_v_layer['top_depth'], curr_v_layer['top_depth'],
curr_v_layer['top_p_velocity'], curr_v_layer['top_p_velocity'],
curr_v_layer['top_s_velocity'], curr_v_layer['top_s_velocity'],
curr_v_layer['top_density'], curr_v_layer['top_density'],
curr_v_layer['top_qp'], curr_v_layer['top_qp'],
curr_v_layer['top_qs'], curr_v_layer['top_qs'])],
dtype=VelocityLayer)
curr_s_layer = create_from_vlayer(
v_layer=curr_v_layer,
is_p_wave=self.s_wave,
radius_of_planet=self.v_mod.radius_of_planet,
is_spherical=self.v_mod.is_spherical)
curr_p_layer = create_from_vlayer(
v_layer=curr_v_layer,
is_p_wave=self.p_wave,
radius_of_planet=self.v_mod.radius_of_planet,
is_spherical=self.v_mod.is_spherical)
# We know that the top is always a critical slowness so add 0
self.critical_depths[cd_count] = (0, 0, 0, 0)
cd_count += 1
# Check to see if starting in fluid zone.
if in_fluid_zone is False and curr_v_layer['top_s_velocity'] == 0:
in_fluid_zone = True
fluid_zone = DepthRange(top_depth=curr_v_layer['top_depth'])
curr_s_layer = curr_p_layer
if min_s_so_far > curr_s_layer['top_p']:
min_s_so_far = curr_s_layer['top_p']
# P is not a typo, it represents slowness, not P-wave speed.
if min_p_so_far > curr_p_layer['top_p']:
min_p_so_far = curr_p_layer['top_p']
for layer_num, layer in enumerate(self.v_mod.layers):
prev_v_layer = curr_v_layer
prev_s_layer = curr_s_layer
prev_p_layer = curr_p_layer
# Could make the following a deep copy, but not necessary.
# Also don't just replace layer here and in the loop
# control with curr_v_layer, or the reference to the first
# zero thickness layer that has been initialised above
# will break.
curr_v_layer = layer
# Check again if in fluid zone
if in_fluid_zone is False and curr_v_layer['top_s_velocity'] == 0:
in_fluid_zone = True
fluid_zone = DepthRange(top_depth=curr_v_layer['top_depth'])
# If already in fluid zone, check if exited (java line 909)
if in_fluid_zone is True and curr_v_layer['top_s_velocity'] != 0:
if prev_v_layer['bot_depth'] > self.v_mod.iocb_depth:
below_outer_core = True
in_fluid_zone = False
fluid_zone.bot_depth = prev_v_layer['bot_depth']
self.fluid_layer_depths.append(fluid_zone)
curr_p_layer = create_from_vlayer(
v_layer=curr_v_layer,
is_p_wave=self.p_wave,
radius_of_planet=self.v_mod.radius_of_planet,
is_spherical=self.v_mod.is_spherical)
# If we are in a fluid zone ( S velocity = 0.0 ) or if we are below
# the outer core and allow_inner_core_s=false then use the P
# velocity structure to look for critical points.
if in_fluid_zone \
or (below_outer_core and self.allow_inner_core_s is False):
curr_s_layer = curr_p_layer
else:
curr_s_layer = create_from_vlayer(
v_layer=curr_v_layer,
is_p_wave=self.s_wave,
radius_of_planet=self.v_mod.radius_of_planet,
is_spherical=self.v_mod.is_spherical)
if prev_s_layer['bot_p'] != curr_s_layer['top_p'] \
or prev_p_layer['bot_p'] != curr_p_layer['top_p']:
# a first order discontinuity
self.critical_depths[cd_count] = (
curr_s_layer['top_depth'],
layer_num,
-1,
-1)
cd_count += 1
if self.debug:
print('First order discontinuity, depth =' +
str(curr_s_layer['top_depth']))
print('between' + str(prev_p_layer), str(curr_p_layer))
if in_high_slowness_zone_s and \
curr_s_layer['top_p'] < min_s_so_far:
if self.debug:
print("Top of current layer is the bottom"
" of a high slowness zone.")
high_slowness_zone_s.bot_depth = curr_s_layer['top_depth']
self.high_slowness_layer_depths_s.append(
high_slowness_zone_s)
in_high_slowness_zone_s = False
if in_high_slowness_zone_p and \
curr_p_layer['top_p'] < min_p_so_far:
if self.debug:
print("Top of current layer is the bottom"
" of a high slowness zone.")
high_slowness_zone_p.bot_depth = curr_s_layer['top_depth']
self.high_slowness_layer_depths_p.append(
high_slowness_zone_p)
in_high_slowness_zone_p = False
# Update min_p_so_far and min_s_so_far as all total reflections
# off of the top of the discontinuity are ok even though below
# the discontinuity could be the start of a high slowness zone.
if min_p_so_far > curr_p_layer['top_p']:
min_p_so_far = curr_p_layer['top_p']
if min_s_so_far > curr_s_layer['top_p']:
min_s_so_far = curr_s_layer['top_p']
if in_high_slowness_zone_s is False and (
prev_s_layer['bot_p'] < curr_s_layer['top_p'] or
curr_s_layer['top_p'] < curr_s_layer['bot_p']):
# start of a high slowness zone S
if self.debug:
print("Found S high slowness at first order " +
"discontinuity, layer = " + str(layer_num))
in_high_slowness_zone_s = True
high_slowness_zone_s = \
DepthRange(top_depth=curr_s_layer['top_depth'])
high_slowness_zone_s.ray_param = min_s_so_far
if in_high_slowness_zone_p is False and (
prev_p_layer['bot_p'] < curr_p_layer['top_p'] or
curr_p_layer['top_p'] < curr_p_layer['bot_p']):
# start of a high slowness zone P
if self.debug:
print("Found P high slowness at first order " +
"discontinuity, layer = " + str(layer_num))
in_high_slowness_zone_p = True
high_slowness_zone_p = \
DepthRange(top_depth=curr_p_layer['top_depth'])
high_slowness_zone_p.ray_param = min_p_so_far
elif ((prev_s_layer['top_p'] - prev_s_layer['bot_p']) *
(prev_s_layer['bot_p'] - curr_s_layer['bot_p']) < 0) or (
(prev_p_layer['top_p'] - prev_p_layer['bot_p']) *
(prev_p_layer['bot_p'] - curr_p_layer['bot_p'])) < 0:
# local slowness extrema, java l 1005
self.critical_depths[cd_count] = (
curr_s_layer['top_depth'],
layer_num,
-1,
-1)
cd_count += 1
if self.debug:
print("local slowness extrema, depth=" +
str(curr_s_layer['top_depth']))
# here is line 1014 of the java src!
if in_high_slowness_zone_p is False \
and curr_p_layer['top_p'] < curr_p_layer['bot_p']:
if self.debug:
print("start of a P high slowness zone, local "
"slowness extrema,min_p_so_far= " +
str(min_p_so_far))
in_high_slowness_zone_p = True
high_slowness_zone_p = \
DepthRange(top_depth=curr_p_layer['top_depth'])
high_slowness_zone_p.ray_param = min_p_so_far
if in_high_slowness_zone_s is False \
and curr_s_layer['top_p'] < curr_s_layer['bot_p']:
if self.debug:
print("start of a S high slowness zone, local "
"slowness extrema, min_s_so_far= " +
str(min_s_so_far))
in_high_slowness_zone_s = True
high_slowness_zone_s = \
DepthRange(top_depth=curr_s_layer['top_depth'])
high_slowness_zone_s.ray_param = min_s_so_far
if in_high_slowness_zone_p and \
curr_p_layer['bot_p'] < min_p_so_far:
# P: layer contains the bottom of a high slowness zone. java
# l 1043
if self.debug:
print("layer contains the bottom of a P " +
"high slowness zone. min_p_so_far=" +
str(min_p_so_far), curr_p_layer)
high_slowness_zone_p.bot_depth = self.find_depth_from_layers(
min_p_so_far, layer_num, layer_num, self.p_wave)
self.high_slowness_layer_depths_p.append(high_slowness_zone_p)
in_high_slowness_zone_p = False
if in_high_slowness_zone_s and \
curr_s_layer['bot_p'] < min_s_so_far:
# S: layer contains the bottom of a high slowness zone. java
# l 1043
if self.debug:
print("layer contains the bottom of a S " +
"high slowness zone. min_s_so_far=" +
str(min_s_so_far), curr_s_layer)
# in fluid layers we want to check p_wave structure
# when looking for S wave critical points
por_s = (self.p_wave
if curr_s_layer == curr_p_layer else self.s_wave)
high_slowness_zone_s.bot_depth = self.find_depth_from_layers(
min_s_so_far, layer_num, layer_num, por_s)
self.high_slowness_layer_depths_s.append(high_slowness_zone_s)
in_high_slowness_zone_s = False
if min_p_so_far > curr_p_layer['bot_p']:
min_p_so_far = curr_p_layer['bot_p']
if min_p_so_far > curr_p_layer['top_p']:
min_p_so_far = curr_p_layer['top_p']
if min_s_so_far > curr_s_layer['bot_p']:
min_s_so_far = curr_s_layer['bot_p']
if min_s_so_far > curr_s_layer['top_p']:
min_s_so_far = curr_s_layer['top_p']
if self.debug and in_high_slowness_zone_s:
print("In S high slowness zone, layer_num = " +
str(layer_num) + " min_s_so_far=" + str(min_s_so_far))
if self.debug and in_high_slowness_zone_p:
print("In P high slowness zone, layer_num = " +
str(layer_num) + " min_p_so_far=" + str(min_p_so_far))
# We know that the bottommost depth is always a critical slowness,
# so we add v_mod.get_num_layers()
# java line 1094
self.critical_depths[cd_count] = (
self.radius_of_planet, len(self.v_mod), -1, -1)
cd_count += 1
# Check if the bottommost depth is contained within a high slowness
# zone, might happen in a flat non-whole-planet model
if in_high_slowness_zone_s:
high_slowness_zone_s.bot_depth = curr_v_layer['bot_depth']
self.high_slowness_layer_depths_s.append(high_slowness_zone_s)
if in_high_slowness_zone_p:
high_slowness_zone_p.bot_depth = curr_v_layer['bot_depth']
self.high_slowness_layer_depths_p.append(high_slowness_zone_p)
# Check if the bottommost depth is contained within a fluid zone, this
# would be the case if we have a non whole planet model with the bottom
# in the outer core or if allow_inner_core_s == false and we want to
# use the P velocity structure in the inner core.
if in_fluid_zone:
fluid_zone.bot_depth = curr_v_layer['bot_depth']
self.fluid_layer_depths.append(fluid_zone)
self.critical_depths = self.critical_depths[:cd_count]
self.validate()
[docs]
def get_num_layers(self, is_p_wave):
"""
Number of slowness layers.
This is meant to return the number of P or S layers.
:param is_p_wave: Return P layer count (``True``) or S layer count
(``False``).
:type is_p_wave: bool
:returns: Number of slowness layers.
:rtype: int
"""
if is_p_wave:
return len(self.p_layers)
else:
return len(self.s_layers)
[docs]
def find_depth_from_depths(self, ray_param, top_depth, bot_depth,
is_p_wave):
"""
Find depth corresponding to a slowness between two given depths.
The given depths are converted to layer numbers before calling
:meth:`find_depth_from_layers`.
:param ray_param: Slowness (aka ray parameter) to find, in s/km.
:type ray_param: float
:param top_depth: Top depth to search, in km.
:type top_depth: float
:param bot_depth: Bottom depth to search, in km.
:type bot_depth: float
:param is_p_wave: ``True`` if P wave or ``False`` for S wave.
:type is_p_wave: bool
:returns: Depth (in km) corresponding to the desired slowness.
:rtype: float
:raises SlownessModelError:
If ``top_critical_layer > bot_critical_layer``
because there are no layers to search, or if there is an increase
in slowness, i.e., a negative velocity gradient, that just balances
the decrease in slowness due to the spherical planet, or if the ray
parameter ``p`` is not contained within the specified layer range.
"""
top_layer_num = self.v_mod.layer_number_below(top_depth)[0]
if self.v_mod.layers[top_layer_num]['bot_depth'] == top_depth:
top_layer_num += 1
bot_layer_num = self.v_mod.layer_number_above(bot_depth)[0]
return self.find_depth_from_layers(ray_param, top_layer_num,
bot_layer_num, is_p_wave)
[docs]
def find_depth_from_layers(self, p, top_critical_layer, bot_critical_layer,
is_p_wave):
"""
Find depth corresponding to a slowness p between two velocity layers.
Here, slowness is defined as ``(radius_of_planet-depth) / velocity``,
and sometimes called ray parameter. Both the top and the bottom
velocity layers are included. We also check to see if the slowness is
less than the bottom slowness of these layers but greater than the top
slowness of the next deeper layer. This corresponds to a total
reflection. In this case a check needs to be made to see if this is an
S wave reflecting off of a fluid layer, use P velocity below in this
case. We assume that slowness is monotonic within these layers and
therefore there is only one depth with the given slowness. This means
we return the first depth that we find.
:param p: Slowness (aka ray parameter) to find, in s/km.
:type p: float
:param top_critical_layer: Top layer number to search.
:type top_critical_layer: int
:param bot_critical_layer: Bottom layer number to search.
:type bot_critical_layer: int
:param is_p_wave: ``True`` if P wave or ``False`` for S wave.
:type is_p_wave: bool
:returns: Depth (in km) corresponding to the desired slowness.
:rtype: float
:raises SlownessModelError: If
``top_critical_layer > bot_critical_layer``
because there are no layers to search, or if there is an increase
in slowness, i.e., a negative velocity gradient, that just balances
the decrease in slowness due to the spherical planet, or if the ray
parameter ``p`` is not contained within the specified layer range.
"""
# top_p = 1.1e300 # dummy numbers
# bot_p = 1.1e300
wave_type = 'P' if is_p_wave else 'S'
if top_critical_layer > bot_critical_layer:
raise SlownessModelError(
"findDepth: no layers to search (wrong layer num?)")
for layer_num in range(top_critical_layer, bot_critical_layer + 1):
vel_layer = self.v_mod.layers[layer_num]
top_velocity = evaluate_velocity_at_top(vel_layer, wave_type)
bot_velocity = evaluate_velocity_at_bottom(vel_layer, wave_type)
top_p = self.to_slowness(top_velocity, vel_layer['top_depth'])
bot_p = self.to_slowness(bot_velocity, vel_layer['bot_depth'])
# Check to see if we are within 'chatter level' (numerical
# error) of the top or bottom and if so then return that depth.
if abs(top_p - p) < self.slowness_tolerance:
return vel_layer['top_depth']
if abs(p - bot_p) < self.slowness_tolerance:
return vel_layer['bot_depth']
if (top_p - p) * (p - bot_p) >= 0:
# Found layer containing p.
# We interpolate assuming that velocity is linear within
# this interval. So slope is the slope for velocity versus
# depth.
slope = (bot_velocity - top_velocity) / (
vel_layer['bot_depth'] - vel_layer['top_depth'])
depth = self.interpolate(p, top_velocity,
vel_layer['top_depth'], slope)
return depth
elif layer_num == top_critical_layer \
and abs(p - top_p) < self.slowness_tolerance:
# Check to see if p is just outside the topmost layer. If so
# then return the top depth.
return vel_layer['top_depth']
# Is p a total reflection? bot_p is the slowness at the bottom
# of the last velocity layer from the previous loop, set top_p
# to be the slowness at the top of the next layer.
if layer_num < len(self.v_mod) - 1:
vel_layer = self.v_mod.layers[layer_num + 1]
top_velocity = evaluate_velocity_at_top(vel_layer, wave_type)
if (is_p_wave is False and
np.any(self.depth_in_fluid(vel_layer['top_depth']))):
# Special case for S waves above a fluid. If top next
# layer is in a fluid then we should set top_velocity to
# be the P velocity at the top of the layer.
top_velocity = evaluate_velocity_at_top(vel_layer, 'P')
top_p = self.to_slowness(top_velocity, vel_layer['top_depth'])
if bot_p >= p >= top_p:
return vel_layer['top_depth']
# noinspection PyUnboundLocalVariable
if abs(p - bot_p) < self.slowness_tolerance:
# java line 1275
# Check to see if p is just outside the bottommost layer. If so
# than return the bottom depth.
print(" p is just outside the bottommost layer. This probably "
"shouldn't be allowed to happen!")
# noinspection PyUnboundLocalVariable
return vel_layer.getBotDepth()
raise SlownessModelError(
"slowness p=" + str(p) +
" is not contained within the specified layers." +
" top_critical_layer=" + str(top_critical_layer) +
" bot_critical_layer=" + str(bot_critical_layer))
[docs]
def to_slowness(self, velocity, depth):
"""
Convert velocity at some depth to slowness.
:param velocity: The velocity to convert, in km/s.
:type velocity: float
:param depth: The depth (in km) at which to perform the calculation.
Must be less than the radius of the planet defined in this
model, or the result is undefined.
:type depth: float
:returns: The slowness, in s/km.
:rtype: float
"""
if np.any(velocity == 0):
raise SlownessModelError(
"to_slowness: velocity can't be zero, at depth" +
str(depth),
"Maybe related to using S velocities in outer core?")
return (self.radius_of_planet - depth) / velocity
[docs]
def interpolate(self, p, top_velocity, top_depth, slope):
"""
Interpolate slowness to depth within a layer.
We interpolate assuming that velocity is linear within
this interval.
All parameters must be of the same shape.
:param p: The slowness to interpolate, in s/km.
:type p: :class:`float` or :class:`~numpy.ndarray`
:param top_velocity: The velocity (in km/s) at the top of the layer.
:type top_velocity: :class:`float` or :class:`~numpy.ndarray`
:param top_depth: The depth (in km) for the top of the layer.
:type top_depth: :class:`float` or :class:`~numpy.ndarray`
:param slope: The slope (in (km/s)/km) for velocity versus depth.
:type slope: :class:`float` or :class:`~numpy.ndarray`
:returns: The depth (in km) of the slowness below the layer boundary.
:rtype: :class:`float` or :class:`~numpy.ndarray`
"""
denominator = p * slope + 1
if np.any(denominator == 0):
raise SlownessModelError(
"Negative velocity gradient that just balances the slowness "
"gradient of the spherical slowness, i.e. planet flattening. "
"Instructions unclear; explode.")
else:
depth = (self.radius_of_planet +
p * (top_depth * slope - top_velocity)) / denominator
return depth
[docs]
def depth_in_fluid(self, depth):
"""
Determine if the given depth is contained within a fluid zone.
The fluid zone includes its upper boundary but not its lower boundary.
The top and bottom of the fluid zone are not returned as a DepthRange,
just like in the Java code, despite its claims to the contrary.
:param depth: The depth to check, in km.
:type depth: :class:`~numpy.ndarray`, dtype = :class:`float`
:returns: ``True`` if the depth is within a fluid zone, ``False``
otherwise.
:rtype: :class:`~numpy.ndarray` (dtype = :class:`bool`)
"""
ret = np.zeros(shape=depth.shape, dtype=np.bool_)
for elem in self.fluid_layer_depths:
ret |= (elem.top_depth <= depth) & (depth < elem.bot_depth)
return ret
[docs]
def coarse_sample(self):
"""
Create a coarse slowness sampling of the velocity model (v_mod).
The resultant slowness layers will satisfy the maximum depth increments
as well as sampling each point specified within the VelocityModel. The
P and S sampling will also be compatible.
"""
self.p_layers = create_from_vlayer(
v_layer=self.v_mod.layers,
is_p_wave=self.p_wave,
radius_of_planet=self.v_mod.radius_of_planet,
is_spherical=self.v_mod.is_spherical)
with np.errstate(divide='ignore'):
self.s_layers = create_from_vlayer(
v_layer=self.v_mod.layers,
is_p_wave=self.s_wave,
radius_of_planet=self.v_mod.radius_of_planet,
is_spherical=self.v_mod.is_spherical)
mask = self.depth_in_fluid(self.v_mod.layers['top_depth'])
if not self.allow_inner_core_s:
mask |= self.v_mod.layers['top_depth'] >= self.v_mod.iocb_depth
self.s_layers[mask] = self.p_layers[mask]
# Check for first order discontinuity. However, we only consider
# S discontinuities in the inner core if allow_inner_core_s is true.
above = self.v_mod.layers[:-1]
below = self.v_mod.layers[1:]
mask = np.logical_or(
above['bot_p_velocity'] != below['top_p_velocity'],
np.logical_and(
above['bot_s_velocity'] != below['top_s_velocity'],
np.logical_or(
self.allow_inner_core_s,
below['top_depth'] < self.v_mod.iocb_depth)))
index = np.where(mask)[0] + 1
above = above[mask]
below = below[mask]
# If we are going from a fluid to a solid or solid to fluid, e.g., core
# mantle or outer core to inner core then we need to use the P velocity
# for determining the S discontinuity.
top_s_vel = np.where(above['bot_s_velocity'] == 0,
above['bot_p_velocity'],
above['bot_s_velocity'])
bot_s_vel = np.where(below['top_s_velocity'] == 0,
below['top_p_velocity'],
below['top_s_velocity'])
# Add the layer, with zero thickness but nonzero slowness step,
# corresponding to the discontinuity.
curr_v_layer = np.empty(shape=above.shape, dtype=VelocityLayer)
curr_v_layer['top_depth'] = above['bot_depth']
curr_v_layer['bot_depth'] = above['bot_depth']
curr_v_layer['top_p_velocity'] = above['bot_p_velocity']
curr_v_layer['bot_p_velocity'] = below['top_p_velocity']
curr_v_layer['top_s_velocity'] = top_s_vel
curr_v_layer['bot_s_velocity'] = bot_s_vel
curr_v_layer['top_density'].fill(_DEFAULT_VALUES["density"])
curr_v_layer['bot_density'].fill(_DEFAULT_VALUES["density"])
curr_v_layer['top_qp'].fill(_DEFAULT_VALUES["qp"])
curr_v_layer['bot_qp'].fill(_DEFAULT_VALUES["qp"])
curr_v_layer['top_qs'].fill(_DEFAULT_VALUES["qs"])
curr_v_layer['bot_qs'].fill(_DEFAULT_VALUES["qs"])
curr_p_layer = create_from_vlayer(
v_layer=curr_v_layer,
is_p_wave=self.p_wave,
radius_of_planet=self.v_mod.radius_of_planet,
is_spherical=self.v_mod.is_spherical)
self.p_layers = np.insert(self.p_layers, index, curr_p_layer)
curr_s_layer = create_from_vlayer(
v_layer=curr_v_layer,
is_p_wave=self.s_wave,
radius_of_planet=self.v_mod.radius_of_planet,
is_spherical=self.v_mod.is_spherical)
mask2 = (above['bot_s_velocity'] == 0) & (below['top_s_velocity'] == 0)
if not self.allow_inner_core_s:
mask2 |= curr_v_layer['top_depth'] >= self.v_mod.iocb_depth
curr_s_layer = np.where(mask2, curr_p_layer, curr_s_layer)
self.s_layers = np.insert(self.s_layers, index, curr_s_layer)
# Make sure that all high slowness layers are sampled exactly
# at their bottom
for high_zone in self.high_slowness_layer_depths_s:
s_layer_num = self.layer_number_above(high_zone.bot_depth,
self.s_wave)
high_s_layer = self.s_layers[s_layer_num]
while high_s_layer['top_depth'] == high_s_layer['bot_depth'] and (
(high_s_layer['top_p'] - high_zone.ray_param) *
(high_zone.ray_param - high_s_layer['bot_p']) < 0):
s_layer_num += 1
high_s_layer = self.s_layers[s_layer_num]
if high_zone.ray_param != high_s_layer['bot_p']:
self.add_slowness(high_zone.ray_param, self.s_wave)
for high_zone in self.high_slowness_layer_depths_p:
s_layer_num = self.layer_number_above(high_zone.bot_depth,
self.p_wave)
high_s_layer = self.p_layers[s_layer_num]
while high_s_layer['top_depth'] == high_s_layer['bot_depth'] and (
(high_s_layer['top_p'] - high_zone.ray_param) *
(high_zone.ray_param - high_s_layer['bot_p']) < 0):
s_layer_num += 1
high_s_layer = self.p_layers[s_layer_num]
if high_zone.ray_param != high_s_layer['bot_p']:
self.add_slowness(high_zone.ray_param, self.p_wave)
# Make sure P and S are consistent by adding discontinuities in one to
# the other.
# Numpy 1.6 compatibility
# _tb = self.p_layers[['top_p', 'bot_p']]
_tb = np.vstack([self.p_layers['top_p'],
self.p_layers['bot_p']]).T.ravel()
uniq = np.unique(_tb)
for p in uniq:
self.add_slowness(p, self.s_wave)
# Numpy 1.6 compatibility
# _tb = self.p_layers[['top_p', 'bot_p']]
_tb = np.vstack([self.s_layers['top_p'],
self.s_layers['bot_p']]).T.ravel()
uniq = np.unique(_tb)
for p in uniq:
self.add_slowness(p, self.p_wave)
[docs]
def layer_number_above(self, depth, is_p_wave):
"""
Find the index of the slowness layer that contains the given depth.
Note that if the depth is a layer boundary, it returns the shallower
of the two or possibly more (since total reflections are zero
thickness layers) layers.
.. seealso:: :meth:`layer_number_below`
:param depth: The depth to find, in km.
:type depth: :class:`float` or :class:`~numpy.ndarray`
:param is_p_wave: Whether to look at P (``True``) velocity or S
(``False``) velocity.
:type is_p_wave: bool
:returns: The slowness layer containing the requested depth.
:rtype: :class:`int` or :class:`~numpy.ndarray` (dtype = :class:`int`,
shape = ``depth.shape``)
:raises SlownessModelError: If no layer in the slowness model contains
the given depth.
"""
if is_p_wave:
layers = self.p_layers
else:
layers = self.s_layers
# Check to make sure depth is within the range available
if np.any(depth < layers[0]['top_depth']) or \
np.any(depth > layers[-1]['bot_depth']):
raise SlownessModelError("No layer contains this depth")
found_layer_num = np.searchsorted(layers['top_depth'], depth)
mask = found_layer_num != 0
if np.isscalar(found_layer_num):
if mask:
found_layer_num -= 1
else:
found_layer_num[mask] -= 1
return found_layer_num
[docs]
def layer_number_below(self, depth, is_p_wave):
"""
Find the index of the slowness layer that contains the given depth.
Note that if the depth is a layer boundary, it returns the deeper of
the two or possibly more (since total reflections are zero thickness
layers) layers.
.. seealso:: :meth:`layer_number_above`
:param depth: The depth to find, in km.
:type depth: :class:`float` or :class:`~numpy.ndarray`
:param is_p_wave: Whether to look at P (``True``) velocity or S
(``False``) velocity.
:type is_p_wave: bool
:returns: The slowness layer containing the requested depth.
:rtype: :class:`int` or :class:`~numpy.ndarray` (dtype = :class:`int`,
shape = ``depth.shape``)
:raises SlownessModelError: If no layer in the slowness model contains
the given depth.
"""
if is_p_wave:
layers = self.p_layers
else:
layers = self.s_layers
# Check to make sure depth is within the range available
if np.any(depth < layers[0]['top_depth']) or \
np.any(depth > layers[-1]['bot_depth']):
raise SlownessModelError("No layer contains this depth")
found_layer_num = np.searchsorted(layers['bot_depth'], depth,
side='right')
mask = found_layer_num == layers.shape[0]
if np.isscalar(found_layer_num):
if mask:
found_layer_num -= 1
else:
found_layer_num[mask] -= 1
return found_layer_num
[docs]
def get_slowness_layer(self, layer, is_p_wave):
"""
Return the Slowness_layer of the requested wave type.
This is not meant to be a clone!
:param layer: The number of the layer(s) to return.
:type layer: :class:`int` or :class:`~numpy.ndarray` (dtype =
:class:`int`)
:param is_p_wave: Whether to return the P layer (``True``) or the S
layer (``False``).
:type is_p_wave: bool
:returns: The slowness layer(s).
:rtype: :class:`~numpy.ndarray`
(dtype = :class:`obspy.taup.helper_classes.SlownessLayer`,
shape = ``layer_num.shape``)
"""
if is_p_wave:
return self.p_layers[layer]
else:
return self.s_layers[layer]
[docs]
def add_slowness(self, p, is_p_wave):
"""
Add a ray parameter to the slowness sampling for the given wave type.
Slowness layers are split as needed and P and S sampling are kept
consistent within fluid layers. Note, this makes use of the velocity
model, so all interpolation is linear in velocity, not in slowness!
:param p: The slowness value to add, in s/km.
:type p: float
:param is_p_wave: Whether to add to the P wave (``True``) or the S wave
(``False``) sampling.
:type is_p_wave: bool
"""
if is_p_wave:
# NB Unlike Java (unfortunately) these are not modified in place!
# NumPy arrays cannot have values inserted in place.
layers = self.p_layers
other_layers = self.s_layers
wave = 'P'
else:
layers = self.s_layers
other_layers = self.p_layers
wave = 'S'
# If depths are the same only need top_velocity, and just to verify we
# are not in a fluid.
nonzero = layers['top_depth'] != layers['bot_depth']
above = self.v_mod.evaluate_above(layers['bot_depth'], wave)
below = self.v_mod.evaluate_below(layers['top_depth'], wave)
top_velocity = np.where(nonzero, below, above)
bot_velocity = np.where(nonzero, above, below)
mask = ((layers['top_p'] - p) * (p - layers['bot_p'])) > 0
# Don't need to check for S waves in a fluid or in inner core if
# allow_inner_core_s is False.
if not is_p_wave:
mask &= top_velocity != 0
if not self.allow_inner_core_s:
iocb_mask = layers['bot_depth'] > self.v_mod.iocb_depth
mask &= ~iocb_mask
index = np.where(mask)[0]
bot_depth = np.copy(layers['bot_depth'])
# Not a zero thickness layer, so calculate the depth for
# the ray parameter.
slope = ((bot_velocity[nonzero] - top_velocity[nonzero]) /
(layers['bot_depth'][nonzero] - layers['top_depth'][nonzero]))
bot_depth[nonzero] = self.interpolate(p, top_velocity[nonzero],
layers['top_depth'][nonzero],
slope)
bot_layer = np.empty(shape=index.shape, dtype=SlownessLayer)
bot_layer['top_p'].fill(p)
bot_layer['top_depth'] = bot_depth[mask]
bot_layer['bot_p'] = layers['bot_p'][mask]
bot_layer['bot_depth'] = layers['bot_depth'][mask]
top_layer = np.empty(shape=index.shape, dtype=SlownessLayer)
top_layer['top_p'] = layers['top_p'][mask]
top_layer['top_depth'] = layers['top_depth'][mask]
top_layer['bot_p'].fill(p)
top_layer['bot_depth'] = bot_depth[mask]
# numpy 1.6 compatibility
other_index = np.where(other_layers.reshape(1, -1) ==
layers[mask].reshape(-1, 1))
layers[index] = bot_layer
layers = np.insert(layers, index, top_layer)
if len(other_index[0]):
other_layers[other_index[1]] = bot_layer[other_index[0]]
other_layers = np.insert(other_layers, other_index[1],
top_layer[other_index[0]])
if is_p_wave:
self.p_layers = layers
self.s_layers = other_layers
else:
self.s_layers = layers
self.p_layers = other_layers
[docs]
def ray_param_inc_check(self):
"""
Check that no slowness layer's ray parameter interval is too large.
The limit is determined by ``self.max_delta_p``.
"""
for wave in [self.s_wave, self.p_wave]:
# These might change with calls to add_slowness, so be sure we have
# the correct copy.
if wave == self.p_wave:
layers = self.p_layers
else:
layers = self.s_layers
diff = layers['top_p'] - layers['bot_p']
absdiff = np.abs(diff)
mask = absdiff > self.max_delta_p
diff = diff[mask]
absdiff = absdiff[mask]
top_p = layers['top_p'][mask]
new_count = np.ceil(absdiff / self.max_delta_p).astype(np.int_)
steps = diff / new_count
for start, n, delta in zip(top_p, new_count, steps):
for j in range(1, n):
newp = start + j * delta
self.add_slowness(newp, self.p_wave)
self.add_slowness(newp, self.s_wave)
[docs]
def depth_inc_check(self):
"""
Check that no slowness layer is too thick.
The maximum is determined by ``self.max_depth_interval``.
"""
for wave in [self.s_wave, self.p_wave]:
# These might change with calls to add_slowness, so be sure we
# have the correct copy.
if wave == self.p_wave:
layers = self.p_layers
else:
layers = self.s_layers
diff = layers['bot_depth'] - layers['top_depth']
mask = diff > self.max_depth_interval
diff = diff[mask]
top_depth = layers['top_depth'][mask]
new_count = np.ceil(diff / self.max_depth_interval).astype(np.int_)
steps = diff / new_count
for start, nd, delta in zip(top_depth, new_count, steps):
new_depth = start + np.arange(1, nd) * delta
if wave == self.s_wave:
velocity = self.v_mod.evaluate_above(new_depth, 'S')
smask = velocity == 0
if not self.allow_inner_core_s:
smask |= new_depth >= self.v_mod.iocb_depth
if np.any(smask):
velocity[smask] = self.v_mod.evaluate_above(
new_depth[smask], 'P')
slowness = self.to_slowness(velocity, new_depth)
else:
slowness = self.to_slowness(
self.v_mod.evaluate_above(new_depth, 'P'),
new_depth)
for p in slowness:
self.add_slowness(p, self.p_wave)
self.add_slowness(p, self.s_wave)
[docs]
def distance_check(self):
"""
Check that no slowness layer is too wide or undersampled.
The width must be less than ``self.max_range_interval`` and the
(estimated) error due to linear interpolation must be less than
``self.max_interp_error``.
"""
for curr_wave_type in [self.s_wave, self.p_wave]:
is_curr_ok = False
is_prev_ok = False
prev_prev_id = None
prev_td = None
curr_td = None
j = 0
s_layer = self.get_slowness_layer(j, curr_wave_type)
while j < self.get_num_layers(curr_wave_type):
prev_s_layer = s_layer
s_layer = self.get_slowness_layer(j, curr_wave_type)
if (self.depth_in_high_slowness(s_layer['bot_depth'],
s_layer['bot_p'],
curr_wave_type) is False and
self.depth_in_high_slowness(s_layer['top_depth'],
s_layer['top_p'],
curr_wave_type) is False):
# Don't calculate prevTD if we can avoid it
if is_curr_ok and curr_td is not None:
if is_prev_ok:
prev_prev_id = prev_td
else:
prev_prev_id = None
prev_td = curr_td
is_prev_ok = True
else:
prev_td = self.approx_distance(j - 1, s_layer['top_p'],
curr_wave_type)
is_prev_ok = True
curr_td = self.approx_distance(j, s_layer['bot_p'],
curr_wave_type)
is_curr_ok = True
# Check for jump of too great distance
if (abs(prev_td['dist'] - curr_td['dist']) >
self.max_range_interval and
abs(s_layer['top_p'] - s_layer['bot_p']) >
2 * self.min_delta_p):
if self.debug:
print("At " + str(j) + " Distance jump too great ("
">max_range_interval " +
str(self.max_range_interval) + "), adding "
"slowness. ")
p = (s_layer['top_p'] + s_layer['bot_p']) / 2
self.add_slowness(p, self.p_wave)
self.add_slowness(p, self.s_wave)
curr_td = prev_td
is_curr_ok = is_prev_ok
prev_td = prev_prev_id
prev_prev_id = None
is_prev_ok = prev_td is not None
else:
# Make guess as to error estimate due to linear
# interpolation if it is not ok, then we split both
# the previous and current slowness layers, this has
# the nice, if unintended, consequence of adding
# extra samples in the neighborhood of poorly
# sampled caustics.
split_ray_param = \
(s_layer['top_p'] + s_layer['bot_p']) / 2
all_but_layer = self.approx_distance(
j - 1, split_ray_param, curr_wave_type)
split_layer = np.array([(
s_layer['top_p'], s_layer['top_depth'],
split_ray_param,
bullen_depth_for(s_layer, split_ray_param,
self.radius_of_planet))],
dtype=SlownessLayer)
just_layer_time, just_layer_dist = \
bullen_radial_slowness(
split_layer, split_ray_param,
self.radius_of_planet)
split_time = \
all_but_layer['time'] + 2 * just_layer_time
split_dist = \
all_but_layer['dist'] + 2 * just_layer_dist
# Python standard division is not IEEE compliant,
# as “The IEEE 754 standard specifies that every
# floating point arithmetic operation, including
# division by zero, has a well-defined result”.
# Use numpy's division instead by using np.array:
with np.errstate(divide='ignore', invalid='ignore'):
diff = (curr_td['time'] -
((split_time - prev_td['time']) *
((curr_td['dist'] - prev_td['dist']) /
(split_dist - prev_td['dist'])) +
prev_td['time']))
if abs(diff) > self.max_interp_error:
p1 = (prev_s_layer['top_p'] +
prev_s_layer['bot_p']) / 2
p2 = (s_layer['top_p'] + s_layer['bot_p']) / 2
self.add_slowness(p1, self.p_wave)
self.add_slowness(p1, self.s_wave)
self.add_slowness(p2, self.p_wave)
self.add_slowness(p2, self.s_wave)
curr_td = prev_prev_id
is_curr_ok = curr_td is not None
is_prev_ok = False
if j > 0:
# Back up one step unless we are at the
# beginning, then stay put.
j -= 1
s_layer = self.get_slowness_layer(
j - 1 if j - 1 >= 0 else 0, curr_wave_type)
# This s_layer will become prevSLayer in the
# next loop.
else:
is_prev_ok = False
is_curr_ok = False
else:
j += 1
if self.debug and j % 10 == 0:
print(j)
else:
prev_prev_id = None
prev_td = None
curr_td = None
is_curr_ok = False
is_prev_ok = False
j += 1
if self.debug and j % 100 == 0:
print(j)
if self.debug:
print("Number of " + ("P" if curr_wave_type else "S") +
" slowness layers: " + str(j))
[docs]
def depth_in_high_slowness(self, depth, ray_param, is_p_wave,
return_depth_range=False):
"""
Determine if depth and slowness are within a high slowness zone.
Whether the high slowness zone includes its upper boundary and its
lower boundaries depends upon the ray parameter. The slowness at the
depth is needed because if depth happens to correspond to a
discontinuity that marks the bottom of the high slowness zone but the
ray is actually a total reflection then it is not part of the high
slowness zone. The ray parameter that delimits the zone, i.e., it can
turn at the top and the bottom, is in the zone at the top, but out of
the zone at the bottom. (?)
:param depth: The depth to check, in km.
:type depth: float
:param ray_param: The slowness to check, in s/km.
:type ray_param: float
:param is_p_wave: Whether to check the P wave (``True``) or the S wave
(``False``).
:type is_p_wave: bool
:param return_depth_range: Whether to also return the DepthRange of
the high slowness zone.
:returns: ``True`` if within a high slowness zone, ``False`` otherwise.
If return_depth_range is ``True``, also returns a DepthRange object
:rtype: bool, or (bool, DepthRange)
"""
if is_p_wave:
high_slowness_layer_depths = self.high_slowness_layer_depths_p
else:
high_slowness_layer_depths = self.high_slowness_layer_depths_s
for temp_range in high_slowness_layer_depths:
if temp_range.top_depth <= depth <= temp_range.bot_depth:
if ray_param > temp_range.ray_param \
or (ray_param == temp_range.ray_param and
depth == temp_range.top_depth):
if return_depth_range:
return True, temp_range
return True
if return_depth_range:
return False, None
return False
[docs]
def approx_distance(self, slowness_turn_layer, p, is_p_wave):
"""
Approximate distance for ray turning at the bottom of a layer.
Generates approximate distance, in radians, for a ray from a surface
source that turns at the bottom of the given slowness layer.
:param slowness_turn_layer: The number of the layer at which the ray
should turn.
:type slowness_turn_layer: int
:param p: The slowness to calculate, in s/km.
:type p: float
:param is_p_wave: Whether to use the P (``True``) or S (``False``)
wave.
:type is_p_wave: bool
:returns: The time (in s) and distance (in rad) the ray travels.
:rtype: :class:`~numpy.ndarray`
(dtype = :class:`obspy.taup.helper_classes.TimeDist`, shape =
(``slowness_turn_layer``, ))
"""
# First, if the slowness model contains less than slowness_turn_layer
# elements we can't calculate a distance.
if slowness_turn_layer >= self.get_num_layers(is_p_wave):
raise SlownessModelError(
"Can't calculate a distance when get_num_layers() is smaller "
"than the given slowness_turn_layer!")
if p < 0:
raise SlownessModelError("Ray parameter must not be negative!")
td = np.zeros(1, dtype=TimeDist)
td['p'] = p
layer_num = np.arange(0, slowness_turn_layer + 1)
if len(layer_num):
time, dist = self.layer_time_dist(p, layer_num, is_p_wave)
# Return 2* distance and time because there is a downgoing as well
# as an upgoing leg, which are equal since this is for a surface
# source.
td['time'] = 2 * np.sum(time)
td['dist'] = 2 * np.sum(dist)
return td
[docs]
def layer_time_dist(self, spherical_ray_param, layer_num, is_p_wave,
check=True, allow_turn=False):
"""
Calculate time and distance for a ray passing through a layer.
Calculates the time and distance increments accumulated by a ray of
spherical ray parameter ``p`` when passing through layer ``layer_num``.
Note that this gives half of the true range and time increments since
there will be both an upgoing and a downgoing path. It also only does
the calculation for the simple cases of the centre of the planet, where
the ray parameter is zero, or for constant velocity layers. Otherwise,
it calls :func:`~.bullen_radial_slowness`.
Either ``spherical_ray_param`` or ``layer_num`` must be 0-D, or they
must have the same shape.
:param spherical_ray_param: The spherical ray parameter of the
ray(s), in s/km.
:type spherical_ray_param: :class:`float` or :class:`~numpy.ndarray`
:param layer_num: The layer(s) in which the calculation should be done.
:type layer_num: :class:`float` or :class:`~numpy.ndarray`
:param is_p_wave: Whether to look at the P (``True``) or S (``False``)
wave.
:type is_p_wave: bool
:param check: Whether to perform checks of input consistency.
:type check: bool
:param allow_turn: Whether to allow the ray to turn in the middle of a
layer.
:type allow_turn: bool
:returns: The time (in s) and distance (in rad) increments for the
specified ray(s) and layer(s).
:rtype: :class:`~numpy.ndarray`
(dtype = :class:`obspy.taup.helper_classes.TimeDist`, shape =
``spherical_ray_param.shape`` or ``layer_num.shape``)
:raises SlownessModelError: If the ray with the given spherical ray
parameter cannot propagate within this layer, or if the ray turns
within this layer but not at the bottom. These checks may be
bypassed by specifying ``check=False``.
"""
spherical_layer = self.get_slowness_layer(layer_num, is_p_wave)
pdim = np.ndim(spherical_ray_param)
ldim = np.ndim(layer_num)
if ldim == 1 and pdim == 0:
time = np.empty(shape=layer_num.shape, dtype=np.float_)
dist = np.empty(shape=layer_num.shape, dtype=np.float_)
elif ldim == 0 and pdim == 1:
time = np.empty(shape=spherical_ray_param.shape, dtype=np.float_)
dist = np.empty(shape=spherical_ray_param.shape, dtype=np.float_)
elif ldim == pdim and (ldim == 0 or
layer_num.shape == spherical_ray_param.shape):
time = np.empty(shape=layer_num.shape, dtype=np.float_)
dist = np.empty(shape=layer_num.shape, dtype=np.float_)
else:
raise TypeError('Either spherical_ray_param or layer_num must be '
'0D, or they must have the same shape.')
# First make sure that a ray with this ray param can propagate
# within this layer and doesn't turn in the middle of the layer. If
# not, raise error.
if check:
if not allow_turn:
minp = np.minimum(spherical_layer['top_p'],
spherical_layer['bot_p'])
if np.any(spherical_ray_param > minp):
raise SlownessModelError(
'Ray turns in the middle of this layer! '
'layer_num = %d' % (layer_num, ))
if np.any(spherical_ray_param > spherical_layer['top_p']):
raise SlownessModelError('Ray cannot propagate within this '
' layer, given ray param too large.')
if np.any(spherical_ray_param < 0):
raise SlownessModelError("Ray parameter must not be negative!")
turning_layers = spherical_ray_param > spherical_layer['bot_p']
if np.any(turning_layers):
if ldim == 1 and pdim == 0:
# Turn in a layer, create temp layers with p at bottom.
tmp_layers = spherical_layer[turning_layers]
turn_depth = bullen_depth_for(
tmp_layers, spherical_ray_param, self.radius_of_planet,
check=False)
spherical_layer['bot_p'][turning_layers] = spherical_ray_param
spherical_layer['bot_depth'][turning_layers] = turn_depth
elif ldim == 0 and pdim == 1:
# Turn in layer, create temp layers with each p at bottom.
# Expand layer array so that each one turns at correct depth.
ldim = 1
new_layers = np.repeat(spherical_layer,
len(spherical_ray_param))
turn_depth = bullen_depth_for(
new_layers, spherical_ray_param, self.radius_of_planet,
check=False)
new_layers['bot_p'][turning_layers] = \
spherical_ray_param[turning_layers]
new_layers['bot_depth'][turning_layers] = \
turn_depth[turning_layers]
spherical_layer = new_layers
elif ldim == pdim == 0:
# Turn in layer, create temp layer with p at bottom.
try:
turn_depth = bullen_depth_for(spherical_layer,
spherical_ray_param,
self.radius_of_planet)
except SlownessModelError:
if check:
raise
else:
turn_depth = np.nan
spherical_layer['bot_p'] = spherical_ray_param
spherical_layer['bot_depth'] = turn_depth
else:
# Turn in layer, create temp layers with each p at bottom.
turn_depth = bullen_depth_for(
spherical_layer, spherical_ray_param,
self.radius_of_planet, check=False)
turning_layers = np.where(turning_layers)
spherical_layer['bot_p'][turning_layers] = \
spherical_ray_param[turning_layers]
spherical_layer['bot_depth'][turning_layers] = \
turn_depth[turning_layers]
if check and np.any(
spherical_ray_param > np.maximum(spherical_layer['top_p'],
spherical_layer['bot_p'])):
raise SlownessModelError("Ray cannot propagate within this layer, "
"given ray param too large.")
# Check to see if this layer has zero thickness, if so then it is
# from a critically reflected slowness sample. That means just
# return 0 for time and distance increments.
zero_thick = \
spherical_layer['top_depth'] == spherical_layer['bot_depth']
if ldim == 0:
if zero_thick:
time.fill(0)
dist.fill(0)
return time, dist
else:
zero_thick = np.zeros(shape=time.shape, dtype=np.bool_)
leftover = ~zero_thick
time[zero_thick] = 0
dist[zero_thick] = 0
# Check to see if this layer contains the centre of the planet. If so
# then the spherical ray parameter should be 0.0 and we calculate the
# range and time increments using a constant velocity layer (sphere).
# See eqns. 43 and 44 of [Buland1983]_, although we implement them
# slightly differently. Note that the distance and time increments are
# for just downgoing or just upgoing, i.e. from the top of the layer
# to the centre of the planet or vice versa but not both. This is in
# keeping with the convention that these are one way distance and time
# increments. We will multiply the result by 2 at the end, or if we are
# doing a 1.5D model, the other direction may be different. The time
# increment for a ray of zero ray parameter passing half way through a
# sphere of constant velocity is just the spherical slowness at the top
# of the sphere. An amazingly simple result!
centre_layer = np.logical_and(leftover, np.logical_and(
spherical_ray_param == 0,
spherical_layer['bot_depth'] == self.radius_of_planet))
leftover &= ~centre_layer
if np.any(layer_num[centre_layer] !=
self.get_num_layers(is_p_wave) - 1):
raise SlownessModelError("There are layers deeper than the "
"centre of the planet!")
time[centre_layer] = spherical_layer['top_p'][centre_layer]
dist[centre_layer] = math.pi / 2
# Now we check to see if this is a constant velocity layer and if so
# than we can do a simple triangle calculation to get the range and
# time increments. To get the time increment we first calculate the
# path length through the layer using the law of cosines, noting
# that the angle at the top of the layer can be obtained from the
# spherical Snell's Law. The time increment is just the path length
# divided by the velocity. To get the distance we first find the
# angular distance traveled, using the law of sines.
top_radius = self.radius_of_planet - spherical_layer['top_depth']
bot_radius = self.radius_of_planet - spherical_layer['bot_depth']
with np.errstate(invalid='ignore'):
vel = bot_radius / spherical_layer['bot_p']
constant_velocity = np.logical_and(
leftover,
np.abs(top_radius / spherical_layer['top_p'] -
vel) < self.slowness_tolerance)
leftover &= ~constant_velocity
top_radius = top_radius[constant_velocity]
bot_radius = bot_radius[constant_velocity]
vel = vel[constant_velocity]
if pdim:
ray_param_const_velocity = spherical_ray_param[constant_velocity]
else:
ray_param_const_velocity = spherical_ray_param
# In cases of a ray turning at the bottom of the layer numerical
# round-off can cause bot_term to be very small (1e-9) but
# negative which causes the sqrt to raise an error. We check for
# values that are within the numerical chatter of zero and just
# set them to zero.
top_term = top_radius ** 2 - (ray_param_const_velocity * vel) ** 2
top_term[np.abs(top_term) < self.slowness_tolerance] = 0
# In this case the ray turns at the bottom of this layer so
# spherical_ray_param*vel == bot_radius, and bot_term should be
# zero. We check for this case specifically because
# numerical chatter can cause small round-off errors that
# lead to bot_term being negative, causing a sqrt error.
bot_term = np.zeros(shape=top_term.shape)
mask = (ray_param_const_velocity !=
spherical_layer['bot_p'][constant_velocity])
if pdim:
bot_term[mask] = bot_radius[mask] ** 2 - (
ray_param_const_velocity[mask] * vel[mask]) ** 2
else:
bot_term[mask] = bot_radius[mask] ** 2 - (
ray_param_const_velocity * vel[mask]) ** 2
b = np.sqrt(top_term) - np.sqrt(bot_term)
time[constant_velocity] = b / vel
dist[constant_velocity] = np.arcsin(
b * ray_param_const_velocity * vel / (top_radius * bot_radius))
# If the layer is not a constant velocity layer or the centre of the
# planet and p is not zero we have to do it the hard way:
time[leftover], dist[leftover] = bullen_radial_slowness(
spherical_layer[leftover] if ldim else spherical_layer,
spherical_ray_param[leftover] if pdim else spherical_ray_param,
self.radius_of_planet,
check=check)
if check and (np.any(time < 0) or np.any(np.isnan(time)) or
np.any(dist < 0) or np.any(np.isnan(dist))):
raise SlownessModelError(
"layer time|dist < 0 or NaN.")
return time, dist
[docs]
def fix_critical_points(self):
"""
Reset the slowness layers that correspond to critical points.
"""
self.critical_depths['p_layer_num'] = self.layer_number_below(
self.critical_depths['depth'],
self.p_wave)
s_layer = self.get_slowness_layer(self.critical_depths['p_layer_num'],
self.p_wave)
# We want the last critical point to be the bottom of the last layer.
mask = (
(self.critical_depths['p_layer_num'] == len(self.p_layers) - 1) &
(s_layer['bot_depth'] == self.critical_depths['depth']))
self.critical_depths['p_layer_num'][mask] += 1
self.critical_depths['s_layer_num'] = self.layer_number_below(
self.critical_depths['depth'],
self.s_wave)
s_layer = self.get_slowness_layer(self.critical_depths['s_layer_num'],
self.s_wave)
# We want the last critical point to be the bottom of the last layer.
mask = (
(self.critical_depths['s_layer_num'] == len(self.s_layers) - 1) &
(s_layer['bot_depth'] == self.critical_depths['depth']))
self.critical_depths['s_layer_num'][mask] += 1
[docs]
def validate(self):
"""
Perform consistency check on the slowness model.
In Java, there is a separate validate method defined in the
SphericalSModel subclass and as such overrides the validate in
SlownessModel, but it itself calls the super method (by
super.validate()), i.e. the code above. Both are merged here (in
fact, it only contained one test).
"""
if self.radius_of_planet <= 0:
raise SlownessModelError("Radius of planet must be positive.")
if self.max_depth_interval <= 0:
raise SlownessModelError(
"max_depth_interval must be positive and non-zero.")
# Check for inconsistencies in high slowness zones.
for is_p_wave in [self.p_wave, self.s_wave]:
if is_p_wave:
high_slowness_layer_depths = self.high_slowness_layer_depths_p
else:
high_slowness_layer_depths = self.high_slowness_layer_depths_s
prev_depth = -1e300
for high_s_zone_depth in high_slowness_layer_depths:
if high_s_zone_depth.top_depth >= high_s_zone_depth.bot_depth:
raise SlownessModelError(
"High Slowness zone has zero or negative thickness!")
if (high_s_zone_depth.top_depth < prev_depth or (
high_s_zone_depth.top_depth == prev_depth and not
self.v_mod.is_discontinuity(
high_s_zone_depth.top_depth))):
raise SlownessModelError(
"High Slowness zone overlaps previous zone.")
prev_depth = high_s_zone_depth.bot_depth
# Check for inconsistencies in fluid zones.
prev_depth = -1e300
for fluid_zone in self.fluid_layer_depths:
if fluid_zone.top_depth >= fluid_zone.bot_depth:
raise SlownessModelError(
"Fluid zone has zero or negative thickness!")
if fluid_zone.top_depth <= prev_depth:
raise SlownessModelError("Fluid zone overlaps previous zone.")
prev_depth = fluid_zone.bot_depth
# Check for inconsistencies in slowness layers.
for layers in [self.p_layers, self.s_layers]:
if layers is None:
continue
if np.any(np.isnan(layers['top_p']) | np.isnan(layers['bot_p'])):
raise SlownessModelError("Slowness layer has NaN values.")
if np.any((layers['top_p'] < 0) | (layers['bot_p'] < 0)):
raise SlownessModelError(
"Slowness layer has negative slowness.")
if np.any(layers['top_p'][1:] != layers['bot_p'][:-1]):
raise SlownessModelError(
"Slowness layer slowness does not agree with "
"previous layer (at same depth)!")
if np.any(layers['top_depth'] > layers['bot_depth']):
raise SlownessModelError(
"Slowness layer has negative thickness.")
if layers['top_depth'][0] > 0:
raise SlownessModelError("Gap between slowness layers!")
if np.any(layers['top_depth'][1:] > layers['bot_depth'][:-1]):
raise SlownessModelError("Gap between slowness layers!")
if layers['top_depth'][0] < 0:
raise SlownessModelError("Slowness layer overlaps previous!")
if np.any(layers['top_depth'][1:] < layers['bot_depth'][:-1]):
raise SlownessModelError("Slowness layer overlaps previous!")
if np.any(np.isnan(layers['top_depth']) |
np.isnan(layers['bot_depth'])):
raise SlownessModelError(
"Slowness layer depth (top or bottom) is NaN!")
if np.any(layers['bot_depth'] > self.radius_of_planet):
raise SlownessModelError(
"Slowness layer has a depth larger than radius of the "
"planet.")
# Everything seems OK.
return True
[docs]
def get_min_turn_ray_param(self, depth, is_p_wave):
"""
Find minimum slowness, turning but not reflected, at or above a depth.
Normally this is the slowness sample at the given depth, but if the
depth is within a high slowness zone, then it may be smaller.
:param depth: The depth to search for, in km.
:type depth: float
:param is_p_wave: Whether to search the P (``True``) or S (``False``)
wave.
:type is_p_wave: bool
:returns: The minimum ray parameter, in s/km.
:rtype: float
"""
min_p_so_far = 1e300
if self.depth_in_high_slowness(depth, 1e300, is_p_wave):
for s_layer in (self.p_layers if is_p_wave else self.s_layers):
if s_layer['bot_depth'] == depth:
min_p_so_far = min(min_p_so_far, s_layer['bot_p'])
return min_p_so_far
elif s_layer['bot_depth'] > depth:
min_p_so_far = min(
min_p_so_far,
evaluate_at_bullen(s_layer, depth,
self.radius_of_planet))
return min_p_so_far
else:
min_p_so_far = min(min_p_so_far, s_layer['bot_p'])
else:
s_layer = self.get_slowness_layer(
self.layer_number_above(depth, is_p_wave), is_p_wave)
if depth == s_layer['bot_depth']:
min_p_so_far = s_layer['bot_p']
else:
min_p_so_far = evaluate_at_bullen(s_layer, depth,
self.radius_of_planet)
return min_p_so_far
[docs]
def get_min_ray_param(self, depth, is_p_wave):
"""
Find minimum slowness, turning or reflected, at or above a depth.
Normally this is the slowness sample at the given depth, but if the
depth is within a high slowness zone, then it may be smaller. Also, at
first order discontinuities, there may be many slowness samples at the
same depth.
:param depth: The depth to search for, in km.
:type depth: float
:param is_p_wave: Whether to search the P (``True``) or S (``False``)
wave.
:type is_p_wave: bool
:returns: The minimum ray parameter, in s/km.
:rtype: float
"""
min_p_so_far = self.get_min_turn_ray_param(depth, is_p_wave)
s_layer_above = self.get_slowness_layer(
self.layer_number_above(depth, is_p_wave), is_p_wave)
s_layer_below = self.get_slowness_layer(
self.layer_number_below(depth, is_p_wave), is_p_wave)
if s_layer_above['bot_depth'] == depth:
min_p_so_far = min(min_p_so_far, s_layer_above['bot_p'],
s_layer_below['top_p'])
return min_p_so_far
[docs]
def split_layer(self, depth, is_p_wave):
"""
Split a slowness layer into two slowness layers.
The interpolation for splitting a layer is a Bullen p=Ar^B and so does
not directly use information from the VelocityModel.
:param depth: The depth at which attempt a split, in km.
:type depth: float
:param is_p_wave: Whether to split based on P (``True``) or S
(``False``) wave.
:type is_p_wave: bool
:returns: Information about the split as (or if) it was performed, such
that:
* ``needed_split=True`` if a layer was actually split;
* ``moved_sample=True`` if a layer was very close, and so moving
the layer's depth is better than making a very thin layer;
* ``ray_param=...``, the new ray parameter (in s/km), if the layer
was split.
:rtype: :class:`~.SplitLayerInfo`
"""
layer_num = self.layer_number_above(depth, is_p_wave)
s_layer = self.get_slowness_layer(layer_num, is_p_wave)
if s_layer['top_depth'] == depth or s_layer['bot_depth'] == depth:
# Depth is already on a slowness layer boundary so no need to
# split any slowness layers.
return SplitLayerInfo(self, False, False, 0)
elif abs(s_layer['top_depth'] - depth) < 0.000001:
# Check for very thin layers, just move the layer to hit the
# boundary.
out = deepcopy(self)
out_layers = out.p_layers if is_p_wave else out.s_layers
out_layers[layer_num] = (s_layer['top_p'], depth,
s_layer['bot_p'], s_layer['bot_depth'])
s_layer = self.get_slowness_layer(layer_num - 1, is_p_wave)
out_layers[layer_num - 1] = (s_layer['top_p'],
s_layer['top_depth'],
s_layer['bot_p'], depth)
return SplitLayerInfo(out, False, True, s_layer['bot_p'])
elif abs(depth - s_layer['bot_depth']) < 0.000001:
# As above.
out = deepcopy(self)
out_layers = out.p_layers if is_p_wave else out.s_layers
out_layers[layer_num] = (s_layer['top_p'], s_layer['top_depth'],
s_layer['bot_p'], depth)
s_layer = self.get_slowness_layer(layer_num + 1, is_p_wave)
out_layers[layer_num + 1] = (s_layer['top_p'], depth,
s_layer['bot_p'],
s_layer['bot_depth'])
return SplitLayerInfo(out, False, True, s_layer['bot_p'])
else:
# Must split properly.
out = deepcopy(self)
p = evaluate_at_bullen(s_layer, depth, self.radius_of_planet)
top_layer = np.array([(s_layer['top_p'], s_layer['top_depth'],
p, depth)],
dtype=SlownessLayer)
bot_layer = (p, depth, s_layer['bot_p'], s_layer['bot_depth'])
out_layers = out.p_layers if is_p_wave else out.s_layers
out_layers[layer_num] = bot_layer
out_layers = np.insert(out_layers, layer_num, top_layer)
# Fix critical layers since we added a slowness layer.
out_critical_depths = self.critical_depths
_fix_critical_depths(out_critical_depths, layer_num, is_p_wave)
if is_p_wave:
out_p_layers = out_layers
out_s_layers = self._fix_other_layers(out.s_layers, p, s_layer,
top_layer, bot_layer,
out_critical_depths,
False)
else:
out_p_layers = self._fix_other_layers(out.p_layers, p, s_layer,
top_layer, bot_layer,
out_critical_depths,
True)
out_s_layers = out_layers
out.critical_depths = out_critical_depths
out.p_layers = out_p_layers
out.s_layers = out_s_layers
return SplitLayerInfo(out, True, False, p)
[docs]
def _fix_other_layers(self, other_layers, p, changed_layer, new_top_layer,
new_bot_layer, critical_depths, is_p_wave):
"""
Fix other wave layers when a split is made.
This performs the second split of the *other* wave type when a split is
made by :meth:`split_layer`.
"""
out = other_layers
# Make sure to keep sampling consistent. If in a fluid, both wave
# types will share a single slowness layer.
other_index = np.where(other_layers == changed_layer)
if len(other_index[0]):
out[other_index[0]] = new_bot_layer
out = np.insert(out, other_index[0], new_top_layer)
number_added = 0
for other_layer_num, s_layer in enumerate(out.copy()):
if (s_layer['top_p'] - p) * (p - s_layer['bot_p']) > 0:
# Found a slowness layer with the other wave type that
# contains the new slowness sample.
top_layer = np.array([(
s_layer['top_p'], s_layer['top_depth'], p,
bullen_depth_for(s_layer, p, self.radius_of_planet))],
dtype=SlownessLayer)
bot_layer = (p, top_layer['bot_depth'],
s_layer['bot_p'], s_layer['bot_depth'])
out[other_layer_num+number_added] = _flat(bot_layer)
out = np.insert(out, other_layer_num+number_added, top_layer)
# Fix critical layers since we have added a slowness layer.
_fix_critical_depths(critical_depths,
other_layer_num, not is_p_wave)
# Skip next layer as it was just added: achieved by slicing
# the list iterator.
number_added += 1
return out