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));
}
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;
}
}
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;
}
}
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();
}
}
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();
}
}
Aggregations