diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py
index d523aee7..6b7657df 100644
--- a/scopesim/effects/spectral_trace_list_utils.py
+++ b/scopesim/effects/spectral_trace_list_utils.py
@@ -772,14 +772,42 @@ def __init__(self, fov, dlam_per_pix):
         self.xi = self.wcs.all_pix2world(self.lam[0], np.arange(n_xi), 0)[1]
         self.npix_xi = n_xi
         self.npix_lam = n_lam
-        # ..todo: cubic spline introduces negative values, linear does not.
-        #  Alternative might be to cubic-spline interpolate on sqrt(image),
-        #  with subsequent squaring of the result. This would require
-        #  wrapping RectBivariateSpline in a new (sub)class.
-        spline_order = (1, 1)
-        self.interp = RectBivariateSpline(self.xi, self.lam, self.image,
-                                          kx=spline_order[0],
-                                          ky=spline_order[1])
+
+    def interp(self, xi, lam, grid=False):
+        """
+        Bilinear interpolation of the xi-lambda image at (xi, lam).
+
+        Direct equivalent of evaluating a ``RectBivariateSpline`` with
+        ``kx=ky=1`` on the regular (xi, lam) grid, including the clamping
+        of out-of-range coordinates to the boundary values -- but avoiding
+        FITPACK's point-evaluation path, which is orders of magnitude
+        slower for large coordinate arrays.
+
+        With ``grid=False``, `xi` and `lam` are coordinate pairs of equal
+        shape. With ``grid=True``, the image is evaluated on the tensor
+        product of the two vectors.
+        """
+        xi = np.asarray(xi, dtype=float)
+        lam = np.asarray(lam, dtype=float)
+
+        i_xi = np.clip(np.searchsorted(self.xi, xi) - 1,
+                       0, self.npix_xi - 2)
+        j_lam = np.clip(np.searchsorted(self.lam, lam) - 1,
+                        0, self.npix_lam - 2)
+        w_xi = np.clip((xi - self.xi[i_xi])
+                       / (self.xi[i_xi + 1] - self.xi[i_xi]), 0., 1.)
+        w_lam = np.clip((lam - self.lam[j_lam])
+                        / (self.lam[j_lam + 1] - self.lam[j_lam]), 0., 1.)
+
+        if grid:
+            i_xi, w_xi = i_xi[:, None], w_xi[:, None]
+            j_lam, w_lam = j_lam[None, :], w_lam[None, :]
+
+        img = self.image
+        return ((1 - w_xi) * (1 - w_lam) * img[i_xi, j_lam]
+                + (1 - w_xi) * w_lam * img[i_xi, j_lam + 1]
+                + w_xi * (1 - w_lam) * img[i_xi + 1, j_lam]
+                + w_xi * w_lam * img[i_xi + 1, j_lam + 1])
 
 
 class Transform2D():
diff --git a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py
index ca2f19d2..cd760ca0 100644
--- a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py
+++ b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py
@@ -161,7 +161,7 @@ def test_grid_false_shape_is_preserved(self, tf2d):
 
 class MockCubeFov:
     """Minimal stand-in for a FieldOfView carrying a spectral cube."""
-    def __init__(self, n_lam=20, n_eta=3, n_xi=11):
+    def __init__(self, n_lam=20, n_eta=3, n_xi=11, data=None):
         hdr = fits.Header()
         hdr["NAXIS"] = 3
         hdr["NAXIS1"], hdr["NAXIS2"], hdr["NAXIS3"] = n_xi, n_eta, n_lam
@@ -170,8 +170,9 @@ def __init__(self, n_lam=20, n_eta=3, n_xi=11):
         hdr["CRVAL1"], hdr["CRVAL2"], hdr["CRVAL3"] = 0., 0., 2.0
         hdr["CDELT1"], hdr["CDELT2"], hdr["CDELT3"] = 0.1, 0.1, 0.001
         hdr["CUNIT1"], hdr["CUNIT2"], hdr["CUNIT3"] = "arcsec", "arcsec", "um"
-        self.cube = fits.ImageHDU(
-            data=np.ones((n_lam, n_eta, n_xi)), header=hdr)
+        if data is None:
+            data = np.ones((n_lam, n_eta, n_xi))
+        self.cube = fits.ImageHDU(data=data, header=hdr)
         self.meta = {"xi_min": -0.5 * u.arcsec, "xi_max": 0.5 * u.arcsec}
 
 
@@ -182,6 +183,34 @@ def test_primary_wcs_keeps_arcsec_cunit(self):
         assert list(xilam.wcs.wcs.cunit) == [u.um, u.arcsec]
         assert list(xilam.wcsa.wcs.cunit) == [u.um, u.dimensionless_unscaled]
 
+    def test_interp_matches_rect_bivariate_spline(self):
+        """XiLamImage.interp must reproduce RectBivariateSpline(kx=ky=1)
+        point evaluation, including clamping outside the grid."""
+        from scipy.interpolate import RectBivariateSpline
+
+        rng = np.random.default_rng(7)
+        n_lam, n_eta, n_xi = 30, 3, 12
+        fov = MockCubeFov(n_lam, n_eta, n_xi,
+                          data=rng.random((n_lam, n_eta, n_xi)))
+        xilam = XiLamImage(fov, dlam_per_pix=0.001)
+
+        spline = RectBivariateSpline(xilam.xi, xilam.lam, xilam.image,
+                                     kx=1, ky=1)
+
+        # scattered points, deliberately extending beyond the grid
+        xi_pts = rng.uniform(xilam.xi[0] - 0.2, xilam.xi[-1] + 0.2, 500)
+        lam_pts = rng.uniform(xilam.lam[0] - 0.01, xilam.lam[-1] + 0.01, 500)
+        np.testing.assert_allclose(
+            xilam.interp(xi_pts, lam_pts, grid=False),
+            spline(xi_pts, lam_pts, grid=False), rtol=1e-10)
+
+        # grid evaluation
+        xi_vec = np.sort(rng.uniform(xilam.xi[0], xilam.xi[-1], 7))
+        lam_vec = np.sort(rng.uniform(xilam.lam[0], xilam.lam[-1], 9))
+        np.testing.assert_allclose(
+            xilam.interp(xi_vec, lam_vec, grid=True),
+            spline(xi_vec, lam_vec, grid=True), rtol=1e-10)
+
 
 class TestImageInterpolations:
     """Tests for function make_image_interpolations"""