moved miri scattering into detector added charge difussion values bas…

…ed on comparisons other small fixes

webbpsf/constants.py

@@ -250,15 +250,16 @@
 # Note, these are parameterized as arcseconds for convenience (and consistency with the jitter paramater)
 # but the underlying physics cares more about detector pixel pitch.
 INSTRUMENT_DETECTOR_CHARGE_DIFFUSION_DEFAULT_PARAMETERS = {
-    'NIRCAM_SW': 0.012,     # Fit by Marcio to WFS TA ePSFs
-    'NIRCAM_LW': 0.024,     # Scaled up by pixel pitch
-    'NIRISS': 0.024,
-    'FGS': 0.024,
+    'NIRCAM_SW': 0.006,     # Fit by Marcio to WFS TA ePSFs
+    'NIRCAM_LW': 0.012,     # Scaled up by pixel pitch
+    'NIRISS': 0.028,         # Fit by Marcio to MIMF-3 F158M (ePSF)
+    'FGS': 0.07,            # Fit by Marcio to FGS_ID images
     'NIRSPEC': 0.036,
-    'MIRI': 0.070,          # Based on user reports, see issue #674
+    'MIRI': 0.05,          # Fit by Marcio to F560W ePSF and single PSF (MIMF-3)
+                           #  0.070 Based on user reports, see issue #674. However, this is before adding IPC effects
 }
 # add Interpixel capacitance (IPC) effects. These are the parameters for each detector kernel
-# For NIRCam we  sue CV#/Flight convolution kernels from Jarron Leisenring, see detectors.apply_detector_ipc for details
+# For NIRCam we  use CV3/Flight convolution kernels from Jarron Leisenring, see detectors.apply_detector_ipc for details
 # NIRISS has different kernels provided by Kevin Volk (STScI), see detectors.apply_detector_ipc for details
 INSTRUMENT_IPC_DEFAULT_KERNEL_PARAMETERS = {
     'MIRI': (0.033, 0.024, 0.013),          # Based on JWST-STScI-002925 by Mike Engesser

webbpsf/detectors.py

@@ -1,17 +1,23 @@
+import copy
+
 import os
 import numpy as np
 import scipy
 import webbpsf
 from webbpsf import utils
 from astropy.convolution.kernels import CustomKernel
 from astropy.convolution import convolve
-import scipy.signal as signal
 from astropy.io import fits
+import astropy.convolution
+import scipy.signal as signal
+
+
+
 
 def apply_detector_ipc(psf_hdulist):
     """Apply a model for interpixel capacitance
     NIRCam: IPC and PPC values derived during ground I&T, primarily ISIM-CV3 from Jarron Leisenring
-    these IPC/PPC kernels will be update it after flight values are available.
+    these IPC/PPC kernels will be update after flight values are available.
     For NIRCam only PPC effects are also included, these are relatively small compared to the IPC contribution
     MIRI: Convolution kernels from JWST-STScI-002925 by Mike Engesser
     NIRISS: Convolution kernels and base code provided by Kevin Volk
@@ -111,12 +117,17 @@ def apply_detector_ipc(psf_hdulist):
         else:
             newimage = numpy.copy(image)
             psf_hdulist[ext].header['IPCINST'] = ('NIRISS', 'Interpixel capacitance (IPC)')
-            psf_hdulist[ext].header['IPCTYPA'] = ('NIRISS', 'no kernel file found')
-            sf_hdulist[ext].header['IPCTYPB'] = ('NIRISS', 'no IPC correction applied')
+            psf_hdulist[ext].header['IPCTYPA'] = ('NIRISS', 'No kernel file found')
+            sf_hdulist[ext].header['IPCTYPB'] = ('NIRISS', 'No IPC correction applied')
             webbpsf.webbpsf_core._log.info("No IPC correction for NIRISS. Check kernel files.")
 
         psf_hdulist[ext].data = newimage
 
+    if inst in ["FGS", "NIRSpec"]:
+        ext = 'DET_DIST'
+        psf_hdulist[ext].header['IPCINST'] = (inst, 'No IPC correction applied')
+        webbpsf.webbpsf_core._log.info("IPC corrections are not implemented yet for {}".format(inst))
+
 
     return psf_hdulist
 
@@ -144,3 +155,175 @@ def apply_detector_charge_diffusion(psf_hdulist, options):
     psf_hdulist[ext].data = out
 
     return psf_hdulist
+
+# Functions for applying MIRI Detector Scattering Effect
+
+def _make_miri_scattering_kernel(image, amplitude, nsamples):
+    """
+    Creates a detector scatter kernel function. For simplicity, we assume a
+    simple exponential dependence. Code is adapted from
+    MIRI-TN-00076-ATC_Imager_PSF_Issue_4.pdf (originally in IDL).
+
+    Parameters
+    ----------
+    image : ndarray
+        PSF array for which to make the kernel
+    amplitude : float
+        Amplitude of the kernel
+    nsamples : int
+        Amount by which the input PSF is oversampled
+
+    Returns
+    -------
+    kernel_x : ndarray
+        1D detector scattering kernel in the x direction
+    """
+
+    # Compute 1d indices
+    x = np.arange(image.shape[1], dtype=float)
+    x -= (image.shape[1]-1)/2
+    x /= nsamples
+
+    # Create 1d kernel
+    kernel_x = amplitude * np.exp(-np.abs(x) / 25)
+
+    # Reshape kernel to 2D image for use in convolution
+    kernel_x.shape = (1, image.shape[1])
+
+    return kernel_x
+
+
+def _apply_miri_scattering_kernel(in_psf, kernel_x, oversample):
+    """
+    Applies the detector scattering kernel created in _make_miri_scattering_kernel
+    function to an input image. Code is adapted from
+    MIRI-TN-00076-ATC_Imager_PSF_Issue_4.pdf
+
+    Parameters
+    ----------
+    in_psf : ndarray
+        PSF array upon which to apply the kernel
+    kernel_x : ndarray
+        The 1D kernel in the x direction, output from _make_miri_scattering_kernel.
+        This will be transposed to create the kernel in the y direction.
+    oversample : int
+        Amount by which the input PSF is oversampled
+
+    Returns
+    -------
+    im_conv_both : ndarray
+        The input image convolved with the input kernel in both the x and
+        y directions
+    """
+
+    # Apply the kernel via convolution in both the X and Y direction
+    # Convolve the input PSF with the kernel for scattering in the X direction
+    im_conv_x = astropy.convolution.convolve_fft(in_psf, kernel_x, boundary='fill', fill_value=0.0,
+                                                 normalize_kernel=False, nan_treatment='fill', allow_huge = True)
+
+    # Transpose to make a kernel for Y and convolve with that too
+    im_conv_y = astropy.convolution.convolve_fft(in_psf, kernel_x.T, boundary='fill', fill_value=0.0,
+                                                 normalize_kernel=False, nan_treatment='fill', allow_huge = True)
+
+    # Sum together both the X and Y scattering.
+    # Note, it appears we do need to correct the amplitude for the sampling factor. Might as well do that here.
+    im_conv_both = (im_conv_x + im_conv_y)/(oversample**2)
+
+    return im_conv_both
+
+
+def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None):
+    """
+    Apply a distortion caused by the MIRI scattering cross artifact effect.
+    In short we convolve a 2D exponentially decaying cross to the PSF where
+    the amplitude of the exponential function is determined by the filter of
+    the PSF. A full description of the distortion and the original code can
+    be found in MIRI-TN-00076-ATC_Imager_PSF_Issue_4.pdf
+
+    Note, this code **edits in place Extension 1 of the supplied HDUlist**. In the typical case where the
+    input PSF is calculated as Extension 0, the calling function must put a copy of that into Extension 1
+    which this will then modify. This happens in webbpsf_core.py/JWInstrument._calc_psf_format_output,
+    which is where this is called from in the usual course of operation.
+
+    Parameters
+    ----------
+    hdulist_or_filename :
+        A PSF from WebbPSF, either as an HDUlist object or as a filename
+    kernel_amp: float
+        Detector scattering kernel amplitude. If set to None,
+        function will pull the value based on best fit analysis
+        using the input PSF's filter. Default = None.
+
+    Returns
+    -------
+    psf : HDUlist object
+        PSF with MIRI detector scattering effect applied
+    """
+
+    # Read in input PSF
+    if isinstance(hdulist_or_filename, str):
+        hdu_list = fits.open(hdulist_or_filename)
+    elif isinstance(hdulist_or_filename, fits.HDUList):
+        hdu_list = hdulist_or_filename
+    else:
+        raise ValueError("input must be a filename or HDUlist")
+
+    # Create a copy of the PSF
+    psf = copy.deepcopy(hdu_list)
+
+    # Log instrument name and filter
+    instrument = hdu_list[0].header["INSTRUME"].upper()
+    filt = hdu_list[0].header["FILTER"].upper()
+
+    if instrument != "MIRI":
+        raise ValueError("MIRI's Scattering Effect should only be applied to MIRI PSFs")
+
+    # Default kernel amplitude values from modeling in MIRI-TN-00076-ATC_Imager_PSF_Issue_4.pdf
+    kernel_amp_dict = {'F560W': 0.00220, 'F770W': 0.00139, 'F1000W': 0.00034,
+                       'F1130W': 0.00007, 'F1280W': 0.00011, 'F1500W': 0.0,
+                       'F1800W': 0.0, 'F2100W': 0.0, 'F2550W': 0.0, 'FND': 0.00087,
+                       'F1065C': 0.00010, 'F1140C': 0.00007, 'F1550C': 0.0,
+                       'F2300C': 0.0}
+
+    # The above values are from that tech report, but empirically we need higher values to
+    # better match the MIRI CDP PSFS. See e.g. MIRI_FM_MIRIMAGE_F560W_PSF_07.02.00.fits
+    # and https://github.com/spacetelescope/webbpsf/issues/415
+    kernel_amp_corrections = {'F560W': 4.05, 'F770W': 4.1, 'F1000W': 3.8,
+                               'F1130W': 2.5, 'F1280W': 2.5, 'F1065C': 2.5, 'F1140C': 2.5}
+
+    # Set values if not already set by a keyword argument
+    if kernel_amp is None:
+        kernel_amp = kernel_amp_dict[filt]
+
+        if filt in kernel_amp_corrections:
+            kernel_amp *= kernel_amp_corrections[filt]
+
+    ext = 1  # edit the oversampled PSF (OVERDIST extension)
+
+    # Set over-sample value
+    oversample = psf[ext].header["DET_SAMP"]
+
+    # Read in PSF
+    in_psf = psf[ext].data
+
+    # Make the kernel
+    kernel_x = _make_miri_scattering_kernel(in_psf, kernel_amp, oversample)
+
+    # Apply the kernel via convolution in both the X and Y direction to produce a 2D output
+    im_conv_both = _apply_miri_scattering_kernel(in_psf, kernel_x, oversample)
+
+    # Add this 2D scattered light output to the PSF
+    psf_new = in_psf + im_conv_both
+
+    # To ensure conservation of intensity, normalize the psf
+    psf_new *= in_psf.sum() / psf_new.sum()
+
+    # Apply data to correct extensions
+    psf[ext].data = psf_new
+
+    # Set new header keywords
+    psf[ext].header["MIR_DIST"] = ("True", "MIRI detector scattering applied")
+    psf[ext].header["KERN_AMP"] = (kernel_amp, "Amplitude (A) in kernel function A*exp(-x/B)")
+    psf[ext].header["KERNFOLD"] = (25, "e-folding length (B) in kernel func A*exp(-x/B)")
+
+    return psf