FieldBrouwerLyddanePropagator.java

  1. /* Copyright 2002-2025 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.FieldKeplerianAnomalyUtility;
  36. import org.orekit.orbits.FieldKeplerianOrbit;
  37. import org.orekit.orbits.FieldOrbit;
  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.propagation.conversion.osc2mean.BrouwerLyddaneTheory;
  44. import org.orekit.propagation.conversion.osc2mean.FixedPointConverter;
  45. import org.orekit.propagation.conversion.osc2mean.MeanTheory;
  46. import org.orekit.propagation.conversion.osc2mean.OsculatingToMeanConverter;
  47. import org.orekit.time.AbsoluteDate;
  48. import org.orekit.time.FieldAbsoluteDate;
  49. import org.orekit.utils.FieldTimeSpanMap;
  50. import org.orekit.utils.ParameterDriver;

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

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

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

  109.     /** Max value for the eccentricity. */
  110.     private static final double MAX_ECC = 0.999999;

  112.     /** Initial Brouwer-Lyddane model. */
  113.     private FieldBLModel<T> initialModel;

  115.     /** All models. */
  116.     private transient FieldTimeSpanMap<FieldBLModel<T>, T> models;

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

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

  124.     /** Un-normalized zonal coefficients. */
  125.     private double[] ck0;

  127.     /** Empirical coefficient used in the drag modeling. */
  128.     private final ParameterDriver M2Driver;

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

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

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

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

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

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

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

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

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


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

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

  416.     /**
  417.      * Private helper constructor.
  418.      * <p>Using this constructor, it is possible to define the initial orbit as
  419.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  420.      * @param initialOrbit initial orbit
  421.      * @param attitude attitude provider
  422.      * @param mass spacecraft mass
  423.      * @param provider for un-normalized zonal coefficients
  424.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  425.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  426.      * @param m2Value value of empirical drag coefficient in rad/s².
  427.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  428.      */
  429.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  430.                                          final AttitudeProvider attitude,
  431.                                          final T mass,
  432.                                          final UnnormalizedSphericalHarmonicsProvider provider,
  433.                                          final UnnormalizedSphericalHarmonics harmonics,
  434.                                          final PropagationType initialType,
  435.                                          final double m2Value) {
  436.         this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
  437.              harmonics.getUnnormalizedCnm(2, 0),
  438.              harmonics.getUnnormalizedCnm(3, 0),
  439.              harmonics.getUnnormalizedCnm(4, 0),
  440.              harmonics.getUnnormalizedCnm(5, 0),
  441.              initialType, m2Value);
  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 m2Value value of empirical drag coefficient in rad/s².
  469.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  470.      */
  471.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  472.                                          final AttitudeProvider attitudeProv,
  473.                                          final T mass,
  474.                                          final double referenceRadius, final T mu,
  475.                                          final double c20, final double c30, final double c40,
  476.                                          final double c50,
  477.                                          final PropagationType initialType,
  478.                                          final double m2Value) {
  479.         this(initialOrbit, attitudeProv, mass, referenceRadius, mu,
  480.              c20, c30, c40, c50, initialType, m2Value,
  481.              new FixedPointConverter(BrouwerLyddanePropagator.EPSILON_DEFAULT,
  482.                                      BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT,
  483.                                      FixedPointConverter.DEFAULT_DAMPING));
  484.     }

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

  534.     /** Build a propagator from orbit, attitude provider, mass and potential.
  535.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  536.      * are related to both the normalized coefficients
  537.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  538.      *  and the J<sub>n</sub> one as follows:</p>
  539.      *
  540.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  541.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  542.      *
  543.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  544.      *
  545.      * <p>Using this constructor, it is possible to define the initial orbit as
  546.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  547.      *
  548.      * @param initialOrbit initial orbit
  549.      * @param attitudeProv attitude provider
  550.      * @param mass spacecraft mass
  551.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  552.      * @param mu central attraction coefficient (m³/s²)
  553.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  554.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  555.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  556.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  557.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  558.      * @param m2Value value of empirical drag coefficient in rad/s².
  559.      *        If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
  560.      * @param converter osculating to mean orbit converter
  561.      * @since 13.0
  562.      */
  563.     public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
  564.                                          final AttitudeProvider attitudeProv,
  565.                                          final T mass,
  566.                                          final double referenceRadius,
  567.                                          final T mu,
  568.                                          final double c20,
  569.                                          final double c30,
  570.                                          final double c40,
  571.                                          final double c50,
  572.                                          final PropagationType initialType,
  573.                                          final double m2Value,
  574.                                          final OsculatingToMeanConverter converter) {

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

  578.         // store model coefficients
  579.         this.referenceRadius = referenceRadius;
  580.         this.mu  = mu;
  581.         this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50};

  583.         // initialize M2 driver
  584.         this.M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, m2Value, SCALE,
  585.                                             Double.NEGATIVE_INFINITY,
  586.                                             Double.POSITIVE_INFINITY);

  588.         // compute mean parameters if needed
  589.         resetInitialState(new FieldSpacecraftState<>(initialOrbit,
  590.                                                      attitudeProv.getAttitude(initialOrbit,
  591.                                                                               initialOrbit.getDate(),
  592.                                                                               initialOrbit.getFrame())).withMass(mass),
  593.                           initialType, converter);

  595.     }

  597.     /** Conversion from osculating to mean orbit.
  598.      * <p>
  599.      * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
  600.      * osculating SpacecraftState in input.
  601.      * </p>
  602.      * <p>
  603.      * Since the osculating orbit is obtained with the computation of
  604.      * short-periodic variation, the resulting output will depend on
  605.      * both the gravity field parameterized in input and the
  606.      * atmospheric drag represented by the {@code m2} parameter.
  607.      * </p>
  608.      * <p>
  609.      * The computation is done through a fixed-point iteration process.
  610.      * </p>
  611.      * @param <T> type of the filed elements
  612.      * @param osculating osculating orbit to convert
  613.      * @param provider for un-normalized zonal coefficients
  614.      * @param harmonics {@code provider.onDate(osculating.getDate())}
  615.      * @param m2Value value of empirical drag coefficient in rad/s².
  616.      *        If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
  617.      * @return mean orbit in a Brouwer-Lyddane sense
  618.      * @since 11.2
  619.      */
  620.     public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
  621.                                                                                               final UnnormalizedSphericalHarmonicsProvider provider,
  622.                                                                                               final UnnormalizedSphericalHarmonics harmonics,
  623.                                                                                               final double m2Value) {
  624.         return computeMeanOrbit(osculating, provider, harmonics, m2Value,
  625.                                 BrouwerLyddanePropagator.EPSILON_DEFAULT,
  626.                                 BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT);
  627.     }

  629.     /** Conversion from osculating to mean orbit.
  630.      * <p>
  631.      * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
  632.      * osculating SpacecraftState in input.
  633.      * </p>
  634.      * <p>
  635.      * Since the osculating orbit is obtained with the computation of
  636.      * short-periodic variation, the resulting output will depend on
  637.      * both the gravity field parameterized in input and the
  638.      * atmospheric drag represented by the {@code m2} parameter.
  639.      * </p>
  640.      * <p>
  641.      * The computation is done through a fixed-point iteration process.
  642.      * </p>
  643.      * @param <T> type of the filed elements
  644.      * @param osculating osculating orbit to convert
  645.      * @param provider for un-normalized zonal coefficients
  646.      * @param harmonics {@code provider.onDate(osculating.getDate())}
  647.      * @param m2Value value of empirical drag coefficient in rad/s².
  648.      *        If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
  649.      * @param epsilon convergence threshold for mean parameters conversion
  650.      * @param maxIterations maximum iterations for mean parameters conversion
  651.      * @return mean orbit in a Brouwer-Lyddane sense
  652.      * @since 11.2
  653.      */
  654.     public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
  655.                                                                                               final UnnormalizedSphericalHarmonicsProvider provider,
  656.                                                                                               final UnnormalizedSphericalHarmonics harmonics,
  657.                                                                                               final double m2Value,
  658.                                                                                               final double epsilon,
  659.                                                                                               final int maxIterations) {
  660.         return computeMeanOrbit(osculating,
  661.                                 provider.getAe(), provider.getMu(),
  662.                                 harmonics.getUnnormalizedCnm(2, 0),
  663.                                 harmonics.getUnnormalizedCnm(3, 0),
  664.                                 harmonics.getUnnormalizedCnm(4, 0),
  665.                                 harmonics.getUnnormalizedCnm(5, 0),
  666.                                 m2Value, epsilon, maxIterations);
  667.     }

  669.     /** Conversion from osculating to mean orbit.
  670.      * <p>
  671.      * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
  672.      * osculating SpacecraftState in input.
  673.      * </p>
  674.      * <p>
  675.      * Since the osculating orbit is obtained with the computation of
  676.      * short-periodic variation, the resulting output will depend on
  677.      * both the gravity field parameterized in input and the
  678.      * atmospheric drag represented by the {@code m2} parameter.
  679.      * </p>
  680.      * <p>
  681.      * The computation is done through a fixed-point iteration process.
  682.      * </p>
  683.      * @param <T> type of the filed elements
  684.      * @param osculating osculating orbit to convert
  685.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  686.      * @param mu central attraction coefficient (m³/s²)
  687.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  688.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  689.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  690.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  691.      * @param m2Value value of empirical drag coefficient in rad/s².
  692.      *        If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
  693.      * @param epsilon convergence threshold for mean parameters conversion
  694.      * @param maxIterations maximum iterations for mean parameters conversion
  695.      * @return mean orbit in a Brouwer-Lyddane sense
  696.      * @since 11.2
  697.      */
  698.     public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
  699.                                                                                               final double referenceRadius,
  700.                                                                                               final double mu,
  701.                                                                                               final double c20,
  702.                                                                                               final double c30,
  703.                                                                                               final double c40,
  704.                                                                                               final double c50,
  705.                                                                                               final double m2Value,
  706.                                                                                               final double epsilon,
  707.                                                                                               final int maxIterations) {
  708.         // Build a fixed-point converter
  709.         final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
  710.                                                                             FixedPointConverter.DEFAULT_DAMPING);
  711.         return computeMeanOrbit(osculating, referenceRadius, mu, c20, c30, c40, c50, m2Value, converter);
  712.     }

  714.     /** Conversion from osculating to mean orbit.
  715.      * <p>
  716.      * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
  717.      * osculating SpacecraftState in input.
  718.      * </p>
  719.      * <p>
  720.      * Since the osculating orbit is obtained with the computation of
  721.      * short-periodic variation, the resulting output will depend on
  722.      * both the gravity field parameterized in input and the
  723.      * atmospheric drag represented by the {@code m2} parameter.
  724.      * </p>
  725.      * <p>
  726.      * The computation is done through the given osculating to mean orbit converter.
  727.      * </p>
  728.      * @param <T> type of the filed elements
  729.      * @param osculating osculating orbit to convert
  730.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  731.      * @param mu central attraction coefficient (m³/s²)
  732.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  733.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  734.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  735.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  736.      * @param m2Value value of empirical drag coefficient in rad/s².
  737.      *        If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
  738.      * @param converter osculating to mean orbit converter
  739.      * @return mean orbit in a Brouwer-Lyddane sense
  740.      * @since 13.0
  741.      */
  742.     public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
  743.                                                                                               final double referenceRadius,
  744.                                                                                               final double mu,
  745.                                                                                               final double c20,
  746.                                                                                               final double c30,
  747.                                                                                               final double c40,
  748.                                                                                               final double c50,
  749.                                                                                               final double m2Value,
  750.                                                                                               final OsculatingToMeanConverter converter) {
  751.         // Set BL as the mean theory for converting
  752.         final MeanTheory theory = new BrouwerLyddaneTheory(referenceRadius, mu, c20, c30, c40, c50, m2Value);
  753.         converter.setMeanTheory(theory);
  754.         return (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(converter.convertToMean(osculating));
  755.     }

  757.     /** {@inheritDoc}
  758.      * <p>The new initial state to consider
  759.      * must be defined with an osculating orbit.</p>
  760.      * @see #resetInitialState(FieldSpacecraftState, PropagationType)
  761.      */
  762.     @Override
  763.     public void resetInitialState(final FieldSpacecraftState<T> state) {
  764.         resetInitialState(state, PropagationType.OSCULATING);
  765.     }

  767.     /** Reset the propagator initial state.
  768.      * @param state new initial state to consider
  769.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  770.      */
  771.     public void resetInitialState(final FieldSpacecraftState<T> state,
  772.                                   final PropagationType stateType) {
  773.         final OsculatingToMeanConverter converter = new FixedPointConverter(BrouwerLyddanePropagator.EPSILON_DEFAULT,
  774.                                                                             BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT,
  775.                                                                             FixedPointConverter.DEFAULT_DAMPING);
  776.         resetInitialState(state, stateType, converter);
  777.     }

  779.     /** Reset the propagator initial state.
  780.      * @param state new initial state to consider
  781.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  782.      * @param epsilon convergence threshold for mean parameters conversion
  783.      * @param maxIterations maximum iterations for mean parameters conversion
  784.      * @since 11.2
  785.      */
  786.     public void resetInitialState(final FieldSpacecraftState<T> state,
  787.                                   final PropagationType stateType,
  788.                                   final double epsilon,
  789.                                   final int maxIterations) {
  790.         final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
  791.                                                                             FixedPointConverter.DEFAULT_DAMPING);
  792.         resetInitialState(state, stateType, converter);
  793.     }

  795.     /** Reset the propagator initial state.
  796.      * @param state     new initial state to consider
  797.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  798.      * @param converter osculating to mean orbit converter
  799.      * @since 13.0
  800.      */
  801.     public void resetInitialState(final FieldSpacecraftState<T> state,
  802.                                   final PropagationType stateType,
  803.                                   final OsculatingToMeanConverter converter) {
  804.         super.resetInitialState(state);
  805.         FieldKeplerianOrbit<T> keplerian = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(state.getOrbit());
  806.         if (stateType == PropagationType.OSCULATING) {
  807.             final MeanTheory theory = new BrouwerLyddaneTheory(referenceRadius, mu.getReal(),
  808.                                                                ck0[2], ck0[3], ck0[4], ck0[5],
  809.                                                                getM2());
  810.             converter.setMeanTheory(theory);
  811.             keplerian = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(converter.convertToMean(keplerian));
  812.         }
  813.         this.initialModel = new FieldBLModel<>(keplerian, state.getMass(), referenceRadius, mu, ck0);
  814.         this.models = new FieldTimeSpanMap<>(initialModel, state.getMass().getField());
  815.     }

  817.     /** {@inheritDoc} */
  818.     @Override
  819.     protected void resetIntermediateState(final FieldSpacecraftState<T> state,
  820.                                           final boolean forward) {
  821.         final OsculatingToMeanConverter converter = new FixedPointConverter(BrouwerLyddanePropagator.EPSILON_DEFAULT,
  822.                                                                             BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT,
  823.                                                                             FixedPointConverter.DEFAULT_DAMPING);
  824.         resetIntermediateState(state, forward, converter);
  825.     }

  827.     /** Reset an intermediate state.
  828.      * @param state new intermediate state to consider
  829.      * @param forward if true, the intermediate state is valid for
  830.      *                propagations after itself
  831.      * @param epsilon convergence threshold for mean parameters conversion
  832.      * @param maxIterations maximum iterations for mean parameters conversion
  833.      * @since 11.2
  834.      */
  835.     protected void resetIntermediateState(final FieldSpacecraftState<T> state,
  836.                                           final boolean forward,
  837.                                           final double epsilon,
  838.                                           final int maxIterations) {
  839.         final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
  840.                                                                             FixedPointConverter.DEFAULT_DAMPING);
  841.         resetIntermediateState(state, forward, converter);
  842.     }

  844.     /** Reset an intermediate state.
  845.      * @param state     new intermediate state to consider
  846.      * @param forward   if true, the intermediate state is valid for
  847.      *                  propagations after itself
  848.      * @param converter osculating to mean orbit converter
  849.      * @since 13.0
  850.      */
  851.     protected void resetIntermediateState(final FieldSpacecraftState<T> state,
  852.                                           final boolean forward,
  853.                                           final OsculatingToMeanConverter converter) {
  854.         final MeanTheory theory = new BrouwerLyddaneTheory(referenceRadius, mu.getReal(),
  855.                                                            ck0[2], ck0[3], ck0[4], ck0[5],
  856.                                                            getM2());
  857.         converter.setMeanTheory(theory);
  858.         final FieldKeplerianOrbit<T> mean = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(converter.convertToMean(state.getOrbit()));
  859.         final FieldBLModel<T> newModel = new FieldBLModel<>(mean, state.getMass(), referenceRadius, mu, ck0);
  860.         if (forward) {
  861.             models.addValidAfter(newModel, state.getDate());
  862.         } else {
  863.             models.addValidBefore(newModel, state.getDate());
  864.         }
  865.         stateChanged(state);
  866.     }

  868.     /** {@inheritDoc} */
  869.     public FieldKeplerianOrbit<T> propagateOrbit(final FieldAbsoluteDate<T> date, final T[] parameters) {
  870.         // compute Cartesian parameters, taking derivatives into account
  871.         final FieldBLModel<T> current = models.get(date);
  872.         return current.propagateParameters(date, parameters);
  873.     }

  875.     /**
  876.      * Get the value of the M2 drag parameter.
  877.      * @return the value of the M2 drag parameter
  878.      */
  879.     public double getM2() {
  880.         return M2Driver.getValue();
  881.     }

  883.     /**
  884.      * Get the value of the M2 drag parameter.
  885.      * @param date date at which the model parameters want to be known
  886.      * @return the value of the M2 drag parameter
  887.      */
  888.     public double getM2(final AbsoluteDate date) {
  889.         return M2Driver.getValue(date);
  890.     }

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

  895.         /** Constant mass. */
  896.         private final T mass;

  898.         /** Central attraction coefficient. */
  899.         private final T mu;

  901.         /** Brouwer-Lyddane mean orbit. */
  902.         private final FieldKeplerianOrbit<T> mean;

  904.         // Preprocessed values

  906.         /** Mean mean motion: n0 = √(μ/a")/a". */
  907.         private final T n0;

  909.         /** η = √(1 - e"²). */
  910.         private final T n;
  911.         /** η². */
  912.         private final T n2;
  913.         /** η³. */
  914.         private final T n3;
  915.         /** η + 1 / (1 + η). */
  916.         private final T t8;

  918.         /** Secular correction for mean anomaly l: &delta;<sub>s</sub>l. */
  919.         private final T dsl;
  920.         /** Secular correction for perigee argument g: &delta;<sub>s</sub>g. */
  921.         private final T dsg;
  922.         /** Secular correction for raan h: &delta;<sub>s</sub>h. */
  923.         private final T dsh;

  925.         /** Secular rate of change of semi-major axis due to drag. */
  926.         private final T aRate;
  927.         /** Secular rate of change of eccentricity due to drag. */
  928.         private final T eRate;

  930.         // CHECKSTYLE: stop JavadocVariable check

  932.         // Storage for speed-up
  933.         private final T yp2;
  934.         private final T ci;
  935.         private final T si;
  936.         private final T oneMci2;
  937.         private final T ci2X3M1;

  939.         // Long periodic corrections factors
  940.         private final T vle1;
  941.         private final T vle2;
  942.         private final T vle3;
  943.         private final T vli1;
  944.         private final T vli2;
  945.         private final T vli3;
  946.         private final T vll2;
  947.         private final T vlh1I;
  948.         private final T vlh2I;
  949.         private final T vlh3I;
  950.         private final T vls1;
  951.         private final T vls2;
  952.         private final T vls3;

  954.         // CHECKSTYLE: resume JavadocVariable check

  956.         /** Create a model for specified mean orbit.
  957.          * @param mean mean Fieldorbit
  958.          * @param mass constant mass
  959.          * @param referenceRadius reference radius of the central body attraction model (m)
  960.          * @param mu central attraction coefficient (m³/s²)
  961.          * @param ck0 un-normalized zonal coefficients
  962.          */
  963.         FieldBLModel(final FieldKeplerianOrbit<T> mean, final T mass,
  964.                      final double referenceRadius, final T mu, final double[] ck0) {

  966.             this.mass = mass;
  967.             this.mu   = mu;

  969.             // mean orbit
  970.             this.mean = mean;

  972.             final T one = mass.getField().getOne();

  974.             // mean eccentricity e"
  975.             final T epp = mean.getE();
  976.             if (epp.getReal() >= 1) {
  977.                 // Only for elliptical (e < 1) orbits
  978.                 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
  979.                                           epp.getReal());
  980.             }
  981.             final T epp2 = epp.square();

  983.             // η
  984.             n2 = one.subtract(epp2);
  985.             n  = n2.sqrt();
  986.             n3 = n2.multiply(n);
  987.             t8 = n.add(one.add(n).reciprocal());

  989.             // mean semi-major axis a"
  990.             final T app = mean.getA();

  992.             // mean mean motion
  993.             n0 = mu.divide(app).sqrt().divide(app);

  995.             // ae/a"
  996.             final T q = app.divide(referenceRadius).reciprocal();

  998.             // γ2'
  999.             T ql = q.square();
  1000.             T nl = n2.square();
  1001.             yp2 = ql.multiply(-0.5 * ck0[2]).divide(nl);
  1002.             final T yp22 = yp2.square();

  1004.             // γ3'
  1005.             ql = ql.multiply(q);
  1006.             nl = nl.multiply(n2);
  1007.             final T yp3 = ql.multiply(ck0[3]).divide(nl);

  1009.             // γ4'
  1010.             ql = ql.multiply(q);
  1011.             nl = nl.multiply(n2);
  1012.             final T yp4 = ql.multiply(0.375 * ck0[4]).divide(nl);

  1014.             // γ5'
  1015.             ql = ql.multiply(q);
  1016.             nl = nl.multiply(n2);
  1017.             final T yp5 = ql.multiply(ck0[5]).divide(nl);

  1019.             // mean inclination I" sin & cos
  1020.             final FieldSinCos<T> sc = FastMath.sinCos(mean.getI());
  1021.             si = sc.sin();
  1022.             ci = sc.cos();
  1023.             final T ci2 = ci.square();
  1024.             oneMci2 = one.subtract(ci2);
  1025.             ci2X3M1 = ci2.multiply(3.).subtract(one);
  1026.             final T ci2X5M1 = ci2.multiply(5.).subtract(one);

  1028.             // secular corrections
  1029.             // true anomaly
  1030.             final T dsl1  = yp2.multiply(n).multiply(1.5);
  1031.             final T dsl2a = n.multiply(n.multiply(25.).add(16.)).subtract(15.);
  1032.             final T dsl2b = n.multiply(n.multiply(90.).add(96.)).negate().add(30.);
  1033.             final T dsl2c = n.multiply(n.multiply(25.).add(144.)).add(105.);
  1034.             final T dsl21 = dsl2a.add(ci2.multiply(dsl2b.add(ci2.multiply(dsl2c))));
  1035.             final T dsl2  = ci2X3M1.add(yp2.multiply(0.0625).multiply(dsl21));
  1036.             final T dsl3  = yp4.multiply(n).multiply(epp2).multiply(0.9375).
  1037.                                 multiply(ci2.multiply(35.0).subtract(30.0).multiply(ci2).add(3.));
  1038.             dsl = dsl1.multiply(dsl2).add(dsl3);

  1040.             // perigee argument
  1041.             final T dsg1  = yp2.multiply(1.5).multiply(ci2X5M1);
  1042.             final T dsg2a = n.multiply(25.).add(24.).multiply(n).add(-35.);
  1043.             final T dsg2b = n.multiply(126.).add(192.).multiply(n).negate().add(90.);
  1044.             final T dsg2c = n.multiply(45.).add(360.).multiply(n).add(385.);
  1045.             final T dsg21 = dsg2a.add(ci2.multiply(dsg2b.add(ci2.multiply(dsg2c))));
  1046.             final T dsg2  = yp22.multiply(0.09375).multiply(dsg21);
  1047.             final T dsg3a = n2.multiply(-9.).add(21.);
  1048.             final T dsg3b = n2.multiply(126.).add(-270.);
  1049.             final T dsg3c = n2.multiply(-189.).add(385.);
  1050.             final T dsg31 = dsg3a.add(ci2.multiply(dsg3b.add(ci2.multiply(dsg3c))));
  1051.             final T dsg3  = yp4.multiply(0.3125).multiply(dsg31);
  1052.             dsg = dsg1.add(dsg2).add(dsg3);

  1054.             // right ascension of ascending node
  1055.             final T dsh1  = yp2.multiply(-3.);
  1056.             final T dsh2a = n.multiply(9.).add(12.).multiply(n).add(-5.);
  1057.             final T dsh2b = n.multiply(5.).add(36.).multiply(n).add(35.);
  1058.             final T dsh21 = dsh2a.subtract(ci2.multiply(dsh2b));
  1059.             final T dsh2  = yp22.multiply(0.375).multiply(dsh21);
  1060.             final T dsh31 = n2.multiply(3.).subtract(5.);
  1061.             final T dsh32 = ci2.multiply(7.).subtract(3.);
  1062.             final T dsh3  = yp4.multiply(1.25).multiply(dsh31).multiply(dsh32);
  1063.             dsh = ci.multiply(dsh1.add(dsh2).add(dsh3));

  1065.             // secular rates of change due to drag
  1066.             // Eq. 2.41 and Eq. 2.45 of Phipps' 1992 thesis
  1067.             final T coef = n0.multiply(one.add(dsl)).multiply(3.).reciprocal().multiply(-4);
  1068.             aRate = coef.multiply(app);
  1069.             eRate = coef.multiply(epp).multiply(n2);

  1071.             // singular term 1/(1 - 5 * cos²(I")) replaced by T2 function
  1072.             final T t2 = T2(ci);

  1074.             // factors for long periodic corrections
  1075.             final T fs12 = yp3.divide(yp2);
  1076.             final T fs13 = yp4.multiply(10).divide(yp2.multiply(3));
  1077.             final T fs14 = yp5.divide(yp2);

  1079.             final T ci2Xt2 = ci2.multiply(t2);
  1080.             final T cA = one.subtract(ci2.multiply(ci2Xt2.multiply(40.) .add(11.)));
  1081.             final T cB = one.subtract(ci2.multiply(ci2Xt2.multiply(8.)  .add(3.)));
  1082.             final T cC = one.subtract(ci2.multiply(ci2Xt2.multiply(24.) .add(9.)));
  1083.             final T cD = one.subtract(ci2.multiply(ci2Xt2.multiply(16.) .add(5.)));
  1084.             final T cE = one.subtract(ci2.multiply(ci2Xt2.multiply(200.).add(33.)));
  1085.             final T cF = one.subtract(ci2.multiply(ci2Xt2.multiply(40.) .add(9.)));

  1087.             final T p5p   = one.add(ci2Xt2.multiply(ci2Xt2.multiply(20.).add(8.)));
  1088.             final T p5p2  = one.add(p5p.multiply(2.));
  1089.             final T p5p4  = one.add(p5p.multiply(4.));
  1090.             final T p5p10 = one.add(p5p.multiply(10.));

  1092.             final T e2X3P4  = epp2.multiply(3.).add(4.);
  1093.             final T ciO1Pci = ci.divide(one.add(ci));
  1094.             final T oneMci  = one.subtract(ci);

  1096.             final T q1 = (yp2.multiply(cA).subtract(fs13.multiply(cB))).
  1097.                             multiply(0.125);
  1098.             final T q2 = (yp2.multiply(p5p10).subtract(fs13.multiply(p5p2))).
  1099.                             multiply(epp2).multiply(ci).multiply(0.125);
  1100.             final T q5 = (fs12.add(e2X3P4.multiply(fs14).multiply(cC).multiply(0.3125))).
  1101.                             multiply(0.25);
  1102.             final T p2 = p5p2.multiply(epp).multiply(ci).multiply(si).multiply(e2X3P4).multiply(fs14).
  1103.                             multiply(0.46875);
  1104.             final T p3 = epp.multiply(si).multiply(fs14).multiply(cC).
  1105.                             multiply(0.15625);
  1106.             final double kf = 35. / 1152.;
  1107.             final T p4 = epp.multiply(fs14).multiply(cD).
  1108.                             multiply(kf);
  1109.             final T p5 = epp.multiply(epp2).multiply(ci).multiply(si).multiply(fs14).multiply(p5p4).
  1110.                             multiply(2. * kf);

  1112.             vle1 = epp.multiply(n2).multiply(q1);
  1113.             vle2 = n2.multiply(si).multiply(q5);
  1114.             vle3 = epp.multiply(n2).multiply(si).multiply(p4).multiply(-3.0);

  1116.             vli1 = epp.multiply(q1).divide(si).negate();
  1117.             vli2 = epp.multiply(ci).multiply(q5).negate();
  1118.             vli3 = epp2.multiply(ci).multiply(p4).multiply(-3.0);

  1120.             vll2 = vle2.add(epp.multiply(n2).multiply(p3).multiply(3.0));

  1122.             vlh1I = si.multiply(q2).negate();
  1123.             vlh2I = epp.multiply(ci).multiply(q5).add(si.multiply(p2));
  1124.             vlh3I = (epp2.multiply(ci).multiply(p4).add(si.multiply(p5))).negate();

  1126.             vls1 = q1.multiply(n3.subtract(one)).
  1127.                    subtract(q2).
  1128.                    add(epp2.multiply(ci2).multiply(ci2Xt2).multiply(ci2Xt2).
  1129.                        multiply(yp2.subtract(fs13.multiply(0.2))).multiply(25.0)).
  1130.                    subtract(epp2.multiply(yp2.multiply(cE).subtract(fs13.multiply(cF))).multiply(0.0625));

  1132.             vls2 = epp.multiply(si).multiply(t8.add(ciO1Pci)).multiply(q5).
  1133.                    add((epp2.subtract(n3).multiply(3.).add(11.)).multiply(p3)).
  1134.                    add(oneMci.multiply(p2));

  1136.             vls3 = si.multiply(p4).multiply(n3.subtract(one).multiply(3.).
  1137.                                             subtract(epp2.multiply(ciO1Pci.add(2.)))).
  1138.                    subtract(oneMci.multiply(p5));
  1139.         }

  1141.         /**
  1142.          * Get true anomaly from mean anomaly.
  1143.          * @param lM  the mean anomaly (rad)
  1144.          * @param ecc the eccentricity
  1145.          * @return the true anomaly (rad)
  1146.          */
  1147.         private FieldUnivariateDerivative1<T> getTrueAnomaly(final FieldUnivariateDerivative1<T> lM,
  1148.                                                              final FieldUnivariateDerivative1<T> ecc) {

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

  1152.             // reduce M to [-PI PI] interval
  1153.             final FieldUnivariateDerivative1<T> reducedM = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(lM.getValue(), zero),
  1154.                                                                                             lM.getFirstDerivative());

  1156.             // compute the true anomaly
  1157.             FieldUnivariateDerivative1<T> lV = FieldKeplerianAnomalyUtility.ellipticMeanToTrue(ecc, lM);

  1159.             // expand the result back to original range
  1160.             lV = lV.add(lM.getValue().subtract(reducedM.getValue()));

  1162.             // Returns the true anomaly
  1163.             return lV;
  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²(i))<sup>-1</sup> (causing the singularity)
  1172.          * by a function, named T2 in the thesis.
  1173.          * </p>
  1174.          * @param cosI cosine of the mean inclination
  1175.          * @return an approximation of (1.0 - 5.0 * cos²(i))<sup>-1</sup> term
  1176.          */
  1177.         private T T2(final T cosI) {

  1179.             // X = (1.0 - 5.0 * cos²(i))
  1180.             final T x  = cosI.square().multiply(-5.0).add(1.0);
  1181.             final T x2 = x.square();
  1182.             final T xb = x2.multiply(BETA);

  1184.             // Eq. 2.48
  1185.             T sum = x.getField().getZero();
  1186.             for (int i = 0; i <= 12; i++) {
  1187.                 final double sign = i % 2 == 0 ? +1.0 : -1.0;
  1188.                 sum = sum.add(FastMath.pow(x2, i).
  1189.                               multiply(FastMath.pow(BETA, i)).
  1190.                               multiply(sign).
  1191.                               divide(CombinatoricsUtils.factorialDouble(i + 1)));
  1192.             }

  1194.             // Right term of equation 2.47
  1195.             final T one = x.getField().getOne();
  1196.             T product = one;
  1197.             for (int i = 0; i <= 10; i++) {
  1198.                 product = product.multiply(one.add(FastMath.exp(xb.multiply(FastMath.scalb(-1.0, i)))));
  1199.             }

  1201.             // Return (Eq. 2.47)
  1202.             return x.multiply(BETA).multiply(sum).multiply(product);
  1203.         }

  1205.         /** Extrapolate an orbit up to a specific target date.
  1206.          * @param date target date for the orbit
  1207.          * @param parameters model parameters
  1208.          * @return propagated parameters
  1209.          */
  1210.         public FieldKeplerianOrbit<T> propagateParameters(final FieldAbsoluteDate<T> date, final T[] parameters) {

  1212.             // Field
  1213.             final Field<T> field = date.getField();
  1214.             final T one  = field.getOne();
  1215.             final T zero = field.getZero();

  1217.             // Empirical drag coefficient M2
  1218.             final T m2 = parameters[0];

  1220.             // Keplerian evolution
  1221.             final FieldUnivariateDerivative1<T> dt  = new FieldUnivariateDerivative1<>(date.durationFrom(mean.getDate()), one);
  1222.             final FieldUnivariateDerivative1<T> not = dt.multiply(n0);

  1224.             final FieldUnivariateDerivative1<T> dtM2  = dt.multiply(m2);
  1225.             final FieldUnivariateDerivative1<T> dt2M2 = dt.multiply(dtM2);

  1227.             // Secular corrections
  1228.             // -------------------

  1230.             // semi-major axis (with drag Eq. 2.41 of Phipps' 1992 thesis)
  1231.             final FieldUnivariateDerivative1<T> app = dtM2.multiply(aRate).add(mean.getA());

  1233.             // eccentricity  (with drag Eq. 2.45 of Phipps' 1992 thesis) reduced to [0, 1[
  1234.             final FieldUnivariateDerivative1<T> tmp = dtM2.multiply(eRate).add(mean.getE());
  1235.             final FieldUnivariateDerivative1<T> epp = FastMath.max(FastMath.min(tmp, MAX_ECC), 0.);

  1237.             // mean argument of perigee
  1238.             final T gp0 = MathUtils.normalizeAngle(mean.getPerigeeArgument().add(dsg.multiply(not.getValue())), zero);
  1239.             final T gp1 = dsg.multiply(n0);
  1240.             final FieldUnivariateDerivative1<T> gpp = new FieldUnivariateDerivative1<>(gp0, gp1);

  1242.             // mean longitude of ascending node
  1243.             final T hp0 = MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode().add(dsh.multiply(not.getValue())), zero);
  1244.             final T hp1 = dsh.multiply(n0);
  1245.             final FieldUnivariateDerivative1<T> hpp = new FieldUnivariateDerivative1<>(hp0, hp1);

  1247.             // mean anomaly (with drag Eq. 2.38 of Phipps' 1992 thesis)
  1248.             final T lp0 = MathUtils.normalizeAngle(mean.getMeanAnomaly().add(dsl.add(one).multiply(not.getValue())).add(dt2M2.getValue()), zero);
  1249.             final T lp1 = dsl.add(one).multiply(n0).add(dtM2.multiply(2.0).getValue());
  1250.             final FieldUnivariateDerivative1<T> lpp = new FieldUnivariateDerivative1<>(lp0, lp1);

  1252.             // Long period corrections
  1253.             //------------------------
  1254.             final FieldSinCos<FieldUnivariateDerivative1<T>> scgpp = gpp.sinCos();
  1255.             final FieldUnivariateDerivative1<T> cgpp = scgpp.cos();
  1256.             final FieldUnivariateDerivative1<T> sgpp = scgpp.sin();
  1257.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc2gpp = gpp.multiply(2).sinCos();
  1258.             final FieldUnivariateDerivative1<T> c2gpp  = sc2gpp.cos();
  1259.             final FieldUnivariateDerivative1<T> s2gpp  = sc2gpp.sin();
  1260.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc3gpp = gpp.multiply(3).sinCos();
  1261.             final FieldUnivariateDerivative1<T> c3gpp  = sc3gpp.cos();
  1262.             final FieldUnivariateDerivative1<T> s3gpp  = sc3gpp.sin();

  1264.             // δ1e
  1265.             final FieldUnivariateDerivative1<T> d1e = c2gpp.multiply(vle1).
  1266.                                                       add(sgpp.multiply(vle2)).
  1267.                                                       add(s3gpp.multiply(vle3));

  1269.             // δ1I
  1270.             FieldUnivariateDerivative1<T> d1I = sgpp.multiply(vli2).
  1271.                                                 add(s3gpp.multiply(vli3));
  1272.             // Pseudo singular term, not to add if I" is zero
  1273.             if (Double.isFinite(vli1.getReal())) {
  1274.                 d1I = d1I.add(c2gpp.multiply(vli1));
  1275.             }

  1277.             // e"δ1l
  1278.             final FieldUnivariateDerivative1<T> eppd1l = s2gpp.multiply(vle1).
  1279.                                                          subtract(cgpp.multiply(vll2)).
  1280.                                                          subtract(c3gpp.multiply(vle3)).
  1281.                                                          multiply(n);

  1283.             // sinI"δ1h
  1284.             final FieldUnivariateDerivative1<T> sIppd1h = s2gpp.multiply(vlh1I).
  1285.                                                           add(cgpp.multiply(vlh2I)).
  1286.                                                           add(c3gpp.multiply(vlh3I));

  1288.             // δ1z = δ1l + δ1g + δ1h
  1289.             final FieldUnivariateDerivative1<T> d1z = s2gpp.multiply(vls1).
  1290.                                                       add(cgpp.multiply(vls2)).
  1291.                                                       add(c3gpp.multiply(vls3));

  1293.             // Short period corrections
  1294.             // ------------------------

  1296.             // true anomaly
  1297.             final FieldUnivariateDerivative1<T> fpp = getTrueAnomaly(lpp, epp);
  1298.             final FieldSinCos<FieldUnivariateDerivative1<T>> scfpp = fpp.sinCos();
  1299.             final FieldUnivariateDerivative1<T> cfpp = scfpp.cos();
  1300.             final FieldUnivariateDerivative1<T> sfpp = scfpp.sin();

  1302.             // e"sin(f')
  1303.             final FieldUnivariateDerivative1<T> eppsfpp = epp.multiply(sfpp);
  1304.             // e"cos(f')
  1305.             final FieldUnivariateDerivative1<T> eppcfpp = epp.multiply(cfpp);
  1306.             // 1 + e"cos(f')
  1307.             final FieldUnivariateDerivative1<T> eppcfppP1 = eppcfpp.add(1);
  1308.             // 2 + e"cos(f')
  1309.             final FieldUnivariateDerivative1<T> eppcfppP2 = eppcfpp.add(2);
  1310.             // 3 + e"cos(f')
  1311.             final FieldUnivariateDerivative1<T> eppcfppP3 = eppcfpp.add(3);
  1312.             // (1 + e"cos(f'))³
  1313.             final FieldUnivariateDerivative1<T> eppcfppP1_3 = eppcfppP1.square().multiply(eppcfppP1);

  1315.             // 2g"
  1316.             final FieldUnivariateDerivative1<T> g2 = gpp.multiply(2);

  1318.             // 2g" + f"
  1319.             final FieldUnivariateDerivative1<T> g2f = g2.add(fpp);
  1320.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc2gf = g2f.sinCos();
  1321.             final FieldUnivariateDerivative1<T> c2gf = sc2gf.cos();
  1322.             final FieldUnivariateDerivative1<T> s2gf = sc2gf.sin();
  1323.             final FieldUnivariateDerivative1<T> eppc2gf = epp.multiply(c2gf);
  1324.             final FieldUnivariateDerivative1<T> epps2gf = epp.multiply(s2gf);

  1326.             // 2g" + 2f"
  1327.             final FieldUnivariateDerivative1<T> g2f2 = g2.add(fpp.multiply(2));
  1328.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc2g2f = g2f2.sinCos();
  1329.             final FieldUnivariateDerivative1<T> c2g2f = sc2g2f.cos();
  1330.             final FieldUnivariateDerivative1<T> s2g2f = sc2g2f.sin();

  1332.             // 2g" + 3f"
  1333.             final FieldUnivariateDerivative1<T> g2f3 = g2.add(fpp.multiply(3));
  1334.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc2g3f = g2f3.sinCos();
  1335.             final FieldUnivariateDerivative1<T> c2g3f = sc2g3f.cos();
  1336.             final FieldUnivariateDerivative1<T> s2g3f = sc2g3f.sin();

  1338.             // e"cos(2g" + 3f")
  1339.             final FieldUnivariateDerivative1<T> eppc2g3f = epp.multiply(c2g3f);
  1340.             // e"sin(2g" + 3f")
  1341.             final FieldUnivariateDerivative1<T> epps2g3f = epp.multiply(s2g3f);

  1343.             // f" + e"sin(f") - l"
  1344.             final FieldUnivariateDerivative1<T> w17 = fpp.add(eppsfpp).subtract(lpp);

  1346.             // ((e"cos(f") + 3)e"cos(f") + 3)cos(f")
  1347.             final FieldUnivariateDerivative1<T> w20 = cfpp.multiply(eppcfppP3.multiply(eppcfpp).add(3.));

  1349.             // 3sin(2g" + 2f") + 3e"sin(2g" + f") + e"sin(2g" + f")
  1350.             final FieldUnivariateDerivative1<T> w21 = s2g2f.add(epps2gf).multiply(3).add(epps2g3f);

  1352.             // (1 + e"cos(f"))(2 + e"cos(f"))/η²
  1353.             final FieldUnivariateDerivative1<T> w22 = eppcfppP1.multiply(eppcfppP2).divide(n2);

  1355.             // sinCos(I"/2)
  1356.             final FieldSinCos<T> sci = FastMath.sinCos(mean.getI().divide(2.));
  1357.             final T siO2 = sci.sin();
  1358.             final T ciO2 = sci.cos();

  1360.             // δ2a
  1361.             final FieldUnivariateDerivative1<T> d2a = app.multiply(yp2).divide(n2).
  1362.                                                       multiply(eppcfppP1_3.subtract(n3).multiply(ci2X3M1).
  1363.                                                                add(c2g2f.multiply(eppcfppP1_3).multiply(oneMci2).multiply(3.)));

  1365.             // δ2e
  1366.             final FieldUnivariateDerivative1<T> d2e = (w20.add(epp.multiply(t8))).multiply(ci2X3M1).
  1367.                                                        add((w20.add(epp.multiply(c2g2f))).multiply(oneMci2.multiply(3))).
  1368.                                                        subtract((eppc2gf.multiply(3).add(eppc2g3f)).multiply(oneMci2.multiply(n2))).
  1369.                                                       multiply(yp2.multiply(0.5));

  1371.             // δ2I
  1372.             final FieldUnivariateDerivative1<T> d2I = ((c2g2f.add(eppc2gf)).multiply(3).add(eppc2g3f)).
  1373.                                                        multiply(yp2.divide(2.).multiply(ci).multiply(si));

  1375.             // e"δ2l
  1376.             final FieldUnivariateDerivative1<T> eppd2l = (w22.add(1).multiply(sfpp).multiply(oneMci2).multiply(2.).
  1377.                                                          add((w22.subtract(1).negate().multiply(s2gf)).
  1378.                                                               add(w22.add(1. / 3.).multiply(s2g3f)).
  1379.                                                              multiply(oneMci2.multiply(3.)))).
  1380.                                                         multiply(yp2.divide(4.).multiply(n3)).negate();

  1382.             // sinI"δ2h
  1383.             final FieldUnivariateDerivative1<T> sIppd2h = (w21.subtract(w17.multiply(6))).
  1384.                                                            multiply(yp2).multiply(ci).multiply(si).divide(2.);

  1386.             // δ2z = δ2l + δ2g + δ2h
  1387.             final T ttt = one.add(ci.multiply(ci.multiply(-5).add(2.)));
  1388.             final FieldUnivariateDerivative1<T> d2z = (epp.multiply(eppd2l).multiply(t8.subtract(one)).divide(n3).
  1389.                                                        add(w17.multiply(ttt).multiply(6).subtract(w21.multiply(ttt.add(2.))).
  1390.                                                            multiply(yp2.divide(4.)))).
  1391.                                                        negate();

  1393.             // Assembling elements
  1394.             // -------------------

  1396.             // e" + δe
  1397.             final FieldUnivariateDerivative1<T> de = epp.add(d1e).add(d2e);

  1399.             // e"δl
  1400.             final FieldUnivariateDerivative1<T> dl = eppd1l.add(eppd2l);

  1402.             // sin(I"/2)δh = sin(I")δh / cos(I"/2) (singular for I" = π, very unlikely)
  1403.             final FieldUnivariateDerivative1<T> dh = sIppd1h.add(sIppd2h).divide(ciO2.multiply(2.));

  1405.             // δI
  1406.             final FieldUnivariateDerivative1<T> di = d1I.add(d2I).multiply(ciO2).divide(2.).add(siO2);

  1408.             // z = l" + g" + h" + δ1z + δ2z
  1409.             final FieldUnivariateDerivative1<T> z = lpp.add(gpp).add(hpp).add(d1z).add(d2z);

  1411.             // Osculating elements
  1412.             // -------------------

  1414.             // Semi-major axis
  1415.             final FieldUnivariateDerivative1<T> a = app.add(d2a);

  1417.             // Eccentricity
  1418.             final FieldUnivariateDerivative1<T> e = FastMath.sqrt(de.square().add(dl.square()));

  1420.             // Mean anomaly
  1421.             final FieldSinCos<FieldUnivariateDerivative1<T>> sclpp = lpp.sinCos();
  1422.             final FieldUnivariateDerivative1<T> clpp = sclpp.cos();
  1423.             final FieldUnivariateDerivative1<T> slpp = sclpp.sin();
  1424.             final FieldUnivariateDerivative1<T> l = FastMath.atan2(de.multiply(slpp).add(dl.multiply(clpp)),
  1425.                                                                    de.multiply(clpp).subtract(dl.multiply(slpp)));

  1427.             // Inclination
  1428.             final FieldUnivariateDerivative1<T> i = FastMath.acos(di.square().add(dh.square()).multiply(2).negate().add(1.));

  1430.             // Longitude of ascending node
  1431.             final FieldSinCos<FieldUnivariateDerivative1<T>> schpp = hpp.sinCos();
  1432.             final FieldUnivariateDerivative1<T> chpp = schpp.cos();
  1433.             final FieldUnivariateDerivative1<T> shpp = schpp.sin();
  1434.             final FieldUnivariateDerivative1<T> h = FastMath.atan2(di.multiply(shpp).add(dh.multiply(chpp)),
  1435.                                                                    di.multiply(chpp).subtract(dh.multiply(shpp)));

  1437.             // Argument of perigee
  1438.             final FieldUnivariateDerivative1<T> g = z.subtract(l).subtract(h);

  1440.             // Return a Keplerian orbit
  1441.             return new FieldKeplerianOrbit<>(a.getValue(), e.getValue(), i.getValue(),
  1442.                                              g.getValue(), h.getValue(), l.getValue(),
  1443.                                              a.getFirstDerivative(), e.getFirstDerivative(), i.getFirstDerivative(),
  1444.                                              g.getFirstDerivative(), h.getFirstDerivative(), l.getFirstDerivative(),
  1445.                                              PositionAngleType.MEAN, mean.getFrame(), date, this.mu);
  1446.         }
  1447.     }

  1449.     /** {@inheritDoc} */
  1450.     @Override
  1451.     protected T getMass(final FieldAbsoluteDate<T> date) {
  1452.         return models.get(date).mass;
  1453.     }

  1455.     /** {@inheritDoc} */
  1456.     @Override
  1457.     public List<ParameterDriver> getParametersDrivers() {
  1458.         return Collections.singletonList(M2Driver);
  1459.     }

  1461. }
Created with JaCoCo 0.8.12.202403310830