Merge pull request #305 from mkelley/vectorial-test-docs

VectorialModel: literature test and documentation

docs/sbpy/activity/gas.rst

@@ -1,6 +1,11 @@
 Gas Comae (`sbpy.activity.gas`)
 ===============================
 
+`sbpy.activity.gas` provides models and reference data for cometary gases and gas comae.
+
+.. toctree::
+   :maxdepth: 2
+
 Photolysis
 ----------
 
@@ -68,7 +73,15 @@ Reference data for fluorescence band emission is available via :func:`~sbpy.acti
 Gas coma models
 ---------------
 
-The Haser (1957) model for parent and daughter species is included, with some calculation enhancements based on Newburn and Johnson (1978).  With `~sbpy.activity.gas.Haser`, we may compute the column density and total number of molecules within aperture:
+Haser Model
+^^^^^^^^^^^
+
+The `Haser (1957)
+<https://ui.adsabs.harvard.edu/abs/1957BSRSL..43..740H/abstract>`_ model is an
+analytical approach to solving for the spatial distribution of parent and
+daughter species.  It is included with some calculation enhancements based on
+Newburn and Johnson (1978).  With `~sbpy.activity.gas.Haser`, we may compute the
+column density and total number of molecules within an aperture:
 
   >>> Q = 1e28 / u.s        # production rate
   >>> v = 0.8 * u.km / u.s  # expansion speed
@@ -87,6 +100,55 @@ The gas coma models work with sbpy's apertures:
   >>> print(coma.total_number(ap))    # doctest: +FLOAT_CMP
   3.8133654170856037e+31
 
+Vectorial Model
+^^^^^^^^^^^^^^^
+
+.. warning::
+
+  Literature tests with the Vectorial model are generally in agreement at the
+  20% level or better.  The cause for the differences with the Festou FORTRAN
+  code are not yet precisely known.  Help testing this feature is appreciated.
+
+The Vectorial model (`Festou 1981
+<https://ui.adsabs.harvard.edu/abs/1981A%26A....95...69F/abstract>`_) describes
+the spatial distribution of coma photolysis products.  Unlike the Haser model,
+daughter products in the Vectorial model may receive an additional velocity
+component from the excess energy after photodissociation.  With the
+`~sbpy.activity.gas.Vectorial` class we may compute the column density and total
+number of molecules in an aperture.  Parent and daughter data is provided via
+`~sbpy.data.Phys` objects, with the following required parameters:
+
++------------------+-----------+------+-------------------------------------------------------+
+|     Species      | Parameter | Unit |        Description (and Festou 1981 variable)         |
++==================+===========+======+=======================================================+
+| parent           | v_outflow | m/s  | outflow velocity (u)                                  |
++------------------+-----------+------+-------------------------------------------------------+
+| parent           | sigma     | m**2 | collisional cross-sectional area                      |
++------------------+-----------+------+-------------------------------------------------------+
+| parent           | tau_d     | s    | photodissociation lifetime                            |
++------------------+-----------+------+-------------------------------------------------------+
+| parent, daughter | tau_T     | s    | total lifetime considering all destruction mechanisms |
++------------------+-----------+------+-------------------------------------------------------+
+| daughter         | v_photo   | m/s  | photodissociation velocity (v_R)                      |
++------------------+-----------+------+-------------------------------------------------------+
+
+  >>> from sbpy.data import Phys
+  >>> water = Phys.from_dict({
+  ...     'tau_T': gas.photo_timescale('H2O') * 0.8,  # approximate
+  ...     'tau_d': gas.photo_timescale('H2O'),
+  ...     'v_outflow': 0.85 * u.km / u.s,
+  ...     'sigma': 3e-16 * u.cm**2
+  ... })
+  >>> hydroxyl = Phys.from_dict({
+  ...     'tau_T': gas.photo_timescale('OH') * 0.8,  # approximate
+  ...     'v_photo': 1.05 * u.km / u.s
+  ... })
+  >>> Q = 1e28 / u.s        # water production rate
+  >>> coma = gas.VectorialModel(Q, water, hydroxyl)
+  >>> print(coma.column_density(10 * u.km))    # doctest: +FLOAT_CMP
+  2.951278139718558e+17 1 / m2
+  >>> print(coma.total_number(1000 * u.km))    # doctest: +FLOAT_CMP
+  6.96687966256294e+29
 
 Production Rate calculations
 ----------------------------

examples/activity/vectorial-model.ipynb

@@ -53,10 +53,29 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 5,
    "id": "narrow-script",
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "<QTable length=1>\n",
+      " tau_T   tau_d  v_outflow  sigma \n",
+      "   s       s      km / s    cm2  \n",
+      "float64 float64  float64  float64\n",
+      "------- ------- --------- -------\n",
+      "41600.0 52000.0      0.85   3e-16\n",
+      "<QTable length=1>\n",
+      " tau_T   v_photo\n",
+      "   s      km / s\n",
+      "float64  float64\n",
+      "-------- -------\n",
+      "128000.0    1.05\n"
+     ]
+    }
+   ],
    "source": [
     "# Rough parameters for a quick example calculation\n",
     "\n",
@@ -79,7 +98,10 @@
     "    'tau_T': sba.photo_timescale('OH') * 0.8,\n",
     "    # Velocity after dissociation, taken from Cochran and Schleicher, 1993\n",
     "    'v_photo': 1.05 * u.km/u.s\n",
-    "    })"
+    "    })\n",
+    "\n",
+    "print(parent_species)\n",
+    "print(fragment_species)"
    ]
   },
   {
@@ -652,7 +674,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -666,7 +688,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.2"
+   "version": "3.9.7"
   }
  },
  "nbformat": 4,