Seite - 86 - in Contributions to GRACE Gravity Field Recovery - Improvements in Dynamic Orbit Integration, Stochastic Modelling of the Antenna Offset Correction, and Co-Estimation of Satellite Orientations

7.2.2 Parametrization of Reference Motion Traditionally, an orbital ellipse in a spherically symmetric potential is parametrised using Kepler parameters ξ= [ a e I ω Ω M ]T . (7.2.30) Numerical tests have shown that this parametrisation for the reference trajectory is not sufficiently stable when using standard double precision arithmetic. This statement is supported by the results later presented in section 7.3.2. This instability could potentially be remedied by computing all parameters relating to the reference trajectory in quadruple precision arithmetic, and then converting the computed state at each epoch to double precision for the further steps. It is undesirable to port the complete orbit integration algorithm to quadruple precision arithmetic, as this is sure to lead to significant performance penalties, with expected slowdowns on the order of a factor of 5 to 10 (Bailey and Borwein, 2015). Another solution is to parametrise the reference motion with a more stable set of orbital elements. Here, the equinoctial elements ξe= [ a h k p q λ ]T (7.2.31) asgiveninBrouckeandCefola (1972)areanattractiveoption. In theseelements, a is the semi-major axis of the ellipse. The elements h and k define the eccentricity and perigee of the orbit. The elements p and q encode the inclination of the orbital plane and the position of the ascending node.λ is the classical mean longitude. The equinoctial elements are a non-canonical set of orbit elements, with the Poincare´ elements their canonical counterpart (Vallado and McClain, 2001). Danielson et al. (1995, Section 2.1) gives a concise but comprehensive introduction to their derivation and use. The transformation from an equinoctial state vector to a Cartesian state vector can mostly be performed without relying on the evaluation of trigonometric functions, making this transformation very numerically stable. In this work all computations relating to the equinoctial elements are performed in double precision arithmetic, with the exception of the computation ofλ, the mean longitude. This is the fast-moving variable defining the position of the satellite along the equinoctial orbit arc.λ is computed and stored in quadruple precision. Asλ is not used in any expensive operations, the impact on overall performance is negligible. With this parametrization and implementation, the conversion from equinoctial elements to Cartesian coordinates shows sufficiently high stability. The results presented in sections 7.3.3 and 7.3.4 illustrate this statement. 7.2.3 Determination of Best-Fit Orbit To determine the best-fit Kepler ellipse from eq. (7.2.27), the partial derivatives of the position and velocity of the satellite in the CRF w.r.t. the equinoctial elements are Chapter7 Numerical Optimization in Orbit Integration86