Repository navigation
Add support for Heliocentric and Barycentric Julian Dates #2244
Description
Activity
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, @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
Timeand, soon,Coordinates). So, I'd suggest to first finish coordinates, then use those for a skeletonEarthLocationclass (#1928), use that inTime, 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's plan sounds promising to me. The remaining question is if we want to change
Timeto use a utility function (withlatandlonkeywords until we haveEarthLocation) for TDB now so that people don't start depending on the currentTime.lonandTime.latapproach...?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 - that looks great! I'd very much like to include something along these lines in
Time. I just checkedepv00and 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...
@StuartLittlefair - just a small comment --
heightis already part ofEarthLocation, so you would not need to override__init__-- but perhaps this is just a leftover from earlier days?@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
Timeobject and possibly anEarthLocationand returns aTimeobject in an appropriate scale. The location attributes for this newTimeobject 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 - I'm good with doing this as utility functions like you have outlined, probably putting this into a new module within the
timesubpackage 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 - 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
EarthLocationinto 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_posmethod. - 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
Timemethod that takes aSkyCoordinstance as import and uses the functions in (1) and (2), returning a time or velocity offset (i.e., aQuantity), just like it is now possible to get a sidereal time. Similarly, or perhaps alternatively, one may want to haveSkyCoordmethods that takes aTimeas an input. - With all the above in place, we might want to consider having
Timesupport more general locations (not onlyEarthLocation). 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
erfaapproximations 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
coordinatesrather thantime. 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...- We need something to convert
@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.
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.