diff --git a/optika/systems/__init__.py b/optika/systems/__init__.py
new file mode 100644
index 0000000..ce20b15
--- /dev/null
+++ b/optika/systems/__init__.py
@@ -0,0 +1,12 @@
+"""
+Optical systems consisting of multiple optical surfaces.
+"""
+
+from ._systems import AbstractSystem
+from ._sequential import AbstractSequentialSystem, SequentialSystem
+
+__all__ = [
+    "AbstractSystem",
+    "AbstractSequentialSystem",
+    "SequentialSystem",
+]
diff --git a/optika/systems.py b/optika/systems/_sequential.py
similarity index 97%
rename from optika/systems.py
rename to optika/systems/_sequential.py
index e9a63db..a39b7cc 100644
--- a/optika/systems.py
+++ b/optika/systems/_sequential.py
@@ -1,12 +1,9 @@
-"""
-Optical systems consisting of multiple optical surfaces.
-"""
-
 from __future__ import annotations
 from typing import Sequence, Callable, Any, ClassVar
 import abc
 import dataclasses
 import functools
+from ezdxf.addons.r12writer import R12FastStreamWriter
 import astropy.units as u
 import astropy.visualization
 import numpy as np
@@ -16,50 +13,14 @@
 import matplotlib.pyplot as plt
 import named_arrays as na
 import optika
-from ezdxf.addons.r12writer import R12FastStreamWriter
+from . import AbstractSystem
 
 __all__ = [
-    "AbstractSystem",
     "AbstractSequentialSystem",
     "SequentialSystem",
 ]
 
 
-@dataclasses.dataclass(eq=False, repr=False)
-class AbstractSystem(
-    optika.mixins.DxfWritable,
-    optika.mixins.Plottable,
-    optika.mixins.Printable,
-    optika.mixins.Transformable,
-    optika.mixins.Shaped,
-):
-    """
-    An interface describing an optical system.
-
-    Could potentially be sequential or non-sequential.
-    """
-
-    @abc.abstractmethod
-    def image(
-        self,
-        scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar],
-        **kwargs: Any,
-    ) -> na.SpectralPositionalVectorArray:
-        """
-        Forward model of the optical system.
-        Maps the given spectral radiance of a scene to detector counts.
-
-        Parameters
-        ----------
-        scene
-            The spectral radiance of the scene as a function of wavelength
-            and field position.
-        kwargs
-            Additional keyword arguments used by subclass implementations
-            of this method.
-        """
-
-
 @dataclasses.dataclass(eq=False, repr=False)
 class AbstractSequentialSystem(
     AbstractSystem,
diff --git a/optika/_tests/test_systems.py b/optika/systems/_sequential_test.py
similarity index 83%
rename from optika/_tests/test_systems.py
rename to optika/systems/_sequential_test.py
index 5f3f985..ee2f3bf 100644
--- a/optika/_tests/test_systems.py
+++ b/optika/systems/_sequential_test.py
@@ -4,54 +4,8 @@
 import astropy.units as u
 import named_arrays as na
 import optika
-from . import test_mixins
-
-
-class AbstractTestAbstractSystem(
-    test_mixins.AbstractTestPlottable,
-    test_mixins.AbstractTestPrintable,
-    test_mixins.AbstractTestTransformable,
-    test_mixins.AbstractTestShaped,
-):
-
-    @pytest.mark.parametrize(
-        argnames="scene",
-        argvalues=[
-            na.FunctionArray(
-                inputs=na.SpectralPositionalVectorArray(
-                    wavelength=na.linspace(
-                        start=530 * u.nm,
-                        stop=531 * u.nm,
-                        axis="wavelength",
-                        num=2,
-                    ),
-                    position=na.Cartesian2dVectorLinearSpace(
-                        start=-1,
-                        stop=+1,
-                        axis=na.Cartesian2dVectorArray("field_x", "field_y"),
-                        num=11,
-                    ),
-                ),
-                outputs=na.random.uniform(
-                    low=0 * u.photon / u.cm**2 / u.arcsec**2 / u.s / u.nm,
-                    high=100 * u.photon / u.cm**2 / u.arcsec**2 / u.s / u.nm,
-                    shape_random=dict(field_x=10, field_y=10),
-                ),
-            )
-        ],
-    )
-    def test_image(
-        self,
-        a: optika.systems.AbstractSystem,
-        scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar],
-    ):
-        result = a.image(scene)
-        assert isinstance(result, na.FunctionArray)
-        assert isinstance(result.inputs, na.SpectralPositionalVectorArray)
-        assert isinstance(result.outputs, na.AbstractScalar)
-        assert np.all(result.inputs.wavelength > 0 * u.nm)
-        assert na.unit_normalized(result.inputs.position).is_equivalent(u.mm)
-        assert result.outputs.sum() != (0 * u.electron)
+from .._tests import test_mixins
+from ._systems_test import AbstractTestAbstractSystem
 
 
 class AbstractTestAbstractSequentialSystem(
diff --git a/optika/systems/_systems.py b/optika/systems/_systems.py
new file mode 100644
index 0000000..5a779ea
--- /dev/null
+++ b/optika/systems/_systems.py
@@ -0,0 +1,45 @@
+from __future__ import annotations
+from typing import Any
+import abc
+import dataclasses
+import named_arrays as na
+import optika
+
+__all__ = [
+    "AbstractSystem",
+]
+
+
+@dataclasses.dataclass(eq=False, repr=False)
+class AbstractSystem(
+    optika.mixins.DxfWritable,
+    optika.mixins.Plottable,
+    optika.mixins.Printable,
+    optika.mixins.Transformable,
+    optika.mixins.Shaped,
+):
+    """
+    An interface describing an optical system.
+
+    Could potentially be sequential or non-sequential.
+    """
+
+    @abc.abstractmethod
+    def image(
+        self,
+        scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar],
+        **kwargs: Any,
+    ) -> na.SpectralPositionalVectorArray:
+        """
+        Forward model of the optical system.
+        Maps the given spectral radiance of a scene to detector counts.
+
+        Parameters
+        ----------
+        scene
+            The spectral radiance of the scene as a function of wavelength
+            and field position.
+        kwargs
+            Additional keyword arguments used by subclass implementations
+            of this method.
+        """
diff --git a/optika/systems/_systems_test.py b/optika/systems/_systems_test.py
new file mode 100644
index 0000000..4f42d24
--- /dev/null
+++ b/optika/systems/_systems_test.py
@@ -0,0 +1,53 @@
+import pytest
+import numpy as np
+import astropy.units as u
+import named_arrays as na
+import optika
+from .._tests import test_mixins
+
+
+class AbstractTestAbstractSystem(
+    test_mixins.AbstractTestPlottable,
+    test_mixins.AbstractTestPrintable,
+    test_mixins.AbstractTestTransformable,
+    test_mixins.AbstractTestShaped,
+):
+
+    @pytest.mark.parametrize(
+        argnames="scene",
+        argvalues=[
+            na.FunctionArray(
+                inputs=na.SpectralPositionalVectorArray(
+                    wavelength=na.linspace(
+                        start=530 * u.nm,
+                        stop=531 * u.nm,
+                        axis="wavelength",
+                        num=2,
+                    ),
+                    position=na.Cartesian2dVectorLinearSpace(
+                        start=-1,
+                        stop=+1,
+                        axis=na.Cartesian2dVectorArray("field_x", "field_y"),
+                        num=11,
+                    ),
+                ),
+                outputs=na.random.uniform(
+                    low=0 * u.photon / u.cm**2 / u.arcsec**2 / u.s / u.nm,
+                    high=100 * u.photon / u.cm**2 / u.arcsec**2 / u.s / u.nm,
+                    shape_random=dict(field_x=10, field_y=10),
+                ),
+            )
+        ],
+    )
+    def test_image(
+        self,
+        a: optika.systems.AbstractSystem,
+        scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar],
+    ):
+        result = a.image(scene)
+        assert isinstance(result, na.FunctionArray)
+        assert isinstance(result.inputs, na.SpectralPositionalVectorArray)
+        assert isinstance(result.outputs, na.AbstractScalar)
+        assert np.all(result.inputs.wavelength > 0 * u.nm)
+        assert na.unit_normalized(result.inputs.position).is_equivalent(u.mm)
+        assert result.outputs.sum() != (0 * u.electron)