diff --git a/webbpsf/distortion.py b/webbpsf/distortion.py
index 5b380c18..0b03b031 100644
--- a/webbpsf/distortion.py
+++ b/webbpsf/distortion.py
@@ -77,7 +77,8 @@ def apply_distortion(HDUlist_or_filename=None, fill_value=0):
     oversamp = psf[ext].header["OVERSAMP"]  # will be 1 for ext=1
     xpix_center = psf[ext].header["DET_X"]  # center x location in pixels
     ypix_center = psf[ext].header["DET_Y"]  # center y location in pixels
-    length = psf[ext].shape[0]  # can assume square
+    len_y = psf[ext].shape[0]
+    len_x = psf[ext].shape[1]
 
     # Convert the PSF center point from pixels to arcseconds using pysiaf
     xarc_center, yarc_center = aper.sci_to_idl(xpix_center, ypix_center)
@@ -85,11 +86,11 @@ def apply_distortion(HDUlist_or_filename=None, fill_value=0):
     # ###############################################
     # Create an array of indices (in pixels) for where the PSF is located on the detector
     # 1) Set up blank indices (in pixels)
-    ypix, xpix = np.indices((length, length), dtype=float)
+    ypix, xpix = np.indices((len_y, len_x), dtype=float)
 
     # 2) Shift indices to be centered on (0,0) (starting to transform into the Ideal frame)
-    ypix -= (length - 1.) / 2.
-    xpix -= (length - 1.) / 2.
+    ypix -= (len_y - 1.) / 2.
+    xpix -= (len_x - 1.) / 2.
 
     # 3) Convert these indices from pixels to arcseconds
     # Note: This also shifts the oversampled indices so they span the same region as the detector-sampled indices
@@ -108,11 +109,11 @@ def apply_distortion(HDUlist_or_filename=None, fill_value=0):
     # ###############################################
     # Create an array of indices (in pixels) that the final data will be interpolated on to
     # 1) Set up blank indices (in pixels)
-    ynew, xnew = np.indices([length, length], dtype=float)
+    ynew, xnew = np.indices([len_y, len_x], dtype=float)
 
     # 2) Shift indices to be in the Ideal frame (centered on 0)
-    xnew -= (length - 1.) / 2.
-    ynew -= (length - 1.) / 2.
+    xnew -= (len_x - 1.) / 2.
+    ynew -= (len_y - 1.) / 2.
 
     # 3) Shift the oversampled indices so they span the same region as the detector-sampled indices
     # Note: the oversampled array is still longer by a factor of the oversample
diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py
index a521ed54..37401256 100644
--- a/webbpsf/webbpsf_core.py
+++ b/webbpsf/webbpsf_core.py
@@ -750,6 +750,9 @@ def calc_psf(self, outfile=None, source=None, nlambda=None, monochromatic=None,
                                                 save_intermediates=save_intermediates,
                                                 return_intermediates=return_intermediates, normalize=normalize)
 
+        if return_intermediates:
+            psf, intermediates = psf
+
         # If chosen to add distortion
         if add_distortion:
 
@@ -779,17 +782,18 @@ def calc_psf(self, outfile=None, source=None, nlambda=None, monochromatic=None,
                 psf_distorted.writeto(outfile, overwrite=True, output_verify='ignore')  # already created in the 1st calc_psf
                 _log.info("Re-saved result with distortion to " + outfile)
 
-            return psf_distorted
+            psf = psf_distorted
+
+        if return_intermediates:
+            return psf, intermediates
 
         else:
             return psf
 
-
     calcPSF = calc_psf
     calc_psf.__doc__ += SpaceTelescopeInstrument.calc_psf.__doc__  # allow users to see poppy calc_psf docstring too
 
 
-
 class MIRI(JWInstrument):
     """ A class modeling the optics of MIRI, the Mid-InfraRed Instrument.