Merge pull request #95 from PaulHuwe/RCAL-748_SynthPhotometry

RCAL-748 Synthetic Photometry

romanisim/bandpass.py

@@ -5,6 +5,7 @@
 the nominal flat AB spectrum of 3631 Jy.  The ultimate source of this
 information is https://roman.gsfc.nasa.gov/science/WFI_technical.html .
 """
+import numpy as np
 from importlib import resources
 from scipy import integrate
 from astropy import constants
@@ -32,6 +33,9 @@
 galsim2roman_bandpass.update(**{k: k for k in roman2galsim_bandpass})
 roman2galsim_bandpass.update(**{k: k for k in galsim_bandpasses})
 
+# AB Zero Spectral Flux Density
+ABZeroSpFluxDens = 3631e-23 * u.erg / (u.s * u.cm**2 * u.hertz)
+
 
 def read_gsfc_effarea(filename=None):
     """Read an effective area file from Roman.
@@ -87,21 +91,10 @@ def compute_abflux(effarea=None):
     # which astropy then gives a default name col12 to.
     filter_names = [x for x in effarea.dtype.names
                     if x != 'Wave' and 'col' not in x]
-    abfv = 3631e-23 * u.erg / (u.s * u.cm**2 * u.hertz)
+    abfv = ABZeroSpFluxDens
     out = dict()
     for bandpass in filter_names:
-        integrand = abfv * constants.c / (
-            effarea['Wave'] * u.micron)**2  # f_lambda
-        integrand /= constants.h * constants.c / (
-            effarea['Wave'] * u.micron)  # hc/lambda
-        integrand *= effarea[bandpass] * u.m**2  # effective area in filter
-        # integrate.simpson looks like it loses units.  So convert to something
-        # we know about.
-        integrand = integrand.to(1 / (u.s * u.micron)).value
-        zpflux = integrate.simpson(integrand, x=effarea['Wave'])
-        # effarea['Wave'] is in microns, so we're left with a number of counts
-        # per second
-        out[bandpass] = zpflux
+        out[bandpass] = compute_count_rate(flux=abfv, bandpass=bandpass, effarea=effarea)
     return out
 
 
@@ -118,7 +111,7 @@ def get_abflux(bandpass):
     Returns
     -------
     float
-        the zero point flux (photons / s)s
+        the zero point flux (photons / s)
     """
     bandpass = galsim2roman_bandpass.get(bandpass, bandpass)
 
@@ -129,3 +122,62 @@ def get_abflux(bandpass):
         abflux = compute_abflux()
         get_abflux.abflux = abflux
     return abflux[bandpass]
+
+
+def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=None):
+    """Compute the count rate in a given filter, for a specified SED.
+
+    How many photons would an object with SED given by
+    flux deposit in Roman's detectors in a second?
+
+    Parameters
+    ----------
+    flux : float or np.ndarray with shape matching wavedist.
+        Spectral flux density in units of ergs per second * hertz * cm^2
+    bandpass : str
+        the name of the bandpass
+    filename : str
+        filename to read in
+    effarea : astropy.Table.table
+        Table from GSFC with effective areas for each filter.
+    wavedist : numpy.ndarray
+        Array of wavelengths along which spectral flux densities are defined in microns
+
+    Returns
+    -------
+    float
+        the total bandpass flux (photons / s)
+    """
+    # Read in default Roman effective areas from Goddard, if areas not supplied
+    if effarea is None:
+        effarea = read_gsfc_effarea(filename)
+
+    # If wavelength distribution is supplied, interpolate flux and area
+    # over it and the effective area table layout
+    if wavedist is not None:
+        # Ensure that wavedist and flux have the same shape
+        if wavedist.shape != flux.shape:
+            raise ValueError('wavedist and flux must have identical shapes!')
+
+        all_wavel = np.unique(np.concatenate((effarea['Wave'], wavedist)))
+        all_flux = np.interp(all_wavel, wavedist, flux)
+        all_effarea = np.interp(all_wavel, effarea['Wave'], effarea[bandpass])
+    else:
+        all_wavel = effarea['Wave']
+        all_flux = flux
+        all_effarea = effarea[bandpass]
+
+    integrand = all_flux * constants.c / (
+        all_wavel * u.micron)**2  # f_lambda
+    integrand /= constants.h * constants.c / (
+        all_wavel * u.micron)  # hc/lambda
+    integrand *= all_effarea * u.m**2  # effective area in filter
+    # integrate.simpson looks like it loses units.  So convert to something
+    # we know about.
+    integrand = integrand.to(1 / (u.s * u.micron)).value
+
+    zpflux = integrate.simpson(integrand, x=all_wavel)
+    # effarea['Wave'] is in microns, so we're left with a number of counts
+    # per second
+
+    return zpflux

romanisim/tests/test_bandpass.py

@@ -5,14 +5,22 @@
 * get_abflux
 """
 
+import os
 import pytest
+import asdf
+from metrics_logger.decorators import metrics_logger
 import numpy as np
 from romanisim import bandpass
 from astropy import constants
 from astropy.table import Table
 from astropy import units as u
-from scipy.stats import norm
+from scipy.stats import norm, argus
+from romanisim import log
 
+# List of all bandpasses for full spectrum tests
+IFILTLIST = ['F062', 'F087', 'F106', 'F129', 'F158', 'F184', 'F213', 'F146']
+
+# List of select filters with calculated AB fluxes for AB Flux test
 FILTERLIST = ['F062', 'F158', 'F213']
 ABVLIST = [4.938e10, 4.0225e10, 2.55e10]
 
@@ -39,7 +47,7 @@ def test_read_gsfc_effarea(tmpdir_factory):
     assert read_table['Dwarf Planet'][2] == 'Makemake'
 
 
-@pytest.mark.parametrize("filter", FILTERLIST)
+@pytest.mark.parametrize("filter", IFILTLIST)
 def test_compute_abflux(filter):
     # Test calculated abfluxes vs analytical values
 
@@ -64,10 +72,176 @@ def test_compute_abflux(filter):
     gauss_flux = bandpass.compute_abflux(data_table)
 
     # Comparing both fluxes as magnitudes
-    assert np.isclose(np.log10(theo_flux.value), np.log10(gauss_flux[filter]), atol=1.0e-6)
+    assert np.isclose(np.log10(theo_flux.value), np.log10(gauss_flux[filter]), rtol=1.0e-6)
 
 
 @pytest.mark.parametrize("filter, value", zip(FILTERLIST, ABVLIST))
 def test_get_abflux(filter, value):
     # Test that proper results (within 10%) are returned for select bands.
     assert np.isclose(bandpass.get_abflux(filter), value, rtol=1.0e-1)
+
+
+@metrics_logger("DMS233")
+@pytest.mark.soctests
+def test_convert_flux_to_counts():
+    # Define dirac delta wavelength
+    dd_wavel = 1.290 * u.micron
+
+    # Define effective area table
+    effarea = bandpass.read_gsfc_effarea()
+
+    # Define wavelength distribution
+    dd_wavedist = effarea['Wave'] * u.micron
+
+    # Check that the wavelength spacing is constant
+    assert np.all(np.isclose(np.diff(dd_wavedist), np.diff(dd_wavedist)[0], rtol=1.0e-6))
+    wave_bin_width = np.diff(dd_wavedist)[0]
+
+    # Create dirac-delta-like distribution
+    flux = norm.pdf(dd_wavedist, dd_wavel.value, 0.001)
+
+    # Add constant flux
+    flux += 100
+
+    # Rescale
+    flux *= 1.0e-35
+
+    # Add spectral flux density units
+    flux *= u.erg / (u.s * u.cm ** 2 * u.hertz)
+
+    theoretical_flux = {}
+    computed_flux = {}
+    flux_distro = {}
+
+    for filter in IFILTLIST:
+        # Define filter area
+        area = bandpass.read_gsfc_effarea()[filter] * u.m ** 2
+
+        # Define pedestal flux
+        flux_AB_ratio = ((100.0e-35 * u.erg / (u.s * u.cm ** 2 * u.hertz))
+                         / (3631e-23 * u.erg / (u.s * u.cm ** 2 * u.hertz)))
+        theoretical_flux[filter] = bandpass.get_abflux(filter) * flux_AB_ratio / u.s
+
+        # Add delta function flux
+        dd_flux = (1.0e-35 * u.erg / (u.s * u.cm ** 2 * u.hertz * constants.h.to(u.erg * u.s))
+                   * np.interp(1.29, bandpass.read_gsfc_effarea()['Wave'], area) * area.unit).to(1 / u.s)
+
+        theoretical_flux[filter] = theoretical_flux[filter] + dd_flux
+
+        # Computed flux
+        computed_flux[filter] = bandpass.compute_count_rate(flux, filter) / u.s
+
+        # Test that proper results (within 0.2%) are returned for select bands.
+        assert np.isclose(theoretical_flux[filter].value, computed_flux[filter].value, rtol=2.0e-03)
+
+        # Flux distribution for artifacts
+        flux_distro[filter] = (wave_bin_width * (np.divide(np.multiply(area, flux),
+                               dd_wavedist) / constants.h.to(u.erg * u.s))).to(1 / u.s)
+
+    # Create log entry and artifacts
+    log.info('DMS233: integrated over an input spectra in physical units to derive the number of photons / s.')
+
+    artifactdir = os.environ.get('TEST_ARTIFACT_DIR', None)
+    if artifactdir is not None:
+        af = asdf.AsdfFile()
+        af.tree = {'flux_distro': flux_distro,
+                   'theoretical_flux': theoretical_flux,
+                   'computed_flux': computed_flux,
+                   'flux': flux}
+        af.write_to(os.path.join(artifactdir, 'dms233.asdf'))
+
+
+@pytest.mark.parametrize("filter", ['F062', 'F087', 'F106', 'F129', 'F158', 'F184', 'F213', 'F146'])
+def test_AB_convert_flux_to_counts(filter):
+    # AB Zero Test
+    abfv = 3631e-23 * u.erg / (u.s * u.cm ** 2 * u.hertz)
+
+    effarea = bandpass.read_gsfc_effarea()
+    wavedist = effarea['Wave'] * u.micron
+
+    flux = abfv * np.ones(len(wavedist))
+    computed_flux = bandpass.compute_count_rate(flux, filter)
+
+    assert np.isclose(bandpass.get_abflux(filter), computed_flux, rtol=1.0e-6)
+
+
+def test_unevenly_sampled_wavelengths_flux_to_counts():
+    # Get filter response table for theoretical curve
+    effarea = bandpass.read_gsfc_effarea()
+
+    # Define default wavelength distribution
+    wavedist = effarea['Wave'] * u.micron
+    wave_bin_width = wavedist[1] - wavedist[0]
+
+    # Define spectral features
+    # Linear slope
+    flux_slope = np.array([1, 0.5, 0])
+    # Curve
+    flux_arg = argus.pdf(x=wavedist, chi=1, loc=1, scale=1)
+    # Pedestal 1
+    flux_flat1 = np.array([0.66, 0.66])
+    # Dirac Delta function
+    wavedist_dd = np.linspace(2.13 - 0.001, 2.13 + 0.001, 1000)
+    flux_dd = norm.pdf(wavedist_dd, 2.13, 0.001)
+    # Pedestal 2
+    flux_flat2 = np.array([0.33, 0.33])
+
+    # Array to store the theoretical spectral flux density
+    an_flux = np.zeros(len(effarea['Wave']))
+
+    # Linear slope from 400nm to 1000nm
+    arg_start = np.where(wavedist.value == 1)[0][0]
+    an_flux[0:arg_start] = np.arange(start=1, stop=0, step=(-1 / (arg_start)))
+
+    # Define spectral flux array and wavelengths (uneven spacing)
+    total_flux = flux_slope.copy()
+    total_wavedist = np.array([0.4, 0.7, 1])
+
+    # Curve from 1000nm to 2000nm
+    arg_stop = np.where(wavedist.value == 2)[0][0]
+    arg_wavedist = wavedist.value[arg_start + 1:arg_stop]
+    total_flux = np.append(total_flux, flux_arg[np.where(wavedist.value == 1)[0][0] + 1:
+                                                np.where(wavedist.value == 2)[0][0]])
+    total_wavedist = np.append(total_wavedist, arg_wavedist)
+    an_flux[arg_start:arg_stop + 1] = flux_arg[arg_start:arg_stop + 1]
+
+    # Pedestal from 2000nm to 2130nm
+    dd_loc = np.where(wavedist.value == 2.13)[0][0] + 1
+    total_flux = np.append(total_flux, flux_flat1)
+    total_wavedist = np.append(total_wavedist, np.array([2.0, 2.13 - 0.001]))
+    an_flux[arg_stop + 1:dd_loc] = flux_flat1[0]
+
+    # Delta function at 2130nm
+    total_flux = np.append(total_flux, flux_dd)
+    total_wavedist = np.append(total_wavedist, np.array(wavedist_dd))
+    an_flux[dd_loc] = 1 / wave_bin_width.value
+
+    # Pedestal from 2130nm to 2600nm
+    total_flux = np.append(total_flux, flux_flat2)
+    total_wavedist = np.append(total_wavedist, np.array([2.13 + 0.001, 2.6]))
+    an_flux[dd_loc + 1:] = flux_flat2[0]
+
+    # Rescale both spectra
+    total_flux *= 1.0e-35
+    an_flux *= 1.0e-35
+
+    # Add spectral flux density units
+    total_flux *= u.erg / (u.s * u.cm ** 2 * u.hertz)
+    an_flux *= u.erg / (u.s * u.cm ** 2 * u.hertz)
+
+    for filter in IFILTLIST:
+        # Define filter area
+        area = effarea[filter] * u.m ** 2
+
+        # Analytical flux
+        an_counts = (wave_bin_width * (np.divide(np.multiply(area, an_flux), wavedist)
+                                       / constants.h.to(u.erg * u.s))).to(1 / u.s)
+
+        # Sum the flux in the filter
+        an_counts_sum = np.sum(an_counts)
+
+        # Computed flux
+        computed_flux = bandpass.compute_count_rate(flux=total_flux, bandpass=filter, wavedist=total_wavedist)
+
+        # Test that proper results (within 4%) are returned for select bands.
+        assert np.isclose(an_counts_sum.value, computed_flux, rtol=4.0e-2)