This notebook shows how to use the METIS LMS mode in Scopesim. The LM Spectrograph (LMS) is the integral-field unit within METIS. Strictly speaking LMS refers to the hardware subsystem and the instrument mode had better been called “IFU”. The use of LMS versus IFU is somewhat inconsistent, it has grown historically and it would be hard to change it now.

from matplotlib import pyplot as plt
from matplotlib.colors import Normalize

import scopesim as sim
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy import units as u

In addition to scopesim we need the configuration and data files that describe the atmosphere, telescope and instrument. These are downloaded from our server with

# Edit this path if you have a custom install directory, otherwise comment it out. [For ReadTheDocs only]
sim.link_irdb("../../../")

# If you want/need to download the instrument packages uncomment the fol
#sim.download_packages(["Armazones", "ELT", "METIS"])

The following command gives the versions of python packages that are used by Scopesim as well as of instrument packages (called IRDB packages). We ask to include the output of this command in any communication with us, in particular in bug reports.

sim.bug_report()

Setting up the instrument

The first step is to set up the instrument in the desired configuration. There are currently two modes that use the LM spectrograph and that are of interest for simulating science observations:

• lms: This is the nominal lms/ifu mode that produces raw two-dimensional detector images in the way that METIS will do. To analyse these data requires a data reduction pipeline, which is not currently available. Scopesim includes a function that allows rectifying these images into a data cube with well-defined spatial and spectral dimensions as we shall demonstrate later in this notebook. This cube corresponds to a single exposure and under-samples the spatial PSF in the across-slice direction.

• lms_cube: This is a convenience mode that directly simulates spectra in a data cube with square spatial pixels that properly sample the spatial PSF. Such a cube will be the result of a procedure that combines a series of dithered and rotated exposures in order to reconstruct the PSF, which is under-sampled in individual exposures. While the simulation contains many of the noise contributions present in real data it cannot include the uncertainties related to an actual data reduction procedure and is therefore rather optimistic.

The extended or spectral IFU mode is not yet available in Scopesim.

We start by setting up a UserCommands object in which we specify the instrument to be used (METIS), the instrument mode (lms) as well as other parameters needed to uniquely specify the instrument setup. For the LMS, a parameter that will always have to be set is the target wavelength. In Scopesim, this can be given simply as a float (units of micrometers).

Using the nominal LMS mode

cmd = sim.UserCommands(use_instrument="METIS", set_modes=["lms"], properties={"!OBS.wavelen": 3.5})

Further parameters can be set either in the command above or separately as in

cmd["!OBS.interp_psf"] = False    # only scale PSF to the central wavelength

The reason for doing this will be explained below.

The OpticalTrain (consisting of atmosphere + telescope + instrument) is then built from the UserCommands object (DeprecationWarnings can be ignored):

metis = sim.OpticalTrain(cmd)

An OpticalTrain consists of a list of Effect objects, which can be listed as follows:

metis.effects.pprint_all()

Observation of blank sky

We start by simply observing a piece of blank sky. This way we don’t have to worry about setting up an astronomical source, which we will however do shortly. IFU simulations are fairly complex and therefore need some time. Much of that time is spent on PSF convolution, which has no effect when the target fills the field of view homogeneously as blank sky does. One could turn off PSF convolution entirely by doing metis['psf'].include = False but one needs to be careful to turn it back on again when observing an actual source. Another way to speed things up is not to scale the PSF by wavelength - this is indeed admissable in LMS observations which cover a short wavelength range (one should not do in long-slit simulations where the PSF at the red end is nearly twice as large as at the blue end). We have already switched off PSF scaling (except to the central wavelength) above.

metis.observe()

The result of observe() is an ImagePlane object. This is a noise-free image of the spectral layout in the detector plane. The units are photons per second (per pixel). A readout() is needed to put this image onto the detector array and add photon noise, dark current and other detector effects to it. We first look at the ImagePlane, then create a readout with exposure time 600 seconds.

sky_implane = metis.image_planes[0].hdu
sky_implane.writeto("sky_implane.fits", overwrite=True)
plt.imshow(sky_implane.data, origin="lower")
plt.colorbar()

The ImagePlane then needs to be read out by the detectors to add all the relevant noise components. This can be repeated as often as desired, for testing different exposure times or for creating exposure sequences with the same exposure time. readout() also applies gain conversion, so the output of a readout is ADU. The gain for the LMS detectors is currently set at 2 electrons/ADU.

sky_read = metis.readout(dit=600, ndit=1)[0]

norm = Normalize(vmin=np.min([hdu.data for hdu in sky_read[1:]]),
                 vmax=np.max([hdu.data for hdu in sky_read[1:]]))
fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(6, 6))
axes[0, 0].imshow(sky_read[1].data, origin='lower', norm=norm)
axes[0, 1].imshow(sky_read[2].data, origin='lower', norm=norm)
axes[1, 0].imshow(sky_read[3].data, origin='lower', norm=norm)
axes[1, 1].imshow(sky_read[4].data, origin='lower', norm=norm)

# Save the sky observation for later use
sky_read.writeto("sky_lms-3.5um-600s.fits", overwrite=True)

Observation of a source

We take the source from the fits file that we created in the previous notebook, METIS_IFU-01-Source_cube.ipynb.

cubefits = "our_first_cube.fits"

with fits.open(cubefits) as hdul:
    src = sim.Source(cube=hdul)

metis.observe(src)

src_implane = metis.image_planes[0].hdu
src_implane.writeto("src_implane.fits", overwrite=True)
plt.imshow(src_implane.data - sky_implane.data, origin='lower')
plt.colorbar()

src_read = metis.readout(dit=600, ndit=1)[0]

# If src and sky had been read out with different dit/ndit settings, sky would have to be scaled
for i in range(1, 5):
    src_read[i].data = src_read[i].data - sky_read[i].data

norm = Normalize(vmin=np.min([hdu.data for hdu in src_read[1:]]),
                 vmax=np.max([hdu.data for hdu in src_read[1:]]))
fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(6, 6))
axes[0, 0].imshow(src_read[1].data, origin='lower', norm=norm)
axes[0, 1].imshow(src_read[2].data, origin='lower', norm=norm)
axes[1, 0].imshow(src_read[3].data, origin='lower', norm=norm)
axes[1, 1].imshow(src_read[4].data, origin='lower', norm=norm)

src_read.writeto("hl_tau_lms_600s.fits", overwrite=True)

rectified = metis['lms_spectral_traces'].rectify_cube(src_read)

We first show an image reconstruction of the cube by summing over the wavelength axis. Since the spatial pixels are rectangular we have to use the aspect keyword to get the correct aspect ratio.

wcs = WCS(rectified.header)
wcs.sub(2)     # celestial part of the WCS (linear for now)
fig, ax = plt.subplots(subplot_kw=dict(projection=wcs.sub(2)), figsize=(6, 3))
# `aspect` is needed because pixels are rectangular!
img = ax.imshow(rectified.data.sum(axis=0), origin="lower", aspect=20.7/8.2)
fig.colorbar(img, label="ADU");

The global spectrum is obtained by summing over the spatial axes. The wavelength values are obtained via the spectral part of the WCS.

swcs = wcs.spectral
wavelength = swcs.all_pix2world(np.arange(rectified.data.shape[0]), 0)[0]
wavelength *= u.Unit(swcs.wcs.cunit[0]).to(u.um)
plt.plot(wavelength, rectified.data.sum(axis=(1, 2)))
plt.xlabel("Wavelength [um]")
plt.ylabel("ADU");

Using the LMS Cube mode

The rectified spectral cube shown above has rectangular spaxels, with detector pixels of 8.2 mas along the x-direction and slices of 20.7 mas along the y-direction. When METIS is on sky a sequence of dithered and rotated exposures will be taken, which will allow reconstruction of the fully sampled field on 8.2 mas pixels in both directions. For convenience, Scopesim provides a mode lms_cube that simulates these reconstructed cubes directly. The noise levels are realistic, but the uncertainties connected to the reconstruction procedure are not.

cmd = sim.UserCommands(use_instrument="METIS", set_modes=['lms_cube'], properties={"!OBS.wavelen": 3.5})

cmd["!SIM.spectral.spectral_bin_width"] = 1e-5
cmd["!OBS.interp_psf"] = False

metis = sim.OpticalTrain(cmd)

Note: `observe()` in cube mode requires a lot of memory for PSF convolution. Since the PSF is not crucial for our extended source we turn it off here. You can include it by changing `False` to `True` in the following cell.

metis["psf"].include = False
metis.observe(src)

metis.image_planes[0].data.shape

src_cube = metis.readout(dit=600, ndit=1)[0][1]

readout() results in a 3D FITS cube with a linear wavelength scale and two spatial dimensions. We will look at the spatial and spectral dimensions, respectively summed over the remaining dimensions.

wcs = WCS(src_cube)
fig, ax = plt.subplots(subplot_kw=dict(projection=wcs.sub(2)))
img = ax.imshow(src_cube.data.sum(axis=0), origin='lower')
fig.colorbar(img, label="ADU");

swcs = wcs.spectral
wavelength_cube = swcs.all_pix2world(np.arange(src_cube.data.shape[0]), 0)[0]
wavelength_cube *= u.Unit(swcs.wcs.cunit[0]).to(u.um)

_, ax = plt.subplots(figsize=(10, 6))
ax.plot(wavelength_cube, src_cube.data.sum(axis=(1, 2)))
ax.set_xlabel("Wavelength [µm]")
ax.set_ylabel("ADU");

Simulate an empty sky observation for background subtraction.

metis.observe()

sky_cube = metis.readout(dit=600, ndit=1)[0][1]

_, ax = plt.subplots(figsize=(12, 6))
ax.plot(wavelength_cube, (src_cube.data - sky_cube.data).sum(axis=(1, 2)), label="Cube mode")
ax.plot(wavelength, rectified.data.sum(axis=(1, 2)), label="Nominal mode")
ax.legend()
ax.set_xlabel("Wavelength [µm]")
ax.set_ylabel("ADU");

The signal in cube mode is higher than that in the nominal mode as the field of view is larger ($0.9\times 0.9$ arcsec$^2$ compared to $0.9\times0.6$ arcsec$^2$. The detector gap at 3.5 um is not taken into account in cube mode. The continuum is modified by the instrumental transmission (in particular the grating efficiency); atmospheric absorption lines are clearly visible throughout. Correction of these effects would require simulation of flux and telluric standard stars.