Oc/grating efficiencies

scopesim/effects/__init__.py

@@ -6,6 +6,7 @@
 
 from .obs_strategies import *
 from .spectral_trace_list import *
+from .spectral_efficiency import *
 from .metis_lms_trace_list import *
 from .surface_list import *
 from .ter_curves import *

scopesim/effects/spectral_efficiency.py

@@ -0,0 +1,127 @@
+"""
+Spectral grating efficiencies
+"""
+import logging
+import numpy as np
+from matplotlib import pyplot as plt
+
+from astropy.io import fits
+from astropy import units as u
+from astropy.wcs import WCS
+from astropy.table import Table
+
+from .effects import Effect
+from .ter_curves import TERCurve
+from .ter_curves_utils import apply_throughput_to_cube
+from ..utils import find_file
+from ..base_classes import FieldOfViewBase, FOVSetupBase
+
+class SpectralEfficiency(Effect):
+    """
+    Applies the grating efficiency (blaze function) for a SpectralTraceList
+
+    Input Data Format
+    -----------------
+    The efficiency curves are taken from a fits file `filename`with a
+    structure similar to the trace definition file (see `SpectralTraceList`).
+    The required extensions are:
+    - 0 : PrimaryHDU [header]
+    - 1 : BinTableHDU or TableHDU[header, data] : Overview table of all traces
+    - 2..N : BinTableHDU or TableHDU : Efficiency curves, one per trace. The
+             tables must have the two columns `wavelength` and `efficiency`
+
+    Note that there must be one extension for each trace defined in the
+    `SpectralTraceList`. Extensions for other traces are ignored.
+
+    EXT 0 : PrimaryHDU
+    ++++++++++++++++++
+    Required header keywords:
+
+    - ECAT : int : Extension number of overview table, normally 1
+    - EDATA : int : Extension number of first Trace table, normally 2
+
+    No data is required in this extension
+
+    EXT 1 : (Bin)TableHDU : Overview of traces
+    ++++++++++++++++++++++++++++++++++++++++++
+    No special header keywords are required in this extension.
+
+    Required Table columns:
+    - description : str : identifier for each trace
+    - extension_id : int : which extension is each trace in
+
+    EXT 2 : (Bin)TableHDU : Efficiencies for individual traces
+    ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    Required header keywords:
+    - EXTNAME : must be identical to the `description` in EXT 1
+
+    Required Table columns:
+    - wavelength : float : [um]
+    - efficiency : float : number [0..1]
+
+    """
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        if "hdulist" in kwargs and isinstance(kwargs["hdulist"], fits.HDUList):
+            self._file = kwargs["hdulist"]
+
+        params = {"z_order": [630]}
+        self.meta.update(params)
+
+        self.efficiencies = self.get_efficiencies()
+
+    def get_efficiencies(self):
+        """Reads effciencies from file, returns a dictionary"""
+        hdul = self._file
+        self.ext_data = hdul[0].header["EDATA"]
+        self.ext_cat = hdul[0].header["ECAT"]
+        self.catalog = Table(hdul[self.ext_cat].data)
+
+        efficiencies = {}
+        for row in self.catalog:
+            params = {col: row[col] for col in row.colnames}
+            params.update(self.meta)
+            hdu = self._file[row["extension_id"]]
+            name = hdu.header['EXTNAME']
+
+            tbl = Table.read(hdu)
+            wavelength = tbl['wavelength'].quantity
+            efficiency = tbl['efficiency'].value
+            effic_curve = TERCurve(wavelength=wavelength,
+                                   transmission=efficiency,
+                                   **params)
+            efficiencies[name] = effic_curve
+
+        hdul.close()
+        return efficiencies
+
+    def apply_to(self, obj, **kwargs):
+        """
+        Interface between FieldOfView and SpectralEfficiency
+        """
+        trace_id = obj.meta['trace_id']
+        try:
+            effic = self.efficiencies[trace_id]
+        except KeyError:
+            logging.warning("No grating efficiency for trace %s" % trace_id)
+            return obj
+
+        wcs = WCS(obj.hdu.header).spectral
+        wave_cube = wcs.all_pix2world(np.arange(obj.hdu.data.shape[0]), 0)[0]
+        wave_cube = (wave_cube * u.Unit(wcs.wcs.cunit[0])).to(u.AA)
+        obj.hdu = apply_throughput_to_cube(obj.hdu, effic.throughput)
+        return obj
+
+    def plot(self):
+        """Plot the grating efficiencies"""
+        for name, effic in self.efficiencies.items():
+            wave = effic.throughput.waveset
+            plt.plot(wave.to(u.um), effic.throughput(wave), label=name)
+
+        plt.xlabel("Wavelength [um]")
+        plt.ylabel("Grating efficiency")
+        plt.title(f"Grating efficiencies from {self.filename}")
+        plt.legend()
+        plt.show()

scopesim/effects/spectral_trace_list.py

@@ -6,6 +6,8 @@
 """
 
 from pathlib import Path
+import logging
+
 import numpy as np
 
 from astropy.io import fits
@@ -61,17 +63,21 @@ class SpectralTraceList(Effect):
 
     Required Table columns:
 
-    - description : str : description of each each trace
+    - description : str : identifier of each trace
     - extension_id : int : which extension is each trace in
     - aperture_id : int : which aperture matches this trace (e.g. MOS / IFU)
     - image_plane_id : int : on which image plane is this trace projected
 
     EXT 2 : BinTableHDU : Individual traces
     +++++++++++++++++++++++++++++++++++++++
-    No special header keywords are required in this extension
+    Required header keywords:
+    - EXTNAME : must be identical to the `description` in EXT 1
 
-    Required Table columns:
+    Recommended header keywords:
+    - DISPDIR : 'x' or 'y' : dispersion axis. If not present, Scopesim tries
+      to determine this automatically; this may be unreliable in some cases.
 
+    Required Table columns:
     - wavelength : float : [um] : wavelength of monochromatic aperture image
     - s : float : [arcsec] : position along aperture perpendicular to trace
     - x : float : [mm] : x position of aperture image on focal plane
@@ -184,6 +190,7 @@ def apply_to(self, obj, **kwargs):
                 # for MAAT
                 pass
             elif obj.hdu is None and obj.cube is None:
+                logging.info("Making cube")
                 obj.cube = obj.make_cube_hdu()
 
             # ..todo: obj will be changed to a single one covering the full field of view

scopesim/effects/spectral_trace_list_utils.py

@@ -378,6 +378,9 @@ def plot(self, wave_min=None, wave_max=None, c="r"):
 
         # Footprint (rectangle enclosing the trace)
         xlim, ylim  = self.footprint(wave_min=wave_min, wave_max=wave_max)
+        if xlim is None:
+            return
+
         xlim.append(xlim[0])
         ylim.append(ylim[0])
         plt.plot(xlim, ylim)
@@ -397,14 +400,14 @@ def plot(self, wave_min=None, wave_max=None, c="r"):
             y = self.table[self.meta["y_colname"]][mask]
             plt.plot(x, y, "o", c=c)
 
-            for wave in np.unique(waves):
-                xx = x[waves==wave]
+            for wave in np.unique(w):
+                xx = x[w==wave]
                 xx.sort()
                 dx = xx[-1] - xx[-2]
-                plt.text(x[waves==wave].max() + 0.5 * dx,
-                         y[waves==wave].mean(),
-                         str(wave), va="center", ha="left")
 
+                plt.text(x[w==wave].max() + 0.5 * dx,
+                         y[w==wave].mean(),
+                         str(wave), va='center', ha='left')
 
             plt.gca().set_aspect("equal")
 

scopesim/tests/tests_effects/test_SpectralEfficiency.py

@@ -0,0 +1,37 @@
+"""Tests for class SpectralEfficiency"""
+import os
+import pytest
+
+from astropy.io import fits
+
+from scopesim import rc
+from scopesim.effects import SpectralEfficiency, TERCurve
+from scopesim.utils import find_file
+
+FILES_PATH = os.path.abspath(os.path.join(os.path.dirname(__file__),
+                                          "../mocks/files/"))
+if FILES_PATH not in rc.__search_path__:
+    rc.__search_path__ += [FILES_PATH]
+
+# pylint: disable=missing-class-docstring
+# pylint: disable=missing-function-docstring
+
+@pytest.fixture(name="speceff", scope="class")
+def fixture_speceff():
+    """Instantiate SpectralEfficiency object"""
+    return SpectralEfficiency(filename="TER_grating.fits")
+
+class TestSpectralEfficiency:
+    def test_initialises_from_file(self, speceff):
+        assert isinstance(speceff, SpectralEfficiency)
+
+    def test_initialises_from_hdulist(self):
+        fitsfile = find_file("TER_grating.fits")
+        hdul = fits.open(fitsfile)
+        speceff = SpectralEfficiency(hdulist=hdul)
+        assert isinstance(speceff, SpectralEfficiency)
+
+    def test_has_efficiencies(self, speceff):
+        efficiencies = speceff.efficiencies
+        assert all(isinstance(effic, TERCurve)
+                   for _, effic in efficiencies.items())