use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.
the class KalmanEstimatorTest method testKeplerianRange.
/**
* Perfect range measurements with a biased start
* Keplerian formalism
* @throws OrekitException
*/
@Test
public void testKeplerianRange() throws OrekitException {
// Create context
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
// Create initial orbit and propagator builder
final OrbitType orbitType = OrbitType.KEPLERIAN;
final PositionAngle positionAngle = PositionAngle.TRUE;
final boolean perfectStart = true;
final double minStep = 1.e-6;
final double maxStep = 60.;
final double dP = 1.;
final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(orbitType, positionAngle, perfectStart, minStep, maxStep, dP);
// Create perfect range measurements
final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(context), 1.0, 4.0, 60.0);
// Reference propagator for estimation performances
final NumericalPropagator referencePropagator = propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
// Reference position/velocity at last measurement date
final Orbit refOrbit = referencePropagator.propagate(measurements.get(measurements.size() - 1).getDate()).getOrbit();
// Change semi-major axis of 1.2m as in the batch test
ParameterDriver aDriver = propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(0);
aDriver.setValue(aDriver.getValue() + 1.2);
aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
// Cartesian covariance matrix initialization
// 100m on position / 1e-2m/s on velocity
final RealMatrix cartesianP = MatrixUtils.createRealDiagonalMatrix(new double[] { 100., 100., 100., 1e-2, 1e-2, 1e-2 });
// Jacobian of the orbital parameters w/r to Cartesian
final Orbit initialOrbit = orbitType.convertType(context.initialOrbit);
final double[][] dYdC = new double[6][6];
initialOrbit.getJacobianWrtCartesian(PositionAngle.TRUE, dYdC);
final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
// Keplerian initial covariance matrix
final RealMatrix initialP = Jac.multiply(cartesianP.multiply(Jac.transpose()));
// Process noise matrix is set to 0 here
RealMatrix Q = MatrixUtils.createRealMatrix(6, 6);
// Build the Kalman filter
final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
kalmanBuilder.builder(propagatorBuilder);
kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
kalmanBuilder.initialCovarianceMatrix(initialP);
kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(Q));
final KalmanEstimator kalman = kalmanBuilder.build();
// Filter the measurements and check the results
final double expectedDeltaPos = 0.;
final double posEps = 1.77e-4;
final double expectedDeltaVel = 0.;
final double velEps = 7.93e-8;
final double[] expectedSigmasPos = { 0.742488, 0.281910, 0.563217 };
final double sigmaPosEps = 1e-6;
final double[] expectedSigmasVel = { 2.206622e-4, 1.306669e-4, 1.293996e-4 };
final double sigmaVelEps = 1e-10;
EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedSigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}
use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.
the class KalmanOrbitDeterminationTest method displayParametersChanges.
/**
* Display parameters changes.
* @param stream output stream
* @param header header message
* @param sort if true, parameters will be sorted lexicographically
* @param parameters parameters list
*/
private void displayParametersChanges(final PrintStream out, final String header, final boolean sort, final int length, final ParameterDriversList parameters) {
List<ParameterDriver> list = new ArrayList<ParameterDriver>(parameters.getDrivers());
if (sort) {
// sort the parameters lexicographically
Collections.sort(list, new Comparator<ParameterDriver>() {
/**
* {@inheritDoc}
*/
@Override
public int compare(final ParameterDriver pd1, final ParameterDriver pd2) {
return pd1.getName().compareTo(pd2.getName());
}
});
}
out.println(header);
int index = 0;
for (final ParameterDriver parameter : list) {
if (parameter.isSelected()) {
final double factor;
if (parameter.getName().endsWith("/az bias") || parameter.getName().endsWith("/el bias")) {
factor = FastMath.toDegrees(1.0);
} else {
factor = 1.0;
}
final double initial = parameter.getReferenceValue();
final double value = parameter.getValue();
out.format(Locale.US, " %2d %s", ++index, parameter.getName());
for (int i = parameter.getName().length(); i < length; ++i) {
out.format(Locale.US, " ");
}
out.format(Locale.US, " %+.12f (final value: % .12f)%n", factor * (value - initial), factor * value);
}
}
}
use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.
the class KalmanOrbitDeterminationTest method displayFinalCovariances.
/**
* Display covariances and sigmas as predicted by a Kalman filter at date t.
* @throws OrekitException
*/
private void displayFinalCovariances(final PrintStream logStream, final KalmanEstimator kalman) throws OrekitException {
// // Get kalman estimated propagator
// final NumericalPropagator kalmanProp = kalman.getProcessModel().getEstimatedPropagator();
//
// // Link the partial derivatives to this propagator
// final String equationName = "kalman-derivatives";
// PartialDerivativesEquations kalmanDerivatives = new PartialDerivativesEquations(equationName, kalmanProp);
//
// // Initialize the derivatives
// final SpacecraftState rawState = kalmanProp.getInitialState();
// final SpacecraftState stateWithDerivatives =
// kalmanDerivatives.setInitialJacobians(rawState);
// kalmanProp.resetInitialState(stateWithDerivatives);
//
// // Propagate to target date
// final SpacecraftState kalmanState = kalmanProp.propagate(targetDate);
//
// // Compute STM
// RealMatrix STM = kalman.getProcessModel().getErrorStateTransitionMatrix(kalmanState, kalmanDerivatives);
//
// // Compute covariance matrix
// RealMatrix P = kalman.getProcessModel().unNormalizeCovarianceMatrix(kalman.predictCovariance(STM,
// kalman.getProcessModel().getProcessNoiseMatrix()));
final RealMatrix P = kalman.getPhysicalEstimatedCovarianceMatrix();
final String[] paramNames = new String[P.getRowDimension()];
int index = 0;
int paramSize = 0;
for (final ParameterDriver driver : kalman.getOrbitalParametersDrivers(true).getDrivers()) {
paramNames[index++] = driver.getName();
paramSize = FastMath.max(paramSize, driver.getName().length());
}
for (final ParameterDriver driver : kalman.getPropagationParametersDrivers(true).getDrivers()) {
paramNames[index++] = driver.getName();
paramSize = FastMath.max(paramSize, driver.getName().length());
}
for (final ParameterDriver driver : kalman.getEstimatedMeasurementsParameters().getDrivers()) {
paramNames[index++] = driver.getName();
paramSize = FastMath.max(paramSize, driver.getName().length());
}
if (paramSize < 20) {
paramSize = 20;
}
// Header
logStream.format("\n%s\n", "Kalman Final Covariances:");
// logStream.format(Locale.US, "\tDate: %-23s UTC\n",
// targetDate.toString(TimeScalesFactory.getUTC()));
logStream.format(Locale.US, "\tDate: %-23s UTC\n", kalman.getCurrentDate().toString(TimeScalesFactory.getUTC()));
// Covariances
String strFormat = String.format("%%%2ds ", paramSize);
logStream.format(strFormat, "Covariances:");
for (int i = 0; i < P.getRowDimension(); i++) {
logStream.format(Locale.US, strFormat, paramNames[i]);
}
logStream.println("");
String numFormat = String.format("%%%2d.6f ", paramSize);
for (int i = 0; i < P.getRowDimension(); i++) {
logStream.format(Locale.US, strFormat, paramNames[i]);
for (int j = 0; j <= i; j++) {
logStream.format(Locale.US, numFormat, P.getEntry(i, j));
}
logStream.println("");
}
// Correlation coeff
final double[] sigmas = new double[P.getRowDimension()];
for (int i = 0; i < P.getRowDimension(); i++) {
sigmas[i] = FastMath.sqrt(P.getEntry(i, i));
}
logStream.format("\n" + strFormat, "Corr coef:");
for (int i = 0; i < P.getRowDimension(); i++) {
logStream.format(Locale.US, strFormat, paramNames[i]);
}
logStream.println("");
for (int i = 0; i < P.getRowDimension(); i++) {
logStream.format(Locale.US, strFormat, paramNames[i]);
for (int j = 0; j <= i; j++) {
logStream.format(Locale.US, numFormat, P.getEntry(i, j) / (sigmas[i] * sigmas[j]));
}
logStream.println("");
}
// Sigmas
logStream.format("\n" + strFormat + "\n", "Sigmas: ");
for (int i = 0; i < P.getRowDimension(); i++) {
logStream.format(Locale.US, strFormat + numFormat + "\n", paramNames[i], sigmas[i]);
}
logStream.println("");
}
use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.
the class KalmanOrbitDeterminationTest method createPropagatorBuilder.
/**
* Create a propagator builder from input parameters
* @param parser input file parser
* @param conventions IERS conventions to use
* @param gravityField gravity field
* @param body central body
* @param orbit first orbit estimate
* @return propagator builder
* @throws NoSuchElementException if input parameters are missing
* @throws OrekitException if body frame cannot be created
*/
private NumericalPropagatorBuilder createPropagatorBuilder(final KeyValueFileParser<ParameterKey> parser, final IERSConventions conventions, final NormalizedSphericalHarmonicsProvider gravityField, final OneAxisEllipsoid body, final Orbit orbit) throws NoSuchElementException, OrekitException {
final double minStep;
if (!parser.containsKey(ParameterKey.PROPAGATOR_MIN_STEP)) {
minStep = 0.001;
} else {
minStep = parser.getDouble(ParameterKey.PROPAGATOR_MIN_STEP);
}
final double maxStep;
if (!parser.containsKey(ParameterKey.PROPAGATOR_MAX_STEP)) {
maxStep = 300;
} else {
maxStep = parser.getDouble(ParameterKey.PROPAGATOR_MAX_STEP);
}
final double dP;
if (!parser.containsKey(ParameterKey.PROPAGATOR_POSITION_ERROR)) {
dP = 10.0;
} else {
dP = parser.getDouble(ParameterKey.PROPAGATOR_POSITION_ERROR);
}
final double positionScale;
if (!parser.containsKey(ParameterKey.ESTIMATOR_ORBITAL_PARAMETERS_POSITION_SCALE)) {
positionScale = dP;
} else {
positionScale = parser.getDouble(ParameterKey.ESTIMATOR_ORBITAL_PARAMETERS_POSITION_SCALE);
}
final NumericalPropagatorBuilder propagatorBuilder = new NumericalPropagatorBuilder(orbit, new DormandPrince853IntegratorBuilder(minStep, maxStep, dP), PositionAngle.MEAN, positionScale);
// initial mass
final double mass;
if (!parser.containsKey(ParameterKey.MASS)) {
mass = 1000.0;
} else {
mass = parser.getDouble(ParameterKey.MASS);
}
propagatorBuilder.setMass(mass);
// gravity field force model
propagatorBuilder.addForceModel(new HolmesFeatherstoneAttractionModel(body.getBodyFrame(), gravityField));
// ocean tides force model
if (parser.containsKey(ParameterKey.OCEAN_TIDES_DEGREE) && parser.containsKey(ParameterKey.OCEAN_TIDES_ORDER)) {
final int degree = parser.getInt(ParameterKey.OCEAN_TIDES_DEGREE);
final int order = parser.getInt(ParameterKey.OCEAN_TIDES_ORDER);
if (degree > 0 && order > 0) {
propagatorBuilder.addForceModel(new OceanTides(body.getBodyFrame(), gravityField.getAe(), gravityField.getMu(), degree, order, conventions, TimeScalesFactory.getUT1(conventions, true)));
}
}
// solid tides force model
List<CelestialBody> solidTidesBodies = new ArrayList<CelestialBody>();
if (parser.containsKey(ParameterKey.SOLID_TIDES_SUN) && parser.getBoolean(ParameterKey.SOLID_TIDES_SUN)) {
solidTidesBodies.add(CelestialBodyFactory.getSun());
}
if (parser.containsKey(ParameterKey.SOLID_TIDES_MOON) && parser.getBoolean(ParameterKey.SOLID_TIDES_MOON)) {
solidTidesBodies.add(CelestialBodyFactory.getMoon());
}
if (!solidTidesBodies.isEmpty()) {
propagatorBuilder.addForceModel(new SolidTides(body.getBodyFrame(), gravityField.getAe(), gravityField.getMu(), gravityField.getTideSystem(), conventions, TimeScalesFactory.getUT1(conventions, true), solidTidesBodies.toArray(new CelestialBody[solidTidesBodies.size()])));
}
// third body attraction
if (parser.containsKey(ParameterKey.THIRD_BODY_SUN) && parser.getBoolean(ParameterKey.THIRD_BODY_SUN)) {
propagatorBuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
}
if (parser.containsKey(ParameterKey.THIRD_BODY_MOON) && parser.getBoolean(ParameterKey.THIRD_BODY_MOON)) {
propagatorBuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
}
// drag
if (parser.containsKey(ParameterKey.DRAG) && parser.getBoolean(ParameterKey.DRAG)) {
final double cd = parser.getDouble(ParameterKey.DRAG_CD);
final double area = parser.getDouble(ParameterKey.DRAG_AREA);
final boolean cdEstimated = parser.getBoolean(ParameterKey.DRAG_CD_ESTIMATED);
MarshallSolarActivityFutureEstimation msafe = new MarshallSolarActivityFutureEstimation("(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}F10\\.(?:txt|TXT)", MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
DataProvidersManager manager = DataProvidersManager.getInstance();
manager.feed(msafe.getSupportedNames(), msafe);
Atmosphere atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), body);
propagatorBuilder.addForceModel(new DragForce(atmosphere, new IsotropicDrag(area, cd)));
if (cdEstimated) {
for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
if (driver.getName().equals(DragSensitive.DRAG_COEFFICIENT)) {
driver.setSelected(true);
}
}
}
}
// solar radiation pressure
if (parser.containsKey(ParameterKey.SOLAR_RADIATION_PRESSURE) && parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE)) {
final double cr = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_CR);
final double area = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_AREA);
final boolean cREstimated = parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE_CR_ESTIMATED);
propagatorBuilder.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), body.getEquatorialRadius(), new IsotropicRadiationSingleCoefficient(area, cr)));
if (cREstimated) {
for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
driver.setSelected(true);
}
}
}
}
// post-Newtonian correction force due to general relativity
if (parser.containsKey(ParameterKey.GENERAL_RELATIVITY) && parser.getBoolean(ParameterKey.GENERAL_RELATIVITY)) {
propagatorBuilder.addForceModel(new Relativity(gravityField.getMu()));
}
// extra polynomial accelerations
if (parser.containsKey(ParameterKey.POLYNOMIAL_ACCELERATION_NAME)) {
final String[] names = parser.getStringArray(ParameterKey.POLYNOMIAL_ACCELERATION_NAME);
final Vector3D[] directions = parser.getVectorArray(ParameterKey.POLYNOMIAL_ACCELERATION_DIRECTION_X, ParameterKey.POLYNOMIAL_ACCELERATION_DIRECTION_Y, ParameterKey.POLYNOMIAL_ACCELERATION_DIRECTION_Z);
final List<String>[] coefficients = parser.getStringsListArray(ParameterKey.POLYNOMIAL_ACCELERATION_COEFFICIENTS, ',');
final boolean[] estimated = parser.getBooleanArray(ParameterKey.POLYNOMIAL_ACCELERATION_ESTIMATED);
for (int i = 0; i < names.length; ++i) {
final PolynomialParametricAcceleration ppa = new PolynomialParametricAcceleration(directions[i], true, names[i], null, coefficients[i].size() - 1);
for (int k = 0; k < coefficients[i].size(); ++k) {
final ParameterDriver driver = ppa.getParameterDriver(names[i] + "[" + k + "]");
driver.setValue(Double.parseDouble(coefficients[i].get(k)));
driver.setSelected(estimated[i]);
}
propagatorBuilder.addForceModel(ppa);
}
}
return propagatorBuilder;
}
use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.
the class KalmanOrbitDeterminationTest method run.
/**
* Function running the Kalman filter estimation.
* @param input Input configuration file
* @param orbitType Orbit type to use (calculation and display)
* @param print Choose whether the results are printed on console or not
* @param cartesianOrbitalP Orbital part of the initial covariance matrix in Cartesian formalism
* @param cartesianOrbitalQ Orbital part of the process noise matrix in Cartesian formalism
* @param propagationP Propagation part of the initial covariance matrix
* @param propagationQ Propagation part of the process noise matrix
* @param measurementP Measurement part of the initial covariance matrix
* @param measurementQ Measurement part of the process noise matrix
*/
private ResultKalman run(final File input, final OrbitType orbitType, final boolean print, final RealMatrix cartesianOrbitalP, final RealMatrix cartesianOrbitalQ, final RealMatrix propagationP, final RealMatrix propagationQ, final RealMatrix measurementP, final RealMatrix measurementQ) throws IOException, IllegalArgumentException, OrekitException, ParseException {
// Read input parameters
KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
parser.parseInput(input.getAbsolutePath(), new FileInputStream(input));
// Log files
final RangeLog rangeLog = new RangeLog();
final RangeRateLog rangeRateLog = new RangeRateLog();
final AzimuthLog azimuthLog = new AzimuthLog();
final ElevationLog elevationLog = new ElevationLog();
final PositionLog positionLog = new PositionLog();
final VelocityLog velocityLog = new VelocityLog();
// Gravity field
GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-5c.gfc", true));
final NormalizedSphericalHarmonicsProvider gravityField = createGravityField(parser);
// Orbit initial guess
Orbit initialGuess = createOrbit(parser, gravityField.getMu());
// Convert to desired orbit type
initialGuess = orbitType.convertType(initialGuess);
// IERS conventions
final IERSConventions conventions;
if (!parser.containsKey(ParameterKey.IERS_CONVENTIONS)) {
conventions = IERSConventions.IERS_2010;
} else {
conventions = IERSConventions.valueOf("IERS_" + parser.getInt(ParameterKey.IERS_CONVENTIONS));
}
// Central body
final OneAxisEllipsoid body = createBody(parser);
// Propagator builder
final NumericalPropagatorBuilder propagatorBuilder = createPropagatorBuilder(parser, conventions, gravityField, body, initialGuess);
// Measurements
final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
for (final String fileName : parser.getStringsList(ParameterKey.MEASUREMENTS_FILES, ',')) {
measurements.addAll(readMeasurements(new File(input.getParentFile(), fileName), createStationsData(parser, body), createPVData(parser), createSatRangeBias(parser), createWeights(parser), createRangeOutliersManager(parser), createRangeRateOutliersManager(parser), createAzElOutliersManager(parser), createPVOutliersManager(parser)));
}
// Building the Kalman filter:
// - Gather the estimated measurement parameters in a list
// - Prepare the initial covariance matrix and the process noise matrix
// - Build the Kalman filter
// --------------------------------------------------------------------
// Build the list of estimated measurements
final ParameterDriversList estimatedMeasurementsParameters = new ParameterDriversList();
for (ObservedMeasurement<?> measurement : measurements) {
final List<ParameterDriver> drivers = measurement.getParametersDrivers();
for (ParameterDriver driver : drivers) {
if (driver.isSelected()) {
// Add the driver
estimatedMeasurementsParameters.add(driver);
}
}
}
// Sort the list lexicographically
estimatedMeasurementsParameters.sort();
// Orbital covariance matrix initialization
// Jacobian of the orbital parameters w/r to Cartesian
final double[][] dYdC = new double[6][6];
initialGuess.getJacobianWrtCartesian(propagatorBuilder.getPositionAngle(), dYdC);
final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
RealMatrix orbitalP = Jac.multiply(cartesianOrbitalP.multiply(Jac.transpose()));
// Orbital process noise matrix
RealMatrix orbitalQ = Jac.multiply(cartesianOrbitalQ.multiply(Jac.transpose()));
// Build the full covariance matrix and process noise matrix
final int nbPropag = (propagationP != null) ? propagationP.getRowDimension() : 0;
final int nbMeas = (measurementP != null) ? measurementP.getRowDimension() : 0;
final RealMatrix initialP = MatrixUtils.createRealMatrix(6 + nbPropag + nbMeas, 6 + nbPropag + nbMeas);
final RealMatrix Q = MatrixUtils.createRealMatrix(6 + nbPropag + nbMeas, 6 + nbPropag + nbMeas);
// Orbital part
initialP.setSubMatrix(orbitalP.getData(), 0, 0);
Q.setSubMatrix(orbitalQ.getData(), 0, 0);
// Propagation part
if (propagationP != null) {
initialP.setSubMatrix(propagationP.getData(), 6, 6);
Q.setSubMatrix(propagationQ.getData(), 6, 6);
}
// Measurement part
if (measurementP != null) {
initialP.setSubMatrix(measurementP.getData(), 6 + nbPropag, 6 + nbPropag);
Q.setSubMatrix(measurementQ.getData(), 6 + nbPropag, 6 + nbPropag);
}
// Build the Kalman
KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
kalmanBuilder.builder(propagatorBuilder);
kalmanBuilder.estimatedMeasurementsParameters(estimatedMeasurementsParameters);
kalmanBuilder.initialCovarianceMatrix(initialP);
kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(Q));
final KalmanEstimator kalman = kalmanBuilder.build();
// Add an observer
kalman.setObserver(new KalmanObserver() {
/**
* Date of the first measurement.
*/
private AbsoluteDate t0;
/**
* {@inheritDoc}
* @throws OrekitException
*/
@Override
@SuppressWarnings("unchecked")
public void evaluationPerformed(final KalmanEstimation estimation) throws OrekitException {
// Current measurement number, date and status
final EstimatedMeasurement<?> estimatedMeasurement = estimation.getCorrectedMeasurement();
final int currentNumber = estimation.getCurrentMeasurementNumber();
final AbsoluteDate currentDate = estimatedMeasurement.getDate();
final EstimatedMeasurement.Status currentStatus = estimatedMeasurement.getStatus();
// Current estimated measurement
final ObservedMeasurement<?> observedMeasurement = estimatedMeasurement.getObservedMeasurement();
// Measurement type & Station name
String measType = "";
String stationName = "";
// Register the measurement in the proper measurement logger
if (observedMeasurement instanceof Range) {
// Add the tuple (estimation, prediction) to the log
rangeLog.add(currentNumber, (EstimatedMeasurement<Range>) estimatedMeasurement);
// Measurement type & Station name
measType = "RANGE";
stationName = ((EstimatedMeasurement<Range>) estimatedMeasurement).getObservedMeasurement().getStation().getBaseFrame().getName();
} else if (observedMeasurement instanceof RangeRate) {
rangeRateLog.add(currentNumber, (EstimatedMeasurement<RangeRate>) estimatedMeasurement);
measType = "RANGE_RATE";
stationName = ((EstimatedMeasurement<RangeRate>) estimatedMeasurement).getObservedMeasurement().getStation().getBaseFrame().getName();
} else if (observedMeasurement instanceof AngularAzEl) {
azimuthLog.add(currentNumber, (EstimatedMeasurement<AngularAzEl>) estimatedMeasurement);
elevationLog.add(currentNumber, (EstimatedMeasurement<AngularAzEl>) estimatedMeasurement);
measType = "AZ_EL";
stationName = ((EstimatedMeasurement<AngularAzEl>) estimatedMeasurement).getObservedMeasurement().getStation().getBaseFrame().getName();
} else if (observedMeasurement instanceof PV) {
positionLog.add(currentNumber, (EstimatedMeasurement<PV>) estimatedMeasurement);
velocityLog.add(currentNumber, (EstimatedMeasurement<PV>) estimatedMeasurement);
measType = "PV";
}
// Header
if (print) {
if (currentNumber == 1) {
// Set t0 to first measurement date
t0 = currentDate;
// Print header
final String formatHeader = "%-4s\t%-25s\t%15s\t%-10s\t%-10s\t%-20s\t%20s\t%20s";
String header = String.format(Locale.US, formatHeader, "Nb", "Epoch", "Dt[s]", "Status", "Type", "Station", "DP Corr", "DV Corr");
// Orbital drivers
for (DelegatingDriver driver : estimation.getEstimatedOrbitalParameters().getDrivers()) {
header += String.format(Locale.US, "\t%20s", driver.getName());
header += String.format(Locale.US, "\t%20s", "D" + driver.getName());
}
// Propagation drivers
for (DelegatingDriver driver : estimation.getEstimatedPropagationParameters().getDrivers()) {
header += String.format(Locale.US, "\t%20s", driver.getName());
header += String.format(Locale.US, "\t%20s", "D" + driver.getName());
}
// Measurements drivers
for (DelegatingDriver driver : estimation.getEstimatedMeasurementsParameters().getDrivers()) {
header += String.format(Locale.US, "\t%20s", driver.getName());
header += String.format(Locale.US, "\t%20s", "D" + driver.getName());
}
// Print header
System.out.println(header);
}
// Print current measurement info in terminal
String line = "";
// Line format
final String lineFormat = "%4d\t%-25s\t%15.3f\t%-10s\t%-10s\t%-20s\t%20.9e\t%20.9e";
// Orbital correction = DP & DV between predicted orbit and estimated orbit
final Vector3D predictedP = estimation.getPredictedSpacecraftStates()[0].getPVCoordinates().getPosition();
final Vector3D predictedV = estimation.getPredictedSpacecraftStates()[0].getPVCoordinates().getVelocity();
final Vector3D estimatedP = estimation.getCorrectedSpacecraftStates()[0].getPVCoordinates().getPosition();
final Vector3D estimatedV = estimation.getCorrectedSpacecraftStates()[0].getPVCoordinates().getVelocity();
final double DPcorr = Vector3D.distance(predictedP, estimatedP);
final double DVcorr = Vector3D.distance(predictedV, estimatedV);
line = String.format(Locale.US, lineFormat, currentNumber, currentDate.toString(), currentDate.durationFrom(t0), currentStatus.toString(), measType, stationName, DPcorr, DVcorr);
// Handle parameters printing (value and error)
int jPar = 0;
final RealMatrix Pest = estimation.getPhysicalEstimatedCovarianceMatrix();
// Orbital drivers
for (DelegatingDriver driver : estimation.getEstimatedOrbitalParameters().getDrivers()) {
line += String.format(Locale.US, "\t%20.9f", driver.getValue());
line += String.format(Locale.US, "\t%20.9e", FastMath.sqrt(Pest.getEntry(jPar, jPar)));
jPar++;
}
// Propagation drivers
for (DelegatingDriver driver : estimation.getEstimatedPropagationParameters().getDrivers()) {
line += String.format(Locale.US, "\t%20.9f", driver.getValue());
line += String.format(Locale.US, "\t%20.9e", FastMath.sqrt(Pest.getEntry(jPar, jPar)));
jPar++;
}
// Measurements drivers
for (DelegatingDriver driver : estimatedMeasurementsParameters.getDrivers()) {
line += String.format(Locale.US, "\t%20.9f", driver.getValue());
line += String.format(Locale.US, "\t%20.9e", FastMath.sqrt(Pest.getEntry(jPar, jPar)));
jPar++;
}
// Print the line
System.out.println(line);
}
}
});
// Process the list measurements
final Orbit estimated = kalman.processMeasurements(measurements).getInitialState().getOrbit();
// Get the last estimated physical covariances
final RealMatrix covarianceMatrix = kalman.getPhysicalEstimatedCovarianceMatrix();
// Parameters and measurements.
final ParameterDriversList propagationParameters = kalman.getPropagationParametersDrivers(true);
final ParameterDriversList measurementsParameters = kalman.getEstimatedMeasurementsParameters();
// Eventually, print parameter changes, statistics and covariances
if (print) {
// Display parameter change for non orbital drivers
int length = 0;
for (final ParameterDriver parameterDriver : propagationParameters.getDrivers()) {
length = FastMath.max(length, parameterDriver.getName().length());
}
for (final ParameterDriver parameterDriver : measurementsParameters.getDrivers()) {
length = FastMath.max(length, parameterDriver.getName().length());
}
if (propagationParameters.getNbParams() > 0) {
displayParametersChanges(System.out, "Estimated propagator parameters changes: ", true, length, propagationParameters);
}
if (measurementsParameters.getNbParams() > 0) {
displayParametersChanges(System.out, "Estimated measurements parameters changes: ", true, length, measurementsParameters);
}
// Measurements statistics summary
System.out.println("");
rangeLog.displaySummary(System.out);
rangeRateLog.displaySummary(System.out);
azimuthLog.displaySummary(System.out);
elevationLog.displaySummary(System.out);
positionLog.displaySummary(System.out);
velocityLog.displaySummary(System.out);
// Covariances and sigmas
displayFinalCovariances(System.out, kalman);
}
// Instantiation of the results
return new ResultKalman(propagationParameters, measurementsParameters, kalman.getCurrentMeasurementNumber(), estimated.getPVCoordinates(), rangeLog.createStatisticsSummary(), rangeRateLog.createStatisticsSummary(), azimuthLog.createStatisticsSummary(), elevationLog.createStatisticsSummary(), positionLog.createStatisticsSummary(), velocityLog.createStatisticsSummary(), covarianceMatrix);
}
Aggregations