From aa4c73b6865d3df3c47885a2779c3560a57b976b Mon Sep 17 00:00:00 2001 From: oczoske Date: Tue, 5 May 2026 13:12:16 +0200 Subject: [PATCH 01/11] Remove offset hack --- scopesim/effects/spectral_trace_list.py | 25 +++++++++++-------- scopesim/effects/spectral_trace_list_utils.py | 11 +++++--- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index d455fc08..0fdeb3a7 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -245,16 +245,21 @@ def apply_to(self, obj, **kwargs): obj.cube = obj.make_hdu() # Check whether an offset slit is used. If so, recompute spectral traces. - offset_x = obj.cube.header["CRVAL1D"] - offset_y = obj.cube.header["CRVAL2D"] - if (offset_x != self.meta["offset_x"] or - offset_y != self.meta["offset_y"]): - logger.debug("Recomputing spectral traces for offset (%.1g, %.1g)", - offset_x, offset_y) - self.meta["offset_x"] = offset_x - self.meta["offset_y"] = offset_y - self.make_spectral_traces() - self.update_meta() + # TODO: this needs to be removed. If transforms are available that + # can be applied to the actual slit position, they ought be used + # without modification. If they are not available but need to be + # hacked, then that hack needs to be defined via the yaml and only + # be used in a particular case, not in general. + #offset_x = obj.cube.header["CRVAL1D"] + #offset_y = obj.cube.header["CRVAL2D"] + #if (offset_x != self.meta["offset_x"] or + # offset_y != self.meta["offset_y"]): + # logger.debug("Recomputing spectral traces for offset (%.1g, %.1g)", + # offset_x, offset_y) + # self.meta["offset_x"] = offset_x + # self.meta["offset_y"] = offset_y + # self.make_spectral_traces() + # self.update_meta() spt = self.spectral_traces[obj.trace_id] obj.hdu = spt.map_spectra_to_focal_plane(obj) diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py index 7541f9cb..34371889 100644 --- a/scopesim/effects/spectral_trace_list_utils.py +++ b/scopesim/effects/spectral_trace_list_utils.py @@ -122,17 +122,20 @@ def compute_interpolation_functions(self): `xi` (spatial coordinate along the slit, in arcsec) and `lam` (wavelength, in um). - The interpolation functions include a shift in the focal-plane - coordinates, determined from the CRVAL of the source FOV. + #The interpolation functions include a shift in the focal-plane + #coordinates, determined from the CRVAL of the source FOV. """ - x_arr = self.table[self.meta["x_colname"]] + self.meta["offset_x"] - y_arr = self.table[self.meta["y_colname"]] + self.meta["offset_y"] + x_arr = self.table[self.meta["x_colname"]] + y_arr = self.table[self.meta["y_colname"]] + xi_arr = self.table[self.meta["s_colname"]] lam_arr = self.table[self.meta["wave_colname"]] self.wave_min = quantify(np.min(lam_arr), u.um).value self.wave_max = quantify(np.max(lam_arr), u.um).value + # TODO There should be an option to include pre- and + # posttransforms. But how should they be defined? 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) From 43dda59f9255efd397bf6f001c889593d11ea1e5 Mon Sep 17 00:00:00 2001 From: oczoske Date: Thu, 7 May 2026 16:15:03 +0200 Subject: [PATCH 02/11] Different approach to offset slits --- scopesim/effects/apertures.py | 8 ++++++++ scopesim/effects/spectral_trace_list.py | 19 ++++++++++++++++--- scopesim/effects/spectral_trace_list_utils.py | 8 ++++++-- 3 files changed, 30 insertions(+), 5 deletions(-) diff --git a/scopesim/effects/apertures.py b/scopesim/effects/apertures.py index fd6984a5..d704265b 100644 --- a/scopesim/effects/apertures.py +++ b/scopesim/effects/apertures.py @@ -129,6 +129,14 @@ def apply_to(self, obj, **kwargs): vol["meta"]["xi_min"] = min(x) * u.arcsec vol["meta"]["xi_max"] = max(x) * u.arcsec + if "offset_type" in self.meta: + print("Yes, got it: ", self.meta["offset_type"]) + vol["meta"]["offset_type"] = self.meta["offset_type"] + vol["meta"]["offset_x"] = self.meta["offset_x"] + vol["meta"]["offset_y"] = self.meta["offset_y"] + else: + vol["meta"]["offset_type"] = "none" + print("What is offset_type?") return obj # Outdated. Remove when removing all old FOVManager code from effects diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index 0fdeb3a7..3ec77f55 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -16,7 +16,8 @@ from .effects import Effect from .ter_curves import FilterCurve -from .spectral_trace_list_utils import SpectralTrace, make_image_interpolations +from .spectral_trace_list_utils import (SpectralTrace, make_image_interpolations, + det_offset) from ..optics.image_plane_utils import header_from_list_of_xy from ..optics.fov import FieldOfView from ..optics.fov_volume_list import FovVolumeList @@ -125,8 +126,6 @@ def __init__(self, **kwargs): "x_colname": "x", "y_colname": "y", "s_colname": "s", - "offset_x": 0, # [mm] in detector plane - "offset_y": 0, # [mm] in detector plane "wave_colname": "wavelength", "center_on_wave_mid": False, "dwave": 0.002, # [um] for finding the best fit dispersion @@ -223,6 +222,7 @@ def apply_to(self, obj, **kwargs): aperture_id=vol["aperture_id"]) for ex_vol in extracted_vols: + print("Extracted volume:", ex_vol["meta"]["offset_type"]) ex_vol["meta"].update(vol) ex_vol["meta"].pop("wave_min") ex_vol["meta"].pop("wave_max") @@ -262,8 +262,21 @@ def apply_to(self, obj, **kwargs): # self.update_meta() spt = self.spectral_traces[obj.trace_id] + remove_pre = False + if ("offset_type" in obj.meta and + obj.meta["offset_type"] == "detector"): + spt.xy2lam.pretransform_x = (det_offset, + {"offset": obj.meta["offset_x"]}) + spt.xy2lam.pretransform_y = (det_offset, + {"offset": obj.meta["offset_y"]}) + remove_pre = True + obj.hdu = spt.map_spectra_to_focal_plane(obj) + if remove_pre: + spt.xy2lam.pretransform_y = None + spt.xy2xi.pretransform_x = None + logger.debug("%s done", self.display_name) return obj diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py index 34371889..1fc8e999 100644 --- a/scopesim/effects/spectral_trace_list_utils.py +++ b/scopesim/effects/spectral_trace_list_utils.py @@ -49,8 +49,6 @@ class SpectralTrace: "x_colname": "x", "y_colname": "y", "s_colname": "s", - "offset_x": 0, - "offset_y": 0, "wave_colname": "wavelength", "dwave": 0.002, "aperture_id": 0, @@ -114,6 +112,7 @@ def fov_grid(self): return {"wave_min": wave_min, "wave_max": wave_max, "trace_id": self.trace_id, "aperture_id": aperture_id} + def compute_interpolation_functions(self): """ Compute various interpolation functions between slit and focal plane. @@ -134,6 +133,7 @@ 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 + # TODO There should be an option to include pre- and # posttransforms. But how should they be defined? self.xy2xi = Transform2D.fit(x_arr, y_arr, xi_arr) @@ -1045,3 +1045,7 @@ def get_affine_parameters(coords): shears = (np.average(shears, axis=0) * rad2deg) - (90 + rotations) return rotations, shears + +def det_offset(x, offset): + """Apply offset to x""" + return x + offset From a5bd2b2dee5317d1a068f63bef4db88a9f5707d0 Mon Sep 17 00:00:00 2001 From: oczoske Date: Sat, 16 May 2026 12:36:06 +0200 Subject: [PATCH 03/11] Remove stray print that caused test failures --- scopesim/effects/spectral_trace_list.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index 3ec77f55..a28973f1 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -222,7 +222,6 @@ def apply_to(self, obj, **kwargs): aperture_id=vol["aperture_id"]) for ex_vol in extracted_vols: - print("Extracted volume:", ex_vol["meta"]["offset_type"]) ex_vol["meta"].update(vol) ex_vol["meta"].pop("wave_min") ex_vol["meta"].pop("wave_max") From 418acadcd43270dd07a87cfcb4a6015c3d70e9cb Mon Sep 17 00:00:00 2001 From: oczoske Date: Sat, 16 May 2026 18:24:18 +0200 Subject: [PATCH 04/11] Start MicadoIFUSpectralTraceList --- scopesim/effects/__init__.py | 1 + scopesim/effects/micado_ifu_trace_list.py | 63 +++++++++++++++++++++++ scopesim/effects/spectral_trace_list.py | 2 + 3 files changed, 66 insertions(+) create mode 100644 scopesim/effects/micado_ifu_trace_list.py diff --git a/scopesim/effects/__init__.py b/scopesim/effects/__init__.py index da03ffc0..170b4192 100644 --- a/scopesim/effects/__init__.py +++ b/scopesim/effects/__init__.py @@ -8,6 +8,7 @@ from .spectral_trace_list import * from .spectral_efficiency import * from .metis_lms_trace_list import * +from .micado_ifu_trace_list import * from .mosaic_trace_list import * from .surface_list import * from .ter_curves import * diff --git a/scopesim/effects/micado_ifu_trace_list.py b/scopesim/effects/micado_ifu_trace_list.py new file mode 100644 index 00000000..826399e1 --- /dev/null +++ b/scopesim/effects/micado_ifu_trace_list.py @@ -0,0 +1,63 @@ +# -*- coding: utf-8 -*- +"""SpectralTraceList and SpectralTrace for the MICADO IFU.""" + +import copy +import warnings + +from tqdm.auto import tqdm +import numpy as np +from scipy.interpolate import RectBivariateSpline + +from astropy.io import fits +from astropy.io import ascii as ioascii +from astropy.table import Table +from astropy.wcs import WCS +from astropy import units as u + +from ..utils import from_currsys, find_file, quantify, get_logger +from .spectral_trace_list import SpectralTraceList +from .spectral_trace_list_utils import SpectralTrace +from .spectral_trace_list_utils import Transform2D +from .spectral_trace_list_utils import make_image_interpolations +from .apertures import ApertureMask +from .ter_curves import TERCurve +from ..optics.fov import FieldOfView, FieldOfView3D +from ..optics.fov_volume_list import FovVolumeList + +logger = get_logger(__name__) + + +class MicadoIFUSpectralTraceList(SpectralTraceList): + """SpectralTraceList for the MICADO IFU.""" + + _class_params = { + "naxis1": 112, + "nslice": 32, + "slicewidth": 0.012, # arcsec + } + + + def __init__(self, **kwargs): + super().__init__(**kwargs) + + # field of view of the instrument + + def make_spectral_traces(self): + """Make a spectral trace for each combination of order and aperture""" + self.ext_data = self._file[0].header["EDATA"] + self.ext_cat = self._file[0].header["ECAT"] + self.catalog = Table(self._file[self.ext_cat].data) + self.slicelist = self._file["Aperture List"].data + spec_traces = {} + for row in self.catalog: + if row["image_plane_id"] == -99: + continue + params = {col: row[col] for col in row.colnames} + params.update(self.meta) + hdu = self._file[row["extension_id"]] + for apid in self.slicelist["id"]: + specid = f"{row['description']}_{apid:02d}" + params["aperture_id"] = apid + spec_traces[specid] = SpectralTrace(hdu, **params) + + self.spectral_traces = spec_traces diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index a28973f1..48e5abef 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -148,6 +148,8 @@ def make_spectral_traces(self): self.catalog = Table(self._file[self.ext_cat].data) spec_traces = {} for row in self.catalog: + if row["image_plane_id"] == -99: + continue params = {col: row[col] for col in row.colnames} params.update(self.meta) hdu = self._file[row["extension_id"]] From 31ce32c66b5ef4d46cb1332db56857d872e9c55e Mon Sep 17 00:00:00 2001 From: oczoske Date: Mon, 18 May 2026 14:23:47 +0200 Subject: [PATCH 05/11] Split apply_to into separate methods for setup and observe --- scopesim/effects/metis_lms_trace_list.py | 173 +++++++++++----------- scopesim/effects/spectral_trace_list.py | 175 +++++++++++++---------- 2 files changed, 187 insertions(+), 161 deletions(-) diff --git a/scopesim/effects/metis_lms_trace_list.py b/scopesim/effects/metis_lms_trace_list.py index bf8c0d88..e6770a90 100644 --- a/scopesim/effects/metis_lms_trace_list.py +++ b/scopesim/effects/metis_lms_trace_list.py @@ -53,93 +53,98 @@ def __init__(self, **kwargs): self.slicelist["top"].max() - self.slicelist["bottom"].min()]) - def apply_to(self, obj, **kwargs): - """See parent docstring.""" - if isinstance(obj, FovVolumeList): - # Create a single volume that covers the aperture and - # the maximum wavelength range of LMS - volumes = [spectral_trace.fov_grid() - for spectral_trace in self.spectral_traces.values()] - - wave_min = min(vol["wave_min"] for vol in volumes) - wave_max = max(vol["wave_max"] for vol in volumes) - extracted_vols = obj.extract(axes=["wave"], - edges=([[wave_min, wave_max]])) - obj.volumes = extracted_vols - - if isinstance(obj, FieldOfView): - # Application to field of view - logger.debug("Executing %s, FoV", self.meta['name']) - if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: - obj.cube = obj.hdu - elif obj.hdu is None and obj.cube is None: - obj.cube = obj.make_hdu() - - fovcube = obj.cube.data - n_z, n_y, n_x = fovcube.shape - fovwcs = WCS(obj.cube.header) - # Make this linear to avoid jump at RA 0 deg - fovwcs.wcs.ctype = ["LINEAR", "LINEAR", fovwcs.wcs.ctype[2]] - fovwcs_spat = fovwcs.sub(2) - ny_slice = self.meta["slice_samples"] - - # Spatial pixel coordinates - xslice, yslice = np.meshgrid(np.arange(n_x), - np.arange(ny_slice)) - - fovimage = np.zeros((obj.detector_header["NAXIS2"], - obj.detector_header["NAXIS1"]), - dtype=np.float32) - - for sptid, spt in tqdm(self.spectral_traces.items(), - desc=" Spectral Traces", position=2): - ymin = spt.meta["fov"]["y_min"] - ymax = spt.meta["fov"]["y_max"] - - slicewcs = fovwcs.deepcopy() - - slicewcs.wcs.ctype = ["LINEAR", "LINEAR", - slicewcs.wcs.ctype[2]] - slicewcs.wcs.crpix[1] = (ny_slice + 1) / 2 - slicewcs.wcs.crval[1] = (ymin + ymax) / 2 / 3600 - slicewcs.wcs.cdelt[1] = (ymax - ymin) / ny_slice / 3600 - slicewcs_spat = slicewcs.sub(2) - - # World coordinates for the slice - xworld, yworld = slicewcs_spat.all_pix2world(xslice, - yslice, 0) - # FOV pixel coordinates for the slice - xfov, yfov = fovwcs_spat.all_world2pix(xworld, yworld, 0) - - slicecube = np.zeros((n_z, ny_slice, n_x)) - for islice in range(n_z): - ifov = RectBivariateSpline(np.arange(n_y), - np.arange(n_x), - fovcube[islice], kx=1, ky=1) - slicecube[islice] = ifov(yfov, xfov, grid=False) - - slicefov = FieldOfView3D(obj.header, - [obj.meta["wave_min"], - obj.meta["wave_max"]]) - slicefov.detector_header = obj.detector_header - slicefov.meta["xi_min"] = obj.meta["xi_min"] - slicefov.meta["xi_max"] = obj.meta["xi_max"] - slicefov.meta["trace_id"] = sptid - slicefov.cube = fits.ImageHDU(header=slicewcs.to_header(), - data=slicecube) - # slicefov.cube.writeto(f"slicefov_{sptid}.fits", overwrite=True) - slicefov.hdu = spt.map_spectra_to_focal_plane(slicefov) - if slicefov.hdu is not None: - sxmin = slicefov.hdu.header["XMIN"] - sxmax = slicefov.hdu.header["XMAX"] - symin = slicefov.hdu.header["YMIN"] - symax = slicefov.hdu.header["YMAX"] - fovimage[symin:symax, sxmin:sxmax] += slicefov.hdu.data - - obj.hdu = fits.ImageHDU(data=fovimage, header=obj.detector_header) + # Overrides for the two apply_to() methods + def apply_to_fovvolumelist(self, obj): + """Setup + Create a single volume that covers the aperture and + the maximum wavelength range of LMS. + """ + logger.debug("Executing %s: Setup", self.meta['name']) + volumes = [spectral_trace.fov_grid() + for spectral_trace in self.spectral_traces.values()] + + wave_min = min(vol["wave_min"] for vol in volumes) + wave_max = max(vol["wave_max"] for vol in volumes) + extracted_vols = obj.extract(axes=["wave"], + edges=([[wave_min, wave_max]])) + obj.volumes = extracted_vols + return obj + + def apply_to_fov(self, obj): + """Application to field of view""" + + logger.debug("Executing %s, FoV", self.meta['name']) + if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: + obj.cube = obj.hdu + elif obj.hdu is None and obj.cube is None: + obj.cube = obj.make_hdu() + + fovcube = obj.cube.data + n_z, n_y, n_x = fovcube.shape + fovwcs = WCS(obj.cube.header) + # Make this linear to avoid jump at RA 0 deg + fovwcs.wcs.ctype = ["LINEAR", "LINEAR", fovwcs.wcs.ctype[2]] + fovwcs_spat = fovwcs.sub(2) + ny_slice = self.meta["slice_samples"] + + # Spatial pixel coordinates + xslice, yslice = np.meshgrid(np.arange(n_x), + np.arange(ny_slice)) + + fovimage = np.zeros((obj.detector_header["NAXIS2"], + obj.detector_header["NAXIS1"]), + dtype=np.float32) + + for sptid, spt in tqdm(self.spectral_traces.items(), + desc=" Spectral Traces", position=2): + ymin = spt.meta["fov"]["y_min"] + ymax = spt.meta["fov"]["y_max"] + + slicewcs = fovwcs.deepcopy() + + slicewcs.wcs.ctype = ["LINEAR", "LINEAR", + slicewcs.wcs.ctype[2]] + slicewcs.wcs.crpix[1] = (ny_slice + 1) / 2 + slicewcs.wcs.crval[1] = (ymin + ymax) / 2 / 3600 + slicewcs.wcs.cdelt[1] = (ymax - ymin) / ny_slice / 3600 + slicewcs_spat = slicewcs.sub(2) + + # World coordinates for the slice + xworld, yworld = slicewcs_spat.all_pix2world(xslice, + yslice, 0) + # FOV pixel coordinates for the slice + xfov, yfov = fovwcs_spat.all_world2pix(xworld, yworld, 0) + + slicecube = np.zeros((n_z, ny_slice, n_x)) + for islice in range(n_z): + ifov = RectBivariateSpline(np.arange(n_y), + np.arange(n_x), + fovcube[islice], kx=1, ky=1) + slicecube[islice] = ifov(yfov, xfov, grid=False) + + slicefov = FieldOfView3D(obj.header, + [obj.meta["wave_min"], + obj.meta["wave_max"]]) + slicefov.detector_header = obj.detector_header + slicefov.meta["xi_min"] = obj.meta["xi_min"] + slicefov.meta["xi_max"] = obj.meta["xi_max"] + slicefov.meta["trace_id"] = sptid + slicefov.cube = fits.ImageHDU(header=slicewcs.to_header(), + data=slicecube) + # slicefov.cube.writeto(f"slicefov_{sptid}.fits", overwrite=True) + slicefov.hdu = spt.map_spectra_to_focal_plane(slicefov) + if slicefov.hdu is not None: + sxmin = slicefov.hdu.header["XMIN"] + sxmax = slicefov.hdu.header["XMAX"] + symin = slicefov.hdu.header["YMIN"] + symax = slicefov.hdu.header["YMAX"] + fovimage[symin:symax, sxmin:sxmax] += slicefov.hdu.data + + obj.hdu = fits.ImageHDU(data=fovimage, header=obj.detector_header) return obj + def make_spectral_traces(self): """Compute the transformations by interpolation.""" # nslice = len(self._file["Aperture List"].data) diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index 48e5abef..c4d3fc92 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -184,6 +184,7 @@ def update_meta(self): self.meta["y_min"] = min(ylim) self.meta["y_max"] = max(ylim) + def apply_to(self, obj, **kwargs): """ Interface between ``FieldOfView`` and ``SpectralTraceList``. @@ -198,89 +199,109 @@ def apply_to(self, obj, **kwargs): list, identified by meta["trace_id"]. """ if isinstance(obj, FovVolumeList): - logger.debug("%s applied to %s", self.display_name, - obj.__class__.__name__) - # Setup of FieldOfView object - # volumes = [spectral_trace.fov_grid() - # for spectral_trace in self.spectral_traces.values()] - - new_vols_list = [] - - # for vol in volumes: - for spt in self.spectral_traces.values(): - vol = spt.fov_grid() - wave_edges = [vol["wave_min"], vol["wave_max"]] - if "x_min" in vol: - x_edges = [vol["x_min"], vol["x_max"]] - y_edges = [vol["y_min"], vol["y_max"]] - extracted_vols = obj.extract( - axes=["wave", "x", "y"], - edges=(wave_edges, x_edges, y_edges), - aperture_id=vol["aperture_id"]) - else: - extracted_vols = obj.extract( - axes=["wave"], - edges=(wave_edges, ), - aperture_id=vol["aperture_id"]) - - for ex_vol in extracted_vols: - ex_vol["meta"].update(vol) - ex_vol["meta"].pop("wave_min") - ex_vol["meta"].pop("wave_max") - new_vols_list.extend(extracted_vols) - - obj.volumes = new_vols_list + obj = self.apply_to_fovvolumelist(obj) if isinstance(obj, FieldOfView): - logger.debug("%s applied to %s", self.display_name, - obj.__class__.__name__) - # Application to field of view - if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: - obj.cube = obj.hdu - elif obj.hdu is not None and obj.hdu.header["NAXIS"] == 2: - # todo: catch the case of obj.hdu.header["NAXIS"] == 2 - # for MAAT - pass - elif obj.hdu is None and obj.cube is None: - logger.info("Making cube") - obj.cube = obj.make_hdu() - - # Check whether an offset slit is used. If so, recompute spectral traces. - # TODO: this needs to be removed. If transforms are available that - # can be applied to the actual slit position, they ought be used - # without modification. If they are not available but need to be - # hacked, then that hack needs to be defined via the yaml and only - # be used in a particular case, not in general. - #offset_x = obj.cube.header["CRVAL1D"] - #offset_y = obj.cube.header["CRVAL2D"] - #if (offset_x != self.meta["offset_x"] or - # offset_y != self.meta["offset_y"]): - # logger.debug("Recomputing spectral traces for offset (%.1g, %.1g)", - # offset_x, offset_y) - # self.meta["offset_x"] = offset_x - # self.meta["offset_y"] = offset_y - # self.make_spectral_traces() - # self.update_meta() - - spt = self.spectral_traces[obj.trace_id] - remove_pre = False - if ("offset_type" in obj.meta and - obj.meta["offset_type"] == "detector"): - spt.xy2lam.pretransform_x = (det_offset, - {"offset": obj.meta["offset_x"]}) - spt.xy2lam.pretransform_y = (det_offset, - {"offset": obj.meta["offset_y"]}) - remove_pre = True - - obj.hdu = spt.map_spectra_to_focal_plane(obj) - - if remove_pre: - spt.xy2lam.pretransform_y = None - spt.xy2xi.pretransform_x = None + obj = self.apply_to_fov(obj) logger.debug("%s done", self.display_name) return obj + def apply_to_fovvolumelist(self, obj): + """ + Setup of SpectralTraceList + + During setup of the required FieldOfView objects, the + SpectralTraceList is asked for the source space volumes that + it requires (spatial limits and wavelength limits). + """ + logger.debug("%s applied to %s", self.display_name, + obj.__class__.__name__) + + new_vols_list = [] + + for spt in self.spectral_traces.values(): + vol = spt.fov_grid() + wave_edges = [vol["wave_min"], vol["wave_max"]] + if "x_min" in vol: + x_edges = [vol["x_min"], vol["x_max"]] + y_edges = [vol["y_min"], vol["y_max"]] + extracted_vols = obj.extract( + axes=["wave", "x", "y"], + edges=(wave_edges, x_edges, y_edges), + aperture_id=vol["aperture_id"]) + else: + extracted_vols = obj.extract( + axes=["wave"], + edges=(wave_edges, ), + aperture_id=vol["aperture_id"]) + + for ex_vol in extracted_vols: + ex_vol["meta"].update(vol) + ex_vol["meta"].pop("wave_min") + ex_vol["meta"].pop("wave_max") + new_vols_list.extend(extracted_vols) + + obj.volumes = new_vols_list + return obj + + + def apply_to_fov(self, obj): + """Apply to FieldOfView + + During "observation" the method is passed a single FieldOfView + object and applies the mapping to the image plane to it. + The FieldOfView object is associated to one SpectralTrace from the + list, identified by meta["trace_id"]. + """ + logger.debug("%s applied to %s %s", self.display_name, + obj.__class__.__name__, obj.trace_id) + + if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: + obj.cube = obj.hdu + elif obj.hdu is not None and obj.hdu.header["NAXIS"] == 2: + # todo: catch the case of obj.hdu.header["NAXIS"] == 2 + # for MAAT + pass + elif obj.hdu is None and obj.cube is None: + logger.info("Making cube") + obj.cube = obj.make_hdu() + + # Check whether an offset slit is used. If so, recompute spectral traces. + # TODO: this needs to be removed. If transforms are available that + # can be applied to the actual slit position, they ought be used + # without modification. If they are not available but need to be + # hacked, then that hack needs to be defined via the yaml and only + # be used in a particular case, not in general. + #offset_x = obj.cube.header["CRVAL1D"] + #offset_y = obj.cube.header["CRVAL2D"] + #if (offset_x != self.meta["offset_x"] or + # offset_y != self.meta["offset_y"]): + # logger.debug("Recomputing spectral traces for offset (%.1g, %.1g)", + # offset_x, offset_y) + # self.meta["offset_x"] = offset_x + # self.meta["offset_y"] = offset_y + # self.make_spectral_traces() + # self.update_meta() + + spt = self.spectral_traces[obj.trace_id] + remove_pre = False + if ("offset_type" in obj.meta and + obj.meta["offset_type"] == "detector"): + spt.xy2lam.pretransform_x = (det_offset, + {"offset": obj.meta["offset_x"]}) + spt.xy2lam.pretransform_y = (det_offset, + {"offset": obj.meta["offset_y"]}) + remove_pre = True + + obj.hdu = spt.map_spectra_to_focal_plane(obj) + + if remove_pre: + spt.xy2lam.pretransform_y = None + spt.xy2xi.pretransform_x = None + + return obj + @property def footprint(self): """Return the footprint of the entire SpectralTraceList.""" From 3375ec53ec2ab2df8bf42fa770d1772a4b6b9ca8 Mon Sep 17 00:00:00 2001 From: oczoske Date: Mon, 18 May 2026 19:01:40 +0200 Subject: [PATCH 06/11] Remove prints --- scopesim/effects/apertures.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/scopesim/effects/apertures.py b/scopesim/effects/apertures.py index d704265b..b8ca2ce2 100644 --- a/scopesim/effects/apertures.py +++ b/scopesim/effects/apertures.py @@ -118,6 +118,7 @@ def apply_to(self, obj, **kwargs): """See parent docstring.""" if isinstance(obj, FovVolumeList): logger.debug("Executing %s, FoV setup", self.meta['name']) + x = quantity_from_table("x", self.table, u.arcsec).to_value(u.arcsec) y = quantity_from_table("y", self.table, @@ -130,13 +131,11 @@ def apply_to(self, obj, **kwargs): vol["meta"]["xi_max"] = max(x) * u.arcsec if "offset_type" in self.meta: - print("Yes, got it: ", self.meta["offset_type"]) vol["meta"]["offset_type"] = self.meta["offset_type"] vol["meta"]["offset_x"] = self.meta["offset_x"] vol["meta"]["offset_y"] = self.meta["offset_y"] else: vol["meta"]["offset_type"] = "none" - print("What is offset_type?") return obj # Outdated. Remove when removing all old FOVManager code from effects From 385cade20f40b9f977aae01704975ded02e5f330 Mon Sep 17 00:00:00 2001 From: oczoske Date: Mon, 18 May 2026 19:02:00 +0200 Subject: [PATCH 07/11] Set ImagePlane to 32bit --- scopesim/optics/image_plane.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scopesim/optics/image_plane.py b/scopesim/optics/image_plane.py index 794c03fb..b58ac50c 100644 --- a/scopesim/optics/image_plane.py +++ b/scopesim/optics/image_plane.py @@ -68,7 +68,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) + image = zeros_from_header(header, dtype=np.float32) self.hdu = fits.ImageHDU(data=image, header=header) self.hdu.header["BUNIT"] = "ph s-1" # photons per second (per pixel) From b2de962a16f2f55e50d97b60110aab46703948e7 Mon Sep 17 00:00:00 2001 From: oczoske Date: Mon, 18 May 2026 19:02:20 +0200 Subject: [PATCH 08/11] MICADO IFU works --- scopesim/effects/micado_ifu_trace_list.py | 153 +++++++++++++++++++++- 1 file changed, 149 insertions(+), 4 deletions(-) diff --git a/scopesim/effects/micado_ifu_trace_list.py b/scopesim/effects/micado_ifu_trace_list.py index 826399e1..c8236e79 100644 --- a/scopesim/effects/micado_ifu_trace_list.py +++ b/scopesim/effects/micado_ifu_trace_list.py @@ -34,6 +34,8 @@ class MicadoIFUSpectralTraceList(SpectralTraceList): "naxis1": 112, "nslice": 32, "slicewidth": 0.012, # arcsec + "wave_min": 1.5, + "wave_max": 2.5, # todo: do not hardcode } @@ -42,6 +44,104 @@ def __init__(self, **kwargs): # field of view of the instrument + + def apply_to_fovvolumelist(self, obj): + """ + Setup + + Create a single volume that covers the aperture and the maximum + wavelength range of the MICADO IFU. + """ + logger.debug("Executing %s: Setup", self.meta['name']) + volumes = [spectral_trace.fov_grid() + for spectral_trace in self.spectral_traces.values()] + + wave_min = min(vol["wave_min"] for vol in volumes) + wave_max = max(vol["wave_max"] for vol in volumes) + extracted_vols = obj.extract(axes=["wave"], + edges=([[wave_min, wave_max]])) + obj.volumes = extracted_vols + self.wave_min = wave_min + self.wave_max = wave_max + return obj + + def apply_to_fov(self, obj): + """Apply to FieldOfView""" + + logger.debug("Executing %s, FoV", self.meta['name']) + if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: + obj.cube = obj.hdu + elif obj.hdu is None and obj.cube is None: + obj.cube = obj.make_hdu() + + fovcube = obj.cube.data + n_z, n_y, n_x = fovcube.shape + fovwcs = WCS(obj.cube.header) + # Make this linear to avoid jump at RA 0 deg + fovwcs.wcs.ctype = ["LINEAR", "LINEAR", fovwcs.wcs.ctype[2]] + fovwcs_spat = fovwcs.sub(2) + ny_slice = self.meta["slice_samples"] + + # Spatial pixel coordinates + xslice, yslice = np.meshgrid(np.arange(n_x), + np.arange(ny_slice)) + + fovimage = np.zeros((obj.detector_header["NAXIS2"], + obj.detector_header["NAXIS1"]), + dtype=np.float32) + + for sptid, spt in tqdm(self.spectral_traces.items(), + desc=" Spectral Traces", position=2): + ymin = spt.meta["fov"]["y_min"] + ymax = spt.meta["fov"]["y_max"] + + + slicewcs = fovwcs.deepcopy() + + slicewcs.wcs.ctype = ["LINEAR", "LINEAR", + slicewcs.wcs.ctype[2]] + slicewcs.wcs.crpix[1] = (ny_slice + 1) / 2 + slicewcs.wcs.crval[1] = (ymin + ymax) / 2 / 3600 + slicewcs.wcs.cdelt[1] = (ymax - ymin) / ny_slice / 3600 + slicewcs_spat = slicewcs.sub(2) + + # World coordinates for the slice + xworld, yworld = slicewcs_spat.all_pix2world(xslice, + yslice, 0) + # FOV pixel coordinates for the slice + xfov, yfov = fovwcs_spat.all_world2pix(xworld, yworld, 0) + + slicecube = np.zeros((n_z, ny_slice, n_x)) + for islice in range(n_z): + ifov = RectBivariateSpline(np.arange(n_y), + np.arange(n_x), + fovcube[islice], kx=1, ky=1) + slicecube[islice] = ifov(yfov, xfov, grid=False) + + slicefov = FieldOfView3D(obj.header, + [obj.meta["wave_min"], + obj.meta["wave_max"]]) + slicefov.detector_header = obj.detector_header + slicefov.meta["xi_min"] = spt.meta["fov"]["xi_min"] + slicefov.meta["xi_max"] = spt.meta["fov"]["xi_max"] + slicefov.meta["trace_id"] = sptid + slicefov.cube = fits.ImageHDU(header=slicewcs.to_header(), + data=slicecube) + # slicefov.cube.writeto(f"slicefov_{sptid}.fits", overwrite=True) + slicefov.hdu = spt.map_spectra_to_focal_plane(slicefov) + if slicefov.hdu is not None: + sxmin = slicefov.hdu.header["XMIN"] + sxmax = slicefov.hdu.header["XMAX"] + symin = slicefov.hdu.header["YMIN"] + symax = slicefov.hdu.header["YMAX"] + fovimage[symin:symax, sxmin:sxmax] += slicefov.hdu.data + + obj.hdu = fits.ImageHDU(data=fovimage, header=obj.detector_header) + return obj + + + + def make_spectral_traces(self): """Make a spectral trace for each combination of order and aperture""" self.ext_data = self._file[0].header["EDATA"] @@ -54,10 +154,55 @@ def make_spectral_traces(self): continue params = {col: row[col] for col in row.colnames} params.update(self.meta) - hdu = self._file[row["extension_id"]] - for apid in self.slicelist["id"]: + trace_hdu = self._file[row["extension_id"]] + + for spslice, apid in enumerate(self.slicelist["id"]): specid = f"{row['description']}_{apid:02d}" - params["aperture_id"] = apid - spec_traces[specid] = SpectralTrace(hdu, **params) + thistrace = MicadoIFUSpectralTrace(trace_hdu, self.slicelist, spslice, + **params) + thistrace.meta["trace_id"] = specid + spec_traces[specid] = thistrace self.spectral_traces = spec_traces + + + +class MicadoIFUSpectralTrace(SpectralTrace): + """SpectralTrace for the MICADO IFU""" + + _class_params = { + "naxis1": 112, + "nslice": 32, + "slicewidth": 0.012, # arcsec + } + + + def __init__(self, trace_tbl, aplist, spslice, cmds=None, **kwargs): + super().__init__(trace_tbl, **kwargs) + # Provisional + self.aplist = aplist + self.meta["slice"] = spslice + self.meta["fov"] = self.fov_grid() + + def fov_grid(self): + """ + Provide information on the source space volume required by the effect. + """ + aperture = self.aplist[self.meta["slice"]] + x_min = aperture["left"] + x_max = aperture["right"] + y_min = aperture["bottom"] + y_max = aperture["top"] + xi_min = aperture["xi1"] + xi_max = aperture["xi2"] + + # ..todo: just a hack - xi and x are the same except xi is a quantity + xi_min = quantify(xi_min, u.arcsec) + xi_max = quantify(xi_max, u.arcsec) + + return {"x_min": x_min, "x_max": x_max, + "y_min": y_min, "y_max": y_max, + "xi_min": xi_min, "xi_max": xi_max, + "wave_min": self.meta["wave_min"], + "wave_max": self.meta["wave_max"], + "trace_id": self.trace_id} From b63315127a42954c33c43baebbbe1457781ddf04 Mon Sep 17 00:00:00 2001 From: oczoske Date: Tue, 19 May 2026 17:45:51 +0200 Subject: [PATCH 09/11] Staggered ifu slices --- scopesim/effects/micado_ifu_trace_list.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/scopesim/effects/micado_ifu_trace_list.py b/scopesim/effects/micado_ifu_trace_list.py index c8236e79..d01ef003 100644 --- a/scopesim/effects/micado_ifu_trace_list.py +++ b/scopesim/effects/micado_ifu_trace_list.py @@ -18,7 +18,7 @@ from .spectral_trace_list import SpectralTraceList from .spectral_trace_list_utils import SpectralTrace from .spectral_trace_list_utils import Transform2D -from .spectral_trace_list_utils import make_image_interpolations +from .spectral_trace_list_utils import make_image_interpolations, det_offset from .apertures import ApertureMask from .ter_curves import TERCurve from ..optics.fov import FieldOfView, FieldOfView3D @@ -42,8 +42,6 @@ class MicadoIFUSpectralTraceList(SpectralTraceList): def __init__(self, **kwargs): super().__init__(**kwargs) - # field of view of the instrument - def apply_to_fovvolumelist(self, obj): """ @@ -174,16 +172,26 @@ class MicadoIFUSpectralTrace(SpectralTrace): "naxis1": 112, "nslice": 32, "slicewidth": 0.012, # arcsec + "pixsize": 0.015, # mm } def __init__(self, trace_tbl, aplist, spslice, cmds=None, **kwargs): super().__init__(trace_tbl, **kwargs) + + self._set_dispersion(self.wave_min, self.wave_max) # Provisional self.aplist = aplist self.meta["slice"] = spslice self.meta["fov"] = self.fov_grid() + lam = np.linspace(self.wave_min, self.wave_max, 1001) * u.um + off_y = aplist['offset'][spslice] + dlam = np.median(self.dlam_per_pix(lam)) / self.meta["pixsize"] * off_y + self.xy2lam.posttransform = (det_offset, {"offset": dlam}) + self.xilam2x.pretransform_y = (det_offset, {"offset": dlam}) + self.xilam2y.pretransform_y = (det_offset, {"offset": dlam}) + def fov_grid(self): """ Provide information on the source space volume required by the effect. From 8b01497d3d3adec269eeffabadeaa5896c023bbc Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 20 May 2026 13:19:03 +0200 Subject: [PATCH 10/11] Flipped sign of offset --- scopesim/effects/micado_ifu_trace_list.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/scopesim/effects/micado_ifu_trace_list.py b/scopesim/effects/micado_ifu_trace_list.py index d01ef003..cd2595b4 100644 --- a/scopesim/effects/micado_ifu_trace_list.py +++ b/scopesim/effects/micado_ifu_trace_list.py @@ -41,7 +41,8 @@ class MicadoIFUSpectralTraceList(SpectralTraceList): def __init__(self, **kwargs): super().__init__(**kwargs) - + self.wave_min = self.meta["wave_min"] + self.wave_max = self.meta["wave_max"] def apply_to_fovvolumelist(self, obj): """ @@ -188,9 +189,9 @@ def __init__(self, trace_tbl, aplist, spslice, cmds=None, **kwargs): lam = np.linspace(self.wave_min, self.wave_max, 1001) * u.um off_y = aplist['offset'][spslice] dlam = np.median(self.dlam_per_pix(lam)) / self.meta["pixsize"] * off_y - self.xy2lam.posttransform = (det_offset, {"offset": dlam}) - self.xilam2x.pretransform_y = (det_offset, {"offset": dlam}) - self.xilam2y.pretransform_y = (det_offset, {"offset": dlam}) + self.xy2lam.posttransform = (det_offset, {"offset": -dlam}) + self.xilam2x.pretransform_y = (det_offset, {"offset": -dlam}) + self.xilam2y.pretransform_y = (det_offset, {"offset": -dlam}) def fov_grid(self): """ From 10a7bdb9609160ee56d6f998fe9e7dcf52b18e2b Mon Sep 17 00:00:00 2001 From: oczoske Date: Mon, 25 May 2026 06:26:30 +0200 Subject: [PATCH 11/11] offset more systematic --- scopesim/effects/spectral_trace_list.py | 40 +++++++++++++++++++++---- 1 file changed, 35 insertions(+), 5 deletions(-) diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index c4d3fc92..4422e183 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -288,10 +288,17 @@ def apply_to_fov(self, obj): remove_pre = False if ("offset_type" in obj.meta and obj.meta["offset_type"] == "detector"): + self.meta["offset_type"] = "detector" + self.meta["offset_x"] = obj.meta["offset_x"] + self.meta["offset_y"] = obj.meta["offset_y"] spt.xy2lam.pretransform_x = (det_offset, {"offset": obj.meta["offset_x"]}) spt.xy2lam.pretransform_y = (det_offset, {"offset": obj.meta["offset_y"]}) + spt.xilam2x.posttransform = (det_offset, + {"offset": -obj.meta["offset_x"]}) + spt.xilam2y.posttransform = (det_offset, + {"offset": -obj.meta["offset_y"]}) remove_pre = True obj.hdu = spt.map_spectra_to_focal_plane(obj) @@ -299,6 +306,8 @@ def apply_to_fov(self, obj): if remove_pre: spt.xy2lam.pretransform_y = None spt.xy2xi.pretransform_x = None + spt.xilam2x.posttransform = None + spt.xilam2y.posttransform = None return obj @@ -401,15 +410,36 @@ def rectify_traces(self, hdulist, xi_min=None, xi_max=None, interps=None, for i, trace_id in tqdm(enumerate(self.spectral_traces, start=1), desc=" Traces", total=len(self.spectral_traces)): - hdu = self[trace_id].rectify(hdulist, - interps=interps, - bin_width=bin_width, - xi_min=xi_min, xi_max=xi_max, - wave_min=wave_min, wave_max=wave_max) + spt = self[trace_id] + remove_pre = False + if ("offset_type" in self.meta and + self.meta["offset_type"] == "detector"): + spt.xy2lam.pretransform_x = (det_offset, + {"offset": self.meta["offset_x"]}) + spt.xy2lam.pretransform_y = (det_offset, + {"offset": self.meta["offset_y"]}) + spt.xilam2x.posttransform = (det_offset, + {"offset": -self.meta["offset_x"]}) + spt.xilam2y.posttransform = (det_offset, + {"offset": -self.meta["offset_y"]}) + remove_pre = True + + + hdu = spt.rectify(hdulist, + interps=interps, + bin_width=bin_width, + xi_min=xi_min, xi_max=xi_max, + wave_min=wave_min, wave_max=wave_max) if hdu is not None: # ..todo: rectify does not do that yet outhdul.append(hdu) outhdul[0].header[f"EXTNAME{i}"] = trace_id + if remove_pre: + spt.xy2lam.pretransform_y = None + spt.xy2xi.pretransform_x = None + spt.xilam2x.posttransform = None + spt.xilam2y.posttransform = None + outhdul[0].header.update(inhdul[0].header) return outhdul