Merge pull request #890 from mperrin/ifu_mode_improvements_part2

IFU mode improvements, continued

webbpsf/constants.py

@@ -401,3 +401,9 @@
 # This is a rough approximation of a detector-position-dependent phenomenon
 MIRI_CRUCIFORM_INNER_RADIUS_PIX = 12
 MIRI_CRUCIFORM_RADIAL_SCALEFACTOR = 0.005   # Brightness factor for the diffuse circular halo
+
+# Parameters for adjusting models of IFU PSFs relative to regular imaging PSFs
+INSTRUMENT_IFU_BROADENING_PARAMETERS = {
+    'NIRSPEC': {'sigma': 0.05},
+    'MIRI': {'sigma': 0.05},
+}

webbpsf/detectors.py

@@ -8,6 +8,8 @@
 from astropy.convolution import convolve
 from astropy.convolution.kernels import CustomKernel
 from astropy.io import fits
+from astropy import units as u
+from astropy.modeling.functional_models import Gaussian1D, GAUSSIAN_SIGMA_TO_FWHM
 
 import webbpsf
 from webbpsf import constants, utils
@@ -195,15 +197,18 @@ def apply_detector_ipc(psf_hdulist, extname='DET_DIST'):
 
             out_ipc_0 = convolve(psf_hdulist[extname].data, ipckernel)
             out_ipc = convolve(out_ipc_0, ppckernel)
+            webbpsf.webbpsf_core._log.info(f'ext {extname}: Added IPC and PPC models for {inst}')
         elif inst.upper() == 'NIRISS':
             # the NIRISS code provided by Kevin Volk was developed for a different convolution function
             if oversample != 1:
                 kernel = oversample_ipc_model(kernel, oversample)
             out_ipc = signal.fftconvolve(psf_hdulist[extname].data, kernel, mode='same')
+            webbpsf.webbpsf_core._log.info(f'ext {extname}: Added IPC model for {inst}')
         else:
             if oversample != 1:
                 kernel = oversample_ipc_model(kernel, oversample)
             out_ipc = convolve(psf_hdulist[extname].data, kernel)
+            webbpsf.webbpsf_core._log.info(f'ext {extname}: Added IPC model for {inst}')
 
         # apply kernel to DET_DIST
         psf_hdulist[extname].data = out_ipc
@@ -214,7 +219,7 @@ def apply_detector_ipc(psf_hdulist, extname='DET_DIST'):
         psf_hdulist[extname].header.add_history('Applied detector interpixel capacitance (IPC) model')
 
     else:
-        webbpsf.webbpsf_core._log.info('IPC corrections are not implemented yet for {}'.format(inst))
+        webbpsf.webbpsf_core._log.info('IPC model not implemented yet for {}'.format(inst))
         psf_hdulist[extname].header['IPCINST'] = (inst, 'No IPC correction applied')
 
     return psf_hdulist
@@ -542,3 +547,129 @@ def _show_miri_cruciform_kernel(filt, npix=101, oversample=4, detector_position=
     ax.plot(0, 0, marker='+', color='yellow')
 
     matplotlib.pyplot.colorbar(mappable=ax.images[0])
+
+# Functions for applying IFU optics systematics models
+#
+# Note, this is not actually a "Detector" effect, but this file is a
+# convenient place to locate that code, because similar to the detector effects
+# it's implemented as a post-processing modification on the output PSF array.
+
+
+def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196):
+    """ Apply a simple empirical model of MIRI IFU broadening to better match observed PSFs
+
+    Parameters
+    -----------
+    hdulist : astropy.io.fits.HDUList
+		PSF calculation output data structure. Will be modified.
+	options : dict
+		Options dict for setting optional behaviors
+    slice_width : float
+		MIRI MRS IFU slice width (across the slice). See MIRI._IFU_pixelscale in webbpsf_core.py
+
+    """
+    # First, check an optional flag to see whether or not to include this effect.
+    # User can set the option to None to disable this step.
+    model_type = options.get('ifu_broadening', 'empirical')
+
+    if model_type is None or model_type.lower() == 'none':
+        return hdulist
+
+    ext = 1  # Apply this effect to the OVERDIST extension, which at this point in the code will be ext 1
+
+    webbpsf.webbpsf_core._log.info(f'Applying MIRI IFU broadening model: {model_type}')
+    hdulist[ext].header.add_history(f"Added broadening model for IFU PSFs: {model_type}")
+
+    hdulist[ext].header['IFUBROAD'] = (True, "IFU PSF broadening model applied")
+    hdulist[ext].header['IFUBTYPE'] = (model_type, "IFU PSF broadening model type")
+
+    if model_type.lower() == 'gaussian':
+        # Very simple model just as a Gaussian convolution kernel
+        sigma = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS['MIRI']['sigma']
+        hdulist[ext].header['IFUBSIGM'] = (sigma, "[arcsec] IFU PSF broadening Gaussian sigma")
+        out = scipy.ndimage.gaussian_filter(hdulist[ext].data, sigma / hdulist[ext].header['PIXELSCL'])
+    elif model_type.lower() == 'empirical':
+        # Model based on empirical PSF properties, Argryiou et al.
+        pixelscl = float(hdulist[ext].header['PIXELSCL'])
+        wavelen = float(hdulist[ext].header['WAVELEN'])
+
+        beta_width = slice_width / pixelscl
+        alpha_width = _miri_mrs_analytical_sigma_alpha_broadening(wavelen * 1e6) / pixelscl
+        out = _miri_mrs_empirical_broadening(psf_model=hdulist[ext].data, alpha_width=alpha_width, beta_width=beta_width)
+
+    hdulist[ext].data = out
+
+    return hdulist
+
+
+def apply_nirspec_ifu_broadening(hdulist, options):
+    """ Apply a simple empirical model of NIRSpec IFU broadening to better match observed PSFs
+
+    """
+    # First, check an optional flag to see whether or not to include this effect
+    model_type = options.get('ifu_broadening', 'gaussian')
+    if model_type is None or model_type.lower() == 'none':
+        return hdulist
+
+    ext = 1  # Apply this effect to the OVERDIST extension, which at this point in the code will be ext 1
+
+    webbpsf.webbpsf_core._log.info(f'Applying NRS IFU broadening model ({model_type}) to '+
+                                   f'ext {hdulist[ext].header["EXTNAME"]}')
+
+    hdulist[ext].header['IFUBROAD'] = (True, "IFU PSF broadening model applied")
+    hdulist[ext].header['IFUBTYPE'] = (model_type, "IFU PSF broadening model type")
+    hdulist[ext].header.add_history(f"Added broadening model for IFU PSFs: {model_type}")
+
+    if model_type.lower() == 'gaussian':
+        sigma = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS['NIRSPEC']['sigma']
+        # currently sigma= 50 mas, half a NIRSpec IFU spaxel. Approximate and loose estimate
+        hdulist[ext].header['IFUBSIGM'] = (sigma, "[arcsec] IFU PSF broadening Gaussian sigma")
+        out = scipy.ndimage.gaussian_filter(hdulist[ext].data, sigma / hdulist[ext].header['PIXELSCL'])
+
+    hdulist[ext].data = out
+
+    return hdulist
+
+
+def _miri_mrs_analytical_sigma_alpha_broadening(wavelength):
+    """
+    Calculate the Gaussian scale of the kernel that broadens the diffraction limited
+    FWHM to the empirically measured FWHM.
+
+    Parameters
+    ----------
+    wavelength : float or ndarray
+		wavelength in MICRONS
+    """
+    empirical_fwhm = 0.033 * wavelength + 0.106  # Law+2023
+    diffraction_fwhm = astropy.coordinates.Angle(1.025*wavelength*1E-6/constants.JWST_CIRCUMSCRIBED_DIAMETER, u.radian).to_value(u.arcsec)
+
+    sigma_emp = empirical_fwhm/GAUSSIAN_SIGMA_TO_FWHM
+    sigma_diffr = diffraction_fwhm/GAUSSIAN_SIGMA_TO_FWHM
+    return np.sqrt(sigma_emp**2 - sigma_diffr**2)  # return kernel width in arcsec
+
+
+def _miri_mrs_empirical_broadening(psf_model, alpha_width, beta_width):
+     """
+     Perform the broadening of a psf model in alpha and beta
+
+     Parameters
+     -----------
+     psf_model : ndarray
+        webbpsf output results, eitehr monochromatic model or datacube
+     alpha_width : float
+        gaussian convolution kernel in pixels, None if no broadening should be performed
+     beta_width : float
+        slice width in pixels
+     """
+     kernel_beta = astropy.convolution.Box1DKernel(beta_width)
+
+    # TODO: extend algorithm to handle the datacube case
+
+     if alpha_width is None:
+         psf_model_alpha_beta = np.apply_along_axis(lambda m: convolve(m, kernel_beta), axis=0, arr=psf_model)
+     else:
+         kernel_alpha = astropy.convolution.Gaussian1DKernel(stddev=alpha_width)
+         psf_model_alpha = np.apply_along_axis(lambda m: convolve(m, kernel_alpha), axis=1, arr=psf_model)
+         psf_model_alpha_beta = np.apply_along_axis(lambda m: convolve(m, kernel_beta), axis=0, arr=psf_model_alpha)
+     return psf_model_alpha_beta

webbpsf/match_data.py

@@ -31,7 +31,11 @@ def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choic
 
     inst = webbpsf.instrument(header['INSTRUME'])
 
-    if inst.name == 'MIRI' and header['FILTER'] == 'P750L':
+    if inst.name == 'MIRI' and header['EXP_TYPE'] == 'MIR_MRS':
+        print("MIRI MRS exposure detected; configuring for IFU mode")
+        inst.mode = 'IFU'
+        # There is no FILTER keyword for MRS, so don't set filter to anything.
+    elif inst.name == 'MIRI' and header['FILTER'] == 'P750L':
         # 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'
@@ -68,7 +72,12 @@ def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choic
             inst.pupil_mask = header['PUPIL']
 
     elif inst.name == 'MIRI':
-        if inst.filter in ['F1065C', 'F1140C', 'F1550C']:
+        if header['EXP_TYPE'] == 'MIR_MRS':
+            ch = header['CHANNEL']
+            band_lookup = {'SHORT': 'A', 'MEDIUM': 'B', 'LONG': 'C'}
+            inst.band = str(ch) + band_lookup[header['BAND']]
+
+        elif inst.filter in ['F1065C', 'F1140C', 'F1550C']:
             inst.image_mask = 'FQPM' + inst.filter[1:5]
         elif inst.filter == 'F2300C':
             inst.image_mask = 'LYOT2300'

webbpsf/tests/test_miri.py

@@ -5,6 +5,7 @@
 import numpy as np
 import pysiaf
 
+import webbpsf
 from .. import webbpsf_core
 from .test_webbpsf import do_test_set_position_from_siaf, do_test_source_offset, generic_output_test
 
@@ -224,3 +225,29 @@ def test_IFU_wavelengths():
     # and test we can specify a reduced wavelength sampling:
     for n in (10, 100):
         assert len(miri.get_IFU_wavelengths(n)) == n
+
+
+def test_miri_ifu_broadening():
+    """ Basic functional test for the code that adjusts PSF outputs to better match empirical IFU PSFs
+    """
+
+    miri = webbpsf_core.MIRI()
+    miri.mode = 'IFU'
+    psf = miri.calc_psf(monochromatic=2.8e-6, fov_pixels=10)
+
+    fwhm_oversamp = webbpsf.measure_fwhm(psf, ext='OVERSAMP')
+    fwhm_overdist = webbpsf.measure_fwhm(psf, ext='OVERDIST')
+    assert fwhm_overdist > fwhm_oversamp, "IFU broadening model should increase the FWHM for the distorted extensions"
+
+    fwhm_detsamp = webbpsf.measure_fwhm(psf, ext='DET_SAMP')
+    fwhm_detdist = webbpsf.measure_fwhm(psf, ext='DET_DIST')
+    assert fwhm_overdist > fwhm_oversamp, "IFU broadening model should increase the FWHM for the distorted extensions"
+
+    # Now test that we can also optionally turn off that effect
+    miri.options['ifu_broadening'] = None
+    psf_nb = miri.calc_psf(monochromatic=2.8e-6, fov_pixels=10)
+
+    fwhm_oversamp = webbpsf.measure_fwhm(psf_nb, ext='OVERSAMP')
+    fwhm_overdist = webbpsf.measure_fwhm(psf_nb, ext='OVERDIST')
+    # The PSF will still be a little broader in this case due to the IPC model, but not by a lot..
+    assert fwhm_oversamp < fwhm_overdist <= 1.1 * fwhm_oversamp, "IFU broadening model should be disabled for this test case"

webbpsf/tests/test_nirspec.py

@@ -4,6 +4,7 @@
 import numpy as np
 import pysiaf
 
+import webbpsf
 from .. import webbpsf_core
 from .test_webbpsf import do_test_set_position_from_siaf, do_test_source_offset, generic_output_test
 
@@ -96,6 +97,7 @@ def test_IFU_wavelengths():
         assert len(nrs.get_IFU_wavelengths(n)) == n
 
 
+
 def test_aperture_rotation_updates():
     """ Test that switching detector or aperture updates the aperture PA
     (this is a tiny detail)"""
@@ -114,3 +116,26 @@ def test_aperture_rotation_updates():
     nrs.detector = 'NRS2'
     assert nrs.aperturename == 'NRS2_FULL'
     assert nrs._rotation != pa_nrs1_full
+
+def test_nrs_ifu_broadening():
+    """ Basic functional test for the code that adjusts PSF outputs to better match empirical IFU PSFs
+    """
+    nrs = webbpsf_core.NIRSpec()
+    nrs.mode = 'IFU'
+    psf = nrs.calc_psf(monochromatic=2.8e-6, fov_pixels=10)
+
+    fwhm_oversamp = webbpsf.measure_fwhm(psf, ext='OVERSAMP')
+    fwhm_overdist = webbpsf.measure_fwhm(psf, ext='OVERDIST')
+    assert fwhm_overdist > fwhm_oversamp, "IFU broadening model should increase the FWHM for the distorted extensions"
+
+    fwhm_detsamp = webbpsf.measure_fwhm(psf, ext='DET_SAMP')
+    fwhm_detdist = webbpsf.measure_fwhm(psf, ext='DET_DIST')
+    assert fwhm_overdist > fwhm_oversamp, "IFU broadening model should increase the FWHM for the distorted extensions"
+
+    # Now test that we can also optionally turn off that effect
+    nrs.options['ifu_broadening'] = None
+    psf_nb = nrs.calc_psf(monochromatic=2.8e-6, fov_pixels=10)
+
+    fwhm_oversamp = webbpsf.measure_fwhm(psf_nb, ext='OVERSAMP')
+    fwhm_overdist = webbpsf.measure_fwhm(psf_nb, ext='OVERDIST')
+    assert fwhm_overdist == fwhm_oversamp, "IFU broadening model should be disabled for this test case"

webbpsf/webbpsf_core.py

@@ -1256,7 +1256,7 @@ def _calc_psf_format_output(self, result, options):
 
         # Add distortion if set in calc_psf
         if add_distortion:
-            _log.debug('Adding PSF distortion(s) and detector effects')
+            _log.info('Adding PSF distortion(s) and detector effects')
 
             # Set up new extensions to add distortion to:
             n_exts = len(result)
@@ -1293,26 +1293,28 @@ def _calc_psf_format_output(self, result, options):
                 else:
                     # there is not yet any distortion calibration for the IFU, and
                     # we don't want to apply charge diffusion directly here
-                    psf_distorted = result
+                    psf_distorted = detectors.apply_miri_ifu_broadening(result, options, slice_width=self._ifu_slice_width)
             elif self.name == 'NIRSpec':
                 # Apply distortion effects to NIRSpec psf: Distortion only
                 # (because applying detector effects would only make sense after simulating spectral dispersion)
                 _log.debug('NIRSpec: Adding optical distortion')
-                if 'IFU' not in self.aperturename:
+                if self.mode != 'IFU':
                     psf_siaf = distortion.apply_distortion(result)  # apply siaf distortion model
+                    psf_distorted = detectors.apply_detector_charge_diffusion(
+                        psf_siaf, options
+                    )  # apply detector charge transfer model
+
                 else:
                     # there is not yet any distortion calibration for the IFU.
-                    psf_siaf = result
-                psf_distorted = detectors.apply_detector_charge_diffusion(
-                    psf_siaf, options
-                )  # apply detector charge transfer model
+                    psf_distorted = detectors.apply_nirspec_ifu_broadening(result, options)
 
             # Edit the variable to match if input didn't request distortion
             # (cannot set result = psf_distorted due to return method)
             [result.append(fits.ImageHDU()) for i in np.arange(len(psf_distorted) - len(result))]
             for ext in np.arange(len(psf_distorted)):
                 result[ext] = psf_distorted[ext]
 
+        _log.info('Formatting output FITS extensions including for sampling.')
         # Rewrite result variable based on output_mode; this includes binning down to detector sampling.
         SpaceTelescopeInstrument._calc_psf_format_output(self, result, options)
 
@@ -2103,10 +2105,11 @@ def __init__(self):
 
         self.monochromatic = 8.0
         self._IFU_pixelscale = {
-            'Ch1': (0.176, 0.196),
-            'Ch2': (0.277, 0.196),
-            'Ch3': (0.387, 0.245),
-            'Ch4': (0.645, 0.273),
+                  # slice width, pixel size.   Values from Argyriou et al. 2023 A&A 675
+            'Ch1': (0.177, 0.196),
+            'Ch2': (0.280, 0.196),
+            'Ch3': (0.390, 0.245),
+            'Ch4': (0.656, 0.273),
         }
         # The above tuples give the pixel resolution (first the 'alpha' direction, perpendicular to the slice,
         # then the 'beta' direction, along the slice).
@@ -2492,7 +2495,7 @@ def band(self, value):
 
         if value in self._IFU_bands_cubepars.keys():
             self._band = value
-            # self._slice_width = self._IFU_pixelscale[f"Ch{self._band[0]}"][0]
+            self._ifu_slice_width = self._IFU_pixelscale[f"Ch{self._band[0]}"][0]
             self.aperturename = 'MIRIFU_CHANNEL' + value
             # setting aperturename will also auto update self._rotation
             # self._rotation = self.MRS_rotation[self._band]