diff --git a/katpoint/__init__.py b/katpoint/__init__.py index 58f1093..6502919 100644 --- a/katpoint/__init__.py +++ b/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 diff --git a/katpoint/target.py b/katpoint/target.py index 31fa3a2..0ae9cd9 100644 --- a/katpoint/target.py +++ b/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,12 +361,11 @@ 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 ---------- @@ -374,6 +373,10 @@ def astrometric_radec(self, timestamp=None, antenna=None): 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,12 +415,12 @@ 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 ---------- @@ -424,6 +428,10 @@ def galactic(self, timestamp=None, antenna=None): 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') diff --git a/katpoint/test/test_target.py b/katpoint/test/test_target.py index 3a2818a..65f9563 100644 --- a/katpoint/test/test_target.py +++ b/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) diff --git a/katpoint/timestamp.py b/katpoint/timestamp.py index da288d0..03cea3f 100644 --- a/katpoint/timestamp.py +++ b/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)