use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.
the class EstimationTestUtils method checkFit.
/**
* Checker for batch LS estimator validation
* @param context Context used for the test
* @param estimator Batch LS estimator
* @param iterations Number of iterations expected
* @param evaluations Number of evaluations expected
* @param expectedRMS Expected RMS value
* @param rmsEps Tolerance on expected RMS
* @param expectedMax Expected weighted residual maximum
* @param maxEps Tolerance on weighted residual maximum
* @param expectedDeltaPos Expected position difference between estimated orbit and initial orbit
* @param posEps Tolerance on expected position difference
* @param expectedDeltaVel Expected velocity difference between estimated orbit and initial orbit
* @param velEps Tolerance on expected velocity difference
* @throws OrekitException
*/
public static void checkFit(final Context context, final BatchLSEstimator estimator, final int iterations, final int evaluations, final double expectedRMS, final double rmsEps, final double expectedMax, final double maxEps, final double expectedDeltaPos, final double posEps, final double expectedDeltaVel, final double velEps) throws OrekitException {
final Orbit estimatedOrbit = estimator.estimate()[0].getInitialState().getOrbit();
final Vector3D estimatedPosition = estimatedOrbit.getPVCoordinates().getPosition();
final Vector3D estimatedVelocity = estimatedOrbit.getPVCoordinates().getVelocity();
Assert.assertEquals(iterations, estimator.getIterationsCount());
Assert.assertEquals(evaluations, estimator.getEvaluationsCount());
Optimum optimum = estimator.getOptimum();
Assert.assertEquals(iterations, optimum.getIterations());
Assert.assertEquals(evaluations, optimum.getEvaluations());
int k = 0;
double sum = 0;
double max = 0;
for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
final ObservedMeasurement<?> m = entry.getKey();
final EstimatedMeasurement<?> e = entry.getValue();
final double[] weight = m.getBaseWeight();
final double[] sigma = m.getTheoreticalStandardDeviation();
final double[] observed = m.getObservedValue();
final double[] theoretical = e.getEstimatedValue();
for (int i = 0; i < m.getDimension(); ++i) {
final double weightedResidual = weight[i] * (theoretical[i] - observed[i]) / sigma[i];
++k;
sum += weightedResidual * weightedResidual;
max = FastMath.max(max, FastMath.abs(weightedResidual));
}
}
Assert.assertEquals(expectedRMS, FastMath.sqrt(sum / k), rmsEps);
Assert.assertEquals(expectedMax, max, maxEps);
Assert.assertEquals(expectedDeltaPos, Vector3D.distance(context.initialOrbit.getPVCoordinates().getPosition(), estimatedPosition), posEps);
Assert.assertEquals(expectedDeltaVel, Vector3D.distance(context.initialOrbit.getPVCoordinates().getVelocity(), estimatedVelocity), velEps);
}
use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.
the class IodGibbsTest method testGibbs1.
@Test
public void testGibbs1() throws OrekitException {
final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
final double mu = context.initialOrbit.getMu();
final Frame frame = context.initialOrbit.getFrame();
final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
// create perfect range measurements
final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new PVMeasurementCreator(), 0.0, 1.0, 60.0);
final Vector3D position1 = new Vector3D(measurements.get(0).getObservedValue()[0], measurements.get(0).getObservedValue()[1], measurements.get(0).getObservedValue()[2]);
final PV pv1 = new PV(measurements.get(0).getDate(), position1, Vector3D.ZERO, 0., 0., 1.);
final Vector3D position2 = new Vector3D(measurements.get(1).getObservedValue()[0], measurements.get(1).getObservedValue()[1], measurements.get(1).getObservedValue()[2]);
final PV pv2 = new PV(measurements.get(1).getDate(), position2, Vector3D.ZERO, 0., 0., 1.);
final Vector3D position3 = new Vector3D(measurements.get(2).getObservedValue()[0], measurements.get(2).getObservedValue()[1], measurements.get(2).getObservedValue()[2]);
final PV pv3 = new PV(measurements.get(2).getDate(), position3, Vector3D.ZERO, 0., 0., 1.);
// instantiate the IOD method
final IodGibbs gibbs = new IodGibbs(mu);
final KeplerianOrbit orbit = gibbs.estimate(frame, pv1, pv2, pv3);
Assert.assertEquals(context.initialOrbit.getA(), orbit.getA(), 1.0e-9 * context.initialOrbit.getA());
Assert.assertEquals(context.initialOrbit.getE(), orbit.getE(), 1.0e-9 * context.initialOrbit.getE());
Assert.assertEquals(context.initialOrbit.getI(), orbit.getI(), 1.0e-9 * context.initialOrbit.getI());
}
use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.
the class IodGibbsTest method testGibbs2.
@Test
public void testGibbs2() throws OrekitException {
// test extracted from "Fundamentals of astrodynamics & applications", D. Vallado, 3rd ed, chap Initial Orbit Determination, Exple 7-3, p457
// extraction of the context.
final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
final double mu = context.initialOrbit.getMu();
// initialisation
final IodGibbs gibbs = new IodGibbs(mu);
// Observation vector (EME2000)
final Vector3D posR1 = new Vector3D(0.0, 0.0, 6378137.0);
final Vector3D posR2 = new Vector3D(0.0, -4464696.0, -5102509.0);
final Vector3D posR3 = new Vector3D(0.0, 5740323.0, 3189068);
// epoch corresponding to the observation vector
AbsoluteDate dateRef = new AbsoluteDate(2000, 01, 01, 0, 0, 0, TimeScalesFactory.getUTC());
AbsoluteDate date2 = dateRef.shiftedBy(76.48);
AbsoluteDate date3 = dateRef.shiftedBy(153.04);
// Reference result (cf. Vallado)
final Vector3D velR2 = new Vector3D(0.0, 5531.148, -5191.806);
// Gibbs IOD
final KeplerianOrbit orbit = gibbs.estimate(FramesFactory.getEME2000(), posR1, dateRef, posR2, date2, posR3, date3);
// test
Assert.assertEquals(0.0, orbit.getPVCoordinates().getVelocity().getNorm() - velR2.getNorm(), 1e-3);
}
use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.
the class IodGibbsTest method testGibbs3.
@Test
public void testGibbs3() throws OrekitException {
// test extracted from "Fundamentals of astrodynamics & applications", D. Vallado, 3rd ed, chap Initial Orbit Determination, Exple 7-4, p463
// Remark: the test value in Vallado is performed with an Herrick-Gibbs methods but results are very close with Gibbs method.
// extraction of context
final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
final double mu = context.initialOrbit.getMu();
// Initialisation
final IodGibbs gibbs = new IodGibbs(mu);
// Observations vector (EME2000)
final Vector3D posR1 = new Vector3D(3419855.64, 6019826.02, 2784600.22);
final Vector3D posR2 = new Vector3D(2935911.95, 6326183.24, 2660595.84);
final Vector3D posR3 = new Vector3D(2434952.02, 6597386.74, 2521523.11);
// epoch corresponding to the observation vector
AbsoluteDate dateRef = new AbsoluteDate(2000, 01, 01, 0, 0, 0, TimeScalesFactory.getUTC());
AbsoluteDate date2 = dateRef.shiftedBy(76.48);
AbsoluteDate date3 = dateRef.shiftedBy(153.04);
// Reference result
final Vector3D velR2 = new Vector3D(-6441.632, 3777.625, -1720.582);
// Gibbs IOD
final KeplerianOrbit orbit = gibbs.estimate(FramesFactory.getEME2000(), posR1, dateRef, posR2, date2, posR3, date3);
// test for the norm of the velocity
Assert.assertEquals(0.0, orbit.getPVCoordinates().getVelocity().getNorm() - velR2.getNorm(), 1e-3);
}
use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.
the class IodGoodingTest method testGooding.
@Test
public void testGooding() throws OrekitException {
final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
final double mu = context.initialOrbit.getMu();
final Frame frame = context.initialOrbit.getFrame();
final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
// create perfect range measurements
final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new PVMeasurementCreator(), 0.0, 1.0, 60.0);
// measurement data 1
final int idMeasure1 = 0;
final AbsoluteDate date1 = measurements.get(idMeasure1).getDate();
final Vector3D stapos1 = Vector3D.ZERO;
/*context.stations.get(0) // FIXME we need to access the station of the measurement
.getBaseFrame()
.getPVCoordinates(date1, frame)
.getPosition();*/
final Vector3D position1 = new Vector3D(measurements.get(idMeasure1).getObservedValue()[0], measurements.get(idMeasure1).getObservedValue()[1], measurements.get(idMeasure1).getObservedValue()[2]);
final double r1 = position1.getNorm();
final Vector3D lineOfSight1 = position1.normalize();
// measurement data 2
final int idMeasure2 = 20;
final AbsoluteDate date2 = measurements.get(idMeasure2).getDate();
final Vector3D stapos2 = Vector3D.ZERO;
/*context.stations.get(0) // FIXME we need to access the station of the measurement
.getBaseFrame()
.getPVCoordinates(date2, frame)
.getPosition();*/
final Vector3D position2 = new Vector3D(measurements.get(idMeasure2).getObservedValue()[0], measurements.get(idMeasure2).getObservedValue()[1], measurements.get(idMeasure2).getObservedValue()[2]);
final Vector3D lineOfSight2 = position2.normalize();
// measurement data 3
final int idMeasure3 = 40;
final AbsoluteDate date3 = measurements.get(idMeasure3).getDate();
final Vector3D stapos3 = Vector3D.ZERO;
/*context.stations.get(0) // FIXME we need to access the station of the measurement
.getBaseFrame()
.getPVCoordinates(date3, frame)
.getPosition();*/
final Vector3D position3 = new Vector3D(measurements.get(idMeasure3).getObservedValue()[0], measurements.get(idMeasure3).getObservedValue()[1], measurements.get(idMeasure3).getObservedValue()[2]);
final double r3 = position3.getNorm();
final Vector3D lineOfSight3 = position3.normalize();
// instantiate the IOD method
final IodGooding iod = new IodGooding(frame, mu);
// the problem is very sensitive, and unless one can provide the exact
// initial range estimate, the estimate may be far off the truth...
final KeplerianOrbit orbit = iod.estimate(stapos1, stapos2, stapos3, lineOfSight1, date1, lineOfSight2, date2, lineOfSight3, date3, r1 * 1.0, r3 * 1.0);
Assert.assertEquals(orbit.getA(), context.initialOrbit.getA(), 1.0e-6 * context.initialOrbit.getA());
Assert.assertEquals(orbit.getE(), context.initialOrbit.getE(), 1.0e-6 * context.initialOrbit.getE());
Assert.assertEquals(orbit.getI(), context.initialOrbit.getI(), 1.0e-6 * context.initialOrbit.getI());
Assert.assertEquals(13127847.99808, iod.getRange1(), 1.0e-3);
Assert.assertEquals(13375711.51931, iod.getRange2(), 1.0e-3);
Assert.assertEquals(13950296.64852, iod.getRange3(), 1.0e-3);
}
Aggregations