/* -Procedure diags2_c ( Diagonalize symmetric 2x2 matrix ) -Abstract Diagonalize a symmetric 2x2 matrix. -Copyright Copyright (1995), California Institute of Technology. U.S. Government sponsorship acknowledged. -Required_Reading ROTATIONS -Keywords ELLIPSE MATRIX ROTATION TRANSFORMATION */ #include "SpiceUsr.h" #include "SpiceZfc.h" #include "SpiceZmc.h" #undef diags2_c void diags2_c ( ConstSpiceDouble symmat [2][2], SpiceDouble diag [2][2], SpiceDouble rotate [2][2] ) /* -Brief_I/O Variable I/O Description -------- --- -------------------------------------------------- symmat I A symmetric 2x2 matrix. diag O A diagonal matrix similar to symmat. rotate O A rotation used as the similarity transformation. -Detailed_Input symmat A symmetric 2x2 matrix. That is, symmat has the form +- -+ | A B | | |. | B C | +- -+ This routine uses only the upper-triangular elements of symmat, that is, the elements symmat[0][0] symmat[0][1] symmat[1][1] to determine the outputs diag and rotate. -Detailed_Output diag, rotate are, respectively, a diagonal matrix and a 2x2 rotation matrix that satisfy the equation T diag = rotate * symmat * rotate. In other words, diag is similar to symmat, and rotate is a change-of-basis matrix that diagonalizes symmat. diags2_c chooses rotate so that its angle of rotation has the smallest possible magnitude. If there are two rotations that meet these criteria (they will be inverses of one another), either rotation may be chosen. -Parameters None. -Files None. -Exceptions Error free. 1) The matrix element symmat[1][0] is not used in this routine's computations, so the condition symmat[0][1] != symmat[1][0] has no effect on this routine's outputs. -Particulars The capability of diagonalizing a 2x2 symmetric matrix is especially useful in a number of geometric applications involving quadratic curves such as ellipses. Such curves are described by expressions of the form 2 2 A x + B xy + C y + D x + E y + F = 0. Diagonalization of the matrix +- -+ | A B/2 | | | | B/2 C | +- -+ allows us to perform a coordinate transformation (a rotation, specifically) such that the equation of the curve becomes 2 2 P u + Q v + R u + S v + T = 0 in the transformed coordinates. This form is much easier to handle. If the quadratic curve in question is an ellipse, we can easily find its center, semi-major axis, and semi-minor axis from the second equation. Ellipses turn up frequently in navigation geometry problems; for example, the limb and terminator (if we treat the Sun as a point source) of a body modelled as a tri-axial ellipsoid are ellipses. A mathematical note: because symmat is symmetric, we can ALWAYS find an orthogonal similarity transformation that diagonalizes symmat, and we can choose the similarity transformation to be a rotation matrix. By `orthogonal' we mean that if the rotate is the matrix in question, then T T rotate rotate = rotate rotate = I. The reasons this routine handles only the 2x2 case are: first, the 2x2 case is much simpler than the general case, in which iterative diagonalization methods must be used, and second, the 2x2 case is adequate for solving problems involving ellipses in 3 dimensional space. Finally, this routine can be used to support a routine that solves the general-dimension diagonalization problem for symmetric matrices. Another feature of the routine that might provoke curiosity is its insistence on choosing the diagonalization matrix that rotates the original basis vectors by the smallest amount. The rotation angle of rotate is of no concern for most applications, but can be important if this routine is used as part of an iterative diagonalization method for higher-dimensional matrices. In that case, it is most undesirable to interchange diagonal matrix elements willy-nilly; the matrix to be diagonalized could get ever closer to being diagonal without converging. Choosing the smallest rotation angle precludes this possibility. -Examples 1) A case that can be verified by hand computation: Suppose symmat is +- -+ | 1.0 4.0 | | | | 4.0 -5.0 | +- -+ Then symmat is similar to the diagonal matrix +- -+ | 3.0 0.0 | | | | 0.0 -7.0 | +- -+ so diag[0][0] = 3. diag[1][0] = 0. diag[0][1] = 0. diag[1][1] = -7. and rotate is +- -+ | 0.89442719099991588 -0.44721359549995794 | | | | 0.44721359549995794 0.89442719099991588 | +- -+ which is an approximation to +- -+ | .4 * 5**(1/2) -.2 * 5**(1/2) | | | | .2 * 5**(1/2) .4 * 5**(1/2) | +- -+ 2) Suppose we want to find the semi-axes of the ellipse defined by 2 2 27 x + 10 xy + 3 y = 1 We can write the above equation as the matrix equation +- -+ +- -+ +- -+ | x y | | 27 5 | | x | = 1 +- -+ | | | | | 5 3 | | y | +- -+ +- -+ Let symmat be the symmetric matrix on the left. The code fragment symmat[0][0] = 27.0; symmat[1][0] = 5.0; symmat[0][1] = 5.0; symmat[1][1] = 3.0; diags2_c ( symmat, diag, rotate ); will return diag, an array containing the eigenvalues of symmat, and rotate, the coordinate transformation required to diagonalize symmat. In this case, diag[0][0] = 28. diag[1][0] = 0. diag[0][1] = 0. diag[1][1] = 2. and rotate[0][0] = 0.980580675690920 rotate[1][0] = 0.196116135138184 rotate[0][1] = -0.196116135138184 rotate[1][1] = 0.980580675690920 The columns of rotate give the ellipse's axes, after scaling them by 1 1 ---------------- and --------------- ____________ ____________ \/ diag[0][0] \/ diag[1][1] respectively. If smajor and sminor are semi-major and semi-minor axes, we can find them as shown below. For brevity, we omit the check for zero or negative eigenvalues. for ( i = 0; i < 2; i++ ) { smajor[i] = rotate[i][0] / sqrt( diag[0][0] ); sminor[i] = rotate[i][1] / sqrt( diag[1][1] ); } -Restrictions None. -Literature_References [1] Calculus, Vol. II. Tom Apostol. John Wiley & Sons, 1969. See Chapter 5, `Eigenvalues of Operators Acting on Euclidean Spaces'. -Author_and_Institution N.J. Bachman (JPL) -Version -CSPICE Version 1.0.0, 13-JUL-1999 (NJB) -Index_Entries diagonalize symmetric 2x2_matrix -& */ { /* Begin diags2_c */ /* Local constants */ /* Static variables */ static SpiceDouble ident [2][2] = { {1., 0.}, {0., 1.} }; /* Local variables */ SpiceDouble a; SpiceDouble b; SpiceDouble c; SpiceDouble eigvec [2]; SpiceDouble root1 [2]; SpiceDouble root2 [2]; SpiceDouble scale; /* Error free. */ /* We check for the case of a diagonal input matrix, since eigenvector determination is simplified by ruling out this case. */ if ( symmat [0][1] == 0. ) { MOVED ( ident, 4, rotate ); MOVED ( symmat, 4, diag ); /* Explicitly zero out the [1][0] entry of diag, since diag is guaranteed to be diagonal. */ diag[1][0] = 0.0; return; } /* Getting here means there's some actual work to do. We start out by scaling our matrix, in order to reduce the chance of overflow. We divide everything by the largest magnitude of any element of symmat. We're guaranteed that scale is non-zero, since the 0 matrix is diagonal. */ scale = MaxAbs ( symmat[0][0], symmat[0][1] ); scale = MaxAbs ( scale, symmat[1][1] ); a = symmat[0][0] / scale; b = symmat[0][1] / scale; c = symmat[1][1] / scale; /* Compute the eigenvalues of the scaled version of symmat. The eigenvalues are roots of the equation det ( (1 / scale) * symmat - x * identity ) = 0, or equivalently, 2 2 x - ( a + c ) x + ( ac - b ) = 0. */ rquad_c ( 1.0, -(a + c), a*c - b*b, root1, root2 ); /* root1 is the root corresponding to the positive discriminant term; this is guaranteed by rquad_c. */ diag[0][0] = root1[0]; diag[1][0] = 0.; diag[0][1] = 0.; diag[1][1] = root2[0]; /* Our next job is to find an eigenvector corresponding to the eigenvalue of smaller magnitude. We can unitize it and choose an orthogonal unit vector so as to create the desired rotation matrix. If our original matrix is +- -+ | a b | | |, | b c | +- -+ then the matrix +- -+ | a - diag[x][x] b | | | | b c - diag[x][x] | +- -+ maps to zero all elements of the eigenspace corresponding to diag[x][x], where x is either 0 or 1. So +- -+ +- -+ | b | | diag[x][x] - c | | | and | | | diag[x][x] - a | | b | +- -+ +- -+ are candidates for eigenvectors for diag[x][x]. To minimize loss of accuracy in our eigenvector due to subtraction of nearly equal quantities, we choose the vector in which the term involving the eigenvalue has the larger magnitude. Note that there is nothing to be gained as far as accuracy is concerned by working with one eigenvalue as opposed to the other: the magnitudes of the quantities diag[x][x] - a and diag[x][x] - c would be interchanged by taking x = 1 instead of x = 0. */ if ( fabs( diag[0][0] - a ) >= fabs( diag[0][0] - c ) ) { /* In this case, the second eigenvector component eigvec[1] should be larger than |b|; we explain why in detail below. We use the MaxVal macro below to guard against reversal of the inequality due to round-off error. */ eigvec[0] = b; eigvec[1] = MaxVal ( diag[0][0] - a, fabs(b) ); /* Recall that diag[0][0] is an eigenvalue of the scaled version of symmat +- -+ | a b | | |. | b c | +- -+ diag[0][0] is the positive-discriminant root of this matrix's characteristic equation. eigvec's y-component diag[0][0] - a is positive and of magnitude at least as large as that of B, since it is the larger of ______________________ / 2 c - a / ( a - c ) 2 diag[0][0] - a = ----- + \ / ---------- + b 2 \/ 4 and ______________________ / 2 a - c / ( a - c ) 2 diag[0][0] - c = ----- + \ / ---------- + b 2 \/ 4 Equality between these expressions can occur only when a is equal to c, in which case both expressions are equal (except for round-off error) to |b|. So the argument of eigvec is in the interval [pi/4, 3*pi/4]. The second eigenvector is eigvec, and the first eigenvector is found by rotating eigvec by -pi/2. Since diag[0][0] is the eigenvalue for the SECOND eigenvector, we must swap the eigenvalues. */ /* Unitize the eigenvector. */ vhatg_c ( eigvec, 2, eigvec ); rotate[0][0] = eigvec[1]; rotate[1][0] = -eigvec[0]; rotate[0][1] = eigvec[0]; rotate[1][1] = eigvec[1]; swapd_ ( &(diag[0][0]), &(diag[1][1]) ); } else { eigvec[0] = MaxVal ( diag[0][0] - c, fabs(b) ); eigvec[1] = b; /* The x-component of eigvec is positive and has magnitude greater than or equal to that of the y-component of eigvec. The argument of eigvec is in [-pi/4, pi/4], and the second eigenvector is found by rotating eigvec by pi/2. */ /* Unitize the eigenvector. */ vhatg_c ( eigvec, 2, eigvec ); rotate[0][0] = eigvec[0]; rotate[1][0] = eigvec[1]; rotate[0][1] = -eigvec[1]; rotate[1][1] = eigvec[0]; } /* We must scale the eigenvalues. */ diag[0][0] *= scale; diag[1][1] *= scale; } /* End diags2_c */