CoefficientsFactory.java

  1. /* Copyright 2002-2022 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 java.util.SortedMap;
  20. import java.util.TreeMap;
  21. import java.util.concurrent.ConcurrentSkipListMap;

  23. import org.hipparchus.Field;
  24. import org.hipparchus.CalculusFieldElement;
  25. import org.hipparchus.util.CombinatoricsUtils;
  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.MathArrays;
  28. import org.orekit.errors.OrekitException;
  29. import org.orekit.errors.OrekitMessages;

  31. /**
  32.  * This class is designed to provide coefficient from the DSST theory.
  33.  *
  34.  * @author Romain Di Costanzo
  35.  */
  36. public class CoefficientsFactory {

  38.     /** Internal storage of the polynomial values. Reused for further computation. */
  39.     private static SortedMap<NSKey, Double> VNS = new ConcurrentSkipListMap<NSKey, Double>();

  41.     /** Last computed order for V<sub>ns</sub> coefficients. */
  42.     private static int         LAST_VNS_ORDER = 2;

  44.     /** Static initialization for the V<sub>ns</sub> coefficient. */
  45.     static {
  46.         // Initialization
  47.         VNS.put(new NSKey(0, 0), 1.);
  48.         VNS.put(new NSKey(1, 0), 0.);
  49.         VNS.put(new NSKey(1, 1), 0.5);
  50.     }

  52.     /** Private constructor as the class is a utility class.
  53.      */
  54.     private CoefficientsFactory() {
  55.     }

  57.     /** Compute the Q<sub>n,s</sub> coefficients evaluated at γ from the recurrence formula 2.8.3-(2).
  58.      *  <p>
  59.      *  Q<sub>n,s</sub> coefficients are computed for n = 0 to nMax
  60.      *  and s = 0 to sMax + 1 in order to also get the derivative dQ<sub>n,s</sub>/dγ = Q(n, s + 1)
  61.      *  </p>
  62.      *  @param gamma γ angle
  63.      *  @param nMax n max value
  64.      *  @param sMax s max value
  65.      *  @return Q<sub>n,s</sub> coefficients array
  66.      */
  67.     public static double[][] computeQns(final double gamma, final int nMax, final int sMax) {

  69.         // Initialization
  70.         final int sDim = FastMath.min(sMax + 1, nMax) + 1;
  71.         final int rows = nMax + 1;
  72.         final double[][] Qns = new double[rows][];
  73.         for (int i = 0; i <= nMax; i++) {
  74.             final int snDim = FastMath.min(i + 1, sDim);
  75.             Qns[i] = new double[snDim];
  76.         }

  78.         // first element
  79.         Qns[0][0] = 1;

  81.         for (int n = 1; n <= nMax; n++) {
  82.             final int snDim = FastMath.min(n + 1, sDim);
  83.             for (int s = 0; s < snDim; s++) {
  84.                 if (n == s) {
  85.                     Qns[n][s] = (2. * s - 1.) * Qns[s - 1][s - 1];
  86.                 } else if (n == (s + 1)) {
  87.                     Qns[n][s] = (2. * s + 1.) * gamma * Qns[s][s];
  88.                 } else {
  89.                     Qns[n][s] = (2. * n - 1.) * gamma * Qns[n - 1][s] - (n + s - 1.) * Qns[n - 2][s];
  90.                     Qns[n][s] /= n - s;
  91.                 }
  92.             }
  93.         }

  95.         return Qns;
  96.     }

  98.     /** Compute the Q<sub>n,s</sub> coefficients evaluated at γ from the recurrence formula 2.8.3-(2).
  99.      *  <p>
  100.      *  Q<sub>n,s</sub> coefficients are computed for n = 0 to nMax
  101.      *  and s = 0 to sMax + 1 in order to also get the derivative dQ<sub>n,s</sub>/dγ = Q(n, s + 1)
  102.      *  </p>
  103.      *  @param gamma γ angle
  104.      *  @param nMax n max value
  105.      *  @param sMax s max value
  106.      *  @param <T> the type of the field elements
  107.      *  @return Q<sub>n,s</sub> coefficients array
  108.      */
  109.     public static <T extends CalculusFieldElement<T>> T[][] computeQns(final T gamma, final int nMax, final int sMax) {

  111.         // Initialization
  112.         final int sDim = FastMath.min(sMax + 1, nMax) + 1;
  113.         final int rows = nMax + 1;
  114.         final T[][] Qns = MathArrays.buildArray(gamma.getField(), rows, FastMath.min(nMax + 1, sDim) - 1);
  115.         for (int i = 0; i <= nMax; i++) {
  116.             final int snDim = FastMath.min(i + 1, sDim);
  117.             Qns[i] = MathArrays.buildArray(gamma.getField(), snDim);
  118.         }

  120.         // first element
  121.         Qns[0][0] = gamma.subtract(gamma).add(1.);

  123.         for (int n = 1; n <= nMax; n++) {
  124.             final int snDim = FastMath.min(n + 1, sDim);
  125.             for (int s = 0; s < snDim; s++) {
  126.                 if (n == s) {
  127.                     Qns[n][s] = Qns[s - 1][s - 1].multiply(2. * s - 1.);
  128.                 } else if (n == (s + 1)) {
  129.                     Qns[n][s] = Qns[s][s].multiply(gamma).multiply(2. * s + 1.);
  130.                 } else {
  131.                     Qns[n][s] = Qns[n - 1][s].multiply(gamma).multiply(2. * n - 1.).subtract(Qns[n - 2][s].multiply(n + s - 1.));
  132.                     Qns[n][s] = Qns[n][s].divide(n - s);
  133.                 }
  134.             }
  135.         }

  137.         return Qns;
  138.     }

  140.     /** Compute recursively G<sub>s</sub> and H<sub>s</sub> polynomials from equation 3.1-(5).
  141.      *  @param k x-component of the eccentricity vector
  142.      *  @param h y-component of the eccentricity vector
  143.      *  @param alpha 1st direction cosine
  144.      *  @param beta 2nd direction cosine
  145.      *  @param order development order
  146.      *  @return Array of G<sub>s</sub> and H<sub>s</sub> polynomials for s from 0 to order.<br>
  147.      *          The 1st column contains the G<sub>s</sub> values.
  148.      *          The 2nd column contains the H<sub>s</sub> values.
  149.      */
  150.     public static double[][] computeGsHs(final double k, final double h,
  151.                                          final double alpha, final double beta,
  152.                                          final int order) {
  153.         // Constant terms
  154.         final double hamkb = h * alpha - k * beta;
  155.         final double kaphb = k * alpha + h * beta;
  156.         // Initialization
  157.         final double[][] GsHs = new double[2][order + 1];
  158.         GsHs[0][0] = 1.;
  159.         GsHs[1][0] = 0.;

  161.         for (int s = 1; s <= order; s++) {
  162.             // Gs coefficient
  163.             GsHs[0][s] = kaphb * GsHs[0][s - 1] - hamkb * GsHs[1][s - 1];
  164.             // Hs coefficient
  165.             GsHs[1][s] = hamkb * GsHs[0][s - 1] + kaphb * GsHs[1][s - 1];
  166.         }

  168.         return GsHs;
  169.     }

  171.     /** Compute recursively G<sub>s</sub> and H<sub>s</sub> polynomials from equation 3.1-(5).
  172.      *  @param k x-component of the eccentricity vector
  173.      *  @param h y-component of the eccentricity vector
  174.      *  @param alpha 1st direction cosine
  175.      *  @param beta 2nd direction cosine
  176.      *  @param order development order
  177.      *  @param field field of elements
  178.      *  @param <T> the type of the field elements
  179.      *  @return Array of G<sub>s</sub> and H<sub>s</sub> polynomials for s from 0 to order.<br>
  180.      *          The 1st column contains the G<sub>s</sub> values.
  181.      *          The 2nd column contains the H<sub>s</sub> values.
  182.      */
  183.     public static <T extends CalculusFieldElement<T>> T[][] computeGsHs(final T k, final T h,
  184.                                          final T alpha, final T beta,
  185.                                          final int order, final Field<T> field) {
  186.         // Zero for initialization
  187.         final T zero = field.getZero();

  189.         // Constant terms
  190.         final T hamkb = h.multiply(alpha).subtract(k.multiply(beta));
  191.         final T kaphb = k.multiply(alpha).add(h.multiply(beta));
  192.         // Initialization
  193.         final T[][] GsHs = MathArrays.buildArray(field, 2, order + 1);
  194.         GsHs[0][0] = zero.add(1.);
  195.         GsHs[1][0] = zero;

  197.         for (int s = 1; s <= order; s++) {
  198.             // Gs coefficient
  199.             GsHs[0][s] = kaphb.multiply(GsHs[0][s - 1]).subtract(hamkb.multiply(GsHs[1][s - 1]));
  200.             // Hs coefficient
  201.             GsHs[1][s] = hamkb.multiply(GsHs[0][s - 1]).add(kaphb.multiply(GsHs[1][s - 1]));
  202.         }

  204.         return GsHs;
  205.     }

  207.     /** Compute the V<sub>n,s</sub> coefficients from 2.8.2-(1)(2).
  208.      * @param order Order of the computation. Computation will be done from 0 to order -1
  209.      * @return Map of the V<sub>n, s</sub> coefficients
  210.      * @deprecated as of 11.3.3, replaced by {@link #computeVnsCoefficients(int)}
  211.      */
  212.     @Deprecated
  213.     public static TreeMap<NSKey, Double> computeVns(final int order) {
  214.         return new TreeMap<>(computeVnsCoefficients(order));
  215.     }

  217.     /** Compute the V<sub>n,s</sub> coefficients from 2.8.2-(1)(2).
  218.      * @param order Order of the computation. Computation will be done from 0 to order -1
  219.      * @return Map of the V<sub>n, s</sub> coefficients
  220.      * @since 11.3.3
  221.      */
  222.     // TODO rename to computeVns in 12.0?
  223.     public static SortedMap<NSKey, Double> computeVnsCoefficients(final int order) {

  225.         if (order > LAST_VNS_ORDER) {
  226.             // Compute coefficient
  227.             // Need previous computation as recurrence relation is done at s + 1 and n + 2
  228.             final int min = (LAST_VNS_ORDER - 2 < 0) ? 0 : (LAST_VNS_ORDER - 2);
  229.             for (int n = min; n < order; n++) {
  230.                 for (int s = 0; s < n + 1; s++) {
  231.                     if ((n - s) % 2 != 0) {
  232.                         VNS.put(new NSKey(n, s), 0.);
  233.                     } else {
  234.                         // s = n
  235.                         if (n == s && (s + 1) < order) {
  236.                             VNS.put(new NSKey(s + 1, s + 1), VNS.get(new NSKey(s, s)) / (2 * s + 2.));
  237.                         }
  238.                         // otherwise
  239.                         if ((n + 2) < order) {
  240.                             VNS.put(new NSKey(n + 2, s), VNS.get(new NSKey(n, s)) * (-n + s - 1.) / (n + s + 2.));
  241.                         }
  242.                     }
  243.                 }
  244.             }
  245.             LAST_VNS_ORDER = order;
  246.         }
  247.         return new ConcurrentSkipListMap<>(VNS);
  248.     }

  250.     /** Get the V<sub>n,s</sub><sup>m</sup> coefficient from V<sub>n,s</sub>.
  251.      *  <br>See § 2.8.2 in Danielson paper.
  252.      * @param m m
  253.      * @param n n
  254.      * @param s s
  255.      * @return The V<sub>n, s</sub> <sup>m</sup> coefficient
  256.      */
  257.     public static double getVmns(final int m, final int n, final int s) {
  258.         if (m > n) {
  259.             throw new OrekitException(OrekitMessages.DSST_VMNS_COEFFICIENT_ERROR_MS, m, n);
  260.         }
  261.         final double fns = CombinatoricsUtils.factorialDouble(n + FastMath.abs(s));
  262.         final double fnm = CombinatoricsUtils.factorialDouble(n  - m);

  264.         double result = 0;
  265.         // If (n - s) is odd, the Vmsn coefficient is null
  266.         if ((n - s) % 2 == 0) {
  267.             // Update the Vns coefficient
  268.             if ((n + 1) > LAST_VNS_ORDER) {
  269.                 computeVnsCoefficients(n + 1);
  270.             }
  271.             if (s >= 0) {
  272.                 result = fns  * VNS.get(new NSKey(n, s)) / fnm;
  273.             } else {
  274.                 // If s < 0 : Vmn-s = (-1)^(-s) Vmns
  275.                 final int mops = (s % 2 == 0) ? 1 : -1;
  276.                 result = mops * fns * VNS.get(new NSKey(n, -s)) / fnm;
  277.             }
  278.         }
  279.         return result;
  280.     }

  282.     /** Key formed by two integer values. */
  283.     public static class NSKey implements Comparable<NSKey> {

  285.         /** n value. */
  286.         private final int n;

  288.         /** s value. */
  289.         private final int s;

  291.         /** Simple constructor.
  292.          * @param n n
  293.          * @param s s
  294.          */
  295.         public NSKey(final int n, final int s) {
  296.             this.n = n;
  297.             this.s = s;
  298.         }

  300.         /** Get n.
  301.          * @return n
  302.          */
  303.         public int getN() {
  304.             return n;
  305.         }

  307.         /** Get s.
  308.          * @return s
  309.          */
  310.         public int getS() {
  311.             return s;
  312.         }

  314.         /** {@inheritDoc} */
  315.         public int compareTo(final NSKey key) {
  316.             int result = 1;
  317.             if (n == key.n) {
  318.                 if (s < key.s) {
  319.                     result = -1;
  320.                 } else if (s == key.s) {
  321.                     result = 0;
  322.                 }
  323.             } else if (n < key.n) {
  324.                 result = -1;
  325.             }
  326.             return result;
  327.         }

  329.         /** {@inheritDoc} */
  330.         public boolean equals(final Object key) {

  332.             if (key == this) {
  333.                 // first fast check
  334.                 return true;
  335.             }

  337.             if (key instanceof NSKey) {
  338.                 return n == ((NSKey) key).n && s == ((NSKey) key).s;
  339.             }

  341.             return false;

  343.         }

  345.         /** {@inheritDoc} */
  346.         public int hashCode() {
  347.             return 0x998493a6 ^ (n << 8) ^ s;
  348.         }

  350.     }

  352. }
Created with JaCoCo 0.8.8.202204050719