Fix some issues with SkyCalcTERCurve

scopesim/effects/ter_curves.py

@@ -6,7 +6,7 @@
 import skycalc_ipy
 from astropy import units as u
 from astropy.io import fits
-from astropy.table import Table
+from astropy.table import Table, vstack
 
 from .effects import Effect
 from .ter_curves_utils import add_edge_zeros
@@ -15,7 +15,7 @@
 from ..base_classes import SourceBase, FOVSetupBase
 from ..optics.surface import SpectralSurface
 from ..source.source import Source
-from ..utils import from_currsys, quantify, check_keys, find_file
+from ..utils import from_currsys, quantify, check_keys, find_file, save_unit_to
 
 
 class TERCurve(Effect):
@@ -326,14 +326,22 @@ def load_skycalc_table(self):
         self.surface.meta.update(tbl.meta)
         self.skycalc_table = tbl
 
-    def query_server(self, **kwargs):
-        self.meta.update(kwargs)
+    def convert_to_nm(self):
+        """The Skycalc server talks in nm, so convert to that.
 
+        It talks back in um though.
+        """
         if "wunit" in self.meta:
-            scale_factor = u.Unit(from_currsys(self.meta["wunit"])).to(u.nm)
+            scale_factor = save_unit_to(u.Unit(from_currsys(self.meta["wunit"])), u.nm)
             for key in ["wmin", "wmax", "wdelta"]:
                 if key in self.meta:
                     self.meta[key] = from_currsys(self.meta[key]) * scale_factor
+            self.meta["wunit"] = "nm"
+
+    def query_server(self, **kwargs):
+        self.meta.update(kwargs)
+
+        self.convert_to_nm()
 
         conn_kwargs = {key: self.meta[key] for key in self.meta
                        if key in self.skycalc_conn.defaults}
@@ -347,6 +355,16 @@ def query_server(self, **kwargs):
             logging.exception(msg)
             raise ValueError(msg)
 
+        if "wunit" in self.meta:
+            # The last row is probably just short of wmax. So just copy
+            # the last row, but then with lam=wmax.
+            tbl_last_row = tbl[-1:].copy()
+            lam_last = tbl_last_row[-1]["lam"] * tbl_last_row["lam"].unit
+            lam_wmax = self.meta["wmax"] * u.Unit(self.meta["wunit"])
+            if lam_last < lam_wmax:
+                tbl_last_row[-1]["lam"] = lam_wmax.to(tbl_last_row["lam"].unit).value
+                tbl = vstack([tbl, tbl_last_row])
+
         return tbl
 
 

scopesim/tests/mocks/files/TC_filter_Q.dat

@@ -0,0 +1,42 @@
+# This table is used in
+# test_SkycalcTERCurve.py::TestRescalingAndUnitConversion::test_rescaling
+# It is essential that this table goes all the way to 2.500 for the
+# test to work properly.
+#
+wavelength      transmission
+     0.780       0.000100000
+     1.895       0.000100000
+     1.932       0.000403475
+     1.940        0.00122141
+     1.948        0.00273684
+     1.949        0.00377545
+     1.954        0.00896850
+     1.955         0.0100071
+     1.959         0.0282328
+     1.960         0.0349870
+     1.961         0.0417412
+     1.964         0.0735057
+     1.965          0.101068
+     1.967          0.187416
+     1.968          0.279336
+     1.969          0.402423
+     1.970          0.580067
+     1.971          0.833265
+     1.972          0.937967
+     2.327          0.940656
+     2.328          0.903457
+     2.329          0.756261
+     2.330          0.600915
+     2.331          0.422942
+     2.333          0.213269
+     2.334          0.145727
+     2.335          0.121079
+     2.336         0.0964298
+     2.343         0.0329677
+     2.349         0.0102813
+     2.350        0.00928342
+     2.351        0.00828549
+     2.380       0.000388406
+     2.402       0.000132468
+     2.435       0.000100000
+     2.500       0.000100000

scopesim/tests/tests_effects/test_SkycalcTERCurve.py

@@ -1,15 +1,16 @@
 import pytest
 import os
 from synphot import SpectralElement, SourceSpectrum
+from astropy import units as u
 
 from scopesim.effects import SkycalcTERCurve
 from scopesim import rc
-from scopesim.utils import from_currsys
+from scopesim.utils import save_unit_to
 
 if rc.__config__["!SIM.tests.run_skycalc_ter_tests"] is False:
     pytestmark = pytest.mark.skip("Ignoring SkyCalc integration tests")
 
-FILES_PATH = os.path.join(os.path.dirname(__file__), "../MOCKS/files/")
+FILES_PATH = os.path.join(os.path.dirname(__file__), "../mocks/files/")
 if FILES_PATH not in rc.__search_path__:
     rc.__search_path__ += [FILES_PATH]
 
@@ -53,3 +54,72 @@ def test_initialises_with_non_skycalc_keys(self):
     def test_initialise_with_local_skycalc_file(self):
         sky_ter = SkycalcTERCurve(use_local_skycalc_file="skycalc_override.fits")
         assert sky_ter.skycalc_table is not None
+
+
+class TestRescalingAndUnitConversion:
+    """Test rescaling of SkyCalcTERCurves and unit conversion."""
+
+    def test_convert_to_nm(self):
+        """Can we idempotent convert to nm?"""
+        tc = SkycalcTERCurve(
+            observatory="armazones",
+            wmin=0.7,
+            wmax=2.5,
+            wunit="um",
+        )
+        tc.convert_to_nm()
+        wmin1 = tc.meta["wmin"]
+        wmax1 = tc.meta["wmax"]
+        wunit1 = tc.meta["wunit"]
+        assert wmin1 == 700
+        assert wmax1 == 2500
+        assert wunit1 == "nm"
+        tc.convert_to_nm()
+        assert wmin1 == 700
+        assert wmax1 == 2500
+        assert wunit1 == "nm"
+
+    def test_save_to(self):
+        """Test whether exact unit conversion works."""
+        for ufrom, uto, fac in [
+            (u.um, u.nm, 1000.),
+            (u.mm, u.um, 1000.),
+            (u.nm, u.um, 0.001),
+            (u.um, u.mm, 0.001),
+        ]:
+            assert save_unit_to(ufrom, uto) == fac
+
+    def test_rescaling(self):
+        """Test rescaling of SkyCalcTERCurves and unit conversion.
+
+        There is a unit conversion problem that prevents rescaling of
+        SkyCalcTERCurves.
+
+        When wunit="um" and wmax=2.5, then this gets converted to nm by astropy,
+        because that is what Skycalc expects. However,
+        u.um.to(u.nm) = 999.9999999999999
+        and therefore 2.5 um is converted to 2499.99 nm.
+
+        2499.99 nm is less than 2.5 um, so there is no overlap anymore between
+        the source spectrum and the to-be-scaled to spectrum, leading to an
+        error.
+        """
+        tc = SkycalcTERCurve(
+            observatory="armazones",
+            wmin=0.7,
+            # It is essential that wmax=2.5 exactly for this test to work properly.
+            wmax=2.5,
+            wunit="um",
+            rescale_emission={
+                # The Q filter also goes to 2.5 exactly
+                "filter_name": "Q",
+                "filename_format": "TC_filter_{}.dat",
+                "value": 17,
+                "unit": "mag",
+            }
+        )
+        # Accessing tc.emission will lead to an exception:
+        # E               synphot.exceptions.PartialOverlap: Source spectrum
+        # and bandpass do not fully overlap. You may use force=[extrap|taper]
+        # to force this Observation anyway.
+        _ = tc.emission

scopesim/utils.py

@@ -1008,3 +1008,14 @@ def return_latest_github_actions_jobs_status(owner_name="AstarVienna", repo_name
         params_list.append(params)
 
     return params_list
+
+
+def save_unit_to(unit_from: u.Unit, unit_to: u.Unit):
+    """Like unit_from.to(unit_to), but u.um to u.nm == 1000, not 999.999.."""
+    to1 = unit_from.to(unit_to)
+    sto1 = str(to1)
+    to2 = 1 / unit_to.to(unit_from)
+    sto2 = str(to2)
+    if len(sto2) < len(sto1):
+        return to2
+    return to1