Power Drude and Gaussian

pahfit/features/features.py

@@ -41,7 +41,7 @@
 PARAM_UNITS = {'temperature': pahfit.units.temperature,
                'wavelength': pahfit.units.wavelength,
                'fwhm': pahfit.units.wavelength,
-               'power': pahfit.units.intensity}
+               'power': pahfit.units.intensity_power}
 
 
 class UniqueKeyLoader(yaml.SafeLoader):

pahfit/fitters/ap_components.py

@@ -3,6 +3,8 @@
 from astropy.modeling.physical_models import Drude1D
 from astropy.modeling import Fittable1DModel
 from astropy.modeling import Parameter
+from astropy import constants
+from pahfit import units
 
 __all__ = ["BlackBody1D", "ModifiedBlackBody1D", "S07_attenuation", "att_Drude1D"]
 
@@ -20,12 +22,11 @@ class BlackBody1D(Fittable1DModel):
 
     @staticmethod
     def evaluate(x, amplitude, temperature):
-        """
-        """
+        """ """
         return (
             amplitude
             * 3.97289e13
-            / x ** 3
+            / x**3
             / (np.exp(1.4387752e4 / x / temperature) - 1.0)
         )
 
@@ -82,7 +83,9 @@ def kvt(in_x):
         # Extend kvt profile to shorter wavelengths
         if min(in_x) < min(kvt_wav):
             kvt_wav_short = in_x[in_x < min(kvt_wav)]
-            kvt_int_short_tmp = min(kvt_int) * np.exp(2.03 * (kvt_wav_short - min(kvt_wav)))
+            kvt_int_short_tmp = min(kvt_int) * np.exp(
+                2.03 * (kvt_wav_short - min(kvt_wav))
+            )
             # Since kvt_int_shoft_tmp does not reach min(kvt_int),
             # we scale it to stitch it.
             kvt_int_short = kvt_int_short_tmp * (kvt_int[0] / max(kvt_int_short_tmp))
@@ -134,3 +137,146 @@ def evaluate(x, tau, x_0, fwhm):
             profile = Drude1D(amplitude=1.0, fwhm=fwhm, x_0=x_0)
             tau_x = tau * profile(x)
             return (1.0 - np.exp(-1.0 * tau_x)) / tau_x
+
+
+class PowerDrude1D(Fittable1DModel):
+    """Drude profile with amplitude determined by power.
+
+    Special attention is needed for the units. An implementation for
+    this is 'unitful' because 'power' is defined as integral over
+    frequency, while the profile is formulated as a function of
+    wavelength. If we assume the output has unit(flux), then without
+    additional conversions the power unit will be unit(flux) * frequency
+    = unit(flux) * unit(c) / unit(lam).
+
+    Example: if x_0 and fwhm are in micron, and the flux is in MJy / sr,
+    then the unit of the fitted power will be MJy sr-1 * unit(c) /
+    micron, which differs by a constant factor from MJy sr-1 Hz,
+    depending on the chosen unit of c.
+
+    For efficiency and to prevent ambiguity, we assume that all units
+    are the internal pahfit units defined in pahfit.units, and
+    precalculate a conversion factor.
+
+    TODO: We need to check if the flux is 'intensity' or 'flux_density',
+    and assume the power parameter has 'intensity_power' or
+    'flux_density_power' units respectively. For now, only intensity is
+    supported.
+
+    """
+
+    power = Parameter(min=0.0)
+    x_0 = Parameter(min=0.0)
+    fwhm = Parameter(default=1, min=0.0)
+
+    # constant factors in the equation to convert power to amplitude of
+    # the profile.
+    intensity_amplitude_factor = (
+        (2 * units.intensity_power * units.wavelength / (constants.c * np.pi))
+        .to(units.intensity)
+        .value
+    )
+
+    @staticmethod
+    def evaluate(x, power, x_0, fwhm):
+        """Smith, et al. (2007) dust features model. Calculation is for
+        a Drude profile (equation in section 4.1.4).
+
+        The intensity profile as a function of
+        wavelength is
+
+        Inu(lambda) = (b * g**2) / ((lambda / x0 - x0 / lambda)**2 + g**2)
+
+        With
+        b = amplitude (has same unit as flux)
+        g = fwhm / x0
+        x0 = central wavelength
+
+        The integrated power (Fnu dnu) of the profile is
+
+        P = (pi * c / 2) * (b * g / x0)
+
+        Which can be solved for the amplitude b.
+
+        b = (P * 2 * x0) / (pi * c * g)
+
+        According to the above equations, without additional
+        conversions, the resulting amplitude unit will be unit(P) *
+        Hz-1. The factor x0 / c = nu0 (Hz-1) will result in very small
+        values for for Inu(lambda), or very large values for P. To avoid
+        numerical problems, we apply a conversion that ensures
+        Inu(lambda) and P are in internal units.
+
+        Parameters
+        ----------
+        power : float
+        fwhm : float
+        central intensity (x_0) : float
+
+        """
+        # The equation and unit conversion for the amplitude:
+        # b = (2 * P * x_0 / (pi * c * g)).to(output_unit).value
+
+        # Use predetermined factor that deals with units and constants.
+        # factor = (2 * unit(power) * unit(wavelength) / (pi * c)).to(unit(intensity))
+
+        g = fwhm / x_0
+        b = power * x_0 / g * PowerDrude1D.intensity_amplitude_factor
+        return b * g**2 / ((x / x_0 - x_0 / x) ** 2 + g**2)
+
+
+class PowerGaussian1D(Fittable1DModel):
+    """Gaussian profile with amplitude derived from power.
+
+    Implementation analogous to PowerDrude1D.
+
+    The amplitude of a gaussian profile given its power P, is P /
+    (stddev sqrt(2 pi)). Since stddev is given in wavelength units, this
+    equation gives the peak density per wavelength interval, Alambda.
+    The profile Flambda is then
+
+    Flambda(lambda) = Alambda * G(lambda; mean, stddev)
+
+    where G is a gaussian profile with amplitude 1.Converting this to
+    Fnu units yields
+
+    Fnu(lambda) = lambda**2 / c * Flambda(lambda)
+
+    Approximating the lambda**2 factor as a constant (the central
+    wavelength = 'mean'), then yields
+
+    Fnu(lambda) = mean**2 / c * Alambda * G(lambda; mean, stddev)
+
+    In other words, for narrow lines, the per-frequency profile is
+    approximately a gaussian with amplitude
+
+    Anu = P * mean**2 / (c * stddev sqrt(2 pi)).
+
+    So the constant factor we can set is
+    (unit(power) * unit(wavelength)**2 / (c * unit(wavelength) * sqrt(2 pi))).to(intensity)
+
+    """
+
+    power = Parameter(min=0.0)
+    mean = Parameter()
+    stddev = Parameter(default=1, min=0.0)
+
+    intensity_amplitude_factor = (
+        (
+            units.intensity_power
+            * (units.wavelength) ** 2
+            / (constants.c * units.wavelength * np.sqrt(2 * np.pi))
+        )
+        .to(units.intensity)
+        .value
+    )
+
+    @staticmethod
+    def evaluate(x, power, mean, stddev):
+        """Evaluate F_nu(lambda) given the power.
+
+        See class description for equations and unit notes."""
+
+        # amplitude in intensity units
+        Anu = power * mean**2 / stddev * PowerGaussian1D.intensity_amplitude_factor
+        return Anu * np.exp(-0.5 * np.square((x - mean) / stddev))

pahfit/fitters/ap_fitter.py

@@ -5,9 +5,9 @@
     ModifiedBlackBody1D,
     S07_attenuation,
     att_Drude1D,
-    Drude1D,
+    PowerDrude1D,
+    PowerGaussian1D,
 )
-from astropy.modeling.functional_models import Gaussian1D
 from astropy.modeling.fitting import LevMarLSQFitter
 import numpy as np
 
@@ -161,10 +161,10 @@ def add_feature_line(self, name, power, wavelength, fwhm):
 
         kwargs = self._astropy_model_kwargs(
             name,
-            ["amplitude", "mean", "stddev"],
+            ["power", "mean", "stddev"],
             [power, wavelength, fwhm / 2.355],
         )
-        self._add_component(Gaussian1D, **kwargs)
+        self._add_component(PowerGaussian1D, **kwargs)
 
     def add_feature_dust_feature(self, name, power, wavelength, fwhm):
         """Register a PowerDrude1D.
@@ -175,9 +175,9 @@ def add_feature_dust_feature(self, name, power, wavelength, fwhm):
         """
         self.feature_types[name] = "dust_feature"
         kwargs = self._astropy_model_kwargs(
-            name, ["amplitude", "x_0", "fwhm"], [power, wavelength, fwhm]
+            name, ["power", "x_0", "fwhm"], [power, wavelength, fwhm]
         )
-        self._add_component(Drude1D, **kwargs)
+        self._add_component(PowerDrude1D, **kwargs)
 
     def add_feature_attenuation(self, name, tau, model="S07", geometry="screen"):
         """Register the S07 attenuation component.
@@ -286,10 +286,6 @@ def get_result(self, component_name):
         (generally the inverse conversions of what was done in the
         register function).
 
-        NOTE: for now, the return units for "power" are (flux unit) x
-        (micron). Still ambiguous, because the flux unit could be flux
-        density or intensity.
-
         Parameters
         ----------
         component_name : str
@@ -300,7 +296,8 @@ def get_result(self, component_name):
         -------
         dict with Parameters according to the PAHFIT definitions.
 
-        e.g. {'power': converted from amplitude, 'fwhm': converted from
+        e.g., for a feature with amplitude, stddev, and mean parameters:
+        {'power': converted from amplitude, 'fwhm': converted from
         stddev, 'mean': wavelength}
 
         """
@@ -322,13 +319,13 @@ def get_result(self, component_name):
             }
         elif c_type == "line":
             return {
-                "power": component.amplitude.value,
+                "power": component.power.value,
                 "wavelength": component.mean.value,
                 "fwhm": component.stddev.value * 2.355,
             }
         elif c_type == "dust_feature":
             return {
-                "power": component.amplitude.value,
+                "power": component.power.value,
                 "wavelength": component.x_0.value,
                 "fwhm": component.fwhm.value,
             }