From b499e1b9a7877c7fe8fa9eaa99928c40184bf68c Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Fri, 26 Jul 2024 13:05:54 -0400 Subject: [PATCH 1/7] start implementation of framework for IFU broadening models --- webbpsf/detectors.py | 62 +++++++++++++++++++++++++++++++++++++++++ webbpsf/webbpsf_core.py | 13 +++++---- 2 files changed, 69 insertions(+), 6 deletions(-) diff --git a/webbpsf/detectors.py b/webbpsf/detectors.py index a3a1e973..2b9b92b0 100644 --- a/webbpsf/detectors.py +++ b/webbpsf/detectors.py @@ -542,3 +542,65 @@ 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, thes are 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): + """ Apply a simple empirical model of MIRI IFU broadening to better match observed PSFs + + """ + # First, check an optional flag to see whether or not to include this effect + perform_this_step = options.get('ifu_broadening', True) + if not perform_this_step: + return hdulist + + ext = 1 # Apply this effect to the OVERDIST extension, which at this point in the code will be ext 1 + + model_type= options.get('ifu_broadening_type', 'gaussian') + webbpsf.webbpsf_core._log.info(f'Applying MIRI IFU broadening model: {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': + sigma = 0.05 # 50 mas, half a NIRSpec IFU spaxel. Approximate and loose estimate + hdulist[ext].header['IFUBSIGM'] = (model_type, "[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 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 + perform_this_step = options.get('ifu_broadening', True) + if not perform_this_step: + return hdulist + + + ext = 1 # Apply this effect to the OVERDIST extension, which at this point in the code will be ext 1 + + model_type= options.get('ifu_broadening_type', 'gaussian') + webbpsf.webbpsf_core._log.info(f'Applying NRS IFU broadening model: {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': + sigma = 0.05 # 50 mas, half a NIRSpec IFU spaxel. Approximate and loose estimate + hdulist[ext].header['IFUBSIGM'] = (model_type, "[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 diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index 292f6f1c..cc06e1da 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -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) 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) From 737835f2058605e021b6ecefacadd472f1371a01 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Fri, 26 Jul 2024 13:06:17 -0400 Subject: [PATCH 2/7] improve setup_sim_to_match_file for MIRI MRS files --- webbpsf/match_data.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/webbpsf/match_data.py b/webbpsf/match_data.py index 8b8e372a..8c74cb53 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' From aad77945f11fb1b8c3c29f61cdba8d8343a89e53 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Fri, 26 Jul 2024 14:19:18 -0400 Subject: [PATCH 3/7] continue IFU effects implementation; add tests; tweak and improve log messages --- webbpsf/constants.py | 6 ++++++ webbpsf/detectors.py | 27 +++++++++++++++------------ webbpsf/tests/test_miri.py | 28 ++++++++++++++++++++++++++++ webbpsf/tests/test_nirspec.py | 25 +++++++++++++++++++++++++ webbpsf/webbpsf_core.py | 3 ++- 5 files changed, 76 insertions(+), 13 deletions(-) diff --git a/webbpsf/constants.py b/webbpsf/constants.py index 4222b398..6d18cd2c 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 2b9b92b0..77308ebb 100644 --- a/webbpsf/detectors.py +++ b/webbpsf/detectors.py @@ -195,15 +195,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 +217,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 @@ -554,21 +557,22 @@ def apply_miri_ifu_broadening(hdulist, options): """ Apply a simple empirical model of MIRI IFU broadening to better match observed PSFs """ - # First, check an optional flag to see whether or not to include this effect - perform_this_step = options.get('ifu_broadening', True) - if not perform_this_step: + # 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', '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 - model_type= options.get('ifu_broadening_type', 'gaussian') webbpsf.webbpsf_core._log.info(f'Applying MIRI IFU broadening model: {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': - sigma = 0.05 # 50 mas, half a NIRSpec IFU spaxel. Approximate and loose estimate + sigma = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS['MIRI']['sigma'] hdulist[ext].header['IFUBSIGM'] = (model_type, "[arcsec] IFU PSF broadening Gaussian sigma") out = scipy.ndimage.gaussian_filter(hdulist[ext].data, sigma / hdulist[ext].header['PIXELSCL']) @@ -582,21 +586,20 @@ def apply_nirspec_ifu_broadening(hdulist, options): """ # First, check an optional flag to see whether or not to include this effect - perform_this_step = options.get('ifu_broadening', True) - if not perform_this_step: + 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 - model_type= options.get('ifu_broadening_type', 'gaussian') - webbpsf.webbpsf_core._log.info(f'Applying NRS IFU broadening model: {model_type}') + webbpsf.webbpsf_core._log.info(f'Applying NRS IFU broadening model ({model_type}) to 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") if model_type.lower() == 'gaussian': - sigma = 0.05 # 50 mas, half a NIRSpec IFU spaxel. Approximate and loose estimate + 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'] = (model_type, "[arcsec] IFU PSF broadening Gaussian sigma") out = scipy.ndimage.gaussian_filter(hdulist[ext].data, sigma / hdulist[ext].header['PIXELSCL']) diff --git a/webbpsf/tests/test_miri.py b/webbpsf/tests/test_miri.py index 12d70fc7..4774a40d 100644 --- a/webbpsf/tests/test_miri.py +++ b/webbpsf/tests/test_miri.py @@ -224,3 +224,31 @@ 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.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" + print(fwhm_oversamp, fwhm_overdist) + + 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') + print(fwhm_oversamp, fwhm_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 53d18ad4..34790f0f 100644 --- a/webbpsf/tests/test_nirspec.py +++ b/webbpsf/tests/test_nirspec.py @@ -91,3 +91,28 @@ def test_IFU_wavelengths(): # and test we can specify a reduced wavelength sampling: for n in (10, 100): assert len(nrs.get_IFU_wavelengths(n)) == n + + +def test_nrs_ifu_broadening(): + """ Basic functional test for the code that adjusts PSF outputs to better match empirical IFU PSFs + """ + nrs = webbpsf.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 cc06e1da..3df10bc4 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) @@ -1299,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) From 5f3aa2957fbfc745ff48c6f4b5264c54fda3a100 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Fri, 26 Jul 2024 14:36:20 -0400 Subject: [PATCH 4/7] code lint cleanup --- webbpsf/constants.py | 2 +- webbpsf/detectors.py | 8 ++++---- webbpsf/match_data.py | 4 ++-- webbpsf/tests/test_miri.py | 6 ++---- 4 files changed, 9 insertions(+), 11 deletions(-) diff --git a/webbpsf/constants.py b/webbpsf/constants.py index 6d18cd2c..1ef44ea4 100644 --- a/webbpsf/constants.py +++ b/webbpsf/constants.py @@ -406,4 +406,4 @@ INSTRUMENT_IFU_BROADENING_PARAMETERS = { 'NIRSPEC': {'sigma': 0.05}, 'MIRI': {'sigma': 0.05}, - } +} diff --git a/webbpsf/detectors.py b/webbpsf/detectors.py index 77308ebb..6d698ea8 100644 --- a/webbpsf/detectors.py +++ b/webbpsf/detectors.py @@ -564,7 +564,7 @@ def apply_miri_ifu_broadening(hdulist, options): 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 + 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}') @@ -590,9 +590,10 @@ def apply_nirspec_ifu_broadening(hdulist, options): 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 + 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 ext {hdulist[ext].header["EXTNAME"]}') + 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") @@ -605,5 +606,4 @@ def apply_nirspec_ifu_broadening(hdulist, options): hdulist[ext].data = out - return hdulist diff --git a/webbpsf/match_data.py b/webbpsf/match_data.py index 8c74cb53..ea983d03 100644 --- a/webbpsf/match_data.py +++ b/webbpsf/match_data.py @@ -72,9 +72,9 @@ def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choic inst.pupil_mask = header['PUPIL'] elif inst.name == 'MIRI': - if header['EXP_TYPE']=='MIR_MRS': + if header['EXP_TYPE'] == 'MIR_MRS': ch = header['CHANNEL'] - band_lookup = {'SHORT':'A', 'MEDIUM': 'B', 'LONG': 'C'} + band_lookup = {'SHORT': 'A', 'MEDIUM': 'B', 'LONG': 'C'} inst.band = str(ch) + band_lookup[header['BAND']] elif inst.filter in ['F1065C', 'F1140C', 'F1550C']: diff --git a/webbpsf/tests/test_miri.py b/webbpsf/tests/test_miri.py index 4774a40d..175ae6e1 100644 --- a/webbpsf/tests/test_miri.py +++ b/webbpsf/tests/test_miri.py @@ -225,6 +225,7 @@ def test_IFU_wavelengths(): 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 """ @@ -236,7 +237,6 @@ def test_miri_ifu_broadening(): 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" - print(fwhm_oversamp, fwhm_overdist) fwhm_detsamp = webbpsf.measure_fwhm(psf, ext='DET_SAMP') fwhm_detdist = webbpsf.measure_fwhm(psf, ext='DET_DIST') @@ -248,7 +248,5 @@ def test_miri_ifu_broadening(): fwhm_oversamp = webbpsf.measure_fwhm(psf_nb, ext='OVERSAMP') fwhm_overdist = webbpsf.measure_fwhm(psf_nb, ext='OVERDIST') - print(fwhm_oversamp, fwhm_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" - + assert fwhm_oversamp < fwhm_overdist <= 1.1 * fwhm_oversamp, "IFU broadening model should be disabled for this test case" From 04e05b9f6a9e43126d92dd27899f1aa3b3f6833f Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Fri, 26 Jul 2024 15:55:46 -0400 Subject: [PATCH 5/7] fix tests --- webbpsf/tests/test_miri.py | 3 ++- webbpsf/tests/test_nirspec.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/webbpsf/tests/test_miri.py b/webbpsf/tests/test_miri.py index 175ae6e1..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 @@ -230,7 +231,7 @@ def test_miri_ifu_broadening(): """ Basic functional test for the code that adjusts PSF outputs to better match empirical IFU PSFs """ - miri = webbpsf.MIRI() + miri = webbpsf_core.MIRI() miri.mode = 'IFU' psf = miri.calc_psf(monochromatic=2.8e-6, fov_pixels=10) diff --git a/webbpsf/tests/test_nirspec.py b/webbpsf/tests/test_nirspec.py index 34790f0f..ef15d4d4 100644 --- a/webbpsf/tests/test_nirspec.py +++ b/webbpsf/tests/test_nirspec.py @@ -96,7 +96,7 @@ def test_IFU_wavelengths(): def test_nrs_ifu_broadening(): """ Basic functional test for the code that adjusts PSF outputs to better match empirical IFU PSFs """ - nrs = webbpsf.NIRSpec() + nrs = webbpsf_core.NIRSpec() nrs.mode = 'IFU' psf = nrs.calc_psf(monochromatic=2.8e-6, fov_pixels=10) From be509065acade35ac866bd46e08b72a40313952a Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Fri, 26 Jul 2024 17:56:01 -0400 Subject: [PATCH 6/7] fix tests --- webbpsf/detectors.py | 70 ++++++++++++++++++++++++++++++++++- webbpsf/tests/test_nirspec.py | 1 + webbpsf/webbpsf_core.py | 13 ++++--- 3 files changed, 76 insertions(+), 8 deletions(-) diff --git a/webbpsf/detectors.py b/webbpsf/detectors.py index 6d698ea8..c9b46713 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 @@ -553,13 +555,22 @@ def _show_miri_cruciform_kernel(filt, npix=101, oversample=4, detector_position= # it's implemented as a post-processing modification on the output PSF array. -def apply_miri_ifu_broadening(hdulist, options): +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', 'gaussian') + model_type = options.get('ifu_broadening', 'empirical') if model_type is None or model_type.lower() == 'none': return hdulist @@ -567,14 +578,24 @@ def apply_miri_ifu_broadening(hdulist, options): 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'] = (model_type, "[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 @@ -597,6 +618,7 @@ def apply_nirspec_ifu_broadening(hdulist, options): 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'] @@ -607,3 +629,47 @@ def apply_nirspec_ifu_broadening(hdulist, options): 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/tests/test_nirspec.py b/webbpsf/tests/test_nirspec.py index ef15d4d4..20f28130 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 diff --git a/webbpsf/webbpsf_core.py b/webbpsf/webbpsf_core.py index 3df10bc4..dac906a0 100644 --- a/webbpsf/webbpsf_core.py +++ b/webbpsf/webbpsf_core.py @@ -1278,7 +1278,7 @@ 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 = detectors.apply_miri_ifu_broadening(result, options) + 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) @@ -2071,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). @@ -2460,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] From e686a22ae2b834e62cc56cb2994f19a96d236afe Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 11 Sep 2024 17:05:07 -0400 Subject: [PATCH 7/7] fixes from review --- webbpsf/detectors.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/webbpsf/detectors.py b/webbpsf/detectors.py index c9b46713..63b86fb6 100644 --- a/webbpsf/detectors.py +++ b/webbpsf/detectors.py @@ -550,7 +550,7 @@ def _show_miri_cruciform_kernel(filt, npix=101, oversample=4, detector_position= # Functions for applying IFU optics systematics models # -# Note, thes are is not actually a "Detector" effect, but this file is a +# 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. @@ -586,7 +586,7 @@ def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196): 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'] = (model_type, "[arcsec] IFU PSF broadening Gaussian 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. @@ -623,7 +623,7 @@ def apply_nirspec_ifu_broadening(hdulist, options): 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'] = (model_type, "[arcsec] IFU PSF broadening Gaussian 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']) hdulist[ext].data = out