Merge pull request spacetelescope#937 from spacetelescope/new_WFSC_fi…

…eld_point

Add support for new sensing point and target phase map

pyproject.toml

@@ -24,7 +24,7 @@ dependencies = [
     "astropy>=5.1.0",
     "photutils>=1.10.0",
     "poppy>=1.0.0",
-    "pysiaf>=0.19.1",
+    "pysiaf>=0.23.3",
     "synphot>=1.0.0",
     "astroquery>=0.4.6",
 ]

requirements.txt

@@ -4,7 +4,7 @@ matplotlib>=3.9.1,<3.10.0
 numpy>=2.1.0,<2.2.0
 photutils>=1.13.0,<2.1.0
 poppy>=1.0.0
-pysiaf>=0.22.0,<0.24.0
+pysiaf>=0.23.3,<=0.24.1
 scipy>=1.14.0,<1.15.0
 synphot>=1.4.0,<1.5.0
 astroquery>=0.4.7,<0.5.0

webbpsf/trending.py

@@ -313,11 +313,13 @@ def wfe_histogram_plot(
             elif ote_only is True:
                 opd_data = fits.getdata(full_file_path, ext=1)
                 mask = opd_data != 0
+                sensing_apername = fits.getheader(full_file_path, ext=0)['APERNAME']
+
+                # Get WSS Target Phase Map for the sensing aperture
+                # Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024.
+                was_targ_file = webbpsf.utils.get_target_phase_map_filename(sensing_apername)
+
 
-                # Get WSS Target Phase Map
-                was_targ_file = os.path.join(
-                    webbpsf.utils.get_webbpsf_data_path(), 'NIRCam', 'OPD', 'wss_target_phase_fp1.fits'
-                )
                 target_1024 = astropy.io.fits.getdata(was_targ_file)
                 target_256 = poppy.utils.krebin(target_1024, (256, 256)) / 16
                 wf_si = target_256 * mask  # Nircam target phase map at FP1
@@ -607,6 +609,7 @@ def single_measurement_trending_plot(
     opd, opdhdu = _read_opd(filename)
     mask = opd != 0
     opd[~mask] = np.nan
+    sensing_apername = opdhdu[0].header['APERNAME']
 
     # hdr_rmswfe = opdhdu[1].header['RMS_WFE']
     visit = opdhdu[0].header['OBS_ID'][0:12]
@@ -628,18 +631,31 @@ def single_measurement_trending_plot(
         print('     Previous OPD is:', prev_filename)
 
     prev_opd, prev_opd_hdu = _read_opd(prev_filename)
+    prev_sensing_apername = prev_opd_hdu[0].header['APERNAME']
 
     if subtract_target:
         if verbose:
             print('     Subtracting NIRCam SI WFE target phase map')
 
-        # Get WSS Target Phase Map
-        was_targ_file = os.path.join(webbpsf.utils.get_webbpsf_data_path(), 'NIRCam', 'OPD', 'wss_target_phase_fp1.fits')
+        # Get WSS Target Phase Map for the sensing aperture
+        # Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024.
+        was_targ_file = webbpsf.utils.get_target_phase_map_filename(sensing_apername)
+        prev_was_targ_file = webbpsf.utils.get_target_phase_map_filename(prev_sensing_apername)
+
+
         target_1024 = astropy.io.fits.getdata(was_targ_file)
         target_256 = poppy.utils.krebin(target_1024, (256, 256)) / 16  # scale factor for rebinning w/out increasing values
 
+        if prev_was_targ_file != was_targ_file:
+            prev_target_1024 = astropy.io.fits.getdata(prev_was_targ_file)
+            prev_target_256 = poppy.utils.krebin(prev_target_1024,
+                                                 (256, 256)) / 16  # scale factor for rebinning w/out increasing values
+        else:
+            prev_target_256 = target_256
+
+
         opd -= target_256
-        prev_opd -= target_256
+        prev_opd -= prev_target_256
 
     # Compute deltas and decompose
     deltatime = get_datetime_utc(opdhdu, return_as='astropy') - get_datetime_utc(prev_opd_hdu, return_as='astropy')
@@ -961,9 +977,13 @@ def vprint(*text):
 
     opd, opdhdu = _read_opd(opdtable[which_opds_mask][0]['fileName'])
     mask = opd != 0
+    sensing_apername = opdhdu[0].header['APERNAME']
+
+    # Get WSS Target Phase Map for the sensing aperture
+    # Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024.
+
+    was_targ_file = webbpsf.utils.get_target_phase_map_filename(sensing_apername)
 
-    # Get WSS Target Phase Map
-    was_targ_file = os.path.join(webbpsf.utils.get_webbpsf_data_path(), 'NIRCam', 'OPD', 'wss_target_phase_fp1.fits')
     target_1024 = astropy.io.fits.getdata(was_targ_file)
     target_256 = poppy.utils.krebin(target_1024, (256, 256))
 

webbpsf/utils.py

@@ -1048,3 +1048,25 @@ def label_wavelength(nwavelengths, wavelength_slices):
     else:
         raise ValueError('Maximum number of wavelengths exceeded. ' 'Cannot be more than 10,000.')
     return label
+
+
+def get_target_phase_map_filename(apername):
+    """Get WSS Target Phase Map for the specified aperture
+    Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024.
+    """
+    path = os.getenv('WEBBPSF_PATH')
+    if apername == 'NRCA3_FP1':
+        fn = 'wss_target_phase_fp1.fits'
+    elif apername == 'NRCA1_FP6':
+        fn = 'wss_target_phase_fp6.fits'
+    else:
+        raise ValueError(f"Target phase map not available for aperture = {apername}")
+
+    was_targ_file = os.path.join(
+        get_webbpsf_data_path(), 'NIRCam', 'OPD', fn)
+
+    if not os.path.exists(was_targ_file):
+        raise ValueError("File wss_target_phase_{}.fits, \
+        not found under {}.".format(apername.split('_')[1].lower(),path))
+
+    return was_targ_file

webbpsf/webbpsf_core.py

@@ -1720,10 +1720,13 @@ def load_wss_opd(self, filename, output_path=None, backout_si_wfe=True, verbose=
                 # note that there is a slight focus offset between the two wavelengths, due to NIRCam's refractive design
             # Set to the sensing aperture, and retrieve the OPD there
             sensing_inst.set_position_from_aperture_name(sensing_apername)
-            # special case: for the main sensing point FP1, we use the official WAS target phase map, rather than the
-            # WebbPSF-internal SI WFE model.
-            was_targ_file = os.path.join(utils.get_webbpsf_data_path(), 'NIRCam', 'OPD', 'wss_target_phase_fp1.fits')
-            if sensing_apername == 'NRCA3_FP1' and os.path.exists(was_targ_file):
+            # special case: for the main sensing points FP1 or FP6, we use the official WAS target phase map,
+            # rather than the WebbPSF-internal SI WFE model.
+
+            # Select correct target phase map based on sensing field point.
+            # Note that the sensing maintenance program changed field point from NRC A3 to A1 around Dec 2024.
+            if sensing_apername in ['NRCA3_FP1', 'NRCA1_FP6']:
+                was_targ_file = utils.get_target_phase_map_filename(sensing_apername)
                 sensing_fp_si_wfe = poppy.FITSOpticalElement(opd=was_targ_file).opd
             else:
                 sensing_fp_si_wfe = sensing_inst.get_wfe('si')