FieldEcksteinHechlerPropagator.java

  1. /* Copyright 2002-2020 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.Field;
  20. import org.hipparchus.RealFieldElement;
  21. import org.hipparchus.analysis.differentiation.FieldUnivariateDerivative2;
  22. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  23. import org.hipparchus.util.FastMath;
  24. import org.hipparchus.util.FieldSinCos;
  25. import org.hipparchus.util.MathArrays;
  26. import org.hipparchus.util.MathUtils;
  27. import org.orekit.annotation.DefaultDataContext;
  28. import org.orekit.attitudes.AttitudeProvider;
  29. import org.orekit.data.DataContext;
  30. import org.orekit.errors.OrekitException;
  31. import org.orekit.errors.OrekitMessages;
  32. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  33. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
  34. import org.orekit.orbits.FieldCartesianOrbit;
  35. import org.orekit.orbits.FieldCircularOrbit;
  36. import org.orekit.orbits.FieldOrbit;
  37. import org.orekit.orbits.OrbitType;
  38. import org.orekit.orbits.PositionAngle;
  39. import org.orekit.propagation.FieldSpacecraftState;
  40. import org.orekit.propagation.PropagationType;
  41. import org.orekit.propagation.Propagator;
  42. import org.orekit.time.FieldAbsoluteDate;
  43. import org.orekit.utils.FieldTimeSpanMap;
  44. import org.orekit.utils.TimeStampedFieldPVCoordinates;

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

  57.     /** Initial Eckstein-Hechler model. */
  58.     private FieldEHModel<T> initialModel;

  60.     /** All models. */
  61.     private transient FieldTimeSpanMap<FieldEHModel<T>, T> models;

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

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

  69.     /** Un-normalized zonal coefficients. */
  70.     private double[] ck0;

  72.     /** Build a propagator from FieldOrbit and potential provider.
  73.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  74.      *
  75.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  76.      *
  77.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  78.      *
  79.      * @param initialOrbit initial FieldOrbit
  80.      * @param provider for un-normalized zonal coefficients
  81.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, AttitudeProvider,
  82.      * UnnormalizedSphericalHarmonicsProvider)
  83.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, UnnormalizedSphericalHarmonicsProvider,
  84.      * PropagationType)
  85.      */
  86.     @DefaultDataContext
  87.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  88.                                           final UnnormalizedSphericalHarmonicsProvider provider) {
  89.         this(initialOrbit, Propagator.getDefaultLaw(DataContext.getDefault().getFrames()),
  90.              initialOrbit.getA().getField().getZero().add(DEFAULT_MASS), provider,
  91.              provider.onDate(initialOrbit.getDate().toAbsoluteDate()));
  92.     }

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

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

  156.     /** Build a propagator from FieldOrbit, mass and potential provider.
  157.      * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
  158.      *
  159.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  160.      *
  161.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  162.      *
  163.      * @param initialOrbit initial FieldOrbit
  164.      * @param mass spacecraft mass
  165.      * @param provider for un-normalized zonal coefficients
  166.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, AttitudeProvider,
  167.      * RealFieldElement, UnnormalizedSphericalHarmonicsProvider)
  168.      */
  169.     @DefaultDataContext
  170.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit, final T mass,
  171.                                           final UnnormalizedSphericalHarmonicsProvider provider) {
  172.         this(initialOrbit, Propagator.getDefaultLaw(DataContext.getDefault().getFrames()),
  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>This constructor uses the {@link DataContext#getDefault() default data context}.
  189.      *
  190.      * <p>Using this constructor, an initial osculating orbit is considered.</p>
  191.      *
  192.      * @param initialOrbit initial FieldOrbit
  193.      * @param mass spacecraft mass
  194.      * @param referenceRadius reference radius of the Earth for the potential model (m)
  195.      * @param mu central attraction coefficient (m³/s²)
  196.      * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
  197.      * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
  198.      * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
  199.      * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
  200.      * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
  201.      * @see #FieldEcksteinHechlerPropagator(FieldOrbit, AttitudeProvider,
  202.      * RealFieldElement, double, RealFieldElement, double, double, double, double, double)
  203.      */
  204.     @DefaultDataContext
  205.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit, final T mass,
  206.                                           final double referenceRadius, final T mu,
  207.                                           final double c20, final double c30, final double c40,
  208.                                           final double c50, final double c60) {
  209.         this(initialOrbit, Propagator.getDefaultLaw(DataContext.getDefault().getFrames()),
  210.                 mass, referenceRadius, mu, c20, c30, c40, c50, c60);
  211.     }

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

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

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

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

  310.     /** Build a propagator from orbit and potential provider.
  311.      * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
  312.      *
  313.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  314.      *
  315.      * <p>Using this constructor, it is possible to define the initial orbit as
  316.      * a mean Eckstein-Hechler orbit or an osculating one.</p>
  317.      *
  318.      * @param initialOrbit initial orbit
  319.      * @param provider for un-normalized zonal coefficients
  320.      * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
  321.      * @since 10.2
  322.      */
  323.     @DefaultDataContext
  324.     public FieldEcksteinHechlerPropagator(final FieldOrbit<T> initialOrbit,
  325.                                           final UnnormalizedSphericalHarmonicsProvider provider,
  326.                                           final PropagationType initialType) {
  327.         this(initialOrbit, Propagator.getDefaultLaw(DataContext.getDefault().getFrames()),
  328.              initialOrbit.getA().getField().getZero().add(DEFAULT_MASS), provider,
  329.              provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType);
  330.     }

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

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

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

  412.         super(mass.getField(), attitudeProv);
  413.         try {

  415.             // store model coefficients
  416.             this.referenceRadius = referenceRadius;
  417.             this.mu  = mu;
  418.             this.ck0 = new double[] {
  419.                 0.0, 0.0, c20, c30, c40, c50, c60
  420.             };

  422.             // compute mean parameters if needed
  423.             // transform into circular adapted parameters used by the Eckstein-Hechler model
  424.             resetInitialState(new FieldSpacecraftState<>(initialOrbit,
  425.                                                          attitudeProv.getAttitude(initialOrbit,
  426.                                                                                   initialOrbit.getDate(),
  427.                                                                                   initialOrbit.getFrame()),
  428.                                                          mass),
  429.                               initialType);

  431.         } catch (OrekitException oe) {
  432.             throw new OrekitException(oe);
  433.         }
  434.     }

  436.     /** {@inheritDoc}
  437.      * <p>The new initial state to consider
  438.      * must be defined with an osculating orbit.</p>
  439.      * @see #resetInitialState(FieldSpacecraftState, PropagationType)
  440.      */
  441.     public void resetInitialState(final FieldSpacecraftState<T> state) {
  442.         resetInitialState(state, PropagationType.OSCULATING);
  443.     }

  445.     /** Reset the propagator initial state.
  446.      * @param state new initial state to consider
  447.      * @param stateType mean Eckstein-Hechler orbit or osculating orbit
  448.      * @since 10.2
  449.      */
  450.     public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType) {
  451.         super.resetInitialState(state);
  452.         final FieldCircularOrbit<T> circular = (FieldCircularOrbit<T>) OrbitType.CIRCULAR.convertType(state.getOrbit());
  453.         this.initialModel = (stateType == PropagationType.MEAN) ?
  454.                              new FieldEHModel<>(circular, state.getMass(), referenceRadius, mu, ck0) :
  455.                              computeMeanParameters(circular, state.getMass());
  456.         this.models = new FieldTimeSpanMap<FieldEHModel<T>, T>(initialModel, state.getA().getField());
  457.     }

  459.     /** {@inheritDoc} */
  460.     protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward) {
  461.         final FieldEHModel<T> newModel = computeMeanParameters((FieldCircularOrbit<T>) OrbitType.CIRCULAR.convertType(state.getOrbit()),
  462.                                                        state.getMass());
  463.         if (forward) {
  464.             models.addValidAfter(newModel, state.getDate());
  465.         } else {
  466.             models.addValidBefore(newModel, state.getDate());
  467.         }
  468.         stateChanged(state);
  469.     }

  471.     /** Compute mean parameters according to the Eckstein-Hechler analytical model.
  472.      * @param osculating osculating FieldOrbit
  473.      * @param mass constant mass
  474.      * @return Eckstein-Hechler mean model
  475.      */
  476.     private FieldEHModel<T> computeMeanParameters(final FieldCircularOrbit<T> osculating, final T mass) {

  478.         // sanity check
  479.         if (osculating.getA().getReal() < referenceRadius) {
  480.             throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
  481.                                            osculating.getA());
  482.         }
  483.         final Field<T> field = mass.getField();
  484.         final T one = field.getOne();
  485.         final T zero = field.getZero();
  486.         // rough initialization of the mean parameters
  487.         FieldEHModel<T> current = new FieldEHModel<>(osculating, mass, referenceRadius, mu, ck0);
  488.         // threshold for each parameter
  489.         final T epsilon         = one .multiply(1.0e-13);
  490.         final T thresholdA      = epsilon.multiply(current.mean.getA().abs().add(1.0));
  491.         final T thresholdE      = epsilon.multiply(current.mean.getE().add(1.0));
  492.         final T thresholdAngles = epsilon.multiply(FastMath.PI);


  495.         int i = 0;
  496.         while (i++ < 100) {

  498.             // recompute the osculating parameters from the current mean parameters
  499.             final FieldUnivariateDerivative2<T>[] parameters = current.propagateParameters(current.mean.getDate());
  500.             // adapted parameters residuals
  501.             final T deltaA      = osculating.getA()         .subtract(parameters[0].getValue());
  502.             final T deltaEx     = osculating.getCircularEx().subtract(parameters[1].getValue());
  503.             final T deltaEy     = osculating.getCircularEy().subtract(parameters[2].getValue());
  504.             final T deltaI      = osculating.getI()         .subtract(parameters[3].getValue());
  505.             final T deltaRAAN   = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode().subtract(
  506.                                                                 parameters[4].getValue()),
  507.                                                                 zero);
  508.             final T deltaAlphaM = MathUtils.normalizeAngle(osculating.getAlphaM().subtract(parameters[5].getValue()), zero);
  509.             // update mean parameters
  510.             current = new FieldEHModel<>(new FieldCircularOrbit<>(current.mean.getA().add(deltaA),
  511.                                                                   current.mean.getCircularEx().add( deltaEx),
  512.                                                                   current.mean.getCircularEy().add( deltaEy),
  513.                                                                   current.mean.getI()         .add( deltaI ),
  514.                                                                   current.mean.getRightAscensionOfAscendingNode().add(deltaRAAN),
  515.                                                                   current.mean.getAlphaM().add(deltaAlphaM),
  516.                                                                   PositionAngle.MEAN,
  517.                                                                   current.mean.getFrame(),
  518.                                                                   current.mean.getDate(), mu),
  519.                                   mass, referenceRadius, mu, ck0);
  520.             // check convergence
  521.             if ((FastMath.abs(deltaA.getReal())      < thresholdA.getReal()) &&
  522.                 (FastMath.abs(deltaEx.getReal())     < thresholdE.getReal()) &&
  523.                 (FastMath.abs(deltaEy.getReal())     < thresholdE.getReal()) &&
  524.                 (FastMath.abs(deltaI.getReal())      < thresholdAngles.getReal()) &&
  525.                 (FastMath.abs(deltaRAAN.getReal())   < thresholdAngles.getReal()) &&
  526.                 (FastMath.abs(deltaAlphaM.getReal()) < thresholdAngles.getReal())) {
  527.                 return current;
  528.             }

  530.         }

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

  534.     }

  536.     /** {@inheritDoc} */
  537.     public FieldCartesianOrbit<T> propagateOrbit(final FieldAbsoluteDate<T> date) {
  538.         // compute Cartesian parameters, taking derivatives into account
  539.         // to make sure velocity and acceleration are consistent
  540.         final FieldEHModel<T> current = models.get(date);
  541.         return new FieldCartesianOrbit<>(toCartesian(date, current.propagateParameters(date)),
  542.                                          current.mean.getFrame(), mu);
  543.     }

  545.     /** Local class for Eckstein-Hechler model, with fixed mean parameters. */
  546.     private static class FieldEHModel<T extends RealFieldElement<T>> {

  548.         /** Mean FieldOrbit. */
  549.         private final FieldCircularOrbit<T> mean;

  551.         /** Constant mass. */
  552.         private final T mass;
  553.         // CHECKSTYLE: stop JavadocVariable check

  555.         // preprocessed values
  556.         private final T xnotDot;
  557.         private final T rdpom;
  558.         private final T rdpomp;
  559.         private final T eps1;
  560.         private final T eps2;
  561.         private final T xim;
  562.         private final T ommD;
  563.         private final T rdl;
  564.         private final T aMD;

  566.         private final T kh;
  567.         private final T kl;

  569.         private final T ax1;
  570.         private final T ay1;
  571.         private final T as1;
  572.         private final T ac2;
  573.         private final T axy3;
  574.         private final T as3;
  575.         private final T ac4;
  576.         private final T as5;
  577.         private final T ac6;

  579.         private final T ex1;
  580.         private final T exx2;
  581.         private final T exy2;
  582.         private final T ex3;
  583.         private final T ex4;

  585.         private final T ey1;
  586.         private final T eyx2;
  587.         private final T eyy2;
  588.         private final T ey3;
  589.         private final T ey4;

  591.         private final T rx1;
  592.         private final T ry1;
  593.         private final T r2;
  594.         private final T r3;
  595.         private final T rl;

  597.         private final T iy1;
  598.         private final T ix1;
  599.         private final T i2;
  600.         private final T i3;
  601.         private final T ih;

  603.         private final T lx1;
  604.         private final T ly1;
  605.         private final T l2;
  606.         private final T l3;
  607.         private final T ll;

  609.         // CHECKSTYLE: resume JavadocVariable check

  611.         /** Create a model for specified mean FieldOrbit.
  612.          * @param mean mean FieldOrbit
  613.          * @param mass constant mass
  614.          * @param referenceRadius reference radius of the central body attraction model (m)
  615.          * @param mu central attraction coefficient (m³/s²)
  616.          * @param ck0 un-normalized zonal coefficients
  617.          */
  618.         FieldEHModel(final FieldCircularOrbit<T> mean, final T mass,
  619.                      final double referenceRadius, final T mu, final double[] ck0) {

  621.             this.mean            = mean;
  622.             this.mass            = mass;
  623.             final T zero = mass.getField().getZero();
  624.             final T one  = mass.getField().getOne();
  625.             // preliminary processing
  626.             T q =  zero.add(referenceRadius).divide(mean.getA());
  627.             T ql = q.multiply(q);
  628.             final T g2 = ql.multiply(ck0[2]);
  629.             ql = ql.multiply(q);
  630.             final T g3 = ql.multiply(ck0[3]);
  631.             ql = ql.multiply(q);
  632.             final T g4 = ql.multiply(ck0[4]);
  633.             ql = ql.multiply(q);
  634.             final T g5 = ql.multiply(ck0[5]);
  635.             ql = ql.multiply(q);
  636.             final T g6 = ql.multiply(ck0[6]);

  638.             final FieldSinCos<T> sc = FastMath.sinCos(mean.getI());
  639.             final T cosI1 = sc.cos();
  640.             final T sinI1 = sc.sin();
  641.             final T sinI2 = sinI1.multiply(sinI1);
  642.             final T sinI4 = sinI2.multiply(sinI2);
  643.             final T sinI6 = sinI2.multiply(sinI4);

  645.             if (sinI2.getReal() < 1.0e-10) {
  646.                 throw new OrekitException(OrekitMessages.ALMOST_EQUATORIAL_ORBIT,
  647.                                                FastMath.toDegrees(mean.getI().getReal()));
  648.             }

  650.             if (FastMath.abs(sinI2.getReal() - 4.0 / 5.0) < 1.0e-3) {
  651.                 throw new OrekitException(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT,
  652.                                                FastMath.toDegrees(mean.getI().getReal()));
  653.             }

  655.             if (mean.getE().getReal() > 0.1) {
  656.                 // if 0.005 < e < 0.1 no error is triggered, but accuracy is poor
  657.                 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
  658.                                                mean.getE());
  659.             }

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

  663.             rdpom = g2.multiply(-0.75).multiply(sinI2.multiply(-5.0).add(4.0));
  664.             rdpomp = g4.multiply(7.5).multiply(sinI2.multiply(-31.0 / 8.0).add(1.0).add( sinI4.multiply(49.0 / 16.0))).subtract(
  665.                     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)) ));


  668.             q = zero.add(3.0).divide(rdpom.multiply(32.0));
  669.             eps1 = q.multiply(g4).multiply(sinI2).multiply(sinI2.multiply(-35.0).add(30.0)).subtract(
  670.                    q.multiply(175.0).multiply(g6).multiply(sinI2).multiply(sinI2.multiply(-3.0).add(sinI4.multiply(2.0625)).add(1.0)));
  671.             q = sinI1.multiply(3.0).divide(rdpom.multiply(8.0));
  672.             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)));

  674.             xim = mean.getI();
  675.             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(
  676.                             g4.multiply(0.9375).multiply(sinI2.multiply(7.0).subtract(4.0))).add(
  677.                             g6.multiply(3.28125).multiply(sinI2.multiply(-9.0).add(2.0).add(sinI4.multiply(8.25)))));

  679.             rdl = g2.multiply(-1.50).multiply(sinI2.multiply(-4.0).add(3.0)).add(1.0);
  680.             aMD = rdl.add(
  681.                     g2.multiply(2.25).multiply(g2.multiply(sinI2.multiply(-263.0 / 12.0 ).add(9.0).add(sinI4.multiply(341.0 / 24.0))))).add(
  682.                     g4.multiply(15.0 / 16.0).multiply(sinI2.multiply(-31.0).add(8.0).add(sinI4.multiply(24.5)))).add(
  683.                     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))));

  685.             final T qq   = g2.divide(rdl).multiply(-1.5);
  686.             final T qA   = g2.multiply(0.75).multiply(g2).multiply(sinI2);
  687.             final T qB   = g4.multiply(0.25).multiply(sinI2);
  688.             final T qC   = g6.multiply(105.0 / 16.0).multiply(sinI2);
  689.             final T qD   = g3.multiply(-0.75).multiply(sinI1);
  690.             final T qE   = g5.multiply(3.75).multiply(sinI1);
  691.             kh = zero.add(0.375).divide(rdpom);
  692.             kl = kh.divide(sinI1);

  694.             ax1 = qq.multiply(sinI2.multiply(-3.5).add(2.0));
  695.             ay1 = qq.multiply(sinI2.multiply(-2.5).add(2.0));
  696.             as1 = qD.multiply(sinI2.multiply(-5.0).add(4.0)).add(
  697.                   qE.multiply(sinI4.multiply(2.625).add(sinI2.multiply(-3.5)).add(1.0)));
  698.             ac2 = qq.multiply(sinI2).add(
  699.                   qA.multiply(7.0).multiply(sinI2.multiply(-3.0).add(2.0))).add(
  700.                   qB.multiply(sinI2.multiply(-17.5).add(15.0))).add(
  701.                   qC.multiply(sinI2.multiply(3.0).subtract(1.0).subtract(sinI4.multiply(33.0 / 16.0))));
  702.             axy3 = qq.multiply(3.5).multiply(sinI2);
  703.             as3 = qD.multiply(5.0 / 3.0).multiply(sinI2).add(
  704.                   qE.multiply(7.0 / 6.0).multiply(sinI2).multiply(sinI2.multiply(-1.125).add(1)));
  705.             ac4 = qA.multiply(sinI2).add(
  706.                   qB.multiply(4.375).multiply(sinI2)).add(
  707.                   qC.multiply(0.75).multiply(sinI4.multiply(1.1).subtract(sinI2)));

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

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

  713.             ex1 = qq.multiply(sinI2.multiply(-1.25).add(1.0));
  714.             exx2 = qq.multiply(0.5).multiply(sinI2.multiply(-5.0).add(3.0));
  715.             exy2 = qq.multiply(sinI2.multiply(-1.5).add(2.0));
  716.             ex3 = qq.multiply(7.0 / 12.0).multiply(sinI2);
  717.             ex4 = qq.multiply(17.0 / 8.0).multiply(sinI2);

  719.             ey1 = qq.multiply(sinI2.multiply(-1.75).add(1.0));
  720.             eyx2 = qq.multiply(sinI2.multiply(-3.0).add(1.0));
  721.             eyy2 = qq.multiply(sinI2.multiply(2.0).subtract(1.5));
  722.             ey3 = qq.multiply(7.0 / 12.0).multiply(sinI2);
  723.             ey4 = qq.multiply(17.0 / 8.0).multiply(sinI2);

  725.             q  = cosI1.multiply(qq).negate();
  726.             rx1 = q.multiply(3.5);
  727.             ry1 = q.multiply(-2.5);
  728.             r2 = q.multiply(-0.5);
  729.             r3 =  q.multiply(7.0 / 6.0);
  730.             rl = g3 .multiply( cosI1).multiply(sinI2.multiply(-15.0).add(4.0)).subtract(
  731.                  g5.multiply(2.5).multiply(cosI1).multiply(sinI2.multiply(-42.0).add(4.0).add(sinI4.multiply(52.5))));

  733.             q = qq.multiply(0.5).multiply(sinI1).multiply(cosI1);
  734.             iy1 =  q;
  735.             ix1 = q.negate();
  736.             i2 =  q;
  737.             i3 =  q.multiply(7.0 / 3.0);
  738.             ih = g3.negate().multiply(cosI1).multiply(sinI2.multiply(-5.0).add(4)).add(
  739.                  g5.multiply(2.5).multiply(cosI1).multiply(sinI2.multiply(-14.0).add(4.0).add(sinI4.multiply(10.5))));
  740.             lx1 = qq.multiply(sinI2.multiply(-77.0 / 8.0).add(7.0));
  741.             ly1 = qq.multiply(sinI2.multiply(55.0 / 8.0).subtract(7.50));
  742.             l2 = qq.multiply(sinI2.multiply(1.25).subtract(0.5));
  743.             l3 = qq.multiply(sinI2.multiply(77.0 / 24.0).subtract(7.0 / 6.0));
  744.             ll = g3.multiply(sinI2.multiply(53.0).subtract(4.0).add(sinI4.multiply(-57.5))).add(
  745.                  g5.multiply(2.5).multiply(sinI2.multiply(-96.0).add(4.0).add(sinI4.multiply(269.5).subtract(sinI6.multiply(183.75)))));

  747.         }

  749.         /** Extrapolate a FieldOrbit up to a specific target date.
  750.          * @param date target date for the FieldOrbit
  751.          * @return propagated parameters
  752.          */
  753.         public FieldUnivariateDerivative2<T>[] propagateParameters(final FieldAbsoluteDate<T> date) {
  754.             final Field<T> field = date.durationFrom(mean.getDate()).getField();
  755.             final T one = field.getOne();
  756.             final T zero = field.getZero();
  757.             // Keplerian evolution
  758.             final FieldUnivariateDerivative2<T> dt = new FieldUnivariateDerivative2<>(date.durationFrom(mean.getDate()), one, zero);
  759.             final FieldUnivariateDerivative2<T> xnot = dt.multiply(xnotDot);

  761.             // secular effects

  763.             // eccentricity
  764.             final FieldUnivariateDerivative2<T> x   = xnot.multiply(rdpom.add(rdpomp));
  765.             final FieldUnivariateDerivative2<T> cx  = x.cos();
  766.             final FieldUnivariateDerivative2<T> sx  = x.sin();
  767.             final FieldUnivariateDerivative2<T> exm = cx.multiply(mean.getCircularEx()).
  768.                                             add(sx.multiply(eps2.subtract(one.subtract(eps1).multiply(mean.getCircularEy()))));
  769.             final FieldUnivariateDerivative2<T> eym = sx.multiply(eps1.add(1.0).multiply(mean.getCircularEx())).
  770.                                             add(cx.multiply(mean.getCircularEy().subtract(eps2))).
  771.                                             add(eps2);
  772.             // no secular effect on inclination

  774.             // right ascension of ascending node
  775.             final FieldUnivariateDerivative2<T> omm =
  776.                            new FieldUnivariateDerivative2<>(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode().add(ommD.multiply(xnot.getValue())),
  777.                                                                                      zero.add(FastMath.PI)),
  778.                                                             ommD.multiply(xnotDot),
  779.                                                             zero);
  780.             // latitude argument
  781.             final FieldUnivariateDerivative2<T> xlm =
  782.                             new FieldUnivariateDerivative2<>(MathUtils.normalizeAngle(mean.getAlphaM().add(aMD.multiply(xnot.getValue())),
  783.                                                                                     zero.add(FastMath.PI)),
  784.                                                            aMD.multiply(xnotDot),
  785.                                                            zero);

  787.             // periodical terms
  788.             final FieldUnivariateDerivative2<T> cl1 = xlm.cos();
  789.             final FieldUnivariateDerivative2<T> sl1 = xlm.sin();
  790.             final FieldUnivariateDerivative2<T> cl2 = cl1.multiply(cl1).subtract(sl1.multiply(sl1));
  791.             final FieldUnivariateDerivative2<T> sl2 = cl1.multiply(sl1).add(sl1.multiply(cl1));
  792.             final FieldUnivariateDerivative2<T> cl3 = cl2.multiply(cl1).subtract(sl2.multiply(sl1));
  793.             final FieldUnivariateDerivative2<T> sl3 = cl2.multiply(sl1).add(sl2.multiply(cl1));
  794.             final FieldUnivariateDerivative2<T> cl4 = cl3.multiply(cl1).subtract(sl3.multiply(sl1));
  795.             final FieldUnivariateDerivative2<T> sl4 = cl3.multiply(sl1).add(sl3.multiply(cl1));
  796.             final FieldUnivariateDerivative2<T> cl5 = cl4.multiply(cl1).subtract(sl4.multiply(sl1));
  797.             final FieldUnivariateDerivative2<T> sl5 = cl4.multiply(sl1).add(sl4.multiply(cl1));
  798.             final FieldUnivariateDerivative2<T> cl6 = cl5.multiply(cl1).subtract(sl5.multiply(sl1));

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

  803.             final FieldUnivariateDerivative2<T> exmCl1 = exm.multiply(cl1);
  804.             final FieldUnivariateDerivative2<T> exmSl1 = exm.multiply(sl1);
  805.             final FieldUnivariateDerivative2<T> eymCl1 = eym.multiply(cl1);
  806.             final FieldUnivariateDerivative2<T> eymSl1 = eym.multiply(sl1);
  807.             final FieldUnivariateDerivative2<T> exmCl2 = exm.multiply(cl2);
  808.             final FieldUnivariateDerivative2<T> exmSl2 = exm.multiply(sl2);
  809.             final FieldUnivariateDerivative2<T> eymCl2 = eym.multiply(cl2);
  810.             final FieldUnivariateDerivative2<T> eymSl2 = eym.multiply(sl2);
  811.             final FieldUnivariateDerivative2<T> exmCl3 = exm.multiply(cl3);
  812.             final FieldUnivariateDerivative2<T> exmSl3 = exm.multiply(sl3);
  813.             final FieldUnivariateDerivative2<T> eymCl3 = eym.multiply(cl3);
  814.             final FieldUnivariateDerivative2<T> eymSl3 = eym.multiply(sl3);
  815.             final FieldUnivariateDerivative2<T> exmCl4 = exm.multiply(cl4);
  816.             final FieldUnivariateDerivative2<T> exmSl4 = exm.multiply(sl4);
  817.             final FieldUnivariateDerivative2<T> eymCl4 = eym.multiply(cl4);
  818.             final FieldUnivariateDerivative2<T> eymSl4 = eym.multiply(sl4);

  820.             // semi major axis
  821.             final FieldUnivariateDerivative2<T> rda = exmCl1.multiply(ax1).
  822.                                             add(eymSl1.multiply(ay1)).
  823.                                             add(sl1.multiply(as1)).
  824.                                             add(cl2.multiply(ac2)).
  825.                                             add(exmCl3.add(eymSl3).multiply(axy3)).
  826.                                             add(sl3.multiply(as3)).
  827.                                             add(cl4.multiply(ac4)).
  828.                                             add(sl5.multiply(as5)).
  829.                                             add(cl6.multiply(ac6));

  831.             // eccentricity
  832.             final FieldUnivariateDerivative2<T> rdex = cl1.multiply(ex1).
  833.                                              add(exmCl2.multiply(exx2)).
  834.                                              add(eymSl2.multiply(exy2)).
  835.                                              add(cl3.multiply(ex3)).
  836.                                              add(exmCl4.add(eymSl4).multiply(ex4));
  837.             final FieldUnivariateDerivative2<T> rdey = sl1.multiply(ey1).
  838.                                              add(exmSl2.multiply(eyx2)).
  839.                                              add(eymCl2.multiply(eyy2)).
  840.                                              add(sl3.multiply(ey3)).
  841.                                              add(exmSl4.subtract(eymCl4).multiply(ey4));

  843.             // ascending node
  844.             final FieldUnivariateDerivative2<T> rdom = exmSl1.multiply(rx1).
  845.                                              add(eymCl1.multiply(ry1)).
  846.                                              add(sl2.multiply(r2)).
  847.                                              add(eymCl3.subtract(exmSl3).multiply(r3)).
  848.                                              add(ql.multiply(rl));

  850.             // inclination
  851.             final FieldUnivariateDerivative2<T> rdxi = eymSl1.multiply(iy1).
  852.                                              add(exmCl1.multiply(ix1)).
  853.                                              add(cl2.multiply(i2)).
  854.                                              add(exmCl3.add(eymSl3).multiply(i3)).
  855.                                              add(qh.multiply(ih));

  857.             // latitude argument
  858.             final FieldUnivariateDerivative2<T> rdxl = exmSl1.multiply(lx1).
  859.                                              add(eymCl1.multiply(ly1)).
  860.                                              add(sl2.multiply(l2)).
  861.                                              add(exmSl3.subtract(eymCl3).multiply(l3)).
  862.                                              add(ql.multiply(ll));
  863.             // osculating parameters
  864.             final FieldUnivariateDerivative2<T>[] FTD = MathArrays.buildArray(rdxl.getField(), 6);

  866.             FTD[0] = rda.add(1.0).multiply(mean.getA());
  867.             FTD[1] = rdex.add(exm);
  868.             FTD[2] = rdey.add(eym);
  869.             FTD[3] = rdxi.add(xim);
  870.             FTD[4] = rdom.add(omm);
  871.             FTD[5] = rdxl.add(xlm);
  872.             return FTD;

  874.         }

  876.     }

  878.     /** Convert circular parameters <em>with derivatives</em> to Cartesian coordinates.
  879.      * @param date date of the parameters
  880.      * @param parameters circular parameters (a, ex, ey, i, raan, alphaM)
  881.      * @return Cartesian coordinates consistent with values and derivatives
  882.      */
  883.     private TimeStampedFieldPVCoordinates<T> toCartesian(final FieldAbsoluteDate<T> date, final FieldUnivariateDerivative2<T>[] parameters) {

  885.         // evaluate coordinates in the FieldOrbit canonical reference frame
  886.         final FieldUnivariateDerivative2<T> cosOmega = parameters[4].cos();
  887.         final FieldUnivariateDerivative2<T> sinOmega = parameters[4].sin();
  888.         final FieldUnivariateDerivative2<T> cosI     = parameters[3].cos();
  889.         final FieldUnivariateDerivative2<T> sinI     = parameters[3].sin();
  890.         final FieldUnivariateDerivative2<T> alphaE   = meanToEccentric(parameters[5], parameters[1], parameters[2]);
  891.         final FieldUnivariateDerivative2<T> cosAE    = alphaE.cos();
  892.         final FieldUnivariateDerivative2<T> sinAE    = alphaE.sin();
  893.         final FieldUnivariateDerivative2<T> ex2      = parameters[1].multiply(parameters[1]);
  894.         final FieldUnivariateDerivative2<T> ey2      = parameters[2].multiply(parameters[2]);
  895.         final FieldUnivariateDerivative2<T> exy      = parameters[1].multiply(parameters[2]);
  896.         final FieldUnivariateDerivative2<T> q        = ex2.add(ey2).subtract(1).negate().sqrt();
  897.         final FieldUnivariateDerivative2<T> beta     = q.add(1).reciprocal();
  898.         final FieldUnivariateDerivative2<T> bx2      = beta.multiply(ex2);
  899.         final FieldUnivariateDerivative2<T> by2      = beta.multiply(ey2);
  900.         final FieldUnivariateDerivative2<T> bxy      = beta.multiply(exy);
  901.         final FieldUnivariateDerivative2<T> u        = bxy.multiply(sinAE).subtract(parameters[1].add(by2.subtract(1).multiply(cosAE)));
  902.         final FieldUnivariateDerivative2<T> v        = bxy.multiply(cosAE).subtract(parameters[2].add(bx2.subtract(1).multiply(sinAE)));
  903.         final FieldUnivariateDerivative2<T> x        = parameters[0].multiply(u);
  904.         final FieldUnivariateDerivative2<T> y        = parameters[0].multiply(v);

  906.         // canonical FieldOrbit reference frame
  907.         final FieldVector3D<FieldUnivariateDerivative2<T>> p =
  908.                 new FieldVector3D<>(x.multiply(cosOmega).subtract(y.multiply(cosI.multiply(sinOmega))),
  909.                                     x.multiply(sinOmega).add(y.multiply(cosI.multiply(cosOmega))),
  910.                                     y.multiply(sinI));

  912.         // dispatch derivatives
  913.         final FieldVector3D<T> p0 = new FieldVector3D<>(p.getX().getValue(),
  914.                                                         p.getY().getValue(),
  915.                                                         p.getZ().getValue());
  916.         final FieldVector3D<T> p1 = new FieldVector3D<>(p.getX().getFirstDerivative(),
  917.                                                         p.getY().getFirstDerivative(),
  918.                                                         p.getZ().getFirstDerivative());
  919.         final FieldVector3D<T> p2 = new FieldVector3D<>(p.getX().getSecondDerivative(),
  920.                                                         p.getY().getSecondDerivative(),
  921.                                                         p.getZ().getSecondDerivative());
  922.         return new TimeStampedFieldPVCoordinates<>(date, p0, p1, p2);

  924.     }

  926.     /** Computes the eccentric latitude argument from the mean latitude argument.
  927.      * @param alphaM = M + Ω mean latitude argument (rad)
  928.      * @param ex e cos(Ω), first component of circular eccentricity vector
  929.      * @param ey e sin(Ω), second component of circular eccentricity vector
  930.      * @return the eccentric latitude argument.
  931.      */
  932.     private FieldUnivariateDerivative2<T> meanToEccentric(final FieldUnivariateDerivative2<T> alphaM,
  933.                                                 final FieldUnivariateDerivative2<T> ex,
  934.                                                 final FieldUnivariateDerivative2<T> ey) {
  935.         // Generalization of Kepler equation to circular parameters
  936.         // with alphaE = PA + E and
  937.         //      alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
  938.         FieldUnivariateDerivative2<T> alphaE        = alphaM;
  939.         FieldUnivariateDerivative2<T> shift         = alphaM.getField().getZero();
  940.         FieldUnivariateDerivative2<T> alphaEMalphaM = alphaM.getField().getZero();
  941.         FieldUnivariateDerivative2<T> cosAlphaE     = alphaE.cos();
  942.         FieldUnivariateDerivative2<T> sinAlphaE     = alphaE.sin();
  943.         int                 iter          = 0;
  944.         do {
  945.             final FieldUnivariateDerivative2<T> f2 = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
  946.             final FieldUnivariateDerivative2<T> f1 = alphaM.getField().getOne().subtract(ex.multiply(cosAlphaE)).subtract(ey.multiply(sinAlphaE));
  947.             final FieldUnivariateDerivative2<T> f0 = alphaEMalphaM.subtract(f2);

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

  952.             alphaEMalphaM  = alphaEMalphaM.subtract(shift);
  953.             alphaE         = alphaM.add(alphaEMalphaM);
  954.             cosAlphaE      = alphaE.cos();
  955.             sinAlphaE      = alphaE.sin();

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

  959.         return alphaE;

  961.     }

  963.     /** {@inheritDoc} */
  964.     protected T getMass(final FieldAbsoluteDate<T> date) {
  965.         return models.get(date).mass;
  966.     }
  967.     /**
  968.      * Normalize an angle in a 2&pi; wide interval around a center value.
  969.      * <p>This method has three main uses:</p>
  970.      * <ul>
  971.      *   <li>normalize an angle between 0 and 2&pi;:<br>
  972.      *       {@code a = MathUtils.normalizeAngle(a, FastMath.PI);}</li>
  973.      *   <li>normalize an angle between -&pi; and +&pi;<br>
  974.      *       {@code a = MathUtils.normalizeAngle(a, 0.0);}</li>
  975.      *   <li>compute the angle between two defining angular positions:<br>
  976.      *       {@code angle = MathUtils.normalizeAngle(end, start) - start;}</li>
  977.      * </ul>
  978.      * <p>Note that due to numerical accuracy and since &pi; cannot be represented
  979.      * exactly, the result interval is <em>closed</em>, it cannot be half-closed
  980.      * as would be more satisfactory in a purely mathematical view.</p>
  981.      * @param a angle to normalize
  982.      * @param center center of the desired 2&pi; interval for the result
  983.      * @param <T> the type of the field elements
  984.      * @return a-2k&pi; with integer k and center-&pi; &lt;= a-2k&pi; &lt;= center+&pi;
  985.      * @since 1.2
  986.      * @deprecated replaced by {@link MathUtils#normalizeAngle(RealFieldElement, RealFieldElement)}
  987.      */
  988.     @Deprecated
  989.     public static <T extends RealFieldElement<T>> T normalizeAngle(final T a, final T center) {
  990.         return a.subtract(2 * FastMath.PI * FastMath.floor((a.getReal() + FastMath.PI - center.getReal()) / (2 * FastMath.PI)));
  991.     }

  993. }
Created with JaCoCo 0.8.6.202009150832