Add arbitrary epoch

katpoint/__init__.py

@@ -14,7 +14,7 @@
 
 from .target import Target, construct_azel_target, construct_radec_target
 from .antenna import Antenna
-from .timestamp import Timestamp
+from .timestamp import Timestamp, epoch_to_ephem
 from .flux import FluxDensityModel
 from .catalogue import Catalogue, specials, _catalogue_completer
 from .ephem_extra import lightspeed, rad2deg, deg2rad, wrap_angle, is_iterable

katpoint/target.py

@@ -3,7 +3,7 @@
 import numpy as np
 import ephem
 
-from .timestamp import Timestamp
+from .timestamp import Timestamp, epoch_to_ephem
 from .flux import FluxDensityModel
 from .ephem_extra import StationaryBody, NullBody, is_iterable, lightspeed, deg2rad, rad2deg
 from .conversion import azel_to_enu
@@ -361,19 +361,22 @@ def _scalar_radec(t):
         else:
             return _scalar_radec(timestamp)
 
-    def astrometric_radec(self, timestamp=None, antenna=None):
+    def astrometric_radec(self, timestamp=None, antenna=None, epoch='J2000'):
         """Calculate target's astrometric (ra, dec) coordinates as seen from antenna at time(s).
 
-        This calculates the J2000 *astrometric geocentric position* of the
-        target, in equatorial coordinates. This is its star atlas position for
-        the epoch of J2000.
+        This calculates the *astrometric geocentric position* of the target, in equatorial coordinates
+        at the specified epoch. This is its star atlas position for the given epoch.
 
         Parameters
         ----------
         timestamp : :class:`Timestamp` object or equivalent, or sequence, optional
             Timestamp(s) in UTC seconds since Unix epoch (defaults to now)
         antenna : :class:`Antenna` object, optional
             Antenna which points at target (defaults to default antenna)
+        epoch : string, optional
+            The desired epoch of the returned coordinates, either B1900 or B1950 
+            or can be of the form 'Jxxxx.xx' which specifies the date in terms
+            of Julian years from J2000 epoch.
 
         Returns
         -------
@@ -389,9 +392,9 @@ def astrometric_radec(self, timestamp=None, antenna=None):
 
         """
         if self.body_type == 'radec':
-            # Convert to J2000 equatorial coordinates
+            # Convert coordinates to specified epoch
             original_radec = ephem.Equatorial(self.body._ra, self.body._dec, epoch=self.body._epoch)
-            ra, dec = ephem.Equatorial(original_radec, epoch=ephem.J2000).get()
+            ra, dec = ephem.Equatorial(original_radec, epoch=epoch_to_ephem(epoch)).get()
             if is_iterable(timestamp):
                 return np.tile(ra, len(timestamp)), np.tile(dec, len(timestamp))
             else:
@@ -400,6 +403,7 @@ def astrometric_radec(self, timestamp=None, antenna=None):
         def _scalar_radec(t):
             """Calculate (ra, dec) coordinates for a single time instant."""
             antenna.observer.date = Timestamp(t).to_ephem_date()
+            antenna.observer.epoch = epoch_to_ephem(epoch)
             self.body.compute(antenna.observer)
             return self.body.a_ra, self.body.a_dec
         if is_iterable(timestamp):
@@ -411,19 +415,23 @@ def _scalar_radec(t):
     # The default (ra, dec) coordinates are the astrometric ones
     radec = astrometric_radec
 
-    def galactic(self, timestamp=None, antenna=None):
+    def galactic(self, timestamp=None, antenna=None, epoch='J2000'):
         """Calculate target's galactic (l, b) coordinates as seen from antenna at time(s).
 
         This calculates the galactic coordinates of the target, based on the
-        J2000 *astrometric* equatorial coordinates. This is its position relative
-        to the Galactic reference frame for the epoch of J2000.
+        *astrometric* equatorial coordinates. This is its position relative
+        to the Galactic reference frame for the given epoch.
 
         Parameters
         ----------
         timestamp : :class:`Timestamp` object or equivalent, or sequence, optional
             Timestamp(s) in UTC seconds since Unix epoch (defaults to now)
         antenna : :class:`Antenna` object, optional
             Antenna which points at target (defaults to default antenna)
+        epoch : string, optional
+            The desired epoch of the returned coordinates, either B1900 or B1950 
+            or can be of the form 'Jxxxx.xx' which specifies the date in terms
+            of Julian years from J2000 epoch.
 
         Returns
         -------
@@ -439,12 +447,15 @@ def galactic(self, timestamp=None, antenna=None):
 
         """
         if self.body_type == 'gal':
-            l, b = ephem.Galactic(ephem.Equatorial(self.body._ra, self.body._dec)).get()
+            #Get J2000 equatorial coordinates
+            original_radec = ephem.Equatorial(self.body._ra, self.body._dec, epoch=self.body._epoch)
+            ra, dec = ephem.Equatorial(original_radec, epoch=epoch_to_ephem(epoch)).get()
+            l, b = ephem.Galactic(ephem.Equatorial(ra, dec)).get()
             if is_iterable(timestamp):
                 return np.tile(l, len(timestamp)), np.tile(b, len(timestamp))
             else:
                 return l, b
-        ra, dec = self.astrometric_radec(timestamp, antenna)
+        ra, dec = self.astrometric_radec(timestamp, antenna, epoch=epoch)
         if is_iterable(ra):
             lb = np.array([ephem.Galactic(ephem.Equatorial(ra[n], dec[n])).get() for n in xrange(len(ra))])
             return lb[:, 0], lb[:, 1]
@@ -851,13 +862,15 @@ def construct_target_params(description):
             body.name = preferred_name
         else:
             body.name = "Ra: %s Dec: %s" % (ra, dec)
-        # Extract epoch info from tags
-        if ('B1900' in tags) or ('b1900' in tags):
-            body._epoch = ephem.B1900
-        elif ('B1950' in tags) or ('b1950' in tags):
-            body._epoch = ephem.B1950
-        else:
-            body._epoch = ephem.J2000
+        # Extract epoch info from tags, default to ephem.J2000
+        epoch = ephem.J2000
+        # Search for a well formed epoch string in tags
+        for tag in tags:
+            try:
+                epoch=epoch_to_ephem(tag)
+            except ValueError:
+                pass
+            body._epoch = epoch
         body._ra = ra
         body._dec = dec
 
@@ -872,7 +885,15 @@ def construct_target_params(description):
             body.name = preferred_name
         else:
             body.name = "Galactic l: %.4f b: %.4f" % (l, b)
-        body._epoch = ephem.J2000
+        # Extract epoch info from tags, default to ephem.J2000
+        epoch = ephem.J2000
+        # Search for a well formed epoch string in tags
+        for tag in tags:
+            try:
+                epoch=epoch_to_ephem(tag)
+            except ValueError:
+                pass
+            body._epoch = epoch
         body._ra = ra
         body._dec = dec
 
@@ -972,11 +993,11 @@ def construct_azel_target(az, el):
 #--- FUNCTION :  construct_radec_target
 #--------------------------------------------------------------------------------------------------
 
-def construct_radec_target(ra, dec):
+def construct_radec_target(ra, dec, epoch=ephem.J2000):
     """Convenience function to create unnamed fixed target (*radec* body type).
 
     The input parameters will also accept :class:`ephem.Angle` objects, as these
-    are floats in radians internally. The epoch is assumed to be J2000.
+    are floats in radians internally.
 
     Parameters
     ----------
@@ -986,6 +1007,8 @@ def construct_radec_target(ra, dec):
     dec : string or float
         Declination, either in 'D:M:S' or decimal degree string format, or as
         a float in radians
+    epoch : string or float or tuple or :class:`ephem.Date` object
+        The epoch of the position, in any format accepted by a call to ephem.Date() 
 
     Returns
     -------
@@ -1002,7 +1025,7 @@ def construct_radec_target(ra, dec):
             pass
     ra, dec = ephem.hours(ra), ephem.degrees(dec)
     body.name = "Ra: %s Dec: %s" % (ra, dec)
-    body._epoch = ephem.J2000
+    body._epoch = epoch
     body._ra = ra
     body._dec = dec
     return Target(body, 'radec')

katpoint/test/test_target.py

@@ -178,3 +178,14 @@ def test_projection(self):
         re_az, re_el = self.target.plane_to_sphere(x, y, self.ts, self.ant1)
         np.testing.assert_almost_equal(re_az, az, decimal=12)
         np.testing.assert_almost_equal(re_el, el, decimal=12)
+
+    def test_epoch(self):
+        """Test epoch conversion."""
+        targ1 = katpoint.Target('Belle Epoch, radec J1900, 18:00:00, -60:00:00',antenna=self.ant1)
+        ra1,dec1 = targ1.astrometric_radec(epoch='J1900')
+        prec_ra1,prec_dec1 = targ1.astrometric_radec(epoch='J1950')
+        targ2 = katpoint.construct_radec_target(prec_ra1,prec_dec1,epoch=katpoint.epoch_to_ephem('J1950'))
+        ra2,dec2 = targ2.astrometric_radec(epoch='J1900')
+        #Check to 9 decimal places should be accurate to better than 100th of an arcsecond at declination 0
+        np.testing.assert_almost_equal(ra1,ra2, decimal=9)
+        np.testing.assert_almost_equal(dec1,dec2, decimal=9)

katpoint/timestamp.py

@@ -4,6 +4,7 @@
 import math
 
 import numpy as np
+import re
 import ephem
 
 
@@ -156,3 +157,25 @@ def to_mjd(self):
         djd = self.to_ephem_date()
         return djd + 2415020 - 2400000.5
 
+def epoch_to_ephem(input_epoch='J2000'):
+    """Convert an epoch string to ephem date object.
+    An epoch can be specified as either B1900 or B1950 to specify the inbuilt Besselian
+    dates in ephem. Or can be of the form 'Jxxxx.xx' which specifies the date in terms
+    of Julian years from J2000 epoch.
+    """
+    # Extract epoch info from tags, default to ephem.J2000
+    if input_epoch in ['B1900','b1900','B1950','b1950','J2000','j2000']:
+        epoch = getattr(ephem, input_epoch.upper())
+    else:
+        #Search for a 'J' epoch string in tags
+        epoch_re = re.compile('^[Jj](\d+(?:\.\d+)?)$')
+        epoch_match = epoch_re.match(input_epoch)
+        if epoch_match:
+            #Convert first epoch string found to Dublin Julian date for input to Ephem
+            #Jxxxx.xx dates are defined (by convention) in Julian years from Gregorian J2000 epoch
+            epoch=(float(epoch_match.groups(1)[0])-2000.0)*365.25 + ephem.J2000
+        else:
+            raise ValueError("Epoch string '%s' not in correct format - " % (input_epoch,) +
+                             "should be Jxxxx.xx, where xxxx.xx specifies the date in terms of " +
+                             "Julian years from J2000 epoch.")
+    return ephem.Date(epoch)