| Ellipses and Ellipsoids Required Reading |
Table of ContentsEllipses and Ellipsoids Required Reading Abstract Introduction References Ellipse Data Type Description Ellipse and ellipsoid routines Constructing ellipses Access to ellipse data elements CGV2EL and EL2CGV are not inverses Triaxial ellipsoid routines Ellipse routines Examples Finding the `limb angle' of an instrument boresight Header examples Use of ellipses with planes Summary of routines Appendix A: Mathematical notes Defining an ellipse parametrically Solving intersection problems Appendix B: Document Revision History 2012 JAN 31, EDW (JPL) 2008 JAN 17, BVS (JPL) 2004 DEC 21, NAIF (JPL) 2002 DEC 12, NAIF (JPL) Ellipses and Ellipsoids Required Reading
Abstract
Introduction
Ellipses turn up frequently in the sort of science analysis problems SPICE is designed to help solve. The shapes of extended bodies--planets, satellites, and the Sun--are frequently modeled by triaxial ellipsoids. The IAU has defined such models for the Sun, all of the planets, and most of their satellites, in the 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 modeling orbits and planetary rings. References
Ellipse Data Type Description
ellipse = CENTER + cos(theta) * V1 + sin(theta) * V2where CENTER, V1, and V2 are 3-vectors, and theta is in the range
(-pi, pi].The set of points "ellipse" is an ellipse (see Appendix A: Mathematical notes). 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 coordinate origin (0, 0, 0,) is completely specified by its generating vectors. Further mention of the center or generating vectors for a particular ellipse, means 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. The internal design of the ellipse data type is not part of its specification. The design is an implementation choice based on the programming language and so the design may change. Users should not write code based on the current implementation; such code might fail when used with a future version of SPICE. NAIF implemented the SPICE ellipse data type in Fortran as arrays of double precision numbers whose elements are set and interpreted by a small set of access routines provided for that purpose. Do not access the ellipse elements in any other way. Arrays used as SPICE planes have length UBEL (`upper bound of ellipse'). Currently, UBEL is 9. As a matter of good programming practice, NAIF strongly recommends declaring ellipses with parameterized lengths as in this example:
INTEGER UBEL PARAMETER ( UBEL = 9 ) DOUBLE PRECISION ELLIPS ( UBEL )The elements of SPICE ellipses are set using CGV2EL (center and generating vectors to ellipse) and accessed using EL2CGV (ellipse to center and generating vectors). Ellipse and ellipsoid routinesConstructing ellipses
Let CENTER, V1, V2, and ELLIPS be declared by:
INTEGER UBEL PARAMETER ( UBEL = 9 ) DOUBLE PRECISION ELLIPS ( UBEL ) DOUBLE PRECISION CENTER ( 3 ) DOUBLE PRECISION V1 ( 3 ) DOUBLE PRECISION V2 ( 3 )After CENTER, V1, and V2 have been assigned values, you can construct a SPICE ellipse using CGV2EL:
CALL CGV2EL ( CENTER, V1, V2, ELLIPS )This call produces the SPICE 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 - Vto
CENTER + V Access to ellipse data elements
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
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
Ellipse routines
ExamplesFinding the `limb angle' of an instrument boresight
We assume that all vectors are given in body-fixed coordinates.
CALL EDLIMB ( A, B, C, OBSERV, LIMB )The ray direction vector is RAYDIR, so the ray is the set of points
OBSERV + t * RAYDIRwhere t is any non-negative real number. The `down' vector is just - OBSERV. The vectors OBSERV and RAYDIR are spanning vectors for the plane we're interested in. We can use PSV2PL to represent this plane by a SPICELIB plane.
CALL PSV2PL ( OBSERV, OBSERV, RAYDIR, PLANE )Find the intersection of the plane defined by OBSERV and RAYDIR with the limb.
CALL INELPL ( LIMB, PLANE, NXPTS, XPT1, XPT2 )We always expect two intersection points, if DOWN is valid. If NXPTS has value less-than two, the user must respond to the error condition. Form the vectors from OBSERV to the intersection points. Find the angular separation between the boresight ray and each vector from OBSERV to the intersection points.
CALL VSUB ( XPT1, OBSERV, VEC1 ) CALL VSUB ( XPT2, OBSERV, VEC2 ) SEP1 = VSEP ( VEC1, RAYDIR ) SEP2 = VSEP ( VEC2, RAYDIR )The angular separation we're after is the minimum of the two separations we've computed.
ANGLE = MIN ( SEP1, SEP2 ) Header examples
Use of ellipses with planes
Summary of routines
CGV2EL Center and generating vectors to ellipse EDLIMB Ellipsoid limb EDTERM Ellipsoid terminator EL2CGV Ellipse to center and generating vectors INEDPL Intersection of ellipsoid and plane INELPL Intersection of ellipse and plane NEARPT Nearest point on ellipsoid to point NPEDLN Nearest point on ellipsoid to line NPELPT Nearest point on ellipse to point PJELPL Projection of ellipse onto plane SAELGV Semi-axes of ellipse from generating vectors SINCPT Surface intercept SURFNM Surface normal on ellipsoid SURFPT Surface intercept point on ellipsoid SURFPV Surface point and velocity Appendix A: Mathematical notesDefining an ellipse parametrically
CENTER + cos(theta) * V1 + sin(theta) * V2where 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) * V2where 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 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 on 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
The distortion map (as it is referred to in SPICE 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:
Appendix B: Document Revision History2012 JAN 31, EDW (JPL)
Removed several obsolete examples. 2008 JAN 17, BVS (JPL)
2004 DEC 21, NAIF (JPL)
2002 DEC 12, NAIF (JPL)
|