BrouwerLyddanePropagator.java

  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.propagation.analytical;

  19. import java.util.ArrayList;
  20. import java.util.Collections;
  21. import java.util.List;

  23. import org.hipparchus.analysis.differentiation.UnivariateDerivative1;
  24. import org.hipparchus.linear.RealMatrix;
  25. import org.hipparchus.util.CombinatoricsUtils;
  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.FieldSinCos;
  28. import org.hipparchus.util.MathUtils;
  29. import org.hipparchus.util.SinCos;
  30. import org.orekit.attitudes.AttitudeProvider;
  31. import org.orekit.attitudes.FrameAlignedProvider;
  32. import org.orekit.errors.OrekitException;
  33. import org.orekit.errors.OrekitMessages;
  34. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  35. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  36. import org.orekit.orbits.FieldKeplerianAnomalyUtility;
  37. import org.orekit.orbits.KeplerianOrbit;
  38. import org.orekit.orbits.Orbit;
  39. import org.orekit.orbits.OrbitType;
  40. import org.orekit.orbits.PositionAngleType;
  41. import org.orekit.propagation.AbstractMatricesHarvester;
  42. import org.orekit.propagation.MatricesHarvester;
  43. import org.orekit.propagation.PropagationType;
  44. import org.orekit.propagation.SpacecraftState;
  45. import org.orekit.propagation.analytical.tle.TLE;
  46. import org.orekit.time.AbsoluteDate;
  47. import org.orekit.utils.DoubleArrayDictionary;
  48. import org.orekit.utils.ParameterDriver;
  49. import org.orekit.utils.ParameterDriversProvider;
  50. import org.orekit.utils.TimeSpanMap;
  51. import org.orekit.utils.TimeSpanMap.Span;

  53. /**
  54.  * This class propagates a {@link org.orekit.propagation.SpacecraftState}
  55.  *  using the analytical Brouwer-Lyddane model (from J2 to J5 zonal harmonics).
  56.  * <p>
  57.  * At the opposite of the {@link EcksteinHechlerPropagator}, the Brouwer-Lyddane model is
  58.  * suited for elliptical orbits, there is no problem having a rather small eccentricity or inclination
  59.  * (Lyddane helped to solve this issue with the Brouwer model). Singularity for the critical
  60.  * inclination i = 63.4° is avoided using the method developed in Warren Phipps' 1992 thesis.
  61.  * <p>
  62.  * By default, Brouwer-Lyddane model considers only the perturbations due to zonal harmonics.
  63.  * However, for low Earth orbits, the magnitude of the perturbative acceleration due to
  64.  * atmospheric drag can be significant. Warren Phipps' 1992 thesis considered the atmospheric
  65.  * drag by time derivatives of the <i>mean</i> mean anomaly using the catch-all coefficient
  66.  * {@link #M2Driver}. Beware that M2Driver must have only 1 span on its TimeSpanMap value.
  67.  *
  68.  * Usually, M2 is adjusted during an orbit determination process and it represents the
  69.  * combination of all unmodeled secular along-track effects (i.e. not just the atmospheric drag).
  70.  * The behavior of M2 is close to the {@link TLE#getBStar()} parameter for the TLE.
  71.  *
  72.  * If the value of M2 is equal to {@link #M2 0.0}, the along-track  secular effects are not
  73.  * considered in the dynamical model. Typical values for M2 are not known. It depends on the
  74.  * orbit type. However, the value of M2 must be very small (e.g. between 1.0e-14 and 1.0e-15).
  75.  * The unit of M2 is rad/s².
  76.  *
  77.  * The along-track effects, represented by the secular rates of the mean semi-major axis
  78.  * and eccentricity, are computed following Eq. 2.38, 2.41, and 2.45 of Warren Phipps' thesis.
  79.  *
  80.  * @see "Brouwer, Dirk. Solution of the problem of artificial satellite theory without drag.
  81.  *       YALE UNIV NEW HAVEN CT NEW HAVEN United States, 1959."
  82.  *
  83.  * @see "Lyddane, R. H. Small eccentricities or inclinations in the Brouwer theory of the
  84.  *       artificial satellite. The Astronomical Journal 68 (1963): 555."
  85.  *
  86.  * @see "Phipps Jr, Warren E. Parallelization of the Navy Space Surveillance Center
  87.  *       (NAVSPASUR) Satellite Model. NAVAL POSTGRADUATE SCHOOL MONTEREY CA, 1992."
  88.  *
  89.  * @author Melina Vanel
  90.  * @author Bryan Cazabonne
  91.  * @since 11.1
  92.  */
  93. public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator implements ParameterDriversProvider {

  95.     /** Parameter name for M2 coefficient. */
  96.     public static final String M2_NAME = "M2";

  98.     /** Default value for M2 coefficient. */
  99.     public static final double M2 = 0.0;

  101.     /** Default convergence threshold for mean parameters conversion. */
  102.     private static final double EPSILON_DEFAULT = 1.0e-13;

  104.     /** Default value for maxIterations. */
  105.     private static final int MAX_ITERATIONS_DEFAULT = 200;

  107.     /** Parameters scaling factor.
  108.      * <p>
  109.      * We use a power of 2 to avoid numeric noise introduction
  110.      * in the multiplications/divisions sequences.
  111.      * </p>
  112.      */
  113.     private static final double SCALE = FastMath.scalb(1.0, -32);

  115.     /** Beta constant used by T2 function. */
  116.     private static final double BETA = FastMath.scalb(100, -11);

  118.     /** Initial Brouwer-Lyddane model. */
  119.     private BLModel initialModel;

  121.     /** All models. */
  122.     private transient TimeSpanMap<BLModel> models;

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

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

  130.     /** Un-normalized zonal coefficients. */
  131.     private double[] ck0;

  133.     /** Empirical coefficient used in the drag modeling. */
  134.     private final ParameterDriver M2Driver;

  136.     /** Build a propagator from orbit and potential provider.
  137.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  138.      *
  139.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  140.      *
  141.      * @param initialOrbit initial orbit
  142.      * @param provider for un-normalized zonal coefficients
  143.      * @param M2 value of empirical drag coefficient in rad/s².
  144.      *        If equal to {@link #M2} drag is not computed
  145.      * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, UnnormalizedSphericalHarmonicsProvider, double)
  146.      * @see #BrouwerLyddanePropagator(Orbit, UnnormalizedSphericalHarmonicsProvider, PropagationType, double)
  147.      */
  148.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  149.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  150.                                     final double M2) {
  151.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  152.              DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), M2);
  153.     }

  155.     /**
  156.      * Private helper constructor.
  157.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  158.      * @param initialOrbit initial orbit
  159.      * @param attitude attitude provider
  160.      * @param mass spacecraft mass
  161.      * @param provider for un-normalized zonal coefficients
  162.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  163.      * @param M2 value of empirical drag coefficient in rad/s².
  164.      *        If equal to {@link #M2} drag is not computed
  165.      * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double,
  166.      *                                 UnnormalizedSphericalHarmonicsProvider,
  167.      *                                 UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics,
  168.      *                                 PropagationType, double)
  169.      */
  170.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  171.                                     final AttitudeProvider attitude,
  172.                                     final double mass,
  173.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  174.                                     final UnnormalizedSphericalHarmonics harmonics,
  175.                                     final double M2) {
  176.         this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
  177.              harmonics.getUnnormalizedCnm(2, 0),
  178.              harmonics.getUnnormalizedCnm(3, 0),
  179.              harmonics.getUnnormalizedCnm(4, 0),
  180.              harmonics.getUnnormalizedCnm(5, 0),
  181.              M2);
  182.     }

  184.     /** Build a propagator from orbit and potential.
  185.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  186.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  187.      * are related to both the normalized coefficients
  188.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  189.      *  and the J<sub>n</sub> one as follows:</p>
  190.      *
  191.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  192.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  193.      *
  194.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  195.      *
  196.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  197.      *
  198.      * @param initialOrbit initial orbit
  199.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  200.      * @param mu central attraction coefficient (m³/s²)
  201.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  202.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  203.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  204.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  205.      * @param M2 value of empirical drag coefficient in rad/s².
  206.      *        If equal to {@link #M2} drag is not computed
  207.      * @see org.orekit.utils.Constants
  208.      * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, double, double,
  209.      * double, double, double, double, double)
  210.      */
  211.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  212.                                     final double referenceRadius, final double mu,
  213.                                     final double c20, final double c30, final double c40,
  214.                                     final double c50, final double M2) {
  215.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  216.              DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, M2);
  217.     }

  219.     /** Build a propagator from orbit, mass and potential provider.
  220.      * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
  221.      *
  222.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  223.      *
  224.      * @param initialOrbit initial orbit
  225.      * @param mass spacecraft mass
  226.      * @param provider for un-normalized zonal coefficients
  227.      * @param M2 value of empirical drag coefficient in rad/s².
  228.      *        If equal to {@link #M2} drag is not computed
  229.      * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, UnnormalizedSphericalHarmonicsProvider, double)
  230.      */
  231.     public BrouwerLyddanePropagator(final Orbit initialOrbit, final double mass,
  232.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  233.                                     final double M2) {
  234.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  235.              mass, provider, provider.onDate(initialOrbit.getDate()), M2);
  236.     }

  238.     /** Build a propagator from orbit, mass and potential.
  239.      * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
  240.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  241.      * are related to both the normalized coefficients
  242.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  243.      *  and the J<sub>n</sub> one as follows:</p>
  244.      *
  245.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  246.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  247.      *
  248.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  249.      *
  250.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  251.      *
  252.      * @param initialOrbit initial orbit
  253.      * @param mass spacecraft mass
  254.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  255.      * @param mu central attraction coefficient (m³/s²)
  256.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  257.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  258.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  259.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  260.      * @param M2 value of empirical drag coefficient in rad/s².
  261.      *        If equal to {@link #M2} drag is not computed
  262.      * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, double, double,
  263.      * double, double, double, double, double)
  264.      */
  265.     public BrouwerLyddanePropagator(final Orbit initialOrbit, final double mass,
  266.                                     final double referenceRadius, final double mu,
  267.                                     final double c20, final double c30, final double c40,
  268.                                     final double c50, final double M2) {
  269.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  270.              mass, referenceRadius, mu, c20, c30, c40, c50, M2);
  271.     }

  273.     /** Build a propagator from orbit, attitude provider and potential provider.
  274.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  275.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  276.      * @param initialOrbit initial orbit
  277.      * @param attitudeProv attitude provider
  278.      * @param provider for un-normalized zonal coefficients
  279.      * @param M2 value of empirical drag coefficient in rad/s².
  280.      *        If equal to {@link #M2} drag is not computed
  281.      */
  282.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  283.                                     final AttitudeProvider attitudeProv,
  284.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  285.                                     final double M2) {
  286.         this(initialOrbit, attitudeProv, DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), M2);
  287.     }

  289.     /** Build a propagator from orbit, attitude provider and potential.
  290.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  291.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  292.      * are related to both the normalized coefficients
  293.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  294.      *  and the J<sub>n</sub> one as follows:</p>
  295.      *
  296.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  297.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  298.      *
  299.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  300.      *
  301.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  302.      *
  303.      * @param initialOrbit initial orbit
  304.      * @param attitudeProv attitude provider
  305.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  306.      * @param mu central attraction coefficient (m³/s²)
  307.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  308.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  309.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  310.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  311.      * @param M2 value of empirical drag coefficient in rad/s².
  312.      *        If equal to {@link #M2} drag is not computed
  313.      */
  314.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  315.                                     final AttitudeProvider attitudeProv,
  316.                                     final double referenceRadius, final double mu,
  317.                                     final double c20, final double c30, final double c40,
  318.                                     final double c50, final double M2) {
  319.         this(initialOrbit, attitudeProv, DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, M2);
  320.     }

  322.     /** Build a propagator from orbit, attitude provider, mass and potential provider.
  323.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  324.      * @param initialOrbit initial orbit
  325.      * @param attitudeProv attitude provider
  326.      * @param mass spacecraft mass
  327.      * @param provider for un-normalized zonal coefficients
  328.      * @param M2 value of empirical drag coefficient in rad/s².
  329.      *        If equal to {@link #M2} drag is not computed
  330.      * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double,
  331.      *                                UnnormalizedSphericalHarmonicsProvider, PropagationType, double)
  332.      */
  333.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  334.                                     final AttitudeProvider attitudeProv,
  335.                                     final double mass,
  336.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  337.                                     final double M2) {
  338.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), M2);
  339.     }

  341.     /** Build a propagator from orbit, attitude provider, mass and potential.
  342.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  343.      * are related to both the normalized coefficients
  344.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  345.      *  and the J<sub>n</sub> one as follows:</p>
  346.      *
  347.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  348.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  349.      *
  350.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  351.      *
  352.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  353.      *
  354.      * @param initialOrbit initial orbit
  355.      * @param attitudeProv attitude provider
  356.      * @param mass spacecraft mass
  357.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  358.      * @param mu central attraction coefficient (m³/s²)
  359.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  360.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  361.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  362.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  363.      * @param M2 value of empirical drag coefficient in rad/s².
  364.      *        If equal to {@link #M2} drag is not computed
  365.      * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, double, double,
  366.      *                                 double, double, double, double, PropagationType, double)
  367.      */
  368.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  369.                                     final AttitudeProvider attitudeProv,
  370.                                     final double mass,
  371.                                     final double referenceRadius, final double mu,
  372.                                     final double c20, final double c30, final double c40,
  373.                                     final double c50, final double M2) {
  374.         this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50,
  375.              PropagationType.OSCULATING, M2);
  376.     }


  379.     /** Build a propagator from orbit and potential provider.
  380.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  381.      *
  382.      * <p>Using this constructor, it is possible to define the initial orbit as
  383.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  384.      *
  385.      * @param initialOrbit initial orbit
  386.      * @param provider for un-normalized zonal coefficients
  387.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  388.      * @param M2 value of empirical drag coefficient in rad/s².
  389.      *        If equal to {@link #M2} drag is not computed
  390.      */
  391.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  392.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  393.                                     final PropagationType initialType, final double M2) {
  394.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  395.              DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), initialType, M2);
  396.     }

  398.     /** Build a propagator from orbit, attitude provider, mass and potential provider.
  399.      * <p>Using this constructor, it is possible to define the initial orbit as
  400.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  401.      * @param initialOrbit initial orbit
  402.      * @param attitudeProv attitude provider
  403.      * @param mass spacecraft mass
  404.      * @param provider for un-normalized zonal coefficients
  405.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  406.      * @param M2 value of empirical drag coefficient in rad/s².
  407.      *        If equal to {@link #M2} drag is not computed
  408.      */
  409.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  410.                                     final AttitudeProvider attitudeProv,
  411.                                     final double mass,
  412.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  413.                                     final PropagationType initialType, final double M2) {
  414.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), initialType, M2);
  415.     }

  417.     /**
  418.      * Private helper constructor.
  419.      * <p>Using this constructor, it is possible to define the initial orbit as
  420.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  421.      * @param initialOrbit initial orbit
  422.      * @param attitude attitude provider
  423.      * @param mass spacecraft mass
  424.      * @param provider for un-normalized zonal coefficients
  425.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  426.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  427.      * @param M2 value of empirical drag coefficient in rad/s².
  428.      *        If equal to {@link #M2} drag is not computed
  429.      */
  430.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  431.                                     final AttitudeProvider attitude,
  432.                                     final double mass,
  433.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  434.                                     final UnnormalizedSphericalHarmonics harmonics,
  435.                                     final PropagationType initialType, final double M2) {
  436.         this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
  437.              harmonics.getUnnormalizedCnm(2, 0),
  438.              harmonics.getUnnormalizedCnm(3, 0),
  439.              harmonics.getUnnormalizedCnm(4, 0),
  440.              harmonics.getUnnormalizedCnm(5, 0),
  441.              initialType, M2);
  442.     }

  444.     /** Build a propagator from orbit, attitude provider, mass and potential.
  445.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  446.      * are related to both the normalized coefficients
  447.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  448.      *  and the J<sub>n</sub> one as follows:</p>
  449.      *
  450.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  451.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  452.      *
  453.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  454.      *
  455.      * <p>Using this constructor, it is possible to define the initial orbit as
  456.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  457.      *
  458.      * @param initialOrbit initial orbit
  459.      * @param attitudeProv attitude provider
  460.      * @param mass spacecraft mass
  461.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  462.      * @param mu central attraction coefficient (m³/s²)
  463.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  464.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  465.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  466.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  467.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  468.      * @param M2 value of empirical drag coefficient in rad/s².
  469.      *        If equal to {@link #M2} drag is not computed
  470.      */
  471.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  472.                                     final AttitudeProvider attitudeProv,
  473.                                     final double mass,
  474.                                     final double referenceRadius, final double mu,
  475.                                     final double c20, final double c30, final double c40,
  476.                                     final double c50,
  477.                                     final PropagationType initialType, final double M2) {
  478.         this(initialOrbit, attitudeProv, mass, referenceRadius, mu,
  479.              c20, c30, c40, c50, initialType, M2, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  480.     }

  482.     /** Build a propagator from orbit, attitude provider, mass and potential.
  483.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  484.      * are related to both the normalized coefficients
  485.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  486.      *  and the J<sub>n</sub> one as follows:</p>
  487.      *
  488.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  489.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  490.      *
  491.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  492.      *
  493.      * <p>Using this constructor, it is possible to define the initial orbit as
  494.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  495.      *
  496.      * @param initialOrbit initial orbit
  497.      * @param attitudeProv attitude provider
  498.      * @param mass spacecraft mass
  499.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  500.      * @param mu central attraction coefficient (m³/s²)
  501.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  502.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  503.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  504.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  505.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  506.      * @param M2 value of empirical drag coefficient in rad/s².
  507.      *        If equal to {@link #M2} drag is not computed
  508.      * @param epsilon convergence threshold for mean parameters conversion
  509.      * @param maxIterations maximum iterations for mean parameters conversion
  510.      * @since 11.2
  511.      */
  512.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  513.                                     final AttitudeProvider attitudeProv,
  514.                                     final double mass,
  515.                                     final double referenceRadius, final double mu,
  516.                                     final double c20, final double c30, final double c40,
  517.                                     final double c50,
  518.                                     final PropagationType initialType, final double M2,
  519.                                     final double epsilon, final int maxIterations) {

  521.         super(attitudeProv);

  523.         // store model coefficients
  524.         this.referenceRadius = referenceRadius;
  525.         this.mu  = mu;
  526.         this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50};

  528.         // initialize M2 driver
  529.         this.M2Driver = new ParameterDriver(M2_NAME, M2, SCALE,
  530.                                             Double.NEGATIVE_INFINITY,
  531.                                             Double.POSITIVE_INFINITY);

  533.         // compute mean parameters if needed
  534.         resetInitialState(new SpacecraftState(initialOrbit,
  535.                                               attitudeProv.getAttitude(initialOrbit,
  536.                                                                        initialOrbit.getDate(),
  537.                                                                        initialOrbit.getFrame()),
  538.                                               mass),
  539.                           initialType, epsilon, maxIterations);

  541.     }

  543.     /** Conversion from osculating to mean orbit.
  544.      * <p>
  545.      * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
  546.      * osculating SpacecraftState in input.
  547.      * </p>
  548.      * <p>
  549.      * Since the osculating orbit is obtained with the computation of
  550.      * short-periodic variation, the resulting output will depend on
  551.      * both the gravity field parameterized in input and the
  552.      * atmospheric drag represented by the {@code m2} parameter.
  553.      * </p>
  554.      * <p>
  555.      * The computation is done through a fixed-point iteration process.
  556.      * </p>
  557.      * @param osculating osculating orbit to convert
  558.      * @param provider for un-normalized zonal coefficients
  559.      * @param harmonics {@code provider.onDate(osculating.getDate())}
  560.      * @param M2Value value of empirical drag coefficient in rad/s².
  561.      *        If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
  562.      * @return mean orbit in a Brouwer-Lyddane sense
  563.      * @since 11.2
  564.      */
  565.     public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
  566.                                                   final UnnormalizedSphericalHarmonicsProvider provider,
  567.                                                   final UnnormalizedSphericalHarmonics harmonics,
  568.                                                   final double M2Value) {
  569.         return computeMeanOrbit(osculating, provider, harmonics, M2Value, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  570.     }

  572.     /** Conversion from osculating to mean orbit.
  573.      * <p>
  574.      * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
  575.      * osculating SpacecraftState in input.
  576.      * </p>
  577.      * <p>
  578.      * Since the osculating orbit is obtained with the computation of
  579.      * short-periodic variation, the resulting output will depend on
  580.      * both the gravity field parameterized in input and the
  581.      * atmospheric drag represented by the {@code m2} parameter.
  582.      * </p>
  583.      * <p>
  584.      * The computation is done through a fixed-point iteration process.
  585.      * </p>
  586.      * @param osculating osculating orbit to convert
  587.      * @param provider for un-normalized zonal coefficients
  588.      * @param harmonics {@code provider.onDate(osculating.getDate())}
  589.      * @param M2Value value of empirical drag coefficient in rad/s².
  590.      *        If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
  591.      * @param epsilon convergence threshold for mean parameters conversion
  592.      * @param maxIterations maximum iterations for mean parameters conversion
  593.      * @return mean orbit in a Brouwer-Lyddane sense
  594.      * @since 11.2
  595.      */
  596.     public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
  597.                                                   final UnnormalizedSphericalHarmonicsProvider provider,
  598.                                                   final UnnormalizedSphericalHarmonics harmonics,
  599.                                                   final double M2Value,
  600.                                                   final double epsilon, final int maxIterations) {
  601.         return computeMeanOrbit(osculating,
  602.                                 provider.getAe(), provider.getMu(),
  603.                                 harmonics.getUnnormalizedCnm(2, 0),
  604.                                 harmonics.getUnnormalizedCnm(3, 0),
  605.                                 harmonics.getUnnormalizedCnm(4, 0),
  606.                                 harmonics.getUnnormalizedCnm(5, 0),
  607.                                 M2Value, epsilon, maxIterations);
  608.     }

  610.     /** Conversion from osculating to mean orbit.
  611.      * <p>
  612.      * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
  613.      * osculating SpacecraftState in input.
  614.      * </p>
  615.      * <p>
  616.      * Since the osculating orbit is obtained with the computation of
  617.      * short-periodic variation, the resulting output will depend on
  618.      * both the gravity field parameterized in input and the
  619.      * atmospheric drag represented by the {@code m2} parameter.
  620.      * </p>
  621.      * <p>
  622.      * The computation is done through a fixed-point iteration process.
  623.      * </p>
  624.      * @param osculating osculating orbit to convert
  625.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  626.      * @param mu central attraction coefficient (m³/s²)
  627.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  628.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  629.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  630.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  631.      * @param M2Value value of empirical drag coefficient in rad/s².
  632.      *        If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
  633.      * @param epsilon convergence threshold for mean parameters conversion
  634.      * @param maxIterations maximum iterations for mean parameters conversion
  635.      * @return mean orbit in a Brouwer-Lyddane sense
  636.      * @since 11.2
  637.      */
  638.     public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
  639.                                                   final double referenceRadius, final double mu,
  640.                                                   final double c20, final double c30, final double c40,
  641.                                                   final double c50, final double M2Value,
  642.                                                   final double epsilon, final int maxIterations) {
  643.         final BrouwerLyddanePropagator propagator =
  644.                         new BrouwerLyddanePropagator(osculating,
  645.                                                      FrameAlignedProvider.of(osculating.getFrame()),
  646.                                                      DEFAULT_MASS,
  647.                                                      referenceRadius, mu, c20, c30, c40, c50,
  648.                                                      PropagationType.OSCULATING, M2Value,
  649.                                                      epsilon, maxIterations);
  650.         return propagator.initialModel.mean;
  651.     }

  653.     /** {@inheritDoc}
  654.      * <p>The new initial state to consider
  655.      * must be defined with an osculating orbit.</p>
  656.      * @see #resetInitialState(SpacecraftState, PropagationType)
  657.      */
  658.     @Override
  659.     public void resetInitialState(final SpacecraftState state) {
  660.         resetInitialState(state, PropagationType.OSCULATING);
  661.     }

  663.     /** Reset the propagator initial state.
  664.      * @param state new initial state to consider
  665.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  666.      */
  667.     public void resetInitialState(final SpacecraftState state, final PropagationType stateType) {
  668.         resetInitialState(state, stateType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  669.     }

  671.     /** Reset the propagator initial state.
  672.      * @param state new initial state to consider
  673.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  674.      * @param epsilon convergence threshold for mean parameters conversion
  675.      * @param maxIterations maximum iterations for mean parameters conversion
  676.      * @since 11.2
  677.      */
  678.     public void resetInitialState(final SpacecraftState state, final PropagationType stateType,
  679.                                   final double epsilon, final int maxIterations) {
  680.         super.resetInitialState(state);
  681.         final KeplerianOrbit keplerian = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state.getOrbit());
  682.         this.initialModel = (stateType == PropagationType.MEAN) ?
  683.                              new BLModel(keplerian, state.getMass(), referenceRadius, mu, ck0) :
  684.                              computeMeanParameters(keplerian, state.getMass(), epsilon, maxIterations);
  685.         this.models = new TimeSpanMap<>(initialModel);
  686.     }

  688.     /** {@inheritDoc} */
  689.     protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
  690.         resetIntermediateState(state, forward, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  691.     }

  693.     /** Reset an intermediate state.
  694.      * @param state new intermediate state to consider
  695.      * @param forward if true, the intermediate state is valid for
  696.      * propagations after itself
  697.      * @param epsilon convergence threshold for mean parameters conversion
  698.      * @param maxIterations maximum iterations for mean parameters conversion
  699.      * @since 11.2
  700.      */
  701.     protected void resetIntermediateState(final SpacecraftState state, final boolean forward,
  702.                                           final double epsilon, final int maxIterations) {
  703.         final BLModel newModel = computeMeanParameters((KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state.getOrbit()),
  704.                                                        state.getMass(), epsilon, maxIterations);
  705.         if (forward) {
  706.             models.addValidAfter(newModel, state.getDate(), false);
  707.         } else {
  708.             models.addValidBefore(newModel, state.getDate(), false);
  709.         }
  710.         stateChanged(state);
  711.     }

  713.     /** Compute mean parameters according to the Brouwer-Lyddane analytical model computation
  714.      * in order to do the propagation.
  715.      * @param osculating osculating orbit
  716.      * @param mass constant mass
  717.      * @param epsilon convergence threshold for mean parameters conversion
  718.      * @param maxIterations maximum iterations for mean parameters conversion
  719.      * @return Brouwer-Lyddane mean model
  720.      */
  721.     private BLModel computeMeanParameters(final KeplerianOrbit osculating, final double mass,
  722.                                           final double epsilon, final int maxIterations) {

  724.         // sanity check
  725.         if (osculating.getA() < referenceRadius) {
  726.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
  727.                                            osculating.getA());
  728.         }

  730.         // rough initialization of the mean parameters
  731.         BLModel current = new BLModel(osculating, mass, referenceRadius, mu, ck0);

  733.         // threshold for each parameter
  734.         final double thresholdA      = epsilon * (1 + FastMath.abs(current.mean.getA()));
  735.         final double thresholdE      = epsilon * (1 + current.mean.getE());
  736.         final double thresholdAngles = epsilon * FastMath.PI;

  738.         int i = 0;
  739.         while (i++ < maxIterations) {

  741.             // recompute the osculating parameters from the current mean parameters
  742.             final KeplerianOrbit parameters = current.propagateParameters(current.mean.getDate());

  744.             // adapted parameters residuals
  745.             final double deltaA     = osculating.getA() - parameters.getA();
  746.             final double deltaE     = osculating.getE() - parameters.getE();
  747.             final double deltaI     = osculating.getI() - parameters.getI();
  748.             final double deltaOmega = MathUtils.normalizeAngle(osculating.getPerigeeArgument() -
  749.                                                                parameters.getPerigeeArgument(),
  750.                                                                0.0);
  751.             final double deltaRAAN  = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode() -
  752.                                                                parameters.getRightAscensionOfAscendingNode(),
  753.                                                                0.0);
  754.             final double deltaAnom = MathUtils.normalizeAngle(osculating.getMeanAnomaly() -
  755.                                                               parameters.getMeanAnomaly(),
  756.                                                               0.0);


  759.             // update mean parameters
  760.             current = new BLModel(new KeplerianOrbit(current.mean.getA() + deltaA,
  761.                                                      FastMath.max(current.mean.getE() + deltaE, 0.0),
  762.                                                      current.mean.getI() + deltaI,
  763.                                                      current.mean.getPerigeeArgument() + deltaOmega,
  764.                                                      current.mean.getRightAscensionOfAscendingNode() + deltaRAAN,
  765.                                                      current.mean.getMeanAnomaly() + deltaAnom,
  766.                                                      PositionAngleType.MEAN, PositionAngleType.MEAN,
  767.                                                      current.mean.getFrame(),
  768.                                                      current.mean.getDate(), mu),
  769.                                   mass, referenceRadius, mu, ck0);
  770.             // check convergence
  771.             if (FastMath.abs(deltaA)     < thresholdA &&
  772.                 FastMath.abs(deltaE)     < thresholdE &&
  773.                 FastMath.abs(deltaI)     < thresholdAngles &&
  774.                 FastMath.abs(deltaOmega) < thresholdAngles &&
  775.                 FastMath.abs(deltaRAAN)  < thresholdAngles &&
  776.                 FastMath.abs(deltaAnom)  < thresholdAngles) {
  777.                 return current;
  778.             }
  779.         }
  780.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_BROUWER_LYDDANE_MEAN_PARAMETERS, i);
  781.     }


  784.     /** {@inheritDoc} */
  785.     public KeplerianOrbit propagateOrbit(final AbsoluteDate date) {
  786.         // compute Keplerian parameters, taking derivatives into account
  787.         final BLModel current = models.get(date);
  788.         return current.propagateParameters(date);
  789.     }

  791.     /**
  792.      * Get the value of the M2 drag parameter. Beware that M2Driver
  793.      * must have only 1 span on its TimeSpanMap value (that is
  794.      * to say setPeriod method should not be called)
  795.      * @return the value of the M2 drag parameter
  796.      */
  797.     public double getM2() {
  798.         // As Brouwer Lyddane is an analytical propagator, for now it is not possible for
  799.         // M2Driver to have several values estimated
  800.         return M2Driver.getValue();
  801.     }

  803.     /**
  804.      * Get the central attraction coefficient μ.
  805.      * @return mu central attraction coefficient (m³/s²)
  806.      */
  807.     public double getMu() {
  808.         return mu;
  809.     }

  811.     /**
  812.      * Get the un-normalized zonal coefficients.
  813.      * @return the un-normalized zonal coefficients
  814.      */
  815.     public double[] getCk0() {
  816.         return ck0.clone();
  817.     }

  819.     /**
  820.      * Get the reference radius of the central body attraction model.
  821.      * @return the reference radius in meters
  822.      */
  823.     public double getReferenceRadius() {
  824.         return referenceRadius;
  825.     }

  827.     /**
  828.      * Get the parameters driver for propagation model.
  829.      * @return drivers for propagation model
  830.      */
  831.     public List<ParameterDriver> getParametersDrivers() {
  832.         return Collections.singletonList(M2Driver);
  833.     }

  835.     /** {@inheritDoc} */
  836.     @Override
  837.     protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
  838.                                                         final DoubleArrayDictionary initialJacobianColumns) {
  839.         // Create the harvester
  840.         final BrouwerLyddaneHarvester harvester = new BrouwerLyddaneHarvester(this, stmName, initialStm, initialJacobianColumns);
  841.         // Update the list of additional state provider
  842.         addAdditionalStateProvider(harvester);
  843.         // Return the configured harvester
  844.         return harvester;
  845.     }

  847.     /**
  848.      * Get the names of the parameters in the matrix returned by {@link MatricesHarvester#getParametersJacobian}.
  849.      * @return names of the parameters (i.e. columns) of the Jacobian matrix
  850.      */
  851.     @Override
  852.     protected List<String> getJacobiansColumnsNames() {
  853.         final List<String> columnsNames = new ArrayList<>();
  854.         for (final ParameterDriver driver : getParametersDrivers()) {
  855.             if (driver.isSelected() && !columnsNames.contains(driver.getNamesSpanMap().getFirstSpan().getData())) {
  856.                 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
  857.                     columnsNames.add(span.getData());
  858.                 }
  859.             }
  860.         }
  861.         Collections.sort(columnsNames);
  862.         return columnsNames;
  863.     }

  865.     /** Local class for Brouwer-Lyddane model. */
  866.     private class BLModel {

  868.         /** Mean orbit. */
  869.         private final KeplerianOrbit mean;

  871.         /** Constant mass. */
  872.         private final double mass;

  874.         // CHECKSTYLE: stop JavadocVariable check

  876.         // preprocessed values

  878.         // Constant for secular terms l'', g'', h''
  879.         // l standing for mean anomaly, g for perigee argument and h for raan
  880.         private final double xnotDot;
  881.         private final double n;
  882.         private final double lt;
  883.         private final double gt;
  884.         private final double ht;


  887.         // Long period terms
  888.         private final double dei3sg;
  889.         private final double de2sg;
  890.         private final double deisg;
  891.         private final double de;


  894.         private final double dlgs2g;
  895.         private final double dlgc3g;
  896.         private final double dlgcg;


  899.         private final double dh2sgcg;
  900.         private final double dhsgcg;
  901.         private final double dhcg;


  904.         // Short period terms
  905.         private final double aC;
  906.         private final double aCbis;
  907.         private final double ac2g2f;


  910.         private final double eC;
  911.         private final double ecf;
  912.         private final double e2cf;
  913.         private final double e3cf;
  914.         private final double ec2f2g;
  915.         private final double ecfc2f2g;
  916.         private final double e2cfc2f2g;
  917.         private final double e3cfc2f2g;
  918.         private final double ec2gf;
  919.         private final double ec2g3f;


  922.         private final double ide;
  923.         private final double isfs2f2g;
  924.         private final double icfc2f2g;
  925.         private final double ic2f2g;


  928.         private final double glf;
  929.         private final double gll;
  930.         private final double glsf;
  931.         private final double glosf;
  932.         private final double gls2f2g;
  933.         private final double gls2gf;
  934.         private final double glos2gf;
  935.         private final double gls2g3f;
  936.         private final double glos2g3f;


  939.         private final double hf;
  940.         private final double hl;
  941.         private final double hsf;
  942.         private final double hcfs2g2f;
  943.         private final double hs2g2f;
  944.         private final double hsfc2g2f;


  947.         private final double edls2g;
  948.         private final double edlcg;
  949.         private final double edlc3g;
  950.         private final double edlsf;
  951.         private final double edls2gf;
  952.         private final double edls2g3f;

  954.         // Drag terms
  955.         private final double aRate;
  956.         private final double eRate;

  958.         // CHECKSTYLE: resume JavadocVariable check

  960.         /** Create a model for specified mean orbit.
  961.          * @param mean mean orbit
  962.          * @param mass constant mass
  963.          * @param referenceRadius reference radius of the central body attraction model (m)
  964.          * @param mu central attraction coefficient (m³/s²)
  965.          * @param ck0 un-normalized zonal coefficients
  966.          */
  967.         BLModel(final KeplerianOrbit mean, final double mass,
  968.                 final double referenceRadius, final double mu, final double[] ck0) {

  970.             this.mean = mean;
  971.             this.mass = mass;

  973.             final double app = mean.getA();
  974.             xnotDot = FastMath.sqrt(mu / app) / app;

  976.             // preliminary processing
  977.             final double q = referenceRadius / app;
  978.             double ql = q * q;
  979.             final double y2 = -0.5 * ck0[2] * ql;

  981.             n = FastMath.sqrt(1 - mean.getE() * mean.getE());
  982.             final double n2 = n * n;
  983.             final double n3 = n2 * n;
  984.             final double n4 = n2 * n2;
  985.             final double n6 = n4 * n2;
  986.             final double n8 = n4 * n4;
  987.             final double n10 = n8 * n2;

  989.             final double yp2 = y2 / n4;
  990.             ql *= q;
  991.             final double yp3 = ck0[3] * ql / n6;
  992.             ql *= q;
  993.             final double yp4 = 0.375 * ck0[4] * ql / n8;
  994.             ql *= q;
  995.             final double yp5 = ck0[5] * ql / n10;


  998.             final SinCos sc    = FastMath.sinCos(mean.getI());
  999.             final double sinI1 = sc.sin();
  1000.             final double sinI2 = sinI1 * sinI1;
  1001.             final double cosI1 = sc.cos();
  1002.             final double cosI2 = cosI1 * cosI1;
  1003.             final double cosI3 = cosI2 * cosI1;
  1004.             final double cosI4 = cosI2 * cosI2;
  1005.             final double cosI6 = cosI4 * cosI2;
  1006.             final double C5c2 = 1.0 / T2(cosI1);
  1007.             final double C3c2 = 3.0 * cosI2 - 1.0;

  1009.             final double epp = mean.getE();
  1010.             final double epp2 = epp * epp;
  1011.             final double epp3 = epp2 * epp;
  1012.             final double epp4 = epp2 * epp2;

  1014.             if (epp >= 1) {
  1015.                 // Only for elliptical (e < 1) orbits
  1016.                 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
  1017.                                           mean.getE());
  1018.             }

  1020.             // secular multiplicative
  1021.             lt = 1 +
  1022.                     1.5 * yp2 * n * C3c2 +
  1023.                     0.09375 * yp2 * yp2 * n * (-15.0 + 16.0 * n + 25.0 * n2 + (30.0 - 96.0 * n - 90.0 * n2) * cosI2 + (105.0 + 144.0 * n + 25.0 * n2) * cosI4) +
  1024.                     0.9375 * yp4 * n * epp2 * (3.0 - 30.0 * cosI2 + 35.0 * cosI4);

  1026.             gt = -1.5 * yp2 * C5c2 +
  1027.                     0.09375 * yp2 * yp2 * (-35.0 + 24.0 * n + 25.0 * n2 + (90.0 - 192.0 * n - 126.0 * n2) * cosI2 + (385.0 + 360.0 * n + 45.0 * n2) * cosI4) +
  1028.                     0.3125 * yp4 * (21.0 - 9.0 * n2 + (-270.0 + 126.0 * n2) * cosI2 + (385.0 - 189.0 * n2) * cosI4 );

  1030.             ht = -3.0 * yp2 * cosI1 +
  1031.                     0.375 * yp2 * yp2 * ((-5.0 + 12.0 * n + 9.0 * n2) * cosI1 + (-35.0 - 36.0 * n - 5.0 * n2) * cosI3) +
  1032.                     1.25 * yp4 * (5.0 - 3.0 * n2) * cosI1 * (3.0 - 7.0 * cosI2);

  1034.             final double cA = 1.0 - 11.0 * cosI2 - 40.0 * cosI4 / C5c2;
  1035.             final double cB = 1.0 - 3.0 * cosI2 - 8.0 * cosI4 / C5c2;
  1036.             final double cC = 1.0 - 9.0 * cosI2 - 24.0 * cosI4 / C5c2;
  1037.             final double cD = 1.0 - 5.0 * cosI2 - 16.0 * cosI4 / C5c2;

  1039.             final double qyp2_4 = 3.0 * yp2 * yp2 * cA -
  1040.                                   10.0 * yp4 * cB;
  1041.             final double qyp52 = epp3 * cosI1 * (0.5 * cD / sinI1 +
  1042.                                  sinI1 * (5.0 + 32.0 * cosI2 / C5c2 + 80.0 * cosI4 / C5c2 / C5c2));
  1043.             final double qyp22 = 2.0 + epp2 - 11.0 * (2.0 + 3.0 * epp2) * cosI2 -
  1044.                                  40.0 * (2.0 + 5.0 * epp2) * cosI4 / C5c2 -
  1045.                                  400.0 * epp2 * cosI6 / C5c2 / C5c2;
  1046.             final double qyp42 = ( qyp22 + 4.0 * (2.0 + epp2 - (2.0 + 3.0 * epp2) * cosI2) ) / 5.0;
  1047.             final double qyp52bis = epp * cosI1 * sinI1 *
  1048.                                     (4.0 + 3.0 * epp2) *
  1049.                                     (3.0 + 16.0 * cosI2 / C5c2 + 40.0 * cosI4 / C5c2 / C5c2);

  1051.             // long periodic multiplicative
  1052.             dei3sg =  35.0 / 96.0 * yp5 / yp2 * epp2 * n2 * cD * sinI1;
  1053.             de2sg = -1.0 / 12.0 * epp * n2 / yp2 * qyp2_4;
  1054.             deisg = ( -35.0 / 128.0 * yp5 / yp2 * epp2 * n2 * cD +
  1055.                     1.0 / 4.0 * n2 / yp2 * (yp3 + 5.0 / 16.0 * yp5 * (4.0 + 3.0 * epp2) * cC)) * sinI1;
  1056.             de = epp2 * n2 / 24.0 / yp2 * qyp2_4;

  1058.             final double qyp52quotient = epp * (-32.0 + 81.0 * epp4) / (4.0 + 3.0 * epp2 + n * (4.0 + 9.0 * epp2));
  1059.             dlgs2g = 1.0 / 48.0 / yp2 * (-3.0 * yp2 * yp2 * qyp22 +
  1060.                      10.0 * yp4 * qyp42 ) +
  1061.                      n3 / yp2 * qyp2_4 / 24.0;
  1062.             dlgc3g = 35.0 / 384.0 * yp5 / yp2 * n3 * epp * cD * sinI1 +
  1063.                      35.0 / 1152.0 * yp5 / yp2 * (2.0 * qyp52 * cosI1 - epp * cD * sinI1 * (3.0 + 2.0 * epp2));
  1064.             dlgcg = -yp3 * epp * cosI2 / ( 4.0 * yp2 * sinI1) +
  1065.                     0.078125 * yp5 / yp2 * (-epp * cosI2 / sinI1 * (4.0 + 3.0 * epp2) + epp2 * sinI1 * (26.0 + 9.0 * epp2)) * cC -
  1066.                     0.46875 * yp5 / yp2 * qyp52bis * cosI1 +
  1067.                     0.25 * yp3 / yp2 * sinI1 * epp / (1.0 + n3) * (3.0 - epp2 * (3.0 - epp2)) +
  1068.                     0.078125 * yp5 / yp2 * n2 * cC * qyp52quotient * sinI1;


  1071.             final double qyp24 = 3.0 * yp2 * yp2 * (11.0 + 80.0 * cosI2 / sinI1 + 200.0 * cosI4 / sinI2) -
  1072.                                  10.0 * yp4 * (3.0 + 16.0 * cosI2 / sinI1 + 40.0 * cosI4 / sinI2);
  1073.             dh2sgcg = 35.0 / 144.0 * yp5 / yp2 * qyp52;
  1074.             dhsgcg = -epp2 * cosI1 / (12.0 * yp2) * qyp24;
  1075.             dhcg = -35.0 / 576.0 * yp5 / yp2 * qyp52 +
  1076.                    epp * cosI1 / (4.0 * yp2 * sinI1) * (yp3 + 0.3125 * yp5 * (4.0 + 3.0 * epp2) * cC) +
  1077.                    1.875 / (4.0 * yp2) * yp5 * qyp52bis;

  1079.             // short periodic multiplicative
  1080.             aC = -yp2 * C3c2 * app / n3;
  1081.             aCbis = y2 * app * C3c2;
  1082.             ac2g2f = y2 * app * 3.0 * sinI2;

  1084.             double qe = 0.5 * n2 * y2 * C3c2 / n6;
  1085.             eC = qe * epp / (1.0 + n3) * (3.0 - epp2 * (3.0 - epp2));
  1086.             ecf = 3.0 * qe;
  1087.             e2cf = 3.0 * epp * qe;
  1088.             e3cf = epp2 * qe;
  1089.             qe = 0.5 * n2 * y2 * 3.0 * (1.0 - cosI2) / n6;
  1090.             ec2f2g = qe * epp;
  1091.             ecfc2f2g = 3.0 * qe;
  1092.             e2cfc2f2g = 3.0 * epp * qe;
  1093.             e3cfc2f2g = epp2 * qe;
  1094.             qe = -0.5 * yp2 * n2 * (1.0 - cosI2);
  1095.             ec2gf = 3.0 * qe;
  1096.             ec2g3f = qe;

  1098.             double qi = epp * yp2 * cosI1 * sinI1;
  1099.             ide = -epp * cosI1 / (n2 * sinI1);
  1100.             isfs2f2g = qi;
  1101.             icfc2f2g = 2.0 * qi;
  1102.             qi = yp2 * cosI1 * sinI1;
  1103.             ic2f2g = 1.5 * qi;

  1105.             double qgl1 = 0.25 * yp2;
  1106.             double qgl2 = 0.25 * yp2 * epp * n2 / (1.0 + n);
  1107.             glf = qgl1 * -6.0 * C5c2;
  1108.             gll = qgl1 * 6.0 * C5c2;
  1109.             glsf = qgl1 * -6.0 * C5c2 * epp +
  1110.                    qgl2 * 2.0 * C3c2;
  1111.             glosf = qgl2 * 2.0 * C3c2;
  1112.             qgl1 = qgl1 * (3.0 - 5.0 * cosI2);
  1113.             qgl2 = qgl2 * 3.0 * (1.0 - cosI2);
  1114.             gls2f2g = 3.0 * qgl1;
  1115.             gls2gf = 3.0 * epp * qgl1 +
  1116.                      qgl2;
  1117.             glos2gf = -1.0 * qgl2;
  1118.             gls2g3f = qgl1 * epp +
  1119.                       1.0 / 3.0 * qgl2;
  1120.             glos2g3f = qgl2;

  1122.             final double qh = 3.0 * yp2 * cosI1;
  1123.             hf = -qh;
  1124.             hl = qh;
  1125.             hsf = -epp * qh;
  1126.             hcfs2g2f = 2.0 * epp * yp2 * cosI1;
  1127.             hs2g2f = 1.5 * yp2 * cosI1;
  1128.             hsfc2g2f = -epp * yp2 * cosI1;

  1130.             final double qedl = -0.25 * yp2 * n3;
  1131.             edls2g = 1.0 / 24.0 * epp * n3 / yp2 * qyp2_4;
  1132.             edlcg = -0.25 * yp3 / yp2 * n3 * sinI1 -
  1133.                     0.078125 * yp5 / yp2 * n3 * sinI1 * (4.0 + 9.0 * epp2) * cC;
  1134.             edlc3g = 35.0 / 384.0 * yp5 / yp2 * n3 * epp2 * cD * sinI1;
  1135.             edlsf = 2.0 * qedl * C3c2;
  1136.             edls2gf = 3.0 * qedl * (1.0 - cosI2);
  1137.             edls2g3f = 1.0 / 3.0 * qedl;

  1139.             // secular rates of the mean semi-major axis and eccentricity
  1140.             // Eq. 2.41 and Eq. 2.45 of Phipps' 1992 thesis
  1141.             aRate = -4.0 * app / (3.0 * xnotDot);
  1142.             eRate = -4.0 * epp * n * n / (3.0 * xnotDot);

  1144.         }

  1146.         /**
  1147.          * Gets eccentric anomaly from mean anomaly.
  1148.          * @param mk the mean anomaly (rad)
  1149.          * @return the eccentric anomaly (rad)
  1150.          */
  1151.         private UnivariateDerivative1 getEccentricAnomaly(final UnivariateDerivative1 mk) {
  1152.             // reduce M to [-PI PI] interval
  1153.             final UnivariateDerivative1 reducedM = new UnivariateDerivative1(MathUtils.normalizeAngle(mk.getValue(), 0.0),
  1154.                                                                              mk.getFirstDerivative());

  1156.             final UnivariateDerivative1 meanE = new UnivariateDerivative1(mean.getE(), 0.);
  1157.             UnivariateDerivative1 ek = FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(meanE, mk);

  1159.             // expand the result back to original range
  1160.             ek = ek.add(mk.getValue() - reducedM.getValue());

  1162.             // Returns the eccentric anomaly
  1163.             return ek;
  1164.         }

  1166.         /**
  1167.          * This method is used in Brouwer-Lyddane model to avoid singularity at the
  1168.          * critical inclination (i = 63.4°).
  1169.          * <p>
  1170.          * This method, based on Warren Phipps's 1992 thesis (Eq. 2.47 and 2.48),
  1171.          * approximate the factor (1.0 - 5.0 * cos²(inc))^-1 (causing the singularity)
  1172.          * by a function, named T2 in the thesis.
  1173.          * </p>
  1174.          * @param cosInc cosine of the mean inclination
  1175.          * @return an approximation of (1.0 - 5.0 * cos²(inc))^-1 term
  1176.          */
  1177.         private double T2(final double cosInc) {

  1179.             // X = (1.0 - 5.0 * cos²(inc))
  1180.             final double x  = 1.0 - 5.0 * cosInc * cosInc;
  1181.             final double x2 = x * x;

  1183.             // Eq. 2.48
  1184.             double sum = 0.0;
  1185.             for (int i = 0; i <= 12; i++) {
  1186.                 final double sign = i % 2 == 0 ? +1.0 : -1.0;
  1187.                 sum += sign * FastMath.pow(BETA, i) * FastMath.pow(x2, i) / CombinatoricsUtils.factorialDouble(i + 1);
  1188.             }

  1190.             // Right term of equation 2.47
  1191.             double product = 1.0;
  1192.             for (int i = 0; i <= 10; i++) {
  1193.                 product *= 1 + FastMath.exp(FastMath.scalb(-1.0, i) * BETA * x2);
  1194.             }

  1196.             // Return (Eq. 2.47)
  1197.             return BETA * x * sum * product;

  1199.         }

  1201.         /** Extrapolate an orbit up to a specific target date.
  1202.          * @param date target date for the orbit
  1203.          * @return propagated parameters
  1204.          */
  1205.         public KeplerianOrbit propagateParameters(final AbsoluteDate date) {

  1207.             // Empirical drag coefficient M2
  1208.             final double m2 = M2Driver.getValue();

  1210.             // Keplerian evolution
  1211.             final UnivariateDerivative1 dt = new UnivariateDerivative1(date.durationFrom(mean.getDate()), 1.0);
  1212.             final UnivariateDerivative1 xnot = dt.multiply(xnotDot);

  1214.             //____________________________________
  1215.             // secular effects

  1217.             // mean mean anomaly (with drag Eq. 2.38 of Phipps' 1992 thesis)
  1218.             final UnivariateDerivative1 dtM2  = dt.multiply(m2);
  1219.             final UnivariateDerivative1 dt2M2 = dt.multiply(dtM2);
  1220.             final UnivariateDerivative1 lpp = new UnivariateDerivative1(MathUtils.normalizeAngle(mean.getMeanAnomaly() + lt * xnot.getValue() + dt2M2.getValue(), 0),
  1221.                                                                       lt * xnotDot + 2.0 * dtM2.getValue());
  1222.             // mean argument of perigee
  1223.             final UnivariateDerivative1 gpp = new UnivariateDerivative1(MathUtils.normalizeAngle(mean.getPerigeeArgument() + gt * xnot.getValue(), 0),
  1224.                                                                       gt * xnotDot);
  1225.             // mean longitude of ascending node
  1226.             final UnivariateDerivative1 hpp = new UnivariateDerivative1(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode() + ht * xnot.getValue(), 0),
  1227.                                                                       ht * xnotDot);

  1229.             // ________________________________________________
  1230.             // secular rates of the mean semi-major axis and eccentricity

  1232.             // semi-major axis
  1233.             final UnivariateDerivative1 appDrag = dt.multiply(aRate * m2);

  1235.             // eccentricity
  1236.             final UnivariateDerivative1 eppDrag = dt.multiply(eRate * m2);

  1238.             //____________________________________
  1239.             // Long periodical terms
  1240.             final FieldSinCos<UnivariateDerivative1> sinCosCg1 = gpp.sinCos();
  1241.             final UnivariateDerivative1 cg1 = sinCosCg1.cos();
  1242.             final UnivariateDerivative1 sg1 = sinCosCg1.sin();
  1243.             final UnivariateDerivative1 c2g = cg1.square().subtract(sg1.square());
  1244.             final UnivariateDerivative1 s2g = cg1.multiply(sg1).add(sg1.multiply(cg1));
  1245.             final UnivariateDerivative1 c3g = c2g.multiply(cg1).subtract(s2g.multiply(sg1));
  1246.             final UnivariateDerivative1 sg2 = sg1.square();
  1247.             final UnivariateDerivative1 sg3 = sg1.multiply(sg2);


  1250.             // de eccentricity
  1251.             final UnivariateDerivative1 d1e = sg3.multiply(dei3sg).
  1252.                                               add(sg1.multiply(deisg)).
  1253.                                               add(sg2.multiply(de2sg)).
  1254.                                               add(de);

  1256.             // l' + g'
  1257.             final UnivariateDerivative1 lpPGp = s2g.multiply(dlgs2g).
  1258.                                                add(c3g.multiply(dlgc3g)).
  1259.                                                add(cg1.multiply(dlgcg)).
  1260.                                                add(lpp).
  1261.                                                add(gpp);

  1263.             // h'
  1264.             final UnivariateDerivative1 hp = sg2.multiply(cg1).multiply(dh2sgcg).
  1265.                                                add(sg1.multiply(cg1).multiply(dhsgcg)).
  1266.                                                add(cg1.multiply(dhcg)).
  1267.                                                add(hpp);

  1269.             // Short periodical terms
  1270.             // eccentric anomaly
  1271.             final UnivariateDerivative1 Ep = getEccentricAnomaly(lpp);
  1272.             final FieldSinCos<UnivariateDerivative1> sinCosEp = Ep.sinCos();
  1273.             final UnivariateDerivative1 cf1 = (sinCosEp.cos().subtract(mean.getE())).divide(sinCosEp.cos().multiply(-mean.getE()).add(1.0));
  1274.             final UnivariateDerivative1 sf1 = (sinCosEp.sin().multiply(n)).divide(sinCosEp.cos().multiply(-mean.getE()).add(1.0));
  1275.             final UnivariateDerivative1 f = FastMath.atan2(sf1, cf1);

  1277.             final UnivariateDerivative1 cf2 = cf1.square();
  1278.             final UnivariateDerivative1 c2f = cf2.subtract(sf1.multiply(sf1));
  1279.             final UnivariateDerivative1 s2f = cf1.multiply(sf1).add(sf1.multiply(cf1));
  1280.             final UnivariateDerivative1 c3f = c2f.multiply(cf1).subtract(s2f.multiply(sf1));
  1281.             final UnivariateDerivative1 s3f = c2f.multiply(sf1).add(s2f.multiply(cf1));
  1282.             final UnivariateDerivative1 cf3 = cf1.multiply(cf2);

  1284.             final UnivariateDerivative1 c2g1f = cf1.multiply(c2g).subtract(sf1.multiply(s2g));
  1285.             final UnivariateDerivative1 c2g2f = c2f.multiply(c2g).subtract(s2f.multiply(s2g));
  1286.             final UnivariateDerivative1 c2g3f = c3f.multiply(c2g).subtract(s3f.multiply(s2g));
  1287.             final UnivariateDerivative1 s2g1f = cf1.multiply(s2g).add(c2g.multiply(sf1));
  1288.             final UnivariateDerivative1 s2g2f = c2f.multiply(s2g).add(c2g.multiply(s2f));
  1289.             final UnivariateDerivative1 s2g3f = c3f.multiply(s2g).add(c2g.multiply(s3f));

  1291.             final UnivariateDerivative1 eE = (sinCosEp.cos().multiply(-mean.getE()).add(1.0)).reciprocal();
  1292.             final UnivariateDerivative1 eE3 = eE.square().multiply(eE);
  1293.             final UnivariateDerivative1 sigma = eE.multiply(n * n).multiply(eE).add(eE);

  1295.             // Semi-major axis (with drag Eq. 2.41 of Phipps' 1992 thesis)
  1296.             final UnivariateDerivative1 a = eE3.multiply(aCbis).add(appDrag.add(mean.getA())).
  1297.                                             add(aC).
  1298.                                             add(eE3.multiply(c2g2f).multiply(ac2g2f));

  1300.             // Eccentricity (with drag Eq. 2.45 of Phipps' 1992 thesis)
  1301.             final UnivariateDerivative1 e = d1e.add(eppDrag.add(mean.getE())).
  1302.                                             add(eC).
  1303.                                             add(cf1.multiply(ecf)).
  1304.                                             add(cf2.multiply(e2cf)).
  1305.                                             add(cf3.multiply(e3cf)).
  1306.                                             add(c2g2f.multiply(ec2f2g)).
  1307.                                             add(c2g2f.multiply(cf1).multiply(ecfc2f2g)).
  1308.                                             add(c2g2f.multiply(cf2).multiply(e2cfc2f2g)).
  1309.                                             add(c2g2f.multiply(cf3).multiply(e3cfc2f2g)).
  1310.                                             add(c2g1f.multiply(ec2gf)).
  1311.                                             add(c2g3f.multiply(ec2g3f));

  1313.             // Inclination
  1314.             final UnivariateDerivative1 i = d1e.multiply(ide).
  1315.                                             add(mean.getI()).
  1316.                                             add(sf1.multiply(s2g2f.multiply(isfs2f2g))).
  1317.                                             add(cf1.multiply(c2g2f.multiply(icfc2f2g))).
  1318.                                             add(c2g2f.multiply(ic2f2g));

  1320.             // Argument of perigee + mean anomaly
  1321.             final UnivariateDerivative1 gPL = lpPGp.add(f.multiply(glf)).
  1322.                                                 add(lpp.multiply(gll)).
  1323.                                                 add(sf1.multiply(glsf)).
  1324.                                                 add(sigma.multiply(sf1).multiply(glosf)).
  1325.                                                 add(s2g2f.multiply(gls2f2g)).
  1326.                                                 add(s2g1f.multiply(gls2gf)).
  1327.                                                 add(sigma.multiply(s2g1f).multiply(glos2gf)).
  1328.                                                 add(s2g3f.multiply(gls2g3f)).
  1329.                                                 add(sigma.multiply(s2g3f).multiply(glos2g3f));


  1332.             // Longitude of ascending node
  1333.             final UnivariateDerivative1 h = hp.add(f.multiply(hf)).
  1334.                                             add(lpp.multiply(hl)).
  1335.                                             add(sf1.multiply(hsf)).
  1336.                                             add(cf1.multiply(s2g2f).multiply(hcfs2g2f)).
  1337.                                             add(s2g2f.multiply(hs2g2f)).
  1338.                                             add(c2g2f.multiply(sf1).multiply(hsfc2g2f));

  1340.             final UnivariateDerivative1 edl = s2g.multiply(edls2g).
  1341.                                             add(cg1.multiply(edlcg)).
  1342.                                             add(c3g.multiply(edlc3g)).
  1343.                                             add(sf1.multiply(edlsf)).
  1344.                                             add(s2g1f.multiply(edls2gf)).
  1345.                                             add(s2g3f.multiply(edls2g3f)).
  1346.                                             add(sf1.multiply(sigma).multiply(edlsf)).
  1347.                                             add(s2g1f.multiply(sigma).multiply(-edls2gf)).
  1348.                                             add(s2g3f.multiply(sigma).multiply(3.0 * edls2g3f));

  1350.             final FieldSinCos<UnivariateDerivative1> sinCosLpp = lpp.sinCos();
  1351.             final UnivariateDerivative1 A = e.multiply(sinCosLpp.cos()).subtract(edl.multiply(sinCosLpp.sin()));
  1352.             final UnivariateDerivative1 B = e.multiply(sinCosLpp.sin()).add(edl.multiply(sinCosLpp.cos()));

  1354.             // Mean anomaly
  1355.             final UnivariateDerivative1 l = FastMath.atan2(B, A);

  1357.             // Argument of perigee
  1358.             final UnivariateDerivative1 g = gPL.subtract(l);

  1360.             // Return a Keplerian orbit
  1361.             return new KeplerianOrbit(a.getValue(), e.getValue(), i.getValue(),
  1362.                                       g.getValue(), h.getValue(), l.getValue(),
  1363.                                       a.getFirstDerivative(), e.getFirstDerivative(), i.getFirstDerivative(),
  1364.                                       g.getFirstDerivative(), h.getFirstDerivative(), l.getFirstDerivative(),
  1365.                                       PositionAngleType.MEAN, mean.getFrame(), date, mu);

  1367.         }

  1369.     }

  1371.     /** {@inheritDoc} */
  1372.     protected double getMass(final AbsoluteDate date) {
  1373.         return models.get(date).mass;
  1374.     }

  1376. }
Created with JaCoCo 0.8.12.202403310830