Gas Comae (sbpy.activity.gas
)¶
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') # doctest: +FLOAT_CMP
<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') # doctest: +FLOAT_CMP
<Quantity 52000. s>
Some sources provide values for the quiet and active Sun (Huebner et al. 1992):
>>> gas.photo_timescale('CN', source='H92') # doctest: +FLOAT_CMP
<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()) # doctest: +REMOTE_DATA
sbpy:
software: sbpy:
Mommert, Kelley, De ValBorro, 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 00 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 00', 1 * u.km / u.s)
>>> print(LN) # doctest: +FLOAT_CMP
[1.54e15] erg / s
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 Haser
, we may compute the column density and total number of molecules within 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)) # doctest: +FLOAT_CMP
7.099280153851781e+17 1 / m2
>>> print(coma.total_number(1000 * u.km)) # doctest: +FLOAT_CMP
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)) # doctest: +FLOAT_CMP
3.8133654170856037e+31
Production Rate calculations¶
productionrate
offers various functions that aid in the calculation
of production rates. 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.591017628812526e05 ... 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 follow the link
for intensity_conversion
under Reference/API section.
>>> 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 follow the link for einstein_coeff
under Reference/API section.
>>> 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.03946054e07 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_nocd
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 follow the link for
beta_factor
under Reference/API section.
>>> from astropy.time import Time
>>> from sbpy.data import Ephem
>>> from sbpy.activity import beta_factor
>>> target = 'C/2016 R2'
>>> time = Time('20171222 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 follow
the link for the function from_Drahus
in LTE
under Reference/API section.
>>> 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 BockeleeMorvan 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 NonLTE 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')
NonLTE Column Density Calculation¶
Once the user has a guess for their column density, the user can then
implement the sbpy.activity
NonLTE function from_pyradex
. This function
calculates the best fitting column density for the integrated flux data using
the python wrapper pyradex of the NonLTE 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 and defaults, search for from_pyradex
in
NonLTE
under the Reference/API section.
>>> from sbpy.activty 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 here
and
here
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, GeorgAugustUniversitä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_nocd
or
any of its alternative names. For more information on the needed parameters
for this function follow the link for total_number_nocd
under Reference/API section.
>>> 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 wellknown 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 follow
the link for the function from_Haser
in LTE
under Reference/API section.
>>> 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¶
sbpy Activity: Gas¶
All things gas coma realated.
Functions¶

Returns beta factor based on timescales from 

Einstein coefficient from molecular data 

Fluorescence band strength. 

Calculate production rate for 

Returns conversion of the integrated line intensity at 300 K (from the JPL molecular spectra catalog) to a chosen temperature 

Photodissociation lengthscale for a gas species. 

Photodissociation timescale for a gas species. 

Equation relating number of molecules with column density, the aperture, and geometry given and accounting for photodissociation, derived from data provided. 
Classes¶

Haser coma model. 
LTE Methods for calculating production_rate 

Class method for non LTE production rate models Not Yet implemented 