DSSTTesseralContext.java

  1. /* Copyright 2002-2022 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;

  18. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  19. import org.hipparchus.util.FastMath;
  20. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  21. import org.orekit.frames.Frame;
  22. import org.orekit.frames.StaticTransform;
  23. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;

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

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

  48.     /** A = sqrt(μ * a). */
  49.     private double A;

  50.     // Common factors for potential computation
  51.     /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
  52.     private double chi;

  53.     /** &Chi;². */
  54.     private double chi2;

  55.     /** Central body rotation angle θ. */
  56.     private double theta;

  57.     // Common factors from equinoctial coefficients
  58.     /** 2 * a / A . */
  59.     private double ax2oA;

  60.     /** 1 / (A * B) . */
  61.     private double ooAB;

  62.     /** B / A . */
  63.     private double BoA;

  64.     /** B / (A * (1 + B)) . */
  65.     private double BoABpo;

  66.     /** C / (2 * A * B) . */
  67.     private double Co2AB;

  68.     /** μ / a . */
  69.     private double moa;

  70.     /** R / a . */
  71.     private double roa;

  72.     /** ecc². */
  73.     private double e2;

  74.     /** Keplerian mean motion. */
  75.     private double n;

  76.     /** Keplerian period. */
  77.     private double period;

  78.     /** Ratio of satellite period to central body rotation period. */
  79.     private double ratio;

  80.     /**
  81.      * Simple constructor.
  82.      *
  83.      * @param auxiliaryElements auxiliary elements related to the current orbit
  84.      * @param centralBodyFrame           rotating body frame
  85.      * @param provider                   provider for spherical harmonics
  86.      * @param maxFrequencyShortPeriodics maximum value for j
  87.      * @param bodyPeriod                 central body rotation period (seconds)
  88.      * @param parameters                 values of the force model parameters
  89.      */
  90.     DSSTTesseralContext(final AuxiliaryElements auxiliaryElements,
  91.                         final Frame centralBodyFrame,
  92.                         final UnnormalizedSphericalHarmonicsProvider provider,
  93.                         final int maxFrequencyShortPeriodics,
  94.                         final double bodyPeriod,
  95.                         final double[] parameters) {

  96.         super(auxiliaryElements);

  97.         final double mu = parameters[0];

  98.         // Keplerian Mean Motion
  99.         final double absA = FastMath.abs(auxiliaryElements.getSma());
  100.         n = FastMath.sqrt(mu / absA) / absA;

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

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

  105.         // Eccentricity square
  106.         e2 = auxiliaryElements.getEcc() * auxiliaryElements.getEcc();

  107.         // Central body rotation angle from equation 2.7.1-(3)(4).
  108.         final StaticTransform t = centralBodyFrame.getStaticTransformTo(
  109.                 auxiliaryElements.getFrame(),
  110.                 auxiliaryElements.getDate());
  111.         final Vector3D xB = t.transformVector(Vector3D.PLUS_I);
  112.         final Vector3D yB = t.transformVector(Vector3D.PLUS_J);
  113.         theta = FastMath.atan2(-auxiliaryElements.getVectorF().dotProduct(yB) + I * auxiliaryElements.getVectorG().dotProduct(xB),
  114.                 auxiliaryElements.getVectorF().dotProduct(xB) + I * auxiliaryElements.getVectorG().dotProduct(yB));

  115.         // Common factors from equinoctial coefficients
  116.         // 2 * a / A
  117.         ax2oA = 2. * auxiliaryElements.getSma() / A;
  118.         // B / A
  119.         BoA = auxiliaryElements.getB() / A;
  120.         // 1 / AB
  121.         ooAB = 1. / (A * auxiliaryElements.getB());
  122.         // C / 2AB
  123.         Co2AB = auxiliaryElements.getC() * ooAB / 2.;
  124.         // B / (A * (1 + B))
  125.         BoABpo = BoA / (1. + auxiliaryElements.getB());
  126.         // &mu / a
  127.         moa = mu / auxiliaryElements.getSma();
  128.         // R / a
  129.         roa = provider.getAe() / auxiliaryElements.getSma();

  130.         // &Chi; = 1 / B
  131.         chi = 1. / auxiliaryElements.getB();
  132.         chi2 = chi * chi;

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

  135.     }

  136.     /** Get ecc².
  137.      * @return e2
  138.      */
  139.     public double getE2() {
  140.         return e2;
  141.     }

  142.     /**
  143.      * Get Central body rotation angle θ.
  144.      * @return theta
  145.      */
  146.     public double getTheta() {
  147.         return theta;
  148.     }

  149.     /**
  150.      * Get ax2oA = 2 * a / A .
  151.      * @return ax2oA
  152.      */
  153.     public double getAx2oA() {
  154.         return ax2oA;
  155.     }

  156.     /**
  157.      * Get &Chi; = 1 / sqrt(1 - e²) = 1 / B.
  158.      * @return chi
  159.      */
  160.     public double getChi() {
  161.         return chi;
  162.     }

  163.     /**
  164.      * Get &Chi;².
  165.      * @return chi2
  166.      */
  167.     public double getChi2() {
  168.         return chi2;
  169.     }

  170.     /**
  171.      * Get B / A.
  172.      * @return BoA
  173.      */
  174.     public double getBoA() {
  175.         return BoA;
  176.     }

  177.     /**
  178.      * Get ooAB = 1 / (A * B).
  179.      * @return ooAB
  180.      */
  181.     public double getOoAB() {
  182.         return ooAB;
  183.     }

  184.     /**
  185.      * Get Co2AB = C / 2AB.
  186.      * @return Co2AB
  187.      */
  188.     public double getCo2AB() {
  189.         return Co2AB;
  190.     }

  191.     /**
  192.      * Get BoABpo = B / A(1 + B).
  193.      * @return BoABpo
  194.      */
  195.     public double getBoABpo() {
  196.         return BoABpo;
  197.     }

  198.     /**
  199.      * Get μ / a .
  200.      * @return moa
  201.      */
  202.     public double getMoa() {
  203.         return moa;
  204.     }

  205.     /**
  206.      * Get roa = R / a.
  207.      * @return roa
  208.      */
  209.     public double getRoa() {
  210.         return roa;
  211.     }

  212.     /**
  213.      * Get the Keplerian period.
  214.      * <p>
  215.      * The Keplerian period is computed directly from semi major axis and central
  216.      * acceleration constant.
  217.      * </p>
  218.      * @return Keplerian period in seconds, or positive infinity for hyperbolic
  219.      *         orbits
  220.      */
  221.     public double getOrbitPeriod() {
  222.         return period;
  223.     }

  224.     /**
  225.      * Get the Keplerian mean motion.
  226.      * <p>
  227.      * The Keplerian mean motion is computed directly from semi major axis and
  228.      * central acceleration constant.
  229.      * </p>
  230.      * @return Keplerian mean motion in radians per second
  231.      */
  232.     public double getMeanMotion() {
  233.         return n;
  234.     }

  235.     /**
  236.      * Get the ratio of satellite period to central body rotation period.
  237.      * @return ratio
  238.      */
  239.     public double getRatio() {
  240.         return ratio;
  241.     }

  242. }