From 030b771c3ff54ebf814d431294c4dba33be895fb Mon Sep 17 00:00:00 2001 From: Bryce Kalmbach Date: Tue, 14 Nov 2023 16:21:58 -0800 Subject: [PATCH 1/2] Add correct zenith angle into imsim configuration files. Add parallactic angle and zenith angle calculation into obsMetadata. --- doc/versionHistory.rst | 7 ++ policy/config/obsVariablesDefault.yaml | 2 +- python/lsst/ts/imsim/closed_loop_task.py | 12 +-- python/lsst/ts/imsim/imsim_cmpt.py | 6 +- python/lsst/ts/imsim/obs_metadata.py | 83 ++++++++++++++++++- python/lsst/ts/imsim/sky_sim.py | 25 ------ .../imsimConfig/imsimConfigLsstCam.yaml | 4 +- tests/test_obs_metadata.py | 53 +++++++++--- tests/test_sky_sim.py | 16 +--- 9 files changed, 141 insertions(+), 67 deletions(-) diff --git a/doc/versionHistory.rst b/doc/versionHistory.rst index 2e441c3..e44b032 100644 --- a/doc/versionHistory.rst +++ b/doc/versionHistory.rst @@ -6,6 +6,13 @@ Version History ################## +------------- +0.8.1 +------------- + +* Add correct zenith angle into imsim configuration files. +* Add parallactic angle and zenith angle calculation into obsMetadata. + ------------- 0.8.0 ------------- diff --git a/policy/config/obsVariablesDefault.yaml b/policy/config/obsVariablesDefault.yaml index b5f63e0..034eb62 100644 --- a/policy/config/obsVariablesDefault.yaml +++ b/policy/config/obsVariablesDefault.yaml @@ -17,7 +17,7 @@ eval_variables: # Zenith angle of the boresight. Note that we don't attempt to set this self consistently # based on boresight+MJD; that's up to the user. You can see down below which computations # explicitly depend on the zenith angle. - azenith: &zenith 0.0 deg + azenith: &zenith 41.407655 deg # Camera rotator angle. AKA RotTelPos artp: &rtp 0.0 deg diff --git a/python/lsst/ts/imsim/closed_loop_task.py b/python/lsst/ts/imsim/closed_loop_task.py index 6d444bb..58692ba 100644 --- a/python/lsst/ts/imsim/closed_loop_task.py +++ b/python/lsst/ts/imsim/closed_loop_task.py @@ -109,10 +109,7 @@ def config_sky_sim( self.sky_sim = SkySim() self.sky_sim.set_camera(inst_name) if path_sky_file == "": - par_angle = self.sky_sim.calc_parallactic_angle(obs_metadata) - self._set_sky_sim_based_on_opd_field_pos( - inst_name, obs_metadata, par_angle, star_mag - ) + self._set_sky_sim_based_on_opd_field_pos(inst_name, obs_metadata, star_mag) else: abs_sky_file_path = os.path.abspath(path_sky_file) self.sky_sim.add_star_by_file(abs_sky_file_path) @@ -121,7 +118,6 @@ def _set_sky_sim_based_on_opd_field_pos( self, inst_name: str, obs_metadata: ObsMetadata, - par_angle: float, star_mag: float, ) -> None: """Set the sky simulator based on the OPD field positions. @@ -134,8 +130,6 @@ def _set_sky_sim_based_on_opd_field_pos( Instrument name. obs_metadata : lsst.ts.imsim.ObsMetadata object Observation metadata. - par_angle : float - Parallactic angle. star_mag : float Star magnitude. This is to pretend there are the stars at OPD field positions. @@ -173,7 +167,9 @@ def _set_sky_sim_based_on_opd_field_pos( # https://lsstc.slack.com/archives/CHXKSF3HC/p1651863987821319?thread_ts=1651863934.274719&cid=CHXKSF3HC # that shows photons farthest from Zenith on sky appear on "top" # of focal plane. - rotation = rotMatrix(obs_metadata.rotator_angle - par_angle + 180) + rotation = rotMatrix( + obs_metadata.rotator_angle - obs_metadata.parallactic_angle + 180 + ) for ra_in_deg, dec_in_deg in zip(ra_in_deg_arr, dec_in_deg_arr): # It is noted that the field position might be < 0. But it is # not the same case for ra (0 <= ra <= 360). diff --git a/python/lsst/ts/imsim/imsim_cmpt.py b/python/lsst/ts/imsim/imsim_cmpt.py index 4a846cf..ee94220 100644 --- a/python/lsst/ts/imsim/imsim_cmpt.py +++ b/python/lsst/ts/imsim/imsim_cmpt.py @@ -230,7 +230,7 @@ def convert_obs_metadata_to_text(self, obs_metadata: ObsMetadata) -> str: obs_variables_text += f" ra: &ra {obs_metadata.ra} deg\n" obs_variables_text += f" dec: &dec {obs_metadata.dec} deg\n" obs_variables_text += f" sband: &band {obs_metadata.band}\n" - obs_variables_text += f" azenith: &zenith {obs_metadata.zenith} deg\n" + obs_variables_text += f" azenith: &zenith {obs_metadata.zenith:.6f} deg\n" obs_variables_text += f" artp: &rtp {obs_metadata.rotator_angle} deg\n" obs_variables_text += f" fexptime: &exptime {obs_metadata.exp_time}\n" obs_variables_text += f" fmjd: &mjd {obs_metadata.mjd}\n" @@ -297,7 +297,9 @@ def add_config_header(self, obs_metadata: ObsMetadata) -> str: header_text += f" fieldRA: {obs_metadata.ra}\n" header_text += f" fieldDec: {obs_metadata.dec}\n" header_text += f" rotTelPos: {obs_metadata.rotator_angle}\n" - header_text += f" airmass: {1.0/np.cos(obs_metadata.zenith)}\n" + header_text += ( + f" airmass: {1.0/np.cos(np.radians(obs_metadata.zenith)):.6f}\n" + ) header_text += f" focusZ: {obs_metadata.focus_z}\n" header_text += f" rawSeeing: {obs_metadata.raw_seeing}\n" diff --git a/python/lsst/ts/imsim/obs_metadata.py b/python/lsst/ts/imsim/obs_metadata.py index 8b45e11..7943ad3 100644 --- a/python/lsst/ts/imsim/obs_metadata.py +++ b/python/lsst/ts/imsim/obs_metadata.py @@ -21,7 +21,10 @@ __all__ = ["ObsMetadata"] -from dataclasses import dataclass +from dataclasses import dataclass, field + +import astropy +from astroplan import Observer @dataclass @@ -31,7 +34,6 @@ class ObsMetadata: ra: float # ra in degrees dec: float # dec in degrees band: str - zenith: float = 0.0 rotator_angle: float = 0.0 exp_time: float = 30.0 mjd: float = 59580.0 @@ -39,3 +41,80 @@ class ObsMetadata: raw_seeing: int = 0.5 obs_id: str = """$f"IM_P_{astropy.time.Time(mjd, format='mjd').strftime('%Y%m%d')}_{seqnum:06d}" """ focus_z: float = 0.0 # Defocal distance in mm + zenith: float = field(init=False) + parallactic_angle: float = field(init=False) + + def __post_init__(self) -> None: + """Populate the zenith and parallactic angles + with values from init. + """ + self.zenith = self.calc_zenith_angle() + self.parallactic_angle = self.calc_parallactic_angle() + + def format_observation_info(self) -> (Observer, astropy.time, astropy.coordinates): + """ + Get the observation info in the data structures needed + to calculate observer information such as parallactic angle + and zenith angle. + + Parameters + ---------- + obs_metadata : lsst.ts.imsim.ObsMetadata + ObsMetadata dataclass object with observation information. + + Returns + ------- + Observer + astroplan.Observer for observer located at Rubin Observatory + astropy.time + Observation time + astropy.coordinates + Observation boresight + """ + + # Observer located at Rubin + rubin_observer = Observer.at_site("cerro pachon") + time = astropy.time.Time(self.mjd, format="mjd") + boresight = astropy.coordinates.SkyCoord( + f"{self.ra}d", f"{self.dec}d", frame="icrs" + ) + + return rubin_observer, time, boresight + + def calc_parallactic_angle(self) -> float: + """Calculate the parallactic angle so we know the + sky rotation angle on alt-az mount for the observation. + + Parameters + ---------- + obs_metadata : lsst.ts.imsim.ObsMetadata + ObsMetadata dataclass object with observation information. + + Returns + ------- + float + Parallactic Angle in degrees. + """ + + observer, time, boresight = self.format_observation_info() + + return observer.parallactic_angle(time, boresight).deg + + def calc_zenith_angle(self) -> float: + """Calculate the zenith angle so we can accurately + populate the imsim configuration. + + Parameters + ---------- + obs_metadata : lsst.ts.imsim.ObsMetadata + ObsMetadata dataclass object with observation information. + + Returns + ------- + float + Zenith Angle in degrees. + """ + + observer, time, boresight = self.format_observation_info() + + return 90.0 - observer.altaz(time, boresight).alt.deg diff --git a/python/lsst/ts/imsim/sky_sim.py b/python/lsst/ts/imsim/sky_sim.py index 73bcc36..e83a096 100644 --- a/python/lsst/ts/imsim/sky_sim.py +++ b/python/lsst/ts/imsim/sky_sim.py @@ -21,10 +21,7 @@ __all__ = ["SkySim"] -import astropy import numpy as np -from astroplan import Observer -from lsst.ts.imsim.obs_metadata import ObsMetadata from lsst.ts.imsim.utils import get_camera @@ -58,28 +55,6 @@ def set_camera(self, inst_name: str) -> None: self._camera = get_camera(inst_name) - def calc_parallactic_angle(self, obs_metadata: ObsMetadata) -> float: - """Calculate the parallactic angle so we know the - sky rotation angle on alt-az mount for the observation. - - Parameters - ---------- - obs_metadata : lsst.ts.imsim.ObsMetadata - ObsMetadata dataclass object with observation information. - - Returns - ------- - float - Parallactic Angle in degrees. - """ - time = astropy.time.Time(obs_metadata.mjd, format="mjd") - rubin = Observer.at_site("cerro pachon") - boresight = astropy.coordinates.SkyCoord( - f"{obs_metadata.ra}d", f"{obs_metadata.dec}d", frame="icrs" - ) - - return rubin.parallactic_angle(time, boresight).deg - def add_star_by_ra_dec_in_deg( self, star_id: int | list[int] | np.ndarray, diff --git a/tests/testData/imsimConfig/imsimConfigLsstCam.yaml b/tests/testData/imsimConfig/imsimConfigLsstCam.yaml index a30ed36..fbe01f3 100644 --- a/tests/testData/imsimConfig/imsimConfigLsstCam.yaml +++ b/tests/testData/imsimConfig/imsimConfigLsstCam.yaml @@ -31,7 +31,7 @@ eval_variables: # Zenith angle of the boresight. Note that we don't attempt to set this self consistently # based on boresight+MJD; that's up to the user. You can see down below which computations # explicitly depend on the zenith angle. - azenith: &zenith 0.0 deg + azenith: &zenith 41.407655 deg # Camera rotator angle. AKA RotTelPos artp: &rtp 0.0 deg @@ -296,7 +296,7 @@ output: # nfiles: 197 # WFS + all science header: - airmass: 1.0 + airmass: 1.333293 band: r fieldDec: 0.0 fieldRA: 0.0 diff --git a/tests/test_obs_metadata.py b/tests/test_obs_metadata.py index 0efe11d..75d83a7 100644 --- a/tests/test_obs_metadata.py +++ b/tests/test_obs_metadata.py @@ -21,25 +21,54 @@ import unittest +from astroplan import FixedTarget, Observer +from astropy.time import Time from lsst.ts.imsim import ObsMetadata class TestObsMetadata(unittest.TestCase): """Test the ObsMetadata dataclass.""" + def setUp(self) -> None: + self.obs_meta_test = ObsMetadata(ra=0.0, dec=0.0, band="r") + def test_obs_metadata(self): - obs_meta_test = ObsMetadata(ra=0.0, dec=0.0, band="r") - self.assertEqual(obs_meta_test.ra, 0.0) - self.assertEqual(obs_meta_test.dec, 0.0) - self.assertEqual(obs_meta_test.band, "r") - self.assertEqual(obs_meta_test.zenith, 0.0) - self.assertEqual(obs_meta_test.rotator_angle, 0.0) - self.assertEqual(obs_meta_test.exp_time, 30.0) - self.assertEqual(obs_meta_test.raw_seeing, 0.5) - self.assertEqual(obs_meta_test.mjd, 59580.0) - self.assertEqual(obs_meta_test.seq_num, 1) + self.assertEqual(self.obs_meta_test.ra, 0.0) + self.assertEqual(self.obs_meta_test.dec, 0.0) + self.assertEqual(self.obs_meta_test.band, "r") + self.assertAlmostEqual(self.obs_meta_test.zenith, 41.407655076) + self.assertEqual(self.obs_meta_test.rotator_angle, 0.0) + self.assertEqual(self.obs_meta_test.exp_time, 30.0) + self.assertEqual(self.obs_meta_test.raw_seeing, 0.5) + self.assertEqual(self.obs_meta_test.mjd, 59580.0) + self.assertEqual(self.obs_meta_test.seq_num, 1) self.assertEqual( - obs_meta_test.obs_id, + self.obs_meta_test.obs_id, "$f\"IM_P_{astropy.time.Time(mjd, format='mjd').strftime('%Y%m%d')}_{seqnum:06d}\" ", ) - self.assertEqual(obs_meta_test.focus_z, 0.0) + self.assertEqual(self.obs_meta_test.focus_z, 0.0) + self.assertAlmostEqual(self.obs_meta_test.parallactic_angle, 139.4727348) + + def test_calc_parallactic_angle(self): + sirius = FixedTarget.from_name("sirius") + t = Time(60000, format="mjd") + obs_metadata = ObsMetadata( + ra=sirius.ra.deg, dec=sirius.dec.deg, band="r", mjd=60000 + ) + rubin = Observer.at_site("cerro pachon") + self.assertAlmostEqual( + obs_metadata.calc_parallactic_angle(), + rubin.parallactic_angle(t, sirius).deg, + ) + + def test_calc_zenith_angle(self): + sirius = FixedTarget.from_name("sirius") + t = Time(59580.0, format="mjd") + obs_metadata = ObsMetadata( + ra=sirius.ra.deg, dec=sirius.dec.deg, band="r", mjd=59580.0 + ) + rubin = Observer.at_site("cerro pachon") + self.assertAlmostEqual( + obs_metadata.calc_zenith_angle(), + 90.0 - rubin.altaz(t, sirius).alt.deg, + ) diff --git a/tests/test_sky_sim.py b/tests/test_sky_sim.py index fa6b4f4..7cb83e9 100644 --- a/tests/test_sky_sim.py +++ b/tests/test_sky_sim.py @@ -22,9 +22,7 @@ import os import unittest -from astroplan import FixedTarget, Observer -from astropy.time import Time -from lsst.ts.imsim import ObsMetadata, SkySim +from lsst.ts.imsim import SkySim from lsst.ts.imsim.utils.utility import get_module_path @@ -36,18 +34,6 @@ def test_set_camera(self): self.sky_sim.set_camera("lsstfam") self.assertEqual(self.sky_sim._camera.getName(), "LSSTCam") - def test_calc_parallactic_angle(self): - sirius = FixedTarget.from_name("sirius") - t = Time(60000, format="mjd") - obs_metadata = ObsMetadata( - ra=sirius.ra.deg, dec=sirius.dec.deg, band="r", mjd=60000 - ) - rubin = Observer.at_site("cerro pachon") - self.assertEqual( - self.sky_sim.calc_parallactic_angle(obs_metadata), - rubin.parallactic_angle(t, sirius).deg, - ) - def test_add_star_by_ra_dec_in_deg(self): self.sky_sim.add_star_by_ra_dec_in_deg(1, 2, 3, 4) self.assertEqual(len(self.sky_sim.star_id), 1) From d86dc0d44414cc1f6075ec931c8d4308edf3926c Mon Sep 17 00:00:00 2001 From: Bryce Kalmbach Date: Wed, 15 Nov 2023 15:05:59 -0800 Subject: [PATCH 2/2] Add comment on conversion from alt to zenith. Change order of attributes in obs_metadata test. --- python/lsst/ts/imsim/obs_metadata.py | 3 +++ tests/test_obs_metadata.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/python/lsst/ts/imsim/obs_metadata.py b/python/lsst/ts/imsim/obs_metadata.py index 7943ad3..ded25d8 100644 --- a/python/lsst/ts/imsim/obs_metadata.py +++ b/python/lsst/ts/imsim/obs_metadata.py @@ -117,4 +117,7 @@ def calc_zenith_angle(self) -> float: observer, time, boresight = self.format_observation_info() + # Altitude refers to elevation angle up from horizon. + # To get zenith we need to convert to the angle down + # from zenith so we subtract altitude from 90 degrees. return 90.0 - observer.altaz(time, boresight).alt.deg diff --git a/tests/test_obs_metadata.py b/tests/test_obs_metadata.py index 75d83a7..cad4f06 100644 --- a/tests/test_obs_metadata.py +++ b/tests/test_obs_metadata.py @@ -36,7 +36,6 @@ def test_obs_metadata(self): self.assertEqual(self.obs_meta_test.ra, 0.0) self.assertEqual(self.obs_meta_test.dec, 0.0) self.assertEqual(self.obs_meta_test.band, "r") - self.assertAlmostEqual(self.obs_meta_test.zenith, 41.407655076) self.assertEqual(self.obs_meta_test.rotator_angle, 0.0) self.assertEqual(self.obs_meta_test.exp_time, 30.0) self.assertEqual(self.obs_meta_test.raw_seeing, 0.5) @@ -47,6 +46,7 @@ def test_obs_metadata(self): "$f\"IM_P_{astropy.time.Time(mjd, format='mjd').strftime('%Y%m%d')}_{seqnum:06d}\" ", ) self.assertEqual(self.obs_meta_test.focus_z, 0.0) + self.assertAlmostEqual(self.obs_meta_test.zenith, 41.407655076) self.assertAlmostEqual(self.obs_meta_test.parallactic_angle, 139.4727348) def test_calc_parallactic_angle(self):