Stellar Velocity Mapping — 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

Filter, throughput, and wavelength grid

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}")

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 emission line sources for simplicity instead of stellar spectrum

Each star gets a Gaussian spectral template centred at the filter wavelength. A small velocity field is imposed by shifting the template peak radially from the field centre.

n_stars = 100
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

# Gaussian spectrum template (normalised to unit integral)
def gauss(x, A, mu, sigma):
    return A * np.exp(-(x - mu)**2 / (2 * sigma**2))

spec = gauss(wave, A=1, mu=filter_info['wavecenter'], sigma=filter_info['wavecenter'] / resolution)
template_flux = spec / spec.sum()
template = (wave, template_flux)

# Radial velocity field: z increases with distance from field centre
r = np.sqrt((xpix - xpix.mean())**2 + (ypix - ypix.mean())**2)
zs = 0.002 * (r / r.max())

print(f"Injected redshifts: {zs.min():.5f} – {zs.max():.5f}")

Load PSFs

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,
)

Sky background

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']

Build source cube and simulate exposure

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,
    z=zs,
)

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}")

SNR map

snr_map = np.nansum(sim['snr']**2, axis=0)**0.5

plt.figure()
im = plt.imshow(snr_map, origin='lower', cmap='viridis', vmin=0)
plt.colorbar(im, label='SNR')
plt.title('SNR integrated over bandpass')
plt.show()

Velocity map recovery

For each star the peak-flux spaxel is used. The flux-weighted centroid of the spectrum within ±6 channels of the peak gives the observed wavelength; converting to redshift recovers the injected velocity field.

zs_rec = np.zeros(n_stars, dtype=np.float32)

for i in range(n_stars):
    ii, jj = int(round(ypix[i])), int(round(xpix[i]))
    spec_i = input_cube_rate[:, ii, jj]
    if spec_i.sum() > 0:
        k = np.argmax(spec_i)
        ki = max(k - 6, 0)
        kf = min(k + 6, len(wave) - 1)
        ws = wave[ki:kf + 1]
        wt = spec_i[ki:kf + 1]
        wc = np.sum(ws * wt) / np.sum(wt)
        zs_rec[i] = (wc - filter_info['wavecenter']) / filter_info['wavecenter']

# Convert to km/s
c_kms = 2.998e5
v_true = zs * c_kms
v_rec  = zs_rec * c_kms

print(f"Residuals (km/s): {(v_rec - v_true)}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

# --- velocity map on the sky ---
sc = axes[0].scatter(xpix, ypix, c=v_rec, cmap='RdBu_r', s=80,
                     vmin=v_true.min(), vmax=v_true.max())
plt.colorbar(sc, ax=axes[0], label='Recovered velocity (km/s)')
axes[0].set_xlim(0, size[1])
axes[0].set_ylim(0, size[0])
axes[0].set_xlabel('x (spaxels)')
axes[0].set_ylabel('y (spaxels)')
axes[0].set_title('Velocity map')

# --- injected vs recovered ---
axes[1].plot(v_true, v_rec, 'o')
vmin, vmax = v_true.min(), v_true.max()
axes[1].plot([vmin, vmax], [vmin, vmax], 'k--', lw=1)
axes[1].set_xlabel('Injected velocity (km/s)')
axes[1].set_ylabel('Recovered velocity (km/s)')
axes[1].set_title('Injected vs recovered')

plt.tight_layout()
plt.show()