Conversion of velocity vectors between ICRF frames with different central objects

Question

Following the very nice answers to this question, I am trying to implement an appropriate procedure for change of central body during orbit propagation.

The very interesting thesis by Mana Vautier linked by @DavidHammen clears up many of the practical problems, including even MATLAB code demonstrating the procedure. Upon detailed reading of the thesis and the included code, I have once again realized that there are more difficulties to this than I originally thought.

In particular, and limiting the problem for now to the Earth-Moon case, I am now dealing with the modifications that must be introduced in the last state vector of the Earth-centered propagation section to convert it into the first state vector for the Moon-centered propagation section. Conversion of the position vector seems very straightforward, since JPL development ephemerides can be used to obtain the position of the Moon from Earth. Subtracting these coordinates from a position vector in an Earth-centered reference frame then results in the position vector in a Moon-centered reference frame.

However, I am struggling to identify the best procedure to convert the velocity in Earth-centered ICRF into Moon-centered ICRF. The MATLAB code in the linked thesis, in particular the fragment at the beginning of page 54, seems to make use of the average Earth-Moon distance (EMdist) and the angular velocity of Moon around Earth (omega) to calculate the velocity of Moon with respect to Earth at a given instant. A perfectly circular orbit of Moon around Earth is assumed, which is then subtracted from the velocity in Earth-centered frame to obtain velocity in Moon-centered frame (note that only X and Y components are calculated because, as stated earlier in the thesis, for simplification Earth, Moon and the spacecraft are assumed to be coplanar).

However, the procedure assumes a perfectly circular orbit of Moon around Earth, which is not the case in reality.

I am therefore wondering, how would such a conversion of velocity vector be performed accurately? (let's say, accurately enough to allow for orbit propagation to be performed with errors to real position in the order of meters-tens of meters) Would it be OK to assume the orbits (of Moon around Earth for the specific case of transitions between the spheres of influence of Earth and Moon, or of the planets around the Solar System barycenter for the cases of interplanetary missions) are perfectly Keplerian elliptical orbits? Or would I need to approximate the velocity of the involved central bodies somehow from the positions calculated from JPL DE?

Answer 1

Assuming Newtonian mechanics are in play, velocities, like displacements, are three dimensional vectors. (There is no reason to go full-bore general relativistic with regard to behaviors in the solar system.) What this assumption means is you need some mechanism to compute the displacement vector and the velocity vector between the Earth and the Moon.

I recommend you use JPL's SPICE package. The Chebyshev polynomial coefficients as functions of time used to calculate displacements are easily differentiable, and SPICE provides a function that does just that, spkgeo_c. In addition to the C/C++ language interface to which I linked, there also are interfaces with other languages such as python.