diff --git a/scopesim/effects/apertures.py b/scopesim/effects/apertures.py
index bc8c50ec..756e7860 100644
--- a/scopesim/effects/apertures.py
+++ b/scopesim/effects/apertures.py
@@ -128,6 +128,7 @@ def apply_to(self, obj, **kwargs):
             for vol in obj.volumes:
                 vol["meta"]["xi_min"] = min(x) * u.arcsec
                 vol["meta"]["xi_max"] = max(x) * u.arcsec
+                vol["meta"]["slit_name"] = self.meta["name"]
 
         return obj
 
diff --git a/scopesim/effects/selector_wheel.py b/scopesim/effects/selector_wheel.py
index 3bebfe17..8eabbcb6 100644
--- a/scopesim/effects/selector_wheel.py
+++ b/scopesim/effects/selector_wheel.py
@@ -49,14 +49,12 @@ class SelectorWheel(Effect):
 
     """
 
-    z_order = (290, 690, 890)
+    z_order = ()
     required_keys = {"selector_key", "wheel"}
 
     def __init__(self, **kwargs):
-        super().__init__(**kwargs)
         check_keys(kwargs, self.required_keys, action="error")
-
-        self.meta.update(kwargs)
+        super().__init__(**kwargs)
 
         self.wheel_effects = {}
         for wheel_entry in self.meta["wheel"]:
@@ -70,9 +68,12 @@ def __init__(self, **kwargs):
             # Instantiate the effect and store it in the wheel_effects dictionary
             if isinstance(selector_value, list):
                 for val in selector_value:
-                    self.wheel_effects[val] = effect_class(**effect_kwargs)
+                    self.wheel_effects[val] = effect_class(cmds=self.cmds, **effect_kwargs)
             else:
-                self.wheel_effects[selector_value] = effect_class(**effect_kwargs)
+                self.wheel_effects[selector_value] = effect_class(cmds=self.cmds, **effect_kwargs)
+
+        # Use the wheel effects' z_order as the z_order of the selector wheel
+        self.z_order = [eff.z_order for eff in self.wheel_effects.values()][0]
 
 
     def apply_to(self, obj, **kwargs):
@@ -95,15 +96,19 @@ def apply_to(self, obj, **kwargs):
 
             for val in unique_selector_values:
                 vols_with_val = [vol for vol in obj.volumes if vol["meta"].get(self.meta["selector_key"], None) == val]
-                effect_to_apply = self.get_effect(val)
-                logger.info(f"Applying effect for selector value: {val} -> {effect_to_apply}, volumes: {len(vols_with_val)}")
 
-                if val is None or effect_to_apply is None:
+                if val is None:
+                    logger.warning(f"Volume(s) with missing selector key {self.meta['selector_key']} value found, "
+                                   f"applying no effect to those volumes.")
                     new_volumes.extend(vols_with_val)
                     continue
 
-                #TODO: If effect_to_apply is a dichroic which reassigns selector_key values, we need to add a check here
-                #TODO: i.e. dichroic.arm_action.keys() should not include items in unique_selector_values other than val
+                effect_to_apply = self.get_effect(val)
+                logger.debug(f"Applying effect for {self.meta['selector_key']}: {val} -> {effect_to_apply}, volumes: {len(vols_with_val)}")
+
+                if effect_to_apply is None:
+                    new_volumes.extend(vols_with_val)
+                    continue
 
                 newvollist = FovVolumeList()
                 newvollist.volumes = vols_with_val
@@ -113,7 +118,7 @@ def apply_to(self, obj, **kwargs):
             obj.volumes = new_volumes
 
         if isinstance(obj, Detector):
-            logger.info("Since passed object is a Detector, selector_key by default is the ID of the Detector object.")
+            logger.debug("Since passed object is a Detector, selector_key by default is the ID of the Detector object.")
             selector_value = obj.meta[real_colname("id", obj.meta)] # Assuming detector ID is the selector
 
             effect_to_apply = self.get_effect(selector_value)
@@ -128,8 +133,9 @@ def apply_to(self, obj, **kwargs):
 
     def get_effect(self, selector_value):
         eff = None
-        if (selector_value is None) or (selector_value not in self.wheel_effects.keys()):
-            logger.warning(f"Either None or Missing value of {self.meta["selector_key"]} requested: {selector_value}, returning None.")
+        if selector_value not in self.wheel_effects.keys():
+            logger.warning(f"Entry for selector value {selector_value} not found in wheel effects. "
+                           f"Assuming no effect to apply for this selector value.")
         else:
             eff = self.wheel_effects[selector_value]
         return eff
diff --git a/scopesim/effects/spectral_efficiency.py b/scopesim/effects/spectral_efficiency.py
index 35fe76fc..7a254785 100644
--- a/scopesim/effects/spectral_efficiency.py
+++ b/scopesim/effects/spectral_efficiency.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 """Spectral grating efficiencies."""
 
-from typing import ClassVar
+from typing import ClassVar, Callable
 
 import numpy as np
 from astropy.io import fits
@@ -12,7 +12,9 @@
 from .effects import Effect
 from .ter_curves import TERCurve
 from .ter_curves_utils import apply_throughput_to_cube
-from ..utils import figure_factory, get_logger
+from ..utils import figure_factory, get_logger, check_keys
+from .data_container import DataContainer
+from ..optics import echelle
 
 
 logger = get_logger(__name__)
@@ -127,3 +129,87 @@ def plot(self):
         axes.legend()
 
         return fig
+
+
+class EchelleSpectralEfficiency(Effect):
+    """
+    Spectral efficiency list from analytical calculations of the blaze function for ZShooter gratings.
+    Requires same input trace parameter table as EchelleSpectralTraceList, supply as kwarg "filename"
+    """
+    z_order: ClassVar[tuple[int, ...]] = (630,)
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self.efficiency_generator = self._generate_efficiency_curve_func()
+        self.efficiencies = {}
+
+    def _generate_efficiency_curve_func(self) -> Callable:
+        trace_params = self.table
+        spectrographs = {}
+        for row in trace_params:
+            prefix = row["prefix"]  # note trance ids are assumed to be prefix_{order}
+            min_order = row['m0'] - row['n']
+            max_order = row['m0']
+            min_wave = row['min_wave'] * u.Unit(trace_params.meta["min_wave_unit"])
+            max_wave = row['max_wave'] * u.Unit(trace_params.meta["max_wave_unit"])
+            design_res = row['design_res']
+            focal_len = row['focal_length'] * u.Unit(trace_params.meta["focal_length_unit"])
+            disp_npix = row['n_disp'] - 2 * row['detector_pad']
+            xdisp_npix = row['n_xdisp']- 2 * row['detector_pad']
+            pix_size = row['pixel_size'] * u.Unit(trace_params.meta["pixel_size_unit"])
+            echelle_angle = np.deg2rad(row['echelle_blaze'])*u.rad
+            xdisp_beta_center = np.deg2rad(row['xbeta_center'])*u.rad
+
+            xdisp_groove_length = u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq']
+            echelle_groove_length = u.Unit(trace_params.meta["disp_freq_unit"]) / row['disp_freq']
+            pix_per_res_elem = row['fwhm']
+
+            spectrograph = echelle.spectrograph_factory(min_wave, max_wave, focal_len,
+                                                        design_res, echelle_angle, min_order, max_order,
+                                                        echelle_groove_length, pix_per_res_elem, disp_npix, xdisp_npix,
+                                                        pix_size, xdisp_groove_length=xdisp_groove_length,
+                                                        xdisp_beta_center=xdisp_beta_center)
+
+            spectrographs[prefix] = spectrograph
+
+        def efficiency_curve(trace_id, wavelength):
+            """Trace ID MUST be in the form prefix_{order}"""
+            prefix, _, order = trace_id.partition('_')
+            order = int(order)
+            spec = spectrographs[prefix]
+            blaze = spec.grating.blaze(spec.grating.beta(wavelength, order), order)
+            xdisp = spec.xdisp_efficiency(wavelength)
+            return blaze*xdisp
+
+        return efficiency_curve
+
+    def apply_to(self, obj, **kwargs):
+        """Interface between FieldOfView and SpectralEfficiency."""
+        trace_id = obj.trace_id
+
+        swcs = WCS(obj.hdu.header).spectral
+        with u.set_enabled_equivalencies(u.spectral()):
+            wave = swcs.pixel_to_world(np.arange(swcs.pixel_shape[0])) << u.um
+
+        efficiency = self.efficiency_generator(trace_id, wave)
+        params = {"description": trace_id}
+        params.update(self.meta)
+        effic_curve = TERCurve(array_dict={"wavelength": wave, "transmission": efficiency}, **params)
+        self.efficiencies[trace_id] = effic_curve
+
+        obj.hdu = apply_throughput_to_cube(obj.hdu, effic_curve.throughput, wave)
+        return obj
+
+    def plot(self):
+        """Plot the grating efficiencies."""
+        fig, axes = figure_factory()
+        for name, effic in self.efficiencies.items():
+            wave = effic.throughput.waveset
+            axes.plot(wave.to(u.um), effic.throughput(wave), label=name)
+
+        axes.set_xlabel("Wavelength [um]")
+        axes.set_ylabel("Grating efficiency")
+        axes.set_title(f"Grating efficiencies {self.display_name}")
+        axes.legend()
+
+        return fig
\ No newline at end of file
diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py
index 816cb260..3c1db629 100644
--- a/scopesim/effects/spectral_trace_list.py
+++ b/scopesim/effects/spectral_trace_list.py
@@ -5,7 +5,6 @@
 The Effect is called `SpectralTraceList`, it applies a list of
 `spectral_trace_list_utils.SpectralTrace` objects to a `FieldOfView`.
 """
-
 from itertools import cycle
 from typing import ClassVar
 
@@ -15,6 +14,7 @@
 from astropy.io import fits
 from astropy.table import Table
 import astropy.units as u
+import os
 
 from .effects import Effect
 from .ter_curves import FilterCurve
@@ -22,7 +22,7 @@
 from ..optics.image_plane_utils import header_from_list_of_xy
 from ..optics.fov import FieldOfView
 from ..optics.fov_volume_list import FovVolumeList
-from ..utils import from_currsys, check_keys, figure_factory, get_logger
+from ..utils import from_currsys, check_keys, figure_factory, get_logger, from_rc_config
 from .data_container import DataContainer
 from ..optics import echelle
 
@@ -108,8 +108,8 @@ class SpectralTraceList(Effect):
     report_plot_include: ClassVar[bool] = True
     report_table_include: ClassVar[bool] = False
 
-    def __init__(self, **kwargs):
-        super().__init__(**kwargs)
+    def __init__(self, cmds=None, **kwargs):
+        super().__init__(cmds=cmds, **kwargs)
 
         if "hdulist" in kwargs and isinstance(kwargs["hdulist"], fits.HDUList):
             self._file = kwargs["hdulist"]
@@ -137,7 +137,6 @@ def __init__(self, **kwargs):
 
         if self._file is not None:
             self.make_spectral_traces()
-
             self.update_meta()
 
     def make_spectral_traces(self):
@@ -524,7 +523,7 @@ class EchelleSpectralTraceList(SpectralTraceList):
     instead of loading them from FITS file. The arguments required to define the echelle traces are supplied through
     a txt file containing a table of parameters using the filename kwarg.
 
-    Below is an example of how to define the echelle trace parameters (see irdb/ZShooter/traces/echelle_trace_parameters.txt):
+    Below is an example of how to define the echelle trace parameters (see irdb/ZShooter_v1/traces/echelle_trace_parameters.txt):
     ----------------------------------------------------------------
     # min_wave_unit : nm
     # max_wave_unit : nm
@@ -539,10 +538,10 @@ class EchelleSpectralTraceList(SpectralTraceList):
     # xdisp_freq_unit : mm
     # slitwidth_unit : arcsec
 
-    prefix    aperture_id    image_plane_id    m0    n    min_wave    max_wave    echelle_blaze    focal_length    fwhm    detector_pad    pixel_size    n_disp    n_xdisp     disp_freq    xdisp_freq    slitwidth    dispdir    plate_scale
-    nIR        0              0                40    24    970         2500        64.2             225             4.7     10              0.015         4096      4096        45           175           10           x          0.159574468085
-    gri        1              1                36    18    490         1020        64.2             225             4.7     10              0.015         4096      4096        100          500           10           x          0.159574468085
-    ub         2              2                29    11    315         515         64.2             225             4.7     10              0.015         4096      4096        200          1000          10           x          0.159574468085
+    prefix    aperture_id    image_plane_id    m0    n    min_wave    max_wave   design_res    echelle_blaze    focal_length    fwhm    detector_pad    pixel_size    n_disp    n_xdisp     disp_freq    xdisp_freq    slitwidth    dispdir
+    ub         0              2                29    11    315         515          20000        64.2             225             4.7     10              0.015         4096      4096        200          1000          10           x
+    gri        1              1                36    18    490         1020         20000        64.2             225             4.7     10              0.015         4096      4096        100          500           10           x
+    nIR        2              0                40    24    970         2500         20000        64.2             225             4.7     10              0.015         4096      4096        45           175           10           x
     ----------------------------------------------------------------
 
     The calculated traces are stored in the same HDUList format as required by SpectralTraceList,
@@ -552,13 +551,19 @@ class EchelleSpectralTraceList(SpectralTraceList):
     required_keys = {"filename"}
     z_order = (71, 271, 671)
 
-    def __init__(self, **kwargs):
+    def __init__(self, cmds=None, **kwargs):
         check_keys(kwargs, self.required_keys, action="error")
+        self.cmds = cmds
 
-        trace_params = DataContainer(filename=kwargs['filename'])
+        trace_param_filename = kwargs.pop("filename")
+        trace_params = DataContainer(filename=trace_param_filename)
         hdulist = self._generate_trace_hdulist(trace_params)
+        hdulist.writeto(f"{from_rc_config('!SIM.file.local_packages_path')}/"
+                        f"{from_currsys('!OBS.instrument', self.cmds)}/"
+                        f"{os.path.dirname(trace_param_filename)}/"
+                        f"analytical_echelle_traces.fits", overwrite=True)
         kwargs["hdulist"] = hdulist
-        super().__init__(**kwargs)
+        super().__init__(cmds=cmds, **kwargs)
 
     def _generate_trace_hdulist(self, trace_params):
         hdul = fits.HDUList()
@@ -588,37 +593,29 @@ def _generate_trace_hdulist(self, trace_params):
             max_order = row['m0']
             min_wave = row['min_wave'] * u.Unit(trace_params.meta["min_wave_unit"])
             max_wave = row['max_wave'] * u.Unit(trace_params.meta["max_wave_unit"])
+            design_res = row['design_res']
             focal_len = row['focal_length'] * u.Unit(trace_params.meta["focal_length_unit"])
-            xdisp_npix = row['n_xdisp']
+            disp_npix = row['n_disp'] - 2 * row['detector_pad']
+            xdisp_npix = row['n_xdisp'] - 2 * row['detector_pad']
             pix_size = row['pixel_size'] * u.Unit(trace_params.meta["pixel_size_unit"])
-            x_disp_len = (xdisp_npix - 2 * row['detector_pad']) * pix_size
-            echelle_angle = np.deg2rad(row['echelle_blaze'])
-            alpha = np.deg2rad(row['alpha'])
-            beta_center = np.deg2rad(row['beta_center'])
-            # cross_disperser = echelle.GratingSetup(
-            #     groove_length=u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq'],
-            #     guess_littrow=(min_wave, max_wave,
-            #                    x_disp_len, focal_len))
-            cross_disperser = echelle.GratingSetup(alpha=alpha, beta_center=beta_center,
-                                                   delta=beta_center,
-                                                   groove_length=u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq'])
-
-            ss = echelle.SpectrographSetup((min_order, max_order),
-                                           max_wave,
-                                           row['fwhm'] * u.Unit(trace_params.meta["fwhm_unit"]),
-                                           focal_len,
-                                           echelle.GratingSetup(alpha=echelle_angle, beta_center=echelle_angle,
-                                                                delta=echelle_angle,
-                                                                groove_length=u.Unit(trace_params.meta["disp_freq_unit"]) / row['disp_freq']),
-                                           echelle.Detector(row['n_disp'], xdisp_npix, pix_size),
-                                           cross_disperser=cross_disperser
-                                           )
+            echelle_angle = np.deg2rad(row['echelle_blaze'])*u.rad
+            xdisp_beta_center = np.deg2rad(row['xbeta_center'])*u.rad
+
+            xdisp_groove_length = u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq']
+            echelle_groove_length = u.Unit(trace_params.meta["disp_freq_unit"]) / row['disp_freq']
+            pix_per_res_elem = row['fwhm']
+
+            ss = echelle.spectrograph_factory(min_wave, max_wave, focal_len,
+                                              design_res, echelle_angle, min_order, max_order,
+                                              echelle_groove_length, pix_per_res_elem, disp_npix, xdisp_npix,
+                                              pix_size, xdisp_groove_length=xdisp_groove_length,
+                                              xdisp_beta_center=xdisp_beta_center)
 
             fsr_edges = ss.edge_wave(fsr=True)
 
-            slit_edge = (row['slitwidth'] / 2) * u.Unit(trace_params.meta["slitwidth_unit"])
+            slit_edge = (row['slitlength'] / 2) * u.Unit(trace_params.meta["slitlength_unit"])
             slit_pos = np.linspace(-slit_edge, slit_edge, num=3)
-            slit_offset_pix = slit_pos / (row['plate_scale'] * u.arcsec)
+            slit_offset_pix = slit_pos / (from_currsys('!INST.pixel_scale', self.cmds) * u.arcsec)
 
             xvals, yvals = [], []
             for i, order in enumerate(ss.orders):
diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py
index 49e394de..200464c1 100644
--- a/scopesim/effects/ter_curves.py
+++ b/scopesim/effects/ter_curves.py
@@ -16,6 +16,7 @@
 
 import skycalc_ipy
 
+from . import DataContainer
 from .effects import Effect
 from .ter_curves_utils import (add_edge_zeros, combine_two_spectra,
                                apply_throughput_to_cube, download_svo_filter,
@@ -1247,80 +1248,15 @@ def _guard_plot_axes(which, axes):
                              "contains only one element"))
 
 
-class Dichroic(TERCurve):
-    """
-    An effect that holds and applies a single dichroic TERCurve.
-    The effect splits each input FoV into two FoVs, one for the transmission arm
-    and one for the reflection arm. "arm_key" kwarg indicates which id key to assign value to (for e.g. aperture_id).
-    "arm_action" kwarg indicates which arm_key value to assign for transmission and reflection arms respectively.
-
-    Example config definition:
-    ```
-    name: dichroic1
-    class: Dichroic
-    kwargs:
-        filename: "TER_dichroic1.dat"
-        arm_key: "aperture_id"
-        arm_actions:
-          0: "transmission"
-          1: "reflection"
-
-    ```
-    """
-
-    z_order: ClassVar[tuple[int, ...]] = (116, 216, 616)
-    required_keys = {"arm_key", "arm_actions"}
-
-    def __init__(self, cmds=None, **kwargs):
-        check_keys(kwargs, self.required_keys, action="error")
-        self.arm_key = kwargs.pop("arm_key")
-        self.arm_actions = kwargs.pop("arm_actions")
-
-        super().__init__(cmds=cmds, **kwargs)
-
-    def apply_to(self, obj, **kwargs):
-        """
-        For FoV volume list, create two FoVs per input FoV for transmission and reflection arms, assign
-        corresponding arm_key values.
-        For single FoV, apply the throughput based on its arm_key value to determine transmission or reflection.
-        """
-        if isinstance(obj, FovVolumeList):
-            logger.debug("Executing %s, FoV setup", self.meta['name'])
-            new_vols = []
-            for idval, action in zip(self.arm_actions):
-                vols = deepcopy(obj)
-                params = {"action": action}
-                self.surface.meta.update(params)
-                vols = super().apply_to(vols, **kwargs)
-                for vol in vols.volumes:
-                    vol["meta"][self.arm_key] = idval
-                new_vols.extend(vols.volumes)
-            obj.volumes = new_vols
-
-        if isinstance(obj, FieldOfView):
-            logger.debug("Executing %s, observe", self.meta['name'])
-            idval = obj.meta.get(self.arm_key, None)
-            if idval is None:
-                raise ValueError(f"This FoV does not have the correct {self.arm_key} to determine "
-                                 f"transmission or reflection arm.")
-            if idval not in self.arm_actions.keys():
-                raise ValueError(f"{self.arm_key} value {idval} not found in arm_action mapping.")
-            action = self.arm_actions[idval]
-            params = {"action": action}
-            self.surface.meta.update(params)
-            super().apply_to(obj, **kwargs)
-
-        return obj
-
-
 class DichroicTree(Effect):
     """
     An effect that holds and applies a tree of dichroics.
 
-    The tree is defined in a table in a txt file where each row is a tree branch and each column is a dichroic.
-    The column names should match the names of the TERCurve file names.
-    The cells contain either a X (for no dichroic at that position), a T (for a dichroic transmission branch)
-    or a R (for a dichroic reflection branch).
+    The tree is defined in a table in a .dat file where each row is a tree branch and columns are dichroics.
+    The first column is an id column to identify the branch. Its name should be the "id" that is set in the FoV meta during setup, e.g. image_plane_id or aperture_id.
+    The other columns are dichroic names and should match the dichroic_names property in the config.
+    The corresponding TER files for the dichroics should contain the dichroic_name in the filename according to the filename_format pattern defined in the config.
+    The cells contain either a X (for no dichroic at that position), a T (for a dichroic transmission branch) or a R (for a dichroic reflection branch).
 
     Example structure of the dichroic tree file with two dichroics:
     ```
@@ -1332,12 +1268,18 @@ class DichroicTree(Effect):
     ```
     The config definition would look like this:
     ```
-    name: dichroic_tree
-    class: DichroicTreeEffect
-    kwargs:
-        tree_filename: "dichroic_tree.dat"
-        dichroic_names: ["dichroic1", "dichroic2"]
-        filename_format: "TER_{}.dat"
+    properties:
+      dichroic_names:
+        - dichroic1
+        - dichroic2
+
+    effects:
+      - name: dichroic_tree
+        class: DichroicTree
+        kwargs:
+          tree_filename: "dichroic_tree.dat"
+          dichroic_names: "!INST.dichroic_names"
+          filename_format: "TER_{}.dat"
     ```
 
     This effect creates a volume of FoVs, one for each row in the dichroic tree table.
@@ -1346,10 +1288,8 @@ class DichroicTree(Effect):
     z_order: ClassVar[tuple[int, ...]] = (117, 217, 617)
 
     def __init__(self, cmds=None, **kwargs):
-        super().__init__(cmds=cmds, **kwargs)
         check_keys(kwargs, self.required_keys, action="error")
-
-        self.meta.update(kwargs)
+        super().__init__(cmds=cmds, **kwargs)
 
         # Load all the dichroic TER curves
         self.dichroics = {}
@@ -1363,7 +1303,8 @@ def __init__(self, cmds=None, **kwargs):
 
         # Load the dichroic tree structure
         tree_path = find_file(from_currsys(self.meta["tree_filename"], cmds=self.cmds))
-        self.table = ioascii.read(tree_path, format="basic")
+        tree = DataContainer(filename=tree_path)
+        self.table = tree.table
 
     def apply_to(self, obj, **kwargs):
         """
@@ -1373,16 +1314,17 @@ def apply_to(self, obj, **kwargs):
         The output should be a list of FoVs, one for each row in the table.
         """
         action_lookup = {"T": "transmission", "R": "reflection", "X": None}
+        id_colname = self.table.colnames[0]
+        dichroic_cols = self.table.colnames[1:]
 
         # 1. During setup of the FieldofViews
         if isinstance(obj, FovVolumeList):
             logger.debug("Executing %s, FoV setup", self.meta['name'])
             new_vols = []
             for row in self.table:
-                arm_id = row["id"]
+                arm_id = row[id_colname]
                 vols = deepcopy(obj)
 
-                dichroic_cols = [col for col in row.colnames if col != "id"]
                 for dichroic_name in dichroic_cols:
                     if dichroic_name not in self.dichroics.keys():
                         raise ValueError(f"Dichroic name {dichroic_name} not found in saved dichroics.")
@@ -1393,25 +1335,24 @@ def apply_to(self, obj, **kwargs):
                         dichroic.surface.meta.update(params)
                         vols = dichroic.apply_to(vols, **kwargs)
 
-                # Set the id of the resulting FoVs. Looks like id field in meta is not being used for anything else.
-                ## TODO: Update to aperture_id after SelectorWheel is in place??
+                # Set the arm_id value of the resulting FoVs
                 for vol in vols.volumes:
-                    vol["meta"]["id"] = arm_id
+                    vol["meta"][id_colname] = arm_id
                 new_vols.extend(vols.volumes)
+                logger.debug(f"Created {len(vols.volumes)} FoVs for {id_colname} {arm_id}")
 
             obj.volumes = new_vols
 
         # 2. During observe
         if isinstance(obj, FieldOfView):
             logger.debug("Executing %s, observe", self.meta['name'])
-            obj_id = obj.meta.get("id", None)
-            row = self.table[self.table["id"] == obj_id]
+            obj_id = obj.meta.get(id_colname, None)
+            row = self.table[self.table[id_colname] == obj_id]
 
             if len(row) != 1:
-                raise ValueError(f"This FoV does not have the correct id to match "
-                                 f"a row in the dichroic tree table: {obj_id}")
+                raise ValueError(f"This FoV's {id_colname}: {obj_id} does not match "
+                                 f"a row in the dichroic tree table.")
 
-            dichroic_cols = [col for col in row.colnames if col != "id"]
             for dichroic_name in dichroic_cols:
                 if dichroic_name not in self.dichroics.keys():
                     raise ValueError(f"Dichroic name {dichroic_name} not found in saved dichroics.")
diff --git a/scopesim/optics/echelle.py b/scopesim/optics/echelle.py
index bdb1100f..e23847a7 100644
--- a/scopesim/optics/echelle.py
+++ b/scopesim/optics/echelle.py
@@ -2,17 +2,109 @@
 
 import numpy as np
 import astropy.units as u
-import logging
+
+from ..utils import get_logger
+logger = get_logger(__name__)
+
+def spectrograph_factory(min_wave: float|u.Quantity, max_wave: float|u.Quantity, focal_len: float|u.Quantity,
+                         design_res: float, echelle_angle: float|u.Quantity, min_order: int, max_order: int,
+                         echelle_groove_length: float|u.Quantity,
+                         pix_per_res_elem: float, disp_npix: int, xdisp_npix: int, pix_size: float|u.Quantity,
+                         xdisp_groove_length: float|u.Quantity = 0.0, xdisp_beta_center: float|u.Quantity = 0.0):
+    """
+    
+    Parameters
+    ----------
+    min_wave: float|u.Quantity
+        Minimum wavelength. If float, assumed to be in nm.
+    max_wave: float|u.Quantity
+        Maximum wavelength. If float, assumed to be in nm.
+    focal_len: float|u.Quantity
+        Focal length. If float, assumed to be in mm.
+    design_res: float|int
+        Design resolution (dimensionless).
+    echelle_angle: float|u.Quantity
+        Echelle angle. If float, assumed to be in degrees.
+    min_order: int
+        Minimum order.
+    max_order: int
+        Maximum order.
+    echelle_groove_length: float|u.Quantity
+        Echelle groove length. If float, assumed to be in mm.
+    pix_per_res_elem: float|int
+        Pixels per resolution element.
+    disp_npix: int
+        Number of dispersion pixels.
+    xdisp_npix: int
+        Number of cross-dispersion pixels.
+    pix_size: float|u.Quantity
+        Pixel size. If float, assumed to be in um.
+    xdisp_groove_length: float|u.Quantity
+        Cross disperser groove length. If float, assumed to be in mm.
+    xdisp_beta_center: float|u.Quantity
+        Cross disperser beta center angle. If float, assumed to be in degrees.
+
+    Returns
+    -------
+
+    """
+    # Convert primitive types to quantities with appropriate units
+    if not isinstance(min_wave, u.Quantity):
+        min_wave = min_wave * u.nm
+    if not isinstance(max_wave, u.Quantity):
+        max_wave = max_wave * u.nm
+    if not isinstance(focal_len, u.Quantity):
+        focal_len = focal_len * u.mm
+    if not isinstance(echelle_angle, u.Quantity):
+        echelle_angle = echelle_angle * u.deg
+    if not isinstance(echelle_groove_length, u.Quantity):
+        echelle_groove_length = echelle_groove_length * u.mm
+    if not isinstance(pix_size, u.Quantity):
+        pix_size = pix_size * u.mm
+    if not isinstance(xdisp_groove_length, u.Quantity):
+        xdisp_groove_length = xdisp_groove_length * u.mm
+    if not isinstance(xdisp_beta_center, u.Quantity):
+        xdisp_beta_center = xdisp_beta_center * u.deg
+
+    x_disp_len = (xdisp_npix*pix_size).to(u.mm)
+
+    if not np.isfinite(xdisp_groove_length):
+        nu, xsdisp_angle = GratingSetup.estimate_xdisp_freq_and_angle(min_wave, max_wave, x_disp_len, focal_len)
+        cross_disperser = GratingSetup(alpha=xsdisp_angle, beta_center=xsdisp_angle, groove_length=1 / nu,
+                                       grating_type='vph', vph_center_wave=(min_wave + max_wave) / 2)
+    else:
+        if xdisp_beta_center != 0:  # TODO probably hsould be a nan, but use of 0 angle incidence isn't likely
+            cross_disperser = GratingSetup(alpha=xdisp_beta_center, beta_center=xdisp_beta_center,
+                                           groove_length=xdisp_groove_length, grating_type='vph',
+                                           vph_center_wave=xdisp_groove_length.to(u.nm) * 2 * np.sin(xdisp_beta_center))
+        else:
+            cross_disperser = GratingSetup(groove_length=xdisp_groove_length,
+                                           guess_littrow=(min_wave, max_wave, x_disp_len, focal_len),
+                                           grating_type='vph', vph_center_wave=(min_wave + max_wave) / 2)
+
+    echelle_grating = GratingSetup(alpha=echelle_angle, beta_center=echelle_angle, delta=echelle_angle,
+                                   groove_length=echelle_groove_length)
+
+    return SpectrographSetup(order_range=(min_order, max_order),
+                                   design_res=design_res,
+                                   pixels_per_res_elem=pix_per_res_elem,
+                                   focal_length=focal_len,
+                                   grating=echelle_grating,
+                                   detector=Detector(disp_npix, xdisp_npix, pix_size),
+                                   cross_disperser=cross_disperser)
+
 
 class GratingSetup:
     def __init__(
             self, *,
-            groove_length:u.Quantity=None,
-            alpha:float=None,
-            delta:float = None,
-            beta_center:float = None,
-            empiric_blaze_factor: float = 1.0,
-            guess_littrow: tuple = None,
+            alpha: u.Quantity = None,
+            delta: u.Quantity = None,
+            beta_center: u.Quantity = None,
+            groove_length: u.Quantity = None,
+            empiric_efficiency_factor: float = 1.0,
+            guess_littrow: tuple[u.Quantity,u.Quantity,u.Quantity,u.Quantity] = None,
+            grating_type: str = 'echelle',
+            vph_center_wave: u.Quantity = None
     ):
         """
         Simulation of an echelle/echellete grating for spectrometer.
@@ -21,34 +113,142 @@ def __init__(
         :param float delta: blaze angle in radians
         :param float beta_center: reflectance angle in radians
         :param u.Quantity groove_length: d, length/groove in units of wavelength (same as schroeder sigma)
-        :param float empiric_blaze_factor: empirically determined blaze factor to correct for grating efficiency relative to ideal peak
-        :param bool guess_littrow: if not None, passed onto estimate_xdisp_angle
+        :param float empiric_efficiency_factor: empirically determined blaze factor to correct for grating efficiency relative to ideal peak
+        :param bool guess_littrow: if not None, passed onto estimate_xdisp_angle_with_groove_len
         """
+        self.grating_type = grating_type.lower()
+        assert self.grating_type in ('echelle', 'vph'), f'grating_type must be echelle or vph, not {self.grating_type}'
+
         if guess_littrow:
-            alpha, beta_center = self.estimate_xdisp_angle(guess_littrow[0], guess_littrow[1], guess_littrow[2],
+            alpha, beta_center = self.estimate_xdisp_angle_with_groove_len(guess_littrow[0], guess_littrow[1], guess_littrow[2],
                                               guess_littrow[3], groove_length)
             self.alpha = alpha
-            self.delta = beta_center
+            self.delta = None if self.grating_type == 'vph' else beta_center
             self.beta_center = beta_center
         else:
-            self.alpha = alpha*u.rad
-            self.delta = delta*u.rad
-            self.beta_center = beta_center*u.rad
-        self.d = groove_length
-        self.empiric_blaze_factor = empiric_blaze_factor
+            self.alpha = alpha.to(u.rad)
+            self.delta = delta.to(u.rad) if delta is not None else None
+            self.beta_center = beta_center.to(u.rad)
+        self.d = groove_length.to(u.mm)
+        self.empiric_efficiency_factor = empiric_efficiency_factor
+
+        # Approximate reality vph_constant = np.pi * deltan * d / bragg_angle
+        # Empircally we would use with peak efficiency near center of range
+        self.vph_constant = np.pi/2*vph_center_wave.to(u.nm) if self.grating_type == 'vph' else None
 
     @staticmethod
-    def estimate_xdisp_angle(l_min, l_max, detector_length, focal_length, d):
+    def estimate_xdisp_angle_with_groove_len(
+            l_min: float | u.Quantity,
+            l_max: float | u.Quantity,
+            detector_length: float | u.Quantity,
+            focal_length: float | u.Quantity,
+            d: u.Quantity) -> tuple[u.Quantity, u.Quantity]:
+        """
+        Estimate cross-disperser angle and beta center for a given groove length.
+        
+        This method calculates the optimal incident angle (alpha) and central reflection 
+        angle (beta) for a cross-disperser grating to match the wavelength range to the 
+        detector length at the focal plane.
+        
+        Parameters
+        ----------
+        l_min : float | u.Quantity
+            Minimum wavelength. If float, assumed to be in nm.
+        l_max : float | u.Quantity
+            Maximum wavelength. If float, assumed to be in nm.
+        detector_length : float | u.Quantity
+            Length of the detector. If float, assumed to be in mm.
+        focal_length : float | u.Quantity
+            Focal length of the system. If float, assumed to be in mm.
+        d : u.Quantity
+            Groove length (spacing) of the grating.
+            
+        Returns
+        -------
+        tuple[u.Quantity, u.Quantity]
+            A tuple containing:
+            - alpha: incident angle in radians
+            - beta_center: central reflection angle in radians
+        """
+        # Convert floats to quantities with appropriate units
+        if not isinstance(l_min, u.Quantity):
+            l_min = l_min * u.nm
+        if not isinstance(l_max, u.Quantity):
+            l_max = l_max * u.nm
+        if not isinstance(detector_length, u.Quantity):
+            detector_length = detector_length * u.mm
+        if not isinstance(focal_length, u.Quantity):
+            focal_length = focal_length * u.mm
+        if not isinstance(d, u.Quantity):
+            d = d * u.mm
+
         # assume we want a dispersion at the midpoint that matches the range to the detector length
-        t = detector_length/focal_length
-        k = (l_max-l_min)/d
-        x = k*(1+np.sqrt(1+t**2))/2/t
-        u=np.arccos(np.sqrt(x))
-        v=np.arcsin(k/2/np.sqrt(x))
+        t = detector_length / focal_length
+        k = (l_max - l_min) / d
+        x = k * (1 + np.sqrt(1 + t ** 2)) / 2 / t
+        u = np.arccos(np.sqrt(x))
+        v = np.arcsin(k / 2 / np.sqrt(x))
         beta_max = u+v
         beta_min = u-v
-        alpha = np.arcsin(l_max/d -np.sin(beta_max))
-        return alpha, (beta_min+beta_max)/2
+        alpha = np.arcsin(l_max / d - np.sin(beta_max))
+        return alpha.to(u.rad), ((beta_min + beta_max) / 2).to(u.rad)
+
+    @staticmethod
+    def estimate_xdisp_freq_and_angle(l_min: float | u.Quantity,
+                                     l_max: float | u.Quantity,
+                                     detector_length: float | u.Quantity,
+                                     focal_length: float | u.Quantity) -> tuple[u.Quantity, u.Quantity]:
+        """
+        Estimate cross-disperser littrow angle and groove frequency.
+
+        This method calculates the analytical solution to the system of equations
+
+        nu l_max == Sin(beta)+Sin(beta+deltabeta) nu l_min == Sin(beta)+Sin(beta-deltabeta)
+        cos[deltabeta] = f/sqrt(f^2 + (d/2)^)
+        sin[deltabeta] = (d/2)/sqrt(f^2 + (d/2)^)
+
+        to match the initial and final wavelengths to the extent of the detector length
+
+
+        Parameters
+        ----------
+        l_min : float | u.Quantity
+            Minimum wavelength. If float, assumed to be in nm.
+        l_max : float | u.Quantity
+            Maximum wavelength. If float, assumed to be in nm.
+        detector_length : float | u.Quantity
+            Length of the detector. If float, assumed to be in mm.
+        focal_length : float | u.Quantity
+            Focal length of the system. If float, assumed to be in mm.
+
+        Returns
+        -------
+        tuple[u.Quantity, u.Quantity]
+            A tuple containing:
+            - ruling frequency
+            - littrow angle
+        """
+        if not isinstance(l_min, u.Quantity):
+            l_min = l_min * u.nm
+        if not isinstance(l_max, u.Quantity):
+            l_max = l_max * u.nm
+        if not isinstance(detector_length, u.Quantity):
+            detector_length = detector_length * u.mm
+        if not isinstance(focal_length, u.Quantity):
+            focal_length = focal_length * u.mm
+
+        f = focal_length
+        p = detector_length
+        li, lf = l_min, l_max
+        a = 4*f**2 + p**2
+        rta = np.sqrt(a)
+        b = lf - li
+        bsq = b**2
+        c = lf**2 + li**2
+        d = 4 * f**2 * bsq + 2 * f * rta * bsq + p**2 * c
+        beta = np.arccos(((2 * f + rta) * b) / np.sqrt(2 * d))
+        nu = (np.sqrt(2) * p * (p**2 + 4 * f * (2 * f + rta))) / ((2 * f + rta) * np.sqrt(a * d))
+        return nu.to('mm-1'), beta.to('rad')
 
     def __str__(self) -> str:
         return (f"alpha={np.rad2deg(self.alpha):.2f}\n"
@@ -64,6 +264,8 @@ def blaze(self, beta, m):
         :param m: order
         :return: grating throughput out of 1
         """
+        if self.grating_type == 'vph':
+            raise ValueError('Cannot calculate blaze for vph gratings.')
         k = np.cos(beta) * np.cos(self.alpha - self.delta) / (np.cos(self.alpha) * np.cos(beta - self.delta))
         k[k > 1] = 1
         # q1 = np.cos(alpha) / np.cos(alpha - delta)
@@ -73,7 +275,17 @@ def blaze(self, beta, m):
         rho = np.cos(self.delta) if self.alpha < self.delta else np.cos(self.alpha) / np.cos(self.alpha - self.delta)
         # 2 different rho depending on whether alpha or delta is larger
         ret = k * np.sinc((m * rho * q4).value) ** 2  # omit np.pi as np.sinc includes it
-        return self.empiric_blaze_factor * ret
+        return self.empiric_efficiency_factor * ret
+
+    def vph_efficiency(self, wavelength):
+        """
+        :param wavelength: wavelength(s) as u.Quantity units of length
+        :return: efficiency of the grating
+        """
+        if self.grating_type != 'vph':
+            raise ValueError('Cannot calculate vph_effiency for echelle.')
+
+        return self.empiric_efficiency_factor * np.sin((self.vph_constant/wavelength).decompose().value)**2
 
     def beta(self, wave, m):
         """
@@ -136,7 +348,7 @@ class SpectrographSetup:
     def __init__(
             self,
             order_range: tuple,
-            final_wave: u.Quantity,
+            design_res: float,
             pixels_per_res_elem: float,
             focal_length: u.Quantity,
             grating: GratingSetup,
@@ -150,7 +362,7 @@ def __init__(
         (and that the slit image is gaussian)
 
         :param tuple order_range: order range of the spectrograph
-        :param u.Quantity final_wave: longest wavelength at the edge of detector
+        # :param u.Quantity final_wave: longest wavelength at the edge of detector
         :param float pixels_per_res_elem: number of pixels per resolution element of spectrometer
         :param u.Quantity focal_length: the focal length of the detector
         :param GratingSetup grating: configured grating
@@ -162,10 +374,11 @@ def __init__(
             order_range = order_range[::-1]
         self.m0 = order_range[0]
         self.m_max = order_range[1]
-        self.l0 = final_wave
+        self.l0 = grating.wave(grating.beta_center, order_range[0])*(1/order_range[0]/2+1)
+        self.design_res = design_res
         self.grating = grating
         self.detector = detector
-        self.focal_length = focal_length
+        self.focal_length = focal_length.to(u.mm)
         self.pixel_scale = np.arctan(self.detector.pixel_size / self.focal_length)
         self.beta_central_pixel = self.grating.beta_center
         self.nord = int(self.m_max - self.m0 + 1)
@@ -174,7 +387,7 @@ def __init__(
         self._orders = None
 
         self.nondimensional_lsf_width = 1 / self.design_res
-        logging.info(f'\nThe spectrograph has been setup with the following properties:'
+        logger.debug(f'\nThe spectrograph has been setup with the following properties:'
                      f'\n\tl0: {self.l0}'
                      # f'\n\tR0: {self.detector.design_R0}'
                      f'\n\tOrders: {self.orders}'
@@ -282,6 +495,13 @@ def blaze(self, wave):
         """
         return self.grating.blaze(self.grating.beta(wave, self.orders[:, None]), self.orders[:, None])
 
+    def xdisp_efficiency(self, wave):
+        """
+        :param wave: wavelength
+        :return: efficiency of the cross disperser
+        """
+        return self.cross_disperser.vph_efficiency(wave)
+
     def mean_blaze_eff_est(self, n=10):
         """
         :param n:
@@ -445,7 +665,7 @@ def angular_dispersion(self):
         return self.grating.angular_dispersion(self.orders[:, None], beta)
 
     @property
-    def design_res(self):
+    def implied_design_res(self):
         """
         :return: design resolution for spectrometer
         Assumes that m0 (the longest wavelength order) FSR fills detector with some sampling
@@ -489,7 +709,7 @@ def plot_echellogram(self, center_orders=True, title='', blaze=False, cross_disp
             oset = waves[1] if center_orders else 0
             plt.plot(waves - oset, [i] * 3, '*', color=f'C{ii}')
             plt.plot(fsr_edges[ii] - oset, [i] * 2, '.', color=f'C{ii}',
-                     label=f'$\lambda=${waves[1]:.0f}')
+                     label=f'$\\lambda=${waves[1]:.0f}')
         plt.legend()
         plt.xlabel('Center relative wavelength (nm)')
         plt.ylabel('Order')