EcksteinHechlerPropagator.java

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

  19. import org.hipparchus.analysis.differentiation.UnivariateDerivative2;
  20. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  21. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  22. import org.hipparchus.linear.RealMatrix;
  23. import org.hipparchus.util.FastMath;
  24. import org.hipparchus.util.MathUtils;
  25. import org.hipparchus.util.SinCos;
  26. import org.orekit.attitudes.AttitudeProvider;
  27. import org.orekit.attitudes.FrameAlignedProvider;
  28. import org.orekit.errors.OrekitException;
  29. import org.orekit.errors.OrekitMessages;
  30. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  31. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  32. import org.orekit.orbits.CartesianOrbit;
  33. import org.orekit.orbits.CircularOrbit;
  34. import org.orekit.orbits.Orbit;
  35. import org.orekit.orbits.OrbitType;
  36. import org.orekit.orbits.PositionAngleType;
  37. import org.orekit.propagation.AbstractMatricesHarvester;
  38. import org.orekit.propagation.PropagationType;
  39. import org.orekit.propagation.SpacecraftState;
  40. import org.orekit.propagation.conversion.osc2mean.EcksteinHechlerTheory;
  41. import org.orekit.propagation.conversion.osc2mean.FixedPointConverter;
  42. import org.orekit.propagation.conversion.osc2mean.MeanTheory;
  43. import org.orekit.propagation.conversion.osc2mean.OsculatingToMeanConverter;
  44. import org.orekit.time.AbsoluteDate;
  45. import org.orekit.utils.DoubleArrayDictionary;
  46. import org.orekit.utils.TimeSpanMap;
  47. import org.orekit.utils.TimeStampedPVCoordinates;

  49. /** This class propagates a {@link org.orekit.propagation.SpacecraftState}
  50.  *  using the analytical Eckstein-Hechler model.
  51.  * <p>The Eckstein-Hechler model is suited for near circular orbits
  52.  * (e &lt; 0.1, with poor accuracy between 0.005 and 0.1) and inclination
  53.  * neither equatorial (direct or retrograde) nor critical (direct or
  54.  * retrograde).</p>
  55.  * <p>
  56.  * Note that before version 7.0, there was a large inconsistency in the generated
  57.  * orbits, and it was fixed as of version 7.0 of Orekit, with a visible side effect.
  58.  * The problems is that if the circular parameters produced by the Eckstein-Hechler
  59.  * model are used to build an orbit considered to be osculating, the velocity deduced
  60.  * from this orbit was <em>inconsistent with the position evolution</em>! The reason is
  61.  * that the model includes non-Keplerian effects but it does not include a corresponding
  62.  * circular/Cartesian conversion. As a consequence, all subsequent computation involving
  63.  * velocity were wrong. This includes attitude modes like yaw compensation and Doppler
  64.  * effect. As this effect was considered serious enough and as accurate velocities were
  65.  * considered important, the propagator now generates {@link CartesianOrbit Cartesian
  66.  * orbits} which are built in a special way to ensure consistency throughout propagation.
  67.  * A side effect is that if circular parameters are rebuilt by user from these propagated
  68.  * Cartesian orbit, the circular parameters will generally <em>not</em> match the initial
  69.  * orbit (differences in semi-major axis can exceed 120 m). The position however <em>will</em>
  70.  * match to sub-micrometer level, and this position will be identical to the positions
  71.  * that were generated by previous versions (in other words, the internals of the models
  72.  * have not been changed, only the output parameters have been changed). The correctness
  73.  * of the initialization has been assessed and is good, as it allows the subsequent orbit
  74.  * to remain close to a numerical reference orbit.
  75.  * </p>
  76.  * <p>
  77.  * If users need a more definitive initialization of an Eckstein-Hechler propagator, they
  78.  * should consider using a {@link org.orekit.propagation.conversion.PropagatorConverter
  79.  * propagator converter} to initialize their Eckstein-Hechler propagator using a complete
  80.  * sample instead of just a single initial orbit.
  81.  * </p>
  82.  * @see Orbit
  83.  * @author Guylaine Prat
  84.  */
  85. public class EcksteinHechlerPropagator extends AbstractAnalyticalPropagator {

  87.     /** Default convergence threshold for mean parameters conversion. */
  88.     private static final double EPSILON_DEFAULT = 1.0e-11;

  90.     /** Default value for maxIterations. */
  91.     private static final int MAX_ITERATIONS_DEFAULT = 100;

  93.     /** Initial Eckstein-Hechler model. */
  94.     private EHModel initialModel;

  96.     /** All models. */
  97.     private TimeSpanMap<EHModel> models;

  99.     /** Reference radius of the central body attraction model (m). */
  100.     private final double referenceRadius;

  102.     /** Central attraction coefficient (m³/s²). */
  103.     private final double mu;

  105.     /** Un-normalized zonal coefficients. */
  106.     private final double[] ck0;

  108.     /** Build a propagator from orbit and potential provider.
  109.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  110.      *
  111.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  112.      *
  113.      * @param initialOrbit initial orbit
  114.      * @param provider for un-normalized zonal coefficients
  115.      * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider,
  116.      * UnnormalizedSphericalHarmonicsProvider)
  117.      * @see #EcksteinHechlerPropagator(Orbit, UnnormalizedSphericalHarmonicsProvider,
  118.      *                                 PropagationType)
  119.      */
  120.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  121.                                      final UnnormalizedSphericalHarmonicsProvider provider) {
  122.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  123.              DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()));
  124.     }

  126.     /**
  127.      * Private helper constructor.
  128.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  129.      * @param initialOrbit initial orbit
  130.      * @param attitude attitude provider
  131.      * @param mass spacecraft mass
  132.      * @param provider for un-normalized zonal coefficients
  133.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  134.      * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
  135.      *                                 UnnormalizedSphericalHarmonicsProvider,
  136.      *                                 UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics,
  137.      *                                 PropagationType)
  138.      */
  139.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  140.                                      final AttitudeProvider attitude,
  141.                                      final double mass,
  142.                                      final UnnormalizedSphericalHarmonicsProvider provider,
  143.                                      final UnnormalizedSphericalHarmonics harmonics) {
  144.         this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
  145.              harmonics.getUnnormalizedCnm(2, 0),
  146.              harmonics.getUnnormalizedCnm(3, 0),
  147.              harmonics.getUnnormalizedCnm(4, 0),
  148.              harmonics.getUnnormalizedCnm(5, 0),
  149.              harmonics.getUnnormalizedCnm(6, 0));
  150.     }

  152.     /** Build a propagator from orbit and potential.
  153.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  154.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  155.      * are related to both the normalized coefficients
  156.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  157.      *  and the J<sub>n</sub> one as follows:</p>
  158.      *
  159.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  160.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  161.      *
  162.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  163.      *
  164.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  165.      *
  166.      * @param initialOrbit initial orbit
  167.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  168.      * @param mu central attraction coefficient (m³/s²)
  169.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  170.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  171.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  172.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  173.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  174.      * @see org.orekit.utils.Constants
  175.      * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
  176.      * double, double, double, double, double)
  177.      */
  178.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  179.                                      final double referenceRadius,
  180.                                      final double mu,
  181.                                      final double c20,
  182.                                      final double c30,
  183.                                      final double c40,
  184.                                      final double c50,
  185.                                      final double c60) {
  186.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  187.              DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
  188.     }

  190.     /** Build a propagator from orbit, mass and potential provider.
  191.      * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
  192.      *
  193.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  194.      *
  195.      * @param initialOrbit initial orbit
  196.      * @param mass spacecraft mass
  197.      * @param provider for un-normalized zonal coefficients
  198.      * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
  199.      * UnnormalizedSphericalHarmonicsProvider)
  200.      */
  201.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  202.                                      final double mass,
  203.                                      final UnnormalizedSphericalHarmonicsProvider provider) {
  204.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  205.              mass, provider, provider.onDate(initialOrbit.getDate()));
  206.     }

  208.     /** Build a propagator from orbit, mass and potential.
  209.      * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
  210.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  211.      * are related to both the normalized coefficients
  212.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  213.      *  and the J<sub>n</sub> one as follows:</p>
  214.      *
  215.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  216.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  217.      *
  218.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  219.      *
  220.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  221.      *
  222.      * @param initialOrbit initial orbit
  223.      * @param mass spacecraft mass
  224.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  225.      * @param mu central attraction coefficient (m³/s²)
  226.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  227.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  228.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  229.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  230.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  231.      * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
  232.      * double, double, double, double, double)
  233.      */
  234.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  235.                                      final double mass,
  236.                                      final double referenceRadius,
  237.                                      final double mu,
  238.                                      final double c20,
  239.                                      final double c30,
  240.                                      final double c40,
  241.                                      final double c50,
  242.                                      final double c60) {
  243.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  244.              mass, referenceRadius, mu, c20, c30, c40, c50, c60);
  245.     }

  247.     /** Build a propagator from orbit, attitude provider and potential provider.
  248.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  249.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  250.      * @param initialOrbit initial orbit
  251.      * @param attitudeProv attitude provider
  252.      * @param provider for un-normalized zonal coefficients
  253.      */
  254.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  255.                                      final AttitudeProvider attitudeProv,
  256.                                      final UnnormalizedSphericalHarmonicsProvider provider) {
  257.         this(initialOrbit, attitudeProv, DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()));
  258.     }

  260.     /** Build a propagator from orbit, attitude provider and potential.
  261.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  262.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  263.      * are related to both the normalized coefficients
  264.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  265.      *  and the J<sub>n</sub> one as follows:</p>
  266.      *
  267.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  268.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  269.      *
  270.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  271.      *
  272.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  273.      *
  274.      * @param initialOrbit initial orbit
  275.      * @param attitudeProv attitude provider
  276.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  277.      * @param mu central attraction coefficient (m³/s²)
  278.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  279.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  280.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  281.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  282.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  283.      */
  284.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  285.                                      final AttitudeProvider attitudeProv,
  286.                                      final double referenceRadius,
  287.                                      final double mu,
  288.                                      final double c20,
  289.                                      final double c30,
  290.                                      final double c40,
  291.                                      final double c50,
  292.                                      final double c60) {
  293.         this(initialOrbit, attitudeProv, DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
  294.     }

  296.     /** Build a propagator from orbit, attitude provider, mass and potential provider.
  297.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  298.      * @param initialOrbit initial orbit
  299.      * @param attitudeProv attitude provider
  300.      * @param mass spacecraft mass
  301.      * @param provider for un-normalized zonal coefficients
  302.      * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
  303.      *                                 UnnormalizedSphericalHarmonicsProvider, PropagationType)
  304.      */
  305.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  306.                                      final AttitudeProvider attitudeProv,
  307.                                      final double mass,
  308.                                      final UnnormalizedSphericalHarmonicsProvider provider) {
  309.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()));
  310.     }

  312.     /** Build a propagator from orbit, attitude provider, mass and potential.
  313.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  314.      * are related to both the normalized coefficients
  315.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  316.      *  and the J<sub>n</sub> one as follows:</p>
  317.      *
  318.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  319.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  320.      *
  321.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  322.      *
  323.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  324.      *
  325.      * @param initialOrbit initial orbit
  326.      * @param attitudeProv attitude provider
  327.      * @param mass spacecraft mass
  328.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  329.      * @param mu central attraction coefficient (m³/s²)
  330.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  331.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  332.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  333.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  334.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  335.      * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
  336.      *                                 double, double, double, double, double, PropagationType)
  337.      */
  338.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  339.                                      final AttitudeProvider attitudeProv,
  340.                                      final double mass,
  341.                                      final double referenceRadius,
  342.                                      final double mu,
  343.                                      final double c20,
  344.                                      final double c30,
  345.                                      final double c40,
  346.                                      final double c50,
  347.                                      final double c60) {
  348.         this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, c60,
  349.              PropagationType.OSCULATING);
  350.     }


  353.     /** Build a propagator from orbit and potential provider.
  354.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  355.      *
  356.      * <p>Using this constructor, it is possible to define the initial orbit as
  357.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  358.      *
  359.      * @param initialOrbit initial orbit
  360.      * @param provider for un-normalized zonal coefficients
  361.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  362.      * @since 10.2
  363.      */
  364.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  365.                                      final UnnormalizedSphericalHarmonicsProvider provider,
  366.                                      final PropagationType initialType) {
  367.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  368.              DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), initialType);
  369.     }

  371.     /** Build a propagator from orbit, attitude provider, mass and potential provider.
  372.      * <p>Using this constructor, it is possible to define the initial orbit as
  373.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  374.      * @param initialOrbit initial orbit
  375.      * @param attitudeProv attitude provider
  376.      * @param mass spacecraft mass
  377.      * @param provider for un-normalized zonal coefficients
  378.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  379.      * @since 10.2
  380.      */
  381.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  382.                                      final AttitudeProvider attitudeProv,
  383.                                      final double mass,
  384.                                      final UnnormalizedSphericalHarmonicsProvider provider,
  385.                                      final PropagationType initialType) {
  386.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), initialType);
  387.     }

  389.     /**
  390.      * Private helper constructor.
  391.      * <p>Using this constructor, it is possible to define the initial orbit as
  392.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  393.      * @param initialOrbit initial orbit
  394.      * @param attitude attitude provider
  395.      * @param mass spacecraft mass
  396.      * @param provider for un-normalized zonal coefficients
  397.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  398.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  399.      * @since 10.2
  400.      */
  401.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  402.                                      final AttitudeProvider attitude,
  403.                                      final double mass,
  404.                                      final UnnormalizedSphericalHarmonicsProvider provider,
  405.                                      final UnnormalizedSphericalHarmonics harmonics,
  406.                                      final PropagationType initialType) {
  407.         this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
  408.              harmonics.getUnnormalizedCnm(2, 0),
  409.              harmonics.getUnnormalizedCnm(3, 0),
  410.              harmonics.getUnnormalizedCnm(4, 0),
  411.              harmonics.getUnnormalizedCnm(5, 0),
  412.              harmonics.getUnnormalizedCnm(6, 0),
  413.              initialType);
  414.     }

  416.     /** Build a propagator from orbit, attitude provider, mass and potential.
  417.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  418.      * are related to both the normalized coefficients
  419.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  420.      *  and the J<sub>n</sub> one as follows:</p>
  421.      *
  422.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  423.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  424.      *
  425.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  426.      *
  427.      * <p>Using this constructor, it is possible to define the initial orbit as
  428.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  429.      *
  430.      * @param initialOrbit initial orbit
  431.      * @param attitudeProv attitude provider
  432.      * @param mass spacecraft mass
  433.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  434.      * @param mu central attraction coefficient (m³/s²)
  435.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  436.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  437.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  438.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  439.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  440.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  441.      * @since 10.2
  442.      */
  443.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  444.                                      final AttitudeProvider attitudeProv,
  445.                                      final double mass,
  446.                                      final double referenceRadius,
  447.                                      final double mu,
  448.                                      final double c20,
  449.                                      final double c30,
  450.                                      final double c40,
  451.                                      final double c50,
  452.                                      final double c60,
  453.                                      final PropagationType initialType) {
  454.         this(initialOrbit, attitudeProv, mass, referenceRadius, mu,
  455.              c20, c30, c40, c50, c60, initialType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  456.     }

  458.     /** Build a propagator from orbit, attitude provider, mass and potential.
  459.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  460.      * are related to both the normalized coefficients
  461.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  462.      *  and the J<sub>n</sub> one as follows:</p>
  463.      *
  464.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  465.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  466.      *
  467.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  468.      *
  469.      * <p>Using this constructor, it is possible to define the initial orbit as
  470.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  471.      *
  472.      * @param initialOrbit initial orbit
  473.      * @param attitudeProv attitude provider
  474.      * @param mass spacecraft mass
  475.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  476.      * @param mu central attraction coefficient (m³/s²)
  477.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  478.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  479.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  480.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  481.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  482.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  483.      * @param epsilon convergence threshold for mean parameters conversion
  484.      * @param maxIterations maximum iterations for mean parameters conversion
  485.      * @since 11.2
  486.      */
  487.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  488.                                      final AttitudeProvider attitudeProv,
  489.                                      final double mass,
  490.                                      final double referenceRadius,
  491.                                      final double mu,
  492.                                      final double c20,
  493.                                      final double c30,
  494.                                      final double c40,
  495.                                      final double c50,
  496.                                      final double c60,
  497.                                      final PropagationType initialType,
  498.                                      final double epsilon,
  499.                                      final int maxIterations) {
  500.         this(initialOrbit, attitudeProv, mass, referenceRadius, mu,
  501.              c20, c30, c40, c50, c60, initialType,
  502.              new FixedPointConverter(epsilon, maxIterations,
  503.                                      FixedPointConverter.DEFAULT_DAMPING));
  504.     }

  506.     /** Build a propagator from orbit, attitude provider, mass and potential.
  507.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  508.      * are related to both the normalized coefficients
  509.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  510.      *  and the J<sub>n</sub> one as follows:</p>
  511.      *
  512.      * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  513.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  514.      *
  515.      * <p> C<sub>n,0</sub> = -J<sub>n</sub>
  516.      *
  517.      * <p>Using this constructor, it is possible to define the initial orbit as
  518.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  519.      *
  520.      * @param initialOrbit initial orbit
  521.      * @param attitudeProv attitude provider
  522.      * @param mass spacecraft mass
  523.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  524.      * @param mu central attraction coefficient (m³/s²)
  525.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  526.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  527.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  528.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  529.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  530.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  531.      * @param converter osculating to mean orbit converter
  532.      * @since 13.0
  533.      */
  534.     public EcksteinHechlerPropagator(final Orbit initialOrbit,
  535.                                      final AttitudeProvider attitudeProv,
  536.                                      final double mass,
  537.                                      final double referenceRadius,
  538.                                      final double mu,
  539.                                      final double c20,
  540.                                      final double c30,
  541.                                      final double c40,
  542.                                      final double c50,
  543.                                      final double c60,
  544.                                      final PropagationType initialType,
  545.                                      final OsculatingToMeanConverter converter) {
  546.         super(attitudeProv);

  548.         // store model coefficients
  549.         this.referenceRadius = referenceRadius;
  550.         this.mu  = mu;
  551.         this.ck0 = new double[] {
  552.             0.0, 0.0, c20, c30, c40, c50, c60
  553.         };

  555.         // compute mean parameters if needed
  556.         // transform into circular adapted parameters used by the Eckstein-Hechler model
  557.         resetInitialState(new SpacecraftState(initialOrbit,
  558.                                               attitudeProv.getAttitude(initialOrbit,
  559.                                                                        initialOrbit.getDate(),
  560.                                                                        initialOrbit.getFrame())).withMass(mass),
  561.                           initialType, converter);

  563.     }

  565.     /** Conversion from osculating to mean orbit.
  566.      * <p>
  567.      * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
  568.      * osculating SpacecraftState in input.
  569.      * </p>
  570.      * <p>
  571.      * Since the osculating orbit is obtained with the computation of
  572.      * short-periodic variation, the resulting output will depend on
  573.      * the gravity field parameterized in input.
  574.      * </p>
  575.      * <p>
  576.      * The computation is done through a fixed-point iteration process.
  577.      * </p>
  578.      * @param osculating osculating orbit to convert
  579.      * @param provider for un-normalized zonal coefficients
  580.      * @param harmonics {@code provider.onDate(osculating.getDate())}
  581.      * @return mean orbit in a Eckstein-Hechler sense
  582.      * @since 11.2
  583.      */
  584.     public static CircularOrbit computeMeanOrbit(final Orbit osculating,
  585.                                                  final UnnormalizedSphericalHarmonicsProvider provider,
  586.                                                  final UnnormalizedSphericalHarmonics harmonics) {
  587.         return computeMeanOrbit(osculating, provider, harmonics, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  588.     }

  590.     /** Conversion from osculating to mean orbit.
  591.      * <p>
  592.      * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
  593.      * osculating SpacecraftState in input.
  594.      * </p>
  595.      * <p>
  596.      * Since the osculating orbit is obtained with the computation of
  597.      * short-periodic variation, the resulting output will depend on
  598.      * the gravity field parameterized in input.
  599.      * </p>
  600.      * <p>
  601.      * The computation is done through a fixed-point iteration process.
  602.      * </p>
  603.      * @param osculating osculating orbit to convert
  604.      * @param provider for un-normalized zonal coefficients
  605.      * @param harmonics {@code provider.onDate(osculating.getDate())}
  606.      * @param epsilon convergence threshold for mean parameters conversion
  607.      * @param maxIterations maximum iterations for mean parameters conversion
  608.      * @return mean orbit in a Eckstein-Hechler sense
  609.      * @since 11.2
  610.      */
  611.     public static CircularOrbit computeMeanOrbit(final Orbit osculating,
  612.                                                  final UnnormalizedSphericalHarmonicsProvider provider,
  613.                                                  final UnnormalizedSphericalHarmonics harmonics,
  614.                                                  final double epsilon, final int maxIterations) {
  615.         return computeMeanOrbit(osculating,
  616.                                 provider.getAe(), provider.getMu(),
  617.                                 harmonics.getUnnormalizedCnm(2, 0),
  618.                                 harmonics.getUnnormalizedCnm(3, 0),
  619.                                 harmonics.getUnnormalizedCnm(4, 0),
  620.                                 harmonics.getUnnormalizedCnm(5, 0),
  621.                                 harmonics.getUnnormalizedCnm(6, 0),
  622.                                 epsilon, maxIterations);
  623.     }

  625.     /** Conversion from osculating to mean orbit.
  626.      * <p>
  627.      * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
  628.      * osculating SpacecraftState in input.
  629.      * </p>
  630.      * <p>
  631.      * Since the osculating orbit is obtained with the computation of
  632.      * short-periodic variation, the resulting output will depend on
  633.      * the gravity field parameterized in input.
  634.      * </p>
  635.      * <p>
  636.      * The computation is done through a fixed-point iteration process.
  637.      * </p>
  638.      * @param osculating osculating orbit to convert
  639.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  640.      * @param mu central attraction coefficient (m³/s²)
  641.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  642.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  643.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  644.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  645.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  646.      * @param epsilon convergence threshold for mean parameters conversion
  647.      * @param maxIterations maximum iterations for mean parameters conversion
  648.      * @return mean orbit in a Eckstein-Hechler sense
  649.      * @since 11.2
  650.      */
  651.     public static CircularOrbit computeMeanOrbit(final Orbit osculating,
  652.                                                  final double referenceRadius,
  653.                                                  final double mu,
  654.                                                  final double c20,
  655.                                                  final double c30,
  656.                                                  final double c40,
  657.                                                  final double c50,
  658.                                                  final double c60,
  659.                                                  final double epsilon,
  660.                                                  final int maxIterations) {
  661.         // Build a fixed-point converter
  662.         final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
  663.                                                                             FixedPointConverter.DEFAULT_DAMPING);
  664.         return computeMeanOrbit(osculating, referenceRadius, mu, c20, c30, c40, c50, c60, converter);
  665.     }

  667.     /** Conversion from osculating to mean orbit.
  668.      * <p>
  669.      * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
  670.      * osculating SpacecraftState in input.
  671.      * </p>
  672.      * <p>
  673.      * Since the osculating orbit is obtained with the computation of
  674.      * short-periodic variation, the resulting output will depend on
  675.      * the gravity field parameterized in input.
  676.      * </p>
  677.      * <p>
  678.      * The computation is done through a fixed-point iteration process.
  679.      * </p>
  680.      * @param osculating osculating orbit to convert
  681.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  682.      * @param mu central attraction coefficient (m³/s²)
  683.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  684.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  685.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  686.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  687.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  688.      * @param converter osculating to mean orbit converter
  689.      * @return mean orbit in a Eckstein-Hechler sense
  690.      * @since 13.0
  691.      */
  692.     public static CircularOrbit computeMeanOrbit(final Orbit osculating,
  693.                                                  final double referenceRadius,
  694.                                                  final double mu,
  695.                                                  final double c20,
  696.                                                  final double c30,
  697.                                                  final double c40,
  698.                                                  final double c50,
  699.                                                  final double c60,
  700.                                                  final OsculatingToMeanConverter converter) {
  701.         // Set EH as the mean theory for converting
  702.         final MeanTheory theory = new EcksteinHechlerTheory(referenceRadius, mu,
  703.                                                             c20, c30, c40, c50, c60);
  704.         converter.setMeanTheory(theory);
  705.         return (CircularOrbit) OrbitType.CIRCULAR.convertType(converter.convertToMean(osculating));
  706.     }

  708.     /** {@inheritDoc}
  709.      * <p>The new initial state to consider
  710.      * must be defined with an osculating orbit.</p>
  711.      * @see #resetInitialState(SpacecraftState, PropagationType)
  712.      */
  713.     @Override
  714.     public void resetInitialState(final SpacecraftState state) {
  715.         resetInitialState(state, PropagationType.OSCULATING);
  716.     }

  718.     /** Reset the propagator initial state.
  719.      * @param state new initial state to consider
  720.      * @param stateType mean Eckstein-Hechler orbit or osculating orbit
  721.      * @since 10.2
  722.      */
  723.     public void resetInitialState(final SpacecraftState state,
  724.                                   final PropagationType stateType) {
  725.         resetInitialState(state, stateType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  726.     }

  728.     /** Reset the propagator initial state.
  729.      * @param state new initial state to consider
  730.      * @param stateType mean Eckstein-Hechler orbit or osculating orbit
  731.      * @param epsilon convergence threshold for mean parameters conversion
  732.      * @param maxIterations maximum iterations for mean parameters conversion
  733.      * @since 11.2
  734.      */
  735.     public void resetInitialState(final SpacecraftState state,
  736.                                   final PropagationType stateType,
  737.                                   final double epsilon,
  738.                                   final int maxIterations) {
  739.         final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
  740.                                                                             FixedPointConverter.DEFAULT_DAMPING);
  741.         resetInitialState(state, stateType, converter);
  742.     }

  744.     /** Reset the propagator initial state.
  745.      * @param state new initial state to consider
  746.      * @param stateType mean Eckstein-Hechler orbit or osculating orbit
  747.      * @param converter osculating to mean orbit converter
  748.      * @since 13.0
  749.      */
  750.     public void resetInitialState(final SpacecraftState state,
  751.                                   final PropagationType stateType,
  752.                                   final OsculatingToMeanConverter converter) {
  753.         super.resetInitialState(state);
  754.         CircularOrbit circular = (CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit());
  755.         if (stateType == PropagationType.OSCULATING) {
  756.             final MeanTheory theory = new EcksteinHechlerTheory(referenceRadius, mu, ck0[2],
  757.                                                                 ck0[3], ck0[4], ck0[5], ck0[6]);
  758.             converter.setMeanTheory(theory);
  759.             circular = (CircularOrbit) OrbitType.CIRCULAR.convertType(converter.convertToMean(circular));
  760.         }
  761.         this.initialModel = new EHModel(circular, state.getMass(), referenceRadius, mu, ck0);
  762.         this.models = new TimeSpanMap<>(initialModel);
  763.     }

  765.     /** {@inheritDoc} */
  766.     protected void resetIntermediateState(final SpacecraftState state,
  767.                                           final boolean forward) {
  768.         resetIntermediateState(state, forward, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  769.     }

  771.     /** Reset an intermediate state.
  772.      * @param state new intermediate state to consider
  773.      * @param forward if true, the intermediate state is valid for
  774.      * propagations after itself
  775.      * @param epsilon convergence threshold for mean parameters conversion
  776.      * @param maxIterations maximum iterations for mean parameters conversion
  777.      * @since 11.2
  778.      */
  779.     protected void resetIntermediateState(final SpacecraftState state,
  780.                                           final boolean forward,
  781.                                           final double epsilon,
  782.                                           final int maxIterations) {
  783.         final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
  784.                                                                             FixedPointConverter.DEFAULT_DAMPING);
  785.         resetIntermediateState(state, forward, converter);
  786.     }

  788.     /** Reset an intermediate state.
  789.      * @param state     new intermediate state to consider
  790.      * @param forward   if true, the intermediate state is valid for
  791.      *                  propagations after itself
  792.      * @param converter osculating to mean orbit converter
  793.      * @since 13.0
  794.      */
  795.     protected void resetIntermediateState(final SpacecraftState state,
  796.                                           final boolean forward,
  797.                                           final OsculatingToMeanConverter converter) {
  798.         final MeanTheory theory = new EcksteinHechlerTheory(referenceRadius, mu, ck0[2],
  799.                                                             ck0[3], ck0[4], ck0[5], ck0[6]);
  800.         converter.setMeanTheory(theory);
  801.         final CircularOrbit mean = (CircularOrbit) OrbitType.CIRCULAR.convertType(converter.convertToMean(state.getOrbit()));
  802.         final EHModel newModel = new EHModel(mean, state.getMass(), referenceRadius, mu, ck0);
  803.         if (forward) {
  804.             models.addValidAfter(newModel, state.getDate(), false);
  805.         } else {
  806.             models.addValidBefore(newModel, state.getDate(), false);
  807.         }
  808.         stateChanged(state);
  809.     }

  811.     /** {@inheritDoc} */
  812.     public CartesianOrbit propagateOrbit(final AbsoluteDate date) {
  813.         // compute Cartesian parameters, taking derivatives into account
  814.         // to make sure velocity and acceleration are consistent
  815.         final EHModel current = models.get(date);
  816.         return new CartesianOrbit(toCartesian(date, current.propagateParameters(date)),
  817.                                   current.mean.getFrame(), mu);
  818.     }

  820.     /**
  821.      * Get the osculating circular orbit from the EH model.
  822.      * <p>This method is only relevant for the conversion from osculating to mean orbit.</p>
  823.      * @param date target date for the orbit
  824.      * @return the osculating circular orbite
  825.      */
  826.     public CircularOrbit getOsculatingCircularOrbit(final AbsoluteDate date) {
  827.         // compute circular parameters from the model
  828.         final EHModel current = models.get(date);
  829.         final UnivariateDerivative2[] parameter = current.propagateParameters(date);
  830.         return new CircularOrbit(parameter[0].getValue(),
  831.                                  parameter[1].getValue(),
  832.                                  parameter[2].getValue(),
  833.                                  parameter[3].getValue(),
  834.                                  parameter[4].getValue(),
  835.                                  parameter[5].getValue(),
  836.                                  PositionAngleType.MEAN,
  837.                                  current.mean.getFrame(),
  838.                                  date, mu);
  839.     }

  841.     /**
  842.      * Get the central attraction coefficient μ.
  843.      * @return mu central attraction coefficient (m³/s²)
  844.      * @since 11.1
  845.      */
  846.     public double getMu() {
  847.         return mu;
  848.     }

  850.     /**
  851.      * Get the un-normalized zonal coefficients.
  852.      * @return the un-normalized zonal coefficients
  853.      * @since 11.1
  854.      */
  855.     public double[] getCk0() {
  856.         return ck0.clone();
  857.     }

  859.     /**
  860.      * Get the reference radius of the central body attraction model.
  861.      * @return the reference radius in meters
  862.      * @since 11.1
  863.      */
  864.     public double getReferenceRadius() {
  865.         return referenceRadius;
  866.     }

  868.     /** {@inheritDoc} */
  869.     @Override
  870.     protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
  871.                                                         final DoubleArrayDictionary initialJacobianColumns) {
  872.         // Create the harvester
  873.         final EcksteinHechlerHarvester harvester = new EcksteinHechlerHarvester(this, stmName, initialStm, initialJacobianColumns);
  874.         // Update the list of additional state provider
  875.         addAdditionalDataProvider(harvester);
  876.         // Return the configured harvester
  877.         return harvester;
  878.     }

  880.     /** Local class for Eckstein-Hechler model, with fixed mean parameters. */
  881.     private static class EHModel {

  883.         /** Mean orbit. */
  884.         private final CircularOrbit mean;

  886.         /** Constant mass. */
  887.         private final double mass;

  889.         // CHECKSTYLE: stop JavadocVariable check

  891.         // preprocessed values
  892.         private final double xnotDot;
  893.         private final double rdpom;
  894.         private final double rdpomp;
  895.         private final double eps1;
  896.         private final double eps2;
  897.         private final double xim;
  898.         private final double ommD;
  899.         private final double rdl;
  900.         private final double aMD;

  902.         private final double kh;
  903.         private final double kl;

  905.         private final double ax1;
  906.         private final double ay1;
  907.         private final double as1;
  908.         private final double ac2;
  909.         private final double axy3;
  910.         private final double as3;
  911.         private final double ac4;
  912.         private final double as5;
  913.         private final double ac6;

  915.         private final double ex1;
  916.         private final double exx2;
  917.         private final double exy2;
  918.         private final double ex3;
  919.         private final double ex4;

  921.         private final double ey1;
  922.         private final double eyx2;
  923.         private final double eyy2;
  924.         private final double ey3;
  925.         private final double ey4;

  927.         private final double rx1;
  928.         private final double ry1;
  929.         private final double r2;
  930.         private final double r3;
  931.         private final double rl;

  933.         private final double iy1;
  934.         private final double ix1;
  935.         private final double i2;
  936.         private final double i3;
  937.         private final double ih;

  939.         private final double lx1;
  940.         private final double ly1;
  941.         private final double l2;
  942.         private final double l3;
  943.         private final double ll;

  945.         // CHECKSTYLE: resume JavadocVariable check

  947.         /** Create a model for specified mean orbit.
  948.          * @param mean mean orbit
  949.          * @param mass constant mass
  950.          * @param referenceRadius reference radius of the central body attraction model (m)
  951.          * @param mu central attraction coefficient (m³/s²)
  952.          * @param ck0 un-normalized zonal coefficients
  953.          */
  954.         EHModel(final CircularOrbit mean, final double mass,
  955.                 final double referenceRadius, final double mu, final double[] ck0) {

  957.             this.mean            = mean;
  958.             this.mass            = mass;

  960.             // preliminary processing
  961.             double q = referenceRadius / mean.getA();
  962.             double ql = q * q;
  963.             final double g2 = ck0[2] * ql;
  964.             ql *= q;
  965.             final double g3 = ck0[3] * ql;
  966.             ql *= q;
  967.             final double g4 = ck0[4] * ql;
  968.             ql *= q;
  969.             final double g5 = ck0[5] * ql;
  970.             ql *= q;
  971.             final double g6 = ck0[6] * ql;

  973.             final SinCos sc    = FastMath.sinCos(mean.getI());
  974.             final double cosI1 = sc.cos();
  975.             final double sinI1 = sc.sin();
  976.             final double sinI2 = sinI1 * sinI1;
  977.             final double sinI4 = sinI2 * sinI2;
  978.             final double sinI6 = sinI2 * sinI4;

  980.             if (sinI2 < 1.0e-10) {
  981.                 throw new OrekitException(OrekitMessages.ALMOST_EQUATORIAL_ORBIT,
  982.                                           FastMath.toDegrees(mean.getI()));
  983.             }

  985.             if (FastMath.abs(sinI2 - 4.0 / 5.0) < 1.0e-3) {
  986.                 throw new OrekitException(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT,
  987.                                           FastMath.toDegrees(mean.getI()));
  988.             }

  990.             if (mean.getE() > 0.1) {
  991.                 // if 0.005 < e < 0.1 no error is triggered, but accuracy is poor
  992.                 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
  993.                                           mean.getE());
  994.             }

  996.             xnotDot = FastMath.sqrt(mu / mean.getA()) / mean.getA();

  998.             rdpom = -0.75 * g2 * (4.0 - 5.0 * sinI2);
  999.             rdpomp = 7.5 * g4 * (1.0 - 31.0 / 8.0 * sinI2 + 49.0 / 16.0 * sinI4) -
  1000.                     13.125 * g6 * (1.0 - 8.0 * sinI2 + 129.0 / 8.0 * sinI4 - 297.0 / 32.0 * sinI6);

  1002.             q = 3.0 / (32.0 * rdpom);
  1003.             eps1 = q * g4 * sinI2 * (30.0 - 35.0 * sinI2) -
  1004.                     175.0 * q * g6 * sinI2 * (1.0 - 3.0 * sinI2 + 2.0625 * sinI4);
  1005.             q = 3.0 * sinI1 / (8.0 * rdpom);
  1006.             eps2 = q * g3 * (4.0 - 5.0 * sinI2) - q * g5 * (10.0 - 35.0 * sinI2 + 26.25 * sinI4);

  1008.             xim = mean.getI();
  1009.             ommD = cosI1 * (1.50    * g2 - 2.25 * g2 * g2 * (2.5 - 19.0 / 6.0 * sinI2) +
  1010.                             0.9375  * g4 * (7.0 * sinI2 - 4.0) +
  1011.                             3.28125 * g6 * (2.0 - 9.0 * sinI2 + 8.25 * sinI4));

  1013.             rdl = 1.0 - 1.50 * g2 * (3.0 - 4.0 * sinI2);
  1014.             aMD = rdl +
  1015.                     2.25 * g2 * g2 * (9.0 - 263.0 / 12.0 * sinI2 + 341.0 / 24.0 * sinI4) +
  1016.                     15.0 / 16.0 * g4 * (8.0 - 31.0 * sinI2 + 24.5 * sinI4) +
  1017.                     105.0 / 32.0 * g6 * (-10.0 / 3.0 + 25.0 * sinI2 - 48.75 * sinI4 + 27.5 * sinI6);

  1019.             final double qq = -1.5 * g2 / rdl;
  1020.             final double qA   = 0.75 * g2 * g2 * sinI2;
  1021.             final double qB   = 0.25 * g4 * sinI2;
  1022.             final double qC   = 105.0 / 16.0 * g6 * sinI2;
  1023.             final double qD   = -0.75 * g3 * sinI1;
  1024.             final double qE   = 3.75 * g5 * sinI1;
  1025.             kh = 0.375 / rdpom;
  1026.             kl = kh / sinI1;

  1028.             ax1 = qq * (2.0 - 3.5 * sinI2);
  1029.             ay1 = qq * (2.0 - 2.5 * sinI2);
  1030.             as1 = qD * (4.0 - 5.0 * sinI2) +
  1031.                   qE * (2.625 * sinI4 - 3.5 * sinI2 + 1.0);
  1032.             ac2 = qq * sinI2 +
  1033.                   qA * 7.0 * (2.0 - 3.0 * sinI2) +
  1034.                   qB * (15.0 - 17.5 * sinI2) +
  1035.                   qC * (3.0 * sinI2 - 1.0 - 33.0 / 16.0 * sinI4);
  1036.             axy3 = qq * 3.5 * sinI2;
  1037.             as3 = qD * 5.0 / 3.0 * sinI2 +
  1038.                   qE * 7.0 / 6.0 * sinI2 * (1.0 - 1.125 * sinI2);
  1039.             ac4 = qA * sinI2 +
  1040.                   qB * 4.375 * sinI2 +
  1041.                   qC * 0.75 * (1.1 * sinI4 - sinI2);

  1043.             as5 = qE * 21.0 / 80.0 * sinI4;

  1045.             ac6 = qC * -11.0 / 80.0 * sinI4;

  1047.             ex1 = qq * (1.0 - 1.25 * sinI2);
  1048.             exx2 = qq * 0.5 * (3.0 - 5.0 * sinI2);
  1049.             exy2 = qq * (2.0 - 1.5 * sinI2);
  1050.             ex3 = qq * 7.0 / 12.0 * sinI2;
  1051.             ex4 = qq * 17.0 / 8.0 * sinI2;

  1053.             ey1 = qq * (1.0 - 1.75 * sinI2);
  1054.             eyx2 = qq * (1.0 - 3.0 * sinI2);
  1055.             eyy2 = qq * (2.0 * sinI2 - 1.5);
  1056.             ey3 = qq * 7.0 / 12.0 * sinI2;
  1057.             ey4 = qq * 17.0 / 8.0 * sinI2;

  1059.             q  = -qq * cosI1;
  1060.             rx1 =  3.5 * q;
  1061.             ry1 = -2.5 * q;
  1062.             r2 = -0.5 * q;
  1063.             r3 =  7.0 / 6.0 * q;
  1064.             rl = g3 * cosI1 * (4.0 - 15.0 * sinI2) -
  1065.                  2.5 * g5 * cosI1 * (4.0 - 42.0 * sinI2 + 52.5 * sinI4);

  1067.             q = 0.5 * qq * sinI1 * cosI1;
  1068.             iy1 =  q;
  1069.             ix1 = -q;
  1070.             i2 =  q;
  1071.             i3 =  q * 7.0 / 3.0;
  1072.             ih = -g3 * cosI1 * (4.0 - 5.0 * sinI2) +
  1073.                  2.5 * g5 * cosI1 * (4.0 - 14.0 * sinI2 + 10.5 * sinI4);

  1075.             lx1 = qq * (7.0 - 77.0 / 8.0 * sinI2);
  1076.             ly1 = qq * (55.0 / 8.0 * sinI2 - 7.50);
  1077.             l2 = qq * (1.25 * sinI2 - 0.5);
  1078.             l3 = qq * (77.0 / 24.0 * sinI2 - 7.0 / 6.0);
  1079.             ll = g3 * (53.0 * sinI2 - 4.0 - 57.5 * sinI4) +
  1080.                  2.5 * g5 * (4.0 - 96.0 * sinI2 + 269.5 * sinI4 - 183.75 * sinI6);

  1082.         }

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

  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.             // secular effects

  1096.             // eccentricity
  1097.             final UnivariateDerivative2 x   = xnot.multiply(rdpom + rdpomp);
  1098.             final UnivariateDerivative2 cx  = x.cos();
  1099.             final UnivariateDerivative2 sx  = x.sin();
  1100.             final UnivariateDerivative2 exm = cx.multiply(mean.getCircularEx()).
  1101.                                               add(sx.multiply(eps2 - (1.0 - eps1) * mean.getCircularEy()));
  1102.             final UnivariateDerivative2 eym = sx.multiply((1.0 + eps1) * mean.getCircularEx()).
  1103.                                               add(cx.multiply(mean.getCircularEy() - eps2)).
  1104.                                               add(eps2);

  1106.             // no secular effect on inclination

  1108.             // right ascension of ascending node
  1109.             final UnivariateDerivative2 omm = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode() + ommD * xnot.getValue(),
  1110.                                                                                                FastMath.PI),
  1111.                                                                       ommD * xnotDot,
  1112.                                                                       0.0);

  1114.             // latitude argument
  1115.             final UnivariateDerivative2 xlm = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getAlphaM() + aMD * xnot.getValue(),
  1116.                                                                                                  FastMath.PI),
  1117.                                                                         aMD * xnotDot,
  1118.                                                                         0.0);

  1120.             // periodical terms
  1121.             final UnivariateDerivative2 cl1 = xlm.cos();
  1122.             final UnivariateDerivative2 sl1 = xlm.sin();
  1123.             final UnivariateDerivative2 cl2 = cl1.multiply(cl1).subtract(sl1.multiply(sl1));
  1124.             final UnivariateDerivative2 sl2 = cl1.multiply(sl1).add(sl1.multiply(cl1));
  1125.             final UnivariateDerivative2 cl3 = cl2.multiply(cl1).subtract(sl2.multiply(sl1));
  1126.             final UnivariateDerivative2 sl3 = cl2.multiply(sl1).add(sl2.multiply(cl1));
  1127.             final UnivariateDerivative2 cl4 = cl3.multiply(cl1).subtract(sl3.multiply(sl1));
  1128.             final UnivariateDerivative2 sl4 = cl3.multiply(sl1).add(sl3.multiply(cl1));
  1129.             final UnivariateDerivative2 cl5 = cl4.multiply(cl1).subtract(sl4.multiply(sl1));
  1130.             final UnivariateDerivative2 sl5 = cl4.multiply(sl1).add(sl4.multiply(cl1));
  1131.             final UnivariateDerivative2 cl6 = cl5.multiply(cl1).subtract(sl5.multiply(sl1));

  1133.             final UnivariateDerivative2 qh  = eym.subtract(eps2).multiply(kh);
  1134.             final UnivariateDerivative2 ql  = exm.multiply(kl);

  1136.             final UnivariateDerivative2 exmCl1 = exm.multiply(cl1);
  1137.             final UnivariateDerivative2 exmSl1 = exm.multiply(sl1);
  1138.             final UnivariateDerivative2 eymCl1 = eym.multiply(cl1);
  1139.             final UnivariateDerivative2 eymSl1 = eym.multiply(sl1);
  1140.             final UnivariateDerivative2 exmCl2 = exm.multiply(cl2);
  1141.             final UnivariateDerivative2 exmSl2 = exm.multiply(sl2);
  1142.             final UnivariateDerivative2 eymCl2 = eym.multiply(cl2);
  1143.             final UnivariateDerivative2 eymSl2 = eym.multiply(sl2);
  1144.             final UnivariateDerivative2 exmCl3 = exm.multiply(cl3);
  1145.             final UnivariateDerivative2 exmSl3 = exm.multiply(sl3);
  1146.             final UnivariateDerivative2 eymCl3 = eym.multiply(cl3);
  1147.             final UnivariateDerivative2 eymSl3 = eym.multiply(sl3);
  1148.             final UnivariateDerivative2 exmCl4 = exm.multiply(cl4);
  1149.             final UnivariateDerivative2 exmSl4 = exm.multiply(sl4);
  1150.             final UnivariateDerivative2 eymCl4 = eym.multiply(cl4);
  1151.             final UnivariateDerivative2 eymSl4 = eym.multiply(sl4);

  1153.             // semi major axis
  1154.             final UnivariateDerivative2 rda = exmCl1.multiply(ax1).
  1155.                                               add(eymSl1.multiply(ay1)).
  1156.                                               add(sl1.multiply(as1)).
  1157.                                               add(cl2.multiply(ac2)).
  1158.                                               add(exmCl3.add(eymSl3).multiply(axy3)).
  1159.                                               add(sl3.multiply(as3)).
  1160.                                               add(cl4.multiply(ac4)).
  1161.                                               add(sl5.multiply(as5)).
  1162.                                               add(cl6.multiply(ac6));

  1164.             // eccentricity
  1165.             final UnivariateDerivative2 rdex = cl1.multiply(ex1).
  1166.                                                add(exmCl2.multiply(exx2)).
  1167.                                                add(eymSl2.multiply(exy2)).
  1168.                                                add(cl3.multiply(ex3)).
  1169.                                                add(exmCl4.add(eymSl4).multiply(ex4));
  1170.             final UnivariateDerivative2 rdey = sl1.multiply(ey1).
  1171.                                                add(exmSl2.multiply(eyx2)).
  1172.                                                add(eymCl2.multiply(eyy2)).
  1173.                                                add(sl3.multiply(ey3)).
  1174.                                                add(exmSl4.subtract(eymCl4).multiply(ey4));

  1176.             // ascending node
  1177.             final UnivariateDerivative2 rdom = exmSl1.multiply(rx1).
  1178.                                                add(eymCl1.multiply(ry1)).
  1179.                                                add(sl2.multiply(r2)).
  1180.                                                add(eymCl3.subtract(exmSl3).multiply(r3)).
  1181.                                                add(ql.multiply(rl));

  1183.             // inclination
  1184.             final UnivariateDerivative2 rdxi = eymSl1.multiply(iy1).
  1185.                                                add(exmCl1.multiply(ix1)).
  1186.                                                add(cl2.multiply(i2)).
  1187.                                                add(exmCl3.add(eymSl3).multiply(i3)).
  1188.                                                add(qh.multiply(ih));

  1190.             // latitude argument
  1191.             final UnivariateDerivative2 rdxl = exmSl1.multiply(lx1).
  1192.                                                add(eymCl1.multiply(ly1)).
  1193.                                                add(sl2.multiply(l2)).
  1194.                                                add(exmSl3.subtract(eymCl3).multiply(l3)).
  1195.                                                add(ql.multiply(ll));

  1197.             // osculating parameters
  1198.             return new UnivariateDerivative2[] {
  1199.                 rda.add(1.0).multiply(mean.getA()),
  1200.                 rdex.add(exm),
  1201.                 rdey.add(eym),
  1202.                 rdxi.add(xim),
  1203.                 rdom.add(omm),
  1204.                 rdxl.add(xlm)
  1205.             };

  1207.         }

  1209.     }

  1211.     /** Convert circular parameters <em>with derivatives</em> to Cartesian coordinates.
  1212.      * @param date date of the orbital parameters
  1213.      * @param parameters circular parameters (a, ex, ey, i, raan, alphaM)
  1214.      * @return Cartesian coordinates consistent with values and derivatives
  1215.      */
  1216.     private TimeStampedPVCoordinates toCartesian(final AbsoluteDate date, final UnivariateDerivative2[] parameters) {

  1218.         // evaluate coordinates in the orbit canonical reference frame
  1219.         final UnivariateDerivative2 cosOmega = parameters[4].cos();
  1220.         final UnivariateDerivative2 sinOmega = parameters[4].sin();
  1221.         final UnivariateDerivative2 cosI     = parameters[3].cos();
  1222.         final UnivariateDerivative2 sinI     = parameters[3].sin();
  1223.         final UnivariateDerivative2 alphaE   = meanToEccentric(parameters[5], parameters[1], parameters[2]);
  1224.         final UnivariateDerivative2 cosAE    = alphaE.cos();
  1225.         final UnivariateDerivative2 sinAE    = alphaE.sin();
  1226.         final UnivariateDerivative2 ex2      = parameters[1].square();
  1227.         final UnivariateDerivative2 ey2      = parameters[2].square();
  1228.         final UnivariateDerivative2 exy      = parameters[1].multiply(parameters[2]);
  1229.         final UnivariateDerivative2 q        = ex2.add(ey2).subtract(1).negate().sqrt();
  1230.         final UnivariateDerivative2 beta     = q.add(1).reciprocal();
  1231.         final UnivariateDerivative2 bx2      = beta.multiply(ex2);
  1232.         final UnivariateDerivative2 by2      = beta.multiply(ey2);
  1233.         final UnivariateDerivative2 bxy      = beta.multiply(exy);
  1234.         final UnivariateDerivative2 u        = bxy.multiply(sinAE).subtract(parameters[1].add(by2.subtract(1).multiply(cosAE)));
  1235.         final UnivariateDerivative2 v        = bxy.multiply(cosAE).subtract(parameters[2].add(bx2.subtract(1).multiply(sinAE)));
  1236.         final UnivariateDerivative2 x        = parameters[0].multiply(u);
  1237.         final UnivariateDerivative2 y        = parameters[0].multiply(v);

  1239.         // canonical orbit reference frame
  1240.         final FieldVector3D<UnivariateDerivative2> p =
  1241.                 new FieldVector3D<>(x.multiply(cosOmega).subtract(y.multiply(cosI.multiply(sinOmega))),
  1242.                                     x.multiply(sinOmega).add(y.multiply(cosI.multiply(cosOmega))),
  1243.                                     y.multiply(sinI));

  1245.         // dispatch derivatives
  1246.         final Vector3D p0 = new Vector3D(p.getX().getValue(),
  1247.                                          p.getY().getValue(),
  1248.                                          p.getZ().getValue());
  1249.         final Vector3D p1 = new Vector3D(p.getX().getFirstDerivative(),
  1250.                                          p.getY().getFirstDerivative(),
  1251.                                          p.getZ().getFirstDerivative());
  1252.         final Vector3D p2 = new Vector3D(p.getX().getSecondDerivative(),
  1253.                                          p.getY().getSecondDerivative(),
  1254.                                          p.getZ().getSecondDerivative());
  1255.         return new TimeStampedPVCoordinates(date, p0, p1, p2);

  1257.     }

  1259.     /** Computes the eccentric latitude argument from the mean latitude argument.
  1260.      * @param alphaM = M + Ω mean latitude argument (rad)
  1261.      * @param ex e cos(Ω), first component of circular eccentricity vector
  1262.      * @param ey e sin(Ω), second component of circular eccentricity vector
  1263.      * @return the eccentric latitude argument.
  1264.      */
  1265.     private UnivariateDerivative2 meanToEccentric(final UnivariateDerivative2 alphaM,
  1266.                                                   final UnivariateDerivative2 ex,
  1267.                                                   final UnivariateDerivative2 ey) {
  1268.         // Generalization of Kepler equation to circular parameters
  1269.         // with alphaE = PA + E and
  1270.         //      alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
  1271.         UnivariateDerivative2 alphaE        = alphaM;
  1272.         UnivariateDerivative2 shift         = alphaM.getField().getZero();
  1273.         UnivariateDerivative2 alphaEMalphaM = alphaM.getField().getZero();
  1274.         UnivariateDerivative2 cosAlphaE     = alphaE.cos();
  1275.         UnivariateDerivative2 sinAlphaE     = alphaE.sin();
  1276.         int                 iter          = 0;
  1277.         do {
  1278.             final UnivariateDerivative2 f2 = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
  1279.             final UnivariateDerivative2 f1 = alphaM.getField().getOne().subtract(ex.multiply(cosAlphaE)).subtract(ey.multiply(sinAlphaE));
  1280.             final UnivariateDerivative2 f0 = alphaEMalphaM.subtract(f2);

  1282.             final UnivariateDerivative2 f12 = f1.multiply(2);
  1283.             shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));

  1285.             alphaEMalphaM  = alphaEMalphaM.subtract(shift);
  1286.             alphaE         = alphaM.add(alphaEMalphaM);
  1287.             cosAlphaE      = alphaE.cos();
  1288.             sinAlphaE      = alphaE.sin();

  1290.         } while (++iter < 50 && FastMath.abs(shift.getValue()) > 1.0e-12);

  1292.         return alphaE;

  1294.     }

  1296.     /** {@inheritDoc} */
  1297.     protected double getMass(final AbsoluteDate date) {
  1298.         return models.get(date).mass;
  1299.     }

  1301. }
Created with JaCoCo 0.8.12.202403310830