From 059ac22f8fe47a7da0903c31011d0035d7450c5d Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Fri, 17 Jan 2014 11:03:00 +0100
Subject: [PATCH] Compute Hansen coefficients using linear transformation.
Previous implementation was based on recursion.
This implementation is based on theoretical work done by Petre Bazavan.
---
pom.xml | 6 +
.../dsst/forces/DSSTThirdBody.java | 191 +-------
.../dsst/forces/TesseralContribution.java | 287 ++++-------
.../dsst/forces/ZonalContribution.java | 173 +------
.../hansen/HansenTesseralLinear.java | 451 ++++++++++++++++++
.../hansen/HansenThirdBodyLinear.java | 322 +++++++++++++
.../utilities/hansen/HansenUtilities.java | 160 +++++++
.../utilities/hansen/HansenZonalLinear.java | 250 ++++++++++
.../hansen/PolynomialFunctionMatrix.java | 146 ++++++
9 files changed, 1451 insertions(+), 535 deletions(-)
create mode 100644 src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenTesseralLinear.java
create mode 100644 src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenThirdBodyLinear.java
create mode 100644 src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenUtilities.java
create mode 100644 src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenZonalLinear.java
create mode 100644 src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/PolynomialFunctionMatrix.java
diff --git a/pom.xml b/pom.xml
index fa37cf8261..06586b2fe0 100644
--- a/pom.xml
+++ b/pom.xml
@@ -124,6 +124,12 @@
<contributor>
<name>Daniel Aguilar Taboada</name>
</contributor>
+ <contributor>
+ <name>Lucian Barbulescu</name>
+ </contributor>
+ <contributor>
+ <name>Petre Bazavan</name>
+ </contributor>
<contributor>
<name>Nicolas Bernard</name>
</contributor>
diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/DSSTThirdBody.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/DSSTThirdBody.java
index 8115e30d01..d1f35cdfcb 100644
--- a/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/DSSTThirdBody.java
+++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/DSSTThirdBody.java
@@ -28,6 +28,7 @@ import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
+import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenThirdBodyLinear;
import org.orekit.time.AbsoluteDate;
/** Third body attraction perturbation to the
@@ -86,6 +87,8 @@ public class DSSTThirdBody implements DSSTForceModel {
private double gamma;
// Common factors for potential computation
+ /** B = sqrt(1 - e<sup>2</sup>). */
+ private double B;
/** Χ = 1 / sqrt(1 - e<sup>2</sup>) = 1 / B. */
private double X;
/** Χ<sup>2</sup>. */
@@ -112,6 +115,10 @@ public class DSSTThirdBody implements DSSTForceModel {
/** Maw power for e in the serie expansion. */
private int maxEccPow;
+ /** An array that contains the objects needed to build the Hansen coefficients. <br/>
+ * The index is s */
+ private HansenThirdBodyLinear[] hansenObjects;
+
/** Complete constructor.
* @param body the 3rd body to consider
* @see org.orekit.bodies.CelestialBodyFactory
@@ -132,6 +139,13 @@ public class DSSTThirdBody implements DSSTForceModel {
for (int i = 1; i < dim; i++) {
fact[i] = i * fact[i - 1];
}
+
+ //Initialise the HansenCoefficient generator
+ this.hansenObjects = new HansenThirdBodyLinear[MAX_POWER + 1];
+ for (int s = 0; s <= MAX_POWER; s++) {
+ this.hansenObjects[s] = new HansenThirdBodyLinear(MAX_POWER, s);
+ }
+
}
/** Get third body.
@@ -247,7 +261,7 @@ public class DSSTThirdBody implements DSSTForceModel {
// Equinoctial coefficients
final double A = aux.getA();
- final double B = aux.getB();
+ B = aux.getB();
final double C = aux.getC();
// Χ
@@ -317,8 +331,6 @@ public class DSSTThirdBody implements DSSTForceModel {
*/
private double[] computeUDerivatives() throws OrekitException {
- // Hansen coefficients
- final HansenThirdBody hansen = new HansenThirdBody();
// Gs and Hs coefficients
final double[][] GsHs = CoefficientsFactory.computeGsHs(k, h, alpha, beta, maxEccPow);
// Qns coefficients
@@ -366,8 +378,9 @@ public class DSSTThirdBody implements DSSTForceModel {
// (n - s) must be even
if ((n - s) % 2 == 0) {
// Extract data from previous computation :
- final double kns = hansen.getValue(n, s);
- final double dkns = hansen.getDerivative(n, s);
+ final double kns = this.hansenObjects[s].getValue(n, B);
+ final double dkns = this.hansenObjects[s].getDerivative(n, B);
+
final double vns = Vns.get(new NSKey(n, s));
final double coef0 = delta0s * aoR3Pow[n] * vns;
final double coef1 = coef0 * Qns[n][s];
@@ -405,172 +418,4 @@ public class DSSTThirdBody implements DSSTForceModel {
};
}
-
- /** Hansen coefficients for 3rd body force model.
- * <p>
- * Hansen coefficients are functions of the eccentricity.
- * For a given eccentricity, the computed elements are stored in a map.
- * </p>
- */
- private class HansenThirdBody {
-
- /** Map to store every Hansen coefficient computed. */
- private TreeMap<NSKey, Double> coefficients;
-
- /** Map to store every Hansen coefficient derivative computed. */
- private TreeMap<NSKey, Double> derivatives;
-
- /** Simple constructor. */
- public HansenThirdBody() {
- coefficients = new TreeMap<CoefficientsFactory.NSKey, Double>();
- derivatives = new TreeMap<CoefficientsFactory.NSKey, Double>();
- initialize();
- }
-
- /** Get the K<sub>0</sub><sup>n,s</sup> coefficient value.
- * @param n n value
- * @param s s value
- * @return K<sub>0</sub><sup>n,s</sup>
- */
- public double getValue(final int n, final int s) {
- if (coefficients.containsKey(new NSKey(n, s))) {
- return coefficients.get(new NSKey(n, s));
- } else {
- // Compute the K0(n,s) coefficients
- return computeValue(n, s);
- }
- }
-
- /** Get the dK<sub>0</sub><sup>n,s</sup> / d&x; coefficient derivative.
- * @param n n value
- * @param s s value
- * @return dK<sub>j</sub><sup>n,s</sup> / d&x;
- */
- public double getDerivative(final int n, final int s) {
- if (derivatives.containsKey(new NSKey(n, s))) {
- return derivatives.get(new NSKey(n, s));
- } else {
- // Compute the dK0(n,s) / dX derivative
- return computeDerivative(n, s);
- }
- }
-
- /** Initialization. */
- private void initialize() {
- final double ec2 = ecc * ecc;
- final double oX3 = 1. / XXX;
- coefficients.put(new NSKey(0, 0), 1.);
- coefficients.put(new NSKey(0, 1), -1.);
- coefficients.put(new NSKey(1, 0), 1. + 0.5 * ec2);
- coefficients.put(new NSKey(1, 1), -1.5);
- coefficients.put(new NSKey(2, 0), 1. + 1.5 * ec2);
- coefficients.put(new NSKey(2, 1), -2. - 0.5 * ec2);
- derivatives.put(new NSKey(0, 0), 0.);
- derivatives.put(new NSKey(1, 0), oX3);
- derivatives.put(new NSKey(2, 0), 3. * oX3);
- derivatives.put(new NSKey(2, 1), -oX3);
- }
-
- /** Compute K<sub>0</sub><sup>n,s</sup> from Equation 2.7.3-(7)(8).
- * @param n n value
- * @param s s value
- * @return K<sub>0</sub><sup>n,s</sup>
- */
- private double computeValue(final int n, final int s) {
- // Initialize return value
- double kns = 0.;
-
- if (n == (s - 1)) {
-
- final NSKey key = new NSKey(s - 2, s - 1);
- if (coefficients.containsKey(key)) {
- kns = coefficients.get(key);
- } else {
- kns = computeValue(s - 2, s - 1);
- }
-
- kns *= -(2. * s - 1.) / s;
-
- } else if (n == s) {
-
- final NSKey key = new NSKey(s - 1, s);
- if (coefficients.containsKey(key)) {
- kns = coefficients.get(key);
- } else {
- kns = computeValue(s - 1, s);
- }
-
- kns *= (2. * s + 1.) / (s + 1.);
-
- } else if (n > s) {
-
- final NSKey key1 = new NSKey(n - 1, s);
- double knM1 = 0.;
- if (coefficients.containsKey(key1)) {
- knM1 = coefficients.get(key1);
- } else {
- knM1 = computeValue(n - 1, s);
- }
-
- final NSKey key2 = new NSKey(n - 2, s);
- double knM2 = 0.;
- if (coefficients.containsKey(key2)) {
- knM2 = coefficients.get(key2);
- } else {
- knM2 = computeValue(n - 2, s);
- }
-
- final double val1 = (2. * n + 1.) / (n + 1.);
- final double val2 = (n + s) * (n - s) / (n * (n + 1.) * XX);
-
- kns = val1 * knM1 - val2 * knM2;
- }
-
- coefficients.put(new NSKey(n, s), kns);
- return kns;
- }
-
- /** Compute dK<sub>0</sub><sup>n,s</sup> / d&x; from Equation 3.2-(3).
- * @param n n value
- * @param s s value
- * @return dK<sub>0</sub><sup>n,s</sup> / d&x;
- */
- private double computeDerivative(final int n, final int s) {
- // Initialize return value
- double dknsdx = 0.;
-
- if (n > s) {
-
- final NSKey keyNm1 = new NSKey(n - 1, s);
- double dKnM1 = 0.;
- if (derivatives.containsKey(keyNm1)) {
- dKnM1 = derivatives.get(keyNm1);
- } else {
- dKnM1 = computeDerivative(n - 1, s);
- }
-
- final NSKey keyNm2 = new NSKey(n - 2, s);
- double dKnM2 = 0.;
- if (derivatives.containsKey(keyNm2)) {
- dKnM2 = derivatives.get(keyNm2);
- } else {
- dKnM2 = computeDerivative(n - 2, s);
- }
-
- final double knM2 = getValue(n - 2, s);
-
- final double val1 = (2. * n + 1.) / (n + 1.);
- final double coef = (n + s) * (n - s) / (n * (n + 1.));
- final double val2 = coef / XX;
- final double val3 = 2. * coef / XXX;
-
- dknsdx = val1 * dKnM1 - val2 * dKnM2 + val3 * knM2;
- }
-
- derivatives.put(new NSKey(n, s), dknsdx);
- return dknsdx;
- }
-
- }
-
}
diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/TesseralContribution.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/TesseralContribution.java
index a1964f21f8..fe1473dc32 100644
--- a/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/TesseralContribution.java
+++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/TesseralContribution.java
@@ -18,7 +18,6 @@ package org.orekit.propagation.semianalytical.dsst.forces;
import java.util.ArrayList;
import java.util.List;
-import java.util.TreeMap;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
@@ -33,11 +32,10 @@ import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
-import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.MNSKey;
import org.orekit.propagation.semianalytical.dsst.utilities.GHmsjPolynomials;
import org.orekit.propagation.semianalytical.dsst.utilities.GammaMnsFunction;
import org.orekit.propagation.semianalytical.dsst.utilities.JacobiPolynomials;
-import org.orekit.propagation.semianalytical.dsst.utilities.NewcombOperators;
+import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenTesseralLinear;
import org.orekit.time.AbsoluteDate;
/** Tesseral contribution to the {@link DSSTCentralBody central body gravitational
@@ -97,23 +95,36 @@ class TesseralContribution implements DSSTForceModel {
// Equinoctial elements (according to DSST notation)
/** a. */
private double a;
+
/** ex. */
private double k;
+
/** ey. */
private double h;
+
/** hx. */
private double q;
+
/** hy. */
private double p;
+
/** lm. */
private double lm;
/** Eccentricity. */
private double ecc;
+ // Common factors for potential computation
+ /** Χ = 1 / sqrt(1 - e<sup>2</sup>) = 1 / B. */
+ private double chi;
+
+ /** Χ<sup>2</sup>. */
+ private double chi2;
+
// Equinoctial reference frame vectors (according to DSST notation)
/** Equinoctial frame f vector. */
private Vector3D f;
+
/** Equinoctial frame g vector. */
private Vector3D g;
@@ -122,27 +133,39 @@ class TesseralContribution implements DSSTForceModel {
/** Direction cosine α. */
private double alpha;
+
/** Direction cosine β. */
private double beta;
+
/** Direction cosine γ. */
private double gamma;
// Common factors from equinoctial coefficients
/** 2 * a / A .*/
private double ax2oA;
+
/** 1 / (A * B) .*/
private double ooAB;
+
/** B / A .*/
private double BoA;
+
/** B / (A * (1 + B)) .*/
private double BoABpo;
+
/** C / (2 * A * B) .*/
private double Co2AB;
+
/** μ / a .*/
private double moa;
+
/** R / a .*/
private double roa;
+ /** A two dimensional array that contains the objects needed to build the Hansen coefficients. <br/>
+ * The indexes are s + maxDegree and j */
+ private HansenTesseralLinear[][] hansenObjects;
+
/** Single constructor.
* @param centralBodyFrame rotating body frame
* @param centralBodyRotationRate central body rotation rate (rad/s)
@@ -174,7 +197,6 @@ class TesseralContribution implements DSSTForceModel {
for (int i = 1; i < maxFact; i++) {
fact[i] = i * fact[i - 1];
}
-
}
/** {@inheritDoc} */
@@ -213,6 +235,45 @@ class TesseralContribution implements DSSTForceModel {
// Set the maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
maxHansen = maxEccPow / 2;
+ //initialise the HansenTesseralLinear objects needed
+ if (resOrders.size() > 0) {
+ initialiseHansenObjects();
+ }
+ }
+
+ /**
+ * Create the objects needed for linear transformation.
+ */
+ private void initialiseHansenObjects() {
+ //compute the j maximum
+ final int jMax = FastMath.max(1, (int) FastMath.round(ratio * resOrders.get(resOrders.size() - 1)));
+
+ //Allocate the two dimensional array
+ this.hansenObjects = new HansenTesseralLinear[2 * maxDegree + 1][jMax + 1];
+
+ //loop through the resonant terms
+ for (int m : resOrders) {
+ //Compute the corresponding j term
+ final int j = FastMath.max(1, (int) FastMath.round(ratio * m));
+
+ //Compute the sMin and sMax values
+ final int sMin = FastMath.min(maxEccPow - j, maxDegree);
+ final int sMax = FastMath.min(maxEccPow + j, maxDegree);
+
+ //loop through the s values
+ for (int s = 0; s <= sMax; s++) {
+ //Compute the n0 value
+ final int n0 = FastMath.max(FastMath.max(2, m), s);
+
+ //Create the object for the pair j,s
+ this.hansenObjects[s + maxDegree][j] = new HansenTesseralLinear(maxDegree, s, j, n0);
+
+ if (s > 0 && s <= sMin) {
+ //Also create the object for the pair j, -s
+ this.hansenObjects[maxDegree - s][j] = new HansenTesseralLinear(maxDegree, -s, j, n0);
+ }
+ }
+ }
}
/** {@inheritDoc} */
@@ -271,6 +332,9 @@ class TesseralContribution implements DSSTForceModel {
// R / a
roa = provider.getAe() / a;
+ // Χ = 1 / B
+ chi = 1. / B;
+ chi2 = chi * chi;
}
/** {@inheritDoc} */
@@ -374,9 +438,6 @@ class TesseralContribution implements DSSTForceModel {
// GAMMAmns function
final GammaMnsFunction gammaMNS = new GammaMnsFunction(fact, gamma, I);
- // Hansen coefficients
- final HansenTesseral hansen = new HansenTesseral(ecc, maxHansen);
-
// R / a up to power degree
final double[] roaPow = new double[maxDegree + 1];
roaPow[0] = 1.;
@@ -418,7 +479,7 @@ class TesseralContribution implements DSSTForceModel {
// n-SUM for s positive
final double[][] nSumSpos = computeNSum(date, j, m, s,
- roaPow, ghMSJ, gammaMNS, hansen);
+ roaPow, ghMSJ, gammaMNS);
dUdaCos += nSumSpos[0][0];
dUdaSin += nSumSpos[0][1];
dUdhCos += nSumSpos[1][0];
@@ -437,7 +498,7 @@ class TesseralContribution implements DSSTForceModel {
// n-SUM for s negative
if (s > 0 && s <= sMin) {
final double[][] nSumSneg = computeNSum(date, j, m, -s,
- roaPow, ghMSJ, gammaMNS, hansen);
+ roaPow, ghMSJ, gammaMNS);
dUdaCos += nSumSneg[0][0];
dUdaSin += nSumSneg[0][1];
dUdhCos += nSumSneg[1][0];
@@ -485,16 +546,13 @@ class TesseralContribution implements DSSTForceModel {
* @param roaPow powers of R/a up to degree <i>n</i>
* @param ghMSJ G<sup>j</sup><sub>m,s</sub> and H<sup>j</sup><sub>m,s</sub> polynomials
* @param gammaMNS Γ<sup>m</sup><sub>n,s</sub>(γ) function
- * @param hansen Hansen coefficients
* @return Components of U<sub>n</sub> derivatives for fixed j, m, s
* @throws OrekitException if some error occurred
*/
private double[][] computeNSum(final AbsoluteDate date,
- final int j, final int m, final int s,
- final double[] roaPow,
- final GHmsjPolynomials ghMSJ,
- final GammaMnsFunction gammaMNS,
- final HansenTesseral hansen) throws OrekitException {
+ final int j, final int m, final int s, final double[] roaPow,
+ final GHmsjPolynomials ghMSJ, final GammaMnsFunction gammaMNS)
+ throws OrekitException {
//spherical harmonics
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
@@ -524,23 +582,15 @@ class TesseralContribution implements DSSTForceModel {
// Initialise lower degree nmin = (Max(2, m, |s|)) for summation over n
final int nmin = FastMath.max(FastMath.max(2, m), FastMath.abs(s));
- // n max value for computing Hansen kernel from Newcomb operators
- final int nmax = nmin + 3;
- // n-SUM from nmin to N
- for (int n = nmin; n <= maxDegree; n++) {
+ //Get the corresponding Hansen object
+ final HansenTesseralLinear hans = this.hansenObjects[maxDegree + s][j];
- // Compute Hansen kernel values and derivatives
- if (n <= nmax) {
- // from Newcomb operators
- hansen.valueFromNewcomb(j, -n - 1, s);
- hansen.derivFromNewcomb(j, -n - 1, s);
- } else {
- // from recurrence relations
- hansen.valueFromRecurrence(j, -n - 1, s);
- hansen.derivFromRecurrence(j, -n - 1, s);
- }
+ //Compute the initial values using newComb operators
+ hans.computeInitValues(ecc * ecc, chi, chi2, maxHansen);
+ // n-SUM from nmin to N
+ for (int n = nmin; n <= maxDegree; n++) {
// If (n - s) is odd, the contribution is null because of Vmns
if ((n - s) % 2 == 0) {
@@ -553,8 +603,8 @@ class TesseralContribution implements DSSTForceModel {
final double dGaMNS = gammaMNS.getDerivative(m, n, s);
// Hansen kernel value and derivative
- final double kJNS = hansen.getValue(j, -n - 1, s);
- final double dkJNS = hansen.getDeriv(j, -n - 1, s);
+ final double kJNS = hans.getValue(-n - 1, chi);
+ final double dkJNS = hans.getDerivative(-n - 1, chi);
// Gjms, Hjms polynomials and derivatives
final double gMSJ = ghMSJ.getGmsj(m, s, j);
@@ -629,179 +679,4 @@ class TesseralContribution implements DSSTForceModel {
{dUdGaCos, dUdGaSin}};
}
- /** Hansen coefficients for tesseral contribution to central body force model.
- * <p>
- * Hansen coefficients are functions of the eccentricity.
- * For a given eccentricity, the computed elements are stored in a map.
- * </p>
- */
- private static class HansenTesseral {
-
- /** Map to store every Hansen kernel value computed. */
- private TreeMap<MNSKey, Double> values;
-
- /** Map to store every Hansen kernel derivative computed. */
- private TreeMap<MNSKey, Double> derivatives;
-
- /** Eccentricity. */
- private final double e2;
-
- /** 1 - e<sup>2</sup>. */
- private final double ome2;
-
- /** χ = 1 / sqrt(1- e<sup>2</sup>). */
- private final double chi;
-
- /** χ<sup>2</sup> = 1 / (1- e<sup>2</sup>). */
- private final double chi2;
-
- /** dχ / de<sup>2</sup> = χ<sup>3</sup> / 2. */
- private final double dchi;
-
- /** dχ<sup>2</sup> / de<sup>2</sup> = χ<sup>4</sup>. */
- private final double dchi2;
-
- /** Max power of e<sup>2</sup> in serie expansion
- * using Newcomb operator for Hansen kernel computation.
- */
- private final int maxNewcomb;
-
- /** Simple constructor.
- * @param ecc eccentricity
- * @param maxHansen maximum power of e<sup>2</sup> to use in series expansion for the Hansen coefficient
- */
- public HansenTesseral(final double ecc, final int maxHansen) {
- this.values = new TreeMap<CoefficientsFactory.MNSKey, Double>();
- this.derivatives = new TreeMap<CoefficientsFactory.MNSKey, Double>();
- this.maxNewcomb = maxHansen;
- this.e2 = ecc * ecc;
- this.ome2 = 1. - e2;
- this.chi = 1. / FastMath.sqrt(ome2);
- this.chi2 = chi * chi;
- this.dchi = 0.5 * chi * chi2;
- this.dchi2 = chi2 * chi2;
- }
-
- /** Get the K<sub>j</sub><sup>-n-1,s</sup> value.
- * @param j j value
- * @param mnm1 -n-1 value
- * @param s s value
- * @return K<sub>j</sub><sup>-n-1,s</sup>
- */
- public double getValue(final int j, final int mnm1, final int s) {
- return values.get(new MNSKey(j, mnm1, s));
- }
-
- /** Get the dK<sub>j</sub><sup>-n-1,s</sup> / de<sup>2</sup> derivative.
- * @param j j value
- * @param mnm1 -n-1 value
- * @param s s value
- * @return dK<sub>j</sub><sup>-n-1,s</sup> / de<sup>2</sup>
- */
- public double getDeriv(final int j, final int mnm1, final int s) {
- return derivatives.get(new MNSKey(j, mnm1, s));
- }
-
- /** Compute the K<sub>j</sub><sup>-n-1,s</sup> from equation 2.7.3-(10).<br>
- * <p>
- * The coefficient value is evaluated from the {@link NewcombOperators} elements.
- * </p>
- * @param j j value
- * @param mnm1 -n-1 value
- * @param s s value
- */
- public void valueFromNewcomb(final int j, final int mnm1, final int s) {
- // Initialization
- final int aHT = FastMath.max(j - s, 0);
- final int bHT = FastMath.max(s - j, 0);
- // Expansion until maxNewcomb, the maximum power in e^2 for the Kernel value
- double sum = NewcombOperators.getValue(maxNewcomb + aHT, maxNewcomb + bHT, mnm1, s);
- for (int alphaHT = maxNewcomb - 1; alphaHT >= 0; alphaHT--) {
- sum *= e2;
- sum += NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
- }
- // Kernel value from equation 2.7.3-(10)
- final double value = FastMath.pow(chi2, -mnm1 - 1) * sum / chi;
- // Storage of the Kernel value in the map
- values.put(new MNSKey(j, mnm1, s), value);
- }
-
- /** Compute dK<sub>j</sub><sup>-n-1,s</sup>/de<sup>2</sup> from equation 3.3-(5).
- * <p>
- * The derivative value is evaluated from the {@link NewcombOperators} elements.
- * </p>
- * @param j j value
- * @param mnm1 -n-1 value
- * @param s s value
- */
- public void derivFromNewcomb(final int j, final int mnm1, final int s) {
- // Initialization
- final int aHT = FastMath.max(j - s, 0);
- final int bHT = FastMath.max(s - j, 0);
- // Expansion until maxNewcomb-1, the maximum power in e^2 for the Kernel derivative
- double sum = maxNewcomb * NewcombOperators.getValue(maxNewcomb + aHT, maxNewcomb + bHT, mnm1, s);
- for (int alphaHT = maxNewcomb - 1; alphaHT >= 1; alphaHT--) {
- sum *= e2;
- sum += alphaHT * NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
- }
- // Kernel derivative from equation 3.3-(5)
- final MNSKey key = new MNSKey(j, mnm1, s);
- final double coef = -(mnm1 + 1.5);
- final double Kjns = values.get(key);
- final double derivative = coef * chi2 * Kjns + FastMath.pow(chi2, -mnm1 - 1) * sum / chi;
- // Storage of the Kernel derivative in the map
- derivatives.put(new MNSKey(j, mnm1, s), derivative);
- }
-
- /** Compute the K<sub>j</sub><sup>-n-1,s</sup> coefficient from equation 2.7.3-(9).
- *
- * @param j j value
- * @param mnm1 -n-1 value
- * @param s s value
- */
- public void valueFromRecurrence(final int j, final int mnm1, final int s) {
- final int n = -mnm1 - 1;
- final double kmn = values.get(new MNSKey(j, -n, s));
- final double kmnp1 = values.get(new MNSKey(j, -n + 1, s));
- final double kmnp3 = values.get(new MNSKey(j, -n + 3, s));
- final double den = (3 - n) * (1 - n + s) * (1 - n - s);
- final double ck = chi2 / den;
- final double ckmn = (3 - n) * (1 - n) * (3 - 2 * n);
- final double ckmnp1 = (2 - n) * ((3 - n) * (1 - n) + (2 * j * s) / chi);
- final double ckmnp3 = j * j * (1 - n);
- final double value = ck * (ckmn * kmn - ckmnp1 * kmnp1 + ckmnp3 * kmnp3);
- // Store value
- values.put(new MNSKey(j, mnm1, s), value);
- }
-
- /** Compute dK<sub>j</sub><sup>-n-1,s</sup>/de<sup>2</sup> from derivation of equation 2.7.3-(9).
- *
- * @param j j value
- * @param mnm1 -n-1 value
- * @param s s value
- */
- public void derivFromRecurrence(final int j, final int mnm1, final int s) {
- final int n = -mnm1 - 1;
- final MNSKey keymn = new MNSKey(j, -n, s);
- final MNSKey keymnp1 = new MNSKey(j, -n + 1, s);
- final MNSKey keymnp3 = new MNSKey(j, -n + 3, s);
- final double kmn = values.get(keymn);
- final double kmnp1 = values.get(keymnp1);
- final double kmnp3 = values.get(keymnp3);
- final double dkmn = derivatives.get(keymn);
- final double dkmnp1 = derivatives.get(keymnp1);
- final double dkmnp3 = derivatives.get(keymnp3);
- final double den = (3 - n) * (1 - n + s) * (1 - n - s);
- final double cdkmn = (3 - n) * (1 - n) * (3 - 2 * n);
- final double c0dkmnp1 = (n - 2) * (3 - n) * (1 - n);
- final double c1dkmnp1 = (n - 2) * (2 * j * s);
- final double cdkmnp3 = j * j * (1 - n);
- final double deriv = chi2 * (cdkmn * dkmn + c0dkmnp1 * dkmnp1 + cdkmnp3 * dkmnp3) +
- chi * c1dkmnp1 * dkmnp1 +
- dchi2 * (cdkmn * kmn + c0dkmnp1 * kmnp1 + cdkmnp3 * kmnp3) +
- dchi * c1dkmnp1 * kmnp1;
- // Store derivative
- derivatives.put(new MNSKey(j, mnm1, s), deriv / den);
- }
- }
}
diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/ZonalContribution.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/ZonalContribution.java
index d276070bcf..d1e1c7f78f 100644
--- a/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/ZonalContribution.java
+++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/forces/ZonalContribution.java
@@ -28,6 +28,7 @@ import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
import org.orekit.propagation.semianalytical.dsst.utilities.UpperBounds;
+import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenZonalLinear;
import org.orekit.time.AbsoluteDate;
/** Zonal contribution to the {@link DSSTCentralBody central body gravitational perturbation}.
@@ -102,6 +103,10 @@ class ZonalContribution implements DSSTForceModel {
/** μ / a .*/
private double muoa;
+ /** An array that contains the objects needed to build the Hansen coefficients. <br/>
+ * The index is s*/
+ private HansenZonalLinear[] hansenObjects;
+
/** Simple constructor.
* @param provider provider for spherical harmonics
*/
@@ -124,6 +129,13 @@ class ZonalContribution implements DSSTForceModel {
// Initialize default values
this.maxEccPow = (maxDegree == 2) ? 0 : Integer.MIN_VALUE;
+
+ //Initialise the HansenCoefficient generator
+ this.hansenObjects = new HansenZonalLinear[maxDegree - 1];
+
+ for (int s = 0; s < maxDegree - 1; s++) {
+ this.hansenObjects[s] = new HansenZonalLinear(maxDegree, s);
+ }
}
/** Get the spherical harmonics provider.
@@ -340,8 +352,6 @@ class ZonalContribution implements DSSTForceModel {
final UnnormalizedSphericalHarmonics harmonics = provider.onDate(date);
- // Hansen coefficients
- final HansenZonal hansen = new HansenZonal();
// Gs and Hs coefficients
final double[][] GsHs = CoefficientsFactory.computeGsHs(k, h, alpha, beta, maxEccPow);
// Qns coefficients
@@ -389,9 +399,11 @@ class ZonalContribution implements DSSTForceModel {
for (int n = s + 2; n <= maxDegree; n++) {
// (n - s) must be even
if ((n - s) % 2 == 0) {
- // Extract data from previous computation :
- final double kns = hansen.getValue(-n - 1, s);
- final double dkns = hansen.getDerivative(-n - 1, s);
+
+ //Extract data from previous computation :
+ final double kns = this.hansenObjects[s].getValue(-n - 1, X);
+ final double dkns = this.hansenObjects[s].getDerivative(-n - 1, X);
+
final double vns = Vns.get(new NSKey(n, s));
final double coef0 = d0s * roaPow[n] * vns * -harmonics.getUnnormalizedCnm(n, 0);
final double coef1 = coef0 * Qns[n][s];
@@ -424,155 +436,4 @@ class ZonalContribution implements DSSTForceModel {
dUdGa * -muoa
};
}
-
- /** Hansen coefficients for zonal contribution to central body force model.
- * <p>
- * Hansen coefficients are functions of the eccentricity.
- * For a given eccentricity, the computed elements are stored in a map.
- * </p>
- */
- private class HansenZonal {
-
- /** Map to store every Hansen coefficient computed. */
- private TreeMap<NSKey, Double> coefficients;
-
- /** Map to store every Hansen coefficient derivative computed. */
- private TreeMap<NSKey, Double> derivatives;
-
- /** Simple constructor. */
- public HansenZonal() {
- coefficients = new TreeMap<CoefficientsFactory.NSKey, Double>();
- derivatives = new TreeMap<CoefficientsFactory.NSKey, Double>();
- initialize();
- }
-
- /** Get the K<sub>0</sub><sup>-n-1,s</sup> coefficient value.
- * @param mnm1 -n-1 value
- * @param s s value
- * @return K<sub>0</sub><sup>-n-1,s</sup>
- */
- public double getValue(final int mnm1, final int s) {
- if (coefficients.containsKey(new NSKey(mnm1, s))) {
- return coefficients.get(new NSKey(mnm1, s));
- } else {
- // Compute the K0(-n-1,s) coefficients
- return computeValue(-mnm1 - 1, s);
- }
- }
-
- /** Get the dK<sub>0</sub><sup>-n-1,s</sup> / dχ coefficient derivative.
- * @param mnm1 -n-1 value
- * @param s s value
- * @return dK<sub>0</sub><sup>-n-1,s</sup> / dχ
- */
- public double getDerivative(final int mnm1, final int s) {
- if (derivatives.containsKey(new NSKey(mnm1, s))) {
- return derivatives.get(new NSKey(mnm1, s));
- } else {
- // Compute the dK0(-n-1,s) / dX derivative
- return computeDerivative(-mnm1 - 1, s);
- }
- }
-
- /** Kernels initialization. */
- private void initialize() {
- // coefficients
-// coefficients.put(new NSKey(0, 0), 1.);
- coefficients.put(new NSKey(-1, 0), 0.);
- coefficients.put(new NSKey(-1, 1), 0.);
- coefficients.put(new NSKey(-2, 0), X);
- coefficients.put(new NSKey(-2, 1), 0.);
- coefficients.put(new NSKey(-3, 0), XXX);
- coefficients.put(new NSKey(-3, 1), 0.5 * XXX);
- // derivatives
- derivatives.put(new NSKey(-1, 0), 0.);
- derivatives.put(new NSKey(-2, 0), 1.);
- derivatives.put(new NSKey(-2, 1), 0.);
- derivatives.put(new NSKey(-3, 1), 1.5 * XX);
- }
-
- /** Compute the K<sub>0</sub><sup>-n-1,s</sup> coefficient from equation 2.7.3-(6).
- * @param n n positive value. For a given 'n', the K<sub>0</sub><sup>-n-1,s</sup> is returned.
- * @param s s value
- * @return K<sub>0</sub><sup>-n-1,s</sup>
- */
- private double computeValue(final int n, final int s) {
- // Initialize return value
- double kns = 0.;
-
- final NSKey key = new NSKey(-n - 1, s);
-
- if (coefficients.containsKey(key)) {
- kns = coefficients.get(key);
- } else {
- if (n == s) {
- kns = 0.;
- } else if (n == (s + 1)) {
- kns = FastMath.pow(X, 1 + 2 * s) / FastMath.pow(2, s);
- } else {
- final NSKey keymNS = new NSKey(-n, s);
- double kmNS = 0.;
- if (coefficients.containsKey(keymNS)) {
- kmNS = coefficients.get(keymNS);
- } else {
- kmNS = computeValue(n - 1, s);
- }
-
- final NSKey keymNp1S = new NSKey(-(n - 1), s);
- double kmNp1S = 0.;
- if (coefficients.containsKey(keymNp1S)) {
- kmNp1S = coefficients.get(keymNp1S);
- } else {
- kmNp1S = computeValue(n - 2, s);
- }
-
- kns = (n - 1.) * XX * ((2. * n - 3.) * kmNS - (n - 2.) * kmNp1S);
- kns /= (n + s - 1.) * (n - s - 1.);
- }
- // Add K(-n-1, s)
- coefficients.put(key, kns);
- }
- return kns;
- }
-
- /** Compute dK<sub>0</sub><sup>-n-1,s</sup> / dχ from Equation 3.1-(7).
- * @param n n positive value. For a given 'n', the dK<sub>0</sub><sup>-n-1,s</sup> / dχ is returned.
- * @param s s value
- * @return dK<sub>0</sub><sup>-n-1,s</sup> / dχ
- */
- private double computeDerivative(final int n, final int s) {
- // Initialize return value
- double dkdxns = 0.;
-
- final NSKey key = new NSKey(-n - 1, s);
- if (n == s) {
- derivatives.put(key, 0.);
- } else if (n == s + 1) {
- dkdxns = (1. + 2. * s) * FastMath.pow(X, 2 * s) / FastMath.pow(2, s);
- } else {
- final NSKey keymN = new NSKey(-n, s);
- double dkmN = 0.;
- if (derivatives.containsKey(keymN)) {
- dkmN = derivatives.get(keymN);
- } else {
- dkmN = computeDerivative(n - 1, s);
- }
- final NSKey keymNp1 = new NSKey(-n + 1, s);
- double dkmNp1 = 0.;
- if (derivatives.containsKey(keymNp1)) {
- dkmNp1 = derivatives.get(keymNp1);
- } else {
- dkmNp1 = computeDerivative(n - 2, s);
- }
- final double kns = getValue(-n - 1, s);
- dkdxns = (n - 1) * XX * ((2. * n - 3.) * dkmN - (n - 2.) * dkmNp1) / ((n + s - 1.) * (n - s - 1.));
- dkdxns += 2. * kns / X;
- }
-
- derivatives.put(key, dkdxns);
- return dkdxns;
- }
-
- }
-
}
diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenTesseralLinear.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenTesseralLinear.java
new file mode 100644
index 0000000000..cc07df5c61
--- /dev/null
+++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenTesseralLinear.java
@@ -0,0 +1,451 @@
+/* Copyright 2002-2014 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.utilities.hansen;
+
+import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
+import org.apache.commons.math3.util.FastMath;
+import org.orekit.propagation.semianalytical.dsst.utilities.NewcombOperators;
+
+/**
+ * Hansen coefficients K(t,n,s) for t!=0 and n < 0.
+ * <p>
+ * Implements Collins 4-236 or Danielson 2.7.3-(9) for Hansen Coefficients and
+ * Collins 4-240 for derivatives. The recursions are transformed into
+ * composition of linear transformations to obtain the associated polynomials
+ * for coefficients and their derivatives - see Petre's paper
+ *
+ * @author Petre Bazavan
+ * @author Lucian Barbulescu
+ */
+public class HansenTesseralLinear {
+ /**
+ * The first vector of polynomials associated to Hansen coefficients and
+ * derivatives.
+ */
+ private PolynomialFunction[][] mpvec;
+
+ /** The second vector of polynomials associated only to derivatives. */
+ private PolynomialFunction[][] mpvecDeriv;
+
+ /** The minimum value for the order. */
+ private int Nmin;
+
+ /** The index of the initial condition, Petre's paper. */
+ private int N0;
+
+ /** The s coefficient. */
+ private int s;
+
+ /** The j coefficient. */
+ private int j;
+
+ /**
+ * The offset used to identify the polynomial that corresponds to a negative.
+ * value of n in the internal array that starts at 0
+ */
+ private int offset;
+
+ /** The objects used to calculate initial data by means of newcomb operators. */
+ private HansenCoefficientsBySeries[] hansenInit;
+
+ /** Initial values of Hansen coefficients. */
+ private double[] hansenValue;
+ /** Initial values of Hansen coefficients derivatives. */
+ private double[] hansenDeriv;
+
+ /**
+ * Constructor.
+ *
+ * @param nMax the maximum (absolute) value of n parameter
+ * @param s s parameter
+ * @param j j parameter
+ * @param n0 the minimum (absolute) value of n
+ */
+ public HansenTesseralLinear(final int nMax, final int s, final int j, final int n0) {
+ //Initialize the fields
+ this.offset = nMax + 1;
+ this.Nmin = -nMax - 1;
+ this.N0 = -n0 - 4;
+ this.s = s;
+ this.j = j;
+
+ // Create the objects used to generate the initial values
+ this.hansenInit = new HansenCoefficientsBySeries[4];
+ for (int i = 0; i < 4; i++) {
+ this.hansenInit[i] = new HansenCoefficientsBySeries(N0 + i, s, j);
+ }
+
+ this.hansenValue = new double[4];
+ this.hansenDeriv = new double[4];
+
+ // The first 4 values are computed with series. No linear combination is needed.
+ final int size = N0 - Nmin;
+ if (size > 0) {
+ mpvec = new PolynomialFunction[size][];
+ mpvecDeriv = new PolynomialFunction[size][];
+
+ // Prepare the database of the associated polynomials
+ generatePolynomials();
+ }
+ }
+
+ /**
+ * Compute polynomial coefficient a.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>j</sub><sup>-n, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
+ * and the coefficient for dK<sub>j</sub><sup>-n, s</sup> / de<sup>2</sup> when computing dK<sub>j</sub><sup>-n-1, s</sup> / de<sup>2</sup>
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
+ * </p>
+ *
+ * @param mnm1 -n-1
+ * @return the polynomial
+ */
+ private PolynomialFunction a(final int mnm1) {
+ // Collins 4-236, Danielson 2.7.3-(9)
+ final double r1 = (mnm1 + 2.) * (2. * mnm1 + 5.);
+ final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
+ return new PolynomialFunction(new double[] {
+ 0.0, 0.0, r1 / r2
+ });
+ }
+
+ /**
+ * Compute polynomial coefficient b.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>j</sub><sup>-n+1, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
+ * and the coefficient for dK<sub>j</sub><sup>-n+1, s</sup> / de<sup>2</sup> when computing dK<sub>j</sub><sup>-n-1, s</sup> / de<sup>2</sup>
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
+ * </p>
+ *
+ * @param mnm1 -n-1
+ * @return the polynomial
+ */
+ private PolynomialFunction b(final int mnm1) {
+ // Collins 4-236, Danielson 2.7.3-(9)
+ final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
+ final double d1 = (mnm1 + 3.) * 2. * j * s / (r2 * (mnm1 + 4.));
+ final double d2 = (mnm1 + 3.) * (mnm1 + 2.) / r2;
+ return new PolynomialFunction(new double[] {
+ 0.0, -d1, -d2
+ });
+ }
+
+ /**
+ * Compute polynomial coefficient c.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>j</sub><sup>-n+3, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
+ * and the coefficient for dK<sub>j</sub><sup>-n+3, s</sup> / de<sup>2</sup> when computing dK<sub>j</sub><sup>-n-1, s</sup> / de<sup>2</sup>
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
+ * </p>
+ *
+ * @param mnm1 -n-1
+ * @return the polynomial
+ */
+ private PolynomialFunction c(final int mnm1) {
+ // Collins 4-236, Danielson 2.7.3-(9)
+ final double r1 = j * j * (mnm1 + 2.);
+ final double r2 = (mnm1 + 4.) * (2. + mnm1 + s) * (2. + mnm1 - s);
+
+ return new PolynomialFunction(new double[] {
+ 0.0, 0.0, r1 / r2
+ });
+ }
+
+ /**
+ * Compute polynomial coefficient d.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>j</sub><sup>-n-1, s</sup> / dχ when computing dK<sub>j</sub><sup>-n-1, s</sup> / de<sup>2</sup>
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
+ * </p>
+ *
+ * @param mnm1 -n-1
+ * @return the polynomial
+ */
+ private PolynomialFunction d(final int mnm1) {
+ // Collins 4-236, Danielson 2.7.3-(9)
+ return new PolynomialFunction(new double[] {
+ 0.0, 0.0, 1.0
+ });
+ }
+
+ /**
+ * Compute polynomial coefficient f.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>j</sub><sup>-n+1, s</sup> / dχ when computing dK<sub>j</sub><sup>-n-1, s</sup> / de<sup>2</sup>
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
+ * </p>
+ *
+ * @param n index
+ * @return the polynomial
+ */
+ private PolynomialFunction f(final int n) {
+ // Collins 4-236, Danielson 2.7.3-(9)
+ final double r1 = (n + 3.0) * j * s;
+ final double r2 = (n + 4.0) * (2.0 + n + s) * (2.0 + n - s);
+ return new PolynomialFunction(new double[] {
+ 0.0, 0.0, 0.0, r1 / r2
+ });
+ }
+
+ /**
+ * Generate the polynomials needed in the linear transformation.
+ *
+ * <p>
+ * See Petre's paper
+ * </p>
+ */
+ private void generatePolynomials() {
+
+
+ // The matrix that contains the coefficients at each step
+ final PolynomialFunctionMatrix a = new PolynomialFunctionMatrix(4);
+
+ // Initialization of the matrices for linear transformations
+ // The final configuration of these matrices are obtained by composition
+ // of linear transformations
+
+ // The matrix of polynomials associated to Hansen coefficients, Petre's
+ // paper
+ PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix4();
+
+ // The matrix of polynomials associated to derivatives, Petre's paper
+ final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix4();
+ PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix4();
+
+ // The matrix of the current linear transformation
+ a.setMatrixLine(0, new PolynomialFunction[] {
+ HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO, HansenUtilities.ZERO
+ });
+ a.setMatrixLine(1, new PolynomialFunction[] {
+ HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO
+ });
+ a.setMatrixLine(2, new PolynomialFunction[] {
+ HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE
+ });
+ // The generation process
+ int index;
+ for (int i = N0 - 1; i > Nmin - 1; i--) {
+ index = i + this.offset;
+ // The matrix of the current linear transformation is updated
+ // Petre's paper
+ a.setMatrixLine(3, new PolynomialFunction[] {
+ c(i), HansenUtilities.ZERO, b(i), a(i)
+ });
+
+ // composition of the linear transformations to calculate
+ // the polynomials associated to Hansen coefficients
+ // Petre's paper
+ A = A.multiply(a);
+ // store the polynomials for Hansen coefficients
+ mpvec[index] = A.getMatrixLine(3);
+ // composition of the linear transformations to calculate
+ // the polynomials associated to derivatives
+ // Petre's paper
+ D = D.multiply(a);
+
+ //Update the B matrix
+ B.setMatrixLine(3, new PolynomialFunction[] {
+ HansenUtilities.ZERO, f(i),
+ HansenUtilities.ZERO, d(i)
+ });
+ D = D.add(A.multiply(B));
+
+ // store the polynomials for Hansen coefficients from the
+ // expressions of derivatives
+ mpvecDeriv[index] = D.getMatrixLine(3);
+ }
+ }
+
+ /**
+ * Compute the values for the first four coefficients and their derivatives by means of series.
+ *
+ * @param e2 e<sup>2</sup>
+ * @param chi Χ
+ * @param chi2 Χ<sup>2</sup>
+ * @param precision the precision that will be used in the series
+ */
+ public void computeInitValues(final double e2, final double chi, final double chi2, final int precision) {
+ for (int i = 0; i < 4; i++) {
+ this.hansenValue[i] = this.hansenInit[i].getValue(e2, chi, chi2, precision);
+ this.hansenDeriv[i] = this.hansenInit[i].getDerivativeValue(e2, chi, chi2,
+ precision, this.hansenValue[i]);
+ }
+ }
+
+ /**
+ * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
+ *
+ * @param mnm1 -n-1
+ * @param chi χ
+ * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
+ */
+ public double getValue(final int mnm1, final double chi) {
+
+ // Check if the required coefficient is one of the initialization values
+ final int i = mnm1 - N0;
+ if (i >= 0 && i < 4) {
+ return hansenValue[i];
+ }
+
+ // Computes the coefficient by linear transformation
+ // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
+ final PolynomialFunction[] vv = mpvec[mnm1 + offset];
+ final double v0 = vv[0].value(chi);
+ final double v1 = vv[1].value(chi);
+ final double v2 = vv[2].value(chi);
+ final double v3 = vv[3].value(chi);
+ return v0 * hansenValue[3] + v1 * hansenValue[2] + v2 * hansenValue[1] + v3 * hansenValue[0];
+
+ }
+
+ /**
+ * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de<sup>2</sup>.
+ *
+ * @param mnm1 -n-1
+ * @param chi χ
+ * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de<sup>2</sup>
+ */
+ public double getDerivative(final int mnm1, final double chi) {
+
+ // Check if the required derivative is one of the initialization values
+ final int i = mnm1 - N0;
+ if (i >= 0 && i < 4) {
+ return hansenDeriv[i];
+ }
+
+ // Compute the derivative by linear transformation
+ // Collins 4-240 and Petre's paper
+ PolynomialFunction[] vv = mpvec[mnm1 + offset];
+ final double dv0 = vv[0].value(chi);
+ final double dv1 = vv[1].value(chi);
+ final double dv2 = vv[2].value(chi);
+ final double dv3 = vv[3].value(chi);
+ double ret = dv0 * hansenDeriv[3] + dv1 * hansenDeriv[2] + dv2 * hansenDeriv[1] + dv3 * hansenDeriv[0];
+
+ vv = mpvecDeriv[mnm1 + offset];
+ final double v0 = vv[0].value(chi);
+ final double v1 = vv[1].value(chi);
+ final double v2 = vv[2].value(chi);
+ final double v3 = vv[3].value(chi);
+ ret += v0 * hansenValue[3] + v1 * hansenValue[2] + v2 * hansenValue[1] + v3 * hansenValue[0];
+
+ return ret;
+
+ }
+
+ /**
+ * Compute a Hansen coefficient with series.
+ *
+ */
+ private class HansenCoefficientsBySeries {
+
+ /** -n-1 coefficient. */
+ private final int mnm1;
+ /** s coefficient. */
+ private final int s;
+ /** j coefficient. */
+ private final int j;
+
+ /**
+ * Class constructor.
+ *
+ * @param mnm1 -n-1 value
+ * @param s s value
+ * @param j j value
+ */
+ public HansenCoefficientsBySeries(final int mnm1, final int s, final int j) {
+ this.mnm1 = mnm1;
+ this.s = s;
+ this.j = j;
+ }
+
+ /**
+ * Compute the value of K<sub>j</sub><sup>-n-1, s</sup> with series.
+ *
+ * @param e2 e<sup>2</sup>
+ * @param chi Χ
+ * @param chi2 Χ<sup>2</sup>
+ * @param maxNewcomb Max power of e<sup>2</sup> in series expansion
+ * @return the value of K<sub>j</sub><sup>-n-1, s</sup>
+ */
+ public double getValue(final double e2, final double chi, final double chi2, final int maxNewcomb) {
+ // Initialization
+ final int a = FastMath.max(j - s, 0);
+ final int b = FastMath.max(s - j, 0);
+ // Expansion until maxNewcomb, the maximum power in e^2 for the
+ // Kernel value
+ double sum = NewcombOperators.getValue(maxNewcomb + a, maxNewcomb + b, mnm1, s);
+ for (int alpha = maxNewcomb - 1; alpha >= 0; alpha--) {
+ sum *= e2;
+ sum += NewcombOperators.getValue(alpha + a, alpha + b, mnm1, s);
+ }
+ // Kernel value from equation 2.7.3-(10)
+ final double value = FastMath.pow(chi2, -mnm1 - 1) * sum / chi;
+
+ return value;
+ }
+
+ /**
+ * Compute the value of dK<sub>j</sub><sup>-n-1, s</sup> / de<sup>2</sup> with series.
+ *
+ * @param e2 e<sup>2</sup>
+ * @param chi Χ
+ * @param chi2 Χ<sup>2</sup>
+ * @param maxNewcomb Max power of e<sup>2</sup> in series expansion
+ * @param Kjmnm1s the value of K<sub>j</sub><sup>-n-1, s</sup>
+ * @return derivative
+ */
+ public double getDerivativeValue(final double e2, final double chi, final double chi2,
+ final int maxNewcomb, final double Kjmnm1s) {
+ // Initialization
+ final int a = FastMath.max(j - s, 0);
+ final int b = FastMath.max(s - j, 0);
+ // Expansion until maxNewcomb-1, the maximum power in e^2 for the
+ // Kernel derivative
+ double sum = maxNewcomb * NewcombOperators.getValue(maxNewcomb + a, maxNewcomb + b, mnm1, s);
+ for (int alpha = maxNewcomb - 1; alpha >= 1; alpha--) {
+ sum *= e2;
+ sum += alpha * NewcombOperators.getValue(alpha + a, alpha + b, mnm1, s);
+ }
+ // Kernel derivative from equation 3.3-(5)
+ final double coef = -(mnm1 + 1.5);
+ final double derivative = coef * chi2 * Kjmnm1s + FastMath.pow(chi2, -mnm1 - 1) * sum / chi;
+
+ return derivative;
+ }
+ }
+}
diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenThirdBodyLinear.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenThirdBodyLinear.java
new file mode 100644
index 0000000000..9f0459d63f
--- /dev/null
+++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenThirdBodyLinear.java
@@ -0,0 +1,322 @@
+/* Copyright 2002-2014 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.utilities.hansen;
+
+import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
+
+/**
+ * Hansen coefficients K(t,n,s) for t=0 and n > 0.
+ * <p>
+ * Implements Collins 4-254 or Danielson 2.7.3-(7) for Hansen Coefficients and
+ * Danielson 3.2-(3) for derivatives. The recursions are transformed into
+ * composition of linear transformations to obtain the associated polynomials
+ * for coefficients and their derivatives - see Petre's paper
+ *
+ * @author Petre Bazavan
+ * @author Lucian Barbulescu
+ */
+public class HansenThirdBodyLinear {
+ /**
+ * The first vector of polynomials associated to Hansen coefficients and
+ * derivatives.
+ */
+ private PolynomialFunction[][] mpvec;
+
+ /** The second vector of polynomials associated only to derivatives. */
+ private PolynomialFunction[][] mpvecDeriv;
+
+ /** The maximum order of n indexes. */
+ private int Nmax;
+
+ /** The index of the initial condition, Petre's paper. */
+ private int N0;
+
+ /** The s index. */
+ private int s;
+
+ /** The value of K<sub>0</sub><sup>s,s</sup> computed with Collins (4-255). */
+ private double hansenS;
+
+ /** The value of K<sub>0</sub><sup>s+1,s</sup> computed with Collins (4-256). */
+ private double hansenSp1;
+ /** Coefficient used to compute K<sub>0</sub><sup>s+1,s</sup>.
+ *
+ * <p>
+ * hansenSp1Coef1 = ( (2*s+1)!! / (s+2)! )
+ * </p>
+ * see Collins (4-256)
+ */
+ private double hansenSp1Coef1;
+ /** Coefficient used to compute K<sub>0</sub><sup>s+1,s</sup>.
+ *
+ * <p>
+ * hansenSp1Coef2 = ( 2*s+3 )
+ * </p>
+ * see Collins (4-256)
+ */
+ private double hansenSp1Coef2;
+
+ /** The value of dK<sub>0</sub><sup>s+1,s</sup> / dΧ computed with Collins (4-259). */
+ private double hansenDerivSp1;
+ /** Coefficient used to compute dK<sub>0</sub><sup>s+1,s</sup> / dΧ.
+ *
+ * <p>
+ * hansenDerivSp1Coef1 = ( 2 * (2*s+1)!! / (s+2)! )
+ * </p>
+ * see Collins (4-256)
+ */
+ private double hansenDerivSp1Coef1;
+
+ /**
+ * Constructor.
+ *
+ * @param Nmax the maximum value of n
+ * @param s the value of s
+ */
+ public HansenThirdBodyLinear(final int Nmax, final int s) {
+
+ this.Nmax = Nmax;
+ N0 = s;
+ this.s = s;
+
+ // initialization of structures for stored data
+ mpvec = new PolynomialFunction[this.Nmax + 1][];
+ mpvecDeriv = new PolynomialFunction[this.Nmax + 1][];
+
+ // Prepare the database of the associated polynomials
+ generatePolynomials();
+ computeInitialValues();
+ }
+
+ /**
+ * Compute polynomial coefficient a.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>0</sub><sup>n-1, s</sup> when computing K<sub>0</sub><sup>n, s</sup>
+ * and the coefficient for dK<sub>0</sub><sup>n-1, s</sup> / dΧ when computing dK<sub>0</sub><sup>n, s</sup> / dΧ
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
+ * </p>
+ *
+ * @param n n value
+ * @return the polynomial
+ */
+ private PolynomialFunction a(final int n) {
+ // from recurrence Danielson 2.7.3-(7c), Collins 4-254
+ final double r1 = 2 * n + 1;
+ final double r2 = n + 1;
+
+ return new PolynomialFunction(new double[] {
+ r1 / r2
+ });
+ }
+
+ /**
+ * Compute polynomial coefficient b.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>0</sub><sup>n-2, s</sup> when computing K<sub>0</sub><sup>n, s</sup>
+ * and the coefficient for dK<sub>0</sub><sup>n-2, s</sup> / dΧ when computing dK<sub>0</sub><sup>n, s</sup> / dΧ
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
+ * </p>
+ *
+ * @param n n value
+ * @return the polynomial
+ */
+ private PolynomialFunction b(final int n) {
+ // from recurrence Danielson 2.7.3-(7c), Collins 4-254
+ final double r1 = (n + s) * (n - s);
+ final double r2 = n * (n + 1);
+
+ return new PolynomialFunction(new double[] {
+ 0.0, 0.0, -r1 / r2
+ });
+ }
+
+ /**
+ * Compute polynomial coefficient d.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>0</sub><sup>n-2, s</sup> when computing dK<sub>0</sub><sup>n, s</sup> / dΧ
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
+ * </p>
+ *
+ * @param n n value
+ * @return the polynomial
+ */
+ private PolynomialFunction d(final int n) {
+ // from Danielson 3.2-(3b)
+ final double r1 = 2 * (n + s) * (n - s);
+ final double r2 = n * (n + 1);
+
+ return new PolynomialFunction(new double[] {
+ 0.0, 0.0, 0.0, r1 / r2
+ });
+ }
+
+ /**
+ * Generate the polynomials needed in the linear transformation.
+ *
+ * <p>
+ * See Petre's paper
+ * </p>
+ */
+ private void generatePolynomials() {
+
+ // Initialization of the matrices for linear transformations
+ // The final configuration of these matrices are obtained by composition
+ // of linear transformations
+
+ // the matrix A for the polynomials associated
+ // to Hansen coefficients, Petre's pater
+ PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix2();
+
+ // the matrix D for the polynomials associated
+ // to derivatives, Petre's paper
+ final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix2();
+ PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix2();
+ PolynomialFunctionMatrix E = HansenUtilities.buildIdentityMatrix2();
+
+ // The matrix that contains the coefficients at each step
+ final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix2();
+ a.setElem(0, 1, HansenUtilities.ONE);
+ // The generation process
+ for (int i = N0 + 2; i <= Nmax; i++) {
+ // Collins 4-254 or Danielson 2.7.3-(7)
+ // Petre's paper
+ // The matrix of the current linear transformation is actualized
+ a.setMatrixLine(1, new PolynomialFunction[] {
+ b(i), a(i)
+ });
+
+ // composition of the linear transformations to calculate
+ // the polynomials associated to Hansen coefficients
+ A = A.multiply(a);
+ // store the polynomials associated to Hansen coefficients
+ this.mpvec[i] = A.getMatrixLine(1);
+ // composition of the linear transformations to calculate
+ // the polynomials associated to derivatives
+ // Danielson 3.2-(3b) and Petre's paper
+ D = D.multiply(a);
+ if (i > N0 + 2) {
+ a.setMatrixLine(1, new PolynomialFunction[] {
+ b(i - 1), a(i - 1)
+ });
+ E = E.multiply(a);
+ }
+
+ B.setElem(1, 0, d(i));
+ // F = E.prod(B);
+ D = D.add(E.multiply(B));
+ // store the polynomials associated to the derivatives
+ this.mpvecDeriv[i] = D.getMatrixLine(1);
+ }
+ }
+
+ /**
+ * Compute the initial values (see Collins, 4-255, 4-256 and 4.259)
+ * <p>
+ * K<sub>0</sub><sup>s, s</sup> = (-1)<sup>s</sup> * ( (2*s+1)!! / (s+1)! )
+ * </p>
+ * <p>
+ * K<sub>0</sub><sup>s+1, s</sup> = (-1)<sup>s</sup> * ( (2*s+1)!! / (s+2)!
+ * ) * (2*s+3 - χ<sup>-2</sup>)
+ * </p>
+ * <p>
+ * dK<sub>0</sub><sup>s+1, s</sup> / dχ = = (-1)<sup>s</sup> * 2 * (
+ * (2*s+1)!! / (s+2)! ) * χ<sup>-3</sup>
+ * </p>
+ */
+ private void computeInitialValues() {
+ // the value for K<sub>0</sub><sup>s, s</s> does not depend of χ
+ hansenS = (s % 2 == 0) ? 1.0 : -1.0;
+ for (int i = s; i >= 1; i--) {
+ hansenS *= (2.0 * i + 1.0) / (i + 1.0);
+ }
+
+ // the value for K<sub>0</sub><sup>s+1, s</s> is hansenSplusoneCoef1 *
+ // (hansenSplusoneCoef2 - χ<sup>-2</sup>)
+ hansenSp1Coef1 = hansenS / (s + 2.);
+ hansenSp1Coef2 = 2. * s + 3.;
+
+ // the value for dK<sub>0</sub><sup>s+1, s</s> is hansenDerivSplusone *
+ // χ<sup>-3</sup>
+ hansenDerivSp1Coef1 = hansenSp1Coef1 * 2.;
+ }
+
+ /**
+ * Compute the value of the Hansen coefficient K<sub>0</sub><sup>n, s</sup>.
+ *
+ * @param n n value
+ * @param chitm1 χ<sup>-1</sup>
+ * @return the coefficient K<sub>0</sub><sup>n, s</sup>
+ */
+ public double getValue(final int n, final double chitm1) {
+ // Danielson 2.7.3-(7a,b)
+ if (n == s) {
+ return hansenS;
+ }
+
+ //Compute K<sub>0</sub><sup>s+1, s</sup>
+ hansenSp1 = hansenSp1Coef1 * (hansenSp1Coef2 - chitm1 * chitm1);
+ if (n == s + 1) {
+ return hansenSp1;
+ }
+ // Computes the coefficient by linear transformation
+ // Danielson 2.7.3-(7c)/Collins 4-254 and Petre's paper
+ final PolynomialFunction[] vv = mpvec[n];
+ return vv[0].value(chitm1) * hansenS + vv[1].value(chitm1) * hansenSp1;
+
+ }
+
+ /**
+ * Compute the value of the Hansen coefficient dK<sub>0</sub><sup>n, s</sup> / dΧ.
+ *
+ * @param n n value
+ * @param chitm1 χ<sup>-1</sup>
+ * @return the coefficient dK<sub>0</sub><sup>n, s</sup> / dΧ
+ */
+ public double getDerivative(final int n, final double chitm1) {
+ if (n == s) {
+ return 0;
+ }
+ //Compute dK<sub>0</sub><sup>s+1, s</sup> / dΧ
+ hansenDerivSp1 = hansenDerivSp1Coef1 * chitm1 * chitm1 * chitm1;
+ if (n == s + 1) {
+ return hansenDerivSp1;
+ }
+ //Compute K<sub>0</sub><sup>s+1, s</sup>
+ hansenSp1 = hansenSp1Coef1 * (hansenSp1Coef2 - chitm1 * chitm1);
+
+ // Computes the coefficient by linear transformation
+ // Danielson 2.7.3-(7c)/Collins 4-254 and Petre's paper
+ final PolynomialFunction[] v = mpvec[n];
+ final PolynomialFunction[] vv = mpvecDeriv[n];
+ return v[1].value(chitm1) * hansenDerivSp1 +
+ vv[0].value(chitm1) * hansenS + vv[1].value(chitm1) * hansenSp1;
+
+ }
+
+}
diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenUtilities.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenUtilities.java
new file mode 100644
index 0000000000..bfd466e7b6
--- /dev/null
+++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenUtilities.java
@@ -0,0 +1,160 @@
+/* Copyright 2002-2014 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.utilities.hansen;
+
+import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
+
+/**
+ * Utilities class.
+ *
+ * @author Lucian Barbulescu
+ */
+public class HansenUtilities {
+
+ /** 1 represented as a polynomial. */
+ public static final PolynomialFunction ONE = new PolynomialFunction(new double[] {
+ 1
+ });
+
+ /** 0 represented as a polynomial. */
+ public static final PolynomialFunction ZERO = new PolynomialFunction(new double[] {
+ 0
+ });
+
+ /** Private constructor as class is a utility.
+ */
+ private HansenUtilities() {
+ }
+
+ /**
+ * Build the identity matrix of order 2.
+ * <p>
+ * <pre>
+ * / 1 0 \
+ * I<sub>2</sub> = | |
+ * \ 0 1 /
+ * </pre>
+ *
+ * @return the identity matrix of order 2
+ */
+ public static final PolynomialFunctionMatrix buildIdentityMatrix2() {
+ final PolynomialFunctionMatrix matrix = new PolynomialFunctionMatrix(2);
+ matrix.setMatrix(new PolynomialFunction[][] {
+ {
+ ONE, ZERO
+ },
+ {
+ ZERO, ONE
+ }
+ });
+ return matrix;
+ }
+
+ /**
+ * Build the empty matrix of order 2.
+ * <p>
+ * <pre>
+ * / 0 0 \
+ * E<sub>2</sub> = | |
+ * \ 0 0 /
+ * </pre>
+ *
+ * @return the identity matrix of order 2
+ */
+ public static final PolynomialFunctionMatrix buildZeroMatrix2() {
+ final PolynomialFunctionMatrix matrix = new PolynomialFunctionMatrix(2);
+ matrix.setMatrix(new PolynomialFunction[][] {
+ {
+ ZERO, ZERO
+ },
+ {
+ ZERO, ZERO
+ }
+ });
+ return matrix;
+ }
+
+
+ /**
+ * Build the identity matrix of order 4.
+ * <p>
+ * <pre>
+ * / 1 0 0 0 \
+ * | |
+ * | 0 1 0 0 |
+ * I<sub>4</sub> = | |
+ * | 0 0 1 0 |
+ * | |
+ * \ 0 0 0 1 /
+ * </pre>
+ *
+ * @return the identity matrix of order 4
+ */
+ public static final PolynomialFunctionMatrix buildIdentityMatrix4() {
+ final PolynomialFunctionMatrix matrix = new PolynomialFunctionMatrix(4);
+ matrix.setMatrix(new PolynomialFunction[][] {
+ {
+ ONE, ZERO, ZERO, ZERO
+ },
+ {
+ ZERO, ONE, ZERO, ZERO
+ },
+ {
+ ZERO, ZERO, ONE, ZERO
+ },
+ {
+ ZERO, ZERO, ZERO, ONE
+ }
+ });
+ return matrix;
+ }
+
+ /**
+ * Build the empty matrix of order 4.
+ * <p>
+ * <pre>
+ * / 0 0 0 0 \
+ * | |
+ * | 0 0 0 0 |
+ * E<sub>4</sub> = | |
+ * | 0 0 0 0 |
+ * | |
+ * \ 0 0 0 0 /
+ * </pre>
+ *
+ * @return the identity matrix of order 4
+ */
+ public static final PolynomialFunctionMatrix buildZeroMatrix4() {
+ final PolynomialFunctionMatrix matrix = new PolynomialFunctionMatrix(4);
+ matrix.setMatrix(new PolynomialFunction[][] {
+ {
+ ZERO, ZERO, ZERO, ZERO
+ },
+ {
+ ZERO, ZERO, ZERO, ZERO
+ },
+ {
+ ZERO, ZERO, ZERO, ZERO
+ },
+ {
+ ZERO, ZERO, ZERO, ZERO
+ }
+ } );
+ return matrix;
+ }
+
+}
diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenZonalLinear.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenZonalLinear.java
new file mode 100644
index 0000000000..446e714819
--- /dev/null
+++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/HansenZonalLinear.java
@@ -0,0 +1,250 @@
+/* Copyright 2002-2014 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.utilities.hansen;
+
+import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * Hansen coefficients K(t,n,s) for t=0 and n < 0.
+ * <p>
+ *Implements Collins 4-242 or echivalently, Danielson 2.7.3-(6) for Hansen Coefficients and
+ * Collins 4-245 or Danielson 3.1-(7) for derivatives. The recursions are transformed into
+ * composition of linear transformations to obtain the associated polynomials
+ * for coefficients and their derivatives - see Petre's paper
+ *
+ * @author Petre Bazavan
+ * @author Lucian Barbulescu
+ */
+public class HansenZonalLinear {
+ /**
+ * The first vector of polynomials associated to Hansen coefficients and
+ * derivatives.
+ */
+ private PolynomialFunction[] mpvec;
+
+ /** The second vector of polynomials associated only to derivatives. */
+ private PolynomialFunction[] mpvecDeriv;
+
+ /** The minimum value for the order. */
+ private int Nmin;
+
+
+ /** The index of the initial condition, Petre's paper. */
+ private int N0;
+
+ /** The s coefficient. */
+ private int s;
+
+ /**
+ * The offset used to identify the polynomial that corresponds to a negative
+ * value of n in the internal array that starts at 0.
+ */
+ private int offset;
+
+ /** 2<sup>s</sup>. */
+ private double twots;
+
+ /** 2<sup>s-1</sup>. */
+ private double twotsm1;
+
+ /** 2*s+1. */
+ private int twosp1;
+
+ /** 2*s-2. */
+ private int twosm2;
+
+ /** 2*s. */
+ private int twos;
+
+ /**
+ * Constructor.
+ *
+ * @param nMax the maximum (absolute) value of n coefficient
+ * @param s s coefficient
+ */
+ public HansenZonalLinear(final int nMax, final int s) {
+
+ //Initialize fields
+ this.offset = nMax + 1;
+ this.Nmin = -nMax - 1;
+ N0 = -(s + 2);
+ this.s = s;
+ this.twots = FastMath.pow(2., s);
+ this.twotsm1 = FastMath.pow(2., s - 1);
+ this.twos = 2 * s;
+ this.twosp1 = this.twos + 1;
+ this.twosm2 = this.twos - 2;
+
+ // prepare structures for stored data
+ final int size = nMax - s - 1;
+ mpvec = new PolynomialFunction[size];
+ mpvecDeriv = new PolynomialFunction[size];
+
+ // Prepare the data base of associated polynomials
+ generatePolynomials();
+ }
+
+ /**
+ * Compute polynomial coefficient a.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>0</sub><sup>-n, s</sup> when computing K<sub>0</sub><sup>-n-1, s</sup>
+ * and the coefficient for dK<sub>0</sub><sup>-n, s</sup> / de<sup>2</sup> when computing dK<sub>0</sub><sup>-n-1, s</sup> / de<sup>2</sup>
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(6) and Collins 4-242 and 4-245
+ * </p>
+ *
+ * @param mnm1 -n-1 value
+ * @return the polynomial
+ */
+ private PolynomialFunction a(final int mnm1) {
+ // from recurence Collins 4-242
+ final double d1 = (mnm1 + 2) * (2 * mnm1 + 5);
+ final double d2 = (mnm1 + 2 - s) * (mnm1 + 2 + s);
+ return new PolynomialFunction(new double[] {
+ 0.0, 0.0, d1 / d2
+ });
+ }
+
+ /**
+ * Compute polynomial coefficient b.
+ *
+ * <p>
+ * It is used to generate the coefficient for K<sub>0</sub><sup>-n+1, s</sup> when computing K<sub>0</sub><sup>-n-1, s</sup>
+ * and the coefficient for dK<sub>0</sub><sup>-n+1, s</sup> / de<sup>2</sup> when computing dK<sub>0</sub><sup>-n-1, s</sup> / de<sup>2</sup>
+ * </p>
+ *
+ * <p>
+ * See Danielson 2.7.3-(6) and Collins 4-242 and 4-245
+ * </p>
+ *
+ * @param mnm1 -n-1 value
+ * @return the polynomial
+ */
+ private PolynomialFunction b(final int mnm1) {
+ // from recurence Collins 4-242
+ final double d1 = (mnm1 + 2) * (mnm1 + 3);
+ final double d2 = (mnm1 + 2 - s) * (mnm1 + 2 + s);
+ return new PolynomialFunction(new double[] {
+ 0.0, 0.0, -d1 / d2
+ });
+ }
+
+ /**
+ * Generate the polynomials needed in the linear transformation.
+ *
+ * <p>
+ * See Petre's paper
+ * </p>
+ */
+ private void generatePolynomials() {
+
+ // Initialisation of matrix for linear transformmations
+ // The final configuration of these matrix are obtained by composition
+ // of linear transformations
+ PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix2();
+ final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix2();
+ PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix2();
+ PolynomialFunctionMatrix E = HansenUtilities.buildIdentityMatrix2();
+
+ // from Collins 4-245 and Petre's paper
+ B.setElem(1, 1, new PolynomialFunction(new double[] {
+ 0, 0, 1
+ }));
+
+ // generation of polynomials associated to Hansen coefficients and to
+ // their derivatives
+ final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix2();
+ a.setElem(0, 1, HansenUtilities.ONE);
+
+ int index;
+ for (int i = N0 - 1; i > Nmin - 1; i--) {
+ index = i + offset;
+ // Matrix of the current linear transformation
+ // Petre's paper
+ a.setMatrixLine(1, new PolynomialFunction[] {
+ b(i), a(i)
+ });
+ // composition of linear transformations
+ // see Petre's paper
+ A = A.multiply(a);
+ // store the polynomials for Hansen coefficients
+ mpvec[index] = A.getElem(1, 1);
+
+ D = D.multiply(a);
+ E = E.multiply(a);
+ D = D.add(E.multiply(B));
+
+ // store the polynomials for Hansen coefficients from the expressions
+ // of derivatives
+ mpvecDeriv[index] = D.getElem(1, 1);
+ }
+ }
+
+ /**
+ * Get the K<sub>0</sub><sup>-n-1,s</sup> coefficient value. <br />
+ * The s value is given in the class constructor <br />
+ *
+ * @param mnm1 (-n-1) coefficient
+ * @param chi The value of χ
+ * @return K<sub>0</sub><sup>-n-1,s</sup>
+ */
+ public double getValue(final int mnm1, final double chi) {
+ // Danielson 2.7.3-(6a,b)
+ if (mnm1 == N0 + 1) {
+ return 0;
+ }
+
+ final double han1 = FastMath.pow(chi, twosp1) / twots;
+ if (mnm1 == N0) {
+ return han1;
+ }
+
+ // Danielson 2.7.3-(6c)/Collins 4-242 and Petre's paper
+ return mpvec[mnm1 + offset].value(chi) * han1;
+ }
+
+ /**
+ * Get the dK<sub>0</sub><sup>-n-1,s</sup> / dΧ coefficient derivative. <br />
+ * The s value is given in the class constructor <br />
+ *
+ * @param mnm1 (-n-1) coefficient
+ * @param chi The value of χ
+ * @return dK<sub>0</sub><sup>-n-1,s</sup> / dΧ
+ */
+ public double getDerivative(final int mnm1, final double chi) {
+ // Danielson 3.1-(7a)
+ if (mnm1 == N0 + 1) {
+ return 0;
+ }
+
+ final double hanDeriv1 = this.twosp1 * FastMath.pow(chi, twos) / twots;
+ if (mnm1 == N0) {
+ return hanDeriv1;
+ }
+
+ // Danielson 3.1-(7c) and Petre's paper
+ double ret = mpvec[mnm1 + offset].value(chi) * hanDeriv1;
+ // Danielson 2.7.3-(6b)
+ final double han1 = FastMath.pow(chi, twosm2) / twotsm1;
+ ret += mpvecDeriv[mnm1 + offset].value(chi) * han1;
+ return ret;
+ }
+}
diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/PolynomialFunctionMatrix.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/PolynomialFunctionMatrix.java
new file mode 100644
index 0000000000..369eac69e0
--- /dev/null
+++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/hansen/PolynomialFunctionMatrix.java
@@ -0,0 +1,146 @@
+/* Copyright 2002-2014 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.utilities.hansen;
+
+import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
+
+/**
+ * A quadratic matrix of
+ * {@link org.apache.commons.math3.analysis.polynomials.PolynomialFunction}.
+ *
+ * @author Petre Bazavan
+ * @author Lucian Barbulescu
+ */
+public class PolynomialFunctionMatrix {
+
+ /** The order of the matrix. */
+ private int order;
+ /** The elements of the matrix. */
+ private PolynomialFunction elements[][];
+
+ /**
+ * Create a matrix.
+ *
+ * <p>
+ * All elements are null
+ *
+ * @param order
+ * the order of the matrix
+ */
+ PolynomialFunctionMatrix(final int order) {
+ this.order = order;
+ this.elements = new PolynomialFunction[order][order];
+ }
+
+ /**
+ * Set an element of the matrix.
+ *
+ * @param line
+ * the line
+ * @param column
+ * the column
+ * @param value
+ * the value
+ */
+ public void setElem(final int line, final int column, final PolynomialFunction value) {
+ elements[line][column] = value;
+ }
+
+ /**
+ * Get the value of an element.
+ *
+ * @param line
+ * the line
+ * @param column
+ * the column
+ * @return the value
+ */
+ public PolynomialFunction getElem(final int line, final int column) {
+ return elements[line][column];
+ }
+
+ /**
+ * Multiply the argument matrix with the current matrix.
+ *
+ * @param matrix
+ * the argument matrix
+ * @return the result of the multiplication
+ */
+ public PolynomialFunctionMatrix multiply(final PolynomialFunctionMatrix matrix) {
+ final PolynomialFunctionMatrix result = new PolynomialFunctionMatrix(order);
+ for (int i = 0; i < order; i++) {
+ for (int j = 0; j < order; j++) {
+ PolynomialFunction cc = HansenUtilities.ZERO;
+ for (int k = 0; k < order; k++) {
+ cc = cc.add(matrix.getElem(i, k).multiply(elements[k][j]));
+ }
+ result.setElem(i, j, cc);
+ }
+ }
+ return result;
+ }
+
+ /**
+ * Set values for all elements.
+ *
+ * @param polynomials
+ * the values that will be used for the matrix
+ */
+ public void setMatrix(final PolynomialFunction[][] polynomials) {
+ elements = polynomials.clone();
+ }
+
+ /**
+ * Set the value of a line of the matrix.
+ *
+ * @param line
+ * the line number
+ * @param polynomials
+ * the values to set
+ */
+ public void setMatrixLine(final int line, final PolynomialFunction[] polynomials) {
+ elements[line] = polynomials;
+ }
+
+ /**
+ * Get a line of the matrix as a {@link PolynomialFunctionVector}.
+ *
+ * @param line
+ * the line number
+ * @return the line of the matrix as a vector
+ */
+ public PolynomialFunction[] getMatrixLine(final int line) {
+ return elements[line].clone();
+ }
+
+ /**
+ * Add the argument matrix with the current matrix.
+ *
+ * @param matrix
+ * the argument matrix
+ * @return the result of the addition
+ */
+ public PolynomialFunctionMatrix add(final PolynomialFunctionMatrix matrix) {
+ final PolynomialFunctionMatrix c = new PolynomialFunctionMatrix(order);
+ for (int i = 0; i < order; i++) {
+ for (int j = 0; j < order; j++) {
+ c.setElem(i, j, elements[i][j].add(matrix.getElem(i, j)));
+ }
+ }
+ return c;
+ }
+}
--
GitLab