Merge branch 'main' of https://github.com/ices-tools-dev/echoSMs

docs/other_software.md

@@ -4,7 +4,7 @@ Other software that provides source code for acoustic scattering models of relev
 
 - [acousticTS](https://github.com/brandynlucca/acousticTS): R code for calculating scattering using the DCM, DWBA, SDWBA, SDWBA_curved, KRM, MSS model, as well as that of calibration spheres.
 - [Coupled BEM acoustic](https://github.com/elavia/coupled_bem_acoustic): Julia code that calculates the TS of three-dimensional shapes with an included object (e.g., a swimbladder).
-- [scatmod](https://github.com/SvenGastauer/scatmod): Open source acoustic scattering models for fisheries acoustics. Python and R code for fluid spheres and calibration spheres.
+- [scatmod](https://github.com/SvenGastauer/scatmod): Open source acoustic scattering models for fisheries acoustics. Python and R code for fluid spheres.
 - [FishAcoustics](https://github.com/gavinmacaulay/FishAcoustics): Contains a Python module that implements the phase-tracking DWBA model.
 - [KRM Model](https://www.fisheries.noaa.gov/data-tools/krm-model): A web page that uses the KRM model to estimate the TS of predefined or user-supplied shapes over a range of input parameters.
 - [KRMr](https://github.com/SvenGastauer/KRMr): KRM model for fish in R. 

src/echosms/__init__.py

@@ -3,9 +3,10 @@
 from .scattermodelbase import ScatterModelBase
 from .benchmarkdata import BenchmarkData
 from .referencemodels import ReferenceModels
+from .esmodel import ESModel
 from .mssmodel import MSSModel
 from .psmsmodel import PSMSModel
 from .dcmmodel import DCMModel
 
 __all__ = ['ScatterModelBase', 'BenchmarkData', 'ReferenceModels', 'MSSModel', 'PSMSModel',
-           'DCMModel', 'k', 'eta', 'h1', 'as_dataframe', 'as_dataarray']
+           'DCMModel', 'ESModel', 'k', 'eta', 'h1', 'as_dataframe', 'as_dataarray']

src/echosms/dcmmodel.py

@@ -78,9 +78,8 @@ def calculate_ts_single(self, medium_c, medium_rho, a, b, theta, f, boundary_typ
         K = k(medium_c, f) * sin(theta_rad)
         Ka = K*a
 
-        m = range(30)  # this needs to vary with f
+        m = range(30)  # TODO this needs to vary with f
 
-        # Some code varies with model type.
         match boundary_type:
             case 'fixed rigid':
                 series = list(map(lambda m: (-1)**m * eta(m)*(jvp(m, Ka) / h1vp(m, Ka)), m))

src/echosms/esmodel.py

@@ -0,0 +1,74 @@
+"""A class that provides the modal series solution scattering model."""
+
+import numpy as np
+from math import log10
+# from mapply.mapply import mapply
+# import swifter
+from scipy.special import spherical_jn, spherical_yn
+from .utils import h1, k
+from .scattermodelbase import ScatterModelBase
+
+
+class ESModel(ScatterModelBase):
+    """Elastic sphere (ES) scattering model.
+
+    This class calculates acoustic scatter from spheres and shells with various
+    boundary conditions, as listed in the ``boundary_types`` class attribute.
+    """
+
+    def __init__(self):
+        super().__init__()
+        self.long_name = 'elastic sphere'
+        self.short_name = 'es'
+        self.analytical_type = 'exact'
+        self.boundary_types = ['elastic sphere']
+        self.shapes = ['sphere']
+        self.max_ka = 20  # [1]
+
+    def calculate_ts_single(self, medium_c, medium_rho, diameter, theta, f, boundary_type,
+                            target_longitudal_c, target_transverse_c, target_rho,
+                            **kwargs) -> float:
+        """
+        Calculate the scatter from an elastic sphere for one set of parameters.
+
+        Parameters
+        ----------
+        medium_c : float
+            Sound speed in the fluid medium surrounding the target [m/s].
+        medium_rho : float
+            Density of the fluid medium surrounding the target [kg/m³].
+        diameter : float
+            Diameter of the sphere [m].
+        theta : float
+            Pitch angle(s) to calculate the scattering at [°]. An angle of 0 is head on,
+            90 is dorsal, and 180 is tail on.
+        f : float
+            Frequencies to calculate the scattering at [Hz].
+        boundary_type : str, optional
+            The boundary type. Supported types are given in the boundary_types class variable.
+        target_longitudal_c : float
+            Longitudal sound speed in the material inside the sphere [m/s].
+        target_transverse_c : float
+            Transverse sound speed in the material inside the sphere [m/s].
+        target_rho : float
+            Density of the material inside the sphere [kg/m³].
+
+        Returns
+        -------
+        : float
+            The target strength (re 1 m²) of the sphere [dB].
+
+        Notes
+        -----
+        The class implements the code in [1].
+
+        References
+        ----------
+        MacLennan, D. N. (1981). The Theory of Solid Spheres as Sonar Calibration Targets
+        Scottish Fisheries Research Report Number 22. Department of Agriculture and Fisheries
+        for Scotland.
+        """
+        k0 = k(medium_c, f)
+        ka = k0*diameter
+
+        return -1.0

src/echosms/mssmodel.py

@@ -87,7 +87,6 @@ def calculate_ts_single(self, medium_c, medium_rho, a, theta, f, boundary_type,
         ka = k0*a
         n = np.arange(0, round(ka+20))
 
-        # Some code varies with model type.
         match boundary_type:
             case 'fixed rigid':
                 A = list(map(lambda x: -spherical_jn(x, ka, True) / h1(x, ka, True), n))

src/echosms/referencemodels.py

@@ -13,6 +13,11 @@ class ReferenceModels:
     ----------
     definitions : dict
         A dict representation of the ``target definitions.toml`` file.
+
+    Raises
+    ------
+    KeyError
+        If the ``target definitions.toml`` file has multiple target entries with the same name.
     """
 
     def __init__(self):
@@ -23,7 +28,12 @@ def __init__(self):
         with open(self.defs_filename, 'rb') as f:
             self.definitions = tomllib.load(f)
 
-        # check that the names parameter is unique across all models
+        # Flag duplicate target names
+        pda = pd.Series(self.names())
+        duplicates = list(pda[pda.duplicated()])
+        if duplicates:
+            raise KeyError(f'The "{self.defs_filename.name}" file has multiple targets '
+                           f'with the same name: '+', '.join(duplicates))
 
         # Substitute parameters names in the target section by the values of those
         # parameters.
@@ -40,7 +50,7 @@ def names(self):
         Returns
         -------
         : iterable of str
-            All model names in the ``target definitions/toml`` file.
+            All model names in the ``target definitions.toml`` file.
         """
         return [n['name'] for n in self.definitions['target']]
 
@@ -55,15 +65,14 @@ def specification(self, name):
         Returns
         -------
         : dict
-            The model definitions for the requested model or ``None`` if no model with that name.
+            The model definitions for the requested model or an empty set if no model
+            with that name.
         """
-        models = pd.DataFrame(self.definitions['target'])
-        m = models.loc[models['name'] == name]
-        if len(m) == 1:
-            m.dropna(axis=1, how='all', inplace=True)
-            return m.iloc[0].to_dict()
+        s = [t for t in self.definitions['target'] if t['name'] == name]
+        if not s:
+            return s
 
-        return None
+        return s[0]
 
     def parameters(self, name):
         """Model parameters for a particular model.
@@ -79,18 +88,19 @@ def parameters(self, name):
         Returns
         -------
         : dict
-            The model parameters for the requested model or ``None`` if no model with that name.
+            The model parameters for the requested model or an empty set if no model with that name.
 
         """
         s = self.specification(name)
 
-        if s is None:
-            return None
+        if not s:
+            return []
 
         # Remove the entries that are not parameters
-        p = s
+        p = s.copy()
         del p['name']
         del p['shape']
         del p['description']
         del p['source']
+        del p['benchmark_model']
         return p

src/echosms/resources/target definitions.toml

@@ -14,6 +14,7 @@ sound_speed = "m/s"
 radius = "m"
 length = "m"
 thickness = "m"
+diameter = "m"
 
 # Model parameters (add more as needed).
 #
@@ -52,8 +53,8 @@ example3_medium_sound_speed_dist = {type="arange", start=1490, stop=1500, step=1
 # For each reference target, include a [[target]] section below.
 #
 # Where sensible, use names from the parameter section above for specifying 
-# target parameters (e.g., the sound speed, density, size, etc). Substitution of the
-# parameter names with their values from the parameters section should be done in the 
+# target parameters (e.g., the sound speed, density, size, etc). # The 'name' entry should be unique.
+# Substitution of the parameter names with their values from the parameters section should be done in the 
 # code that reads this file.
 #
 # The source item in each section should be used to identify the source of the model configuration.
@@ -69,6 +70,7 @@ a = "sphere_radius"
 medium_rho = "medium_density"
 medium_c = "medium_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "mss"
 
 [[target]]
 name = "pressure release sphere"
@@ -79,6 +81,7 @@ a = "sphere_radius"
 medium_rho = "medium_density"
 medium_c = "medium_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "mss"
 
 [[target]]
 name = "gas filled sphere"
@@ -91,6 +94,7 @@ medium_c = "medium_sound_speed"
 target_rho = "target_gas_density"
 target_c = "target_gas_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "mss"
 
 [[target]]
 name = "weakly scattering sphere"
@@ -103,6 +107,7 @@ medium_c = "medium_sound_speed"
 target_rho = "target_weakly_density"
 target_c = "target_weakly_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "mss"
 
 [[target]]
 name = "spherical fluid shell with pressure release interior"
@@ -116,6 +121,7 @@ medium_c = "medium_sound_speed"
 shell_rho = "target_shell_pressure_release_density"
 shell_c = "target_shell_pressure_release_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "mss"
 
 [[target]]
 name = "spherical fluid shell with gas interior"
@@ -131,6 +137,7 @@ shell_c = "target_shell_gas_filled_sound_speed"
 target_rho = "target_gas_density"
 target_c = "target_gas_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "mss"
 
 [[target]]
 name = "spherical fluid shell with weakly scattering interior"
@@ -146,6 +153,7 @@ shell_c = "target_shell_weakly_sound_speed"
 target_rho = "target_interior_weakly_density"
 target_c = "target_interior_weakly_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "mss"
 
 [[target]]
 name = "fixed rigid prolate spheroid"
@@ -157,6 +165,7 @@ b = "prolate_spheroid_minor_axis_radius"
 medium_rho = "medium_density"
 medium_c = "medium_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "psms"
 
 [[target]]
 name = "pressure release prolate spheroid"
@@ -168,6 +177,7 @@ b = "prolate_spheroid_minor_axis_radius"
 medium_rho = "medium_density"
 medium_c = "medium_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "psms"
 
 [[target]]
 name = "gas filled prolate spheroid"
@@ -181,6 +191,7 @@ medium_c = "medium_sound_speed"
 target_rho = "target_gas_density"
 target_c = "target_gas_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "psms"
 
 [[target]]
 name = "weakly scattering prolate spheroid"
@@ -194,6 +205,7 @@ medium_c = "medium_sound_speed"
 target_rho = "target_weakly_density"
 target_c = "target_weakly_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "psms"
 
 [[target]]
 name = "fixed rigid finite cylinder"
@@ -205,6 +217,7 @@ b = "finite_cylinder_length"
 medium_rho = "medium_density"
 medium_c = "medium_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "dcm"
 
 [[target]]
 name = "pressure release finite cylinder"
@@ -216,6 +229,7 @@ b = "finite_cylinder_length"
 medium_rho = "medium_density"
 medium_c = "medium_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "dcm"
 
 [[target]]
 name = "gas filled finite cylinder"
@@ -229,6 +243,7 @@ medium_c = "medium_sound_speed"
 target_rho = "target_gas_density"
 target_c = "target_gas_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "dcm"
 
 [[target]]
 name = "weakly scattering finite cylinder"
@@ -242,3 +257,32 @@ medium_c = "medium_sound_speed"
 target_rho = "target_weakly_density"
 target_c = "target_weakly_sound_speed"
 source = "https://doi.org/10.1121/1.4937607"
+benchmark_model = "dcm"
+
+[[target]]
+name = "WC38.1 calibration sphere"
+shape = "sphere"
+boundary_type = "elastic sphere"
+description = "A 38.1 mm diameter tungsten carbide sphere with 6% cobalt binder"
+diameter = 0.0381 # diameter
+medium_rho = "medium_density"
+medium_c = "medium_c"
+target_rho = 14900
+target_longitudal_c = 6853
+target_transverse_c = 4171
+source = " "
+benchmark_model = "es"
+
+[[target]]
+name = "Cu60 calibration sphere"
+shape = "sphere"
+boundary_type = "elastic sphere"
+description = "A 60 mm diameter copper sphere"
+diameter = 0.060 # diameter
+medium_rho = "medium_density"
+medium_c = "medium_sound_speed"
+target_rho = 8947
+target_longitudal_c = 4760
+target_transverse_c = 2288.5
+source = " "
+benchmark_model = "es"

src/echosms/scattermodelbase.py

@@ -68,7 +68,6 @@ def calculate_ts(self, data, multiprocess=False):
         elif isinstance(data, pd.DataFrame):
             pass
         elif isinstance(data, xr.DataArray):
-            # For the moment just convert DataArrays into DataFrames
             data = data.to_dataframe().reset_index()
         else:
             raise ValueError(f'Data type of {type(data)} is not supported'
@@ -81,7 +80,6 @@ def calculate_ts(self, data, multiprocess=False):
             # ts = df.swifter.apply(self.__ts_helper, axis=1)
             ts = data.apply(self.__ts_helper, axis=1)
         else:  # this uses just one CPU
-            # ts = data.apply(self.__ts_helper, axis=1)
             ts = data.apply(self.__ts_helper, axis=1)
 
         return ts.to_numpy()  # TODO - return data type that matches the input data type

src/echosms/utils.py

@@ -44,14 +44,12 @@ def as_dataarray(params: dict) -> xr.DataArray:
         in the input dict.
 
     """
-    # Convert scalars to iterables so xarray is happier later on
+    # Convert scalars to iterables so xarray is happy
     for k, v in params.items():
         if not isinstance(v, Iterable) or isinstance(v, str):
             params[k] = [v]
 
-    # Lengths of each parameter array
     sz = [len(v) for k, v in params.items()]
-    # Create the DataArray
     return xr.DataArray(data=np.full(sz, np.nan), coords=params, name='ts')