# Dust comae (sbpy.activity.dust)¶

## Afρ and εfρ¶

sbpy has two classes to support observations and models of a coma continuum: Afrho and Efrho.

The Afρ parameter of A’Hearn et al (1984) is based on observations of idealized cometary dust comae. It is proportional to the observed flux density within a circular aperture. The quantity is the product of dust albedo, dust filling factor, and the radius of the aperture at the distance of the comet. It carries the units of ρ (length), and under certain assumptions is proportional to the dust production rate of the comet:

$Afρ = \frac{4 Δ^2 r_h^2}{ρ}\frac{F_c}{F_⊙}$

where Δ and ρ have the same (linear) units, but $$r_h$$ is in units of au. $$F_c$$ * is the flux density of the comet in the aperture, and $$F_⊙$$ is that of the Sun at 1 au in the same units. See A’Hearn et al. (1984) and Fink & Rubin (2012) for more details.

The εfρ parameter is the thermal emission counterpart to Afρ, replacing albedo with IR emissivity, ε, and the solar spectrum with the Planck function, B:

$εfρ = \frac{F_c Δ^2}{π ρ B(T_c)}$

where $$T_c$$ is the spectral temperature of the continuum (Kelley et al. 2013).

### Afρ and εfρ are quantities¶

Afrho and Efrho are subclasses of astropy’s Quantity and carry units of length.

>>> import numpy as np
>>> import astropy.units as u
>>> from sbpy.activity import Afrho, Efrho
>>>
>>> afrho = Afrho(100 * u.cm)
>>> print(afrho)    # doctest: +FLOAT_CMP
100.0 cm
>>> efrho = Efrho(afrho * 3.5)
>>> print(efrho)    # doctest: +FLOAT_CMP
350.0 cm


They may be converted to other units of length just like any Quantity:

>>> afrho.to('m')    # doctest: +FLOAT_CMP
<Afrho 1. m>


### Convert to/from flux density¶

The quantities may be initialized from flux densities. Here, we reproduce one of the calculations from the original A’Hearn et al. (1984) work:

• $$r_h$$ = 4.785 au

• Δ = 3.822 au

• phase angle = 3.3°

• aperture radius = 9.8” (27200 km)

• $$\log{F_λ}$$ = -13.99 erg/(cm2 s Å) at λ = 5240 Å

The solar flux density at 1 au is also needed. We use 1868 W/(m2 μm).

>>> from sbpy.data import Ephem
>>> from sbpy.calib import solar_fluxd
>>>
>>> solar_fluxd.set({
...     'λ5240': 1868 * u.W / u.m**2 / u.um,
...     'λ5240(lambda pivot)': 5240 * u.AA
... })              # doctest: +IGNORE_OUTPUT
>>>
>>> flam = 10**-13.99 * u.Unit('erg/(s cm2 AA)')
>>> aper = 27200 * u.km
>>>
>>> eph = Ephem.from_dict({'rh': 4.785 * u.au, 'delta': 3.822 * u.au})
>>>
>>> afrho = Afrho.from_fluxd('λ5240', flam, aper, eph)
>>> print(afrho)    # doctest: +FLOAT_CMP
6029.90248952895 cm


Which is within a few percent of 6160 cm computed by A’Hearn et al.. The difference is likely due to the assumed solar flux density in the bandpass.

The Afrho class may be converted to a flux density, and the original value is recovered.

>>> f = afrho.to_fluxd('λ5240', aper, eph).to('erg/(s cm2 AA)')
>>> print(np.log10(f.value))    # doctest: +FLOAT_CMP
-13.99


Afrho works seamlessly with sbpy’s spectral calibration framework (Spectral Standards and Photometric Calibration (sbpy.calib)) when the astropy affiliated package synphot is installed. The solar flux density (via solar_fluxd) is not required, but instead the spectral wavelengths or the system transmission of the instrument and filter:

>>> wave = [0.4, 0.5, 0.6] * u.um
>>> print(afrho.to_fluxd(wave, aper, eph))
[7.76770018e-14 1.05542873e-13 9.57978939e-14] W / (m2 um)

>>> from synphot import SpectralElement, Box1D
>>> # bandpass width is a guess
>>> bp = SpectralElement(Box1D, x_0=5240 * u.AA, width=50 * u.AA)
>>> print(Afrho.from_fluxd(bp, flam, aper, eph))
5994.110239075767 cm


### Thermal emission with εfρ¶

The Efrho class has the same functionality as the Afrho class. The most important difference is that εfρ is calculated using a Planck function and temperature. sbpy follows common practice and parameterizes the temperature as a constant scale factor of $$T_{BB} = 278\,r_h^{1/2}$$ K, the equilibrium temperature of a large blackbody sphere at a distance $$r_h$$ from the Sun.

Reproduce the εfρ of 246P/NEAT from Kelley et al. (2013).

>>> wave = [15.8, 22.3] * u.um
>>> fluxd = [25.75, 59.2] * u.mJy
>>> aper = 11.1 * u.arcsec
>>> eph = Ephem.from_dict({'rh': 4.28 * u.au, 'delta': 3.71 * u.au})
>>> efrho = Efrho.from_fluxd(wave, fluxd, aper, eph)
>>> for i in range(len(wave)):
...     print('{:5.1f} at {:.1f}'.format(efrho[i], wave[i]))
406.2 cm at 15.8 um
427.9 cm at 22.3 um


Compare to 397.0 cm and 424.6 cm listed in Kelley et al. (2013).

### To/from magnitudes¶

Afrho and Efrho work with astropy’s magnitude units. If the conversion between Vega-based magnitudes is required, then sbpy’s calibration framework (Spectral Standards and Photometric Calibration (sbpy.calib)) will be used.

Estimate the Afρ of comet C/2012 S1 (ISON) based on Pan-STARRS 1 photometry in the r band (Meech et al. 2013)

>>> w = 0.617 * u.um
>>> m = 16.02 * u.ABmag
>>> aper = 5 * u.arcsec
>>> eph = {'rh': 5.234 * u.au, 'delta': 4.275 * u.au, 'phase': 2.6 * u.deg}
>>> afrho = Afrho.from_fluxd(w, m, aper, eph)
>>> print(afrho)
1948.496075629444 cm
>>> m2 = afrho.to_fluxd(w, aper, eph, unit=u.ABmag)
>>> print(m2)
16.02 mag(AB)


### Phase angles and functions¶

Phase angle was not used in the previous section. In the Afρ formalism, “albedo” includes the scattering phase function, and is more precisely written A(θ), where θ is the phase angle. The default behavior for Afrho is to compute A(θ)fρ as opposed to A(0°)fρ. Returning to the A’Hearn et al. data, we scale Afρ to 0° from 3.3° phase using the to_phase() method:

>>> afrho = Afrho(6029.9 * u.cm)
>>> print(afrho.to_phase(0 * u.deg, 3.3 * u.deg))    # doctest: +FLOAT_CMP
6886.825981017757 cm


The default phase function is the Halley-Marcus composite phase function (phase_HalleyMarcus()). Any function or callable object that accepts an angle as a Quantity and returns a scalar value may be used:

>>> Phi = lambda phase: 10**(-0.016 / u.deg * phase.to('deg'))
>>> print(afrho.to_phase(0 * u.deg, 3.3 * u.deg, Phi=Phi))    # doctest: +FLOAT_CMP
6809.419810008357 cm


To correct an observed flux density for the phase function, use the phasecor option of to_fluxd() and from_fluxd() methods:

>>> flam = 10**-13.99 * u.Unit('erg/(s cm2 AA)')
>>> aper = 27200 * u.km
>>> eph = Ephem.from_dict({
...     'rh': 4.785 * u.au,
...     'delta': 3.822 * u.au,
...     'phase': 3.3 * u.deg
... })
>>>
>>> afrho = Afrho.from_fluxd('λ5240', flam, aper, eph, phasecor=True)
>>> print(afrho)    # doctest: +FLOAT_CMP
6886.828824340642 cm


### Apertures¶

Other apertures may be used, as long as they can be converted into an equivalent radius, assuming a coma with a 1/ρ surface brightness distribution. activity has a collection of useful geometries.

>>> from sbpy.activity import CircularAperture, AnnularAperture, RectangularAperture, GaussianAperture
>>> apertures = (
...   ( '10" radius circle', CircularAperture(10 * u.arcsec)),
...   (    '5"–10" annulus', AnnularAperture([5, 10] * u.arcsec)),
...   (       '2"x10" slit', RectangularAperture([2, 10] * u.arcsec)),
...   ('σ=5" Gaussian beam', GaussianAperture(5 * u.arcsec))
... )
>>> for name, aper in apertures:
...     afrho = Afrho.from_fluxd('λ5240', flam, aper, eph)
...     print('{:18s} = {:5.0f}'.format(name, afrho))    # doctest: +FLOAT_CMP
10" radius circle =  5917 cm
5"–10" annulus = 11834 cm
2"x10" slit = 28114 cm
σ=5" Gaussian beam =  9442 cm


## Reference/API¶

### sbpy Activity: Dust¶

All things dust coma related.

### Functions¶

 phase_HalleyMarcus(phase) Halley-Marcus composite dust phase function.

### Classes¶

 Afrho Coma dust quantity for scattered light. Efrho Coma dust quantity for thermal emission.