FieldAuxiliaryElements.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.utilities;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.RealFieldElement;
  21. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.MathUtils;
  24. import org.orekit.frames.Frame;
  25. import org.orekit.orbits.FieldOrbit;
  26. import org.orekit.time.FieldAbsoluteDate;

  28. /** Container class for common parameters used by all DSST forces.
  29.  *  <p>
  30.  *  Most of them are defined in Danielson paper at § 2.1.
  31.  *  </p>
  32.  */
  33. public class FieldAuxiliaryElements<T extends RealFieldElement<T>> {

  35.     /** \(2\pi\) . */
  36.     public static final double TWO_PI = 2 * FastMath.PI;

  38.     /** Orbit date. */
  39.     private final FieldAbsoluteDate<T> date;

  41.     /** Orbit frame. */
  42.     private final Frame frame;

  44.     /** Eccentricity. */
  45.     private final T ecc;

  47.     /** Keplerian mean motion. */
  48.     private final T n;

  50.     /** Keplerian period. */
  51.     private final T period;

  53.     /** Semi-major axis. */
  54.     private final T sma;

  56.     /** x component of eccentricity vector. */
  57.     private final T k;

  59.     /** y component of eccentricity vector. */
  60.     private final T h;

  62.     /** x component of inclination vector. */
  63.     private final T q;

  65.     /** y component of inclination vector. */
  66.     private final T p;

  68.     /** Mean longitude. */
  69.     private final T lm;

  71.     /** True longitude. */
  72.     private final T lv;

  74.     /** Eccentric longitude. */
  75.     private final T le;

  77.     /** Retrograde factor I.
  78.      *  <p>
  79.      *  DSST model needs equinoctial orbit as internal representation.
  80.      *  Classical equinoctial elements have discontinuities when inclination
  81.      *  is close to zero. In this representation, I = +1. <br>
  82.      *  To avoid this discontinuity, another representation exists and equinoctial
  83.      *  elements can be expressed in a different way, called "retrograde" orbit.
  84.      *  This implies I = -1. <br>
  85.      *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  86.      *  But for the sake of consistency with the theory, the retrograde factor
  87.      *  has been kept in the formulas.
  88.      *  </p>
  89.      */
  90.     private final int    I;

  92.     /** B = sqrt(1 - h² - k²). */
  93.     private final T B;

  95.     /** C = 1 + p² + q². */
  96.     private final T C;

  98.     /** Equinoctial frame f vector. */
  99.     private final FieldVector3D<T> f;

  101.     /** Equinoctial frame g vector. */
  102.     private final FieldVector3D<T> g;

  104.     /** Equinoctial frame w vector. */
  105.     private final FieldVector3D<T> w;

  107.     /** Direction cosine α. */
  108.     private final T alpha;

  110.     /** Direction cosine β. */
  111.     private final T beta;

  113.     /** Direction cosine γ. */
  114.     private final T gamma;

  116.     /** Simple constructor.
  117.      * @param orbit related mean orbit for auxiliary elements
  118.      * @param retrogradeFactor retrograde factor I [Eq. 2.1.2-(2)]
  119.      */
  120.     public FieldAuxiliaryElements(final FieldOrbit<T> orbit, final int retrogradeFactor) {

  122.         final Field<T> field = orbit.getDate().getField();
  123.         final T zero = field.getZero();
  124.         final T pi = zero.add(FastMath.PI);

  126.         // Date of the orbit
  127.         date = orbit.getDate();

  129.         // Orbit definition frame
  130.         frame = orbit.getFrame();

  132.         // Eccentricity
  133.         ecc = orbit.getE();

  135.         // Keplerian mean motion
  136.         n = orbit.getKeplerianMeanMotion();

  138.         // Keplerian period
  139.         period = orbit.getKeplerianPeriod();

  141.         // Equinoctial elements [Eq. 2.1.2-(1)]
  142.         sma = orbit.getA();
  143.         k   = orbit.getEquinoctialEx();
  144.         h   = orbit.getEquinoctialEy();
  145.         q   = orbit.getHx();
  146.         p   = orbit.getHy();
  147.         lm  = MathUtils.normalizeAngle(orbit.getLM(), pi);
  148.         lv  = MathUtils.normalizeAngle(orbit.getLv(), pi);
  149.         le  = MathUtils.normalizeAngle(orbit.getLE(), pi);

  151.         // Retrograde factor [Eq. 2.1.2-(2)]
  152.         I = retrogradeFactor;

  154.         final T k2 = k.multiply(k);
  155.         final T h2 = h.multiply(h);
  156.         final T q2 = q.multiply(q);
  157.         final T p2 = p.multiply(p);

  159.         // A, B, C parameters [Eq. 2.1.6-(1)]
  160.         B = FastMath.sqrt(k2.add(h2).negate().add(1.));
  161.         C = q2.add(p2).add(1);

  163.         // Equinoctial reference frame [Eq. 2.1.4-(1)]
  164.         final T ooC = C.reciprocal();
  165.         final T px2 = p.multiply(2.);
  166.         final T qx2 = q.multiply(2.);
  167.         final T pq2 = px2.multiply(q);
  168.         f = new FieldVector3D<>(ooC, new FieldVector3D<>(p2.negate().add(1.).add(q2), pq2, px2.multiply(I).negate()));
  169.         g = new FieldVector3D<>(ooC, new FieldVector3D<>(pq2.multiply(I), (p2.add(1.).subtract(q2)).multiply(I), qx2));
  170.         w = new FieldVector3D<>(ooC, new FieldVector3D<>(px2, qx2.negate(), (p2.add(q2).negate().add(1.)).multiply(I)));

  172.         // Direction cosines for central body [Eq. 2.1.9-(1)]
  173.         alpha = (T) f.getZ();
  174.         beta  = (T) g.getZ();
  175.         gamma = (T) w.getZ();
  176.     }

  178.     /**
  179.      * Normalize an angle in a 2&pi; wide interval around a center value.
  180.      * <p>This method has three main uses:</p>
  181.      * <ul>
  182.      *   <li>normalize an angle between 0 and 2&pi;:<br>
  183.      *       {@code a = MathUtils.normalizeAngle(a, FastMath.PI);}</li>
  184.      *   <li>normalize an angle between -&pi; and +&pi;<br>
  185.      *       {@code a = MathUtils.normalizeAngle(a, 0.0);}</li>
  186.      *   <li>compute the angle between two defining angular positions:<br>
  187.      *       {@code angle = MathUtils.normalizeAngle(end, start) - start;}</li>
  188.      * </ul>
  189.      * <p>Note that due to numerical accuracy and since &pi; cannot be represented
  190.      * exactly, the result interval is <em>closed</em>, it cannot be half-closed
  191.      * as would be more satisfactory in a purely mathematical view.</p>
  192.      * @param a angle to normalize
  193.      * @param center center of the desired 2&pi; interval for the result
  194.      * @return a-2k&pi; with integer k and center-&pi; &lt;= a-2k&pi; &lt;= center+&pi;
  195.      * @deprecated replaced by {@link MathUtils#normalizeAngle(RealFieldElement, RealFieldElement)}
  196.      */
  197.     @Deprecated
  198.     public T normalizeAngle(final T a, final T center) {
  199.         return a.subtract(FastMath.floor((a.add(FastMath.PI).subtract(center)).divide(TWO_PI)).multiply(TWO_PI));
  200.     }

  202.     /** Get the date of the orbit.
  203.      * @return the date
  204.      */
  205.     public FieldAbsoluteDate<T> getDate() {
  206.         return date;
  207.     }

  209.     /** Get the definition frame of the orbit.
  210.      * @return the definition frame
  211.      */
  212.     public Frame getFrame() {
  213.         return frame;
  214.     }

  216.     /** Get the eccentricity.
  217.      * @return ecc
  218.      */
  219.     public T getEcc() {
  220.         return ecc;
  221.     }

  223.     /** Get the Keplerian mean motion.
  224.      * @return n
  225.      */
  226.     public T getMeanMotion() {
  227.         return n;
  228.     }

  230.     /** Get the Keplerian period.
  231.      * @return period
  232.      */
  233.     public T getKeplerianPeriod() {
  234.         return period;
  235.     }

  237.     /** Get the semi-major axis.
  238.      * @return the semi-major axis a
  239.      */
  240.     public T getSma() {
  241.         return sma;
  242.     }

  244.     /** Get the x component of eccentricity vector.
  245.      * <p>
  246.      * This element called k in DSST corresponds to ex
  247.      * for the {@link org.orekit.orbits.EquinoctialOrbit}
  248.      * </p>
  249.      * @return k
  250.      */
  251.     public T getK() {
  252.         return k;
  253.     }

  255.     /** Get the y component of eccentricity vector.
  256.      * <p>
  257.      * This element called h in DSST corresponds to ey
  258.      * for the {@link org.orekit.orbits.EquinoctialOrbit}
  259.      * </p>
  260.      * @return h
  261.      */
  262.     public T getH() {
  263.         return h;
  264.     }

  266.     /** Get the x component of inclination vector.
  267.      * <p>
  268.      * This element called q in DSST corresponds to hx
  269.      * for the {@link org.orekit.orbits.EquinoctialOrbit}
  270.      * </p>
  271.      * @return q
  272.      */
  273.     public T getQ() {
  274.         return q;
  275.     }

  277.     /** Get the y component of inclination vector.
  278.      * <p>
  279.      * This element called p in DSST corresponds to hy
  280.      * for the {@link org.orekit.orbits.EquinoctialOrbit}
  281.      * </p>
  282.      * @return p
  283.      */
  284.     public T getP() {
  285.         return p;
  286.     }

  288.     /** Get the mean longitude.
  289.      * @return lm
  290.      */
  291.     public T getLM() {
  292.         return lm;
  293.     }

  295.     /** Get the true longitude.
  296.      * @return lv
  297.      */
  298.     public T getLv() {
  299.         return lv;
  300.     }

  302.     /** Get the eccentric longitude.
  303.      * @return le
  304.      */
  305.     public T getLe() {
  306.         return le;
  307.     }

  309.     /** Get the retrograde factor.
  310.      * @return the retrograde factor I
  311.      */
  312.     public int getRetrogradeFactor() {
  313.         return I;
  314.     }

  316.     /** Get B = sqrt(1 - e²).
  317.      * @return B
  318.      */
  319.     public T getB() {
  320.         return B;
  321.     }

  323.     /** Get C = 1 + p² + q².
  324.      * @return C
  325.      */
  326.     public T getC() {
  327.         return C;
  328.     }

  330.     /** Get equinoctial frame vector f.
  331.      * @return f vector
  332.      */
  333.     public FieldVector3D<T> getVectorF() {
  334.         return f;
  335.     }

  337.     /** Get equinoctial frame vector g.
  338.      * @return g vector
  339.      */
  340.     public FieldVector3D<T> getVectorG() {
  341.         return g;
  342.     }

  344.     /** Get equinoctial frame vector w.
  345.      * @return w vector
  346.      */
  347.     public FieldVector3D<T> getVectorW() {
  348.         return w;
  349.     }

  351.     /** Get direction cosine α for central body.
  352.      * @return α
  353.      */
  354.     public T getAlpha() {
  355.         return alpha;
  356.     }

  358.     /** Get direction cosine β for central body.
  359.      * @return β
  360.      */
  361.     public T getBeta() {
  362.         return beta;
  363.     }

  365.     /** Get direction cosine γ for central body.
  366.      * @return γ
  367.      */
  368.     public T getGamma() {
  369.         return gamma;
  370.     }

  372. }
Created with JaCoCo 0.8.6.202009150832