diff --git a/scopesim/effects/atmo_dispersion.py b/scopesim/effects/atmo_dispersion.py index 268cd7fe..3fd85fd1 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/electronic/__init__.py b/scopesim/effects/electronic/__init__.py index e086896c..991d1aa3 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 f1123737..328eefb9 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 8c6aef4f..e9d9198a 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 65b5b818..5179aa7d 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/psfs/analytical.py b/scopesim/effects/psfs/analytical.py index dc39f520..1e8bb0ff 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 64d852ef..8b5c3335 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/selector_wheel.py b/scopesim/effects/selector_wheel.py index 7ee03f24..c5869257 100644 --- a/scopesim/effects/selector_wheel.py +++ b/scopesim/effects/selector_wheel.py @@ -13,12 +13,13 @@ 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 from ..optics.image_plane import ImagePlane from ..detector.detector import Detector +from ..source.source import Source logger = get_logger(__name__) @@ -62,8 +63,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") @@ -78,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/spectral_efficiency.py b/scopesim/effects/spectral_efficiency.py index 2baae5bd..b551b72f 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 3e7fbc71..09e0af97 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 c0d68eb2..f4a6e421 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 @@ -232,7 +233,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 @@ -1556,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) @@ -1563,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 = {} @@ -1579,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. @@ -1617,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/fov.py b/scopesim/optics/fov.py index a2b2ffa9..3f56c4e2 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 4e7a894f..261f48f8 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 3ecba4d4..f65c86f7 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 15f91daf..d64c995a 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 a6031f2a..ebad53f7 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__ @@ -35,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): @@ -532,9 +541,26 @@ 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 + # 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 @@ -546,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, @@ -562,14 +588,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)) + extract_source(fov, source, branch_effects, spatial_header=spatial_header, _spectral_wavesets=spectral_wavesets) + else: + extract_source(fov, source, branch_effects, _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 +655,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 +868,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 a578674a..a72cc531 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 27706b5d..cfbbd337 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 1119df8d..03b7bbcc 100644 --- a/scopesim/tests/tests_effects/test_MoffatPSF.py +++ b/scopesim/tests/tests_effects/test_MoffatPSF.py @@ -1,14 +1,43 @@ import logging +from unittest.mock import Mock 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 +45,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 +120,165 @@ 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. + + +@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_calibrated_spectral_trace.py b/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py index 00114bce..e263f9ca 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_detector_operating_points.py b/scopesim/tests/tests_effects/test_detector_operating_points.py new file mode 100644 index 00000000..1763f1d2 --- /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 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 00000000..4f74407b --- /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_imaging_branches.py b/scopesim/tests/tests_effects/test_imaging_branches.py new file mode 100644 index 00000000..ce308f39 --- /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_effects/test_shift_planes.py b/scopesim/tests/tests_effects/test_shift_planes.py new file mode 100644 index 00000000..5f4bcabe --- /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 00000000..e8ef5a18 --- /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 61ced191..bab799ad 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 c941e903..2866f345 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 1c6b77ca..cd055f5d 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 00000000..3df01129 --- /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], 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) + 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 00000000..2a12e4ba --- /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.)