BrouwerLyddanePropagator.java

  1. /* Copyright 2002-2022 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.propagation.analytical;

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

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

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

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

  95.     /** Default value for M2 coefficient. */
  96.     public static final double M2 = 0.0;

  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.     /** Initial Brouwer-Lyddane model. */
  110.     private BLModel initialModel;

  112.     /** All models. */
  113.     private transient TimeSpanMap<BLModel> models;

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

  118.     /** Central attraction coefficient (m³/s²). */
  119.     private double 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 #M2} drag is not computed
  136.      * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, UnnormalizedSphericalHarmonicsProvider, double)
  137.      * @see #BrouwerLyddanePropagator(Orbit, UnnormalizedSphericalHarmonicsProvider, PropagationType, double)
  138.      */
  139.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  140.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  141.                                     final double M2) {
  142.         this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
  143.              DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), 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 #M2} drag is not computed
  156.      * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double,
  157.      *                                 UnnormalizedSphericalHarmonicsProvider,
  158.      *                                 UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics,
  159.      *                                 PropagationType, double)
  160.      */
  161.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  162.                                     final AttitudeProvider attitude,
  163.                                     final double mass,
  164.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  165.                                     final UnnormalizedSphericalHarmonics harmonics,
  166.                                     final double M2) {
  167.         this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
  168.              harmonics.getUnnormalizedCnm(2, 0),
  169.              harmonics.getUnnormalizedCnm(3, 0),
  170.              harmonics.getUnnormalizedCnm(4, 0),
  171.              harmonics.getUnnormalizedCnm(5, 0),
  172.              M2);
  173.     }

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

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

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

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

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

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

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


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

  389.     /** Build a propagator from orbit, attitude provider, mass and potential provider.
  390.      * <p>Using this constructor, it is possible to define the initial orbit as
  391.      * a mean Brouwer-Lyddane orbit or an osculating one.</p>
  392.      * @param initialOrbit initial orbit
  393.      * @param attitudeProv attitude provider
  394.      * @param mass spacecraft mass
  395.      * @param provider for un-normalized zonal coefficients
  396.      * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
  397.      * @param M2 value of empirical drag coefficient in rad/s².
  398.      *        If equal to {@link #M2} drag is not computed
  399.      */
  400.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  401.                                     final AttitudeProvider attitudeProv,
  402.                                     final double mass,
  403.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  404.                                     final PropagationType initialType, final double M2) {
  405.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), 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 #M2} drag is not computed
  420.      */
  421.     public BrouwerLyddanePropagator(final Orbit initialOrbit,
  422.                                     final AttitudeProvider attitude,
  423.                                     final double mass,
  424.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  425.                                     final UnnormalizedSphericalHarmonics harmonics,
  426.                                     final PropagationType initialType, final double M2) {
  427.         this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
  428.              harmonics.getUnnormalizedCnm(2, 0),
  429.              harmonics.getUnnormalizedCnm(3, 0),
  430.              harmonics.getUnnormalizedCnm(4, 0),
  431.              harmonics.getUnnormalizedCnm(5, 0),
  432.              initialType, M2);
  433.     }

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

  470.         super(attitudeProv);

  472.         // store model coefficients
  473.         this.referenceRadius = referenceRadius;
  474.         this.mu  = mu;
  475.         this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50};

  477.         // initialize M2 driver
  478.         this.M2Driver = new ParameterDriver(M2_NAME, M2, SCALE,
  479.                                             Double.NEGATIVE_INFINITY,
  480.                                             Double.POSITIVE_INFINITY);

  482.         // compute mean parameters if needed
  483.         resetInitialState(new SpacecraftState(initialOrbit,
  484.                                               attitudeProv.getAttitude(initialOrbit,
  485.                                                                        initialOrbit.getDate(),
  486.                                                                        initialOrbit.getFrame()),
  487.                                               mass),
  488.                                               initialType);

  490.     }

  492.     /** {@inheritDoc}
  493.      * <p>The new initial state to consider
  494.      * must be defined with an osculating orbit.</p>
  495.      * @see #resetInitialState(SpacecraftState, PropagationType)
  496.      */
  497.     public void resetInitialState(final SpacecraftState state) {
  498.         resetInitialState(state, PropagationType.OSCULATING);
  499.     }

  501.     /** Reset the propagator initial state.
  502.      * @param state new initial state to consider
  503.      * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
  504.      */
  505.     public void resetInitialState(final SpacecraftState state, final PropagationType stateType) {
  506.         super.resetInitialState(state);
  507.         final KeplerianOrbit keplerian = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state.getOrbit());
  508.         this.initialModel = (stateType == PropagationType.MEAN) ?
  509.                              new BLModel(keplerian, state.getMass(), referenceRadius, mu, ck0) :
  510.                              computeMeanParameters(keplerian, state.getMass());
  511.         this.models = new TimeSpanMap<BLModel>(initialModel);
  512.     }

  514.     /** {@inheritDoc} */
  515.     protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
  516.         final BLModel newModel = computeMeanParameters((KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state.getOrbit()),
  517.                                                        state.getMass());
  518.         if (forward) {
  519.             models.addValidAfter(newModel, state.getDate(), false);
  520.         } else {
  521.             models.addValidBefore(newModel, state.getDate(), false);
  522.         }
  523.         stateChanged(state);
  524.     }

  526.     /** Compute mean parameters according to the Brouwer-Lyddane analytical model computation
  527.      * in order to do the propagation.
  528.      * @param osculating osculating orbit
  529.      * @param mass constant mass
  530.      * @return Brouwer-Lyddane mean model
  531.      */
  532.     private BLModel computeMeanParameters(final KeplerianOrbit osculating, final double mass) {

  534.         // sanity check
  535.         if (osculating.getA() < referenceRadius) {
  536.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
  537.                                            osculating.getA());
  538.         }

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

  543.         // threshold for each parameter
  544.         final double epsilon         = 1.0e-13;
  545.         final double thresholdA      = epsilon * (1 + FastMath.abs(current.mean.getA()));
  546.         final double thresholdE      = epsilon * (1 + current.mean.getE());
  547.         final double thresholdAngles = epsilon * FastMath.PI;

  549.         int i = 0;
  550.         while (i++ < 200) {

  552.             // recompute the osculating parameters from the current mean parameters
  553.             final UnivariateDerivative2[] parameters = current.propagateParameters(current.mean.getDate());

  555.             // adapted parameters residuals
  556.             final double deltaA     = osculating.getA() - parameters[0].getValue();
  557.             final double deltaE     = osculating.getE() - parameters[1].getValue();
  558.             final double deltaI     = osculating.getI() - parameters[2].getValue();
  559.             final double deltaOmega = MathUtils.normalizeAngle(osculating.getPerigeeArgument() -
  560.                                       parameters[3].getValue(),
  561.                                       0.0);
  562.             final double deltaRAAN  = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode() -
  563.                                       parameters[4].getValue(),
  564.                                       0.0);
  565.             final double deltaAnom = MathUtils.normalizeAngle(osculating.getMeanAnomaly() -
  566.                                      parameters[5].getValue(),
  567.                                      0.0);


  570.             // update mean parameters
  571.             current = new BLModel(new KeplerianOrbit(current.mean.getA() + deltaA,
  572.                                                      current.mean.getE() + deltaE,
  573.                                                      current.mean.getI() + deltaI,
  574.                                                      current.mean.getPerigeeArgument() + deltaOmega,
  575.                                                      current.mean.getRightAscensionOfAscendingNode() + deltaRAAN,
  576.                                                      current.mean.getMeanAnomaly() + deltaAnom,
  577.                                                      PositionAngle.MEAN,
  578.                                                      current.mean.getFrame(),
  579.                                                      current.mean.getDate(), mu),
  580.                                   mass, referenceRadius, mu, ck0);
  581.             // check convergence
  582.             if (FastMath.abs(deltaA)     < thresholdA &&
  583.                 FastMath.abs(deltaE)     < thresholdE &&
  584.                 FastMath.abs(deltaI)     < thresholdAngles &&
  585.                 FastMath.abs(deltaOmega) < thresholdAngles &&
  586.                 FastMath.abs(deltaRAAN)  < thresholdAngles &&
  587.                 FastMath.abs(deltaAnom)  < thresholdAngles) {
  588.                 return current;
  589.             }
  590.         }
  591.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_BROUWER_LYDDANE_MEAN_PARAMETERS, i);
  592.     }


  595.     /** {@inheritDoc} */
  596.     public KeplerianOrbit propagateOrbit(final AbsoluteDate date) {
  597.         // compute Cartesian parameters, taking derivatives into account
  598.         // to make sure velocity and acceleration are consistent
  599.         final BLModel current = models.get(date);
  600.         final UnivariateDerivative2[] propOrb_parameters = current.propagateParameters(date);
  601.         return new KeplerianOrbit(propOrb_parameters[0].getValue(), propOrb_parameters[1].getValue(),
  602.                                   propOrb_parameters[2].getValue(), propOrb_parameters[3].getValue(),
  603.                                   propOrb_parameters[4].getValue(), propOrb_parameters[5].getValue(),
  604.                                   PositionAngle.MEAN, current.mean.getFrame(), date, mu);
  605.     }

  607.     /**
  608.      * Get the value of the M2 drag parameter.
  609.      * @return the value of the M2 drag parameter
  610.      */
  611.     public double getM2() {
  612.         return M2Driver.getValue();
  613.     }

  615.     /**
  616.      * Get the central attraction coefficient μ.
  617.      * @return mu central attraction coefficient (m³/s²)
  618.      */
  619.     public double getMu() {
  620.         return mu;
  621.     }

  623.     /**
  624.      * Get the un-normalized zonal coefficients.
  625.      * @return the un-normalized zonal coefficients
  626.      */
  627.     public double[] getCk0() {
  628.         return ck0.clone();
  629.     }

  631.     /**
  632.      * Get the reference radius of the central body attraction model.
  633.      * @return the reference radius in meters
  634.      */
  635.     public double getReferenceRadius() {
  636.         return referenceRadius;
  637.     }

  639.     /**
  640.      * Get the parameters driver for propagation model.
  641.      * @return drivers for propagation model
  642.      */
  643.     public List<ParameterDriver> getParametersDrivers() {
  644.         return Collections.singletonList(M2Driver);
  645.     }

  647.     /** {@inheritDoc} */
  648.     @Override
  649.     protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
  650.                                                         final DoubleArrayDictionary initialJacobianColumns) {
  651.         return new BrouwerLyddaneHarvester(this, stmName, initialStm, initialJacobianColumns);
  652.     }

  654.     /**
  655.      * Get the names of the parameters in the matrix returned by {@link MatricesHarvester#getParametersJacobian}.
  656.      * @return names of the parameters (i.e. columns) of the Jacobian matrix
  657.      */
  658.     protected List<String> getJacobiansColumnsNames() {
  659.         final List<String> columnsNames = new ArrayList<>();
  660.         for (final ParameterDriver driver : getParametersDrivers()) {
  661.             if (driver.isSelected() && !columnsNames.contains(driver.getName())) {
  662.                 columnsNames.add(driver.getName());
  663.             }
  664.         }
  665.         Collections.sort(columnsNames);
  666.         return columnsNames;
  667.     }

  669.     /** Local class for Brouwer-Lyddane model. */
  670.     private class BLModel {

  672.         /** Mean orbit. */
  673.         private final KeplerianOrbit mean;

  675.         /** Constant mass. */
  676.         private final double mass;

  678.         // CHECKSTYLE: stop JavadocVariable check

  680.         // preprocessed values

  682.         // Constant for secular terms l'', g'', h''
  683.         // l standing for true anomaly, g for perigee argument and h for raan
  684.         private final double xnotDot;
  685.         private final double n;
  686.         private final double lt;
  687.         private final double gt;
  688.         private final double ht;


  691.         // Long period terms
  692.         private final double dei3sg;
  693.         private final double de2sg;
  694.         private final double deisg;
  695.         private final double de;


  698.         private final double dlgs2g;
  699.         private final double dlgc3g;
  700.         private final double dlgcg;


  703.         private final double dh2sgcg;
  704.         private final double dhsgcg;
  705.         private final double dhcg;


  708.         // Short period terms
  709.         private final double aC;
  710.         private final double aCbis;
  711.         private final double ac2g2f;


  714.         private final double eC;
  715.         private final double ecf;
  716.         private final double e2cf;
  717.         private final double e3cf;
  718.         private final double ec2f2g;
  719.         private final double ecfc2f2g;
  720.         private final double e2cfc2f2g;
  721.         private final double e3cfc2f2g;
  722.         private final double ec2gf;
  723.         private final double ec2g3f;


  726.         private final double ide;
  727.         private final double isfs2f2g;
  728.         private final double icfc2f2g;
  729.         private final double ic2f2g;


  732.         private final double glf;
  733.         private final double gll;
  734.         private final double glsf;
  735.         private final double glosf;
  736.         private final double gls2f2g;
  737.         private final double gls2gf;
  738.         private final double glos2gf;
  739.         private final double gls2g3f;
  740.         private final double glos2g3f;


  743.         private final double hf;
  744.         private final double hl;
  745.         private final double hsf;
  746.         private final double hcfs2g2f;
  747.         private final double hs2g2f;
  748.         private final double hsfc2g2f;


  751.         private final double edls2g;
  752.         private final double edlcg;
  753.         private final double edlc3g;
  754.         private final double edlsf;
  755.         private final double edls2gf;
  756.         private final double edls2g3f;

  758.         // Drag terms
  759.         private final double aRate;
  760.         private final double eRate;

  762.         // CHECKSTYLE: resume JavadocVariable check

  764.         /** Create a model for specified mean orbit.
  765.          * @param mean mean orbit
  766.          * @param mass constant mass
  767.          * @param referenceRadius reference radius of the central body attraction model (m)
  768.          * @param mu central attraction coefficient (m³/s²)
  769.          * @param ck0 un-normalized zonal coefficients
  770.          */
  771.         BLModel(final KeplerianOrbit mean, final double mass,
  772.                 final double referenceRadius, final double mu, final double[] ck0) {

  774.             this.mean = mean;
  775.             this.mass = mass;

  777.             final double app = mean.getA();
  778.             xnotDot = FastMath.sqrt(mu / app) / app;

  780.             // preliminary processing
  781.             final double q = referenceRadius / app;
  782.             double ql = q * q;
  783.             final double y2 = -0.5 * ck0[2] * ql;

  785.             n = FastMath.sqrt(1 - mean.getE() * mean.getE());
  786.             final double n2 = n * n;
  787.             final double n3 = n2 * n;
  788.             final double n4 = n2 * n2;
  789.             final double n6 = n4 * n2;
  790.             final double n8 = n4 * n4;
  791.             final double n10 = n8 * n2;

  793.             final double yp2 = y2 / n4;
  794.             ql *= q;
  795.             final double yp3 = ck0[3] * ql / n6;
  796.             ql *= q;
  797.             final double yp4 = 0.375 * ck0[4] * ql / n8;
  798.             ql *= q;
  799.             final double yp5 = ck0[5] * ql / n10;


  802.             final SinCos sc    = FastMath.sinCos(mean.getI());
  803.             final double sinI1 = sc.sin();
  804.             final double sinI2 = sinI1 * sinI1;
  805.             final double cosI1 = sc.cos();
  806.             final double cosI2 = cosI1 * cosI1;
  807.             final double cosI3 = cosI2 * cosI1;
  808.             final double cosI4 = cosI2 * cosI2;
  809.             final double cosI6 = cosI4 * cosI2;
  810.             final double C5c2 = 1.0 / T2(cosI1);
  811.             final double C3c2 = 3.0 * cosI2 - 1.0;

  813.             final double epp = mean.getE();
  814.             final double epp2 = epp * epp;
  815.             final double epp3 = epp2 * epp;
  816.             final double epp4 = epp2 * epp2;

  818.             if (epp >= 1) {
  819.                 // Only for elliptical (e < 1) orbits
  820.                 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
  821.                                           mean.getE());
  822.             }

  824.             // secular multiplicative
  825.             lt = 1 +
  826.                     1.5 * yp2 * n * C3c2 +
  827.                     0.09375 * yp2 * yp2 * n * (-15.0 + 16.0 * n + 25.0 * n2 + (30.0 - 96.0 * n - 90.0 * n2) * cosI2 + (105.0 + 144.0 * n + 25.0 * n2) * cosI4) +
  828.                     0.9375 * yp4 * n * epp2 * (3.0 - 30.0 * cosI2 + 35.0 * cosI4);

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

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

  838.             final double cA = 1.0 - 11.0 * cosI2 - 40.0 * cosI4 / C5c2;
  839.             final double cB = 1.0 - 3.0 * cosI2 - 8.0 * cosI4 / C5c2;
  840.             final double cC = 1.0 - 9.0 * cosI2 - 24.0 * cosI4 / C5c2;
  841.             final double cD = 1.0 - 5.0 * cosI2 - 16.0 * cosI4 / C5c2;

  843.             final double qyp2_4 = 3.0 * yp2 * yp2 * cA -
  844.                                   10.0 * yp4 * cB;
  845.             final double qyp52 = epp3 * cosI1 * (0.5 * cD / sinI1 +
  846.                                  sinI1 * (5.0 + 32.0 * cosI2 / C5c2 + 80.0 * cosI4 / C5c2 / C5c2));
  847.             final double qyp22 = 2.0 + epp2 - 11.0 * (2.0 + 3.0 * epp2) * cosI2 -
  848.                                  40.0 * (2.0 + 5.0 * epp2) * cosI4 / C5c2 -
  849.                                  400.0 * epp2 * cosI6 / C5c2 / C5c2;
  850.             final double qyp42 = ( qyp22 + 4.0 * (2.0 + epp2 - (2.0 + 3.0 * epp2) * cosI2) ) / 5.0;
  851.             final double qyp52bis = epp * cosI1 * sinI1 *
  852.                                     (4.0 + 3.0 * epp2) *
  853.                                     (3.0 + 16.0 * cosI2 / C5c2 + 40.0 * cosI4 / C5c2 / C5c2);

  855.             // long periodic multiplicative
  856.             dei3sg =  35.0 / 96.0 * yp5 / yp2 * epp2 * n2 * cD * sinI1;
  857.             de2sg = -1.0 / 12.0 * epp * n2 / yp2 * qyp2_4;
  858.             deisg = ( -35.0 / 128.0 * yp5 / yp2 * epp2 * n2 * cD +
  859.                     1.0 / 4.0 * n2 / yp2 * (yp3 + 5.0 / 16.0 * yp5 * (4.0 + 3.0 * epp2) * cC)) * sinI1;
  860.             de = epp2 * n2 / 24.0 / yp2 * qyp2_4;

  862.             final double qyp52quotient = epp * (-32.0 + 81.0 * epp4) / (4.0 + 3.0 * epp2 + n * (4.0 + 9.0 * epp2));
  863.             dlgs2g = 1.0 / 48.0 / yp2 * (-3.0 * yp2 * yp2 * qyp22 +
  864.                      10.0 * yp4 * qyp42 ) +
  865.                      n3 / yp2 * qyp2_4 / 24.0;
  866.             dlgc3g = 35.0 / 384.0 * yp5 / yp2 * n3 * epp * cD * sinI1 +
  867.                      35.0 / 1152.0 * yp5 / yp2 * (2.0 * qyp52 * cosI1 - epp * cD * sinI1 * (3.0 + 2.0 * epp2));
  868.             dlgcg = -yp3 * epp * cosI2 / ( 4.0 * yp2 * sinI1) +
  869.                     0.078125 * yp5 / yp2 * (-epp * cosI2 / sinI1 * (4.0 + 3.0 * epp2) + epp2 * sinI1 * (26.0 + 9.0 * epp2)) * cC -
  870.                     0.46875 * yp5 / yp2 * qyp52bis * cosI1 +
  871.                     0.25 * yp3 / yp2 * sinI1 * epp / (1.0 + n3) * (3.0 - epp2 * (3.0 - epp2)) +
  872.                     0.078125 * yp5 / yp2 * n2 * cC * qyp52quotient * sinI1;


  875.             final double qyp24 = 3.0 * yp2 * yp2 * (11.0 + 80.0 * cosI2 / sinI1 + 200.0 * cosI4 / sinI2) -
  876.                                  10.0 * yp4 * (3.0 + 16.0 * cosI2 / sinI1 + 40.0 * cosI4 / sinI2);
  877.             dh2sgcg = 35.0 / 144.0 * yp5 / yp2 * qyp52;
  878.             dhsgcg = -epp2 * cosI1 / (12.0 * yp2) * qyp24;
  879.             dhcg = -35.0 / 576.0 * yp5 / yp2 * qyp52 +
  880.                    epp * cosI1 / (4.0 * yp2 * sinI1) * (yp3 + 0.3125 * yp5 * (4.0 + 3.0 * epp2) * cC) +
  881.                    1.875 / (4.0 * yp2) * yp5 * qyp52bis;

  883.             // short periodic multiplicative
  884.             aC = -yp2 * C3c2 * app / n3;
  885.             aCbis = y2 * app * C3c2;
  886.             ac2g2f = y2 * app * 3.0 * sinI2;

  888.             double qe = 0.5 * n2 * y2 * C3c2 / n6;
  889.             eC = qe * epp / (1.0 + n3) * (3.0 - epp2 * (3.0 - epp2));
  890.             ecf = 3.0 * qe;
  891.             e2cf = 3.0 * epp * qe;
  892.             e3cf = epp2 * qe;
  893.             qe = 0.5 * n2 * y2 * 3.0 * (1.0 - cosI2) / n6;
  894.             ec2f2g = qe * epp;
  895.             ecfc2f2g = 3.0 * qe;
  896.             e2cfc2f2g = 3.0 * epp * qe;
  897.             e3cfc2f2g = epp2 * qe;
  898.             qe = -0.5 * yp2 * n2 * (1.0 - cosI2);
  899.             ec2gf = 3.0 * qe;
  900.             ec2g3f = qe;

  902.             double qi = epp * yp2 * cosI1 * sinI1;
  903.             ide = -epp * cosI1 / (n2 * sinI1);
  904.             isfs2f2g = qi;
  905.             icfc2f2g = 2.0 * qi;
  906.             qi = yp2 * cosI1 * sinI1;
  907.             ic2f2g = 1.5 * qi;

  909.             double qgl1 = 0.25 * yp2;
  910.             double qgl2 = 0.25 * yp2 * epp * n2 / (1.0 + n);
  911.             glf = qgl1 * -6.0 * C5c2;
  912.             gll = qgl1 * 6.0 * C5c2;
  913.             glsf = qgl1 * -6.0 * C5c2 * epp +
  914.                    qgl2 * 2.0 * C3c2;
  915.             glosf = qgl2 * 2.0 * C3c2;
  916.             qgl1 = qgl1 * (3.0 - 5.0 * cosI2);
  917.             qgl2 = qgl2 * 3.0 * (1.0 - cosI2);
  918.             gls2f2g = 3.0 * qgl1;
  919.             gls2gf = 3.0 * epp * qgl1 +
  920.                      qgl2;
  921.             glos2gf = -1.0 * qgl2;
  922.             gls2g3f = qgl1 * epp +
  923.                       1.0 / 3.0 * qgl2;
  924.             glos2g3f = qgl2;

  926.             final double qh = 3.0 * yp2 * cosI1;
  927.             hf = -qh;
  928.             hl = qh;
  929.             hsf = -epp * qh;
  930.             hcfs2g2f = 2.0 * epp * yp2 * cosI1;
  931.             hs2g2f = 1.5 * yp2 * cosI1;
  932.             hsfc2g2f = -epp * yp2 * cosI1;

  934.             final double qedl = -0.25 * yp2 * n3;
  935.             edls2g = 1.0 / 24.0 * epp * n3 / yp2 * qyp2_4;
  936.             edlcg = -0.25 * yp3 / yp2 * n3 * sinI1 -
  937.                     0.078125 * yp5 / yp2 * n3 * sinI1 * (4.0 + 9.0 * epp2) * cC;
  938.             edlc3g = 35.0 / 384.0 * yp5 / yp2 * n3 * epp2 * cD * sinI1;
  939.             edlsf = 2.0 * qedl * C3c2;
  940.             edls2gf = 3.0 * qedl * (1.0 - cosI2);
  941.             edls2g3f = 1.0 / 3.0 * qedl;

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

  948.         }

  950.         /**
  951.          * Accurate computation of E - e sin(E).
  952.          *
  953.          * @param E eccentric anomaly
  954.          * @return E - e sin(E)
  955.          */
  956.         private UnivariateDerivative2 eMeSinE(final UnivariateDerivative2 E) {
  957.             UnivariateDerivative2 x = E.sin().multiply(1 - mean.getE());
  958.             final UnivariateDerivative2 mE2 = E.negate().multiply(E);
  959.             UnivariateDerivative2 term = E;
  960.             UnivariateDerivative2 d    = E.getField().getZero();
  961.             // the inequality test below IS intentional and should NOT be replaced by a check with a small tolerance
  962.             for (UnivariateDerivative2 x0 = d.add(Double.NaN); !Double.valueOf(x.getValue()).equals(Double.valueOf(x0.getValue()));) {
  963.                 d = d.add(2);
  964.                 term = term.multiply(mE2.divide(d.multiply(d.add(1))));
  965.                 x0 = x;
  966.                 x = x.subtract(term);
  967.             }
  968.             return x;
  969.         }

  971.         /**
  972.          * Gets eccentric anomaly from mean anomaly.
  973.          * <p>The algorithm used to solve the Kepler equation has been published in:
  974.          * "Procedures for  solving Kepler's Equation", A. W. Odell and R. H. Gooding,
  975.          * Celestial Mechanics 38 (1986) 307-334</p>
  976.          * <p>It has been copied from the OREKIT library (KeplerianOrbit class).</p>
  977.          *
  978.          * @param mk the mean anomaly (rad)
  979.          * @return the eccentric anomaly (rad)
  980.          */
  981.         private UnivariateDerivative2 getEccentricAnomaly(final UnivariateDerivative2 mk) {
  982.             final double k1 = 3 * FastMath.PI + 2;
  983.             final double k2 = FastMath.PI - 1;
  984.             final double k3 = 6 * FastMath.PI - 1;
  985.             final double A  = 3.0 * k2 * k2 / k1;
  986.             final double B  = k3 * k3 / (6.0 * k1);
  987.             // reduce M to [-PI PI] interval
  988.             final UnivariateDerivative2 reducedM = new UnivariateDerivative2(MathUtils.normalizeAngle(mk.getValue(), 0.0),
  989.                                                                              mk.getFirstDerivative(),
  990.                                                                              mk.getSecondDerivative());

  992.             // compute start value according to A. W. Odell and R. H. Gooding S12 starter
  993.             UnivariateDerivative2 ek;
  994.             if (FastMath.abs(reducedM.getValue()) < 1.0 / 6.0) {
  995.                 if (FastMath.abs(reducedM.getValue()) < Precision.SAFE_MIN) {
  996.                     // this is an Orekit change to the S12 starter.
  997.                     // If reducedM is 0.0, the derivative of cbrt is infinite which induces NaN appearing later in
  998.                     // the computation. As in this case E and M are almost equal, we initialize ek with reducedM
  999.                     ek = reducedM;
  1000.                 } else {
  1001.                     // this is the standard S12 starter
  1002.                     ek = reducedM.add(reducedM.multiply(6).cbrt().subtract(reducedM).multiply(mean.getE()));
  1003.                 }
  1004.             } else {
  1005.                 if (reducedM.getValue() < 0) {
  1006.                     final UnivariateDerivative2 w = reducedM.add(FastMath.PI);
  1007.                     ek = reducedM.add(w.multiply(-A).divide(w.subtract(B)).subtract(FastMath.PI).subtract(reducedM).multiply(mean.getE()));
  1008.                 } else {
  1009.                     final UnivariateDerivative2 minusW = reducedM.subtract(FastMath.PI);
  1010.                     ek = reducedM.add(minusW.multiply(A).divide(minusW.add(B)).add(FastMath.PI).subtract(reducedM).multiply(mean.getE()));
  1011.                 }
  1012.             }

  1014.             final double e1 = 1 - mean.getE();
  1015.             final boolean noCancellationRisk = (e1 + ek.getValue() * ek.getValue() / 6) >= 0.1;

  1017.             // perform two iterations, each consisting of one Halley step and one Newton-Raphson step
  1018.             for (int j = 0; j < 2; ++j) {
  1019.                 final UnivariateDerivative2 f;
  1020.                 UnivariateDerivative2 fd;
  1021.                 final UnivariateDerivative2 fdd  = ek.sin().multiply(mean.getE());
  1022.                 final UnivariateDerivative2 fddd = ek.cos().multiply(mean.getE());
  1023.                 if (noCancellationRisk) {
  1024.                     f  = ek.subtract(fdd).subtract(reducedM);
  1025.                     fd = fddd.subtract(1).negate();
  1026.                 } else {
  1027.                     f  = eMeSinE(ek).subtract(reducedM);
  1028.                     final UnivariateDerivative2 s = ek.multiply(0.5).sin();
  1029.                     fd = s.multiply(s).multiply(2 * mean.getE()).add(e1);
  1030.                 }
  1031.                 final UnivariateDerivative2 dee = f.multiply(fd).divide(f.multiply(0.5).multiply(fdd).subtract(fd.multiply(fd)));

  1033.                 // update eccentric anomaly, using expressions that limit underflow problems
  1034.                 final UnivariateDerivative2 w = fd.add(dee.multiply(0.5).multiply(fdd.add(dee.multiply(fdd).divide(3))));
  1035.                 fd = fd.add(dee.multiply(fdd.add(dee.multiply(0.5).multiply(fdd))));
  1036.                 ek = ek.subtract(f.subtract(dee.multiply(fd.subtract(w))).divide(fd));
  1037.             }

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

  1042.             // Returns the eccentric anomaly
  1043.             return ek;
  1044.         }

  1046.         /**
  1047.          * This method is used in Brouwer-Lyddane model to avoid singularity at the
  1048.          * critical inclination (i = 63.4°).
  1049.          * <p>
  1050.          * This method, based on Warren Phipps's 1992 thesis (Eq. 2.47 and 2.48),
  1051.          * approximate the factor (1.0 - 5.0 * cos²(inc))^-1 (causing the singularity)
  1052.          * by a function, named T2 in the thesis.
  1053.          * </p>
  1054.          * @param cosInc cosine of the mean inclination
  1055.          * @return an approximation of (1.0 - 5.0 * cos²(inc))^-1 term
  1056.          */
  1057.         private double T2(final double cosInc) {

  1059.             // X = (1.0 - 5.0 * cos²(inc))
  1060.             final double x  = 1.0 - 5.0 * cosInc * cosInc;
  1061.             final double x2 = x * x;

  1063.             // Eq. 2.48
  1064.             double sum = 0.0;
  1065.             for (int i = 0; i <= 12; i++) {
  1066.                 final double sign = i % 2 == 0 ? +1.0 : -1.0;
  1067.                 sum += sign * FastMath.pow(BETA, i) * FastMath.pow(x2, i) / CombinatoricsUtils.factorialDouble(i + 1);
  1068.             }

  1070.             // Right term of equation 2.47
  1071.             double product = 1.0;
  1072.             for (int i = 0; i <= 10; i++) {
  1073.                 product *= 1 + FastMath.exp(FastMath.scalb(-1.0, i) * BETA * x2);
  1074.             }

  1076.             // Return (Eq. 2.47)
  1077.             return BETA * x * sum * product;

  1079.         }

  1081.         /** Extrapolate an orbit up to a specific target date.
  1082.          * @param date target date for the orbit
  1083.          * @return propagated parameters
  1084.          */
  1085.         public UnivariateDerivative2[] propagateParameters(final AbsoluteDate date) {

  1087.             // Empirical drag coefficient M2
  1088.             final double m2 = M2Driver.getValue();

  1090.             // Keplerian evolution
  1091.             final UnivariateDerivative2 dt = new UnivariateDerivative2(date.durationFrom(mean.getDate()), 1.0, 0.0);
  1092.             final UnivariateDerivative2 xnot = dt.multiply(xnotDot);

  1094.             //____________________________________
  1095.             // secular effects

  1097.             // mean mean anomaly (with drag Eq. 2.38 of Phipps' 1992 thesis)
  1098.             final UnivariateDerivative2 dtM2  = dt.multiply(m2);
  1099.             final UnivariateDerivative2 dt2M2 = dt.multiply(dtM2);
  1100.             final UnivariateDerivative2 lpp = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getMeanAnomaly() + lt * xnot.getValue() + dt2M2.getValue(),
  1101.                                                                                                FastMath.PI),
  1102.                                                                       lt * xnotDot + 2.0 * dtM2.getValue(),
  1103.                                                                       2.0 * m2);
  1104.             // mean argument of perigee
  1105.             final UnivariateDerivative2 gpp = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getPerigeeArgument() + gt * xnot.getValue(),
  1106.                                                                                                FastMath.PI),
  1107.                                                                       gt * xnotDot,
  1108.                                                                       0.0);
  1109.             // mean longitude of ascending node
  1110.             final UnivariateDerivative2 hpp = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode() + ht * xnot.getValue(),
  1111.                                                                                                FastMath.PI),
  1112.                                                                       ht * xnotDot,
  1113.                                                                       0.0);

  1115.             // ________________________________________________
  1116.             // secular rates of the mean semi-major axis and eccentricity

  1118.             // semi-major axis
  1119.             final UnivariateDerivative2 appDrag = dt.multiply(aRate * m2);

  1121.             // eccentricity
  1122.             final UnivariateDerivative2 eppDrag = dt.multiply(eRate * m2);

  1124.             //____________________________________
  1125.             // Long periodical terms
  1126.             final UnivariateDerivative2 cg1 = gpp.cos();
  1127.             final UnivariateDerivative2 sg1 = gpp.sin();
  1128.             final UnivariateDerivative2 c2g = cg1.multiply(cg1).subtract(sg1.multiply(sg1));
  1129.             final UnivariateDerivative2 s2g = cg1.multiply(sg1).add(sg1.multiply(cg1));
  1130.             final UnivariateDerivative2 c3g = c2g.multiply(cg1).subtract(s2g.multiply(sg1));
  1131.             final UnivariateDerivative2 sg2 = sg1.multiply(sg1);
  1132.             final UnivariateDerivative2 sg3 = sg1.multiply(sg2);


  1135.             // de eccentricity
  1136.             final UnivariateDerivative2 d1e = sg3.multiply(dei3sg).
  1137.                                               add(sg1.multiply(deisg)).
  1138.                                               add(sg2.multiply(de2sg)).
  1139.                                               add(de);

  1141.             // l' + g'
  1142.             final UnivariateDerivative2 lp_p_gp = s2g.multiply(dlgs2g).
  1143.                                                add(c3g.multiply(dlgc3g)).
  1144.                                                add(cg1.multiply(dlgcg)).
  1145.                                                add(lpp).
  1146.                                                add(gpp);

  1148.             // h'
  1149.             final UnivariateDerivative2 hp = sg2.multiply(cg1).multiply(dh2sgcg).
  1150.                                                add(sg1.multiply(cg1).multiply(dhsgcg)).
  1151.                                                add(cg1.multiply(dhcg)).
  1152.                                                add(hpp);

  1154.             // Short periodical terms
  1155.             // eccentric anomaly
  1156.             final UnivariateDerivative2 Ep = getEccentricAnomaly(lpp);
  1157.             final UnivariateDerivative2 cf1 = (Ep.cos().subtract(mean.getE())).divide(Ep.cos().multiply(-mean.getE()).add(1.0));
  1158.             final UnivariateDerivative2 sf1 = (Ep.sin().multiply(n)).divide(Ep.cos().multiply(-mean.getE()).add(1.0));
  1159.             final UnivariateDerivative2 f = FastMath.atan2(sf1, cf1);

  1161.             final UnivariateDerivative2 c2f = cf1.multiply(cf1).subtract(sf1.multiply(sf1));
  1162.             final UnivariateDerivative2 s2f = cf1.multiply(sf1).add(sf1.multiply(cf1));
  1163.             final UnivariateDerivative2 c3f = c2f.multiply(cf1).subtract(s2f.multiply(sf1));
  1164.             final UnivariateDerivative2 s3f = c2f.multiply(sf1).add(s2f.multiply(cf1));
  1165.             final UnivariateDerivative2 cf2 = cf1.multiply(cf1);
  1166.             final UnivariateDerivative2 cf3 = cf1.multiply(cf2);

  1168.             final UnivariateDerivative2 c2g1f = cf1.multiply(c2g).subtract(sf1.multiply(s2g));
  1169.             final UnivariateDerivative2 c2g2f = c2f.multiply(c2g).subtract(s2f.multiply(s2g));
  1170.             final UnivariateDerivative2 c2g3f = c3f.multiply(c2g).subtract(s3f.multiply(s2g));
  1171.             final UnivariateDerivative2 s2g1f = cf1.multiply(s2g).add(c2g.multiply(sf1));
  1172.             final UnivariateDerivative2 s2g2f = c2f.multiply(s2g).add(c2g.multiply(s2f));
  1173.             final UnivariateDerivative2 s2g3f = c3f.multiply(s2g).add(c2g.multiply(s3f));

  1175.             final UnivariateDerivative2 eE = (Ep.cos().multiply(-mean.getE()).add(1.0)).reciprocal();
  1176.             final UnivariateDerivative2 eE3 = eE.multiply(eE).multiply(eE);
  1177.             final UnivariateDerivative2 sigma = eE.multiply(n * n).multiply(eE).add(eE);

  1179.             // Semi-major axis (with drag Eq. 2.41 of Phipps' 1992 thesis)
  1180.             final UnivariateDerivative2 a = eE3.multiply(aCbis).add(appDrag.add(mean.getA())).
  1181.                                             add(aC).
  1182.                                             add(eE3.multiply(c2g2f).multiply(ac2g2f));

  1184.             // Eccentricity (with drag Eq. 2.45 of Phipps' 1992 thesis)
  1185.             final UnivariateDerivative2 e = d1e.add(eppDrag.add(mean.getE())).
  1186.                                             add(eC).
  1187.                                             add(cf1.multiply(ecf)).
  1188.                                             add(cf2.multiply(e2cf)).
  1189.                                             add(cf3.multiply(e3cf)).
  1190.                                             add(c2g2f.multiply(ec2f2g)).
  1191.                                             add(c2g2f.multiply(cf1).multiply(ecfc2f2g)).
  1192.                                             add(c2g2f.multiply(cf2).multiply(e2cfc2f2g)).
  1193.                                             add(c2g2f.multiply(cf3).multiply(e3cfc2f2g)).
  1194.                                             add(c2g1f.multiply(ec2gf)).
  1195.                                             add(c2g3f.multiply(ec2g3f));

  1197.             // Inclination
  1198.             final UnivariateDerivative2 i = d1e.multiply(ide).
  1199.                                             add(mean.getI()).
  1200.                                             add(sf1.multiply(s2g2f.multiply(isfs2f2g))).
  1201.                                             add(cf1.multiply(c2g2f.multiply(icfc2f2g))).
  1202.                                             add(c2g2f.multiply(ic2f2g));

  1204.             // Argument of perigee + True anomaly
  1205.             final UnivariateDerivative2 g_p_l = lp_p_gp.add(f.multiply(glf)).
  1206.                                                 add(lpp.multiply(gll)).
  1207.                                                 add(sf1.multiply(glsf)).
  1208.                                                 add(sigma.multiply(sf1).multiply(glosf)).
  1209.                                                 add(s2g2f.multiply(gls2f2g)).
  1210.                                                 add(s2g1f.multiply(gls2gf)).
  1211.                                                 add(sigma.multiply(s2g1f).multiply(glos2gf)).
  1212.                                                 add(s2g3f.multiply(gls2g3f)).
  1213.                                                 add(sigma.multiply(s2g3f).multiply(glos2g3f));


  1216.             // Longitude of ascending node
  1217.             final UnivariateDerivative2 h = hp.add(f.multiply(hf)).
  1218.                                             add(lpp.multiply(hl)).
  1219.                                             add(sf1.multiply(hsf)).
  1220.                                             add(cf1.multiply(s2g2f).multiply(hcfs2g2f)).
  1221.                                             add(s2g2f.multiply(hs2g2f)).
  1222.                                             add(c2g2f.multiply(sf1).multiply(hsfc2g2f));

  1224.             final UnivariateDerivative2 edl = s2g.multiply(edls2g).
  1225.                                             add(cg1.multiply(edlcg)).
  1226.                                             add(c3g.multiply(edlc3g)).
  1227.                                             add(sf1.multiply(edlsf)).
  1228.                                             add(s2g1f.multiply(edls2gf)).
  1229.                                             add(s2g3f.multiply(edls2g3f)).
  1230.                                             add(sf1.multiply(sigma).multiply(edlsf)).
  1231.                                             add(s2g1f.multiply(sigma).multiply(-edls2gf)).
  1232.                                             add(s2g3f.multiply(sigma).multiply(3.0 * edls2g3f));

  1234.             final UnivariateDerivative2 A = e.multiply(lpp.cos()).subtract(edl.multiply(lpp.sin()));
  1235.             final UnivariateDerivative2 B = e.multiply(lpp.sin()).add(edl.multiply(lpp.cos()));

  1237.             // True anomaly
  1238.             final UnivariateDerivative2 l = FastMath.atan2(B, A);

  1240.             // Argument of perigee
  1241.             final UnivariateDerivative2 g = g_p_l.subtract(l);

  1243.             return new UnivariateDerivative2[] { a, e, i, g, h, l };

  1245.         }

  1247.     }

  1249.     /** {@inheritDoc} */
  1250.     protected double getMass(final AbsoluteDate date) {
  1251.         return models.get(date).mass;
  1252.     }

  1254. }
Created with JaCoCo 0.8.7.202105040129