DSSTThirdBodyContext.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 org.hipparchus.geometry.euclidean.threed.Vector3D;
  20. import org.hipparchus.util.CombinatoricsUtils;
  21. import org.hipparchus.util.FastMath;
  22. import org.orekit.bodies.CelestialBody;
  23. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  24. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  25. import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;

  27. /**
  28.  * This class is a container for the common parameters used in {@link DSSTThirdBody}.
  29.  * <p>
  30.  * It performs parameters initialization at each integration step for the third
  31.  * body attraction perturbation.
  32.  * <p>
  33.  * @author Bryan Cazabonne
  34.  * @since 10.0
  35.  */
  36. class DSSTThirdBodyContext extends ForceModelContext {

  38.     /** Max power for summation. */
  39.     private static final int    MAX_POWER = 22;

  41.     /** Truncation tolerance for big, eccentric  orbits. */
  42.     private static final double BIG_TRUNCATION_TOLERANCE = 1.e-1;

  44.     /** Truncation tolerance for small orbits. */
  45.     private static final double SMALL_TRUNCATION_TOLERANCE = 1.9e-6;

  47.     /** Maximum power for eccentricity used in short periodic computation. */
  48.     private static final int    MAX_ECCPOWER_SP = 4;

  50.     /** Max power for a/R3 in the serie expansion. */
  51.     private int    maxAR3Pow;

  53.     /** Max power for e in the serie expansion. */
  54.     private int    maxEccPow;

  56.     /** a / R3 up to power maxAR3Pow. */
  57.     private double[] aoR3Pow;

  59.     /** Max power for e in the serie expansion (for short periodics). */
  60.     private int    maxEccPowShort;

  62.     /** Max frequency of F. */
  63.     private int    maxFreqF;

  65.     /** Qns coefficients. */
  66.     private double[][] Qns;

  68.     /** Standard gravitational parameter μ for the body in m³/s². */
  69.     private final double gm;

  71.     /** Distance from center of mass of the central body to the 3rd body. */
  72.     private double R3;

  74.     /** A = sqrt(μ * a). */
  75.     private final double A;

  77.     // Direction cosines of the symmetry axis
  78.     /** α. */
  79.     private double alpha;
  80.     /** β. */
  81.     private double beta;
  82.     /** γ. */
  83.     private double gamma;

  85.     /** B². */
  86.     private double BB;
  87.     /** B³. */
  88.     private double BBB;

  90.     /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
  91.     private double X;
  92.     /** &Chi;². */
  93.     private double XX;
  94.     /** &Chi;³. */
  95.     private double XXX;
  96.     /** -2 * a / A. */
  97.     private double m2aoA;
  98.     /** B / A. */
  99.     private double BoA;
  100.     /** 1 / (A * B). */
  101.     private double ooAB;
  102.     /** -C / (2 * A * B). */
  103.     private double mCo2AB;
  104.     /** B / A(1 + B). */
  105.     private double BoABpo;

  107.     /** mu3 / R3. */
  108.     private double muoR3;

  110.     /** b = 1 / (1 + sqrt(1 - e²)) = 1 / (1 + B).*/
  111.     private double b;

  113.     /** h * &Chi;³. */
  114.     private double hXXX;
  115.     /** k * &Chi;³. */
  116.     private double kXXX;

  118.     /** Keplerian mean motion. */
  119.     private final double motion;

  121.     /**
  122.      * Simple constructor.
  123.      *
  124.      * @param auxiliaryElements auxiliary elements related to the current orbit
  125.      * @param thirdBody body the 3rd body to consider
  126.      * @param parameters values of the force model parameters
  127.      */
  128.     DSSTThirdBodyContext(final AuxiliaryElements auxiliaryElements, final CelestialBody thirdBody, final double[] parameters) {

  130.         super(auxiliaryElements);

  132.         final double mu = parameters[1];
  133.         A = FastMath.sqrt(mu * auxiliaryElements.getSma());
  134.         this.gm = parameters[0];

  136.         // Keplerian Mean Motion
  137.         final double absA = FastMath.abs(auxiliaryElements.getSma());
  138.         motion = FastMath.sqrt(mu / absA) / absA;

  140.         // Distance from center of mass of the central body to the 3rd body
  141.         final Vector3D bodyPos = thirdBody.getPVCoordinates(auxiliaryElements.getDate(), auxiliaryElements.getFrame()).getPosition();
  142.         R3 = bodyPos.getNorm();

  144.         // Direction cosines
  145.         final Vector3D bodyDir = bodyPos.normalize();
  146.         alpha = bodyDir.dotProduct(auxiliaryElements.getVectorF());
  147.         beta  = bodyDir.dotProduct(auxiliaryElements.getVectorG());
  148.         gamma = bodyDir.dotProduct(auxiliaryElements.getVectorW());

  150.         //&Chi;<sup>-2</sup>.
  151.         BB = auxiliaryElements.getB() * auxiliaryElements.getB();
  152.         //&Chi;<sup>-3</sup>.
  153.         BBB = BB * auxiliaryElements.getB();

  155.         //b = 1 / (1 + B)
  156.         b = 1. / (1. + auxiliaryElements.getB());

  158.         // &Chi;
  159.         X = 1. / auxiliaryElements.getB();
  160.         XX = X * X;
  161.         XXX = X * XX;
  162.         // -2 * a / A
  163.         m2aoA = -2. * auxiliaryElements.getSma() / A;
  164.         // B / A
  165.         BoA = auxiliaryElements.getB() / A;
  166.         // 1 / AB
  167.         ooAB = 1. / (A * auxiliaryElements.getB());
  168.         // -C / 2AB
  169.         mCo2AB = -auxiliaryElements.getC() * ooAB / 2.;
  170.         // B / A(1 + B)
  171.         BoABpo = BoA / (1. + auxiliaryElements.getB());

  173.         // mu3 / R3
  174.         muoR3 = gm / R3;

  176.         //h * &Chi;³
  177.         hXXX = auxiliaryElements.getH() * XXX;
  178.         //k * &Chi;³
  179.         kXXX = auxiliaryElements.getK() * XXX;

  181.         // Truncation tolerance.
  182.         final double aoR3 = auxiliaryElements.getSma() / R3;
  183.         final double tol = ( aoR3 > .3 || (aoR3 > .15  && auxiliaryElements.getEcc() > .25) ) ? BIG_TRUNCATION_TOLERANCE : SMALL_TRUNCATION_TOLERANCE;

  185.         // Utilities for truncation
  186.         // Set a lower bound for eccentricity
  187.         final double eo2  = FastMath.max(0.0025, auxiliaryElements.getEcc() / 2.);
  188.         final double x2o2 = XX / 2.;
  189.         final double[] eccPwr = new double[MAX_POWER];
  190.         final double[] chiPwr = new double[MAX_POWER];
  191.         eccPwr[0] = 1.;
  192.         chiPwr[0] = X;
  193.         for (int i = 1; i < MAX_POWER; i++) {
  194.             eccPwr[i] = eccPwr[i - 1] * eo2;
  195.             chiPwr[i] = chiPwr[i - 1] * x2o2;
  196.         }

  198.         // Auxiliary quantities.
  199.         final double ao2rxx = aoR3 / (2. * XX);
  200.         double xmuarn       = ao2rxx * ao2rxx * gm / (X * R3);
  201.         double term         = 0.;

  203.         // Compute max power for a/R3 and e.
  204.         maxAR3Pow = 2;
  205.         maxEccPow = 0;
  206.         int n     = 2;
  207.         int m     = 2;
  208.         int nsmd2 = 0;

  210.         do {
  211.             // Upper bound for Tnm.
  212.             term =  xmuarn *
  213.                             (CombinatoricsUtils.factorialDouble(n + m) / (CombinatoricsUtils.factorialDouble(nsmd2) * CombinatoricsUtils.factorialDouble(nsmd2 + m))) *
  214.                             (CombinatoricsUtils.factorialDouble(n + m + 1) / (CombinatoricsUtils.factorialDouble(m) * CombinatoricsUtils.factorialDouble(n + 1))) *
  215.                             (CombinatoricsUtils.factorialDouble(n - m + 1) / CombinatoricsUtils.factorialDouble(n + 1)) *
  216.                             eccPwr[m] * UpperBounds.getDnl(XX, chiPwr[m], n + 2, m);

  218.             if (term < tol) {
  219.                 if (m == 0) {
  220.                     break;
  221.                 } else if (m < 2) {
  222.                     xmuarn *= ao2rxx;
  223.                     m = 0;
  224.                     n++;
  225.                     nsmd2++;
  226.                 } else {
  227.                     m -= 2;
  228.                     nsmd2++;
  229.                 }
  230.             } else {
  231.                 maxAR3Pow = n;
  232.                 maxEccPow = FastMath.max(m, maxEccPow);
  233.                 xmuarn *= ao2rxx;
  234.                 m++;
  235.                 n++;
  236.             }
  237.         } while (n < MAX_POWER);

  239.         maxEccPow = FastMath.min(maxAR3Pow, maxEccPow);

  241.         // allocate the array aoR3Pow
  242.         aoR3Pow = new double[maxAR3Pow + 1];

  244.         aoR3Pow[0] = 1.;
  245.         for (int i = 1; i <= maxAR3Pow; i++) {
  246.             aoR3Pow[i] = aoR3 * aoR3Pow[i - 1];
  247.         }

  249.         maxFreqF = maxAR3Pow + 1;
  250.         maxEccPowShort = MAX_ECCPOWER_SP;

  252.         Qns = CoefficientsFactory.computeQns(gamma, maxAR3Pow, FastMath.max(maxEccPow, maxEccPowShort));

  254.     }

  256.     /** Get A = sqrt(μ * a).
  257.      * @return A
  258.      */
  259.     public double getA() {
  260.         return A;
  261.     }

  263.     /** Get direction cosine α for central body.
  264.      * @return α
  265.      */
  266.     public double getAlpha() {
  267.         return alpha;
  268.     }

  270.     /** Get direction cosine β for central body.
  271.      * @return β
  272.      */
  273.     public double getBeta() {
  274.         return beta;
  275.     }

  277.     /** Get direction cosine γ for central body.
  278.      * @return γ
  279.      */
  280.     public double getGamma() {
  281.         return gamma;
  282.     }

  284.     /** Get B².
  285.      * @return B²
  286.      */
  287.     public double getBB() {
  288.         return BB;
  289.     }

  291.     /** Get B³.
  292.      * @return B³
  293.      */
  294.     public double getBBB() {
  295.         return BBB;
  296.     }

  298.     /** Get b = 1 / (1 + sqrt(1 - e²)) = 1 / (1 + B).
  299.      * @return b
  300.      */
  301.     public double getb() {
  302.         return b;
  303.     }

  305.     /** Get &Chi; = 1 / sqrt(1 - e²) = 1 / B.
  306.      * @return &Chi;
  307.      */
  308.     public double getX() {
  309.         return X;
  310.     }

  312.     /** Get m2aoA = -2 * a / A.
  313.      * @return m2aoA
  314.      */
  315.     public double getM2aoA() {
  316.         return m2aoA;
  317.     }

  319.     /** Get B / A.
  320.      * @return BoA
  321.      */
  322.     public double getBoA() {
  323.         return BoA;
  324.     }

  326.     /** Get ooAB = 1 / (A * B).
  327.      * @return ooAB
  328.      */
  329.     public double getOoAB() {
  330.         return ooAB;
  331.     }

  333.     /** Get mCo2AB = -C / 2AB.
  334.      * @return mCo2AB
  335.      */
  336.     public double getMCo2AB() {
  337.         return mCo2AB;
  338.     }

  340.     /** Get BoABpo = B / A(1 + B).
  341.      * @return BoABpo
  342.      */
  343.     public double getBoABpo() {
  344.         return BoABpo;
  345.     }

  347.     /** Get muoR3 = mu3 / R3.
  348.      * @return muoR3
  349.      */
  350.     public double getMuoR3() {
  351.         return muoR3;
  352.     }

  354.     /** Get hXXX = h * &Chi;³.
  355.      * @return hXXX
  356.      */
  357.     public double getHXXX() {
  358.         return hXXX;
  359.     }

  361.     /** Get kXXX = h * &Chi;³.
  362.      * @return kXXX
  363.      */
  364.     public double getKXXX() {
  365.         return kXXX;
  366.     }

  368.     /** Get the value of max power for a/R3 in the serie expansion.
  369.      * @return maxAR3Pow
  370.      */
  371.     public int getMaxAR3Pow() {
  372.         return maxAR3Pow;
  373.     }

  375.     /** Get the value of max power for e in the serie expansion.
  376.      * @return maxEccPow
  377.      */
  378.     public int getMaxEccPow() {
  379.         return maxEccPow;
  380.     }

  382.     /** Get the value of a / R3 up to power maxAR3Pow.
  383.      * @return aoR3Pow
  384.      */
  385.     public double[] getAoR3Pow() {
  386.         return aoR3Pow;
  387.     }

  389.    /** Get the value of max frequency of F.
  390.      * @return maxFreqF
  391.      */
  392.     public int getMaxFreqF() {
  393.         return maxFreqF;
  394.     }

  396.     /** Get the Keplerian mean motion.
  397.      * <p>The Keplerian mean motion is computed directly from semi major axis
  398.      * and central acceleration constant.</p>
  399.      * @return Keplerian mean motion in radians per second
  400.      */
  401.     public double getMeanMotion() {
  402.         return motion;
  403.     }

  405.     /** Get the value of Qns coefficients.
  406.      * @return Qns
  407.      */
  408.     public double[][] getQns() {
  409.         return Qns;
  410.     }

  412. }
Created with JaCoCo 0.8.6.202009150832