FieldBrouwerLyddanePropagator.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.Collections;
  20. import java.util.List;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.Field;
  24. import org.hipparchus.analysis.differentiation.FieldUnivariateDerivative1;
  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.orekit.attitudes.AttitudeProvider;
  30. import org.orekit.attitudes.FrameAlignedProvider;
  31. import org.orekit.errors.OrekitException;
  32. import org.orekit.errors.OrekitMessages;
  33. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  34. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  35. import org.orekit.orbits.FieldKeplerianOrbit;
  36. import org.orekit.orbits.FieldOrbit;
  37. import org.orekit.orbits.FieldKeplerianAnomalyUtility;
  38. import org.orekit.orbits.OrbitType;
  39. import org.orekit.orbits.PositionAngleType;
  40. import org.orekit.propagation.FieldSpacecraftState;
  41. import org.orekit.propagation.PropagationType;
  42. import org.orekit.propagation.analytical.tle.FieldTLE;
  43. import org.orekit.time.AbsoluteDate;
  44. import org.orekit.time.FieldAbsoluteDate;
  45. import org.orekit.utils.FieldTimeSpanMap;
  46. import org.orekit.utils.ParameterDriver;

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

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

  93.     /** Default value for maxIterations. */
  94.     private static final int MAX_ITERATIONS_DEFAULT = 200;

  96.     /** Parameters scaling factor.
  97.      * <p>
  98.      * We use a power of 2 to avoid numeric noise introduction
  99.      * in the multiplications/divisions sequences.
  100.      * </p>
  101.      */
  102.     private static final double SCALE = FastMath.scalb(1.0, -20);

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

  107.     /** Initial Brouwer-Lyddane model. */
  108.     private FieldBLModel<T> initialModel;

  110.     /** All models. */
  111.     private transient FieldTimeSpanMap<FieldBLModel<T>, T> models;

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

  116.     /** Central attraction coefficient (m³/s²). */
  117.     private T mu;

  119.     /** Un-normalized zonal coefficients. */
  120.     private double[] ck0;

  122.     /** Empirical coefficient used in the drag modeling. */
  123.     private final ParameterDriver M2Driver;

  125.     /** Build a propagator from orbit and potential provider.
  126.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  127.      *
  128.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  129.      *
  130.      * @param initialOrbit initial orbit
  131.      * @param provider for un-normalized zonal coefficients
  132.      * @param M2 value of empirical drag coefficient in rad/s².
  133.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  134.      * @see #FieldBrouwerLyddanePropagator(FieldOrbit, UnnormalizedSphericalHarmonicsProvider, PropagationType, double)
  135.      */
  136.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  137.                                          final UnnormalizedSphericalHarmonicsProvider provider,
  138.                                          final double M2) {
  139.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  140.              initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
  141.              provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
  142.     }

  144.     /**
  145.      * Private helper constructor.
  146.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  147.      * @param initialOrbit initial orbit
  148.      * @param attitude attitude provider
  149.      * @param mass spacecraft mass
  150.      * @param provider for un-normalized zonal coefficients
  151.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  152.      * @param M2 value of empirical drag coefficient in rad/s².
  153.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  154.      * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement,
  155.      * UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics, PropagationType, double)
  156.      */
  157.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  158.                                          final AttitudeProvider attitude,
  159.                                          final T mass,
  160.                                          final UnnormalizedSphericalHarmonicsProvider provider,
  161.                                          final UnnormalizedSphericalHarmonics harmonics,
  162.                                          final double M2) {
  163.         this(initialOrbit, attitude,  mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
  164.              harmonics.getUnnormalizedCnm(2, 0),
  165.              harmonics.getUnnormalizedCnm(3, 0),
  166.              harmonics.getUnnormalizedCnm(4, 0),
  167.              harmonics.getUnnormalizedCnm(5, 0),
  168.              M2);
  169.     }

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

  206.     /** Build a propagator from orbit, mass and potential provider.
  207.      * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
  208.      *
  209.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  210.      *
  211.      * @param initialOrbit initial orbit
  212.      * @param mass spacecraft mass
  213.      * @param provider for un-normalized zonal coefficients
  214.      * @param M2 value of empirical drag coefficient in rad/s².
  215.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  216.      * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider, double)
  217.      */
  218.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit, final T mass,
  219.                                          final UnnormalizedSphericalHarmonicsProvider provider,
  220.                                          final double M2) {
  221.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  222.              mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
  223.     }

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

  259.     /** Build a propagator from orbit, attitude provider and potential provider.
  260.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  261.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  262.      * @param initialOrbit initial orbit
  263.      * @param attitudeProv attitude provider
  264.      * @param provider for un-normalized zonal coefficients
  265.      * @param M2 value of empirical drag coefficient in rad/s².
  266.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  267.      */
  268.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  269.                                          final AttitudeProvider attitudeProv,
  270.                                          final UnnormalizedSphericalHarmonicsProvider provider,
  271.                                          final double M2) {
  272.         this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
  273.              provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
  274.     }

  276.     /** Build a propagator from orbit, attitude provider and potential.
  277.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  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.      *
  283.      * <p> 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.      *
  286.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  287.      *
  288.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  289.      *
  290.      * @param initialOrbit initial orbit
  291.      * @param attitudeProv attitude provider
  292.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  293.      * @param mu central attraction coefficient (m³/s²)
  294.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  295.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  296.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  297.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  298.      * @param M2 value of empirical drag coefficient in rad/s².
  299.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  300.      */
  301.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  302.                                          final AttitudeProvider attitudeProv,
  303.                                          final double referenceRadius, final T mu,
  304.                                          final double c20, final double c30, final double c40,
  305.                                          final double c50, final double M2) {
  306.         this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS),
  307.              referenceRadius, mu, c20, c30, c40, c50, M2);
  308.     }

  310.     /** Build a propagator from orbit, attitude provider, mass and potential provider.
  311.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  312.      * @param initialOrbit initial orbit
  313.      * @param attitudeProv attitude provider
  314.      * @param mass spacecraft mass
  315.      * @param provider for un-normalized zonal coefficients
  316.      * @param M2 value of empirical drag coefficient in rad/s².
  317.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  318.      * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider, PropagationType, double)
  319.      */
  320.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  321.                                          final AttitudeProvider attitudeProv,
  322.                                          final T mass,
  323.                                          final UnnormalizedSphericalHarmonicsProvider provider,
  324.                                          final double M2) {
  325.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
  326.     }

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


  364.     /** Build a propagator from orbit and potential provider.
  365.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  366.      *
  367.      * <p>Using this constructor, it is possible to define the initial orbit as
  368.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  369.      *
  370.      * @param initialOrbit initial orbit
  371.      * @param provider for un-normalized zonal coefficients
  372.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  373.      * @param M2 value of empirical drag coefficient in rad/s².
  374.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  375.      */
  376.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  377.                                          final UnnormalizedSphericalHarmonicsProvider provider,
  378.                                          final PropagationType initialType,
  379.                                          final double M2) {
  380.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  381.              initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
  382.              provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType, M2);
  383.     }

  385.     /** Build a propagator from orbit, attitude provider, mass and potential provider.
  386.      * <p>Using this constructor, it is possible to define the initial orbit as
  387.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  388.      * @param initialOrbit initial orbit
  389.      * @param attitudeProv attitude provider
  390.      * @param mass spacecraft mass
  391.      * @param provider for un-normalized zonal coefficients
  392.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  393.      * @param M2 value of empirical drag coefficient in rad/s².
  394.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  395.      */
  396.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  397.                                          final AttitudeProvider attitudeProv,
  398.                                          final T mass,
  399.                                          final UnnormalizedSphericalHarmonicsProvider provider,
  400.                                          final PropagationType initialType,
  401.                                          final double M2) {
  402.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType, M2);
  403.     }

  405.     /**
  406.      * Private helper constructor.
  407.      * <p>Using this constructor, it is possible to define the initial orbit as
  408.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  409.      * @param initialOrbit initial orbit
  410.      * @param attitude attitude provider
  411.      * @param mass spacecraft mass
  412.      * @param provider for un-normalized zonal coefficients
  413.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  414.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  415.      * @param M2 value of empirical drag coefficient in rad/s².
  416.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  417.      */
  418.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  419.                                          final AttitudeProvider attitude,
  420.                                          final T mass,
  421.                                          final UnnormalizedSphericalHarmonicsProvider provider,
  422.                                          final UnnormalizedSphericalHarmonics harmonics,
  423.                                          final PropagationType initialType,
  424.                                          final double M2) {
  425.         this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
  426.              harmonics.getUnnormalizedCnm(2, 0),
  427.              harmonics.getUnnormalizedCnm(3, 0),
  428.              harmonics.getUnnormalizedCnm(4, 0),
  429.              harmonics.getUnnormalizedCnm(5, 0),
  430.              initialType, M2);
  431.     }

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

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

  511.         super(mass.getField(), attitudeProv);

  513.         // store model coefficients
  514.         this.referenceRadius = referenceRadius;
  515.         this.mu  = mu;
  516.         this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50};

  518.         // initialize M2 driver
  519.         this.M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, M2, SCALE,
  520.                                             Double.NEGATIVE_INFINITY,
  521.                                             Double.POSITIVE_INFINITY);

  523.         // compute mean parameters if needed
  524.         resetInitialState(new FieldSpacecraftState<>(initialOrbit,
  525.                                                  attitudeProv.getAttitude(initialOrbit,
  526.                                                                           initialOrbit.getDate(),
  527.                                                                           initialOrbit.getFrame()),
  528.                                                  mass),
  529.                                                  initialType, epsilon, maxIterations);

  531.     }

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

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

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

  647.     /** {@inheritDoc}
  648.      * <p>The new initial state to consider
  649.      * must be defined with an osculating orbit.</p>
  650.      * @see #resetInitialState(FieldSpacecraftState, PropagationType)
  651.      */
  652.     @Override
  653.     public void resetInitialState(final FieldSpacecraftState<T> state) {
  654.         resetInitialState(state, PropagationType.OSCULATING);
  655.     }

  657.     /** Reset the propagator initial state.
  658.      * @param state new initial state to consider
  659.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  660.      */
  661.     public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType) {
  662.         resetInitialState(state, stateType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  663.     }

  665.     /** Reset the propagator initial state.
  666.      * @param state new initial state to consider
  667.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  668.      * @param epsilon convergence threshold for mean parameters conversion
  669.      * @param maxIterations maximum iterations for mean parameters conversion
  670.      * @since 11.2
  671.      */
  672.     public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType,
  673.                                   final double epsilon, final int maxIterations) {
  674.         super.resetInitialState(state);
  675.         final FieldKeplerianOrbit<T> keplerian = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(state.getOrbit());
  676.         this.initialModel = (stateType == PropagationType.MEAN) ?
  677.                              new FieldBLModel<>(keplerian, state.getMass(), referenceRadius, mu, ck0) :
  678.                              computeMeanParameters(keplerian, state.getMass(), epsilon, maxIterations);
  679.         this.models = new FieldTimeSpanMap<>(initialModel, state.getA().getField());
  680.     }

  682.     /** {@inheritDoc} */
  683.     @Override
  684.     protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward) {
  685.         resetIntermediateState(state, forward, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  686.     }

  688.     /** Reset an intermediate state.
  689.      * @param state new intermediate state to consider
  690.      * @param forward if true, the intermediate state is valid for
  691.      * propagations after itself
  692.      * @param epsilon convergence threshold for mean parameters conversion
  693.      * @param maxIterations maximum iterations for mean parameters conversion
  694.      * @since 11.2
  695.      */
  696.     protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward,
  697.                                           final double epsilon, final int maxIterations) {
  698.         final FieldBLModel<T> newModel = computeMeanParameters((FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(state.getOrbit()),
  699.                                                                state.getMass(), epsilon, maxIterations);
  700.         if (forward) {
  701.             models.addValidAfter(newModel, state.getDate());
  702.         } else {
  703.             models.addValidBefore(newModel, state.getDate());
  704.         }
  705.         stateChanged(state);
  706.     }

  708.     /** Compute mean parameters according to the Brouwer-Lyddane analytical model computation
  709.      * in order to do the propagation.
  710.      * @param osculating osculating orbit
  711.      * @param mass constant mass
  712.      * @param epsilon convergence threshold for mean parameters conversion
  713.      * @param maxIterations maximum iterations for mean parameters conversion
  714.      * @return Brouwer-Lyddane mean model
  715.      */
  716.     private FieldBLModel<T> computeMeanParameters(final FieldKeplerianOrbit<T> osculating, final T mass,
  717.                                                   final double epsilon, final int maxIterations) {

  719.         // sanity check
  720.         if (osculating.getA().getReal() < referenceRadius) {
  721.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
  722.                                            osculating.getA());
  723.         }

  725.         final Field<T> field = mass.getField();
  726.         final T one = field.getOne();
  727.         final T zero = field.getZero();
  728.         // rough initialization of the mean parameters
  729.         FieldBLModel<T> current = new FieldBLModel<>(osculating, mass, referenceRadius, mu, ck0);

  731.         // threshold for each parameter
  732.         final T thresholdA      = current.mean.getA().abs().add(1.0).multiply(epsilon);
  733.         final T thresholdE      = current.mean.getE().add(1.0).multiply(epsilon);
  734.         final T thresholdAngles = one.getPi().multiply(epsilon);

  736.         int i = 0;
  737.         while (i++ < maxIterations) {

  739.             // recompute the osculating parameters from the current mean parameters
  740.             final FieldKeplerianOrbit<T> parameters = current.propagateParameters(current.mean.getDate(), getParameters(mass.getField(), current.mean.getDate()));

  742.             // adapted parameters residuals
  743.             final T deltaA     = osculating.getA()  .subtract(parameters.getA());
  744.             final T deltaE     = osculating.getE()  .subtract(parameters.getE());
  745.             final T deltaI     = osculating.getI()  .subtract(parameters.getI());
  746.             final T deltaOmega = MathUtils.normalizeAngle(osculating.getPerigeeArgument().subtract(parameters.getPerigeeArgument()), zero);
  747.             final T deltaRAAN  = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode().subtract(parameters.getRightAscensionOfAscendingNode()), zero);
  748.             final T deltaAnom  = MathUtils.normalizeAngle(osculating.getMeanAnomaly().subtract(parameters.getMeanAnomaly()), zero);


  751.             // update mean parameters
  752.             current = new FieldBLModel<>(new FieldKeplerianOrbit<>(current.mean.getA()            .add(deltaA),
  753.                                                      FastMath.max(current.mean.getE().add(deltaE), zero),
  754.                                                      current.mean.getI()                            .add(deltaI),
  755.                                                      current.mean.getPerigeeArgument()              .add(deltaOmega),
  756.                                                      current.mean.getRightAscensionOfAscendingNode().add(deltaRAAN),
  757.                                                      current.mean.getMeanAnomaly()                  .add(deltaAnom),
  758.                                                      PositionAngleType.MEAN, PositionAngleType.MEAN,
  759.                                                      current.mean.getFrame(),
  760.                                                      current.mean.getDate(), mu),
  761.                                   mass, referenceRadius, mu, ck0);
  762.             // check convergence
  763.             if (FastMath.abs(deltaA.getReal())     < thresholdA.getReal() &&
  764.                 FastMath.abs(deltaE.getReal())     < thresholdE.getReal() &&
  765.                 FastMath.abs(deltaI.getReal())     < thresholdAngles.getReal() &&
  766.                 FastMath.abs(deltaOmega.getReal()) < thresholdAngles.getReal() &&
  767.                 FastMath.abs(deltaRAAN.getReal())  < thresholdAngles.getReal() &&
  768.                 FastMath.abs(deltaAnom.getReal())  < thresholdAngles.getReal()) {
  769.                 return current;
  770.             }
  771.         }
  772.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_BROUWER_LYDDANE_MEAN_PARAMETERS, i);
  773.     }


  776.     /** {@inheritDoc} */
  777.     public FieldKeplerianOrbit<T> propagateOrbit(final FieldAbsoluteDate<T> date, final T[] parameters) {
  778.         // compute Cartesian parameters, taking derivatives into account
  779.         final FieldBLModel<T> current = models.get(date);
  780.         return current.propagateParameters(date, parameters);
  781.     }

  783.     /**
  784.      * Get the value of the M2 drag parameter.
  785.      * @return the value of the M2 drag parameter
  786.      */
  787.     public double getM2() {
  788.         return M2Driver.getValue();
  789.     }

  791.     /**
  792.      * Get the value of the M2 drag parameter.
  793.      * @param date date at which the model parameters want to be known
  794.      * @return the value of the M2 drag parameter
  795.      */
  796.     public double getM2(final AbsoluteDate date) {
  797.         return M2Driver.getValue(date);
  798.     }

  800.     /** Local class for Brouwer-Lyddane model. */
  801.     private static class FieldBLModel<T extends CalculusFieldElement<T>> {

  803.         /** Mean orbit. */
  804.         private final FieldKeplerianOrbit<T> mean;

  806.         /** Constant mass. */
  807.         private final T mass;

  809.         /** Central attraction coefficient. */
  810.         private final T mu;

  812.         // CHECKSTYLE: stop JavadocVariable check

  814.         // preprocessed values

  816.         // Constant for secular terms l'', g'', h''
  817.         // l standing for mean anomaly, g for perigee argument and h for raan
  818.         private final T xnotDot;
  819.         private final T n;
  820.         private final T lt;
  821.         private final T gt;
  822.         private final T ht;


  825.         // Long periodT
  826.         private final T dei3sg;
  827.         private final T de2sg;
  828.         private final T deisg;
  829.         private final T de;


  832.         private final T dlgs2g;
  833.         private final T dlgc3g;
  834.         private final T dlgcg;


  837.         private final T dh2sgcg;
  838.         private final T dhsgcg;
  839.         private final T dhcg;


  842.         // Short perioTs
  843.         private final T aC;
  844.         private final T aCbis;
  845.         private final T ac2g2f;


  848.         private final T eC;
  849.         private final T ecf;
  850.         private final T e2cf;
  851.         private final T e3cf;
  852.         private final T ec2f2g;
  853.         private final T ecfc2f2g;
  854.         private final T e2cfc2f2g;
  855.         private final T e3cfc2f2g;
  856.         private final T ec2gf;
  857.         private final T ec2g3f;


  860.         private final T ide;
  861.         private final T isfs2f2g;
  862.         private final T icfc2f2g;
  863.         private final T ic2f2g;


  866.         private final T glf;
  867.         private final T gll;
  868.         private final T glsf;
  869.         private final T glosf;
  870.         private final T gls2f2g;
  871.         private final T gls2gf;
  872.         private final T glos2gf;
  873.         private final T gls2g3f;
  874.         private final T glos2g3f;


  877.         private final T hf;
  878.         private final T hl;
  879.         private final T hsf;
  880.         private final T hcfs2g2f;
  881.         private final T hs2g2f;
  882.         private final T hsfc2g2f;


  885.         private final T edls2g;
  886.         private final T edlcg;
  887.         private final T edlc3g;
  888.         private final T edlsf;
  889.         private final T edls2gf;
  890.         private final T edls2g3f;

  892.         // Drag terms
  893.         private final T aRate;
  894.         private final T eRate;

  896.         // CHECKSTYLE: resume JavadocVariable check

  898.         /** Create a model for specified mean orbit.
  899.          * @param mean mean Fieldorbit
  900.          * @param mass constant mass
  901.          * @param referenceRadius reference radius of the central body attraction model (m)
  902.          * @param mu central attraction coefficient (m³/s²)
  903.          * @param ck0 un-normalized zonal coefficients
  904.          */
  905.         FieldBLModel(final FieldKeplerianOrbit<T> mean, final T mass,
  906.                 final double referenceRadius, final T mu, final double[] ck0) {

  908.             this.mean = mean;
  909.             this.mass = mass;
  910.             this.mu   = mu;
  911.             final T one  = mass.getField().getOne();

  913.             final T app = mean.getA();
  914.             xnotDot = mu.divide(app).sqrt().divide(app);

  916.             // preliminary processing
  917.             final T q = app.divide(referenceRadius).reciprocal();
  918.             T ql = q.square();
  919.             final T y2 = ql.multiply(-0.5 * ck0[2]);

  921.             n = ((mean.getE().square().negate()).add(1.0)).sqrt();
  922.             final T n2 = n.square();
  923.             final T n3 = n2.multiply(n);
  924.             final T n4 = n2.square();
  925.             final T n6 = n4.multiply(n2);
  926.             final T n8 = n4.square();
  927.             final T n10 = n8.multiply(n2);

  929.             final T yp2 = y2.divide(n4);
  930.             ql = ql.multiply(q);
  931.             final T yp3 = ql.multiply(ck0[3]).divide(n6);
  932.             ql = ql.multiply(q);
  933.             final T yp4 = ql.multiply(0.375 * ck0[4]).divide(n8);
  934.             ql = ql.multiply(q);
  935.             final T yp5 = ql.multiply(ck0[5]).divide(n10);


  938.             final FieldSinCos<T> sc = FastMath.sinCos(mean.getI());
  939.             final T sinI1 = sc.sin();
  940.             final T sinI2 = sinI1.square();
  941.             final T cosI1 = sc.cos();
  942.             final T cosI2 = cosI1.square();
  943.             final T cosI3 = cosI2.multiply(cosI1);
  944.             final T cosI4 = cosI2.square();
  945.             final T cosI6 = cosI4.multiply(cosI2);
  946.             final T C5c2 = T2(cosI1).reciprocal();
  947.             final T C3c2 = cosI2.multiply(3.0).subtract(1.0);

  949.             final T epp = mean.getE();
  950.             final T epp2 = epp.square();
  951.             final T epp3 = epp2.multiply(epp);
  952.             final T epp4 = epp2.square();

  954.             if (epp.getReal() >= 1) {
  955.                 // Only for elliptical (e < 1) orbits
  956.                 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
  957.                                           mean.getE().getReal());
  958.             }

  960.             // secular multiplicative
  961.             lt = one.add(yp2.multiply(n).multiply(C3c2).multiply(1.5)).
  962.                      add(yp2.multiply(0.09375).multiply(yp2).multiply(n).multiply(n2.multiply(25.0).add(n.multiply(16.0)).add(-15.0).
  963.                              add(cosI2.multiply(n2.multiply(-90.0).add(n.multiply(-96.0)).add(30.0))).
  964.                              add(cosI4.multiply(n2.multiply(25.0).add(n.multiply(144.0)).add(105.0))))).
  965.                      add(yp4.multiply(0.9375).multiply(n).multiply(epp2).multiply(cosI4.multiply(35.0).add(cosI2.multiply(-30.0)).add(3.0)));

  967.             gt = yp2.multiply(-1.5).multiply(C5c2).
  968.                      add(yp2.multiply(0.09375).multiply(yp2).multiply(n2.multiply(25.0).add(n.multiply(24.0)).add(-35.0).
  969.                             add(cosI2.multiply(n2.multiply(-126.0).add(n.multiply(-192.0)).add(90.0))).
  970.                             add(cosI4.multiply(n2.multiply(45.0).add(n.multiply(360.0)).add(385.0))))).
  971.                      add(yp4.multiply(0.3125).multiply(n2.multiply(-9.0).add(21.0).
  972.                             add(cosI2.multiply(n2.multiply(126.0).add(-270.0))).
  973.                             add(cosI4.multiply(n2.multiply(-189.0).add(385.0)))));

  975.             ht = yp2.multiply(-3.0).multiply(cosI1).
  976.                      add(yp2.multiply(0.375).multiply(yp2).multiply(cosI1.multiply(n2.multiply(9.0).add(n.multiply(12.0)).add(-5.0)).
  977.                                                                     add(cosI3.multiply(n2.multiply(-5.0).add(n.multiply(-36.0)).add(-35.0))))).
  978.                      add(yp4.multiply(1.25).multiply(cosI1).multiply(n2.multiply(-3.0).add(5.0)).multiply(cosI2.multiply(-7.0).add(3.0)));

  980.             final T cA = one.subtract(cosI2.multiply(11.0)).subtract(cosI4.multiply(40.0).divide(C5c2));
  981.             final T cB = one.subtract(cosI2.multiply(3.0)).subtract(cosI4.multiply(8.0).divide(C5c2));
  982.             final T cC = one.subtract(cosI2.multiply(9.0)).subtract(cosI4.multiply(24.0).divide(C5c2));
  983.             final T cD = one.subtract(cosI2.multiply(5.0)).subtract(cosI4.multiply(16.0).divide(C5c2));

  985.             final T qyp2_4 = yp2.multiply(3.0).multiply(yp2).multiply(cA).
  986.                              subtract(yp4.multiply(10.0).multiply(cB));
  987.             final T qyp52 = cosI1.multiply(epp3).multiply(cD.multiply(0.5).divide(sinI1).
  988.                                                           add(sinI1.multiply(cosI4.divide(C5c2).divide(C5c2).multiply(80.0).
  989.                                                                              add(cosI2.divide(C5c2).multiply(32.0).
  990.                                                                              add(5.0)))));
  991.             final T qyp22 = epp2.add(2.0).subtract(epp2.multiply(3.0).add(2.0).multiply(11.0).multiply(cosI2)).
  992.                             subtract(epp2.multiply(5.0).add(2.0).multiply(40.0).multiply(cosI4.divide(C5c2))).
  993.                             subtract(epp2.multiply(400.0).multiply(cosI6).divide(C5c2).divide(C5c2));
  994.             final T qyp42 = one.divide(5.0).multiply(qyp22.
  995.                                                      add(one.newInstance(4.0).multiply(epp2.
  996.                                                                                     add(2.0).
  997.                                                                                     subtract(cosI2.multiply(epp2.multiply(3.0).add(2.0))))));
  998.             final T qyp52bis = cosI1.multiply(sinI1).multiply(epp).multiply(epp2.multiply(3.0).add(4.0)).
  999.                                                                    multiply(cosI4.divide(C5c2).divide(C5c2).multiply(40.0).
  1000.                                                                             add(cosI2.divide(C5c2).multiply(16.0)).
  1001.                                                                             add(3.0));
  1002.            // long periodic multiplicative
  1003.             dei3sg =  yp5.divide(yp2).multiply(35.0 / 96.0).multiply(epp2).multiply(n2).multiply(cD).multiply(sinI1);
  1004.             de2sg = qyp2_4.divide(yp2).multiply(epp).multiply(n2).multiply(-1.0 / 12.0);
  1005.             deisg = sinI1.multiply(yp5.multiply(-35.0 / 128.0).divide(yp2).multiply(epp2).multiply(n2).multiply(cD).
  1006.                             add(n2.multiply(0.25).divide(yp2).multiply(yp3.add(yp5.multiply(5.0 / 16.0).multiply(epp2.multiply(3.0).add(4.0)).multiply(cC)))));
  1007.             de = epp2.multiply(n2).multiply(qyp2_4).divide(24.0).divide(yp2);

  1009.             final T qyp52quotient = epp.multiply(epp4.multiply(81.0).add(-32.0)).divide(n.multiply(epp2.multiply(9.0).add(4.0)).add(epp2.multiply(3.0)).add(4.0));
  1010.             dlgs2g = yp2.multiply(48.0).reciprocal().multiply(yp2.multiply(-3.0).multiply(yp2).multiply(qyp22).
  1011.                             add(yp4.multiply(10.0).multiply(qyp42))).
  1012.                             add(n3.divide(yp2).multiply(qyp2_4).divide(24.0));
  1013.             dlgc3g =  yp5.multiply(35.0 / 384.0).divide(yp2).multiply(n3).multiply(epp).multiply(cD).multiply(sinI1).
  1014.                             add(yp5.multiply(35.0 / 1152.0).divide(yp2).multiply(qyp52.multiply(2.0).multiply(cosI1).subtract(epp.multiply(cD).multiply(sinI1).multiply(epp2.multiply(2.0).add(3.0)))));
  1015.             dlgcg = yp3.negate().multiply(cosI2).multiply(epp).divide(yp2.multiply(sinI1).multiply(4.0)).
  1016.                     add(yp5.divide(yp2).multiply(0.078125).multiply(cC).multiply(cosI2.divide(sinI1).multiply(epp.negate()).multiply(epp2.multiply(3.0).add(4.0)).
  1017.                                                                     add(sinI1.multiply(epp2).multiply(epp2.multiply(9.0).add(26.0)))).
  1018.                     subtract(yp5.divide(yp2).multiply(0.46875).multiply(qyp52bis).multiply(cosI1)).
  1019.                     add(yp3.divide(yp2).multiply(0.25).multiply(sinI1).multiply(epp).divide(n3.add(1.0)).multiply(epp2.multiply(epp2.subtract(3.0)).add(3.0))).
  1020.                     add(yp5.divide(yp2).multiply(0.078125).multiply(n2).multiply(cC).multiply(qyp52quotient).multiply(sinI1)));
  1021.             final T qyp24 = yp2.multiply(yp2).multiply(3.0).multiply(cosI4.divide(sinI2).multiply(200.0).add(cosI2.divide(sinI1).multiply(80.0)).add(11.0)).
  1022.                             subtract(yp4.multiply(10.0).multiply(cosI4.divide(sinI2).multiply(40.0).add(cosI2.divide(sinI1).multiply(16.0)).add(3.0)));
  1023.             dh2sgcg = yp5.divide(yp2).multiply(qyp52).multiply(35.0 / 144.0);
  1024.             dhsgcg = qyp24.multiply(cosI1).multiply(epp2.negate()).divide(yp2.multiply(12.0));
  1025.             dhcg = yp5.divide(yp2).multiply(qyp52).multiply(-35.0 / 576.0).
  1026.                    add(cosI1.multiply(epp).divide(yp2.multiply(sinI1).multiply(4.0)).multiply(yp3.add(yp5.multiply(0.3125).multiply(cC).multiply(epp2.multiply(3.0).add(4.0))))).
  1027.                    add(yp5.multiply(qyp52bis).multiply(1.875).divide(yp2.multiply(4.0)));

  1029.             // short periodic multiplicative
  1030.             aC = yp2.negate().multiply(C3c2).multiply(app).divide(n3);
  1031.             aCbis = y2.multiply(app).multiply(C3c2);
  1032.             ac2g2f = y2.multiply(app).multiply(sinI2).multiply(3.0);

  1034.             T qe = y2.multiply(C3c2).multiply(0.5).multiply(n2).divide(n6);
  1035.             eC = qe.multiply(epp).divide(n3.add(1.0)).multiply(epp2.multiply(epp2.subtract(3.0)).add(3.0));
  1036.             ecf = qe.multiply(3.0);
  1037.             e2cf = qe.multiply(3.0).multiply(epp);
  1038.             e3cf = qe.multiply(epp2);
  1039.             qe = y2.multiply(0.5).multiply(n2).multiply(3.0).multiply(cosI2.negate().add(1.0)).divide(n6);
  1040.             ec2f2g = qe.multiply(epp);
  1041.             ecfc2f2g = qe.multiply(3.0);
  1042.             e2cfc2f2g = qe.multiply(3.0).multiply(epp);
  1043.             e3cfc2f2g = qe.multiply(epp2);
  1044.             qe = yp2.multiply(-0.5).multiply(n2).multiply(cosI2.negate().add(1.0));
  1045.             ec2gf = qe.multiply(3.0);
  1046.             ec2g3f = qe;

  1048.             T qi = yp2.multiply(epp).multiply(cosI1).multiply(sinI1);
  1049.             ide = cosI1.multiply(epp.negate()).divide(sinI1.multiply(n2));
  1050.             isfs2f2g = qi;
  1051.             icfc2f2g = qi.multiply(2.0);
  1052.             qi = yp2.multiply(cosI1).multiply(sinI1);
  1053.             ic2f2g = qi.multiply(1.5);

  1055.             T qgl1 = yp2.multiply(0.25);
  1056.             T qgl2 =  yp2.multiply(epp).multiply(n2).multiply(0.25).divide(n.add(1.0));
  1057.             glf = qgl1.multiply(C5c2).multiply(-6.0);
  1058.             gll = qgl1.multiply(C5c2).multiply(6.0);
  1059.             glsf = qgl1.multiply(C5c2).multiply(-6.0).multiply(epp).
  1060.                    add(qgl2.multiply(C3c2).multiply(2.0));
  1061.             glosf = qgl2.multiply(C3c2).multiply(2.0);
  1062.             qgl1 = qgl1.multiply(cosI2.multiply(-5.0).add(3.0));
  1063.             qgl2 = qgl2.multiply(3.0).multiply(cosI2.negate().add(1.0));
  1064.             gls2f2g = qgl1.multiply(3.0);
  1065.             gls2gf = qgl1.multiply(3.0).multiply(epp).
  1066.                      add(qgl2);
  1067.             glos2gf = qgl2.negate();
  1068.             gls2g3f = qgl1.multiply(epp).
  1069.                       add(qgl2.multiply(1.0 / 3.0));
  1070.             glos2g3f = qgl2;

  1072.             final T qh = yp2.multiply(cosI1).multiply(3.0);
  1073.             hf = qh.negate();
  1074.             hl = qh;
  1075.             hsf = qh.multiply(epp).negate();
  1076.             hcfs2g2f = yp2.multiply(cosI1).multiply(epp).multiply(2.0);
  1077.             hs2g2f = yp2.multiply(cosI1).multiply(1.5);
  1078.             hsfc2g2f = yp2.multiply(cosI1).multiply(epp).negate();

  1080.             final T qedl = yp2.multiply(n3).multiply(-0.25);
  1081.             edls2g = qyp2_4.multiply(1.0 / 24.0).multiply(epp).multiply(n3).divide(yp2);
  1082.             edlcg = yp3.divide(yp2).multiply(-0.25).multiply(n3).multiply(sinI1).
  1083.                     subtract(yp5.divide(yp2).multiply(0.078125).multiply(n3).multiply(sinI1).multiply(cC).multiply(epp2.multiply(9.0).add(4.0)));
  1084.             edlc3g = yp5.divide(yp2).multiply(n3).multiply(epp2).multiply(cD).multiply(sinI1).multiply(35.0 / 384.0);
  1085.             edlsf = qedl.multiply(C3c2).multiply(2.0);
  1086.             edls2gf = qedl.multiply(3.0).multiply(cosI2.negate().add(1.0));
  1087.             edls2g3f = qedl.multiply(1.0 / 3.0);

  1089.             // secular rates of the mean semi-major axis and eccentricity
  1090.             // Eq. 2.41 and Eq. 2.45 of Phipps' 1992 thesis
  1091.             aRate = app.multiply(-4.0).divide(xnotDot.multiply(3.0));
  1092.             eRate = epp.multiply(n).multiply(n).multiply(-4.0).divide(xnotDot.multiply(3.0));

  1094.         }

  1096.         /**
  1097.          * Gets eccentric anomaly from mean anomaly.
  1098.          * @param mk the mean anomaly (rad)
  1099.          * @return the eccentric anomaly (rad)
  1100.          */
  1101.         private FieldUnivariateDerivative1<T> getEccentricAnomaly(final FieldUnivariateDerivative1<T> mk) {

  1103.             final T zero = mean.getE().getField().getZero();

  1105.             // reduce M to [-PI PI] interval
  1106.             final FieldUnivariateDerivative1<T> reducedM = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(mk.getValue(), zero ),
  1107.                                                                              mk.getFirstDerivative());

  1109.             final FieldUnivariateDerivative1<T> meanE = new FieldUnivariateDerivative1<>(mean.getE(), zero);
  1110.             FieldUnivariateDerivative1<T> ek = FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(meanE, mk);

  1112.             // expand the result back to original range
  1113.             ek = ek.add(mk.getValue().subtract(reducedM.getValue()));

  1115.             // Returns the eccentric anomaly
  1116.             return ek;
  1117.         }

  1119.         /**
  1120.          * This method is used in Brouwer-Lyddane model to avoid singularity at the
  1121.          * critical inclination (i = 63.4°).
  1122.          * <p>
  1123.          * This method, based on Warren Phipps's 1992 thesis (Eq. 2.47 and 2.48),
  1124.          * approximate the factor (1.0 - 5.0 * cos²(inc))^-1 (causing the singularity)
  1125.          * by a function, named T2 in the thesis.
  1126.          * </p>
  1127.          * @param cosInc cosine of the mean inclination
  1128.          * @return an approximation of (1.0 - 5.0 * cos²(inc))^-1 term
  1129.          */
  1130.         private T T2(final T cosInc) {

  1132.             // X = (1.0 - 5.0 * cos²(inc))
  1133.             final T x  = cosInc.square().multiply(-5.0).add(1.0);
  1134.             final T x2 = x.square();

  1136.             // Eq. 2.48
  1137.             T sum = x.getField().getZero();
  1138.             for (int i = 0; i <= 12; i++) {
  1139.                 final double sign = i % 2 == 0 ? +1.0 : -1.0;
  1140.                 sum = sum.add(FastMath.pow(x2, i).multiply(FastMath.pow(BETA, i)).multiply(sign).divide(CombinatoricsUtils.factorialDouble(i + 1)));
  1141.             }

  1143.             // Right term of equation 2.47
  1144.             T product = x.getField().getOne();
  1145.             for (int i = 0; i <= 10; i++) {
  1146.                 product = product.multiply(FastMath.exp(x2.multiply(BETA).multiply(FastMath.scalb(-1.0, i))).add(1.0));
  1147.             }

  1149.             // Return (Eq. 2.47)
  1150.             return x.multiply(BETA).multiply(sum).multiply(product);

  1152.         }

  1154.         /** Extrapolate an orbit up to a specific target date.
  1155.          * @param date target date for the orbit
  1156.          * @param parameters model parameters
  1157.          * @return propagated parameters
  1158.          */
  1159.         public FieldKeplerianOrbit<T> propagateParameters(final FieldAbsoluteDate<T> date, final T[] parameters) {

  1161.             // Field
  1162.             final Field<T> field = date.getField();
  1163.             final T one  = field.getOne();
  1164.             final T zero = field.getZero();

  1166.             // Empirical drag coefficient M2
  1167.             final T m2 = parameters[0];

  1169.             // Keplerian evolution
  1170.             final FieldUnivariateDerivative1<T> dt = new FieldUnivariateDerivative1<>(date.durationFrom(mean.getDate()), one);
  1171.             final FieldUnivariateDerivative1<T> xnot = dt.multiply(xnotDot);

  1173.             //____________________________________
  1174.             // secular effects

  1176.             // mean mean anomaly
  1177.             final FieldUnivariateDerivative1<T> dtM2  = dt.multiply(m2);
  1178.             final FieldUnivariateDerivative1<T> dt2M2 = dt.multiply(dtM2);
  1179.             final FieldUnivariateDerivative1<T> lpp = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(mean.getMeanAnomaly().add(lt.multiply(xnot.getValue())).add(dt2M2.getValue()), zero),
  1180.                                                                                         lt.multiply(xnotDot).add(dtM2.multiply(2.0).getValue()));
  1181.             // mean argument of perigee
  1182.             final FieldUnivariateDerivative1<T> gpp = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(mean.getPerigeeArgument().add(gt.multiply(xnot.getValue())), zero),
  1183.                                                                                         gt.multiply(xnotDot));
  1184.             // mean longitude of ascending node
  1185.             final FieldUnivariateDerivative1<T> hpp = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode().add(ht.multiply(xnot.getValue())), zero),
  1186.                                                                                         ht.multiply(xnotDot));

  1188.             // ________________________________________________
  1189.             // secular rates of the mean semi-major axis and eccentricity

  1191.             // semi-major axis
  1192.             final FieldUnivariateDerivative1<T> appDrag = dt.multiply(aRate.multiply(m2));

  1194.             // eccentricity
  1195.             final FieldUnivariateDerivative1<T> eppDrag = dt.multiply(eRate.multiply(m2));

  1197.             //____________________________________
  1198.             // Long periodical terms
  1199.             final FieldSinCos<FieldUnivariateDerivative1<T>> sinCosGpp = gpp.sinCos();
  1200.             final FieldUnivariateDerivative1<T> cg1 = sinCosGpp.cos();
  1201.             final FieldUnivariateDerivative1<T> sg1 = sinCosGpp.sin();
  1202.             final FieldUnivariateDerivative1<T> c2g = cg1.multiply(cg1).subtract(sg1.multiply(sg1));
  1203.             final FieldUnivariateDerivative1<T> s2g = cg1.multiply(sg1).add(sg1.multiply(cg1));
  1204.             final FieldUnivariateDerivative1<T> c3g = c2g.multiply(cg1).subtract(s2g.multiply(sg1));
  1205.             final FieldUnivariateDerivative1<T> sg2 = sg1.square();
  1206.             final FieldUnivariateDerivative1<T> sg3 = sg1.multiply(sg2);


  1209.             // de eccentricity
  1210.             final FieldUnivariateDerivative1<T> d1e = sg3.multiply(dei3sg).
  1211.                                                add(sg1.multiply(deisg)).
  1212.                                                add(sg2.multiply(de2sg)).
  1213.                                                add(de);

  1215.             // l' + g'
  1216.             final FieldUnivariateDerivative1<T> lpPGp = s2g.multiply(dlgs2g).
  1217.                                                add(c3g.multiply(dlgc3g)).
  1218.                                                add(cg1.multiply(dlgcg)).
  1219.                                                add(lpp).
  1220.                                                add(gpp);

  1222.             // h'
  1223.             final FieldUnivariateDerivative1<T> hp = sg2.multiply(cg1).multiply(dh2sgcg).
  1224.                                                add(sg1.multiply(cg1).multiply(dhsgcg)).
  1225.                                                add(cg1.multiply(dhcg)).
  1226.                                                add(hpp);

  1228.             // Short periodical terms
  1229.             // eccentric anomaly
  1230.             final FieldUnivariateDerivative1<T> Ep = getEccentricAnomaly(lpp);
  1231.             final FieldSinCos<FieldUnivariateDerivative1<T>> sinCosEp = Ep.sinCos();
  1232.             final FieldUnivariateDerivative1<T> cf1 = (sinCosEp.cos().subtract(mean.getE())).divide(sinCosEp.cos().multiply(mean.getE().negate()).add(1.0));
  1233.             final FieldUnivariateDerivative1<T> sf1 = (sinCosEp.sin().multiply(n)).divide(sinCosEp.cos().multiply(mean.getE().negate()).add(1.0));
  1234.             final FieldUnivariateDerivative1<T> f = FastMath.atan2(sf1, cf1);

  1236.             final FieldUnivariateDerivative1<T> cf2 = cf1.square();
  1237.             final FieldUnivariateDerivative1<T> c2f = cf2.subtract(sf1.multiply(sf1));
  1238.             final FieldUnivariateDerivative1<T> s2f = cf1.multiply(sf1).add(sf1.multiply(cf1));
  1239.             final FieldUnivariateDerivative1<T> c3f = c2f.multiply(cf1).subtract(s2f.multiply(sf1));
  1240.             final FieldUnivariateDerivative1<T> s3f = c2f.multiply(sf1).add(s2f.multiply(cf1));
  1241.             final FieldUnivariateDerivative1<T> cf3 = cf1.multiply(cf2);

  1243.             final FieldUnivariateDerivative1<T> c2g1f = cf1.multiply(c2g).subtract(sf1.multiply(s2g));
  1244.             final FieldUnivariateDerivative1<T> c2g2f = c2f.multiply(c2g).subtract(s2f.multiply(s2g));
  1245.             final FieldUnivariateDerivative1<T> c2g3f = c3f.multiply(c2g).subtract(s3f.multiply(s2g));
  1246.             final FieldUnivariateDerivative1<T> s2g1f = cf1.multiply(s2g).add(c2g.multiply(sf1));
  1247.             final FieldUnivariateDerivative1<T> s2g2f = c2f.multiply(s2g).add(c2g.multiply(s2f));
  1248.             final FieldUnivariateDerivative1<T> s2g3f = c3f.multiply(s2g).add(c2g.multiply(s3f));

  1250.             final FieldUnivariateDerivative1<T> eE = (sinCosEp.cos().multiply(mean.getE().negate()).add(1.0)).reciprocal();
  1251.             final FieldUnivariateDerivative1<T> eE3 = eE.square().multiply(eE);
  1252.             final FieldUnivariateDerivative1<T> sigma = eE.multiply(n.square()).multiply(eE).add(eE);

  1254.             // Semi-major axis
  1255.             final FieldUnivariateDerivative1<T> a = eE3.multiply(aCbis).add(appDrag.add(mean.getA())).
  1256.                                             add(aC).
  1257.                                             add(eE3.multiply(c2g2f).multiply(ac2g2f));

  1259.             // Eccentricity
  1260.             final FieldUnivariateDerivative1<T> e = d1e.add(eppDrag.add(mean.getE())).
  1261.                                             add(eC).
  1262.                                             add(cf1.multiply(ecf)).
  1263.                                             add(cf2.multiply(e2cf)).
  1264.                                             add(cf3.multiply(e3cf)).
  1265.                                             add(c2g2f.multiply(ec2f2g)).
  1266.                                             add(c2g2f.multiply(cf1).multiply(ecfc2f2g)).
  1267.                                             add(c2g2f.multiply(cf2).multiply(e2cfc2f2g)).
  1268.                                             add(c2g2f.multiply(cf3).multiply(e3cfc2f2g)).
  1269.                                             add(c2g1f.multiply(ec2gf)).
  1270.                                             add(c2g3f.multiply(ec2g3f));

  1272.             // Inclination
  1273.             final FieldUnivariateDerivative1<T> i = d1e.multiply(ide).
  1274.                                             add(mean.getI()).
  1275.                                             add(sf1.multiply(s2g2f.multiply(isfs2f2g))).
  1276.                                             add(cf1.multiply(c2g2f.multiply(icfc2f2g))).
  1277.                                             add(c2g2f.multiply(ic2f2g));

  1279.             // Argument of perigee + mean anomaly
  1280.             final FieldUnivariateDerivative1<T> gPL = lpPGp.add(f.multiply(glf)).
  1281.                                              add(lpp.multiply(gll)).
  1282.                                              add(sf1.multiply(glsf)).
  1283.                                              add(sigma.multiply(sf1).multiply(glosf)).
  1284.                                              add(s2g2f.multiply(gls2f2g)).
  1285.                                              add(s2g1f.multiply(gls2gf)).
  1286.                                              add(sigma.multiply(s2g1f).multiply(glos2gf)).
  1287.                                              add(s2g3f.multiply(gls2g3f)).
  1288.                                              add(sigma.multiply(s2g3f).multiply(glos2g3f));

  1290.             // Longitude of ascending node
  1291.             final FieldUnivariateDerivative1<T> h = hp.add(f.multiply(hf)).
  1292.                                             add(lpp.multiply(hl)).
  1293.                                             add(sf1.multiply(hsf)).
  1294.                                             add(cf1.multiply(s2g2f).multiply(hcfs2g2f)).
  1295.                                             add(s2g2f.multiply(hs2g2f)).
  1296.                                             add(c2g2f.multiply(sf1).multiply(hsfc2g2f));

  1298.             final FieldUnivariateDerivative1<T> edl = s2g.multiply(edls2g).
  1299.                                             add(cg1.multiply(edlcg)).
  1300.                                             add(c3g.multiply(edlc3g)).
  1301.                                             add(sf1.multiply(edlsf)).
  1302.                                             add(s2g1f.multiply(edls2gf)).
  1303.                                             add(s2g3f.multiply(edls2g3f)).
  1304.                                             add(sf1.multiply(sigma).multiply(edlsf)).
  1305.                                             add(s2g1f.multiply(sigma).multiply(edls2gf.negate())).
  1306.                                             add(s2g3f.multiply(sigma).multiply(edls2g3f.multiply(3.0)));

  1308.             final FieldUnivariateDerivative1<T> A = e.multiply(lpp.cos()).subtract(edl.multiply(lpp.sin()));
  1309.             final FieldUnivariateDerivative1<T> B = e.multiply(lpp.sin()).add(edl.multiply(lpp.cos()));

  1311.             // Mean anomaly
  1312.             final FieldUnivariateDerivative1<T> l = FastMath.atan2(B, A);

  1314.             // Argument of perigee
  1315.             final FieldUnivariateDerivative1<T> g = gPL.subtract(l);

  1317.             // Return a Keplerian orbit
  1318.             return new FieldKeplerianOrbit<>(a.getValue(), e.getValue(), i.getValue(),
  1319.                                              g.getValue(), h.getValue(), l.getValue(),
  1320.                                              a.getFirstDerivative(), e.getFirstDerivative(), i.getFirstDerivative(),
  1321.                                              g.getFirstDerivative(), h.getFirstDerivative(), l.getFirstDerivative(),
  1322.                                              PositionAngleType.MEAN, mean.getFrame(), date, this.mu);

  1324.         }
  1325.     }

  1327.     /** {@inheritDoc} */
  1328.     @Override
  1329.     protected T getMass(final FieldAbsoluteDate<T> date) {
  1330.         return models.get(date).mass;
  1331.     }

  1333.     /** {@inheritDoc} */
  1334.     @Override
  1335.     public List<ParameterDriver> getParametersDrivers() {
  1336.         return Collections.singletonList(M2Driver);
  1337.     }

  1339. }
Created with JaCoCo 0.8.12.202403310830