Use hardware bandpass for sky level

imsim/bandpass.py

@@ -96,6 +96,7 @@ def RubinBandpass(
         case _:
             camera = camera
     tp_path = Path(os.getenv("RUBIN_SIM_DATA_DIR")) / "throughputs"
+    bp_hardware = None
     if airmass is None and camera is None:
         file_name = tp_path / "baseline" / f"total_{band}.dat"
         if not file_name.is_file():
@@ -183,6 +184,12 @@ def RubinBandpass(
     bp = bp.truncate(relative_throughput=1.e-3)
     bp = bp.thin()
     bp = bp.withZeropoint('AB')
+    if bp_hardware is not None:
+        bp_hardware.truncate(relative_throughput=1.e-3)
+        bp_hardware.thin()
+        bp_hardware.withZeropoint('AB')
+        bp.bp_hardware = bp_hardware
+
     return bp
 
 

imsim/sky_model.py

@@ -3,21 +3,21 @@
 import warnings
 import numpy as np
 import galsim
-import pickle
 from galsim.config import InputLoader, RegisterInputType, RegisterValueType
 from scipy.interpolate import RegularGridInterpolator, RectBivariateSpline
 import os
 from .meta_data import data_dir
 
-RUBIN_AREA = 0.25 * np.pi * 649**2  # cm^2
+
+RUBIN_AREA = 0.25 * np.pi * 642.3**2  # cm^2
 
 
 __all__ = ['SkyModel', 'SkyGradient']
 
 
 class SkyModel:
     """Interface to rubin_sim.skybrightness model."""
-    def __init__(self, exptime, mjd, bandpass, eff_area=RUBIN_AREA):
+    def __init__(self, exptime, mjd, bandpass, eff_area=RUBIN_AREA, logger=None):
         """
         Parameters
         ----------
@@ -29,13 +29,20 @@ def __init__(self, exptime, mjd, bandpass, eff_area=RUBIN_AREA):
             Bandpass to use for flux calculation.
         eff_area : `float`
             Collecting area of telescope in cm^2. Default: Rubin value from
-            https://confluence.lsstcorp.org/display/LKB/LSST+Key+Numbers
+            ls.st/SMTN-002
+        logger : Logger object.
         """
         from rubin_sim import skybrightness
+        logger = galsim.config.LoggerWrapper(logger)
         self.exptime = exptime
         self.mjd = mjd
         self.eff_area = eff_area
-        self.bandpass = bandpass
+        if hasattr(bandpass, "bp_hardware"):
+            self.bandpass = bandpass.bp_hardware
+        else:
+            logger.info("A separate hardware bandpass is not available. "
+                        "Using the total bandpass for the sky level.")
+            self.bandpass = bandpass
         self._rubin_sim_sky_model = skybrightness.SkyModel()
 
     def get_sky_level(self, skyCoord):
@@ -112,6 +119,7 @@ def __call__(self, x, y):
         """
         return (self.a*x + self.b*y + self.c)/self.sky_level_center
 
+
 class CCD_Fringing:
     """
     Class generates normalized fringing map.
@@ -164,7 +172,6 @@ def generate_heightfield(self, fractal_dimension=2.5, n=4096):
 
         return np.fft.ifft2(A)
 
-
     def simulate_fringes(self, amp=0.002):
         '''
         # Generate random fringing pattern from a heightfield
@@ -183,7 +190,6 @@ def simulate_fringes(self, amp=0.002):
         #Z  += 1
         return(Z)
 
-
     def fringe_variation_level(self):
         '''
         Function implementing temporal and spatial variation of fringing.
@@ -237,6 +243,7 @@ def calculate_fringe_amplitude(self,x,y,amplitude = 0.002):
 
         return (interp_func((x,y)))
 
+
 class SkyModelLoader(InputLoader):
     """
     Class to load a SkyModel object.
@@ -249,6 +256,7 @@ def getKwargs(self, config, base, logger):
         bandpass, safe1 = galsim.config.BuildBandpass(base['image'], 'bandpass', base, logger)
         safe = safe and safe1
         kwargs['bandpass'] = bandpass
+        kwargs['logger'] = logger
         return kwargs, safe
 
 
@@ -261,5 +269,6 @@ def SkyLevel(config, base, value_type):
     value = sky_model.get_sky_level(base['world_center'])
     return value, False
 
-RegisterInputType('sky_model', SkyModelLoader(SkyModel))
+
+RegisterInputType('sky_model', SkyModelLoader(SkyModel, takes_logger=True))
 RegisterValueType('SkyLevel', SkyLevel, [float], input_type='sky_model')

tests/data/reference_sky_levels.json

@@ -1 +1 @@
-{"u": 966.8740383242452, "g": 9565.16844049534, "r": 20335.64330578381, "i": 39770.05348346839, "z": 66820.68250072759, "y": 62387.30799288845}
+{"u": 1769.5583744933522, "g": 12497.156387235353, "r": 25110.19597937431, "i": 44332.949518546884, "z": 70823.00493285632, "y": 59993.94756134916}

tests/sky_level_reference_values.py

@@ -7,7 +7,7 @@
 import imsim
 
 
-RUBIN_AREA = 0.25 * np.pi * 649**2  # cm^2
+RUBIN_AREA = 0.25 * np.pi * 642.3**2  # cm^2
 
 ra = 54.9348753510528
 dec = -35.8385705255579
@@ -24,7 +24,7 @@
 
 sky_levels = {}
 for band in 'ugrizy':
-    bandpass = imsim.RubinBandpass(band)
+    bandpass = imsim.RubinBandpass(band, camera='LsstCam', det_name='R22_S11').bp_hardware
     wave, spec = sky_model.return_wave_spec()
     lut = galsim.LookupTable(wave, spec[0])
     sed = galsim.SED(lut, wave_type='nm', flux_type='flambda')

tests/test_sky_model.py

@@ -30,7 +30,7 @@ def test_sky_model():
         expected_sky_levels = json.load(fobj)
 
     for band in 'ugrizy':
-        bandpass = RubinBandpass(band)
+        bandpass = RubinBandpass(band, camera='LsstCam', det_name='R22_S11')
         sky_model = SkyModel(exptime, mjd, bandpass)
         sky_level = sky_model.get_sky_level(skyCoord)
         np.testing.assert_approx_equal(sky_level, expected_sky_levels[band],
@@ -41,7 +41,7 @@ def test_sky_model():
     # presumably innocuous reasons (makes a difference ~part per 10000).
     # Test still passes at significant=3.
     for band in 'ugrizy':
-        bandpass = RubinBandpass(band, airmass=1.2)
+        bandpass = RubinBandpass(band, airmass=1.2, camera='LsstCam', det_name='R22_S11')
         sky_model = SkyModel(exptime, mjd, bandpass)
         sky_level = sky_model.get_sky_level(skyCoord)
         np.testing.assert_approx_equal(sky_level, expected_sky_levels[band],