GammaMnsFunction.java
/* Copyright 2002-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation.semianalytical.dsst.utilities;
import java.util.Arrays;
import org.hipparchus.fraction.BigFraction;
import org.hipparchus.util.FastMath;
/** Compute the Γ<sup>m</sup><sub>n,s</sub>(γ) function from equation 2.7.1-(13).
*
* @author Romain Di Costanzo
*/
public class GammaMnsFunction {
/** Factorial ratios. */
private static double[] PRECOMPUTED_RATIOS;
/** Factorial ratios. */
private final double[] ratios;
/** Storage array. */
private final double[] values;
/** 1 + I * γ. */
private final double opIg;
/** I = +1 for a prograde orbit, -1 otherwise. */
private final int I;
/** Simple constructor.
* @param nMax max value for n
* @param gamma γ
* @param I retrograde factor
*/
public GammaMnsFunction(final int nMax, final double gamma, final int I) {
final int size = (nMax + 1) * (nMax + 2) * (4 * nMax + 3) / 6;
this.values = new double[size];
this.ratios = getRatios(nMax, size);
Arrays.fill(values, Double.NaN);
this.opIg = 1. + I * gamma;
this.I = I;
}
/** Compute the array index.
* @param m m
* @param n n
* @param s s
* @return index for element m, n, s
*/
private static int index(final int m, final int n, final int s) {
return n * (n + 1) * (4 * n - 1) / 6 + // index for 0, n, 0
m * (2 * n + 1) + // index for m, n, 0
s + n; // index for m, n, s
}
/** Get the ratios for the given size.
* @param nMax max value for n
* @param size ratio size array
* @return factorial ratios
*/
private static double[] getRatios(final int nMax, final int size) {
synchronized (GammaMnsFunction.class) {
if (PRECOMPUTED_RATIOS == null || PRECOMPUTED_RATIOS.length < size) {
// we need to compute a larger reference array
final BigFraction[] bF = new BigFraction[size];
for (int n = 0; n <= nMax; ++n) {
// populate ratios for s = 0
bF[index(0, n, 0)] = BigFraction.ONE;
for (int m = 1; m <= n; ++m) {
bF[index(m, n, 0)] = bF[index(m - 1, n, 0)].multiply(n + m).divide(n - (m - 1));
}
// populate ratios for s != 0
for (int absS = 1; absS <= n; ++absS) {
for (int m = 0; m <= n; ++m) {
bF[index(m, n, +absS)] = bF[index(m, n, absS - 1)].divide(n + absS).multiply(n - (absS - 1));
bF[index(m, n, -absS)] = bF[index(m, n, absS)];
}
}
}
// convert to double
PRECOMPUTED_RATIOS = new double[size];
for (int i = 0; i < bF.length; ++i) {
PRECOMPUTED_RATIOS[i] = bF[i].doubleValue();
}
}
return PRECOMPUTED_RATIOS;
}
}
/** Get Γ function value.
* @param m m
* @param n n
* @param s s
* @return Γ<sup>m</sup><sub>n, s</sub>(γ)
*/
public double getValue(final int m, final int n, final int s) {
final int i = index(m, n, s);
if (Double.isNaN(values[i])) {
if (s <= -m) {
values[i] = (((m - s) & 0x1) == 0 ? +1 : -1) * FastMath.scalb(FastMath.pow(opIg, -I * m), s);
} else if (s <= m) {
values[i] = (((m - s) & 0x1) == 0 ? +1 : -1) * FastMath.scalb(FastMath.pow(opIg, I * s), -m) * ratios[i];
} else {
values[i] = FastMath.scalb(FastMath.pow(opIg, I * m), -s);
}
}
return values[i];
}
/** Get Γ function derivative.
* @param m m
* @param n n
* @param s s
* @return dΓ<sup>m</sup><sub>n,s</sub>(γ)/dγ
*/
public double getDerivative(final int m, final int n, final int s) {
double res = 0.;
if (s <= -m) {
res = -m * I * getValue(m, n, s) / opIg;
} else if (s >= m) {
res = m * I * getValue(m, n, s) / opIg;
} else {
res = s * I * getValue(m, n, s) / opIg;
}
return res;
}
}