Simulating IFU mode and Datacubes

This page describes how to simulate IFU mode data and spectral datacubes for JWST’s NIRSpec IFU and MIRI MRS.

IFU support is a relatively recent addition, and efforts are still in progress to refine and improve the fidelity of IFU mode simulations.

Selecting IFU mode simulations

These instruments have a mode attribute, which can be set to either ‘IFU’ or ‘imaging’. Selecting IFU mode has the following effects:

• The normal filter attribute for selecting spectral bandpass is superceded by a band attribute for selected IFU bands, the specific details of which differ for NIRSpec and MIRI. A get_IFU_wavelengths() function is added, which allows looking up the wavelength range for each band.

• The PSF simulation gets an added step for adding “IFU broadening” effects, which are empirical models for slightly broadening/blurring the simulated PSF to better match the observed PSF FWHMs. Physically this is a simplified model for optical blurring effects due to imperfect wavefront quality in the IFU image slicer optics, for example.

• For NIRSpec IFU simulations only, the PSF output is rotated by an additional 90 degrees to match the orientation of JWST pipeline output dataproducts created with the “orient=’IFUalign’” option in the Cube Build step.

Note there are three options for computing PSFs in IFU mode.

1. If you only need one wavelength, use regular calc_psf() with the monochromatic keyword to specify a wavelength.

2. If you want a datacube at all wavelengths, use calc_datacube() with a list or array of wavelengths. This is the recommended path for many typical use cases.

3. If you want a datacube at all wavelengths, you can also use calc_datacube_fast() which achieves a much-faster calculation runtime by making a simplifying assumption in the optical calculation, and currently by not including the detector effects or distortion modeling steps.

   • Specifically, it assumes that the wavefront optical path difference in the IFU exit pupil is independent of wavelength. This assumption is reasonably true for both MIRI and NIRSpec IFU modes within the current level of fidelity.

Datacube Coordinate Frames for IFU PSFs: SkyAlign vs IFUAlign

STPSF is intended to simulate IFU PSFs in IFUAlign orientation only.

For IFU observations, the JWST pipeline cube_build step (see pipeline docs) can produce output datacubes in two different reference frames or orientations:

• coord_system='skyalign'. In these cubes, north is up and east is left. This is the default orientation.

• coord_system='ifualign'. These cubes are build on the local IFU slicer plane’s native coordinate system. The X and Y coordinates corresponds to the along-slice and across-slice directions, respectively (which are often labeled as the \(\alpha\) and \(\beta\) axes for MIRI MRS; see JDox.)

PSF simulations are always produced in the ifualign orientation only. For PSF modeling work, we recommend reducing your science data using the ifualign frame. This allows PSF models and orientations to be consistent between multiple datasets, independent of position angles on the sky.

If your science data has been reduced in skyalign frame, you can rotate the PSF based on the aperture position angle to match your data. However this will introduce numerical interpolation noise, particularly for spatially undersampled IFU data. That can be avoided by working in IFUalign frame, as we recommend here.

Pixelscales: By default, IFU PSFs are made using the same default pixel scales as the pipeline defaults (0.1 arcsec/pixel for NIRSpec, and variously 0.13 - 0.35 arcsec/pixel for MIRI depending on which MRS channel). If you have reduced your science data using a different pixelscale in the JWST pipeline, e.g. 0.123 arcsec/pix, simply set e.g. miri.pixelscale = 0.123 to the same value for creating the PSF model.

NIRSpec IFU example

For NIRSpec, the band attribute is derived from the prism and disperser elements.

nrs = stpsf.NIRSpec()
nrs.mode = 'IFU'

nrs.disperser = 'PRISM'
nrs.filter = 'CLEAR'
print("Band is", nrs.band)

Band is PRISM/CLEAR

The wavelength sampling can be obtained using the get_IFU_wavelengths() function. By default this returns the same wavelength sampling as the pipeline uses. But if desired you can also specify some other number of wavelengths nlambda, for instance to reduce simulation runtimes when the full spectral resolution is not needed.

allwaves = nrs.get_IFU_wavelengths()
print(f"{nrs.band} default wavelength sampling uses {len(allwaves)} wavelengths")
allwaves

PRISM/CLEAR default wavelength sampling uses 940 wavelengths

$[0.6,~0.605,~0.61,~\dots,~5.285,~5.29,~5.295] \; \mathrm{\mu m}$

# let's get a subset for a faster PSF sim runtime
waves = nrs.get_IFU_wavelengths(nlambda=50)

# convert waves from  microns to meters
# (this is a work around for a current issue with the poppy library upstream)
waves = waves.to_value(u.meter)

Compute a datacube:

cube = nrs.calc_datacube(waves)

The output FITS file has the same extensions as a regular PSF calculation, but each extension contains a 3D datacube rather than a 2D image:

cube.info()

Filename: (No file associated with this HDUList)
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  OVERSAMP      1 PrimaryHDU    1288   (192, 192, 50)   float64
  1  DET_SAMP      1 ImageHDU      1290   (48, 48, 50)   float64
  2  OVERDIST      1 ImageHDU      1343   (192, 192, 50)   float64
  3  DET_DIST      1 ImageHDU      1344   (48, 48, 50)   float64

The calc_datacube_fast routine does a simplified calculation, much faster:

quickcube = nrs.calc_datacube_fast(waves)

Note that in this case, the output FITS file only contains the first oversampled extension:

quickcube.info()

Filename: (No file associated with this HDUList)
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  OVERSAMP      1 PrimaryHDU     159   (192, 192, 50)   float64

The display_psf function works with datacubes, but you have to specify which slice of the cube to display.

index = 20
stpsf.display_psf(cube, ext=3, cube_slice=index,
                    # Note that currently the default plot title isn't very informative for datacube modes
                    # so we can specify a better title directly here:
                    title=f'NRS IFU cube slice {index}, $\\lambda$={cube[0].header["WAVELN"+str(index)]*1e6:.4} micron')

MIRI MRS example

This is mostly the same, except the selection of the IFU bands is done via setting band directly.

miri = stpsf.MIRI()
miri.mode = 'IFU'
miri.band= '2A'
print("Band is", miri.band)

Band is 2A

Note that selecting an IFU band automatically configures the simulated pixelscale to match the default scale used in pipeline output datacubes:

print(f"Pixelscale for band {miri.band} is {miri.pixelscale} arcsec/pix")

miri.band = '3B'
print(f"Pixelscale for band {miri.band} is {miri.pixelscale} arcsec/pix")

Pixelscale for band 2A is 0.17 arcsec/pix
Pixelscale for band 3B is 0.2 arcsec/pix

# let's get a subset for a faster PSF sim runtime
waves = miri.get_IFU_wavelengths(nlambda=50)

# convert waves from  microns to meters
# (this is a work around for a current issue with the poppy library upstream)
waves = waves.to_value(u.meter)

Compute a datacube:

(Note, for MIRI MRS you may see a warning message about the valid region of the field dependence model; this is a benign warning and can be mostly ignored.)

cube = miri.calc_datacube(waves)

/Users/mperrin/Dropbox (Personal)/Documents/software/stpsf/stpsf/opds.py:1759: UserWarning: For (V2,V3) = [-8.40143167 -5.31995333] arcmin, Field point -8.254199999999999 arcmin, -2.4800466666666674 arcmin not within valid region for field dependence model of OTE WFE for MIRI: -8.254199999999999 arcmin--6.21738 arcmin, -2.557224 arcmin--0.5632056 arcmin.  Clipping to closest available valid location, 0.14723166666666643 arcmin away from the requested coordinates.
  warnings.warn(warning_message)

The output FITS file has the same extensions as a regular PSF calculation, but each extension contains a 3D datacube rather than a 2D image:

cube.info()

Filename: (No file associated with this HDUList)
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  OVERSAMP      1 PrimaryHDU     990   (240, 240, 50)   float64
  1  DET_SAMP      1 ImageHDU       992   (60, 60, 50)   float64
  2  OVERDIST      1 ImageHDU      1097   (240, 240, 50)   float64
  3  DET_DIST      1 ImageHDU      1098   (60, 60, 50)   float64

The display_psf function works with datacubes, but you have to specify which slice of the cube to display.

index = 20
stpsf.display_psf(cube, ext=3, cube_slice=index,
                    # Note that currently the default plot title isn't very informative for datacube modes
                    # so we can specify a better title directly here:
                    title=f'MIRI MRS band {miri.band}, cube slice {index}, $\\lambda$={cube[0].header["WAVELN"+str(index)]*1e6:.4} micron')