1   /* Copyright 2002-2024 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.orbits;
18  
19  import java.io.Serializable;
20  
21  import org.hipparchus.analysis.differentiation.UnivariateDerivative1;
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.util.FastMath;
24  import org.hipparchus.util.SinCos;
25  import org.orekit.annotation.DefaultDataContext;
26  import org.orekit.data.DataContext;
27  import org.orekit.errors.OrekitIllegalArgumentException;
28  import org.orekit.errors.OrekitInternalError;
29  import org.orekit.errors.OrekitMessages;
30  import org.orekit.frames.Frame;
31  import org.orekit.time.AbsoluteDate;
32  import org.orekit.utils.PVCoordinates;
33  import org.orekit.utils.TimeStampedPVCoordinates;
34  
35  
36  /**
37   * This class handles circular orbital parameters.
38  
39   * <p>
40   * The parameters used internally are the circular elements which can be
41   * related to Keplerian elements as follows:
42   *   <ul>
43   *     <li>a</li>
44   *     <li>e<sub>x</sub> = e cos(ω)</li>
45   *     <li>e<sub>y</sub> = e sin(ω)</li>
46   *     <li>i</li>
47   *     <li>Ω</li>
48   *     <li>α<sub>v</sub> = v + ω</li>
49   *   </ul>
50   * where Ω stands for the Right Ascension of the Ascending Node and
51   * α<sub>v</sub> stands for the true latitude argument
52   *
53   * <p>
54   * The conversion equations from and to Keplerian elements given above hold only
55   * when both sides are unambiguously defined, i.e. when orbit is neither equatorial
56   * nor circular. When orbit is circular (but not equatorial), the circular
57   * parameters are still unambiguously defined whereas some Keplerian elements
58   * (more precisely ω and Ω) become ambiguous. When orbit is equatorial,
59   * neither the Keplerian nor the circular parameters can be defined unambiguously.
60   * {@link EquinoctialOrbit equinoctial orbits} is the recommended way to represent
61   * orbits.
62   * </p>
63   * <p>
64   * The instance <code>CircularOrbit</code> is guaranteed to be immutable.
65   * </p>
66   * @see    Orbit
67   * @see    KeplerianOrbit
68   * @see    CartesianOrbit
69   * @see    EquinoctialOrbit
70   * @author Luc Maisonobe
71   * @author Fabien Maussion
72   * @author V&eacute;ronique Pommier-Maurussane
73   */
74  
75  public class CircularOrbit extends Orbit implements PositionAngleBased {
76  
77      /** Serializable UID. */
78      private static final long serialVersionUID = 20231217L;
79  
80      /** Semi-major axis (m). */
81      private final double a;
82  
83      /** First component of the circular eccentricity vector. */
84      private final double ex;
85  
86      /** Second component of the circular eccentricity vector. */
87      private final double ey;
88  
89      /** Inclination (rad). */
90      private final double i;
91  
92      /** Right Ascension of Ascending Node (rad). */
93      private final double raan;
94  
95      /** Cached latitude argument (rad). */
96      private final double cachedAlpha;
97  
98      /** Type of cached position angle (latitude argument). */
99      private final PositionAngleType cachedPositionAngleType;
100 
101     /** Semi-major axis derivative (m/s). */
102     private final double aDot;
103 
104     /** First component of the circular eccentricity vector derivative. */
105     private final double exDot;
106 
107     /** Second component of the circular eccentricity vector derivative. */
108     private final double eyDot;
109 
110     /** Inclination derivative (rad/s). */
111     private final double iDot;
112 
113     /** Right Ascension of Ascending Node derivative (rad/s). */
114     private final double raanDot;
115 
116     /** True latitude argument derivative (rad/s). */
117     private final double cachedAlphaDot;
118 
119     /** Indicator for {@link PVCoordinates} serialization. */
120     private final boolean serializePV;
121 
122     /** Partial Cartesian coordinates (position and velocity are valid, acceleration may be missing). */
123     private transient PVCoordinates partialPV;
124 
125     /** Creates a new instance.
126      * @param a  semi-major axis (m)
127      * @param ex e cos(ω), first component of circular eccentricity vector
128      * @param ey e sin(ω), second component of circular eccentricity vector
129      * @param i inclination (rad)
130      * @param raan right ascension of ascending node (Ω, rad)
131      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
132      * @param type type of latitude argument
133      * @param cachedPositionAngleType type of cached latitude argument
134      * @param frame the frame in which are defined the parameters
135      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
136      * @param date date of the orbital parameters
137      * @param mu central attraction coefficient (m³/s²)
138      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
139      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
140      * @since 12.1
141      */
142     public CircularOrbit(final double a, final double ex, final double ey,
143                          final double i, final double raan, final double alpha,
144                          final PositionAngleType type, final PositionAngleType cachedPositionAngleType,
145                          final Frame frame, final AbsoluteDate date, final double mu)
146         throws IllegalArgumentException {
147         this(a, ex, ey, i, raan, alpha,
148              Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN,
149              type, cachedPositionAngleType, frame, date, mu);
150     }
151 
152     /** Creates a new instance without derivatives and with cached position angle same as value inputted.
153      * @param a  semi-major axis (m)
154      * @param ex e cos(ω), first component of circular eccentricity vector
155      * @param ey e sin(ω), second component of circular eccentricity vector
156      * @param i inclination (rad)
157      * @param raan right ascension of ascending node (Ω, rad)
158      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
159      * @param type type of latitude argument
160      * @param frame the frame in which are defined the parameters
161      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
162      * @param date date of the orbital parameters
163      * @param mu central attraction coefficient (m³/s²)
164      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
165      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
166      */
167     public CircularOrbit(final double a, final double ex, final double ey,
168                          final double i, final double raan, final double alpha,
169                          final PositionAngleType type,
170                          final Frame frame, final AbsoluteDate date, final double mu)
171             throws IllegalArgumentException {
172         this(a, ex, ey, i, raan, alpha, type, type, frame, date, mu);
173     }
174 
175     /** Creates a new instance.
176      * @param a  semi-major axis (m)
177      * @param ex e cos(ω), first component of circular eccentricity vector
178      * @param ey e sin(ω), second component of circular eccentricity vector
179      * @param i inclination (rad)
180      * @param raan right ascension of ascending node (Ω, rad)
181      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
182      * @param aDot  semi-major axis derivative (m/s)
183      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
184      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
185      * @param iDot inclination  derivative(rad/s)
186      * @param raanDot right ascension of ascending node derivative (rad/s)
187      * @param alphaDot  d(an + ω), mean, eccentric or true latitude argument derivative (rad/s)
188      * @param type type of latitude argument
189      * @param cachedPositionAngleType type of cached latitude argument
190      * @param frame the frame in which are defined the parameters
191      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
192      * @param date date of the orbital parameters
193      * @param mu central attraction coefficient (m³/s²)
194      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
195      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
196      * @since 12.1
197      */
198     public CircularOrbit(final double a, final double ex, final double ey,
199                          final double i, final double raan, final double alpha,
200                          final double aDot, final double exDot, final double eyDot,
201                          final double iDot, final double raanDot, final double alphaDot,
202                          final PositionAngleType type, final PositionAngleType cachedPositionAngleType,
203                          final Frame frame, final AbsoluteDate date, final double mu)
204         throws IllegalArgumentException {
205         super(frame, date, mu);
206         if (ex * ex + ey * ey >= 1.0) {
207             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
208                                                      getClass().getName());
209         }
210         this.a       =  a;
211         this.aDot    =  aDot;
212         this.ex      = ex;
213         this.exDot   = exDot;
214         this.ey      = ey;
215         this.eyDot   = eyDot;
216         this.i       = i;
217         this.iDot    = iDot;
218         this.raan    = raan;
219         this.raanDot = raanDot;
220         this.cachedPositionAngleType = cachedPositionAngleType;
221 
222         if (hasDerivatives()) {
223             final UnivariateDerivative1 alphaUD = initializeCachedAlpha(alpha, alphaDot, type);
224             this.cachedAlpha = alphaUD.getValue();
225             this.cachedAlphaDot = alphaUD.getFirstDerivative();
226         } else {
227             this.cachedAlpha = initializeCachedAlpha(alpha, type);
228             this.cachedAlphaDot = Double.NaN;
229         }
230 
231         serializePV = false;
232         partialPV   = null;
233 
234     }
235 
236     /** Creates a new instance with derivatives and with cached position angle same as value inputted.
237      * @param a  semi-major axis (m)
238      * @param ex e cos(ω), first component of circular eccentricity vector
239      * @param ey e sin(ω), second component of circular eccentricity vector
240      * @param i inclination (rad)
241      * @param raan right ascension of ascending node (Ω, rad)
242      * @param alpha  an + ω, mean, eccentric or true latitude argument (rad)
243      * @param aDot  semi-major axis derivative (m/s)
244      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
245      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
246      * @param iDot inclination  derivative(rad/s)
247      * @param raanDot right ascension of ascending node derivative (rad/s)
248      * @param alphaDot  d(an + ω), mean, eccentric or true latitude argument derivative (rad/s)
249      * @param type type of latitude argument
250      * @param frame the frame in which are defined the parameters
251      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
252      * @param date date of the orbital parameters
253      * @param mu central attraction coefficient (m³/s²)
254      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
255      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
256      */
257     public CircularOrbit(final double a, final double ex, final double ey,
258                          final double i, final double raan, final double alpha,
259                          final double aDot, final double exDot, final double eyDot,
260                          final double iDot, final double raanDot, final double alphaDot,
261                          final PositionAngleType type,
262                          final Frame frame, final AbsoluteDate date, final double mu)
263             throws IllegalArgumentException {
264         this(a, ex, ey, i, raan, alpha, aDot, exDot, eyDot, iDot, raanDot, alphaDot, type, type,
265                 frame, date, mu);
266     }
267 
268     /** Creates a new instance.
269      * @param a  semi-major axis (m)
270      * @param ex e cos(ω), first component of circular eccentricity vector
271      * @param ey e sin(ω), second component of circular eccentricity vector
272      * @param i inclination (rad)
273      * @param raan right ascension of ascending node (Ω, rad)
274      * @param alpha  input latitude argument (rad)
275      * @param aDot  semi-major axis derivative (m/s)
276      * @param exDot d(e cos(ω))/dt, first component of circular eccentricity vector derivative
277      * @param eyDot d(e sin(ω))/dt, second component of circular eccentricity vector derivative
278      * @param iDot inclination  derivative(rad/s)
279      * @param raanDot right ascension of ascending node derivative (rad/s)
280      * @param alphaDot  input latitude argument derivative (rad/s)
281      * @param pvCoordinates the {@link PVCoordinates} in inertial frame
282      * @param positionAngleType type of position angle
283      * @param frame the frame in which are defined the parameters
284      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
285      * @param mu central attraction coefficient (m³/s²)
286      * @exception IllegalArgumentException if eccentricity is equal to 1 or larger or
287      * if frame is not a {@link Frame#isPseudoInertial pseudo-inertial frame}
288      */
289     private CircularOrbit(final double a, final double ex, final double ey,
290                           final double i, final double raan, final double alpha,
291                           final double aDot, final double exDot, final double eyDot,
292                           final double iDot, final double raanDot, final double alphaDot,
293                           final TimeStampedPVCoordinates pvCoordinates,
294                           final PositionAngleType positionAngleType, final Frame frame, final double mu)
295         throws IllegalArgumentException {
296         super(pvCoordinates, frame, mu);
297         this.a           =  a;
298         this.aDot        =  aDot;
299         this.ex          = ex;
300         this.exDot       = exDot;
301         this.ey          = ey;
302         this.eyDot       = eyDot;
303         this.i           = i;
304         this.iDot        = iDot;
305         this.raan        = raan;
306         this.raanDot     = raanDot;
307         this.cachedAlpha = alpha;
308         this.cachedAlphaDot = alphaDot;
309         this.cachedPositionAngleType = positionAngleType;
310         this.serializePV = true;
311         this.partialPV   = null;
312     }
313 
314     /** Constructor from Cartesian parameters.
315      *
316      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
317      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
318      * use {@code mu} and the position to compute the acceleration, including
319      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
320      *
321      * @param pvCoordinates the {@link PVCoordinates} in inertial frame
322      * @param frame the frame in which are defined the {@link PVCoordinates}
323      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
324      * @param mu central attraction coefficient (m³/s²)
325      * @exception IllegalArgumentException if frame is not a {@link
326      * Frame#isPseudoInertial pseudo-inertial frame}
327      */
328     public CircularOrbit(final TimeStampedPVCoordinates pvCoordinates, final Frame frame, final double mu)
329         throws IllegalArgumentException {
330         super(pvCoordinates, frame, mu);
331         this.cachedPositionAngleType = PositionAngleType.TRUE;
332 
333         // compute semi-major axis
334         final Vector3D pvP = pvCoordinates.getPosition();
335         final Vector3D pvV = pvCoordinates.getVelocity();
336         final Vector3D pvA = pvCoordinates.getAcceleration();
337         final double r2 = pvP.getNormSq();
338         final double r  = FastMath.sqrt(r2);
339         final double V2 = pvV.getNormSq();
340         final double rV2OnMu = r * V2 / mu;
341         a = r / (2 - rV2OnMu);
342 
343         if (!isElliptical()) {
344             throw new OrekitIllegalArgumentException(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS,
345                                                      getClass().getName());
346         }
347 
348         // compute inclination
349         final Vector3D momentum = pvCoordinates.getMomentum();
350         i = Vector3D.angle(momentum, Vector3D.PLUS_K);
351 
352         // compute right ascension of ascending node
353         final Vector3D node  = Vector3D.crossProduct(Vector3D.PLUS_K, momentum);
354         raan = FastMath.atan2(node.getY(), node.getX());
355 
356         // 2D-coordinates in the canonical frame
357         final SinCos scRaan = FastMath.sinCos(raan);
358         final SinCos scI    = FastMath.sinCos(i);
359         final double xP     = pvP.getX();
360         final double yP     = pvP.getY();
361         final double zP     = pvP.getZ();
362         final double x2     = (xP * scRaan.cos() + yP * scRaan.sin()) / a;
363         final double y2     = ((yP * scRaan.cos() - xP * scRaan.sin()) * scI.cos() + zP * scI.sin()) / a;
364 
365         // compute eccentricity vector
366         final double eSE    = Vector3D.dotProduct(pvP, pvV) / FastMath.sqrt(mu * a);
367         final double eCE    = rV2OnMu - 1;
368         final double e2     = eCE * eCE + eSE * eSE;
369         final double f      = eCE - e2;
370         final double g      = FastMath.sqrt(1 - e2) * eSE;
371         final double aOnR   = a / r;
372         final double a2OnR2 = aOnR * aOnR;
373         ex = a2OnR2 * (f * x2 + g * y2);
374         ey = a2OnR2 * (f * y2 - g * x2);
375 
376         // compute latitude argument
377         final double beta = 1 / (1 + FastMath.sqrt(1 - ex * ex - ey * ey));
378         cachedAlpha = CircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, FastMath.atan2(y2 + ey + eSE * beta * ex, x2 + ex - eSE * beta * ey));
379 
380         partialPV   = pvCoordinates;
381 
382         if (hasNonKeplerianAcceleration(pvCoordinates, mu)) {
383             // we have a relevant acceleration, we can compute derivatives
384 
385             final double[][] jacobian = new double[6][6];
386             getJacobianWrtCartesian(PositionAngleType.MEAN, jacobian);
387 
388             final Vector3D keplerianAcceleration    = new Vector3D(-mu / (r * r2), pvP);
389             final Vector3D nonKeplerianAcceleration = pvA.subtract(keplerianAcceleration);
390             final double   aX                       = nonKeplerianAcceleration.getX();
391             final double   aY                       = nonKeplerianAcceleration.getY();
392             final double   aZ                       = nonKeplerianAcceleration.getZ();
393             aDot    = jacobian[0][3] * aX + jacobian[0][4] * aY + jacobian[0][5] * aZ;
394             exDot   = jacobian[1][3] * aX + jacobian[1][4] * aY + jacobian[1][5] * aZ;
395             eyDot   = jacobian[2][3] * aX + jacobian[2][4] * aY + jacobian[2][5] * aZ;
396             iDot    = jacobian[3][3] * aX + jacobian[3][4] * aY + jacobian[3][5] * aZ;
397             raanDot = jacobian[4][3] * aX + jacobian[4][4] * aY + jacobian[4][5] * aZ;
398 
399             // in order to compute latitude argument derivative, we must compute
400             // mean latitude argument derivative including Keplerian motion and convert to true latitude argument
401             final double alphaMDot = getKeplerianMeanMotion() +
402                                      jacobian[5][3] * aX + jacobian[5][4] * aY + jacobian[5][5] * aZ;
403             final UnivariateDerivative1 exUD     = new UnivariateDerivative1(ex, exDot);
404             final UnivariateDerivative1 eyUD     = new UnivariateDerivative1(ey, eyDot);
405             final UnivariateDerivative1 alphaMUD = new UnivariateDerivative1(getAlphaM(), alphaMDot);
406             final UnivariateDerivative1 alphavUD = FieldCircularLatitudeArgumentUtility.meanToTrue(exUD, eyUD, alphaMUD);
407             cachedAlphaDot = alphavUD.getFirstDerivative();
408 
409         } else {
410             // acceleration is either almost zero or NaN,
411             // we assume acceleration was not known
412             // we don't set up derivatives
413             aDot      = Double.NaN;
414             exDot     = Double.NaN;
415             eyDot     = Double.NaN;
416             iDot      = Double.NaN;
417             raanDot   = Double.NaN;
418             cachedAlphaDot = Double.NaN;
419         }
420 
421         serializePV = true;
422 
423     }
424 
425     /** Constructor from Cartesian parameters.
426      *
427      * <p> The acceleration provided in {@code pvCoordinates} is accessible using
428      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
429      * use {@code mu} and the position to compute the acceleration, including
430      * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
431      *
432      * @param pvCoordinates the {@link PVCoordinates} in inertial frame
433      * @param frame the frame in which are defined the {@link PVCoordinates}
434      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
435      * @param date date of the orbital parameters
436      * @param mu central attraction coefficient (m³/s²)
437      * @exception IllegalArgumentException if frame is not a {@link
438      * Frame#isPseudoInertial pseudo-inertial frame}
439      */
440     public CircularOrbit(final PVCoordinates pvCoordinates, final Frame frame,
441                          final AbsoluteDate date, final double mu)
442         throws IllegalArgumentException {
443         this(new TimeStampedPVCoordinates(date, pvCoordinates), frame, mu);
444     }
445 
446     /** Constructor from any kind of orbital parameters.
447      * @param op orbital parameters to copy
448      */
449     public CircularOrbit(final Orbit op) {
450 
451         super(op.getFrame(), op.getDate(), op.getMu());
452 
453         a    = op.getA();
454         i    = op.getI();
455         final double hx = op.getHx();
456         final double hy = op.getHy();
457         final double h2 = hx * hx + hy * hy;
458         final double h  = FastMath.sqrt(h2);
459         raan = FastMath.atan2(hy, hx);
460         final SinCos scRaan  = FastMath.sinCos(raan);
461         final double cosRaan = h == 0 ? scRaan.cos() : hx / h;
462         final double sinRaan = h == 0 ? scRaan.sin() : hy / h;
463         final double equiEx = op.getEquinoctialEx();
464         final double equiEy = op.getEquinoctialEy();
465         ex     = equiEx * cosRaan + equiEy * sinRaan;
466         ey     = equiEy * cosRaan - equiEx * sinRaan;
467         cachedPositionAngleType = PositionAngleType.TRUE;
468         cachedAlpha = op.getLv() - raan;
469 
470         if (op.hasDerivatives()) {
471             aDot    = op.getADot();
472             final double hxDot = op.getHxDot();
473             final double hyDot = op.getHyDot();
474             iDot    = 2 * (cosRaan * hxDot + sinRaan * hyDot) / (1 + h2);
475             raanDot = (hx * hyDot - hy * hxDot) / h2;
476             final double equiExDot = op.getEquinoctialExDot();
477             final double equiEyDot = op.getEquinoctialEyDot();
478             exDot   = (equiExDot + equiEy * raanDot) * cosRaan +
479                       (equiEyDot - equiEx * raanDot) * sinRaan;
480             eyDot   = (equiEyDot - equiEx * raanDot) * cosRaan -
481                       (equiExDot + equiEy * raanDot) * sinRaan;
482             cachedAlphaDot = op.getLvDot() - raanDot;
483         } else {
484             aDot      = Double.NaN;
485             exDot     = Double.NaN;
486             eyDot     = Double.NaN;
487             iDot      = Double.NaN;
488             raanDot   = Double.NaN;
489             cachedAlphaDot = Double.NaN;
490         }
491 
492         serializePV = false;
493         partialPV   = null;
494 
495     }
496 
497     /** {@inheritDoc} */
498     @Override
499     public OrbitType getType() {
500         return OrbitType.CIRCULAR;
501     }
502 
503     /** {@inheritDoc} */
504     @Override
505     public double getA() {
506         return a;
507     }
508 
509     /** {@inheritDoc} */
510     @Override
511     public double getADot() {
512         return aDot;
513     }
514 
515     /** {@inheritDoc} */
516     @Override
517     public double getEquinoctialEx() {
518         final SinCos sc = FastMath.sinCos(raan);
519         return ex * sc.cos() - ey * sc.sin();
520     }
521 
522     /** {@inheritDoc} */
523     @Override
524     public double getEquinoctialExDot() {
525         final SinCos sc = FastMath.sinCos(raan);
526         return (exDot - ey * raanDot) * sc.cos() - (eyDot + ex * raanDot) * sc.sin();
527     }
528 
529     /** {@inheritDoc} */
530     @Override
531     public double getEquinoctialEy() {
532         final SinCos sc = FastMath.sinCos(raan);
533         return ey * sc.cos() + ex * sc.sin();
534     }
535 
536     /** {@inheritDoc} */
537     @Override
538     public double getEquinoctialEyDot() {
539         final SinCos sc = FastMath.sinCos(raan);
540         return (eyDot + ex * raanDot) * sc.cos() + (exDot - ey * raanDot) * sc.sin();
541     }
542 
543     /** Get the first component of the circular eccentricity vector.
544      * @return ex = e cos(ω), first component of the circular eccentricity vector
545      */
546     public double getCircularEx() {
547         return ex;
548     }
549 
550     /** Get the first component of the circular eccentricity vector derivative.
551      * @return ex = e cos(ω), first component of the circular eccentricity vector derivative
552      * @since 9.0
553      */
554     public double getCircularExDot() {
555         return exDot;
556     }
557 
558     /** Get the second component of the circular eccentricity vector.
559      * @return ey = e sin(ω), second component of the circular eccentricity vector
560      */
561     public double getCircularEy() {
562         return ey;
563     }
564 
565     /** Get the second component of the circular eccentricity vector derivative.
566      * @return ey = e sin(ω), second component of the circular eccentricity vector derivative
567      */
568     public double getCircularEyDot() {
569         return eyDot;
570     }
571 
572     /** {@inheritDoc} */
573     @Override
574     public double getHx() {
575         // Check for equatorial retrograde orbit
576         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
577             return Double.NaN;
578         }
579         return  FastMath.cos(raan) * FastMath.tan(i / 2);
580     }
581 
582     /** {@inheritDoc} */
583     @Override
584     public double getHxDot() {
585         // Check for equatorial retrograde orbit
586         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
587             return Double.NaN;
588         }
589         final SinCos sc  = FastMath.sinCos(raan);
590         final double tan = FastMath.tan(0.5 * i);
591         return 0.5 * sc.cos() * (1 + tan * tan) * iDot - sc.sin() * tan * raanDot;
592     }
593 
594     /** {@inheritDoc} */
595     @Override
596     public double getHy() {
597         // Check for equatorial retrograde orbit
598         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
599             return Double.NaN;
600         }
601         return  FastMath.sin(raan) * FastMath.tan(i / 2);
602     }
603 
604     /** {@inheritDoc} */
605     @Override
606     public double getHyDot() {
607         // Check for equatorial retrograde orbit
608         if (FastMath.abs(i - FastMath.PI) < 1.0e-10) {
609             return Double.NaN;
610         }
611         final SinCos sc  = FastMath.sinCos(raan);
612         final double tan = FastMath.tan(0.5 * i);
613         return 0.5 * sc.sin() * (1 + tan * tan) * iDot + sc.cos() * tan * raanDot;
614     }
615 
616     /** Get the true latitude argument.
617      * @return v + ω true latitude argument (rad)
618      */
619     public double getAlphaV() {
620         switch (cachedPositionAngleType) {
621             case TRUE:
622                 return cachedAlpha;
623 
624             case ECCENTRIC:
625                 return CircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, cachedAlpha);
626 
627             case MEAN:
628                 return CircularLatitudeArgumentUtility.meanToTrue(ex, ey, cachedAlpha);
629 
630             default:
631                 throw new OrekitInternalError(null);
632         }
633     }
634 
635     /** Get the true latitude argument derivative.
636      * <p>
637      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
638      * </p>
639      * @return v + ω true latitude argument derivative (rad/s)
640      * @since 9.0
641      */
642     public double getAlphaVDot() {
643         switch (cachedPositionAngleType) {
644             case ECCENTRIC:
645                 final UnivariateDerivative1 alphaEUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
646                 final UnivariateDerivative1 exUD     = new UnivariateDerivative1(ex,     exDot);
647                 final UnivariateDerivative1 eyUD     = new UnivariateDerivative1(ey,     eyDot);
648                 final UnivariateDerivative1 alphaVUD = FieldCircularLatitudeArgumentUtility.eccentricToTrue(exUD, eyUD,
649                         alphaEUD);
650                 return alphaVUD.getFirstDerivative();
651 
652             case TRUE:
653                 return cachedAlphaDot;
654 
655             case MEAN:
656                 final UnivariateDerivative1 alphaMUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
657                 final UnivariateDerivative1 exUD2    = new UnivariateDerivative1(ex,     exDot);
658                 final UnivariateDerivative1 eyUD2    = new UnivariateDerivative1(ey,     eyDot);
659                 final UnivariateDerivative1 alphaVUD2 = FieldCircularLatitudeArgumentUtility.meanToTrue(exUD2,
660                         eyUD2, alphaMUD);
661                 return alphaVUD2.getFirstDerivative();
662 
663             default:
664                 throw new OrekitInternalError(null);
665         }
666     }
667 
668     /** Get the eccentric latitude argument.
669      * @return E + ω eccentric latitude argument (rad)
670      */
671     public double getAlphaE() {
672         switch (cachedPositionAngleType) {
673             case TRUE:
674                 return CircularLatitudeArgumentUtility.trueToEccentric(ex, ey, cachedAlpha);
675 
676             case ECCENTRIC:
677                 return cachedAlpha;
678 
679             case MEAN:
680                 return CircularLatitudeArgumentUtility.meanToEccentric(ex, ey, cachedAlpha);
681 
682             default:
683                 throw new OrekitInternalError(null);
684         }
685     }
686 
687     /** Get the eccentric latitude argument derivative.
688      * <p>
689      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
690      * </p>
691      * @return d(E + ω)/dt eccentric latitude argument derivative (rad/s)
692      * @since 9.0
693      */
694     public double getAlphaEDot() {
695         switch (cachedPositionAngleType) {
696             case TRUE:
697                 final UnivariateDerivative1 alphaVUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
698                 final UnivariateDerivative1 exUD     = new UnivariateDerivative1(ex,     exDot);
699                 final UnivariateDerivative1 eyUD     = new UnivariateDerivative1(ey,     eyDot);
700                 final UnivariateDerivative1 alphaEUD = FieldCircularLatitudeArgumentUtility.trueToEccentric(exUD, eyUD,
701                         alphaVUD);
702                 return alphaEUD.getFirstDerivative();
703 
704             case ECCENTRIC:
705                 return cachedAlphaDot;
706 
707             case MEAN:
708                 final UnivariateDerivative1 alphaMUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
709                 final UnivariateDerivative1 exUD2    = new UnivariateDerivative1(ex,     exDot);
710                 final UnivariateDerivative1 eyUD2    = new UnivariateDerivative1(ey,     eyDot);
711                 final UnivariateDerivative1 alphaVUD2 = FieldCircularLatitudeArgumentUtility.meanToEccentric(exUD2,
712                         eyUD2, alphaMUD);
713                 return alphaVUD2.getFirstDerivative();
714 
715             default:
716                 throw new OrekitInternalError(null);
717         }
718     }
719 
720     /** Get the mean latitude argument.
721      * @return M + ω mean latitude argument (rad)
722      */
723     public double getAlphaM() {
724         switch (cachedPositionAngleType) {
725             case TRUE:
726                 return CircularLatitudeArgumentUtility.trueToMean(ex, ey, cachedAlpha);
727 
728             case MEAN:
729                 return cachedAlpha;
730 
731             case ECCENTRIC:
732                 return CircularLatitudeArgumentUtility.eccentricToMean(ex, ey, cachedAlpha);
733 
734             default:
735                 throw new OrekitInternalError(null);
736         }
737     }
738 
739     /** Get the mean latitude argument derivative.
740      * <p>
741      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
742      * </p>
743      * @return d(M + ω)/dt mean latitude argument derivative (rad/s)
744      * @since 9.0
745      */
746     public double getAlphaMDot() {
747         switch (cachedPositionAngleType) {
748             case TRUE:
749                 final UnivariateDerivative1 alphaVUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
750                 final UnivariateDerivative1 exUD     = new UnivariateDerivative1(ex,     exDot);
751                 final UnivariateDerivative1 eyUD     = new UnivariateDerivative1(ey,     eyDot);
752                 final UnivariateDerivative1 alphaMUD = FieldCircularLatitudeArgumentUtility.trueToMean(exUD, eyUD,
753                         alphaVUD);
754                 return alphaMUD.getFirstDerivative();
755 
756             case MEAN:
757                 return cachedAlphaDot;
758 
759             case ECCENTRIC:
760                 final UnivariateDerivative1 alphaEUD = new UnivariateDerivative1(cachedAlpha, cachedAlphaDot);
761                 final UnivariateDerivative1 exUD2    = new UnivariateDerivative1(ex,     exDot);
762                 final UnivariateDerivative1 eyUD2    = new UnivariateDerivative1(ey,     eyDot);
763                 final UnivariateDerivative1 alphaMUD2 = FieldCircularLatitudeArgumentUtility.eccentricToMean(exUD2,
764                         eyUD2, alphaEUD);
765                 return alphaMUD2.getFirstDerivative();
766 
767             default:
768                 throw new OrekitInternalError(null);
769         }
770     }
771 
772     /** Get the latitude argument.
773      * @param type type of the angle
774      * @return latitude argument (rad)
775      */
776     public double getAlpha(final PositionAngleType type) {
777         return (type == PositionAngleType.MEAN) ? getAlphaM() :
778                                               ((type == PositionAngleType.ECCENTRIC) ? getAlphaE() :
779                                                                                    getAlphaV());
780     }
781 
782     /** Get the latitude argument derivative.
783      * <p>
784      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
785      * </p>
786      * @param type type of the angle
787      * @return latitude argument derivative (rad/s)
788      * @since 9.0
789      */
790     public double getAlphaDot(final PositionAngleType type) {
791         return (type == PositionAngleType.MEAN) ? getAlphaMDot() :
792                                               ((type == PositionAngleType.ECCENTRIC) ? getAlphaEDot() :
793                                                                                    getAlphaVDot());
794     }
795 
796     /** Computes the true latitude argument from the eccentric latitude argument.
797      * @param alphaE = E + ω eccentric latitude argument (rad)
798      * @param ex e cos(ω), first component of circular eccentricity vector
799      * @param ey e sin(ω), second component of circular eccentricity vector
800      * @return the true latitude argument.
801      */
802     @Deprecated
803     public static double eccentricToTrue(final double alphaE, final double ex, final double ey) {
804         return CircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, alphaE);
805     }
806 
807     /** Computes the eccentric latitude argument from the true latitude argument.
808      * @param alphaV = V + ω true latitude argument (rad)
809      * @param ex e cos(ω), first component of circular eccentricity vector
810      * @param ey e sin(ω), second component of circular eccentricity vector
811      * @return the eccentric latitude argument.
812      */
813     @Deprecated
814     public static double trueToEccentric(final double alphaV, final double ex, final double ey) {
815         return CircularLatitudeArgumentUtility.trueToEccentric(ex, ey, alphaV);
816     }
817 
818     /** Computes the eccentric latitude argument from the mean latitude argument.
819      * @param alphaM = M + ω  mean latitude argument (rad)
820      * @param ex e cos(ω), first component of circular eccentricity vector
821      * @param ey e sin(ω), second component of circular eccentricity vector
822      * @return the eccentric latitude argument.
823      */
824     @Deprecated
825     public static double meanToEccentric(final double alphaM, final double ex, final double ey) {
826         return CircularLatitudeArgumentUtility.meanToEccentric(ex, ey, alphaM);
827     }
828 
829     /** Computes the mean latitude argument from the eccentric latitude argument.
830      * @param alphaE = E + ω  mean latitude argument (rad)
831      * @param ex e cos(ω), first component of circular eccentricity vector
832      * @param ey e sin(ω), second component of circular eccentricity vector
833      * @return the mean latitude argument.
834      */
835     @Deprecated
836     public static double eccentricToMean(final double alphaE, final double ex, final double ey) {
837         return CircularLatitudeArgumentUtility.eccentricToMean(ex, ey, alphaE);
838     }
839 
840     /** {@inheritDoc} */
841     @Override
842     public double getE() {
843         return FastMath.sqrt(ex * ex + ey * ey);
844     }
845 
846     /** {@inheritDoc} */
847     @Override
848     public double getEDot() {
849         return (ex * exDot + ey * eyDot) / getE();
850     }
851 
852     /** {@inheritDoc} */
853     @Override
854     public double getI() {
855         return i;
856     }
857 
858     /** {@inheritDoc} */
859     @Override
860     public double getIDot() {
861         return iDot;
862     }
863 
864     /** Get the right ascension of the ascending node.
865      * @return right ascension of the ascending node (rad)
866      */
867     public double getRightAscensionOfAscendingNode() {
868         return raan;
869     }
870 
871     /** Get the right ascension of the ascending node derivative.
872      * <p>
873      * If the orbit was created without derivatives, the value returned is {@link Double#NaN}.
874      * </p>
875      * @return right ascension of the ascending node derivative (rad/s)
876      * @since 9.0
877      */
878     public double getRightAscensionOfAscendingNodeDot() {
879         return raanDot;
880     }
881 
882     /** {@inheritDoc} */
883     @Override
884     public double getLv() {
885         return getAlphaV() + raan;
886     }
887 
888     /** {@inheritDoc} */
889     @Override
890     public double getLvDot() {
891         return getAlphaVDot() + raanDot;
892     }
893 
894     /** {@inheritDoc} */
895     @Override
896     public double getLE() {
897         return getAlphaE() + raan;
898     }
899 
900     /** {@inheritDoc} */
901     @Override
902     public double getLEDot() {
903         return getAlphaEDot() + raanDot;
904     }
905 
906     /** {@inheritDoc} */
907     @Override
908     public double getLM() {
909         return getAlphaM() + raan;
910     }
911 
912     /** {@inheritDoc} */
913     @Override
914     public double getLMDot() {
915         return getAlphaMDot() + raanDot;
916     }
917 
918     /** Compute position and velocity but not acceleration.
919      */
920     private void computePVWithoutA() {
921 
922         if (partialPV != null) {
923             // already computed
924             return;
925         }
926 
927         // get equinoctial parameters
928         final double equEx = getEquinoctialEx();
929         final double equEy = getEquinoctialEy();
930         final double hx = getHx();
931         final double hy = getHy();
932         final double lE = getLE();
933 
934         // inclination-related intermediate parameters
935         final double hx2   = hx * hx;
936         final double hy2   = hy * hy;
937         final double factH = 1. / (1 + hx2 + hy2);
938 
939         // reference axes defining the orbital plane
940         final double ux = (1 + hx2 - hy2) * factH;
941         final double uy =  2 * hx * hy * factH;
942         final double uz = -2 * hy * factH;
943 
944         final double vx = uy;
945         final double vy = (1 - hx2 + hy2) * factH;
946         final double vz =  2 * hx * factH;
947 
948         // eccentricity-related intermediate parameters
949         final double exey = equEx * equEy;
950         final double ex2  = equEx * equEx;
951         final double ey2  = equEy * equEy;
952         final double e2   = ex2 + ey2;
953         final double eta  = 1 + FastMath.sqrt(1 - e2);
954         final double beta = 1. / eta;
955 
956         // eccentric latitude argument
957         final SinCos scLe   = FastMath.sinCos(lE);
958         final double cLe    = scLe.cos();
959         final double sLe    = scLe.sin();
960         final double exCeyS = equEx * cLe + equEy * sLe;
961 
962         // coordinates of position and velocity in the orbital plane
963         final double x      = a * ((1 - beta * ey2) * cLe + beta * exey * sLe - equEx);
964         final double y      = a * ((1 - beta * ex2) * sLe + beta * exey * cLe - equEy);
965 
966         final double factor = FastMath.sqrt(getMu() / a) / (1 - exCeyS);
967         final double xdot   = factor * (-sLe + beta * equEy * exCeyS);
968         final double ydot   = factor * ( cLe - beta * equEx * exCeyS);
969 
970         final Vector3D position =
971                         new Vector3D(x * ux + y * vx, x * uy + y * vy, x * uz + y * vz);
972         final Vector3D velocity =
973                         new Vector3D(xdot * ux + ydot * vx, xdot * uy + ydot * vy, xdot * uz + ydot * vz);
974 
975         partialPV = new PVCoordinates(position, velocity);
976 
977     }
978 
979     /** Initialize cached alpha with rate.
980      * @param alpha input alpha
981      * @param alphaDot rate of input alpha
982      * @param inputType position angle type passed as input
983      * @return alpha to cache with rate
984      * @since 12.1
985      */
986     private UnivariateDerivative1 initializeCachedAlpha(final double alpha, final double alphaDot,
987                                                         final PositionAngleType inputType) {
988         if (cachedPositionAngleType == inputType) {
989             return new UnivariateDerivative1(alpha, alphaDot);
990 
991         } else {
992             final UnivariateDerivative1 exUD = new UnivariateDerivative1(ex, exDot);
993             final UnivariateDerivative1 eyUD = new UnivariateDerivative1(ey, eyDot);
994             final UnivariateDerivative1 alphaUD = new UnivariateDerivative1(alpha, alphaDot);
995 
996             switch (cachedPositionAngleType) {
997 
998                 case ECCENTRIC:
999                     if (inputType == PositionAngleType.MEAN) {
1000                         return FieldCircularLatitudeArgumentUtility.meanToEccentric(exUD, eyUD, alphaUD);
1001                     } else {
1002                         return FieldCircularLatitudeArgumentUtility.trueToEccentric(exUD, eyUD, alphaUD);
1003                     }
1004 
1005                 case TRUE:
1006                     if (inputType == PositionAngleType.MEAN) {
1007                         return FieldCircularLatitudeArgumentUtility.meanToTrue(exUD, eyUD, alphaUD);
1008                     } else {
1009                         return FieldCircularLatitudeArgumentUtility.eccentricToTrue(exUD, eyUD, alphaUD);
1010                     }
1011 
1012                 case MEAN:
1013                     if (inputType == PositionAngleType.TRUE) {
1014                         return FieldCircularLatitudeArgumentUtility.trueToMean(exUD, eyUD, alphaUD);
1015                     } else {
1016                         return FieldCircularLatitudeArgumentUtility.eccentricToMean(exUD, eyUD, alphaUD);
1017                     }
1018 
1019                 default:
1020                     throw new OrekitInternalError(null);
1021 
1022             }
1023 
1024         }
1025 
1026     }
1027 
1028     /** Initialize cached alpha.
1029      * @param alpha input alpha
1030      * @param positionAngleType position angle type passed as input
1031      * @return alpha to cache
1032      * @since 12.1
1033      */
1034     private double initializeCachedAlpha(final double alpha, final PositionAngleType positionAngleType) {
1035         if (positionAngleType == cachedPositionAngleType) {
1036             return alpha;
1037 
1038         } else {
1039             switch (cachedPositionAngleType) {
1040 
1041                 case ECCENTRIC:
1042                     if (positionAngleType == PositionAngleType.MEAN) {
1043                         return CircularLatitudeArgumentUtility.meanToEccentric(ex, ey, alpha);
1044                     } else {
1045                         return CircularLatitudeArgumentUtility.trueToEccentric(ex, ey, alpha);
1046                     }
1047 
1048                 case MEAN:
1049                     if (positionAngleType == PositionAngleType.TRUE) {
1050                         return CircularLatitudeArgumentUtility.trueToMean(ex, ey, alpha);
1051                     } else {
1052                         return CircularLatitudeArgumentUtility.eccentricToMean(ex, ey, alpha);
1053                     }
1054 
1055                 case TRUE:
1056                     if (positionAngleType == PositionAngleType.MEAN) {
1057                         return CircularLatitudeArgumentUtility.meanToTrue(ex, ey, alpha);
1058                     } else {
1059                         return CircularLatitudeArgumentUtility.eccentricToTrue(ex, ey, alpha);
1060                     }
1061 
1062                 default:
1063                     throw new OrekitInternalError(null);
1064             }
1065         }
1066     }
1067 
1068     /** Compute non-Keplerian part of the acceleration from first time derivatives.
1069      * <p>
1070      * This method should be called only when {@link #hasDerivatives()} returns true.
1071      * </p>
1072      * @return non-Keplerian part of the acceleration
1073      */
1074     private Vector3D nonKeplerianAcceleration() {
1075 
1076         final double[][] dCdP = new double[6][6];
1077         getJacobianWrtParameters(PositionAngleType.MEAN, dCdP);
1078 
1079         final double nonKeplerianMeanMotion = getAlphaMDot() - getKeplerianMeanMotion();
1080         final double nonKeplerianAx = dCdP[3][0] * aDot    + dCdP[3][1] * exDot   + dCdP[3][2] * eyDot   +
1081                                       dCdP[3][3] * iDot    + dCdP[3][4] * raanDot + dCdP[3][5] * nonKeplerianMeanMotion;
1082         final double nonKeplerianAy = dCdP[4][0] * aDot    + dCdP[4][1] * exDot   + dCdP[4][2] * eyDot   +
1083                                       dCdP[4][3] * iDot    + dCdP[4][4] * raanDot + dCdP[4][5] * nonKeplerianMeanMotion;
1084         final double nonKeplerianAz = dCdP[5][0] * aDot    + dCdP[5][1] * exDot   + dCdP[5][2] * eyDot   +
1085                                       dCdP[5][3] * iDot    + dCdP[5][4] * raanDot + dCdP[5][5] * nonKeplerianMeanMotion;
1086 
1087         return new Vector3D(nonKeplerianAx, nonKeplerianAy, nonKeplerianAz);
1088 
1089     }
1090 
1091     /** {@inheritDoc} */
1092     @Override
1093     protected Vector3D initPosition() {
1094 
1095         // get equinoctial parameters
1096         final double equEx = getEquinoctialEx();
1097         final double equEy = getEquinoctialEy();
1098         final double hx = getHx();
1099         final double hy = getHy();
1100         final double lE = getLE();
1101 
1102         // inclination-related intermediate parameters
1103         final double hx2   = hx * hx;
1104         final double hy2   = hy * hy;
1105         final double factH = 1. / (1 + hx2 + hy2);
1106 
1107         // reference axes defining the orbital plane
1108         final double ux = (1 + hx2 - hy2) * factH;
1109         final double uy =  2 * hx * hy * factH;
1110         final double uz = -2 * hy * factH;
1111 
1112         final double vx = uy;
1113         final double vy = (1 - hx2 + hy2) * factH;
1114         final double vz =  2 * hx * factH;
1115 
1116         // eccentricity-related intermediate parameters
1117         final double exey = equEx * equEy;
1118         final double ex2  = equEx * equEx;
1119         final double ey2  = equEy * equEy;
1120         final double e2   = ex2 + ey2;
1121         final double eta  = 1 + FastMath.sqrt(1 - e2);
1122         final double beta = 1. / eta;
1123 
1124         // eccentric latitude argument
1125         final SinCos scLe   = FastMath.sinCos(lE);
1126         final double cLe    = scLe.cos();
1127         final double sLe    = scLe.sin();
1128 
1129         // coordinates of position and velocity in the orbital plane
1130         final double x      = a * ((1 - beta * ey2) * cLe + beta * exey * sLe - equEx);
1131         final double y      = a * ((1 - beta * ex2) * sLe + beta * exey * cLe - equEy);
1132 
1133         return new Vector3D(x * ux + y * vx, x * uy + y * vy, x * uz + y * vz);
1134 
1135     }
1136 
1137     /** {@inheritDoc} */
1138     @Override
1139     protected TimeStampedPVCoordinates initPVCoordinates() {
1140 
1141         // position and velocity
1142         computePVWithoutA();
1143 
1144         // acceleration
1145         final double r2 = partialPV.getPosition().getNormSq();
1146         final Vector3D keplerianAcceleration = new Vector3D(-getMu() / (r2 * FastMath.sqrt(r2)), partialPV.getPosition());
1147         final Vector3D acceleration = hasDerivatives() ?
1148                                       keplerianAcceleration.add(nonKeplerianAcceleration()) :
1149                                       keplerianAcceleration;
1150 
1151         return new TimeStampedPVCoordinates(getDate(), partialPV.getPosition(), partialPV.getVelocity(), acceleration);
1152 
1153     }
1154 
1155     /** {@inheritDoc} */
1156     @Override
1157     public CircularOrbit shiftedBy(final double dt) {
1158 
1159         // use Keplerian-only motion
1160         final CircularOrbit keplerianShifted = new CircularOrbit(a, ex, ey, i, raan,
1161                                                                  getAlphaM() + getKeplerianMeanMotion() * dt,
1162                                                                  PositionAngleType.MEAN, cachedPositionAngleType,
1163                                                                  getFrame(), getDate().shiftedBy(dt), getMu());
1164 
1165         if (hasDerivatives()) {
1166 
1167             // extract non-Keplerian acceleration from first time derivatives
1168             final Vector3D nonKeplerianAcceleration = nonKeplerianAcceleration();
1169 
1170             // add quadratic effect of non-Keplerian acceleration to Keplerian-only shift
1171             keplerianShifted.computePVWithoutA();
1172             final Vector3D fixedP   = new Vector3D(1, keplerianShifted.partialPV.getPosition(),
1173                                                    0.5 * dt * dt, nonKeplerianAcceleration);
1174             final double   fixedR2 = fixedP.getNormSq();
1175             final double   fixedR  = FastMath.sqrt(fixedR2);
1176             final Vector3D fixedV  = new Vector3D(1, keplerianShifted.partialPV.getVelocity(),
1177                                                   dt, nonKeplerianAcceleration);
1178             final Vector3D fixedA  = new Vector3D(-getMu() / (fixedR2 * fixedR), keplerianShifted.partialPV.getPosition(),
1179                                                   1, nonKeplerianAcceleration);
1180 
1181             // build a new orbit, taking non-Keplerian acceleration into account
1182             return new CircularOrbit(new TimeStampedPVCoordinates(keplerianShifted.getDate(),
1183                                                                   fixedP, fixedV, fixedA),
1184                                      keplerianShifted.getFrame(), keplerianShifted.getMu());
1185 
1186         } else {
1187             // Keplerian-only motion is all we can do
1188             return keplerianShifted;
1189         }
1190 
1191     }
1192 
1193     /** {@inheritDoc} */
1194     @Override
1195     protected double[][] computeJacobianMeanWrtCartesian() {
1196 
1197 
1198         final double[][] jacobian = new double[6][6];
1199 
1200         computePVWithoutA();
1201         final Vector3D position = partialPV.getPosition();
1202         final Vector3D velocity = partialPV.getVelocity();
1203         final double x          = position.getX();
1204         final double y          = position.getY();
1205         final double z          = position.getZ();
1206         final double vx         = velocity.getX();
1207         final double vy         = velocity.getY();
1208         final double vz         = velocity.getZ();
1209         final double pv         = Vector3D.dotProduct(position, velocity);
1210         final double r2         = position.getNormSq();
1211         final double r          = FastMath.sqrt(r2);
1212         final double v2         = velocity.getNormSq();
1213 
1214         final double mu         = getMu();
1215         final double oOsqrtMuA  = 1 / FastMath.sqrt(mu * a);
1216         final double rOa        = r / a;
1217         final double aOr        = a / r;
1218         final double aOr2       = a / r2;
1219         final double a2         = a * a;
1220 
1221         final double ex2        = ex * ex;
1222         final double ey2        = ey * ey;
1223         final double e2         = ex2 + ey2;
1224         final double epsilon    = FastMath.sqrt(1 - e2);
1225         final double beta       = 1 / (1 + epsilon);
1226 
1227         final double eCosE      = 1 - rOa;
1228         final double eSinE      = pv * oOsqrtMuA;
1229 
1230         final SinCos scI    = FastMath.sinCos(i);
1231         final SinCos scRaan = FastMath.sinCos(raan);
1232         final double cosI       = scI.cos();
1233         final double sinI       = scI.sin();
1234         final double cosRaan    = scRaan.cos();
1235         final double sinRaan    = scRaan.sin();
1236 
1237         // da
1238         fillHalfRow(2 * aOr * aOr2, position, jacobian[0], 0);
1239         fillHalfRow(2 * a2 / mu, velocity, jacobian[0], 3);
1240 
1241         // differentials of the normalized momentum
1242         final Vector3D danP = new Vector3D(v2, position, -pv, velocity);
1243         final Vector3D danV = new Vector3D(r2, velocity, -pv, position);
1244         final double recip  = 1 / partialPV.getMomentum().getNorm();
1245         final double recip2 = recip * recip;
1246         final Vector3D dwXP = new Vector3D(recip, new Vector3D(  0,  vz, -vy), -recip2 * sinRaan * sinI, danP);
1247         final Vector3D dwYP = new Vector3D(recip, new Vector3D(-vz,   0,  vx),  recip2 * cosRaan * sinI, danP);
1248         final Vector3D dwZP = new Vector3D(recip, new Vector3D( vy, -vx,   0), -recip2 * cosI,           danP);
1249         final Vector3D dwXV = new Vector3D(recip, new Vector3D(  0,  -z,   y), -recip2 * sinRaan * sinI, danV);
1250         final Vector3D dwYV = new Vector3D(recip, new Vector3D(  z,   0,  -x),  recip2 * cosRaan * sinI, danV);
1251         final Vector3D dwZV = new Vector3D(recip, new Vector3D( -y,   x,   0), -recip2 * cosI,           danV);
1252 
1253         // di
1254         fillHalfRow(sinRaan * cosI, dwXP, -cosRaan * cosI, dwYP, -sinI, dwZP, jacobian[3], 0);
1255         fillHalfRow(sinRaan * cosI, dwXV, -cosRaan * cosI, dwYV, -sinI, dwZV, jacobian[3], 3);
1256 
1257         // dRaan
1258         fillHalfRow(sinRaan / sinI, dwYP, cosRaan / sinI, dwXP, jacobian[4], 0);
1259         fillHalfRow(sinRaan / sinI, dwYV, cosRaan / sinI, dwXV, jacobian[4], 3);
1260 
1261         // orbital frame: (p, q, w) p along ascending node, w along momentum
1262         // the coordinates of the spacecraft in this frame are: (u, v, 0)
1263         final double u     =  x * cosRaan + y * sinRaan;
1264         final double cv    = -x * sinRaan + y * cosRaan;
1265         final double v     = cv * cosI + z * sinI;
1266 
1267         // du
1268         final Vector3D duP = new Vector3D(cv * cosRaan / sinI, dwXP,
1269                                           cv * sinRaan / sinI, dwYP,
1270                                           1, new Vector3D(cosRaan, sinRaan, 0));
1271         final Vector3D duV = new Vector3D(cv * cosRaan / sinI, dwXV,
1272                                           cv * sinRaan / sinI, dwYV);
1273 
1274         // dv
1275         final Vector3D dvP = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXP,
1276                                           -u * sinRaan * cosI / sinI - cosRaan * z, dwYP,
1277                                           cv, dwZP,
1278                                           1, new Vector3D(-sinRaan * cosI, cosRaan * cosI, sinI));
1279         final Vector3D dvV = new Vector3D(-u * cosRaan * cosI / sinI + sinRaan * z, dwXV,
1280                                           -u * sinRaan * cosI / sinI - cosRaan * z, dwYV,
1281                                           cv, dwZV);
1282 
1283         final Vector3D dc1P = new Vector3D(aOr2 * (2 * eSinE * eSinE + 1 - eCosE) / r2, position,
1284                                             -2 * aOr2 * eSinE * oOsqrtMuA, velocity);
1285         final Vector3D dc1V = new Vector3D(-2 * aOr2 * eSinE * oOsqrtMuA, position,
1286                                             2 / mu, velocity);
1287         final Vector3D dc2P = new Vector3D(aOr2 * eSinE * (eSinE * eSinE - (1 - e2)) / (r2 * epsilon), position,
1288                                             aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, velocity);
1289         final Vector3D dc2V = new Vector3D(aOr2 * (1 - e2 - eSinE * eSinE) * oOsqrtMuA / epsilon, position,
1290                                             eSinE / (mu * epsilon), velocity);
1291 
1292         final double cof1   = aOr2 * (eCosE - e2);
1293         final double cof2   = aOr2 * epsilon * eSinE;
1294         final Vector3D dexP = new Vector3D(u, dc1P,  v, dc2P, cof1, duP,  cof2, dvP);
1295         final Vector3D dexV = new Vector3D(u, dc1V,  v, dc2V, cof1, duV,  cof2, dvV);
1296         final Vector3D deyP = new Vector3D(v, dc1P, -u, dc2P, cof1, dvP, -cof2, duP);
1297         final Vector3D deyV = new Vector3D(v, dc1V, -u, dc2V, cof1, dvV, -cof2, duV);
1298         fillHalfRow(1, dexP, jacobian[1], 0);
1299         fillHalfRow(1, dexV, jacobian[1], 3);
1300         fillHalfRow(1, deyP, jacobian[2], 0);
1301         fillHalfRow(1, deyV, jacobian[2], 3);
1302 
1303         final double cle = u / a + ex - eSinE * beta * ey;
1304         final double sle = v / a + ey + eSinE * beta * ex;
1305         final double m1  = beta * eCosE;
1306         final double m2  = 1 - m1 * eCosE;
1307         final double m3  = (u * ey - v * ex) + eSinE * beta * (u * ex + v * ey);
1308         final double m4  = -sle + cle * eSinE * beta;
1309         final double m5  = cle + sle * eSinE * beta;
1310         fillHalfRow((2 * m3 / r + aOr * eSinE + m1 * eSinE * (1 + m1 - (1 + aOr) * m2) / epsilon) / r2, position,
1311                     (m1 * m2 / epsilon - 1) * oOsqrtMuA, velocity,
1312                     m4, dexP, m5, deyP, -sle / a, duP, cle / a, dvP,
1313                     jacobian[5], 0);
1314         fillHalfRow((m1 * m2 / epsilon - 1) * oOsqrtMuA, position,
1315                     (2 * m3 + eSinE * a + m1 * eSinE * r * (eCosE * beta * 2 - aOr * m2) / epsilon) / mu, velocity,
1316                     m4, dexV, m5, deyV, -sle / a, duV, cle / a, dvV,
1317                     jacobian[5], 3);
1318 
1319         return jacobian;
1320 
1321     }
1322 
1323     /** {@inheritDoc} */
1324     @Override
1325     protected double[][] computeJacobianEccentricWrtCartesian() {
1326 
1327         // start by computing the Jacobian with mean angle
1328         final double[][] jacobian = computeJacobianMeanWrtCartesian();
1329 
1330         // Differentiating the Kepler equation aM = aE - ex sin aE + ey cos aE leads to:
1331         // daM = (1 - ex cos aE - ey sin aE) daE - sin aE dex + cos aE dey
1332         // which is inverted and rewritten as:
1333         // daE = a/r daM + sin aE a/r dex - cos aE a/r dey
1334         final double alphaE = getAlphaE();
1335         final SinCos scAe   = FastMath.sinCos(alphaE);
1336         final double cosAe  = scAe.cos();
1337         final double sinAe  = scAe.sin();
1338         final double aOr    = 1 / (1 - ex * cosAe - ey * sinAe);
1339 
1340         // update longitude row
1341         final double[] rowEx = jacobian[1];
1342         final double[] rowEy = jacobian[2];
1343         final double[] rowL  = jacobian[5];
1344         for (int j = 0; j < 6; ++j) {
1345             rowL[j] = aOr * (rowL[j] + sinAe * rowEx[j] - cosAe * rowEy[j]);
1346         }
1347 
1348         return jacobian;
1349 
1350     }
1351 
1352     /** {@inheritDoc} */
1353     @Override
1354     protected double[][] computeJacobianTrueWrtCartesian() {
1355 
1356         // start by computing the Jacobian with eccentric angle
1357         final double[][] jacobian = computeJacobianEccentricWrtCartesian();
1358 
1359         // Differentiating the eccentric latitude equation
1360         // tan((aV - aE)/2) = [ex sin aE - ey cos aE] / [sqrt(1-ex^2-ey^2) + 1 - ex cos aE - ey sin aE]
1361         // leads to
1362         // cT (daV - daE) = cE daE + cX dex + cY dey
1363         // with
1364         // cT = [d^2 + (ex sin aE - ey cos aE)^2] / 2
1365         // d  = 1 + sqrt(1-ex^2-ey^2) - ex cos aE - ey sin aE
1366         // cE = (ex cos aE + ey sin aE) (sqrt(1-ex^2-ey^2) + 1) - ex^2 - ey^2
1367         // cX =  sin aE (sqrt(1-ex^2-ey^2) + 1) - ey + ex (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
1368         // cY = -cos aE (sqrt(1-ex^2-ey^2) + 1) + ex + ey (ex sin aE - ey cos aE) / sqrt(1-ex^2-ey^2)
1369         // which can be solved to find the differential of the true latitude
1370         // daV = (cT + cE) / cT daE + cX / cT deX + cY / cT deX
1371         final double alphaE    = getAlphaE();
1372         final SinCos scAe      = FastMath.sinCos(alphaE);
1373         final double cosAe     = scAe.cos();
1374         final double sinAe     = scAe.sin();
1375         final double eSinE     = ex * sinAe - ey * cosAe;
1376         final double ecosE     = ex * cosAe + ey * sinAe;
1377         final double e2        = ex * ex + ey * ey;
1378         final double epsilon   = FastMath.sqrt(1 - e2);
1379         final double onePeps   = 1 + epsilon;
1380         final double d         = onePeps - ecosE;
1381         final double cT        = (d * d + eSinE * eSinE) / 2;
1382         final double cE        = ecosE * onePeps - e2;
1383         final double cX        = ex * eSinE / epsilon - ey + sinAe * onePeps;
1384         final double cY        = ey * eSinE / epsilon + ex - cosAe * onePeps;
1385         final double factorLe  = (cT + cE) / cT;
1386         final double factorEx  = cX / cT;
1387         final double factorEy  = cY / cT;
1388 
1389         // update latitude row
1390         final double[] rowEx = jacobian[1];
1391         final double[] rowEy = jacobian[2];
1392         final double[] rowA = jacobian[5];
1393         for (int j = 0; j < 6; ++j) {
1394             rowA[j] = factorLe * rowA[j] + factorEx * rowEx[j] + factorEy * rowEy[j];
1395         }
1396 
1397         return jacobian;
1398 
1399     }
1400 
1401     /** {@inheritDoc} */
1402     @Override
1403     public void addKeplerContribution(final PositionAngleType type, final double gm,
1404                                       final double[] pDot) {
1405         final double oMe2;
1406         final double ksi;
1407         final double n  = FastMath.sqrt(gm / a) / a;
1408         final SinCos sc;
1409         switch (type) {
1410             case MEAN :
1411                 pDot[5] += n;
1412                 break;
1413             case ECCENTRIC :
1414                 sc = FastMath.sinCos(getAlphaE());
1415                 ksi   = 1. / (1 - ex * sc.cos() - ey * sc.sin());
1416                 pDot[5] += n * ksi;
1417                 break;
1418             case TRUE :
1419                 sc = FastMath.sinCos(getAlphaV());
1420                 oMe2  = 1 - ex * ex - ey * ey;
1421                 ksi   = 1 + ex * sc.cos() + ey * sc.sin();
1422                 pDot[5] += n * ksi * ksi / (oMe2 * FastMath.sqrt(oMe2));
1423                 break;
1424             default :
1425                 throw new OrekitInternalError(null);
1426         }
1427     }
1428 
1429     /**  Returns a string representation of this Orbit object.
1430      * @return a string representation of this object
1431      */
1432     public String toString() {
1433         return new StringBuilder().append("circular parameters: ").append('{').
1434                                   append("a: ").append(a).
1435                                   append(", ex: ").append(ex).append(", ey: ").append(ey).
1436                                   append(", i: ").append(FastMath.toDegrees(i)).
1437                                   append(", raan: ").append(FastMath.toDegrees(raan)).
1438                                   append(", alphaV: ").append(FastMath.toDegrees(getAlphaV())).
1439                                   append(";}").toString();
1440     }
1441 
1442     /** {@inheritDoc} */
1443     @Override
1444     public PositionAngleType getCachedPositionAngleType() {
1445         return cachedPositionAngleType;
1446     }
1447 
1448     /** {@inheritDoc} */
1449     @Override
1450     public boolean hasRates() {
1451         return hasDerivatives();
1452     }
1453 
1454     /** {@inheritDoc} */
1455     @Override
1456     public CircularOrbit removeRates() {
1457         final PositionAngleType positionAngleType = getCachedPositionAngleType();
1458         return new CircularOrbit(a, ex, ey, i, raan, cachedAlpha, positionAngleType, positionAngleType,
1459                 getFrame(), getDate(), getMu());
1460     }
1461 
1462     /** Replace the instance with a data transfer object for serialization.
1463      * @return data transfer object that will be serialized
1464      */
1465     @DefaultDataContext
1466     private Object writeReplace() {
1467         return new DTO(this);
1468     }
1469 
1470     /** Internal class used only for serialization. */
1471     @DefaultDataContext
1472     private static class DTO implements Serializable {
1473 
1474         /** Serializable UID. */
1475         private static final long serialVersionUID = 20231217L;
1476 
1477         /** Double values. */
1478         private final double[] d;
1479 
1480         /** Frame in which are defined the orbital parameters. */
1481         private final Frame frame;
1482 
1483         /** Type of cached position angle. */
1484         private final PositionAngleType positionAngleType;
1485 
1486         /** Simple constructor.
1487          * @param orbit instance to serialize
1488          */
1489         private DTO(final CircularOrbit orbit) {
1490 
1491             final AbsoluteDate date = orbit.getDate();
1492             positionAngleType = orbit.getCachedPositionAngleType();
1493 
1494             // decompose date
1495             final AbsoluteDate j2000Epoch =
1496                     DataContext.getDefault().getTimeScales().getJ2000Epoch();
1497             final double epoch  = FastMath.floor(date.durationFrom(j2000Epoch));
1498             final double offset = date.durationFrom(j2000Epoch.shiftedBy(epoch));
1499 
1500             if (orbit.serializePV) {
1501                 final TimeStampedPVCoordinates pv = orbit.getPVCoordinates();
1502                 if (orbit.hasDerivatives()) {
1503                     this.d = new double[] {
1504                         // date + mu + orbit + derivatives + Cartesian : 24 parameters
1505                         epoch, offset, orbit.getMu(),
1506                         orbit.a, orbit.ex, orbit.ey,
1507                         orbit.i, orbit.raan, orbit.cachedAlpha,
1508                         orbit.aDot, orbit.exDot, orbit.eyDot,
1509                         orbit.iDot, orbit.raanDot, orbit.cachedAlphaDot,
1510                         pv.getPosition().getX(),     pv.getPosition().getY(),     pv.getPosition().getZ(),
1511                         pv.getVelocity().getX(),     pv.getVelocity().getY(),     pv.getVelocity().getZ(),
1512                         pv.getAcceleration().getX(), pv.getAcceleration().getY(), pv.getAcceleration().getZ(),
1513                     };
1514                 } else {
1515                     this.d = new double[] {
1516                         // date + mu + orbit + Cartesian : 18 parameters
1517                         epoch, offset, orbit.getMu(),
1518                         orbit.a, orbit.ex, orbit.ey,
1519                         orbit.i, orbit.raan, orbit.cachedAlpha,
1520                         pv.getPosition().getX(),     pv.getPosition().getY(),     pv.getPosition().getZ(),
1521                         pv.getVelocity().getX(),     pv.getVelocity().getY(),     pv.getVelocity().getZ(),
1522                         pv.getAcceleration().getX(), pv.getAcceleration().getY(), pv.getAcceleration().getZ(),
1523                     };
1524                 }
1525             } else {
1526                 if (orbit.hasDerivatives()) {
1527                     // date + mu + orbit + derivatives: 15 parameters
1528                     this.d = new double[] {
1529                         epoch, offset, orbit.getMu(),
1530                         orbit.a, orbit.ex, orbit.ey,
1531                         orbit.i, orbit.raan, orbit.cachedAlpha,
1532                         orbit.aDot, orbit.exDot, orbit.eyDot,
1533                         orbit.iDot, orbit.raanDot, orbit.cachedAlphaDot
1534                     };
1535                 } else {
1536                     // date + mu + orbit: 9 parameters
1537                     this.d = new double[] {
1538                         epoch, offset, orbit.getMu(),
1539                         orbit.a, orbit.ex, orbit.ey,
1540                         orbit.i, orbit.raan, orbit.cachedAlpha
1541                     };
1542                 }
1543             }
1544 
1545             this.frame = orbit.getFrame();
1546 
1547         }
1548 
1549         /** Replace the deserialized data transfer object with a {@link CircularOrbit}.
1550          * @return replacement {@link CircularOrbit}
1551          */
1552         private Object readResolve() {
1553             final AbsoluteDate j2000Epoch =
1554                     DataContext.getDefault().getTimeScales().getJ2000Epoch();
1555             switch (d.length) {
1556                 case 24 : // date + mu + orbit + derivatives + Cartesian
1557                     return new CircularOrbit(d[ 3], d[ 4], d[ 5], d[ 6], d[ 7], d[ 8],
1558                                              d[ 9], d[10], d[11], d[12], d[13], d[14],
1559                                              new TimeStampedPVCoordinates(j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
1560                                                                           new Vector3D(d[15], d[16], d[17]),
1561                                                                           new Vector3D(d[18], d[19], d[20]),
1562                                                                           new Vector3D(d[21], d[22], d[23])),
1563                                              positionAngleType, frame,
1564                                              d[2]);
1565                 case 18 : // date + mu + orbit + Cartesian
1566                     return new CircularOrbit(d[3], d[4], d[5], d[6], d[7], d[8],
1567                                              Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN, Double.NaN,
1568                                              new TimeStampedPVCoordinates(j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
1569                                                                           new Vector3D(d[ 9], d[10], d[11]),
1570                                                                           new Vector3D(d[12], d[13], d[14]),
1571                                                                           new Vector3D(d[15], d[16], d[17])),
1572                                              positionAngleType, frame,
1573                                              d[2]);
1574                 case 15 : // date + mu + orbit + derivatives
1575                     return new CircularOrbit(d[ 3], d[ 4], d[ 5], d[ 6], d[ 7], d[ 8],
1576                                              d[ 9], d[10], d[11], d[12], d[13], d[14],
1577                                              positionAngleType, positionAngleType,
1578                                              frame, j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
1579                                              d[2]);
1580                 default : // date + mu + orbit
1581                     return new CircularOrbit(d[3], d[4], d[5], d[6], d[7], d[8], positionAngleType, positionAngleType,
1582                                              frame, j2000Epoch.shiftedBy(d[0]).shiftedBy(d[1]),
1583                                              d[2]);
1584 
1585             }
1586         }
1587 
1588     }
1589 
1590 }