# -*- coding: utf-8 -*-
"""
Functions acting on slowness layers.
"""
import math
import numpy as np
from .c_wrappers import clibtau
from .helper_classes import SlownessLayer, SlownessModelError
from .velocity_layer import (evaluate_velocity_at_bottom,
evaluate_velocity_at_top)
[docs]
def bullen_radial_slowness(layer, p, radius_of_planet, check=True):
"""
Calculate time and distance increments of a spherical ray.
The time and distance (in radians) increments accumulated by a ray of
spherical ray parameter p when passing through this layer. Note that this
gives half of the true range and time increments since there will be both
an upgoing and a downgoing path. Here we use the Mohorovicic or Bullen
law: p=A*r^B
The ``layer`` and ``p`` parameters must be either 0-D, or both of the same
shape.
:param layer: The layer(s) in which to calculate the increments.
:type layer: :class:`~numpy.ndarray`,
dtype = :class:`obspy.taup.helper_classes.SlownessLayer`
:param p: The spherical ray paramater to use for calculation, in s/km.
:type p: :class:`~numpy.ndarray`, dtype = :class:`float`
:param radius_of_planet: The radius of the planet to use, in km.
:type radius_of_planet: float
:param check: Check that the calculated results are not invalid. This
check may be disabled if the layers requested are expected not to
include the specified ray.
:type check: bool
:returns: Time (in s) and distance (in rad) increments.
:rtype: tuple(:class:`~numpy.ndarray`)
"""
ldim = np.ndim(layer)
pdim = np.ndim(p)
if ldim == 1 and pdim == 0:
time = np.zeros(shape=layer.shape, dtype=np.float_)
dist = np.zeros(shape=layer.shape, dtype=np.float_)
elif ldim == 0 and pdim == 1:
time = np.zeros(shape=p.shape, dtype=np.float_)
dist = np.zeros(shape=p.shape, dtype=np.float_)
elif ldim == pdim and (ldim == 0 or layer.shape == p.shape):
time = np.zeros(shape=layer.shape, dtype=np.float_)
dist = np.zeros(shape=layer.shape, dtype=np.float_)
else:
raise TypeError('Either layer or p must be 0D, or they must have '
'the same shape.')
length = len(time)
if isinstance(p, np.float_):
p = p * np.ones(length, dtype=np.float_)
clibtau.bullen_radial_slowness_inner_loop(
layer, p, time, dist, radius_of_planet, length)
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("timedist.time or .dist < 0 or Nan")
return time, dist
[docs]
def bullen_depth_for(layer, ray_param, radius_of_planet, check=True):
"""
Finds the depth for a ray parameter within this layer.
Uses a Bullen interpolant, Ar^B. Special case for ``bot_p == 0`` or
``bot_depth == radius_of_planet`` as these cause division by 0; use linear
interpolation in this case.
The ``layer`` and ``ray_param`` parameters must be either 0-D, or both of
the same shape.
:param layer: The layer(s) to check.
:type layer: :class:`~numpy.ndarray` (shape = (1, ), dtype =
:class:`obspy.taup.helper_classes.SlownessLayer`)
:param ray_param: The ray parameter(s) to use for calculation, in s/km.
:type ray_param: float
:param radius_of_planet: The radius (in km) of the planet to use.
:type radius_of_planet: float
:returns: The depth (in km) for the specified ray parameter.
:rtype: float
"""
ldim = np.ndim(layer)
pdim = np.ndim(ray_param)
if ldim == 1 and pdim == 0:
ray_param = ray_param * np.ones(layer.shape, dtype=np.float_)
depth = np.zeros(shape=layer.shape, dtype=np.float_)
elif ldim == 0 and pdim == 1:
layer = layer * np.ones(ray_param.shape, dtype=SlownessLayer)
depth = np.zeros(shape=ray_param.shape, dtype=np.float_)
elif ldim == pdim and (ldim == 0 or layer.shape == ray_param.shape):
if ldim == 0:
# Make array-like to work with NumPy < 1.9.
layer = np.array([layer], dtype=SlownessLayer)
ray_param = np.array([ray_param])
depth = np.zeros(shape=layer.shape, dtype=np.float_)
else:
raise TypeError('Either layer or ray_param must be 0D, or they must '
'have the same shape.')
valid = (layer['top_p'] - ray_param) * (ray_param - layer['bot_p']) >= 0
if not check or np.all(valid):
leftover = np.ones_like(depth, dtype=np.bool_)
# Easy cases for 0 thickness layer, or ray parameter found at
# top or bottom.
mask = layer['top_depth'] == layer['bot_depth']
depth[mask] = layer['bot_depth'][mask]
leftover &= ~mask
mask = leftover & (layer['top_p'] == ray_param)
depth[mask] = layer['top_depth'][mask]
leftover &= ~mask
mask = leftover & (layer['bot_p'] == ray_param)
depth[mask] = layer['bot_depth'][mask]
leftover &= ~mask
mask = leftover & (
(layer['bot_p'] != 0) & (layer['bot_depth'] != radius_of_planet))
if np.any(mask):
top_p_mask = layer['top_p'][mask]
bot_p_mask = layer['bot_p'][mask]
top_depth_mask = layer['top_depth'][mask]
bot_depth_mask = layer['bot_depth'][mask]
ray_param_mask = ray_param[mask]
b = np.divide(np.log(top_p_mask / bot_p_mask),
np.log((radius_of_planet - top_depth_mask) /
(radius_of_planet - bot_depth_mask)))
with np.errstate(over='ignore'):
denom = np.power(radius_of_planet - top_depth_mask, b)
a = np.divide(top_p_mask, denom)
temp_depth = np.empty_like(a)
mask2 = (a != 0) & (b != 0)
temp_depth[mask2] = radius_of_planet - np.exp(
1.0 / b[mask2] * np.log(np.divide(ray_param_mask[mask2],
a[mask2])))
# or equivalent (maybe better stability?):
# tempDepth = radius_of_planet - math.pow(ray_param_mask/A, 1/B)
# Overflow. Use linear interpolation.
temp_depth[~mask2] = (
(bot_depth_mask[~mask2] - top_depth_mask[~mask2]) /
(bot_p_mask[~mask2] - top_p_mask[~mask2]) *
(ray_param_mask[~mask2] - top_p_mask[~mask2])
) + top_depth_mask[~mask2]
# Check if slightly outside layer due to rounding or
# numerical instability:
mask2 = ((top_depth_mask > temp_depth) &
(temp_depth > top_depth_mask - 0.000001))
temp_depth[mask2] = top_depth_mask[mask2]
mask2 = ((bot_depth_mask < temp_depth) &
(temp_depth < bot_depth_mask + 0.000001))
temp_depth[mask2] = bot_depth_mask[mask2]
mask2 = ((temp_depth < 0) | np.isnan(temp_depth) |
np.isinf(temp_depth) |
(temp_depth < top_depth_mask) |
(temp_depth > bot_depth_mask))
# Numerical instability in power law calculation? Try a
# linear interpolation if the layer is small (<5km).
small_layer = bot_depth_mask[mask2] - top_depth_mask[mask2] > 5
if np.any(small_layer):
if check:
raise SlownessModelError(
"Calculated depth is outside layer, negative, or NaN.")
else:
temp_depth[mask2][small_layer] = np.nan
linear = (
(bot_depth_mask[mask2] - top_depth_mask[mask2]) /
(bot_p_mask[mask2] - top_p_mask[mask2]) *
(ray_param_mask[mask2] - top_p_mask[mask2])
) + top_depth_mask[mask2]
outside_layer = small_layer & (
linear < 0 | np.isnan(linear) | np.isinf(linear))
if np.any(outside_layer):
if check:
raise SlownessModelError(
"Calculated depth is outside layer, negative, or NaN.")
else:
temp_depth[mask2][outside_layer] = np.nan
temp_depth[mask2] = linear
# Check for tempDepth just above top_depth or below bottomDepth.
mask2 = ((temp_depth < top_depth_mask) &
(top_depth_mask - temp_depth < 1e-10))
temp_depth[mask2] = top_depth_mask[mask2]
mask2 = ((temp_depth > bot_depth_mask) &
(temp_depth - bot_depth_mask < 1e-10))
temp_depth[mask2] = bot_depth_mask[mask2]
depth[mask] = temp_depth
leftover &= ~mask
# Special case for the centre of the planet, since Ar^B might
# blow up at r = 0.
mask = leftover & (layer['top_p'] != layer['bot_p'])
depth[mask] = (layer['bot_depth'][mask] +
(ray_param[mask] - layer['bot_p'][mask]) *
(layer['top_depth'][mask] - layer['bot_depth'][mask]) /
(layer['top_p'][mask] - layer['bot_p'][mask]))
leftover &= ~mask
# weird case, return bot_depth??
depth[leftover] = layer['bot_depth'][leftover]
# Make sure invalid cases are left out
depth[~valid] = np.nan
# NumPy < 1.9 compatibility.
if ldim == 0 and pdim == 0:
return depth[0]
else:
return depth
else:
raise SlownessModelError(
"Ray parameter is not contained within this slowness layer.")
[docs]
def evaluate_at_bullen(layer, depth, radius_of_planet):
"""
Find the slowness at the given depth.
Note that this method assumes a Bullen type of slowness interpolation,
i.e., p(r) = a*r^b. This will produce results consistent with a tau model
that uses this interpolant, but it may differ slightly from going directly
to the velocity model. Also, if the tau model is generated using another
interpolant, linear for instance, then the result may not be consistent
with the tau model.
:param layer: The layer to use for the calculation.
:type layer: :class:`numpy.ndarray`,
dtype = :class:`obspy.taup.helper_classes.SlownessLayer`
:param depth: The depth (in km) to use for the calculation. It must be
contained within the provided ``layer`` or else results are undefined.
:type depth: float
:param radius_of_planet: The radius of the planet to use, in km.
:type radius_of_planet: float
"""
top_p = layer['top_p']
bot_p = layer['bot_p']
top_depth = layer['top_depth']
bot_depth = layer['bot_depth']
# Could do some safeguard asserts...
assert not bot_depth > radius_of_planet
assert not (top_depth - depth) * (depth - bot_depth) < 0
if depth == top_depth:
return top_p
elif depth == bot_depth:
return bot_p
else:
try:
# The power law calculation has some stability issues with very
# small layers. Catch them here to trigger the fallback
# computations later on.
with np.errstate(all="raise"):
b = math.log(top_p / bot_p) / \
math.log((radius_of_planet - top_depth) /
(radius_of_planet - bot_depth))
a_denominator = pow((radius_of_planet - top_depth), b)
a = top_p / a_denominator
answer = a * pow((radius_of_planet - depth), b)
except FloatingPointError:
answer = np.nan
if answer < 0 or math.isnan(answer) or math.isinf(answer):
# numerical instability in power law calculation???
# try a linear interpolation if the layer is small ( <2 km)
# or if denominator of A is infinity as we probably overflowed
# the double in that case.
if bot_depth - top_depth < 2 \
or math.isinf(a_denominator) \
or bot_p == 0:
linear = (bot_p - top_p) / (bot_depth - top_depth) * \
(depth - top_depth) + top_p
if linear < 0 \
or math.isinf(linear) \
or math.isnan(linear):
pass
else:
return linear
raise SlownessModelError(
"Calculated Slowness is NaN or negative!")
return answer
[docs]
def create_from_vlayer(v_layer, is_p_wave, radius_of_planet,
is_spherical=True):
"""
Compute the slowness layer from a velocity layer.
:param v_layer: The velocity layer to convert.
:type v_layer: :class:`numpy.ndarray`,
dtype = :class:`obspy.taup.velocity_layer.VelocityLayer`
:param is_p_wave: Whether this velocity layer is for compressional/P
(``True``) or shear/S (``False``) waves.
:type is_p_wave: bool
:param radius_of_planet: The radius of the planet to use, in km.
:type radius_of_planet: float
:param is_spherical: Whether the model is spherical. Non-spherical models
are not currently supported.
:type is_spherical: bool
"""
ret = np.empty(shape=v_layer.shape, dtype=SlownessLayer)
ret['top_depth'] = v_layer['top_depth']
ret['bot_depth'] = v_layer['bot_depth']
wave_type = ('p' if is_p_wave else 's')
if is_spherical:
ret['top_p'] = (radius_of_planet - ret['top_depth']) / \
evaluate_velocity_at_top(v_layer, wave_type)
bot_depth = ret["bot_depth"]
bot_vel = evaluate_velocity_at_bottom(v_layer, wave_type)
if bot_depth.shape and bot_depth.size:
if bot_depth[-1] == radius_of_planet and bot_vel[-1] == 0.0:
bot_depth[-1] = 1.0
ret['bot_p'] = (radius_of_planet - bot_depth) / bot_vel
else:
raise NotImplementedError("no flat models yet")
return ret