FieldGammaMnsFunction.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.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.  * @param <T> type of the field elements
  26.  */
  27. public class FieldGammaMnsFunction <T extends CalculusFieldElement<T>> {

  28.     /** Factorial ratios. */
  29.     private static double[] PRECOMPUTED_RATIOS;

  30.     /** Field element. */
  31.     private final Field<T> field;

  32.     /** Factorial ratios. */
  33.     private final double[] ratios;

  34.     /** Storage array. */
  35.     private final T[] values;

  36.     /** 1 + I * γ. */
  37.     private final T opIg;

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

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

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

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

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

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

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

  90.                 }

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

  96.             }
  97.             return PRECOMPUTED_RATIOS;
  98.         }
  99.     }

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

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

  137. }