Add support for Heliocentric and Barycentric Julian Dates

From a discussion on astropy-dev:
https://groups.google.com/forum/#!topic/astropy-dev/5KH5Ny1m6WY

The tentative conclusion there was to support Heliocentric and Barycentric Julian Dates via utility functions that accept a Time and Coordinate object.

@eteq

If we're willing to simply make this a utility function which accepts a Time and Coordinate and returns a corrected Time, then this should be relatively straightforward. I agree with @eteq that using some of the abstractions in coordinates should allow decoupling the Time and Coordinate. Basically BJD and HJD appear to be low-precision standards, so you should just be able to transform the coordinate to some reasonable frame (ICRS) and compute a correction.

I guess we need an Earth ephemeris though.

Yep, Earth location is definitely necessary here. But that's also going to be necessary at least in part for implementing the IAU2000 transformation stack, so it might be possible to use the ERFA pieces to do that. Even if not then, we'll eventually want a more complete interface to solar system ephemerides anyway.

But I'm tagging this as "future" for now, and we can re-tag it something firm once we've got pre-requisites in place.

@taldcroft

@taldcroft, @eteq - just saw this and the astropy-dev. I like the idea of providing utility functions since that doesn't commit us to an API for the classes yet (which needs some thought). I'm slightly wary of the "low precision" aspect, though, I think we should at least do (something like) the Stumpff (1980, A&AS 41, 1 [1]) algorithm, which also includes barycentric velocity. Ideally, of course, we do better and use the JPL ephemeris... It should be possible to lift at least the time aspect out of PINT [2] and/or jplephem [3](cc @scottransom, @demorest).

Overall, my sense would be for now to refer people to jplephem (or some lower precision equivalent) and for astropy aim to do it right from the get go, at full precision (as for Time and, soon, Coordinates). So, I'd suggest to first finish coordinates, then use those for a skeleton EarthLocation class (#1928), use that in Time, and then implement the proper JPL ephemeris.

[1] http://adsabs.harvard.edu/abs/1980A%26AS...41....1S
[2] https://github.com/nanograv/PINT
[3] https://github.com/brandon-rhodes/python-jplephem/

@mhvk

@mhvk's plan sounds promising to me. The remaining question is if we want to change Time to use a utility function (with lat and lon keywords until we have EarthLocation) for TDB now so that people don't start depending on the current Time.lon and Time.lat approach...?

I rated this "package-intermediate" if we go the utility function route, but it might be "package-expert" if it ends up being better done as a method

Probably, a function initially is best. I fear it is more than an intermediate amount of work, though, at least for the DE450 type of precision.

I don't know what the status of this issue is, or what the best way of contributing code to astropy is. However, I had to do this for my own work, and it was relatively easy following the work that has been done on erfa. I proceeded by sub-classing time.Time as shown in the code below to add helper functions that accept an EarthLocation object. In my tests of random points of the sky, it seems to agree with TEMPO2 answers to a few tens of microseconds.

I imagine it will need much work to include in the astropy project, but I thought I'd paste it here so you can use it if you wish

from astropy import time
from astropy import constants as const
from astropy import units as u
from astropy.utils.iers import IERS
from astropy.utils.iers import IERS_A, IERS_A_URL
from astropy.utils.data import download_file
from astropy.io import ascii
from astropy import coordinates as coord
from astropy import erfa
import numpy as np
import warnings

# mean sidereal rate (at J2000) in radians per (UT1) second
SR = 7.292115855306589e-5

''' time class: inherits astropy time object, and adds heliocentric, barycentric
    correction utilities.'''
class Time(time.Time):

    def __init__(self,*args,**kwargs):
        super(Time, self).__init__(*args, **kwargs)
        self.height = kwargs.get('height',0.0)

    def _pvobs(self):
        '''calculates position and velocity of the observatory
           returns position/velocity in AU and AU/d in GCRS reference frame
        '''

        # convert obs position from WGS84 (lat long) to ITRF geocentric coords in AU
        xyz = self.location.to(u.AU).value

        # now we need to convert this position to Celestial Coords
        # specifically, the GCRS coords.
        # conversion from celestial to terrestrial coords given by
        # [TRS] = RPOM * R_3(ERA) * RC2I * [CRS]
        # where:
        # [CRS] is vector in GCRS (geocentric celestial system)
        # [TRS] is vector in ITRS (International Terrestrial Ref System)
        # ERA is earth rotation angle
        # RPOM = polar motion matrix

        tt = self.tt
        mjd = self.utc.mjd

        # we need the IERS values to correct for the precession/nutation of the Earth
        iers_tab = IERS.open()

        # Find UT1, which is needed to calculate ERA
        # uses IERS_B by default , for more recent times use IERS_A download
        try:      
            ut1 = self.ut1 
        except:
            try:
                iers_a_file = download_file(IERS_A_URL, cache=True)
                iers_a = IERS_A.open(iers_a_file)
                self.delta_ut1_utc = self.get_delta_ut1_utc(iers_a)
                ut1 = self.ut1
            except:
                # fall back to UTC with degraded accuracy
                warnings.warn('Cannot calculate UT1: using UTC with degraded accuracy') 
                ut1 = self.utc

        # Gets x,y coords of Celestial Intermediate Pole (CIP) and CIO locator s
        # CIO = Celestial Intermediate Origin
        # Both in GCRS
        X,Y,S = erfa.xys00a(tt.jd1,tt.jd2)

        # Get dX and dY from IERS B
        dX = np.interp(mjd, iers_tab['MJD'], iers_tab['dX_2000A']) * u.arcsec 
        dY = np.interp(mjd, iers_tab['MJD'], iers_tab['dY_2000A']) * u.arcsec

        # Get GCRS to CIRS matrix
        # can be used to convert to Celestial Intermediate Ref Sys
        # from GCRS.
        rc2i = erfa.c2ixys(X+dX.to(u.rad).value, Y+dY.to(u.rad).value, S)

        # Gets the Terrestrial Intermediate Origin (TIO) locator s'
        # Terrestrial Intermediate Ref Sys (TIRS) defined by TIO and CIP.
        # TIRS related to to CIRS by Earth Rotation Angle
        sp = erfa.sp00(tt.jd1,tt.jd2)

        # Get X and Y from IERS B
        # X and Y are
        xp = np.interp(mjd, iers_tab['MJD'], iers_tab['PM_x']) * u.arcsec
        yp = np.interp(mjd, iers_tab['MJD'], iers_tab['PM_y']) * u.arcsec 

        # Get the polar motion matrix. Relates ITRF to TIRS.
        rpm = erfa.pom00(xp.to(u.rad).value, yp.to(u.rad).value, sp)

        # multiply ITRF position of obs by transpose of polar motion matrix
        # Gives Intermediate Ref Frame position of obs
        x,y,z = np.array([rpmMat.T.dot(xyz) for rpmMat in rpm]).T

        # Functions of Earth Rotation Angle, theta
        # Theta is angle bewtween TIO and CIO (along CIP)
        # USE UT1 here.
        theta = erfa.era00(ut1.jd1,ut1.jd2)
        S,C = np.sin(theta),np.cos(theta)

        # Position #GOT HERE
        pos = np.asarray([C*x - S*y, S*x + C*y, z]).T

        # multiply by inverse of GCRS to CIRS matrix
        # different methods for scalar times vs arrays
        if pos.ndim > 1:
            pos = np.array([np.dot(rc2i[j].T,pos[j]) for j in range(len(pos))])
        else:   
            pos = np.dot(rc2i.T,pos)

        # Velocity
        vel = np.asarray([SR*(-S*x - C*y), SR*(C*x-S*y), np.zeros_like(x)]).T
        # multiply by inverse of GCRS to CIRS matrix
        if vel.ndim > 1:
            vel = np.array([np.dot(rc2i[j].T,vel[j]) for j in range(len(pos))])
        else:        
            vel = np.dot(rc2i.T,vel)

        #return position and velocity
        return pos,vel

    def _obs_pos(self):
        '''calculates heliocentric and barycentric position of the earth in AU and AU/d'''
        tdb = self.tdb

        # get heliocentric and barycentric position and velocity of Earth
        # BCRS reference frame
        h_pv,b_pv = erfa.epv00(tdb.jd1,tdb.jd2)

        # h_pv etc can be shape (ntimes,2,3) or (2,3) if given a scalar time
        if h_pv.ndim == 2:
            h_pv = h_pv[np.newaxis,:]
        if b_pv.ndim == 2:
            b_pv = b_pv[np.newaxis,:]

        # unpack into position and velocity arrays
        h_pos = h_pv[:,0,:]
        h_vel = h_pv[:,1,:]

        # unpack into position and velocity arrays
        b_pos = b_pv[:,0,:]
        b_vel = b_pv[:,1,:]

        #now need position and velocity of observing station
        pos_obs, vel_obs = self._pvobs()

        #add this to heliocentric and barycentric position of center of Earth
        h_pos += pos_obs
        b_pos += pos_obs
        h_vel += vel_obs        
        b_vel += vel_obs        
        return (h_pos,h_vel,b_pos,b_vel)

    def _vect(self,coord):
        '''get unit vector pointing to star, and modulus of vector, in AU
           coordinate of star supplied as astropy.coordinate object

           assume zero proper motion, parallax and radial velocity'''
        pmra = pmdec = px = rv = 0.0

        rar  = coord.ra.radian
        decr = coord.dec.radian
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            # ignore warnings about 0 parallax
            pos,vel = erfa.starpv(rar,decr,pmra,pmdec,px,rv)

        modulus = np.sqrt(pos.dot(pos))
        unit = pos/modulus
        modulus /= const.au.value
        return modulus, unit

    def hcor(self,coord):
        mod, spos = self._vect(coord)
        # get helio/bary-centric position and velocity of observatory, in AU, AU/d
        h_pos,h_vel,b_pos,b_vel = self._obs_pos()

        # heliocentric light travel time, s
        tcor_hel = const.au.value * np.array([np.dot(spos,hpos) for hpos in h_pos]) / const.c.value
        #print 'Correction to add to get time at heliocentre = %.7f s' % tcor_hel
        dt = time.TimeDelta(tcor_hel, format='sec', scale='tdb')
        return self.utc + dt

    def bcor(self,coord):
        mod, spos = self._vect(coord)
        # get helio/bary-centric position and velocity of observatory, in AU, AU/d
        h_pos,h_vel,b_pos,b_vel  = self._obs_pos()

        # barycentric light travel time, s
        tcor_bar = const.au.value *  np.array([np.dot(spos,bpos) for bpos in b_pos]) / const.c.value
        #print 'Correction to add to get time at barycentre  = %.7f s' % tcor_bar
        dt = time.TimeDelta(tcor_bar, format='sec', scale='tdb')
        return self.tdb + dt

@StuartLittlefair

@StuartLittlefair - that looks great! I'd very much like to include something along these lines in Time. I just checked epv00 and it is consistent with the JPL ephemeris within 11.2 km = 4micro-sec and 5 mm/s. I think this suffices for most applications, including RV planet work, and for those for which it is not good enough (e.g., pulsar timing), we can later include the option to work with the JPL ephemerides (without having those huge files be part of astropy!).

@eteq - how far along are you for the transformation to the various intermediate schemes? Presumably, your work includes quite a bit of this...

@dpshelio

ping @dpshelio

@StuartLittlefair

@StuartLittlefair - just a small comment -- height is already part of EarthLocation, so you would not need to override __init__ -- but perhaps this is just a leftover from earlier days?

@mhvk

@mhvk indeed, this is just a hangover. Thanks for pointing this out - the height property is unused anyway, so there's no need to override init. The code should probably have some assertions in there to raise errors when the time object does not have a location property though.

Summary:

High priority

Has anyone yet formed a strong opinion about the need to develop a coherent framework for this? The problem is that a barycentric/heliocentric corrected time really has a position attached so it might be considered ideal implement this together with a framework that supports velocities, proper motions etc.

However, it's also an option to add a utility function to the time package that simply accepts a Time object and possibly an EarthLocation and returns a Time object in an appropriate scale. The location attributes for this new Time object would be nonsense of course.

If people are happy with the latter approach I can have a go at implementing this. It will be my first contribution to a big collaborative project though, so go easy on me...

@StuartLittlefair

@StuartLittlefair - I'm good with doing this as utility functions like you have outlined, probably putting this into a new module within the time subpackage and exporting the new public function(s). If it ever comes to pass that a grander framework is developed that includes this functionality then we can slowly deprecate the functional interface.

@StuartLittlefair

@StuartLittlefair - Help would be great! I think the first thing would be to split the problem up into pieces -- as your sample code above does already:

We need something to convert EarthLocation into a position relative to the solar system barycentre or the Sun. I think this may already be possible within coordinates, since the same would seem to be needed to support heliocentric ecliptic coordinates. @eteq? Otherwise, this might for now be a utility function that replaces your _obs_pos method.
Same for unit vector to object. My sense would be that coordinates supports this already too... @eteq?
The remainder is quite trivial, of course, and easiest may be to have a Time method that takes a SkyCoord instance as import and uses the functions in (1) and (2), returning a time or velocity offset (i.e., a Quantity), just like it is now possible to get a sidereal time. Similarly, or perhaps alternatively, one may want to have SkyCoord methods that takes a Time as an input.
With all the above in place, we might want to consider having Time support more general locations (not only EarthLocation). Since the corrections depend on the source one looked at, this doesn't seem trivial, and probably is of much lower priority.

The nice aspect of the separation that you already have, is that it should be relatively straightforward to substitute erfa approximations in (1) and (2) for functions that use the full JPL ephemerides.

Anyway, in the end I think @taldcroft is right in that the best way to start is with some utility functions, though arguably best within coordinates rather than time. Maybe @eteq can point us to what exists already -- by prodding around in the transformations, we may at the same time be able to help address #4268...

@mhvk

@mhvk - I agree that looks like a very sensible way forward.

With a few changes to the coordinates package I think the point 1 is quite straightforward. I'm not sure about 2, but as you say from that point on it is quite trivial, and just a question of API design.

Navigation Menu

Repository navigation

Add support for Heliocentric and Barycentric Julian Dates #2244

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Add support for Heliocentric and Barycentric Julian Dates #2244

Description

Activity

taldcroft commented on Mar 27, 2014

eteq commented on Mar 27, 2014

mhvk commented on Mar 28, 2014

eteq commented on Mar 28, 2014

eteq commented on Jul 11, 2014

mhvk commented on Jul 11, 2014

StuartLittlefair commented on Nov 28, 2014

mhvk commented on Nov 28, 2014

Cadair commented on Nov 28, 2014

mhvk commented on Nov 28, 2014

StuartLittlefair commented on Nov 28, 2014

taldcroft commented on Nov 2, 2015

StuartLittlefair commented on Nov 3, 2015

taldcroft commented on Nov 3, 2015

mhvk commented on Nov 3, 2015

StuartLittlefair commented on Nov 3, 2015

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions