HB filter set characteristics and data reduction tools

docs/sbpy/calib.rst

@@ -41,6 +41,9 @@ The names of the built-in sources are stored as an internal array.  They can be
     >>> print(sun)
     <Sun: E490-00a (2014) low resolution reference solar spectrum (Table 4)>
 
+
+.. _sbpy-spectroscopic-calibration:
+
 Controlling the default spectra
 -------------------------------
 
@@ -92,6 +95,8 @@ An example showing how to change the default Vega spectrum:
     ...   # vega.txt in effect
 
 
+.. _sbpy-photometric-calibration:
+
 Photometric calibration (without spectra or `synphot`)
 ------------------------------------------------------
 

docs/sbpy/photometry.rst

@@ -181,6 +181,13 @@ For other bandpasses, obtain the photon-counting relative spectral response curv
 
 See `synphot.spectrum.SpectralElement` for other options and file format details.
 
+
+HB cometary narrow-band filter set
+----------------------------------
+
+The HB filter set is a set of narrow-band filters designed for the characterization of comets at near-UV and optical wavelengths (`Farnham, Schleicher, & A'Hearn 2000 <https://scixplorer.org/abs/2000Icar..147..180F/abstract>`_).  Most filter bandpass and photometric properties are 
+
+
 Reference/API
 -------------
 .. automodapi:: sbpy.photometry

sbpy/photometry/__init__.py

@@ -5,3 +5,4 @@
 from .core import *
 from .iau import *
 from .bandpass import *
+from . import hb

sbpy/photometry/hb/__init__.py

@@ -0,0 +1,10 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+"""
+HB cometary narrow-band filter set.
+
+Filter set characteristics and data reduction tools.
+
+"""
+
+from .filter import *

sbpy/photometry/hb/data/filter.ecsv

@@ -0,0 +1,28 @@
+# %ECSV 1.0
+# ---
+# datatype:
+# - {name: name, datatype: string}
+# - {name: designation, datatype: string}
+# - {name: cwave, unit: AA, datatype: float64, description: center wavelength}
+# - {name: width, unit: AA, datatype: string, subtype: 'float64[4]', description: full-width power points at 80%, 50%, 10%, and 1% of the peak transmission}
+# - {name: Flam0, unit: erg / (cm2 s AA), datatype: float64, description: flux of 0 magnitude star}
+# - {name: XXmBCsun, datatype: float64, description: "solar color, XX-BC"}
+# - {name: gamma_XX_XX, unit: 1 / AA, datatype: float64, description: "fraction of gas emission band captured by filter, normalized by the filter's equivalent width and transmission properties"}
+# - {name: gamma_prime_XX_XX, datatype: float64, description: fraction of gas emission band captured by filter}
+# meta: !!omap
+# - {title: Representative filter characteristics and photometric properties.}
+# - {source: "Farnham, Schleicher, & A'Hearn 2000. The HB Narrowband Comet Filters: Standard Stars and Calibrations.  Icarus 147, 180\u2013\
+#     204."}
+# schema: astropy-2.0
+name designation cwave width Flam0 XXmBCsun gamma_XX_XX gamma_prime_XX_XX
+OH 3090/62 3097 [52,58,68,87] 1.056e-08 1.791 1.698e-2 0.98
+NH 3362/58 3361 [47,54,64,81] 8.42e-09 1.188 1.907e-2 0.99
+UC 3448/84 3449 [72,79,93,116] 7.802e-09 1.101 0 0
+CN 3870/62 3869 [50,56,65,82] 8.6e-09 1.031 1.812e-2 0.99
+C3 4062/62 4063 [43,58,69,84] 8.16e-09 0.497 3.352e-3 0.19
+CO+ 4266/64 4266 [58,64,74,90] 7.323e-09 0.338 1.549e-2 0.99
+BC 4450/67 4453 [55,61,71,86] 6.21e-09 0.000 0 0
+C2 5141/118 5135 [109,119,140,171] 3.887e-09 -0.423 5.433e-3 0.66
+GC 5260/56 5259 [52,56,65,79] 3.616e-09 -0.507 0 0
+H2O+ 7020/170 7028 [148,164,193,239] 1.38e-09 -1.249 5.424e-3 1.00
+RC 7128/58 7133 [53,58,71,92] 1.316e-09 -1.276 0 0

sbpy/photometry/hb/filter.py

@@ -0,0 +1,255 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+"""
+HB cometary narrow-band filter set.
+
+Filter set characteristics and data reduction tools.
+
+"""
+
+__all__ = ["filter_data", "HBFilterSet", "gamma", "gamma_prime"]
+
+from enum import Enum
+from typing import TypeVar
+
+import astropy.units as u
+from astropy.io import ascii
+from astropy.table import QTable
+from astropy.utils.data import get_pkg_data_path
+
+HBFilterSetType = TypeVar("HBFilterSetType", bound="HBFilterSet")
+
+"""HB filter set characteristics and photometric properties."""
+filter_data = {
+    row["name"]: row
+    for row in QTable(ascii.read(get_pkg_data_path("data", "filter.ecsv")))
+}
+
+
+class HBFilterSet(Enum):
+    """Representative filter characteristics and photometric properties.
+
+
+    Examples
+    --------
+
+    Get the wavelength of the BC filter:
+
+    >>> from sbpy.photometry.hb import HBFilterSet
+    >>> HBFilterSet.BC.wavelength
+    <Quantity 4453. Angstrom>
+
+
+    References
+    ----------
+
+    Farnham, Schleicher, & A'Hearn 2000. The HB Narrowband Comet Filters:
+    Standard Stars and Calibrations.  Icarus 147, 180–204.
+
+    """
+
+    OH = "OH"
+    NH = "NH"
+    UC = "UC"
+    CN = "CN"
+    C3 = "C3"
+    COplus = "CO+"
+    BC = "BC"
+    C2 = "C2"
+    GC = "GC"
+    H2Oplus = "H2O+"
+    RC = "RC"
+
+    def __str__(self) -> str:
+        """Filter ID as a string."""
+        return self.value
+
+    def __lt__(self, filter: HBFilterSetType) -> bool:
+        """Compare wavelengths."""
+        return self.wavelength < filter.wavelength
+
+    def __gt__(self, filter: HBFilterSetType) -> bool:
+        """Compare wavelengths."""
+        return self.wavelength > filter.wavelength
+
+    @property
+    def designation(self) -> str:
+        """Filter designation."""
+
+        return filter_data[self.value]["designation"]
+
+    @property
+    def wavelength(self) -> u.Quantity[u.AA]:
+        """Center wavelength."""
+
+        return filter_data[self.value]["cwave"]
+
+    @property
+    def widths(self) -> dict[int, u.Quantity[u.AA]]:
+        """Power point width.
+
+        Returns
+        -------
+        widths : dict
+            Measured full-width power points at 80%, 50%, 10%, and 1% of the
+            peak transmission.
+
+        """
+
+        power_percent = [80, 50, 10, 1]
+
+        return {power_percent[i]: filter_data[self.value]["width"][i] for i in range(4)}
+
+    @property
+    def fluxd0(self) -> u.Quantity[u.Unit("erg / (cm2 s AA)")]:
+        """Spectral flux density of a 0 magnitude star."""
+        return filter_data[self.value]["Flam0"]
+
+    @property
+    def solar_color(self) -> u.Magnitude:
+        """Solar color: [filter]-BC."""
+        return u.Magnitude(filter_data[self.value]["XXmBCsun"])
+
+    def gamma(self, molecule: str) -> u.Quantity[1 / u.AA]:
+        """Fraction of molecule present in this filter, normalized by filter equivalent width.
+
+        Tables VI and VII of Farnham, Schleicher, & A'Hearn 2000.
+
+
+        Parameters
+        ----------
+        molecule : str
+            The molecule name: OH, NH, CN, C3, CO+, C2, or H2O+.
+
+
+        Returns
+        -------
+        gamma : |Quantity|
+
+        """
+
+        if self.value in ["UC", "BC", "GC", "RC"]:
+            raise ValueError("gamma is not defined for the continuum filters")
+        elif self.value == molecule:
+            # Table VI of Farnham, Schleicher, & A'Hearn 2000
+            g = filter_data[self.value]["gamma_XX_XX"]
+        elif molecule == "C3":
+            # Table VII of Farnham, Schleicher, & A'Hearn 2000
+            table_vii = {
+                "NH": 1.433e-5 / u.AA,
+                "CN": 1.427e-3 / u.AA,
+                "CO+": 4.607e-4 / u.AA,
+            }
+            g = table_vii.get(self.value, 0 / u.AA)
+        else:
+            g = 0 / u.AA
+
+        return g.to("1/AA")
+
+    def gamma_prime(self, molecule: str) -> u.Quantity[u.dimensionless_unscaled]:
+        """Fraction of molecule present in this filter.
+
+        Tables VI and VII of Farnham, Schleicher, & A'Hearn 2000.
+
+
+        Parameters
+        ----------
+        molecule : str
+            The molecule name: OH, NH, CN, C3, CO+, C2, or H2O+.
+
+
+        Returns
+        -------
+        gamma_prime : |Quantity|
+
+        """
+
+        if self.value in ["UC", "BC", "GC", "RC"]:
+            raise ValueError("gamma is not defined for the continuum filters")
+        elif self.value == molecule:
+            # Table VI of Farnham, Schleicher, & A'Hearn 2000
+            gp = filter_data[self.value]["gamma_prime_XX_XX"]
+        elif molecule == "C3":
+            data = filter_data[self.value]
+            equivalent_width = data["gamma_prime_XX_XX"] / data["gamma_XX_XX"]
+            gp = self.gamma(molecule) * equivalent_width
+        else:
+            gp = 0
+
+        return u.Quantity(gp, u.dimensionless_unscaled)
+
+
+def gamma(filter: str, molecule: str):
+    """Fraction of molecule present filter, normalized by filter equivalent width.
+
+    Tables VI and VII of Farnham, Schleicher, & A'Hearn 2000.
+
+
+    Parameters
+    ----------
+    filter : str
+        The filter name: OH, NH, CN, C3, CO+, C2, or H2O+.
+
+    molecule : str
+        The molecule name: OH, NH, CN, C3, CO+, C2, or H2O+.
+
+
+    Returns
+    -------
+    gamma : |Quantity|
+
+
+    Examples
+    --------
+
+    >>> from sbpy.photometry import hb
+    >>> hb.gamma("CN", "CN")
+    <Quantity 0.01812 1 / Angstrom>
+
+    Using this function is equivalent to using :func:`HBFilterSet.gamma`:
+
+    >>> from sbpy.photometry.hb import HBFilterSet
+    >>> HBFilterSet.CN.gamma("CN")
+    <Quantity 0.01812 1 / Angstrom>
+
+    """
+
+    return HBFilterSet(filter).gamma(molecule)
+
+
+def gamma_prime(filter: str, molecule: str) -> u.Quantity[u.dimensionless_unscaled]:
+    """Fraction of molecule present in filter.
+
+    Tables VI and VII of Farnham, Schleicher, & A'Hearn 2000.
+
+
+    Parameters
+    ----------
+    filter : str
+        The filter name: OH, NH, CN, C3, CO+, C2, or H2O+.
+
+    molecule : str
+        The molecule name: OH, NH, CN, C3, CO+, C2, or H2O+.
+
+
+    Returns
+    -------
+    gamma_prime : |Quantity|
+
+
+    Examples
+    --------
+
+    >>> from sbpy.photometry import hb
+    >>> hb.gamma_prime("OH", "OH")
+    <Quantity 0.98>
+
+    Using this function is equivalent to using :func:`HBFilterSet.gamma_prime`:
+
+    >>> from sbpy.photometry.hb import HBFilterSet
+    >>> HBFilterSet.OH.gamma_prime("OH")
+    <Quantity 0.98>
+
+    """
+
+    return HBFilterSet(filter).gamma_prime(molecule)