diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 0efed64..b563bbe 100644 --- a/pahfit/features/features.py +++ b/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): diff --git a/pahfit/fitters/ap_components.py b/pahfit/fitters/ap_components.py index db71949..921e23b 100644 --- a/pahfit/fitters/ap_components.py +++ b/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,147 @@ 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 + ) + + def evaluate(self, 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 = 1 / 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 * self.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 + ) + + def evaluate(self, 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 * self.intensity_amplitude_factor + return Anu * np.exp(-0.5 * np.square((x - mean) / stddev)) diff --git a/pahfit/fitters/ap_fitter.py b/pahfit/fitters/ap_fitter.py index 92760bf..af8be31 100644 --- a/pahfit/fitters/ap_fitter.py +++ b/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, } diff --git a/pahfit/model.py b/pahfit/model.py index 83380b0..70b55ae 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -1,5 +1,6 @@ from specutils import Spectrum1D from astropy import units as u +from astropy import constants import copy import matplotlib as mpl from matplotlib import pyplot as plt @@ -187,8 +188,13 @@ def guess( # parse spectral data self.features.meta["user_unit"]["flux"] = spec.flux.unit _, _, _, lam, flux, _ = self._convert_spec_data(spec, z) - wmin = min(lam) - wmax = max(lam) + lam_min = min(lam) + lam_max = max(lam) + + # Some useful quantities for guessing + median_flux = np.median(flux) + Flambda = flux * units.intensity * (lam * units.wavelength) ** -2 * constants.c + total_power = integrate.trapezoid(Flambda, lam * units.wavelength) # simple linear interpolation function for spectrum sp = interpolate.interp1d(lam, flux) @@ -205,7 +211,7 @@ def loop_over_non_fixed(kind, parameter, estimate_function, force=False): # guess starting point of bb def starlight_guess(row): bb = BlackBody1D(1, row["temperature"][0]) - w = wmin + 0.1 # the wavelength used to compare + w = lam_min + 0.1 # the wavelength used to compare if w < 5: # wavelength is short enough to not have numerical # issues. Evaluate both at w. @@ -227,12 +233,12 @@ def dust_continuum_guess(row): temp = row["temperature"][0] fmax_lam = 2898.0 / temp bb = BlackBody1D(1, temp) - if fmax_lam >= wmin and fmax_lam <= wmax: + if fmax_lam >= lam_min and fmax_lam <= lam_max: lam_ref = fmax_lam - elif fmax_lam > wmax: - lam_ref = wmax + elif fmax_lam > lam_max: + lam_ref = lam_max else: - lam_ref = wmin + lam_ref = lam_min flux_ref = np.median(flux[(lam > lam_ref - 0.2) & (lam < lam_ref + 0.2)]) amp_guess = flux_ref / bb(lam_ref) @@ -241,47 +247,59 @@ def dust_continuum_guess(row): loop_over_non_fixed("dust_continuum", "tau", dust_continuum_guess) def line_fwhm_guess(row): - w = row["wavelength"][0] - if not instrument.within_segment(w, inst): + lam_line = row["wavelength"][0] + if not instrument.within_segment(lam_line, inst): return 0 - fwhm = instrument.fwhm(inst, w, as_bounded=True)[0][0] + fwhm = instrument.fwhm(inst, lam_line, as_bounded=True)[0][0] return fwhm - def amp_guess(row, fwhm): - w = row["wavelength"][0] - if not instrument.within_segment(w, inst): + def power_guess(row, fwhm): + # local integration for the lines + lam_line = row["wavelength"][0] + if not instrument.within_segment(lam_line, inst): return 0 factor = 1.5 - wmin = w - factor * fwhm - wmax = w + factor * fwhm - lam_window = (lam > wmin) & (lam < wmax) + lam_min = lam_line - factor * fwhm + lam_max = lam_line + factor * fwhm + lam_window = (lam > lam_min) & (lam < lam_max) xpoints = lam[lam_window] ypoints = flux[lam_window] if np.count_nonzero(lam_window) >= 2: # difference between flux in window and flux around it - power_guess = integrate.trapezoid(flux[lam_window], lam[lam_window]) + Fnu_dlambda = integrate.trapezoid(flux[lam_window], lam[lam_window]) # subtract continuum estimate, but make sure we don't go negative continuum = (ypoints[0] + ypoints[-1]) / 2 * (xpoints[-1] - xpoints[0]) - if continuum < power_guess: - power_guess -= continuum + if continuum < Fnu_dlambda: + Fnu_dlambda -= continuum else: - power_guess = 0 - return power_guess / fwhm + Fnu_dlambda = 0 + + # this is an unphysical power (Fnu * dlambda), but we + # convert to Fnu dnu = Fnu dnu/dlambda dlambda = Fnu c / + # lambda **2 dlambda + Fnu_dlambda *= units.intensity * units.wavelength + Fnu_dnu = Fnu_dlambda * constants.c / (lam_line * units.wavelength) ** 2 + return Fnu_dnu.to(units.intensity_power).value - # Same logic as in the old function: just use same amp for all - # dust features. - some_flux = 0.5 * np.median(flux) - loop_over_non_fixed("dust_feature", "power", lambda row: some_flux) + def drude_power_guess(row): + # multiply total power by some fraction to guess Drude power + fwhm = row["fwhm"][0] * units.wavelength + delta_w = spec.spectral_axis[-1] - spec.spectral_axis[0] + return (total_power * fwhm / delta_w).to(units.intensity_power).value + + loop_over_non_fixed("dust_feature", "power", drude_power_guess) if integrate_line_flux: - # calc line amplitude using instrumental fwhm and integral over data + # calc line power using instrumental fwhm and integral over data loop_over_non_fixed( - "line", "power", lambda row: amp_guess(row, line_fwhm_guess(row)) + "line", "power", lambda row: power_guess(row, line_fwhm_guess(row)) ) else: - loop_over_non_fixed("line", "power", lambda row: some_flux) + loop_over_non_fixed( + "line", "power", lambda row: median_flux * line_fwhm_guess(row) + ) # Set the fwhms in the features table. Slightly different logic, # as the fwhm for lines are masked by default. TODO: leave FWHM @@ -745,14 +763,16 @@ def tabulate( if wavelengths is None: ranges = instrument.wave_range(instrumentname) if isinstance(ranges[0], float): - wmin, wmax = ranges + lam_min, lam_max = ranges else: # In case of multiple ranges (multiple segments), choose # the min and max instead - wmin = min(r[0] for r in ranges) - wmax = max(r[1] for r in ranges) - wfwhm = instrument.fwhm(instrumentname, wmin, as_bounded=True)[0, 0] - lam = np.arange(wmin, wmax, wfwhm / 2) * u.micron + lam_min = min(r[0] for r in ranges) + lam_max = max(r[1] for r in ranges) + + wfwhm = instrument.fwhm(instrumentname, lam_min, as_bounded=True)[0, 0] + lam = np.arange(lam_min, lam_max, wfwhm / 2) * u.micron + elif isinstance(wavelengths, Spectrum1D): lam = wavelengths.spectral_axis.to(u.micron) / (1 + z) else: @@ -902,13 +922,12 @@ def cleaned(features_tuple3): elif kind == "line": # be careful with lines that have masked FWHM values here - if use_instrument_fwhm or row['fwhm'] is np.ma.masked: + if use_instrument_fwhm or row["fwhm"] is np.ma.masked: # One caveat here: redshift. We do the necessary # adjustment as follows : 1. shift to observed wav; # 2. evaluate fwhm at oberved wav; 3. shift back to # rest frame wav (width in rest frame will be # narrower than observed width) - lam_obs = row["wavelength"]["val"] * (1.0 + redshift) # returned value is tuple (value, min, max). And # min/max are already masked in case of fixed value