use of org.hipparchus.analysis.RealFieldUnivariateFunction 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;
}
}
Aggregations