Merge pull request #182 from AstarVienna/dev_master

Prepare for 0.5.5 release

scopesim/effects/__init__.py

@@ -16,4 +16,6 @@
 from .obs_strategies import *
 from .fits_headers import *
 
+from .rotation import *
+
 # from . import effects_utils

scopesim/effects/apertures.py

@@ -443,12 +443,10 @@ def __init__(self, **kwargs):
 
         self.table = self.get_table()
 
-
     def apply_to(self, obj, **kwargs):
         """Use apply_to of current_slit"""
         return self.current_slit.apply_to(obj, **kwargs)
 
-
     def fov_grid(self, which="edges", **kwargs):
         return self.current_slit.fov_grid(which=which, **kwargs)
 
@@ -460,6 +458,21 @@ def change_slit(self, slitname=None):
         else:
             raise ValueError("Unknown slit requested: " + slitname)
 
+    def add_slit(self, newslit, name=None):
+        """
+        Add a slit to the SlitWheel
+
+        Parameters
+        ==========
+        newslit : Slit
+        name : string
+           Name to be used for the new slit. If `None`, a name from
+           the newslit object is used.
+        """
+        if name is None:
+            name = newslit.display_name
+        self.slits[name] = newslit
+
     @property
     def current_slit(self):
         """Return the currently used slit"""

scopesim/effects/detector_list.py

@@ -87,7 +87,7 @@ class : DetectorList
             active_detectors : [1, 3]
             image_plane_id : 0
 
-    where the file detecotr_list.dat contains the following information
+    where the file detector_list.dat contains the following information
     ::
 
         # x_cen_unit : mm

scopesim/effects/electronic.py

@@ -12,6 +12,7 @@
 - LinearityCurve - apply detector (non-)linearity and saturation
 - ReferencePixelBorder
 - BinnedImage
+- UnequalBinnedImage
 - Bias - adds constant bias level to readout
 
 Functions:
@@ -557,6 +558,24 @@ def apply_to(self, det, **kwargs):
 
         return det
 
+class UnequalBinnedImage(Effect):
+    def __init__(self, **kwargs):
+        super(UnequalBinnedImage, self).__init__(**kwargs)
+        self.meta["z_order"] = [870]
+
+        self.required_keys = ["binx","biny"]
+        utils.check_keys(self.meta, self.required_keys, action="error")
+
+    def apply_to(self, det, **kwargs):
+        if isinstance(det, DetectorBase):
+            bx = from_currsys(self.meta["binx"])
+            by = from_currsys(self.meta["biny"])
+            image = det._hdu.data
+            h, w = image.shape
+            new_image = image.reshape((h//bx, bx, w//by, by))
+            det._hdu.data = new_image.sum(axis=3).sum(axis=1)
+
+        return det
 
 ################################################################################
 

scopesim/effects/fits_headers.py

@@ -107,7 +107,7 @@ class ExtraFitsKeywords(Effect):
               ESO:
                 DET:
                   DIT: [5, '[s] exposure length']   # example of adding a comment
-            EXTNAME: "DET++.DATA"                   # example of extension specific qualifier
+            EXTNAME: "DET§.DATA"                   # example of extension specific qualifier
 
     The keywords can be added to one or more extensions, based on one of the
     following ``ext_`` qualifiers: ``ext_name``, ``ext_number``, ``ext_type``
@@ -116,11 +116,11 @@ class ExtraFitsKeywords(Effect):
     For a list, ScopeSim will add the keywords to all extensions matching the
     specified type/name/number
 
-    The number of the extension can be used in a value by using the "++" string.
-    That is, keyword values with "++" with have the extension number inserted
-    where the "++" is.
+    The number of the extension can be used in a value by using the "§" string.
+    That is, keyword values with "§" with have the extension number inserted
+    where the "§" is.
 
-    The above example (``EXTNAME: "DET++.DATA"``) will result in the following
+    The above example (``EXTNAME: "DET§.DATA"``) will result in the following
     keyword added only to extensions 1 and 2:
 
     - PrimaryHDU (ext 0)::
@@ -269,9 +269,10 @@ def apply_to(self, hdul, **kwargs):
                 unresolved = flatten_dict(dic.get("unresolved_keywords", {}))
                 exts = get_relevant_extensions(dic, hdul)
                 for i in exts:
+                    # On windows machines  appears in the string when using §
                     resolved_with_counters = {
-                        k: v.replace("++", str(i)) if isinstance(v, str) else v
-                        for k, v in resolved.items()
+                        k: v.replace("Â", "").replace("§", str(i)).replace("++", str(i))
+                        if isinstance(v, str) else v for k, v in resolved.items()
                     }
                     hdul[i].header.update(resolved_with_counters)
                     hdul[i].header.update(unresolved)
@@ -349,6 +350,11 @@ def flatten_dict(dic, base_key="", flat_dict=None,
                 comment = f"[{str(value.unit)}] " + comment
                 value = value.value
 
+            # Convert e.g.  Unit(mag) to just "mag". Not sure how this will
+            # work when deserializing though.
+            if isinstance(value, u.Unit):
+                value = str(value)
+
             if isinstance(value, (list, np.ndarray)):
                 value = f"{value.__class__.__name__}:{str(list(value))}"
                 max_len = 80 - len(flat_key)

scopesim/effects/psf_utils.py

@@ -122,11 +122,11 @@ def rescale_kernel(image, scale_factor, spline_order=None):
     if dy > 0:
         image = image[2*dy:, :]
     elif dy < 0:
-        image = image[:2*abs(dy), :]
+        image = image[:2*dy, :]
     if dx > 0:
         image = image[:, 2*dx:]
     elif dx < 0:
-        image = image[:, 2*abs(dx):]
+        image = image[:, :2*dx]
 
     sum_new_image = np.sum(image)
     image *= sum_image / sum_new_image

scopesim/effects/rotation.py

@@ -0,0 +1,36 @@
+"""Effects related to rotation of the field/CCD
+
+Classes:
+- RotateCCD - Rotates CCD by integer multiples of 90 degrees
+"""
+
+import logging
+import numpy as np
+
+from . import Effect
+from .. import utils
+from ..utils import from_currsys
+from ..base_classes import DetectorBase
+
+
+class Rotate90CCD(Effect):
+    """
+    Rotates CCD by integer multiples of 90 degrees
+    rotations kwarg is number of counter-clockwise rotations
+
+    """
+    def __init__(self, **kwargs):
+        super(Rotate90CCD, self).__init__(**kwargs)
+        params = {"z_order": [809]}
+        self.meta.update(params)
+        self.meta.update(kwargs)
+
+        required_keys = ["rotations"]
+        utils.check_keys(self.meta, required_keys, action="error")
+
+    def apply_to(self, obj, **kwargs):
+        if isinstance(obj, DetectorBase):
+            rotations = from_currsys(self.meta["rotations"])
+            obj._hdu.data = np.rot90(obj._hdu.data,rotations)
+
+        return obj

scopesim/effects/spectral_trace_list_utils.py

@@ -109,6 +109,11 @@ def compute_interpolation_functions(self):
         xi_arr = self.table[self.meta['s_colname']]
         lam_arr = self.table[self.meta['wave_colname']]
 
+        wi0, wi1 = lam_arr.argmin(), lam_arr.argmax()
+        x_disp_length = np.diff([x_arr[wi0], x_arr[wi1]])
+        y_disp_length = np.diff([y_arr[wi0], y_arr[wi1]])
+        self.dispersion_axis = "x" if x_disp_length > y_disp_length else "y"
+
         self.wave_min = quantify(np.min(lam_arr), u.um).value
         self.wave_max = quantify(np.max(lam_arr), u.um).value
 
@@ -203,8 +208,11 @@ def map_spectra_to_focal_plane(self, fov):
         #         - The dispersion direction is selected by the direction of the
         #           gradient of lam(x, y). This works if the lam-axis is well
         #           aligned with x or y. Needs to be tested for MICADO.
-        dlam_by_dx, dlam_by_dy = self.xy2lam.gradient()
-        if np.abs(dlam_by_dx(0, 0)) > np.abs(dlam_by_dy(0, 0)):
+
+
+        # dlam_by_dx, dlam_by_dy = self.xy2lam.gradient()
+        # if np.abs(dlam_by_dx(0, 0)) > np.abs(dlam_by_dy(0, 0)):
+        if self.dispersion_axis == "x":
             avg_dlam_per_pix = (wave_max - wave_min) / sub_naxis1
         else:
             avg_dlam_per_pix = (wave_max - wave_min) / sub_naxis2
@@ -259,6 +267,7 @@ def map_spectra_to_focal_plane(self, fov):
         image = xilam.interp(xi_fpa, lam_fpa, grid=False) * ijmask
 
         # Scale to ph / s / pixel
+        dlam_by_dx, dlam_by_dy = self.xy2lam.gradient()
         dlam_per_pix = pixsize * np.sqrt(dlam_by_dx(ximg_fpa, yimg_fpa)**2 +
                                          dlam_by_dy(ximg_fpa, yimg_fpa)**2)
         image *= pixscale * dlam_per_pix        # [arcsec/pix] * [um/pix]
@@ -414,7 +423,6 @@ def __repr__(self):
         return msg
 
 
-
 class XiLamImage():
     """
     Class to compute a rectified 2D spectrum
@@ -459,7 +467,8 @@ def __init__(self, fov, dlam_per_pix):
             # overlaps with the wavelength range covered by the cube
             if lam0.min() < cube_lam.max() and lam0.max() > cube_lam.min():
                 plane = fov.cube.data[:, i, :].T
-                plane_interp = RectBivariateSpline(cube_xi, cube_lam, plane)
+                plane_interp = RectBivariateSpline(cube_xi, cube_lam, plane,
+                                                   kx=1, ky=1)
                 self.image += plane_interp(cube_xi, lam0)
 
         self.image *= d_eta     # ph/s/um/arcsec2 --> ph/s/um/arcsec

scopesim/effects/ter_curves.py

@@ -1,3 +1,4 @@
+'''Transmission, emissivity, reflection curves'''
 import numpy as np
 from astropy import units as u
 from os import path as pth
@@ -355,7 +356,7 @@ class FilterCurve(TERCurve):
         ``TC_filter_{}.dat``
 
     Can either be created using the standard 3 options:
-    - ``filename``: direct filename of the filer curve
+    - ``filename``: direct filename of the filter curve
     - ``table``: an ``astropy.Table``
     - ``array_dict``: a dictionary version of a table: ``{col_name1: values, }``
 
@@ -487,10 +488,10 @@ def __init__(self, **kwargs):
         red = kwargs["red_cutoff"]
         peak = kwargs["transmission"]
         wing = kwargs.get("wing_transmission", 0)
-        
+
         waveset = [wave_min, 0.999*blue, blue, red, red*1.001, wave_max]
         transmission = [wing, wing, peak, peak, wing, wing]
-        
+
         tbl = Table(names=["wavelength", "transmission"],
                     data=[waveset, transmission])
         super(TopHatFilterCurve, self).__init__(table=tbl,
@@ -601,6 +602,21 @@ def change_filter(self, filtername=None):
         else:
             raise ValueError("Unknown filter requested: " + filtername)
 
+    def add_filter(self, newfilter, name=None):
+        """
+        Add a filter to the FilterWheel
+
+        Parameters
+        ==========
+        newfilter : FilterCurve
+        name : string
+           Name to be used for the new filter. If `None` a name from
+           the newfilter object is used.
+        """
+        if name is None:
+            name = newfilter.display_name
+        self.filters[name] = newfilter
+
     @property
     def current_filter(self):
         filter_eff = None

scopesim/source/source.py

@@ -144,6 +144,9 @@ def __init__(self, filename=None, cube=None, ext=0,
 
         self.meta = {}
         self.meta.update(kwargs)
+        # ._meta_dicts contains a meta for each of the .fields. It is primarily
+        # used to set proper FITS header keywords for each field so the source
+        # can be reconstructed from the FITS headers.
         self._meta_dicts = [self.meta]
 
         self.fields = []
@@ -490,7 +493,6 @@ def shift(self, dx=0, dy=0, layers=None):
             which .fields entries to shift
 
         """
-
         if layers is None:
             layers = np.arange(len(self.fields))
 
@@ -559,6 +561,13 @@ def make_copy(self):
 
     def append(self, source_to_add):
         new_source = source_to_add.make_copy()
+        # If there is no field yet, then self._meta_dicts contains a
+        # reference to self.meta, which is empty. This ensures that both are
+        # updated at the same time. However, it is important that the fields
+        # and _meta_dicts match when appending sources.
+        if len(self.fields) == 0:
+            assert self._meta_dicts == [{}]
+            self._meta_dicts = []
         if isinstance(source_to_add, Source):
             for field in new_source.fields:
                 if isinstance(field, Table):
@@ -601,4 +610,4 @@ def __repr__(self):
                     msg += f", referencing spectrum {num_spec}"
                 msg += "\n"
 
-        return msg
+        return msg

scopesim/tests/tests_effects/test_TERCurve.py

@@ -133,6 +133,16 @@ def test_change_to_unknown_filter(self, fwheel):
         with pytest.raises(ValueError):
             fwheel.change_filter('X')
 
+    def test_add_filter_to_wheel(self, fwheel):
+        num_filt_old = len(fwheel.filters)
+        newfilter = tc.TopHatFilterCurve(transmission=0.9, blue_cutoff=1.,
+                                         red_cutoff=2., name='blank')
+        fwheel.add_filter(newfilter, name='blank')
+        assert len(fwheel.filters) == num_filt_old + 1
+
+        fwheel.change_filter("blank")
+        assert fwheel.current_filter.display_name == "blank"
+
     def test_plots_all_filters(self, fwheel):
         if PLOTS:
             fwheel.plot()

scopesim/tests/tests_effects/test_fits_headers.py

@@ -176,21 +176,30 @@ def test_works_for_ext_type(self, comb_hdul, ext_type, answer):
 
 class TestFlattenDict:
     def test_works(self):
-        dic = {"HIERARCH":
-                   {"ESO":
-                        {"ATM":
-                             {"PWV": 1.0, "AIRMASS": 2.0},
-                         "DPR": {"TYPE": "DARK"}},
-                    "SIM": {
-                        "area": ("!TEL.area", "area")}
+        dic = {
+            "HIERARCH": {
+                "ESO": {
+                    "ATM": {
+                        "PWV": 1.0,
+                        "AIRMASS": 2.0,
+                    },
+                    "DPR": {
+                        "TYPE": "DARK",
                     }
-               }
+                },
+                "SIM": {
+                    "area": ("!TEL.area", "area"),
+                    "SRC0": {"scaling_unit": u.mag},
+                },
+           },
+        }
         flat_dict = fh.flatten_dict(dic)
         assert flat_dict["HIERARCH ESO ATM PWV"] == 1.0
         assert flat_dict["HIERARCH ESO ATM AIRMASS"] == 2.0
         assert flat_dict["HIERARCH ESO DPR TYPE"] == "DARK"
         assert flat_dict["HIERARCH SIM area"][0] == "!TEL.area"
         assert flat_dict["HIERARCH SIM area"][1] == "area"
+        assert flat_dict["HIERARCH SIM SRC0 scaling_unit"] == "mag"
 
     def test_resolves_bang_strings(self):
         dic = {"SIM": {"random_seed": "!SIM.random.seed"}}

scopesim/tests/tests_effects/test_slitwheel.py

@@ -49,3 +49,17 @@ def test_reports_current_slit_false(self):
     def test_changes_to_false(self, swheel):
         swheel.change_slit(False)
         assert not swheel.current_slit
+
+    def test_add_slit_to_wheel(self, swheel):
+        num_slit_old = len(swheel.slits)
+        kwargs = {"array_dict": {"x": [-2, -1, 1, 2],
+                                 "y": [-1, -2, 2, 1]},
+                  "x_unit": "arcsec",
+                  "y_unit": "arcsec"}
+        newslit = ApertureMask(name="newslit", **kwargs)
+
+        swheel.add_slit(newslit, name='newslit')
+        assert len(swheel.slits) == num_slit_old + 1
+
+        swheel.change_slit('newslit')
+        assert swheel.current_slit.display_name == "newslit"