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. Ellipses =========================================================================== Revisions -------------------------------------------------------- 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. 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. 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. 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 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. 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 ) 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 ) 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 `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. 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. 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 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. 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. 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. 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 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. 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. 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 ) 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 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. 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. 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. Examples =========================================================================== 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, TIME, and KERNEL 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 LDPOOL ( LEAP ) WRITE (*,*) 'Enter name of P_constants kernel' READ (*,FMT='(A)') PCK CALL LDPOOL ( 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 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 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 and KERNEL 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 LDPOOL ( 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 required reading for more information about C P_constants kernels). C WRITE (*,*) 'Enter name of P_constants kernel' READ (*,FMT='(A)') PCK CALL LDPOOL ( 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 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 ) 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 Mathematical notes =========================================================================== 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. 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.