from __future__ import annotations
import math
from collections.abc import Sequence
from dataclasses import dataclass
import numpy as np
[docs]
@dataclass(frozen=True)
class BioLayer:
"""Biological layer for Fortran ``AttenMod`` volume ``B`` (dB/km resonance parameters)."""
z1: float
z2: float
f0: float
q: float
a0: float
@dataclass(frozen=True)
class FrancoisGarrisonParams:
"""Parameters for Francois-Garrison volume attenuation."""
T: float
S: float
pH: float
z_bar: float
def thorp_db_per_km(f_khz: float | np.ndarray) -> float | np.ndarray:
"""Thorp volume attenuation (dB/km), JKPS Eq. 1.34.
Port of Matlab `thorp.m`.
"""
f2 = np.asarray(f_khz) ** 2
alpha = 0.0033 + 0.11 * f2 / (1 + f2) + 44 * f2 / (4100 + f2) + 0.0003 * f2
return float(alpha) if np.isscalar(f_khz) else alpha
def franc_garr_db_per_km(
f_khz: float | np.ndarray,
*,
T: float,
S: float,
pH: float,
z: float,
) -> float | np.ndarray:
"""Francois-Garrison attenuation (dB/km).
Port of Matlab `franc_garr.m`.
"""
f = np.asarray(f_khz, dtype=float)
c = 1412 + 3.21 * T + 1.19 * S + 0.0167 * z
# Boric acid contribution
A1 = 8.86 / c * 10 ** (0.78 * pH - 5)
P1 = 1.0
f1 = 2.8 * math.sqrt(S / 35) * 10 ** (4 - 1245 / (T + 273))
# Magnesium sulfate contribution
A2 = 21.44 * S / c * (1 + 0.025 * T)
P2 = 1 - 1.37e-4 * z + 6.2e-9 * z**2
f2 = 8.17 * 10 ** (8 - 1990 / (T + 273)) / (1 + 0.0018 * (S - 35))
# Viscosity
P3 = 1 - 3.83e-5 * z + 4.9e-10 * z**2
if T < 20:
A3 = 4.937e-4 - 2.59e-5 * T + 9.11e-7 * T**2 - 1.5e-8 * T**3
else:
A3 = 3.964e-4 - 1.146e-5 * T + 1.45e-7 * T**2 - 6.5e-10 * T**3
alpha = (
A1 * P1 * (f1 * f**2) / (f1**2 + f**2)
+ A2 * P2 * (f2 * f**2) / (f2**2 + f**2)
+ A3 * P3 * f**2
)
return float(alpha) if np.isscalar(f_khz) else alpha
def crci(
c: float,
alpha: float,
*,
freq_hz: float,
atten_unit: str,
fg: FrancoisGarrisonParams | None = None,
depth_m: float = 0.0,
freq0_hz: float | None = None,
power_law_beta: float | None = None,
freq_transition_hz: float | None = None,
bio_layers: Sequence[BioLayer] | None = None,
) -> complex:
"""Convert real wave speed and attenuation to a complex wave speed.
Port of Matlab `crci.m`, extended to match Fortran ``AttenMod:CRCI``.
- First character of ``atten_unit``: ``N``, ``M``, ``m``, ``F``, ``W``, ``Q``, ``L``
(see Fortran ``AttenMod.f90``; ``m`` is dB/m with frequency power law).
- Second character: ``T``, ``F``, ``B``, or blank (Thorp, Francois-Garrison, biological).
``depth_m`` is the profile depth (m) used only for biological ``B`` volume loss.
For ``m``, supply ``freq0_hz``, ``power_law_beta``, and ``freq_transition_hz`` (Fortran
``freq0``, ``beta``, ``fT``).
"""
if len(atten_unit) < 1:
raise ValueError("atten_unit must be at least 1 char")
omega = 2.0 * math.pi * freq_hz
# Convert to nepers/m
unit = atten_unit[0:1]
alphaT = 0.0
if unit == "N":
alphaT = alpha
elif unit == "M":
alphaT = alpha / 8.6858896
elif unit == "m":
if freq0_hz is None or power_law_beta is None or freq_transition_hz is None:
raise ValueError(
"attenuation unit 'm' requires freq0_hz, power_law_beta, freq_transition_hz"
)
if freq0_hz == 0.0:
raise ValueError("freq0_hz must be non-zero for attenuation unit 'm'")
alphaT = alpha / 8.6858896
if freq_hz < freq_transition_hz:
alphaT *= (freq_hz / freq0_hz) ** power_law_beta
else:
alphaT *= (freq_hz / freq0_hz) * (freq_transition_hz / freq0_hz) ** (
power_law_beta - 1.0
)
elif unit == "F":
alphaT = alpha * freq_hz / 8685.8896
elif unit == "W":
alphaT = 0.0 if c == 0.0 else alpha * freq_hz / (8.6858896 * c)
elif unit == "Q":
alphaT = 0.0 if (c * alpha) == 0.0 else omega / (2.0 * c * alpha)
elif unit == "L":
# Fortran AttenMod.f90: loss parameter -> Nepers/m before c^2/omega scaling
alphaT = 0.0 if c == 0.0 else alpha * omega / c
else:
raise ValueError(f"unknown attenuation unit: {unit!r}")
# Optional volume attenuation
vol = atten_unit[1:2] if len(atten_unit) >= 2 else ""
if vol == "T":
th = thorp_db_per_km(freq_hz / 1000.0) / 8685.8896 # nepers/m
alphaT += float(th)
elif vol == "F":
if fg is None:
raise ValueError("Francois-Garrison enabled but fg params not provided")
fg_db_km = franc_garr_db_per_km(freq_hz / 1000.0, T=fg.T, S=fg.S, pH=fg.pH, z=fg.z_bar)
alphaT += float(fg_db_km) / 8685.8896 # nepers/m
elif vol == "B":
if bio_layers is None:
raise ValueError("biological volume attenuation (B) requires bio_layers")
f_sq = max(freq_hz * freq_hz, 1e-30)
for layer in bio_layers:
if layer.z1 <= depth_m <= layer.z2:
q2 = max(layer.q * layer.q, 1e-30)
den = (1.0 - layer.f0**2 / f_sq) ** 2 + 1.0 / q2
a_db_km = layer.a0 / den
alphaT += float(a_db_km) / 8685.8896
# Convert nepers/m to equivalent imaginary sound speed
alphaT = alphaT * (c**2) / omega
return complex(c, alphaT)