<jupyter_start><jupyter_text>sbpy.activity: gas[sbpy.activity](https://sbpy.readthedocs.io/en/latest/sbpy/activity.html) has classes and functions for models of cometary gas comae.  Comet brightness can be estimated for observations of scattered light or themral emission.  Haser coma modelThe [Haser (1957)](https://ui.adsabs.harvard.edu/abs/1957BSRSL..43..740H/abstract) model presents the distribution of coma gas models destroyed and/or produced in the coma. Column density versus radial distanceTo model the distribution of CO$_2$ in a coma at 1 au from the Sun, use `Haser.column_density`:<jupyter_code>import numpy as np
import astropy.units as u
from sbpy.activity import Haser

v = 0.8 * u.km / u.s  # gas radial expansion speed
tau = 2.0e-6 / u.s * u.au**2 # CO2 photodissociation rate (Crovisier & Encrenaz 1983)

Q = 1e28 / u.s  # CO2 production rate
rh = 1.0 * u.au  # heliocentric distance
gamma = v / tau * rh**2  # coma lengthscale

co2 = Haser(Q, v, gamma)
print('coma lengthscale = {:.0g}\n'.format(gamma))

for rho in 10**np.arange(0, 5) * u.km:  # projected radial distance from the nucleus
    print('{:10.1f} {:10.2e}'.format(rho, co2.column_density(rho)))<jupyter_output>coma lengthscale = 4e+05 km

       1.0 km   3.12e+21 1 / m2
      10.0 km   3.12e+20 1 / m2
     100.0 km   3.12e+19 1 / m2
    1000.0 km   3.09e+18 1 / m2
   10000.0 km   2.89e+17 1 / m2<jupyter_text>Plot the distribution of CO$_2$ as a function of heliocentric distance.<jupyter_code>import matplotlib.pyplot as plt
%matplotlib notebook

rho = np.logspace(4, 7) * u.km

fig = plt.figure(1)
fig.clear()
ax = fig.gca()

for rh in np.arange(1, 11, 2) * u.au:
    co2 = Haser(Q, v, v / tau * rh**2)
    ax.plot(rho, co2.column_density(rho), label='{:.1f}'.format(rh))
    
plt.setp(ax, xlabel='Radial distance (km)', xscale='log',
         ylabel='Column density (m$^{-2}$)', ylim=(1e12, 1e18), yscale='log')
ax.legend()
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Fragment species and production rateTo model a fragment species produced in the coma, provide the lengthscales of the parent and daughter to `Haser`.  `sbpy` will have a library of effective lengthscales based on the literature.  Water and OH from [Cochran & Schleicher (1993)](https://ui.adsabs.harvard.edu/?abs/1993Icar..105..235C) are currently included.<jupyter_code>from sbpy.activity import photo_lengthscale
from sbpy import bib

with bib.Tracking():
    gamma_H2O = photo_lengthscale('H2O')
    gamma_OH = photo_lengthscale('OH')
print(bib.to_text())
print(gamma_H2O, gamma_OH)<jupyter_output>activity.gas.photo_lengthscale:
  H2O photodissociation lengthscale: 1993Icar..105..235C
  OH photodissociation lengthscale: 1993Icar..105..235C

24000.0 km 160000.0 km<jupyter_text>Reproduce the OH production rates for comet 46P/Wirtanen from [Farnham & Schleicher 1998](https://ui.adsabs.harvard.edu/abs/1998A&A...335L..50F/abstract) with `Haser.total_number`:<jupyter_code># OH production rates from Farnham & Schleicher 1998
f98 = 10**np.array((27.89, 27.87, 27.84)) / u.s

eph = {
    'rh': [1.12, 1.14, 1.07] * u.au,
    'delta': [1.41, 1.60, 1.54] * u.au
}
rho = 10**np.array((4.29, 4.64, 4.49)) * u.km  # aperture radius
fluxd = 10**np.array((-10.49, -10.34, -10.47)) * u.erg / u.cm**2 / u.s  # observed flux density

# luminosity per molecule at 1 au for the 0-0 band and comet heliocentric velocity (F&S table 2)
L_N = 10**np.array((-14.576, -14.765, -14.773)) * u.erg / u.s
L_N /= (eph['rh'] / u.au)**2  # scale to observed rh

Q0 = 1 / u.s
v = 1.0 * u.km / u.s

Q = np.zeros_like(f98)
for i in range(len(Q)):
    OH = Haser(Q0, v, gamma_H2O * (eph['rh'][i] / u.au)**2, gamma_OH * (eph['rh'][i] / u.au)**2)
    N = OH.total_number(rho[i])

    Q[i] = Q0 * fluxd[i] / (L_N[i] * N / 4 / np.pi / eph['delta'][i]**2)
    Q[i] = Q[i].to(f98.unit)

print('Mean percent difference: {:.1%}'.format(np.mean(np.abs(Q - f98) / f98)))<jupyter_output>Mean percent difference: 1.1%
{'fraction_non_alphanumeric': 0.15347661188369152, 'fraction_numerical': 0.07509481668773704, 'mean_word_length': 3.8009708737864076, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 14, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '474869', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1761, 'n_tokens_neox': 1539, 'n_words': 412, 'path': 'notebooks/activity/gas.ipynb', 'repo_name': 'jianyangli/sbpy-tutorial', 'size': '3955'}