FieldDSSTThirdBodyContext.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.Field;
  20. import org.hipparchus.RealFieldElement;
  21. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  22. import org.hipparchus.util.CombinatoricsUtils;
  23. import org.hipparchus.util.FastMath;
  24. import org.hipparchus.util.MathArrays;
  25. import org.orekit.bodies.CelestialBody;
  26. import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
  27. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  28. import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;

  30. /**
  31.  * This class is a container for the common "field" parameters used in {@link DSSTThirdBody}.
  32.  * <p>
  33.  * It performs parameters initialization at each integration step for the third
  34.  * body attraction perturbation.
  35.  * <p>
  36.  * @author Bryan Cazabonne
  37.  * @since 10.0
  38.  */
  39. class FieldDSSTThirdBodyContext<T extends RealFieldElement <T>> extends FieldForceModelContext<T> {

  41.     /** Max power for summation. */
  42.     private static final int    MAX_POWER = 22;

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

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

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

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

  56.     /** Max power for e in the serie expansion. */
  57.     private int maxEccPow;

  59.     /** a / R3 up to power maxAR3Pow. */
  60.     private T[] aoR3Pow;

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

  65.     /** Max frequency of F. */
  66.     private int maxFreqF;

  68.     /** Qns coefficients. */
  69.     private T[][] Qns;

  71.     /** Standard gravitational parameter μ for the body in m³/s². */
  72.     private final T gm;

  74.     /** Distance from center of mass of the central body to the 3rd body. */
  75.     private T R3;

  77.     /** A = sqrt(μ * a). */
  78.     private final T A;

  80.     // Direction cosines of the symmetry axis
  81.     /** α. */
  82.     private final T alpha;
  83.     /** β. */
  84.     private final T beta;
  85.     /** γ. */
  86.     private final T gamma;

  88.     /** B². */
  89.     private final T BB;
  90.     /** B³. */
  91.     private final T BBB;

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

  110.     /** mu3 / R3. */
  111.     private final T muoR3;

  113.     /** b = 1 / (1 + sqrt(1 - e²)) = 1 / (1 + B).*/
  114.     private final T b;

  116.     /** h * &Chi;³. */
  117.     private final T hXXX;
  118.     /** k * &Chi;³. */
  119.     private final T kXXX;

  121.     /** Keplerian mean motion. */
  122.     private final T motion;

  124.     /**
  125.      * Simple constructor.
  126.      *
  127.      * @param auxiliaryElements auxiliary elements related to the current orbit
  128.      * @param thirdBody body the 3rd body to consider
  129.      * @param parameters values of the force model parameters
  130.      */
  131.     FieldDSSTThirdBodyContext(final FieldAuxiliaryElements<T> auxiliaryElements,
  132.                                      final CelestialBody thirdBody,
  133.                                      final T[] parameters) {

  135.         super(auxiliaryElements);

  137.         // Field for array building
  138.         final Field<T> field = auxiliaryElements.getDate().getField();
  139.         final T zero = field.getZero();

  141.         final T mu = parameters[1];
  142.         A = FastMath.sqrt(mu.multiply(auxiliaryElements.getSma()));

  144.         this.gm = parameters[0];

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

  150.         // Distance from center of mass of the central body to the 3rd body
  151.         final FieldVector3D<T> bodyPos = thirdBody.getPVCoordinates(auxiliaryElements.getDate(), auxiliaryElements.getFrame()).getPosition();
  152.         R3 = bodyPos.getNorm();

  154.         // Direction cosines
  155.         final FieldVector3D<T> bodyDir = bodyPos.normalize();
  156.         alpha = (T) bodyDir.dotProduct(auxiliaryElements.getVectorF());
  157.         beta  = (T) bodyDir.dotProduct(auxiliaryElements.getVectorG());
  158.         gamma = (T) bodyDir.dotProduct(auxiliaryElements.getVectorW());

  160.         //&Chi;<sup>-2</sup>.
  161.         BB = auxiliaryElements.getB().multiply(auxiliaryElements.getB());
  162.         //&Chi;<sup>-3</sup>.
  163.         BBB = BB.multiply(auxiliaryElements.getB());

  165.         //b = 1 / (1 + B)
  166.         b = auxiliaryElements.getB().add(1.).reciprocal();

  168.         // &Chi;
  169.         X = auxiliaryElements.getB().reciprocal();
  170.         XX = X.multiply(X);
  171.         XXX = X.multiply(XX);

  173.         // -2 * a / A
  174.         m2aoA = auxiliaryElements.getSma().multiply(-2.).divide(A);
  175.         // B / A
  176.         BoA = auxiliaryElements.getB().divide(A);
  177.         // 1 / AB
  178.         ooAB = (A.multiply(auxiliaryElements.getB())).reciprocal();
  179.         // -C / 2AB
  180.         mCo2AB = auxiliaryElements.getC().multiply(ooAB).divide(2.).negate();
  181.         // B / A(1 + B)
  182.         BoABpo = BoA.divide(auxiliaryElements.getB().add(1.));
  183.         // mu3 / R3
  184.         muoR3 = R3.divide(gm).reciprocal();
  185.         //h * &Chi;³
  186.         hXXX = XXX.multiply(auxiliaryElements.getH());
  187.         //k * &Chi;³
  188.         kXXX = XXX.multiply(auxiliaryElements.getK());

  190.         // Truncation tolerance.
  191.         final T aoR3 = auxiliaryElements.getSma().divide(R3);
  192.         final double tol = ( aoR3.getReal() > .3 || (aoR3.getReal() > .15  && auxiliaryElements.getEcc().getReal() > .25) ) ? BIG_TRUNCATION_TOLERANCE : SMALL_TRUNCATION_TOLERANCE;

  194.         // Utilities for truncation
  195.         // Set a lower bound for eccentricity
  196.         final T eo2  = FastMath.max(zero.add(0.0025), auxiliaryElements.getEcc().divide(2.));
  197.         final T x2o2 = XX.divide(2.);
  198.         final T[] eccPwr = MathArrays.buildArray(field, MAX_POWER);
  199.         final T[] chiPwr = MathArrays.buildArray(field, MAX_POWER);
  200.         eccPwr[0] = zero.add(1.);
  201.         chiPwr[0] = X;
  202.         for (int i = 1; i < MAX_POWER; i++) {
  203.             eccPwr[i] = eccPwr[i - 1].multiply(eo2);
  204.             chiPwr[i] = chiPwr[i - 1].multiply(x2o2);
  205.         }

  207.         // Auxiliary quantities.
  208.         final T ao2rxx = aoR3.divide(XX.multiply(2.));
  209.         T xmuarn       = ao2rxx.multiply(ao2rxx).multiply(gm).divide(X.multiply(R3));
  210.         T term         = zero;

  212.         // Compute max power for a/R3 and e.
  213.         maxAR3Pow = 2;
  214.         maxEccPow = 0;
  215.         int n     = 2;
  216.         int m     = 2;
  217.         int nsmd2 = 0;

  219.         do {
  220.             term =  xmuarn.multiply((CombinatoricsUtils.factorialDouble(n + m) / (CombinatoricsUtils.factorialDouble(nsmd2) * CombinatoricsUtils.factorialDouble(nsmd2 + m))) *
  221.                             (CombinatoricsUtils.factorialDouble(n + m + 1) / (CombinatoricsUtils.factorialDouble(m) * CombinatoricsUtils.factorialDouble(n + 1))) *
  222.                             (CombinatoricsUtils.factorialDouble(n - m + 1) / CombinatoricsUtils.factorialDouble(n + 1))).
  223.                             multiply(eccPwr[m]).multiply(UpperBounds.getDnl(XX, chiPwr[m], n + 2, m));

  225.             if (term.getReal() < tol) {
  226.                 if (m == 0) {
  227.                     break;
  228.                 } else if (m < 2) {
  229.                     xmuarn = xmuarn.multiply(ao2rxx);
  230.                     m = 0;
  231.                     n++;
  232.                     nsmd2++;
  233.                 } else {
  234.                     m -= 2;
  235.                     nsmd2++;
  236.                 }
  237.             } else {
  238.                 maxAR3Pow = n;
  239.                 maxEccPow = FastMath.max(m, maxEccPow);
  240.                 xmuarn = xmuarn.multiply(ao2rxx);
  241.                 m++;
  242.                 n++;
  243.             }
  244.         } while (n < MAX_POWER);

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

  248.         // allocate the array aoR3Pow
  249.         aoR3Pow = MathArrays.buildArray(field, maxAR3Pow + 1);

  251.         aoR3Pow[0] = field.getOne();
  252.         for (int i = 1; i <= maxAR3Pow; i++) {
  253.             aoR3Pow[i] = aoR3.multiply(aoR3Pow[i - 1]);
  254.         }

  256.         maxFreqF = maxAR3Pow + 1;
  257.         maxEccPowShort = MAX_ECCPOWER_SP;

  259.         Qns = CoefficientsFactory.computeQns(gamma, maxAR3Pow, FastMath.max(maxEccPow, maxEccPowShort));
  260.     }

  262.     /** Get A = sqrt(μ * a).
  263.      * @return A
  264.      */
  265.     public T getA() {
  266.         return A;
  267.     }

  269.     /** Get direction cosine α for central body.
  270.      * @return α
  271.      */
  272.     public T getAlpha() {
  273.         return alpha;
  274.     }

  276.     /** Get direction cosine β for central body.
  277.      * @return β
  278.      */
  279.     public T getBeta() {
  280.         return beta;
  281.     }

  283.     /** Get direction cosine γ for central body.
  284.      * @return γ
  285.      */
  286.     public T getGamma() {
  287.         return gamma;
  288.     }

  290.     /** Get B².
  291.      * @return B²
  292.      */
  293.     public T getBB() {
  294.         return BB;
  295.     }

  297.     /** Get B³.
  298.      * @return B³
  299.      */
  300.     public T getBBB() {
  301.         return BBB;
  302.     }

  304.     /** Get b = 1 / (1 + sqrt(1 - e²)) = 1 / (1 + B).
  305.      * @return b
  306.      */
  307.     public T getb() {
  308.         return b;
  309.     }

  311.     /** Get &Chi; = 1 / sqrt(1 - e²) = 1 / B.
  312.      * @return &Chi;
  313.      */
  314.     public T getX() {
  315.         return X;
  316.     }

  318.     /** Get m2aoA = -2 * a / A.
  319.      * @return m2aoA
  320.      */
  321.     public T getM2aoA() {
  322.         return m2aoA;
  323.     }

  325.     /** Get B / A.
  326.      * @return BoA
  327.      */
  328.     public T getBoA() {
  329.         return BoA;
  330.     }

  332.     /** Get ooAB = 1 / (A * B).
  333.      * @return ooAB
  334.      */
  335.     public T getOoAB() {
  336.         return ooAB;
  337.     }

  339.     /** Get mCo2AB = -C / 2AB.
  340.      * @return mCo2AB
  341.      */
  342.     public T getMCo2AB() {
  343.         return mCo2AB;
  344.     }

  346.     /** Get BoABpo = B / A(1 + B).
  347.      * @return BoABpo
  348.      */
  349.     public T getBoABpo() {
  350.         return BoABpo;
  351.     }

  353.     /** Get muoR3 = mu3 / R3.
  354.      * @return muoR3
  355.      */
  356.     public T getMuoR3() {
  357.         return muoR3;
  358.     }

  360.     /** Get hXXX = h * &Chi;³.
  361.      * @return hXXX
  362.      */
  363.     public T getHXXX() {
  364.         return hXXX;
  365.     }

  367.     /** Get kXXX = h * &Chi;³.
  368.      * @return kXXX
  369.      */
  370.     public T getKXXX() {
  371.         return kXXX;
  372.     }

  374.     /** Get the value of max power for a/R3 in the serie expansion.
  375.      * @return maxAR3Pow
  376.      */
  377.     public int getMaxAR3Pow() {
  378.         return maxAR3Pow;
  379.     }

  381.     /** Get the value of max power for e in the serie expansion.
  382.      * @return maxEccPow
  383.      */
  384.     public int getMaxEccPow() {
  385.         return maxEccPow;
  386.     }

  388.     /** Get the value of a / R3 up to power maxAR3Pow.
  389.      * @return aoR3Pow
  390.      */
  391.     public T[] getAoR3Pow() {
  392.         return aoR3Pow;
  393.     }

  395.    /** Get the value of max frequency of F.
  396.      * @return maxFreqF
  397.      */
  398.     public int getMaxFreqF() {
  399.         return maxFreqF;
  400.     }

  402.     /** Get the Keplerian mean motion.
  403.      * <p>The Keplerian mean motion is computed directly from semi major axis
  404.      * and central acceleration constant.</p>
  405.      * @return Keplerian mean motion in radians per second
  406.      */
  407.     public T getMeanMotion() {
  408.         return motion;
  409.     }

  411.     /** Get the value of Qns coefficients.
  412.      * @return Qns
  413.      */
  414.     public T[][] getQns() {
  415.         return Qns;
  416.     }

  418. }
Created with JaCoCo 0.8.6.202009150832