Repository navigation
Add Heliocentric and Barycentric frames to coordinates #4288
Description
With reference to #2244. Specifically addressing @mhvk's point 1 in #2244 (comment). It would be very useful if we could add the above frames to the coordinates package. I don't think this should be too hard - I've included a first stab below.
The transformation framework allows you to calculate the location of the observatory in GCRS coordinates (see below). @eteq's transformation agrees with the code I posted in #2244 (comment) to within 70m, so that's a tiny part of #4268 ticked off.
from astropy.coordinates import sites
import astropy.coordinates as coord
from astropy.time import Time
from astropy import units
site_list = sites.get_builtin_sites()
WHT = site_list['lapalma']
now = Time.now()
itrs = coord.ITRS(WHT.itrs.cartesian,obstime=now)
gcrs = itrs.transform_to(coord.GCRS(obstime=now))
The position and velocity of the heliocenter and barycenter can be found using erfa.epv00, which is is consistent with the JPL ephemeris within 11.2 km = 4micro-sec and 5 mm/s, so the addition of a Heliocentric and Barycentric frame to the coordinates package doesn't look hard, e.g
def gcrs_to_bcrs(gcrs_loc):
tdb = gcrs_loc.obstime.tdb
#returns BCRS position and velocity of Earth in AU and AU/d
h_pv,b_pv = erfa.epv00(tdb.jd1,tdb.jd2)
# assume scalar obstime! and unpack position array
b_pos = b_pv[0,:] * units.AU
new_rep = CartesianRepresentation(gcrs.cartesian.xyz + b_pos)
return BCRS(new_rep,obstime=gcrs_loc.obstime)
One snag with the example above is that it assumes the obstime of the frame is a scalar.
Activity
mhvk commentedon Nov 4, 2015
What makes it fail for vectors? Since erfa is supposed to deal well with vectors, this might be a bug...
StuartLittlefair commentedon Nov 4, 2015
It's more that I don't know how to properly handle the coordinates. I can deal with the shape of the position-velocity vector returned by
erfa.epv00. I think the code below works, but it's different to how it's implemented in theicrs_to_helioecliptictransform.mhvk commentedon Nov 4, 2015
@StuartLittlefair - I think it should be possible to write it such that you don't have to worry about scalar versus vector at all. In particular, with #4123, it will be possible for coordinates to have more than 1 dimension, so we may as well future-proof it... To get the last axis to be the first one, one would need to do
b_pos = np.rollaxis(b_pos, -1)instead of doingb_pos.T(this works for scalars too, so the whole if statement could go).StuartLittlefair commentedon Nov 4, 2015
There ought to be a like button for comments @mhvk. I didn't know about
np.rollaxis.mhvk commentedon Nov 4, 2015
Nor did I until last night...
embray commentedon Nov 16, 2015
Indeed, rollaxis is often what you want over transpose for more than 2 dimensions.