diags2_c

Index Page diags2_c
A B C D E F G H I J K L M N O P Q R S T U V W X
Procedure Abstract Required_Reading Keywords Brief_I/O Detailed_Input | Detailed_Output Parameters Exceptions Files Particulars Examples | Restrictions Literature_References Author_and_Institution Version Index_Entries
#### Procedure
       void diags2_c ( ConstSpiceDouble    symmat [2][2],
                       SpiceDouble         diag   [2][2],
                       SpiceDouble         rotate [2][2]  )
#### Abstract
       Diagonalize a symmetric 2x2 matrix.
#### Required_Reading
       ROTATION
#### Keywords
       ELLIPSE
       MATRIX
       ROTATION
       TRANSFORMATION
#### 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.
#### 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.
#### Files
       None.
#### 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
#### Link to routine diags2_c source file diags2_c.c

    Thu Apr 16 08:57:44 2009