use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.
the class BatchLSEstimator method estimate.
/**
* Estimate the orbital, propagation and measurements parameters.
* <p>
* The initial guess for all parameters must have been set before calling this method
* using {@link #getOrbitalParametersDrivers(boolean)}, {@link #getPropagatorParametersDrivers(boolean)},
* and {@link #getMeasurementsParametersDrivers(boolean)} and then {@link ParameterDriver#setValue(double)
* setting the values} of the parameters.
* </p>
* <p>
* For parameters whose reference date has not been set to a non-null date beforehand (i.e.
* the parameters for which {@link ParameterDriver#getReferenceDate()} returns {@code null},
* a default reference date will be set automatically at the start of the estimation to the
* {@link NumericalPropagatorBuilder#getInitialOrbitDate() initial orbit date} of the first
* propagator builder. For parameters whose reference date has been set to a non-null date,
* this reference date is untouched.
* </p>
* <p>
* After this method returns, the estimated parameters can be retrieved using
* {@link #getOrbitalParametersDrivers(boolean)}, {@link #getPropagatorParametersDrivers(boolean)},
* and {@link #getMeasurementsParametersDrivers(boolean)} and then {@link ParameterDriver#getValue()
* getting the values} of the parameters.
* </p>
* <p>
* As a convenience, the method also returns a fully configured and ready to use
* propagator set up with all the estimated values.
* </p>
* <p>
* For even more in-depth information, the {@link #getOptimum()} method provides detailed
* elements (covariance matrix, estimated parameters standard deviation, weighted Jacobian, RMS,
* χ², residuals and more).
* </p>
* @return propagators configured with estimated orbits as initial states, and all
* propagators estimated parameters also set
* @exception OrekitException if there is a conflict in parameters names
* or if orbit cannot be determined
*/
public NumericalPropagator[] estimate() throws OrekitException {
// set reference date for all parameters that lack one (including the not estimated parameters)
for (final ParameterDriver driver : getOrbitalParametersDrivers(false).getDrivers()) {
if (driver.getReferenceDate() == null) {
driver.setReferenceDate(builders[0].getInitialOrbitDate());
}
}
for (final ParameterDriver driver : getPropagatorParametersDrivers(false).getDrivers()) {
if (driver.getReferenceDate() == null) {
driver.setReferenceDate(builders[0].getInitialOrbitDate());
}
}
for (final ParameterDriver driver : getMeasurementsParametersDrivers(false).getDrivers()) {
if (driver.getReferenceDate() == null) {
driver.setReferenceDate(builders[0].getInitialOrbitDate());
}
}
// get all estimated parameters
final ParameterDriversList estimatedOrbitalParameters = getOrbitalParametersDrivers(true);
final ParameterDriversList estimatedPropagatorParameters = getPropagatorParametersDrivers(true);
final ParameterDriversList estimatedMeasurementsParameters = getMeasurementsParametersDrivers(true);
// create start point
final double[] start = new double[estimatedOrbitalParameters.getNbParams() + estimatedPropagatorParameters.getNbParams() + estimatedMeasurementsParameters.getNbParams()];
int iStart = 0;
for (final ParameterDriver driver : estimatedOrbitalParameters.getDrivers()) {
start[iStart++] = driver.getNormalizedValue();
}
for (final ParameterDriver driver : estimatedPropagatorParameters.getDrivers()) {
start[iStart++] = driver.getNormalizedValue();
}
for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
start[iStart++] = driver.getNormalizedValue();
}
lsBuilder.start(start);
// create target (which is an array set to 0, as we compute weighted residuals ourselves)
int p = 0;
for (final ObservedMeasurement<?> measurement : measurements) {
if (measurement.isEnabled()) {
p += measurement.getDimension();
}
}
final double[] target = new double[p];
lsBuilder.target(target);
// set up the model
final ModelObserver modelObserver = new ModelObserver() {
/**
* {@inheritDoc}
*/
@Override
public void modelCalled(final Orbit[] newOrbits, final Map<ObservedMeasurement<?>, EstimatedMeasurement<?>> newEstimations) {
BatchLSEstimator.this.orbits = newOrbits;
BatchLSEstimator.this.estimations = newEstimations;
}
};
final Model model = new Model(builders, measurements, estimatedMeasurementsParameters, modelObserver);
lsBuilder.model(model);
// add a validator for orbital parameters
lsBuilder.parameterValidator(new Validator(estimatedOrbitalParameters, estimatedPropagatorParameters, estimatedMeasurementsParameters));
lsBuilder.checker(new ConvergenceChecker<LeastSquaresProblem.Evaluation>() {
/**
* {@inheritDoc}
*/
@Override
public boolean converged(final int iteration, final LeastSquaresProblem.Evaluation previous, final LeastSquaresProblem.Evaluation current) {
final double lInf = current.getPoint().getLInfDistance(previous.getPoint());
return lInf <= parametersConvergenceThreshold;
}
});
// set up the problem to solve
final LeastSquaresProblem problem = new TappedLSProblem(lsBuilder.build(), model, estimatedOrbitalParameters, estimatedPropagatorParameters, estimatedMeasurementsParameters);
try {
// solve the problem
optimum = optimizer.optimize(problem);
// create a new configured propagator with all estimated parameters
return model.createPropagators(optimum.getPoint());
} catch (MathRuntimeException mrte) {
throw new OrekitException(mrte);
} catch (OrekitExceptionWrapper oew) {
throw oew.getException();
}
}
use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.
the class Differentiation method differentiate.
/**
* Differentiate a scalar function using finite differences.
* @param function function to differentiate
* @param driver driver for the parameter
* @param nbPoints number of points used for finite differences
* @param step step for finite differences
* @return scalar function evaluating to the derivative of the original function
*/
public static ParameterFunction differentiate(final ParameterFunction function, final ParameterDriver driver, final int nbPoints, final double step) {
final UnivariateFunction uf = new UnivariateFunction() {
/**
* {@inheritDoc}
*/
@Override
public double value(final double normalizedValue) throws OrekitExceptionWrapper {
try {
final double saved = driver.getNormalizedValue();
driver.setNormalizedValue(normalizedValue);
final double functionValue = function.value(driver);
driver.setNormalizedValue(saved);
return functionValue;
} catch (OrekitException oe) {
throw new OrekitExceptionWrapper(oe);
}
}
};
final FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(nbPoints, step);
final UnivariateDifferentiableFunction differentiated = differentiator.differentiate(uf);
return new ParameterFunction() {
/**
* {@inheritDoc}
*/
@Override
public double value(final ParameterDriver parameterDriver) throws OrekitException {
if (!parameterDriver.getName().equals(driver.getName())) {
throw new OrekitException(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, parameterDriver.getName(), driver.getName());
}
try {
final DerivativeStructure dsParam = FACTORY.variable(0, parameterDriver.getNormalizedValue());
final DerivativeStructure dsValue = differentiated.value(dsParam);
return dsValue.getPartialDerivative(1);
} catch (OrekitExceptionWrapper oew) {
throw oew.getException();
}
}
};
}
use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.
the class TopocentricFrame method computeLimitVisibilityPoint.
/**
* Compute the limit visibility point for a satellite in a given direction.
* <p>
* This method can be used to compute visibility circles around ground stations
* for example, using a simple loop on azimuth, with either a fixed elevation
* or an elevation that depends on azimuth to take ground masks into account.
* </p>
* @param radius satellite distance to Earth center
* @param azimuth pointing azimuth from station
* @param elevation pointing elevation from station
* @return limit visibility point for the satellite
* @throws OrekitException if point cannot be found
*/
public GeodeticPoint computeLimitVisibilityPoint(final double radius, final double azimuth, final double elevation) throws OrekitException {
try {
// convergence threshold on point position: 1mm
final double deltaP = 0.001;
final UnivariateSolver solver = new BracketingNthOrderBrentSolver(deltaP / Constants.WGS84_EARTH_EQUATORIAL_RADIUS, deltaP, deltaP, 5);
// find the distance such that a point in the specified direction and at the solved-for
// distance is exactly at the specified radius
final double distance = solver.solve(1000, new UnivariateFunction() {
/**
* {@inheritDoc}
*/
public double value(final double x) {
try {
final GeodeticPoint gp = pointAtDistance(azimuth, elevation, x);
return parentShape.transform(gp).getNorm() - radius;
} catch (OrekitException oe) {
throw new OrekitExceptionWrapper(oe);
}
}
}, 0, 2 * radius);
// return the limit point
return pointAtDistance(azimuth, elevation, distance);
} catch (MathRuntimeException mrte) {
throw new OrekitException(mrte);
} catch (OrekitExceptionWrapper lwe) {
throw lwe.getException();
}
}
use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.
the class FieldTransformGenerator method generate.
/**
* {@inheritDoc}
*/
public List<FieldTransform<T>> generate(final AbsoluteDate existingDate, final AbsoluteDate date) {
try {
final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
final List<FieldTransform<T>> generated = new ArrayList<>();
if (existingDate == null) {
// no prior existing transforms, just generate a first set
for (int i = 0; i < neighborsSize; ++i) {
generated.add(provider.getTransform(fieldDate.shiftedBy(i * step)));
}
} else {
// some transforms have already been generated
// add the missing ones up to specified date
FieldAbsoluteDate<T> t = new FieldAbsoluteDate<>(field, existingDate);
if (date.compareTo(t.toAbsoluteDate()) > 0) {
// forward generation
do {
t = t.shiftedBy(step);
generated.add(generated.size(), provider.getTransform(t));
} while (t.compareTo(fieldDate) <= 0);
} else {
// backward generation
do {
t = t.shiftedBy(-step);
generated.add(0, provider.getTransform(t));
} while (t.compareTo(fieldDate) >= 0);
}
}
// return the generated transforms
return generated;
} catch (OrekitException oe) {
throw new OrekitExceptionWrapper(oe);
}
}
use of org.orekit.errors.OrekitExceptionWrapper in project Orekit by CS-SI.
the class EventState method findRoot.
/**
* Find a root in a bracketing interval.
*
* <p> When calling this method one of the following must be true. Either ga == 0, gb
* == 0, (ga < 0 and gb > 0), or (ga > 0 and gb < 0).
*
* @param interpolator that covers the interval.
* @param ta earliest possible time for root.
* @param ga g(ta).
* @param tb latest possible time for root.
* @param gb g(tb).
* @return if a zero crossing was found.
* @throws OrekitException if the event detector throws one
*/
private boolean findRoot(final OrekitStepInterpolator interpolator, final AbsoluteDate ta, final double ga, final AbsoluteDate tb, final double gb) throws OrekitException {
// check there appears to be a root in [ta, tb]
check(ga == 0.0 || gb == 0.0 || (ga > 0.0 && gb < 0.0) || (ga < 0.0 && gb > 0.0));
final double convergence = detector.getThreshold();
final int maxIterationCount = detector.getMaxIterationCount();
final BracketedUnivariateSolver<UnivariateFunction> solver = new BracketingNthOrderBrentSolver(0, convergence, 0, 5);
// event time, just at or before the actual root.
AbsoluteDate beforeRootT = null;
double beforeRootG = Double.NaN;
// time on the other side of the root.
// Initialized the the loop below executes once.
AbsoluteDate afterRootT = ta;
double afterRootG = 0.0;
// the ga == 0.0 case is handled by the loop below
if (ta.equals(tb)) {
// both non-zero but times are the same. Probably due to reset state
beforeRootT = ta;
beforeRootG = ga;
afterRootT = shiftedBy(beforeRootT, convergence);
afterRootG = g(interpolator.getInterpolatedState(afterRootT));
} else if (ga != 0.0 && gb == 0.0) {
// hard: ga != 0.0 and gb == 0.0
// look past gb by up to convergence to find next sign
// throw an exception if g(t) = 0.0 in [tb, tb + convergence]
beforeRootT = tb;
beforeRootG = gb;
afterRootT = shiftedBy(beforeRootT, convergence);
afterRootG = g(interpolator.getInterpolatedState(afterRootT));
} else if (ga != 0.0) {
final double newGa = g(interpolator.getInterpolatedState(ta));
if (ga > 0 != newGa > 0) {
// both non-zero, step sign change at ta, possibly due to reset state
beforeRootT = ta;
beforeRootG = newGa;
afterRootT = minTime(shiftedBy(beforeRootT, convergence), tb);
afterRootG = g(interpolator.getInterpolatedState(afterRootT));
}
}
// loop to skip through "fake" roots, i.e. where g(t) = g'(t) = 0.0
// executed once if we didn't hit a special case above
AbsoluteDate loopT = ta;
double loopG = ga;
while ((afterRootG == 0.0 || afterRootG > 0.0 == g0Positive) && strictlyAfter(afterRootT, tb)) {
if (loopG == 0.0) {
// ga == 0.0 and gb may or may not be 0.0
// handle the root at ta first
beforeRootT = loopT;
beforeRootG = loopG;
afterRootT = minTime(shiftedBy(beforeRootT, convergence), tb);
afterRootG = g(interpolator.getInterpolatedState(afterRootT));
} else {
// both non-zero, the usual case, use a root finder.
try {
// time zero for evaluating the function f. Needs to be final
final AbsoluteDate fT0 = loopT;
final UnivariateFunction f = dt -> {
try {
return g(interpolator.getInterpolatedState(fT0.shiftedBy(dt)));
} catch (OrekitException oe) {
throw new OrekitExceptionWrapper(oe);
}
};
// tb as a double for use in f
final double tbDouble = tb.durationFrom(fT0);
if (forward) {
final Interval interval = solver.solveInterval(maxIterationCount, f, 0, tbDouble);
beforeRootT = fT0.shiftedBy(interval.getLeftAbscissa());
beforeRootG = interval.getLeftValue();
afterRootT = fT0.shiftedBy(interval.getRightAbscissa());
afterRootG = interval.getRightValue();
} else {
final Interval interval = solver.solveInterval(maxIterationCount, f, tbDouble, 0);
beforeRootT = fT0.shiftedBy(interval.getRightAbscissa());
beforeRootG = interval.getRightValue();
afterRootT = fT0.shiftedBy(interval.getLeftAbscissa());
afterRootG = interval.getLeftValue();
}
} catch (OrekitExceptionWrapper oew) {
throw oew.getException();
}
}
// assume tolerance is 1 ulp
if (beforeRootT.equals(afterRootT)) {
afterRootT = nextAfter(afterRootT);
afterRootG = g(interpolator.getInterpolatedState(afterRootT));
}
// check loop is making some progress
check((forward && afterRootT.compareTo(beforeRootT) > 0) || (!forward && afterRootT.compareTo(beforeRootT) < 0));
// setup next iteration
loopT = afterRootT;
loopG = afterRootG;
}
// figure out the result of root finding, and return accordingly
if (afterRootG == 0.0 || afterRootG > 0.0 == g0Positive) {
// loop gave up and didn't find any crossing within this step
return false;
} else {
// real crossing
check(beforeRootT != null && !Double.isNaN(beforeRootG));
// variation direction, with respect to the integration direction
increasing = !g0Positive;
pendingEventTime = beforeRootT;
stopTime = beforeRootG == 0.0 ? beforeRootT : afterRootT;
pendingEvent = true;
afterEvent = afterRootT;
afterG = afterRootG;
// check increasing set correctly
check(afterG > 0 == increasing);
check(increasing == gb >= ga);
return true;
}
}
Aggregations