diff --git a/webbpsf/constants.py b/webbpsf/constants.py index 4222b398..1ef44ea4 100644 --- a/webbpsf/constants.py +++ b/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}, +} diff --git a/webbpsf/detectors.py b/webbpsf/detectors.py index a3a1e973..63b86fb6 100644 --- a/webbpsf/detectors.py +++ b/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 diff --git a/webbpsf/match_data.py b/webbpsf/match_data.py index 2f0fb098..960a110d 100644 --- a/webbpsf/match_data.py +++ b/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' diff --git a/webbpsf/tests/test_miri.py b/webbpsf/tests/test_miri.py index 12d70fc7..298a7d53 100644 --- a/webbpsf/tests/test_miri.py +++ b/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" diff --git a/webbpsf/tests/test_nirspec.py b/webbpsf/tests/test_nirspec.py index 58c6b7ac..8ba17f69 100644 --- a/webbpsf/tests/test_nirspec.py +++ b/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" diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index 31c895d6..73ab432e 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -1241,7 +1241,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) @@ -1278,19 +1278,20 @@ 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) @@ -1298,6 +1299,7 @@ def _calc_psf_format_output(self, result, options): 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) @@ -2069,10 +2071,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). @@ -2458,7 +2461,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]