Geoid.java

/* Contributed in the public domain.
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.models.earth;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
import org.hipparchus.analysis.UnivariateFunction;
import org.hipparchus.analysis.solvers.AllowedSolution;
import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
import org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver;
import org.hipparchus.analysis.solvers.UnivariateSolver;
import org.hipparchus.exception.MathRuntimeException;
import org.hipparchus.geometry.euclidean.threed.FieldLine;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Line;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.errors.OrekitException;
import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.potential.GravityFields;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.forces.gravity.potential.TideSystem;
import org.orekit.frames.FieldStaticTransform;
import org.orekit.frames.Frame;
import org.orekit.frames.StaticTransform;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
import org.orekit.utils.TimeStampedPVCoordinates;

/**
 * A geoid is a level surface of the gravity potential of a body. The gravity
 * potential, W, is split so W = U + T, where U is the normal potential (defined
 * by the ellipsoid) and T is the anomalous potential.[3](eq. 2-137)
 *
 * <p> The {@link #getIntersectionPoint(Line, Vector3D, Frame, AbsoluteDate)}
 * method is tailored specifically for Earth's geoid. All of the other methods
 * in this class are general and will work for an arbitrary body.
 *
 * <p> There are several components that are needed to define a geoid[1]:
 *
 * <ul> <li>Geopotential field. These are the coefficients of the spherical
 * harmonics: S<sub>n,m</sub> and C<sub>n,m</sub></li>
 *
 * <li>Reference Ellipsoid. The ellipsoid is used to define the undulation of
 * the geoid (distance between ellipsoid and geoid) and U<sub>0</sub> the value
 * of the normal gravity potential at the surface of the ellipsoid.</li>
 *
 * <li>W<sub>0</sub>, the potential at the geoid. The value of the potential on
 * the level surface. This is taken to be U<sub>0</sub>, the normal gravity
 * potential at the surface of the {@link ReferenceEllipsoid}.</li>
 *
 * <li>Permanent Tide System. This implementation assumes that the geopotential
 * field and the reference ellipsoid use the same permanent tide system. If the
 * assumption is false it will produce errors of about 0.5 m. Conversion between
 * tide systems is a possible improvement.[1,2]</li>
 *
 * <li>Topographic Masses. That is mass outside of the geoid, e.g. mountains.
 * This implementation ignores topographic masses, which causes up to 3m error
 * in the Himalayas, and ~ 1.5m error in the Rockies. This could be improved
 * through the use of DTED and calculating height anomalies or using the
 * correction coefficients.[1]</li> </ul>
 *
 * <p> This implementation also assumes that the normal to the reference
 * ellipsoid is the same as the normal to the geoid. This assumption enables the
 * equation: (height above geoid) = (height above ellipsoid) - (undulation),
 * which is used in {@link #transform(GeodeticPoint)} and {@link
 * #transform(Vector3D, Frame, AbsoluteDate)}.
 *
 * <p> In testing, the error in the undulations calculated by this class were
 * off by less than 3 meters, which matches the assumptions outlined above.
 *
 * <p> References:
 *
 * <ol> <li>Dru A. Smith. There is no such thing as "The" EGM96 geoid: Subtle
 * points on the use of a global geopotential model. IGeS Bulletin No. 8:17-28,
 * 1998. <a href= "http://www.ngs.noaa.gov/PUBS_LIB/EGM96_GEOID_PAPER/egm96_geoid_paper.html"
 * >http://www.ngs.noaa.gov/PUBS_LIB/EGM96_GEOID_PAPER/egm96_geoid_paper.html</a></li>
 *
 * <li> Martin Losch, Verena Seufer. How to Compute Geoid Undulations (Geoid
 * Height Relative to a Given Reference Ellipsoid) from Spherical Harmonic
 * Coefficients for Satellite Altimetry Applications. , 2003. <a
 * href="http://mitgcm.org/~mlosch/geoidcookbook.pdf">mitgcm.org/~mlosch/geoidcookbook.pdf</a>
 * </li>
 *
 * <li>Weikko A. Heiskanen, Helmut Moritz. Physical Geodesy. W. H. Freeman and
 * Company, 1967. (especially sections 2.13 and equation 2-144 Bruns
 * Formula)</li>
 *
 * <li>S. A. Holmes, W. E. Featherstone. A unified approach to the Clenshaw
 * summation and the recursive computation of very high degree and order
 * normalised associated Legendre functions. Journal of Geodesy, 76(5):279,
 * 2002.</li>
 *
 * <li>DMA TR 8350.2. 1984.</li>
 *
 * <li>Department of Defense World Geodetic System 1984. 2000. NIMA TR 8350.2
 * Third Edition, Amendment 1.</li> </ol>
 *
 * @author Evan Ward
 */
public class Geoid implements EarthShape {

    /**
     * uid is date of last modification.
     */
    private static final long serialVersionUID = 20150312L;

    /**
     * A number larger than the largest undulation. Wikipedia says the geoid
     * height is in [-106, 85]. I chose 100 to be safe.
     */
    private static final double MAX_UNDULATION = 100;
    /**
     * A number smaller than the smallest undulation. Wikipedia says the geoid
     * height is in [-106, 85]. I chose -150 to be safe.
     */
    private static final double MIN_UNDULATION = -150;
    /**
     * the maximum number of evaluations for the line search in {@link
     * #getIntersectionPoint(Line, Vector3D, Frame, AbsoluteDate)}.
     */
    private static final int MAX_EVALUATIONS = 100;

    /**
     * the default date to use when evaluating the {@link #harmonics}. Used when
     * no other dates are available. Should be removed in a future release.
     */
    private final AbsoluteDate defaultDate;
    /**
     * the reference ellipsoid.
     */
    private final ReferenceEllipsoid referenceEllipsoid;
    /**
     * the geo-potential combined with an algorithm for evaluating the spherical
     * harmonics. The Holmes and Featherstone method is very robust.
     */
    private final transient HolmesFeatherstoneAttractionModel harmonics;

    /**
     * Creates a geoid from the given geopotential, reference ellipsoid and the
     * assumptions in the comment for {@link Geoid}.
     *
     * @param geopotential       the gravity potential. Only the anomalous
     *                           potential will be used. It is assumed that the
     *                           {@code geopotential} and the {@code
     *                           referenceEllipsoid} are defined in the same
     *                           frame. Usually a {@link GravityFields#getConstantNormalizedProvider(int,
     *                           int, AbsoluteDate) constant geopotential} is used to define a
     *                           time-invariant Geoid.
     * @param referenceEllipsoid the normal gravity potential.
     * @throws NullPointerException if {@code geopotential == null ||
     *                              referenceEllipsoid == null}
     */
    public Geoid(final NormalizedSphericalHarmonicsProvider geopotential,
                 final ReferenceEllipsoid referenceEllipsoid) {
        // parameter check
        if (geopotential == null || referenceEllipsoid == null) {
            throw new NullPointerException();
        }

        // subtract the ellipsoid from the geopotential
        final SubtractEllipsoid potential = new SubtractEllipsoid(geopotential,
                referenceEllipsoid);

        // set instance parameters
        this.referenceEllipsoid = referenceEllipsoid;
        this.harmonics = new HolmesFeatherstoneAttractionModel(
                referenceEllipsoid.getBodyFrame(), potential);
        this.defaultDate = AbsoluteDate.ARBITRARY_EPOCH;
    }

    @Override
    public Frame getBodyFrame() {
        // same as for reference ellipsoid.
        return this.getEllipsoid().getBodyFrame();
    }

    /**
     * Gets the Undulation of the Geoid, N at the given position. N is the
     * distance between the {@link #getEllipsoid() reference ellipsoid} and the
     * geoid. The latitude and longitude parameters are both defined with
     * respect to the reference ellipsoid. For EGM96 and the WGS84 ellipsoid the
     * undulation is between -107m and +86m.
     *
     * <p> NOTE: Restrictions are not put on the range of the arguments {@code
     * geodeticLatitude} and {@code longitude}.
     *
     * @param geodeticLatitude geodetic latitude (angle between the local normal
     *                         and the equatorial plane on the reference
     *                         ellipsoid), in radians.
     * @param longitude        on the reference ellipsoid, in radians.
     * @param date             of evaluation. Used for time varying geopotential
     *                         fields.
     * @return the undulation in m, positive means the geoid is higher than the
     * ellipsoid.
     * @see Geoid
     * @see <a href="http://en.wikipedia.org/wiki/Geoid">Geoid on Wikipedia</a>
     */
    public double getUndulation(final double geodeticLatitude,
                                final double longitude,
                                final AbsoluteDate date) {
        /*
         * equations references are to the algorithm printed in the geoid
         * cookbook[2]. See comment for Geoid.
         */
        // reference ellipsoid
        final ReferenceEllipsoid ellipsoid = this.getEllipsoid();

        // position in geodetic coordinates
        final GeodeticPoint gp = new GeodeticPoint(geodeticLatitude, longitude, 0);
        // position in Cartesian coordinates, is converted to geocentric lat and
        // lon in the Holmes and Featherstone class
        final Vector3D position = ellipsoid.transform(gp);

        // get normal gravity from ellipsoid, eq 15
        final double normalGravity = ellipsoid
                .getNormalGravity(geodeticLatitude);

        // calculate disturbing potential, T, eq 30.
        final double mu = this.harmonics.getMu(date);
        final double T  = this.harmonics.nonCentralPart(date, position, mu);
        // calculate undulation, eq 30
        return T / normalGravity;
    }

    @Override
    public ReferenceEllipsoid getEllipsoid() {
        return this.referenceEllipsoid;
    }

    /**
     * This class implements equations 20-24 in the geoid cook book.(Losch and
     * Seufer) It modifies C<sub>2n,0</sub> where n = 1,2,...,5.
     *
     * @see "DMA TR 8350.2. 1984."
     */
    private static final class SubtractEllipsoid implements
            NormalizedSphericalHarmonicsProvider {
        /**
         * provider of the fully normalized coefficients, includes the reference
         * ellipsoid.
         */
        private final NormalizedSphericalHarmonicsProvider provider;
        /**
         * the reference ellipsoid to subtract from {@link #provider}.
         */
        private final ReferenceEllipsoid ellipsoid;

        /**
         * @param provider  potential used for GM<sub>g</sub> and a<sub>g</sub>,
         *                  and of course the coefficients Cnm, and Snm.
         * @param ellipsoid Used to calculate the fully normalized
         *                  J<sub>2n</sub>
         */
        private SubtractEllipsoid(
                final NormalizedSphericalHarmonicsProvider provider,
                final ReferenceEllipsoid ellipsoid) {
            super();
            this.provider = provider;
            this.ellipsoid = ellipsoid;
        }

        @Override
        public int getMaxDegree() {
            return this.provider.getMaxDegree();
        }

        @Override
        public int getMaxOrder() {
            return this.provider.getMaxOrder();
        }

        @Override
        public double getMu() {
            return this.provider.getMu();
        }

        @Override
        public double getAe() {
            return this.provider.getAe();
        }

        @Override
        public AbsoluteDate getReferenceDate() {
            return this.provider.getReferenceDate();
        }

        @Override
        public NormalizedSphericalHarmonics onDate(final AbsoluteDate date) {
            return new NormalizedSphericalHarmonics() {

                /** the original harmonics */
                private final NormalizedSphericalHarmonics delegate = provider.onDate(date);

                @Override
                public double getNormalizedCnm(final int n, final int m) {
                    return getCorrectedCnm(n, m, this.delegate.getNormalizedCnm(n, m));
                }

                @Override
                public double getNormalizedSnm(final int n, final int m) {
                    return this.delegate.getNormalizedSnm(n, m);
                }

                @Override
                public AbsoluteDate getDate() {
                    return date;
                }
            };
        }

        /**
         * Get the corrected Cnm for different GM or a values.
         *
         * @param n              degree
         * @param m              order
         * @param uncorrectedCnm uncorrected Cnm coefficient
         * @return the corrected Cnm coefficient.
         */
        private double getCorrectedCnm(final int n,
                                       final int m,
                                       final double uncorrectedCnm) {
            double Cnm = uncorrectedCnm;
            // n = 2,4,6,8, or 10 and m = 0
            if (m == 0 && n <= 10 && n % 2 == 0 && n > 0) {
                // correction factor for different GM and a, 1 if no difference
                final double gmRatio = this.ellipsoid.getGM() / this.getMu();
                final double aRatio = this.ellipsoid.getEquatorialRadius() /
                        this.getAe();
                /*
                 * eq 20 in the geoid cook book[2], with eq 3-61 in chapter 3 of
                 * DMA TR 8350.2
                 */
                // halfN = 1,2,3,4,5 for n = 2,4,6,8,10, respectively
                final int halfN = n / 2;
                Cnm = Cnm - gmRatio * FastMath.pow(aRatio, halfN) *
                        this.ellipsoid.getC2n0(halfN);
            }
            // return is a modified Cnm
            return Cnm;
        }

        @Override
        public TideSystem getTideSystem() {
            return this.provider.getTideSystem();
        }

    }

    /**
     * {@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) {
        /*
         * 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 StaticTransform frameToBody =
                frame.getStaticTransformTo(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;
        }
    }

    @Override
    public Vector3D projectToGround(final Vector3D point,
                                    final AbsoluteDate date,
                                    final Frame frame) {
        final GeodeticPoint gp = this.transform(point, frame, date);
        final GeodeticPoint gpZero =
                new GeodeticPoint(gp.getLatitude(), gp.getLongitude(), 0);
        final StaticTransform bodyToFrame =
                this.getBodyFrame().getStaticTransformTo(frame, date);
        return bodyToFrame.transformPosition(this.transform(gpZero));
    }

    /**
     * {@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 CalculusFieldElement<T>> FieldGeodeticPoint<T> getIntersectionPoint(final FieldLine<T> lineInFrame,
                                                                                      final FieldVector3D<T> closeInFrame,
                                                                                      final Frame frame,
                                                                                      final FieldAbsoluteDate<T> date) {

        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 FieldStaticTransform<T> frameToBody = frame.getStaticTransformTo(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 CalculusFieldUnivariateFunction<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;
        }
    }

    @Override
    public TimeStampedPVCoordinates projectToGround(
            final TimeStampedPVCoordinates pv,
            final Frame frame) {
        throw new UnsupportedOperationException();
    }

    /**
     * {@inheritDoc}
     *
     * @param date date of the conversion. Used for computing frame
     *             transformations and for time dependent geopotential.
     * @return The surface relative point at the same location. Altitude is
     * orthometric height, that is height above the {@link Geoid}. Latitude and
     * longitude are both geodetic and defined with respect to the {@link
     * #getEllipsoid() reference ellipsoid}.
     * @see #transform(GeodeticPoint)
     * @see <a href="http://en.wikipedia.org/wiki/Orthometric_height">Orthometric_height</a>
     */
    @Override
    public GeodeticPoint transform(final Vector3D point, final Frame frame,
                                   final AbsoluteDate date) {
        // convert using reference ellipsoid, altitude referenced to ellipsoid
        final GeodeticPoint ellipsoidal = this.getEllipsoid().transform(
                point, frame, date);
        // convert altitude to orthometric using the undulation.
        final double undulation = this.getUndulation(ellipsoidal.getLatitude(),
                ellipsoidal.getLongitude(), date);
        // add undulation to the altitude
        return new GeodeticPoint(
                ellipsoidal.getLatitude(),
                ellipsoidal.getLongitude(),
                ellipsoidal.getAltitude() - undulation
        );
    }

    /**
     * {@inheritDoc}
     *
     * @param date date of the conversion. Used for computing frame
     *             transformations and for time dependent geopotential.
     * @return The surface relative point at the same location. Altitude is
     * orthometric height, that is height above the {@link Geoid}. Latitude and
     * longitude are both geodetic and defined with respect to the {@link
     * #getEllipsoid() reference ellipsoid}.
     * @see #transform(GeodeticPoint)
     * @see <a href="http://en.wikipedia.org/wiki/Orthometric_height">Orthometric_height</a>
     */
    @Override
    public <T extends CalculusFieldElement<T>> FieldGeodeticPoint<T> transform(final FieldVector3D<T> point, final Frame frame,
                                                                           final FieldAbsoluteDate<T> date) {
        // convert using reference ellipsoid, altitude referenced to ellipsoid
        final FieldGeodeticPoint<T> ellipsoidal = this.getEllipsoid().transform(
                point, frame, date);
        // convert altitude to orthometric using the undulation.
        final double undulation = this.getUndulation(ellipsoidal.getLatitude().getReal(),
                                                     ellipsoidal.getLongitude().getReal(),
                                                     date.toAbsoluteDate());
        // add undulation to the altitude
        return new FieldGeodeticPoint<>(
                ellipsoidal.getLatitude(),
                ellipsoidal.getLongitude(),
                ellipsoidal.getAltitude().subtract(undulation)
        );
    }

    /**
     * {@inheritDoc}
     *
     * @param point The surface relative point to transform. Altitude is
     *              orthometric height, that is height above the {@link Geoid}.
     *              Latitude and longitude are both geodetic and defined with
     *              respect to the {@link #getEllipsoid() reference ellipsoid}.
     * @return point at the same location but as a Cartesian point in the {@link
     * #getBodyFrame() body frame}.
     * @see #transform(Vector3D, Frame, AbsoluteDate)
     */
    @Override
    public Vector3D transform(final GeodeticPoint point) {
        try {
            // convert orthometric height to height above ellipsoid using undulation
            // TODO pass in date to allow user to specify
            final double undulation = this.getUndulation(
                    point.getLatitude(),
                    point.getLongitude(),
                    this.defaultDate
            );
            final GeodeticPoint ellipsoidal = new GeodeticPoint(
                    point.getLatitude(),
                    point.getLongitude(),
                    point.getAltitude() + undulation
            );
            // transform using reference ellipsoid
            return this.getEllipsoid().transform(ellipsoidal);
        } catch (OrekitException e) {
            //this method, as defined in BodyShape, is not permitted to throw
            //an OrekitException, so wrap in an exception we can throw.
            throw new RuntimeException(e);
        }
    }

    /**
     * {@inheritDoc}
     *
     * @param point The surface relative point to transform. Altitude is
     *              orthometric height, that is height above the {@link Geoid}.
     *              Latitude and longitude are both geodetic and defined with
     *              respect to the {@link #getEllipsoid() reference ellipsoid}.
     * @param <T> type of the field elements
     * @return point at the same location but as a Cartesian point in the {@link
     * #getBodyFrame() body frame}.
     * @see #transform(Vector3D, Frame, AbsoluteDate)
     * @since 9.0
     */
    @Override
    public <T extends CalculusFieldElement<T>> FieldVector3D<T> transform(final FieldGeodeticPoint<T> point) {
        try {
            // convert orthometric height to height above ellipsoid using undulation
            // TODO pass in date to allow user to specify
            final double undulation = this.getUndulation(
                    point.getLatitude().getReal(),
                    point.getLongitude().getReal(),
                    this.defaultDate
            );
            final FieldGeodeticPoint<T> ellipsoidal = new FieldGeodeticPoint<>(
                    point.getLatitude(),
                    point.getLongitude(),
                    point.getAltitude().add(undulation)
            );
            // transform using reference ellipsoid
            return this.getEllipsoid().transform(ellipsoidal);
        } catch (OrekitException e) {
            //this method, as defined in BodyShape, is not permitted to throw
            //an OrekitException, so wrap in an exception we can throw.
            throw new RuntimeException(e);
        }
    }

}