FieldGammaMnsFunction.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;

  18. import java.util.Arrays;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.CalculusFieldElement;
  21. import org.hipparchus.fraction.BigFraction;
  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.MathArrays;

  24. /** Compute the &Gamma;<sup>m</sup><sub>n,s</sub>(γ) function from equation 2.7.1-(13). */
  25. public class FieldGammaMnsFunction <T extends CalculusFieldElement<T>> {

  26.     /** Factorial ratios. */
  27.     private static double[] PRECOMPUTED_RATIOS;

  28.     /** Field element. */
  29.     private final Field<T> field;

  30.     /** Factorial ratios. */
  31.     private final double[] ratios;

  32.     /** Storage array. */
  33.     private final T[] values;

  34.     /** 1 + I * γ. */
  35.     private final T opIg;

  36.     /** I = +1 for a prograde orbit, -1 otherwise. */
  37.     private final int    I;

  38.     /** Simple constructor.
  39.      *  @param nMax max value for n
  40.      *  @param gamma γ
  41.      *  @param I retrograde factor
  42.      *  @param field field element
  43.      */
  44.     public FieldGammaMnsFunction(final int nMax, final T gamma, final int I, final Field<T> field) {
  45.         this.field = field;
  46.         final T zero = field.getZero();
  47.         final int size = (nMax + 1) * (nMax + 2) * (4 * nMax + 3) / 6;
  48.         this.values = MathArrays.buildArray(field, size);
  49.         this.ratios = getRatios(nMax, size);
  50.         Arrays.fill(values, zero.add(Double.NaN));
  51.         this.opIg   = gamma.multiply(I).add(1.);
  52.         this.I      = I;
  53.     }

  54.     /** Compute the array index.
  55.      *  @param m m
  56.      *  @param n n
  57.      *  @param s s
  58.      *  @return index for element m, n, s
  59.      */
  60.     private static int index(final int m, final int n, final int s) {
  61.         return n * (n + 1) * (4 * n - 1) / 6 + // index for 0, n, 0
  62.                m * (2 * n + 1) +               // index for m, n, 0
  63.                s + n;                          // index for m, n, s
  64.     }

  65.     /** Get the ratios for the given size.
  66.      * @param nMax max value for n
  67.      * @param size ratio size array
  68.      * @return factorial ratios
  69.      */
  70.     private static double[] getRatios(final int nMax, final int size) {
  71.         synchronized (GammaMnsFunction.class) {
  72.             if (PRECOMPUTED_RATIOS == null || PRECOMPUTED_RATIOS.length < size) {
  73.                 // we need to compute a larger reference array

  74.                 final BigFraction[] bF = new BigFraction[size];
  75.                 for (int n = 0; n <= nMax; ++n) {

  76.                     // populate ratios for s = 0
  77.                     bF[index(0, n, 0)] = BigFraction.ONE;
  78.                     for (int m = 1; m <= n; ++m) {
  79.                         bF[index(m, n, 0)] = bF[index(m - 1, n, 0)].multiply(n + m).divide(n - (m - 1));
  80.                     }

  81.                     // populate ratios for s != 0
  82.                     for (int absS = 1; absS <= n; ++absS) {
  83.                         for (int m = 0; m <= n; ++m) {
  84.                             bF[index(m, n, +absS)] = bF[index(m, n, absS - 1)].divide(n + absS).multiply(n - (absS - 1));
  85.                             bF[index(m, n, -absS)] = bF[index(m, n, absS)];
  86.                         }
  87.                     }

  88.                 }

  89.                 // convert to double
  90.                 PRECOMPUTED_RATIOS = new double[size];
  91.                 for (int i = 0; i < bF.length; ++i) {
  92.                     PRECOMPUTED_RATIOS[i] = bF[i].doubleValue();
  93.                 }

  94.             }
  95.             return PRECOMPUTED_RATIOS;
  96.         }
  97.     }

  98.     /** Get &Gamma; function value.
  99.      *  @param m m
  100.      *  @param n n
  101.      *  @param s s
  102.      *  @return &Gamma;<sup>m</sup><sub>n, s</sub>(γ)
  103.      */
  104.     public T getValue(final int m, final int n, final int s) {
  105.         final int i = index(m, n, s);
  106.         if (Double.isNaN(values[i].getReal())) {
  107.             if (s <= -m) {
  108.                 values[i] = FastMath.scalb(FastMath.pow(opIg, -I * m), s).multiply(((m - s) & 0x1) == 0 ? +1 : -1);
  109.             } else if (s <= m) {
  110.                 values[i] = FastMath.scalb(FastMath.pow(opIg, I * s), -m).multiply(ratios[i]).multiply(((m - s) & 0x1) == 0 ? +1 : -1);
  111.             } else {
  112.                 values[i] = FastMath.scalb(FastMath.pow(opIg, I * m), -s);
  113.             }
  114.         }
  115.         return values[i];
  116.     }

  117.     /** Get &Gamma; function derivative.
  118.      * @param m m
  119.      * @param n n
  120.      * @param s s
  121.      * @return d&Gamma;<sup>m</sup><sub>n,s</sub>(γ)/dγ
  122.      */
  123.     public T getDerivative(final int m, final int n, final int s) {
  124.         final T zero = field.getZero();
  125.         T res = zero;
  126.         if (s <= -m) {
  127.             res = getValue(m, n, s).multiply(I).multiply(-m).divide(opIg);
  128.         } else if (s >= m) {
  129.             res =  getValue(m, n, s).multiply(I).multiply(m).divide(opIg);;
  130.         } else {
  131.             res =  getValue(m, n, s).multiply(I).multiply(s).divide(opIg);;
  132.         }
  133.         return res;
  134.     }

  135. }