DSSTZonal.java
/* Copyright 2002-2024 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation.semianalytical.dsst.forces;
import java.lang.reflect.Array;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.SortedMap;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.util.CombinatoricsUtils;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.FieldSinCos;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.SinCos;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitInternalError;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
import org.orekit.frames.Frame;
import org.orekit.orbits.FieldOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.PropagationType;
import org.orekit.propagation.SpacecraftState;
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.FieldAuxiliaryElements;
import org.orekit.propagation.semianalytical.dsst.utilities.FieldCjSjCoefficient;
import org.orekit.propagation.semianalytical.dsst.utilities.FieldGHIJjsPolynomials;
import org.orekit.propagation.semianalytical.dsst.utilities.FieldLnsCoefficients;
import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
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.FieldHansenZonalLinear;
import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenZonalLinear;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.utils.FieldTimeSpanMap;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.TimeSpanMap;
/** Zonal contribution to the central body gravitational perturbation.
*
* @author Romain Di Costanzo
* @author Pascal Parraud
* @author Bryan Cazabonne (field translation)
*/
public class DSSTZonal implements DSSTForceModel {
/** Name of the prefix for short period coefficients keys. */
public static final String SHORT_PERIOD_PREFIX = "DSST-central-body-zonal-";
/** 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;
/** Number of points for interpolation. */
private static final int INTERPOLATION_POINTS = 3;
/** Truncation tolerance. */
private static final double TRUNCATION_TOLERANCE = 1e-4;
/** Central attraction scaling factor.
* <p>
* We use a power of 2 to avoid numeric noise introduction
* in the multiplications/divisions sequences.
* </p>
*/
private static final double MU_SCALE = FastMath.scalb(1.0, 32);
/** Coefficient used to define the mean disturbing function V<sub>ns</sub> coefficient. */
private final SortedMap<NSKey, Double> Vns;
/** Provider for spherical harmonics. */
private final UnnormalizedSphericalHarmonicsProvider provider;
/** Central body rotating frame (fixed with respect to body). */
private final Frame bodyFixedFrame;
/** Maximal degree to consider for harmonics potential. */
private final int maxDegree;
/** Maximal degree to consider for harmonics potential. */
private final int maxDegreeShortPeriodics;
/** 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;
/** Highest power of the eccentricity to be used in mean elements computations. */
private int maxEccPowMeanElements;
/** Highest power of the eccentricity. */
private int maxEccPow;
/** Short period terms. */
private ZonalShortPeriodicCoefficients zonalSPCoefs;
/** Short period terms. */
private Map<Field<?>, FieldZonalShortPeriodicCoefficients<?>> zonalFieldSPCoefs;
/** Driver for gravitational parameter. */
private final ParameterDriver gmParameterDriver;
/** Hansen objects. */
private HansenObjects hansen;
/** Hansen objects for field elements. */
private Map<Field<?>, FieldHansenObjects<?>> fieldHansen;
/** Constructor with default reference values.
* <p>
* When this constructor is used, maximum allowed values are used
* for the short periodic coefficients:
* </p>
* <ul>
* <li> {@link #maxDegreeShortPeriodics} is set to {@code provider.getMaxDegree()} </li>
* <li> {@link #maxEccPowShortPeriodics} is set to {@code min(provider.getMaxDegree() - 1, 4)}.
* This parameter should not exceed 4 as higher values will exceed computer capacity </li>
* <li> {@link #maxFrequencyShortPeriodics} is set to {@code 2 * provider.getMaxDegree() + 1} </li>
* </ul>
* @param bodyFixedFrame rotating body frame
* <p>If set, this frame will be used to compute the Earth pole direction and the zonal contribution.
* If not, the inertial frame from the propagator will be used.<br>
* Ideally this frame should be always set to ITRF.
* However, using TOD is advised since it's generally a good trade-off between performance and precision.
*
* @param provider provider for spherical harmonics
* @since 12.1
*/
public DSSTZonal(final Frame bodyFixedFrame, final UnnormalizedSphericalHarmonicsProvider provider) {
this(bodyFixedFrame, provider, provider.getMaxDegree(), FastMath.min(4, provider.getMaxDegree() - 1), 2 * provider.getMaxDegree() + 1);
}
/**
* Constructor with bodyFixedFrame = orbit frame, and default reference values.
* <p>
* Kept for compatibility with anterior versions.
* <p>
* Setting bodyFixedFrame to null will lead to large errors if the orbit frame is far from Earth rotating frame (GCRF, EME2000...).
* The error gets smaller as the orbit frame gets closer to Earth rotating frame (MOD, then TOD).
* @see <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1104">issue-1104 on the forge</a>
* <p>
* When this constructor is used, maximum allowed values are used
* for the short periodic coefficients:
* </p>
* <ul>
* <li> {@link #maxDegreeShortPeriodics} is set to {@code provider.getMaxDegree()} </li>
* <li> {@link #maxEccPowShortPeriodics} is set to {@code min(provider.getMaxDegree() - 1, 4)}.
* This parameter should not exceed 4 as higher values will exceed computer capacity </li>
* <li> {@link #maxFrequencyShortPeriodics} is set to {@code 2 * provider.getMaxDegree() + 1} </li>
* </ul>
* @param provider provider for spherical harmonics
* @since 10.1
*/
public DSSTZonal(final UnnormalizedSphericalHarmonicsProvider provider) {
this(null, provider);
}
/** Constructor.
* @param bodyFixedFrame rotating body frame
* <p>If set, this frame will be used to compute the Earth pole direction and the zonal contribution.
* If not, the inertial frame from the propagator will be used.<br>
* Ideally this frame should be always set to ITRF.
* However, using TOD is advised since it's a good trade-off between performance and precision.
*
* @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})
* @since 12.1
*/
public DSSTZonal(final Frame bodyFixedFrame,
final UnnormalizedSphericalHarmonicsProvider provider,
final int maxDegreeShortPeriodics,
final int maxEccPowShortPeriodics,
final int maxFrequencyShortPeriodics) {
gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
provider.getMu(), MU_SCALE,
0.0, Double.POSITIVE_INFINITY);
// Central body rotating frame
this.bodyFixedFrame = bodyFixedFrame;
// Vns coefficients
this.Vns = CoefficientsFactory.computeVns(provider.getMaxDegree() + 1);
this.provider = provider;
this.maxDegree = provider.getMaxDegree();
checkIndexRange(maxDegreeShortPeriodics, 2, provider.getMaxDegree());
this.maxDegreeShortPeriodics = maxDegreeShortPeriodics;
checkIndexRange(maxEccPowShortPeriodics, 0, maxDegreeShortPeriodics - 1);
this.maxEccPowShortPeriodics = maxEccPowShortPeriodics;
checkIndexRange(maxFrequencyShortPeriodics, 1, 2 * maxDegreeShortPeriodics + 1);
this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
// Initialize default values
this.maxEccPowMeanElements = (maxDegree == 2) ? 0 : Integer.MIN_VALUE;
zonalFieldSPCoefs = new HashMap<>();
fieldHansen = new HashMap<>();
}
/**
* Constructor with bodyFixedFrame = orbit frame.
* <p>
* Added for retro-compatibility with anterior versions for initialization.
* <p>
* Setting bodyFixedFrame to null will lead to large errors if the orbit frame is far from Earth rotating frame (GCRF, EME2000...).
* The error gets smaller as the orbit frame gets closer to Earth rotating frame (MOD, then TOD).
* @see <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1104">issue-1104 on the forge</a>
* @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})
* @since 7.2
*/
public DSSTZonal(final UnnormalizedSphericalHarmonicsProvider provider,
final int maxDegreeShortPeriodics,
final int maxEccPowShortPeriodics,
final int maxFrequencyShortPeriodics) {
this(null, provider, maxDegreeShortPeriodics, maxEccPowShortPeriodics, maxFrequencyShortPeriodics);
}
/** Check an index range.
* @param index index value
* @param min minimum value for index
* @param max maximum value for index
*/
private void checkIndexRange(final int index, final int min, final int max) {
if (index < min || index > max) {
throw new OrekitException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, index, min, 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> initializeShortPeriodTerms(final AuxiliaryElements auxiliaryElements,
final PropagationType type,
final double[] parameters) {
computeMeanElementsTruncations(auxiliaryElements, parameters);
switch (type) {
case MEAN:
maxEccPow = maxEccPowMeanElements;
break;
case OSCULATING:
maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
break;
default:
throw new OrekitInternalError(null);
}
hansen = new HansenObjects();
final List<ShortPeriodTerms> list = new ArrayList<>();
zonalSPCoefs = new ZonalShortPeriodicCoefficients(maxFrequencyShortPeriodics,
INTERPOLATION_POINTS,
new TimeSpanMap<>(new Slot(maxFrequencyShortPeriodics, INTERPOLATION_POINTS)));
list.add(zonalSPCoefs);
return list;
}
/** {@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 <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>> initializeShortPeriodTerms(final FieldAuxiliaryElements<T> auxiliaryElements,
final PropagationType type,
final T[] parameters) {
// Field used by default
final Field<T> field = auxiliaryElements.getDate().getField();
computeMeanElementsTruncations(auxiliaryElements, parameters, field);
switch (type) {
case MEAN:
maxEccPow = maxEccPowMeanElements;
break;
case OSCULATING:
maxEccPow = FastMath.max(maxEccPowMeanElements, maxEccPowShortPeriodics);
break;
default:
throw new OrekitInternalError(null);
}
fieldHansen.put(field, new FieldHansenObjects<>(field));
final FieldZonalShortPeriodicCoefficients<T> fzspc =
new FieldZonalShortPeriodicCoefficients<>(maxFrequencyShortPeriodics,
INTERPOLATION_POINTS,
new FieldTimeSpanMap<>(new FieldSlot<>(maxFrequencyShortPeriodics,
INTERPOLATION_POINTS),
field));
zonalFieldSPCoefs.put(field, fzspc);
return Collections.singletonList(fzspc);
}
/** Compute indices truncations for mean elements computations.
* @param auxiliaryElements auxiliary elements
* @param parameters values of the force model parameters for state date (only 1 value for each parameter)
*/
private void computeMeanElementsTruncations(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
final DSSTZonalContext context = new DSSTZonalContext(auxiliaryElements, bodyFixedFrame, provider, parameters);
//Compute the max eccentricity power for the mean element rate expansion
if (maxDegree == 2) {
maxEccPowMeanElements = 0;
} else {
// Initializes specific parameters.
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(auxiliaryElements.getDate());
// Utilities for truncation
final double ax2or = 2. * auxiliaryElements.getSma() / provider.getAe();
double xmuran = parameters[0] / auxiliaryElements.getSma();
// Set a lower bound for eccentricity
final double eo2 = FastMath.max(0.0025, auxiliaryElements.getEcc() / 2.);
final double x2o2 = context.getChi2() / 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] = context.getChi();
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 maxDeg = 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(maxDeg, m);
final double snm = harmonics.getUnnormalizedSnm(maxDeg, 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 = (maxDeg - maxEccPowMeanElements - 1) / 2;
int l = maxDeg - 2 * nsld2;
// Loop over l
double term = 0.;
do {
// Compute magnitude of current spherical harmonic term.
if (m < l) {
term = csnm * xmuran *
(CombinatoricsUtils.factorialDouble(maxDeg - l) / (CombinatoricsUtils.factorialDouble(maxDeg - m))) *
(CombinatoricsUtils.factorialDouble(maxDeg + l) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l))) *
eccPwr[l] * UpperBounds.getDnl(context.getChi2(), chiPwr[l], maxDeg, l) *
(UpperBounds.getRnml(context.getGamma(), maxDeg, l, m, 1, I) + UpperBounds.getRnml(context.getGamma(), maxDeg, l, m, -1, I));
} else {
term = csnm * xmuran *
(CombinatoricsUtils.factorialDouble(maxDeg + m) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l))) *
eccPwr[l] * hafPwr[m - l] * UpperBounds.getDnl(context.getChi2(), chiPwr[l], maxDeg, l) *
(UpperBounds.getRnml(context.getGamma(), maxDeg, m, l, 1, I) + UpperBounds.getRnml(context.getGamma(), maxDeg, 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 < maxDeg);
// 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(maxDeg, provider.getMaxOrder()));
// Proceed to next degree.
xmuran *= ax2or;
maxDeg--;
} while (maxDeg > maxEccPowMeanElements + 2);
maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
}
}
/** Compute indices truncations for mean elements computations.
* @param <T> type of the elements
* @param auxiliaryElements auxiliary elements
* @param parameters values of the force model parameters
* @param field field used by default
*/
private <T extends CalculusFieldElement<T>> void computeMeanElementsTruncations(final FieldAuxiliaryElements<T> auxiliaryElements,
final T[] parameters,
final Field<T> field) {
final T zero = field.getZero();
final FieldDSSTZonalContext<T> context = new FieldDSSTZonalContext<>(auxiliaryElements, bodyFixedFrame, provider, parameters);
//Compute the max eccentricity power for the mean element rate expansion
if (maxDegree == 2) {
maxEccPowMeanElements = 0;
} else {
// Initializes specific parameters.
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(auxiliaryElements.getDate().toAbsoluteDate());
// Utilities for truncation
final T ax2or = auxiliaryElements.getSma().multiply(2.).divide(provider.getAe());
T xmuran = parameters[0].divide(auxiliaryElements.getSma());
// Set a lower bound for eccentricity
final T eo2 = FastMath.max(zero.newInstance(0.0025), auxiliaryElements.getEcc().divide(2.));
final T x2o2 = context.getChi2().divide(2.);
final T[] eccPwr = MathArrays.buildArray(field, maxDegree + 1);
final T[] chiPwr = MathArrays.buildArray(field, maxDegree + 1);
final T[] hafPwr = MathArrays.buildArray(field, maxDegree + 1);
eccPwr[0] = zero.newInstance(1.);
chiPwr[0] = context.getChi();
hafPwr[0] = zero.newInstance(1.);
for (int i = 1; i <= maxDegree; i++) {
eccPwr[i] = eccPwr[i - 1].multiply(eo2);
chiPwr[i] = chiPwr[i - 1].multiply(x2o2);
hafPwr[i] = hafPwr[i - 1].multiply(0.5);
xmuran = xmuran.divide(ax2or);
}
// Set highest power of e and degree of current spherical harmonic.
maxEccPowMeanElements = 0;
int maxDeg = 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 T cnm = zero.newInstance(harmonics.getUnnormalizedCnm(maxDeg, m));
final T snm = zero.newInstance(harmonics.getUnnormalizedSnm(maxDeg, m));
final T csnm = FastMath.hypot(cnm, snm);
if (csnm.getReal() == 0.) {
break;
}
// Set magnitude of last spherical harmonic term.
T lastTerm = zero;
// Set current power of e and related indices.
int nsld2 = (maxDeg - maxEccPowMeanElements - 1) / 2;
int l = maxDeg - 2 * nsld2;
// Loop over l
T term = zero;
do {
// Compute magnitude of current spherical harmonic term.
if (m < l) {
term = csnm.multiply(xmuran).
multiply((CombinatoricsUtils.factorialDouble(maxDeg - l) / (CombinatoricsUtils.factorialDouble(maxDeg - m))) *
(CombinatoricsUtils.factorialDouble(maxDeg + l) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l)))).
multiply(eccPwr[l]).multiply(UpperBounds.getDnl(context.getChi2(), chiPwr[l], maxDeg, l)).
multiply(UpperBounds.getRnml(context.getGamma(), maxDeg, l, m, 1, I).add(UpperBounds.getRnml(context.getGamma(), maxDeg, l, m, -1, I)));
} else {
term = csnm.multiply(xmuran).
multiply(CombinatoricsUtils.factorialDouble(maxDeg + m) / (CombinatoricsUtils.factorialDouble(nsld2) * CombinatoricsUtils.factorialDouble(nsld2 + l))).
multiply(eccPwr[l]).multiply(hafPwr[m - l]).multiply(UpperBounds.getDnl(context.getChi2(), chiPwr[l], maxDeg, l)).
multiply(UpperBounds.getRnml(context.getGamma(), maxDeg, m, l, 1, I).add(UpperBounds.getRnml(context.getGamma(), maxDeg, m, l, -1, I)));
}
// Is the current spherical harmonic term bigger than the truncation tolerance ?
if (term.getReal() >= TRUNCATION_TOLERANCE) {
maxEccPowMeanElements = l;
} else {
// Is the current term smaller than the last term ?
if (term.getReal() < lastTerm.getReal()) {
break;
}
}
// Proceed to next power of e.
lastTerm = term;
l += 2;
nsld2--;
} while (l < maxDeg);
// Is the current spherical harmonic term bigger than the truncation tolerance ?
if (term.getReal() >= TRUNCATION_TOLERANCE) {
maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
return;
}
// Proceed to next order.
m++;
} while (m <= FastMath.min(maxDeg, provider.getMaxOrder()));
// Proceed to next degree.
xmuran = xmuran.multiply(ax2or);
maxDeg--;
} while (maxDeg > maxEccPowMeanElements + 2);
maxEccPowMeanElements = FastMath.min(maxDegree - 2, maxEccPowMeanElements);
}
}
/** Performs initialization at each integration step for the current force model.
* <p>
* This method aims at being called before mean elements rates computation.
* </p>
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param parameters values of the force model parameters
* @return new force model context
*/
private DSSTZonalContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
return new DSSTZonalContext(auxiliaryElements, bodyFixedFrame, provider, parameters);
}
/** Performs initialization at each integration step for the current force model.
* <p>
* This method aims at being called before mean elements rates computation.
* </p>
* @param <T> type of the elements
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param parameters values of the force model parameters
* @return new force model context
*/
private <T extends CalculusFieldElement<T>> FieldDSSTZonalContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
final T[] parameters) {
return new FieldDSSTZonalContext<>(auxiliaryElements, bodyFixedFrame, provider, parameters);
}
/** {@inheritDoc} */
@Override
public double[] getMeanElementRate(final SpacecraftState spacecraftState,
final AuxiliaryElements auxiliaryElements, final double[] parameters) {
// Container of attributes
final DSSTZonalContext context = initializeStep(auxiliaryElements, parameters);
// Access to potential U derivatives
final UAnddU udu = new UAnddU(spacecraftState.getDate(), context, auxiliaryElements, hansen);
return computeMeanElementRates(context, udu);
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> spacecraftState,
final FieldAuxiliaryElements<T> auxiliaryElements,
final T[] parameters) {
// Field used by default
final Field<T> field = auxiliaryElements.getDate().getField();
// Container of attributes
final FieldDSSTZonalContext<T> context = initializeStep(auxiliaryElements, parameters);
@SuppressWarnings("unchecked")
final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
// Access to potential U derivatives
final FieldUAnddU<T> udu = new FieldUAnddU<>(spacecraftState.getDate(), context, auxiliaryElements, fho);
return computeMeanElementRates(spacecraftState.getDate(), context, udu);
}
/** Compute the mean element rates.
* @param context container for attributes
* @param udu derivatives of the gravitational potential U
* @return the mean element rates
*/
private double[] computeMeanElementRates(final DSSTZonalContext context,
final UAnddU udu) {
// Auxiliary elements related to the current orbit
final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
// Compute cross derivatives [Eq. 2.2-(8)]
// U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
final double UAlphaGamma = context.getAlpha() * udu.getdUdGa() - context.getGamma() * udu.getdUdAl();
// U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
final double UBetaGamma = context.getBeta() * udu.getdUdGa() - context.getGamma() * udu.getdUdBe();
// Common factor
final double pUAGmIqUBGoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();
// Compute mean elements rates [Eq. 3.1-(1)]
final double da = 0.;
final double dh = context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUAGmIqUBGoAB;
final double dk = -context.getBoA() * udu.getdUdh() - auxiliaryElements.getH() * pUAGmIqUBGoAB;
final double dp = -context.getCo2AB() * UBetaGamma;
final double dq = -context.getCo2AB() * UAlphaGamma * I;
final double dM = -context.getAx2oA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUAGmIqUBGoAB;
return new double[] {da, dk, dh, dq, dp, dM};
}
/** Compute the mean element rates.
* @param <T> type of the elements
* @param date current date
* @param context container for attributes
* @param udu derivatives of the gravitational potential U
* @return the mean element rates
*/
private <T extends CalculusFieldElement<T>> T[] computeMeanElementRates(final FieldAbsoluteDate<T> date,
final FieldDSSTZonalContext<T> context,
final FieldUAnddU<T> udu) {
// Auxiliary elements related to the current orbit
final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
// Parameter for array building
final Field<T> field = date.getField();
// Compute cross derivatives [Eq. 2.2-(8)]
// U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
final T UAlphaGamma = udu.getdUdGa().multiply(context.getAlpha()).subtract(udu.getdUdAl().multiply(context.getGamma()));
// U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
final T UBetaGamma = udu.getdUdGa().multiply(context.getBeta()).subtract(udu.getdUdBe().multiply(context.getGamma()));
// Common factor
final T pUAGmIqUBGoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(I).multiply(auxiliaryElements.getQ()))).multiply(context.getOoAB());
// Compute mean elements rates [Eq. 3.1-(1)]
final T da = field.getZero();
final T dh = udu.getdUdk().multiply(context.getBoA()).add(pUAGmIqUBGoAB.multiply(auxiliaryElements.getK()));
final T dk = (udu.getdUdh().multiply(context.getBoA()).negate()).subtract(pUAGmIqUBGoAB.multiply(auxiliaryElements.getH()));
final T dp = UBetaGamma.multiply(context.getCo2AB().negate());
final T dq = UAlphaGamma.multiply(context.getCo2AB().negate()).multiply(I);
final T dM = pUAGmIqUBGoAB.add(udu.getdUda().multiply(context.getAx2oA().negate())).add((udu.getdUdh().multiply(auxiliaryElements.getH()).add(udu.getdUdk().multiply(auxiliaryElements.getK()))).multiply(context.getBoABpo()));
final T[] elements = MathArrays.buildArray(field, 6);
elements[0] = da;
elements[1] = dk;
elements[2] = dh;
elements[3] = dq;
elements[4] = dp;
elements[5] = dM;
return elements;
}
/** {@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 double[] parameters, final SpacecraftState... meanStates) {
final Slot slot = zonalSPCoefs.createSlot(meanStates);
for (final SpacecraftState meanState : meanStates) {
// Auxiliary elements related to the current orbit
final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(meanState.getOrbit(), I);
// Container of attributes
// Extract the proper parameters valid for the corresponding meanState date from the input array
final double[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());
final DSSTZonalContext context = initializeStep(auxiliaryElements, extractedParameters);
// Access to potential U derivatives
final UAnddU udu = new UAnddU(meanState.getDate(), context, auxiliaryElements, hansen);
// Compute rhoj and sigmaj
final double[][] rhoSigma = computeRhoSigmaCoefficients(slot, auxiliaryElements);
// Compute Di
computeDiCoefficients(meanState.getDate(), slot, context, udu);
// generate the Cij and Sij coefficients
final FourierCjSjCoefficients cjsj = new FourierCjSjCoefficients(meanState.getDate(),
maxDegreeShortPeriodics,
maxEccPowShortPeriodics,
maxFrequencyShortPeriodics,
context,
hansen);
computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma, context, auxiliaryElements, udu);
}
}
/** {@inheritDoc} */
@Override
@SuppressWarnings("unchecked")
public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
final FieldSpacecraftState<T>... meanStates) {
// Field used by default
final Field<T> field = meanStates[0].getDate().getField();
final FieldZonalShortPeriodicCoefficients<T> fzspc = (FieldZonalShortPeriodicCoefficients<T>) zonalFieldSPCoefs.get(field);
final FieldSlot<T> slot = fzspc.createSlot(meanStates);
for (final FieldSpacecraftState<T> meanState : meanStates) {
// Auxiliary elements related to the current orbit
final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(meanState.getOrbit(), I);
// Container of attributes
// Extract the proper parameters valid for the corresponding meanState date from the input array
final T[] extractedParameters = this.extractParameters(parameters, auxiliaryElements.getDate());
final FieldDSSTZonalContext<T> context = initializeStep(auxiliaryElements, extractedParameters);
final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
// Access to potential U derivatives
final FieldUAnddU<T> udu = new FieldUAnddU<>(meanState.getDate(), context, auxiliaryElements, fho);
// Compute rhoj and sigmaj
final T[][] rhoSigma = computeRhoSigmaCoefficients(slot, auxiliaryElements, field);
// Compute Di
computeDiCoefficients(meanState.getDate(), slot, context, field, udu);
// generate the Cij and Sij coefficients
final FieldFourierCjSjCoefficients<T> cjsj = new FieldFourierCjSjCoefficients<>(meanState.getDate(),
maxDegreeShortPeriodics,
maxEccPowShortPeriodics,
maxFrequencyShortPeriodics,
context,
fho);
computeCijSijCoefficients(meanState.getDate(), slot, cjsj, rhoSigma, context, auxiliaryElements, field, udu);
}
}
/** {@inheritDoc} */
public List<ParameterDriver> getParametersDrivers() {
return Collections.singletonList(gmParameterDriver);
}
/** Generate the values for the D<sub>i</sub> coefficients.
* @param date target date
* @param slot slot to which the coefficients belong
* @param context container for attributes
* @param udu derivatives of the gravitational potential U
*/
private void computeDiCoefficients(final AbsoluteDate date,
final Slot slot,
final DSSTZonalContext context,
final UAnddU udu) {
final double[] meanElementRates = computeMeanElementRates(context, udu);
final double[] currentDi = new double[6];
// Add the coefficients to the interpolation grid
for (int i = 0; i < 6; i++) {
currentDi[i] = meanElementRates[i] / context.getMeanMotion();
if (i == 5) {
currentDi[i] += -1.5 * 2 * udu.getU() * context.getOON2A2();
}
}
slot.di.addGridPoint(date, currentDi);
}
/** Generate the values for the D<sub>i</sub> coefficients.
* @param <T> type of the elements
* @param date target date
* @param slot slot to which the coefficients belong
* @param context container for attributes
* @param field field used by default
* @param udu derivatives of the gravitational potential U
*/
private <T extends CalculusFieldElement<T>> void computeDiCoefficients(final FieldAbsoluteDate<T> date,
final FieldSlot<T> slot,
final FieldDSSTZonalContext<T> context,
final Field<T> field,
final FieldUAnddU<T> udu) {
final T[] meanElementRates = computeMeanElementRates(date, context, udu);
final T[] currentDi = MathArrays.buildArray(field, 6);
// Add the coefficients to the interpolation grid
for (int i = 0; i < 6; i++) {
currentDi[i] = meanElementRates[i].divide(context.getMeanMotion());
if (i == 5) {
currentDi[i] = currentDi[i].add(context.getOON2A2().multiply(udu.getU()).multiply(2.).multiply(-1.5));
}
}
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>
* @param context container for attributes
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param udu derivatives of the gravitational potential U
*/
private void computeCijSijCoefficients(final AbsoluteDate date, final Slot slot,
final FourierCjSjCoefficients cjsj,
final double[][] rhoSigma, final DSSTZonalContext context,
final AuxiliaryElements auxiliaryElements,
final UAnddU udu) {
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 * auxiliaryElements.getK() * udu.getU() - context.getHK() * cjsj.getCj(1) + context.getK2MH2O2() * cjsj.getSj(1);
final double coef2 = 4 * auxiliaryElements.getH() * udu.getU() + context.getK2MH2O2() * cjsj.getCj(1) + context.getHK() * cjsj.getSj(1);
final double coef3 = (auxiliaryElements.getK() * cjsj.getCj(1) + auxiliaryElements.getH() * cjsj.getSj(1)) / 4.;
final double coef4 = (8 * udu.getU() - auxiliaryElements.getH() * cjsj.getCj(1) + auxiliaryElements.getK() * 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 = context.getK2MH2() * udu.getU();
final double coef2 = 2 * context.getHK() * udu.getU();
final double coef3 = auxiliaryElements.getH() * udu.getU() / 2;
final double coef4 = auxiliaryElements.getK() * udu.getU() / 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) && j + 2 < cjsj.jMax) {
final double coef1 = ( j + 2 ) * (-context.getHK() * cjsj.getCj(j + 2) + context.getK2MH2O2() * cjsj.getSj(j + 2));
final double coef2 = ( j + 2 ) * (context.getK2MH2O2() * cjsj.getCj(j + 2) + context.getHK() * cjsj.getSj(j + 2));
final double coef3 = ( j + 2 ) * (auxiliaryElements.getK() * cjsj.getCj(j + 2) + auxiliaryElements.getH() * cjsj.getSj(j + 2)) / 4;
final double coef4 = ( j + 2 ) * (auxiliaryElements.getH() * cjsj.getCj(j + 2) - auxiliaryElements.getK() * 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) && j + 1 < cjsj.jMax) {
final double coef1 = 2 * ( j + 1 ) * (-auxiliaryElements.getH() * cjsj.getCj(j + 1) + auxiliaryElements.getK() * cjsj.getSj(j + 1));
final double coef2 = 2 * ( j + 1 ) * (auxiliaryElements.getK() * cjsj.getCj(j + 1) + auxiliaryElements.getH() * 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) && j - 1 < cjsj.jMax) {
final double coef1 = 2 * ( j - 1 ) * (auxiliaryElements.getH() * cjsj.getCj(j - 1) + auxiliaryElements.getK() * cjsj.getSj(j - 1));
final double coef2 = 2 * ( j - 1 ) * (auxiliaryElements.getK() * cjsj.getCj(j - 1) - auxiliaryElements.getH() * 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) && j - 2 < cjsj.jMax) {
final double coef1 = ( j - 2 ) * (context.getHK() * cjsj.getCj(j - 2) + context.getK2MH2O2() * cjsj.getSj(j - 2));
final double coef2 = ( j - 2 ) * (-context.getK2MH2O2() * cjsj.getCj(j - 2) + context.getHK() * cjsj.getSj(j - 2));
final double coef3 = ( j - 2 ) * (auxiliaryElements.getK() * cjsj.getCj(j - 2) - auxiliaryElements.getH() * cjsj.getSj(j - 2)) / 4;
final double coef4 = ( j - 2 ) * (auxiliaryElements.getH() * cjsj.getCj(j - 2) + auxiliaryElements.getK() * cjsj.getSj(j - 2)) / 4;
final double coef5 = ( j - 2 ) * (context.getK2MH2O2() * cjsj.getCj(j - 2) - context.getHK() * 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] *= context.getX3ON2A();
currentSij[0] *= context.getX3ON2A();
//for k (i == 1) -> χ / (n² * a²)
currentCij[1] *= context.getXON2A2();
currentSij[1] *= context.getXON2A2();
//for h (i == 2) -> χ / (n² * a²)
currentCij[2] *= context.getXON2A2();
currentSij[2] *= context.getXON2A2();
//for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
currentCij[5] *= context.getX2ON2A2XP1();
currentSij[5] *= context.getX2ON2A2XP1();
// j is between 1 and 2 * N - 1
if (isBetween(j, 1, 2 * nMax - 1) && j < cjsj.jMax) {
// Compute cross derivatives
// Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
final double CjAlphaGamma = context.getAlpha() * cjsj.getdCjdGamma(j) - context.getGamma() * cjsj.getdCjdAlpha(j);
// Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
final double CjAlphaBeta = context.getAlpha() * cjsj.getdCjdBeta(j) - context.getBeta() * cjsj.getdCjdAlpha(j);
// Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
final double CjBetaGamma = context.getBeta() * cjsj.getdCjdGamma(j) - context.getGamma() * cjsj.getdCjdBeta(j);
// Cj(h,k) = h * dC/dk - k * dC/dh
final double CjHK = auxiliaryElements.getH() * cjsj.getdCjdK(j) - auxiliaryElements.getK() * cjsj.getdCjdH(j);
// Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
final double SjAlphaGamma = context.getAlpha() * cjsj.getdSjdGamma(j) - context.getGamma() * cjsj.getdSjdAlpha(j);
// Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
final double SjAlphaBeta = context.getAlpha() * cjsj.getdSjdBeta(j) - context.getBeta() * cjsj.getdSjdAlpha(j);
// Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
final double SjBetaGamma = context.getBeta() * cjsj.getdSjdGamma(j) - context.getGamma() * cjsj.getdSjdBeta(j);
// Sj(h,k) = h * dS/dk - k * dS/dh
final double SjHK = auxiliaryElements.getH() * cjsj.getdSjdK(j) - auxiliaryElements.getK() * cjsj.getdSjdH(j);
//Coefficients for a
final double coef1 = context.getX3ON2A() * (3 - context.getBB()) * j;
currentCij[0] += coef1 * cjsj.getSj(j);
currentSij[0] -= coef1 * cjsj.getCj(j);
//Coefficients for k and h
final double coef2 = auxiliaryElements.getP() * CjAlphaGamma - I * auxiliaryElements.getQ() * CjBetaGamma;
final double coef3 = auxiliaryElements.getP() * SjAlphaGamma - I * auxiliaryElements.getQ() * SjBetaGamma;
currentCij[1] -= context.getXON2A2() * (auxiliaryElements.getH() * coef2 + context.getBB() * cjsj.getdCjdH(j) - 1.5 * auxiliaryElements.getK() * j * cjsj.getSj(j));
currentSij[1] -= context.getXON2A2() * (auxiliaryElements.getH() * coef3 + context.getBB() * cjsj.getdSjdH(j) + 1.5 * auxiliaryElements.getK() * j * cjsj.getCj(j));
currentCij[2] += context.getXON2A2() * (auxiliaryElements.getK() * coef2 + context.getBB() * cjsj.getdCjdK(j) + 1.5 * auxiliaryElements.getH() * j * cjsj.getSj(j));
currentSij[2] += context.getXON2A2() * (auxiliaryElements.getK() * coef3 + context.getBB() * cjsj.getdSjdK(j) - 1.5 * auxiliaryElements.getH() * 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] = context.getCXO2N2A2() * (-I * CjAlphaGamma + auxiliaryElements.getQ() * coef4);
currentSij[3] = context.getCXO2N2A2() * (-I * SjAlphaGamma + auxiliaryElements.getQ() * coef5);
currentCij[4] = context.getCXO2N2A2() * (-CjBetaGamma + auxiliaryElements.getP() * coef4);
currentSij[4] = context.getCXO2N2A2() * (-SjBetaGamma + auxiliaryElements.getP() * coef5);
//Coefficients for λ
final double coef6 = auxiliaryElements.getH() * cjsj.getdCjdH(j) + auxiliaryElements.getK() * cjsj.getdCjdK(j);
final double coef7 = auxiliaryElements.getH() * cjsj.getdSjdH(j) + auxiliaryElements.getK() * cjsj.getdSjdK(j);
currentCij[5] += context.getOON2A2() * (-2 * auxiliaryElements.getSma() * cjsj.getdCjdA(j) + coef6 / (context.getChi() + 1) + context.getChi() * coef2 - 3 * cjsj.getCj(j));
currentSij[5] += context.getOON2A2() * (-2 * auxiliaryElements.getSma() * cjsj.getdSjdA(j) + coef7 / (context.getChi() + 1) + context.getChi() * 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);
}
/**
* Generate the values for the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
* @param <T> type of the elements
* @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>
* @param context container for attributes
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param field field used by default
* @param udu derivatives of the gravitational potential U
*/
private <T extends CalculusFieldElement<T>> void computeCijSijCoefficients(final FieldAbsoluteDate<T> date,
final FieldSlot<T> slot,
final FieldFourierCjSjCoefficients<T> cjsj,
final T[][] rhoSigma,
final FieldDSSTZonalContext<T> context,
final FieldAuxiliaryElements<T> auxiliaryElements,
final Field<T> field,
final FieldUAnddU<T> udu) {
// Zero
final T zero = field.getZero();
final int nMax = maxDegreeShortPeriodics;
// The C<sub>i</sub>⁰ coefficients
final T[] currentCi0 = MathArrays.buildArray(field, 6);
Arrays.fill(currentCi0, zero);
for (int j = 1; j < slot.cij.length; j++) {
// Create local arrays
final T[] currentCij = MathArrays.buildArray(field, 6);
final T[] currentSij = MathArrays.buildArray(field, 6);
Arrays.fill(currentCij, zero);
Arrays.fill(currentSij, zero);
// j == 1
if (j == 1) {
final T coef1 = auxiliaryElements.getK().multiply(udu.getU()).multiply(4.).subtract(context.getHK().multiply(cjsj.getCj(1))).add(context.getK2MH2O2().multiply(cjsj.getSj(1)));
final T coef2 = auxiliaryElements.getH().multiply(udu.getU()).multiply(4.).add(context.getK2MH2O2().multiply(cjsj.getCj(1))).add(context.getHK().multiply(cjsj.getSj(1)));
final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(1)).add(auxiliaryElements.getH().multiply(cjsj.getSj(1))).divide(4.);
final T coef4 = udu.getU().multiply(8.).subtract(auxiliaryElements.getH().multiply(cjsj.getCj(1))).add(auxiliaryElements.getK().multiply(cjsj.getSj(1))).divide(4.);
//Coefficients for a
currentCij[0] = currentCij[0].add(coef1);
currentSij[0] = currentSij[0].add(coef2);
//Coefficients for k
currentCij[1] = currentCij[1].add(coef4);
currentSij[1] = currentSij[1].add(coef3);
//Coefficients for h
currentCij[2] = currentCij[2].subtract(coef3);
currentSij[2] = currentSij[2].add(coef4);
//Coefficients for λ
currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
currentSij[5] = currentSij[5].add(coef1.divide(2.));
}
// j == 2
if (j == 2) {
final T coef1 = context.getK2MH2().multiply(udu.getU());
final T coef2 = context.getHK().multiply(udu.getU()).multiply(2.);
final T coef3 = auxiliaryElements.getH().multiply(udu.getU()).divide(2.);
final T coef4 = auxiliaryElements.getK().multiply(udu.getU()).divide(2.);
//Coefficients for a
currentCij[0] = currentCij[0].add(coef1);
currentSij[0] = currentSij[0].add(coef2);
//Coefficients for k
currentCij[1] = currentCij[1].add(coef4);
currentSij[1] = currentSij[1].add(coef3);
//Coefficients for h
currentCij[2] = currentCij[2].subtract(coef3);
currentSij[2] = currentSij[2].add(coef4);
//Coefficients for λ
currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
currentSij[5] = currentSij[5].add(coef1.divide(2.));
}
// j between 1 and 2N-3
if (isBetween(j, 1, 2 * nMax - 3) && j + 2 < cjsj.jMax) {
final T coef1 = context.getHK().negate().multiply(cjsj.getCj(j + 2)).add(context.getK2MH2O2().multiply(cjsj.getSj(j + 2))).multiply(j + 2);
final T coef2 = context.getK2MH2O2().multiply(cjsj.getCj(j + 2)).add(context.getHK().multiply(cjsj.getSj(j + 2))).multiply(j + 2);
final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(j + 2)).add(auxiliaryElements.getH().multiply(cjsj.getSj(j + 2))).multiply(j + 2).divide(4.);
final T coef4 = auxiliaryElements.getH().multiply(cjsj.getCj(j + 2)).subtract(auxiliaryElements.getK().multiply(cjsj.getSj(j + 2))).multiply(j + 2).divide(4.);
//Coefficients for a
currentCij[0] = currentCij[0].add(coef1);
currentSij[0] = currentSij[0].subtract(coef2);
//Coefficients for k
currentCij[1] = currentCij[1].add(coef4.negate());
currentSij[1] = currentSij[1].subtract(coef3);
//Coefficients for h
currentCij[2] = currentCij[2].subtract(coef3);
currentSij[2] = currentSij[2].add(coef4);
//Coefficients for λ
currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
currentSij[5] = currentSij[5].add(coef1.negate().divide(2.));
}
// j between 1 and 2N-2
if (isBetween(j, 1, 2 * nMax - 2) && j + 1 < cjsj.jMax) {
final T coef1 = auxiliaryElements.getH().negate().multiply(cjsj.getCj(j + 1)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j + 1))).multiply(2. * (j + 1));
final T coef2 = auxiliaryElements.getK().multiply(cjsj.getCj(j + 1)).add(auxiliaryElements.getH().multiply(cjsj.getSj(j + 1))).multiply(2. * (j + 1));
final T coef3 = cjsj.getCj(j + 1).multiply(j + 1);
final T coef4 = cjsj.getSj(j + 1).multiply(j + 1);
//Coefficients for a
currentCij[0] = currentCij[0].add(coef1);
currentSij[0] = currentSij[0].subtract(coef2);
//Coefficients for k
currentCij[1] = currentCij[1].add(coef4);
currentSij[1] = currentSij[1].subtract(coef3);
//Coefficients for h
currentCij[2] = currentCij[2].subtract(coef3);
currentSij[2] = currentSij[2].subtract(coef4);
//Coefficients for λ
currentCij[5] = currentCij[5].subtract(coef2.divide(2.));
currentSij[5] = currentSij[5].add(coef1.negate().divide(2.));
}
// j between 2 and 2N
if (isBetween(j, 2, 2 * nMax) && j - 1 < cjsj.jMax) {
final T coef1 = auxiliaryElements.getH().multiply(cjsj.getCj(j - 1)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j - 1))).multiply(2 * ( j - 1 ));
final T coef2 = auxiliaryElements.getK().multiply(cjsj.getCj(j - 1)).subtract(auxiliaryElements.getH().multiply(cjsj.getSj(j - 1))).multiply(2 * ( j - 1 ));
final T coef3 = cjsj.getCj(j - 1).multiply(j - 1);
final T coef4 = cjsj.getSj(j - 1).multiply(j - 1);
//Coefficients for a
currentCij[0] = currentCij[0].add(coef1);
currentSij[0] = currentSij[0].subtract(coef2);
//Coefficients for k
currentCij[1] = currentCij[1].add(coef4);
currentSij[1] = currentSij[1].subtract(coef3);
//Coefficients for h
currentCij[2] = currentCij[2].add(coef3);
currentSij[2] = currentSij[2].add(coef4);
//Coefficients for λ
currentCij[5] = currentCij[5].add(coef2.divide(2.));
currentSij[5] = currentSij[5].add(coef1.divide(2.));
}
// j between 3 and 2N + 1
if (isBetween(j, 3, 2 * nMax + 1) && j - 2 < cjsj.jMax) {
final T coef1 = context.getHK().multiply(cjsj.getCj(j - 2)).add(context.getK2MH2O2().multiply(cjsj.getSj(j - 2))).multiply(j - 2);
final T coef2 = context.getK2MH2O2().negate().multiply(cjsj.getCj(j - 2)).add(context.getHK().multiply(cjsj.getSj(j - 2))).multiply(j - 2);
final T coef3 = auxiliaryElements.getK().multiply(cjsj.getCj(j - 2)).subtract(auxiliaryElements.getH().multiply(cjsj.getSj(j - 2))).multiply(j - 2).divide(4.);
final T coef4 = auxiliaryElements.getH().multiply(cjsj.getCj(j - 2)).add(auxiliaryElements.getK().multiply(cjsj.getSj(j - 2))).multiply(j - 2).divide(4.);
final T coef5 = context.getK2MH2O2().multiply(cjsj.getCj(j - 2)).subtract(context.getHK().multiply(cjsj.getSj(j - 2))).multiply(j - 2);
//Coefficients for a
currentCij[0] = currentCij[0].add(coef1);
currentSij[0] = currentSij[0].add(coef2);
//Coefficients for k
currentCij[1] = currentCij[1].add(coef4);
currentSij[1] = currentSij[1].add(coef3.negate());
//Coefficients for h
currentCij[2] = currentCij[2].add(coef3);
currentSij[2] = currentSij[2].add(coef4);
//Coefficients for λ
currentCij[5] = currentCij[5].add(coef5.divide(2.));
currentSij[5] = currentSij[5].add(coef1.divide(2.));
}
//multiply by the common factor
//for a (i == 0) -> χ³ / (n² * a)
currentCij[0] = currentCij[0].multiply(context.getX3ON2A());
currentSij[0] = currentSij[0].multiply(context.getX3ON2A());
//for k (i == 1) -> χ / (n² * a²)
currentCij[1] = currentCij[1].multiply(context.getXON2A2());
currentSij[1] = currentSij[1].multiply(context.getXON2A2());
//for h (i == 2) -> χ / (n² * a²)
currentCij[2] = currentCij[2].multiply(context.getXON2A2());
currentSij[2] = currentSij[2].multiply(context.getXON2A2());
//for λ (i == 5) -> (χ²) / (n² * a² * (χ + 1 ) )
currentCij[5] = currentCij[5].multiply(context.getX2ON2A2XP1());
currentSij[5] = currentSij[5].multiply(context.getX2ON2A2XP1());
// j is between 1 and 2 * N - 1
if (isBetween(j, 1, 2 * nMax - 1) && j < cjsj.jMax) {
// Compute cross derivatives
// Cj(alpha,gamma) = alpha * dC/dgamma - gamma * dC/dalpha
final T CjAlphaGamma = context.getAlpha().multiply(cjsj.getdCjdGamma(j)).subtract(context.getGamma().multiply(cjsj.getdCjdAlpha(j)));
// Cj(alpha,beta) = alpha * dC/dbeta - beta * dC/dalpha
final T CjAlphaBeta = context.getAlpha().multiply(cjsj.getdCjdBeta(j)).subtract(context.getBeta().multiply(cjsj.getdCjdAlpha(j)));
// Cj(beta,gamma) = beta * dC/dgamma - gamma * dC/dbeta
final T CjBetaGamma = context.getBeta().multiply(cjsj.getdCjdGamma(j)).subtract(context.getGamma().multiply(cjsj.getdCjdBeta(j)));
// Cj(h,k) = h * dC/dk - k * dC/dh
final T CjHK = auxiliaryElements.getH().multiply(cjsj.getdCjdK(j)).subtract(auxiliaryElements.getK().multiply(cjsj.getdCjdH(j)));
// Sj(alpha,gamma) = alpha * dS/dgamma - gamma * dS/dalpha
final T SjAlphaGamma = context.getAlpha().multiply(cjsj.getdSjdGamma(j)).subtract(context.getGamma().multiply(cjsj.getdSjdAlpha(j)));
// Sj(alpha,beta) = alpha * dS/dbeta - beta * dS/dalpha
final T SjAlphaBeta = context.getAlpha().multiply(cjsj.getdSjdBeta(j)).subtract(context.getBeta().multiply(cjsj.getdSjdAlpha(j)));
// Sj(beta,gamma) = beta * dS/dgamma - gamma * dS/dbeta
final T SjBetaGamma = context.getBeta().multiply(cjsj.getdSjdGamma(j)).subtract(context.getGamma().multiply(cjsj.getdSjdBeta(j)));
// Sj(h,k) = h * dS/dk - k * dS/dh
final T SjHK = auxiliaryElements.getH().multiply(cjsj.getdSjdK(j)).subtract(auxiliaryElements.getK().multiply(cjsj.getdSjdH(j)));
//Coefficients for a
final T coef1 = context.getX3ON2A().multiply(context.getBB().negate().add(3.)).multiply(j);
currentCij[0] = currentCij[0].add(coef1.multiply(cjsj.getSj(j)));
currentSij[0] = currentSij[0].subtract(coef1.multiply(cjsj.getCj(j)));
//Coefficients for k and h
final T coef2 = auxiliaryElements.getP().multiply(CjAlphaGamma).subtract(auxiliaryElements.getQ().multiply(CjBetaGamma).multiply(I));
final T coef3 = auxiliaryElements.getP().multiply(SjAlphaGamma).subtract(auxiliaryElements.getQ().multiply(SjBetaGamma).multiply(I));
currentCij[1] = currentCij[1].subtract(context.getXON2A2().multiply(auxiliaryElements.getH().multiply(coef2).add(context.getBB().multiply(cjsj.getdCjdH(j))).subtract(auxiliaryElements.getK().multiply(1.5).multiply(j).multiply(cjsj.getSj(j)))));
currentSij[1] = currentSij[1].subtract(context.getXON2A2().multiply(auxiliaryElements.getH().multiply(coef3).add(context.getBB().multiply(cjsj.getdSjdH(j))).add(auxiliaryElements.getK().multiply(1.5).multiply(j).multiply(cjsj.getCj(j)))));
currentCij[2] = currentCij[2].add(context.getXON2A2().multiply(auxiliaryElements.getK().multiply(coef2).add(context.getBB().multiply(cjsj.getdCjdK(j))).add(auxiliaryElements.getH().multiply(1.5).multiply(j).multiply(cjsj.getSj(j)))));
currentSij[2] = currentSij[2].add(context.getXON2A2().multiply(auxiliaryElements.getK().multiply(coef3).add(context.getBB().multiply(cjsj.getdSjdK(j))).subtract(auxiliaryElements.getH().multiply(1.5).multiply(j).multiply(cjsj.getCj(j)))));
//Coefficients for q and p
final T coef4 = CjHK.subtract(CjAlphaBeta).subtract(cjsj.getSj(j).multiply(j));
final T coef5 = SjHK.subtract(SjAlphaBeta).add(cjsj.getCj(j).multiply(j));
currentCij[3] = context.getCXO2N2A2().multiply(CjAlphaGamma.multiply(-I).add(auxiliaryElements.getQ().multiply(coef4)));
currentSij[3] = context.getCXO2N2A2().multiply(SjAlphaGamma.multiply(-I).add(auxiliaryElements.getQ().multiply(coef5)));
currentCij[4] = context.getCXO2N2A2().multiply(CjBetaGamma.negate().add(auxiliaryElements.getP().multiply(coef4)));
currentSij[4] = context.getCXO2N2A2().multiply(SjBetaGamma.negate().add(auxiliaryElements.getP().multiply(coef5)));
//Coefficients for λ
final T coef6 = auxiliaryElements.getH().multiply(cjsj.getdCjdH(j)).add(auxiliaryElements.getK().multiply(cjsj.getdCjdK(j)));
final T coef7 = auxiliaryElements.getH().multiply(cjsj.getdSjdH(j)).add(auxiliaryElements.getK().multiply(cjsj.getdSjdK(j)));
currentCij[5] = currentCij[5].add(context.getOON2A2().multiply(auxiliaryElements.getSma().multiply(-2.).multiply(cjsj.getdCjdA(j)).add(coef6.divide(context.getChi().add(1.))).add(context.getChi().multiply(coef2)).subtract(cjsj.getCj(j).multiply(3.))));
currentSij[5] = currentSij[5].add(context.getOON2A2().multiply(auxiliaryElements.getSma().multiply(-2.).multiply(cjsj.getdSjdA(j)).add(coef7.divide(context.getChi().add(1.))).add(context.getChi().multiply(coef3)).subtract(cjsj.getSj(j).multiply(3.))));
}
for (int i = 0; i < 6; i++) {
//Add the current coefficients contribution to C<sub>i</sub>⁰
currentCi0[i] = currentCi0[i].subtract(currentCij[i].multiply(rhoSigma[j][0]).add(currentSij[i].multiply(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 slot slot to which the coefficients belong
* @param auxiliaryElements auxiliary elements related to the current orbit
* @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
*/
private double[][] computeRhoSigmaCoefficients(final Slot slot,
final AuxiliaryElements auxiliaryElements) {
final CjSjCoefficient cjsjKH = new CjSjCoefficient(auxiliaryElements.getK(), auxiliaryElements.getH());
final double b = 1. / (1 + auxiliaryElements.getB());
// (-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 * auxiliaryElements.getB()) * 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;
}
/**
* 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 <T> type of the elements
* @param slot slot to which the coefficients belong
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param field field used by default
* @return array containing ρ<sub>j</sub> and σ<sub>j</sub>
*/
private <T extends CalculusFieldElement<T>> T[][] computeRhoSigmaCoefficients(final FieldSlot<T> slot,
final FieldAuxiliaryElements<T> auxiliaryElements,
final Field<T> field) {
final T zero = field.getZero();
final FieldCjSjCoefficient<T> cjsjKH = new FieldCjSjCoefficient<>(auxiliaryElements.getK(), auxiliaryElements.getH(), field);
final T b = auxiliaryElements.getB().add(1.).reciprocal();
// (-b)<sup>j</sup>
T mbtj = zero.newInstance(1.);
final T[][] rhoSigma = MathArrays.buildArray(field, slot.cij.length, 2);
for (int j = 1; j < rhoSigma.length; j++) {
//Compute current rho and sigma;
mbtj = mbtj.multiply(b.negate());
final T coef = mbtj.multiply(auxiliaryElements.getB().multiply(j).add(1.));
final T rho = coef.multiply(cjsjKH.getCj(j));
final T sigma = coef.multiply(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 {
/** 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 maxFrequencyShortPeriodics maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
* @param slots all coefficients slots
*/
ZonalShortPeriodicCoefficients(final int maxFrequencyShortPeriodics, final int interpolationPoints,
final TimeSpanMap<Slot> slots) {
// Save parameters
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();
final int compare = first.compareTo(last);
if (compare < 0) {
slots.addValidAfter(slot, first, false);
} else if (compare > 0) {
slots.addValidBefore(slot, first, false);
} else {
// single date, valid for all time
slots.addValidAfter(slot, AbsoluteDate.PAST_INFINITY, false);
}
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();
// 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 <= maxFrequencyShortPeriodics; j++) {
final double[] c = slot.cij[j].value(meanOrbit.getDate());
final double[] s = slot.sij[j].value(meanOrbit.getDate());
final SinCos sc = FastMath.sinCos(j * L);
for (int i = 0; i < 6; i++) {
// add corresponding term to the short periodic variation
shortPeriodicVariation[i] += c[i] * sc.cos();
shortPeriodicVariation[i] += s[i] * sc.sin();
}
}
return shortPeriodicVariation;
}
/** {@inheritDoc} */
@Override
public String getCoefficientsKeyPrefix() {
return DSSTZonal.SHORT_PERIOD_PREFIX;
}
/** {@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) {
// select the coefficients slot
final Slot slot = slots.get(date);
final Map<String, double[]> coefficients = new HashMap<>(2 * maxFrequencyShortPeriodics + 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 <= maxFrequencyShortPeriodics; 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);
}
}
}
/** 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 FieldZonalShortPeriodicCoefficients <T extends CalculusFieldElement<T>> implements FieldShortPeriodTerms<T> {
/** 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 FieldTimeSpanMap<FieldSlot<T>, T> slots;
/** Constructor.
* @param maxFrequencyShortPeriodics maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
* @param slots all coefficients slots
*/
FieldZonalShortPeriodicCoefficients(final int maxFrequencyShortPeriodics, final int interpolationPoints,
final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
// Save parameters
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
*/
@SuppressWarnings("unchecked")
public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
final FieldSlot<T> slot = new FieldSlot<>(maxFrequencyShortPeriodics, interpolationPoints);
final FieldAbsoluteDate<T> first = meanStates[0].getDate();
final FieldAbsoluteDate<T> 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 T[] value(final FieldOrbit<T> meanOrbit) {
// select the coefficients slot
final FieldSlot<T> slot = slots.get(meanOrbit.getDate());
// Get the True longitude L
final T L = meanOrbit.getLv();
// Compute the center
final T center = L.subtract(meanOrbit.getLM());
// Initialize short periodic variations
final T[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
final T[] d = slot.di.value(meanOrbit.getDate());
for (int i = 0; i < 6; i++) {
shortPeriodicVariation[i] = shortPeriodicVariation[i].add(center.multiply(d[i]));
}
for (int j = 1; j <= maxFrequencyShortPeriodics; j++) {
final T[] c = slot.cij[j].value(meanOrbit.getDate());
final T[] s = slot.sij[j].value(meanOrbit.getDate());
final FieldSinCos<T> sc = FastMath.sinCos(L.multiply(j));
for (int i = 0; i < 6; i++) {
// add corresponding term to the short periodic variation
shortPeriodicVariation[i] = shortPeriodicVariation[i].add(c[i].multiply(sc.cos()));
shortPeriodicVariation[i] = shortPeriodicVariation[i].add(s[i].multiply(sc.sin()));
}
}
return shortPeriodicVariation;
}
/** {@inheritDoc} */
@Override
public String getCoefficientsKeyPrefix() {
return DSSTZonal.SHORT_PERIOD_PREFIX;
}
/** {@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, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {
// select the coefficients slot
final FieldSlot<T> slot = slots.get(date);
final Map<String, T[]> coefficients = new HashMap<>(2 * maxFrequencyShortPeriodics + 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 <= maxFrequencyShortPeriodics; 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, T[]> map, final Set<String> selected,
final T[] 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);
}
}
}
/** 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
* @param context container for attributes
* @param hansenObjects initialization of hansen objects
*/
FourierCjSjCoefficients(final AbsoluteDate date,
final int nMax, final int sMax, final int jMax, final DSSTZonalContext context,
final HansenObjects hansenObjects) {
final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
this.ghijCoef = new GHIJjsPolynomials(auxiliaryElements.getK(), auxiliaryElements.getH(), context.getAlpha(), context.getBeta());
// Qns coefficients
final double[][] Qns = CoefficientsFactory.computeQns(context.getGamma(), nMax, nMax);
this.lnsCoef = new LnsCoefficients(nMax, nMax, Qns, Vns, context.getRoa());
this.nMax = nMax;
this.sMax = sMax;
this.jMax = jMax;
// compute the common factors that depends on the mean elements
this.hXXX = auxiliaryElements.getH() * context.getChi3();
this.kXXX = auxiliaryElements.getK() * context.getChi3();
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.computeHansenObjectsInitValues(context, s);
}
generateCoefficients(date, context, auxiliaryElements, hansenObjects);
}
/** Generate all coefficients.
* @param date the current date
* @param context container for attributes
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param hansenObjects initialization of hansen objects
*/
private void generateCoefficients(final AbsoluteDate date,
final DSSTZonalContext context,
final AuxiliaryElements auxiliaryElements,
final HansenObjects hansenObjects) {
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.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final double dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
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.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final double dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
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.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final double dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
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.getHansenObjects()[0].getValue(-j - 1, context.getChi());
double dkns = hansenObjects.getHansenObjects()[0].getDerivative(-j - 1, context.getChi());
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.getHansenObjects()[s].getValue(-j - 1, context.getChi());
dkns = hansenObjects.getHansenObjects()[s].getDerivative(-j - 1, context.getChi());
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.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final double dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
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.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final double dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
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.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final double dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
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.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final double dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
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] *= -context.getMuoa() / j;
cCoef[1][j] *= context.getMuoa() / ( j * auxiliaryElements.getSma() );
cCoef[2][j] *= -context.getMuoa() / j;
cCoef[3][j] *= -context.getMuoa() / j;
cCoef[4][j] *= -context.getMuoa() / j;
cCoef[5][j] *= -context.getMuoa() / j;
cCoef[6][j] *= -context.getMuoa() / j;
sCoef[0][j] *= -context.getMuoa() / j;
sCoef[1][j] *= context.getMuoa() / ( j * auxiliaryElements.getSma() );
sCoef[2][j] *= -context.getMuoa() / j;
sCoef[3][j] *= -context.getMuoa() / j;
sCoef[4][j] *= -context.getMuoa() / j;
sCoef[5][j] *= -context.getMuoa() / j;
sCoef[6][j] *= -context.getMuoa() / 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];
}
}
/** 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 FieldFourierCjSjCoefficients <T extends CalculusFieldElement<T>> {
/** The G<sub>js</sub>, H<sub>js</sub>, I<sub>js</sub> and J<sub>js</sub> polynomials. */
private final FieldGHIJjsPolynomials<T> ghijCoef;
/** L<sub>n</sub><sup>s</sup>(γ). */
private final FieldLnsCoefficients<T> 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 T[][] 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 T[][] sCoef;
/** h * Χ³. */
private final T hXXX;
/** k * Χ³. */
private final T 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
* @param context container for attributes
* @param hansenObjects initialization of hansen objects
*/
FieldFourierCjSjCoefficients(final FieldAbsoluteDate<T> date,
final int nMax, final int sMax, final int jMax,
final FieldDSSTZonalContext<T> context,
final FieldHansenObjects<T> hansenObjects) {
//Field used by default
final Field<T> field = date.getField();
final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
this.ghijCoef = new FieldGHIJjsPolynomials<>(auxiliaryElements.getK(), auxiliaryElements.getH(), context.getAlpha(), context.getBeta());
// Qns coefficients
final T[][] Qns = CoefficientsFactory.computeQns(context.getGamma(), nMax, nMax);
this.lnsCoef = new FieldLnsCoefficients<>(nMax, nMax, Qns, Vns, context.getRoa(), field);
this.nMax = nMax;
this.sMax = sMax;
this.jMax = jMax;
// compute the common factors that depends on the mean elements
this.hXXX = auxiliaryElements.getH().multiply(context.getChi3());
this.kXXX = auxiliaryElements.getK().multiply(context.getChi3());
this.cCoef = MathArrays.buildArray(field, 7, jMax + 1);
this.sCoef = MathArrays.buildArray(field, 7, jMax + 1);
for (int s = 0; s <= sMax; s++) {
//Initialise the Hansen roots
hansenObjects.computeHansenObjectsInitValues(context, s);
}
generateCoefficients(date, context, auxiliaryElements, hansenObjects, field);
}
/** Generate all coefficients.
* @param date the current date
* @param context container for attributes
* @param hansenObjects initialization of hansen objects
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param field field used by default
*/
private void generateCoefficients(final FieldAbsoluteDate<T> date,
final FieldDSSTZonalContext<T> context,
final FieldAuxiliaryElements<T> auxiliaryElements,
final FieldHansenObjects<T> hansenObjects,
final Field<T> field) {
//Zero
final T zero = field.getZero();
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());
for (int j = 1; j <= jMax; j++) {
//init arrays
for (int i = 0; i <= 6; i++) {
cCoef[i][j] = zero;
sCoef[i][j] = zero;
}
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 T lns = lnsCoef.getLns(n, -jms);
final T dlns = lnsCoef.getdLnsdGamma(n, -jms);
final T hjs = ghijCoef.getHjs(s, -jms);
final T dHjsdh = ghijCoef.getdHjsdh(s, -jms);
final T dHjsdk = ghijCoef.getdHjsdk(s, -jms);
final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, -jms);
final T dHjsdBeta = ghijCoef.getdHjsdBeta(s, -jms);
final T gjs = ghijCoef.getGjs(s, -jms);
final T dGjsdh = ghijCoef.getdGjsdh(s, -jms);
final T dGjsdk = ghijCoef.getdGjsdk(s, -jms);
final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, -jms);
final T dGjsdBeta = ghijCoef.getdGjsdBeta(s, -jms);
// J<sub>n</sub>
final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));
// K₀<sup>-n-1,s</sup>
final T kns = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final T dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
final T coef0 = jn.multiply(d0smj);
final T coef1 = coef0.multiply(lns);
final T coef2 = coef1.multiply(kns);
final T coef3 = coef2.multiply(hjs);
final T coef4 = coef2.multiply(gjs);
// Add the term to the coefficients
cCoef[0][j] = cCoef[0][j].add(coef3);
cCoef[1][j] = cCoef[1][j].add(coef3.multiply(n + 1));
cCoef[2][j] = cCoef[2][j].add(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
cCoef[3][j] = cCoef[3][j].add(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
cCoef[4][j] = cCoef[4][j].add(coef2.multiply(dHjsdAlpha));
cCoef[5][j] = cCoef[5][j].add(coef2.multiply(dHjsdBeta));
cCoef[6][j] = cCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(hjs));
sCoef[0][j] = sCoef[0][j].add(coef4);
sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(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 T lns = lnsCoef.getLns(n, jps);
final T dlns = lnsCoef.getdLnsdGamma(n, jps);
final T hjs = ghijCoef.getHjs(s, jps);
final T dHjsdh = ghijCoef.getdHjsdh(s, jps);
final T dHjsdk = ghijCoef.getdHjsdk(s, jps);
final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(s, jps);
final T dHjsdBeta = ghijCoef.getdHjsdBeta(s, jps);
final T gjs = ghijCoef.getGjs(s, jps);
final T dGjsdh = ghijCoef.getdGjsdh(s, jps);
final T dGjsdk = ghijCoef.getdGjsdk(s, jps);
final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(s, jps);
final T dGjsdBeta = ghijCoef.getdGjsdBeta(s, jps);
// J<sub>n</sub>
final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));
// K₀<sup>-n-1,s</sup>
final T kns = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final T dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
final T coef0 = jn.multiply(d0spj);
final T coef1 = coef0.multiply(lns);
final T coef2 = coef1.multiply(kns);
final T coef3 = coef2.multiply(hjs);
final T coef4 = coef2.multiply(gjs);
// Add the term to the coefficients
cCoef[0][j] = cCoef[0][j].subtract(coef3);
cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dHjsdAlpha));
cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dHjsdBeta));
cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(hjs));
sCoef[0][j] = sCoef[0][j].add(coef4);
sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dGjsdBeta));
sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(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 T lns = lnsCoef.getLns(n, jms);
final T dlns = lnsCoef.getdLnsdGamma(n, jms);
final T ijs = ghijCoef.getIjs(s, jms);
final T dIjsdh = ghijCoef.getdIjsdh(s, jms);
final T dIjsdk = ghijCoef.getdIjsdk(s, jms);
final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final T dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final T jjs = ghijCoef.getJjs(s, jms);
final T dJjsdh = ghijCoef.getdJjsdh(s, jms);
final T dJjsdk = ghijCoef.getdJjsdk(s, jms);
final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final T dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));
// K₀<sup>-n-1,s</sup>
final T kns = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final T dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
final T coef0 = jn.multiply(d0smj);
final T coef1 = coef0.multiply(lns);
final T coef2 = coef1.multiply(kns);
final T coef3 = coef2.multiply(ijs);
final T coef4 = coef2.multiply(jjs);
// Add the term to the coefficients
cCoef[0][j] = cCoef[0][j].subtract(coef3);
cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));
sCoef[0][j] = sCoef[0][j].add(coef4);
sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
}
}
}
}
if (isBetween(j, 2, nMax)) {
//add first term
// J<sub>j</sub>
final T jj = zero.subtract(harmonics.getUnnormalizedCnm(j, 0));
T kns = hansenObjects.getHansenObjects()[0].getValue(-j - 1, context.getChi());
T dkns = hansenObjects.getHansenObjects()[0].getDerivative(-j - 1, context.getChi());
T lns = lnsCoef.getLns(j, j);
//dlns is 0 because n == s == j
final T hjs = ghijCoef.getHjs(0, j);
final T dHjsdh = ghijCoef.getdHjsdh(0, j);
final T dHjsdk = ghijCoef.getdHjsdk(0, j);
final T dHjsdAlpha = ghijCoef.getdHjsdAlpha(0, j);
final T dHjsdBeta = ghijCoef.getdHjsdBeta(0, j);
final T gjs = ghijCoef.getGjs(0, j);
final T dGjsdh = ghijCoef.getdGjsdh(0, j);
final T dGjsdk = ghijCoef.getdGjsdk(0, j);
final T dGjsdAlpha = ghijCoef.getdGjsdAlpha(0, j);
final T dGjsdBeta = ghijCoef.getdGjsdBeta(0, j);
// 2 * J<sub>j</sub> * K₀<sup>-j-1,0</sup> * L<sub>j</sub><sup>j</sup>
T coef0 = jj.multiply(2.);
T coef1 = coef0.multiply(lns);
T coef2 = coef1.multiply(kns);
T coef3 = coef2.multiply(hjs);
T coef4 = coef2.multiply(gjs);
// Add the term to the coefficients
cCoef[0][j] = cCoef[0][j].subtract(coef3);
cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(j + 1));
cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dHjsdh).add(hjs.multiply(hXXX).multiply(dkns))));
cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dHjsdk).add(hjs.multiply(kXXX).multiply(dkns))));
cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dHjsdAlpha));
cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dHjsdBeta));
//no contribution to cCoef[6][j] because dlns is 0
sCoef[0][j] = sCoef[0][j].add(coef4);
sCoef[1][j] = sCoef[1][j].add(coef4.multiply(j + 1));
sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dGjsdh).add(gjs.multiply(hXXX).multiply(dkns))));
sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dGjsdk).add(gjs.multiply(kXXX).multiply(dkns))));
sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dGjsdAlpha));
sCoef[5][j] = sCoef[5][j].add(coef2.multiply(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.getHansenObjects()[s].getValue(-j - 1, context.getChi());
dkns = hansenObjects.getHansenObjects()[s].getDerivative(-j - 1, context.getChi());
lns = lnsCoef.getLns(j, jms);
final T dlns = lnsCoef.getdLnsdGamma(j, jms);
final T ijs = ghijCoef.getIjs(s, jms);
final T dIjsdh = ghijCoef.getdIjsdh(s, jms);
final T dIjsdk = ghijCoef.getdIjsdk(s, jms);
final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final T dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final T jjs = ghijCoef.getJjs(s, jms);
final T dJjsdh = ghijCoef.getdJjsdh(s, jms);
final T dJjsdk = ghijCoef.getdJjsdk(s, jms);
final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final T dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
coef0 = jj.multiply(d0smj);
coef1 = coef0.multiply(lns);
coef2 = coef1.multiply(kns);
coef3 = coef2.multiply(ijs);
coef4 = coef2.multiply(jjs);
// Add the term to the coefficients
cCoef[0][j] = cCoef[0][j].subtract(coef3);
cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(j + 1));
cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));
sCoef[0][j] = sCoef[0][j].add(coef4);
sCoef[1][j] = sCoef[1][j].add(coef4.multiply(j + 1));
sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(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 T lns = lnsCoef.getLns(n, jms);
final T dlns = lnsCoef.getdLnsdGamma(n, jms);
final T ijs = ghijCoef.getIjs(s, jms);
final T dIjsdh = ghijCoef.getdIjsdh(s, jms);
final T dIjsdk = ghijCoef.getdIjsdk(s, jms);
final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final T dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final T jjs = ghijCoef.getJjs(s, jms);
final T dJjsdh = ghijCoef.getdJjsdh(s, jms);
final T dJjsdk = ghijCoef.getdJjsdk(s, jms);
final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final T dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));
// K₀<sup>-n-1,s</sup>
final T kns = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final T dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
final T coef0 = jn.multiply(d0smj);
final T coef1 = coef0.multiply(lns);
final T coef2 = coef1.multiply(kns);
final T coef3 = coef2.multiply(ijs);
final T coef4 = coef2.multiply(jjs);
// Add the term to the coefficients
cCoef[0][j] = cCoef[0][j].subtract(coef3);
cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));
sCoef[0][j] = sCoef[0][j].add(coef4);
sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(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 T lns = lnsCoef.getLns(n, jms);
final T dlns = lnsCoef.getdLnsdGamma(n, jms);
final T ijs = ghijCoef.getIjs(s, jms);
final T dIjsdh = ghijCoef.getdIjsdh(s, jms);
final T dIjsdk = ghijCoef.getdIjsdk(s, jms);
final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final T dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final T jjs = ghijCoef.getJjs(s, jms);
final T dJjsdh = ghijCoef.getdJjsdh(s, jms);
final T dJjsdk = ghijCoef.getdJjsdk(s, jms);
final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final T dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));
// K₀<sup>-n-1,s</sup>
final T kns = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final T dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
final T coef0 = jn.multiply(d0smj);
final T coef1 = coef0.multiply(lns);
final T coef2 = coef1.multiply(kns);
final T coef3 = coef2.multiply(ijs);
final T coef4 = coef2.multiply(jjs);
// Add the term to the coefficients
cCoef[0][j] = cCoef[0][j].subtract(coef3);
cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));
sCoef[0][j] = sCoef[0][j].add(coef4);
sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(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 T lns = lnsCoef.getLns(n, jms);
final T dlns = lnsCoef.getdLnsdGamma(n, jms);
final T ijs = ghijCoef.getIjs(s, jms);
final T dIjsdh = ghijCoef.getdIjsdh(s, jms);
final T dIjsdk = ghijCoef.getdIjsdk(s, jms);
final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final T dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final T jjs = ghijCoef.getJjs(s, jms);
final T dJjsdh = ghijCoef.getdJjsdh(s, jms);
final T dJjsdk = ghijCoef.getdJjsdk(s, jms);
final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final T dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));
// K₀<sup>-n-1,s</sup>
final T kns = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final T dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
final T coef0 = jn.multiply(d0smj);
final T coef1 = coef0.multiply(lns);
final T coef2 = coef1.multiply(kns);
final T coef3 = coef2.multiply(ijs);
final T coef4 = coef2.multiply(jjs);
// Add the term to the coefficients
cCoef[0][j] = cCoef[0][j].subtract(coef3);
cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));
sCoef[0][j] = sCoef[0][j].add(coef4);
sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(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 T lns = lnsCoef.getLns(n, jms);
final T dlns = lnsCoef.getdLnsdGamma(n, jms);
final T ijs = ghijCoef.getIjs(s, jms);
final T dIjsdh = ghijCoef.getdIjsdh(s, jms);
final T dIjsdk = ghijCoef.getdIjsdk(s, jms);
final T dIjsdAlpha = ghijCoef.getdIjsdAlpha(s, jms);
final T dIjsdBeta = ghijCoef.getdIjsdBeta(s, jms);
final T jjs = ghijCoef.getJjs(s, jms);
final T dJjsdh = ghijCoef.getdJjsdh(s, jms);
final T dJjsdk = ghijCoef.getdJjsdk(s, jms);
final T dJjsdAlpha = ghijCoef.getdJjsdAlpha(s, jms);
final T dJjsdBeta = ghijCoef.getdJjsdBeta(s, jms);
// J<sub>n</sub>
final T jn = zero.subtract(harmonics.getUnnormalizedCnm(n, 0));
// K₀<sup>-n-1,s</sup>
final T kns = hansenObjects.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final T dkns = hansenObjects.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
final T coef0 = jn.multiply(d0smj);
final T coef1 = coef0.multiply(lns);
final T coef2 = coef1.multiply(kns);
final T coef3 = coef2.multiply(ijs);
final T coef4 = coef2.multiply(jjs);
// Add the term to the coefficients
cCoef[0][j] = cCoef[0][j].subtract(coef3);
cCoef[1][j] = cCoef[1][j].subtract(coef3.multiply(n + 1));
cCoef[2][j] = cCoef[2][j].subtract(coef1.multiply(kns.multiply(dIjsdh).add(ijs.multiply(hXXX).multiply(dkns))));
cCoef[3][j] = cCoef[3][j].subtract(coef1.multiply(kns.multiply(dIjsdk).add(ijs.multiply(kXXX).multiply(dkns))));
cCoef[4][j] = cCoef[4][j].subtract(coef2.multiply(dIjsdAlpha));
cCoef[5][j] = cCoef[5][j].subtract(coef2.multiply(dIjsdBeta));
cCoef[6][j] = cCoef[6][j].subtract(coef0.multiply(dlns).multiply(kns).multiply(ijs));
sCoef[0][j] = sCoef[0][j].add(coef4);
sCoef[1][j] = sCoef[1][j].add(coef4.multiply(n + 1));
sCoef[2][j] = sCoef[2][j].add(coef1.multiply(kns.multiply(dJjsdh).add(jjs.multiply(hXXX).multiply(dkns))));
sCoef[3][j] = sCoef[3][j].add(coef1.multiply(kns.multiply(dJjsdk).add(jjs.multiply(kXXX).multiply(dkns))));
sCoef[4][j] = sCoef[4][j].add(coef2.multiply(dJjsdAlpha));
sCoef[5][j] = sCoef[5][j].add(coef2.multiply(dJjsdBeta));
sCoef[6][j] = sCoef[6][j].add(coef0.multiply(dlns).multiply(kns).multiply(jjs));
}
}
}
}
}
}
cCoef[0][j] = cCoef[0][j].multiply(context.getMuoa().divide(j).negate());
cCoef[1][j] = cCoef[1][j].multiply(context.getMuoa().divide(auxiliaryElements.getSma().multiply(j)));
cCoef[2][j] = cCoef[2][j].multiply(context.getMuoa().divide(j).negate());
cCoef[3][j] = cCoef[3][j].multiply(context.getMuoa().divide(j).negate());
cCoef[4][j] = cCoef[4][j].multiply(context.getMuoa().divide(j).negate());
cCoef[5][j] = cCoef[5][j].multiply(context.getMuoa().divide(j).negate());
cCoef[6][j] = cCoef[6][j].multiply(context.getMuoa().divide(j).negate());
sCoef[0][j] = sCoef[0][j].multiply(context.getMuoa().divide(j).negate());
sCoef[1][j] = sCoef[1][j].multiply(context.getMuoa().divide(auxiliaryElements.getSma().multiply(j)));
sCoef[2][j] = sCoef[2][j].multiply(context.getMuoa().divide(j).negate());
sCoef[3][j] = sCoef[3][j].multiply(context.getMuoa().divide(j).negate());
sCoef[4][j] = sCoef[4][j].multiply(context.getMuoa().divide(j).negate());
sCoef[5][j] = sCoef[5][j].multiply(context.getMuoa().divide(j).negate());
sCoef[6][j] = sCoef[6][j].multiply(context.getMuoa().divide(j).negate());
}
}
/** 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 T 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 T 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 T 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 T 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 T 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 T 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 T 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 T 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 T 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 T 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 T 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 T 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 T 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 T getdSjdGamma(final int j) {
return sCoef[6][j];
}
}
/** Coefficients valid for one time slot. */
private static class Slot {
/**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);
}
}
}
/** Coefficients valid for one time slot. */
private static class FieldSlot <T extends CalculusFieldElement<T>> {
/**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 FieldShortPeriodicsInterpolatedCoefficient<T> 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 FieldShortPeriodicsInterpolatedCoefficient<T>[] 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 FieldShortPeriodicsInterpolatedCoefficient<T>[] sij;
/** Simple constructor.
* @param maxFrequencyShortPeriodics maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
*/
@SuppressWarnings("unchecked")
FieldSlot(final int maxFrequencyShortPeriodics, final int interpolationPoints) {
final int rows = maxFrequencyShortPeriodics + 1;
di = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows);
sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, rows);
//Initialize the arrays
for (int j = 0; j <= maxFrequencyShortPeriodics; j++) {
cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
}
}
}
/** Compute potential and potential derivatives with respect to orbital parameters. */
private class UAnddU {
/** The current value of the U function. <br/>
* Needed for the short periodic contribution */
private double U;
/** dU / da. */
private double dUda;
/** dU / dk. */
private double dUdk;
/** dU / dh. */
private double dUdh;
/** dU / dAlpha. */
private double dUdAl;
/** dU / dBeta. */
private double dUdBe;
/** dU / dGamma. */
private double dUdGa;
/** Simple constuctor.
* @param date current date
* @param context container for attributes
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param hansen initialization of hansen objects
*/
UAnddU(final AbsoluteDate date,
final DSSTZonalContext context,
final AuxiliaryElements auxiliaryElements,
final HansenObjects hansen) {
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
//Reset U
U = 0.;
// Gs and Hs coefficients
final double[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), context.getAlpha(), context.getBeta(), maxEccPowMeanElements);
// Qns coefficients
final double[][] Qns = CoefficientsFactory.computeQns(context.getGamma(), maxDegree, maxEccPowMeanElements);
final double[] roaPow = new double[maxDegree + 1];
roaPow[0] = 1.;
for (int i = 1; i <= maxDegree; i++) {
roaPow[i] = context.getRoa() * roaPow[i - 1];
}
// Potential derivatives
dUda = 0.;
dUdk = 0.;
dUdh = 0.;
dUdAl = 0.;
dUdBe = 0.;
dUdGa = 0.;
for (int s = 0; s <= maxEccPowMeanElements; s++) {
//Initialize the Hansen roots
hansen.computeHansenObjectsInitValues(context, s);
// 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 = context.getBeta() * sxgsm1 - context.getAlpha() * sxhsm1;
dGsdk = context.getAlpha() * sxgsm1 + context.getBeta() * sxhsm1;
dGsdAl = auxiliaryElements.getK() * sxgsm1 - auxiliaryElements.getH() * sxhsm1;
dGsdBe = auxiliaryElements.getH() * sxgsm1 + auxiliaryElements.getK() * 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 = hansen.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final double dkns = hansen.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
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 + auxiliaryElements.getK() * context.getChi3() * gs * dkns);
// Compute dU / dEy
dUdh += coef1 * (kns * dGsdh + auxiliaryElements.getH() * context.getChi3() * 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)
this.U = -context.getMuoa() * U;
this.dUda = dUda * context.getMuoa() / auxiliaryElements.getSma();
this.dUdk = dUdk * -context.getMuoa();
this.dUdh = dUdh * -context.getMuoa();
this.dUdAl = dUdAl * -context.getMuoa();
this.dUdBe = dUdBe * -context.getMuoa();
this.dUdGa = dUdGa * -context.getMuoa();
}
/** Return value of U.
* @return U
*/
public double getU() {
return U;
}
/** Return value of dU / da.
* @return dUda
*/
public double getdUda() {
return dUda;
}
/** Return value of dU / dk.
* @return dUdk
*/
public double getdUdk() {
return dUdk;
}
/** Return value of dU / dh.
* @return dUdh
*/
public double getdUdh() {
return dUdh;
}
/** Return value of dU / dAlpha.
* @return dUdAl
*/
public double getdUdAl() {
return dUdAl;
}
/** Return value of dU / dBeta.
* @return dUdBe
*/
public double getdUdBe() {
return dUdBe;
}
/** Return value of dU / dGamma.
* @return dUdGa
*/
public double getdUdGa() {
return dUdGa;
}
}
/** 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>
*/
private class FieldUAnddU <T extends CalculusFieldElement<T>> {
/** The current value of the U function. <br/>
* Needed for the short periodic contribution */
private T U;
/** dU / da. */
private T dUda;
/** dU / dk. */
private T dUdk;
/** dU / dh. */
private T dUdh;
/** dU / dAlpha. */
private T dUdAl;
/** dU / dBeta. */
private T dUdBe;
/** dU / dGamma. */
private T dUdGa;
/** Simple constuctor.
* @param date current date
* @param context container for attributes
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param hansen initialization of hansen objects
*/
FieldUAnddU(final FieldAbsoluteDate<T> date,
final FieldDSSTZonalContext<T> context,
final FieldAuxiliaryElements<T> auxiliaryElements,
final FieldHansenObjects<T> hansen) {
// Zero for initialization
final Field<T> field = date.getField();
final T zero = field.getZero();
//spherical harmonics
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date.toAbsoluteDate());
//Reset U
U = zero;
// Gs and Hs coefficients
final T[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(),
context.getAlpha(), context.getBeta(),
maxEccPowMeanElements, field);
// Qns coefficients
final T[][] Qns = CoefficientsFactory.computeQns(context.getGamma(), maxDegree, maxEccPowMeanElements);
final T[] roaPow = MathArrays.buildArray(field, maxDegree + 1);
roaPow[0] = zero.newInstance(1.);
for (int i = 1; i <= maxDegree; i++) {
roaPow[i] = roaPow[i - 1].multiply(context.getRoa());
}
// Potential derivatives
dUda = zero;
dUdk = zero;
dUdh = zero;
dUdAl = zero;
dUdBe = zero;
dUdGa = zero;
for (int s = 0; s <= maxEccPowMeanElements; s++) {
//Initialize the Hansen roots
hansen.computeHansenObjectsInitValues(context, s);
// Get the current Gs coefficient
final T gs = GsHs[0][s];
// Compute Gs partial derivatives from 3.1-(9)
T dGsdh = zero;
T dGsdk = zero;
T dGsdAl = zero;
T dGsdBe = zero;
if (s > 0) {
// First get the G(s-1) and the H(s-1) coefficients
final T sxgsm1 = GsHs[0][s - 1].multiply(s);
final T sxhsm1 = GsHs[1][s - 1].multiply(s);
// Then compute derivatives
dGsdh = sxgsm1.multiply(context.getBeta()).subtract(sxhsm1.multiply(context.getAlpha()));
dGsdk = sxgsm1.multiply(context.getAlpha()).add(sxhsm1.multiply(context.getBeta()));
dGsdAl = sxgsm1.multiply(auxiliaryElements.getK()).subtract(sxhsm1.multiply(auxiliaryElements.getH()));
dGsdBe = sxgsm1.multiply(auxiliaryElements.getH()).add(sxhsm1.multiply(auxiliaryElements.getK()));
}
// Kronecker symbol (2 - delta(0,s))
final T d0s = zero.newInstance((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 T kns = hansen.getHansenObjects()[s].getValue(-n - 1, context.getChi());
final T dkns = hansen.getHansenObjects()[s].getDerivative(-n - 1, context.getChi());
final double vns = Vns.get(new NSKey(n, s));
final T coef0 = d0s.multiply(roaPow[n]).multiply(vns).multiply(-harmonics.getUnnormalizedCnm(n, 0));
final T coef1 = coef0.multiply(Qns[n][s]);
final T coef2 = coef1.multiply(kns);
final T coef3 = coef2.multiply(gs);
// dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
final T dqns = Qns[n][s + 1];
// Compute U
U = U.add(coef3);
// Compute dU / da :
dUda = dUda.add(coef3.multiply(n + 1));
// Compute dU / dEx
dUdk = dUdk.add(coef1.multiply(dGsdk.multiply(kns).add(auxiliaryElements.getK().multiply(context.getChi3()).multiply(dkns).multiply(gs))));
// Compute dU / dEy
dUdh = dUdh.add(coef1.multiply(dGsdh.multiply(kns).add(auxiliaryElements.getH().multiply(context.getChi3()).multiply(dkns).multiply(gs))));
// Compute dU / dAlpha
dUdAl = dUdAl.add(coef2.multiply(dGsdAl));
// Compute dU / dBeta
dUdBe = dUdBe.add(coef2.multiply(dGsdBe));
// Compute dU / dGamma
dUdGa = dUdGa.add(coef0.multiply(kns).multiply(dqns).multiply(gs));
}
}
}
// Multiply by -(μ / a)
U = U.multiply(context.getMuoa().negate());
dUda = dUda.multiply(context.getMuoa().divide(auxiliaryElements.getSma()));
dUdk = dUdk.multiply(context.getMuoa()).negate();
dUdh = dUdh.multiply(context.getMuoa()).negate();
dUdAl = dUdAl.multiply(context.getMuoa()).negate();
dUdBe = dUdBe.multiply(context.getMuoa()).negate();
dUdGa = dUdGa.multiply(context.getMuoa()).negate();
}
/** Return value of U.
* @return U
*/
public T getU() {
return U;
}
/** Return value of dU / da.
* @return dUda
*/
public T getdUda() {
return dUda;
}
/** Return value of dU / dk.
* @return dUdk
*/
public T getdUdk() {
return dUdk;
}
/** Return value of dU / dh.
* @return dUdh
*/
public T getdUdh() {
return dUdh;
}
/** Return value of dU / dAlpha.
* @return dUdAl
*/
public T getdUdAl() {
return dUdAl;
}
/** Return value of dU / dBeta.
* @return dUdBe
*/
public T getdUdBe() {
return dUdBe;
}
/** Return value of dU / dGamma.
* @return dUdGa
*/
public T getdUdGa() {
return dUdGa;
}
}
/** Computes init values of the Hansen Objects. */
private class HansenObjects {
/** An array that contains the objects needed to build the Hansen coefficients. <br/>
* The index is s*/
private HansenZonalLinear[] hansenObjects;
/** Simple constructor. */
HansenObjects() {
this.hansenObjects = new HansenZonalLinear[maxEccPow + 1];
for (int s = 0; s <= maxEccPow; s++) {
this.hansenObjects[s] = new HansenZonalLinear(maxDegree, s);
}
}
/** Compute init values for hansen objects.
* @param context container for attributes
* @param element element of the array to compute the init values
*/
public void computeHansenObjectsInitValues(final DSSTZonalContext context, final int element) {
hansenObjects[element].computeInitValues(context.getChi());
}
/** Get the Hansen Objects.
* @return hansenObjects
*/
public HansenZonalLinear[] getHansenObjects() {
return hansenObjects;
}
}
/** Computes init values of the Hansen Objects.
* @param <T> type of the elements
*/
private class FieldHansenObjects<T extends CalculusFieldElement<T>> {
/** An array that contains the objects needed to build the Hansen coefficients. <br/>
* The index is s*/
private FieldHansenZonalLinear<T>[] hansenObjects;
/** Simple constructor.
* @param field field used by default
*/
@SuppressWarnings("unchecked")
FieldHansenObjects(final Field<T> field) {
this.hansenObjects = (FieldHansenZonalLinear<T>[]) Array.newInstance(FieldHansenZonalLinear.class, maxEccPow + 1);
for (int s = 0; s <= maxEccPow; s++) {
this.hansenObjects[s] = new FieldHansenZonalLinear<>(maxDegree, s, field);
}
}
/** Compute init values for hansen objects.
* @param context container for attributes
* @param element element of the array to compute the init values
*/
public void computeHansenObjectsInitValues(final FieldDSSTZonalContext<T> context, final int element) {
hansenObjects[element].computeInitValues(context.getChi());
}
/** Get the Hansen Objects.
* @return hansenObjects
*/
public FieldHansenZonalLinear<T>[] getHansenObjects() {
return hansenObjects;
}
}
}