Activity Sketches

A collection of comet comae models for gas and dust activity.

Table of contents

• Haser coma model
• Vectorial coma model
• Syndynes and synchrones

Haser coma model

Create a Haser coma model (Haser 1957) for OH at comet Encke and image it from Earth on 2017 Jan 1.

Code

import aplpy
from astropy.time import Time
import astropy.units as u
from sbpy import activity, utils
from sbpy.data import Ephemeris

# Molecular scale lengths at 1 au, Cochran and Schleicher 1993
gamma_h2o = 2.4e4 * u.km
gamma_oh = 1.6e5 * u.km

# Fluorescence band efficiency for OH 0-0 band at 3090 Å (r-dot = 0 km/s) (Schleicher and A'Hearn 1988)
L_N = 2e-15 * u.erg / u.s / u.molecule * u.au**2

# Create the coma
Q_oh = 1e28 * u.molecule / u.s
oh_coma = activity.Haser(Q_oh, gamma_h2o, gamma_oh)

# Get observing geometry for comet Encke on 2017 Jan 1
comet = '2P'
observer = 'G37'  # Lowell/DCT
date = Time('2017-01-01', scale='utc')
eph = Ephemeris.from_horizons(comet, observer, epoch=date)

# Region of interest a 1000 km radius circle: CircularAperture(1000 * u.km),
# or RectangularAperture(100 * u.km, 1000 * u.km) for long-slit spectroscopy,
# or GaussianAperture(1000 * u.km) for radio telescopes.
aperture = activity.CircularAperture(1000 * u.km)

# Show some calculations.
print('Comet Encke observed from the DCT: rh={:.2f}, Delta={:.2f}'.format(eph.rh, eph.delta))
print('Volume density at 1000 km: {:.4g}'.format(oh_coma.volume_density(aperture, eph)))
print('Column density at 1000 km: {:.4g}'.format(oh_coma.column_density(aperture, eph)))
print('Total number within 1000 km: {:.4g}'.format(oh_coma.total_number(aperture, eph)))

F = L_N / eph.rh**2 * oh_coma.total_number(aperture, eph)
print('3090 Å OH band flux within 1000 km radius aperture: {:.4g}'.format(F))

# Simulate an image of the 0-0 band and display it.
fov = (100, 100) * u.pixel
scale = 1 * u.arcsec / u.pixel
rho = utils.rarray(fov) * scale
im = L_N / eph.rh**2 * oh_coma.column_density(rho, eph)
fig = aplpy.FITSFigure(im.value)
fig.show_colorscale(cmap='viridis')
fig.colorbar.set_axis_label_text('Flux ({})'.format(im.unit))

Notes

• astropy.units does not have the molecule unit. We should create one, though. It does have mole.
• There probably should be some close relationship between Ephemeris and astropy's SkyCoord.

Vectorial coma model

For the vectorial model of Festou (1981), we create and instance of the Vectorial class, but otherwise the script is the same. Following the Web Vectorial Model inputs:

v_h2o = 0.85 * u.km / u.s  # speed of the parent
tau_h2o = 86430 * u.s      # total lifetime of the parent
tau_h2o_d = 101730 * u.s   # parent disassociative lifetime
v_oh = 1.05 * u.km / u.s   # speed of daughter
tau_oh = 1.29e5 * u.s      # total lifetime of daughter
oh_coma = activity.Vectorial(Q_oh, v_h2o, tau_h2o, tau_h2o_d, v_oh, tau_oh)

Syndynes and synchrones

Create a set of syndynes for comet Encke as seen from Earth on 2017 Jan 1.

Code

import matplotlib.pyplot as plt
from astropy.time import Time
import astropy.units as u
import astropy.coordinates as coords
from sbpy import activity
from sbpy.data import Ephemeris, OrbElements

# To generate syndynes, we need heliocentric rectangular coordinates of our target at
# every time step.  Here, we will compute these from the orbital elements, but they
# could be computed from any class with a `state` method.  
comet = '2P'
date = Time('2017-01-01', scale='utc')
orbit = OrbElements.from_horizons(comet, epoch=date)

print("""Encke's position and velocity on {date}:
  r = {rv[0]}
  v = {rv[1]}
""".format(date=date, rv=orbit.state(date)))

syn = activity.Syndynes(orbit, date)

# Syndynes are parameterized with the unitless `beta`, the ratio of the force from solar
# radiation pressure to the force from solar gravity.
betas = np.logspace(0, -3, 13)
ages = np.linspace(0, 100, 101)  # time steps for each beta

# Set up figure
fig = plt.figure(1)
fig.clear()
ax = fig.add_subplot(111, projection='polar')

dct = coords.EarthLocation.of_site('DCT')
syn.plot_syndynes(betas, ages, location=dct)
syn.plot_synchrones(betas, ages[::10], location=dct)  # only plot every 10th time step
plt.show()

The plot_* convenience functions call syn.syndynes and syn.synchrones. Each call to syn.syndynes returns an array of SkyCoord, one for each time step. The coordinates are computed as Cartesian coordinates in a heliocentric reference frame. If an observer location is provided, they are transformed into a coordinate system at the observer.

# syndynes
for beta in betas:
    c = syn.syndyne(beta, ages, location=location)
    r = c[0].separation(c)
    pa = c[0].position_angle(c)
    ax.plot(r.arcsec, pa.radian)

# synchrones
for age in ages:
    c = syn.synchrone(betas, age, location=location)
    r = c[0].separation(c)
    pa = c[0].position_angle(c)
    ax.plot(r.arcsec, pa.radian)

Notes

• We should add IAU observatory codes to astropy, or otherwise figure out how to use them in the EarthLocation.of_site() call.