Field of Stars — Liger Imager

from liger_iris_sim.sources import make_point_source_image
from liger_iris_sim.expose import compute_throughput, expose_imager
from liger_iris_sim.utils import LIGER_PROPS, compute_filter_photon_flux, get_psf
from liger_iris_sim.sky import get_maunakea_sky_background
from liger_iris_drp_resources.filters import load_filters_summary

import numpy as np
import matplotlib.pyplot as plt

Instrument and exposure parameters

instrument_name = 'Liger'
instrument_mode = 'img'          # imager mode
filt = 'J'
size = (2048, 2048)     # detector size in pixels
read_noise = 9        # e- RMS
resolution = 4000
dark_current = 0.025  # e- / sec / pixel
scale = 0.01          # arcsec / pixel
itime = 300           # integration time, sec
n_frames = 1          # number of coadded frames
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'],
    tel=0.8, ao=0.65, filt=filt
)
print(f"Total throughput for {instrument_mode} mode at {filter_info['wavecenter']:.3f} µm: {tput:.3f}")

Get the PSF

get_psf returns the PSF at a given wavelength and field position binned on the detector or lenslet/slicer grid.

psf, psf_info = get_psf(
    instrument_name=instrument_name,
    instrument_mode=instrument_mode,
    wave=filter_info['wavecenter'],
    xs=0, ys=0,
    output_plate_scale=scale
)
print(f"PSF shape after rebinning: {psf.shape}")

Create a field of stars

compute_filter_photon_flux converts a magnitude to photons / sec / m² using the filter zero-point. make_point_source_image places each star into the source rate image, convolved with the PSF.

np.random.seed(1)
n_stars = 10
mag_range = (18, 19)

source_rate = np.zeros(size, dtype=np.float32)

for i in range(n_stars):
    mag = np.random.uniform(mag_range[0], mag_range[1])
    xpix = int(np.random.uniform(5, size[1] - 5))
    ypix = int(np.random.uniform(5, size[0] - 5))
    photon_flux = compute_filter_photon_flux(mag, zp=filter_info['zpphot'])
    print(f"Star {i+1:2d}: mag={mag:.2f}, x={xpix}, y={ypix}, flux={photon_flux:.2e} phot/s/m²")
    make_point_source_image(
        xpix, ypix, photon_flux,
        psf=psf,
        image_out=source_rate,
    )

Sky background

The per-pixel sky emission rate is derived from the filter background magnitude and integrated over the pixel solid angle.

sky_data = get_maunakea_sky_background(
    filter_info=filter_info,
    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_rate_bandpass_tot']
sky_trans = sky_data['sky_trans_bandpass_mean']

Run the exposure simulation

expose_imager adds Poisson noise, dark current, and read noise, and returns a dict of simulated images and noise components.

sim = expose_imager(
    source_rate,
    itime=itime, n_frames=n_frames, collarea=collarea,
    sky_emission_rate=sky_em_rate,
    tput=tput, read_noise=read_noise, dark_current=dark_current,
)
print(f"Peak SNR: {np.nanmax(sim['snr']):.1f}")

Visualize results

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
rate_map = sim['observed_rate'].copy()
rate_map[rate_map <= 0] = np.nan  # Mask non-positive values for
im0 = axes[0].imshow(
    np.log10(rate_map),
    origin='lower', cmap='inferno',
)
axes[0].set_title('log10(Observed rate (e-/s))')
plt.colorbar(im0, ax=axes[0])

im1 = axes[1].imshow(sim['snr'], origin='lower', cmap='viridis', vmin=0)
axes[1].set_title('SNR per pixel')
plt.colorbar(im1, ax=axes[1])

plt.tight_layout()
plt.show()