Add new photolysis rates; reformat photo_timescale exception; light code reformatting.

Author: mkelley

CHANGES.rst

@@ -10,23 +10,35 @@ Update supported versions [#419]:
 New Features
 ------------
 
+sbpy.activity.gas
+^^^^^^^^^^^^^^^^^
+
+- Added active and quiet sun photolysis rates for H2O and CO2 from Huebner &
+  Mukherjee 2015 to `sbpy.activity.gas.photo_timescale`. [#425]
+
+
 sbpy.names
 ^^^^^^^^^^
+
 - Added functionality to `sbpy.Names.from_packed()` and
   `sbpy.Names.to_packed()` to handle new extended provisional designations
   to be implemented by the MPC in anticipation of higher asteroid discovery
   rates in the LSST survey era [#406]
 
+
 API Changes
 -----------
+
 - Deprecated `sbpy.ginga_plugins` in favor of using `sbpy-ginga` at
   https://github.com/NASA-Planetary-Science/sbpy-ginga [#413]
 
+
 sbpy.data.ephem
 ^^^^^^^^^^^^^^^
 
 - New command-line script: ``sbpy-ephem``. [#396]
 
+
 Other Changes and Additions
 ---------------------------
 
@@ -35,6 +47,7 @@ sbpy.activity.gas
 - Replaced calls to the deprecated function `scipy.integrate.romberg` with
   `scipy.integrate.quad`.  [#412]
 
+
 sbpy.names
 ^^^^^^^^^^
 - Fixed `sbpy.Names.to_packed()' to raise an error in cases of invalid
@@ -53,6 +66,7 @@ sbpy.names
   "P2015XN77"), which was previously being interpreted as a packed
   asteroid designation, but now raises a TargetNameParseError [#422]
 
+
 0.5.0 (2024-08-28)
 ==================
 

docs/sbpy/activity/gas.rst

@@ -26,21 +26,42 @@ Two functions provide reference data for the photolysis of gas molecules in opti
 Use :func:`~sbpy.activity.gas.photo_timescale` to retrieve photolysis timescales:
 
   >>> gas.photo_timescale(None)
-  Traceback (most recent call last):
-      ...
-  ValueError: Invalid species None.  Choose from:
-  CH3OH [C94]
-  CN [H92]
-  CO [CE83]
-  CO2 [CE83]
-  H2CO [C94]
-  H2O [CS93]
-  HCN [C94]
-  OH [CS93]
+  species source default         tau               bibcode      
+  ------- ------ ------- ------------------- -------------------
+      H2O   CS93    True           52000.0 s 1993Icar..105..235C
+      H2O   HM15   False   [82940. 49390.] s 2015P&SS..106...11H
+       OH   CS93    True          160000.0 s 1993Icar..105..235C
+      HCN    C94    True           67000.0 s 1994JGR....99.3777C
+    CH3OH    C94    True           77000.0 s 1994JGR....99.3777C
+     H2CO    C94    True            5000.0 s 1994JGR....99.3777C
+       CO   CE83    True         1500000.0 s 1983A&A...126..170C
+      CO2   CE83    True          500000.0 s 1983A&A...126..170C
+      CO2   HM15   False [494800. 210100.] s 2015P&SS..106...11H
+       CN    H92    True [315000. 135000.] s 1992Ap&SS.195....1H
+  -------------------------------------------------------------------------
+  ValueError                                Traceback (most recent call last)
+  Cell In[1], line 1
+  ----> 1 photo_timescale()
+  
+  File /disks/data0/astro/Projects/sbpy/sbpy/activity/gas/core.py:177, in photo_timescale(species, source)
+      167             rows.append(
+      168                 {
+      169                     "species": sp,
+     (...)
+      174                 }
+      175             )
+      176     Table(rows).pprint_all()
+  --> 177     raise ValueError("Invalid species {}".format(species))
+      179 gas = data[species]
+      180 source = default_sources[species] if source is None else source
+  
+  ValueError: Invalid species None
+  
   >>> gas.photo_timescale('H2O')  # doctest: +FLOAT_CMP
   <Quantity 52000. s>
 
-Some sources provide values for the quiet and active Sun (Huebner et al. 1992):
+Some sources provide values for the quiet and active Sun (e.g., Huebner et al.
+1992):
 
   >>> gas.photo_timescale('CN', source='H92')  # doctest: +FLOAT_CMP
   <Quantity [315000., 135000.] s>

sbpy/activity/gas/core.py

@@ -1,7 +1,5 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
-"""activity.gas core
-
-"""
+"""activity.gas core"""
 
 __all__ = [
     "photo_lengthscale",
@@ -18,6 +16,7 @@
 
 import numpy as np
 import astropy.units as u
+from astropy.table import Table
 
 try:
     import scipy
@@ -86,8 +85,7 @@ def photo_lengthscale(species, source=None):
         for k, v in sorted(data.items()):
             summary += "\n{} [{}]".format(k, ", ".join(v.keys()))
 
-        raise ValueError(
-            "Invalid species {}.  Choose from:{}".format(species, summary))
+        raise ValueError("Invalid species {}.  Choose from:{}".format(species, summary))
 
     gas = data[species]
     source = default_sources[species] if source is None else source
@@ -105,24 +103,25 @@ def photo_lengthscale(species, source=None):
     return gamma
 
 
-def photo_timescale(species, source=None):
+def photo_timescale(species=None, source=None):
     """Photodissociation timescale for a gas species.
 
 
     Parameters
     ----------
-    species : string
-        Species to look up.
+    species : string, optional
+        Species to look up.  Use ``None`` to show a list of all species and
+        sources.
 
     source : string, optional
-        Retrieve values from this source.  See references for keys.
+        Retrieve values from this source.  See References for keys.
 
 
     Returns
     -------
     tau : `~astropy.units.Quantity`
-      The timescale at 1 au.  May be a two-element array: (quiet Sun,
-      active Sun).
+      The timescale at 1 au.  May be a two-element array: (quiet Sun, active
+      Sun).
 
 
     Examples
@@ -133,15 +132,17 @@ def photo_timescale(species, source=None):
 
     References
     ----------
-    [CS93] Table IV of Cochran & Schleicher 1993, Icarus 105, 235-253.
-    Quoted for intermediate solar activity.
+    [CS93] Table IV of Cochran & Schleicher 1993, Icarus 105, 235-253. Quoted
+    for intermediate solar activity.
 
     [C94] Crovisier 1994, JGR 99, 3777-3781.
 
     [CE83] Crovisier & Encrenaz 1983, A&A 126, 170-182.
 
     [H92] Huebner et al. 1992, Astroph. & Space Sci. 195, 1-294.
 
+    [HM15] Huebner and Mukherjee 2015, Plan. & Space Sci. 106, 11-45.
+
     """
 
     from .data import photo_timescale as data
@@ -158,12 +159,20 @@ def photo_timescale(species, source=None):
     }
 
     if species not in data:
-        summary = ""
-        for k, v in sorted(data.items()):
-            summary += "\n{} [{}]".format(k, ", ".join(v.keys()))
-
-        raise ValueError(
-            "Invalid species {}.  Choose from:{}".format(species, summary))
+        rows = []
+        for sp, sources in data.items():
+            for source, (tau, bibcode) in sources.items():
+                rows.append(
+                    {
+                        "species": sp,
+                        "source": source,
+                        "default": source == default_sources[sp],
+                        "tau": tau,
+                        "bibcode": ",".join(bibcode.values()),
+                    }
+                )
+        Table(rows).pprint_all()
+        raise ValueError("Invalid species {}".format(species))
 
     gas = data[species]
     source = default_sources[species] if source is None else source
@@ -549,9 +558,7 @@ def f(rho, sigma):
 
                 """
                 return (
-                    rho
-                    * np.exp(-(rho**2) / sigma**2 / 2)
-                    * self._column_density(rho)
+                    rho * np.exp(-(rho**2) / sigma**2 / 2) * self._column_density(rho)
                 )
 
             sigma = aper.sigma.to_value("m")
@@ -662,12 +669,7 @@ def _total_number(self, aper):
             daughter = self.daughter.to(rho.unit)
             y = (rho / daughter).to_value("")
             N *= (daughter / (parent - daughter)).to_value("") * (
-                self._iK0(y)
-                - self._iK0(x)
-                + x**-1
-                - y**-1
-                + self._K1(y)
-                - self._K1(x)
+                self._iK0(y) - self._iK0(x) + x**-1 - y**-1 + self._K1(y) - self._K1(x)
             )
 
         return N
@@ -692,6 +694,7 @@ class VMParent:
     sigma : float
         See VectorialModel documentation.
     """
+
     v_outflow: float
     tau_d: float
     tau_T: float
@@ -711,6 +714,7 @@ class VMFragment:
     tau_T : float
         See VectorialModel documentation.
     """
+
     v_photo: float
     tau_T: float
 
@@ -733,6 +737,7 @@ class VMGridParams:
     radial_substeps : int
         See VectorialModel documentation.
     """
+
     radial_points: int
     angular_points: int
     radial_substeps: int
@@ -755,6 +760,7 @@ class VMParams:
     max_fragment_lifetimes : float
         See VectorialModel documentation.
     """
+
     parent_destruction_level: float
     fragment_destruction_level: float
     max_fragment_lifetimes: float
@@ -780,6 +786,7 @@ class VMFragmentSputterPolar:
         List of fragment densities at the corresponding ``rs[i]`` and
         ``thetas[j]``.
     """
+
     rs: np.ndarray
     thetas: np.ndarray
     fragment_density: np.ndarray
@@ -861,6 +868,7 @@ class VMResult:
         will also measure how long the effects of a single outburst can affect
         the density.
     """
+
     volume_density_grid: np.ndarray = None
     volume_density: np.ndarray = None
     column_density_grid: np.ndarray = None
@@ -1275,13 +1283,11 @@ def _setup_calculations(self) -> None:
         # occupy)
         # NOTE: Equation (16) of Festou 1981 where alpha is the percent
         # destruction of molecules
-        parent_beta_r = - \
-            np.log(1.0 - self.model_params.parent_destruction_level)
+        parent_beta_r = -np.log(1.0 - self.model_params.parent_destruction_level)
         parent_r = parent_beta_r * vp * self.parent.tau_T
         self.vmr.coma_radius = parent_r * u.m
 
-        fragment_beta_r = - \
-            np.log(1.0 - self.model_params.fragment_destruction_level)
+        fragment_beta_r = -np.log(1.0 - self.model_params.fragment_destruction_level)
         # Calculate the time needed to hit a steady, permanent production of
         # fragments
         perm_flow_radius = self.vmr.coma_radius.value + (
@@ -1510,8 +1516,7 @@ def _fragment_sputter(self, r: np.float64, theta: np.float64) -> np.float64:
 
                 # differential addition to the density integrating along dr,
                 # similar to eq. (36) Festou 1981
-                n_r = (p_extinction * f_extinction *
-                       q_r_eps) / (sep_dist**2 * v)
+                n_r = (p_extinction * f_extinction * q_r_eps) / (sep_dist**2 * v)
 
                 sputter += n_r * dr
 
@@ -1567,8 +1572,7 @@ def _compute_fragment_density(self) -> None:
         # Equivalent to summing over j for sin(theta[j]) *
         # fragment_sputter[i][j] with numpy magic
         self.solid_angle_sputter = np.sin(thetas) * self.fragment_sputter
-        self.fast_voldens = 2.0 * np.pi * \
-            np.sum(self.solid_angle_sputter, axis=1)
+        self.fast_voldens = 2.0 * np.pi * np.sum(self.solid_angle_sputter, axis=1)
 
         # Tag with proper units
         self.vmr.volume_density = self.fast_voldens / (u.m**3)
@@ -1633,8 +1637,7 @@ def _column_density_at_rho(self, rho: np.float64) -> np.float64:
         def column_density_integrand(z):
             return self._volume_density(np.sqrt(z**2 + rhosq))
 
-        c_dens = 2 * quad(column_density_integrand, 0,
-                          z_max, epsrel=0.0001)[0]
+        c_dens = 2 * quad(column_density_integrand, 0, z_max, epsrel=0.0001)[0]
 
         # result is in 1/m^2
         return c_dens
@@ -1709,8 +1712,7 @@ def _calc_num_fragments_theory(self) -> np.float64:
         for i, t in enumerate(time_slices[:-1]):
             extinction_one = t / p_tau_T
             extinction_two = time_slices[i + 1] / p_tau_T
-            mult_factor = -np.e ** (-extinction_one) + \
-                np.e ** (-extinction_two)
+            mult_factor = -np.e ** (-extinction_one) + np.e ** (-extinction_two)
             theory_total += self.production_at_time(t) * mult_factor
 
         return theory_total * alpha - edge_adjust