Commit 94b45f8b authored by Luc Maisonobe's avatar Luc Maisonobe

Merge branch 'develop' into bug-fix-setmu

Conflicts:
	src/test/java/org/orekit/estimation/leastsquares/BatchLSEstimatorTest.java
parents aaf623c5 3708181b
Bruno Revelin <bruno.revelin@c-s.fr> <Bruno.Revelin@c-s.fr>
Evan Ward <evan.ward@nrl.navy.mil> <evan.ward@nrl.navy.mil>
Evan Ward <evan.ward@nrl.navy.mil> <nomail@nowhere.org>
Fabien Maussion <fabien@orekit.org> <fabien@orekit.org>
Francesco Rocca <francesco.rocca@telespazio.com> <>
Luc Maisonobe <luc@orekit.org> <Luc.Maisonobe@free.fr>
Luc Maisonobe <luc@orekit.org> <Luc.Maisonobe@orekit.org>
Luc Maisonobe <luc@orekit.org> <maisonobe@users.noreply.github.com>
Pascal Parraud <pascal@orekit.org> <pascal@orekit.org>
Romain di Costanzo <romain.di-costanzo@c-s.fr> <romain.di-costanzo@c-s.fr>
Steven <sljkwsk@yahoo.com> <sljkwsk@yahoo.com>
Thierry Ceolin <thierry@orekit.org> <thierry.ceolin@c-s.fr>
Thierry Ceolin <thierry@orekit.org> <thierry@orekit.org>
Thomas Neidhart <thomas.neidhart@gmail.com> <thomas.neidhart@gmail.com>
Véronique Pommier-Maurussane <vero@orekit.org> <vero@orekit.org>
Maxime Journot <maxime.journot@c-s.fr> <Maxime.Journot@c-s.fr>
Andrea Antolino <andrea.antolino@c-s.fr> <Andrea.Antolino@c-s.fr>
Joris Olympio <joris.olympio@c-s.fr> <joris.olympio@c-s.fr>
Nicolas Bernard <nicolas.bernard@c-s.fr> <nicolas.bernard@c-s.fr>
Bruno Revelin <bruno.revelin@c-s.fr> Bruno Revelin <Bruno.Revelin@c-s.fr>
Evan Ward <evan.ward@nrl.navy.mil> Evan Ward <evan.ward@nrl.navy.mil>
Evan Ward <evan.ward@nrl.navy.mil> Evan Ward <nomail@nowhere.org>
Fabien Maussion <fabien@orekit.org> Fabien Maussion <fabien@orekit.org>
Francesco Rocca <francesco.rocca@telespazio.com> Francesco Rocca <>
Luc Maisonobe <luc@orekit.org> Luc Maisonobe <Luc.Maisonobe@free.fr>
Luc Maisonobe <luc@orekit.org> Luc Maisonobe <Luc.Maisonobe@orekit.org>
Luc Maisonobe <luc@orekit.org> Luc Maisonobe <luc@orekit.org>
Luc Maisonobe <luc@orekit.org> maisonobe <maisonobe@users.noreply.github.com>
Pascal Parraud <pascal@orekit.org> Pascal Parraud <pascal@orekit.org>
Romain di Costanzo <romain.di-costanzo@c-s.fr> rdicosta <romain.di-costanzo@c-s.fr>
Steven <sljkwsk@yahoo.com> Steven <sljkwsk@yahoo.com>
Thierry Ceolin <thierry.ceolin@c-s.fr> Thierry Ceolin <thierry@orekit.org>
Thierry Ceolin <thierry.ceolin@c-s.fr> thierry ceolin <thierry@orekit.org>
Thierry Ceolin <thierry.ceolin@c-s.fr> thierry ceolin <thierry.ceolin@c-s.fr>
Thomas Neidhart <thomas.neidhart@gmail.com> Thomas Neidhart <thomas.neidhart@gmail.com>
Véronique Pommier-Maurussane <vero@orekit.org> Véronique Pommier-Maurussane <vero@orekit.org>
Nicolas Bernard <nicolas.bernard@c-s.fr> Nicola Bernard <nicolas.bernard@c-s.fr>
Romain di Costanzo <romain.di-costanzo@c-s.fr> Romain Di Costanzo <romain.di-costanzo@c-s.fr>
Joris Olympio <joris.olympio@c-s.fr> jolympio <joris.olympio@c-s.fr>
Piotr Listkiewicz <piotr.listkiewicz@gmail.com> liscju <piotr.listkiewicz@gmail.com>
......@@ -23,7 +23,9 @@ import java.util.List;
import java.util.Map;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.exception.MathRuntimeException;
import org.hipparchus.linear.RealMatrix;
import org.hipparchus.linear.RealVector;
import org.hipparchus.optim.ConvergenceChecker;
import org.hipparchus.optim.nonlinear.vector.leastsquares.EvaluationRmsChecker;
......@@ -447,12 +449,66 @@ public class BatchLSEstimator {
* The {@link Optimum} object contains detailed elements (covariance matrix, estimated
* parameters standard deviation, weighted Jacobian, RMS, χ², residuals and more).
* </p>
* <p>
* Beware that the returned object is the raw view from the underlying mathematical
* library. At this ral level, parameters have {@link ParameterDriver#getNormalizedValue()
* normalized values} whereas the space flight parameters have {@link ParameterDriver#getValue()
* physical values} with their units. So there are {@link ParameterDriver#getScale() scaling
* factors} to apply when using these elements.
* </p>
* @return optimum found after last call to {@link #estimate()}
*/
public Optimum getOptimum() {
return optimum;
}
/** Get the covariances matrix in space flight dynamics physical units.
* <p>
* This method retrieve the {@link
* org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation#getCovariances(double)
* covariances} from the [@link {@link #getOptimum() optimum} and applies the scaling factors
* to it in order to convert it from raw normalized values back to physical values.
* </p>
* @param threshold threshold to identify matrix singularity
* @return covariances matrix in space flight dynamics physical units
* @exception OrekitException if the covariance matrix cannot be computed (singular problem).
* @since 9.1
*/
public RealMatrix getPhysicalCovariances(final double threshold)
throws OrekitException {
final RealMatrix covariances;
try {
// get the normalized matrix
covariances = optimum.getCovariances(threshold).copy();
} catch (MathIllegalArgumentException miae) {
// the problem is singular
throw new OrekitException(miae);
}
// retrieve the scaling factors
final double[] scale = new double[covariances.getRowDimension()];
int index = 0;
for (final ParameterDriver driver : getOrbitalParametersDrivers(true).getDrivers()) {
scale[index++] = driver.getScale();
}
for (final ParameterDriver driver : getPropagatorParametersDrivers(true).getDrivers()) {
scale[index++] = driver.getScale();
}
for (final ParameterDriver driver : getMeasurementsParametersDrivers(true).getDrivers()) {
scale[index++] = driver.getScale();
}
// unnormalize the matrix, to retrieve physical covariances
for (int i = 0; i < covariances.getRowDimension(); ++i) {
for (int j = 0; j < covariances.getColumnDimension(); ++j) {
covariances.setEntry(i, j, scale[i] * scale[j] * covariances.getEntry(i, j));
}
}
return covariances;
}
/** Get the number of iterations used for last estimation.
* @return number of iterations used for last estimation
* @see #setMaxIterations(int)
......
......@@ -314,7 +314,7 @@ public abstract class AbstractMeasurement<T extends ObservedMeasurement<T>>
* {@code firstDerivative + 1} and {@code firstDerivative + 2}.
* The velocity will correspond to variables {@code firstDerivative + 3},
* {@code firstDerivative + 4} and {@code firstDerivative + 5}.
* The acceleration will corredpon to constants.
* The acceleration will correspond to constants.
* </p>
* @param state state of the satellite considered
* @param firstDerivative index of the first derivative
......
......@@ -16,6 +16,8 @@
*/
package org.orekit.forces.maneuvers;
import java.util.Map;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.attitudes.Attitude;
......@@ -230,8 +232,13 @@ public class ImpulseManeuver<T extends EventDetector> extends AbstractDetector<I
final double newMass = oldState.getMass() * FastMath.exp(-sign * deltaV.getNorm() / im.vExhaust);
// pack everything in a new state
return new SpacecraftState(oldState.getOrbit().getType().convertType(newOrbit),
attitude, newMass);
SpacecraftState newState = new SpacecraftState(oldState.getOrbit().getType().convertType(newOrbit),
attitude, newMass);
for (final Map.Entry<String, double[]> entry : oldState.getAdditionalStates().entrySet()) {
newState = newState.addAdditionalState(entry.getKey(), entry.getValue());
}
return newState;
}
......
......@@ -84,7 +84,7 @@ public class MariniMurrayModel implements TroposphericModel {
@Override
public double pathDelay(final double elevation, final double height) {
final double A = 0.002357 * P0 + 0.000141 * e0;
final double K = 1.163 * 0.00968 * FastMath.cos(2 * latitude) - 0.00104 * T0 + 0.00001435 * P0;
final double K = 1.163 - 0.00968 * FastMath.cos(2 * latitude) - 0.00104 * T0 + 0.00001435 * P0;
final double B = (1.084 * 1e-8) * P0 * T0 * K + (4.734 * 1e-8) * P0 * (P0 / T0) * (2 * K) / (3 * K - 1);
final double flambda = getLaserFrequencyParameter();
......@@ -134,7 +134,7 @@ public class MariniMurrayModel implements TroposphericModel {
// enhancement factor, equation (4) of reference paper
final double fw = 1.00062 + (3.14 * 1e-6) * P0 + (5.6 * 1e-7) * FastMath.pow(T0 - 273.15, 2);
final double e = rh * fw * es;
final double e = (rh / 100.0) * fw * es;
return e;
}
}
......@@ -234,6 +234,7 @@ discrete events. Here is a short list of the features offered by the library:</p
<ul>
<li>tropospheric delay (modified Saastamoinen)</li>
<li>tropospheric refraction correction angle (Recommendation ITU-R P.834-7 and Saemundssen's formula quoted by Meeus)</li>
<li>tropospheric model for laser ranging (Marini-Murray)</li>
<li>Klobuchar ionospheric model (including parsing α and β coefficients from University of Bern Astronomical Institute files)</li>
<li>geomagnetic field (WMM, IGRF)</li>
<li>geoid model from any gravity field</li>
......
......@@ -19,7 +19,9 @@ package org.orekit.propagation.analytical;
import java.io.NotSerializableException;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.Map;
import java.util.SortedSet;
import org.orekit.attitudes.Attitude;
......@@ -120,7 +122,7 @@ public class KeplerianPropagator extends AbstractAnalyticalPropagator implements
getAttitudeProvider().getAttitude(initialOrbit,
initialOrbit.getDate(),
initialOrbit.getFrame()),
mass, mu);
mass, mu, Collections.emptyMap());
states = new TimeSpanMap<SpacecraftState>(initialState);
super.resetInitialState(initialState);
......@@ -135,16 +137,21 @@ public class KeplerianPropagator extends AbstractAnalyticalPropagator implements
* @param attitude current attitude
* @param mass current mass
* @param mu gravity coefficient to use
* @param additionalStates additional states
* @return fixed orbit
*/
private SpacecraftState fixState(final Orbit orbit, final Attitude attitude, final double mass,
final double mu) {
final double mu, final Map<String, double[]> additionalStates) {
final OrbitType type = orbit.getType();
final double[] stateVector = new double[6];
type.mapOrbitToArray(orbit, PositionAngle.TRUE, stateVector, null);
final Orbit fixedOrbit = type.mapArrayToOrbit(stateVector, null, PositionAngle.TRUE,
orbit.getDate(), mu, orbit.getFrame());
return new SpacecraftState(fixedOrbit, attitude, mass);
SpacecraftState fixedState = new SpacecraftState(fixedOrbit, attitude, mass);
for (final Map.Entry<String, double[]> entry : additionalStates.entrySet()) {
fixedState = fixedState.addAdditionalState(entry.getKey(), entry.getValue());
}
return fixedState;
}
/** {@inheritDoc} */
......@@ -156,7 +163,8 @@ public class KeplerianPropagator extends AbstractAnalyticalPropagator implements
final SpacecraftState fixedState = fixState(state.getOrbit(),
state.getAttitude(),
state.getMass(),
mu);
mu,
state.getAdditionalStates());
initialState = fixedState;
states = new TimeSpanMap<SpacecraftState>(initialState);
......
......@@ -35,7 +35,7 @@ public class DormandPrince853IntegratorBuilder implements ODEIntegratorBuilder {
/** Maximum step size (s). */
private final double maxStep;
/** Minimum step size (s). */
/** Position error (m). */
private final double dP;
/** Build a new instance.
......
......@@ -18,6 +18,7 @@ package org.orekit.propagation.numerical;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import org.hipparchus.Field;
......@@ -29,6 +30,7 @@ import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.FieldAttitude;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitIllegalArgumentException;
import org.orekit.errors.OrekitInternalError;
import org.orekit.errors.OrekitMessages;
import org.orekit.forces.ForceModel;
import org.orekit.forces.gravity.NewtonianAttraction;
......@@ -44,6 +46,8 @@ import org.orekit.propagation.integration.FieldStateMapper;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.utils.FieldPVCoordinates;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.ParameterObserver;
import org.orekit.utils.TimeStampedFieldPVCoordinates;
/** This class propagates {@link org.orekit.orbits.FieldOrbit orbits} using
......@@ -141,9 +145,6 @@ import org.orekit.utils.TimeStampedFieldPVCoordinates;
*/
public class FieldNumericalPropagator<T extends RealFieldElement<T>> extends FieldAbstractIntegratedPropagator<T> {
/** Central body attraction. */
private NewtonianAttraction newtonianAttraction;
/** Force models used during the extrapolation of the FieldOrbit<T>, without Jacobians. */
private final List<ForceModel> forceModels;
......@@ -175,8 +176,25 @@ public class FieldNumericalPropagator<T extends RealFieldElement<T>> extends Fie
* @see #addForceModel(ForceModel)
*/
public void setMu(final double mu) {
addForceModel(new NewtonianAttraction(mu));
}
/** Set the central attraction coefficient μ only in upper class.
* @param mu central attraction coefficient (m³/s²)
*/
private void superSetMu(final double mu) {
super.setMu(mu);
newtonianAttraction = new NewtonianAttraction(mu);
}
/** Check if Newtonian attraction force model is available.
* <p>
* Newtonian attraction is always the last force model in the list.
* </p>
* @return true if Newtonian attraction force model is available
*/
private boolean hasNewtonianAttraction() {
final int last = forceModels.size() - 1;
return last >= 0 && forceModels.get(last) instanceof NewtonianAttraction;
}
/** Add a force model to the global perturbation model.
......@@ -187,7 +205,43 @@ public class FieldNumericalPropagator<T extends RealFieldElement<T>> extends Fie
* @see #setMu(double)
*/
public void addForceModel(final ForceModel model) {
forceModels.add(model);
if (model instanceof NewtonianAttraction) {
// we want to add the central attraction force model
try {
// ensure we are notified of any mu change
model.getParametersDrivers()[0].addObserver(new ParameterObserver() {
/** {@inheritDoc} */
@Override
public void valueChanged(final double previousValue, final ParameterDriver driver) {
superSetMu(driver.getValue());
}
});
} catch (OrekitException oe) {
// this should never happen
throw new OrekitInternalError(oe);
}
if (hasNewtonianAttraction()) {
// there is already a central attraction model, replace it
forceModels.set(forceModels.size() - 1, model);
} else {
// there are no central attraction model yet, add it at the end of the list
forceModels.add(model);
}
} else {
// we want to add a perturbing force model
if (hasNewtonianAttraction()) {
// insert the new force model before Newtonian attraction,
// which should always be the last one in the list
forceModels.add(forceModels.size() - 1, model);
} else {
// we only have perturbing force models up to now, just append at the end of the list
forceModels.add(model);
}
}
}
/** Remove all perturbing force models from the global perturbation model.
......@@ -200,22 +254,41 @@ public class FieldNumericalPropagator<T extends RealFieldElement<T>> extends Fie
forceModels.clear();
}
/** Get all the force models, perturbing forces and Newtonian attraction included.
* @return list of perturbing force models, with Newtonian attraction being the
* last one
* @see #addForceModel(ForceModel)
* @see #setMu(double)
* @since 9.1
*/
public List<ForceModel> getAllForceModels() {
return Collections.unmodifiableList(forceModels);
}
/** Get perturbing force models list.
* @return list of perturbing force models
* @see #addForceModel(ForceModel)
* @see #getNewtonianAttractionForceModel()
* @deprecated as of 9.1, this method is deprecated, the perturbing
* force models are retrieved together with the Newtonian attraction
* by calling {@link #getAllForceModels()}
*/
@Deprecated
public List<ForceModel> getForceModels() {
return forceModels;
return hasNewtonianAttraction() ? forceModels.subList(0, forceModels.size() - 1) : forceModels;
}
/** Get the Newtonian attraction from the central body force model.
* @return Newtonian attraction force model
* @see #setMu(double)
* @see #getForceModels()
* @deprecated as of 9.1, this method is deprecated, the Newtonian
* attraction force model (if any) is the last in the {@link #getAllForceModels()}
*/
@Deprecated
public NewtonianAttraction getNewtonianAttractionForceModel() {
return newtonianAttraction;
final int last = forceModels.size() - 1;
if (last >= 0 && forceModels.get(last) instanceof NewtonianAttraction) {
return (NewtonianAttraction) forceModels.get(last);
} else {
return null;
}
}
/** Set propagation orbit type.
......@@ -263,7 +336,7 @@ public class FieldNumericalPropagator<T extends RealFieldElement<T>> extends Fie
/** {@inheritDoc} */
public void resetInitialState(final FieldSpacecraftState<T> state) throws OrekitException {
super.resetInitialState(state);
if (newtonianAttraction == null) {
if (!hasNewtonianAttraction()) {
setMu(state.getMu());
}
setStartDate(state.getDate());
......@@ -376,12 +449,14 @@ public class FieldNumericalPropagator<T extends RealFieldElement<T>> extends Fie
orbit = state.getOrbit();
Arrays.fill(yDot, zero);
orbit.getJacobianWrtCartesian(getPositionAngleType(), jacobian);
// compute the contributions of all perturbing forces
// compute the contributions of all perturbing forces,
// using the Kepler contribution at the end since
// NewtonianAttraction is always the last instance in the list
for (final ForceModel forceModel : forceModels) {
forceModel.addContribution(state, this);
}
// finalize derivatives by adding the Kepler contribution
newtonianAttraction.addContribution(state, this);
return yDot.clone();
}
......
......@@ -194,6 +194,7 @@
* tropospheric delay (modified Saastamoinen)
* tropospheric refraction correction angle (Recommendation ITU-R P.834-7 and Saemundssen's formula quoted by Meeus)
* tropospheric model for laser ranging (Marini-Murray)
* Klobuchar ionospheric model (including parsing α and β coefficients from University of Bern Astronomical Institute files)
* geomagnetic field (WMM, IGRF)
* geoid model from any gravity field
......
......@@ -21,6 +21,34 @@
</properties>
<body>
<release version="9.1" date="TBD" description="TBD">
<action dev="luc" type="fix">
Avoid intermixed ChangeForwarder instances calling each other.
Fixes issue #360.
</action>
<action dev="maxime" type="fix">
Modified the way the propagation parameter drivers are mapped in the
Jacobian matrix in class "Model".
Added a test for multi-sat orbit determination with estimated
propagation parameters (µ and SRP coefficients).
Fixes issue #354.
</action>
<action dev="luc" type="fix">
Added a convenience method to retrieve covariance matrix in
physical units in orbit determination.
Fixes issue #353.
</action>
<action dev="luc" type="fix" due-to="Rongwang Li">
Fixed two errors in Marini-Murray model implementation.
Fixes issue #352.
</action>
<action dev="luc" type="fix">
Prevent duplicated Newtonian attraction in FieldNumericalPropagator.
Fixes issue #350.
</action>
<action dev="luc" type="fix">
Copy additional states through impulse maneuvers.
Fixes issue #349.
</action>
<action dev="luc" type="fix">
Removed unused construction parameters in ShiftingTransformProvider
and InterpolatingTransformProvider.
......
......@@ -21,6 +21,7 @@ import java.util.List;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.linear.RealMatrix;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
import org.junit.Assert;
......@@ -91,6 +92,14 @@ public class BatchLSEstimatorTest {
0.0, 1.4e-8,
0.0, 6.3e-12);
RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
Assert.assertEquals(6, normalizedCovariances.getRowDimension());
Assert.assertEquals(6, normalizedCovariances.getColumnDimension());
Assert.assertEquals(6, physicalCovariances.getRowDimension());
Assert.assertEquals(6, physicalCovariances.getColumnDimension());
Assert.assertEquals(0.00258, physicalCovariances.getEntry(0, 0), 1.0e-5);
}
@Test
......@@ -541,7 +550,7 @@ public class BatchLSEstimatorTest {
for (final ObservedMeasurement<?> range : r1) {
estimator.addMeasurement(range);
}
estimator.setParametersConvergenceThreshold(1.0e-2);
estimator.setParametersConvergenceThreshold(1.0e-3);
estimator.setMaxIterations(10);
estimator.setMaxEvaluations(20);
estimator.setObserver(new BatchLSObserver() {
......@@ -623,7 +632,7 @@ public class BatchLSEstimatorTest {
Assert.assertEquals(0.0,
Vector3D.distance(closeOrbit.getPVCoordinates().getPosition(),
determined.getPVCoordinates().getPosition()),
5.8e-6);
2.0e-6);
Assert.assertEquals(0.0,
Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(),
determined.getPVCoordinates().getVelocity()),
......
......@@ -37,6 +37,7 @@ import java.util.TreeSet;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.linear.RealMatrix;
import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer;
import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.Decomposition;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer;
......@@ -295,6 +296,22 @@ public class OrbitDeterminationTest {
Assert.assertEquals(RefStatEle[2], odsatW3.getElevStat().getMean(), angleAccuracy);
Assert.assertEquals(RefStatEle[3], odsatW3.getElevStat().getStandardDeviation(), angleAccuracy);
RealMatrix covariances = odsatW3.getCovariances();
Assert.assertEquals(28, covariances.getRowDimension());
Assert.assertEquals(28, covariances.getColumnDimension());
// drag coefficient variance
Assert.assertEquals(0.77754, covariances.getEntry(6, 6), 1.0e-5);
// leak-X constant term variance
Assert.assertEquals(2.1016e-12, covariances.getEntry(7, 7), 1.0e-16);
// leak-Y constant term variance
Assert.assertEquals(2.8869e-11, covariances.getEntry(9, 9), 1.0e-15);
// leak-Z constant term variance
Assert.assertEquals(9.0489e-11, covariances.getEntry(11, 11), 1.0e-15);
}
private class ResultOD {
......@@ -304,14 +321,16 @@ public class OrbitDeterminationTest {
private StreamingStatistics rangeStat;
private StreamingStatistics azimStat;
private StreamingStatistics elevStat;
private ParameterDriversList propagatorParameters ;
private ParameterDriversList measurementsParameters;
private ParameterDriversList propagatorParameters ;
private ParameterDriversList measurementsParameters;
private RealMatrix covariances;
ResultOD(ParameterDriversList propagatorParameters,
ParameterDriversList measurementsParameters,
int numberOfIteration, int numberOfEvaluation, TimeStampedPVCoordinates estimatedPV,
StreamingStatistics rangeStat, StreamingStatistics rangeRateStat,
StreamingStatistics azimStat, StreamingStatistics elevStat,
StreamingStatistics posStat, StreamingStatistics velStat) {
StreamingStatistics posStat, StreamingStatistics velStat,
RealMatrix covariances) {
this.propagatorParameters = propagatorParameters;
this.measurementsParameters = measurementsParameters;
......@@ -321,39 +340,36 @@ public class OrbitDeterminationTest {
this.rangeStat = rangeStat;
this.azimStat = azimStat;
this.elevStat = elevStat;
this.covariances = covariances;
}
public int getNumberOfIteration() {
return numberOfIteration;
}
public int getNumberOfIteration() {
return numberOfIteration;
}
public int getNumberOfEvaluation() {
return numberOfEvaluation;
}
public PVCoordinates getEstimatedPV() {
return estimatedPV;
}
public StreamingStatistics getRangeStat() {
return rangeStat;
}
public int getNumberOfEvaluation() {
return numberOfEvaluation;
}
public PVCoordinates getEstimatedPV() {
return estimatedPV;
}
public StreamingStatistics getAzimStat() {
return azimStat;
}
public StreamingStatistics getRangeStat() {
return rangeStat;
}
public StreamingStatistics getAzimStat() {
return azimStat;
}
public StreamingStatistics getElevStat() {
return elevStat;
}
public StreamingStatistics getElevStat() {
return elevStat;
}
public RealMatrix getCovariances() {
return covariances;
}
}
......@@ -514,7 +530,8 @@ public class OrbitDeterminationTest {
estimator.getIterationsCount(), estimator.getEvaluationsCount(), estimated.getPVCoordinates(),
rangeLog.createStatisticsSummary(), rangeRateLog.createStatisticsSummary(),
azimuthLog.createStatisticsSummary(), elevationLog.createStatisticsSummary(),
positionLog.createStatisticsSummary(), velocityLog.createStatisticsSummary());
positionLog.createStatisticsSummary(), velocityLog.createStatisticsSummary(),
estimator.getPhysicalCovariances(1.0e-10));
}
/** Sort parameters changes.
......
......@@ -28,6 +28,7 @@ import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.linear.RealMatrix;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
import org.hipparchus.random.RandomGenerator;
import org.hipparchus.random.Well19937a;
......@@ -134,6 +135,17 @@ public class GroundStationTest {
Assert.assertEquals(reference.getLatitude(), result.getLatitude(), 6.5e-15);
Assert.assertEquals(reference.getLongitude(), result.getLongitude(), 1.5e-14);
Assert.assertEquals(reference.getAltitude(), result.getAltitude(), 1.6e-7);
RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
Assert.assertEquals(9, normalizedCovariances.getRowDimension());
Assert.assertEquals(9, normalizedCovariances.getColumnDimension());
Assert.assertEquals(9, physicalCovariances.getRowDimension());
Assert.assertEquals(9, physicalCovariances.getColumnDimension());
Assert.assertEquals(0.55431, physicalCovariances.getEntry(6, 6), 1.0e-5);
Assert.assertEquals(0.22694, physicalCovariances.getEntry(7, 7), 1.0e-5);
Assert.assertEquals(0.13106, physicalCovariances.getEntry(8, 8), 1.0e-5);
}
@Test
......
......@@ -21,11 +21,13 @@ import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.RotationConvention;
import org.hipparchus.geometry.euclidean.threed.RotationOrder;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.hipparchus.util.FastMath;
import org.junit.Assert;
import org.junit.Before;
import org.junit.Test;
import org.orekit.Utils;
import org.orekit.attitudes.Attitude;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.InertialProvider;
import org.orekit.attitudes.LofOffset;
......@@ -41,6 +43,8 @@ import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.analytical.KeplerianPropagator;
import org.orekit.propagation.events.DateDetector;
import org.orekit.propagation.events.NodeDetector;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.propagation.numerical.PartialDerivativesEquations;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.DateComponents;
import org.orekit.time.TimeComponents;
......@@ -190,6 +194,115 @@ public class ImpulseManeuverTest {
}
@Test
public void testAdditionalStateKeplerian() throws OrekitException {
final double mu = CelestialBodyFactory.getEarth().getGM();
final double initialX = 7100e3;
final double initialY = 0.0;
final double initialZ = 1300e3;
final double initialVx = 0;
final double initialVy = 8000;
final double initialVz = 1000;