Hi everyone,
I am still experiencing discrepancies between the estimator object itself and what I recompute with a separate propagator. I can’t identify the source as one of things I mentioned previously, but it seems linked to drag.
I’ve created a test case in Java that hopefully you guys can run.
Imports:
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.nonstiff.DormandPrince54Integrator;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.estimation.leastsquares.BatchLSEstimator;
import org.orekit.estimation.measurements.EstimatedMeasurement;
import org.orekit.estimation.measurements.ObservableSatellite;
import org.orekit.estimation.measurements.ObservedMeasurement;
import org.orekit.estimation.measurements.PV;
import org.orekit.forces.drag.DragForce;
import org.orekit.forces.drag.IsotropicDrag;
import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.ICGEMFormatReader;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.ITRFVersion;
import org.orekit.models.earth.Geoid;
import org.orekit.models.earth.ReferenceEllipsoid;
import org.orekit.models.earth.atmosphere.NRLMSISE00;
import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation;
import org.orekit.orbits.CartesianOrbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.conversion.DormandPrince54IntegratorBuilder;
import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
import java.io.File;
import java.util.Locale;
import java.util.Map;
Main:
public static void main(final String[] args) {
try {
// configure Orekit
final File home = new File(System.getProperty("user.home"));
final File orekitData = new File(home, "Dev//orekit-data");
if (!orekitData.exists()) {
System.err.format(Locale.US, "Failed to find %s folder%n",
orekitData.getAbsolutePath());
System.err.format(Locale.US, "You need to download %s from %s, unzip it in %s and rename it 'orekit-data' for this tutorial to work%n",
"orekit-data-master.zip", "https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip",
home.getAbsolutePath());
System.exit(1);
}
final DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
manager.addProvider(new DirectoryCrawler(orekitData));
// Initial date in UTC time scale
final TimeScale utc = TimeScalesFactory.getUTC();
final AbsoluteDate initialDate = new AbsoluteDate(2021, 10, 2, 0, 0, 00.000, utc);
// Frames
final Frame inertialFrame = FramesFactory.getGCRF();
final Frame itrf = FramesFactory.getITRF(ITRFVersion.ITRF_2014, IERSConventions.IERS_2010, false);
// Initial conditions
final double r0 = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 500e3;
final double v0 = FastMath.sqrt(Constants.EGM96_EARTH_MU / r0);
final double mass = 10.;
final double cs = 1.;
final double dragCoef = 2.;
final Vector3D pos = new Vector3D(r0, 0., 0.);
final Vector3D vel = new Vector3D(0., v0, 0.);
final PVCoordinates pv = new PVCoordinates(pos, vel);
// Forces
GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader(GravityFieldFactory.ICGEM_FILENAME, false));
final int degree = 2;
final int order = 1;
final NormalizedSphericalHarmonicsProvider gravityFields = GravityFieldFactory.getNormalizedProvider(degree, order);
final Geoid geoid = new Geoid(gravityFields, ReferenceEllipsoid.getWgs84(itrf));
final HolmesFeatherstoneAttractionModel holmes = new HolmesFeatherstoneAttractionModel(geoid.getBodyFrame(),
gravityFields);
final MarshallSolarActivityFutureEstimation atmParams = new MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
final NRLMSISE00 atmo = new NRLMSISE00(atmParams, CelestialBodyFactory.getSun(), geoid);
final double coeffMin = 0.;
final double coeffMax = 1000.;
final IsotropicDrag drag = new IsotropicDrag(cs, dragCoef, coeffMin, coeffMax);
final DragForce dragForce = new DragForce(atmo, drag);
final double mu = holmes.getMu();
final CartesianOrbit orbit = new CartesianOrbit(pv, inertialFrame, initialDate, mu);
final SpacecraftState spacecraftState = new SpacecraftState(orbit, mass);
// Integrators
final OrbitType orbitType = OrbitType.CARTESIAN;
final PositionAngle positionAngle = PositionAngle.TRUE;
final double minStep = 0.1;
final double maxStep = 60.;
final double positionTolerance = 0.1;
final double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, orbit, orbitType);
final DormandPrince54Integrator dp = new DormandPrince54Integrator(minStep, maxStep,
tolerances[0], tolerances[1]);
// Create Cartesian pseudo-measurements by propagating without drag
final double dt = 10000.;
final AbsoluteDate date2 = initialDate.shiftedBy(dt);
final AbsoluteDate date3 = date2.shiftedBy(dt);
final NumericalPropagator propReal = new NumericalPropagator(dp);
propReal.setInitialState(spacecraftState);
propReal.setOrbitType(orbitType);
propReal.addForceModel(holmes);
propReal.setResetAtEnd(false);
final PVCoordinates pv2 = propReal.propagate(date2).getPVCoordinates(inertialFrame);
final PVCoordinates pv3 = propReal.propagate(date3).getPVCoordinates(inertialFrame);
// Guess for estimator (real orbit with error introduced)
final PVCoordinates modifiedPv = new PVCoordinates(pos.add(new Vector3D(1e3, -2e3, 4e2)),
vel.add(new Vector3D(-3e-4, 5e-5, 1e-4)));
final CartesianOrbit modifiedOrbit = new CartesianOrbit(modifiedPv, inertialFrame, initialDate, mu);
// Propagator builder for estimator
final DormandPrince54IntegratorBuilder dpBuilder = new DormandPrince54IntegratorBuilder(minStep, maxStep,
positionTolerance);
final double positionScale = 1.;
final NumericalPropagatorBuilder propBuilder = new NumericalPropagatorBuilder(modifiedOrbit, dpBuilder,
positionAngle, positionScale);
propBuilder.addForceModel(holmes);
propBuilder.addForceModel(dragForce);
// Create estimator
final double weight = 1.;
final double sigmaPos = 1e-3;
final double sigmaVel = 1e-6;
final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
final BatchLSEstimator estimator = new BatchLSEstimator(optimizer, propBuilder);
estimator.setMaxEvaluations(60);
estimator.setMaxIterations(40);
estimator.setParametersConvergenceThreshold(1e-4);
final ObservableSatellite sat = new ObservableSatellite(0);
estimator.addMeasurement(new PV(initialDate, pv.getPosition(), pv.getVelocity(),
sigmaPos, sigmaVel, weight, sat));
estimator.addMeasurement(new PV(date2, pv2.getPosition(), pv2.getVelocity(),
sigmaPos, sigmaVel, weight, sat));
estimator.addMeasurement(new PV(date3, pv3.getPosition(), pv3.getVelocity(),
sigmaPos, sigmaVel, weight, sat));
final String name = "drag coefficient";
estimator.getPropagatorParametersDrivers(false).findByName(name).setSelected(true);
estimator.getPropagatorParametersDrivers(false).findByName(name).setReferenceValue(dragCoef);
// perform OD and post process it
final NumericalPropagator estimatedProp = (NumericalPropagator) estimator.estimate()[0];
estimatedProp.setResetAtEnd(false);
final SpacecraftState estimatedState = estimatedProp.getInitialState();
final PVCoordinates estimatedPV = estimatedState.getPVCoordinates(inertialFrame);
final PVCoordinates pv2Fitted = estimatedProp.propagate(date2).getPVCoordinates(inertialFrame);
final PVCoordinates pv3Fitted = estimatedProp.propagate(date3).getPVCoordinates(inertialFrame);
// Create independent propagator and recompute estimated measurements
final double[][] tolerancesBis = NumericalPropagator.tolerances(positionTolerance,
estimatedState.getOrbit(), estimatedState.getOrbit().getType());
final DormandPrince54Integrator dpBis = new DormandPrince54Integrator(minStep, maxStep,
tolerancesBis[0], tolerancesBis[1]);
final NumericalPropagator propBis = new NumericalPropagator(dpBis);
final CartesianOrbit orbitBis = new CartesianOrbit(new PVCoordinates(estimatedPV.getPosition(),
estimatedPV.getVelocity()), inertialFrame, initialDate, estimatedState.getMu());
final SpacecraftState scBis = new SpacecraftState(orbitBis, spacecraftState.getMass());
propBis.setInitialState(scBis);
propBis.setOrbitType(orbitType);
final Geoid geoid2 = new Geoid(gravityFields, ReferenceEllipsoid.getWgs84(itrf));
propBis.addForceModel(new HolmesFeatherstoneAttractionModel(geoid2.getBodyFrame(), gravityFields));
final double fittedCoef = estimator.getPropagatorParametersDrivers(true).findByName(name).getValue();
final IsotropicDrag newDrag = new IsotropicDrag(cs, fittedCoef, coeffMin, coeffMax);
final NRLMSISE00 atmo2 = new NRLMSISE00(atmParams, CelestialBodyFactory.getSun(), geoid2);
propBis.addForceModel(new DragForce(atmo2, newDrag));
propBis.setResetAtEnd(false);
final PVCoordinates pv2Repropagated = propBis.propagate(date2).getPVCoordinates(inertialFrame);
final PVCoordinates pv3Repropagated = propBis.propagate(date3).getPVCoordinates(inertialFrame);
// assess differences between propagators
final double[] diffPos2 = (pv2Fitted.getPosition().subtract(pv2Repropagated.getPosition())).toArray();
final double[] diffVel2 = (pv2Fitted.getVelocity().subtract(pv2Repropagated.getVelocity())).toArray();
final double[] diffPos3 = (pv3Fitted.getPosition().subtract(pv3Repropagated.getPosition())).toArray();
final double[] diffVel3 = (pv3Fitted.getVelocity().subtract(pv3Repropagated.getVelocity())).toArray();
// look at estimated measurements (for debug)
final Map<ObservedMeasurement<?>, EstimatedMeasurement<?>> estimatedMeasurement = estimator.getLastEstimations();
for (EstimatedMeasurement meas: estimatedMeasurement.values()) {
final double[] estimatedValue = meas.getEstimatedValue();
final double[] observedValue = meas.getObservedValue();
}
} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
}
}
Specifically, the issue is that diffPositionX
is non-negligible.
N.B.: if some objects are re-instantiated rather than reused, it is because this is what my original Python code does and I’ve tried to change things as little as possible when transcribing it into Java.
Cheers,
Romain.