combining issues 726 and 736

webbpsf/__init__.py

@@ -128,4 +128,4 @@ class Conf(_config.ConfigNamespace):
 
 from .jupyter_gui import show_notebook_interface
 
-from .match_data import setup_sim_to_match_file
+from .match_data import setup_sim_to_match_file

webbpsf/detectors.py

@@ -170,11 +170,15 @@ def apply_detector_ipc(psf_hdulist, extname = 'DET_DIST'):
         webbpsf.webbpsf_core._log.debug(f"Skipping IPC simulation since ext {extname} is not found")
         return
 
+    # This avoid applying IPC effect simulations twice
+    keyword = 'IPCINST'
+    if keyword in psf_hdulist[extname].header._keyword_indices:
+        return
+
     inst = psf_hdulist[extname].header['INSTRUME'].upper()
     oversample = psf_hdulist[extname].header['OVERSAMP']
 
     kernel, meta = get_detector_ipc_model(inst, psf_hdulist[extname].header)
-
     if kernel is not None:
 
         if inst.upper()=='NIRCAM':

webbpsf/gridded_library.py

@@ -291,6 +291,12 @@ def create_grid(self):
                 if self.verbose is True:
                     print("    Position {}/{}: {} pixels".format(i+1, len(self.location_list), loc))
 
+                add_ipc_gridded = False   # add variable to keep track of if this function locally changed the user's IPC input
+                # Deactivate IPC corrections, if any, before calc_psf as we are applying them later
+                if self.webb.options.get('add_ipc', True):
+                    self.webb.options['add_ipc'] = False
+                    add_ipc_gridded = True   # we should reset the IPC flag below, after computing the PSF
+
                 # Create PSF
                 psf = self.webb.calc_psf(**self._kwargs)
                 if self.verbose is True:
@@ -302,10 +308,9 @@ def create_grid(self):
                 psf[ext].data = astropy.convolution.convolve(psf[ext].data, kernel)
 
                 # Convolve PSF with a model for interpixel capacitance
-                # note, normally this is applied in calc_psf to the detector-sampled data;
-                # here we specially apply this to the oversampled data
-                if self.add_distortion and self.webb.options.get('add_ipc', True):
+                if self.add_distortion and add_ipc_gridded:
                     webbpsf.detectors.apply_detector_ipc(psf, extname=ext)
+                    self.webb.options['add_ipc'] = True  # restore the user's value for the IPC option
 
                 # Add PSF to 5D array
                 psf_arr[i, :, :] = psf[ext].data
@@ -390,7 +395,7 @@ def create_grid(self):
             # copy all the jitter-related keys (the exact set of keywords varies based on jitter type)
             # Also copy charge diffusion and IPC
             for k in psf[ext].header.keys():   # do the rest
-                if k.startswith('JITR') or k.startswith('IPC') or k.startswith('CHDF'):
+                if k.startswith('JITR') or k.startswith('IPC') or k.startswith('PPC') or k.startswith('CHDF'):
                     meta[k] = (psf[ext].header[k], psf[ext].header.comments[k])
 
             meta["DATE"] = (psf[ext].header["DATE"], "Date of calculation")

webbpsf/match_data.py

@@ -6,6 +6,7 @@
 import scipy.optimize
 
 import webbpsf
+import pysiaf
 
 
 def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choice='closest'):
@@ -37,6 +38,10 @@
         # webbpsf doesn't model the MIRI LRS prism spectral response
         print("Please note, webbpsf does not currently model the LRS spectral response. Setting filter to F770W instead.")
         inst.filter='F770W'
+    elif (inst.name == 'NIRCam') and (header['PUPIL'][0] == 'F') and (header['PUPIL'][-1] in ['N', 'M']):
+        # These NIRCam filters are physically in the pupil wheel, but still act as filters.
+        # Grab the filter name from the PUPIL keyword in this case.
+        inst.filter = header['PUPIL']
     else:
         inst.filter=header['filter']
     inst.set_position_from_aperture_name(header['APERNAME'])
@@ -48,10 +53,22 @@
     # per-instrument specializations
     if inst.name == 'NIRCam':
         if header['PUPIL'].startswith('MASK'):
-            inst.pupil_mask = header['PUPIL']
+            if header['PUPIL'] == 'MASKBAR':
+                inst.pupil_mask = header['CORONMSK'].replace('MASKA', 'MASK')
+            else:
+                inst.pupil_mask = header['PUPIL']
             inst.image_mask = header['CORONMSK'].replace('MASKA', 'MASK')  # note, have to modify the value slightly for
                                                                            # consistency with the labels used in webbpsf
-            inst.set_position_from_aperture_name(header['APERNAME'])   # Redo this, in case the image_mask setting auto switched it to something else
+            # The apername keyword is not always correct for cases with dual-channel coronagraphy
+            # in some such cases, APERNAME != PPS_APER. Let's ensure we have the proper apername for this channel:
+            apername = get_nrc_coron_apname(header)
+            inst.set_position_from_aperture_name(apername)
+
+        elif header['PUPIL'].startswith('F'):
+            inst.filter = header['PUPIL']
+        else:
+            inst.pupil_mask = header['PUPIL']
+
     elif inst.name == 'MIRI':
         if inst.filter in ['F1065C', 'F1140C', 'F1550C']:
             inst.image_mask = 'FQPM'+inst.filter[1:5]
@@ -78,3 +95,119 @@
 
     return inst
 
+
+
+def get_nrc_coron_apname(input):
+    """Get NIRCam coronagraph aperture name from header or data model
+
+    Handles edge cases for dual-channel coronagraphy.
+
+    By Jarron Leisenring originally in webbpsf_ext, copied here by permission
+
+    Parameters
+    ==========
+    input : fits.header.Header or datamodels.DataModel
+        Input header or data model
+    """
+
+    if isinstance(input, (fits.header.Header)):
+        # Aperture names
+        apname = input['APERNAME']
+        apname_pps = input['PPS_APER']
+        subarray = input['SUBARRAY']
+    else:
+        # Data model meta info
+        meta = input.meta
+
+        # Aperture names
+        apname = meta.aperture.name
+        apname_pps = meta.aperture.pps_name
+        subarray = meta.subarray.name
+
+    # print(apname, apname_pps, subarray)
+
+    # No need to do anything if the aperture names are the same
+    # Also skip if MASK not in apname_pps
+    if ((apname==apname_pps) or ('MASK' not in apname_pps)) and ('400X256' not in subarray):
+        apname_new = apname
+    else:
+        # Should only get here if coron mask and apname doesn't match PPS
+        apname_str_split = apname.split('_')
+        sca = apname_str_split[0]
+        image_mask = get_nrc_coron_mask_from_pps_apername(apname_pps)
+
+        # Get subarray info
+        # Sometimes apname erroneously has 'FULL' in it
+        # So, first for subarray info in apname_pps 
+        if ('400X256' in apname_pps) or ('400X256' in subarray):
+            apn0 = f'{sca}_400X256'
+        elif ('FULL' in apname_pps):
+            apn0 = f'{sca}_FULL'
+        else:
+            apn0 = sca
+
+        apname_new = f'{apn0}_{image_mask}'
+
+        # Append filter or NARROW if needed
+        pps_str_arr = apname_pps.split('_')
+        last_str = pps_str_arr[-1]
+        # Look for filter specified in PPS aperture name
+        if ('_F1' in apname_pps) or ('_F2' in apname_pps) or ('_F3' in apname_pps) or ('_F4' in apname_pps):
+            # Find all instances of "_"
+            inds = [pos for pos, char in enumerate(apname_pps) if char == '_']
+            # Filter is always appended to end, but can have different string sizes (F322W2)
+            filter = apname_pps[inds[-1]+1:]
+            apname_new += f'_{filter}'
+        elif last_str=='NARROW':
+            apname_new += '_NARROW'
+        elif ('TAMASK' in apname_pps) and ('WB' in apname_pps[-1]):
+            apname_new += '_WEDGE_BAR'
+        elif ('TAMASK' in apname_pps) and (apname_pps[-1]=='R'):
+            apname_new += '_WEDGE_RND'
+
+    # print(apname_new)
+
+    # If apname_new doesn't exit, we need to fall back to apname
+    # even if it may not fully make sense.
+    if apname_new in pysiaf.Siaf('NIRCam').apernames:
+        return apname_new
+    else:
+        return apname
+
+
+def get_nrc_coron_mask_from_pps_apername(apname_pps):
+    """Get NIRCam coronagraph mask name from PPS aperture name
+
+    The PPS aperture name is of the form:
+        NRC[A/B][1-5]_[FULL]_[TA][MASK]
+    where MASK is the name of the coronagraphic mask used.
+
+    For target acquisition apertures the mask name can be
+    prependend with "TA" (eg., TAMASK335R).
+
+    Return '' if MASK not in input aperture name.
+    """
+
+    if 'MASK' not in apname_pps:
+        return ''
+
+    pps_str_arr = apname_pps.split('_')
+    for s in pps_str_arr:
+        if 'MASK' in s:
+            image_mask = s
+            break
+
+    # Special case for TA apertures
+    if 'TA' in image_mask:
+        # Remove TA from mask name
+        image_mask = image_mask.replace('TA', '')
+
+        # Remove FS from mask name
+        if 'FS' in image_mask:
+            image_mask = image_mask.replace('FS', '')
+
+        # Remove trailing S or L from LWB and SWB TA apertures
+        if ('WB' in image_mask) and (image_mask[-1]=='S' or image_mask[-1]=='L'):
+            image_mask = image_mask[:-1]
+
+    return image_mask

webbpsf/tests/test_detectors.py

@@ -169,3 +169,27 @@ def test_ipc_oversampling_equivalence(oversamp = 2):
     psf_detdist_v2 = poppy.utils.rebin_array(testpsf['OVERDIST'].data, (oversamp,oversamp))
 
     assert np.allclose(psf_detdist, psf_detdist_v2), "PSFs calculated should be equivalent for IPC convolution and binning in either order"
+
+
+def test_ipc_basic_effect_on_psf_fwhm():
+    """A basic test that the IPC model has the expected effect: making PSFs slightly broader.
+
+    Tests that (a) there is no change to the first two extensions, which are 'pure optical PSF'
+               (b) that the FWHM increases for the other two extensions, which are distortion+detector effects
+    """
+    nrc = webbpsf_core.NIRCam()
+    psf_withipc = nrc.calc_psf(nlambda=1, fov_pixels=101)
+    nrc.options['add_ipc'] = False
+    psf_noipc = nrc.calc_psf(nlambda=1, fov_pixels=101)
+
+    for extname in ['OVERSAMP', 'DET_SAMP']:
+        fwhm_ipc = poppy.measure_fwhm(psf_withipc, ext=extname)
+        fwhm_noipc = poppy.measure_fwhm(psf_noipc, ext=extname)
+        assert fwhm_ipc==fwhm_noipc, f'Adding IPC should not have any effect on the {extname} data.'
+        print(f'test ok for {extname}')
+
+    for extname in ['OVERDIST', 'DET_DIST']:
+        fwhm_ipc = poppy.measure_fwhm(psf_withipc, ext=extname)
+        fwhm_noipc = poppy.measure_fwhm(psf_noipc, ext=extname)
+        assert fwhm_ipc>fwhm_noipc, f'Adding IPC should not blur and increase the FWHM in the {extname} data.'
+        print(f'test ok for {extname}')

webbpsf/webbpsf_core.py

@@ -1118,6 +1118,8 @@
         # Pull values from options dictionary
         add_distortion = options.get('add_distortion', True)
         crop_psf = options.get('crop_psf', True)
+        # you can turn on/off IPC corrections via the add_ipc option, default True.
+        add_ipc = options.get('add_ipc', True)
 
         # Add distortion if set in calc_psf
         if add_distortion:
@@ -1166,11 +1168,11 @@
 
         # Rewrite result variable based on output_mode; this includes binning down to detector sampling.
         SpaceTelescopeInstrument._calc_psf_format_output(self, result, options)
-        # you can turn on/off IPC corrections via the add_ipc option, default True.
-        add_ipc = options.get('add_ipc', True)
-        if add_ipc and add_distortion:
-            result = detectors.apply_detector_ipc(result)  # apply detector IPC model (after binning to detector sampling)
 
+        if add_ipc and add_distortion and ('DET_DIST' in result):
+            result = detectors.apply_detector_ipc(result)  # apply detector IPC model (after binning to detector sampling)
+        if add_ipc and add_distortion and ('OVERDIST' in result):
+            result = detectors.apply_detector_ipc(result, extname = 'OVERDIST')  # apply detector IPC model to oversampled PSF
 
     def interpolate_was_opd(self, array, newdim):
         """ Interpolates an input 2D  array to any given size.
@@ -2434,13 +2436,31 @@
             SAM_box_size = 5.0
         elif ((self.image_mask == 'MASKSWB') or (self.image_mask == 'MASKLWB')):
             bar_offset = self.options.get('bar_offset', None)
-            # If the bar offset is not provided, use the filter name to lookup the default
+            # If the bar offset is not provided, use the SIAF aperture name, or else the filter name to lookup the default
             # position. If an offset is provided and is a floating point value, use that
             # directly as the offset. Otherwise assume it's a filter name and try passing
             # that in to the auto offset. (that allows for selecting the narrow position, or
             # for simulating using a given filter at some other filter's position.)
+            # This code is somewhat convoluted, for historical reasons and back-compatibility
             if bar_offset is None:
-                auto_offset = self.filter
+                # Try to use the SIAF aperture name to determine the offset
+                # This can help better automate simulations matching data, since match_data.py will
+                # have copied the aperturename from the header, like NRCA5_MASKLWB_NARROW or similar
+                if 'MASK' in self.aperturename:
+                    apname_last_part = self.aperturename.split('_')[-1]
+                    if apname_last_part=='NARROW':
+                        auto_offset = 'narrow'   # set to lower case for consistency with existing code in optics.py
+                        _log.info(f"Set bar offset to {auto_offset} based on current aperture name {self.aperturename}")
+                    elif apname_last_part.startswith('F'):
+                        auto_offset = apname_last_part
+                        _log.info(f"Set bar offset to {auto_offset} based on current aperture name {self.aperturename}")
+                    else:
+                        auto_offset = self.filter
+                        _log.info(f"Set bar offset to {auto_offset} based on current filter {self.filter}")
+                else:
+                    auto_offset = self.filter
+                    _log.info(f"Set bar offset to {auto_offset} based on current filter {self.filter}")
+
             else:
                 try:
                     _ = float(bar_offset)