Merge pull request #908 from spacetelescope/webbpsf_load_wss_pupil_size

Enhance capabilities to simulate PSFs with larger JWST pupil

docs/index.rst

@@ -60,6 +60,7 @@ Contents
    jwst_measured_opds.ipynb
    jwst_detector_effects.ipynb
    jwst_matching_psfs_to_data.ipynb
+   jwst_large_psf.ipynb
    jwst_optical_budgets.ipynb
    jwst_psf_subtraction.ipynb
    jwst_wavefront_deltas.ipynb

webbpsf/optics.py

@@ -1223,7 +1223,7 @@ def _get_initial_pupil_sampling(instrument):
         if isinstance(instrument.pupil, fits.HDUList):
             pupilheader = instrument.pupil[0].header
         else:
-            pupilfile = os.path.join(instrument._datapath, 'OPD', instrument.pupil)
+            pupilfile = os.path.join(utils.get_webbpsf_data_path(), instrument.pupil)
             pupilheader = fits.getheader(pupilfile)
 
         npix = pupilheader['NAXIS1']

webbpsf/webbpsf_core.py

@@ -1627,6 +1627,7 @@ def load_wss_opd(self, filename, output_path=None, backout_si_wfe=True, verbose=
                 an estimate of the OTE wavefront nominally at the master chief ray location (between the NIRCams).
                 WebbPSF will automatically add back on top of this the OTE field dependent WFE for the appropriate
                 field point. as usual.
+            - Scale the OPD to match the same size of the user provide pupil file
 
         Parameters
         ----------
@@ -1654,6 +1655,15 @@ def load_wss_opd(self, filename, output_path=None, backout_si_wfe=True, verbose=
             which is most of use for some other tool.
 
         """
+        # We use the size of the user supplied name of the JWST pupil in order to create the matching size OPD
+        # The code assume the naming convention for the JWST pupil file: jwst_pupil_RevW_npix<size in pixels>.fits.gz
+        npix_out = int(self.pupil[self.pupil.find('npix') + len('npix'):self.pupil.find('.fits')])
+
+        if verbose and npix_out != 1024:
+            print(
+                  f'The size of the JWST pupil is different than nominal (1024px), {self.pupil}. '
+                  f'The OPD will be scaled accordingly'
+            )
 
         # If the provided filename doesn't exist on the local disk, try retrieving it from MAST
         # Note, this will automatically use cached versions downloaded previously, if present
@@ -1662,12 +1672,12 @@ def load_wss_opd(self, filename, output_path=None, backout_si_wfe=True, verbose=
 
         if verbose:
             print(f'Importing and format-converting OPD from {filename}')
-        opdhdu = webbpsf.mast_wss.import_wss_opd(filename)
+        opdhdu = webbpsf.mast_wss.import_wss_opd(filename, npix_out=npix_out)
 
         # Mask out any pixels in the OPD array which are outside the OTE pupil.
         # This is mostly cosmetic, and helps mask out some edge effects from the extrapolation + interpolation in
         # resizing the OPDs
-        ote_pupil_mask = utils.get_pupil_mask() != 0
+        ote_pupil_mask = utils.get_pupil_mask(npix=npix_out) != 0
         opdhdu[0].data *= ote_pupil_mask
 
         # opdhdu[0].header['RMS_OBS'] = (webbpsf.utils.rms(opdhdu[0].data, mask=ote_pupil_mask)*1e9,
@@ -1698,6 +1708,10 @@ def load_wss_opd(self, filename, output_path=None, backout_si_wfe=True, verbose=
             sensing_inst.pupil = (
                 self.pupil
             )  # handle the case if the user has selected a different NPIX other than the default 1024
+            sensing_inst.pupilopd = (
+                opdhdu
+            )  # handle the case if the user has selected a different NPIX other than the default 1024
+
             if sensing_inst.name == 'NRC':
                 sensing_inst.filter = 'F212N'
                 # TODO: optionally check for the edge case in which the sensing was done in F187N
@@ -1712,6 +1726,11 @@ def load_wss_opd(self, filename, output_path=None, backout_si_wfe=True, verbose=
             else:
                 sensing_fp_si_wfe = sensing_inst.get_wfe('si')
 
+            if npix_out != 1024:   # handle the case if the user has selected a different NPIX other than the default
+                # the results from the zoom function preserve the STD between both phase maps and
+                # the total sum between the phase maps is proportional to the zoom value
+                sensing_fp_si_wfe = scipy.ndimage.zoom(sensing_fp_si_wfe, npix_out / 1024)
+
             sensing_fp_ote_wfe = sensing_inst.get_wfe('ote_field_dep')
 
             sihdu = fits.ImageHDU(sensing_fp_si_wfe)
@@ -1746,7 +1765,7 @@ def load_wss_opd(self, filename, output_path=None, backout_si_wfe=True, verbose=
             if plot or save_ote_wfe:
                 # Either of these options will need the total OTE WFE.
                 # Under normal circumstances webbpsf will compute this later automatically, but if needed we do it here too
-                selected_fp_ote_wfe = self.get_wfe('ote_field_dep')
+                selected_fp_ote_wfe = sensing_inst.get_wfe('ote_field_dep')
                 total_ote_wfe_at_fp = opdhdu[0].data + (selected_fp_ote_wfe * ote_pupil_mask)
 
             if plot: