Skip to content

Tweak echelle code to support different parameterization #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Yashvi-Sharma wants to merge 8 commits into
base: main
Choose a base branch
from
Draft

Tweak echelle code to support different parameterization #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
2365eb7
Tweak echelle code to support different parameterization
baileyji Feb 6, 2026
f3abb8a
ZSHTR-27 analytical echelle trace calculation, updated spectraltracel…
Yashvi-Sharma Feb 24, 2026
8254953
Stub function for analytical spectral efficiency effect, updated ZSho…
Yashvi-Sharma Mar 5, 2026
94a8367
Add spectrograph factory and enhance echelle parameter handling. add …
baileyji Mar 6, 2026
9b80cc3
Fix logging format for selector wheel warning message.
baileyji Mar 6, 2026
d5a107f
Implement initial computation of echelle efficiencies from analytical…
baileyji Mar 6, 2026
ba9ed12
Fix incorrect dictionary method usage in spectral trace initialization
baileyji Mar 6, 2026
45a699c
Refactor spectral trace initialization and parameter handling for ech…
baileyji Mar 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion scopesim/effects/selector_wheel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,8 @@ def apply_to(self, obj, **kwargs):
def get_effect(self, selector_value):
eff = None
if (selector_value is None) or (selector_value not in self.wheel_effects.keys()):
logger.warning(f"Either None or Missing value of {self.meta["selector_key"]} requested: {selector_value}, returning None.")
logger.warning(f"Either None or Missing value of {self.meta['selector_key']} "
f"requested: {selector_value}, returning None.")
else:
eff = self.wheel_effects[selector_value]
return eff
Expand Down
93 changes: 91 additions & 2 deletions scopesim/effects/spectral_efficiency.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# -*- coding: utf-8 -*-
"""Spectral grating efficiencies."""

from typing import ClassVar
from typing import ClassVar, Callable

import numpy as np
from astropy.io import fits
Expand All @@ -12,7 +12,9 @@
from .effects import Effect
from .ter_curves import TERCurve
from .ter_curves_utils import apply_throughput_to_cube
from ..utils import figure_factory, get_logger
from ..utils import figure_factory, get_logger, check_keys
from .data_container import DataContainer
from ..optics import echelle


logger = get_logger(__name__)
Expand Down Expand Up @@ -127,3 +129,90 @@ def plot(self):
axes.legend()

return fig


class EchelleSpectralEfficiency(Effect):
"""
Spectral efficiency list from analytical calculations of the blaze function for ZShooter gratings.
Requires same input trace parameter table as EchelleSpectralTraceList, supply as kwarg "filename"
"""
z_order: ClassVar[tuple[int, ...]] = (630,)

def __init__(self, **kwargs):
super().__init__(**kwargs)
self.efficiency_generator = self._generate_efficiency_curve_func(DataContainer(filename=kwargs.pop('filename')))
Copy link
Collaborator Author

@Yashvi-Sharma Yashvi-Sharma Mar 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need to call DataContainer here, the super().init already handles the loading of files specified under 'filename' kwarg. So the trace_params table can be accessed with self.table

self.efficiencies = {}

def _generate_efficiency_curve_func(self, trace_params) -> Callable:
spectrographs = {}
for row in trace_params.table:
prefix = row["prefix"] # note trance ids are assumed to be prefix_{order}
min_order = row['m0'] - row['n']
max_order = row['m0']
min_wave = row['min_wave'] * u.Unit(trace_params.meta["min_wave_unit"])
max_wave = row['max_wave'] * u.Unit(trace_params.meta["max_wave_unit"])
design_res = row['design_res']
focal_len = row['focal_length'] * u.Unit(trace_params.meta["focal_length_unit"])
disp_npix = row['n_disp'] - 2 * row['detector_pad']
xdisp_npix = row['n_xdisp']- 2 * row['detector_pad']
pix_size = row['pixel_size'] * u.Unit(trace_params.meta["pixel_size_unit"])
echelle_angle = np.deg2rad(row['echelle_blaze'])*u.rad
xdisp_beta_center = np.deg2rad(row['xbeta_center'])*u.rad

xdisp_groove_length = u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq']
echelle_groove_length = u.Unit(trace_params.meta["disp_freq_unit"]) / row['disp_freq']
pix_per_res_elem = row['fwhm']

spectrograph = echelle.spectrograph_factory(min_wave, max_wave, focal_len,
design_res, echelle_angle, min_order, max_order,
echelle_groove_length, pix_per_res_elem, disp_npix, xdisp_npix,
pix_size, xdisp_groove_length=xdisp_groove_length,
xdisp_beta_center=xdisp_beta_center)

spectrographs[prefix] = spectrograph

def efficiency_curve(trace_id, wavelength):
"""Trace ID MUST be in the form prefix_{order}"""
prefix, _, order = trace_id.partition('_')
order = int(order)
spec = spectrograph[prefix]
# blaze = spec.blaze(wavelength)[trace_id] # compute all of them and then subscript
blaze = spec.grating.blaze(spec.grating.beta(wavelength, order), order)
xdisp = spec.xdisp_efficiency(wavelength)
return blaze*xdisp

return efficiency_curve

def apply_to(self, obj, **kwargs):
"""Interface between FieldOfView and SpectralEfficiency."""
trace_id = obj.trace_id

swcs = WCS(obj.hdu.header).spectral
with u.set_enabled_equivalencies(u.spectral()):
wave = swcs.pixel_to_world(np.arange(swcs.pixel_shape[0])) << u.um

try:
efficiency = self.efficiency_generator(trace_id, wave)
params = {"description": trace_id}
params.update(self.meta)
effic_curve = TERCurve(array_dict={"wavelength": wave, "transmission": efficiency}, **params)
self.efficiencies[trace_id] = effic_curve
except:
raise ValueError(f"Error generating efficiency curve for trace {trace_id} with wavelength range {wave.min()} - {wave.max()}")

obj.hdu = apply_throughput_to_cube(obj.hdu, effic_curve.throughput, wave)
return obj

def plot(self):
"""Plot the grating efficiencies."""
fig, axes = figure_factory()
for name, effic in self.efficiencies.items():
wave = effic.throughput.waveset
axes.plot(wave.to(u.um), effic.throughput(wave), label=name)

axes.set_xlabel("Wavelength [um]")
axes.set_ylabel("Grating efficiency")
axes.set_title(f"Grating efficiencies {self.display_name}")
axes.legend()

return fig
53 changes: 22 additions & 31 deletions scopesim/effects/spectral_trace_list.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,7 +137,6 @@ def __init__(self, **kwargs):

if self._file is not None:
self.make_spectral_traces()

self.update_meta()

def make_spectral_traces(self):
Expand Down Expand Up @@ -524,7 +523,7 @@ class EchelleSpectralTraceList(SpectralTraceList):
instead of loading them from FITS file. The arguments required to define the echelle traces are supplied through
a txt file containing a table of parameters using the filename kwarg.

Below is an example of how to define the echelle trace parameters (see irdb/ZShooter/traces/echelle_trace_parameters.txt):
Below is an example of how to define the echelle trace parameters (see irdb/ZShooter_v1/traces/echelle_trace_parameters.txt):
----------------------------------------------------------------
# min_wave_unit : nm
# max_wave_unit : nm
Expand All @@ -539,10 +538,10 @@ class EchelleSpectralTraceList(SpectralTraceList):
# xdisp_freq_unit : mm
# slitwidth_unit : arcsec

prefix aperture_id image_plane_id m0 n min_wave max_wave echelle_blaze focal_length fwhm detector_pad pixel_size n_disp n_xdisp disp_freq xdisp_freq slitwidth dispdir plate_scale
nIR 0 0 40 24 970 2500 64.2 225 4.7 10 0.015 4096 4096 45 175 10 x 0.159574468085
gri 1 1 36 18 490 1020 64.2 225 4.7 10 0.015 4096 4096 100 500 10 x 0.159574468085
ub 2 2 29 11 315 515 64.2 225 4.7 10 0.015 4096 4096 200 1000 10 x 0.159574468085
prefix aperture_id image_plane_id m0 n min_wave max_wave design_res echelle_blaze focal_length fwhm detector_pad pixel_size n_disp n_xdisp disp_freq xdisp_freq slitwidth dispdir plate_scale xbeta_center
ub 2 2 29 11 315 515 20000 64.2 225 4.7 10 0.015 4096 4096 200 1000 10 x 0.159574468085
nIR 0 0 40 24 970 2500 20000 64.2 225 4.7 10 0.015 4096 4096 45 175 10 x 0.159574468085
gri 1 1 36 18 490 1020 20000 64.2 225 4.7 10 0.015 4096 4096 100 500 10 x 0.159574468085
----------------------------------------------------------------

The calculated traces are stored in the same HDUList format as required by SpectralTraceList,
Expand All @@ -555,7 +554,7 @@ class EchelleSpectralTraceList(SpectralTraceList):
def __init__(self, **kwargs):
check_keys(kwargs, self.required_keys, action="error")

trace_params = DataContainer(filename=kwargs['filename'])
trace_params = DataContainer(filename=kwargs.pop('filename'))
hdulist = self._generate_trace_hdulist(trace_params)
kwargs["hdulist"] = hdulist
super().__init__(**kwargs)
Expand Down Expand Up @@ -588,35 +587,27 @@ def _generate_trace_hdulist(self, trace_params):
max_order = row['m0']
min_wave = row['min_wave'] * u.Unit(trace_params.meta["min_wave_unit"])
max_wave = row['max_wave'] * u.Unit(trace_params.meta["max_wave_unit"])
design_res = row['design_res']
focal_len = row['focal_length'] * u.Unit(trace_params.meta["focal_length_unit"])
xdisp_npix = row['n_xdisp']
disp_npix = row['n_disp'] - 2 * row['detector_pad']
xdisp_npix = row['n_xdisp'] - 2 * row['detector_pad']
pix_size = row['pixel_size'] * u.Unit(trace_params.meta["pixel_size_unit"])
x_disp_len = (xdisp_npix - 2 * row['detector_pad']) * pix_size
echelle_angle = np.deg2rad(row['echelle_blaze'])
alpha = np.deg2rad(row['alpha'])
beta_center = np.deg2rad(row['beta_center'])
# cross_disperser = echelle.GratingSetup(
# groove_length=u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq'],
# guess_littrow=(min_wave, max_wave,
# x_disp_len, focal_len))
cross_disperser = echelle.GratingSetup(alpha=alpha, beta_center=beta_center,
delta=beta_center,
groove_length=u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq'])

ss = echelle.SpectrographSetup((min_order, max_order),
max_wave,
row['fwhm'] * u.Unit(trace_params.meta["fwhm_unit"]),
focal_len,
echelle.GratingSetup(alpha=echelle_angle, beta_center=echelle_angle,
delta=echelle_angle,
groove_length=u.Unit(trace_params.meta["disp_freq_unit"]) / row['disp_freq']),
echelle.Detector(row['n_disp'], xdisp_npix, pix_size),
cross_disperser=cross_disperser
)
echelle_angle = np.deg2rad(row['echelle_blaze'])*u.rad
xdisp_beta_center = np.deg2rad(row['xbeta_center'])*u.rad

xdisp_groove_length = u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq']
echelle_groove_length = u.Unit(trace_params.meta["disp_freq_unit"]) / row['disp_freq']
pix_per_res_elem = row['fwhm']

ss = echelle.spectrograph_factory(min_wave, max_wave, focal_len,
design_res, echelle_angle, min_order, max_order,
echelle_groove_length, pix_per_res_elem, disp_npix, xdisp_npix,
pix_size, xdisp_groove_length=xdisp_groove_length,
xdisp_beta_center=xdisp_beta_center)

fsr_edges = ss.edge_wave(fsr=True)

slit_edge = (row['slitwidth'] / 2) * u.Unit(trace_params.meta["slitwidth_unit"])
slit_edge = (row['slitlength'] / 2) * u.Unit(trace_params.meta["slitlength_unit"])
slit_pos = np.linspace(-slit_edge, slit_edge, num=3)
slit_offset_pix = slit_pos / (row['plate_scale'] * u.arcsec)

Expand Down
Loading