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.FieldEquinoctialOrbit;
  36. import org.orekit.orbits.FieldKeplerianAnomalyUtility;
  37. import org.orekit.orbits.FieldKeplerianOrbit;
  38. import org.orekit.orbits.FieldOrbit;
  39. import org.orekit.orbits.OrbitType;
  40. import org.orekit.orbits.PositionAngleType;
  41. import org.orekit.propagation.FieldSpacecraftState;
  42. import org.orekit.propagation.PropagationType;
  43. import org.orekit.propagation.analytical.tle.FieldTLE;
  44. import org.orekit.time.AbsoluteDate;
  45. import org.orekit.time.FieldAbsoluteDate;
  46. import org.orekit.utils.FieldTimeSpanMap;
  47. import org.orekit.utils.ParameterDriver;

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

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

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

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

  109.     /** Initial Brouwer-Lyddane model. */
  110.     private FieldBLModel<T> initialModel;

  112.     /** All models. */
  113.     private transient FieldTimeSpanMap<FieldBLModel<T>, T> models;

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

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

  121.     /** Un-normalized zonal coefficients. */
  122.     private double[] ck0;

  124.     /** Empirical coefficient used in the drag modeling. */
  125.     private final ParameterDriver M2Driver;

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

  518.         // store model coefficients
  519.         this.referenceRadius = referenceRadius;
  520.         this.mu  = mu;
  521.         this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50};

  523.         // initialize M2 driver
  524.         this.M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, M2, SCALE,
  525.                                             Double.NEGATIVE_INFINITY,
  526.                                             Double.POSITIVE_INFINITY);

  528.         // compute mean parameters if needed
  529.         resetInitialState(new FieldSpacecraftState<>(initialOrbit,
  530.                                                      attitudeProv.getAttitude(initialOrbit,
  531.                                                                               initialOrbit.getDate(),
  532.                                                                               initialOrbit.getFrame()),
  533.                                                      mass),
  534.                                                      initialType, epsilon, maxIterations);

  536.     }

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

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

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

  659.     /** {@inheritDoc}
  660.      * <p>The new initial state to consider
  661.      * must be defined with an osculating orbit.</p>
  662.      * @see #resetInitialState(FieldSpacecraftState, PropagationType)
  663.      */
  664.     @Override
  665.     public void resetInitialState(final FieldSpacecraftState<T> state) {
  666.         resetInitialState(state, PropagationType.OSCULATING);
  667.     }

  669.     /** Reset the propagator initial state.
  670.      * @param state new initial state to consider
  671.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  672.      */
  673.     public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType) {
  674.         resetInitialState(state, stateType,
  675.                           BrouwerLyddanePropagator.EPSILON_DEFAULT,
  676.                           BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT);
  677.     }

  679.     /** Reset the propagator initial state.
  680.      * @param state new initial state to consider
  681.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  682.      * @param epsilon convergence threshold for mean parameters conversion
  683.      * @param maxIterations maximum iterations for mean parameters conversion
  684.      * @since 11.2
  685.      */
  686.     public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType,
  687.                                   final double epsilon, final int maxIterations) {
  688.         super.resetInitialState(state);
  689.         final FieldKeplerianOrbit<T> keplerian = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(state.getOrbit());
  690.         this.initialModel = (stateType == PropagationType.MEAN) ?
  691.                              new FieldBLModel<>(keplerian, state.getMass(), referenceRadius, mu, ck0) :
  692.                              computeMeanParameters(keplerian, state.getMass(), epsilon, maxIterations);
  693.         this.models = new FieldTimeSpanMap<>(initialModel, state.getA().getField());
  694.     }

  696.     /** {@inheritDoc} */
  697.     @Override
  698.     protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward) {
  699.         resetIntermediateState(state, forward,
  700.                                BrouwerLyddanePropagator.EPSILON_DEFAULT,
  701.                                BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT);
  702.     }

  704.     /** Reset an intermediate state.
  705.      * @param state new intermediate state to consider
  706.      * @param forward if true, the intermediate state is valid for
  707.      * propagations after itself
  708.      * @param epsilon convergence threshold for mean parameters conversion
  709.      * @param maxIterations maximum iterations for mean parameters conversion
  710.      * @since 11.2
  711.      */
  712.     protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward,
  713.                                           final double epsilon, final int maxIterations) {
  714.         final FieldBLModel<T> newModel = computeMeanParameters((FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(state.getOrbit()),
  715.                                                                state.getMass(), epsilon, maxIterations);
  716.         if (forward) {
  717.             models.addValidAfter(newModel, state.getDate());
  718.         } else {
  719.             models.addValidBefore(newModel, state.getDate());
  720.         }
  721.         stateChanged(state);
  722.     }

  724.     /** Compute mean parameters according to the Brouwer-Lyddane analytical model computation
  725.      * in order to do the propagation.
  726.      * @param osculating osculating orbit
  727.      * @param mass constant mass
  728.      * @param epsilon convergence threshold for mean parameters conversion
  729.      * @param maxIterations maximum iterations for mean parameters conversion
  730.      * @return Brouwer-Lyddane mean model
  731.      */
  732.     private FieldBLModel<T> computeMeanParameters(final FieldKeplerianOrbit<T> osculating, final T mass,
  733.                                                   final double epsilon, final int maxIterations) {

  735.         // damping factor
  736.         final double damping = BrouwerLyddanePropagator.DAMPING_DEFAULT;

  738.         // sanity check
  739.         if (osculating.getA().getReal() < referenceRadius) {
  740.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
  741.                                            osculating.getA());
  742.         }

  744.         final Field<T> field = osculating.getDate().getField();
  745.         final T zero = field.getZero();
  746.         final T pi   = zero.getPi();

  748.         // rough initialization of the mean parameters
  749.         FieldBLModel<T> current = new FieldBLModel<>(osculating, mass, referenceRadius, mu, ck0);

  751.         // Get equinoctial parameters
  752.         T sma = osculating.getA();
  753.         T ex  = osculating.getEquinoctialEx();
  754.         T ey  = osculating.getEquinoctialEy();
  755.         T hx  = osculating.getHx();
  756.         T hy  = osculating.getHy();
  757.         T lv  = osculating.getLv();

  759.         // threshold for each parameter
  760.         final T thresholdA  = osculating.getA().abs().add(1).multiply(epsilon);
  761.         final T thresholdE  = FastMath.hypot(ex, ey).add(1).multiply(epsilon);
  762.         final T thresholdH  = FastMath.hypot(hx, hy).add(1).multiply(epsilon);
  763.         final T thresholdLv = pi.multiply(epsilon);

  765.         int i = 0;
  766.         while (i++ < maxIterations) {

  768.             // recompute the osculating parameters from the current mean parameters
  769.             final FieldKeplerianOrbit<T> parameters = current.propagateParameters(current.mean.getDate(),
  770.                                                                                   getParameters(field,
  771.                                                                                                 current.mean.getDate()));

  773.             // adapted parameters residuals
  774.             final T deltaA  = osculating.getA().subtract(parameters.getA());
  775.             final T deltaEx = osculating.getEquinoctialEx().subtract(parameters.getEquinoctialEx());
  776.             final T deltaEy = osculating.getEquinoctialEy().subtract(parameters.getEquinoctialEy());
  777.             final T deltaHx = osculating.getHx().subtract(parameters.getHx());
  778.             final T deltaHy = osculating.getHy().subtract(parameters.getHy());
  779.             final T deltaLv = MathUtils.normalizeAngle(osculating.getLv().subtract(parameters.getLv()), zero);

  781.             // update state
  782.             sma = sma.add(deltaA.multiply(damping));
  783.             ex  = ex.add(deltaEx.multiply(damping));
  784.             ey  = ey.add(deltaEy.multiply(damping));
  785.             hx  = hx.add(deltaHx.multiply(damping));
  786.             hy  = hy.add(deltaHy.multiply(damping));
  787.             lv  = lv.add(deltaLv.multiply(damping));

  789.             // Update mean orbit
  790.             final FieldEquinoctialOrbit<T> meanEquinoctial =
  791.                             new FieldEquinoctialOrbit<>(sma, ex, ey, hx, hy, lv, PositionAngleType.TRUE,
  792.                                                         osculating.getFrame(), osculating.getDate(),
  793.                                                         osculating.getMu());
  794.             final FieldKeplerianOrbit<T> meanKeplerian =
  795.                             (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(meanEquinoctial);

  797.             // update mean parameters
  798.             current = new FieldBLModel<>(meanKeplerian, mass, referenceRadius, mu, ck0);

  800.             // check convergence
  801.             if (FastMath.abs(deltaA.getReal())  < thresholdA.getReal() &&
  802.                 FastMath.abs(deltaEx.getReal()) < thresholdE.getReal() &&
  803.                 FastMath.abs(deltaEy.getReal()) < thresholdE.getReal() &&
  804.                 FastMath.abs(deltaHx.getReal()) < thresholdH.getReal() &&
  805.                 FastMath.abs(deltaHy.getReal()) < thresholdH.getReal() &&
  806.                 FastMath.abs(deltaLv.getReal()) < thresholdLv.getReal()) {
  807.                 return current;
  808.             }
  809.         }
  810.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_BROUWER_LYDDANE_MEAN_PARAMETERS, i);
  811.     }


  814.     /** {@inheritDoc} */
  815.     public FieldKeplerianOrbit<T> propagateOrbit(final FieldAbsoluteDate<T> date, final T[] parameters) {
  816.         // compute Cartesian parameters, taking derivatives into account
  817.         final FieldBLModel<T> current = models.get(date);
  818.         return current.propagateParameters(date, parameters);
  819.     }

  821.     /**
  822.      * Get the value of the M2 drag parameter.
  823.      * @return the value of the M2 drag parameter
  824.      */
  825.     public double getM2() {
  826.         return M2Driver.getValue();
  827.     }

  829.     /**
  830.      * Get the value of the M2 drag parameter.
  831.      * @param date date at which the model parameters want to be known
  832.      * @return the value of the M2 drag parameter
  833.      */
  834.     public double getM2(final AbsoluteDate date) {
  835.         return M2Driver.getValue(date);
  836.     }

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

  841.         /** Constant mass. */
  842.         private final T mass;

  844.         /** Central attraction coefficient. */
  845.         private final T mu;

  847.         /** Brouwer-Lyddane mean orbit. */
  848.         private final FieldKeplerianOrbit<T> mean;

  850.         // Preprocessed values

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

  855.         /** η = √(1 - e"²). */
  856.         private final T n;
  857.         /** η². */
  858.         private final T n2;
  859.         /** η³. */
  860.         private final T n3;
  861.         /** η + 1 / (1 + η). */
  862.         private final T t8;

  864.         /** Secular correction for mean anomaly l: &delta;<sub>s</sub>l. */
  865.         private final T dsl;
  866.         /** Secular correction for perigee argument g: &delta;<sub>s</sub>g. */
  867.         private final T dsg;
  868.         /** Secular correction for raan h: &delta;<sub>s</sub>h. */
  869.         private final T dsh;

  871.         /** Secular rate of change of semi-major axis due to drag. */
  872.         private final T aRate;
  873.         /** Secular rate of change of eccentricity due to drag. */
  874.         private final T eRate;

  876.         // CHECKSTYLE: stop JavadocVariable check

  878.         // Storage for speed-up
  879.         private final T yp2;
  880.         private final T ci;
  881.         private final T si;
  882.         private final T oneMci2;
  883.         private final T ci2X3M1;

  885.         // Long periodic corrections factors
  886.         private final T vle1;
  887.         private final T vle2;
  888.         private final T vle3;
  889.         private final T vli1;
  890.         private final T vli2;
  891.         private final T vli3;
  892.         private final T vll2;
  893.         private final T vlh1I;
  894.         private final T vlh2I;
  895.         private final T vlh3I;
  896.         private final T vls1;
  897.         private final T vls2;
  898.         private final T vls3;

  900.         // CHECKSTYLE: resume JavadocVariable check

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

  912.             this.mass = mass;
  913.             this.mu   = mu;

  915.             // mean orbit
  916.             this.mean = mean;

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

  920.             // mean eccentricity e"
  921.             final T epp = mean.getE();
  922.             if (epp.getReal() >= 1) {
  923.                 // Only for elliptical (e < 1) orbits
  924.                 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
  925.                                           epp.getReal());
  926.             }
  927.             final T epp2 = epp.square();

  929.             // η
  930.             n2 = one.subtract(epp2);
  931.             n  = n2.sqrt();
  932.             n3 = n2.multiply(n);
  933.             t8 = n.add(one.add(n).reciprocal());

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

  938.             // mean mean motion
  939.             n0 = mu.divide(app).sqrt().divide(app);

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

  944.             // γ2'
  945.             T ql = q.square();
  946.             T nl = n2.square();
  947.             yp2 = ql.multiply(-0.5 * ck0[2]).divide(nl);
  948.             final T yp22 = yp2.square();

  950.             // γ3'
  951.             ql = ql.multiply(q);
  952.             nl = nl.multiply(n2);
  953.             final T yp3 = ql.multiply(ck0[3]).divide(nl);

  955.             // γ4'
  956.             ql = ql.multiply(q);
  957.             nl = nl.multiply(n2);
  958.             final T yp4 = ql.multiply(0.375 * ck0[4]).divide(nl);

  960.             // γ5'
  961.             ql = ql.multiply(q);
  962.             nl = nl.multiply(n2);
  963.             final T yp5 = ql.multiply(ck0[5]).divide(nl);

  965.             // mean inclination I" sin & cos
  966.             final FieldSinCos<T> sc = FastMath.sinCos(mean.getI());
  967.             si = sc.sin();
  968.             ci = sc.cos();
  969.             final T ci2 = ci.square();
  970.             oneMci2 = one.subtract(ci2);
  971.             ci2X3M1 = ci2.multiply(3.).subtract(one);
  972.             final T ci2X5M1 = ci2.multiply(5.).subtract(one);

  974.             // secular corrections
  975.             // true anomaly
  976.             final T dsl1  = yp2.multiply(n).multiply(1.5);
  977.             final T dsl2a = n.multiply(n.multiply(25.).add(16.)).subtract(15.);
  978.             final T dsl2b = n.multiply(n.multiply(90.).add(96.)).negate().add(30.);
  979.             final T dsl2c = n.multiply(n.multiply(25.).add(144.)).add(105.);
  980.             final T dsl21 = dsl2a.add(ci2.multiply(dsl2b.add(ci2.multiply(dsl2c))));
  981.             final T dsl2  = ci2X3M1.add(yp2.multiply(0.0625).multiply(dsl21));
  982.             final T dsl3  = yp4.multiply(n).multiply(epp2).multiply(0.9375).
  983.                                 multiply(ci2.multiply(35.0).subtract(30.0).multiply(ci2).add(3.));
  984.             dsl = dsl1.multiply(dsl2).add(dsl3);

  986.             // perigee argument
  987.             final T dsg1  = yp2.multiply(1.5).multiply(ci2X5M1);
  988.             final T dsg2a = n.multiply(25.).add(24.).multiply(n).add(-35.);
  989.             final T dsg2b = n.multiply(126.).add(192.).multiply(n).negate().add(90.);
  990.             final T dsg2c = n.multiply(45.).add(360.).multiply(n).add(385.);
  991.             final T dsg21 = dsg2a.add(ci2.multiply(dsg2b.add(ci2.multiply(dsg2c))));
  992.             final T dsg2  = yp22.multiply(0.09375).multiply(dsg21);
  993.             final T dsg3a = n2.multiply(-9.).add(21.);
  994.             final T dsg3b = n2.multiply(126.).add(-270.);
  995.             final T dsg3c = n2.multiply(-189.).add(385.);
  996.             final T dsg31 = dsg3a.add(ci2.multiply(dsg3b.add(ci2.multiply(dsg3c))));
  997.             final T dsg3  = yp4.multiply(0.3125).multiply(dsg31);
  998.             dsg = dsg1.add(dsg2).add(dsg3);

  1000.             // right ascension of ascending node
  1001.             final T dsh1  = yp2.multiply(-3.);
  1002.             final T dsh2a = n.multiply(9.).add(12.).multiply(n).add(-5.);
  1003.             final T dsh2b = n.multiply(5.).add(36.).multiply(n).add(35.);
  1004.             final T dsh21 = dsh2a.subtract(ci2.multiply(dsh2b));
  1005.             final T dsh2  = yp22.multiply(0.375).multiply(dsh21);
  1006.             final T dsh31 = n2.multiply(3.).subtract(5.);
  1007.             final T dsh32 = ci2.multiply(7.).subtract(3.);
  1008.             final T dsh3  = yp4.multiply(1.25).multiply(dsh31).multiply(dsh32);
  1009.             dsh = ci.multiply(dsh1.add(dsh2).add(dsh3));

  1011.             // secular rates of change due to drag
  1012.             // Eq. 2.41 and Eq. 2.45 of Phipps' 1992 thesis
  1013.             final T coef = n0.multiply(one.add(dsl)).multiply(3.).reciprocal().multiply(-4);
  1014.             aRate = coef.multiply(app);
  1015.             eRate = coef.multiply(epp).multiply(n2);

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

  1020.             // factors for long periodic corrections
  1021.             final T fs12 = yp3.divide(yp2);
  1022.             final T fs13 = yp4.multiply(10).divide(yp2.multiply(3));
  1023.             final T fs14 = yp5.divide(yp2);

  1025.             final T ci2Xt2 = ci2.multiply(t2);
  1026.             final T cA = one.subtract(ci2.multiply(ci2Xt2.multiply(40.) .add(11.)));
  1027.             final T cB = one.subtract(ci2.multiply(ci2Xt2.multiply(8.)  .add(3.)));
  1028.             final T cC = one.subtract(ci2.multiply(ci2Xt2.multiply(24.) .add(9.)));
  1029.             final T cD = one.subtract(ci2.multiply(ci2Xt2.multiply(16.) .add(5.)));
  1030.             final T cE = one.subtract(ci2.multiply(ci2Xt2.multiply(200.).add(33.)));
  1031.             final T cF = one.subtract(ci2.multiply(ci2Xt2.multiply(40.) .add(9.)));

  1033.             final T p5p   = one.add(ci2Xt2.multiply(ci2Xt2.multiply(20.).add(8.)));
  1034.             final T p5p2  = one.add(p5p.multiply(2.));
  1035.             final T p5p4  = one.add(p5p.multiply(4.));
  1036.             final T p5p10 = one.add(p5p.multiply(10.));

  1038.             final T e2X3P4  = epp2.multiply(3.).add(4.);
  1039.             final T ciO1Pci = ci.divide(one.add(ci));
  1040.             final T oneMci  = one.subtract(ci);

  1042.             final T q1 = (yp2.multiply(cA).subtract(fs13.multiply(cB))).
  1043.                             multiply(0.125);
  1044.             final T q2 = (yp2.multiply(p5p10).subtract(fs13.multiply(p5p2))).
  1045.                             multiply(epp2).multiply(ci).multiply(0.125);
  1046.             final T q5 = (fs12.add(e2X3P4.multiply(fs14).multiply(cC).multiply(0.3125))).
  1047.                             multiply(0.25);
  1048.             final T p2 = p5p2.multiply(epp).multiply(ci).multiply(si).multiply(e2X3P4).multiply(fs14).
  1049.                             multiply(0.46875);
  1050.             final T p3 = epp.multiply(si).multiply(fs14).multiply(cC).
  1051.                             multiply(0.15625);
  1052.             final double kf = 35. / 1152.;
  1053.             final T p4 = epp.multiply(fs14).multiply(cD).
  1054.                             multiply(kf);
  1055.             final T p5 = epp.multiply(epp2).multiply(ci).multiply(si).multiply(fs14).multiply(p5p4).
  1056.                             multiply(2. * kf);

  1058.             vle1 = epp.multiply(n2).multiply(q1);
  1059.             vle2 = n2.multiply(si).multiply(q5);
  1060.             vle3 = epp.multiply(n2).multiply(si).multiply(p4).multiply(-3.0);

  1062.             vli1 = epp.multiply(q1).divide(si).negate();
  1063.             vli2 = epp.multiply(ci).multiply(q5).negate();
  1064.             vli3 = epp2.multiply(ci).multiply(p4).multiply(-3.0);

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

  1068.             vlh1I = si.multiply(q2).negate();
  1069.             vlh2I = epp.multiply(ci).multiply(q5).add(si.multiply(p2));
  1070.             vlh3I = (epp2.multiply(ci).multiply(p4).add(si.multiply(p5))).negate();

  1072.             vls1 = q1.multiply(n3.subtract(one)).
  1073.                    subtract(q2).
  1074.                    add(epp2.multiply(ci2).multiply(ci2Xt2).multiply(ci2Xt2).
  1075.                        multiply(yp2.subtract(fs13.multiply(0.2))).multiply(25.0)).
  1076.                    subtract(epp2.multiply(yp2.multiply(cE).subtract(fs13.multiply(cF))).multiply(0.0625));

  1078.             vls2 = epp.multiply(si).multiply(t8.add(ciO1Pci)).multiply(q5).
  1079.                    add((epp2.subtract(n3).multiply(3.).add(11.)).multiply(p3)).
  1080.                    add(oneMci.multiply(p2));

  1082.             vls3 = si.multiply(p4).multiply(n3.subtract(one).multiply(3.).
  1083.                                             subtract(epp2.multiply(ciO1Pci.add(2.)))).
  1084.                    subtract(oneMci.multiply(p5));
  1085.         }

  1087.         /**
  1088.          * Get true anomaly from mean anomaly.
  1089.          * @param lM  the mean anomaly (rad)
  1090.          * @param ecc the eccentricity
  1091.          * @return the true anomaly (rad)
  1092.          */
  1093.         private FieldUnivariateDerivative1<T> getTrueAnomaly(final FieldUnivariateDerivative1<T> lM,
  1094.                                                              final FieldUnivariateDerivative1<T> ecc) {

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

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

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

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

  1108.             // Returns the true anomaly
  1109.             return lV;
  1110.         }

  1112.         /**
  1113.          * This method is used in Brouwer-Lyddane model to avoid singularity at the
  1114.          * critical inclination (i = 63.4°).
  1115.          * <p>
  1116.          * This method, based on Warren Phipps's 1992 thesis (Eq. 2.47 and 2.48),
  1117.          * approximate the factor (1.0 - 5.0 * cos²(i))<sup>-1</sup> (causing the singularity)
  1118.          * by a function, named T2 in the thesis.
  1119.          * </p>
  1120.          * @param cosI cosine of the mean inclination
  1121.          * @return an approximation of (1.0 - 5.0 * cos²(i))<sup>-1</sup> term
  1122.          */
  1123.         private T T2(final T cosI) {

  1125.             // X = (1.0 - 5.0 * cos²(i))
  1126.             final T x  = cosI.square().multiply(-5.0).add(1.0);
  1127.             final T x2 = x.square();
  1128.             final T xb = x2.multiply(BETA);

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

  1140.             // Right term of equation 2.47
  1141.             final T one = x.getField().getOne();
  1142.             T product = one;
  1143.             for (int i = 0; i <= 10; i++) {
  1144.                 product = product.multiply(one.add(FastMath.exp(xb.multiply(FastMath.scalb(-1.0, i)))));
  1145.             }

  1147.             // Return (Eq. 2.47)
  1148.             return x.multiply(BETA).multiply(sum).multiply(product);
  1149.         }

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

  1158.             // Field
  1159.             final Field<T> field = date.getField();
  1160.             final T one  = field.getOne();
  1161.             final T zero = field.getZero();

  1163.             // Empirical drag coefficient M2
  1164.             final T m2 = parameters[0];

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

  1170.             final FieldUnivariateDerivative1<T> dtM2  = dt.multiply(m2);
  1171.             final FieldUnivariateDerivative1<T> dt2M2 = dt.multiply(dtM2);

  1173.             // Secular corrections
  1174.             // -------------------

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

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

  1183.             // mean argument of perigee
  1184.             final T gp0 = MathUtils.normalizeAngle(mean.getPerigeeArgument().add(dsg.multiply(not.getValue())), zero);
  1185.             final T gp1 = dsg.multiply(n0);
  1186.             final FieldUnivariateDerivative1<T> gpp = new FieldUnivariateDerivative1<>(gp0, gp1);

  1188.             // mean longitude of ascending node
  1189.             final T hp0 = MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode().add(dsh.multiply(not.getValue())), zero);
  1190.             final T hp1 = dsh.multiply(n0);
  1191.             final FieldUnivariateDerivative1<T> hpp = new FieldUnivariateDerivative1<>(hp0, hp1);

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

  1198.             // Long period corrections
  1199.             //------------------------
  1200.             final FieldSinCos<FieldUnivariateDerivative1<T>> scgpp = gpp.sinCos();
  1201.             final FieldUnivariateDerivative1<T> cgpp = scgpp.cos();
  1202.             final FieldUnivariateDerivative1<T> sgpp = scgpp.sin();
  1203.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc2gpp = gpp.multiply(2).sinCos();
  1204.             final FieldUnivariateDerivative1<T> c2gpp  = sc2gpp.cos();
  1205.             final FieldUnivariateDerivative1<T> s2gpp  = sc2gpp.sin();
  1206.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc3gpp = gpp.multiply(3).sinCos();
  1207.             final FieldUnivariateDerivative1<T> c3gpp  = sc3gpp.cos();
  1208.             final FieldUnivariateDerivative1<T> s3gpp  = sc3gpp.sin();

  1210.             // δ1e
  1211.             final FieldUnivariateDerivative1<T> d1e = c2gpp.multiply(vle1).
  1212.                                                       add(sgpp.multiply(vle2)).
  1213.                                                       add(s3gpp.multiply(vle3));

  1215.             // δ1I
  1216.             FieldUnivariateDerivative1<T> d1I = sgpp.multiply(vli2).
  1217.                                                 add(s3gpp.multiply(vli3));
  1218.             // Pseudo singular term, not to add if I" is zero
  1219.             if (Double.isFinite(vli1.getReal())) {
  1220.                 d1I = d1I.add(c2gpp.multiply(vli1));
  1221.             }

  1223.             // e"δ1l
  1224.             final FieldUnivariateDerivative1<T> eppd1l = s2gpp.multiply(vle1).
  1225.                                                          subtract(cgpp.multiply(vll2)).
  1226.                                                          subtract(c3gpp.multiply(vle3)).
  1227.                                                          multiply(n);

  1229.             // sinI"δ1h
  1230.             final FieldUnivariateDerivative1<T> sIppd1h = s2gpp.multiply(vlh1I).
  1231.                                                           add(cgpp.multiply(vlh2I)).
  1232.                                                           add(c3gpp.multiply(vlh3I));

  1234.             // δ1z = δ1l + δ1g + δ1h
  1235.             final FieldUnivariateDerivative1<T> d1z = s2gpp.multiply(vls1).
  1236.                                                       add(cgpp.multiply(vls2)).
  1237.                                                       add(c3gpp.multiply(vls3));

  1239.             // Short period corrections
  1240.             // ------------------------

  1242.             // true anomaly
  1243.             final FieldUnivariateDerivative1<T> fpp = getTrueAnomaly(lpp, epp);
  1244.             final FieldSinCos<FieldUnivariateDerivative1<T>> scfpp = fpp.sinCos();
  1245.             final FieldUnivariateDerivative1<T> cfpp = scfpp.cos();
  1246.             final FieldUnivariateDerivative1<T> sfpp = scfpp.sin();

  1248.             // e"sin(f')
  1249.             final FieldUnivariateDerivative1<T> eppsfpp = epp.multiply(sfpp);
  1250.             // e"cos(f')
  1251.             final FieldUnivariateDerivative1<T> eppcfpp = epp.multiply(cfpp);
  1252.             // 1 + e"cos(f')
  1253.             final FieldUnivariateDerivative1<T> eppcfppP1 = eppcfpp.add(1);
  1254.             // 2 + e"cos(f')
  1255.             final FieldUnivariateDerivative1<T> eppcfppP2 = eppcfpp.add(2);
  1256.             // 3 + e"cos(f')
  1257.             final FieldUnivariateDerivative1<T> eppcfppP3 = eppcfpp.add(3);
  1258.             // (1 + e"cos(f'))³
  1259.             final FieldUnivariateDerivative1<T> eppcfppP1_3 = eppcfppP1.square().multiply(eppcfppP1);

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

  1264.             // 2g" + f"
  1265.             final FieldUnivariateDerivative1<T> g2f = g2.add(fpp);
  1266.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc2gf = g2f.sinCos();
  1267.             final FieldUnivariateDerivative1<T> c2gf = sc2gf.cos();
  1268.             final FieldUnivariateDerivative1<T> s2gf = sc2gf.sin();
  1269.             final FieldUnivariateDerivative1<T> eppc2gf = epp.multiply(c2gf);
  1270.             final FieldUnivariateDerivative1<T> epps2gf = epp.multiply(s2gf);

  1272.             // 2g" + 2f"
  1273.             final FieldUnivariateDerivative1<T> g2f2 = g2.add(fpp.multiply(2));
  1274.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc2g2f = g2f2.sinCos();
  1275.             final FieldUnivariateDerivative1<T> c2g2f = sc2g2f.cos();
  1276.             final FieldUnivariateDerivative1<T> s2g2f = sc2g2f.sin();

  1278.             // 2g" + 3f"
  1279.             final FieldUnivariateDerivative1<T> g2f3 = g2.add(fpp.multiply(3));
  1280.             final FieldSinCos<FieldUnivariateDerivative1<T>> sc2g3f = g2f3.sinCos();
  1281.             final FieldUnivariateDerivative1<T> c2g3f = sc2g3f.cos();
  1282.             final FieldUnivariateDerivative1<T> s2g3f = sc2g3f.sin();

  1284.             // e"cos(2g" + 3f")
  1285.             final FieldUnivariateDerivative1<T> eppc2g3f = epp.multiply(c2g3f);
  1286.             // e"sin(2g" + 3f")
  1287.             final FieldUnivariateDerivative1<T> epps2g3f = epp.multiply(s2g3f);

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

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

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

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

  1301.             // sinCos(I"/2)
  1302.             final FieldSinCos<T> sci = FastMath.sinCos(mean.getI().divide(2.));
  1303.             final T siO2 = sci.sin();
  1304.             final T ciO2 = sci.cos();

  1306.             // δ2a
  1307.             final FieldUnivariateDerivative1<T> d2a = app.multiply(yp2).divide(n2).
  1308.                                                       multiply(eppcfppP1_3.subtract(n3).multiply(ci2X3M1).
  1309.                                                                add(c2g2f.multiply(eppcfppP1_3).multiply(oneMci2).multiply(3.)));

  1311.             // δ2e
  1312.             final FieldUnivariateDerivative1<T> d2e = (w20.add(epp.multiply(t8))).multiply(ci2X3M1).
  1313.                                                        add((w20.add(epp.multiply(c2g2f))).multiply(oneMci2.multiply(3))).
  1314.                                                        subtract((eppc2gf.multiply(3).add(eppc2g3f)).multiply(oneMci2.multiply(n2))).
  1315.                                                       multiply(yp2.multiply(0.5));

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

  1321.             // e"δ2l
  1322.             final FieldUnivariateDerivative1<T> eppd2l = (w22.add(1).multiply(sfpp).multiply(oneMci2).multiply(2.).
  1323.                                                          add((w22.subtract(1).negate().multiply(s2gf)).
  1324.                                                               add(w22.add(1. / 3.).multiply(s2g3f)).
  1325.                                                              multiply(oneMci2.multiply(3.)))).
  1326.                                                         multiply(yp2.divide(4.).multiply(n3)).negate();

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

  1332.             // δ2z = δ2l + δ2g + δ2h
  1333.             final T ttt = one.add(ci.multiply(ci.multiply(-5).add(2.)));
  1334.             final FieldUnivariateDerivative1<T> d2z = (epp.multiply(eppd2l).multiply(t8.subtract(one)).divide(n3).
  1335.                                                        add(w17.multiply(ttt).multiply(6).subtract(w21.multiply(ttt.add(2.))).
  1336.                                                            multiply(yp2.divide(4.)))).
  1337.                                                        negate();

  1339.             // Assembling elements
  1340.             // -------------------

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

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

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

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

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

  1357.             // Osculating elements
  1358.             // -------------------

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

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

  1366.             // Mean anomaly
  1367.             final FieldSinCos<FieldUnivariateDerivative1<T>> sclpp = lpp.sinCos();
  1368.             final FieldUnivariateDerivative1<T> clpp = sclpp.cos();
  1369.             final FieldUnivariateDerivative1<T> slpp = sclpp.sin();
  1370.             final FieldUnivariateDerivative1<T> l = FastMath.atan2(de.multiply(slpp).add(dl.multiply(clpp)),
  1371.                                                                    de.multiply(clpp).subtract(dl.multiply(slpp)));

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

  1376.             // Longitude of ascending node
  1377.             final FieldSinCos<FieldUnivariateDerivative1<T>> schpp = hpp.sinCos();
  1378.             final FieldUnivariateDerivative1<T> chpp = schpp.cos();
  1379.             final FieldUnivariateDerivative1<T> shpp = schpp.sin();
  1380.             final FieldUnivariateDerivative1<T> h = FastMath.atan2(di.multiply(shpp).add(dh.multiply(chpp)),
  1381.                                                                    di.multiply(chpp).subtract(dh.multiply(shpp)));

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

  1386.             // Return a Keplerian orbit
  1387.             return new FieldKeplerianOrbit<>(a.getValue(), e.getValue(), i.getValue(),
  1388.                                              g.getValue(), h.getValue(), l.getValue(),
  1389.                                              a.getFirstDerivative(), e.getFirstDerivative(), i.getFirstDerivative(),
  1390.                                              g.getFirstDerivative(), h.getFirstDerivative(), l.getFirstDerivative(),
  1391.                                              PositionAngleType.MEAN, mean.getFrame(), date, this.mu);
  1392.         }
  1393.     }

  1395.     /** {@inheritDoc} */
  1396.     @Override
  1397.     protected T getMass(final FieldAbsoluteDate<T> date) {
  1398.         return models.get(date).mass;
  1399.     }

  1401.     /** {@inheritDoc} */
  1402.     @Override
  1403.     public List<ParameterDriver> getParametersDrivers() {
  1404.         return Collections.singletonList(M2Driver);
  1405.     }

  1407. }
Created with JaCoCo 0.8.12.202403310830