diff --git a/docs/relnotes.rst b/docs/relnotes.rst index 1129039b..7e227a32 100644 --- a/docs/relnotes.rst +++ b/docs/relnotes.rst @@ -28,12 +28,14 @@ Version History and Change Log ------------------------------- Version 1.3.0 -============= + *Date TBD* This release comes with new features and improvements including but not limited to: X. Improved support for NIRSpec and MIRI IFU PSF calculations, including addition of a ``mode`` attribute for toggling between imaging mode and IFU mode simulations; an option for much faster (but slightly simplified) IFU datacube PSF calculations. Spectral bandpass information added for the IFU bands for both NIRSpec and MIRI. For NIRSpec, IFU mode PSF outputs are rotated by an additional 90 degrees to match the convention used in pipeline-output s3d datacubes made using the IFUAlign orientation. This extra rotation can be optionally disabled if so desired; see the NIRSpec class docstring. +X. Improved PSF models for MIRI imager, in particular including an empirical model for the field-dependent shifts of the cruciform artifact seen at short wavelengths. + Version 1.2.1 diff --git a/webbpsf/constants.py b/webbpsf/constants.py index 214c88af..27785ba5 100644 --- a/webbpsf/constants.py +++ b/webbpsf/constants.py @@ -400,3 +400,4 @@ # See Gaspar et al. 2021 for illustrative figures. # 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 \ No newline at end of file diff --git a/webbpsf/detectors.py b/webbpsf/detectors.py index b687deaa..e7d04e33 100644 --- a/webbpsf/detectors.py +++ b/webbpsf/detectors.py @@ -284,46 +284,73 @@ def oversample_ipc_model(kernel, oversample): # Functions for applying MIRI Detector Scattering Effect +# Lookup tables of shifts of the cruciform, estimated roughly from F560W ePSFs (ePSFs by Libralatto, shift estimate by Perrin) +cruciform_xshifts = scipy.interpolate.interp1d([0, 357, 1031], [1.5,0.5,-0.9], kind='linear', fill_value='extrapolate') +cruciform_yshifts = scipy.interpolate.interp1d([0, 511, 1031], [1.6,0,-1.6], kind='linear', fill_value='extrapolate') -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). +def _make_miri_scattering_kernel_2d(in_psf, kernel_amp, oversample=1, wavelength=5.5, detector_position=(0, 0)): + """Improved / more complex model of the MIRI cruciform, with parameterization to model + additional features as seen in the GimMIRI models of Gaspar et al. 2021, PASP 133 + See in particular their Figure 12. + + Note, this contains a moderate amount of ad-hoc parameter fitting for scale factors to match observed PSFs from flight. Parameters ---------- - image : ndarray + in_psf : ndarray PSF array for which to make the kernel - amplitude : float - Amplitude of the kernel - nsamples : int + kernel_amp : float + Amplitude scale factor of the kernel + oversample : int Amount by which the input PSF is oversampled + wavelength : float + Wavelength in microns, for use in adding wavelength-dependent effects + detector_position : tuple of floats + X, Y position, for use in adding wavelength-dependent effects - Returns - ------- - kernel_x : ndarray - 1D detector scattering kernel in the x direction """ + # make output array + npix = in_psf.shape[0] + cen = (npix-1) // 2 + kernel_2d = np.zeros( (npix, npix), float) + ### make 1d kernels for the main cruciform bright lines # Compute 1d indices - x = np.arange(image.shape[1], dtype=float) - x -= (image.shape[1] - 1) / 2 - x /= nsamples + x = np.arange(npix, dtype=float) + x -= (npix-1)/2 + x /= oversample + y = x # we're working in 1d in this part, but clarify let's have separate coords for each axis - # Create 1d kernel - kernel_x = amplitude * np.exp(-np.abs(x) / 25) + # Create 1d kernels + kernel_x = kernel_amp * np.exp(-np.abs(x) / 25) + kernel_y = kernel_amp * np.exp(-np.abs(x) / 25) # reduce intensity in the inner part, since the cruciform is suppressed at small radii kernel_x[np.abs(x) < constants.MIRI_CRUCIFORM_INNER_RADIUS_PIX] *= 0.5 + kernel_y[np.abs(y) < constants.MIRI_CRUCIFORM_INNER_RADIUS_PIX] *= 0.5 + + # Add in the offset copies of the main 1d kernels + # Empirically, the 'center' of the cruciform shifts inwards towards the center of the detector + # i.e. for the upper right corner, the cruciform shifts down and left a bit, etc. + yshift = cruciform_yshifts(detector_position[1]) + xshift = cruciform_xshifts(detector_position[0]) + kernel_2d[cen + int(round(yshift*oversample))] = kernel_x + kernel_2d[:, cen + int(round(xshift*oversample))] = kernel_y + + ### create and add in the more diffuse radial term + # Model this as an expoential falloff outside the inner radius, times some scale factor relative to the above + y, x = np.indices(kernel_2d.shape) + r = np.sqrt((x-cen)**2 + (y-cen)**2) / oversample + radial_term = np.exp(-r/2/webbpsf.constants.MIRI_CRUCIFORM_INNER_RADIUS_PIX) * kernel_amp \ + * (r > webbpsf.constants.MIRI_CRUCIFORM_INNER_RADIUS_PIX) \ + * webbpsf.constants.MIRI_CRUCIFORM_RADIAL_SCALEFACTOR - # Reshape kernel to 2D image for use in convolution - kernel_x.shape = (1, image.shape[1]) + kernel_2d += radial_term - return kernel_x + return kernel_2d -def _apply_miri_scattering_kernel(in_psf, kernel_x, oversample): +def _apply_miri_scattering_kernel_2d(in_psf, kernel_2d, oversample): """ Applies the detector scattering kernel created in _make_miri_scattering_kernel function to an input image. Code is adapted from @@ -346,25 +373,65 @@ def _apply_miri_scattering_kernel(in_psf, kernel_x, oversample): 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 - ) + # Convolve the input PSF with the kernel for scattering + im_conv = astropy.convolution.convolve_fft(in_psf, kernel_2d, 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. + # Normalize. # 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) + im_conv_both = im_conv / oversample**2 return im_conv_both -def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None): +def get_miri_cruciform_amplitude(filt): + # 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, + 'FND': 3.0, + } + # FND value is a WAG, interpolating between the F1000W and F1130W values; in reality it varies over that + # huge bandpass, but we can't compute it per-wavelength here. + + # In-flight correction based on measured cycle 1 ePSFs, coarsely + for k in kernel_amp_corrections: + kernel_amp_corrections[k] *= 0.5 + + kernel_amp = kernel_amp_dict[filt] + + if filt in kernel_amp_corrections: + kernel_amp *= kernel_amp_corrections[filt] + return kernel_amp + + +def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None, old_method=False): """ Apply a distortion caused by the MIRI scattering cross artifact effect. In short we convolve a 2D exponentially decaying cross to the PSF where @@ -410,47 +477,9 @@ def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None): 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, - 'FND': 3.0, - } # FND value is a WAG, interpolating between the F1000W and F1130W values; in reality it varies over that huge bandpass, but we can't compute it per-wavelengthhere. - # In-flight correction based on measured cycle 1 ePSFs, coarsely - for k in kernel_amp_corrections: - kernel_amp_corrections[k] *= 0.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] + kernel_amp = get_miri_cruciform_amplitude(filt) ext = 1 # edit the oversampled PSF (OVERDIST extension) @@ -460,11 +489,11 @@ def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None): # 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) + # create cruciform model using improved method using a 2d convolution kernel, attempting to model more physics. + kernel_2d = _make_miri_scattering_kernel_2d(in_psf, kernel_amp, oversample, + detector_position= (hdu_list[0].header['DET_X'], hdu_list[0].header['DET_Y']), + wavelength = hdu_list[0].header['WAVELEN']*1e6 ) + im_conv_both = _apply_miri_scattering_kernel_2d(in_psf, kernel_2d, oversample) # Add this 2D scattered light output to the PSF psf_new = in_psf + im_conv_both @@ -481,3 +510,27 @@ def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None): psf[ext].header['KERNFOLD'] = (25, 'e-folding length (B) in kernel func A*exp(-x/B)') return psf + + +def _show_miri_cruciform_kernel(filt, npix=101, oversample=4, detector_position=(512,512), ax=None): + """ utility function for viewing/visualizing the cruciform kernel + """ + import matplotlib + + placeholder = np.zeros((npix*oversample, npix*oversample)) + kernel_amp = get_miri_cruciform_amplitude(filt) + extent =[-npix/2, npix/2, -npix/2, npix/2] + + kernel_2d = _make_miri_scattering_kernel_2d(placeholder, kernel_amp, oversample, + detector_position= detector_position) + norm = matplotlib.colors.LogNorm(1e-6, 1) + cmap = matplotlib.cm.viridis + cmap.set_bad(cmap(0)) + if ax is None: + ax = matplotlib.pyplot.gca() + ax.imshow(kernel_2d, norm=norm, cmap=cmap, extent=extent, origin='lower') + ax.set_title(f"MIRI cruciform model for {filt}, position {detector_position}, oversample {oversample}") + ax.plot(0,0,marker='+', color='yellow') + + matplotlib.pyplot.colorbar(mappable=ax.images[0]) + diff --git a/webbpsf/tests/test_detectors.py b/webbpsf/tests/test_detectors.py index f8952404..bff7e9fd 100644 --- a/webbpsf/tests/test_detectors.py +++ b/webbpsf/tests/test_detectors.py @@ -60,10 +60,12 @@ def test_apply_miri_scattering(): # Test that the 4 corners of the box contain very small (close to 0) values ylen, xlen = diff.shape - # Choose the start/stop points for these squares (each will take up 1/3 of the total array) + # Choose the start/stop points for these squares (each will take up 1/10th of the total array per side) + # Following the addition of the radial term to the cruciform model, we restrict to just the very corners + # to avoid including much from that circular halo term. first = 0 - second = int(0.33 * xlen) - third = int(0.67 * xlen) + second = int(0.1 * xlen) + third = int(0.9 * xlen) fourth = xlen - 1 # Pull these squares out of the data