AbstractGaussianContribution.java
/* Copyright 2002-2019 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation.semianalytical.dsst.forces;
import java.lang.reflect.Array;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
import org.hipparchus.analysis.RealFieldUnivariateVectorFunction;
import org.hipparchus.analysis.UnivariateVectorFunction;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;
import org.orekit.attitudes.Attitude;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.FieldAttitude;
import org.orekit.forces.ForceModel;
import org.orekit.orbits.EquinoctialOrbit;
import org.orekit.orbits.FieldEquinoctialOrbit;
import org.orekit.orbits.FieldOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
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.FieldAuxiliaryElements;
import org.orekit.propagation.semianalytical.dsst.utilities.FieldCjSjCoefficient;
import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
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;
/** Common handling of {@link DSSTForceModel} methods for Gaussian contributions to DSST propagation.
* <p>
* This abstract class allows to provide easily a subset of {@link DSSTForceModel} methods
* for specific Gaussian contributions.
* </p><p>
* This class implements the notion of numerical averaging of the DSST theory.
* Numerical averaging is mainly used for non-conservative disturbing forces such as
* atmospheric drag and solar radiation pressure.
* </p><p>
* Gaussian contributions can be expressed as: da<sub>i</sub>/dt = δa<sub>i</sub>/δv . q<br>
* where:
* <ul>
* <li>a<sub>i</sub> are the six equinoctial elements</li>
* <li>v is the velocity vector</li>
* <li>q is the perturbing acceleration due to the considered force</li>
* </ul>
*
* <p> The averaging process and other considerations lead to integrate this contribution
* over the true longitude L possibly taking into account some limits.
*
* <p> To create a numerically averaged contribution, one needs only to provide a
* {@link ForceModel} and to implement in the derived class the methods:
* {@link #getLLimits(SpacecraftState, AuxiliaryElements)}
* and {@link #getParametersDriversWithoutMu()}.
* </p>
* @author Pascal Parraud
* @author Bryan Cazabonne (field translation)
*/
public abstract class AbstractGaussianContribution implements DSSTForceModel {
/** 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;
/** 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);
/** Available orders for Gauss quadrature. */
private static final int[] GAUSS_ORDER = {12, 16, 20, 24, 32, 40, 48};
/** Max rank in Gauss quadrature orders array. */
private static final int MAX_ORDER_RANK = GAUSS_ORDER.length - 1;
/** Number of points for interpolation. */
private static final int INTERPOLATION_POINTS = 3;
/** Maximum value for j index. */
private static final int JMAX = 12;
/** Contribution to be numerically averaged. */
private final ForceModel contribution;
/** Gauss integrator. */
private final double threshold;
/** Gauss integrator. */
private GaussQuadrature integrator;
/** Flag for Gauss order computation. */
private boolean isDirty;
/** Attitude provider. */
private AttitudeProvider attitudeProvider;
/** Prefix for coefficients keys. */
private final String coefficientsKeyPrefix;
/** Short period terms. */
private GaussianShortPeriodicCoefficients gaussianSPCoefs;
/** Short period terms. */
private Map<Field<?>, FieldGaussianShortPeriodicCoefficients<?>> gaussianFieldSPCoefs;
/** Driver for gravitational parameter. */
private final ParameterDriver gmParameterDriver;
/** Build a new instance.
* @param coefficientsKeyPrefix prefix for coefficients keys
* @param threshold tolerance for the choice of the Gauss quadrature order
* @param contribution the {@link ForceModel} to be numerically averaged
* @param mu central attraction coefficient
*/
protected AbstractGaussianContribution(final String coefficientsKeyPrefix,
final double threshold,
final ForceModel contribution,
final double mu) {
gmParameterDriver = new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
mu, MU_SCALE,
0.0, Double.POSITIVE_INFINITY);
this.coefficientsKeyPrefix = coefficientsKeyPrefix;
this.contribution = contribution;
this.threshold = threshold;
this.integrator = new GaussQuadrature(GAUSS_ORDER[MAX_ORDER_RANK]);
this.isDirty = true;
gaussianFieldSPCoefs = new HashMap<>();
}
/** {@inheritDoc} */
@Override
public ParameterDriver[] getParametersDrivers() {
final ParameterDriver[] driversWithoutMu = getParametersDriversWithoutMu();
// + 1 for central attraction coefficient driver
final ParameterDriver[] drivers = new ParameterDriver[driversWithoutMu.length + 1];
int index = 0;
for (final ParameterDriver driver : driversWithoutMu) {
drivers[index] = driver;
index++;
}
// We put central attraction coefficient driver at the end of the array
drivers[driversWithoutMu.length] = gmParameterDriver;
return drivers;
}
/** Get the drivers for force model parameters except the one for the central attraction coefficient.
* <p>
* The driver for central attraction coefficient is automatically
* added at the last element of the {@link ParameterDriver} array
* into {@link #getParametersDrivers()} method.
* </p>
* @return drivers for force model parameters
*/
protected abstract ParameterDriver[] getParametersDriversWithoutMu();
/** {@inheritDoc} */
@Override
public List<ShortPeriodTerms> initialize(final AuxiliaryElements auxiliaryElements,
final PropagationType type,
final double[] parameters) {
final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
gaussianSPCoefs = new GaussianShortPeriodicCoefficients(coefficientsKeyPrefix,
JMAX, INTERPOLATION_POINTS,
new TimeSpanMap<Slot>(new Slot(JMAX, INTERPOLATION_POINTS)));
list.add(gaussianSPCoefs);
return list;
}
/** {@inheritDoc} */
@Override
public <T extends RealFieldElement<T>> List<FieldShortPeriodTerms<T>> initialize(final FieldAuxiliaryElements<T> auxiliaryElements,
final PropagationType type,
final T[] parameters) {
final Field<T> field = auxiliaryElements.getDate().getField();
final FieldGaussianShortPeriodicCoefficients<T> fgspc =
new FieldGaussianShortPeriodicCoefficients<>(coefficientsKeyPrefix,
JMAX, INTERPOLATION_POINTS,
new FieldTimeSpanMap<>(new FieldSlot<>(JMAX,
INTERPOLATION_POINTS),
field));
gaussianFieldSPCoefs.put(field, fgspc);
return Collections.singletonList(fgspc);
}
/** 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 parameters values of the force model parameters
* @return new force model context
*/
private AbstractGaussianContributionContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
return new AbstractGaussianContributionContext(auxiliaryElements, 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 parameters values of the force model parameters
* @return new force model context
*/
private <T extends RealFieldElement<T>> FieldAbstractGaussianContributionContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
final T[] parameters) {
return new FieldAbstractGaussianContributionContext<>(auxiliaryElements, parameters);
}
/** {@inheritDoc} */
@Override
public double[] getMeanElementRate(final SpacecraftState state,
final AuxiliaryElements auxiliaryElements, final double[] parameters) {
// Container for attributes
final AbstractGaussianContributionContext context = initializeStep(auxiliaryElements, parameters);
double[] meanElementRate = new double[6];
// Computes the limits for the integral
final double[] ll = getLLimits(state, auxiliaryElements);
// Computes integrated mean element rates if Llow < Lhigh
if (ll[0] < ll[1]) {
meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1], context, parameters);
if (isDirty) {
boolean next = true;
for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
final double[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0], ll[1], context, parameters);
if (getRatesDiff(meanElementRate, meanRates, context) < threshold) {
integrator = new GaussQuadrature(GAUSS_ORDER[i]);
next = false;
}
}
isDirty = false;
}
}
return meanElementRate;
}
/** {@inheritDoc} */
@Override
public <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> state,
final FieldAuxiliaryElements<T> auxiliaryElements,
final T[] parameters) {
// Container for attributes
final FieldAbstractGaussianContributionContext<T> context = initializeStep(auxiliaryElements, parameters);
final Field<T> field = state.getDate().getField();
T[] meanElementRate = MathArrays.buildArray(field, 6);
// Computes the limits for the integral
final T[] ll = getLLimits(state, auxiliaryElements);
// Computes integrated mean element rates if Llow < Lhigh
if (ll[0].getReal() < ll[1].getReal()) {
meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1], context, parameters);
if (isDirty) {
boolean next = true;
for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
final T[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0], ll[1], context, parameters);
if (getRatesDiff(meanElementRate, meanRates, context).getReal() < threshold) {
integrator = new GaussQuadrature(GAUSS_ORDER[i]);
next = false;
}
}
isDirty = false;
}
}
return meanElementRate;
}
/** Compute the limits in L, the true longitude, for integration.
*
* @param state current state information: date, kinematics, attitude
* @param auxiliaryElements auxiliary elements related to the current orbit
* @return the integration limits in L
*/
protected abstract double[] getLLimits(SpacecraftState state, AuxiliaryElements auxiliaryElements);
/** Compute the limits in L, the true longitude, for integration.
*
* @param <T> type of the elements
* @param state current state information: date, kinematics, attitude
* @param auxiliaryElements auxiliary elements related to the current orbit
* @return the integration limits in L
*/
protected abstract <T extends RealFieldElement<T>> T[] getLLimits(FieldSpacecraftState<T> state,
FieldAuxiliaryElements<T> auxiliaryElements);
/** Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
*
* @param state current state
* @param gauss Gauss quadrature
* @param low lower bound of the integral interval
* @param high upper bound of the integral interval
* @param context container for attributes
* @param parameters values of the force model parameters
* @return the mean element rates
*/
private double[] getMeanElementRate(final SpacecraftState state,
final GaussQuadrature gauss,
final double low,
final double high,
final AbstractGaussianContributionContext context,
final double[] parameters) {
// Auxiliary elements related to the current orbit
final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
final double[] meanElementRate = gauss.integrate(new IntegrableFunction(state, true, 0, parameters), low, high);
// Constant multiplier for integral
final double coef = 1. / (2. * FastMath.PI * auxiliaryElements.getB());
// Corrects mean element rates
for (int i = 0; i < 6; i++) {
meanElementRate[i] *= coef;
}
return meanElementRate;
}
/** Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
*
* @param <T> type of the elements
* @param state current state
* @param gauss Gauss quadrature
* @param low lower bound of the integral interval
* @param high upper bound of the integral interval
* @param context container for attributes
* @param parameters values of the force model parameters
* @return the mean element rates
*/
private <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> state,
final GaussQuadrature gauss,
final T low,
final T high,
final FieldAbstractGaussianContributionContext<T> context,
final T[] parameters) {
// Field
final Field<T> field = context.getA().getField();
// Auxiliary elements related to the current orbit
final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
final T[] meanElementRate = gauss.integrate(new FieldIntegrableFunction<>(state, true, 0, parameters, field), low, high, field);
// Constant multiplier for integral
final T coef = auxiliaryElements.getB().multiply(FastMath.PI).multiply(2.).reciprocal();
// Corrects mean element rates
for (int i = 0; i < 6; i++) {
meanElementRate[i] = meanElementRate[i].multiply(coef);
}
return meanElementRate;
}
/** Estimates the weighted magnitude of the difference between 2 sets of equinoctial elements rates.
*
* @param meanRef reference rates
* @param meanCur current rates
* @param context container for attributes
* @return estimated magnitude of weighted differences
*/
private double getRatesDiff(final double[] meanRef, final double[] meanCur, final AbstractGaussianContributionContext context) {
// Auxiliary elements related to the current orbit
final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
double maxDiff = FastMath.abs(meanRef[0] - meanCur[0]) / auxiliaryElements.getSma();
// Corrects mean element rates
for (int i = 1; i < meanRef.length; i++) {
final double diff = FastMath.abs(meanRef[i] - meanCur[i]);
if (maxDiff < diff) maxDiff = diff;
}
return maxDiff;
}
/** Estimates the weighted magnitude of the difference between 2 sets of equinoctial elements rates.
*
* @param <T> type of the elements
* @param meanRef reference rates
* @param meanCur current rates
* @param context container for attributes
* @return estimated magnitude of weighted differences
*/
private <T extends RealFieldElement<T>> T getRatesDiff(final T[] meanRef, final T[] meanCur, final FieldAbstractGaussianContributionContext<T> context) {
// Auxiliary elements related to the current orbit
final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
T maxDiff = FastMath.abs(meanRef[0].subtract(meanCur[0])).divide(auxiliaryElements.getSma());;
// Corrects mean element rates
for (int i = 1; i < meanRef.length; i++) {
final T diff = FastMath.abs(meanRef[i].subtract(meanCur[i]));
if (maxDiff.getReal() < diff.getReal()) maxDiff = diff;
}
return maxDiff;
}
/** {@inheritDoc} */
@Override
public void registerAttitudeProvider(final AttitudeProvider provider) {
this.attitudeProvider = provider;
}
/** {@inheritDoc} */
@Override
public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {
final Slot slot = gaussianSPCoefs.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
final AbstractGaussianContributionContext context = initializeStep(auxiliaryElements, parameters);
// Compute rhoj and sigmaj
final double[][] currentRhoSigmaj = computeRhoSigmaCoefficients(meanState.getDate(), auxiliaryElements);
// Generate the Cij and Sij coefficients
final FourierCjSjCoefficients fourierCjSj = new FourierCjSjCoefficients(meanState, JMAX, auxiliaryElements, parameters);
// Generate the Uij and Vij coefficients
final UijVijCoefficients uijvij = new UijVijCoefficients(currentRhoSigmaj, fourierCjSj, JMAX);
gaussianSPCoefs.computeCoefficients(meanState, slot, fourierCjSj, uijvij, context.getMeanMotion(), auxiliaryElements.getSma());
}
}
/** {@inheritDoc} */
@Override
@SuppressWarnings("unchecked")
public <T extends RealFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
final FieldSpacecraftState<T>... meanStates) {
// Field used by default
final Field<T> field = meanStates[0].getDate().getField();
final FieldGaussianShortPeriodicCoefficients<T> fgspc = (FieldGaussianShortPeriodicCoefficients<T>) gaussianFieldSPCoefs.get(field);
final FieldSlot<T> slot = fgspc.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
final FieldAbstractGaussianContributionContext<T> context = initializeStep(auxiliaryElements, parameters);
// Compute rhoj and sigmaj
final T[][] currentRhoSigmaj = computeRhoSigmaCoefficients(meanState.getDate(), context, field);
// Generate the Cij and Sij coefficients
final FieldFourierCjSjCoefficients<T> fourierCjSj = new FieldFourierCjSjCoefficients<>(meanState, JMAX, auxiliaryElements, parameters, field);
// Generate the Uij and Vij coefficients
final FieldUijVijCoefficients<T> uijvij = new FieldUijVijCoefficients<>(currentRhoSigmaj, fourierCjSj, JMAX, field);
fgspc.computeCoefficients(meanState, slot, fourierCjSj, uijvij, context.getMeanMotion(), auxiliaryElements.getSma(), field);
}
}
/**
* Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
* <p>
* The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
* ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
* σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
* </p>
* @param date current date
* @param auxiliaryElements auxiliary elements related to the current orbit
* @return computed coefficients
*/
private double[][] computeRhoSigmaCoefficients(final AbsoluteDate date, final AuxiliaryElements auxiliaryElements) {
final double[][] currentRhoSigmaj = new double[2][3 * JMAX + 1];
final CjSjCoefficient cjsjKH = new CjSjCoefficient(auxiliaryElements.getK(), auxiliaryElements.getH());
final double b = 1. / (1 + auxiliaryElements.getB());
// (-b)<sup>j</sup>
double mbtj = 1;
for (int j = 1; j <= 3 * JMAX; j++) {
//Compute current rho and sigma;
mbtj *= -b;
final double coef = (1 + j * auxiliaryElements.getB()) * mbtj;
currentRhoSigmaj[0][j] = coef * cjsjKH.getCj(j);
currentRhoSigmaj[1][j] = coef * cjsjKH.getSj(j);
}
return currentRhoSigmaj;
}
/**
* 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 date current date
* @param context container for attributes
* @param field field used by default
* @return computed coefficients
*/
private <T extends RealFieldElement<T>> T[][] computeRhoSigmaCoefficients(final FieldAbsoluteDate<T> date,
final FieldAbstractGaussianContributionContext<T> context,
final Field<T> field) {
// zero
final T zero = field.getZero();
final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
final T[][] currentRhoSigmaj = MathArrays.buildArray(field, 2, 3 * JMAX + 1);
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.add(1.);
for (int j = 1; j <= 3 * JMAX; j++) {
//Compute current rho and sigma;
mbtj = mbtj.multiply(b.negate());
final T coef = mbtj.multiply(auxiliaryElements.getB().multiply(j).add(1.));
currentRhoSigmaj[0][j] = coef.multiply(cjsjKH.getCj(j));
currentRhoSigmaj[1][j] = coef.multiply(cjsjKH.getSj(j));
}
return currentRhoSigmaj;
}
/** Internal class for numerical quadrature.
* <p>
* This class is a rewrite of {@link IntegrableFunction} for field elements
* </p>
*/
private class FieldIntegrableFunction <T extends RealFieldElement<T>> implements RealFieldUnivariateVectorFunction<T> {
/** Current state. */
private final FieldSpacecraftState<T> state;
/** Signal that this class is used to compute the values required by the mean element variations
* or by the short periodic element variations. */
private final boolean meanMode;
/** The j index.
* <p>
* Used only for short periodic variation. Ignored for mean elements variation.
* </p> */
private final int j;
/** Container for attributes. */
private final FieldAbstractGaussianContributionContext<T> context;
/** Auxiliary Elements. */
private final FieldAuxiliaryElements<T> auxiliaryElements;
/** Drivers for solar radiation and atmospheric drag forces. */
private final T[] parameters;
/** Build a new instance with a new field.
* @param state current state information: date, kinematics, attitude
* @param meanMode if true return the value associated to the mean elements variation,
* if false return the values associated to the short periodic elements variation
* @param j the j index. used only for short periodic variation. Ignored for mean elements variation.
* @param parameters values of the force model parameters
* @param field field utilized by default
*/
FieldIntegrableFunction(final FieldSpacecraftState<T> state,
final boolean meanMode,
final int j,
final T[] parameters,
final Field<T> field) {
this.meanMode = meanMode;
this.j = j;
this.parameters = parameters;
this.auxiliaryElements = new FieldAuxiliaryElements<>(state.getOrbit(), I);
this.context = new FieldAbstractGaussianContributionContext<>(auxiliaryElements, parameters);
// remove derivatives from state
final T[] stateVector = MathArrays.buildArray(field, 6);
OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngle.TRUE, stateVector, null);
final FieldOrbit<T> fixedOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(stateVector, null, PositionAngle.TRUE,
state.getDate(),
context.getMu(),
state.getFrame());
this.state = new FieldSpacecraftState<>(fixedOrbit, state.getAttitude(), state.getMass());
}
/** {@inheritDoc} */
@Override
public T[] value(final T x) {
// Parameters for array building
final Field<T> field = auxiliaryElements.getDate().getField();
final int dimension = 6;
//Compute the time difference from the true longitude difference
final T shiftedLm = trueToMean(x);
final T dLm = shiftedLm.subtract(auxiliaryElements.getLM());
final T dt = dLm.divide(context.getMeanMotion());
final T cosL = FastMath.cos(x);
final T sinL = FastMath.sin(x);
final T roa = auxiliaryElements.getB().multiply(auxiliaryElements.getB()).divide(auxiliaryElements.getH().multiply(sinL).add(auxiliaryElements.getK().multiply(cosL)).add(1.));
final T roa2 = roa.multiply(roa);
final T r = auxiliaryElements.getSma().multiply(roa);
final T X = r.multiply(cosL);
final T Y = r.multiply(sinL);
final T naob = context.getMeanMotion().multiply(auxiliaryElements.getSma()).divide(auxiliaryElements.getB());
final T Xdot = naob.multiply(auxiliaryElements.getH().add(sinL)).negate();
final T Ydot = naob.multiply(auxiliaryElements.getK().add(cosL));
final FieldVector3D<T> vel = new FieldVector3D<>(Xdot, auxiliaryElements.getVectorF(), Ydot, auxiliaryElements.getVectorG());
// Compute acceleration
FieldVector3D<T> acc = FieldVector3D.getZero(field);
// shift the orbit to dt
final FieldOrbit<T> shiftedOrbit = state.getOrbit().shiftedBy(dt);
// Recompose an orbit with time held fixed to be compliant with DSST theory
final FieldOrbit<T> recomposedOrbit =
new FieldEquinoctialOrbit<>(shiftedOrbit.getA(),
shiftedOrbit.getEquinoctialEx(),
shiftedOrbit.getEquinoctialEy(),
shiftedOrbit.getHx(),
shiftedOrbit.getHy(),
shiftedOrbit.getLv(),
PositionAngle.TRUE,
shiftedOrbit.getFrame(),
state.getDate(),
context.getMu());
// Get the corresponding attitude
final FieldAttitude<T> recomposedAttitude =
attitudeProvider.getAttitude(recomposedOrbit,
recomposedOrbit.getDate(),
recomposedOrbit.getFrame());
// create shifted SpacecraftState with attitude at specified time
final FieldSpacecraftState<T> shiftedState =
new FieldSpacecraftState<>(recomposedOrbit, recomposedAttitude, state.getMass());
acc = contribution.acceleration(shiftedState, parameters);
//Compute the derivatives of the elements by the speed
final T[] deriv = MathArrays.buildArray(field, dimension);
// da/dv
deriv[0] = getAoV(vel).dotProduct(acc);
// dex/dv
deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
// dey/dv
deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
// dhx/dv
deriv[3] = getQoV(X).dotProduct(acc);
// dhy/dv
deriv[4] = getPoV(Y).dotProduct(acc);
// dλ/dv
deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);
// Compute mean elements rates
T[] val = null;
if (meanMode) {
val = MathArrays.buildArray(field, dimension);
for (int i = 0; i < 6; i++) {
// da<sub>i</sub>/dt
val[i] = deriv[i].multiply(roa2);
}
} else {
val = MathArrays.buildArray(field, dimension * 2);
//Compute cos(j*L) and sin(j*L);
final T cosjL = j == 1 ? cosL : FastMath.cos(x.multiply(j));
final T sinjL = j == 1 ? sinL : FastMath.sin(x.multiply(j));
for (int i = 0; i < 6; i++) {
// da<sub>i</sub>/dv * cos(jL)
val[i] = deriv[i].multiply(cosjL);
// da<sub>i</sub>/dv * sin(jL)
val[i + 6] = deriv[i].multiply(sinjL);
}
}
return val;
}
/** Converts true longitude to mean longitude.
* @param x True longitude
* @return Eccentric longitude
*/
private T trueToMean (final T x) {
return eccentricToMean(trueToEccentric(x));
}
/** Converts true longitude to eccentric longitude.
* @param lv True longitude
* @return Eccentric longitude
*/
private T trueToEccentric (final T lv) {
final T cosLv = FastMath.cos(lv);
final T sinLv = FastMath.sin(lv);
final T num = auxiliaryElements.getH().multiply(cosLv).subtract(auxiliaryElements.getK().multiply(sinLv));
final T den = auxiliaryElements.getB().add(auxiliaryElements.getK().multiply(cosLv)).add(auxiliaryElements.getH().multiply(sinLv)).add(1.);
return FastMath.atan(num.divide(den)).multiply(2.).add(lv);
}
/** Converts eccentric longitude to mean longitude.
* @param le Eccentric longitude
* @return Mean longitude
*/
private T eccentricToMean (final T le) {
return le.subtract(auxiliaryElements.getK().multiply(FastMath.sin(le))).add(auxiliaryElements.getH().multiply(FastMath.cos(le)));
}
/** Compute δa/δv.
* @param vel satellite velocity
* @return δa/δv
*/
private FieldVector3D<T> getAoV(final FieldVector3D<T> vel) {
return new FieldVector3D<>(context.getTon2a(), vel);
}
/** Compute δh/δv.
* @param X satellite position component along f, equinoctial reference frame 1st vector
* @param Y satellite position component along g, equinoctial reference frame 2nd vector
* @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
* @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
* @return δh/δv
*/
private FieldVector3D<T> getHoV(final T X, final T Y, final T Xdot, final T Ydot) {
final T kf = (Xdot.multiply(Y).multiply(2.).subtract(X.multiply(Ydot))).multiply(context.getOoMU());
final T kg = X.multiply(Xdot).multiply(context.getOoMU());
final T kw = auxiliaryElements.getK().multiply(auxiliaryElements.getQ().multiply(Y).multiply(I).subtract(auxiliaryElements.getP().multiply(X))).multiply(context.getOOAB());
return new FieldVector3D<>(kf, auxiliaryElements.getVectorF(), kg.negate(), auxiliaryElements.getVectorG(), kw, auxiliaryElements.getVectorW());
}
/** Compute δk/δv.
* @param X satellite position component along f, equinoctial reference frame 1st vector
* @param Y satellite position component along g, equinoctial reference frame 2nd vector
* @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
* @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
* @return δk/δv
*/
private FieldVector3D<T> getKoV(final T X, final T Y, final T Xdot, final T Ydot) {
final T kf = Y.multiply(Ydot).multiply(context.getOoMU());
final T kg = (X.multiply(Ydot).multiply(2.).subtract(Xdot.multiply(Y))).multiply(context.getOoMU());
final T kw = auxiliaryElements.getH().multiply(auxiliaryElements.getQ().multiply(Y).multiply(I).subtract(auxiliaryElements.getP().multiply(X))).multiply(context.getOOAB());
return new FieldVector3D<>(kf.negate(), auxiliaryElements.getVectorF(), kg, auxiliaryElements.getVectorG(), kw.negate(), auxiliaryElements.getVectorW());
}
/** Compute δp/δv.
* @param Y satellite position component along g, equinoctial reference frame 2nd vector
* @return δp/δv
*/
private FieldVector3D<T> getPoV(final T Y) {
return new FieldVector3D<>(context.getCo2AB().multiply(Y), auxiliaryElements.getVectorW());
}
/** Compute δq/δv.
* @param X satellite position component along f, equinoctial reference frame 1st vector
* @return δq/δv
*/
private FieldVector3D<T> getQoV(final T X) {
return new FieldVector3D<>(context.getCo2AB().multiply(X).multiply(I), auxiliaryElements.getVectorW());
}
/** Compute δλ/δv.
* @param X satellite position component along f, equinoctial reference frame 1st vector
* @param Y satellite position component along g, equinoctial reference frame 2nd vector
* @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
* @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
* @return δλ/δv
*/
private FieldVector3D<T> getLoV(final T X, final T Y, final T Xdot, final T Ydot) {
final FieldVector3D<T> pos = new FieldVector3D<>(X, auxiliaryElements.getVectorF(), Y, auxiliaryElements.getVectorG());
final FieldVector3D<T> v2 = new FieldVector3D<>(auxiliaryElements.getK(), getHoV(X, Y, Xdot, Ydot), auxiliaryElements.getH().negate(), getKoV(X, Y, Xdot, Ydot));
return new FieldVector3D<>(context.getOOA().multiply(-2.), pos, context.getOoBpo(), v2, context.getOOA().multiply(auxiliaryElements.getQ().multiply(Y).multiply(I).subtract(auxiliaryElements.getP().multiply(X))), auxiliaryElements.getVectorW());
}
}
/** Internal class for numerical quadrature. */
private class IntegrableFunction implements UnivariateVectorFunction {
/** Current state. */
private final SpacecraftState state;
/** Signal that this class is used to compute the values required by the mean element variations
* or by the short periodic element variations. */
private final boolean meanMode;
/** The j index.
* <p>
* Used only for short periodic variation. Ignored for mean elements variation.
* </p> */
private final int j;
/** Container for attributes. */
private final AbstractGaussianContributionContext context;
/** Auxiliary Elements. */
private final AuxiliaryElements auxiliaryElements;
/** Drivers for solar radiation and atmospheric drag forces. */
private final double[] parameters;
/** Build a new instance.
* @param state current state information: date, kinematics, attitude
* @param meanMode if true return the value associated to the mean elements variation,
* if false return the values associated to the short periodic elements variation
* @param j the j index. used only for short periodic variation. Ignored for mean elements variation.
* @param parameters values of the force model parameters
*/
IntegrableFunction(final SpacecraftState state,
final boolean meanMode,
final int j,
final double[] parameters) {
this.meanMode = meanMode;
this.j = j;
this.parameters = parameters;
this.auxiliaryElements = new AuxiliaryElements(state.getOrbit(), I);
this.context = new AbstractGaussianContributionContext(auxiliaryElements, parameters);
// remove derivatives from state
final double[] stateVector = new double[6];
OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngle.TRUE, stateVector, null);
final Orbit fixedOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(stateVector, null, PositionAngle.TRUE,
state.getDate(),
context.getMu(),
state.getFrame());
this.state = new SpacecraftState(fixedOrbit, state.getAttitude(), state.getMass());
}
/** {@inheritDoc} */
@Override
public double[] value(final double x) {
//Compute the time difference from the true longitude difference
final double shiftedLm = trueToMean(x);
final double dLm = shiftedLm - auxiliaryElements.getLM();
final double dt = dLm / context.getMeanMotion();
final double cosL = FastMath.cos(x);
final double sinL = FastMath.sin(x);
final double roa = auxiliaryElements.getB() * auxiliaryElements.getB() / (1. + auxiliaryElements.getH() * sinL + auxiliaryElements.getK() * cosL);
final double roa2 = roa * roa;
final double r = auxiliaryElements.getSma() * roa;
final double X = r * cosL;
final double Y = r * sinL;
final double naob = context.getMeanMotion() * auxiliaryElements.getSma() / auxiliaryElements.getB();
final double Xdot = -naob * (auxiliaryElements.getH() + sinL);
final double Ydot = naob * (auxiliaryElements.getK() + cosL);
final Vector3D vel = new Vector3D(Xdot, auxiliaryElements.getVectorF(), Ydot, auxiliaryElements.getVectorG());
// Compute acceleration
Vector3D acc = Vector3D.ZERO;
// shift the orbit to dt
final Orbit shiftedOrbit = state.getOrbit().shiftedBy(dt);
// Recompose an orbit with time held fixed to be compliant with DSST theory
final Orbit recomposedOrbit =
new EquinoctialOrbit(shiftedOrbit.getA(),
shiftedOrbit.getEquinoctialEx(),
shiftedOrbit.getEquinoctialEy(),
shiftedOrbit.getHx(),
shiftedOrbit.getHy(),
shiftedOrbit.getLv(),
PositionAngle.TRUE,
shiftedOrbit.getFrame(),
state.getDate(),
context.getMu());
// Get the corresponding attitude
final Attitude recomposedAttitude =
attitudeProvider.getAttitude(recomposedOrbit,
recomposedOrbit.getDate(),
recomposedOrbit.getFrame());
// create shifted SpacecraftState with attitude at specified time
final SpacecraftState shiftedState =
new SpacecraftState(recomposedOrbit, recomposedAttitude, state.getMass());
acc = contribution.acceleration(shiftedState, parameters);
//Compute the derivatives of the elements by the speed
final double[] deriv = new double[6];
// da/dv
deriv[0] = getAoV(vel).dotProduct(acc);
// dex/dv
deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
// dey/dv
deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
// dhx/dv
deriv[3] = getQoV(X).dotProduct(acc);
// dhy/dv
deriv[4] = getPoV(Y).dotProduct(acc);
// dλ/dv
deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);
// Compute mean elements rates
double[] val = null;
if (meanMode) {
val = new double[6];
for (int i = 0; i < 6; i++) {
// da<sub>i</sub>/dt
val[i] = roa2 * deriv[i];
}
} else {
val = new double[12];
//Compute cos(j*L) and sin(j*L);
final double cosjL = j == 1 ? cosL : FastMath.cos(j * x);
final double sinjL = j == 1 ? sinL : FastMath.sin(j * x);
for (int i = 0; i < 6; i++) {
// da<sub>i</sub>/dv * cos(jL)
val[i] = cosjL * deriv[i];
// da<sub>i</sub>/dv * sin(jL)
val[i + 6] = sinjL * deriv[i];
}
}
return val;
}
/** Converts true longitude to eccentric longitude.
* @param lv True longitude
* @return Eccentric longitude
*/
private double trueToEccentric (final double lv) {
final double cosLv = FastMath.cos(lv);
final double sinLv = FastMath.sin(lv);
final double num = auxiliaryElements.getH() * cosLv - auxiliaryElements.getK() * sinLv;
final double den = auxiliaryElements.getB() + 1. + auxiliaryElements.getK() * cosLv + auxiliaryElements.getH() * sinLv;
return lv + 2. * FastMath.atan(num / den);
}
/** Converts eccentric longitude to mean longitude.
* @param le Eccentric longitude
* @return Mean longitude
*/
private double eccentricToMean (final double le) {
return le - auxiliaryElements.getK() * FastMath.sin(le) + auxiliaryElements.getH() * FastMath.cos(le);
}
/** Converts true longitude to mean longitude.
* @param lv True longitude
* @return Eccentric longitude
*/
private double trueToMean (final double lv) {
return eccentricToMean(trueToEccentric(lv));
}
/** Compute δa/δv.
* @param vel satellite velocity
* @return δa/δv
*/
private Vector3D getAoV(final Vector3D vel) {
return new Vector3D(context.getTon2a(), vel);
}
/** Compute δh/δv.
* @param X satellite position component along f, equinoctial reference frame 1st vector
* @param Y satellite position component along g, equinoctial reference frame 2nd vector
* @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
* @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
* @return δh/δv
*/
private Vector3D getHoV(final double X, final double Y, final double Xdot, final double Ydot) {
final double kf = (2. * Xdot * Y - X * Ydot) * context.getOoMU();
final double kg = X * Xdot * context.getOoMU();
final double kw = auxiliaryElements.getK() * (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOAB();
return new Vector3D(kf, auxiliaryElements.getVectorF(), -kg, auxiliaryElements.getVectorG(), kw, auxiliaryElements.getVectorW());
}
/** Compute δk/δv.
* @param X satellite position component along f, equinoctial reference frame 1st vector
* @param Y satellite position component along g, equinoctial reference frame 2nd vector
* @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
* @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
* @return δk/δv
*/
private Vector3D getKoV(final double X, final double Y, final double Xdot, final double Ydot) {
final double kf = Y * Ydot * context.getOoMU();
final double kg = (2. * X * Ydot - Xdot * Y) * context.getOoMU();
final double kw = auxiliaryElements.getH() * (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOAB();
return new Vector3D(-kf, auxiliaryElements.getVectorF(), kg, auxiliaryElements.getVectorG(), -kw, auxiliaryElements.getVectorW());
}
/** Compute δp/δv.
* @param Y satellite position component along g, equinoctial reference frame 2nd vector
* @return δp/δv
*/
private Vector3D getPoV(final double Y) {
return new Vector3D(context.getCo2AB() * Y, auxiliaryElements.getVectorW());
}
/** Compute δq/δv.
* @param X satellite position component along f, equinoctial reference frame 1st vector
* @return δq/δv
*/
private Vector3D getQoV(final double X) {
return new Vector3D(I * context.getCo2AB() * X, auxiliaryElements.getVectorW());
}
/** Compute δλ/δv.
* @param X satellite position component along f, equinoctial reference frame 1st vector
* @param Y satellite position component along g, equinoctial reference frame 2nd vector
* @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
* @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
* @return δλ/δv
*/
private Vector3D getLoV(final double X, final double Y, final double Xdot, final double Ydot) {
final Vector3D pos = new Vector3D(X, auxiliaryElements.getVectorF(), Y, auxiliaryElements.getVectorG());
final Vector3D v2 = new Vector3D(auxiliaryElements.getK(), getHoV(X, Y, Xdot, Ydot), -auxiliaryElements.getH(), getKoV(X, Y, Xdot, Ydot));
return new Vector3D(-2. * context.getOOA(), pos, context.getOoBpo(), v2, (I * auxiliaryElements.getQ() * Y - auxiliaryElements.getP() * X) * context.getOOA(), auxiliaryElements.getVectorW());
}
}
/** Class used to {@link #integrate(UnivariateVectorFunction, double, double) integrate}
* a {@link org.hipparchus.analysis.UnivariateVectorFunction function}
* of the orbital elements using the Gaussian quadrature rule to get the acceleration.
*/
private static class GaussQuadrature {
// Points and weights for the available quadrature orders
/** Points for quadrature of order 12. */
private static final double[] P_12 = {
-0.98156063424671910000,
-0.90411725637047490000,
-0.76990267419430470000,
-0.58731795428661740000,
-0.36783149899818024000,
-0.12523340851146890000,
0.12523340851146890000,
0.36783149899818024000,
0.58731795428661740000,
0.76990267419430470000,
0.90411725637047490000,
0.98156063424671910000
};
/** Weights for quadrature of order 12. */
private static final double[] W_12 = {
0.04717533638651220000,
0.10693932599531830000,
0.16007832854334633000,
0.20316742672306584000,
0.23349253653835478000,
0.24914704581340286000,
0.24914704581340286000,
0.23349253653835478000,
0.20316742672306584000,
0.16007832854334633000,
0.10693932599531830000,
0.04717533638651220000
};
/** Points for quadrature of order 16. */
private static final double[] P_16 = {
-0.98940093499164990000,
-0.94457502307323260000,
-0.86563120238783160000,
-0.75540440835500310000,
-0.61787624440264380000,
-0.45801677765722737000,
-0.28160355077925890000,
-0.09501250983763745000,
0.09501250983763745000,
0.28160355077925890000,
0.45801677765722737000,
0.61787624440264380000,
0.75540440835500310000,
0.86563120238783160000,
0.94457502307323260000,
0.98940093499164990000
};
/** Weights for quadrature of order 16. */
private static final double[] W_16 = {
0.02715245941175405800,
0.06225352393864777000,
0.09515851168249283000,
0.12462897125553388000,
0.14959598881657685000,
0.16915651939500256000,
0.18260341504492360000,
0.18945061045506847000,
0.18945061045506847000,
0.18260341504492360000,
0.16915651939500256000,
0.14959598881657685000,
0.12462897125553388000,
0.09515851168249283000,
0.06225352393864777000,
0.02715245941175405800
};
/** Points for quadrature of order 20. */
private static final double[] P_20 = {
-0.99312859918509490000,
-0.96397192727791390000,
-0.91223442825132600000,
-0.83911697182221890000,
-0.74633190646015080000,
-0.63605368072651510000,
-0.51086700195082700000,
-0.37370608871541955000,
-0.22778585114164507000,
-0.07652652113349734000,
0.07652652113349734000,
0.22778585114164507000,
0.37370608871541955000,
0.51086700195082700000,
0.63605368072651510000,
0.74633190646015080000,
0.83911697182221890000,
0.91223442825132600000,
0.96397192727791390000,
0.99312859918509490000
};
/** Weights for quadrature of order 20. */
private static final double[] W_20 = {
0.01761400713915226400,
0.04060142980038684000,
0.06267204833410904000,
0.08327674157670477000,
0.10193011981724048000,
0.11819453196151844000,
0.13168863844917678000,
0.14209610931838212000,
0.14917298647260380000,
0.15275338713072600000,
0.15275338713072600000,
0.14917298647260380000,
0.14209610931838212000,
0.13168863844917678000,
0.11819453196151844000,
0.10193011981724048000,
0.08327674157670477000,
0.06267204833410904000,
0.04060142980038684000,
0.01761400713915226400
};
/** Points for quadrature of order 24. */
private static final double[] P_24 = {
-0.99518721999702130000,
-0.97472855597130950000,
-0.93827455200273270000,
-0.88641552700440100000,
-0.82000198597390300000,
-0.74012419157855440000,
-0.64809365193697550000,
-0.54542147138883950000,
-0.43379350762604520000,
-0.31504267969616340000,
-0.19111886747361634000,
-0.06405689286260563000,
0.06405689286260563000,
0.19111886747361634000,
0.31504267969616340000,
0.43379350762604520000,
0.54542147138883950000,
0.64809365193697550000,
0.74012419157855440000,
0.82000198597390300000,
0.88641552700440100000,
0.93827455200273270000,
0.97472855597130950000,
0.99518721999702130000
};
/** Weights for quadrature of order 24. */
private static final double[] W_24 = {
0.01234122979998733500,
0.02853138862893380600,
0.04427743881741981000,
0.05929858491543691500,
0.07334648141108027000,
0.08619016153195320000,
0.09761865210411391000,
0.10744427011596558000,
0.11550566805372553000,
0.12167047292780335000,
0.12583745634682825000,
0.12793819534675221000,
0.12793819534675221000,
0.12583745634682825000,
0.12167047292780335000,
0.11550566805372553000,
0.10744427011596558000,
0.09761865210411391000,
0.08619016153195320000,
0.07334648141108027000,
0.05929858491543691500,
0.04427743881741981000,
0.02853138862893380600,
0.01234122979998733500
};
/** Points for quadrature of order 32. */
private static final double[] P_32 = {
-0.99726386184948160000,
-0.98561151154526840000,
-0.96476225558750640000,
-0.93490607593773970000,
-0.89632115576605220000,
-0.84936761373256990000,
-0.79448379596794250000,
-0.73218211874028970000,
-0.66304426693021520000,
-0.58771575724076230000,
-0.50689990893222950000,
-0.42135127613063540000,
-0.33186860228212767000,
-0.23928736225213710000,
-0.14447196158279646000,
-0.04830766568773831000,
0.04830766568773831000,
0.14447196158279646000,
0.23928736225213710000,
0.33186860228212767000,
0.42135127613063540000,
0.50689990893222950000,
0.58771575724076230000,
0.66304426693021520000,
0.73218211874028970000,
0.79448379596794250000,
0.84936761373256990000,
0.89632115576605220000,
0.93490607593773970000,
0.96476225558750640000,
0.98561151154526840000,
0.99726386184948160000
};
/** Weights for quadrature of order 32. */
private static final double[] W_32 = {
0.00701861000947013600,
0.01627439473090571200,
0.02539206530926214200,
0.03427386291302141000,
0.04283589802222658600,
0.05099805926237621600,
0.05868409347853559000,
0.06582222277636193000,
0.07234579410884862000,
0.07819389578707042000,
0.08331192422694673000,
0.08765209300440380000,
0.09117387869576390000,
0.09384439908080441000,
0.09563872007927487000,
0.09654008851472784000,
0.09654008851472784000,
0.09563872007927487000,
0.09384439908080441000,
0.09117387869576390000,
0.08765209300440380000,
0.08331192422694673000,
0.07819389578707042000,
0.07234579410884862000,
0.06582222277636193000,
0.05868409347853559000,
0.05099805926237621600,
0.04283589802222658600,
0.03427386291302141000,
0.02539206530926214200,
0.01627439473090571200,
0.00701861000947013600
};
/** Points for quadrature of order 40. */
private static final double[] P_40 = {
-0.99823770971055930000,
-0.99072623869945710000,
-0.97725994998377420000,
-0.95791681921379170000,
-0.93281280827867660000,
-0.90209880696887420000,
-0.86595950321225960000,
-0.82461223083331170000,
-0.77830565142651940000,
-0.72731825518992710000,
-0.67195668461417960000,
-0.61255388966798030000,
-0.54946712509512820000,
-0.48307580168617870000,
-0.41377920437160500000,
-0.34199409082575850000,
-0.26815218500725370000,
-0.19269758070137110000,
-0.11608407067525522000,
-0.03877241750605081600,
0.03877241750605081600,
0.11608407067525522000,
0.19269758070137110000,
0.26815218500725370000,
0.34199409082575850000,
0.41377920437160500000,
0.48307580168617870000,
0.54946712509512820000,
0.61255388966798030000,
0.67195668461417960000,
0.72731825518992710000,
0.77830565142651940000,
0.82461223083331170000,
0.86595950321225960000,
0.90209880696887420000,
0.93281280827867660000,
0.95791681921379170000,
0.97725994998377420000,
0.99072623869945710000,
0.99823770971055930000
};
/** Weights for quadrature of order 40. */
private static final double[] W_40 = {
0.00452127709853309800,
0.01049828453115270400,
0.01642105838190797300,
0.02224584919416689000,
0.02793700698002338000,
0.03346019528254786500,
0.03878216797447199000,
0.04387090818567333000,
0.04869580763507221000,
0.05322784698393679000,
0.05743976909939157000,
0.06130624249292891000,
0.06480401345660108000,
0.06791204581523394000,
0.07061164739128681000,
0.07288658239580408000,
0.07472316905796833000,
0.07611036190062619000,
0.07703981816424793000,
0.07750594797842482000,
0.07750594797842482000,
0.07703981816424793000,
0.07611036190062619000,
0.07472316905796833000,
0.07288658239580408000,
0.07061164739128681000,
0.06791204581523394000,
0.06480401345660108000,
0.06130624249292891000,
0.05743976909939157000,
0.05322784698393679000,
0.04869580763507221000,
0.04387090818567333000,
0.03878216797447199000,
0.03346019528254786500,
0.02793700698002338000,
0.02224584919416689000,
0.01642105838190797300,
0.01049828453115270400,
0.00452127709853309800
};
/** Points for quadrature of order 48. */
private static final double[] P_48 = {
-0.99877100725242610000,
-0.99353017226635080000,
-0.98412458372282700000,
-0.97059159254624720000,
-0.95298770316043080000,
-0.93138669070655440000,
-0.90587913671556960000,
-0.87657202027424800000,
-0.84358826162439350000,
-0.80706620402944250000,
-0.76715903251574020000,
-0.72403413092381470000,
-0.67787237963266400000,
-0.62886739677651370000,
-0.57722472608397270000,
-0.52316097472223300000,
-0.46690290475095840000,
-0.40868648199071680000,
-0.34875588629216070000,
-0.28736248735545555000,
-0.22476379039468908000,
-0.16122235606889174000,
-0.09700469920946270000,
-0.03238017096286937000,
0.03238017096286937000,
0.09700469920946270000,
0.16122235606889174000,
0.22476379039468908000,
0.28736248735545555000,
0.34875588629216070000,
0.40868648199071680000,
0.46690290475095840000,
0.52316097472223300000,
0.57722472608397270000,
0.62886739677651370000,
0.67787237963266400000,
0.72403413092381470000,
0.76715903251574020000,
0.80706620402944250000,
0.84358826162439350000,
0.87657202027424800000,
0.90587913671556960000,
0.93138669070655440000,
0.95298770316043080000,
0.97059159254624720000,
0.98412458372282700000,
0.99353017226635080000,
0.99877100725242610000
};
/** Weights for quadrature of order 48. */
private static final double[] W_48 = {
0.00315334605230596250,
0.00732755390127620800,
0.01147723457923446900,
0.01557931572294386600,
0.01961616045735556700,
0.02357076083932435600,
0.02742650970835688000,
0.03116722783279807000,
0.03477722256477045000,
0.03824135106583080600,
0.04154508294346483000,
0.04467456085669424000,
0.04761665849249054000,
0.05035903555385448000,
0.05289018948519365000,
0.05519950369998416500,
0.05727729210040315000,
0.05911483969839566000,
0.06070443916589384000,
0.06203942315989268000,
0.06311419228625403000,
0.06392423858464817000,
0.06446616443595010000,
0.06473769681268386000,
0.06473769681268386000,
0.06446616443595010000,
0.06392423858464817000,
0.06311419228625403000,
0.06203942315989268000,
0.06070443916589384000,
0.05911483969839566000,
0.05727729210040315000,
0.05519950369998416500,
0.05289018948519365000,
0.05035903555385448000,
0.04761665849249054000,
0.04467456085669424000,
0.04154508294346483000,
0.03824135106583080600,
0.03477722256477045000,
0.03116722783279807000,
0.02742650970835688000,
0.02357076083932435600,
0.01961616045735556700,
0.01557931572294386600,
0.01147723457923446900,
0.00732755390127620800,
0.00315334605230596250
};
/** Node points. */
private final double[] nodePoints;
/** Node weights. */
private final double[] nodeWeights;
/** Number of points. */
private final int numberOfPoints;
/** Creates a Gauss integrator of the given order.
*
* @param numberOfPoints Order of the integration rule.
*/
GaussQuadrature(final int numberOfPoints) {
this.numberOfPoints = numberOfPoints;
switch(numberOfPoints) {
case 12 :
this.nodePoints = P_12.clone();
this.nodeWeights = W_12.clone();
break;
case 16 :
this.nodePoints = P_16.clone();
this.nodeWeights = W_16.clone();
break;
case 20 :
this.nodePoints = P_20.clone();
this.nodeWeights = W_20.clone();
break;
case 24 :
this.nodePoints = P_24.clone();
this.nodeWeights = W_24.clone();
break;
case 32 :
this.nodePoints = P_32.clone();
this.nodeWeights = W_32.clone();
break;
case 40 :
this.nodePoints = P_40.clone();
this.nodeWeights = W_40.clone();
break;
case 48 :
default :
this.nodePoints = P_48.clone();
this.nodeWeights = W_48.clone();
break;
}
}
/** Integrates a given function on the given interval.
*
* @param f Function to integrate.
* @param lowerBound Lower bound of the integration interval.
* @param upperBound Upper bound of the integration interval.
* @return the integral of the weighted function.
*/
public double[] integrate(final UnivariateVectorFunction f,
final double lowerBound, final double upperBound) {
final double[] adaptedPoints = nodePoints.clone();
final double[] adaptedWeights = nodeWeights.clone();
transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
return basicIntegrate(f, adaptedPoints, adaptedWeights);
}
/** Integrates a given function on the given interval.
*
* @param <T> the type of the field elements
* @param f Function to integrate.
* @param lowerBound Lower bound of the integration interval.
* @param upperBound Upper bound of the integration interval.
* @param field field utilized by default
* @return the integral of the weighted function.
*/
public <T extends RealFieldElement<T>> T[] integrate(final RealFieldUnivariateVectorFunction<T> f,
final T lowerBound, final T upperBound,
final Field<T> field) {
final T zero = field.getZero();
final T[] adaptedPoints = MathArrays.buildArray(field, numberOfPoints);
final T[] adaptedWeights = MathArrays.buildArray(field, numberOfPoints);;
for (int i = 0; i < numberOfPoints; i++) {
adaptedPoints[i] = zero.add(nodePoints[i]);
adaptedWeights[i] = zero.add(nodeWeights[i]);
}
transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
return basicIntegrate(f, adaptedPoints, adaptedWeights, field);
}
/** Performs a change of variable so that the integration
* can be performed on an arbitrary interval {@code [a, b]}.
* <p>
* It is assumed that the natural interval is {@code [-1, 1]}.
* </p>
*
* @param points Points to adapt to the new interval.
* @param weights Weights to adapt to the new interval.
* @param a Lower bound of the integration interval.
* @param b Lower bound of the integration interval.
*/
private void transform(final double[] points, final double[] weights,
final double a, final double b) {
// Scaling
final double scale = (b - a) / 2;
final double shift = a + scale;
for (int i = 0; i < points.length; i++) {
points[i] = points[i] * scale + shift;
weights[i] *= scale;
}
}
/** Performs a change of variable so that the integration
* can be performed on an arbitrary interval {@code [a, b]}.
* <p>
* It is assumed that the natural interval is {@code [-1, 1]}.
* </p>
* @param <T> the type of the field elements
* @param points Points to adapt to the new interval.
* @param weights Weights to adapt to the new interval.
* @param a Lower bound of the integration interval.
* @param b Lower bound of the integration interval
*/
private <T extends RealFieldElement<T>> void transform(final T[] points, final T[] weights,
final T a, final T b) {
// Scaling
final T scale = (b.subtract(a)).divide(2.);
final T shift = a.add(scale);
for (int i = 0; i < points.length; i++) {
points[i] = scale.multiply(points[i]).add(shift);
weights[i] = scale.multiply(weights[i]);
}
}
/** Returns an estimate of the integral of {@code f(x) * w(x)},
* where {@code w} is a weight function that depends on the actual
* flavor of the Gauss integration scheme.
*
* @param f Function to integrate.
* @param points Nodes.
* @param weights Nodes weights.
* @return the integral of the weighted function.
*/
private double[] basicIntegrate(final UnivariateVectorFunction f,
final double[] points,
final double[] weights) {
double x = points[0];
double w = weights[0];
double[] v = f.value(x);
final double[] y = new double[v.length];
for (int j = 0; j < v.length; j++) {
y[j] = w * v[j];
}
final double[] t = y.clone();
final double[] c = new double[v.length];
final double[] s = t.clone();
for (int i = 1; i < points.length; i++) {
x = points[i];
w = weights[i];
v = f.value(x);
for (int j = 0; j < v.length; j++) {
y[j] = w * v[j] - c[j];
t[j] = s[j] + y[j];
c[j] = (t[j] - s[j]) - y[j];
s[j] = t[j];
}
}
return s;
}
/** Returns an estimate of the integral of {@code f(x) * w(x)},
* where {@code w} is a weight function that depends on the actual
* flavor of the Gauss integration scheme.
*
* @param <T> the type of the field elements.
* @param f Function to integrate.
* @param points Nodes.
* @param weights Nodes weight
* @param field field utilized by default
* @return the integral of the weighted function.
*/
private <T extends RealFieldElement<T>> T[] basicIntegrate(final RealFieldUnivariateVectorFunction<T> f,
final T[] points,
final T[] weights,
final Field<T> field) {
T x = points[0];
T w = weights[0];
T[] v = f.value(x);
final T[] y = MathArrays.buildArray(field, v.length);
for (int j = 0; j < v.length; j++) {
y[j] = v[j].multiply(w);
}
final T[] t = y.clone();
final T[] c = MathArrays.buildArray(field, v.length);;
final T[] s = t.clone();
for (int i = 1; i < points.length; i++) {
x = points[i];
w = weights[i];
v = f.value(x);
for (int j = 0; j < v.length; j++) {
y[j] = v[j].multiply(w).subtract(c[j]);
t[j] = y[j].add(s[j]);
c[j] = (t[j].subtract(s[j])).subtract(y[j]);
s[j] = t[j];
}
}
return s;
}
}
/** Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
* <p>
* Those coefficients are given in Danielson paper by expression 4.4-(6)
* </p>
* @author Petre Bazavan
* @author Lucian Barbulescu
*/
private class FourierCjSjCoefficients {
/** Maximum possible value for j. */
private final int jMax;
/** The C<sub>i</sub><sup>j</sup> coefficients.
* <p>
* the index i corresponds to the following elements: <br/>
* - 0 for a <br>
* - 1 for k <br>
* - 2 for h <br>
* - 3 for q <br>
* - 4 for p <br>
* - 5 for λ <br>
* </p>
*/
private final double[][] cCoef;
/** The C<sub>i</sub><sup>j</sup> coefficients.
* <p>
* the index i corresponds to the following elements: <br/>
* - 0 for a <br>
* - 1 for k <br>
* - 2 for h <br>
* - 3 for q <br>
* - 4 for p <br>
* - 5 for λ <br>
* </p>
*/
private final double[][] sCoef;
/** Standard constructor.
* @param state the current state
* @param jMax maximum value for j
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param parameters values of the force model parameters
*/
FourierCjSjCoefficients(final SpacecraftState state, final int jMax,
final AuxiliaryElements auxiliaryElements, final double[] parameters) {
//Initialise the fields
this.jMax = jMax;
//Allocate the arrays
final int rows = jMax + 1;
cCoef = new double[rows][6];
sCoef = new double[rows][6];
//Compute the coefficients
computeCoefficients(state, auxiliaryElements, parameters);
}
/**
* Compute the Fourrier coefficients.
* <p>
* Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients need to be computed
* as D<sub>i</sub><sup>m</sup> is always 0.
* </p>
* @param state the current state
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param parameters values of the force model parameters
*/
private void computeCoefficients(final SpacecraftState state,
final AuxiliaryElements auxiliaryElements,
final double[] parameters) {
// Computes the limits for the integral
final double[] ll = getLLimits(state, auxiliaryElements);
// Computes integrated mean element rates if Llow < Lhigh
if (ll[0] < ll[1]) {
//Compute 1 / PI
final double ooPI = 1 / FastMath.PI;
// loop through all values of j
for (int j = 0; j <= jMax; j++) {
final double[] curentCoefficients =
integrator.integrate(new IntegrableFunction(state, false, j, parameters), ll[0], ll[1]);
//divide by PI and set the values for the coefficients
for (int i = 0; i < 6; i++) {
cCoef[j][i] = ooPI * curentCoefficients[i];
sCoef[j][i] = ooPI * curentCoefficients[i + 6];
}
}
}
}
/** Get the coefficient C<sub>i</sub><sup>j</sup>.
* @param i i index - corresponds to the required variation
* @param j j index
* @return the coefficient C<sub>i</sub><sup>j</sup>
*/
public double getCij(final int i, final int j) {
return cCoef[j][i];
}
/** Get the coefficient S<sub>i</sub><sup>j</sup>.
* @param i i index - corresponds to the required variation
* @param j j index
* @return the coefficient S<sub>i</sub><sup>j</sup>
*/
public double getSij(final int i, final int j) {
return sCoef[j][i];
}
}
/** Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients with field elements.
* <p>
* Those coefficients are given in Danielson paper by expression 4.4-(6)
* </p>
* @author Petre Bazavan
* @author Lucian Barbulescu
*/
private class FieldFourierCjSjCoefficients <T extends RealFieldElement<T>> {
/** Maximum possible value for j. */
private final int jMax;
/** The C<sub>i</sub><sup>j</sup> coefficients.
* <p>
* the index i corresponds to the following elements: <br/>
* - 0 for a <br>
* - 1 for k <br>
* - 2 for h <br>
* - 3 for q <br>
* - 4 for p <br>
* - 5 for λ <br>
* </p>
*/
private final T[][] cCoef;
/** The C<sub>i</sub><sup>j</sup> coefficients.
* <p>
* the index i corresponds to the following elements: <br/>
* - 0 for a <br>
* - 1 for k <br>
* - 2 for h <br>
* - 3 for q <br>
* - 4 for p <br>
* - 5 for λ <br>
* </p>
*/
private final T[][] sCoef;
/** Standard constructor.
* @param state the current state
* @param jMax maximum value for j
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param parameters values of the force model parameters
* @param field field used by default
*/
FieldFourierCjSjCoefficients(final FieldSpacecraftState<T> state, final int jMax,
final FieldAuxiliaryElements<T> auxiliaryElements,
final T[] parameters,
final Field<T> field) {
//Initialise the fields
this.jMax = jMax;
//Allocate the arrays
final int rows = jMax + 1;
cCoef = MathArrays.buildArray(field, rows, 6);
sCoef = MathArrays.buildArray(field, rows, 6);
//Compute the coefficients
computeCoefficients(state, auxiliaryElements, parameters, field);
}
/**
* Compute the Fourrier coefficients.
* <p>
* Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients need to be computed
* as D<sub>i</sub><sup>m</sup> is always 0.
* </p>
* @param state the current state
* @param auxiliaryElements auxiliary elements related to the current orbit
* @param parameters values of the force model parameters
* @param field field used by default
*/
private void computeCoefficients(final FieldSpacecraftState<T> state,
final FieldAuxiliaryElements<T> auxiliaryElements,
final T[] parameters,
final Field<T> field) {
// Zero
final T zero = field.getZero();
// Computes the limits for the integral
final T[] ll = getLLimits(state, auxiliaryElements);
// Computes integrated mean element rates if Llow < Lhigh
if (ll[0].getReal() < ll[1].getReal()) {
//Compute 1 / PI
final T ooPI = zero.add(FastMath.PI).reciprocal();
// loop through all values of j
for (int j = 0; j <= jMax; j++) {
final T[] curentCoefficients =
integrator.integrate(new FieldIntegrableFunction<>(state, false, j, parameters, field), ll[0], ll[1], field);
//divide by PI and set the values for the coefficients
for (int i = 0; i < 6; i++) {
cCoef[j][i] = curentCoefficients[i].multiply(ooPI);
sCoef[j][i] = curentCoefficients[i + 6].multiply(ooPI);
}
}
}
}
/** Get the coefficient C<sub>i</sub><sup>j</sup>.
* @param i i index - corresponds to the required variation
* @param j j index
* @return the coefficient C<sub>i</sub><sup>j</sup>
*/
public T getCij(final int i, final int j) {
return cCoef[j][i];
}
/** Get the coefficient S<sub>i</sub><sup>j</sup>.
* @param i i index - corresponds to the required variation
* @param j j index
* @return the coefficient S<sub>i</sub><sup>j</sup>
*/
public T getSij(final int i, final int j) {
return sCoef[j][i];
}
}
/** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
*
* <p>
* The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
* then the values of the coefficients D<sub>i</sub><sup>m</sup> for m > 2 are also 0.
* </p>
* @author Petre Bazavan
* @author Lucian Barbulescu
*
*/
private static class GaussianShortPeriodicCoefficients implements ShortPeriodTerms {
/** Maximum value for j index. */
private final int jMax;
/** Number of points used in the interpolation process. */
private final int interpolationPoints;
/** Prefix for coefficients keys. */
private final String coefficientsKeyPrefix;
/** All coefficients slots. */
private final transient TimeSpanMap<Slot> slots;
/** Constructor.
* @param coefficientsKeyPrefix prefix for coefficients keys
* @param jMax maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
* @param slots all coefficients slots
*/
GaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix,
final int jMax, final int interpolationPoints,
final TimeSpanMap<Slot> slots) {
//Initialize fields
this.jMax = jMax;
this.interpolationPoints = interpolationPoints;
this.coefficientsKeyPrefix = coefficientsKeyPrefix;
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(jMax, interpolationPoints);
final AbsoluteDate first = meanStates[0].getDate();
final AbsoluteDate last = meanStates[meanStates.length - 1].getDate();
if (first.compareTo(last) <= 0) {
slots.addValidAfter(slot, first);
} else {
slots.addValidBefore(slot, first);
}
return slot;
}
/** Compute the short periodic coefficients.
*
* @param state current state information: date, kinematics, attitude
* @param slot coefficients slot
* @param fourierCjSj Fourier coefficients
* @param uijvij U and V coefficients
* @param n Keplerian mean motion
* @param a semi major axis
*/
private void computeCoefficients(final SpacecraftState state, final Slot slot,
final FourierCjSjCoefficients fourierCjSj,
final UijVijCoefficients uijvij,
final double n, final double a) {
// get the current date
final AbsoluteDate date = state.getDate();
// compute the k₂⁰ coefficient
final double k20 = computeK20(jMax, uijvij.currentRhoSigmaj);
// 1. / n
final double oon = 1. / n;
// 3. / (2 * a * n)
final double to2an = 1.5 * oon / a;
// 3. / (4 * a * n)
final double to4an = to2an / 2;
// Compute the coefficients for each element
final int size = jMax + 1;
final double[] di1 = new double[6];
final double[] di2 = new double[6];
final double[][] currentCij = new double[size][6];
final double[][] currentSij = new double[size][6];
for (int i = 0; i < 6; i++) {
// compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
di1[i] = -oon * fourierCjSj.getCij(i, 0);
if (i == 5) {
di1[i] += to2an * uijvij.getU1(0, 0);
}
di2[i] = 0.;
if (i == 5) {
di2[i] += -to4an * fourierCjSj.getCij(0, 0);
}
//the C<sub>i</sub>⁰ is computed based on all others
currentCij[0][i] = -di2[i] * k20;
for (int j = 1; j <= jMax; j++) {
// compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
currentCij[j][i] = oon * uijvij.getU1(j, i);
if (i == 5) {
currentCij[j][i] += -to2an * uijvij.getU2(j);
}
currentSij[j][i] = oon * uijvij.getV1(j, i);
if (i == 5) {
currentSij[j][i] += -to2an * uijvij.getV2(j);
}
// add the computed coefficients to C<sub>i</sub>⁰
currentCij[0][i] += -(currentCij[j][i] * uijvij.currentRhoSigmaj[0][j] + currentSij[j][i] * uijvij.currentRhoSigmaj[1][j]);
}
}
// add the values to the interpolators
slot.cij[0].addGridPoint(date, currentCij[0]);
slot.dij[1].addGridPoint(date, di1);
slot.dij[2].addGridPoint(date, di2);
for (int j = 1; j <= jMax; j++) {
slot.cij[j].addGridPoint(date, currentCij[j]);
slot.sij[j].addGridPoint(date, currentSij[j]);
}
}
/** Compute the coefficient k₂⁰ by using the equation
* 2.5.3-(9a) from Danielson.
* <p>
* After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
* k₂⁰ = Σ<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
* </p>
* @param kMax max value fot k index
* @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
* @return the coefficient k₂⁰
*/
private double computeK20(final int kMax, final double[][] currentRhoSigmaj) {
double k20 = 0.;
for (int kIndex = 1; kIndex <= kMax; kIndex++) {
// After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:
//k₂⁰ = Σ<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
double currentTerm = currentRhoSigmaj[1][kIndex] * currentRhoSigmaj[1][kIndex] +
currentRhoSigmaj[0][kIndex] * currentRhoSigmaj[0][kIndex];
//multiply by 2 / k²
currentTerm *= 2. / (kIndex * kIndex);
// add the term to the result
k20 += currentTerm;
}
return k20;
}
/** {@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 (l - λ)
final double center = L - meanOrbit.getLM();
// Compute (l - λ)²
final double center2 = center * center;
// Initialize short periodic variations
final double[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
final double[] d1 = slot.dij[1].value(meanOrbit.getDate());
final double[] d2 = slot.dij[2].value(meanOrbit.getDate());
for (int i = 0; i < 6; i++) {
shortPeriodicVariation[i] += center * d1[i] + center2 * d2[i];
}
for (int j = 1; j <= JMAX; j++) {
final double[] c = slot.cij[j].value(meanOrbit.getDate());
final double[] s = slot.sij[j].value(meanOrbit.getDate());
final double cos = FastMath.cos(j * L);
final double sin = FastMath.sin(j * L);
for (int i = 0; i < 6; i++) {
// add corresponding term to the short periodic variation
shortPeriodicVariation[i] += c[i] * cos;
shortPeriodicVariation[i] += s[i] * sin;
}
}
return shortPeriodicVariation;
}
/** {@inheritDoc} */
public String getCoefficientsKeyPrefix() {
return coefficientsKeyPrefix;
}
/** {@inheritDoc}
* <p>
* For Gaussian forces, there are JMAX cj coefficients,
* JMAX sj coefficients and 3 dj coefficients. As JMAX = 12,
* this sums up to 27 coefficients. The j index is the integer
* multiplier for the true longitude argument in the cj and sj
* coefficients and to 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<String, double[]>(2 * JMAX + 3);
storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
for (int j = 1; j <= JMAX; 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);
}
}
}
/** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
*
* <p>
* The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
* then the values of the coefficients D<sub>i</sub><sup>m</sup> for m > 2 are also 0.
* </p>
* @author Petre Bazavan
* @author Lucian Barbulescu
*
*/
private static class FieldGaussianShortPeriodicCoefficients <T extends RealFieldElement<T>> implements FieldShortPeriodTerms<T> {
/** Maximum value for j index. */
private final int jMax;
/** Number of points used in the interpolation process. */
private final int interpolationPoints;
/** Prefix for coefficients keys. */
private final String coefficientsKeyPrefix;
/** All coefficients slots. */
private final transient FieldTimeSpanMap<FieldSlot<T>, T> slots;
/** Constructor.
* @param coefficientsKeyPrefix prefix for coefficients keys
* @param jMax maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
* @param slots all coefficients slots
*/
FieldGaussianShortPeriodicCoefficients(final String coefficientsKeyPrefix,
final int jMax, final int interpolationPoints,
final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
//Initialize fields
this.jMax = jMax;
this.interpolationPoints = interpolationPoints;
this.coefficientsKeyPrefix = coefficientsKeyPrefix;
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<>(jMax, 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;
}
/** Compute the short periodic coefficients.
*
* @param state current state information: date, kinematics, attitude
* @param slot coefficients slot
* @param fourierCjSj Fourier coefficients
* @param uijvij U and V coefficients
* @param n Keplerian mean motion
* @param a semi major axis
* @param field field used by default
*/
private void computeCoefficients(final FieldSpacecraftState<T> state, final FieldSlot<T> slot,
final FieldFourierCjSjCoefficients<T> fourierCjSj,
final FieldUijVijCoefficients<T> uijvij,
final T n, final T a,
final Field<T> field) {
// Zero
final T zero = field.getZero();
// get the current date
final FieldAbsoluteDate<T> date = state.getDate();
// compute the k₂⁰ coefficient
final T k20 = computeK20(jMax, uijvij.currentRhoSigmaj, field);
// 1. / n
final T oon = n.reciprocal();
// 3. / (2 * a * n)
final T to2an = oon.multiply(1.5).divide(a);
// 3. / (4 * a * n)
final T to4an = to2an.divide(2.);
// Compute the coefficients for each element
final int size = jMax + 1;
final T[] di1 = MathArrays.buildArray(field, 6);
final T[] di2 = MathArrays.buildArray(field, 6);
final T[][] currentCij = MathArrays.buildArray(field, size, 6);
final T[][] currentSij = MathArrays.buildArray(field, size, 6);
for (int i = 0; i < 6; i++) {
// compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
di1[i] = oon.negate().multiply(fourierCjSj.getCij(i, 0));
if (i == 5) {
di1[i] = di1[i].add(to2an.multiply(uijvij.getU1(0, 0)));
}
di2[i] = zero;
if (i == 5) {
di2[i] = di2[i].add(to4an.negate().multiply(fourierCjSj.getCij(0, 0)));
}
//the C<sub>i</sub>⁰ is computed based on all others
currentCij[0][i] = di2[i].negate().multiply(k20);
for (int j = 1; j <= jMax; j++) {
// compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
currentCij[j][i] = oon.multiply(uijvij.getU1(j, i));
if (i == 5) {
currentCij[j][i] = currentCij[j][i].add(to2an.negate().multiply(uijvij.getU2(j)));
}
currentSij[j][i] = oon.multiply(uijvij.getV1(j, i));
if (i == 5) {
currentSij[j][i] = currentSij[j][i].add(to2an.negate().multiply(uijvij.getV2(j)));
}
// add the computed coefficients to C<sub>i</sub>⁰
currentCij[0][i] = currentCij[0][i].add(currentCij[j][i].multiply(uijvij.currentRhoSigmaj[0][j]).add(currentSij[j][i].multiply(uijvij.currentRhoSigmaj[1][j])).negate());
}
}
// add the values to the interpolators
slot.cij[0].addGridPoint(date, currentCij[0]);
slot.dij[1].addGridPoint(date, di1);
slot.dij[2].addGridPoint(date, di2);
for (int j = 1; j <= jMax; j++) {
slot.cij[j].addGridPoint(date, currentCij[j]);
slot.sij[j].addGridPoint(date, currentSij[j]);
}
}
/** Compute the coefficient k₂⁰ by using the equation
* 2.5.3-(9a) from Danielson.
* <p>
* After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
* k₂⁰ = Σ<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
* </p>
* @param kMax max value fot k index
* @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
* @param field field used by default
* @return the coefficient k₂⁰
*/
private T computeK20(final int kMax, final T[][] currentRhoSigmaj, final Field<T> field) {
final T zero = field.getZero();
T k20 = zero;
for (int kIndex = 1; kIndex <= kMax; kIndex++) {
// After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:
//k₂⁰ = Σ<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
T currentTerm = currentRhoSigmaj[1][kIndex].multiply(currentRhoSigmaj[1][kIndex]).
add(currentRhoSigmaj[0][kIndex].multiply(currentRhoSigmaj[0][kIndex]));
//multiply by 2 / k²
currentTerm = currentTerm.multiply(2. / (kIndex * kIndex));
// add the term to the result
k20 = k20.add(currentTerm);
}
return k20;
}
/** {@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 (l - λ)
final T center = L.subtract(meanOrbit.getLM());
// Compute (l - λ)²
final T center2 = center.multiply(center);
// Initialize short periodic variations
final T[] shortPeriodicVariation = slot.cij[0].value(meanOrbit.getDate());
final T[] d1 = slot.dij[1].value(meanOrbit.getDate());
final T[] d2 = slot.dij[2].value(meanOrbit.getDate());
for (int i = 0; i < 6; i++) {
shortPeriodicVariation[i] = shortPeriodicVariation[i].add(center.multiply(d1[i]).add(center2.multiply(d2[i])));
}
for (int j = 1; j <= JMAX; j++) {
final T[] c = slot.cij[j].value(meanOrbit.getDate());
final T[] s = slot.sij[j].value(meanOrbit.getDate());
final T cos = FastMath.cos(L.multiply(j));
final T sin = FastMath.sin(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(cos));
shortPeriodicVariation[i] = shortPeriodicVariation[i].add(s[i].multiply(sin));
}
}
return shortPeriodicVariation;
}
/** {@inheritDoc} */
public String getCoefficientsKeyPrefix() {
return coefficientsKeyPrefix;
}
/** {@inheritDoc}
* <p>
* For Gaussian forces, there are JMAX cj coefficients,
* JMAX sj coefficients and 3 dj coefficients. As JMAX = 12,
* this sums up to 27 coefficients. The j index is the integer
* multiplier for the true longitude argument in the cj and sj
* coefficients and to 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<String, T[]>(2 * JMAX + 3);
storeIfSelected(coefficients, selected, slot.cij[0].value(date), "d", 0);
storeIfSelected(coefficients, selected, slot.dij[1].value(date), "d", 1);
storeIfSelected(coefficients, selected, slot.dij[2].value(date), "d", 2);
for (int j = 1; j <= JMAX; 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);
}
}
}
/** The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients described by
* equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
* <p>
* The index i takes only the values 1 and 2<br>
* For U only the index 0 for j is used.
* </p>
*
* @author Petre Bazavan
* @author Lucian Barbulescu
*/
private static class UijVijCoefficients {
/** The U₁<sup>j</sup> coefficients.
* <p>
* The first index identifies the Fourier coefficients used<br>
* Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
* The only exception is when j = 0 when only the coefficient for fourier index = 1 (i == 0) is needed.<br>
* Also, for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
* to compute the coefficients U₂<sup>j</sup>
* </p>
*/
private final double[][] u1ij;
/** The V₁<sup>j</sup> coefficients.
* <p>
* The first index identifies the Fourier coefficients used<br>
* Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
* for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
* to compute the coefficients V₂<sup>j</sup>
* </p>
*/
private final double[][] v1ij;
/** The U₂<sup>j</sup> coefficients.
* <p>
* Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
* </p>
*/
private final double[] u2ij;
/** The V₂<sup>j</sup> coefficients.
* <p>
* Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
* </p>
*/
private final double[] v2ij;
/** The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients. */
private final double[][] currentRhoSigmaj;
/** The C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> Fourier coefficients. */
private final FourierCjSjCoefficients fourierCjSj;
/** The maximum value for j index. */
private final int jMax;
/** Constructor.
* @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
* @param fourierCjSj the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
* @param jMax maximum value for j index
*/
UijVijCoefficients(final double[][] currentRhoSigmaj, final FourierCjSjCoefficients fourierCjSj, final int jMax) {
this.currentRhoSigmaj = currentRhoSigmaj;
this.fourierCjSj = fourierCjSj;
this.jMax = jMax;
// initialize the internal arrays.
this.u1ij = new double[6][2 * jMax + 1];
this.v1ij = new double[6][2 * jMax + 1];
this.u2ij = new double[jMax + 1];
this.v2ij = new double[jMax + 1];
//compute the coefficients
computeU1V1Coefficients();
computeU2V2Coefficients();
}
/** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients. */
private void computeU1V1Coefficients() {
// generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
// for j >= 1
// also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
u1ij[0][0] = 0;
for (int j = 1; j <= jMax; j++) {
// compute 1 / j
final double ooj = 1. / j;
for (int i = 0; i < 6; i++) {
//j is aready between 1 and J
u1ij[i][j] = fourierCjSj.getSij(i, j);
v1ij[i][j] = fourierCjSj.getCij(i, j);
// 1 - δ<sub>1j</sub> is 1 for all j > 1
if (j > 1) {
// k starts with 1 because j-J is less than or equal to 0
for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
// C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
// S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
u1ij[i][j] += fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[1][kIndex] +
fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[0][kIndex];
// C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
// S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
v1ij[i][j] += fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[0][kIndex] -
fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[1][kIndex];
}
}
// since j must be between 1 and J-1 and is already between 1 and J
// the following sum is skiped only for j = jMax
if (j != jMax) {
for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
// -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
// S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
u1ij[i][j] += -fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[1][kIndex] +
fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[0][kIndex];
// C<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub> +
// S<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub>
v1ij[i][j] += fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[0][kIndex] +
fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[1][kIndex];
}
}
for (int kIndex = 1; kIndex <= jMax; kIndex++) {
// C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
// S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
u1ij[i][j] += -fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[1][j + kIndex] -
fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[0][j + kIndex];
// C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
// S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
v1ij[i][j] += fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[0][j + kIndex] +
fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[1][j + kIndex];
}
// divide by 1 / j
u1ij[i][j] *= -ooj;
v1ij[i][j] *= ooj;
// if index = 1 (i == 0) add the computed terms to U₁⁰
if (i == 0) {
//- (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
u1ij[0][0] += -u1ij[0][j] * currentRhoSigmaj[0][j] - v1ij[0][j] * currentRhoSigmaj[1][j];
}
}
}
// Terms with j > jMax are required only when computing the coefficients
// U₂<sup>j</sup> and V₂<sup>j</sup>
// and those coefficients are only required for Fourier index = 1 (i == 0).
for (int j = jMax + 1; j <= 2 * jMax; j++) {
// compute 1 / j
final double ooj = 1. / j;
//the value of i is 0
u1ij[0][j] = 0.;
v1ij[0][j] = 0.;
//k starts from j-J as it is always greater than or equal to 1
for (int kIndex = j - jMax; kIndex <= j - 1; kIndex++) {
// C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
// S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
u1ij[0][j] += fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[1][kIndex] +
fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[0][kIndex];
// C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
// S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
v1ij[0][j] += fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[0][kIndex] -
fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[1][kIndex];
}
for (int kIndex = 1; kIndex <= jMax; kIndex++) {
// C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
// S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
u1ij[0][j] += -fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[1][j + kIndex] -
fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[0][j + kIndex];
// C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
// S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
v1ij[0][j] += fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[0][j + kIndex] +
fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[1][j + kIndex];
}
// divide by 1 / j
u1ij[0][j] *= -ooj;
v1ij[0][j] *= ooj;
}
}
/** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
* <p>
* Only the coefficients for Fourier index = 1 (i == 0) are required.
* </p>
*/
private void computeU2V2Coefficients() {
for (int j = 1; j <= jMax; j++) {
// compute 1 / j
final double ooj = 1. / j;
// only the values for i == 0 are computed
u2ij[j] = v1ij[0][j];
v2ij[j] = u1ij[0][j];
// 1 - δ<sub>1j</sub> is 1 for all j > 1
if (j > 1) {
for (int l = 1; l <= j - 1; l++) {
// U₁<sup>j-l</sup> * σ<sub>l</sub> +
// V₁<sup>j-l</sup> * ρ<sub>l</sub>
u2ij[j] += u1ij[0][j - l] * currentRhoSigmaj[1][l] +
v1ij[0][j - l] * currentRhoSigmaj[0][l];
// U₁<sup>j-l</sup> * ρ<sub>l</sub> -
// V₁<sup>j-l</sup> * σ<sub>l</sub>
v2ij[j] += u1ij[0][j - l] * currentRhoSigmaj[0][l] -
v1ij[0][j - l] * currentRhoSigmaj[1][l];
}
}
for (int l = 1; l <= jMax; l++) {
// -U₁<sup>j+l</sup> * σ<sub>l</sub> +
// U₁<sup>l</sup> * σ<sub>j+l</sub> +
// V₁<sup>j+l</sup> * ρ<sub>l</sub> -
// V₁<sup>l</sup> * ρ<sub>j+l</sub>
u2ij[j] += -u1ij[0][j + l] * currentRhoSigmaj[1][l] +
u1ij[0][l] * currentRhoSigmaj[1][j + l] +
v1ij[0][j + l] * currentRhoSigmaj[0][l] -
v1ij[0][l] * currentRhoSigmaj[0][j + l];
// U₁<sup>j+l</sup> * ρ<sub>l</sub> +
// U₁<sup>l</sup> * ρ<sub>j+l</sub> +
// V₁<sup>j+l</sup> * σ<sub>l</sub> +
// V₁<sup>l</sup> * σ<sub>j+l</sub>
u2ij[j] += u1ij[0][j + l] * currentRhoSigmaj[0][l] +
u1ij[0][l] * currentRhoSigmaj[0][j + l] +
v1ij[0][j + l] * currentRhoSigmaj[1][l] +
v1ij[0][l] * currentRhoSigmaj[1][j + l];
}
// divide by 1 / j
u2ij[j] *= -ooj;
v2ij[j] *= ooj;
}
}
/** Get the coefficient U₁<sup>j</sup> for Fourier index i.
*
* @param j j index
* @param i Fourier index (starts at 0)
* @return the coefficient U₁<sup>j</sup> for the given Fourier index i
*/
public double getU1(final int j, final int i) {
return u1ij[i][j];
}
/** Get the coefficient V₁<sup>j</sup> for Fourier index i.
*
* @param j j index
* @param i Fourier index (starts at 0)
* @return the coefficient V₁<sup>j</sup> for the given Fourier index i
*/
public double getV1(final int j, final int i) {
return v1ij[i][j];
}
/** Get the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0).
*
* @param j j index
* @return the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0)
*/
public double getU2(final int j) {
return u2ij[j];
}
/** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
*
* @param j j index
* @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
*/
public double getV2(final int j) {
return v2ij[j];
}
}
/** The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients described by
* equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
* <p>
* The index i takes only the values 1 and 2<br>
* For U only the index 0 for j is used.
* </p>
*
* @author Petre Bazavan
* @author Lucian Barbulescu
*/
private static class FieldUijVijCoefficients <T extends RealFieldElement<T>> {
/** The U₁<sup>j</sup> coefficients.
* <p>
* The first index identifies the Fourier coefficients used<br>
* Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
* The only exception is when j = 0 when only the coefficient for fourier index = 1 (i == 0) is needed.<br>
* Also, for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
* to compute the coefficients U₂<sup>j</sup>
* </p>
*/
private final T[][] u1ij;
/** The V₁<sup>j</sup> coefficients.
* <p>
* The first index identifies the Fourier coefficients used<br>
* Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
* for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
* to compute the coefficients V₂<sup>j</sup>
* </p>
*/
private final T[][] v1ij;
/** The U₂<sup>j</sup> coefficients.
* <p>
* Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
* </p>
*/
private final T[] u2ij;
/** The V₂<sup>j</sup> coefficients.
* <p>
* Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
* </p>
*/
private final T[] v2ij;
/** The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients. */
private final T[][] currentRhoSigmaj;
/** The C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> Fourier coefficients. */
private final FieldFourierCjSjCoefficients<T> fourierCjSj;
/** The maximum value for j index. */
private final int jMax;
/** Constructor.
* @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
* @param fourierCjSj the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
* @param jMax maximum value for j index
* @param field field used by default
*/
FieldUijVijCoefficients(final T[][] currentRhoSigmaj,
final FieldFourierCjSjCoefficients<T> fourierCjSj,
final int jMax,
final Field<T> field) {
this.currentRhoSigmaj = currentRhoSigmaj;
this.fourierCjSj = fourierCjSj;
this.jMax = jMax;
// initialize the internal arrays.
this.u1ij = MathArrays.buildArray(field, 6, 2 * jMax + 1);
this.v1ij = MathArrays.buildArray(field, 6, 2 * jMax + 1);
this.u2ij = MathArrays.buildArray(field, jMax + 1);
this.v2ij = MathArrays.buildArray(field, jMax + 1);
//compute the coefficients
computeU1V1Coefficients(field);
computeU2V2Coefficients(field);
}
/** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
* @param field field used by default
*/
private void computeU1V1Coefficients(final Field<T> field) {
// Zero
final T zero = field.getZero();
// generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
// for j >= 1
// also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
u1ij[0][0] = zero;
for (int j = 1; j <= jMax; j++) {
// compute 1 / j
final double ooj = 1. / j;
for (int i = 0; i < 6; i++) {
//j is aready between 1 and J
u1ij[i][j] = fourierCjSj.getSij(i, j);
v1ij[i][j] = fourierCjSj.getCij(i, j);
// 1 - δ<sub>1j</sub> is 1 for all j > 1
if (j > 1) {
// k starts with 1 because j-J is less than or equal to 0
for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
// C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
// S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, j - kIndex).multiply(currentRhoSigmaj[1][kIndex]).
add(fourierCjSj.getSij(i, j - kIndex).multiply(currentRhoSigmaj[0][kIndex])));
// C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
// S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
v1ij[i][j] = v1ij[i][j].add(fourierCjSj.getCij(i, j - kIndex).multiply(currentRhoSigmaj[0][kIndex]).
subtract(fourierCjSj.getSij(i, j - kIndex).multiply(currentRhoSigmaj[1][kIndex])));
}
}
// since j must be between 1 and J-1 and is already between 1 and J
// the following sum is skiped only for j = jMax
if (j != jMax) {
for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
// -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
// S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, j + kIndex).negate().multiply(currentRhoSigmaj[1][kIndex]).
add(fourierCjSj.getSij(i, j + kIndex).multiply(currentRhoSigmaj[0][kIndex])));
// C<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub> +
// S<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub>
v1ij[i][j] = v1ij[i][j].add(fourierCjSj.getCij(i, j + kIndex).multiply(currentRhoSigmaj[0][kIndex]).
add(fourierCjSj.getSij(i, j + kIndex).multiply(currentRhoSigmaj[1][kIndex])));
}
}
for (int kIndex = 1; kIndex <= jMax; kIndex++) {
// C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
// S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
u1ij[i][j] = u1ij[i][j].add(fourierCjSj.getCij(i, kIndex).negate().multiply(currentRhoSigmaj[1][j + kIndex]).
subtract(fourierCjSj.getSij(i, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])));
// C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
// S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
v1ij[i][j] = v1ij[i][j].add(fourierCjSj.getCij(i, kIndex).multiply(currentRhoSigmaj[0][j + kIndex]).
add(fourierCjSj.getSij(i, kIndex).multiply(currentRhoSigmaj[1][j + kIndex])));
}
// divide by 1 / j
u1ij[i][j] = u1ij[i][j].multiply(-ooj);
v1ij[i][j] = v1ij[i][j].multiply(ooj);
// if index = 1 (i == 0) add the computed terms to U₁⁰
if (i == 0) {
//- (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
u1ij[0][0] = u1ij[0][0].add(u1ij[0][j].negate().multiply(currentRhoSigmaj[0][j]).subtract(v1ij[0][j].multiply(currentRhoSigmaj[1][j])));
}
}
}
// Terms with j > jMax are required only when computing the coefficients
// U₂<sup>j</sup> and V₂<sup>j</sup>
// and those coefficients are only required for Fourier index = 1 (i == 0).
for (int j = jMax + 1; j <= 2 * jMax; j++) {
// compute 1 / j
final double ooj = 1. / j;
//the value of i is 0
u1ij[0][j] = zero;
v1ij[0][j] = zero;
//k starts from j-J as it is always greater than or equal to 1
for (int kIndex = j - jMax; kIndex <= j - 1; kIndex++) {
// C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
// S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
u1ij[0][j] = u1ij[0][j].add(fourierCjSj.getCij(0, j - kIndex).multiply(currentRhoSigmaj[1][kIndex]).
add(fourierCjSj.getSij(0, j - kIndex).multiply(currentRhoSigmaj[0][kIndex])));
// C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
// S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
v1ij[0][j] = v1ij[0][j].add(fourierCjSj.getCij(0, j - kIndex).multiply(currentRhoSigmaj[0][kIndex]).
subtract(fourierCjSj.getSij(0, j - kIndex).multiply(currentRhoSigmaj[1][kIndex])));
}
for (int kIndex = 1; kIndex <= jMax; kIndex++) {
// C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
// S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
u1ij[0][j] = u1ij[0][j].add(fourierCjSj.getCij(0, kIndex).negate().multiply(currentRhoSigmaj[1][j + kIndex]).
subtract(fourierCjSj.getSij(0, kIndex).multiply(currentRhoSigmaj[0][j + kIndex])));
// C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
// S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
v1ij[0][j] = v1ij[0][j].add(fourierCjSj.getCij(0, kIndex).multiply(currentRhoSigmaj[0][j + kIndex]).add(
fourierCjSj.getSij(0, kIndex).multiply(currentRhoSigmaj[1][j + kIndex])));
}
// divide by 1 / j
u1ij[0][j] = u1ij[0][j].multiply(-ooj);
v1ij[0][j] = v1ij[0][j].multiply(ooj);
}
}
/** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
* <p>
* Only the coefficients for Fourier index = 1 (i == 0) are required.
* </p>
* @param field field used by default
*/
private void computeU2V2Coefficients(final Field<T> field) {
for (int j = 1; j <= jMax; j++) {
// compute 1 / j
final double ooj = 1. / j;
// only the values for i == 0 are computed
u2ij[j] = v1ij[0][j];
v2ij[j] = u1ij[0][j];
// 1 - δ<sub>1j</sub> is 1 for all j > 1
if (j > 1) {
for (int l = 1; l <= j - 1; l++) {
// U₁<sup>j-l</sup> * σ<sub>l</sub> +
// V₁<sup>j-l</sup> * ρ<sub>l</sub>
u2ij[j] = u2ij[j].add(u1ij[0][j - l].multiply(currentRhoSigmaj[1][l]).
add(v1ij[0][j - l].multiply(currentRhoSigmaj[0][l])));
// U₁<sup>j-l</sup> * ρ<sub>l</sub> -
// V₁<sup>j-l</sup> * σ<sub>l</sub>
v2ij[j] = v2ij[j].add(u1ij[0][j - l].multiply(currentRhoSigmaj[0][l]).
subtract(v1ij[0][j - l].multiply(currentRhoSigmaj[1][l])));
}
}
for (int l = 1; l <= jMax; l++) {
// -U₁<sup>j+l</sup> * σ<sub>l</sub> +
// U₁<sup>l</sup> * σ<sub>j+l</sub> +
// V₁<sup>j+l</sup> * ρ<sub>l</sub> -
// V₁<sup>l</sup> * ρ<sub>j+l</sub>
u2ij[j] = u2ij[j].add(u1ij[0][j + l].negate().multiply(currentRhoSigmaj[1][l]).
add(u1ij[0][l].multiply(currentRhoSigmaj[1][j + l])).
add(v1ij[0][j + l].multiply(currentRhoSigmaj[0][l])).
subtract(v1ij[0][l].multiply(currentRhoSigmaj[0][j + l])));
// U₁<sup>j+l</sup> * ρ<sub>l</sub> +
// U₁<sup>l</sup> * ρ<sub>j+l</sub> +
// V₁<sup>j+l</sup> * σ<sub>l</sub> +
// V₁<sup>l</sup> * σ<sub>j+l</sub>
u2ij[j] = u2ij[j].add(u1ij[0][j + l].multiply(currentRhoSigmaj[0][l]).
add(u1ij[0][l].multiply(currentRhoSigmaj[0][j + l])).
add(v1ij[0][j + l].multiply(currentRhoSigmaj[1][l])).
add(v1ij[0][l].multiply(currentRhoSigmaj[1][j + l])));
}
// divide by 1 / j
u2ij[j] = u2ij[j].multiply(-ooj);
v2ij[j] = v2ij[j].multiply(ooj);
}
}
/** Get the coefficient U₁<sup>j</sup> for Fourier index i.
*
* @param j j index
* @param i Fourier index (starts at 0)
* @return the coefficient U₁<sup>j</sup> for the given Fourier index i
*/
public T getU1(final int j, final int i) {
return u1ij[i][j];
}
/** Get the coefficient V₁<sup>j</sup> for Fourier index i.
*
* @param j j index
* @param i Fourier index (starts at 0)
* @return the coefficient V₁<sup>j</sup> for the given Fourier index i
*/
public T getV1(final int j, final int i) {
return v1ij[i][j];
}
/** Get the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0).
*
* @param j j index
* @return the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0)
*/
public T getU2(final int j) {
return u2ij[j];
}
/** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
*
* @param j j index
* @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
*/
public T getV2(final int j) {
return v2ij[j];
}
}
/** Coefficients valid for one time slot. */
private static class Slot {
/**The coefficients D<sub>i</sub><sup>j</sup>.
* <p>
* Only for j = 1 and j = 2 the coefficients are not 0. <br>
* 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[] dij;
/** The coefficients C<sub>i</sub><sup>j</sup>.
* <p>
* The index order is cij[j][i] <br/>
* i corresponds to the equinoctial element, as follows: <br/>
* - i=0 for a <br/>
* - i=1 for k <br/>
* - i=2 for h <br/>
* - i=3 for q <br/>
* - i=4 for p <br/>
* - i=5 for λ <br/>
* </p>
*/
private final ShortPeriodicsInterpolatedCoefficient[] cij;
/** The coefficients S<sub>i</sub><sup>j</sup>.
* <p>
* The index order is sij[j][i] <br/>
* i corresponds to the equinoctial element, as follows: <br/>
* - i=0 for a <br/>
* - i=1 for k <br/>
* - i=2 for h <br/>
* - i=3 for q <br/>
* - i=4 for p <br/>
* - i=5 for λ <br/>
* </p>
*/
private final ShortPeriodicsInterpolatedCoefficient[] sij;
/** Simple constructor.
* @param jMax maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
*/
Slot(final int jMax, final int interpolationPoints) {
dij = new ShortPeriodicsInterpolatedCoefficient[3];
cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
// Initialize the C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup> and D<sub>i</sub><sup>j</sup> coefficients
for (int j = 0; j <= jMax; j++) {
cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
if (j > 0) {
sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
}
// Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
if (j == 1 || j == 2) {
dij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
}
}
}
}
/** Coefficients valid for one time slot. */
private static class FieldSlot <T extends RealFieldElement<T>> {
/**The coefficients D<sub>i</sub><sup>j</sup>.
* <p>
* Only for j = 1 and j = 2 the coefficients are not 0. <br>
* 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>[] dij;
/** The coefficients C<sub>i</sub><sup>j</sup>.
* <p>
* The index order is cij[j][i] <br/>
* i corresponds to the equinoctial element, as follows: <br/>
* - i=0 for a <br/>
* - i=1 for k <br/>
* - i=2 for h <br/>
* - i=3 for q <br/>
* - i=4 for p <br/>
* - i=5 for λ <br/>
* </p>
*/
private final 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 jMax maximum value for j index
* @param interpolationPoints number of points used in the interpolation process
*/
@SuppressWarnings("unchecked")
FieldSlot(final int jMax, final int interpolationPoints) {
dij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, 3);
cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
// Initialize the C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup> and D<sub>i</sub><sup>j</sup> coefficients
for (int j = 0; j <= jMax; j++) {
cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
if (j > 0) {
sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
}
// Initialize only the non-zero D<sub>i</sub><sup>j</sup> coefficients
if (j == 1 || j == 2) {
dij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
}
}
}
}
}