of Jupiter’s axis in the ecliptic reference frame: $$\label{kaxis}
\begin{array}{rcl}
\hat{k}_x&=&-0.01460\; ,\\
\noalign{\smallskip}
\hat{k}_y&=&-0.03582\; ,\\
\noalign{\smallskip}
\hat{k}_z&=&0.99925\; .
\end{array}$$ Finally, from the vector in Eq. (\[kaxis\]) and the spacecraft’s position vector we can determine the colatitude angle in Eqs. (\[Fgeor\]) and (\[Fgeot\]) in order to compute the force arising from the zonal harmonics.

We will also point out that we have used an iterative procedure to solve the equations of motion in Eq. (\[eqmotionr\])-(\[eqmotionv\]) because its simplicity and stability in comparison with the alternative approach in which the unknown functions, ${\bf r}(t)$ and ${\bf v}(t)$, appear also in the perturbation terms. The algorithm proceeds as follows:

-   We select a given timestep (in minutes in the ephemeris for Juno [@Horizons]) as close to the perijove as possible. This would be our initial condition. The fact that it does not coincide with the perijove is not relevant for our purpose.

    The main reason for this backwards and forward integration procedure is to reduce the propagation of errors in the numerical method. Alternatively, we can start from an instant $180$ minutes before the perijove and integrate throughout the perijove to another instant $180$ minutes after the perijove but we have found that 