Repository: LouisDesdoigts/dLux (language: python)

Issue:
Nd pixel positions -> Plane summaries



Discussion / hints:

Fix patch:
diff --git a/dLux/__init__.py b/dLux/__init__.py
index fcd54c6c..2f93ccc5 100644
--- a/dLux/__init__.py
+++ b/dLux/__init__.py
@@ -1,5 +1,5 @@
 name = "dLux"
-__version__ = "0.10.0"
+__version__ = "0.10.1"
 
 # Import as modules
 from . import core
@@ -32,4 +32,6 @@
 # Check for 64-bit
 from jax import config
 if not config.x64_enabled:
-    print("dLux: Jax is running in 32-bit, to enable 64-bit visit:\nhttps://jax.readthedocs.io/en/latest/notebooks/Common_Gotchas_in_JAX.html#double-64bit-precision")
+    print("dLux: Jax is running in 32-bit, to enable 64-bit visit: "
+          "https://jax.readthedocs.io/en/latest/notebooks/"
+          "Common_Gotchas_in_JAX.html#double-64bit-precision")
diff --git a/dLux/apertures.py b/dLux/apertures.py
index 75737b91..71c37b31 100644
--- a/dLux/apertures.py
+++ b/dLux/apertures.py
@@ -4,8 +4,10 @@
 from jax import numpy as np, lax, vmap
 from jax.tree_util import tree_map, tree_flatten
 from equinox import filter, static_field
-from dLux.utils import get_pixel_coordinates, coordinates as c, opd_to_phase, \
+from dLux.utils import get_pixel_positions, coordinates as c, opd_to_phase, \
     factorial, cartesian_to_polar, list_to_dictionary
+from dLux.utils.helpers import two_image_plot
+from dLux.utils.units import convert_angular, convert_cartesian
 
 
 Array = np.ndarray
@@ -87,7 +89,9 @@ def get_aperture(self     : ApertureLayer,
         aperture : Array 
             The array representing the transmission of the aperture.
         """
-        coordinates = get_pixel_coordinates(npixels, diameter / npixels)
+        npixels_in = (npixels, npixels)
+        pixel_scales = (diameter / npixels, diameter / npixels)
+        coordinates = get_pixel_positions(npixels_in, pixel_scales)
         return self._aperture(coordinates)
 
 
@@ -447,10 +451,6 @@ def __init__(self        : ApertureLayer,
                  occulting   : bool = False, 
                  softening   : Array = np.array(1.),
                  name        : str = "CircularAperture",
-
-                 static = False,
-                 npixels = None,
-                 diameter = None,
                  ) -> Array:
         """
         Constructor for the CircularAperture class.
@@ -532,6 +532,47 @@ def _extent(self : ApertureLayer) -> Array:
             The maximum distance from the centre to edge of aperture.
         """
         return self.radius
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        radius = convert_cartesian(self.radius, 'meters', cartesian_units)
+        center = convert_cartesian(self.center, 'meters', cartesian_units)
+        transmissive = "transmissive" if not self.occulting else "occulting"
+
+        summary = (f"Applies a {transmissive} Circular Aperture with radius "
+                   f"{radius} {cartesian_units}")
+        
+        if self.softening != np.array(0):
+            summary += f" softened by ~{self.softening} pixels"
+        if self.center != np.array([0., 0.]):
+            summary += f" centred at {center}"
+        if self.shear != np.array([0., 0.]):
+            summary += f" sheared by {self.shear}"
+        if self.compression != np.array([1., 1.]):
+            summary += f" compressed by {self.compression}"
+        return summary + "."
 
 
 class AnnularAperture(DynamicAperture):
@@ -659,6 +700,49 @@ def _extent(self : ApertureLayer) -> Array:
             The maximum distance from the centre to edge of aperture.
         """
         return self.rmax
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        rmin = convert_cartesian(self.rmin, 'meters', cartesian_units)
+        rmax = convert_cartesian(self.rmax, 'meters', cartesian_units)
+        center = convert_cartesian(self.center, 'meters', cartesian_units)
+        transmissive = "transmissive" if not self.occulting else "occulting"
+
+        summary = (f"Applies a {transmissive} Annular Aperture with inner "
+                   f"radius {rmin} {cartesian_units} and outer radius {rmax} "
+                   f"{cartesian_units}")
+        
+        if self.softening != np.array(0):
+            summary += f" softened by ~{self.softening} pixels"
+        if self.center != np.array([0., 0.]):
+            summary += f" centred at {center}"
+        if self.shear != np.array([0., 0.]):
+            summary += f" sheared by {self.shear}"
+        if self.compression != np.array([1., 1.]):
+            summary += f" compressed by {self.compression}"
+        return summary + "."
 
 
 class RectangularAperture(DynamicAperture):
@@ -792,6 +876,52 @@ def _extent(self : ApertureLayer) -> Array:
             The maximum distance from the centre to edge of aperture.
         """
         return np.hypot(self.height / 2., self.width / 2.)
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        height = convert_cartesian(self.height, 'meters', cartesian_units)
+        width = convert_cartesian(self.width, 'meters', cartesian_units)
+        center = convert_cartesian(self.center, 'meters', cartesian_units)
+        rotation = convert_angular(self.rotation, 'radians', angular_units)
+        transmissive = "transmissive" if not self.occulting else "occulting"
+
+        summary = (f"Applies a {transmissive} Rectangular Aperture with height "
+                   f"{height} {cartesian_units} and width {width} "
+                   f"{cartesian_units}")
+        
+        if self.softening != np.array(0):
+            summary += f" softened by ~{self.softening} pixels"
+        if self.center != np.array([0., 0.]):
+            summary += f" centred at {center}"
+        if self.rotation != np.array(0.):
+            summary += f" rotated by {rotation} {angular_units}"
+        if self.shear != np.array([0., 0.]):
+            summary += f" sheared by {self.shear}"
+        if self.compression != np.array([1., 1.]):
+            summary += f" compressed by {self.compression}"
+        return summary + "."
 
 
 class SquareAperture(DynamicAperture):
@@ -917,6 +1047,50 @@ def _extent(self : ApertureLayer) -> Array:
             The maximum distance from the centre to edge of aperture.
         """
         return np.sqrt(2) * self.width / 2.
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        width = convert_cartesian(self.width, 'meters', cartesian_units)
+        center = convert_cartesian(self.center, 'meters', cartesian_units)
+        rotation = convert_angular(self.rotation, 'radians', angular_units)
+        transmissive = "transmissive" if not self.occulting else "occulting"
+
+        summary = (f"Applies a {transmissive} Rectangular Aperture with width "
+                   f"{width} {cartesian_units}")
+        
+        if self.softening != np.array(0):
+            summary += f" softened by ~{self.softening} pixels"
+        if self.center != np.array([0., 0.]):
+            summary += f" centred at {center}"
+        if self.rotation != np.array(0.):
+            summary += f" rotated by {rotation} {angular_units}"
+        if self.shear != np.array([0., 0.]):
+            summary += f" sheared by {self.shear}"
+        if self.compression != np.array([1., 1.]):
+            summary += f" compressed by {self.compression}"
+        return summary + "."
 
 
 class PolygonalAperture(DynamicAperture, ABC):
@@ -1301,6 +1475,48 @@ def _hard_edged(self : ApertureLayer, coordinates : Array) -> Array:
         edges = (dist_from_edges * dist_sgn) > 0.
 
         return (edges).prod(axis=0)
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        center = convert_cartesian(self.center, 'meters', cartesian_units)
+        rotation = convert_angular(self.rotation, 'radians', angular_units)
+        transmissive = "transmissive" if not self.occulting else "occulting"
+
+        summary = f"Applies a {transmissive} Irregular Polygonal Aperture"
+        
+        if self.softening != np.array(0):
+            summary += f" softened by ~{self.softening} pixels"
+        if self.center != np.array([0., 0.]):
+            summary += f" centred at {center}"
+        if self.rotation != np.array(0.):
+            summary += f" rotated by {rotation} {angular_units}"
+        if self.shear != np.array([0., 0.]):
+            summary += f" sheared by {self.shear}"
+        if self.compression != np.array([1., 1.]):
+            summary += f" compressed by {self.compression}"
+        return summary + "."
 
 
 class RegularPolygonalAperture(PolygonalAperture):
@@ -1465,6 +1681,51 @@ def _hard_edged(self : ApertureLayer, coordinates : Array) -> Array:
         return dist.prod(axis=0)
 
 
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        rmax = convert_cartesian(self.rmax, 'meters', cartesian_units)
+        center = convert_cartesian(self.center, 'meters', cartesian_units)
+        rotation = convert_angular(self.rotation, 'radians', angular_units)
+        transmissive = "transmissive" if not self.occulting else "occulting"
+
+        summary = (f"Applies a {transmissive} {self.nsides} sided Regular "
+                   f"Polygonal Aperture of max radius {rmax:.{sigfigs}} "
+                   f"{cartesian_units}")
+        
+        if self.softening != np.array(0):
+            summary += f" softened by ~{self.softening} pixels"
+        if self.center != np.array([0., 0.]):
+            summary += f" centred at {center}"
+        if self.rotation != np.array(0.):
+            summary += f" rotated by {rotation} {angular_units}"
+        if self.shear != np.array([0., 0.]):
+            summary += f" sheared by {self.shear}"
+        if self.compression != np.array([1., 1.]):
+            summary += f" compressed by {self.compression}"
+        return summary + "."
+
+
 class HexagonalAperture(RegularPolygonalAperture):
     """
     A hexagonal aperture parameterised by the maximum radius to the vertices 
@@ -1538,6 +1799,48 @@ def __init__(self        : ApertureLayer,
                          name = name)
 
 
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        rmax = convert_cartesian(self.rmax, 'meters', cartesian_units)
+        center = convert_cartesian(self.center, 'meters', cartesian_units)
+        rotation = convert_angular(self.rotation, 'radians', angular_units)
+        transmissive = "transmissive" if not self.occulting else "occulting"
+
+        summary = (f"Applies a {transmissive} Hexagonal Aperture of max radius "
+                   f"{rmax:.{sigfigs}} {cartesian_units}")
+        
+        if self.softening != np.array(0):
+            summary += f" softened by ~{self.softening} pixels"
+        if self.center != np.array([0., 0.]):
+            summary += f" centred at {center}"
+        if self.rotation != np.array(0.):
+            summary += f" rotated by {rotation} {angular_units}"
+        if self.shear != np.array([0., 0.]):
+            summary += f" sheared by {self.shear}"
+        if self.compression != np.array([1., 1.]):
+            summary += f" compressed by {self.compression}"
+        return summary + "."
 
 
 ###############
@@ -1778,6 +2081,48 @@ def _hard_edged(self : ApertureLayer, coordinates : Array) -> Array:
         return struts.prod(axis=0)
 
 
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        width = convert_cartesian(self.strut_width, 'meters', cartesian_units)
+        center = convert_cartesian(self.center, 'meters', cartesian_units)
+        rotation = convert_angular(self.rotation, 'radians', angular_units)
+
+        summary = (f"Applies a {self.nstrut} strut spider with widths {width} "
+                   f"{cartesian_units}")
+        
+        if self.softening != np.array(0):
+            summary += f" softened by ~{self.softening} pixels"
+        if self.center != np.array([0., 0.]):
+            summary += f" centred at {center}"
+        if self.rotation != np.array(0.):
+            summary += f" rotated by {rotation} {angular_units}"
+        if self.shear != np.array([0., 0.]):
+            summary += f" sheared by {self.shear}"
+        if self.compression != np.array([1., 1.]):
+            summary += f" compressed by {self.compression}"
+        return summary + "."
+
 
 
 ###################
@@ -2041,7 +2386,9 @@ def get_aperture(self     : ApertureLayer,
         aperture : Array 
             The array representing the transmission of the aperture.
         """
-        coordinates = get_pixel_coordinates(npixels, diameter / npixels)
+        npixels_in = (npixels, npixels)
+        pixel_scales = (diameter / npixels, diameter / npixels)
+        coordinates = get_pixel_positions(npixels_in, pixel_scales)
         return self.aperture._aperture(coordinates)
 
 
@@ -2096,7 +2443,9 @@ def get_basis(self     : ApertureLayer,
         basis : Array 
             The array of the basis vectors of the aperture aberrations.
         """
-        coordinates = get_pixel_coordinates(npixels, diameter / npixels)
+        npixels_in = (npixels, npixels)
+        pixel_scales = (diameter / npixels, diameter / npixels)
+        coordinates = get_pixel_positions(npixels_in, pixel_scales)
         return self._basis(coordinates)
  
 
@@ -2137,7 +2486,9 @@ def get_opd(self : ApertureLayer, npixels : int, diameter : float) -> Array:
         basis : Array 
             The array of the total opd of the aperture aberrations.
         """
-        coordinates = get_pixel_coordinates(npixels, diameter / npixels)
+        npixels_in = (npixels, npixels)
+        pixel_scales = (diameter / npixels, diameter / npixels)
+        coordinates = get_pixel_positions(npixels_in, pixel_scales)
         return self._opd(coordinates)
 
 
@@ -2409,7 +2760,33 @@ def _orthonormalise(self     : ApertureLayer,
                     np.sqrt((intermediate ** 2).sum() / pixel_area))
         
         return basis
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
 
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        summary = super().summary(angular_units, cartesian_units, sigfigs)
+        return summary[:-1] + f" with {len(self.coefficients)} aberrations."
 
 
 
@@ -2548,7 +2925,9 @@ def get_aperture(self     : ApertureLayer,
         aperture : Array 
             The array representing the transmission of the aperture.
         """
-        coordinates = get_pixel_coordinates(npixels, diameter/npixels)
+        npixels_in = (npixels, npixels)
+        pixel_scales = (diameter / npixels, diameter / npixels)
+        coordinates = get_pixel_positions(npixels_in, pixel_scales)
         return self._aperture(coordinates)
 
 
@@ -2594,7 +2973,9 @@ def get_basis(self     : ApertureLayer,
         basis : Array 
             The array of the basis vectors of the aperture aberrations.
         """
-        coordinates = get_pixel_coordinates(npixels, diameter/npixels)
+        npixels_in = (npixels, npixels)
+        pixel_scales = (diameter / npixels, diameter / npixels)
+        coordinates = get_pixel_positions(npixels_in, pixel_scales)
         return self._basis(coordinates)
 
 
@@ -2685,7 +3066,9 @@ def get_opd(self : ApertureLayer, npixels : int, diameter : float) -> Array:
         basis : Array 
             The array of the total opd of the aperture aberrations.
         """
-        coordinates = get_pixel_coordinates(npixels, diameter/npixels)
+        npixels_in = (npixels, npixels)
+        pixel_scales = (diameter / npixels, diameter / npixels)
+        coordinates = get_pixel_positions(npixels_in, pixel_scales)
         return self._opd(coordinates)
 
 
@@ -2824,6 +3207,37 @@ def _aperture(self : ApertureLayer, coordinates : Array) -> Array:
         """
         aps = self._stacked_apertures(coordinates)
         return aps.prod(axis=0)
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        summary = ("Applies a Compound Aperture with the following "
+                   "sub-apertures: \n")
+        for ap in self.apertures:
+            ap_summary = ap.summary(angular_units, cartesian_units, sigfigs)
+            summary += ap_summary + "\n"
+        return summary
 
 
 class MultiAperture(CompositeAperture):
@@ -2932,6 +3346,35 @@ def is_aberrated(leaf):
         return tree_flatten(aberrated, is_leaf=is_aberrated)[0]
 
 
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        summary = ("Applies a Multi Aperture with the following "
+                   "sub-apertures: \n")
+        for ap in self.apertures:
+            ap_summary = ap.summary(angular_units, cartesian_units, sigfigs)
+            summary += ap_summary + "\n"
+        return summary
 
 
 ########################
@@ -2987,7 +3430,9 @@ def __init__(self        : ApertureLayer,
         
         # Generate coordinates if not provided
         if coordinates is None:
-            coordinates = get_pixel_coordinates(npixels, diameter/npixels)
+            npixels_in = (npixels, npixels)
+            pixel_scales = (diameter / npixels, diameter / npixels)
+            coordinates = get_pixel_positions(npixels_in, pixel_scales)
 
         super().__init__(name = name)
         self.aperture = aperture._aperture(coordinates)
@@ -3098,6 +3543,33 @@ def __init__(self        : ApertureLayer,
                          diameter = diameter, 
                          coordinates = coordinates, 
                          name = name)
+        
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return "Applies a pre-calculated Static Aperture."
+
 
 
 class StaticAberratedAperture(AbstractAberratedAperture, AbstractStaticAperture):
@@ -3159,7 +3631,9 @@ def __init__(self        : ApertureLayer,
         
         # Generate coordinates if not provided
         if coordinates is None:
-            coordinates = get_pixel_coordinates(npixels, diameter/npixels)
+            npixels_in = (npixels, npixels)
+        pixel_scales = (diameter / npixels, diameter / npixels)
+        coordinates = get_pixel_positions(npixels_in, pixel_scales)
 
         super().__init__(aperture=aperture, coordinates=coordinates, 
             coefficients=aperture.coefficients, name=name)
@@ -3274,6 +3748,33 @@ def opd(self : ApertureLayer) -> Array:
             The array of the total opd of the aperture aberrations.
         """
         return self._opd()
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"Applies a pre-calculated Static Aperture with "
+                f"{len(self.coefficients)} aberrations.")
 
 
 
diff --git a/dLux/core.py b/dLux/core.py
index 4a11411c..d0b6d48f 100644
--- a/dLux/core.py
+++ b/dLux/core.py
@@ -8,6 +8,7 @@
 from collections import OrderedDict
 from copy import deepcopy
 from functools import partial
+from typing import Union
 import dLux
 
 
@@ -438,6 +439,39 @@ def normalise(self : Instrument) -> Instrument:
         """
         return tree_at(lambda instrument: instrument.scene, self, \
                            self.scene.normalise())
+    
+
+    def summarise(self : Instrument) -> None:
+        """
+        Prints a summary of all instrument
+        """
+        if self.scene is not None:
+            print("Scene summary:")
+            self.scene.summarise()
+        print("Optics summary:")
+        self.optics.summarise()
+        if self.detector is not None:
+            print("Detector summary:")
+            self.detector.summarise()
+
+
+    def plot(self       : Optics, 
+             wavelength : Array, 
+             offset     : Array = np.zeros(2)) -> None:
+        """
+        Prints the summary of all the planes and then plots a wavefront as it
+        propagates through the optics.
+
+        Parameters
+        ----------
+        wavelength : Array, meters
+            The wavelength of the wavefront to propagate through the optical
+            layers.
+        offset : Array, radians = np.zeros(2)
+            The (x, y) offset from the optical axis of the source.
+        """
+        wf = self.optics.plot(wavelength, offset)
+        self.detector.plot(wf.psf)
 
 
     def model(self : Instrument, **kwargs):
@@ -720,6 +754,80 @@ def debug_prop(self       : Optics,
 
         return params_dict["Wavefront"].wavefront_to_psf(), \
                                 intermediate_dicts, intermediate_layers
+    
+
+    def get_planes(self : Optics) -> list:
+        """
+        Breaks the optical layers into planes, where each plane is a list of
+        layers.
+
+        Returns
+        -------
+        planes : list
+            A list of lists, with the inner lists being optical layers, and the
+            outer list being planes.
+        """
+        planes = []
+        plane = []
+        keys = self.layers.keys()
+        for key in keys:
+            layer = self.layers[key]
+            plane.append(layer)
+            if isinstance(layer, dLux.propagators.Propagator):
+                planes.append(plane)
+                plane = []
+        return planes
+
+
+    def summarise(self : Optics) -> None:
+        """
+        Prints a summary of all the planes in the optical system.
+        """
+        planes = self.get_planes()
+        # TODO: Add plane type (Plane 0: Pupil)
+        print("Text summary:")
+        for i in range(len(planes)):
+            print(f'Plane {i}')
+            for layer in planes[i]:
+                print(f"  {layer.summary(angular_units='arcseconds')}")
+        print('\n')
+
+
+    def plot(self       : Optics, 
+             wavelength : Array, 
+             offset     : Array = np.zeros(2)) -> None:
+        """
+        Prints the summary of all the planes and then plots a wavefront as it
+        propagates through the optics.
+
+        Parameters
+        ----------
+        wavelength : Array, meters
+            The wavelength of the wavefront to propagate through the optical
+            layers.
+        offset : Array, radians = np.zeros(2)
+            The (x, y) offset from the optical axis of the source.
+        
+        Returns
+        -------
+        wf : Wavefront
+            The final wavefront after being propagated through the optical
+            layers.
+        """
+        planes = self.get_planes()
+        self.summarise()
+
+        for i in range(len(planes)):
+            print(f'Plane {i}')
+            for layer in planes[i]:
+                print(f"  {layer.summary()}")
+                if isinstance(layer, dLux.CreateWavefront):
+                    wf, _ = layer(None, 
+                                  {"wavelength": wavelength, 'offset': offset})
+                else:
+                    wf = layer(wf)
+                layer.display(wf)
+        return wf
 
 
     def model(self : Optics, **kwargs):
@@ -839,6 +947,32 @@ def normalise(self : Scene) -> Scene:
                 is_leaf = lambda leaf: isinstance(leaf, dLux.sources.Source))
 
         return tree_at(lambda scene: scene.sources, self, normalised_sources)
+    
+
+    def summarise(self : Scene) -> None:
+        """
+        Prints a summary of all the Scene.
+        """
+        for source in self.sources.values():
+            print(source.summary())
+
+
+    def plot(self : Scene) -> None:
+        """
+        Prints the summary of all the planes and then plots a wavefront as it
+        propagates through the optics.
+
+        Parameters
+        ----------
+        wavelength : Array, meters
+            The wavelength of the wavefront to propagate through the optical
+            layers.
+        offset : Array, radians = np.zeros(2)
+            The (x, y) offset from the optical axis of the source.
+        """
+        for source in self.sources.values():
+            print(source.summary())
+            source.display()
 
 
     def model(self : Scene, optics : Optics, **kwargs):
@@ -997,6 +1131,37 @@ def debug_apply_detector(self  : Instrument,
             intermediate_images.append(image.copy())
             intermediate_layers.append(deepcopy(layer))
         return image, intermediate_images, intermediate_layers
+    
+
+    def summarise(self : Detector) -> None:
+        """
+        Prints a summary of all the layers in the detector.
+        """
+        print("Text summary:")
+        keys = self.layers.keys()
+        for key in keys:
+            layer = self.layers[key]
+            print(f"  {layer.summary()}")
+        print('\n')
+
+
+    def plot(self : Optics, image : Array) -> None:
+        """
+        Prints the summary of all the layers and then plots a image as it
+        propagates through the detector layer.
+
+        Parameters
+        ----------
+        iamge : Array
+            The image to propagate through the detector.
+        """
+        self.summarise()
+        keys = self.layers.keys()
+        for key in keys:
+            layer = self.layers[key]
+            print(f"{layer.summary()}")
+            image = layer(image)
+            layer.display(image)
 
 
     def model(self : Detector, optics : Optics, **kwargs):
diff --git a/dLux/detectors.py b/dLux/detectors.py
index a625d97c..1f24035d 100644
--- a/dLux/detectors.py
+++ b/dLux/detectors.py
@@ -5,8 +5,10 @@
 from jax.scipy.stats import norm
 from equinox import tree_at, static_field
 from zodiax import ExtendedBase
-import dLux
 from dLux.utils.interpolation import rotate, fourier_rotate
+from dLux.utils.units import convert_angular, convert_cartesian
+from dLux.utils.helpers import single_image_plot
+import dLux
 
 
 __all__ = ["ApplyPixelResponse", "ApplyJitter", "ApplySaturation",
@@ -48,6 +50,65 @@ def __call__(self : DetectorLayer, image : Array) -> Array: # pragma: no cover
         Abstract method for Detector Layers
         """
         return
+    
+
+    def summary(self            : DetectorLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return f"{self.name} layer has no summary method yet."
+    
+
+    def display(self            : DetectorLayer, 
+                image           : Array,
+                figsize         : tuple = (5, 4),
+                dpi             : int = 120,
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> None:
+        """
+        Displays a plot of the image propagating through the layer.
+
+        Parameters
+        ----------
+        image : Array
+            The dummy image to propagate though the detector.
+        figsize : tuple = (10, 4)
+            The size of the figure to display.
+        cmap : str = 'inferno'
+            The colour map to use.
+        dpi : int = 120
+            The resolution of the figure.
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+        """
+        single_image_plot(image, figsize=figsize, title="PSF", 
+                          cbar_label="Counts", cmap='inferno', dpi=dpi)
 
 
 class ApplyPixelResponse(DetectorLayer):
@@ -99,6 +160,32 @@ def __call__(self : DetectorLayer, image) -> Array:
             The image with the pixel_response applied.
         """
         return image * self.pixel_response
+    
+
+    def summary(self            : DetectorLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return f"{self.name}: Applies the pixel response to the image."
 
 
 class ApplyJitter(DetectorLayer):
@@ -173,6 +260,33 @@ def __call__(self : DetectorLayer, image : Array) -> Array:
         """
         kernel = self.generate_kernel()
         return convolve(image, kernel, mode='same')
+    
+
+    def summary(self            : DetectorLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Applies a jitter of {self.jitter:.{sigfigs}} "
+                "pixels.")
 
 
 class ApplySaturation(DetectorLayer):
@@ -224,6 +338,33 @@ def __call__(self : DetectorLayer, image : Array) -> Array:
             The image with the saturation applied.
         """
         return np.minimum(image, self.saturation)
+    
+
+    def summary(self            : DetectorLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Applies a saturation of "
+                f"{self.saturation:.{sigfigs}} counts.")
 
 
 class AddConstant(DetectorLayer):
@@ -274,6 +415,33 @@ def __call__(self : DetectorLayer, image : Array) -> Array:
             The image with the value added.
         """
         return image + self.value
+    
+
+    def summary(self            : DetectorLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Adds a constant value of "
+                f"{self.value:.{sigfigs}} counts.")
 
 
 class IntegerDownsample(DetectorLayer):
@@ -353,6 +521,33 @@ def __call__(self, image):
             The downsampled image.
         """
         return self.downsample(image, self.kernel_size)
+    
+
+    def summary(self            : DetectorLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Downsamples the image by {self.kernel_size} "
+                f"pixels.")
 
 
 class Rotate(DetectorLayer):
@@ -420,4 +615,38 @@ def __call__(self : DetectorLayer, image : Array) -> Array:
         if self.fourier:
             return fourier_rotate(image, self.angle, self.padding)
         else:
-            return rotate(image, self.angle)
\ No newline at end of file
+            return rotate(image, self.angle)
+        
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        angle = convert_angular(self.angle, 'radians', angular_units)
+        
+        if self.fourier:
+            method = f"a Fourier method with padding of {self.padding}"
+        else:
+            method = "an Interpolation method of order 1"
+
+        return (f"{self.name}: Applies a {angle:.{sigfigs}} {angular_units} "
+                f"rotation to the image using {method}.")
\ No newline at end of file
diff --git a/dLux/optics.py b/dLux/optics.py
index 440e46c5..58823f10 100644
--- a/dLux/optics.py
+++ b/dLux/optics.py
@@ -8,13 +8,14 @@
 from zodiax import ExtendedBase
 from abc import ABC, abstractmethod
 from inspect import signature
+from dLux.utils.helpers import two_image_plot
+from dLux.utils.units import convert_angular, convert_cartesian
 import dLux
 
 
-__all__ = ["CreateWavefront", "TiltWavefront",
-           "NormaliseWavefront", "ApplyBasisOPD", "AddPhase", "AddOPD",
-           "TransmissiveOptic", "CompoundAperture", "ApplyBasisCLIMB",
-           "Rotate"]
+__all__ = ["CreateWavefront", "TiltWavefront", "NormaliseWavefront", 
+           "ApplyBasisOPD", "AddPhase", "AddOPD", "TransmissiveOptic", 
+           "ApplyBasisCLIMB", "Rotate", "CompoundAperture"]
 
 
 Array = np.ndarray
@@ -103,6 +104,69 @@ def apply(self : OpticalLayer, parameters : dict) -> dict:
 
         # Return updated parameters dictionary
         return parameters
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return f"{self.name} layer has no summary method yet."
+
+
+    def display(self            : OpticalLayer, 
+                wavefront       : Wavefront,
+                figsize         : tuple = (10, 4),
+                dpi             : int = 120,
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> None:
+        """
+        Displays a plot of the wavefront amplitude and opd or phase.
+
+        Parameters
+        ----------
+        wavefront : Wavefront
+            The dummy wavefront to propagate though the optics.
+        figsize : tuple = (10, 4)
+            The size of the figure to display.
+        cmap : str = 'inferno'
+            The colour map to use.
+        dpi : int = 120
+            The resolution of the figure.
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+        """
+        amplitude = wavefront.amplitude
+        phase = convert_angular(wavefront.phase, "radians", angular_units)
+        two_image_plot(amplitude[0], phase[0], figsize=figsize, 
+            titles=("Amplitude", "Phase"), cbar_labels=("Intensity", 
+            f"Phase ({angular_units})"), cmaps=('inferno', 'twilight'), 
+            bounds=(None, 2*np.pi), dpi=dpi)
 
 
 class CreateWavefront(OpticalLayer):
@@ -224,6 +288,33 @@ def __call__(self               : OpticalLayer,
             return wavefront, parameters
         else:
             return wavefront
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Creates a {self.wavefront_type} wavefront of " 
+                f"size {self.npixels} pixels and diameter {self.diameter} m.")
 
 
 class TiltWavefront(OpticalLayer):
@@ -279,6 +370,34 @@ def __call__(self : OpticalLayer, wavefront : Wavefront) -> Wavefront:
         return wavefront.tilt_wavefront(self.tilt_angles)
 
 
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        angles = convert_angular(self.tilt_angles, 'radians', angular_units)
+        return (f"{self.name}: Tilts the wavefront by {angles:.{sigfigs}} "
+                f"{angular_units} in the (x, y) dimension.")
+
+
 class NormaliseWavefront(OpticalLayer):
     """
     Normalises the input wavefront using the in-built wavefront normalisation
@@ -320,6 +439,32 @@ def __call__(self : OpticalLayer, wavefront : Wavefront) -> Wavefront:
             The wavefront with the wavefront normalisation method applied.
         """
         return wavefront.normalise()
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return f"{self.name}: Normalises the wavefront to unity power."
 
 
 class ApplyBasisOPD(OpticalLayer):
@@ -406,6 +551,33 @@ def __call__(self : OpticalLayer, wavefront : Wavefront) -> Wavefront:
             The wavefront with the appropriate phase applied.
         """
         return wavefront.add_opd(self.get_total_opd())
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Applies an Optical Path Difference (OPD) to the "
+            f"wavefront calculated from the basis vectors and coefficients.")
 
 
 class AddPhase(OpticalLayer):
@@ -432,8 +604,8 @@ def __init__(self  : OpticalLayer,
         ----------
         phase : Array, radians
             Array of phase values to be applied to the input wavefront. This
-            must a 0, 2 or 3 dimensional array with equal to that of the wavefront
-            at time of aplication.
+            must a 0, 2 or 3 dimensional array with equal to that of the 
+            wavefront at time of aplication.
         name : str = 'AddPhase'
             The name of the layer, which is used to index the layers dictionary.
         """
@@ -460,6 +632,32 @@ def __call__(self : OpticalLayer, wavefront : Wavefront) -> Wavefront:
             The wavefront with the phase added.
         """
         return wavefront.add_phase(self.phase)
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return f"{self.name}: Add an array of phase values to the wavefront."
 
 
 class AddOPD(OpticalLayer):
@@ -514,6 +712,33 @@ def __call__(self : OpticalLayer, wavefront : Wavefront) -> Wavefront:
             The wavefront with the OPD added.
         """
         return wavefront.add_opd(self.opd)
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Add an array of Optical Path Differences (OPD) "
+                "to the wavefront.")
 
 
 class TransmissiveOptic(OpticalLayer):
@@ -571,6 +796,33 @@ def __call__(self : OpticalLayer, wavefront : Wavefront) -> Wavefront:
             The wavefront with the tranmission applied.
         """
         return wavefront.multiply_amplitude(self.transmission)
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Applies an array of tranmission values to the "
+                "Wavefront.")
 
 
 class CompoundAperture(OpticalLayer):
@@ -751,7 +1003,8 @@ def construct_combined_aperture(self     : OpticalLayer,
 
         # Generate coordinate grid
         pixel_scale = diameter/npixels
-        xycoords = dLux.utils.get_pixel_coordinates(npixels, pixel_scale)
+        xycoords = dLux.utils.get_pixel_positions((npixels, npixels), 
+                                                  (pixel_scale, pixel_scale))
 
         # Generate aperture/occulters
         outer_apers = mapped_aperture(self.aperture_radii, \
@@ -934,7 +1187,8 @@ def get_binary_phase(self):
 
 
     def lsq_params(self, img):
-        xx, yy = np.meshgrid(np.linspace(0,1,img.shape[0]),np.linspace(0,1,img.shape[1]))
+        xx, yy = np.meshgrid(np.linspace(0,1,img.shape[0]),
+                             np.linspace(0,1,img.shape[1]))
         A = np.vstack([xx.ravel(), yy.ravel(), np.ones_like(xx).ravel()]).T
         matrix = np.linalg.inv(np.dot(A.T,A)).dot(A.T)
         return matrix, xx, yy, A
@@ -947,7 +1201,9 @@ def lsq(self, img):
 
     def area(self, img, epsilon = 1e-15):
         a,b,c = self.lsq(img)
-        a, b, c = np.where(a==0,epsilon,a), np.where(b==0,epsilon,b), np.where(c==0,epsilon,c)
+        a = np.where(a==0,epsilon,a)
+        b = np.where(b==0,epsilon,b)
+        c = np.where(c==0,epsilon,c)
         x1 = (-b-c)/(a) # don't divide by zero
         x2 = -c/(a) # don't divide by zero
         x1, x2 = np.min(np.array([x1,x2])), np.max(np.array([x1,x2]))
@@ -986,6 +1242,33 @@ def CLIMB(self, wf, ppsz = 256):
         soft_bin = vmap_mask(flat).reshape(ppsz, ppsz)
 
         return soft_bin
+    
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Applies a binary OPD to the Wavefront using the "
+                "CLIMB algorithm.")
 
 
 class Rotate(OpticalLayer):
@@ -1001,6 +1284,9 @@ class Rotate(OpticalLayer):
         or the real and imaginary arrays.
     fourier : bool
         Should the rotation be done using fourier methods or interpolation.
+    order : int = 1
+        The order of the interpolation to use. Only applies if fourier is
+        False. Must be 0, 1, or 3.
     padding : int
         The amount of padding to use if the fourier method is used.
     name : str
@@ -1009,6 +1295,7 @@ class Rotate(OpticalLayer):
     angle          : Array
     real_imaginary : bool
     fourier        : bool
+    order          : int
     padding        : int
 
 
@@ -1016,6 +1303,7 @@ def __init__(self           : OpticalLayer,
                  angle          : Array,
                  real_imaginary : bool = False,
                  fourier        : bool = False,
+                 order          : int  = 1,
                  padding        : int  = None,
                  name           : str  = 'Rotate') -> OpticalLayer:
         """
@@ -1024,13 +1312,17 @@ def __init__(self           : OpticalLayer,
         Parameters
         ----------
         angle: float, radians
-            The angle by which to rotate the wavefront in the clockwise direction.
+            The angle by which to rotate the wavefront in the clockwise 
+            direction.
         real_imaginary : bool = False
             Should the rotation be performed on the amplitude and phase array
             or the real and imaginary arrays.
         fourier : bool = False
             Should the fourier rotation method be used (True), or regular
             interpolation method be used (False).
+        order : int = 1
+            The order of the interpolation to use. Only applies if fourier is
+            False. Must be 0, 1, or 3.
         padding : int = None
             The amount of fourier padding to use. Only applies if fourier is
             True.
@@ -1040,6 +1332,9 @@ def __init__(self           : OpticalLayer,
         super().__init__(name)
         self.angle          = np.asarray(angle, dtype=float)
         self.real_imaginary = bool(real_imaginary)
+        if order not in (0, 1, 3):
+            raise ValueError("Order must be 0, 1, or 3.")
+        self.order = int(order)
         self.fourier        = bool(fourier)
         self.padding = padding if padding is None else int(padding)
         assert self.angle.ndim == 0, ("angle must be scalar array.")
@@ -1059,6 +1354,45 @@ def __call__(self : OpticalLayer, wavefront : Wavefront) -> Wavefront:
         wavefront : Wavefront
             The rotated wavefront.
         """
-        args = [self.angle, self.real_imaginary, self.fourier]
+        args = [self.angle, self.real_imaginary, self.fourier, self.order]
         args += [self.padding] if self.padding is not None else []
         return wavefront.rotate(*args)
+
+
+    def summary(self            : OpticalLayer, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        angle = convert_angular(self.angle, 'radians', angular_units)
+        
+        if self.fourier:
+            method = f"a Fourier method with padding of {self.padding}"
+        else:
+            method = "an Interpolation method of order 1"
+
+        if self.real_imaginary:
+            wf_type = "real and imaginary arrays"
+        else:
+            wf_type = "amplitude and phase arrays"
+
+        return (f"{self.name}: Applies a {angle:.{sigfigs}} {angular_units} "
+                f"rotation to the wavefront {wf_type} using {method}.")
\ No newline at end of file
diff --git a/dLux/propagators.py b/dLux/propagators.py
index bbf49fcd..6fe7c1c7 100644
--- a/dLux/propagators.py
+++ b/dLux/propagators.py
@@ -2,7 +2,9 @@
 import jax.numpy as np
 from equinox import tree_at
 from abc import ABC, abstractmethod
-from dLux.utils.coordinates import get_pixel_coordinates, get_coordinates_vector
+from dLux.utils.coordinates import get_pixel_positions
+from dLux.utils.helpers import two_image_plot
+from dLux.utils.units import convert_angular, convert_cartesian
 import dLux
 
 
@@ -193,11 +195,11 @@ def _generate_transfer_matrices(self         : Propagator,
         pixels_input, npixels_out = npixels
         sign = 1 if self.inverse else -1
 
-        input_coordinates = get_coordinates_vector(pixels_input, input_scale,
-                                                   pixel_offset/input_scale)
+        input_coordinates = get_pixel_positions(pixels_input, input_scale,
+                                                    pixel_offset * input_scale)
 
-        output_coordinates = get_coordinates_vector(npixels_out, output_scale,
-                                                    pixel_offset/output_scale)
+        output_coordinates = get_pixel_positions(npixels_out, output_scale,
+                                                    pixel_offset * output_scale)
 
         input_to_output = np.outer(input_coordinates, output_coordinates)
 
@@ -569,7 +571,56 @@ def get_nfringes(self      : Propagator,
         size_in = wavefront.diameter
         size_out = self.pixel_scale_out * self.npixels_out
         return size_in * size_out / self.focal_length / wavefront.wavelength
+    
 
+    def summary(self            : Propagator, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        plane_out = 'pupil' if self.inverse else 'focal'
+
+        pixel_scale = convert_cartesian(self.pixel_scale_out, 'meters', 
+                                        cartesian_units)
+        focal_length = convert_cartesian(self.focal_length, 'meters', 
+                                        cartesian_units)
+        shift = convert_cartesian(self.shift, 'meters', cartesian_units) \
+            if not self.pixel_shift else self.shift
+
+        prop_string = (f"Propagates the Wavefront {focal_length} "
+                       f"{cartesian_units} to a {plane_out} plane "
+                       f"with {self.npixels_out}x{self.npixels_out} pixels of "
+                       f"size {pixel_scale:.{sigfigs}} {cartesian_units} ")
+        
+        if (self.shift == np.zeros(2)).all():
+            return prop_string +  "using a 2-sided MFT."
+        else:
+            if self.pixel_shift:
+                shift_string = (f"shifted by {shift} pixels in the (x, y) "
+                                "dimension ")
+            else:
+                shift_string = (f"shifted by {shift} {cartesian_units} in the "
+                                "(x, y) dimension ")
+            return prop_string + shift_string + "using a 2-sided MFT."
+        
 
 class AngularMFT(AngularPropagator, VariableSamplingPropagator):
     """
@@ -652,6 +703,52 @@ def get_nfringes(self      : Propagator,
         fringe_size = wavefront.wavelength / wavefront.diameter
         detector_size = self.npixels_out * self.pixel_scale_out
         return detector_size / fringe_size
+    
+
+    def summary(self            : Propagator, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        plane_out = 'pupil' if self.inverse else 'focal'
+            
+        pixel_scale = convert_angular(self.pixel_scale_out, 'radians', 
+                                        angular_units)
+        shift = convert_angular(self.shift, 'radians', angular_units) \
+            if not self.pixel_shift else self.shift
+
+        prop_string = (f"Propagates the Wavefront to a {plane_out} plane "
+                       f"with {self.npixels_out}x{self.npixels_out} pixels of "
+                       f"size {pixel_scale:.{sigfigs}} {angular_units} ")
+        
+        if (self.shift == np.zeros(2)).all():
+            return prop_string +  "using a 2-sided MFT."
+        else:
+            if self.pixel_shift:
+                shift_string = (f"shifted by {shift} pixels in the (x, y) "
+                                "dimension ")
+            else:
+                shift_string = (f"shifted by {shift} {angular_units} in the "
+                                "(x, y) dimension ")
+            return prop_string + shift_string + "using a 2-sided MFT."
 
 
 class CartesianFFT(CartesianPropagator, FixedSamplingPropagator):
@@ -708,6 +805,37 @@ def get_pixel_scale_out(self : Propagator, wavefront : Wavefront) -> Array:
             The pixel scale in the output plane.
         """
         return self.focal_length * wavefront.wavelength / wavefront.diameter
+    
+
+    def summary(self            : Propagator, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        plane_out = 'pupil' if self.inverse else 'focal'
+        focal_length = convert_cartesian(self.focal_length, 'meters', 
+                                        cartesian_units)
+
+        return (f"Propagates the Wavefront {focal_length} {cartesian_units} to "
+                f"a {plane_out} plane using a FFT.")
 
 
 class AngularFFT(AngularPropagator, FixedSamplingPropagator):
@@ -761,6 +889,34 @@ def get_pixel_scale_out(self : Propagator, wavefront : Wavefront) -> Array:
             The pixel scale in the output plane.
         """
         return wavefront.wavelength / wavefront.diameter
+    
+
+    def summary(self            : Propagator, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        plane_out = 'pupil' if self.inverse else 'focal'
+
+        return f"Propagates the Wavefront to a {plane_out} plane using a FFT."
 
 
 class CartesianFresnel(FarFieldFresnel, CartesianMFT):
@@ -935,8 +1091,9 @@ def propagate(self : Propagator, wavefront : Wavefront) -> Array:
         offsets = self.get_shift()
 
         input_positions = wavefront.pixel_coordinates
-        output_positions = get_pixel_coordinates(self.npixels_out,
-                                                 self.pixel_scale_out)
+        output_positions = get_pixel_positions(
+                            (self.npixels_out, self.npixels_out), 
+                            (self.pixel_scale_out, self.pixel_scale_out))
 
         propagation_distance = self.focal_length + self.propagation_shift
 
@@ -952,6 +1109,61 @@ def propagate(self : Propagator, wavefront : Wavefront) -> Array:
         field *= self.quadratic_phase(*output_positions,
             wavefront.wavelength, propagation_distance)
         return field
+    
+
+    def summary(self            : Propagator, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        plane_out = 'pupil' if self.inverse else 'focal'
+
+        pixel_scale = convert_cartesian(self.pixel_scale_out, 'meters', 
+                                        cartesian_units)
+        prop_distance = convert_cartesian(self.focal_length, 'meters', 
+                                        cartesian_units)
+        defocus = convert_cartesian(self.propagation_shift, 'meters', 
+                                        cartesian_units)
+        shift = convert_cartesian(self.shift, 'meters', cartesian_units) \
+            if not self.pixel_shift else self.shift
+
+
+        prop_string = (f"Propagates the Wavefront {prop_distance} "
+                       f"{cartesian_units} to a {plane_out} plane "
+                       f"with {self.npixels_out}x{self.npixels_out} pixels of "
+                       f"size {pixel_scale:.{sigfigs}} {cartesian_units} ")
+        
+        defocus_string = (f"with a defocus of {defocus} {cartesian_units} "
+                           "using a 2-sided MFT.")
+        
+        if (self.shift == np.zeros(2)).all():
+            return prop_string +  defocus_string
+        else:
+            if self.pixel_shift:
+                shift_string = (f"shifted by {shift} pixels in the (x, y) "
+                                "dimension ")
+            else:
+                shift_string = (f"shifted by {shift} {cartesian_units} in the "
+                                "(x, y) dimension ")
+            return prop_string + shift_string + defocus_string
 
 
 # TODO: Implement eventually
diff --git a/dLux/sources.py b/dLux/sources.py
index a4371532..d730a522 100644
--- a/dLux/sources.py
+++ b/dLux/sources.py
@@ -5,6 +5,8 @@
 from equinox import tree_at, static_field
 from zodiax import ExtendedBase
 from abc import ABC, abstractmethod
+from dLux.utils.units import convert_angular, convert_cartesian
+from dLux.utils.helpers import spectrum_plot
 import dLux
 
 
@@ -163,6 +165,62 @@ def normalise(self : Source) -> Source:
         normalised_spectrum = self.spectrum.normalise()
         return tree_at(
             lambda source : source.spectrum, self, normalised_spectrum)
+    
+
+    def summary(self            : Source, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name} has no summary method yet.")
+    
+
+    def display(self            : Source, 
+                figsize         : tuple = (6, 3),
+                dpi             : int = 120,
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> None:
+        """
+        Displays a plot of the wavefront amplitude and opd or phase.
+
+        Parameters
+        ----------
+        figsize : tuple = (6, 3)
+            The size of the figure to display.
+        cmap : str = 'inferno'
+            The colour map to use.
+        dpi : int = 120
+            The resolution of the figure.
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+        """
+        self.spectrum.normalise().display(figsize=figsize, dpi=dpi, 
+            angular_units=angular_units, cartesian_units=cartesian_units)
 
 
     def model(self      : Source,
@@ -458,6 +516,39 @@ def __init__(self        : Source,
         super().__init__(position, flux, spectrum, name=name)
 
 
+    def summary(self            : Source, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        position = np.array_str(convert_angular(self.position, 'radians', 
+                                            angular_units), precision=sigfigs)
+        spectrum_str = self.spectrum.summary(angular_units, cartesian_units,
+                                                sigfigs)
+        flux = self.flux
+        return (f"{self.name}: A Point Source with flux {flux:.{sigfigs}} "
+                f"photons at offset from the optical axis by {position} "
+                f"{angular_units} with spectrum\n  {spectrum_str}")
+
+
 class MultiPointSource(Source):
     """
     Concrete Class for multiple unresolved point source objects.
@@ -575,6 +666,40 @@ def model(self      : Source,
 
         # Apply detector if supplied
         return psf if detector is None else detector.apply_detector(psf)
+    
+
+    def summary(self            : Source, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        position = np.array_str(convert_angular(self.position, 'radians', 
+                                            angular_units), precision=sigfigs)
+        spectrum_str = self.spectrum.summary(angular_units, cartesian_units,
+                                                sigfigs)
+        n = len(self.flux)
+        return (f"{self.name}: {n} Point Sources with fluxes "
+                f"{self.flux:.{sigfigs}} photons at offsets from the optical "
+                f"axis by {position} {angular_units} with spectrum\n  "
+                f"{spectrum_str}")
 
 
 class ArrayDistribution(ResolvedSource):
@@ -673,6 +798,39 @@ def normalise(self : Source) -> Source:
         return tree_at(
             lambda source : (source.spectrum, source.distribution), self, \
                             (normalised_spectrum, normalised_distribution))
+    
+
+    def summary(self            : Source, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        position = np.array_str(convert_angular(self.position, 'radians', 
+                                            angular_units), precision=sigfigs)
+        spectrum_str = self.spectrum.summary(angular_units, cartesian_units,
+                                                sigfigs)
+        return (f"{self.name}: A Resolved Source with flux "
+                f"{self.flux:.{sigfigs}} photons at offset from the optical "
+                f"axis by {position} {angular_units} with spectrum\n  "
+                f"{spectrum_str}")
 
 
 class BinarySource(RelativePositionSource, RelativeFluxSource):
@@ -795,6 +953,44 @@ def model(self      : Source,
 
         # Apply detector if supplied
         return psf if detector is None else detector.apply_detector(psf)
+    
+
+    def summary(self            : Source, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        position = np.array_str(convert_angular(self.position, 'radians', 
+                                            angular_units), precision=sigfigs)
+        separation = convert_angular(self.separation, 'radians', angular_units)
+        position_angle = convert_angular(self.position_angle, 'radians', 
+                                         angular_units)
+        spectrum_str = self.spectrum.summary(angular_units, cartesian_units,
+                                                sigfigs)
+        return (f"{self.name}: A Binary Source with mean flux "
+                f"{self.flux:.{sigfigs}} photons and constrast {self.contrast} "
+                f"offset from the optical axis by {position} "
+                f"{angular_units} with separation {separation:.{sigfigs}} "
+                f"{angular_units}, position angle {position_angle:.{sigfigs}} "
+                f"{angular_units} and spectrum\n  {spectrum_str}")
 
 
 class PointExtendedSource(RelativeFluxSource, ArrayDistribution):
@@ -915,7 +1111,40 @@ def model(self      : Source,
 
         # Apply detector if supplied
         return psf if detector is None else detector.apply_detector(psf)
+    
 
+    def summary(self            : Source, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        position = np.array_str(convert_angular(self.position, 'radians', 
+                                            angular_units), precision=sigfigs)
+        spectrum_str = self.spectrum.summary(angular_units, cartesian_units,
+                                                sigfigs)
+        return (f"{self.name}: A Point and Resolved Source with mean flux "
+                f"{self.flux:.{sigfigs}} photons and constrast "
+                f"{self.contrast:.{sigfigs}} offset from the optical axis by "
+                f"{position} {angular_units} and spectrum\n  "
+                f"{spectrum_str}")
 
 
 class PointAndExtendedSource(RelativeFluxSource, ArrayDistribution):
@@ -1038,3 +1267,36 @@ def model(self      : Source,
 
         # Apply detector if supplied
         return psf if detector is None else detector.apply_detector(psf)
+    
+
+    def summary(self            : Source, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        position = np.array_str(convert_angular(self.position, 'radians', 
+                                            angular_units), precision=sigfigs)
+        spectrum_str = self.spectrum.summary(angular_units, cartesian_units,
+                                                sigfigs)
+        return (f"{self.name}: A Point and Resolved Source with mean flux "
+                f"{self.flux} photons and constrast {self.contrast} offset "
+                f"from the optical axis by {position} "
+                f"{angular_units} and spectrum\n  {spectrum_str}")
diff --git a/dLux/spectrums.py b/dLux/spectrums.py
index 438b69bf..dfbcf944 100644
--- a/dLux/spectrums.py
+++ b/dLux/spectrums.py
@@ -4,6 +4,8 @@
 from equinox import tree_at
 from zodiax import ExtendedBase
 from jax import vmap
+from dLux.utils.units import convert_cartesian
+from dLux.utils.helpers import spectrum_plot
 import dLux
 
 
@@ -21,12 +23,16 @@ class Spectrum(ExtendedBase, ABC):
     ----------
     wavelengths : Array, meters
         The array of wavelengths at which the spectrum is defined.
+    name : str
+        The name of the spectrum.
     """
+    name        : str
     wavelengths : Array
 
 
     def __init__(self        : Spectrum,
-                 wavelengths : Array) -> Spectrum:
+                 wavelengths : Array,
+                 name        : str = 'Spectrum') -> Spectrum:
         """
         Constructor for the Spectrum class.
 
@@ -34,7 +40,10 @@ def __init__(self        : Spectrum,
         ----------
         wavelengths : Array, meters
             The array of wavelengths at which the spectrum is defined.
+        name : str = 'Spectrum'
+            The name of the spectrum.
         """
+        self.name = name
         self.wavelengths = np.asarray(wavelengths, dtype=float)
 
         # Input checking
@@ -72,6 +81,64 @@ def get_weights(self : Spectrum) -> Array: # pragma: no cover
         classes.
         """
         return
+    
+
+    def summary(self            : Spectrum, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return f"{self.name} layer has no summary method yet."
+    
+
+    def display(self            : Spectrum, 
+                figsize         : tuple = (6, 3),
+                dpi             : int = 120,
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> None:
+        """
+        Displays a plot of the wavefront amplitude and opd or phase.
+
+        Parameters
+        ----------
+        figsize : tuple = (6, 3)
+            The size of the figure to display.
+        cmap : str = 'inferno'
+            The colour map to use.
+        dpi : int = 120
+            The resolution of the figure.
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+        """
+        wavelengths = convert_cartesian(self.wavelengths, "meters", 
+                                        cartesian_units)
+        spectrum_plot(wavelengths, self.get_weights(), figsize=figsize,
+                      cartesian_units=cartesian_units, dpi=dpi)
 
 
 class ArraySpectrum(Spectrum):
@@ -85,13 +152,16 @@ class ArraySpectrum(Spectrum):
         The array of wavelengths at which the spectrum is defined.
     weights : Array
         The relative weights of each wavelength.
+    name : str
+        The name of the spectrum.
     """
     weights : Array
 
 
     def __init__(self        : Spectrum,
                  wavelengths : Array,
-                 weights     : Array = None) -> Spectrum:
+                 weights     : Array = None,
+                 name        : str = 'ArraySpectrum') -> Spectrum:
         """
         Constructor for the ArraySpectrum class.
 
@@ -102,8 +172,10 @@ def __init__(self        : Spectrum,
         weights : Array = None
             The relative weights of each wavelength. Defaults to uniform
             spectrum. Weights are automatically normalised to a sum of 1.
+        name : str = 'ArraySpectrum'
+            The name of the spectrum.
         """
-        super().__init__(wavelengths)
+        super().__init__(wavelengths, name)
         weights = np.ones(len(self.wavelengths))/len(wavelengths) \
                                     if weights is None else weights
         weights = np.asarray(weights, dtype=float)
@@ -143,6 +215,33 @@ def normalise(self : Spectrum) -> Spectrum:
         total_power = self.get_weights().sum()
         norm_weights = self.get_weights()/total_power
         return tree_at(lambda spectrum: spectrum.weights, self, norm_weights)
+    
+
+    def summary(self            : Spectrum, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Applies a spectrum with weights stored as an "
+                "array.")
 
 
 class PolynomialSpectrum(Spectrum):
@@ -161,6 +260,8 @@ class PolynomialSpectrum(Spectrum):
         The degree of the polynomial.
     coefficients : Array
         The array of polynomial coefficient values.
+    name : str
+        The name of the spectrum.
     """
     degree       : int # Just a helper
     coefficients : Array
@@ -168,7 +269,8 @@ class PolynomialSpectrum(Spectrum):
 
     def __init__(self         : Spectrum,
                  wavelengths  : Array,
-                 coefficients : Array) -> Spectrum:
+                 coefficients : Array,
+                 name        : str = 'PolynomialSpectrum') -> Spectrum:
         """
         Constructor for the PolynomialSpectrum class.
 
@@ -178,8 +280,10 @@ def __init__(self         : Spectrum,
             The array of wavelengths at which the spectrum is defined.
         coefficients : Array
             The array of polynomial coefficient values.
+        name : str = 'PolynomialSpectrum'
+            The name of the spectrum.
         """
-        super().__init__(wavelengths)
+        super().__init__(wavelengths, name)
 
         self.coefficients = np.asarray(coefficients, dtype=float)
 
@@ -227,6 +331,33 @@ def normalise(self : Spectrum) -> Spectrum:
             The unmodified spectrum object
         """
         return self
+    
+
+    def summary(self            : Spectrum, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Applies a spectrum with weights parametried by "
+                f"a polynomial of degree {self.degree}.")
 
 
 class CombinedSpectrum(ArraySpectrum):
@@ -241,17 +372,20 @@ class CombinedSpectrum(ArraySpectrum):
     Attributes
     ----------
     wavelengths : Array, meters
-        The (2, n) array of wavelengths at which the spectrum is defined.
-        Input can also be a 1d
+        The (nspectra, nwavels) array of wavelengths at which the spectrum is 
+        defined. Input can also be a 1d
     weights : Array
-        The (2, n) relative weights of each wavelength. Defaults to uniform
-        throughput.
+        The (nspectra, nwavels) relative weights of each wavelength. Defaults 
+        to uniform throughput.
+    name : str
+        The name of the spectrum.
     """
 
 
     def __init__(self        : Spectrum,
                  wavelengths : Array,
-                 weights     : Array = None) -> Spectrum:
+                 weights     : Array = None,
+                 name        : str = 'CombinedSpectrum') -> Spectrum:
         """
         Constructor for the CombinedSpectrum class. Expects wavelengths and
         weights to have the same dimensionality, ie (nsources, nwavelengths).
@@ -259,13 +393,16 @@ def __init__(self        : Spectrum,
         Parameters
         ----------
         wavelengths : Array, meters
-            The (2, n) array of wavelengths at which the spectrum is defined.
-            Input can also be a 1d
+            The (nspectra, nwavels) array of wavelengths at which the spectrum 
+            is defined. Input can also be a 1d
         weights : Array (optional)
-            The (2, n) relative weights of each wavelength. Defaults to uniform
-            throughput.
+            The (nspectra, nwavels) relative weights of each wavelength. 
+            Defaults to uniform throughput.
+        name : str = 'CombinedSpectrum'
+            The name of the spectrum.
         """
         super() # Access methods but don't instatiate attributes
+        self.name = str(name)
         self.wavelengths = np.asarray(wavelengths, dtype=float)
 
         # Wavelengths
@@ -312,4 +449,31 @@ def normalise(self : Spectrum) -> Spectrum:
         weights = self.get_weights()
         total_power = weights.sum(1).reshape((len(weights), 1))
         norm_weights = weights/total_power
-        return tree_at(lambda spectrum: spectrum.weights, self, norm_weights)
\ No newline at end of file
+        return tree_at(lambda spectrum: spectrum.weights, self, norm_weights)
+    
+
+    def summary(self            : Spectrum, 
+                angular_units   : str = 'radians', 
+                cartesian_units : str = 'meters', 
+                sigfigs         : int = 4) -> str:
+        """
+        Returns a summary of the class.
+
+        Parameters
+        ----------
+        angular_units : str = 'radians'
+            The angular units to use in the summary. Options are 'radians', 
+            'degrees', 'arcseconds' and 'arcminutes'.
+        cartesian_units : str = 'meters'
+            The cartesian units to use in the summary. Options are 'meters',
+            'millimeters' and 'microns'.
+        sigfigs : int = 4
+            The number of significant figures to use in the summary.
+
+        Returns
+        -------
+        summary : str
+            A summary of the class.
+        """
+        return (f"{self.name}: Applies a two spectrums with weights stored as "
+                "arrays.")
\ No newline at end of file
diff --git a/dLux/utils/coordinates.py b/dLux/utils/coordinates.py
index 5ca4e92c..e718965a 100644
--- a/dLux/utils/coordinates.py
+++ b/dLux/utils/coordinates.py
@@ -1,11 +1,10 @@
 import jax.numpy as np
+from jax.tree_util import tree_map
+from typing import Union
 import dLux
 
 __all__ = ["cartesian_to_polar", "polar_to_cartesian",
-           "get_positions_vector",  "get_pixel_positions",
-           "get_polar_positions",   "get_coordinates_vector",
-           "get_pixel_coordinates", "get_polar_coordinates",
-           "rotate", "translate", "shear", "compress"]
+           "get_pixel_positions", "rotate", "translate", "shear", "compress"]
 
 
 Array = np.ndarray
@@ -53,166 +52,116 @@ def polar_to_cartesian(coordinates : Array) -> Array:
 
 
 ### Positions Calculations ###
-def get_positions_vector(npixels      : int,
-                         pixel_offset : Array = np.array(0.)) -> Array:
-    """
-    Generate the vector pixel positions relative to the center.
-
-    Parameters
-    ----------
-    npixels : int
-        The number of pixels in the vector.
-    pixel_offset : Array = np.array(0.)
-        The number of pixels to offset from the center.
-
-    Returns
-    -------
-    positions : Array
-        The vector of pixel positions.
-    """
-    return np.arange(npixels) - \
-            (npixels - 1) / 2. - pixel_offset
-
-
-def get_pixel_positions(npixels        : int,
-                        x_pixel_offset : Array = np.array(0.),
-                        y_pixel_offset : Array = np.array(0.)) -> Array:
-    """
-    Returns arrays of pixel positions relative to the (central) optical axis.
-
-    Parameters
-    ----------
-    npixels : int
-        The number of pixels along the side of the output array.
-    x_pixel_offset : Array, pixels = np.array(0.)
-        The x offset of the centre of the coordinate system in the
-        square output array.
-    y_pixel_offset : Array, pixels = np.array(0.)
-        The y offset of the centre of the coordinate system in the
-        square output array.
-
-    Returns
-    -------
-    positions : Array
-        The (x, y) array of pixel positions.
-    """
-    x = get_positions_vector(npixels, x_pixel_offset)
-    y = get_positions_vector(npixels, y_pixel_offset)
-    return np.array(np.meshgrid(x, y))
-
-
-def get_polar_positions(npixels        : int,
-                        x_pixel_offset : Array = np.array(0.),
-                        y_pixel_offset : Array = np.array(0.)) -> Array:
-    """
-    Generate the polar positions of each pixel.
-
+def get_pixel_positions(npixels      : Union[int, tuple], 
+                        pixel_scales : Union[tuple, float, None] = None,
+                        offsets      : Union[tuple, float, None] = None,
+                        polar        : bool = False,
+                        indexing     : str = 'xy') -> Array:
+    """
+    Calculates the positions of the pixel centers for the given input. All 
+    coordinates are output in units of meters. 
+    
+    The indexing argument is the same
+    as in numpy.meshgrid., ie:  Giving the string ‘ij’ returns a meshgrid with 
+    matrix indexing, while ‘xy’ returns a meshgrid with Cartesian indexing. In 
+    the 2-D case with inputs of length M and N, the outputs are of shape (N, M) 
+    for ‘xy’ indexing and (M, N) for ‘ij’ indexing. In the 3-D case with inputs 
+    of length M, N and P, outputs are of shape (N, M, P) for ‘xy’ indexing and 
+    (M, N, P) for ‘ij’ indexing. If the output is in polar coordainates, 
+    indexing is set to 'xy' and the input must be 2d
+    
     Parameters
     ----------
-    npixels : int
-        The number of pixels along one edge of the square array.
-    x_pixel_offset : Array = np.array(0.)
-        The x offset of the centre of the coordinate system in the
-        square output array.
-    y_pixel_offset : Array = np.array(0.)
-        The y offset of the centre of the coordinate system in the
-        square output array.
-
+    npixels : Union[int, tuple]
+        The number of pixels in each dimension.
+    pixel_scales : Union[tuple, float, None] = None
+        The pixel scales in each dimension. If a tuple, the length
+        of the tuple must match the number of dimensions. If a float, the same
+        scale is applied to all dimensions. If None, the scale is set to 1.
+    offsets : Union[tuple, float, None] = None
+        The offset of the pixel centers in each dimension. If a tuple, the 
+        length of the tuple must match the number of dimensions. If a float, 
+        the same offset is applied to all dimensions. If None, the offset is 
+        set to 0.
+    polar : bool = False
+        If True, the output is in polar coordinates. If False, the output is in
+        cartesian coordinates. Default is False.
+    indexing : str = 'xy'
+        The indexing of the output. Default is 'xy'. See numpy.meshgrid for more
+        details.
+    
     Returns
     -------
     positions : Array
-        The (r, phi) array of the radial coordinates.
-    """
-    positions = get_pixel_positions(npixels, x_pixel_offset, y_pixel_offset)
-    return cartesian_to_polar(positions)
-
-
-### Coordinate Calculations ###
-def get_coordinates_vector(npixels     : int,
-                           pixel_scale : Array,
-                           offset      : Array = np.array(0.)) -> Array:
-    """
-    Generate the vector pixel coordinates relative to the center.
-
-    Parameters
-    ----------
-    npixels : int
-        The number of pixels in the vector.
-    offset : Array, meters = np.array(0.)
-        The offset from the center.
-
-    Returns
-    -------
-    x : Array
-        The vector of pixel coordinates.
-    """
-    return pixel_scale * get_positions_vector(npixels, pixel_scale * offset)
-
-
-def get_pixel_coordinates(npixels     : int,
-                          pixel_scale : Array,
-                          x_offset    : Array = np.array(0.),
-                          y_offset    : Array = np.array(0.)) -> Array:
-    """
-    Returns the physical (x, y) coordinate array of pixel
-    positions relative to the (central) optical axis.
-
-    Parameters
-    ----------
-    npixels : int
-        The number of pixels along the side of the output array.
-    pixel_scale : Array, meters
-        The physical size of each pixel
-    x_offset : Array = np.array(0.), meters
-        The x offset of the centre of the coordinate system in the
-        square output array.
-    y_offset : Array = np.array(0.), meters
-        The y offset of the centre of the coordinate system in the
-        square output array.
-
-    Returns
-    -------
-    pixel_positions : Array
-        The (x, y) coordinate arrays in the square output array with the
-        correct offsets.
-    """
-    return pixel_scale * get_pixel_positions(npixels,
-                                             pixel_scale * x_offset,
-                                             pixel_scale * y_offset)
-
-
-def get_polar_coordinates(npixels     : int,
-                          pixel_scale : Array,
-                          x_offset    : Array = np.array(0.),
-                          y_offset    : Array = np.array(0.)) -> Array:
-    """
-    Generate the (r, phi) polar coordinates of each pixel.
-
-    Parameters
-    ----------
-    npixels : int
-        The number of pixels along one edge of the square array.
-    pixel_scale : Array, meters
-        The physical size of each pixel
-    x_offset : Array = np.array(0.).
-        The x offset of the centre of the coordinate system in the
-        square output array.
-    y_offset : Array = np.array(0.)
-        The y offset of the centre of the coordinate system in the
-        square output array.
-
-    Returns
-    -------
-    positions : Array
-        The (r, phi) coordinate arrays in the square output array with the
-        correct offsets
-    """
-    pixel_scaler = np.array([pixel_scale, 1.])[:, None, None]
-    return pixel_scaler * get_polar_positions(npixels,
-                                              x_offset / pixel_scale,
-                                              y_offset / pixel_scale)
-
-
+        The positions of the pixel centers in the given dimensions.
+    """
+    if indexing not in ['xy', 'ij']:
+        raise ValueError("indexing must be either 'xy' or 'ij'.")
+    
+    if polar and indexing == 'ij':
+        indexing = 'xy'
+
+    # Turn inputs into tuples
+    if isinstance(npixels, int):
+        npixels = (npixels,)
+
+        if offsets is None:
+            offsets = (0.,)
+        elif not isinstance(offsets, (float, Array)):
+            raise ValueError("offset must be an be a float or Array if npixels "
+                             "is an int.")
+        else:
+            offsets = (offsets,)
+
+        if pixel_scales is None:
+            pixel_scales = (1.,)
+        elif not isinstance(pixel_scales, (float, Array)):
+            raise ValueError("pixel_scales must be an be a float if npixels is an int.")
+        else:
+            pixel_scales = (pixel_scales,)
+        
+    # Check input 
+    else:
+        if offsets is None:
+            offsets = tuple([0.]*len(npixels))
+        elif not isinstance(offsets, tuple):
+            raise ValueError("offset must be an be a float or Array if npixels "
+                             "is an int.")
+        else:
+            if len(offsets) != len(npixels):
+                raise ValueError("offset must have the same length as npixels.")
+            
+        if pixel_scales is None:
+            pixel_scales = tuple([1.]*len(npixels))
+        elif isinstance(pixel_scales, float):
+            pixel_scales = tuple([pixel_scales]*len(npixels))
+        elif not isinstance(pixel_scales, tuple):
+            raise ValueError("pixel_scales must be a tuple if npixels is a tuple.")
+        else:
+            if len(pixel_scales) != len(npixels):
+                raise ValueError("pixel_scales must have the same length as npixels.")
+    
+    def pixel_fn(n, offset, scale):
+        pix = np.arange(n) - (n - 1) / 2.
+        pix *= scale
+        pix -= offset
+        return pix
+    
+    pixels = tree_map(pixel_fn, npixels, offsets, pixel_scales)
+
+    # ouput (x, y) for 2d, else in order
+    positions = np.array(np.meshgrid(*pixels, indexing=indexing))
+
+    if polar:
+        if len(npixels) != 2:
+            raise ValueError("polar coordinates are only defined for 2D arrays.")
+        return cartesian_to_polar(positions)
+
+    # Squeeze for empty axis removal with 1d
+    return np.squeeze(positions)
+
+
+### Coordinates Transformations ###
 def rotate(coordinates: Array, rotation: Array) -> Array:
     """
     Rotate the coordinate system by a pre-specified amount.
diff --git a/dLux/utils/gradient_energy.py b/dLux/utils/gradient_energy.py
index 7f9b5978..74b5eba7 100644
--- a/dLux/utils/gradient_energy.py
+++ b/dLux/utils/gradient_energy.py
@@ -42,7 +42,8 @@ def get_RGE(array : Array, epsilon : float = 1e-8) -> Array:
     array : Array
         The array of radial gradient energies.
     """
-    positions = dLux.utils.coordinates.get_pixel_positions(array.shape[0])
+    npix = array.shape[0]
+    positions = dLux.utils.coordinates.get_pixel_positions((npix, npix))
     grads_vec = np.gradient(array)
 
     xnorm = positions[1]*grads_vec[0]
@@ -67,8 +68,10 @@ def get_RWGE(array : Array, epsilon : float = 1e-8) -> Array:
     array : Array
         The array of radial radially weighted energies.
     """
-    positions = dLux.utils.coordinates.get_pixel_positions(array.shape[0])
-    radii = dLux.utils.coordinates.get_polar_positions(array.shape[0])[0]
+    npix = array.shape[0]
+    positions = dLux.utils.coordinates.get_pixel_positions((npix, npix))
+    radii = dLux.utils.coordinates.get_pixel_positions((npix, npix), 
+                                                        polar=True)[0]
     radii_norm = positions/(radii + epsilon)
     grads_vec = np.gradient(array)
 
@@ -98,5 +101,6 @@ def get_radial_mask(npixels : int,
     mask: Array
         A mask with the the values below rmin and above rmax masked out.
     """
-    radii = dLux.utils.coordinates.get_polar_positions(npixels)[0]
+    radii = dLux.utils.coordinates.get_pixel_positions((npixels, npixels), 
+                                                        polar=True)[0]
     return np.asarray((radii < rmax) & (radii > rmin), dtype=float)
\ No newline at end of file
diff --git a/dLux/utils/helpers.py b/dLux/utils/helpers.py
index 921572b8..b2a7af51 100644
--- a/dLux/utils/helpers.py
+++ b/dLux/utils/helpers.py
@@ -1,8 +1,10 @@
 import jax.numpy as np
 from collections import OrderedDict
+import matplotlib.pyplot as plt
 
 
-__all__ = ["list_to_dictionary"]
+__all__ = ["list_to_dictionary", "single_image_plot", "two_image_plot",
+           "spectrum_plot"]
 
 
 Array = np.ndarray
@@ -69,3 +71,129 @@ def list_to_dictionary(list_in : list, ordered : bool = True) -> dict:
         # dict_out[names[i]] = list_in[i].set_name(names[i])
 
     return dict_out
+
+
+def single_image_plot(array       : Array, 
+                      figsize     : tuple = (5, 4), 
+                      title       : str   = "Array",
+                      cbar_label  : str   = None,
+                      cmap        : str   = "inferno",
+                      bound       : float = None,
+                      dpi         : int   = 120) -> None:
+    """
+    Plots a  single image.
+
+    Parameters
+    ----------
+    array : Array
+        The first array to plot.
+    figsize : tuple = (5, 4)
+        The size of the figure to display.
+    title : str = "Array"
+        The title of the array.
+    cmap : str = "inferno"
+        The colour map to use.
+    bound : float = None
+        The bound of the colour map.
+    dpi : int = 120
+        The resolution of the figure.
+    """
+    plt.figure(figsize=figsize)
+    plt.title(title)
+    if bound is not None:
+        array = array % bound
+    plt.imshow(array, cmap=cmap)
+    cbar = plt.colorbar()
+    if cbar_label is not None:
+        cbar.set_label(cbar_label)
+    plt.show()
+
+
+def two_image_plot(array1      : Array, 
+                   array2      : Array, 
+                   figsize     : tuple = (10, 4), 
+                   titles      : tuple = ("Array 1", "Array 2"),
+                   cbar_labels : tuple = (None, None),
+                   cmaps       : tuple = ("inferno", "inferno"),
+                   bounds      : tuple = (None, None),
+                   dpi         : int   = 120) -> None:
+    """
+    Plots two images side by side.
+
+    Parameters
+    ----------
+    array1 : Array
+        The first array to plot.
+    array2 : Array
+        The second array to plot.
+    figsize : tuple = (10, 4)
+        The size of the figure to display.
+    titles : tuple = ("Array 1", "Array 2")
+        The titles of the arrays.
+    cmaps : tuple = ("inferno", "inferno")
+        The colour maps to use.
+    bounds : tuple = (None, None)
+        The bounds of the colour maps.
+    dpi : int = 120
+        The resolution of the figure.
+    """
+    plt.figure(figsize=figsize)
+    plt.subplot(1, 2, 1)
+    plt.title(titles[0])
+    if bounds[0] is not None:
+        array1 = array1 % bounds[0]
+    plt.imshow(array1, cmap=cmaps[0])
+    cbar = plt.colorbar()
+    if cbar_labels[0] is not None:
+        cbar.set_label(cbar_labels[0])
+
+    plt.subplot(1, 2, 2)
+    plt.title(titles[1])
+    if bounds[1] is not None:
+        array2 = array2 % bounds[1]
+    plt.imshow(array2, cmap=cmaps[1])
+    cbar = plt.colorbar()
+    if cbar_labels[1] is not None:
+        cbar.set_label(cbar_labels[1])
+    plt.show()
+
+
+def spectrum_plot(wavelengths     : Array, 
+                  weights         : Array, 
+                  figsize         : tuple = (6, 3),
+                  labels          : tuple = None, 
+                  cartesian_units : str = 'meters',
+                  dpi             : int = 120) -> None:
+    """
+    Plots a spectrum based on wavelgths and weights.
+
+    Parameters
+    ----------
+    wavelengths : Array, meters
+        The wavelengths of the spectrum.
+    weights : Array
+        The weights of the spectrum.
+    figsize : tuple = (6, 3)
+        The size of the figure to display.
+    labels : tuple = None
+        The labels of the spectra.
+    cartesian_units : str = 'meters'
+        The units of the wavelengths.
+    dpi : int = 120
+        The resolution of the figure.
+    """
+    nspectra = 1 if wavelengths.ndim == 1 else len(wavelengths)
+    if labels is None:
+        labels = ['Spectrum {}'.format(i) for i in range(nspectra)]
+    
+    plt.figure(figsize=figsize, dpi=dpi)
+    if nspectra == 1:
+        plt.scatter(wavelengths, weights, label=labels[0])
+    else:
+        for i in range(nspectra):
+            plt.scatter(wavelengths[i], weights[i], label=labels[i])
+    plt.legend()
+    plt.xlabel(f"Wavelengths {cartesian_units}")
+    plt.ylabel("Weights")
+    plt.ylim(-0.01)
+    plt.show()
\ No newline at end of file
diff --git a/dLux/utils/interpolation.py b/dLux/utils/interpolation.py
index 1d64225b..6d448117 100644
--- a/dLux/utils/interpolation.py
+++ b/dLux/utils/interpolation.py
@@ -2,7 +2,8 @@
 from jax import vmap
 from jax.scipy.ndimage import map_coordinates
 from functools import partial
-from dLux.utils.coordinates import get_polar_positions, polar_to_cartesian
+import dLux.utils.coordinates as c
+
 
 
 __all__ = ["scale_array", "generate_coordinates", "interpolate_field",
@@ -175,6 +176,7 @@ def rotate_field(field          : Array,
                  angle          : Array,
                  real_imaginary : bool = False,
                  fourier        : bool = False,
+                 order          : int  = 1,
                  padding        : int  = 2) -> Array:
     """
     Paraxially rotates a wavefront field (either in ampltude and phase, or
@@ -195,6 +197,8 @@ def rotate_field(field          : Array,
     fourier : bool = False
         Should the fourier rotation method be used (True), or regular
         interpolation method be used (False).
+    order : int = 2
+        The interpolation order to use. Must be 0, 1, or 3.
     padding : int = 2
         The amount of fourier padding to use. Only applies if fourier is True.
 
@@ -208,7 +212,8 @@ def rotate_field(field          : Array,
         padded_rotate = partial(fourier_rotate, padding=padding)
         rotator = vmap(padded_rotate, in_axes=(0, None))
     else:
-        rotator = vmap(rotate, in_axes=(0, None))
+        order_rotate = partial(rotate, order=order)
+        rotator = vmap(order_rotate, in_axes=(0, None))
 
     # Rotate
     rotated_field = rotator(field, angle)
@@ -223,7 +228,7 @@ def rotate_field(field          : Array,
     return np.array([amplitude, phase])
 
 
-def rotate(array : Array, angle : Array) -> Array:
+def rotate(array : Array, angle : Array, order : int = 1) -> Array:
     """
     Rotates an array by the angle, using linear interpolation.
 
@@ -239,16 +244,16 @@ def rotate(array : Array, angle : Array) -> Array:
     array : Array
         The rotated array.
     """
+    if order not in (0, 1, 3):
+        raise ValueError("Order must be 0, 1, or 3.")
     # Get coordinates
     npixels = array.shape[0]
     centre = (npixels - 1) / 2
-    coordinates = get_polar_positions(npixels)
-    coordinates += np.array([0., angle])[:, None, None]
-    coordinates_rotated = np.roll(polar_to_cartesian(coordinates) + centre,
-                                  shift=1, axis=0)
+    coordinates = c.get_pixel_positions((npixels, npixels), indexing='ij')
+    coordinates_rotated = c.rotate(coordinates, angle) + centre
 
     # Interpolate
-    return map_coordinates(array, coordinates_rotated, order=1)
+    return map_coordinates(array, coordinates_rotated, order=order)
 
 
 def fourier_rotate(array   : Array,
diff --git a/dLux/utils/units.py b/dLux/utils/units.py
index 283e253e..fe16c48c 100644
--- a/dLux/utils/units.py
+++ b/dLux/utils/units.py
@@ -1,15 +1,134 @@
 import jax.numpy as np
 
-__all__ = ["radians_to_arcseconds", "arcseconds_to_radians"]
-
+__all__ = ["convert_cartesian", "convert_angular",
+           "radians_to_arcseconds", "radians_to_degrees", 
+           "radians_to_arcminutes", "degrees_to_radians", 
+           "degrees_to_arcminutes", "degrees_to_arcseconds", 
+           "arcminutes_to_radians", "arcminutes_to_degrees", 
+           "arcminutes_to_arcseconds", "arcseconds_to_radians", 
+           "arcseconds_to_degrees", "arcseconds_to_arcminutes",
+           "r2s", "r2d", "r2m", "d2r", "d2m", "d2s", 
+           "m2r", "m2d", "m2s", "s2r", "s2d", "s2m"]
 
 Array = np.ndarray
 
 
+# General conversion classes:
+def convert_cartesian(values : Array, 
+                      input  : str = 'meters', 
+                      output : str = 'meters') -> Array:
+    """
+    Converts the input values from one unit to another.
+
+    Parameters
+    ----------
+    values : Array
+        The input values to be converted.
+    input : str = 'meters'
+        The input units. Must be one of 'meters', 'millimeters', or 'microns'.
+    output : str = 'meters'
+        The output units. Must be one of 'meters', 'millimeters', or 'microns'.
+    
+    Returns
+    -------
+    values : Array
+        The input values converted into the output units.
+    """
+    if input not in ('meters', 'millimeters', 'microns'):
+        raise ValueError("input must be 'meters', 'millimeters', or 'microns'.")
+    if output not in ('meters', 'millimeters', 'microns'):
+        raise ValueError("output must be 'meters', 'millimeters', or "
+                         "'microns'.")
+    
+    if input == output:
+        factor = 1
+    elif input == 'meters':
+        if output == 'millimeters':
+            factor = 1e-3
+        elif output == 'microns':
+            factor = 1e-6
+    elif input == 'millimeter':
+        if output == 'meters':
+            factor = 1e3
+        elif output == 'microns':
+            factor = 1e-3
+    elif input == 'microns':
+        if output == 'meters':
+            factor = 1e6
+        elif output == 'millimeters':
+            factor = 1e3
+    return values * factor
+
+
+def convert_angular(values : Array,
+                    input  : str = 'radians',
+                    output : str = 'radians') -> Array:
+    """
+    Converts the input values from one unit to another.
+
+    Parameters
+    ----------
+    values : Array
+        The input values to be converted.
+    input : str = 'radians'
+        The input units. Must be one of 'radians', 'degrees', 'arcseconds', or
+        'arcminutes'.
+    output : str = 'radians'
+        The output units. Must be one of 'radians', 'degrees', 'arcseconds', or
+        'arcminutes'.
+    
+    Returns
+    -------
+    values : Array
+        The input values converted into the output units.
+    """
+    if input not in ('radians', 'degrees', 'arcseconds', 'arcminutes'):
+        raise ValueError(f"input must be one of 'radians', 'degrees', "
+                         f"'arcseconds' or 'arcminutes'.")
+    if output not in ('radians', 'degrees', 'arcseconds', 'arcminutes'):
+        raise ValueError(f"output must be one of 'radians', 'degrees', "
+                         f"'arcseconds' or 'arcminutes'.")
+
+    if input == output:
+        factor = 1
+    elif input == 'radians':
+        if output == 'degrees':
+            factor = r2d(1)
+        elif output == 'arcminutes':
+            factor = r2m(1)
+        elif output == 'arcseconds':
+            factor = r2s(1)
+    elif input == 'degrees':
+        if output == 'radians':
+            factor = d2r(1)
+        elif output == 'arcminutes':
+            factor = d2m(1)
+        elif output == 'arcseconds':
+            factor = d2s(1)
+    elif input == 'arcminutes':
+        if output == 'radians':
+            factor = m2r(1)
+        elif output == 'degrees':
+            factor = m2d(1)
+        elif output == 'arcseconds':
+            factor = m2s(1)
+    elif input == 'arcseconds':
+        if output == 'radians':
+            factor = s2r(1)
+        elif output == 'degrees':
+            factor = s2d(1)
+        elif output == 'arcminutes':
+            factor = s2m(1)
+    return values * factor
+
+
+# Radians to:
 def radians_to_arcseconds(values : Array) -> Array:
     """
     Converts the inputs values from radians to arcseconds.
 
+    Can also be imported as r2s.
+
     Parameters
     ----------
     values : Array, radians
@@ -23,10 +142,166 @@ def radians_to_arcseconds(values : Array) -> Array:
     return values * 3600 * 180 / np.pi
 
 
+def radians_to_degrees(values : Array) -> Array:
+    """
+    Converts the inputs values from radians to degrees.
+
+    Can also be imported as r2d.
+
+    Parameters
+    ----------
+    values : Array, radians
+        The input values in units of radians to be converted into degrees.
+
+    Returns
+    -------
+    values : Array, degrees
+        The input values converted into degrees.
+    """
+    return values * 180 / np.pi
+
+
+def radians_to_arcminutes(values : Array) -> Array:
+    """
+    Converts the inputs values from radians to arcminutes.
+
+    Can also be imported as r2m.
+
+    Parameters
+    ----------
+    values : Array, radians
+        The input values in units of radians to be converted into arcminutes.
+
+    Returns
+    -------
+    values : Array, arcminutes
+        The input values converted into arcminutes.
+    """
+    return values * 60 * 180 / np.pi
+
+
+# Degrees to:
+def degrees_to_radians(values : Array) -> Array:
+    """
+    Converts the inputs values from degrees to radians.
+
+    Can also be imported as d2r.
+
+    Parameters
+    ----------
+    values : Array, degrees
+        The input values in units of degrees to be converted into radians.
+
+    Returns
+    -------
+    values : Array, radians
+        The input values converted into radians.
+    """
+    return values * np.pi / 180
+
+
+def degrees_to_arcminutes(values : Array) -> Array:
+    """
+    Converts the inputs values from degrees to arcminutes.
+
+    Can also be imported as d2m.
+
+    Parameters
+    ----------
+    values : Array, degrees
+        The input values in units of degrees to be converted into arcminutes.
+
+    Returns
+    -------
+    values : Array, arcminutes
+        The input values converted into arcminutes.
+    """
+    return values * 60
+
+
+def degrees_to_arcseconds(values : Array) -> Array:
+    """
+    Converts the inputs values from degrees to arcseconds.
+
+    Can also be imported as d2s.
+
+    Parameters
+    ----------
+    values : Array, degrees
+        The input values in units of degrees to be converted into arcseconds.
+
+    Returns
+    -------
+    values : Array, arcseconds
+        The input values converted into arcseconds.
+    """
+    return values * 3600
+
+# Arcminutes to:
+def arcminutes_to_radians(values : Array) -> Array:
+    """
+    Converts the inputs values from arcminutes to radians.
+
+    Can also be imported as m2r.
+
+    Parameters
+    ----------
+    values : Array, arcminutes
+        The input values in units of arcminutes to be converted into radians.
+
+    Returns
+    -------
+    values : Array, radians
+        The input values converted into radians.
+    """
+    return values * np.pi / (60 * 180)
+
+
+def arcminutes_to_degrees(values : Array) -> Array:
+    """
+    Converts the inputs values from arcminutes to degrees.
+
+    Can also be imported as m2d.
+
+    Parameters
+    ----------
+    values : Array, arcminutes
+        The input values in units of arcminutes to be converted into degrees.
+
+    Returns
+    -------
+    values : Array, degrees
+        The input values converted into degrees.
+    """
+    return values / 60
+
+
+def arcminutes_to_arcseconds(values : Array) -> Array:
+    """
+    Converts the inputs values from arcminutes to arcseconds.
+
+    Can also be imported as m2s.
+
+    Parameters
+    ----------
+    values : Array, arcminutes
+        The input values in units of arcminutes to be converted into arcseconds.
+
+    Returns
+    -------
+    values : Array, arcseconds
+        The input values converted into arcseconds.
+    """
+    return values * 60
+
+
+# Arcseconds to:
 def arcseconds_to_radians(values : Array) -> Array:
     """
     Converts the inputs values from arcseconds to radians.
 
+    Can also be imported as s2r.
+
     Parameters
     ----------
     values : Array, arcseconds
@@ -40,4 +315,53 @@ def arcseconds_to_radians(values : Array) -> Array:
     return values * np.pi / (3600 * 180)
 
 
+def arcseconds_to_degrees(values : Array) -> Array:
+    """
+    Converts the inputs values from arcseconds to degrees.
+
+    Can also be imported as s2d.
+
+    Parameters
+    ----------
+    values : Array, arcseconds
+        The input values in units of arcseconds to be converted into degrees.
+
+    Returns
+    -------
+    values : Array, degrees
+        The input values converted into degrees.
+    """
+    return values / 3600
+
+
+def arcseconds_to_arcminutes(values : Array) -> Array:
+    """
+    Converts the inputs values from arcseconds to arcminutes.
+
+    Can also be imported as s2m.
+
+    Parameters
+    ----------
+    values : Array, arcseconds
+        The input values in units of arcseconds to be converted into arcminutes.
+
+    Returns
+    -------
+    values : Array, arcminutes
+        The input values converted into arcminutes.
+    """
+    return values / 60
 
+# Alias to simpler names
+r2s = radians_to_arcseconds
+r2d = radians_to_degrees
+r2m = radians_to_arcminutes
+d2r = degrees_to_radians
+d2m = degrees_to_arcminutes
+d2s = degrees_to_arcseconds
+m2r = arcminutes_to_radians
+m2d = arcminutes_to_degrees
+m2s = arcminutes_to_arcseconds
+s2r = arcseconds_to_radians
+s2d = arcseconds_to_degrees
+s2m = arcseconds_to_arcminutes
\ No newline at end of file
diff --git a/dLux/wavefronts.py b/dLux/wavefronts.py
index da431a71..54cf7831 100644
--- a/dLux/wavefronts.py
+++ b/dLux/wavefronts.py
@@ -5,7 +5,7 @@
 from zodiax import ExtendedBase
 from enum import IntEnum
 from abc import ABC
-from dLux.utils.coordinates import get_pixel_coordinates
+from dLux.utils.coordinates import get_pixel_positions
 from dLux.utils.interpolation import interpolate_field, rotate_field
 import dLux
 
@@ -214,7 +214,9 @@ def pixel_coordinates(self : Wavefront) -> Array:
             The coordinates of the centers of each pixel representing the
             wavefront.
         """
-        return get_pixel_coordinates(self.npixels, self.pixel_scale)
+        # return get_pixel_coordinates(self.npixels, self.pixel_scale)
+        return get_pixel_positions((self.npixels, self.npixels), 
+                                   (self.pixel_scale, self.pixel_scale))
 
 
     @property
@@ -571,6 +573,7 @@ def rotate(self           : Wavefront,
                angle          : Array,
                real_imaginary : bool = False,
                fourier        : bool = False,
+               order          : int  = 1,
                padding        : int  = 2) -> Wavefront:
         """
         Performs a paraxial rotation on the wavefront, determined by the
@@ -591,6 +594,8 @@ def rotate(self           : Wavefront,
         fourier : bool = False
             Should the fourier rotation method be used (True), or regular
             interpolation method be used (False).
+        order : int = 2
+            The interpolation order to use. Must be 0, 1, or 3.
         padding : int = 2
             The amount of fourier padding to use. Only applies if fourier is
             True.
@@ -611,11 +616,11 @@ def rotate(self           : Wavefront,
         if field.shape[1] == 1:
             new_amplitude, new_phase = \
             rotate_field(field[:, 0], angle, fourier=fourier,
-                         real_imaginary=real_imaginary)[:, None, :, :]
+                    real_imaginary=real_imaginary, order=order)[:, None, :, :]
         else:
             rotator = vmap(rotate_field, in_axes=(1, None))
             new_amplitude, new_phase = rotator(field, angle, fourier=fourier,
-                                               real_imaginary=real_imaginary)
+                                    real_imaginary=real_imaginary, order=order)
 
         # Update parameters
         return tree_at(lambda wavefront: (wavefront.amplitude, wavefront.phase),
diff --git a/docs/utils/coordinates/get_coordinates_vector.md b/docs/utils/coordinates/get_coordinates_vector.md
deleted file mode 100644
index 2202a1e2..00000000
--- a/docs/utils/coordinates/get_coordinates_vector.md
+++ /dev/null
@@ -1,1 +0,0 @@
-::: dLux.utils.coordinates.get_coordinates_vector
diff --git a/docs/utils/coordinates/get_pixel_coordinates.md b/docs/utils/coordinates/get_pixel_coordinates.md
deleted file mode 100644
index 52e69e72..00000000
--- a/docs/utils/coordinates/get_pixel_coordinates.md
+++ /dev/null
@@ -1,1 +0,0 @@
-::: dLux.utils.coordinates.get_pixel_coordinates
diff --git a/docs/utils/coordinates/get_polar_coordinates.md b/docs/utils/coordinates/get_polar_coordinates.md
deleted file mode 100644
index 1a23696b..00000000
--- a/docs/utils/coordinates/get_polar_coordinates.md
+++ /dev/null
@@ -1,1 +0,0 @@
-::: dLux.utils.coordinates.get_polar_coordinates
diff --git a/docs/utils/coordinates/get_polar_positions.md b/docs/utils/coordinates/get_polar_positions.md
deleted file mode 100644
index a478cad4..00000000
--- a/docs/utils/coordinates/get_polar_positions.md
+++ /dev/null
@@ -1,1 +0,0 @@
-::: dLux.utils.coordinates.get_polar_positions
diff --git a/docs/utils/coordinates/get_positions_vector.md b/docs/utils/coordinates/get_positions_vector.md
deleted file mode 100644
index 9f495d38..00000000
--- a/docs/utils/coordinates/get_positions_vector.md
+++ /dev/null
@@ -1,1 +0,0 @@
-::: dLux.utils.coordinates.get_positions_vector
diff --git a/docs/utils/helpers/single_image_plot.md b/docs/utils/helpers/single_image_plot.md
new file mode 100644
index 00000000..b66e5ee5
--- /dev/null
+++ b/docs/utils/helpers/single_image_plot.md
@@ -0,0 +1,1 @@
+::: dLux.utils.helpers.single_image_plot
\ No newline at end of file
diff --git a/docs/utils/helpers/spectrum_plot.md b/docs/utils/helpers/spectrum_plot.md
new file mode 100644
index 00000000..fae5b1ff
--- /dev/null
+++ b/docs/utils/helpers/spectrum_plot.md
@@ -0,0 +1,1 @@
+::: dLux.utils.helpers.spectrum_plot
\ No newline at end of file
diff --git a/docs/utils/helpers/two_image_plot.md b/docs/utils/helpers/two_image_plot.md
new file mode 100644
index 00000000..ded13c40
--- /dev/null
+++ b/docs/utils/helpers/two_image_plot.md
@@ -0,0 +1,1 @@
+::: dLux.utils.helpers.two_image_plot
\ No newline at end of file
diff --git a/docs/utils/units/arcminutes_to_arcseconds.md b/docs/utils/units/arcminutes_to_arcseconds.md
new file mode 100644
index 00000000..77c5041d
--- /dev/null
+++ b/docs/utils/units/arcminutes_to_arcseconds.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.arcminutes_to_arcseconds
\ No newline at end of file
diff --git a/docs/utils/units/arcminutes_to_degrees.md b/docs/utils/units/arcminutes_to_degrees.md
new file mode 100644
index 00000000..0e974457
--- /dev/null
+++ b/docs/utils/units/arcminutes_to_degrees.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.arcminutes_to_degrees
\ No newline at end of file
diff --git a/docs/utils/units/arcminutes_to_radians.md b/docs/utils/units/arcminutes_to_radians.md
new file mode 100644
index 00000000..8c388fc5
--- /dev/null
+++ b/docs/utils/units/arcminutes_to_radians.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.arcminutes_to_radians
\ No newline at end of file
diff --git a/docs/utils/units/arcseconds_to_arcminutes.md b/docs/utils/units/arcseconds_to_arcminutes.md
new file mode 100644
index 00000000..ad43d12f
--- /dev/null
+++ b/docs/utils/units/arcseconds_to_arcminutes.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.arcseconds_to_arcminutes
\ No newline at end of file
diff --git a/docs/utils/units/arcseconds_to_degrees.md b/docs/utils/units/arcseconds_to_degrees.md
new file mode 100644
index 00000000..96bb9cf0
--- /dev/null
+++ b/docs/utils/units/arcseconds_to_degrees.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.arcseconds_to_degrees
\ No newline at end of file
diff --git a/docs/utils/units/convert_angular.md b/docs/utils/units/convert_angular.md
new file mode 100644
index 00000000..b3dd6962
--- /dev/null
+++ b/docs/utils/units/convert_angular.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.convert_angular
\ No newline at end of file
diff --git a/docs/utils/units/convert_cartesian.md b/docs/utils/units/convert_cartesian.md
new file mode 100644
index 00000000..77618c8f
--- /dev/null
+++ b/docs/utils/units/convert_cartesian.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.convert_cartesian
\ No newline at end of file
diff --git a/docs/utils/units/degrees_to_arcminutes.md b/docs/utils/units/degrees_to_arcminutes.md
new file mode 100644
index 00000000..72635aea
--- /dev/null
+++ b/docs/utils/units/degrees_to_arcminutes.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.degrees_to_arcminutes
\ No newline at end of file
diff --git a/docs/utils/units/degrees_to_arcseconds.md b/docs/utils/units/degrees_to_arcseconds.md
new file mode 100644
index 00000000..b3148438
--- /dev/null
+++ b/docs/utils/units/degrees_to_arcseconds.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.degrees_to_arcseconds
\ No newline at end of file
diff --git a/docs/utils/units/degrees_to_radians.md b/docs/utils/units/degrees_to_radians.md
new file mode 100644
index 00000000..8e9a7fdb
--- /dev/null
+++ b/docs/utils/units/degrees_to_radians.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.degrees_to_radians
\ No newline at end of file
diff --git a/docs/utils/units/radians_to_arcminutes.md b/docs/utils/units/radians_to_arcminutes.md
new file mode 100644
index 00000000..a8fc7294
--- /dev/null
+++ b/docs/utils/units/radians_to_arcminutes.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.radians_to_arcminutes
\ No newline at end of file
diff --git a/docs/utils/units/radians_to_degrees.md b/docs/utils/units/radians_to_degrees.md
new file mode 100644
index 00000000..1ec82ed0
--- /dev/null
+++ b/docs/utils/units/radians_to_degrees.md
@@ -0,0 +1,1 @@
+::: dLux.utils.units.radians_to_degrees
\ No newline at end of file
diff --git a/mkdocs.yml b/mkdocs.yml
index bd9f1b08..f9f040ab 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -133,8 +133,20 @@ nav:
             - simple optical system: utils/models/simple_optical_system.md
             - toliman:               utils/models/toliman.md
         - Units:
-            - radians to arcseconds: utils/units/radians_to_arcseconds.md
-            - arcseconds to radians: utils/units/arcseconds_to_radians.md
+            - convert cartesian:        utils/units/convert_cartesian.md
+            - convert angular:          utils/units/convert_angular.md
+            - radians to arcseconds:    utils/units/radians_to_arcseconds.md
+            - radians to degrees:       utils/units/radians_to_degrees.md
+            - radians to arcminutes:    utils/units/radians_to_arcminutes.md
+            - degrees to radians:       utils/units/degrees_to_radians.md
+            - degrees to arcminutes:    utils/units/degrees_to_arcminutes.md
+            - degrees to arcseconds:    utils/units/degrees_to_arcseconds.md
+            - arcminutes to radians:    utils/units/arcminutes_to_radians.md
+            - arcminutes to degrees:    utils/units/arcminutes_to_degrees.md
+            - arcminutes to arcseconds: utils/units/arcminutes_to_arcseconds.md
+            - arcseconds to radians:    utils/units/arcseconds_to_radians.md
+            - arcseconds to degrees:    utils/units/arcseconds_to_degrees.md
+            - arcseconds to arcminutes: utils/units/arcseconds_to_arcminutes.md
         - Math:
             - factorial: utils/math/factorial.md
         - Coordinates:
@@ -144,12 +156,7 @@ nav:
             - shear:                  utils/coordinates/shear.md
             - cartesian to polar:     utils/coordinates/cartesian_to_polar.md
             - polar to cartesian:     utils/coordinates/polar_to_cartesian.md
-            - get coordinates vector: utils/coordinates/get_coordinates_vector.md
-            - get positions vector:   utils/coordinates/get_positions_vector.md
             - get pixel positions:    utils/coordinates/get_pixel_positions.md
-            - get pixel coordinates:  utils/coordinates/get_pixel_coordinates.md
-            - get polar positions:    utils/coordinates/get_polar_positions.md
-            - get polar coordinates:  utils/coordinates/get_polar_coordinates.md
         - Interpolation:
             - scale array:          utils/interpolation/scale_array.md
             - generate coordinates: utils/interpolation/generate_coordinates.md
@@ -178,3 +185,4 @@ nav:
             - get radial mask: utils/gradient_energy/get_radial_mask.md
         - Helpers:
             - list to dictionary: utils/helpers/list_to_dictionary.md
+            - two image plot:     utils/helpers/two_image_plot.md
diff --git a/setup.py b/setup.py
index 215b3d2c..022dc9d4 100644
--- a/setup.py
+++ b/setup.py
@@ -26,7 +26,7 @@ def find_version(*file_paths):
 
 
 setuptools.setup(
-    python_requires='>=3.7,<3.11',
+    python_requires='>=3.7,<4.0',
     name="dLux",
     version=find_version("dLux", "__init__.py"),
     description="A fully differentiable optical simulator build in Jax",


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index aca8ca06..d35008c0 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -6,7 +6,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.7", "3.8", "3.9", "3.10"]
+        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"]
 
     steps:
       - name: checkout
diff --git a/tests/conftest.py b/tests/conftest.py
index 3a7e0db2..398590f8 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -166,8 +166,9 @@ def _create_rotate(
             angle: Array = np.array(np.pi),
             real_imaginary: bool = False,
             fourier: bool = False,
+            order: int = 1,
             padding: int = 2) -> OpticalLayer:
-        return dLux.optics.Rotate(angle, real_imaginary, fourier, padding)
+        return dLux.optics.Rotate(angle, real_imaginary, fourier, order, padding)
     return _create_rotate
 
 
diff --git a/tests/test_apertures.py b/tests/test_apertures.py
index 16061f40..f08e3bec 100644
--- a/tests/test_apertures.py
+++ b/tests/test_apertures.py
@@ -41,8 +41,7 @@ def test_all_apertures(self, create_square_aperture : callable,
             self._test_single_aperture_class(ctor)
     
     def _test_single_aperture_class(self, aperture_fixture):
-        
-        coordinates = dLux.utils.get_pixel_coordinates(512, 2./512)
+        coordinates = dLux.utils.get_pixel_positions((16, 16), (2./16, 2./16))
         
         x_offset = 1.
         y_offset = 1.
@@ -101,9 +100,11 @@ def test_range_hard(self, create_uniform_spider: callable) -> None:
         """
         Checks that the apertures fall into the correct range.
         """
-        npix = 128
+        npix = 16
         width = 2.
-        coordinates = dLux.utils.get_pixel_coordinates(npix, width / npix)
+        pixel_scale = width / npix
+        coordinates = dLux.utils.get_pixel_positions((npix,npix), 
+                                                    (pixel_scale, pixel_scale))
 
         # Case Translated 
         spider = create_uniform_spider(centre=[1., 1.], softening=0.)
@@ -130,9 +131,11 @@ def test_range_soft(self, create_uniform_spider: callable) -> None:
         """
         Checks that the aperture falls into the correct range.
         """
-        npix = 128
+        npix = 16
         width = 2.
-        coordinates = dLux.utils.get_pixel_coordinates(npix, width / npix)
+        pixel_scale = width / npix
+        coordinates = dLux.utils.get_pixel_positions((npix,npix), 
+                                                    (pixel_scale, pixel_scale))
 
         # Case Translated 
         spider = create_uniform_spider(centre=[1., 1.], softening=0.)
@@ -185,9 +188,11 @@ def test_on_aperture(self: object,
         """
         jax.config.update("jax_debug_nans", True)
         
+        npix = 16
         width = 2.
-        npix = 128
-        coordinates = dLux.utils.get_pixel_coordinates(npix, width / npix)
+        pixel_scale = width / npix
+        coordinates = dLux.utils.get_pixel_positions((npix,npix), 
+                                                    (pixel_scale, pixel_scale))
 
         ap = create_circular_aperture()
 
diff --git a/tests/test_wavefronts.py b/tests/test_wavefronts.py
index f807639a..d29ed0da 100644
--- a/tests/test_wavefronts.py
+++ b/tests/test_wavefronts.py
@@ -123,8 +123,8 @@ def test_pixel_coordinates(self, create_wavefront: callable) -> None:
         """
         wf = create_wavefront()
         assert (wf.pixel_coordinates == \
-        dLux.utils.coordinates.get_pixel_coordinates(wf.npixels,
-                                                     wf.pixel_scale)).all()
+        dLux.utils.coordinates.get_pixel_positions((wf.npixels, wf.npixels,),
+                                        (wf.pixel_scale, wf.pixel_scale))).all()
 
 
     def test_set_amplitude(self, create_wavefront: callable) -> None:
@@ -434,16 +434,18 @@ def test_rotate(self, create_wavefront: callable) -> None:
         """
         wf = create_wavefront()
         wf = dLux.CircularAperture(1.)(wf)
-        flipped_amplitude = np.flipud(wf.amplitude)
-        flipped_phase = np.flipud(wf.phase)
+        flipped_amplitude = np.flip(wf.amplitude, axis=(-1, -2))
+        flipped_phase = np.flip(wf.phase, axis=(-1, -2))
 
-        new_wf = wf.rotate(np.pi)
+        new_wf = wf.rotate(np.pi, order=1)
         assert np.allclose(new_wf.amplitude, flipped_amplitude, atol=1e-5)
         assert np.allclose(new_wf.phase, flipped_phase)
 
-        new_wf = wf.rotate(np.pi, real_imaginary=True)
+        new_wf = wf.rotate(np.pi, real_imaginary=True, order=1)
+
         assert np.allclose(new_wf.amplitude, flipped_amplitude, atol=1e-5)
-        assert np.allclose(new_wf.phase, flipped_phase)
+        # Add small remainer to fix 0-pi instability
+        assert np.allclose((new_wf.phase+1e-6)%np.pi, flipped_phase, atol=1e-5)
 
         with pytest.raises(NotImplementedError):
             wf.rotate(np.pi, fourier=True)