Much improved PSMSModel calculations

src/echosms/psmsmodel.py

@@ -1,7 +1,7 @@
 """A class that provides the prolate spheroidal modal series scattering model."""
 
 import numpy as np
-from scipy.integrate import quad
+from math import factorial
 from .scattermodelbase import ScatterModelBase
 from .utils import pro_rad1, pro_rad2, pro_ang1, wavenumber, Neumann
 
@@ -65,26 +65,20 @@ def calculate_ts_single(self, medium_c, medium_rho, a, b, theta, f, boundary_typ
             “Prediction of krill target strength by liquid prolate spheroid
             model,” Fish. Sci., 60, 261–265.
         """
-        match boundary_type:
-            case 'pressure release' | 'fluid filled':
-                pass
-            case 'fixed rigid':
-                raise ValueError(f'Model type "{boundary_type}" has not yet been implemented '
-                                 f'for the {self.long_name} model.')
-            case _:
-                raise ValueError(f'The {self.long_name} model does not support '
-                                 f'a model type of "{boundary_type}".')
-
-        if boundary_type == 'fluid filled':
-            hc = target_c / medium_c
-            rh = target_rho / medium_rho
+        if boundary_type not in self.boundary_types:
+            raise ValueError(f'The {self.long_name} model does not support '
+                             f'a model type of "{boundary_type}".')
 
         xiw = (1.0 - (b/a)**2)**(-.5)
         q = a/xiw  # semi-focal length
 
         kw = wavenumber(medium_c, f)
         hw = kw*q
 
+        if boundary_type == 'fluid filled':
+            g = target_rho / medium_rho
+            ht = wavenumber(target_c, f)*q
+
         # Phi, the port/starboard angle is fixed for this code
         phi_inc = np.pi  # incident direction
         phi_sca = np.pi + phi_inc  # scattered direction
@@ -100,46 +94,44 @@ def calculate_ts_single(self, medium_c, medium_rho, a, b, theta, f, boundary_typ
         for m in range(m_max+1):
             epsilon_m = Neumann(m)
             for n in range(m, n_max+1):
-                Smn_inc, _ = pro_ang1(m, n, hw, np.cos(theta_inc))
-                Smn_sca, _ = pro_ang1(m, n, hw, np.cos(theta_sca))
+                Smn_inc, _ = pro_ang1(m, n, hw, np.cos(theta_inc), norm=True)
+                Smn_sca, _ = pro_ang1(m, n, hw, np.cos(theta_sca), norm=True)
+                # The Meixner-Schäfke normalisation scheme for the angular function of the first
+                # kind. Refer to eqn 21.7.11 in Abramowitz, M., and Stegun, I. A. (1964).
+                # Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables
+                # (Dover, New York), 10th ed.
+                n_mn = 2/(2*n+1) * factorial(n+m) / factorial(n-m)
+                # since n_mn and the angular functions can be very large for large m,
+                # structure things to reduce roundoff errors.
+                ss = (Smn_inc / n_mn) * Smn_sca
+
                 match boundary_type:
                     case 'fluid filled':
-                        r_type1A, dr_type1A = pro_rad1(m, n, hw, xiw)
-                        r_type2A, dr_type2A = pro_rad2(m, n, hw, xiw)
-                        r_type1B, dr_type1B = pro_rad1(m, n, hw/hc, xiw)
-
-                        eeA = r_type1A - rh*r_type1B/dr_type1B*dr_type1A
-                        eeB = eeA + 1j*(r_type2A - rh*r_type1B/dr_type1B*dr_type2A)
-                        Amn = -eeA/eeB  # Furusawa (1988) Eq. 5 p 15
+                        # Note: we can implement the simplier equations if hw is similar to ht,
+                        # but that only applies to weakly scattering conditions. The gas-filled
+                        # condition does not meet that, so we have two paths here. The simplified
+                        # equations are quicker, so it is worth to do.
+                        R1w, dR1w = pro_rad1(m, n, hw, xiw)
+                        R1t, dR1t = pro_rad1(m, n, ht, xiw)
+                        R2w, dR2w = pro_rad2(m, n, hw, xiw)
+
+                        R3w = R1w + 1j*R2w
+                        dR3w = dR1w + 1j*dR2w
+                        if abs((hw - ht)/ht) <= 0.1:
+                            E1 = R1w - g * R1t / dR1t * dR1w
+                            E3 = R3w - g * R1t / dR1t * dR3w
+                            Amn = -E1/E3
+                        else:
+                            Amn = 1.0
                     case 'pressure release':
-                        r_type1, _ = pro_rad1(m, n, hw, xiw)
-                        r_type2, _ = pro_rad2(m, n, hw, xiw)
-                        Amn = -r_type1/(r_type1 + 1j*r_type2)
+                        R1w, _ = pro_rad1(m, n, hw, xiw)
+                        R2w, _ = pro_rad2(m, n, hw, xiw)
+                        Amn = -R1w/(R1w + 1j*R2w)
                     case 'fixed rigid':
-                        pass  # see eqn of (3) of Furusawa, 1988
-
-                # This definition of the norm of S is in Yeh (1967), and is equation
-                # 21.7.11 in Abramowitz & Stegun (10th printing) as the
-                # Meixner-Schãfke normalisation scheme. Note that the RHS of
-                # 21.7.11 doesn't give correct results compared to doing the actual
-                # integration.
-                #
-                # Yeh, C. (1967). "Scattering of Acoustic Waves by a Penetrable
-                #   Prolate Spheroid I Liquid Prolate Spheroid," J. Acoust. Soc.
-                #   Am. 42, 518-521.
-                #
-                # Abramowitz, M., and Stegun, I. A. (1964). Handbook of
-                #   Mathematical Functions with Formulas, Graphs, and Mathematical
-                #   Tables (Dover, New York), 10th ed.
-
-                n_mn = quad(PSMSModel._aswfa2, -1.0, 1.0, args=(m, n, hw), epsrel=1e-5)
-
-                f_sc += epsilon_m / n_mn[0]\
-                    * Smn_inc * Amn * Smn_sca * np.cos(m*(phi_sca - phi_inc))
+                        _, dR1w = pro_rad1(m, n, hw, xiw)
+                        _, dR2w = pro_rad2(m, n, hw, xiw)
+                        Amn = -dR1w/(dR1w + 1j*dR2w)
 
-        return 20*np.log10(np.abs(-2j / kw * f_sc))
+                f_sc += epsilon_m * ss * Amn * np.cos(m*(phi_sca - phi_inc))
 
-    @staticmethod
-    def _aswfa2(x, m, n, h0):
-        """Just pro_ang1()**2, but with x as the first parameter for use with `quad()`."""
-        return pro_ang1(m, n, h0, x)[0]**2
+        return 20*np.log10(np.abs(-2j / kw * f_sc))

src/example_code.py

@@ -57,9 +57,10 @@ def plot_compare(f1, ts1, label1, f2, ts2, label2, title):
           ('pressure release finite cylinder', 'Cylinder_PressureRelease'),
           ('gas filled finite cylinder', 'Cylinder_Gas'),
           ('weakly scattering finite cylinder', 'Cylinder_WeaklyScattering'),
-          # ('fixed rigid prolate spheroid', 'ProlateSpheroid_Rigid'),
-          # ('pressure release prolate spheroid', 'ProlateSpheroid_PressureRelease'),
-          # ('gas filled prolate spheroid', 'ProlateSpheroid_Gas'),
+          ('fixed rigid prolate spheroid', 'ProlateSpheroid_Rigid'),
+          ('pressure release prolate spheroid', 'ProlateSpheroid_PressureRelease'),
+          ('gas filled prolate spheroid', 'ProlateSpheroid_Gas'),
+          # weakly scattering takes a while to run, so leave it out for the moment
           # ('weakly scattering prolate spheroid', 'ProlateSpheroid_WeaklyScattering'),
           ]
 
@@ -79,24 +80,37 @@ def plot_compare(f1, ts1, label1, f2, ts2, label2, title):
         case _:
             pass
 
-    # Add frequencies that have finite TS values and incident angle to the model parameters
+    # Add frequencies that have non nan benchmark TS values to the model parameters
     m['f'] = bmf['Frequency_kHz'][~np.isnan(bmf[bm_name])]*1e3  # [Hz]
+
+    # No benchmark TS is available for this model, so add some sensible frequencies in
+    if bm_name == 'ProlateSpheroid_Gas':
+        m['f'] = np.arange(12, 82, 2)*1e3
+
     m['theta'] = 90.0
 
-    # and run these
+    # and run the models
     ts = mod.calculate_ts(m)
 
-    # The finite benchmark TS values
+    # Get the benchmark TS values
     bm_ts = bmf[bm_name][~np.isnan(bmf[bm_name])]
 
+    # Cope with there being no benchmark TS for this model
+    if bm_name == 'ProlateSpheroid_Gas':
+        bm_ts = m['f'] * np.nan
+
     plot_compare(m['f'], ts, s['benchmark_model'], m['f'], bm_ts, 'Benchmark', name)
 
 # %% ###############################################################################################
 # 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'),
           ('gas filled finite cylinder', 'Cylinder_Gas'),
-          ('weakly scattering finite cylinder', 'Cylinder_WeaklyScattering')]
+          ('weakly scattering finite cylinder', 'Cylinder_WeaklyScattering'),
+          ('fixed rigid prolate spheroid', 'ProlateSpheroid_Rigid'),
+          ('pressure release prolate spheroid', 'ProlateSpheroid_PressureRelease'),
+          ('gas filled prolate spheroid', 'ProlateSpheroid_Gas'),
+          ('weakly scattering prolate spheroid', 'ProlateSpheroid_WeaklyScattering'),]
 
 for name, bm_name in models:
     # Get the model parameters used in Jech et al. (2015) for a particular model.