diff --git a/stpsf/constants.py b/stpsf/constants.py index 1ef44ea..4eb448c 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}, + 'MIRI': {'sigma': 0.05, + 'fwhm_cruciform': 15, + "offset_cruciform": -3}, } diff --git a/stpsf/detectors.py b/stpsf/detectors.py index d05d97f..841c0b2 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) @@ -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)) @@ -548,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 @@ -558,6 +563,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 @@ -570,10 +578,12 @@ 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': + 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 @@ -596,6 +606,25 @@ 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() == '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']) + + 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"]["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) + # enforce strict flux conservation + out *= hdulist[ext].data.sum() / out.sum() hdulist[ext].data = out @@ -631,10 +660,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 ---------- @@ -650,26 +681,94 @@ 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): + """ 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 + else: + return i + + +def _Lorentz1DKernel(amp, fwhm, x_0): + """ 1D Lorentz model as an astropy 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 for the cruciform. + This implements the additional effect of the cruciform on MRS data. + See also _miri_mrs_empirical_broadening + + 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) + 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 flux_conv_factor * (psf_model + amp * psf_model_cruciform) diff --git a/stpsf/match_data.py b/stpsf/match_data.py index 69ec4e5..3af53e2 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,24 @@ 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 +241,38 @@ 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 diff --git a/stpsf/opds.py b/stpsf/opds.py index 3a952eb..e2d2f6b 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} ' diff --git a/stpsf/stpsf_core.py b/stpsf/stpsf_core.py index ad6a1e7..dc4e4e2 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: @@ -1278,33 +1286,36 @@ 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. 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 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_detectors.py b/stpsf/tests/test_detectors.py index 08d6138..8766570 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' diff --git a/stpsf/tests/test_miri.py b/stpsf/tests/test_miri.py index 0874c14..5ee99f1 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,11 +243,21 @@ 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" + # 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()), "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" + + diff --git a/stpsf/utils.py b/stpsf/utils.py index 1f83c25..97b4a95 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)