Source code for disba._sensitivity

from collections import namedtuple

import numpy as np

from ._base import BaseSensitivity
from ._common import ifunc, ipar
from ._cps import srfker96, swegn96

__all__ = [
    "SensitivityKernel",
    "PhaseSensitivity",
    "GroupSensitivity",
    "EllipticitySensitivity",
]


SensitivityKernel = namedtuple(
    "SensitivityKernel",
    ("depth", "kernel", "period", "velocity", "mode", "wave", "type", "parameter"),
)


[docs]class PhaseSensitivity(BaseSensitivity): def __init__( self, thickness, velocity_p, velocity_s, density, algorithm="dunkin", dc=0.005, dp=0.025, ): """ Phase velocity sensitivity kernel class. Parameters ---------- thickness : array_like Layer thickness (in km). velocity_p : array_like Layer P-wave velocity (in km/s). velocity_s : array_like Layer S-wave velocity (in km/s). density : array_like Layer density (in g/cm3). algorithm : str {'dunkin', 'fast-delta'}, optional, default 'dunkin' Algorithm to use for computation of Rayleigh-wave dispersion: - 'dunkin': Dunkin's matrix (adapted from surf96), - 'fast-delta': fast delta matrix (after Buchen and Ben-Hador, 1996). dc : scalar, optional, default 0.005 Phase velocity increment for root finding. dp : scalar, optional, default 0.025 Parameter increment (%) for numerical partial derivatives. """ super().__init__(thickness, velocity_p, velocity_s, density, algorithm, dc, dp)
[docs] def __call__(self, t, mode=0, wave="rayleigh", parameter="velocity_s"): """ Calculate phase velocity sensitivity kernel for a given period and parameter. Parameters ---------- t : scalar Period (in s). mode : int, optional, default 0 Mode number (0 if fundamental). wave : str {'love', 'rayleigh'}, optional, default 'rayleigh' Wave type. parameter : str {'thickness', 'velocity_p', 'velocity_s', 'density'}, optional, default 'velocity_s' Parameter with respect to which sensitivity kernel is calculated. Returns ------- :class:`disba.SensitivityKernel` Sensitivity kernel as a namedtuple (depth, kernel, period, velocity, mode, wave, type, parameter). """ c1, kernel = srfker96( t, self._thickness, self._velocity_p, self._velocity_s, self._density, mode=mode, itype=0, ifunc=ifunc[self._algorithm][wave], ipar=ipar[parameter], dc=self._dc, dp=self._dp, ) return SensitivityKernel( np.insert(self._thickness.cumsum()[:-1], 0, 0.0), kernel, t, c1, mode, wave, "phase", parameter, )
[docs]class GroupSensitivity(BaseSensitivity): def __init__( self, thickness, velocity_p, velocity_s, density, algorithm="dunkin", dc=0.005, dt=0.025, dp=0.025, ): """ Phase velocity sensitivity kernel class. Parameters ---------- thickness : array_like Layer thickness (in km). velocity_p : array_like Layer P-wave velocity (in km/s). velocity_s : array_like Layer S-wave velocity (in km/s). density : array_like Layer density (in g/cm3). algorithm : str {'dunkin', 'fast-delta'}, optional, default 'dunkin' Algorithm to use for computation of Rayleigh-wave dispersion: - 'dunkin': Dunkin's matrix (adapted from surf96), - 'fast-delta': fast delta matrix (after Buchen and Ben-Hador, 1996). dc : scalar, optional, default 0.005 Phase velocity increment for root finding. dt : scalar, optional, default 0.025 Frequency increment (%) for calculating group velocity. dp : scalar, optional, default 0.025 Parameter increment (%) for numerical partial derivatives. """ if not isinstance(dt, float): raise TypeError() super().__init__(thickness, velocity_p, velocity_s, density, algorithm, dc, dp) self._dt = dt
[docs] def __call__(self, t, mode=0, wave="rayleigh", parameter="velocity_s"): """ Calculate group velocity sensitivity kernel for a given period and parameter. Parameters ---------- t : scalar Period (in s). mode : int, optional, default 0 Mode number (0 if fundamental). wave : str {'love', 'rayleigh'}, optional, default 'rayleigh' Wave type. parameter : str {'thickness', 'velocity_p', 'velocity_s', 'density'}, optional, default 'velocity_s' Parameter with respect to which sensitivity kernel is calculated. Returns ------- :class:`disba.SensitivityKernel` Sensitivity kernel as a namedtuple (depth, kernel, period, velocity, mode, wave, type, parameter). """ c1, kernel = srfker96( t, self._thickness, self._velocity_p, self._velocity_s, self._density, mode=mode, itype=1, ifunc=ifunc[self._algorithm][wave], ipar=ipar[parameter], dc=self._dc, dt=self._dt, dp=self._dp, ) return SensitivityKernel( np.insert(self._thickness.cumsum()[:-1], 0, 0.0), kernel, t, c1, mode, wave, "group", parameter, )
@property def dt(self): """Return frequency increment (%) for calculating group velocity.""" return self._dt
[docs]class EllipticitySensitivity(BaseSensitivity): def __init__( self, thickness, velocity_p, velocity_s, density, algorithm="dunkin", dc=0.005, dp=0.025, ): """ Rayleigh-wave ellipticity sensitivity kernel class. Parameters ---------- thickness : array_like Layer thickness (in km). velocity_p : array_like Layer P-wave velocity (in km/s). velocity_s : array_like Layer S-wave velocity (in km/s). density : array_like Layer density (in g/cm3). algorithm : str {'dunkin', 'fast-delta'}, optional, default 'dunkin' Algorithm to use for computation of Rayleigh-wave dispersion: - 'dunkin': Dunkin's matrix (adapted from surf96), - 'fast-delta': fast delta matrix (after Buchen and Ben-Hador, 1996). dc : scalar, optional, default 0.005 Phase velocity increment for root finding. dp : scalar, optional, default 0.025 Parameter increment (%) for numerical partial derivatives. """ super().__init__(thickness, velocity_p, velocity_s, density, algorithm, dc, dp)
[docs] def __call__(self, t, mode=0, parameter="velocity_s"): """ Calculate Rayleigh-wave ellipticity sensitivity kernel for a given period and parameter. Parameters ---------- t : scalar Period (in s). mode : int, optional, default 0 Mode number (0 if fundamental). parameter : str {'thickness', 'velocity_p', 'velocity_s', 'density'}, optional, default 'velocity_s' Parameter with respect to which sensitivity kernel is calculated. Returns ------- :class:`disba.SensitivityKernel` Sensitivity kernel as a namedtuple (depth, kernel, period, velocity, mode, wave, type, parameter). """ # Reference ellipticity ell1 = self._ellipticity(t, mode) # Initialize kernel mmax = len(self._thickness) kernel = np.empty(mmax) # Loop over layers fac = 1.0 + self._dp par = getattr(self, parameter) for i in range(mmax): tmp = par[i] par[i] /= fac ell2 = self._ellipticity(t, mode) kernel[i] = (ell2 - ell1) / (par[i] - tmp) par[i] *= fac return SensitivityKernel( np.insert(self._thickness.cumsum()[:-1], 0, 0.0), kernel, t, None, mode, "rayleigh", "ellipticity", parameter, )
def _ellipticity(self, t, mode): """Compute Rayleigh-wave ellipticity for input period and mode.""" eig = swegn96( t, self._thickness, self._velocity_p, self._velocity_s, self._density, mode, ifunc[self._algorithm]["rayleigh"], self._dc, )[:, :2] return eig[0, 0] / eig[0, 1]