Implement basic uvw-by-target tracking.

TO DO:
 - Flatten uvw-map, producing matched utc / lha arrays.
 - Output all as table.

Author: timstaley

src/fastimgproto/coords.py

@@ -2,7 +2,7 @@
 import numpy as np
 
 
-def time_of_next_transit(start_time, observer_longitude, target_ra,
+def time_of_next_transit(observer_longitude, target_ra, start_time,
                          tolerance=0.02 * u.arcsec):
     """
     Calculate time, ``t`` when the local-hour-angle of target will be 0.
@@ -15,11 +15,11 @@ def time_of_next_transit(start_time, observer_longitude, target_ra,
     time-period.
 
     Args:
-        start_time (astropy.time.Time): Time to start from.
         observer_longitude (astropy.coordinates.Longitude): Longitude
             of position on Earth
         target_ra (astropy.coordinates.Longitude): Right ascension of
             sky-target.
+        start_time (astropy.time.Time): Time to start from.
         tolerance (astropy.units.Quantity): If target's LHA is within
             ``tolerance`` of zero at ``start_time``, simply return
             ``start_time``. Otherwise look for the next transit.

src/fastimgproto/telescope/base.py

@@ -1,5 +1,6 @@
 from __future__ import print_function
 
+from collections import OrderedDict
 from itertools import combinations
 
 import astropy.io.ascii as ascii
@@ -9,10 +10,10 @@
 from astropy.coordinates import (Latitude, Longitude, EarthLocation, )
 from astropy.table import Table
 from attr import attrib, attrs
-
 from fastimgproto.coords import (
+    time_of_next_transit,
     xyz_to_uvw_rotation_matrix,
-    z_rotation_matrix
+    z_rotation_matrix,
 )
 
 _validator_optional_ndarray = attr.validators.optional(
@@ -98,6 +99,19 @@ def from_itrf(ant_itrf_xyz, ant_labels):
             ant_local_xyz=ant_local_xyz,
         )
 
+    def lha(self, ra, time):
+        """
+        Calculate the local hour angle of a target-RA at a given time
+
+        Args:
+            ra (astropy.coordinates.Longitude): Right ascension
+            time (astropy.time.Time): Timestamp
+
+        Returns:
+            astropy.coordinates.Longitude: Local Hour Angle of target-RA.
+        """
+        return self.lst(time) - ra
+
     def lst(self, time):
         """
         Calculate the local sidereal time at the telescope
@@ -106,8 +120,29 @@ def lst(self, time):
             time (astropy.time.Time): Global timestamp
 
         Returns:
+            astropy.coordinates.Longitude: Local sidereal time expressed as
+            an angle-Quantity.
+        """
+        return time.sidereal_time('apparent',
+                                  longitude=self.longitude)
+
+    def next_transit(self, target_ra, start_time, ):
+        """
+        Wrapper around :py:func:`.time_of_next_transit`
+
+        See :py:func:`.time_of_next_transit` for details.
+
+        Args:
+            target_ra (astropy.coordinates.Longitude): Right ascension of
+                sky-target.
+            start_time (astropy.time.Time): Time to start from.
+        Returns:
+            astropy.time.Time: Approximate time of the next transit
 
         """
+        return time_of_next_transit(observer_longitude=self.longitude,
+                                    target_ra=target_ra,
+                                    start_time=start_time)
 
     def uvw_at_local_hour_angle(self, local_hour_angle, dec):
         """
@@ -128,8 +163,14 @@ def uvw_at_local_hour_angle(self, local_hour_angle, dec):
         rotation = xyz_to_uvw_rotation_matrix(local_hour_angle, dec)
         return np.dot(rotation, self.baseline_local_xyz.T).T
 
-    def uvw_at_skycoord_and_time(self, pointing_centre, utc_time):
-        pass
+    def uvw_tracking_skycoord(self, pointing_centre, obs_times):
+        lha_uvw_map = OrderedDict()
+        for time in obs_times:
+            lha = self.lha(pointing_centre.ra, time)
+            lha_uvw_map[lha] = self.uvw_at_local_hour_angle(
+                local_hour_angle=lha, dec=pointing_centre.dec
+            )
+        return lha_uvw_map
 
 
 def generate_baselines_and_labels(antenna_positions, antenna_labels):
@@ -238,4 +279,3 @@ def parse_itrf_ascii_to_xyz_and_labels(tabledata):
     xyz = hstack_table_columns_as_ndarray(tbl.columns[:3]) * u.m
     labels = tbl['label'].data
     return xyz, labels
-

tests/test_coords/test_transit_calcs.py

@@ -17,27 +17,23 @@ def test_basic_functionality():
     tol = 0.02 * u.arcsec
 
     # Next transit-time at lon=0
-    tt_lon0 = time_of_next_transit(start_time=t0,
-                                   observer_longitude=lon0,
-                                   target_ra=ra,
-                                   tolerance=tol)
+    tt_lon0 = time_of_next_transit(observer_longitude=lon0, target_ra=ra,
+                                   start_time=t0, tolerance=tol)
 
     tt1_lst = tt_lon0.sidereal_time('apparent', longitude=lon0)
     assert np.fabs(tt1_lst) < tol
 
     # If we provide a close-enough-to-transit time, just return it, don't
     # endlessly cycle by a sidereal day:
-    tt_lon0_tnt = time_of_next_transit(tt_lon0, lon0, ra, tol)
+    tt_lon0_tnt = time_of_next_transit(lon0, ra, tt_lon0, tol)
     assert tt_lon0_tnt == tt_lon0
 
     # Longitude is positive east. So a negative longitude observer (West) should
     # see a later transit:
 
     lon_m30 = Longitude(-30 * u.deg)
-    tt_lon30 = time_of_next_transit(start_time=t0,
-                                    observer_longitude=lon_m30,
-                                    target_ra=ra,
-                                    tolerance=tol)
+    tt_lon30 = time_of_next_transit(observer_longitude=lon_m30, target_ra=ra,
+                                    start_time=t0, tolerance=tol)
     assert tt_lon30 > tt_lon0
     # Difference should be 30 / 360 = 2 sidereal hours
     shour = 1./24.*u.sday # Sidereal hour

tests/test_telescope/test_telescope.py

@@ -1,10 +1,12 @@
 from __future__ import print_function
 import numpy as np
 import astropy.units as u
-from astropy.coordinates import Latitude, Longitude
+from astropy.coordinates import Latitude, Longitude, SkyCoord
+from astropy.time import Time
 from fastimgproto.telescope import Telescope
 from fastimgproto.telescope.base import parse_itrf_ascii_to_xyz_and_labels
 from fastimgproto.telescope.data import meerkat_itrf_filepath
+from fastimgproto.telescope.readymade import Meerkat
 
 
 def test_telescope_from_itrf():
@@ -38,9 +40,9 @@ def test_xyz_antennae_to_uvw_baselines():
     ]
 
     expected_uvw_baselines = np.array([
-        [1., 0., 0.], #origin,1-east
-        [0., 1., 0.], #origin,1-north
-        [-1., 1., 0.]]) #1-east,1-north
+        [1., 0., 0.],  # origin,1-east
+        [0., 1., 0.],  # origin,1-north
+        [-1., 1., 0.]])  # 1-east,1-north
 
     tel = Telescope(latitude=Latitude(0 * u.deg),
                     longitude=Longitude(0 * u.deg),
@@ -53,7 +55,30 @@ def test_xyz_antennae_to_uvw_baselines():
     lha = 0 * u.deg
     dec = 0 * u.deg
     uvw = tel.uvw_at_local_hour_angle(lha, dec)
-    assert (uvw==expected_uvw_baselines).all()
+    assert (uvw == expected_uvw_baselines).all()
     # print("UVW:")
     # print(repr(uvw))
     # print(tel.baseline_labels)
+
+
+def test_uvw_tracking_skyposn():
+    """
+    Check UVW calculations for Time+RA+Dec vs LHA+Dec produce consistent results
+    """
+    telescope = Meerkat()
+    azimuth_target = SkyCoord(ra=0 * u.deg, dec=telescope.latitude)
+
+    t0 = Time('2017-01-01 12:00:00', format='iso', scale='utc')
+    tt = telescope.next_transit(azimuth_target.ra, t0)
+    tt_lha = telescope.lha(azimuth_target.ra, tt)
+    assert np.fabs(tt_lha) < 0.2 * u.arcsec
+
+    obs_times = [tt - 10*u.second, tt, tt + 10*u.second]
+    uvw_by_tracking = telescope.uvw_tracking_skycoord(
+        azimuth_target,
+        obs_times=obs_times)
+    assert len(uvw_by_tracking) == len(obs_times)
+
+    uvw_by_lha = telescope.uvw_at_local_hour_angle(tt_lha, azimuth_target.dec)
+    assert uvw_by_tracking.keys()[1] == tt_lha
+    assert (uvw_by_lha == uvw_by_tracking.values()[1]).all()