MOSAIC Demo

MOSAIC Demo#

This is a simple demo of the “ETC” mode for MOSAIC. To use this you need the oc/mosaic branch of Scopesim and the oc/mosaic branch of the irdb. The implementation creates one-dimensional spectra for the output (summed over the fibres in a bundle for the MOS modes). The first sections demonstrate how to use Scopesim for the MOS and mIFU modes. The final section provides a more detailed description of the implementation and instrument definition that is currently used.

import scopesim as sim
import numpy as np
from matplotlib import pyplot as plt
from astropy import units as u

Change the path below to your local irdb copy.

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

To simulate observations with MOSAIC, the instrument packages for MOSAIC, the ELT, and Armazones are required. The packages are downloaded from the server and installed into sub-directory inst_pkgs in the current working directory. If you have not got the packages yet, uncomment the following cell.

#sim.download_packages(["MOSAIC", "ELT", "Armazones"])

The following command gives an overview of installed packages used by Scopesim. Please include the output in every bug report or question.

sim.bug_report()

Instrument modes#

cmd = sim.UserCommands(use_instrument="MOSAIC")

cmd.modes

  • The visual modes map spectra to a pseudo-detector of length 13000 pixels; the gap between the two 6k detectors in the real instrument is not simulated.

  • The near-infrared modes map spectra to a 4k x 4k detector.

  • As stated above, the MOS modes collapse the 7 (low-res modes) or 19 (high-res modes) fibres in a bundle to a single one-dimensional spectrum, the output is a FITS binary table with wavelength and flux (in ADU).

  • The mIFU modes return a 4k x 4k image with each of the 4xx fibre spectra mapped onto a detector row each. Rearrangement of the spectra into a cube has not been implemented.

MOS modes#

All the MOS modes work identically, and we will only demonstrate the low-resolution R band mode.

cmd = sim.UserCommands(use_instrument="MOSAIC", set_modes=['MOS-LR-R'])

We shall use an exposure time of 10 seconds throughout. This could be set in the Usercommands object (cmd["!OBS.exptime"]) but we prefer to set it explicitely in the readout commands.

t_exp = 10   # seconds

We can now build the instrument model as an instance of the OpticalTrain class.

mosaic = sim.OpticalTrain(cmd)

The following gives an overview of the effects that are included in the optical train.

mosaic.effects

We start by observing blank sky; this will be useful for background subtraction when we observe a source later.

mosaic.observe()

This step creates an ImagePlane object, a noise-less image of the detector that records the expected number of electrons per second per pixel. In the current implementation for the visual branch of MOSAIC this is a pseudo-detector with 13000 pixels in the dispersion direction (x-axis), and just enough pixels in the y-direction to hold the 19 fibres of the HR modes (only 7 for LR modes, as in the example). Each fibre is mapped to a single row in the ImagePlane. Normally, the ImagePlane is not of interest to the user; we save it here for later when we will use to estimate signal-to-noise ratios.

sky_implane = mosaic.image_planes[0].data
print(f"Size of imageplane: {sky_implane.shape}")
plt.imshow(sky_implane[:, 2000:3000], origin='lower', vmax=0.1);

The actual detector “image” including dark current, shot noise and readout noise is created in the next step, where we have to specify the exposure time.

sky = mosaic.readout(exptime=t_exp)[0]

The result is a binary table with columns wavelength and spectrum:

sky[1].data

plt.plot(sky[1].data['wavelength'], sky[1].data['spectrum'])

As a source, we’ll use a star of 15 mag:

src = sim.source.source_templates.star(flux=15 * u.mag)

mosaic.observe(src)

Again, we save the ImagePlane for later. Note that this contains the expected electron counts (per second) from both the source and the background. We then proceed to create the actual detector readout.

star_implane = mosaic.image_planes[0].data

star = mosaic.readout(exptime=t_exp)[0]      # same exposure time as for the background simulation

plt.plot(star[1].data['wavelength'], star[1].data['spectrum'], label="raw spectrum", lw=1)
plt.plot(star[1].data['wavelength'], star[1].data['spectrum'] - sky[1].data['spectrum'], label="background subtracted", lw=1)
plt.plot(sky[1].data['wavelength'], sky[1].data['spectrum'], label="background spectrum", lw=1)
plt.legend()
plt.text(0.656281, 1700, r"H$\alpha$", horizontalalignment="center")
plt.text(0.686719, 1700, r"O$_2$", horizontalalignment="center")
plt.text(0.72, 2500, r"H$_2$O", horizontalalignment="center")
plt.text(0.759370, 450, r"O$_2$", horizontalalignment="center")
plt.text(0.822696, 1300, r"O$_2$", horizontalalignment="center")
plt.xlabel("Wavelength [um]")
plt.ylabel("Electrons");

Strictly speaking, the flux units is ADU, but the gain is currently set to 1, so that 1 ADU corresponds to 1 electron collected over the exposure time (60 s). The spectrum is the sum over the seven fibres (19 for the high-resolution modes), each covering a solid angle of 0.0118 arcsec$^2$.

Note that the mosaic.readout() method above has been set by the exptime keyword, and the exposure is a single readout of that integration time. This is appropriate for the visual modes. The near-infrared modes can also receive a simple exptime. In this case, the exposure time is automatically split into ndit subexposures of integration time dit each (cf. the readout of the mIFU example below). It is also possible to set dit and ndit directly.

Signal-to-noise estimation from the simulated data#

Scopesim does not have a built-in function to compute a signal-to-noise ratio for the simulation. Using the simulations for the background and star+background, it can be easily estimated.

To start with we have to convert the spectra from ADU to electrons by multiplying by the gain (in the visual mode, this is currently set to 1; in the near-infrared it is set to 2.5 from our previous experience with HxRG detectors).

gain = mosaic.cmds["!DET.gain"]
print(f"Gain: {gain}")

noise = np.sqrt(star[1].data['spectrum'] / gain)
signal = (star[1].data['spectrum'] - sky[1].data['spectrum'])/ gain
SNR_emp = signal/noise

plt.plot(star[1].data['wavelength'], SNR_emp)
plt.ylim(0, 80)
plt.xlabel("Wavelength [um]")
plt.ylabel("Signal-to-noise ration")
plt.title("Empirical signal-to-noise ratio");

Signal-to-noise ration determination from expected values#

We can also break down the noise into its constituent components as they are currently implemented in Scopesim. For the signal and the Poisson noise due to source and background we go back to the expected values that are provided in the ImagePlane of the respective simulations. We use the following prescription for the SNR computation (ELT Spectroscopy ETC document): $$ \frac{S}{N} = \frac{\sqrt{n_\mathrm{exp}}\cdot N_\mathrm{obj}}{\sqrt{N_\mathrm{obj} + N_\mathrm{sky} + n_\mathrm{pix} R^2 + n_\mathrm{pix} D t_\mathrm{exp}}} $$ For the electron numbers $N_\mathrm{obj}$ and $N_\mathrm{sky}$ from the source and the background, respectively, we recur to the ImagePlanes that we saved earlier. These give the expected number of electrons per second per pixels.

N_obj = (star_implane - sky_implane).sum(axis=0) * t_exp    # sum over all fibres (non-fibre rows are zero, so summing over y-axis is fine)
N_sky = sky_implane.sum(axis=0) * t_exp

Scopesim uses $R=7$ electrons and $D=0.005$ electrons per second. These values were taken from E-MOS-SYS-ANR-0063-2_0, where they are given as $R\approx 2.5$ electrons and $D=3$ electrons per hour, respectively. Our adopted values take into account that each row in the ImagePlane is the sum over a fibre, which we take to be around 6 pixels wide in a (real) detector image ($D$ scales linearly with that number, while $R$ scales with its square root). In the formula above $n_\mathrm{pix}$ corresponds to the number of fibres that have been summed over, which is 7 in our case. The number of exposures is $n_\mathrm{exp} = 1$ in the visual.

rdnoise = mosaic.cmds["!DET.readout_noise"]
dark = mosaic.cmds["!DET.dark_current"]
print(f"Readout noise: {rdnoise} electrons")
print(f"Dark current:  {dark} electrons/second")

SNR_exp = N_obj / np.sqrt(N_obj + N_sky + 7 * rdnoise**2 + 7 * dark * t_exp)

plt.plot(star[1].data['wavelength'], SNR_exp)
plt.xlabel("Wavelength [um]")
plt.ylabel("Signal-to-noise ration")
plt.title("Expected signal-to-noise ratio")
plt.ylim(0, 80);

Using a galaxy as a source#

The scopesim_templates package includes models for a variety of astronomical source. We will use it to simulate an observation of an elliptical galaxy at some redshift.

import scopesim_templates as sim_tp

gal = sim_tp.extragalactic.galaxy("kc96/elliptical", z=0.1, amplitude=15, filter_curve="g", pixel_scale=0.01, r_eff=0.3, 
                                  n=4, ellip=0.5, theta=45, extend=3)

plt.imshow(gal.fields[0].data, origin="lower", norm="log")

gal.fields[0].header

cmd_b = sim.UserCommands(use_instrument="MOSAIC", set_modes=["MOS-LR-B"])
mosaic_b = sim.OpticalTrain(cmd_b)

mosaic_b.observe(gal)

galobs = mosaic_b.readout(exptime=t_exp)[0]

plt.plot(galobs[1].data['wavelength'], galobs[1].data['spectrum'],  label="galaxy")
plt.plot(star[1].data['wavelength'], star[1].data['spectrum'] - sky[1].data['spectrum'], label="star")
plt.xlabel("Wavelength [um]")
plt.legend();

mIFU modes#

cmd = sim.UserCommands(use_instrument="MOSAIC", set_modes=["mIFU-LR-J"])

mosaic = sim.OpticalTrain(cmd)

mosaic.observe()

hdul = mosaic.readout(exptime=3600)[0]

plt.imshow(hdul[1].data[900:1100, 1900:2100], norm='log', origin='lower')
plt.colorbar();

plt.plot(hdul[1].data[991, :])

The mIFU mode is incomplete - it is currently not possible to deduce the wavelength vector and the spatial arrangement of the fibres from the readout alone. A table format might be better suited for that purpose.

Instrument package#

In this section we describe how the instrument is defined for Scopesim and what data files and parameters are used. The instrument configuration is found in irdb/MOSAIC (a release package will be provided once the code and configuration has been validated). The main file is default.yaml where in particular the available instrument modes are defined (see above for a list of these modes). The OpticalTrain, i.e. the software representation of the entire system of atmosphere, telescope, instrument optics and detector, is built as a series of Effect objects, which are set along with relevant parameters in further yaml files, as appropriate for each instrument mode. All modes use the additional Armazones and ELT packages:

  • irdb/Armazones/Armazones.yaml sets a single effect, skycalc_atmosphere, which provides atmospheric transmission and emission spectra as taken from ESO’s skycalc server. The two parameters spectral_resolution and spectral_bin_width have been carefully chosen for each mode and should not be changed (unless for experimental purposes).

  • irdb/ELT/ELT.yaml defines telescope_reflection, which provides reflectivity and thermal emission from the telescope mirrors. The parameter !TEL.ter_curve.filename is by default set to TER_ELT_5_mirror.dat. The irdb contains files for 5- and 6-mirror configurations (the latter including a flat mirror to deflect light to a side port), the latter with exposed (field-tracking) or hidden (pupil-tracking) telescope spiders, and for clean and not so clean mirror segments.

The instrument is split into a visual and a near-infrared arm, described by MOSAIC_VIS.yaml and MOSAIC_NIR.yaml, respectively.

  • The transmission of the visual arm is given as a single combined system transmission effect, which uses the file named in the parameter “!OBS.ter_file” (TER = transmission, emission, reflection). For the near-infrared arm, the transmission is split into fibre transmission and spectrograph transmission (parameters !OBS.fibre_ter_file and !OBS.spec_ter_file). The fibre transmission files are named TER_mos_* to distinguish them from the IFU fibres, TER_mIFU_*. The available data are rather crude at present.

  • The point-spread function (effect psf) is currently taken to be a SeeingPSF, with FWHM settable via !OBS.psf_fwhm. This parameter is currently set to a default of 0.2 arcsec for all modes, probably not always a good value.

  • The line-spread function (effect lsf) is a convolution of a top-hat of width 4.5 spectral bins in the visual and 2.75 spectral bins in the near-infrared (the effect parameter lsfwidth is not currently user-settable, except by editing the yaml file; this should change in the future) with a Gaussian with standard deviation of 1 spectral bin. A spectral bin is given by the dispersion (um per pixel) on the detector.

  • The fibre_bundle effect determines the outline of the fibre bundle on the sky, and the spectral_traces effect is responsible for mapping the spectra into the two-dimensional ImagePlane (i.e. the detector focal plane). Both effects use the same trace file. Trace files are available for all instrument modes as multi-extension FITS files. The spatial arrangement of the fibres in a bundle is given as a table in the second extension of the file. Note that the fibre apertures are modelled as squares rather than hexagons (this makes it easier to sample the source cube); of course the squares cover the same solid angle on sky as the hexagons.

The remaining traces map the spectra into the detector. This is done in a strictly linear manner, i.e. with constant dispersion (different for each mode). Each fibre is mapped as a one-dimensional spectrum onto a single row of the ImagePlane.

The ImagePlane is produced by the mosaic.observe() command and is an image of the expected flux in photon/s in the detector plane, without noise. It is accessible by mosaic.image_plane.hdu.

The mosaic.readout() command subsequently invokes the detector configuration to create a noisy detector readout from the ImagePlane. The detectors are defined in MOSAIC_DET_VIS.yaml and MOSAIC_DET_NIR.yaml, respectively. The geometric layout of the detectors is given by the effect detector_array, which reads FPA_mosaic_VIS_layout.dat and FPA_mosaic_NIR_layout.dat, respectively. The visual detector is currently described as a single 13000 x 160 array rather than the four 6k x 6k detectors foreseen for the actual instrument (hence the “s(i)mpl(e)” in the yaml file name) . This is because Scopesim only simulates a single fibre bundle, so a narrow effective detector is sufficient. Masking the gap between the two detectors is left to the user. The NIR detector is modeled as a single 4k x 4k detector (as in the actual instrument) to have enough room for the mIFU mode.

The QE is handled by the effect quantum_efficiency, which uses the files QE_VIS_ML2_depleted.dat and QE_detector_NIR.dat, respectively.

Finally, both detector yamls include the effect collapse_1d, which sums up the separate spectra from the fibres and produces a single output spectrum as described above. It is possible to switch this off with mosaic['collapse_1d'].include=False to see the “detector image”. However, due to the way the spectra are currently formed, this is probably not of much interest.