EcksteinHechlerPropagator.java

  1. /* Copyright 2002-2016 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.propagation.analytical;

  19. import java.io.NotSerializableException;
  20. import java.io.Serializable;
  21. import java.util.ArrayList;
  22. import java.util.List;
  23. import java.util.SortedSet;

  25. import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
  26. import org.apache.commons.math3.geometry.euclidean.threed.FieldVector3D;
  27. import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
  28. import org.apache.commons.math3.util.FastMath;
  29. import org.apache.commons.math3.util.MathUtils;
  30. import org.orekit.attitudes.AttitudeProvider;
  31. import org.orekit.errors.OrekitException;
  32. import org.orekit.errors.OrekitInternalError;
  33. import org.orekit.errors.OrekitMessages;
  34. import org.orekit.errors.PropagationException;
  35. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  36. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  37. import org.orekit.orbits.CartesianOrbit;
  38. import org.orekit.orbits.CircularOrbit;
  39. import org.orekit.orbits.Orbit;
  40. import org.orekit.orbits.OrbitType;
  41. import org.orekit.orbits.PositionAngle;
  42. import org.orekit.propagation.AdditionalStateProvider;
  43. import org.orekit.propagation.SpacecraftState;
  44. import org.orekit.time.AbsoluteDate;
  45. import org.orekit.utils.TimeSpanMap;
  46. import org.orekit.utils.TimeStampedPVCoordinates;

  48. /** This class propagates a {@link org.orekit.propagation.SpacecraftState}
  49.  *  using the analytical Eckstein-Hechler model.
  50.  * <p>The Eckstein-Hechler model is suited for near circular orbits
  51.  * (e < 0.1, with poor accuracy between 0.005 and 0.1) and inclination
  52.  * neither equatorial (direct or retrograde) nor critical (direct or
  53.  * retrograde).</p>
  54.  * <p>
  55.  * Note that before version 7.0, there was a large inconsistency in the generated
  56.  * orbits, and it was fixed as of version 7.0 of Orekit, with a visible side effect.
  57.  * The problems is that if the circular parameters produced by the Eckstein-Hechler
  58.  * model are used to build an orbit considered to be osculating, the velocity deduced
  59.  * from this orbit was <em>inconsistent with the position evolution</em>! The reason is
  60.  * that the model includes non-Keplerian effects but it does not include a corresponding
  61.  * circular/Cartesian conversion. As a consequence, all subsequent computation involving
  62.  * velocity were wrong. This includes attitude modes like yaw compensation and Doppler
  63.  * effect. As this effect was considered serious enough and as accurate velocities were
  64.  * considered important, the propagator now generates {@link CartesianOrbit Cartesian
  65.  * orbits} which are built in a special way to ensure consistency throughout propagation.
  66.  * A side effect is that if circular parameters are rebuilt by user from these propagated
  67.  * Cartesian orbit, the circular parameters will generally <em>not</em> match the initial
  68.  * orbit (differences in semi-major axis can exceed 120 m). The position however <em>will</em>
  69.  * match to sub-micrometer level, and this position will be identical to the positions
  70.  * that were generated by previous versions (in other words, the internals of the models
  71.  * have not been changed, only the output parameters have been changed). The correctness
  72.  * of the initialization has been assessed and is good, as it allows the subsequent orbit
  73.  * to remain close to a numerical reference orbit.
  74.  * </p>
  75.  * <p>
  76.  * If users need a more definitive initialization of an Eckstein-Hechler propagator, they
  77.  * should consider using a {@link org.orekit.propagation.conversion.PropagatorConverter
  78.  * propagator converter} to initialize their Eckstein-Hechler propagator using a complete
  79.  * sample instead of just a single initial orbit.
  80.  * </p>
  81.  * @see Orbit
  82.  * @author Guylaine Prat
  83.  */
  84. public class EcksteinHechlerPropagator extends AbstractAnalyticalPropagator implements Serializable {

  86.     /** Serializable UID. */
  87.     private static final long serialVersionUID = 20151202L;

  89.     /** Initial Eckstein-Hechler model. */
  90.     private EHModel initialModel;

  92.     /** All models. */
  93.     private transient TimeSpanMap<EHModel> models;

  95.     /** Reference radius of the central body attraction model (m). */
  96.     private double referenceRadius;

  98.     /** Central attraction coefficient (m³/s²). */
  99.     private double mu;

  101.     /** Un-normalized zonal coefficients. */
  102.     private double[] ck0;

  104.     /** Build a propagator from orbit and potential provider.
  105.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  106.      * @param initialOrbit initial orbit
  107.      * @param provider for un-normalized zonal coefficients
  108.      * @exception OrekitException if the zonal coefficients cannot be retrieved
  109.      * @exception PropagationException if the mean parameters cannot be computed
  110.      */
  111.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  112.                                      final UnnormalizedSphericalHarmonicsProvider provider)
  113.         throws PropagationException, OrekitException {
  114.         this(initialOrbit, DEFAULT_LAW, DEFAULT_MASS, provider,
  115.              provider.onDate(initialOrbit.getDate()));
  116.     }

  118.     /**
  119.      * Private helper constructor.
  120.      * @param initialOrbit initial orbit
  121.      * @param attitude attitude provider
  122.      * @param mass spacecraft mass
  123.      * @param provider for un-normalized zonal coefficients
  124.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  125.      * @exception OrekitException if the zonal coefficients cannot be retrieved
  126.      */
  127.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  128.                                      final AttitudeProvider attitude,
  129.                                      final double mass,
  130.                                      final UnnormalizedSphericalHarmonicsProvider provider,
  131.                                      final UnnormalizedSphericalHarmonics harmonics)
  132.         throws OrekitException {
  133.         this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
  134.              harmonics.getUnnormalizedCnm(2, 0),
  135.              harmonics.getUnnormalizedCnm(3, 0),
  136.              harmonics.getUnnormalizedCnm(4, 0),
  137.              harmonics.getUnnormalizedCnm(5, 0),
  138.              harmonics.getUnnormalizedCnm(6, 0));
  139.     }

  141.     /** Build a propagator from orbit and potential.
  142.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  143.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  144.      * are related to both the normalized coefficients
  145.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  146.      *  and the J<sub>n</sub> one as follows:</p>
  147.      * <pre>
  148.      *   C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  149.      *                      <span style="text-decoration: overline">C</span><sub>n,0</sub>
  150.      *   C<sub>n,0</sub> = -J<sub>n</sub>
  151.      * </pre>
  152.      * @param initialOrbit initial orbit
  153.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  154.      * @param mu central attraction coefficient (m³/s²)
  155.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  156.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  157.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  158.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  159.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  160.      * @exception PropagationException if the mean parameters cannot be computed
  161.      * @see org.orekit.utils.Constants
  162.      */
  163.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  164.                                      final double referenceRadius, final double mu,
  165.                                      final double c20, final double c30, final double c40,
  166.                                      final double c50, final double c60)
  167.         throws PropagationException {
  168.         this(initialOrbit, DEFAULT_LAW, DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
  169.     }

  171.     /** Build a propagator from orbit, mass and potential provider.
  172.      * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
  173.      * @param initialOrbit initial orbit
  174.      * @param mass spacecraft mass
  175.      * @param provider for un-normalized zonal coefficients
  176.      * @exception OrekitException if the zonal coefficients cannot be retrieved
  177.      * @exception PropagationException if the mean parameters cannot be computed
  178.      */
  179.     public EcksteinHechlerPropagator(final Orbit initialOrbit, final double mass,
  180.                                      final UnnormalizedSphericalHarmonicsProvider provider)
  181.         throws PropagationException, OrekitException {
  182.         this(initialOrbit, DEFAULT_LAW, mass, provider, provider.onDate(initialOrbit.getDate()));
  183.     }

  185.     /** Build a propagator from orbit, mass and potential.
  186.      * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
  187.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  188.      * are related to both the normalized coefficients
  189.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  190.      *  and the J<sub>n</sub> one as follows:</p>
  191.      * <pre>
  192.      *   C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  193.      *                      <span style="text-decoration: overline">C</span><sub>n,0</sub>
  194.      *   C<sub>n,0</sub> = -J<sub>n</sub>
  195.      * </pre>
  196.      * @param initialOrbit initial orbit
  197.      * @param mass spacecraft mass
  198.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  199.      * @param mu central attraction coefficient (m³/s²)
  200.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  201.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  202.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  203.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  204.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  205.      * @exception PropagationException if the mean parameters cannot be computed
  206.      */
  207.     public EcksteinHechlerPropagator(final Orbit initialOrbit, final double mass,
  208.                                      final double referenceRadius, final double mu,
  209.                                      final double c20, final double c30, final double c40,
  210.                                      final double c50, final double c60)
  211.         throws PropagationException {
  212.         this(initialOrbit, DEFAULT_LAW, mass, referenceRadius, mu, c20, c30, c40, c50, c60);
  213.     }

  215.     /** Build a propagator from orbit, attitude provider and potential provider.
  216.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  217.      * @param initialOrbit initial orbit
  218.      * @param attitudeProv attitude provider
  219.      * @param provider for un-normalized zonal coefficients
  220.      * @exception OrekitException if the zonal coefficients cannot be retrieved
  221.      * @exception PropagationException if the mean parameters cannot be computed
  222.      */
  223.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  224.                                      final AttitudeProvider attitudeProv,
  225.                                      final UnnormalizedSphericalHarmonicsProvider provider)
  226.         throws PropagationException, OrekitException {
  227.         this(initialOrbit, attitudeProv, DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()));
  228.     }

  230.     /** Build a propagator from orbit, attitude provider and potential.
  231.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  232.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  233.      * are related to both the normalized coefficients
  234.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  235.      *  and the J<sub>n</sub> one as follows:</p>
  236.      * <pre>
  237.      *   C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  238.      *                     <span style="text-decoration: overline">C</span><sub>n,0</sub>
  239.      *   C<sub>n,0</sub> = -J<sub>n</sub>
  240.      * </pre>
  241.      * @param initialOrbit initial orbit
  242.      * @param attitudeProv attitude provider
  243.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  244.      * @param mu central attraction coefficient (m³/s²)
  245.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  246.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  247.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  248.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  249.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  250.      * @exception PropagationException if the mean parameters cannot be computed
  251.      */
  252.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  253.                                      final AttitudeProvider attitudeProv,
  254.                                      final double referenceRadius, final double mu,
  255.                                      final double c20, final double c30, final double c40,
  256.                                      final double c50, final double c60)
  257.         throws PropagationException {
  258.         this(initialOrbit, attitudeProv, DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
  259.     }

  261.     /** Build a propagator from orbit, attitude provider, mass and potential provider.
  262.      * @param initialOrbit initial orbit
  263.      * @param attitudeProv attitude provider
  264.      * @param mass spacecraft mass
  265.      * @param provider for un-normalized zonal coefficients
  266.      * @exception OrekitException if the zonal coefficients cannot be retrieved
  267.      * @exception PropagationException if the mean parameters cannot be computed
  268.      */
  269.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  270.                                      final AttitudeProvider attitudeProv,
  271.                                      final double mass,
  272.                                      final UnnormalizedSphericalHarmonicsProvider provider)
  273.         throws PropagationException, OrekitException {
  274.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()));
  275.     }

  277.     /** Build a propagator from orbit, attitude provider, mass and potential.
  278.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  279.      * are related to both the normalized coefficients
  280.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  281.      *  and the J<sub>n</sub> one as follows:</p>
  282.      * <pre>
  283.      *   C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  284.      *                      <span style="text-decoration: overline">C</span><sub>n,0</sub>
  285.      *   C<sub>n,0</sub> = -J<sub>n</sub>
  286.      * </pre>
  287.      * @param initialOrbit initial orbit
  288.      * @param attitudeProv attitude provider
  289.      * @param mass spacecraft mass
  290.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  291.      * @param mu central attraction coefficient (m³/s²)
  292.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  293.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  294.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  295.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  296.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  297.      * @exception PropagationException if the mean parameters cannot be computed
  298.      */
  299.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  300.                                      final AttitudeProvider attitudeProv,
  301.                                      final double mass,
  302.                                      final double referenceRadius, final double mu,
  303.                                      final double c20, final double c30, final double c40,
  304.                                      final double c50, final double c60)
  305.         throws PropagationException {

  307.         super(attitudeProv);

  309.         try {

  311.             // store model coefficients
  312.             this.referenceRadius = referenceRadius;
  313.             this.mu  = mu;
  314.             this.ck0 = new double[] {
  315.                 0.0, 0.0, c20, c30, c40, c50, c60
  316.             };

  318.             // compute mean parameters
  319.             // transform into circular adapted parameters used by the Eckstein-Hechler model
  320.             resetInitialState(new SpacecraftState(initialOrbit,
  321.                                                   attitudeProv.getAttitude(initialOrbit,
  322.                                                                            initialOrbit.getDate(),
  323.                                                                            initialOrbit.getFrame()),
  324.                                                   mass));

  326.         } catch (OrekitException oe) {
  327.             throw new PropagationException(oe);
  328.         }
  329.     }

  331.     /** {@inheritDoc} */
  332.     public void resetInitialState(final SpacecraftState state)
  333.         throws PropagationException {
  334.         super.resetInitialState(state);
  335.         this.initialModel = computeMeanParameters((CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit()),
  336.                                                   state.getMass());
  337.         this.models       = new TimeSpanMap<EHModel>(initialModel);
  338.     }

  340.     /** {@inheritDoc} */
  341.     protected void resetIntermediateState(final SpacecraftState state, final boolean forward)
  342.         throws PropagationException {
  343.         super.resetIntermediateState(state, forward);
  344.         final EHModel newModel = computeMeanParameters((CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit()),
  345.                                                        state.getMass());
  346.         if (forward) {
  347.             models.addValidAfter(newModel, state.getDate());
  348.         } else {
  349.             models.addValidBefore(newModel, state.getDate());
  350.         }
  351.     }

  353.     /** Compute mean parameters according to the Eckstein-Hechler analytical model.
  354.      * @param osculating osculating orbit
  355.      * @param mass constant mass
  356.      * @return Eckstein-Hechler mean model
  357.      * @exception PropagationException if orbit goes outside of supported range
  358.      * (trajectory inside the Brillouin sphere, too eccentric, equatorial, critical
  359.      * inclination) or if convergence cannot be reached
  360.      */
  361.     private EHModel computeMeanParameters(final CircularOrbit osculating, final double mass)
  362.         throws PropagationException {

  364.         // sanity check
  365.         if (osculating.getA() < referenceRadius) {
  366.             throw new PropagationException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
  367.                                            osculating.getA());
  368.         }

  370.         // rough initialization of the mean parameters
  371.         EHModel current = new EHModel(osculating, mass, referenceRadius, mu, ck0);

  373.         // threshold for each parameter
  374.         final double epsilon         = 1.0e-13;
  375.         final double thresholdA      = epsilon * (1 + FastMath.abs(current.mean.getA()));
  376.         final double thresholdE      = epsilon * (1 + current.mean.getE());
  377.         final double thresholdAngles = epsilon * FastMath.PI;

  379.         int i = 0;
  380.         while (i++ < 100) {

  382.             // recompute the osculating parameters from the current mean parameters
  383.             final DerivativeStructure[] parameters = current.propagateParameters(current.mean.getDate());

  385.             // adapted parameters residuals
  386.             final double deltaA      = osculating.getA()          - parameters[0].getValue();
  387.             final double deltaEx     = osculating.getCircularEx() - parameters[1].getValue();
  388.             final double deltaEy     = osculating.getCircularEy() - parameters[2].getValue();
  389.             final double deltaI      = osculating.getI()          - parameters[3].getValue();
  390.             final double deltaRAAN   = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode() -
  391.                                                                 parameters[4].getValue(),
  392.                                                                 0.0);
  393.             final double deltaAlphaM = MathUtils.normalizeAngle(osculating.getAlphaM() - parameters[5].getValue(), 0.0);

  395.             // update mean parameters
  396.             current = new EHModel(new CircularOrbit(current.mean.getA()          + deltaA,
  397.                                                     current.mean.getCircularEx() + deltaEx,
  398.                                                     current.mean.getCircularEy() + deltaEy,
  399.                                                     current.mean.getI()          + deltaI,
  400.                                                     current.mean.getRightAscensionOfAscendingNode() + deltaRAAN,
  401.                                                     current.mean.getAlphaM()     + deltaAlphaM,
  402.                                                     PositionAngle.MEAN,
  403.                                                     current.mean.getFrame(),
  404.                                                     current.mean.getDate(), mu),
  405.                                   mass, referenceRadius, mu, ck0);

  407.             // check convergence
  408.             if ((FastMath.abs(deltaA)      < thresholdA) &&
  409.                 (FastMath.abs(deltaEx)     < thresholdE) &&
  410.                 (FastMath.abs(deltaEy)     < thresholdE) &&
  411.                 (FastMath.abs(deltaI)      < thresholdAngles) &&
  412.                 (FastMath.abs(deltaRAAN)   < thresholdAngles) &&
  413.                 (FastMath.abs(deltaAlphaM) < thresholdAngles)) {
  414.                 return current;
  415.             }

  417.         }

  419.         throw new PropagationException(OrekitMessages.UNABLE_TO_COMPUTE_ECKSTEIN_HECHLER_MEAN_PARAMETERS, i);

  421.     }

  423.     /** {@inheritDoc} */
  424.     public CartesianOrbit propagateOrbit(final AbsoluteDate date)
  425.         throws PropagationException {
  426.         // compute Cartesian parameters, taking derivatives into account
  427.         // to make sure velocity and acceleration are consistent
  428.         final EHModel current = models.get(date);
  429.         return new CartesianOrbit(toCartesian(date, current.propagateParameters(date)),
  430.                                   current.mean.getFrame(), mu);
  431.     }

  433.     /** Local class for Eckstein-Hechler model, with fixed mean parameters. */
  434.     private static class EHModel implements Serializable {

  436.         /** Serializable UID. */
  437.         private static final long serialVersionUID = 20160115L;

  439.         /** Mean orbit. */
  440.         private final CircularOrbit mean;

  442.         /** Constant mass. */
  443.         private final double mass;

  445.         // CHECKSTYLE: stop JavadocVariable check

  447.         // preprocessed values
  448.         private final double xnotDot;
  449.         private final double rdpom;
  450.         private final double rdpomp;
  451.         private final double eps1;
  452.         private final double eps2;
  453.         private final double xim;
  454.         private final double ommD;
  455.         private final double rdl;
  456.         private final double aMD;

  458.         private final double kh;
  459.         private final double kl;

  461.         private final double ax1;
  462.         private final double ay1;
  463.         private final double as1;
  464.         private final double ac2;
  465.         private final double axy3;
  466.         private final double as3;
  467.         private final double ac4;
  468.         private final double as5;
  469.         private final double ac6;

  471.         private final double ex1;
  472.         private final double exx2;
  473.         private final double exy2;
  474.         private final double ex3;
  475.         private final double ex4;

  477.         private final double ey1;
  478.         private final double eyx2;
  479.         private final double eyy2;
  480.         private final double ey3;
  481.         private final double ey4;

  483.         private final double rx1;
  484.         private final double ry1;
  485.         private final double r2;
  486.         private final double r3;
  487.         private final double rl;

  489.         private final double iy1;
  490.         private final double ix1;
  491.         private final double i2;
  492.         private final double i3;
  493.         private final double ih;

  495.         private final double lx1;
  496.         private final double ly1;
  497.         private final double l2;
  498.         private final double l3;
  499.         private final double ll;

  501.         // CHECKSTYLE: resume JavadocVariable check

  503.         /** Create a model for specified mean orbit.
  504.          * @param mean mean orbit
  505.          * @param mass constant mass
  506.          * @param referenceRadius reference radius of the central body attraction model (m)
  507.          * @param mu central attraction coefficient (m³/s²)
  508.          * @param ck0 un-normalized zonal coefficients
  509.          * @exception PropagationException if mean orbit is not within model supported domain
  510.          */
  511.         EHModel(final CircularOrbit mean, final double mass,
  512.                 final double referenceRadius, final double mu, final double[] ck0)
  513.             throws PropagationException {

  515.             this.mean            = mean;
  516.             this.mass            = mass;

  518.             // preliminary processing
  519.             double q = referenceRadius / mean.getA();
  520.             double ql = q * q;
  521.             final double g2 = ck0[2] * ql;
  522.             ql *= q;
  523.             final double g3 = ck0[3] * ql;
  524.             ql *= q;
  525.             final double g4 = ck0[4] * ql;
  526.             ql *= q;
  527.             final double g5 = ck0[5] * ql;
  528.             ql *= q;
  529.             final double g6 = ck0[6] * ql;

  531.             final double cosI1 = FastMath.cos(mean.getI());
  532.             final double sinI1 = FastMath.sin(mean.getI());
  533.             final double sinI2 = sinI1 * sinI1;
  534.             final double sinI4 = sinI2 * sinI2;
  535.             final double sinI6 = sinI2 * sinI4;

  537.             if (sinI2 < 1.0e-10) {
  538.                 throw new PropagationException(OrekitMessages.ALMOST_EQUATORIAL_ORBIT,
  539.                                                FastMath.toDegrees(mean.getI()));
  540.             }

  542.             if (FastMath.abs(sinI2 - 4.0 / 5.0) < 1.0e-3) {
  543.                 throw new PropagationException(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT,
  544.                                                FastMath.toDegrees(mean.getI()));
  545.             }

  547.             if (mean.getE() > 0.1) {
  548.                 // if 0.005 < e < 0.1 no error is triggered, but accuracy is poor
  549.                 throw new PropagationException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
  550.                                                mean.getE());
  551.             }

  553.             xnotDot = FastMath.sqrt(mu / mean.getA()) / mean.getA();

  555.             rdpom = -0.75 * g2 * (4.0 - 5.0 * sinI2);
  556.             rdpomp = 7.5 * g4 * (1.0 - 31.0 / 8.0 * sinI2 + 49.0 / 16.0 * sinI4) -
  557.                     13.125 * g6 * (1.0 - 8.0 * sinI2 + 129.0 / 8.0 * sinI4 - 297.0 / 32.0 * sinI6);

  559.             q = 3.0 / (32.0 * rdpom);
  560.             eps1 = q * g4 * sinI2 * (30.0 - 35.0 * sinI2) -
  561.                     175.0 * q * g6 * sinI2 * (1.0 - 3.0 * sinI2 + 2.0625 * sinI4);
  562.             q = 3.0 * sinI1 / (8.0 * rdpom);
  563.             eps2 = q * g3 * (4.0 - 5.0 * sinI2) - q * g5 * (10.0 - 35.0 * sinI2 + 26.25 * sinI4);

  565.             xim = mean.getI();
  566.             ommD = cosI1 * (1.50    * g2 - 2.25 * g2 * g2 * (2.5 - 19.0 / 6.0 * sinI2) +
  567.                             0.9375  * g4 * (7.0 * sinI2 - 4.0) +
  568.                             3.28125 * g6 * (2.0 - 9.0 * sinI2 + 8.25 * sinI4));

  570.             rdl = 1.0 - 1.50 * g2 * (3.0 - 4.0 * sinI2);
  571.             aMD = rdl +
  572.                     2.25 * g2 * g2 * (9.0 - 263.0 / 12.0 * sinI2 + 341.0 / 24.0 * sinI4) +
  573.                     15.0 / 16.0 * g4 * (8.0 - 31.0 * sinI2 + 24.5 * sinI4) +
  574.                     105.0 / 32.0 * g6 * (-10.0 / 3.0 + 25.0 * sinI2 - 48.75 * sinI4 + 27.5 * sinI6);

  576.             final double qq = -1.5 * g2 / rdl;
  577.             final double qA   = 0.75 * g2 * g2 * sinI2;
  578.             final double qB   = 0.25 * g4 * sinI2;
  579.             final double qC   = 105.0 / 16.0 * g6 * sinI2;
  580.             final double qD   = -0.75 * g3 * sinI1;
  581.             final double qE   = 3.75 * g5 * sinI1;
  582.             kh = 0.375 / rdpom;
  583.             kl = kh / sinI1;

  585.             ax1 = qq * (2.0 - 3.5 * sinI2);
  586.             ay1 = qq * (2.0 - 2.5 * sinI2);
  587.             as1 = qD * (4.0 - 5.0 * sinI2) +
  588.                   qE * (2.625 * sinI4 - 3.5 * sinI2 + 1.0);
  589.             ac2 = qq * sinI2 +
  590.                   qA * 7.0 * (2.0 - 3.0 * sinI2) +
  591.                   qB * (15.0 - 17.5 * sinI2) +
  592.                   qC * (3.0 * sinI2 - 1.0 - 33.0 / 16.0 * sinI4);
  593.             axy3 = qq * 3.5 * sinI2;
  594.             as3 = qD * 5.0 / 3.0 * sinI2 +
  595.                   qE * 7.0 / 6.0 * sinI2 * (1.0 - 1.125 * sinI2);
  596.             ac4 = qA * sinI2 +
  597.                   qB * 4.375 * sinI2 +
  598.                   qC * 0.75 * (1.1 * sinI4 - sinI2);

  600.             as5 = qE * 21.0 / 80.0 * sinI4;

  602.             ac6 = qC * -11.0 / 80.0 * sinI4;

  604.             ex1 = qq * (1.0 - 1.25 * sinI2);
  605.             exx2 = qq * 0.5 * (3.0 - 5.0 * sinI2);
  606.             exy2 = qq * (2.0 - 1.5 * sinI2);
  607.             ex3 = qq * 7.0 / 12.0 * sinI2;
  608.             ex4 = qq * 17.0 / 8.0 * sinI2;

  610.             ey1 = qq * (1.0 - 1.75 * sinI2);
  611.             eyx2 = qq * (1.0 - 3.0 * sinI2);
  612.             eyy2 = qq * (2.0 * sinI2 - 1.5);
  613.             ey3 = qq * 7.0 / 12.0 * sinI2;
  614.             ey4 = qq * 17.0 / 8.0 * sinI2;

  616.             q  = -qq * cosI1;
  617.             rx1 =  3.5 * q;
  618.             ry1 = -2.5 * q;
  619.             r2 = -0.5 * q;
  620.             r3 =  7.0 / 6.0 * q;
  621.             rl = g3 * cosI1 * (4.0 - 15.0 * sinI2) -
  622.                  2.5 * g5 * cosI1 * (4.0 - 42.0 * sinI2 + 52.5 * sinI4);

  624.             q = 0.5 * qq * sinI1 * cosI1;
  625.             iy1 =  q;
  626.             ix1 = -q;
  627.             i2 =  q;
  628.             i3 =  q * 7.0 / 3.0;
  629.             ih = -g3 * cosI1 * (4.0 - 5.0 * sinI2) +
  630.                  2.5 * g5 * cosI1 * (4.0 - 14.0 * sinI2 + 10.5 * sinI4);

  632.             lx1 = qq * (7.0 - 77.0 / 8.0 * sinI2);
  633.             ly1 = qq * (55.0 / 8.0 * sinI2 - 7.50);
  634.             l2 = qq * (1.25 * sinI2 - 0.5);
  635.             l3 = qq * (77.0 / 24.0 * sinI2 - 7.0 / 6.0);
  636.             ll = g3 * (53.0 * sinI2 - 4.0 - 57.5 * sinI4) +
  637.                  2.5 * g5 * (4.0 - 96.0 * sinI2 + 269.5 * sinI4 - 183.75 * sinI6);

  639.         }

  641.         /** Extrapolate an orbit up to a specific target date.
  642.          * @param date target date for the orbit
  643.          * @return propagated parameters
  644.          * @exception PropagationException if some parameters are out of bounds
  645.          */
  646.         public DerivativeStructure[] propagateParameters(final AbsoluteDate date)
  647.             throws PropagationException {

  649.             // keplerian evolution
  650.             final DerivativeStructure dt =
  651.                     new DerivativeStructure(1, 2, 0, date.durationFrom(mean.getDate()));
  652.             final DerivativeStructure xnot = dt.multiply(xnotDot);

  654.             // secular effects

  656.             // eccentricity
  657.             final DerivativeStructure x   = xnot.multiply(rdpom + rdpomp);
  658.             final DerivativeStructure cx  = x.cos();
  659.             final DerivativeStructure sx  = x.sin();
  660.             final DerivativeStructure exm = cx.multiply(mean.getCircularEx()).
  661.                                             add(sx.multiply(eps2 - (1.0 - eps1) * mean.getCircularEy()));
  662.             final DerivativeStructure eym = sx.multiply((1.0 + eps1) * mean.getCircularEx()).
  663.                                             add(cx.multiply(mean.getCircularEy() - eps2)).
  664.                                             add(eps2);

  666.             // no secular effect on inclination

  668.             // right ascension of ascending node
  669.             final DerivativeStructure omm =
  670.                     new DerivativeStructure(1, 2,
  671.                                             MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode() + ommD * xnot.getValue(),
  672.                                                                      FastMath.PI),
  673.                                             ommD * xnotDot,
  674.                                             0.0);

  676.             // latitude argument
  677.             final DerivativeStructure xlm =
  678.                     new DerivativeStructure(1, 2,
  679.                                             MathUtils.normalizeAngle(mean.getAlphaM() + aMD * xnot.getValue(), FastMath.PI),
  680.                                             aMD * xnotDot,
  681.                                             0.0);

  683.             // periodical terms
  684.             final DerivativeStructure cl1 = xlm.cos();
  685.             final DerivativeStructure sl1 = xlm.sin();
  686.             final DerivativeStructure cl2 = cl1.multiply(cl1).subtract(sl1.multiply(sl1));
  687.             final DerivativeStructure sl2 = cl1.multiply(sl1).add(sl1.multiply(cl1));
  688.             final DerivativeStructure cl3 = cl2.multiply(cl1).subtract(sl2.multiply(sl1));
  689.             final DerivativeStructure sl3 = cl2.multiply(sl1).add(sl2.multiply(cl1));
  690.             final DerivativeStructure cl4 = cl3.multiply(cl1).subtract(sl3.multiply(sl1));
  691.             final DerivativeStructure sl4 = cl3.multiply(sl1).add(sl3.multiply(cl1));
  692.             final DerivativeStructure cl5 = cl4.multiply(cl1).subtract(sl4.multiply(sl1));
  693.             final DerivativeStructure sl5 = cl4.multiply(sl1).add(sl4.multiply(cl1));
  694.             final DerivativeStructure cl6 = cl5.multiply(cl1).subtract(sl5.multiply(sl1));

  696.             final DerivativeStructure qh  = eym.subtract(eps2).multiply(kh);
  697.             final DerivativeStructure ql  = exm.multiply(kl);

  699.             final DerivativeStructure exmCl1 = exm.multiply(cl1);
  700.             final DerivativeStructure exmSl1 = exm.multiply(sl1);
  701.             final DerivativeStructure eymCl1 = eym.multiply(cl1);
  702.             final DerivativeStructure eymSl1 = eym.multiply(sl1);
  703.             final DerivativeStructure exmCl2 = exm.multiply(cl2);
  704.             final DerivativeStructure exmSl2 = exm.multiply(sl2);
  705.             final DerivativeStructure eymCl2 = eym.multiply(cl2);
  706.             final DerivativeStructure eymSl2 = eym.multiply(sl2);
  707.             final DerivativeStructure exmCl3 = exm.multiply(cl3);
  708.             final DerivativeStructure exmSl3 = exm.multiply(sl3);
  709.             final DerivativeStructure eymCl3 = eym.multiply(cl3);
  710.             final DerivativeStructure eymSl3 = eym.multiply(sl3);
  711.             final DerivativeStructure exmCl4 = exm.multiply(cl4);
  712.             final DerivativeStructure exmSl4 = exm.multiply(sl4);
  713.             final DerivativeStructure eymCl4 = eym.multiply(cl4);
  714.             final DerivativeStructure eymSl4 = eym.multiply(sl4);

  716.             // semi major axis
  717.             final DerivativeStructure rda = exmCl1.multiply(ax1).
  718.                                             add(eymSl1.multiply(ay1)).
  719.                                             add(sl1.multiply(as1)).
  720.                                             add(cl2.multiply(ac2)).
  721.                                             add(exmCl3.add(eymSl3).multiply(axy3)).
  722.                                             add(sl3.multiply(as3)).
  723.                                             add(cl4.multiply(ac4)).
  724.                                             add(sl5.multiply(as5)).
  725.                                             add(cl6.multiply(ac6));

  727.             // eccentricity
  728.             final DerivativeStructure rdex = cl1.multiply(ex1).
  729.                                              add(exmCl2.multiply(exx2)).
  730.                                              add(eymSl2.multiply(exy2)).
  731.                                              add(cl3.multiply(ex3)).
  732.                                              add(exmCl4.add(eymSl4).multiply(ex4));
  733.             final DerivativeStructure rdey = sl1.multiply(ey1).
  734.                                              add(exmSl2.multiply(eyx2)).
  735.                                              add(eymCl2.multiply(eyy2)).
  736.                                              add(sl3.multiply(ey3)).
  737.                                              add(exmSl4.subtract(eymCl4).multiply(ey4));

  739.             // ascending node
  740.             final DerivativeStructure rdom = exmSl1.multiply(rx1).
  741.                                              add(eymCl1.multiply(ry1)).
  742.                                              add(sl2.multiply(r2)).
  743.                                              add(eymCl3.subtract(exmSl3).multiply(r3)).
  744.                                              add(ql.multiply(rl));

  746.             // inclination
  747.             final DerivativeStructure rdxi = eymSl1.multiply(iy1).
  748.                                              add(exmCl1.multiply(ix1)).
  749.                                              add(cl2.multiply(i2)).
  750.                                              add(exmCl3.add(eymSl3).multiply(i3)).
  751.                                              add(qh.multiply(ih));

  753.             // latitude argument
  754.             final DerivativeStructure rdxl = exmSl1.multiply(lx1).
  755.                                              add(eymCl1.multiply(ly1)).
  756.                                              add(sl2.multiply(l2)).
  757.                                              add(exmSl3.subtract(eymCl3).multiply(l3)).
  758.                                              add(ql.multiply(ll));

  760.             // osculating parameters
  761.             return new DerivativeStructure[] {
  762.                 rda.add(1.0).multiply(mean.getA()),
  763.                 rdex.add(exm),
  764.                 rdey.add(eym),
  765.                 rdxi.add(xim),
  766.                 rdom.add(omm),
  767.                 rdxl.add(xlm)
  768.             };

  770.         }

  772.     }

  774.     /** Convert circular parameters <em>with derivatives</em> to Cartesian coordinates.
  775.      * @param date date of the orbital parameters
  776.      * @param parameters circular parameters (a, ex, ey, i, raan, alphaM)
  777.      * @return Cartesian coordinates consistent with values and derivatives
  778.      */
  779.     private TimeStampedPVCoordinates toCartesian(final AbsoluteDate date, final DerivativeStructure[] parameters) {

  781.         // evaluate coordinates in the orbit canonical reference frame
  782.         final DerivativeStructure cosOmega = parameters[4].cos();
  783.         final DerivativeStructure sinOmega = parameters[4].sin();
  784.         final DerivativeStructure cosI     = parameters[3].cos();
  785.         final DerivativeStructure sinI     = parameters[3].sin();
  786.         final DerivativeStructure alphaE   = meanToEccentric(parameters[5], parameters[1], parameters[2]);
  787.         final DerivativeStructure cosAE    = alphaE.cos();
  788.         final DerivativeStructure sinAE    = alphaE.sin();
  789.         final DerivativeStructure ex2      = parameters[1].multiply(parameters[1]);
  790.         final DerivativeStructure ey2      = parameters[2].multiply(parameters[2]);
  791.         final DerivativeStructure exy      = parameters[1].multiply(parameters[2]);
  792.         final DerivativeStructure q        = ex2.add(ey2).subtract(1).negate().sqrt();
  793.         final DerivativeStructure beta     = q.add(1).reciprocal();
  794.         final DerivativeStructure bx2      = beta.multiply(ex2);
  795.         final DerivativeStructure by2      = beta.multiply(ey2);
  796.         final DerivativeStructure bxy      = beta.multiply(exy);
  797.         final DerivativeStructure u        = bxy.multiply(sinAE).subtract(parameters[1].add(by2.subtract(1).multiply(cosAE)));
  798.         final DerivativeStructure v        = bxy.multiply(cosAE).subtract(parameters[2].add(bx2.subtract(1).multiply(sinAE)));
  799.         final DerivativeStructure x        = parameters[0].multiply(u);
  800.         final DerivativeStructure y        = parameters[0].multiply(v);

  802.         // canonical orbit reference frame
  803.         final FieldVector3D<DerivativeStructure> p =
  804.                 new FieldVector3D<DerivativeStructure>(x.multiply(cosOmega).subtract(y.multiply(cosI.multiply(sinOmega))),
  805.                                                        x.multiply(sinOmega).add(y.multiply(cosI.multiply(cosOmega))),
  806.                                                        y.multiply(sinI));

  808.         // dispatch derivatives
  809.         final Vector3D p0 = new Vector3D(p.getX().getValue(),
  810.                                          p.getY().getValue(),
  811.                                          p.getZ().getValue());
  812.         final Vector3D p1 = new Vector3D(p.getX().getPartialDerivative(1),
  813.                                          p.getY().getPartialDerivative(1),
  814.                                          p.getZ().getPartialDerivative(1));
  815.         final Vector3D p2 = new Vector3D(p.getX().getPartialDerivative(2),
  816.                                          p.getY().getPartialDerivative(2),
  817.                                          p.getZ().getPartialDerivative(2));
  818.         return new TimeStampedPVCoordinates(date, p0, p1, p2);

  820.     }

  822.     /** Computes the eccentric latitude argument from the mean latitude argument.
  823.      * @param alphaM = M + Ω mean latitude argument (rad)
  824.      * @param ex e cos(Ω), first component of circular eccentricity vector
  825.      * @param ey e sin(Ω), second component of circular eccentricity vector
  826.      * @return the eccentric latitude argument.
  827.      */
  828.     private DerivativeStructure meanToEccentric(final DerivativeStructure alphaM,
  829.                                                 final DerivativeStructure ex,
  830.                                                 final DerivativeStructure ey) {
  831.         // Generalization of Kepler equation to circular parameters
  832.         // with alphaE = PA + E and
  833.         //      alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
  834.         DerivativeStructure alphaE        = alphaM;
  835.         DerivativeStructure shift         = alphaM.getField().getZero();
  836.         DerivativeStructure alphaEMalphaM = alphaM.getField().getZero();
  837.         DerivativeStructure cosAlphaE     = alphaE.cos();
  838.         DerivativeStructure sinAlphaE     = alphaE.sin();
  839.         int                 iter          = 0;
  840.         do {
  841.             final DerivativeStructure f2 = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
  842.             final DerivativeStructure f1 = alphaM.getField().getOne().subtract(ex.multiply(cosAlphaE)).subtract(ey.multiply(sinAlphaE));
  843.             final DerivativeStructure f0 = alphaEMalphaM.subtract(f2);

  845.             final DerivativeStructure f12 = f1.multiply(2);
  846.             shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));

  848.             alphaEMalphaM  = alphaEMalphaM.subtract(shift);
  849.             alphaE         = alphaM.add(alphaEMalphaM);
  850.             cosAlphaE      = alphaE.cos();
  851.             sinAlphaE      = alphaE.sin();

  853.         } while ((++iter < 50) && (FastMath.abs(shift.getValue()) > 1.0e-12));

  855.         return alphaE;

  857.     }

  859.     /** {@inheritDoc} */
  860.     protected double getMass(final AbsoluteDate date) {
  861.         return models.get(date).mass;
  862.     }

  864.     /** Replace the instance with a data transfer object for serialization.
  865.      * @return data transfer object that will be serialized
  866.      * @exception NotSerializableException if an additional state provider is not serializable
  867.      */
  868.     private Object writeReplace() throws NotSerializableException {
  869.         try {
  870.             // managed states providers
  871.             final List<AdditionalStateProvider> serializableProviders = new ArrayList<AdditionalStateProvider>();
  872.             for (final AdditionalStateProvider provider : getAdditionalStateProviders()) {
  873.                 if (provider instanceof Serializable) {
  874.                     serializableProviders.add(provider);
  875.                 } else {
  876.                     throw new NotSerializableException(provider.getClass().getName());
  877.                 }
  878.             }

  880.             // states transitions
  881.             final AbsoluteDate[]  transitionDates;
  882.             final CircularOrbit[] allOrbits;
  883.             final double[]        allMasses;
  884.             final SortedSet<TimeSpanMap.Transition<EHModel>> transitions = models.getTransitions();
  885.             if (transitions.size() == 1  && transitions.first().getBefore() == transitions.first().getAfter()) {
  886.                 // the single entry is a dummy one, without a real transition
  887.                 // we ignore it completely
  888.                 transitionDates = null;
  889.                 allOrbits       = null;
  890.                 allMasses       = null;
  891.             } else {
  892.                 transitionDates = new AbsoluteDate[transitions.size()];
  893.                 allOrbits       = new CircularOrbit[transitions.size() + 1];
  894.                 allMasses       = new double[transitions.size() + 1];
  895.                 int i = 0;
  896.                 for (final TimeSpanMap.Transition<EHModel> transition : transitions) {
  897.                     if (i == 0) {
  898.                         // model before the first transition
  899.                         allOrbits[i] = transition.getBefore().mean;
  900.                         allMasses[i] = transition.getBefore().mass;
  901.                     }
  902.                     transitionDates[i] = transition.getDate();
  903.                     allOrbits[++i]     = transition.getAfter().mean;
  904.                     allMasses[i]       = transition.getAfter().mass;
  905.                 }
  906.             }

  908.             return new DataTransferObject(getInitialState().getOrbit(), initialModel.mass,
  909.                                           referenceRadius, mu, ck0, getAttitudeProvider(),
  910.                                           transitionDates, allOrbits, allMasses,
  911.                                           serializableProviders.toArray(new AdditionalStateProvider[serializableProviders.size()]));
  912.         } catch (OrekitException orekitException) {
  913.             // this should never happen
  914.             throw new OrekitInternalError(null);
  915.         }

  917.     }

  919.     /** Internal class used only for serialization. */
  920.     private static class DataTransferObject implements Serializable {

  922.         /** Serializable UID. */
  923.         private static final long serialVersionUID = 20151202L;

  925.         /** Initial orbit. */
  926.         private final Orbit orbit;

  928.         /** Attitude provider. */
  929.         private final AttitudeProvider attitudeProvider;

  931.         /** Mass and gravity field. */
  932.         private double[] g;

  934.         /** Transition dates (may be null). */
  935.         private final AbsoluteDate[] transitionDates;

  937.         /** Orbits before and after transitions (may be null). */
  938.         private final CircularOrbit[] allOrbits;

  940.         /** Masses before and after transitions (may be null). */
  941.         private final double[] allMasses;

  943.         /** Providers for additional states. */
  944.         private final AdditionalStateProvider[] providers;

  946.         /** Simple constructor.
  947.          * @param orbit initial orbit
  948.          * @param mass spacecraft mass
  949.          * @param referenceRadius reference radius of the Earth for the potential model (m)
  950.          * @param mu central attraction coefficient (m³/s²)
  951.          * @param ck0 un-normalized zonal coefficients
  952.          * @param attitudeProvider attitude provider
  953.          * @param transitionDates transition dates (may be null)
  954.          * @param allOrbits orbits before and after transitions (may be null)
  955.          * @param allMasses masses before and after transitions (may be null)
  956.          * @param providers providers for additional states
  957.          */
  958.         DataTransferObject(final Orbit orbit, final double mass,
  959.                            final double referenceRadius, final double mu,
  960.                            final double[] ck0,
  961.                            final AttitudeProvider attitudeProvider,
  962.                            final AbsoluteDate[] transitionDates,
  963.                            final CircularOrbit[] allOrbits,
  964.                            final double[] allMasses,
  965.                            final AdditionalStateProvider[] providers) {
  966.             this.orbit            = orbit;
  967.             this.attitudeProvider = attitudeProvider;
  968.             this.g = new double[] {
  969.                 mass, referenceRadius, mu,
  970.                 ck0[2], ck0[3], ck0[4], ck0[5], ck0[6] // ck0[0] and ck0[1] are both zero so not serialized
  971.             };
  972.             this.transitionDates  = transitionDates;
  973.             this.allOrbits        = allOrbits;
  974.             this.allMasses        = allMasses;
  975.             this.providers        = providers;
  976.         }

  978.         /** Replace the deserialized data transfer object with a {@link EcksteinHechlerPropagator}.
  979.          * @return replacement {@link EcksteinHechlerPropagator}
  980.          */
  981.         private Object readResolve() {
  982.             try {
  983.                 final EcksteinHechlerPropagator propagator =
  984.                                 new EcksteinHechlerPropagator(orbit, attitudeProvider,
  985.                                                               g[0], g[1], g[2],              // mass, referenceRadius, mu
  986.                                                               g[3], g[4], g[5], g[6], g[7]); // c20, c30, c40, c50, c60
  987.                 for (final AdditionalStateProvider provider : providers) {
  988.                     propagator.addAdditionalStateProvider(provider);

  990.                 }
  991.                 if (transitionDates != null) {
  992.                     // override the state transitions
  993.                     final double[] ck0 = new double[] {
  994.                         0, 0, g[3], g[4], g[5], g[6], g[7]
  995.                     };
  996.                     propagator.models = new TimeSpanMap<EHModel>(new EHModel(allOrbits[0], allMasses[0],
  997.                                                                              g[1], g[2], ck0));
  998.                     for (int i = 0; i < transitionDates.length; ++i) {
  999.                         propagator.models.addValidAfter(new EHModel(allOrbits[i + 1], allMasses[i + 1],
  1000.                                                                     g[1], g[2], ck0),
  1001.                                                         transitionDates[i]);
  1002.                     }
  1003.                 }

  1005.                 return propagator;
  1006.             } catch (OrekitException oe) {
  1007.                 throw new OrekitInternalError(oe);
  1008.             }
  1009.         }

  1011.     }

  1013. }
Created with JaCoCo 0.7.4.201502262128