# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
sbpy Activity: Dust
===================
All things dust coma related.
"""
__all__ = ["phase_HalleyMarcus", "Afrho", "Efrho"]
__doctest_requires__ = {"Afrho.to_fluxd": ["astropy>=5.3"]}
from warnings import warn
import abc
import numpy as np
import astropy.units as u
import astropy.version as astropy_version
try:
from astropy.utils.misc import InheritDocstrings
except ImportError:
InheritDocstrings = None
from .. import bib
from ..calib import Sun
from ..spectroscopy import BlackbodySource
from .. import data as sbd
from .. import exceptions as sbe
from .. import units as sbu
from ..spectroscopy.sources import SinglePointSpectrumError
from .core import Aperture
[docs]@bib.cite(
{
"Halley-Marcus phase function": "2011AJ....141..177S",
"Halley phase function": "1998Icar..132..397S",
"Marcus phase function": "2007ICQ....29...39M",
}
)
def phase_HalleyMarcus(phase):
"""Halley-Marcus composite dust phase function.
Uses `~scipy.interpolate` for spline interpolation, otherwise uses
linear interpolation from `~numpy.interp`.
Parameters
----------
phase : `~astropy.units.Quantity`, `~astropy.coordinate.Angle`
Phase angle.
Returns
-------
Phi : float, `~np.ndarray`
Notes
-----
The Halley-Marcus phase function was first used by Schleicher and
Bair (2011), but only described in detail by Schleicher and Marcus
(May 2010) online at:
https://asteroid.lowell.edu/comet/dustphase.html
"To distinguish this curve from others, we designate this as
the HM phase function, for the sources of the two components:
Halley and Marcus, where the Halley curve for smaller phase
angles comes from our previous work (Schleicher et al. 1998)
while Joe Marcus has fit a Henyey-Greenstein function to a
variety of mid- and large-phase angle data sets (Marcus 2007);
see here for details. Note that we do not consider our
composite curve to be a definitive result, but rather
appropriate for performing first-order adjustments to dust
measurements for changing phase angle."
References
----------
Schleicher & Bair 2011, AJ 141, 177.
Schleicher, Millis, & Birch 1998, Icarus 132, 397-417.
Marcus 2007, International Comets Quarterly 29, 39-66.
Examples
--------
>>> from sbpy.activity import phase_HalleyMarcus
>>> import astropy.units as u
>>> phase_HalleyMarcus(0 * u.deg) # doctest: +FLOAT_CMP
1.0
>>> phase_HalleyMarcus(15 * u.deg) # doctest: +FLOAT_CMP
5.8720e-01
"""
try:
import scipy
from scipy.interpolate import splrep, splev
except ImportError:
scipy = None
th = np.arange(181)
ph = np.array(
[
1.0000e00,
9.5960e-01,
9.2170e-01,
8.8590e-01,
8.5220e-01,
8.2050e-01,
7.9060e-01,
7.6240e-01,
7.3580e-01,
7.1070e-01,
6.8710e-01,
6.6470e-01,
6.4360e-01,
6.2370e-01,
6.0490e-01,
5.8720e-01,
5.7040e-01,
5.5460e-01,
5.3960e-01,
5.2550e-01,
5.1220e-01,
4.9960e-01,
4.8770e-01,
4.7650e-01,
4.6590e-01,
4.5590e-01,
4.4650e-01,
4.3770e-01,
4.2930e-01,
4.2150e-01,
4.1420e-01,
4.0730e-01,
4.0090e-01,
3.9490e-01,
3.8930e-01,
3.8400e-01,
3.7920e-01,
3.7470e-01,
3.7060e-01,
3.6680e-01,
3.6340e-01,
3.6030e-01,
3.5750e-01,
3.5400e-01,
3.5090e-01,
3.4820e-01,
3.4580e-01,
3.4380e-01,
3.4210e-01,
3.4070e-01,
3.3970e-01,
3.3890e-01,
3.3850e-01,
3.3830e-01,
3.3850e-01,
3.3890e-01,
3.3960e-01,
3.4050e-01,
3.4180e-01,
3.4320e-01,
3.4500e-01,
3.4700e-01,
3.4930e-01,
3.5180e-01,
3.5460e-01,
3.5760e-01,
3.6090e-01,
3.6450e-01,
3.6830e-01,
3.7240e-01,
3.7680e-01,
3.8150e-01,
3.8650e-01,
3.9170e-01,
3.9730e-01,
4.0320e-01,
4.0940e-01,
4.1590e-01,
4.2280e-01,
4.3000e-01,
4.3760e-01,
4.4560e-01,
4.5400e-01,
4.6270e-01,
4.7200e-01,
4.8160e-01,
4.9180e-01,
5.0240e-01,
5.1360e-01,
5.2530e-01,
5.3750e-01,
5.5040e-01,
5.6380e-01,
5.7800e-01,
5.9280e-01,
6.0840e-01,
6.2470e-01,
6.4190e-01,
6.5990e-01,
6.7880e-01,
6.9870e-01,
7.1960e-01,
7.4160e-01,
7.6480e-01,
7.8920e-01,
8.1490e-01,
8.4200e-01,
8.7060e-01,
9.0080e-01,
9.3270e-01,
9.6640e-01,
1.0021e00,
1.0399e00,
1.0799e00,
1.1223e00,
1.1673e00,
1.2151e00,
1.2659e00,
1.3200e00,
1.3776e00,
1.4389e00,
1.5045e00,
1.5744e00,
1.6493e00,
1.7294e00,
1.8153e00,
1.9075e00,
2.0066e00,
2.1132e00,
2.2281e00,
2.3521e00,
2.4861e00,
2.6312e00,
2.7884e00,
2.9592e00,
3.1450e00,
3.3474e00,
3.5685e00,
3.8104e00,
4.0755e00,
4.3669e00,
4.6877e00,
5.0418e00,
5.4336e00,
5.8682e00,
6.3518e00,
6.8912e00,
7.4948e00,
8.1724e00,
8.9355e00,
9.7981e00,
1.0777e01,
1.1891e01,
1.3166e01,
1.4631e01,
1.6322e01,
1.8283e01,
2.0570e01,
2.3252e01,
2.6418e01,
3.0177e01,
3.4672e01,
4.0086e01,
4.6659e01,
5.4704e01,
6.4637e01,
7.7015e01,
9.2587e01,
1.1237e02,
1.3775e02,
1.7060e02,
2.1348e02,
2.6973e02,
3.4359e02,
4.3989e02,
5.6292e02,
7.1363e02,
8.8448e02,
1.0533e03,
1.1822e03,
1.2312e03,
]
)
_phase = np.abs(u.Quantity(phase, "deg").value)
if scipy:
Phi = splev(_phase, splrep(th, ph))
else:
warn(
sbe.OptionalPackageUnavailable(
"scipy is not present, using linear interpolation."
)
)
Phi = np.interp(_phase, th, ph)
if np.iterable(_phase):
Phi = np.array(Phi).reshape(np.shape(_phase))
else:
Phi = float(Phi)
return Phi
# checks if the installed astropy version is up-to-date
if astropy_version.major < 4:
class DustComaQuantityMeta(InheritDocstrings, abc.ABCMeta):
pass
else:
# astropy >= 4.0
DustComaQuantityMeta = abc.ABCMeta
class DustComaQuantity(u.SpecificTypeQuantity, metaclass=DustComaQuantityMeta):
"""Abstract base class for dust coma photometric models: Afrho, Efrho."""
_equivalent_unit = u.meter
_include_easy_conversion_members = False
def __new__(cls, value, unit=None, dtype=None, copy=True):
return super().__new__(cls, value, unit=unit, dtype=dtype, copy=copy)
@classmethod
def from_fluxd(cls, wfb, fluxd, aper, eph, **kwargs):
"""Initialize from spectral flux density.
Parameters
----------
wfb : `~astropy.units.Quantity`, `~synphot.SpectralElement`, list
Wavelengths, frequencies, bandpass, or list of
bandpasses of the observation. Bandpasses require
`~synphot`.
fluxd : `~astropy.units.Quantity`
Flux density per unit wavelength or frequency.
aper : `~astropy.units.Quantity` or `~sbpy.activity.Aperture`
Aperture of the observation as a circular radius (length
or angular units), or as an `~sbpy.activity.Aperture`.
eph: dictionary-like, `~sbpy.data.Ephem`
Ephemerides of the comet. Required fields: 'rh', 'delta'.
Optional: 'phase'.
**kwargs
Keyword arguments for `~to_fluxd`.
"""
fluxd1cm = cls(1 * u.cm).to_fluxd(wfb, aper,
eph, unit=fluxd.unit, **kwargs)
if isinstance(fluxd1cm, u.Magnitude):
coma = cls((fluxd - fluxd1cm).physical * u.cm)
else:
coma = cls((fluxd / fluxd1cm).decompose() * u.cm)
return coma
@sbd.dataclass_input(eph=sbd.Ephem)
def to_fluxd(self, wfb, aper, eph, unit=None, **kwargs):
"""Express as spectral flux density in an observation.
Assumes the small angle approximation.
Parameters
----------
wfb : `~astropy.units.Quantity`, `~synphot.SpectralElement`, list
Wavelengths, frequencies, bandpass, or list of
bandpasses of the observation. Bandpasses require
`~synphot`. Ignored if ``S`` is provided.
aper: `~astropy.units.Quantity`, `~sbpy.activity.Aperture`
Aperture of the observation as a circular radius (length
or angular units), or as an sbpy `~sbpy.activity.Aperture`.
eph: dictionary-like, `~sbpy.data.Ephem`
Ephemerides of the comet. Required fields: 'rh', 'delta'.
Optional: 'phase'.
unit : `~astropy.units.Unit`, string, optional
The flux density unit for the output.
"""
# This method handles the geometric quantities. Sub-classes
# will handle the photometric quantities in `_source_fluxd`.
# rho = effective circular aperture radius at the distance of
# the comet. Keep track of array dimensionality as Ephem
# objects can needlessly increase the number of dimensions.
if isinstance(aper, Aperture):
rho = aper.coma_equivalent_radius()
ndim = np.ndim(rho)
else:
rho = aper
ndim = np.ndim(rho)
rho = rho.to("km", sbu.projected_size(eph))
ndim = max(ndim, np.ndim(self))
# validate unit
if unit is not None:
unit = u.Unit(unit)
# get source spectral flux density
# * sunlight for Afrho,
# * blackbody emission for Efrho
# quantity = (delta**2 * F / rho) / source
# must have spectral flux density units
source = self._source_fluxd(wfb, eph, unit=unit, **kwargs)
if isinstance(source, u.Magnitude):
_source = source.physical
else:
_source = source
fluxd = self * rho / eph["delta"] ** 2 * _source
# using Ephem can unnecessarily promote fluxd to an array
if np.ndim(fluxd) > ndim:
fluxd = np.squeeze(fluxd)
# and back to magnitudes, as needed
return fluxd.to(source.unit)
@abc.abstractmethod
def _source_fluxd(self, wfb, eph, unit=None, **kwargs):
"""Photometric calibration of dust coma quantity.
quantity = delta**2 * F / rho / source
delta - observer-comet distance
F - observed spectral flux density
rho - photometric aperture radius at the distance of the comet
source - source function flux density (this method)
For Afrho, source = S / rh**2 / 4 * Phi(phase).
For Efrho, source = 1 / pi / B(T).
Must respect requested units.
"""
[docs]class Afrho(DustComaQuantity):
"""Coma dust quantity for scattered light.
``Afrho`` objects behave like `~astropy.units.Quantity` objects
with units of length.
Parameters
----------
value : number, `~astropy.units.Quantity`
The value(s).
unit : string, `~astropy.units.Unit`, optional
The unit of the input value. Strings must be parseable by
:mod:`~astropy.units` package.
dtype : `~numpy.dtype`, optional
See `~astropy.units.Quantity`.
copy : bool, optional
See `~astropy.units.Quantity`.
Notes
-----
Afρ is the product of dust albedo, dust filling factor, and
circular aperture radius. It is nominally a constant for a
steady-state coma in free expansion. See A'Hearn et al. (1984)
for details.
References
----------
A'Hearn et al. 1984, AJ 89, 579-591.
Examples
--------
>>> from sbpy.activity import Afrho
>>> print(Afrho(1000.0, 'cm'))
1000.0 cm
"""
[docs] @classmethod
def from_fluxd(cls, wfb, fluxd, aper, eph, **kwargs):
return super().from_fluxd(wfb, fluxd, aper, eph, **kwargs)
from_fluxd.__doc__ = (
DustComaQuantity.from_fluxd.__doc__
+ """
Examples
--------
Convert observed flux density to Afρ, with a user-provided
solar flux density for the V-band:
>>> from sbpy.activity import Afrho
>>> import astropy.units as u
>>> from sbpy.calib import solar_fluxd
>>>
>>> solar_fluxd.set({'V': 1869 * u.W / u.m**2 / u.um})
>>>
>>> fluxd = 6.730018324465526e-14 * u.W / u.m**2 / u.um
>>> aper = 1 * u.arcsec
>>> eph = dict(rh=1.5 * u.au, delta=1.0 * u.au)
>>> afrho = Afrho.from_fluxd('V', fluxd, aper, eph=eph)
>>> print(afrho) # doctest: +FLOAT_CMP
999.9999999999999 cm
"""
)
[docs] def to_fluxd(self, wfb, aper, eph, unit=None, phasecor=False, Phi=None):
return super().to_fluxd(wfb, aper, eph, unit=unit, phasecor=phasecor, Phi=Phi)
to_fluxd.__doc__ = (
DustComaQuantity.to_fluxd.__doc__
+ """
phasecor: bool, optional
Scale the result by the phase function ``Phi``, assuming
``Afrho`` is quoted for 0° phase. Requires phase angle in
``eph``.
Phi : callable, optional
Phase function, see :func:`~Afrho.to_phase`.
**kwargs
Keyword arguments for `~Sun.observe`.
Returns
-------
fluxd : `~astropy.units.Quantity`
Spectral flux density.
Examples
--------
>>> from sbpy.activity import Afrho
>>> import astropy.units as u
>>> afrho = Afrho(1000 * u.cm)
>>> wave = 0.55 * u.um
>>> aper = 1 * u.arcsec
>>> eph = dict(rh=1.5 * u.au, delta=1.0 * u.au)
>>> fluxd = afrho.to_fluxd(wave, aper, eph)
>>> print(fluxd) # doctest: +FLOAT_CMP
6.730018324465526e-14 W / (um m2)
With a phase correction:
>>> eph['phase'] = 30 * u.deg
>>> fluxd = afrho.to_fluxd(wave, aper, eph, phasecor=True)
>>> print(fluxd) # doctest: +FLOAT_CMP
2.8017202649540757e-14 W / (um m2)
In magnitudes through the Johnson V filter:
>>> import sbpy.units as sbu
>>> from sbpy.photometry import bandpass
>>> bp = bandpass('Johnson V')
>>> fluxd = afrho.to_fluxd(bp, aper, eph, unit=sbu.JMmag,
... phasecor=True)
>>> print(fluxd) # doctest: +FLOAT_CMP
15.321242371548918 mag(JM)
"""
)
@bib.cite({"model": "1984AJ.....89..579A"})
def _source_fluxd(self, wfb, eph, unit=None, phasecor=False, Phi=None, **kwargs):
# get solar flux density
sun = Sun.from_default()
try:
S = sun.observe(wfb, unit=unit, **kwargs)
except SinglePointSpectrumError:
S = sun(wfb, unit=unit)
if not (
S.unit.is_equivalent(u.W / u.m**2 / u.um)
or S.unit.is_equivalent(u.W / u.m**2 / u.Hz)
or isinstance(S, u.Magnitude)
):
raise ValueError(
"Solar flux density must have units of spectral flux "
"density, e.g., W/m2/μm or W/m2/Hz"
)
if phasecor:
Phi = phase_HalleyMarcus if Phi is None else Phi
_Phi = Phi(eph["phase"]) / Phi(0 * u.deg)
else:
_Phi = 1
# compute
_S = S.physical if isinstance(S, u.Magnitude) else S
source = _S * _Phi / 4 * u.au**2 / eph["rh"] ** 2
return source.to(S.unit)
[docs] def to_phase(self, to_phase, from_phase, Phi=None):
"""Scale to another phase angle.
Parameters
----------
to_phase : `~astropy.units.Quantity`
New target phase angle.
from_phase : `~astropy.units.Quantity`
Current target phase angle.
Phi : callable, optional
Phase function, a callable object that takes a single
parameter, phase angle as a `~astropy.units.Quantity`, and
returns a scale factor. Default is `~phase_HalleyMarcus`.
Returns
-------
afrho : `~Afrho`
The scaled Afρ quantity.
Examples
--------
>>> from sbpy.activity import Afrho
>>> afrho = Afrho(10 * u.cm).to_phase(15 * u.deg, 0 * u.deg)
>>> print(afrho) # doctest: +FLOAT_CMP
5.87201 cm
"""
if Phi is None:
Phi = phase_HalleyMarcus
return self * Phi(to_phase) / Phi(from_phase)
[docs]class Efrho(DustComaQuantity):
"""Coma dust quantity for thermal emission.
``Efrho`` behave like `~astropy.units.Quantity` objects with units
of length.
Parameters
----------
value : number, `~astropy.units.Quantity`
The value(s).
unit : str, `~astropy.units.Unit`, optional
The unit of the input value. Strings must be parseable by
:mod:`~astropy.units` package.
dtype : `~numpy.dtype`, optional
See `~astropy.units.Quantity`.
copy : bool, optional
See `~astropy.units.Quantity`.
Notes
-----
εfρ is the product of dust emissivity, dust filling factor, and
circular aperture radius. It is nominally a constant for a
steady-state coma in free expansion, and is the thermal emission
equivalent for the Afρ quantity. See Kelley et al. (2013) for
details.
References
----------
A'Hearn et al. 1984, AJ 89, 579-591.
Kelley et al. 2013, Icarus 225, 475-494.
Examples
--------
>>> from sbpy.activity import Efrho
>>> print(Efrho(1000.0, 'cm'))
1000.0 cm
"""
[docs] @classmethod
def from_fluxd(cls, wfb, fluxd, aper, eph, **kwargs):
return super().from_fluxd(wfb, fluxd, aper, eph, **kwargs)
from_fluxd.__doc__ = (
DustComaQuantity.from_fluxd.__doc__
+ """
Examples
--------
>>> from sbpy.activity import Efrho
>>> import astropy.units as u
>>> wave = 15.8 * u.um
>>> fluxd = 6.52 * u.mJy
>>> aper = 11.1 * u.arcsec
>>> eph = {'rh': 4.42 * u.au, 'delta': 4.01 * u.au}
>>> efrho = Efrho.from_fluxd(wave, fluxd, aper, eph=eph)
>>> print(efrho) # doctest: +FLOAT_CMP
120.00836963059808 cm
"""
)
[docs] def to_fluxd(self, wfb, aper, eph, unit=None, Tscale=1.1, T=None, B=None):
return super().to_fluxd(wfb, aper, eph, unit=unit, Tscale=Tscale, T=T, B=B)
to_fluxd.__doc__ = (
DustComaQuantity.to_fluxd.__doc__
+ """
Tscale : float, optional
Scale factor for blackbody in LTE with sunlight. Ignored
if ``T`` or ``B`` is provided.
T : `~astropy.units.Quantity`, optional
Blackbody temperature. Ignored if ``B`` is provided.
B : `~astropy.units.Quantity`, optional
Observed spectral flux density from a blackbody sphere,
i.e., pi * Planck function. Overrides ``T`` and
``Tscale``.
Returns
-------
fluxd : `~astropy.units.Quantity`
Spectral flux density.
Examples
--------
>>> from sbpy.activity import Efrho
>>> import astropy.units as u
>>> efrho = Efrho(120.0, 'cm')
>>> freq = 15.8 * u.um
>>> aper = 11.1 * u.arcsec
>>> eph = {'rh': 4.42 * u.au, 'delta': 4.01 * u.au}
>>> fluxd = efrho.to_fluxd(freq, aper, eph=eph, unit='Jy')
>>> print(fluxd) # doctest: +FLOAT_CMP
0.006519545281786034 Jy
"""
)
@bib.cite({"model": "2013Icar..225..475K"})
def _source_fluxd(self, wfb, eph, unit=None, Tscale=1.1, T=None, B=None):
if T is None:
T = Tscale * 278 / np.sqrt(eph["rh"].to("au").value)
if B is None:
BB = BlackbodySource(T)
try:
B = BB.observe(wfb, unit=unit)
except SinglePointSpectrumError:
B = BB(wfb, unit=unit)
else:
if not (
B.unit.is_equivalent(u.W / u.m**2 / u.um)
or B.unit.is_equivalent(u.W / u.m**2 / u.Hz)
or isinstance(B, u.Magnitude)
):
raise ValueError(
"B must be a magnitude or have units of spectral "
"flux density, e.g., W/m2/μm or W/m2/Hz"
)
return B