Kilonova Simulation

Kilonovae evolve rapidly, the SED moves hard to the red, and ZShooter will observe early and late in the sequence.

See e.g. Kasen et al. 2017 synthetic kilonova spectra, plus one observed sequence panel from GW170817-era follow-up.

Set up Simulator

1. UserCommands(use_instrument=..., set_modes=[...])
        └─ parse default.yaml + mode yamls, build rc config

import matplotlib.pyplot as plt
from astropy import units as u
import scopesim as sim
import scopesim_templates as sim_tp
import importlib.resources
import pathlib
import numpy as np

from spextra import Spextrum
import synphot
import astropy.units as u
from scopesim.source.source import Source

def plot_source(source):
    wave = np.linspace(0.3, 2.5, 1001) * u.um
    num_fields = len(source.fields)
    fig, axs = plt.subplots(figsize=(6,2*num_fields), nrows=num_fields, ncols=2, width_ratios=[1,2], constrained_layout=True)
    if num_fields == 1:
        axs = np.array([axs])
    for i, field in enumerate(source.fields):
        ax1, ax2 = axs[i]
        ax1.imshow(field.data, origin='lower', cmap='viridis')
        ax2.plot(wave, field.spectrum(wave))
        ax1.set_title('Spatial Profile')
        ax1.set_xlabel('X (arcsec)')
        ax1.set_ylabel('Y (arcsec)')
        ax2.set_title('Spectrum')
        ax2.set_xlabel('Wavelength (um)')
        ax2.set_ylabel(r'Flux (ph s$^{-1}$ cm$^{-2}$ $\AA^{-1}$)')

## configure IRDB path
try:
    irdb_path = str(pathlib.Path(importlib.resources.files('irdb')).parent)
except ModuleNotFoundError:
    raise ModuleNotFoundError("Please install irdb package or comment this error and set the line below to your irdb path")
    # e.g. "/Users/yashvi/Desktop/ZShooter/irdb/"
    irdb_path = ""

sim.rc.__config__["!SIM.file.local_packages_path"] = irdb_path
sim.rc.__config__["!SIM.file.search_path"].append(irdb_path)

## configure logging level
sim.utils.set_console_log_level("DEBUG")  ## INFO, DEBUG, WARNING, ERROR, CRITICAL

## user commands
cmd = sim.UserCommands(use_instrument="ZShooter_v2", set_modes=["SPEC"])

Set Settings

cmd['!OBS.sky_brightnesss']='bright'
cmd['!OBS.airmass']=2.0

Create ZShooter’s OpticalTrain

2. OpticalTrain(cmd)
        ├─ OpticsManager loads all yaml effects into optical elements
        ├─ FOVManager runs z_order 200–299 effects → FovVolumeList
        ├─ ImagePlane objects initialised from z_order 100–199 / 400–499
        └─ DetectorManager creates Detector objects from z_order 400–499 effects

zs = sim.OpticalTrain(cmd)

Set Other Settings

Create a scene

spectrum_file=''


# Start with a Galaxy off to the side
scene = sim_tp.extragalactic.galaxies.spiral_two_component(extent=200*u.arcsec, fluxes=(20, 20)*u.mag)  # loaded from saved fits template of NGC1232L
scene.shift(dy=-30*0.004)


# download spectrum from spextra library

sed = Spextrum('sne/sn1a').redshift(0.02) # redshift it to match the host galaxy
spec = sed.scale_to_magnitude(24*u.ABmag, 'g')


# add sn to scene 3" above center
scene.append(Source(spectra=spec, x=[0.0], y=[3.0], ref=[0.], weight=[1.0]))

# add user file to scene -3" above center
scene.append(Source(spectra=synphot.SourceSpectrumn.from_file(spectrum_file), x=[0.0], y=[-3.0], ref=[0.], weight=[1.0]))

plot_source(scene)
plt.figure()
spec.plot()
plt.figure()
plot_source(scene)

Observe It

3. observe(source)
        ├─ z_order 500–599: source-level spectral transmission applied
        ├─ For each FoV:
        │     ├─ Extract overlapping source data → attach to FoV
        │     └─ z_order 600–699: 3D FOV effects (PSF, trace, ADC, shifts…)
        ├─ Project FoV HDU → correct ImagePlane
        └─ z_order 700–799: 2D image plane effects

Note that some final settings can be set here. The locations of settings varies due to ….. Best practice is to …., think of them as settings of category ….

kwargs = {'exposure_time': 1000*u.ms, 'n_exposures': 1}  ## Final sewttings

zs.observe(scene, **kwargs)

Read It Out

4. readout()
        ├─ Map ImagePlane → Detector
        ├─ z_order 800–899: detector effects (noise, dark, readout)
        ├─ z_order 900–999: detector array effects
        └─ z_order 1000–1099: FITS header effects → return list of HDUs

import datetime
file_root = f"{sim.instrument}_{sim.mode}-{sim.scene_name}_{datetime.datetime.now().strftime('%Y%m%d_%H%M%S')}"
file = f"{file_root}.fits"

hdus = zs.readout(file)

from astropy.visualization import ZScaleInterval

# reduce space between subplots
fig, axes = plt.subplots(2,3, figsize=(12,8), gridspec_kw={'wspace':0.2}, tight_layout=True)
titles = ['B', 'G', 'R', 'YJ', 'H', 'K']
i = 0
for ax, imhdu in zip(axes.flat, hdus):
    interval = ZScaleInterval()
    vmin, vmax = interval.get_limits(imhdu[1].data)
    ax.axis("off")
    im = ax.imshow(imhdu[1].data, origin='lower', vmin=vmin, vmax=vmax, cmap='Grays_r')
    ax.set_title(titles[i])
    plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
    i += 1

plt.subplots_adjust(wspace=0.0, hspace=0.0)