Merge pull request #20 from lsst-ts/tickets/DM-41714

Add correct zenith angle into imsim configuration files.

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
 -------------

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
 

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).

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"
 

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,11 +34,90 @@ 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
     seq_num: int = 1
     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()
+
+        # 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

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,

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

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.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.zenith, 41.407655076)
+        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,
+        )

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)