diff --git a/cleo/light/light_dependence.py b/cleo/light/light_dependence.py
index 5acd14fa..d995daf8 100644
--- a/cleo/light/light_dependence.py
+++ b/cleo/light/light_dependence.py
@@ -55,6 +55,11 @@ def log_makima_interpolator(lambdas_nm, epsilons, lambda_new_nm):
     return _log_(makima_interpolator, lambdas_nm, epsilons, lambda_new_nm)
 
 
+def _sorted_spectrum(spectrum):
+    spectrum_array = np.array(spectrum)
+    return spectrum_array[spectrum_array[:, 0].argsort()].T
+
+
 # hacky MRO stuff...multiple inheritance only works because slots=False,
 # and must be placed *before* SynapseDevice to work right
 @define(eq=False, slots=False)
@@ -122,8 +127,7 @@ def _get_source_for_synapse(
     def epsilon(self, lambda_new: Quantity) -> float:
         """Returns the :math:`\\varepsilon` value for a given lambda (including units)
         representing the relative sensitivity of the opsin to that wavelength."""
-        lam_eps_array = np.array(self.spectrum)
-        lambdas, epsilons = lam_eps_array[lam_eps_array[:, 0].argsort()].T
+        lambdas, epsilons = _sorted_spectrum(self.spectrum)
         lambda_new /= nmeter
         eps_new = self.spectrum_interpolator(lambdas, epsilons, lambda_new)
 
@@ -205,7 +209,7 @@ def plot_spectra(
 
     fig, ax = plt.subplots()
     for ldd in ldds:
-        lambdas, epsilons = np.array(ldd.spectrum).T
+        lambdas, epsilons = _sorted_spectrum(ldd.spectrum)
         i_data_in_range = (lambdas > xlim[0]) & (lambdas < xlim[1])
         lambdas, epsilons = lambdas[i_data_in_range], epsilons[i_data_in_range]
         if not extrapolate:
diff --git a/tests/light/test_light.py b/tests/light/test_light.py
index 2bf57783..caffe9da 100644
--- a/tests/light/test_light.py
+++ b/tests/light/test_light.py
@@ -1,10 +1,19 @@
+import matplotlib.pyplot as plt
 import neo
 import pytest
 import quantities as pq
 from brian2 import Network, NeuronGroup, asarray, mm, mm2, ms, mwatt, nmeter, np, um
 
 import cleo
-from cleo.light import GaussianEllipsoid, KoehlerBeam, Light, LightModel, fiber473nm
+from cleo.light import (
+    GaussianEllipsoid,
+    KoehlerBeam,
+    Light,
+    LightModel,
+    fiber473nm,
+    plot_spectra,
+)
+from cleo.opto import chr2_4s
 from cleo.utilities import normalize_coords, unit_safe_allclose
 
 
@@ -106,6 +115,27 @@ def test_reset():
     assert light.value == 0
 
 
+def test_plot_spectra_sorts_wavelengths():
+    sorted_opsin = chr2_4s()
+    unsorted_opsin = chr2_4s()
+    unsorted_opsin.spectrum = unsorted_opsin.spectrum[::-1]
+    original_spectrum = unsorted_opsin.spectrum.copy()
+
+    sorted_fig, sorted_ax = plot_spectra(sorted_opsin)
+    unsorted_fig, unsorted_ax = plot_spectra(unsorted_opsin)
+
+    np.testing.assert_allclose(
+        unsorted_ax.lines[0].get_xdata(), sorted_ax.lines[0].get_xdata()
+    )
+    np.testing.assert_allclose(
+        unsorted_ax.lines[0].get_ydata(), sorted_ax.lines[0].get_ydata()
+    )
+    assert unsorted_opsin.spectrum == original_spectrum
+
+    plt.close(sorted_fig)
+    plt.close(unsorted_fig)
+
+
 def test_coords():
     light = Light(light_model=fiber473nm(), coords=[[0, 0, 0], [1, 1, 1]] * mm)
     assert light.n == 2