diff --git a/README.md b/README.md
index 02a5c75..64d9a61 100644
--- a/README.md
+++ b/README.md
@@ -1,25 +1,30 @@
-# EyeSimulator
+# Visisipy: accessible vision simulations in Python
 
-This project aims to provide a standardized method to perform optical (ray tracing) simulations on eye models. `EyeSimulator` is just a working name and should preferrably be replaced with something better before releasing this library.
+Visisipy (pronounced `/ˌvɪsəˈsɪpi/`, like Mississippi but with a V) is a Python library for optical simulations of the eye.
+It provides an easy-to-use interface to define and build eye models, and to perform common ophthalmic analyses on these models.
 
 ## Goals
 
 1. Provide a uniform interface to define, build and analyze various types of eye models, using abstractions that make sense in a clinical context;
-2. Provide an accessible interface to the most common analyses on these models;
-3. Decouple the model definition interface as much as possible from the simulation software, i.e. Zemax OpticStudio;
-4. Introduce this as a standardized method for optical simulations of the eye, that can be used by the broader physiological optics community.
+2. Provide a collection of ready-to-use eye models, such as the Navarro model[^navarro], that can be customized at need;
+3. Provide an accessible interface to clinically relevant analyses with these models.
+
+All calculations are currently performed in OpticStudio through the [ZOSPy][zospy] library[^zospy], but visisipy is designed in a modular fashion to allow for other backends in the future.
 
 ## Contributing
 
-Visisipy aims to be a community-driven project and warmly accepts contributions. If you want to contribute, please email us (visisipy@mreye.nl) or [open a new discussion](https://github.com/MREYE-LUMC/visisipy/discussions).
+Visisipy aims to be a community-driven project and warmly accepts contributions.
+If you want to contribute, please email us (visisipy@mreye.nl) or [open a new discussion](https://github.com/MREYE-LUMC/visisipy/discussions).
+
+## Installation
 
-## Current implementation
+Visisipy can be installed through `pip`:
 
-In its current state, `EyeSimulator` provides
+```bash
+pip install git+https://github.com/MREYE-LUMC/visisipy.git
+```
 
-- Classes to define models in terms of clinically relevant and measurable parameters (`EyeGeometry`);
-- An abstract base class `BaseEye` and a simple implementation `Eye`, which allows to build and analyze these eye models;
-- A `Surface` class, which acts as a bridge between the abstractions used in this library and Zemax OpticStudio. This allows for seamless integration of OpticStudio models in our `Eye` class.
+Visisipy will be made available through PyPI and Conda as soon as possible.
 
 ## Example
 
@@ -35,7 +40,11 @@ model = visisipy.EyeModel()
 model.build()
 
 # Perform a raytrace analysis
-raytrace = visisipy.analysis.raytrace(coordinates=zip([0] * 5, range(0, 60, 10)))
+coordinates = [(0, 0), (0, 10), (0, 20), (0, 30), (0, 40)]
+raytrace = visisipy.analysis.raytrace(coordinates=coordinates)
+
+# Alternatively, the model can be built and analyzed in one go:
+# raytrace = visisipy.analysis.raytrace(model, coordinates=zip([0] * 5, range(0, 60, 10)))
 
 # Visualize the model
 fig, ax = plt.subplots()
@@ -49,42 +58,91 @@ sns.lineplot(raytrace, x="z", y="y", hue="field", ax=ax)
 plt.show()
 ```
 
+### Configure the backend
+
+Visisipy uses OpticStudio as a backend for calculations; this is currently the only supported backend.
+This backend is automatically started and managed in the background, but can also be configured manually.
+
 ```python
-import zospy as zp
-from visisipy import EyeModel, NavarroGeometry
-from visisipy.opticstudio import OpticStudioEye
+import visisipy
+
+# Use OpticStudio in standalone mode (default)
+visisipy.set_backend("opticstudio")
+
+# Use OpticStudio in extension mode
+visisipy.set_backend("opticstudio", mode="extension")
+
+# Use OpticStudio in extension mode with ray aiming enabled
+visisipy.set_backend("opticstudio", mode="extension", ray_aiming="real")
+
+# Get the OpticStudioSystem from visisipy to interact with it manually
+# This only works when the backend is set to "opticstudio"
+# See https://zospy.readthedocs.io/en/latest/api/zospy.zpcore.OpticStudioSystem.html for documentation of this object
+oss = visisipy.backend.get_oss()
+```
 
-# Initialize ZOSPy
-zos = zp.ZOS()
-oss = zos.connect()
+### Create a custom eye model from clinical parameters
 
-# Initialize an eye, using a slightly modified Navarro model
-eye_model = EyeModel(NavarroGeometry(iris_radius=0.5))
-eye = OpticStudioEye(eye_model)
+An eye model in visispy consists of two parts: the geometry and the material properties.
+The geometry is defined by `visisipy.models.EyeGeometry`, and the material properties are defined by `visisipy.models.Materials`.
+They are combined in `visisipy.EyeModel` to constitute a complete eye model.
 
-# Update a parameter of one of the eye's surfaces
-eye.lens_front.refractive_index += 1.0
+```python
+import visisipy
+
+geometry = visisipy.models.create_geometry(
+    axial_length=20,
+    cornea_thickness=0.5,
+    anterior_chamber_depth=3,
+    lens_thickness=4,
+    cornea_front_radius=7,
+    cornea_front_asphericity=0,
+    cornea_back_radius=6,
+    cornea_back_asphericity=0,
+    lens_front_radius=10,
+    lens_front_asphericity=0,
+    lens_back_radius=-6,
+    lens_back_asphericity=0,
+    retina_radius=-12,
+    retina_asphericity=0,
+    pupil_radius=1.0,
+)
+
+# Use this geometry together with the refractive indices of the Navarro model
+model = visisipy.EyeModel(geometry=geometry, materials=visisipy.models.materials.NavarroMaterials())
+
+# NavarroMaterials is the default, so this is equivalent:
+model = visisipy.EyeModel(geometry=geometry)
+```
 
-# Build the eye in OpticStudio
-eye.build(oss)
+### Interact with the eye model in OpticStudio
 
-# Change the lens's refractive index back
-eye.lens_front.refractive_index -= 1.0
+```python
+import visisipy
 
-# Insert a new (dummy) surface in OpticStudio
-oss.LDE.InsertNewSurfaceAt(1).Comment = "input beam"
+# Just use the default Navarro model
+model = visisipy.EyeModel()
 
-# Relink the surfaces of the eye, because a new surface was added
-eye.relink_surfaces(oss)
+# Build the model in OpticStudio
+built_model: visisipy.opticstudio.OpticStudioEye = model.build()
 
-# Check if the relinking succeeded
-assert eye.lens_back.radius == eye_model.geometry.lens_back_curvature
+# Update the lens front radius
+built_model.lens_front.radius = 10.5
 ```
 
+## Planned functions
+
+- Generation of realistic randomized eye models using the method proposed by Rozema et al.[^rozema]
+
 ## Future ideas
 
-- Provide (customizable) geometry definitions for various standard eye models, e.g. `NavarroGeometry`, `GullstrandGeometry`.
-    - Most likely as child classes of `EyeGeometry`. The default values currently used in `EyeGeometry` will then be removed.
-    - The same applies for the `EyeMaterials` class.
-- Add more implementations of `BaseEye`, e.g. `ReverseEye`, a reversed version of `Eye`.
-- Integrate the eye plot functions defined in [chaasjes/utilities](https://git.lumc.nl/chaasjes/utilities/-/tree/main/utilities/plots/eye) in this library for easy visualization of eye models.
\ No newline at end of file
+- Provide (customizable) geometry definitions for other standard eye models, e.g. `GullstrandGeometry`.
+- Add support for reversed eyes.
+- Add support for other (open source) ray tracing backends.
+
+[zospy]: https://zospy.readthedocs.io/
+
+[//]: # (References)
+[^navarro]: Escudero-Sanz, I., & Navarro, R. (1999). Off-axis aberrations of a wide-angle schematic eye model. JOSA A, 16(8), 1881–1891. https://doi.org/10.1364/JOSAA.16.001881
+[^rozema]: Rozema, J. J., Rodriguez, P., Navarro, R., & Tassignon, M.-J. (2016). SyntEyes: A Higher-Order Statistical Eye Model for Healthy Eyes. Investigative Ophthalmology & Visual Science, 57(2), 683–691. https://doi.org/10.1167/iovs.15-18067
+[^zospy]: Vught, L. van, Haasjes, C., & Beenakker, J.-W. M. (2024). ZOSPy: Optical ray tracing in Python through OpticStudio. Journal of Open Source Software, 9(96), 5756. https://doi.org/10.21105/joss.05756
diff --git a/examples/Patient-specific mapping of fundus photographs to three-dimensional ocular imaging/helpers.py b/examples/Patient-specific mapping of fundus photographs to three-dimensional ocular imaging/helpers.py
index f21cac5..3c87054 100644
--- a/examples/Patient-specific mapping of fundus photographs to three-dimensional ocular imaging/helpers.py	
+++ b/examples/Patient-specific mapping of fundus photographs to three-dimensional ocular imaging/helpers.py	
@@ -1,15 +1,17 @@
 from __future__ import annotations
 
-from typing import NamedTuple
+from typing import TYPE_CHECKING, NamedTuple
 
 import numpy as np
-import pandas as pd
+
+if TYPE_CHECKING:
+    import pandas as pd
 
 
 def get_ray_output_angle(
-        df: pd.DataFrame,
-        reference_point: tuple[float, float] = (0, 0),
-        coordinate="y",
+    df: pd.DataFrame,
+    reference_point: tuple[float, float] = (0, 0),
+    coordinate="y",
 ):
     """Calculate the output angle of a ray with respect to the optical axis and a reference point."""
     x0, y0 = reference_point
@@ -30,22 +32,18 @@ class InputOutputAngles(NamedTuple):
 
     @classmethod
     def from_ray_trace_result(
-            cls,
-            raytrace_result: pd.DataFrame,
-            np2: float,
-            np2_navarro: float = None,
-            retina_center: float = None,
-            patient: int = None,
-            coordinate="y",
-    ) -> "InputOutputAngles":
+        cls,
+        raytrace_result: pd.DataFrame,
+        np2: float,
+        np2_navarro: float | None = None,
+        retina_center: float | None = None,
+        patient: int | None = None,
+        coordinate="y",
+    ) -> InputOutputAngles:
         return cls(
             input_angle_field=raytrace_result.field[0][1],
-            output_angle_pupil=get_ray_output_angle(
-                raytrace_result, reference_point=(0, 0), coordinate=coordinate
-            ),
-            output_angle_np2=get_ray_output_angle(
-                raytrace_result, reference_point=(np2, 0), coordinate=coordinate
-            ),
+            output_angle_pupil=get_ray_output_angle(raytrace_result, reference_point=(0, 0), coordinate=coordinate),
+            output_angle_np2=get_ray_output_angle(raytrace_result, reference_point=(np2, 0), coordinate=coordinate),
             output_angle_retina_center=(
                 get_ray_output_angle(
                     raytrace_result,
diff --git a/pyproject.toml b/pyproject.toml
index 5d7d957..4e7d48c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -57,6 +57,7 @@ extend-ignore = [
 ]
 
 [tool.ruff.lint.extend-per-file-ignores]
+"examples/*" = [ "INP001" ]
 "tests/*" = [ "ARG001", "ARG002" ]
 "tests/opticstudio/*" = [ "SLF001" ]
 
diff --git a/tests/opticstudio/test_analysis.py b/tests/opticstudio/test_analysis.py
index e16fe5b..8bec3cf 100644
--- a/tests/opticstudio/test_analysis.py
+++ b/tests/opticstudio/test_analysis.py
@@ -13,15 +13,32 @@ def opticstudio_analysis(opticstudio_backend):
 
 
 class TestCardinalPointsAnalysis:
-    @pytest.mark.parametrize("surface_1,surface_2,expectation", [
-        (None, None, does_not_raise()),
-        (1, 6, does_not_raise()),
-        (2, 4, does_not_raise()),
-        (-1, 7, pytest.raises(ValueError, match="surface_1 and surface_2 must be between 1 and the number of surfaces in the system")),
-        (1, 8, pytest.raises(ValueError, match="surface_1 and surface_2 must be between 1 and the number of surfaces in the system")),
-        (3, 3, pytest.raises(ValueError, match="surface_1 must be less than surface_2")),
-        (4, 2, pytest.raises(ValueError, match="surface_1 must be less than surface_2")),
- ])
+    @pytest.mark.parametrize(
+        "surface_1,surface_2,expectation",
+        [
+            (None, None, does_not_raise()),
+            (1, 6, does_not_raise()),
+            (2, 4, does_not_raise()),
+            (
+                -1,
+                7,
+                pytest.raises(
+                    ValueError,
+                    match="surface_1 and surface_2 must be between 1 and the number of surfaces in the system",
+                ),
+            ),
+            (
+                1,
+                8,
+                pytest.raises(
+                    ValueError,
+                    match="surface_1 and surface_2 must be between 1 and the number of surfaces in the system",
+                ),
+            ),
+            (3, 3, pytest.raises(ValueError, match="surface_1 must be less than surface_2")),
+            (4, 2, pytest.raises(ValueError, match="surface_1 must be less than surface_2")),
+        ],
+    )
     def test_cardinal_points(self, opticstudio_backend, opticstudio_analysis, surface_1, surface_2, expectation):
         opticstudio_backend.build_model(EyeModel())
 
@@ -96,12 +113,12 @@ class TestRefractionAnalysis:
         ],
     )
     def test_refraction(
-            self,
-            opticstudio_analysis,
-            use_higher_order_aberrations,
-            wavelength,
-            pupil_diameter,
-            monkeypatch,
+        self,
+        opticstudio_analysis,
+        use_higher_order_aberrations,
+        wavelength,
+        pupil_diameter,
+        monkeypatch,
     ):
         monkeypatch.setattr(opticstudio_analysis._backend, "model", MockOpticstudioModel())
 
diff --git a/tests/opticstudio/test_surfaces.py b/tests/opticstudio/test_surfaces.py
index 458a90d..68dc09d 100644
--- a/tests/opticstudio/test_surfaces.py
+++ b/tests/opticstudio/test_surfaces.py
@@ -254,6 +254,7 @@ class MockOpticStudioZernikeSurface(BaseOpticStudioZernikeSurface):
 
         This class is necessary because BaseOpticStudioZernikeSurface cannot be instantiated directly.
         """
+
         _TYPE = "ZernikeStandardSag"
 
         def __init__(self, *args, **kwargs):
@@ -281,13 +282,19 @@ def test_init_negative_zernike_coefficient_raises_valueerror(self, key):
             )
 
     @pytest.mark.parametrize(
-        "n,value,maximum_term,expectation", [
+        "n,value,maximum_term,expectation",
+        [
             (0, 1.234, 1, pytest.raises(ValueError, match="Zernike coefficient must be larger than 0")),
             (1, 1.234, 2, does_not_raise()),
             (2, 1.234, 2, does_not_raise()),
-            (2, 1.234, 1,
-             pytest.raises(ValueError, match="Zernike coefficient must be smaller than the maximum term 1")),
-        ])
+            (
+                2,
+                1.234,
+                1,
+                pytest.raises(ValueError, match="Zernike coefficient must be smaller than the maximum term 1"),
+            ),
+        ],
+    )
     def test_set_zernike_coefficient(self, new_oss, n: int, value: float, maximum_term: int, expectation):
         surface = self.MockOpticStudioZernikeSurface(
             comment="Test",
@@ -301,13 +308,16 @@ def test_set_zernike_coefficient(self, new_oss, n: int, value: float, maximum_te
             assert surface.surface.SurfaceData.GetNthZernikeCoefficient(n) == value
 
     @pytest.mark.parametrize(
-        "n,coefficients,expectation", [
+        "n,coefficients,expectation",
+        [
             (2, {1: 1.234, 2: 3.456}, does_not_raise()),
-            (4, {1: 1.234, 2: 3.456}, pytest.raises(
-                ValueError, match="Zernike coefficient must be smaller than the "
-                                  "maximum term 3")),
+            (
+                4,
+                {1: 1.234, 2: 3.456},
+                pytest.raises(ValueError, match="Zernike coefficient must be smaller than the " "maximum term 3"),
+            ),
             (0, {1: 1.234, 2: 3.456}, pytest.raises(ValueError, match="Zernike coefficient must be larger than 0")),
-        ]
+        ],
     )
     def test_get_zernike_coefficient(self, new_oss, n, coefficients, expectation):
         surface = OpticStudioZernikeStandardSagSurface(
@@ -324,13 +334,15 @@ def test_get_zernike_coefficient(self, new_oss, n, coefficients, expectation):
 
 class TestOpticStudioZernikeStandardSagSurface:
     @pytest.mark.parametrize(
-        "maximum_term,coefficients,extrapolate,decenter_x,decenter_y", [
+        "maximum_term,coefficients,extrapolate,decenter_x,decenter_y",
+        [
             (3, {1: 1.0, 2: 2.0, 3: 3.0}, 1, 0.0, 0.0),
             (3, {1: 1.0, 2: 2.0, 3: 3.0}, 0, -1.0, 1.0),
             (3, {1: 1.0, 2: 2.0, 3: 3.0}, 1, 1.0, -1.0),
             (231, {}, 0, 0.0, 0.0),
             (231, None, 1, 0.0, 0.0),
-        ])
+        ],
+    )
     def test_build(self, new_oss, maximum_term, coefficients, extrapolate, decenter_x, decenter_y):
         surface = OpticStudioZernikeStandardSagSurface(
             comment="Test comment",
@@ -397,13 +409,15 @@ def test_extrapolate(self, new_oss, extrapolate):
 
 class TestOpticStudioZernikeStandardPhaseSurface:
     @pytest.mark.parametrize(
-        "maximum_term,coefficients,extrapolate,diffract_order", [
+        "maximum_term,coefficients,extrapolate,diffract_order",
+        [
             (3, {1: 1.0, 2: 2.0, 3: 3.0}, 1, 1),
             (3, {1: 1.0, 2: 2.0, 3: 3.0}, 0, 2.3),
             (3, {1: 1.0, 2: 2.0, 3: 3.0}, 1, 3.4),
             (231, {}, 0, 4.5),
             (231, None, 1, 5.6),
-        ])
+        ],
+    )
     def test_build(self, new_oss, maximum_term, coefficients, extrapolate, diffract_order):
         surface = OpticStudioZernikeStandardPhaseSurface(
             comment="Test comment",
@@ -467,11 +481,11 @@ def test_make_surface(self):
         [
             (1, 2, 3, 4, "BK7"),
             (
-                    5,
-                    6,
-                    7,
-                    8,
-                    MaterialModel(refractive_index=1.5, abbe_number=50, partial_dispersion=0.67),
+                5,
+                6,
+                7,
+                8,
+                MaterialModel(refractive_index=1.5, abbe_number=50, partial_dispersion=0.67),
             ),
         ],
     )
@@ -496,9 +510,9 @@ def test_make_standard_surface(self, radius, thickness, semi_diameter, aspherici
         [
             (1, 2, "BK7"),
             (
-                    3,
-                    4,
-                    MaterialModel(refractive_index=1.5, abbe_number=50, partial_dispersion=0.67),
+                3,
+                4,
+                MaterialModel(refractive_index=1.5, abbe_number=50, partial_dispersion=0.67),
             ),
         ],
     )
diff --git a/tests/test_wavefront.py b/tests/test_wavefront.py
index 2db1f00..dbfd592 100644
--- a/tests/test_wavefront.py
+++ b/tests/test_wavefront.py
@@ -6,28 +6,34 @@
 
 
 class TestZernikeCoefficients:
-    @pytest.mark.parametrize("terms,result,expectation", [
-        (None, {}, does_not_raise()),
-        ({}, {}, does_not_raise()),
-        ({1: 0.5, 2: 0.9}, {1: 0.5, 2: 0.9}, does_not_raise()),
-        ({1: 0.5, "x": 0.9}, {}, pytest.raises(TypeError, match="All keys must be integers")),
-        ({1: 0.5, 2.8: 0.9}, {}, pytest.raises(TypeError, match="All keys must be integers")),
-        ({1: 0.5, 0: 0.9}, {}, pytest.raises(ValueError, match="The Zernike coefficients must be larger than 0")),
-        ({-1: 0.5, 2: 0.9}, {}, pytest.raises(ValueError, match="The Zernike coefficients must be larger than 0")),
-    ])
+    @pytest.mark.parametrize(
+        "terms,result,expectation",
+        [
+            (None, {}, does_not_raise()),
+            ({}, {}, does_not_raise()),
+            ({1: 0.5, 2: 0.9}, {1: 0.5, 2: 0.9}, does_not_raise()),
+            ({1: 0.5, "x": 0.9}, {}, pytest.raises(TypeError, match="All keys must be integers")),
+            ({1: 0.5, 2.8: 0.9}, {}, pytest.raises(TypeError, match="All keys must be integers")),
+            ({1: 0.5, 0: 0.9}, {}, pytest.raises(ValueError, match="The Zernike coefficients must be larger than 0")),
+            ({-1: 0.5, 2: 0.9}, {}, pytest.raises(ValueError, match="The Zernike coefficients must be larger than 0")),
+        ],
+    )
     def test_init(self, terms, result, expectation):
         with expectation:
             zernike = ZernikeCoefficients(terms)
             assert zernike == result
 
-    @pytest.mark.parametrize("n,value,expectation", [
-        (1, 1.2, does_not_raise()),
-        (1, 0.0, does_not_raise()),
-        (0, 0.0, pytest.raises(ValueError, match="The coefficient must be larger than 0")),
-        (-1, 0.0, pytest.raises(ValueError, match="The coefficient must be larger than 0")),
-        ("x", 0.0, pytest.raises(TypeError, match="The key must be an integer")),
-        (1.5, 0.0, pytest.raises(TypeError, match="The key must be an integer")),
-    ])
+    @pytest.mark.parametrize(
+        "n,value,expectation",
+        [
+            (1, 1.2, does_not_raise()),
+            (1, 0.0, does_not_raise()),
+            (0, 0.0, pytest.raises(ValueError, match="The coefficient must be larger than 0")),
+            (-1, 0.0, pytest.raises(ValueError, match="The coefficient must be larger than 0")),
+            ("x", 0.0, pytest.raises(TypeError, match="The key must be an integer")),
+            (1.5, 0.0, pytest.raises(TypeError, match="The key must be an integer")),
+        ],
+    )
     def test_setitem(self, n: int, value: float, expectation):
         zernike = ZernikeCoefficients()
 
@@ -35,15 +41,18 @@ def test_setitem(self, n: int, value: float, expectation):
             zernike[n] = value
             assert zernike[n] == value
 
-    @pytest.mark.parametrize("n,expected_value,expectation", [
-        (1, 0.5, does_not_raise()),
-        (2, 0.9, does_not_raise()),
-        (3, 0.0, does_not_raise()),
-        (0, 0.0, pytest.raises(ValueError, match="The coefficient must be larger than 0")),
-        (-1, 0.0, pytest.raises(ValueError, match="The coefficient must be larger than 0")),
-        ("x", 0.0, pytest.raises(TypeError, match="The key must be an integer")),
-        (1.5, 0.0, pytest.raises(TypeError, match="The key must be an integer")),
-    ])
+    @pytest.mark.parametrize(
+        "n,expected_value,expectation",
+        [
+            (1, 0.5, does_not_raise()),
+            (2, 0.9, does_not_raise()),
+            (3, 0.0, does_not_raise()),
+            (0, 0.0, pytest.raises(ValueError, match="The coefficient must be larger than 0")),
+            (-1, 0.0, pytest.raises(ValueError, match="The coefficient must be larger than 0")),
+            ("x", 0.0, pytest.raises(TypeError, match="The key must be an integer")),
+            (1.5, 0.0, pytest.raises(TypeError, match="The key must be an integer")),
+        ],
+    )
     def test_getitem(self, n, expected_value, expectation):
         zernike = ZernikeCoefficients({1: 0.5, 2: 0.9})
 
diff --git a/visisipy/analysis/cardinal_points.py b/visisipy/analysis/cardinal_points.py
index c38e359..3fa1332 100644
--- a/visisipy/analysis/cardinal_points.py
+++ b/visisipy/analysis/cardinal_points.py
@@ -7,7 +7,6 @@
 from visisipy.backend import get_backend
 
 if TYPE_CHECKING:
-    from pandas import DataFrame
     from visisipy import EyeModel
 
 
@@ -18,6 +17,7 @@ class CardinalPoints(NamedTuple):
     """
     The cardinal points of a system in object and image space.
     """
+
     object: float
     image: float
 
@@ -42,6 +42,7 @@ class CardinalPointsResult:
     anti_nodal_points : CardinalPoints
         The anti-nodal points of the system.
     """
+
     focal_lengths: CardinalPoints
     focal_points: CardinalPoints
     principal_points: CardinalPoints
diff --git a/visisipy/analysis/refraction.py b/visisipy/analysis/refraction.py
index 2303944..6e32d92 100644
--- a/visisipy/analysis/refraction.py
+++ b/visisipy/analysis/refraction.py
@@ -15,14 +15,14 @@
 
 @analysis
 def refraction(
-        model: EyeModel | None,  # noqa: ARG001
-        field_coordinate: tuple[float, float] | None = None,
-        wavelength: float | None = None,
-        pupil_diameter: float | None = None,
-        field_type: Literal["angle", "object_height"] = "angle",
-        *,
-        use_higher_order_aberrations: bool = True,
-        return_raw_result: bool = False,  # noqa: ARG001
+    model: EyeModel | None,  # noqa: ARG001
+    field_coordinate: tuple[float, float] | None = None,
+    wavelength: float | None = None,
+    pupil_diameter: float | None = None,
+    field_type: Literal["angle", "object_height"] = "angle",
+    *,
+    use_higher_order_aberrations: bool = True,
+    return_raw_result: bool = False,  # noqa: ARG001
 ) -> FourierPowerVectorRefraction | tuple[FourierPowerVectorRefraction, Any]:
     """Calculates the ocular refraction.
 
diff --git a/visisipy/backend.py b/visisipy/backend.py
index 958a9b2..2b35aad 100644
--- a/visisipy/backend.py
+++ b/visisipy/backend.py
@@ -124,6 +124,6 @@ def get_oss() -> OpticStudioSystem | None:
     os_backend = importlib.import_module("visisipy.opticstudio.backend")
 
     if _BACKEND is os_backend.OpicStudioBackend:
-        return os_backend._oss
+        return os_backend._oss  # noqa: SLF001
 
     return None
diff --git a/visisipy/models/geometry.py b/visisipy/models/geometry.py
index d157acb..250a6a8 100644
--- a/visisipy/models/geometry.py
+++ b/visisipy/models/geometry.py
@@ -152,6 +152,7 @@ class ZernikeStandardSagSurface(StandardSurface):
     ValueError
         If the Zernike coefficients contain terms that are greater than the maximum term.
     """
+
     zernike_coefficients: ZernikeCoefficients | dict[int, float] = field(default_factory=dict)
     extrapolate: bool = True
     zernike_decenter_x: float = 0
@@ -191,6 +192,7 @@ class ZernikeStandardPhaseSurface(StandardSurface):
     ValueError
         If the Zernike coefficients contain terms that are greater than the maximum term.
     """
+
     zernike_coefficients: ZernikeCoefficients | dict[int, float] = field(default_factory=dict)
     extrapolate: bool = True
     diffraction_order: float = 1
@@ -342,10 +344,10 @@ def retina(self, value: StandardSurface) -> None:
     def axial_length(self) -> float:
         """Axial length of the eye, in mm."""
         return (
-                self.cornea_front.thickness
-                + self.cornea_back.thickness
-                + self.lens_front.thickness
-                + self.lens_back.thickness
+            self.cornea_front.thickness
+            + self.cornea_back.thickness
+            + self.lens_front.thickness
+            + self.lens_back.thickness
         )
 
     @property
@@ -434,11 +436,11 @@ def _update_attribute_if_specified(obj: Surface, attribute: str, value: Any):
 
 
 def _calculate_vitreous_thickness(
-        geometry: EyeGeometry,
-        axial_length: float | None = None,
-        cornea_thickness: float | None = None,
-        anterior_chamber_depth: float | None = None,
-        lens_thickness: float | None = None,
+    geometry: EyeGeometry,
+    axial_length: float | None = None,
+    cornea_thickness: float | None = None,
+    anterior_chamber_depth: float | None = None,
+    lens_thickness: float | None = None,
 ) -> float:
     """Calculate the thickness of the vitreous body for a partially initialized eye geometry."""
     _axial_length = geometry.axial_length if axial_length is None else axial_length
@@ -449,10 +451,10 @@ def _calculate_vitreous_thickness(
     _lens_thickness = geometry.lens_thickness if lens_thickness is None else lens_thickness
 
     if None in (
-            _axial_length,
-            _cornea_thickness,
-            _anterior_chamber_depth,
-            _lens_thickness,
+        _axial_length,
+        _cornea_thickness,
+        _anterior_chamber_depth,
+        _lens_thickness,
     ):
         raise ValueError("Cannot calculate vitreous thickness from the supplied parameters.")
 
@@ -463,24 +465,24 @@ def _calculate_vitreous_thickness(
 
 
 def create_geometry(
-        base: type[GeometryType] = NavarroGeometry,
-        axial_length: float | None = None,
-        cornea_thickness: float | None = None,
-        cornea_front_radius: float | None = None,
-        cornea_front_asphericity: float | None = None,
-        cornea_back_radius: float | None = None,
-        cornea_back_asphericity: float | None = None,
-        anterior_chamber_depth: float | None = None,
-        pupil_radius: float | None = None,
-        lens_thickness: float | None = None,
-        lens_back_radius: float | None = None,
-        lens_back_asphericity: float | None = None,
-        lens_front_radius: float | None = None,
-        lens_front_asphericity: float | None = None,
-        retina_radius: float | None = None,
-        retina_asphericity: float | None = None,
-        *,
-        estimate_cornea_back: bool = False,
+    base: type[GeometryType] = NavarroGeometry,
+    axial_length: float | None = None,
+    cornea_thickness: float | None = None,
+    cornea_front_radius: float | None = None,
+    cornea_front_asphericity: float | None = None,
+    cornea_back_radius: float | None = None,
+    cornea_back_asphericity: float | None = None,
+    anterior_chamber_depth: float | None = None,
+    pupil_radius: float | None = None,
+    lens_thickness: float | None = None,
+    lens_back_radius: float | None = None,
+    lens_back_asphericity: float | None = None,
+    lens_front_radius: float | None = None,
+    lens_front_asphericity: float | None = None,
+    retina_radius: float | None = None,
+    retina_asphericity: float | None = None,
+    *,
+    estimate_cornea_back: bool = False,
 ) -> GeometryType:
     """Create a geometry instance from clinically used parameters.
 
diff --git a/visisipy/opticstudio/analysis.py b/visisipy/opticstudio/analysis.py
index d805ae6..7f78b5d 100644
--- a/visisipy/opticstudio/analysis.py
+++ b/visisipy/opticstudio/analysis.py
@@ -1,6 +1,6 @@
 from __future__ import annotations
 
-from typing import Literal, TYPE_CHECKING
+from typing import TYPE_CHECKING, Literal
 
 import numpy as np
 import pandas as pd
@@ -81,11 +81,11 @@ def _get_zernike_coefficient(zernike_result: zp.analyses.base.AttrDict, coeffici
 
 
 def _zernike_data_to_refraction(
-        zernike_data: zp.analyses.base.AttrDict,
-        pupil_data: zp.functions.lde.PupilData,
-        wavelength: float,
-        *,
-        use_higher_order_aberrations: bool = True,
+    zernike_data: zp.analyses.base.AttrDict,
+    pupil_data: zp.functions.lde.PupilData,
+    wavelength: float,
+    *,
+    use_higher_order_aberrations: bool = True,
 ) -> FourierPowerVectorRefraction:
     z4 = _get_zernike_coefficient(zernike_data, 4) * wavelength * 4 * np.sqrt(3)
     z11 = _get_zernike_coefficient(zernike_data, 11) * wavelength * 12 * np.sqrt(5)
@@ -106,15 +106,15 @@ def _zernike_data_to_refraction(
 
     if use_higher_order_aberrations:
         return FourierPowerVectorRefraction(
-            M=(-z4 + z11 - z22 + z37) / (exit_pupil_radius ** 2),
-            J0=(-z6 + z12 - z24 + z38) / (exit_pupil_radius ** 2),
-            J45=(-z5 + z13 - z23 + z39) / (exit_pupil_radius ** 2),
+            M=(-z4 + z11 - z22 + z37) / (exit_pupil_radius**2),
+            J0=(-z6 + z12 - z24 + z38) / (exit_pupil_radius**2),
+            J45=(-z5 + z13 - z23 + z39) / (exit_pupil_radius**2),
         )
 
     return FourierPowerVectorRefraction(
-        M=(-z4) / (exit_pupil_radius ** 2),
-        J0=(-z6) / (exit_pupil_radius ** 2),
-        J45=(-z5) / (exit_pupil_radius ** 2),
+        M=(-z4) / (exit_pupil_radius**2),
+        J0=(-z6) / (exit_pupil_radius**2),
+        J45=(-z5) / (exit_pupil_radius**2),
     )
 
 
@@ -127,8 +127,9 @@ def __init__(self, backend: OpticStudioBackend):
         self._backend = backend
         self._oss = backend.oss
 
-    def cardinal_points(self, surface_1: int | None = None, surface_2: int | None = None) -> tuple[
-        CardinalPointsResult, zp.analyses.base.AnalysisResult]:
+    def cardinal_points(
+        self, surface_1: int | None = None, surface_2: int | None = None
+    ) -> tuple[CardinalPointsResult, zp.analyses.base.AnalysisResult]:
         """
         Get the cardinal points of the system between `surface_1` and `surface_2`.
 
@@ -170,11 +171,11 @@ def cardinal_points(self, surface_1: int | None = None, surface_2: int | None =
         return _build_cardinal_points_result(cardinal_points_result), cardinal_points_result
 
     def raytrace(
-            self,
-            coordinates: Iterable[tuple[float, float]],
-            wavelengths: Iterable[float] = (0.543,),
-            field_type: Literal["angle", "object_height"] = "angle",
-            pupil: tuple[float, float] = (0, 0),
+        self,
+        coordinates: Iterable[tuple[float, float]],
+        wavelengths: Iterable[float] = (0.543,),
+        field_type: Literal["angle", "object_height"] = "angle",
+        pupil: tuple[float, float] = (0, 0),
     ) -> tuple[DataFrame, list[zp.analyses.base.AnalysisResult]]:
         """
         Perform a ray trace analysis using the given parameters.
@@ -235,12 +236,12 @@ def raytrace(
         return _build_raytrace_result(raytrace_results), raytrace_results
 
     def zernike_standard_coefficients(
-            self,
-            field_coordinate: tuple[float, float] | None = None,
-            wavelength: float | None = None,
-            field_type: Literal["angle", "object_height"] = "angle",
-            sampling: str = "512x512",
-            maximum_term: int = 45,
+        self,
+        field_coordinate: tuple[float, float] | None = None,
+        wavelength: float | None = None,
+        field_type: Literal["angle", "object_height"] = "angle",
+        sampling: str = "512x512",
+        maximum_term: int = 45,
     ) -> tuple[zp.analyses.base.AttrDict, zp.analyses.base.AnalysisResult]:
         """
         Calculates the Zernike standard coefficients at the retina surface.
@@ -291,13 +292,13 @@ def zernike_standard_coefficients(
         return zernike_result.Data, zernike_result
 
     def refraction(
-            self,
-            field_coordinate: tuple[float, float] | None = None,
-            wavelength: float | None = None,
-            pupil_diameter: float | None = None,
-            field_type: Literal["angle", "object_height"] = "angle",
-            *,
-            use_higher_order_aberrations: bool = True,
+        self,
+        field_coordinate: tuple[float, float] | None = None,
+        wavelength: float | None = None,
+        pupil_diameter: float | None = None,
+        field_type: Literal["angle", "object_height"] = "angle",
+        *,
+        use_higher_order_aberrations: bool = True,
     ) -> tuple[FourierPowerVectorRefraction, zp.analyses.base.AnalysisResult]:
         """Calculates the ocular refraction.
 
diff --git a/visisipy/opticstudio/surfaces.py b/visisipy/opticstudio/surfaces.py
index 1ea920e..c5104cb 100644
--- a/visisipy/opticstudio/surfaces.py
+++ b/visisipy/opticstudio/surfaces.py
@@ -71,15 +71,15 @@ class OpticStudioSurface(BaseSurface):
     _TYPE: str = "Standard"
 
     def __init__(
-            self,
-            comment: str,
-            *,
-            radius: float = float("inf"),
-            thickness: float = 0.0,
-            semi_diameter: float | None = None,
-            conic: float = 0.0,
-            material: MaterialModel | str | None = None,
-            is_stop: bool | None = None,
+        self,
+        comment: str,
+        *,
+        radius: float = float("inf"),
+        thickness: float = 0.0,
+        semi_diameter: float | None = None,
+        conic: float = 0.0,
+        material: MaterialModel | str | None = None,
+        is_stop: bool | None = None,
     ):
         self._comment = comment
         self._radius = radius
@@ -237,18 +237,18 @@ def __new__(cls, *args, **kwargs):  # noqa: ARG003
         return super().__new__(cls)
 
     def __init__(
-            self,
-            comment: str,
-            *,
-            radius: float = float("inf"),
-            thickness: float = 0.0,
-            semi_diameter: float | None = None,
-            conic: float = 0.0,
-            material: MaterialModel | str | None = None,
-            is_stop: bool | None = None,
-            number_of_terms: int = 0,
-            norm_radius: float = 100,
-            zernike_coefficients: ZernikeCoefficients | None = None,
+        self,
+        comment: str,
+        *,
+        radius: float = float("inf"),
+        thickness: float = 0.0,
+        semi_diameter: float | None = None,
+        conic: float = 0.0,
+        material: MaterialModel | str | None = None,
+        is_stop: bool | None = None,
+        number_of_terms: int = 0,
+        norm_radius: float = 100,
+        zernike_coefficients: ZernikeCoefficients | None = None,
     ):
         super().__init__(
             comment=comment,
@@ -347,24 +347,25 @@ def build(self, oss: OpticStudioSystem, *, position: int, replace_existing: bool
 
 class OpticStudioZernikeStandardSagSurface(BaseOpticStudioZernikeSurface):
     """Zernike Standard Sag surface in OpticStudio."""
+
     _TYPE = "ZernikeStandardSag"
 
     def __init__(
-            self,
-            comment: str,
-            *,
-            radius: float = float("inf"),
-            thickness: float = 0.0,
-            semi_diameter: float | None = None,
-            conic: float = 0.0,
-            material: MaterialModel | str | None = None,
-            is_stop: bool | None = None,
-            extrapolate: int = 0,
-            zernike_decenter_x: float = 0.0,
-            zernike_decenter_y: float = 0.0,
-            number_of_terms: int = 0,
-            norm_radius: float = 100,
-            zernike_coefficients: ZernikeCoefficients | None = None,
+        self,
+        comment: str,
+        *,
+        radius: float = float("inf"),
+        thickness: float = 0.0,
+        semi_diameter: float | None = None,
+        conic: float = 0.0,
+        material: MaterialModel | str | None = None,
+        is_stop: bool | None = None,
+        extrapolate: int = 0,
+        zernike_decenter_x: float = 0.0,
+        zernike_decenter_y: float = 0.0,
+        number_of_terms: int = 0,
+        norm_radius: float = 100,
+        zernike_coefficients: ZernikeCoefficients | None = None,
     ):
         super().__init__(
             comment=comment,
@@ -412,23 +413,24 @@ def build(self, oss: OpticStudioSystem, *, position: int, replace_existing: bool
 
 class OpticStudioZernikeStandardPhaseSurface(BaseOpticStudioZernikeSurface):
     """Zernike Standard Phase surface in OpticStudio."""
+
     _TYPE = "ZernikeStandardPhase"
 
     def __init__(
-            self,
-            comment: str,
-            *,
-            radius: float = float("inf"),
-            thickness: float = 0.0,
-            semi_diameter: float | None = None,
-            conic: float = 0.0,
-            material: MaterialModel | str | None = None,
-            is_stop: bool | None = None,
-            extrapolate: int = 0,
-            diffract_order: float = 0.0,
-            number_of_terms: int = 0,
-            norm_radius: float = 100,
-            zernike_coefficients: ZernikeCoefficients | None = None,
+        self,
+        comment: str,
+        *,
+        radius: float = float("inf"),
+        thickness: float = 0.0,
+        semi_diameter: float | None = None,
+        conic: float = 0.0,
+        material: MaterialModel | str | None = None,
+        is_stop: bool | None = None,
+        extrapolate: int = 0,
+        diffract_order: float = 0.0,
+        number_of_terms: int = 0,
+        norm_radius: float = 100,
+        zernike_coefficients: ZernikeCoefficients | None = None,
     ):
         super().__init__(
             comment=comment,
@@ -495,9 +497,9 @@ def make_surface(surface: Surface, material: str | MaterialModel, comment: str =
 
 @make_surface.register
 def _make_surface(
-        surface: StandardSurface,
-        material: Union[str, MaterialModel],  # noqa: UP007
-        comment: str = "",
+    surface: StandardSurface,
+    material: Union[str, MaterialModel],  # noqa: UP007
+    comment: str = "",
 ) -> OpticStudioSurface:
     return OpticStudioSurface(
         comment=comment,
@@ -511,9 +513,10 @@ def _make_surface(
 
 @make_surface.register
 def _make_surface(
-        surface: Stop,
-        material: Union[str, MaterialModel] = "",  # noqa: UP007
-        comment: str = "") -> OpticStudioSurface:
+    surface: Stop,
+    material: Union[str, MaterialModel] = "",  # noqa: UP007
+    comment: str = "",
+) -> OpticStudioSurface:
     return OpticStudioSurface(
         comment=comment,
         thickness=surface.thickness,
@@ -525,9 +528,9 @@ def _make_surface(
 
 @make_surface.register
 def _make_surface(
-        surface: ZernikeStandardSagSurface,
-        material: Union[str, MaterialModel] = "",  # noqa: UP007
-        comment: str = ""
+    surface: ZernikeStandardSagSurface,
+    material: Union[str, MaterialModel] = "",  # noqa: UP007
+    comment: str = "",
 ) -> OpticStudioZernikeStandardSagSurface:
     return OpticStudioZernikeStandardSagSurface(
         comment=comment,
@@ -547,9 +550,9 @@ def _make_surface(
 
 @make_surface.register
 def _make_surface(
-        surface: ZernikeStandardPhaseSurface,
-        material: Union[str, MaterialModel] = "",  # noqa: UP007
-        comment: str = ""
+    surface: ZernikeStandardPhaseSurface,
+    material: Union[str, MaterialModel] = "",  # noqa: UP007
+    comment: str = "",
 ) -> OpticStudioZernikeStandardPhaseSurface:
     return OpticStudioZernikeStandardPhaseSurface(
         comment=comment,
diff --git a/visisipy/refraction.py b/visisipy/refraction.py
index 530b856..a150c93 100644
--- a/visisipy/refraction.py
+++ b/visisipy/refraction.py
@@ -57,12 +57,12 @@ def to_polar_power_vectors(self) -> PolarPowerVectorRefraction:
         """
         return PolarPowerVectorRefraction(
             M=self.M,
-            J=np.sqrt(self.J0 ** 2 + self.J45 ** 2),
+            J=np.sqrt(self.J0**2 + self.J45**2),
             axis=np.rad2deg(np.arctan2(self.J45, self.J0) / 2),
         )
 
     def to_sphero_cylindrical(
-            self, cylinder_form: Literal["positive", "negative"] = "negative"
+        self, cylinder_form: Literal["positive", "negative"] = "negative"
     ) -> SpheroCylindricalRefraction:
         """Converts the refraction to sphero-cylindrical form.
 
@@ -80,8 +80,8 @@ def to_sphero_cylindrical(
             raise ValueError("cylinder_form must be either 'positive' or 'negative'.")
 
         sphero_cylinder = SpheroCylindricalRefraction(
-            sphere=self.M + np.sqrt(self.J0 ** 2 + self.J45 ** 2),
-            cylinder=-2 * np.sqrt(self.J0 ** 2 + self.J45 ** 2),
+            sphere=self.M + np.sqrt(self.J0**2 + self.J45**2),
+            cylinder=-2 * np.sqrt(self.J0**2 + self.J45**2),
             axis=np.rad2deg(np.arctan2(self.J45, self.J0) / 2),
         )