From 245b7cea587f25753beefc3bdf77f54cc12897a7 Mon Sep 17 00:00:00 2001
From: PaulHuwe <phuwe@stsci.edu>
Date: Sun, 4 Feb 2024 17:51:33 -0500
Subject: [PATCH 1/7] Initial tests.

---
 romanisim/bandpass.py            | 51 ++++++++++++++-----
 romanisim/tests/test_bandpass.py | 85 +++++++++++++++++++++++++++++++-
 2 files changed, 123 insertions(+), 13 deletions(-)

diff --git a/romanisim/bandpass.py b/romanisim/bandpass.py
index 93d6e210..ee5679cc 100644
--- a/romanisim/bandpass.py
+++ b/romanisim/bandpass.py
@@ -32,6 +32,8 @@
 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.
@@ -81,20 +83,22 @@ 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
+        zpflux = compute_count_rate(thru=abfv, bandpass=bandpass, effarea=effarea)
+
+        # 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
     return out
 
@@ -123,3 +127,26 @@ def get_abflux(bandpass):
         abflux = compute_abflux()
         get_abflux.abflux = abflux
     return abflux[bandpass]
+
+def compute_count_rate(thru=None, bandpass=None, filename=None, effarea=None):
+    # ab_zero = get_abflux(bandpass)
+
+    if not effarea:
+        effarea = read_gsfc_effarea(filename)
+
+    integrand = thru * 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
+    # integrand /= compute_photflam(bandpass).value
+    zpflux = integrate.simpson(integrand, x=effarea['Wave'])
+    # effarea['Wave'] is in microns, so we're left with a number of counts
+    # per second
+
+    return zpflux
+
+
diff --git a/romanisim/tests/test_bandpass.py b/romanisim/tests/test_bandpass.py
index 4d0fc835..42d373c9 100644
--- a/romanisim/tests/test_bandpass.py
+++ b/romanisim/tests/test_bandpass.py
@@ -5,17 +5,20 @@
 * get_abflux
 """
 
+import os
 import pytest
+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 romanisim import log
 
 FILTERLIST = ['F062', 'F158', 'F213']
 ABVLIST = [4.938e10, 4.0225e10, 2.55e10]
-
+IFILTLIST = ['F062', 'F087', 'F106', 'F129', 'F158', 'F184', 'F213'] #, 'F146']
 
 def test_read_gsfc_effarea(tmpdir_factory):
     table_file = str(tmpdir_factory.mktemp("ndata").join("table.csv"))
@@ -71,3 +74,83 @@ def test_compute_abflux(filter):
 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 AB zero flux, and dirac delta wavelength
+    # abfv = 3631e-23 * u.erg / (u.s * u.cm ** 2 * u.hertz)
+    dd_wavel = 1.290 * u.micron
+    effarea = bandpass.read_gsfc_effarea()
+
+    # # wavedist = np.linspace(0.4, 2.6, len(area))
+    # # thru = norm.pdf(wavedist, wavel.value, 0.01)
+    # # thru += 100
+
+    # Create dirac-delta-like distribution for filter response
+    # dd_wavedist = np.linspace(dd_wavel.value - 0.001, dd_wavel.value + 0.001, 1000)
+    # dd_wavedist = np.linspace(0.4, 2.6, len(bandpass.read_gsfc_effarea()['F129'])) * u.micron
+    dd_wavedist = effarea['Wave'] * u.micron
+    wave_bin_width = dd_wavedist[1] - dd_wavedist[0]
+    # thru = norm.pdf(dd_wavedist, dd_wavel.value, 1)
+    thru = norm.pdf(dd_wavedist, dd_wavel.value, 0.001)# * u.erg / (u.s * u.cm ** 2 * u.hertz)
+
+    # Add constant flux
+    thru += 100
+
+    # Rescale
+    thru *= 1.0e-35
+
+    # Add flux units
+    thru *= u.erg / (u.s * u.cm ** 2 * u.hertz)
+
+    theo_flux={}
+    theo_flux_sum={}
+    gauss_flux={}
+    for filter in IFILTLIST:
+        # Define filter area
+        area = bandpass.read_gsfc_effarea()[filter] * u.m ** 2
+
+        # Analytical flux
+        theo_flux[filter] = (wave_bin_width * (np.divide(np.multiply(area, thru), dd_wavedist) / constants.h.to(u.erg * u.s))).to(1 / u.s)
+
+        # Sum the flux in the filter
+        theo_flux_sum[filter] = np.sum(theo_flux[filter])
+        
+        # Computed flux
+        gauss_flux[filter] = bandpass.compute_count_rate(thru, filter)
+
+        # Comparing both fluxes
+        if filter != 'F129':
+            assert np.isclose(theo_flux_sum[filter].value, gauss_flux[filter], atol=1.0e-6)
+        else:
+            # The granualrity of the binning causes greater difference in the
+            # narrow filter containing the dirac delta function
+            assert np.isclose(theo_flux_sum[filter].value, gauss_flux[filter], atol=1.0e-4)
+
+    # Create log entry and artifacts
+    log.info(f'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 = {'theo_flux': theo_flux,
+                   'theo_flux_sum': theo_flux_sum,
+                   'gauss_flux': gauss_flux,
+                   'thru': thru}
+        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
+
+    thru = abfv * np.ones(len(wavedist))
+    gauss_flux = bandpass.compute_count_rate(thru, filter)
+
+    assert np.isclose(bandpass.get_abflux(filter), gauss_flux, atol=1.0e-6)

From ec0cb922850d64dda1293db680c13df7f89990c6 Mon Sep 17 00:00:00 2001
From: PaulHuwe <phuwe@stsci.edu>
Date: Mon, 12 Feb 2024 14:47:59 -0500
Subject: [PATCH 2/7] Added code and test for custom wavelength spacing.

---
 romanisim/bandpass.py            |  73 ++++++++++------
 romanisim/tests/test_bandpass.py | 140 ++++++++++++++++++++++++-------
 2 files changed, 155 insertions(+), 58 deletions(-)

diff --git a/romanisim/bandpass.py b/romanisim/bandpass.py
index 5d1fddbf..2095ec17 100644
--- a/romanisim/bandpass.py
+++ b/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
@@ -35,6 +36,7 @@
 # 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.
 
@@ -92,20 +94,7 @@ def compute_abflux(effarea=None):
     abfv = ABZeroSpFluxDens
     out = dict()
     for bandpass in filter_names:
-        zpflux = compute_count_rate(thru=abfv, bandpass=bandpass, effarea=effarea)
-
-        # 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
 
 
@@ -122,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)
 
@@ -134,25 +123,57 @@ def get_abflux(bandpass):
         get_abflux.abflux = abflux
     return abflux[bandpass]
 
-def compute_count_rate(thru=None, bandpass=None, filename=None, effarea=None):
-    # ab_zero = get_abflux(bandpass)
 
-    if not effarea:
+def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=None):
+    """Compute the AB zero point fluxes for each filter.
+
+    How many photons would a zeroth magnitude AB star deposit in
+    Roman's detectors in a second?
+
+    Parameters
+    ----------
+    flux : float
+        Spectral flux density
+    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
+
+    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)
 
-    integrand = thru * constants.c / (
-        effarea['Wave'] * u.micron)**2  # f_lambda
+    # If wavelength distribution is supplied, interpolate flux and area
+    # over it and the effective area table layout
+    if wavedist is not None:
+        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 / (
-        effarea['Wave'] * u.micron)  # hc/lambda
-    integrand *= effarea[bandpass] * u.m**2  # effective area in filter
+        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
-    # integrand /= compute_photflam(bandpass).value
-    zpflux = integrate.simpson(integrand, x=effarea['Wave'])
+
+    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
-
-
diff --git a/romanisim/tests/test_bandpass.py b/romanisim/tests/test_bandpass.py
index 42d373c9..fe8c2940 100644
--- a/romanisim/tests/test_bandpass.py
+++ b/romanisim/tests/test_bandpass.py
@@ -7,18 +7,23 @@
 
 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]
-IFILTLIST = ['F062', 'F087', 'F106', 'F129', 'F158', 'F184', 'F213'] #, 'F146']
+
 
 def test_read_gsfc_effarea(tmpdir_factory):
     table_file = str(tmpdir_factory.mktemp("ndata").join("table.csv"))
@@ -42,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
 
@@ -67,7 +72,7 @@ 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))
@@ -76,25 +81,17 @@ def test_get_abflux(filter, value):
     assert np.isclose(bandpass.get_abflux(filter), value, rtol=1.0e-1)
 
 
-#metrics_logger("DMS233")
-#pytest.mark.soctests
+@metrics_logger("DMS233")
+@pytest.mark.soctests
 def test_convert_flux_to_counts():
     # Define AB zero flux, and dirac delta wavelength
-    # abfv = 3631e-23 * u.erg / (u.s * u.cm ** 2 * u.hertz)
     dd_wavel = 1.290 * u.micron
     effarea = bandpass.read_gsfc_effarea()
 
-    # # wavedist = np.linspace(0.4, 2.6, len(area))
-    # # thru = norm.pdf(wavedist, wavel.value, 0.01)
-    # # thru += 100
-
     # Create dirac-delta-like distribution for filter response
-    # dd_wavedist = np.linspace(dd_wavel.value - 0.001, dd_wavel.value + 0.001, 1000)
-    # dd_wavedist = np.linspace(0.4, 2.6, len(bandpass.read_gsfc_effarea()['F129'])) * u.micron
     dd_wavedist = effarea['Wave'] * u.micron
     wave_bin_width = dd_wavedist[1] - dd_wavedist[0]
-    # thru = norm.pdf(dd_wavedist, dd_wavel.value, 1)
-    thru = norm.pdf(dd_wavedist, dd_wavel.value, 0.001)# * u.erg / (u.s * u.cm ** 2 * u.hertz)
+    thru = norm.pdf(dd_wavedist, dd_wavel.value, 0.001)
 
     # Add constant flux
     thru += 100
@@ -105,32 +102,29 @@ def test_convert_flux_to_counts():
     # Add flux units
     thru *= u.erg / (u.s * u.cm ** 2 * u.hertz)
 
-    theo_flux={}
-    theo_flux_sum={}
-    gauss_flux={}
+    theo_flux = {}
+    theo_flux_sum = {}
+    gauss_flux = {}
+
     for filter in IFILTLIST:
         # Define filter area
         area = bandpass.read_gsfc_effarea()[filter] * u.m ** 2
 
         # Analytical flux
-        theo_flux[filter] = (wave_bin_width * (np.divide(np.multiply(area, thru), dd_wavedist) / constants.h.to(u.erg * u.s))).to(1 / u.s)
+        theo_flux[filter] = (wave_bin_width * (np.divide(np.multiply(area, thru),
+                             dd_wavedist) / constants.h.to(u.erg * u.s))).to(1 / u.s)
 
         # Sum the flux in the filter
         theo_flux_sum[filter] = np.sum(theo_flux[filter])
-        
+
         # Computed flux
         gauss_flux[filter] = bandpass.compute_count_rate(thru, filter)
 
-        # Comparing both fluxes
-        if filter != 'F129':
-            assert np.isclose(theo_flux_sum[filter].value, gauss_flux[filter], atol=1.0e-6)
-        else:
-            # The granualrity of the binning causes greater difference in the
-            # narrow filter containing the dirac delta function
-            assert np.isclose(theo_flux_sum[filter].value, gauss_flux[filter], atol=1.0e-4)
+        # Test that proper results (within 2%) are returned for select bands.
+        assert np.isclose(theo_flux_sum[filter].value, gauss_flux[filter], rtol=2.0e-02)
 
     # Create log entry and artifacts
-    log.info(f'DMS233: integrated over an input spectra in physical units to derive the number of photons / s.')
+    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:
@@ -140,7 +134,7 @@ def test_convert_flux_to_counts():
                    'gauss_flux': gauss_flux,
                    'thru': thru}
         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):
@@ -150,7 +144,89 @@ def test_AB_convert_flux_to_counts(filter):
     effarea = bandpass.read_gsfc_effarea()
     wavedist = effarea['Wave'] * u.micron
 
-    thru = abfv * np.ones(len(wavedist))
-    gauss_flux = bandpass.compute_count_rate(thru, filter)
+    flux = abfv * np.ones(len(wavedist))
+    gauss_flux = bandpass.compute_count_rate(flux, filter)
+
+    assert np.isclose(bandpass.get_abflux(filter), gauss_flux, rtol=1.0e-6)
+
+
+def test_uneven_area_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
+    theo_flux = np.zeros(len(effarea['Wave']))
+
+    # Linear slope from 400nm to 1000nm
+    arg_start = np.where(wavedist.value == 1)[0][0]
+    theo_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)
+    theo_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]))
+    theo_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))
+    theo_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]))
+    theo_flux[dd_loc + 1:] = flux_flat2[0]
+
+    # Rescale both spectra
+    total_flux *= 1.0e-35
+    theo_flux *= 1.0e-35
+
+    # Add spectral flux density units
+    total_flux *= u.erg / (u.s * u.cm ** 2 * u.hertz)
+    theo_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
+        theo_counts = (wave_bin_width * (np.divide(np.multiply(area, theo_flux), wavedist)
+                                         / constants.h.to(u.erg * u.s))).to(1 / u.s)
+
+        # Sum the flux in the filter
+        theo_counts_sum = np.sum(theo_counts)
+
+        # Computed flux
+        gauss_flux = bandpass.compute_count_rate(flux=total_flux, bandpass=filter, wavedist=total_wavedist)
 
-    assert np.isclose(bandpass.get_abflux(filter), gauss_flux, atol=1.0e-6)
+        # Test that proper results (within 4%) are returned for select bands.
+        assert np.isclose(theo_counts_sum.value, gauss_flux, rtol=4.0e-2)

From c95726c12b9a1b6550702fc210597f53800b1c6a Mon Sep 17 00:00:00 2001
From: PaulHuwe <phuwe@stsci.edu>
Date: Wed, 14 Feb 2024 10:02:30 -0500
Subject: [PATCH 3/7] Addressed PR comments

---
 romanisim/bandpass.py            | 12 ++++--
 romanisim/tests/test_bandpass.py | 72 +++++++++++++++++---------------
 2 files changed, 46 insertions(+), 38 deletions(-)

diff --git a/romanisim/bandpass.py b/romanisim/bandpass.py
index 2095ec17..0e8d0d1f 100644
--- a/romanisim/bandpass.py
+++ b/romanisim/bandpass.py
@@ -127,12 +127,12 @@ def get_abflux(bandpass):
 def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=None):
     """Compute the AB zero point fluxes for each filter.
 
-    How many photons would a zeroth magnitude AB star deposit in
-    Roman's detectors in a second?
+    How many photons would an object with SED given by
+    flux deposit in Roman's detectors in a second.
 
     Parameters
     ----------
-    flux : float
+    flux : float or np.ndarray with shape matching wavedist.
         Spectral flux density
     bandpass : str
         the name of the bandpass
@@ -141,7 +141,7 @@ def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=Non
     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
+        Array of wavelengths along which spectral flux densities are defined in microns
 
     Returns
     -------
@@ -155,6 +155,10 @@ def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=Non
     # 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])
diff --git a/romanisim/tests/test_bandpass.py b/romanisim/tests/test_bandpass.py
index fe8c2940..807e321c 100644
--- a/romanisim/tests/test_bandpass.py
+++ b/romanisim/tests/test_bandpass.py
@@ -58,21 +58,21 @@ def test_compute_abflux(filter):
 
     # Create dirac-delta-like distribution for filter response
     wavedist = np.linspace(wavel.value - 0.001, wavel.value + 0.001, 1000)
-    thru = norm.pdf(wavedist, wavel.value, 0.0001)
+    flux = norm.pdf(wavedist, wavel.value, 0.0001)
 
     # Analytical flux
-    theo_flux = (abfv * area / (constants.h.to(u.erg * u.s) * wavel)).to(1 / (u.s * u.micron))
+    an_flux = (abfv * area / (constants.h.to(u.erg * u.s) * wavel)).to(1 / (u.s * u.micron))
 
     # Table for filter data storage
     data_table = Table()
     data_table['Wave'] = wavedist
-    data_table[filter] = thru
+    data_table[filter] = flux
 
     # Computed flux
     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]), rtol=1.0e-6)
+    assert np.isclose(np.log10(an_flux.value), np.log10(gauss_flux[filter]), rtol=1.0e-6)
 
 
 @pytest.mark.parametrize("filter, value", zip(FILTERLIST, ABVLIST))
@@ -90,38 +90,42 @@ def test_convert_flux_to_counts():
 
     # Create dirac-delta-like distribution for filter response
     dd_wavedist = effarea['Wave'] * u.micron
-    wave_bin_width = dd_wavedist[1] - dd_wavedist[0]
-    thru = norm.pdf(dd_wavedist, dd_wavel.value, 0.001)
+
+    # 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]
+
+    flux = norm.pdf(dd_wavedist, dd_wavel.value, 0.001)
 
     # Add constant flux
-    thru += 100
+    flux += 100
 
     # Rescale
-    thru *= 1.0e-35
+    flux *= 1.0e-35
 
     # Add flux units
-    thru *= u.erg / (u.s * u.cm ** 2 * u.hertz)
+    flux *= u.erg / (u.s * u.cm ** 2 * u.hertz)
 
-    theo_flux = {}
-    theo_flux_sum = {}
-    gauss_flux = {}
+    an_flux = {}
+    an_flux_sum = {}
+    computed_flux = {}
 
     for filter in IFILTLIST:
         # Define filter area
         area = bandpass.read_gsfc_effarea()[filter] * u.m ** 2
 
         # Analytical flux
-        theo_flux[filter] = (wave_bin_width * (np.divide(np.multiply(area, thru),
+        an_flux[filter] = (wave_bin_width * (np.divide(np.multiply(area, flux),
                              dd_wavedist) / constants.h.to(u.erg * u.s))).to(1 / u.s)
 
         # Sum the flux in the filter
-        theo_flux_sum[filter] = np.sum(theo_flux[filter])
+        an_flux_sum[filter] = np.sum(an_flux[filter])
 
         # Computed flux
-        gauss_flux[filter] = bandpass.compute_count_rate(thru, filter)
+        computed_flux[filter] = bandpass.compute_count_rate(flux, filter)
 
         # Test that proper results (within 2%) are returned for select bands.
-        assert np.isclose(theo_flux_sum[filter].value, gauss_flux[filter], rtol=2.0e-02)
+        assert np.isclose(an_flux_sum[filter].value, computed_flux[filter], rtol=2.0e-02)
 
     # Create log entry and artifacts
     log.info('DMS233: integrated over an input spectra in physical units to derive the number of photons / s.')
@@ -129,10 +133,10 @@ def test_convert_flux_to_counts():
     artifactdir = os.environ.get('TEST_ARTIFACT_DIR', None)
     if artifactdir is not None:
         af = asdf.AsdfFile()
-        af.tree = {'theo_flux': theo_flux,
-                   'theo_flux_sum': theo_flux_sum,
-                   'gauss_flux': gauss_flux,
-                   'thru': thru}
+        af.tree = {'an_flux': an_flux,
+                   'an_flux_sum': an_flux_sum,
+                   'computed_flux': computed_flux,
+                   'flux': flux}
         af.write_to(os.path.join(artifactdir, 'dms233.asdf'))
 
 
@@ -145,12 +149,12 @@ def test_AB_convert_flux_to_counts(filter):
     wavedist = effarea['Wave'] * u.micron
 
     flux = abfv * np.ones(len(wavedist))
-    gauss_flux = bandpass.compute_count_rate(flux, filter)
+    computed_flux = bandpass.compute_count_rate(flux, filter)
 
-    assert np.isclose(bandpass.get_abflux(filter), gauss_flux, rtol=1.0e-6)
+    assert np.isclose(bandpass.get_abflux(filter), computed_flux, rtol=1.0e-6)
 
 
-def test_uneven_area_flux_to_counts():
+def test_unevenly_sampled_wavelengths_flux_to_counts():
     # Get filter response table for theoretical curve
     effarea = bandpass.read_gsfc_effarea()
 
@@ -172,11 +176,11 @@ def test_uneven_area_flux_to_counts():
     flux_flat2 = np.array([0.33, 0.33])
 
     # Array to store the theoretical spectral flux density
-    theo_flux = np.zeros(len(effarea['Wave']))
+    an_flux = np.zeros(len(effarea['Wave']))
 
     # Linear slope from 400nm to 1000nm
     arg_start = np.where(wavedist.value == 1)[0][0]
-    theo_flux[0:arg_start] = np.arange(start=1, stop=0, step=(-1 / (arg_start)))
+    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()
@@ -188,45 +192,45 @@ def test_uneven_area_flux_to_counts():
     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)
-    theo_flux[arg_start:arg_stop + 1] = flux_arg[arg_start:arg_stop + 1]
+    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]))
-    theo_flux[arg_stop + 1:dd_loc] = flux_flat1[0]
+    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))
-    theo_flux[dd_loc] = 1 / wave_bin_width.value
+    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]))
-    theo_flux[dd_loc + 1:] = flux_flat2[0]
+    an_flux[dd_loc + 1:] = flux_flat2[0]
 
     # Rescale both spectra
     total_flux *= 1.0e-35
-    theo_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)
-    theo_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
-        theo_counts = (wave_bin_width * (np.divide(np.multiply(area, theo_flux), wavedist)
+        theo_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
         theo_counts_sum = np.sum(theo_counts)
 
         # Computed flux
-        gauss_flux = bandpass.compute_count_rate(flux=total_flux, bandpass=filter, wavedist=total_wavedist)
+        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(theo_counts_sum.value, gauss_flux, rtol=4.0e-2)
+        assert np.isclose(theo_counts_sum.value, computed_flux, rtol=4.0e-2)

From 6c57d231c68477a7372b4a777fdb82e89c39e62a Mon Sep 17 00:00:00 2001
From: PaulHuwe <phuwe@stsci.edu>
Date: Wed, 14 Feb 2024 10:13:03 -0500
Subject: [PATCH 4/7] Addressed PR comments

---
 romanisim/tests/test_bandpass.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/romanisim/tests/test_bandpass.py b/romanisim/tests/test_bandpass.py
index 807e321c..412835d3 100644
--- a/romanisim/tests/test_bandpass.py
+++ b/romanisim/tests/test_bandpass.py
@@ -223,14 +223,14 @@ def test_unevenly_sampled_wavelengths_flux_to_counts():
         area = effarea[filter] * u.m ** 2
 
         # Analytical flux
-        theo_counts = (wave_bin_width * (np.divide(np.multiply(area, an_flux), wavedist)
+        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
-        theo_counts_sum = np.sum(theo_counts)
+        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(theo_counts_sum.value, computed_flux, rtol=4.0e-2)
+        assert np.isclose(an_counts_sum.value, computed_flux, rtol=4.0e-2)

From b1b5139443ce4a12ef07d93ac08665f7ac365bad Mon Sep 17 00:00:00 2001
From: PaulHuwe <phuwe@stsci.edu>
Date: Wed, 14 Feb 2024 10:21:34 -0500
Subject: [PATCH 5/7] Specified expected spectral flux density units.

---
 romanisim/bandpass.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/romanisim/bandpass.py b/romanisim/bandpass.py
index 0e8d0d1f..4a12dc8f 100644
--- a/romanisim/bandpass.py
+++ b/romanisim/bandpass.py
@@ -133,7 +133,7 @@ def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=Non
     Parameters
     ----------
     flux : float or np.ndarray with shape matching wavedist.
-        Spectral flux density
+        Spectral flux density in units of ergs per second * hertz * cm^2
     bandpass : str
         the name of the bandpass
     filename : str

From ee21dc8c91b7a80b24e88121e2ed7097f2383087 Mon Sep 17 00:00:00 2001
From: PaulHuwe <phuwe@stsci.edu>
Date: Wed, 14 Feb 2024 13:41:19 -0500
Subject: [PATCH 6/7] Addressed PR comments. Improved theoretical flux in a
 test.

---
 romanisim/bandpass.py            |  4 +--
 romanisim/tests/test_bandpass.py | 50 +++++++++++++++++++-------------
 2 files changed, 32 insertions(+), 22 deletions(-)

diff --git a/romanisim/bandpass.py b/romanisim/bandpass.py
index 4a12dc8f..d146b613 100644
--- a/romanisim/bandpass.py
+++ b/romanisim/bandpass.py
@@ -125,10 +125,10 @@ def get_abflux(bandpass):
 
 
 def compute_count_rate(flux, bandpass, filename=None, effarea=None, wavedist=None):
-    """Compute the AB zero point fluxes for each filter.
+    """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.
+    flux deposit in Roman's detectors in a second?
 
     Parameters
     ----------
diff --git a/romanisim/tests/test_bandpass.py b/romanisim/tests/test_bandpass.py
index 412835d3..a42c6f8c 100644
--- a/romanisim/tests/test_bandpass.py
+++ b/romanisim/tests/test_bandpass.py
@@ -58,21 +58,21 @@ def test_compute_abflux(filter):
 
     # Create dirac-delta-like distribution for filter response
     wavedist = np.linspace(wavel.value - 0.001, wavel.value + 0.001, 1000)
-    flux = norm.pdf(wavedist, wavel.value, 0.0001)
+    thru = norm.pdf(wavedist, wavel.value, 0.0001)
 
     # Analytical flux
-    an_flux = (abfv * area / (constants.h.to(u.erg * u.s) * wavel)).to(1 / (u.s * u.micron))
+    theo_flux = (abfv * area / (constants.h.to(u.erg * u.s) * wavel)).to(1 / (u.s * u.micron))
 
     # Table for filter data storage
     data_table = Table()
     data_table['Wave'] = wavedist
-    data_table[filter] = flux
+    data_table[filter] = thru
 
     # Computed flux
     gauss_flux = bandpass.compute_abflux(data_table)
 
     # Comparing both fluxes as magnitudes
-    assert np.isclose(np.log10(an_flux.value), np.log10(gauss_flux[filter]), rtol=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))
@@ -84,17 +84,20 @@ def test_get_abflux(filter, value):
 @metrics_logger("DMS233")
 @pytest.mark.soctests
 def test_convert_flux_to_counts():
-    # Define AB zero flux, and dirac delta wavelength
+    # Define dirac delta wavelength
     dd_wavel = 1.290 * u.micron
+
+    # Define effective area table
     effarea = bandpass.read_gsfc_effarea()
 
-    # Create dirac-delta-like distribution for filter response
+    # 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
@@ -103,29 +106,36 @@ def test_convert_flux_to_counts():
     # Rescale
     flux *= 1.0e-35
 
-    # Add flux units
+    # Add spectral flux density units
     flux *= u.erg / (u.s * u.cm ** 2 * u.hertz)
 
-    an_flux = {}
-    an_flux_sum = {}
+    theoretical_flux = {}
     computed_flux = {}
+    flux_distro = {}
 
     for filter in IFILTLIST:
         # Define filter area
         area = bandpass.read_gsfc_effarea()[filter] * u.m ** 2
 
-        # Analytical flux
-        an_flux[filter] = (wave_bin_width * (np.divide(np.multiply(area, flux),
-                             dd_wavedist) / constants.h.to(u.erg * u.s))).to(1 / u.s)
+        # 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
 
-        # Sum the flux in the filter
-        an_flux_sum[filter] = np.sum(an_flux[filter])
+        # 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))
+                   * area[bandpass.read_gsfc_effarea()['Wave'] == 1.29]).to(1 / u.s)
+        theoretical_flux[filter] = theoretical_flux[filter] + dd_flux
 
         # Computed flux
-        computed_flux[filter] = bandpass.compute_count_rate(flux, filter)
+        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)
 
-        # Test that proper results (within 2%) are returned for select bands.
-        assert np.isclose(an_flux_sum[filter].value, computed_flux[filter], rtol=2.0e-02)
+        # 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.')
@@ -133,8 +143,8 @@ def test_convert_flux_to_counts():
     artifactdir = os.environ.get('TEST_ARTIFACT_DIR', None)
     if artifactdir is not None:
         af = asdf.AsdfFile()
-        af.tree = {'an_flux': an_flux,
-                   'an_flux_sum': an_flux_sum,
+        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'))
@@ -224,7 +234,7 @@ def test_unevenly_sampled_wavelengths_flux_to_counts():
 
         # 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)
+                                       / constants.h.to(u.erg * u.s))).to(1 / u.s)
 
         # Sum the flux in the filter
         an_counts_sum = np.sum(an_counts)

From e9483b3f7e76c1412c1578d701ff8bbb1614de9a Mon Sep 17 00:00:00 2001
From: PaulHuwe <phuwe@stsci.edu>
Date: Wed, 14 Feb 2024 15:25:55 -0500
Subject: [PATCH 7/7] Softened wavelength link to table formatting.

---
 romanisim/tests/test_bandpass.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/romanisim/tests/test_bandpass.py b/romanisim/tests/test_bandpass.py
index a42c6f8c..bada80ef 100644
--- a/romanisim/tests/test_bandpass.py
+++ b/romanisim/tests/test_bandpass.py
@@ -124,7 +124,8 @@ def test_convert_flux_to_counts():
 
         # 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))
-                   * area[bandpass.read_gsfc_effarea()['Wave'] == 1.29]).to(1 / 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