From 3b877481284e8bb43763ec64847883b16ff64689 Mon Sep 17 00:00:00 2001 From: Hugo Buddelmeijer Date: Tue, 4 Jul 2023 19:55:44 +0200 Subject: [PATCH 1/4] Allow rescaling of SkycalcTERCurves --- scopesim/effects/ter_curves.py | 27 +++++-- scopesim/tests/mocks/files/TC_filter_Q.dat | 42 +++++++++++ .../tests_effects/test_SkycalcTERCurve.py | 74 ++++++++++++++++++- scopesim/utils.py | 11 +++ 4 files changed, 147 insertions(+), 7 deletions(-) create mode 100644 scopesim/tests/mocks/files/TC_filter_Q.dat diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py index 5d0e3335..55483354 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,15 @@ def query_server(self, **kwargs): logging.exception(msg) raise ValueError(msg) + # 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..ed97f78d --- /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 al 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 From 7d0378f1ee0b02be4e957fe699d925c1f529c072 Mon Sep 17 00:00:00 2001 From: Hugo Buddelmeijer Date: Wed, 5 Jul 2023 10:07:11 +0200 Subject: [PATCH 2/4] Add check whether wunit in self.meta Why is this even necessary?? The class should either have a wunit or not, but not like have it half of the time in meta --- scopesim/effects/ter_curves.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/scopesim/effects/ter_curves.py b/scopesim/effects/ter_curves.py index 55483354..e36d0d38 100644 --- a/scopesim/effects/ter_curves.py +++ b/scopesim/effects/ter_curves.py @@ -355,14 +355,15 @@ def query_server(self, **kwargs): logging.exception(msg) raise ValueError(msg) - # 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]) + 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 From c2d8a6e1066b9b0e3723f8cabe4c9b2c51d679aa Mon Sep 17 00:00:00 2001 From: Hugo Buddelmeijer Date: Wed, 5 Jul 2023 11:33:38 +0200 Subject: [PATCH 3/4] Empty commit to rerender the PR From 81bb307dc3085032d51dfb895a8028aa489ed861 Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 5 Jul 2023 21:35:32 +0200 Subject: [PATCH 4/4] Update TC_filter_Q.dat: typo --- scopesim/tests/mocks/files/TC_filter_Q.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scopesim/tests/mocks/files/TC_filter_Q.dat b/scopesim/tests/mocks/files/TC_filter_Q.dat index ed97f78d..9eea4d52 100644 --- a/scopesim/tests/mocks/files/TC_filter_Q.dat +++ b/scopesim/tests/mocks/files/TC_filter_Q.dat @@ -1,6 +1,6 @@ # This table is used in # test_SkycalcTERCurve.py::TestRescalingAndUnitConversion::test_rescaling -# It is essential that this table goes al the way to 2.500 for the +# It is essential that this table goes all the way to 2.500 for the # test to work properly. # wavelength transmission