diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py index 5d0e3335..e36d0d38 100644 --- a/scopesim/effects/ter_curves.py +++ b/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 diff --git a/scopesim/tests/mocks/files/TC_filter_Q.dat b/scopesim/tests/mocks/files/TC_filter_Q.dat new file mode 100644 index 00000000..9eea4d52 --- /dev/null +++ b/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 diff --git a/scopesim/tests/tests_effects/test_SkycalcTERCurve.py b/scopesim/tests/tests_effects/test_SkycalcTERCurve.py index ca696de6..61f0423c 100644 --- a/scopesim/tests/tests_effects/test_SkycalcTERCurve.py +++ b/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 diff --git a/scopesim/utils.py b/scopesim/utils.py index eafb6f82..0bc2c8b4 100644 --- a/scopesim/utils.py +++ b/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