Procedure

     NEARPT ( Nearest point on an ellipsoid )

     SUBROUTINE NEARPT ( POSITN, A, B, C, NPOINT, ALT )

Abstract

     Locate the point on the surface of an ellipsoid that is nearest
     to a specified position. Also return the altitude of the position
     above the ellipsoid.

Required_Reading

     None.

Keywords

     ALTITUDE
     ELLIPSOID
     GEOMETRY

Declarations

     IMPLICIT NONE

     DOUBLE PRECISION      POSITN ( 3 )
     DOUBLE PRECISION      A
     DOUBLE PRECISION      B
     DOUBLE PRECISION      C
     DOUBLE PRECISION      NPOINT ( 3 )
     DOUBLE PRECISION      ALT

Brief_I/O

     VARIABLE  I/O  DESCRIPTION
     --------  ---  --------------------------------------------------
     POSITN     I   Position of a point in body-fixed frame.
     A          I   Length of semi-axis parallel to x-axis.
     B          I   Length of semi-axis parallel to y-axis.
     C          I   Length on semi-axis parallel to z-axis.
     NPOINT     O   Point on the ellipsoid closest to POSITN.
     ALT        O   Altitude of POSITN above the ellipsoid.

Detailed_Input

     POSITN   is a 3-vector giving the position of a point with respect
              to the center of an ellipsoid. The vector is expressed
              in a body-fixed reference frame. The semi-axes of the
              ellipsoid are aligned with the x, y, and z-axes of the
              body-fixed frame.

     A        is the length of the semi-axis of the ellipsoid that is
              parallel to the x-axis of the body-fixed reference frame.

     B        is the length of the semi-axis of the ellipsoid that is
              parallel to the y-axis of the body-fixed reference frame.

     C        is the length of the semi-axis of the ellipsoid that is
              parallel to the z-axis of the body-fixed reference frame.

Detailed_Output

     NPOINT   is the nearest point on the ellipsoid to POSITN.
              NPOINT is a 3-vector expressed in the body-fixed
              reference frame.

     ALT      is the altitude of POSITN above the ellipsoid. If
              POSITN is inside the ellipsoid, ALT will be negative
              and have magnitude equal to the distance between
              NPOINT and POSITN.

Parameters

     None.

Exceptions

     1)  If any of the axis lengths A, B or C are non-positive, the
         error SPICE(BADAXISLENGTH) is signaled.

     2)  If the ratio of the longest to the shortest ellipsoid axis
         is large enough so that arithmetic expressions involving its
         squared value may overflow, the error SPICE(BADAXISLENGTH)
         is signaled.

     3)  If any of the expressions

            A * ABS( POSITN(1) ) / m**2
            B * ABS( POSITN(2) ) / m**2
            C * ABS( POSITN(3) ) / m**2

         where `m' is the minimum of { A, B, C }, is large enough so
         that arithmetic expressions involving these sub-expressions
         may overflow, the error SPICE(INPUTSTOOLARGE) is signaled.

     4)  If the axes of the ellipsoid have radically different
         magnitudes, for example if the ratios of the axis lengths vary
         by 10 orders of magnitude, the results may have poor
         precision. No error checks are done to identify this problem.

     5)  If the axes of the ellipsoid and the input point POSITN have
         radically different magnitudes, for example if the ratio of
         the magnitude of POSITN to the length of the shortest axis is
         1.E25, the results may have poor precision. No error checks
         are done to identify this problem.

Files

     None.

Particulars

     Many applications of this routine are more easily performed
     using the higher-level SPICELIB routine SUBPNT. This routine
     is the mathematical workhorse on which SUBPNT relies.

     This routine solves for the location, N, on the surface of an
     ellipsoid nearest to an arbitrary location, P, relative to that
     ellipsoid.

Examples

     Example 1.

     The code fragment below illustrates how you can use SPICELIB to
     compute the apparent sub-earth point on the moon.

     C
     C     Load the ephemeris, leapseconds and physical constants
     C     files first. We assume the names of these files are
     C     stored in the character variables SPK, LSK and
     C     PCK.
     C
           CALL FURNSH ( SPK )
           CALL FURNSH ( LSK )
           CALL FURNSH ( PCK )

     C
     C     Get the apparent position of the moon as seen from the
     C     earth. Look up this position vector in the moon
     C     body-fixed frame IAU_MOON. The orientation of the
     C     IAU_MOON frame will be computed at epoch ET-LT.
     C
           CALL SPKPOS ( 'moon',  ET,    'IAU_MOON', 'LT+S',
          .              'earth', TRGPOS, LT                 )

     C
     C     Negate the moon's apparent position to obtain the
     C     position of the earth in the moon's body-fixed frame.
     C
           CALL VMINUS ( TRGPOS, EVEC )

     C
     C     Get the lengths of the principal axes of the moon.
     C     Transfer the elements of the array RADII to the
     C     variables A, B, C to enhance readability.
     C
           CALL BODVRD (  'MOON',    'RADII', DIM, RADII )
           CALL VUPACK (  RADII,  A,       B,   C     )

     C
     C     Finally get the point SUBPNT on the surface of the
     C     moon closest to the earth --- the sub-earth point.
     C     SUBPNT is expressed in the IAU_MOON reference frame.
     C
           CALL NEARPT ( EVEC, A, B, C, SUBPNT, ALT )


     Example 2.

     One can use this routine to define a generalization of GEODETIC
     coordinates called GAUSSIAN coordinates of a triaxial body. (The
     name is derived from the famous Gauss-map of classical
     differential geometry).  The coordinates are longitude,
     latitude, and altitude.

     We let the x-axis of the body fixed coordinate system point
     along the longest axis of the triaxial body. The y-axis points
     along the middle axis and the z-axis points along the shortest
     axis.

     Given a point P, there is a point on the ellipsoid that is
     closest to P, call it Q. The latitude and longitude of P
     are determined by constructing the outward pointing unit normal
     to the ellipsoid at Q. Latitude of P is the latitude that the
     normal points toward in the body-fixed frame. Longitude is
     the longitude the normal points to in the body-fixed frame.
     The altitude is the signed distance from P to Q, positive if P
     is outside the ellipsoid, negative if P is inside.
     (the mapping of the point Q to the unit normal at Q is the
     Gauss-map of Q).

     To obtain the Gaussian coordinates of a point whose position
     in body-fixed rectangular coordinates is given by a vector P,
     the code fragment below will suffice.

        CALL NEARPT ( P,    A, B, C, Q, ALT  )
        CALL SURFNM (       A, B, C  Q, NRML )

        CALL RECLAT ( NRML, R, LONG, LAT     )

     The Gaussian coordinates are LONG, LAT, and ALT.

Restrictions

     See $Exceptions section.

Literature_References

     None.

Author_and_Institution

     C.H. Acton         (JPL)
     N.J. Bachman       (JPL)
     J. Diaz del Rio    (ODC Space)
     W.L. Taber         (JPL)
     E.D. Wright        (JPL)

Version

    SPICELIB Version 2.0.0, 26-OCT-2021 (NJB) (JDR) (EDW)

        Edit to logic to reduce unneeded operations when
        error or projection vectors equal zero. Addition
        of details concerning the "ellipsoid near point"
        problem and solution.

        Added IMPLICIT NONE statement.

        Edited the header to comply with NAIF standard.

    SPICELIB Version 1.4.0, 27-JUN-2013 (NJB)

        Updated in-line comments.

     Last update was 04-MAR-2013 (NJB)

        Bug fix: now correctly computes off-axis solution for
        the case of a prolate ellipsoid and a viewing point
        on the interior long axis.

    SPICELIB Version 1.3.1, 07-FEB-2008 (NJB)

        Header update: header now refers to SUBPNT rather
        than deprecated routine SUBPT.

    SPICELIB Version 1.3.0, 07-AUG-2006 (NJB)

        Bug fix: added initialization of variable SNGLPT to support
                  operation under the Macintosh Intel Fortran
                  compiler. Note that this bug did not affect
                  operation of this routine on other platforms.

    SPICELIB Version 1.2.0, 15-NOV-2005 (EDW) (NJB)

        Various changes were made to ensure that all loops terminate.

        Bug fix: scale of transverse component of error vector
        was corrected for the exterior point case.

        Bug fix: non-standard use of duplicate arguments in VSCL
        calls was corrected.

        Error checking was added to screen out inputs that might
        cause numeric overflow.

        Replaced BODVAR call in examples to BODVRD.

    SPICELIB Version 1.1.1, 28-JUL-2003 (NJB) (CHA)

        Various header corrections were made.

    SPICELIB Version 1.1.0, 27-NOV-1990 (WLT)

        The routine was substantially rewritten to achieve
        more robust behavior and document the mathematics
        of the routine.

    SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)

        Comment section for permuted index source lines was added
        following the header.

    SPICELIB Version 1.0.0, 31-JAN-1990 (WLT)