# -*- coding: utf-8 -*-
"""
Misc functionality.
"""
import inspect
from pathlib import Path
import numpy as np
ROOT = Path(Path(inspect.getfile(inspect.currentframe())).resolve()).parent
[docs]
def parse_phase_list(phase_list):
"""
Takes a list of phases, returns a list of individual phases. Performs e.g.
replacing e.g. ``"ttall"`` with the relevant phases.
"""
phase_names = []
for phase_name in phase_list:
phase_names += get_phase_names(phase_name)
# Remove duplicates.
return sorted(list(set(phase_names)))
[docs]
def get_phase_names(phase_name):
"""
Called by parse_phase_list to replace e.g. ttall with the relevant phases.
"""
lphase = phase_name.lower()
names = []
if lphase in ("ttp", "tts", "ttbasic", "tts+", "ttp+", "ttall"):
if lphase in ("ttp", "ttp+", "ttbasic", "ttall"):
names.extend(["p", "P", "Pn", "Pdiff", "PKP", "PKiKP", "PKIKP"])
if lphase in ("tts", "tts+", "ttbasic", "ttall"):
names.extend(["s", "S", "Sn", "Sdiff", "SKS", "SKIKS"])
if lphase in ("ttp+", "ttbasic", "ttall"):
names.extend(["PcP", "pP", "pPdiff", "pPKP", "pPKIKP", "pPKiKP",
"sP", "sPdiff", "sPKP", "sPKIKP", "sPKiKP"])
if lphase in ("tts+", "ttbasic", "ttall"):
names.extend(["sS", "sSdiff", "sSKS", "sSKIKS", "ScS", "pS",
"pSdiff", "pSKS", "pSKIKS"])
if lphase in ("ttbasic", "ttall"):
names.extend(["ScP", "SKP", "SKIKP", "PKKP", "PKIKKIKP", "SKKP",
"SKIKKIKP", "PP", "PKPPKP", "PKIKPPKIKP"])
if lphase == "ttall":
names.extend(["SKiKP", "PP", "ScS", "PcS", "PKS", "PKIKS", "PKKS",
"PKIKKIKS", "SKKS", "SKIKKIKS", "SKSSKS",
"SKIKSSKIKS", "SS", "SP", "PS"])
else:
names.append(phase_name)
return names
[docs]
def split_ray_path(path, model):
"""
Split and label ray path according to type of wave.
:param path: Path taken by ray.
:type path: :class:`~numpy.ndarray`
(dtype = :class:`~obspy.taup.helper_classes.TimeDist`)
:param model: The model used to calculate the ray path.
:type model: :class:`obspy.taup.tau_model.TauModel`
:returns: A list of paths and a list of wave types.
Wave types are either "p", "s", or "diff".
:rtype: list[:class:`~numpy.ndarray`]
and list[str]
The ray path is split at all discontinuities in the model
and labelled according to wave type.
"""
# Get discontinuity depths in the model in km
discs = model.s_mod.v_mod.get_discontinuity_depths()[:-1]
# Split path at discontinuity depths
depths = path["depth"]
is_disc = np.isin(depths, discs)
is_disc[0] = False # Don't split on first point
idx = np.where(is_disc)[0]
# Split ray path, including start and end points in each segment
splitted = np.split(path, idx)
paths = [
np.append(s, splitted[i + 1][0]) for i, s in enumerate(splitted[:-1])
]
# Classify the waves as p, s, or diff
wave_types = [_classify_path(p, model) for p in paths]
return paths, wave_types
[docs]
def _expected_delay_time(ray_param, depth0, depth1, wave, model):
"""
Expected delay time between two depths for a given wave type (p or s).
"""
# Convert depths to radii
radius0 = model.radius_of_planet - depth0
radius1 = model.radius_of_planet - depth1
# Velocity model from TauModel
v_mod = model.s_mod.v_mod
# Get velocities
if depth1 >= depth0:
v0 = v_mod.evaluate_below(depth0, wave)[0]
v1 = v_mod.evaluate_above(depth1, wave)[0]
else:
v0 = v_mod.evaluate_above(depth0, wave)[0]
v1 = v_mod.evaluate_below(depth1, wave)[0]
# Calculate time for segment if velocity non-zero
# - if velocity zero then return zero time
if v0 > 0.0:
eta0 = radius0 / v0
eta1 = radius1 / v1
def vertical_slowness(eta, p):
y = eta**2 - p**2
return np.sqrt(y * (y > 0)) # in s
n0 = vertical_slowness(eta0, ray_param)
n1 = vertical_slowness(eta1, ray_param)
if ray_param == 0.0:
return 0.5 * ((1.0 / v0) + (1.0 / v1)) * abs(radius1 - radius0)
return 0.5 * (n0 + n1) * abs(np.log(radius1 / radius0))
return 0.0
[docs]
def _classify_path(path, model):
"""
Determine whether we have a p or s-wave path by comparing delay times.
"""
# Examine just the first two points near the shallowest part of the path
if path[0]["depth"] < path[-1]["depth"]:
point0 = path[0]
point1 = path[1]
else:
point0 = path[-2]
point1 = path[-1]
# Ray parameter
ray_param = point0["p"]
# Depths
depth0 = point0["depth"]
depth1 = point1["depth"]
# If no change in depth then this is a diffracted/head wave segment
if depth0 == depth1:
return "diff"
# Delay time for this segment from ray path
travel_time = point1["time"] - point0["time"]
distance = abs(point0["dist"] - point1["dist"])
delay_time = travel_time - ray_param * distance
# Get the expected delay time for each wave type
delay_p = _expected_delay_time(ray_param, depth0, depth1, "p", model)
delay_s = _expected_delay_time(ray_param, depth0, depth1, "s", model)
# Difference between predictions and given delay times
error_p = (delay_p / delay_time) - 1.0
error_s = (delay_s / delay_time) - 1.0
# Check which wave type matches the given delay time the best
if abs(error_p) < abs(error_s):
return "p"
return "s"
[docs]
def _flat(list_or_tuple):
"""
Return a flat version of the input tuple. Intention is to avoid numpy
single-item arrays with ndim >0. Reason being that newer numpy has
deprecated automatic conversion to scalars from these types of arrays and
this would result in errors in some future numpy releases.
Seems like we rely on this automatic conversion to scalars a lot. This very
likely could be done faster and more elegant, not sure if this is impacting
execution speed of our taup routines.
"""
def _item_or_return(x):
try:
return x.item()
except AttributeError:
return x
return tuple(_item_or_return(x) for x in list_or_tuple)