diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py index dfb28193..2569ec96 100644 --- a/scopesim/effects/spectral_trace_list_utils.py +++ b/scopesim/effects/spectral_trace_list_utils.py @@ -15,7 +15,9 @@ import numpy as np from scipy.interpolate import RectBivariateSpline -from scipy.interpolate import interp1d +from scipy.interpolate import InterpolatedUnivariateSpline, interp1d +from scipy.ndimage import zoom +from matplotlib import pyplot as plt from astropy.table import Table, vstack from astropy.modeling import fitting @@ -24,8 +26,11 @@ from astropy.wcs import WCS from astropy.modeling.models import Polynomial2D -from ..utils import (power_vector, quantify, from_currsys, close_loop, - figure_factory, get_logger) +import pandas as pd + +from ..optics import image_plane_utils as imp_utils +from ..utils import (deriv_polynomial2d, power_vector, interp2, check_keys, + from_currsys, quantify, close_loop, figure_factory, get_logger) logger = get_logger(__name__) @@ -132,7 +137,10 @@ def compute_interpolation_functions(self): 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._xix2y = Transform2D.fit(xi_arr, x_arr, y_arr) self._xiy2x = Transform2D.fit(xi_arr, y_arr, x_arr) + self._xix2lam = Transform2D.fit(xi_arr, x_arr, lam_arr) self._xiy2lam = Transform2D.fit(xi_arr, y_arr, lam_arr) if self.dispersion_axis == "unknown": @@ -423,6 +431,12 @@ def footprint(self, wave_min=None, wave_max=None, xi_min=None, xi_max=None): Minimum and maximum slit position on the sky. If `None`, use the full range that spectral trace is defined on. Float values are interpreted as arcsec. + + Returns + ------- + xy_edges : tuple of lists + (x_edges, y_edges) in [mm] + """ # Define the wavelength range of the footprint. This is a compromise # between the requested range (by method args) and the definition @@ -673,8 +687,12 @@ def __init__(self, fov, dlam_per_pix): # dimensions are set by the slit aperture (n_lam, n_eta, n_xi) = fov.cube.data.shape + m_xi = int(np.round((fov.meta['xi_max'].to(u.arcsec).value - + fov.meta['xi_min'].to(u.arcsec).value) / d_xi)) + assert abs(m_xi - n_xi) <= 1 + # arrays of cube coordinates - cube_xi = d_xi * np.arange(n_xi) + fov.meta["xi_min"].value + cube_xi = d_xi * np.arange(n_xi) + fov.meta["xi_min"].to(u.arcsec).value cube_eta = d_eta * (np.arange(n_eta) - (n_eta - 1) / 2) cube_lam = wcs_lam.all_pix2world(np.arange(n_lam), 1)[0] cube_lam *= u.Unit(wcs_lam.wcs.cunit[0]).to(u.um) @@ -983,8 +1001,6 @@ def make_image_interpolations(hdulist, **kwargs): return interps # ..todo: Check whether the following functions are actually used - - def rolling_median(x, n): """Calculate the rolling median of a sequence for +/- n entries.""" y = [np.median(x[max(0, i-n):min(len(x), i+n+1)]) for i in range(len(x))] @@ -1038,3 +1054,152 @@ def get_affine_parameters(coords): shears = (np.average(shears, axis=0) * rad2deg) - (90 + rotations) return rotations, shears + + +class TracesHDUListGenerator: + def __init__(self, trace_hdus, aperture_ids, image_plane_ids): + """ + + Parameters + ---------- + trace_hdus + aperture_ids + image_plane_ids + """ + self.trace_hdus = trace_hdus + self.aperture_ids = aperture_ids + self.image_plane_ids = image_plane_ids + + def make_hdulist(self, filename=None, overwrite=False): + trace_hdulist = fits.HDUList([self.pri_hdu, self.cat_hdu] + + self.trace_hdus) + + if filename: + trace_hdulist.writeto(filename, overwrite=overwrite) + + return trace_hdulist + + @property + def hdulist(self): + return self.make_hdulist() + + @property + def pri_hdu(self): + pri_hdu = fits.PrimaryHDU() + pri_hdu.header["ECAT"] = 1 + pri_hdu.header["EDATA"] = 2 + + meta = {"author": "", + "source": "", + "descript": "Spectral Trace List", + "date-cre": "", + "date-mod": "", + } + pri_hdu.header.update(meta) + + return pri_hdu + + @property + def cat_hdu(self): + trace_names = [f"TRACE_{i}" for i in range(len(self.trace_hdus))] + cat_table = Table(names=["description", "extension_id", + "aperture_id", "image_plane_id"], + data=[trace_names, + 2 + np.arange(len(self.trace_hdus)), + self.aperture_ids, + self.image_plane_ids + ]) + cat_hdu = fits.table_to_hdu(cat_table) + cat_hdu.header["EXTNAME"] = "TOC" + + return cat_hdu + + +def make_trace_hdu(const_wave_coords_dicts: list[dict], + n_extra_points: int = 0, + spline_order: int = 1): + """ + Creates the BinTableHDU objects needed to make a spectral trace HDUList + + Parameters + ---------- + const_wave_coords_dicts : list of dicts + A list containing dictionaries for every line of constant wavelength + [{"wave": w [um], "x": list(xs) [mm],"y": list(ys) [mm]}, + {..}, ..] + n_extra_points : int, 2-tuple + Default = 0. How many extra points to put in the spectral table + If a 2-tuple is passed, then n extra points will be added such: + (wavelength, spatial) + spline_order : int + Default = 1 + + Examples + -------- + :: + dict_list = [{"wave": 1.9, "x": [0, 1], "y": [-2, -2]}, + {"wave": 2.4, "x": [0, 1], "y": [2, 2]}] + tg = TraceHDUGenerator(dict_list, n_extra_points=3) + tg.trace_hdu + + """ + wave, x, y = coords_from_lines_of_const_wavelength(const_wave_coords_dicts, + n_extra_points, + spline_order) + tbl = Table(names=["wavelength", "x", "y"], + units=[u.um, u.mm, u.mm], + data=[wave, x, y]) + + hdu = fits.table_to_hdu(tbl) + hdu.header[""] + return + + +def coords_from_lines_of_const_wavelength(const_wave_coords_dicts, + n_extra_points=0, spline_order=1): + """ + Returns expanded coordinates for constant wavelength in a spectral trace + + Parameters + ---------- + const_wave_coords_dicts : list of dicts + A list containing dictionaries for every line of constant wavelength + [{"wave": w [um], "x": list(xs) [mm],"y": list(ys) [mm]}, + {..}, ..] + n_extra_points : int, 2-tuple + Default = 0. How many extra points to put in the spectral table + If a 2-tuple is passed, then n extra points will be added such: + (wavelength, spatial) + spline_order : int + Default = 1 + + Examples + -------- + The following code creates a vertically bent trace with horizontal lines of + constant wavelength, with 20 additional points along the wavelength + direction, but 0 additional points along the slit direction + :: + dict_list = [{"wave": 1.9, "x": [0, 1, 2], "y": [-2.5, -2, -1.5]}, + {"wave": 2.15, "x": [1, 2, 3], "y": [0., 0, 0.]}, + {"wave": 2.4, "x": [0, 1, 2], "y": [1.5, 2, 2.5]}] + w, x, y = coords_from_lines_of_const_wavelength(dict_list, + n_extra_points=(20, 0), + spline_order=2) + + """ + # assert spline_order < len(const_wave_coords_dicts) + # assert spline_order < len(const_wave_coords_dicts[0]["x"]) + + x_im = np.array([dic["x"] for dic in const_wave_coords_dicts]) + y_im = np.array([dic["y"] for dic in const_wave_coords_dicts]) + w_im = np.array([[dic["wave"]]*len(dic["x"]) for dic in const_wave_coords_dicts]) + + cube = np.dstack([w_im, x_im, y_im]) + + xy_scale = 1 + np.array(n_extra_points) / np.array(x_im.shape) + if any(s != 1. for s in xy_scale): + cube = zoom(cube.astype(float), zoom=(*xy_scale, 1), + order=spline_order) + + ws, xs, ys = cube.reshape(-1, cube.shape[-1]).T + return ws, xs, ys diff --git a/scopesim/optics/image_plane_utils.py b/scopesim/optics/image_plane_utils.py index d0742856..a771623f 100644 --- a/scopesim/optics/image_plane_utils.py +++ b/scopesim/optics/image_plane_utils.py @@ -9,6 +9,7 @@ from astropy.io import fits from astropy.table import Table from astropy.utils.exceptions import AstropyWarning +from astropy.wcs import wcs from scipy import ndimage as ndi from ..utils import (unit_from_table, quantity_from_table, has_needed_keywords, @@ -345,6 +346,22 @@ def _add_intpixel_sources_to_canvas(canvas_hdu, xpix, ypix, flux, mask): def _add_subpixel_sources_to_canvas(canvas_hdu, xpix, ypix, flux, mask): + """ + + Parameters + ---------- + canvas_hdu : fits.ImageHDU + xpix, ypix : list of floats + [pix] Coordinates of sources in units of images pixel + flux : list of floats + mask : list of bools + flags for adding or ignoring a source + + Returns + ------- + canvas_hdu : fits.ImageHDU + + """ canvas_hdu.header["comment"] = f"Adding {len(flux)} sub-pixel files" canvas_shape = canvas_hdu.data.shape for xpx, ypx, flx, msk in zip(xpix, ypix, flux, mask): @@ -1029,6 +1046,51 @@ def split_header(hdr, chunk_size, wcs_suffix=""): return hdr_list +def extract_region_from_imagehdu(hdu, x_edges, y_edges, wcs_suffix=""): + """ + Returns an ImageHDU for a subsection of the + + Parameters + ---------- + hdu + x_edges : list of floats + [deg, mm] + y_edges : list of floats + [deg, mm] + + Returns + ------- + + """ + + s = "D" if wcs_suffix == "D" else " " + w = wcs.WCS(hdu.header, key=s, naxis=2) + xps, yps = w.all_world2pix([min(x_edges), max(x_edges)], + [min(y_edges), max(y_edges)], 1) + #xps, yps = np.round(xps).astype(int), np.round(yps).astype(int) + xps = [int(np.floor(min(xps))), int(np.floor(max(xps)))] + yps = [int(np.floor(min(yps))), int(np.floor(max(yps)))] + + w, h = hdu.header["NAXIS1"], hdu.header["NAXIS2"] + assert xps[0] >= 0 and xps[1] <= w and yps[0] >= 0 and yps[1] <= h, \ + f"Pixel edges ({xps, yps}) were not within the bounds " \ + f"of the array shape ({0, hdu.data.shape[-1], 0, hdu.data.shape[-2]})" + + if hdu.header["NAXIS"] == 2: + new_data = hdu.data[yps[0]:yps[1]+1, xps[0]:xps[1]+1] + elif hdu.header["NAXIS"] == 3: + new_data = hdu.data[:, yps[0]:yps[1]+1, xps[0]:xps[1]+1] + + new_hdu = fits.ImageHDU(data=new_data) + new_hdu.header.update(hdu.header) + new_hdu.header.update({"CRVAL1": np.average(x_edges), + "CRVAL2": np.average(y_edges), + "CRPIX1": np.average(xps), + "CRPIX2": np.average(yps)}) + + return new_hdu + + def _fix_360(arr): """Fix the "full circle overflow" that occurs with deg.""" if isinstance(arr, u.Quantity): diff --git a/scopesim/rc.py b/scopesim/rc.py index cf486656..3ea843e4 100644 --- a/scopesim/rc.py +++ b/scopesim/rc.py @@ -34,3 +34,9 @@ Path(__pkg_dir__).absolute(), *[Path(pth).absolute() for pth in __config__["!SIM.file.search_path"]], ]) + +# Needed by kl/mos_branch +# Can refactor this later to remove it once things work again +import os +__basic_inst_path__ = os.path.abspath(os.path.join(__pkg_dir__, + "tests/mocks/basic_instrument/")) diff --git a/scopesim/tests/conftest.py b/scopesim/tests/conftest.py index 0be2a072..05a659cc 100644 --- a/scopesim/tests/conftest.py +++ b/scopesim/tests/conftest.py @@ -91,6 +91,12 @@ def patch_mock_path_micado(mock_path_micado): yield +@pytest.fixture(scope="package") +def mock_path_basic_instrument(): + """Path to MICADO mock files.""" + return MOCK_DIR / "basic_instrument" + + @pytest.fixture(scope="function") def no_file_error(): """Patch currsys to avoid missing file error.""" diff --git a/scopesim/tests/mocks/basic_instrument/INS_mos_apertures.dat b/scopesim/tests/mocks/basic_instrument/INS_mos_apertures.dat new file mode 100644 index 00000000..ffe52115 --- /dev/null +++ b/scopesim/tests/mocks/basic_instrument/INS_mos_apertures.dat @@ -0,0 +1,11 @@ +# x_unit : arcsec +# y_unit : arcsec +# r_unit : arcsec +# angle_unit : deg + +id x y r angle conserve_image shape +0 -5 -5 0.4 0 False round +1 5 -5 0.6 0 False rect +2 0 0 0.8 0 False hex +3 -5 5 1 0 False oct +4 5 5 1.2 0 False 3 \ No newline at end of file diff --git a/scopesim/tests/mocks/basic_instrument/INS_mos_bundle_ter_corrections.dat b/scopesim/tests/mocks/basic_instrument/INS_mos_bundle_ter_corrections.dat new file mode 100644 index 00000000..b7a6769d --- /dev/null +++ b/scopesim/tests/mocks/basic_instrument/INS_mos_bundle_ter_corrections.dat @@ -0,0 +1,20 @@ +# source : MOSAIC MVP inst pkg +# zenith_dist_unit : deg +# wavelength_unit : um + +wavelength zenith_dist TRACE_ALL TRACE_0 TRACE_1 TRACE_2 TRACE_3 TRACE_4 TRACE_5 TRACE_6 +0.920 0. 0.590 0.165 0.068 0.076 0.069 0.069 0.075 0.068 +1.202 0. 0.638 0.213 0.067 0.078 0.069 0.068 0.076 0.066 +1.638 0. 0.688 0.290 0.062 0.076 0.064 0.063 0.073 0.060 +0.920 15 0.577 0.157 0.068 0.075 0.069 0.068 0.074 0.067 +1.202 15 0.626 0.203 0.067 0.078 0.069 0.068 0.076 0.066 +1.638 15 0.678 0.277 0.063 0.076 0.065 0.064 0.073 0.061 +0.920 30 0.537 0.135 0.065 0.071 0.066 0.066 0.070 0.064 +1.202 30 0.588 0.174 0.066 0.075 0.067 0.067 0.073 0.065 +1.638 30 0.645 0.239 0.064 0.076 0.066 0.065 0.073 0.062 +0.920 45 0.461 0.102 0.058 0.063 0.059 0.059 0.062 0.058 +1.202 45 0.515 0.131 0.062 0.068 0.063 0.063 0.067 0.061 +1.638 45 0.578 0.180 0.063 0.073 0.065 0.064 0.071 0.062 +0.920 60 0.334 0.063 0.045 0.046 0.045 0.045 0.046 0.045 +1.202 60 0.381 0.077 0.050 0.052 0.050 0.050 0.052 0.049 +1.638 60 0.443 0.104 0.055 0.060 0.056 0.055 0.059 0.054 \ No newline at end of file diff --git a/scopesim/tests/mocks/basic_instrument/INS_mos_traces.fits b/scopesim/tests/mocks/basic_instrument/INS_mos_traces.fits new file mode 100644 index 00000000..d0dff37e Binary files /dev/null and b/scopesim/tests/mocks/basic_instrument/INS_mos_traces.fits differ diff --git a/scopesim/tests/mocks/basic_instrument/YAML_basic_instrument.yaml b/scopesim/tests/mocks/basic_instrument/YAML_basic_instrument.yaml index aaa74f83..fb4557fe 100644 --- a/scopesim/tests/mocks/basic_instrument/YAML_basic_instrument.yaml +++ b/scopesim/tests/mocks/basic_instrument/YAML_basic_instrument.yaml @@ -51,3 +51,11 @@ effects : include: "!OBS.include_slicer" kwargs: filename: "INS_ifu_apertures.dat" + extend_fov_beyond_slit: 5 + +- name: fibre_list + description: collection of apertures corresponding to the MOS fibre heads + class: FibreApertureList + include: "!OBS.include_fibres" + kwargs: + filename: "INS_mos_apertures.dat" diff --git a/scopesim/tests/mocks/basic_instrument/YAML_mode_mos.yaml b/scopesim/tests/mocks/basic_instrument/YAML_mode_mos.yaml new file mode 100644 index 00000000..44a8791f --- /dev/null +++ b/scopesim/tests/mocks/basic_instrument/YAML_mode_mos.yaml @@ -0,0 +1,26 @@ +# MOS Spectroscopy mode specific Effects and properties +name : mode_mos +alias : INST + +properties : + decouple_detector_from_sky_headers: True + +effects : +- name : mos_spectral_traces + description : list of mos trace geometries on the focal plane + class : UnresolvedSpectralTraceList + kwargs : + filename : "INS_mos_traces.fits" + +--- + +name: config_overrides +alias: SIM +description: SIM config param overrides + +properties: + spectral: + wave_min: 1.7 + wave_mid: 2.1 + wave_max: 2.5 + spectral_bin_width: !!float 9E-4 \ No newline at end of file diff --git a/scopesim/tests/mocks/basic_instrument/code/make_spectral_trace.py b/scopesim/tests/mocks/basic_instrument/code/make_spectral_trace.py index 388c2124..bf088ebd 100644 --- a/scopesim/tests/mocks/basic_instrument/code/make_spectral_trace.py +++ b/scopesim/tests/mocks/basic_instrument/code/make_spectral_trace.py @@ -75,12 +75,13 @@ def make_lss_trace_file(): # make_lss_trace_file() + def make_ifu_trace_file(): names = [f"Ap{i}" for i in range(5)] wave_min = 1.75 wave_max = 2.5 - s_cens = np.array([-64., -32., 0., 32., 64.]) # across detector - dss = np.array([-14., -7., 0., 7., 14.]) # across slit + s_cens = np.array([-64., -32., 0., 32., 64.]) # across detector [arcsec] + dss = np.array([-14., -7., 0., 7., 14.]) # across slit [arcsec] hdus = [] for s_cen, name in zip(s_cens, names): @@ -118,3 +119,45 @@ def make_ifu_trace_file(): hdul.writeto("../INS_ifu_traces.fits", overwrite=True) # make_ifu_trace_file() + + +def make_mos_trace_file(): + names = [f"Ap{i}" for i in range(5)] + wave_min = 1.75 + wave_max = 2.5 + # detector edges (x,y) +/-5.12 [mm] + traces = [([-4, -4, -4], [-4, 0, 4]), # vertical line + ([-2.5, -2, -1.5], [-4, 0, 4]), # tilted line + ([0, 0.5, 0.7, 0.5, 0], [-4, -2, -0, 2, 4]), # curve + ([2, 2.5, 2, 1.5, 2], [-4, -2, -0, 2, 4]), # S-shape + ([3, 4, 5], [-0.1, 0, 0.1])] # horizontal line + + hdus = [] + for trace, name in zip(traces, names): + xs, ys = trace + waves = np.geomspace(wave_min, wave_max, len(xs)) + + tbl = Table(names=["wavelength", "x", "y"], + data=[waves, xs, ys], + units=[u.um, u.mm, u.mm]) + hdu = fits.table_to_hdu(tbl) + hdu.header["EXTNAME"] = f"TRACE_{name}" + hdus += [hdu] + + tbl = Table( + names=["description", "extension_id", "aperture_id", "image_plane_id"], + data=[[hdu.header["EXTNAME"] for hdu in hdus], + [2, 3, 4, 5, 6], + [0, 1, 2, 3, 4], + [0, 0, 0, 0, 0]]) + cat_hdu = fits.table_to_hdu(tbl) + cat_hdu.header["EXTNAME"] = "TOC" + + pri_hdr = fits.PrimaryHDU() + pri_hdr.header["ECAT"] = 1 # where to find the catalogue table + pri_hdr.header["EDATA"] = 2 # where the data tables start + + hdul = fits.HDUList([pri_hdr, cat_hdu] + hdus) + hdul.writeto("../INS_mos_traces.fits", overwrite=True) + +# make_mos_trace_file() diff --git a/scopesim/tests/mocks/load_basic_instrument.py b/scopesim/tests/mocks/load_basic_instrument.py index f750645a..d0487554 100644 --- a/scopesim/tests/mocks/load_basic_instrument.py +++ b/scopesim/tests/mocks/load_basic_instrument.py @@ -8,6 +8,8 @@ def load_example_optical_train(**kwargs): """ Return an basic example ``OpticalTrain`` object with IMG and SPEC modes. + Available modes: ["imaging", "spectroscopy", "ifu"] + Parameters ---------- Any of the additional parameters taken by UserCommands. diff --git a/scopesim/tests/mocks/py_objects/header_objects.py b/scopesim/tests/mocks/py_objects/header_objects.py index 748998a1..cb7f4260 100644 --- a/scopesim/tests/mocks/py_objects/header_objects.py +++ b/scopesim/tests/mocks/py_objects/header_objects.py @@ -5,19 +5,18 @@ from scopesim.optics.image_plane_utils import header_from_list_of_xy -def _basic_fov_header(): - w, h = 150, 150 +def _basic_fov_header(w=150, h=150, cdelt=0.1, cdelt_d=1): skywcs = wcs.WCS(naxis=2) # skywcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] skywcs.wcs.ctype = ["LINEAR", "LINEAR"] - skywcs.wcs.cdelt = [0.1, 0.1] + skywcs.wcs.cdelt = [cdelt, cdelt] skywcs.wcs.cunit = ["arcsec", "arcsec"] skywcs.wcs.crval = [0, 0] skywcs.wcs.crpix = [(w + 1) / 2, (h + 1) / 2] detwcs = wcs.WCS(naxis=2, key="D") detwcs.wcs.ctype = ["LINEAR", "LINEAR"] - detwcs.wcs.cdelt = [1, 1] + detwcs.wcs.cdelt = [cdelt_d, cdelt_d] detwcs.wcs.cunit = ["mm", "mm"] detwcs.wcs.crval = [0, 0] detwcs.wcs.crpix = [(w + 1) / 2, (h + 1) / 2] diff --git a/scopesim/tests/tests_effects/test_ApertureList.py b/scopesim/tests/tests_effects/test_ApertureList.py index 8a077e09..3a6adc16 100644 --- a/scopesim/tests/tests_effects/test_ApertureList.py +++ b/scopesim/tests/tests_effects/test_ApertureList.py @@ -1,8 +1,9 @@ - import numpy as np from astropy.io import fits +from scopesim import rc from scopesim.effects import ApertureList, ApertureMask +from scopesim.optics.fov_manager import FovVolumeList import matplotlib.pyplot as plt @@ -10,18 +11,28 @@ PLOTS = False +def basic_aperture_list(**kwargs): + params = {"array_dict": {"id": [0, 1, 2], + "left": [-3, -1, -1], + "right": [3, 1, 1], + "top": [-1, 1, 3], + "bottom": [-3, -1, 1], + "angle": [0, 0, 0], + "conserve_image": [True]*3, + "shape": ["rect", "hex", 7]}, + "x_unit": "arcsec", + "y_unit": "arcsec", + "pixel_scale": 0.01} + params.update(kwargs) + return ApertureList(**params) + + class TestInit: def test_initialises_with_nothing(self): isinstance(ApertureList(), ApertureList) def test_initialises_with_array_dict(self): - params = {"array_dict": {"id": [0], "left": [-1], "right": [1], - "top": [1], "bottom": [1], "angle": [0], - "conserve_image": [True], "shape": [7]}, - "x_unit": "arcsec", - "y_unit": "arcsec", - "pixel_scale": 0.01} - apl = ApertureList(**params) + apl = basic_aperture_list() assert isinstance(apl, ApertureList) def test_initialises_with_filename(self, mock_path): @@ -29,6 +40,42 @@ def test_initialises_with_filename(self, mock_path): assert isinstance(apl, ApertureList) +class TestApplyTo: + def test_extracts_volumes_for_all_apertures(self): + apl = basic_aperture_list() + fvl = FovVolumeList() + fvl = apl.apply_to(fvl) + assert len(fvl.volumes) == 3 + assert fvl.volumes[0]["x_min"] == -3 + assert fvl.volumes[1]["x_min"] == -1 + + def test_extracts_single_volume_for_all_apertures(self): + apl = basic_aperture_list(fov_for_each_aperture=False) + fvl = FovVolumeList() + fvl = apl.apply_to(fvl) + assert len(fvl.volumes) == 1 + assert fvl.volumes[0]["x_min"] == -3 + assert fvl.volumes[0]["y_max"] == 3 + + def test_increases_all_volume_sizes_for_extend_beyond_slit(self): + apl = basic_aperture_list(fov_for_each_aperture=True, + extend_fov_beyond_slit=2) + fvl = FovVolumeList() + fvl = apl.apply_to(fvl) + assert len(fvl.volumes) == 3 + assert fvl.volumes[0]["x_min"] == -5 + assert fvl.volumes[1]["x_min"] == -3 + + def test_increases_single_volume_for_extend_beyond_slit(self): + apl = basic_aperture_list(fov_for_each_aperture=False, + extend_fov_beyond_slit=3) + fvl = FovVolumeList() + fvl = apl.apply_to(fvl) + assert len(fvl.volumes) == 1 + assert fvl.volumes[0]["x_min"] == -6 + assert fvl.volumes[0]["y_max"] == 6 + + class TestApertures: def test_returns_list_of_aperture_masks(self, mock_path): apl = ApertureList(filename=str(mock_path / "test_aperture_list.dat"), @@ -41,30 +88,3 @@ def test_returns_list_of_aperture_masks(self, mock_path): plt.subplot(2, 2, ii+1) plt.imshow(apertures[ii].mask.T) plt.show() - - -class TestFovGrid: - def test_returns_headers(self, mock_path): - apl = ApertureList(filename=str(mock_path / "test_aperture_list.dat"), - no_mask=False, pixel_scale=0.01) - hdrs = apl.fov_grid() - assert all([isinstance(hdr, fits.Header) for hdr in hdrs]) - - if PLOTS: - from scopesim.optics.image_plane import ImagePlane - from scopesim.optics.image_plane_utils import header_from_list_of_xy - from astropy.io.fits import ImageHDU - - x = np.array([-3, 3]) / 3600. - y = np.array([-3, 3]) / 3600. - pixel_scale = 0.01 / 3600 - hdr = header_from_list_of_xy(x, y, pixel_scale) - implane = ImagePlane(hdr) - - for i in range(4): - hdu = ImageHDU(data=apl[i].mask.astype(int)+1, - header=apl[i].header) - implane.add(hdu) - - plt.imshow(implane.data, origin="lower") - plt.show() diff --git a/scopesim/tests/tests_effects/test_ApertureMask.py b/scopesim/tests/tests_effects/test_ApertureMask.py index 8c282b39..be4582f5 100644 --- a/scopesim/tests/tests_effects/test_ApertureMask.py +++ b/scopesim/tests/tests_effects/test_ApertureMask.py @@ -1,12 +1,13 @@ - import pytest from pytest import approx import numpy as np from astropy.io import fits +from scopesim import rc from scopesim.effects import ApertureMask from scopesim.effects.apertures import points_on_a_circle, make_aperture_polygon +from scopesim.optics.fov_manager import FovVolumeList import matplotlib.pyplot as plt @@ -14,16 +15,20 @@ PLOTS = False +def basic_aperture_mask(x=(-2, -1, 1, 2), y=(-1, -2, 2, 1), **kwargs): + params = {"array_dict": {"x": list(x), "y": list(y)}, + "x_unit": "arcsec", + "y_unit": "arcsec"} + params.update(kwargs) + return ApertureMask(**params) + + class TestInit: def test_initialises_with_nothing(self): assert isinstance(ApertureMask(), ApertureMask) def test_initialises_with_x_y_array_dict(self): - kwargs = {"array_dict": {"x": [-2, -1, 1, 2], - "y": [-1, -2, 2, 1]}, - "x_unit": "arcsec", - "y_unit": "arcsec"} - apm = ApertureMask(**kwargs) + apm = basic_aperture_mask() assert isinstance(apm, ApertureMask) assert "x" in apm.table.colnames @@ -33,15 +38,28 @@ def test_initialises_from_file(self, mock_path): assert "y" in apm.table.colnames +class TestApplyTo: + def test_shrinks_fov_volume_list_as_expected(self): + apm = basic_aperture_mask(x=[-2, 2, 2, -2], y=[-1, -1, 1, 1], + pixel_scale=0.1) + fvl = FovVolumeList() + fvl = apm.apply_to(fvl) + assert fvl.volumes[0]["x_min"] == -2 + + def test_shrinks_fov_volume_list_with_extended_region(self): + apm = basic_aperture_mask(x=[-2, 2, 2, -2], y=[-1, -1, 1, 1], + extend_fov_beyond_slit=2, pixel_scale=0.1) + fvl = FovVolumeList() + fvl = apm.apply_to(fvl) + assert fvl.volumes[0]["x_min"] == -4 + + class TestHeader: def test_returns_header_surrounding_tilted_aperture(self): e = 1e-7 - kwargs = {"array_dict": {"x": [-2, -1, 1+e, 2+e], - "y": [-1, -2, 2+e, 1+e]}, - "x_unit": "arcsec", - "y_unit": "arcsec", - "pixel_scale": 0.1} - apm = ApertureMask(**kwargs) + apm = basic_aperture_mask(x=[-2, -1, 1+e, 2+e], + y=[-1, -2, 2+e, 1+e], + pixel_scale=0.1) hdr = apm.header assert isinstance(hdr, fits.Header) assert hdr["NAXIS1"] == hdr["NAXIS2"] == 40 @@ -77,11 +95,6 @@ def test_returns_mask_for_super_wierd_aperture(self): plt.show() -class TestFovGrid: - # not needed because all it does is return the header - pass - - class TestMakeAperturePolygon: @pytest.mark.parametrize("shape, n_corners", [("rect", 4), ("hex", 6), ("round", 32), (7, 7)]) diff --git a/scopesim/tests/tests_effects/test_MosaicBundleTracePSF.py b/scopesim/tests/tests_effects/test_MosaicBundleTracePSF.py new file mode 100644 index 00000000..3827ff2e --- /dev/null +++ b/scopesim/tests/tests_effects/test_MosaicBundleTracePSF.py @@ -0,0 +1,35 @@ +import pytest +from pytest import approx + +from matplotlib import pyplot as plt +import numpy as np + +from scopesim.effects import MosaicBundleTracePSF, AnalyticalPSF + +PLOTS = False + + +class TestInit: + def test_initialised_with_nothing(self): + assert isinstance(MosaicBundleTracePSF(), AnalyticalPSF) + + def test_initalises_with_radial_profile(self): + mbt_psf = MosaicBundleTracePSF(r=[0, 3.5, 4.2], + z=[1, 1, 0]) + assert isinstance(mbt_psf, MosaicBundleTracePSF) + + +class TestGetKernel: + def test_returned_kernel_is_zipped_version_of_single_kernel(self): + mbt_psf = MosaicBundleTracePSF(n_traces_per_bundle=2, + trace_spacing=2.1, + trace_flux_weights=[0.5, 2, 1]) + kernel = mbt_psf.get_kernel(None) + + if PLOTS: + plt.imshow(kernel) + plt.show() + + ws = mbt_psf.meta["trace_flux_weights"] + assert np.sum(kernel) == np.sum(ws) * np.sum(mbt_psf.single_kernel) + assert kernel.shape[0] > ws * kernel.shape[0] diff --git a/scopesim/tests/tests_effects/test_RadialProfilePSF.py b/scopesim/tests/tests_effects/test_RadialProfilePSF.py new file mode 100644 index 00000000..96c08524 --- /dev/null +++ b/scopesim/tests/tests_effects/test_RadialProfilePSF.py @@ -0,0 +1,56 @@ +import pytest +from pytest import approx + +from matplotlib import pyplot as plt +import numpy as np + +from scopesim.effects import RadialProfilePSF + +PLOTS = False + + +class TestInit: + def test_throws_error_for_no_profile(self): + with pytest.raises(ValueError): + RadialProfilePSF() + + def test_initalises_with_radial_profile(self): + rppsf = RadialProfilePSF(r=[0, 3.5, 4.2], + z=[1, 1, 0]) + assert isinstance(rppsf, RadialProfilePSF) + +class TestApplyTo: + def test_kernel_sums_to_unity(self): + rppsf = RadialProfilePSF(r=[0, 3.5, 4.2], + z=[1, 1, 0]) + kernel = rppsf.get_kernel(None) + assert np.sum(kernel) == approx(1) + + def test_kernel_size_for_unit_pixel(self): + rppsf = RadialProfilePSF(r=[0, 3.5, 4.2], + z=[1, 1, 0]) + kernel = rppsf.get_kernel(None) + assert kernel.shape[0] == 11 and kernel.shape[1] == 11 + + def test_kernel_size_for_unit_mm(self): + rppsf = RadialProfilePSF(r=[0, 3.5, 4.2], + z=[1, 1, 0], + unit="mm", + pixel_scale=0.25, + plate_scale=1) + kernel = rppsf.get_kernel(None) + assert kernel.shape[0] == int(np.ceil(4.2/0.25)*2+1) # 35 + + if PLOTS: + plt.imshow(kernel) + plt.show() + + def test_really_small_kernel_size(self): + rppsf = RadialProfilePSF(r=[0, 0.9, 1.28, 1.5], + z=[1, 1, 0.5, 0]) + kernel = rppsf.get_kernel(None) + assert kernel.shape[0] == 5 + + if PLOTS: + plt.imshow(kernel) + plt.show() diff --git a/scopesim/tests/tests_effects/test_SpectralTraceList.py b/scopesim/tests/tests_effects/test_SpectralTraceList.py index fbc460ab..05074d80 100644 --- a/scopesim/tests/tests_effects/test_SpectralTraceList.py +++ b/scopesim/tests/tests_effects/test_SpectralTraceList.py @@ -1,17 +1,34 @@ """Tests for module spectral_trace_list.py""" - +import os import pytest from unittest.mock import patch +import numpy as np from astropy.io import fits +from astropy.wcs import WCS +from astropy import units as u +from matplotlib import pyplot as plt +from matplotlib.colors import LogNorm +from synphot import Empirical1D, SourceSpectrum +from scopesim.effects import spectral_trace_list as spt +from scopesim.optics.fov_manager import FovVolumeList from scopesim.effects.spectral_trace_list import SpectralTraceList, \ SpectralTraceListWheel from scopesim.effects.spectral_trace_list_utils import SpectralTrace from scopesim.tests.mocks.py_objects import trace_list_objects as tlo from scopesim.tests.mocks.py_objects import header_objects as ho +from scopesim.base_classes import PoorMansHeader +from scopesim import rc + +MOCK_PATH = os.path.abspath(os.path.join(os.path.dirname(__file__), + "../mocks/MICADO_SPEC/")) +if MOCK_PATH not in rc.__search_path__: + rc.__search_path__ += [MOCK_PATH] +if rc.__basic_inst_path__ not in rc.__search_path__: + rc.__search_path__.append(rc.__basic_inst_path__) PLOTS = False @@ -39,6 +56,7 @@ class TestInit: def test_initialises_with_nothing(self): assert isinstance(SpectralTraceList(), SpectralTraceList) + # @pytest.mark.usefixtures("full_trace_list") def test_initialises_with_a_hdulist(self, full_trace_list): spt = SpectralTraceList(hdulist=full_trace_list) assert isinstance(spt, SpectralTraceList) @@ -67,11 +85,120 @@ def test_setitem_appends_correctly(self, full_trace_list): assert len(slist.spectral_traces) == n_trace + 1 +@pytest.mark.skip(reason="Ignoring old Spectroscopy integration tests") +class TestGetFOVHeaders: + @pytest.mark.usefixtures("full_trace_list", "slit_header") + def test_gets_the_headers(self, full_trace_list, slit_header): + sptl = spt.SpectralTraceList(hdulist=full_trace_list) + params = {"pixel_scale": 0.015, "plate_scale": 0.26666, + "wave_min": 0.7, "wave_max": 2.5} + hdrs = sptl.get_fov_headers(slit_header, **params) + + # assert all([isinstance(hdr, fits.Header) for hdr in hdrs]) + assert all([isinstance(hdr, PoorMansHeader) for hdr in hdrs]) + # ..todo:: add in some better test of correctness + + if PLOTS: + # pixel coords + for hdr in hdrs[::50]: + xp = [0, hdr["NAXIS1"], hdr["NAXIS1"], 0] + yp = [0, 0, hdr["NAXIS2"], hdr["NAXIS2"]] + wcs = WCS(hdr, key="D") + # world coords + xw, yw = wcs.all_pix2world(xp, yp, 1) + plt.fill(xw / hdr["CDELT1D"], yw / hdr["CDELT2D"], alpha=0.2) + plt.show() + + def test_gets_headers_from_real_file(self): + slit_hdr = ho._long_micado_slit_header() + # slit_hdr = ho._short_micado_slit_header() + wave_min = 1.0 + wave_max = 1.3 + sptl = spt.SpectralTraceList(filename="TRACE_15arcsec.fits", + s_colname="xi", + wave_colname="lam", + spline_order=1) + params = {"wave_min": wave_min, "wave_max": wave_max, + "pixel_scale": 0.004, "plate_scale": 0.266666667} + hdrs = sptl.get_fov_headers(slit_hdr, **params) + assert isinstance(sptl, spt.SpectralTraceList) + + print(len(hdrs)) + + if PLOTS: + sptl.plot(wave_min, wave_max) + + # pixel coords + for hdr in hdrs[::300]: + xp = [0, hdr["NAXIS1"], hdr["NAXIS1"], 0] + yp = [0, 0, hdr["NAXIS2"], hdr["NAXIS2"]] + wcs = WCS(hdr, key="D") + # world coords + xw, yw = wcs.all_pix2world(xp, yp, 1) + plt.plot(xw, yw, alpha=0.2) + plt.show() + + +# class TestApplyTo: +# def test_fov_setup_base_returns_only_extracted_fov_limits(self): +# fname = r"F:\Work\irdb\MICADO\TRACE_MICADO.fits" +# sptl = spt.SpectralTraceList(filename=fname, s_colname='xi') +# +# fvl = FovVolumeList() +# fvl = spt.apply_to(fvl) +# +# assert len(fvl) == 17 + + @pytest.fixture(name="spectral_trace_list", scope="class") def fixture_spectral_trace_list(): """Instantiate a SpectralTraceList""" return SpectralTraceList(hdulist=tlo.make_trace_hdulist()) + +def test_set_pc_matrix(rotation_ang=0, shear_ang=10): + n = 100 + im = np.arange(n**2).reshape(n, n) + hdu = fits.ImageHDU(im) + hdr_dict = {"CTYPE1": "LINEAR", + "CTYPE2": "LINEAR", + "CUNIT1": "deg", + "CUNIT2": "deg", + "CDELT1": 1, + "CDELT2": 1, + "CRVAL1": 0, + "CRVAL2": 0, + "CRPIX1": 0, + "CRPIX2": 0} + hdu.header.update(hdr_dict) + + c = np.cos(rotation_ang / 57.29578) * 2 + s = np.sin(rotation_ang / 57.29578) * 2 + t = np.tan(shear_ang / 57.29578) + + n = 5 + pc_dict = {"PC1_1": c + t*s, + "PC1_2": -s + t*c, + "PC2_1": s, + "PC2_2": c} + det = np.sqrt(np.abs(pc_dict["PC1_1"] * pc_dict["PC2_2"] - \ + pc_dict["PC1_2"] * pc_dict["PC2_1"])) + for key in pc_dict: + pc_dict[key] /= det + hdu.header.update(pc_dict) + w = WCS(hdu) + + xd = np.array([0, 10, 10, 0]) + yd = np.array([0, 0, 10, 10]) + xs, ys = w.all_pix2world(xd, yd, 1) + + if PLOTS: + plt.figure(figsize=(6, 6)) + plt.plot(xd, yd, "o-") + plt.plot(xs, ys, "o-") + plt.show() + + class TestRectification: def test_rectify_cube_not_implemented(self, spectral_trace_list): hdulist = fits.HDUList() @@ -101,3 +228,70 @@ def test_basic_init(self): assert stw.meta["trace_list_names"] == ["foo"] assert isinstance(stw.trace_lists["foo"], SpectralTraceList) assert stw.trace_lists["foo"].meta["filename"] == "bogus_foo" + + +class TestUnresolvedSpectralTraceListInit: + def test_init(self): + sptl = spt.UnresolvedSpectralTraceList(filename="INS_mos_traces.fits") + assert len(sptl.spectral_traces) == 5 + assert sptl.spectral_traces["TRACE_Ap0"].meta["trace_id"] == "TRACE_Ap0" + + if PLOTS: + for trace in list(sptl.spectral_traces.values()): + plt.plot(trace.table["x"], trace.table["y"]) + plt.show() + + def test_get_xyz_for_trace(self): + sptl = spt.UnresolvedSpectralTraceList(filename="INS_mos_traces.fits", + plate_scale=20, + pixel_scale=0.2) + waves = np.linspace(1.75, 2.5, 101) * u.um + flux = np.zeros_like(waves.value) + flux[::10] = 1 + spec = SourceSpectrum(Empirical1D, points=waves, lookup_table=flux) + + for trace in sptl.spectral_traces.values(): + x, y, z = sptl.get_xyz_for_trace(trace, spec) + plt.plot(x, y) + + if PLOTS: + plt.show() + plt.pause(0) + + def test_trace_to_image(self): + sptl = spt.UnresolvedSpectralTraceList(filename="INS_mos_traces.fits", + plate_scale=20, + pixel_scale=0.2) + waves = np.linspace(1.75, 2.5, 101) * u.um + flux = np.ones_like(waves.value) + flux[::10] = 10 + spec = SourceSpectrum(Empirical1D, points=waves, lookup_table=flux) + + fig, axes = plt.subplots(2, 3) + + for ax, trace in zip(axes.flatten(), sptl.spectral_traces.values()): + x, y, z = sptl.get_xyz_for_trace(trace, spec) + image = sptl.trace_to_image(x, y, z) + ax.imshow(image, origin="lower") + + if PLOTS: + fig.show() + + +class TestMosaicSpectralTraceList(): + def test_intit(self): + sptl = spt.MosaicSpectralTraceList( + plate_scale=6.33333, + pixel_scale=0.095, + wave_min=1.420, + wave_mid=1.6385, + wave_max=1.857, + spectral_bin_width=0.0001 + ) + + if PLOTS: + sptl.plot("proj") + plt.show() + + # print(sptl.spectral_traces["TRACE_0"].table) + assert len(sptl.spectral_traces) == 14 diff --git a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py index 21dadc9a..5a27c222 100644 --- a/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py +++ b/scopesim/tests/tests_effects/test_SpectralTraceListUtils.py @@ -6,14 +6,21 @@ import pytest import numpy as np - from astropy.io import fits +from astropy.table import Table +import matplotlib.pyplot as plt from scopesim.effects.spectral_trace_list_utils import SpectralTrace from scopesim.effects.spectral_trace_list_utils import Transform2D, power_vector from scopesim.effects.spectral_trace_list_utils import make_image_interpolations +from scopesim.effects.spectral_trace_list_utils import coords_from_lines_of_const_wavelength as const_wave +from scopesim.effects.spectral_trace_list_utils import TracesHDUListGenerator +from scopesim.effects.spectral_trace_list_utils import make_trace_hdu from scopesim.tests.mocks.py_objects import trace_list_objects as tlo + +PLOTS = False + class TestSpectralTrace: """Tests not covered in test_SpectralTraceList.py""" def test_initialises_with_table(self): @@ -153,3 +160,100 @@ def test_interpolation_is_accurate(self): interps = make_image_interpolations(hdul, kx=1, ky=1) imginterp = interps[0](yy, xx, grid=False) assert np.allclose(imginterp, img) + + +class TestTracesHDUListGenerator: + def test_returns_hdulist_with_4_ext_for_2_traces(self): + dict_list = [{"wave": 1.9, "x": [0, 1], "y": [-2, -2]}, + {"wave": 2.4, "x": [0, 1], "y": [2, 2]}] + trace_hdu = make_trace_hdu(dict_list, n_extra_points=3) + sptl = TracesHDUListGenerator([trace_hdu, trace_hdu], + [0, 0], [0, 0]) + + assert isinstance(sptl.hdulist, fits.HDUList) + assert len(sptl.hdulist) == 4 + + +class TestMakeTraceHDU: + def test_returns_tablehdu_with_straight_trace(self): + dict_list = [{"wave": 1.9, "x": [0, 1], "y": [-2, -2]}, + {"wave": 2.4, "x": [0, 1], "y": [2, 2]}] + + trace_hdu = make_trace_hdu(dict_list, n_extra_points=3) + assert isinstance(trace_hdu, fits.BinTableHDU) + assert len(Table(trace_hdu.data)) == (2 + 3) ** 2 + + +class TestCoordsFromLinesOfConstWavelength: + def test_returns_input_for_no_expansion(self): + dict_list = [{"wave": 1.9, "x": [0, 1], "y": [-2, -2]}, + {"wave": 2.4, "x": [0, 1], "y": [2, 2]}] + w, x, y = const_wave(dict_list) + + assert min(w) == 1.9 and max(w) == 2.4 + assert min(x) == 0. and max(x) == 1. + assert min(y) == -2 and max(y) == 2. + + def test_returns_same_edges_when_expanded(self): + dict_list = [{"wave": 1.9, "x": [0, 1], "y": [-2, -2]}, + {"wave": 2.4, "x": [0, 1], "y": [2, 2]}] + w, x, y = const_wave(dict_list, n_extra_points=2) + + assert min(w) == 1.9 and max(w) == 2.4 + assert min(x) == 0. and max(x) == 1. + assert min(y) == -2 and max(y) == 2. + + assert len(x) == 4 + + def test_linear_interpolation_for_spline_1(self): + dict_list = [{"wave": 1.9, "x": [0, 1], "y": [-2, -2]}, + {"wave": 2.15, "x": [1, 2], "y": [0, 0]}, + {"wave": 2.4, "x": [0, 1], "y": [2, 2]}] + w, x, y = const_wave(dict_list, n_extra_points=2) + + if PLOTS: + plt.plot(x, y, "o") + plt.show() + + assert x[1] == 0.5 and x[2] == 1 + + def test_cubic_interpolation_for_spline_2(self): + dict_list = [{"wave": 1.9, "x": [0, 1, 2], "y": [-2, -2, -2]}, + {"wave": 2.15, "x": [1, 2, 3], "y": [0, 0, 0]}, + {"wave": 2.4, "x": [0, 1, 2], "y": [2, 2, 2]}] + w, x, y = const_wave(dict_list, n_extra_points=20, spline_order=2) + + if not PLOTS: + plt.scatter(x, y, c=w) + plt.show() + + def test_echelle_type_expansion(self): + dict_list = [{"wave": 1.9, "x": [0, 1, 2], "y": [-2.5, -2, -1.5]}, + {"wave": 2.15, "x": [1, 2, 3], "y": [0., 0, 0.]}, + {"wave": 2.4, "x": [0, 1, 2], "y": [1.5, 2, 2.5]}] + w, x, y = const_wave(dict_list, n_extra_points=20, spline_order=2) + + if PLOTS: + plt.scatter(x, y, c=w) + plt.show() + + def test_super_funky_expansions(self): + dict_list = [{"wave": 2.0, "x": [0, 1, 2], "y": [0, 0, 0]}, + {"wave": 2.1, "x": [0.5, 1, 1.5], "y": [2., 2, 2.]}, + {"wave": 2.2, "x": [0, 1, 2], "y": [4, 4, 4]}, + {"wave": 2.3, "x": [0.5, 1.5, 2.5], "y": [5.5, 6, 5.5]}, + {"wave": 2.4, "x": [0, 1, 2], "y": [8, 8, 8]}] + w, x, y = const_wave(dict_list, n_extra_points=(100, 10), spline_order=3) + + if PLOTS: + plt.scatter(x, y, c=w) + plt.show() + + def test_expands_only_wave_direction_with_tuple_extra_points(self): + dict_list = [{"wave": 1.9, "x": [0, 1], "y": [-2, -2]}, + {"wave": 2.4, "x": [0, 1], "y": [2, 2]}] + w, x, y = const_wave(dict_list, n_extra_points=(3, 0)) + + assert len(w) == 5 * 2 + assert len(set(w)) == 5 + assert len(set(x)) == 2 diff --git a/scopesim/tests/tests_effects/test_TERCurve.py b/scopesim/tests/tests_effects/test_TERCurve.py index 8e9a1847..dc3f7b62 100644 --- a/scopesim/tests/tests_effects/test_TERCurve.py +++ b/scopesim/tests/tests_effects/test_TERCurve.py @@ -1,12 +1,19 @@ import pytest - import numpy as np from matplotlib import pyplot as plt + from astropy import units as u +from astropy.wcs import wcs +from astropy.table import Table + +from synphot import Empirical1D, SourceSpectrum +from synphot.units import PHOTLAM +from scopesim.optics import FieldOfView from scopesim.effects import ter_curves as tc from scopesim.tests.mocks.py_objects import source_objects as so from scopesim.tests.mocks.py_objects import effects_objects as eo +from scopesim.tests.mocks.py_objects import header_objects as ho PLOTS = False @@ -20,6 +27,37 @@ def _filter_wheel(mock_path_micado): "filename_format": fname, "current_filter": "Br-gamma"}) +@pytest.fixture(scope="function") +def mos_bundle_tc(mock_path_basic_instrument): + fname = str(mock_path_basic_instrument / "INS_mos_bundle_ter_corrections.dat") + return tc.AirmassDependantFibreBundleTransmission(filename=fname, + zenith_dist=15, + n_fibres=7) + + +def mock_fov(trace_id=0, aperture_id=0, src_type="both"): + skyhdr = ho._basic_fov_header(w=7, h=7, cdelt=0.1) + waverange = (0.92, 1.638) * u.um + fov = FieldOfView(header=skyhdr, + waverange=waverange, + area=1*u.m**2, + sub_pixel=False, + conserve_image=False, + aperture_id=aperture_id, + trace_id=trace_id) + + fields = [so._unity_source(), + Table(names=["x", "y", "ref", "weight"], + data=[[0], [0], [0], [1]]) + ] + if "point" in src_type: fields = fields[1:] + elif "ext" in src_type: fields = fields[0:1] + fov.fields += fields + fov.spectra[0] = SourceSpectrum(Empirical1D, points=waverange, + lookup_table=[1, 1]*PHOTLAM) + + return fov + class TestTERCurveApplyTo: def test_adds_bg_to_source_if_source_has_no_bg(self): @@ -202,3 +240,48 @@ def test_initialises_with_correct_input(self): assert filt_wheel.filters["J"].throughput(1.15*u.um) == 0.9 assert filt_wheel.filters["J"].throughput(1.13*u.um) == 0. assert filt_wheel.meta["current_filter"] == "K" + + +class TestAirmassDependantFibreBundleTransmission: + def test_init_with_real_input(self, mos_bundle_tc): + assert isinstance(mos_bundle_tc, + tc.AirmassDependantFibreBundleTransmission) + assert mos_bundle_tc.ter_cube.shape == (5, 3) + + def test_get_table_for_zenith_dist(self, mos_bundle_tc): + tbl15 = mos_bundle_tc.get_table_for_zenith_dist(15) + tbl20 = mos_bundle_tc.get_table_for_zenith_dist(20) + + assert tbl15["TRACE_ALL"][0] > tbl20["TRACE_ALL"][0] + assert tbl15["TRACE_1"][1] > tbl20["TRACE_1"][1] + + def test_apply_to_correction_factor_on_single_bundle(self, mos_bundle_tc): + wave = 0.92 * u.um + fovs = [mock_fov(trace_id=i) for i in range(7)] + flux_0 = [fov.spectra[0](wave) for fov in fovs] + for fov in fovs: + mos_bundle_tc.apply_to(fov) + flux_1 = [fov.spectra[0](wave) for fov in fovs] + + # print(f"{flux_0=}") + # print(f"{flux_1=}") + assert all([f0 > f1 for f0, f1 in zip(flux_0, flux_1)]) + + def test_apply_to_correction_factor_on_three_bundle(self, mos_bundle_tc): + wave = 1.638 * u.um + fovs = [mock_fov(trace_id=i) for i in range(7 * 3)] + for fov in fovs: + mos_bundle_tc.apply_to(fov) + flux = [fov.spectra[0](wave) for fov in fovs] + + assert all([flux[i] == flux[i+7] == flux[i+7*2] for i in range(7)]) + + def test_apply_to_extended_source(self, mos_bundle_tc): + wave = 1.638 * u.um + fovs = [mock_fov(trace_id=i, src_type="ext") for i in range(7 * 3)] + flux_0 = [fov.spectra[0](wave) for fov in fovs] + for fov in fovs: + mos_bundle_tc.apply_to(fov) + flux_1 = [fov.spectra[0](wave) for fov in fovs] + + assert all([f0 == f1 for f0, f1 in zip(flux_0, flux_1)]) diff --git a/scopesim/tests/tests_optics/test_FieldOfView.py b/scopesim/tests/tests_optics/test_FieldOfView.py index 380325ab..bba88069 100644 --- a/scopesim/tests/tests_optics/test_FieldOfView.py +++ b/scopesim/tests/tests_optics/test_FieldOfView.py @@ -488,7 +488,7 @@ def test_makes_spectrum_from_all_types_of_source_object(self): if PLOTS: waves = fov.waveset - plt.plot(waves, spec(waves)) + plt.plot(waves, spec(waves), "r--") for spectrum in src_all.spectra: plt.plot(waves, spectrum(waves)) plt.show() diff --git a/scopesim/tests/tests_optics/test_FovVolumeList.py b/scopesim/tests/tests_optics/test_FovVolumeList.py index a4534a8b..4d94f907 100644 --- a/scopesim/tests/tests_optics/test_FovVolumeList.py +++ b/scopesim/tests/tests_optics/test_FovVolumeList.py @@ -19,6 +19,7 @@ def test_init_with_nothing(self): assert isinstance(fvl, FovVolumeList) assert len(fvl.volumes) == 1 + class TestSplit: @pytest.mark.parametrize("axis, value", [("wave", 3.0), ("x", 0.), diff --git a/setup.py b/setup.py new file mode 100644 index 00000000..e69de29b