diff --git a/abel/__init__.py b/abel/__init__.py
index 553baebf..1dd91641 100644
--- a/abel/__init__.py
+++ b/abel/__init__.py
@@ -42,6 +42,7 @@
 from .classes.source.impl.source_basic import *
 from .classes.source.impl.source_trapezoid import SourceTrapezoid
 from .classes.source.impl.source_combiner import *
+from .classes.source.impl.source_energy_comb import *
 from .classes.source.impl.source_from_file import *
 from .classes.source.impl.source_flattop import *
 from .classes.source.impl.source_capsule import *
@@ -72,6 +73,8 @@
 from .classes.bds.impl.bds_fbt import *
 from .classes.spectrometer.impl.spectrometer_facet_ocelot import *
 from .classes.spectrometer.impl.spectrometer_basic_clear import *
+from .classes.spectrometer.impl.achromatic.spectrometer_achromatic import SpectrometerAchromatic
+from .classes.spectrometer.impl.achromatic.impl.spectrometer_achromatic_basic import SpectrometerAchromaticBasic
 from .classes.ip.impl.ip_basic import *
 from .classes.ip.impl.ip_guineapig import *
 
diff --git a/abel/classes/beam.py b/abel/classes/beam.py
index a51a273e..ff260b5f 100644
--- a/abel/classes/beam.py
+++ b/abel/classes/beam.py
@@ -1364,7 +1364,7 @@ def plot_trace_space_y(self):
         cb = fig.colorbar(p)
         cb.ax.set_ylabel('Charge density (pC/um/mrad)')
 
-    def plot_transverse_profile(self):
+    def plot_transverse_profile(self, xlims=(None, None), ylims=(None, None), title=''):
         dQdxdy, xs, ys = self.phase_space_density(self.xs, self.ys)
 
         fig, ax = plt.subplots()
@@ -1379,6 +1379,9 @@ def plot_transverse_profile(self):
         cb = fig.colorbar(p)
         cb.ax.set_ylabel('Charge density (pC/um^2)')
 
+        ax.set_xlim(xlims)
+        ax.set_ylim(ylims)
+        ax.set_title(title)
     
     # TODO: unfinished!
     # def plot_bunch_pattern(self):
diff --git a/abel/classes/beamline/impl/experiment/experiment.py b/abel/classes/beamline/impl/experiment/experiment.py
index b4dff5bc..08cdf4e1 100644
--- a/abel/classes/beamline/impl/experiment/experiment.py
+++ b/abel/classes/beamline/impl/experiment/experiment.py
@@ -9,33 +9,34 @@
 import warnings
 
 class Experiment(Beamline):
-    
+
     def __init__(self, linac=None, component=None, spectrometer=None, num_bunches_in_train=1, bunch_separation=0, rep_rate_trains=10):
         super().__init__(num_bunches_in_train, bunch_separation, rep_rate_trains)
-        
+
         self.linac = linac
         self.component = component
         self.spectrometer = spectrometer
-        
-    
-    
+
+
+
     # assemble the trackables
     def assemble_trackables(self):
 
         self.trackables = []
-        
+
         # add the linac/source
         assert(isinstance(self.linac, Linac) or isinstance(self.linac, Source))
         self.trackables.append(self.linac)
 
         # add the experimental component
-        self.trackables.append(self.component)
+        if self.component is not None:
+            self.trackables.append(self.component)
 
         # add the spectrometer (if it exists
         if self.spectrometer is not None:
             assert(isinstance(self.spectrometer, Spectrometer))
             self.trackables.append(self.spectrometer)
-        
+
         # set the bunch train pattern etc.
         super().assemble_trackables()
 
@@ -45,10 +46,10 @@ def energy_usage(self):
 
     def get_cost_breakdown():
         return ('Experiment', 0)
-    
+
     # density plots
     def plot_spectrometer_screen(self, xlims=None, ylims=None, E_calib = False, diverg = None, plot_m12 = False, savefig = None):
-        
+
         # load phase space
         beam = self.final_beam
 
@@ -60,15 +61,15 @@ def plot_spectrometer_screen(self, xlims=None, ylims=None, E_calib = False, dive
         if ylims is not None:
             ybins = np.linspace(min(ylims), max(ylims), num_bins)
         dQdxdy, xedges, yedges = beam.density_transverse(hbins=xbins, vbins=ybins)
-        
+
         # prepare figure
         fig, ax = plt.subplots(2, 1, gridspec_kw={'height_ratios': [2, 1]})
         fig.set_figwidth(8)
         fig.set_figheight(8)
-        
+
         # current profile
         c0 = ax[0].pcolor(xedges*1e3, yedges * 1e3, abs(dQdxdy) * 1e3, cmap=CONFIG.default_cmap, shading='auto')
-        
+
         # make plot
         ax[0].set_xlabel('x (mm)')
         ax[0].set_ylabel('y (mm)')
@@ -77,25 +78,25 @@ def plot_spectrometer_screen(self, xlims=None, ylims=None, E_calib = False, dive
         cbar0 = fig.colorbar(c0, ax=ax[0], pad = 0.015*fig.get_figwidth())
         cbar0.ax.set_ylabel('Charge density (nC ' r'$\mathrm{mm^{-2})}$')
         ax[0].set_ylim(min(yedges * 1e3), max(yedges * 1e3))
-        
+
         # calculate energy axis (E = E0*Dy/y)
         Dy_img = self.spectrometer.get_dispersion(energy=self.spectrometer.img_energy)
         E_times_Dy = self.spectrometer.img_energy*Dy_img
-        
+
         # add imaging energies
         y_img = -E_times_Dy/self.spectrometer.img_energy
         y_img_y = -E_times_Dy/self.spectrometer.img_energy_y
         ax[0].axhline(y_img*1e3, color = 'black', linestyle = '--', linewidth = 1, label = 'Imaging energy x (GeV)')
         ax[0].axhline(y_img_y*1e3, color = 'black', linestyle = ':', linewidth = 1, label = 'Imaging energy y (GeV)')
-        
+
         # add energy axis
         warnings.simplefilter('ignore', category=RuntimeWarning)
         ax2 = ax[0].secondary_yaxis('right', functions=(lambda y: -E_times_Dy/(y/1e3)/1e9, lambda y: -E_times_Dy/(y/1e3)/1e9))
         ax2.set_ylabel('Energy (GeV)')
-        
+
         if diverg is not None:
             energies = -E_times_Dy/yedges
-            m12s = self.spectrometer.get_m12(energies)        
+            m12s = self.spectrometer.get_m12(energies)
             x_pos = m12s*diverg
             x_neg = -m12s*diverg
             ax[0].plot(x_pos*1e3, yedges*1e3, color = 'orange', label = str(diverg*1e3) + ' mrad butterfly', alpha = 0.7)
@@ -104,12 +105,11 @@ def plot_spectrometer_screen(self, xlims=None, ylims=None, E_calib = False, dive
             ax[1].set_ylabel(r'$\mathrm{m_{12}}$')
             ax[1].set_xlabel('E (GeV)')
             ax[1].grid(False, which='major')
-            
+
         ax[0].set_xlim(min(xedges*1e3), max(xedges*1e3))
         ax[0].legend(loc = 'center right', fontsize = 6)
-        
+
         if plot_m12 == False:
             fig.delaxes(ax[1])
         if savefig is not None:
             plt.savefig(str(savefig), dpi = 700)
-    
\ No newline at end of file
diff --git a/abel/classes/plasma_lens/impl/plasma_lens_nonlinear_thin.py b/abel/classes/plasma_lens/impl/plasma_lens_nonlinear_thin.py
index a0a996dc..9aa5ccb9 100644
--- a/abel/classes/plasma_lens/impl/plasma_lens_nonlinear_thin.py
+++ b/abel/classes/plasma_lens/impl/plasma_lens_nonlinear_thin.py
@@ -2,28 +2,52 @@
 from matplotlib import patches
 import numpy as np
 import scipy.constants as SI
+import matplotlib.pyplot as plt
 
 class PlasmaLensNonlinearThin(PlasmaLens):
-    
-    def __init__(self, length=None, radius=None, current=None, rel_nonlinearity=0):
+
+    def __init__(self, length=None, radius=None, current=None, rel_nonlinearity=0, nonlinearity_in_x=True):
 
         super().__init__(length, radius, current)
-        
+
         # set nonlinearity (defined as R/Dx)
         self.rel_nonlinearity = rel_nonlinearity
-        
-    
+        self.nonlinearity_in_x = nonlinearity_in_x
+
+    def remove_charge_outside(self, beam, plot=False):
+
+        mask = np.sqrt((beam.xs()-self.offset_x)**2 + (beam.ys()-self.offset_y)**2) > self.radius
+
+        if plot:
+
+            fig, ax = plt.subplots()
+            ax.add_patch(patches.Circle((self.offset_x, self.offset_y), self.radius, fill=False, linestyle='--'))
+
+            ax.scatter(beam.xs(), beam.ys(), s=1, alpha=0.1, zorder=0)
+            ax.set_aspect('equal', 'box')
+            plt.show()
+
+            rel_loss = mask.sum() / mask.size
+            print(f"Removing {rel_loss*100:.2f}% of the charge outside the lens")
+
+        # remove charge outside the lens
+        del beam[mask]
+
+        return beam
+
     def track(self, beam, savedepth=0, runnable=None, verbose=False):
 
         # remove charge outside the lens (start)
-        del beam[np.sqrt((beam.xs()-self.offset_x)**2 + (beam.ys()-self.offset_y)**2) > self.radius]
-        
+        # del beam[np.sqrt((beam.xs()-self.offset_x)**2 + (beam.ys()-self.offset_y)**2) > self.radius]
+        self.remove_charge_outside(beam, plot=False)
+
         # drift half the distance
         beam.transport(self.length/2)
-        
+
         # remove charge outside the lens (middle)
-        del beam[np.sqrt((beam.xs()-self.offset_x)**2 + (beam.ys()-self.offset_y)**2) > self.radius]
-        
+        # del beam[np.sqrt((beam.xs()-self.offset_x)**2 + (beam.ys()-self.offset_y)**2) > self.radius]
+        self.remove_charge_outside(beam)
+
         # get particles
         xs = beam.xs()-self.offset_x
         xps = beam.xps()
@@ -34,16 +58,21 @@ def track(self, beam, savedepth=0, runnable=None, verbose=False):
         g0 = self.get_focusing_gradient()
 
         # calculate the nonlinearity
-        inv_Dx = self.rel_nonlinearity/self.radius
-        
+        if self.nonlinearity_in_x:
+            inv_Dx = self.rel_nonlinearity/self.radius
+            inv_Dy = 0.0
+        else:
+            inv_Dx = 0.0
+            inv_Dy = self.rel_nonlinearity/self.radius
+
         # thin lens kick
-        Bx = g0*(ys + xs*ys*inv_Dx)
-        By = -g0*(xs + (xs**2 + ys**2)/2*inv_Dx)
+        Bx = g0*(ys + xs*ys*inv_Dx  + (xs**2 + ys**2)/2*inv_Dy)
+        By = -g0*(xs + (xs**2 + ys**2)/2*inv_Dx + xs*ys*inv_Dy)
 
         # calculate the angular kicks
         delta_xp = self.length*(By*beam.charge_sign()*SI.c/beam.Es())
         delta_yp = -self.length*(Bx*beam.charge_sign()*SI.c/beam.Es())
-        
+
         # set new beam positions and angles (shift back plasma-lens offsets)
         beam.set_xps(xps + delta_xp)
         beam.set_yps(yps + delta_yp)
@@ -52,11 +81,11 @@ def track(self, beam, savedepth=0, runnable=None, verbose=False):
         beam.transport(self.length/2)
 
         # remove charge outside the lens (end)
-        del beam[np.sqrt((beam.xs()-self.offset_x)**2 + (beam.ys()-self.offset_y)**2) > self.radius]
-        
-        return super().track(beam, savedepth, runnable, verbose)   
+        # del beam[np.sqrt((beam.xs()-self.offset_x)**2 + (beam.ys()-self.offset_y)**2) > self.radius]
+        self.remove_charge_outside(beam)
+
+        return super().track(beam, savedepth, runnable, verbose)
 
 
     def get_focusing_gradient(self):
         return SI.mu_0 * self.current / (2*np.pi * self.radius**2)
-    
\ No newline at end of file
diff --git a/abel/classes/plasma_lens/plasma_lens.py b/abel/classes/plasma_lens/plasma_lens.py
index 26d55134..cc36c92e 100644
--- a/abel/classes/plasma_lens/plasma_lens.py
+++ b/abel/classes/plasma_lens/plasma_lens.py
@@ -4,12 +4,12 @@
 import numpy as np
 
 class PlasmaLens(Trackable):
-    
+
     @abstractmethod
     def __init__(self, length, radius, current, offset_x=0, offset_y=0):
 
         super().__init__()
-        
+
         # common variables
         self.length = length
         self.radius = radius
@@ -17,8 +17,8 @@ def __init__(self, length, radius, current, offset_x=0, offset_y=0):
 
         self.offset_x = offset_x
         self.offset_y = offset_y
-        
-    
+
+
     @abstractmethod
     def track(self, beam, savedepth=0, runnable=None, verbose=False):
         return super().track(beam, savedepth, runnable, verbose)
@@ -26,10 +26,13 @@ def track(self, beam, savedepth=0, runnable=None, verbose=False):
     @abstractmethod
     def get_focusing_gradient(self):
         pass
-        
+
     def get_length(self):
         return self.length
 
     def survey_object(self):
         return patches.Rectangle((0, -1), self.get_length(), 2)
-    
\ No newline at end of file
+
+    @property
+    def l(self):
+        return self.length
diff --git a/abel/classes/source/impl/source_energy_comb.py b/abel/classes/source/impl/source_energy_comb.py
new file mode 100644
index 00000000..c59bb94c
--- /dev/null
+++ b/abel/classes/source/impl/source_energy_comb.py
@@ -0,0 +1,317 @@
+"""
+Module: source_energy_comb
+
+This module defines the `SourceEnergyComb` class, which is a subclass
+of `Source`. It provides methods to initialize the source and track
+the beam with comb-like energy distributions.
+
+Classes
+-------
+`SourceEnergyComb` : `Source`
+    A source class for creating comb-like energy distributions.
+
+Imports
+-------
+import time
+import numpy as np
+import scipy.constants as SI
+from abel import Source, Beam
+from abel.utilities.beam_physics import generate_trace_space_xy, generate_symm_trace_space_xyz
+from abel.utilities.relativity import energy2gamma
+"""
+
+import time
+import numpy as np
+import scipy.constants as SI
+from abel import Source, Beam
+from abel.utilities.beam_physics import generate_trace_space_xy, generate_symm_trace_space_xyz
+from abel.utilities.relativity import energy2gamma
+
+class SourceEnergyComb(Source):
+    """
+    A source class that extends the functionality of `SourceBasic` by distributing
+    particle energies within a given list of energies. This is useful for creating
+comb-like energy distributions for resolution estimates.
+
+    Parameters
+    ----------
+    length : [m] float, optional
+        Length of the source. Default is 0.
+
+    num_particles : int, optional
+        Number of particles in the beam. Default is 1000.
+
+    energy : [eV] float, optional
+        Central energy of the beam. Default is ``None``.
+
+    charge : [C] float, optional
+        Total charge of the beam. Default is 0.
+
+    rel_energy_spread : float, optional
+        Relative energy spread (fraction). Default is ``None``.
+
+    energy_spread : [eV] float, optional
+        Absolute energy spread. Default is ``None``.
+
+    bunch_length : [m] float, optional
+        RMS bunch length. Default is ``None``.
+
+    z_offset : [m] float, optional
+        Longitudinal offset of the beam. Default is 0.
+
+    x_offset, y_offset : [m] float, optional
+        Transverse offsets of the beam. Default is 0.
+
+    x_angle, y_angle : [rad] float, optional
+        Transverse angles of the beam. Default is 0.
+
+    emit_nx, emit_ny : [m rad] float, optional
+        Normalized emittances in x and y. Default is 0.
+
+    beta_x, beta_y : [m] float, optional
+        Beta functions in x and y. Default is ``None``.
+
+    alpha_x, alpha_y : float, optional
+        Alpha functions in x and y. Default is 0.
+
+    angular_momentum : float, optional
+        Angular momentum of the beam. Default is 0.
+
+    wallplug_efficiency : float, optional
+        Efficiency of the source. Default is 1.
+
+    accel_gradient : [V/m] float, optional
+        Acceleration gradient. Default is ``None``.
+
+    symmetrize : bool, optional
+        Whether to symmetrize the transverse phase space. Default is False.
+
+    symmetrize_6d : bool, optional
+        Whether to symmetrize the full 6D phase space. Default is False.
+
+    z_cutoff : [m] float, optional
+        Longitudinal cutoff for filtering particles. Default is ``None``.
+
+    energy_comb_delta : float, optional
+        Relative energy step size for creating comb-like energy distributions. Default is ``None``.
+    """
+
+    def __init__(self, length=0, num_particles=1000, energy=None, charge=0, rel_energy_spread=None, energy_spread=None, bunch_length=None, z_offset=0, x_offset=0, y_offset=0, x_angle=0, y_angle=0, emit_nx=0, emit_ny=0, beta_x=None, beta_y=None, alpha_x=0, alpha_y=0, angular_momentum=0, wallplug_efficiency=1, accel_gradient=None, symmetrize=False, symmetrize_6d=False, z_cutoff=None, energy_comb_delta=None):
+
+        super().__init__(length, charge, energy, accel_gradient, wallplug_efficiency, x_offset, y_offset, x_angle, y_angle)
+
+        self.rel_energy_spread = rel_energy_spread # [eV]
+        self.energy_spread = energy_spread
+        self.bunch_length = bunch_length # [m]
+        self.z_offset = z_offset # [m]
+        self.num_particles = num_particles
+
+        self.emit_nx = emit_nx # [m rad]
+        self.emit_ny = emit_ny # [m rad]
+        self.beta_x = beta_x # [m]
+        self.beta_y = beta_y # [m]
+        self.alpha_x = alpha_x # [m]
+        self.alpha_y = alpha_y # [m]
+
+        self.angular_momentum = angular_momentum
+        self.symmetrize = symmetrize
+        self.symmetrize_6d = symmetrize_6d
+        self.z_cutoff = z_cutoff
+
+        # TODO: Make this an energy_diff or delta instead of a linspace
+        # the diff should be in percentage
+        # this will be easier for scans
+        self.energy_comb_delta = energy_comb_delta
+
+
+    def track(self, _=None, savedepth=0, runnable=None, verbose=False):
+        """
+        Tracks the beam through the source, generating its phase space and applying
+        energy comb distributions or filters if specified.
+
+        Parameters
+        ----------
+        _ : ``None``, optional
+            Placeholder parameter. Default is ``None``.
+
+        savedepth : int, optional
+            Depth of saved tracking data. Default is 0.
+
+        runnable : object, optional
+            Runnable object for tracking. Default is ``None``.
+
+        verbose : bool, optional
+            Whether to print verbose output. Default is False.
+
+        Returns
+        -------
+        Beam
+            The tracked beam object.
+        """
+
+        # make empty beam
+        beam = Beam()
+
+        # Lorentz gamma
+        gamma = energy2gamma(self.energy)
+
+        # generate relative/absolute energy spreads
+        if self.rel_energy_spread is not None:
+            if self.energy_spread is None:
+                self.energy_spread = self.energy * self.rel_energy_spread
+            elif abs(self.energy_spread - self.energy * self.rel_energy_spread) > 0:
+                raise Exception("Both absolute and relative energy spread defined.")
+
+        if self.symmetrize_6d is False:
+
+            # horizontal and vertical phase spaces
+            xs, xps, ys, yps = generate_trace_space_xy(self.emit_nx/gamma, self.beta_x, self.alpha_x, self.emit_ny/gamma, self.beta_y, self.alpha_y, self.num_particles, self.angular_momentum/gamma, symmetrize=self.symmetrize)
+
+            # longitudinal phase space
+            if self.symmetrize:
+                num_tiling = 4
+                num_particles_actual = round(self.num_particles/num_tiling)
+            else:
+                num_particles_actual = self.num_particles
+            zs = np.random.normal(loc=self.z_offset, scale=self.bunch_length, size=num_particles_actual)
+            # Es = np.random.normal(loc=self.energy, scale=self.energy_spread, size=num_particles_actual)
+            # this should yield a uniform energy distribution with length 2*energy_spread around self.energy
+            Es = (np.random.random(size=num_particles_actual) - 0.5) * 2 * self.energy_spread + self.energy
+            if self.symmetrize:
+                zs = np.tile(zs, num_tiling)
+                Es = np.tile(Es, num_tiling)
+
+        else:
+            xs, xps, ys, yps, zs, Es = generate_symm_trace_space_xyz(self.emit_nx/gamma, self.beta_x, self.alpha_x, self.emit_ny/gamma, self.beta_y, self.alpha_y, self.num_particles, self.bunch_length, self.energy_spread, self.angular_momentum/gamma)
+
+            # add longitudinal offsets
+            zs += self.z_offset
+            Es += self.energy
+
+        # create phase space
+        beam.set_phase_space(xs=xs, ys=ys, zs=zs, xps=xps, yps=yps, Es=Es, Q=self.charge)
+
+        # Apply filter(s) if desired
+        if self.z_cutoff is not None:
+            beam = self.z_filter(beam)
+
+        if self.energy_comb_delta is not None:
+            self.reorder_particle_energies(beam)
+
+        # add jitters and offsets in super function
+        return super().track(beam, savedepth, runnable, verbose)
+
+    def reorder_particle_energies(self, beam, plot=False):
+        """
+        Reorders the particle energies in the beam to create a comb-like energy
+        distribution based on `self.energy_comb_delta`.
+
+        Parameters
+        ----------
+        beam : Beam
+            The beam object whose particle energies will be reordered.
+
+        plot : bool, optional
+            Whether to plot the longitudinal phase space after reordering. Default is False.
+
+        Returns
+        -------
+        ``None``
+        """
+
+        # takes in the beam and uses the attribute self.energy_comb_delta
+        # to create a list of energies that will be given to each beam particle
+        # to create a comb-like energy distribution for resolution estimates
+
+        # the energies of all bunch particles
+        vals = beam.Es()
+
+        energy_beam = beam.energy()
+
+        # the energy delta is given as relative, thus multiply times the average
+        # particle energy
+        energy_delta = self.energy_comb_delta * energy_beam
+
+        # half of the number of steps depending on the step size and the energy
+        # range of the beam particles
+        num_steps_half = int((vals.max() - vals.min()) * 0.5
+                             / energy_delta  + 1)
+
+        # go half the number of steps in both directions, with the bunch
+        # average energy in the middle
+        energy_steps = np.arange(-num_steps_half, num_steps_half+1)
+        energies_comb = energy_beam + energy_steps * energy_delta
+
+        # find the index for each value that shows the closest value in energies
+        # this function requires "left" or "right" to be specified
+        # use "left"
+        idx = np.searchsorted(energies_comb, vals, side="left")
+
+        # this is required as otherwise the energies_comb[idx] in the next line
+        # will throw an error
+        idx[idx == len(energies_comb)] = len(energies_comb) - 1
+
+        # find out if the value with the left-index or the value with the
+        # right-index value is closer and accordingly subtract the index value
+        idx = idx - (np.abs(vals - energies_comb[idx-1])
+                     < np.abs(vals - energies_comb[idx]))
+
+        # energies_comb[idx] has all the energies that were closest for each
+        # particle.  Set the particle energies to this new energies
+        beam.set_Es(energies_comb[idx])
+
+        if plot:
+            # print(np.diff(energies_comb)*1E-9)
+            beam.plot_lps()
+
+        return
+
+
+    # ==================================================
+    # Filter out particles whose z < z_cutoff for testing instability etc.
+    def z_filter(self, beam):
+        """
+        Filters out particles whose longitudinal positions (z) are below the
+        specified cutoff (`self.z_cutoff`).
+
+        Parameters
+        ----------
+        beam : Beam
+            The beam object to be filtered.
+
+        Returns
+        -------
+        Beam
+            A new beam object with particles filtered based on the z cutoff.
+        """
+        xs = beam.xs()
+        ys = beam.ys()
+        zs = beam.zs()
+        pxs = beam.pxs()
+        pys = beam.pys()
+        pzs = beam.pzs()
+        weights = beam.weightings()
+
+        # Apply the filter
+        bool_indices = (zs > self.z_cutoff)
+        zs_filtered = zs[bool_indices]
+        xs_filtered = xs[bool_indices]
+        ys_filtered = ys[bool_indices]
+        pxs_filtered = pxs[bool_indices]
+        pys_filtered = pys[bool_indices]
+        pzs_filtered = pzs[bool_indices]
+        weights_filtered = weights[bool_indices]
+
+        # Initialise ABEL Beam object
+        beam_out = Beam()
+
+        # Set the phase space of the ABEL beam
+        beam_out.set_phase_space(Q=np.sum(weights_filtered)*-SI.e,
+                             xs=xs_filtered,
+                             ys=ys_filtered,
+                             zs=zs_filtered,
+                             pxs=pxs_filtered,  # Always use single particle momenta?
+                             pys=pys_filtered,
+                             pzs=pzs_filtered)
+
+        return beam_out
diff --git a/abel/classes/spectrometer/impl/achromatic/impl/spectrometer_achromatic_basic.py b/abel/classes/spectrometer/impl/achromatic/impl/spectrometer_achromatic_basic.py
new file mode 100644
index 00000000..74595b97
--- /dev/null
+++ b/abel/classes/spectrometer/impl/achromatic/impl/spectrometer_achromatic_basic.py
@@ -0,0 +1,313 @@
+"""
+Module: spectrometer_achromatic_basic
+
+This module defines the `SpectrometerAchromaticBasic` class, which is a
+subclass of `SpectrometerAchromatic`. It provides methods to initialize
+the spectrometer and track the beam.
+
+Classes
+-------
+`SpectrometerAchromaticBasic` : `SpectrometerAchromatic`
+    A basic implementation of an achromatic spectrometer.
+
+Imports
+-------
+from abel import SpectrometerAchromatic, PlasmaLensNonlinearThin
+import scipy.constants as SI
+import numpy as np
+"""
+
+from abel import SpectrometerAchromatic, PlasmaLensNonlinearThin
+import scipy.constants as SI
+import numpy as np
+
+class SpectrometerAchromaticBasic(SpectrometerAchromatic):
+    """
+    A basic implementation of an achromatic spectrometer.
+
+    Parameters
+    ----------
+    length_drift_to_dipole : [m] float, optional
+        Length of the drift space to the dipole magnet. Default is 0.2.
+
+    field_dipole : [T] float, optional
+        Magnetic field strength of the dipole magnet. Negative values
+        bend the beam downwards. Default is -0.85.
+
+    length_dipole : [m] float, optional
+        Length of the dipole magnet. Default is 0.15.
+
+    length_drift_dipole_to_lens : [m] float, optional
+        Length of the drift space from the dipole magnet to the lens.
+        Default is 0.15.
+
+    length_plasma_lens : [m] float, optional
+        Length of the plasma lens. Default is 20E-3.
+
+    radius_plasma_lens : [m] float, optional
+        Radius of the plasma lens. Default is 1E-3.
+
+    length_drift_lens_to_screen : [m] float, optional
+        Length of the drift space from the lens to the screen. Default
+        is 0.4.
+
+    imaging_energy : [eV] float or None, optional
+        Imaging energy. Default is None.
+
+    disable_lens_nonlinearity : bool, optional
+        Flag to disable lens nonlinearity. Default is `True`.
+    """
+
+    def __init__(self,
+                 length_drift_to_dipole=0.2,
+                 field_dipole=-0.85,
+                 length_dipole=0.15,
+                 length_drift_dipole_to_lens=0.15,
+                 length_plasma_lens=20E-3,
+                 radius_plasma_lens=1E-3,
+                 length_drift_lens_to_screen=0.4,
+                 imaging_energy=None,
+                 disable_lens_nonlinearity=True):
+        """
+        Initializes the `SpectrometerAchromaticBasic` with the given
+        parameters.
+
+        Parameters
+        ----------
+        length_drift_to_dipole : [m] float, optional
+            Length of the drift space to the dipole magnet. Default is
+            0.2.
+
+        field_dipole : [T] float, optional
+            Magnetic field strength of the dipole magnet. Negative
+            values bend the beam downwards. Default is -0.85.
+
+        length_dipole : [m] float, optional
+            Length of the dipole magnet. Default is 0.15.
+
+        length_drift_dipole_to_lens : [m] float, optional
+            Length of the drift space from the dipole magnet to the lens.
+            Default is 0.15.
+
+        length_plasma_lens : [m] float, optional
+            Length of the plasma lens. Default is 20E-3.
+
+        radius_plasma_lens : [m] float, optional
+            Radius of the plasma lens. Default is 1E-3.
+
+        length_drift_lens_to_screen : [m] float, optional
+            Length of the drift space from the lens to the screen.
+            Default is 0.4.
+
+        imaging_energy : [eV] float or None, optional
+            Imaging energy. Default is None.
+
+        disable_lens_nonlinearity : bool, optional
+            Flag to disable lens nonlinearity. Default is `True`.
+        """
+
+        super().__init__(length_drift_to_dipole=length_drift_to_dipole,
+                         field_dipole=field_dipole,
+                         length_dipole=length_dipole,
+                         length_drift_dipole_to_lens=length_drift_dipole_to_lens,
+                         length_plasma_lens=length_plasma_lens,
+                         radius_plasma_lens=radius_plasma_lens,
+                         length_drift_lens_to_screen=length_drift_lens_to_screen,
+                         imaging_energy=imaging_energy,
+                         disable_lens_nonlinearity=disable_lens_nonlinearity)
+
+    def track(self, beam0, savedepth=0, runnable=None, verbose=False):
+        """
+        Tracks the beam through the spectrometer.
+
+        Parameters
+        ----------
+        beam0 : `Beam` object
+            The initial beam to be tracked.
+
+        savedepth : int, optional
+            The depth at which to save the beam state. Default is 0.
+
+        runnable : callable, optional
+            A callable to be run during tracking. Default is ``None``.
+
+        verbose : bool, optional
+            Whether to print verbose output. Default is `False`.
+
+        Returns
+        -------
+        `Beam` object
+            The tracked beam object.
+        """
+
+        # transport the beam to the center of the dipole
+        beam0.transport(self.length_drift_to_dipole + self.length_dipole/2)
+
+        # this function allows to plot the trace space for debugging
+        # beam0.plot_trace_space_y()
+
+        # computing the bending angle of the dipole for each particle
+        # sin(bend_angle/2) = length_dipole*0.5 / bending_radius
+        # bending_radius = momentum / (B * e)
+        # for small angles sin(bend_angle/2) \approx bend_angle/2
+        # bend_angle = length_dipole * B * e / momentum
+        # momentum [eV/c]
+        bend_angles = self.get_bending_angle(beam0.Es())
+        # bend_angles = self.length_dipole * self.field_dipole * SI.c / beam0.Es()
+
+        # if self.imaging_energy is None:
+        #     energy_nom = beam0.energy()
+        # else:
+        energy_nom = beam0.energy()
+
+        # the bending angle for the nominal beam energy
+        bend_angle_nom = self.get_bending_angle(energy_nom)
+        # bend_angle_nom = self.length_dipole * self.field_dipole * SI.c / energy_nom
+
+        # the bending angle of each particle relative to the bending of
+        # the beam at nominal energy
+        bend_angles_rel = bend_angles - bend_angle_nom
+
+        # updating the angles of the particles with the dipole kick
+        beam0.set_yps(beam0.yps() + bend_angles_rel)
+
+        # transport the beam to the lens with the new angles of the particles
+        beam0.transport(self.length_dipole/2 + self.length_drift_dipole_to_lens)
+
+        # calculate the focal length
+        length_object_to_lens = (self.length_drift_to_dipole
+                                 + self.length_dipole
+                                 + self.length_drift_dipole_to_lens
+                                 + self.length_plasma_lens/2)
+        length_lens_to_screen = (self.length_plasma_lens/2
+                                 + self.length_drift_lens_to_screen)
+        focal_length = 1 / (1 / length_object_to_lens
+                            + 1 / length_lens_to_screen)
+
+        # calculate the lens strength based on the focal length
+        k = 1 / (focal_length * self.length_plasma_lens)
+
+        if self.imaging_energy is None:
+            imaging_energy = beam0.energy()
+        else:
+            imaging_energy = self.imaging_energy
+
+        # calculate the magnetic field gradient with the imaging energy
+        g = k / (SI.c/ imaging_energy)
+
+        # calculate the current based on the magnetic field gradient
+        current = -g * 2 * np.pi * self.radius_plasma_lens**2 / (SI.mu_0)
+
+        # calculate the dispersion introduced by the dipole at the position of
+        # the plasma lens
+        length_middle_dipole_to_lens = (self.length_dipole/2
+                                        + self.length_drift_dipole_to_lens
+                                        + self.length_plasma_lens/2)
+        dispersion = self.get_dispersion(bend_angle=bend_angle_nom,
+                                         length=length_middle_dipole_to_lens)
+
+        # calculate the nonlinearity of the plasma lens, which is
+        # defined by R/Dy (negative sign found by trial and error)
+        if not self.disable_lens_nonlinearity:
+            rel_nonlinearity = - self.radius_plasma_lens / dispersion
+        else:
+            rel_nonlinearity = 0.0
+
+        # use the already implemented PlasmaLensNonlinearThin class to calculate
+        # the kick produced by the plasma lens. The kicks are updated in the
+        # class beam
+        # nonlinearity_in_x=False as the dipole here disperses vertically
+        plasma_lens = PlasmaLensNonlinearThin(length=self.length_plasma_lens,
+                                              radius=self.radius_plasma_lens,
+                                              current=current,
+                                              rel_nonlinearity=rel_nonlinearity,
+                                              nonlinearity_in_x=False)
+
+        beam = plasma_lens.track(beam=beam0)
+
+        # transport the beam with the updated angles
+        length_middle_lens_to_screen = (self.length_plasma_lens/2
+                                        + self.length_drift_lens_to_screen)
+        beam.transport(length_middle_lens_to_screen)
+
+        # shift beam by the dispersion (special case for spectrometers)
+        # this will allow to plot the beam as if there was a screen
+        length_middle_dipole_to_screen = (self.length_dipole / 2
+                                           + self.length_drift_dipole_to_lens
+                                           + self.length_plasma_lens
+                                           + self.length_drift_lens_to_screen)
+        # dispersion_on_screen = bend_angle_nom * length_middle_dipole_to_screen
+        dispersion_on_screen = self.get_dispersion(bend_angle=bend_angle_nom,
+                                                   length= length_middle_dipole_to_screen)
+        beam.set_ys(beam.ys() + dispersion_on_screen)
+
+        return super().track(beam, savedepth, runnable, verbose)
+
+    def get_bending_angle(self, energy):
+        """
+        Calculates the bending angle of the spectrometer for a given energy.
+
+        Parameters
+        ----------
+        energy : [eV] float
+            The energy of the particle beam.
+
+        Returns
+        -------
+        float
+            The calculated bending angle [rad].
+        """
+
+        # computing the bending angle of the dipole for each particle
+        # sin(bend_angle/2) = length_dipole*0.5 / bending_radius
+        # bending_radius = momentum / (B * e)
+        # for small angles sin(bend_angle/2) \approx bend_angle/2
+        # bend_angle = length_dipole * B * e / momentum
+        # momentum [eV/c]
+        bend_angle = self.length_dipole * self.field_dipole * SI.c / energy
+
+        return bend_angle
+
+    def get_dispersion(self, energy=None, bend_angle=None, length=None):
+        """
+        Calculates the dispersion of the spectrometer for a given energy or
+        bending angle.
+
+        Parameters
+        ----------
+        energy : [eV] float, optional
+            The energy of the particle beam. Either `energy` or
+            `bend_angle` must be provided.
+
+        bend_angle : [rad] float, optional
+            Bending angle of the dipole. Either `bend_angle` or
+            `energy` must be provided.
+
+        length : [m] float, optional
+            The distance to use for the calculation. If not provided,
+            it defaults to the distance from the center of the dipole to
+            the screen.
+
+        Returns
+        -------
+        float
+            The calculated dispersion value [m].
+        """
+
+        text = 'Please give either the bending angle *or* the energy.'
+        assert ((bend_angle is None and energy is not None)
+                or (bend_angle is not None and energy is None)), text
+
+        if bend_angle is None:
+            bend_angle = self.get_bending_angle(energy)
+
+        if length is None:
+            # assume lenth from center of dipole to the screen
+            length = (self.length_dipole / 2
+                      + self.length_drift_dipole_to_lens
+                      + self.length_plasma_lens
+                      + self.length_drift_lens_to_screen)
+
+        dispersion = bend_angle * length
+
+        return dispersion
diff --git a/abel/classes/spectrometer/impl/achromatic/spectrometer_achromatic.py b/abel/classes/spectrometer/impl/achromatic/spectrometer_achromatic.py
new file mode 100644
index 00000000..e0e50bdc
--- /dev/null
+++ b/abel/classes/spectrometer/impl/achromatic/spectrometer_achromatic.py
@@ -0,0 +1,202 @@
+"""
+Module: spectrometer_achromatic
+
+This module defines the `SpectrometerAchromatic` class, which is an
+abstract base class for achromatic spectrometers. It provides methods
+and properties to model the spectrometer's behavior.
+
+Classes
+-------
+`SpectrometerAchromatic` : `Spectrometer`
+    Abstract base class for achromatic spectrometers.
+
+Imports
+-------
+from abc import abstractmethod
+from abel import Spectrometer
+import numpy as np
+import matplotlib.pyplot as plt
+from abel import CONFIG
+"""
+
+from abc import abstractmethod
+from abel import Spectrometer
+import numpy as np
+import matplotlib.pyplot as plt
+from abel import CONFIG
+
+
+class SpectrometerAchromatic(Spectrometer):
+    """
+    Abstract base class for achromatic spectrometers.
+
+    This class models an achromatic spectrometer, which is used to analyze
+    particle beams by bending them through a magnetic field and focusing
+    them using a plasma lens.
+
+    Parameters
+    ----------
+    length_drift_to_dipole : [m] float
+        Length of the drift space to the dipole magnet. Default is 0.2.
+
+    field_dipole : [T] float
+        Magnetic field strength of the dipole magnet. Negative values
+        indicate bending the beam downwards. Default is -0.85.
+
+    length_dipole : [m] float
+        Length of the dipole magnet. Default is 0.8.
+
+    length_drift_dipole_to_lens : [m] float
+        Length of the drift space from the dipole magnet to the lens.
+        Default is 0.2.
+
+    length_plasma_lens : [m] float
+        Length of the plasma lens. Default is 20E-3.
+
+    radius_plasma_lens : [m] float
+        Radius of the plasma lens. Default is 1E-3.
+
+    length_drift_lens_to_screen : [m] float
+        Length of the drift space from the lens to the screen. Default
+        is 0.4.
+
+    imaging_energy : [eV] float or None
+        Imaging energy. Default is None.
+
+    disable_lens_nonlinearity : bool
+        Flag to disable lens nonlinearity. Default is ``True``.
+    """
+
+    @abstractmethod
+    def __init__(self,
+                 length_drift_to_dipole,
+                 field_dipole,
+                 length_dipole,
+                 length_drift_dipole_to_lens,
+                 length_plasma_lens,
+                 radius_plasma_lens,
+                 length_drift_lens_to_screen,
+                 imaging_energy,
+                 disable_lens_nonlinearity):
+        """
+        Initializes the `SpectrometerAchromatic` object with the specified
+        parameters.
+
+        Parameters
+        ----------
+        length_drift_to_dipole : [m] float
+            Length of the drift space to the dipole magnet. Default is 0.2.
+
+        field_dipole : [T] float
+            Magnetic field strength of the dipole magnet. Negative values
+            indicate bending the beam downwards. Default is -0.85.
+
+        length_dipole : [m] float
+            Length of the dipole magnet. Default is 0.8.
+
+        length_drift_dipole_to_lens : [m] float
+            Length of the drift space from the dipole magnet to the lens.
+            Default is 0.2.
+
+        length_plasma_lens : [m] float
+            Length of the plasma lens. Default is 20E-3.
+
+        radius_plasma_lens : [m] float
+            Radius of the plasma lens. Default is 1E-3.
+
+        length_drift_lens_to_screen : [m] float
+            Length of the drift space from the lens to the screen.
+            Default is 0.4.
+
+        imaging_energy : [eV] float or None
+            Imaging energy. Default is None.
+
+        disable_lens_nonlinearity : bool
+            Flag to disable lens nonlinearity. Default is ``True``.
+        """
+        self.length_drift_to_dipole = length_drift_to_dipole  # [m]
+        self.field_dipole = field_dipole  # [T]
+        self.length_dipole = length_dipole  # [m]
+
+        self.length_drift_dipole_to_lens = length_drift_dipole_to_lens  # [m]
+        self.length_plasma_lens = length_plasma_lens  # [m]
+
+        self.imaging_energy = imaging_energy  # [eV]
+        self.length_drift_lens_to_screen = length_drift_lens_to_screen  # [m]
+
+        self._radius_plasma_lens = radius_plasma_lens  # [m]
+        self.disable_lens_nonlinearity = disable_lens_nonlinearity  # [bool]
+
+        super().__init__()
+
+    def get_length(self):
+        """
+        Calculates the total length of the spectrometer lattice.
+
+        Returns
+        -------
+        total_length : [m] float
+            Total length of the spectrometer lattice [m].
+        """
+        total_length = (self.length_drift_to_dipole
+                        + self.length_dipole
+                        + self.length_drift_dipole_to_lens
+                        + self.length_plasma_lens
+                        + self.length_drift_lens_to_screen)
+        return total_length
+
+    @property
+    def radius_plasma_lens(self):
+        """
+        Gets the radius of the plasma lens.
+
+        Returns
+        -------
+        self._radius_plasma_lens : [m] float
+            Radius of the plasma lens [m].
+        """
+        return self._radius_plasma_lens
+
+    @radius_plasma_lens.setter
+    def radius_plasma_lens(self, value):
+        """
+        Sets the radius of the plasma lens.
+
+        Parameters
+        ----------
+        value : [m] float
+            Radius of the plasma lens. Must be greater than or equal to zero.
+
+        Raises
+        ------
+        AssertionError
+            If the radius is less than zero.
+        """
+        assert value >= 0, 'Radius of the plasma lens must be greater than or equal to zero.'
+        self._radius_plasma_lens = value
+
+    @abstractmethod
+    def track(self, beam, savedepth=0, runnable=None, verbose=False):
+        """
+        Tracks the beam through the spectrometer.
+
+        Parameters
+        ----------
+        beam : object
+            The beam object to be tracked.
+
+        savedepth : int, optional
+            Depth of saved states during tracking. Default is 0.
+
+        runnable : callable or None, optional
+            Runnable object for tracking. Default is ``None``.
+
+        verbose : bool, optional
+            Flag to enable verbose output. Default is ``False``.
+
+        Returns
+        -------
+        object
+            Result of the tracking process.
+        """
+        return super().track(beam, savedepth, runnable, verbose)