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?


1 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.

  • $\begingroup$ Great! I am in fact currently using already the Chebyshev coefficients to obtain positions, from which I calculate displacements. So by differentiating these polynomials, we can obtain directly the velocities of the bodies also from JPL DE? really interesting! $\endgroup$
    – Rafa
    Commented Mar 1, 2022 at 11:06
  • $\begingroup$ I guess chbint_c might be the algorithm I'm looking for? $\endgroup$
    – Rafa
    Commented Mar 1, 2022 at 11:23
  • $\begingroup$ Ok in fact it seems the Clenshaw algorithm is easily applicable for the n-th derivatives as well as explained in this answer . So it would seem that makes it very straightforward to get velocities and accelerations from JPLDE! $\endgroup$
    – Rafa
    Commented Mar 1, 2022 at 11:41
  • 2
    $\begingroup$ The derivative calculation is fairly straightforward. That said, just use what SPICE already provides. I'm not fond of the 50+ year old API (pure ugliness), but it does work. $\endgroup$ Commented Mar 1, 2022 at 13:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.