Ellipses Required Reading

Index Page Ellipses Required Reading
## Table of Contents
       Ellipses Required Reading
          Abstract
          Note on FORTRAN and C Versions
          Introduction
          References
       The ellipse data type
       Ellipse and ellipsoid routines
          Ellipse access routines
             Creating a SPICELIB ellipse
             `Breaking apart' a SPICELIB ellipse
             CGV2EL and EL2CGV are not inverses
          Triaxial ellipsoid routines
             Finding the nearest point on an ellipsoid to a point
             Finding the intersection of a ray and an ellipsoid
             Finding the outward normal at a surface point
             Finding the limb of an ellipsoid
             Finding the nearest point on an ellipsoid to a line
             Finding the intersection of an ellipsoid and a plane
          Ellipse routines
             Finding the intersection of an ellipse and plane
             Finding the nearest point on an ellipse to a point
             Finding the semi-axes of an ellipse from generating vectors
             Finding the orthogonal projection of an ellipse onto a plane
       Examples
             Finding the sub-observer point on an extended body
             Testing for occultation of a star
             Finding the `limb angle' of an instrument boresight
             Finding the instrument boresight---limb distance
       Mathematical notes
          Defining an ellipse parametrically
          Solving intersection problems
       Appendix: Document Revision History
             December 21, 2004
             December 12, 2002.
Top # Ellipses Required Reading Last revised on 2008 JAN 17 by B. V. Semenov. Top ## Abstract SPICELIB contains a substantial set of subroutines that solve common mathematical problems involving ellipses and triaxial ellipsoids. Top ## Note on FORTRAN and C Versions This document covers the FORTRAN version of the interfaces of this subsystem. CSPICE provides f2c translated equivalents for all, and native C wrappers for some of them. If you wish to use the C versions of the interfaces described in this document, refer to the CSPICE Required Reading, cspice.req, for more information on naming conventions, locations, and usage of the f2c'ed routines and native C wrappers. Top ## Introduction Ellipses turn up frequently in the sort of science analysis problems SPICELIB is designed to help solve. The shapes of extended bodies--planets, satellites, and the Sun--are frequently modelled by triaxial ellipsoids. The IAU has defined such models for the Sun, all of the planets, and most of their satellites, in the 1988 IAU/IAG/COSPAR working group report [1]. Many geometry problems involving triaxial ellipsoids give rise to ellipses as `mathematical byproducts'. Ellipses are also used in modelling orbits and planetary rings. SPICELIB contains a substantial set of subroutines that solve common mathematical problems involving ellipses and triaxial ellipsoids. This required reading file documents those routines, gives examples of their use, and presents some of the mathematical background required to understand the routines. This required reading file is divided into four chapters: `--` The ellipse data type `--` Ellipse and ellipsoid routines `--` Examples `--` Mathematical notes The first chapter discusses a structured data type in SPICELIB called `ellipses'. The first chapter is a prerequisite for the next two, since many SPICELIB subroutines make use of the ellipse data type in their calling sequences. The second chapter lists the SPICELIB routines that deal with ellipses and ellipsoids. The third chapter presents some extended examples that illustrate the use of SPICELIB ellipse and ellipsoid routines to solve practical geometry problems. The fourth chapter presents some of the mathematical background necessary to fully understand the SPICELIB ellipse and ellipsoid routines. Top ## References [1] `Report of the IAU/IAG/COSPAR Working Group on Cartographic Coordinates and Rotational Elements of the Planets and Satellites: 1988', January 4, 1989. [2] `Calculus, Vol. II'. Tom Apostol. John Wiley and Sons, 1969. See Chapter 5, `Eigenvalues of Operators Acting on Euclidean Spaces'. [3] PLANES required reading (planes.req). Top # The ellipse data type The `ellipse' is a structured data type used in SPICELIB to represent ellipses in three-dimensional space. SPICELIB ellipses are used to simplify calling sequences of routines that output or accept as input data that defines ellipses. We will not dwell on the argument justifying this choice; see the introduction of the PLANES required reading, planes.req, for a discussion of the topic. SPICELIB ellipses are declared as double precision arrays of length UBEL; the value of the parameter UBEL is 9. As a matter of good programming practice, we recommend that you declare ellipses using the parameter UBEL, as in the example below.
       INTEGER               UBEL
       PARAMETER           ( UBEL = 9 )

       DOUBLE PRECISION      ELLIPS ( UBEL )
The values of SPICELIB ellipses are set and retrieved via two access routines, CGV2EL (Center and generating vectors to ellipse) and EL2CGV (Ellipse to center and generating vectors). You should never directly assign values to the elements of a SPICELIB ellipse or make an assignment from an element of a SPICELIB ellipse. The contents of SPICELIB ellipses are considered an implementation choice, and are not described by the specification of any SPICELIB routine. In the future, the definition of the contents of SPICELIB ellipses may be changed in any way that does not affect the user-visible behavior of the access routines. The following representation of an ellipse is used throughout SPICELIB, and in particular by the ellipse access routines: An ellipse is the set of points
       CENTER    +    cos(theta) * V1    +    sin(theta) * V2
where CENTER, V1, and V2 are three vectors, and theta is in the range
       (-pi, pi].
It is a fact, which we will prove in the chapter titled `Mathematical notes', that this set of points is an ellipse. The ellipse defined by this parametric representation is non-degenerate if and only if V1 and V2 are linearly independent. We call CENTER the `center' of the ellipse, and we refer to V1 and V2 as `generating vectors'. Note that an ellipse centered at the origin is completely specified by its generating vectors. In the rest of this document, when we refer to the center or generating vectors for a particular ellipse, we mean vectors that play the role of CENTER or V1 and V2 in defining that ellipse. This representation of ellipses has the particularly convenient property that it allows easy computation of the image of an ellipse under a linear transformation. If M is a matrix representing a linear transformation, and E is the ellipse
       CENTER    +    cos(theta) * V1    +    sin(theta) * V2,
then the image of E under the transformation represented by M is
       M*CENTER    +    cos(theta) * M*V1    +    sin(theta) * M*V2.
If we accept that the first set of points is an ellipse, then we can see that the image of an ellipse under a linear transformation is always another (possibly degenerate) ellipse. Since many geometric computations involving ellipses and ellipsoids may be greatly simplified by judicious application of linear transformations to ellipses, it is useful to have a representation for ellipses that allows ready computation of their images under such mappings. Top # Ellipse and ellipsoid routines The SPICELIB routines that deal with ellipses and ellipsoids are: CGV2EL ( Center and generating vectors to ellipse ) EL2CGV ( Ellipse to center and generating vectors ) NEARPT ( Nearest point on ellipsoid to point ) SURFPT ( Surface intercept point on ellipsoid ) SURFNM ( Surface normal on ellipsoid ) EDLIMB ( Ellipsoid limb ) NPEDLN ( Nearest point on ellipsoid to line ) INEDPL ( Intersection of ellipsoid and plane) INELPL ( Intersection of ellipse and plane ) NPELPT ( Nearest point on ellipse to point ) PJELPL ( Projection of ellipse onto plane ) SAELGV ( Semi-axes of ellipse from generating vectors ) In the code examples that follow, we will repeatedly use the variable names CENTER, V1, V2, and ELLIPS. Wherever these appear, they refer to the center and a pair of generating vectors for an ellipse, and a SPICELIB ellipse. These variables can always be assumed to have been declared as follows:
       INTEGER               UBEL
       PARAMETER           ( UBEL = 9 )

       DOUBLE PRECISION      ELLIPS ( UBEL )
       DOUBLE PRECISION      CENTER ( 3 )
       DOUBLE PRECISION      V1     ( 3 )
       DOUBLE PRECISION      V2     ( 3 )
Top ## Ellipse access routines The SPICELIB routines used for accessing ellipses are CGV2EL ( Center and generating vectors to ellipse ) EL2CGV ( Ellipse to center and generating vectors ) Top ### Creating a SPICELIB ellipse Let CENTER, V1, and V2 be a center vector and two generating vectors for an ellipse. To create a SPICELIB ellipse from these vectors, we could make the declarations shown above for CENTER, V1, and V2, and use the subroutine call
       CALL CGV2EL ( CENTER, V1, V2, ELLIPS )
This call produces the SPICELIB ellipse ELLIPS, which represents the same mathematical ellipse as do CENTER, V1, and V2. The generating vectors need not be linearly independent. If they are not, the resulting ellipse will be degenerate. Specifically, if the generating vectors are both zero, the ellipse will be the single point represented by CENTER, and if just one of the semi-axis vectors (call it V) is non-zero, the ellipse will be the line segment extending from
       CENTER - V
to
       CENTER + V
Top ### `Breaking apart' a SPICELIB ellipse Let ELLIPS be a SPICELIB ellipse. To produce the center and two generating vectors for ELLIPS, we can make the call
       CALL EL2CGV ( ELLIPS, CENTER, V1, V2 )
On output, V1 will be a semi-major axis vector for the ellipse represented by ELLIPS, and V2 will be a semi-minor axis vector. Semi-axis vectors are never unique; if X is a semi-axis vector; then so is -X. V1 is a vector of maximum norm extending from the ellipse's center to the ellipse itself; V2 is an analogous vector of minimum norm. V1 and V2 are orthogonal vectors. Top ### CGV2EL and EL2CGV are not inverses Because the routine EL2CGV always returns semi-axes as generating vectors, if V1 and V2 are not semi-axes on input to CGV2EL, the sequence of calls
       CALL CGV2EL ( CENTER,  V1,      V2,  ELLIPS )
       CALL EL2CGV ( ELLIPS,  CENTER,  V1,  V2     )
will certainly modify V1 and V2. Even if V1 and V2 are semi-axes to start out with, because of the non-uniqueness of semi-axes, one or both of these vectors could be negated on output from EL2CGV. There is a sense in which CGV2EL and EL2CGV are inverses, though: the above sequence of calls returns a center and generating vectors that define the same ellipse as the input center and generating vectors. Top ## Triaxial ellipsoid routines The SPICELIB routines that deal with triaxial ellipsoids are NEARPT ( Nearest point on ellipsoid to point ) SURFPT ( Surface intercept point on ellipsoid ) SURFNM ( Surface normal on ellipsoid ) EDLIMB ( Ellipsoid limb ) NPEDLN ( Nearest point on ellipsoid to line ) INEDPL ( Intersection of ellipsoid and plane) In the rest of this section, the variables A, B, and C will be used to denote the semi-axis lengths of a triaxial ellipsoid. The declarations of these variables are
       DOUBLE PRECISION      A
       DOUBLE PRECISION      B
       DOUBLE PRECISION      C
Top ### Finding the nearest point on an ellipsoid to a point The routine NEARPT finds the nearest point on a triaxial ellipsoid to a specified point. The point may be outside, on, or inside the ellipsoid. If the variable POINT and output arguments NRPT and ALT are declared by
       DOUBLE PRECISION      POINT ( 3 )
       DOUBLE PRECISION      NRPT  ( 3 )
       DOUBLE PRECISION      ALT
then the code fragment
       CALL NEARPT ( POINT, A, B, C, NRPT, ALT )
will return NRPT as the nearest point on the ellipsoid to POINT and ALT as the altitude of POINT above the ellipsoid. If POINT is inside the ellipsoid, ALT will be negative. NEARPT is useful for, among other things, finding the sub-spacecraft point or sub-solar point on a target body. It also is used as a utility routine by the SPICELIB routines that convert between geodetic and other coordinate systems. Top ### Finding the intersection of a ray and an ellipsoid The routine SURFPT finds the intersection of a ray and a triaxial ellipsoid. Let a ray be defined by two variables, VERTEX and DIR, that represent the vertex and direction vector. Let XPT represent the intersection point and let FOUND be a logical flag indicating whether the ray and ellipsoid intersect. VERTEX, DIR, XPT, and FOUND are declared as follows:
       DOUBLE PRECISION      VERTEX ( 3 )
       DOUBLE PRECISION      DIR    ( 3 )
       DOUBLE PRECISION      XPT    ( 3 )
       LOGICAL               FOUND
Then if A, B, and C are the semi-axis lengths of the ellipsoid, the call
       CALL SURFPT ( VERTEX, DIR, A, B, C, XPT, FOUND )
returns the location of the intercept point in XPT and sets FOUND to .TRUE. if the ray does intersect the ellipsoid. Otherwise, FOUND is set to .FALSE. and XPT is set to the zero vector. A typical use of SURFPT would be finding the surface intercept of an instrument boresight ray on an extended body. SURFPT can also be used to for a low-accuracy determination of whether a star is occulted by an extended body. Such a test would ignore factors such as topography of the body, differential stellar aberration, differential light-time corrections, and light-bending. Top ### Finding the outward normal at a surface point The routine SURFNM finds the unit outward normal vector at a specified point on a triaxial ellipsoid. Let A, B, C be the lengths of the semi-axes of the ellipsoid, let the variable POINT represent a surface point on the ellipsoid, and let NORMAL represent the unit outward normal at POINT. The variables POINT and NORMAL are declared
       DOUBLE PRECISION      POINT  ( 3 )
       DOUBLE PRECISION      NORMAL ( 3 )
The call
       CALL SURFNM ( A, B, C, POINT, NORMAL )
returns NORMAL as the unit outward normal vector at the surface point represented by POINT. SURFNM is useful in any application requiring a local surface coordinate system with one coordinate plane tangent to the surface at a specified point. It is also useful for computing emission and solar incidence angles at a surface point. SURFNM is also used as a utility by SPICELIB routines that convert between geodetic and other coordinate systems. Top ### Finding the limb of an ellipsoid The routine EDLIMB finds the limb of a triaxial ellipsoid as viewed from a specified location. (The limb is the boundary of the visible portion of the ellipsoid; this boundary is an ellipse as long as the viewing location is outside of the ellipsoid.) Let A, B, C be the lengths of the semi-axes of the ellipsoid, and let VIEWPT be the viewing location. The output argument LIMB is declared as a SPICELIB ellipse. We also declare variables that represent the center and semi-axes of the limb.
       INTEGER               UBEL
       PARAMETER           ( UBEL = 9 )

       DOUBLE PRECISION      VIEWPT ( 3 )
       DOUBLE PRECISION      LIMB   ( UBEL )

       DOUBLE PRECISION      CENTER ( 3 )
       DOUBLE PRECISION      SMAJOR ( 3 )
       DOUBLE PRECISION      SMINOR ( 3 )
After the call
       CALL EDLIMB ( A, B, C, VIEWPT, LIMB )
LIMB represents the limb of the ellipsoid as viewed from VIEWPT. The routine EL2CGV can be used to extract the center and semi-axes of the limb:
       CALL EL2CGV ( LIMB, CENTER, SMAJOR, SMINOR )
Several uses of the calculated limb are: `--` Predicting the appearance of an extended body in an instrument's field of view `--` Creating graphical overlays for limb-fitting `--` Finding the terminator of a body (using the vertex of the shadow cone as the viewing point) `--` Determining whether one ellipsoidal body occults another Top ### Finding the nearest point on an ellipsoid to a line The routine NPEDLN finds the nearest point on a triaxial ellipsoid to a specified line. Let A, B, C be the lengths of the semi-axes of the ellipsoid, and let the line be specified by a point and a direction vector represented by the variables POINT and DIR respectively. NPEDLN returns the nearest point on the ellipsoid to the line and the distance between that point and the line. The variables NRPT and DIST serve as these arguments.
       DOUBLE PRECISION      POINT  ( 3 )
       DOUBLE PRECISION      DIR    ( 3 )
       DOUBLE PRECISION      NRPT   ( 3 )
       DOUBLE PRECISION      DIST
The call
       CALL NPEDLN ( A, B, C, POINT, DIR, NRPT, DIST )
sets the values of NRPT and DIST as described. The capability of finding the nearest surface point to a line is useful in constructing maps that associate the nearest surface point to the boresight of a remote-sensing instrument at a given epoch with a measurement made by that instrument at that epoch. Note that this capability is NOT useful for determining whether the boresight of a remote sensing instrument intersects the surface of a body when topography is taken into account. Top ### Finding the intersection of an ellipsoid and a plane The routine INEDPL finds the intersection of a triaxial ellipsoid and a plane. Let A, B, C be the lengths of the semi-axes of the ellipsoid, and let PLANE and ELLIPS be declared as a SPICELIB plane and SPICELIB ellipse respectively:
       INTEGER               UBEL
       PARAMETER           ( UBEL = 9 )

       INTEGER               UBPL
       PARAMETER           ( UBPL = 4 )

       DOUBLE PRECISION      PLANE  ( UBPL )
       DOUBLE PRECISION      ELLIPS ( UBEL )
       LOGICAL               FOUND
Then the call
       CALL INEDPL ( A, B, C, PLANE, ELLIPS, FOUND )
returns ELLIPS as the ellipse of intersection of the ellipsoid and plane, if the intersection is non-empty. FOUND indicates whether the plane and ellipsoid intersect. INEDPL is used as a utility routine by the SPICELIB routines EDLIMB and NPEDLN. Other uses include determining whether an ellipsoidal body is partially within the field of view of a remote sensing instrument when the boundary of the field of view is polygonal, and finding the curve on the surface of an ellipsoid that corresponds to a line segment in an image. Top ## Ellipse routines The SPICELIB routines that deal with ellipses and do not fall into any of the above categories are: INELPL ( Intersection of ellipse and plane ) NPELPT ( Nearest point on ellipse to point ) PJELPL ( Projection of ellipse onto plane ) SAELGV ( Semi-axes of ellipse from generating vectors ) Top ### Finding the intersection of an ellipse and plane The routine INELPL finds the intersection of an ellipse and plane. Let ELLIPS and PLANE be declared as a SPICELIB ellipse and SPICELIB plane, respectively:
       INTEGER               UBEL
       PARAMETER           ( UBEL = 9 )

       INTEGER               UBPL
       PARAMETER           ( UBPL = 4 )

       DOUBLE PRECISION      ELLIPS ( UBEL )
       DOUBLE PRECISION      PLANE  ( UBPL )
       DOUBLE PRECISION      XPT1   ( 3 )
       DOUBLE PRECISION      XPT2   ( 3 )
       INTEGER               NXPTS
Then the call
       CALL INELPL ( ELLIPS, PLANE, NXPTS, XPT1, XPT2 )
sets the output argument NXPTS equal to the number of intersection points of the ELLIPS and PLANE. The possible values of NXPTS are -1, 0, 1, and 2. The value -1 is returned when the input ellipse lies in the input plane; the number of intersection points is infinite in that case. When NXPTS is equal to 1, both XPT1 and XPT2 are set to the value of the single point of intersection. When NXPTS is equal to 2, XPT1 and XPT2 are the points of intersection. Among the uses of INELPL are: `--` Finding the boundary of the illuminated portion of the limb of a body (the intersection of the terminator plane and the limb) `--` Finding the intersection of the limb of an ellipsoidal body with a plane defined by the focal point of a camera and an edge of the field of view of the camera--this calculation is useful for determining whether any part of the body appears in the field of view `--` Finding the limb points that lie in the plane defined by the body-centered position and velocity vectors of a spacecraft orbiting a body Top ### Finding the nearest point on an ellipse to a point The routine NPELPT finds the nearest point on an ellipse to a specified point. Let ELLIPS be declared as a SPICELIB ellipse; POINT is just a point in three-dimensional space:
       INTEGER               UBEL
       PARAMETER           ( UBEL = 9 )

       DOUBLE PRECISION      POINT  ( 3 )
       DOUBLE PRECISION      ELLIPS ( UBEL )
       DOUBLE PRECISION      NRPT   ( 3 )
       DOUBLE PRECISION      DIST
Then the call
       CALL NPELPT ( POINT, ELLIPS, NRPT, DIST )
returns NRPT as the point on ELLIPS closest to POINT, and DIST as the distance between NRPT and POINT. NPELPT is used as a utility routine by the SPICELIB routine NPEDLN. It is also useful for testing ellipses for equality: given a parametric representation of one ellipse, a number of points on that ellipse can be generated, and the distance of each from a second ellipse can be found. Top ### Finding the semi-axes of an ellipse from generating vectors The routine SAELGV finds the semi-axes of an ellipse in three-dimensional space, given two generating vectors. The generating vectors are not required to be linearly independent. If V1, V2, SMAJOR, and SMINOR are declared as follows,
       DOUBLE PRECISION      V1     ( 3 )
       DOUBLE PRECISION      V2     ( 3 )
       DOUBLE PRECISION      SMAJOR ( 3 )
       DOUBLE PRECISION      SMINOR ( 3 )
then the call
       CALL SAELGV ( V1, V2, SMAJOR, SMINOR )
returns SMAJOR as a semi-major axis of the ellipse and SMINOR as a semi-minor axis. SMAJOR and SMINOR are orthogonal; SMAJOR is a point on the ellipse of largest norm, and SMINOR is a point of smallest norm. The semi-axis vectors are not unique: the inverse of any semi-axis vector is also a semi-axis vector. SAELGV is used within SPICELIB as a utility routine. Top ### Finding the orthogonal projection of an ellipse onto a plane The routine PJELPL finds the orthogonal projection of an ellipse onto a specified plane. The orthogonal projection of an ellipse is another, possibly degenerate, ellipse. The projection is just the union of the orthogonal projections onto the plane of each point in the ellipse. The orthogonal projection of a point X onto a plane is the closest point in the plane to X. Let ELLIN and ELLOUT be declared as SPICELIB ellipses, and let PLANE be declared as a SPICELIB plane:
       INTEGER               UBEL
       PARAMETER           ( UBEL = 9 )

       INTEGER               UBPL
       PARAMETER           ( UBPL = 4 )

       DOUBLE PRECISION      ELLIN  ( UBEL )
       DOUBLE PRECISION      PLANE  ( UBPL )
       DOUBLE PRECISION      ELLOUT ( UBEL )
Then the call
       CALL PJELPL ( ELLIN, PLANE, ELLOUT )
returns the orthogonal projection of ELLIN onto PLANE as the SPICELIB ellipse ELLOUT. PJELPL is used as a utility routine in the SPICELIB routine INEDPL. Another application is finding the altitude of an instrument boresight above the limb of an ellipsoidal body; this procedure is demonstrated in the `Examples' chapter below. Top # Examples Top ### Finding the sub-observer point on an extended body In this example, we present a small program that solves a realistic science analysis geometry problem: finding the sub-observer point on an extended body at a specified time. We will assume that we have an SPK file available that contains ephemeris data for observer and target body at the time of interest. If you are not familar with SPICE kernel files, consult the SPK.REQ, TIME.REQ, and KERNEL.REQ required reading.
             PROGRAM SUBOBS
       C
       C     Compute the planetocentric latitude and longitude of the
       C     sub-observer point on an extended target body at a
       C     specified UTC epoch.  The observer may be a spacecraft,
       C     planet, or any body for which ephemeris data is available
       C     in an SPK file.
       C

       C
       C     SPICELIB functions
       C
             DOUBLE PRECISION      DPR

       C
       C     Local parameters
       C
             INTEGER               FILEN
             PARAMETER           ( FILEN = 128 )

       C
       C     Local variables
       C
             CHARACTER*(FILEN)     LEAP
             CHARACTER*(FILEN)     PCK
             CHARACTER*(FILEN)     SPK
             CHARACTER*(30)        UTC
             CHARACTER*(5)         CORR

             DOUBLE PRECISION      ALT
             DOUBLE PRECISION      DIST
             DOUBLE PRECISION      EPOCH
             DOUBLE PRECISION      ET
             DOUBLE PRECISION      LT
             DOUBLE PRECISION      OBSPOS ( 3 )
             DOUBLE PRECISION      RADII  ( 3 )
             DOUBLE PRECISION      STATE  ( 6 )
             DOUBLE PRECISION      SUBLAT
             DOUBLE PRECISION      SUBLON
             DOUBLE PRECISION      SUBPT  ( 3 )
             DOUBLE PRECISION      TIPM   ( 3, 3 )

             INTEGER               HANDLE
             INTEGER               N
             INTEGER               OBSRVR
             INTEGER               TARGET

             LOGICAL               CONT

             DATA CONT / .TRUE. /

       C
       C     Start out by prompting for the names of kernel files.
       C     Load each kernel as the name is supplied.
       C
             WRITE (*,*) ' '
             WRITE (*,*) 'Enter name of SPK file'
             READ  (*,FMT='(A)') SPK

             CALL SPKLEF ( SPK, HANDLE )

             WRITE (*,*) 'Enter name of leapseconds kernel'
             READ  (*,FMT='(A)') LEAP

             CALL FURNSH ( LEAP )

             WRITE (*,*) 'Enter name of P_constants kernel'
             READ  (*,FMT='(A)') PCK

             CALL FURNSH ( PCK )

             DO WHILE ( CONT )

                WRITE (*,*) ' '
                WRITE (*,*) 'Enter NAIF integer code of target body'
                READ   *,    TARGET

                WRITE (*,*) 'Enter NAIF integer code of observer body'
                READ   *,    OBSRVR

                WRITE (*,*) 'Enter UTC epoch of observation'
                READ  (*,FMT='(A)') UTC

                WRITE (*,*) 'Enter aberration correction to use:'
                WRITE (*,*) 'NONE, LT, or LT+S'
                READ  (*,FMT='(A)') CORR

       C
       C        Convert the UTC epoch to ET.
       C
                CALL UTC2ET ( UTC, ET )

       C
       C        Find the state of the target as seen from the observer at
       C        ET.
       C
                CALL SPKEZ ( TARGET, ET, 'J2000', CORR, OBSRVR, STATE,
            .                LT                                        )

       C
       C        The pure mathematical problem embedded in this application
       C        is finding the nearest point on a triaxial ellipsoid to a
       C        specified point.  The SPICELIB routine NEARPT solves this
       C        problem.  NEARPT deals with an ellipsoid that is centered
       C        at the origin and whose axes are lined up with the x, y,
       C        and z coordinate axes.  So to use NEARPT, we'll have to
       C        convert the problem at hand into a form that NEARPT can
       C        accept.
       C
       C        The way to do this is to convert all of our vectors into
       C        body equator and prime meridian coordinates.  The body
       C        equator and prime meridian system is one in which the
       C        coordinate axes are lined up with the body's geometric
       C        axes (there may be a few exceptions to the rule--check the
       C        PCK file to be sure about the body you're working with).
       C

       C
       C        Find the epoch for which the inertial to body equator and
       C        prime meridian transformation should be computed.
       C
                IF ( CORR .EQ. 'NONE' ) THEN
                   EPOCH = ET
                ELSE
                   EPOCH = ET - LT
                END IF

       C
       C        Find the transformation from inertial to body equator and
       C        prime meridian coordinates at the appropriate epoch.
       C
                CALL BODMAT ( TARGET, EPOCH, TIPM )

       C
       C        Negate the observer-target vector and rotate the result
       C        into body-fixed coordinates.
       C
                CALL VMINUS ( STATE, OBSPOS )

                CALL MXV    ( TIPM,  OBSPOS,  OBSPOS )

       C
       C        Look up the radii of the target body from the PCK file.
       C
                CALL BODVAR ( TARGET, 'RADII', N, RADII )

       C
       C        Find the nearest point on the body to the observer.
       C
                CALL NEARPT ( OBSPOS, RADII(1), RADII(2), RADII(3),
            .                 SUBPT,  ALT                           )

       C
       C        Find the latitude and longitude of the sub-observer point.
       C
                CALL RECLAT ( SUBPT,  DIST,  SUBLON,  SUBLAT )

       C
       C        Write out the results, converting radians to degrees.
       C
                WRITE(*,*) ' '
                WRITE(*,*) 'Sub-observer longitude (deg) ', DPR() * SUBLON
                WRITE(*,*) 'Sub-observer latitude  (deg) ', DPR() * SUBLAT
                WRITE(*,*) ' '

                WRITE (*,*) 'Do you wish to continue? (T/F)'
                READ   *,    CONT

             END DO

             END
Top ### Testing for occultation of a star We've mentioned that SURFPT can be used for a low-accuracy determination of whether a star is occulted by an extended body. Such a test would ignore factors such as topography of the body, differential stellar aberration, differential light-time corrections, and light-bending. Nonetheless, such a determination can be useful, especially for deciding which objects to draw in a graphical depiction of the sky as seen in the field of view of an instrument. We demonstrate this procedure, assuming that the observing body is the Earth, and the occulting body is Jupiter. The NAIF integer code for the Earth is 399, and that of Jupiter is 599 (see the SPK required reading, spk.req, for the complete set of NAIF integer body codes). The program below could be easily generalized to deal with occulting bodies other than Jupiter: the necessary change would be to make the ID code of the occulting body a variable. We're going to take the offset of the observer from the Earth's center into account because doing so gives us an opportunity to demonstrate the use of a number of SPICELIB routines. The actual parallax involved is not necessarily important.
             PROGRAM STAROC
       C
       C     Determine whether a star is occulted by Jupiter as seen from
       C     Earth at a specified surface location and specified UTC time.
       C     This program does not determine whether the Earth blocks the
       C     line of sight to the star.
       C
       C     Inputs to the program are:
       C
       C        -- Name of a leapseconds kernel file
       C
       C        -- Name of a PCK file containing radii and orientation
       C           models for Earth and Jupiter
       C
       C        -- Name of an SPK file containing ephemeris data for Earth
       C           and Jupiter at the time of observation
       C
       C        -- UTC time of observation
       C
       C        -- RA and Dec of star, relative to the J2000 inertial
       C           reference frame
       C
       C        -- Planetocentric latitude and longitude of the observer.
       C
       C     The program prints out the result to the user's terminal
       C     screen.
       C
       C     Test data:  the inputs
       C
       C        UTC:    1990 DEC 12 18:02.01
       C        RA:     136.10735
       C        DEC:     17.37258
       C        OBSLAT:  12.0
       C        OBSLON:   1.0
       C
       C     should cause this program to find that the star at the
       C     specified RA and Dec is occulted.
       C

       C
       C     SPICELIB functions
       C
             DOUBLE PRECISION      DPR
             DOUBLE PRECISION      RPD
       C
       C     Local parameters
       C
             INTEGER               FILEN
             PARAMETER           ( FILEN = 128 )

       C
       C     Local variables
       C
             CHARACTER*(30)        UTC
             CHARACTER*(FILEN)     LEAP
             CHARACTER*(FILEN)     PCK
             CHARACTER*(FILEN)     SPK

             DOUBLE PRECISION      DEC
             DOUBLE PRECISION      ET
             DOUBLE PRECISION      ETIPM  ( 3, 3 )
             DOUBLE PRECISION      LT
             DOUBLE PRECISION      OBSLAT
             DOUBLE PRECISION      OBSLON
             DOUBLE PRECISION      OBSPOS ( 3 )
             DOUBLE PRECISION      OBSPT  ( 3 )
             DOUBLE PRECISION      ORIGIN ( 3 )
             DOUBLE PRECISION      JTIPM  ( 3, 3 )
             DOUBLE PRECISION      RA
             DOUBLE PRECISION      RANGE
             DOUBLE PRECISION      REARTH ( 3 )
             DOUBLE PRECISION      RJUP   ( 3 )
             DOUBLE PRECISION      JSTATE ( 6 )
             DOUBLE PRECISION      UVEC   ( 3 )
             DOUBLE PRECISION      VSTAR  ( 3 )
             DOUBLE PRECISION      XPT    ( 3 )

             INTEGER               HANDLE
             INTEGER               N

             LOGICAL               FOUND

             DATA ORIGIN  /  3 * 0.D0                /

       C
       C     Convert the UTC epoch to ephemeris time.  Before doing this,
       C     we must load a leapseconds kernel into the kernel pool.
       C     Consult the TIME.REQ and KERNEL.REQ required reading for more
       C     information about NAIF time conversion routines.
       C
       C     (In your own program, you would substitute the name of an
       C     actual leapseconds kernel accessible to your program for the
       C     name shown below.)
       C
             WRITE (*,*) 'Enter name of leapseconds kernel'
             READ  (*,FMT='(A)') LEAP

             CALL FURNSH ( LEAP )

       C
       C     We're also going to need the radii of Earth and Jupiter, and
       C     the orientation of Earth and Jupiter at the desired epoch.
       C     In order to make the required information available to our
       C     code, we must load a P_constants kernel into the kernel pool
       C     (in your own program, you would have to supply the name of a
       C     P_constants kernel accessible to your own program.  See the
       C     KERNEL.REQ required reading for more information about
       C     P_constants kernels).
       C
             WRITE (*,*) 'Enter name of P_constants kernel'
             READ  (*,FMT='(A)') PCK

             CALL FURNSH ( PCK )

       C
       C     We'll need to load an SPK file containing ephemeris data for
       C     both Earth and Jupiter.
       C
             WRITE (*,*) 'Enter name of SPK file'
             READ  (*,FMT='(A)') SPK

             CALL SPKLEF ( SPK, HANDLE )

       C
       C     Now obtain the UTC time and convert UTC to ET.
       C
             WRITE (*,*) 'Enter time of observation, UTC'
             READ  (*,FMT='(A)') UTC

             CALL UTC2ET ( UTC, ET )

       C
       C     Obtain the RA and Dec of the star.
       C
             WRITE (*,*) 'Enter RA of star, relative to '               //
            .            'J2000 reference frame (degrees)'
             READ   *,    RA
             WRITE (*,*) 'Enter Dec of star, relative to '              //
            .            'J2000 reference frame (degrees)'
             READ   *,    DEC

       C
       C     Convert input RA and Dec to radians.
       C
             RA  = RA  * RPD()
             DEC = DEC * RPD()

       C
       C     Obtain the observation longitude and latitude.
       C
             WRITE (*,*) 'Enter planetocentric longitude of the '       //
            .            'observer (degrees)'
             READ   *,    OBSLON

             WRITE (*,*) 'Enter planetocentric latitude of the '        //
            .            'observer (degrees)'
             READ   *,    OBSLAT

             OBSLON = RPD() * OBSLON
             OBSLAT = RPD() * OBSLAT

       C
       C     Look up the radii used in the ellipsoidal models of the
       C     Earth and Jupiter.
       C
             CALL BODVAR ( 399, 'RADII', N, REARTH )
             CALL BODVAR ( 599, 'RADII', N, RJUP   )

       C
       C     Find the position of the observer relative to the Earth's
       C     center in cartesian coordinates.  Although it's a little
       C     tedious to do so, we'll find the observation point's exact
       C     altitude (according to our ellipsoidal Earth model).  To do
       C     this, we find the unit vector pointing from the Earth's
       C     center toward the observer's surface location, and compute
       C     the cartesian coordinates of the intercept point generated
       C     by extending this unit vector.
       C

       C
       C     Find the cartesian unit vector:
       C
             CALL LATREC ( 1.D0, OBSLON, OBSLAT, UVEC )

       C
       C     Extend the unit vector from the origin and find the surface
       C     intercept.  The output argument FOUND indicates whether the
       C     intercept point was found; there's no need to check this
       C     when the intercepting ray has its vertex inside the
       C     ellipsoid.  OBSPT is the surface point on output.
       C
             CALL SURFPT ( ORIGIN, UVEC, REARTH(1), REARTH(2), REARTH(3),
            .              OBSPT,  FOUND                                 )

       C
       C     Now find the vector from the Earth's center to the
       C     observation point in J2000 coordinates.  We'll need the
       C     transformation from J2000 to Earth body equator and prime
       C     meridian coordinates.  This transformation, which we'll name
       C     ETIPM, can be obtained from BODMAT.
       C
             CALL BODMAT ( 399, ET, ETIPM )

       C
       C     To convert OBSPT to J2000 coordinates, use the transpose
       C     of ETIPM.  The routine MTXV left-multiplies a vector by
       C     the transpose of a specified matrix.
       C
             CALL MTXV ( ETIPM, OBSPT, OBSPT )

       C
       C     Find the state of Jupiter, as seen from Earth, in J2000
       C     coordinates.  Use light time correction.  The output
       C     argument JSTATE is a state vector--the first three
       C     components are position; the next three are velocity.
       C
             CALL SPKEZ ( 599, ET, 'J2000', 'LT', 399, JSTATE, LT )

       C
       C     Find the J2000 to body equator and prime meridian
       C     transformation for Jupiter.  Note:  the epoch at which to
       C     evaluation this transformation must be adjusted for
       C     light time.  Anticipating our need, SPKEZ (called above)
       C     has returned this light time value in the argument LT.
       C
             CALL BODMAT ( 599, ET-LT, JTIPM )

       C
       C     Find the position of the Earth's center in Jupiter equator
       C     and prime meridian coordinates.  To do this, negate the
       C     Earth center-Jupiter center vector and rotate the result
       C     into Jupiter equator and prime meridian coordinates.  VMINUS
       C     negates a three-dimensional vector.
       C
             CALL VMINUS ( JSTATE, OBSPOS )
             CALL MXV    ( JTIPM,  OBSPOS, OBSPOS )

       C
       C     While we're at it, transform the observer's offset vector
       C     from J2000 to Jupiter equator and prime meridian coordinates.
       C
             CALL MXV  ( JTIPM, OBSPT, OBSPT )

       C
       C     Finally, we can find the position of the observer in
       C     Jupiter equator and prime meridian coordinates.
       C
             CALL VADD ( OBSPOS, OBSPT, OBSPOS )

       C
       C     Convert the RA and Dec of the star into a unit vector
       C     specified in J2000 coordinates.
       C
             CALL RADREC ( 1.D0, RA, DEC, VSTAR )

       C
       C     Rotate the star's direction vector into
       C     Jupiter equator and prime meridian coordinates.
       C
             CALL MXV  ( JTIPM, VSTAR, VSTAR )

       C
       C     Now that the observer's position and the star's direction
       C     vector are in Jupiter equator and prime meridian coordinates,
       C     we can apply SURFPT to find out whether the line segment
       C     from the observer to the star intersects Jupiter.  The point
       C     of working in this coordinate system is that our Jupiter
       C     axes are lined up with the coordinate axes of this system.
       C
       C     If SURFPT sets the output argument FOUND to .TRUE., the
       C     star is occulted.  Otherwise, it is visible.
       C
       C
             CALL SURFPT ( OBSPOS, VSTAR, RJUP(1), RJUP(2), RJUP(3),
            .              XPT,    FOUND                                )

             WRITE (*,*) ' '

             IF ( FOUND ) THEN
                WRITE (*,*) 'Star is occulted by Jupiter at time ', UTC
             ELSE
                WRITE (*,*) 'Star is not occulted by Jupiter at time ',UTC
             END IF

             END
Top ### Finding the `limb angle' of an instrument boresight If we want to find the angle of some ray above the limb of an ellipsoid, where the angle is measured in a plane containing the ray and a `down' vector, we can follow the procedure given below. We assume the ray does not intersect the ellipsoid. The result we seek is called ANGLE, imaginatively enough. We assume that all vectors are given in body-fixed coordinates.
       C
       C     The body-center--observer vector is OBSERV; RAYDIR is
       C     the boresight ray's direction vector in body-fixed
       C     coordinates.
       C
       C     Find the limb of the ellipsoid as seen from the
       C     point OBSERV.  Here A, B, and C are the lengths of
       C     the semi-axes of the ellipsoid.
       C
             CALL EDLIMB ( A, B, C, OBSERV, LIMB )

       C
       C     The ray direction vector is RAYDIR, so the ray is the
       C     set of points
       C
       C        OBSERV  +  t * RAYDIR
       C
       C     where t is any non-negative real number.
       C
       C     The `down' vector is just -OBSERV.  The vectors
       C     OBSERV and RAYDIR are spanning vectors for the plane
       C     we're interested in.  We can use PSV2PL to represent
       C     this plane by a SPICELIB plane.
       C
             CALL PSV2PL ( OBSERV, OBSERV, RAYDIR, PLANE )
       C
       C
       C     Find the intersection of the plane defined by OBSERV
       C     and RAYDIR with the limb.
       C
             CALL INELPL ( LIMB, PLANE, NXPTS, XPT1, XPT2 )
       C
       C
       C     We always expect two intersection points, if DOWN
       C     is valid.
       C
       C                  IF ( NXPTS .LT. 2 ) THEN
       C
       C                     [ do something about the error ]
       C
       C                  ENDIF
       C
       C
       C     Form the vectors from OBSERV to the intersection
       C     points.  Find the angular separation between the
       C     boresight ray and each vector from OBSERV to the
       C     intersection points.
       C
             CALL VSUB ( XPT1, OBSERV, VEC1 )
             CALL VSUB ( XPT2, OBSERV, VEC2 )

             SEP1 = VSEP ( VEC1, RAYDIR )
             SEP2 = VSEP ( VEC2, RAYDIR )

       C
       C     The angular separation we're after is the minimum of
       C     the two separations we've computed.
       C
             ANGLE = MIN ( SEP1, SEP2 )
Top ### Finding the instrument boresight---limb distance The technique demonstrated here can be used to find the altitude of a spacecraft-mounted remote sensing instrument's boresight ray above a the limb of a body modelled by a triaxial ellipsoid. We assume that the angular separation of the boresight ray and the `down' direction (the direction of the normal to the limb plane that points away from the observer) is less than pi/2 radians. In practice, this criterion will usually be met when the limb of the body is in the instrument field of view. For example, in the case of a camera with a field of view of eight milliradians (roughly the field of view of the Galileo SSI camera and slightly more than the field of view of the Voyager 2 narrow angle ISS camera), this technique will work for observations of the Earth as long as the limb of the Earth is in the field of view, and as long as the altitude of the observer above the Earth is at least 52 meters (assuming a spherical Earth with all radii equal to the equatorial radius). The following subroutine demonstrates the computation.
             SUBROUTINE RAYALT  ( TARG,   ET,  CORR,  SC,  RAYDIR,
            .                     LNEAR,  ALT                      )

       C
       C     Find the altitude of an instrument boresight ray above the
       C     limb of a specified target body, and find the nearest point
       C     on the limb to the ray.
       C
             INTEGER               TARG
             DOUBLE PRECISION      ET
             CHARACTER*(*)         CORR
             INTEGER               SC
             DOUBLE PRECISION      RAYDIR ( 3 )
             DOUBLE PRECISION      LNEAR  ( 3 )
             DOUBLE PRECISION      ALT

       C
       C     SPICELIB functions
       C
             DOUBLE PRECISION      DPR

       C
       C     Local parameters
       C
             INTEGER               UBEL
             PARAMETER           ( UBEL  =   9 )

             INTEGER               UBPL
             PARAMETER           ( UBPL  =   4 )

       C
       C     Local variables
       C
             DOUBLE PRECISION      CENTER ( 3 )
             DOUBLE PRECISION      DIST
             DOUBLE PRECISION      EPOCH
             DOUBLE PRECISION      LIMB   ( UBEL )
             DOUBLE PRECISION      NEAR   ( 3 )
             DOUBLE PRECISION      LPLANE ( UBPL )
             DOUBLE PRECISION      LT
             DOUBLE PRECISION      PJLIMB ( UBEL )
             DOUBLE PRECISION      PJPOS  ( 3 )
             DOUBLE PRECISION      PLANE  ( UBPL )
             DOUBLE PRECISION      SCPOS  ( 3 )
             DOUBLE PRECISION      SMAJOR ( 3 )
             DOUBLE PRECISION      SMINOR ( 3 )
             DOUBLE PRECISION      RADII  ( 3 )
             DOUBLE PRECISION      STATE  ( 6 )
             DOUBLE PRECISION      SUBPT  ( 3 )
             DOUBLE PRECISION      TIPM   ( 3, 3 )

             INTEGER               N

             LOGICAL               FOUND

       C
       C                  Glossary of variables
       C
       C
       C       Name                            Meaning
       C    ----------        ------------------------------------------
       C
       C     SC                NAIF ID code of a spacecraft.
       C
       C     TARG              NAIF ID code of a target body.
       C
       C     ET                Ephemeris time of the observation.
       C
       C     STATE             State (position and velocity), light-time
       C                       corrected, of the target body as seen from
       C                       the spacecraft at ET.
       C
       C     LT                Light time from target body to the
       C                       spacecraft at ET.
       C
       C     SCPOS             Spacecraft position relative to the
       C                       light-time corrected body position, in body
       C                       centered, inertial coordinates.
       C
       C     TIPM              Transformation from inertial (J2000)
       C                       coordinates to light-time corrected body
       C                       equator and prime meridian coordinates.
       C
       C     RADII             Radii of the triaxial ellipsoid used to
       C                       model the body.
       C
       C     RAYDIR            `Ray direction'--the instrument boresight
       C                       vector.
       C
       C     LIMB              A SPICELIB ellipse that represents the
       C                       body's limb. (This is an array of length
       C                       9.)
       C
       C     LPLANE            The limb plane.
       C
       C     PLANE             A plane orthogonal to the vector RAYDIR.
       C                       (This is an array of length 4.)
       C
       C     PJLIMB            The orthogonal projection of the limb onto
       C                       PLANE.
       C
       C     PJPOS             The orthogonal projection of the
       C                       spacecraft position onto PLANE.
       C
       C     NEAR              The nearest point on the projected limb
       C                       PJLIMB to the projected vertex PJPOS.
       C
       C     LNEAR             The inverse orthogonal projection of PJPOS
       C                       back to the limb plane. This is the point
       C                       on the limb closest to the boresight ray.
       C
       C     ALT               The altitude of the boresight above the
       C                       limb of the target body.

       C
       C
       C     Step 1:  Find the light-time corrected position of the body
       C              as seen from the spacecraft.  We will request the
       C              state vector in J2000 coordinates.
       C
       C              Also find the light-time corrected inertial to body
       C              equator and prime meridian transformation matrix.
       C
       C
             CALL SPKEZ  ( TARG, ET, 'J2000', CORR, SC, STATE, LT )

             CALL BODMAT ( TARG, ET-LT, TIPM )

       C
       C     Step 2:  Find the position of the spacecraft as seen from
       C              the light-time corrected body position.  Convert
       C              this position vector to equator and prime meridian
       C              coordinates.  Also transform the boresight ray.
       C
             CALL VMINUS ( STATE, SCPOS          )
             CALL MXV    ( TIPM,  SCPOS,  SCPOS  )
             CALL MXV    ( TIPM,  RAYDIR, RAYDIR )

       C
       C     Step 3:  Look up the radii of the body and find the limb.
       C              Find the limb plane as well, since we'll use it
       C              soon.
       C
             CALL BODVAR ( TARG, 'RADII', N, RADII )

             CALL EDLIMB ( RADII(1), RADII(2), RADII(3), SCPOS,  LIMB   )
             CALL EL2CGV ( LIMB,     CENTER,   SMAJOR,   SMINOR         )
             CALL PSV2PL (           CENTER,   SMAJOR,   SMINOR, LPLANE )

       C
       C     Step 4.  Create a plane that is orthogonal to the instrument
       C              boresight vector.  Orthogonally project the limb
       C              onto this plane.  Under this projection, every
       C              point on the limb maps to a point that is the same
       C              distance from the boresight as the pre-image of
       C              that point.  The intersection of
       C              the boresight itself with the plane is just the
       C              orthogonal projection of the boresight ray's vertex
       C              onto the plane.
       C
             CALL NVC2PL ( RAYDIR, 0.D0,   PLANE  )

             CALL PJELPL ( LIMB,   PLANE,  PJLIMB )

             CALL VPRJP  ( SCPOS,  PLANE,  PJPOS  )

       C
       C     Step 5.  Find the nearest point (NEAR) on the projected
       C              limb to the projected vertex of the boresight
       C              ray.  The distance between the projected ellipse
       C              PJLIMB and the point PJPOS is the altitude of the
       C              boresight ray above the limb.  To find the point
       C              on the limb closest to the boresight, just find
       C              the inverse orthogonal projection of NEAR onto the
       C              limb plane.
       C

             CALL NPELPT ( PJPOS,   PJLIMB,   NEAR,    ALT           )

             CALL VPRJPI ( NEAR,    PLANE,   LPLANE,  LNEAR,  FOUND )

       C
       C     If FOUND is false, the inverse projection could not be
       C     performed.  This indicates that the limb plane is too
       C     close to parallel to the boresight ray for this method to
       C     work.  This problem is very unlikely to occur.
       C
             IF ( .NOT. FOUND ) THEN
                CALL SETMSG ( 'Could not compute limb altitude.' )
                CALL SIGERR ( 'SPICE(DEGENERATECASE)'            )
                RETURN
             END IF

       C
       C     LNEAR and ALT are the limb point nearest to the boresight
       C     ray, and the altitude of the boresight ray above the limb of
       C     the target body, respectively.
       C
             END
Top # Mathematical notes Top ## Defining an ellipse parametrically Our aim is to show that the set of points
       CENTER    +    cos(theta) * V1    +    sin(theta) * V2
where CENTER, V1, and V2 are specified vectors in three-dimensional space, and where theta is a real number in the interval (-pi, pi], is in fact an ellipse as we've claimed. Since the vector CENTER simply translates the set, we may assume without loss of generality that it is the zero vector. So we'll re-write our expression for the alleged ellipse as
       cos(theta) * V1    +    sin(theta) * V2
where theta is a real number in the interval (-pi, pi]. We'll give the name S to the above set of vectors. Without loss of generality, we can assume that V1 and V2 lie in the x-y plane. Therefore, we can treat V1 and V2 as two-dimensional vectors. If V1 and V2 are linearly dependent, S is obviously a line segment or a point, so there is nothing to prove. We'll assume from now on that V1 and V2 are linearly independent. Every point in S has coordinates ( cos(theta), sin(theta) ) relative to the basis
       {V1, V2}.
Define the change-of-basis matrix C by setting the first and second columns of C equal to V1 and V2, respectively. If (x,y) are the coordinates of a point P of S relative to the standard basis
       { (1,0), (0,1) },
then the coordinates of P relative to the basis
       {V1, V2}
are
                  +- -+
             -1   | x |
            C     |   |
                  | y |
                  +- -+

                +-          -+
                | cos(theta) |
       =        |            |
                | sin(theta) |
                +-          -+
Taking inner products, we find
            +-    -+      -1 T     -1   +- -+
            | x  y |   ( C  )     C     | x |
            +-    -+                    |   |
                                        | y |
                                        +- -+

            +-                      -+  +-          -+
       =    | cos(theta)  sin(theta) |  | cos(theta) |
            +-                      -+  |            |
                                        | sin(theta) |
                                        +-          -+

       =    1
The matrix
          -1  T   -1
       ( C   )   C
is symmetric; let's say that it has entries
       +-          -+
       |   a   b/2  |
       |            |.
       |  b/2   c   |
       +-          -+
We know that a and c are positive because they are squares of norms of the columns of
        -1
       C
which is a non-singular matrix. Then the equation above reduces to
          2                2
       a x   +  b xy  + c y   =  1,     a, c  >  0.
We can find a new orthogonal basis such that this equation transforms to
           2           2
       d1 u    +   d2 v
with respect to this new basis. Let's give the name SYM to the matrix
       +-          -+
       |   a   b/2  |
       |            |;
       |  b/2   c   |
       +-          -+
since SYM is symmetric, there exists an orthogonal matrix M that diagonalizes SYM. That is, we can find an orthogonal matrix M such that
                        +-      -+
        T               | d1   0 |
       M  SYM  M    =   |        |.
                        | 0   d2 |
                        +-      -+
The existence of such a matrix M will not be proved here; see reference [2]. The columns of M are the elements of the basis we're looking for: if we define the variables (u,v) by the transformation
       +- -+        +- -+
       | u |      T | x |
       |   |  =  M  |   |,
       | v |        | y |
       +- -+        +- -+
then our equation in x and y transforms to the equation
           2           2
       d1 u    +   d2 v
since
            2                 2
           a x   +  b xy  +  c y

            +-    -+              +- -+
       =    | x  y |      SYM     | x |
            +-    -+              |   |
                                  | y |
                                  +- -+

            +-    -+   T          +- -+
       =    | u  v |  M   SYM  M  | u |
            +-    -+              |   |
                                  | v |
                                  +- -+

            +-    -+  +-      -+  +- -+
       =    | u  v |  | d1   0 |  | u |
            +-    -+  |        |  |   |
                      | 0   d2 |  | v |
                      +-      -+  +- -+

                2            2
       =    d1 u    +    d2 v
This last equation is that of an ellipse, as long as d1 and d2 are positive. To verify that they are, note that d1 and d2 are the eigenvalues of the matrix SYM, and SYM is the product
          -1  T   -1
       ( C   )   C,
which is of the form
        T
       M   M,
so SYM is positive semi-definite (its eigenvalues are non-negative). Furthermore, since the product
          -1  T   -1
       ( C   )   C
is non-singular if C is non-singular, and since the columns of C are V1 and V2, SYM exists and is non-singular precisely when V1 and V2 are linearly independent, a condition that we have assumed. So the eigenvalues of SYM can't be zero. They're not negative either. We conclude they're positive. Top ## Solving intersection problems There is one problem solving technique used in SPICELIB ellipse and ellipsoid routines that is so useful that it deserves special mention: using a `distortion map' to solve intersection problems. The distortion map (as it is referred to in SPICELIB routines) is simply a linear transformation that maps an ellipsoid to the unit sphere. The distortion map defined by an ellipsoid whose semi-axes are A, B, and C is represented by the matrix
       +-                -+
       |  1/A   0    0    |
       |   0   1/B   0    |.
       |   0    0    1/C  |
       +-                -+
The distortion map is (as is clear from examining the matrix) one-to-one and onto, and in particular is invertible, so it preserves set operations such as intersection. That is, if M is a distortion map and X, Y are two sets, then
       M( X intersect Y ) = M(X) intersect M(Y).
The same is true of the inverse of the distortion map. The utility of these facts is that frequently it's easier to find the intersection of the images under the distortion map of two sets than it is to find the intersection of the original two sets. Having found the intersection of the `distorted' sets, we apply the inverse distortion map to arrive at the intersection of the original sets. Some examples: `--` To find the intersection of a ray and an ellipsoid, apply the distortion map to both. Because the distortion map is linear, the ray maps to another ray, and the ellipsoid maps to the unit sphere. The resulting intersection problem is easy to solve. Having found the points of intersection of the new ray and the unit sphere, if any, we apply the inverse distortion map to these points, and we're done. `--` To find the intersection of a plane and an ellipsoid, apply the distortion map to both. The linearity of the distortion map ensures that the original plane maps to a second plane (whose formula is easily calculated). The ellipsoid maps to the unit sphere. The intersection of a plane and a unit sphere is easily found. The inverse distortion map is then applied to the circle of intersection (when the intersection is non-trivial), and the ellipse of intersection of the original plane and ellipsoid results. This procedure is used in the SPICELIB routine INEDPL. `--` To find the image under gnomonic projection onto a plane (camera projection) of an ellipsoid, given a focal point, we must find the intersection of the plane and the cone generated by ellipsoid and the focal point. Applying the distortion map to the ellipsoid, plane, and focal point, the problem is transformed into that of finding the intersection of the transformed plane with the cone generated by a unit sphere and the transformed focal point. This `transformed' problem is much easier to solve. The resulting intersection ellipse is then mapped back to the original intersection ellipse by the inverse distortion mapping. Top # Appendix: Document Revision History Top ### December 21, 2004 LDPOOL was replaced with FURNSH in all examples. Top ### December 12, 2002. Corrections were made to comments in the code example that computes the altitude of a ray above the limb of an ellipsoid. Previously, the quantity computed was incorrectly described as the altitude of a ray above an ellipsoid.