Commit 29045223 authored by Bryan Cazabonne's avatar Bryan Cazabonne

Merge branch 'issue-736' into 'develop'

Fixed initialization of maxEccPow, maxHansen, and resOrders terms.

Closes #736

See merge request orekit/orekit!120
parents 7c1bd047 951a150d
Pipeline #868 passed with stages
in 28 minutes and 29 seconds
......@@ -21,6 +21,9 @@
</properties>
<body>
<release version="11.0" date="TBD" description="TBD">
<action dev="bryan" type="fix" issue="736">
Fixed NullPointerException in DSSTTesseral Hansen object.
</action>
<action dev="bryan" type="update" issue="601">
Changed getPVInPZ90() method to private.
</action>
......
......@@ -169,12 +169,21 @@ public class DSSTTesseral implements DSSTForceModel {
/** Maximum value for j. */
private final int maxFrequencyShortPeriodics;
/** Maximum power of the eccentricity to use in summation over s. */
private int maxEccPow;
/** Maximum power of the eccentricity to use in Hansen coefficient Kernel expansion. */
private int maxHansen;
/** Maximum value between maxOrderMdailyTesseralSP and maxOrderTesseralSP. */
private int mMax;
/** List of non resonant orders with j != 0. */
private final SortedMap<Integer, List<Integer> > nonResOrders;
/** List of resonant orders. */
private final List<Integer> resOrders;
/** Short period terms. */
private TesseralShortPeriodicCoefficients shortPeriodTerms;
......@@ -293,12 +302,16 @@ public class DSSTTesseral implements DSSTForceModel {
this.maxFrequencyShortPeriodics = maxFrequencyShortPeriodics;
// Initialize default values
this.nonResOrders = new TreeMap<Integer, List <Integer> >();
this.resOrders = new ArrayList<Integer>();
this.nonResOrders = new TreeMap<Integer, List <Integer>>();
pendingInitialization = true;
fieldShortPeriodTerms = new HashMap<>();
fieldHansen = new HashMap<>();
// Initialize default values
this.fieldShortPeriodTerms = new HashMap<>();
this.fieldHansen = new HashMap<>();
this.maxEccPow = 0;
this.maxHansen = 0;
}
......@@ -322,15 +335,21 @@ public class DSSTTesseral implements DSSTForceModel {
// Initializes specific parameters.
final DSSTTesseralContext context = initializeStep(auxiliaryElements, parameters);
// Set the highest power of the eccentricity in the analytical power
// series expansion for the averaged high order resonant central body
// spherical harmonic perturbation
maxEccPow = getMaxEccPow(auxiliaryElements.getEcc());
// Set the maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
maxHansen = maxEccPow / 2;
// The following terms are only used for hansen objects initialization
final double ratio = context.getRatio();
final int maxEccPow = context.getMaxEccPow();
final List<Integer> resOrders = context.getResOrders();
final double ratio = context.getRatio();
// Compute the non resonant tesseral harmonic terms if not set by the user
getNonResonantTerms(type, context);
getResonantAndNonResonantTerms(type, context.getOrbitPeriod(), ratio);
hansen = new HansenObjects(ratio, maxEccPow, resOrders, type);
hansen = new HansenObjects(ratio, type);
mMax = FastMath.max(maxOrderTesseralSP, maxOrderMdailyTesseralSP);
......@@ -360,18 +379,24 @@ public class DSSTTesseral implements DSSTForceModel {
// Initializes specific parameters.
final FieldDSSTTesseralContext<T> context = initializeStep(auxiliaryElements, parameters);
// Compute the non resonant tesseral harmonic terms if not set by the user
getNonResonantTerms(type, context, field);
// Set the highest power of the eccentricity in the analytical power
// series expansion for the averaged high order resonant central body
// spherical harmonic perturbation
maxEccPow = getMaxEccPow(auxiliaryElements.getEcc().getReal());
// Set the maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
maxHansen = maxEccPow / 2;
// The following terms are only used for hansen objects initialization
final T ratio = context.getRatio();
final int maxEccPow = context.getMaxEccPow();
final List<Integer> resOrders = context.getResOrders();
final T ratio = context.getRatio();
// Compute the non resonant tesseral harmonic terms if not set by the user
// Field information is not important here
getResonantAndNonResonantTerms(type, context.getOrbitPeriod().getReal(), ratio.getReal());
mMax = FastMath.max(maxOrderTesseralSP, maxOrderMdailyTesseralSP);
fieldHansen.put(field, new FieldHansenObjects<>(ratio, maxEccPow, resOrders, type, field));
fieldHansen.put(field, new FieldHansenObjects<>(ratio, type));
pendingInitialization = false;
}
......@@ -390,6 +415,30 @@ public class DSSTTesseral implements DSSTForceModel {
}
/**
* Get the maximum power of the eccentricity to use in summation over s.
* @param e eccentricity
* @return the maximum power of the eccentricity
*/
private int getMaxEccPow(final double e) {
// maxEccPow depends on satellite eccentricity
if (e <= 0.005) {
return 3;
} else if (e <= 0.02) {
return 4;
} else if (e <= 0.1) {
return 7;
} else if (e <= 0.2) {
return 10;
} else if (e <= 0.3) {
return 12;
} else if (e <= 0.4) {
return 15;
} else {
return 20;
}
}
/** Performs initialization at each integration step for the current force model.
* <p>
* This method aims at being called before mean elements rates computation.
......@@ -723,66 +772,29 @@ public class DSSTTesseral implements DSSTForceModel {
}
/**
* Get the non-resonant tesseral terms in the central body spherical harmonic field.
* Get the resonant and non-resonant tesseral terms in the central body spherical harmonic field.
*
* @param type type of the elements used during the propagation
* @param context container for attributes
* @param orbitPeriod Keplerian period
* @param ratio ratio of satellite period to central body rotation period
*/
private void getNonResonantTerms(final PropagationType type, final DSSTTesseralContext context) {
private void getResonantAndNonResonantTerms(final PropagationType type, final double orbitPeriod,
final double ratio) {
// Compute natural resonant terms
final double tolerance = 1. / FastMath.max(MIN_PERIOD_IN_SAT_REV,
MIN_PERIOD_IN_SECONDS / context.getOrbitPeriod());
MIN_PERIOD_IN_SECONDS / orbitPeriod);
// Search the resonant orders in the tesseral harmonic field
resOrders.clear();
nonResOrders.clear();
for (int m = 1; m <= maxOrder; m++) {
final double resonance = context.getRatio() * m;
final double resonance = ratio * m;
int jRes = 0;
final int jComputedRes = (int) FastMath.round(resonance);
if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance - jComputedRes) <= tolerance) {
// Store each resonant index and order
jRes = jComputedRes;
}
if (type == PropagationType.OSCULATING && maxDegreeTesseralSP >= 0 && m <= maxOrderTesseralSP) {
//compute non resonant orders in the tesseral harmonic field
final List<Integer> listJofM = new ArrayList<Integer>();
//for the moment we take only the pairs (j,m) with |j| <= maxDegree + maxEccPow (from |s-j| <= maxEccPow and |s| <= maxDegree)
for (int j = -maxFrequencyShortPeriodics; j <= maxFrequencyShortPeriodics; j++) {
if (j != 0 && j != jRes) {
listJofM.add(j);
}
}
nonResOrders.put(m, listJofM);
}
}
}
/**
* Get the non-resonant tesseral terms in the central body spherical harmonic field.
*
* @param <T> type of the elements
* @param type type of the elements used during the propagation
* @param context container for attributes
* @param field field used by default
*/
private <T extends RealFieldElement<T>> void getNonResonantTerms(final PropagationType type,
final FieldDSSTTesseralContext<T> context,
final Field<T> field) {
final T zero = field.getZero();
// Compute natural resonant terms
final T tolerance = FastMath.max(zero.add(MIN_PERIOD_IN_SAT_REV),
context.getOrbitPeriod().divide(MIN_PERIOD_IN_SECONDS).reciprocal()).reciprocal();
nonResOrders.clear();
for (int m = 1; m <= maxOrder; m++) {
final T resonance = context.getRatio().multiply(m);
int jRes = 0;
final int jComputedRes = (int) FastMath.round(resonance);
if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance.subtract(jComputedRes)).getReal() <= tolerance.getReal()) {
// Store each resonant index and order
resOrders.add(m);
jRes = jComputedRes;
}
......@@ -2320,7 +2332,7 @@ public class DSSTTesseral implements DSSTForceModel {
dUdGa = 0.;
// Compute only if there is at least one resonant tesseral
if (!context.getResOrders().isEmpty()) {
if (!resOrders.isEmpty()) {
// Gmsj and Hmsj polynomials
final GHmsjPolynomials ghMSJ = new GHmsjPolynomials(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), I);
......@@ -2335,7 +2347,7 @@ public class DSSTTesseral implements DSSTForceModel {
}
// SUM over resonant terms {j,m}
for (int m : context.getResOrders()) {
for (int m : resOrders) {
// Resonant index for the current resonant order
final int j = FastMath.max(1, (int) FastMath.round(context.getRatio() * m));
......@@ -2363,8 +2375,8 @@ public class DSSTTesseral implements DSSTForceModel {
double dUdGaSin = 0.;
// s-SUM from -sMin to sMax
final int sMin = FastMath.min(context.getMaxEccPow() - j, maxDegree);
final int sMax = FastMath.min(context.getMaxEccPow() + j, maxDegree);
final int sMin = FastMath.min(maxEccPow - j, maxDegree);
final int sMax = FastMath.min(maxEccPow + j, maxDegree);
for (int s = 0; s <= sMax; s++) {
//Compute the initial values for Hansen coefficients using newComb operators
......@@ -2545,7 +2557,7 @@ public class DSSTTesseral implements DSSTForceModel {
dUdGa = zero;
// Compute only if there is at least one resonant tesseral
if (!context.getResOrders().isEmpty()) {
if (!resOrders.isEmpty()) {
// Gmsj and Hmsj polynomials
final FieldGHmsjPolynomials<T> ghMSJ = new FieldGHmsjPolynomials<>(auxiliaryElements.getK(), auxiliaryElements.getH(), auxiliaryElements.getAlpha(), auxiliaryElements.getBeta(), I, field);
......@@ -2560,7 +2572,7 @@ public class DSSTTesseral implements DSSTForceModel {
}
// SUM over resonant terms {j,m}
for (int m : context.getResOrders()) {
for (int m : resOrders) {
// Resonant index for the current resonant order
final int j = FastMath.max(1, (int) FastMath.round(context.getRatio().multiply(m)));
......@@ -2588,8 +2600,8 @@ public class DSSTTesseral implements DSSTForceModel {
T dUdGaSin = zero;
// s-SUM from -sMin to sMax
final int sMin = FastMath.min(context.getMaxEccPow() - j, maxDegree);
final int sMax = FastMath.min(context.getMaxEccPow() + j, maxDegree);
final int sMin = FastMath.min(maxEccPow - j, maxDegree);
final int sMax = FastMath.min(maxEccPow + j, maxDegree);
for (int s = 0; s <= sMax; s++) {
//Compute the initial values for Hansen coefficients using newComb operators
......@@ -2713,27 +2725,17 @@ public class DSSTTesseral implements DSSTForceModel {
/** Computes init values of the Hansen Objects. */
private class HansenObjects {
/** Maximum power of the eccentricity to use in Hansen coefficient Kernel expansion. */
private int maxHansen;
/** 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;
/** Simple constructor.
* @param ratio Ratio of satellite period to central body rotation period
* @param maxEccPow Maximum power of the eccentricity to use in summation over s.
* @param resOrders List of resonant orders
* @param type type of the elements used during the propagation
*/
HansenObjects(final double ratio,
final int maxEccPow,
final List<Integer> resOrders,
final PropagationType type) {
// Set the maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
maxHansen = maxEccPow / 2;
//Allocate the two dimensional array
final int rows = 2 * maxDegree + 1;
final int columns = maxFrequencyShortPeriodics + 1;
......@@ -2804,26 +2806,17 @@ public class DSSTTesseral implements DSSTForceModel {
/** Computes init values of the Hansen Objects. */
private class FieldHansenObjects<T extends RealFieldElement<T>> {
/** Maximum power of the eccentricity to use in Hansen coefficient Kernel expansion. */
private int maxHansen;
/** A two dimensional array that contains the objects needed to build the Hansen coefficients. <br/>
* The indexes are s + maxDegree and j */
private FieldHansenTesseralLinear<T>[][] hansenObjects;
/** Simple constructor.
* @param ratio Ratio of satellite period to central body rotation period
* @param maxEccPow Maximum power of the eccentricity to use in summation over s.
* @param resOrders List of resonant orders
* @param type type of the elements used during the propagation
* @param field field used by default.
*/
@SuppressWarnings("unchecked")
FieldHansenObjects(final T ratio,
final int maxEccPow,
final List<Integer> resOrders,
final PropagationType type,
final Field<T> field) {
final PropagationType type) {
// Set the maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
maxHansen = maxEccPow / 2;
......@@ -2850,11 +2843,11 @@ public class DSSTTesseral implements DSSTForceModel {
final int n0 = FastMath.max(FastMath.max(2, m), s);
//Create the object for the pair j, s
this.hansenObjects[s + maxDegree][j] = new FieldHansenTesseralLinear<>(maxDegree, s, j, n0, maxHansen, field);
this.hansenObjects[s + maxDegree][j] = new FieldHansenTesseralLinear<>(maxDegree, s, j, n0, maxHansen, ratio.getField());
if (s > 0 && s <= sMin) {
//Also create the object for the pair j, -s
this.hansenObjects[maxDegree - s][j] = new FieldHansenTesseralLinear<>(maxDegree, -s, j, n0, maxHansen, field);
this.hansenObjects[maxDegree - s][j] = new FieldHansenTesseralLinear<>(maxDegree, -s, j, n0, maxHansen, ratio.getField());
}
}
}
......@@ -2866,7 +2859,7 @@ public class DSSTTesseral implements DSSTForceModel {
for (int s = -maxDegree; s <= maxDegree; s++) {
//Compute the n0 value
final int n0 = FastMath.max(2, FastMath.abs(s));
this.hansenObjects[s + maxDegree][j] = new FieldHansenTesseralLinear<>(maxDegree, s, j, n0, maxHansen, field);
this.hansenObjects[s + maxDegree][j] = new FieldHansenTesseralLinear<>(maxDegree, s, j, n0, maxHansen, ratio.getField());
}
}
break;
......
......@@ -16,9 +16,6 @@
*/
package org.orekit.propagation.semianalytical.dsst.forces;
import java.util.ArrayList;
import java.util.List;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
......@@ -52,16 +49,6 @@ class DSSTTesseralContext extends ForceModelContext {
*/
private static final int I = 1;
/** Minimum period for analytically averaged high-order resonant
* central body spherical harmonics in seconds.
*/
private static final double MIN_PERIOD_IN_SECONDS = 864000.;
/** Minimum period for analytically averaged high-order resonant
* central body spherical harmonics in satellite revolutions.
*/
private static final double MIN_PERIOD_IN_SAT_REV = 10.;
/** A = sqrt(μ * a). */
private double A;
......@@ -100,9 +87,6 @@ class DSSTTesseralContext extends ForceModelContext {
/** ecc². */
private double e2;
/** Maximum power of the eccentricity to use in summation over s. */
private int maxEccPow;
/** Keplerian mean motion. */
private double n;
......@@ -112,9 +96,6 @@ class DSSTTesseralContext extends ForceModelContext {
/** Ratio of satellite period to central body rotation period. */
private double ratio;
/** List of resonant orders. */
private final List<Integer> resOrders;
/**
* Simple constructor.
*
......@@ -134,9 +115,6 @@ class DSSTTesseralContext extends ForceModelContext {
super(auxiliaryElements);
this.maxEccPow = 0;
this.resOrders = new ArrayList<Integer>();
final double mu = parameters[0];
// Keplerian Mean Motion
......@@ -179,51 +157,9 @@ class DSSTTesseralContext extends ForceModelContext {
chi = 1. / auxiliaryElements.getB();
chi2 = chi * chi;
// Set the highest power of the eccentricity in the analytical power
// series expansion for the averaged high order resonant central body
// spherical harmonic perturbation
final double e = auxiliaryElements.getEcc();
if (e <= 0.005) {
maxEccPow = 3;
} else if (e <= 0.02) {
maxEccPow = 4;
} else if (e <= 0.1) {
maxEccPow = 7;
} else if (e <= 0.2) {
maxEccPow = 10;
} else if (e <= 0.3) {
maxEccPow = 12;
} else if (e <= 0.4) {
maxEccPow = 15;
} else {
maxEccPow = 20;
}
// Ratio of satellite to central body periods to define resonant terms
ratio = period / bodyPeriod;
// Compute natural resonant terms
final double tolerance = 1. / FastMath.max(MIN_PERIOD_IN_SAT_REV,
MIN_PERIOD_IN_SECONDS / period);
// Search the resonant orders in the tesseral harmonic field
resOrders.clear();
for (int m = 1; m <= provider.getMaxOrder(); m++) {
final double resonance = ratio * m;
final int jComputedRes = (int) FastMath.round(resonance);
if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance - jComputedRes) <= tolerance) {
// Store each resonant index and order
this.resOrders.add(m);
}
}
}
/** Get the list of resonant orders.
* @return resOrders
*/
public List<Integer> getResOrders() {
return resOrders;
}
/** Get ecc².
......@@ -303,13 +239,6 @@ class DSSTTesseralContext extends ForceModelContext {
return roa;
}
/** Get the maximum power of the eccentricity to use in summation over s.
* @return maxEccPow
*/
public int getMaxEccPow() {
return maxEccPow;
}
/** Get the Keplerian period.
* <p>The Keplerian period is computed directly from semi major axis
* and central acceleration constant.</p>
......
......@@ -16,9 +16,6 @@
*/
package org.orekit.propagation.semianalytical.dsst.forces;
import java.util.ArrayList;
import java.util.List;
import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
......@@ -54,16 +51,6 @@ class FieldDSSTTesseralContext<T extends RealFieldElement<T>> extends FieldForce
*/
private static final int I = 1;
/** Minimum period for analytically averaged high-order resonant
* central body spherical harmonics in seconds.
*/
private static final double MIN_PERIOD_IN_SECONDS = 864000.;
/** Minimum period for analytically averaged high-order resonant
* central body spherical harmonics in satellite revolutions.
*/
private static final double MIN_PERIOD_IN_SAT_REV = 10.;
/** A = sqrt(μ * a). */
private T A;
......@@ -108,15 +95,9 @@ class FieldDSSTTesseralContext<T extends RealFieldElement<T>> extends FieldForce
/** Keplerian period. */
private T period;
/** Maximum power of the eccentricity to use in summation over s. */
private int maxEccPow;
/** Ratio of satellite period to central body rotation period. */
private T ratio;
/** List of resonant orders. */
private final List<Integer> resOrders;
/**
* Simple constructor.
*
......@@ -139,9 +120,6 @@ class FieldDSSTTesseralContext<T extends RealFieldElement<T>> extends FieldForce
final Field<T> field = auxiliaryElements.getDate().getField();
final T zero = field.getZero();
this.maxEccPow = 0;
this.resOrders = new ArrayList<Integer>();
final T mu = parameters[0];
// Keplerian mean motion
......@@ -184,51 +162,9 @@ class FieldDSSTTesseralContext<T extends RealFieldElement<T>> extends FieldForce
chi = auxiliaryElements.getB().reciprocal();
chi2 = chi.multiply(chi);
// Set the highest power of the eccentricity in the analytical power
// series expansion for the averaged high order resonant central body
// spherical harmonic perturbation
final T e = auxiliaryElements.getEcc();
if (e.getReal() <= 0.005) {
maxEccPow = 3;
} else if (e.getReal() <= 0.02) {
maxEccPow = 4;
} else if (e.getReal() <= 0.1) {
maxEccPow = 7;
} else if (e.getReal() <= 0.2) {
maxEccPow = 10;
} else if (e.getReal() <= 0.3) {
maxEccPow = 12;
} else if (e.getReal() <= 0.4) {
maxEccPow = 15;
} else {
maxEccPow = 20;
}
// Ratio of satellite to central body periods to define resonant terms
ratio = period.divide(bodyPeriod);
// Compute natural resonant terms
final T tolerance = FastMath.max(zero.add(MIN_PERIOD_IN_SAT_REV),
period.divide(MIN_PERIOD_IN_SECONDS).reciprocal()).reciprocal();
// Search the resonant orders in the tesseral harmonic field
resOrders.clear();
for (int m = 1; m <= provider.getMaxOrder(); m++) {
final T resonance = ratio.multiply(m);
final int jComputedRes = (int) FastMath.round(resonance);
if (jComputedRes > 0 && jComputedRes <= maxFrequencyShortPeriodics && FastMath.abs(resonance.subtract(jComputedRes)).getReal() <= tolerance.getReal()) {
// Store each resonant index and order
this.resOrders.add(m);
}
}
}
/** Get the list of resonant orders.
* @return resOrders
*/
public List<Integer> getResOrders() {
return resOrders;
}
/** Get ecc².
......@@ -308,13 +244,6 @@ class FieldDSSTTesseralContext<T extends RealFieldElement<T>> extends FieldForce
return roa;