FieldEcksteinHechlerPropagator.java

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

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

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.Field;
  24. import org.hipparchus.analysis.differentiation.FieldUnivariateDerivative2;
  25. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.FieldSinCos;
  28. import org.hipparchus.util.MathArrays;
  29. import org.hipparchus.util.MathUtils;
  30. import org.orekit.attitudes.AttitudeProvider;
  31. import org.orekit.attitudes.FrameAlignedProvider;
  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.FieldCartesianOrbit;
  37. import org.orekit.orbits.FieldCircularOrbit;
  38. import org.orekit.orbits.FieldOrbit;
  39. import org.orekit.orbits.OrbitType;
  40. import org.orekit.orbits.PositionAngleType;
  41. import org.orekit.propagation.FieldSpacecraftState;
  42. import org.orekit.propagation.PropagationType;
  43. import org.orekit.time.FieldAbsoluteDate;
  44. import org.orekit.utils.FieldTimeSpanMap;
  45. import org.orekit.utils.ParameterDriver;
  46. import org.orekit.utils.TimeStampedFieldPVCoordinates;

  48. /** This class propagates a {@link org.orekit.propagation.FieldSpacecraftState}
  49.  *  using the analytical Eckstein-Hechler model.
  50.  * <p>The Eckstein-Hechler model is suited for near circular orbits
  51.  * (e &lt; 0.1, with poor accuracy between 0.005 and 0.1) and inclination
  52.  * neither equatorial (direct or retrograde) nor critical (direct or
  53.  * retrograde).</p>
  54.  * @see FieldOrbit
  55.  * @author Guylaine Prat
  56.  * @param <T> type of the field elements
  57.  */
  58. public class FieldEcksteinHechlerPropagator<T extends CalculusFieldElement<T>> extends FieldAbstractAnalyticalPropagator<T> {

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

  63.     /** Default value for maxIterations. */
  64.     private static final int MAX_ITERATIONS_DEFAULT = 100;

  66.     /** Initial Eckstein-Hechler model. */
  67.     private FieldEHModel<T> initialModel;

  69.     /** All models. */
  70.     private transient FieldTimeSpanMap<FieldEHModel<T>, T> models;

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

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

  78.     /** Un-normalized zonal coefficients. */
  79.     private double[] ck0;

  81.     /** Build a propagator from FieldOrbit and potential provider.
  82.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  83.      *
  84.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  85.      *
  86.      * @param initialOrbit initial FieldOrbit
  87.      * @param provider for un-normalized zonal coefficients
  88.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, AttitudeProvider,
  89.      * UnnormalizedSphericalHarmonicsProvider)
  90.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, UnnormalizedSphericalHarmonicsProvider,
  91.      * PropagationType)
  92.      */
  93.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  94.                                           final UnnormalizedSphericalHarmonicsProvider provider) {
  95.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  96.              initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
  97.              provider.onDate(initialOrbit.getDate().toAbsoluteDate()));
  98.     }

  100.     /**
  101.      * Private helper constructor.
  102.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  103.      * @param initialOrbit initial FieldOrbit
  104.      * @param attitude attitude provider
  105.      * @param mass spacecraft mass
  106.      * @param provider for un-normalized zonal coefficients
  107.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  108.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement,
  109.      * UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics, PropagationType)
  110.      */
  111.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  112.                                           final AttitudeProvider attitude,
  113.                                           final T mass,
  114.                                           final UnnormalizedSphericalHarmonicsProvider provider,
  115.                                           final UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics harmonics) {
  116.         this(initialOrbit, attitude,  mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
  117.              harmonics.getUnnormalizedCnm(2, 0),
  118.              harmonics.getUnnormalizedCnm(3, 0),
  119.              harmonics.getUnnormalizedCnm(4, 0),
  120.              harmonics.getUnnormalizedCnm(5, 0),
  121.              harmonics.getUnnormalizedCnm(6, 0));
  122.     }

  124.     /** Build a propagator from FieldOrbit and potential.
  125.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  126.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  127.      * are related to both the normalized coefficients
  128.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  129.      *  and the J<sub>n</sub> one as follows:
  130.      * <p>
  131.      *   C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  132.      *                      <span style="text-decoration: overline">C</span><sub>n,0</sub>
  133.      * <p>
  134.      *   C<sub>n,0</sub> = -J<sub>n</sub>
  135.      *
  136.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  137.      *
  138.      * @param initialOrbit initial FieldOrbit
  139.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  140.      * @param mu central attraction coefficient (m³/s²)
  141.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  142.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  143.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  144.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  145.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  146.      * @see org.orekit.utils.Constants
  147.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, AttitudeProvider, double,
  148.      * CalculusFieldElement, double, double, double, double, double)
  149.      */
  150.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  151.                                           final double referenceRadius, final T mu,
  152.                                           final double c20, final double c30, final double c40,
  153.                                           final double c50, final double c60) {
  154.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  155.              initialOrbit.getMu().newInstance(DEFAULT_MASS),
  156.              referenceRadius, mu, c20, c30, c40, c50, c60);
  157.     }

  159.     /** Build a propagator from FieldOrbit, mass and potential provider.
  160.      * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
  161.      *
  162.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  163.      *
  164.      * @param initialOrbit initial FieldOrbit
  165.      * @param mass spacecraft mass
  166.      * @param provider for un-normalized zonal coefficients
  167.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, AttitudeProvider,
  168.      * CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider)
  169.      */
  170.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit, final T mass,
  171.                                           final UnnormalizedSphericalHarmonicsProvider provider) {
  172.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  173.              mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()));
  174.     }

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

  210.     /** Build a propagator from FieldOrbit, attitude provider and potential provider.
  211.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  212.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  213.      * @param initialOrbit initial FieldOrbit
  214.      * @param attitudeProv attitude provider
  215.      * @param provider for un-normalized zonal coefficients
  216.      */
  217.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  218.                                           final AttitudeProvider attitudeProv,
  219.                                           final UnnormalizedSphericalHarmonicsProvider provider) {
  220.         this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS), provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()));
  221.     }

  223.     /** Build a propagator from FieldOrbit, attitude provider and potential.
  224.      * <p>Mass is set to an unspecified non-null arbitrary value.</p>
  225.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  226.      * are related to both the normalized coefficients
  227.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  228.      *  and the J<sub>n</sub> one as follows:</p>
  229.      * <p>
  230.      *   C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  231.      *                     <span style="text-decoration: overline">C</span><sub>n,0</sub>
  232.      * <p>
  233.      *   C<sub>n,0</sub> = -J<sub>n</sub>
  234.      *
  235.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  236.      *
  237.      * @param initialOrbit initial FieldOrbit
  238.      * @param attitudeProv attitude provider
  239.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  240.      * @param mu central attraction coefficient (m³/s²)
  241.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  242.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  243.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  244.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  245.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  246.      */
  247.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  248.                                           final AttitudeProvider attitudeProv,
  249.                                           final double referenceRadius, final T mu,
  250.                                           final double c20, final double c30, final double c40,
  251.                                           final double c50, final double c60) {
  252.         this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS),
  253.              referenceRadius, mu, c20, c30, c40, c50, c60);
  254.     }

  256.     /** Build a propagator from FieldOrbit, attitude provider, mass and potential provider.
  257.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  258.      * @param initialOrbit initial FieldOrbit
  259.      * @param attitudeProv attitude provider
  260.      * @param mass spacecraft mass
  261.      * @param provider for un-normalized zonal coefficients
  262.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement,
  263.      * UnnormalizedSphericalHarmonicsProvider, PropagationType)
  264.      */
  265.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  266.                                           final AttitudeProvider attitudeProv,
  267.                                           final T mass,
  268.                                           final UnnormalizedSphericalHarmonicsProvider provider) {
  269.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()));
  270.     }

  272.     /** Build a propagator from FieldOrbit, attitude provider, mass and potential.
  273.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  274.      * are related to both the normalized coefficients
  275.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  276.      *  and the J<sub>n</sub> one as follows:</p>
  277.      * <p>
  278.      *   C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  279.      *                      <span style="text-decoration: overline">C</span><sub>n,0</sub>
  280.      * <p>
  281.      *   C<sub>n,0</sub> = -J<sub>n</sub>
  282.      *
  283.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  284.      *
  285.      * @param initialOrbit initial FieldOrbit
  286.      * @param attitudeProv attitude provider
  287.      * @param mass spacecraft mass
  288.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  289.      * @param mu central attraction coefficient (m³/s²)
  290.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  291.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  292.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  293.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  294.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  295.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, double,
  296.      *                                      CalculusFieldElement, double, double, double, double, double, PropagationType)
  297.      */
  298.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  299.                                           final AttitudeProvider attitudeProv,
  300.                                           final T mass,
  301.                                           final double referenceRadius, final T mu,
  302.                                           final double c20, final double c30, final double c40,
  303.                                           final double c50, final double c60) {
  304.         this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, c60, PropagationType.OSCULATING);
  305.     }

  307.     /** Build a propagator from orbit and potential provider.
  308.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  309.      *
  310.      * <p>Using this constructor, it is possible to define the initial orbit as
  311.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  312.      *
  313.      * @param initialOrbit initial orbit
  314.      * @param provider for un-normalized zonal coefficients
  315.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  316.      * @since 10.2
  317.      */
  318.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  319.                                           final UnnormalizedSphericalHarmonicsProvider provider,
  320.                                           final PropagationType initialType) {
  321.         this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
  322.              initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
  323.              provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType);
  324.     }

  326.     /** Build a propagator from orbit, attitude provider, mass and potential provider.
  327.      * <p>Using this constructor, it is possible to define the initial orbit as
  328.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  329.      * @param initialOrbit initial orbit
  330.      * @param attitudeProv attitude provider
  331.      * @param mass spacecraft mass
  332.      * @param provider for un-normalized zonal coefficients
  333.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  334.      * @since 10.2
  335.      */
  336.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  337.                                           final AttitudeProvider attitudeProv,
  338.                                           final T mass,
  339.                                           final UnnormalizedSphericalHarmonicsProvider provider,
  340.                                           final PropagationType initialType) {
  341.         this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType);
  342.     }

  344.     /**
  345.      * Private helper constructor.
  346.      * <p>Using this constructor, it is possible to define the initial orbit as
  347.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  348.      * @param initialOrbit initial orbit
  349.      * @param attitude attitude provider
  350.      * @param mass spacecraft mass
  351.      * @param provider for un-normalized zonal coefficients
  352.      * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
  353.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  354.      * @since 10.2
  355.      */
  356.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  357.                                           final AttitudeProvider attitude,
  358.                                           final T mass,
  359.                                           final UnnormalizedSphericalHarmonicsProvider provider,
  360.                                           final UnnormalizedSphericalHarmonics harmonics,
  361.                                           final PropagationType initialType) {
  362.         this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
  363.              harmonics.getUnnormalizedCnm(2, 0),
  364.              harmonics.getUnnormalizedCnm(3, 0),
  365.              harmonics.getUnnormalizedCnm(4, 0),
  366.              harmonics.getUnnormalizedCnm(5, 0),
  367.              harmonics.getUnnormalizedCnm(6, 0),
  368.              initialType);
  369.     }

  371.     /** Build a propagator from FieldOrbit, attitude provider, mass and potential.
  372.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  373.      * are related to both the normalized coefficients
  374.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  375.      *  and the J<sub>n</sub> one as follows:</p>
  376.      * <p>
  377.      *   C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  378.      *                      <span style="text-decoration: overline">C</span><sub>n,0</sub>
  379.      * <p>
  380.      *   C<sub>n,0</sub> = -J<sub>n</sub>
  381.      *
  382.      * <p>Using this constructor, it is possible to define the initial orbit as
  383.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  384.      *
  385.      * @param initialOrbit initial FieldOrbit
  386.      * @param attitudeProv attitude provider
  387.      * @param mass spacecraft mass
  388.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  389.      * @param mu central attraction coefficient (m³/s²)
  390.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  391.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  392.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  393.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  394.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  395.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  396.      * @since 10.2
  397.      */
  398.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  399.                                           final AttitudeProvider attitudeProv,
  400.                                           final T mass,
  401.                                           final double referenceRadius, final T mu,
  402.                                           final double c20, final double c30, final double c40,
  403.                                           final double c50, final double c60,
  404.                                           final PropagationType initialType) {
  405.         this(initialOrbit, attitudeProv, mass, referenceRadius, mu,
  406.              c20, c30, c40, c50, c60, initialType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  407.     }

  409.     /** Build a propagator from FieldOrbit, attitude provider, mass and potential.
  410.      * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
  411.      * are related to both the normalized coefficients
  412.      * <span style="text-decoration: overline">C</span><sub>n,0</sub>
  413.      *  and the J<sub>n</sub> one as follows:</p>
  414.      * <p>
  415.      *   C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
  416.      *                      <span style="text-decoration: overline">C</span><sub>n,0</sub>
  417.      * <p>
  418.      *   C<sub>n,0</sub> = -J<sub>n</sub>
  419.      *
  420.      * <p>Using this constructor, it is possible to define the initial orbit as
  421.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  422.      *
  423.      * @param initialOrbit initial FieldOrbit
  424.      * @param attitudeProv attitude provider
  425.      * @param mass spacecraft mass
  426.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  427.      * @param mu central attraction coefficient (m³/s²)
  428.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  429.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  430.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  431.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  432.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  433.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  434.      * @param epsilon convergence threshold for mean parameters conversion
  435.      * @param maxIterations maximum iterations for mean parameters conversion
  436.      * @since 11.2
  437.      */
  438.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  439.                                           final AttitudeProvider attitudeProv,
  440.                                           final T mass,
  441.                                           final double referenceRadius, final T mu,
  442.                                           final double c20, final double c30, final double c40,
  443.                                           final double c50, final double c60,
  444.                                           final PropagationType initialType,
  445.                                           final double epsilon, final int maxIterations) {
  446.         super(mass.getField(), attitudeProv);
  447.         try {

  449.             // store model coefficients
  450.             this.referenceRadius = referenceRadius;
  451.             this.mu  = mu;
  452.             this.ck0 = new double[] {
  453.                 0.0, 0.0, c20, c30, c40, c50, c60
  454.             };

  456.             // compute mean parameters if needed
  457.             // transform into circular adapted parameters used by the Eckstein-Hechler model
  458.             resetInitialState(new FieldSpacecraftState<>(initialOrbit,
  459.                                                          attitudeProv.getAttitude(initialOrbit,
  460.                                                                                   initialOrbit.getDate(),
  461.                                                                                   initialOrbit.getFrame()),
  462.                                                          mass),
  463.                               initialType, epsilon, maxIterations);

  465.         } catch (OrekitException oe) {
  466.             throw new OrekitException(oe);
  467.         }
  468.     }

  470.     /** Conversion from osculating to mean orbit.
  471.      * <p>
  472.      * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
  473.      * osculating SpacecraftState in input.
  474.      * </p>
  475.      * <p>
  476.      * Since the osculating orbit is obtained with the computation of
  477.      * short-periodic variation, the resulting output will depend on
  478.      * the gravity field parameterized in input.
  479.      * </p>
  480.      * <p>
  481.      * The computation is done through a fixed-point iteration process.
  482.      * </p>
  483.      * @param <T> type of the filed elements
  484.      * @param osculating osculating orbit to convert
  485.      * @param provider for un-normalized zonal coefficients
  486.      * @param harmonics {@code provider.onDate(osculating.getDate())}
  487.      * @return mean orbit in a Eckstein-Hechler sense
  488.      * @since 11.2
  489.      */
  490.     public static <T extends CalculusFieldElement<T>> FieldCircularOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
  491.                                                                                              final UnnormalizedSphericalHarmonicsProvider provider,
  492.                                                                                              final UnnormalizedSphericalHarmonics harmonics) {
  493.         return computeMeanOrbit(osculating, provider, harmonics, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  494.     }

  496.     /** Conversion from osculating to mean orbit.
  497.      * <p>
  498.      * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
  499.      * osculating SpacecraftState in input.
  500.      * </p>
  501.      * <p>
  502.      * Since the osculating orbit is obtained with the computation of
  503.      * short-periodic variation, the resulting output will depend on
  504.      * the gravity field parameterized in input.
  505.      * </p>
  506.      * <p>
  507.      * The computation is done through a fixed-point iteration process.
  508.      * </p>
  509.      * @param <T> type of the filed elements
  510.      * @param osculating osculating orbit to convert
  511.      * @param provider for un-normalized zonal coefficients
  512.      * @param harmonics {@code provider.onDate(osculating.getDate())}
  513.      * @param epsilon convergence threshold for mean parameters conversion
  514.      * @param maxIterations maximum iterations for mean parameters conversion
  515.      * @return mean orbit in a Eckstein-Hechler sense
  516.      * @since 11.2
  517.      */
  518.     public static <T extends CalculusFieldElement<T>> FieldCircularOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
  519.                                                                                              final UnnormalizedSphericalHarmonicsProvider provider,
  520.                                                                                              final UnnormalizedSphericalHarmonics harmonics,
  521.                                                                                              final double epsilon, final int maxIterations) {
  522.         return computeMeanOrbit(osculating,
  523.                                 provider.getAe(), provider.getMu(),
  524.                                 harmonics.getUnnormalizedCnm(2, 0),
  525.                                 harmonics.getUnnormalizedCnm(3, 0),
  526.                                 harmonics.getUnnormalizedCnm(4, 0),
  527.                                 harmonics.getUnnormalizedCnm(5, 0),
  528.                                 harmonics.getUnnormalizedCnm(6, 0),
  529.                                 epsilon, maxIterations);
  530.     }

  532.     /** Conversion from osculating to mean orbit.
  533.      * <p>
  534.      * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
  535.      * osculating SpacecraftState in input.
  536.      * </p>
  537.      * <p>
  538.      * Since the osculating orbit is obtained with the computation of
  539.      * short-periodic variation, the resulting output will depend on
  540.      * the gravity field parameterized in input.
  541.      * </p>
  542.      * <p>
  543.      * The computation is done through a fixed-point iteration process.
  544.      * </p>
  545.      * @param <T> type of the filed elements
  546.      * @param osculating osculating orbit to convert
  547.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  548.      * @param mu central attraction coefficient (m³/s²)
  549.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  550.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  551.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  552.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  553.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  554.      * @param epsilon convergence threshold for mean parameters conversion
  555.      * @param maxIterations maximum iterations for mean parameters conversion
  556.      * @return mean orbit in a Eckstein-Hechler sense
  557.      * @since 11.2
  558.      */
  559.     public static <T extends CalculusFieldElement<T>> FieldCircularOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
  560.                                                                                              final double referenceRadius, final double mu,
  561.                                                                                              final double c20, final double c30, final double c40,
  562.                                                                                              final double c50, final double c60,
  563.                                                                                              final double epsilon, final int maxIterations) {
  564.         final FieldEcksteinHechlerPropagator<T> propagator =
  565.                         new FieldEcksteinHechlerPropagator<>(osculating,
  566.                                                              FrameAlignedProvider.of(osculating.getFrame()),
  567.                                                              osculating.getMu().newInstance(DEFAULT_MASS),
  568.                                                              referenceRadius, osculating.getMu().newInstance(mu),
  569.                                                              c20, c30, c40, c50, c60,
  570.                                                              PropagationType.OSCULATING, epsilon, maxIterations);
  571.         return propagator.initialModel.mean;
  572.     }

  574.     /** {@inheritDoc}
  575.      * <p>The new initial state to consider
  576.      * must be defined with an osculating orbit.</p>
  577.      * @see #resetInitialState(FieldSpacecraftState, PropagationType)
  578.      */
  579.     @Override
  580.     public void resetInitialState(final FieldSpacecraftState<T> state) {
  581.         resetInitialState(state, PropagationType.OSCULATING);
  582.     }

  584.     /** Reset the propagator initial state.
  585.      * @param state new initial state to consider
  586.      * @param stateType mean Eckstein-Hechler orbit or osculating orbit
  587.      * @since 10.2
  588.      */
  589.     public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType) {
  590.         resetInitialState(state, stateType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  591.     }

  593.     /** Reset the propagator initial state.
  594.      * @param state new initial state to consider
  595.      * @param stateType mean Eckstein-Hechler orbit or osculating orbit
  596.      * @param epsilon convergence threshold for mean parameters conversion
  597.      * @param maxIterations maximum iterations for mean parameters conversion
  598.      * @since 11.2
  599.      */
  600.     public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType,
  601.                                   final double epsilon, final int maxIterations) {
  602.         super.resetInitialState(state);
  603.         final FieldCircularOrbit<T> circular = (FieldCircularOrbit<T>) OrbitType.CIRCULAR.convertType(state.getOrbit());
  604.         this.initialModel = (stateType == PropagationType.MEAN) ?
  605.                              new FieldEHModel<>(circular, state.getMass(), referenceRadius, mu, ck0) :
  606.                              computeMeanParameters(circular, state.getMass(), epsilon, maxIterations);
  607.         this.models = new FieldTimeSpanMap<FieldEHModel<T>, T>(initialModel, state.getA().getField());
  608.     }

  610.     /** {@inheritDoc} */
  611.     @Override
  612.     protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward) {
  613.         resetIntermediateState(state, forward, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  614.     }

  616.     /** Reset an intermediate state.
  617.      * @param state new intermediate state to consider
  618.      * @param forward if true, the intermediate state is valid for
  619.      * propagations after itself
  620.      * @param epsilon convergence threshold for mean parameters conversion
  621.      * @param maxIterations maximum iterations for mean parameters conversion
  622.      * @since 11.2
  623.      */
  624.     protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward,
  625.                                           final double epsilon, final int maxIterations) {
  626.         final FieldEHModel<T> newModel = computeMeanParameters((FieldCircularOrbit<T>) OrbitType.CIRCULAR.convertType(state.getOrbit()),
  627.                                                                state.getMass(), epsilon, maxIterations);
  628.         if (forward) {
  629.             models.addValidAfter(newModel, state.getDate());
  630.         } else {
  631.             models.addValidBefore(newModel, state.getDate());
  632.         }
  633.         stateChanged(state);
  634.     }

  636.     /** Compute mean parameters according to the Eckstein-Hechler analytical model.
  637.      * @param osculating osculating FieldOrbit
  638.      * @param mass constant mass
  639.      * @param epsilon convergence threshold for mean parameters conversion
  640.      * @param maxIterations maximum iterations for mean parameters conversion
  641.      * @return Eckstein-Hechler mean model
  642.      */
  643.     private FieldEHModel<T> computeMeanParameters(final FieldCircularOrbit<T> osculating, final T mass,
  644.                                                   final double epsilon, final int maxIterations) {

  646.         // sanity check
  647.         if (osculating.getA().getReal() < referenceRadius) {
  648.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
  649.                                            osculating.getA());
  650.         }
  651.         final Field<T> field = mass.getField();
  652.         final T one = field.getOne();
  653.         final T zero = field.getZero();
  654.         // rough initialization of the mean parameters
  655.         FieldEHModel<T> current = new FieldEHModel<>(osculating, mass, referenceRadius, mu, ck0);
  656.         // threshold for each parameter
  657.         final T thresholdA      = current.mean.getA().abs().add(1.0).multiply(epsilon);
  658.         final T thresholdE      = current.mean.getE().add(1.0).multiply(epsilon);
  659.         final T thresholdAngles = one.getPi().multiply(2).multiply(epsilon);


  662.         int i = 0;
  663.         while (i++ < maxIterations) {

  665.             // recompute the osculating parameters from the current mean parameters
  666.             final FieldUnivariateDerivative2<T>[] parameters = current.propagateParameters(current.mean.getDate());
  667.             // adapted parameters residuals
  668.             final T deltaA      = osculating.getA()         .subtract(parameters[0].getValue());
  669.             final T deltaEx     = osculating.getCircularEx().subtract(parameters[1].getValue());
  670.             final T deltaEy     = osculating.getCircularEy().subtract(parameters[2].getValue());
  671.             final T deltaI      = osculating.getI()         .subtract(parameters[3].getValue());
  672.             final T deltaRAAN   = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode().subtract(
  673.                                                                 parameters[4].getValue()),
  674.                                                                 zero);
  675.             final T deltaAlphaM = MathUtils.normalizeAngle(osculating.getAlphaM().subtract(parameters[5].getValue()), zero);
  676.             // update mean parameters
  677.             current = new FieldEHModel<>(new FieldCircularOrbit<>(current.mean.getA().add(deltaA),
  678.                                                                   current.mean.getCircularEx().add( deltaEx),
  679.                                                                   current.mean.getCircularEy().add( deltaEy),
  680.                                                                   current.mean.getI()         .add( deltaI ),
  681.                                                                   current.mean.getRightAscensionOfAscendingNode().add(deltaRAAN),
  682.                                                                   current.mean.getAlphaM().add(deltaAlphaM),
  683.                                                                   PositionAngleType.MEAN,
  684.                                                                   current.mean.getFrame(),
  685.                                                                   current.mean.getDate(), mu),
  686.                                   mass, referenceRadius, mu, ck0);
  687.             // check convergence
  688.             if (FastMath.abs(deltaA.getReal())      < thresholdA.getReal() &&
  689.                 FastMath.abs(deltaEx.getReal())     < thresholdE.getReal() &&
  690.                 FastMath.abs(deltaEy.getReal())     < thresholdE.getReal() &&
  691.                 FastMath.abs(deltaI.getReal())      < thresholdAngles.getReal() &&
  692.                 FastMath.abs(deltaRAAN.getReal())   < thresholdAngles.getReal() &&
  693.                 FastMath.abs(deltaAlphaM.getReal()) < thresholdAngles.getReal()) {
  694.                 return current;
  695.             }

  697.         }

  699.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_ECKSTEIN_HECHLER_MEAN_PARAMETERS, i);

  701.     }

  703.     /** {@inheritDoc} */
  704.     @Override
  705.     public FieldCartesianOrbit<T> propagateOrbit(final FieldAbsoluteDate<T> date, final T[] parameters) {
  706.         // compute Cartesian parameters, taking derivatives into account
  707.         // to make sure velocity and acceleration are consistent
  708.         final FieldEHModel<T> current = models.get(date);
  709.         return new FieldCartesianOrbit<>(toCartesian(date, current.propagateParameters(date)),
  710.                                          current.mean.getFrame(), mu);
  711.     }

  713.     /** Local class for Eckstein-Hechler model, with fixed mean parameters. */
  714.     private static class FieldEHModel<T extends CalculusFieldElement<T>> {

  716.         /** Mean FieldOrbit. */
  717.         private final FieldCircularOrbit<T> mean;

  719.         /** Constant mass. */
  720.         private final T mass;
  721.         // CHECKSTYLE: stop JavadocVariable check

  723.         // preprocessed values
  724.         private final T xnotDot;
  725.         private final T rdpom;
  726.         private final T rdpomp;
  727.         private final T eps1;
  728.         private final T eps2;
  729.         private final T xim;
  730.         private final T ommD;
  731.         private final T rdl;
  732.         private final T aMD;

  734.         private final T kh;
  735.         private final T kl;

  737.         private final T ax1;
  738.         private final T ay1;
  739.         private final T as1;
  740.         private final T ac2;
  741.         private final T axy3;
  742.         private final T as3;
  743.         private final T ac4;
  744.         private final T as5;
  745.         private final T ac6;

  747.         private final T ex1;
  748.         private final T exx2;
  749.         private final T exy2;
  750.         private final T ex3;
  751.         private final T ex4;

  753.         private final T ey1;
  754.         private final T eyx2;
  755.         private final T eyy2;
  756.         private final T ey3;
  757.         private final T ey4;

  759.         private final T rx1;
  760.         private final T ry1;
  761.         private final T r2;
  762.         private final T r3;
  763.         private final T rl;

  765.         private final T iy1;
  766.         private final T ix1;
  767.         private final T i2;
  768.         private final T i3;
  769.         private final T ih;

  771.         private final T lx1;
  772.         private final T ly1;
  773.         private final T l2;
  774.         private final T l3;
  775.         private final T ll;

  777.         // CHECKSTYLE: resume JavadocVariable check

  779.         /** Create a model for specified mean FieldOrbit.
  780.          * @param mean mean FieldOrbit
  781.          * @param mass constant mass
  782.          * @param referenceRadius reference radius of the central body attraction model (m)
  783.          * @param mu central attraction coefficient (m³/s²)
  784.          * @param ck0 un-normalized zonal coefficients
  785.          */
  786.         FieldEHModel(final FieldCircularOrbit<T> mean, final T mass,
  787.                      final double referenceRadius, final T mu, final double[] ck0) {

  789.             this.mean            = mean;
  790.             this.mass            = mass;
  791.             final T zero = mass.getField().getZero();
  792.             final T one  = mass.getField().getOne();
  793.             // preliminary processing
  794.             T q =  zero.newInstance(referenceRadius).divide(mean.getA());
  795.             T ql = q.multiply(q);
  796.             final T g2 = ql.multiply(ck0[2]);
  797.             ql = ql.multiply(q);
  798.             final T g3 = ql.multiply(ck0[3]);
  799.             ql = ql.multiply(q);
  800.             final T g4 = ql.multiply(ck0[4]);
  801.             ql = ql.multiply(q);
  802.             final T g5 = ql.multiply(ck0[5]);
  803.             ql = ql.multiply(q);
  804.             final T g6 = ql.multiply(ck0[6]);

  806.             final FieldSinCos<T> sc = FastMath.sinCos(mean.getI());
  807.             final T cosI1 = sc.cos();
  808.             final T sinI1 = sc.sin();
  809.             final T sinI2 = sinI1.multiply(sinI1);
  810.             final T sinI4 = sinI2.multiply(sinI2);
  811.             final T sinI6 = sinI2.multiply(sinI4);

  813.             if (sinI2.getReal() < 1.0e-10) {
  814.                 throw new OrekitException(OrekitMessages.ALMOST_EQUATORIAL_ORBIT,
  815.                                                FastMath.toDegrees(mean.getI().getReal()));
  816.             }

  818.             if (FastMath.abs(sinI2.getReal() - 4.0 / 5.0) < 1.0e-3) {
  819.                 throw new OrekitException(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT,
  820.                                                FastMath.toDegrees(mean.getI().getReal()));
  821.             }

  823.             if (mean.getE().getReal() > 0.1) {
  824.                 // if 0.005 < e < 0.1 no error is triggered, but accuracy is poor
  825.                 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
  826.                                                mean.getE());
  827.             }

  829.             xnotDot = mu.divide(mean.getA()).sqrt().divide(mean.getA());

  831.             rdpom = g2.multiply(-0.75).multiply(sinI2.multiply(-5.0).add(4.0));
  832.             rdpomp = g4.multiply(7.5).multiply(sinI2.multiply(-31.0 / 8.0).add(1.0).add( sinI4.multiply(49.0 / 16.0))).subtract(
  833.                     g6.multiply(13.125).multiply(one.subtract(sinI2.multiply(8.0)).add(sinI4.multiply(129.0 / 8.0)).subtract(sinI6.multiply(297.0 / 32.0)) ));


  836.             q = zero.newInstance(3.0).divide(rdpom.multiply(32.0));
  837.             eps1 = q.multiply(g4).multiply(sinI2).multiply(sinI2.multiply(-35.0).add(30.0)).subtract(
  838.                    q.multiply(175.0).multiply(g6).multiply(sinI2).multiply(sinI2.multiply(-3.0).add(sinI4.multiply(2.0625)).add(1.0)));
  839.             q = sinI1.multiply(3.0).divide(rdpom.multiply(8.0));
  840.             eps2 = q.multiply(g3).multiply(sinI2.multiply(-5.0).add(4.0)).subtract(q.multiply(g5).multiply(sinI2.multiply(-35.0).add(sinI4.multiply(26.25)).add(10.0)));

  842.             xim = mean.getI();
  843.             ommD = cosI1.multiply(g2.multiply(1.50).subtract(g2.multiply(2.25).multiply(g2).multiply(sinI2.multiply(-19.0 / 6.0).add(2.5))).add(
  844.                             g4.multiply(0.9375).multiply(sinI2.multiply(7.0).subtract(4.0))).add(
  845.                             g6.multiply(3.28125).multiply(sinI2.multiply(-9.0).add(2.0).add(sinI4.multiply(8.25)))));

  847.             rdl = g2.multiply(-1.50).multiply(sinI2.multiply(-4.0).add(3.0)).add(1.0);
  848.             aMD = rdl.add(
  849.                     g2.multiply(2.25).multiply(g2.multiply(sinI2.multiply(-263.0 / 12.0 ).add(9.0).add(sinI4.multiply(341.0 / 24.0))))).add(
  850.                     g4.multiply(15.0 / 16.0).multiply(sinI2.multiply(-31.0).add(8.0).add(sinI4.multiply(24.5)))).add(
  851.                     g6.multiply(105.0 / 32.0).multiply(sinI2.multiply(25.0).add(-10.0 / 3.0).subtract(sinI4.multiply(48.75)).add(sinI6.multiply(27.5))));

  853.             final T qq   = g2.divide(rdl).multiply(-1.5);
  854.             final T qA   = g2.multiply(0.75).multiply(g2).multiply(sinI2);
  855.             final T qB   = g4.multiply(0.25).multiply(sinI2);
  856.             final T qC   = g6.multiply(105.0 / 16.0).multiply(sinI2);
  857.             final T qD   = g3.multiply(-0.75).multiply(sinI1);
  858.             final T qE   = g5.multiply(3.75).multiply(sinI1);
  859.             kh = zero.newInstance(0.375).divide(rdpom);
  860.             kl = kh.divide(sinI1);

  862.             ax1 = qq.multiply(sinI2.multiply(-3.5).add(2.0));
  863.             ay1 = qq.multiply(sinI2.multiply(-2.5).add(2.0));
  864.             as1 = qD.multiply(sinI2.multiply(-5.0).add(4.0)).add(
  865.                   qE.multiply(sinI4.multiply(2.625).add(sinI2.multiply(-3.5)).add(1.0)));
  866.             ac2 = qq.multiply(sinI2).add(
  867.                   qA.multiply(7.0).multiply(sinI2.multiply(-3.0).add(2.0))).add(
  868.                   qB.multiply(sinI2.multiply(-17.5).add(15.0))).add(
  869.                   qC.multiply(sinI2.multiply(3.0).subtract(1.0).subtract(sinI4.multiply(33.0 / 16.0))));
  870.             axy3 = qq.multiply(3.5).multiply(sinI2);
  871.             as3 = qD.multiply(5.0 / 3.0).multiply(sinI2).add(
  872.                   qE.multiply(7.0 / 6.0).multiply(sinI2).multiply(sinI2.multiply(-1.125).add(1)));
  873.             ac4 = qA.multiply(sinI2).add(
  874.                   qB.multiply(4.375).multiply(sinI2)).add(
  875.                   qC.multiply(0.75).multiply(sinI4.multiply(1.1).subtract(sinI2)));

  877.             as5 = qE.multiply(21.0 / 80.0).multiply(sinI4);

  879.             ac6 = qC.multiply(-11.0 / 80.0).multiply(sinI4);

  881.             ex1 = qq.multiply(sinI2.multiply(-1.25).add(1.0));
  882.             exx2 = qq.multiply(0.5).multiply(sinI2.multiply(-5.0).add(3.0));
  883.             exy2 = qq.multiply(sinI2.multiply(-1.5).add(2.0));
  884.             ex3 = qq.multiply(7.0 / 12.0).multiply(sinI2);
  885.             ex4 = qq.multiply(17.0 / 8.0).multiply(sinI2);

  887.             ey1 = qq.multiply(sinI2.multiply(-1.75).add(1.0));
  888.             eyx2 = qq.multiply(sinI2.multiply(-3.0).add(1.0));
  889.             eyy2 = qq.multiply(sinI2.multiply(2.0).subtract(1.5));
  890.             ey3 = qq.multiply(7.0 / 12.0).multiply(sinI2);
  891.             ey4 = qq.multiply(17.0 / 8.0).multiply(sinI2);

  893.             q  = cosI1.multiply(qq).negate();
  894.             rx1 = q.multiply(3.5);
  895.             ry1 = q.multiply(-2.5);
  896.             r2 = q.multiply(-0.5);
  897.             r3 =  q.multiply(7.0 / 6.0);
  898.             rl = g3 .multiply( cosI1).multiply(sinI2.multiply(-15.0).add(4.0)).subtract(
  899.                  g5.multiply(2.5).multiply(cosI1).multiply(sinI2.multiply(-42.0).add(4.0).add(sinI4.multiply(52.5))));

  901.             q = qq.multiply(0.5).multiply(sinI1).multiply(cosI1);
  902.             iy1 =  q;
  903.             ix1 = q.negate();
  904.             i2 =  q;
  905.             i3 =  q.multiply(7.0 / 3.0);
  906.             ih = g3.negate().multiply(cosI1).multiply(sinI2.multiply(-5.0).add(4)).add(
  907.                  g5.multiply(2.5).multiply(cosI1).multiply(sinI2.multiply(-14.0).add(4.0).add(sinI4.multiply(10.5))));
  908.             lx1 = qq.multiply(sinI2.multiply(-77.0 / 8.0).add(7.0));
  909.             ly1 = qq.multiply(sinI2.multiply(55.0 / 8.0).subtract(7.50));
  910.             l2 = qq.multiply(sinI2.multiply(1.25).subtract(0.5));
  911.             l3 = qq.multiply(sinI2.multiply(77.0 / 24.0).subtract(7.0 / 6.0));
  912.             ll = g3.multiply(sinI2.multiply(53.0).subtract(4.0).add(sinI4.multiply(-57.5))).add(
  913.                  g5.multiply(2.5).multiply(sinI2.multiply(-96.0).add(4.0).add(sinI4.multiply(269.5).subtract(sinI6.multiply(183.75)))));

  915.         }

  917.         /** Extrapolate a FieldOrbit up to a specific target date.
  918.          * @param date target date for the FieldOrbit
  919.          * @return propagated parameters
  920.          */
  921.         public FieldUnivariateDerivative2<T>[] propagateParameters(final FieldAbsoluteDate<T> date) {
  922.             final Field<T> field = date.durationFrom(mean.getDate()).getField();
  923.             final T one = field.getOne();
  924.             final T zero = field.getZero();
  925.             // Keplerian evolution
  926.             final FieldUnivariateDerivative2<T> dt = new FieldUnivariateDerivative2<>(date.durationFrom(mean.getDate()), one, zero);
  927.             final FieldUnivariateDerivative2<T> xnot = dt.multiply(xnotDot);

  929.             // secular effects

  931.             // eccentricity
  932.             final FieldUnivariateDerivative2<T> x   = xnot.multiply(rdpom.add(rdpomp));
  933.             final FieldUnivariateDerivative2<T> cx  = x.cos();
  934.             final FieldUnivariateDerivative2<T> sx  = x.sin();
  935.             final FieldUnivariateDerivative2<T> exm = cx.multiply(mean.getCircularEx()).
  936.                                             add(sx.multiply(eps2.subtract(one.subtract(eps1).multiply(mean.getCircularEy()))));
  937.             final FieldUnivariateDerivative2<T> eym = sx.multiply(eps1.add(1.0).multiply(mean.getCircularEx())).
  938.                                             add(cx.multiply(mean.getCircularEy().subtract(eps2))).
  939.                                             add(eps2);
  940.             // no secular effect on inclination

  942.             // right ascension of ascending node
  943.             final FieldUnivariateDerivative2<T> omm =
  944.                            new FieldUnivariateDerivative2<>(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode().add(ommD.multiply(xnot.getValue())),
  945.                                                                                      one.getPi()),
  946.                                                             ommD.multiply(xnotDot),
  947.                                                             zero);
  948.             // latitude argument
  949.             final FieldUnivariateDerivative2<T> xlm =
  950.                             new FieldUnivariateDerivative2<>(MathUtils.normalizeAngle(mean.getAlphaM().add(aMD.multiply(xnot.getValue())),
  951.                                                                                       one.getPi()),
  952.                                                            aMD.multiply(xnotDot),
  953.                                                            zero);

  955.             // periodical terms
  956.             final FieldUnivariateDerivative2<T> cl1 = xlm.cos();
  957.             final FieldUnivariateDerivative2<T> sl1 = xlm.sin();
  958.             final FieldUnivariateDerivative2<T> cl2 = cl1.multiply(cl1).subtract(sl1.multiply(sl1));
  959.             final FieldUnivariateDerivative2<T> sl2 = cl1.multiply(sl1).add(sl1.multiply(cl1));
  960.             final FieldUnivariateDerivative2<T> cl3 = cl2.multiply(cl1).subtract(sl2.multiply(sl1));
  961.             final FieldUnivariateDerivative2<T> sl3 = cl2.multiply(sl1).add(sl2.multiply(cl1));
  962.             final FieldUnivariateDerivative2<T> cl4 = cl3.multiply(cl1).subtract(sl3.multiply(sl1));
  963.             final FieldUnivariateDerivative2<T> sl4 = cl3.multiply(sl1).add(sl3.multiply(cl1));
  964.             final FieldUnivariateDerivative2<T> cl5 = cl4.multiply(cl1).subtract(sl4.multiply(sl1));
  965.             final FieldUnivariateDerivative2<T> sl5 = cl4.multiply(sl1).add(sl4.multiply(cl1));
  966.             final FieldUnivariateDerivative2<T> cl6 = cl5.multiply(cl1).subtract(sl5.multiply(sl1));

  968.             final FieldUnivariateDerivative2<T> qh  = eym.subtract(eps2).multiply(kh);
  969.             final FieldUnivariateDerivative2<T> ql  = exm.multiply(kl);

  971.             final FieldUnivariateDerivative2<T> exmCl1 = exm.multiply(cl1);
  972.             final FieldUnivariateDerivative2<T> exmSl1 = exm.multiply(sl1);
  973.             final FieldUnivariateDerivative2<T> eymCl1 = eym.multiply(cl1);
  974.             final FieldUnivariateDerivative2<T> eymSl1 = eym.multiply(sl1);
  975.             final FieldUnivariateDerivative2<T> exmCl2 = exm.multiply(cl2);
  976.             final FieldUnivariateDerivative2<T> exmSl2 = exm.multiply(sl2);
  977.             final FieldUnivariateDerivative2<T> eymCl2 = eym.multiply(cl2);
  978.             final FieldUnivariateDerivative2<T> eymSl2 = eym.multiply(sl2);
  979.             final FieldUnivariateDerivative2<T> exmCl3 = exm.multiply(cl3);
  980.             final FieldUnivariateDerivative2<T> exmSl3 = exm.multiply(sl3);
  981.             final FieldUnivariateDerivative2<T> eymCl3 = eym.multiply(cl3);
  982.             final FieldUnivariateDerivative2<T> eymSl3 = eym.multiply(sl3);
  983.             final FieldUnivariateDerivative2<T> exmCl4 = exm.multiply(cl4);
  984.             final FieldUnivariateDerivative2<T> exmSl4 = exm.multiply(sl4);
  985.             final FieldUnivariateDerivative2<T> eymCl4 = eym.multiply(cl4);
  986.             final FieldUnivariateDerivative2<T> eymSl4 = eym.multiply(sl4);

  988.             // semi major axis
  989.             final FieldUnivariateDerivative2<T> rda = exmCl1.multiply(ax1).
  990.                                             add(eymSl1.multiply(ay1)).
  991.                                             add(sl1.multiply(as1)).
  992.                                             add(cl2.multiply(ac2)).
  993.                                             add(exmCl3.add(eymSl3).multiply(axy3)).
  994.                                             add(sl3.multiply(as3)).
  995.                                             add(cl4.multiply(ac4)).
  996.                                             add(sl5.multiply(as5)).
  997.                                             add(cl6.multiply(ac6));

  999.             // eccentricity
  1000.             final FieldUnivariateDerivative2<T> rdex = cl1.multiply(ex1).
  1001.                                              add(exmCl2.multiply(exx2)).
  1002.                                              add(eymSl2.multiply(exy2)).
  1003.                                              add(cl3.multiply(ex3)).
  1004.                                              add(exmCl4.add(eymSl4).multiply(ex4));
  1005.             final FieldUnivariateDerivative2<T> rdey = sl1.multiply(ey1).
  1006.                                              add(exmSl2.multiply(eyx2)).
  1007.                                              add(eymCl2.multiply(eyy2)).
  1008.                                              add(sl3.multiply(ey3)).
  1009.                                              add(exmSl4.subtract(eymCl4).multiply(ey4));

  1011.             // ascending node
  1012.             final FieldUnivariateDerivative2<T> rdom = exmSl1.multiply(rx1).
  1013.                                              add(eymCl1.multiply(ry1)).
  1014.                                              add(sl2.multiply(r2)).
  1015.                                              add(eymCl3.subtract(exmSl3).multiply(r3)).
  1016.                                              add(ql.multiply(rl));

  1018.             // inclination
  1019.             final FieldUnivariateDerivative2<T> rdxi = eymSl1.multiply(iy1).
  1020.                                              add(exmCl1.multiply(ix1)).
  1021.                                              add(cl2.multiply(i2)).
  1022.                                              add(exmCl3.add(eymSl3).multiply(i3)).
  1023.                                              add(qh.multiply(ih));

  1025.             // latitude argument
  1026.             final FieldUnivariateDerivative2<T> rdxl = exmSl1.multiply(lx1).
  1027.                                              add(eymCl1.multiply(ly1)).
  1028.                                              add(sl2.multiply(l2)).
  1029.                                              add(exmSl3.subtract(eymCl3).multiply(l3)).
  1030.                                              add(ql.multiply(ll));
  1031.             // osculating parameters
  1032.             final FieldUnivariateDerivative2<T>[] FTD = MathArrays.buildArray(rdxl.getField(), 6);

  1034.             FTD[0] = rda.add(1.0).multiply(mean.getA());
  1035.             FTD[1] = rdex.add(exm);
  1036.             FTD[2] = rdey.add(eym);
  1037.             FTD[3] = rdxi.add(xim);
  1038.             FTD[4] = rdom.add(omm);
  1039.             FTD[5] = rdxl.add(xlm);
  1040.             return FTD;

  1042.         }

  1044.     }

  1046.     /** Convert circular parameters <em>with derivatives</em> to Cartesian coordinates.
  1047.      * @param date date of the parameters
  1048.      * @param parameters circular parameters (a, ex, ey, i, raan, alphaM)
  1049.      * @return Cartesian coordinates consistent with values and derivatives
  1050.      */
  1051.     private TimeStampedFieldPVCoordinates<T> toCartesian(final FieldAbsoluteDate<T> date, final FieldUnivariateDerivative2<T>[] parameters) {

  1053.         // evaluate coordinates in the FieldOrbit canonical reference frame
  1054.         final FieldUnivariateDerivative2<T> cosOmega = parameters[4].cos();
  1055.         final FieldUnivariateDerivative2<T> sinOmega = parameters[4].sin();
  1056.         final FieldUnivariateDerivative2<T> cosI     = parameters[3].cos();
  1057.         final FieldUnivariateDerivative2<T> sinI     = parameters[3].sin();
  1058.         final FieldUnivariateDerivative2<T> alphaE   = meanToEccentric(parameters[5], parameters[1], parameters[2]);
  1059.         final FieldUnivariateDerivative2<T> cosAE    = alphaE.cos();
  1060.         final FieldUnivariateDerivative2<T> sinAE    = alphaE.sin();
  1061.         final FieldUnivariateDerivative2<T> ex2      = parameters[1].square();
  1062.         final FieldUnivariateDerivative2<T> ey2      = parameters[2].square();
  1063.         final FieldUnivariateDerivative2<T> exy      = parameters[1].multiply(parameters[2]);
  1064.         final FieldUnivariateDerivative2<T> q        = ex2.add(ey2).subtract(1).negate().sqrt();
  1065.         final FieldUnivariateDerivative2<T> beta     = q.add(1).reciprocal();
  1066.         final FieldUnivariateDerivative2<T> bx2      = beta.multiply(ex2);
  1067.         final FieldUnivariateDerivative2<T> by2      = beta.multiply(ey2);
  1068.         final FieldUnivariateDerivative2<T> bxy      = beta.multiply(exy);
  1069.         final FieldUnivariateDerivative2<T> u        = bxy.multiply(sinAE).subtract(parameters[1].add(by2.subtract(1).multiply(cosAE)));
  1070.         final FieldUnivariateDerivative2<T> v        = bxy.multiply(cosAE).subtract(parameters[2].add(bx2.subtract(1).multiply(sinAE)));
  1071.         final FieldUnivariateDerivative2<T> x        = parameters[0].multiply(u);
  1072.         final FieldUnivariateDerivative2<T> y        = parameters[0].multiply(v);

  1074.         // canonical FieldOrbit reference frame
  1075.         final FieldVector3D<FieldUnivariateDerivative2<T>> p =
  1076.                 new FieldVector3D<>(x.multiply(cosOmega).subtract(y.multiply(cosI.multiply(sinOmega))),
  1077.                                     x.multiply(sinOmega).add(y.multiply(cosI.multiply(cosOmega))),
  1078.                                     y.multiply(sinI));

  1080.         // dispatch derivatives
  1081.         final FieldVector3D<T> p0 = new FieldVector3D<>(p.getX().getValue(),
  1082.                                                         p.getY().getValue(),
  1083.                                                         p.getZ().getValue());
  1084.         final FieldVector3D<T> p1 = new FieldVector3D<>(p.getX().getFirstDerivative(),
  1085.                                                         p.getY().getFirstDerivative(),
  1086.                                                         p.getZ().getFirstDerivative());
  1087.         final FieldVector3D<T> p2 = new FieldVector3D<>(p.getX().getSecondDerivative(),
  1088.                                                         p.getY().getSecondDerivative(),
  1089.                                                         p.getZ().getSecondDerivative());
  1090.         return new TimeStampedFieldPVCoordinates<>(date, p0, p1, p2);

  1092.     }

  1094.     /** Computes the eccentric latitude argument from the mean latitude argument.
  1095.      * @param alphaM = M + Ω mean latitude argument (rad)
  1096.      * @param ex e cos(Ω), first component of circular eccentricity vector
  1097.      * @param ey e sin(Ω), second component of circular eccentricity vector
  1098.      * @return the eccentric latitude argument.
  1099.      */
  1100.     private FieldUnivariateDerivative2<T> meanToEccentric(final FieldUnivariateDerivative2<T> alphaM,
  1101.                                                 final FieldUnivariateDerivative2<T> ex,
  1102.                                                 final FieldUnivariateDerivative2<T> ey) {
  1103.         // Generalization of Kepler equation to circular parameters
  1104.         // with alphaE = PA + E and
  1105.         //      alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
  1106.         FieldUnivariateDerivative2<T> alphaE        = alphaM;
  1107.         FieldUnivariateDerivative2<T> shift         = alphaM.getField().getZero();
  1108.         FieldUnivariateDerivative2<T> alphaEMalphaM = alphaM.getField().getZero();
  1109.         FieldUnivariateDerivative2<T> cosAlphaE     = alphaE.cos();
  1110.         FieldUnivariateDerivative2<T> sinAlphaE     = alphaE.sin();
  1111.         int                 iter          = 0;
  1112.         do {
  1113.             final FieldUnivariateDerivative2<T> f2 = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
  1114.             final FieldUnivariateDerivative2<T> f1 = alphaM.getField().getOne().subtract(ex.multiply(cosAlphaE)).subtract(ey.multiply(sinAlphaE));
  1115.             final FieldUnivariateDerivative2<T> f0 = alphaEMalphaM.subtract(f2);

  1117.             final FieldUnivariateDerivative2<T> f12 = f1.multiply(2);
  1118.             shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));

  1120.             alphaEMalphaM  = alphaEMalphaM.subtract(shift);
  1121.             alphaE         = alphaM.add(alphaEMalphaM);
  1122.             cosAlphaE      = alphaE.cos();
  1123.             sinAlphaE      = alphaE.sin();

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

  1127.         return alphaE;

  1129.     }

  1131.     /** {@inheritDoc} */
  1132.     @Override
  1133.     protected T getMass(final FieldAbsoluteDate<T> date) {
  1134.         return models.get(date).mass;
  1135.     }

  1137.     /** {@inheritDoc} */
  1138.     @Override
  1139.     public List<ParameterDriver> getParametersDrivers() {
  1140.         // Eckstein Hechler propagation model does not have parameter drivers.
  1141.         return Collections.emptyList();
  1142.     }

  1144. }
Created with JaCoCo 0.8.12.202403310830