diff --git a/scopesim/optics/fov.py b/scopesim/optics/fov.py
index 503214f5..b92dd1ec 100644
--- a/scopesim/optics/fov.py
+++ b/scopesim/optics/fov.py
@@ -1165,8 +1165,10 @@ def extract_range_from_spectrum(spectrum, waverange):
     assert isinstance(spectrum, SourceSpectrum), (
         f"spectrum must be of type synphot.SourceSpectrum: {type(spectrum)}")
 
-    wave_min, wave_max = quantify(waverange, u.um).to(u.AA).value
-    spec_waveset = spectrum.waveset.to(u.AA).value
+    # Round to 8 decimal places, because the conversion to u.AA will cause
+    # irrelevant floating point differences.
+    wave_min, wave_max = quantify(waverange, u.um).to(u.AA).value.round(8)
+    spec_waveset = spectrum.waveset.to(u.AA).value.round(8)
     mask = (spec_waveset > wave_min) * (spec_waveset < wave_max)
 
     # FIXME: Why did I comment this out in 2023? Seems useful to have...
diff --git a/scopesim/optics/optical_train.py b/scopesim/optics/optical_train.py
index 4b3448cb..075ae6af 100644
--- a/scopesim/optics/optical_train.py
+++ b/scopesim/optics/optical_train.py
@@ -340,7 +340,12 @@ def prepare_source(self, source):
                 wave_max = max(fov.meta["wave_max"] for fov in self.fov_manager.fovs)
                 wave_unit = u.Unit(from_currsys("!SIM.spectral.wave_unit", self.cmds))
                 dwave = from_currsys("!SIM.spectral.spectral_bin_width", self.cmds)  # Not a quantity
-                fov_waveset = np.arange(wave_min.value, wave_max.value, dwave) * wave_unit
+                # Include the last point by adding dwave, because arange excludes the last point.
+                fov_waveset = np.arange(wave_min.value, wave_max.value + 0.9 * dwave, dwave)
+                # Ensure the start and end values are exact.
+                fov_waveset[0] = wave_min.value
+                fov_waveset[-1] = wave_max.value
+                fov_waveset *= wave_unit
                 fov_waveset = fov_waveset.to(u.um)
 
                 field.spectra[ispec] = SourceSpectrum(Empirical1D,
@@ -382,11 +387,17 @@ def prepare_source(self, source):
             cube.header["CUNIT2"] = "deg"
 
             # Put on fov wavegrid
+            # TODO: This is the same code as earlier in the function; refactor.
             wave_min = min(fov.meta["wave_min"] for fov in self.fov_manager.fovs)
             wave_max = max(fov.meta["wave_max"] for fov in self.fov_manager.fovs)
             wave_unit = u.Unit(from_currsys("!SIM.spectral.wave_unit", self.cmds))
             dwave = from_currsys("!SIM.spectral.spectral_bin_width", self.cmds)  # Not a quantity
-            fov_waveset = np.arange(wave_min.value, wave_max.value, dwave) * wave_unit
+            # Include the last point by adding dwave, because arange excludes the last point.
+            fov_waveset = np.arange(wave_min.value, wave_max.value + 0.9 * dwave, dwave)
+            # Ensure the start and end values are exact.
+            fov_waveset[0] = wave_min.value
+            fov_waveset[-1] = wave_max.value
+            fov_waveset *= wave_unit
             fov_waveset = fov_waveset.to(u.um)
 
             if (wave.to(u.um).min() > fov_waveset.max() or
diff --git a/scopesim/tests/tests_optics/test_fov_utls.py b/scopesim/tests/tests_optics/test_fov_utls.py
index ba08540d..f7252e31 100644
--- a/scopesim/tests/tests_optics/test_fov_utls.py
+++ b/scopesim/tests/tests_optics/test_fov_utls.py
@@ -112,7 +112,10 @@ def test_extracts_the_wave_range_needed(self):
         new_spec = extract_range_from_spectrum(spec, waverange)
 
         assert len(new_spec.waverange) == 2
-        assert new_spec.waverange[0] == 1.98 * u.um
+        # Cannot test for equality, because this does not hold in astropy.units:
+        # 1.98 * u.um == 19800.0 * u.AA
+        # assert new_spec.waverange[0] == 1.98 * u.um
+        assert abs(new_spec.waverange[0] == 1.98 * u.um) < 1e-10 * u.AA
         assert new_spec(1.98 * u.um).value == approx(12.8)
 
     @pytest.mark.parametrize(("endpoint", "msg"),