Repository: AstarVienna/ScopeSim

Issue:
Allow `Source` to be created directly from existing `SourceField`
The current API of the `Source` class has a rather convoluted constructor, to keep it very flexible for many different use cases. This is useful in general, and before the introduction of the `SourceField` (sub)class(es), there were not many other options. Direct supply of tables and image HDUs plus spectra was always possible, and is now redirected to do little more than create the appropriate field class instance.

It would be useful to have an equally easy way to simply supply a `field` argument to `Source()`, using a preexisting `SourceField` instance to initialize the `Source` object. This should fit neatly in the existing constructor logic for the class.

Discussion:

Solution patch:
diff --git a/.github/workflows/publish_pypi.yml b/.github/workflows/publish_pypi.yml
index 2dad0848..0370131b 100644
--- a/.github/workflows/publish_pypi.yml
+++ b/.github/workflows/publish_pypi.yml
@@ -67,7 +67,7 @@ jobs:
         run: echo "## Publishing $(poetry version) to PyPI" >> $GITHUB_STEP_SUMMARY
 
       - name: Retrieve distribution files
-        uses: actions/download-artifact@v4.2.1
+        uses: actions/download-artifact@v4.3.0
         with:
           name: distribution
           path: dist/
diff --git a/scopesim/effects/detector_list.py b/scopesim/effects/detector_list.py
index 1052852d..03c23099 100644
--- a/scopesim/effects/detector_list.py
+++ b/scopesim/effects/detector_list.py
@@ -20,7 +20,7 @@ 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,
+from ..utils import (from_currsys, close_loop, figure_factory, array_minmax,
                      quantity_from_table, unit_from_table, get_logger)
 
 logger = get_logger(__name__)
@@ -164,7 +164,7 @@ class DetectorList(Effect):
                    values=(tuple(zip(xy_sky.min(axis=0),
                                      xy_sky.max(axis=0)))))
 
-        lims = np.array((xy_mm.min(axis=0), xy_mm.max(axis=0)))
+        lims = array_minmax(xy_mm)
         keys = ["xd_min", "xd_max", "yd_min", "yd_max"]
         obj.detector_limits = dict(zip(keys, lims.T.flatten()))
 
diff --git a/scopesim/effects/metis_wcu/fpmask.py b/scopesim/effects/metis_wcu/fpmask.py
index d559ae29..ec40efd4 100644
--- a/scopesim/effects/metis_wcu/fpmask.py
+++ b/scopesim/effects/metis_wcu/fpmask.py
@@ -1,45 +1,54 @@
-"""A class for the METIS WCU focal-plane mask"""
+# -*- coding: utf-8 -*-
+"""A class for the METIS WCU focal-plane mask."""
 
 from pathlib import Path
+
 import numpy as np
 from astropy.io import fits
 from astropy import units as u
+
 from ..data_container import DataContainer
 from ...utils import find_file, from_currsys, get_logger, figure_factory
 from ...optics.image_plane_utils import sub_pixel_fractions
 
 logger = get_logger(__name__)
 
+
 class FPMask:
-    """Focal-plane mask for the METIS WCU
+    """
+    Focal-plane mask for the METIS WCU.
 
     Parameters
     ----------
     See :class:`DataContainer` for input parameters
 
     """
-    hdr = {"BG_SRC": True,
-           "BG_SURF": "WCU focal plane mask",   # TODO more specific?
-           "CTYPE1": "LINEAR",
-           "CTYPE2": "LINEAR",
-           "CRPIX1": 1024.5,
-           "CRPIX2": 1024.5,
-           "CRVAL1": 0.,
-           "CRVAL2": 0.,
-           "CUNIT1": "arcsec",
-           "CUNIT2": "arcsec",
-           "CDELT1": 0.00547,
-           "CDELT2": 0.00547,
-           "BUNIT": "PHOTLAM arcsec-2",
-           "SOLIDANG": "arcsec-2"}
-
-    def __init__(self,
-                 maskname: Path | str | None = None,
-                 fpmask_filename_format: str | None = None,
-                 angle: float = 0,
-                 shift: tuple = (0, 0),
-                 **kwargs
-                 ):
+
+    hdr = {
+        "BG_SRC": True,
+        "BG_SURF": "WCU focal plane mask",  # TODO: more specific?
+        "CTYPE1": "LINEAR",
+        "CTYPE2": "LINEAR",
+        "CRPIX1": 1024.5,
+        "CRPIX2": 1024.5,
+        "CRVAL1": 0.,
+        "CRVAL2": 0.,
+        "CUNIT1": "arcsec",
+        "CUNIT2": "arcsec",
+        "CDELT1": 0.00547,
+        "CDELT2": 0.00547,
+        "BUNIT": "PHOTLAM arcsec-2",
+        "SOLIDANG": "arcsec-2",
+    }
+
+    def __init__(
+        self,
+        maskname: Path | str | None = None,
+        fpmask_filename_format: str | None = None,
+        angle: float = 0,
+        shift: tuple[float, float] = (0, 0),
+        **kwargs
+    ):
         logger.debug("Initialising FPMask with %s", maskname)
         self.name = maskname
         self.angle = angle
@@ -59,16 +68,17 @@ class FPMask:
                 self.filename = maskname
 
             self.data_container = DataContainer(filename=self.filename, **kwargs)
-            self.pixarea = (self.hdr['CDELT1'] * u.Unit(self.hdr['CUNIT1'])
-                            * self.hdr['CDELT2'] * u.Unit(self.hdr['CUNIT2']))
+            self.pixarea = (self.hdr["CDELT1"] * u.Unit(self.hdr["CUNIT1"])
+                            * self.hdr["CDELT2"] * u.Unit(self.hdr["CUNIT2"]))
             self.make_hdus(header=self.hdr)
 
-
     def make_hdus(self, header):
-        """Create an hdu for the holes in fpmask
+        """
+        Create an hdu for the holes in fpmask.
 
-        The holes are assumed to be unresolved. They therefore cover one pixel and have
-        a value corresponding to the actual solid angle covered by the hole.
+        The holes are assumed to be unresolved. They therefore cover one pixel
+        and have a value corresponding to the actual solid angle covered by
+        the hole.
         """
         holehdu = fits.ImageHDU()
         holehdu.header.update(header)
@@ -80,9 +90,9 @@ class FPMask:
 
         # Hole locations
         tab = self.data_container.table
-        xhole = tab['x'].data
-        yhole = tab['y'].data
-        diam = tab['diam'].data
+        xhole = tab["x"].data
+        yhole = tab["y"].data
+        diam = tab["diam"].data
 
         if self.angle != 0:
             rangle = np.deg2rad(self.angle)
@@ -94,33 +104,34 @@ class FPMask:
         xhole += self.shift[0]
         yhole += self.shift[1]
 
-        xpix = (xhole - header['CRVAL1']) / header['CDELT1'] + header['CRPIX1'] - 1
-        ypix = (yhole - header['CRVAL2']) / header['CDELT2'] + header['CRPIX2'] - 1
+        xpix = (xhole - header["CRVAL1"]) / header["CDELT1"] + header["CRPIX1"] - 1
+        ypix = (yhole - header["CRVAL2"]) / header["CDELT2"] + header["CRPIX2"] - 1
         in_field = (xpix > 0) * (xpix < 2047) * (ypix > 0) * (ypix < 2047)
 
-
         for x, y, d in zip(xpix[in_field], ypix[in_field], diam[in_field]):
             holearea = (d/2)**2 * np.pi
             xint, yint, fracs = sub_pixel_fractions(x, y)
             holehdu.data[yint, xint] = np.array(fracs) * holearea
             opaquehdu.data[yint, xint] = 0
+
         self.xpix = xpix
         self.ypix = ypix
         self.holehdu = holehdu
         self.opaquehdu = opaquehdu
 
-
     def plot(self):
-        """Plot the location of the holes"""
+        """Plot the location of the holes."""
         _, axes = figure_factory()
-        axes.plot(self.xpix, self.ypix, 'o')
+        axes.plot(self.xpix, self.ypix, "o")
         axes.set_xlim(0, 2048)
         axes.set_ylim(0, 2048)
-        axes.set_aspect('equal')
+        axes.set_aspect("equal")
 
         return axes
 
     def __str__(self) -> str:
-        return f"""{self.__class__.__name__}: "{self.name}"
-    Angle:        {self.angle} deg
-    Shift:        {self.shift} arcsec"""
+        """Return str(self)."""
+        msg = (f"{self.__class__.__name__}: \"{self.name}\"\n"
+               f"Angle:        {self.angle} deg\n"
+               f"Shift:        {self.shift} arcsec\n")
+        return msg
diff --git a/scopesim/effects/metis_wcu/metis_wcu.py b/scopesim/effects/metis_wcu/metis_wcu.py
index 105dd640..7029c151 100644
--- a/scopesim/effects/metis_wcu/metis_wcu.py
+++ b/scopesim/effects/metis_wcu/metis_wcu.py
@@ -1,4 +1,5 @@
-"""Classes for the METIS Warm Calibration Unit"""
+# -*- coding: utf-8 -*-
+"""Classes for the METIS Warm Calibration Unit."""
 
 from typing import ClassVar
 
@@ -14,16 +15,18 @@ import yaml
 
 from .fpmask import FPMask
 from ..ter_curves import TERCurve, FilterCurve
-from ...utils import get_logger, seq, find_file,\
-    from_currsys
+from ...utils import get_logger, seq, find_file, from_currsys
 from ...source.source import Source
+from ...source.source_fields import BackgroundSourceField
 from ...optics.surface import SpectralSurface
 
 
 logger = get_logger(__name__)
 
+
 class WCUSource(TERCurve):
-    """Warm Calibration Unit Source
+    """
+    Warm Calibration Unit Source.
 
     This class returns a TERCurve that describes the output emission
     from the integrating sphere of the METIS WCU. The calculations include
@@ -37,29 +40,36 @@ class WCUSource(TERCurve):
     The default configuration file for the METIS WCU is provided in the irdb.
     To modify the configuration, copy the file <irdb>/METIS/metis_wcu_config.yaml
     to the working directory and edit it there. To use the file do e.g.
-    >>> cmds = sim.UserCommands(use_instrument="METIS", set_modes=["wcu_img_lm"]
+
+    >>> cmds = sim.UserCommands(use_instrument="METIS", set_modes=["wcu_img_lm"])
     >>> cmds["!WCU.config_file"] = "my_wcu_config.yaml"
     >>> metis = sim.OpticalTrain(cmds)
 
     Changing important parameters
     -----------------------------
-    The temperatures of the black-body source, the integrating sphere as well as the
-    ambient temperature of the WCU (which is the temperature of the focal-plane mask)
-    can be set by calling
-    >>> metis['wcu_source'].set_temperature(bb_temp=1200*u.K, is_temp=320*u.K,
+    The temperatures of the black-body source, the integrating sphere as well
+    as the ambient temperature of the WCU (which is the temperature of the
+    focal-plane mask) can be set by calling
+
+    >>> metis["wcu_source"].set_temperature(bb_temp=1200*u.K, is_temp=320*u.K,
                                             wcu_temp=295*u.K)
 
-    The flux-controlling masks are currently implemented as a float (in [0, 1]) that
-    gives the fraction of flux transmitted. It can be change with
-    >>> metis['wcu_source'].set_bb_aperture(0.8)
+    The flux-controlling masks are currently implemented as a float (in [0, 1])
+    that gives the fraction of flux transmitted. It can be changed with
 
-    The focal-plane mask and/or its orientation and position can be changed by calling
-    >>> metis['wcu_source'].set_fpmask("pinhole", angle=10, shift=(0.1, 0))
+    >>> metis["wcu_source"].set_bb_aperture(0.8)
+
+    The focal-plane mask and/or its orientation and position can be changed by
+    calling
+
+    >>> metis["wcu_source"].set_fpmask("pinhole", angle=10, shift=(0.1, 0))
 
     To use the lasers instead of the black-body source, do
-    >>> metis['wcu_source'].set_lamp("laser")
-    Which laser is seen depends on the wavelength of the observation. Note that the
-    tunable laser currently cannot be tuned.
+
+    >>> metis["wcu_source"].set_lamp("laser")
+
+    Which laser is seen depends on the wavelength of the observation. Note that
+    the tunable laser currently cannot be tuned.
     """
 
     z_order: ClassVar[tuple[int, ...]] = (113, 513)
@@ -71,8 +81,8 @@ class WCUSource(TERCurve):
             "position": 0,  # position in surface table
         }
         self.meta.update(params)
-        if 'config_file' in self.meta:
-            config_file = from_currsys(self.meta['config_file'], self.cmds)
+        if "config_file" in self.meta:
+            config_file = from_currsys(self.meta["config_file"], self.cmds)
             with open(find_file(config_file), encoding="utf-8") as fd:
                 config = yaml.safe_load(fd)
                 self.meta.update(config)
@@ -82,15 +92,15 @@ class WCUSource(TERCurve):
         # Check on the presence of one vital parameter
         if "rho_tube" not in self.meta:
             raise ValueError(
-                "Parameters not present: please provide config file or parameter values"
-                "for WCUSource")
+                "Parameters not present: please provide config file or "
+                "parameter values for WCUSource")
 
-        self.set_fpmask(maskname=self.meta['current_fpmask'],
+        self.set_fpmask(maskname=self.meta["current_fpmask"],
                         angle=self.meta["fpmask_angle"],
                         shift=self.meta["fpmask_shift"])
         self._background_source = None
 
-        self.bb_aperture = self.meta['bb_aperture']
+        self.bb_aperture = self.meta["bb_aperture"]
 
         # Define surfaces (the IS surface is currently taken from parent class)
         # we add another one here (TODO: make this nicer)
@@ -101,9 +111,9 @@ class WCUSource(TERCurve):
         self.get_wavelength()
         self.bb_scale = 1 * u.ph / (u.s * u.m**2 * u.arcsec**2 * u.um)
         self.bb_to_is = self.bb_to_is_throughput()
-        self.rho_tube = get_reflectivity(self.meta['rho_tube'])
-        self.rho_is = get_reflectivity(self.meta['rho_is'])
-        self.emiss_mask = self.meta['emiss_mask']
+        self.rho_tube = get_reflectivity(self.meta["rho_tube"])
+        self.rho_is = get_reflectivity(self.meta["rho_is"])
+        self.emiss_mask = self.meta["emiss_mask"]
 
         # Compute the emission components
         self.compute_lamp_emission()
@@ -119,69 +129,101 @@ class WCUSource(TERCurve):
 
     @property
     def background_source(self):
-        """Define a source field for the FP mask"""
-
+        """Define a source field for the FP mask."""
         bb_flux = self.emission
         mask_flux = self.mask_emission
         self._background_source = []
 
-        self._background_source.append(Source(image_hdu=self.fpmask.holehdu,
-                                              spectra=bb_flux))
-        if self.fpmask.opaquehdu is not None:   # not the open mask
-            self._background_source.append(Source(image_hdu=self.fpmask.opaquehdu,
-                                                  spectra=mask_flux))
+        # TODO: Turn this into one Source object with multiple fields...
+        # TODO: Not sure how this is generally, if this check is really needed.
+        hole_hdu = self.fpmask.holehdu
+        if hole_hdu.data is None:
+            hole_spec = {0: bb_flux}
+            hole_hdu.header.update({"SPEC_REF": 0})
+            hole_fld = BackgroundSourceField(field=None, spectra=hole_spec,
+                                             header=hole_hdu.header)
+            hole_src = Source(field=hole_fld)
+        else:
+            hole_src = Source(image_hdu=self.fpmask.holehdu, spectra=bb_flux)
+
+        self._background_source.append(hole_src)
+
+        opaque_hdu = self.fpmask.opaquehdu
+        if opaque_hdu is not None:   # not the open mask
+            # TODO: same as above...
+            if opaque_hdu.data is None:
+                opaque_spec = {0: mask_flux}
+                opaque_hdu.header.update({"SPEC_REF": 0})
+                opaque_fld = BackgroundSourceField(field=None,
+                                                   spectra=opaque_spec,
+                                                   header=opaque_hdu.header)
+                opaque_src = Source(field=opaque_fld)
+            else:
+                opaque_src = Source(image_hdu=self.fpmask.opaquehdu,
+                                    spectra=mask_flux)
+            self._background_source.append(opaque_src)
+
         return self._background_source
 
-    def set_lamp(self, lamp='bb'):
-        """Change the WCU lamp
+    def set_lamp(self, lamp="bb"):
+        """
+        Change the WCU lamp.
 
         The list of available lamps can be obtained with
-        >>> metis['wcu_source'].lamps
+
+        >>> metis["wcu_source"].lamps
+
         The lamp currently in use is
-        >>> metis['wcu_source'].current_lamp
+
+        >>> metis["wcu_source"].current_lamp
+
         It is changed with, e.g.,
-        >>> metis['wcu_source'].set_lamp('laser')
+
+        >>> metis["wcu_source"].set_lamp("laser")
         """
         if lamp not in self.lamps:
             raise ValueError(f"'lamp' needs to be one of {self.lamps}")
 
-        self.meta['current_lamp'] = lamp
+        self.meta["current_lamp"] = lamp
 
         self.compute_lamp_emission()
         self.compute_fp_emission()
 
-
     def get_wavelength(self):
-        """Try to set the appropriate wavelength vector for the mode and filter"""
-        if 'wcu_lms' in self.cmds['!OBS.modes']:     ## Need to provide for wcu_lms_extended
-            lamc = self.cmds['!OBS.wavelen']
-            dlam = self.cmds['!SIM.spectral.spectral_bin_width']
+        """Try to set the appropriate wavelength vector for mode and filter."""
+        # Need to provide for wcu_lms_extended
+        if "wcu_lms" in self.cmds["!OBS.modes"]:
+            lamc = self.cmds["!OBS.wavelen"]
+            dlam = self.cmds["!SIM.spectral.spectral_bin_width"]
             lam = seq(lamc - 3000 * dlam, lamc + 3000 * dlam, dlam) * u.um
         else:
-            filter_name = self.cmds['!OBS.filter_name']
-            filename_format = self.cmds['!INST.filter_file_format']
+            filter_name = self.cmds["!OBS.filter_name"]
+            filename_format = self.cmds["!INST.filter_file_format"]
             tempfilter = FilterCurve(filter_name=filter_name,
                                      filename_format=filename_format)
             lammin, lammax = tempfilter.throughput.waverange << u.um
-            dlam = self.cmds['!SIM.spectral.spectral_bin_width']
+            dlam = self.cmds["!SIM.spectral.spectral_bin_width"]
             lam = seq(lammin.value, lammax.value, dlam) * u.um
         self.wavelength = lam
 
 
     @property
     def lamps(self):
-        """List of available lamps"""
-        return self.meta['lamps']
+        """Return list of available lamps."""
+        return self.meta["lamps"]
 
     @property
     def current_lamp(self):
-        """Name of the lamp currently in use"""
-        return self.meta['current_lamp']
-
-    def set_temperature(self, bb_temp: [float | u.Quantity]=None,
-                        is_temp: [float | u.Quantity]=None,
-                        wcu_temp: [float | u.Quantity]=None):
-        """Change the black-body temperature
+        """Return name of the lamp currently in use."""
+        return self.meta["current_lamp"]
+
+    def set_temperature(
+        self, bb_temp: float | u.Quantity | None = None,
+        is_temp: float | u.Quantity | None = None,
+        wcu_temp: float | u.Quantity | None = None,
+    ) -> None:
+        """
+        Change the black-body temperature.
 
         Parameters
         ----------
@@ -217,8 +259,9 @@ class WCUSource(TERCurve):
         self.compute_lamp_emission()
         self.compute_fp_emission()
 
-    def set_bb_aperture(self, value: float):
-        """Change the flux-controlling aperture for the black-body source
+    def set_bb_aperture(self, value: float) -> None:
+        """
+        Change the flux-controlling aperture for the black-body source.
 
         Currently, this accepts a float value between 0 and 1 that gives the
         fraction of black-body emission that exits the integrating sphere
@@ -231,12 +274,18 @@ class WCUSource(TERCurve):
         self.bb_aperture = value
         self.compute_lamp_emission()
 
-    def set_fpmask(self, maskname: str = None, angle: float = None, shift: tuple = None):
-        """Change the focal-plane mask
+    def set_fpmask(
+        self,
+        maskname: str | None = None,
+        angle: float | None = None,
+        shift: tuple[float, float] | None = None,
+    ) -> None:
+        """
+        Change the focal-plane mask.
 
-        If `maskname` is not given, the currently inserted mask is rotated to `angle`
-        and shifted to `shift`. If `maskname` is given, angle and shift are reset to
-        zero, unless explicitely specified.
+        If `maskname` is not given, the currently inserted mask is rotated to
+        `angle` and shifted to `shift`. If `maskname` is given, angle and shift
+        are reset to zero, unless explicitely specified.
 
         See also :class:`FPMask`.
 
@@ -259,45 +308,46 @@ class WCUSource(TERCurve):
         if shift is not None:
             self.meta["fpmask_shift"] = shift
 
-        self.fpmask = FPMask(maskname=self.meta["current_fpmask"],
-                             fpmask_filename_format=self.meta['fpmask_filename_format'],
-                             angle=self.meta["fpmask_angle"],
-                             shift=self.meta["fpmask_shift"])
-
+        self.fpmask = FPMask(
+            maskname=self.meta["current_fpmask"],
+            fpmask_filename_format=self.meta["fpmask_filename_format"],
+            angle=self.meta["fpmask_angle"],
+            shift=self.meta["fpmask_shift"],
+        )
 
     @property
     def current_fpmask(self):
-        """Name of the mask currently in use"""
-        return self.meta['current_fpmask']
-
+        """Return name of the mask currently in use."""
+        return self.meta["current_fpmask"]
 
     def bb_to_is_throughput(self):
-        """Load throughput and return interpolation function
+        """
+        Load throughput and return interpolation function.
 
         This loads a lookup table for the transmission from the black-body
-        source to the entrance port of the integrating sphere. It returns
-        an interpolation function that can be evaluated on a reflectivity for the tube.
+        source to the entrance port of the integrating sphere. It returns an
+        interpolation function that can be evaluated on a reflectivity for the
+        tube.
 
         The model used is the general model, currently with no gaps.
         """
-
-        path = find_file(self.meta['bb_to_is'])
+        path = find_file(self.meta["bb_to_is"])
         if path is None:
             return lambda x: 1
         with fits.open(path) as hdul:
-            rho_tube = hdul[1].data['rho_tube']
-            throughput = hdul[1].data['t_gen_no_gap']    # TODO: update
+            rho_tube = hdul[1].data["rho_tube"]
+            throughput = hdul[1].data["t_gen_no_gap"]  # TODO: update
         return interp1d(rho_tube, throughput)
 
     def compute_lamp_emission(self):
-        """Compute the emission at the exit of the integrating sphere"""
+        """Compute the emission at the exit of the integrating sphere."""
         self.d_is = self.meta["diam_is"] << u.mm
         self.d_is_in = self.meta["diam_is_in"] << u.mm
         self.d_is_out = self.meta["diam_is_out"] << u.mm
         self.bb_temp = self.meta["bb_temp"] << u.K
         self.is_temp = self.meta["is_temp"] << u.K
         self.wcu_temp = self.meta["wcu_temp"] << u.K
-        self.emiss_bb = self.meta["emiss_bb"]      # <<<<<< that needs to be a function
+        self.emiss_bb = self.meta["emiss_bb"]  # that needs to be a function
 
         lam = self.wavelength
 
@@ -306,10 +356,12 @@ class WCUSource(TERCurve):
         # continuum black-body source
         if self.current_lamp == "bb":
             self.is_lamp = BlackBody(self.bb_temp, scale=self.bb_scale)
-            self.flux_lamp = (self.emiss_bb * self.is_lamp(lam) * self.bb_aperture
-                            * (np.pi * self.d_is_in**2 / 4) * (np.pi * u.sr))
+            self.flux_lamp = (self.emiss_bb * self.is_lamp(lam)
+                              * self.bb_aperture
+                              * (np.pi * self.d_is_in**2 / 4) * (np.pi * u.sr))
             self.flux_lamp *= self.bb_to_is(self.rho_tube(lam))
-            self.intens_lamp = self.flux_lamp / (np.pi * self.d_is**2) * mult_is / (np.pi * u.sr)
+            self.intens_lamp = (self.flux_lamp / (np.pi * self.d_is**2)
+                                * mult_is / (np.pi * u.sr))
         elif self.current_lamp == "laser":
             self.intens_lamp = self._laser_intensity()
         elif self.current_lamp == "none":
@@ -319,7 +371,7 @@ class WCUSource(TERCurve):
 
         # background emission from integrating sphere
         self.is_bg = BlackBody(self.is_temp, scale=self.bb_scale)
-        self.intens_bg  = self.rho_is(lam) * self.is_bg(lam)
+        self.intens_bg = self.rho_is(lam) * self.is_bg(lam)
 
         self.intensity = self.intens_lamp + self.intens_bg
 
@@ -334,7 +386,8 @@ class WCUSource(TERCurve):
         self.surface.meta.update(tbl.meta)
 
     def _laser_intensity(self):
-        """Compute the intensity for the single-line lasers
+        """
+        Compute the intensity for the single-line lasers.
 
         The function computes both lasers at once. This is possible because the
         lines are so far apart that there is no band that sees them both.
@@ -342,7 +395,7 @@ class WCUSource(TERCurve):
         lam = self.wavelength
         dlam = lam[1] - lam[0]
 
-        # Laser 1 (L band)     ## TODO move to yaml
+        # Laser 1 (L band)  # TODO move to yaml
         lamc_l = 3.39 * u.um
         power_l = 5e-3 * u.W / (c.c * c.h / lamc_l) * u.ph
         # Laser 2 (tunable), power divided among multiple lines
@@ -360,13 +413,13 @@ class WCUSource(TERCurve):
         line_m = Gaussian1D(amplitude=amp, mean=lamc_m, stddev=sigma)
         list_t = [Gaussian1D(amplitude=amp, mean=ll, stddev=sigma) for ll in lam_t]
         line_t = sum(list_t[1:], start=list_t[0])
-        flux = (power_l * line_l(lam) + power_m * line_m(lam) + power_t * line_t(lam)) \
-            / (np.pi * self.d_is_in**2 / 4)
+        flux = ((power_l * line_l(lam) + power_m * line_m(lam) + power_t
+                 * line_t(lam)) / (np.pi * self.d_is_in**2 / 4))
         intens = flux / (np.pi * u.sr)
         return intens.to(self.bb_scale)
 
     def compute_fp_emission(self):
-        """Compute the emission spectrum from the opaque part of the focal-plane mask"""
+        """Compute emission spectrum from opaque part of focal-plane mask."""
         self.wcu_temp = self.meta["wcu_temp"] << u.K
 
         lam = self.wavelength
@@ -386,9 +439,9 @@ class WCUSource(TERCurve):
         self.mask_surf.table = tbl
         self.mask_surf.meta.update(tbl.meta)
 
-
     def is_multiplication(self, wavelength, nport=2):
-        """Intensity multiplication factor for the integrating sphere
+        """
+        Intensity multiplication factor for the integrating sphere.
 
         Parameters
         ----------
@@ -402,7 +455,7 @@ class WCUSource(TERCurve):
         """
         rho = self.rho_is(wavelength)
         area_in = (self.d_is_in**2
-                   * ( 1 - np.cos(np.arcsin(self.d_is_in / self.d_is))))
+                   * (1 - np.cos(np.arcsin(self.d_is_in / self.d_is))))
         area_out = (self.d_is_out**2
                     * (1 - np.cos(np.arcsin(self.d_is_out / self.d_is))))
         area_is = np.pi * self.d_is**2
@@ -414,25 +467,29 @@ class WCUSource(TERCurve):
         return rho / (1 - rho * (1 - frac))
 
     def __str__(self) -> str:
-        return f"""{self.__class__.__name__}: "{self.display_name}"
-    Current lamp:            {self.current_lamp}     {self.lamps}
-    BlackBody temperature:   {self.meta['bb_temp'] << u.K}
-    Integrating sphere temp: {self.meta['is_temp'] << u.K}
-    WCU temperature:         {self.meta['wcu_temp'] << u.K}
-    BlackBody aperture:      {self.bb_aperture}
-    Focal-plane mask:        {self.meta['current_fpmask']}  {self.meta['fpmasks']}"""
-
+        """Return str(self)."""
+        msg = (
+            f"{self.__class__.__name__}: \"{self.display_name}\"\n"
+            f"Current lamp:            {self.current_lamp}     {self.lamps}\n"
+            f"BlackBody temperature:   {self.meta['bb_temp'] << u.K}\n"
+            f"Integrating sphere temp: {self.meta['is_temp'] << u.K}\n"
+            f"WCU temperature:         {self.meta['wcu_temp'] << u.K}\n"
+            f"BlackBody aperture:      {self.bb_aperture}\n"
+            f"Focal-plane mask:        {self.meta['current_fpmask']}  "
+            f"{self.meta['fpmasks']}\n"
+        )
+        return msg
 
 
 def get_reflectivity(file_or_number):
     """
-    Get a reflectivity from either a file or a number
+    Get a reflectivity from either a file or a number.
 
     This returns a function of wavelength.
     Does not work with Quantity, assumes wavelength in um.
     """
     if isinstance(file_or_number, str):
         tbl = ioascii.read(find_file(file_or_number))
-        return interp1d(tbl['wavelength'], tbl['reflection'])
+        return interp1d(tbl["wavelength"], tbl["reflection"])
 
     return lambda x: file_or_number * np.ones_like(x.value)
diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py
index fdb14d26..bcd681bd 100644
--- a/scopesim/effects/ter_curves.py
+++ b/scopesim/effects/ter_curves.py
@@ -18,7 +18,8 @@ from .ter_curves_utils import (add_edge_zeros, combine_two_spectra,
 from ..optics.fov_volume_list import FovVolumeList
 from ..optics.surface import SpectralSurface
 from ..source.source import Source
-from ..source.source_fields import CubeSourceField, SpectrumSourceField
+from ..source.source_fields import (CubeSourceField, SpectrumSourceField,
+                                    BackgroundSourceField)
 from ..utils import (from_currsys, quantify, check_keys, find_file,
                      figure_factory, get_logger)
 
@@ -160,17 +161,28 @@ class TERCurve(Effect):
     def background_source(self):
         if self._background_source is None:
             flux = self.emission
-            bg_hdu = fits.ImageHDU()
-
-            bg_hdu.header.update({"BG_SRC": True,
-                                  "BG_SURF": self.display_name,
-                                  "CUNIT1": "ARCSEC",
-                                  "CUNIT2": "ARCSEC",
-                                  "CDELT1": 0,
-                                  "CDELT2": 0,
-                                  "BUNIT": "PHOTLAM arcsec-2",
-                                  "SOLIDANG": "arcsec-2"})
-            self._background_source = Source(image_hdu=bg_hdu, spectra=flux)
+            bkg_spec = {0: flux}
+            bkg_hdr = fits.Header()
+
+            bkg_hdr.update({
+                "BG_SRC": True,
+                "BG_SURF": self.display_name,
+                "SPEC_REF": 0,
+                # Seem those are not needed...
+                # "CUNIT1": "ARCSEC",
+                # "CUNIT2": "ARCSEC",
+                # "CDELT1": 0,
+                # "CDELT2": 0,
+                "BUNIT": "PHOTLAM arcsec-2",
+                "SOLIDANG": "arcsec-2",
+            })
+            bkg_fld = BackgroundSourceField(field=None, spectra=bkg_spec, header=bkg_hdr)
+            bkg_src = Source(field=bkg_fld)
+
+            # Before BackgroundSourceField:
+            # bkg_hdu = fits.ImageHDU(header=bkg_hdr)
+            # bkg_src = Source(image_hdu=bkg_hdu, spectra=flux)
+            self._background_source = bkg_src
 
         return [self._background_source]
 
diff --git a/scopesim/optics/fov.py b/scopesim/optics/fov.py
index fefdd6b9..58c8b8dd 100644
--- a/scopesim/optics/fov.py
+++ b/scopesim/optics/fov.py
@@ -1,26 +1,35 @@
+# -*- coding: utf-8 -*-
 """Defines FieldOfView class."""
 
 from copy import deepcopy
 from itertools import chain
-from collections.abc import Iterable
+from collections.abc import Iterable, Generator
 
 import numpy as np
 from scipy.interpolate import interp1d
 
 from astropy import units as u
 from astropy.io import fits
-from astropy.table import Table
+from astropy.wcs import WCS
 from synphot import Empirical1D, SourceSpectrum
 from synphot.units import PHOTLAM
 
-from . import fov_utils as fu
 from . import image_plane_utils as imp_utils
+from ..source.source_fields import (
+    SourceField,
+    ImageSourceField,
+    CubeSourceField,
+    TableSourceField,
+    BackgroundSourceField,
+)
 
 from ..utils import (
     from_currsys,
     quantify,
     has_needed_keywords,
     get_logger,
+    array_minmax,
+    close_loop,
     unit_includes_per_physical_type,
 )
 from ..source.source import Source
@@ -36,16 +45,18 @@ class FieldOfView:
     Flux units after extracting the fields from the Source are in ph/s/pixel
 
     The initial header should contain an on-sky WCS description:
+
     - CDELT-, CUNIT-, NAXIS- : for pixel scale and size (assumed CUNIT in deg)
     - CRVAL-, CRPIX- : for positioning the final image
     - CTYPE- : is assumed to be "RA---TAN", "DEC--TAN"
 
     and an image-plane WCS description
+
     - CDELT-D, CUNIT-D, NAXISn : for pixel scale and size (assumed CUNIT in mm)
     - CRVAL-D, CRPIX-D : for positioning the final image
     - CTYPE-D : is assumed to be "LINEAR", "LINEAR"
 
-    The wavelength range is given by waverange
+    The wavelength range is given by waverange.
 
     """
 
@@ -60,11 +71,13 @@ class FieldOfView:
             "area": 0 * u.m**2,
             "pixel_area": None,  # [arcsec]
             "sub_pixel": "!SIM.sub_pixel.flag",
-            "distortion": {"scale": [1, 1],
-                           "offset": [0, 0],
-                           "shear": [1, 1],
-                           "rotation": 0,
-                           "radius_of_curvature": None},
+            "distortion": {
+                "scale": [1, 1],
+                "offset": [0, 0],
+                "shear": [1, 1],
+                "rotation": 0,
+                "radius_of_curvature": None,
+            },
             "conserve_image": True,
             "trace_id": None,
             "aperture_id": None,
@@ -73,7 +86,7 @@ class FieldOfView:
 
         self.cmds = cmds
 
-        if not any((has_needed_keywords(header, s) for s in ["", "S"])):
+        if not any(has_needed_keywords(header, s) for s in {"", "S"}):
             raise ValueError(
                 f"Header must contain a valid sky-plane WCS: {dict(header)}")
         if not has_needed_keywords(header, "D"):
@@ -90,7 +103,7 @@ class FieldOfView:
         self.hdu = None
 
         self.image_plane_id = 0
-        self.fields = []
+        self.fields: list[SourceField] = []
         self.spectra = {}
 
         # These are apparently not supposed to be used?
@@ -102,7 +115,8 @@ class FieldOfView:
         self._wavelength = None
         self._volume = None
 
-    def _pixarea(self, hdr):
+    @staticmethod
+    def _pixarea(hdr):
         return (hdr["CDELT1"] * u.Unit(hdr["CUNIT1"]) *
                 hdr["CDELT2"] * u.Unit(hdr["CUNIT2"])).to(u.arcsec ** 2)
 
@@ -113,70 +127,118 @@ class FieldOfView:
 
     @property
     def pixel_area(self):
+        """Return the area in arcsec**2 covered by one pixel."""
         if self.meta["pixel_area"] is None:
             # [arcsec] (really?)
             self.meta["pixel_area"] = self._pixarea(self.header).value
         return self.meta["pixel_area"]
 
-    def extract_from(self, src):
-        """..assumption: Bandpass has been applied.
+    def extract_from(self, src) -> None:
+        """
+        Extract relevent fields from source object.
 
-        .. note:: Spectra are cut and copied from the original Source object.
-            They are in original units. ph/s/pix comes in the make_**** methods
+        Parameters
+        ----------
+        src : Source
+            Input Source object to be "observed".
+
+        Returns
+        -------
+        None
+
+        Notes
+        -----
+        Spectra are cut and copied from the original Source object.
+        They are in original units. ph/s/pix comes in the make_**** methods.
+
+        This method assumes that Bandpass has been applied.
 
         """
         assert isinstance(src, Source), f"expected Source: {type(src)}"
 
-        fields_in_fov = [field.field for field in src.fields
-                         if fu.is_field_in_fov(self.header, field)]
+        fields_in_fov = list(self.get_fields_in_fov(src.fields))
+
         if not fields_in_fov:
             logger.warning("No fields in FOV.")
+        else:
+            logger.debug("%d fields in FOV", len(fields_in_fov))
+
+        corners_arcsec, _ = self.get_corners("arcsec")
+        corners_deg, _ = self.get_corners("deg")
+        minmax = array_minmax(corners_arcsec) * u.arcsec
+
+        # TODO: Not sure why it's necessary to recreate the fields here, but
+        #       perhasps for multi-fov instruments?
+        # TODO: Perhaps already split into cube, image and table fields here,
+        #       i.e. don't even bother with the total fields list and thus also
+        #       eleminate the need for the properties later on? It doesn't
+        #       seem like either the .fields attribute nor the various
+        #       fields properties are actually used outside the class itself...
+        # TODO: Perhaps remove .spectra and just keep spectra as source field
+        #       attributes, removing the need for SPEC_REF?
+        for field in fields_in_fov:
+            if isinstance(field, TableSourceField):
+                extracted = self.extract_area_from_table(field.field, minmax)
+                if not len(extracted):
+                    continue  # Can happen for small FOV and sparse table
+                # TODO: Rework extract_area_from_table to also affect spectra
+                #       and just return new copy of field.
+                new_fld = TableSourceField(
+                    field=extracted,
+                    spectra={
+                        ref: spec for ref, spec in field.spectra.items()
+                        if ref in extracted["ref"]
+                    },
+                )
+                self.fields.append(new_fld)
+                for ref, spec in new_fld.spectra.items():
+                    extracted = extract_range_from_spectrum(spec, self.waverange)
+                    self.spectra[ref] = extracted
+
+            elif isinstance(field, ImageSourceField):
+                assert field.header["NAXIS"] == 2, "Invalid image HDU"
+                extracted = self.extract_area_from_imagehdu(field.field, corners_deg)
+                replace_nans(extracted, self.cmds)
+                new_fld = ImageSourceField(
+                    field=extracted,
+                    spectra=field.spectra,
+                )
+                self.fields.append(new_fld)
+                ref, spec = new_fld.spectrum  # ImageField has only 1 spectrum
+                self.spectra[ref] = extract_range_from_spectrum(spec, self.waverange)
+
+            elif isinstance(field, CubeSourceField):
+                assert field.header["NAXIS"] == 3, "Invalid cube HDU"
+                extracted = self.extract_area_from_imagehdu(field.field, corners_deg)
+                replace_nans(extracted, self.cmds)
+                new_fld = CubeSourceField(field=extracted)
+                self.fields.append(new_fld)
+
+            elif isinstance(field, BackgroundSourceField):
+                self.fields.append(field)
+                ref, spec = field.spectrum  # BkgField has only 1 spectrum
+                self.spectra[ref] = extract_range_from_spectrum(spec, self.waverange)
 
-        spec_refs = set()
-        volume = self.volume()
-        for ifld, fld in enumerate(fields_in_fov):
-            if isinstance(fld, Table):
-                extracted_field = fu.extract_area_from_table(fld, volume)
-                spec_refs.update(extracted_field["ref"])
-                fields_in_fov[ifld] = extracted_field
-
-            elif isinstance(fld, fits.ImageHDU):
-                if fld.header["NAXIS"] in (2, 3):
-                    extracted = fu.extract_area_from_imagehdu(fld, volume)
-                    try:
-                        fill_value = self.cmds["!SIM.computing.nan_fill_value"]
-                    except TypeError:
-                        # This occurs in testing, not sure why
-                        fill_value = 0.
-                    if np.ma.fix_invalid(extracted.data, copy=False,
-                                         fill_value=fill_value).mask.any():
-                        logger.warning("Source contained NaN values, which "
-                                       "were replaced by %f.", fill_value)
-                        logger.info(
-                            "The replacement value for NaNs in sources can be "
-                            "set in '!SIM.computing.nan_fill_value'.")
-                    fields_in_fov[ifld] = extracted
-                if fld.header["NAXIS"] == 2 or fld.header.get("BG_SRC"):
-                    ref = fld.header.get("SPEC_REF")
-                    if ref is not None:
-                        spec_refs.add(ref)
-
-        waves = volume["waves"] * u.Unit(volume["wave_unit"])
-        spectra = {ref: fu.extract_range_from_spectrum(src.spectra[int(ref)], waves)
-                   for ref in spec_refs}
-
-        self.fields = fields_in_fov
-        self.spectra = spectra
-
-    def view(self, hdu_type="image", sub_pixel=None, use_photlam=None):
+            else:
+                raise TypeError(f"Unexpected source field type {type(field)}.")
+
+    def view(
+        self,
+        hdu_type: str = "image",
+        sub_pixel: bool | None = None,
+        use_photlam: bool | None = None,
+    ):
         """
         Force the self.fields to be viewed as a single object.
 
         Parameters
         ----------
-        sub_pixel : bool
-        hdu_type : str
-            ["cube", "image", "spectrum"]
+        hdu_type : {"image", "cube", "spectrum"}
+            DESCRIPTION.
+        sub_pixel : bool | None, optional
+            If None (the default), use value from meta.
+        use_photlam : bool | None, optional
+            If None (the default), assume False. Only used in imaging (why?).
 
         Returns
         -------
@@ -203,6 +265,184 @@ class FieldOfView:
             image = np.sum(self.hdu.data, axis=0)
             self.hdu.data = image
 
+    def get_fields_in_fov(self, fields: Iterable[SourceField]) -> Generator:
+        """Return True if Source.field footprint is inside FOV footprint."""
+        fov_corners, _ = self.get_corners("arcsec")
+
+        for field in fields:
+            field_corners = field.get_corners("arcsec")
+            is_inside_fov = (
+                (field_corners.max(axis=0) > fov_corners.min(axis=0)).all() and
+                (field_corners.min(axis=0) < fov_corners.max(axis=0)).all()
+            )
+            if is_inside_fov:
+                yield field
+
+    @staticmethod
+    def extract_area_from_table(table, minmax):
+        """
+        Extract the entries of a ``Table`` that fit inside the FOV volume.
+
+        Parameters
+        ----------
+        table : table.Table
+            The field table.
+        minmax : quantity
+            From FOV corners in the form of [[xmin, ymin], [xmax, ymax]].
+
+        Returns
+        -------
+        cut_table : table.Table
+            Table reduced to sources inside the FOV.
+
+        """
+        mask = ((table["x"].quantity >= minmax[0, 0]) *
+                (table["x"].quantity < minmax[1, 0]) *
+                (table["y"].quantity >= minmax[0, 1]) *
+                (table["y"].quantity < minmax[1, 1]))
+        return table[mask]
+
+    def extract_area_from_imagehdu(self, imagehdu, corners):
+        """
+        Extract the part of a ``ImageHDU`` that fits inside the FOV volume.
+
+        Parameters
+        ----------
+        imagehdu : fits.ImageHDU
+            The field ImageHDU, either an image or a cube with wavelength [um].
+        corners : quantity
+            From FOV corners in deg.
+
+        Returns
+        -------
+        new_imagehdu : fits.ImageHDU
+
+        """
+        # TODO: At some point streamline the debug logging here...
+        hdr = imagehdu.header
+        image_wcs = WCS(hdr, naxis=2)
+        naxis1, naxis2 = hdr["NAXIS1"], hdr["NAXIS2"]
+        logger.debug("old naxis: %s", [naxis1, naxis2])
+        xy_hdu = image_wcs.calc_footprint(center=False, axes=(naxis1, naxis2))
+
+        if image_wcs.wcs.cunit[0] == "deg":
+            logger.debug("Found 'deg' in image WCS, applying 360 fix.")
+            imp_utils._fix_360(xy_hdu)
+        elif image_wcs.wcs.cunit[0] == "arcsec":
+            logger.debug("Found 'arcsec' in image WCS, converting to deg.")
+            xy_hdu *= u.arcsec.to(u.deg)
+
+        logger.debug("XY HDU:\n%s", xy_hdu)
+        logger.debug("XY FOV:\n%s", corners)
+
+        xy0s = np.array((xy_hdu.min(axis=0), corners.min(axis=0))).max(axis=0)
+        xy1s = np.array((xy_hdu.max(axis=0), corners.max(axis=0))).min(axis=0)
+        logger.debug("xy0s: %s; xy1s: %s", xy0s, xy1s)
+
+        # Round to avoid floating point madness
+        xyp = image_wcs.wcs_world2pix(np.array([xy0s, xy1s]), 0).round(7)
+
+        # To deal with negative CDELTs
+        logger.debug("xyp:\n%s", xyp)
+        xyp.sort(axis=0)
+        logger.debug("xyp:\n%s", xyp)
+
+        xy0p = np.max(((0, 0), np.floor(xyp[0]).astype(int)), axis=0)
+        xy1p = np.min(((naxis1, naxis2), np.ceil(xyp[1]).astype(int)), axis=0)
+        logger.debug("xy0p: %s; xy1p: %s", xy0p, xy1p)
+
+        # Add 1 if the same
+        xy1p += (xy0p == xy1p)
+        logger.debug("xy0p: %s; xy1p: %s", xy0p, xy1p)
+
+        new_wcs, new_naxis = imp_utils.create_wcs_from_points(
+            np.array([xy0s, xy1s]).round(11), pixel_scale=hdr["CDELT1"])
+
+        # TODO: Come back at some point and figure out if the failing tests
+        #       here are relevant or can be ignored...
+        # FIXME: Commented out for now because it appears too often IRL...
+        # try:
+        #     roundtrip = new_wcs.wcs_world2pix(
+        #         np.array([xy0s, xy1s - .5*image_wcs.wcs.cdelt]), 0).round(5)
+        #     np.testing.assert_array_equal(roundtrip[0], [0, 0])
+        #     np.testing.assert_array_equal(roundtrip[1], new_naxis - [1, 1])
+        # except AssertionError:
+        #     logger.exception("WCS roundtrip assertion failed.")
+
+        new_hdr = new_wcs.to_header()
+        new_hdr.update({"NAXIS1": new_naxis[0], "NAXIS2": new_naxis[1]})
+
+        if hdr["NAXIS"] == 3:
+            new_hdr, data = self._extract_volume_from_cube(imagehdu, new_hdr,
+                                                           xy0p, xy1p)
+        else:
+            data = imagehdu.data[xy0p[1]:xy1p[1],
+                                 xy0p[0]:xy1p[0]]
+            new_hdr["SPEC_REF"] = hdr.get("SPEC_REF")
+
+        if not data.size:
+            logger.warning("Empty image HDU.")
+
+        return fits.ImageHDU(header=new_hdr, data=data)
+
+    def _extract_volume_from_cube(self, cubehdu, new_hdr, xy0p, xy1p):
+        hdr = cubehdu.header
+        # Look 0.5*wdel past the fov edges in each direction to catch any
+        # slices where the middle wavelength value doesn't fall inside the
+        # fov waverange, but up to 50% of the slice is actually inside the
+        # fov waverange:
+        # E.g. FOV: [1.92, 2.095], HDU bin centres: [1.9, 2.0, 2.1]
+        # CDELT3 = 0.1, and HDU bin edges: [1.85, 1.95, 2.05, 2.15]
+        # So 1.9 slice needs to be multiplied by 0.3, and 2.1 slice should be
+        # multipled by 0.45 to reach the scaled contribution of the edge slices
+        # This scaling factor is:
+        # f = ((hdu_bin_centre - fov_edge [+/-] 0.5 * cdelt3) % cdelt3) / cdelt3
+
+        hdu_waves = get_cube_waveset(hdr)
+        wdel = hdr["CDELT3"]
+        wunit = u.Unit(hdr.get("CUNIT3", "AA"))
+        fov_wmin, fov_wmax = (self.waverange).to(wunit).value
+        mask = ((hdu_waves > fov_wmin - 0.5 * wdel) *
+                (hdu_waves <= fov_wmax + 0.5 * wdel))  # need to go [+/-] half a bin
+
+        # OC [2021-12-14] if fov range is not covered by the source return nothing
+        if not np.any(mask):
+            logger.warning("FOV %s um - %s um: not covered by Source",
+                           fov_wmin, fov_wmax)
+            # FIXME: returning None here breaks the principle that a function
+            #        should always return the same type. Maybe this should
+            #        instead raise an exception that's caught higher up...
+            return None
+
+        i0p, i1p = np.where(mask)[0][0], np.where(mask)[0][-1]
+        f0 = (abs(hdu_waves[i0p] - fov_wmin + 0.5 * wdel) % wdel) / wdel    # blue edge
+        f1 = (abs(hdu_waves[i1p] - fov_wmax - 0.5 * wdel) % wdel) / wdel    # red edge
+        data = cubehdu.data[i0p:i1p+1,
+                            xy0p[1]:xy1p[1],
+                            xy0p[0]:xy1p[0]]
+        data[0, :, :] *= f0
+        if i1p > i0p:
+            data[-1, :, :] *= f1
+
+        # w0, w1 : the closest cube wavelengths outside the fov edge wavelengths
+        # fov_waves : the fov edge wavelengths
+        # f0, f1 : the scaling factors for the blue and red edge cube slices
+        #
+        # w0, w1 = hdu_waves[i0p], hdu_waves[i1p]
+
+        new_hdr.update({
+            "NAXIS": 3,
+            "NAXIS3": data.shape[0],
+            "CRVAL3": round(hdu_waves[i0p], 11),
+            "CRPIX3": 0,
+            "CDELT3": hdr["CDELT3"],
+            "CUNIT3": hdr["CUNIT3"],
+            "CTYPE3": hdr["CTYPE3"],
+            "BUNIT":  hdr["BUNIT"],
+        })
+
+        return new_hdr, data
+
     def _evaluate_spectrum_with_weight(self, ref, waveset, weight):
         return self.spectra[int(ref)](waveset).value * weight
 
@@ -222,8 +462,10 @@ class FieldOfView:
         * yield scaled flux to be added to canvas flux
         """
         for field in self.cube_fields:
-            hdu_waveset = fu.get_cube_waveset(field.header,
-                                              return_quantity=True)
+            # TODO: if SourceFields were kept in the _fields lists, the wave
+            #       attribute of CubeSourceField might be used directly (but
+            #       check if extraction limits the waves accordingly!!).
+            hdu_waveset = get_cube_waveset(field.header, return_quantity=True)
             fluxes = field.data.sum(axis=2).sum(axis=1)
             fov_waveset_fluxes = np.interp(fov_waveset, hdu_waveset, fluxes)
 
@@ -341,7 +583,7 @@ class FieldOfView:
             xpix, ypix = imp_utils.val2pix(self.header,
                                            field["x"] / 3600,
                                            field["y"] / 3600)
-            if from_currsys(self.meta["sub_pixel"], self.cmds):
+            if self.sub_pixel:
                 for idx, row in enumerate(field):
                     xs, ys, fracs = imp_utils.sub_pixel_fractions(xpix[idx],
                                                                   ypix[idx])
@@ -385,16 +627,14 @@ class FieldOfView:
             [ph s-1 pixel-1] or PHOTLAM (if use_photlam=True)
 
         """
-        spline_order = from_currsys("!SIM.computing.spline_order", self.cmds)
         # Make waveset and canvas image
         fov_waveset = self.waveset
         bin_widths = np.diff(fov_waveset)       # u.um
         bin_widths = 0.5 * (np.r_[0, bin_widths] + np.r_[bin_widths, 0])
-        area = from_currsys(self.meta["area"], self.cmds) << u.m**2   # u.m2
 
         # PHOTLAM * u.um * u.m2 --> ph / s
         specs = {ref: spec(fov_waveset) if use_photlam
-                 else (spec(fov_waveset) * bin_widths * area).to(u.ph / u.s)
+                 else (spec(fov_waveset) * bin_widths * self.area).to(u.ph / u.s)
                  for ref, spec in self.spectra.items()}
 
         fluxes = {ref: np.sum(spec.value) for ref, spec in specs.items()}
@@ -403,16 +643,16 @@ class FieldOfView:
             data=np.zeros((self.header["NAXIS2"], self.header["NAXIS1"])),
             header=self.header)
 
-        for tmp_hdu in chain(self._make_image_cubefields(area),
+        for tmp_hdu in chain(self._make_image_cubefields(self.area),
                              self._make_image_imagefields(fluxes)):
             canvas_image_hdu = imp_utils.add_imagehdu_to_imagehdu(
                 tmp_hdu,
                 canvas_image_hdu,
                 conserve_flux=True,
-                spline_order=spline_order)
+                spline_order=self.spline_order)
 
         for flux, weight, x, y in self._make_image_tablefields(fluxes):
-            if from_currsys(self.meta["sub_pixel"], self.cmds):
+            if self.sub_pixel:
                 # These x and y should not be arrays when sub_pixel is
                 # enabled, it is therefore not necessary to deploy the fix
                 # below in the else-branch.
@@ -454,8 +694,7 @@ class FieldOfView:
         for field in self.cube_fields:
             # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping
             # Assumption is that ImageHDUs have units of PHOTLAM arcsec-2
-            field_waveset = fu.get_cube_waveset(field.header,
-                                                return_quantity=True)
+            field_waveset = get_cube_waveset(field.header, return_quantity=True)
 
             # ..todo: Deal with this bounds_error in a more elegant way
             field_interp = interp1d(field_waveset.to(u.um).value,
@@ -521,7 +760,7 @@ class FieldOfView:
                 xpix, ypix = imp_utils.val2pix(self.header, xsky, ysky)
                 flux_vector = specs[row["ref"]].value * row["weight"] / self.pixel_area
 
-                if from_currsys(self.meta["sub_pixel"], self.cmds):
+                if self.sub_pixel:
                     xs, ys, fracs = imp_utils.sub_pixel_fractions(xpix, ypix)
                     for i, j, k in zip(xs, ys, fracs):
                         yield flux_vector * k, i, j
@@ -587,8 +826,6 @@ class FieldOfView:
             [ph s-1 AA-1 arcsec-2]      # as needed by SpectralTrace
 
         """
-        spline_order = from_currsys("!SIM.computing.spline_order", self.cmds)
-
         # 1. Make waveset and canvas cube (area, bin_width are applied at end)
         # TODO: Why is this not self.waveset? What's different?
         # -> For non-cube input but cube output (e.g. 2D image in spec mode),
@@ -626,11 +863,11 @@ class FieldOfView:
             canvas_cube_hdu = imp_utils.add_imagehdu_to_imagehdu(
                 field_hdu,
                 canvas_cube_hdu,
-                spline_order=spline_order)
+                spline_order=self.spline_order)
 
-        canvas_cube_hdu.data = sum(self._make_cube_imagefields(specs,
-                                                               spline_order),
-                                   start=canvas_cube_hdu.data)
+        canvas_cube_hdu.data = sum(self._make_cube_imagefields(
+            specs, self.spline_order),
+            start=canvas_cube_hdu.data)
 
         for flux, x, y in self._make_cube_tablefields(specs):
             # To prevent adding array values in this manner.
@@ -645,8 +882,7 @@ class FieldOfView:
         #    PHOTLAM = ph/s/cm-2/AA
         #    area = m2, fov_waveset = um
         # SpectralTrace wants ph/s/um/arcsec2 --> get rid of m2, leave um
-        area = from_currsys(self.meta["area"], self.cmds)  # u.m2
-        canvas_cube_hdu.data *= area.to(u.cm ** 2).value
+        canvas_cube_hdu.data *= self.area.to(u.cm ** 2).value
         canvas_cube_hdu.data *= 1e4       # ph/s/AA/arcsec2 --> ph/s/um/arcsec2
 
         # TODO: what's with this code??
@@ -660,22 +896,9 @@ class FieldOfView:
                                        "CRPIX3": 1,
                                        "CUNIT3": "um",
                                        "CTYPE3": "WAVE"})
-        # ..todo: Add the log wavelength keyword here, if log scale is needed
+        # TODO: Add the log wavelength keyword here, if log scale is needed
         return canvas_cube_hdu      # [ph s-1 AA-1 (arcsec-2)]
 
-    def volume(self, wcs_prefix=""):
-        xy = imp_utils.calc_footprint(self.header, wcs_suffix=wcs_prefix)
-        unit = self.header[f"CUNIT1{wcs_prefix}"].lower()
-        # FIXME: This is unused!!
-        # wave_corners = self.waverange
-        minmax = np.array((xy.min(axis=0), xy.max(axis=0)))
-        self._volume = {"xs": minmax[:, 0],
-                        "ys": minmax[:, 1],
-                        "waves": self.waverange,
-                        "xy_unit": unit,
-                        "wave_unit": "um"}
-        return self._volume
-
     @property
     def data(self):
         """Return either hdu.data, image, cube, spectrum or None."""
@@ -689,14 +912,9 @@ class FieldOfView:
             return self.spectrum
         return None
 
-    @property
-    def corners(self):
+    def get_corners(self, new_unit: str = None):
         """Return sky footprint, image plane footprint."""
-        # Couldn't find where this is ever used, put warning here just in case
-        logger.warning("calc_footprint has been updated, any code that "
-                       "relies on this .corners property must likely be "
-                       "adapted as well.")
-        sky_corners = imp_utils.calc_footprint(self.header)
+        sky_corners = imp_utils.calc_footprint(self.header, new_unit=new_unit)
         imp_corners = imp_utils.calc_footprint(self.header, "D")
         return sky_corners, imp_corners
 
@@ -704,17 +922,17 @@ class FieldOfView:
     def waverange(self):
         """Return wavelength range in um [wave_min, wave_max]."""
         if self._waverange is None:
-            wave_min = quantify(self.meta["wave_min"], u.um).value
-            wave_max = quantify(self.meta["wave_max"], u.um).value
-            self._waverange = [wave_min, wave_max]
+            wave_min = self.meta["wave_min"]
+            wave_max = self.meta["wave_max"]
+            self._waverange = [wave_min, wave_max] << u.um
         return self._waverange
 
     @property
     def wavelength(self):
         """Return central wavelength in um."""
         if self._wavelength is None:
-            self._wavelength = np.average(self.waverange)
-        return quantify(self._wavelength, u.um)
+            self._wavelength = self.waverange.mean() << u.um
+        return self._wavelength
 
     @property
     def waveset(self):
@@ -722,13 +940,15 @@ class FieldOfView:
         if field_cubes := self.cube_fields:
             naxis3_max = np.argmax([cube.header["NAXIS3"]
                                     for cube in field_cubes])
-            _waveset = fu.get_cube_waveset(field_cubes[naxis3_max].header,
-                                           return_quantity=True)
+            _waveset = get_cube_waveset(
+                field_cubes[naxis3_max].header,
+                return_quantity=True,
+            )
         elif self.spectra:
             wavesets = [spec.waveset for spec in self.spectra.values()]
             _waveset = np.concatenate(wavesets)
         else:
-            _waveset = self.waverange * u.um
+            _waveset = self.waverange << u.um
 
         # TODO: tie the round to a global precision setting (this is defined
         #       better in another TODO somewhere...)
@@ -740,34 +960,48 @@ class FieldOfView:
         _waveset = np.unique(_waveset.to(u.um).round(10))
         return _waveset
 
+    @property
+    def area(self) -> float:
+        """Return meta["area"] from cmds."""
+        # TODO: Make this an attribute once meta resolving is implemented.
+        return from_currsys(self.meta["area"], self.cmds) << u.m**2
+
+    @property
+    def sub_pixel(self) -> bool:
+        """Return meta["sub_pixel"] from cmds."""
+        # TODO: Make this an attribute once meta resolving is implemented.
+        return from_currsys(self.meta["sub_pixel"], self.cmds)
+
+    @property
+    def spline_order(self) -> int:
+        """Return "!SIM.computing.spline_order" from cmds."""
+        # TODO: Make this an attribute once meta resolving is implemented.
+        return from_currsys("!SIM.computing.spline_order", self.cmds)
+
     @property
     def cube_fields(self):
         """Return list of non-BG_SRC ImageHDU fields with NAXIS=3."""
-        return [field for field in self.fields
-                if isinstance(field, fits.ImageHDU)
-                and field.header["NAXIS"] == 3
-                and not field.header.get("BG_SRC", False)]
+        return [field.field for field in self.fields
+                if isinstance(field, CubeSourceField)]
 
     @property
     def image_fields(self):
         """Return list of non-BG_SRC ImageHDU fields with NAXIS=2."""
-        return [field for field in self.fields
-                if isinstance(field, fits.ImageHDU)
-                and field.header["NAXIS"] == 2
-                and not field.header.get("BG_SRC", False)]
+        return [field.field for field in self.fields
+                if isinstance(field, ImageSourceField)]
 
     @property
     def table_fields(self):
         """Return list of Table fields."""
-        return [field for field in self.fields
-                if isinstance(field, Table)]
+        return [field.field for field in self.fields
+                if isinstance(field, TableSourceField)]
 
     @property
     def background_fields(self):
         """Return list of BG_SRC ImageHDU fields."""
-        return [field for field in self.fields
-                if isinstance(field, fits.ImageHDU)
-                and field.header.get("BG_SRC", False)]
+        # HACK: While Sourcefields are not fully supported.
+        return [fits.ImageHDU(header=field.header) for field in self.fields
+                if isinstance(field, BackgroundSourceField)]
 
     def _ensure_deg_header(self):
         cunit = u.Unit(self.header["CUNIT1"].lower())
@@ -779,19 +1013,88 @@ class FieldOfView:
         self.header["CUNIT1"] = "deg"
         self.header["CUNIT2"] = "deg"
 
-    def __repr__(self):
-        waverange = [self.meta["wave_min"].value, self.meta["wave_max"].value]
-        msg = (f"{self.__class__.__name__}({self.header!r}, {waverange!r}, "
-               f"{self.detector_header!r}, **{self.meta!r})")
+    def plot(self, axes, units: str = "arcsec") -> None:
+        """Plot FOV footprint."""
+        outline = np.array(list(close_loop(self.get_corners(units)[0])))
+        axes.plot(*outline.T, label=f"FOV id: {self.meta['id']}")
+
+    def __repr__(self) -> str:
+        """Return repr(self)."""
+        msg = (f"{self.__class__.__name__}({self.header!r}, "
+               f"{self.waverange!r}, {self.detector_header!r}, "
+               f"**{self.meta!r})")
         return msg
 
-    def __str__(self):
+    def __str__(self) -> str:
+        """Return str(self)."""
         msg = (f"FOV id: {self.meta['id']}, with dimensions "
                f"({self.header['NAXIS1']}, {self.header['NAXIS2']})\n"
                f"Sky centre: ({self.header['CRVAL1']}, "
                f"{self.header['CRVAL2']})\n"
                f"Image centre: ({self.header['CRVAL1D']}, "
                f"{self.header['CRVAL2D']})\n"
-               f"Wavelength range: ({self.meta['wave_min']}, "
-               f"{self.meta['wave_max']})um\n")
+               f"Wavelength range: {self.waverange!s}\n")
         return msg
+
+
+def replace_nans(field, cmds) -> None:
+    """Replace NaN values in image field."""
+    try:
+        fill_value = cmds["!SIM.computing.nan_fill_value"]
+    except TypeError:
+        # This occurs in testing, not sure why
+        fill_value = 0.
+    if np.ma.fix_invalid(field.data, copy=False,
+                         fill_value=fill_value).mask.any():
+        logger.warning("Source contained NaN values, which "
+                       "were replaced by %f.", fill_value)
+        logger.info(
+            "The replacement value for NaNs in sources can be "
+            "set in '!SIM.computing.nan_fill_value'.")
+
+
+def get_cube_waveset(hdr, return_quantity=False):
+    wval, wdel, wpix, wlen, = [hdr[kw] for kw in ["CRVAL3", "CDELT3",
+                                                  "CRPIX3", "NAXIS3"]]
+    # ASSUMPTION - cube wavelength is in regularly spaced units of um
+    wmin = wval - wdel * wpix
+    wmax = wmin + wdel * (wlen - 1)
+    wunit = u.Unit(hdr.get("CUNIT3", "AA"))
+
+    if "LOG" in hdr.get("CTYPE3", ""):
+        hdu_waves = np.logspace(wmin, wmax, wlen)
+    else:
+        hdu_waves = np.linspace(wmin, wmax, wlen)
+
+    if return_quantity:
+        hdu_waves = hdu_waves << wunit
+        hdu_waves.to(u.um)
+
+    return hdu_waves
+
+
+def extract_range_from_spectrum(spectrum, waverange):
+    assert isinstance(spectrum, SourceSpectrum), (
+        f"spectrum must be of type synphot.SourceSpectrum: {type(spectrum)}")
+
+    wave_min, wave_max = quantify(waverange, u.um).to(u.AA).value
+    spec_waveset = spectrum.waveset.to(u.AA).value
+    mask = (spec_waveset > wave_min) * (spec_waveset < wave_max)
+
+    # FIXME: Why did I comment this out in 2023? Seems useful to have...
+    # if sum(mask) == 0:
+    #     logger.info(
+    #         "Waverange does not overlap with Spectrum waveset: %s <> %s for "
+    #         "spectrum %s", [wave_min, wave_max], spec_waveset, spectrum)
+    if wave_min < min(spec_waveset) or wave_max > max(spec_waveset):
+        logger.info(("Waverange only partially overlaps with Spectrum waveset: "
+                     "%s <> %s for spectrum %s"),
+                     [wave_min, wave_max], spec_waveset, spectrum)
+
+    wave = np.r_[wave_min, spec_waveset[mask], wave_max]
+    flux = spectrum(wave)
+
+    new_spectrum = SourceSpectrum(Empirical1D, points=wave, lookup_table=flux)
+    new_spectrum.meta.update(spectrum.meta)
+
+    return new_spectrum
diff --git a/scopesim/optics/fov_utils.py b/scopesim/optics/fov_utils.py
index a69a91c9..75a1f788 100644
--- a/scopesim/optics/fov_utils.py
+++ b/scopesim/optics/fov_utils.py
@@ -1,77 +1,15 @@
+# -*- coding: utf-8 -*-
 
 import numpy as np
 from astropy import units as u
 from astropy.io import fits
 from astropy.table import Table, Column
-from astropy.wcs import WCS
-from synphot import SourceSpectrum, Empirical1D
 
 from . import image_plane_utils as imp_utils
-from ..source.source_fields import SourceField
-from ..utils import from_currsys, quantify, quantity_from_table, get_logger
-
-
-logger = get_logger(__name__)
-
-
-def is_field_in_fov(fov_header, field, wcs_suffix=""):
-    """
-    Return True if Source.field footprint is inside the FieldOfView footprint.
-
-    Parameters
-    ----------
-    fov_header : fits.Header
-        Header from a FieldOfView object
-    field : [SourceField, astropy.Table, astropy.ImageHDU]
-        Field object from a Source object. Should now be SourceField, but bare
-        Table and ImageHDU still supported.
-    wcs_suffix : str
-        ["S", "D"] Coordinate system: Sky or Detector
-
-    Returns
-    -------
-    is_inside_fov : bool
-
-    """
-    if isinstance(field, SourceField):
-        field = field.field
-    # TODO: Check if this can always get a SourceField, if yes then just do
-    #       that and remove this.
-
-    if isinstance(field, fits.ImageHDU) and \
-            field.header.get("BG_SRC") is not None:
-        is_inside_fov = True
-    else:
-        if isinstance(field, Table):
-            x = list(quantity_from_table("x", field,
-                                               u.arcsec).to(u.deg).value)
-            y = list(quantity_from_table("y", field,
-                                               u.arcsec).to(u.deg).value)
-            s = wcs_suffix
-            # cdelt = quantify(fov_header["CDELT1" + s], u.deg).value
-            cdelt = fov_header[f"CDELT1{s}"] * u.Unit(fov_header[f"CUNIT1{s}"]).to(u.deg)
-            field_header = imp_utils.header_from_list_of_xy(x, y, cdelt, s)
-        elif isinstance(field, (fits.ImageHDU, fits.PrimaryHDU)):
-            field_header = field.header
-        else:
-            logger.warning("Input was neither Table nor ImageHDU: %s", field)
-            return False
-
-        xy = imp_utils.calc_footprint(field_header, wcs_suffix)
-        ext_xsky, ext_ysky = xy[:, 0], xy[:, 1]
-        xy = imp_utils.calc_footprint(fov_header, wcs_suffix)
-        fov_xsky, fov_ysky = xy[:, 0], xy[:, 1]
-        fov_xsky *= u.Unit(fov_header["CUNIT1"].lower()).to(u.deg)
-        fov_ysky *= u.Unit(fov_header["CUNIT2"].lower()).to(u.deg)
-
-        is_inside_fov = (min(ext_xsky) < max(fov_xsky) and
-                         max(ext_xsky) > min(fov_xsky) and
-                         min(ext_ysky) < max(fov_ysky) and
-                         max(ext_ysky) > min(fov_ysky))
-
-    return is_inside_fov
+from ..utils import from_currsys, quantity_from_table
 
 
+# TODO: Unused function. Remove?
 def make_flux_table(source_tbl, src, wave_min, wave_max, area):
     fluxes = np.zeros(len(src.spectra))
     ref_set = list(set(source_tbl["ref"]))
@@ -90,6 +28,7 @@ def make_flux_table(source_tbl, src, wave_min, wave_max, area):
     return tbl
 
 
+# TODO: Unused function. Remove?
 def combine_table_fields(fov_header, src, field_indexes):
     """
     Combine list of Table objects into a single one bounded by the Header WCS.
@@ -138,6 +77,7 @@ def combine_table_fields(fov_header, src, field_indexes):
     return tbl
 
 
+# TODO: Unused functionv (except in a comment somewhere). Remove?
 def combine_imagehdu_fields(fov_header, src, fields_indexes, wave_min, wave_max,
                             area, wcs_suffix="", cmds=None):
     """
@@ -196,264 +136,7 @@ def combine_imagehdu_fields(fov_header, src, fields_indexes, wave_min, wave_max,
     return canvas_hdu
 
 
-def sky2fp(header, xsky, ysky):
-    """
-    Convert sky coordinates to image plane coordinated.
-
-    Parameters
-    ----------
-    header : Header
-        Header of a FieldOfView object which contains two sets of WCS keywords
-    xsky, ysky : float, array
-        [deg] The on-sky coordinated
-
-    Returns
-    -------
-    xdet, ydet : float, array
-        [mm] The coordinated on the image plane
-
-    """
-    xpix, ypix = imp_utils.val2pix(header, xsky, ysky)
-    xdet, ydet = imp_utils.pix2val(header, xpix, ypix, "D")
-    return xdet, ydet
-
-
-def extract_common_field(field, fov_volume):
-    """
-    Extract the overlapping parts of a field within a FOV volume.
-
-    Parameters
-    ----------
-    field : Table or ImageHDU
-    fov_volume : dict
-        Contains {"xs": [xmin, xmax], "ys": [ymin, ymax],
-                  "waves": [wave_min, wave_max],
-                  "xy_unit": "deg" or "mm", "wave_unit": "um"}
-
-    Returns
-    -------
-    field_new : Table or ImageHDU
-
-    """
-    if isinstance(field, Table):
-        field_new = extract_area_from_table(field, fov_volume)
-    elif isinstance(field, fits.ImageHDU):
-        field_new = extract_area_from_imagehdu(field, fov_volume)
-    else:
-        raise ValueError("field must be either Table or ImageHDU, but is "
-                         f"{type(field)}")
-
-    return field_new
-
-
-def extract_area_from_table(table, fov_volume):
-    """
-    Extract the entries of a ``Table`` that fits inside the `fov_volume`.
-
-    Parameters
-    ----------
-    table : fits.ImageHDU
-        The field ImageHDU, either an image of a wavelength [um] cube
-    fov_volume : dict
-        Contains {"xs": [xmin, xmax], "ys": [ymin, ymax],
-                  "waves": [wave_min, wave_max],
-                  "xy_unit": "deg" or "mm", "wave_unit": "um"}
-
-    Returns
-    -------
-    new_imagehdu : fits.ImageHDU
-
-    """
-    fov_unit = u.Unit(fov_volume["xy_unit"])
-    fov_xs = (fov_volume["xs"] * fov_unit).to(table["x"].unit)
-    fov_ys = (fov_volume["ys"] * fov_unit).to(table["y"].unit)
-
-    mask = ((table["x"].data >= fov_xs[0].value) *
-            (table["x"].data < fov_xs[1].value) *
-            (table["y"].data >= fov_ys[0].value) *
-            (table["y"].data < fov_ys[1].value))
-    table_new = table[mask]
-
-    return table_new
-
-
-def extract_area_from_imagehdu(imagehdu, fov_volume):
-    """
-    Extract the part of a ``ImageHDU`` that fits inside the `fov_volume`.
-
-    Parameters
-    ----------
-    imagehdu : fits.ImageHDU
-        The field ImageHDU, either an image or a cube with wavelength [um]
-    fov_volume : dict
-        Contains {"xs": [xmin, xmax], "ys": [ymin, ymax],
-                  "waves": [wave_min, wave_max],
-                  "xy_unit": "deg" or "mm", "wave_unit": "um"}
-
-    Returns
-    -------
-    new_imagehdu : fits.ImageHDU
-
-    """
-    hdr = imagehdu.header
-    image_wcs = WCS(hdr, naxis=2)
-    naxis1, naxis2 = hdr["NAXIS1"], hdr["NAXIS2"]
-    xy_hdu = image_wcs.calc_footprint(center=False, axes=(naxis1, naxis2))
-
-    if image_wcs.wcs.cunit[0] == "deg":
-        logger.debug("Found 'deg' in image WCS, applying 360 fix.")
-        imp_utils._fix_360(xy_hdu)
-    elif image_wcs.wcs.cunit[0] == "arcsec":
-        logger.debug("Found 'arcsec' in image WCS, converting to deg.")
-        xy_hdu *= u.arcsec.to(u.deg)
-    logger.debug("XY HDU:\n%s", xy_hdu)
-
-    xy_fov = np.array([fov_volume["xs"], fov_volume["ys"]]).T
-
-    if fov_volume["xy_unit"] == "arcsec":
-        xy_fov *= u.arcsec.to(u.deg)
-
-    logger.debug("XY FOV:\n%s", xy_fov)
-
-    xy0s = np.array((xy_hdu.min(axis=0), xy_fov.min(axis=0))).max(axis=0)
-    xy1s = np.array((xy_hdu.max(axis=0), xy_fov.max(axis=0))).min(axis=0)
-
-    # Round to avoid floating point madness
-    xyp = image_wcs.wcs_world2pix(np.array([xy0s, xy1s]), 0).round(7)
-
-    # To deal with negative CDELTs
-    logger.debug("xyp:\n%s", xyp)
-    xyp.sort(axis=0)
-    logger.debug("xyp:\n%s", xyp)
-
-    xy0p = np.max(((0, 0), np.floor(xyp[0]).astype(int)), axis=0)
-    xy1p = np.min(((naxis1, naxis2), np.ceil(xyp[1]).astype(int)), axis=0)
-
-    # Add 1 if the same
-    xy1p += (xy0p == xy1p)
-
-    logger.debug("xy0p: %s; xy1p: %s", xy0p, xy1p)
-
-    new_wcs, new_naxis = imp_utils.create_wcs_from_points(
-        np.array([xy0s, xy1s]), pixel_scale=hdr["CDELT1"])
-    new_hdr = new_wcs.to_header()
-    new_hdr.update({"NAXIS1": new_naxis[0], "NAXIS2": new_naxis[1]})
-
-    if hdr["NAXIS"] == 3:
-
-        # Look 0.5*wdel past the fov edges in each direction to catch any
-        # slices where the middle wavelength value doesn't fall inside the
-        # fov waverange, but up to 50% of the slice is actually inside the
-        # fov waverange:
-        # E.g. FOV: [1.92, 2.095], HDU bin centres: [1.9, 2.0, 2.1]
-        # CDELT3 = 0.1, and HDU bin edges: [1.85, 1.95, 2.05, 2.15]
-        # So 1.9 slice needs to be multiplied by 0.3, and 2.1 slice should be
-        # multipled by 0.45 to reach the scaled contribution of the edge slices
-        # This scaling factor is:
-        # f = ((hdu_bin_centre - fov_edge [+/-] 0.5 * cdelt3) % cdelt3) / cdelt3
-
-        hdu_waves = get_cube_waveset(hdr)
-        wdel = hdr["CDELT3"]
-        wunit = u.Unit(hdr.get("CUNIT3", "AA"))
-        fov_waves = quantify(fov_volume["waves"], u.um).to(wunit).value
-        mask = ((hdu_waves > fov_waves[0] - 0.5 * wdel) *
-                (hdu_waves <= fov_waves[1] + 0.5 * wdel))  # need to go [+/-] half a bin
-
-        # OC [2021-12-14] if fov range is not covered by the source return nothing
-        if not np.any(mask):
-            logger.warning("FOV %s um - %s um: not covered by Source",
-                           fov_waves[0], fov_waves[1])
-            # FIXME: returning None here breaks the principle that a function
-            #        should always return the same type. Maybe this should
-            #        instead raise an exception that's caught higher up...
-            return None
-
-        i0p, i1p = np.where(mask)[0][0], np.where(mask)[0][-1]
-        f0 = (abs(hdu_waves[i0p] - fov_waves[0] + 0.5 * wdel) % wdel) / wdel    # blue edge
-        f1 = (abs(hdu_waves[i1p] - fov_waves[1] - 0.5 * wdel) % wdel) / wdel    # red edge
-        data = imagehdu.data[i0p:i1p+1,
-                             xy0p[1]:xy1p[1],
-                             xy0p[0]:xy1p[0]]
-        data[0, :, :] *= f0
-        if i1p > i0p:
-            data[-1, :, :] *= f1
-
-        # w0, w1 : the closest cube wavelengths outside the fov edge wavelengths
-        # fov_waves : the fov edge wavelengths
-        # f0, f1 : the scaling factors for the blue and red edge cube slices
-        #
-        # w0, w1 = hdu_waves[i0p], hdu_waves[i1p]
-
-        new_hdr.update({"NAXIS": 3,
-                        "NAXIS3": data.shape[0],
-                        "CRVAL3": hdu_waves[i0p],
-                        "CRPIX3": 0,
-                        "CDELT3": hdr["CDELT3"],
-                        "CUNIT3": hdr["CUNIT3"],
-                        "CTYPE3": hdr["CTYPE3"],
-                        "BUNIT":  hdr["BUNIT"]})
-
-    else:
-        data = imagehdu.data[xy0p[1]:xy1p[1],
-                             xy0p[0]:xy1p[0]]
-        new_hdr["SPEC_REF"] = hdr.get("SPEC_REF")
-
-    if not data.size:
-        logger.warning("Empty image HDU.")
-
-    new_imagehdu = fits.ImageHDU(data=data)
-    new_imagehdu.header.update(new_hdr)
-
-    return new_imagehdu
-
-
-def get_cube_waveset(hdr, return_quantity=False):
-    wval, wdel, wpix, wlen, = [hdr[kw] for kw in ["CRVAL3", "CDELT3",
-                                                  "CRPIX3", "NAXIS3"]]
-    # ASSUMPTION - cube wavelength is in regularly spaced units of um
-    wmin = wval - wdel * wpix
-    wmax = wmin + wdel * (wlen - 1)
-    wunit = u.Unit(hdr.get("CUNIT3", "AA"))
-
-    if "LOG" in hdr.get("CTYPE3", ""):
-        hdu_waves = np.logspace(wmin, wmax, wlen)
-    else:
-        hdu_waves = np.linspace(wmin, wmax, wlen)
-
-    if return_quantity:
-        hdu_waves = hdu_waves << wunit
-        hdu_waves.to(u.um)
-
-    return hdu_waves
-
-
-def extract_range_from_spectrum(spectrum, waverange):
-    if not isinstance(spectrum, SourceSpectrum):
-        raise ValueError(f"spectrum must be of type synphot.SourceSpectrum: "
-                         f"{type(spectrum)}")
-
-    wave_min, wave_max = quantify(waverange, u.um).to(u.AA).value
-    spec_waveset = spectrum.waveset.to(u.AA).value
-    mask = (spec_waveset > wave_min) * (spec_waveset < wave_max)
-
-    # if sum(mask) == 0:
-    #     logger.info(
-    #         "Waverange does not overlap with Spectrum waveset: %s <> %s for "
-    #         "spectrum %s", [wave_min, wave_max], spec_waveset, spectrum)
-    if wave_min < min(spec_waveset) or wave_max > max(spec_waveset):
-        logger.info(("Waverange only partially overlaps with Spectrum waveset: "
-                     "%s <> %s for spectrum %s"),
-                     [wave_min, wave_max], spec_waveset, spectrum)
-
-    wave = np.r_[wave_min, spec_waveset[mask], wave_max]
-    flux = spectrum(wave)
-
-    new_spectrum = SourceSpectrum(Empirical1D, points=wave, lookup_table=flux)
-    new_spectrum.meta.update(spectrum.meta)
-
-    return new_spectrum
-
-
+# TODO: Unused function. Remove?
 def make_cube_from_table(table, spectra, waveset, fov_header, sub_pixel=False):
     """
 
diff --git a/scopesim/optics/image_plane_utils.py b/scopesim/optics/image_plane_utils.py
index a188641d..d0742856 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)
+    offset = (extent - naxis) * pixel_scale  # compensate rounding
 
     # FIXME: Woule be nice to have D headers referenced at (1, 1), but that
     #        breaks some things, likely down to rescaling...
@@ -202,13 +203,14 @@ def create_wcs_from_points(points: np.ndarray,
                                         f"pix, got '{naxis.unit}' instead.")
         naxis = naxis.value
 
-    if (naxis == 0).any():
+    if (zeroaxis := (naxis == 0)).any():
         # Ensure at least one pixel in each axis.
         logger.debug("NAXISn == 0, using minimum of 1.")
-        naxis[naxis == 0] = 1
+        naxis[zeroaxis] = 1
+        offset[zeroaxis] = 0
 
     crpix = (naxis + 1) / 2
-    crval = (points.min(axis=0) + points.max(axis=0)) / 2
+    crval = (points.min(axis=0) + points.max(axis=0) + offset) / 2
 
     # Cannot do `in "DX"` here because that would also match the empty string.
     linsuff = {"D", "X"}
diff --git a/scopesim/source/source.py b/scopesim/source/source.py
index a0f11399..c040208c 100644
--- a/scopesim/source/source.py
+++ b/scopesim/source/source.py
@@ -155,7 +155,7 @@ class Source:
 
     def __init__(self, filename=None, cube=None, ext=0,
                  lam=None, spectra=None, x=None, y=None, ref=None, weight=None,
-                 table=None, image_hdu=None, flux=None, **kwargs):
+                 table=None, image_hdu=None, flux=None, field=None, **kwargs):
 
         self._meta = {}
         self.fields: list[SourceField] = []
@@ -200,6 +200,9 @@ class Source:
                 ref is not None and spectra is not None:
             self._from_arrays(x, y, ref, weight, spectra, **kwargs)
 
+        elif field is not None:
+            self._from_field(field)
+
         self.meta.update(kwargs)
 
     def _from_file(self, filename, spectra, flux, **kwargs):
@@ -374,6 +377,11 @@ class Source:
 
         self.fields = [CubeSourceField(cube_hdu, meta=kwargs)]
 
+    def _from_field(self, field: SourceField):
+        """Create from single source field."""
+        assert not self.fields, "Constructor method must act on empty instance!"
+        self.fields = [field]
+
     def _get_fields(self, subclass):
         """Yield fields of specific subclass."""
         for field in self.fields:
@@ -601,10 +609,10 @@ class Source:
                 field.field["ref"] += specrefoffset
                 self.fields.append(field)
 
-            elif isinstance(field, HDUSourceField):
+            else:
                 if ("SPEC_REF" in field.header and
                         isinstance(field.header["SPEC_REF"], int)):
-                    field.field.header["SPEC_REF"] += specrefoffset
+                    field.header["SPEC_REF"] += specrefoffset
                 self.fields.append(field)
 
             field.spectra = {k + specrefoffset: v for k, v in
diff --git a/scopesim/source/source_fields.py b/scopesim/source/source_fields.py
index 3f1cd4ef..ae7282de 100644
--- a/scopesim/source/source_fields.py
+++ b/scopesim/source/source_fields.py
@@ -112,6 +112,10 @@ class SourceField:
         """Name of the object (if set)."""
         return self.meta.get("object", "<unknown>")
 
+    def get_corners(self, unit: u.Unit | str = "arcsec") -> np.ndarray:
+        """Calculate and return footprint corner points in `unit`."""
+        raise NotImplementedError("Subclasses should implement this.")
+
 
 @dataclass
 class SpectrumSourceField(SourceField):
@@ -119,6 +123,14 @@ class SpectrumSourceField(SourceField):
 
     spectra: dict
 
+    @property
+    def spectrum(self) -> SourceSpectrum:
+        """Return single spectrum and ref if only one spectrum in spectra."""
+        if len(self.spectra) > 1:
+            raise TypeError("More than one spectrum in field -> use spectra!")
+        spec, = self.spectra.items()
+        return spec
+
 
 @dataclass
 class TableSourceField(SpectrumSourceField):
@@ -181,6 +193,20 @@ class TableSourceField(SpectrumSourceField):
         stream.write(f"Table with {len(self)} rows, referencing "
                      f"spectra {set(self.spectra)}")
 
+    def get_corners(self, unit: u.Unit | str = "arcsec") -> np.ndarray:
+        """Calculate and return footprint corner points in `unit`."""
+        x_qty = self.field["x"].quantity.to(unit).value
+        y_qty = self.field["y"].quantity.to(unit).value
+        x_min, x_max = x_qty.min(), x_qty.max()
+        y_min, y_max = y_qty.min(), y_qty.max()
+        corners = np.array([
+            [x_min, y_min],
+            [x_min, y_max],
+            [x_max, y_min],
+            [x_max, y_max],
+        ])
+        return corners
+
     def plot(self, axes, color) -> None:
         """Plot source."""
         axes.plot(self.field["x"], self.field["y"], color+".", label=self.name)
@@ -243,11 +269,13 @@ class HDUSourceField(SourceField):
         stream.write(f"ImageHDU with size {self.img_size}, referencing "
                      f"spectrum {self.field.header.get('SPEC_REF', '-')}")
 
+    def get_corners(self, unit: u.Unit | str = "arcsec") -> np.ndarray:
+        """Calculate and return footprint corner points in `unit`."""
+        return imp_utils.calc_footprint(self.header, new_unit=unit)
+
     def plot(self, axes, color) -> None:
         """Plot source."""
-        xypts = imp_utils.calc_footprint(self.header)
-        convf = u.Unit(self.header["CUNIT1"]).to(u.arcsec)
-        outline = np.array(list(close_loop(xypts))) * convf
+        outline = np.array(list(close_loop(self.get_corners())))
         axes.plot(*outline.T, color, label=self.name)
 
     def shift(self, dx, dy) -> None:
@@ -286,7 +314,10 @@ class CubeSourceField(HDUSourceField):
             self.wcs = WCS(self.field)
 
         try:
-            bunit = u.Unit(self.header["BUNIT"])
+            bunit = str(self.header["BUNIT"])
+            # Can't just do .lower because some units are uppercase (e.g. J)
+            bunit = bunit.replace("PHOTLAM", "photlam")
+            bunit = u.Unit(bunit)
         except KeyError:
             bunit = u.erg / u.s / u.cm**2 / u.arcsec**2
             logger.warning(
@@ -318,3 +349,17 @@ class CubeSourceField(HDUSourceField):
         with u.set_enabled_equivalencies(u.spectral()):
             wave = swcs.pixel_to_world(np.arange(swcs.pixel_shape[0])) << u.um
         return wave
+
+
+@dataclass
+class BackgroundSourceField(SpectrumSourceField):
+    """Source field with spectrum only, for TER curve emissions."""
+
+    header: fits.Header
+
+    def get_corners(self, unit: u.Unit | str = "arcsec") -> np.ndarray:
+        """Return imaginary corner from + to - infinity."""
+        return np.array([-np.inf, np.inf])
+
+    def _write_stream(self, stream: TextIO) -> None:
+        stream.write(f"Background field ref. spectrum {self.spectrum[0]}")
diff --git a/scopesim/utils.py b/scopesim/utils.py
index 78ed0aff..7421676f 100644
--- a/scopesim/utils.py
+++ b/scopesim/utils.py
@@ -134,6 +134,11 @@ def deriv_polynomial2d(poly):
     return dpoly_dx, dpoly_dy
 
 
+def array_minmax(arr):
+    """Create min and max from array of points, useful for footprint."""
+    return np.array((arr.min(axis=0), arr.max(axis=0)))
+
+
 def _get_required_packages():
     reqs = metadata.requires(__package__)
     # metadata.requires can return None if the package metadata cannot be found


Test patch:
diff --git a/scopesim/tests/tests_integrations/test_flux_is_conserved_through_full_system.py b/scopesim/tests/tests_integrations/test_flux_is_conserved_through_full_system.py
index 946ef2f0..3aea50f7 100644
--- a/scopesim/tests/tests_integrations/test_flux_is_conserved_through_full_system.py
+++ b/scopesim/tests/tests_integrations/test_flux_is_conserved_through_full_system.py
@@ -60,7 +60,7 @@ class TestObserve:
         print(src_flux, bg_flux)
         area = opt.optics_manager.area.value  # u.m**2
         assert src_flux == approx(1)          # u.Unit("ph s-1")
-        assert np.sum(im) == approx(src_flux * area, rel=2e-3)
+        assert np.sum(im) == approx(src_flux * area)
 
     def test_flux_is_conserved_for_yes_bg_emission(self, non_unity_cmds, tbl_src):
         """
@@ -80,8 +80,8 @@ class TestObserve:
 
         # given a 1 um bandpass
         assert src_flux == approx(1)          # u.Unit("ph s-1")
-        assert np.sum(im - np.median(im)) == approx(0.45, rel=1e-2)
-        assert np.median(im) == approx(1.5, abs=1e-2)
+        assert np.sum(im - np.median(im)) == approx(0.45)
+        assert np.median(im) == approx(1.5, abs=2e-3)
 
 
 @pytest.mark.usefixtures("protect_currsys", "patch_all_mock_paths")
diff --git a/scopesim/tests/tests_optics/test_FOVManager.py b/scopesim/tests/tests_optics/test_FOVManager.py
index b00502a3..a8e2c93b 100644
--- a/scopesim/tests/tests_optics/test_FOVManager.py
+++ b/scopesim/tests/tests_optics/test_FOVManager.py
@@ -1,6 +1,7 @@
 import pytest
 from pytest import approx
 import numpy as np
+from astropy import units as u
 
 from scopesim.optics.fov_manager import FOVManager
 from scopesim.tests.mocks.py_objects import effects_objects as eo
@@ -37,11 +38,11 @@ class TestGenerateFovList:
         effects = eo._mvs_effects_list()
         fov_man = FOVManager(effects=effects, pixel_scale=1, plate_scale=1)
         fovs = list(fov_man.generate_fovs_list())
-        fov_volume = fovs[0].volume()
+        fov_skycorners, _ = fovs[0].get_corners("deg")
 
         assert len(fovs) == 1
-        assert fov_volume["xs"][0] == approx(-1024 / 3600)      # [deg] 2k detector / pixel_scale
-        assert fov_volume["waves"][0] == 0.6            # [um] filter blue edge
+        assert fov_skycorners.min(axis=0)[0] == approx(-1024 / 3600)  # [deg] 2k detector / pixel_scale
+        assert fovs[0].waverange[0] == 0.6 * u.um  # filter blue edge
 
     @pytest.mark.parametrize("chunk_size, n_fovs",
                              [(500, 25), (512, 16), (1000, 9), (1024, 4)])
@@ -50,11 +51,11 @@ class TestGenerateFovList:
         fov_man = FOVManager(effects=effects, pixel_scale=1, plate_scale=1,
                              max_segment_size=1024**2, chunk_size=1024)
         fovs = list(fov_man.generate_fovs_list())
-        fov_volume = fovs[0].volume()
+        fov_skycorners, _ = fovs[0].get_corners("deg")
 
         assert len(fovs) == 4
-        assert fov_volume["xs"][0] == approx(-1024 / 3600)     # [deg] 2k detector / pixel_scale
-        assert fov_volume["waves"][0] == 0.6            # [um] filter blue edge
+        assert fov_skycorners.min(axis=0)[0] == approx(-1024 / 3600)  # [deg] 2k detector / pixel_scale
+        assert fovs[0].waverange[0] == 0.6 * u.um  # filter blue edge
 
     def test_fov_volumes_have_detector_dimensions_from_detector_list(self):
         effects = eo._mvs_effects_list()
diff --git a/scopesim/tests/tests_optics/test_FieldOfView.py b/scopesim/tests/tests_optics/test_FieldOfView.py
index 5ce12043..07f43d21 100644
--- a/scopesim/tests/tests_optics/test_FieldOfView.py
+++ b/scopesim/tests/tests_optics/test_FieldOfView.py
@@ -9,8 +9,7 @@ from matplotlib.colors import LogNorm
 
 from scopesim.tests.mocks.py_objects import header_objects as ho
 from scopesim.tests.mocks.py_objects import source_objects as so
-from scopesim.optics.fov import FieldOfView
-from scopesim.optics.fov_utils import get_cube_waveset
+from scopesim.optics.fov import FieldOfView, get_cube_waveset
 
 PLOTS = False
 
@@ -60,8 +59,8 @@ class TestExtractFrom:
     #                            "is extracted..."))
     def test_extract_point_sources_from_table(self):
         src = so._table_source()
-        src.fields[0]["x"] = [-15,-5,0,0] * u.arcsec
-        src.fields[0]["y"] = [0,0,5,15] * u.arcsec
+        src.fields[0]["x"] = [-15, -5, 0, 0] * u.arcsec
+        src.fields[0]["y"] = [0, 0, 5, 15] * u.arcsec
         fov = _fov_190_210_um()
         fov.extract_from(src)
 
@@ -94,15 +93,15 @@ class TestExtractFrom:
         fov = _fov_197_202_um()
         fov.extract_from(src)
 
-        assert fov.fields[0].shape == (3, 51, 25)
+        assert fov.fields[0].field.shape == (3, 51, 25)
 
     # @pytest.mark.xfail(reason=("is_field_in_fov drops table if anything is "
     #                            "outside fov volume, therefore no point source "
     #                            "is extracted..."))
     def test_extract_one_of_each_type_from_source_object(self):
         src_table = so._table_source()              # 4 sources, put two outside of FOV
-        src_table.fields[0]["x"] = [-15,-5,0,0] * u.arcsec
-        src_table.fields[0]["y"] = [0,0,5,15] * u.arcsec
+        src_table.fields[0]["x"] = [-15, -5, 0, 0] * u.arcsec
+        src_table.fields[0]["y"] = [0, 0, 5, 15] * u.arcsec
         src_image = so._image_source(dx=10)         # 10x10" @ 0.2"/pix
         src_cube = so._cube_source()                # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm
         src = src_cube + src_image + src_table
@@ -110,9 +109,9 @@ class TestExtractFrom:
         fov = _fov_197_202_um()
         fov.extract_from(src)
 
-        assert fov.fields[0].shape == (3, 51, 51)
-        assert fov.fields[1].shape == (51, 25)
-        assert len(fov.fields[2]) == 2
+        assert fov.fields[0].field.shape == (3, 51, 51)
+        assert fov.fields[1].field.shape == (51, 25)
+        assert len(fov.fields[2].field) == 2
 
         assert len(fov.spectra) == 3
         assert fov.fields[1].header["SPEC_REF"] == 0
@@ -147,10 +146,10 @@ class TestExtractFrom:
                               area=1 * u.m ** 2)
         the_fov.extract_from(src)
         assert len(the_fov.fields) == 3
-        assert isinstance(the_fov.fields[0], fits.ImageHDU)
-        assert isinstance(the_fov.fields[1], fits.ImageHDU)
+        assert isinstance(the_fov.fields[0].field, fits.ImageHDU)
+        assert isinstance(the_fov.fields[1].field, fits.ImageHDU)
         assert the_fov.fields[1].header["NAXIS"] == 3
-        assert isinstance(the_fov.fields[2], Table)
+        assert isinstance(the_fov.fields[2].field, Table)
 
     def test_handles_nans(self):
         src = so._image_source()
diff --git a/scopesim/tests/tests_optics/test_fov_utls.py b/scopesim/tests/tests_optics/test_fov_utls.py
index 2009718c..267e1f9d 100644
--- a/scopesim/tests/tests_optics/test_fov_utls.py
+++ b/scopesim/tests/tests_optics/test_fov_utls.py
@@ -8,7 +8,8 @@ from synphot.units import PHOTLAM
 from astropy import units as u
 from astropy.io import fits
 
-from scopesim.optics import FieldOfView, fov_utils
+from scopesim.optics.fov import FieldOfView, extract_range_from_spectrum
+from scopesim.optics.fov_utils import make_cube_from_table
 from scopesim.optics import image_plane_utils as imp_utils
 
 from scopesim.tests.mocks.py_objects import header_objects as ho
@@ -32,19 +33,19 @@ class TestExtractAreaFromImageHDU:
                                              basic_fov_header):
         fov = FieldOfView(basic_fov_header, [0.5, 2.5])
         field = cube_source.fields[0]
-        new_field = fov_utils.extract_area_from_imagehdu(field, fov.volume())
+        new_field = fov.extract_area_from_imagehdu(field, fov.get_corners("deg")[0])
 
         if PLOTS:
             xy = imp_utils.calc_footprint(basic_fov_header)
-            x, y = xy[:, 0], xy[:, 1]
-            plt.fill(x, y, c="r")
+            plt.fill(xy[:, 0], xy[:, 1], c="r", label="FOV")
             xy = imp_utils.calc_footprint(field.header)
-            x, y = xy[:, 0], xy[:, 1]
-            plt.fill(x, y, c="y")
+            plt.fill(xy[:, 0], xy[:, 1], c="y", label="Source")
             xy = imp_utils.calc_footprint(new_field.header)
-            x, y = xy[:, 0], xy[:, 1]
-            plt.fill(x, y, c="g")
+            plt.fill(xy[:, 0], xy[:, 1], c="g", label="Extracted")
 
+            plt.grid()
+            plt.legend()
+            plt.gca().set_aspect("equal")
             plt.show()
 
         assert new_field.header["NAXIS1"] == field.header["NAXIS1"]
@@ -55,19 +56,19 @@ class TestExtractAreaFromImageHDU:
                                                       basic_fov_header):
         fov = FieldOfView(basic_fov_header, [1.3, 1.7])
         field = cube_source.fields[0]
-        new_field = fov_utils.extract_area_from_imagehdu(field, fov.volume())
+        new_field = fov.extract_area_from_imagehdu(field, fov.get_corners("deg")[0])
 
         if PLOTS:
             xy = imp_utils.calc_footprint(basic_fov_header)
-            x, y = xy[:, 0], xy[:, 1]
-            plt.fill(x, y, c="r")
+            plt.fill(xy[:, 0], xy[:, 1], c="r", label="FOV")
             xy = imp_utils.calc_footprint(field.header)
-            x, y = xy[:, 0], xy[:, 1]
-            plt.fill(x, y, c="y")
+            plt.fill(xy[:, 0], xy[:, 1], c="y", label="Source")
             xy = imp_utils.calc_footprint(new_field.header)
-            x, y = xy[:, 0], xy[:, 1]
-            plt.fill(x, y, c="g")
+            plt.fill(xy[:, 0], xy[:, 1], c="g", label="Extracted")
 
+            plt.grid()
+            plt.legend()
+            plt.gca().set_aspect("equal")
             plt.show()
 
         assert new_field.header["NAXIS1"] == field.header["NAXIS1"]
@@ -75,29 +76,31 @@ class TestExtractAreaFromImageHDU:
         assert new_field.header["NAXIS3"] == 21
 
     def test_returns_eigth_cube_for_3d_offset_fov(self, cube_source,
-                                                         basic_fov_header):
+                                                  basic_fov_header):
         hdr = basic_fov_header
         hdr["CRVAL1"] += 75 * hdr["CDELT1"]
         hdr["CRVAL2"] += 75 * hdr["CDELT2"]
         fov = FieldOfView(hdr, [1.5, 2.5])
         field = cube_source.fields[0]
-        new_field = fov_utils.extract_area_from_imagehdu(field, fov.volume())
+        new_field = fov.extract_area_from_imagehdu(field, fov.get_corners("deg")[0])
 
         if PLOTS:
             xy = imp_utils.calc_footprint(basic_fov_header)
-            x, y = xy[:, 0], xy[:, 1]
-            plt.fill(x, y, c="r")
+            plt.fill(xy[:, 0], xy[:, 1], c="r", label="FOV")
             xy = imp_utils.calc_footprint(field.header)
-            x, y = xy[:, 0], xy[:, 1]
-            plt.fill(x, y, c="y")
+            plt.fill(xy[:, 0], xy[:, 1], c="y", label="Source")
             xy = imp_utils.calc_footprint(new_field.header)
-            x, y = xy[:, 0], xy[:, 1]
-            plt.fill(x, y, c="g")
+            plt.fill(xy[:, 0], xy[:, 1], c="g", label="Extracted")
 
+            plt.grid()
+            plt.legend()
+            plt.gca().set_aspect("equal")
             plt.show()
 
-        assert new_field.header["NAXIS1"] == 25
-        assert new_field.header["NAXIS2"] == 25
+        # Note: 26 is correct because there are actually 25.5 source pixels in
+        #       the FOV, but the cutout is "generous".
+        assert new_field.header["NAXIS1"] == 26
+        assert new_field.header["NAXIS2"] == 26
         assert new_field.header["NAXIS3"] == 51
 
 
@@ -108,38 +111,32 @@ class TestExtractRangeFromSpectrum:
         spec = SourceSpectrum(Empirical1D, points=wave, lookup_table=flux)
 
         waverange = [1.98, 2.12] * u.um
-        new_spec = fov_utils.extract_range_from_spectrum(spec, waverange)
+        new_spec = extract_range_from_spectrum(spec, waverange)
 
         assert len(new_spec.waverange) == 2
         assert new_spec.waverange[0] == 1.98 * u.um
         assert new_spec(1.98 * u.um).value == approx(12.8)
 
-    @pytest.mark.skip(reason="Kicking the can down the road")
-    def test_throws_error_if_no_overlap_between_waverange_and_waveset(self):
-        wave = np.arange(0.7, 1.5, 0.1) * u.um
+    @pytest.mark.parametrize(("endpoint", "msg"),
+                             [pytest.param(1.5, "Waverange does not overlap", marks=pytest.mark.xfail(reason="Check was disabled in function, dunno why.")),
+                              (2.05, "Waverange only partially overlaps")])
+    def test_logs_msg_for_waverang_overlap_mismatch(
+            self, endpoint, msg, caplog):
+        wave = np.arange(0.7, endpoint, 0.1) * u.um
         flux = np.arange(len(wave)) * PHOTLAM
         spec = SourceSpectrum(Empirical1D, points=wave, lookup_table=flux)
 
-        with pytest.raises(ValueError):
-            waverange = [1.98, 2.12] * u.um
-            new_spec = fov_utils.extract_range_from_spectrum(spec, waverange)
-
-    @pytest.mark.skip(reason="Kicking the can down the road")
-    def test_throws_error_if_only_partial_overlap_exists(self):
-        wave = np.arange(0.7, 2.05, 0.1) * u.um
-        flux = np.arange(len(wave)) * PHOTLAM
-        spec = SourceSpectrum(Empirical1D, points=wave, lookup_table=flux)
+        waverange = [1.98, 2.12] * u.um
+        extract_range_from_spectrum(spec, waverange)
 
-        with pytest.raises(ValueError):
-            waverange = [1.98, 2.12] * u.um
-            new_spec = fov_utils.extract_range_from_spectrum(spec, waverange)
+        assert msg in caplog.text
 
 
 class TestMakeCubeFromTable():
     def test_returns_an_imagehdu(self):
         src_table = so._table_source()
-        src_table.fields[0]["x"] = [-15,-5,0,0] * u.arcsec
-        src_table.fields[0]["y"] = [0,0,5,15] * u.arcsec
+        src_table.fields[0]["x"] = [-15, -5, 0, 0] * u.arcsec
+        src_table.fields[0]["y"] = [0, 0, 5, 15] * u.arcsec
 
         hdr = ho._fov_header()  # 20x20" @ 0.2" --> [-10, 10]"
         wav = [1.9, 2.1] * u.um
@@ -148,7 +145,7 @@ class TestMakeCubeFromTable():
         fov.extract_from(src_table)
 
         waveset = np.linspace(wav[0], wav[1], 51)
-        hdu = fov_utils.make_cube_from_table(fov.fields[0], fov.spectra,
-                                             waveset, fov.header)
+        hdu = make_cube_from_table(fov.fields[0], fov.spectra,
+                                   waveset, fov.header)
 
         assert isinstance(hdu, fits.ImageHDU)