Repository navigation
Add Heliocentric and Barycentric frames to coordinates #4288
Description
Activity
One snag with the example above is that it assumes the obstime of the frame is a scalar.
What makes it fail for vectors? Since erfa is supposed to deal well with vectors, this might be a bug...
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.def gcrs_to_bcrs(gcrs_loc): tdb = gcrs_loc.obstime.tdb #returns BCRS position and velocity of Earth in AU and AU/d pvh,pvb = erfa.epv00(tdb.jd1,tdb.jd2) # extract barycentric position from pos-vel vectors b_pos = pvb[..., 0, :] * units.AU # different paths for scalar and vector obstimes if b_pos.ndim == 1: #scalar obstime new_rep = CartesianRepresentation(gcrs.cartesian.xyz + b_pos) else: #vector obstime. b_pv has shape (N,2,3) new_rep = CartesianRepresentation(gcrs.cartesian.xyz + b_pos.T) return BCRS(new_rep,obstime=gcrs_loc.obstime)@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).There ought to be a like button for comments @mhvk. I didn't know about
np.rollaxis.Nor did I until last night...
Indeed, rollaxis is often what you want over transpose for more than 2 dimensions.
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.
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.gOne snag with the example above is that it assumes the obstime of the frame is a scalar.