Search in sources :

Example 11 with OrekitExceptionWrapper

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();
    }
}
Also used : MathRuntimeException(org.hipparchus.exception.MathRuntimeException) OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) ParameterDriver(org.orekit.utils.ParameterDriver) ParameterDriversList(org.orekit.utils.ParameterDriversList) LeastSquaresProblem(org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem) OrekitException(org.orekit.errors.OrekitException) Map(java.util.Map) ParameterValidator(org.hipparchus.optim.nonlinear.vector.leastsquares.ParameterValidator)

Example 12 with OrekitExceptionWrapper

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();
            }
        }
    };
}
Also used : OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) UnivariateDifferentiableFunction(org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction) OrekitException(org.orekit.errors.OrekitException) FiniteDifferencesDifferentiator(org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator)

Example 13 with OrekitExceptionWrapper

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();
    }
}
Also used : OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) MathRuntimeException(org.hipparchus.exception.MathRuntimeException) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) UnivariateSolver(org.hipparchus.analysis.solvers.UnivariateSolver) OrekitException(org.orekit.errors.OrekitException) GeodeticPoint(org.orekit.bodies.GeodeticPoint) BracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver)

Example 14 with OrekitExceptionWrapper

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);
    }
}
Also used : OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) ArrayList(java.util.ArrayList) OrekitException(org.orekit.errors.OrekitException) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate)

Example 15 with OrekitExceptionWrapper

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;
    }
}
Also used : BracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) MathRuntimeException(org.hipparchus.exception.MathRuntimeException) OrekitInternalError(org.orekit.errors.OrekitInternalError) OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) Interval(org.hipparchus.analysis.solvers.BracketedUnivariateSolver.Interval) Serializable(java.io.Serializable) Precision(org.hipparchus.util.Precision) EventHandler(org.orekit.propagation.events.handlers.EventHandler) OrekitException(org.orekit.errors.OrekitException) SpacecraftState(org.orekit.propagation.SpacecraftState) OrekitStepInterpolator(org.orekit.propagation.sampling.OrekitStepInterpolator) BracketedUnivariateSolver(org.hipparchus.analysis.solvers.BracketedUnivariateSolver) FastMath(org.hipparchus.util.FastMath) AbsoluteDate(org.orekit.time.AbsoluteDate) OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) OrekitException(org.orekit.errors.OrekitException) BracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver) AbsoluteDate(org.orekit.time.AbsoluteDate) Interval(org.hipparchus.analysis.solvers.BracketedUnivariateSolver.Interval)

Aggregations

OrekitExceptionWrapper (org.orekit.errors.OrekitExceptionWrapper)27 OrekitException (org.orekit.errors.OrekitException)25 UnivariateFunction (org.hipparchus.analysis.UnivariateFunction)10 BracketingNthOrderBrentSolver (org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver)9 SpacecraftState (org.orekit.propagation.SpacecraftState)9 AbsoluteDate (org.orekit.time.AbsoluteDate)8 Map (java.util.Map)7 UnivariateSolver (org.hipparchus.analysis.solvers.UnivariateSolver)7 MathRuntimeException (org.hipparchus.exception.MathRuntimeException)6 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)6 Frame (org.orekit.frames.Frame)5 Transform (org.orekit.frames.Transform)5 ArrayList (java.util.ArrayList)4 HashMap (java.util.HashMap)3 DerivativeStructure (org.hipparchus.analysis.differentiation.DerivativeStructure)3 FiniteDifferencesDifferentiator (org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator)3 NumericalPropagator (org.orekit.propagation.numerical.NumericalPropagator)3 UnivariateVectorFunction (org.hipparchus.analysis.UnivariateVectorFunction)2 UnivariateDifferentiableVectorFunction (org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction)2 BracketedUnivariateSolver (org.hipparchus.analysis.solvers.BracketedUnivariateSolver)2