DSSTThirdBody.java
/* Copyright 2002-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (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.forces;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.TreeMap;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.util.FastMath;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.bodies.CelestialBody;
import org.orekit.errors.OrekitException;
import org.orekit.orbits.Orbit;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
import org.orekit.propagation.semianalytical.dsst.utilities.JacobiPolynomials;
import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenThirdBodyLinear;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.TimeSpanMap;
/** Third body attraction perturbation to the
* {@link org.orekit.propagation.semianalytical.dsst.DSSTPropagator DSSTPropagator}.
*
* @author Romain Di Costanzo
* @author Pascal Parraud
*/
public class DSSTThirdBody implements DSSTForceModel {
/** Max power for summation. */
private static final int MAX_POWER = 22;
/** Truncation tolerance for big, eccentric orbits. */
private static final double BIG_TRUNCATION_TOLERANCE = 1.e-1;
/** Truncation tolerance for small orbits. */
private static final double SMALL_TRUNCATION_TOLERANCE = 1.9e-6;
/** Number of points for interpolation. */
private static final int INTERPOLATION_POINTS = 3;
/** Maximum power for eccentricity used in short periodic computation. */
private static final int MAX_ECCPOWER_SP = 4;
/** Retrograde factor I.
* <p>
* DSST model needs equinoctial orbit as internal representation.
* Classical equinoctial elements have discontinuities when inclination
* is close to zero. In this representation, I = +1. <br>
* To avoid this discontinuity, another representation exists and equinoctial
* elements can be expressed in a different way, called "retrograde" orbit.
* This implies I = -1. <br>
* As Orekit doesn't implement the retrograde orbit, I is always set to +1.
* But for the sake of consistency with the theory, the retrograde factor
* has been kept in the formulas.
* </p>
*/
private static final int I = 1;
/** The 3rd body to consider. */
private final CelestialBody body;
/** Standard gravitational parameter μ for the body in m³/s². */
private final double gm;
/** Factorial. */
private final double[] fact;
/** V<sub>ns</sub> coefficients. */
private final TreeMap<NSKey, Double> Vns;
/** Distance from center of mass of the central body to the 3rd body. */
private double R3;
/** Short period terms. */
private ThirdBodyShortPeriodicCoefficients shortPeriods;
// Equinoctial elements (according to DSST notation)
/** a. */
private double a;
/** e<sub>x</sub>. */
private double k;
/** e<sub>y</sub>. */
private double h;
/** h<sub>x</sub>. */
private double q;
/** h<sub>y</sub>. */
private double p;
/** Eccentricity. */
private double ecc;
// Direction cosines of the symmetry axis
/** α. */
private double alpha;
/** β. */
private double beta;
/** γ. */
private double gamma;
// Common factors for potential computation
/** A = n * a². */
private double A;
/** B = sqrt(1 - e²). */
private double B;
/** C = 1 + p² + q². */
private double C;
/** B². */
private double BB;
/** B³. */
private double BBB;
/** The mean motion (n). */
private double meanMotion;
/** Χ = 1 / sqrt(1 - e²) = 1 / B. */
private double X;
/** Χ². */
private double XX;
/** Χ³. */
private double XXX;
/** -2 * a / A. */
private double m2aoA;
/** B / A. */
private double BoA;
/** 1 / (A * B). */
private double ooAB;
/** -C / (2 * A * B). */
private double mCo2AB;
/** B / A(1 + B). */
private double BoABpo;
/** Max power for a/R3 in the serie expansion. */
private int maxAR3Pow;
/** Max power for e in the serie expansion. */
private int maxEccPow;
/** Max power for e in the serie expansion (for short periodics). */
private int maxEccPowShort;
/** Max frequency of F. */
private int maxFreqF;
/** An array that contains the objects needed to build the Hansen coefficients. <br/>
* The index is s */
private final HansenThirdBodyLinear[] hansenObjects;
/** The current value of the U function. <br/>
* Needed for the short periodic contribution */
private double U;
/** Qns coefficients. */
private double[][] Qns;
/** a / R3 up to power maxAR3Pow. */
private double[] aoR3Pow;
/** mu3 / R3. */
private double muoR3;
/** b = 1 / (1 + sqrt(1 - e²)) = 1 / (1 + B).*/
private double b;
/** h * Χ³. */
private double hXXX;
/** k * Χ³. */
private double kXXX;
/** Complete constructor.
* @param body the 3rd body to consider
* @see org.orekit.bodies.CelestialBodyFactory
*/
public DSSTThirdBody(final CelestialBody body) {
this.body = body;
this.gm = body.getGM();
this.maxAR3Pow = Integer.MIN_VALUE;
this.maxEccPow = Integer.MIN_VALUE;
this.Vns = CoefficientsFactory.computeVns(MAX_POWER);
// Factorials computation
final int dim = 2 * MAX_POWER;
this.fact = new double[dim];
fact[0] = 1.;
for (int i = 1; i < dim; i++) {
fact[i] = i * fact[i - 1];
}
//Initialise the HansenCoefficient generator
this.hansenObjects = new HansenThirdBodyLinear[MAX_POWER + 1];
for (int s = 0; s <= MAX_POWER; s++) {
this.hansenObjects[s] = new HansenThirdBodyLinear(MAX_POWER, s);
}
}
/** Get third body.
* @return third body
*/
public CelestialBody getBody() {
return body;
}
/** Computes the highest power of the eccentricity and the highest power
* of a/R3 to appear in the truncated analytical power series expansion.
* <p>
* This method computes the upper value for the 3rd body potential and
* determines the maximal powers for the eccentricity and a/R3 producing
* potential terms bigger than a defined tolerance.
* </p>
* @param aux auxiliary elements related to the current orbit
* @param meanOnly only mean elements will be used for the propagation
* @throws OrekitException if some specific error occurs
*/
@Override
public List<ShortPeriodTerms> initialize(final AuxiliaryElements aux, final boolean meanOnly)
throws OrekitException {
// Initializes specific parameters.
initializeStep(aux);
// Truncation tolerance.
final double aor = a / R3;
final double tol = ( aor > .3 || (aor > .15 && ecc > .25) ) ? BIG_TRUNCATION_TOLERANCE : SMALL_TRUNCATION_TOLERANCE;
// Utilities for truncation
// Set a lower bound for eccentricity
final double eo2 = FastMath.max(0.0025, ecc / 2.);
final double x2o2 = XX / 2.;
final double[] eccPwr = new double[MAX_POWER];
final double[] chiPwr = new double[MAX_POWER];
eccPwr[0] = 1.;
chiPwr[0] = X;
for (int i = 1; i < MAX_POWER; i++) {
eccPwr[i] = eccPwr[i - 1] * eo2;
chiPwr[i] = chiPwr[i - 1] * x2o2;
}
// Auxiliary quantities.
final double ao2rxx = aor / (2. * XX);
double xmuarn = ao2rxx * ao2rxx * gm / (X * R3);
double term = 0.;
// Compute max power for a/R3 and e.
maxAR3Pow = 2;
maxEccPow = 0;
int n = 2;
int m = 2;
int nsmd2 = 0;
do {
// Upper bound for Tnm.
term = xmuarn *
(fact[n + m] / (fact[nsmd2] * fact[nsmd2 + m])) *
(fact[n + m + 1] / (fact[m] * fact[n + 1])) *
(fact[n - m + 1] / fact[n + 1]) *
eccPwr[m] * UpperBounds.getDnl(XX, chiPwr[m], n + 2, m);
if (term < tol) {
if (m == 0) {
break;
} else if (m < 2) {
xmuarn *= ao2rxx;
m = 0;
n++;
nsmd2++;
} else {
m -= 2;
nsmd2++;
}
} else {
maxAR3Pow = n;
maxEccPow = FastMath.max(m, maxEccPow);
xmuarn *= ao2rxx;
m++;
n++;
}
} while (n < MAX_POWER);
maxEccPow = FastMath.min(maxAR3Pow, maxEccPow);
// allocate the array aoR3Pow
aoR3Pow = new double[maxAR3Pow + 1];
maxFreqF = maxAR3Pow + 1;
maxEccPowShort = MAX_ECCPOWER_SP;
Qns = CoefficientsFactory.computeQns(gamma, maxAR3Pow, FastMath.max(maxEccPow, maxEccPowShort));
shortPeriods = new ThirdBodyShortPeriodicCoefficients(maxAR3Pow + 1, INTERPOLATION_POINTS,
maxFreqF, body.getName());
final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
list.add(shortPeriods);
return list;
}
/** {@inheritDoc} */
@Override
public void initializeStep(final AuxiliaryElements aux) throws OrekitException {
// Equinoctial elements
a = aux.getSma();
k = aux.getK();
h = aux.getH();
q = aux.getQ();
p = aux.getP();
// Eccentricity
ecc = aux.getEcc();
// Distance from center of mass of the central body to the 3rd body
final Vector3D bodyPos = body.getPVCoordinates(aux.getDate(), aux.getFrame()).getPosition();
R3 = bodyPos.getNorm();
// Direction cosines
final Vector3D bodyDir = bodyPos.normalize();
alpha = bodyDir.dotProduct(aux.getVectorF());
beta = bodyDir.dotProduct(aux.getVectorG());
gamma = bodyDir.dotProduct(aux.getVectorW());
// Equinoctial coefficients
A = aux.getA();
B = aux.getB();
C = aux.getC();
meanMotion = aux.getMeanMotion();
//Χ<sup>-2</sup>.
BB = B * B;
//Χ<sup>-3</sup>.
BBB = BB * B;
//b = 1 / (1 + B)
b = 1. / (1. + B);
// Χ
X = 1. / B;
XX = X * X;
XXX = X * XX;
// -2 * a / A
m2aoA = -2. * a / A;
// B / A
BoA = B / A;
// 1 / AB
ooAB = 1. / (A * B);
// -C / 2AB
mCo2AB = -C * ooAB / 2.;
// B / A(1 + B)
BoABpo = BoA / (1. + B);
// mu3 / R3
muoR3 = gm / R3;
//h * Χ³
hXXX = h * XXX;
//k * Χ³
kXXX = k * XXX;
}
/** {@inheritDoc} */
@Override
public double[] getMeanElementRate(final SpacecraftState currentState) {
// Qns coefficients
Qns = CoefficientsFactory.computeQns(gamma, maxAR3Pow, maxEccPow);
// a / R3 up to power maxAR3Pow
final double aoR3 = a / R3;
aoR3Pow[0] = 1.;
for (int i = 1; i <= maxAR3Pow; i++) {
aoR3Pow[i] = aoR3 * aoR3Pow[i - 1];
}
// Compute potential U derivatives
final double[] dU = computeUDerivatives();
final double dUda = dU[0];
final double dUdk = dU[1];
final double dUdh = dU[2];
final double dUdAl = dU[3];
final double dUdBe = dU[4];
final double dUdGa = dU[5];
// Compute cross derivatives [Eq. 2.2-(8)]
// U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
final double UAlphaGamma = alpha * dUdGa - gamma * dUdAl;
// U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
final double UBetaGamma = beta * dUdGa - gamma * dUdBe;
// Common factor
final double pUAGmIqUBGoAB = (p * UAlphaGamma - I * q * UBetaGamma) * ooAB;
// Compute mean elements rates [Eq. 3.1-(1)]
final double da = 0.;
final double dh = BoA * dUdk + k * pUAGmIqUBGoAB;
final double dk = -BoA * dUdh - h * pUAGmIqUBGoAB;
final double dp = mCo2AB * UBetaGamma;
final double dq = mCo2AB * UAlphaGamma * I;
final double dM = m2aoA * dUda + BoABpo * (h * dUdh + k * dUdk) + pUAGmIqUBGoAB;
return new double[] {da, dk, dh, dq, dp, dM};
}
/** {@inheritDoc} */
@Override
public void updateShortPeriodTerms(final SpacecraftState ... meanStates)
throws OrekitException {
final Slot slot = shortPeriods.createSlot(meanStates);
for (final SpacecraftState meanState : meanStates) {
initializeStep(new AuxiliaryElements(meanState.getOrbit(), I));
// Qns coefficients
Qns = CoefficientsFactory.computeQns(gamma, maxAR3Pow, FastMath.max(maxEccPow, maxEccPowShort));
final GeneratingFunctionCoefficients gfCoefs =
new GeneratingFunctionCoefficients(maxAR3Pow, MAX_ECCPOWER_SP, maxAR3Pow + 1);
// a / R3 up to power maxAR3Pow
final double aoR3 = a / R3;
aoR3Pow[0] = 1.;
for (int i = 1; i <= maxAR3Pow; i++) {
aoR3Pow[i] = aoR3 * aoR3Pow[i - 1];
}
//Compute additional quantities
// 2 * a / An
final double ax2oAn = -m2aoA / meanMotion;
// B / An
final double BoAn = BoA / meanMotion;
// 1 / ABn
final double ooABn = ooAB / meanMotion;
// C / 2ABn
final double Co2ABn = -mCo2AB / meanMotion;
// B / (A * (1 + B) * n)
final double BoABpon = BoABpo / meanMotion;
// -3 / n²a² = -3 / nA
final double m3onA = -3 / (A * meanMotion);
//Compute the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients.
for (int j = 1; j < slot.cij.length; j++) {
// First compute the C<sub>i</sub><sup>j</sup> coefficients
final double[] currentCij = new double[6];
// Compute the cross derivatives operator :
final double SAlphaGammaCj = alpha * gfCoefs.getdSdgammaCj(j) - gamma * gfCoefs.getdSdalphaCj(j);
final double SAlphaBetaCj = alpha * gfCoefs.getdSdbetaCj(j) - beta * gfCoefs.getdSdalphaCj(j);
final double SBetaGammaCj = beta * gfCoefs.getdSdgammaCj(j) - gamma * gfCoefs.getdSdbetaCj(j);
final double ShkCj = h * gfCoefs.getdSdkCj(j) - k * gfCoefs.getdSdhCj(j);
final double pSagmIqSbgoABnCj = (p * SAlphaGammaCj - I * q * SBetaGammaCj) * ooABn;
final double ShkmSabmdSdlCj = ShkCj - SAlphaBetaCj - gfCoefs.getdSdlambdaCj(j);
currentCij[0] = ax2oAn * gfCoefs.getdSdlambdaCj(j);
currentCij[1] = -(BoAn * gfCoefs.getdSdhCj(j) + h * pSagmIqSbgoABnCj + k * BoABpon * gfCoefs.getdSdlambdaCj(j));
currentCij[2] = BoAn * gfCoefs.getdSdkCj(j) + k * pSagmIqSbgoABnCj - h * BoABpon * gfCoefs.getdSdlambdaCj(j);
currentCij[3] = Co2ABn * (q * ShkmSabmdSdlCj - I * SAlphaGammaCj);
currentCij[4] = Co2ABn * (p * ShkmSabmdSdlCj - SBetaGammaCj);
currentCij[5] = -ax2oAn * gfCoefs.getdSdaCj(j) + BoABpon * (h * gfCoefs.getdSdhCj(j) + k * gfCoefs.getdSdkCj(j)) + pSagmIqSbgoABnCj + m3onA * gfCoefs.getSCj(j);
// add the computed coefficients to the interpolators
slot.cij[j].addGridPoint(meanState.getDate(), currentCij);
// Compute the S<sub>i</sub><sup>j</sup> coefficients
final double[] currentSij = new double[6];
// Compute the cross derivatives operator :
final double SAlphaGammaSj = alpha * gfCoefs.getdSdgammaSj(j) - gamma * gfCoefs.getdSdalphaSj(j);
final double SAlphaBetaSj = alpha * gfCoefs.getdSdbetaSj(j) - beta * gfCoefs.getdSdalphaSj(j);
final double SBetaGammaSj = beta * gfCoefs.getdSdgammaSj(j) - gamma * gfCoefs.getdSdbetaSj(j);
final double ShkSj = h * gfCoefs.getdSdkSj(j) - k * gfCoefs.getdSdhSj(j);
final double pSagmIqSbgoABnSj = (p * SAlphaGammaSj - I * q * SBetaGammaSj) * ooABn;
final double ShkmSabmdSdlSj = ShkSj - SAlphaBetaSj - gfCoefs.getdSdlambdaSj(j);
currentSij[0] = ax2oAn * gfCoefs.getdSdlambdaSj(j);
currentSij[1] = -(BoAn * gfCoefs.getdSdhSj(j) + h * pSagmIqSbgoABnSj + k * BoABpon * gfCoefs.getdSdlambdaSj(j));
currentSij[2] = BoAn * gfCoefs.getdSdkSj(j) + k * pSagmIqSbgoABnSj - h * BoABpon * gfCoefs.getdSdlambdaSj(j);
currentSij[3] = Co2ABn * (q * ShkmSabmdSdlSj - I * SAlphaGammaSj);
currentSij[4] = Co2ABn * (p * ShkmSabmdSdlSj - SBetaGammaSj);
currentSij[5] = -ax2oAn * gfCoefs.getdSdaSj(j) + BoABpon * (h * gfCoefs.getdSdhSj(j) + k * gfCoefs.getdSdkSj(j)) + pSagmIqSbgoABnSj + m3onA * gfCoefs.getSSj(j);
// add the computed coefficients to the interpolators
slot.sij[j].addGridPoint(meanState.getDate(), currentSij);
if (j == 1) {
//Compute the C⁰ coefficients using Danielson 2.5.2-15a.
final double[] value = new double[6];
for (int i = 0; i < 6; ++i) {
value[i] = currentCij[i] * k / 2. + currentSij[i] * h / 2.;
}
slot.cij[0].addGridPoint(meanState.getDate(), value);
}
}
}
}
/** {@inheritDoc} */
@Override
public EventDetector[] getEventsDetectors() {
return null;
}
/** Compute potential derivatives.
* @return derivatives of the potential with respect to orbital parameters
*/
private double[] computeUDerivatives() {
// Gs and Hs coefficients
final double[][] GsHs = CoefficientsFactory.computeGsHs(k, h, alpha, beta, maxEccPow);
// Initialise U.
U = 0.;
// Potential derivatives
double dUda = 0.;
double dUdk = 0.;
double dUdh = 0.;
double dUdAl = 0.;
double dUdBe = 0.;
double dUdGa = 0.;
for (int s = 0; s <= maxEccPow; s++) {
// initialise the Hansen roots
this.hansenObjects[s].computeInitValues(B, BB, BBB);
// Get the current Gs coefficient
final double gs = GsHs[0][s];
// Compute Gs partial derivatives from 3.1-(9)
double dGsdh = 0.;
double dGsdk = 0.;
double dGsdAl = 0.;
double dGsdBe = 0.;
if (s > 0) {
// First get the G(s-1) and the H(s-1) coefficients
final double sxGsm1 = s * GsHs[0][s - 1];
final double sxHsm1 = s * GsHs[1][s - 1];
// Then compute derivatives
dGsdh = beta * sxGsm1 - alpha * sxHsm1;
dGsdk = alpha * sxGsm1 + beta * sxHsm1;
dGsdAl = k * sxGsm1 - h * sxHsm1;
dGsdBe = h * sxGsm1 + k * sxHsm1;
}
// Kronecker symbol (2 - delta(0,s))
final double delta0s = (s == 0) ? 1. : 2.;
for (int n = FastMath.max(2, s); n <= maxAR3Pow; n++) {
// (n - s) must be even
if ((n - s) % 2 == 0) {
// Extract data from previous computation :
final double kns = this.hansenObjects[s].getValue(n, B);
final double dkns = this.hansenObjects[s].getDerivative(n, B);
final double vns = Vns.get(new NSKey(n, s));
final double coef0 = delta0s * aoR3Pow[n] * vns;
final double coef1 = coef0 * Qns[n][s];
final double coef2 = coef1 * kns;
// dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
// for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
final double dqns = (n == s) ? 0. : Qns[n][s + 1];
//Compute U:
U += coef2 * gs;
// Compute dU / da :
dUda += coef2 * n * gs;
// Compute dU / dh
dUdh += coef1 * (kns * dGsdh + hXXX * gs * dkns);
// Compute dU / dk
dUdk += coef1 * (kns * dGsdk + kXXX * gs * dkns);
// Compute dU / dAlpha
dUdAl += coef2 * dGsdAl;
// Compute dU / dBeta
dUdBe += coef2 * dGsdBe;
// Compute dU / dGamma
dUdGa += coef0 * kns * dqns * gs;
}
}
}
// multiply by mu3 / R3
U *= muoR3;
return new double[] {
dUda * muoR3 / a,
dUdk * muoR3,
dUdh * muoR3,
dUdAl * muoR3,
dUdBe * muoR3,
dUdGa * muoR3
};
}
/** {@inheritDoc} */
@Override
public void registerAttitudeProvider(final AttitudeProvider provider) {
//nothing is done since this contribution is not sensitive to attitude
}
/** Computes the C<sup>j</sup> and S<sup>j</sup> coefficients Danielson 4.2-(15,16)
* and their derivatives.
* <p>
* CS Mathematical Report $3.5.3.2
* </p>
*/
private class FourierCjSjCoefficients {
/** The coefficients G<sub>n, s</sub> and their derivatives. */
private final GnsCoefficients gns;
/** the coefficients e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and their derivatives by h and k. */
private final WnsjEtomjmsCoefficient wnsjEtomjmsCoefficient;
/** The terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h). */
private final CjSjAlphaBetaKH ABDECoefficients;
/** The Fourier coefficients C<sup>j</sup> and their derivatives.
* <p>
* Each column of the matrix contains the following values: <br/>
* - C<sup>j</sup> <br/>
* - dC<sup>j</sup> / da <br/>
* - dC<sup>j</sup> / dk <br/>
* - dC<sup>j</sup> / dh <br/>
* - dC<sup>j</sup> / dα <br/>
* - dC<sup>j</sup> / dβ <br/>
* - dC<sup>j</sup> / dγ <br/>
* </p>
*/
private final double[][] cj;
/** The S<sup>j</sup> coefficients and their derivatives.
* <p>
* Each column of the matrix contains the following values: <br/>
* - S<sup>j</sup> <br/>
* - dS<sup>j</sup> / da <br/>
* - dS<sup>j</sup> / dk <br/>
* - dS<sup>j</sup> / dh <br/>
* - dS<sup>j</sup> / dα <br/>
* - dS<sup>j</sup> / dβ <br/>
* - dS<sup>j</sup> / dγ <br/>
* </p>
*/
private final double[][] sj;
/** The Coefficients C<sup>j</sup><sub>,λ</sub>.
* <p>
* See Danielson 4.2-21
* </p>
*/
private final double[] cjlambda;
/** The Coefficients S<sup>j</sup><sub>,λ</sub>.
* <p>
* See Danielson 4.2-21
* </p>
*/
private final double[] sjlambda;
/** Maximum value for n. */
private final int nMax;
/** Maximum value for s. */
private final int sMax;
/** Maximum value for j. */
private final int jMax;
/**
* Private constructor.
*
* @param nMax maximum value for n index
* @param sMax maximum value for s index
* @param jMax maximum value for j index
*/
FourierCjSjCoefficients(final int nMax, final int sMax, final int jMax) {
//Save parameters
this.nMax = nMax;
this.sMax = sMax;
this.jMax = jMax;
//Create objects
wnsjEtomjmsCoefficient = new WnsjEtomjmsCoefficient();
ABDECoefficients = new CjSjAlphaBetaKH();
gns = new GnsCoefficients(nMax, sMax);
//create arays
this.cj = new double[7][jMax + 1];
this.sj = new double[7][jMax + 1];
this.cjlambda = new double[jMax];
this.sjlambda = new double[jMax];
computeCoefficients();
}
/**
* Compute all coefficients.
*/
private void computeCoefficients() {
for (int j = 1; j <= jMax; j++) {
// initialise the coefficients
for (int i = 0; i <= 6; i++) {
cj[i][j] = 0.;
sj[i][j] = 0.;
}
if (j < jMax) {
// initialise the C<sup>j</sup><sub>,λ</sub> and S<sup>j</sup><sub>,λ</sub> coefficients
cjlambda[j] = 0.;
sjlambda[j] = 0.;
}
for (int s = 0; s <= sMax; s++) {
// Compute the coefficients A<sub>j, s</sub>, B<sub>j, s</sub>, D<sub>j, s</sub> and E<sub>j, s</sub>
ABDECoefficients.computeCoefficients(j, s);
// compute starting value for n
final int minN = FastMath.max(2, FastMath.max(j - 1, s));
for (int n = minN; n <= nMax; n++) {
// check if n-s is even
if ((n - s) % 2 == 0) {
// compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n+1, s</sup> and its derivatives
final double[] wjnp1semjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n + 1);
// compute the coefficient e<sup>-|j-s|</sup>*w<sub>-j</sub><sup>n+1, s</sup> and its derivatives
final double[] wmjnp1semjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(-j, s, n + 1);
// compute common factors
final double coef1 = -(wjnp1semjms[0] * ABDECoefficients.getCoefA() + wmjnp1semjms[0] * ABDECoefficients.getCoefB());
final double coef2 = wjnp1semjms[0] * ABDECoefficients.getCoefD() + wmjnp1semjms[0] * ABDECoefficients.getCoefE();
//Compute C<sup>j</sup>
cj[0][j] += gns.getGns(n, s) * coef1;
//Compute dC<sup>j</sup> / da
cj[1][j] += gns.getdGnsda(n, s) * coef1;
//Compute dC<sup>j</sup> / dk
cj[2][j] += -gns.getGns(n, s) *
(
wjnp1semjms[1] * ABDECoefficients.getCoefA() +
wjnp1semjms[0] * ABDECoefficients.getdCoefAdk() +
wmjnp1semjms[1] * ABDECoefficients.getCoefB() +
wmjnp1semjms[0] * ABDECoefficients.getdCoefBdk()
);
//Compute dC<sup>j</sup> / dh
cj[3][j] += -gns.getGns(n, s) *
(
wjnp1semjms[2] * ABDECoefficients.getCoefA() +
wjnp1semjms[0] * ABDECoefficients.getdCoefAdh() +
wmjnp1semjms[2] * ABDECoefficients.getCoefB() +
wmjnp1semjms[0] * ABDECoefficients.getdCoefBdh()
);
//Compute dC<sup>j</sup> / dα
cj[4][j] += -gns.getGns(n, s) *
(
wjnp1semjms[0] * ABDECoefficients.getdCoefAdalpha() +
wmjnp1semjms[0] * ABDECoefficients.getdCoefBdalpha()
);
//Compute dC<sup>j</sup> / dβ
cj[5][j] += -gns.getGns(n, s) *
(
wjnp1semjms[0] * ABDECoefficients.getdCoefAdbeta() +
wmjnp1semjms[0] * ABDECoefficients.getdCoefBdbeta()
);
//Compute dC<sup>j</sup> / dγ
cj[6][j] += gns.getdGnsdgamma(n, s) * coef1;
//Compute S<sup>j</sup>
sj[0][j] += gns.getGns(n, s) * coef2;
//Compute dS<sup>j</sup> / da
sj[1][j] += gns.getdGnsda(n, s) * coef2;
//Compute dS<sup>j</sup> / dk
sj[2][j] += gns.getGns(n, s) *
(
wjnp1semjms[1] * ABDECoefficients.getCoefD() +
wjnp1semjms[0] * ABDECoefficients.getdCoefDdk() +
wmjnp1semjms[1] * ABDECoefficients.getCoefE() +
wmjnp1semjms[0] * ABDECoefficients.getdCoefEdk()
);
//Compute dS<sup>j</sup> / dh
sj[3][j] += gns.getGns(n, s) *
(
wjnp1semjms[2] * ABDECoefficients.getCoefD() +
wjnp1semjms[0] * ABDECoefficients.getdCoefDdh() +
wmjnp1semjms[2] * ABDECoefficients.getCoefE() +
wmjnp1semjms[0] * ABDECoefficients.getdCoefEdh()
);
//Compute dS<sup>j</sup> / dα
sj[4][j] += gns.getGns(n, s) *
(
wjnp1semjms[0] * ABDECoefficients.getdCoefDdalpha() +
wmjnp1semjms[0] * ABDECoefficients.getdCoefEdalpha()
);
//Compute dS<sup>j</sup> / dβ
sj[5][j] += gns.getGns(n, s) *
(
wjnp1semjms[0] * ABDECoefficients.getdCoefDdbeta() +
wmjnp1semjms[0] * ABDECoefficients.getdCoefEdbeta()
);
//Compute dS<sup>j</sup> / dγ
sj[6][j] += gns.getdGnsdgamma(n, s) * coef2;
//Check if n is greater or equal to j and j is at most jMax-1
if (n >= j && j < jMax) {
// compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives
final double[] wjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n);
// compute the coefficient e<sup>-|j-s|</sup>*w<sub>-j</sub><sup>n, s</sup> and its derivatives
final double[] wmjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(-j, s, n);
//Compute C<sup>j</sup><sub>,λ</sub>
cjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefD() + wmjnsemjms[0] * ABDECoefficients.getCoefE());
//Compute S<sup>j</sup><sub>,λ</sub>
sjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefA() + wmjnsemjms[0] * ABDECoefficients.getCoefB());
}
}
}
}
// Divide by j
for (int i = 0; i <= 6; i++) {
cj[i][j] /= j;
sj[i][j] /= j;
}
}
//The C⁰ coefficients are not computed here.
//They are evaluated at the final point.
//C⁰<sub>,λ</sub>
cjlambda[0] = k * cjlambda[1] / 2. + h * sjlambda[1] / 2.;
}
/** Get the Fourier coefficient C<sup>j</sup>.
* @param j j index
* @return C<sup>j</sup>
*/
public double getCj(final int j) {
return cj[0][j];
}
/** Get the derivative dC<sup>j</sup>/da.
* @param j j index
* @return dC<sup>j</sup>/da
*/
public double getdCjda(final int j) {
return cj[1][j];
}
/** Get the derivative dC<sup>j</sup>/dk.
* @param j j index
* @return dC<sup>j</sup>/dk
*/
public double getdCjdk(final int j) {
return cj[2][j];
}
/** Get the derivative dC<sup>j</sup>/dh.
* @param j j index
* @return dC<sup>j</sup>/dh
*/
public double getdCjdh(final int j) {
return cj[3][j];
}
/** Get the derivative dC<sup>j</sup>/dα.
* @param j j index
* @return dC<sup>j</sup>/dα
*/
public double getdCjdalpha(final int j) {
return cj[4][j];
}
/** Get the derivative dC<sup>j</sup>/dβ.
* @param j j index
* @return dC<sup>j</sup>/dβ
*/
public double getdCjdbeta(final int j) {
return cj[5][j];
}
/** Get the derivative dC<sup>j</sup>/dγ.
* @param j j index
* @return dC<sup>j</sup>/dγ
*/
public double getdCjdgamma(final int j) {
return cj[6][j];
}
/** Get the Fourier coefficient S<sup>j</sup>.
* @param j j index
* @return S<sup>j</sup>
*/
public double getSj(final int j) {
return sj[0][j];
}
/** Get the derivative dS<sup>j</sup>/da.
* @param j j index
* @return dS<sup>j</sup>/da
*/
public double getdSjda(final int j) {
return sj[1][j];
}
/** Get the derivative dS<sup>j</sup>/dk.
* @param j j index
* @return dS<sup>j</sup>/dk
*/
public double getdSjdk(final int j) {
return sj[2][j];
}
/** Get the derivative dS<sup>j</sup>/dh.
* @param j j index
* @return dS<sup>j</sup>/dh
*/
public double getdSjdh(final int j) {
return sj[3][j];
}
/** Get the derivative dS<sup>j</sup>/dα.
* @param j j index
* @return dS<sup>j</sup>/dα
*/
public double getdSjdalpha(final int j) {
return sj[4][j];
}
/** Get the derivative dS<sup>j</sup>/dβ.
* @param j j index
* @return dS<sup>j</sup>/dβ
*/
public double getdSjdbeta(final int j) {
return sj[5][j];
}
/** Get the derivative dS<sup>j</sup>/dγ.
* @param j j index
* @return dS<sup>j</sup>/dγ
*/
public double getdSjdgamma(final int j) {
return sj[6][j];
}
/** Get the coefficient C⁰<sub>,λ</sub>.
* @return C⁰<sub>,λ</sub>
*/
public double getC0Lambda() {
return cjlambda[0];
}
/** Get the coefficient C<sup>j</sup><sub>,λ</sub>.
* @param j j index
* @return C<sup>j</sup><sub>,λ</sub>
*/
public double getCjLambda(final int j) {
if (j < 1 || j >= jMax) {
return 0.;
}
return cjlambda[j];
}
/** Get the coefficient S<sup>j</sup><sub>,λ</sub>.
* @param j j index
* @return S<sup>j</sup><sub>,λ</sub>
*/
public double getSjLambda(final int j) {
if (j < 1 || j >= jMax) {
return 0.;
}
return sjlambda[j];
}
}
/** This class covers the coefficients e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and their derivatives by h and k.
*
* <p>
* Starting from Danielson 4.2-9,10,11 and taking into account that fact that: <br />
* c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
* the expression e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup>
* can be written as: <br/ >
* - for |s| > |j| <br />
* e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
* (((n + s)!(n - s)!)/((n + j)!(n - j)!)) *
* (-b)<sup>|j-s|</sup> *
* ((1 - c²)<sup>n-|s|</sup>/(1 + c²)<sup>n</sup>) *
* P<sub>n-|s|</sub><sup>|j-s|, |j+s|</sup>(χ) <br />
* <br />
* - for |s| <= |j| <br />
* e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
* (-b)<sup>|j-s|</sup> *
* ((1 - c²)<sup>n-|j|</sup>/(1 + c²)<sup>n</sup>) *
* P<sub>n-|j|</sub><sup>|j-s|, |j+s|</sup>(χ)
* </p>
*
* @author Lucian Barbulescu
*/
private class WnsjEtomjmsCoefficient {
/** The value c.
* <p>
* c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
* </p>
* */
private final double c;
/** c².*/
private final double c2;
/** db / dh. */
private final double dbdh;
/** db / dk. */
private final double dbdk;
/** de / dh. */
private final double dedh;
/** de / dk. */
private final double dedk;
/** dc / dh = e * db/dh + b * de/dh. */
private final double dcdh;
/** dc / dk = e * db/dk + b * de/dk. */
private final double dcdk;
/** The values (1 - c²)<sup>n</sup>. <br />
* The maximum possible value for the power is N + 1 */
private final double[] omc2tn;
/** The values (1 + c²)<sup>n</sup>. <br />
* The maximum possible value for the power is N + 1 */
private final double[] opc2tn;
/** The values b<sup>|j-s|</sup>. */
private final double[] btjms;
/**
* Standard constructor.
*/
WnsjEtomjmsCoefficient() {
//initialise fields
c = ecc * b;
c2 = c * c;
//b² * χ
final double b2Chi = b * b * X;
//Compute derivatives of b
dbdh = h * b2Chi;
dbdk = k * b2Chi;
//Compute derivatives of e
dedh = h / ecc;
dedk = k / ecc;
//Compute derivatives of c
dcdh = ecc * dbdh + b * dedh;
dcdk = ecc * dbdk + b * dedk;
//Compute the powers (1 - c²)<sup>n</sup> and (1 + c²)<sup>n</sup>
omc2tn = new double[maxAR3Pow + maxFreqF + 2];
opc2tn = new double[maxAR3Pow + maxFreqF + 2];
final double omc2 = 1. - c2;
final double opc2 = 1. + c2;
omc2tn[0] = 1.;
opc2tn[0] = 1.;
for (int i = 1; i <= maxAR3Pow + maxFreqF + 1; i++) {
omc2tn[i] = omc2tn[i - 1] * omc2;
opc2tn[i] = opc2tn[i - 1] * opc2;
}
//Compute the powers of b
btjms = new double[maxAR3Pow + maxFreqF + 1];
btjms[0] = 1.;
for (int i = 1; i <= maxAR3Pow + maxFreqF; i++) {
btjms[i] = btjms[i - 1] * b;
}
}
/** Compute the value of the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives by h and k. <br />
*
* @param j j index
* @param s s index
* @param n n index
* @return an array containing the value of the coefficient at index 0, the derivative by k at index 1 and the derivative by h at index 2
*/
public double[] computeWjnsEmjmsAndDeriv(final int j, final int s, final int n) {
final double[] wjnsemjms = new double[] {0., 0., 0.};
// |j|
final int absJ = FastMath.abs(j);
// |s|
final int absS = FastMath.abs(s);
// |j - s|
final int absJmS = FastMath.abs(j - s);
// |j + s|
final int absJpS = FastMath.abs(j + s);
//The lower index of P. Also the power of (1 - c²)
int l;
// the factorial ratio coefficient or 1. if |s| <= |j|
double factCoef;
if (absS > absJ) {
factCoef = (fact[n + s] / fact[n + j]) * (fact[n - s] / fact [n - j]);
l = n - absS;
} else {
factCoef = 1.;
l = n - absJ;
}
// (-1)<sup>|j-s|</sup>
final double sign = absJmS % 2 != 0 ? -1. : 1.;
//(1 - c²)<sup>n-|s|</sup> / (1 + c²)<sup>n</sup>
final double coef1 = omc2tn[l] / opc2tn[n];
//-b<sup>|j-s|</sup>
final double coef2 = sign * btjms[absJmS];
// P<sub>l</sub><sup>|j-s|, |j+s|</sup>(χ)
final DerivativeStructure jac =
JacobiPolynomials.getValue(l, absJmS, absJpS, new DerivativeStructure(1, 1, 0, X));
// the derivative of coef1 by c
final double dcoef1dc = -coef1 * 2. * c * (((double) n) / opc2tn[1] + ((double) l) / omc2tn[1]);
// the derivative of coef1 by h
final double dcoef1dh = dcoef1dc * dcdh;
// the derivative of coef1 by k
final double dcoef1dk = dcoef1dc * dcdk;
// the derivative of coef2 by b
final double dcoef2db = absJmS == 0 ? 0 : sign * (double) absJmS * btjms[absJmS - 1];
// the derivative of coef2 by h
final double dcoef2dh = dcoef2db * dbdh;
// the derivative of coef2 by k
final double dcoef2dk = dcoef2db * dbdk;
// the jacobi polinomial value
final double jacobi = jac.getValue();
// the derivative of the Jacobi polynomial by h
final double djacobidh = jac.getPartialDerivative(1) * hXXX;
// the derivative of the Jacobi polynomial by k
final double djacobidk = jac.getPartialDerivative(1) * kXXX;
//group the above coefficients to limit the mathematical operations
final double term1 = factCoef * coef1 * coef2;
final double term2 = factCoef * coef1 * jacobi;
final double term3 = factCoef * coef2 * jacobi;
//compute e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives by k and h
wjnsemjms[0] = term1 * jacobi;
wjnsemjms[1] = dcoef1dk * term3 + dcoef2dk * term2 + djacobidk * term1;
wjnsemjms[2] = dcoef1dh * term3 + dcoef2dh * term2 + djacobidh * term1;
return wjnsemjms;
}
}
/** The G<sub>n,s</sub> coefficients and their derivatives.
* <p>
* See Danielson, 4.2-17
*
* @author Lucian Barbulescu
*/
private class GnsCoefficients {
/** Maximum value for n index. */
private final int nMax;
/** Maximum value for s index. */
private final int sMax;
/** The coefficients G<sub>n,s</sub>. */
private final double gns[][];
/** The derivatives of the coefficients G<sub>n,s</sub> by a. */
private final double dgnsda[][];
/** The derivatives of the coefficients G<sub>n,s</sub> by γ. */
private final double dgnsdgamma[][];
/** Standard constructor.
*
* @param nMax maximum value for n indes
* @param sMax maximum value for s index
*/
GnsCoefficients(final int nMax, final int sMax) {
this.nMax = nMax;
this.sMax = sMax;
final int rows = nMax + 1;
final int columns = sMax + 1;
this.gns = new double[rows][columns];
this.dgnsda = new double[rows][columns];
this.dgnsdgamma = new double[rows][columns];
// Generate the coefficients
generateCoefficients();
}
/**
* Compute the coefficient G<sub>n,s</sub> and its derivatives.
* <p>
* Only the derivatives by a and γ are computed as all others are 0
* </p>
*/
private void generateCoefficients() {
for (int s = 0; s <= sMax; s++) {
// The n index is always at least the maximum between 2 and s
final int minN = FastMath.max(2, s);
for (int n = minN; n <= nMax; n++) {
// compute the coefficients only if (n - s) % 2 == 0
if ( (n - s) % 2 == 0 ) {
// Kronecker symbol (2 - delta(0,s))
final double delta0s = (s == 0) ? 1. : 2.;
final double vns = Vns.get(new NSKey(n, s));
final double coef0 = delta0s * aoR3Pow[n] * vns * muoR3;
final double coef1 = coef0 * Qns[n][s];
// dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
// for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
final double dqns = (n == s) ? 0. : Qns[n][s + 1];
//Compute the coefficient and its derivatives.
this.gns[n][s] = coef1;
this.dgnsda[n][s] = coef1 * n / a;
this.dgnsdgamma[n][s] = coef0 * dqns;
} else {
// the coefficient and its derivatives is 0
this.gns[n][s] = 0.;
this.dgnsda[n][s] = 0.;
this.dgnsdgamma[n][s] = 0.;
}
}
}
}
/** Get the coefficient G<sub>n,s</sub>.
*
* @param n n index
* @param s s index
* @return the coefficient G<sub>n,s</sub>
*/
public double getGns(final int n, final int s) {
return this.gns[n][s];
}
/** Get the derivative dG<sub>n,s</sub> / da.
*
* @param n n index
* @param s s index
* @return the derivative dG<sub>n,s</sub> / da
*/
public double getdGnsda(final int n, final int s) {
return this.dgnsda[n][s];
}
/** Get the derivative dG<sub>n,s</sub> / dγ.
*
* @param n n index
* @param s s index
* @return the derivative dG<sub>n,s</sub> / dγ
*/
public double getdGnsdgamma(final int n, final int s) {
return this.dgnsdgamma[n][s];
}
}
/** This class computes the terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h).
*
* <p>
* The following terms and their derivatives by k, h, alpha and beta are considered: <br/ >
* - sign(j-s) * C<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) + S<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) <br />
* - C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) <br />
* - C<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) - sign(j-s) * S<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) <br />
* - C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) <br />
* For the ease of usage the above terms are renamed A<sub>js</sub>, B<sub>js</sub>, D<sub>js</sub> and E<sub>js</sub> respectively <br />
* See the CS Mathematical report $3.5.3.2 for more details
* </p>
* @author Lucian Barbulescu
*/
private class CjSjAlphaBetaKH {
/** The C<sub>j</sub>(k, h) and the S<sub>j</sub>(k, h) series. */
private final CjSjCoefficient cjsjkh;
/** The C<sub>j</sub>(α, β) and the S<sub>j</sub>(α, β) series. */
private final CjSjCoefficient cjsjalbe;
/** The coeficient sign(j-s) * C<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) + S<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h)
* and its derivative by k, h, α and β. */
private final double coefAandDeriv[];
/** The coeficient C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h)
* and its derivative by k, h, α and β. */
private final double coefBandDeriv[];
/** The coeficient C<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) - sign(j-s) * S<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h)
* and its derivative by k, h, α and β. */
private final double coefDandDeriv[];
/** The coeficient C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h)
* and its derivative by k, h, α and β. */
private final double coefEandDeriv[];
/**
* Standard constructor.
*/
CjSjAlphaBetaKH() {
cjsjkh = new CjSjCoefficient(k, h);
cjsjalbe = new CjSjCoefficient(alpha, beta);
coefAandDeriv = new double[5];
coefBandDeriv = new double[5];
coefDandDeriv = new double[5];
coefEandDeriv = new double[5];
}
/** Compute the coefficients and their derivatives for a given (j,s) pair.
* @param j j index
* @param s s index
*/
public void computeCoefficients(final int j, final int s) {
// sign of j-s
final int sign = j < s ? -1 : 1;
//|j-s|
final int absJmS = FastMath.abs(j - s);
//j+s
final int jps = j + s;
//Compute the coefficient A and its derivatives
coefAandDeriv[0] = sign * cjsjalbe.getCj(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getCj(absJmS);
coefAandDeriv[1] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDk(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDk(absJmS);
coefAandDeriv[2] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDh(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDh(absJmS);
coefAandDeriv[3] = sign * cjsjalbe.getDcjDk(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDk(s) * cjsjkh.getCj(absJmS);
coefAandDeriv[4] = sign * cjsjalbe.getDcjDh(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDh(s) * cjsjkh.getCj(absJmS);
//Compute the coefficient B and its derivatives
coefBandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getSj(jps) - cjsjalbe.getSj(s) * cjsjkh.getCj(jps);
coefBandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDsjDk(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDk(jps);
coefBandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDsjDh(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDh(jps);
coefBandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDk(s) * cjsjkh.getCj(jps);
coefBandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDh(s) * cjsjkh.getCj(jps);
//Compute the coefficient D and its derivatives
coefDandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getSj(absJmS);
coefDandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDk(absJmS);
coefDandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDh(absJmS);
coefDandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDk(s) * cjsjkh.getSj(absJmS);
coefDandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDh(s) * cjsjkh.getSj(absJmS);
//Compute the coefficient E and its derivatives
coefEandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(jps) + cjsjalbe.getSj(s) * cjsjkh.getSj(jps);
coefEandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDk(jps);
coefEandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDh(jps);
coefEandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDk(s) * cjsjkh.getSj(jps);
coefEandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDh(s) * cjsjkh.getSj(jps);
}
/** Get the value of coefficient A<sub>j,s</sub>.
*
* @return the coefficient A<sub>j,s</sub>
*/
public double getCoefA() {
return coefAandDeriv[0];
}
/** Get the value of coefficient dA<sub>j,s</sub>/dk.
*
* @return the coefficient dA<sub>j,s</sub>/dk
*/
public double getdCoefAdk() {
return coefAandDeriv[1];
}
/** Get the value of coefficient dA<sub>j,s</sub>/dh.
*
* @return the coefficient dA<sub>j,s</sub>/dh
*/
public double getdCoefAdh() {
return coefAandDeriv[2];
}
/** Get the value of coefficient dA<sub>j,s</sub>/dα.
*
* @return the coefficient dA<sub>j,s</sub>/dα
*/
public double getdCoefAdalpha() {
return coefAandDeriv[3];
}
/** Get the value of coefficient dA<sub>j,s</sub>/dβ.
*
* @return the coefficient dA<sub>j,s</sub>/dβ
*/
public double getdCoefAdbeta() {
return coefAandDeriv[4];
}
/** Get the value of coefficient B<sub>j,s</sub>.
*
* @return the coefficient B<sub>j,s</sub>
*/
public double getCoefB() {
return coefBandDeriv[0];
}
/** Get the value of coefficient dB<sub>j,s</sub>/dk.
*
* @return the coefficient dB<sub>j,s</sub>/dk
*/
public double getdCoefBdk() {
return coefBandDeriv[1];
}
/** Get the value of coefficient dB<sub>j,s</sub>/dh.
*
* @return the coefficient dB<sub>j,s</sub>/dh
*/
public double getdCoefBdh() {
return coefBandDeriv[2];
}
/** Get the value of coefficient dB<sub>j,s</sub>/dα.
*
* @return the coefficient dB<sub>j,s</sub>/dα
*/
public double getdCoefBdalpha() {
return coefBandDeriv[3];
}
/** Get the value of coefficient dB<sub>j,s</sub>/dβ.
*
* @return the coefficient dB<sub>j,s</sub>/dβ
*/
public double getdCoefBdbeta() {
return coefBandDeriv[4];
}
/** Get the value of coefficient D<sub>j,s</sub>.
*
* @return the coefficient D<sub>j,s</sub>
*/
public double getCoefD() {
return coefDandDeriv[0];
}
/** Get the value of coefficient dD<sub>j,s</sub>/dk.
*
* @return the coefficient dD<sub>j,s</sub>/dk
*/
public double getdCoefDdk() {
return coefDandDeriv[1];
}
/** Get the value of coefficient dD<sub>j,s</sub>/dh.
*
* @return the coefficient dD<sub>j,s</sub>/dh
*/
public double getdCoefDdh() {
return coefDandDeriv[2];
}
/** Get the value of coefficient dD<sub>j,s</sub>/dα.
*
* @return the coefficient dD<sub>j,s</sub>/dα
*/
public double getdCoefDdalpha() {
return coefDandDeriv[3];
}
/** Get the value of coefficient dD<sub>j,s</sub>/dβ.
*
* @return the coefficient dD<sub>j,s</sub>/dβ
*/
public double getdCoefDdbeta() {
return coefDandDeriv[4];
}
/** Get the value of coefficient E<sub>j,s</sub>.
*
* @return the coefficient E<sub>j,s</sub>
*/
public double getCoefE() {
return coefEandDeriv[0];
}
/** Get the value of coefficient dE<sub>j,s</sub>/dk.
*
* @return the coefficient dE<sub>j,s</sub>/dk
*/
public double getdCoefEdk() {
return coefEandDeriv[1];
}
/** Get the value of coefficient dE<sub>j,s</sub>/dh.
*
* @return the coefficient dE<sub>j,s</sub>/dh
*/
public double getdCoefEdh() {
return coefEandDeriv[2];
}
/** Get the value of coefficient dE<sub>j,s</sub>/dα.
*
* @return the coefficient dE<sub>j,s</sub>/dα
*/
public double getdCoefEdalpha() {
return coefEandDeriv[3];
}
/** Get the value of coefficient dE<sub>j,s</sub>/dβ.
*
* @return the coefficient dE<sub>j,s</sub>/dβ
*/
public double getdCoefEdbeta() {
return coefEandDeriv[4];
}
}
/** This class computes the coefficients for the generating function S and its derivatives.
* <p>
* The form of the generating functions is: <br>
* S = C⁰ + Σ<sub>j=1</sub><sup>N+1</sup>(C<sup>j</sup> * cos(jF) + S<sup>j</sup> * sin(jF)) <br>
* The coefficients C⁰, C<sup>j</sup>, S<sup>j</sup> are the Fourrier coefficients
* presented in Danielson 4.2-14,15 except for the case j=1 where
* C¹ = C¹<sub>Fourier</sub> - hU and
* S¹ = S¹<sub>Fourier</sub> + kU <br>
* Also the coefficients of the derivatives of S by a, k, h, α, β, γ and λ
* are computed end expressed in a similar manner. The formulas used are 4.2-19, 20, 23, 24
* </p>
* @author Lucian Barbulescu
*/
private class GeneratingFunctionCoefficients {
/** The Fourier coefficients as presented in Danielson 4.2-14,15. */
private final FourierCjSjCoefficients cjsjFourier;
/** Maximum value of j index. */
private final int jMax;
/** The coefficients C<sup>j</sup> of the function S and its derivatives.
* <p>
* The index j belongs to the interval [0,jMax]. The coefficient C⁰ is the free coefficient.<br>
* Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
* - S <br/>
* - dS / da <br/>
* - dS / dk <br/>
* - dS / dh <br/>
* - dS / dα <br/>
* - dS / dβ <br/>
* - dS / dγ <br/>
* - dS / dλ
* </p>
*/
private final double[][] cjCoefs;
/** The coefficients S<sup>j</sup> of the function S and its derivatives.
* <p>
* The index j belongs to the interval [0,jMax].<br>
* Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
* - S <br/>
* - dS / da <br/>
* - dS / dk <br/>
* - dS / dh <br/>
* - dS / dα <br/>
* - dS / dβ <br/>
* - dS / dγ <br/>
* - dS / dλ
* </p>
*/
private final double[][] sjCoefs;
/**
* Standard constructor.
*
* @param nMax maximum value of n index
* @param sMax maximum value of s index
* @param jMax maximum value of j index
*/
GeneratingFunctionCoefficients(final int nMax, final int sMax, final int jMax) {
this.jMax = jMax;
this.cjsjFourier = new FourierCjSjCoefficients(nMax, sMax, jMax);
this.cjCoefs = new double[8][jMax + 1];
this.sjCoefs = new double[8][jMax + 1];
computeGeneratingFunctionCoefficients();
}
/**
* Compute the coefficients for the generating function S and its derivatives.
*/
private void computeGeneratingFunctionCoefficients() {
// Compute potential U and derivatives
final double[] dU = computeUDerivatives();
//Compute the C<sup>j</sup> coefficients
for (int j = 1; j <= jMax; j++) {
//Compute the C<sup>j</sup> coefficients
cjCoefs[0][j] = cjsjFourier.getCj(j);
cjCoefs[1][j] = cjsjFourier.getdCjda(j);
cjCoefs[2][j] = cjsjFourier.getdCjdk(j) - (cjsjFourier.getSjLambda(j - 1) - cjsjFourier.getSjLambda(j + 1)) / 2;
cjCoefs[3][j] = cjsjFourier.getdCjdh(j) - (cjsjFourier.getCjLambda(j - 1) + cjsjFourier.getCjLambda(j + 1)) / 2;
cjCoefs[4][j] = cjsjFourier.getdCjdalpha(j);
cjCoefs[5][j] = cjsjFourier.getdCjdbeta(j);
cjCoefs[6][j] = cjsjFourier.getdCjdgamma(j);
cjCoefs[7][j] = cjsjFourier.getCjLambda(j);
//Compute the S<sup>j</sup> coefficients
sjCoefs[0][j] = cjsjFourier.getSj(j);
sjCoefs[1][j] = cjsjFourier.getdSjda(j);
sjCoefs[2][j] = cjsjFourier.getdSjdk(j) + (cjsjFourier.getCjLambda(j - 1) - cjsjFourier.getCjLambda(j + 1)) / 2;
sjCoefs[3][j] = cjsjFourier.getdSjdh(j) - (cjsjFourier.getSjLambda(j - 1) + cjsjFourier.getSjLambda(j + 1)) / 2;
sjCoefs[4][j] = cjsjFourier.getdSjdalpha(j);
sjCoefs[5][j] = cjsjFourier.getdSjdbeta(j);
sjCoefs[6][j] = cjsjFourier.getdSjdgamma(j);
sjCoefs[7][j] = cjsjFourier.getSjLambda(j);
//In the special case j == 1 there are some additional terms to be added
if (j == 1) {
//Additional terms for C<sup>j</sup> coefficients
cjCoefs[0][j] += -h * U;
cjCoefs[1][j] += -h * dU[0];
cjCoefs[2][j] += -h * dU[1];
cjCoefs[3][j] += -(h * dU[2] + U + cjsjFourier.getC0Lambda());
cjCoefs[4][j] += -h * dU[3];
cjCoefs[5][j] += -h * dU[4];
cjCoefs[6][j] += -h * dU[5];
//Additional terms for S<sup>j</sup> coefficients
sjCoefs[0][j] += k * U;
sjCoefs[1][j] += k * dU[0];
sjCoefs[2][j] += k * dU[1] + U + cjsjFourier.getC0Lambda();
sjCoefs[3][j] += k * dU[2];
sjCoefs[4][j] += k * dU[3];
sjCoefs[5][j] += k * dU[4];
sjCoefs[6][j] += k * dU[5];
}
}
}
/** Get the coefficient C<sup>j</sup> for the function S.
* <br>
* Possible values for j are within the interval [0,jMax].
* The value 0 is used to obtain the free coefficient C⁰
* @param j j index
* @return C<sup>j</sup> for the function S
*/
public double getSCj(final int j) {
return cjCoefs[0][j];
}
/** Get the coefficient S<sup>j</sup> for the function S.
* <br>
* Possible values for j are within the interval [1,jMax].
* @param j j index
* @return S<sup>j</sup> for the function S
*/
public double getSSj(final int j) {
return sjCoefs[0][j];
}
/** Get the coefficient C<sup>j</sup> for the derivative dS/da.
* <br>
* Possible values for j are within the interval [0,jMax].
* The value 0 is used to obtain the free coefficient C⁰
* @param j j index
* @return C<sup>j</sup> for the function dS/da
*/
public double getdSdaCj(final int j) {
return cjCoefs[1][j];
}
/** Get the coefficient S<sup>j</sup> for the derivative dS/da.
* <br>
* Possible values for j are within the interval [1,jMax].
* @param j j index
* @return S<sup>j</sup> for the derivative dS/da
*/
public double getdSdaSj(final int j) {
return sjCoefs[1][j];
}
/** Get the coefficient C<sup>j</sup> for the derivative dS/dk
* <br>
* Possible values for j are within the interval [0,jMax].
* The value 0 is used to obtain the free coefficient C⁰
* @param j j index
* @return C<sup>j</sup> for the function dS/dk
*/
public double getdSdkCj(final int j) {
return cjCoefs[2][j];
}
/** Get the coefficient S<sup>j</sup> for the derivative dS/dk.
* <br>
* Possible values for j are within the interval [1,jMax].
* @param j j index
* @return S<sup>j</sup> for the derivative dS/dk
*/
public double getdSdkSj(final int j) {
return sjCoefs[2][j];
}
/** Get the coefficient C<sup>j</sup> for the derivative dS/dh
* <br>
* Possible values for j are within the interval [0,jMax].
* The value 0 is used to obtain the free coefficient C⁰
* @param j j index
* @return C<sup>j</sup> for the function dS/dh
*/
public double getdSdhCj(final int j) {
return cjCoefs[3][j];
}
/** Get the coefficient S<sup>j</sup> for the derivative dS/dh.
* <br>
* Possible values for j are within the interval [1,jMax].
* @param j j index
* @return S<sup>j</sup> for the derivative dS/dh
*/
public double getdSdhSj(final int j) {
return sjCoefs[3][j];
}
/** Get the coefficient C<sup>j</sup> for the derivative dS/dα
* <br>
* Possible values for j are within the interval [0,jMax].
* The value 0 is used to obtain the free coefficient C⁰
* @param j j index
* @return C<sup>j</sup> for the function dS/dα
*/
public double getdSdalphaCj(final int j) {
return cjCoefs[4][j];
}
/** Get the coefficient S<sup>j</sup> for the derivative dS/dα.
* <br>
* Possible values for j are within the interval [1,jMax].
* @param j j index
* @return S<sup>j</sup> for the derivative dS/dα
*/
public double getdSdalphaSj(final int j) {
return sjCoefs[4][j];
}
/** Get the coefficient C<sup>j</sup> for the derivative dS/dβ
* <br>
* Possible values for j are within the interval [0,jMax].
* The value 0 is used to obtain the free coefficient C⁰
* @param j j index
* @return C<sup>j</sup> for the function dS/dβ
*/
public double getdSdbetaCj(final int j) {
return cjCoefs[5][j];
}
/** Get the coefficient S<sup>j</sup> for the derivative dS/dβ.
* <br>
* Possible values for j are within the interval [1,jMax].
* @param j j index
* @return S<sup>j</sup> for the derivative dS/dβ
*/
public double getdSdbetaSj(final int j) {
return sjCoefs[5][j];
}
/** Get the coefficient C<sup>j</sup> for the derivative dS/dγ
* <br>
* Possible values for j are within the interval [0,jMax].
* The value 0 is used to obtain the free coefficient C⁰
* @param j j index
* @return C<sup>j</sup> for the function dS/dγ
*/
public double getdSdgammaCj(final int j) {
return cjCoefs[6][j];
}
/** Get the coefficient S<sup>j</sup> for the derivative dS/dγ.
* <br>
* Possible values for j are within the interval [1,jMax].
* @param j j index
* @return S<sup>j</sup> for the derivative dS/dγ
*/
public double getdSdgammaSj(final int j) {
return sjCoefs[6][j];
}
/** Get the coefficient C<sup>j</sup> for the derivative dS/dλ
* <br>
* Possible values for j are within the interval [0,jMax].
* The value 0 is used to obtain the free coefficient C⁰
* @param j j index
* @return C<sup>j</sup> for the function dS/dλ
*/
public double getdSdlambdaCj(final int j) {
return cjCoefs[7][j];
}
/** Get the coefficient S<sup>j</sup> for the derivative dS/dλ.
* <br>
* Possible values for j are within the interval [1,jMax].
* @param j j index
* @return S<sup>j</sup> for the derivative dS/dλ
*/
public double getdSdlambdaSj(final int j) {
return sjCoefs[7][j];
}
}
/**
* The coefficients used to compute the short periodic contribution for the Third body perturbation.
* <p>
* The short periodic contribution for the Third Body is expressed in Danielson 4.2-25.<br>
* The coefficients C<sub>i</sub>⁰, C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup>
* are computed by replacing the corresponding values in formula 2.5.5-10.
* </p>
* @author Lucian Barbulescu
*/
private static class ThirdBodyShortPeriodicCoefficients implements ShortPeriodTerms {
/** Serializable UID. */
private static final long serialVersionUID = 20151119L;
/** Maximal value for j. */
private final int jMax;
/** Number of points used in the interpolation process. */
private final int interpolationPoints;
/** Max frequency of F. */
private final int maxFreqF;
/** Coefficients prefix. */
private final String prefix;
/** All coefficients slots. */
private final TimeSpanMap<Slot> slots;
/**
* Standard constructor.
* @param interpolationPoints number of points used in the interpolation process
* @param jMax maximal value for j
* @param maxFreqF Max frequency of F
* @param bodyName third body name
*/
ThirdBodyShortPeriodicCoefficients(final int jMax, final int interpolationPoints,
final int maxFreqF, final String bodyName) {
this.jMax = jMax;
this.interpolationPoints = interpolationPoints;
this.maxFreqF = maxFreqF;
this.prefix = "DSST-3rd-body-" + bodyName + "-";
this.slots = new TimeSpanMap<Slot>(new Slot(jMax, interpolationPoints));
}
/** Get the slot valid for some date.
* @param meanStates mean states defining the slot
* @return slot valid at the specified date
*/
public Slot createSlot(final SpacecraftState ... meanStates) {
final Slot slot = new Slot(jMax, interpolationPoints);
final AbsoluteDate first = meanStates[0].getDate();
final AbsoluteDate last = meanStates[meanStates.length - 1].getDate();
if (first.compareTo(last) <= 0) {
slots.addValidAfter(slot, first);
} else {
slots.addValidBefore(slot, first);
}
return slot;
}
/** {@inheritDoc} */
@Override
public double[] value(final Orbit meanOrbit) {
// select the coefficients slot
final Slot slot = slots.get(meanOrbit.getDate());
// the current eccentric longitude
final double F = meanOrbit.getLE();
//initialize the short periodic contribution with the corresponding C⁰ coeficient
final double[] shortPeriodic = slot.cij[0].value(meanOrbit.getDate());
// Add the cos and sin dependent terms
for (int j = 1; j <= maxFreqF; j++) {
//compute cos and sin
final double cosjF = FastMath.cos(j * F);
final double sinjF = FastMath.sin(j * F);
final double[] c = slot.cij[j].value(meanOrbit.getDate());
final double[] s = slot.sij[j].value(meanOrbit.getDate());
for (int i = 0; i < 6; i++) {
shortPeriodic[i] += c[i] * cosjF + s[i] * sinjF;
}
}
return shortPeriodic;
}
/** {@inheritDoc} */
@Override
public String getCoefficientsKeyPrefix() {
return prefix;
}
/** {@inheritDoc}
* <p>
* For third body attraction forces,there are maxFreqF + 1 cj coefficients,
* maxFreqF sj coefficients where maxFreqF depends on the orbit.
* The j index is the integer multiplier for the eccentric longitude argument
* in the cj and sj coefficients.
* </p>
*/
@Override
public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected)
throws OrekitException {
// select the coefficients slot
final Slot slot = slots.get(date);
final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * maxFreqF + 1);
storeIfSelected(coefficients, selected, slot.cij[0].value(date), "c", 0);
for (int j = 1; j <= maxFreqF; j++) {
storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
}
return coefficients;
}
/** Put a coefficient in a map if selected.
* @param map map to populate
* @param selected set of coefficients that should be put in the map
* (empty set means all coefficients are selected)
* @param value coefficient value
* @param id coefficient identifier
* @param indices list of coefficient indices
*/
private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
final double[] value, final String id, final int ... indices) {
final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
keyBuilder.append(id);
for (int index : indices) {
keyBuilder.append('[').append(index).append(']');
}
final String key = keyBuilder.toString();
if (selected.isEmpty() || selected.contains(key)) {
map.put(key, value);
}
}
}
/** Coefficients valid for one time slot. */
private static class Slot {
/** The coefficients C<sub>i</sub><sup>j</sup>.
* <p>
* The index order is cij[j][i] <br/>
* i corresponds to the equinoctial element, as follows: <br/>
* - i=0 for a <br/>
* - i=1 for k <br/>
* - i=2 for h <br/>
* - i=3 for q <br/>
* - i=4 for p <br/>
* - i=5 for λ <br/>
* </p>
*/
private final ShortPeriodicsInterpolatedCoefficient[] cij;
/** The coefficients S<sub>i</sub><sup>j</sup>.
* <p>
* The index order is sij[j][i] <br/>
* i corresponds to the equinoctial element, as follows: <br/>
* - i=0 for a <br/>
* - i=1 for k <br/>
* - i=2 for h <br/>
* - i=3 for q <br/>
* - i=4 for p <br/>
* - i=5 for λ <br/>
* </p>
*/
private final ShortPeriodicsInterpolatedCoefficient[] sij;
/** Simple constructor.
* @param jMax maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
*/
Slot(final int jMax, final int interpolationPoints) {
// allocate the coefficients arrays
cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
for (int j = 0; j <= jMax; j++) {
cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
}
}
}
}