Search in sources :

Example 1 with MathRuntimeException

use of org.hipparchus.exception.MathRuntimeException in project Orekit by CS-SI.

the class Phasing method findLatitudeCrossing.

/**
 * Find the state at which the reference latitude is crossed.
 * @param latitude latitude to search for
 * @param guessDate guess date for the crossing
 * @param endDate maximal date not to overtake
 * @param shift shift value used to evaluate the latitude function bracketing around the guess date
 * @param maxShift maximum value that the shift value can take
 * @param propagator propagator used
 * @return state at latitude crossing time
 * @throws OrekitException if state cannot be propagated
 * @throws MathRuntimeException if latitude cannot be bracketed in the search interval
 */
private SpacecraftState findLatitudeCrossing(final double latitude, final AbsoluteDate guessDate, final AbsoluteDate endDate, final double shift, final double maxShift, final Propagator propagator) throws OrekitException, MathRuntimeException {
    // function evaluating to 0 at latitude crossings
    final UnivariateFunction latitudeFunction = new UnivariateFunction() {

        /**
         * {@inheritDoc}
         */
        public double value(double x) {
            try {
                final SpacecraftState state = propagator.propagate(guessDate.shiftedBy(x));
                final Vector3D position = state.getPVCoordinates(earth.getBodyFrame()).getPosition();
                final GeodeticPoint point = earth.transform(position, earth.getBodyFrame(), state.getDate());
                return point.getLatitude() - latitude;
            } catch (OrekitException oe) {
                throw new RuntimeException(oe);
            }
        }
    };
    // try to bracket the encounter
    double span;
    if (guessDate.shiftedBy(shift).compareTo(endDate) > 0) {
        // Take a 1e-3 security margin
        span = endDate.durationFrom(guessDate) - 1e-3;
    } else {
        span = shift;
    }
    while (!UnivariateSolverUtils.isBracketing(latitudeFunction, -span, span)) {
        if (2 * span > maxShift) {
            // let the Hipparchus exception be thrown
            UnivariateSolverUtils.verifyBracketing(latitudeFunction, -span, span);
        } else if (guessDate.shiftedBy(2 * span).compareTo(endDate) > 0) {
            // Out of range :
            return null;
        }
        // expand the search interval
        span *= 2;
    }
    // find the encounter in the bracketed interval
    final BaseUnivariateSolver<UnivariateFunction> solver = new BracketingNthOrderBrentSolver(0.1, 5);
    final double dt = solver.solve(1000, latitudeFunction, -span, span);
    return propagator.propagate(guessDate.shiftedBy(dt));
}
Also used : SpacecraftState(org.orekit.propagation.SpacecraftState) MathRuntimeException(org.hipparchus.exception.MathRuntimeException) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) OrekitException(org.orekit.errors.OrekitException) GeodeticPoint(org.orekit.bodies.GeodeticPoint) BracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver)

Example 2 with MathRuntimeException

use of org.hipparchus.exception.MathRuntimeException in project Orekit by CS-SI.

the class Geoid method getIntersectionPoint.

/**
 * {@inheritDoc}
 *
 * <p> The intersection point is computed using a line search along the
 * specified line. This is accurate when the geoid is slowly varying.
 */
@Override
public GeodeticPoint getIntersectionPoint(final Line lineInFrame, final Vector3D closeInFrame, final Frame frame, final AbsoluteDate date) throws OrekitException {
    /*
         * It is assumed that the geoid is slowly varying over it's entire
         * surface. Therefore there will one local intersection.
         */
    // transform to body frame
    final Frame bodyFrame = this.getBodyFrame();
    final Transform frameToBody = frame.getTransformTo(bodyFrame, date);
    final Vector3D close = frameToBody.transformPosition(closeInFrame);
    final Line lineInBodyFrame = frameToBody.transformLine(lineInFrame);
    // set the line's direction so the solved for value is always positive
    final Line line;
    if (lineInBodyFrame.getAbscissa(close) < 0) {
        line = lineInBodyFrame.revert();
    } else {
        line = lineInBodyFrame;
    }
    final ReferenceEllipsoid ellipsoid = this.getEllipsoid();
    // calculate end points
    // distance from line to center of earth, squared
    final double d2 = line.pointAt(0.0).getNormSq();
    // the minimum abscissa, squared
    final double n = ellipsoid.getPolarRadius() + MIN_UNDULATION;
    final double minAbscissa2 = n * n - d2;
    // smaller end point of the interval = 0.0 or intersection with
    // min_undulation sphere
    final double lowPoint = FastMath.sqrt(FastMath.max(minAbscissa2, 0.0));
    // the maximum abscissa, squared
    final double x = ellipsoid.getEquatorialRadius() + MAX_UNDULATION;
    final double maxAbscissa2 = x * x - d2;
    // larger end point of the interval
    final double highPoint = FastMath.sqrt(maxAbscissa2);
    // line search function
    final UnivariateFunction heightFunction = new UnivariateFunction() {

        @Override
        public double value(final double x) {
            try {
                final GeodeticPoint geodetic = transform(line.pointAt(x), bodyFrame, date);
                return geodetic.getAltitude();
            } catch (OrekitException e) {
                // due to frame transform -> re-throw
                throw new RuntimeException(e);
            }
        }
    };
    // compute answer
    if (maxAbscissa2 < 0) {
        // ray does not pierce bounding sphere -> no possible intersection
        return null;
    }
    // solve line search problem to find the intersection
    final UnivariateSolver solver = new BracketingNthOrderBrentSolver();
    try {
        final double abscissa = solver.solve(MAX_EVALUATIONS, heightFunction, lowPoint, highPoint);
        // return intersection point
        return this.transform(line.pointAt(abscissa), bodyFrame, date);
    } catch (MathRuntimeException e) {
        // no intersection
        return null;
    }
}
Also used : Frame(org.orekit.frames.Frame) MathRuntimeException(org.hipparchus.exception.MathRuntimeException) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) RealFieldUnivariateFunction(org.hipparchus.analysis.RealFieldUnivariateFunction) UnivariateSolver(org.hipparchus.analysis.solvers.UnivariateSolver) FieldBracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver) BracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver) Line(org.hipparchus.geometry.euclidean.threed.Line) FieldLine(org.hipparchus.geometry.euclidean.threed.FieldLine) MathRuntimeException(org.hipparchus.exception.MathRuntimeException) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) OrekitException(org.orekit.errors.OrekitException) Transform(org.orekit.frames.Transform) FieldTransform(org.orekit.frames.FieldTransform) GeodeticPoint(org.orekit.bodies.GeodeticPoint) FieldGeodeticPoint(org.orekit.bodies.FieldGeodeticPoint)

Example 3 with MathRuntimeException

use of org.hipparchus.exception.MathRuntimeException in project Orekit by CS-SI.

the class Geoid method getIntersectionPoint.

/**
 * {@inheritDoc}
 *
 * <p> The intersection point is computed using a line search along the
 * specified line. This is accurate when the geoid is slowly varying.
 */
@Override
public <T extends RealFieldElement<T>> FieldGeodeticPoint<T> getIntersectionPoint(final FieldLine<T> lineInFrame, final FieldVector3D<T> closeInFrame, final Frame frame, final FieldAbsoluteDate<T> date) throws OrekitException {
    final Field<T> field = date.getField();
    /*
         * It is assumed that the geoid is slowly varying over it's entire
         * surface. Therefore there will one local intersection.
         */
    // transform to body frame
    final Frame bodyFrame = this.getBodyFrame();
    final FieldTransform<T> frameToBody = frame.getTransformTo(bodyFrame, date);
    final FieldVector3D<T> close = frameToBody.transformPosition(closeInFrame);
    final FieldLine<T> lineInBodyFrame = frameToBody.transformLine(lineInFrame);
    // set the line's direction so the solved for value is always positive
    final FieldLine<T> line;
    if (lineInBodyFrame.getAbscissa(close).getReal() < 0) {
        line = lineInBodyFrame.revert();
    } else {
        line = lineInBodyFrame;
    }
    final ReferenceEllipsoid ellipsoid = this.getEllipsoid();
    // calculate end points
    // distance from line to center of earth, squared
    final T d2 = line.pointAt(0.0).getNormSq();
    // the minimum abscissa, squared
    final double n = ellipsoid.getPolarRadius() + MIN_UNDULATION;
    final T minAbscissa2 = d2.negate().add(n * n);
    // smaller end point of the interval = 0.0 or intersection with
    // min_undulation sphere
    final T lowPoint = minAbscissa2.getReal() < 0 ? field.getZero() : minAbscissa2.sqrt();
    // the maximum abscissa, squared
    final double x = ellipsoid.getEquatorialRadius() + MAX_UNDULATION;
    final T maxAbscissa2 = d2.negate().add(x * x);
    // larger end point of the interval
    final T highPoint = maxAbscissa2.sqrt();
    // line search function
    final RealFieldUnivariateFunction<T> heightFunction = z -> {
        try {
            final FieldGeodeticPoint<T> geodetic = transform(line.pointAt(z), bodyFrame, date);
            return geodetic.getAltitude();
        } catch (OrekitException e) {
            // due to frame transform -> re-throw
            throw new RuntimeException(e);
        }
    };
    // compute answer
    if (maxAbscissa2.getReal() < 0) {
        // ray does not pierce bounding sphere -> no possible intersection
        return null;
    }
    // solve line search problem to find the intersection
    final FieldBracketingNthOrderBrentSolver<T> solver = new FieldBracketingNthOrderBrentSolver<>(field.getZero().add(1.0e-14), field.getZero().add(1.0e-6), field.getZero().add(1.0e-15), 5);
    try {
        final T abscissa = solver.solve(MAX_EVALUATIONS, heightFunction, lowPoint, highPoint, AllowedSolution.ANY_SIDE);
        // return intersection point
        return this.transform(line.pointAt(abscissa), bodyFrame, date);
    } catch (MathRuntimeException e) {
        // no intersection
        return null;
    }
}
Also used : AllowedSolution(org.hipparchus.analysis.solvers.AllowedSolution) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) GeodeticPoint(org.orekit.bodies.GeodeticPoint) FieldGeodeticPoint(org.orekit.bodies.FieldGeodeticPoint) MathRuntimeException(org.hipparchus.exception.MathRuntimeException) Frame(org.orekit.frames.Frame) NormalizedSphericalHarmonicsProvider(org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) FastMath(org.hipparchus.util.FastMath) FieldBracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver) BracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) UnivariateFunction(org.hipparchus.analysis.UnivariateFunction) TideSystem(org.orekit.forces.gravity.potential.TideSystem) UnivariateSolver(org.hipparchus.analysis.solvers.UnivariateSolver) Line(org.hipparchus.geometry.euclidean.threed.Line) Field(org.hipparchus.Field) RealFieldUnivariateFunction(org.hipparchus.analysis.RealFieldUnivariateFunction) OrekitException(org.orekit.errors.OrekitException) RealFieldElement(org.hipparchus.RealFieldElement) FieldLine(org.hipparchus.geometry.euclidean.threed.FieldLine) Transform(org.orekit.frames.Transform) HolmesFeatherstoneAttractionModel(org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel) FieldTransform(org.orekit.frames.FieldTransform) AbsoluteDate(org.orekit.time.AbsoluteDate) Frame(org.orekit.frames.Frame) FieldGeodeticPoint(org.orekit.bodies.FieldGeodeticPoint) MathRuntimeException(org.hipparchus.exception.MathRuntimeException) FieldBracketingNthOrderBrentSolver(org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver) MathRuntimeException(org.hipparchus.exception.MathRuntimeException) OrekitException(org.orekit.errors.OrekitException)

Example 4 with MathRuntimeException

use of org.hipparchus.exception.MathRuntimeException 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 5 with MathRuntimeException

use of org.hipparchus.exception.MathRuntimeException 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)

Aggregations

MathRuntimeException (org.hipparchus.exception.MathRuntimeException)11 OrekitException (org.orekit.errors.OrekitException)8 UnivariateFunction (org.hipparchus.analysis.UnivariateFunction)5 BracketingNthOrderBrentSolver (org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver)5 GeodeticPoint (org.orekit.bodies.GeodeticPoint)5 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)4 OrekitExceptionWrapper (org.orekit.errors.OrekitExceptionWrapper)4 SpacecraftState (org.orekit.propagation.SpacecraftState)4 UnivariateSolver (org.hipparchus.analysis.solvers.UnivariateSolver)3 RealFieldUnivariateFunction (org.hipparchus.analysis.RealFieldUnivariateFunction)2 FieldBracketingNthOrderBrentSolver (org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver)2 FieldLine (org.hipparchus.geometry.euclidean.threed.FieldLine)2 FieldVector3D (org.hipparchus.geometry.euclidean.threed.FieldVector3D)2 Line (org.hipparchus.geometry.euclidean.threed.Line)2 Test (org.junit.Test)2 FieldGeodeticPoint (org.orekit.bodies.FieldGeodeticPoint)2 FieldTransform (org.orekit.frames.FieldTransform)2 Frame (org.orekit.frames.Frame)2 Transform (org.orekit.frames.Transform)2 AbsoluteDate (org.orekit.time.AbsoluteDate)2