Fully working elastic sphere model (#24)

* Fully working elastic sphere model

* typo and spelling fixes

src/echosms/__init__.py

@@ -1,5 +1,5 @@
 """Setup the public API for echoSMs."""
-from .utils import k, eta, h1, as_dataframe, as_dataarray
+from .utils import k, eta, h1, spherical_jnpp, as_dataframe, as_dataarray
 from .scattermodelbase import ScatterModelBase
 from .benchmarkdata import BenchmarkData
 from .referencemodels import ReferenceModels
@@ -9,4 +9,5 @@
 from .dcmmodel import DCMModel
 
 __all__ = ['ScatterModelBase', 'BenchmarkData', 'ReferenceModels', 'MSSModel', 'PSMSModel',
-           'DCMModel', 'ESModel', 'k', 'eta', 'h1', 'as_dataframe', 'as_dataarray']
+           'DCMModel', 'ESModel', 'k', 'eta', 'h1', 'spherical_jnpp',
+           'as_dataframe', 'as_dataarray']

src/echosms/esmodel.py

@@ -1,18 +1,17 @@
-"""A class that provides the modal series solution scattering model."""
+"""A class that provides the elastic scattering model."""
 
 import numpy as np
-from math import log10
-# from mapply.mapply import mapply
-# import swifter
+from math import log10, sin, atan
+from cmath import exp
 from scipy.special import spherical_jn, spherical_yn
-from .utils import h1, k
+from .utils import k, spherical_jnpp
 from .scattermodelbase import ScatterModelBase
 
 
 class ESModel(ScatterModelBase):
     """Elastic sphere (ES) scattering model.
 
-    This class calculates acoustic scatter from elastic spheres.
+    This class calculates acoustic backscatter from elastic spheres.
     """
 
     def __init__(self):
@@ -24,27 +23,24 @@ def __init__(self):
         self.shapes = ['sphere']
         self.max_ka = 20  # [1]
 
-    def calculate_ts_single(self, medium_c, medium_rho, diameter, theta, f,
-                            target_longitudal_c, target_transverse_c, target_rho,
+    def calculate_ts_single(self, medium_c, medium_rho, a, f,
+                            target_longitudinal_c, target_transverse_c, target_rho,
                             **kwargs) -> float:
         """
-        Calculate the scatter from an elastic sphere for one set of parameters.
+        Calculate the backscatter from an elastic sphere for one set of parameters.
 
         Parameters
         ----------
         medium_c : float
             Sound speed in the fluid medium surrounding the sphere [m/s].
         medium_rho : float
             Density of the fluid medium surrounding the sphere [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.
+        a : float
+            Radius of the sphere [m].
         f : float
-            Frequencies to calculate the scattering at [Hz].
-        target_longitudal_c : float
-            Longitudal sound speed in the material inside the sphere [m/s].
+            Frequency to calculate the scattering at [Hz].
+        target_longitudinal_c : float
+            Longitudinal 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
@@ -65,7 +61,26 @@ def calculate_ts_single(self, medium_c, medium_rho, diameter, theta, f,
         Scottish Fisheries Research Report Number 22. Department of Agriculture and Fisheries
         for Scotland.
         """
-        k0 = k(medium_c, f)
-        ka = k0*diameter
+        q = k(medium_c, f)*a
+        q1 = q*medium_c/target_longitudinal_c
+        q2 = q*medium_c/target_transverse_c
+        alpha = 2. * (target_rho/medium_rho) * (target_transverse_c/medium_c)**2
+        beta = (target_rho/medium_rho) * (target_longitudinal_c/medium_c)**2 - alpha
 
-        return -1.0
+        n = range(20)
+
+        # Use n instead of l (ell) because l looks like 1.
+        def S(n):
+            A2 = (n**2 + n-2) * spherical_jn(n, q2) + q2**2 * spherical_jnpp(n, q2)
+            A1 = 2*n*(n+1) * (q1*spherical_jn(n, q1, True) - spherical_jn(n, q1))
+            B2 = A2*q1**2 * (beta*spherical_jn(n, q1) - alpha*spherical_jnpp(n, q1))\
+                - A1*alpha * (spherical_jn(n, q2) - q2*spherical_jn(n, q2, True))
+            B1 = q * (A2*q1*spherical_jn(n, q1, True) - A1*spherical_jn(n, q2))
+            eta_n = atan(-(B2*spherical_jn(n, q, True) - B1*spherical_jn(n, q))
+                         / (B2*spherical_yn(n, q, True) - B1*spherical_yn(n, q)))
+
+            return (-1)**n * (2*n+1) * sin(eta_n) * exp(1j*eta_n)
+
+        f_inf = -2.0/q * np.sum(list(map(S, n)))
+
+        return 10*log10(a**2 * abs(f_inf)**2 / 4.0)

src/echosms/resources/target definitions.toml

@@ -264,11 +264,11 @@ 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
+a = 0.01905 # radius
 medium_rho = "medium_density"
-medium_c = "medium_c"
+medium_c = "medium_sound_speed"
 target_rho = 14900
-target_longitudal_c = 6853
+target_longitudinal_c = 6853
 target_transverse_c = 4171
 source = " "
 benchmark_model = "es"
@@ -278,11 +278,11 @@ name = "Cu60 calibration sphere"
 shape = "sphere"
 boundary_type = "elastic sphere"
 description = "A 60 mm diameter copper sphere"
-diameter = 0.060 # diameter
+a = 0.030 # radius
 medium_rho = "medium_density"
 medium_c = "medium_sound_speed"
 target_rho = 8947
-target_longitudal_c = 4760
+target_longitudinal_c = 4760
 target_transverse_c = 2288.5
 source = " "
 benchmark_model = "es"

src/echosms/utils.py

@@ -132,3 +132,22 @@ def h1(n: int, z: float, derivative=False) -> complex:
         return spherical_jn(n, z) + 1j*spherical_yn(n, z)
     else:
         return -h1(n+1, z) + (n/z) * h1(n, z)
+
+
+def spherical_jnpp(n: int, z: float) -> float:
+    """Second derivative of the spherical Bessel function.
+
+    Parameters
+    ----------
+    n :
+        Order (n ≥ 0)
+    z :
+        Argument of the Bessel function.
+
+    Returns
+    -------
+    :
+        The second derivative of the spherical Bessel function.
+
+    """
+    return 1./z**2 * ((n**2-n-z**2)*spherical_jn(n, z) + 2.*z*spherical_jn(n+1, z))

src/example_code.py

@@ -83,7 +83,7 @@
     plt.suptitle(name[0])
     plt.show()
 
-    # %% ###############################################################################################
+# %% ###############################################################################################
 # Run the benchmark models and compare to the angle-varying benchmark results.
 models = [('fixed rigid finite cylinder', 'Cylinder_Rigid'),
           ('pressure release finite cylinder', 'Cylinder_PressureRelease'),
@@ -134,6 +134,30 @@
     plt.suptitle(name[0])
     plt.show()
 
+# %% ###############################################################################################
+# Use the ES model on a calibration sphere
+name = 'WC38.1 calibration sphere'
+p = rm.parameters(name)
+p['f'] = np.arange(10, 100, 0.05) * 1e3  # [kHz]
+
+es = ESModel()
+ts = es.calculate_ts(p)
+
+plt.plot(p['f']*1e-3, ts)
+plt.xlabel('Freq [kHz]')
+plt.ylabel('TS re 1 m$^2$ [dB] ')
+plt.title(name)
+plt.show()
+
+# Can readily modify the parameters for a different sphere
+p['a'] = 0.012/2
+ts = es.calculate_ts(p)
+plt.plot(p['f']*1e-3, ts)
+plt.xlabel('Freq [kHz]')
+plt.ylabel('TS re 1 m$^2$ [dB] ')
+plt.title('WC12 calibration sphere')
+plt.show()
+
 # %% ###############################################################################################
 # Some other ways to run models.
 
@@ -189,18 +213,3 @@
 
 # and it can be inserted into params_xa
 # TODO once the data is returned in an appropriate form
-
-# %% ###############################################################################################
-# Use the ES model on a calibration sphere
-p = rm.parameters('Cu60 calibration sphere')
-p['theta'] = 90
-p['f'] = np.arange(10, 400, 0.5) * 1e3  # [kHz]
-
-es = ESModel()
-ts = es.calculate_ts(p)
-
-plt.plot(p['f']*1e-3, ts)
-
-# Can readily modify the parameters for a different sphere
-p['a'] = 0.063
-ts = es.calculate_ts(p)