FieldDSSTTesseralContext.java

  1. /* Copyright 2002-2023 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.CalculusFieldElement;
  20. import org.hipparchus.Field;
  21. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  22. import org.hipparchus.util.FastMath;
  23. import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
  24. import org.orekit.frames.FieldStaticTransform;
  25. import org.orekit.frames.Frame;
  26. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  27. import org.orekit.time.AbsoluteDate;

  29. /**
  30.  * This class is a container for the common "field" 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.  * @param <T> type of the field elements
  38.  */
  39. public class FieldDSSTTesseralContext<T extends CalculusFieldElement<T>> extends FieldForceModelContext<T> {

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

  56.     /** A = sqrt(μ * a). */
  57.     private T A;

  59.     // Common factors for potential computation
  60.     /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
  61.     private T chi;

  63.     /** &Chi;². */
  64.     private T chi2;

  66.     /** Central body rotation angle θ. */
  67.     private T theta;

  69.     // Common factors from equinoctial coefficients
  70.     /** 2 * a / A .*/
  71.     private T ax2oA;

  73.     /** 1 / (A * B) .*/
  74.     private T ooAB;

  76.     /** B / A .*/
  77.     private T BoA;

  79.     /** B / (A * (1 + B)) .*/
  80.     private T BoABpo;

  82.     /** C / (2 * A * B) .*/
  83.     private T Co2AB;

  85.     /** μ / a .*/
  86.     private T moa;

  88.     /** R / a .*/
  89.     private T roa;

  91.     /** ecc². */
  92.     private T e2;

  94.     /** Keplerian mean motion. */
  95.     private T n;

  97.     /** Keplerian period. */
  98.     private T period;

  100.     /** Ratio of satellite period to central body rotation period. */
  101.     private T ratio;

  103.     /**
  104.      * Simple constructor.
  105.      *
  106.      * @param auxiliaryElements auxiliary elements related to the current orbit
  107.      * @param centralBodyFrame rotating body frame
  108.      * @param provider provider for spherical harmonics
  109.      * @param maxFrequencyShortPeriodics maximum value for j
  110.      * @param bodyPeriod central body rotation period (seconds)
  111.      * @param parameters values of the force model parameters (only 1 values
  112.      * for each parameters corresponding to state date) obtained by calling
  113.      * the extract parameter method {@link #extractParameters(double[], AbsoluteDate)}
  114.      * to selected the right value for state date or by getting the parameters for a specific date
  115.      */
  116.     FieldDSSTTesseralContext(final FieldAuxiliaryElements<T> auxiliaryElements,
  117.                                     final Frame centralBodyFrame,
  118.                                     final UnnormalizedSphericalHarmonicsProvider provider,
  119.                                     final int maxFrequencyShortPeriodics,
  120.                                     final double bodyPeriod,
  121.                                     final T[] parameters) {

  123.         super(auxiliaryElements);

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

  128.         final T mu = parameters[0];

  130.         // Keplerian mean motion
  131.         final T absA = FastMath.abs(auxiliaryElements.getSma());
  132.         n = FastMath.sqrt(mu.divide(absA)).divide(absA);

  134.         // Keplerian period
  135.         final T a = auxiliaryElements.getSma();
  136.         period = (a.getReal() < 0) ? zero.add(Double.POSITIVE_INFINITY) : a.multiply(a.getPi().multiply(2.0)).multiply(a.divide(mu).sqrt());

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

  140.         // Eccentricity square
  141.         e2 = auxiliaryElements.getEcc().multiply(auxiliaryElements.getEcc());

  143.         // Central body rotation angle from equation 2.7.1-(3)(4).
  144.         final FieldStaticTransform<T> t = centralBodyFrame.getStaticTransformTo(auxiliaryElements.getFrame(), auxiliaryElements.getDate());
  145.         final FieldVector3D<T> xB = t.transformVector(FieldVector3D.getPlusI(field));
  146.         final FieldVector3D<T> yB = t.transformVector(FieldVector3D.getPlusJ(field));
  147.         theta = FastMath.atan2(auxiliaryElements.getVectorF().dotProduct(yB).negate().add((auxiliaryElements.getVectorG().dotProduct(xB)).multiply(I)),
  148.                                auxiliaryElements.getVectorF().dotProduct(xB).add(auxiliaryElements.getVectorG().dotProduct(yB).multiply(I)));

  150.         // Common factors from equinoctial coefficients
  151.         // 2 * a / A
  152.         ax2oA  = auxiliaryElements.getSma().divide(A).multiply(2.);
  153.         // B / A
  154.         BoA    = auxiliaryElements.getB().divide(A);
  155.         // 1 / AB
  156.         ooAB   = A.multiply(auxiliaryElements.getB()).reciprocal();
  157.         // C / 2AB
  158.         Co2AB  = auxiliaryElements.getC().multiply(ooAB).divide(2.);
  159.         // B / (A * (1 + B))
  160.         BoABpo = BoA.divide(auxiliaryElements.getB().add(1.));
  161.         // &mu / a
  162.         moa    = mu.divide(auxiliaryElements.getSma());
  163.         // R / a
  164.         roa    = auxiliaryElements.getSma().divide(provider.getAe()).reciprocal();

  166.         // &Chi; = 1 / B
  167.         chi  = auxiliaryElements.getB().reciprocal();
  168.         chi2 = chi.multiply(chi);

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

  173.     }

  175.     /** Get ecc².
  176.      * @return e2
  177.      */
  178.     public T getE2() {
  179.         return e2;
  180.     }

  182.     /** Get Central body rotation angle θ.
  183.      * @return theta
  184.      */
  185.     public T getTheta() {
  186.         return theta;
  187.     }

  189.     /** Get ax2oA = 2 * a / A .
  190.      * @return ax2oA
  191.      */
  192.     public T getAx2oA() {
  193.         return ax2oA;
  194.     }

  196.     /** Get &Chi; = 1 / sqrt(1 - e²) = 1 / B.
  197.      * @return chi
  198.      */
  199.     public T getChi() {
  200.         return chi;
  201.     }

  203.     /** Get &Chi;².
  204.      * @return chi2
  205.      */
  206.     public T getChi2() {
  207.         return chi2;
  208.     }

  210.     /** Get B / A.
  211.      * @return BoA
  212.      */
  213.     public T getBoA() {
  214.         return BoA;
  215.     }

  217.     /** Get ooAB = 1 / (A * B).
  218.      * @return ooAB
  219.      */
  220.     public T getOoAB() {
  221.         return ooAB;
  222.     }

  224.     /** Get Co2AB = C / 2AB.
  225.      * @return Co2AB
  226.      */
  227.     public T getCo2AB() {
  228.         return Co2AB;
  229.     }

  231.     /** Get BoABpo = B / A(1 + B).
  232.      * @return BoABpo
  233.      */
  234.     public T getBoABpo() {
  235.         return BoABpo;
  236.     }

  238.     /** Get μ / a .
  239.      * @return moa
  240.      */
  241.     public T getMoa() {
  242.         return moa;
  243.     }

  245.     /** Get roa = R / a.
  246.      * @return roa
  247.      */
  248.     public T getRoa() {
  249.         return roa;
  250.     }

  252.     /** Get the Keplerian period.
  253.      * <p>The Keplerian period is computed directly from semi major axis
  254.      * and central acceleration constant.</p>
  255.      * @return Keplerian period in seconds, or positive infinity for hyperbolic orbits
  256.      */
  257.     public T getOrbitPeriod() {
  258.         return period;
  259.     }

  261.     /** Get the Keplerian mean motion.
  262.      * <p>The Keplerian mean motion is computed directly from semi major axis
  263.      * and central acceleration constant.</p>
  264.      * @return Keplerian mean motion in radians per second
  265.      */
  266.     public T getMeanMotion() {
  267.         return n;
  268.     }

  270.     /** Get the ratio of satellite period to central body rotation period.
  271.      * @return ratio
  272.      */
  273.     public T getRatio() {
  274.         return ratio;
  275.     }

  277. }
Created with JaCoCo 0.8.11.202310140853