Implement simplified IFU mode with cube readout

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "ScopeSim"
-version = "0.10.0a0"
+version = "0.10.0a1"
 description = "Generalised telescope observation simulator"
 license = "GPL-3.0-or-later"
 authors = ["Kieran Leschinski <[email protected]>"]

scopesim/detector/detector.py

@@ -6,16 +6,15 @@
 from ..optics import ImagePlane
 from ..optics.image_plane_utils import (add_imagehdu_to_imagehdu,
                                         sky_wcs_from_det_wcs)
-from ..utils import get_logger, from_currsys, stringify_dict
+from ..utils import get_logger, from_currsys, stringify_dict, zeros_from_header
 
 
 logger = get_logger(__name__)
 
 
 class Detector:
     def __init__(self, header, cmds=None, **kwargs):
-        image = np.zeros((header["NAXIS2"], header["NAXIS1"]))
-        self._hdu = ImageHDU(header=header, data=image)
+        self._hdu = ImageHDU(header=header, data=zeros_from_header(header))
         self.meta = {}
         self.meta.update(header)
         self.meta.update(kwargs)
@@ -32,6 +31,14 @@ def extract_from(self, image_plane, spline_order=1, reset=True):
         self._hdu = add_imagehdu_to_imagehdu(
             image_plane.hdu, self.hdu, spline_order, wcs_suffix="D")
 
+        # HACK: to get sky WCS from image plane into detector...
+        if self._hdu.header["NAXIS"] == 3:
+            sky_wcs = WCS(image_plane.hdu)
+            # HACK: hardcode force back to celestial, this isn't great but will
+            #       have to do for now
+            sky_wcs.wcs.ctype = ["RA---TAN", "DEC--TAN", "WAVE"]
+            self.header.update(sky_wcs.to_header())
+
     def reset(self):
         """Reset internal HDU data to all-zeros-array."""
         # The detector might have been converted to integers by the
@@ -45,7 +52,8 @@ def hdu(self):
 
         pixel_scale = from_currsys("!INST.pixel_scale", self.cmds)
         plate_scale = from_currsys("!INST.plate_scale", self.cmds)
-        if pixel_scale == 0 or plate_scale == 0:
+        if (pixel_scale == 0 or plate_scale == 0
+                or self._hdu.header["NAXIS"] == 3):
             logger.warning("Could not create sky WCS.")
         else:
             sky_wcs, _ = sky_wcs_from_det_wcs(

scopesim/effects/__init__.py

@@ -15,6 +15,7 @@
 from .detector_list import *
 from .electronic import *
 from .shutter import *
+from .binning_3d import *
 from .obs_strategies import *
 from .fits_headers import *
 

scopesim/effects/binning_3d.py

@@ -0,0 +1,57 @@
+# -*- coding: utf-8 -*-
+"""Where else to put this?."""
+
+from typing import ClassVar
+
+from astropy import units as u
+
+from .effects import Effect
+from ..optics.fov import FieldOfView
+from ..utils import get_logger, from_currsys, unit_includes_per_physical_type
+
+
+logger = get_logger(__name__)
+
+
+class FluxBinning3D(Effect):
+    """Takes care of cube flux conversion in absence of a SpectralTraceList."""
+
+    z_order: ClassVar[tuple[int, ...]] = (690,)
+
+    def apply_to(self, fov, **kwargs):
+        """See parent docstring."""
+        if not isinstance(fov, FieldOfView):
+            return fov
+
+        if fov.hdu is None or fov.hdu.header["NAXIS"] != 3:
+            logger.error("Cannot apply cube flux binning.")
+            return fov
+
+        bunit = u.Unit(fov.hdu.header["BUNIT"].lower())
+        pixel_area = fov.pixel_area << u.arcsec**2
+
+        # Spatial binning
+        if unit_includes_per_physical_type(bunit, "solid angle"):
+            fov.hdu.data *= pixel_area.value
+            bunit *= pixel_area.unit
+        else:
+            logger.warning("Cube is already binned spatially.")
+
+        # Spectral binning
+        if unit_includes_per_physical_type(bunit, "length"):
+            fov.hdu.data *= self.dwave.value
+            bunit *= self.dwave.unit
+        else:
+            logger.warning("Cube is already binned spectrally.")
+
+        fov.hdu.header["BUNIT"] = bunit.to_string("fits")
+
+        # This is done in SpectralTraceList as well, idk...
+        fov.cube = fov.hdu
+        return fov
+
+    @property
+    def dwave(self) -> u.Quantity[u.um]:
+        # TODO: turn into class attribute once cmds lookup works...
+        dwave = from_currsys("!SIM.spectral.spectral_bin_width", self.cmds)
+        return dwave << u.um