# 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 auphase angle = 3.3°

aperture radius = 9.8” (27200 km)

\(\log{F_λ}\) = -13.99 erg/(cm

^{2}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
```