End-to-end Kilonova Simulation

Kilonovae are high-energy astrophysical transients formed by the merger of a two neutron stars. They are short-lived (a few weeks) and evolve rapidly (i.e. drop in brightness), with the SED quickly shifting to the redder wavelengths in a matter of days after the merger. Which is why it is important to follow-up kilonovae ASAP, and get spectra across visible and NIR bands.

ZShooter will be uniquely suited with its rapid ToO capability, to observe KNe early and late in the sequence.

Imports

# basic
from copy import deepcopy
from pathlib import Path
from typing import Optional, Literal
import logging
import importlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import astropy.units as u
from astropy.io import fits
from astropy.table import Table
from scipy.ndimage import median_filter

# zs_scopesim_tools convenience routines
import zs_scopesim_tools.helpers
from zs_scopesim_tools.helpers import (disable_scopesim_top_level_catch, warning_prevent_sync_alt_ra_dec,
                                       ignore_warnings, set_scopesim_log_level, make_local_output_dir, configure_irdb_path, timed, list_effects, simulate, save_zshooter_readout, subtract_hduls)
# scopesim and scopesim source creation related
import scopesim as sim
import scopesim_templates as sim_tp
from scopesim.source.source import Source
from scopesim.source import source_templates
from zs_scopesim_tools.sources import lamp_flat, transient, kilonova_spectra
from zs_scopesim_tools.plots import PlotScene, plot_readout_overview
from synphot.units import convert_flux, FLAM

# DRP
from zsdrp.ZSHOOTER import ZSHOOTER
import zsdrp.pipe, zsdrp.steps


configure_irdb_path()  # REQUIRED the first time after installation, sets local IRDB directory path in ScopeSim
set_scopesim_log_level("warning")  # set ScopeSim log toggles
disable_scopesim_top_level_catch()
ignore_warnings()
logging.disable(logging.WARNING)

Setup

Set up OUTPUT directory, remember to change it if running different simulations from same notebook to prevent accidental overrides.

OUTDIR = make_local_output_dir(dirname="01_Kilonovae/")

Load the latest ZShooter instrument model configurations (from local IRDB/ZShooter_v2).

zs_config = sim.UserCommands(use_instrument="ZShooter_v2", set_modes=["SPEC"])

list_effects() is a helper to make a dataframe of all the optical effects in the model with some basic information about each.

df = list_effects(zs_config)
df.head()

	included	effect	class	config	selector	data file
0	True	atmo_transmission	TAPASTERCurve	MaunaKea (ATMO)		TER_maunakea_tapas.dat
1	True	continuum_emission	SkyBackgroundTERCurve	MaunaKea (ATMO)
2	True	airglow_and_interline_continuum	PalaceAirglowEmission	MaunaKea (ATMO)
3	True	seeing_psf	AOEnhanceablePSF	MaunaKea (ATMO)
4	True	adc_residuals	ADCShift	MaunaKea (ATMO)		adc_residuals.dat

Helper function to update observing parameters and run simulation in one-shot

Example usage:

hdul, train = update_params_and_observe(
                        object_name = "object_name_identifier", # string name
                        imagetype = "SCIENCE", # one of BIAS, DOMEFLAT, ARCS, SCIENCE
                        source = Source, # scopesim Source object instance
                        exptimes = [300, 300, 300, 100, 100, 100], # list of 6 numbers, one per channel
                        ndit = 2, # Number of exposures to take at above exptime, same for all 6 channels
                        output_dir = "/path/to/output/directory",
                        disable_effects = ["name_of_effect_to_disable", "..."], # see list_effects()
                        disable_backgrounds = ["sky_lines"], # from "sky", "sky_lines", "sky_continuum"
                        airmass = 1.3,
                        seeing = 0.7,
                        brightness = "grey", # options: "bright", "grey", "dark"
                        vis_slit = 0.7, # options: 1.25, 0.70, 0.57, 0.50, 0.42, 0.37, 0.33
                        nir_slit = 0.7, # options: 1.25, 0.70, 0.57, 0.50, 0.42, 0.37, 0.33
                        pixel_scale = None, # uses instrument default
                        pupil_angle = None, # uses instrument default
                        ao_on = False, # AO effect ON/OFF
                        diffuse_emissivity_on = False
)

Returns a tuple of:

• list of size ndit, each item of which is a list of 6 fits.HDUList (one per channel), the raw images from scopesim readout

• The scopesim OpticalTrain object

@timed
def update_params_and_observe(object_name: str,
                              imagetype: str,
                              source: Source,
                              exptimes: Optional[list | float] = None,
                              ndit: int = 1,
                              output_dir: Optional[str | Path]=None,
                              *,
                              disable_effects: Optional[list[str]]=None,
                              disable_backgrounds: Optional[list[Literal["sky", "sky_continuum", "sky_lines"]]]=None,
                              airmass: Optional[float]=None,
                              seeing: Optional[float]=None,
                              brightness: Optional[Literal["bright", "grey", "dark"]]=None,
                              vis_slit: Optional[float]=None,
                              nir_slit: Optional[float]=None,
                              pixel_scale: Optional[float]=None,
                              pupil_angle: Optional[float]=None,
                              ao_on: bool=False,
                              diffuse_emissivity_on: bool=False,
                              params: Optional[dict]=None
                              ) -> tuple[list, sim.optics.OpticalTrain]:
    """
    :param object_name: Name of the Source object
    :param imagetype: IMAGETYPE keyword (BIAS, DOMEFLAT, ARCS, SKY, SCIENCE)
    :param source: Source instance
    :param exptimes: List of exposure times in channel order (blue, green, red, yj, h, k) in seconds, or single value.
                If None, it is set to 0.0.
    :param ndit: number of exposures
    :param output_dir: Directory to save observed data in, FITS files are saved with {object_name}_{CHANNEL}.fits names
                        where CHANNEL is in [BLUE, GREEN, RED, YJ, H, K] (upper case)
    :param disable_effects: List of effect names to disable
    :param disable_backgrounds: List of background names to disable, options are "sky", "sky_continuum", or "sky_lines".
    :param airmass: Airmass
    :param seeing: Seeing (natural)
    :param brightness: Brightness of sky ("bright", "grey", or "dark")
    :param vis_slit: Slit width for Visible channel
    :param nir_slit: Slit width for NIR channel
    :param pixel_scale: Instrument pixel scale
    :param pupil_angle: Instrument pupil angle
    :param ao_on: True for AO correction on, defaults to False
    :param diffuse_emissivity_on: True for diffuse emissivity effect on, defaults to False
    :param params: Dict of any other UserCommands alias keys and values to set

    :return: tuple(HDUList, OpticalTrain)
    """
    if imagetype not in ["BIAS", "DOMEFLAT", "ARCS", "SKY", "SCIENCE"]:
        raise ValueError(f"Imagetype {imagetype} is incorrect! Choose from BIAS, DOMEFLAT, ARCS, SCIENCE.")
    # load ZShooter_v2 irdb model config
    config = sim.UserCommands(use_instrument="ZShooter_v2", set_modes=["SPEC"])
    # update with passed params
    config["!OBS.airmass"] = airmass if airmass else config["!OBS.airmass"]
    config["!OBS.seeing"] = seeing if seeing else config["!OBS.seeing"]
    config["!OBS.brightness"] = brightness if brightness else config["!OBS.brightness"]
    config["!OBS.pupil_angle"] = pupil_angle if pupil_angle else config["!OBS.pupil_angle"]
    config["!INST.vis_curr_slit"] = vis_slit if vis_slit else config["!INST.vis_curr_slit"]
    config["!INST.nir_curr_slit"] = nir_slit if nir_slit else config["!INST.nir_curr_slit"]
    config["!INST.pixel_scale"] = pixel_scale if pixel_scale else config["!INST.pixel_scale"]
    config["!INST.include_diffuse_emissivity"] = diffuse_emissivity_on
    config["!ATMO.ao_enabled"] = ao_on

    dits = [0.]*6 if exptimes is None else [exptimes]*6 if not isinstance(exptimes, list) else exptimes
    for dit, expkey in zip(dits, ['blue', 'green', 'red', 'yj', 'h', 'k']):
        config[f"!OBS.dit_{expkey}"] = dit
        config[f"!OBS.ndit_{expkey}"] = 1.

    if params is not None and isinstance(params, dict):
        params = {key:val for key, val in params.items() if key in config}
        # USE THIS TO UPDATE THE CONFIG WITH NEW PROPERTIES, IT IS NOT A NORMAL DICT
        config.update_alias(config, params)
    # make sure airmass, ra, dec are in sync if all are provided
    warning_prevent_sync_alt_ra_dec(config)

    # call simulate(), it inits optical train, runs observe() and readout()
    print(f'Simulating {object_name}')
    n_hduls = []
    for i in range(ndit):
        zs, hdul = simulate(config, source=source, disable_effects=disable_effects, disable_backgrounds=disable_backgrounds, hide_progress_bars=True)
        # check desired counts
        print(f'Max counts per channel: {[int(np.max(hdu[1].data)) for hdu in hdul]}')
        # save the simulated readout to fits files
        if output_dir is not None:
            _ = save_zshooter_readout(hdul, output_dir=output_dir, object_name=object_name,
                                          imagetype=imagetype.upper(), ndit=i, cmds=config, train=zs)
            print(f"Saved readout #{i}/{ndit} in {output_dir}")
        n_hduls.append(hdul)

    return n_hduls, zs

Create sources

Helper function to simulate kilonova on top of a galaxy.

def make_default_kilonova(
                    z: float = 0.0098, # redshift, default is for GW170817
                    host_mag: float = 14, # host galaxy surface brightness in AB mag, default is 14 for GW170817 host
                    pixel_scale: float = 0.01,
                    model_files: list[str] = ["knova_d1_n10_m0.025_vk0.30_Xlan1e-4.0.h5",
                                              "knova_d1_n10_m0.040_vk0.15_Xlan1e-1.5.h5"], # combo best fit to GW170817
                    phases: list[float] = [1.5, 7.5], # phases in kilonova evolution to create the sources for
                    plot: bool = True # plot each created source
):
    ...

List the desired kilonova model spectra filenames in the argument model_files and they will be pulled automatically from this repo (https://github.com/dnkasen/Kasen_Kilonova_Models_2017/tree/master/kilonova_models) to generate the source. The currently listed names are of the models that produced the combined best fit to GW170817 KN (Kasen 2017, https://www.nature.com/articles/nature24453)

List all the phases (in days) at which you want the model spectra in phases.

The output of the helper function is a dictionary of Source objects, keyed by the phase.

def make_default_kilonova(z: float = 0.0098,
                          host_mag: float = 14,
                          pixel_scale: float = 0.01,
                          model_files: list[str] = ["knova_d1_n10_m0.025_vk0.30_Xlan1e-4.0.h5",
                                                    "knova_d1_n10_m0.040_vk0.15_Xlan1e-1.5.h5"],
                          phases: list[float] = [1.5, 7.5, 20.0],
                          plot: bool = True):
    """
    :param z - redshift
    :param host_mag - apparent magnitude of the host galaxy
    :param pixel_scale - pixel scale of the image
    :param model_files - list of file names corresponding to the model files
    :param phases - list of phases of spectra needed
    :param plot - Plot the source
    """
    # base galaxy centered at 0, 0 arcsec
    gal = sim_tp.galaxy(sed='brown/NGC3870', z=z, amplitude=host_mag, pixel_scale=pixel_scale,
                        r_eff=2.5, n=1., ellip=0.3, theta=90, extend=5)
    # shift upwards by 2"
    gal.shift(dy=10.5 * u.arcsec)
    # kilonova spectra from external Kasen spectrum models
    kn_spectra = kilonova_spectra(model_files, phases=phases, redshift=z, combine_models=True)
    # init source dictionary
    kns = dict()
    for phase in phases:
        source = gal + transient(x=0, y=0, spextrum_template=kn_spectra[phase], name=f'KN {phase}d')
        kns[phase] = source
        if plot:
            _ = PlotScene().from_source(source, pixel_scale=pixel_scale)

    return kns

Make the kilonova sources at 1.5 d, 7.5 d and 20 d from explosion, and simulate!

KNs = make_default_kilonova(plot=True)

Simulating cal+science at one observing configuration

Set your observing parameters in cell below.

If you change any instrument related parameters for a science object simulation (like vis_slit, nir_slit, ao_on or diffuse_emissivity_on)

• Make sure to simulate a complementary set of cals with the same parameters!

• Make sure you CHANGE the output_dir so that your data doesn’t get overwritten!

airmass = 1.3
seeing = 0.7
brightness = "dark"
vis_slit = 0.7
nir_slit = 0.7
ao_on = False
diffuse_emissivity_on = False
params = {"!SIM.spectral.spectral_bin_width": 1e-5} # in um, (== 0.1 Angstrom, use 0.05 Angstrom for smaller slit widths)

args = {'airmass':airmass, 'seeing':seeing, 'brightness':brightness, 'vis_slit':vis_slit, 'nir_slit':nir_slit,
        'ao_on': ao_on, 'diffuse_emissivity_on': diffuse_emissivity_on, 'params': params}

# creating a sub-directory in OUTDIR for this run config, naming it based on config params for uniqueness
output_dir = OUTDIR / f"kilonova_X-{airmass}_NS-{seeing}_{brightness}_slit-{vis_slit}_ao-{int(ao_on)}"
output_dir.mkdir(parents=True, exist_ok=True)

overwrite_cals = True  # set True to redo cals
plot = True  # Set True for plotting

Cell below simulates bias, domeflat, arcs and a standard star

# for cals
disable_effects_cals = ["atmo_transmission", "continuum_emission", "airglow_and_interline_continuum",
                        "seeing_psf", "adc_residuals"]

if overwrite_cals or len(sorted(output_dir.glob("*bias*")))!=6:
    # bias
    bias_hduls, _ = update_params_and_observe(object_name="bias", imagetype="BIAS", source=None, exptimes=0,
                                             output_dir=output_dir, disable_effects=disable_effects_cals, **args)
    # arcs
    LAMP = lamp_flat()
    arcs_hduls, _ = update_params_and_observe(object_name="arcs", imagetype="ARCS", source=LAMP,
                                             exptimes=[40, 40, 3, 3, 20, 10], output_dir=output_dir,
                                             disable_effects=disable_effects_cals, **args)
    # domeflat
    DOMEFLAT = sim_tp.calibration.flat_field(temperature=3000*u.K, amplitude=8*u.ABmag, filter_curve='V', extend=60)
    domeflat_hduls, _ = update_params_and_observe(object_name="domeflat", imagetype="DOMEFLAT", source=DOMEFLAT,
                                                 exptimes=[120, 10, 3, 3, 3, 3], output_dir=output_dir,
                                                 disable_effects=disable_effects_cals, **args)
    # standard star, sky disabled to save time doing one less sim for nodding
    STANDARD = sim_tp.misc.point_source(sed=source_templates.ab_spectrum(mag=16))
    std_hduls, _ = update_params_and_observe(object_name="standard", imagetype="SCIENCE", source=STANDARD,
                                            exptimes=100, disable_backgrounds=["sky"], output_dir=output_dir, **args)
    if plot:
        _, _ = plot_readout_overview(bias_hduls[0],
                                        titles=[f"bias {chan}" for chan in ["BLUE", "GREEN", "RED", "YJ", "H", "K"]])
        _, _ = PlotScene().from_source(LAMP)
        _, _ = plot_readout_overview(arcs_hduls[0],
                                        titles=[f"arcs {chan}" for chan in ["BLUE", "GREEN", "RED", "YJ", "H", "K"]])
        _, _ = PlotScene().from_source(DOMEFLAT)
        _, _ = plot_readout_overview(domeflat_hduls[0],
                                        titles=[f"domeflat {chan}" for chan in ["BLUE", "GREEN", "RED", "YJ", "H", "K"]])
        _, _ = PlotScene().from_source(STANDARD)
        _, _ = plot_readout_overview(std_hduls[0],
                                        titles=[f"standard {chan}" for chan in ["BLUE", "GREEN", "RED", "YJ", "H", "K"]])
# for the precious RAM
del bias_hduls, domeflat_hduls, arcs_hduls, std_hduls

Simulating bias
Max counts per channel: [1040, 1040, 1040, 2, 2, 2]
Saved readout #0/1 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
update_params_and_observe(bias): 28.2 s
Simulating arcs
Max counts per channel: [24788, 33360, 33281, 32000, 35132, 37771]
Saved readout #0/1 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
update_params_and_observe(arcs): 29.0 s
Simulating domeflat
Max counts per channel: [39807, 41802, 49388, 68750, 76798, 72284]
Saved readout #0/1 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
update_params_and_observe(domeflat): 32.6 s
Simulating standard
Max counts per channel: [1138, 1173, 1203, 184, 198, 859]
Saved readout #0/1 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
update_params_and_observe(standard): 48.8 s

Observe each KN, and one sky exposure for each unique frame exposure time.

importlib.reload(zs_scopesim_tools.helpers)
from zs_scopesim_tools.helpers import subtract_hduls, save_zshooter_readout

sky_hduls, _ = update_params_and_observe(object_name=f"sky", imagetype="SKY", source=None,
                                     exptimes=1800, output_dir=output_dir, **args)

# 2 frames x 1800s = 3600s = 1hr
for phase, source in KNs.items():
    # with full sky background
    objname = f"KN_{phase}d"
    hduls, _ = update_params_and_observe(object_name=objname, imagetype="SCIENCE", source=source,
                                            exptimes=1800, ndit=2, output_dir=output_dir, **args)
    # subtract sky
    for i, hdul in enumerate(hduls):
        sub_hdul = subtract_hduls(hdul, sky_hduls[0])
        _ = save_zshooter_readout(sub_hdul, output_dir=output_dir, object_name=f"{objname}_SUB", imagetype="SCIENCE", ndit=i, extra_headers={"SUBTRACTED": True})

_ = plot_readout_overview(sky_hduls[0], titles=[f"Sky, {chan}" for chan in ["BLUE", "GREEN", "RED", "YJ", "H", "K"]])
_ = plot_readout_overview(hdul, titles=[f"{objname}, {chan}" for chan in ["BLUE", "GREEN", "RED", "YJ", "H", "K"]])
_ = plot_readout_overview(sub_hdul, titles=[f"{objname}-Sky, {chan}" for chan in ["BLUE", "GREEN", "RED", "YJ", "H", "K"]])

Simulating sky
Max counts per channel: [1073, 6234, 17905, 115365, 583790, 332603]
Saved readout #0/1 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
update_params_and_observe(sky): 58.1 s
Simulating KN_1.5d
Max counts per channel: [1079, 6245, 18243, 114992, 583922, 333386]
Saved readout #0/2 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
Max counts per channel: [1081, 6225, 18372, 115252, 584787, 332740]
Saved readout #1/2 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
update_params_and_observe(KN_1.5d): 115.2 s
Simulating KN_7.5d
Max counts per channel: [1068, 6148, 17928, 115038, 585589, 333439]
Saved readout #0/2 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
Max counts per channel: [1070, 6210, 17982, 115352, 584312, 333295]
Saved readout #1/2 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
update_params_and_observe(KN_7.5d): 108.6 s
Simulating KN_20.0d
Max counts per channel: [1070, 6190, 17941, 114954, 584203, 333314]
Saved readout #0/2 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
Max counts per channel: [1069, 6182, 17875, 114887, 584024, 333977]
Saved readout #1/2 in /Users/yashvi/Desktop/ZShooter/zs-scopesim/notebooks/01_Kilonovae/kilonova_X-1.3_NS-0.7_dark_slit-0.7_ao-0
update_params_and_observe(KN_20.0d): 122.6 s

Reducing the simulated readouts with Pyreduce

Set REDUCED_DIR - where outputs of reduction will be saved (simply set to reduced/ folder inside simulation output) directory.

REDUCED_DIR = output_dir / "reduced"
REDUCED_DIR.mkdir(exist_ok=True)

Since pyreduce reduction steps are designed to run per channel, load the raw file paths for each channel and object and store them in a dictionary for easy access using in-built method of pyreduce’s Instrument subclass - ZSHOOTER

from zsdrp.ZSHOOTER import ZSHOOTER

zspy = ZSHOOTER()
files = {}
for chan in zspy.config.channels:
    stepnames = ['bias','flat','wavecal_master','science']
    file_groups = zspy.sort_files(str(output_dir), target='', night='', channel=chan, steps=stepnames)
    _, fls = file_groups[0]
    files[chan] = {name: fls[name] for name in stepnames}
# files[chan]

100%|██████████| 102/102 [00:00<00:00, 1148.22it/s]
100%|██████████| 102/102 [00:00<00:00, 1210.48it/s]
100%|██████████| 102/102 [00:00<00:00, 1383.60it/s]
100%|██████████| 102/102 [00:00<00:00, 1540.57it/s]
100%|██████████| 102/102 [00:00<00:00, 1517.31it/s]
100%|██████████| 102/102 [00:00<00:00, 1250.84it/s]

Commented cell below shows how to construct files nested dict of filenames manually, if needed.

# files = {}
# for chan in ['BLUE', 'GREEN', 'RED', 'YJ', 'H', 'K']:
#     channel_files = dict()
#     channel_files["flat"] = sorted(output_dir.glob(f"*domeflat_{chan}*")) # pattern match filenames with 'domeflat' in them inside OUTDIR folder
#     channel_files["wavecal_master"] = sorted(output_dir.glob(f"*arcs_{chan}*"))
#     channel_files["bias"] = sorted(output_dir.glob(f"*bias_{chan}*"))
#
#     # similarly as above, pattern match science simulation files created earlier in the notebook
#     channel_files["science"] = []
#     for obj in ["KN_1.5d", "KN_7.5d"]:
#         channel_files["science"] += sorted(output_dir.glob(f"*{obj}_{chan}*"))
#     print(f"Found {len(channel_files['flat'])} domeflats, {len(channel_files['wavecal_master'])} arcs, {len(channel_files['bias'])} bias and {len(channel_files['science'])} science images for channel {chan}.")
#     print('--------------------------------')

Reduce

The output of zsdrp.pipe.run_reduction() is a dictionary of per channel calibration products and intermediate spectra (up to blaze normalization). Flux correction is applied after.

import zsdrp.pipe, zsdrp.steps
importlib.reload(zsdrp.pipe)
importlib.reload(zsdrp.steps)

results = zsdrp.pipe.run_reduction(files, channels=None, output_dir=REDUCED_DIR, disable_steps=['mask'], print_params=False, disable_all_plots=True, reuse_cals=True)

Running reduction for channel BLUE...
Running reduction for channel GREEN...
Running reduction for channel RED...
Running reduction for channel YJ...
Running reduction for channel H...
Running reduction for channel K...

# import zsdrp.utils
# importlib.reload(zsdrp.utils)
# from zsdrp.utils import plot_spectra_object
# for objname, spectra in results['RED']['spectra'].items():
#     _ = plot_spectra_object(spectra, title=f'{objname}-{'RED'}', xlabel='x-pixel', ylabel='Scaled counts + offset')
#     plt.show()

Apply flux normalization. The standard source simulated above is a simple ab_spectrum at 16th mag so that’s our ground truth, i.e. reference.

from scopesim.source.source_templates import ab_spectrum
from synphot.units import convert_flux

reference = ab_spectrum(mag=16)
reference_wave = reference.waveset
reference_flux = convert_flux(reference_wave, reference(reference_wave), "FLAM")

final_spectra = zsdrp.pipe.run_flux_calibration_and_stitch(results, reference_wave=reference_wave, reference_flux=reference_flux, plot=False, output_dir=str(REDUCED_DIR))

# for the precious RAM
del results
del final_spectra

final_spectra is a dictionary of stitched and flux normalized spectra keyed by object name (science objects only).

Plot

fig, ax = plt.subplots(2,1, figsize=(20,8))
ax[0].grid(True)
ax[1].grid(True)
ax[0].tick_params(axis='y', labelsize=12)
ax[1].tick_params(axis='y', labelsize=12)
ax[1].tick_params(axis='x', labelsize=16)
ax[0].set_ylabel('Flux (erg/s/cm$^2/\AA$)', fontsize=18)
ax[1].set_ylabel('SNR', fontsize=18)
ax[1].set_xlabel('Wavelength ($\AA$)', fontsize=18)

for i,obj in enumerate(["KN_1.5D", "KN_7.5D", "KN_20.0D"]):
    objname = f"{obj}_SUB" # set pattern string prefix
    flname = sorted(REDUCED_DIR.glob(f'*{objname.upper()}*stitched.txt'))[0] # set pattern string ACCURATELY here
    with open(flname, 'r', encoding='utf-8') as f:
        tbl = pd.read_csv(f, sep='\t', comment='# ', engine='python')
        wave = tbl['wavelength']
        spec = tbl['flux']
        sig = tbl['eflux']
        snr = spec/sig
        good = np.isfinite(snr)
        print(f"{objname}: mean SNR = {np.mean(snr[good]):.2f}, "
              f"median SNR = {np.median(snr[good]):.2f}")

        ax[0].plot(wave, spec, alpha=0.2, color=f"C{i}")
        ax[0].plot(wave, median_filter(spec, 51), label=obj, color=f"C{i}")
        ax[1].plot(wave, snr, alpha=0.2, color=f"C{i}")
        ax[1].plot(wave, median_filter(snr , 51), color=f"C{i}")


ax[0].set_title('Exptime = 3600 s', fontsize=24)
ax[0].set_xlim(3070, 25000)
ax[1].set_xlim(3070, 25000)
ax[0].set_ylim(-0.01e-15, 0.2e-15)
ax[1].set_ylim(-10, 100)
ax[0].legend(ncol=5)
plt.show()

KN_1.5D_SUB: mean SNR = 37.18, median SNR = 34.54
KN_7.5D_SUB: mean SNR = 12.60, median SNR = 7.71
KN_20.0D_SUB: mean SNR = 0.45, median SNR = 0.22

Checking stitching

from astropy.table import Table
from scipy.ndimage import median_filter

obj = "KN_1.5D_SUB"
specfiles = sorted((output_dir / "reduced").glob(f'*{obj.upper()}_*_fluxcal.txt'))
stitched = sorted((output_dir / f"reduced").glob(f'*{obj.upper()}_*stitched.txt'))[0]

fig, ax = plt.subplots(2,1, figsize=(20,8), tight_layout=True, sharex=True)
ax[0].grid(True)
ax[1].grid(True)
ax[0].tick_params(axis='y', labelsize=12)
ax[1].tick_params(axis='y', labelsize=12)
ax[1].tick_params(axis='x', labelsize=16)
ax[0].set_ylabel('Flux (erg/s/cm$^2/\AA$)', fontsize=18)
ax[1].set_ylabel('SNR', fontsize=18)
ax[1].set_xlabel('Wavelength ($\AA$)', fontsize=18)
ax[0].set_title(f"{obj.upper()} 3600s, slit {vis_slit}", fontsize=18)

colors = {'BLUE': 'cornflowerblue', 'GREEN': 'green', 'RED': 'orange', 'YJ': 'orangered', 'H': 'red', 'K': 'maroon'}
for fl in specfiles:
    chan = str(fl).split('/')[-1].split('_')[-2].upper()
    tbl = Table.read(str(fl), format='ascii', delimiter='\t')
    snr = tbl['flux']/tbl['eflux']
    ax[0].plot(tbl['wavelength'], tbl['flux'], color=colors[chan], alpha=0.3)
    ax[0].plot(tbl['wavelength'], median_filter(tbl['flux'], 51), label=chan, color=colors[chan])
    ax[1].plot(tbl['wavelength'], snr, alpha=0.3, color=colors[chan])
    ax[1].plot(tbl['wavelength'], median_filter(snr, 51), color=colors[chan])
stbl = Table.read(str(stitched), format='ascii', delimiter='\t')
ax[0].plot(stbl['wavelength'], median_filter(stbl['flux'], 51), color='gray', alpha=0.7, label='Stitched')
ax[1].plot(stbl['wavelength'], median_filter(stbl['flux']/stbl['eflux'], 51), color='white', alpha=0.7)
model = convert_flux(stbl['wavelength'], KNs[1.5].spectra[1](stbl['wavelength']), "FLAM")
ax[0].plot(stbl['wavelength'], model.value, color='white', label='Model', ls='--')

ax[0].set_xlim(3000, 25000)
ax[0].set_ylim(-0.1e-15, 0.2e-15)
ax[1].set_ylim(-10, 100)
ax[0].legend(ncol=7, fontsize=16)
plt.show()

model

$[9.9355329 \times 10^{-15},~9.9354824 \times 10^{-15},~9.935432 \times 10^{-15},~\dots,~6.0472525 \times 10^{-16},~6.0471365 \times 10^{-16},~6.0470205 \times 10^{-16}] \; \mathrm{flam}$