Polish KRM model code and fix bug with non-swimbladdered fish example datasets

Author: gavinmacaulay

src/echosms/krmdata.py

@@ -14,7 +14,7 @@
 
 @dataclass
 class KRMshape():
-    """KRM shape and property storage.
+    """KRM shape and property class.
 
     Attributes
     ----------
@@ -25,36 +25,58 @@ class KRMshape():
     w :
         Width of the shape [m].
     z_U :
-        Distance from the fish axis to the upper surface of the shape [m].
+        Distance from the axis to the upper surface of the shape [m].
     z_L :
-        Distance from the fish axis to the lower surface of the shape [m].
+        Distance from the axis to the lower surface of the shape [m].
     c :
         Sound speed in the shape [m/s].
     rho :
-        Density in the shape [kg/m³].
+        Density of the shape material [kg/m³].
     """
 
     boundary: str  # 'soft' (aka swimbladder) or 'fluid' (aka body)
     x: np.ndarray
     w: np.ndarray
     z_U: np.ndarray
     z_L: np.ndarray
-    c: float = None
-    rho: float = None
+    c: float
+    rho: float
+
+    def volume(self) -> float:
+        """Volume of the shape.
+
+        Returns
+        -------
+        :
+            The volume of the shape [m³].
+        """
+        thickness = np.diff(self.x)
+        thickness = np.append(thickness, thickness[1])
+        return np.sum(np.pi * (self.z_U - self.z_L) * self.w * thickness)
+
+    def length(self) -> float:
+        """Length of the shape.
+
+        Returns
+        -------
+        :
+            The length of the shape [m].
+        """
+        return self.x[-1] - self.x[0]
 
 
 @dataclass
-class KRMfish():
-    """KRM fish and swimbladder model shapes.
+class KRMorganism():
+    """KRM body and swimbladder model shape.
 
     Attributes
     ----------
     name :
-        A name for the fish.
+        A name for the organism.
     source :
-        A link to or description of the source of the fish data.
+        A link to or description of the source of the organism data.
     shapes :
-        The shapes that make up the fish.
+        The shapes that make up the organism.
     """
 
     name: str
@@ -63,7 +85,7 @@ class KRMfish():
 
 
 class KRMdata():
-    """Provides example fish datasets for the KRM model."""
+    """Example datasets for the KRM model."""
 
     def __init__(self):
         # Load in the NOAA KRM shapes data
@@ -74,15 +96,17 @@ def __init__(self):
             except tomllib.TOMLDecodeError as e:
                 raise SyntaxError(f'Error while parsing file "{self.defs_filename.name}"') from e
 
-        # Put the shapes into a dict of KRMfish(). Use some default values for sound speed and
+        # Put the shapes into a dict of KRMorganism(). Use some default values for sound speed and
         # density
         self.krm_models = {}
         for s in shapes['shape']:
             body = KRMshape('fluid', np.array(s['x_b']), np.array(s['w_b']),
-                            np.array(s['z_bU']), np.array(s['z_bL']), 1570, 1070)
+                            np.array(s['z_bU']), np.array(s['z_bL']),
+                            s['body_c'], s['body_rho'])
             swimbladder = KRMshape('soft', np.array(s['x_sb']), np.array(s['w_sb']),
-                                   np.array(s['z_sbU']), np.array(s['z_sbL']), 345, 1.24)
-            self.krm_models[s['name']] = KRMfish(s['name'], s['source'], [body, swimbladder])
+                                   np.array(s['z_sbU']), np.array(s['z_sbL']),
+                                   s['swimbladder_c'], s['swimbladder_rho'])
+            self.krm_models[s['name']] = KRMorganism(s['name'], s['source'], [body, swimbladder])
 
     def names(self):
         """Available KRM model names."""
@@ -95,12 +119,12 @@ def as_dict(self) -> dict:
         -------
         :
             All the KRM model shapes. The dataset name is the dict key and the value is an instance
-            of `KRMfish`.
+            of `KRMorganism`.
 
         """
         return self.krm_models
 
-    def model(self, name: str) -> KRMfish:
+    def model(self, name: str) -> KRMorganism:
         """KRM model shape with requested name.
 
         Parameters
@@ -111,7 +135,7 @@ def model(self, name: str) -> KRMfish:
         Returns
         -------
         :
-            An instance of `KRMfish` or None if there is no model with `name`.
+            An instance of `KRMorganism` or None if there is no model with `name`.
 
         """
         try:
@@ -121,7 +145,7 @@ def model(self, name: str) -> KRMfish:
 
     @staticmethod
     def ts(name: str) -> np.ndarray:
-        """KRM model ts with requested name.
+        """KRM model ts from model `name`.
 
         Parameters
         ----------
@@ -131,7 +155,8 @@ def ts(name: str) -> np.ndarray:
         Returns
         -------
         :
-            The TS (re 1 m²) for some default model parameters [dB] or None of no TS data exist.
+            The TS (re 1 m²) for some default model parameters [dB] or None if no TS data
+            are available.
 
         """
         # Sometimes there will be TS results for the model (available for testing of the

src/echosms/krmmodel.py

@@ -1,10 +1,12 @@
 """A class that provides the Kirchhoff ray mode scattering model."""
 
 from math import log10, pi, sqrt, cos, sin, radians
+from cmath import exp
 import numpy as np
+from scipy.special import j0, y0, jvp, yvp
 from .utils import wavenumber, as_dict
 from .scattermodelbase import ScatterModelBase
-from .krmdata import KRMfish, KRMshape
+from .krmdata import KRMshape
 
 
 def _u(x, z, theta):
@@ -54,7 +56,11 @@ def validate_parameters(self, params):
     def calculate_ts_single(self, medium_c, medium_rho, theta, f, bodies,
                             validate_parameters=True, **kwargs) -> float:
         """
-        Calculate the scatter using the kirchhoff ray mode model for one set of parameters.
+        Calculate the scatter using the Kirchhoff ray mode model for one set of parameters.
+
+        Warning
+        --------
+        The low _ka_ part of this model has not yet been verified to give correct results.
 
         Parameters
         ----------
@@ -68,11 +74,11 @@ def calculate_ts_single(self, medium_c, medium_rho, theta, f, bodies,
             conventions/#coordinate-systems) [°].
         f : float
             Frequency to calculate the scattering at [Hz].
-        bodies: KRMfish
+        bodies: KRMorganism
             The body shapes that make up the model. Currently, `bodies` should contain only two
             shapes, one of which should have a boundary of `fluid` (aka, the fish body) and the
-            other a boundary of `soft` (aka, the swimbladder). KRMfish.shapes[0] should be the
-            fish body and KRMfish.shapes[1] the swimbladder.
+            other a boundary of `soft` (aka, the swimbladder). KRMorganism.shapes[0] should be the
+            fish body and KRMorganism.shapes[1] the swimbladder.
         validate_parameters : bool
             Whether to validate the model parameters.
 
@@ -83,7 +89,7 @@ def calculate_ts_single(self, medium_c, medium_rho, theta, f, bodies,
 
         Notes
         -----
-        The class implements the code in Clay & Horne (1994) and when ka < 0.15 as per Clay (1992).
+        The class implements the code in Clay & Horne (1994) and when ka < 0.15 uses Clay (1992).
 
         References
         ----------
@@ -126,14 +132,15 @@ def calculate_ts_single(self, medium_c, medium_rho, theta, f, bodies,
         R_bc = (gp*hp-1) / (gp*hp+1)  # Eqn (9)
 
         # Equivalent radius of swimbladder (as per Part A of paper)
-        a_e = sqrt(self._volume(swimbladder)
-                   / (pi * (np.max(swimbladder.x) - np.min(swimbladder.x))))
-        a_e = 10
+        a_e = sqrt(swimbladder.volume() / (pi * swimbladder.length()))
 
         # Choose which modelling approach to use
         if k*a_e < 0.15:
             # Do the mode solution for the swimbladder (and ignore the body?)
-            mode_sl = self._mode_solution(swimbladder)
+            # TODO: need to check if it should be target or medium for the numerators
+            g = target_rho / swimbladder_rho
+            h = target_c / target_c
+            mode_sl = self._mode_solution(swimbladder, g, h, k, a_e, swimbladder.length(), theta)
             return 20*log10(abs(mode_sl))
 
         # Do the Kirchhoff-ray approximation for the swimbladder and body
@@ -142,13 +149,53 @@ def calculate_ts_single(self, medium_c, medium_rho, theta, f, bodies,
 
         return 20*log10(abs(soft_sl + fluid_sl))
 
-    def _volume(self, shape):
-        """Volume of the object."""
-        return 1.0
+    def _mode_solution(self, swimbladder: KRMshape, g: float, h: float,
+                       k: float, a: float, L_e: float, theta: float) -> float:
+        """Backscatter from a soft swimbladder at low ka.
+
+        Parameters
+        ----------
+        swimbladder :
+            The shape.
+        g :
+            Ratio of medium density over swimbladder density.
+        h :
+            Ratio of medium sound speed over swimbladder sound speed.
+        k :
+            The wavenumber in the medium surrounding the swimbladder.
+        a :
+            Equivalent radius of swimbladder [m].
+        L_e :
+            Equivalent length of swimbladder [m].
+        theta :
+            Pitch angle to calculate the scattering at, as per the echoSMs
+            [coordinate system](https://ices-tools-dev.github.io/echoSMs/
+            conventions/#coordinate-systems) [°].
+
+        Returns
+        -------
+        :
+            The scattering length [m].
+        """
+        # Note: equation references in this function are to Clay (1992)
+        if h == 0.0:
+            raise ValueError('Ratio of sound speeds (h) cannot be zero for low ka solution.')
+
+        # Chi is approximately this. More accurate equations are in Appendix B of Clay (1992)
+        chi = -pi/4  # Eqn (B10) and paragraph below that equation
+
+        ka = k*a
+        kca = ka/h
+
+        C_0 = (jvp(0, kca)*y0(ka) - g*h*yvp(0, ka)*j0(kca))\
+            / (jvp(0, kca)*j0(ka) - g*h*jvp(0, ka)*j0(kca))  # Eqn (A1) with m=0
+        b_0 = -1 / (1+1j*C_0)  # Also Eqn (A1)
+
+        delta = k*L_e*cos(theta)  # Eqn (4)
+
+        S_M = (exp(1j*(chi - pi/4)) * L_e)/pi * sin(delta)/delta * b_0  # Eqn (15)
 
-    def _mode_solution(self, swimbladder):
-        """Backscatter from a soft swimbladder at low ka."""
-        return -1.
+        return S_M
 
     def _soft_KA(self, swimbladder, k, k_b, R_bc, TwbTbw, theta):
         """Backscatter from a soft object using the Kirchhoff approximation."""