From 4296266bc2aadd12cee2ad4ba23940d23585d212 Mon Sep 17 00:00:00 2001 From: Jeb Bailey Date: Fri, 23 Oct 2026 12:39:12 -0700 Subject: [PATCH] Use calibrated traces and physical slit and detector coordinates Source commits: fe052a1d, c76a4c21, 93237d85, c0e0dfdc, 7ffffb87, fa2598b7 (command arrays), f9dccd8f (spectral padding), 3264a6e6, ce5c8cd5, 06f64b09 --- scopesim/effects/apertures.py | 4 + scopesim/effects/illumination.py | 11 +- scopesim/effects/metis_lms_trace_list.py | 5 +- scopesim/effects/spectral_trace_list_utils.py | 180 +++++++++++------- scopesim/effects/ter_curves.py | 7 +- scopesim/tests/test_utils_functions.py | 9 + .../tests_effects/test_SkycalcTERCurve.py | 4 +- .../tests_effects/test_SpectralTraceList.py | 101 +++++++++- .../test_SpectralTraceListUtils.py | 8 +- .../test_calibrated_spectral_trace.py | 87 +++++++++ .../tests/tests_effects/test_local_skycalc.py | 24 +++ .../test_registered_background.py | 32 ++++ scopesim/utils.py | 5 +- 13 files changed, 392 insertions(+), 85 deletions(-) create mode 100644 scopesim/tests/tests_effects/test_calibrated_spectral_trace.py create mode 100644 scopesim/tests/tests_effects/test_local_skycalc.py create mode 100644 scopesim/tests/tests_effects/test_registered_background.py diff --git a/scopesim/effects/apertures.py b/scopesim/effects/apertures.py index aca858141..199914541 100644 --- a/scopesim/effects/apertures.py +++ b/scopesim/effects/apertures.py @@ -128,6 +128,8 @@ def apply_to(self, obj, **kwargs): for vol in obj.volumes: vol["meta"]["xi_min"] = min(x) * u.arcsec vol["meta"]["xi_max"] = max(x) * u.arcsec + vol["meta"]["eta_min"] = min(y) * u.arcsec + vol["meta"]["eta_max"] = max(y) * u.arcsec vol["meta"]["slit_name"] = self.meta["name"] if "offset_type" in self.meta: @@ -305,6 +307,8 @@ def apply_to(self, obj, **kwargs): # ..todo: HUGE HACK - Get rid of this! vol["meta"]["xi_min"] = row["left"] * u.arcsec vol["meta"]["xi_max"] = row["right"] * u.arcsec + vol["meta"]["eta_min"] = row["bottom"] * u.arcsec + vol["meta"]["eta_max"] = row["top"] * u.arcsec vol["conserve_image"] = row["conserve_image"] vol["shape"] = row["shape"] vol["angle"] = row["angle"] diff --git a/scopesim/effects/illumination.py b/scopesim/effects/illumination.py index d4ff12bfe..65b5b8181 100644 --- a/scopesim/effects/illumination.py +++ b/scopesim/effects/illumination.py @@ -538,6 +538,7 @@ class ImagePlaneBackground(Effect): trace list. The supplied value or model should already include downstream throughput and detector QE. For tapered QE maps, use an average or otherwise representative positional QE when deriving this non-dispersed background. + A FITS ``filename`` and ``extname`` can supply an already registered map. """ z_order: ClassVar[tuple[int, ...]] = (760,) @@ -562,7 +563,7 @@ def apply_to(self, obj, **kwargs): return obj shape = obj.hdu.data.shape - if self._map is None or shape != self._map_shape: + if self._file is not None or self._map is None or shape != self._map_shape: self._map = self._make_map(shape) self._map_shape = shape @@ -570,7 +571,13 @@ def apply_to(self, obj, **kwargs): return obj def _make_map(self, shape): - if self._model is None: + if self._file is not None: + # The calibration image must be registered to this image plane. + extname = str(from_currsys(self.meta["extname"], self.cmds)) + background = self._file[extname].data + if background.shape != shape: + raise ValueError("Background calibration shape differs from the image plane") + elif self._model is None: value = from_currsys(self.meta["value"], self.cmds) background = np.full(shape, value, dtype=np.float32) else: diff --git a/scopesim/effects/metis_lms_trace_list.py b/scopesim/effects/metis_lms_trace_list.py index 92780daca..64f6f460d 100644 --- a/scopesim/effects/metis_lms_trace_list.py +++ b/scopesim/effects/metis_lms_trace_list.py @@ -144,6 +144,8 @@ def apply_to_fov(self, obj): slicefov.detector_header = obj.detector_header slicefov.meta["xi_min"] = obj.meta["xi_min"] slicefov.meta["xi_max"] = obj.meta["xi_max"] + slicefov.meta["eta_min"] = ymin * u.arcsec + slicefov.meta["eta_max"] = ymax * u.arcsec slicefov.meta["trace_id"] = sptid slicefov.cube = fits.ImageHDU(header=slicewcs.to_header(), data=slicecube) @@ -322,6 +324,7 @@ class MetisLMSSpectralTrace(SpectralTrace): """SpectralTrace for the METIS LM spectrograph.""" _class_params = { + **SpectralTrace._class_params, "naxis1": 122, "nslice": 28, "slicewidth": 0.0207, # arcsec @@ -393,7 +396,7 @@ def get_waverange(self, det_mm_lims): lam0 = from_currsys(self.meta["wavelen"], self.cmds) xi0 = 0. - ymid = self.xilam2y(xi0, lam0)[0] # estimate y level of trace + ymid = self.xilam2y(xi0, lam0) # scalar input gives a scalar trace coordinate waverange = self.xy2lam(np.array([xmin, xmax]), np.array([ymid, ymid]), grid=False) diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py index d5635bc71..3e7fbc71b 100644 --- a/scopesim/effects/spectral_trace_list_utils.py +++ b/scopesim/effects/spectral_trace_list_utils.py @@ -14,8 +14,8 @@ import numpy as np -from scipy.interpolate import RectBivariateSpline -from scipy.interpolate import interp1d +from scipy.interpolate import interp1d, RectBivariateSpline +from scipy.ndimage import gaussian_filter1d, map_coordinates from astropy.table import Table, vstack from astropy.modeling import fitting @@ -103,6 +103,7 @@ class SpectralTrace: "pixel_size": None, "description": "", "trace_flux_jacobian": "!SIM.spectral.trace_flux_jacobian", + "fit_degree": 4, } def __init__(self, trace_tbl, cmds=None, **kwargs): @@ -128,6 +129,9 @@ def __init__(self, trace_tbl, cmds=None, **kwargs): ("PLTSCALE", "plate_scale"), ("DETPAD", "detector_pad"), ("DETANG", "detector_angle"), + ("FITDEG", "fit_degree"), + ("LSFRES", "optical_resolving_power"), + ("SLITCOL", "slit_dispersion_colname"), ): if header_key in trace_tbl.header: self.meta[meta_key] = trace_tbl.header[header_key] @@ -194,12 +198,19 @@ def compute_interpolation_functions(self): self.wave_min = quantify(np.min(lam_arr), u.um).value self.wave_max = quantify(np.max(lam_arr), u.um).value - self.xy2xi = Transform2D.fit(x_arr, y_arr, xi_arr) - self.xy2lam = Transform2D.fit(x_arr, y_arr, lam_arr) - self.xilam2x = Transform2D.fit(xi_arr, lam_arr, x_arr) - self.xilam2y = Transform2D.fit(xi_arr, lam_arr, y_arr) - self._xiy2x = Transform2D.fit(xi_arr, y_arr, x_arr) - self._xiy2lam = Transform2D.fit(xi_arr, y_arr, lam_arr) + degree = self.meta["fit_degree"] + self.xy2xi = Transform2D.fit(x_arr, y_arr, xi_arr, degree=degree) + self.xy2lam = Transform2D.fit(x_arr, y_arr, lam_arr, degree=degree) + self.xilam2x = Transform2D.fit(xi_arr, lam_arr, x_arr, degree=degree) + self.xilam2y = Transform2D.fit(xi_arr, lam_arr, y_arr, degree=degree) + self._xiy2x = Transform2D.fit(xi_arr, y_arr, x_arr, degree=degree) + self._xiy2lam = Transform2D.fit(xi_arr, y_arr, lam_arr, degree=degree) + if "slit_dispersion_colname" in self.meta: + center = xi_arr == 0 + derivative = self.table[self.meta["slit_dispersion_colname"]].quantity.to_value(u.um / u.arcsec) + # FOVs are restricted to this trace; allow unit-conversion roundoff + # at the wavelength endpoints when interpolating the derivative. + self.slit_dispersion = interp1d(lam_arr[center], derivative[center], bounds_error=False, fill_value="extrapolate") if self.dispersion_axis == "unknown": dlam_dx, dlam_dy = self.xy2lam.gradient() @@ -261,6 +272,24 @@ def _focal_plane_grid(self, fov): y_mm = np.linspace(ymin_mm, ymax_mm, sub_naxis2, endpoint=False, dtype=np.float32) return xmin, ymin, x_mm, y_mm, det_wcs + def spectral_padding(self, fov): + """Input wavelength support for slit displacement and optical blur.""" + lower, upper = fov.meta['wave_min'].to_value(u.um), fov.meta['wave_max'].to_value(u.um) + wave = np.unique(u.Quantity(self.table[self.meta['wave_colname']]).to_value(u.um)) + wave = np.r_[lower, wave[(wave > lower) & (wave < upper)], upper] + width = (fov.meta['eta_max'] - fov.meta['eta_min']).to_value(u.arcsec) + if 'slit_dispersion_colname' in self.meta: + derivative = self.slit_dispersion(wave) + else: + pixsize = (fov.detector_header['CDELT1D']*u.Unit(fov.detector_header['CUNIT1D'])).to_value(u.mm) + self._set_dispersion(lower, upper, pixsize=pixsize) + step = abs(fov.header['CDELT2'])*3600 + derivative = self.dlam_per_pix(wave)/step + padding = np.max(abs(derivative))*width/2 + if 'optical_resolving_power' in self.meta: + padding += 5*upper/(self.meta['optical_resolving_power']*np.sqrt(8*np.log(2))) + return padding*u.um + def map_spectra_to_focal_plane(self, fov): """ Apply the spectral trace mapping to a spectral cube. @@ -294,14 +323,13 @@ def map_spectra_to_focal_plane(self, fov): image = np.zeros((sub_naxis2, sub_naxis1), dtype=np.float32) self._set_dispersion(wave_min, wave_max, pixsize=pixsize) - try: - xilam = XiLamImage(fov, self.dlam_per_pix) - self._xilamimg = xilam # ..todo: remove or make available with a debug flag? - except ValueError: - logger.warning(" ---> %s gave ValueError", self.trace_id) - - npix_xi, npix_lam = xilam.npix_xi, xilam.npix_lam - xilam_wcs = xilam.wcs + slit_dispersion = None + if "slit_dispersion_colname" in self.meta: + slit_dispersion = self.slit_dispersion + xilam = XiLamImage(fov, self.dlam_per_pix, slit_dispersion=slit_dispersion) + if "optical_resolving_power" in self.meta: + xilam.convolve_lsf(self.meta["optical_resolving_power"]) + self._xilamimg = xilam # focal-plane coordinate images ximg_fpa, yimg_fpa = np.meshgrid(x_mm, y_mm) @@ -310,33 +338,11 @@ def map_spectra_to_focal_plane(self, fov): xi_fpa = self.xy2xi(ximg_fpa, yimg_fpa).astype(np.float32) lam_fpa = self.xy2lam(ximg_fpa, yimg_fpa).astype(np.float32) - # mask everything outside the wavelength range - mask = (xi_fpa >= xi_min) & (xi_fpa <= xi_max) - xi_fpa *= mask - lam_fpa *= mask - - # mask everything outside the wavelength range - # mask = (lam_fpa >= wave_min) & (lam_fpa <= wave_max) - # xi_fpa *= mask - # lam_fpa *= mask - - # Convert to pixel images - # These are the pixel coordinates in the image corresponding to - # xi, lambda - # It is much quicker to do the linear transformation by hand - # than to use the astropy.wcs functions for conversion. - i_img = ((lam_fpa - xilam_wcs.wcs.crval[0]) - / xilam_wcs.wcs.cdelt[0]).astype(int) - j_img = ((xi_fpa - xilam_wcs.wcs.crval[1]) - / xilam_wcs.wcs.cdelt[1]).astype(int) - - # truncate images to remove pixel coordinates outside the image - ijmask = ((i_img >= 0) * (i_img < npix_lam) - * (j_img >= 0) * (j_img < npix_xi)) - - # do the actual interpolation - # image is in [ph/s/um/arcsec] - image = xilam.interp(xi_fpa, lam_fpa, grid=False) * ijmask + # The physical slit ends include the outer half-cells. A truncated + # integer index incorrectly admits pixels below the lower boundary. + inside = (xi_fpa >= xi_min) & (xi_fpa < xi_max) + inside &= (lam_fpa >= wave_min) & (lam_fpa < wave_max) + image = xilam.interp(xi_fpa, lam_fpa, grid=False)*inside # Scale to ph / s / pixel image *= self._trace_flux_scale( @@ -380,11 +386,7 @@ def map_spectra_to_focal_plane(self, fov): def _trace_flux_scale(self, x_mm, y_mm, pixsize, pixscale, det_header): """Return the local [arcsec um] per detector pixel flux scale.""" - try: - use_jacobian = bool(from_currsys( - self.meta["trace_flux_jacobian"], self.cmds)) - except ValueError: - use_jacobian = False + use_jacobian = from_currsys(self.meta["trace_flux_jacobian"], self.cmds) if not use_jacobian: dlam_by_dx, dlam_by_dy = self.xy2lam.gradient() @@ -401,13 +403,10 @@ def _trace_flux_scale(self, x_mm, y_mm, pixsize, pixscale, det_header): @staticmethod def _detector_pixel_area(det_header): """Return detector WCS pixel area in mm2.""" - try: - pixel_scale_matrix = WCS(det_header, key="D").pixel_scale_matrix - return np.abs(np.linalg.det(pixel_scale_matrix)) - except Exception: # pragma: no cover - FITS fallback for odd headers - cdelt1 = det_header["CDELT1D"] * u.Unit(det_header["CUNIT1D"]) - cdelt2 = det_header["CDELT2D"] * u.Unit(det_header["CUNIT2D"]) - return np.abs((cdelt1 * cdelt2).to_value(u.mm**2)) + detector_wcs = WCS(det_header, key="D") + pixel_area = np.abs(np.linalg.det(detector_wcs.pixel_scale_matrix)) + pixel_area *= u.Unit(detector_wcs.wcs.cunit[0]) * u.Unit(detector_wcs.wcs.cunit[1]) + return pixel_area.to_value(u.mm**2) def rectify( self, @@ -851,9 +850,13 @@ class XiLamImage(): dlam_per_pix : a 1-D interpolation function from wavelength (in um) to dispersion (in um/pixel); alternatively a number giving an average dispersion + slit_dispersion : callable, optional + Physical cross-slit derivative d(lambda)/d(eta) in um/arcsec as a + function of wavelength in um. When supplied, this replaces the + legacy coupling between sky-grid sampling and detector dispersion. """ - def __init__(self, fov, dlam_per_pix): + def __init__(self, fov, dlam_per_pix, slit_dispersion=None): # ..todo: we assume that we always have a cube. We use SpecCADO's # add_cube_layer method cube_wcs = WCS(fov.cube.header, key=" ") @@ -871,8 +874,16 @@ def __init__(self, fov, dlam_per_pix): (n_lam, n_eta, n_xi) = fov.cube.data.shape # arrays of cube coordinates - cube_xi = d_xi * np.arange(n_xi) + fov.meta["xi_min"].value - cube_eta = d_eta * (np.arange(n_eta) - (n_eta - 1) / 2) + cube_xi = cube_wcs.all_pix2world(np.arange(n_xi), np.zeros(n_xi), np.zeros(n_xi), 0)[0] + cube_eta = cube_wcs.all_pix2world(np.zeros(n_eta), np.arange(n_eta), np.zeros(n_eta), 0)[1] + if cube_wcs.has_celestial: + # Slit coordinates are offsets about zero; RA wraps at 360 degrees. + cube_xi = (cube_xi + 180) % 360 - 180 + cube_xi *= u.Unit(cube_wcs.wcs.cunit[0]).to(u.arcsec) + cube_eta *= u.Unit(cube_wcs.wcs.cunit[1]).to(u.arcsec) + # Each IFU slice has its own optical axis, independent of sky position. + eta_center = (fov.meta["eta_min"] + fov.meta["eta_max"]) / 2 + cube_eta -= eta_center.to_value(u.arcsec) cube_lam = wcs_lam.all_pix2world(np.arange(n_lam), 0)[0] cube_lam *= u.Unit(wcs_lam.wcs.cunit[0]).to(u.um) @@ -886,24 +897,31 @@ def __init__(self, fov, dlam_per_pix): logger.warning("Using scalar dlam_per_pix = %.2g", dlam_per_pix_val) + dlam_per_arcsec = dlam_per_pix_val / d_eta if slit_dispersion is None else slit_dispersion(self.lam) for i, eta in enumerate(cube_eta): - lam0 = self.lam + dlam_per_pix_val * eta / d_eta + # A positive calibrated derivative moves an input line to a larger + # apparent wavelength: sample the source at lambda_out - shift. + lam0 = self.lam + dlam_per_arcsec * eta if slit_dispersion is None else self.lam - dlam_per_arcsec * eta # lam0 is the target wavelength. We need to check that this # overlaps with the wavelength range covered by the cube if lam0.min() < cube_lam.max() and lam0.max() > cube_lam.min(): - plane = fov.cube.data[:, i, :].T - plane_interp = RectBivariateSpline(cube_xi, cube_lam, plane, - kx=1, ky=1) - self.image += plane_interp(cube_xi, lam0) + # Only wavelength changes here; xi is already on its native + # grid. Building a 2D spline for every slit cell is redundant. + index = np.searchsorted(cube_lam, lam0, side='right') - 1 + index = np.clip(index, 0, n_lam - 2) + fraction = (lam0 - cube_lam[index])/d_lam + valid = (lam0 >= cube_lam[0]) & (lam0 <= cube_lam[-1]) + plane = fov.cube.data[:, i, :] + self.image += ((1 - fraction[:, None])*plane[index] + fraction[:, None]*plane[index + 1]).T*valid - self.image *= d_eta # ph/s/um/arcsec2 --> ph/s/um/arcsec + self.image *= abs(d_eta) # ph/s/um/arcsec2 --> ph/s/um/arcsec # WCS for the xi-lambda image, i.e. the rectified 2D spectrum # Default WCS with xi in arcsec self.wcs = WCS(naxis=2) self.wcs.wcs.crpix = [1, 1] - self.wcs.wcs.crval = [self.lam[0], fov.meta["xi_min"].value] + self.wcs.wcs.crval = [self.lam[0], cube_xi[0]] self.wcs.wcs.pc = [[1, 0], [0, 1]] self.wcs.wcs.cdelt = [d_lam, d_xi] self.wcs.wcs.ctype = ["LINEAR", "LINEAR"] @@ -923,14 +941,30 @@ def __init__(self, fov, dlam_per_pix): self.xi = self.wcs.all_pix2world(self.lam[0], np.arange(n_xi), 0)[1] self.npix_xi = n_xi self.npix_lam = n_lam - # ..todo: cubic spline introduces negative values, linear does not. - # Alternative might be to cubic-spline interpolate on sqrt(image), - # with subsequent squaring of the result. This would require - # wrapping RectBivariateSpline in a new (sub)class. - spline_order = (1, 1) - self.interp = RectBivariateSpline(self.xi, self.lam, self.image, - kx=spline_order[0], - ky=spline_order[1]) + + def interp(self, xi, wavelength, grid=True): + """Bilinear density interpolation on the existing uniform WCS grid.""" + xi, wavelength = np.meshgrid(xi, wavelength, indexing='ij') if grid else np.broadcast_arrays(xi, wavelength) + row = (xi - self.xi[0])/self.wcs.wcs.cdelt[1] + column = (wavelength - self.lam[0])/self.wcs.wcs.cdelt[0] + return map_coordinates(self.image, [row, column], order=1, mode='nearest', prefilter=False) + + def convolve_lsf(self, resolving_power): + """Apply an optical Gaussian LSF at constant R, preserving photon density. + + This is the optical blur alone; slit illumination is integrated above. + Convolving lambda * I_lambda on a uniform log(lambda) grid preserves + integrated flux rather than treating I_lambda as a log-bin count. + """ + if resolving_power <= 0: + raise ValueError("Optical resolving power must be positive") + log_wave = np.log(self.lam) + n = int(np.ceil((log_wave[-1] - log_wave[0]) / np.min(np.diff(log_wave)))) + 1 + log_grid = np.linspace(log_wave[0], log_wave[-1], n) + sigma = 1 / (resolving_power * np.sqrt(8 * np.log(2)) * (log_grid[1] - log_grid[0])) + density = interp1d(log_wave, self.image * self.lam, axis=1)(log_grid) + blurred = gaussian_filter1d(density, sigma, axis=1, mode="nearest") + self.image = interp1d(log_grid, blurred, axis=1)(log_wave) / self.lam class Transform2D(): diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py index 6efd90f7b..c0d68eb2e 100644 --- a/scopesim/effects/ter_curves.py +++ b/scopesim/effects/ter_curves.py @@ -547,9 +547,14 @@ def load_skycalc_table(self): else: path = find_file(use_local_file) - fits_tbl = fits.getdata(path, ext=1) + with fits.open(path) as sky_file: + fits_tbl = sky_file[1].data + # Older cached files require an explicit unit if TUNIT is absent. + unit = self.meta["local_skycalc_wavelength_unit"] if "local_skycalc_wavelength_unit" in self.meta else sky_file[1].columns["lam"].unit + wave_unit = u.Unit(unit) fits_hdr = fits.getheader(path, ext=0) tbl = Table(fits_tbl) + tbl["lam"] = tbl["lam"] * wave_unit.to(u.um) tbl["lam"].unit = u.um for colname in tbl.colnames: if "flux" in colname: diff --git a/scopesim/tests/test_utils_functions.py b/scopesim/tests/test_utils_functions.py index a1df0f1b3..1469b32d9 100644 --- a/scopesim/tests/test_utils_functions.py +++ b/scopesim/tests/test_utils_functions.py @@ -199,6 +199,15 @@ def test_converts_list(self): def test_converts_numpy_array(self): assert utils.from_currsys(np.array(["!SIM.random.seed"]*2))[1] is None + @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32]) + def test_numeric_calibration_array_keeps_values_and_precision(self, dtype): + values = np.arange(12, dtype=dtype).reshape(3, 4) + resolved = utils.from_currsys(values) + np.testing.assert_array_equal(resolved, values) + assert resolved.dtype == dtype + resolved[0, 0] = 20 + assert values[0, 0] == 0 + def test_converts_dict(self): assert utils.from_currsys({"seed": "!SIM.random.seed"})["seed"] is None diff --git a/scopesim/tests/tests_effects/test_SkycalcTERCurve.py b/scopesim/tests/tests_effects/test_SkycalcTERCurve.py index 689621b10..b87995c81 100644 --- a/scopesim/tests/tests_effects/test_SkycalcTERCurve.py +++ b/scopesim/tests/tests_effects/test_SkycalcTERCurve.py @@ -59,7 +59,9 @@ def test_initialises_with_non_skycalc_keys(self): def test_initialise_with_local_skycalc_file(self, mock_path): sky_ter = SkycalcTERCurve( - use_local_skycalc_file=str(mock_path / "skycalc_override.fits")) + use_local_skycalc_file=str(mock_path / "skycalc_override.fits"), + local_skycalc_wavelength_unit=u.um) + # This legacy fixture predates wavelength units in its FITS columns. assert sky_ter.skycalc_table is not None diff --git a/scopesim/tests/tests_effects/test_SpectralTraceList.py b/scopesim/tests/tests_effects/test_SpectralTraceList.py index 8461857a2..5ff0633aa 100644 --- a/scopesim/tests/tests_effects/test_SpectralTraceList.py +++ b/scopesim/tests/tests_effects/test_SpectralTraceList.py @@ -6,6 +6,7 @@ import numpy as np from astropy import units as u from astropy.io import fits +from astropy.table import Table from scopesim.effects.spectral_trace_list import SpectralTraceList, \ @@ -137,6 +138,93 @@ def _echelle_cmds(): class TestEchelleSpectralTraceList: + @pytest.mark.parametrize('angle', [-25, 0, 25]) + @pytest.mark.parametrize('offset', [(0., 0.), (1.5, -1.2)]) + def test_fixed_placement_and_detector_intersections(self, tmp_path, angle, offset): + params_file = tmp_path / 'echelle_trace_parameters.dat' + _write_echelle_trace_params(params_file, detector_angle=angle) + params = Table.read(params_file, format='ascii') + for key, value in {'m0': 75, 'n': 6, 'min_wave': 960, 'max_wave': 1340, 'n_disp': 2048, 'n_xdisp': 1024, 'disp_freq': 24., 'xdisp_freq': 350.}.items(): + params[key][0] = value + params['xbeta_center'] = [10.9] + params['detector_offset_x'], params['detector_offset_y'] = [offset[0]], [offset[1]] + params.meta['comments'].extend(['detector_offset_x_unit : mm', 'detector_offset_y_unit : mm']) + params.write(params_file, format='ascii', overwrite=True) + effect = EchelleSpectralTraceList(cmds=_echelle_cmds(), filename=str(params_file), wave_colname='wavelength', s_colname='s') + model = effect._spectrographs['b'] + nx, ny = model.detector.n_pix_x, model.detector.n_pix_y + raw_grid = np.linspace(-2*nx, 3*nx, 50001) + for order, edges in zip(model.orders, model.edge_wave(fsr=True)): + wave = np.linspace(edges[0], edges[1], 4*nx + 1) + x, y = model._detector_pixels(wave, order) + # Rotate about the camera origin, then translate in detector axes. + xraw = model.wavelength_to_x_pixel(wave, order).value - nx/2 + yraw = model.wavelength_to_y_pixel(wave).value[None, :] - ny/2 + model._slit_offset_pix[:, None] + c, s = np.cos(np.deg2rad(angle)), np.sin(np.deg2rad(angle)) + expected_x = c*xraw - s*yraw + nx/2 + offset[0]/.015 + expected_y = s*xraw + c*yraw + ny/2 + offset[1]/.015 + np.testing.assert_allclose(x, expected_x, atol=1e-9, rtol=0) + np.testing.assert_allclose(y, expected_y, atol=1e-9, rtol=0) + # Translation changes placement, not dispersion or slit length. + xrot, yrot = model._rotated_detector_pixels(wave, order) + np.testing.assert_allclose(np.diff(x), np.diff(xrot), atol=1e-9, rtol=0) + np.testing.assert_allclose(np.diff(y), np.diff(yrot), atol=1e-9, rtol=0) + slit_length = np.hypot(x[2] - x[0], y[2] - y[0]) + np.testing.assert_allclose(slit_length, model._slit_offset_pix[2] - model._slit_offset_pix[0], atol=1e-9, rtol=0) + + # A deliberately oversized wavelength search independently checks + # that the translated clipping search has not lost either edge. + wave = model.x_pixel_to_wavelength(raw_grid, order) + x, y = model._detector_pixels(wave, order) + inside = np.any((x >= 0) & (x <= nx) & (y >= 0) & (y <= ny), axis=0) + trace_id = f'b_{order}' + if not inside.any(): + assert trace_id not in effect.spectral_traces + continue + trace = effect.spectral_traces[trace_id] + actual = model.wavelength_to_x_pixel(u.Quantity(trace.table['wavelength']), order).value + expected = raw_grid[inside][[0, -1]] + # The generator locates each edge on a roughly one-pixel grid. + np.testing.assert_allclose([actual.min(), actual.max()], expected, atol=3, rtol=0) + + @pytest.mark.parametrize('missing', ['detector_offset_x', 'detector_offset_y', 'detector_offset_x_unit', 'detector_offset_y_unit']) + def test_incomplete_placement_fails(self, tmp_path, missing): + path = tmp_path / 'echelle_trace_parameters.dat' + _write_echelle_trace_params(path) + params = Table.read(path, format='ascii') + for axis in ['x', 'y']: + key = f'detector_offset_{axis}' + if missing != key: + params[key] = [0.] + if missing != f'{key}_unit': + params.meta['comments'].append(f'{key}_unit : mm') + params.write(path, format='ascii', overwrite=True) + with pytest.raises(KeyError, match=missing): + EchelleSpectralTraceList(cmds=_echelle_cmds(), filename=str(path)) + + @pytest.mark.parametrize('missing_value', ['detector_offset_x', 'detector_offset_y']) + def test_blank_placement_value_does_not_become_zero(self, tmp_path, missing_value): + path = tmp_path / 'echelle_trace_parameters.dat' + _write_echelle_trace_params(path) + params = Table.read(path, format='ascii') + for axis in ['x', 'y']: + key = f'detector_offset_{axis}' + params[key] = np.ma.array([0.], mask=[key == missing_value]) + params.meta['comments'].append(f'{key}_unit : mm') + params.write(path, format='ascii', overwrite=True) + with pytest.raises(ValueError, match='Both detector offsets must have values'): + EchelleSpectralTraceList(cmds=_echelle_cmds(), filename=str(path)) + + def test_placement_requires_length_units(self, tmp_path): + path = tmp_path / 'echelle_trace_parameters.dat' + _write_echelle_trace_params(path) + params = Table.read(path, format='ascii') + params['detector_offset_x'], params['detector_offset_y'] = [0.], [0.] + params.meta['comments'].extend(['detector_offset_x_unit : arcsec', 'detector_offset_y_unit : mm']) + params.write(path, format='ascii', overwrite=True) + with pytest.raises(u.UnitConversionError): + EchelleSpectralTraceList(cmds=_echelle_cmds(), filename=str(path)) + def test_does_not_write_generated_hdulist_by_default(self, tmp_path): params_file = tmp_path / "echelle_trace_parameters.dat" _write_echelle_trace_params(params_file) @@ -197,6 +285,13 @@ def test_detector_angle_clips_to_detector_not_padding(self, tmp_path): assert trace.meta["detector_angle"] == 25 assert trace.meta["detector_pad"] == 10 + model = spt._spectrographs['b'] + order = int(trace.meta['trace_id'].split('_')[1]) + wave = u.Quantity(trace.table['wavelength'][:n_wave]) + model_x, model_y = model._detector_pixels(wave, order) + np.testing.assert_allclose(model_x, xpix, atol=1e-8, rtol=0) + np.testing.assert_allclose(model_y, ypix, atol=1e-8, rtol=0) + padded_file = tmp_path / "echelle_trace_parameters_more_padding.dat" _write_echelle_trace_params( padded_file, @@ -224,7 +319,9 @@ def test_detector_angle_clips_to_detector_not_padding(self, tmp_path): ) unrotated_trace = next(iter(unrotated.spectral_traces.values())) - assert trace.wave_min == pytest.approx(unrotated_trace.wave_min) - assert trace.wave_max == pytest.approx(unrotated_trace.wave_max) + # Rotation changes which wavelengths reach the rectangular detector; + # clipping must occur after rotation, not at the old unrotated edges. + assert trace.wave_min < trace.wave_max + assert not np.allclose([trace.wave_min, trace.wave_max], [unrotated_trace.wave_min, unrotated_trace.wave_max]) assert not np.allclose(trace.table["x"], unrotated_trace.table["x"]) assert not np.allclose(trace.table["y"], unrotated_trace.table["y"]) diff --git a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py index dfaa2d7f6..04be6beef 100644 --- a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py +++ b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py @@ -196,7 +196,7 @@ def test_call_gives_correct_result(self, quadratic): # matrix and explicit function td2d = Transform2D(quadratic['matrix']) - assert td2d(x, y) == quadratic['function'](x, y) + assert td2d(x, y) == pytest.approx(quadratic['function'](x, y)) def test_gradient_gives_correct_result(self, quadratic): x = np.random.randn() @@ -205,8 +205,8 @@ def test_gradient_gives_correct_result(self, quadratic): tf2d = Transform2D(quadratic['matrix']) tf2d_grad = tf2d.gradient() - assert tf2d_grad[0](x, y) == quadratic['gradient'][0](x, y) - assert tf2d_grad[1](x, y) == quadratic['gradient'][1](x, y) + assert tf2d_grad[0](x, y) == pytest.approx(quadratic['gradient'][0](x, y)) + assert tf2d_grad[1](x, y) == pytest.approx(quadratic['gradient'][1](x, y)) def test_grid_true_gives_correct_shape(self, tf2d): n_x, n_y = 12, 3 @@ -303,7 +303,7 @@ def __init__(self, n_lam=20, n_eta=3, n_xi=11): hdr["CUNIT1"], hdr["CUNIT2"], hdr["CUNIT3"] = "arcsec", "arcsec", "um" self.cube = fits.ImageHDU( data=np.ones((n_lam, n_eta, n_xi)), header=hdr) - self.meta = {"xi_min": -0.5 * u.arcsec, "xi_max": 0.5 * u.arcsec} + self.meta = {"xi_min": -0.5 * u.arcsec, "xi_max": 0.5 * u.arcsec, "eta_min": -0.5 * u.arcsec, "eta_max": 0.5 * u.arcsec} class TestXiLamImage: diff --git a/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py b/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py new file mode 100644 index 000000000..00114bcee --- /dev/null +++ b/scopesim/tests/tests_effects/test_calibrated_spectral_trace.py @@ -0,0 +1,87 @@ +from types import SimpleNamespace +import numpy as np +import pytest +from astropy import units as u +from astropy.io import fits +from scopesim.effects.spectral_trace_list_utils import XiLamImage, SpectralTrace + + +def line_cube(step, offset=0): + wave = np.linspace(.995, 1.005, 2001) + eta = np.arange(-.5+step/2, .5, step) + spatial = np.exp(-.5*((eta-offset)/.12)**2) + profile = np.exp(-.5*((wave-1)/.00002)**2) + cube = profile[:, None, None]*spatial[None, :, None]*np.ones((1, 1, 3)) + header = fits.Header({'CTYPE1': 'LINEAR', 'CTYPE2': 'LINEAR', 'CTYPE3': 'WAVE', 'CUNIT1': 'arcsec', 'CUNIT2': 'arcsec', 'CUNIT3': 'um', 'CDELT1': .1, 'CDELT2': step, 'CDELT3': wave[1]-wave[0], 'CRPIX1': 1., 'CRPIX2': 1., 'CRPIX3': 1., 'CRVAL1': 0., 'CRVAL2': eta[0], 'CRVAL3': wave[0]}) + return SimpleNamespace(cube=fits.ImageHDU(cube, header), meta={'xi_min': 0*u.arcsec, 'eta_min': -.5*u.arcsec, 'eta_max': .5*u.arcsec}) + + + +def test_slit_dispersion_is_independent_of_numerical_sky_grid(): + spectra = [XiLamImage(line_cube(step), lambda w: .00001, slit_dispersion=lambda w: .0002) for step in [.1, .05, .025]] + # Finite spatial quadrature changes, while the physical slit dispersion stays fixed. + for image in spectra[:-1]: + np.testing.assert_allclose(image.image, spectra[-1].image, atol=.001, rtol=.02) + + + +def test_calibrated_slit_shift_has_forward_model_sign(): + image = XiLamImage(line_cube(.025, offset=.1), lambda w: .00001, slit_dispersion=lambda w: .002) + centroid = np.sum(image.image[1]*image.lam)/image.image[1].sum() + assert centroid == pytest.approx(1.0002, abs=2e-6) + + + +def test_optical_lsf_preserves_photons_and_adds_variance(): + image = XiLamImage(line_cube(.05), lambda w: 0., slit_dispersion=lambda w: 0.) + before = np.trapezoid(image.image[1], image.lam) + image.convolve_lsf(10000) + after = np.trapezoid(image.image[1], image.lam) + assert after == pytest.approx(before, rel=.001) + center = np.sum(image.lam*image.image[1])/image.image[1].sum() + variance = np.sum((image.lam-center)**2*image.image[1])/image.image[1].sum() + assert variance == pytest.approx(.00002**2 + (1/10000/np.sqrt(8*np.log(2)))**2, rel=.02) + + +def test_slit_coordinates_follow_cube_wcs(): + fov = line_cube(.025) + fov.cube.header['CRVAL1'] = 2. + fov.cube.header['CRPIX1'] = 2. + fov.cube.header['CRVAL2'] += .1 + image = XiLamImage(fov, lambda wave: .00001, slit_dispersion=lambda wave: .002) + np.testing.assert_allclose(image.xi, [1.9, 2., 2.1]) + centroid = np.sum(image.image[1] * image.lam) / image.image[1].sum() + assert centroid == pytest.approx(1.0002, abs=2e-6) + + +def test_translating_an_ifu_slice_preserves_its_local_spectrum(): + fov = line_cube(.025) + centered = XiLamImage(fov, lambda wave: .00001, slit_dispersion=lambda wave: .002) + fov.cube.header['CRVAL2'] += 2. + fov.meta['eta_min'] += 2*u.arcsec + fov.meta['eta_max'] += 2*u.arcsec + translated = XiLamImage(fov, lambda wave: .00001, slit_dispersion=lambda wave: .002) + np.testing.assert_allclose(translated.image, centered.image, atol=1e-7) + + +def test_missing_aperture_geometry_fails(): + fov = line_cube(.025) + del fov.meta['eta_min'] + with pytest.raises(KeyError, match='eta_min'): + XiLamImage(fov, lambda wave: .00001, slit_dispersion=lambda wave: .002) + + +def test_detector_jacobian_uses_physical_pixel_area(): + header = fits.Header({'CTYPE1D': 'LINEAR', 'CTYPE2D': 'LINEAR', 'CUNIT1D': 'm', 'CUNIT2D': 'm', 'CDELT1D': .000015, 'CDELT2D': .000030, 'CRPIX1D': 1., 'CRPIX2D': 1., 'CRVAL1D': 0., 'CRVAL2D': 0.}) + assert SpectralTrace._detector_pixel_area(header) == pytest.approx(.015*.030) + + +def test_rectified_interpolation_reproduces_a_linear_density_field(): + fov = line_cube(.1) + wave = .995 + np.arange(2001)*fov.cube.header['CDELT3'] + xi = np.arange(3)*.1 + fov.cube.data[:] = 2 + 3*xi[None, None, :] + 5*wave[:, None, None] + image = XiLamImage(fov, lambda w: 0., slit_dispersion=lambda w: 0.) + positions, wavelengths = np.array([.03, .17]), np.array([.998, 1.003]) + expected = 2 + 3*positions + 5*wavelengths + np.testing.assert_allclose(image.interp(positions, wavelengths, grid=False), expected, rtol=2e-6) diff --git a/scopesim/tests/tests_effects/test_local_skycalc.py b/scopesim/tests/tests_effects/test_local_skycalc.py new file mode 100644 index 000000000..8baab5ba1 --- /dev/null +++ b/scopesim/tests/tests_effects/test_local_skycalc.py @@ -0,0 +1,24 @@ +import numpy as np +import pytest +from astropy import units as u +from astropy.table import Table +from scopesim.effects import SkycalcTERCurve + + +@pytest.mark.parametrize('wave_unit', [u.nm, u.um]) +def test_local_skycalc_respects_wavelength_unit(tmp_path, wave_unit): + wave = np.array([.5, 1., 2.]) * u.um + table = Table([wave.to(wave_unit), [.5, .8, .9], [100., 200., 300.]], names=['lam', 'trans', 'flux']) + path = tmp_path / 'sky.fits' + table.write(path) + sky = SkycalcTERCurve(use_local_skycalc_file=str(path)) + np.testing.assert_allclose(sky.skycalc_table['wavelength'], wave.value) + assert sky.surface.transmission(1*u.um).value == pytest.approx(.8) + + + +def test_skycalc_missing_wavelength_unit_fails(tmp_path): + path = tmp_path / 'sky.fits' + Table([[500., 1000., 2000.], [.5, .8, .9], [100., 200., 300.]], names=['lam', 'trans', 'flux']).write(path) + with pytest.raises(TypeError, match='Unit'): + SkycalcTERCurve(use_local_skycalc_file=str(path)) diff --git a/scopesim/tests/tests_effects/test_registered_background.py b/scopesim/tests/tests_effects/test_registered_background.py new file mode 100644 index 000000000..54606b9a5 --- /dev/null +++ b/scopesim/tests/tests_effects/test_registered_background.py @@ -0,0 +1,32 @@ +import numpy as np +import pytest +from scopesim.detector import Detector +from scopesim.optics.image_plane_utils import header_from_list_of_xy +from astropy.io import fits +from scopesim.effects import ImagePlaneBackground +from scopesim.optics.image_plane import ImagePlane + + +def detector(nx=120, ny=80): + header = header_from_list_of_xy([-.015*nx/2, .015*nx/2], [-.015*ny/2, .015*ny/2], .015, 'D') + return Detector(header) + + + +def test_background_uses_current_named_calibration(tmp_path): + path = tmp_path / 'background.fits' + fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(np.full((80, 120), .15), name='0.9'), fits.ImageHDU(np.full((80, 120), .04), name='0.9JH')]).writeto(path) + cmds = {'!INST.slit': '0.9'} + effect = ImagePlaneBackground(filename=str(path), extname='!INST.slit', cmds=cmds) + plane = ImagePlane(detector().header) + effect.apply_to(plane) + np.testing.assert_allclose(plane.data, .15) + cmds['!INST.slit'] = '0.9JH' + plane.data[:] = 0 + effect.apply_to(plane) + np.testing.assert_allclose(plane.data, .04) + with pytest.raises(ValueError, match='shape'): + effect.apply_to(ImagePlane(detector(60, 40).header)) + cmds['!INST.slit'] = 'missing' + with pytest.raises(KeyError): + effect.apply_to(plane) diff --git a/scopesim/utils.py b/scopesim/utils.py index 36b2e5080..a6c60801d 100644 --- a/scopesim/utils.py +++ b/scopesim/utils.py @@ -610,6 +610,9 @@ def from_currsys(item, cmds=None): item.meta = item_meta if isinstance(item, np.ndarray) and not isinstance(item, u.Quantity): + # Numeric calibration arrays cannot contain command references. + if item.dtype.kind in "biufc": + return item.copy() item = np.array([from_currsys(x, cmds) for x in item]) if isinstance(item, list): @@ -1261,4 +1264,4 @@ def _get_next_moon_cached(moontype="full", today=None, location_key=None): nextm = times[mask][np.argmax(moonalts[mask])] else: raise ValueError(f"Invalid moon type: {moontype}, should be 'full', 'new' or 'half'.") - return nextm \ No newline at end of file + return nextm