FieldAuxiliaryElements.java

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

  18. import org.hipparchus.CalculusFieldElement;
  19. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  20. import org.hipparchus.util.FastMath;
  21. import org.hipparchus.util.MathUtils;
  22. import org.orekit.frames.Frame;
  23. import org.orekit.orbits.FieldOrbit;
  24. import org.orekit.propagation.semianalytical.dsst.forces.FieldDSSTGravityContext;
  25. import org.orekit.time.FieldAbsoluteDate;

  26. /** Container class for common parameters used by all DSST forces.
  27.  *  <p>
  28.  *  Most of them are defined in Danielson paper at § 2.1.
  29.  *  </p>
  30.  *  @author Bryan Cazabonne
  31.  * @param <T> type of the field elements
  32.  */
  33. public class FieldAuxiliaryElements<T extends CalculusFieldElement<T>> {

  34.     /** Orbit. */
  35.     private final FieldOrbit<T> orbit;

  36.     /** Orbit date. */
  37.     private final FieldAbsoluteDate<T> date;

  38.     /** Orbit frame. */
  39.     private final Frame frame;

  40.     /** Eccentricity. */
  41.     private final T ecc;

  42.     /** Keplerian mean motion. */
  43.     private final T n;

  44.     /** Keplerian period. */
  45.     private final T period;

  46.     /** Semi-major axis. */
  47.     private final T sma;

  48.     /** x component of eccentricity vector. */
  49.     private final T k;

  50.     /** y component of eccentricity vector. */
  51.     private final T h;

  52.     /** x component of inclination vector. */
  53.     private final T q;

  54.     /** y component of inclination vector. */
  55.     private final T p;

  56.     /** Mean longitude. */
  57.     private final T lm;

  58.     /** True longitude. */
  59.     private final T lv;

  60.     /** Eccentric longitude. */
  61.     private final T le;

  62.     /** Retrograde factor I.
  63.      *  <p>
  64.      *  DSST model needs equinoctial orbit as internal representation.
  65.      *  Classical equinoctial elements have discontinuities when inclination
  66.      *  is close to zero. In this representation, I = +1. <br>
  67.      *  To avoid this discontinuity, another representation exists and equinoctial
  68.      *  elements can be expressed in a different way, called "retrograde" orbit.
  69.      *  This implies I = -1. <br>
  70.      *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  71.      *  But for the sake of consistency with the theory, the retrograde factor
  72.      *  has been kept in the formulas.
  73.      *  </p>
  74.      */
  75.     private final int    I;

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

  78.     /** C = 1 + p² + q². */
  79.     private final T C;

  80.     /** Equinoctial frame f vector. */
  81.     private final FieldVector3D<T> f;

  82.     /** Equinoctial frame g vector. */
  83.     private final FieldVector3D<T> g;

  84.     /** Equinoctial frame w vector. */
  85.     private final FieldVector3D<T> w;

  86.     /** Simple constructor.
  87.      * @param orbit related mean orbit for auxiliary elements
  88.      * @param retrogradeFactor retrograde factor I [Eq. 2.1.2-(2)]
  89.      */
  90.     public FieldAuxiliaryElements(final FieldOrbit<T> orbit, final int retrogradeFactor) {

  91.         final T pi = orbit.getDate().getField().getZero().getPi();

  92.         // Orbit
  93.         this.orbit = orbit;

  94.         // Date of the orbit
  95.         date = orbit.getDate();

  96.         // Orbit definition frame
  97.         frame = orbit.getFrame();

  98.         // Eccentricity
  99.         ecc = orbit.getE();

  100.         // Keplerian mean motion
  101.         n = orbit.getKeplerianMeanMotion();

  102.         // Keplerian period
  103.         period = orbit.getKeplerianPeriod();

  104.         // Equinoctial elements [Eq. 2.1.2-(1)]
  105.         sma = orbit.getA();
  106.         k   = orbit.getEquinoctialEx();
  107.         h   = orbit.getEquinoctialEy();
  108.         q   = orbit.getHx();
  109.         p   = orbit.getHy();
  110.         lm  = MathUtils.normalizeAngle(orbit.getLM(), pi);
  111.         lv  = MathUtils.normalizeAngle(orbit.getLv(), pi);
  112.         le  = MathUtils.normalizeAngle(orbit.getLE(), pi);

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

  115.         final T k2 = k.multiply(k);
  116.         final T h2 = h.multiply(h);
  117.         final T q2 = q.multiply(q);
  118.         final T p2 = p.multiply(p);

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

  122.         // Equinoctial reference frame [Eq. 2.1.4-(1)]
  123.         final T ooC = C.reciprocal();
  124.         final T px2 = p.multiply(2.);
  125.         final T qx2 = q.multiply(2.);
  126.         final T pq2 = px2.multiply(q);
  127.         f = new FieldVector3D<>(ooC, new FieldVector3D<>(p2.negate().add(1.).add(q2), pq2, px2.multiply(I).negate()));
  128.         g = new FieldVector3D<>(ooC, new FieldVector3D<>(pq2.multiply(I), (p2.add(1.).subtract(q2)).multiply(I), qx2));
  129.         w = new FieldVector3D<>(ooC, new FieldVector3D<>(px2, qx2.negate(), (p2.add(q2).negate().add(1.)).multiply(I)));
  130.     }

  131.     /** Get the orbit.
  132.      * @return the orbit
  133.      */
  134.     public FieldOrbit<T> getOrbit() {
  135.         return orbit;
  136.     }

  137.     /** Get the date of the orbit.
  138.      * @return the date
  139.      */
  140.     public FieldAbsoluteDate<T> getDate() {
  141.         return date;
  142.     }

  143.     /** Get the definition frame of the orbit.
  144.      * @return the definition frame
  145.      */
  146.     public Frame getFrame() {
  147.         return frame;
  148.     }

  149.     /** Get the eccentricity.
  150.      * @return ecc
  151.      */
  152.     public T getEcc() {
  153.         return ecc;
  154.     }

  155.     /** Get the Keplerian mean motion.
  156.      * @return n
  157.      */
  158.     public T getMeanMotion() {
  159.         return n;
  160.     }

  161.     /** Get the Keplerian period.
  162.      * @return period
  163.      */
  164.     public T getKeplerianPeriod() {
  165.         return period;
  166.     }

  167.     /** Get the semi-major axis.
  168.      * @return the semi-major axis a
  169.      */
  170.     public T getSma() {
  171.         return sma;
  172.     }

  173.     /** Get the x component of eccentricity vector.
  174.      * <p>
  175.      * This element called k in DSST corresponds to ex
  176.      * for the {@link org.orekit.orbits.EquinoctialOrbit}
  177.      * </p>
  178.      * @return k
  179.      */
  180.     public T getK() {
  181.         return k;
  182.     }

  183.     /** Get the y component of eccentricity vector.
  184.      * <p>
  185.      * This element called h in DSST corresponds to ey
  186.      * for the {@link org.orekit.orbits.EquinoctialOrbit}
  187.      * </p>
  188.      * @return h
  189.      */
  190.     public T getH() {
  191.         return h;
  192.     }

  193.     /** Get the x component of inclination vector.
  194.      * <p>
  195.      * This element called q in DSST corresponds to hx
  196.      * for the {@link org.orekit.orbits.EquinoctialOrbit}
  197.      * </p>
  198.      * @return q
  199.      */
  200.     public T getQ() {
  201.         return q;
  202.     }

  203.     /** Get the y component of inclination vector.
  204.      * <p>
  205.      * This element called p in DSST corresponds to hy
  206.      * for the {@link org.orekit.orbits.EquinoctialOrbit}
  207.      * </p>
  208.      * @return p
  209.      */
  210.     public T getP() {
  211.         return p;
  212.     }

  213.     /** Get the mean longitude.
  214.      * @return lm
  215.      */
  216.     public T getLM() {
  217.         return lm;
  218.     }

  219.     /** Get the true longitude.
  220.      * @return lv
  221.      */
  222.     public T getLv() {
  223.         return lv;
  224.     }

  225.     /** Get the eccentric longitude.
  226.      * @return le
  227.      */
  228.     public T getLe() {
  229.         return le;
  230.     }

  231.     /** Get the retrograde factor.
  232.      * @return the retrograde factor I
  233.      */
  234.     public int getRetrogradeFactor() {
  235.         return I;
  236.     }

  237.     /** Get B = sqrt(1 - e²).
  238.      * @return B
  239.      */
  240.     public T getB() {
  241.         return B;
  242.     }

  243.     /** Get C = 1 + p² + q².
  244.      * @return C
  245.      */
  246.     public T getC() {
  247.         return C;
  248.     }

  249.     /** Get equinoctial frame vector f.
  250.      * @return f vector
  251.      */
  252.     public FieldVector3D<T> getVectorF() {
  253.         return f;
  254.     }

  255.     /** Get equinoctial frame vector g.
  256.      * @return g vector
  257.      */
  258.     public FieldVector3D<T> getVectorG() {
  259.         return g;
  260.     }

  261.     /** Get equinoctial frame vector w.
  262.      * @return w vector
  263.      */
  264.     public FieldVector3D<T> getVectorW() {
  265.         return w;
  266.     }

  267.     /** Transforms the FieldAuxiliaryElements instance into an AuxiliaryElements instance.
  268.      * @return the AuxiliaryElements instance
  269.      * @since 11.3.3
  270.      */
  271.     public AuxiliaryElements toAuxiliaryElements() {
  272.         return new AuxiliaryElements(orbit.toOrbit(), getRetrogradeFactor());
  273.     }

  274.     /** Get direction cosine α for central body.
  275.      * @return α
  276.      * @deprecated since 12.2, use {@link FieldDSSTGravityContext#getAlpha()} instead
  277.      */
  278.     @Deprecated
  279.     public T getAlpha() {
  280.         return (T) f.getZ();
  281.     }

  282.     /** Get direction cosine β for central body.
  283.      * @return β
  284.      * @deprecated since 12.2, use {@link FieldDSSTGravityContext#getBeta()} instead
  285.      */
  286.     @Deprecated
  287.     public T getBeta() {
  288.         return (T) g.getZ();
  289.     }

  290.     /** Get direction cosine γ for central body.
  291.      * @return γ
  292.      * @deprecated since 12.2, use {@link FieldDSSTGravityContext#getGamma()} instead
  293.      */
  294.     @Deprecated
  295.     public T getGamma() {
  296.         return (T) w.getZ();
  297.     }
  298. }