From 6037ec134e44743af235208bb2d081f163c1dd4e Mon Sep 17 00:00:00 2001 From: "Michael S. P. Kelley" Date: Tue, 18 Apr 2023 11:18:08 -0400 Subject: [PATCH] New coma continuum models; convert Afrho/Efrho to/from cross-sectional area. --- sbpy/activity/dust.py | 205 ++++++++++++++++++++++++- sbpy/spectroscopy/coma.py | 220 +++++++++++++++++++++++++++ sbpy/spectroscopy/sources.py | 3 +- sbpy/spectroscopy/tests/test_coma.py | 119 +++++++++++++++ 4 files changed, 544 insertions(+), 3 deletions(-) create mode 100644 sbpy/spectroscopy/coma.py create mode 100644 sbpy/spectroscopy/tests/test_coma.py diff --git a/sbpy/activity/dust.py b/sbpy/activity/dust.py index 424289a5..49df4e6f 100644 --- a/sbpy/activity/dust.py +++ b/sbpy/activity/dust.py @@ -648,6 +648,110 @@ def to_phase(self, to_phase, from_phase, Phi=None): return self * Phi(to_phase) / Phi(from_phase) + @sbd.dataclass_input(eph=sbd.Ephem) + @sbd.quantity_to_dataclass(eph=(sbd.Ephem, "delta")) + def to_cross_section(self, Ap, aper, eph): + """Convert from Afρ to dust geometric cross-sectional area. + + + Parameters + ---------- + Ap : float + Grain geometric albedo. Note that A in the Afρ quantity is ``4 * + Ap`` (A'Hearn et al. 1984). + + 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: `~astropy.units.Quantity`, dictionary-like, `~sbpy.data.Ephem` + Observer-comet distance, or ephemeris data object with "delta" + defined. + + + Returns + ------- + G : `~astropy.units.Quantity` + + + Examples + -------- + >>> afrho = Afrho(1000 * u.cm) + >>> eph = {"rh": 1 * u.au, "delta": 1 * u.au, "phase": 60 * u.deg} + >>> aper = 1 * u.arcsec + >>> afrho.to_cross_section(0.05, aper, eph) + + To account for phase angle, first use ``to_phase``: + >>> afrho.to_phase(0 * u.deg, eph["phase"]).to_cross_section(0.1, aper, eph) + + """ + + # G = pi (rh delta)**2 F_comet / (Ap F_sun) + # Ap = A / 4 + # G = pi 4 * (rh delta)**2 F_comet / (A F_sun) + # G A / pi = 4 * (rh delta)**2 F_comet / F_sun + # + # Afrho = 4 (rh delta)**2 F_comet / (rho F_sun) + # Afrho rho = 4 * (rh delta)**2 F_comet / F_sun = G A / pi + # + # G = Afrho rho pi / A + # G = Afrho rho pi / (4 Ap) + + # 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)) + + G = (self * rho * np.pi / (4 * Ap)).to("km2") + return G + + @classmethod + @sbd.dataclass_input(eph=sbd.Ephem) + @sbd.quantity_to_dataclass(eph=(sbd.Ephem, "delta")) + def from_cross_section(cls, G, Ap, aper, eph): + """Initialize from dust geometric cross-sectional area. + + + Parameters + ---------- + G : `~astropy.units.Quantity` + Dust geometric cross-sectional area. + + Ap : float + Grain geometric albedo. Note that A in the Afρ quantity is ``4 * + Ap`` (A'Hearn et al. 1984). + + 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: `~astropy.units.Quantity`, dictionary-like, `~sbpy.data.Ephem` + Observer-comet distance, or ephemeris data object with "delta" + defined. + + + Returns + ------- + G : `~astropy.units.Quantity` + + + Examples + -------- + >>> eph = {"rh": 1 * u.au, "delta": 1 * u.au} + >>> aper = 1 * u.arcsec + >>> Afrho.from_cross_section(1 * u.km**2, 0.05, aper, eph) + + """ + + G1 = cls(1 * u.cm).to_cross_section(Ap, aper, eph) + return cls((G / G1).to("") * u.cm) + class Efrho(DustComaQuantity): """Coma dust quantity for thermal emission. @@ -662,7 +766,7 @@ class Efrho(DustComaQuantity): The value(s). unit : str, `~astropy.units.Unit`, optional - The unit of the input value. Strings must be parseable by + The unit of the input value. Strings must be parsable by :mod:`~astropy.units` package. dtype : `~numpy.dtype`, optional @@ -780,3 +884,102 @@ def _source_fluxd(self, wfb, eph, unit=None, Tscale=1.1, T=None, B=None): ) return B + + @sbd.dataclass_input(eph=sbd.Ephem) + @sbd.quantity_to_dataclass(eph=(sbd.Ephem, "delta")) + def to_cross_section(self, epsilon, aper, eph): + """Convert from εfρ to dust geometric cross-sectional area. + + + Parameters + ---------- + epsilon : float + Grain emissivity. + + 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: `~astropy.units.Quantity`, dictionary-like, `~sbpy.data.Ephem` + Observer-comet distance, or ephemeris data object with "delta" + defined. + + + Returns + ------- + G : `~astropy.units.Quantity` + + + Examples + -------- + >>> efrho = Efrho(1000 * u.cm) + >>> eph = {"rh": 1 * u.au, "delta": 1 * u.au} + >>> aper = 1 * u.arcsec + >>> efrho.to_cross_section(0.95, aper, eph) + + """ + + # F_comet = G epsilon pi B(T) / delta**2 + # G = F_comet delta**2 / (epsilon pi B(T)) + # G epsilon = F_comet delta**2 / (pi B(T)) + # + # efrho = F_comet rho / (pi rho_angular**2 B(T)) + # = F_comet rho / (pi (rho**2 / delta**2) B(T)) + # = F_comet delta**2 / (pi rho B(T)) + # efrho rho = F_comet delta**2 / (pi B(T)) = G epsilon + # + # G = efrho rho / epsilon + + # 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)) + + G = (self * rho / epsilon).to("km2") + return G + + @classmethod + @sbd.dataclass_input(eph=sbd.Ephem) + @sbd.quantity_to_dataclass(eph=(sbd.Ephem, "delta")) + def from_cross_section(cls, G, epsilon, aper, eph): + """Initialize from dust geometric cross-sectional area. + + + Parameters + ---------- + G : `~astropy.units.Quantity` + Dust geometric cross-sectional area. + + epsilon : float + Grain emissivity. + + 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: `~astropy.units.Quantity`, dictionary-like, `~sbpy.data.Ephem` + Observer-comet distance, or ephemeris data object with "delta" + defined. + + + Returns + ------- + G : `~astropy.units.Quantity` + + + Examples + -------- + >>> eph = {"rh": 1 * u.au, "delta": 1 * u.au} + >>> aper = 1 * u.arcsec + >>> Efrho.from_cross_section(1 * u.km**2, 0.95, aper, eph) + + """ + + G1 = Efrho(1 * u.cm).to_cross_section(epsilon, aper, eph) + return Efrho((G / G1).to("") * u.cm) diff --git a/sbpy/spectroscopy/coma.py b/sbpy/spectroscopy/coma.py new file mode 100644 index 00000000..cccec0ab --- /dev/null +++ b/sbpy/spectroscopy/coma.py @@ -0,0 +1,220 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +"""sbpy Coma Spectroscopy + +Spectrophotometric classes for cometary comae. + +Uses the astropy modeling framework. + +Requires synphot. + +""" + +import numpy as np +import astropy.units as u +from astropy.modeling import Fittable1DModel, Parameter +from .core import SpectralGradient +from .sources import BlackbodySource +from ..activity.dust import Afrho, Efrho +from .. import data as sbd +from .. import units as sbu +from ..calib import Sun + +__all__ = ["Scattered", "Thermal"] + + +class Scattered(Fittable1DModel): + """Light scattered from coma dust. + + Does not account for geometric albedo or phase function: + + F_lambda = cross_section * F_sun * R_lambda(S) / pi / rh**2 / delta**2 + + where R is the spectral reddening function, F_sun is the flux density of the + Sun at 1 au, rh is the heliocentric distance scale factor for the solar + flux density, and delta is the observer-coma distance. + + + Parameters + ---------- + eph: dictionary-like, `~sbpy.data.Ephem` + Ephemerides of the comet. Required fields: 'rh', 'delta'. + + unit : string or `astropy.units.Unit` + Unit of the resultant spectrum (spectral flux density). + + cross_section : float or `~astropy.units.Quantity`, optional + Cross-sectional area of the dust. For float input, km**2 is assumed. + + S : float or `~astropy.units.Quantity`, optional + Spectral gradient at 550 nm. For float input, %/100 nm is assumed. + + """ + + n_inputs = 1 + n_outputs = 1 + input_units = {"cross_section": u.km**2, "S": u.percent / sbu.hundred_nm} + + cross_section = Parameter( + description="cross-sectional area of dust", default=1, unit=u.km**2 + ) + S = Parameter( + description="spectral gradient", default=5, unit=u.percent / sbu.hundred_nm + ) + + @sbd.dataclass_input(eph=sbd.Ephem) + def __init__(self, eph, unit, *args, **kwargs): + self.sun = Sun.from_default() + self.eph = eph + self.unit = u.Unit(unit) + super().__init__(*args, **kwargs) + + def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): + return { + "cross_section": u.km**2, + "S": u.percent / sbu.hundred_nm, + } + + @property + def rh(self): + return self.eph["rh"][0] + + @property + def delta(self): + return self.eph["delta"][0] + + def evaluate(self, wave, cross_section, S): + _wave = u.Quantity(wave, u.um) + _cross_section = u.Quantity(cross_section, u.km**2) + # Keep S as a scalar value + _S = np.atleast_1d(u.Quantity(S, u.percent / sbu.hundred_nm))[0] + + sun = self.sun.redden(SpectralGradient(_S, wave0=550 * u.nm)) + if np.size(_wave) == 1: + F_sun = sun(_wave, unit=self.unit) + else: + F_sun = sun.observe(_wave, unit=self.unit) + + spec = ( + _cross_section + * F_sun + / np.pi + / self.delta**2 + / self.rh.to_value(u.au) ** 2 + ).to(self.unit) + return spec if u.Quantity(wave).unit.is_equivalent(u.um) else spec.value + + # def fit_deriv(self, wave, cross_section, S): + # spec = self.evaluate(wave, cross_section, S) + # dspec_dcross_section = spec / cross_section + # dspec_dS = spec * (wave - 0.55 * u.um) + # return [dspec_dcross_section, dspec_dS] + + @u.quantity_input(wave=u.m) + def afrho(self, wave, aper): + """Convert spectrum to Afρ quantity. + + + Parameters + ---------- + wave : `~astropy.units.Quantity` + Spectral wavelengths. + + 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`. + + """ + + return Afrho.from_fluxd(wave, self(wave), aper, self.eph) + + +class Thermal(Fittable1DModel): + """Thermal emission from coma dust. + + Simple blackbody model: + + F_lambda = cross_section * B_lambda(T) / delta**2 + + where B_lambda is the Planck function, and ``T = Tscale * 278 / sqrt(rh)`` + is the dust temperature. + + + Parameters + ---------- + eph: dictionary-like, `~sbpy.data.Ephem` + Ephemerides of the comet. Required fields: 'rh', 'delta' + + unit : string or `astropy.units.Unit` + Unit of the resultant spectrum (spectral flux density). + + cross_section : float or `~astropy.units.Quantity`, optional + Cross-sectional area of the dust. For float input, km**2 is assumed. + + Tscale : float or `~astropy.units.Quantity`, optional + LTE blackbody sphere temperature scale factor: ``T = Tscale * 278 / + sqrt(rh)``. + + """ + + n_inputs = 1 + n_outputs = 1 + input_units = {"cross_section": u.km**2} + + cross_section = Parameter( + description="cross-sectional area of dust", default=1, unit=u.km**2 + ) + Tscale = Parameter( + description="LTE blackbody sphere temperature scale factor", default=1.0 + ) + + @sbd.dataclass_input(eph=sbd.Ephem) + def __init__(self, eph, unit, *args, **kwargs): + self.eph = eph + self.unit = unit + super().__init__(*args, **kwargs) + + def _parameter_units_for_data_units(self, inputs_unit, outputs_unit): + return { + "cross_section": u.km**2, + "Tscale": "", + } + + @property + def rh(self): + return self.eph["rh"][0] + + @property + def delta(self): + return self.eph["delta"][0] + + def evaluate(self, wave, cross_section, Tscale): + _wave = u.Quantity(wave, u.um) + _cross_section = u.Quantity(cross_section, u.km**2) + T = Tscale * u.Quantity(278, u.K) / np.sqrt(self.rh.to_value("au")) + + B = BlackbodySource(T) + if np.size(_wave) == 1: + F_bb = B(_wave, unit=self.unit) + else: + F_bb = B.observe(_wave, unit=self.unit) + + spec = (_cross_section * F_bb / self.delta**2).to(self.unit) + return spec if u.Quantity(wave).unit.is_equivalent(u.um) else spec.value + + @u.quantity_input(wave=u.m) + def efrho(self, wave, aper): + """Convert spectrum to εfρ quantity. + + + Parameters + ---------- + wave : `~astropy.units.Quantity` + Spectral wavelengths. + + 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`. + + """ + + return Efrho.from_fluxd(wave, self(wave), aper, self.eph, Tscale=self.Tscale) diff --git a/sbpy/spectroscopy/sources.py b/sbpy/spectroscopy/sources.py index 868e069d..feb8bb82 100644 --- a/sbpy/spectroscopy/sources.py +++ b/sbpy/spectroscopy/sources.py @@ -533,8 +533,7 @@ class Reddening(BaseUnitlessSpectrum): @u.quantity_input(S=u.percent / u.um) def __init__(self, S): if getattr(S, 'wave0', None) is None: - raise ValueError("Normalization wavelength in `S` (.wave0) is " - "required by not available.") + raise ValueError("Normalization wavelength in `S` (.wave0) is required.") wv = [1, 2] * S.wave0 df = (S.wave0 * S).to('').value super().__init__( diff --git a/sbpy/spectroscopy/tests/test_coma.py b/sbpy/spectroscopy/tests/test_coma.py new file mode 100644 index 00000000..537cc129 --- /dev/null +++ b/sbpy/spectroscopy/tests/test_coma.py @@ -0,0 +1,119 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import pytest +import numpy as np +import astropy.units as u +from astropy.modeling.models import Scale +import astropy.constants as const +import synphot +from ..sources import BlackbodySource, SpectralSource, SynphotRequired, Reddening +from ..coma import Scattered, Thermal +from ..core import SpectralGradient +from ...activity.dust import Afrho, Efrho +from ...calib import Sun +from ...units import hundred_nm + + +class TestScattered: + def test_init(self): + s = Scattered({"rh": 1 * u.au, "delta": 1 * u.au}, "W/(m2 um)") + assert s.rh == 1 * u.au + assert s.delta == 1 * u.au + assert s.unit == u.Unit("W/(m2 um)") + + def test_evaluate(self): + cross_section = 1 * u.km**2 + S = 0 * u.percent / hundred_nm + + w = [2.0, 2.2] * u.um + eph = {"rh": 1 * u.au, "delta": 1 * u.au} + unit = u.Unit("W/(m2 um)") + + # expected value + sun = Sun.from_default() + expected = ( + cross_section + * sun.observe(w, unit=unit) + / np.pi + / eph["rh"].to_value("au") ** 2 + / eph["delta"] ** 2 + ).to(unit) + + # test multiple wavelengths and no reddening + s = Scattered(eph, unit, cross_section=cross_section, S=S) + assert u.allclose(s(w), expected) + + # test single wavelength and add reddening + w = 2 * u.um + s.S = 1 * u.percent / hundred_nm + R = 1 + ((w - 0.55 * u.um) * s.S).to_value("") + expected = ( + cross_section + * sun(w, unit=unit) + * R + / np.pi + / eph["rh"].to_value("au") ** 2 + / eph["delta"] ** 2 + ).to(unit) + assert u.allclose(s(w), expected) + + def test_afrho(self): + cross_section = 1 * u.km**2 + S = 0 * u.percent / hundred_nm + w = 2.0 * u.um + eph = {"rh": 1 * u.au, "delta": 1 * u.au} + rap = 1 * u.arcsec + unit = u.Unit("W/(m2 um)") + + s = Scattered(eph, unit, cross_section=cross_section, S=S) + a = s.afrho(w, rap) + b = Afrho.from_cross_section(cross_section, 1, rap, eph) + assert u.isclose(a, b) + + +class TestThermal: + def test_init(self): + s = Thermal({"rh": 1 * u.au, "delta": 1 * u.au}, "W/(m2 um)") + assert s.rh == 1 * u.au + assert s.delta == 1 * u.au + assert s.unit == u.Unit("W/(m2 um)") + + def test_evaluate(self): + cross_section = 1 * u.km**2 + Tscale = 1.0 + + w = [20.0, 22.0] * u.um + eph = {"rh": 1 * u.au, "delta": 1 * u.au} + unit = u.Unit("W/(m2 um)") + + # expected value + B = BlackbodySource(278 / np.sqrt(eph["rh"].to_value("au"))) + expected = (cross_section * B.observe(w, unit=unit) / eph["delta"] ** 2).to( + unit + ) + + # test multiple wavelengths and no reddening + t = Thermal(eph, unit, cross_section=cross_section, Tscale=Tscale) + assert u.allclose(t(w), expected) + + # test single wavelength and hotter temperatures + w = 2 * u.um + Tscale = 1.1 + eph = {"rh": 0.8 * u.au, "delta": 1 * u.au} + t = Thermal(eph, unit, cross_section=cross_section, Tscale=Tscale) + B = BlackbodySource(Tscale * 278 / np.sqrt(eph["rh"].to_value("au"))) + expected = (cross_section * B(w, unit=unit) / eph["delta"] ** 2).to(unit) + assert u.allclose(t(w), expected) + + def test_efrho(self): + cross_section = 1 * u.km**2 + Tscale = 1.1 + w = 20.0 * u.um + eph = {"rh": 1 * u.au, "delta": 1 * u.au} + rap = 1 * u.arcsec + unit = u.Unit("W/(m2 um)") + + t = Thermal(eph, unit, cross_section=cross_section, Tscale=Tscale) + a = t.efrho(w, rap) + b = Efrho.from_cross_section(cross_section, 1, rap, eph) + assert u.isclose(a, b)