diff --git a/wmpl/MetSim/MetSimErosionCyTools.pyx b/wmpl/MetSim/MetSimErosionCyTools.pyx index d7ec84bb..097f2f0d 100644 --- a/wmpl/MetSim/MetSimErosionCyTools.pyx +++ b/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: diff --git a/wmpl/Utils/Physics.py b/wmpl/Utils/Physics.py index 1b28c0f4..8beb7751 100644 --- a/wmpl/Utils/Physics.py +++ b/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.