Add RC2001 deceleration correction and calcMass support

wmpl/MetSim/MetSimErosionCyTools.pyx

@@ -376,7 +376,7 @@ cpdef double decelerationRK4(double dt, double K, double m, double rho_atm, doub
 
 
 @cython.cdivision(True) 
-cpdef double luminousEfficiency(int lum_eff_type, double lum_eff, double vel, double mass):
+cpdef double luminousEfficiency(int lum_eff_type, double lum_eff, double vel, double mass, double v_init=-1.0):
     """ Compute the luminous efficiency of the given type, velocity, and mass.
 
     Arguments:
@@ -393,13 +393,19 @@ cpdef double luminousEfficiency(int lum_eff_type, double lum_eff, double vel, do
         lum_eff: [double] Value of the constant luminous efficiency (percent).
         vel: [double] Velocity (m/s).
         mass: [double] Mass (kg).
+        v_init: [double] Initial velocity (m/s). If <= 0, the
+            velocity-dependent deceleration correction term from
+            Revelle & Ceplecha (2001) is not applied. Default is -1.
+            The Revelle & Ceplecha (2001) deceleration correction contains
+            log(v_init - vel). For v_init = vel, the correction is set to zero.
+            Values with v_init < vel are inconsistent and will fail.
 
     Return:
         tau: [double] Luminous efficiency (ratio).
 
     """
 
-    cdef double c1, c2
+    cdef double c1, c2, lv, decel, dv
 
     # Constant luminous efficiency
     if lum_eff_type == 0:
@@ -423,15 +429,41 @@ cpdef double luminousEfficiency(int lum_eff_type, double lum_eff, double vel, do
             c1 = -2.670
             c2 = -4.674
 
+        lv = log(vel/1000.0)
+
+        decel = 0.0
+
+        if v_init > 0:
+            # The Revelle & Ceplecha (2001) deceleration correction contains
+            # log(v_init - vel). For v_init = vel, the correction is set to zero.
+            # Values with v_init < vel are physically inconsistent and will fail.
+            dv = (v_init - vel)/1000.0
+
+            if dv == 0:
+                decel = 0.0
+            else:
+                decel = 0.26*log(dv) + 0.0042*log(dv)**3
 
         # Slow meteoroids
         if vel < 25372:
-            return (exp(c1 - 10.307*log(vel/1000.0) + 9.781*log(vel/1000.0)**2 - 3.0414*log(vel/1000.0)**3 \
-                + 0.3213*log(vel/1000.0)**4 + 1.15*tanh(0.38*log(mass))))/100.0
+            return exp(
+                c1
+                - 10.307*lv
+                + 9.781*lv**2
+                - 3.0414*lv**3
+                + 0.3213*lv**4
+                + 1.15*tanh(0.38*log(mass))
+                + decel
+            )/100.0
 
         # Fast meteoroids
         else:
-            return (exp(c2 + log(vel/1000.0) + 1.15*tanh(0.38*log(mass))))/100.0
+            return exp(
+                c2
+                + lv
+                + 1.15*tanh(0.38*log(mass))
+                + decel
+            )/100.0
 
     # Borovicka et al. (2013) - Kosice
     elif lum_eff_type == 4:

wmpl/Utils/Physics.py

@@ -65,7 +65,6 @@ def dynamicPressure(lat, lon, height, jd, velocity, gamma=1.0):
     return dyn_pressure
 
 
-
 def dynamicMass(bulk_density, lat, lon, height, jd, velocity, decel, gamma=1.0, shape_factor=1.21):
     """ Calculate dynamic mass at the given point on meteor's trajectory. 
 
@@ -104,7 +103,6 @@ def dynamicMass(bulk_density, lat, lon, height, jd, velocity, decel, gamma=1.0,
     return dyn_mass
 
 
-
 def calcIntensity(mag_abs, P_0m=840.0):
     """ Calculate the radiated power (intensity) from absolute magnitudes.
 
@@ -158,8 +156,7 @@ def calcRadiatedEnergy(time, mag_abs, P_0m=840.0):
     return intens_int
 
 
-
-def calcMass(time, mag_abs, velocity, tau=0.007, P_0m=840.0, lum_eff_mass=-1):
+def calcMass(time, mag_abs, velocity, tau=0.007, P_0m=840.0, lum_eff_mass=-1, v_init=None):
     """ Calculates the mass of a meteoroid from the time and absolute magnitude.
 
     Arguments:
@@ -186,7 +183,15 @@ def calcMass(time, mag_abs, velocity, tau=0.007, P_0m=840.0, lum_eff_mass=-1):
         lum_eff_mass: [float] Meteoroid mass (kg) used to evaluate the mass-dependent luminous efficiency
             models. Ignored for a float tau. If -1 (default), the mass is iterated to self-consistency
             (compute mass -> recompute tau -> repeat until convergence).
-
+        v_init: [float or None] Pre-atmospheric velocity (m/s) used by the ReVelle & Ceplecha (2001)
+            deceleration correction in the luminous efficiency calculation. Default is None.
+            - None: the deceleration correction is not applied.
+            - -1: use the maximum velocity value in the velocity array as the pre-atmospheric velocity.
+              If v_init=-1 and velocity is scalar, the deceleration correction is not applied.
+            - >0: use the supplied value as the pre-atmospheric velocity.
+            NOTE: the velocity array should represent a fitted deceleration curve (i.e. a physically
+            consistent monotonic velocity evolution) rather than individual noisy measurements.
+
     Return:
         mass: [float] Photometric mass of the meteoroid in kg.
 
@@ -200,14 +205,43 @@ def calcMass(time, mag_abs, velocity, tau=0.007, P_0m=840.0, lum_eff_mass=-1):
     #
     # For constant velocity and luminous efficiency this reduces to:
     #   M = (2/(tau*v^2))*integral(I dt)
+    #
+    # If tau is constant but velocity varies along the trajectory:
+    #   M = (2/tau)*integral(I/v^2 dt)
 
-    # A fixed luminous efficiency was given: compute the mass using that constant value.
-    if not isinstance(tau, str) and not isinstance(tau, (bool, int)):
-
-        intens_int = calcRadiatedEnergy(time, mag_abs, P_0m=P_0m)
-
-        return (2.0/(tau*velocity**2))*intens_int
+    if len(time) != len(mag_abs):
+        raise ValueError(
+            "time and mag_abs must have the same length."
+        )
+
+    # If a velocity array is given, it must have the same length as the
+    # photometric measurements.
+    if not np.isscalar(velocity) and len(velocity) != len(time):
+        raise ValueError(
+            "Velocity array length ({:d}) does not match the number of "
+            "time/magnitude measurements ({:d}).".format(
+                len(velocity), len(time)
+            )
+        )
+
+    if lum_eff_mass < 0 and lum_eff_mass != -1:
+        raise ValueError(
+            "lum_eff_mass must be >= 0 or -1 for self-consistent iteration."
+        )
 
+    # Constant luminous efficiency (tau given directly as a ratio).
+    if not isinstance(tau, str) and not isinstance(tau, (bool, int)):
+        # Constant velocity
+        if np.isscalar(velocity):
+            intens_int = calcRadiatedEnergy(time, mag_abs, P_0m=P_0m)
+            return (2.0/(tau*velocity**2))*intens_int
+
+        # Velocity array
+        else:
+            intens = calcIntensity(mag_abs, P_0m=P_0m)
+            vel_arr = np.asarray(velocity, dtype=float)
+            dm_dt = 2.0*intens/(tau*vel_arr**2)
+            return simpson(dm_dt, x=time)
 
     # Velocity-dependent model: resolve the lum_eff_type code.
     if isinstance(tau, str):
@@ -226,14 +260,26 @@ def calcMass(time, mag_abs, velocity, tau=0.007, P_0m=840.0, lum_eff_mass=-1):
     # Radiated power at every measurement point
     intens = calcIntensity(mag_abs, P_0m=P_0m)
 
-    # Work with array velocity so a single code path handles scalar and array input
+
+    # Convert velocity to an array for the velocity-dependent luminous efficiency calculations.
     vel_arr = np.atleast_1d(np.asarray(velocity, dtype=float))
 
+    if v_init is None:
+        # Disable the ReVelle & Ceplecha deceleration correction
+        v_init_eff = -1.0
+    elif v_init == -1:
+        if np.isscalar(velocity):
+            v_init_eff = -1.0
+        else:
+            v_init_eff = float(np.max(vel_arr))
+    else:
+        v_init_eff = float(v_init)
+
     def _mass_for(mass_assumed):
         """ Compute the photometric mass given an assumed meteoroid mass for the lum. eff. models. """
 
         # Vectorize the scalar Cython call over the velocity points (lum_eff arg only used for type 0)
-        tau_arr = np.array([luminousEfficiency(lum_eff_type, 0.7, v, mass_assumed) for v in vel_arr])
+        tau_arr = np.array([luminousEfficiency(lum_eff_type, 0.7, v, mass_assumed, v_init_eff) for v in vel_arr])
 
         # Integrate the instantaneous mass-loss rate dm/dt = 2*I/(tau*v^2)
         dm_dt = 2.0*intens/(tau_arr*vel_arr**2)
@@ -262,11 +308,15 @@ def _mass_for(mass_assumed):
             mass = mass_new
             break
         mass = mass_new
+    else:
+        print(
+            "WARNING: luminous efficiency self-consistency iteration "
+            "did not converge after 50 iterations."
+        )
 
     return mass
 
 
-
 def calcKB(lat, lon, ht_beg, jd, v, zenith_angle, gmn_correction=False):
     """ Calculate the Celpecha (1958) KB parameter.