Simulate spectrograms#

This example simulates a single emission line galaxy in one grism exposure, then extracts its spectrum and builds an emission line map. You only need the SIR calibration files (see Core concepts); no Euclid observations are used.

Set up a frame#

Create a gelsa.Gelsa instance and an empty spectroscopic frame at a chosen pointing:

import numpy as np
from matplotlib import pyplot as plt

import gelsa
from gelsa import analysis, galaxy, visu, utils

G = gelsa.Gelsa(calibdir="path/to/SIR/calibration")

frame = G.new_spec_frame(ra=0, dec=0, pa=0, grism_name="RGS000")

Choose the source position by asking where pixel (1000, 1000) on detector 5 sees light at 1.5 µm:

ra, dec = frame.pixel_to_radec(1000, 1000, 5, 15000)

Define a galaxy#

A Galaxy has a position, a redshift, a light profile, a continuum and a set of emission lines. Here it has a Gaussian profile and an Hα line (with [NII]) at redshift 1.2:

redshift = 1.2

gal = galaxy.Galaxy(
    ra=ra, dec=dec,
    redshift=redshift,
    fwhm_arcsec=0.5,
    axis_ratio=1,
)
gal.set_flux_Ha(2e-15, NII_Ha_ratio=0.2)   # erg/s/cm²

Render the galaxy into the frame. With noise=True, detector noise is added:

frame.add_sources([gal], noise=True)

Look at the spectrogram#

cutout() returns a dictionary of SpecCrop objects, one for each detector the trace falls on. gelsa.visu.show() displays a crop:

crops = frame.cutout(ra, dec, redshift=redshift, padx=5, pady=25)
crop = crops[5]
visu.show(crop)

Extract the spectrum#

Analysis works the same way on one frame as on many. Add the frame, cut out the source, and resample it onto a wavelength grid:

A = analysis.Analysis(G)
A.add_frame(frame)
A.crop(ra, dec, redshift=redshift, padx=5, pady=25)

spec = A.build_stacked_2D_spectrum(wave_range=(12000, 19000))

spec is a dictionary. Its main entries are:

  • wavelength: the wavelength grid (Å),

  • spec2d, spec2d_var: the 2D spectrum and its variance,

  • spec1d, spec1d_var: the extracted 1D spectrum and its variance,

  • pixel_perp: the offset from the trace, in pixels.

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)

extent = (spec["wavelength"][0], spec["wavelength"][-1],
          spec["pixel_perp"][0], spec["pixel_perp"][-1])
spec2d = np.ma.array(spec["spec2d"], mask=spec["spec2d_norm"] == 0)
ax1.imshow(visu.normalize_image(spec2d), extent=extent, aspect="auto")

ax2.plot(spec["wavelength"], spec["spec1d"])
ax2.set_xlabel("Wavelength (Å)")

The Hα line appears at 6564.6 × (1 + z) ≈ 14442 Å.

Build an emission line map#

A line map resamples the dispersed light at one wavelength back onto the sky. Define the output grid with an astropy.wcs.WCS, here 41 × 41 pixels of 0.3″ centred on the source:

wcs = utils.make_wcs(ra, dec, pixel_size_arcsec=0.3, shape=(41, 41))

line_wavelength = 6564.61 * (1 + redshift)
line_map, line_var, norm, count = A.build_stacked_line_map(
    line_wavelength, wcs, width=50)

plt.subplot(projection=wcs)
plt.imshow(line_map, origin="lower")

Many exposures#

To simulate a dithered observation, create several frames with slightly different pointings or grisms, add the same galaxy to each, and pass all of them to the same Analysis before calling crop:

A = analysis.Analysis(G)
for grism in ("RGS000", "RGS180"):
    f = G.new_spec_frame(ra=0, dec=0, pa=0, grism_name=grism)
    f.add_sources([gal], noise=True)
    A.add_frame(f)

A.crop(ra, dec, redshift=redshift, padx=5, pady=25)
spec = A.build_stacked_2D_spectrum(wave_range=(12000, 19000))