1   /* Copyright 2002-2025 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  
19  import org.hipparchus.fraction.BigFraction;
20  import org.hipparchus.util.FastMath;
21  
22  import java.util.Arrays;
23  
24  /** Compute the &Gamma;<sup>m</sup><sub>n,s</sub>(γ) function from equation 2.7.1-(13).
25   *
26   *  @author Romain Di Costanzo
27   */
28  public class GammaMnsFunction {
29  
30      /** Factorial ratios. */
31      private static double[] PRECOMPUTED_RATIOS;
32  
33      /** Factorial ratios. */
34      private final double[] ratios;
35  
36      /** Storage array. */
37      private final double[] values;
38  
39      /** 1 + I * γ. */
40      private final double opIg;
41  
42      /** I = +1 for a prograde orbit, -1 otherwise. */
43      private final int    I;
44  
45      /** Simple constructor.
46       *  @param nMax max value for n
47       *  @param gamma γ
48       *  @param I retrograde factor
49       */
50      public GammaMnsFunction(final int nMax, final double gamma, final int I) {
51          final int size = (nMax + 1) * (nMax + 2) * (4 * nMax + 3) / 6;
52          this.values = new double[size];
53          this.ratios = getRatios(nMax, size);
54          Arrays.fill(values, Double.NaN);
55          this.opIg   = 1. + I * gamma;
56          this.I      = I;
57      }
58  
59      /** Compute the array index.
60       *  @param m m
61       *  @param n n
62       *  @param s s
63       *  @return index for element m, n, s
64       */
65      private static int index(final int m, final int n, final int s) {
66          return n * (n + 1) * (4 * n - 1) / 6 + // index for 0, n, 0
67                 m * (2 * n + 1) +               // index for m, n, 0
68                 s + n;                          // index for m, n, s
69      }
70  
71      /** Get the ratios for the given size.
72       * @param nMax max value for n
73       * @param size ratio size array
74       * @return factorial ratios
75       */
76      private static double[] getRatios(final int nMax, final int size) {
77          synchronized (GammaMnsFunction.class) {
78              if (PRECOMPUTED_RATIOS == null || PRECOMPUTED_RATIOS.length < size) {
79                  // we need to compute a larger reference array
80  
81                  final BigFraction[] bF = new BigFraction[size];
82                  for (int n = 0; n <= nMax; ++n) {
83  
84                      // populate ratios for s = 0
85                      bF[index(0, n, 0)] = BigFraction.ONE;
86                      for (int m = 1; m <= n; ++m) {
87                          bF[index(m, n, 0)] = bF[index(m - 1, n, 0)].multiply(n + m).divide(n - (m - 1));
88                      }
89  
90                      // populate ratios for s != 0
91                      for (int absS = 1; absS <= n; ++absS) {
92                          for (int m = 0; m <= n; ++m) {
93                              bF[index(m, n, +absS)] = bF[index(m, n, absS - 1)].divide(n + absS).multiply(n - (absS - 1));
94                              bF[index(m, n, -absS)] = bF[index(m, n, absS)];
95                          }
96                      }
97  
98                  }
99  
100                 // convert to double
101                 PRECOMPUTED_RATIOS = new double[size];
102                 for (int i = 0; i < bF.length; ++i) {
103                     PRECOMPUTED_RATIOS[i] = bF[i].doubleValue();
104                 }
105 
106             }
107             return PRECOMPUTED_RATIOS;
108         }
109     }
110 
111     /** Get &Gamma; function value.
112      *  @param m m
113      *  @param n n
114      *  @param s s
115      *  @return &Gamma;<sup>m</sup><sub>n, s</sub>(γ)
116      */
117     public double getValue(final int m, final int n, final int s) {
118         final int i = index(m, n, s);
119         if (Double.isNaN(values[i])) {
120             if (s <= -m) {
121                 values[i] = (((m - s) & 0x1) == 0 ? +1 : -1) * FastMath.scalb(FastMath.pow(opIg, -I * m), s);
122             } else if (s <= m) {
123                 values[i] = (((m - s) & 0x1) == 0 ? +1 : -1) * FastMath.scalb(FastMath.pow(opIg, I * s), -m) * ratios[i];
124             } else {
125                 values[i] = FastMath.scalb(FastMath.pow(opIg, I * m), -s);
126             }
127         }
128         return values[i];
129     }
130 
131     /** Get &Gamma; function derivative.
132      * @param m m
133      * @param n n
134      * @param s s
135      * @return d&Gamma;<sup>m</sup><sub>n,s</sub>(γ)/dγ
136      */
137     public double getDerivative(final int m, final int n, final int s) {
138         if (s <= -m) {
139             return -m * I * getValue(m, n, s) / opIg;
140         } else if (s >= m) {
141             return m * I * getValue(m, n, s) / opIg;
142         } else {
143             return s * I * getValue(m, n, s) / opIg;
144         }
145     }
146 
147 }