use of org.orekit.propagation.conversion.NumericalPropagatorBuilder in project Orekit by CS-SI.
the class AngularAzElTest method testParameterDerivatives.
@Test
public void testParameterDerivatives() throws OrekitException {
Context context = EstimationTestUtils.geoStationnaryContext("regular-data:potential:tides");
final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, false, 1.0e-6, 60.0, 0.001);
// create perfect azimuth-elevation measurements
for (final GroundStation station : context.stations) {
station.getEastOffsetDriver().setSelected(true);
station.getNorthOffsetDriver().setSelected(true);
station.getZenithOffsetDriver().setSelected(true);
}
final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new AngularAzElMeasurementCreator(context), 0.25, 3.0, 600.0);
propagator.setSlaveMode();
for (final ObservedMeasurement<?> measurement : measurements) {
// parameter corresponding to station position offset
final GroundStation stationParameter = ((AngularAzEl) measurement).getStation();
// We intentionally propagate to a date which is close to the
// real spacecraft state but is *not* the accurate date, by
// compensating only part of the downlink delay. This is done
// in order to validate the partial derivatives with respect
// to velocity. If we had chosen the proper state date, the
// angular would have depended only on the current position but
// not on the current velocity.
final AbsoluteDate datemeas = measurement.getDate();
final SpacecraftState stateini = propagator.propagate(datemeas);
final Vector3D stationP = stationParameter.getOffsetToInertial(stateini.getFrame(), datemeas).transformPosition(Vector3D.ZERO);
final double meanDelay = AbstractMeasurement.signalTimeOfFlight(stateini.getPVCoordinates(), stationP, datemeas);
final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
final SpacecraftState state = propagator.propagate(date);
final ParameterDriver[] drivers = new ParameterDriver[] { stationParameter.getEastOffsetDriver(), stationParameter.getNorthOffsetDriver(), stationParameter.getZenithOffsetDriver() };
for (int i = 0; i < 3; ++i) {
final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
Assert.assertEquals(2, measurement.getDimension());
Assert.assertEquals(2, gradient.length);
for (final int k : new int[] { 0, 1 }) {
final ParameterFunction dMkdP = Differentiation.differentiate(new ParameterFunction() {
/**
* {@inheritDoc}
*/
@Override
public double value(final ParameterDriver parameterDriver) throws OrekitException {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[k];
}
}, drivers[i], 3, 50.0);
final double ref = dMkdP.value(drivers[i]);
if (ref > 1.e-12) {
Assert.assertEquals(ref, gradient[k], 3e-10 * FastMath.abs(ref));
}
}
}
}
}
use of org.orekit.propagation.conversion.NumericalPropagatorBuilder in project Orekit by CS-SI.
the class GroundStationTest method testEstimateEOP.
@Test
public void testEstimateEOP() throws OrekitException {
Context linearEOPContext = EstimationTestUtils.eccentricContext("linear-EOP:regular-data/de431-ephemerides:potential:tides");
final AbsoluteDate refDate = new AbsoluteDate(2000, 2, 24, linearEOPContext.utc);
final double dut10 = 0.3079738;
final double lod = 0.0011000;
final double xp0 = 68450.0e-6;
final double xpDot = -50.0e-6;
final double yp0 = 60.0e-6;
final double ypDot = 2.0e-6;
for (double dt = -2 * Constants.JULIAN_DAY; dt < 2 * Constants.JULIAN_DAY; dt += 300.0) {
AbsoluteDate date = refDate.shiftedBy(dt);
Assert.assertEquals(dut10 - dt * lod / Constants.JULIAN_DAY, linearEOPContext.ut1.getEOPHistory().getUT1MinusUTC(date), 1.0e-15);
Assert.assertEquals(lod, linearEOPContext.ut1.getEOPHistory().getLOD(date), 1.0e-15);
Assert.assertEquals((xp0 + xpDot * dt / Constants.JULIAN_DAY) * Constants.ARC_SECONDS_TO_RADIANS, linearEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getXp(), 1.0e-15);
Assert.assertEquals((yp0 + ypDot * dt / Constants.JULIAN_DAY) * Constants.ARC_SECONDS_TO_RADIANS, linearEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getYp(), 1.0e-15);
}
final NumericalPropagatorBuilder linearPropagatorBuilder = linearEOPContext.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
// create perfect range measurements
final Propagator propagator = EstimationTestUtils.createPropagator(linearEOPContext.initialOrbit, linearPropagatorBuilder);
final List<ObservedMeasurement<?>> linearMeasurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(linearEOPContext), 1.0, 5.0, 60.0);
Utils.clearFactories();
Context zeroEOPContext = EstimationTestUtils.eccentricContext("zero-EOP:regular-data/de431-ephemerides:potential:potential:tides");
for (double dt = -2 * Constants.JULIAN_DAY; dt < 2 * Constants.JULIAN_DAY; dt += 300.0) {
AbsoluteDate date = refDate.shiftedBy(dt);
Assert.assertEquals(0.0, zeroEOPContext.ut1.getEOPHistory().getUT1MinusUTC(date), 1.0e-15);
Assert.assertEquals(0.0, zeroEOPContext.ut1.getEOPHistory().getLOD(date), 1.0e-15);
Assert.assertEquals(0.0, zeroEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getXp(), 1.0e-15);
Assert.assertEquals(0.0, zeroEOPContext.ut1.getEOPHistory().getPoleCorrection(date).getYp(), 1.0e-15);
}
// create orbit estimator
final NumericalPropagatorBuilder zeroPropagatorBuilder = linearEOPContext.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), zeroPropagatorBuilder);
for (final ObservedMeasurement<?> linearMeasurement : linearMeasurements) {
Range linearRange = (Range) linearMeasurement;
for (final GroundStation station : zeroEOPContext.stations) {
if (station.getBaseFrame().getName().equals(linearRange.getStation().getBaseFrame().getName())) {
Range zeroRange = new Range(station, linearRange.getDate(), linearRange.getObservedValue()[0], linearRange.getTheoreticalStandardDeviation()[0], linearRange.getBaseWeight()[0]);
estimator.addMeasurement(zeroRange);
}
}
}
estimator.setParametersConvergenceThreshold(1.0e-3);
estimator.setMaxIterations(100);
estimator.setMaxEvaluations(200);
// we want to estimate pole and prime meridian
GroundStation station = zeroEOPContext.stations.get(0);
station.getPrimeMeridianOffsetDriver().setReferenceDate(refDate);
station.getPrimeMeridianOffsetDriver().setSelected(true);
station.getPrimeMeridianDriftDriver().setSelected(true);
station.getPolarOffsetXDriver().setReferenceDate(refDate);
station.getPolarOffsetXDriver().setSelected(true);
station.getPolarDriftXDriver().setSelected(true);
station.getPolarOffsetYDriver().setReferenceDate(refDate);
station.getPolarOffsetYDriver().setSelected(true);
station.getPolarDriftYDriver().setSelected(true);
// just for the fun and to speed up test, we will use orbit determination, *without* estimating orbit
for (final ParameterDriver driver : zeroPropagatorBuilder.getOrbitalParametersDrivers().getDrivers()) {
driver.setSelected(false);
}
estimator.estimate();
final double computedDut1 = station.getPrimeMeridianOffsetDriver().getValue() / EstimatedEarthFrameProvider.EARTH_ANGULAR_VELOCITY;
final double computedLOD = station.getPrimeMeridianDriftDriver().getValue() * (-Constants.JULIAN_DAY / EstimatedEarthFrameProvider.EARTH_ANGULAR_VELOCITY);
final double computedXp = station.getPolarOffsetXDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS;
final double computedXpDot = station.getPolarDriftXDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS * Constants.JULIAN_DAY;
final double computedYp = station.getPolarOffsetYDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS;
final double computedYpDot = station.getPolarDriftYDriver().getValue() / Constants.ARC_SECONDS_TO_RADIANS * Constants.JULIAN_DAY;
Assert.assertEquals(dut10, computedDut1, 4.3e-10);
Assert.assertEquals(lod, computedLOD, 4.9e-10);
Assert.assertEquals(xp0, computedXp, 5.6e-9);
Assert.assertEquals(xpDot, computedXpDot, 7.2e-9);
Assert.assertEquals(yp0, computedYp, 1.1e-9);
Assert.assertEquals(ypDot, computedYpDot, 2.8e-11);
// thresholds to use if orbit is estimated
// (i.e. when commenting out the loop above that sets orbital parameters drivers to "not selected")
// Assert.assertEquals(dut10, computedDut1, 6.6e-3);
// Assert.assertEquals(lod, computedLOD, 1.1e-9);
// Assert.assertEquals(xp0, computedXp, 3.3e-8);
// Assert.assertEquals(xpDot, computedXpDot, 2.2e-8);
// Assert.assertEquals(yp0, computedYp, 3.3e-8);
// Assert.assertEquals(ypDot, computedYpDot, 3.8e-8);
}
use of org.orekit.propagation.conversion.NumericalPropagatorBuilder in project Orekit by CS-SI.
the class TurnAroundRangeAnalyticTest method genericTestStateDerivatives.
/**
* Generic test function for derivatives with respect to state
* @param isModifier Use of atmospheric modifiers
* @param isFiniteDifferences Finite differences reference calculation if true, TurnAroundRange class otherwise
* @param printResults Print the results ?
* @throws OrekitException
*/
void genericTestStateDerivatives(final boolean isModifier, final boolean isFiniteDifferences, final boolean printResults, final double refErrorsPMedian, final double refErrorsPMean, final double refErrorsPMax, final double refErrorsVMedian, final double refErrorsVMean, final double refErrorsVMax) throws OrekitException {
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
// Context context = EstimationTestUtils.geoStationnaryContext();
final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
// create perfect range2 measurements
final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new TurnAroundRangeMeasurementCreator(context), 1.0, 3.0, 300.0);
propagator.setSlaveMode();
double[] errorsP = new double[3 * measurements.size()];
double[] errorsV = new double[3 * measurements.size()];
int indexP = 0;
int indexV = 0;
// Print the results ? Header
if (printResults) {
System.out.format(Locale.US, "%-15s %-15s %-23s %-23s " + "%10s %10s %10s " + "%10s %10s %10s " + "%10s %10s %10s " + "%10s %10s %10s%n", "Master Station", "Slave Station", "Measurement Date", "State Date", "ΔdPx", "ΔdPy", "ΔdPz", "ΔdVx", "ΔdVy", "ΔdVz", "rel ΔdPx", "rel ΔdPy", "rel ΔdPz", "rel ΔdVx", "rel ΔdVy", "rel ΔdVz");
}
// Loop on the measurements
for (final ObservedMeasurement<?> measurement : measurements) {
// Add modifiers if test implies it
final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
if (isModifier) {
((TurnAroundRange) measurement).addModifier(modifier);
}
// We intentionally propagate to a date which is close to the
// real spacecraft state but is *not* the accurate date, by
// compensating only part of the downlink delay. This is done
// in order to validate the partial derivatives with respect
// to velocity. If we had chosen the proper state date, the
// range would have depended only on the current position but
// not on the current velocity.
final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
final SpacecraftState state = propagator.propagate(date);
final EstimatedMeasurement<TurnAroundRange> TAR = new TurnAroundRangeAnalytic((TurnAroundRange) measurement).theoreticalEvaluationAnalytic(0, 0, propagator.getInitialState(), state);
if (isModifier) {
modifier.modify(TAR);
}
final double[][] jacobian = TAR.getStateDerivatives(0);
// Jacobian reference value
final double[][] jacobianRef;
if (isFiniteDifferences) {
// Compute a reference value using finite differences
jacobianRef = Differentiation.differentiate(new StateFunction() {
public double[] value(final SpacecraftState state) throws OrekitException {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue();
}
}, measurement.getDimension(), propagator.getAttitudeProvider(), OrbitType.CARTESIAN, PositionAngle.TRUE, 2.0, 3).value(state);
} else {
// Compute a reference value using TurnAroundRange class function
jacobianRef = ((TurnAroundRange) measurement).theoreticalEvaluation(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
}
// //Test: Test point by point with the debugger
// if (!isFiniteDifferences && !isModifier) {
// final EstimatedMeasurement<TurnAroundRange> test =
// new TurnAroundRangeAnalytic((TurnAroundRange)measurement).theoreticalEvaluationValidation(0, 0, state);
// }
// //Test
Assert.assertEquals(jacobianRef.length, jacobian.length);
Assert.assertEquals(jacobianRef[0].length, jacobian[0].length);
double[][] dJacobian = new double[jacobian.length][jacobian[0].length];
double[][] dJacobianRelative = new double[jacobian.length][jacobian[0].length];
for (int i = 0; i < jacobian.length; ++i) {
for (int j = 0; j < jacobian[i].length; ++j) {
dJacobian[i][j] = jacobian[i][j] - jacobianRef[i][j];
dJacobianRelative[i][j] = FastMath.abs(dJacobian[i][j] / jacobianRef[i][j]);
if (j < 3) {
errorsP[indexP++] = dJacobianRelative[i][j];
} else {
errorsV[indexV++] = dJacobianRelative[i][j];
}
}
}
// Print results on the console ? Print the Jacobian
if (printResults) {
String masterStationName = ((TurnAroundRange) measurement).getMasterStation().getBaseFrame().getName();
String slaveStationName = ((TurnAroundRange) measurement).getSlaveStation().getBaseFrame().getName();
System.out.format(Locale.US, "%-15s %-15s %-23s %-23s " + "%10.3e %10.3e %10.3e " + "%10.3e %10.3e %10.3e " + "%10.3e %10.3e %10.3e " + "%10.3e %10.3e %10.3e%n", masterStationName, slaveStationName, measurement.getDate(), date, dJacobian[0][0], dJacobian[0][1], dJacobian[0][2], dJacobian[0][3], dJacobian[0][4], dJacobian[0][5], dJacobianRelative[0][0], dJacobianRelative[0][1], dJacobianRelative[0][2], dJacobianRelative[0][3], dJacobianRelative[0][4], dJacobianRelative[0][5]);
}
}
// End loop on the measurements
// Compute some statistics
final double errorsPMedian = new Median().evaluate(errorsP);
final double errorsPMean = new Mean().evaluate(errorsP);
final double errorsPMax = new Max().evaluate(errorsP);
final double errorsVMedian = new Median().evaluate(errorsV);
final double errorsVMean = new Mean().evaluate(errorsV);
final double errorsVMax = new Max().evaluate(errorsV);
// Print the results on console ? Final results
if (printResults) {
System.out.println();
System.out.format(Locale.US, "Relative errors dR/dP -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", errorsPMedian, errorsPMean, errorsPMax);
System.out.format(Locale.US, "Relative errors dR/dV -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", errorsVMedian, errorsVMean, errorsVMax);
}
// Assert the results / max values depend on the test
Assert.assertEquals(0.0, errorsPMedian, refErrorsPMedian);
Assert.assertEquals(0.0, errorsPMean, refErrorsPMean);
Assert.assertEquals(0.0, errorsPMax, refErrorsPMax);
Assert.assertEquals(0.0, errorsVMedian, refErrorsVMedian);
Assert.assertEquals(0.0, errorsVMean, refErrorsVMean);
Assert.assertEquals(0.0, errorsVMax, refErrorsVMax);
}
use of org.orekit.propagation.conversion.NumericalPropagatorBuilder in project Orekit by CS-SI.
the class TurnAroundRangeTest method genericTestParameterDerivatives.
void genericTestParameterDerivatives(final boolean isModifier, final boolean printResults, final double refErrorQMMedian, final double refErrorQMMean, final double refErrorQMMax, final double refErrorQSMedian, final double refErrorQSMean, final double refErrorQSMax) throws OrekitException {
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
// Create perfect TAR measurements
for (Map.Entry<GroundStation, GroundStation> entry : context.TARstations.entrySet()) {
final GroundStation masterStation = entry.getKey();
final GroundStation slaveStation = entry.getValue();
masterStation.getEastOffsetDriver().setSelected(true);
masterStation.getNorthOffsetDriver().setSelected(true);
masterStation.getZenithOffsetDriver().setSelected(true);
slaveStation.getEastOffsetDriver().setSelected(true);
slaveStation.getNorthOffsetDriver().setSelected(true);
slaveStation.getZenithOffsetDriver().setSelected(true);
}
final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new TurnAroundRangeMeasurementCreator(context), 1.0, 3.0, 300.0);
propagator.setSlaveMode();
// Print results on console ? Header
if (printResults) {
System.out.format(Locale.US, "%-15s %-15s %-23s %-23s " + "%10s %10s %10s " + "%10s %10s %10s " + "%10s %10s %10s " + "%10s %10s %10s%n", "Master Station", "Slave Station", "Measurement Date", "State Date", "ΔdQMx", "rel ΔdQMx", "ΔdQMy", "rel ΔdQMy", "ΔdQMz", "rel ΔdQMz", "ΔdQSx", "rel ΔdQSx", "ΔdQSy", "rel ΔdQSy", "ΔdQSz", "rel ΔdQSz");
}
// List to store the results for master and slave station
final List<Double> relErrorQMList = new ArrayList<Double>();
final List<Double> relErrorQSList = new ArrayList<Double>();
// Loop on the measurements
for (final ObservedMeasurement<?> measurement : measurements) {
// Add modifiers if test implies it
final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
if (isModifier) {
((TurnAroundRange) measurement).addModifier(modifier);
}
// parameter corresponding to station position offset
final GroundStation masterStationParameter = ((TurnAroundRange) measurement).getMasterStation();
final GroundStation slaveStationParameter = ((TurnAroundRange) measurement).getSlaveStation();
// We intentionally propagate to a date which is close to the
// real spacecraft state but is *not* the accurate date, by
// compensating only part of the downlink delay. This is done
// in order to validate the partial derivatives with respect
// to velocity. If we had chosen the proper state date, the
// range would have depended only on the current position but
// not on the current velocity.
final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
final SpacecraftState state = propagator.propagate(date);
final ParameterDriver[] drivers = new ParameterDriver[] { masterStationParameter.getEastOffsetDriver(), masterStationParameter.getNorthOffsetDriver(), masterStationParameter.getZenithOffsetDriver(), slaveStationParameter.getEastOffsetDriver(), slaveStationParameter.getNorthOffsetDriver(), slaveStationParameter.getZenithOffsetDriver() };
// Print results on console ? Stations' names
if (printResults) {
String masterStationName = masterStationParameter.getBaseFrame().getName();
String slaveStationName = slaveStationParameter.getBaseFrame().getName();
System.out.format(Locale.US, "%-15s %-15s %-23s %-23s ", masterStationName, slaveStationName, measurement.getDate(), date);
}
// Loop on the parameters
for (int i = 0; i < 6; ++i) {
final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
Assert.assertEquals(1, measurement.getDimension());
Assert.assertEquals(1, gradient.length);
// Compute a reference value using finite differences
final ParameterFunction dMkdP = Differentiation.differentiate(new ParameterFunction() {
/**
* {@inheritDoc}
*/
@Override
public double value(final ParameterDriver parameterDriver) throws OrekitException {
return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
}
}, drivers[i], 3, 20.0);
final double ref = dMkdP.value(drivers[i]);
// Deltas
double dGradient = gradient[0] - ref;
double dGradientRelative = FastMath.abs(dGradient / ref);
// Print results on console ? Gradient difference
if (printResults) {
System.out.format(Locale.US, "%10.3e %10.3e ", dGradient, dGradientRelative);
}
// Add relative error to the list
if (i < 3) {
relErrorQMList.add(dGradientRelative);
} else {
relErrorQSList.add(dGradientRelative);
}
}
// End for loop on the parameters
if (printResults) {
System.out.format(Locale.US, "%n");
}
}
// End for loop on the measurements
// Convert error list to double[]
final double[] relErrorQM = relErrorQMList.stream().mapToDouble(Double::doubleValue).toArray();
final double[] relErrorQS = relErrorQSList.stream().mapToDouble(Double::doubleValue).toArray();
// Compute statistics
final double relErrorsQMMedian = new Median().evaluate(relErrorQM);
final double relErrorsQMMean = new Mean().evaluate(relErrorQM);
final double relErrorsQMMax = new Max().evaluate(relErrorQM);
final double relErrorsQSMedian = new Median().evaluate(relErrorQS);
final double relErrorsQSMean = new Mean().evaluate(relErrorQS);
final double relErrorsQSMax = new Max().evaluate(relErrorQS);
// Print the results on console ?
if (printResults) {
System.out.println();
System.out.format(Locale.US, "Relative errors dR/dQ master station -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", relErrorsQMMedian, relErrorsQMMean, relErrorsQMMax);
System.out.format(Locale.US, "Relative errors dR/dQ slave station -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", relErrorsQSMedian, relErrorsQSMean, relErrorsQSMax);
}
// Check values
Assert.assertEquals(0.0, relErrorsQMMedian, refErrorQMMedian);
Assert.assertEquals(0.0, relErrorsQMMean, refErrorQMMean);
Assert.assertEquals(0.0, relErrorsQMMax, refErrorQMMax);
Assert.assertEquals(0.0, relErrorsQSMedian, refErrorQSMedian);
Assert.assertEquals(0.0, relErrorsQSMean, refErrorQSMean);
Assert.assertEquals(0.0, relErrorsQSMax, refErrorQSMax);
}
use of org.orekit.propagation.conversion.NumericalPropagatorBuilder in project Orekit by CS-SI.
the class BiasTest method testEstimateBias.
@SuppressWarnings("unchecked")
@Test
public void testEstimateBias() throws OrekitException {
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
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 RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
// create range biases: one bias for each station
final RandomGenerator random = new Well19937a(0x0c4b69da5d64b35al);
final Bias<?>[] stationsRangeBiases = new Bias<?>[context.stations.size()];
final double[] realStationsBiases = new double[context.stations.size()];
for (int i = 0; i < context.stations.size(); ++i) {
final TopocentricFrame base = context.stations.get(i).getBaseFrame();
stationsRangeBiases[i] = new Bias<Range>(new String[] { base.getName() + " range bias" }, new double[] { 0.0 }, new double[] { 1.0 }, new double[] { Double.NEGATIVE_INFINITY }, new double[] { Double.POSITIVE_INFINITY });
realStationsBiases[i] = 2 * random.nextDouble() - 1;
}
// create orbit estimator
final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
// add the measurements, with both spacecraft and stations biases
for (final ObservedMeasurement<?> measurement : measurements) {
final Range range = (Range) measurement;
for (int i = 0; i < context.stations.size(); ++i) {
if (range.getStation() == context.stations.get(i)) {
double biasedRange = range.getObservedValue()[0] + realStationsBiases[i];
final Range m = new Range(range.getStation(), range.getDate(), biasedRange, range.getTheoreticalStandardDeviation()[0], range.getBaseWeight()[0]);
m.addModifier((Bias<Range>) stationsRangeBiases[i]);
estimator.addMeasurement(m);
}
}
}
estimator.setParametersConvergenceThreshold(1.0e-3);
estimator.setMaxIterations(10);
estimator.setMaxEvaluations(20);
// we want to estimate the biases
for (Bias<?> bias : stationsRangeBiases) {
for (final ParameterDriver driver : bias.getParametersDrivers()) {
driver.setSelected(true);
}
}
EstimationTestUtils.checkFit(context, estimator, 2, 3, 0.0, 7.2e-7, 0.0, 2.1e-6, 0.0, 3.7e-7, 0.0, 1.7e-10);
for (int i = 0; i < stationsRangeBiases.length; ++i) {
Assert.assertEquals(realStationsBiases[i], stationsRangeBiases[i].getParametersDrivers().get(0).getValue(), 3.3e-6);
}
}
Aggregations