From a744f1ce7ff2c35b47e0d0ddb01941bba19b6c7c Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:37:10 -0700 Subject: [PATCH 1/7] 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) From 4296266bc2aadd12cee2ad4ba23940d23585d212 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:39:12 -0700 Subject: [PATCH 2/7] Use calibrated traces and physical slit and detector coordinates Source commits: fe052a1d, c76a4c21, 93237d85, c0e0dfdc, 7ffffb87, fa2598b7 (command arrays), f9dccd8f (spectral padding), 3264a6e6, ce5c8cd5, 06f64b09 --- scopesim/effects/apertures.py | 4 + scopesim/effects/illumination.py | 11 +- scopesim/effects/metis_lms_trace_list.py | 5 +- scopesim/effects/spectral_trace_list_utils.py | 180 +++++++++++------- scopesim/effects/ter_curves.py | 7 +- scopesim/tests/test_utils_functions.py | 9 + .../tests_effects/test_SkycalcTERCurve.py | 4 +- .../tests_effects/test_SpectralTraceList.py | 101 +++++++++- .../test_SpectralTraceListUtils.py | 8 +- .../test_calibrated_spectral_trace.py | 87 +++++++++ .../tests/tests_effects/test_local_skycalc.py | 24 +++ .../test_registered_background.py | 32 ++++ scopesim/utils.py | 5 +- 13 files changed, 392 insertions(+), 85 deletions(-) create mode 100644 scopesim/tests/tests_effects/test_calibrated_spectral_trace.py create mode 100644 scopesim/tests/tests_effects/test_local_skycalc.py create mode 100644 scopesim/tests/tests_effects/test_registered_background.py diff --git a/scopesim/effects/apertures.py b/scopesim/effects/apertures.py index aca858141..199914541 100644 --- a/scopesim/effects/apertures.py +++ b/scopesim/effects/apertures.py @@ -128,6 +128,8 @@ def apply_to(self, obj, **kwargs): for vol in obj.volumes: vol["meta"]["xi_min"] = min(x) * u.arcsec vol["meta"]["xi_max"] = max(x) * u.arcsec + vol["meta"]["eta_min"] = min(y) * u.arcsec + vol["meta"]["eta_max"] = max(y) * u.arcsec vol["meta"]["slit_name"] = self.meta["name"] if "offset_type" in self.meta: @@ -305,6 +307,8 @@ def apply_to(self, obj, **kwargs): # ..todo: HUGE HACK - Get rid of this! vol["meta"]["xi_min"] = row["left"] * u.arcsec vol["meta"]["xi_max"] = row["right"] * u.arcsec + vol["meta"]["eta_min"] = row["bottom"] * u.arcsec + vol["meta"]["eta_max"] = row["top"] * u.arcsec vol["conserve_image"] = row["conserve_image"] vol["shape"] = row["shape"] vol["angle"] = row["angle"] diff --git a/scopesim/effects/illumination.py b/scopesim/effects/illumination.py index d4ff12bfe..65b5b8181 100644 --- a/scopesim/effects/illumination.py +++ b/scopesim/effects/illumination.py @@ -538,6 +538,7 @@ class ImagePlaneBackground(Effect): trace list. The supplied value or model should already include downstream throughput and detector QE. For tapered QE maps, use an average or otherwise representative positional QE when deriving this non-dispersed background. + A FITS ``filename`` and ``extname`` can supply an already registered map. """ z_order: ClassVar[tuple[int, ...]] = (760,) @@ -562,7 +563,7 @@ def apply_to(self, obj, **kwargs): return obj shape = obj.hdu.data.shape - if self._map is None or shape != self._map_shape: + if self._file is not None or self._map is None or shape != self._map_shape: self._map = self._make_map(shape) self._map_shape = shape @@ -570,7 +571,13 @@ def apply_to(self, obj, **kwargs): return obj def _make_map(self, shape): - if self._model is None: + if self._file is not None: + # The calibration image must be registered to this image plane. + extname = str(from_currsys(self.meta["extname"], self.cmds)) + background = self._file[extname].data + if background.shape != shape: + raise ValueError("Background calibration shape differs from the image plane") + elif self._model is None: value = from_currsys(self.meta["value"], self.cmds) background = np.full(shape, value, dtype=np.float32) else: diff --git a/scopesim/effects/metis_lms_trace_list.py b/scopesim/effects/metis_lms_trace_list.py index 92780daca..64f6f460d 100644 --- a/scopesim/effects/metis_lms_trace_list.py +++ b/scopesim/effects/metis_lms_trace_list.py @@ -144,6 +144,8 @@ def apply_to_fov(self, obj): slicefov.detector_header = obj.detector_header slicefov.meta["xi_min"] = obj.meta["xi_min"] slicefov.meta["xi_max"] = obj.meta["xi_max"] + slicefov.meta["eta_min"] = ymin * u.arcsec + slicefov.meta["eta_max"] = ymax * u.arcsec slicefov.meta["trace_id"] = sptid slicefov.cube = fits.ImageHDU(header=slicewcs.to_header(), data=slicecube) @@ -322,6 +324,7 @@ class MetisLMSSpectralTrace(SpectralTrace): """SpectralTrace for the METIS LM spectrograph.""" _class_params = { + **SpectralTrace._class_params, "naxis1": 122, "nslice": 28, "slicewidth": 0.0207, # arcsec @@ -393,7 +396,7 @@ def get_waverange(self, det_mm_lims): lam0 = from_currsys(self.meta["wavelen"], self.cmds) xi0 = 0. - ymid = self.xilam2y(xi0, lam0)[0] # estimate y level of trace + ymid = self.xilam2y(xi0, lam0) # scalar input gives a scalar trace coordinate waverange = self.xy2lam(np.array([xmin, xmax]), np.array([ymid, ymid]), grid=False) diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py index d5635bc71..3e7fbc71b 100644 --- a/scopesim/effects/spectral_trace_list_utils.py +++ b/scopesim/effects/spectral_trace_list_utils.py @@ -14,8 +14,8 @@ import numpy as np -from scipy.interpolate import RectBivariateSpline -from scipy.interpolate import interp1d +from scipy.interpolate import interp1d, RectBivariateSpline +from scipy.ndimage import gaussian_filter1d, map_coordinates from astropy.table import Table, vstack from astropy.modeling import fitting @@ -103,6 +103,7 @@ class SpectralTrace: "pixel_size": None, "description": "", "trace_flux_jacobian": "!SIM.spectral.trace_flux_jacobian", + "fit_degree": 4, } def __init__(self, trace_tbl, cmds=None, **kwargs): @@ -128,6 +129,9 @@ def __init__(self, trace_tbl, cmds=None, **kwargs): ("PLTSCALE", "plate_scale"), ("DETPAD", "detector_pad"), ("DETANG", "detector_angle"), + ("FITDEG", "fit_degree"), + ("LSFRES", "optical_resolving_power"), + ("SLITCOL", "slit_dispersion_colname"), ): if header_key in trace_tbl.header: self.meta[meta_key] = trace_tbl.header[header_key] @@ -194,12 +198,19 @@ def compute_interpolation_functions(self): self.wave_min = quantify(np.min(lam_arr), u.um).value self.wave_max = quantify(np.max(lam_arr), u.um).value - self.xy2xi = Transform2D.fit(x_arr, y_arr, xi_arr) - self.xy2lam = Transform2D.fit(x_arr, y_arr, lam_arr) - self.xilam2x = Transform2D.fit(xi_arr, lam_arr, x_arr) - self.xilam2y = Transform2D.fit(xi_arr, lam_arr, y_arr) - self._xiy2x = Transform2D.fit(xi_arr, y_arr, x_arr) - self._xiy2lam = Transform2D.fit(xi_arr, y_arr, lam_arr) + degree = self.meta["fit_degree"] + self.xy2xi = Transform2D.fit(x_arr, y_arr, xi_arr, degree=degree) + self.xy2lam = Transform2D.fit(x_arr, y_arr, lam_arr, degree=degree) + self.xilam2x = Transform2D.fit(xi_arr, lam_arr, x_arr, degree=degree) + self.xilam2y = Transform2D.fit(xi_arr, lam_arr, y_arr, degree=degree) + self._xiy2x = Transform2D.fit(xi_arr, y_arr, x_arr, degree=degree) + self._xiy2lam = Transform2D.fit(xi_arr, y_arr, lam_arr, degree=degree) + if "slit_dispersion_colname" in self.meta: + center = xi_arr == 0 + derivative = self.table[self.meta["slit_dispersion_colname"]].quantity.to_value(u.um / u.arcsec) + # FOVs are restricted to this trace; allow unit-conversion roundoff + # at the wavelength endpoints when interpolating the derivative. + self.slit_dispersion = interp1d(lam_arr[center], derivative[center], bounds_error=False, fill_value="extrapolate") if self.dispersion_axis == "unknown": dlam_dx, dlam_dy = self.xy2lam.gradient() @@ -261,6 +272,24 @@ def _focal_plane_grid(self, fov): y_mm = np.linspace(ymin_mm, ymax_mm, sub_naxis2, endpoint=False, dtype=np.float32) return xmin, ymin, x_mm, y_mm, det_wcs + def spectral_padding(self, fov): + """Input wavelength support for slit displacement and optical blur.""" + lower, upper = fov.meta['wave_min'].to_value(u.um), fov.meta['wave_max'].to_value(u.um) + wave = np.unique(u.Quantity(self.table[self.meta['wave_colname']]).to_value(u.um)) + wave = np.r_[lower, wave[(wave > lower) & (wave < upper)], upper] + width = (fov.meta['eta_max'] - fov.meta['eta_min']).to_value(u.arcsec) + if 'slit_dispersion_colname' in self.meta: + derivative = self.slit_dispersion(wave) + else: + pixsize = (fov.detector_header['CDELT1D']*u.Unit(fov.detector_header['CUNIT1D'])).to_value(u.mm) + self._set_dispersion(lower, upper, pixsize=pixsize) + step = abs(fov.header['CDELT2'])*3600 + derivative = self.dlam_per_pix(wave)/step + padding = np.max(abs(derivative))*width/2 + if 'optical_resolving_power' in self.meta: + padding += 5*upper/(self.meta['optical_resolving_power']*np.sqrt(8*np.log(2))) + return padding*u.um + def map_spectra_to_focal_plane(self, fov): """ Apply the spectral trace mapping to a spectral cube. @@ -294,14 +323,13 @@ def map_spectra_to_focal_plane(self, fov): image = np.zeros((sub_naxis2, sub_naxis1), dtype=np.float32) self._set_dispersion(wave_min, wave_max, pixsize=pixsize) - try: - xilam = XiLamImage(fov, self.dlam_per_pix) - self._xilamimg = xilam # ..todo: remove or make available with a debug flag? - except ValueError: - logger.warning(" ---> %s gave ValueError", self.trace_id) - - npix_xi, npix_lam = xilam.npix_xi, xilam.npix_lam - xilam_wcs = xilam.wcs + slit_dispersion = None + if "slit_dispersion_colname" in self.meta: + slit_dispersion = self.slit_dispersion + xilam = XiLamImage(fov, self.dlam_per_pix, slit_dispersion=slit_dispersion) + if "optical_resolving_power" in self.meta: + xilam.convolve_lsf(self.meta["optical_resolving_power"]) + self._xilamimg = xilam # focal-plane coordinate images ximg_fpa, yimg_fpa = np.meshgrid(x_mm, y_mm) @@ -310,33 +338,11 @@ def map_spectra_to_focal_plane(self, fov): xi_fpa = self.xy2xi(ximg_fpa, yimg_fpa).astype(np.float32) lam_fpa = self.xy2lam(ximg_fpa, yimg_fpa).astype(np.float32) - # mask everything outside the wavelength range - mask = (xi_fpa >= xi_min) & (xi_fpa <= xi_max) - xi_fpa *= mask - lam_fpa *= mask - - # mask everything outside the wavelength range - # mask = (lam_fpa >= wave_min) & (lam_fpa <= wave_max) - # xi_fpa *= mask - # lam_fpa *= mask - - # Convert to pixel images - # These are the pixel coordinates in the image corresponding to - # xi, lambda - # It is much quicker to do the linear transformation by hand - # than to use the astropy.wcs functions for conversion. - i_img = ((lam_fpa - xilam_wcs.wcs.crval[0]) - / xilam_wcs.wcs.cdelt[0]).astype(int) - j_img = ((xi_fpa - xilam_wcs.wcs.crval[1]) - / xilam_wcs.wcs.cdelt[1]).astype(int) - - # truncate images to remove pixel coordinates outside the image - ijmask = ((i_img >= 0) * (i_img < npix_lam) - * (j_img >= 0) * (j_img < npix_xi)) - - # do the actual interpolation - # image is in [ph/s/um/arcsec] - image = xilam.interp(xi_fpa, lam_fpa, grid=False) * ijmask + # The physical slit ends include the outer half-cells. A truncated + # integer index incorrectly admits pixels below the lower boundary. + inside = (xi_fpa >= xi_min) & (xi_fpa < xi_max) + inside &= (lam_fpa >= wave_min) & (lam_fpa < wave_max) + image = xilam.interp(xi_fpa, lam_fpa, grid=False)*inside # Scale to ph / s / pixel image *= self._trace_flux_scale( @@ -380,11 +386,7 @@ def map_spectra_to_focal_plane(self, fov): def _trace_flux_scale(self, x_mm, y_mm, pixsize, pixscale, det_header): """Return the local [arcsec um] per detector pixel flux scale.""" - try: - use_jacobian = bool(from_currsys( - self.meta["trace_flux_jacobian"], self.cmds)) - except ValueError: - use_jacobian = False + use_jacobian = from_currsys(self.meta["trace_flux_jacobian"], self.cmds) if not use_jacobian: dlam_by_dx, dlam_by_dy = self.xy2lam.gradient() @@ -401,13 +403,10 @@ def _trace_flux_scale(self, x_mm, y_mm, pixsize, pixscale, det_header): @staticmethod def _detector_pixel_area(det_header): """Return detector WCS pixel area in mm2.""" - try: - pixel_scale_matrix = WCS(det_header, key="D").pixel_scale_matrix - return np.abs(np.linalg.det(pixel_scale_matrix)) - except Exception: # pragma: no cover - FITS fallback for odd headers - cdelt1 = det_header["CDELT1D"] * u.Unit(det_header["CUNIT1D"]) - cdelt2 = det_header["CDELT2D"] * u.Unit(det_header["CUNIT2D"]) - return np.abs((cdelt1 * cdelt2).to_value(u.mm**2)) + detector_wcs = WCS(det_header, key="D") + pixel_area = np.abs(np.linalg.det(detector_wcs.pixel_scale_matrix)) + pixel_area *= u.Unit(detector_wcs.wcs.cunit[0]) * u.Unit(detector_wcs.wcs.cunit[1]) + return pixel_area.to_value(u.mm**2) def rectify( self, @@ -851,9 +850,13 @@ class XiLamImage(): dlam_per_pix : a 1-D interpolation function from wavelength (in um) to dispersion (in um/pixel); alternatively a number giving an average dispersion + slit_dispersion : callable, optional + Physical cross-slit derivative d(lambda)/d(eta) in um/arcsec as a + function of wavelength in um. When supplied, this replaces the + legacy coupling between sky-grid sampling and detector dispersion. """ - def __init__(self, fov, dlam_per_pix): + def __init__(self, fov, dlam_per_pix, slit_dispersion=None): # ..todo: we assume that we always have a cube. We use SpecCADO's # add_cube_layer method cube_wcs = WCS(fov.cube.header, key=" ") @@ -871,8 +874,16 @@ def __init__(self, fov, dlam_per_pix): (n_lam, n_eta, n_xi) = fov.cube.data.shape # arrays of cube coordinates - cube_xi = d_xi * np.arange(n_xi) + fov.meta["xi_min"].value - cube_eta = d_eta * (np.arange(n_eta) - (n_eta - 1) / 2) + cube_xi = cube_wcs.all_pix2world(np.arange(n_xi), np.zeros(n_xi), np.zeros(n_xi), 0)[0] + cube_eta = cube_wcs.all_pix2world(np.zeros(n_eta), np.arange(n_eta), np.zeros(n_eta), 0)[1] + if cube_wcs.has_celestial: + # Slit coordinates are offsets about zero; RA wraps at 360 degrees. + cube_xi = (cube_xi + 180) % 360 - 180 + cube_xi *= u.Unit(cube_wcs.wcs.cunit[0]).to(u.arcsec) + cube_eta *= u.Unit(cube_wcs.wcs.cunit[1]).to(u.arcsec) + # Each IFU slice has its own optical axis, independent of sky position. + eta_center = (fov.meta["eta_min"] + fov.meta["eta_max"]) / 2 + cube_eta -= eta_center.to_value(u.arcsec) cube_lam = wcs_lam.all_pix2world(np.arange(n_lam), 0)[0] cube_lam *= u.Unit(wcs_lam.wcs.cunit[0]).to(u.um) @@ -886,24 +897,31 @@ def __init__(self, fov, dlam_per_pix): logger.warning("Using scalar dlam_per_pix = %.2g", dlam_per_pix_val) + dlam_per_arcsec = dlam_per_pix_val / d_eta if slit_dispersion is None else slit_dispersion(self.lam) for i, eta in enumerate(cube_eta): - lam0 = self.lam + dlam_per_pix_val * eta / d_eta + # A positive calibrated derivative moves an input line to a larger + # apparent wavelength: sample the source at lambda_out - shift. + lam0 = self.lam + dlam_per_arcsec * eta if slit_dispersion is None else self.lam - dlam_per_arcsec * eta # lam0 is the target wavelength. We need to check that this # overlaps with the wavelength range covered by the cube if lam0.min() < cube_lam.max() and lam0.max() > cube_lam.min(): - plane = fov.cube.data[:, i, :].T - plane_interp = RectBivariateSpline(cube_xi, cube_lam, plane, - kx=1, ky=1) - self.image += plane_interp(cube_xi, lam0) + # Only wavelength changes here; xi is already on its native + # grid. Building a 2D spline for every slit cell is redundant. + index = np.searchsorted(cube_lam, lam0, side='right') - 1 + index = np.clip(index, 0, n_lam - 2) + fraction = (lam0 - cube_lam[index])/d_lam + valid = (lam0 >= cube_lam[0]) & (lam0 <= cube_lam[-1]) + plane = fov.cube.data[:, i, :] + self.image += ((1 - fraction[:, None])*plane[index] + fraction[:, None]*plane[index + 1]).T*valid - self.image *= d_eta # ph/s/um/arcsec2 --> ph/s/um/arcsec + self.image *= abs(d_eta) # ph/s/um/arcsec2 --> ph/s/um/arcsec # WCS for the xi-lambda image, i.e. the rectified 2D spectrum # Default WCS with xi in arcsec self.wcs = WCS(naxis=2) self.wcs.wcs.crpix = [1, 1] - self.wcs.wcs.crval = [self.lam[0], fov.meta["xi_min"].value] + self.wcs.wcs.crval = [self.lam[0], cube_xi[0]] self.wcs.wcs.pc = [[1, 0], [0, 1]] self.wcs.wcs.cdelt = [d_lam, d_xi] self.wcs.wcs.ctype = ["LINEAR", "LINEAR"] @@ -923,14 +941,30 @@ def __init__(self, fov, dlam_per_pix): self.xi = self.wcs.all_pix2world(self.lam[0], np.arange(n_xi), 0)[1] self.npix_xi = n_xi self.npix_lam = n_lam - # ..todo: cubic spline introduces negative values, linear does not. - # Alternative might be to cubic-spline interpolate on sqrt(image), - # with subsequent squaring of the result. This would require - # wrapping RectBivariateSpline in a new (sub)class. - spline_order = (1, 1) - self.interp = RectBivariateSpline(self.xi, self.lam, self.image, - kx=spline_order[0], - ky=spline_order[1]) + + def interp(self, xi, wavelength, grid=True): + """Bilinear density interpolation on the existing uniform WCS grid.""" + xi, wavelength = np.meshgrid(xi, wavelength, indexing='ij') if grid else np.broadcast_arrays(xi, wavelength) + row = (xi - self.xi[0])/self.wcs.wcs.cdelt[1] + column = (wavelength - self.lam[0])/self.wcs.wcs.cdelt[0] + return map_coordinates(self.image, [row, column], order=1, mode='nearest', prefilter=False) + + def convolve_lsf(self, resolving_power): + """Apply an optical Gaussian LSF at constant R, preserving photon density. + + This is the optical blur alone; slit illumination is integrated above. + Convolving lambda * I_lambda on a uniform log(lambda) grid preserves + integrated flux rather than treating I_lambda as a log-bin count. + """ + if resolving_power <= 0: + raise ValueError("Optical resolving power must be positive") + log_wave = np.log(self.lam) + n = int(np.ceil((log_wave[-1] - log_wave[0]) / np.min(np.diff(log_wave)))) + 1 + log_grid = np.linspace(log_wave[0], log_wave[-1], n) + sigma = 1 / (resolving_power * np.sqrt(8 * np.log(2)) * (log_grid[1] - log_grid[0])) + density = interp1d(log_wave, self.image * self.lam, axis=1)(log_grid) + blurred = gaussian_filter1d(density, sigma, axis=1, mode="nearest") + self.image = interp1d(log_grid, blurred, axis=1)(log_wave) / self.lam class Transform2D(): diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py index 6efd90f7b..c0d68eb2e 100644 --- a/scopesim/effects/ter_curves.py +++ b/scopesim/effects/ter_curves.py @@ -547,9 +547,14 @@ def load_skycalc_table(self): else: path = find_file(use_local_file) - fits_tbl = fits.getdata(path, ext=1) + with fits.open(path) as sky_file: + fits_tbl = sky_file[1].data + # Older cached files require an explicit unit if TUNIT is absent. + unit = self.meta["local_skycalc_wavelength_unit"] if "local_skycalc_wavelength_unit" in self.meta else sky_file[1].columns["lam"].unit + wave_unit = u.Unit(unit) fits_hdr = fits.getheader(path, ext=0) tbl = Table(fits_tbl) + tbl["lam"] = tbl["lam"] * wave_unit.to(u.um) tbl["lam"].unit = u.um for colname in tbl.colnames: if "flux" in colname: diff --git a/scopesim/tests/test_utils_functions.py b/scopesim/tests/test_utils_functions.py index a1df0f1b3..1469b32d9 100644 --- a/scopesim/tests/test_utils_functions.py +++ b/scopesim/tests/test_utils_functions.py @@ -199,6 +199,15 @@ def test_converts_list(self): def test_converts_numpy_array(self): assert utils.from_currsys(np.array(["!SIM.random.seed"]*2))[1] is None + @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32]) + def test_numeric_calibration_array_keeps_values_and_precision(self, dtype): + values = np.arange(12, dtype=dtype).reshape(3, 4) + resolved = utils.from_currsys(values) + np.testing.assert_array_equal(resolved, values) + assert resolved.dtype == dtype + resolved[0, 0] = 20 + assert values[0, 0] == 0 + def test_converts_dict(self): assert utils.from_currsys({"seed": "!SIM.random.seed"})["seed"] is None diff --git a/scopesim/tests/tests_effects/test_SkycalcTERCurve.py b/scopesim/tests/tests_effects/test_SkycalcTERCurve.py index 689621b10..b87995c81 100644 --- a/scopesim/tests/tests_effects/test_SkycalcTERCurve.py +++ b/scopesim/tests/tests_effects/test_SkycalcTERCurve.py @@ -59,7 +59,9 @@ def test_initialises_with_non_skycalc_keys(self): def test_initialise_with_local_skycalc_file(self, mock_path): sky_ter = SkycalcTERCurve( - use_local_skycalc_file=str(mock_path / "skycalc_override.fits")) + use_local_skycalc_file=str(mock_path / "skycalc_override.fits"), + local_skycalc_wavelength_unit=u.um) + # This legacy fixture predates wavelength units in its FITS columns. assert sky_ter.skycalc_table is not None diff --git a/scopesim/tests/tests_effects/test_SpectralTraceList.py b/scopesim/tests/tests_effects/test_SpectralTraceList.py index 8461857a2..5ff0633aa 100644 --- a/scopesim/tests/tests_effects/test_SpectralTraceList.py +++ b/scopesim/tests/tests_effects/test_SpectralTraceList.py @@ -6,6 +6,7 @@ import numpy as np from astropy import units as u from astropy.io import fits +from astropy.table import Table from scopesim.effects.spectral_trace_list import SpectralTraceList, \ @@ -137,6 +138,93 @@ def _echelle_cmds(): class TestEchelleSpectralTraceList: + @pytest.mark.parametrize('angle', [-25, 0, 25]) + @pytest.mark.parametrize('offset', [(0., 0.), (1.5, -1.2)]) + def test_fixed_placement_and_detector_intersections(self, tmp_path, angle, offset): + params_file = tmp_path / 'echelle_trace_parameters.dat' + _write_echelle_trace_params(params_file, detector_angle=angle) + params = Table.read(params_file, format='ascii') + for key, value in {'m0': 75, 'n': 6, 'min_wave': 960, 'max_wave': 1340, 'n_disp': 2048, 'n_xdisp': 1024, 'disp_freq': 24., 'xdisp_freq': 350.}.items(): + params[key][0] = value + params['xbeta_center'] = [10.9] + params['detector_offset_x'], params['detector_offset_y'] = [offset[0]], [offset[1]] + params.meta['comments'].extend(['detector_offset_x_unit : mm', 'detector_offset_y_unit : mm']) + params.write(params_file, format='ascii', overwrite=True) + effect = EchelleSpectralTraceList(cmds=_echelle_cmds(), filename=str(params_file), wave_colname='wavelength', s_colname='s') + model = effect._spectrographs['b'] + nx, ny = model.detector.n_pix_x, model.detector.n_pix_y + raw_grid = np.linspace(-2*nx, 3*nx, 50001) + for order, edges in zip(model.orders, model.edge_wave(fsr=True)): + wave = np.linspace(edges[0], edges[1], 4*nx + 1) + x, y = model._detector_pixels(wave, order) + # Rotate about the camera origin, then translate in detector axes. + xraw = model.wavelength_to_x_pixel(wave, order).value - nx/2 + yraw = model.wavelength_to_y_pixel(wave).value[None, :] - ny/2 + model._slit_offset_pix[:, None] + c, s = np.cos(np.deg2rad(angle)), np.sin(np.deg2rad(angle)) + expected_x = c*xraw - s*yraw + nx/2 + offset[0]/.015 + expected_y = s*xraw + c*yraw + ny/2 + offset[1]/.015 + np.testing.assert_allclose(x, expected_x, atol=1e-9, rtol=0) + np.testing.assert_allclose(y, expected_y, atol=1e-9, rtol=0) + # Translation changes placement, not dispersion or slit length. + xrot, yrot = model._rotated_detector_pixels(wave, order) + np.testing.assert_allclose(np.diff(x), np.diff(xrot), atol=1e-9, rtol=0) + np.testing.assert_allclose(np.diff(y), np.diff(yrot), atol=1e-9, rtol=0) + slit_length = np.hypot(x[2] - x[0], y[2] - y[0]) + np.testing.assert_allclose(slit_length, model._slit_offset_pix[2] - model._slit_offset_pix[0], atol=1e-9, rtol=0) + + # A deliberately oversized wavelength search independently checks + # that the translated clipping search has not lost either edge. + wave = model.x_pixel_to_wavelength(raw_grid, order) + x, y = model._detector_pixels(wave, order) + inside = np.any((x >= 0) & (x <= nx) & (y >= 0) & (y <= ny), axis=0) + trace_id = f'b_{order}' + if not inside.any(): + assert trace_id not in effect.spectral_traces + continue + trace = effect.spectral_traces[trace_id] + actual = model.wavelength_to_x_pixel(u.Quantity(trace.table['wavelength']), order).value + expected = raw_grid[inside][[0, -1]] + # The generator locates each edge on a roughly one-pixel grid. + np.testing.assert_allclose([actual.min(), actual.max()], expected, atol=3, rtol=0) + + @pytest.mark.parametrize('missing', ['detector_offset_x', 'detector_offset_y', 'detector_offset_x_unit', 'detector_offset_y_unit']) + def test_incomplete_placement_fails(self, tmp_path, missing): + path = tmp_path / 'echelle_trace_parameters.dat' + _write_echelle_trace_params(path) + params = Table.read(path, format='ascii') + for axis in ['x', 'y']: + key = f'detector_offset_{axis}' + if missing != key: + params[key] = [0.] + if missing != f'{key}_unit': + params.meta['comments'].append(f'{key}_unit : mm') + params.write(path, format='ascii', overwrite=True) + with pytest.raises(KeyError, match=missing): + EchelleSpectralTraceList(cmds=_echelle_cmds(), filename=str(path)) + + @pytest.mark.parametrize('missing_value', ['detector_offset_x', 'detector_offset_y']) + def test_blank_placement_value_does_not_become_zero(self, tmp_path, missing_value): + path = tmp_path / 'echelle_trace_parameters.dat' + _write_echelle_trace_params(path) + params = Table.read(path, format='ascii') + for axis in ['x', 'y']: + key = f'detector_offset_{axis}' + params[key] = np.ma.array([0.], mask=[key == missing_value]) + params.meta['comments'].append(f'{key}_unit : mm') + params.write(path, format='ascii', overwrite=True) + with pytest.raises(ValueError, match='Both detector offsets must have values'): + EchelleSpectralTraceList(cmds=_echelle_cmds(), filename=str(path)) + + def test_placement_requires_length_units(self, tmp_path): + path = tmp_path / 'echelle_trace_parameters.dat' + _write_echelle_trace_params(path) + params = Table.read(path, format='ascii') + params['detector_offset_x'], params['detector_offset_y'] = [0.], [0.] + params.meta['comments'].extend(['detector_offset_x_unit : arcsec', 'detector_offset_y_unit : mm']) + params.write(path, format='ascii', overwrite=True) + with pytest.raises(u.UnitConversionError): + EchelleSpectralTraceList(cmds=_echelle_cmds(), filename=str(path)) + def test_does_not_write_generated_hdulist_by_default(self, tmp_path): params_file = tmp_path / "echelle_trace_parameters.dat" _write_echelle_trace_params(params_file) @@ -197,6 +285,13 @@ def test_detector_angle_clips_to_detector_not_padding(self, tmp_path): assert trace.meta["detector_angle"] == 25 assert trace.meta["detector_pad"] == 10 + model = spt._spectrographs['b'] + order = int(trace.meta['trace_id'].split('_')[1]) + wave = u.Quantity(trace.table['wavelength'][:n_wave]) + model_x, model_y = model._detector_pixels(wave, order) + np.testing.assert_allclose(model_x, xpix, atol=1e-8, rtol=0) + np.testing.assert_allclose(model_y, ypix, atol=1e-8, rtol=0) + padded_file = tmp_path / "echelle_trace_parameters_more_padding.dat" _write_echelle_trace_params( padded_file, @@ -224,7 +319,9 @@ def test_detector_angle_clips_to_detector_not_padding(self, tmp_path): ) unrotated_trace = next(iter(unrotated.spectral_traces.values())) - assert trace.wave_min == pytest.approx(unrotated_trace.wave_min) - assert trace.wave_max == pytest.approx(unrotated_trace.wave_max) + # Rotation changes which wavelengths reach the rectangular detector; + # clipping must occur after rotation, not at the old unrotated edges. + assert trace.wave_min < trace.wave_max + assert not np.allclose([trace.wave_min, trace.wave_max], [unrotated_trace.wave_min, unrotated_trace.wave_max]) assert not np.allclose(trace.table["x"], unrotated_trace.table["x"]) assert not np.allclose(trace.table["y"], unrotated_trace.table["y"]) diff --git a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py index dfaa2d7f6..04be6beef 100644 --- a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py +++ b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py @@ -196,7 +196,7 @@ def test_call_gives_correct_result(self, quadratic): # matrix and explicit function td2d = Transform2D(quadratic['matrix']) - assert td2d(x, y) == quadratic['function'](x, y) + assert td2d(x, y) == pytest.approx(quadratic['function'](x, y)) def test_gradient_gives_correct_result(self, quadratic): x = np.random.randn() @@ -205,8 +205,8 @@ def test_gradient_gives_correct_result(self, quadratic): tf2d = Transform2D(quadratic['matrix']) tf2d_grad = tf2d.gradient() - assert tf2d_grad[0](x, y) == quadratic['gradient'][0](x, y) - assert tf2d_grad[1](x, y) == quadratic['gradient'][1](x, y) + assert tf2d_grad[0](x, y) == pytest.approx(quadratic['gradient'][0](x, y)) + assert tf2d_grad[1](x, y) == pytest.approx(quadratic['gradient'][1](x, y)) def test_grid_true_gives_correct_shape(self, tf2d): n_x, n_y = 12, 3 @@ -303,7 +303,7 @@ def __init__(self, n_lam=20, n_eta=3, n_xi=11): hdr["CUNIT1"], hdr["CUNIT2"], hdr["CUNIT3"] = "arcsec", "arcsec", "um" self.cube = fits.ImageHDU( data=np.ones((n_lam, n_eta, n_xi)), header=hdr) - self.meta = {"xi_min": -0.5 * u.arcsec, "xi_max": 0.5 * u.arcsec} + self.meta = {"xi_min": -0.5 * u.arcsec, "xi_max": 0.5 * u.arcsec, "eta_min": -0.5 * u.arcsec, "eta_max": 0.5 * u.arcsec} class TestXiLamImage: diff --git a/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py b/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py new file mode 100644 index 000000000..00114bcee --- /dev/null +++ b/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py @@ -0,0 +1,87 @@ +from types import SimpleNamespace +import numpy as np +import pytest +from astropy import units as u +from astropy.io import fits +from scopesim.effects.spectral_trace_list_utils import XiLamImage, SpectralTrace + + +def line_cube(step, offset=0): + wave = np.linspace(.995, 1.005, 2001) + eta = np.arange(-.5+step/2, .5, step) + spatial = np.exp(-.5*((eta-offset)/.12)**2) + profile = np.exp(-.5*((wave-1)/.00002)**2) + cube = profile[:, None, None]*spatial[None, :, None]*np.ones((1, 1, 3)) + header = fits.Header({'CTYPE1': 'LINEAR', 'CTYPE2': 'LINEAR', 'CTYPE3': 'WAVE', 'CUNIT1': 'arcsec', 'CUNIT2': 'arcsec', 'CUNIT3': 'um', 'CDELT1': .1, 'CDELT2': step, 'CDELT3': wave[1]-wave[0], 'CRPIX1': 1., 'CRPIX2': 1., 'CRPIX3': 1., 'CRVAL1': 0., 'CRVAL2': eta[0], 'CRVAL3': wave[0]}) + return SimpleNamespace(cube=fits.ImageHDU(cube, header), meta={'xi_min': 0*u.arcsec, 'eta_min': -.5*u.arcsec, 'eta_max': .5*u.arcsec}) + + + +def test_slit_dispersion_is_independent_of_numerical_sky_grid(): + spectra = [XiLamImage(line_cube(step), lambda w: .00001, slit_dispersion=lambda w: .0002) for step in [.1, .05, .025]] + # Finite spatial quadrature changes, while the physical slit dispersion stays fixed. + for image in spectra[:-1]: + np.testing.assert_allclose(image.image, spectra[-1].image, atol=.001, rtol=.02) + + + +def test_calibrated_slit_shift_has_forward_model_sign(): + image = XiLamImage(line_cube(.025, offset=.1), lambda w: .00001, slit_dispersion=lambda w: .002) + centroid = np.sum(image.image[1]*image.lam)/image.image[1].sum() + assert centroid == pytest.approx(1.0002, abs=2e-6) + + + +def test_optical_lsf_preserves_photons_and_adds_variance(): + image = XiLamImage(line_cube(.05), lambda w: 0., slit_dispersion=lambda w: 0.) + before = np.trapezoid(image.image[1], image.lam) + image.convolve_lsf(10000) + after = np.trapezoid(image.image[1], image.lam) + assert after == pytest.approx(before, rel=.001) + center = np.sum(image.lam*image.image[1])/image.image[1].sum() + variance = np.sum((image.lam-center)**2*image.image[1])/image.image[1].sum() + assert variance == pytest.approx(.00002**2 + (1/10000/np.sqrt(8*np.log(2)))**2, rel=.02) + + +def test_slit_coordinates_follow_cube_wcs(): + fov = line_cube(.025) + fov.cube.header['CRVAL1'] = 2. + fov.cube.header['CRPIX1'] = 2. + fov.cube.header['CRVAL2'] += .1 + image = XiLamImage(fov, lambda wave: .00001, slit_dispersion=lambda wave: .002) + np.testing.assert_allclose(image.xi, [1.9, 2., 2.1]) + centroid = np.sum(image.image[1] * image.lam) / image.image[1].sum() + assert centroid == pytest.approx(1.0002, abs=2e-6) + + +def test_translating_an_ifu_slice_preserves_its_local_spectrum(): + fov = line_cube(.025) + centered = XiLamImage(fov, lambda wave: .00001, slit_dispersion=lambda wave: .002) + fov.cube.header['CRVAL2'] += 2. + fov.meta['eta_min'] += 2*u.arcsec + fov.meta['eta_max'] += 2*u.arcsec + translated = XiLamImage(fov, lambda wave: .00001, slit_dispersion=lambda wave: .002) + np.testing.assert_allclose(translated.image, centered.image, atol=1e-7) + + +def test_missing_aperture_geometry_fails(): + fov = line_cube(.025) + del fov.meta['eta_min'] + with pytest.raises(KeyError, match='eta_min'): + XiLamImage(fov, lambda wave: .00001, slit_dispersion=lambda wave: .002) + + +def test_detector_jacobian_uses_physical_pixel_area(): + header = fits.Header({'CTYPE1D': 'LINEAR', 'CTYPE2D': 'LINEAR', 'CUNIT1D': 'm', 'CUNIT2D': 'm', 'CDELT1D': .000015, 'CDELT2D': .000030, 'CRPIX1D': 1., 'CRPIX2D': 1., 'CRVAL1D': 0., 'CRVAL2D': 0.}) + assert SpectralTrace._detector_pixel_area(header) == pytest.approx(.015*.030) + + +def test_rectified_interpolation_reproduces_a_linear_density_field(): + fov = line_cube(.1) + wave = .995 + np.arange(2001)*fov.cube.header['CDELT3'] + xi = np.arange(3)*.1 + fov.cube.data[:] = 2 + 3*xi[None, None, :] + 5*wave[:, None, None] + image = XiLamImage(fov, lambda w: 0., slit_dispersion=lambda w: 0.) + positions, wavelengths = np.array([.03, .17]), np.array([.998, 1.003]) + expected = 2 + 3*positions + 5*wavelengths + np.testing.assert_allclose(image.interp(positions, wavelengths, grid=False), expected, rtol=2e-6) diff --git a/scopesim/tests/tests_effects/test_local_skycalc.py b/scopesim/tests/tests_effects/test_local_skycalc.py new file mode 100644 index 000000000..8baab5ba1 --- /dev/null +++ b/scopesim/tests/tests_effects/test_local_skycalc.py @@ -0,0 +1,24 @@ +import numpy as np +import pytest +from astropy import units as u +from astropy.table import Table +from scopesim.effects import SkycalcTERCurve + + +@pytest.mark.parametrize('wave_unit', [u.nm, u.um]) +def test_local_skycalc_respects_wavelength_unit(tmp_path, wave_unit): + wave = np.array([.5, 1., 2.]) * u.um + table = Table([wave.to(wave_unit), [.5, .8, .9], [100., 200., 300.]], names=['lam', 'trans', 'flux']) + path = tmp_path / 'sky.fits' + table.write(path) + sky = SkycalcTERCurve(use_local_skycalc_file=str(path)) + np.testing.assert_allclose(sky.skycalc_table['wavelength'], wave.value) + assert sky.surface.transmission(1*u.um).value == pytest.approx(.8) + + + +def test_skycalc_missing_wavelength_unit_fails(tmp_path): + path = tmp_path / 'sky.fits' + Table([[500., 1000., 2000.], [.5, .8, .9], [100., 200., 300.]], names=['lam', 'trans', 'flux']).write(path) + with pytest.raises(TypeError, match='Unit'): + SkycalcTERCurve(use_local_skycalc_file=str(path)) diff --git a/scopesim/tests/tests_effects/test_registered_background.py b/scopesim/tests/tests_effects/test_registered_background.py new file mode 100644 index 000000000..54606b9a5 --- /dev/null +++ b/scopesim/tests/tests_effects/test_registered_background.py @@ -0,0 +1,32 @@ +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.effects import ImagePlaneBackground +from scopesim.optics.image_plane import ImagePlane + + +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_background_uses_current_named_calibration(tmp_path): + path = tmp_path / 'background.fits' + fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(np.full((80, 120), .15), name='0.9'), fits.ImageHDU(np.full((80, 120), .04), name='0.9JH')]).writeto(path) + cmds = {'!INST.slit': '0.9'} + effect = ImagePlaneBackground(filename=str(path), extname='!INST.slit', cmds=cmds) + plane = ImagePlane(detector().header) + effect.apply_to(plane) + np.testing.assert_allclose(plane.data, .15) + cmds['!INST.slit'] = '0.9JH' + plane.data[:] = 0 + effect.apply_to(plane) + np.testing.assert_allclose(plane.data, .04) + with pytest.raises(ValueError, match='shape'): + effect.apply_to(ImagePlane(detector(60, 40).header)) + cmds['!INST.slit'] = 'missing' + with pytest.raises(KeyError): + effect.apply_to(plane) diff --git a/scopesim/utils.py b/scopesim/utils.py index 36b2e5080..a6c60801d 100644 --- a/scopesim/utils.py +++ b/scopesim/utils.py @@ -610,6 +610,9 @@ def from_currsys(item, cmds=None): item.meta = item_meta if isinstance(item, np.ndarray) and not isinstance(item, u.Quantity): + # Numeric calibration arrays cannot contain command references. + if item.dtype.kind in "biufc": + return item.copy() item = np.array([from_currsys(x, cmds) for x in item]) if isinstance(item, list): @@ -1261,4 +1264,4 @@ def _get_next_moon_cached(moontype="full", today=None, location_key=None): nextm = times[mask][np.argmax(moonalts[mask])] else: raise ValueError(f"Invalid moon type: {moontype}, should be 'full', 'new' or 'half'.") - return nextm \ No newline at end of file + return nextm From 008b9761ab911fc2417b9c151a309e7ecec36a46 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:39:12 -0700 Subject: [PATCH 3/7] 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 f9d727383dfc8ae1745c635b6d899356ad706a24 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:41:42 -0700 Subject: [PATCH 4/7] Preserve native spectral photon integrals and source representations Source commits: a23b80c8, cafce81d, 8a43ec43 (source precision), a0495f72 (mutable spectrum sampling) --- scopesim/optics/fov.py | 90 +++++--- scopesim/optics/optical_train.py | 54 +---- .../tests/tests_optics/test_FieldOfView.py | 208 ++++-------------- .../test_spectral_bin_integration.py | 78 +++++++ 4 files changed, 184 insertions(+), 246 deletions(-) create mode 100644 scopesim/tests/tests_optics/test_spectral_bin_integration.py diff --git a/scopesim/optics/fov.py b/scopesim/optics/fov.py index d38f4abe5..a2b2ffa9b 100644 --- a/scopesim/optics/fov.py +++ b/scopesim/optics/fov.py @@ -4,10 +4,9 @@ from warnings import warn from itertools import chain from collections.abc import Iterable, Generator -from weakref import WeakKeyDictionary import numpy as np -from scipy.interpolate import interp1d +from scipy.integrate import cumulative_trapezoid from astropy import units as u from astropy.io import fits @@ -43,7 +42,6 @@ logger = get_logger(__name__) -_SPECTRUM_WAVESET_AA_CACHE = WeakKeyDictionary() class FieldOfView: @@ -490,10 +488,11 @@ def waveset(self): for cube in cube_fields]) _waveset = cube_fields[naxis3_max].waveset elif specfields := self._get_fields(SpectrumSourceField): - _waveset = np.concatenate([ + _waveset = np.concatenate([self.waverange << u.um] + [ spec.waveset.to(u.um) for field in specfields for spec in field.spectra.values() + if spec.waveset is not None ]) else: _waveset = self.waverange << u.um @@ -641,14 +640,10 @@ def _make_cubefields(self): * yield scaled flux to be added to canvas flux """ for field in self._get_cube_fields(): - fluxes = field.data.sum(axis=(1, 2)) - fov_waveset_fluxes = np.interp(self.waveset, field.waveset, fluxes) - - # .lower() is needed because astropy doesn't recognize PHOTLAM - field_unit = field.header.get("BUNIT", PHOTLAM).lower() - flux_scale_factor = u.Unit(field_unit).to(PHOTLAM) - - yield fov_waveset_fluxes * flux_scale_factor + fluxes = field.data.sum(axis=(1, 2), dtype=np.float64)*u.Unit(field.header['BUNIT'].lower()) + if field.is_bunit_spatially_differential: + fluxes *= field.pixel_area + yield np.interp(self.waveset, field.waveset, fluxes.to_value(PHOTLAM)) def _make_imagefields(self): """ @@ -659,7 +654,7 @@ def _make_imagefields(self): * yield spectrum multiply by sum to be added to canvas flux """ for field in self._get_image_fields(): - weight = np.sum(field.data) # Shouldn't that be 1 by convention? + weight = np.sum(field.data, dtype=np.float64) # Shouldn't that be 1 by convention? yield field.spectrum(self.waveset).value * weight def _make_tablefields(self): @@ -928,7 +923,7 @@ def _make_cubefields(self, fov_waveset): Find Cube fields. * rescale and reorient cubes - * interp1d smaller cubes with waveset + * integrate native spectral samples into the FOV wavelength bins * yield cubes to be added to cavas cube """ for field in self._get_cube_fields(): @@ -936,17 +931,9 @@ def _make_cubefields(self, fov_waveset): # Assumption is that ImageHDUs have units of PHOTLAM arcsec-2 # TODO: Deal with this bounds_error in a more elegant way - field_interp = interp1d( - field.waveset.to(u.um).value, - field.data, - axis=0, - kind="linear", - bounds_error=False, - fill_value=0, - ) - - field_data = field_interp(fov_waveset.value) + field_data = rebin_spectral_density(field.waveset.to_value(u.um), field.data, spectral_bin_edges(fov_waveset.to_value(u.um))) field_hdu = fits.ImageHDU(data=field_data, header=field.header) + field_hdu.header.update(CRPIX3=1., CRVAL3=fov_waveset[0].to_value(u.um), CDELT3=np.diff(fov_waveset[:2])[0].to_value(u.um), CUNIT3='um', CTYPE3='WAVE') yield field_hdu def _make_imagefields(self, fov_waveset, spline_order): @@ -975,7 +962,7 @@ def _make_imagefields(self, fov_waveset, spline_order): differential=True, ) - spec = field.spectrum(fov_waveset) + spec = average_spectrum(field.spectrum, fov_waveset) # 2D * 1D -> 3D field_cube = canvas_image_hdu.data[None, :, :] * spec[:, None, None] yield field_cube.value @@ -998,7 +985,7 @@ def _make_tablefields(self, fov_waveset): # if this is the case and if yes, remove the condition here!! # Fluxes (yes that's the correct plural) in PHOTLAM fluxes = { - ref: spec(fov_waveset) + ref: average_spectrum(spec, fov_waveset) for ref, spec in field.spectra.items() if ref in field.field["ref"] } @@ -1025,7 +1012,7 @@ def _make_backfields(self, fov_waveset): # Remember, cube output needs PHOTLAM arcsec-2 ! # So if SOLIDANG or BUNIT or whatever is not in arcsec-2, # the spectrum shoule be scaled accordingly! - spec = field.spectrum(fov_waveset) + spec = average_spectrum(field.spectrum, fov_waveset) yield spec[:, None, None].value def make_hdu(self): @@ -1183,12 +1170,14 @@ def replace_nans(field, cmds) -> None: "set in '!SIM.computing.nan_fill_value'.") -def extract_range_from_spectrum(spectrum, waverange): +def extract_range_from_spectrum(spectrum, waverange, _waveset=None): assert isinstance(spectrum, SourceSpectrum), ( f"spectrum must be of type synphot.SourceSpectrum: {type(spectrum)}") wave_min, wave_max = quantify(waverange, u.um).to(u.AA).value - spec_waveset = _spectrum_waveset_aa_value(spectrum) + spec_waveset = _spectrum_waveset_aa_value(spectrum) if _waveset is None else _waveset + if spec_waveset is None: + return spectrum mask = (spec_waveset > wave_min) & (spec_waveset < wave_max) # FIXME: Why did I comment this out in 2023? Seems useful to have... @@ -1211,9 +1200,40 @@ def extract_range_from_spectrum(spectrum, waverange): def _spectrum_waveset_aa_value(spectrum): - """Return cached spectrum waveset values in Angstrom.""" - spec_waveset = _SPECTRUM_WAVESET_AA_CACHE.get(spectrum) - if spec_waveset is None: - spec_waveset = spectrum.waveset.to(u.AA).value - _SPECTRUM_WAVESET_AA_CACHE[spectrum] = spec_waveset - return spec_waveset + # A spectrum can change its sampling in place, including through redshift. + waveset = spectrum.waveset + return None if waveset is None else waveset.to_value(u.AA) + + +def spectral_bin_edges(centers): + """Edges of bins whose centers are the supplied wavelength samples.""" + half_steps = np.diff(centers) / 2 + return np.r_[centers[0] - half_steps[0], centers[:-1] + half_steps, centers[-1] + half_steps[-1]] + + +def rebin_spectral_density(wave, density, edges): + """Integrate a piecewise-linear density, with zero flux beyond its support. + + Wavelength is the first axis; density may also carry spatial dimensions. + The result remains a density, averaged over each requested bin. + """ + if wave[1] < wave[0]: + wave, density = wave[::-1], density[::-1] + integral = cumulative_trapezoid(density, wave, axis=0, initial=0) + clipped = np.clip(edges, wave[0], wave[-1]) + index = np.clip(np.searchsorted(wave, clipped, side='right') - 1, 0, len(wave) - 2) + dimensions = (len(edges),) + (1,)*(density.ndim - 1) + offset = (clipped - wave[index]).reshape(dimensions) + step = (wave[index + 1] - wave[index]).reshape(dimensions) + at_edges = integral[index] + density[index]*offset + (density[index + 1] - density[index])*offset**2/(2*step) + widths = np.diff(edges).reshape((len(edges) - 1,) + (1,)*(density.ndim - 1)) + return np.diff(at_edges, axis=0) / widths + + +def average_spectrum(spectrum, wavelengths): + """Average PHOTLAM over numerical bins without discarding native line samples.""" + edges = spectral_bin_edges(wavelengths.to_value(u.AA)) + native = _spectrum_waveset_aa_value(spectrum) + wave = edges if native is None else np.union1d(edges, native[(native > edges[0]) & (native < edges[-1])]) + flux = spectrum(wave*u.AA).to_value(PHOTLAM) + return rebin_spectral_density(wave, flux, edges) * PHOTLAM diff --git a/scopesim/optics/optical_train.py b/scopesim/optics/optical_train.py index 6974bbe9f..a6031f2a7 100644 --- a/scopesim/optics/optical_train.py +++ b/scopesim/optics/optical_train.py @@ -6,7 +6,6 @@ from pathlib import Path import numpy as np -from scipy.interpolate import interp1d from astropy import units as u from astropy.table import QTable from astropy.wcs import WCS @@ -14,7 +13,6 @@ from tqdm.auto import tqdm -from synphot import SourceSpectrum, Empirical1D from synphot.units import PHOTLAM from astar_utils.nested_mapping import RecursiveNestedMapping, is_bangkey, recursive_update @@ -27,7 +25,6 @@ from ..effects import ExtraFitsKeywords from ..utils import from_currsys, top_level_catch, get_logger from ..source.source_templates import empty_sky -from ..source.source_fields import SpectrumSourceField from .. import __version__ @@ -620,28 +617,12 @@ def prepare_source(self, source): # TODO: this is not sufficiently general # TODO: Maybe move this to source_fields?? - for field in source.fields: - if not isinstance(field, SpectrumSourceField): - continue - for ispec, spec in field.spectra.items(): - # Put on fov wavegrid - wave_min = min(fov.meta["wave_min"] for fov in self.fov_manager.fovs) - wave_max = max(fov.meta["wave_max"] for fov in self.fov_manager.fovs) - wave_unit = u.Unit(from_currsys("!SIM.spectral.wave_unit", self.cmds)) - dwave = from_currsys("!SIM.spectral.spectral_bin_width", self.cmds) # Not a quantity - fov_waveset = np.arange(wave_min.value, wave_max.value, dwave) * wave_unit - fov_waveset = fov_waveset.to(u.um) - - field.spectra[ispec] = SourceSpectrum(Empirical1D, - points=fov_waveset, - lookup_table=spec(fov_waveset)) - for cube in source.cube_fields: header, data, wave = cube.header, cube.data, cube.waveset # Need to check whether BUNIT is per arcsec2 or per pixel inunit = u.Unit(header["BUNIT"]) - data = data.astype(np.float32) * inunit + data = data.astype(np.float64, copy=False) * inunit factor = 1 for base, power in zip(inunit.bases, inunit.powers): if (base**power).is_equivalent(u.arcsec**(-2)): @@ -670,35 +651,10 @@ def prepare_source(self, source): cube.header["CUNIT1"] = "deg" cube.header["CUNIT2"] = "deg" - # Put on fov wavegrid - wave_min = min(fov.meta["wave_min"] for fov in self.fov_manager.fovs) - wave_max = max(fov.meta["wave_max"] for fov in self.fov_manager.fovs) - wave_unit = u.Unit(from_currsys("!SIM.spectral.wave_unit", self.cmds)) - dwave = from_currsys("!SIM.spectral.spectral_bin_width", self.cmds) # Not a quantity - fov_waveset = np.arange(wave_min.value, wave_max.value, dwave) * wave_unit - fov_waveset = fov_waveset.to(u.um) - - if (wave.to(u.um).min() > fov_waveset.max() or - wave.to(u.um).max() < fov_waveset.min()): - logger.warning("Source and FOV wavesets disjoint.") - - # Interpolate into new data cube. - # This is done layer by layer for memory reasons. - new_data = np.zeros((fov_waveset.shape[0], data.shape[1], data.shape[2]), - dtype=np.float32) - for j in range(data.shape[1]): - cube_interp = interp1d(wave.to_value(u.um), data[:, j, :], - axis=0, kind="linear", - bounds_error=False, fill_value=0) - new_data[:, j, :] = cube_interp(fov_waveset.value) - - cube.data = new_data - cube.header["CTYPE3"] = "WAVE" - cube.header["CRPIX3"] = 1 - cube.header["CRVAL3"] = wave_min.value - cube.header["CDELT3"] = dwave - cube.header["CUNIT3"] = wave_unit.name - + # Keep the native wavelength samples until the spatial FOV is cut. + # Rebinning the entire source here loses narrow lines and expands + # resolved cubes over the full instrument band unnecessarily. + cube.data = data cube.wcs = WCS(cube.field) return source diff --git a/scopesim/tests/tests_optics/test_FieldOfView.py b/scopesim/tests/tests_optics/test_FieldOfView.py index de961eb0f..61ced1915 100644 --- a/scopesim/tests/tests_optics/test_FieldOfView.py +++ b/scopesim/tests/tests_optics/test_FieldOfView.py @@ -164,123 +164,6 @@ def test_handles_nans(self): class TestMakeCube: - @pytest.mark.xfail(reason="cube flux is broken") - def test_makes_cube_from_table(self): - src_table = so._table_source() # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm - fov = _fov_190_210_um(ndim=3) - fov.extract_from(src_table) - - cube = fov.make_hdu() - - in_sum = 0 - waveset = fov.spectra[0].waveset - for x, y, ref, weight in src_table.fields[0]: - flux = src_table.spectra[ref](waveset).to(u.ph/u.s/u.m**2/u.um).value - in_sum += np.sum(flux) * weight - - out_sum = np.sum(cube.data) - - assert out_sum == approx(in_sum, rel=0.01) - - if PLOTS: - plt.imshow(cube.data[0, :, :], origin="lower") - plt.show() - - @pytest.mark.xfail(reason="cube flux is broken") - def test_makes_cube_from_imagehdu(self): - src_image = so._image_source() # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm - fov = _fov_190_210_um(ndim=3) - fov.extract_from(src_image) - - cube = fov.make_hdu() - - waveset = np.linspace(1.9, 2.1, np.shape(cube)[0]) * u.um - spec = fov.spectra[0](waveset).to(u.ph/u.s/u.m**2/u.um).value - in_sum = np.sum(src_image.fields[0].data) * np.sum(spec) - out_sum = np.sum(cube.data) - - assert out_sum == approx(in_sum, rel=0.01) - - if PLOTS: - plt.imshow(cube.data[0, :, :], origin="lower") - plt.show() - - @pytest.mark.xfail(reason="cube flux is broken") - def test_makes_cube_from_other_cube_imagehdu(self): - import scopesim as sim - sim.rc.__currsys__["!SIM.spectral.spectral_bin_width"] = 0.01 - src_cube = so._cube_source() # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm - fov = _fov_197_202_um(ndim=3) - fov.extract_from(src_cube) - - cube = fov.make_hdu() - - # layer 74 to 77 are extracted by FOV - in_sum = np.sum(src_cube.fields[0].data[74:77, :, :]) - - out_sum = np.sum(cube.data) - - assert out_sum == approx(in_sum) - - if PLOTS: - plt.imshow(cube.data[0, :, :], origin="lower") - plt.show() - - @pytest.mark.xfail(reason="cube flux is broken") - def test_makes_cube_from_two_similar_cube_imagehdus(self): - src_cube = so._cube_source() + so._cube_source(dx=1) # 2 cubes 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm - fov = _fov_197_202_um(ndim=3) - fov.extract_from(src_cube) - - cube = fov.make_hdu() - - # layer 74 to 77 are extracted by FOV - bin_widths = np.array([0.01, 0.02, 0.01])[:, None, None] * 1e4 # um -> AA - in_sum = 2 * np.sum(src_cube.fields[0].data[74:77, :, :] * bin_widths) - - out_sum = np.sum(cube.data) - - assert out_sum == approx(in_sum, rel=1e-3) - - if PLOTS: - plt.imshow(cube.data[0, :, :], origin="lower") - plt.show() - - @pytest.mark.xfail(reason="cube flux is broken") - def test_makes_cube_from_all_types_of_source_object(self): - src_all = so._table_source() + \ - so._image_source(dx=-4, dy=-4) + \ - so._cube_source(weight=1e-8, dx=4) - - fov = _fov_190_210_um(ndim=3) - fov.extract_from(src_all) - - cube = fov.make_hdu() - - # sum up the expected flux in the output cube - # bin_width * half width edge bin * PHOTLAM -> SI - waveset = fov.spectra[0].waveset - - table_sum = 0 - for x, y, ref, weight in src_all.fields[0]: - flux = src_all.spectra[ref](waveset).value - table_sum += np.sum(flux) * weight - - ref = src_all.fields[1].header["SPEC_REF"] - spec = fov.spectra[ref](waveset).value - image_sum = np.sum(src_all.fields[1].data) * np.sum(spec) - - cube_sum = np.sum(src_all.fields[2].data[70:81, :, :]) * 0.02 * 0.95 - - in_sum = table_sum + image_sum + cube_sum - out_sum = np.sum(cube.data) - - assert out_sum == approx(in_sum, rel=0.02) - - if PLOTS: - im = cube.data[0, :, :] - plt.imshow(im, origin="lower", norm="log", vmin=1e-8) - plt.show() @pytest.mark.parametrize("src", [so._table_source(), so._image_source(), @@ -494,49 +377,50 @@ def test_makes_spectrum_from_all_types_of_source_object(self): plt.show() -@pytest.mark.skip(reason="needs split FOVs after subclass refactor") class TestMakeSpectrumImageCubeAllPlayNicely: - def test_make_cube_and_make_spectrum_return_the_same_fluxes(self): - src_all = so._table_source() + \ - so._image_source(dx=-4, dy=-4) + \ - so._cube_source(weight=1e-8, dx=4) - - fov = _fov_190_210_um(ndim=3) - fov.extract_from(src_all) - - # if photlam, units of ph / s / cm2 / AA, else units of ph / s / voxel - cube = fov.make_hdu() - # cube_waves = get_cube_waveset(cube.header) - cube_spectrum = cube.data.sum(axis=2).sum(axis=1) - - # always units of ph / s / cm-2 / AA-1 - waves = fov.waveset - spectrum = fov.make_hdu()(waves).value - - bin_width = 0.02 # um - photlam_to_si = 1e8 # cm-2 AA-1 --> m-2 um-1 - edge_compensation = 0.91 - - spectrum *= bin_width * photlam_to_si * edge_compensation - - if PLOTS: - plt.plot(waves, spectrum, "k") - # plt.plot(cube_waves, cube_spectrum, "r") - plt.show() - - assert np.sum(cube.data) == approx(np.sum(spectrum), rel=0.001) - - def test_make_cube_and_make_image_return_the_same_fluxes(self): - src_all = so._cube_source(weight=1e-8, dx=4) - # so._image_source(dx=-4, dy=-4) # - # so._table_source()# + \ - - fov = _fov_190_210_um() - fov.extract_from(src_all) - - # if photlam, units of ph / s / cm2 / AA, else units of ph / s / voxel - cube_sum = np.sum(fov.make_hdu().data) - # if photlam, units of ph / s / cm2 / AA, else units of ph / s / pixel - image_sum = np.sum(fov.make_hdu().data) - - assert cube_sum == approx(image_sum, rel=0.05) + @pytest.mark.parametrize("kind", ["table", "image", "cube", "two_cubes", "mixed"]) + @pytest.mark.parametrize("spectral_step", [.01, .005]) + def test_photon_integrals_agree_between_representations(self, kind, spectral_step): + from synphot import SourceSpectrum, Empirical1D + from scopesim.commands import UserCommands + from scopesim.optics.optical_train import OpticalTrain + + # A line fully inside the band isolates representation and grid errors + # from partial end bins. Its integral is the native trapezoidal integral. + wave = np.linspace(1.5, 2.5, 101)*u.um + profile = np.exp(-.5*((wave.to_value(u.um)-2.)/.02)**2) + spectrum = SourceSpectrum(Empirical1D, points=wave, lookup_table=profile*u.ph/u.s/u.m**2/u.um) + table = so._table_source() + table.fields[0].spectra = {ref: spectrum for ref in table.fields[0].spectra} + image = so._image_source() + image.fields[0].spectra = {0: spectrum} + cube = so._cube_source() + cube.fields[0].data = cube.fields[0].data[-1][None, :, :]/4 * profile[:, None, None] + cube.fields[0].header['CRVAL3'] = 2. + cube.fields[0].header['CDELT3'] = .01 + cube.fields[0].header['CUNIT3'] = 'um' + sources = {'table': table, 'image': image, 'cube': cube, 'two_cubes': cube+cube, 'mixed': table+image+cube} + source = sources[kind] + spatial_weight = {'table': table.fields[0].field['weight'].sum(), 'image': image.fields[0].data.sum(), 'cube': cube.fields[0].data[50].sum()} + spatial_weight['two_cubes'] = 2*spatial_weight['cube'] + spatial_weight['mixed'] = spatial_weight['table']+spatial_weight['image']+spatial_weight['cube'] + expected = spatial_weight[kind]*np.trapezoid(profile, wave.to_value(u.um)) + source = OpticalTrain().prepare_source(source) + + commands = UserCommands(properties={'!SIM.spectral.spectral_bin_width': spectral_step}) + rates = [] + for cls in (FieldOfView1D, FieldOfView2D, FieldOfView3D): + fov = cls(ho._fov_header(), [1.9, 2.1]*u.um, area=1*u.m**2, cmds=commands) + fov.extract_from(source) + product = fov.make_hdu() + if cls is FieldOfView1D: + samples = product(fov.waveset).to(u.ph/u.s/u.m**2/u.um) + rate = np.trapezoid(samples, fov.waveset.to(u.um))*fov.area + rates.append(rate.to_value(u.ph/u.s)) + elif cls is FieldOfView2D: + rates.append(product.data.sum()) + else: + rate = product.data.sum()*u.Unit(product.header['BUNIT'])*fov.pixel_area*product.header['CDELT3']*u.um + rates.append(rate.to_value(u.ph/u.s)) + # The line tails outside the band contain less than 10^-6 of its flux. + assert rates == approx([expected]*3, rel=2e-5) diff --git a/scopesim/tests/tests_optics/test_spectral_bin_integration.py b/scopesim/tests/tests_optics/test_spectral_bin_integration.py new file mode 100644 index 000000000..a037bb5b5 --- /dev/null +++ b/scopesim/tests/tests_optics/test_spectral_bin_integration.py @@ -0,0 +1,78 @@ +import numpy as np +import pytest +from astropy import units as u +from synphot import SourceSpectrum, Empirical1D +from synphot.units import PHOTLAM +from scopesim.optics.fov import average_spectrum, spectral_bin_edges, rebin_spectral_density +from scopesim.optics.optical_train import OpticalTrain +from scopesim.source.source import Source +from scopesim.source.source_templates import ab_spectrum + + +@pytest.mark.parametrize('resolving_power', [20000, 200000]) +@pytest.mark.parametrize('phase', [0, .2, .5, .8]) +def test_narrow_line_photons_survive_coarse_grid_phase(resolving_power, phase): + center = 5000 + phase*.1 + sigma = center / resolving_power / np.sqrt(8*np.log(2)) + native = center + np.linspace(-8, 8, 193)*sigma + flux = np.exp(-.5*((native - center)/sigma)**2) / (sigma*np.sqrt(2*np.pi)) + spectrum = SourceSpectrum(Empirical1D, points=native*u.AA, lookup_table=flux*PHOTLAM) + wave = np.arange(4998, 5002, .1)*u.AA + density = average_spectrum(spectrum, wave) + photons = (density * np.diff(spectral_bin_edges(wave.to_value(u.AA)))*u.AA).sum() + assert photons.to_value(u.ph/u.s/u.cm**2) == pytest.approx(1, rel=.001) + + +def test_flat_ab_continuum_bin_integral(): + spectrum = ab_spectrum(mag=20) + wave = np.linspace(5000, 5100, 101)*u.AA + edges = spectral_bin_edges(wave.to_value(u.AA)) + # Constant f_nu gives photon density proportional to 1/lambda. + scale = spectrum(5000*u.AA).to_value(PHOTLAM)*5000 + expected = scale*np.log(edges[1:]/edges[:-1])/np.diff(edges) + np.testing.assert_allclose(average_spectrum(spectrum, wave).to_value(PHOTLAM), expected, rtol=1e-6) + + +def test_mutating_spectrum_redshift_updates_native_sampling(): + points = np.array([4999., 4999.99, 5000., 5000.01, 5001.])*u.AA + spectrum = SourceSpectrum(Empirical1D, points=points, lookup_table=[0., 0., 1., 0., 0.]) + average_spectrum(spectrum, np.arange(4998., 5002., .1)*u.AA) + spectrum.z = .01 + fresh = SourceSpectrum(Empirical1D, points=points, lookup_table=[0., 0., 1., 0., 0.], z=.01) + wave = np.arange(5048., 5052., .1)*u.AA + actual = average_spectrum(spectrum, wave) + np.testing.assert_allclose(actual, average_spectrum(fresh, wave), rtol=1e-12) + assert actual.sum() > 0 + + +@pytest.mark.parametrize('descending', [False, True]) +def test_cube_rebinning_preserves_flux_spatial_weights_and_support(descending): + wave = np.array([1., 2., 3.]) + weights = np.array([[1., 2.], [3., 4.]]) + density = wave[:, None, None]*weights + if descending: + wave, density = wave[::-1], density[::-1] + rebinned = rebin_spectral_density(wave, density, np.array([.5, 1.5, 2.5, 3.5])) + np.testing.assert_allclose(rebinned, np.array([.625, 2., 1.375])[:, None, None]*weights) + np.testing.assert_allclose(rebinned.sum(axis=0), 4*weights) + + +def test_prepare_source_keeps_native_line_sampling(): + wave = np.linspace(4999.9, 5000.1, 101)*u.AA + spectrum = SourceSpectrum(Empirical1D, points=wave, lookup_table=np.exp(-((wave.value - 5000)/.01)**2)*PHOTLAM) + source = Source(x=[0.], y=[0.], ref=[0], weight=[1.], spectra=[spectrum]) + train = OpticalTrain() + prepared = train.prepare_source(source) + np.testing.assert_array_equal(prepared.spectra[0].waveset, spectrum.waveset) + + +def test_prepare_cube_retains_faint_contrast_in_double_precision(): + from scopesim.tests.mocks.py_objects.source_objects import _cube_source + source = _cube_source() + field = source.cube_fields[0] + field.data.fill(1e8) + field.data[:, 0, 0] += 1. + prepared = OpticalTrain().prepare_source(source) + data = prepared.cube_fields[0].data + contrast = (data[:, 0, 0] - data[:, 0, 1])/data[:, 0, 1] + np.testing.assert_allclose(contrast, 1e-8, rtol=1e-6) From 4e2b76dde36347fe4da22848cd636fd0d4edabb6 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:42:47 -0700 Subject: [PATCH 5/7] Model configured avalanche noise and warn at detector exposure limits Source commits: e2cefad8 --- scopesim/effects/electronic/__init__.py | 4 +- scopesim/effects/electronic/exposure.py | 30 +++++++ scopesim/effects/electronic/noise.py | 46 ++++++++-- scopesim/effects/illumination.py | 1 + scopesim/effects/selector_wheel.py | 6 +- .../test_detector_operating_points.py | 85 +++++++++++++++++++ 6 files changed, 163 insertions(+), 9 deletions(-) create mode 100644 scopesim/tests/tests_effects/test_detector_operating_points.py diff --git a/scopesim/effects/electronic/__init__.py b/scopesim/effects/electronic/__init__.py index e086896cf..991d1aa33 100644 --- a/scopesim/effects/electronic/__init__.py +++ b/scopesim/effects/electronic/__init__.py @@ -52,8 +52,8 @@ def apply_to(self, obj, **kwargs): from .electrons import LinearityCurve, ADConversion, InterPixelCapacitance from .noise import (Bias, PoorMansHxRGReadoutNoise, BasicReadoutNoise, - PixelResponseNonUniformity, ShotNoise, DarkCurrent, ReadoutNoiseFromTable) -from .exposure import AutoExposure, ExposureIntegration, ExposureOutput + PixelResponseNonUniformity, ShotNoise, AvalancheShotNoise, DarkCurrent, ReadoutNoiseFromTable) +from .exposure import AutoExposure, ExposureIntegration, ExposureOutput, ExposureLimit from .pixels import ReferencePixelBorder, BinnedImage, UnequalBinnedImage, OnChipBinning from .dmps import DetectorModePropertiesSetter diff --git a/scopesim/effects/electronic/exposure.py b/scopesim/effects/electronic/exposure.py index f11237375..328eefb90 100644 --- a/scopesim/effects/electronic/exposure.py +++ b/scopesim/effects/electronic/exposure.py @@ -2,6 +2,7 @@ """Exposure actions.""" from typing import ClassVar +from warnings import warn import numpy as np @@ -293,3 +294,32 @@ def apply_to(self, obj, **kwargs): obj.data = obj.data * dit * ndit return obj + + + +class ExposureLimit(Effect): + """Warn about expected charge per integration without changing the exposure. + + Apply after integration, dark current and on-chip binning, before noise. + The configured usable well refers to a charge packet at this readout stage. + Software binning is deliberately later and does not change this limit. + """ + + required_keys = {"usable_well", "ndit"} + z_order = (819,) + + def __init__(self, **kwargs): + super().__init__(**kwargs) + check_keys(self.meta, self.required_keys, action="error") + + def apply_to(self, obj, **kwargs): + if isinstance(obj, Detector): + well = from_currsys(self.meta["usable_well"], self.cmds) + ndit = from_currsys(self.meta["ndit"], self.cmds) + peak = obj.data.max() / ndit + obj.header["PEAKDIT"] = (peak, "Expected peak charge per integration [electron]") + obj.meta["WELLDEP"] = well + obj.header["WELLDEP"] = (well, "Usable charge per integration [input electron]") + if peak > well: + warn(f"Detector {obj.det_id}: peak {peak:.1f} e per integration exceeds usable well {well:.1f} e; DIT/NDIT are unchanged.", UserWarning, stacklevel=2) + return obj diff --git a/scopesim/effects/electronic/noise.py b/scopesim/effects/electronic/noise.py index 8c6aef4fa..e9d9198a0 100644 --- a/scopesim/effects/electronic/noise.py +++ b/scopesim/effects/electronic/noise.py @@ -141,10 +141,9 @@ def create_rng(self, det_id: int | None = None) -> np.random.Generator: New Generator instance. """ - if det_id is not None: - seed = [det_id, *self.random_seed] - else: - seed = self.random_seed + seed = self.random_seed + if det_id is not None and seed is not None: + seed = [det_id, *seed] return np.random.default_rng(seed) def plot(self, det): @@ -194,7 +193,7 @@ def _apply_to_det(self, det: Detector) -> None: logger.debug("Apply %s to %s", self.display_name, det) ndit = from_currsys(self.meta["ndit"], self.cmds) - det.data = self(det.data, ndit) + det.data = self(det.data, ndit, det_id=det.det_id) class PoorMansHxRGReadoutNoise(BasicReadoutNoise): @@ -427,6 +426,9 @@ class ShotNoise(RandomEffect): """ z_order: ClassVar[tuple[int, ...]] = (820,) + def _apply_to_det(self, det): + det.data = self(det.data, det_id=det.det_id) + def __call__(self, data: ArrayLike, det_id: int | None = None) -> NDArray: rng = self.create_rng(det_id) @@ -449,3 +451,37 @@ def __call__(self, data: ArrayLike, det_id: int | None = None) -> NDArray: ) return data + + + +class AvalancheShotNoise(ShotNoise): + """Input-referred photon multiplication noise, with independent dark charge. + + Conditional on N detected photons, a gamma variate of mean N and variance + N(F-1) gives total photon variance F*mu. F is a variance multiplier, not + an RMS factor. The adopted input-referred dark term remains Poisson. + Apply after ExposureIntegration and DarkCurrent, before read noise. + """ + + required_keys = {"excess_noise_factor", "dark_current", "dit", "ndit"} + + @property + def excess_noise_factor(self): + return from_currsys(self.meta["excess_noise_factor"], self.cmds) + + def __call__(self, data, dark_charge=0., det_id=None): + factor = self.excess_noise_factor + if factor < 1: + raise ValueError("The avalanche variance multiplier must be at least one.") + if factor == 1: + return super().__call__(data, det_id=det_id) + rng = self.create_rng(det_id) + photons = rng.poisson(data - dark_charge) + signal = rng.gamma(photons / (factor - 1), factor - 1) + return signal + rng.poisson(dark_charge, size=data.shape) + + def _apply_to_det(self, det): + dark = from_currsys(self.meta["dark_current"], self.cmds) + dit = from_currsys(self.meta["dit"], self.cmds) + ndit = from_currsys(self.meta["ndit"], self.cmds) + det.data = self(det.data, dark * dit * ndit, det_id=det.det_id) diff --git a/scopesim/effects/illumination.py b/scopesim/effects/illumination.py index 65b5b8181..5179aa7db 100644 --- a/scopesim/effects/illumination.py +++ b/scopesim/effects/illumination.py @@ -212,6 +212,7 @@ def _effect_from_spec(spec, cmds=None): def _detector_qe_from_spec(spec, filename=None, cmds=None): + spec = from_currsys(spec, cmds) if spec is not None: if isinstance(spec, str): return SpectralQuantumEfficiency(filename=spec, cmds=cmds) diff --git a/scopesim/effects/selector_wheel.py b/scopesim/effects/selector_wheel.py index 7ee03f242..333635bb1 100644 --- a/scopesim/effects/selector_wheel.py +++ b/scopesim/effects/selector_wheel.py @@ -13,7 +13,7 @@ import importlib from numbers import Integral -from ..utils import (check_keys, get_logger, real_colname) +from ..utils import from_currsys, check_keys, get_logger, real_colname from .effects import Effect from ..optics.fov_volume_list import FovVolumeList from ..optics.fov import FieldOfView @@ -62,8 +62,10 @@ def __init__(self, **kwargs): self.wheel_effects = {} for wheel_entry in self.meta["wheel"]: selector_value = wheel_entry["selector_value"] # can be single value or list of values - effect_class_name = wheel_entry["effect_class"] + effect_class_name = from_currsys(wheel_entry["effect_class"], self.cmds) effect_kwargs = wheel_entry.get("effect_kwargs", {}) + if isinstance(effect_kwargs, str): + effect_kwargs = from_currsys(effect_kwargs, self.cmds) # Dynamically import the effect class effect_module = importlib.import_module("scopesim.effects") diff --git a/scopesim/tests/tests_effects/test_detector_operating_points.py b/scopesim/tests/tests_effects/test_detector_operating_points.py new file mode 100644 index 000000000..1763f1d26 --- /dev/null +++ b/scopesim/tests/tests_effects/test_detector_operating_points.py @@ -0,0 +1,85 @@ +import warnings + +import numpy as np +import pytest + +from scopesim import UserCommands +from scopesim.detector import Detector +from scopesim.effects.electronic import AvalancheShotNoise, DarkCurrent, ExposureIntegration, ExposureLimit, OnChipBinning, BasicReadoutNoise +from scopesim.optics.image_plane_utils import header_from_list_of_xy + + +@pytest.mark.parametrize('factor', [1., 1.2]) +@pytest.mark.parametrize('photons,dark', [(0., 4.), (.2, .01), (50., 3.)]) +def test_avalanche_preserves_physical_mean_and_only_amplifies_photon_variance(factor, photons, dark): + cmds = UserCommands(properties={'!SIM.random.seed': 7143}) + effect = AvalancheShotNoise(excess_noise_factor=factor, dark_current=dark, dit=1., ndit=1, cmds=cmds) + data = effect(np.full(400000, photons + dark), dark_charge=dark, det_id=2) + mean, variance = photons + dark, factor * photons + dark + assert data.mean() == pytest.approx(mean, abs=5*np.sqrt(variance/data.size)) + assert data.var() == pytest.approx(variance, rel=.025) + + +def test_independent_integrations_and_read_noise_after_avalanche(): + cmds = UserCommands(properties={'!SIM.random.seed': 736}) + count, photons, dark, ron, factor = 12, 7., .4, .5, 1.2 + effect = AvalancheShotNoise(excess_noise_factor=factor, dark_current=dark, dit=1., ndit=count, cmds=cmds) + data = effect(np.full(400000, count*(photons+dark)), dark_charge=count*dark, det_id=3) + data = BasicReadoutNoise(noise_std=ron, ndit=count, cmds=cmds)(data, ndit=count, det_id=3) + assert data.mean() == pytest.approx(count*(photons+dark), abs=.1) + assert data.var() == pytest.approx(count*(factor*photons + dark + ron**2), rel=.02) + cmds['!OBS.roid'] = 2 + other = effect(np.full(400000, count*(photons+dark)), dark_charge=count*dark, det_id=3) + assert not np.array_equal(data, other) + + +def test_exposure_limit_warns_after_charge_binning_without_changing_counts_or_commands(): + cmds = UserCommands(properties={'!OBS.dit': 10., '!OBS.ndit': 4}) + header = header_from_list_of_xy([-2, 2], [-2, 2], 1., 'D') + header['ID'] = 0 + detector = Detector(header, cmds=cmds) + detector.data[:] = 400. + limit = ExposureLimit(usable_well=15000., ndit='!OBS.ndit', cmds=cmds) + ExposureIntegration(dit='!OBS.dit', ndit='!OBS.ndit', cmds=cmds).apply_to(detector) + DarkCurrent(value=0., dit='!OBS.dit', ndit='!OBS.ndit', cmds=cmds).apply_to(detector) + with warnings.catch_warnings(record=True) as caught: + limit.apply_to(detector) + assert not caught # Summed 16000 e is only 4000 e per integration. + OnChipBinning(binx=2, biny=2).apply_to(detector) + before = detector.data.copy() + with pytest.warns(UserWarning, match='peak 16000.0 e.*DIT/NDIT are unchanged'): + limit.apply_to(detector) + np.testing.assert_array_equal(detector.data, before) + assert detector.header['PEAKDIT'] == 16000 + assert cmds['!OBS.dit'] == 10 and cmds['!OBS.ndit'] == 4 + + + +def test_avalanche_supports_unseeded_readouts(): + cmds = UserCommands(properties={'!SIM.random.seed': None}) + effect = AvalancheShotNoise(excess_noise_factor=1.2, dark_current=0., dit=1., ndit=1, cmds=cmds) + first = effect(np.full(1000, 10.), det_id=3) + second = effect(np.full(1000, 10.), det_id=3) + assert not np.array_equal(first, second) + + +def test_detector_ids_produce_independent_read_noise(): + cmds = UserCommands(properties={'!SIM.random.seed': 7143}) + effect = BasicReadoutNoise(noise_std=2., ndit=1, cmds=cmds) + header = header_from_list_of_xy([-50, 50], [-50, 50], 1., 'D') + header['ID'] = 1 + first = Detector(header.copy(), cmds=cmds) + header['ID'] = 2 + second = Detector(header.copy(), cmds=cmds) + effect.apply_to(first) + effect.apply_to(second) + assert abs(np.corrcoef(first.data.ravel(), second.data.ravel())[0, 1]) < .05 + + +def test_selector_detector_properties_keep_exposure_references_live(): + from scopesim.effects import SelectorWheel + from scopesim.utils import from_currsys + cmds = UserCommands(properties={'!DET.noise_class': 'BasicReadoutNoise', '!DET.noise_kwargs': {'noise_std': 2., 'ndit': '!OBS.ndit'}, '!OBS.ndit': 1}) + effect = SelectorWheel(selector_key='detector_id', wheel=[{'selector_value': 0, 'effect_class': '!DET.noise_class', 'effect_kwargs': '!DET.noise_kwargs'}], cmds=cmds) + cmds['!OBS.ndit'] = 7 + assert from_currsys(effect.get_effect(0).meta['ndit'], cmds) == 7 From 8335e57e26d85a4989722a433ad814dc9f986c3b Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:42:47 -0700 Subject: [PATCH 6/7] Integrate slit illumination and detector pixels with pre-slit ADC shifts Source commits: b6f31f6a, 0701bb57, dcb05181, f9dccd8f, fa2598b7, ee4e4e67, e9f4b77a, 0b15a616, 089852a6, 8a43ec43, a0495f72, 8727eafe, aef3a53d --- scopesim/effects/atmo_dispersion.py | 141 ++++++---- scopesim/effects/psfs/analytical.py | 144 ++++++++-- scopesim/effects/psfs/discrete.py | 7 +- scopesim/effects/spectral_efficiency.py | 7 +- scopesim/effects/spectral_trace_list_utils.py | 145 +++++++--- scopesim/effects/ter_curves.py | 10 +- scopesim/optics/fov.py | 251 +++++++++--------- scopesim/optics/fov_manager.py | 29 +- scopesim/optics/image_plane.py | 2 +- scopesim/optics/image_plane_utils.py | 33 +++ scopesim/optics/optical_train.py | 83 +++++- .../test_AtmosphericDispersion.py | 15 ++ scopesim/tests/tests_effects/test_ConstPSF.py | 34 +++ .../tests/tests_effects/test_MoffatPSF.py | 170 +++++++++++- .../test_calibrated_spectral_trace.py | 60 +++++ .../test_finite_spectral_cells.py | 36 +++ .../tests/tests_effects/test_shift_planes.py | 14 + .../tests_effects/test_spatial_cube_shifts.py | 30 +++ .../tests/tests_optics/test_FieldOfView.py | 6 +- .../tests/tests_optics/test_ImagePlane.py | 9 + .../tests/tests_optics/test_OpticsManager.py | 16 +- .../test_adc_slit_illumination.py | 190 +++++++++++++ .../test_slit_spatial_integration.py | 146 ++++++++++ 23 files changed, 1281 insertions(+), 297 deletions(-) create mode 100644 scopesim/tests/tests_effects/test_finite_spectral_cells.py create mode 100644 scopesim/tests/tests_effects/test_shift_planes.py create mode 100644 scopesim/tests/tests_effects/test_spatial_cube_shifts.py create mode 100644 scopesim/tests/tests_optics/test_adc_slit_illumination.py create mode 100644 scopesim/tests/tests_optics/test_slit_spatial_integration.py diff --git a/scopesim/effects/atmo_dispersion.py b/scopesim/effects/atmo_dispersion.py index 268cd7fee..3fd85fd18 100644 --- a/scopesim/effects/atmo_dispersion.py +++ b/scopesim/effects/atmo_dispersion.py @@ -1,5 +1,5 @@ from typing import ClassVar, Any -from functools import partial +from functools import partial, lru_cache import numpy as np from astropy import units as u from astropy.coordinates import SkyCoord, EarthLocation @@ -9,6 +9,7 @@ from astropy.coordinates import SkyCoord, HADec from scipy.interpolate import make_interp_spline, RegularGridInterpolator from scipy.ndimage import map_coordinates +from scipy.sparse import csr_array from .effects import Effect from ..optics.fov import FieldOfView3D @@ -61,15 +62,17 @@ class ShiftFoV3D(Effect): def __init__(self, **kwargs): super().__init__(**kwargs) + self.meta.setdefault("method", None) - if self.meta.get("method", None) not in [None, "roll", "map_coords"]: + if self.meta["method"] not in [None, "roll", "map_coords"]: raise ValueError("Method must be 'roll' or 'map_coords' if provided.") def get_shifts(self, obj): """ Shifts calculated from obj (FoV3D) properties or from init data. Overwrite this function when sub-classing. - Should return interp spline values for dx and dy at data wavelength bins in arcsec + Return dx and dy at cube wavelength bins in arcsec. The FOV can be + unrendered when its incoming illumination needs spatial padding. """ if self.data is not None: if isinstance(self.data, Table): @@ -81,14 +84,28 @@ def get_shifts(self, obj): logger.warning("No shifts as a function of wavelength provided in init.") raise NotImplementedError - 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 - return dx(wave), dy(wave) + wave = _cube_wavelengths(obj) + return dx(wave), dy(wave) + + def _padded_header(self, header, shifts, rounded=False): + """Spatial support needed before applying this displacement.""" + header = header.copy() + # Cubic interpolation needs neighboring samples as well as the shift. + guard = 0 if self.meta['method'] == 'roll' else 4 + for axis, displacement in enumerate(shifts, start=1): + cells = (displacement/(abs(header[f'CDELT{axis}'])*u.deg)).decompose().value + if rounded: + cells = np.round(cells, 1) + if self.meta['method'] == 'roll': + cells = np.round(np.asarray(cells, dtype=np.float32)) + padding = int(np.ceil(np.max(np.abs(cells)))) + guard if np.any(cells) else 0 + header[f'NAXIS{axis}'] += 2*padding + header[f'CRPIX{axis}'] += padding + return header def apply_to(self, obj, **kwargs): if isinstance(obj, FieldOfView3D) and obj.hdu is not None and obj.hdu.header['NAXIS'] == 3: - dxs, dys = self.get_shifts(obj) + dxs, dys = kwargs['_shifts'] if '_shifts' in kwargs else self.get_shifts(obj) # get bin widths in spatial axis and divide to get dxs and dys in steps dxs = np.round((dxs / (obj.hdu.header["CDELT1"]*u.deg)).decompose().value, 1) @@ -119,19 +136,29 @@ def roll_shift(data, dx, dy): raise ValueError("dx and dy must have shape (nwave,)") out = np.empty_like(data) - for w in range(nw): - out[w] = np.roll(data[w], (dy[w], dx[w]), axis=(0, 1)) - # zero out the rolled in part, select the correct slice based on sign of dxint, dyint - xslice = slice(0, dx[w]) if dx[w] > 0 else slice(nx-dx[w], nx) - yslice = slice(0, dy[w]) if dy[w] > 0 else slice(ny-dy[w], ny) - out[w][yslice, xslice] = 0 + boundaries = np.r_[0, np.flatnonzero((np.diff(dx) != 0) | (np.diff(dy) != 0))+1, nw] + for start, stop in zip(boundaries[:-1], boundaries[1:]): + xshift, yshift = dx[start], dy[start] + if abs(xshift) >= nx or abs(yshift) >= ny: + out[start:stop] = 0 + continue + if xshift == 0 and yshift == 0: + out[start:stop] = data[start:stop] + continue + out[start:stop] = np.roll(data[start:stop], (yshift, xshift), axis=(1, 2)) + if xshift != 0: + xslice = slice(0, xshift) if xshift > 0 else slice(nx + xshift, nx) + out[start:stop, :, xslice] = 0 + if yshift != 0: + yslice = slice(0, yshift) if yshift > 0 else slice(ny + yshift, ny) + out[start:stop, yslice, :] = 0 return out @staticmethod - def map_coords_shift(data, dx, dy, chunk_w: int=8, order: int=3): + def map_coords_shift(data, dx, dy, chunk_w: int=64, order: int=3): """ Shift each wavelength plane data[w, :, :] by (dx[w], dy[w]) using - scipy.ndimage.map_coordinates with chunked processing. + scipy.ndimage.map_coordinates independently on each spatial plane. Parameters ---------- @@ -142,7 +169,7 @@ def map_coords_shift(data, dx, dy, chunk_w: int=8, order: int=3): dy : np.ndarray, shape (nwave,) Shift along y (axis=0) in pixels. Positive moves content toward +y. chunk_w : int - Number of wavelength planes to process per chunk. Default 8. + Number of independent wavelength planes processed together. order : int Spline interpolation order. Default 3. @@ -156,36 +183,34 @@ def map_coords_shift(data, dx, dy, chunk_w: int=8, order: int=3): if dx.shape != (nw,) or dy.shape != (nw,): raise ValueError("dx and dy must have shape (nwave,)") - # Base grids for a single chunk (will be broadcast) - yi = np.arange(ny, dtype=np.float32)[None, :, None] # (1, ny, 1) - xi = np.arange(nx, dtype=np.float32)[None, None, :] # (1, 1, nx) - - out = np.empty_like(data) - prefilter = (order > 1) - - for w0 in range(0, nw, chunk_w): - w1 = min(w0 + chunk_w, nw) - k = w1 - w0 - - # wave index grid for this chunk: (k, 1, 1) - wi = np.arange(w0, w1, dtype=np.float32)[:, None, None] - - # To shift content by +dx/+dy, sample input at (x - dx, y - dy) - xcoords = xi - dx[w0:w1][:, None, None] # (k, 1, nx) -> broadcast to (k, ny, nx) - ycoords = yi - dy[w0:w1][:, None, None] # (k, ny, 1) -> broadcast to (k, ny, nx) - - # Broadcast to full chunk volume - xcoords = np.broadcast_to(xcoords, (k, ny, nx)) - ycoords = np.broadcast_to(ycoords, (k, ny, nx)) - wcoords = np.broadcast_to(wi, (k, ny, nx)) - - # Axis order in coords must match axis order of data: (wave, y, x) - coords = np.array([wcoords, ycoords, xcoords], dtype=np.float32) # (3, k, ny, nx) - - out[w0:w1, :, :] = map_coordinates( - data, coords, order=order, mode='nearest', prefilter=prefilter - ) - return out + # The same rounded displacement recurs in many wavelength planes. + # Reuse its spatial interpolation operator, without ever interpolating + # along wavelength or filtering the same profile thousands of times. + shifted = data + for axis, displacement, length in ((2, dx, nx), (1, dy, ny)): + if not np.any(displacement): + continue + out = np.empty(data.shape, dtype=float) + for offset in np.unique(displacement): + operator = _spatial_shift_operator(length, float(offset), order) + planes = np.flatnonzero(displacement == offset) + for start in range(0, len(planes), chunk_w): + selected = planes[start:start + chunk_w] + block = np.moveaxis(shifted[selected], axis, 0) + result = operator @ block.reshape(length, -1) + out[selected] = np.moveaxis(result.reshape(block.shape), 0, axis) + shifted = out + return shifted.astype(data.dtype, copy=False) + + +@lru_cache(maxsize=128) +def _spatial_shift_operator(length, shift, order): + # Apply scipy's interpolation to unit basis profiles once. Cubic spline + # coefficients decay exponentially; omit only weights below 1e-13. + row, column = np.indices((length, length), dtype=float) + response = map_coordinates(np.eye(length), [row, column - shift], order=order, mode='nearest', prefilter=order > 1) + response[abs(response) < 1e-13] = 0 + return csr_array(response.T) class ADShift(ShiftFoV3D): @@ -236,7 +261,7 @@ def apply_to(self, obj, **kwargs): def get_shifts(self, obj: FieldOfView3D): shift_arcsec = self._get_shifts_arcsec(obj) - pos_angle_y = field_rotation_pa_y(obj.hdu.header) + pos_angle_y = field_rotation_pa_y(obj.header if obj.hdu is None else obj.hdu.header) par_angle = get_parallactic_angle(self.target, self.location, self.time) theta = np.deg2rad(par_angle - pos_angle_y) # angle of zenith direction from +y-axis counter-clockwise logger.info(f'Direction of zenith, measured counter-clockwise from +y axis is {int(par_angle - pos_angle_y)}.') @@ -249,9 +274,7 @@ def _get_shifts_arcsec(self, obj: FieldOfView3D): wave_ref = from_currsys(self.meta["wave_ref_um"], self.cmds) * u.um n_ref = self.n(wave_ref) - 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 + wave = _cube_wavelengths(obj) n_wave = self.n(wave) delta_n = n_wave - n_ref return 206265 * delta_n * np.tan(self.zenith_angle.to(u.rad).value) * u.arcsec @@ -309,9 +332,7 @@ def get_shifts(self, obj: FieldOfView3D): """ Computes an AD correction residual for a given zenith angle. """ - 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 + wave = _cube_wavelengths(obj) shifts = np.zeros_like(wave.value) * u.arcsec if self.data is not None and isinstance(self.data, Table): @@ -354,7 +375,7 @@ def get_shifts(self, obj: FieldOfView3D): shifts += ad_shift ad.zenith_angle = self.zenith_angle - pos_angle_y = field_rotation_pa_y(obj.hdu.header) + pos_angle_y = field_rotation_pa_y(obj.header if obj.hdu is None else obj.hdu.header) par_angle = get_parallactic_angle(self.target, self.location, self.time) theta = np.deg2rad(par_angle - pos_angle_y) # angle of zenith direction from +y-axis counter-clockwise logger.info(f'Direction of zenith, measured counter-clockwise from +y axis is {int(par_angle - pos_angle_y)}.') @@ -401,12 +422,20 @@ def _get_ad_shift(self, ad_kwargs): ########################### AD utils ############################### +def _cube_wavelengths(obj): + if obj.hdu is None: + return obj.cube_waveset.to(u.um) + swcs = WCS(obj.hdu.header).spectral + with u.set_enabled_equivalencies(u.spectral()): + return swcs.pixel_to_world(np.arange(swcs.pixel_shape[0])).to(u.um) + + def field_rotation_pa_y(header) -> float: """ Returns PA of detector +Y axis on sky, in deg (N->E convention). 0 deg = Y points north, 90 deg = Y points east """ - w = WCS(header).dropaxis(-1) + w = WCS(header, naxis=2) # local anchor (image center) x0, y0 = w.wcs.crpix lon0, lat0 = w.pixel_to_world(x0, y0) diff --git a/scopesim/effects/psfs/analytical.py b/scopesim/effects/psfs/analytical.py index dc39f5207..1e8bb0ff4 100644 --- a/scopesim/effects/psfs/analytical.py +++ b/scopesim/effects/psfs/analytical.py @@ -2,23 +2,31 @@ """Contains simple Vibration, NCPA, Seeing and Diffraction PSFs.""" from typing import ClassVar, Any -from functools import partial +from functools import partial, lru_cache +from pathlib import Path import numpy as np from astropy import units as u from astropy.convolution import Gaussian2DKernel -from astropy.modeling.models import Moffat2D from scipy.interpolate import make_interp_spline +from scipy.special import gammaln, stdtr +from astropy.wcs import WCS from .. import DataContainer from ...optics import ImagePlane -from ...optics.fov import FieldOfView +from ...optics.fov import FieldOfView, FieldOfView3D from ...utils import (from_currsys, quantify, quantity_from_table, figure_factory, check_keys, get_logger, find_file, get_zenith_angle, get_observation_info_from_cmds) from . import PSF, PoorMansFOV logger = get_logger(__name__) + +def _psf_table_key(filename, cmds): + path = Path(find_file(from_currsys(filename, cmds))) + return str(path), path.read_bytes() + + class AnalyticalPSF(PSF): """Base class for analytical PSFs.""" @@ -261,9 +269,13 @@ def fwhm(self, wavelengths): fwhm = self.meta["fwhm"] if isinstance(fwhm, dict) and {"pivot_wave", "pivot_wave_unit"}.issubset(fwhm): target, location, time, _ = get_observation_info_from_cmds(self.cmds) - key = tuple(fwhm.items()), target, location, time + # The seeing law depends on zenith angle, not the coordinate frame + # instance (SkyCoord equality raises for different observation times). + key = tuple(fwhm.items()), get_zenith_angle(target, location, time) elif isinstance(fwhm, dict): key = tuple(fwhm.items()) + elif isinstance(fwhm, str): + key = _psf_table_key(fwhm, self.cmds) else: key = fwhm @@ -311,8 +323,9 @@ def get_fwhm_interp(self): raise TypeError("fwhm kwarg must be of type dict or float or str") def get_kernel(self, fov): - pixel_scale = fov.header["CDELT1"] * u.deg.to(u.arcsec) + pixel_scale = abs(fov.header["CDELT1"]) * u.deg.to(u.arcsec) pixel_scale = quantify(pixel_scale, u.arcsec) + ystep = abs(fov.header["CDELT2"] / fov.header["CDELT1"]) # Sample the wavelength-dependent seeing law with a bounded number of # representative points. Kernel sizing is driven by the broadest PSF @@ -324,14 +337,14 @@ def get_kernel(self, fov): target = self._target_enclosed_energy() max_ksize = self._max_kernel_size() - ksize = self._minimum_kernel_size(np.max(fwhms)) + ksize = self._minimum_kernel_size(np.max(fwhms) / min(1., ystep)) if max_ksize is not None: ksize = min(ksize, max_ksize) - kernel, norm = self._make_moffat_kernel(gammas, alpha, ksize) + kernel, norm = self._make_moffat_kernel(gammas, alpha, ksize, ystep) while norm < target and (max_ksize is None or ksize < max_ksize): ksize = self._next_kernel_size(ksize, max_ksize) - kernel, norm = self._make_moffat_kernel(gammas, alpha, ksize) + kernel, norm = self._make_moffat_kernel(gammas, alpha, ksize, ystep) if norm < target: logger.warning( @@ -354,6 +367,65 @@ def get_kernel(self, fov): kernel /= norm return kernel + def point_source_image(self, fov, x, y): + """Pixel-integrated effective PSF at an actual sky position [arcsec].""" + dx, dy = abs(fov.header['CDELT1'])*3600, abs(fov.header['CDELT2'])*3600 + xp, yp = WCS(fov.header, naxis=2).all_world2pix(x/3600, y/3600, 0) + xx = (np.arange(fov.header['NAXIS1']) - xp)[None, :]*dx + yy = (np.arange(fov.header['NAXIS2']) - yp)[:, None]*dy + wave = self._sample_kernel_wavelengths(fov.waveset) + gammas = self.fwhm2gamma(self.fwhm(wave), self.alpha).to_value(u.arcsec) + return moffat_cell_integrals(xx, yy, dx, dy, gammas, self.alpha) + + def _apply_convolution(self, obj, spatial_shifts=()): + if isinstance(obj, FieldOfView3D) and obj.hdu is None: + if obj.trace_id is not None: + # Finite cross-slit cells carry a uniform illumination density. + # Resolve the effective PSF before shearing those cells in + # wavelength. Refine only as needed to resolve this effective PSF. + wave = self._sample_kernel_wavelengths(obj.waveset) + width = (obj.meta['eta_max'] - obj.meta['eta_min']).to_value(u.arcsec) + fwhm = self.fwhm(wave).to_value(u.arcsec).min() + obj.meta['effective_psf_fwhm'] = fwhm*u.arcsec + step = fwhm / 10 + rows = max(obj.header['NAXIS2'], int(np.ceil(width / step))) + obj.header['NAXIS2'] = rows + obj.header['CDELT2'] = np.copysign(width / rows / 3600, obj.header['CDELT2']) + obj.header['CRPIX2'] = (rows + 1) / 2 + obj.header['CRVAL2'] = ((obj.meta['eta_min'] + obj.meta['eta_max']) / 2).to_value(u.deg) + length = (obj.meta['xi_max'] - obj.meta['xi_min']).to_value(u.arcsec) + columns = max(obj.header['NAXIS1'], int(np.ceil(6*length/fwhm))) + obj.header['NAXIS1'] = columns + obj.header['CDELT1'] = np.copysign(length/columns/3600, obj.header['CDELT1']) + obj.header['CRPIX1'] = (columns+1)/2 + obj.header['CRVAL1'] = ((obj.meta['xi_min'] + obj.meta['xi_max'])/2).to_value(u.deg) + obj.meta['pixel_area'] = None + obj._volume = None + if 'slit_spectral_width' in obj.meta: + obj.meta['slit_spectral_width'] *= min(1., fwhm / width) + obj.hdu = obj.make_hdu(spatial_psf=self, spatial_shifts=spatial_shifts) + else: + super()._apply_convolution(obj) + + def padded_header(self, fov, point_sources=False): + """Spatial support needed to illuminate this FOV through the PSF.""" + if point_sources: + # Points use analytic cell integrals, so constructing a convolution + # kernel merely to determine source-selection bounds is redundant. + wave = self._sample_kernel_wavelengths(fov.waveset) + gamma = self.fwhm2gamma(self.fwhm(wave), self.alpha).to_value(u.arcsec).max() + radius = gamma*np.sqrt((1-self._target_enclosed_energy())**(1/(1-self.alpha))-1) + padding_x = int(np.ceil(radius/(abs(fov.header['CDELT1'])*3600))) + padding_y = int(np.ceil(radius/(abs(fov.header['CDELT2'])*3600))) + else: + kernel = self.get_kernel(fov) + padding_y, padding_x = kernel.shape[0]//2, kernel.shape[1]//2 + header = fov.header.copy() + for axis, padding in ((1, padding_x), (2, padding_y)): + header[f'NAXIS{axis}'] += 2*padding + header[f'CRPIX{axis}'] += padding + return header + def _sample_kernel_wavelengths(self, waveset: u.Quantity) -> u.Quantity: max_samples = max(2, int(from_currsys( self.meta.get("max_kernel_wavelength_samples", 16), self.cmds))) @@ -387,22 +459,9 @@ def _minimum_kernel_size(self, max_fwhm_pix: float) -> int: kx = float(from_currsys(self.meta.get("kernel_size", 4.0), self.cmds)) return self._ensure_odd_int(kx * max_fwhm_pix) - def _make_moffat_kernel(self, gammas: np.ndarray, alpha:float, ksize: int) -> tuple[np.ndarray, float]: - amplitude = (alpha - 1) / (np.pi * gammas**2) - x, y = np.meshgrid( - np.arange(ksize) - ksize // 2, - np.arange(ksize) - ksize // 2, - ) - cube = Moffat2D.evaluate( - x=x[None, ...], - y=y[None, ...], - amplitude=amplitude[:, None, None], - x_0=0, - y_0=0, - gamma=gammas[:, None, None], - alpha=alpha, - ) - kernel = np.mean(cube, axis=0) + def _make_moffat_kernel(self, gammas, alpha, ksize, ystep=1.): + positions = np.arange(ksize) - ksize//2 + kernel = moffat_cell_integrals(positions[None, :], positions[:, None]*ystep, 1., ystep, gammas, alpha) if from_currsys(self.meta.get("rounded_edges", False), self.cmds): kernel = self._round_kernel_edges(kernel) return kernel, float(np.sum(kernel)) @@ -451,6 +510,33 @@ def fwhm_from_table(table): return lambda wavelengths: interp(wavelengths.to_value(u.um)) +@lru_cache(maxsize=16) +def _legendre_rule(order): + nodes, weights = np.polynomial.legendre.leggauss(order) + nodes.flags.writeable = weights.flags.writeable = False + return nodes, weights + + +def moffat_cell_integrals(x, y, dx, dy, gammas, alpha): + """Integrate normalized Moffat profiles over rectangles, averaging wavelengths. + + Integrate the conditional Student-t profile analytically along y and use + Gauss-Legendre quadrature along x. This resolves narrow profiles without + oversampling a spectral cube. All coordinates share one angular unit. + """ + nodes, weights = _legendre_rule(max(8, int(np.ceil(8*dx/np.min(gammas))))) + xx = x[..., None] + dx*nodes/2 + lower, upper = y[..., None] - dy/2, y[..., None] + dy/2 + integral = np.zeros(np.broadcast_shapes(x.shape, y.shape)) + normalization = np.exp(gammaln(alpha - .5) - gammaln(alpha - 1)) / np.sqrt(np.pi) + for gamma in gammas: + marginal = normalization/gamma * (1 + (xx/gamma)**2)**(.5 - alpha) + scale = np.sqrt((2*alpha - 1)/(gamma**2 + xx**2)) + probability = stdtr(2*alpha - 1, upper*scale) - stdtr(2*alpha - 1, lower*scale) + integral += np.sum(marginal*probability*weights, axis=-1)*dx/2 + return integral/len(gammas) + + class AOEnhanceablePSF(MoffatPSF): """ Wavelength dependent Moffat PSF but with added AO scaling. @@ -505,9 +591,10 @@ def natural_alpha(self): @property def ao_alpha(self): - if self._ao_alpha is None or self._ao_alpha[1] != self.meta['ao_table']: + key = _psf_table_key(self.meta['ao_table'], self.cmds), self.meta['alpha'] + if self._ao_alpha is None or self._ao_alpha[1] != key: _, alpha = self.ao_table_data() - self._ao_alpha = alpha, self.meta['ao_table'] + self._ao_alpha = alpha, key return self._ao_alpha[0] @property @@ -530,9 +617,10 @@ def ao_table_data(self): @property def ao_scale(self): - if self._ao_interp is None or self._ao_interp[1] != self.meta['ao_table']: + key = _psf_table_key(self.meta['ao_table'], self.cmds) + if self._ao_interp is None or self._ao_interp[1] != key: ao_table, _ = self.ao_table_data() - self._ao_interp = self.fwhm_from_table(ao_table), self.meta['ao_table'] + self._ao_interp = self.fwhm_from_table(ao_table), key return self._ao_interp[0] def ao_fwhm(self, wavelengths): diff --git a/scopesim/effects/psfs/discrete.py b/scopesim/effects/psfs/discrete.py index 64d852efc..8b5c33355 100644 --- a/scopesim/effects/psfs/discrete.py +++ b/scopesim/effects/psfs/discrete.py @@ -197,8 +197,6 @@ def get_kernel(self, fov): """Find nearest wavelength and build PSF kernel from file.""" idx = _nearest_index(fov.wavelength, self._waveset) ext = self.kernel_indices[idx] - if ext == self.current_layer_id: - return self.kernel try: spline_order = from_currsys(self.meta["interp_order"], cmds=self.cmds) @@ -207,13 +205,13 @@ def get_kernel(self, fov): logger.warning("PSF parameter interp_order not found, using 1 (linear interpolation).") spline_order = 1 - if fov.hdu.header["NAXIS"] == 3: self.current_layer_id = ext self.make_psf_cube(fov, spline_order) return self.kernel - self.kernel = self._file[ext].data + # Keep the loaded reference independent of callers modifying a kernel. + self.kernel = self._file[ext].data.copy() self.current_layer_id = ext hdr = self._file[ext].header @@ -265,6 +263,7 @@ def make_psf_cube(self, fov, spline_order=1): - fov.hdu.header["CRPIX3"]) + fov.hdu.header["CRVAL3"] ) + lam *= u.Unit(fov.hdu.header['CUNIT3']).to(u.um) if not self.meta.get("interp_psf", True): lam = np.array([lam[len(lam)//2]]) diff --git a/scopesim/effects/spectral_efficiency.py b/scopesim/effects/spectral_efficiency.py index 2baae5bd7..b551b72fc 100644 --- a/scopesim/effects/spectral_efficiency.py +++ b/scopesim/effects/spectral_efficiency.py @@ -112,12 +112,7 @@ def apply_to(self, obj, **kwargs): logger.warning("No grating efficiency for trace %s", trace_id) return obj - 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 - thru = effic.throughput(wave) - obj.hdu.data = obj.hdu.data * thru.value[:, None, None] - return obj + return effic.apply_to(obj) def plot(self): """Plot the grating efficiencies.""" diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py index 3e7fbc71b..09e0af973 100644 --- a/scopesim/effects/spectral_trace_list_utils.py +++ b/scopesim/effects/spectral_trace_list_utils.py @@ -16,6 +16,7 @@ from scipy.interpolate import interp1d, RectBivariateSpline from scipy.ndimage import gaussian_filter1d, map_coordinates +from scipy.sparse import csr_matrix from astropy.table import Table, vstack from astropy.modeling import fitting @@ -90,6 +91,7 @@ class SpectralTrace: - x, y : [mm] """ + _pixel_quadrature_order = 4 _class_params = { "x_colname": "x", "y_colname": "y", @@ -146,6 +148,7 @@ def __init__(self, trace_tbl, cmds=None, **kwargs): # Declaration of other attributes self._xilamimg = None + self._slit_cache = {} self.dlam_per_pix = None self.fit_inverse = True @@ -268,8 +271,8 @@ def _focal_plane_grid(self, fov): xmin_mm, ymin_mm = fpa_wcsd.all_pix2world(xmin, ymin, 0) xmax_mm, ymax_mm = fpa_wcsd.all_pix2world(xmax, ymax, 0) - x_mm = np.linspace(xmin_mm, xmax_mm, sub_naxis1, endpoint=False, dtype=np.float32) - y_mm = np.linspace(ymin_mm, ymax_mm, sub_naxis2, endpoint=False, dtype=np.float32) + x_mm = np.linspace(xmin_mm, xmax_mm, sub_naxis1, endpoint=False) + y_mm = np.linspace(ymin_mm, ymax_mm, sub_naxis2, endpoint=False) return xmin, ymin, x_mm, y_mm, det_wcs def spectral_padding(self, fov): @@ -286,6 +289,9 @@ def spectral_padding(self, fov): step = abs(fov.header['CDELT2'])*3600 derivative = self.dlam_per_pix(wave)/step padding = np.max(abs(derivative))*width/2 + # Resolve the physical slit profile rather than using one fixed + # wavelength step for all orders. The PSF may request a finer step. + fov.meta['slit_spectral_width'] = np.min(abs(derivative))*width*u.um if 'optical_resolving_power' in self.meta: padding += 5*upper/(self.meta['optical_resolving_power']*np.sqrt(8*np.log(2))) return padding*u.um @@ -320,39 +326,57 @@ def map_spectra_to_focal_plane(self, fov): sub_naxis2 = len(y_mm) xmax = xmin + sub_naxis1 ymax = ymin + sub_naxis2 - image = np.zeros((sub_naxis2, sub_naxis1), dtype=np.float32) + image = np.zeros((sub_naxis2, sub_naxis1), dtype=np.float64) self._set_dispersion(wave_min, wave_max, pixsize=pixsize) slit_dispersion = None if "slit_dispersion_colname" in self.meta: slit_dispersion = self.slit_dispersion - xilam = XiLamImage(fov, self.dlam_per_pix, slit_dispersion=slit_dispersion) + xilam = XiLamImage(fov, self.dlam_per_pix, slit_dispersion=slit_dispersion, _cache=self._slit_cache) if "optical_resolving_power" in self.meta: xilam.convolve_lsf(self.meta["optical_resolving_power"]) self._xilamimg = xilam - # focal-plane coordinate images - ximg_fpa, yimg_fpa = np.meshgrid(x_mm, y_mm) - - # Image mapping (xi, lambda) on the focal plane - xi_fpa = self.xy2xi(ximg_fpa, yimg_fpa).astype(np.float32) - lam_fpa = self.xy2lam(ximg_fpa, yimg_fpa).astype(np.float32) - - # The physical slit ends include the outer half-cells. A truncated - # integer index incorrectly admits pixels below the lower boundary. - inside = (xi_fpa >= xi_min) & (xi_fpa < xi_max) - inside &= (lam_fpa >= wave_min) & (lam_fpa < wave_max) - image = xilam.interp(xi_fpa, lam_fpa, grid=False)*inside - - # Scale to ph / s / pixel - image *= self._trace_flux_scale( - ximg_fpa, yimg_fpa, pixsize, pixscale, det_header) + # Integrate the density over each detector pixel. The inverse optical + # map is locally affine within a pixel. Recompute centers and gradients: + # configured transforms can be edited in place between observations. + xx, yy = np.meshgrid(np.arange(sub_naxis1), np.arange(sub_naxis2)) + ximg_fpa, yimg_fpa = det_wcs.all_pix2world(xx, yy, 0) + lam_fpa = self.xy2lam(ximg_fpa, yimg_fpa) + matrix = det_wcs.pixel_scale_matrix + dxidx, dxidy = self.xy2xi.gradient() + dlamdx, dlamdy = self.xy2lam.gradient() + ax, ay = dxidx(ximg_fpa, yimg_fpa), dxidy(ximg_fpa, yimg_fpa) + bx, by = dlamdx(ximg_fpa, yimg_fpa), dlamdy(ximg_fpa, yimg_fpa) + xi_fpa = self.xy2xi(ximg_fpa, yimg_fpa) + xi_x, xi_y = ax*matrix[0, 0] + ay*matrix[1, 0], ax*matrix[0, 1] + ay*matrix[1, 1] + lam_x, lam_y = bx*matrix[0, 0] + by*matrix[1, 0], bx*matrix[0, 1] + by*matrix[1, 1] + scale = self._trace_flux_scale(ximg_fpa, yimg_fpa, pixsize, pixscale, det_header) + quadrature = np.full(2, self._pixel_quadrature_order) + if 'effective_psf_fwhm' in fov.meta: + middle = (wave_min + wave_max)/2 + xm, ym = self.xilam2x(0., middle), self.xilam2y(0., middle) + dxidx, dxidy = self.xy2xi.gradient() + dlamdx, dlamdy = self.xy2lam.gradient() + angular_dispersion = self.slit_dispersion(middle) if slit_dispersion is not None else -self.dlam_per_pix(middle)/pixscale + angular_map = np.array([[dxidx(xm, ym), dxidy(xm, ym)], [dlamdx(xm, ym)/angular_dispersion, dlamdy(xm, ym)/angular_dispersion]]) + angular_steps = np.linalg.norm(angular_map @ det_wcs.pixel_scale_matrix, axis=0) + psf_width = fov.meta['effective_psf_fwhm'].to_value(u.arcsec) + quadrature = np.maximum(quadrature, np.ceil(8*angular_steps/psf_width).astype(int)) + xnodes, xweights = np.polynomial.legendre.leggauss(quadrature[0]) + ynodes, yweights = np.polynomial.legendre.leggauss(quadrature[1]) + for nx, wx in zip(xnodes/2, xweights/2): + for ny, wy in zip(ynodes/2, yweights/2): + xi = xi_fpa + nx*xi_x + ny*xi_y + lam = lam_fpa + nx*lam_x + ny*lam_y + inside = (xi >= xi_min) & (xi < xi_max) & (lam >= wave_min) & (lam < wave_max) + image[inside] += xilam.interp(xi[inside], lam[inside], grid=False) * scale[inside] * wx * wy detector_qe = fov.meta.get("detector_qe") if detector_qe is not None: xpix_img, ypix_img = np.meshgrid( - np.arange(xmin, xmax, dtype=np.float32), - np.arange(ymin, ymax, dtype=np.float32), + np.arange(xmin, xmax, dtype=np.float64), + np.arange(ymin, ymax, dtype=np.float64), ) image = apply_detector_qe_to_trace_image( image, @@ -837,7 +861,41 @@ def __str__(self): return msg -class XiLamImage(): +def _uniform_sum_cdf(position, width_a, width_b): + """CDF of two centered uniform intervals, including a zero-width interval.""" + wide, narrow = np.maximum(width_a, width_b), np.minimum(width_a, width_b) + probability = np.clip(.5 + position / wide, 0, 1) + finite = narrow > 0 + a, b, x = wide[finite], narrow[finite], position[:, finite] + tail = np.clip((a + b)/2 - abs(x), 0, b)**2 / (2*a*b) + probability[:, finite] = np.where(abs(x) <= (a - b)/2, .5 + x/a, np.where(x < 0, tail, 1 - tail)) + return probability + + +def _slit_cell_operator(wavelength, derivative, eta, deta): + """Conservative overlap of wavelength bins sheared across one finite slit cell. + + Wavelength and its derivative use um and um/arcsec; eta uses arcsec. + Locally affine dispersion maps each rectangular (lambda, eta) cell to a + trapezoidal wavelength profile. Integrating its CDF at output edges keeps + the input photon integral, including lines narrower than a wavelength bin. + """ + step = wavelength[1] - wavelength[0] + center = np.arange(len(wavelength)) + derivative * eta / step + spectral_width = 1 + eta * np.gradient(derivative, wavelength) + slit_width = abs(derivative * deta / step) + assert np.all(spectral_width > 0), 'Slit displacement must preserve wavelength ordering.' + span = spectral_width + slit_width + first = np.floor(center - span/2 + .5).astype(int) + output = first[None, :] + np.arange(int(np.ceil(span.max())) + 1)[:, None] + lower = output - .5 - center + probability = _uniform_sum_cdf(lower + 1, spectral_width, slit_width) - _uniform_sum_cdf(lower, spectral_width, slit_width) + source = np.broadcast_to(np.arange(len(wavelength)), output.shape) + valid = (output >= 0) & (output < len(wavelength)) & (probability > 0) + return csr_matrix((probability[valid], (output[valid], source[valid])), shape=(len(wavelength), len(wavelength))) + + +class XiLamImage: """ Class to compute a rectified 2D spectrum. @@ -856,7 +914,7 @@ class XiLamImage(): legacy coupling between sky-grid sampling and detector dispersion. """ - def __init__(self, fov, dlam_per_pix, slit_dispersion=None): + def __init__(self, fov, dlam_per_pix, slit_dispersion=None, _cache=None): # ..todo: we assume that we always have a cube. We use SpecCADO's # add_cube_layer method cube_wcs = WCS(fov.cube.header, key=" ") @@ -888,7 +946,6 @@ def __init__(self, fov, dlam_per_pix, slit_dispersion=None): cube_lam *= u.Unit(wcs_lam.wcs.cunit[0]).to(u.um) # Initialise the array to hold the xi-lambda image - self.image = np.zeros((n_xi, n_lam), dtype=np.float32) self.lam = cube_lam try: dlam_per_pix_val = dlam_per_pix(np.asarray(self.lam)) @@ -897,25 +954,25 @@ def __init__(self, fov, dlam_per_pix, slit_dispersion=None): logger.warning("Using scalar dlam_per_pix = %.2g", dlam_per_pix_val) - dlam_per_arcsec = dlam_per_pix_val / d_eta if slit_dispersion is None else slit_dispersion(self.lam) - for i, eta in enumerate(cube_eta): - # A positive calibrated derivative moves an input line to a larger - # apparent wavelength: sample the source at lambda_out - shift. - lam0 = self.lam + dlam_per_arcsec * eta if slit_dispersion is None else self.lam - dlam_per_arcsec * eta - - # lam0 is the target wavelength. We need to check that this - # overlaps with the wavelength range covered by the cube - if lam0.min() < cube_lam.max() and lam0.max() > cube_lam.min(): - # Only wavelength changes here; xi is already on its native - # grid. Building a 2D spline for every slit cell is redundant. - index = np.searchsorted(cube_lam, lam0, side='right') - 1 - index = np.clip(index, 0, n_lam - 2) - fraction = (lam0 - cube_lam[index])/d_lam - valid = (lam0 >= cube_lam[0]) & (lam0 <= cube_lam[-1]) - plane = fov.cube.data[:, i, :] - self.image += ((1 - fraction[:, None])*plane[index] + fraction[:, None]*plane[index + 1]).T*valid - - self.image *= abs(d_eta) # ph/s/um/arcsec2 --> ph/s/um/arcsec + derivative = -dlam_per_pix_val / d_eta if slit_dispersion is None else slit_dispersion(self.lam) + derivative = np.broadcast_to(derivative, self.lam.shape) + # The operator depends only on these evaluated physical arrays. Include + # every element, so direct WCS/table/callable changes cannot alias it. + key = (cube_lam.tobytes(), cube_eta.tobytes(), d_eta, derivative.tobytes()) + if _cache is None: + _cache = {} + if _cache.get('key') != key: + _cache.clear() + operators = [_slit_cell_operator(cube_lam, derivative, eta, d_eta).tocoo() for eta in cube_eta] + rows = np.concatenate([operator.row for operator in operators]) + columns = np.concatenate([operator.col*n_eta+i for i, operator in enumerate(operators)]) + weights = np.concatenate([operator.data for operator in operators])*abs(d_eta) + # Interleaved wavelength/slit cells match the contiguous cube layout, + # avoiding a separate strided dense copy for every across-slit row. + _cache['operator'] = csr_matrix((weights, (rows, columns)), shape=(n_lam, n_lam*n_eta)) + _cache['key'] = key + cells = fov.cube.data.reshape(n_lam*n_eta, n_xi) + self.image = (_cache['operator'] @ cells).T # WCS for the xi-lambda image, i.e. the rectified 2D spectrum # Default WCS with xi in arcsec diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py index c0d68eb2e..e35d2129f 100644 --- a/scopesim/effects/ter_curves.py +++ b/scopesim/effects/ter_curves.py @@ -232,7 +232,15 @@ def apply_to(self, obj, **kwargs): ## If obj is a FieldOfView, apply the throughput if isinstance(obj, FieldOfView): thru = self.throughput - if obj.hdu is not None and obj.hdu.header['NAXIS'] == 3: + if obj.hdu is None: + # Retain native spectral knots until spatial rendering, so + # narrow lines survive multiplication by the throughput. + for field in obj.fields: + if isinstance(field, CubeSourceField): + field.data *= thru(field.waveset).value[:, None, None] + elif isinstance(field, SpectrumSourceField): + field.spectra = {ref: spec * thru for ref, spec in field.spectra.items()} + elif obj.hdu.header['NAXIS'] == 3: 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 diff --git a/scopesim/optics/fov.py b/scopesim/optics/fov.py index a2b2ffa9b..3f56c4e28 100644 --- a/scopesim/optics/fov.py +++ b/scopesim/optics/fov.py @@ -1,12 +1,14 @@ # -*- coding: utf-8 -*- """Defines FieldOfView class.""" +from copy import copy from warnings import warn from itertools import chain from collections.abc import Iterable, Generator import numpy as np from scipy.integrate import cumulative_trapezoid +from scipy.signal import fftconvolve from astropy import units as u from astropy.io import fits @@ -73,6 +75,7 @@ def __init__(self, header, waverange, detector_header=None, cmds=None, **kwargs) "wave_max": quantify(waverange[1], u.um), "wave_bin_n": 1, "wave_bin_type": "linear", + 'spectral_padding': 0*u.um, "area": 0 * u.m**2, "pixel_area": None, # [arcsec] @@ -132,7 +135,7 @@ def pixel_area(self): self.meta["pixel_area"] = pixel_area(self.header) return self.meta["pixel_area"] - def extract_from(self, src) -> None: + def extract_from(self, src, spatial_header=None, _spectral_wavesets=None) -> None: """ Extract relevent fields from source object. @@ -154,17 +157,22 @@ def extract_from(self, src) -> None: """ assert isinstance(src, Source), f"expected Source: {type(src)}" + self.fields = [] + self.hdu = None - fields_in_fov = list(self.get_fields_in_fov(src.fields)) + corners_arcsec, _ = self.get_corners("arcsec") + if spatial_header is not None: + corners_arcsec = imp_utils.calc_footprint(spatial_header, new_unit='arcsec') + corners_deg = corners_arcsec / 3600 + fields_in_fov = list(self.get_fields_in_fov(src.fields, corners_arcsec)) if not fields_in_fov: logger.warning("No fields in FOV.") else: logger.debug("%d fields in FOV", len(fields_in_fov)) - corners_arcsec, _ = self.get_corners("arcsec") - corners_deg, _ = self.get_corners("deg") minmax = array_minmax(corners_arcsec) * u.arcsec + source_range = self.waverange + np.array([-1, 1])*self.meta['spectral_padding'] for field in fields_in_fov: if isinstance(field, TableSourceField): @@ -176,7 +184,7 @@ def extract_from(self, src) -> None: new_fld = TableSourceField( field=extracted, spectra={ - ref: extract_range_from_spectrum(spec, self.waverange) + ref: extract_range_from_spectrum(spec, source_range, _waveset=None if _spectral_wavesets is None else _spectral_wavesets[spec]) for ref, spec in field.spectra.items() if ref in extracted["ref"] }, @@ -190,7 +198,7 @@ def extract_from(self, src) -> None: new_fld = ImageSourceField( field=extracted, spectra={ - ref: extract_range_from_spectrum(spec, self.waverange) + ref: extract_range_from_spectrum(spec, source_range, _waveset=None if _spectral_wavesets is None else _spectral_wavesets[spec]) for ref, spec in field.spectra.items() }, ) @@ -208,7 +216,7 @@ def extract_from(self, src) -> None: field=None, header=field.header, spectra={ - ref: extract_range_from_spectrum(spec, self.waverange) + ref: extract_range_from_spectrum(spec, source_range, _waveset=None if _spectral_wavesets is None else _spectral_wavesets[spec]) for ref, spec in field.spectra.items() }, ) @@ -246,12 +254,13 @@ def view(self, sub_pixel: bool | None = None): def flatten(self): """If cube, collapse along first axis.""" if self.hdu and self.hdu.header["NAXIS"] == 3: - image = np.sum(self.hdu.data, axis=0) + image = np.sum(self.hdu.data, axis=0, dtype=np.float64) self.hdu.data = image - def get_fields_in_fov(self, fields: Iterable[SourceField]) -> Generator: + def get_fields_in_fov(self, fields: Iterable[SourceField], fov_corners=None) -> Generator: """Return True if Source.field footprint is inside FOV footprint.""" - fov_corners, _ = self.get_corners("arcsec") + if fov_corners is None: + fov_corners, _ = self.get_corners("arcsec") for field in fields: field_corners = field.get_corners("arcsec") @@ -356,7 +365,6 @@ def extract_area_from_imagehdu(self, imagehdu, corners): # up a duplicated or a dropped row/column. new_wcs = image_wcs.deepcopy() new_wcs.wcs.crpix = image_wcs.wcs.crpix - xy0p - new_wcs.wcs.cdelt = np.abs(new_wcs.wcs.cdelt) new_naxis = xy1p - xy0p logger.debug("orig image wcs: %s", image_wcs) logger.debug("new cutout wcs: %s", new_wcs) @@ -400,10 +408,14 @@ def _extract_volume_from_cube(self, cubehdu, new_hdr, xy0p, xy1p): hdu_waves = hdu_waves.quantity.to(u.um) wdel = hdr["CDELT3"] * u.Unit(hdr.get("CUNIT3", "AA")) - fov_wmin, fov_wmax = self.waverange + fov_wmin, fov_wmax = self.waverange + np.array([-1, 1])*self.meta['spectral_padding'] # need to go [+/-] half a bin - mask = ((hdu_waves > fov_wmin - 0.5 * wdel) * - (hdu_waves <= fov_wmax + 0.5 * wdel)) + margin = abs(wdel)/2 + if isinstance(self, FieldOfView3D): + # Spectroscopic inputs remain densities until bin integration; + # retain neighbours on both sides of the output bin edges. + margin = abs(wdel) + np.diff(self.cube_waveset[:2])[0]/2 + mask = (hdu_waves > fov_wmin - margin) & (hdu_waves <= fov_wmax + margin) # OC [2021-12-14] if fov range is not covered by the source return nothing if not np.any(mask): @@ -418,12 +430,11 @@ def _extract_volume_from_cube(self, cubehdu, new_hdr, xy0p, xy1p): f0 = (abs(hdu_waves[i0p] - fov_wmin + 0.5 * wdel) % wdel) / wdel # blue edge f1 = (abs(hdu_waves[i1p] - fov_wmax - 0.5 * wdel) % wdel) / wdel # red edge - data = cubehdu.data[i0p:i1p+1, - xy0p[1]:xy1p[1], - xy0p[0]:xy1p[0]] - data[0, :, :] *= f0.to(1).round(11).value # decompose dimensionless - if i1p > i0p: - data[-1, :, :] *= f1.to(1).round(11).value + data = cubehdu.data[i0p:i1p+1, xy0p[1]:xy1p[1], xy0p[0]:xy1p[0]].copy() + if not isinstance(self, FieldOfView3D): + data[0] *= f0.to_value(1).round(11) + if i1p > i0p: + data[-1] *= f1.to_value(1).round(11) # w0, w1 : the closest cube wavelengths outside the fov edge wavelengths # fov_waves : the fov edge wavelengths @@ -745,7 +756,7 @@ def _make_cubefields(self): u.Unit(field.header["CUNIT3"]) ).to(u.Angstrom) # First collapse to image, then convert units - image = np.sum(field.data, axis=0) * PHOTLAM/u.arcsec**2 + image = np.sum(field.data, axis=0, dtype=np.float64) * PHOTLAM/u.arcsec**2 image = (image * self.area * field.pixel_area * spectral_bin_width).to(u.ph/u.s) logger.debug("2D FOV make_cubefields: image.mean() = %f %s", @@ -811,8 +822,8 @@ def _make_tablefields(self, fov_waveset, bin_widths, use_photlam=False): else: # Note: these had x/ypix+0.5 until a06ab75 # TODO: could these be something more numpythonic grid-ish? - x = np.array(xpix).astype(int) - y = np.array(ypix).astype(int) # quickest way to round + x = np.floor(xpix + .5).astype(int) + y = np.floor(ypix + .5).astype(int) flux = np.array([fluxes[int(ref)] for ref in field.field["ref"]]) yield flux, np.array(field.field["weight"]), x, y @@ -870,14 +881,11 @@ def make_hdu(self): # enabled, it is therefore not necessary to deploy the fix # below in the else-branch. assert not isinstance(x, Iterable), "x must be an integer" - try: + if 0 <= x < canvas_image_hdu.data.shape[1] and 0 <= y < canvas_image_hdu.data.shape[0]: canvas_image_hdu.data[y, x] += flux * weight - except IndexError: - continue else: # Mask out any stars that were pushed out of the fov by rounding - mask = ((x < canvas_image_hdu.data.shape[1]) * - (y < canvas_image_hdu.data.shape[0])) + mask = (x >= 0) & (y >= 0) & (x < canvas_image_hdu.data.shape[1]) & (y < canvas_image_hdu.data.shape[0]) # This used to contain this line: # canvas_image_hdu.data[y[mask], x[mask]] += flux[mask] * weight[mask] @@ -918,6 +926,16 @@ class FieldOfView3D(FieldOfView): """ + @property + def cube_waveset(self): + wave_unit = u.Unit(from_currsys('!SIM.spectral.wave_unit', self.cmds)) + step = from_currsys('!SIM.spectral.spectral_bin_width', self.cmds) + if 'slit_spectral_width' in self.meta: + step = min(step, self.meta['slit_spectral_width'].to_value(wave_unit)/20) + start, stop = self.meta['wave_min'].value, self.meta['wave_max'].value + padding = int(np.ceil(self.meta['spectral_padding'].to_value(wave_unit)/step)) + return (start + np.arange(-padding, int((stop - start)/step) + padding)*step)*wave_unit + def _make_cubefields(self, fov_waveset): """ Find Cube fields. @@ -936,7 +954,7 @@ def _make_cubefields(self, fov_waveset): field_hdu.header.update(CRPIX3=1., CRVAL3=fov_waveset[0].to_value(u.um), CDELT3=np.diff(fov_waveset[:2])[0].to_value(u.um), CUNIT3='um', CTYPE3='WAVE') yield field_hdu - def _make_imagefields(self, fov_waveset, spline_order): + def _make_imagefields(self, fov_waveset, spline_order, spatial_psf=None): """ Find Image fields. @@ -951,20 +969,19 @@ def _make_imagefields(self, fov_waveset, spline_order): # ImageHDUs have photons/second/pixel. field_hdu = self._make_scaled_hdu(field) - canvas_image_hdu = fits.ImageHDU( - data=zeros_from_header(self.header, ndims=2), - header=self.header) - canvas_image_hdu = imp_utils.add_imagehdu_to_imagehdu( - field_hdu, - canvas_image_hdu, - spline_order=spline_order, - conserve_flux=True, - differential=True, - ) + if spatial_psf is None: + image = imp_utils.project_spatial_density(field_hdu, self.header) + else: + padded = spatial_psf.padded_header(self) + image = imp_utils.project_spatial_density(field_hdu, padded) + image = fftconvolve(image, spatial_psf.get_kernel(self), mode='same') + px = (padded['NAXIS1'] - self.header['NAXIS1'])//2 + py = (padded['NAXIS2'] - self.header['NAXIS2'])//2 + image = image[py:py + self.header['NAXIS2'], px:px + self.header['NAXIS1']] spec = average_spectrum(field.spectrum, fov_waveset) # 2D * 1D -> 3D - field_cube = canvas_image_hdu.data[None, :, :] * spec[:, None, None] + field_cube = image[None, :, :] * spec[:, None, None] yield field_cube.value def _make_tablefields(self, fov_waveset): @@ -1003,7 +1020,7 @@ def _make_tablefields(self, fov_waveset): for i, j, k in zip(xs, ys, fracs): yield flux_vector * k, i, j else: - yield flux_vector, int(xpix), int(ypix) + yield flux_vector, int(np.floor(xpix + .5)), int(np.floor(ypix + .5)) def _make_backfields(self, fov_waveset): for field in self._get_background_fields(): @@ -1015,81 +1032,47 @@ def _make_backfields(self, fov_waveset): spec = average_spectrum(field.spectrum, fov_waveset) yield spec[:, None, None].value - def make_hdu(self): - """ - TBA. - - Used for IFUs, slit spectrographs, and coherent MOSs (e.g.KMOS) - - Returned cube units are ``ph s-1 voxel-1`` - - .. note:: ``self.make_cube()`` does NOT store anything in ``self.cube`` - - self.cube and self.make_cube() are deliberately kept seperately - so that self.cube will not be accidently overwritten by a rogue - call from an Effect object. - - All Effect objects should specifically test whether - ``self.cube is None`` before assigning a new cube it - - The cube is made with these steps: - - 1. Make waveset and canvas cube:: - - if at least one cube: - set waveset to equal largest cube waveset - else: - make waveset from self.meta values - make canvas cube based on waveset of largest cube and NAXIS1,2 from fov.header - - 2. Find Cube fields (see ``FieldOfView._make_cube_cubefields()``). - 3. Find Image fields (see ``FieldOfView._make_cube_imagefields()``). - 4. Find Table fields (see ``FieldOfView._make_cube_tablefields()``). - - ``PHOTLAM = ph / (cm2 * s * AA)``. - Original source fields are in units of: - - - tables: (PHOTLAM in spectrum) - - images: arcsec-2 (PHOTLAM in spectrum) - - cubes: PHOTLAM arcsec-2 - - .. warning:: Input Images and Cubes should have units of PHOTLAM arcsec-2 - - Returns - ------- - canvas_cube_hdu : fits.ImageHDU - [ph s-1 um-1 arcsec-2] # as needed by SpectralTrace + def make_hdu(self, spatial_psf=None, spatial_shifts=()): + """Render spectral surface-brightness density in ph/s/um/arcsec2. + Native spectra are integrated into the cube wavelength cells. Point + sources use cell-integrated PSF illumination; resolved fields are + convolved before the physical aperture is clipped. Collecting area is + applied here. Slit and detector cell integrals are applied by the trace. """ + if spatial_shifts: + # Keep incoming illumination outside the aperture until every + # atmospheric displacement is applied. The physical slit metadata + # and cell sizes stay unchanged; only temporary spatial support grows. + shifts = [(effect, effect.get_shifts(self)) for effect in spatial_shifts] + padded = copy(self) + padded.header = self.header.copy() + for effect, displacement in shifts: + # This is the final PSF grid. Earlier source selection keeps + # conservative support because its coarser grid can still change. + padded.header = effect._padded_header(padded.header, displacement, rounded=True) + padded.hdu = padded.make_hdu(spatial_psf=spatial_psf) + for effect, displacement in shifts: + effect.apply_to(padded, _shifts=displacement) + px = (padded.header['NAXIS1'] - self.header['NAXIS1'])//2 + py = (padded.header['NAXIS2'] - self.header['NAXIS2'])//2 + header = padded.hdu.header.copy() + header.update(CRPIX1=self.header['CRPIX1'], CRPIX2=self.header['CRPIX2']) + data = padded.hdu.data[:, py:py + self.header['NAXIS2'], px:px + self.header['NAXIS1']].copy() + return fits.ImageHDU(data=data, header=header) + # 1. Make waveset and canvas cube (area, bin_width are applied at end) # TODO: Why is this not self.waveset? What's different? # -> For non-cube input but cube output (e.g. 2D image in spec mode), # waverange needs to be resampled (I guess) to spectral_bin_width, # but self.waverange can only access the fields. # -> Perhaps change that?? - wave_unit = u.Unit(from_currsys("!SIM.spectral.wave_unit", self.cmds)) - fov_waveset = np.arange( - self.meta["wave_min"].value, self.meta["wave_max"].value, - from_currsys("!SIM.spectral.spectral_bin_width", self.cmds)) * wave_unit - - # Note: There is an edge case where some floating point madness - # results in the arange ticking over by another bin by just - # .000000000005 (yeah...), which creates and off-by-one error - # further down the line. - # TODO: Find a better way to solve this, perhaps with linspace... - n_wave = int( - (self.meta["wave_max"].value - self.meta["wave_min"].value) / - from_currsys("!SIM.spectral.spectral_bin_width", self.cmds) - ) - fov_waveset = fov_waveset[:n_wave].to(u.um) - - # make canvas cube based on waveset of largest cube and NAXIS1,2 - # from fov.header - canvas_cube_hdu = fits.ImageHDU( - data=np.zeros((len(fov_waveset), - self.header["NAXIS2"], - self.header["NAXIS1"])), - header=self.header) + fov_waveset = self.cube_waveset.to(u.um) + + # Retain double precision when combining faint sources with sky and + # thermal backgrounds; subtraction can expose storage-rounding error. + shape = (len(fov_waveset), self.header['NAXIS2'], self.header['NAXIS1']) + canvas_cube_hdu = fits.ImageHDU(data=np.zeros(shape, dtype=np.float64), header=self.header) # set BUNIT initially to PHOTLAM / arcsec**2 canvas_cube_hdu.header["BUNIT"] = "ph cm-2 s-1 AA-1 arcsec-2" @@ -1112,26 +1095,38 @@ def make_hdu(self): assert canvas_cube_hdu.header["NAXIS3"] == self.detector_header["NAXIS3"] for field_hdu in self._make_cubefields(fov_waveset): - canvas_cube_hdu = imp_utils.add_imagehdu_to_imagehdu( - field_hdu, - canvas_cube_hdu, - spline_order=self.spline_order, - conserve_flux=True, - differential=True, - ) - - canvas_cube_hdu.data = sum( - self._make_imagefields(fov_waveset, self.spline_order), - start=canvas_cube_hdu.data, - ) - - for flux, x, y in self._make_tablefields(fov_waveset): - # To prevent adding array values in this manner. - assert not isinstance(x, Iterable), "x should be integer" - canvas_cube_hdu.data[:, y, x] += flux - - canvas_cube_hdu.data = sum(self._make_backfields(fov_waveset), - start=canvas_cube_hdu.data) + if spatial_psf is None: + canvas_cube_hdu.data += imp_utils.project_spatial_density(field_hdu, self.header) + else: + padded = spatial_psf.padded_header(self) + kernel = spatial_psf.get_kernel(self)[None, :, :] + px = (padded['NAXIS1'] - self.header['NAXIS1'])//2 + py = (padded['NAXIS2'] - self.header['NAXIS2'])//2 + # Bound the temporary padded cube, even for long input spectra. + for start in range(0, len(fov_waveset), 64): + chunk = fits.ImageHDU(data=field_hdu.data[start:start + 64], header=field_hdu.header) + data = imp_utils.project_spatial_density(chunk, padded) + data = fftconvolve(data, kernel, mode='same', axes=(-2, -1)) + canvas_cube_hdu.data[start:start + 64] += data[:, py:py + self.header['NAXIS2'], px:px + self.header['NAXIS1']] + + for image in self._make_imagefields(fov_waveset, self.spline_order, spatial_psf): + canvas_cube_hdu.data += image + + if spatial_psf is None: + for flux, x, y in self._make_tablefields(fov_waveset): + if 0 <= x < self.header['NAXIS1'] and 0 <= y < self.header['NAXIS2']: + canvas_cube_hdu.data[:, y, x] += flux + else: + for field in self._get_table_fields(): + spectra = {ref: average_spectrum(spec, fov_waveset).value for ref, spec in field.spectra.items()} + for row in field.field: + image = spatial_psf.point_source_image(self, row['x'], row['y'])/self.pixel_area.value + for start in range(0, len(fov_waveset), 64): + spectrum = spectra[row['ref']][start:start + 64] + canvas_cube_hdu.data[start:start + 64] += spectrum[:, None, None]*image[None, :, :]*row['weight'] + + for background in self._make_backfields(fov_waveset): + canvas_cube_hdu.data += background # SpectralTrace wants ph/s/um/arcsec2 --> get rid of m2, leave um canvas_cube_hdu.data *= self.area.to(u.cm ** 2).value diff --git a/scopesim/optics/fov_manager.py b/scopesim/optics/fov_manager.py index 4e7a894fe..261f48f86 100644 --- a/scopesim/optics/fov_manager.py +++ b/scopesim/optics/fov_manager.py @@ -189,7 +189,10 @@ def generate_fovs_list(self) -> Iterator[FieldOfView]: for vol in self.volumes_list: image_plane_id = vol["meta"].get("image_plane_id") det_eff = None - if image_plane_id is not None: + if image_plane_id is None and len(detector_effects) == 1: + # Single-plane fibre/imaging packages need no arm-routing key. + det_eff, = detector_effects + elif image_plane_id is not None: for candidate in detector_effects: if (from_currsys(candidate.meta["image_plane_id"], self.cmds) == image_plane_id): @@ -201,10 +204,9 @@ def generate_fovs_list(self) -> Iterator[FieldOfView]: "plate_scale": self.meta["plate_scale"], } if det_eff is not None: - scale_meta["pixel_scale"] = det_eff.meta.get( - "pixel_scale", scale_meta["pixel_scale"]) - scale_meta["plate_scale"] = det_eff.meta.get( - "plate_scale", scale_meta["plate_scale"]) + for key in ('pixel_scale', 'plate_scale'): + if det_eff.meta[key] != f'!INST.{key}': + scale_meta[key] = det_eff.meta[key] scales = from_currsys(scale_meta, self.cmds) pixel_scale = scales["pixel_scale"] plate_scale = scales["plate_scale"] @@ -215,14 +217,27 @@ def generate_fovs_list(self) -> Iterator[FieldOfView]: [ys_min, ys_max], pixel_scale=pixel_scale / 3600.) + if self.is_spectroscope and "eta_min" in vol["meta"]: + # Tile the physical slit exactly. Detector sampling sets the + # approximate integration step, never the admitted slit area. + for axis, lower, upper in ((1, xs_min, xs_max), (2, ys_min, ys_max)): + count = max(1, int(np.ceil((upper - lower)*3600/pixel_scale))) + skyhdr[f'NAXIS{axis}'] = count + skyhdr[f'CDELT{axis}'] = (upper - lower)/count + skyhdr[f'CRPIX{axis}'] = (count + 1)/2 + skyhdr[f'CRVAL{axis}'] = (lower + upper)/2 + fov_dethdr, _ = ipu.det_wcs_from_sky_wcs( WCS(skyhdr), pixel_scale, plate_scale) skyhdr.update(fov_dethdr.to_header()) - dethdr = fov_dethdr + dethdr = fov_dethdr.to_header() + dethdr.update(NAXIS=2, NAXIS1=skyhdr['NAXIS1'], NAXIS2=skyhdr['NAXIS2']) # useful for spectroscopy mode where slit dimensions is not the same # as detector dimensions - if decouple and det_eff is not None: + if decouple: + if det_eff is None: + raise ValueError('A decoupled FOV needs an image-plane ID or one unambiguous detector layout.') dethdr = det_eff.image_plane_header # TODO: Why is this .image_plane_header and not # .detector_headers()[0] or something? diff --git a/scopesim/optics/image_plane.py b/scopesim/optics/image_plane.py index 3ecba4d4f..f65c86f73 100644 --- a/scopesim/optics/image_plane.py +++ b/scopesim/optics/image_plane.py @@ -64,7 +64,7 @@ def __init__(self, header, cmds=None, **kwargs): raise ValueError(f"header must have a valid image-plane WCS: " f"{dict(header)}") - image = zeros_from_header(header, dtype=np.float32) + image = zeros_from_header(header, dtype=np.float64) self.hdu = fits.ImageHDU(data=image, header=header) self.hdu.header["BUNIT"] = "ph s-1" # photons per second (per pixel) diff --git a/scopesim/optics/image_plane_utils.py b/scopesim/optics/image_plane_utils.py index 15f91daf5..d64c995a7 100644 --- a/scopesim/optics/image_plane_utils.py +++ b/scopesim/optics/image_plane_utils.py @@ -7,6 +7,7 @@ from astropy.io import fits from astropy.utils.exceptions import AstropyWarning from scipy import ndimage as ndi +from scipy.sparse import csr_array from ..utils import get_logger @@ -625,6 +626,38 @@ def affine_map( return output +def project_spatial_density(image_hdu, canvas_header): + """Average a surface brightness image/cube over destination pixel areas. + + Overlaps are exact for aligned rectangular pixels, including signed and + unequal scales. Rotated inputs retain the existing affine interpolation + before this area integration. No global flux renormalization is applied + after clipping to the destination footprint. + """ + source = reorient_imagehdu(image_hdu.copy()) + source_wcs = WCS(source.header, naxis=2) + canvas_wcs = WCS(canvas_header, naxis=2) + nx, ny = source.data.shape[-1], source.data.shape[-2] + x = np.arange(canvas_header['NAXIS1'] + 1) - .5 + y = np.arange(canvas_header['NAXIS2'] + 1) - .5 + xedges, _ = source_wcs.all_world2pix(*canvas_wcs.all_pix2world(x, np.zeros_like(x), 0), 0) + _, yedges = source_wcs.all_world2pix(*canvas_wcs.all_pix2world(np.zeros_like(y), y, 0), 0) + + weights = [] + for count, edges in ((nx, xedges), (ny, yedges)): + lower = np.minimum(edges[:-1], edges[1:])[:, None] + upper = np.maximum(edges[:-1], edges[1:])[:, None] + input_lower = np.arange(count)[None, :] - .5 + overlap = np.maximum(0, np.minimum(upper, input_lower + 1) - np.maximum(lower, input_lower)) + weights.append(csr_array(overlap/(upper - lower))) + wx, wy = weights + planes = source.data.reshape(-1, ny, nx) + rows = wy @ planes.transpose(1, 0, 2).reshape(ny, -1) + rows = rows.reshape(wy.shape[0], len(planes), nx).transpose(1, 0, 2) + mapped = (wx @ rows.reshape(-1, nx).T).T + return mapped.reshape(source.data.shape[:-2] + (wy.shape[0], wx.shape[0])) + + def add_imagehdu_to_imagehdu( image_hdu: fits.ImageHDU, canvas_hdu: fits.ImageHDU, diff --git a/scopesim/optics/optical_train.py b/scopesim/optics/optical_train.py index a6031f2a7..9aa230d28 100644 --- a/scopesim/optics/optical_train.py +++ b/scopesim/optics/optical_train.py @@ -20,9 +20,12 @@ from .optics_manager import OpticsManager from .fov_manager import FOVManager from .image_plane import ImagePlane +from .fov import FieldOfView3D from ..commands.user_commands import UserCommands from ..detector import DetectorManager -from ..effects import ExtraFitsKeywords +from ..effects import ExtraFitsKeywords, SelectorWheel, MoffatPSF, TERCurve, SpectralEfficiency, EchelleSpectralEfficiency, SpectralTraceList +from ..effects.ter_curves import DichroicTree, FilterWheelBase, TaperedQuantumEfficiency +from ..effects.atmo_dispersion import ShiftFoV3D from ..utils import from_currsys, top_level_catch, get_logger from ..source.source_templates import empty_sky from .. import __version__ @@ -535,6 +538,15 @@ def observe(self, orig_source=None, update=True, **kwargs): for effect in self.optics_manager.source_effects: source = effect.apply_to(source) + # This source is a private copy, after all source effects. FOV extraction + # copies native-sampled spectra before FOV effects can modify them. Evaluate + # native sampling once for this observation; retain no cache on the + # source, spectrum, FOV or train, and never reuse it on the next call. + spectral_wavesets = {} + for spectrum in source.spectra.values(): + waveset = spectrum.waveset + spectral_wavesets[spectrum] = None if waveset is None else waveset.to_value(u.AA) + # [3D - Atmospheric shifts, PSF, NCPAs, Grating shift/distortion] # START OF MULTIPROCESSING @@ -562,14 +574,59 @@ def observe(self, orig_source=None, update=True, **kwargs): fovs = self.fov_manager.fovs nobar = len(fovs) <= 1 for fov in tqdm(fovs, desc=" FOVs", position=0, disable=nobar): - fov.extract_from(source) - - fov.view() - foveffs = self.optics_manager.fov_effects + # FOVManager can reuse a FOV whose previous observation is + # already mapped to detector coordinates. Prepare new sky + # illumination from its physical header, not that old HDU. + fov.hdu = None + foveffs = [] + for effect in self.optics_manager.fov_effects: + if isinstance(effect, SelectorWheel): + effect = effect.get_effect(fov.meta[effect.meta['selector_key']]) + if isinstance(effect, FilterWheelBase): + effect = effect.current_filter + if effect is not None: + foveffs.append(effect) + + for effect in foveffs: + # Slit shear needs spectral support before extraction. + # Fibre bundles have no single FOV trace and map separately. + if isinstance(effect, SpectralTraceList) and isinstance(fov, FieldOfView3D) and fov.trace_id is not None: + fov.meta['spectral_padding'] = effect.spectral_traces[fov.trace_id].spectral_padding(fov) + + # Render at the PSF stage when preceding effects only multiply + # spectra. Point positions then stay continuous, and resolved + # images are convolved before clipping to the aperture. + psf_indices = [i for i, effect in enumerate(foveffs) if isinstance(effect, MoffatPSF)] + defer_render = isinstance(fov, FieldOfView3D) and len(psf_indices) == 1 + if defer_render: + psf_index = psf_indices[0] + defer_render = all(isinstance(effect, (TERCurve, SpectralEfficiency, EchelleSpectralEfficiency, TaperedQuantumEfficiency)) for effect in foveffs[:psf_index]) + spatial_shifts = [] + if defer_render: + for effect in foveffs[psf_index + 1:]: + if isinstance(effect, ShiftFoV3D): + spatial_shifts.append(effect) + elif not isinstance(effect, (TERCurve, DichroicTree, SpectralEfficiency, EchelleSpectralEfficiency, TaperedQuantumEfficiency)): + break + # Wavelength-only multipliers commute with a displacement + # within each spectral plane; stop before other spatial effects. + point_sources = not (source.image_fields or source.cube_fields) + spatial_header = foveffs[psf_index].padded_header(fov, point_sources=point_sources) + for effect in spatial_shifts: + spatial_header = effect._padded_header(spatial_header, effect.get_shifts(fov)) + fov.extract_from(source, spatial_header=spatial_header, _spectral_wavesets=spectral_wavesets) + else: + fov.extract_from(source, _spectral_wavesets=spectral_wavesets) + fov.view() nobar = len(foveffs) <= 1 for effect in tqdm(foveffs, disable=nobar, desc=" FOV effects", position=1):#, leave=False): - fov = effect.apply_to(fov) + if effect in spatial_shifts: + continue # Applied to incoming illumination at the PSF stage. + if defer_render and effect is foveffs[psf_index] and spatial_shifts: + effect._apply_convolution(fov, spatial_shifts=spatial_shifts) + else: + fov = effect.apply_to(fov) if fov.hdu is None: logger.info(" Skipping empty FOV") @@ -584,6 +641,10 @@ def observe(self, orig_source=None, update=True, **kwargs): self.image_planes[fov.image_plane_id].header.update( WCS(fov.hdu).to_header()) + # Mapping is complete. Keep the rendered FOV and trace's + # rectified spectrum, without retaining every input cube. + fov.cube = None + # ..todo: finish off the multiple image plane stuff # END OF MULTIPROCESSING @@ -793,10 +854,12 @@ def write_header(self, hdulist, detector_array): # Image HDUs for hdu in hdulist[1:]: - hdu.header["EXPTIME"] = from_currsys("!OBS.exptime", self.cmds), "[s]" - hdu.header["DIT"] = from_currsys("!OBS.dit", self.cmds), "[s]" - hdu.header["NDIT"] = from_currsys("!OBS.ndit", self.cmds) - hdu.header["DETECTOR"] = detector_array._detector_list.meta["detector"] + # Arm-specific integration effects already recorded their actual + # exposure. Global command values are only for minimal setups. + for key in ('EXPTIME', 'DIT', 'NDIT'): + if key not in hdu.header: + hdu.header[key] = from_currsys(f'!OBS.{key.lower()}', self.cmds) + hdu.header['DETECTOR'] = detector_array._detector_list.display_name return hdulist def shutdown(self): diff --git a/scopesim/tests/tests_effects/test_AtmosphericDispersion.py b/scopesim/tests/tests_effects/test_AtmosphericDispersion.py index a578674a1..a72cc531f 100644 --- a/scopesim/tests/tests_effects/test_AtmosphericDispersion.py +++ b/scopesim/tests/tests_effects/test_AtmosphericDispersion.py @@ -3,12 +3,27 @@ import numpy as np +from scopesim.effects.atmo_dispersion import ADCShift + from scopesim.effects.shifts import AtmosphericDispersion, \ get_pixel_border_waves_from_atmo_disp from ..mocks.py_objects.yaml_objects import _atmo_yaml_dict +def test_repeated_integer_adc_shifts_preserve_positions_and_clip_edges(): + cube = np.zeros((10, 5, 7)) + cube[:, 2, 3] = 1 + dx = np.array([-1.2, -1.2, 0., 0., .9, .9, 8., -8., 0., 0.]) + dy = np.array([1., 1., 0., 0., -1., -1., 0., 0., 6., -6.]) + shifted = ADCShift.roll_shift(cube, dx, dy) + expected = np.zeros_like(cube) + expected[:2, 3, 2] = 1 + expected[2:4, 2, 3] = 1 + expected[4:6, 1, 4] = 1 + np.testing.assert_array_equal(shifted, expected) + + @pytest.fixture(scope="function") def atmo_yaml_dict(): atmo_dict = _atmo_yaml_dict() diff --git a/scopesim/tests/tests_effects/test_ConstPSF.py b/scopesim/tests/tests_effects/test_ConstPSF.py index 27706b5d9..cfbbd337f 100644 --- a/scopesim/tests/tests_effects/test_ConstPSF.py +++ b/scopesim/tests/tests_effects/test_ConstPSF.py @@ -16,9 +16,12 @@ """ import pytest +from types import SimpleNamespace import numpy as np from numpy import testing as npt +from astropy import units as u +from astropy.io import fits import matplotlib.pyplot as plt @@ -76,6 +79,37 @@ def test_initialised_with_psf_name_and_file_format(self, mock_path): class TestGetKernel: + def test_returned_kernel_does_not_alias_reference_or_next_call(self, const_psf, centre_fov): + centre_fov.view() + initial = const_psf.get_kernel(centre_fov) + expected = initial.copy() + initial.fill(0.) + npt.assert_array_equal(const_psf.get_kernel(centre_fov), expected) + + def test_reference_image_edits_reach_next_kernel(self, const_psf, centre_fov): + centre_fov.view() + first = const_psf.get_kernel(centre_fov).copy() + ext = const_psf.current_layer_id + const_psf._file[ext].data[:] = np.roll(const_psf._file[ext].data, 3, axis=0) + second = const_psf.get_kernel(centre_fov) + assert not np.allclose(first, second) + + def test_spectral_kernel_follows_current_grid_and_units(self, const_psf): + const_psf.cmds = {'!SIM.computing.psf_maxsize': 21} + header = fits.Header({'CDELT1': .004/3600, 'CDELT2': .004/3600, 'CUNIT1': 'deg', 'CUNIT2': 'deg', 'CDELT3': .01, 'CRPIX3': 1., 'CRVAL3': 1.55, 'CUNIT3': 'um'}) + first = fits.ImageHDU(np.zeros((3, 15, 15)), header) + fov = SimpleNamespace(header=first.header, hdu=first, wavelength=1.6*u.um) + initial = const_psf.get_kernel(fov).copy() + second = fits.ImageHDU(np.zeros((5, 15, 15)), header) + fov.header, fov.hdu = second.header, second + expanded = const_psf.get_kernel(fov).copy() + assert expanded.shape[0] == 5 + npt.assert_allclose(expanded[:3], initial) + second.header['CUNIT3'] = 'nm' + second.header['CDELT3'] *= 1000 + second.header['CRVAL3'] *= 1000 + npt.assert_allclose(const_psf.get_kernel(fov), expanded, rtol=1e-6) + @pytest.mark.parametrize("waves, max_pixel", [([1.1, 1.3], 1 / 3.), ([1.5, 1.7], 1 / 3.), diff --git a/scopesim/tests/tests_effects/test_MoffatPSF.py b/scopesim/tests/tests_effects/test_MoffatPSF.py index 1119df8dd..6fe78db4b 100644 --- a/scopesim/tests/tests_effects/test_MoffatPSF.py +++ b/scopesim/tests/tests_effects/test_MoffatPSF.py @@ -3,12 +3,40 @@ import numpy as np import pytest from astropy import units as u +from astropy.coordinates import AltAz, EarthLocation, SkyCoord +from astropy.table import Table +from astropy.time import Time +from astropy.utils import iers -from scopesim.effects.psfs.analytical import MoffatPSF +from scopesim.commands import UserCommands +from scopesim.effects.psfs.analytical import AOEnhanceablePSF, MoffatPSF + + +def test_seeing_cache_uses_physical_zenith_angle(monkeypatch): + from astropy.coordinates import SkyCoord, AltAz, EarthLocation + from astropy.time import Time + from scopesim.effects.psfs import analytical + + location = EarthLocation.from_geodetic(0*u.deg, 30*u.deg) + altitude = 60*u.deg + times = iter(Time(60310 + np.arange(10)/86400, format='mjd')) + + def observing_info(cmds): + time = next(times) + target = SkyCoord(alt=altitude, az=0*u.deg, frame=AltAz(obstime=time, location=location)) + return target, location, time, 'dark' + + monkeypatch.setattr(analytical, 'get_observation_info_from_cmds', observing_info) + psf = MoffatPSF(alpha=4.765, fwhm={'seeing': .7, 'seeing_unit': 'arcsec', 'pivot_wave': 500, 'pivot_wave_unit': 'nm'}) + first = psf.fwhm(np.array([.5])*u.um) + np.testing.assert_allclose(psf.fwhm(np.array([.5])*u.um), first) + altitude = 30*u.deg + expected = .7 / np.sin(altitude.to_value(u.rad))**.6 + assert psf.fwhm(np.array([.5])*u.um)[0].to_value(u.arcsec) == pytest.approx(expected) class SimpleFov: - header = {"CDELT1": 0.1 / 3600} + header = {"CDELT1": 0.1 / 3600, "CDELT2": 0.1 / 3600} meta = {"wave_min": 1.0 * u.um, "wave_max": 1.1 * u.um} waveset = np.linspace(1.0, 1.1, 5) * u.um @@ -16,7 +44,11 @@ class SimpleFov: def test_constant_fwhm_does_not_require_observation_commands(): psf = MoffatPSF(alpha=4.765, fwhm=0.7) - assert isinstance(psf, MoffatPSF) + wave = np.array([500., 1000.])*u.nm + np.testing.assert_allclose(psf.fwhm(wave).to_value(u.arcsec), [.7, .7]) + cached = psf._fwhm_interp[0] + np.testing.assert_allclose(psf.fwhm(wave.to(u.um)).to_value(u.arcsec), [.7, .7]) + assert psf._fwhm_interp[0] is cached def test_kernel_grows_to_flux_accuracy_without_renormalizing(): @@ -87,3 +119,135 @@ def test_legacy_renormalization_can_be_requested(): kernel = psf.get_kernel(SimpleFov()) assert np.sum(kernel) == pytest.approx(1.0) + + +@pytest.fixture +def seeing_cmds(monkeypatch): + monkeypatch.setattr(iers.conf, 'auto_download', False) + # The regression uses historical timestamps, independent of table freshness. + monkeypatch.setattr(iers.conf, 'auto_max_age', None) + monkeypatch.setattr(Time, 'now', lambda: Time('2026-10-18T10:50:25.709309')) + return UserCommands(properties={ + '!ATMO.longitude': -155.47833, '!ATMO.latitude': 19.82833, '!ATMO.altitude': 4160., + '!OBS.alt': 65.38002267134289, '!OBS.az': 0., '!OBS.airmass': None, + '!OBS.ra': None, '!OBS.dec': None, '!OBS.brightness': 'dark', + }) + + +@pytest.fixture +def natural_psf(seeing_cmds): + return MoffatPSF(cmds=seeing_cmds, alpha=4.765, + fwhm={'seeing': .65, 'seeing_unit': 'arcsec', 'pivot_wave': 500., 'pivot_wave_unit': 'nm'}) + + +def test_nighttime_timestamp_change_reuses_the_same_seeing(natural_psf, monkeypatch): + wave = np.array([500., 1000.])*u.nm + first = natural_psf.fwhm(wave) + cached = natural_psf._fwhm_interp[0] + # These are the two times in the performance notebook's original traceback. + monkeypatch.setattr(Time, 'now', lambda: Time('2026-10-18T10:50:25.735796')) + second = natural_psf.fwhm(wave.to(u.um)) + np.testing.assert_allclose(first.to_value(u.arcsec), [.68825435, .59916022], atol=1e-8, rtol=0) + np.testing.assert_array_equal(second, first) + assert natural_psf._fwhm_interp[0] is cached + + +def test_seeing_and_pivot_changes_refresh_the_cache(natural_psf): + wave = np.array([500., 1000.])*u.nm + first = natural_psf.fwhm(wave) + natural_psf.meta['fwhm']['seeing'] = 1.3 + np.testing.assert_allclose(natural_psf.fwhm(wave).to_value(u.arcsec), 2*first.to_value(u.arcsec)) + natural_psf.meta['fwhm']['pivot_wave'] = 1000. + np.testing.assert_allclose(natural_psf.fwhm(wave).to_value(u.arcsec), 2*2**.2*first.to_value(u.arcsec)) + + +def test_altitude_change_refreshes_the_seeing(natural_psf): + wave = np.array([500.])*u.nm + natural_psf.cmds['!OBS.alt'] = 90. + np.testing.assert_allclose(natural_psf.fwhm(wave).to_value(u.arcsec), [.65]) + natural_psf.cmds['!OBS.alt'] = 30. + np.testing.assert_allclose(natural_psf.fwhm(wave).to_value(u.arcsec), [.65*2**.6]) + + +def test_radec_observation_time_changes_the_seeing(natural_psf): + natural_psf.cmds['!OBS.ra'], natural_psf.cmds['!OBS.dec'] = 30., 20. + target = SkyCoord(ra=30*u.deg, dec=20*u.deg, distance=1*u.AU) + location = EarthLocation.from_geodetic(-155.47833*u.deg, 19.82833*u.deg, 4160*u.m) + widths = [] + for time in ['2026-10-18T10:50:25', '2026-10-18T12:50:25']: + natural_psf.cmds['!OBS.mjdobs'] = time + alt = target.transform_to(AltAz(obstime=Time(time), location=location)).alt + assert alt > 0*u.deg + width = natural_psf.fwhm(np.array([500.])*u.nm).to_value(u.arcsec)[0] + np.testing.assert_allclose(width, .65/np.sin(alt)**.6) + widths.append(width) + assert not np.isclose(widths[0], widths[1]) + + +def test_relative_ao_uses_the_same_natural_seeing_cache(seeing_cmds, tmp_path, monkeypatch): + path = tmp_path / 'ao.dat' + Table({'wavelength': [.4, .5, 1., 2.5], 'fwhm': [.5, .5, .5, .5]}, meta={'comments': ['wavelength_unit : um']}).write(path, format='ascii') + psf = AOEnhanceablePSF(cmds=seeing_cmds, ao_table=str(path), is_absolute=False, enable_ao=True, alpha=4.765, + fwhm={'seeing': .65, 'seeing_unit': 'arcsec', 'pivot_wave': 500., 'pivot_wave_unit': 'nm'}) + wave = np.array([500., 1000.])*u.nm + first = psf.fwhm(wave) + cached = psf._fwhm_interp[0] + monkeypatch.setattr(Time, 'now', lambda: Time('2026-10-18T10:50:25.735796')) + np.testing.assert_array_equal(psf.fwhm(wave), first) + assert psf._fwhm_interp[0] is cached + psf.meta['enable_ao'] = False + np.testing.assert_allclose(psf.fwhm(wave).to_value(u.arcsec), 2*first.to_value(u.arcsec)) + + +def test_table_fwhm_reuses_interpolator_without_observation_commands(tmp_path): + path = tmp_path / 'fwhm.dat' + Table({'wavelength': [.4, .5, 1., 2.5], 'fwhm': [.4, .5, 1., 2.5]}, meta={'comments': ['wavelength_unit : um']}).write(path, format='ascii') + psf = MoffatPSF(alpha=4.765, fwhm=str(path)) + wave = np.array([500., 1000.])*u.nm + np.testing.assert_allclose(psf.fwhm(wave).to_value(u.arcsec), [.5, 1.]) + cached = psf._fwhm_interp[0] + np.testing.assert_allclose(psf.fwhm(wave.to(u.um)).to_value(u.arcsec), [.5, 1.]) + assert psf._fwhm_interp[0] is cached + + +def test_returned_kernel_mutation_does_not_change_later_psf(): + psf = MoffatPSF(alpha=4.765, fwhm=.7) + first = psf.get_kernel(SimpleFov()) + expected = first.copy() + first.fill(0.) + np.testing.assert_array_equal(psf.get_kernel(SimpleFov()), expected) + + +def test_overwriting_fwhm_table_refreshes_interpolator(tmp_path): + path = tmp_path/'fwhm.dat' + table = Table({'wavelength': [.4, .5, 1., 2.5], 'fwhm': [.4, .5, 1., 2.5]}, meta={'comments': ['wavelength_unit : um']}) + table.write(path, format='ascii') + psf = MoffatPSF(alpha=4.765, fwhm=str(path)) + wave = np.array([.5, 1.])*u.um + first = psf.fwhm(wave) + table['fwhm'] *= 2 + table.write(path, format='ascii', overwrite=True) + np.testing.assert_allclose(psf.fwhm(wave), 2*first) + + +def test_ao_cache_tracks_resolved_path_contents_and_alpha(tmp_path): + first, second = tmp_path/'first.dat', tmp_path/'second.dat' + table = Table({'wavelength': [.4, .5, 1., 2.5], 'fwhm': [.3, .3, .3, .3]}, meta={'comments': ['wavelength_unit : um']}) + table.write(first, format='ascii') + table['fwhm'] *= 2 + table.write(second, format='ascii') + cmds = UserCommands(properties={'!INST.ao_file': str(first)}) + psf = AOEnhanceablePSF(cmds=cmds, ao_table='!INST.ao_file', is_absolute=True, enable_ao=True, alpha=4.765) + # PSF construction resolves command references. Also exercise a reference + # assigned directly to the configured metadata after construction. + psf.meta['ao_table'] = '!INST.ao_file' + wave = np.array([.5, 1.])*u.um + np.testing.assert_allclose(psf.fwhm(wave).value, .3) + assert psf.alpha == 4.765 + cmds['!INST.ao_file'] = str(second) + np.testing.assert_allclose(psf.fwhm(wave).value, .6) + table['fwhm'] /= 3 + table.write(second, format='ascii', overwrite=True) + np.testing.assert_allclose(psf.fwhm(wave).value, .2) + psf.meta['alpha'] = 3. + assert psf.alpha == 3. diff --git a/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py b/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py index 00114bcee..e263f9ca6 100644 --- a/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py +++ b/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py @@ -4,6 +4,20 @@ from astropy import units as u from astropy.io import fits from scopesim.effects.spectral_trace_list_utils import XiLamImage, SpectralTrace +from astropy.table import Table +from scopesim.optics.image_plane_utils import header_from_list_of_xy + + +def test_single_precision_cube_is_integrated_in_double_precision(): + fov = line_cube(.05) + fov.cube.data = np.zeros(fov.cube.data.shape, dtype=np.float32) + fov.cube.data[1000, :, :] = 1. + fov.cube.data[1000, 0, :] = 1e8 + image = XiLamImage(fov, lambda wave: .00001, slit_dispersion=lambda wave: .002) + # Integration across eta must retain the 19 faint cells beside the bright + # cell. The remapped line lies well inside the spectral boundaries. + expected = (1e8 + 19)*.05 + np.testing.assert_allclose(image.image.sum(axis=1), expected, rtol=1e-12) def line_cube(step, offset=0): @@ -85,3 +99,49 @@ def test_rectified_interpolation_reproduces_a_linear_density_field(): positions, wavelengths = np.array([.03, .17]), np.array([.998, 1.003]) expected = 2 + 3*positions + 5*wavelengths np.testing.assert_allclose(image.interp(positions, wavelengths, grid=False), expected, rtol=2e-6) + + +@pytest.mark.parametrize('change', ['wavelength_origin', 'wavelength_step', 'wavelength_pc', 'slit_origin', 'slit_step', 'slit_pc', 'slice_center', 'dispersion', 'source', 'units']) +def test_reused_slit_operator_matches_fresh_after_direct_changes(change): + fov = line_cube(.1) + cache = {} + slope = [.002] + dispersion = lambda wave: slope[0]*wave + XiLamImage(fov, lambda wave: .00001, slit_dispersion=dispersion, _cache=cache) + if change == 'wavelength_origin': fov.cube.header['CRVAL3'] += .001 + if change == 'wavelength_step': fov.cube.header['CDELT3'] *= 1.1 + if change == 'wavelength_pc': fov.cube.header['PC3_3'] = 1.1 + if change == 'slit_origin': fov.cube.header['CRVAL2'] += .04 + if change == 'slit_step': fov.cube.header['CDELT2'] *= 1.1 + if change == 'slit_pc': fov.cube.header['PC2_2'] = 1.1 + if change == 'slice_center': fov.meta['eta_min'] += .04*u.arcsec + if change == 'dispersion': slope[0] *= 1.3 + if change == 'source': fov.cube.data *= 2. + if change == 'units': + fov.cube.header['CUNIT2'], fov.cube.header['CUNIT3'] = 'deg', 'nm' + for keyword in ('CDELT2', 'CRVAL2'): fov.cube.header[keyword] /= 3600 + for keyword in ('CDELT3', 'CRVAL3'): fov.cube.header[keyword] *= 1000 + reused = XiLamImage(fov, lambda wave: .00001, slit_dispersion=dispersion, _cache=cache) + fresh = XiLamImage(fov, lambda wave: .00001, slit_dispersion=dispersion) + np.testing.assert_array_equal(reused.image, fresh.image) + + +@pytest.mark.parametrize('coordinate', ['xy2lam', 'xy2xi']) +def test_direct_optical_transform_edits_reach_next_projection(coordinate): + xi, wave = np.meshgrid(np.linspace(-1., 1., 5), np.linspace(.99, 1.01, 9)) + table = Table({'wavelength': wave.ravel()*u.um, 's': xi.ravel()*u.arcsec, 'x': (wave.ravel()-1)*1000*u.mm, 'y': xi.ravel()*.15*u.mm}) + reused, fresh = [SpectralTrace(table.copy(), trace_flux_jacobian=True) for _ in range(2)] + fov = line_cube(.1) + fov.cube.data *= np.array([1., 2., 4.])[None, None, :] + fov.header = fov.cube.header + fov.detector_header = header_from_list_of_xy([-10., 10.], [-2., 2.], .015, 'D') + fov.trace_id = 'linear' + fov.meta.update(wave_min=.995*u.um, wave_max=1.005*u.um, xi_min=-.05*u.arcsec, xi_max=.25*u.arcsec) + before = reused.map_spectra_to_focal_plane(fov).data.copy() + for trace in (reused, fresh): + if coordinate == 'xy2lam': trace.xy2lam.matrix[0, 0] += .0002 + if coordinate == 'xy2xi': trace.xy2xi.matrix[0, 0] += .04 + actual = reused.map_spectra_to_focal_plane(fov).data + expected = fresh.map_spectra_to_focal_plane(fov).data + np.testing.assert_array_equal(actual, expected) + assert np.max(abs(actual-before)) > .001*before.max() diff --git a/scopesim/tests/tests_effects/test_finite_spectral_cells.py b/scopesim/tests/tests_effects/test_finite_spectral_cells.py new file mode 100644 index 000000000..4f74407bc --- /dev/null +++ b/scopesim/tests/tests_effects/test_finite_spectral_cells.py @@ -0,0 +1,36 @@ +"""Photon integrals and line moments across finite slit and wavelength cells.""" +import numpy as np +import pytest +from scopesim.effects.spectral_trace_list_utils import _slit_cell_operator + + +@pytest.mark.parametrize('phase', [0., .25, .5, .75]) +@pytest.mark.parametrize('slope', [0., .02]) +def test_slit_remap_preserves_photons_and_illumination_centroid(phase, slope): + wave = np.linspace(.99, 1.01, 2001) + derivative = -.001 + slope * (wave - 1.) + step = wave[1] - wave[0] + eta, deta = .0825 + phase*step/abs(derivative[1000]), .03 + photons = np.zeros(len(wave)) + photons[1000] = 17. + output = _slit_cell_operator(wave, derivative, eta, deta) @ photons + assert output.sum() == pytest.approx(17., rel=1e-12) + center = np.sum(output * wave) / output.sum() + assert abs(center - (1. + eta*derivative[1000])) < .01*step + + +def test_zero_slit_dispersion_is_identity(): + wave = np.linspace(1., 1.1, 101) + photons = np.arange(len(wave)) + 1. + actual = _slit_cell_operator(wave, np.zeros_like(wave), .1, .1) @ photons + np.testing.assert_allclose(actual, photons, rtol=1e-12) + + +def test_finite_slit_cell_width_has_uniform_variance(): + wave = np.linspace(.99, 1.01, 10001) + derivative = np.full_like(wave, -.001) + photons = np.zeros(len(wave)) + photons[5000] = 1. + output = _slit_cell_operator(wave, derivative, 0., .1) @ photons + variance = np.sum(output * (wave - 1.)**2) + assert variance == pytest.approx((.001*.1)**2 / 12, rel=.002) diff --git a/scopesim/tests/tests_effects/test_shift_planes.py b/scopesim/tests/tests_effects/test_shift_planes.py new file mode 100644 index 000000000..5f4bcabe1 --- /dev/null +++ b/scopesim/tests/tests_effects/test_shift_planes.py @@ -0,0 +1,14 @@ +import numpy as np +from scipy.ndimage import map_coordinates +from scopesim.effects.atmo_dispersion import ShiftFoV3D + + +def test_dispersion_does_not_mix_wavelengths(): + data = np.zeros((19, 15, 17)) + data[:, 7, 8] = np.arange(19) + dx, dy = np.linspace(-.2, .2, 19), np.linspace(.4, .6, 19) + shifted = ShiftFoV3D.map_coords_shift(data, dx, dy) + yy, xx = np.indices(data.shape[1:]) + for wave in range(19): + expected = map_coordinates(data[wave], [yy-dy[wave], xx-dx[wave]], order=3, mode='nearest') + np.testing.assert_allclose(shifted[wave], expected, atol=2e-5) diff --git a/scopesim/tests/tests_effects/test_spatial_cube_shifts.py b/scopesim/tests/tests_effects/test_spatial_cube_shifts.py new file mode 100644 index 000000000..e8ef5a18b --- /dev/null +++ b/scopesim/tests/tests_effects/test_spatial_cube_shifts.py @@ -0,0 +1,30 @@ +import numpy as np +import pytest +from scipy.ndimage import map_coordinates +from scopesim.effects.atmo_dispersion import ShiftFoV3D + + +@pytest.mark.parametrize('shape', [(13, 5, 7), (9, 1, 4), (9, 8, 1)]) +def test_chunked_cubic_shift_matches_independent_spatial_planes(shape): + data = np.random.default_rng(127).uniform(size=shape) + dx = np.linspace(-2.3, 3.7, shape[0], dtype=np.float32) + dy = np.linspace(1.8, -.4, shape[0], dtype=np.float32) + yy, xx = np.indices(shape[1:], dtype=np.float32) + expected = np.stack([map_coordinates(plane, [yy - sy, xx - sx], order=3, mode='nearest') for plane, sx, sy in zip(data, dx, dy)]) + actual = ShiftFoV3D.map_coords_shift(data, dx, dy, chunk_w=4) + np.testing.assert_allclose(actual, expected, rtol=2e-6, atol=1e-6) + + +def test_spatial_shift_does_not_broaden_a_spectral_line(): + cube = np.zeros((11, 5, 7)) + cube[5, 2, 3] = 1. + shifted = ShiftFoV3D.map_coords_shift(cube, np.full(11, .3), np.full(11, -.2)) + np.testing.assert_array_equal(shifted[:5], 0.) + np.testing.assert_array_equal(shifted[6:], 0.) + + +@pytest.mark.parametrize('dx,dy', [(1, 0), (-1, 0), (0, 1), (0, -1), (2, -2)]) +def test_integer_roll_does_not_wrap_light_at_edges(dx, dy): + cube = np.ones((1, 5, 7)) + shifted = ShiftFoV3D.roll_shift(cube, [dx], [dy]) + assert shifted.sum() == (7 - abs(dx))*(5 - abs(dy)) diff --git a/scopesim/tests/tests_optics/test_FieldOfView.py b/scopesim/tests/tests_optics/test_FieldOfView.py index 61ced1915..bab799ad9 100644 --- a/scopesim/tests/tests_optics/test_FieldOfView.py +++ b/scopesim/tests/tests_optics/test_FieldOfView.py @@ -336,14 +336,18 @@ def test_make_spectrum_from_image(self): def test_make_spectrum_from_cube(self): src_cube = so._cube_source() # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm fov = _fov_197_202_um(ndim=1) + original = src_cube.fields[0].data.copy() fov.extract_from(src_cube) spec = fov.make_hdu() - in_sum = np.sum(src_cube.fields[0].data[74:77, :, :]) * 1e-8 + # 1.98 and 2.00 um cells are complete; the 2.02 um cell is half + # inside this FOV. Cropping must not change the input cube in place. + in_sum = (original[74:76].sum() + original[76].sum()/2)*1e-8 out_sum = np.sum(spec(fov.waveset).value) assert in_sum == approx(out_sum) + np.testing.assert_array_equal(src_cube.fields[0].data, original) def test_makes_spectrum_from_all_types_of_source_object(self): src_table = so._table_source() diff --git a/scopesim/tests/tests_optics/test_ImagePlane.py b/scopesim/tests/tests_optics/test_ImagePlane.py index c941e903e..2866f3451 100644 --- a/scopesim/tests/tests_optics/test_ImagePlane.py +++ b/scopesim/tests/tests_optics/test_ImagePlane.py @@ -284,6 +284,15 @@ def test_mm_images_on_large_canvas( class TestImagePlaneAdd: + def test_faint_signal_is_retained_on_bright_background(self, image_hdu_square): + plane = opt_imp.ImagePlane(image_hdu_square.header) + background = fits.ImageHDU(np.full(image_hdu_square.data.shape, 1e8), image_hdu_square.header) + signal = fits.ImageHDU(np.ones(image_hdu_square.data.shape), image_hdu_square.header) + plane.add(background) + before = plane.data.copy() + plane.add(signal) + np.testing.assert_allclose(plane.data - before, 1., atol=1e-7) + def test_simple_add_imagehdu_conserves_flux( self, image_hdu_square, image_hdu_rect, ): diff --git a/scopesim/tests/tests_optics/test_OpticsManager.py b/scopesim/tests/tests_optics/test_OpticsManager.py index 1c6b77cac..cd055f5d2 100644 --- a/scopesim/tests/tests_optics/test_OpticsManager.py +++ b/scopesim/tests/tests_optics/test_OpticsManager.py @@ -59,7 +59,7 @@ def test_effect_order_table_uses_runtime_phase_order(self, detector_setup_rows = tbl[tbl["phase"] == "detector_setup"] assert source_rows["effect"].tolist() == ["detector_qe_curve"] - assert fov_rows["effect"].tolist() == ["detector_qe_curve"] + assert len(fov_rows) == 0 # Ordinary QE is applied once, at the source. assert detector_setup_rows["effect"].tolist() == [ "micado_detector_geometry"] @@ -76,7 +76,7 @@ def test_effect_order_table_summarises_selector_wheels(self): "selector_key": "aperture_id", "wheel": [{ "selector_value": [1, 2], - "effect_class": "TERCurve", + "effect_class": "SpectralQuantumEfficiency", "effect_kwargs": {"filename": "TER_blank.dat"}, }], }, @@ -86,14 +86,14 @@ def test_effect_order_table_summarises_selector_wheels(self): opt_man = opt_mgr.OpticsManager([yaml_dict]) tbl = opt_man.effect_order_table() - source_row = tbl[tbl["phase"] == "source"][0] fov_row = tbl[tbl["phase"] == "fov"][0] - assert source_row["effect"] == "selected_qe" - assert source_row["class"] == "SelectorWheel" - assert source_row["selector_key"] == "aperture_id" - assert source_row["selector_values"] == "1, 2" - assert source_row["selected_classes"] == "TERCurve" + assert len(tbl[tbl['phase'] == 'source']) == 0 + assert fov_row['effect'] == 'selected_qe' + assert fov_row['class'] == 'SelectorWheel' + assert fov_row['selector_key'] == 'aperture_id' + assert fov_row['selector_values'] == '1, 2' + assert fov_row['selected_classes'] == 'SpectralQuantumEfficiency' assert fov_row["phase_z_orders"] == "610" diff --git a/scopesim/tests/tests_optics/test_adc_slit_illumination.py b/scopesim/tests/tests_optics/test_adc_slit_illumination.py new file mode 100644 index 000000000..4dae1e4d3 --- /dev/null +++ b/scopesim/tests/tests_optics/test_adc_slit_illumination.py @@ -0,0 +1,190 @@ +import numpy as np +import pytest +from astropy import units as u +from astropy.io import fits +from synphot import SourceSpectrum, ConstFlux1D +from synphot.units import PHOTLAM + +from scopesim import Source +from scopesim.effects.atmo_dispersion import ShiftFoV3D, _cube_wavelengths +from scopesim.effects.psfs.analytical import MoffatPSF +from scopesim.source.source_fields import BackgroundSourceField +from scopesim.optics.image_plane_utils import header_from_list_of_xy +from scopesim.tests.tests_optics.test_slit_spatial_integration import physical_fov + + +class Displacement(ShiftFoV3D): + def __init__(self, dx, dy, **kwargs): + super().__init__(**kwargs) + self.dx, self.dy = dx*u.arcsec, dy*u.arcsec + + def get_shifts(self, obj): + wave = _cube_wavelengths(obj) + return np.broadcast_to(self.dx, wave.shape, subok=True), np.broadcast_to(self.dy, wave.shape, subok=True) + + +def render(fov, source, psf, shifts): + support = psf.padded_header(fov, point_sources=not (source.image_fields or source.cube_fields)) + for effect in shifts: + support = effect._padded_header(support, effect.get_shifts(fov)) + fov.extract_from(source, spatial_header=support) + psf._apply_convolution(fov, spatial_shifts=shifts) + return fov.hdu.data*fov.pixel_area.value/1e4 + + +@pytest.mark.parametrize('width', [.33, .7, .9]) +@pytest.mark.parametrize('position,steps', [((0., 0.), (0, 1)), ((.037, .4), (0, -2)), ((5.6, .1), (-2, 1)), ((7., 0.), (-15, 0))]) +@pytest.mark.parametrize('sign', [-1, 1]) +def test_shifted_point_uses_incoming_psf_before_slit_and_ends(width, position, steps, sign): + fov = physical_fov(width) + fov.header['CDELT2'] *= sign + original = fov.header.copy() + cell = np.array([fov.header['CDELT1'], fov.header['CDELT2']])*3600 + displacement = np.array(steps)*cell + shift = Displacement(*displacement, method='roll') + psf = MoffatPSF(alpha=4.765, fwhm=.7) + source = Source(x=[position[0]], y=[position[1]], ref=[0], weight=[1.], spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + actual = render(fov, source, psf, [shift]) + expected = psf.point_source_image(fov, *(position + displacement)) + np.testing.assert_allclose(actual, np.broadcast_to(expected, actual.shape), rtol=2e-6, atol=1e-12) + for key in ('NAXIS1', 'NAXIS2', 'CRPIX1', 'CRPIX2', 'CDELT1', 'CDELT2', 'CRPIX1D', 'CRPIX2D'): + assert fov.header[key] == original[key] + assert fov.hdu.header[key] == original[key] + + +@pytest.mark.parametrize('method', ['roll', 'map_coords']) +def test_uniform_sky_is_not_lost_at_shifted_slit_edges(method): + fov = physical_fov(.7) + source = Source() + source.fields.append(BackgroundSourceField(field=None, header=fits.Header(), spectra={0: SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)})) + psf = MoffatPSF(alpha=4.765, fwhm=.7) + shift = Displacement(.23, -.17, method=method) + actual = render(fov, source, psf, [shift]) + np.testing.assert_allclose(actual, fov.pixel_area.value, rtol=1e-12) + + +def test_wavelength_dependent_shift_keeps_spectral_planes_independent(): + fov = physical_fov(.7) + count = len(fov.cube_waveset) + cell_y = fov.header['CDELT2']*3600 + displacement = np.resize([-1., 0., 1.], count)*cell_y + shift = Displacement(0., displacement, method='roll') + psf = MoffatPSF(alpha=4.765, fwhm=.7) + source = Source(x=[0], y=[.2], ref=[0], weight=[1.], spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + actual = render(fov, source, psf, [shift]) + expected = np.stack([psf.point_source_image(fov, 0., .2 + dy) for dy in displacement]) + np.testing.assert_allclose(actual, expected, rtol=1e-6) + + +def test_two_displacements_keep_light_that_reenters_slit(): + fov = physical_fov(.7) + step = fov.header['CDELT2']*3600 + shifts = [Displacement(0., 4*step, method='roll'), Displacement(0., -4*step, method='roll')] + psf = MoffatPSF(alpha=4.765, fwhm=.7) + source = Source(x=[0], y=[0], ref=[0], weight=[1.], spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + actual = render(fov, source, psf, shifts) + expected = psf.point_source_image(fov, 0., 0.) + np.testing.assert_allclose(actual, np.broadcast_to(expected, actual.shape), rtol=1e-6) + + +def test_shifted_resolved_image_and_cube_agree(): + header = header_from_list_of_xy([-1/3600, 1/3600], [-1/3600, 1/3600], .05/3600) + yy, xx = np.indices((header['NAXIS2'], header['NAXIS1'])) + density = np.exp(-((xx-20)**2 + (yy-28)**2)/8) + image_source = Source(image_hdu=fits.ImageHDU(data=density*.05**2, header=header), spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + cube_header = header.copy() + cube_header.update(CRPIX3=1., CRVAL3=.49, CDELT3=.01, CUNIT3='um', CTYPE3='WAVE', BUNIT='photlam arcsec-2') + cube_source = Source(cube=fits.ImageHDU(data=np.broadcast_to(density, (13, *density.shape)).copy(), header=cube_header)) + psf = MoffatPSF(alpha=4.765, fwhm=.7) + shift = Displacement(.15, -.2, method='roll') + images = [render(physical_fov(.7), source, psf, [shift]) for source in (image_source, cube_source)] + np.testing.assert_allclose(images[0], images[1], rtol=2e-6, atol=1e-10) + fov = physical_fov(.7) + cell = np.array([fov.header['CDELT1'], fov.header['CDELT2']])*3600 + displacement = np.rint(np.round(np.array([.15, -.2])/cell, 1))*cell + moved = image_source.make_copy() + moved.shift(*displacement) + expected = render(fov, moved, psf, []) + np.testing.assert_allclose(images[0], expected, rtol=2e-6, atol=1e-10) + + +@pytest.mark.parametrize('method', ['roll', 'map_coords']) +def test_zero_shift_reproduces_unshifted_render(method): + psf = MoffatPSF(alpha=4.765, fwhm=.7) + source = Source(x=[.1], y=[.2], ref=[0], weight=[1.], spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + unshifted = render(physical_fov(.7), source, psf, []) + shift = Displacement(0., 0., method=method) + fov = physical_fov(.7) + assert shift._padded_header(fov.header, shift.get_shifts(fov)) == fov.header + shifted = render(fov, source, psf, [shift]) + np.testing.assert_allclose(unshifted, shifted, rtol=1e-14, atol=1e-14) + + +@pytest.mark.parametrize('position', [(0., 0.), (.1, .4), (5.6, 0.)]) +def test_fractional_shift_agrees_with_pre_slit_psf_flux_and_centroid(position): + fov = physical_fov(.7) + cell = np.array([fov.header['CDELT1'], fov.header['CDELT2']])*3600 + shift = Displacement(*(np.array([.3, -.4])*cell), method='map_coords') + psf = MoffatPSF(alpha=4.765, fwhm=.7) + source = Source(x=[position[0]], y=[position[1]], ref=[0], weight=[1.], spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + actual = render(fov, source, psf, [shift])[0] + expected = psf.point_source_image(fov, *(position + np.array([.3, -.4])*cell)) + assert actual.sum() == pytest.approx(expected.sum(), rel=1e-3) + yy, xx = np.indices(actual.shape) + for coordinate in (xx, yy): + measured = np.sum(actual*coordinate)/actual.sum() + reference = np.sum(expected*coordinate)/expected.sum() + assert measured == pytest.approx(reference, abs=.01) + + +@pytest.mark.parametrize('intervening_throughput', [False, True]) +def test_observe_applies_shift_before_slit_through_wavelength_throughput(intervening_throughput): + from types import SimpleNamespace + from scopesim import OpticalTrain + from scopesim.effects import TERCurve + + psf = MoffatPSF(alpha=4.765, fwhm=.7) + images = [] + for move_source in (False, True): + fov = physical_fov(.7) + displacement = -2*fov.header['CDELT2']*3600 + source_y = .4 + displacement if move_source else .4 + source = Source(x=[0], y=[source_y], ref=[0], weight=[1.], spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + effects = [psf] + if intervening_throughput: + effects.append(TERCurve(array_dict={'wavelength': [.4, .7], 'transmission': [.2, .8]}, wavelength_unit='um')) + if not move_source: + effects.append(Displacement(0., displacement, method='roll')) + train = OpticalTrain() + train.cmds = fov.cmds + train.fov_manager = SimpleNamespace(fovs=[fov]) + train.optics_manager = SimpleNamespace(source_effects=[], fov_effects=effects, image_plane_effects=[]) + train.image_planes = [SimpleNamespace(add=lambda hdu, **kwargs: images.append(hdu.data.copy()))] + train.observe(source, update=False) + train.observe(source, update=False) + for image in images[1:]: + np.testing.assert_allclose(images[0], image, rtol=2e-6, atol=1e-12) + + +@pytest.mark.parametrize('method', ['roll', 'map_coords']) +@pytest.mark.parametrize('fraction', [1e-12, .04, .54, .56]) +def test_final_padding_uses_applied_displacement_but_source_support_is_conservative(method, fraction): + fov = physical_fov(.7) + cell_y = fov.header['CDELT2']*3600 + shift = Displacement(1e-14, fraction*cell_y, method=method) + shifts = shift.get_shifts(fov) + source_support = shift._padded_header(fov.header, shifts) + final_support = shift._padded_header(fov.header, shifts, rounded=True) + assert source_support['NAXIS2'] > fov.header['NAXIS2'] + assert final_support['NAXIS1'] == fov.header['NAXIS1'] + applied = np.round(fraction, 1) + if method == 'roll': + applied = np.round(np.float32(applied)) + padding = int(np.ceil(abs(applied))) + (4 if method == 'map_coords' else 0) if applied else 0 + assert final_support['NAXIS2'] == fov.header['NAXIS2'] + 2*padding + if applied == 0: + psf = MoffatPSF(alpha=4.765, fwhm=.7) + source = Source(x=[.1], y=[.2], ref=[0], weight=[1.], spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + actual = render(fov, source, psf, [shift]) + expected = render(physical_fov(.7), source, psf, []) + np.testing.assert_array_equal(actual, expected) diff --git a/scopesim/tests/tests_optics/test_slit_spatial_integration.py b/scopesim/tests/tests_optics/test_slit_spatial_integration.py new file mode 100644 index 000000000..2a12e4ba5 --- /dev/null +++ b/scopesim/tests/tests_optics/test_slit_spatial_integration.py @@ -0,0 +1,146 @@ +from types import SimpleNamespace +import numpy as np +import pytest +from astropy import units as u +from astropy.wcs import WCS +from astropy.io import fits +from synphot import SourceSpectrum, ConstFlux1D +from synphot.units import PHOTLAM +from scipy.special import stdtr +from scopesim import Source, UserCommands +from scopesim.effects import TERCurve, DetectorWindow +from scopesim.effects.psfs.analytical import MoffatPSF +from scopesim.optics.fov import FieldOfView3D +from scopesim.optics.fov_manager import FOVManager +from scopesim.optics.image_plane_utils import header_from_list_of_xy, det_wcs_from_sky_wcs, project_spatial_density + + +def slit_fov(width, step=.1570545): + nx, ny = int(np.ceil(11/step)), int(np.ceil(width/step)) + header = header_from_list_of_xy([-5.5/3600, 5.5/3600], [-width/7200, width/7200], step/3600) + header.update(NAXIS1=nx, NAXIS2=ny, CRPIX1=(nx+1)/2, CRPIX2=(ny+1)/2, CRVAL1=0., CRVAL2=0., CDELT1=11/nx/3600, CDELT2=width/ny/3600) + return SimpleNamespace(header=header, waveset=np.linspace(.5, .6, 16)*u.um) + + +@pytest.mark.parametrize('width', [.30, .33, .70]) +@pytest.mark.parametrize('fwhm', [.10, .14, .7]) +@pytest.mark.parametrize('offset', [0, .13, .4]) +def test_integrated_point_psf_matches_analytic_slit_loss(width, fwhm, offset): + alpha = 4.765 + psf = MoffatPSF(alpha=alpha, fwhm=fwhm, rounded_edges=False) + image = psf.point_source_image(slit_fov(width), .037, offset) + gamma = fwhm/(2*np.sqrt(2**(1/alpha) - 1)) + scale = np.sqrt(2*alpha - 2)/gamma + expected = stdtr(2*alpha - 2, (width/2 - offset)*scale) - stdtr(2*alpha - 2, (-width/2 - offset)*scale) + # The analytic comparison has an infinite slit length; the model has 11 arcsec. + assert image.sum() == pytest.approx(expected, abs=2e-5) + + +def test_off_slit_source_remains_in_psf_extraction_support(): + fov = slit_fov(.33) + psf = MoffatPSF(alpha=4.765, fwhm=.7) + padded = psf.padded_header(fov) + x, y = WCS(padded).all_world2pix(0, .4/3600, 0) + assert -.5 < x < padded['NAXIS1'] - .5 + assert -.5 < y < padded['NAXIS2'] - .5 + assert psf.point_source_image(fov, 0, .4).sum() > .1 + + +def test_kernel_cache_responds_to_seeing_and_pixel_scale(): + fov = slit_fov(.33) + psf = MoffatPSF(alpha=4.765, fwhm=.7, rounded_edges=False) + first = psf.get_kernel(fov) + np.testing.assert_array_equal(psf.get_kernel(fov), first) + psf.meta['fwhm'] = .2 + narrower = psf.get_kernel(fov) + assert narrower.max() > first.max() + fov.header['CDELT2'] /= 2 + finer = psf.get_kernel(fov) + assert finer.max() < narrower.max() + + +def physical_fov(width=.33): + header = slit_fov(width).header + detector_wcs, _ = det_wcs_from_sky_wcs(WCS(header), .1570545, 10.) + header.update(detector_wcs.to_header()) + cmds = UserCommands(properties={'!SIM.spectral.spectral_bin_width': .01}) + return FieldOfView3D(header, [.5, .6], cmds=cmds, area=1*u.cm**2) + + +def test_throughput_then_psf_admits_off_slit_point_light(): + fov = physical_fov() + psf = MoffatPSF(alpha=4.765, fwhm=.7, rounded_edges=False) + source = Source(x=[.037], y=[.4], ref=[0], weight=[1.], spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + fov.extract_from(source, spatial_header=psf.padded_header(fov)) + transmission = TERCurve(array_dict={'wavelength': [.4, .7], 'transmission': [.5, .5]}, wavelength_unit='um') + transmission.apply_to(fov) + psf.apply_to(fov) + expected = psf.point_source_image(fov, .037, .4).sum()*.5 + actual = fov.hdu.data.sum(axis=(1, 2))*fov.pixel_area.value/1e4 + np.testing.assert_allclose(actual, expected, rtol=1e-6) + + +@pytest.mark.parametrize('sign', [-1, 1]) +def test_surface_brightness_uses_signed_wcs_and_pixel_overlaps(sign): + header = header_from_list_of_xy([-.8/3600, .8/3600], [-.6/3600, .6/3600], .1/3600) + header['CDELT1'] *= sign + field = fits.ImageHDU(data=np.full((header['NAXIS2'], header['NAXIS1']), 7.), header=header) + target = slit_fov(.33).header + projected = project_spatial_density(field, target) + # Only the small input's length illuminates this 11-arcsec slit. + expected_area = header['NAXIS1']*.1*.33 + assert projected.sum()*abs(target['CDELT1']*target['CDELT2'])*3600**2 == pytest.approx(7*expected_area, rel=1e-6) + + +def test_resolved_image_and_cube_have_same_psf_and_spectral_normalization(): + header = header_from_list_of_xy([-1/3600, 1/3600], [-1/3600, 1/3600], .05/3600) + yy, xx = np.indices((header['NAXIS2'], header['NAXIS1'])) + density = np.exp(-((xx - 20)**2 + (yy - 26)**2)/8) + image = fits.ImageHDU(data=density*.05**2, header=header) + spectrum = SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM) + image_source = Source(image_hdu=image, spectra=[spectrum]) + cube_header = header.copy() + cube_header.update(CRPIX3=1., CRVAL3=.49, CDELT3=.01, CUNIT3='um', CTYPE3='WAVE', BUNIT='photlam arcsec-2') + cube = fits.ImageHDU(data=np.broadcast_to(density, (13, *density.shape)).copy(), header=cube_header) + cube_source = Source(cube=cube) + psf = MoffatPSF(alpha=4.765, fwhm=.7, rounded_edges=False) + rendered = [] + for source in (image_source, cube_source): + fov = physical_fov() + fov.extract_from(source, spatial_header=psf.padded_header(fov)) + psf.apply_to(fov) + rendered.append(fov.hdu.data) + np.testing.assert_allclose(rendered[0], rendered[1], rtol=2e-6, atol=1e-10) + + +@pytest.mark.parametrize('width', [.30, .33, .70]) +def test_fov_manager_tiles_physical_slit_width(width): + class SlitVolume: + def apply_to(self, volumes, **kwargs): + for volume in volumes: + volume.update(x_min=-5.5, x_max=5.5, y_min=-width/2, y_max=width/2) + volume['meta'].update(eta_min=-width/2*u.arcsec, eta_max=width/2*u.arcsec, image_plane_id=0) + return volumes + + detector = DetectorWindow(pixel_size=.015, x=0, y=0, width=100, height=100, units='pixel', image_plane_id=0, pixel_scale=.1570545, plate_scale=10.) + manager = FOVManager([SlitVolume(), detector], pixel_scale=.1570545, plate_scale=10., preload_fovs=False) + manager.is_spectroscope = True + fov = next(manager.generate_fovs_list()) + assert fov.header['NAXIS2']*fov.header['CDELT2']*3600 == pytest.approx(width) + assert fov.header['NAXIS1']*fov.header['CDELT1']*3600 == pytest.approx(11.) + + +@pytest.mark.parametrize('edge', ['left', 'right']) +def test_subpixel_deposition_does_not_wrap_or_index_beyond_slit(edge): + fov = physical_fov() + fov.meta['sub_pixel'] = True + xpixel = -.25 if edge == 'left' else fov.header['NAXIS1'] - .75 + ypixel = (fov.header['NAXIS2'] - 1)/2 + x, y = WCS(fov.header, naxis=2).all_pix2world(xpixel, ypixel, 0) + source = Source(x=[x*3600], y=[y*3600], ref=[0], weight=[1.], spectra=[SourceSpectrum(ConstFlux1D, amplitude=1*PHOTLAM)]) + fov.extract_from(source) + cube = fov.make_hdu() + admitted = cube.data.sum(axis=(1, 2))*fov.pixel_area.value/1e4 + np.testing.assert_allclose(admitted, .75, rtol=1e-6) + opposite = -1 if edge == 'left' else 0 + np.testing.assert_array_equal(cube.data[:, :, opposite], 0.) From fd14b0906d5dc0e3693221cfa1a3f5cd7a8afc13 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:55:53 -0700 Subject: [PATCH 7/7] Apply imaging branch throughput before wavelength integration Source commits: dbbba294, 38bf6c59, efa913a2, 71321def, b9ba36e6 --- scopesim/effects/selector_wheel.py | 11 +- scopesim/effects/ter_curves.py | 13 +- scopesim/optics/optical_train.py | 26 +++- .../tests/tests_effects/test_MoffatPSF.py | 31 +++++ .../tests_effects/test_imaging_branches.py | 115 ++++++++++++++++++ .../test_adc_slit_illumination.py | 2 +- 6 files changed, 186 insertions(+), 12 deletions(-) create mode 100644 scopesim/tests/tests_effects/test_imaging_branches.py diff --git a/scopesim/effects/selector_wheel.py b/scopesim/effects/selector_wheel.py index 333635bb1..c58692572 100644 --- a/scopesim/effects/selector_wheel.py +++ b/scopesim/effects/selector_wheel.py @@ -19,6 +19,7 @@ from ..optics.fov import FieldOfView from ..optics.image_plane import ImagePlane from ..detector.detector import Detector +from ..source.source import Source logger = get_logger(__name__) @@ -80,8 +81,14 @@ def __init__(self, **kwargs): self.z_order = self._resolve_z_order() - def apply_to(self, obj, **kwargs): - """Based on the selector_key's value in obj.meta, apply the corresponding effect from the wheel.""" + def apply_to(self, obj, fov=None, **kwargs): + """Select from obj.meta, or fov.meta when applying branch effects to a Source.""" + if isinstance(obj, Source) and fov is not None: + selector_value = fov.meta[self.meta["selector_key"]] + effect_to_apply = self.get_effect(selector_value) + if effect_to_apply is not None: + obj = effect_to_apply.apply_to(obj, fov=fov, **kwargs) + if isinstance(obj, FieldOfView): if self.meta['selector_key'] not in obj.meta.keys(): raise ValueError(f"Selector key {self.meta['selector_key']} not found in FieldOfView meta.") diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py index e35d2129f..f4a6e4219 100644 --- a/scopesim/effects/ter_curves.py +++ b/scopesim/effects/ter_curves.py @@ -4,6 +4,7 @@ import warnings from typing import ClassVar from collections.abc import Collection, Iterable +from numbers import Integral import numpy as np from numpy.typing import ArrayLike, NDArray from scipy import integrate @@ -1564,6 +1565,7 @@ class DichroicTree(Effect): ``` This effect creates a volume of FoVs, one for each row in the dichroic tree table. + For imaging, pass a Source with its fov context to apply the branch throughput before rendering. """ required_keys = {"tree_filename", "dichroic_names", "filename_format"} z_order: ClassVar[tuple[int, ...]] = (117, 217, 617) @@ -1571,6 +1573,9 @@ class DichroicTree(Effect): def __init__(self, cmds=None, **kwargs): check_keys(kwargs, self.required_keys, action="error") super().__init__(cmds=cmds, **kwargs) + if "z_order" in self.meta: + z_order = self.meta["z_order"] + self.z_order = (int(z_order),) if isinstance(z_order, Integral) else tuple(z_order) # Load all the dichroic TER curves self.dichroics = {} @@ -1587,7 +1592,7 @@ def __init__(self, cmds=None, **kwargs): tree = DataContainer(filename=tree_path) self.table = tree.table - def apply_to(self, obj, **kwargs): + def apply_to(self, obj, fov=None, **kwargs): """ To the object, loop through each row of the table and apply the dichroics (defined by columns) with the corresponding action defined in the table cell. @@ -1625,9 +1630,11 @@ def apply_to(self, obj, **kwargs): obj.volumes = new_vols # 2. During observe - if isinstance(obj, FieldOfView): + if isinstance(obj, FieldOfView) or (isinstance(obj, Source) and fov is not None): logger.debug("Executing %s, observe", self.meta['name']) - obj_id = obj.meta.get(id_colname, None) + branch_fov = obj if isinstance(obj, FieldOfView) else fov + obj_id = branch_fov.meta[id_colname] + setattr(branch_fov, id_colname, obj_id) row = self.table[self.table[id_colname] == obj_id] if len(row) != 1: diff --git a/scopesim/optics/optical_train.py b/scopesim/optics/optical_train.py index 9aa230d28..ebad53f79 100644 --- a/scopesim/optics/optical_train.py +++ b/scopesim/optics/optical_train.py @@ -38,8 +38,14 @@ N_PROCESSES = mp.cpu_count() - 1 USE_MULTIPROCESSING = False -def extract_source(fov, source): - fov.extract_from(source) +def extract_source(fov, source, source_effects=(), spatial_header=None, _spectral_wavesets=None): + if source_effects: + source = source.make_copy() + for effect in source_effects: + source = effect.apply_to(source, fov=fov) + # Branch filtering changes the spectrum and its native wavelength grid. + _spectral_wavesets = None + fov.extract_from(source, spatial_header=spatial_header, _spectral_wavesets=_spectral_wavesets) return fov def view_fov(fov, hdu_type): @@ -535,7 +541,15 @@ def observe(self, orig_source=None, update=True, **kwargs): source = self.prepare_source(source) # [1D - transmission curves] - for effect in self.optics_manager.source_effects: + source_effects = self.optics_manager.source_effects + branch_effects = [] + if not self.fov_manager.is_spectroscope: + for ii, effect in enumerate(source_effects): + if isinstance(effect, (DichroicTree, SelectorWheel)): + # Attenuation and emission downstream of a split belong to each branch. + source_effects, branch_effects = source_effects[:ii], source_effects[ii:] + break + for effect in source_effects: source = effect.apply_to(source) # This source is a private copy, after all source effects. FOV extraction @@ -558,7 +572,7 @@ def observe(self, orig_source=None, update=True, **kwargs): with mp.Pool(processes=N_PROCESSES) as pool: fovs = pool.starmap(extract_source, - zip(fovs, [source] * len(fovs))) + zip(fovs, [source] * len(fovs), [branch_effects] * len(fovs))) with mp.Pool(processes=N_PROCESSES) as pool: fovs = pool.starmap(view_fov, @@ -614,9 +628,9 @@ def observe(self, orig_source=None, update=True, **kwargs): spatial_header = foveffs[psf_index].padded_header(fov, point_sources=point_sources) for effect in spatial_shifts: spatial_header = effect._padded_header(spatial_header, effect.get_shifts(fov)) - fov.extract_from(source, spatial_header=spatial_header, _spectral_wavesets=spectral_wavesets) + extract_source(fov, source, branch_effects, spatial_header=spatial_header, _spectral_wavesets=spectral_wavesets) else: - fov.extract_from(source, _spectral_wavesets=spectral_wavesets) + extract_source(fov, source, branch_effects, _spectral_wavesets=spectral_wavesets) fov.view() nobar = len(foveffs) <= 1 for effect in tqdm(foveffs, disable=nobar, diff --git a/scopesim/tests/tests_effects/test_MoffatPSF.py b/scopesim/tests/tests_effects/test_MoffatPSF.py index 6fe78db4b..03b7bbccf 100644 --- a/scopesim/tests/tests_effects/test_MoffatPSF.py +++ b/scopesim/tests/tests_effects/test_MoffatPSF.py @@ -1,4 +1,5 @@ import logging +from unittest.mock import Mock import numpy as np import pytest @@ -251,3 +252,33 @@ def test_ao_cache_tracks_resolved_path_contents_and_alpha(tmp_path): np.testing.assert_allclose(psf.fwhm(wave).value, .2) psf.meta['alpha'] = 3. assert psf.alpha == 3. + + +@pytest.mark.parametrize("psf_class", [MoffatPSF, AOEnhanceablePSF]) +def test_observation_configuration_is_only_required_when_evaluating(psf_class, monkeypatch): + observation_info = Mock(side_effect=KeyError("!OBS.airmass")) + monkeypatch.setattr("scopesim.effects.psfs.analytical.get_observation_info_from_cmds", observation_info) + kwargs = {"alpha": 4.765, "fwhm": {"seeing": 0.7, "seeing_unit": "arcsec", "pivot_wave": 500, "pivot_wave_unit": "nm"}} + if psf_class is AOEnhanceablePSF: + kwargs.update(ao_table="unused.dat", enable_ao=False, is_absolute=False) + psf = psf_class(**kwargs) + observation_info.assert_not_called() + with pytest.raises(KeyError, match="airmass"): + psf.fwhm(np.array([0.5]) * u.um) + + + +def test_ao_toggle_preserves_natural_and_absolute_widths(tmp_path): + filename = tmp_path / "ao.ecsv" + table = Table({"wavelength": np.array([0.5, 1.0, 1.5, 2.0]) * u.um, "fwhm": [0.5, 0.5, 0.5, 0.5]}) + table.write(filename, format="ascii.ecsv") + psf = AOEnhanceablePSF(alpha=4.765, fwhm=0.7, ao_table=str(filename), enable_ao=False, is_absolute=False) + wavelengths = np.array([1.0]) * u.um + np.testing.assert_allclose(psf.fwhm(wavelengths).to_value(u.arcsec), [0.7]) + psf.meta["enable_ao"] = True + np.testing.assert_allclose(psf.fwhm(wavelengths).to_value(u.arcsec), [0.35]) + psf.meta["is_absolute"] = True + np.testing.assert_allclose(psf.fwhm(wavelengths).to_value(u.arcsec), [0.5]) + psf.meta["fwhm"] = 0.9 + psf.meta["enable_ao"] = False + np.testing.assert_allclose(psf.fwhm(wavelengths).to_value(u.arcsec), [0.9]) diff --git a/scopesim/tests/tests_effects/test_imaging_branches.py b/scopesim/tests/tests_effects/test_imaging_branches.py new file mode 100644 index 000000000..ce308f392 --- /dev/null +++ b/scopesim/tests/tests_effects/test_imaging_branches.py @@ -0,0 +1,115 @@ +"""Branch throughput must be applied before integrating an imaging spectrum.""" + +from types import SimpleNamespace + +import numpy as np +import pytest + +from scopesim.commands import UserCommands +from scopesim.effects import DichroicTree, SelectorWheel, TERCurve +from scopesim.optics.image_plane import ImagePlane +from scopesim.optics.optical_train import OpticalTrain, extract_source +from scopesim.source.source import Source +from scopesim.tests.mocks.py_objects.fov_objects import _centre_fov +from scopesim.tests.mocks.py_objects.source_objects import _table_source + + +@pytest.fixture +def tree(tmp_path): + (tmp_path / 'tree.dat').write_text('image_plane_id first second\n0 R X\n1 T R\n2 T T\n') + for name, transmission in [('first', 0.75), ('second', 0.6)]: + curve = '# wavelength_unit: um\nwavelength transmission reflection emissivity\n' + curve += f'0.5 {transmission} {1 - transmission} 0\n2.5 {transmission} {1 - transmission} 0\n' + (tmp_path / f'{name}.dat').write_text(curve) + return DichroicTree(name='tree', tree_filename=str(tmp_path / 'tree.dat'), dichroic_names=['first', 'second'], filename_format=str(tmp_path / '{}.dat')) + + +@pytest.fixture +def selector(): + wheel = [] + for arm, transmission in enumerate([0.25, 0.75]): + curve = {'wavelength': [0.5, 2.5], 'transmission': [transmission, transmission], 'emissivity': [0, 0]} + wheel.append({'selector_value': arm, 'effect_class': 'QuantumEfficiencyCurve', 'effect_kwargs': {'array_dict': curve, 'wavelength_unit': 'um'}}) + return SelectorWheel(selector_key='image_plane_id', wheel=wheel, allowed_missing_selector_values=[2]) + + +@pytest.fixture +def source(): + source = _table_source() + background = TERCurve(array_dict={'wavelength': [0.5, 2.5], 'transmission': [1, 1], 'emission': [1e-8, 1e-8]}, wavelength_unit='um', emission_unit='photlam') + source.append(background.background_source[0]) + return source + + +@pytest.mark.parametrize('effect_name, transmissions', [('tree', [0.25, 0.3, 0.45]), ('selector', [0.25, 0.75, 1.0])]) +def test_branch_attenuates_coloured_sources_and_background_independently(request, effect_name, transmissions, source): + effect = request.getfixturevalue(effect_name) + reference = extract_source(_centre_fov(n=10), source).make_hdu().data + for arm, transmission in enumerate(transmissions): + fov = _centre_fov(n=10) + fov.meta['image_plane_id'] = arm + extract_source(fov, source, [effect]) + np.testing.assert_allclose(fov.make_hdu().data, transmission * reference) + # Rendering each branch must leave the shared source, including its sky, intact. + np.testing.assert_array_equal(extract_source(_centre_fov(n=10), source).make_hdu().data, reference) + + +@pytest.mark.parametrize('effect_name', ['tree', 'selector']) +def test_source_without_fov_context_is_unchanged(request, effect_name, source): + effect = request.getfixturevalue(effect_name) + reference = extract_source(_centre_fov(n=10), source).make_hdu().data + effect.apply_to(source) + np.testing.assert_array_equal(extract_source(_centre_fov(n=10), source).make_hdu().data, reference) + + +@pytest.mark.parametrize('effect_name', ['tree', 'selector']) +def test_missing_branch_id_fails(request, effect_name, source): + effect = request.getfixturevalue(effect_name) + fov = _centre_fov(n=10) + fov.meta.clear() + with pytest.raises(KeyError, match='image_plane_id'): + effect.apply_to(source, fov=fov) + + +def test_observe_keeps_downstream_emission_after_branch_split(tree, source): + common = TERCurve(array_dict={'wavelength': [0.5, 2.5], 'transmission': [0.8, 0.8], 'emissivity': [0, 0]}, wavelength_unit='um') + downstream = TERCurve(array_dict={'wavelength': [0.5, 2.5], 'transmission': [0.5, 0.5], 'emission': [2e-8, 2e-8]}, wavelength_unit='um', emission_unit='photlam') + train = OpticalTrain() + train.cmds = UserCommands() + fovs = [_centre_fov(n=10) for arm in range(3)] + for arm, fov in enumerate(fovs): + fov.meta['image_plane_id'] = arm + train.fov_manager = SimpleNamespace(fovs=fovs, is_spectroscope=False) + train.optics_manager = SimpleNamespace(source_effects=[common, tree, downstream], fov_effects=[], image_plane_effects=[]) + train.image_planes = [ImagePlane(fov.header) for fov in fovs] + + prepared = train.prepare_source(source.make_copy()) + reference = extract_source(_centre_fov(n=10), prepared).make_hdu().data + emission = extract_source(_centre_fov(n=10), downstream.background_source[0]).make_hdu().data + train.observe(source, update=False) + for arm, transmission in enumerate([0.25, 0.3, 0.45]): + expected = reference * 0.8 * transmission * 0.5 + emission + np.testing.assert_allclose(train._last_fovs[arm].hdu.data, expected) + assert train._last_fovs[arm].image_plane_id == arm + + +@pytest.mark.parametrize('spectrum_id, expected', [(0, 0.5), (1, 0.5 + 1 / 24), (2, 0.5 - 1 / 24)]) +def test_wavelength_dependent_branch_preserves_source_colour(spectrum_id, expected): + spectrum = _table_source().fields[0].spectra[spectrum_id] + source = Source(x=[0], y=[0], ref=[0], weight=[1], spectra=[spectrum]) + curve = {'wavelength': [0.5, 2.5], 'transmission': [0, 1], 'emissivity': [0, 0]} + selector = SelectorWheel(selector_key='image_plane_id', wheel=[{'selector_value': 0, 'effect_class': 'QuantumEfficiencyCurve', 'effect_kwargs': {'array_dict': curve, 'wavelength_unit': 'um'}}]) + reference = extract_source(_centre_fov(n=10), source).make_hdu().data.sum() + fov = _centre_fov(n=10) + fov.meta['image_plane_id'] = 0 + measured = extract_source(fov, source, [selector]).make_hdu().data.sum() / reference + # Linear spectra and throughput give these analytic integrals over 1--2 um. + # The input spectra are sampled every 0.02 um. + np.testing.assert_allclose(measured, expected, rtol=0, atol=1e-4) + + +@pytest.mark.parametrize('z_order, expected', [([117, 217, 517], (117, 217, 517)), (517, (517,))]) +def test_dichroic_source_stage_is_opt_in(tree, z_order, expected): + assert tree.z_order == (117, 217, 617) + configured = DichroicTree(z_order=z_order, **{key: tree.meta[key] for key in tree.required_keys}) + assert configured.z_order == expected diff --git a/scopesim/tests/tests_optics/test_adc_slit_illumination.py b/scopesim/tests/tests_optics/test_adc_slit_illumination.py index 4dae1e4d3..3df011292 100644 --- a/scopesim/tests/tests_optics/test_adc_slit_illumination.py +++ b/scopesim/tests/tests_optics/test_adc_slit_illumination.py @@ -157,7 +157,7 @@ def test_observe_applies_shift_before_slit_through_wavelength_throughput(interve effects.append(Displacement(0., displacement, method='roll')) train = OpticalTrain() train.cmds = fov.cmds - train.fov_manager = SimpleNamespace(fovs=[fov]) + train.fov_manager = SimpleNamespace(fovs=[fov], is_spectroscope=True) train.optics_manager = SimpleNamespace(source_effects=[], fov_effects=effects, image_plane_effects=[]) train.image_planes = [SimpleNamespace(add=lambda hdu, **kwargs: images.append(hdu.data.copy()))] train.observe(source, update=False)