Gas Comae (sbpy.activity.gas)¶
sbpy.activity.gas provides models and reference data for cometary gases and gas comae.
Photolysis¶
Two functions provide reference data for the photolysis of gas molecules in optically thin comae: photo_lengthscale() and photo_timescale(). The source data is stored in sbpy.activity.gas.data.
photo_lengthscale() provides empirical comae lengthscales (defaults to Cochran and Schleicher 1993)):
>>> from sbpy.activity import gas
>>> gas.photo_lengthscale(None)
Traceback (most recent call last):
...
ValueError: Invalid species None. Choose from:
H2O [CS93]
OH [CS93]
>>> gas.photo_lengthscale('H2O')
<Quantity 24000. km>
Use 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]
>>> gas.photo_timescale('H2O')
<Quantity 52000. s>
Some sources provide values for the quiet and active Sun (Huebner et al. 1992):
>>> gas.photo_timescale('CN', source='H92')
<Quantity [315000., 135000.] s>
With the Bibliography Tracking Module (sbpy.bib), the citation may be discovered:
>>> from sbpy import bib
>>> bib.reset() # clear any old citations
>>> with bib.Tracking():
... tau = gas.photo_timescale('H2O')
>>> print(bib.to_text())
sbpy:
software: sbpy:
Mommert, Kelley, de Val-Borro, Li et al. 2019, The Journal of Open Source Software, Vol 4, 38, 1426
sbpy.activity.gas.core.photo_timescale:
H2O photodissociation timescale:
Cochran & Schleicher 1993, Icarus, Vol 105, 1, 235
Fluorescence¶
Reference data for fluorescence band emission is available via fluorescence_band_strength(). Compute the fluorescence band strength (luminosity per molecule) of the OH 0-0 band at 1 au from the Sun, moving towards the Sun at 1 km/s (defaults to Schleicher and A’Hearn 1988):
>>> import astropy.units as u
>>> LN = gas.fluorescence_band_strength('OH 0-0', -1 * u.km / u.s)
>>> print(LN)
[1.54e-15] erg / s
Gas coma models¶
Haser Model¶
The Haser (1957) 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 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
>>> parent = gas.photo_lengthscale('H2O')
>>> daughter = gas.photo_lengthscale('OH')
>>> coma = gas.Haser(Q, v, parent, daughter)
>>> print(coma.column_density(10 * u.km))
7.099280153851781e+17 1 / m2
>>> print(coma.total_number(1000 * u.km))
1.161357452192558e+30
The gas coma models work with sbpy’s apertures:
>>> from sbpy.activity import AnnularAperture
>>> ap = AnnularAperture((5000, 10000) * u.km)
>>> print(coma.total_number(ap))
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) 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
Vectorial class we may compute the column density and total
number of molecules in an aperture. Parent and daughter data is provided via
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))
2.951278139718558e+17 1 / m2
>>> print(coma.total_number(1000 * u.km))
6.96687966256294e+29
Production Rate calculations¶
Various functions that aid in the calculation of production rates are offered.
Phys has a function called from_jplspec
which takes care of querying the JPL Molecular Spectral Catalog through the use of
jplspec and calculates all the necessary constants needed for
production rate calculations in this module. Yet, the option for the user to
provide their own molecular data is possible through the use of an Phys
object, as long as it has the required information. It is imperative to read
the documentation of the functions in this section to understand what is needed
for each. If the user does not have the necessary data, they can build an object
using JPLSpec:
>>> from sbpy.data.phys import Phys
>>> import astropy.units as u
>>> temp_estimate = 47. * u.K
>>> transition_freq = (230.53799 * u.GHz).to('MHz')
>>> integrated_flux = 0.26 * u.K * u.km / u.s
>>> mol_tag = '^CO$'
>>> mol_data = Phys.from_jplspec(temp_estimate, transition_freq, mol_tag)
>>> mol_data
<QTable length=1>
t_freq temp lgint300 ... degfreedom mol_tag
MHz K MHz nm2 ...
float64 float64 float64 ... int64 int64
-------- ------- --------------------- ... ---------- -------
230538.0 47.0 7.591017628812526e-05 ... 2 28001
Having this information, we can move forward towards the calculation of production rate. The functions that sbpy currently provides to calculate production rates are listed below.
Integrated Line Intensity Conversion¶
The JPL Molecular Spectroscopy Catalog offers the integrated line intensity
at 300 K for a molecule. Yet, in order to calculate production rate, we need
to know the integrated line intensity at a given temperature. This function
takes care of converting the integrated line intensity at 300 K to its equivalent
in the desired temperature using equations provided by the JPLSpec documentation.
For more information on the needed parameters for this function see
intensity_conversion.
>>> from sbpy.activity import intensity_conversion
>>> intl = intensity_conversion(mol_data)
>>> mol_data.apply([intl.value] * intl.unit, name='intl')
11
>>> intl
<Quantity 0.00280051 MHz nm2>
Einstein Coefficient Calculation¶
Einstein coefficients give us insight into the molecule’s probability of spontaneous
absorption, which is useful for production rate calculations.
Unlike catalogs like LAMDA, JPLSpec does not offer the Eistein coefficient
and it must be calculated using equations provided by the JPL Molecular
Spectroscopy Catalog. These equations have been compared to established LAMDA
values of the Einstein Coefficient for HCN and CO, and no more than
a 24% difference has been found between the calculation from JPLSpec and the
LAMDA catalog value. Since JPLSpec and LAMDA are two very different
catalogs with different data, the difference is expected, and the user is
allowed to provide their own Einstein Coefficient if they want. If the user
does want to provide their own Einstein Coefficient, they may do so simply
by appending their value with the unit 1/s to the Phys object, called
mol_data in these examples. For more information on the needed parameters
for this function see einstein_coeff.
>>> from sbpy.activity import einstein_coeff
>>> au = einstein_coeff(mol_data)
>>> mol_data.apply([au.value] * au.unit, name = 'Einstein Coefficient')
12
>>> au
<Quantity 7.03946054e-07 1 / s>
Beta Factor Calculation¶
Returns beta factor based on timescales from gas and distance
from the Sun using an Ephem object. The calculation is
parent photodissociation timescale * (distance from comet to Sun)**2
and it accounts for certain photodissociation and geometric factors needed
in the calculation of total number of molecules total_number
If you wish to provide your own beta factor, you can calculate the equation
expressed in units of AU**2 * s , all that is needed is the timescale
of the molecule and the distance of the comet from the Sun. Once you
have the beta factor you can append it to your mol_data phys object
with the name ‘beta’ or any of its alternative names. For more information on
the needed parameters for this function see beta_factor.
>>> from astropy.time import Time
>>> from sbpy.data import Ephem
>>> from sbpy.activity import beta_factor
>>> target = 'C/2016 R2'
>>> time = Time('2017-12-22 05:24:20', format = 'iso')
>>> ephemobj = Ephem.from_horizons(target, epochs=time.jd)
>>> beta = beta_factor(mol_data, ephemobj)
>>> mol_data.apply([beta.value] * beta.unit, name='beta')
13
>>> beta
<Quantity [13333365.25745597] AU2 s>
Simplified Model for Production Rate¶
activity provides several models to calculate production rates in comets.
One of the models followed by this module is based on the following paper:
The following example shows the usage of the function. This LTE model does not
include photodissociation, but it does serve as way to obtain educated
first guesses for other models within sbpy. For more information on the
parameters that are needed for the function see
from_Drahus.
>>> from sbpy.activity import LTE
>>> vgas = 0.5 * u.km / u.s # gas velocity
>>> lte = LTE()
>>> q = lte.from_Drahus(integrated_flux, mol_data, ephemobj, vgas, aper, b=b)
>>> q
<Quantity 3.59397119e+28 1 / s>
LTE Column Density Calculation¶
To calculate a column density with no previous column density information, we can use equation 10 from Bockelee-Morvan et al. 2004. This function is very useful to obtain a column density with no previous guess for it, and also useful to provide a first guess for the more involved Non-LTE model for column density explained in the next section.
>>> cdensity = lte.cdensity_Bockelee(integrated_flux, mol_data)
>>> mol_data.apply([cdensity.value] * cdensity.unit, name='cdensity')
Non-LTE Column Density Calculation¶
Once the user has a guess for their column density, the user can then
implement the sbpy.activity NonLTE function sbpy.activity.NonLTE.from_pyradex.
This function calculates the best fitting column density for the integrated
flux data using the python wrapper pyradex of the Non-LTE iterative code RADEX.
The code utilizes the LAMDA catalog collection of molecular data files,
presently this is the only functionality available, yet in the future a function
will be provided by sbpy to build your own molecular data file from JPLSpec
for use in this function. The code will look for a ‘cdensity’ column value
within mol_data to use as its first guess. For a more detailed look at the
input parameters, please see from_pyradex.
>>> from sbpy.activity import NonLTE
>>> nonlte = NonLTE()
>>> cdensity = nonlte.from_pyradex(integrated_flux, mol_data, iter=500)
>>> mol_data.apply([cdensity.value] * cdensity.unit, name='cdensity')
Note that for this calculation the installation of pyradex is needed. Pyradex
is a python wrapper for the RADEX fortran code. See pyradex installation and
README file
for installation instruction and tips as well as a briefing of how pyradex
works and what common errors might arise. You need to make sure you have a
fortran compiler installed in order for pyradex to work (gfortran works and can
be installed with homebrew for easier management).
Total Number¶
In order to obtain our total number of molecules from flux data, we use the millimeter/submillimeter spectroscopy beam factors explained and detailed in equation 1.3 from:
Drahus, M. (2010). Microwave observations and modeling of the molecularcoma in comets. PhD Thesis, Georg-August-Universität Göttingen.
If the user prefers to give the total number, they may do so by appending
to the mol_data Phys object with the name total_number or
any of its alternative names. For more information on the needed parameters
for this function see total_number.
>>> from sbpy.activity import total_number
>>> integrated_flux = 0.26 * u.K * u.km / u.s
>>> b = 0.74
>>> aper = 10 * u.m
>>> tnum = total_number(integrated_flux, mol_data, aper, b)
>>> mol_data.apply([tnum], name='total_number')
14
>>> tnum
<Quantity [2.93988826e+26]>
Haser Model for Production Rate¶
Another model included in the module is based off of the model in the following literature:
This model is well-known as the Haser model. In the case of our implementation
the function takes in an initial guess for the production rate, and uses the
module found in gas to find a ratio between the model
total number of molecules and the number of molecules calculated from the data
to scale the model Q and output the new production rate from the result. This
LTE model does account for the effects of photolysis. For more information
on the parameters that are needed for the function see
from_Haser().
>>> from sbpy.activity import Haser, photo_timescale, from_Haser
>>> Q_estimate = 3.5939*10**(28) / u.s
>>> parent = photo_timescale('CO') * vgas
>>> coma = Haser(Q_estimate, vgas, parent)
>>> Q = from_Haser(coma, mol_data, aper=aper)
>>> Q
<Quantity [[9.35795579e+27]] 1 / s>
For more involved examples and literature comparison for any of the production rate modules, please see notebook examples.
Reference/API¶
activity.gas core
Functions¶
|
Photodissociation lengthscale for a gas species. |
|
Photodissociation timescale for a gas species. |
|
Fluorescence band strength. |
Classes¶
|
Haser coma model. |
|
Vectorial model for fragments in a coma produced |