DSSTZonal.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.io.NotSerializableException;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.SortedSet;
import java.util.TreeMap;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.FastMath;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.errors.OrekitException;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
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.GHIJjsPolynomials;
import org.orekit.propagation.semianalytical.dsst.utilities.LnsCoefficients;
import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenZonalLinear;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.TimeSpanMap;
/** Zonal contribution to the {@link DSSTCentralBody central body gravitational perturbation}.
*
* @author Romain Di Costanzo
* @author Pascal Parraud
*/
public class DSSTZonal implements DSSTForceModel {
/** Truncation tolerance. */
private static final double TRUNCATION_TOLERANCE = 1e-4;
/** Number of points for interpolation. */
private static final int INTERPOLATION_POINTS = 3;
/** 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;
/** Provider for spherical harmonics. */
private final UnnormalizedSphericalHarmonicsProvider provider;
/** Maximal degree to consider for harmonics potential. */
private final int maxDegree;
/** Maximal degree to consider for harmonics potential. */
private final int maxDegreeShortPeriodics;
/** Maximal degree to consider for harmonics potential in short periodic computations. */
private final int maxOrder;
/** Factorial. */
private final double[] fact;
/** Coefficient used to define the mean disturbing function V<sub>ns</sub> coefficient. */
private final TreeMap<NSKey, Double> Vns;
/** Highest power of the eccentricity to be used in mean elements computations. */
private int maxEccPowMeanElements;
/** Highest power of the eccentricity to be used in short periodic computations. */
private final int maxEccPowShortPeriodics;
/** Maximum frequency in true longitude for short periodic computations. */
private final int maxFrequencyShortPeriodics;
/** Short period terms. */
private ZonalShortPeriodicCoefficients zonalSPCoefs;
// Equinoctial elements (according to DSST notation)
/** a. */
private double a;
/** ex. */
private double k;
/** ey. */
private double h;
/** hx. */
private double q;
/** hy. */
private double p;
/** Eccentricity. */
private double ecc;
/** Direction cosine &alpha. */
private double alpha;
/** Direction cosine &beta. */
private double beta;
/** Direction cosine &gamma. */
private double gamma;
// Common factors for potential computation
/** Χ = 1 / sqrt(1 - e²) = 1 / B. */
private double X;
/** Χ². */
private double XX;
/** Χ³. */
private double XXX;
/** 1 / (A * B) .*/
private double ooAB;
/** B / A .*/
private double BoA;
/** B / A(1 + B) .*/
private double BoABpo;
/** -C / (2 * A * B) .*/
private double mCo2AB;
/** -2 * a / A .*/
private double m2aoA;
/** μ / a .*/
private double muoa;
/** R / a .*/
private double roa;
/** An array that contains the objects needed to build the Hansen coefficients. <br/>
* The index is s*/
private HansenZonalLinear[] hansenObjects;
/** The current value of the U function. <br/>
* Needed when computed the short periodic contribution */
private double U;
/** A = sqrt( μ * a ) = n * a². */
private double A;
/** B = sqrt( 1 - k² - h² ). */
private double B;
/** C = 1 + p² + Q². */
private double C;
/** The mean motion (n). */
private double meanMotion;
/** h * k. */
private double hk;
/** k² - h². */
private double k2mh2;
/** (k² - h²) / 2. */
private double k2mh2o2;
/** 1 / (n² * a²). */
private double oon2a2;
/** 1 / (n² * a) . */
private double oon2a;
/** χ³ / (n² * a). */
private double x3on2a;
/** χ / (n² * a²). */
private double xon2a2;
/** (C * χ) / ( 2 * n² * a² ). */
private double cxo2n2a2;
/** (χ²) / (n² * a² * (χ + 1 ) ). */
private double x2on2a2xp1;
/** B * B.*/
private double BB;
/** Simple constructor.
* @param provider provider for spherical harmonics
* @param maxDegreeShortPeriodics maximum degree to consider for short periodics zonal harmonics potential
* (must be between 2 and {@code provider.getMaxDegree()})
* @param maxEccPowShortPeriodics maximum power of the eccentricity to be used in short periodic computations
* (must be between 0 and {@code maxDegreeShortPeriodics - 1}, but should typically not exceed 4 as higher
* values will exceed computer capacity)
* @param maxFrequencyShortPeriodics maximum frequency in true longitude for short periodic computations
* (must be between 1 and {@code 2 * maxDegreeShortPeriodics + 1})
* @exception OrekitException if degrees or powers are out of range
* @since 7.2
*/
public DSSTZonal(final UnnormalizedSphericalHarmonicsProvider provider,
final int maxDegreeShortPeriodics,
final int maxEccPowShortPeriodics,
final int maxFrequencyShortPeriodics)
throws OrekitException {
this.provider = provider;
this.maxDegree = provider.getMaxDegree();
this.maxOrder = provider.getMaxOrder();
checkIndexRange(maxDegreeShortPeriodics, 2, maxDegree);
this.maxDegreeShortPeriodics = maxDegreeShortPeriodics;
checkIndexRange(maxEccPowShortPeriodics, 0, maxDegreeShortPeriodics - 1);
this.maxEccPowShortPeriodics = maxEccPowShortPeriodics;
checkIndexRange(maxFrequencyShortPeriodics, 1, 2 * maxDegreeShortPeriodics + 1);
this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
// Vns coefficients
this.Vns = CoefficientsFactory.computeVns(maxDegree + 1);
// Factorials computation
final int maxFact = 2 * maxDegree + 1;
this.fact = new double[maxFact];
fact[0] = 1.;
for (int i = 1; i < maxFact; i++) {
fact[i] = i * fact[i - 1];
}
// Initialize default values
this.maxEccPowMeanElements = (maxDegree == 2) ? 0 : Integer.MIN_VALUE;
}
/** Check an index range.
* @param index index value
* @param min minimum value for index
* @param max maximum value for index
* @exception OrekitException if index is out of range
*/
private void checkIndexRange(final int index, final int min, final int max)
throws OrekitException {
if (index < min) {
throw new OrekitException(LocalizedFormats.NUMBER_TOO_SMALL, index, min);
}
if (index > max) {
throw new OrekitException(LocalizedFormats.NUMBER_TOO_LARGE, index, max);
}
}
/** Get the spherical harmonics provider.
* @return the spherical harmonics provider
*/
public UnnormalizedSphericalHarmonicsProvider getProvider() {
return provider;
}
/** {@inheritDoc}
* <p>
* Computes the highest power of the eccentricity to appear in the truncated
* analytical power series expansion.
* </p>
* <p>
* This method computes the upper value for the central body potential and
* determines the maximal power for the eccentricity producing potential
* terms bigger than a defined tolerance.
* </p>
*/
@Override
public List<ShortPeriodTerms> initialize(final AuxiliaryElements aux, final boolean meanOnly)
throws OrekitException {
computeMeanElementsTruncations(aux);
final int maxEccPow;
if (!meanOnly) {
maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
} else {
maxEccPow = maxEccPowMeanElements;
}
//Initialize the HansenCoefficient generator
this.hansenObjects = new HansenZonalLinear[maxEccPow + 1];
for (int s = 0; s <= maxEccPow; s++) {
this.hansenObjects[s] = new HansenZonalLinear(maxDegree, s);
}
final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
zonalSPCoefs = new ZonalShortPeriodicCoefficients(maxDegreeShortPeriodics, maxFrequencyShortPeriodics,
INTERPOLATION_POINTS,
new TimeSpanMap<Slot>(new Slot(maxFrequencyShortPeriodics,
INTERPOLATION_POINTS)));
list.add(zonalSPCoefs);
return list;
}
/** Compute indices truncations for mean elements computations.
* @param aux auxiliary elements
* @throws OrekitException if an error occurs
*/
private void computeMeanElementsTruncations(final AuxiliaryElements aux) throws OrekitException {
//Compute the max eccentricity power for the mean element rate expansion
if (maxDegree == 2) {
maxEccPowMeanElements = 0;
} else {
// Initializes specific parameters.
initializeStep(aux);
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(aux.getDate());
// Utilities for truncation
final double ax2or = 2. * a / provider.getAe();
double xmuran = provider.getMu() / a;
// 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[maxDegree + 1];
final double[] chiPwr = new double[maxDegree + 1];
final double[] hafPwr = new double[maxDegree + 1];
eccPwr[0] = 1.;
chiPwr[0] = X;
hafPwr[0] = 1.;
for (int i = 1; i <= maxDegree; i++) {
eccPwr[i] = eccPwr[i - 1] * eo2;
chiPwr[i] = chiPwr[i - 1] * x2o2;
hafPwr[i] = hafPwr[i - 1] * 0.5;
xmuran /= ax2or;
}
// Set highest power of e and degree of current spherical harmonic.
maxEccPowMeanElements = 0;
int n = maxDegree;
// Loop over n
do {
// Set order of current spherical harmonic.
int m = 0;
// Loop over m
do {
// Compute magnitude of current spherical harmonic coefficient.
final double cnm = harmonics.getUnnormalizedCnm(n, m);
final double snm = harmonics.getUnnormalizedSnm(n, m);
final double csnm = FastMath.hypot(cnm, snm);
if (csnm == 0.) break;
// Set magnitude of last spherical harmonic term.
double lastTerm = 0.;
// Set current power of e and related indices.
int nsld2 = (n - maxEccPowMeanElements - 1) / 2;
int l = n - 2 * nsld2;
// Loop over l
double term = 0.;
do {
// Compute magnitude of current spherical harmonic term.
if (m < l) {
term = csnm * xmuran *
(fact[n - l] / (fact[n - m])) *
(fact[n + l] / (fact[nsld2] * fact[nsld2 + l])) *
eccPwr[l] * UpperBounds.getDnl(XX, chiPwr[l], n, l) *
(UpperBounds.getRnml(gamma, n, l, m, 1, I) + UpperBounds.getRnml(gamma, n, l, m, -1, I));
} else {
term = csnm * xmuran *
(fact[n + m] / (fact[nsld2] * fact[nsld2 + l])) *
eccPwr[l] * hafPwr[m - l] * UpperBounds.getDnl(XX, chiPwr[l], n, l) *
(UpperBounds.getRnml(gamma, n, m, l, 1, I) + UpperBounds.getRnml(gamma, n, m, l, -1, I));
}
// Is the current spherical harmonic term bigger than the truncation tolerance ?
if (term >= TRUNCATION_TOLERANCE) {
maxEccPowMeanElements = l;
} else {
// Is the current term smaller than the last term ?
if (term < lastTerm) {
break;
}
}
// Proceed to next power of e.
lastTerm = term;
l += 2;
nsld2--;
} while (l < n);
// Is the current spherical harmonic term bigger than the truncation tolerance ?
if (term >= TRUNCATION_TOLERANCE) {
maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
return;
}
// Proceed to next order.
m++;
} while (m <= FastMath.min(n, maxOrder));
// Proceed to next degree.
xmuran *= ax2or;
n--;
} while (n > maxEccPowMeanElements + 2);
maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
}
}
/** {@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();
// Direction cosines
alpha = aux.getAlpha();
beta = aux.getBeta();
gamma = aux.getGamma();
// Equinoctial coefficients
A = aux.getA();
B = aux.getB();
C = aux.getC();
// Χ = 1 / B
X = 1. / B;
XX = X * X;
XXX = X * XX;
// 1 / AB
ooAB = 1. / (A * B);
// B / A
BoA = B / A;
// -C / 2AB
mCo2AB = -C * ooAB / 2.;
// B / A(1 + B)
BoABpo = BoA / (1. + B);
// -2 * a / A
m2aoA = -2 * a / A;
// μ / a
muoa = provider.getMu() / a;
// R / a
roa = provider.getAe() / a;
// Mean motion
meanMotion = aux.getMeanMotion();
}
/** {@inheritDoc} */
@Override
public double[] getMeanElementRate(final SpacecraftState spacecraftState) throws OrekitException {
return computeMeanElementRates(spacecraftState.getDate());
}
/** {@inheritDoc} */
@Override
public EventDetector[] getEventsDetectors() {
return null;
}
/** Compute the mean element rates.
* @param date current date
* @return the mean element rates
* @throws OrekitException if an error occurs in hansen computation
*/
private double[] computeMeanElementRates(final AbsoluteDate date) throws OrekitException {
// Compute potential derivative
final double[] dU = computeUDerivatives(date);
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};
}
/** Compute the derivatives of the gravitational potential U [Eq. 3.1-(6)].
* <p>
* The result is the array
* [dU/da, dU/dk, dU/dh, dU/dα, dU/dβ, dU/dγ]
* </p>
* @param date current date
* @return potential derivatives
* @throws OrekitException if an error occurs in hansen computation
*/
private double[] computeUDerivatives(final AbsoluteDate date) throws OrekitException {
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
//Reset U
U = 0.;
// Gs and Hs coefficients
final double[][] GsHs = CoefficientsFactory.computeGsHs(k, h, alpha, beta, maxEccPowMeanElements);
// Qns coefficients
final double[][] Qns = CoefficientsFactory.computeQns(gamma, maxDegree, maxEccPowMeanElements);
final double[] roaPow = new double[maxDegree + 1];
roaPow[0] = 1.;
for (int i = 1; i <= maxDegree; i++) {
roaPow[i] = roa * roaPow[i - 1];
}
// 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 <= maxEccPowMeanElements; s++) {
//Initialize the Hansen roots
this.hansenObjects[s].computeInitValues(X);
// 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 d0s = (s == 0) ? 1 : 2;
for (int n = s + 2; n <= maxDegree; n++) {
// (n - s) must be even
if ((n - s) % 2 == 0) {
//Extract data from previous computation :
final double kns = this.hansenObjects[s].getValue(-n - 1, X);
final double dkns = this.hansenObjects[s].getDerivative(-n - 1, X);
final double vns = Vns.get(new NSKey(n, s));
final double coef0 = d0s * roaPow[n] * vns * -harmonics.getUnnormalizedCnm(n, 0);
final double coef1 = coef0 * Qns[n][s];
final double coef2 = coef1 * kns;
final double coef3 = coef2 * gs;
// dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
final double dqns = Qns[n][s + 1];
// Compute U
U += coef3;
// Compute dU / da :
dUda += coef3 * (n + 1);
// Compute dU / dEx
dUdk += coef1 * (kns * dGsdk + k * XXX * gs * dkns);
// Compute dU / dEy
dUdh += coef1 * (kns * dGsdh + h * XXX * gs * dkns);
// Compute dU / dAlpha
dUdAl += coef2 * dGsdAl;
// Compute dU / dBeta
dUdBe += coef2 * dGsdBe;
// Compute dU / dGamma
dUdGa += coef0 * kns * dqns * gs;
}
}
}
// Multiply by -(μ / a)
U *= -muoa;
return new double[] {
dUda * muoa / a,
dUdk * -muoa,
dUdh * -muoa,
dUdAl * -muoa,
dUdBe * -muoa,
dUdGa * -muoa
};
}
/** {@inheritDoc} */
@Override
public void registerAttitudeProvider(final AttitudeProvider attitudeProvider) {
//nothing is done since this contribution is not sensitive to attitude
}
/** Check if an index is within the accepted interval.
*
* @param index the index to check
* @param lowerBound the lower bound of the interval
* @param upperBound the upper bound of the interval
* @return true if the index is between the lower and upper bounds, false otherwise
*/
private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
return index >= lowerBound && index <= upperBound;
}
/** {@inheritDoc} */
@Override
public void updateShortPeriodTerms(final SpacecraftState ... meanStates)
throws OrekitException {
final Slot slot = zonalSPCoefs.createSlot(meanStates);
for (final SpacecraftState meanState : meanStates) {
initializeStep(new AuxiliaryElements(meanState.getOrbit(), I));
// h * k.
this.hk = h * k;
// k² - h².
this.k2mh2 = k * k - h * h;
// (k² - h²) / 2.
this.k2mh2o2 = k2mh2 / 2.;
// 1 / (n² * a²) = 1 / (n * A)
this.oon2a2 = 1 / (A * meanMotion);
// 1 / (n² * a) = a / (n * A)
this.oon2a = a * oon2a2;
// χ³ / (n² * a)
this.x3on2a = XXX * oon2a;
// χ / (n² * a²)
this.xon2a2 = X * oon2a2;
// (C * χ) / ( 2 * n² * a² )
this.cxo2n2a2 = xon2a2 * C / 2;
// (χ²) / (n² * a² * (χ + 1 ) )
this.x2on2a2xp1 = xon2a2 * X / (X + 1);
// B * B
this.BB = B * B;
// Compute rhoj and sigmaj
final double[][] rhoSigma = computeRhoSigmaCoefficients(meanState.getDate(), slot);
// Compute Di
computeDiCoefficients(meanState.getDate(), slot);
// generate the Cij and Sij coefficients
final FourierCjSjCoefficients cjsj = new FourierCjSjCoefficients(meanState.getDate(),
maxDegreeShortPeriodics,
maxEccPowShortPeriodics,
maxFrequencyShortPeriodics);
computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma);
}
}
/** Generate the values for the D<sub>i</sub> coefficients.
* @param date target date
* @param slot slot to which the coefficients belong
* @throws OrekitException if an error occurs during the coefficient computation
*/
private void computeDiCoefficients(final AbsoluteDate date, final Slot slot)
throws OrekitException {
final double[] meanElementRates = computeMeanElementRates(date);
final double[] currentDi = new double[6];
// Add the coefficients to the interpolation grid
for (int i = 0; i < 6; i++) {
currentDi[i] = meanElementRates[i] / meanMotion;
if (i == 5) {
currentDi[i] += -1.5 * 2 * U * oon2a2;
}
}
slot.di.addGridPoint(date, currentDi);
}
/**
* Generate the values for the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
* @param date date of computation
* @param slot slot to which the coefficients belong
* @param cjsj Fourier coefficients
* @param rhoSigma ρ<sub>j</sub> and σ<sub>j</sub>
*/
private void computeCijSijCoefficients(final AbsoluteDate date, final Slot slot,
final FourierCjSjCoefficients cjsj,
final double[][] rhoSigma) {
final int nMax = maxDegreeShortPeriodics;
// The C<sub>i</sub>⁰ coefficients
final double[] currentCi0 = new double[] {0., 0., 0., 0., 0., 0.};
for (int j = 1; j < slot.cij.length; j++) {
// Create local arrays
final double[] currentCij = new double[] {0., 0., 0., 0., 0., 0.};
final double[] currentSij = new double[] {0., 0., 0., 0., 0., 0.};
// j == 1
if (j == 1) {
final double coef1 = 4 * k * U - hk * cjsj.getCj(1) + k2mh2o2 * cjsj.getSj(1);
final double coef2 = 4 * h * U + k2mh2o2 * cjsj.getCj(1) + hk * cjsj.getSj(1);
final double coef3 = (k * cjsj.getCj(1) + h * cjsj.getSj(1)) / 4.;
final double coef4 = (8 * U - h * cjsj.getCj(1) + k * cjsj.getSj(1)) / 4.;
//Coefficients for a
currentCij[0] += coef1;
currentSij[0] += coef2;
//Coefficients for k
currentCij[1] += coef4;
currentSij[1] += coef3;
//Coefficients for h
currentCij[2] -= coef3;
currentSij[2] += coef4;
//Coefficients for λ
currentCij[5] -= coef2 / 2;
currentSij[5] += coef1 / 2;
}
// j == 2
if (j == 2) {
final double coef1 = k2mh2 * U;
final double coef2 = 2 * hk * U;
final double coef3 = h * U / 2;
final double coef4 = k * U / 2;
//Coefficients for a
currentCij[0] += coef1;
currentSij[0] += coef2;
//Coefficients for k
currentCij[1] += coef4;
currentSij[1] += coef3;
//Coefficients for h
currentCij[2] -= coef3;
currentSij[2] += coef4;
//Coefficients for λ
currentCij[5] -= coef2 / 2;
currentSij[5] += coef1 / 2;
}
// j between 1 and 2N-3
if (isBetween(j, 1, 2 * nMax - 3)) {
final double coef1 = ( j + 2 ) * (-hk * cjsj.getCj(j + 2) + k2mh2o2 * cjsj.getSj(j + 2));
final double coef2 = ( j + 2 ) * (k2mh2o2 * cjsj.getCj(j + 2) + hk * cjsj.getSj(j + 2));
final double coef3 = ( j + 2 ) * (k * cjsj.getCj(j + 2) + h * cjsj.getSj(j + 2)) / 4;
final double coef4 = ( j + 2 ) * (h * cjsj.getCj(j + 2) - k * cjsj.getSj(j + 2)) / 4;
//Coefficients for a
currentCij[0] += coef1;
currentSij[0] -= coef2;
//Coefficients for k
currentCij[1] += -coef4;
currentSij[1] -= coef3;
//Coefficients for h
currentCij[2] -= coef3;
currentSij[2] += coef4;
//Coefficients for λ
currentCij[5] -= coef2 / 2;
currentSij[5] += -coef1 / 2;
}
// j between 1 and 2N-2
if (isBetween(j, 1, 2 * nMax - 2)) {
final double coef1 = 2 * ( j + 1 ) * (-h * cjsj.getCj(j + 1) + k * cjsj.getSj(j + 1));
final double coef2 = 2 * ( j + 1 ) * (k * cjsj.getCj(j + 1) + h * cjsj.getSj(j + 1));
final double coef3 = ( j + 1 ) * cjsj.getCj(j + 1);
final double coef4 = ( j + 1 ) * cjsj.getSj(j + 1);
//Coefficients for a
currentCij[0] += coef1;
currentSij[0] -= coef2;
//Coefficients for k
currentCij[1] += coef4;
currentSij[1] -= coef3;
//Coefficients for h
currentCij[2] -= coef3;
currentSij[2] -= coef4;
//Coefficients for λ
currentCij[5] -= coef2 / 2;
currentSij[5] += -coef1 / 2;
}
// j between 2 and 2N
if (isBetween(j, 2, 2 * nMax)) {
final double coef1 = 2 * ( j - 1 ) * (h * cjsj.getCj(j - 1) + k * cjsj.getSj(j - 1));
final double coef2 = 2 * ( j - 1 ) * (k * cjsj.getCj(j - 1) - h * cjsj.getSj(j - 1));
final double coef3 = ( j - 1 ) * cjsj.getCj(j - 1);
final double coef4 = ( j - 1 ) * cjsj.getSj(j - 1);
//Coefficients for a
currentCij[0] += coef1;
currentSij[0] -= coef2;
//Coefficients for k
currentCij[1] += coef4;
currentSij[1] -= coef3;
//Coefficients for h
currentCij[2] += coef3;
currentSij[2] += coef4;
//Coefficients for λ
currentCij[5] += coef2 / 2;
currentSij[5] += coef1 / 2;
}
// j between 3 and 2N + 1
if (isBetween(j, 3, 2 * nMax + 1)) {
final double coef1 = ( j - 2 ) * (hk * cjsj.getCj(j - 2) + k2mh2o2 * cjsj.getSj(j - 2));
final double coef2 = ( j - 2 ) * (-k2mh2o2 * cjsj.getCj(j - 2) + hk * cjsj.getSj(j - 2));
final double coef3 = ( j - 2 ) * (k * cjsj.getCj(j - 2) - h * cjsj.getSj(j - 2)) / 4;
final double coef4 = ( j - 2 ) * (h * cjsj.getCj(j - 2) + k * cjsj.getSj(j - 2)) / 4;
final double coef5 = ( j - 2 ) * (k2mh2o2 * cjsj.getCj(j - 2) - hk * cjsj.getSj(j - 2));
//Coefficients for a
currentCij[0] += coef1;
currentSij[0] += coef2;
//Coefficients for k
currentCij[1] += coef4;
currentSij[1] += -coef3;
//Coefficients for h
currentCij[2] += coef3;
currentSij[2] += coef4;
//Coefficients for λ
currentCij[5] += coef5 / 2;
currentSij[5] += coef1 / 2;
}
//multiply by the common factor
//for a (i == 0) -> χ³ / (n² * a)
currentCij[0] *= this.x3on2a;
currentSij[0] *= this.x3on2a;
//for k (i == 1) -> χ / (n² * a²)
currentCij[1] *= this.xon2a2;
currentSij[1] *= this.xon2a2;
//for h (i == 2) -> χ / (n² * a²)
currentCij[2] *= this.xon2a2;
currentSij[2] *= this.xon2a2;
//for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
currentCij[5] *= this.x2on2a2xp1;
currentSij[5] *= this.x2on2a2xp1;
// j is between 1 and 2 * N - 1
if (isBetween(j, 1, 2 * nMax - 1)) {
// Compute cross derivatives
// Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
final double CjAlphaGamma = alpha * cjsj.getdCjdGamma(j) - gamma * cjsj.getdCjdAlpha(j);
// Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
final double CjAlphaBeta = alpha * cjsj.getdCjdBeta(j) - beta * cjsj.getdCjdAlpha(j);
// Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
final double CjBetaGamma = beta * cjsj.getdCjdGamma(j) - gamma * cjsj.getdCjdBeta(j);
// Cj(h,k) = h * dC/dk - k * dC/dh
final double CjHK = h * cjsj.getdCjdK(j) - k * cjsj.getdCjdH(j);
// Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
final double SjAlphaGamma = alpha * cjsj.getdSjdGamma(j) - gamma * cjsj.getdSjdAlpha(j);
// Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
final double SjAlphaBeta = alpha * cjsj.getdSjdBeta(j) - beta * cjsj.getdSjdAlpha(j);
// Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
final double SjBetaGamma = beta * cjsj.getdSjdGamma(j) - gamma * cjsj.getdSjdBeta(j);
// Sj(h,k) = h * dS/dk - k * dS/dh
final double SjHK = h * cjsj.getdSjdK(j) - k * cjsj.getdSjdH(j);
//Coefficients for a
final double coef1 = this.x3on2a * (3 - BB) * j;
currentCij[0] += coef1 * cjsj.getSj(j);
currentSij[0] -= coef1 * cjsj.getCj(j);
//Coefficients for k and h
final double coef2 = p * CjAlphaGamma - I * q * CjBetaGamma;
final double coef3 = p * SjAlphaGamma - I * q * SjBetaGamma;
currentCij[1] -= this.xon2a2 * (h * coef2 + BB * cjsj.getdCjdH(j) - 1.5 * k * j * cjsj.getSj(j));
currentSij[1] -= this.xon2a2 * (h * coef3 + BB * cjsj.getdSjdH(j) + 1.5 * k * j * cjsj.getCj(j));
currentCij[2] += this.xon2a2 * (k * coef2 + BB * cjsj.getdCjdK(j) + 1.5 * h * j * cjsj.getSj(j));
currentSij[2] += this.xon2a2 * (k * coef3 + BB * cjsj.getdSjdK(j) - 1.5 * h * j * cjsj.getCj(j));
//Coefficients for q and p
final double coef4 = CjHK - CjAlphaBeta - j * cjsj.getSj(j);
final double coef5 = SjHK - SjAlphaBeta + j * cjsj.getCj(j);
currentCij[3] = this.cxo2n2a2 * (-I * CjAlphaGamma + q * coef4);
currentSij[3] = this.cxo2n2a2 * (-I * SjAlphaGamma + q * coef5);
currentCij[4] = this.cxo2n2a2 * (-CjBetaGamma + p * coef4);
currentSij[4] = this.cxo2n2a2 * (-SjBetaGamma + p * coef5);
//Coefficients for λ
final double coef6 = h * cjsj.getdCjdH(j) + k * cjsj.getdCjdK(j);
final double coef7 = h * cjsj.getdSjdH(j) + k * cjsj.getdSjdK(j);
currentCij[5] += this.oon2a2 * (-2 * a * cjsj.getdCjdA(j) + coef6 / (X + 1) + X * coef2 - 3 * cjsj.getCj(j));
currentSij[5] += this.oon2a2 * (-2 * a * cjsj.getdSjdA(j) + coef7 / (X + 1) + X * coef3 - 3 * cjsj.getSj(j));
}
for (int i = 0; i < 6; i++) {
//Add the current coefficients contribution to C<sub>i</sub>⁰
currentCi0[i] -= currentCij[i] * rhoSigma[j][0] + currentSij[i] * rhoSigma[j][1];
}
// Add the coefficients to the interpolation grid
slot.cij[j].addGridPoint(date, currentCij);
slot.sij[j].addGridPoint(date, currentSij);
}
//Add C<sub>i</sub>⁰ to the interpolation grid
slot.cij[0].addGridPoint(date, currentCi0);
}
/**
* Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
* <p>
* The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
* ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
* σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
* </p>
* @param date target date
* @param slot slot to which the coefficients belong
* @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
*/
private double[][] computeRhoSigmaCoefficients(final AbsoluteDate date, final Slot slot) {
final CjSjCoefficient cjsjKH = new CjSjCoefficient(k, h);
final double b = 1. / (1 + B);
// (-b)<sup>j</sup>
double mbtj = 1;
final double[][] rhoSigma = new double[slot.cij.length][2];
for (int j = 1; j < rhoSigma.length; j++) {
//Compute current rho and sigma;
mbtj *= -b;
final double coef = (1 + j * B) * mbtj;
final double rho = coef * cjsjKH.getCj(j);
final double sigma = coef * cjsjKH.getSj(j);
// Add the coefficients to the interpolation grid
rhoSigma[j][0] = rho;
rhoSigma[j][1] = sigma;
}
return rhoSigma;
}
/** The coefficients used to compute the short-periodic zonal contribution.
*
* <p>
* Those coefficients are given in Danielson paper by expressions 4.1-(20) to 4.1.-(25)
* </p>
* <p>
* The coefficients are: <br/>
* - C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> <br/>
* - ρ<sub>j</sub> and σ<sub>j</sub> <br/>
* - C<sub>i</sub>⁰
* </p>
*
* @author Lucian Barbulescu
*/
private static class ZonalShortPeriodicCoefficients implements ShortPeriodTerms {
/** Serializable UID. */
private static final long serialVersionUID = 20151118L;
/** Maximal degree to consider for harmonics potential. */
private final int maxDegreeShortPeriodics;
/** Maximum value for j index. */
private final int maxFrequencyShortPeriodics;
/** Number of points used in the interpolation process. */
private final int interpolationPoints;
/** All coefficients slots. */
private final transient TimeSpanMap<Slot> slots;
/** Constructor.
* @param maxDegreeShortPeriodics maximal degree to consider for harmonics potential
* @param maxFrequencyShortPeriodics maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
* @param slots all coefficients slots
*/
ZonalShortPeriodicCoefficients(final int maxDegreeShortPeriodics,
final int maxFrequencyShortPeriodics, final int interpolationPoints,
final TimeSpanMap<Slot> slots) {
// Save parameters
this.maxDegreeShortPeriodics = maxDegreeShortPeriodics;
this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
this.interpolationPoints = interpolationPoints;
this.slots = slots;
}
/** 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(maxFrequencyShortPeriodics, 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());
// Get the True longitude L
final double L = meanOrbit.getLv();
// Define maxJ
final int maxJ = 2 * maxDegreeShortPeriodics + 1;
// Compute the center
final double center = L - meanOrbit.getLM();
// Initialize short periodic variations
final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
final double[] d = slot.di.value(meanOrbit.getDate());
for (int i = 0; i < 6; i++) {
shortPeriodicVariation[i] += center * d[i];
}
for (int j = 1; j <= maxJ; j++) {
final double[] c = slot.cij[j].value(meanOrbit.getDate());
final double[] s = slot.sij[j].value(meanOrbit.getDate());
final double cos = FastMath.cos(j * L);
final double sin = FastMath.sin(j * L);
for (int i = 0; i < 6; i++) {
// add corresponding term to the short periodic variation
shortPeriodicVariation[i] += c[i] * cos;
shortPeriodicVariation[i] += s[i] * sin;
}
}
return shortPeriodicVariation;
}
/** {@inheritDoc} */
@Override
public String getCoefficientsKeyPrefix() {
return "DSST-central-body-zonal-";
}
/** {@inheritDoc}
* <p>
* For zonal terms contributions,there are maxJ cj coefficients,
* maxJ sj coefficients and 2 dj coefficients, where maxJ depends
* on the orbit. The j index is the integer multiplier for the true
* longitude argument in the cj and sj coefficients and the degree
* in the polynomial dj 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 int maxJ = 2 * maxDegreeShortPeriodics + 1;
final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * maxJ + 2);
storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
storeIfSelected(coefficients, selected, slot.di.value(date), "d", 1);
for (int j = 1; j <= maxJ; 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);
}
}
/** Replace the instance with a data transfer object for serialization.
* @return data transfer object that will be serialized
* @exception NotSerializableException if an additional state provider is not serializable
*/
private Object writeReplace() throws NotSerializableException {
// slots transitions
final SortedSet<TimeSpanMap.Transition<Slot>> transitions = slots.getTransitions();
final AbsoluteDate[] transitionDates = new AbsoluteDate[transitions.size()];
final Slot[] allSlots = new Slot[transitions.size() + 1];
int i = 0;
for (final TimeSpanMap.Transition<Slot> transition : transitions) {
if (i == 0) {
// slot before the first transition
allSlots[i] = transition.getBefore();
}
if (i < transitionDates.length) {
transitionDates[i] = transition.getDate();
allSlots[++i] = transition.getAfter();
}
}
return new DataTransferObject(maxDegreeShortPeriodics,
maxFrequencyShortPeriodics, interpolationPoints,
transitionDates, allSlots);
}
/** Internal class used only for serialization. */
private static class DataTransferObject implements Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20160319L;
/** Maximal degree to consider for harmonics potential. */
private final int maxDegreeShortPeriodics;
/** Maximum value for j index. */
private final int maxFrequencyShortPeriodics;
/** Number of points used in the interpolation process. */
private final int interpolationPoints;
/** Transitions dates. */
private final AbsoluteDate[] transitionDates;
/** All slots. */
private final Slot[] allSlots;
/** Simple constructor.
* @param maxDegreeShortPeriodics maximal degree to consider for harmonics potential
* @param maxFrequencyShortPeriodics maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
* @param transitionDates transitions dates
* @param allSlots all slots
*/
DataTransferObject(final int maxDegreeShortPeriodics,
final int maxFrequencyShortPeriodics, final int interpolationPoints,
final AbsoluteDate[] transitionDates, final Slot[] allSlots) {
this.maxDegreeShortPeriodics = maxDegreeShortPeriodics;
this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
this.interpolationPoints = interpolationPoints;
this.transitionDates = transitionDates;
this.allSlots = allSlots;
}
/** Replace the deserialized data transfer object with a {@link ZonalShortPeriodicCoefficients}.
* @return replacement {@link ZonalShortPeriodicCoefficients}
*/
private Object readResolve() {
final TimeSpanMap<Slot> slots = new TimeSpanMap<Slot>(allSlots[0]);
for (int i = 0; i < transitionDates.length; ++i) {
slots.addValidAfter(allSlots[i + 1], transitionDates[i]);
}
return new ZonalShortPeriodicCoefficients(maxDegreeShortPeriodics,
maxFrequencyShortPeriodics,
interpolationPoints,
slots);
}
}
}
/** Compute the C<sup>j</sup> and the S<sup>j</sup> coefficients.
* <p>
* Those coefficients are given in Danielson paper by expressions 4.1-(13) to 4.1.-(16b)
* </p>
*/
private class FourierCjSjCoefficients {
/** The G<sub>js</sub>, H<sub>js</sub>, I<sub>js</sub> and J<sub>js</sub> polynomials. */
private final GHIJjsPolynomials ghijCoef;
/** L<sub>n</sub><sup>s</sup>(γ). */
private final LnsCoefficients lnsCoef;
/** Maximum possible value for n. */
private final int nMax;
/** Maximum possible value for s. */
private final int sMax;
/** Maximum possible value for j. */
private final int jMax;
/** The C<sup>j</sup> coefficients 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> / dh <br/>
* - dC<sup>j</sup> / dk <br/>
* - dC<sup>j</sup> / dα <br/>
* - dC<sup>j</sup> / dβ <br/>
* - dC<sup>j</sup> / dγ <br/>
* </p>
*/
private final double[][] cCoef;
/** 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> / dh <br/>
* - dS<sup>j</sup> / dk <br/>
* - dS<sup>j</sup> / dα <br/>
* - dS<sup>j</sup> / dβ <br/>
* - dS<sup>j</sup> / dγ <br/>
* </p>
*/
private final double[][] sCoef;
/** h * Χ³. */
private final double hXXX;
/** k * Χ³. */
private final double kXXX;
/** Create a set of C<sup>j</sup> and the S<sup>j</sup> coefficients.
* @param date the current date
* @param nMax maximum possible value for n
* @param sMax maximum possible value for s
* @param jMax maximum possible value for j
* @throws OrekitException if an error occurs while generating the coefficients
*/
FourierCjSjCoefficients(final AbsoluteDate date,
final int nMax, final int sMax, final int jMax)
throws OrekitException {
this.ghijCoef = new GHIJjsPolynomials(k, h, alpha, beta);
// Qns coefficients
final double[][] Qns = CoefficientsFactory.computeQns(gamma, nMax, nMax);
this.lnsCoef = new LnsCoefficients(nMax, nMax, Qns, Vns, roa);
this.nMax = nMax;
this.sMax = sMax;
this.jMax = jMax;
// compute the common factors that depends on the mean elements
this.hXXX = h * XXX;
this.kXXX = k * XXX;
this.cCoef = new double[7][jMax + 1];
this.sCoef = new double[7][jMax + 1];
for (int s = 0; s <= sMax; s++) {
//Initialise the Hansen roots
hansenObjects[s].computeInitValues(X);
}
generateCoefficients(date);
}
/** Generate all coefficients.
* @param date the current date
* @throws OrekitException if an error occurs while generating the coefficients
*/
private void generateCoefficients(final AbsoluteDate date) throws OrekitException {
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
for (int j = 1; j <= jMax; j++) {
//init arrays
for (int i = 0; i <= 6; i++) {
cCoef[i][j] = 0.;
sCoef[i][j] = 0.;
}
if (isBetween(j, 1, nMax - 1)) {
//compute first double sum where s: j -> N-1 and n: s+1 -> N
for (int s = j; s <= FastMath.min(nMax - 1, sMax); s++) {
// j - s
final int jms = j - s;
// Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
final int d0smj = (s == j) ? 1 : 2;
for (int n = s + 1; n <= nMax; n++) {
// if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
if ((n + jms) % 2 == 0) {
// (2 - delta(0,s-j)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
final double lns = lnsCoef.getLns(n, -jms);
final double dlns = lnsCoef.getdLnsdGamma(n, -jms);
final double hjs = ghijCoef.getHjs(s, -jms);
final double dHjsdh = ghijCoef.getdHjsdh(s, -jms);
final double dHjsdk = ghijCoef.getdHjsdk(s, -jms);
final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, -jms);
final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, -jms);
final double gjs = ghijCoef.getGjs(s, -jms);
final double dGjsdh = ghijCoef.getdGjsdh(s, -jms);
final double dGjsdk = ghijCoef.getdGjsdk(s, -jms);
final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, -jms);
final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, -jms);
// J<sub>n</sub>
final double jn = -harmonics.getUnnormalizedCnm(n, 0);
// K₀<sup>-n-1,s</sup>
final double kns = hansenObjects[s].getValue(-n - 1, X);
final double dkns = hansenObjects[s].getDerivative(-n - 1, X);
final double coef0 = d0smj * jn;
final double coef1 = coef0 * lns;
final double coef2 = coef1 * kns;
final double coef3 = coef2 * hjs;
final double coef4 = coef2 * gjs;
// Add the term to the coefficients
cCoef[0][j] += coef3;
cCoef[1][j] += coef3 * (n + 1);
cCoef[2][j] += coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
cCoef[3][j] += coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
cCoef[4][j] += coef2 * dHjsdAlpha;
cCoef[5][j] += coef2 * dHjsdBeta;
cCoef[6][j] += coef0 * dlns * kns * hjs;
sCoef[0][j] += coef4;
sCoef[1][j] += coef4 * (n + 1);
sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
sCoef[4][j] += coef2 * dGjsdAlpha;
sCoef[5][j] += coef2 * dGjsdBeta;
sCoef[6][j] += coef0 * dlns * kns * gjs;
}
}
}
//compute second double sum where s: 0 -> N-j and n: max(j+s, j+1) -> N
for (int s = 0; s <= FastMath.min(nMax - j, sMax); s++) {
// j + s
final int jps = j + s;
// Kronecker symbols (2 - delta(0,j+s))
final double d0spj = (s == -j) ? 1 : 2;
for (int n = FastMath.max(j + s, j + 1); n <= nMax; n++) {
// if n + (j+s) is odd, then the term is equal to zero due to the factor Vn,s+j
if ((n + jps) % 2 == 0) {
// (2 - delta(0,s+j)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j+s</sup>
final double lns = lnsCoef.getLns(n, jps);
final double dlns = lnsCoef.getdLnsdGamma(n, jps);
final double hjs = ghijCoef.getHjs(s, jps);
final double dHjsdh = ghijCoef.getdHjsdh(s, jps);
final double dHjsdk = ghijCoef.getdHjsdk(s, jps);
final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, jps);
final double dHjsdBeta = ghijCoef.getdHjsdBeta(s, jps);
final double gjs = ghijCoef.getGjs(s, jps);
final double dGjsdh = ghijCoef.getdGjsdh(s, jps);
final double dGjsdk = ghijCoef.getdGjsdk(s, jps);
final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, jps);
final double dGjsdBeta = ghijCoef.getdGjsdBeta(s, jps);
// J<sub>n</sub>
final double jn = -harmonics.getUnnormalizedCnm(n, 0);
// K₀<sup>-n-1,s</sup>
final double kns = hansenObjects[s].getValue(-n - 1, X);
final double dkns = hansenObjects[s].getDerivative(-n - 1, X);
final double coef0 = d0spj * jn;
final double coef1 = coef0 * lns;
final double coef2 = coef1 * kns;
final double coef3 = coef2 * hjs;
final double coef4 = coef2 * gjs;
// Add the term to the coefficients
cCoef[0][j] -= coef3;
cCoef[1][j] -= coef3 * (n + 1);
cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
cCoef[4][j] -= coef2 * dHjsdAlpha;
cCoef[5][j] -= coef2 * dHjsdBeta;
cCoef[6][j] -= coef0 * dlns * kns * hjs;
sCoef[0][j] += coef4;
sCoef[1][j] += coef4 * (n + 1);
sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
sCoef[4][j] += coef2 * dGjsdAlpha;
sCoef[5][j] += coef2 * dGjsdBeta;
sCoef[6][j] += coef0 * dlns * kns * gjs;
}
}
}
//compute third double sum where s: 1 -> j and n: j+1 -> N
for (int s = 1; s <= FastMath.min(j, sMax); s++) {
// j - s
final int jms = j - s;
// Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
final int d0smj = (s == j) ? 1 : 2;
for (int n = j + 1; n <= nMax; n++) {
// if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
if ((n + jms) % 2 == 0) {
// (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
final double lns = lnsCoef.getLns(n, jms);
final double dlns = lnsCoef.getdLnsdGamma(n, jms);
final double ijs = ghijCoef.getIjs(s, jms);
final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final double jjs = ghijCoef.getJjs(s, jms);
final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final double jn = -harmonics.getUnnormalizedCnm(n, 0);
// K₀<sup>-n-1,s</sup>
final double kns = hansenObjects[s].getValue(-n - 1, X);
final double dkns = hansenObjects[s].getDerivative(-n - 1, X);
final double coef0 = d0smj * jn;
final double coef1 = coef0 * lns;
final double coef2 = coef1 * kns;
final double coef3 = coef2 * ijs;
final double coef4 = coef2 * jjs;
// Add the term to the coefficients
cCoef[0][j] -= coef3;
cCoef[1][j] -= coef3 * (n + 1);
cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
cCoef[4][j] -= coef2 * dIjsdAlpha;
cCoef[5][j] -= coef2 * dIjsdBeta;
cCoef[6][j] -= coef0 * dlns * kns * ijs;
sCoef[0][j] += coef4;
sCoef[1][j] += coef4 * (n + 1);
sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
sCoef[4][j] += coef2 * dJjsdAlpha;
sCoef[5][j] += coef2 * dJjsdBeta;
sCoef[6][j] += coef0 * dlns * kns * jjs;
}
}
}
}
if (isBetween(j, 2, nMax)) {
//add first term
// J<sub>j</sub>
final double jj = -harmonics.getUnnormalizedCnm(j, 0);
double kns = hansenObjects[0].getValue(-j - 1, X);
double dkns = hansenObjects[0].getDerivative(-j - 1, X);
double lns = lnsCoef.getLns(j, j);
//dlns is 0 because n == s == j
final double hjs = ghijCoef.getHjs(0, j);
final double dHjsdh = ghijCoef.getdHjsdh(0, j);
final double dHjsdk = ghijCoef.getdHjsdk(0, j);
final double dHjsdAlpha = ghijCoef.getdHjsdAlpha(0, j);
final double dHjsdBeta = ghijCoef.getdHjsdBeta(0, j);
final double gjs = ghijCoef.getGjs(0, j);
final double dGjsdh = ghijCoef.getdGjsdh(0, j);
final double dGjsdk = ghijCoef.getdGjsdk(0, j);
final double dGjsdAlpha = ghijCoef.getdGjsdAlpha(0, j);
final double dGjsdBeta = ghijCoef.getdGjsdBeta(0, j);
// 2 * J<sub>j</sub> * K₀<sup>-j-1,0</sup> * L<sub>j</sub><sup>j</sup>
double coef0 = 2 * jj;
double coef1 = coef0 * lns;
double coef2 = coef1 * kns;
double coef3 = coef2 * hjs;
double coef4 = coef2 * gjs;
// Add the term to the coefficients
cCoef[0][j] -= coef3;
cCoef[1][j] -= coef3 * (j + 1);
cCoef[2][j] -= coef1 * (kns * dHjsdh + hjs * hXXX * dkns);
cCoef[3][j] -= coef1 * (kns * dHjsdk + hjs * kXXX * dkns);
cCoef[4][j] -= coef2 * dHjsdAlpha;
cCoef[5][j] -= coef2 * dHjsdBeta;
//no contribution to cCoef[6][j] because dlns is 0
sCoef[0][j] += coef4;
sCoef[1][j] += coef4 * (j + 1);
sCoef[2][j] += coef1 * (kns * dGjsdh + gjs * hXXX * dkns);
sCoef[3][j] += coef1 * (kns * dGjsdk + gjs * kXXX * dkns);
sCoef[4][j] += coef2 * dGjsdAlpha;
sCoef[5][j] += coef2 * dGjsdBeta;
//no contribution to sCoef[6][j] because dlns is 0
//compute simple sum where s: 1 -> j-1
for (int s = 1; s <= FastMath.min(j - 1, sMax); s++) {
// j - s
final int jms = j - s;
// Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
final int d0smj = (s == j) ? 1 : 2;
// if s is odd, then the term is equal to zero due to the factor Vj,s-j
if (s % 2 == 0) {
// (2 - delta(0,j-s)) * J<sub>j</sub> * K₀<sup>-j-1,s</sup> * L<sub>j</sub><sup>j-s</sup>
kns = hansenObjects[s].getValue(-j - 1, X);
dkns = hansenObjects[s].getDerivative(-j - 1, X);
lns = lnsCoef.getLns(j, jms);
final double dlns = lnsCoef.getdLnsdGamma(j, jms);
final double ijs = ghijCoef.getIjs(s, jms);
final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final double jjs = ghijCoef.getJjs(s, jms);
final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
coef0 = d0smj * jj;
coef1 = coef0 * lns;
coef2 = coef1 * kns;
coef3 = coef2 * ijs;
coef4 = coef2 * jjs;
// Add the term to the coefficients
cCoef[0][j] -= coef3;
cCoef[1][j] -= coef3 * (j + 1);
cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
cCoef[4][j] -= coef2 * dIjsdAlpha;
cCoef[5][j] -= coef2 * dIjsdBeta;
cCoef[6][j] -= coef0 * dlns * kns * ijs;
sCoef[0][j] += coef4;
sCoef[1][j] += coef4 * (j + 1);
sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
sCoef[4][j] += coef2 * dJjsdAlpha;
sCoef[5][j] += coef2 * dJjsdBeta;
sCoef[6][j] += coef0 * dlns * kns * jjs;
}
}
}
if (isBetween(j, 3, 2 * nMax - 1)) {
//compute uppercase sigma expressions
//min(j-1,N)
final int minjm1on = FastMath.min(j - 1, nMax);
//if j is even
if (j % 2 == 0) {
//compute first double sum where s: j-min(j-1,N) -> j/2-1 and n: j-s -> min(j-1,N)
for (int s = j - minjm1on; s <= FastMath.min(j / 2 - 1, sMax); s++) {
// j - s
final int jms = j - s;
// Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
final int d0smj = (s == j) ? 1 : 2;
for (int n = j - s; n <= minjm1on; n++) {
// if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
if ((n + jms) % 2 == 0) {
// (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
final double lns = lnsCoef.getLns(n, jms);
final double dlns = lnsCoef.getdLnsdGamma(n, jms);
final double ijs = ghijCoef.getIjs(s, jms);
final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final double jjs = ghijCoef.getJjs(s, jms);
final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final double jn = -harmonics.getUnnormalizedCnm(n, 0);
// K₀<sup>-n-1,s</sup>
final double kns = hansenObjects[s].getValue(-n - 1, X);
final double dkns = hansenObjects[s].getDerivative(-n - 1, X);
final double coef0 = d0smj * jn;
final double coef1 = coef0 * lns;
final double coef2 = coef1 * kns;
final double coef3 = coef2 * ijs;
final double coef4 = coef2 * jjs;
// Add the term to the coefficients
cCoef[0][j] -= coef3;
cCoef[1][j] -= coef3 * (n + 1);
cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
cCoef[4][j] -= coef2 * dIjsdAlpha;
cCoef[5][j] -= coef2 * dIjsdBeta;
cCoef[6][j] -= coef0 * dlns * kns * ijs;
sCoef[0][j] += coef4;
sCoef[1][j] += coef4 * (n + 1);
sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
sCoef[4][j] += coef2 * dJjsdAlpha;
sCoef[5][j] += coef2 * dJjsdBeta;
sCoef[6][j] += coef0 * dlns * kns * jjs;
}
}
}
//compute second double sum where s: j/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
for (int s = j / 2; s <= FastMath.min(minjm1on - 1, sMax); s++) {
// j - s
final int jms = j - s;
// Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
final int d0smj = (s == j) ? 1 : 2;
for (int n = s + 1; n <= minjm1on; n++) {
// if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
if ((n + jms) % 2 == 0) {
// (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
final double lns = lnsCoef.getLns(n, jms);
final double dlns = lnsCoef.getdLnsdGamma(n, jms);
final double ijs = ghijCoef.getIjs(s, jms);
final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final double jjs = ghijCoef.getJjs(s, jms);
final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final double jn = -harmonics.getUnnormalizedCnm(n, 0);
// K₀<sup>-n-1,s</sup>
final double kns = hansenObjects[s].getValue(-n - 1, X);
final double dkns = hansenObjects[s].getDerivative(-n - 1, X);
final double coef0 = d0smj * jn;
final double coef1 = coef0 * lns;
final double coef2 = coef1 * kns;
final double coef3 = coef2 * ijs;
final double coef4 = coef2 * jjs;
// Add the term to the coefficients
cCoef[0][j] -= coef3;
cCoef[1][j] -= coef3 * (n + 1);
cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
cCoef[4][j] -= coef2 * dIjsdAlpha;
cCoef[5][j] -= coef2 * dIjsdBeta;
cCoef[6][j] -= coef0 * dlns * kns * ijs;
sCoef[0][j] += coef4;
sCoef[1][j] += coef4 * (n + 1);
sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
sCoef[4][j] += coef2 * dJjsdAlpha;
sCoef[5][j] += coef2 * dJjsdBeta;
sCoef[6][j] += coef0 * dlns * kns * jjs;
}
}
}
}
//if j is odd
else {
//compute first double sum where s: (j-1)/2 -> min(j-1,N)-1 and n: s+1 -> min(j-1,N)
for (int s = (j - 1) / 2; s <= FastMath.min(minjm1on - 1, sMax); s++) {
// j - s
final int jms = j - s;
// Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
final int d0smj = (s == j) ? 1 : 2;
for (int n = s + 1; n <= minjm1on; n++) {
// if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
if ((n + jms) % 2 == 0) {
// (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
final double lns = lnsCoef.getLns(n, jms);
final double dlns = lnsCoef.getdLnsdGamma(n, jms);
final double ijs = ghijCoef.getIjs(s, jms);
final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final double jjs = ghijCoef.getJjs(s, jms);
final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final double jn = -harmonics.getUnnormalizedCnm(n, 0);
// K₀<sup>-n-1,s</sup>
final double kns = hansenObjects[s].getValue(-n - 1, X);
final double dkns = hansenObjects[s].getDerivative(-n - 1, X);
final double coef0 = d0smj * jn;
final double coef1 = coef0 * lns;
final double coef2 = coef1 * kns;
final double coef3 = coef2 * ijs;
final double coef4 = coef2 * jjs;
// Add the term to the coefficients
cCoef[0][j] -= coef3;
cCoef[1][j] -= coef3 * (n + 1);
cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
cCoef[4][j] -= coef2 * dIjsdAlpha;
cCoef[5][j] -= coef2 * dIjsdBeta;
cCoef[6][j] -= coef0 * dlns * kns * ijs;
sCoef[0][j] += coef4;
sCoef[1][j] += coef4 * (n + 1);
sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
sCoef[4][j] += coef2 * dJjsdAlpha;
sCoef[5][j] += coef2 * dJjsdBeta;
sCoef[6][j] += coef0 * dlns * kns * jjs;
}
}
}
//the second double sum is added only if N >= 4 and j between 5 and 2*N-3
if (nMax >= 4 && isBetween(j, 5, 2 * nMax - 3)) {
//compute second double sum where s: j-min(j-1,N) -> (j-3)/2 and n: j-s -> min(j-1,N)
for (int s = j - minjm1on; s <= FastMath.min((j - 3) / 2, sMax); s++) {
// j - s
final int jms = j - s;
// Kronecker symbols (2 - delta(0,s-j)) and (2 - delta(0,j-s))
final int d0smj = (s == j) ? 1 : 2;
for (int n = j - s; n <= minjm1on; n++) {
// if n + (j-s) is odd, then the term is equal to zero due to the factor Vn,s-j
if ((n + jms) % 2 == 0) {
// (2 - delta(0,j-s)) * J<sub>n</sub> * K₀<sup>-n-1,s</sup> * L<sub>n</sub><sup>j-s</sup>
final double lns = lnsCoef.getLns(n, jms);
final double dlns = lnsCoef.getdLnsdGamma(n, jms);
final double ijs = ghijCoef.getIjs(s, jms);
final double dIjsdh = ghijCoef.getdIjsdh(s, jms);
final double dIjsdk = ghijCoef.getdIjsdk(s, jms);
final double dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final double dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final double jjs = ghijCoef.getJjs(s, jms);
final double dJjsdh = ghijCoef.getdJjsdh(s, jms);
final double dJjsdk = ghijCoef.getdJjsdk(s, jms);
final double dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final double dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final double jn = -harmonics.getUnnormalizedCnm(n, 0);
// K₀<sup>-n-1,s</sup>
final double kns = hansenObjects[s].getValue(-n - 1, X);
final double dkns = hansenObjects[s].getDerivative(-n - 1, X);
final double coef0 = d0smj * jn;
final double coef1 = coef0 * lns;
final double coef2 = coef1 * kns;
final double coef3 = coef2 * ijs;
final double coef4 = coef2 * jjs;
// Add the term to the coefficients
cCoef[0][j] -= coef3;
cCoef[1][j] -= coef3 * (n + 1);
cCoef[2][j] -= coef1 * (kns * dIjsdh + ijs * hXXX * dkns);
cCoef[3][j] -= coef1 * (kns * dIjsdk + ijs * kXXX * dkns);
cCoef[4][j] -= coef2 * dIjsdAlpha;
cCoef[5][j] -= coef2 * dIjsdBeta;
cCoef[6][j] -= coef0 * dlns * kns * ijs;
sCoef[0][j] += coef4;
sCoef[1][j] += coef4 * (n + 1);
sCoef[2][j] += coef1 * (kns * dJjsdh + jjs * hXXX * dkns);
sCoef[3][j] += coef1 * (kns * dJjsdk + jjs * kXXX * dkns);
sCoef[4][j] += coef2 * dJjsdAlpha;
sCoef[5][j] += coef2 * dJjsdBeta;
sCoef[6][j] += coef0 * dlns * kns * jjs;
}
}
}
}
}
}
cCoef[0][j] *= -muoa / j;
cCoef[1][j] *= muoa / ( j * a );
cCoef[2][j] *= -muoa / j;
cCoef[3][j] *= -muoa / j;
cCoef[4][j] *= -muoa / j;
cCoef[5][j] *= -muoa / j;
cCoef[6][j] *= -muoa / j;
sCoef[0][j] *= -muoa / j;
sCoef[1][j] *= muoa / ( j * a );
sCoef[2][j] *= -muoa / j;
sCoef[3][j] *= -muoa / j;
sCoef[4][j] *= -muoa / j;
sCoef[5][j] *= -muoa / j;
sCoef[6][j] *= -muoa / j;
}
}
/** Check if an index is within the accepted interval.
*
* @param index the index to check
* @param lowerBound the lower bound of the interval
* @param upperBound the upper bound of the interval
* @return true if the index is between the lower and upper bounds, false otherwise
*/
private boolean isBetween(final int index, final int lowerBound, final int upperBound) {
return index >= lowerBound && index <= upperBound;
}
/**Get the value of C<sup>j</sup>.
*
* @param j j index
* @return C<sup>j</sup>
*/
public double getCj(final int j) {
return cCoef[0][j];
}
/**Get the value of dC<sup>j</sup> / da.
*
* @param j j index
* @return dC<sup>j</sup> / da
*/
public double getdCjdA(final int j) {
return cCoef[1][j];
}
/**Get the value of dC<sup>j</sup> / dh.
*
* @param j j index
* @return dC<sup>j</sup> / dh
*/
public double getdCjdH(final int j) {
return cCoef[2][j];
}
/**Get the value of dC<sup>j</sup> / dk.
*
* @param j j index
* @return dC<sup>j</sup> / dk
*/
public double getdCjdK(final int j) {
return cCoef[3][j];
}
/**Get the value of dC<sup>j</sup> / dα.
*
* @param j j index
* @return dC<sup>j</sup> / dα
*/
public double getdCjdAlpha(final int j) {
return cCoef[4][j];
}
/**Get the value of dC<sup>j</sup> / dβ.
*
* @param j j index
* @return dC<sup>j</sup> / dβ
*/
public double getdCjdBeta(final int j) {
return cCoef[5][j];
}
/**Get the value of dC<sup>j</sup> / dγ.
*
* @param j j index
* @return dC<sup>j</sup> / dγ
*/
public double getdCjdGamma(final int j) {
return cCoef[6][j];
}
/**Get the value of S<sup>j</sup>.
*
* @param j j index
* @return S<sup>j</sup>
*/
public double getSj(final int j) {
return sCoef[0][j];
}
/**Get the value of dS<sup>j</sup> / da.
*
* @param j j index
* @return dS<sup>j</sup> / da
*/
public double getdSjdA(final int j) {
return sCoef[1][j];
}
/**Get the value of dS<sup>j</sup> / dh.
*
* @param j j index
* @return dS<sup>j</sup> / dh
*/
public double getdSjdH(final int j) {
return sCoef[2][j];
}
/**Get the value of dS<sup>j</sup> / dk.
*
* @param j j index
* @return dS<sup>j</sup> / dk
*/
public double getdSjdK(final int j) {
return sCoef[3][j];
}
/**Get the value of dS<sup>j</sup> / dα.
*
* @param j j index
* @return dS<sup>j</sup> / dα
*/
public double getdSjdAlpha(final int j) {
return sCoef[4][j];
}
/**Get the value of dS<sup>j</sup> / dβ.
*
* @param j j index
* @return dS<sup>j</sup> / dβ
*/
public double getdSjdBeta(final int j) {
return sCoef[5][j];
}
/**Get the value of dS<sup>j</sup> / dγ.
*
* @param j j index
* @return dS<sup>j</sup> / dγ
*/
public double getdSjdGamma(final int j) {
return sCoef[6][j];
}
}
/** Coefficients valid for one time slot. */
private static class Slot implements Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20160319L;
/**The coefficients D<sub>i</sub>.
* <p>
* i corresponds to the equinoctial element, as follows:
* - 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 di;
/** The coefficients C<sub>i</sub><sup>j</sup>.
* <p>
* The constant term C<sub>i</sub>⁰ is also stored in this variable at index j = 0 <br>
* 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 maxFrequencyShortPeriodics maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
*/
Slot(final int maxFrequencyShortPeriodics, final int interpolationPoints) {
final int rows = maxFrequencyShortPeriodics + 1;
di = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
cij = new ShortPeriodicsInterpolatedCoefficient[rows];
sij = new ShortPeriodicsInterpolatedCoefficient[rows];
//Initialize the arrays
for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
}
}
}
}