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)
100.0 cm
>>> efrho = Efrho(afrho * 3.5)
>>> print(efrho)
350.0 cm
They may be converted to other units of length just like any Quantity
:
>>> afrho.to('m')
<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/(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
... })
>>>
>>> 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)
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))
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.76770018e14 1.05542873e13 9.57978939e14] 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 Vegabased 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 PanSTARRS 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))
6886.825981017757 cm
The default phase function is the HalleyMarcus 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))
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)
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))
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¶

HalleyMarcus composite dust phase function. 
Classes¶
Coma dust quantity for scattered light. 

Coma dust quantity for thermal emission. 