From 4e2b76dde36347fe4da22848cd636fd0d4edabb6 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:42:47 -0700 Subject: [PATCH 1/2] 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 fd14b0906d5dc0e3693221cfa1a3f5cd7a8afc13 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:55:53 -0700 Subject: [PATCH 2/2] 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)