diff --git a/docs/_img/python/jitter.py b/docs/_img/python/jitter.py
index 3e48386..bf34c95 100644
--- a/docs/_img/python/jitter.py
+++ b/docs/_img/python/jitter.py
@@ -5,8 +5,8 @@
 mpl.rcParams['figure.figsize'] = (4.5, 4.5)
 
 mask = lentil.circlemask((256, 256), 120)
-phase = lentil.zernike_compose(mask, [0, 0, 0, -300e-9, 50e-9, -100e-9, 50e-9])
-pupil = lentil.Pupil(amplitude=mask, phase=phase, focal_length=10, pixelscale=1 / 240)
+opd = lentil.zernike_compose(mask, [0, 0, 0, -300e-9, 50e-9, -100e-9, 50e-9])
+pupil = lentil.Pupil(amplitude=mask, opd=opd, focal_length=10, pixelscale=1 / 240)
 w = lentil.Wavefront(650e-9)
 w *= pupil
 w = lentil.propagate_dft(w, pixelscale=5e-6, shape=64, oversample=5)
diff --git a/docs/_img/python/pixelate.py b/docs/_img/python/pixelate.py
index b30ad69..b8f83f5 100644
--- a/docs/_img/python/pixelate.py
+++ b/docs/_img/python/pixelate.py
@@ -5,8 +5,8 @@
 mpl.rcParams['figure.figsize'] = (4.5, 4.5)
 
 mask = lentil.circlemask((256, 256), 120)
-phase = lentil.zernike_compose(mask, [0, 0, 0, -300e-9, 50e-9, -100e-9, 50e-9])
-pupil = lentil.Pupil(amplitude=mask, phase=phase, focal_length=10, pixelscale=1 / 240)
+opd = lentil.zernike_compose(mask, [0, 0, 0, -300e-9, 50e-9, -100e-9, 50e-9])
+pupil = lentil.Pupil(amplitude=mask, opd=opd, focal_length=10, pixelscale=1 / 240)
 w = lentil.Wavefront(650e-9)
 w *= pupil
 w = lentil.propagate_dft(w, pixelscale=5e-6, shape=64, oversample=5)
diff --git a/docs/_img/python/plane_transformations.py b/docs/_img/python/plane_transformations.py
index 0b11d49..df3804a 100644
--- a/docs/_img/python/plane_transformations.py
+++ b/docs/_img/python/plane_transformations.py
@@ -3,7 +3,7 @@
 
 mask = lentil.util.circle((256, 256), 128)
 opd = lentil.zernike.zernike(mask, 8) * 500e-9  # 500 nm of coma
-pupil = lentil.Pupil(amplitude=mask, phase=opd,
+pupil = lentil.Pupil(amplitude=mask, opd=opd,
                      focal_length=10, pixelscale=1/256)
 rotation = lentil.Rotate(angle=30, unit='degrees')
 flip = lentil.Flip(1)
diff --git a/docs/_img/python/simple_example.py b/docs/_img/python/simple_example.py
index 10658f1..deb78f8 100644
--- a/docs/_img/python/simple_example.py
+++ b/docs/_img/python/simple_example.py
@@ -4,7 +4,7 @@
 mask = lentil.util.circle((256, 256), 128) - lentil.util.circle((256, 256), 128/3)
 opd = lentil.zernike.zernike_compose(mask, coeffs=[0, 0, 0, 300e-9, 50e-9, -100e-9, 50e-9])
 
-pupil = lentil.Pupil(amplitude=mask, phase=opd, focal_length=10,
+pupil = lentil.Pupil(amplitude=mask, opd=opd, focal_length=10,
                      pixelscale=1/256)
 detector = lentil.Image(pixelscale=5e-6)
 
diff --git a/docs/_img/python/smear.py b/docs/_img/python/smear.py
index 31d5f6f..d7b510e 100644
--- a/docs/_img/python/smear.py
+++ b/docs/_img/python/smear.py
@@ -5,8 +5,8 @@
 mpl.rcParams['figure.figsize'] = (4.5, 4.5)
 
 mask = lentil.circlemask((256, 256), 120)
-phase = lentil.zernike_compose(mask, [0, 0, 0, -300e-9, 50e-9, -100e-9, 50e-9])
-pupil = lentil.Pupil(amplitude=mask, phase=phase, focal_length=10, pixelscale=1 / 240)
+opd = lentil.zernike_compose(mask, [0, 0, 0, -300e-9, 50e-9, -100e-9, 50e-9])
+pupil = lentil.Pupil(amplitude=mask, opd=opd, focal_length=10, pixelscale=1 / 240)
 w = lentil.Wavefront(650e-9)
 w *= pupil
 w = lentil.propagate_dft(w, pixelscale=5e-6, shape=64, oversample=5)
diff --git a/docs/_img/python/smear_directional.py b/docs/_img/python/smear_directional.py
index fbe5281..c1b7f31 100644
--- a/docs/_img/python/smear_directional.py
+++ b/docs/_img/python/smear_directional.py
@@ -5,8 +5,8 @@
 mpl.rcParams['figure.figsize'] = (4.5, 4.5)
 
 mask = lentil.circlemask((256, 256), 120)
-phase = lentil.zernike_compose(mask, [0, 0, 0, -300e-9, 50e-9, -100e-9, 50e-9])
-pupil = lentil.Pupil(amplitude=mask, phase=phase, focal_length=10, pixelscale=1 / 240)
+opd = lentil.zernike_compose(mask, [0, 0, 0, -300e-9, 50e-9, -100e-9, 50e-9])
+pupil = lentil.Pupil(amplitude=mask, opd=opd, focal_length=10, pixelscale=1 / 240)
 w = lentil.Wavefront(650e-9)
 w *= pupil
 w = lentil.propagate_dft(w, pixelscale=5e-6, shape=64, oversample=5)
diff --git a/docs/dev/tech_notes/prop_algorithm.rst b/docs/dev/tech_notes/prop_algorithm.rst
index 709028a..88e9e5c 100644
--- a/docs/dev/tech_notes/prop_algorithm.rst
+++ b/docs/dev/tech_notes/prop_algorithm.rst
@@ -9,9 +9,9 @@ The general propagation algorithm implemented by :class:`prop.Propagate` is as f
 .. code::
 
     for each plane in planes
-        cache plane amplitude, phase, ptt_vector
+        cache plane amplitude, opd, ptt_vector
 
 
 
     for each plane in planes
-        clear plane amplitude, phase, ptt_vector cache
\ No newline at end of file
+        clear plane amplitude, opd, ptt_vector cache
\ No newline at end of file
diff --git a/docs/examples/planes/rb_element.rst b/docs/examples/planes/rb_element.rst
index 9cf76e3..1ea23ea 100644
--- a/docs/examples/planes/rb_element.rst
+++ b/docs/examples/planes/rb_element.rst
@@ -16,7 +16,7 @@ The example below assumes a ``dwdx`` influence function matrix is available.
             self._x = np.zeros(6)
 
         @property
-        def phase(self):
+        def opd(self):
             return np.einsum('ijk,i->jk',DWDX, self.x)
 
         @property
@@ -62,7 +62,7 @@ influence function matrix.
             self._uerror = 1e-6
 
         @property
-        def phase(self):
+        def opd(self):
             return np.einsum('ijk,i->jk',DWDX, self.x)
 
         @property
diff --git a/docs/user/basics.optical_systems.rst b/docs/user/basics.optical_systems.rst
index aa2cd7f..2bb3524 100644
--- a/docs/user/basics.optical_systems.rst
+++ b/docs/user/basics.optical_systems.rst
@@ -50,7 +50,7 @@ multiplied in two ways:
         >>> w1 = plane * w0
 
 The :func:`~Plane.multiply` method constructs a complex phasor from the plane's
-:attr:`~lentil.Plane.amplitude` and :attr:`~lentil.Plane.phase` attributes and the
+:attr:`~lentil.Plane.amplitude` and :attr:`~lentil.Plane.opd` attributes and the
 |Wavefront| wavelength. The plane complex phasor is then multiplied element-wise with
 the wavefront's complex data array:
 
@@ -77,15 +77,15 @@ represented by a single |Pupil| plane:
     >>> amplitude = lentil.circle(shape=(256, 256), radius=120)
     >>> opd = lentil.zernike_compose(mask=amplitude,
     ...                              coeffs=[0, 0, 0, 100e-9, 300e-9, 0, -100e-9])
-    >>> pupil = lentil.Pupil(amplitude=amplitude, phase=opd, focal_length=10,
+    >>> pupil = lentil.Pupil(amplitude=amplitude, opd=opd, focal_length=10,
     ...                      pixelscale=1/240)
-    >>> plt.imshow(pupil.phase, origin='lower')
+    >>> plt.imshow(pupil.opd, origin='lower')
 
 Segmented optical systems
 =========================
 Creating a model of a segmented aperture optical system in Lentil doesn't require any
 special treatment. The |Plane| and |Pupil| objects work the same with sparse or
-segmented amplitude, phase, and mask attributes as with monolithic ones.
+segmented amplitude, opd, and mask attributes as with monolithic ones.
 
 That being said, it is advantageous from a performance point of view to supply a
 3-dimensional `segment mask` when specifying a Plane's :attr:`~lentil.Plane.mask`
diff --git a/docs/user/basics.planes.rst b/docs/user/basics.planes.rst
index 8251202..3087046 100644
--- a/docs/user/basics.planes.rst
+++ b/docs/user/basics.planes.rst
@@ -52,7 +52,7 @@ plane. A plane is defined by the following parameters:
 
 * :attr:`~lentil.Plane.amplitude` - Defines the relative electric field amplitude
   transmission through the plane
-* :attr:`~lentil.Plane.phase` - Defines the electric field phase shift that a wavefront
+* :attr:`~lentil.Plane.opd` - Defines the optical path difference that a wavefront
   experiences when propagating through the plane
 * :attr:`~lentil.Plane.mask` - Defines the binary mask over which the plane data is
   valid. If `mask` is 2-dimensional, the plane is assumed to be monolithic. If `mask`
@@ -93,15 +93,15 @@ Once a Plane is defined, its attributes can be modified at any time:
     :scale: 50
 
     >>> p = lentil.Plane(amplitude=lentil.util.circle((256,256), 120))
-    >>> p.phase = 2e-6 * lentil.zernike(p.mask, index=4)
-    >>> plt.imshow(p.phase, origin='lower')
+    >>> p.opd = 2e-6 * lentil.zernike(p.mask, index=4)
+    >>> plt.imshow(p.opd, origin='lower')
 
 
 Resampling or rescaling a Plane
 -------------------------------
 It is possible to resample a plane using either the :func:`~lentil.Plane.resample`
 or :func:`~lentil.Plane.rescale` methods. Both methods use intrepolation to
-resample the amplitude, phase, and mask attributes and readjust the pixelscale
+resample the amplitude, opd, and mask attributes and readjust the pixelscale
 attribute as necessary. The default behavior is to perform this interpolation
 on a copy of the plane, but it is possible to operate in-place by setting
 ``inplace=True``.
@@ -111,7 +111,7 @@ on a copy of the plane, but it is possible to operate in-place by setting
 Fitting and removing Plane tilt
 -------------------------------
 The plane's :func:`~lentil.Plane.fit_tilt` method performs a least squares fit to
-estimate and remove tilt from the phase attribute. The tilt removed from the phase
+estimate and remove tilt from the opd attribute. The tilt removed from the opd
 attribute is accounted for by appending an equivalent :class:`~lentil.Tilt` object
 to the plane's :attr:`~lentil.Plane.tilt` attribute. The default behavior is to
 perform this operation on a copy of the plane, but it is possible to operate
@@ -147,9 +147,8 @@ Discretely sampled pupil attributes can also be specified:
 
 * :attr:`~lentil.Pupil.amplitude` - Defines the relative electric field amplitude
   transmission through the pupil
-* :attr:`~lentil.Pupil.phase` - Defines the electric field phase shift that a wavefront
-  experiences when propagating through the pupil. This term is commonly known as the
-  optical path difference (OPD).
+* :attr:`~lentil.Pupil.opd` - Defines the optical path difference that a wavefront
+  experiences when propagating through the pupil.
 * :attr:`~lentil.Pupil.mask` - Defines the binary mask over which the pupil data is
   valid. If `mask` is 2-dimensional, the pupil is assumed to be monolithic. If `mask`
   is 3-dimensional, the pupil is assumed to be segmented with the segment masks
@@ -165,7 +164,7 @@ Create a pupil with:
 
 .. code-block:: pycon
 
-    >>> p = lentil.Pupil(focal_length=10, pixelscale=1/100, amplitude=1, phase=0)
+    >>> p = lentil.Pupil(focal_length=10, pixelscale=1/100, amplitude=1, opd=0)
 
 Image
 =====
@@ -180,7 +179,7 @@ of the following can be specified:
   the image plane will grow as necessary to capture all data.
 * :attr:`~lentil.Image.amplitude` - Definers the relative electric field amplitude
   transmission through the image plane.
-* :attr:`~lentil.Image.phase` - Defines the electric field phase shift that a wavefront
+* :attr:`~lentil.Image.opd` - Defines the optical path difference that a wavefront
   experiences when propagating through the image plane.
 
 Detector
@@ -318,12 +317,12 @@ efficiently modeling a grism.
 
 Active optics and deformable mirrors
 ====================================
-Active optics and deformable mirrors are easily represented by defining a phase that
+Active optics and deformable mirrors are easily represented by defining an OPD that
 depends on some parameterized state. Because there is no standard architecture for these
 types of optical elements, Lentil does not provide a concrete implementation. Instead,
 a custom subclass of either |Plane| or |Pupil| should be defined. The exact
 implementation details will vary by application, but a simple example of a tip-tilt
-mirror where the plane's phase is computed dynamically based on the state `x` is
+mirror where the plane's OPD is computed dynamically based on the state `x` is
 provided below. Additional examples can be found in Model Patterns under
 :ref:`patterns.planes`.
 
@@ -346,14 +345,14 @@ provided below. Additional examples can be found in Model Patterns under
                                                        normalize=False)
 
         @property
-        def phase(self):
+        def opd(self):
             return np.einsum('ijk,i->jk', self._infl_fn, self.x)
 
 .. code-block:: pycon
 
     >>> tt = TipTiltMirror()
     >>> tt.x = [1e-6, 3e-6]
-    >>> plt.imshow(tt.phase)
+    >>> plt.imshow(tt.opd)
     >>> plt.colorbar()
 
 .. plot::
@@ -363,9 +362,9 @@ provided below. Additional examples can be found in Model Patterns under
     import lentil
 
     mask = lentil.circlemask((256,256), 120)
-    phase = lentil.zernike_compose(mask, [0, 1e-6, 3e-6], normalize=False)
+    opd = lentil.zernike_compose(mask, [0, 1e-6, 3e-6], normalize=False)
 
-    im = plt.imshow(phase, origin='lower')
+    im = plt.imshow(opd, origin='lower')
     plt.colorbar(im, fraction=0.046, pad=0.04)
 
 Customizing Plane
@@ -373,7 +372,7 @@ Customizing Plane
 The Plane class or any of the classes derived from Plane can be subclassed to modify
 any of the default behavior. Reasons to do this may include but are not limited to:
 
-* Dynamically computing the :attr:`~lentil.Plane.phase` attribute
+* Dynamically computing the :attr:`~lentil.Plane.opd` attribute
 * Changing the Plane-Wavefront interaction by redefining the `Plane.multiply()` method
 * Modifying the way a Plane is resampled or rescaled
 
@@ -390,9 +389,9 @@ Some general guidance for how to safely subclass Plane is provided below.
     passing these attributes along to the ``super().__init__()`` call to ensure they are
     properly set.
 
-Redefining the amplitude, phase, or mask attributes
+Redefining the amplitude, OPD, or mask attributes
 ---------------------------------------------------
-Plane :attr:`~lentil.Plane.amplitude`, :attr:`~lentil.Plane.phase`, and
+Plane :attr:`~lentil.Plane.amplitude`, :attr:`~lentil.Plane.opd`, and
 :attr:`~lentil.Plane.mask` are all defined as properties, but Python allows you to
 redefine them as class attributes without issue:
 
@@ -403,10 +402,10 @@ redefine them as class attributes without issue:
     class CustomPlane(le.Plane):
         def __init__(self):
             self.amplitude = lentil.circle((256,256), 128)
-            self.phase = lentil.zernike(lentil.circlemask((256,256),128), 4)
+            self.opd = lentil.zernike(lentil.circlemask((256,256),128), 4)
 
 If more dynamic behavior is required, the property can be redefined. For example, to
-return a new random phase each time the :attr:`~lentil.Plane.phase` attribute is
+return a new random OPD each time the :attr:`~lentil.Plane.opd` attribute is
 accessed:
 
 .. code-block:: python3
@@ -420,11 +419,11 @@ accessed:
             self.amplitude = lentil.circle((256,256), 128)
 
         @property
-        def phase(self):
+        def phaopdse(self):
             return lentil.zernike_compose(self.mask, np.random.random(10))
 
-It is also straightforward to implement a custom :attr:`~lentil.Plane.phase` property to
-provide a stateful phase attribute:
+It is also straightforward to implement a custom :attr:`~lentil.Plane.opd` property to
+provide a stateful OPD attribute:
 
 .. code-block:: python3
 
@@ -438,12 +437,12 @@ provide a stateful phase attribute:
             self.x = x
 
         @property
-        def phase(self):
+        def opd(self):
             return lentil.zernike_compose(self.mask, self.x)
 
 .. note::
 
-    Polychromatic or broadband diffraction propagations access the phase, amplitude,
+    Polychromatic or broadband diffraction propagations access the OPD, amplitude,
     and mask attributes for each propagatioon wavelength. Because these attributes
     remain fixed during a propagation, it is inefficient to repeatedly recompute
     them. To mitigate this, it can be very useful to provide a mechanism for freezing
@@ -462,14 +461,14 @@ provide a stateful phase attribute:
                 self.amplitude = lentil.circle((256,256), 128)
 
             @property
-            def phase(self):
+            def opd(self):
                 return lentil.zernike_compose(self.mask, np.random.random(10))
 
             def freeze(self):
-                # Return a copy of CustomPlane with the phase attribute redefined
-                # to be a static copy of the phase when freeze() is called
+                # Return a copy of CustomPlane with the OPD attribute redefined
+                # to be a static copy of the OPD when freeze() is called
                 out = copy.deepcopy(self)
-                out.phase = self.phase.copy()
+                out.opd = self.opd.copy()
                 return out
 
 
diff --git a/docs/user/basics.wavefront_error.rst b/docs/user/basics.wavefront_error.rst
index 6d85c37..9a44f36 100644
--- a/docs/user/basics.wavefront_error.rst
+++ b/docs/user/basics.wavefront_error.rst
@@ -5,7 +5,7 @@ Representing wavefront error
 ****************************
 
 Wavefront error is represented in a :class:`~lentil.Plane` by specifying its
-:attr:`~lentil.Plane.phase` attribute. For static errors, a simple wavefront error map is
+:attr:`~lentil.Plane.opd` attribute. For static errors, a simple wavefront error map is
 sufficient. For more complicated errors that are random or time-varying in nature, a
 more dynamic and/or state-based approach is required.
 
@@ -14,7 +14,7 @@ Lentil to a Numpy array.
 
 .. note::
 
-    For :class:`~lentil.Pupil` planes, the :attr:`~lentil.Pupil.phase` attribute represents the optical
+    For :class:`~lentil.Pupil` planes, the :attr:`~lentil.Pupil.opd` attribute represents the optical
     path difference (OPD) relative to the pupil's reference sphere.
 
 .. _user.wavefront_error.sign:
@@ -71,8 +71,8 @@ containing the JWST NITCam static wavefront error:
 
     >>> import numpy as np
     >>> import lentil
-    >>> phase = np.load('path/to/nircam_wfe.npy')
-    >>> pupil = lentil.Pupil(focal_length=119.77, pixelscale=6.6035/1024, phase=phase)
+    >>> opd = np.load('path/to/nircam_wfe.npy')
+    >>> pupil = lentil.Pupil(focal_length=119.77, pixelscale=6.6035/1024, opd=opd)
 
 .. image:: /_static/img/nircam.png
     :scale: 50
diff --git a/docs/user/performance.rst b/docs/user/performance.rst
index ae67e02..2e6ce74 100644
--- a/docs/user/performance.rst
+++ b/docs/user/performance.rst
@@ -109,8 +109,8 @@ to clear any cached values.
 
 The cached attributes are defined in a list in each Plane's
 :attr:`~lentil.Plane.cache_attrs` attribute. This list is user-settable but the only
-valid values are `'amplitude'` and `'phase'`. The default behavior is to cache both
-amplitude and phase attributes.
+valid values are `'amplitude'` and `'opd'`. The default behavior is to cache both
+amplitude and OPD attributes.
 
 DFT matrices
 ------------
diff --git a/docs/user/plots/focus_images.py b/docs/user/plots/focus_images.py
index e8492d7..b6982fb 100644
--- a/docs/user/plots/focus_images.py
+++ b/docs/user/plots/focus_images.py
@@ -13,13 +13,13 @@
 focus = lentil.zernike(mask=np.ones((256, 256)), index=4)
 
 pupil_neg = lentil.Pupil(amplitude=amp, pixelscale=1/240, focal_length=10)
-pupil_neg.phase = -6e-6 * focus
+pupil_neg.opd = -6e-6 * focus
 w_neg = lentil.Wavefront(650e-9)
 w_neg *= pupil_neg
 w_neg = lentil.propagate_dft(w_neg, pixelscale=5e-6, shape=200, oversample=2)
 
 pupil_pos = lentil.Pupil(amplitude=amp, pixelscale=1/240, focal_length=10)
-pupil_pos.phase = 6e-6 * focus
+pupil_pos.opd = 6e-6 * focus
 w_pos = lentil.Wavefront(650e-9)
 w_pos *= pupil_pos
 w_pos = lentil.propagate_dft(w_pos, pixelscale=5e-6, shape=200, oversample=2)
diff --git a/docs/user/plots/npix_prop.py b/docs/user/plots/npix_prop.py
index 69a062a..5666f41 100644
--- a/docs/user/plots/npix_prop.py
+++ b/docs/user/plots/npix_prop.py
@@ -6,7 +6,7 @@
 
 amp = lentil.circle((256, 256), 120)
 opd = lentil.zernike(amp, 4) * 1e-6
-pupil = lentil.Pupil(amplitude=amp, phase=opd, pixelscale=1 / 240, focal_length=10)
+pupil = lentil.Pupil(amplitude=amp, opd=opd, pixelscale=1 / 240, focal_length=10)
 
 w1 = lentil.Wavefront(wavelength=500e-9)
 w1 *= pupil
diff --git a/docs/user/plots/tilt_images.py b/docs/user/plots/tilt_images.py
index a116a31..ae7494d 100644
--- a/docs/user/plots/tilt_images.py
+++ b/docs/user/plots/tilt_images.py
@@ -6,12 +6,12 @@
 x_tilt = 8e-6 * lentil.zernike(amp, 3)  # +x tilt
 y_tilt = 8e-6 * lentil.zernike(amp, 2)  # +y tilt
 
-py = lentil.Pupil(focal_length=10, pixelscale=1 / 240, amplitude=amp, phase=y_tilt)
+py = lentil.Pupil(focal_length=10, pixelscale=1 / 240, amplitude=amp, opd=y_tilt)
 wy = lentil.Wavefront(650e-9)
 wy *= py
 wy = lentil.propagate_dft(wy, pixelscale=5e-6, shape=200, oversample=5)
 
-px = lentil.Pupil(focal_length=10, pixelscale=1 /240, amplitude=amp, phase=x_tilt)
+px = lentil.Pupil(focal_length=10, pixelscale=1 /240, amplitude=amp, opd=x_tilt)
 wx = lentil.Wavefront(650e-9)
 wx *= px
 wx = lentil.propagate_dft(wx, pixelscale=5e-6, shape=200, oversample=5)
diff --git a/docs/user/quickstart.rst b/docs/user/quickstart.rst
index 32a477d..d84d210 100644
--- a/docs/user/quickstart.rst
+++ b/docs/user/quickstart.rst
@@ -52,7 +52,7 @@ a focal length of 20 meters and a diameter of 1 meter:
     :context: close-figs
     :include-source:
 
-    >>> pupil = lentil.Pupil(amplitude=amp, phase=opd, pixelscale=1/240, 
+    >>> pupil = lentil.Pupil(amplitude=amp, opd=opd, pixelscale=1/240, 
     ...                      focal_length=20)
 
 Note the diameter is implicitly defined via the 
diff --git a/lentil/plane.py b/lentil/plane.py
index b16efa3..ca739d4 100644
--- a/lentil/plane.py
+++ b/lentil/plane.py
@@ -23,9 +23,9 @@ class Plane:
         through diffraction propagation is required. If not specified (default),
         amplitude is created which has no effect on wavefront propagation. Can
         also be specified using the ``amp`` keyword.
-    phase : array_like, optional
-        Phase change caused by plane. If not specified (default), phase is
-        created which has no effect on wavefront propagation.
+    opd : array_like, optional
+        Optical path difference (OPD) induced by plane. If not specified (default), 
+        zero OPD is created which has no effect on wavefront propagation.
     mask : array_like, optional
         Binary mask. If not specified, a mask is created from the amplitude.
         If ``mask`` has 2 dimensions, the plane is assumed to be monolithic. If
@@ -46,7 +46,7 @@ class Plane:
         Plane type
     """
 
-    def __init__(self, amplitude=1, phase=0, mask=None, pixelscale=None, diameter=None,
+    def __init__(self, amplitude=1, opd=0, mask=None, pixelscale=None, diameter=None,
                  ptype=None, **kwargs):
         
         if 'amp' in kwargs.keys():
@@ -56,7 +56,7 @@ def __init__(self, amplitude=1, phase=0, mask=None, pixelscale=None, diameter=No
             amplitude = kwargs['amp']
 
         self.amplitude = np.asarray(amplitude)
-        self.phase = np.asarray(phase)
+        self.opd = np.asarray(opd)
         self.mask = mask
         self.pixelscale = pixelscale
         self.diameter = diameter
@@ -84,19 +84,19 @@ def amplitude(self, value):
         self._amplitude = np.asarray(value)
     
     @property
-    def phase(self):
-        """Electric field phase
+    def opd(self):
+        """Optical path difference
 
         Returns
         -------
         ndarray
 
         """
-        return self._phase
+        return self._opd
     
-    @phase.setter
-    def phase(self, value):
-        self._phase = np.asarray(value)
+    @opd.setter
+    def opd(self, value):
+        self._opd = np.asarray(value)
 
     @property
     def mask(self):
@@ -262,7 +262,7 @@ def copy(self):
 
     def fit_tilt(self, inplace=False):
         """
-        Fit and remove tilt from Plane :attr:`phase` via least squares. The
+        Fit and remove tilt from Plane :attr:`opd` via least squares. The
         equivalent angular tilt is bookkept in Plane :attr:`tilt`.
 
         Parameters
@@ -284,27 +284,27 @@ def fit_tilt(self, inplace=False):
 
         # There are a couple of cases where we don't have enough information to remove the
         # tilt, so we just return the Plane as-is
-        if ptt_vector is None or plane.phase.size == 1:
+        if ptt_vector is None or plane.opd.size == 1:
             return plane
 
         if self.size == 1:
-            t = np.linalg.lstsq(ptt_vector.T, plane.phase.ravel(), rcond=None)[0]
-            phase_tilt = np.einsum('ij,i->j', ptt_vector[1:3], t[1:3])
-            plane.phase -= phase_tilt.reshape(plane.phase.shape)
+            t = np.linalg.lstsq(ptt_vector.T, plane.opd.ravel(), rcond=None)[0]
+            opd_tilt = np.einsum('ij,i->j', ptt_vector[1:3], t[1:3])
+            plane.opd -= opd_tilt.reshape(plane.opd.shape)
             plane.tilt.append(Tilt(x=t[1], y=t[2]))
 
         else:
             t = np.empty((self.size, 3))
-            phase_no_tilt = np.empty((self.size, plane.phase.shape[0], plane.phase.shape[1]))
+            opd_no_tilt = np.empty((self.size, plane.opd.shape[0], plane.opd.shape[1]))
 
             # iterate over the segments and compute the tilt term
             for seg in np.arange(self.size):
-                t[seg] = np.linalg.lstsq(ptt_vector[3 * seg:3 * seg + 3].T, plane.phase.ravel(),
+                t[seg] = np.linalg.lstsq(ptt_vector[3 * seg:3 * seg + 3].T, plane.opd.ravel(),
                                          rcond=None)[0]
                 seg_tilt = np.einsum('ij,i->j', ptt_vector[3 * seg + 1:3 * seg + 3], t[seg, 1:3])
-                phase_no_tilt[seg] = (plane.phase - seg_tilt.reshape(plane.phase.shape)) * self.mask[seg]
+                opd_no_tilt[seg] = (plane.opd - seg_tilt.reshape(plane.opd.shape)) * self.mask[seg]
 
-            plane.phase = np.sum(phase_no_tilt, axis=0)
+            plane.opd = np.sum(opd_no_tilt, axis=0)
             plane.tilt.extend([Tilt(x=t[seg, 1], y=t[seg, 2]) for seg in range(self.size)])
 
         return plane
@@ -316,7 +316,7 @@ def rescale(self, scale):
 
         * `Plane.amplitude` is rescaled via 3rd order spline interpolation. The
             result is scaled to preserve total power.
-        * `Plane.phase` is rescaled via 3rd order spline interpolation.
+        * `Plane.opd` is rescaled via 3rd order spline interpolation.
         * `Plane.mask` is rescaled via 0-order nearest neighbor interpolation.
         * `Plane.pixelscale` is adjusted appropriately.
 
@@ -346,9 +346,9 @@ def rescale(self, scale):
                                                 mask=None, order=3, mode='nearest',
                                                 unitary=False)/scale
 
-        if plane.phase.ndim > 1:
-            plane.phase = lentil.rescale(plane.phase, scale=scale, shape=None, mask=None,
-                                         order=3, mode='nearest', unitary=False)
+        if plane.opd.ndim > 1:
+            plane.opd = lentil.rescale(plane.opd, scale=scale, shape=None, mask=None,
+                                       order=3, mode='nearest', unitary=False)
 
         # because plane.mask is automatically computed from amplitude if it is not
         # provided, we really only want to reinterpolate if a mask was provided (stored
@@ -377,7 +377,7 @@ def resample(self, pixelscale):
 
         * `Plane.amplitude` is resampled via 3rd order spline interpolation. The
           result is scaled to preserve total power.
-        * `Plane.phase` is resampled via 3rd order spline interpolation.
+        * `Plane.opd` is resampled via 3rd order spline interpolation.
         * `Plane.mask` is resampled via 0-order nearest neighbor interpolation.
         * `Plane.pixelscale` is adjusted appropriately.
 
@@ -449,10 +449,10 @@ def multiply(self, wavefront):
                 # slice of the amplitude array may contain parts of adjacent segments
                 mask = self.mask if self.size == 1 else self.mask[n]
                 amp = self.amplitude if self.amplitude.size == 1 else self.amplitude[s] * mask[s]
-                phase = self.phase if self.phase.size == 1 else self.phase[s]
+                opd = self.opd if self.opd.size == 1 else self.opd[s]
 
                 # construct complex phasor
-                phasor = Field(data=amp*np.exp(2*np.pi*1j*phase/wavefront.wavelength),
+                phasor = Field(data=amp*np.exp(2*np.pi*1j*opd/wavefront.wavelength),
                                pixelscale=self.pixelscale,
                                offset=lentil.helper.slice_offset(s, self.shape),
                                tilt=[self.tilt[n]] if self.tilt else [])
@@ -569,9 +569,9 @@ class Pupil(Plane):
         through a diffraction propagation is required. If not specified, a
         default amplitude is created which has no effect on wavefront
         propagation. Can also be specified using the ``amp`` keyword.
-    phase : array_like, optional
-        Phase change caused by plane. If not specified, a default phase is
-        created which has no effect on wavefront propagation.
+    opd : array_like, optional
+        Optical path difference (OPD) induced by plane. If not specified (default), 
+        zero OPD is created which has no effect on wavefront propagation.
     mask : array_like, optional
         Binary mask. If not specified, a mask is created from the amplitude.
         If ``mask`` has 2 dimensions, the plane is assumed to be monolithic. If
@@ -590,9 +590,9 @@ class Pupil(Plane):
     """
 
     def __init__(self, focal_length=None, pixelscale=None, amplitude=1,
-                 phase=0, mask=None, **kwargs):
+                 opd=0, mask=None, **kwargs):
 
-        super().__init__(pixelscale=pixelscale, amplitude=amplitude, phase=phase,
+        super().__init__(pixelscale=pixelscale, amplitude=amplitude, opd=opd,
                          mask=mask, ptype=lentil.pupil, **kwargs)
 
         self.focal_length = focal_length
@@ -645,12 +645,13 @@ class Image(Plane):
     Detector
 
     """
-
-    def __init__(self, shape=None, pixelscale=None, amplitude=1, mask=None,
-                 **kwargs):
-        super().__init__(amplitude=amplitude, mask=mask,
+    def __init__(self, shape=None, pixelscale=None, amplitude=1, opd=0, 
+                 mask=None, **kwargs):
+        
+        super().__init__(amplitude=amplitude, opd=opd, mask=mask, 
                          pixelscale=pixelscale, ptype=lentil.image,
                          **kwargs)
+
         self.shape = shape
 
     @property
@@ -962,13 +963,13 @@ def _arc_len(dist_func, a, b):
         return scipy.integrate.quad(dist_func, a, b)[0]
         
 
-class DispersivePhase(Plane):
+class DispersiveAberration(Plane):
 
     def multiply(self, wavefront):
-        # NOTE: we can handle wavelength-dependent phase terms here (e.g. chromatic
-        # aberrations). Since the phase will vary by wavelength, we can't fit out the
+        # NOTE: we can handle wavelength-dependent OPD terms here (e.g. chromatic
+        # aberrations). Since the OPD will vary by wavelength, we can't fit out the
         # tilt pre-propagation and apply the same tilt for each wavelength like we can
-        # with run of the mill tilt
+        # with run-of-the-mill tilt
         raise NotImplementedError
 
 
diff --git a/lentil/wfe.py b/lentil/wfe.py
index d8af20f..49c285a 100644
--- a/lentil/wfe.py
+++ b/lentil/wfe.py
@@ -63,12 +63,12 @@ def power_spectrum(mask, pixelscale, rms, half_power_freq, exp, seed=None):
     # Generate noise, filter it to realize the requested PSD, and enforce
     # the pupil mask
     noise = rng.normal(size=[n, m])
-    phase = np.real(np.fft.ifft2(np.fft.fft2(noise) * H)) * np.sqrt(m * n)
+    opd = np.real(np.fft.ifft2(np.fft.fft2(noise) * H)) * np.sqrt(m * n)
 
-    phase *= mask
-    phase = phase * np.sqrt(np.count_nonzero(phase)/np.sum(np.abs(phase)**2)) * rms
+    opd *= mask
+    opd = opd * np.sqrt(np.count_nonzero(opd)/np.sum(np.abs(opd)**2)) * rms
 
-    return phase
+    return opd
 
 
 def translation_defocus(mask, f_number, translation):
diff --git a/tests/fixtures/pupil.py b/tests/fixtures/pupil.py
index e81e290..5fa2381 100644
--- a/tests/fixtures/pupil.py
+++ b/tests/fixtures/pupil.py
@@ -16,7 +16,7 @@ def _pupil(focal_length, diameter, shape, radius, coeffs=None):
             opd = 0
 
         p = lentil.Pupil(amplitude=amplitude,
-                         phase=opd,
+                         opd=opd,
                          pixelscale=pixelscale,
                          focal_length=focal_length)
         
diff --git a/tests/test_plane.py b/tests/test_plane.py
index 3d80de5..d607f06 100644
--- a/tests/test_plane.py
+++ b/tests/test_plane.py
@@ -12,18 +12,18 @@ class RandomPlane(lentil.Plane):
     def __init__(self):
         super().__init__(pixelscale=1,
                          amplitude=np.random.uniform(size=size),
-                         phase=np.random.uniform(size=size))
+                         opd=np.random.uniform(size=size))
 
 
 def test_default_plane():
     # Ensure that a default Plane creates an object that won't have any
-    # impact on an optical system (a perfect optic with no phase change and
-    # perfect optical and spectral transmission).
+    # impact on an optical system (a perfect optic with no wavefront error 
+    # and perfect optical and spectral transmission).
 
     p = lentil.Plane()
     assert p.pixelscale == None
     assert np.all(p.amplitude == 1)
-    assert np.all(p.phase == 0)
+    assert np.all(p.opd == 0)
     assert p.mask == p.amplitude
 
 
@@ -52,7 +52,7 @@ def test_wavefront_plane_multiply():
     w1 = p.multiply(w)
 
     slc = lentil.helper.boundary_slice(p.mask)
-    phasor = p.amplitude[slc] * np.exp(2*np.pi*1j*p.phase[slc]/w.wavelength)
+    phasor = p.amplitude[slc] * np.exp(2*np.pi*1j*p.opd[slc]/w.wavelength)
 
     assert np.array_equal(w1.data[0].data, phasor)
 
@@ -79,14 +79,14 @@ def __init__(self):
         super().__init__(focal_length=10,
                          pixelscale=2/256,
                          amplitude=lentil.util.circle((256, 256), 128),
-                         phase=np.zeros((256, 256)))
+                         opd=np.zeros((256, 256)))
 
 
 def test_wavefront_pupil_multiply():
     p = CircularPupil()
     w = lentil.wavefront.Wavefront(650e-9)
     w = p.multiply(w)
-    phasor = p.amplitude * np.exp(1j*p.phase * 2 * np.pi / w.wavelength)
+    phasor = p.amplitude * np.exp(1j*p.opd * 2 * np.pi / w.wavelength)
 
     assert np.array_equal(w.data[0].data, phasor)
     assert w.focal_length == p.focal_length
diff --git a/tests/test_propagate.py b/tests/test_propagate.py
index 3f7beb0..ec37d15 100644
--- a/tests/test_propagate.py
+++ b/tests/test_propagate.py
@@ -135,7 +135,7 @@ def __init__(self, npix, diameter=1, coeffs=None):
         super().__init__(focal_length=10,
                          pixelscale=1/npix,
                          amplitude=amplitude,
-                         phase=opd,
+                         opd=opd,
                          mask=mask)
 
         self.diameter = diameter
@@ -179,30 +179,30 @@ def test_propagate_airy():
     assert np.all(np.isclose(psf, psf_airy, atol=1e-3))
 
 
-def test_propagate_tilt_angle():
+def test_propagate_fit_tilt():
 
     p = TiltPupil(npix=256)
 
-    w_phase = lentil.Wavefront(650e-9)
-    w_phase = p.multiply(w_phase)
-    w_phase = lentil.propagate_dft(w_phase, shape=128, pixelscale=5e-6, oversample=2)
-    psf_phase = w_phase.intensity
+    w = lentil.Wavefront(650e-9)
+    w = p.multiply(w)
+    w = lentil.propagate_dft(w, shape=128, pixelscale=5e-6, oversample=2)
+    psf = w.intensity
 
     p.fit_tilt(inplace=True)
-    w_angle = lentil.Wavefront(650e-9)
-    w_angle = w_angle * p
-    w_angle = lentil.propagate_dft(w_angle, shape=128, pixelscale=5e-6, oversample=2)
-    psf_angle = w_angle.intensity
+    w_fit_tilt = lentil.Wavefront(650e-9)
+    w_fit_tilt = w_fit_tilt * p
+    w_fit_tilt = lentil.propagate_dft(w_fit_tilt, shape=128, pixelscale=5e-6, oversample=2)
+    psf_fit_tilt = w_fit_tilt.intensity
 
     # threshold the PSFs so that the centroiding is consistent
-    psf_phase[psf_phase < 1e-5] = 0
-    psf_angle[psf_angle < 1e-5] = 0
+    psf[psf < 1e-5] = 0
+    psf_fit_tilt[psf_fit_tilt < 1e-5] = 0
 
-    delta = np.abs(np.asarray(lentil.util.centroid(psf_phase)) - np.asarray(lentil.util.centroid(psf_angle)))
+    delta = np.abs(np.asarray(lentil.util.centroid(psf)) - np.asarray(lentil.util.centroid(psf_fit_tilt)))
     assert np.all(delta <= 1e-10)
 
 
-def test_propagate_tilt_phase_analytic():
+def test_propagate_tilt_analytic():
     oversample = 10
     pixelscale = 5e-6
     npix = np.array([64, 64])
@@ -231,7 +231,7 @@ def test_propagate_tilt_phase_analytic():
     assert np.all((np.abs(shift - analytic_shift)/oversample) < 0.2)
 
 
-def test_propagate_tilt_angle_analytic():
+def test_propagate_tilt_fit_tilt_analytic():
     oversample = 10
     pixelscale = 5e-6
     npix = np.array([64, 64])
@@ -267,7 +267,7 @@ def test_propagate_resample():
     coeffs = np.random.uniform(low=-200e-9, high=200e-9, size=6)
     opd = lentil.zernike_compose(amp, coeffs)
 
-    p = lentil.Pupil(focal_length=10, pixelscale=1 / 240, amplitude=amp, phase=opd)
+    p = lentil.Pupil(focal_length=10, pixelscale=1 / 240, amplitude=amp, opd=opd)
     w = lentil.Wavefront(650e-9)
     w *= p
     wi = lentil.propagate_dft(w, shape=(64,64), pixelscale=5e-6, oversample=10)
diff --git a/tests/test_propagate_segmented.py b/tests/test_propagate_segmented.py
index f2d1b71..ecb1ad6 100644
--- a/tests/test_propagate_segmented.py
+++ b/tests/test_propagate_segmented.py
@@ -20,7 +20,7 @@ def __init__(self, npix=256):
         opd = opdp + opd1 + opd2 + opd3
 
         super().__init__(focal_length=10, pixelscale=1/npix,
-                         amplitude=lentil.normalize_power(global_mask), phase=opd,
+                         amplitude=lentil.normalize_power(global_mask), opd=opd,
                          mask=mask)
 
 
@@ -31,20 +31,20 @@ def test_propagate_tilt_angle_mono():
     w1 = lentil.Wavefront(wavelength=650e-9)
     w1 *= p
     w1 = lentil.propagate_dft(w1, shape=(128,128), pixelscale=5e-6)
-    psf_phase = w1.intensity
+    psf = w1.intensity
 
     p.fit_tilt(inplace=True)
     w2 = lentil.Wavefront(wavelength=650e-9)
     w2 *= p
     w2 = lentil.propagate_dft(w2, shape=(128,128), pixelscale=5e-6)
-    psf_angle = w2.intensity
+    psf_fit_tilt = w2.intensity
 
     # Normalize and threshold the PSFs so that the centroiding is consistent
-    psf_phase /= np.max(psf_phase)
-    psf_phase[psf_phase < 0.2] = 0
+    psf /= np.max(psf)
+    psf[psf < 0.2] = 0
 
-    psf_angle /= np.max(psf_angle)
-    psf_angle[psf_angle < 0.2] = 0
+    psf_fit_tilt /= np.max(psf_fit_tilt)
+    psf_fit_tilt[psf_fit_tilt < 0.2] = 0
 
-    delta = np.abs(np.asarray(lentil.util.centroid(psf_phase)) - np.asarray(lentil.util.centroid(psf_angle)))
+    delta = np.abs(np.asarray(lentil.util.centroid(psf)) - np.asarray(lentil.util.centroid(psf_fit_tilt)))
     assert np.all(delta <= 0.2)
diff --git a/tests/test_propagate_slice.py b/tests/test_propagate_slice.py
index 630f804..fd317a0 100644
--- a/tests/test_propagate_slice.py
+++ b/tests/test_propagate_slice.py
@@ -21,20 +21,20 @@ def test_propagate_slice_one():
     rho, theta = lentil.zernike_coordinates(amp, shift=shift)
     coeffs = np.random.uniform(low=-1, high=1, size=11)*100e-9
     coeffs[0:3] = 0 # no ptt
-    phase = lentil.zernike_compose(amp, coeffs, rho=rho, theta=theta)
+    opd = lentil.zernike_compose(amp, coeffs, rho=rho, theta=theta)
 
     alpha = (dx*du)/(wavelength*focal_length)
     oversample = 5
     shape = 64
 
-    phasor = amp*np.exp(-2*1j*np.pi*phase/wavelength)
+    phasor = amp*np.exp(-2*1j*np.pi*opd/wavelength)
     F = lentil.fourier.dft2(phasor, alpha=alpha/oversample, shape=shape*oversample)
 
 
     slc = lentil.helper.boundary_slice(amp)
     ofst = lentil.helper.slice_offset(slc, shape=amp.shape)
 
-    phasor = amp[slc]*np.exp(-2*1j*np.pi*phase[slc]/wavelength)
+    phasor = amp[slc]*np.exp(-2*1j*np.pi*opd[slc]/wavelength)
     F_slc = lentil.fourier.dft2(phasor, alpha=alpha/oversample, shape=shape*oversample, offset=ofst)
 
     assert np.allclose(F, F_slc)
@@ -57,7 +57,7 @@ def test_propagate_slice_multi():
     rho, theta = lentil.zernike_coordinates(amp1, shift=(0, -0.3*n))
     coeffs = np.random.uniform(low=-1, high=1, size=11)*100e-9
     coeffs[0:3] = 0 # no ptt
-    phase1 = lentil.zernike_compose(amp1, coeffs, rho=rho, theta=theta)
+    opd1 = lentil.zernike_compose(amp1, coeffs, rho=rho, theta=theta)
 
 
     amp2 = lentil.util.circle((n,n), n//5, shift=(0, .3*n))
@@ -67,22 +67,22 @@ def test_propagate_slice_multi():
     rho, theta = lentil.zernike_coordinates(amp2, shift=(0, 0.3*n))
     coeffs = np.random.uniform(low=-1, high=1, size=11)*100e-9
     coeffs[0:3] = 0 # no ptt
-    phase2 = lentil.zernike_compose(amp2, coeffs, rho=rho, theta=theta)
+    opd2 = lentil.zernike_compose(amp2, coeffs, rho=rho, theta=theta)
 
     alpha = (dx*du)/(wavelength*focal_length)
     oversample = 5
     shape = 64
 
     amp = amp1 + amp2
-    phase = phase1 + phase2
+    opd = opd1 + opd2
 
-    phasor = amp*np.exp(-2*1j*np.pi*phase/wavelength)
+    phasor = amp*np.exp(-2*1j*np.pi*opd/wavelength)
     F = lentil.fourier.dft2(phasor, alpha=alpha/oversample, shape=shape*oversample)
 
-    phasor1 = amp1[slc1]*np.exp(-2*1j*np.pi*phase1[slc1]/wavelength)
+    phasor1 = amp1[slc1]*np.exp(-2*1j*np.pi*opd1[slc1]/wavelength)
     F1 = lentil.fourier.dft2(phasor1, alpha=alpha/oversample, shape=shape*oversample, offset=ofst1)
 
-    phasor2 = amp2[slc2]*np.exp(-2*1j*np.pi*phase2[slc2]/wavelength)
+    phasor2 = amp2[slc2]*np.exp(-2*1j*np.pi*opd2[slc2]/wavelength)
     F2 = lentil.fourier.dft2(phasor2, alpha=alpha/oversample, shape=shape*oversample, offset=ofst2)
 
     F_slc = F1 + F2
diff --git a/tests/test_wfe.py b/tests/test_wfe.py
index 04727f8..53fda78 100644
--- a/tests/test_wfe.py
+++ b/tests/test_wfe.py
@@ -5,10 +5,10 @@
 def test_power_spectrum_wfe():
     mask = lentil.util.circlemask((256, 256), 128)
     rms = 50e-9
-    phase = lentil.wfe.power_spectrum(mask, pixelscale=1/256, rms=rms,
+    opd = lentil.wfe.power_spectrum(mask, pixelscale=1/256, rms=rms,
                                       half_power_freq=5, exp=3)
 
-    assert np.std(phase[np.nonzero(phase)])/rms >= 0.8
+    assert np.std(opd[np.nonzero(opd)])/rms >= 0.8
 
 
 def test_translation_defocus():
@@ -16,7 +16,7 @@ def test_translation_defocus():
     f_number = np.random.normal(10)
     translation = np.random.uniform(low=-0.5e-3, high=0.5e-3)
 
-    phase = lentil.wfe.translation_defocus(mask, f_number, translation)
+    opd = lentil.wfe.translation_defocus(mask, f_number, translation)
     pv_defocus = translation/(8*f_number**2)
 
-    assert np.isclose(np.abs(pv_defocus), np.abs(np.max(phase) - np.min(phase)))
+    assert np.isclose(np.abs(pv_defocus), np.abs(np.max(opd) - np.min(opd)))
diff --git a/tests/test_zernike.py b/tests/test_zernike.py
index 828c497..4dfd88f 100644
--- a/tests/test_zernike.py
+++ b/tests/test_zernike.py
@@ -21,16 +21,16 @@ def test_zernike_basis_vectorize():
 def test_zernike_fit():
     mask = lentil.util.circlemask((256, 256), 128)
     coeffs = np.random.rand(4)*100e-9
-    phase = lentil.zernike_compose(mask, coeffs)
-    fit_coeffs = lentil.zernike_fit(phase, mask, np.arange(1, 5))
+    opd = lentil.zernike_compose(mask, coeffs)
+    fit_coeffs = lentil.zernike_fit(opd, mask, np.arange(1, 5))
     assert np.all(np.isclose(coeffs, fit_coeffs))
 
 
 def test_zernike_remove():
     mask = lentil.util.circlemask((256, 256), 128)
     coeffs = np.random.rand(4)*100e-9
-    phase = lentil.zernike_compose(mask, coeffs)
-    residual = lentil.zernike_remove(phase, mask, np.arange(1, 5))
+    opd = lentil.zernike_compose(mask, coeffs)
+    residual = lentil.zernike_remove(opd, mask, np.arange(1, 5))
     assert np.all(np.isclose(residual, np.zeros_like(residual)))
 
 def test_zernike_random_mask():