Repository: AstarVienna/ScopeSim

Issue:
"Science" mode for the METIS LMS
The METIS Science Team request a short-cut mode for Scopesim that produces "reduced and reconstructed" data cubes. This requires a few enhancements within Scopesim, but a new mode and new datafiles in the irdb package.
- a `Detector3D` class that can readout a three-dimensional `ImagePlane` (this is produced as an intermediate step in the standard LMS mode), essentially identical to the existing `Detector` object, but with three axes. The pixel scale should be LMS detector scale (8.2 mas) in both spatial directions, the wavelength bin can be taken as the one by the existing `rectify_traces` method. 
- photon noise is trivial (as long as we treat it as independent between voxels). Dark current and readout noise requires information on the number of exposures that should go into the output cube (minimum 6, but should depend on the requested total exposure time), but also on the number of pixels  that go into a voxel.
- It is unclear what PSF to use, but the existing `FieldConstantPSF` class should be sufficient.
- We need a new effect for the line-spread function (LSF). In the extisting LMS mode, the LSF is created automatically when the 3D `ImagePlane` is collapsed to the `XiLamImage`. Here, we should use a pre-defined LMS kernel that is convolved onto the wavelength direction of the input cube (or 3D ImagePlane).
- Most of the transmission curves can be applied as before.  The grating efficiencies might be a bit trickier, as they depend on the settings of prism and echelle grating; however, these can be calculated from the target wavelength as before.

Discussion:
@oczoske: Actually, there is no 3D ImagePlane at all. `MetisLMSSpectralTraceList` operates on `FOVSetupBase` to determine the covered volume, and then on `FieldOfViewBase` to create a cube (unless the `Source` is already a cube), which is then mapped onto the detector. This cube is available as `opttrain.fov_manager.fovs[0].cube`. In any case, the functionality to compute the volume and to create a cube needs to be carried over into the `sci_lms` mode. 

Solution patch:
diff --git a/README.md b/README.md
index 93b97fa0..44a7b845 100644
--- a/README.md
+++ b/README.md
@@ -10,6 +10,7 @@
 [![codecov](https://codecov.io/gh/AstarVienna/ScopeSim/graph/badge.svg)](https://codecov.io/gh/AstarVienna/ScopeSim)
 [![PyPI - Version](https://img.shields.io/pypi/v/ScopeSim)](https://pypi.org/project/ScopeSim/)
 ![PyPI - Python Version](https://img.shields.io/pypi/pyversions/ScopeSim)
+![GitHub Release Date](https://img.shields.io/github/release-date/AstarVienna/ScopeSim)
 
 [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
 [![Citation](https://img.shields.io/badge/DOI-10.1117%2F12.2559784-blue)](https://doi.org/10.1117/12.2559784)
diff --git a/pyproject.toml b/pyproject.toml
index 1b568787..f08d09ce 100644
--- a/pyproject.toml
+++ b/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 <kieran.leschinski@unive.ac.at>"]
diff --git a/scopesim/detector/detector.py b/scopesim/detector/detector.py
index 71eebc49..10b43ca9 100644
--- a/scopesim/detector/detector.py
+++ b/scopesim/detector/detector.py
@@ -6,7 +6,7 @@ from astropy.wcs import WCS
 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__)
@@ -14,8 +14,7 @@ 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 @@ class Detector:
         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 @@ class Detector:
 
         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(
diff --git a/scopesim/effects/__init__.py b/scopesim/effects/__init__.py
index d5772a81..91f471cf 100644
--- a/scopesim/effects/__init__.py
+++ b/scopesim/effects/__init__.py
@@ -15,6 +15,7 @@ from . import ter_curves_utils
 from .detector_list import *
 from .electronic import *
 from .shutter import *
+from .binning_3d import *
 from .obs_strategies import *
 from .fits_headers import *
 
diff --git a/scopesim/effects/binning_3d.py b/scopesim/effects/binning_3d.py
new file mode 100644
index 00000000..d355340a
--- /dev/null
+++ b/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
diff --git a/scopesim/effects/detector_list.py b/scopesim/effects/detector_list.py
index fbc22abc..8ef64a47 100644
--- a/scopesim/effects/detector_list.py
+++ b/scopesim/effects/detector_list.py
@@ -9,23 +9,32 @@ above. The effects here (or one of them) are what ScopeSim "actually needs" to
 have defined in order to run.
 """
 
-import warnings
 from typing import ClassVar
 
 import numpy as np
 from astropy import units as u
+from astropy.io import fits
 from astropy.table import Table
 
 from ..optics.fov_volume_list import FovVolumeList
 from .effects import Effect
-from .apertures import ApertureMask
-from ..optics.image_plane_utils import header_from_list_of_xy, calc_footprint
-from ..utils import (from_currsys, close_loop, figure_factory, array_minmax,
-                     quantity_from_table, unit_from_table, get_logger)
+from ..optics.image_plane_utils import (
+    header_from_list_of_xy,
+    calc_footprint,
+    create_wcs_from_points,
+)
+from ..utils import (
+    from_currsys,
+    close_loop,
+    figure_factory,
+    quantity_from_table,
+    unit_from_table,
+    get_logger,
+)
 
 logger = get_logger(__name__)
 
-__all__ = ["DetectorList", "DetectorWindow"]
+__all__ = ["DetectorList", "DetectorWindow", "DetectorList3D"]
 
 
 class DetectorList(Effect):
@@ -121,14 +130,17 @@ class DetectorList(Effect):
 
     """
 
+    dims = "xy"  # 2 spatial dimensions
     z_order: ClassVar[tuple[int, ...]] = (90, 290, 390, 490)
     report_plot_include: ClassVar[bool] = True
     report_table_include: ClassVar[bool] = True
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
-        params = {"pixel_scale": "!INST.pixel_scale",      # arcsec
-                  "active_detectors": "all",}
+        params = {
+            "pixel_scale": "!INST.pixel_scale",  # arcsec
+            "active_detectors": "all",
+        }
         self.meta.update(params)
         self.meta.update(kwargs)
 
@@ -151,48 +163,115 @@ class DetectorList(Effect):
         if not isinstance(obj, FovVolumeList):
             return obj
 
-        hdr = self.image_plane_header
-        xy_mm = calc_footprint(hdr, "D")
-        pixel_size = hdr["CDELT1D"]  # mm
-        pixel_scale = kwargs.get("pixel_scale", self.meta["pixel_scale"])  # ["]
-        pixel_scale = from_currsys(pixel_scale, self.cmds)
+        logger.debug("apply_to got kwargs: %s", kwargs)
+        shrink_axis, shrink_values = self._get_fov_limits(
+            pixel_scale=kwargs.get("pixel_scale", None))
+        obj.shrink(axis=shrink_axis, values=shrink_values)
+
+        return obj
 
-        # x["] = x[mm] * ["] / [mm]
-        xy_sky = xy_mm * pixel_scale / pixel_size
+    @property
+    def pixel_size(self) -> u.Quantity[u.mm]:
+        """Return size of one pixel in mm."""
+        pixel_size = np.min(
+            quantity_from_table("pixel_size", self.active_table, u.mm))
+        return pixel_size
 
-        obj.shrink(axis=["x", "y"],
-                   values=(tuple(zip(xy_sky.min(axis=0),
-                                     xy_sky.max(axis=0)))))
+    @property
+    def pixel_scale_mm(self) -> u.Equivalency:
+        """Return pixel scale (mm / pix) as equivalency."""
+        return u.pixel_scale(self.pixel_size / u.pixel)
+
+    def _get_pixel_scale_arcsec(self, pixel_scale=None) -> u.Equivalency:
+        """To allow overriding defaults, but use defaults in property."""
+        pixel_scale = u.pixel_scale(
+            from_currsys(
+                pixel_scale if pixel_scale is not None else self.meta["pixel_scale"],
+                self.cmds,
+            ) << u.arcsec / u.pixel
+        )
+        return pixel_scale
 
-        return obj
+    @property
+    def pixel_scale_arcsec(self) -> u.Equivalency:
+        """Return pixel scale (arcsec / pix) as equivalency."""
+        return self._get_pixel_scale_arcsec()
+
+    def _get_fov_limits(self, pixel_scale=None):
+        mm_points = self._get_corner_points()[:, :2]
+        pixel_scale = self._get_pixel_scale_arcsec(pixel_scale)
+
+        # The line below combines the pixesl <=> arcsec scale with the
+        # pixel <=> mm scale (both are astropy equivalencies, which can be
+        # simply added together), to allow an indirect conversion from mm via
+        # pixel to arcsec, converting the detector's mm footprint to an on-sky
+        # FOV limit. Conceptually, this should involve the plate scale, but
+        # ScopeSim is still somewhat inconsistent / redundant in defining
+        # pixel size, pixel scale and plate scale.
+        # TODO: Consider using plate scale here once that's generally been
+        #       refactored to be more consistent everywhere.
+        with u.set_enabled_equivalencies(pixel_scale + self.pixel_scale_mm):
+            sky_points = (mm_points << u.pixel).to_value(u.arcsec)
+
+        axis = ["x", "y"]
+        values = (tuple(zip(sky_points.min(axis=0), sky_points.max(axis=0))))
+        return axis, values
+
+    def _get_corner_points(self):
+        tbl = self.active_table
+        with u.set_enabled_equivalencies(self.pixel_scale_mm):
+            cens = {
+                dim: (tbl[f"{dim}_cen"].data.astype(float) *
+                      unit_from_table(f"{dim}_cen", tbl, u.mm) << u.mm)
+                for dim in self.dims
+            }
+
+            ds = {
+                dim: (0.5 * tbl[f"{dim}_size"].data.astype(float) *
+                      unit_from_table(f"{dim}_size", tbl, u.mm) << u.mm)
+                for dim in self.dims
+            }
+
+        det_mins = {dim: np.min(cens[dim] - ds[dim]) for dim in self.dims}
+        det_maxs = {dim: np.max(cens[dim] + ds[dim]) for dim in self.dims}
+
+        points = np.array([
+            [det_mins[dim].to_value(u.mm) for dim in self.dims],
+            [det_maxs[dim].to_value(u.mm) for dim in self.dims],
+        ])
+
+        return points * u.mm
 
     def fov_grid(self, which="edges", **kwargs):
         """Return an ApertureMask object. kwargs are "pixel_scale" [arcsec]."""
-        warnings.warn("The fov_grid method is deprecated and will be removed "
-                      "in a future release.", DeprecationWarning, stacklevel=2)
-        aperture_mask = None
-        if which == "edges":
-            self.meta.update(kwargs)
-            self.meta = from_currsys(self.meta, self.cmds)
-
-            hdr = self.image_plane_header
-            xy_mm = calc_footprint(hdr, "D")
-            pixel_size = hdr["CDELT1D"]              # mm
-            pixel_scale = self.meta["pixel_scale"]   # ["]
+        # This has really been taken care of by apply_to now
+        # TODO v1.0: finally rm this completely
+        raise AttributeError("The DetectorList.fov_grid() method has been "
+                             "removed in vPLACEHOLDER_NEXT_RELEASE_VERSION.")
+        # aperture_mask = None
+        # if which == "edges":
+        #     self.meta.update(kwargs)
+        #     self.meta = from_currsys(self.meta, self.cmds)
+
+        #     hdr = self.image_plane_header
+        #     xy_mm = calc_footprint(hdr, "D")
+        #     pixel_size = hdr["CDELT1D"]              # mm
+        #     pixel_scale = self.meta["pixel_scale"]   # ["]
 
-            # x["] = x[mm] * ["] / [mm]
-            xy_sky = xy_mm * pixel_scale / pixel_size
+        #     # x["] = x[mm] * ["] / [mm]
+        #     xy_sky = xy_mm * pixel_scale / pixel_size
 
-            aperture_mask = ApertureMask(array_dict={"x": xy_sky[:, 0],
-                                                     "y": xy_sky[:, 1]},
-                                         pixel_scale=pixel_scale)
+        #     aperture_mask = ApertureMask(array_dict={"x": xy_sky[:, 0],
+        #                                              "y": xy_sky[:, 1]},
+        #                                  pixel_scale=pixel_scale)
 
-        return aperture_mask
+        # return aperture_mask
 
     @property
     def image_plane_header(self):
         """Create and return the Image Plane Header."""
         # FIXME: Heavy property.....
+        # TODO: refactor this using (parts of) 3D implementation below
         tbl = self.active_table
         pixel_size = np.min(quantity_from_table("pixel_size", tbl, u.mm))
         x_size_unit = unit_from_table("x_size", tbl, u.mm)
@@ -246,6 +325,8 @@ class DetectorList(Effect):
 
     def detector_headers(self, ids=None):
         """Create detector headers from active detectors or given IDs."""
+        # Note: the ids argument is not used anywhere in ScopeSim.
+
         if ids is not None and all(isinstance(ii, int) for ii in ids):
             self.meta["active_detectors"] = list(ids)
 
@@ -360,3 +441,73 @@ class DetectorWindow(DetectorList):
         }
 
         super().__init__(array_dict=array_dict, **params)
+
+
+class DetectorList3D(DetectorList):
+    """Pseudo-detector for simple IFU mode.
+
+    .. versionadded:: PLACEHOLDER_NEXT_RELEASE_VERSION
+    """
+
+    dims = "xyz"  # 2 spatial dimensions, 1 spectral dimension
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        params = {
+            "pixel_scale": "!INST.pixel_scale",  # arcsec
+            "active_detectors": "all",
+            "wavelen": "!OBS.wavelen",
+            "dwave": "!SIM.spectral.spectral_bin_width",
+        }
+        self.meta.update(params)
+        self.meta.update(kwargs)
+
+    @property
+    def waverange(self) -> u.Quantity[u.um]:
+        """Return wavelength range in um [wave_min, wave_max]."""
+        # TODO: dynamic wave unit ?
+        wave_points = self._get_corner_points()[:, 2]
+        wave_mid = from_currsys(self.meta["wavelen"], self.cmds) * u.um
+        dwave = from_currsys(self.meta["dwave"], self.cmds) * u.um / u.pixel
+
+        with u.set_enabled_equivalencies(self.pixel_scale_mm):
+            waverange = wave_points.to(u.pixel) * dwave + wave_mid
+        return waverange
+
+    def _get_fov_limits(self, pixel_scale=None):
+        mm_points = self._get_corner_points()[:, :2]
+        pixel_scale = self._get_pixel_scale_arcsec(pixel_scale)
+
+        # See parent class for explanation of this operation.
+        with u.set_enabled_equivalencies(pixel_scale + self.pixel_scale_mm):
+            sky_points = (mm_points << u.pixel).to_value(u.arcsec)
+
+        axis = ["x", "y", "wave"]
+        values = (
+            *zip(sky_points.min(axis=0), sky_points.max(axis=0)),
+            tuple(self.waverange.to_value(u.um).round(7))
+        )
+
+        return axis, values  # [arcsec, arcsec, um]
+
+    @property
+    def image_plane_header(self):
+        """Create and return the Image Plane Header."""
+        # FIXME: Heavy property.....
+        points = self._get_corner_points().value
+        new_wcs, naxis = create_wcs_from_points(
+            points, self.pixel_size.to_value(u.mm), "D")
+
+        hdr = fits.Header()
+        hdr["NAXIS"] = 3
+        hdr["NAXIS1"] = int(naxis[0])
+        hdr["NAXIS2"] = int(naxis[1])
+        hdr["NAXIS3"] = int(naxis[2])
+        hdr.update(new_wcs.to_header())
+        hdr["IMGPLANE"] = self.image_plane_id
+
+        return hdr
+
+    def detector_headers(self, ids=None):
+        """Override for simplified 3D case."""
+        return [self.image_plane_header]
diff --git a/scopesim/effects/effects_utils.py b/scopesim/effects/effects_utils.py
index ae3e556e..151bc055 100644
--- a/scopesim/effects/effects_utils.py
+++ b/scopesim/effects/effects_utils.py
@@ -69,8 +69,11 @@ def make_effect(effect_dict, cmds=None, **properties):
 
 
 def is_spectroscope(effects):
-    spec_classes = (efs.SpectralTraceList,
-                    efs.SpectralTraceListWheel)
+    spec_classes = (
+        efs.SpectralTraceList,
+        efs.SpectralTraceListWheel,
+        efs.DetectorList3D,
+    )
     return any(isinstance(eff, spec_classes) for eff in effects)
 
 
diff --git a/scopesim/effects/psfs/discrete.py b/scopesim/effects/psfs/discrete.py
index f759714a..1939cf24 100644
--- a/scopesim/effects/psfs/discrete.py
+++ b/scopesim/effects/psfs/discrete.py
@@ -179,6 +179,8 @@ class FieldConstantPSF(DiscretePSF):
 
         xworld, yworld = cubewcs.all_pix2world(xcube, ycube, 1)
         outcube = np.zeros((lam.shape[0], nypsf, nxpsf), dtype=np.float32)
+        logger.info("Interpolating PSF onto %s cube", outcube.shape)
+
         for i, wave in enumerate(lam):
             psf_wave_pixscale = ref_pixel_scale * wave / refwave
             psfwcs.wcs.cdelt = [psf_wave_pixscale,
diff --git a/scopesim/optics/fov.py b/scopesim/optics/fov.py
index c9c7dc65..742dba0c 100644
--- a/scopesim/optics/fov.py
+++ b/scopesim/optics/fov.py
@@ -366,6 +366,14 @@ class FieldOfView:
         #     np.testing.assert_array_equal(roundtrip[1], new_naxis - [1, 1])
         # except AssertionError:
         #     logger.exception("WCS roundtrip assertion failed.")
+        # FIXME: Related to the above, this sometimes fails:
+        # np.testing.assert_equal(xy1p - xy0p, new_naxis)
+        # This occurs when the floor and ceil in xy0s and xy1s produce an
+        # off-by-one error. Using .round instead for both would solve things
+        # in those cases, but breaks in other cases. Find a proper solution!
+        # Note: This is not super fatal, because the resulting projections
+        #       will trim off that extra pixel later on, but this should still
+        #       be addressed.
 
         new_hdr = new_wcs.to_header()
         new_hdr.update({"NAXIS1": new_naxis[0], "NAXIS2": new_naxis[1]})
@@ -432,7 +440,7 @@ class FieldOfView:
             "NAXIS": 3,
             "NAXIS3": data.shape[0],
             "CRVAL3": round(hdu_waves[i0p], 11),
-            "CRPIX3": 0,
+            "CRPIX3": 1,  # 1 because FITS
             "CDELT3": hdr["CDELT3"],
             "CUNIT3": hdr["CUNIT3"],
             "CTYPE3": hdr["CTYPE3"],
@@ -833,7 +841,7 @@ class FieldOfView:
         Returns
         -------
         canvas_cube_hdu : fits.ImageHDU
-            [ph s-1 AA-1 arcsec-2]      # as needed by SpectralTrace
+            [ph s-1 um-1 arcsec-2]      # as needed by SpectralTrace
 
         """
         # 1. Make waveset and canvas cube (area, bin_width are applied at end)
@@ -846,7 +854,17 @@ class FieldOfView:
         fov_waveset = np.arange(
             self.meta["wave_min"].value, self.meta["wave_max"].value,
             from_currsys("!SIM.spectral.spectral_bin_width", self.cmds)) * wave_unit
-        fov_waveset = fov_waveset.to(u.um)
+
+        # Note: There is an edge case where some floating point madness
+        #       results in the arange ticking over by another bin by just
+        #       .000000000005 (yeah...), which creates and off-by-one error
+        #       further down the line.
+        # TODO: Find a better way to solve this, perhaps with linspace...
+        n_wave = int(
+            (self.meta["wave_max"].value - self.meta["wave_min"].value) /
+            from_currsys("!SIM.spectral.spectral_bin_width", self.cmds)
+        )
+        fov_waveset = fov_waveset[:n_wave].to(u.um)
 
         # TODO: what's with this code??
         # fov_waveset = self.waveset
@@ -864,7 +882,26 @@ class FieldOfView:
                            self.header["NAXIS2"],
                            self.header["NAXIS1"])),
             header=self.header)
-        # canvas_cube_hdu.header["BUNIT"] = "ph s-1 cm-2 AA-1"
+        # set BUNIT initially to PHOTLAM / arcsec**2
+        canvas_cube_hdu.header["BUNIT"] = "ph cm-2 s-1 AA-1 arcsec-2"
+
+        canvas_cube_hdu.header.update({
+            "CDELT3": np.diff(fov_waveset[:2])[0].to_value(u.um),
+            "CRVAL3": fov_waveset[0].value,
+            "CRPIX3": 1,
+            "CUNIT3": "um",
+            "CTYPE3": "WAVE",
+        })
+        # TODO: Add the log wavelength keyword here, if log scale is needed
+
+        if (self.detector_header is not None and
+                self.detector_header["NAXIS"] == 3):  # cube ifu mode
+            # TODO: Find out why not always do that, probably related to this
+            #       "!INST.decouple_detector_from_sky_headers" thingy
+            new_det_wcs = WCS(self.detector_header, key="D")
+            canvas_cube_hdu.header.update(new_det_wcs.to_header())
+
+            assert canvas_cube_hdu.header["NAXIS3"] == self.detector_header["NAXIS3"]
 
         for field_hdu in self._make_cube_cubefields(fov_waveset):
             canvas_cube_hdu = imp_utils.add_imagehdu_to_imagehdu(
@@ -891,20 +928,13 @@ class FieldOfView:
         # SpectralTrace wants ph/s/um/arcsec2 --> get rid of m2, leave um
         canvas_cube_hdu.data *= self.area.to(u.cm ** 2).value
         canvas_cube_hdu.data *= 1e4       # ph/s/AA/arcsec2 --> ph/s/um/arcsec2
+        canvas_cube_hdu.header["BUNIT"] = "ph s-1 um-1 arcsec-2"
 
         # TODO: what's with this code??
         # bin_widths = np.diff(fov_waveset).to(u.AA).value
         # bin_widths = 0.5 * (np.r_[0, bin_widths] + np.r_[bin_widths, 0])
         # canvas_cube_hdu.data *= bin_widths[:, None, None]
 
-        cdelt3 = np.diff(fov_waveset[:2])[0]
-        canvas_cube_hdu.header.update({"CDELT3": cdelt3.to(u.um).value,
-                                       "CRVAL3": fov_waveset[0].value,
-                                       "CRPIX3": 1,
-                                       "CUNIT3": "um",
-                                       "CTYPE3": "WAVE"})
-        canvas_cube_hdu.header["BUNIT"] = "ph s-1 um-1 arcsec-2"
-        # TODO: Add the log wavelength keyword here, if log scale is needed
         return canvas_cube_hdu      # [ph s-1 um-1 (arcsec-2)]
 
     @property
diff --git a/scopesim/optics/fov_manager.py b/scopesim/optics/fov_manager.py
index 2842db09..bf163f84 100644
--- a/scopesim/optics/fov_manager.py
+++ b/scopesim/optics/fov_manager.py
@@ -169,10 +169,12 @@ class FOVManager:
 
             # useful for spectroscopy mode where slit dimensions is not the same
             # as detector dimensions
-            # ..todo: Make sure this changes for multiple image planes
+            # TODO: Make sure this changes for multiple image planes
             if from_currsys(self.meta["decouple_sky_det_hdrs"], self.cmds):
                 det_eff = eu.get_all_effects(self.effects, DetectorList)[0]
                 dethdr = det_eff.image_plane_header
+                # TODO: Why is this .image_plane_header and not
+                #       .detector_headers()[0] or something?
 
             yield FieldOfView(skyhdr,
                               waverange,
diff --git a/scopesim/optics/image_plane.py b/scopesim/optics/image_plane.py
index 8d361e61..194a34ee 100644
--- a/scopesim/optics/image_plane.py
+++ b/scopesim/optics/image_plane.py
@@ -11,7 +11,13 @@ from astropy.wcs import WCS
 
 from .image_plane_utils import add_table_to_imagehdu, add_imagehdu_to_imagehdu
 
-from ..utils import from_currsys, has_needed_keywords, get_logger
+from ..utils import (
+    from_currsys,
+    has_needed_keywords,
+    get_logger,
+    zeros_from_header,
+)
+
 
 logger = get_logger(__name__)
 
@@ -53,16 +59,16 @@ class ImagePlane:
         self.cmds = cmds
         self.meta = {}
         self.meta.update(kwargs)
-        self.id = header["IMGPLANE"] if "IMGPLANE" in header else 0
+        self.id = header.get("IMGPLANE", 0)
 
         if not any(has_needed_keywords(header, s)
                    for s in ["", "D", "S"]):
             raise ValueError(f"header must have a valid image-plane WCS: "
                              f"{dict(header)}")
 
-        # image = np.zeros((header["NAXIS2"]+1, header["NAXIS1"]+1))
-        image = np.zeros((header["NAXIS2"], header["NAXIS1"]))
+        image = zeros_from_header(header)
         self.hdu = fits.ImageHDU(data=image, header=header)
+        self.hdu.header["BUNIT"] = "ph s-1"  # photons per second (per pixel)
 
         self._det_wcs = self._get_wcs(header, "D")
         logger.debug("det %s", self._det_wcs)
@@ -139,6 +145,12 @@ class ImagePlane:
                              hdus_or_tables.data.shape)
                 self.hdu = add_imagehdu_to_imagehdu(hdus_or_tables, self.hdu,
                                                     spline_order, wcs_suffix)
+                if (img_bunit := hdus_or_tables.header.get("BUNIT")) is not None:
+                    if img_bunit != (imp_bunit := self.hdu.header["BUNIT"]):
+                        logger.warning("Added mismatched BUNIT %s to %s.",
+                                       img_bunit, imp_bunit)
+                else:
+                    logger.info("No BUNIT found in added HDU.")
 
     @property
     def header(self):
diff --git a/scopesim/optics/image_plane_utils.py b/scopesim/optics/image_plane_utils.py
index d0742856..90315a0c 100644
--- a/scopesim/optics/image_plane_utils.py
+++ b/scopesim/optics/image_plane_utils.py
@@ -187,6 +187,7 @@ def create_wcs_from_points(points: np.ndarray,
     pixel_scale = abs(pixel_scale)
     extent = points.ptp(axis=0) / pixel_scale
     naxis = extent.round().astype(int)
+    ndims = len(naxis)
     offset = (extent - naxis) * pixel_scale  # compensate rounding
 
     # FIXME: Woule be nice to have D headers referenced at (1, 1), but that
@@ -214,10 +215,15 @@ def create_wcs_from_points(points: np.ndarray,
 
     # Cannot do `in "DX"` here because that would also match the empty string.
     linsuff = {"D", "X"}
-    lintype = ["LINEAR", "LINEAR"]
+    lintype = ndims * ["LINEAR"]
     # skytype = ["RA---TAN", "DEC--TAN"]  # ScopeSim can't handle the truth yet
     skytype = ["LINEAR", "LINEAR"]
-    ctype = lintype if wcs_suffix in linsuff else skytype
+    if wcs_suffix in linsuff:
+        ctype = lintype
+    else:
+        if ndims != 2:
+            raise ValueError("The sky must be two-dimensional!")
+        ctype = skytype
 
     if isinstance(points, u.Quantity):
         cunit = points.unit
@@ -227,10 +233,10 @@ def create_wcs_from_points(points: np.ndarray,
     if isinstance(pixel_scale, u.Quantity):
         pixel_scale = pixel_scale.value
 
-    new_wcs = WCS(key=wcs_suffix)
+    new_wcs = WCS(key=wcs_suffix, naxis=ndims)
     new_wcs.wcs.ctype = ctype
-    new_wcs.wcs.cunit = 2 * [cunit]
-    new_wcs.wcs.cdelt = np.array(2 * [pixel_scale])
+    new_wcs.wcs.cunit = ndims * [cunit]
+    new_wcs.wcs.cdelt = np.array(ndims * [pixel_scale])
     new_wcs.wcs.crval = crval
     new_wcs.wcs.crpix = crpix
 
@@ -424,40 +430,47 @@ def overlay_image(small_im, big_im, coords, mask=None, sub_pixel=False):
     if sub_pixel:
         raise NotImplementedError
 
+    if len(coords) != small_im.ndim:
+        coords = np.array([*coords, (big_im.shape[0] - 1) / 2])
+
     # FIXME: this would not be necessary if we used WCS instead of manual 2pix
     coords = np.ceil(np.asarray(coords).round(10)).astype(int)
-    y, x = coords.astype(int)[::-1] - np.array(small_im.shape[-2:]) // 2
+    idx = coords.astype(int)[::-1] - np.array(small_im.shape) // 2
 
     # Image ranges
-    x1, x2 = max(0, x), min(big_im.shape[-1], x + small_im.shape[-1])
-    y1, y2 = max(0, y), min(big_im.shape[-2], y + small_im.shape[-2])
+    idx1 = np.maximum(0, idx)
+    idx2 = np.minimum(big_im.shape, idx + small_im.shape)
 
     # Overlay ranges
-    x1o, x2o = max(0, -x), min(small_im.shape[-1], big_im.shape[-1] - x)
-    y1o, y2o = max(0, -y), min(small_im.shape[-2], big_im.shape[-2] - y)
+    idx1o = np.maximum(0, -idx)
+    idx2o = np.minimum(small_im.shape, big_im.shape - idx)
+
+    # Convert to tuples of slices for indexing
+    idx_slcs = tuple(slice(i1, i2) for i1, i2 in zip(idx1, idx2))
+    idx_slcso = tuple(slice(i1, i2) for i1, i2 in zip(idx1o, idx2o))
 
     # Exit if nothing to do
-    if y1 >= y2 or x1 >= x2 or y1o >= y2o or x1o >= x2o:
+    if (idx1 >= idx2).any() or (idx1o >= idx2o).any():
         return big_im
 
-    if small_im.ndim == 2 and big_im.ndim == 2:
-        small_im_3 = small_im[None, :, :]
-        big_im_3 = big_im[None, :, :]
-    elif small_im.ndim == 3 and big_im.ndim == 3:
-        small_im_3 = small_im
-        big_im_3 = big_im
-    else:
+    if mask is not None:
+        # The mask option was never used anywhere and got complicated.
+        # Keeping the previous code around just in case.
+        # mask = mask[None, idx_slcso[-2:]] * np.ones(small_im_3.shape[-3])
+        # mask = mask.astype(bool)
+        # big_im_3[idx_slcs][mask] = big_im_3[idx_slcs][mask] + \
+        #     small_im_3[idx_slcso][mask]
+        raise NotImplementedError()
+
+    if small_im.ndim == 2 != big_im.ndim == 2:
         raise ValueError(f"Dimensions mismatch between big_im and small_im: "
                          f"{big_im.ndim} : {small_im.ndim}")
 
-    if mask is None:
-        big_im_3[:, y1:y2, x1:x2] = big_im_3[:, y1:y2, x1:x2] + \
-                                    small_im_3[:, y1o:y2o, x1o:x2o]
-    else:
-        mask = mask[None, y1o:y2o, x1o:x2o] * np.ones(small_im_3.shape[-3])
-        mask = mask.astype(bool)
-        big_im_3[:, y1:y2, x1:x2][mask] = big_im_3[:, y1:y2, x1:x2][mask] + \
-                                          small_im_3[:, y1o:y2o, x1o:x2o][mask]
+    # Note: Could also do np.add(big_im, idx_slcs, small_im),
+    #       IF small_im[idx_slcso] == small_im (is that always?)
+    #       Both options modify the original as far as I can tell.
+
+    big_im[idx_slcs] = big_im[idx_slcs] + small_im[idx_slcso]
 
     return big_im
 
@@ -744,7 +757,7 @@ def add_imagehdu_to_imagehdu(image_hdu: fits.ImageHDU,
         canvas_wcs = wcs_suffix
     else:
         wcs_suffix = wcs_suffix or " "
-        canvas_wcs = WCS(canvas_hdu.header, key=wcs_suffix, naxis=2)
+        canvas_wcs = WCS(canvas_hdu.header, key=wcs_suffix)
 
     if isinstance(image_hdu.data, u.Quantity):
         image_hdu.data = image_hdu.data.value
@@ -766,11 +779,16 @@ def add_imagehdu_to_imagehdu(image_hdu: fits.ImageHDU,
                                 spline_order=spline_order,
                                 conserve_flux=conserve_flux)
 
-    img_center = np.array([[new_hdu.header["NAXIS1"],
-                            new_hdu.header["NAXIS2"]]])
+    img_center = np.array([[new_hdu.header[f"NAXIS{i+1}"]
+                           for i in range(new_hdu.header["NAXIS"])]])
     img_center = (img_center - 1) / 2
 
-    new_wcs = WCS(new_hdu.header, key=canvas_wcs.wcs.alt, naxis=2)
+    new_wcs = WCS(
+        new_hdu.header,
+        key=canvas_wcs.wcs.alt,
+        naxis=canvas_wcs.naxis,
+    )
+
     sky_center = new_wcs.wcs_pix2world(img_center, 0)
     if new_wcs.wcs.cunit[0] == "deg":
         sky_center = _fix_360(sky_center)
diff --git a/scopesim/optics/optical_train.py b/scopesim/optics/optical_train.py
index 91ed75a3..a6ff266a 100644
--- a/scopesim/optics/optical_train.py
+++ b/scopesim/optics/optical_train.py
@@ -283,8 +283,15 @@ class OpticalTrain:
                                    desc=" FOV effects", position=1):#, leave=False):
                     fov = effect.apply_to(fov)
 
-                fov.flatten()
-                self.image_planes[fov.image_plane_id].add(fov.hdu, wcs_suffix="D")
+                if self.cmds.get("!INST.flatten", True):
+                    fov.flatten()
+                    self.image_planes[fov.image_plane_id].add(fov.hdu, wcs_suffix="D")
+                else:  # cube output
+                    self.image_planes[fov.image_plane_id].add(fov.hdu, wcs_suffix="D")
+                    # HACK: to get sky WCS from FOV into image plane...
+                    self.image_planes[fov.image_plane_id].header.update(
+                        WCS(fov.hdu).to_header())
+
                 # ..todo: finish off the multiple image plane stuff
 
         # END OF MULTIPROCESSING
@@ -372,6 +379,10 @@ class OpticalTrain:
             fov_waveset = np.arange(wave_min.value, wave_max.value, dwave) * wave_unit
             fov_waveset = fov_waveset.to(u.um)
 
+            if (wave.to(u.um).min() > fov_waveset.max() or
+                    wave.to(u.um).max() < fov_waveset.min()):
+                logger.warning("Source and FOV wavesets disjoint.")
+
             # Interpolate into new data cube.
             # This is done layer by layer for memory reasons.
             new_data = np.zeros((fov_waveset.shape[0], data.shape[1], data.shape[2]),
@@ -630,6 +641,19 @@ class OpticalTrain:
             with p.indent(2):
                 p.pretty(self.effects)
 
+    def _repr_html_(self) -> str:
+        """For notebooks."""
+        html = (
+            "<div>"
+            f"<h4>{self.__class__.__name__} "
+            f"for {self.cmds['!OBS.instrument']} "
+            f"@ {self.cmds['!TEL.telescope']}</h4>"
+            f"<div><h5>Settings</h5>{self.cmds._repr_html_()}</div>"
+            f"<div><h5>Effects</h5>{self.effects._repr_html_()}</div>"
+            "</div>"
+        )
+        return html
+
     def __getitem__(self, item):
         return self.optics_manager[item]
 


Test patch:
diff --git a/scopesim/tests/mocks/basic_instrument/YAML_detector.yaml b/scopesim/tests/mocks/basic_instrument/YAML_detector.yaml
index c38ff61c..87da7b22 100644
--- a/scopesim/tests/mocks/basic_instrument/YAML_detector.yaml
+++ b/scopesim/tests/mocks/basic_instrument/YAML_detector.yaml
@@ -18,6 +18,7 @@ effects:
 - name: detector_window
   class: DetectorWindow
   description: Cut-out of the focal plane image with custom dimensions and coordinates
+  include: "!OBS.include_det_window"
   kwargs:
     image_plane_id: 0
     pixel_size: 0.01
@@ -27,6 +28,28 @@ effects:
     height: "!DET.height"
     units: pixel
 
+- name: detector_3d
+  class: DetectorList3D
+  description: Virtual 3D detector for simulating reduced ifu cubes
+  include: "!OBS.include_det_3d"
+  kwargs:
+    x_size_unit: pixel
+    y_size_unit: pixel
+    z_size_unit: pixel
+    image_plane_id: 0
+    pixel_scale: 0.5
+    array_dict:
+      id: [0]
+      x_cen: [0]
+      y_cen: [0]
+      z_cen: [0]
+      x_size: [15]
+      y_size: [25]
+      z_size: [5]
+      pixel_size: [0.01]
+      angle: [0]
+      gain: [1]
+
 - name : qe_curve
   description : Quantum efficiency curve (pseudo transmission curve)
   class : QuantumEfficiencyCurve
@@ -84,3 +107,12 @@ effects:
   class : ExtraFitsKeywords
   kwargs :
     filename: FITS_extra_keywords.yaml
+
+---
+
+name: config_overrides
+alias: OBS
+description: OBS config param overrides
+
+properties:
+  wavelen: 2.1
diff --git a/scopesim/tests/mocks/basic_instrument/YAML_mode_simple_ifu.yaml b/scopesim/tests/mocks/basic_instrument/YAML_mode_simple_ifu.yaml
new file mode 100644
index 00000000..163ff73f
--- /dev/null
+++ b/scopesim/tests/mocks/basic_instrument/YAML_mode_simple_ifu.yaml
@@ -0,0 +1,33 @@
+# IFU Spectroscopy mode specific Effects and properties
+name : mode_simple_ifu
+alias : INST
+
+properties :
+  decouple_detector_from_sky_headers: True
+
+effects :
+- name : flux_binning
+  description : Binning sampled flux for virtual 3D detector
+  class : FluxBinning3D
+
+---
+
+name: config_overrides
+alias: SIM
+description: SIM config param overrides
+
+properties:
+  spectral:
+    wave_min: 1.7
+    wave_mid: 2.1
+    wave_max: 2.5
+    spectral_bin_width: !!float 9E-4
+
+---
+
+name: config_overrides
+alias: INST
+description: INST config param overrides
+
+properties:
+  flatten: False
diff --git a/scopesim/tests/mocks/basic_instrument/default.yaml b/scopesim/tests/mocks/basic_instrument/default.yaml
index 1f270daa..fb4cec99 100644
--- a/scopesim/tests/mocks/basic_instrument/default.yaml
+++ b/scopesim/tests/mocks/basic_instrument/default.yaml
@@ -30,6 +30,8 @@ mode_yamls :
   properties :
     include_slit : False
     include_slicer: False
+    include_det_window: True
+    include_det_3d: False
     filter_name : "J"
 
 - name : spectroscopy
@@ -39,6 +41,8 @@ mode_yamls :
   properties :
     include_slit : True
     include_slicer: False
+    include_det_window: True
+    include_det_3d: False
     filter_name : "open"
   yamls :
     - YAML_mode_spectroscopy.yaml   # mode specific effects and properties
@@ -50,10 +54,25 @@ mode_yamls :
   properties:
     include_slit: False
     include_slicer: True
+    include_det_window: True
+    include_det_3d: False
     filter_name: "open"
   yamls:
     - YAML_mode_ifu.yaml   # mode specific effects and properties
 
+- name : simple_ifu
+  description: For testing the virtual 3D Detector with direct cube output
+  status: development
+  alias: OBS
+  properties:
+    include_slit: False
+    include_slicer: False
+    include_det_window: False
+    include_det_3d: True
+    filter_name: "open"
+  yamls:
+    - YAML_mode_simple_ifu.yaml   # mode specific effects and properties
+
 - name: mock_concept_mode
   description: Dummy mode to test concept status.
   status: concept
diff --git a/scopesim/tests/test_basic_instrument/test_basic_instrument.py b/scopesim/tests/test_basic_instrument/test_basic_instrument.py
index aa660b8f..f0847b2c 100644
--- a/scopesim/tests/test_basic_instrument/test_basic_instrument.py
+++ b/scopesim/tests/test_basic_instrument/test_basic_instrument.py
@@ -218,6 +218,35 @@ class TestObserveIfuMode:
             plt.show()
 
 
+@pytest.mark.slow
+@pytest.mark.usefixtures("protect_currsys", "patch_all_mock_paths")
+class TestObserveSimpleIfuMode:
+    def test_runs(self):
+        wave = np.arange(0.7, 2.5, 0.001)
+        spec = np.zeros_like(wave)
+        spec[25::50] += 100      # every 0.05µm, offset by 0.025µm
+        x = np.tile(np.arange(-4, 5, 2), 5)
+        y = np.repeat(np.arange(-4, 5, 2), 5)
+        src = sim.Source(lam=wave*u.um, spectra=spec,
+                         x=x, y=y, ref=[0]*len(x), weight=[1e-3]*len(x))
+
+        cmd = sim.UserCommands(
+            use_instrument="basic_instrument",
+            set_modes=["simple_ifu"],
+            ignore_effects=SWITCHOFF,
+        )
+        opt = sim.OpticalTrain(cmd)
+
+        opt.observe(src)
+        hdul = opt.readout()[0]
+
+        imp_im = opt.image_planes[0].data
+        det_im = hdul[1].data
+
+        assert imp_im.ndim == 3
+        assert det_im.ndim == 3
+
+
 @pytest.mark.usefixtures("protect_currsys", "patch_all_mock_paths")
 class TestFitsHeader:
     def test_source_keywords_in_header(self):
@@ -230,7 +259,7 @@ class TestFitsHeader:
         hdr = hdul[0].header
 
         assert hdr["SIM SRC0 object"] == 'star'
-        assert hdr["SIM EFF14 class"] == 'SourceDescriptionFitsKeywords'
+        assert hdr["SIM EFF15 class"] == 'SourceDescriptionFitsKeywords'
         assert hdr["SIM CONFIG OBS filter_name"] == 'J'
         assert hdr["ESO ATM SEEING"] == opt.cmds["!OBS.psf_fwhm"]
 
diff --git a/scopesim/tests/tests_effects/test_DetectorList.py b/scopesim/tests/tests_effects/test_DetectorList.py
index bd9f6faf..a9570167 100644
--- a/scopesim/tests/tests_effects/test_DetectorList.py
+++ b/scopesim/tests/tests_effects/test_DetectorList.py
@@ -2,8 +2,9 @@ import pytest
 from unittest.mock import patch
 
 from astropy.table import Table
+from astropy import units as u
 
-from scopesim.effects import DetectorList, DetectorWindow, ApertureMask
+from scopesim.effects import DetectorList, DetectorWindow, DetectorList3D
 
 
 @pytest.fixture(scope="class")
@@ -154,24 +155,6 @@ class TestDetecotrHeaders:
             assert hdr["NAXIS2"] == 32
 
 
-@pytest.mark.usefixtures("patch_mock_path_micado")
-class TestFovGrid:
-    def test_returns_aperture_mask_object(self):
-        det_list = DetectorList(filename="FPA_array_layout.dat",
-                                image_plane_id=0)
-        apm = det_list.fov_grid(pixel_scale=0.004)
-        assert isinstance(apm, ApertureMask)
-
-    def test_aperture_mask_header_covers_all_of_detector_header(self):
-        det_list = DetectorList(filename="FPA_array_layout.dat",
-                                image_plane_id=0)
-        apm = det_list.fov_grid(pixel_scale=0.004)
-        apm_hdr = apm.header
-        det_hdr = det_list.image_plane_header
-        assert apm_hdr["NAXIS1"] == det_hdr["NAXIS1"]
-        assert apm_hdr["NAXIS2"] == det_hdr["NAXIS2"]
-
-
 class TestDetectorWindowInit:
     def test_throws_error_when_initialises_with_nothing(self):
         with pytest.raises(TypeError):
@@ -213,3 +196,44 @@ class TestDetectorWindowInit:
                                  width="!DET.width", units="pixel")
             assert det.image_plane_header["NAXIS1"] == 42
             assert det.detector_headers()[0]["NAXIS1"] == 42
+
+
+@pytest.fixture(name="det_list")
+def basic_detectorlist3d():
+    det_list = DetectorList3D(array_dict={
+        "id": [0],
+        "x_cen": [0],
+        "y_cen": [0],
+        "z_cen": [0],
+        "x_size": [15],
+        "y_size": [25],
+        "z_size": [5],
+        "pixel_size": [0.01],
+        "angle": [0],
+        "gain": [1]},
+        x_size_unit="pixel",
+        y_size_unit="pixel",
+        z_size_unit="pixel",
+        image_plane_id=0,
+        pixel_scale=0.5,
+    )
+    return det_list
+
+
+@pytest.mark.usefixtures("patch_all_mock_paths")
+class TestDetectorList3D:
+    def test_image_plane_header_is_3d(self, det_list):
+        impl_hdr = det_list.image_plane_header
+
+        assert impl_hdr["NAXIS"] == 3
+        assert impl_hdr["NAXIS1"] == 15
+        assert impl_hdr["NAXIS2"] == 25
+        assert impl_hdr["NAXIS3"] == 5
+
+    def test_pixel_scale_arcsec(self, det_list):
+        equiv = det_list.pixel_scale_arcsec
+        assert (1 * u.pixel).to(u.arcsec, equiv) == .5 * u.arcsec
+
+    def test_pixel_scale_mm(self, det_list):
+        equiv = det_list.pixel_scale_mm
+        assert (1 * u.pixel).to(u.mm, equiv) == .01 * u.mm
diff --git a/scopesim/tests/tests_optics/test_FieldOfView.py b/scopesim/tests/tests_optics/test_FieldOfView.py
index 380325ab..77dd9fad 100644
--- a/scopesim/tests/tests_optics/test_FieldOfView.py
+++ b/scopesim/tests/tests_optics/test_FieldOfView.py
@@ -464,6 +464,7 @@ class TestMakeSpectrum:
 
         assert in_sum == approx(out_sum)
 
+    @pytest.mark.xfail(reason="Fails after fixing spectral off-by-one for cubes.")
     def test_makes_spectrum_from_all_types_of_source_object(self):
         src_table = so._table_source()
         src_image = so._image_source(dx=-4, dy=-4)
diff --git a/scopesim/tests/tests_optics/test_fov_manager_utils.py b/scopesim/tests/tests_optics/test_fov_manager_utils.py
index 1384ebf4..5673d251 100644
--- a/scopesim/tests/tests_optics/test_fov_manager_utils.py
+++ b/scopesim/tests/tests_optics/test_fov_manager_utils.py
@@ -141,6 +141,8 @@ class TestGetImagingHeaders:
         area_sum = np.sum([hdr["NAXIS1"] * hdr["NAXIS2"] for hdr in hdrs])
         assert area_sum == 228 * 328
 
+    @pytest.mark.skip(reason="this now fails because it's testing an unused "
+                      "function that uses an otherwise unused method...")
     def test_returns_set_of_headers_from_detector_list_effect(self):
         # det = eo._full_detector_list()
         det = eo._detector_list()