From 008b9761ab911fc2417b9c151a309e7ecec36a46 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:39:12 -0700 Subject: [PATCH 1/2] Correct analytical echelle angles, placement and physical slit dispersion Source commits: acfa5a71, aefbe60d, 87d6f581, 4c64fe55, f2c9c01c, 4c10ecdf, e2cefad8 (configured optical-table values) --- scopesim/effects/spectral_efficiency.py | 22 +-- scopesim/effects/spectral_trace_list.py | 142 ++++++++++---------- scopesim/optics/echelle.py | 91 ++++++++++--- scopesim/tests/tests_optics/test_echelle.py | 102 ++++++++++++++ 4 files changed, 263 insertions(+), 94 deletions(-) create mode 100644 scopesim/tests/tests_optics/test_echelle.py diff --git a/scopesim/effects/spectral_efficiency.py b/scopesim/effects/spectral_efficiency.py index 472389889..2baae5bd7 100644 --- a/scopesim/effects/spectral_efficiency.py +++ b/scopesim/effects/spectral_efficiency.py @@ -11,7 +11,7 @@ from .effects import Effect from .ter_curves import TERCurve -from ..utils import figure_factory, get_logger, check_keys +from ..utils import figure_factory, get_logger, check_keys, from_currsys from .data_container import DataContainer from ..optics import echelle @@ -149,7 +149,7 @@ def __init__(self, **kwargs): self.efficiency_generator = self._generate_efficiency_curve_func() def _generate_efficiency_curve_func(self) -> Callable: - trace_params = self.table + trace_params = from_currsys(self.table, self.cmds) spectrographs = {} for row in trace_params: prefix = row["prefix"] # note trance ids are assumed to be prefix_{order} @@ -173,13 +173,17 @@ def _generate_efficiency_curve_func(self) -> Callable: 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'] + prescription = {} + if 'collimator1_focal_length' in trace_params.colnames: + for key in ('collimator1_focal_length', 'collimator2_focal_length', 'entrance_slit_scale'): + prescription[key] = row[key] * u.Unit(trace_params.meta[key + '_unit']) spectrographs[prefix] = 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, - dispersion_focal_len=dispersion_focal_len) + dispersion_focal_len=dispersion_focal_len, **prescription) self._spectrographs = spectrographs def efficiency_curve(trace_id, wavelength): @@ -197,9 +201,12 @@ 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 + if obj.hdu is None: + wave = obj.cube_waveset + else: + 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 efficiency = self.efficiency_generator(trace_id, wave) params = {"description": trace_id} @@ -207,8 +214,7 @@ def apply_to(self, obj, **kwargs): effic_curve = TERCurve(array_dict={"wavelength": wave, "transmission": efficiency}, **params) self.efficiencies[trace_id] = effic_curve - obj.hdu.data *= effic_curve.throughput(wave).value[:, None, None] - return obj + return effic_curve.apply_to(obj) def plot(self): """Plot the grating efficiencies.""" diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index 6597ebb20..b10cbb49d 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -678,6 +678,19 @@ class EchelleSpectralTraceList(SpectralTraceList): nIR 2 0 40 24 970 2500 20000 64.2 225 225 4.7 0.7 10.0 10 0.015 4096 4096 45 175 10 x ---------------------------------------------------------------- + Optional ``detector_offset_x`` and ``detector_offset_y`` columns specify + the rotated analytical camera origin relative to detector center, positive + right/up. Both columns and their ``*_unit`` metadata are required together; + values must have length units. Explicit zeros place that origin at detector + center. Offsets are applied after rotation, before detector clipping, with + no automatic recentering. Omitting both columns retains automatic placement. + + With ``physical_slit_dispersion=True``, supply ``collimator1_focal_length``, + ``collimator2_focal_length`` and ``entrance_slit_scale`` columns with their + ``*_unit`` metadata. The entrance scale is length per on-sky arcsecond. + These determine the illumination angle and asymmetric pupil relay; + the along-slit plate scale is not used for across-slit dispersion. + The calculated traces are stored in the same HDUList format as required by SpectralTraceList, and supplied to the parent class through hdulist kwarg. @@ -690,10 +703,13 @@ def __init__(self, cmds=None, **kwargs): self.cmds = cmds trace_param_filename = kwargs.pop("filename") + self.physical_slit_dispersion = kwargs.pop("physical_slit_dispersion", False) + self.slit_samples = kwargs.pop("slit_samples", 3) save_generated_hdulist = kwargs.pop("save_generated_hdulist", False) generated_hdulist_filename = kwargs.pop( "generated_hdulist_filename", "analytical_echelle_traces.fits") trace_params = DataContainer(filename=trace_param_filename) + trace_params.table = from_currsys(trace_params.table, self.cmds) hdulist = self._generate_trace_hdulist(trace_params) if save_generated_hdulist: output_path = Path(generated_hdulist_filename).expanduser() @@ -712,6 +728,7 @@ def _generate_trace_hdulist(self, trace_params): hdul[0].header["EDATA"] = 2 trace_ids, ap_ids, im_ids = [], [], [] + self._spectrographs = {} for row in trace_params.table: prefix = row["prefix"] min_order = row['m0'] - row['n'] + 1 @@ -740,13 +757,19 @@ def _generate_trace_hdulist(self, trace_params): 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'] + prescription = {} + if 'collimator1_focal_length' in trace_params.table.colnames: + for key in ('collimator1_focal_length', 'collimator2_focal_length', 'entrance_slit_scale'): + prescription[key] = row[key] * u.Unit(trace_params.meta[key + '_unit']) + nominal_width = row['nominal_slit_width'] * u.Unit(trace_params.meta['nominal_slit_width_unit']) + pix_per_res_elem = (nominal_width * prescription['entrance_slit_scale'] * focal_len / prescription['collimator2_focal_length'] / pix_size).to_value(u.dimensionless_unscaled) 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, - dispersion_focal_len=dispersion_focal_len) + dispersion_focal_len=dispersion_focal_len, **prescription) _warn_if_echelle_design_res_inconsistent( prefix, design_res, @@ -757,7 +780,7 @@ def _generate_trace_hdulist(self, trace_params): ) slit_edge = (row['slitlength'] / 2) * u.Unit(trace_params.meta["slitlength_unit"]) - slit_pos = np.linspace(-slit_edge, slit_edge, num=3) + slit_pos = np.linspace(-slit_edge, slit_edge, num=self.slit_samples) if "plate_scale" in trace_params.table.colnames: plate_scale = row["plate_scale"] * u.Unit( trace_params.meta["plate_scale_unit"]) @@ -776,74 +799,52 @@ def _generate_trace_hdulist(self, trace_params): ).to_value(u.deg) cang = np.cos(np.deg2rad(detector_angle)) sang = np.sin(np.deg2rad(detector_angle)) - detector_x_padding = max( - 0.0, - ( - disp_npix * abs(cang) - + xdisp_npix * abs(sang) - - disp_npix - ) / 2, - ) - slit_x_padding = ( - abs(np.tan(np.deg2rad(detector_angle))) - * float(np.nanmax(np.abs(slit_offset_pix))) - ) - raw_x_min = 0.5 - detector_x_padding - slit_x_padding - raw_x_max = disp_npix - 0.5 + detector_x_padding + slit_x_padding - - def raw_detector_pixels(wave, order): - x = ss.wavelength_to_x_pixel(wave, order) - y = ss.wavelength_to_y_pixel(wave) - x = u.Quantity(x, copy=False).to_value( - u.dimensionless_unscaled) - y = u.Quantity(y, copy=False).to_value( - u.dimensionless_unscaled) - xpix = np.broadcast_to(x, (slit_offset_pix.size, wave.size)) - ypix = y[None, :] + slit_offset_pix[:, None] - return xpix, ypix - - rotated_x_min = rotated_y_min = np.inf - rotated_x_max = rotated_y_max = -np.inf - # The analytical echelle layout is relative, not yet detector - # placed. Rotate the footprint first, then place that rotated - # footprint on the detector. Detector padding is not part of the - # optical trace geometry; downstream detector/FOV code is - # responsible for clipping padded display or extraction regions. - for i, order in enumerate(ss.orders): - raw_x = np.linspace( - raw_x_min, raw_x_max, num=max(disp_npix, 2)) - wave = ss.x_pixel_to_wavelength(raw_x, order) - xpix, ypix = raw_detector_pixels(wave, order) - xrot = cang * xpix - sang * ypix - yrot = sang * xpix + cang * ypix - rotated_x_min = min(rotated_x_min, float(np.nanmin(xrot))) - rotated_x_max = max(rotated_x_max, float(np.nanmax(xrot))) - rotated_y_min = min(rotated_y_min, float(np.nanmin(yrot))) - rotated_y_max = max(rotated_y_max, float(np.nanmax(yrot))) - rotated_x_center = rotated_x_min + (rotated_x_max - rotated_x_min) / 2 - rotated_y_center = rotated_y_min + (rotated_y_max - rotated_y_min) / 2 - - def rotated_detector_pixels(wave, order): - xpix, ypix = raw_detector_pixels(wave, order) - xrot = cang * xpix - sang * ypix - yrot = sang * xpix + cang * ypix - return ( - xrot - rotated_x_center + disp_npix / 2, - yrot - rotated_y_center + xdisp_npix / 2, - ) + ss._detector_cos, ss._detector_sin = cang, sang + ss._slit_offset_pix = slit_offset_pix + self._spectrographs[prefix] = ss + if 'detector_offset_x' in trace_params.table.colnames or 'detector_offset_y' in trace_params.table.colnames: + # Multiplying an ASCII table's masked cell by a Unit yields + # zero, so reject missing values before attaching units. + if np.ma.is_masked(row['detector_offset_x']) or np.ma.is_masked(row['detector_offset_y']): + raise ValueError('Both detector offsets must have values.') + offset_x = row['detector_offset_x'] * u.Unit(trace_params.meta['detector_offset_x_unit']) + offset_y = row['detector_offset_y'] * u.Unit(trace_params.meta['detector_offset_y_unit']) + # Offsets locate the rotated camera origin relative to detector + # center. The raw mapping includes the unrotated pixel origin. + rotated_x_center = cang*disp_npix/2 - sang*xdisp_npix/2 - (offset_x/pix_size).to_value(u.dimensionless_unscaled) + rotated_y_center = sang*disp_npix/2 + cang*xdisp_npix/2 - (offset_y/pix_size).to_value(u.dimensionless_unscaled) + + # Project the translated detector rectangle back onto model x; + # slit offsets are in model y and cancel in this projection. + raw_x_center = cang*rotated_x_center + sang*rotated_y_center + raw_x_halfspan = (abs(cang)*disp_npix + abs(sang)*xdisp_npix)/2 + candidate_raw_x = np.linspace(raw_x_center - raw_x_halfspan, raw_x_center + raw_x_halfspan, disp_npix + 1) + else: + # Preserve the original automatic placement and sampling for + # parameter files that do not specify a detector position. + detector_x_padding = max(0., (disp_npix*abs(cang) + xdisp_npix*abs(sang) - disp_npix)/2) + slit_x_padding = abs(np.tan(np.deg2rad(detector_angle))) * np.max(np.abs(slit_offset_pix)) + raw_x_min = .5 - detector_x_padding - slit_x_padding + raw_x_max = disp_npix - .5 + detector_x_padding + slit_x_padding + candidate_raw_x = np.linspace(raw_x_min, raw_x_max, num=max(disp_npix, 2)) + rotated_x_min = rotated_y_min = np.inf + rotated_x_max = rotated_y_max = -np.inf + for order in ss.orders: + wave = ss.x_pixel_to_wavelength(candidate_raw_x, order) + xrot, yrot = ss._rotated_detector_pixels(wave, order) + rotated_x_min = min(rotated_x_min, xrot.min()) + rotated_x_max = max(rotated_x_max, xrot.max()) + rotated_y_min = min(rotated_y_min, yrot.min()) + rotated_y_max = max(rotated_y_max, yrot.max()) + rotated_x_center = rotated_x_min + (rotated_x_max - rotated_x_min)/2 + rotated_y_center = rotated_y_min + (rotated_y_max - rotated_y_min)/2 + ss._detector_center = (rotated_x_center, rotated_y_center) def on_detector(wave, order): - xpix, ypix = rotated_detector_pixels(wave, order) - return ( - (xpix >= 0) - & (xpix <= disp_npix) - & (ypix >= 0) - & (ypix <= xdisp_npix) - ) + xpix, ypix = ss._detector_pixels(wave, order) + return (xpix >= 0) & (xpix <= disp_npix) & (ypix >= 0) & (ypix <= xdisp_npix) - for i, order in enumerate(ss.orders): - candidate_raw_x = np.linspace( - raw_x_min, raw_x_max, num=max(disp_npix, 2)) + for order in ss.orders: candidate_wave = ss.x_pixel_to_wavelength(candidate_raw_x, order) valid_wave = np.any( on_detector(candidate_wave, order), axis=0) @@ -872,7 +873,7 @@ def on_detector(wave, order): continue s = np.tile(slit_pos, wave.size).reshape(wave.size, slit_pos.size).T.ravel() w = np.tile(wave, slit_offset_pix.size) - xpix, ypix = rotated_detector_pixels(wave, order) + xpix, ypix = ss._detector_pixels(wave, order) xval = (xpix.ravel() - disp_npix / 2) * pix_size.to(u.mm) yval = (ypix.ravel() - xdisp_npix / 2) * pix_size.to(u.mm) @@ -881,7 +882,12 @@ def on_detector(wave, order): 'x': xval, 'y': yval}) + if self.physical_slit_dispersion: + order_table['dlam_deta'] = ss.slit_dispersion(w, order) + trace_hdu = fits.BinTableHDU(order_table) + if self.physical_slit_dispersion: + trace_hdu.header['SLITCOL'] = 'dlam_deta' trace_hdu.header['DISPDIR'] = row['dispdir'] trace_hdu.header["EXTNAME"] = f'{prefix}_{order:d}' trace_hdu.header["DESIGNR"] = ( diff --git a/scopesim/optics/echelle.py b/scopesim/optics/echelle.py index 4c7614e63..709b2fd2e 100644 --- a/scopesim/optics/echelle.py +++ b/scopesim/optics/echelle.py @@ -11,7 +11,10 @@ def spectrograph_factory(min_wave: float|u.Quantity, max_wave: float|u.Quantity, 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, - dispersion_focal_len: float|u.Quantity = None): + dispersion_focal_len: float|u.Quantity = None, + collimator1_focal_length: u.Quantity = None, + collimator2_focal_length: u.Quantity = None, + entrance_slit_scale: u.Quantity = None): """ Parameters @@ -47,6 +50,11 @@ def spectrograph_factory(min_wave: float|u.Quantity, max_wave: float|u.Quantity, dispersion_focal_len: float|u.Quantity Effective focal length for the echelle/main-dispersion detector coordinate. If float, assumed to be in mm. Defaults to focal_len. + collimator1_focal_length, collimator2_focal_length: u.Quantity + Physical collimator lengths. When supplied together, derive the + dispersion focal length from focal_len times their ratio. + entrance_slit_scale: u.Quantity + Entrance-slit length per on-sky angle, required with the collimators. Returns ------- @@ -98,6 +106,9 @@ def spectrograph_factory(min_wave: float|u.Quantity, max_wave: float|u.Quantity, pixels_per_res_elem=pix_per_res_elem, focal_length=focal_len, dispersion_focal_length=dispersion_focal_len, + collimator1_focal_length=collimator1_focal_length, + collimator2_focal_length=collimator2_focal_length, + entrance_slit_scale=entrance_slit_scale, grating=echelle_grating, detector=Detector(disp_npix, xdisp_npix, pix_size), cross_disperser=cross_disperser) @@ -191,16 +202,21 @@ def estimate_xdisp_angle_with_groove_len( 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 - beta_center_guess = np.arccos(np.sqrt(x)) - beta_half_range = np.arcsin(k / 2 / np.sqrt(x)) - beta_max = beta_center_guess + beta_half_range - beta_min = beta_center_guess - beta_half_range - alpha = np.arcsin(l_max / d - np.sin(beta_max)) - return alpha.to(u.rad), ((beta_min + beta_max) / 2).to(u.rad) + # The camera maps beta_center +/- half_angle to +/- detector_length/2. + # Subtracting the two grating equations eliminates the incidence angle: + # (l_max - l_min)/d = 2*cos(beta_center)*sin(half_angle). + if l_max <= l_min or detector_length <= 0 * u.mm or focal_length <= 0 * u.mm or d <= 0 * u.mm: + raise ValueError("Wavelength range, detector length, focal length and groove spacing must be positive.") + half_angle = np.arctan(detector_length / (2 * focal_length)) + cos_center = (l_max - l_min) / (2 * d * np.sin(half_angle)) + if not 0 < cos_center <= 1: + raise ValueError("The wavelength range cannot fit this detector with the specified groove spacing.") + beta_center = np.arccos(cos_center) + sin_alpha = l_min / d - np.sin(beta_center - half_angle) + if beta_center + half_angle >= 90 * u.deg or not -1 <= sin_alpha <= 1: + raise ValueError("The requested detector span requires a non-propagating grating angle.") + alpha = np.arcsin(sin_alpha) + return alpha.to(u.rad), beta_center.to(u.rad) @staticmethod def estimate_xdisp_freq_and_angle(l_min: float | u.Quantity, @@ -364,6 +380,9 @@ def __init__( detector: Detector, cross_disperser: GratingSetup = None, dispersion_focal_length: u.Quantity = None, + collimator1_focal_length: u.Quantity = None, + collimator2_focal_length: u.Quantity = None, + entrance_slit_scale: u.Quantity = None, ): """ @@ -398,6 +417,12 @@ def __init__( if dispersion_focal_length is None else dispersion_focal_length.to(u.mm) ) + if collimator1_focal_length is not None: + self.collimator1_focal_length = collimator1_focal_length.to(u.mm) + self.collimator2_focal_length = collimator2_focal_length.to(u.mm) + self.entrance_slit_scale = entrance_slit_scale.to(u.mm / u.arcsec) + # Angular dispersion is amplified by the asymmetric pupil relay. + self.dispersion_focal_length = self.focal_length * self.collimator1_focal_length / self.collimator2_focal_length self.pixel_scale = np.arctan( self.detector.pixel_size / self.dispersion_focal_length) self.beta_central_pixel = self.grating.beta_center @@ -432,6 +457,22 @@ def set_beta_center(self, beta, littrow: bool = False): if littrow: self.grating.alpha = beta + def slit_dispersion(self, wavelength, order): + """Apparent wavelength shift per across-slit arcsecond at fixed incident wavelength. + + m*lambda = d*(sin(alpha) + sin(beta)); increasing entrance-slit eta + increases alpha and reduces the wavelength inferred at nominal alpha. + The first collimator sets this angle; the pupil relay is already in + the camera's dispersion focal length and must not enter a second time. + """ + derivative = -self.grating.d * np.cos(self.grating.alpha) * self.entrance_slit_scale / (order * self.collimator1_focal_length) + return np.ones(np.shape(wavelength)) * derivative.to(u.um / u.arcsec) + + def projected_slit_width(self, wavelength, order, slit_width): + """Unrotated echelle-coordinate slit width, including local anamorphism.""" + half_shift = self.slit_dispersion(wavelength, order) * slit_width / 2 + return abs(self.wavelength_to_x_pixel(wavelength + half_shift, order) - self.wavelength_to_x_pixel(wavelength - half_shift, order)) * u.pixel + @property def orders(self): try: @@ -558,17 +599,15 @@ def order_mask(self, wave, fsr_edge: bool = False): def edge_wave(self, fsr=True): """ :param fsr: True to return the FSR edge wavelengths, False to return the detector edges - :return: wavelengths at the edges of the detector (or FSR) for each order, always include central pixel wavelength + :return: the two endpoint wavelengths for each order """ - pix = np.array([0, self.detector.n_pix_x // 2, self.detector.n_pix_x-1]) + .5 - fiducial_waves = self.x_pixel_to_wavelength(pix, self.orders[:, None]) + pixels = np.array([0, self.detector.n_pix_x//2, self.detector.n_pix_x - 1]) + .5 + fiducial_waves = self.x_pixel_to_wavelength(pixels, self.orders[:, None]) if not fsr: return fiducial_waves[:, [0, -1]] - + # Retain the original central-pixel wavelength convention for the FSR. central_fsr = fiducial_waves[:, 1] / self.orders - fsr_edges = (u.Quantity([-central_fsr / 2, central_fsr / 2]) + fiducial_waves[:, 1]).T - - return fsr_edges + return (u.Quantity([-central_fsr/2, central_fsr/2]) + fiducial_waves[:, 1]).T def central_wave(self, order): """ @@ -616,6 +655,22 @@ def x_pixel_wavelengths(self, pixel_phase=.5): def wavelength_to_y_pixel(self, wave): return self.y_pixel_for_beta(self.cross_disperser.beta(wave, 1)) + def _rotated_detector_pixels(self, wave, order): + x = self.wavelength_to_x_pixel(wave, order).to_value(u.dimensionless_unscaled) + y = self.wavelength_to_y_pixel(wave).to_value(u.dimensionless_unscaled) + xpix = np.broadcast_to(x, (self._slit_offset_pix.size, wave.size)) + ypix = y[None, :] + self._slit_offset_pix[:, None] + return self._detector_cos*xpix - self._detector_sin*ypix, self._detector_sin*xpix + self._detector_cos*ypix + + def _detector_pixels(self, wave, order): + """The trace generator's placed coordinates for all configured slit nodes. + + Unlike the clipped trace table, this mapping can also inspect the full + FSR outside the detector. It uses the same rotation and placement. + """ + x, y = self._rotated_detector_pixels(wave, order) + return x - self._detector_center[0] + self.detector.n_pix_x/2, y - self._detector_center[1] + self.detector.n_pix_y/2 + def y_pixel_to_wavelength(self, pixel, m): """ :param pixel: pixel diff --git a/scopesim/tests/tests_optics/test_echelle.py b/scopesim/tests/tests_optics/test_echelle.py new file mode 100644 index 000000000..2370f03dd --- /dev/null +++ b/scopesim/tests/tests_optics/test_echelle.py @@ -0,0 +1,102 @@ +"""Planar grating geometry checked through the camera's forward projection.""" + +import numpy as np +import pytest +from astropy import units as u + +from scopesim.optics.echelle import GratingSetup, spectrograph_factory + + +# Wavelength endpoints [nm], cross-disperser density [mm^-1], detector height [pix]. +CHANNELS = [(310, 430, 2000, 4096), (430, 650, 1100, 4096), (650, 980, 625, 4096), + (960, 1340, 350, 2048), (1470, 1830, 300, 2048), (1990, 2450, 270, 2048)] + + +def test_asymmetric_pupil_slit_dispersion_and_anamorphism(): + instrument = spectrograph_factory(960, 1340, 225, 18000, 65.6, 58, 80, 1/24, 4.46, + 2048, 2048, .015, xdisp_groove_length=1/350, xdisp_beta_center=10.9, + collimator1_focal_length=1500*u.mm, collimator2_focal_length=1250*u.mm, entrance_slit_scale=.5306*u.mm/u.arcsec) + assert instrument.dispersion_focal_length == 270*u.mm + center = instrument.central_wave(65) + derivative = instrument.slit_dispersion(center, 65) + assert derivative < 0*u.um/u.arcsec + for slit, expected in ((.7, 17800), (.33, 37800)): + resolution = (center / abs(derivative * slit*u.arcsec)).to_value(u.dimensionless_unscaled) + assert resolution == pytest.approx(expected, rel=.01) + width = instrument.projected_slit_width(center, 65, .7*u.arcsec) + assert width.to_value(u.pixel) == pytest.approx(4.45704, rel=1e-5) + wavelengths = center * np.array([.99, 1., 1.01]) + widths = instrument.projected_slit_width(wavelengths, 65, .7*u.arcsec) + assert np.all(np.diff(widths) > 0*u.pixel) + half_angle = (.35*u.arcsec * instrument.entrance_slit_scale / instrument.collimator1_focal_length).value * u.rad + beta_low = np.arcsin(65*wavelengths/instrument.grating.d - np.sin(instrument.grating.alpha + half_angle)) + beta_high = np.arcsin(65*wavelengths/instrument.grating.d - np.sin(instrument.grating.alpha - half_angle)) + edge_width = abs(instrument.x_pixel_for_beta(beta_high) - instrument.x_pixel_for_beta(beta_low)) + # Independent finite slit edges retain the same anamorphism as dispersion. + np.testing.assert_allclose(widths.to_value(u.pixel), edge_width.value, rtol=2e-7) + step = center*1e-7 + dispersion = (instrument.wavelength_to_x_pixel(wavelengths+step, 65)-instrument.wavelength_to_x_pixel(wavelengths-step, 65))/(2*step) + measured_delta = widths.to_value(u.pixel)/dispersion + np.testing.assert_allclose(measured_delta.to_value(u.nm), abs(derivative*.7*u.arcsec).to_value(u.nm), rtol=2e-7) + eta = .0825*u.arcsec + dalpha = (eta * instrument.entrance_slit_scale / instrument.collimator1_focal_length).to_value(u.dimensionless_unscaled) * u.rad + exact = instrument.grating.d / 65 * (np.sin(instrument.grating.alpha) - np.sin(instrument.grating.alpha + dalpha)) + assert (derivative * eta).to_value(u.nm) == pytest.approx(exact.to_value(u.nm), rel=1e-4) + + +def test_incomplete_physical_prescription_fails(): + with pytest.raises(AttributeError): + spectrograph_factory(960, 1340, 225, 18000, 65.6, 58, 80, 1/24, 4.46, 2048, 2048, .015, + xdisp_groove_length=1/350, collimator1_focal_length=1500*u.mm) + + +@pytest.mark.parametrize("low,high,density,height", CHANNELS) +def test_fixed_groove_endpoints_fill_camera(low, high, density, height): + instrument = spectrograph_factory(low, high, 225, 18000, 65.6, 40, 50, 1/24, 4.46, + 2048, height, .015, xdisp_groove_length=1/density, dispersion_focal_len=270) + wavelength = np.array([low, high]) * u.nm + pixels = instrument.wavelength_to_y_pixel(wavelength) + # Algebraic inverse/forward evaluation; 1e-7 pixel allows floating-point cancellation. + np.testing.assert_allclose(pixels.value, [0, height], atol=1e-7, rtol=0) + recovered = instrument.y_pixel_to_wavelength(pixels, 1) + np.testing.assert_allclose(recovered.to_value(u.nm), wavelength.value, atol=1e-9, rtol=0) + + +def test_units_and_signed_incidence(): + numeric = GratingSetup.estimate_xdisp_angle_with_groove_len(1470, 1830, 30.72, 225, 1/300) + physical = GratingSetup.estimate_xdisp_angle_with_groove_len(1.47*u.um, 1.83*u.um, 3.072*u.cm, .225*u.m, (1/300)*u.mm) + assert physical[0] < 0 * u.rad + for plain, quantity in zip(numeric, physical): + assert plain.unit == u.rad + assert quantity.unit == u.rad + assert plain.to_value(u.rad) == pytest.approx(quantity.to_value(u.rad), abs=1e-12) + + +@pytest.mark.parametrize("inputs", [(310, 430, 1, 225, .0005), (3000, 3010, 61.44, 225, .0001), + (430, 310, 61.44, 225, .0005), (310, 430, 0, 225, .0005)]) +def test_impossible_geometry_fails(inputs): + with pytest.raises(ValueError): + GratingSetup.estimate_xdisp_angle_with_groove_len(*inputs) + + +def test_explicit_angle_keeps_configured_geometry(): + instrument = spectrograph_factory(960, 1340, 225, 18000, 65.6, 57, 80, 1/24, 4.46, + 2048, 2048, .015, xdisp_groove_length=1/350, xdisp_beta_center=12) + assert instrument.cross_disperser.alpha.to_value(u.deg) == pytest.approx(12) + assert instrument.cross_disperser.beta_center.to_value(u.deg) == pytest.approx(12) + + +def test_free_frequency_endpoints_fill_camera(): + instrument = spectrograph_factory(960, 1340, 225, 18000, 65.6, 57, 80, 1/24, 4.46, + 2048, 2048, .015, xdisp_groove_length=np.nan) + pixels = instrument.wavelength_to_y_pixel(np.array([960, 1340]) * u.nm) + np.testing.assert_allclose(pixels.value, [0, 2048], atol=1e-7, rtol=0) + + +def test_fsr_retains_the_central_pixel_wavelength_convention(): + instrument = spectrograph_factory(1990, 2450, 225, 18000, 65.6, 31, 38, 1/24, 4.46, + 2048, 2048, .015, xdisp_groove_length=1/270, dispersion_focal_len=270) + edges = instrument.edge_wave(fsr=True) + center = instrument.x_pixel_to_wavelength(instrument.detector.n_pix_x//2 + .5, instrument.orders) + np.testing.assert_allclose(edges.mean(axis=1).to_value(u.nm), center.to_value(u.nm), rtol=1e-12) + np.testing.assert_allclose(np.diff(edges, axis=1).ravel().to_value(u.nm), (center/instrument.orders).to_value(u.nm), rtol=1e-12) From a744f1ce7ff2c35b47e0d0ddb01941bba19b6c7c Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:37:10 -0700 Subject: [PATCH 2/2] Preserve arm exposures, charge binning and independent detector readouts Source commits: 2a68fbc0, 0ac61b83, ef55346d, b3c91073, c4f0ec17, 307faba5, 2780412e, 1f528b4a, 7e30187e --- scopesim/detector/detector.py | 8 +- scopesim/effects/detector_list.py | 11 ++- scopesim/effects/electronic/__init__.py | 4 +- scopesim/effects/electronic/dmps.py | 14 +-- scopesim/effects/electronic/electrons.py | 13 +-- scopesim/effects/electronic/exposure.py | 3 + scopesim/effects/electronic/noise.py | 19 +++- scopesim/effects/electronic/pixels.py | 77 +++++++++++---- scopesim/effects/fits_headers.py | 9 ++ scopesim/effects/mosaic_trace_list.py | 4 + scopesim/optics/optical_train.py | 95 +++++++------------ .../tests/tests_effects/test_OnChipBinning.py | 75 +++++++++++++++ .../test_ReadoutNoiseFromTable.py | 25 +++++ .../test_detector_arm_configuration.py | 53 +++++++++++ .../tests_effects/test_mosaic_trace_list.py | 23 +++++ .../tests/tests_optics/test_OpticalTrain.py | 73 ++++++++++++-- 16 files changed, 398 insertions(+), 108 deletions(-) create mode 100644 scopesim/tests/tests_effects/test_OnChipBinning.py create mode 100644 scopesim/tests/tests_effects/test_ReadoutNoiseFromTable.py create mode 100644 scopesim/tests/tests_effects/test_detector_arm_configuration.py diff --git a/scopesim/detector/detector.py b/scopesim/detector/detector.py index b99c5dba1..9e1336217 100644 --- a/scopesim/detector/detector.py +++ b/scopesim/detector/detector.py @@ -25,6 +25,7 @@ class Detector: def __init__(self, header, cmds=None, **kwargs): + self._native_header = header.copy() self._hdu = ImageHDU(header=header, data=zeros_from_header(header)) self.meta = {} self.meta.update(header) @@ -51,10 +52,9 @@ def extract_from(self, image_plane, spline_order=1, reset=True): self.header.update(sky_wcs.to_header()) def reset(self): - """Reset internal HDU data to all-zeros-array.""" - # The detector might have been converted to integers by the - # Quantization effect, so it is not possible to use zeros_like. - self._hdu.data = np.zeros(self._hdu.data.shape) + """Reset charge and restore the configured, unbinned detector geometry.""" + self._hdu = ImageHDU(header=self._native_header, data=zeros_from_header(self._native_header)) + self.meta.update(self._native_header) @property def hdu(self): diff --git a/scopesim/effects/detector_list.py b/scopesim/effects/detector_list.py index 7cb075146..4a3c5ea55 100644 --- a/scopesim/effects/detector_list.py +++ b/scopesim/effects/detector_list.py @@ -330,16 +330,17 @@ def detector_headers(self, ids=None): hdr["PC2_1"] = -sang hdr["PC2_2"] = cang - # hdr["GAIN"] = row["gain"] - if "id" in row: - hdr["DET_ID"] = row["id"] - hdr["EXTNAME"] = f"DET_{row['id']}" - row_dict = { f"HIERARCH {col.upper()}" if len(col) > 8 else col: row[col] for col in row.colnames } hdr.update(row_dict) + if "id" in row.colnames: + # Keep routing identifiers numeric when detector metadata is + # serialized to FITS (NumPy integers otherwise become strings). + hdr["ID"] = int(row["id"]) + hdr["DET_ID"] = int(row["id"]) + hdr["EXTNAME"] = f"DET_{row['id']}" hdrs.append(hdr) return hdrs diff --git a/scopesim/effects/electronic/__init__.py b/scopesim/effects/electronic/__init__.py index c35dfe06d..e086896cf 100644 --- a/scopesim/effects/electronic/__init__.py +++ b/scopesim/effects/electronic/__init__.py @@ -52,9 +52,9 @@ def apply_to(self, obj, **kwargs): from .electrons import LinearityCurve, ADConversion, InterPixelCapacitance from .noise import (Bias, PoorMansHxRGReadoutNoise, BasicReadoutNoise, - PixelResponseNonUniformity, ShotNoise, DarkCurrent) + PixelResponseNonUniformity, ShotNoise, DarkCurrent, ReadoutNoiseFromTable) from .exposure import AutoExposure, ExposureIntegration, ExposureOutput -from .pixels import ReferencePixelBorder, BinnedImage, UnequalBinnedImage +from .pixels import ReferencePixelBorder, BinnedImage, UnequalBinnedImage, OnChipBinning from .dmps import DetectorModePropertiesSetter diff --git a/scopesim/effects/electronic/dmps.py b/scopesim/effects/electronic/dmps.py index 296703a88..2cf3e2130 100644 --- a/scopesim/effects/electronic/dmps.py +++ b/scopesim/effects/electronic/dmps.py @@ -23,6 +23,9 @@ class DetectorModePropertiesSetter(Effect): mode_properties : dict A dictionary containing the DET parameters to be changed for each mode. See below for an example yaml entry. + mode_key : str + Command selecting the mode, by default ``!OBS.detector_readout_mode``. + Separate keys allow independent detector arms to select their modes. Examples -------- @@ -70,21 +73,20 @@ def __init__(self, **kwargs): check_keys(self.meta, self.required_keys, action="error") self.mode_properties = kwargs["mode_properties"] + self.mode_key = kwargs.get("mode_key", "!OBS.detector_readout_mode") def apply_to(self, obj, **kwargs): logger.debug("Entering DetectorModePropertiesSetter.apply_to") - mode_name = kwargs.get( - "detector_readout_mode", - from_currsys("!OBS.detector_readout_mode", self.cmds), - ) - + mode_name = kwargs.get("detector_readout_mode", + from_currsys(self.mode_key, + self.cmds)) if isinstance(obj, ImagePlane) and mode_name == "auto": mode_name = self.select_mode(obj, **kwargs) logger.info("Detector mode set to %s", mode_name) self.meta["detector_readout_mode"] = mode_name props_dict = self.mode_properties[mode_name] - self.cmds["!OBS.detector_readout_mode"] = mode_name + self.cmds[self.mode_key] = mode_name for key, value in props_dict.items(): self.cmds[key] = value logger.debug("DMPS setting %s: %s", key, value) diff --git a/scopesim/effects/electronic/electrons.py b/scopesim/effects/electronic/electrons.py index c2bee25b2..ac3e13236 100644 --- a/scopesim/effects/electronic/electrons.py +++ b/scopesim/effects/electronic/electrons.py @@ -263,7 +263,7 @@ def __init__(self, **kwargs): super().__init__(**kwargs) params = { "dtype": "uint16", - "gain": 1. # default, usually overridden from yaml + "gain": "!DET.gain" } self.meta.update(params) self.meta.update(kwargs) @@ -297,17 +297,18 @@ def _get_gain(self, det_id: int) -> Real: # Apply the gain value (copy from DarkCurrent) # Note that this does not cater for the case where the gain is given # as a plain dictionary. Should we implement that? - if hasattr(self.cmds["!DET.gain"], "dic"): - gain = self.cmds["!DET.gain"].dic[det_id] + configured_gain = from_currsys(self.meta["gain"], self.cmds) + if hasattr(configured_gain, "dic"): + gain = configured_gain.dic[det_id] logger.info(f"Detector {det_id}: applying gain {gain}") return gain - if isinstance(self.cmds["!DET.gain"], Real): - gain = self.cmds["!DET.gain"] + if isinstance(configured_gain, Real): + gain = configured_gain logger.info(f"Applying gain {gain}") return gain raise ValueError( f"{self.__class__.__name__}.meta['gain'] must be either " - f"dict or float, but is {self.cmds['!DET.gain']}") + f"dict or float, but is {configured_gain}") def _apply_to_det(self, det: Detector) -> None: """Subclasses can override if more params needed in call.""" diff --git a/scopesim/effects/electronic/exposure.py b/scopesim/effects/electronic/exposure.py index 342e9ac67..f11237375 100644 --- a/scopesim/effects/electronic/exposure.py +++ b/scopesim/effects/electronic/exposure.py @@ -287,6 +287,9 @@ def apply_to(self, obj, **kwargs): "readout call." ) + obj.header["EXPTIME"] = (dit * ndit, "Total accumulated integration time [s]") + obj.header['DIT'] = (dit, 'Integration time per exposure [s]') + obj.header['NDIT'] = ndit obj.data = obj.data * dit * ndit return obj diff --git a/scopesim/effects/electronic/noise.py b/scopesim/effects/electronic/noise.py index 555218549..8c6aef4fa 100644 --- a/scopesim/effects/electronic/noise.py +++ b/scopesim/effects/electronic/noise.py @@ -54,7 +54,7 @@ def __call__( return data + dark_level * dit * ndit def _get_dark_level(self, det_id: int) -> float: - dark_level = float(from_currsys(self.meta["value"], self.cmds)) + dark_level = from_currsys(self.meta["value"], self.cmds) if isinstance(dark_level, Real): return dark_level if isinstance(dark_level, Mapping): @@ -286,6 +286,23 @@ def _pseudo_random_field( return image +class ReadoutNoiseFromTable(BasicReadoutNoise): + """Interpolate measured read noise in electrons as a function of DIT [s].""" + + required_keys = {"filename", "dit", "ndit"} + + def apply_to(self, det, **kwargs): + if not isinstance(det, Detector): + return det + dit = from_currsys(self.meta["dit"], self.cmds) + times = self.table["dit"] + if not times[0] <= dit <= times[-1]: + raise ValueError(f"DIT {dit} s is outside the measured read-noise table") + self.meta["noise_std"] = np.interp(dit, times, self.table["read_noise"]) + det.header["RDNOISE"] = (self.meta["noise_std"], "Read noise per integration [electron]") + return super().apply_to(det, **kwargs) + + class PixelResponseNonUniformity(RandomEffect): """Pixel Response Non-Uniformity (PRNU). diff --git a/scopesim/effects/electronic/pixels.py b/scopesim/effects/electronic/pixels.py index e26b26e96..efb966cf2 100644 --- a/scopesim/effects/electronic/pixels.py +++ b/scopesim/effects/electronic/pixels.py @@ -93,19 +93,55 @@ class BinnedImageBase(ElectronicEffect): z_order: ClassVar[tuple[int, ...]] = (870,) + def __call__(self, data: ArrayLike) -> NDArray: + bx, by = self.bin_factors + ny, nx = data.shape + return data.reshape(ny//by, by, nx//bx, bx).sum(axis=(1, 3)) + + def bin_header(self, native_header): + """WCS of the summed pixel blocks, also used by trace coordinates.""" + bx, by = self.bin_factors + if not isinstance(bx, int) or not isinstance(by, int) or min(bx, by) < 1: + raise ValueError('Bin factors must be positive integers') + if 'allowed_bins' in self.meta and [bx, by] not in self.meta['allowed_bins']: + raise ValueError(f'Unsupported binning: {bx} x {by}') + header = native_header.copy() + for axis, factor in ((1, bx), (2, by)): + if header[f'NAXIS{axis}'] % factor: + raise ValueError('Binning must divide the full detector dimensions') + header[f'NAXIS{axis}'] //= factor + for suffix in [''] + list('ABCDEFGHIJKLMNOPQRSTUVWXYZ'): + if f'CRPIX1{suffix}' not in header: + continue + for axis, factor in ((1, bx), (2, by)): + header[f'CRPIX{axis}{suffix}'] = (header[f'CRPIX{axis}{suffix}'] - .5)/factor + .5 + if f'CD1_1{suffix}' in header: + for i in (1, 2): + for j, factor in ((1, bx), (2, by)): + header[f'CD{i}_{j}{suffix}'] *= factor + else: + header[f'CDELT1{suffix}'] *= bx + header[f'CDELT2{suffix}'] *= by + if f'PC1_2{suffix}' in header: + header[f'PC1_2{suffix}'] *= by/bx + if f'PC2_1{suffix}' in header: + header[f'PC2_1{suffix}'] *= bx/by + header['BINX'], header['BINY'] = bx, by + return header + + def _apply_to_det(self, det): + header = self.bin_header(det.header) + det.data = self(det.data) + det.header.update(header) + det.meta.update(header) + class BinnedImage(BinnedImageBase): required_keys = {"bin_size"} - def __call__(self, data: ArrayLike) -> NDArray: - height, width = data.shape - binned_data = data.reshape(( - height//self.bin_size, - self.bin_size, - width//self.bin_size, - self.bin_size - )).sum(axis=3).sum(axis=1) - return binned_data + @property + def bin_factors(self): + return self.bin_size, self.bin_size @property def bin_size(self) -> int: @@ -115,15 +151,9 @@ def bin_size(self) -> int: class UnequalBinnedImage(BinnedImageBase): required_keys = {"binx","biny"} - def __call__(self, data: ArrayLike) -> NDArray: - height, width = data.shape - binned_data = data.reshape(( - height//self.biny, - self.biny, - width//self.binx, - self.binx - )).sum(axis=3).sum(axis=1) - return binned_data + @property + def bin_factors(self): + return self.binx, self.biny @property def binx(self) -> int: @@ -132,3 +162,14 @@ def binx(self) -> int: @property def biny(self) -> int: return from_currsys(self.meta["biny"], self.cmds) + + +class OnChipBinning(UnequalBinnedImage): + """Sum charge in FITS x/y pixel blocks before adding read noise. + + Unlike a software sum after readout, each output pixel receives one read + noise draw. Configure this effect between charge integration and read noise. + WCS reference pixels refer to the centers of the new pixel blocks. + """ + + z_order = (815,) diff --git a/scopesim/effects/fits_headers.py b/scopesim/effects/fits_headers.py index 9110f3331..987117edf 100644 --- a/scopesim/effects/fits_headers.py +++ b/scopesim/effects/fits_headers.py @@ -153,6 +153,9 @@ class ExtraFitsKeywords(Effect): header['HIERARCH ESO DET DIT'] = (5, '[s] exposure length') header['EXTNAME'] = "DET2.DATA" + A group may additionally specify ``detector_id``. It then applies only to + HDULists containing that ``DET_ID``, including their primary headers. + Resolved and un-resolved keywords --------------------------------- ScopeSim uses bang-strings to resolve global parameters. @@ -285,6 +288,12 @@ def apply_to(self, hdul, **kwargs): return hdul for dic in self.dict_list: + # Each detector array produces its own HDUList. A detector selector + # also permits arm-specific metadata in that array's primary HDU. + if "detector_id" in dic: + detector_ids = [hdu.header["DET_ID"] for hdu in hdul if "DET_ID" in hdu.header] + if dic["detector_id"] not in detector_ids: + continue resolved = flatten_dict(dic.get("keywords", {}), resolve=True, optics_manager=opt_train) unresolved = flatten_dict(dic.get("unresolved_keywords", {})) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index 34e6bb66e..c74c9d214 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -255,6 +255,9 @@ def apply_to(self, det, **kwargs): image = det.data detwcs = WCS(det.header, key="D") + # ExposureIntegration has already recorded the actual accumulated time. + # New table HDUs need those cards, but not the detector-image WCS. + exposure_header = fits.Header(card for card in det.header.cards if card.keyword in ("DIT", "NDIT", "EXPTIME")) output_format = from_currsys(self.meta["format"], self.cmds) @@ -309,6 +312,7 @@ def apply_to(self, det, **kwargs): else: logger.error("Unknown output format: %s", output_format) + det.header.update(exposure_header) return det diff --git a/scopesim/optics/optical_train.py b/scopesim/optics/optical_train.py index cd8ab65bd..6974bbe9f 100644 --- a/scopesim/optics/optical_train.py +++ b/scopesim/optics/optical_train.py @@ -122,6 +122,7 @@ def __init__(self, cmds=None): self.detector_managers = [] self.yaml_dicts = None self._last_source = None + self._readout_id = 0 if cmds is not None: self.load(cmds) @@ -234,12 +235,13 @@ def trace_detector_coordinates(self, *, xi=0 * u.arcsec, explicit definition of which output image owns each trace coordinate. """ from ..effects import ( - BinnedImage, + DetectorModePropertiesSetter, + ImagePlaneBackground, Rotate90CCD, SelectorWheel, SpectralTraceList, - UnequalBinnedImage, ) + from ..effects.electronic.pixels import BinnedImageBase trace_lists = self.optics_manager.get_all(SpectralTraceList) if len(trace_lists) != 1: @@ -248,14 +250,18 @@ def trace_detector_coordinates(self, *, xi=0 * u.arcsec, "SpectralTraceList.") trace_list = trace_lists[0] - if self.optics_manager.image_plane_effects: - raise NotImplementedError( - "Trace detector coordinates do not yet support image-plane " - "effects.") - if self.optics_manager.detector_array_effects: - raise NotImplementedError( - "Trace detector coordinates do not yet support detector-array " - "effects.") + for effect in self.optics_manager.image_plane_effects: + selected = effect.wheel_effects.values() if isinstance(effect, SelectorWheel) else [effect] + if any(not isinstance(item, ImagePlaneBackground) for item in selected): + raise NotImplementedError('Trace detector coordinates only support additive image-plane backgrounds.') + for effect in self.optics_manager.detector_array_effects: + if not isinstance(effect, DetectorModePropertiesSetter): + raise NotImplementedError('Trace detector coordinates do not yet support this detector-array effect.') + # Mode setters update commands, not image coordinates. Resolve an + # explicit mode before reading any binning command it can set. + if from_currsys(effect.mode_key, self.cmds) == 'auto': + raise NotImplementedError('Select an explicit detector readout mode before inspecting coordinates.') + effect.apply_to(None) fovs_by_trace = {} for fov in self.fov_manager.fovs: @@ -326,31 +332,14 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): "detectors on one image plane.") detector = detector_manager[0] detector_id = detector.header["ID"] + readout_header = detector._native_header.copy() for effect in self.optics_manager.detector_effects: - selected_effect = ( - effect.wheel_effects.get(detector_id) - if isinstance(effect, SelectorWheel) - else effect - ) + selected_effect = effect.get_effect(detector_id) if isinstance(effect, SelectorWheel) else effect if selected_effect is None: continue - if isinstance(selected_effect, BinnedImage): - bin_size = from_currsys( - selected_effect.meta["bin_size"], self.cmds) - if bin_size != 1: - raise NotImplementedError( - "Trace detector coordinates do not yet support " - "post-extraction detector binning.") - if isinstance(selected_effect, UnequalBinnedImage): - binx = from_currsys( - selected_effect.meta["binx"], self.cmds) - biny = from_currsys( - selected_effect.meta["biny"], self.cmds) - if binx != 1 or biny != 1: - raise NotImplementedError( - "Trace detector coordinates do not yet support " - "post-extraction detector binning.") + if isinstance(selected_effect, BinnedImageBase): + readout_header = selected_effect.bin_header(readout_header) if isinstance(selected_effect, Rotate90CCD): rotations = from_currsys( selected_effect.meta["rotations"], self.cmds) @@ -359,7 +348,7 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): "Trace detector coordinates do not yet support " "post-extraction detector rotation.") - detector_wcs = WCS(detector.header, key="D") + detector_wcs = WCS(detector._native_header, key="D") image_plane_wcs = image_plane._det_wcs if ( not np.allclose(image_plane_wcs.wcs.pc, np.eye(2)) @@ -429,13 +418,7 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): f"Trace {trace_name!r} has overlapping FOV wavelength " "ranges.") - ( - _, - _, - x_grid_mm, - y_grid_mm, - subimage_wcs, - ) = trace._focal_plane_grid(fov) + _, _, x_grid_mm, y_grid_mm, subimage_wcs = trace._focal_plane_grid(fov) if ( not np.allclose(subimage_wcs.wcs.pc, np.eye(2)) or not np.allclose( @@ -450,16 +433,7 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): x_mm = trace.xilam2x(xi_flat[in_fov], wave_flat[in_fov]) y_mm = trace.xilam2y(xi_flat[in_fov], wave_flat[in_fov]) - subimage_x = ( - (x_mm - x_grid_mm[0]) - / (x_grid_mm[-1] - x_grid_mm[0]) - * (len(x_grid_mm) - 1) - ) - subimage_y = ( - (y_mm - y_grid_mm[0]) - / (y_grid_mm[-1] - y_grid_mm[0]) - * (len(y_grid_mm) - 1) - ) + subimage_x, subimage_y = subimage_wcs.all_world2pix(x_mm, y_mm, 0) subimage_shape = np.array( [len(x_grid_mm), len(y_grid_mm)]) @@ -477,6 +451,8 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): f"Requested wavelengths for trace {trace_name!r} " "do not map through a configured FOV.") + focal_x, focal_y = detector_wcs.all_pix2world(x_pix, y_pix, 0) + x_pix, y_pix = WCS(readout_header, key='D').all_world2pix(focal_x, focal_y, 0) point_count = wave_grid.size columns["readout_index"].extend( np.full(point_count, readout_index)) @@ -534,6 +510,7 @@ def observe(self, orig_source=None, update=True, **kwargs): """ if update: self.update(**kwargs) + self._readout_id = 0 if self.cmds["!SIM.random.seed"] is None: # Concretize seed from system entropy to keep constant during the @@ -731,7 +708,7 @@ def readout( self, filename: Path | str | None = None, reset: bool = True, - roid: int = 0, + roid: int | None = None, **kwargs ) -> list[HDUList]: """ @@ -746,13 +723,10 @@ def readout( reset : bool, optional If true (the default), do not keep any parameters set during any previous readout. This is usually what you want. - roid : int, optional - Readout ID, used for seeding randomness. If you run multiple - readouts on the same observation and want the random effects to - produce different results, set this parameter with e.g. the loop - counter. Simple integers are fine (0, 1, 2, ...), because this gets - wrangled into a proper random seed downstream. If you need to - reproduce a specific readout in a sequence, set e.g. `roid=4`. + roid : int or None, optional + Readout ID for reproducible noise. By default, successive calls + use 0, 1, 2, ... within each observation. Supply an explicit ID to + replay that noise realization with the same root seed. **kwargs : Any other kwargs passed to the detector effects. @@ -764,6 +738,9 @@ def readout( """ if reset: self.cmds.clear() + if roid is None: + roid = self._readout_id + self._readout_id += 1 # Hack to make sure AutoExposure and Quantization work properly. # Should probably be removed once #428 is fixed properly. @@ -774,7 +751,7 @@ def readout( # self.cmds["!OBS.ndit"] = None # TODO: This is still hacky but seems to work for now... params = {"exptime": None, "roid": roid} - params.update(kwargs) + params.update({key.removeprefix('!OBS.'): value for key, value in kwargs.items()}) if params["exptime"] is not None and params.get("dit") is None: params.update( dit=None, @@ -789,8 +766,8 @@ def readout( for key, value in params.items(): if is_bangkey(key): - # skip if already explicitly !ABC.xyz logger.debug("%s: %s given in kwargs", key, value) + self.cmds[key] = value continue # otherwise fallback to !OBS logger.debug("%s: %s given in kwargs, put in !OBS", key, value) diff --git a/scopesim/tests/tests_effects/test_OnChipBinning.py b/scopesim/tests/tests_effects/test_OnChipBinning.py new file mode 100644 index 000000000..ec3347f70 --- /dev/null +++ b/scopesim/tests/tests_effects/test_OnChipBinning.py @@ -0,0 +1,75 @@ +import numpy as np +import pytest +from scopesim.detector import Detector +from scopesim.optics.image_plane_utils import header_from_list_of_xy +from astropy.wcs import WCS +from scopesim.effects.electronic import OnChipBinning, UnequalBinnedImage, BasicReadoutNoise + + +def detector(nx=120, ny=80): + header = header_from_list_of_xy([-.015*nx/2, .015*nx/2], [-.015*ny/2, .015*ny/2], .015, 'D') + return Detector(header) + + + +@pytest.mark.parametrize('bx,by', [(1, 1), (1, 2), (2, 2)]) +def test_binning_charge_and_focal_coordinates(bx, by): + det = detector() + det.data[:] = 7 + original = WCS(det.header, key='D') + OnChipBinning(binx=bx, biny=by).apply_to(det) + assert det.data.shape == (80//by, 120//bx) + assert np.all(det.data == 7*bx*by) + assert det.data.sum() == 7*120*80 + binned = WCS(det.hdu.header, key='D') + expected = original.all_pix2world((np.array([10, 20])+.5)*bx-.5, (np.array([12, 23])+.5)*by-.5, 0) + np.testing.assert_allclose(binned.all_pix2world([10, 20], [12, 23], 0), expected, atol=1e-12) + det.reset() + assert det.data.shape == (80, 120) + assert det.header['CDELT2D'] == .015 + OnChipBinning(binx=bx, biny=by).apply_to(det) + assert det.data.shape == (80//by, 120//bx) + + + +def test_binning_rotated_cd_matrix(): + det = detector() + det.header['CD1_1D'], det.header['CD1_2D'] = .012, -.009 + det.header['CD2_1D'], det.header['CD2_2D'] = .009, .012 + wcs = WCS(det.header, key='D') + OnChipBinning(binx=1, biny=2).apply_to(det) + actual = WCS(det.header, key='D').all_pix2world(9, 15, 0) + np.testing.assert_allclose(actual, wcs.all_pix2world(9, 30.5, 0), atol=1e-12) + + + +def test_binning_rejects_unsupported_modes(): + with pytest.raises(ValueError, match='Unsupported'): + OnChipBinning(binx=2, biny=1, allowed_bins=[[1, 1], [1, 2], [2, 2]]).apply_to(detector()) + with pytest.raises(ValueError): + OnChipBinning(binx=1, biny=3).apply_to(detector()) + + + +def test_ccd_bin_does_not_multiply_amplifier_noise(): + det = detector(600, 400) + OnChipBinning(binx=2, biny=2).apply_to(det) + BasicReadoutNoise(noise_std=3, ndit=1, random_seed=17).apply_to(det) + assert det.data.std() == pytest.approx(3, rel=.02) + + +def test_software_sum_accumulates_independent_read_noise(): + det = detector(600, 400) + BasicReadoutNoise(noise_std=3, ndit=1, random_seed=17).apply_to(det) + UnequalBinnedImage(binx=2, biny=2).apply_to(det) + assert det.data.std() == pytest.approx(6, rel=.02) + assert det.header['CDELT1D'] == pytest.approx(.03) + assert det.header['CDELT2D'] == pytest.approx(.03) + + +def test_header_projection_matches_applied_binning(): + det = detector() + binning = OnChipBinning(binx=1, biny=2) + expected = binning.bin_header(det.header) + binning.apply_to(det) + assert det.header == expected diff --git a/scopesim/tests/tests_effects/test_ReadoutNoiseFromTable.py b/scopesim/tests/tests_effects/test_ReadoutNoiseFromTable.py new file mode 100644 index 000000000..ef372b486 --- /dev/null +++ b/scopesim/tests/tests_effects/test_ReadoutNoiseFromTable.py @@ -0,0 +1,25 @@ +import numpy as np +import pytest +from scopesim.detector import Detector +from scopesim.optics.image_plane_utils import header_from_list_of_xy +from scopesim.effects.electronic import ReadoutNoiseFromTable, ExposureOutput + + +def detector(nx=120, ny=80): + header = header_from_list_of_xy([-.015*nx/2, .015*nx/2], [-.015*ny/2, .015*ny/2], .015, 'D') + return Detector(header) + + + +@pytest.mark.parametrize('ndit', [1, 4, 9]) +def test_nir_averaging_preserves_signal_and_reduces_noise(tmp_path, ndit): + noise_file = tmp_path / 'noise.dat' + noise_file.write_text('dit read_noise\n10 12\n100 6\n') + det = detector(500, 400) + det.data[:] = 30*ndit + ReadoutNoiseFromTable(filename=str(noise_file), dit=55, ndit=ndit, random_seed=21).apply_to(det) + ExposureOutput(dit=55, ndit=ndit, mode='average').apply_to(det) + assert det.data.mean() == pytest.approx(30, abs=.1) + assert det.data.std() == pytest.approx(9/np.sqrt(ndit), rel=.02) + with pytest.raises(ValueError, match='outside'): + ReadoutNoiseFromTable(filename=str(noise_file), dit=101, ndit=ndit).apply_to(det) diff --git a/scopesim/tests/tests_effects/test_detector_arm_configuration.py b/scopesim/tests/tests_effects/test_detector_arm_configuration.py new file mode 100644 index 000000000..2b81134c4 --- /dev/null +++ b/scopesim/tests/tests_effects/test_detector_arm_configuration.py @@ -0,0 +1,53 @@ +import numpy as np +import pytest +from scopesim.detector import Detector +from scopesim.optics.image_plane_utils import header_from_list_of_xy +from astropy.io import fits +from scopesim.commands import UserCommands +from scopesim.effects import DetectorList, ExtraFitsKeywords, DetectorModePropertiesSetter +from scopesim.effects.electronic import ADConversion + + +def detector(nx=120, ny=80): + header = header_from_list_of_xy([-.015*nx/2, .015*nx/2], [-.015*ny/2, .015*ny/2], .015, 'D') + return Detector(header) + + + +def test_gain_uses_selected_arm_and_missing_gain_fails(): + det = detector() + det.data[:] = 12 + ADConversion(gain='!DET.gain_uvb', dtype='float32', cmds={'!DET.gain_uvb': .6}).apply_to(det) + np.testing.assert_allclose(det.data, 20) + with pytest.raises(ValueError, match='gain_uvb'): + ADConversion(gain='!DET.gain_uvb', dtype='float32', cmds={}).apply_to(det) + + + +def test_detector_specific_primary_headers(): + eff = ExtraFitsKeywords(header_dict=[{'detector_id': 1, 'ext_number': [0, 1], 'keywords': {'ARM': 'VIS'}}]) + for detector_id in [0, 1, 2]: + hdul = fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(header=fits.Header({'DET_ID': detector_id}))]) + eff.apply_to(hdul) + assert ('ARM' in hdul[0].header) == (detector_id == 1) + + + +def test_detector_layout_sets_real_id_header(): + effect = DetectorList(array_dict={'id': [2], 'x_cen': [0.], 'y_cen': [0.], 'x_size': [16], 'y_size': [8], 'pixel_size': [.018], 'angle': [0.], 'gain': [2.29]}, x_cen_unit='mm', y_cen_unit='mm', x_size_unit='pixel', y_size_unit='pixel', pixel_size_unit='mm', angle_unit='deg', gain_unit='electron/adu') + header = effect.detector_headers()[0] + assert header['DET_ID'] == 2 + assert header['EXTNAME'] == 'DET_2' + det = Detector(header) + assert det.meta['ID'] == 2 + assert det.hdu.header['DET_ID'] == 2 + + + +def test_independent_detector_readout_modes(): + cmds = UserCommands(properties={'!DET.mode_blue': 'slow', '!DET.mode_red': 'fast'}) + for arm in ['blue', 'red']: + modes = {'slow': {f'!DET.gain_{arm}': .6}, 'fast': {f'!DET.gain_{arm}': 1.8}} + DetectorModePropertiesSetter(mode_properties=modes, mode_key=f'!DET.mode_{arm}', cmds=cmds).apply_to(None) + assert cmds['!DET.gain_blue'] == .6 + assert cmds['!DET.gain_red'] == 1.8 diff --git a/scopesim/tests/tests_effects/test_mosaic_trace_list.py b/scopesim/tests/tests_effects/test_mosaic_trace_list.py index 39433af43..f5bef8e8d 100644 --- a/scopesim/tests/tests_effects/test_mosaic_trace_list.py +++ b/scopesim/tests/tests_effects/test_mosaic_trace_list.py @@ -12,6 +12,7 @@ MosaicSpectralTraceList, MosaicSpectralTrace, MosaicOutputFormat) +from scopesim.effects.electronic.exposure import ExposureIntegration from scopesim.optics.image_plane_utils import header_from_list_of_xy from scopesim.detector import Detector @@ -201,3 +202,25 @@ def test_image_returns_detector(self, apply_image): def test_image_detector_has_image(self, apply_image): assert isinstance(apply_image._hdu, fits.ImageHDU) + + +@pytest.mark.usefixtures("patch_mock_path_mosaic") +@pytest.mark.parametrize("output_format", ["collapse1d", "table", "image"]) +def test_output_format_preserves_accumulated_counts_and_exposure(output_format): + header = header_from_list_of_xy([-25, 25], [-2048, 2048], 1, "D") + detector = Detector(header) + rows = [2023, 2031, 2039, 2047, 2055, 2063, 2071] + detector.data[rows, :] = 1 + rate_sum = detector.data.sum() + ExposureIntegration(dit=17., ndit=3).apply_to(detector) + output = MosaicOutputFormat(filename="TRACE_MOS-LR-J.fits", format=output_format).apply_to(detector) + assert output.header["DIT"] == 17. + assert output.header["NDIT"] == 3 + assert output.header["EXPTIME"] == 51. + if output_format == "image": + counts = output.data.sum() + assert "CTYPE1D" in output.header + else: + counts = output.data["spectrum"].sum() + assert "CTYPE1D" not in output.header + assert counts == pytest.approx(rate_sum * 51.) diff --git a/scopesim/tests/tests_optics/test_OpticalTrain.py b/scopesim/tests/tests_optics/test_OpticalTrain.py index c1931499a..1585e350d 100644 --- a/scopesim/tests/tests_optics/test_OpticalTrain.py +++ b/scopesim/tests/tests_optics/test_OpticalTrain.py @@ -23,6 +23,11 @@ Effect, SpectralTraceList, UnequalBinnedImage, + OnChipBinning, + BasicReadoutNoise, + ImagePlaneBackground, + DetectorModePropertiesSetter, + SelectorWheel, ) from scopesim.utils import find_file @@ -197,6 +202,33 @@ def test_ignore_effects_works(self, cmds_with_ignore): class TestTraceDetectorCoordinates: + def test_explicit_detector_mode_resolves_its_binning_before_mapping(self): + train, _ = spectral_geometry_train() + before = train.trace_detector_coordinates() + train.cmds = UserCommands(properties={'!OBS.detector_readout_mode': 'slow'}) + setter = DetectorModePropertiesSetter(mode_properties={'slow': {'!DET.binx': 2}}, cmds=train.cmds) + train.optics_manager.detector_array_effects = [setter] + train.optics_manager.detector_effects = [OnChipBinning(binx='!DET.binx', biny=1, cmds=train.cmds)] + after = train.trace_detector_coordinates() + np.testing.assert_allclose(after['detector_x'], (before['detector_x'] + .5*u.pixel)/2 - .5*u.pixel) + np.testing.assert_array_equal(after['detector_y'], before['detector_y']) + + @pytest.mark.parametrize('selected', [False, True]) + def test_additive_background_preserves_trace_coordinates(self, selected): + train, _ = spectral_geometry_train() + before = train.trace_detector_coordinates() + background = SelectorWheel(selector_key='image_plane_id', wheel=[{'selector_value': [0], 'effect_class': 'ImagePlaneBackground', 'effect_kwargs': {'value': 1.}}]) if selected else ImagePlaneBackground(value=1.) + train.optics_manager.image_plane_effects = [background] + after = train.trace_detector_coordinates() + np.testing.assert_array_equal(before['detector_x'], after['detector_x']) + np.testing.assert_array_equal(before['detector_y'], after['detector_y']) + + def test_unsupported_image_plane_transform_is_not_silently_ignored(self): + train, _ = spectral_geometry_train() + train.optics_manager.image_plane_effects = [Effect()] + with pytest.raises(NotImplementedError, match='additive'): + train.trace_detector_coordinates() + def test_maps_live_trace_to_zero_origin_detector_pixels(self): train, _ = spectral_geometry_train() @@ -258,14 +290,19 @@ def test_rejects_multiple_detectors_on_one_image_plane(self): with pytest.raises(NotImplementedError, match="multiple detectors"): train.trace_detector_coordinates() - def test_rejects_post_extraction_binning(self): + @pytest.mark.parametrize('binning_class', [UnequalBinnedImage, OnChipBinning]) + def test_binned_coordinates_match_before_and_after_readout_geometry(self, binning_class): train, _ = spectral_geometry_train() - train.optics_manager.detector_effects = [ - UnequalBinnedImage(binx=2, biny=1), - ] - - with pytest.raises(NotImplementedError, match="binning"): - train.trace_detector_coordinates() + native = train.trace_detector_coordinates() + binning = binning_class(binx=2, biny=1) + train.optics_manager.detector_effects = [binning] + before = train.trace_detector_coordinates() + np.testing.assert_allclose(before['detector_x'], (native['detector_x'] + .5*u.pixel)/2 - .5*u.pixel) + np.testing.assert_allclose(before['detector_y'], native['detector_y']) + binning.apply_to(train.detector_managers[0][0]) + after = train.trace_detector_coordinates() + np.testing.assert_array_equal(before['detector_x'], after['detector_x']) + np.testing.assert_array_equal(before['detector_y'], after['detector_y']) @pytest.mark.slow @@ -378,6 +415,28 @@ def test_works_with_a_pointer_to_fits_imagehdu(self, cmds): @pytest.mark.usefixtures("patch_mock_path") class TestReadout: + def test_explicit_readout_keys_and_reset(self, unity_cmds): + opt = OpticalTrain(deepcopy(unity_cmds)) + for seconds in (3600., 300., 3600.): + opt.readout(**{'!OBS.dit': seconds, '!OBS.ndit': 1, '!OBS.dit_blue': seconds / 2}) + assert opt.cmds['!OBS.dit'] == seconds + assert opt.cmds['!OBS.dit_blue'] == seconds / 2 + opt.readout(dit=17., ndit=2) + assert opt.cmds['!OBS.dit'] == 17. + assert opt.cmds['!OBS.ndit'] == 2 + + def test_successive_noise_realizations_can_be_replayed(self, unity_cmds, unity_src): + commands = deepcopy(unity_cmds) + commands['!SIM.random.seed'] = 17 + opt = OpticalTrain(commands) + opt.optics_manager.add_effect(BasicReadoutNoise(noise_std=3., ndit=1, cmds=opt.cmds)) + opt.observe(unity_src) + first = opt.readout()[0][1].data + second = opt.readout()[0][1].data + replay = opt.readout(roid=0)[0][1].data + assert not np.array_equal(first, second) + np.testing.assert_array_equal(first, replay) + def test_readout_works_when_source_observed(self, unity_cmds, unity_src): opt = OpticalTrain(unity_cmds)