Repository: AstarVienna/ScopeSim (language: python)

Issue:
Preparing for 0.6.1
Release notes in version.py and version number in pyproject.toml


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/notebooks_with_irdb_download.yml b/.github/workflows/notebooks_with_irdb_download.yml
index 26193abc..3b440c89 100644
--- a/.github/workflows/notebooks_with_irdb_download.yml
+++ b/.github/workflows/notebooks_with_irdb_download.yml
@@ -18,8 +18,15 @@ jobs:
         # that users have with ScopeSim / IRDB.
         # However, only use minimum and maximum supported Python version,
         # as the IRDB download often fails.
-        os: [ubuntu-latest, windows-latest, macos-latest]
-        python-version: ['3.8', '3.11']
+        # As of 2023/07/19, this test still often fails. See
+        # https://github.com/AstarVienna/ScopeSim/issues/254
+        # For now, only run on one configuration, as having one green test
+        # is more useful than 5 green and one red. The others can be enabled
+        # once we fix the problem properly.
+        # os: [ubuntu-latest, windows-latest, macos-latest]
+        # python-version: ['3.8', '3.11']
+        os: [ubuntu-latest]
+        python-version: ['3.11']
 
     steps:
       - uses: actions/checkout@v3
diff --git a/pyproject.toml b/pyproject.toml
index bed6f163..f05624ad 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "ScopeSim"
-version = "0.6.0"
+version = "0.6.1"
 description = "Generalised telescope observation simulator"
 readme = "README.md"
 requires-python = ">=3.8"
@@ -58,6 +58,7 @@ docs = [
     "sphinxcontrib-apidoc",
     "nbsphinx",
     "numpydoc",
+    "scopesim_templates>=0.4.4",
 ]
 
 [project.urls]
diff --git a/scopesim/commands/user_commands.py b/scopesim/commands/user_commands.py
index 793c3bbf..10dff1e9 100644
--- a/scopesim/commands/user_commands.py
+++ b/scopesim/commands/user_commands.py
@@ -231,6 +231,8 @@ def set_modes(self, modes=None):
                 for mode in modes:
                     if mode in self.modes_dict:
                         defyam["properties"]["modes"].append(mode)
+                        if "deprecate" in self.modes_dict[mode]:
+                            logging.warning(self.modes_dict[mode]["deprecate"])
                     else:
                         raise ValueError(f"mode '{mode}' was not recognised")
 
@@ -243,6 +245,8 @@ def list_modes(self):
                 dic = self.modes_dict[mode_name]
                 desc = dic["description"] if "description" in dic else "<None>"
                 modes[mode_name] = desc
+                if "deprecate" in dic:
+                    modes[mode_name] += " (deprecated)"
 
             msg = "\n".join([f"{key}: {value}" for key, value in modes.items()])
         else:
@@ -265,6 +269,15 @@ def __contains__(self, item):
     def __repr__(self):
         return f"{self.__class__.__name__}(**{self.kwargs!r})"
 
+    def __str__(self):
+        return str(self.cmds)
+
+    def _repr_pretty_(self, p, cycle):
+        """For ipython"""
+        if cycle:
+            p.text("UserCommands(...)")
+        else:
+            p.text(str(self))
 
 def check_for_updates(package_name):
     """
diff --git a/scopesim/detector/detector_array.py b/scopesim/detector/detector_array.py
index 10321eb9..b53df817 100644
--- a/scopesim/detector/detector_array.py
+++ b/scopesim/detector/detector_array.py
@@ -19,7 +19,7 @@ def __init__(self, detector_list=None, **kwargs):
         self.detectors = []
         self.latest_exposure = None
 
-    def readout(self, image_planes, array_effects=[], dtcr_effects=[], **kwargs):
+    def readout(self, image_planes, array_effects=None, dtcr_effects=None, **kwargs):
         """
         Read out the detector array into a FITS file
 
@@ -54,8 +54,8 @@ def readout(self, image_planes, array_effects=[], dtcr_effects=[], **kwargs):
         # - add ImageHDUs
         # - add ASCIITableHDU with Effects meta data in final table extension
 
-        self.array_effects = array_effects
-        self.dtcr_effects = dtcr_effects
+        self.array_effects = array_effects or []
+        self.dtcr_effects = dtcr_effects or []
         self.meta.update(kwargs)
 
         # 0. Get the image plane that corresponds to this detector array
@@ -101,6 +101,21 @@ def readout(self, image_planes, array_effects=[], dtcr_effects=[], **kwargs):
 
         return self.latest_exposure
 
+    def __repr__(self):
+        msg = (f"{self.__class__.__name__}"
+               f"({self.detector_list!r}, **{self.meta!r})")
+        return msg
+
+    def __str__(self):
+        return f"{self.__class__.__name__} with {self.detector_list!s}"
+
+    def _repr_pretty_(self, p, cycle):
+        """For ipython"""
+        if cycle:
+            p.text(f"{self.__class__.__name__}(...)")
+        else:
+            p.text(str(self))
+
 
 def make_primary_hdu(meta):
     """Create the primary header from meta data"""
diff --git a/scopesim/effects/apertures.py b/scopesim/effects/apertures.py
index 872d0894..4953b7d4 100644
--- a/scopesim/effects/apertures.py
+++ b/scopesim/effects/apertures.py
@@ -14,7 +14,8 @@
 from ..optics import image_plane_utils as imp_utils
 from ..base_classes import FOVSetupBase
 
-from ..utils import quantify, quantity_from_table, from_currsys, check_keys
+from ..utils import quantify, quantity_from_table, from_currsys, check_keys, \
+    figure_factory
 
 
 class ApertureMask(Effect):
@@ -83,7 +84,7 @@ def __init__(self, **kwargs):
                 w, h = kwargs["width"], kwargs["height"]
                 kwargs["filename"] = kwargs["filename_format"].format(w, h)
 
-        super(ApertureMask, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"pixel_scale": "!INST.pixel_scale",
                   "no_mask": True,
                   "angle": 0,
@@ -175,22 +176,23 @@ def get_mask(self):
 
         return mask
 
-    def plot(self, new_figure=True):
-        import matplotlib.pyplot as plt
-        if new_figure:
-            plt.gcf().clf()
+    def plot(self, axes=None):
+        if axes is None:
+            fig, ax = figure_factory()
+        else:
+            fig = axes.figure
 
         x = list(self.table["x"].data)
         y = list(self.table["y"].data)
-        plt.plot(x + [x[0]], y + [y[0]])
-        plt.gca().set_aspect("equal")
+        ax.plot(x + [x[0]], y + [y[0]])
+        ax.set_aspect("equal")
 
-        return plt.gcf()
+        return fig
 
 
 class RectangularApertureMask(ApertureMask):
     def __init__(self, **kwargs):
-        super(RectangularApertureMask, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"x_unit" : "arcsec",
                   "y_unit" : "arcsec"}
         self.meta.update(params)
@@ -268,7 +270,7 @@ class ApertureList(Effect):
 
     """
     def __init__(self, **kwargs):
-        super(ApertureList, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"pixel_scale": "!INST.pixel_scale",
                   "n_round_corners": 32,        # number of corners use to estimate ellipse
                   "no_mask": False,             # .. todo:: is this necessary when we have conserve_image?
@@ -342,26 +344,27 @@ def get_apertures(self, row_ids):
         return apertures_list
 
     def plot(self):
-        import matplotlib.pyplot as plt
-        plt.gcf().clf()
+        fig, ax = figure_factory()
 
         for ap in self.apertures:
-            ap.plot(new_figure=False)
+            ap.plot(ax)
 
-        return plt.gcf()
+        return fig
 
     def plot_masks(self):
-        import matplotlib.pyplot as plt
-
         aps = self.apertures
         n = len(aps)
         w = np.ceil(n ** 0.5).astype(int)
+        assert int(n ** 0.5) == w + 1
         h = np.ceil(n / w).astype(int)
+        assert int(n / w) == h + 1
+        # TODO: change these?
 
-        for ii, ap in enumerate(aps):
-            plt.subplot(w, h, ii + 1)
-            plt.imshow(ap.mask.T)
-        plt.show()
+        fig, axes = figure_factory(w, h)
+        for ap, ax in zip(aps, axes):
+            ax.imshow(ap.mask.T)
+        fig.show()
+        return fig
 
     def __add__(self, other):
         if isinstance(other, ApertureList):
@@ -419,11 +422,13 @@ class SlitWheel(Effect):
             current_slit: "C"
 
     """
-    def __init__(self, **kwargs):
-        required_keys = ["slit_names", "filename_format", "current_slit"]
-        check_keys(kwargs, required_keys, action="error")
 
-        super(SlitWheel, self).__init__(**kwargs)
+    required_keys = {"slit_names", "filename_format", "current_slit"}
+    _current_str = "current_slit"
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        check_keys(kwargs, self.required_keys, action="error")
 
         params = {"z_order": [80, 280, 580],
                   "path": "",
@@ -433,7 +438,7 @@ def __init__(self, **kwargs):
         self.meta.update(params)
         self.meta.update(kwargs)
 
-        path = Path(self.meta["path"], from_currsys(self.meta["filename_format"]))
+        path = self._get_path()
         self.slits = {}
         for name in from_currsys(self.meta["slit_names"]):
             kwargs["name"] = name
@@ -480,12 +485,6 @@ def current_slit(self):
             return False
         return self.slits[currslit]
 
-    @property
-    def display_name(self):
-        return f"{self.meta['name']} : " \
-               f"[{from_currsys(self.meta['current_slit'])}]"
-
-
     def __getattr__(self, item):
         return getattr(self.current_slit, item)
 
diff --git a/scopesim/effects/detector_list.py b/scopesim/effects/detector_list.py
index 405e21e6..347078b5 100644
--- a/scopesim/effects/detector_list.py
+++ b/scopesim/effects/detector_list.py
@@ -8,6 +8,7 @@
 from .effects import Effect
 from .apertures import ApertureMask
 from .. import utils
+from ..utils import close_loop, figure_factory
 from ..optics.image_plane_utils import header_from_list_of_xy, calc_footprint
 
 __all__ = ["DetectorList", "DetectorWindow"]
@@ -105,7 +106,7 @@ class : DetectorList
 
     """
     def __init__(self, **kwargs):
-        super(DetectorList, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [90, 290, 390, 490],
                   "pixel_scale": "!INST.pixel_scale",      # arcsec
                   "active_detectors": "all",
@@ -251,22 +252,21 @@ def detector_headers(self, ids=None):
 
         return hdrs
 
-    def plot(self):
-        import matplotlib.pyplot as plt
-        plt.gcf().clf()
+    def plot(self, axes=None):
+        if axes is None:
+            _, axes = figure_factory()
 
         for hdr in self.detector_headers():
             x_mm, y_mm = calc_footprint(hdr, "D")
-            x_cen, y_cen = np.average(x_mm), np.average(y_mm)
-            x_mm = list(x_mm) + [x_mm[0]]
-            y_mm = list(y_mm) + [y_mm[0]]
-            plt.gca().plot(x_mm, y_mm)
-            plt.gca().text(x_cen, y_cen, hdr["ID"])
+            axes.plot(list(close_loop(x_mm)), list(close_loop(y_mm)))
+            axes.text(*np.mean((x_mm, y_mm), axis=1), hdr["ID"],
+                      ha="center", va="center")
 
-        plt.gca().set_aspect("equal")
-        plt.ylabel("Size [mm]")
+        axes.set_aspect("equal")
+        axes.set_xlabel("Size [mm]")
+        axes.set_ylabel("Size [mm]")
 
-        return plt.gcf()
+        return axes
 
 
 class DetectorWindow(DetectorList):
@@ -314,4 +314,4 @@ def __init__(self, pixel_size, x, y, width, height=None, angle=0, gain=1,
                            "angle", "gain", "pixel_size"])
         tbl.meta.update(params)
 
-        super(DetectorWindow, self).__init__(table=tbl, **params)
+        super().__init__(table=tbl, **params)
diff --git a/scopesim/effects/effects.py b/scopesim/effects/effects.py
index 00dbe9e9..50ed67be 100644
--- a/scopesim/effects/effects.py
+++ b/scopesim/effects/effects.py
@@ -36,7 +36,7 @@ class Effect(DataContainer):
     """
 
     def __init__(self, **kwargs):
-        super(Effect, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = []
         self.meta["include"] = True
         self.meta.update(kwargs)
@@ -109,19 +109,19 @@ def include(self, item):
 
     @property
     def display_name(self):
-        return self.meta.get("name", self.meta.get("filename", "<empty>"))
+        name = self.meta.get("name", self.meta.get("filename", "<untitled>"))
+        if not hasattr(self, "_current_str"):
+            return name
+        return f"{name} : [{from_currsys(self.meta[self._current_str])}]"
 
     @property
     def meta_string(self):
-        meta_str = ""
-        max_key_len = max(len(key) for key in self.meta.keys())
-        padlen = max_key_len + 4
-        for key in self.meta:
-            if key not in {"comments", "changes", "description", "history",
-                           "report_table_caption", "report_plot_caption",
-                           "table"}:
-                meta_str += f"{key:>{padlen}} : {self.meta[key]}\n"
-
+        padlen = 4 + len(max(self.meta, key=len))
+        exclude = {"comments", "changes", "description", "history",
+                   "report_table_caption", "report_plot_caption", "table"}
+        meta_str = "\n".join(f"{key:>{padlen}} : {value}"
+                             for key, value in self.meta.items()
+                             if key not in exclude)
         return meta_str
 
     def report(self, filename=None, output="rst", rst_title_chars="*+",
@@ -199,7 +199,7 @@ def report(self, filename=None, output="rst", rst_title_chars="*+",
 
         """
         changes = self.meta.get("changes", [])
-        changes_str = "- " + "\n- ".join([str(entry) for entry in changes])
+        changes_str = "- " + "\n- ".join(str(entry) for entry in changes)
         cls_doc = self.__doc__ if self.__doc__ is not None else "<no docstring>"
         cls_descr = cls_doc.lstrip().splitlines()[0]
 
@@ -239,7 +239,12 @@ def report(self, filename=None, output="rst", rst_title_chars="*+",
 """
 
         if params["report_plot_include"] and hasattr(self, "plot"):
+            from matplotlib.figure import Figure
             fig = self.plot()
+            # HACK: plot methods should always return the same, while this is
+            #       not sorted out, deal with both fig and ax
+            if not isinstance(fig, Figure):
+                fig = fig.figure
 
             if fig is not None:
                 path = params["report_image_path"]
@@ -257,6 +262,9 @@ def report(self, filename=None, output="rst", rst_title_chars="*+",
                 #                            params["report_rst_path"])
                 # rel_file_path = os.path.join(rel_path, fname)
 
+                # TODO: fname is set in a loop above, so using it here in the
+                #       fstring will only access the last value from the loop,
+                #       is that intended?
                 rst_str += f"""
 .. figure:: {fname}
     :name: {"fig:" + params.get("name", "<unknown Effect>")}
@@ -287,16 +295,11 @@ def report(self, filename=None, output="rst", rst_title_chars="*+",
         return rst_str
 
     def info(self):
-        """
-        Prints basic information on the effect, notably the description
-        """
-        text = str(self)
-
-        desc = self.meta.get("description")
-        if desc is not None:
-            text += f"\nDescription: {desc}"
-
-        print(text)
+        """Print basic information on the effect, notably the description."""
+        if (desc := self.meta.get("description")) is not None:
+            print(f"{self}\nDescription: {desc}")
+        else:
+            print(self)
 
     def __repr__(self):
         return f"{self.__class__.__name__}(**{self.meta!r})"
@@ -321,3 +324,9 @@ def __getitem__(self, item):
             raise ValueError(f"__getitem__ calls must start with '#': {item}")
 
         return value
+
+    def _get_path(self):
+        if any(key not in self.meta for key in ("path", "filename_format")):
+            return None
+        return Path(self.meta["path"],
+                    from_currsys(self.meta["filename_format"]))
diff --git a/scopesim/effects/electronic.py b/scopesim/effects/electronic.py
index 59601575..1ecf7c9d 100644
--- a/scopesim/effects/electronic.py
+++ b/scopesim/effects/electronic.py
@@ -22,14 +22,13 @@
 import logging
 
 import numpy as np
-import matplotlib.pyplot as plt
 
 from astropy.io import fits
 
 from .. import rc
 from . import Effect
 from ..base_classes import DetectorBase, ImagePlaneBase
-from ..utils import from_currsys
+from ..utils import from_currsys, figure_factory
 from .. import utils
 
 
@@ -231,7 +230,7 @@ class SummedExposure(Effect):
 
     """
     def __init__(self, **kwargs):
-        super(SummedExposure, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [860]}
         self.meta.update(params)
         self.meta.update(kwargs)
@@ -255,7 +254,7 @@ class Bias(Effect):
 
     """
     def __init__(self, **kwargs):
-        super(Bias, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [855]}
         self.meta.update(params)
         self.meta.update(kwargs)
@@ -272,7 +271,7 @@ def apply_to(self, obj, **kwargs):
 
 class PoorMansHxRGReadoutNoise(Effect):
     def __init__(self, **kwargs):
-        super(PoorMansHxRGReadoutNoise, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [811],
                   "pedestal_fraction": 0.3,
                   "read_fraction": 0.4,
@@ -311,19 +310,21 @@ def apply_to(self, det, **kwargs):
 
         return det
 
-    def plot(self, det, new_figure=False, **kwargs):
+    def plot(self, det, **kwargs):
         dtcr = self.apply_to(det)
-        plt.imshow(dtcr.data, origin="lower")
+        fig, ax = figure_factory()
+        ax.imshow(dtcr.data, origin="lower")
 
     def plot_hist(self, det, **kwargs):
         dtcr = self.apply_to(det)
-        plt.hist(dtcr.data.flatten())
+        fig, ax = figure_factory()
+        ax.hist(dtcr.data.flatten())
 
 
 class BasicReadoutNoise(Effect):
     """Readout noise computed as: ron * sqrt(NDIT)"""
     def __init__(self, **kwargs):
-        super(BasicReadoutNoise, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = [811]
         self.meta["random_seed"] = "!SIM.random.seed"
         self.meta.update(kwargs)
@@ -347,16 +348,18 @@ def apply_to(self, det, **kwargs):
 
     def plot(self, det):
         dtcr = self.apply_to(det)
-        plt.imshow(dtcr.data)
+        fig, ax = figure_factory()
+        ax.imshow(dtcr.data)
 
     def plot_hist(self, det, **kwargs):
         dtcr = self.apply_to(det)
-        plt.hist(dtcr.data.flatten())
+        fig, ax = figure_factory()
+        ax.hist(dtcr.data.flatten())
 
 
 class ShotNoise(Effect):
     def __init__(self, **kwargs):
-        super(ShotNoise, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = [820]
         self.meta["random_seed"] = "!SIM.random.seed"
         self.meta.update(kwargs)
@@ -392,11 +395,13 @@ def apply_to(self, det, **kwargs):
 
     def plot(self, det):
         dtcr = self.apply_to(det)
-        plt.imshow(dtcr.data)
+        fig, ax = figure_factory()
+        ax.imshow(dtcr.data)
 
     def plot_hist(self, det, **kwargs):
         dtcr = self.apply_to(det)
-        plt.hist(dtcr.data.flatten())
+        fig, ax = figure_factory()
+        ax.hist(dtcr.data.flatten())
 
 
 class DarkCurrent(Effect):
@@ -404,7 +409,7 @@ class DarkCurrent(Effect):
     required: dit, ndit, value
     """
     def __init__(self, **kwargs):
-        super(DarkCurrent, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = [830]
 
         required_keys = ["value", "dit", "ndit"]
@@ -436,9 +441,10 @@ def plot(self, det, **kwargs):
         dtcr = self.apply_to(det)
         dark_level = dtcr.data[0, 0] / total_time  # just read one pixel
         levels = dark_level * times
-        plt.plot(times, levels, **kwargs)
-        plt.xlabel("time")
-        plt.ylabel("dark level")
+        fig, ax = figure_factory()
+        ax.plot(times, levels, **kwargs)
+        ax.set_xlabel("time")
+        ax.set_ylabel("dark level")
 
 
 class LinearityCurve(Effect):
@@ -496,22 +502,22 @@ def apply_to(self, obj, **kwargs):
         return obj
 
     def plot(self, **kwargs):
-        plt.gcf().clf()
+        fig, ax = figure_factory()
 
         ndit = from_currsys(self.meta["ndit"])
         incident = self.table["incident"] * ndit
         measured = self.table["measured"] * ndit
 
-        plt.loglog(incident, measured, **kwargs)
-        plt.xlabel("Incident [ph s$^-1$]")
-        plt.ylabel("Measured [e- s$^-1$]")
+        ax.loglog(incident, measured, **kwargs)
+        ax.set_xlabel("Incident [ph s$^-1$]")
+        ax.set_ylabel("Measured [e- s$^-1$]")
 
-        return plt.gcf()
+        return fig
 
 
 class ReferencePixelBorder(Effect):
     def __init__(self, **kwargs):
-        super(ReferencePixelBorder, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = [780]
         val = 0
         if "all" in kwargs:
@@ -536,13 +542,14 @@ def apply_to(self, implane, **kwargs):
 
     def plot(self, implane, **kwargs):
         implane = self.apply_to(implane)
-        plt.imshow(implane.data, origin="bottom", **kwargs)
-        plt.show()
+        fig, ax = figure_factory()
+        ax.imshow(implane.data, origin="bottom", **kwargs)
+        # fig.show()
 
 
 class BinnedImage(Effect):
     def __init__(self, **kwargs):
-        super(BinnedImage, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = [870]
 
         self.required_keys = ["bin_size"]
@@ -560,7 +567,7 @@ def apply_to(self, det, **kwargs):
 
 class UnequalBinnedImage(Effect):
     def __init__(self, **kwargs):
-        super(UnequalBinnedImage, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = [870]
 
         self.required_keys = ["binx","biny"]
diff --git a/scopesim/effects/fits_headers.py b/scopesim/effects/fits_headers.py
index db858b20..e92d0a2a 100644
--- a/scopesim/effects/fits_headers.py
+++ b/scopesim/effects/fits_headers.py
@@ -216,7 +216,7 @@ class ExtraFitsKeywords(Effect):
     """
     def __init__(self, **kwargs):
         # don't pass kwargs, as DataContainer can't handle yaml files
-        super(ExtraFitsKeywords, self).__init__()
+        super().__init__()
         params = {"name": "extra_fits_keywords",
                   "description": "Extra FITS headers",
                   "z_order": [999],
@@ -326,7 +326,7 @@ def flatten_dict(dic, base_key="", flat_dict=None,
     if flat_dict is None:
         flat_dict = {}
     for key, val in dic.items():
-        flat_key = base_key + f"{key} "
+        flat_key = f"{base_key}{key} "
         if isinstance(val, dict):
             flatten_dict(val, flat_key, flat_dict, resolve, optics_manager)
         else:
@@ -340,17 +340,17 @@ def flatten_dict(dic, base_key="", flat_dict=None,
                 value = deepcopy(val)
 
             # resolve any bang or hash strings
-            if resolve and isinstance(value, str) and len(value) > 0:
-                if value[0] == "!":
+            if resolve and isinstance(value, str):
+                if value.startswith("!"):
                     value = from_currsys(value)
-                elif value[0] == "#":
+                elif value.startswith("#"):
                     if optics_manager is None:
                         raise ValueError("An OpticsManager object must be "
                                          "passed in order to resolve #-strings")
                     value = optics_manager[value]
 
             if isinstance(value, u.Quantity):
-                comment = f"[{str(value.unit)}] " + comment
+                comment = f"[{str(value.unit)}] {comment}"
                 value = value.value
 
             # Convert e.g.  Unit(mag) to just "mag". Not sure how this will
@@ -362,13 +362,13 @@ def flatten_dict(dic, base_key="", flat_dict=None,
                 value = f"{value.__class__.__name__}:{str(list(value))}"
                 max_len = 80 - len(flat_key)
                 if len(value) > max_len:
-                    value = value[:max_len-4] + " ..."
+                    value = f"{value[:max_len-4]} ..."
 
             if isinstance(value, (datetime.time, datetime.date, datetime.datetime)):
                 value = value.isoformat()
 
             # Add the flattened KEYWORD = (value, comment) to the header dict
-            if len(comment) > 0:
+            if comment:
                 flat_dict[flat_key] = (value, str(comment))
             else:
                 flat_dict[flat_key] = value
@@ -449,13 +449,18 @@ def apply_to(self, hdul, **kwargs):
                 # add effect under the EFFn keyword
                 prefix = self.meta["keyword_prefix"]
                 class_name = opt_train[eff_name].__class__.__name__
-                self.dict_list = [{"ext_number": self.meta["ext_number"],
-                                   "keywords": {
-                                       f"{prefix} EFF{i} class": [class_name, "ScopeSim class name"],
-                                       f"{prefix} EFF{i}": eff_meta}
-                                   }]
-                super_apply_to = super(EffectsMetaKeywords, self).apply_to
-                hdul = super_apply_to(hdul=hdul, optical_train=opt_train)
+                self.dict_list = [
+                    {"ext_number": self.meta["ext_number"],
+                     "keywords": {
+                         f"{prefix} EFF{i} class": [
+                             class_name,
+                             "ScopeSim class name"
+                             ],
+                         f"{prefix} EFF{i}": eff_meta
+                         }
+                     }
+                    ]
+                hdul = super().apply_to(hdul=hdul, optical_train=opt_train)
 
         return hdul
 
@@ -592,7 +597,6 @@ def apply_to(self, hdul, **kwargs):
                                f"{resolve_prefix}keywords": {
                                    f"{sim_prefix} CONFIG": cmds}
                                }]
-            super_apply_to = super(SimulationConfigFitsKeywords, self).apply_to
-            hdul = super_apply_to(hdul=hdul, optical_train=opt_train)
+            hdul = super().apply_to(hdul=hdul, optical_train=opt_train)
 
         return hdul
diff --git a/scopesim/effects/metis_lms_trace_list.py b/scopesim/effects/metis_lms_trace_list.py
index 040b6967..7a101dc7 100644
--- a/scopesim/effects/metis_lms_trace_list.py
+++ b/scopesim/effects/metis_lms_trace_list.py
@@ -217,7 +217,6 @@ def fov_grid(self):
         x_max = aperture["right"]
         y_min = aperture["bottom"]
         y_max = aperture["top"]
-        trace_id = self.meta["trace_id"]
 
         layout = ioascii.read(find_file("!DET.layout.file_name"))
         det_lims = {}
@@ -237,7 +236,7 @@ def fov_grid(self):
                 "y_min": y_min, "y_max": y_max,
                 "xi_min": xi_min, "xi_max": xi_max,
                 "wave_min": wave_min, "wave_max": wave_max,
-                "trace_id": trace_id}
+                "trace_id": self.trace_id}
 
     def get_waverange(self, det_mm_lims):
         """Determine wavelength range that spectral trace covers on image plane"""
diff --git a/scopesim/effects/obs_strategies.py b/scopesim/effects/obs_strategies.py
index 3b7fbd16..97fda51f 100644
--- a/scopesim/effects/obs_strategies.py
+++ b/scopesim/effects/obs_strategies.py
@@ -3,7 +3,6 @@
 - ChopNodCombiner: simulate chop-nod cycle
 """
 import numpy as np
-from matplotlib import pyplot as plt
 
 from scopesim.base_classes import DetectorBase
 from scopesim.effects import Effect
diff --git a/scopesim/effects/psf_utils.py b/scopesim/effects/psf_utils.py
index abe63773..dfe1324b 100644
--- a/scopesim/effects/psf_utils.py
+++ b/scopesim/effects/psf_utils.py
@@ -5,11 +5,11 @@
 from astropy import units as u
 from astropy.convolution import Gaussian2DKernel
 from astropy.io import fits
-import matplotlib.pyplot as plt
 from matplotlib.colors import LogNorm
 
 from .. import rc, utils
 from ..optics import image_plane_utils as imp_utils
+from ..utils import figure_factory
 
 
 def round_kernel_edges(kernel):
@@ -74,9 +74,10 @@ def nmrms_from_strehl_and_wavelength(strehl, wavelength, strehl_hdu,
         nm = np.interp(strehl, strehls[::-1], nms[::-1])
 
     if plot:
-        plt.plot(nms, strehls)
-        plt.plot(nm, strehl, "ro")
-        plt.show()
+        fig, ax = figure_factory()
+        ax.plot(nms, strehls)
+        ax.plot(nm, strehl, "ro")
+        fig.show()
 
     return nm
 
diff --git a/scopesim/effects/psfs.py b/scopesim/effects/psfs.py
index fe26049f..ca788c24 100644
--- a/scopesim/effects/psfs.py
+++ b/scopesim/effects/psfs.py
@@ -13,6 +13,9 @@
 from . import psf_utils as pu
 from ..base_classes import ImagePlaneBase, FieldOfViewBase, FOVSetupBase
 from .. import utils
+from ..utils import figure_grid_factory
+
+# TODO: directly import currsys stuff, replace utils.
 
 
 class PoorMansFOV:
@@ -133,14 +136,13 @@ def get_kernel(self, obj):
         return self.kernel
 
     def plot(self, obj=None, **kwargs):
-        import matplotlib.pyplot as plt
         from matplotlib.colors import LogNorm
-        plt.gcf().clf()
+        fig, ax = figure_factory()
 
         kernel = self.get_kernel(obj)
-        plt.imshow(kernel, norm=LogNorm(), origin="lower", **kwargs)
+        ax.imshow(kernel, norm=LogNorm(), origin="lower", **kwargs)
 
-        return plt.gcf()
+        return fig
 
 
 
@@ -247,8 +249,7 @@ def total_wfe(self):
         return self._total_wfe
 
     def plot(self):
-        import matplotlib.pyplot as plt
-        plt.gcf().clf()
+        fig, ax = figure_factory()
 
         wave_min, wave_max = utils.from_currsys([self.meta["wave_min"],
                                                  self.meta["wave_max"]])
@@ -256,11 +257,11 @@ def plot(self):
         wfe = self.total_wfe
         strehl = pu.wfe2strehl(wfe=wfe, wave=waves)
 
-        plt.plot(waves, strehl)
-        plt.xlabel(f"Wavelength [{waves.unit}]")
-        plt.ylabel(f"Strehl Ratio \n[Total WFE = {wfe}]")
+        ax.plot(waves, strehl)
+        ax.set_xlabel(f"Wavelength [{waves.unit}]")
+        ax.set_ylabel(f"Strehl Ratio \n[Total WFE = {wfe}]")
 
-        return plt.gcf()
+        return fig
 
 
 class SeeingPSF(AnalyticalPSF):
@@ -427,7 +428,7 @@ class AnisocadoConstPSF(SemiAnalyticalPSF):
 
     """
     def __init__(self, **kwargs):
-        super(AnisocadoConstPSF, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [42, 652],
                   "psf_side_length": 512,
                   "offset": (0, 0),
@@ -507,51 +508,58 @@ def nmRms(self):
         return nm_rms
 
     def plot(self, obj=None, **kwargs):
-        import matplotlib.pyplot as plt
         from matplotlib.colors import LogNorm
-        plt.figure(figsize=(10, 10))
+        fig, gs = figure_grid_factory(2, 2, height_ratios=(3, 2),
+            left=0.3, right=0.7, bottom=0.15, top=0.85,
+            wspace=0.05, hspace=0.05)
+        # or no height_ratios and bottom=0.1, top=0.9
 
         pixel_scale = utils.from_currsys("!INST.pixel_scale")
-
         kernel = self.get_kernel(pixel_scale)
-        # subplot2grid((n_vert, n_horiz), (from_top, from_left), colspan=1, rowspan=1)
-        plt.subplot2grid((2, 2), (0, 0))
+
+        ax = fig.add_subplot(gs[0, 0])
         im = kernel
         r_sky = pixel_scale * im.shape[0]
-        plt.imshow(im, norm=LogNorm(), origin="lower",
-                   extent= [-r_sky, r_sky, -r_sky, r_sky], **kwargs)
-        plt.ylabel("[arcsec]")
-
-        plt.subplot2grid((2, 2), (0, 1))
+        ax.imshow(im, norm=LogNorm(), origin="lower",
+                  extent= [-r_sky, r_sky, -r_sky, r_sky], **kwargs)
+        ax.set_aspect("equal")
+        ax.set_xlabel("[arcsec]")
+        ax.set_ylabel("[arcsec]")
+        ax.xaxis.set_ticks_position("top")
+        ax.xaxis.set_label_position("top")
+
+        ax = fig.add_subplot(gs[0, 1])
         x = kernel.shape[1] // 2
         y = kernel.shape[0] // 2
         r = 16
         im = kernel[y-r:y+r, x-r:x+r]
         r_sky = pixel_scale * im.shape[0]
-        plt.imshow(im, norm=LogNorm(), origin="lower",
-                   extent= [-r_sky, r_sky, -r_sky, r_sky], **kwargs)
-        plt.ylabel("[arcsec]")
-        plt.gca().yaxis.set_label_position("right")
-
-        plt.subplot2grid((2, 2), (1, 0), colspan=2)
+        ax.imshow(im, norm=LogNorm(), origin="lower",
+                  extent= [-r_sky, r_sky, -r_sky, r_sky], **kwargs)
+        ax.set_aspect("equal")
+        ax.set_xlabel("[arcsec]")
+        ax.set_ylabel("[arcsec]")
+        ax.xaxis.set_ticks_position("top")
+        ax.xaxis.set_label_position("top")
+        ax.yaxis.set_ticks_position("right")
+        ax.yaxis.set_label_position("right")
+
+        ax = fig.add_subplot(gs[1, :])
         hdr = self._file[0].header
         data = self._file[0].data
-
-        hdr = self._file[0].header
-        data = self._file[0].data
-
         wfes = np.arange(hdr["NAXIS1"]) * hdr["CDELT1"] + hdr["CRVAL1"]
         waves = np.arange(hdr["NAXIS2"]) * hdr["CDELT2"] + hdr["CRVAL2"]
-        for i in np.arange(len(waves))[::-1]:
-            plt.plot(wfes, data[i, :],
-                     label=f"{waves[i]:.3f} " + r"$\mu m$")
-
-        plt.xlabel("RMS Wavefront Error [um]")
-        plt.ylabel("Strehl Ratio")
-        plt.legend()
 
-        return plt.gcf()
+        # TODO: Get unit dynamically? Then again, it's hardcoded elsewhere in this module...
+        unit_str = u.Unit("um").to_string("latex")
+        for strehl, wav in reversed(list(zip(data, waves))):
+            ax.plot(wfes, strehl, label=f"{wav:.3f} {unit_str}")
 
+        ax.set_xlabel(f"RMS Wavefront Error [{unit_str}]")
+        ax.set_ylabel("Strehl Ratio")
+        ax.legend()
+        fig.align_labels()
+        return fig
 
 
 ################################################################################
@@ -701,7 +709,7 @@ class FieldVaryingPSF(DiscretePSF):
     """
     def __init__(self, **kwargs):
         # sub_pixel_flag and flux_accuracy are taken care of in PSF base class
-        super(FieldVaryingPSF, self).__init__(**kwargs)
+        super().__init__(**kwargs)
 
         self.required_keys = ["filename"]
         utils.check_keys(self.meta, self.required_keys, action="error")
diff --git a/scopesim/effects/rotation.py b/scopesim/effects/rotation.py
index 28ac85ef..9c347b99 100644
--- a/scopesim/effects/rotation.py
+++ b/scopesim/effects/rotation.py
@@ -22,7 +22,7 @@ class Rotate90CCD(Effect):
 
     """
     def __init__(self, **kwargs):
-        super(Rotate90CCD, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [809]}
         self.meta.update(params)
         self.meta.update(kwargs)
diff --git a/scopesim/effects/shifts.py b/scopesim/effects/shifts.py
index 84db3961..99e0c3d3 100644
--- a/scopesim/effects/shifts.py
+++ b/scopesim/effects/shifts.py
@@ -3,13 +3,14 @@
 from astropy.table import Table
 
 from .effects import Effect
-from ..utils import airmass2zendist, from_currsys, check_keys, quantify
+from ..utils import airmass2zendist, from_currsys, check_keys, quantify, \
+    figure_factory
 from ..base_classes import FieldOfViewBase
 
 
 class Shift3D(Effect):
     def __init__(self, **kwargs):
-        super(Shift3D, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [30, 230],
                   "report_plot_include": True,
                   "report_table_include": False,}
@@ -38,19 +39,18 @@ def get_table(self, **kwargs):
         return tbl
 
     def plot(self):
-        import matplotlib.pyplot as plt
-        plt.gcf().clf()
+        fig, ax = figure_factory()
 
         tbl = self.get_table()
-        plt.scatter(x=tbl["dx"], y=tbl["dy"], c=tbl["wavelength"])
-        plt.colorbar()
-        plt.xlabel(f"dx [{quantify(tbl['dx'], u.arcsec).unit}]")
-        plt.ylabel(f"dy [{quantify(tbl['dy'], u.arcsec).unit}]")
-        plt.axvline(0, ls=":")
-        plt.axhline(0, ls=":")
-        # plt.gca().set_aspect("equal")
+        pnts = ax.scatter(x=tbl["dx"], y=tbl["dy"], c=tbl["wavelength"])
+        fig.colorbar(pnts)
+        ax.set_xlabel(f"dx [{quantify(tbl['dx'], u.arcsec).unit}]")
+        ax.set_ylabel(f"dy [{quantify(tbl['dy'], u.arcsec).unit}]")
+        ax.axvline(0, ls=":")
+        ax.axhline(0, ls=":")
+        # ax.set_aspect("equal")
 
-        return plt.gcf()
+        return fig
 
 
 class AtmosphericDispersion(Shift3D):
@@ -93,7 +93,7 @@ class AtmosphericDispersion(Shift3D):
 
     """
     def __init__(self, **kwargs):
-        super(AtmosphericDispersion, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [231],
                   "wave_min": "!SIM.spectral.wave_min",
                   "wave_mid": "!SIM.spectral.wave_mid",
@@ -161,7 +161,7 @@ def __init__(self, **kwargs):
         ----------
         kwargs
         """
-        super(AtmosphericDispersionCorrection, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = [632]
         if "quick_adc" in self.meta and self.meta["quick_adc"] is True:
             self.meta["z_order"] += [232]
diff --git a/scopesim/effects/spectral_efficiency.py b/scopesim/effects/spectral_efficiency.py
index 6b255eec..c81f73d6 100644
--- a/scopesim/effects/spectral_efficiency.py
+++ b/scopesim/effects/spectral_efficiency.py
@@ -3,7 +3,6 @@
 """
 import logging
 import numpy as np
-from matplotlib import pyplot as plt
 
 from astropy.io import fits
 from astropy import units as u
@@ -13,7 +12,7 @@
 from .effects import Effect
 from .ter_curves import TERCurve
 from .ter_curves_utils import apply_throughput_to_cube
-from ..utils import find_file
+from ..utils import find_file, figure_factory
 from ..base_classes import FieldOfViewBase, FOVSetupBase
 
 class SpectralEfficiency(Effect):
@@ -101,11 +100,11 @@ def apply_to(self, obj, **kwargs):
         """
         Interface between FieldOfView and SpectralEfficiency
         """
-        trace_id = obj.meta['trace_id']
+        trace_id = obj.trace_id
         try:
             effic = self.efficiencies[trace_id]
         except KeyError:
-            logging.warning("No grating efficiency for trace %s" % trace_id)
+            logging.warning("No grating efficiency for trace %s", trace_id)
             return obj
 
         wcs = WCS(obj.hdu.header).spectral
@@ -116,12 +115,14 @@ def apply_to(self, obj, **kwargs):
 
     def plot(self):
         """Plot the grating efficiencies"""
+        fig, ax = figure_factory()
         for name, effic in self.efficiencies.items():
             wave = effic.throughput.waveset
-            plt.plot(wave.to(u.um), effic.throughput(wave), label=name)
+            ax.plot(wave.to(u.um), effic.throughput(wave), label=name)
 
-        plt.xlabel("Wavelength [um]")
-        plt.ylabel("Grating efficiency")
-        plt.title(f"Grating efficiencies from {self.filename}")
-        plt.legend()
-        plt.show()
+        ax.set_xlabel("Wavelength [um]")
+        ax.set_ylabel("Grating efficiency")
+        ax.set_title(f"Grating efficiencies {self.display_name}")
+        ax.legend()
+
+        return fig
diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py
index 553c7879..456fa043 100644
--- a/scopesim/effects/spectral_trace_list.py
+++ b/scopesim/effects/spectral_trace_list.py
@@ -7,6 +7,7 @@
 
 from pathlib import Path
 import logging
+from itertools import cycle
 
 import numpy as np
 
@@ -16,7 +17,7 @@
 from .effects import Effect
 from .ter_curves import FilterCurve
 from .spectral_trace_list_utils import SpectralTrace, make_image_interpolations
-from ..utils import from_currsys, check_keys
+from ..utils import from_currsys, check_keys, figure_factory
 from ..optics.image_plane_utils import header_from_list_of_xy
 from ..base_classes import FieldOfViewBase, FOVSetupBase
 
@@ -179,7 +180,7 @@ def apply_to(self, obj, **kwargs):
                     ex_vol["meta"].update(vol)
                     ex_vol["meta"].pop("wave_min")
                     ex_vol["meta"].pop("wave_max")
-                new_vols_list += extracted_vols
+                new_vols_list.extend(extracted_vols)
 
             obj.volumes = new_vols_list
 
@@ -195,8 +196,7 @@ def apply_to(self, obj, **kwargs):
                 logging.info("Making cube")
                 obj.cube = obj.make_cube_hdu()
 
-            trace_id = obj.meta["trace_id"]
-            spt = self.spectral_traces[trace_id]
+            spt = self.spectral_traces[obj.trace_id]
             obj.hdu = spt.map_spectra_to_focal_plane(obj)
 
         return obj
@@ -208,11 +208,11 @@ def footprint(self):
         xfoot, yfoot = [], []
         for spt in self.spectral_traces.values():
             xtrace, ytrace = spt.footprint()
-            xfoot += xtrace
-            yfoot += ytrace
+            xfoot.extend(xtrace)
+            yfoot.extend(ytrace)
 
-        xfoot = [np.min(xfoot), np.max(xfoot), np.max(xfoot), np.min(xfoot)]
-        yfoot = [np.min(yfoot), np.min(yfoot), np.max(yfoot), np.max(yfoot)]
+        xfoot = [min(xfoot), max(xfoot), max(xfoot), min(xfoot)]
+        yfoot = [min(yfoot), min(yfoot), max(yfoot), max(yfoot)]
 
         return xfoot, yfoot
 
@@ -299,7 +299,7 @@ def rectify_traces(self, hdulist, xi_min=None, xi_max=None, interps=None,
         #pdu.header['FILTER'] = from_currsys("!OBS.filter_name_fw1")
         outhdul = fits.HDUList([pdu])
 
-        for i, trace_id in enumerate(self.spectral_traces):
+        for i, trace_id in enumerate(self.spectral_traces, start=1):
             hdu = self[trace_id].rectify(hdulist,
                                          interps=interps,
                                          bin_width=bin_width,
@@ -307,34 +307,53 @@ def rectify_traces(self, hdulist, xi_min=None, xi_max=None, interps=None,
                                          wave_min=wave_min, wave_max=wave_max)
             if hdu is not None:   # ..todo: rectify does not do that yet
                 outhdul.append(hdu)
-                outhdul[0].header[f"EXTNAME{i+1}"] = trace_id
+                outhdul[0].header[f"EXTNAME{i}"] = trace_id
 
         outhdul[0].header.update(inhdul[0].header)
 
         return outhdul
 
-
     def rectify_cube(self, hdulist):
         """Rectify traces and combine into a cube"""
         raise(NotImplementedError)
 
-    def plot(self, wave_min=None, wave_max=None, **kwargs):
+    def plot(self, wave_min=None, wave_max=None, axes=None, **kwargs):
+        """Plot every spectral trace in the spectral trace list.
+
+        Parameters
+        ----------
+        wave_min : float, optional
+            Minimum wavelength, if any. If None, value from_currsys is used.
+        wave_max : float, optional
+            Maximum wavelength, if any. If None, value from_currsys is used.
+        axes : matplotlib axes, optional
+            The axes object to use for the plot. If None (default), a new
+            figure with one axes will be created.
+        **kwargs : dict
+            Any other parameters passed along to the plot method of the
+            individual spectral traces.
+
+        Returns
+        -------
+        fig : matplotlib figure
+            DESCRIPTION.
+
+        """
         if wave_min is None:
             wave_min = from_currsys("!SIM.spectral.wave_min")
         if wave_max is None:
             wave_max = from_currsys("!SIM.spectral.wave_max")
 
-        from matplotlib import pyplot as plt
-        from matplotlib._pylab_helpers import Gcf
-        if len(Gcf.figs) == 0:
-            plt.figure(figsize=(12, 12))
+        if axes is None:
+            fig, axes = figure_factory()
+        else:
+            fig = axes.figure
 
         if self.spectral_traces is not None:
-            clrs = "rgbcymk" * (1 + len(self.spectral_traces) // 7)
-            for spt, c in zip(self.spectral_traces.values(), clrs):
-                spt.plot(wave_min, wave_max, c=c)
+            for spt, c in zip(self.spectral_traces.values(), cycle("rgbcymk")):
+                spt.plot(wave_min, wave_max, c=c, axes=axes, **kwargs)
 
-        return plt.gcf()
+        return fig
 
     def __repr__(self):
         # "\n".join([spt.__repr__() for spt in self.spectral_traces])
@@ -409,13 +428,14 @@ class SpectralTraceListWheel(Effect):
               filter_names: "!INST.grism_names"
 
     """
-    def __init__(self, **kwargs):
-        required_keys = ["trace_list_names",
-                         "filename_format",
-                         "current_trace_list"]
-        check_keys(kwargs, required_keys, action="error")
 
-        super(SpectralTraceListWheel, self).__init__(**kwargs)
+    required_keys = {"trace_list_names", "filename_format",
+                     "current_trace_list"}
+    _current_str = "current_trace_list"
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        check_keys(kwargs, self.required_keys, action="error")
 
         params = {"z_order": [70, 270, 670],
                   "path": "",
@@ -425,7 +445,7 @@ def __init__(self, **kwargs):
         self.meta.update(params)
         self.meta.update(kwargs)
 
-        path = Path(self.meta["path"], from_currsys(self.meta["filename_format"]))
+        path = self._get_path()
         self.trace_lists = {}
         for name in from_currsys(self.meta["trace_list_names"]):
             kwargs["name"] = name
@@ -443,8 +463,3 @@ def current_trace_list(self):
         if trace_list_name is not None:
             trace_list_eff = self.trace_lists[trace_list_name]
         return trace_list_eff
-
-    @property
-    def display_name(self):
-        name = self.meta.get("name", self.meta.get("filename", "<untitled>"))
-        return f"{name} : [{from_currsys(self.meta['current_trace_list'])}]"
diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py
index d43df124..da1241b3 100644
--- a/scopesim/effects/spectral_trace_list_utils.py
+++ b/scopesim/effects/spectral_trace_list_utils.py
@@ -15,16 +15,16 @@
 
 from scipy.interpolate import RectBivariateSpline
 from scipy.interpolate import interp1d
-from matplotlib import pyplot as plt
 
-from astropy.table import Table
+from astropy.table import Table, vstack
 from astropy.modeling import fitting
 from astropy.io import fits
 from astropy import units as u
 from astropy.wcs import WCS
 from astropy.modeling.models import Polynomial2D
 
-from ..utils import power_vector, quantify, from_currsys
+from ..utils import power_vector, quantify, from_currsys, close_loop, \
+    figure_factory
 
 
 class SpectralTrace:
@@ -76,6 +76,11 @@ def __init__(self, trace_tbl, **kwargs):
         self._xilamimg = None
         self.dlam_per_pix = None
 
+    @property
+    def trace_id(self):
+        """Return the name of the trace"""
+        return self.meta["trace_id"]
+
     def fov_grid(self):
         """
         Provide information on the source space volume required by the effect
@@ -86,7 +91,6 @@ def fov_grid(self):
         Spatial limits are determined by the `ApertureMask` effect
         and are not returned here.
         """
-        trace_id = self.meta["trace_id"]
         aperture_id = self.meta["aperture_id"]
         lam_arr = self.table[self.meta["wave_colname"]]
 
@@ -94,7 +98,7 @@ def fov_grid(self):
         wave_min = np.min(lam_arr)
 
         return {"wave_min": wave_min, "wave_max": wave_max,
-                "trace_id": trace_id, "aperture_id": aperture_id}
+                "trace_id": self.trace_id, "aperture_id": aperture_id}
 
     def compute_interpolation_functions(self):
         """
@@ -118,7 +122,7 @@ def compute_interpolation_functions(self):
         self._xiy2x = Transform2D.fit(xi_arr, y_arr, x_arr)
         self._xiy2lam = Transform2D.fit(xi_arr, y_arr, lam_arr)
 
-        if self.dispersion_axis == 'unknown':
+        if self.dispersion_axis == "unknown":
             dlam_dx, dlam_dy = self.xy2lam.gradient()
             wave_mid = 0.5 * (self.wave_min + self.wave_max)
             xi_mid = np.mean(xi_arr)
@@ -142,7 +146,7 @@ def map_spectra_to_focal_plane(self, fov):
         The method returns a section of the fov image along with info on
         where this image lies in the focal plane.
         """
-        logging.info("Mapping %s", fov.meta['trace_id'])
+        logging.info("Mapping %s", fov.trace_id)
         # Initialise the image based on the footprint of the spectral
         # trace and the focal plane WCS
         wave_min = fov.meta["wave_min"].value       # [um]
@@ -176,7 +180,7 @@ def map_spectra_to_focal_plane(self, fov):
         ## Check if spectral trace footprint is outside FoV
         if xmax < 0 or xmin > naxis1d or ymax < 0 or ymin > naxis2d:
             logging.info("Spectral trace %s: footprint is outside FoV",
-                         fov.meta['trace_id'])
+                         fov.trace_id)
             return None
 
         # Only work on parts within the FoV
@@ -207,7 +211,7 @@ def map_spectra_to_focal_plane(self, fov):
             xilam = XiLamImage(fov, self.dlam_per_pix)
             self._xilamimg = xilam   # ..todo: remove or make available with a debug flag?
         except ValueError:
-            print(f" ---> {self.meta['trace_id']} gave ValueError")
+            print(f" ---> {self.trace_id} gave ValueError")
 
         npix_xi, npix_lam = xilam.npix_xi, xilam.npix_lam
         xilam_wcs = xilam.wcs
@@ -321,7 +325,7 @@ def rectify(self, hdulist, interps=None, wcs=None, **kwargs):
             bin_width = np.abs(self.dlam_per_pix.y).min()
         logging.info("   Bin width %.02g um", bin_width)
 
-        pixscale = from_currsys(self.meta['pixel_scale'])
+        pixscale = from_currsys(self.meta["pixel_scale"])
 
         # Temporary solution to get slit length
         xi_min = kwargs.get("xi_min", None)
@@ -341,8 +345,8 @@ def rectify(self, hdulist, interps=None, wcs=None, **kwargs):
 
         if wcs is None:
             wcs = WCS(naxis=2)
-            wcs.wcs.ctype = ['WAVE', 'LINEAR']
-            wcs.wcs.cunit = ['um', 'arcsec']
+            wcs.wcs.ctype = ["WAVE", "LINEAR"]
+            wcs.wcs.cunit = ["um", "arcsec"]
             wcs.wcs.crpix = [1, 1]
             wcs.wcs.cdelt = [bin_width, pixscale] # PIXSCALE
 
@@ -372,13 +376,14 @@ def rectify(self, hdulist, interps=None, wcs=None, **kwargs):
         rect_spec = np.zeros_like(Xarr, dtype=np.float32)
 
         ihdu = 0
+        # TODO: make this more iteratory
         for hdu in hdulist:
             if not isinstance(hdu, fits.ImageHDU):
                 continue
 
             wcs_fp = WCS(hdu.header, key="D")
-            n_x = hdu.header['NAXIS1']
-            n_y = hdu.header['NAXIS2']
+            n_x = hdu.header["NAXIS1"]
+            n_y = hdu.header["NAXIS2"]
             iarr, jarr = wcs_fp.all_world2pix(Xarr, Yarr, 0)
             mask = (iarr > 0) * (iarr < n_x) * (jarr > 0) * (jarr < n_y)
             if np.any(mask):
@@ -388,7 +393,7 @@ def rectify(self, hdulist, interps=None, wcs=None, **kwargs):
             ihdu += 1
 
         header = wcs.to_header()
-        header['EXTNAME'] = self.trace_id
+        header["EXTNAME"] = self.trace_id
         return fits.ImageHDU(data=rect_spec, header=header)
 
 
@@ -481,48 +486,117 @@ def footprint(self, wave_min=None, wave_max=None, xi_min=None, xi_max=None):
         return ([x_edge.min(), x_edge.max(), x_edge.max(), x_edge.min()],
                 [y_edge.min(), y_edge.min(), y_edge.max(), y_edge.max()])
 
-    def plot(self, wave_min=None, wave_max=None, c="r"):
-        """Plot control points of the SpectralTrace"""
+    def plot(self, wave_min=None, wave_max=None, xi_min=None, xi_max=None, *,
+             c="r", axes=None, plot_footprint=True, plot_wave=True,
+             plot_ctrlpnts=True, plot_outline=False, plot_trace_id=False):
+        """Plot control points (and/or footprint) of the SpectralTrace.
+
+        Parameters
+        ----------
+        wave_min : float, optional
+            Minimum wavelength, if any.
+        wave_max : float, optional
+            Maximum wavelength, if any.
+        xi_min : float, optional
+            Minimum slit, if any.
+        xi_max : float, optional
+            Maximum slit, if any.
+        c : str, optional
+            Colour, any valid matplotlib colour string. The default is "r".
+        axes : matplotlib axes, optional
+            The axes object to use for the plot. If None (default), a new
+            figure with one axes will be created.
+
+        Returns
+        -------
+        axes : matplotlib axes
+            The axes object containing the plot.
+
+        Other Parameters
+        ----------------
+        plot_footprint : bool, optional
+            Plot a rectangle encompassing all control points, which may be
+            larger than the area actually covered by the trace, if the trace is
+            not exactly perpendicular to the detector. The default is True.
+        plot_wave : bool, optional
+            Annotate the wavelength points. The default is True.
+        plot_ctrlpnts : bool, optional
+            Plot the individual control points as makers. The default is True.
+        plot_outline : bool, optional
+            Plot the smallest tetragon encompassing all control points.
+            The default is False.
+        plot_trace_id : bool, optional
+            Write the trace ID in the middle of the trace.
+            The default is False.
+        """
+        if axes is None:
+            _, axes = figure_factory()
 
         # Footprint (rectangle enclosing the trace)
         xlim, ylim  = self.footprint(wave_min=wave_min, wave_max=wave_max)
         if xlim is None:
-            return
+            return axes
+
+        if plot_footprint:
+            axes.plot(list(close_loop(xlim)), list(close_loop(ylim)))
 
-        xlim.append(xlim[0])
-        ylim.append(ylim[0])
-        plt.plot(xlim, ylim)
+        # for convenience...
+        xname = self.meta["x_colname"]
+        yname = self.meta["y_colname"]
+        wname = self.meta["wave_colname"]
+        sname = self.meta["s_colname"]
 
         # Control points
-        waves = self.table[self.meta["wave_colname"]]
+        waves = self.table[wname]
         if wave_min is None:
             wave_min = waves.min()
         if wave_max is None:
             wave_max = waves.max()
-
-        mask = (waves >= wave_min) * (waves <= wave_max)
-        if sum(mask) > 2:
-            w = waves[mask]
-
-            x = self.table[self.meta["x_colname"]][mask]
-            y = self.table[self.meta["y_colname"]][mask]
-            plt.plot(x, y, "o", c=c)
-
-            for wave in np.unique(w):
-                xx = x[w==wave]
-                xx.sort()
-                dx = xx[-1] - xx[-2]
-
-                plt.text(x[w==wave].max() + 0.5 * dx,
-                         y[w==wave].mean(),
-                         str(wave), va='center', ha='left')
-
-            plt.gca().set_aspect("equal")
-
-    @property
-    def trace_id(self):
-        """Return the name of the trace"""
-        return self.meta['trace_id']
+        xis = self.table[sname]
+        if xi_min is None:
+            xi_min = xis.min()
+        if xi_max is None:
+            xi_max = xis.max()
+
+        mask = ((waves >= wave_min) & (waves <= wave_max)
+                & (xis >= xi_min) & (xis <= xi_max))
+        if sum(mask) <= 2:
+            return axes
+
+        w = waves[mask]
+
+        x = self.table[xname][mask]
+        y = self.table[yname][mask]
+        if plot_ctrlpnts:
+            axes.plot(x, y, "o", c=c)
+
+        if plot_outline:
+            blue_end = self.table[mask][w == w.min()]
+            red_end = self.table[mask][w == w.max()]
+            blue_end.sort(sname)
+            red_end.sort(sname)
+            corners = vstack([blue_end[[0, -1]][xname, yname],
+                              red_end[[-1, 0]][xname, yname],
+                              blue_end[0][xname, yname]])
+            axes.plot(corners[xname], corners[yname], c=c)
+
+        if plot_trace_id:
+            axes.text(corners[xname][:-1].mean(), corners[yname][:-1].mean(),
+                      self.trace_id, c=c, rotation="vertical",
+                      ha="center", va="center")
+
+        for wave in np.unique(w):
+            xx = x[w==wave]
+            xx.sort()
+            dx = xx[-1] - xx[-2]
+
+            if plot_wave:
+                axes.text(x[w==wave].max() + 0.5 * dx,
+                          y[w==wave].mean(),
+                          str(wave), va="center", ha="left")
+
+        axes.set_aspect("equal")
+        return axes
 
     def _set_dispersion(self, wave_min, wave_max, pixsize=None):
         """Computation of dispersion dlam_per_pix along xi=0
@@ -538,8 +612,8 @@ def _set_dispersion(self, wave_min, wave_max, pixsize=None):
             dlam_grad = self.xy2lam.gradient()[0]  # dlam_by_dx
         else:
             dlam_grad = self.xy2lam.gradient()[1]  # dlam_by_dy
-        pixsize = (from_currsys(self.meta['pixel_scale']) /
-                   from_currsys(self.meta['plate_scale']))
+        pixsize = (from_currsys(self.meta["pixel_scale"]) /
+                   from_currsys(self.meta["plate_scale"]))
         self.dlam_per_pix = interp1d(lam,
                                      dlam_grad(x_mm, y_mm) * pixsize,
                                      fill_value="extrapolate")
@@ -548,7 +622,7 @@ def __repr__(self):
         return f"{self.__class__.__name__}({self.table!r}, **{self.meta!r})"
 
     def __str__(self):
-        msg = (f"<SpectralTrace> \"{self.meta['trace_id']}\" : "
+        msg = (f"<SpectralTrace> \"{self.trace_id}\" : "
                f"[{self.wave_min:.4f}, {self.wave_max:.4f}]um : "
                f"Ext {self.meta['extension_id']} : "
                f"Aperture {self.meta['aperture_id']} : "
@@ -892,8 +966,8 @@ def make_image_interpolations(hdulist, **kwargs):
     for hdu in hdulist:
         if isinstance(hdu, fits.ImageHDU):
             interps.append(
-                RectBivariateSpline(np.arange(hdu.header['NAXIS1']),
-                                    np.arange(hdu.header['NAXIS2']),
+                RectBivariateSpline(np.arange(hdu.header["NAXIS1"]),
+                                    np.arange(hdu.header["NAXIS2"]),
                                     hdu.data, **kwargs)
             )
     return interps
diff --git a/scopesim/effects/surface_list.py b/scopesim/effects/surface_list.py
index 7c5b2c44..252c5203 100644
--- a/scopesim/effects/surface_list.py
+++ b/scopesim/effects/surface_list.py
@@ -1,16 +1,17 @@
+from collections.abc import Collection
+
 import numpy as np
 from astropy import units as u
-import matplotlib.pyplot as plt
 
 from .ter_curves import TERCurve
 from ..optics import radiometry_utils as rad_utils
 from ..optics.surface import PoorMansSurface, SpectralSurface
-from ..utils import quantify, from_currsys
+from ..utils import quantify, from_currsys, figure_factory
 
 
 class SurfaceList(TERCurve):
     def __init__(self, **kwargs):
-        super(SurfaceList, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [20, 120, 520],
                   "minimum_throughput": "!SIM.spectral.minimum_throughput",
                   "etendue": "!TEL.etendue",
@@ -135,41 +136,48 @@ def add_surface_list(self, surface_list, prepend=False):
             self.table = rad_utils.combine_tables(surface_list.table,
                                                   self.table, prepend)
 
-    def plot(self, which="x", wavelength=None, ax=None, **kwargs):
-        """
+    def plot(self, which="x", wavelength=None, *, axes=None, **kwargs):
+        """Plot TER curves.
 
         Parameters
         ----------
-        which : str
-            "x" plots throughput. "t","e","r" plot trans/emission/refl
-        wavelength
-        kwargs
+        which : {"x", "t", "e", "r"}, optional
+            "x" plots throughput. "t","e","r" plot trans/emission/refl.
+            Can be a combination, e.g. "tr" or "tex" to plot each.
+        wavelength : array_like, optional
+            Passed to TERCurve.plot() for each surface. The default is None.
+        axes : matplotlib axes, optional
+            If given, plot into existing axes. The default is None.
 
         Returns
         -------
+        fig : matplotlib figure
+            Figure containing plots.
 
         """
-        import matplotlib.pyplot as plt
-        plt.gcf().clf()
-
-        for ii, which_part in enumerate(which):
-            ax = plt.subplot(len(which), 1, ii+1)
+        if axes is None:
+            fig, axes = figure_factory(len(which), 1, iterable_axes=True)
+        else:
+            fig = axes.figure
+            self._axes_guard(which, axes)
 
+        for ter, ax in zip(which, axes):
             # Plot the individual surfaces
-            for key in self.surfaces:
-                ter = TERCurve(**self.surfaces[key].meta)
-                ter.surface = self.surfaces[key]
-                ter.plot(which=which_part, wavelength=None, ax=None,
-                         new_figure=False,
-                         plot_kwargs={"ls": "-", "label": key}, **kwargs)
+            # TODO: do we want separate plots for these? if yes, how (row/col)?
+            for key, surface in self.surfaces.items():
+                curve = TERCurve(**surface.meta)
+                curve.surface = surface
+                kwargs.update(plot_kwargs={"ls": "-", "label": key})
+                curve.plot(ter, wavelength, axes=ax, **kwargs)
 
             # Plot the system surface
-            ter = TERCurve(**self.meta)
-            ter.surface = self.surface
-            ter.plot(which=which_part, wavelength=None, ax=ax, new_figure=False,
-                     plot_kwargs={"ls": "-.", "label": "System Throughput"},
-                     **kwargs)
+            # TODO: self is a subclass of TERCurve, why create again??
+            curve = TERCurve(**self.meta)
+            curve.surface = self.surface
+            kwargs.update(plot_kwargs={"ls": "-.",
+                                       "label": "System Throughput"})
+            curve.plot(ter, axes=ax, **kwargs)
 
-        plt.legend()
+        fig.legend()
 
-        return plt.gcf()
+        return fig
diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py
index 5d0e3335..cfd9a548 100644
--- a/scopesim/effects/ter_curves.py
+++ b/scopesim/effects/ter_curves.py
@@ -1,6 +1,7 @@
 """Transmission, emissivity, reflection curves"""
 import logging
 from pathlib import Path
+from collections.abc import Collection
 
 import numpy as np
 import skycalc_ipy
@@ -15,7 +16,8 @@
 from ..base_classes import SourceBase, FOVSetupBase
 from ..optics.surface import SpectralSurface
 from ..source.source import Source
-from ..utils import from_currsys, quantify, check_keys, find_file
+from ..utils import from_currsys, quantify, check_keys, find_file, \
+    figure_factory
 
 
 class TERCurve(Effect):
@@ -66,7 +68,7 @@ class TERCurve
 
     """
     def __init__(self, **kwargs):
-        super(TERCurve, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         params = {"z_order": [10, 110, 510],
                   "ignore_wings": False,
                   "wave_min": "!SIM.spectral.wave_min",
@@ -92,8 +94,6 @@ def __init__(self, **kwargs):
             self.surface.table = data
             self.surface.table.meta.update(self.meta)
 
-    # ####### added in new branch
-
     def apply_to(self, obj, **kwargs):
         if isinstance(obj, SourceBase):
             assert isinstance(obj, Source), "Only Source supported."
@@ -166,74 +166,63 @@ def background_source(self):
 
         return self._background_source
 
-    # #######
-
-    def plot(self, which="x", wavelength=None, ax=None, new_figure=True,
-             label=None, **kwargs):
-        """
+    def plot(self, which="x", wavelength=None, *, axes=None, **kwargs):
+        """Plot TER curves.
 
         Parameters
         ----------
-        which : str
-            "x" plots throughput. "t","e","r" plot trans/emission/refl
-        wavelength : list, np.ndarray
-        ax : matplotlib.Axis
-        new_figure : start a new figure (or add to the existing one)
-        label : the label to use (ignored)
-        kwargs
+        which : {"x", "t", "e", "r"}, optional
+            "x" plots throughput. "t","e","r" plot trans/emission/refl.
+            Can be a combination, e.g. "tr" or "tex" to plot each.
+        wavelength : array_like, optional
+            Wavelength on x-axis, taken from currsys if None (default).
+        axes : matplotlib axes, optional
+            If given, plot into existing axes. The default is None.
 
         Returns
         -------
+        fig : matplotlib figure
+            Figure containing plots.
 
         """
-        import matplotlib.pyplot as plt
-        if new_figure:
-            # from matplotlib import rcParams
-            # rcParams["axes.formatter.useoffset"] = False
-            plt.figure(figsize=(10, 5))
+        if axes is None:
+            fig, axes = figure_factory(len(which), 1, iterable_axes=True)
+        else:
+            fig = axes.figure
+            _guard_plot_axes(which, axes)
 
         self.meta.update(kwargs)
         params = from_currsys(self.meta)
 
-        for ii, ter in enumerate(which):
-            if ax is None:
-                plt.subplot(len(which), 1, ii+1)
-
-            if wavelength is None:
-                wunit = params["wave_unit"]
-                wave = np.arange(quantify(params["wave_min"], wunit).value,
-                                 quantify(params["wave_max"], wunit).value,
-                                 quantify(params["wave_bin"], wunit).value)
-                wave *= u.Unit(wunit)
-            else:
-                wave = wavelength
-
-            plot_kwargs = self.meta.get("plot_kwargs", {})
-            surf = self.surface
-
-            if "t" in ter:
-                y = surf.transmission(wave)
-            elif "e" in ter:
-                y = surf.emission(wave)
-            elif "r" in ter:
-                y = surf.reflection(wave)
-            else:
-                y = surf.throughput(wave)
+        wave_unit = self.meta.get("wavelength_unit")
+        if wavelength is None:
+            wunit = params["wave_unit"]
+            # TODO: shouldn't need both, make sure they're equal
+            assert wunit == wave_unit
+            wave = np.arange(quantify(params["wave_min"], wunit).value,
+                             quantify(params["wave_max"], wunit).value,
+                             quantify(params["wave_bin"], wunit).value)
+            wave *= u.Unit(wunit)
+        else:
+            wave = wavelength
 
-            plt.plot(wave, y, **plot_kwargs)
+        plot_kwargs = self.meta.get("plot_kwargs", {})
+        abbrs = {"t": "transmission", "e": "emission",
+                 "r": "reflection", "x": "throughput"}
 
-            wave_unit = self.meta.get("wavelength_unit")
-            plt.xlabel(f"Wavelength [{wave_unit}]")
-            y_str = {"t": "Transmission", "e": "Emission",
-                     "r": "Reflectivity", "x": "Throughput"}
-            plt.ylabel(f"{y_str[ter]} [{y.unit}]")
+        for ter, ax in zip(which, axes):
+            y_name = abbrs.get(ter, "throughput")
+            y = getattr(self.surface, y_name)
+            ax.plot(wave, y, **plot_kwargs)
+            ax.set_xlabel(f"Wavelength [{wave_unit}]")
+            ax.set_ylabel(f"{y_name.title()} [{y.unit}]")
 
-        return plt.gcf()
+        return fig
 
 
 class AtmosphericTERCurve(TERCurve):
     def __init__(self, **kwargs):
-        super(AtmosphericTERCurve, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = [111, 511]
         self.meta["action"] = "transmission"
         self.meta["position"] = 0       # position in surface table
@@ -271,7 +260,7 @@ class : SkycalcTERCurve
                 outer_unit : "m"
 
         """
-        super(SkycalcTERCurve, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["z_order"] = [112, 512]
         self.meta["use_local_skycalc_file"] = False
         self.meta.update(kwargs)
@@ -352,7 +341,7 @@ def query_server(self, **kwargs):
 
 class QuantumEfficiencyCurve(TERCurve):
     def __init__(self, **kwargs):
-        super(QuantumEfficiencyCurve, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         self.meta["action"] = "transmission"
         self.meta["z_order"] = [113, 513]
         self.meta["position"] = -1          # position in surface table
@@ -391,7 +380,7 @@ def __init__(self, **kwargs):
                                  " `array_dict`, `table`) or both "
                                  f"(`filter_name`, `filename_format`): {kwargs}")
 
-        super(FilterCurve, self).__init__(**kwargs)
+        super().__init__(**kwargs)
         if self.table is None:
             raise ValueError("Could not initialise filter. Either filename not "
                              "found, or array are not compatible")
@@ -508,9 +497,9 @@ def __init__(self, **kwargs):
 
         tbl = Table(names=["wavelength", "transmission"],
                     data=[waveset, transmission])
-        super(TopHatFilterCurve, self).__init__(table=tbl,
-                                                wavelength_unit="um",
-                                                action="transmission")
+        super().__init__(table=tbl,
+                         wavelength_unit="um",
+                         action="transmission")
         self.meta.update(kwargs)
 
 
@@ -520,7 +509,7 @@ def __init__(self, **kwargs):
         check_keys(kwargs, required_keys, action="error")
         filt_str = kwargs["filename_format"].format(kwargs["filter_name"])
         tbl = download_svo_filter(filt_str, return_style="table")
-        super(DownloadableFilterCurve, self).__init__(table=tbl, **kwargs)
+        super().__init__(table=tbl, **kwargs)
 
 
 class SpanishVOFilterCurve(FilterCurve):
@@ -555,48 +544,27 @@ def __init__(self, **kwargs):
         kwargs["svo_id"] = filt_str
 
         tbl = download_svo_filter(filt_str, return_style="table")
-        super(SpanishVOFilterCurve, self).__init__(table=tbl, **kwargs)
+        super().__init__(table=tbl, **kwargs)
 
 
-class FilterWheel(Effect):
-    """
-    This wheel holds a selection of predefined filters.
+class FilterWheelBase(Effect):
+    """Base class for Filter Wheels"""
 
-    Examples
-    --------
-    ::
-
-        name: filter_wheel
-        class: FilterWheel
-        kwargs:
-            filter_names: []
-            filename_format: "filters/{}.
-            current_filter: "Ks"
-
-    """
+    required_keys = set()
+    _current_str = "current_filter"
 
     def __init__(self, **kwargs):
-        required_keys = ["filter_names", "filename_format", "current_filter"]
-        check_keys(kwargs, required_keys, action="error")
-
-        super(FilterWheel, self).__init__(**kwargs)
+        super().__init__(**kwargs)
+        check_keys(kwargs, self.required_keys, action="error")
 
         params = {"z_order": [124, 224, 524],
-                  "path": "",
                   "report_plot_include": True,
                   "report_table_include": True,
                   "report_table_rounding": 4}
         self.meta.update(params)
         self.meta.update(kwargs)
 
-        path = Path(self.meta["path"], from_currsys(self.meta["filename_format"]))
         self.filters = {}
-        for name in from_currsys(self.meta["filter_names"]):
-            kwargs["name"] = name
-            self.filters[name] = FilterCurve(filename=str(path).format(name),
-                                             **kwargs)
-
-        self.table = self.get_table()
 
     def apply_to(self, obj, **kwargs):
         """Use apply_to of current filter"""
@@ -635,42 +603,42 @@ def current_filter(self):
             filter_eff = self.filters[filt_name]
         return filter_eff
 
-    @property
-    def display_name(self):
-        return (f"{self.meta['name']} : "
-                f"[{from_currsys(self.meta['current_filter'])}]")
-
     def __getattr__(self, item):
         return getattr(self.current_filter, item)
 
-    def plot(self, which="x", wavelength=None, **kwargs):
-        """
+    def plot(self, which="x", wavelength=None, *, axes=None, **kwargs):
+        """Plot TER curves.
 
         Parameters
         ----------
-        which : str
-            "x" plots throughput. "t","e","r" plot trans/emission/refl
-        wavelength
-        kwargs
+        which : {"x", "t", "e", "r"}, optional
+            "x" plots throughput. "t","e","r" plot trans/emission/refl.
+            Can be a combination, e.g. "tr" or "tex" to plot each.
+        wavelength : array_like, optional
+            DESCRIPTION. The default is None.
+        axes : matplotlib axes, optional
+            If given, plot into existing axes. The default is None.
 
         Returns
         -------
+        fig : matplotlib figure
+            Figure containing plots.
 
         """
-        import matplotlib.pyplot as plt
-        plt.figure(figsize=(10, 5))
+        if axes is None:
+            fig, axes = figure_factory(len(which), 1, iterable_axes=True)
+        else:
+            fig = axes.figure
+            _guard_plot_axes(which, axes)
 
-        for ii, ter in enumerate(which):
-            ax = plt.subplot(len(which), 1, ii+1)
-            for name, _filter in self.filters.items():
-                _filter.plot(which=ter, wavelength=wavelength, ax=ax,
-                             new_figure=False, plot_kwargs={"label": name},
-                             **kwargs)
 
-        # plt.semilogy()
-        plt.legend()
+        for ter, ax in zip(which, axes):
+            for name, _filter in self.filters.items():
+                _filter.plot(which=ter, wavelength=wavelength, axes=ax,
+                             plot_kwargs={"label": name}, **kwargs)
 
-        return plt.gcf()
+        fig.legend()
+        return fig
 
     def get_table(self):
         names = list(self.filters.keys())
@@ -686,7 +654,42 @@ def get_table(self):
         return tbl
 
 
-class TopHatFilterWheel(FilterWheel):
+class FilterWheel(FilterWheelBase):
+    """
+    This wheel holds a selection of predefined filters.
+
+    Examples
+    --------
+    ::
+
+        name: filter_wheel
+        class: FilterWheel
+        kwargs:
+            filter_names: []
+            filename_format: "filters/{}.
+            current_filter: "Ks"
+
+    """
+
+    required_keys = {"filter_names", "filename_format", "current_filter"}
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        params = {"path": ""}
+        self.meta.update(params)
+        self.meta.update(kwargs)
+
+        path = self._get_path()
+        for name in from_currsys(self.meta["filter_names"]):
+            kwargs["name"] = name
+            self.filters[name] = FilterCurve(filename=str(path).format(name),
+                                             **kwargs)
+
+        self.table = self.get_table()
+
+
+class TopHatFilterWheel(FilterWheelBase):
     """
     A selection of top-hat filter curves as defined in the input lists
 
@@ -708,7 +711,7 @@ class TopHatFilterWheel(FilterWheel):
 
     current_filter: str, optional
         Name of current filter at initialisation. If no name is given, the
-        first entry in ``filter_names`` is used by default.
+        first entry in `filter_names` is used by default.
 
     Examples
     --------
@@ -725,35 +728,30 @@ class TopHatFilterWheel(FilterWheel):
             current_filter: "K"
 
     """
-    def __init__(self, **kwargs):
-        required_keys = ["filter_names", "transmissions", "wing_transmissions",
-                         "blue_cutoffs", "red_cutoffs"]
-        check_keys(kwargs, required_keys, action="error")
 
-        super(FilterWheel, self).__init__(**kwargs)
+    required_keys = {"filter_names", "transmissions", "wing_transmissions",
+                     "blue_cutoffs", "red_cutoffs"}
 
-        params = {"z_order": [124, 224, 524],
-                  "report_plot_include": True,
-                  "report_table_include": True,
-                  "report_table_rounding": 4,
-                  "current_filter": kwargs.get("current_filter",
-                                               kwargs["filter_names"][0])
-                  }
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        current_filter = kwargs.get("current_filter",
+                                    kwargs["filter_names"][0])
+        params = {"current_filter": current_filter}
         self.meta.update(params)
         self.meta.update(kwargs)
 
-        self.filters = {}
-        for i in range(len(self.meta["filter_names"])):
-            name = self.meta["filter_names"][i]
-            effect_kwargs = {"name": name,
-                             "transmission": self.meta["transmissions"][i],
-                             "wing_transmission": self.meta["wing_transmissions"][i],
-                             "blue_cutoff": self.meta["blue_cutoffs"][i],
-                             "red_cutoff": self.meta["red_cutoffs"][i]}
+        for i_filt, name in enumerate(self.meta["filter_names"]):
+            effect_kwargs = {
+                "name": name,
+                "transmission": self.meta["transmissions"][i_filt],
+                "wing_transmission": self.meta["wing_transmissions"][i_filt],
+                "blue_cutoff": self.meta["blue_cutoffs"][i_filt],
+                "red_cutoff": self.meta["red_cutoffs"][i_filt]}
             self.filters[name] = TopHatFilterCurve(**effect_kwargs)
 
 
-class SpanishVOFilterWheel(FilterWheel):
+class SpanishVOFilterWheel(FilterWheelBase):
     """
     A FilterWheel that loads all the filters from the Spanish VO service
 
@@ -795,37 +793,32 @@ class SpanishVOFilterWheel(FilterWheel):
             include_str: "_filter"
 
     """
-    def __init__(self, **kwargs):
-        required_keys = ["observatory", "instrument", "current_filter"]
-        check_keys(kwargs, required_keys, action="error")
 
-        # Call Effect.init, *NOT* FilterWheel.init --> different required_keys
-        super(FilterWheel, self).__init__(**kwargs)
+    required_keys = {"observatory", "instrument", "current_filter"}
 
-        params = {"z_order": [124, 224, 524],
-                  "report_plot_include": True,
-                  "report_table_include": True,
-                  "report_table_rounding": 4,
-                  "include_str": None,         # passed to
+    def __init__(self, **kwargs):        
+        super().__init__(**kwargs)
+
+        params = {"include_str": None,         # passed to
                   "exclude_str": None,
                   }
         self.meta.update(params)
         self.meta.update(kwargs)
 
         obs, inst = self.meta["observatory"], self.meta["instrument"]
-        inc, exc = self.meta["include_str"], self.meta["exclude_str"]
-        filter_names = download_svo_filter_list(obs, inst, short_names=True,
-                                                include=inc, exclude=exc)
+        filter_names = download_svo_filter_list(
+            obs, inst, short_names=True,
+            include=self.meta["include_str"], exclude=self.meta["exclude_str"])
 
         self.meta["filter_names"] = filter_names
-        self.filters = {name: SpanishVOFilterCurve(observatory=obs,
-                                                   instrument=inst,
-                                                   filter_name=name)
-                        for name in filter_names}
-        self.filters["open"] = FilterCurve(array_dict={"wavelength": [0.3, 3.0],
-                                                       "transmission": [1., 1.]},
-                                           wavelength_unit="um",
-                                           name="unity transmission")
+        for name in filter_names:
+            self.filters[name] = SpanishVOFilterCurve(observatory=obs,
+                                                      instrument=inst,
+                                                      filter_name=name)
+
+        self.filters["open"] = FilterCurve(
+            array_dict={"wavelength": [0.3, 3.0], "transmission": [1., 1.]},
+            wavelength_unit="um", name="unity transmission")
 
         self.table = self.get_table()
 
@@ -870,11 +863,13 @@ class ADCWheel(Effect):
            filename_format: "TER_ADC_{}.dat"
            current_adc: "const_90"
     """
-    def __init__(self, **kwargs):
-        required_keys = ["adc_names", "filename_format", "current_adc"]
-        check_keys(kwargs, required_keys, action="error")
 
+    required_keys = {"adc_names", "filename_format", "current_adc"}
+    _current_str = "current_adc"
+
+    def __init__(self, **kwargs):
         super().__init__(**kwargs)
+        check_keys(kwargs, self.required_keys, action="error")
 
         params = {"z_order": [125, 225, 525],
                   "path": "",
@@ -884,7 +879,7 @@ def __init__(self, **kwargs):
         self.meta.update(params)
         self.meta.update(kwargs)
 
-        path = Path(self.meta["path"], from_currsys(self.meta["filename_format"]))
+        path = self._get_path()
         self.adcs = {}
         for name in from_currsys(self.meta["adc_names"]):
             kwargs["name"] = name
@@ -913,11 +908,6 @@ def current_adc(self):
             return False
         return self.adcs[curradc]
 
-    @property
-    def display_name(self):
-        return (f"{self.meta['name']} : "
-                f"[{from_currsys(self.meta['current_adc'])}]")
-
     def __getattr__(self, item):
         return getattr(self.current_adc, item)
 
@@ -930,3 +920,16 @@ def get_table(self):
         tbl = Table(names=["name", "max_transmission"],
                     data=[names, tmax])
         return tbl
+
+
+def _guard_plot_axes(which, axes):
+    if len(which) > 1:
+        if not isinstance(axes, Collection):
+            raise TypeError(("axes must be collection of axes if which "
+                             "contains more than one element"))
+        if not len(axes) == len(which):
+            raise ValueError("len of which and axes must match")
+    else:
+        if isinstance(axes, Collection):
+            raise TypeError(("axes must be a single axes object if which "
+                             "contains only one element"))
diff --git a/scopesim/optics/fov.py b/scopesim/optics/fov.py
index 73af7cd3..6a87d70f 100644
--- a/scopesim/optics/fov.py
+++ b/scopesim/optics/fov.py
@@ -1,5 +1,7 @@
 """Defines FieldOfView class"""
+import logging
 from copy import deepcopy
+from itertools import chain
 
 import numpy as np
 from scipy.interpolate import interp1d
@@ -83,14 +85,19 @@ def __init__(self, header, waverange, detector_header=None, **kwargs):
         self._wavelength = None
         self._volume = None
 
+    def _pixarea(self, hdr):
+        return (hdr["CDELT1"] * u.Unit(hdr["CUNIT1"]) *
+                hdr["CDELT2"] * u.Unit(hdr["CUNIT2"])).to(u.arcsec ** 2)
+
+    @property
+    def trace_id(self):
+        """Return the name of the trace"""
+        return self.meta["trace_id"]
+
     @property
     def pixel_area(self):
         if self.meta["pixel_area"] is None:
-            hdr = self.header
-            pixarea = (hdr["CDELT1"] * u.Unit(hdr["CUNIT1"]) *
-                       hdr["CDELT2"] * u.Unit(hdr["CUNIT2"])).to(u.arcsec ** 2)
-            self.meta["pixel_area"] = pixarea.value     # [arcsec]
-
+            self.meta["pixel_area"] = self._pixarea(self.header).value  # [arcsec]
         return self.meta["pixel_area"]
 
     def extract_from(self, src):
@@ -106,13 +113,16 @@ def extract_from(self, src):
 
         fields_in_fov = [field for field in src.fields
                          if fu.is_field_in_fov(self.header, field)]
+        if not fields_in_fov:
+            logging.warning("No fields in FOV.")
 
-        spec_refs = []
+        spec_refs = set()
         volume = self.volume()
         for ifld, fld in enumerate(fields_in_fov):
             if isinstance(fld, Table):
-                fields_in_fov[ifld] = fu.extract_area_from_table(fld, volume)
-                spec_refs += list(np.unique(fields_in_fov[ifld] ["ref"]))
+                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):
@@ -120,11 +130,11 @@ def extract_from(self, src):
                 if fld.header["NAXIS"] == 2 or fld.header.get("BG_SRC"):
                     ref = fld.header.get("SPEC_REF")
                     if ref is not None:
-                        spec_refs += [ref]
+                        spec_refs.add(ref)
 
         waves = volume["waves"] * u.Unit(volume["wave_unit"])
         spectra = {ref: fu.extract_range_from_spectrum(src.spectra[ref], waves)
-                   for ref in np.unique(spec_refs)}
+                   for ref in spec_refs}
 
         self.fields = fields_in_fov
         self.spectra = spectra
@@ -148,6 +158,7 @@ def view(self, hdu_type="image", sub_pixel=None, use_photlam=None):
             self.meta["sub_pixel"] = sub_pixel
 
         if hdu_type == "image":
+            # FIXME: why not just make False the default value??
             use_photlam = False if use_photlam is None else use_photlam
             self.hdu = self.make_image_hdu(use_photlam=use_photlam)
         elif hdu_type == "cube":
@@ -158,82 +169,163 @@ def view(self, hdu_type="image", sub_pixel=None, use_photlam=None):
         return self.hdu
 
     def flatten(self):
+        """If cube, collapse along first axis."""
         if self.hdu and self.hdu.header["NAXIS"] == 3:
             image = np.sum(self.hdu.data, axis=0)
             self.hdu.data = image
 
-    def make_spectrum(self):
-        """
-        This is needed for when we do incoherent MOS instruments.
-        Each fibre doesn't care about the spatial information.
-
-        1. Make waveset and zeros flux
-            make dict of spectra evaluated at waveset
-
-        2. Find Cube fields
-            collapse cube along spatial dimensions --> spectrum vector
-            convert vector to PHOTLAM
-            interpolate at waveset
-            add to zeros flux vector
+    def _evaluate_spectrum_with_weight(self, ref, waveset, weight):
+        return self.spectra[ref](waveset).value * weight
 
-        3. Find Image fields
-            sum image over both dimensions
-            evaluate SPEC_REF spectum at waveset
-            multiply by sum
-            add to zeros flux vector
-
-        4. Find Table fields
-            evaluate all spectra at waveset
-            for each unique ref, sum the weights
-            add each spectra * sum of weights to zeros flux vector
-
-        Returns
-        -------
-        spec : SourceSpectrum
-            [PHOTLAM]
+    def _calc_area_factor(self, field):
+        bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2)
+        area_factor = self.pixel_area * bg_solid_angle  # arcsec**2 * arcsec**-2
+        return area_factor
 
+    def _make_spectrum_cubefields(self, fov_waveset):
         """
-        fov_waveset = self.waveset
-        canvas_flux = np.zeros(len(fov_waveset))
-
-        specs = {ref: spec(fov_waveset) for ref, spec in self.spectra.items()}
+        Find Cube fields.
 
+        * collapse cube along spatial dimensions --> spectrum vector
+        * convert vector to PHOTLAM
+        * interpolate at waveset
+        * 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)
+            hdu_waveset = fu.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)
 
             field_unit = field.header.get("BUNIT", PHOTLAM)
             flux_scale_factor = u.Unit(field_unit).to(PHOTLAM)
 
-            canvas_flux += fov_waveset_fluxes * flux_scale_factor
+            yield fov_waveset_fluxes * flux_scale_factor
 
+    def _make_spectrum_imagefields(self, waveset):
+        """
+        Find Image fields.
+
+        * sum image over both dimensions
+        * evaluate SPEC_REF spectum at waveset
+        * yield spectrum multiply by sum to be added to canvas flux
+        """
         for field in self.image_fields:
-            ref = field.header["SPEC_REF"]
             weight = np.sum(field.data)
-            canvas_flux += self.spectra[ref](fov_waveset).value * weight
+            yield self._evaluate_spectrum_with_weight(field.header["SPEC_REF"],
+                                                      waveset, weight)
 
+    def _make_spectrum_tablefields(self, waveset):
+        """
+        Find Table fields.
 
+        * evaluate all spectra at waveset
+        * for each unique ref, sum the weights
+        * yield each spectrum * sum of weights to be added to canvas flux
+        """
         for field in self.table_fields:
             refs = np.array(field["ref"])
             weights = np.array(field["weight"])
-            weight_sums = {ref: np.sum(weights[refs == ref])
-                           for ref in np.unique(refs)}
-            for ref, weight in weight_sums.items():
-                canvas_flux += self.spectra[ref](fov_waveset).value * weight
+            # TODO: could do grouping of table with both columns??
+            for ref in set(refs):
+                weight = np.sum(weights, where=refs==ref)
+                yield self._evaluate_spectrum_with_weight(ref, waveset, weight)
 
+    def _make_spectrum_backfields(self, waveset):
         for field in self.background_fields:
-            bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2)
-            area_factor = self.pixel_area * bg_solid_angle       # arcsec**2 * arcsec**-2
+            yield self._evaluate_spectrum_with_weight(
+                field.header["SPEC_REF"], waveset,
+                self._calc_area_factor(field))
+
+    def make_spectrum(self):
+        """
+        This is needed for when we do incoherent MOS instruments.
+        Each fibre doesn't care about the spatial information.
 
-            ref = field.header["SPEC_REF"]
-            canvas_flux += self.spectra[ref](fov_waveset).value * area_factor
+        Returns
+        -------
+        spec : SourceSpectrum
+            [PHOTLAM]
+
+        """
+        fov_waveset = self.waveset
+        # Start with zero flux no ensure correct array shape even if none of
+        # the sub-functions yield anything.
+        canvas_flux = sum(chain(
+            self._make_spectrum_cubefields(fov_waveset),
+            self._make_spectrum_imagefields(fov_waveset),
+            self._make_spectrum_tablefields(fov_waveset),
+            self._make_spectrum_backfields(fov_waveset),
+            ), start=np.zeros_like(fov_waveset.value))
 
         spectrum = SourceSpectrum(Empirical1D, points=fov_waveset,
                                   lookup_table=canvas_flux)
-
         return spectrum
 
+    def _make_image_cubefields(self, area):
+        """
+        Find Cube fields.
+
+        * collapse cube along wavelength axis
+        * rescale and reorient image
+        * yield cube image  to be added to canvas image
+        """
+        for field in self.cube_fields:
+            # cube_fields come in with units of photlam/arcsec2, need to convert to ph/s
+            # We need to the voxel volume (spectral and solid angle) for that.
+            # ..todo: implement branch for use_photlam is True
+            spectral_bin_width = (field.header["CDELT3"] *
+                                  u.Unit(field.header["CUNIT3"])).to(u.Angstrom)
+            # First collapse to image, then convert units
+            image = np.sum(field.data, axis=0) * PHOTLAM/u.arcsec**2
+            image = (image * self._pixarea(field.header) * area *
+                     spectral_bin_width).to(u.ph/u.s)
+            yield fits.ImageHDU(data=image, header=field.header)
+
+    def _make_image_imagefields(self, fluxes):
+        """
+        Find Image fields.
+
+        * sum spectra between wavelength edges
+        * multiply image by summed spectrum
+        * yield image  to be added to canvas image
+        """
+        for field in self.image_fields:
+            image = deepcopy(field.data)
+            image *= fluxes[field.header["SPEC_REF"]]  # ph / s
+            yield fits.ImageHDU(data=image, header=field.header)
+
+    def _make_image_tablefields(self, fluxes):
+        """
+        Find Table fields.
+
+        * sum spectra between wavelength edges
+        * yield summed flux at x,y position to be added to canvas image
+        """
+        for field in self.table_fields:
+            # x, y are ALWAYS in arcsec - crval is in deg
+            xpix, ypix = imp_utils.val2pix(self.header,
+                                           field["x"] / 3600,
+                                           field["y"] / 3600)
+            if utils.from_currsys(self.meta["sub_pixel"]):
+                for idx, row in enumerate(field):
+                    xs, ys, fracs = imp_utils.sub_pixel_fractions(xpix[idx],
+                                                                  ypix[idx])
+                    for x, y, frac in zip(xs, ys, fracs):
+                        yield fluxes[row["ref"]] * frac, row["weight"], x, y
+            else:
+                # Note: these had x/ypix+0.5 until a06ab75
+                # TODO: could these be something more numpythonic grid-ish?
+                x = np.array(xpix).astype(int)
+                y = np.array(ypix).astype(int)     # quickest way to round
+                flux = np.array([fluxes[ref] for ref in field["ref"]])
+                yield flux, np.array(field["weight"]), x, y
+
+    def _make_image_backfields(self, fluxes):
+        for field in self.background_fields:
+            flux = fluxes[field.header["SPEC_REF"]]
+            yield flux * self._calc_area_factor(field)
+
     def make_image_hdu(self, use_photlam=False):
         """
         Used for imaging
@@ -244,23 +336,7 @@ def make_image_hdu(self, use_photlam=False):
 
             See make_cube for an explanation
 
-        1. Make waveset and canvas image
-            make canvas image from NAXIS1,2 from fov.header
-
-        2. Find Cube fields
-            collapse cube along wavelength axis
-            rescale and reorient image
-            add cube image to canvas image
-
-        3. Find Image fields
-            rescale and reorient images
-            sum spectra between wavelength edges
-            multiply image by summed spectrum
-            add image to canvas image
-
-        4. Find Table fields
-            sum spectra between wavelength edges
-            add summed flux at x,y position in canvas image
+        Make canvas image from NAXIS1,2 from fov.header
 
         Parameters
         ----------
@@ -275,91 +351,131 @@ def make_image_hdu(self, use_photlam=False):
         """
         spline_order = utils.from_currsys("!SIM.computing.spline_order")
 
-        # 1. Make waveset and canvas image
-        fov_waveset = np.unique(self.waveset)
+        # 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 = utils.from_currsys(self.meta["area"])    # u.m2
 
         # PHOTLAM * u.um * u.m2 --> ph / s
-        specs = {ref: spec(fov_waveset) for ref, spec in self.spectra.items()}
-        if use_photlam is False:
-            for key in specs:
-                specs[key] = (specs[key] * bin_widths * area).to(u.ph / u.s)
-
-        fluxes = {ref: np.sum(spec).value for ref, spec in specs.items()}
-        naxis1, naxis2 = self.header["NAXIS1"], self.header["NAXIS2"]
-        canvas_image_hdu = fits.ImageHDU(data=np.zeros((naxis2, naxis1)),
-                                         header=self.header)
-
-        # 2. Find Cube fields
-        for field in self.cube_fields:
-            # cube_fields come in with units of photlam/arcsec2, need to convert to ph/s
-            # We need to the voxel volume (spectral and solid angle) for that.
-            # ..todo: implement branch for use_photlam is True
-            spectral_bin_width = (field.header["CDELT3"] *
-                                  u.Unit(field.header["CUNIT3"])).to(u.Angstrom)
-            pixarea = (field.header["CDELT1"] * u.Unit(field.header["CUNIT1"]) *
-                       field.header["CDELT2"] * u.Unit(field.header["CUNIT2"])).to(u.arcsec**2)
-
-            # First collapse to image, then convert units
-            image = np.sum(field.data, axis=0) * PHOTLAM/u.arcsec**2
-            image = (image * pixarea * area * spectral_bin_width).to(u.ph/u.s)
+        specs = {ref: spec(fov_waveset) if use_photlam
+                 else (spec(fov_waveset) * bin_widths * 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()}
+        canvas_image_hdu = fits.ImageHDU(
+            data=np.zeros((self.header["NAXIS2"], self.header["NAXIS1"])),
+            header=self.header)
 
-            tmp_hdu = fits.ImageHDU(data=image, header=field.header)
+        for tmp_hdu in chain(self._make_image_cubefields(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)
 
-        # 2. Find Image fields
+        for flux, weight, x, y in self._make_image_tablefields(fluxes):
+            if utils.from_currsys(self.meta["sub_pixel"]):
+                canvas_image_hdu.data[y, x] += flux * weight
+            else:
+                # Mask out any stars that were pushed out of the fov by rounding
+                mask = ((x < canvas_image_hdu.data.shape[1]) *
+                        (y < canvas_image_hdu.data.shape[0]))
+                canvas_image_hdu.data[y[mask], x[mask]] += flux[mask] * weight[mask]
+
+        canvas_image_hdu.data = sum(self._make_image_backfields(fluxes),
+                                    start=canvas_image_hdu.data)
+        return canvas_image_hdu  # [ph s-1 pixel-1]
+
+    def _make_cube_cubefields(self, fov_waveset):
+        """
+        Find Cube fields.
+
+        * rescale and reorient cubes
+        * interp1d smaller cubes with waveset
+        * yield cubes to be added to cavas cube
+        """
+        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)
+
+            # ..todo: Deal with this bounds_error in a more elegant way
+            field_interp = interp1d(field_waveset.to(u.um).value,
+                                    field.data, axis=0, kind="linear",
+                                    bounds_error=False, fill_value=0)
+
+            field_data = field_interp(fov_waveset.value)
+
+            # Pixel scale conversion
+            field_data *= self._pixarea(field.header).value / self.pixel_area
+            field_hdu = fits.ImageHDU(data=field_data, header=field.header)
+            yield field_hdu
+
+    def _make_cube_imagefields(self, specs, spline_order):
+        """
+        Find Image fields.
+
+        * rescale and reorient images
+        * evaluate spectra at waveset
+        * expand image by spectra to 3D form
+        * yield image cubes to be added to cavas cube
+        """
         for field in self.image_fields:
-            image = deepcopy(field.data)
-            tmp_hdu = fits.ImageHDU(data=image, header=field.header)
-            tmp_hdu.data *= fluxes[field.header["SPEC_REF"]]  # ph / s
+            # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping
+            # Assumption is that ImageHDUs have units of PHOTLAM arcsec-2
+            # ImageHDUs have photons/second/pixel.
+            # ..todo: Add a catch to get ImageHDU with BUNITs
+            canvas_image_hdu = fits.ImageHDU(
+                data=np.zeros((self.header["NAXIS2"], self.header["NAXIS1"])),
+                header=self.header)
+            # FIXME: Why here not "Pixel scale conversion" as above?
+            field.data /= self.pixel_area
             canvas_image_hdu = imp_utils.add_imagehdu_to_imagehdu(
-                tmp_hdu,
+                field,
                 canvas_image_hdu,
-                conserve_flux=True,
                 spline_order=spline_order)
+            spec = specs[field.header["SPEC_REF"]]
+            field_cube = canvas_image_hdu.data[None, :, :] * spec[:, None, None]  # 2D * 1D -> 3D
+            yield field_cube.value
 
-        # 3. Find Table fields
+    def _make_cube_tablefields(self, specs):
+        """
+        Find Table fields.
+
+        * evaluate spectra at waveset
+        * yield spectrum at x,y position to be added cavas to cube
+        """
         for field in self.table_fields:
-            # x, y are ALWAYS in arcsec - crval is in deg
-            xpix, ypix = imp_utils.val2pix(self.header,
-                                           field["x"] / 3600,
-                                           field["y"] / 3600)
-            if utils.from_currsys(self.meta["sub_pixel"]):
-                for i, row in enumerate(field):
-                    xs, ys, fracs = imp_utils.sub_pixel_fractions(xpix[i],
-                                                                  ypix[i])
-                    ref, weight = row["ref"], row["weight"]
-                    for x, y, f in zip(xs, ys, fracs):
-                        canvas_image_hdu.data[y, x] += fluxes[ref] * weight * f
-            else:
-                # Note: these had x/ypix+0.5 until a06ab75
-                x = np.array(xpix).astype(int)
-                y = np.array(ypix).astype(int)     # quickest way to round
-                f = np.array([fluxes[ref] for ref in field["ref"]])
-                weight = np.array(field["weight"])
+            # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping
+            # Point sources are in PHOTLAM per pixel
+            # Point sources need to be scaled up by inverse pixel_area
+            for row in field:
+                xsky, ysky = row["x"] / 3600, row["y"] / 3600
+                # x, y are ALWAYS in arcsec - crval is in deg
+                xpix, ypix = imp_utils.val2pix(self.header, xsky, ysky)
+                flux_vector = specs[row["ref"]].value * row["weight"] / self.pixel_area
 
-                # Mask out any stars that were pushed out of the fov by rounding
-                m = (x < canvas_image_hdu.data.shape[1]) * \
-                    (y < canvas_image_hdu.data.shape[0])
-                canvas_image_hdu.data[y[m], x[m]] += f[m] * weight[m]
+                if utils.from_currsys(self.meta["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
+                else:
+                    yield flux_vector, int(xpix), int(ypix)
 
-        # 4. Find Background fields
+    def _make_cube_backfields(self, specs):
         for field in self.background_fields:
-            bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2)
-            area_factor = self.pixel_area * bg_solid_angle       # arcsec**2 * arcsec**-2
-
-            flux = fluxes[field.header["SPEC_REF"]] * area_factor
-            canvas_image_hdu.data += flux
-
-        image_hdu = canvas_image_hdu        # [ph s-1 pixel-1]
+            # TODO: The following would have been identical to the other two
+            #       make methods, but was commented out. Why?
+            # bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2)  # float [arcsec-2]
+            # pixel_area = utils.from_currsys(self.meta["pixel_scale"]) ** 2      # float [arcsec2]
+            # area_factor = pixel_area * bg_solid_angle                           # float [arcsec2 * arcsec-2]
 
-        return image_hdu
+            # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping
+            # spec = specs[field.header["SPEC_REF"]] * area_factor
+            spec = specs[field.header["SPEC_REF"]]
+            yield spec[:, None, None].value
 
     def make_cube_hdu(self):
         """
@@ -386,23 +502,9 @@ def make_cube_hdu(self):
                 make waveset from self.meta values
             make canvas cube based on waveset of largest cube and NAXIS1,2 from fov.header
 
-        2. Find Cube fields::
-
-            rescale and reorient cubes
-            interp1d smaller cubes with waveset
-            add cubes to cavas cube
-
-        3. Find Image fields::
-
-            rescale and reorient images
-            evaluate spectra at waveset
-            expand image by spectra to 3D form
-            add image cubes to canvas cube
-
-        4. Find Table fields::
-
-            evaluate spectra at waveset
-            add spectrum at x,y position in canvas cube
+        2. Find Cube fields (see ``FieldOfView._make_cube_cubefields()``).
+        3. Find Image fields (see ``FieldOfView._make_cube_imagefields()``).
+        4. Find Table fields (see ``FieldOfView._make_cube_tablefields()``).
 
         ``PHOTLAM = ph/s/m2/um``.
         Original source fields are in units of:
@@ -423,15 +525,14 @@ def make_cube_hdu(self):
 
 
         # 1. Make waveset and canvas cube (area, bin_width are applied at end)
-
-        wave_min = self.meta["wave_min"]        # Quantity [um]
-        wave_max = self.meta["wave_max"]
-
+        # TODO: Why is this not self.waveset? What's different?
         wave_unit = u.Unit(utils.from_currsys("!SIM.spectral.wave_unit"))
-        dwave = utils.from_currsys("!SIM.spectral.spectral_bin_width")  # Not a quantity
-        fov_waveset = np.arange(wave_min.value, wave_max.value, dwave) * wave_unit
+        fov_waveset = np.arange(
+            self.meta["wave_min"].value, self.meta["wave_max"].value,
+            utils.from_currsys("!SIM.spectral.spectral_bin_width")) * wave_unit
         fov_waveset = fov_waveset.to(u.um)
 
+        # TODO: what's with this code??
         # fov_waveset = self.waveset
         # wave_bin_n = len(fov_waveset)
         # if "lin" in self.meta["wave_bin_type"]:
@@ -444,90 +545,30 @@ def make_cube_hdu(self):
                  for ref, spec in self.spectra.items()}
 
         # make canvas cube based on waveset of largest cube and NAXIS1,2 from fov.header
-        naxis1, naxis2 = self.header["NAXIS1"], self.header["NAXIS2"]
-        naxis3 = len(fov_waveset)
-        canvas_cube_hdu = fits.ImageHDU(data=np.zeros((naxis3, naxis2, naxis1)),
-                                        header=self.header)
+        canvas_cube_hdu = fits.ImageHDU(
+            data=np.zeros((len(fov_waveset),
+                           self.header["NAXIS2"],
+                           self.header["NAXIS1"])),
+            header=self.header)
         canvas_cube_hdu.header["BUNIT"] = "ph s-1 cm-2 AA-1"
 
-        # 2. Add Cube fields
-        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)
-
-            # ..todo: Deal with this bounds_error in a more elegant way
-            field_interp = interp1d(field_waveset.to(u.um).value,
-                                    field.data, axis=0, kind="linear",
-                                    bounds_error=False, fill_value=0)
-
-            field_data = field_interp(fov_waveset.value)
-
-            # Pixel scale conversion
-            field_pixarea = (field.header["CDELT1"]
-                             * field.header["CDELT2"]
-                             * u.Unit(field.header["CUNIT1"])
-                             * u.Unit(field.header["CUNIT2"])).to(u.arcsec**2)
-            field_pixarea = field_pixarea.value
-            field_data *= field_pixarea / self.pixel_area
-            field_hdu = fits.ImageHDU(data=field_data, header=field.header)
+        for field_hdu in self._make_cube_cubefields(fov_waveset):
             canvas_cube_hdu = imp_utils.add_imagehdu_to_imagehdu(
                 field_hdu,
                 canvas_cube_hdu,
                 spline_order=spline_order)
 
-        # 3. Find Image fields
-        for field in self.image_fields:
-            # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping
-            # Assumption is that ImageHDUs have units of PHOTLAM arcsec-2
-            # ImageHDUs have photons/second/pixel.
-            # ..todo: Add a catch to get ImageHDU with BUNITs
-            canvas_image_hdu = fits.ImageHDU(data=np.zeros((naxis2, naxis1)),
-                                             header=self.header)
-            pixarea = (field.header["CDELT1"] * u.Unit(field.header["CUNIT1"]) *
-                       field.header["CDELT2"] * u.Unit(field.header["CUNIT2"])).to(u.arcsec**2)
-
-            field.data = field.data / self.pixel_area
-            canvas_image_hdu = imp_utils.add_imagehdu_to_imagehdu(field,
-                                                    canvas_image_hdu,
-                                                    spline_order=spline_order)
-            spec = specs[field.header["SPEC_REF"]]
-            field_cube = canvas_image_hdu.data[None, :, :] * spec[:, None, None]  # 2D * 1D -> 3D
-            canvas_cube_hdu.data += field_cube.value
+        canvas_cube_hdu.data = sum(self._make_cube_imagefields(specs,
+                                                               spline_order),
+                                   start=canvas_cube_hdu.data)
 
-        # 4. Find Table fields
-        for field in self.table_fields:
-            # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping
-            # Point sources are in PHOTLAM per pixel
-            # Point sources need to be scaled up by inverse pixel_area
-            pixel_area = self.pixel_area
-            for row in field:
-                xsky, ysky = row["x"], row["y"]
-                ref, weight = row["ref"], row["weight"]
-                # x, y are ALWAYS in arcsec - crval is in deg
-                xpix, ypix = imp_utils.val2pix(self.header, xsky / 3600, ysky / 3600)
-                if utils.from_currsys(self.meta["sub_pixel"]):
-                    xs, ys, fracs = imp_utils.sub_pixel_fractions(xpix, ypix)
-                    for i, j, k in zip(xs, ys, fracs):
-                        flux_vector = specs[ref].value * weight * k / pixel_area
-                        canvas_cube_hdu.data[:, j, i] +=  flux_vector
-                else:
-                    x, y = int(xpix), int(ypix)
-                    flux_vector = specs[ref].value * weight / pixel_area
-                    canvas_cube_hdu.data[:, y, x] += flux_vector
+        for flux, x, y in self._make_cube_tablefields(specs):
+            canvas_cube_hdu.data[:, y, x] += flux
 
-        # 5. Add Background fields
-        for field in self.background_fields:
-            # bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2)  # float [arcsec-2]
-            # pixel_area = utils.from_currsys(self.meta["pixel_scale"]) ** 2      # float [arcsec2]
-            # area_factor = pixel_area * bg_solid_angle                           # float [arcsec2 * arcsec-2]
-
-            # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping
-            # spec = specs[field.header["SPEC_REF"]] * area_factor
-            spec = specs[field.header["SPEC_REF"]]
-            canvas_cube_hdu.data += spec[:, None, None].value
+        canvas_cube_hdu.data = sum(self._make_cube_backfields(specs),
+                                   start=canvas_cube_hdu.data)
 
+        # TODO: what's with this code??
         # 6. Convert from PHOTLAM to ph/s/voxel
         #    PHOTLAM = ph/s/cm-2/AA
         #    area = m2, fov_waveset = um
@@ -536,6 +577,7 @@ def make_cube_hdu(self):
         canvas_cube_hdu.data *= 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??
         # bin_widths = np.diff(fov_waveset).to(u.AA).value
         # bin_widths = 0.5 * (np.r_[0, bin_widths] + np.r_[bin_widths, 0])
         # canvas_cube_hdu.data *= bin_widths[:, None, None]
@@ -547,14 +589,12 @@ def make_cube_hdu(self):
                                        "CUNIT3": "um",
                                        "CTYPE3": "WAVE"})
         # ..todo: Add the log wavelength keyword here, if log scale is needed
-
-        cube_hdu = canvas_cube_hdu      # [ph s-1 AA-1 (arcsec-2)]
-
-        return cube_hdu
+        return canvas_cube_hdu      # [ph s-1 AA-1 (arcsec-2)]
 
     def volume(self, wcs_prefix=""):
         xs, ys = imp_utils.calc_footprint(self.header, wcs_suffix=wcs_prefix)
-        wave_corners = self.waverange
+        # FIXME: This is unused!!
+        # wave_corners = self.waverange
         self._volume = {"xs": [min(xs), max(xs)],
                         "ys": [min(ys), max(ys)],
                         "waves": self.waverange,
@@ -564,28 +604,27 @@ def volume(self, wcs_prefix=""):
 
     @property
     def data(self):
+        """Return either hdu.data, image, cube, spectrum or None."""
         if self.hdu is not None:
-            data = self.hdu.data
-        elif self.image is not None:
-            data = self.image
-        elif self.cube is not None:
-            data = self.cube
-        elif self.spectrum is not None:
-            data = self.spectrum
-        else:
-            data = None
-
-        return data
+            return self.hdu.data
+        if self.image is not None:
+            return self.image
+        if self.cube is not None:
+            return self.cube
+        if self.spectrum is not None:
+            return self.spectrum
+        return None
 
     @property
     def corners(self):
+        """Return sky footprint, image plane footprint."""
         sky_corners = imp_utils.calc_footprint(self.header)
         imp_corners = imp_utils.calc_footprint(self.header, "D")
         return sky_corners, imp_corners
 
     @property
     def waverange(self):
-        """Returns wavelength range in um [wave_min, wave_max]"""
+        """Return wavelength range in um [wave_min, wave_max]."""
         if self._waverange is None:
             wave_min = utils.quantify(self.meta["wave_min"], u.um).value
             wave_max = utils.quantify(self.meta["wave_max"], u.um).value
@@ -594,55 +633,63 @@ def waverange(self):
 
     @property
     def wavelength(self):
-        """Returns central wavelength in um"""
+        """Return central wavelength in um."""
         if self._wavelength is None:
             self._wavelength = np.average(self.waverange)
         return utils.quantify(self._wavelength, u.um)
 
     @property
     def waveset(self):
-        """Returns a wavelength vector in um"""
-
-        field_cubes = self.cube_fields
-        if len(field_cubes) > 0:
-            i = np.argmax([cube.header["NAXIS3"] for cube in field_cubes])
-            _waveset = fu.get_cube_waveset(field_cubes[i].header,
+        """Return a wavelength vector in um."""
+        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)
-        elif len(self.spectra) > 0:
+        elif self.spectra:
             wavesets = [spec.waveset for spec in self.spectra.values()]
             _waveset = np.concatenate(wavesets)
         else:
             _waveset = self.waverange * u.um
 
-        _waveset = np.unique(_waveset)
-
-        return _waveset.to(u.um)
+        # TODO: tie the round to a global precision setting (this is defined
+        #       better in another TODO somewhere...)
+        # The rounding is necessary to not end up with things like:
+        #   0.7 um
+        #   0.7000000000000001 um
+        #   0.7000000000000002 um
+        # and yes, that actually happend...
+        _waveset = np.unique(_waveset.to(u.um).round(10))
+        return _waveset
 
     @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 field.header.get("BG_SRC", False) is False]
+                and not field.header.get("BG_SRC", False)]
 
     @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 field.header.get("BG_SRC", False) is False]
+                and not field.header.get("BG_SRC", False)]
 
     @property
     def table_fields(self):
+        """Return list of Table fields."""
         return [field for field in self.fields
                 if isinstance(field, Table)]
 
-
     @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) is True]
+                and field.header.get("BG_SRC", False)]
 
     def __repr__(self):
         waverange = [self.meta["wave_min"].value, self.meta["wave_max"].value]
diff --git a/scopesim/optics/fov_manager.py b/scopesim/optics/fov_manager.py
index 819385e9..6c361af9 100644
--- a/scopesim/optics/fov_manager.py
+++ b/scopesim/optics/fov_manager.py
@@ -43,10 +43,11 @@
 # """
 
 from copy import deepcopy
-import numpy as np
 from typing import TextIO
 from io import StringIO
+from collections.abc import Iterable, MutableSequence
 
+import numpy as np
 from astropy import units as u
 
 from . import image_plane_utils as ipu
@@ -72,7 +73,7 @@ class FOVManager:
     All observation parameters as passed from UserCommands
 
     """
-    def __init__(self, effects=[], **kwargs):
+    def __init__(self, effects=None, **kwargs):
         self.meta = {"area": "!TEL.area",
                      "pixel_scale": "!INST.pixel_scale",
                      "plate_scale": "!INST.plate_scale",
@@ -94,9 +95,9 @@ def __init__(self, effects=[], **kwargs):
         params["meta"] = from_currsys({key: self.meta[key] for key in fvl_meta})
         self.volumes_list = FovVolumeList(initial_volume=params)
 
-        self.effects = effects
+        self.effects = effects or []
         self._fovs_list = []
-        self.is_spectroscope = eu.is_spectroscope(effects)
+        self.is_spectroscope = eu.is_spectroscope(self.effects)
 
         if from_currsys(self.meta["preload_fovs"]) is True:
             self._fovs_list = self.generate_fovs_list()
@@ -173,11 +174,11 @@ def fovs(self):
         return self._fovs_list
 
     @property
-    def fov_footprints(self, which="both"):
+    def fov_footprints(self):
         return None
 
 
-class FovVolumeList(FOVSetupBase):
+class FovVolumeList(FOVSetupBase, MutableSequence):
     """
     List of FOV volumes for FOVManager
 
@@ -211,9 +212,9 @@ def __init__(self, initial_volume=None):
                                 "yd_min": 0,
                                 "yd_max": 0}
 
-    def split(self, axis, value, aperture_id=None):
+    def split(self, axis, value, aperture_id=None) -> None:
         """
-        Splits the all volumes that include axis=value into two.
+        Split the all volumes that include axis=value into two.
 
         - Loop through all volume dict
         - Find any entries where min < value < max
@@ -221,8 +222,8 @@ def split(self, axis, value, aperture_id=None):
 
         Parameters
         ----------
-        axis : str, list of str
-            "wave", "x", "y"
+        axis : {"wave", "x", "y"}, or list thereof
+            Which axis (``str``) or axes (``list[str]``) to use.
         value : float, list of floats
         aperture_id : int, optional
             Default None. If ``None``, split all volumes. If ``int``, only split
@@ -244,20 +245,25 @@ def split(self, axis, value, aperture_id=None):
         if isinstance(axis, (tuple, list)):
             for ax, val in zip(axis, value):
                 self.split(ax, val)
-        elif isinstance(value, (tuple, list, np.ndarray)):
+            return
+
+        if isinstance(value, Iterable):
             for val in value:
                 self.split(axis, val)
-        else:
-            for i, vol_old in enumerate(self.volumes):
-                if aperture_id in (vol_old["meta"]["aperture_id"], None)  \
-                        and vol_old[f"{axis}_min"] < value \
-                        and vol_old[f"{axis}_max"] > value:
-                    vol_new = deepcopy(vol_old)
-                    vol_new[f"{axis}_min"] = value
-                    vol_old[f"{axis}_max"] = value
-                    self.volumes.insert(i+1, vol_new)
-
-    def shrink(self, axis, values, aperture_id=None):
+            return
+
+        for vol in self:
+            if (aperture_id is not None and
+                aperture_id != vol["meta"]["aperture_id"]):
+                continue
+            if vol[f"{axis}_min"] >= value or vol[f"{axis}_max"] <= value:
+                continue
+            new_vol = deepcopy(vol)
+            new_vol[f"{axis}_min"] = value
+            vol[f"{axis}_max"] = value
+            self.insert(self.index(vol) + 1, new_vol)
+
+    def shrink(self, axis, values, aperture_id=None) -> None:
         """
         - Loop through all volume dict
         - Replace any entries where min < values.min
@@ -265,8 +271,8 @@ def shrink(self, axis, values, aperture_id=None):
 
         Parameters
         ----------
-        axis : str
-            "wave", "x", "y"
+        axis : {"wave", "x", "y"} or list thereof
+            Which axis (``str``) or axes (``list[str]``) to use.
         values : list of 2 floats
             [min, max], [min, None], [None, max]
         aperture_id : int, optional
@@ -283,47 +289,52 @@ def shrink(self, axis, values, aperture_id=None):
 
 
         """
+        # FIXME: Isn't this method just the same as setting self.volumes to the
+        #        output list of self.extract()?? Except the None values.
         if isinstance(axis, (tuple, list)):
             for ax, val in zip(axis, values):
                 self.shrink(ax, val)
-        else:
-            to_pop = []
+            return
+
+        to_pop = []
+
+        for vol in self:
+            if (aperture_id is not None and
+                aperture_id != vol["meta"]["aperture_id"]):
+                continue
 
             if values[0] is not None:
-                for i, vol in enumerate(self.volumes):
-                    if aperture_id in (vol["meta"]["aperture_id"], None):
-                        if vol[f"{axis}_max"] <= values[0]:
-                            to_pop.append(i)
-                        elif vol[f"{axis}_min"] < values[0]:
-                            vol[f"{axis}_min"] = values[0]
+                if vol[f"{axis}_max"] <= values[0]:
+                    to_pop.append(self.index(vol))
+                    continue
+                vol[f"{axis}_min"] = max(values[0], vol[f"{axis}_min"])
 
             if values[1] is not None:
-                for i, vol in enumerate(self.volumes):
-                    if aperture_id in (vol["meta"]["aperture_id"], None):
-                        if vol[f"{axis}_min"] >= values[1]:
-                            to_pop.append(i)
-                        if vol[f"{axis}_max"] > values[1]:
-                            vol[f"{axis}_max"] = values[1]
+                if vol[f"{axis}_min"] >= values[1]:
+                    to_pop.append(self.index(vol))
+                    continue
+                vol[f"{axis}_max"] = min(values[1], vol[f"{axis}_max"])
 
-            for i in sorted(to_pop)[::-1]:
-                self.volumes.pop(i)
+        for idx in reversed(sorted(to_pop)):
+            self.pop(idx)
 
     def extract(self, axes, edges, aperture_id=None):
         """
-        Returns new volumes from within all existing volumes
+        Return new volumes from within all existing volumes.
 
         This method DOES NOT alter the existing self.volumes list
         To include the returned volumes, add them to the self.volumes list
 
         Parameters
         ----------
-        axes : str, list of str
-            "wave", "x", "y"
+        axes : list of either {"wave", "x", "y"}
+            Which axis (``list`` of single ``str``) or axes (``list[str]``)
+            to use. Must be ``list`` in either case.
         edges : list, tuple of lists
             Edge points for each axes listed
         aperture_id : int, optional
             Default None. If ``None``, extract from all volumes. If ``int``,
-            only extract from volumes with this ``aperture_id`` in the meta dict
+            only extract from volumes with this `aperture_id` in the meta dict
 
         Examples
         --------
@@ -345,38 +356,36 @@ def extract(self, axes, edges, aperture_id=None):
             A list of all new volumes extracted from existing volumes
 
         """
-        new_vols = []
-        for old_vol in self.volumes:
-            if aperture_id in (old_vol["meta"]["aperture_id"], None):
-                add_flag = True
-                new_vol = deepcopy(old_vol)
+        def _get_new_vols():
+            for vol in self:
+                if (aperture_id is not None and
+                    aperture_id != vol["meta"]["aperture_id"]):
+                    continue
+                if not all(_volume_in_range(vol, axis, edge) for axis, edge
+                           in zip(axes, edges)):
+                    continue
+                new_vol = deepcopy(vol)
                 for axis, edge in zip(axes, edges):
-                    if edge[0] <= old_vol[f"{axis}_max"] and \
-                       edge[1] >= old_vol[f"{axis}_min"]:
-                        new_vol[f"{axis}_min"] = max(edge[0], old_vol[f"{axis}_min"])
-                        new_vol[f"{axis}_max"] = min(edge[1], old_vol[f"{axis}_max"])
-                    else:
-                        add_flag = False
+                    new_vol[f"{axis}_min"] = max(edge[0], vol[f"{axis}_min"])
+                    new_vol[f"{axis}_max"] = min(edge[1], vol[f"{axis}_max"])
+                yield new_vol
 
-                if add_flag is True:
-                    new_vols.append(new_vol)
-
-        return new_vols
+        return list(_get_new_vols())
 
     def __len__(self):
         return len(self.volumes)
 
-    def __iter__(self):
-        return iter(self.volumes)
+    def __getitem__(self, index):
+        return self.volumes[index]
 
-    def __getitem__(self, key):
-        return self.volumes[key]
+    def __setitem__(self, index, value):
+        self.volumes[index] = value
 
-    def __setitem__(self, key, value):
-        self.volumes[item] = value
+    def __delitem__(self, index):
+        del self.volumes[index]
 
-    def __delitem__(self, key):
-        del self.volumes[key]
+    def insert(self, index, value):
+        self.volumes.insert(index, value)
 
     def write_string(self, stream: TextIO) -> None:
         """Write formatted string representation to I/O stream"""
@@ -403,21 +412,18 @@ def __str__(self) -> str:
             output = str_stream.getvalue()
         return output
 
-    def __iadd__(self, other):
-        if isinstance(other, list):
-            self.volumes += other
-        else:
-            raise ValueError(f"Can only add lists of volumes: {other}")
-
-        return self
-
     def __add__(self, other):
+        # TODO: Is this functionality actually used anywhere?
         new_self = deepcopy(self)
         new_self += other
 
         return new_self
 
 
+def _volume_in_range(vol: dict, axis: str, edge) -> bool:
+    return edge[0] <= vol[f"{axis}_max"] and edge[1] >= vol[f"{axis}_min"]
+
+
 # Spectroscopy FOV setup
 # ======================
 #
diff --git a/scopesim/optics/fov_manager_utils.py b/scopesim/optics/fov_manager_utils.py
index 05165939..f97db5aa 100644
--- a/scopesim/optics/fov_manager_utils.py
+++ b/scopesim/optics/fov_manager_utils.py
@@ -1,4 +1,7 @@
+import logging
 from copy import deepcopy
+from itertools import product
+from more_itertools import pairwise
 
 import numpy as np
 from astropy import units as u
@@ -9,6 +12,8 @@
 from ..effects.effects_utils import get_all_effects
 from ..utils import check_keys
 
+# TODO: Where are all these functions used??
+
 
 def get_3d_shifts(effects, **kwargs):
     """
@@ -41,9 +46,11 @@ def get_3d_shifts(effects, **kwargs):
         shifts = [eff.fov_grid(which="shifts", **kwargs) for eff in effects]
 
         old_bin_edges = [shift[0] for shift in shifts if len(shift[0]) >= 2]
+        # TODO: could this use combine_wavesets?
         new_bin_edges = np.unique(np.sort(np.concatenate(old_bin_edges),
                                           kind="stable"))
 
+        # TODO: could this be zeros_like?
         x_shifts = np.zeros(len(new_bin_edges))
         y_shifts = np.zeros(len(new_bin_edges))
         # .. todo:: replace the 1e-7 with a variable in !SIM
@@ -99,25 +106,27 @@ def get_imaging_waveset(effects_list, **kwargs):
 
     wave_bin_edges = [filt.fov_grid(which="waveset", **kwargs)
                       for filt in filters]
-    if len(wave_bin_edges) > 0:
-        kwargs["wave_min"] = np.max([w[0].value for w in wave_bin_edges])
-        kwargs["wave_max"] = np.min([w[1].value for w in wave_bin_edges])
+    if wave_bin_edges:
+        kwargs["wave_min"] = max(wave[0].value for wave in wave_bin_edges)
+        kwargs["wave_max"] = min(wave[1].value for wave in wave_bin_edges)
+    # Bit confusing...
     wave_bin_edges = [[kwargs["wave_min"], kwargs["wave_max"]]]
 
     if kwargs["wave_min"] > kwargs["wave_max"]:
-        raise ValueError(f"Filter wavelength ranges do not overlap: {wave_bin_edges}")
+        raise ValueError("Filter wavelength ranges do not overlap: "
+                         f"{wave_bin_edges}.")
 
     # ..todo: add in Atmospheric dispersion and ADC here
     for effect_class in [efs.PSF]:
-        effects = get_all_effects(effects_list, effect_class)
-        for eff in effects:
+        for eff in get_all_effects(effects_list, effect_class):
             waveset = eff.fov_grid(which="waveset", **kwargs)
             if waveset is not None:
-                wave_bin_edges += [waveset]
+                wave_bin_edges.append(waveset)
 
-    wave_bin_edges = combine_wavesets(wave_bin_edges)
+    wave_bin_edges = combine_wavesets(*wave_bin_edges)
 
-    if len(wave_bin_edges) == 0:
+    if not wave_bin_edges:
+        # This is already set at the top, why again here?
         wave_bin_edges = [kwargs["wave_min"], kwargs["wave_max"]]
 
     return wave_bin_edges
@@ -125,7 +134,7 @@ def get_imaging_waveset(effects_list, **kwargs):
 
 def get_imaging_headers(effects, **kwargs):
     """
-    Returns a list of Header objects for each of the FieldOfVIew objects
+    Return a generator of Header objects for each of the FieldOfVIew objects.
 
     Parameters
     ----------
@@ -135,7 +144,7 @@ def get_imaging_headers(effects, **kwargs):
 
     Returns
     -------
-    hdrs : list of Header objects
+    hdrs : generator of Header objects
 
     Notes
     -----
@@ -163,33 +172,39 @@ def get_imaging_headers(effects, **kwargs):
     aperture_effects = get_all_effects(effects, (efs.ApertureMask,
                                                  efs.SlitWheel,
                                                  efs.ApertureList))
-    if len(aperture_effects) == 0:
+    if not aperture_effects:
         detector_arrays = get_all_effects(effects, efs.DetectorList)
-        if len(detector_arrays) > 0:
-            aperture_effects += [detarr.fov_grid(which="edges",
-                                                 pixel_scale=pixel_scale)
-                                 for detarr in detector_arrays]
-        else:
-            raise ValueError("No ApertureMask or DetectorList was provided. At "
-                             "least one must be passed to make an ImagePlane: "
-                             f"{effects}")
-
+        if not detector_arrays:
+            raise ValueError("No ApertureMask or DetectorList was provided. "
+                             "At least one must be passed to make an "
+                             f"ImagePlane: {effects}.")
+        aperture_effects.extend(
+            detarr.fov_grid(which="edges", pixel_scale=pixel_scale)
+            for detarr in detector_arrays)
+
+    # FIXME: all of this is a bit inconsistent; fov_grid(which="edges" is
+    #        called afterwards, but when looking in detector_arrays, the same
+    #        is called immediately; does that even work? is this all tested??
     # get aperture headers from fov_grid()
     # - for-loop catches mutliple headers from ApertureList.fov_grid()
-    hdrs = []
-    for ap_eff in aperture_effects:
-        # ..todo:: add this functionality to ApertureList effect
-        hdr = ap_eff.fov_grid(which="edges", pixel_scale=pixel_scale)
-        hdrs += hdr if isinstance(hdr, (list, tuple)) else [hdr]
+    def _get_hdrs(ap_effs):
+        for ap_eff in ap_effs:
+            # ..todo:: add this functionality to ApertureList effect
+            hdr = ap_eff.fov_grid(which="edges", pixel_scale=pixel_scale)
+            if isinstance(hdr, (list, tuple)):
+                yield from hdr
+            else:
+                yield hdr
+    headers = _get_hdrs(aperture_effects)
 
     # check size of aperture in pixels - split if necessary
-    sky_hdrs = []
-    for hdr in hdrs:
-        if hdr["NAXIS1"] * hdr["NAXIS2"] > kwargs["max_segment_size"]:
-            sky_hdrs += imp_utils.split_header(hdr, kwargs["chunk_size"])
-        else:
-            sky_hdrs += [hdr]
-
+    def _get_sky_hdrs(hdrs):
+        for hdr in hdrs:
+            if hdr["NAXIS1"] * hdr["NAXIS2"] > kwargs["max_segment_size"]:
+                yield from imp_utils.split_header(hdr, kwargs["chunk_size"])
+            else:
+                yield hdr
+    sky_hdrs = _get_sky_hdrs(headers)
     # ..todo:: Deal with the case that two or more ApertureMasks overlap
 
     # map the on-sky apertures directly to the image plane using plate_scale
@@ -208,13 +223,12 @@ def get_imaging_headers(effects, **kwargs):
 
         dethdr = imp_utils.header_from_list_of_xy(x_det, y_det, pixel_size, "D")
         skyhdr.update(dethdr)
-
-    return sky_hdrs
+        yield skyhdr
 
 
 def get_imaging_fovs(headers, waveset, shifts, **kwargs):
     """
-    Returns a list of ``FieldOfView`` objects
+    Return a generator of ``FieldOfView`` objects.
 
     Parameters
     ----------
@@ -224,54 +238,50 @@ def get_imaging_fovs(headers, waveset, shifts, **kwargs):
     waveset : list of floats
         [um] N+1 wavelengths for N spectral layers
 
-    shifts : list of tuples
+    shifts : list of tuples (or actually arrays?)
         [deg] x,y shifts w.r.t to the optical axis plane. N shifts for N
         spectral layers
 
     Returns
     -------
-    fovs : list of FieldOfView objects
+    fovs : generator of ``FieldOfView`` objects
 
     """
-
-    shift_waves = shifts["wavelengths"]     # in [um]
+    # Ensure array for later indexing
+    shift_waves = np.array(shifts["wavelengths"])     # in [um]
     shift_dx = shifts["x_shifts"]           # in [deg]
     shift_dy = shifts["y_shifts"]
 
     # combine the wavelength bins from 1D spectral effects and 3D shift effects
-    if len(shifts["wavelengths"]) > 0:
-        mask = (shift_waves > np.min(waveset)) * (shift_waves < np.max(waveset))
-        waveset = combine_wavesets([waveset, shift_waves[mask]])
-
-    counter = 0
-    fovs = []
+    if shift_waves.size:
+        mask = (shift_waves > min(waveset)) * (shift_waves < max(waveset))
+        waveset = combine_wavesets(waveset, shift_waves[mask])
 
-    print(f"Preparing {(len(waveset)-1)*len(headers)} FieldOfViews", flush=True)
+    # Actually evaluating the generators here is only necessary for the log msg
+    waveset = list(waveset)
+    headers = list(headers)
+    logging.info("Preparing %d FieldOfViews", (len(waveset) - 1) * len(headers))
 
-    for ii in range(len(waveset) - 1):
-        for hdr in headers:
-            # add any pre-instrument shifts to the FOV sky coords
-            wave_mid = 0.5 * (waveset[ii] + waveset[ii+1])
-            x_shift = np.interp(wave_mid, shift_waves, shift_dx)
-            y_shift = np.interp(wave_mid, shift_waves, shift_dy)
-
-            fov_hdr = deepcopy(hdr)
-            fov_hdr["CRVAL1"] += x_shift      # headers are in [deg]
-            fov_hdr["CRVAL2"] += y_shift
+    combos = product(pairwise(waveset), headers)
+    for fov_id, ((wave_min, wave_max), hdr) in enumerate(combos):
+        # add any pre-instrument shifts to the FOV sky coords
+        wave_mid = 0.5 * (wave_min + wave_max)
+        x_shift = np.interp(wave_mid, shift_waves, shift_dx)
+        y_shift = np.interp(wave_mid, shift_waves, shift_dy)
 
-            # define the wavelength range for the FOV
-            waverange = [waveset[ii], waveset[ii + 1]]
+        fov_hdr = deepcopy(hdr)
+        fov_hdr["CRVAL1"] += x_shift      # headers are in [deg]
+        fov_hdr["CRVAL2"] += y_shift
 
-            # Make the FOV
-            fov = FieldOfView(fov_hdr, waverange, id=counter, **kwargs)
-            fovs += [fov]
-            counter += 1
+        # define the wavelength range for the FOV
+        waverange = [wave_min, wave_max]
 
-    return fovs
+        # Make the FOV
+        yield FieldOfView(fov_hdr, waverange, id=fov_id, **kwargs)
 
 
 def get_spectroscopy_headers(effects, **kwargs):
-
+    """Return generator of Header objects."""
     required_keys = ["pixel_scale", "plate_scale",
                      "wave_min", "wave_max"]
     check_keys(kwargs, required_keys, action="error")
@@ -285,15 +295,15 @@ def get_spectroscopy_headers(effects, **kwargs):
                                                  efs.SlitWheel,
                                                  efs.ApertureMask))
 
-    if len(surface_list_effects) > 0:
+    if surface_list_effects:
         waves = surface_list_effects[0].fov_grid(which="waveset")
         if len(waves) == 2:
             kwargs["wave_min"] = np.max([waves[0].value, kwargs["wave_min"]])
             kwargs["wave_max"] = np.min([waves[1].value, kwargs["wave_max"]])
 
-    if len(detector_list_effects) > 0:
+    if detector_list_effects:
         implane_hdr = detector_list_effects[0].image_plane_header
-    elif len(spec_trace_effects) > 0:
+    elif spec_trace_effects:
         implane_hdr = spec_trace_effects[0].image_plane_header
     else:
         raise ValueError("Missing a way to determine the image plane size")
@@ -304,6 +314,7 @@ def get_spectroscopy_headers(effects, **kwargs):
                          f"{spec_trace_effects}")
     spec_trace = spec_trace_effects[0]
 
+    # TODO: The following is WET with the code in get_imaging_headers
     sky_hdrs = []
     for ap_eff in aperture_effects:
         # if ApertureList, a list of ApertureMask headers is returned
@@ -320,19 +331,23 @@ def get_spectroscopy_headers(effects, **kwargs):
                                        plate_scale=kwargs["plate_scale"]
                                        )
                    for sky_hdr in sky_hdrs]
-    fov_headers = [hdr for hdr_list in fov_headers for hdr in hdr_list]
-    # ..todo: check that each header is not larger than chunk_size
-
-    return fov_headers
+    for hdr_list in fov_headers:
+        for hdr in hdr_list:
+            yield hdr
+    # TODO: check that each header is not larger than chunk_size
+    # that's already done in get_imaging_headers, isn't it?
 
 
-def get_spectroscopy_fovs(headers, shifts, effects=[], **kwargs):
+def get_spectroscopy_fovs(headers, shifts, effects=None, **kwargs):
+    """Return a generator of ``FieldOfView`` objects."""
+    if effects is None:
+        effects = []
 
     shift_waves = shifts["wavelengths"]     # in [um]
     shift_dx = shifts["x_shifts"]           # in [deg]
     shift_dy = shifts["y_shifts"]
 
-    print(f"Preparing {len(headers)} FieldOfViews", flush=True)
+    logging.info("Preparing %d FieldOfViews", len(headers))
 
     apertures = get_all_effects(effects, (efs.ApertureList, efs.ApertureMask))
     masks = [ap.fov_grid(which="masks") for ap in apertures]
@@ -343,8 +358,7 @@ def get_spectroscopy_fovs(headers, shifts, effects=[], **kwargs):
         elif isinstance(mask, np.ndarray):
             mask_dict[len(mask_dict)] = mask
 
-    fovs = []
-    for ii, hdr in enumerate(headers):
+    for fov_id, hdr in enumerate(headers):
         # add any pre-instrument shifts to the FOV sky coords
         wave_mid = hdr["WAVE_MID"]
         x_shift = np.interp(wave_mid, shift_waves, shift_dx)
@@ -355,29 +369,30 @@ def get_spectroscopy_fovs(headers, shifts, effects=[], **kwargs):
         fov_hdr["CRVAL2"] += y_shift
 
         # Make the FOV
-        fov = FieldOfView(fov_hdr, waverange=[hdr["WAVE_MIN"], hdr["WAVE_MAX"]],
-                          **kwargs)
+        waverange = [hdr["WAVE_MIN"], hdr["WAVE_MAX"]]
+        fov = FieldOfView(fov_hdr, waverange=waverange, **kwargs)
         fov.meta["distortion"]["rotation"] = hdr["ROTANGD"]
         fov.meta["distortion"]["shear"] = hdr["SKEWANGD"]
         fov.meta["conserve_image"] = hdr["IMG_CONS"]
-        fov.meta["fov_id"] = ii
+        # TODO: In the other function, the id is set via the contructor.
+        #       What's the difference?
+        fov.meta["fov_id"] = fov_id
         fov.meta["aperture_id"] = hdr["APERTURE"]
 
         # .. todo: get these masks working
         # there needs to be fov_grid(which="mask") in ApertureList/Mask
         # fov.mask = mask_dict[hdr["APERTURE"]]
-        fovs += [fov]
+        yield fov
 
-    return fovs
 
-
-def combine_wavesets(waveset_list):
+# FIXME: This functions doesn't seem to be covered by any separate unit test.
+def combine_wavesets(*wavesets):
     """
-    Joins and sorts several sets of wavelengths into a single 1D array
+    Join and sorts several sets of wavelengths into a single 1D array.
 
     Parameters
     ----------
-    waveset_list : list
+    wavesets : one or more iterables
         A group of wavelength arrays or lists
 
     Returns
@@ -385,12 +400,22 @@ def combine_wavesets(waveset_list):
     wave_set : np.ndarray
         Combined set of wavelengths
 
+    Note
+    ----
+    This assumes that all wavesets are given in the same unit!
     """
-    wave_set = []
-    for wbe in waveset_list:        # wbe = waveset bin edges
-        wbe = wbe.value if isinstance(wbe, u.Quantity) else wbe
-        wave_set += list(wbe)
-    # ..todo:: set variable in !SIM.computing for rounding to the 7th decimal
-    wave_set = np.unique(np.round(np.sort(wave_set, kind="stable"), 7))
-
+    # TODO: set variable in !SIM.computing for rounding to the 7th decimal
+    decimals = 7
+
+    def _get_waves(waves):
+        for wave in waves:
+            if isinstance(wave, u.Quantity):
+                round_wave = wave.round(decimals).value
+            else:
+                round_wave = np.round(wave, decimals)
+            yield from round_wave
+
+    # NOTE: This function previously used np.sort(wave_set, kind="stable").
+    #       If any issues occur with the buitin sorted, go back to that!
+    wave_set = sorted(set(_get_waves(wavesets)))
     return wave_set
diff --git a/scopesim/optics/fov_utils.py b/scopesim/optics/fov_utils.py
index 8c554b49..86747eb4 100644
--- a/scopesim/optics/fov_utils.py
+++ b/scopesim/optics/fov_utils.py
@@ -329,7 +329,11 @@ def extract_area_from_imagehdu(imagehdu, fov_volume):
 
         # OC [2021-12-14] if fov range is not covered by the source return nothing
         if not np.any(mask):
-            print(f"FOV {fov_waves[0]} um - {fov_waves[1]} um: not covered by Source")
+            logging.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]
@@ -360,6 +364,9 @@ def extract_area_from_imagehdu(imagehdu, fov_volume):
         data = imagehdu.data[y0p:y1p, x0p:x1p]
         new_hdr["SPEC_REF"] = hdr.get("SPEC_REF")
 
+    if not data.size:
+        logging.warning("Empty image HDU.")
+
     new_imagehdu = fits.ImageHDU(data=data)
     new_imagehdu.header.update(new_hdr)
 
diff --git a/scopesim/optics/image_plane_utils.py b/scopesim/optics/image_plane_utils.py
index 116c6c18..6944175f 100644
--- a/scopesim/optics/image_plane_utils.py
+++ b/scopesim/optics/image_plane_utils.py
@@ -756,14 +756,15 @@ def val2pix(header, a, b, wcs_suffix=""):
     if isinstance(header, fits.ImageHDU):
         header = header.header
 
-    if "PC1_1"+s in header:
-        pc11 = header["PC1_1"+s]
-        pc12 = header["PC1_2"+s]
-        pc21 = header["PC2_1"+s]
-        pc22 = header["PC2_2"+s]
+    pckeys = [key + s for key in ["PC1_1", "PC1_2", "PC2_1", "PC2_2"]]
+    if all(key in header for key in pckeys):
+        pc11, pc12, pc21, pc22 = (header[key] for key in pckeys)
     else:
         pc11, pc12, pc21, pc22 = 1, 0, 0, 1
 
+    if (pc11 * pc22 + pc12 * pc21) != 1.0:
+        logging.error("PC matrix det != 1.0")
+
     da = float(header["CDELT1"+s])
     db = float(header["CDELT2"+s])
     x0 = float(header["CRPIX1"+s])
diff --git a/scopesim/optics/optical_element.py b/scopesim/optics/optical_element.py
index f2b46545..7f630c49 100644
--- a/scopesim/optics/optical_element.py
+++ b/scopesim/optics/optical_element.py
@@ -215,6 +215,13 @@ def __repr__(self):
     def __str__(self):
         return f"{self.__class__.__name__}: \"{self.display_name}\""
 
+    def _repr_pretty_(self, p, cycle):
+        """For ipython"""
+        if cycle:
+            p.text(f"{self.__class__.__name__}(...)")
+        else:
+            p.text(str(self))
+
     @property
     def properties_str(self):
         prop_str = ""
diff --git a/scopesim/optics/optical_train.py b/scopesim/optics/optical_train.py
index d0ed34dc..65178279 100644
--- a/scopesim/optics/optical_train.py
+++ b/scopesim/optics/optical_train.py
@@ -497,6 +497,36 @@ def __repr__(self):
     def __str__(self):
         return self._description
 
+    def _repr_pretty_(self, p, cycle):
+        """For ipython"""
+        if cycle:
+            p.text(f"{self.__class__.__name__}(...)")
+        else:
+            p.text(f"{self.__class__.__name__} ")
+            p.text(f"for {self.cmds['!OBS.instrument']} ")
+            p.text(f"@ {self.cmds['!TEL.telescope']}:")
+            p.breakable()
+            p.text("UserCommands:")
+            p.breakable()
+            p.pretty(self.cmds)
+            p.breakable()
+            p.text("OpticalElements:")
+            with p.indent(2):
+                for item in self:
+                    p.breakable()
+                    p.pretty(item)
+            p.breakable()
+            p.text("DetectorArrays:")
+            with p.indent(2):
+                for item in self.detector_arrays:
+                    p.breakable()
+                    p.pretty(item)
+            p.breakable()
+            p.text("Effects:")
+            p.breakable()
+            with p.indent(2):
+                p.pretty(self.effects)
+
     def __getitem__(self, item):
         return self.optics_manager[item]
 
diff --git a/scopesim/source/source.py b/scopesim/source/source.py
index 23ac7a23..de9248e3 100644
--- a/scopesim/source/source.py
+++ b/scopesim/source/source.py
@@ -55,6 +55,8 @@
 
 from ..base_classes import SourceBase
 from .. import utils
+from ..utils import close_loop, figure_factory
+# TODO: explicit util imports above...
 
 
 class Source(SourceBase):
@@ -529,23 +531,21 @@ def plot(self):
         Source components instantiated from 2d or 3d ImageHDUs are represented by their
         spatial footprint. Source components instantiated from tables are shown as points.
         """
-        # pylint: disable=import-outside-toplevel
-        import matplotlib.pyplot as plt
+        _, axes = figure_factory()
 
         colours = "rgbcymk" * (len(self.fields) // 7 + 1)
         for col, field in zip(colours, self.fields):
             if isinstance(field, Table):
-                plt.plot(field["x"], field["y"], col+".")
+                axes.plot(field["x"], field["y"], col+".")
             elif isinstance(field, (fits.ImageHDU, fits.PrimaryHDU)):
                 xpts, ypts = imp_utils.calc_footprint(field.header)
-                xpts *= 3600   # Because ImageHDUs are always in CUNIT=DEG
-                ypts *= 3600
-                xpts = list(xpts) + [xpts[0]]
-                ypts = list(ypts) + [ypts[0]]
-                plt.plot(xpts, ypts, col)
-                plt.xlabel("x [arcsec]")
-                plt.ylabel("y [arcsec]")
-        plt.gca().set_aspect("equal")
+                # * 3600, because ImageHDUs are always in CUNIT=DEG
+                xpts = list(close_loop(xpts * 3600))
+                ypts = list(close_loop(ypts * 3600))
+                axes.plot(xpts, ypts, col)
+                axes.set_xlabel("x [arcsec]")
+                axes.set_ylabel("y [arcsec]")
+        axes.set_aspect("equal")
 
     def make_copy(self):
         new_source = Source()
diff --git a/scopesim/utils.py b/scopesim/utils.py
index eafb6f82..d653b40c 100644
--- a/scopesim/utils.py
+++ b/scopesim/utils.py
@@ -6,12 +6,15 @@
 import sys
 import logging
 from collections import OrderedDict
-from docutils.core import publish_string
+from collections.abc import Iterable, Generator
+from itertools import chain
 from copy import deepcopy
 
+from docutils.core import publish_string
 import requests
 import yaml
 import numpy as np
+from matplotlib import pyplot as plt
 from astropy import units as u
 from astropy.io import fits
 from astropy.table import Column, Table
@@ -178,7 +181,7 @@ def deriv_polynomial2d(poly):
     degree = poly.degree
     dpoly_dx = Polynomial2D(degree=degree - 1)
     dpoly_dy = Polynomial2D(degree=degree - 1)
-    regexp = re.compile(r'c(\d+)_(\d+)')
+    regexp = re.compile(r"c(\d+)_(\d+)")
     for pname in poly.param_names:
         # analyse the name
         match = regexp.match(pname)
@@ -261,8 +264,7 @@ def seq(start, stop, step=1):
     except ZeroDivisionError:
         if step == 0 and delta == 0:
             return start
-        else:
-            raise ValueError("invalid '(stop - start) / step'")
+        raise ValueError("invalid '(stop - start) / step'")
 
     if npts < 0:
         raise ValueError("wrong sign in 'step' argument")
@@ -338,7 +340,8 @@ def telescope_diffraction_limit(aperture_size, wavelength, distance=None):
 
     """
 
-    diff_limit = (((wavelength * u.um) / (aperture_size * u.m)) * u.rad).to(u.arcsec).value
+    diff_limit = (((wavelength * u.um) / (aperture_size * u.m))
+                  * u.rad).to(u.arcsec).value
 
     if distance is not None:
         diff_limit *= distance / u.pc.to(u.AU)
@@ -598,9 +601,9 @@ def change_table_entry(tbl, col_name, new_val, old_val=None, position=None):
 
 
 def real_colname(name, colnames, silent=True):
-    names = [name.lower(), name.upper(), name[0].upper() + name[1:].lower()]
+    names = [name.lower(), name.upper(), name.capitalize()]
     real_name = [name for name in names if name in colnames]
-    if len(real_name) == 0:
+    if not real_name:
         real_name = None
         if not silent:
             logging.warning("None of %s were found in %s", names, colnames)
@@ -654,13 +657,13 @@ def get_meta_quantity(meta_dict, name, fallback_unit=""):
 
     """
 
-    if isinstance(meta_dict[name], str) and meta_dict[name][0] == "!":
+    if isinstance(meta_dict[name], str) and meta_dict[name].startswith("!"):
         meta_dict[name] = from_currsys(meta_dict[name])
 
     if isinstance(meta_dict[name], u.Quantity):
         unit = meta_dict[name].unit
-    elif name + "_unit" in meta_dict:
-        unit = meta_dict[name + "_unit"]
+    elif f"{name}_unit" in meta_dict:
+        unit = meta_dict[f"{name}_unit"]
     else:
         unit = u.Unit(fallback_unit)
 
@@ -684,7 +687,7 @@ def quantify(item, unit):
 
     """
 
-    if isinstance(item, str) and item[0] == "!":
+    if isinstance(item, str) and item.startswith("!"):
         item = from_currsys(item)
     if isinstance(item, u.Quantity):
         quant = item.to(u.Unit(unit))
@@ -752,13 +755,9 @@ def extract_base_from_unit(unit, base_unit):
     return new_unit, extracted_units
 
 
-def is_fits(filename):
-    flag = False
-    if filename is not None:
-        if filename.split(".")[-1].lower() in "fits":
-            flag = True
-
-    return flag
+def is_fits(filename) -> bool:
+    # Using 'in ".fits"' to also catch ".fit", which exists sometimes...
+    return (filename is not None and Path(filename).suffix.lower() in ".fits")
 
 
 def get_fits_type(filename):
@@ -779,7 +778,7 @@ def quantity_from_table(colname, table, default_unit=""):
         else:
             col = col.data << col.unit
     else:
-        colname_u = colname + "_unit"
+        colname_u = f"{colname}_unit"
         if colname_u in table.meta:
             col = col * u.Unit(table.meta[colname_u])
         else:
@@ -802,7 +801,7 @@ def unit_from_table(colname, table, default_unit=""):
     """
     Looks for the unit for a column based on the meta dict keyword "<col>_unit"
     """
-    colname_u = colname + "_unit"
+    colname_u = f"{colname}_unit"
     col = table[colname]
     if col.unit is not None:
         unit = col.unit
@@ -833,16 +832,16 @@ def has_needed_keywords(header, suffix=""):
     """
     Check to see if the WCS keywords are in the header
     """
-    keys = ["CDELT1", "CRVAL1", "CRPIX1"]
-    return (sum(key + suffix in header.keys() for key in keys) == 3 and
-            "NAXIS1" in header.keys())
+    keys = {"CDELT1", "CRVAL1", "CRPIX1"}
+    keys = {key + suffix for key in keys}
+    keys.add("NAXIS1")
+    return all(key in header.keys() for key in keys)
 
 
 def stringify_dict(dic, ignore_types=(str, int, float)):
     """
     Turns a dict entries into strings for addition to FITS headers
     """
-    from copy import deepcopy
     dic_new = deepcopy(dic)
     for key in dic_new:
         if not isinstance(dic_new[key], ignore_types):
@@ -925,7 +924,8 @@ def check_keys(input_dict, required_keys, action="error", all_any="all"):
     if not keys_present:
         if "error" in action:
             raise ValueError("One or more of the following keys missing from "
-                             f"input_dict: \n{required_keys} \n{input_dict.keys()}")
+                             f"input_dict: \n{required_keys} "
+                             f"\n{input_dict.keys()}")
         if "warn" in action:
             logging.warning(("One or more of the following keys missing "
                              "from input_dict: \n%s \n%s"), required_keys,
@@ -976,35 +976,71 @@ def pretty_print_dict(dic, indent=0):
     return text
 
 
-def return_latest_github_actions_jobs_status(owner_name="AstarVienna", repo_name="ScopeSim",
-                                             branch="dev_master", actions_yaml_name="tests.yml"):
+def return_latest_github_actions_jobs_status(
+        owner_name="AstarVienna",
+        repo_name="ScopeSim",
+        branch="dev_master",
+        actions_yaml_name="tests.yml",
+    ):
     """
     Gets the status of the latest test run
     """
-    response = requests.get(f"https://api.github.com/repos/{owner_name}/{repo_name}/"
-                            f"actions/workflows/{actions_yaml_name}/runs?branch={branch}&per_page=1")
+    response = requests.get(
+        f"https://api.github.com/repos/{owner_name}/{repo_name}/actions/"
+        f"workflows/{actions_yaml_name}/runs?branch={branch}&per_page=1"
+    )
     dic = response.json()
     run_id = dic["workflow_runs"][0]["id"]
 
-    response = requests.get(f"https://api.github.com/repos/{owner_name}/{repo_name}/actions/runs/{run_id}/jobs")
+    response = requests.get(
+        f"https://api.github.com/repos/{owner_name}/{repo_name}/actions/runs/"
+        f"{run_id}/jobs"
+    )
     dic = response.json()
     params_list = []
     for job in dic["jobs"]:
         params = {
-            "name": job['name'],
-            "status": job['status'],
-            "conclusion": job['conclusion'],
-            "started_at": job['started_at'],
-            "completed_at": job['completed_at'],
-            "url": job['html_url'],
+            "name": job["name"],
+            "status": job["status"],
+            "conclusion": job["conclusion"],
+            "started_at": job["started_at"],
+            "completed_at": job["completed_at"],
+            "url": job["html_url"],
             "badge_url": None
         }
 
-        key = "Python_" + job['name'].split()[-1][:-1]
-        value = "passing" if job['conclusion'] == "success" else "failing"
-        colour = "brightgreen" if job['conclusion'] == "success" else "red"
+        # TODO: this could use the new badges from IRDB, once that's in
+        #       scopesim_core...
+        key = "Python_" + job["name"].split()[-1][:-1]
+        value = "passing" if job["conclusion"] == "success" else "failing"
+        colour = "brightgreen" if job["conclusion"] == "success" else "red"
         badge_url = f"https://img.shields.io/badge/{key}-{value}-{colour}"
         params["badge_url"] = badge_url
         params_list.append(params)
 
     return params_list
+
+
+def close_loop(iterable: Iterable) -> Generator:
+    """x, y = zip(*close_loop(zip(x, y)))"""
+    iterator = iter(iterable)
+    first = next(iterator)
+    yield first
+    yield from iterator
+    yield first
+
+
+def figure_factory(nrows=1, ncols=1, **kwargs):
+    """Default way to init fig and ax, to easily modify later."""
+    iterable_axes = kwargs.pop("iterable_axes", False)
+    fig, ax = plt.subplots(nrows, ncols, **kwargs)
+    if iterable_axes and not isinstance(ax, Iterable):
+        ax = (ax,)
+    return fig, ax
+
+
+def figure_grid_factory(nrows=1, ncols=1, **kwargs):
+    """Gridspec variant."""
+    fig = plt.figure()
+    gs = fig.add_gridspec(nrows, ncols, **kwargs)
+    return fig, gs
diff --git a/scopesim/version.py b/scopesim/version.py
index 844b2eba..03f73be0 100644
--- a/scopesim/version.py
+++ b/scopesim/version.py
@@ -1,7 +1,15 @@
 from importlib import metadata
 version = metadata.version(__package__)
-date = '2023-07-10 10:00:00 GMT'
+date = '2023-09-06 13:00:00 GMT'
 yaml_descriptions = """
+- version : 0.6.1
+  date : 2023-09-06
+  comment : Patch with visualisation improvements and general refactoring
+  changes :
+  - Improvements to console representation and plot methods of various classes: #252, #260, #263, #266
+  - Refactoring and changes to the inheritance of some classes: #261, #264, #265
+  - Changes to the CI configuration to reduce fails caused by web request timeouts #255, #262
+
 - version : 0.6.0
   date : 2023-07-10
   comment : Summer 2023


Test patch:
diff --git a/scopesim/tests/tests_commands/test_UserCommands.py b/scopesim/tests/tests_commands/test_UserCommands.py
index aeeaafa6..5f90bfe5 100644
--- a/scopesim/tests/tests_commands/test_UserCommands.py
+++ b/scopesim/tests/tests_commands/test_UserCommands.py
@@ -87,6 +87,14 @@ def test_mode_yamls_read_from_file(self):
         assert cmd["!OBS.airmass"] == 2
         assert cmd.yaml_dicts[-1]["effects"][0]["kwargs"]["meaning_of_life"] == 42
 
+    def test_init_through_repr(self):
+        """Check whether we can recreate a UserCommand by evaluating its __repr__."""
+        cmd1 = UserCommands(use_instrument="test_package")
+        cmd2 = eval(repr(cmd1))
+        # TODO: Create a proper __eq__ so we can assert cmd1 == cmd2
+        assert str(cmd1) == str(cmd2)
+        assert cmd1.cmds == cmd2.cmds
+
 
 class TestMiscFeatures:
     def test_updates_with_yaml_dict(self):
@@ -101,6 +109,11 @@ def test_update_works_via_setitem(self):
         cmd["!TEL.gigawatts"] = 1.21
         assert cmd["!TEL.gigawatts"] == 1.21
 
+    def test_str(self):
+        """Test whether __str__ gives a pretty result."""
+        cmd = UserCommands(use_instrument="test_package")
+        assert "├─" in str(cmd)
+
 
 class TestListLocalPackages:
     def test_all_packages_listed(self):
diff --git a/scopesim/tests/tests_effects/test_SpectralTraceList.py b/scopesim/tests/tests_effects/test_SpectralTraceList.py
index 535e610a..17aa38cc 100644
--- a/scopesim/tests/tests_effects/test_SpectralTraceList.py
+++ b/scopesim/tests/tests_effects/test_SpectralTraceList.py
@@ -5,7 +5,8 @@
 from astropy.io import fits
 
 
-from scopesim.effects.spectral_trace_list import SpectralTraceList
+from scopesim.effects.spectral_trace_list import SpectralTraceList, \
+    SpectralTraceListWheel
 from scopesim.effects.spectral_trace_list_utils import SpectralTrace
 from scopesim.tests.mocks.py_objects import trace_list_objects as tlo
 from scopesim.tests.mocks.py_objects import header_objects as ho
@@ -21,19 +22,23 @@
 # pylint: disable=missing-class-docstring,
 # pylint: disable=missing-function-docstring
 
+
 @pytest.fixture(name="slit_header", scope="class")
 def fixture_slit_header():
     return ho._short_micado_slit_header()
 
+
 @pytest.fixture(name="long_slit_header", scope="class")
 def fixture_long_slit_header():
     return ho._long_micado_slit_header()
 
+
 @pytest.fixture(name="full_trace_list", scope="class")
 def fixture_full_trace_list():
     """Instantiate a trace definition hdu list"""
     return tlo.make_trace_hdulist()
 
+
 class TestInit:
     def test_initialises_with_nothing(self):
         assert isinstance(SpectralTraceList(), SpectralTraceList)
@@ -70,13 +75,32 @@ def fixture_spectral_trace_list():
     """Instantiate a SpectralTraceList"""
     return SpectralTraceList(hdulist=tlo.make_trace_hdulist())
 
+
 class TestRectification:
     def test_rectify_cube_not_implemented(self, spectral_trace_list):
         hdulist = fits.HDUList()
         with pytest.raises(NotImplementedError):
             spectral_trace_list.rectify_cube(hdulist)
 
-    #def test_rectify_traces_needs_ximin_and_ximax(self, spectral_trace_list):
+    # def test_rectify_traces_needs_ximin_and_ximax(self, spectral_trace_list):
     #    hdulist = fits.HDUList([fits.PrimaryHDU()])
     #    with pytest.raises(KeyError):
     #        spectral_trace_list.rectify_traces(hdulist)
+
+
+class TestSpectralTraceListWheel:
+    def test_basic_init(self):
+        """
+        This is a super basic test just to see the thing basically works and
+        parameters are passed correctly. Please feel free to improve this!!
+        """
+        kwargs = {"current_trace_list": "bogus",
+                  "filename_format": "bogus_{}",
+                  "trace_list_names": ["foo"]}
+        stw = SpectralTraceListWheel(**kwargs)
+        assert isinstance(stw, SpectralTraceListWheel)
+        assert stw.meta["current_trace_list"] == "bogus"
+        assert stw.meta["filename_format"] == "bogus_{}"
+        assert stw.meta["trace_list_names"] == ["foo"]
+        assert isinstance(stw.trace_lists["foo"], SpectralTraceList)
+        assert stw.trace_lists["foo"].meta["filename"] == "bogus_foo"
diff --git a/scopesim/tests/tests_effects/test_TERCurve.py b/scopesim/tests/tests_effects/test_TERCurve.py
index 14fe65e7..95ba57bc 100644
--- a/scopesim/tests/tests_effects/test_TERCurve.py
+++ b/scopesim/tests/tests_effects/test_TERCurve.py
@@ -158,7 +158,7 @@ def test_returns_filter_as_wanted(self, observatory, instrument,
                                              current_filter=default_filter,
                                              name="test_svo_wheel")
 
-        assert isinstance(filt_wheel, tc.FilterWheel)
+        assert isinstance(filt_wheel, tc.FilterWheelBase)
         assert default_filter in filt_wheel.filters
 
     def test_returns_filters_with_include_str(self):
diff --git a/scopesim/tests/tests_optics/test_FieldOfView.py b/scopesim/tests/tests_optics/test_FieldOfView.py
index 99ff7152..45e571f4 100644
--- a/scopesim/tests/tests_optics/test_FieldOfView.py
+++ b/scopesim/tests/tests_optics/test_FieldOfView.py
@@ -140,8 +140,9 @@ def test_contains_all_fields_inside_fov(self, basic_fov_header,
             assert isinstance(the_fov.fields[2], Table)
 
 
-@pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ")
+# @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ")
 class TestMakeCube:
+    @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ")
     def test_makes_cube_from_table(self):
         src_table = so._table_source()            # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm
         fov = _fov_190_210_um()
@@ -163,6 +164,7 @@ def test_makes_cube_from_table(self):
             plt.imshow(cube.data[0, :, :], origin="lower")
             plt.show()
 
+    @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ")
     def test_makes_cube_from_imagehdu(self):
         src_image = so._image_source()            # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm
         fov = _fov_190_210_um()
@@ -181,6 +183,7 @@ def test_makes_cube_from_imagehdu(self):
             plt.imshow(cube.data[0, :, :], origin="lower")
             plt.show()
 
+    @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ")
     def test_makes_cube_from_other_cube_imagehdu(self):
         import scopesim as sim
         sim.rc.__currsys__["!SIM.spectral.spectral_bin_width"] = 0.01
@@ -201,6 +204,7 @@ def test_makes_cube_from_other_cube_imagehdu(self):
             plt.imshow(cube.data[0, :, :], origin="lower")
             plt.show()
 
+    @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ")
     def test_makes_cube_from_two_similar_cube_imagehdus(self):
         src_cube = so._cube_source() + so._cube_source(dx=1)            # 2 cubes 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm
         fov = _fov_197_202_um()
@@ -220,6 +224,7 @@ def test_makes_cube_from_two_similar_cube_imagehdus(self):
             plt.imshow(cube.data[0, :, :], origin="lower")
             plt.show()
 
+    @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ")
     def test_makes_cube_from_all_types_of_source_object(self):
         src_all = so._table_source() + \
                   so._image_source(dx=-4, dy=-4) + \
@@ -271,7 +276,7 @@ def test_cube_has_full_wcs(self, src):
         assert "CTYPE3" in cube.header
 
 
-@pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset")
+# @pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset")
 class TestMakeImage:
     def test_makes_image_from_table(self):
         src_table = so._table_source()            # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm
@@ -317,6 +322,7 @@ def test_makes_image_from_image(self):
             plt.imshow(img.data, origin="lower")
             plt.show()
 
+    @pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset")
     def test_makes_image_from_cube(self):
         src_cube = so._cube_source()  # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm
         fov = _fov_197_202_um()
@@ -334,6 +340,7 @@ def test_makes_image_from_cube(self):
 
         assert out_sum == approx(in_sum)
 
+    @pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset")
     def test_makes_image_from_all_types_of_source_object(self):
         src_all = so._table_source() + \
                   so._image_source(dx=-4, dy=-4) + \
@@ -370,7 +377,7 @@ def test_makes_image_from_all_types_of_source_object(self):
             plt.show()
 
 
-@pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset")
+# @pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset")
 class TestMakeSpectrum:
     def test_make_spectrum_from_table(self):
         src_table = so._table_source()            # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm
diff --git a/scopesim/tests/tests_optics/test_fov_manager_utils.py b/scopesim/tests/tests_optics/test_fov_manager_utils.py
index 1007af7e..63468858 100644
--- a/scopesim/tests/tests_optics/test_fov_manager_utils.py
+++ b/scopesim/tests/tests_optics/test_fov_manager_utils.py
@@ -13,11 +13,31 @@
 from scopesim.tests.mocks.py_objects import trace_list_objects as tlo
 
 
+@pytest.fixture(scope="function")
+def wave_kwargs():
+    return {"wave_min": 0.5, "wave_max": 2.5}
+
+
+@pytest.fixture(scope="function")
+def th_filt():
+    return eo._filter_tophat_curve()
+
+
 @pytest.fixture(scope="function")
 def full_trace_list():
     return tlo.make_trace_hdulist()
 
 
+@pytest.fixture(scope="function")
+def spec_hdrs(full_trace_list):
+    params = {"pixel_scale": 0.1, "plate_scale": 0.1,
+              "wave_min": 0.7, "wave_max": 2.5}
+    spt = SpectralTraceList(hdulist=full_trace_list, **params)
+    apm = apo._basic_aperture()
+    hdrs = fm_utils.get_spectroscopy_headers(effects=[spt, apm], **params)
+    return list(hdrs)
+
+
 PLOTS = False
 
 
@@ -73,35 +93,35 @@ def test_combined_shifts_reduced_to_usable_number(self, sub_pix_frac):
 
 
 class TestGetImagingWaveset:
-    def test_returns_default_wave_range_when_passed_no_effects(self):
-        kwargs = {"wave_min": 0.5, "wave_max": 2.5}
-        wave_bin_edges = fm_utils.get_imaging_waveset([], **kwargs)
+    @pytest.mark.usefixtures("wave_kwargs")
+    def test_returns_default_wave_range_when_passed_no_effects(self, wave_kwargs):
+        wave_bin_edges = fm_utils.get_imaging_waveset([], **wave_kwargs)
         assert len(wave_bin_edges) == 2
 
-    def test_returns_waveset_of_filter(self):
-        filt = eo._filter_tophat_curve()
-        kwargs = {"wave_min": 0.5, "wave_max": 2.5}
-        wave_bin_edges = fm_utils.get_imaging_waveset([filt], **kwargs)
+    @pytest.mark.usefixtures("wave_kwargs", "th_filt")
+    def test_returns_waveset_of_filter(self, wave_kwargs, th_filt):
+        wave_bin_edges = fm_utils.get_imaging_waveset([th_filt], **wave_kwargs)
         assert len(wave_bin_edges) == 2
 
-    def test_returns_waveset_of_psf(self):
+    @pytest.mark.usefixtures("wave_kwargs")
+    def test_returns_waveset_of_psf(self, wave_kwargs):
         psf = eo._const_psf()
-        kwargs = {"wave_min": 0.5, "wave_max": 2.5}
-        wave_bin_edges = fm_utils.get_imaging_waveset([psf], **kwargs)
+        wave_bin_edges = fm_utils.get_imaging_waveset([psf], **wave_kwargs)
         assert len(wave_bin_edges) == 4
 
-    def test_returns_waveset_of_psf_and_filter(self):
-        filt = eo._filter_tophat_curve()
+    @pytest.mark.usefixtures("wave_kwargs", "th_filt")
+    def test_returns_waveset_of_psf_and_filter(self, wave_kwargs, th_filt):
         psf = eo._const_psf()
-        kwargs = {"wave_min": 0.5, "wave_max": 2.5}
-        wave_bin_edges = fm_utils.get_imaging_waveset([filt, psf], **kwargs)
+        wave_bin_edges = fm_utils.get_imaging_waveset([th_filt, psf],
+                                                      **wave_kwargs)
         assert len(wave_bin_edges) == 4
 
-    def test_returns_waveset_of_ncpa_psf_inside_filter_edges(self):
-        filt = eo._filter_tophat_curve()
+    @pytest.mark.usefixtures("wave_kwargs", "th_filt")
+    def test_returns_waveset_of_ncpa_psf_inside_filter_edges(self, wave_kwargs,
+                                                             th_filt):
         psf = eo._ncpa_psf()
-        kwargs = {"wave_min": 0.5, "wave_max": 2.5}
-        wave_bin_edges = fm_utils.get_imaging_waveset([psf, filt], **kwargs)
+        wave_bin_edges = fm_utils.get_imaging_waveset([psf, th_filt],
+                                                      **wave_kwargs)
         assert min(wave_bin_edges) == 1.
         assert max(wave_bin_edges) == 2.
         assert len(wave_bin_edges) == 9
@@ -111,7 +131,7 @@ class TestGetImagingHeaders:
     def test_throws_error_if_not_all_keywords_are_passed(self):
         apm = eo._img_aperture_mask()
         with pytest.raises(ValueError):
-            fm_utils.get_imaging_headers([apm])
+            list(fm_utils.get_imaging_headers([apm]))
 
     def test_returns_set_of_headers_from_aperture_effects(self):
         apm = eo._img_aperture_mask(array_dict={"x": [-1.28, 1., 1., -1.28],
@@ -129,6 +149,8 @@ def test_returns_set_of_headers_from_detector_list_effect(self):
         kwargs = {"pixel_scale": 0.004, "plate_scale": 0.26666666666,
                   "max_segment_size": 2048 ** 2, "chunk_size": 1024}
         hdrs = fm_utils.get_imaging_headers([det], **kwargs)
+        # Evaluate generator for testing
+        hdrs = list(hdrs)
 
         area_sum = np.sum([hdr["NAXIS1"] * hdr["NAXIS2"] for hdr in hdrs])
         assert area_sum == 4096**2
@@ -161,6 +183,8 @@ def test_returns_fov_objects_for_basic_input(self):
                   "max_segment_size": 100 ** 2, "chunk_size": 100}
 
         hdrs = fm_utils.get_imaging_headers([apm], **kwargs)
+        # Evaluate generator for testing
+        hdrs = list(hdrs)
         waveset = np.linspace(1, 2, 6)
         shifts = {"wavelengths": np.array([1, 2]),
                   "x_shifts": np.zeros(2),
@@ -168,7 +192,10 @@ def test_returns_fov_objects_for_basic_input(self):
         fovs = fm_utils.get_imaging_fovs(headers=hdrs, waveset=waveset,
                                          shifts=shifts)
 
+        # Evaluate generator for testing
+        fovs = list(fovs)
         assert len(fovs) == (len(waveset)-1) * len(hdrs)
+        assert fovs
 
         if PLOTS:
             from scopesim.optics.image_plane_utils import calc_footprint
@@ -189,19 +216,13 @@ def test_returns_fov_objects_for_basic_input(self):
             plt.show()
 
 
-@pytest.mark.usefixtures("full_trace_list")
+@pytest.mark.usefixtures("spec_hdrs")
 class TestGetSpectroscopyHeaders:
-    def test_returns_headers(self, full_trace_list):
-        params = {"pixel_scale": 0.1, "plate_scale": 0.1,
-                  "wave_min": 0.7, "wave_max": 2.5}
-        spt = SpectralTraceList(hdulist=full_trace_list, **params)
-        apm = apo._basic_aperture()
-
-        hdrs = fm_utils.get_spectroscopy_headers(effects=[spt, apm], **params)
-        assert all([isinstance(hdr, PoorMansHeader) for hdr in hdrs])
+    def test_returns_headers(self, spec_hdrs):
+        assert all([isinstance(hdr, PoorMansHeader) for hdr in spec_hdrs])
 
         if PLOTS:
-            for hdr in hdrs:
+            for hdr in spec_hdrs:
                 x = np.array([0, hdr["NAXIS1"], hdr["NAXIS1"], 0])
                 y = np.array([0, 0, hdr["NAXIS2"], hdr["NAXIS2"]])
                 xw, yw = pix2val(hdr, x, y, "D")
@@ -210,19 +231,15 @@ def test_returns_headers(self, full_trace_list):
             plt.show()
 
 
-@pytest.mark.usefixtures("full_trace_list")
+@pytest.mark.usefixtures("spec_hdrs")
 class TestGetSpectroscopyFOVs:
-    def test_returns_fovs(self, full_trace_list):
-        params = {"pixel_scale": 0.1, "plate_scale": 0.1,
-                  "wave_min": 0.7, "wave_max": 2.5}
-        spt = SpectralTraceList(hdulist=full_trace_list, **params)
-        apm = apo._basic_aperture()
-
+    def test_returns_fovs(self, spec_hdrs):
         shifts = {"wavelengths": np.array([0.7, 2.5]),
                   "x_shifts": np.array([0, 0]),
                   "y_shifts": np.array([0, 1/3600.])}
-        hdrs = fm_utils.get_spectroscopy_headers(effects=[spt, apm], **params)
-        fovs = fm_utils.get_spectroscopy_fovs(hdrs, shifts)
+        fovs = fm_utils.get_spectroscopy_fovs(spec_hdrs, shifts)
+        # Evaluate generator for testing
+        fovs = list(fovs)
 
         assert all([isinstance(fov, FieldOfView) for fov in fovs])
 
diff --git a/scopesim/tests/tests_server/test_database.py b/scopesim/tests/tests_server/test_database.py
index b0bbbc46..73e0fa9a 100644
--- a/scopesim/tests/tests_server/test_database.py
+++ b/scopesim/tests/tests_server/test_database.py
@@ -178,6 +178,7 @@ def test_downloads_specific_package(self):
 
             assert version_dict["version"] == release
 
+    @pytest.mark.skip(reason="fails too often with timeout")
     @pytest.mark.webtest
     def test_downloads_github_version_of_package_with_semicolon(self):
         release = "github:728761fc76adb548696205139e4e9a4260401dfc"
@@ -188,6 +189,7 @@ def test_downloads_github_version_of_package_with_semicolon(self):
 
             assert os.path.exists(filename)
 
+    @pytest.mark.skip(reason="fails too often with timeout")
     @pytest.mark.webtest
     def test_downloads_github_version_of_package_with_at_symbol(self):
         release = "github@728761fc76adb548696205139e4e9a4260401dfc"
@@ -199,6 +201,7 @@ def test_downloads_github_version_of_package_with_at_symbol(self):
             assert os.path.exists(filename)
 
 
+@pytest.mark.skip(reason="fails too often with timeout")
 class TestDownloadGithubFolder:
     @pytest.mark.webtest
     def test_downloads_current_package(self):