FieldDSSTTesseralContext.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.semianalytical.dsst.forces;

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

  22. import org.hipparchus.Field;
  23. import org.hipparchus.RealFieldElement;
  24. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  25. import org.hipparchus.util.FastMath;
  26. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  27. import org.orekit.frames.FieldTransform;
  28. import org.orekit.frames.Frame;
  29. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;

  31. /**
  32.  * This class is a container for the common "field" parameters used in {@link DSSTTesseral}.
  33.  * <p>
  34.  * It performs parameters initialization at each integration step for the Tesseral contribution
  35.  * to the central body gravitational perturbation.
  36.  * <p>
  37.  * @author Bryan Cazabonne
  38.  * @since 10.0
  39.  */
  40. class FieldDSSTTesseralContext<T extends RealFieldElement<T>> extends FieldForceModelContext<T> {

  42.     /** Retrograde factor I.
  43.      *  <p>
  44.      *  DSST model needs equinoctial orbit as internal representation.
  45.      *  Classical equinoctial elements have discontinuities when inclination
  46.      *  is close to zero. In this representation, I = +1. <br>
  47.      *  To avoid this discontinuity, another representation exists and equinoctial
  48.      *  elements can be expressed in a different way, called "retrograde" orbit.
  49.      *  This implies I = -1. <br>
  50.      *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  51.      *  But for the sake of consistency with the theory, the retrograde factor
  52.      *  has been kept in the formulas.
  53.      *  </p>
  54.      */
  55.     private static final int I = 1;

  57.     /** Minimum period for analytically averaged high-order resonant
  58.      *  central body spherical harmonics in seconds.
  59.      */
  60.     private static final double MIN_PERIOD_IN_SECONDS = 864000.;

  62.     /** Minimum period for analytically averaged high-order resonant
  63.      *  central body spherical harmonics in satellite revolutions.
  64.      */
  65.     private static final double MIN_PERIOD_IN_SAT_REV = 10.;

  67.     /** A = sqrt(μ * a). */
  68.     private T A;

  70.     // Common factors for potential computation
  71.     /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
  72.     private T chi;

  74.     /** &Chi;². */
  75.     private T chi2;

  77.     /** Central body rotation angle θ. */
  78.     private T theta;

  80.     // Common factors from equinoctial coefficients
  81.     /** 2 * a / A .*/
  82.     private T ax2oA;

  84.     /** 1 / (A * B) .*/
  85.     private T ooAB;

  87.     /** B / A .*/
  88.     private T BoA;

  90.     /** B / (A * (1 + B)) .*/
  91.     private T BoABpo;

  93.     /** C / (2 * A * B) .*/
  94.     private T Co2AB;

  96.     /** μ / a .*/
  97.     private T moa;

  99.     /** R / a .*/
  100.     private T roa;

  102.     /** ecc². */
  103.     private T e2;

  105.     /** Keplerian mean motion. */
  106.     private T n;

  108.     /** Keplerian period. */
  109.     private T period;

  111.     /** Maximum power of the eccentricity to use in summation over s. */
  112.     private int maxEccPow;

  114.     /** Ratio of satellite period to central body rotation period. */
  115.     private T ratio;

  117.     /** List of resonant orders. */
  118.     private final List<Integer> resOrders;

  120.     /**
  121.      * Simple constructor.
  122.      *
  123.      * @param auxiliaryElements auxiliary elements related to the current orbit
  124.      * @param centralBodyFrame rotating body frame
  125.      * @param provider provider for spherical harmonics
  126.      * @param maxFrequencyShortPeriodics maximum value for j
  127.      * @param bodyPeriod central body rotation period (seconds)
  128.      * @param parameters values of the force model parameters
  129.      */
  130.     FieldDSSTTesseralContext(final FieldAuxiliaryElements<T> auxiliaryElements,
  131.                                     final Frame centralBodyFrame,
  132.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  133.                                     final int maxFrequencyShortPeriodics,
  134.                                     final double bodyPeriod,
  135.                                     final T[] parameters) {

  137.         super(auxiliaryElements);

  139.         final Field<T> field = auxiliaryElements.getDate().getField();
  140.         final T zero = field.getZero();

  142.         this.maxEccPow = 0;
  143.         this.resOrders = new ArrayList<Integer>();

  145.         final T mu = parameters[0];

  147.         // Keplerian mean motion
  148.         final T absA = FastMath.abs(auxiliaryElements.getSma());
  149.         n = FastMath.sqrt(mu.divide(absA)).divide(absA);

  151.         // Keplerian period
  152.         final T a = auxiliaryElements.getSma();
  153.         period = (a.getReal() < 0) ? zero.add(Double.POSITIVE_INFINITY) : a.multiply(2.0 * FastMath.PI).multiply(a.divide(mu).sqrt());

  155.         A = FastMath.sqrt(mu.multiply(auxiliaryElements.getSma()));

  157.         // Eccentricity square
  158.         e2 = auxiliaryElements.getEcc().multiply(auxiliaryElements.getEcc());

  160.         // Central body rotation angle from equation 2.7.1-(3)(4).
  161.         final FieldTransform<T> t = centralBodyFrame.getTransformTo(auxiliaryElements.getFrame(), auxiliaryElements.getDate());
  162.         final FieldVector3D<T> xB = t.transformVector(FieldVector3D.getPlusI(field));
  163.         final FieldVector3D<T> yB = t.transformVector(FieldVector3D.getPlusJ(field));
  164.         theta = FastMath.atan2(auxiliaryElements.getVectorF().dotProduct(yB).negate().add((auxiliaryElements.getVectorG().dotProduct(xB)).multiply(I)),
  165.                                auxiliaryElements.getVectorF().dotProduct(xB).add(auxiliaryElements.getVectorG().dotProduct(yB).multiply(I)));

  167.         // Common factors from equinoctial coefficients
  168.         // 2 * a / A
  169.         ax2oA  = auxiliaryElements.getSma().divide(A).multiply(2.);
  170.         // B / A
  171.         BoA    = auxiliaryElements.getB().divide(A);
  172.         // 1 / AB
  173.         ooAB   = A.multiply(auxiliaryElements.getB()).reciprocal();
  174.         // C / 2AB
  175.         Co2AB  = auxiliaryElements.getC().multiply(ooAB).divide(2.);
  176.         // B / (A * (1 + B))
  177.         BoABpo = BoA.divide(auxiliaryElements.getB().add(1.));
  178.         // &mu / a
  179.         moa    = mu.divide(auxiliaryElements.getSma());
  180.         // R / a
  181.         roa    = auxiliaryElements.getSma().divide(provider.getAe()).reciprocal();

  183.         // &Chi; = 1 / B
  184.         chi  = auxiliaryElements.getB().reciprocal();
  185.         chi2 = chi.multiply(chi);

  187.         // Set the highest power of the eccentricity in the analytical power
  188.         // series expansion for the averaged high order resonant central body
  189.         // spherical harmonic perturbation
  190.         final T e = auxiliaryElements.getEcc();
  191.         if (e.getReal() <= 0.005) {
  192.             maxEccPow = 3;
  193.         } else if (e.getReal() <= 0.02) {
  194.             maxEccPow = 4;
  195.         } else if (e.getReal() <= 0.1) {
  196.             maxEccPow = 7;
  197.         } else if (e.getReal() <= 0.2) {
  198.             maxEccPow = 10;
  199.         } else if (e.getReal() <= 0.3) {
  200.             maxEccPow = 12;
  201.         } else if (e.getReal() <= 0.4) {
  202.             maxEccPow = 15;
  203.         } else {
  204.             maxEccPow = 20;
  205.         }

  207.         // Ratio of satellite to central body periods to define resonant terms
  208.         ratio = period.divide(bodyPeriod);

  210.         // Compute natural resonant terms
  211.         final T tolerance = FastMath.max(zero.add(MIN_PERIOD_IN_SAT_REV),
  212.                                                    period.divide(MIN_PERIOD_IN_SECONDS).reciprocal()).reciprocal();

  214.         // Search the resonant orders in the tesseral harmonic field
  215.         resOrders.clear();
  216.         for (int m = 1; m <= provider.getMaxOrder(); m++) {
  217.             final T resonance = ratio.multiply(m);
  218.             final int jComputedRes = (int) FastMath.round(resonance);
  219.             if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance.subtract(jComputedRes)).getReal() <= tolerance.getReal()) {
  220.                 // Store each resonant index and order
  221.                 this.resOrders.add(m);
  222.             }
  223.         }

  225.     }

  227.     /** Get the list of resonant orders.
  228.      * @return resOrders
  229.      */
  230.     public List<Integer> getResOrders() {
  231.         return resOrders;
  232.     }

  234.     /** Get ecc².
  235.      * @return e2
  236.      */
  237.     public T getE2() {
  238.         return e2;
  239.     }

  241.     /** Get Central body rotation angle θ.
  242.      * @return theta
  243.      */
  244.     public T getTheta() {
  245.         return theta;
  246.     }

  248.     /** Get ax2oA = 2 * a / A .
  249.      * @return ax2oA
  250.      */
  251.     public T getAx2oA() {
  252.         return ax2oA;
  253.     }

  255.     /** Get &Chi; = 1 / sqrt(1 - e²) = 1 / B.
  256.      * @return chi
  257.      */
  258.     public T getChi() {
  259.         return chi;
  260.     }

  262.     /** Get &Chi;².
  263.      * @return chi2
  264.      */
  265.     public T getChi2() {
  266.         return chi2;
  267.     }

  269.     /** Get B / A.
  270.      * @return BoA
  271.      */
  272.     public T getBoA() {
  273.         return BoA;
  274.     }

  276.     /** Get ooAB = 1 / (A * B).
  277.      * @return ooAB
  278.      */
  279.     public T getOoAB() {
  280.         return ooAB;
  281.     }

  283.     /** Get Co2AB = C / 2AB.
  284.      * @return Co2AB
  285.      */
  286.     public T getCo2AB() {
  287.         return Co2AB;
  288.     }

  290.     /** Get BoABpo = B / A(1 + B).
  291.      * @return BoABpo
  292.      */
  293.     public T getBoABpo() {
  294.         return BoABpo;
  295.     }

  297.     /** Get μ / a .
  298.      * @return moa
  299.      */
  300.     public T getMoa() {
  301.         return moa;
  302.     }

  304.     /** Get roa = R / a.
  305.      * @return roa
  306.      */
  307.     public T getRoa() {
  308.         return roa;
  309.     }

  311.     /** Get the maximum power of the eccentricity to use in summation over s.
  312.      * @return roa
  313.      */
  314.     public int getMaxEccPow() {
  315.         return maxEccPow;
  316.     }

  318.     /** Get the Keplerian period.
  319.      * <p>The Keplerian period is computed directly from semi major axis
  320.      * and central acceleration constant.</p>
  321.      * @return Keplerian period in seconds, or positive infinity for hyperbolic orbits
  322.      */
  323.     public T getOrbitPeriod() {
  324.         return period;
  325.     }

  327.     /** Get the Keplerian mean motion.
  328.      * <p>The Keplerian mean motion is computed directly from semi major axis
  329.      * and central acceleration constant.</p>
  330.      * @return Keplerian mean motion in radians per second
  331.      */
  332.     public T getMeanMotion() {
  333.         return n;
  334.     }

  336.     /** Get the ratio of satellite period to central body rotation period.
  337.      * @return ratio
  338.      */
  339.     public T getRatio() {
  340.         return ratio;
  341.     }

  343. }
Created with JaCoCo 0.8.6.202009150832