LSF works

Author: oczoske

scopesim/effects/metis_ifu_simple/ifu_simple.py

@@ -3,12 +3,16 @@
 
 from typing import ClassVar
 
+import numpy as np
+from astropy.convolution import convolve, Box1DKernel, Gaussian1DKernel
+
 from ..effects import Effect
 from ...optics.fov import FieldOfView
 from ...utils import from_currsys, get_logger
 
 logger = get_logger(__name__)
 
+
 class LineSpreadFunction(Effect):
     """
     Compute and apply line spread function to IFU cube
@@ -26,13 +30,52 @@ def __init__(self, **kwargs):
         self.meta.update(kwargs)
         self.meta = from_currsys(self.meta, self.cmds)
 
-        self.lsfwidth = None #self.get_width()
-        self.kernel = None   #self.get_kernel()
+        if "lsfwidth" not in self.meta:
+            self.lsfwidth = self.get_lsf_width()
+        else:
+            self.lsfwidth = self.meta["lsfwidth"]
+        self.kernel = self.get_kernel()
 
-    def apply_to(self, obj, **kwargs):
+    def apply_to(self, fov, **kwargs):
         """Apply the LSF"""
-        if not isinstance(obj, FieldOfView):
-            return obj
+        if not isinstance(fov, FieldOfView):
+            return fov
+
+        if fov.hdu is None or fov.hdu.header["NAXIS"] != 3:
+            logger.error("Cannot apply LSF convolution")
+            return fov
+
+
+        print("Happily convolving with LSF")
+        fov.hdu.data = convolve(fov.hdu.data, self.kernel)
+        fov.cube = fov.hdu
+        return fov
+
+    def get_lsf_width(self):
+        """Determine width of the LSF kernel at central wavelength"""
+
+        slope = self.meta['fit_slope']
+        intercept = self.meta['fit_intercept']
+        lamc = self.meta["wavelen"]
+        dlam_per_pix = slope * lamc + intercept
+
+        slice_width = self.meta['slice_width']
+        pixel_scale = self.meta['pixel_scale']
+
+        dlam_per_slice = dlam_per_pix * slice_width / pixel_scale
+
+        spec_binwidth = self.meta["spec_binwidth"]
+
+        return dlam_per_slice / spec_binwidth
+
+
+    def get_kernel(self):
+        """Build LSF kernel: box kernel smoothed with narrow Gauss"""
 
-        lamc = from_currsys("!OBS.wavelen", self.cmds)
-        print(">>>>>>>>> Central wavelength:", lamc)
+        box = Box1DKernel(width=self.lsfwidth)
+        gauss = Gaussian1DKernel(1)
+        if box.shape > gauss.shape:
+            kernel = convolve(box.array, gauss)[np.newaxis, np.newaxis, :]
+        else:
+            kernel = convolve(gauss.array, box)[np.newaxis, np.newaxis, :]
+        return kernel

scopesim/tests/tests_effects/test_MetisIfuSimple.py

@@ -4,8 +4,26 @@
 # pylint: disable=missing-class-docstring
 # pylint: disable=missing-function-docstring
 
+import pytest
+import numpy as np
 from scopesim.effects.metis_ifu_simple import LineSpreadFunction
 
 class TestLineSpreadFunction:
     def test_initialises_correctly(self):
-        assert isinstance(LineSpreadFunction(), LineSpreadFunction)
+        assert isinstance(LineSpreadFunction(lsfwidth=2.7), LineSpreadFunction)
+
+    @pytest.mark.parametrize("width", [1.5, 2.7, 3.0, 5.2])
+    def test_kernel_is_normalised(self, width):
+        lsf = LineSpreadFunction(lsfwidth=width)
+        assert np.isclose(lsf.kernel.sum(), 1, 0.02)
+
+
+    def test_kernel_is_normalised_with_params(self):
+        lsf = LineSpreadFunction(wavelen=3.0,
+                                 fit_slope=3.795e-6,
+                                 fit_intercept=-4.659e-7,
+                                 slice_width=0.0207,
+                                 pixel_scale=0.0082,
+                                 spec_binwidth=1e-5)
+        assert isinstance(lsf, LineSpreadFunction)
+        assert np.isclose(lsf.kernel.sum(), 1, 1e-4)