Field of Stars — Liger IFS

from liger_iris_sim.sources import make_point_source_ifs_cube
from liger_iris_sim.expose import compute_throughput, expose_ifs
from liger_iris_sim.sky import get_maunakea_sky_background
from liger_iris_sim.utils import LIGER_PROPS, compute_filter_photon_flux, generate_wave_grid_for_filter, get_psfs

from liger_iris_drp_resources.filters import load_filters_summary

import numpy as np
import matplotlib.pyplot as plt

Instrument and exposure parameters

np.random.seed(1)

instrument_name = 'Liger'
instrument_mode = 'ifs'
ifs_mode = 'lenslet'  # 'lenslet' or 'slicer'
filt = 'J'
size = (128, 128)       # IFS spatial dimensions (y, x) in spaxels
read_noise = 9        # e- RMS
dark_current = 0.025  # e- / sec / pixel
scale = 0.014         # arcsec / spaxel (lenslet)
itime = 30           # integration time, sec
n_frames = 1          # number of coadded frames
resolution = 4000     # spectral resolution R = lambda / delta_lambda
collarea = LIGER_PROPS['collarea']  # m^2

Load filter data and compute throughput

load_filters_summary returns filter metadata including the central wavelength, zero-point flux, and background magnitude.

compute_throughput folds together telescope, AO, filter, and instrument throughput.

filter_info = load_filters_summary(filter_name=filt)

tput = compute_throughput(
    instrument_name=instrument_name,
    instrument_mode=instrument_mode,
    wave=filter_info['wavecenter'],
    ifs_mode=ifs_mode,
    filt=filt,
)
print(f"Total throughput for {ifs_mode} mode at {filter_info['wavecenter']:.3f} µm: {tput:.3f}")

Build the wavelength grid

generate_wave_grid_for_filter returns a Nyquist-sampled wavelength array covering the filter bandpass at the requested spectral resolution.

wave = generate_wave_grid_for_filter(filter_info, resolution=resolution)
print(f"Wavelength grid: {wave[0]:.4f} – {wave[-1]:.4f} µm, {len(wave)} channels")

Define field of stars and spectrum

n_stars = 10
mag_range = (16, 17)

mag = np.random.uniform(mag_range[0], mag_range[1], size=n_stars)
xpix = np.random.uniform(5, size[1] - 5, size=n_stars)
ypix = np.random.uniform(5, size[0] - 5, size=n_stars)
photon_flux = compute_filter_photon_flux(mag, zp=filter_info['zpphot'])  # photons / sec / m^2

# Flat spectrum normalised to unit integral
flat_spec = np.ones(len(wave), dtype=np.float32)
template = (wave, flat_spec / flat_spec.sum())

Load PSF

psfs, psf_infos, psf_indices = get_psfs(
    instrument_name=instrument_name,
    instrument_mode=instrument_mode,
    wave=filter_info['wavecenter'],
    xdet=xpix, ydet=ypix,
    output_plate_scale=scale,
)

Compute sky emission and transmission

sky_data = get_maunakea_sky_background(
    wave, resolution=resolution,
    T_tel=275, T_atm=258, T_aos=243,
    Em_tel=0.09, Em_atm=0.2, Em_aos=0.01,
    plate_scale=scale,
)

sky_em_rate = sky_data['sky_em']
sky_trans = sky_data['sky_trans']

Create a field of stars

input_cube_rate = np.zeros((len(wave), *size), dtype=np.float32)

make_point_source_ifs_cube(
    xdet=xpix, ydet=ypix, wave=wave,
    template=template,
    flux_int=photon_flux,
    psf=psfs,
    psf_indices=psf_indices,
    cube_out=input_cube_rate,
)

Run the IFS exposure simulation

expose_ifs applies throughput, sky emission, sky transmission, Poisson noise, dark current, and read noise to produce a simulated data cube.

sim = expose_ifs(
    input_cube_rate,
    itime=itime, n_frames=n_frames, collarea=collarea,
    sky_emission_rate=sky_em_rate,
    sky_transmission=sky_trans,
    tput=tput, read_noise=read_noise, dark_current=dark_current,
)

print(f"Output cube shape: {sim['observed_tot'].shape}")
print(f"Peak SNR (per voxel): {np.nanmax(sim['snr']):.1f}")

Visualize results

Collapse the cube along the wavelength axis to produce a white-light image, and show the peak-SNR spectrum.

snr_map = np.nansum(sim['snr']**2, axis=0)**0.5
im = plt.imshow(snr_map, origin='lower', cmap='viridis', vmin=0)
plt.colorbar(im, label='SNR')
plt.title('SNR integrated over bandpass')
plt.show()

# Extract and plot the spectrum of the brightest spaxel
peak_idx = np.unravel_index(np.nanargmax(snr_map), snr_map.shape)
ypeak, xpeak = peak_idx

fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

axes[0].plot(wave, sim['observed_rate'][:, ypeak, xpeak], lw=1, label='Observed')
axes[0].plot(wave, sim['source_rate'][:, ypeak, xpeak], lw=1, ls='--', label='Source (noiseless)')
axes[0].set_ylabel('Signal (e-/sec)')
axes[0].legend()
axes[0].set_title(f'Spectrum at brightest spaxel (x={xpeak}, y={ypeak})')

axes[1].plot(wave, sim['snr'][:, ypeak, xpeak], color='C2', lw=1)
axes[1].set_ylabel('SNR')
axes[1].set_xlabel('Wavelength (µm)')
plt.show()