Search in sources :

Example 1 with Line

use of org.hipparchus.geometry.euclidean.threed.Line 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 2 with Line

use of org.hipparchus.geometry.euclidean.threed.Line 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 3 with Line

use of org.hipparchus.geometry.euclidean.threed.Line in project Orekit by CS-SI.

the class Transform method transformLine.

/**
 * Transform a line.
 * @param line to transform
 * @return transformed line
 */
public Line transformLine(final Line line) {
    final Vector3D transformedP0 = transformPosition(line.getOrigin());
    final Vector3D transformedP1 = transformPosition(line.pointAt(1.0e6));
    return new Line(transformedP0, transformedP1, 1.0e-10);
}
Also used : Line(org.hipparchus.geometry.euclidean.threed.Line) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D)

Example 4 with Line

use of org.hipparchus.geometry.euclidean.threed.Line in project Orekit by CS-SI.

the class OneAxisEllipsoid method getIntersectionPoint.

/**
 * {@inheritDoc}
 */
public GeodeticPoint getIntersectionPoint(final Line line, final Vector3D close, final Frame frame, final AbsoluteDate date) throws OrekitException {
    // transform line and close to body frame
    final Transform frameToBodyFrame = frame.getTransformTo(bodyFrame, date);
    final Line lineInBodyFrame = frameToBodyFrame.transformLine(line);
    final Vector3D closeInBodyFrame = frameToBodyFrame.transformPosition(close);
    final double closeAbscissa = lineInBodyFrame.getAbscissa(closeInBodyFrame);
    // compute some miscellaneous variables outside of the loop
    final Vector3D point = lineInBodyFrame.getOrigin();
    final double x = point.getX();
    final double y = point.getY();
    final double z = point.getZ();
    final double z2 = z * z;
    final double r2 = x * x + y * y;
    final Vector3D direction = lineInBodyFrame.getDirection();
    final double dx = direction.getX();
    final double dy = direction.getY();
    final double dz = direction.getZ();
    final double cz2 = dx * dx + dy * dy;
    // abscissa of the intersection as a root of a 2nd degree polynomial :
    // a k^2 - 2 b k + c = 0
    final double a = 1.0 - e2 * cz2;
    final double b = -(g2 * (x * dx + y * dy) + z * dz);
    final double c = g2 * (r2 - ae2) + z2;
    final double b2 = b * b;
    final double ac = a * c;
    if (b2 < ac) {
        return null;
    }
    final double s = FastMath.sqrt(b2 - ac);
    final double k1 = (b < 0) ? (b - s) / a : c / (b + s);
    final double k2 = c / (a * k1);
    // select the right point
    final double k = (FastMath.abs(k1 - closeAbscissa) < FastMath.abs(k2 - closeAbscissa)) ? k1 : k2;
    final Vector3D intersection = lineInBodyFrame.pointAt(k);
    final double ix = intersection.getX();
    final double iy = intersection.getY();
    final double iz = intersection.getZ();
    final double lambda = FastMath.atan2(iy, ix);
    final double phi = FastMath.atan2(iz, g2 * FastMath.sqrt(ix * ix + iy * iy));
    return new GeodeticPoint(phi, lambda, 0.0);
}
Also used : Line(org.hipparchus.geometry.euclidean.threed.Line) FieldLine(org.hipparchus.geometry.euclidean.threed.FieldLine) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Transform(org.orekit.frames.Transform) FieldTransform(org.orekit.frames.FieldTransform)

Example 5 with Line

use of org.hipparchus.geometry.euclidean.threed.Line in project Orekit by CS-SI.

the class TransformTest method testLine.

@Test
public void testLine() {
    RandomGenerator random = new Well19937a(0x4a5ff67426c5731fl);
    for (int i = 0; i < 100; ++i) {
        Transform transform = randomTransform(random);
        for (int j = 0; j < 20; ++j) {
            Vector3D p0 = randomVector(1.0e3, random);
            Vector3D p1 = randomVector(1.0e3, random);
            Line l = new Line(p0, p1, 1.0e-10);
            Line transformed = transform.transformLine(l);
            for (int k = 0; k < 10; ++k) {
                Vector3D p = l.pointAt(random.nextDouble() * 1.0e6);
                Assert.assertEquals(0.0, transformed.distance(transform.transformPosition(p)), 1.0e-9);
            }
        }
    }
}
Also used : Line(org.hipparchus.geometry.euclidean.threed.Line) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) Well19937a(org.hipparchus.random.Well19937a) RandomGenerator(org.hipparchus.random.RandomGenerator) Test(org.junit.Test)

Aggregations

Line (org.hipparchus.geometry.euclidean.threed.Line)20 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)19 FieldVector3D (org.hipparchus.geometry.euclidean.threed.FieldVector3D)12 Test (org.junit.Test)11 GeodeticPoint (org.orekit.bodies.GeodeticPoint)10 Frame (org.orekit.frames.Frame)9 Transform (org.orekit.frames.Transform)6 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)6 TimeStampedPVCoordinates (org.orekit.utils.TimeStampedPVCoordinates)6 FieldLine (org.hipparchus.geometry.euclidean.threed.FieldLine)5 AbsoluteDate (org.orekit.time.AbsoluteDate)5 OrekitException (org.orekit.errors.OrekitException)4 PVCoordinates (org.orekit.utils.PVCoordinates)4 Rotation (org.hipparchus.geometry.euclidean.threed.Rotation)3 RandomGenerator (org.hipparchus.random.RandomGenerator)3 Well19937a (org.hipparchus.random.Well19937a)3 FieldGeodeticPoint (org.orekit.bodies.FieldGeodeticPoint)3 FieldTransform (org.orekit.frames.FieldTransform)3 DateComponents (org.orekit.time.DateComponents)3 TimeStampedFieldPVCoordinates (org.orekit.utils.TimeStampedFieldPVCoordinates)3