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