From 9a11864f523cfbfbeff96bedd37ace4995dfb52a Mon Sep 17 00:00:00 2001 From: Polychronis Patapis Date: Sun, 3 Nov 2024 19:30:58 +0100 Subject: [PATCH 01/12] Implemeted MRI IFU empirical cruciform --- stpsf/constants.py | 2 +- stpsf/detectors.py | 66 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 67 insertions(+), 1 deletion(-) diff --git a/stpsf/constants.py b/stpsf/constants.py index 1ef44ea4..46249c74 100644 --- a/stpsf/constants.py +++ b/stpsf/constants.py @@ -405,5 +405,5 @@ # Parameters for adjusting models of IFU PSFs relative to regular imaging PSFs INSTRUMENT_IFU_BROADENING_PARAMETERS = { 'NIRSPEC': {'sigma': 0.05}, - 'MIRI': {'sigma': 0.05}, + 'MIRI': {'sigma': 0.05, 'fhwm_cruciform': 15, "offset_cruciform": -3}, } diff --git a/stpsf/detectors.py b/stpsf/detectors.py index d05d97f3..e2dc307b 100644 --- a/stpsf/detectors.py +++ b/stpsf/detectors.py @@ -596,6 +596,21 @@ def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196): 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) + elif model_type.lower() == 'cruciform': + # 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) + if wavelen * 1e6 <= 7.5: + amplitude_cruciform = get_mrs_cruciform_amplitude(wavelen * 1e6) + oversample_factor = hdulist[ext].header['DET_SAMP']/7 # optimised parameters with oversampling = 7 + fwhm_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["fhwm_cruciform"]*oversample_factor + offset_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["offset_cruciform"]*oversample_factor + out = _miri_mrs_empirical_cruciform(psf_model=out, amp=amplitude_cruciform, + fwhm=fwhm_cruciform, x_0=offset_cruciform) hdulist[ext].data = out @@ -673,3 +688,54 @@ def _miri_mrs_empirical_broadening(psf_model, alpha_width, beta_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 + + +def get_mrs_cruciform_amplitude(wavelen): + """ + Empirical amplitude of additional cruciform component in MIRI IFU data + wavelen: wavelength (only applicable if wavelength < 7.5 um - see apply_miri_ifu_broadening) + """ + return -0.16765378*wavelen + 1.23632423 # Patapis 2025 PSF paper + + +def _round_up_to_odd_integer(value): + i = np.ceil(value) + if i % 2 == 0: + return i + 1 + else: + return i + + +def _Lorentz1DKernel(amp, fwhm, x_0): + """ + 1D Lorentz model as a convolution kernel + x_size : size of kernel, should be odd + """ + if amp is None: + amp = 2 / (np.pi * fwhm) + + fwhm_int = _round_up_to_odd_integer(fwhm) + return astropy.convolution.Model1DKernel(astropy.modeling.models.Lorentz1D(amplitude=amp, fwhm=fwhm, x_0=x_0), + x_size=int(8 * fwhm_int + 1)) + + +def _miri_mrs_empirical_cruciform(psf_model, amp, fwhm, x_0): + """ + Perform the broadening of a psf model in alpha and beta + + Parameters + ----------- + psf_model : ndarray + webbpsf output results, either monochromatic model or datacube + amp : float + amplitude of Lorentzian in pixels + fwhm : float + Full Width Half Max of Lorentzian in pixels + x_0 : float + offset of Lorentzian in pixels + """ + kernel_cruciform = _Lorentz1DKernel(1.0, fwhm, x_0) + + # TODO: extend algorithm to handle the datacube case + psf_model_cruciform = np.apply_along_axis(lambda m: convolve(m, kernel_cruciform), axis=1, arr=psf_model) + return psf_model+amp*psf_model_cruciform From f98158fde853aa6244f6234a824109daa2fa1e88 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Thu, 7 Nov 2024 14:48:52 -0500 Subject: [PATCH 02/12] fix typo --- stpsf/constants.py | 4 +++- stpsf/detectors.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/stpsf/constants.py b/stpsf/constants.py index 46249c74..4eb448c4 100644 --- a/stpsf/constants.py +++ b/stpsf/constants.py @@ -405,5 +405,7 @@ # Parameters for adjusting models of IFU PSFs relative to regular imaging PSFs INSTRUMENT_IFU_BROADENING_PARAMETERS = { 'NIRSPEC': {'sigma': 0.05}, - 'MIRI': {'sigma': 0.05, 'fhwm_cruciform': 15, "offset_cruciform": -3}, + 'MIRI': {'sigma': 0.05, + 'fwhm_cruciform': 15, + "offset_cruciform": -3}, } diff --git a/stpsf/detectors.py b/stpsf/detectors.py index e2dc307b..beefc0b5 100644 --- a/stpsf/detectors.py +++ b/stpsf/detectors.py @@ -607,7 +607,7 @@ def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196): if wavelen * 1e6 <= 7.5: amplitude_cruciform = get_mrs_cruciform_amplitude(wavelen * 1e6) oversample_factor = hdulist[ext].header['DET_SAMP']/7 # optimised parameters with oversampling = 7 - fwhm_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["fhwm_cruciform"]*oversample_factor + fwhm_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["fwhm_cruciform"]*oversample_factor offset_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["offset_cruciform"]*oversample_factor out = _miri_mrs_empirical_cruciform(psf_model=out, amp=amplitude_cruciform, fwhm=fwhm_cruciform, x_0=offset_cruciform) From 598c005c8d600f1b7412f990b3eb1840a7367720 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Thu, 7 Nov 2024 14:49:55 -0500 Subject: [PATCH 03/12] don't warn a field point warning for MRS field point --- stpsf/opds.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/stpsf/opds.py b/stpsf/opds.py index 3a952ebe..6ab2798a 100644 --- a/stpsf/opds.py +++ b/stpsf/opds.py @@ -1748,7 +1748,14 @@ def _get_zernikes_for_ote_field_dep( y_field_pt = np.clip(y_field_pt0, min_y_field, max_y_field) clip_dist = np.sqrt((x_field_pt - x_field_pt0) ** 2 + (y_field_pt - y_field_pt0) ** 2) - if clip_dist > 0.1 * u.arcsec: + + # special case for MIRI: don't do the following warning specifically for the MRS field point since it's + # a known issue that # point is slightly outside of the valid region. This specific case is benign and + # it's not helpful to emit this warning for every MRS calculation ever + mrs_v2v3 = [-8.39108833, -5.32144667]*u.arcmin + is_mrs_fieldpoint = np.allclose(v2v3, mrs_v2v3, atol=0.01) + + if (clip_dist > 0.1 * u.arcsec) and not is_mrs_fieldpoint: # warn the user we're making an adjustment here (but no need to do so if the distance is trivially small) warning_message = ( f'For (V2,V3) = {v2v3}, Field point {x_field_pt}, {y_field_pt} ' From b643763c96e53acc35cfc6306e7f84d8eec891d6 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Fri, 8 Nov 2024 17:21:07 -0500 Subject: [PATCH 04/12] add warnings for MRS cubes other than ifualign and single-band --- stpsf/match_data.py | 51 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 50 insertions(+), 1 deletion(-) diff --git a/stpsf/match_data.py b/stpsf/match_data.py index 69ec4e54..335f6551 100644 --- a/stpsf/match_data.py +++ b/stpsf/match_data.py @@ -1,4 +1,5 @@ # Functions to match or fit PSFs to observed JWST data +import warnings import astropy import astropy.io.fits as fits import pysiaf @@ -74,8 +75,23 @@ def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choic elif inst.name == 'MIRI': if header['EXP_TYPE'] == 'MIR_MRS': ch = header['CHANNEL'] + band = header['BAND'] + inferred_output_type, inferred_coord_system = infer_mrs_cube_type(filename_or_HDUList) + if band == 'MULTIPLE' or inferred_output_type != 'band': + warnings.warn(f"** The input file seems to be an MRS datacube with output_type='{inferred_output_type}', " + "combining multiple bands. Note that PSF models can be computed for only 1 band at a time. " + "You will need to make multiple PSF simulations to model the PSF in this dataset. For high " + "precision work, be aware of the small (<1 deg) rotation differences between the individual bands. **") + band = 'SHORT' # just pick one, arbitrarily + if inferred_coord_system != 'ifualign': + warnings.warn(f"** The input file seems to be an MRS datacube with coord_system='skyalign'. Note that PSF " + "models can be computed only for coord_system='ifualign'.. You will need to either re-reduce " + "your data using the ifualign coord_system (preferred) or rotate the PSF model based on the " + "position angle (less preferred, due to numerical interpolation noise). **") band_lookup = {'SHORT': 'A', 'MEDIUM': 'B', 'LONG': 'C'} - inst.band = str(ch) + band_lookup[header['BAND']] + inst.band = str(ch) + band_lookup[band] + + elif inst.filter in ['F1065C', 'F1140C', 'F1550C']: inst.image_mask = 'FQPM' + inst.filter[1:5] @@ -224,3 +240,36 @@ def get_nrc_coron_mask_from_pps_apername(apname_pps): image_mask = image_mask[:-1] return image_mask + +def infer_mrs_cube_type(filename, verbose=False): + """attempt to infer cube coordinate system and output type from header metadata. + The cube_build step doesn't record in metadata several of its key input parameters; this function attempts to infer what their values were. + + Returns (guessed_output_type, guessed_coord_system) + """ + with fits.open(filename) as hdul: + + if hdul[0].header['INSTRUME'] != 'MIRI': + raise RuntimeError("This function only applies to MIRI MRS data") + + naxis3 = hdul['SCI'].header['NAXIS3'] + # Infer cube output type, based on the number of wavelengths in this cube. + # This is inelegant but seems to work sufficiently + if naxis3 < 2000: + output_type = 'band' + elif naxis3 > 2000 and naxis3 < 4000: + output_type = 'channel' + elif naxis3 > 4000: + output_type = 'multi' + + # Infer coordinate system, based on PC rotation matrix (?) + pc11 = hdul['SCI'].header['PC1_1'] + pc22 = hdul['SCI'].header['PC2_2'] + if pc11 == -1 and pc22 == 1: + coord_system = 'skyalign' + else: + coord_system = 'ifualign' + + if verbose: + print(filename, naxis3, output_type, pc11, pc22, coord_system) + return output_type, coord_system From 0c53fab9ded5a2dd02d9177fd24a7766a5e2c01a Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Fri, 8 Nov 2024 17:32:01 -0500 Subject: [PATCH 05/12] pep8 lint --- stpsf/detectors.py | 10 +++++----- stpsf/match_data.py | 19 +++++++++++-------- stpsf/opds.py | 2 +- 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/stpsf/detectors.py b/stpsf/detectors.py index beefc0b5..02591dac 100644 --- a/stpsf/detectors.py +++ b/stpsf/detectors.py @@ -606,9 +606,9 @@ def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196): out = _miri_mrs_empirical_broadening(psf_model=hdulist[ext].data, alpha_width=alpha_width, beta_width=beta_width) if wavelen * 1e6 <= 7.5: amplitude_cruciform = get_mrs_cruciform_amplitude(wavelen * 1e6) - oversample_factor = hdulist[ext].header['DET_SAMP']/7 # optimised parameters with oversampling = 7 - fwhm_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["fwhm_cruciform"]*oversample_factor - offset_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["offset_cruciform"]*oversample_factor + oversample_factor = hdulist[ext].header['DET_SAMP'] / 7 # optimised parameters with oversampling = 7 + fwhm_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["fwhm_cruciform"] * oversample_factor + offset_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["offset_cruciform"] * oversample_factor out = _miri_mrs_empirical_cruciform(psf_model=out, amp=amplitude_cruciform, fwhm=fwhm_cruciform, x_0=offset_cruciform) @@ -695,7 +695,7 @@ def get_mrs_cruciform_amplitude(wavelen): Empirical amplitude of additional cruciform component in MIRI IFU data wavelen: wavelength (only applicable if wavelength < 7.5 um - see apply_miri_ifu_broadening) """ - return -0.16765378*wavelen + 1.23632423 # Patapis 2025 PSF paper + return -0.16765378 * wavelen + 1.23632423 # Patapis 2025 PSF paper def _round_up_to_odd_integer(value): @@ -738,4 +738,4 @@ def _miri_mrs_empirical_cruciform(psf_model, amp, fwhm, x_0): # TODO: extend algorithm to handle the datacube case psf_model_cruciform = np.apply_along_axis(lambda m: convolve(m, kernel_cruciform), axis=1, arr=psf_model) - return psf_model+amp*psf_model_cruciform + return psf_model + amp * psf_model_cruciform diff --git a/stpsf/match_data.py b/stpsf/match_data.py index 335f6551..3af53e25 100644 --- a/stpsf/match_data.py +++ b/stpsf/match_data.py @@ -79,15 +79,16 @@ def setup_sim_to_match_file(filename_or_HDUList, verbose=True, plot=False, choic inferred_output_type, inferred_coord_system = infer_mrs_cube_type(filename_or_HDUList) if band == 'MULTIPLE' or inferred_output_type != 'band': warnings.warn(f"** The input file seems to be an MRS datacube with output_type='{inferred_output_type}', " - "combining multiple bands. Note that PSF models can be computed for only 1 band at a time. " - "You will need to make multiple PSF simulations to model the PSF in this dataset. For high " - "precision work, be aware of the small (<1 deg) rotation differences between the individual bands. **") - band = 'SHORT' # just pick one, arbitrarily + "combining multiple bands. Note that PSF models can be computed for only 1 band at a time. " + "You will need to make multiple PSF simulations to model the PSF in this dataset. For high " + "precision work, be aware of the small (<1 deg) rotation differences between the individual" + "bands. **") + band = 'SHORT' # just pick one, arbitrarily if inferred_coord_system != 'ifualign': warnings.warn(f"** The input file seems to be an MRS datacube with coord_system='skyalign'. Note that PSF " - "models can be computed only for coord_system='ifualign'.. You will need to either re-reduce " - "your data using the ifualign coord_system (preferred) or rotate the PSF model based on the " - "position angle (less preferred, due to numerical interpolation noise). **") + "models can be computed only for coord_system='ifualign'.. You will need to either re-reduce " + "your data using the ifualign coord_system (preferred) or rotate the PSF model based on the " + "position angle (less preferred, due to numerical interpolation noise). **") band_lookup = {'SHORT': 'A', 'MEDIUM': 'B', 'LONG': 'C'} inst.band = str(ch) + band_lookup[band] @@ -241,9 +242,11 @@ def get_nrc_coron_mask_from_pps_apername(apname_pps): return image_mask + def infer_mrs_cube_type(filename, verbose=False): """attempt to infer cube coordinate system and output type from header metadata. - The cube_build step doesn't record in metadata several of its key input parameters; this function attempts to infer what their values were. + The cube_build step doesn't record in metadata several of its key input parameters; + this function attempts to infer what their values were. Returns (guessed_output_type, guessed_coord_system) """ diff --git a/stpsf/opds.py b/stpsf/opds.py index 6ab2798a..e2d2f6bb 100644 --- a/stpsf/opds.py +++ b/stpsf/opds.py @@ -1752,7 +1752,7 @@ def _get_zernikes_for_ote_field_dep( # special case for MIRI: don't do the following warning specifically for the MRS field point since it's # a known issue that # point is slightly outside of the valid region. This specific case is benign and # it's not helpful to emit this warning for every MRS calculation ever - mrs_v2v3 = [-8.39108833, -5.32144667]*u.arcmin + mrs_v2v3 = [-8.39108833, -5.32144667] * u.arcmin is_mrs_fieldpoint = np.allclose(v2v3, mrs_v2v3, atol=0.01) if (clip_dist > 0.1 * u.arcsec) and not is_mrs_fieldpoint: From c3809bedd398427d1d6cc7db2e93698010d90fc1 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Sat, 30 Nov 2024 08:48:11 -0500 Subject: [PATCH 06/12] rename functions for MIRI cruciform, for clarity --- stpsf/detectors.py | 29 ++++++++++++++++++----------- stpsf/stpsf_core.py | 2 +- stpsf/tests/test_detectors.py | 6 +++--- 3 files changed, 22 insertions(+), 15 deletions(-) diff --git a/stpsf/detectors.py b/stpsf/detectors.py index 02591dac..7d6d30c5 100644 --- a/stpsf/detectors.py +++ b/stpsf/detectors.py @@ -297,7 +297,7 @@ def oversample_ipc_model(kernel, oversample): cruciform_yshifts = scipy.interpolate.interp1d([0, 511, 1031], [1.6, 0, -1.6], kind='linear', fill_value='extrapolate') -def _make_miri_scattering_kernel_2d(in_psf, kernel_amp, oversample=1, wavelength=5.5, detector_position=(0, 0)): +def _make_miri_cruciform_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. @@ -362,9 +362,9 @@ def _make_miri_scattering_kernel_2d(in_psf, kernel_amp, oversample=1, wavelength return kernel_2d -def _apply_miri_scattering_kernel_2d(in_psf, kernel_2d, oversample): +def _apply_miri_cruciform_kernel_2d(in_psf, kernel_2d, oversample): """ - Applies the detector scattering kernel created in _make_miri_scattering_kernel + Applies the detector scattering kernel created in _make_miri_cruciform_kernel_2d function to an input image. Code is adapted from MIRI-TN-00076-ATC_Imager_PSF_Issue_4.pdf @@ -372,9 +372,8 @@ def _apply_miri_scattering_kernel_2d(in_psf, kernel_2d, oversample): ---------- 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. + kernel_2d : ndarray + The 2D kernel, output from _make_miri_cruciform_kernel_2d. oversample : int Amount by which the input PSF is oversampled @@ -443,7 +442,7 @@ def get_miri_cruciform_amplitude(filt): return kernel_amp -def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None, old_method=False): +def apply_miri_imager_cruciform(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 @@ -456,6 +455,9 @@ def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None, old_method= which this will then modify. This happens in stpsf_core.py/JWInstrument._calc_psf_format_output, which is where this is called from in the usual course of operation. + Note - this function is called only for **MIRI IMAGER** simulations. + For MRS, see instead the apply_miri_ifu_broadening function, below. + Parameters ---------- hdulist_or_filename : @@ -502,14 +504,14 @@ def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None, old_method= in_psf = psf[ext].data # create cruciform model using improved method using a 2d convolution kernel, attempting to model more physics. - kernel_2d = _make_miri_scattering_kernel_2d( + kernel_2d = _make_miri_cruciform_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) + im_conv_both = _apply_miri_cruciform_kernel_2d(in_psf, kernel_2d, oversample) # Add this 2D scattered light output to the PSF psf_new = in_psf + im_conv_both @@ -529,14 +531,16 @@ def apply_miri_scattering(hdulist_or_filename=None, kernel_amp=None, old_method= def _show_miri_cruciform_kernel(filt, npix=101, oversample=4, detector_position=(512, 512), ax=None): - """utility function for viewing/visualizing the cruciform kernel""" + """utility function for viewing/visualizing the cruciform kernel. + This displays a plot on screen, as a user convenience. + This display function is Not used as part of PSF calculations.""" 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) + kernel_2d = _make_miri_cruciform_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)) @@ -558,6 +562,9 @@ def _show_miri_cruciform_kernel(filt, npix=101, oversample=4, detector_position= 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 + Note - this function is called only for **MIRI MRS** simulations. + For MRS, see instead the apply_miri_imager_cruciform function, above. + Parameters ----------- hdulist : astropy.io.fits.HDUList diff --git a/stpsf/stpsf_core.py b/stpsf/stpsf_core.py index ad6a1e7f..7a7acafd 100644 --- a/stpsf/stpsf_core.py +++ b/stpsf/stpsf_core.py @@ -1286,7 +1286,7 @@ def _calc_psf_format_output(self, result, options): # slit type aperture, specifically LRS SLIT, does not have distortion polynomials # therefore omit apply_distortion if a SLIT aperture is selected. psf_siaf = result - psf_siaf_rot = detectors.apply_miri_scattering(psf_siaf) # apply scattering effect + psf_siaf_rot = detectors.apply_miri_imager_cruciform(psf_siaf) # apply scattering effect psf_distorted = detectors.apply_detector_charge_diffusion( psf_siaf_rot, options ) # apply detector charge transfer model diff --git a/stpsf/tests/test_detectors.py b/stpsf/tests/test_detectors.py index 08d6138c..87665704 100644 --- a/stpsf/tests/test_detectors.py +++ b/stpsf/tests/test_detectors.py @@ -17,7 +17,7 @@ def test_apply_miri_scattering_error(): # Test that running this function will raise a ValueError with pytest.raises(ValueError) as excinfo: - detectors.apply_miri_scattering(psf) + detectors.apply_miri_imager_cruciform(psf) assert 'ValueError' in str(excinfo), 'Non-MIRI PSFs should not be able to run through apply_miri_scattering' @@ -45,7 +45,7 @@ def test_apply_miri_scattering(): psf[ext_new].header['EXTNAME'] = psf[ext].header['EXTNAME'][0:4] + 'DIST' # change extension name # Run it through just the apply_miri_scattering function - psf_cross = detectors.apply_miri_scattering(psf) + psf_cross = detectors.apply_miri_imager_cruciform(psf) # Rebin data to get 3rd extension mir.options['output_mode'] = 'Both extensions' @@ -137,7 +137,7 @@ def test_miri_conservation_energy(): psf[ext_new].header['EXTNAME'] = psf[ext].header['EXTNAME'][0:4] + 'DIST' # change extension name # Run it through just the apply_miri_scattering function - psf_cross = detectors.apply_miri_scattering(psf) + psf_cross = detectors.apply_miri_imager_cruciform(psf) # Rebin data to get 3rd extension mir.options['output_mode'] = 'Both extensions' From b3d47d08c1608aeef03bc773c33b51c5bc2f7c11 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Sat, 30 Nov 2024 08:48:44 -0500 Subject: [PATCH 07/12] Add util function to display all extensions in a PSF calc. Should've done this long ago! --- stpsf/utils.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/stpsf/utils.py b/stpsf/utils.py index 1f83c25a..97b4a953 100644 --- a/stpsf/utils.py +++ b/stpsf/utils.py @@ -1069,4 +1069,18 @@ def get_target_phase_map_filename(apername): raise ValueError("File wss_target_phase_{}.fits, \ not found under {}.".format(apername.split('_')[1].lower(),path)) - return was_targ_file \ No newline at end of file + return was_targ_file + + +def display_psf_exts(psf, figsize=(12,3), imagecrop=5, colorbar=False, **kwargs): + """ Display all extensions of a given PSF. + + See display_psf for additional information on parameters. + """ + import poppy + n_exts = len(psf) + + fig, axes = plt.subplots(figsize=figsize, ncols=4) + for ext in range(n_exts): + poppy.display_psf(psf, ext=ext, ax=axes[ext], title=f'Ext {ext}: {psf[ext].header["EXTNAME"]}', + imagecrop=imagecrop, colorbar=colorbar, **kwargs) From e3239c1133fe0e784e9c5b06405ebe405a0ffbb0 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Sat, 30 Nov 2024 09:26:11 -0500 Subject: [PATCH 08/12] improve doc strings; fix mis-indented code --- stpsf/detectors.py | 70 ++++++++++++++++++++++++++++------------------ 1 file changed, 43 insertions(+), 27 deletions(-) diff --git a/stpsf/detectors.py b/stpsf/detectors.py index 7d6d30c5..2e307fd4 100644 --- a/stpsf/detectors.py +++ b/stpsf/detectors.py @@ -289,7 +289,7 @@ def oversample_ipc_model(kernel, oversample): return kernel_oversample -# Functions for applying MIRI Detector Scattering Effect +# Functions for applying MIRI Cruciform Detector "Scattering" Diffraction Effect # Lookup tables of shifts of the cruciform # Estimated roughly from F560W ePSFs (ePSFs by Libralatto, shift estimate by Perrin) @@ -552,6 +552,7 @@ def _show_miri_cruciform_kernel(filt, npix=101, oversample=4, detector_position= 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 @@ -611,6 +612,7 @@ def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196): 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) + if wavelen * 1e6 <= 7.5: amplitude_cruciform = get_mrs_cruciform_amplitude(wavelen * 1e6) oversample_factor = hdulist[ext].header['DET_SAMP'] / 7 # optimised parameters with oversampling = 7 @@ -653,10 +655,12 @@ def apply_nirspec_ifu_broadening(hdulist, options): return hdulist +# Functions for modeling MIRI MRS effects, including both the cruciform (as seen in the MRS) and other IFU broadening. +# Invoked from apply_miri_ifu_broadening, above. + 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. + """ Calculate the Gaussian scale of the kernel that broadens the MRS IFU PSF, + from the diffraction limited FWHM to the empirically measured FWHM. Parameters ---------- @@ -672,40 +676,51 @@ def _miri_mrs_analytical_sigma_alpha_broadening(wavelength): def _miri_mrs_empirical_broadening(psf_model, alpha_width, beta_width): - """ - Perform the broadening of a psf model in alpha and beta + """ Perform the broadening of a MRS PSF model in alpha and beta. + This only includes the IFU broadening. See also _miri_mrs_empirical_cruciform + + The beta (across-slice) convolution is implemented as a Box1D kernel + The alpha (along-slice) convolution is implemented as a Gaussian kernel - Parameters - ----------- - psf_model : ndarray + Parameters + ----------- + psf_model : ndarray stpsf output results, eitehr monochromatic model or datacube - alpha_width : float + alpha_width : float gaussian convolution kernel in pixels, None if no broadening should be performed - beta_width : float + beta_width : float slice width in pixels - """ - kernel_beta = astropy.convolution.Box1DKernel(beta_width) + """ + kernel_beta = astropy.convolution.Box1DKernel(beta_width) + webbpsf.webbpsf_core._log.info(f' MRS empirical broadening: alpha width {alpha_width}, beta width {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 + 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 def get_mrs_cruciform_amplitude(wavelen): - """ - Empirical amplitude of additional cruciform component in MIRI IFU data - wavelen: wavelength (only applicable if wavelength < 7.5 um - see apply_miri_ifu_broadening) + """ Calculate empirical amplitude of additional cruciform component in MIRI IFU data + See Patapis et al. 2025 + + Parameters + ---------- + wavelen : float + wavelength (only applicable if wavelength < 7.5 um - see apply_miri_ifu_broadening) """ return -0.16765378 * wavelen + 1.23632423 # Patapis 2025 PSF paper def _round_up_to_odd_integer(value): + """ ensure an integer is odd. + Minor utility function used by convolution kernel creation""" i = np.ceil(value) if i % 2 == 0: return i + 1 @@ -714,8 +729,8 @@ def _round_up_to_odd_integer(value): def _Lorentz1DKernel(amp, fwhm, x_0): - """ - 1D Lorentz model as a convolution kernel + """ 1D Lorentz model as an astropy convolution kernel + x_size : size of kernel, should be odd """ if amp is None: @@ -727,8 +742,9 @@ def _Lorentz1DKernel(amp, fwhm, x_0): def _miri_mrs_empirical_cruciform(psf_model, amp, fwhm, x_0): - """ - Perform the broadening of a psf model in alpha and beta + """ Perform the broadening of a psf model in alpha and beta for the cruciform. + This implements the additional effect of the cruciform on MRS data. + See also _miri_mrs_empirical_broadening Parameters ----------- From 4cc7d40d1a50ff7b890d203cdce9d623a5eee31e Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Sat, 30 Nov 2024 09:26:37 -0500 Subject: [PATCH 09/12] Conserve flux in MIRI MRS IFU cruciform --- stpsf/detectors.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/stpsf/detectors.py b/stpsf/detectors.py index 2e307fd4..f9b64de3 100644 --- a/stpsf/detectors.py +++ b/stpsf/detectors.py @@ -754,11 +754,16 @@ def _miri_mrs_empirical_cruciform(psf_model, amp, fwhm, x_0): amplitude of Lorentzian in pixels fwhm : float Full Width Half Max of Lorentzian in pixels - x_0 : float + x_0 : float offset of Lorentzian in pixels """ kernel_cruciform = _Lorentz1DKernel(1.0, fwhm, x_0) + webbpsf.webbpsf_core._log.info(f' MRS empirical cruciform: amp {amp} fwhm {fwhm} x_0 {x_0}') + + # Flux conservation: the integral of the Lorentz kernel is the same as the amplitude + # therefore the following will keep the total flux conserved in the summed output PSF. + flux_conv_factor = 1 / (1 + amp) # TODO: extend algorithm to handle the datacube case psf_model_cruciform = np.apply_along_axis(lambda m: convolve(m, kernel_cruciform), axis=1, arr=psf_model) - return psf_model + amp * psf_model_cruciform + return flux_conv_factor * (psf_model + amp * psf_model_cruciform) From 8b9f8d95a121775b4913fdc9910fdbf83062d23b Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Sat, 30 Nov 2024 09:39:59 -0500 Subject: [PATCH 10/12] make empirical_cruciform the default for MRS sims. enhance test case. --- stpsf/detectors.py | 5 +++-- stpsf/tests/test_miri.py | 4 +++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/stpsf/detectors.py b/stpsf/detectors.py index f9b64de3..b61bda5a 100644 --- a/stpsf/detectors.py +++ b/stpsf/detectors.py @@ -578,7 +578,7 @@ def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196): """ # 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') + model_type = options.get('ifu_broadening', 'empirical_cruciform') if model_type is None or model_type.lower() == 'none': return hdulist @@ -604,7 +604,8 @@ def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196): 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) - elif model_type.lower() == 'cruciform': + elif model_type.lower() == 'empirical_cruciform': + # The above empirical mode, plus an additional cruciform effect at < 7.5 microns # Model based on empirical PSF properties, Argryiou et al. pixelscl = float(hdulist[ext].header['PIXELSCL']) wavelen = float(hdulist[ext].header['WAVELEN']) diff --git a/stpsf/tests/test_miri.py b/stpsf/tests/test_miri.py index 0874c141..9213374c 100644 --- a/stpsf/tests/test_miri.py +++ b/stpsf/tests/test_miri.py @@ -233,7 +233,7 @@ def test_miri_ifu_broadening(): miri = stpsf_core.MIRI() miri.mode = 'IFU' - psf = miri.calc_psf(monochromatic=2.8e-6, fov_pixels=10) + psf = miri.calc_psf(monochromatic=6.8e-6, fov_pixels=20) fwhm_oversamp = stpsf.measure_fwhm(psf, ext='OVERSAMP') fwhm_overdist = stpsf.measure_fwhm(psf, ext='OVERDIST') @@ -243,6 +243,8 @@ def test_miri_ifu_broadening(): fwhm_detdist = stpsf.measure_fwhm(psf, ext='DET_DIST') assert fwhm_overdist > fwhm_oversamp, "IFU broadening model should increase the FWHM for the distorted extensions" + assert np.isclose(psf['DET_SAMP'].data.sum(), psf['DET_DIST'].data.sum(), rtol=5e-3), "IFU broadening should not change total flux much" + # 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) From 3408e31da500f148bdb0d80811a263a4339e78f2 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 4 Dec 2024 13:45:06 -0500 Subject: [PATCH 11/12] improve flux conservation and testing for MRS IFU sims --- stpsf/detectors.py | 2 ++ stpsf/stpsf_core.py | 7 +++++-- stpsf/tests/test_miri.py | 10 +++++++++- 3 files changed, 16 insertions(+), 3 deletions(-) diff --git a/stpsf/detectors.py b/stpsf/detectors.py index b61bda5a..f331bc2f 100644 --- a/stpsf/detectors.py +++ b/stpsf/detectors.py @@ -621,6 +621,8 @@ def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196): offset_cruciform = constants.INSTRUMENT_IFU_BROADENING_PARAMETERS["MIRI"]["offset_cruciform"] * oversample_factor out = _miri_mrs_empirical_cruciform(psf_model=out, amp=amplitude_cruciform, fwhm=fwhm_cruciform, x_0=offset_cruciform) + # enforce strict flux conservation + out *= hdulist[ext].data.sum() / out.sum() hdulist[ext].data = out diff --git a/stpsf/stpsf_core.py b/stpsf/stpsf_core.py index 7a7acafd..c1af5b4e 100644 --- a/stpsf/stpsf_core.py +++ b/stpsf/stpsf_core.py @@ -1278,10 +1278,11 @@ def _calc_psf_format_output(self, result, options): ) # apply detector charge transfer model elif self.name == 'MIRI': # Apply distortion effects to MIRI psf: Distortion and MIRI Scattering - _log.debug('MIRI: Adding optical distortion and Si:As detector internal scattering') if self.mode != 'IFU': + _log.debug('MIRI imager: Adding optical distortion and Si:As detector internal scattering') if self._detector_geom_info.aperture.AperType != 'SLIT': psf_siaf = distortion.apply_distortion(result) # apply siaf distortion + _log.debug('MIRI: Applied optical distortion based on SIAF parameters') else: # slit type aperture, specifically LRS SLIT, does not have distortion polynomials # therefore omit apply_distortion if a SLIT aperture is selected. @@ -1291,20 +1292,22 @@ def _calc_psf_format_output(self, result, options): psf_siaf_rot, options ) # apply detector charge transfer model else: + _log.debug('MIRI MRS: Adding IFU PSF broadening effects.') # 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, 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 self.mode != 'IFU': + _log.debug('NIRSpec: Adding optical distortion and detector charge transfer') 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: + _log.debug('NIRSpec IFU: Adding IFU PSF broadening') # there is not yet any distortion calibration for the IFU. psf_distorted = detectors.apply_nirspec_ifu_broadening(result, options) diff --git a/stpsf/tests/test_miri.py b/stpsf/tests/test_miri.py index 9213374c..757f546c 100644 --- a/stpsf/tests/test_miri.py +++ b/stpsf/tests/test_miri.py @@ -243,7 +243,10 @@ def test_miri_ifu_broadening(): fwhm_detdist = stpsf.measure_fwhm(psf, ext='DET_DIST') assert fwhm_overdist > fwhm_oversamp, "IFU broadening model should increase the FWHM for the distorted extensions" - assert np.isclose(psf['DET_SAMP'].data.sum(), psf['DET_DIST'].data.sum(), rtol=5e-3), "IFU broadening should not change total flux much" + # test flux conservation. THis is close but not exact, due to the optical distortion part which is distinct from but + # happens at same point in the calculation as the IFU broadening effect. + assert np.isclose(psf['DET_SAMP'].data.sum(), psf['DET_DIST'].data.sum(), rtol=2e-4), "IFU broadening should not change total flux much" + # Now test that we can also optionally turn off that effect miri.options['ifu_broadening'] = None @@ -253,3 +256,8 @@ def test_miri_ifu_broadening(): fwhm_overdist = stpsf.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" + + # test flux conservation with and without the IFU broadening. This should be a more exact match + assert np.isclose(psf_nb['DET_DIST'].data.sum(), psf['DET_DIST'].data.sum() ), "IFU broadening should not change total flux much" + + From 95dad0b39d02ed0cf1f6112693a46aaf9e80ada1 Mon Sep 17 00:00:00 2001 From: Marshall Perrin Date: Wed, 4 Dec 2024 17:11:59 -0500 Subject: [PATCH 12/12] Avoid adding detector IPC to MIRI MRS IFU sims. Improve testing rigor & consistency --- stpsf/detectors.py | 2 ++ stpsf/stpsf_core.py | 10 +++++++++- stpsf/tests/test_miri.py | 6 +++--- 3 files changed, 14 insertions(+), 4 deletions(-) diff --git a/stpsf/detectors.py b/stpsf/detectors.py index f331bc2f..841c0b28 100644 --- a/stpsf/detectors.py +++ b/stpsf/detectors.py @@ -581,7 +581,9 @@ def apply_miri_ifu_broadening(hdulist, options, slice_width=0.196): model_type = options.get('ifu_broadening', 'empirical_cruciform') if model_type is None or model_type.lower() == 'none': + webbpsf.webbpsf_core._log.debug('MIRI MRS: IFU broadening option is set to None, skipping IFU PSF broadening effects.') return hdulist + webbpsf.webbpsf_core._log.debug('MIRI MRS: Adding IFU PSF broadening effects.') ext = 1 # Apply this effect to the OVERDIST extension, which at this point in the code will be ext 1 diff --git a/stpsf/stpsf_core.py b/stpsf/stpsf_core.py index c1af5b4e..dc4e4e25 100644 --- a/stpsf/stpsf_core.py +++ b/stpsf/stpsf_core.py @@ -1147,6 +1147,12 @@ def _get_pixelscale_from_apername(self, apername): # The slight departures from this are handled in the distortion model; see distortion.py return (ap.XSciScale + ap.YSciScale) / 2 + @property + def mode(self): + # This exists just for API consistency with the subclasses that have an imaging vs IFU mode toggle, + # so that all JWST instrument classes have a mode attribute, consistently, whether used or not + return "imaging" + def _get_fits_header(self, result, options): """populate FITS Header keywords""" super(JWInstrument, self)._get_fits_header(result, options) @@ -1252,7 +1258,9 @@ def _calc_psf_format_output(self, result, options): add_distortion = options.get('add_distortion', True) crop_psf = options.get('crop_psf', True) # you can turn on/off IPC corrections via the add_ipc option, default True. - add_ipc = options.get('add_ipc', True) + # except for IFU mode simulations, it doesn't make sense to add the regular detector IPC + # instead the more complex IFU broadening models should be applied + add_ipc = options.get('add_ipc', True if self.mode != 'IFU' else False) # Add distortion if set in calc_psf if add_distortion: diff --git a/stpsf/tests/test_miri.py b/stpsf/tests/test_miri.py index 757f546c..5ee99f1c 100644 --- a/stpsf/tests/test_miri.py +++ b/stpsf/tests/test_miri.py @@ -245,17 +245,17 @@ def test_miri_ifu_broadening(): # test flux conservation. THis is close but not exact, due to the optical distortion part which is distinct from but # happens at same point in the calculation as the IFU broadening effect. - assert np.isclose(psf['DET_SAMP'].data.sum(), psf['DET_DIST'].data.sum(), rtol=2e-4), "IFU broadening should not change total flux much" + assert np.isclose(psf['DET_SAMP'].data.sum(), psf['DET_DIST'].data.sum()), "IFU broadening should not change total flux much" # 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) + psf_nb = miri.calc_psf(monochromatic=6.8e-6, fov_pixels=20) fwhm_oversamp = stpsf.measure_fwhm(psf_nb, ext='OVERSAMP') fwhm_overdist = stpsf.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" + assert fwhm_overdist == fwhm_oversamp, "IFU broadening model should be disabled for this test case" # test flux conservation with and without the IFU broadening. This should be a more exact match assert np.isclose(psf_nb['DET_DIST'].data.sum(), psf['DET_DIST'].data.sum() ), "IFU broadening should not change total flux much"