<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%
StarcoderdataJupyterStructuredCleanDedup