DSSTTesseralContext.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.geometry.euclidean.threed.Vector3D;
  23. import org.hipparchus.util.FastMath;
  24. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  25. import org.orekit.frames.Frame;
  26. import org.orekit.frames.Transform;
  27. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;

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

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

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

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

  65.     /** A = sqrt(μ * a). */
  66.     private double A;

  68.     // Common factors for potential computation
  69.     /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
  70.     private double chi;

  72.     /** &Chi;². */
  73.     private double chi2;

  75.     /** Central body rotation angle θ. */
  76.     private double theta;

  78.     // Common factors from equinoctial coefficients
  79.     /** 2 * a / A .*/
  80.     private double ax2oA;

  82.     /** 1 / (A * B) .*/
  83.     private double ooAB;

  85.     /** B / A .*/
  86.     private double BoA;

  88.     /** B / (A * (1 + B)) .*/
  89.     private double BoABpo;

  91.     /** C / (2 * A * B) .*/
  92.     private double Co2AB;

  94.     /** μ / a .*/
  95.     private double moa;

  97.     /** R / a .*/
  98.     private double roa;

  100.     /** ecc². */
  101.     private double e2;

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

  106.     /** Keplerian mean motion. */
  107.     private double n;

  109.     /** Keplerian period. */
  110.     private double period;

  112.     /** Ratio of satellite period to central body rotation period. */
  113.     private double ratio;

  115.     /** List of resonant orders. */
  116.     private final List<Integer> resOrders;

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

  135.         super(auxiliaryElements);

  137.         this.maxEccPow = 0;
  138.         this.resOrders = new ArrayList<Integer>();

  140.         final double mu = parameters[0];

  142.         // Keplerian Mean Motion
  143.         final double absA = FastMath.abs(auxiliaryElements.getSma());
  144.         n = FastMath.sqrt(mu / absA) / absA;

  146.         // Keplerian period
  147.         final double a = auxiliaryElements.getSma();
  148.         period = (a < 0) ? Double.POSITIVE_INFINITY : 2.0 * FastMath.PI * a * FastMath.sqrt(a / mu);

  150.         A = FastMath.sqrt(mu * auxiliaryElements.getSma());

  152.         // Eccentricity square
  153.         e2 = auxiliaryElements.getEcc() * auxiliaryElements.getEcc();

  155.         // Central body rotation angle from equation 2.7.1-(3)(4).
  156.         final Transform t = centralBodyFrame.getTransformTo(auxiliaryElements.getFrame(), auxiliaryElements.getDate());
  157.         final Vector3D xB = t.transformVector(Vector3D.PLUS_I);
  158.         final Vector3D yB = t.transformVector(Vector3D.PLUS_J);
  159.         theta = FastMath.atan2(-auxiliaryElements.getVectorF().dotProduct(yB) + I * auxiliaryElements.getVectorG().dotProduct(xB),
  160.                                auxiliaryElements.getVectorF().dotProduct(xB) + I * auxiliaryElements.getVectorG().dotProduct(yB));

  162.         // Common factors from equinoctial coefficients
  163.         // 2 * a / A
  164.         ax2oA  = 2. * auxiliaryElements.getSma() / A;
  165.         // B / A
  166.         BoA    = auxiliaryElements.getB() / A;
  167.         // 1 / AB
  168.         ooAB   = 1. / (A * auxiliaryElements.getB());
  169.         // C / 2AB
  170.         Co2AB  = auxiliaryElements.getC() * ooAB / 2.;
  171.         // B / (A * (1 + B))
  172.         BoABpo = BoA / (1. + auxiliaryElements.getB());
  173.         // &mu / a
  174.         moa    = mu / auxiliaryElements.getSma();
  175.         // R / a
  176.         roa    = provider.getAe() / auxiliaryElements.getSma();

  178.         // &Chi; = 1 / B
  179.         chi  = 1. / auxiliaryElements.getB();
  180.         chi2 = chi * chi;

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

  202.         // Ratio of satellite to central body periods to define resonant terms
  203.         ratio = period / bodyPeriod;

  205.         // Compute natural resonant terms
  206.         final double tolerance = 1. / FastMath.max(MIN_PERIOD_IN_SAT_REV,
  207.                                                    MIN_PERIOD_IN_SECONDS / period);

  209.         // Search the resonant orders in the tesseral harmonic field
  210.         resOrders.clear();
  211.         for (int m = 1; m <= provider.getMaxOrder(); m++) {
  212.             final double resonance = ratio * m;
  213.             final int jComputedRes = (int) FastMath.round(resonance);
  214.             if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance - jComputedRes) <= tolerance) {
  215.                 // Store each resonant index and order
  216.                 this.resOrders.add(m);
  217.             }
  218.         }

  220.     }

  222.     /** Get the list of resonant orders.
  223.      * @return resOrders
  224.      */
  225.     public List<Integer> getResOrders() {
  226.         return resOrders;
  227.     }

  229.     /** Get ecc².
  230.      * @return e2
  231.      */
  232.     public double getE2() {
  233.         return e2;
  234.     }

  236.     /** Get Central body rotation angle θ.
  237.      * @return theta
  238.      */
  239.     public double getTheta() {
  240.         return theta;
  241.     }

  243.     /** Get ax2oA = 2 * a / A .
  244.      * @return ax2oA
  245.      */
  246.     public double getAx2oA() {
  247.         return ax2oA;
  248.     }

  250.     /** Get &Chi; = 1 / sqrt(1 - e²) = 1 / B.
  251.      * @return chi
  252.      */
  253.     public double getChi() {
  254.         return chi;
  255.     }

  257.     /** Get &Chi;².
  258.      * @return chi2
  259.      */
  260.     public double getChi2() {
  261.         return chi2;
  262.     }

  264.     /** Get B / A.
  265.      * @return BoA
  266.      */
  267.     public double getBoA() {
  268.         return BoA;
  269.     }

  271.     /** Get ooAB = 1 / (A * B).
  272.      * @return ooAB
  273.      */
  274.     public double getOoAB() {
  275.         return ooAB;
  276.     }

  278.     /** Get Co2AB = C / 2AB.
  279.      * @return Co2AB
  280.      */
  281.     public double getCo2AB() {
  282.         return Co2AB;
  283.     }

  285.     /** Get BoABpo = B / A(1 + B).
  286.      * @return BoABpo
  287.      */
  288.     public double getBoABpo() {
  289.         return BoABpo;
  290.     }

  292.     /** Get μ / a .
  293.      * @return moa
  294.      */
  295.     public double getMoa() {
  296.         return moa;
  297.     }

  299.     /** Get roa = R / a.
  300.      * @return roa
  301.      */
  302.     public double getRoa() {
  303.         return roa;
  304.     }

  306.     /** Get the maximum power of the eccentricity to use in summation over s.
  307.      * @return maxEccPow
  308.      */
  309.     public int getMaxEccPow() {
  310.         return maxEccPow;
  311.     }

  313.     /** Get the Keplerian period.
  314.      * <p>The Keplerian period is computed directly from semi major axis
  315.      * and central acceleration constant.</p>
  316.      * @return Keplerian period in seconds, or positive infinity for hyperbolic orbits
  317.      */
  318.     public double getOrbitPeriod() {
  319.         return period;
  320.     }

  322.     /** Get the Keplerian mean motion.
  323.      * <p>The Keplerian mean motion is computed directly from semi major axis
  324.      * and central acceleration constant.</p>
  325.      * @return Keplerian mean motion in radians per second
  326.      */
  327.     public double getMeanMotion() {
  328.         return n;
  329.     }

  331.     /** Get the ratio of satellite period to central body rotation period.
  332.      * @return ratio
  333.      */
  334.     public double getRatio() {
  335.         return ratio;
  336.     }

  338. }
Created with JaCoCo 0.8.6.202009150832