Skip to content

Commit

Permalink
Merge pull request #837 from mperrin/better_miri_cruciform
Browse files Browse the repository at this point in the history 
Improved model for MIRI cruciform artifact, part 1
  • Loading branch information
@mperrin
mperrin authored May 17, 2024
2 parents 6cf62b0 + a8889bf commit 0957cd7
Show file tree
Hide file tree
Showing 4 changed files with 140 additions and 82 deletions.
4 changes: 3 additions & 1 deletion docs/relnotes.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions webbpsf/constants.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
209 changes: 131 additions & 78 deletions webbpsf/detectors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand All @@ -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
Expand All @@ -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])

8 changes: 5 additions & 3 deletions webbpsf/tests/test_detectors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 0957cd7

Please sign in to comment.