From 2365eb7069642530a6a9e86b20741705ba38fa53 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Thu, 5 Feb 2026 16:35:04 -0800 Subject: [PATCH 1/8] Tweak echelle code to support different parameterization --- scopesim/optics/echelle.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/scopesim/optics/echelle.py b/scopesim/optics/echelle.py index bdb1100f..27a57831 100644 --- a/scopesim/optics/echelle.py +++ b/scopesim/optics/echelle.py @@ -136,7 +136,7 @@ class SpectrographSetup: def __init__( self, order_range: tuple, - final_wave: u.Quantity, + design_res: u.Quantity, pixels_per_res_elem: float, focal_length: u.Quantity, grating: GratingSetup, @@ -150,7 +150,7 @@ def __init__( (and that the slit image is gaussian) :param tuple order_range: order range of the spectrograph - :param u.Quantity final_wave: longest wavelength at the edge of detector + # :param u.Quantity final_wave: longest wavelength at the edge of detector :param float pixels_per_res_elem: number of pixels per resolution element of spectrometer :param u.Quantity focal_length: the focal length of the detector :param GratingSetup grating: configured grating @@ -162,7 +162,8 @@ def __init__( order_range = order_range[::-1] self.m0 = order_range[0] self.m_max = order_range[1] - self.l0 = final_wave + self.l0 = grating.wave(grating.beta_center, order_range[0])*(1/order_range[0]/2+1) + self.design_res = design_res self.grating = grating self.detector = detector self.focal_length = focal_length @@ -444,14 +445,14 @@ def angular_dispersion(self): beta = self.beta_for_x_pixel(np.arange(self.detector.n_pix_x) + .5) return self.grating.angular_dispersion(self.orders[:, None], beta) - @property - def design_res(self): - """ - :return: design resolution for spectrometer - Assumes that m0 (the longest wavelength order) FSR fills detector with some sampling - """ - dlambda = self.fsr(self.m0) / self.detector.n_pix_x * self.nominal_pixels_per_res_elem - return self.l0 / dlambda + # @property + # def design_res(self): + # """ + # :return: design resolution for spectrometer + # Assumes that m0 (the longest wavelength order) FSR fills detector with some sampling + # """ + # dlambda = self.fsr(self.m0) / self.detector.n_pix_x * self.nominal_pixels_per_res_elem + # return self.l0 / dlambda @property def average_res(self): @@ -489,7 +490,7 @@ def plot_echellogram(self, center_orders=True, title='', blaze=False, cross_disp oset = waves[1] if center_orders else 0 plt.plot(waves - oset, [i] * 3, '*', color=f'C{ii}') plt.plot(fsr_edges[ii] - oset, [i] * 2, '.', color=f'C{ii}', - label=f'$\lambda=${waves[1]:.0f}') + label=f'$\\lambda=${waves[1]:.0f}') plt.legend() plt.xlabel('Center relative wavelength (nm)') plt.ylabel('Order') From f3abb8ac98770ba7cc9cd0ceb520ef52a7538098 Mon Sep 17 00:00:00 2001 From: Yashvi-Sharma Date: Tue, 24 Feb 2026 11:10:08 -0800 Subject: [PATCH 2/8] ZSHTR-27 analytical echelle trace calculation, updated spectraltracelist effect to work with tweaks --- scopesim/effects/spectral_trace_list.py | 33 +++++++++++++------------ scopesim/optics/echelle.py | 2 +- 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index 816cb260..2a298de1 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -539,10 +539,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 + 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, @@ -588,6 +588,7 @@ 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'] pix_size = row['pixel_size'] * u.Unit(trace_params.meta["pixel_size_unit"]) @@ -599,18 +600,18 @@ def _generate_trace_hdulist(self, trace_params): # 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 = echelle.GratingSetup(alpha=alpha, beta_center=beta_center, delta=beta_center, + groove_length=u.Unit(trace_params.meta["xdisp_freq_unit"]) / row['xdisp_freq']) + + echelle_grating = 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']) + + ss = echelle.SpectrographSetup(order_range=(min_order, max_order), + design_res=design_res, + pixels_per_res_elem=row['fwhm'] * u.Unit(trace_params.meta["fwhm_unit"]), + focal_length=focal_len, + grating=echelle_grating, + detector=echelle.Detector(row['n_disp'], xdisp_npix, pix_size), cross_disperser=cross_disperser ) diff --git a/scopesim/optics/echelle.py b/scopesim/optics/echelle.py index 27a57831..bca69a90 100644 --- a/scopesim/optics/echelle.py +++ b/scopesim/optics/echelle.py @@ -7,10 +7,10 @@ class GratingSetup: def __init__( self, *, - groove_length:u.Quantity=None, alpha:float=None, delta:float = None, beta_center:float = None, + groove_length: u.Quantity = None, empiric_blaze_factor: float = 1.0, guess_littrow: tuple = None, ): From 8254953a21300adb2b8e1d68fba3c44ddd7d1d2b Mon Sep 17 00:00:00 2001 From: Yashvi-Sharma Date: Thu, 5 Mar 2026 15:53:06 -0800 Subject: [PATCH 3/8] Stub function for analytical spectral efficiency effect, updated ZShooter_v2 spec mode config with Echelle trace effects --- scopesim/effects/spectral_efficiency.py | 61 ++++++++++++++++++++++++- scopesim/effects/spectral_trace_list.py | 4 +- 2 files changed, 61 insertions(+), 4 deletions(-) diff --git a/scopesim/effects/spectral_efficiency.py b/scopesim/effects/spectral_efficiency.py index 35fe76fc..35962675 100644 --- a/scopesim/effects/spectral_efficiency.py +++ b/scopesim/effects/spectral_efficiency.py @@ -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 @@ -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__) @@ -127,3 +129,58 @@ def plot(self): axes.legend() return fig + + +class EchelletteSpectralEfficiency(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() + self.efficiencies = {} + + def _generate_efficiency_curve_func(self) -> Callable: + trace_params_table = self.table + # TODO: Call spectrograph and grating setup, return a callable? that takes wavelength as input and returns trace efficiency array? Or return a TERCurve object directly? + # For now, just return a dummy function that returns 1 for all wavelengths + def dummy_efficiency_curve(trace_id, wavelength): + return np.ones_like(wavelength) + return dummy_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 \ No newline at end of file diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index 2a298de1..4e84383b 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -524,7 +524,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 @@ -555,7 +555,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) From 94a8367471978eaa63265e6d203840def6c9fa5f Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 6 Mar 2026 10:23:04 -0800 Subject: [PATCH 4/8] Add spectrograph factory and enhance echelle parameter handling. add greater ability to guess cross-disperser. --- scopesim/optics/echelle.py | 284 ++++++++++++++++++++++++++++++++----- 1 file changed, 251 insertions(+), 33 deletions(-) diff --git a/scopesim/optics/echelle.py b/scopesim/optics/echelle.py index bca69a90..758e1767 100644 --- a/scopesim/optics/echelle.py +++ b/scopesim/optics/echelle.py @@ -4,15 +4,106 @@ import astropy.units as u import logging + +def spectrograph_factory(min_wave: float|u.Quantity, max_wave: float|u.Quantity, focal_len: float|u.Quantity, + design_res: float, echelle_angle: float|u.Quantity, min_order: int, max_order: int, + echelle_groove_length: float|u.Quantity, + pix_per_res_elem: float, disp_npix: int, xdisp_npix: int, pix_size: float|u.Quantity, + xdisp_groove_length: float|u.Quantity = 0.0, xdisp_beta_center: float|u.Quantity = 0.0): + """ + + Parameters + ---------- + min_wave: float|u.Quantity + Minimum wavelength. If float, assumed to be in nm. + max_wave: float|u.Quantity + Maximum wavelength. If float, assumed to be in nm. + focal_len: float|u.Quantity + Focal length. If float, assumed to be in mm. + design_res: float|int + Design resolution (dimensionless). + echelle_angle: float|u.Quantity + Echelle angle. If float, assumed to be in degrees. + min_order: int + Minimum order. + max_order: int + Maximum order. + echelle_groove_length: float|u.Quantity + Echelle groove length. If float, assumed to be in mm. + pix_per_res_elem: float|int + Pixels per resolution element. + disp_npix: int + Number of dispersion pixels. + xdisp_npix: int + Number of cross-dispersion pixels. + pix_size: float|u.Quantity + Pixel size. If float, assumed to be in um. + xdisp_groove_length: float|u.Quantity + Cross disperser groove length. If float, assumed to be in mm. + xdisp_beta_center: float|u.Quantity + Cross disperser beta center angle. If float, assumed to be in degrees. + + Returns + ------- + + """ + # Convert primitive types to quantities with appropriate units + if not isinstance(min_wave, u.Quantity): + min_wave = min_wave * u.nm + if not isinstance(max_wave, u.Quantity): + max_wave = max_wave * u.nm + if not isinstance(focal_len, u.Quantity): + focal_len = focal_len * u.mm + if not isinstance(echelle_angle, u.Quantity): + echelle_angle = echelle_angle * u.deg + if not isinstance(echelle_groove_length, u.Quantity): + echelle_groove_length = echelle_groove_length * u.mm + if not isinstance(pix_size, u.Quantity): + pix_size = pix_size * u.mm + if not isinstance(xdisp_groove_length, u.Quantity): + xdisp_groove_length = xdisp_groove_length * u.mm + if not isinstance(xdisp_beta_center, u.Quantity): + xdisp_beta_center = xdisp_beta_center * u.deg + + x_disp_len = (xdisp_npix*pix_size).to(u.mm) + + if not np.isfinite(xdisp_groove_length): + nu, xsdisp_angle = GratingSetup.estimate_xdisp_freq_and_angle(min_wave, max_wave, x_disp_len, focal_len) + cross_disperser = GratingSetup(alpha=xsdisp_angle, beta_center=xsdisp_angle, groove_length=1 / nu, + grating_type='vph', vph_center_wave=(min_wave + max_wave) / 2) + else: + if xdisp_beta_center != 0: # TODO probably hsould be a nan, but use of 0 angle incidence isn't likely + cross_disperser = GratingSetup(alpha=xdisp_beta_center, beta_center=xdisp_beta_center, + groove_length=xdisp_groove_length, grating_type='vph', + vph_center_wave=xdisp_groove_length.to(u.nm) * 2 * np.sin(xdisp_beta_center)) + else: + cross_disperser = GratingSetup(groove_length=xdisp_groove_length, + guess_littrow=(min_wave, max_wave, x_disp_len, focal_len), + grating_type='vph', vph_center_wave=(min_wave + max_wave) / 2) + + echelle_grating = GratingSetup(alpha=echelle_angle, beta_center=echelle_angle, delta=echelle_angle, + groove_length=echelle_groove_length) + + return SpectrographSetup(order_range=(min_order, max_order), + design_res=design_res, + pixels_per_res_elem=pix_per_res_elem, + focal_length=focal_len, + grating=echelle_grating, + detector=Detector(disp_npix, xdisp_npix, pix_size), + cross_disperser=cross_disperser) + + class GratingSetup: def __init__( self, *, - alpha:float=None, - delta:float = None, - beta_center:float = None, + alpha: u.Quantity = None, + delta: u.Quantity = None, + beta_center: u.Quantity = None, groove_length: u.Quantity = None, - empiric_blaze_factor: float = 1.0, - guess_littrow: tuple = None, + empiric_efficiency_factor: float = 1.0, + guess_littrow: tuple[u.Quantity,u.Quantity,u.Quantity,u.Quantity] = None, + grating_type: str = 'echelle', + vph_center_wave: u.Quantity = None ): """ Simulation of an echelle/echellete grating for spectrometer. @@ -21,34 +112,142 @@ def __init__( :param float delta: blaze angle in radians :param float beta_center: reflectance angle in radians :param u.Quantity groove_length: d, length/groove in units of wavelength (same as schroeder sigma) - :param float empiric_blaze_factor: empirically determined blaze factor to correct for grating efficiency relative to ideal peak - :param bool guess_littrow: if not None, passed onto estimate_xdisp_angle + :param float empiric_efficiency_factor: empirically determined blaze factor to correct for grating efficiency relative to ideal peak + :param bool guess_littrow: if not None, passed onto estimate_xdisp_angle_with_groove_len """ + self.grating_type = grating_type.lower() + assert self.grating_type in ('echelle', 'vph'), f'grating_type must be echelle or vph, not {self.grating_type}' + if guess_littrow: - alpha, beta_center = self.estimate_xdisp_angle(guess_littrow[0], guess_littrow[1], guess_littrow[2], + alpha, beta_center = self.estimate_xdisp_angle_with_groove_len(guess_littrow[0], guess_littrow[1], guess_littrow[2], guess_littrow[3], groove_length) self.alpha = alpha - self.delta = beta_center + self.delta = None if self.grating_type == 'vph' else beta_center self.beta_center = beta_center else: - self.alpha = alpha*u.rad - self.delta = delta*u.rad - self.beta_center = beta_center*u.rad - self.d = groove_length - self.empiric_blaze_factor = empiric_blaze_factor + self.alpha = alpha.to(u.rad) + self.delta = delta.to(u.rad) if delta is not None else None + self.beta_center = beta_center.to(u.rad) + self.d = groove_length.to(u.mm) + self.empiric_efficiency_factor = empiric_efficiency_factor + + # Approximate reality vph_constant = np.pi * deltan * d / bragg_angle + # Empircally we would use with peak efficiency near center of range + self.vph_constant = np.pi/2*vph_center_wave.to(u.nm) if self.grating_type == 'vph' else None @staticmethod - def estimate_xdisp_angle(l_min, l_max, detector_length, focal_length, d): + def estimate_xdisp_angle_with_groove_len( + l_min: float | u.Quantity, + l_max: float | u.Quantity, + detector_length: float | u.Quantity, + focal_length: float | u.Quantity, + d: u.Quantity) -> tuple[u.Quantity, u.Quantity]: + """ + Estimate cross-disperser angle and beta center for a given groove length. + + This method calculates the optimal incident angle (alpha) and central reflection + angle (beta) for a cross-disperser grating to match the wavelength range to the + detector length at the focal plane. + + Parameters + ---------- + l_min : float | u.Quantity + Minimum wavelength. If float, assumed to be in nm. + l_max : float | u.Quantity + Maximum wavelength. If float, assumed to be in nm. + detector_length : float | u.Quantity + Length of the detector. If float, assumed to be in mm. + focal_length : float | u.Quantity + Focal length of the system. If float, assumed to be in mm. + d : u.Quantity + Groove length (spacing) of the grating. + + Returns + ------- + tuple[u.Quantity, u.Quantity] + A tuple containing: + - alpha: incident angle in radians + - beta_center: central reflection angle in radians + """ + # Convert floats to quantities with appropriate units + if not isinstance(l_min, u.Quantity): + l_min = l_min * u.nm + if not isinstance(l_max, u.Quantity): + l_max = l_max * u.nm + if not isinstance(detector_length, u.Quantity): + detector_length = detector_length * u.mm + if not isinstance(focal_length, u.Quantity): + focal_length = focal_length * u.mm + if not isinstance(d, u.Quantity): + d = d * u.mm + # assume we want a dispersion at the midpoint that matches the range to the detector length - t = detector_length/focal_length - k = (l_max-l_min)/d - x = k*(1+np.sqrt(1+t**2))/2/t - u=np.arccos(np.sqrt(x)) - v=np.arcsin(k/2/np.sqrt(x)) + t = detector_length / focal_length + k = (l_max - l_min) / d + x = k * (1 + np.sqrt(1 + t ** 2)) / 2 / t + u = np.arccos(np.sqrt(x)) + v = np.arcsin(k / 2 / np.sqrt(x)) beta_max = u+v beta_min = u-v - alpha = np.arcsin(l_max/d -np.sin(beta_max)) - return alpha, (beta_min+beta_max)/2 + alpha = np.arcsin(l_max / d - np.sin(beta_max)) + return alpha.to(u.rad), ((beta_min + beta_max) / 2).to(u.rad) + + @staticmethod + def estimate_xdisp_freq_and_angle(l_min: float | u.Quantity, + l_max: float | u.Quantity, + detector_length: float | u.Quantity, + focal_length: float | u.Quantity) -> tuple[u.Quantity, u.Quantity]: + """ + Estimate cross-disperser littrow angle and groove frequency. + + This method calculates the analytical solution to the system of equations + + nu l_max == Sin(beta)+Sin(beta+deltabeta) nu l_min == Sin(beta)+Sin(beta-deltabeta) + cos[deltabeta] = f/sqrt(f^2 + (d/2)^) + sin[deltabeta] = (d/2)/sqrt(f^2 + (d/2)^) + + to match the initial and final wavelengths to the extent of the detector length + + + Parameters + ---------- + l_min : float | u.Quantity + Minimum wavelength. If float, assumed to be in nm. + l_max : float | u.Quantity + Maximum wavelength. If float, assumed to be in nm. + detector_length : float | u.Quantity + Length of the detector. If float, assumed to be in mm. + focal_length : float | u.Quantity + Focal length of the system. If float, assumed to be in mm. + + Returns + ------- + tuple[u.Quantity, u.Quantity] + A tuple containing: + - ruling frequency + - littrow angle + """ + if not isinstance(l_min, u.Quantity): + l_min = l_min * u.nm + if not isinstance(l_max, u.Quantity): + l_max = l_max * u.nm + if not isinstance(detector_length, u.Quantity): + detector_length = detector_length * u.mm + if not isinstance(focal_length, u.Quantity): + focal_length = focal_length * u.mm + + f = focal_length + p = detector_length + li, lf = l_min, l_max + a = 4*f**2 + p**2 + rta = np.sqrt(a) + b = lf - li + bsq = b**2 + c = lf**2 + li**2 + d = 4 * f**2 * bsq + 2 * f * rta * bsq + p**2 * c + beta = np.arccos(((2 * f + rta) * b) / np.sqrt(2 * d)) + nu = (np.sqrt(2) * p * (p**2 + 4 * f * (2 * f + rta))) / ((2 * f + rta) * np.sqrt(a * d)) + return nu.to('mm-1'), beta.to('rad') def __str__(self) -> str: return (f"alpha={np.rad2deg(self.alpha):.2f}\n" @@ -64,6 +263,8 @@ def blaze(self, beta, m): :param m: order :return: grating throughput out of 1 """ + if self.grating_type == 'vph': + raise ValueError('Cannot calculate blaze for vph gratings.') k = np.cos(beta) * np.cos(self.alpha - self.delta) / (np.cos(self.alpha) * np.cos(beta - self.delta)) k[k > 1] = 1 # q1 = np.cos(alpha) / np.cos(alpha - delta) @@ -73,7 +274,17 @@ def blaze(self, beta, m): rho = np.cos(self.delta) if self.alpha < self.delta else np.cos(self.alpha) / np.cos(self.alpha - self.delta) # 2 different rho depending on whether alpha or delta is larger ret = k * np.sinc((m * rho * q4).value) ** 2 # omit np.pi as np.sinc includes it - return self.empiric_blaze_factor * ret + return self.empiric_efficiency_factor * ret + + def vph_efficiency(self, wavelength): + """ + :param wavelength: wavelength(s) as u.Quantity units of length + :return: efficiency of the grating + """ + if self.grating_type != 'vph': + raise ValueError('Cannot calculate vph_effiency for echelle.') + + return self.empiric_efficiency_factor * np.sin(self.vph_constant/wavelength)**2 def beta(self, wave, m): """ @@ -136,7 +347,7 @@ class SpectrographSetup: def __init__( self, order_range: tuple, - design_res: u.Quantity, + design_res: float, pixels_per_res_elem: float, focal_length: u.Quantity, grating: GratingSetup, @@ -166,7 +377,7 @@ def __init__( self.design_res = design_res self.grating = grating self.detector = detector - self.focal_length = focal_length + self.focal_length = focal_length.to(u.mm) self.pixel_scale = np.arctan(self.detector.pixel_size / self.focal_length) self.beta_central_pixel = self.grating.beta_center self.nord = int(self.m_max - self.m0 + 1) @@ -283,6 +494,13 @@ def blaze(self, wave): """ return self.grating.blaze(self.grating.beta(wave, self.orders[:, None]), self.orders[:, None]) + def xdisp_efficiency(self, wave): + """ + :param wave: wavelength + :return: efficiency of the cross disperser + """ + return self.grating.vph_efficiency(wave) + def mean_blaze_eff_est(self, n=10): """ :param n: @@ -445,14 +663,14 @@ def angular_dispersion(self): beta = self.beta_for_x_pixel(np.arange(self.detector.n_pix_x) + .5) return self.grating.angular_dispersion(self.orders[:, None], beta) - # @property - # def design_res(self): - # """ - # :return: design resolution for spectrometer - # Assumes that m0 (the longest wavelength order) FSR fills detector with some sampling - # """ - # dlambda = self.fsr(self.m0) / self.detector.n_pix_x * self.nominal_pixels_per_res_elem - # return self.l0 / dlambda + @property + def implied_design_res(self): + """ + :return: design resolution for spectrometer + Assumes that m0 (the longest wavelength order) FSR fills detector with some sampling + """ + dlambda = self.fsr(self.m0) / self.detector.n_pix_x * self.nominal_pixels_per_res_elem + return self.l0 / dlambda @property def average_res(self): From 9b80cc3af3a19c33b58521e7f35596eb52fb30e2 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 6 Mar 2026 10:23:15 -0800 Subject: [PATCH 5/8] Fix logging format for selector wheel warning message. --- scopesim/effects/selector_wheel.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scopesim/effects/selector_wheel.py b/scopesim/effects/selector_wheel.py index 3bebfe17..dd4c11a8 100644 --- a/scopesim/effects/selector_wheel.py +++ b/scopesim/effects/selector_wheel.py @@ -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 From d5a107f1533952095885f9b325b1157900b8cb5d Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 6 Mar 2026 10:24:19 -0800 Subject: [PATCH 6/8] Implement initial computation of echelle efficiencies from analytical calculations. Echellette->Echelle --- scopesim/effects/spectral_efficiency.py | 50 ++++++++++++++++++++----- 1 file changed, 41 insertions(+), 9 deletions(-) diff --git a/scopesim/effects/spectral_efficiency.py b/scopesim/effects/spectral_efficiency.py index 35962675..dc498d5d 100644 --- a/scopesim/effects/spectral_efficiency.py +++ b/scopesim/effects/spectral_efficiency.py @@ -131,7 +131,7 @@ def plot(self): return fig -class EchelletteSpectralEfficiency(Effect): +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" @@ -140,16 +140,48 @@ class EchelletteSpectralEfficiency(Effect): def __init__(self, **kwargs): super().__init__(**kwargs) - self.efficiency_generator = self._generate_efficiency_curve_func() + self.efficiency_generator = self._generate_efficiency_curve_func(DataContainer(filename=kwargs.pop('filename'))) self.efficiencies = {} - def _generate_efficiency_curve_func(self) -> Callable: - trace_params_table = self.table - # TODO: Call spectrograph and grating setup, return a callable? that takes wavelength as input and returns trace efficiency array? Or return a TERCurve object directly? - # For now, just return a dummy function that returns 1 for all wavelengths - def dummy_efficiency_curve(trace_id, wavelength): - return np.ones_like(wavelength) - return dummy_efficiency_curve + 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.""" From ba9ed1293acee593f44bfdcc0dba189e891ba8a4 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 6 Mar 2026 10:25:17 -0800 Subject: [PATCH 7/8] Fix incorrect dictionary method usage in spectral trace initialization --- scopesim/effects/spectral_trace_list.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index 4e84383b..f684e06f 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -555,7 +555,7 @@ class EchelleSpectralTraceList(SpectralTraceList): def __init__(self, **kwargs): check_keys(kwargs, self.required_keys, action="error") - trace_params = DataContainer(filename=kwargs.pop['filename']) + trace_params = DataContainer(filename=kwargs.pop('filename')) hdulist = self._generate_trace_hdulist(trace_params) kwargs["hdulist"] = hdulist super().__init__(**kwargs) From 45a699c65e2397b11498c7fad4940f450460a24e Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 6 Mar 2026 10:25:30 -0800 Subject: [PATCH 8/8] Refactor spectral trace initialization and parameter handling for echelle spectrograph factory. --- scopesim/effects/spectral_trace_list.py | 42 ++++++++++--------------- 1 file changed, 16 insertions(+), 26 deletions(-) diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index f684e06f..32c91746 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -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): @@ -539,7 +538,7 @@ class EchelleSpectralTraceList(SpectralTraceList): # xdisp_freq_unit : mm # slitwidth_unit : arcsec - 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 + 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 @@ -590,34 +589,25 @@ def _generate_trace_hdulist(self, trace_params): 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']) - - echelle_grating = 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']) - - ss = echelle.SpectrographSetup(order_range=(min_order, max_order), - design_res=design_res, - pixels_per_res_elem=row['fwhm'] * u.Unit(trace_params.meta["fwhm_unit"]), - focal_length=focal_len, - grating=echelle_grating, - detector=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)