Procedure

     DIAGS2   ( Diagonalize symmetric 2x2 matrix )

     SUBROUTINE DIAGS2 ( SYMMAT, DIAG, ROTATE )

Abstract

     Diagonalize a symmetric 2x2 matrix.

Required_Reading

     ROTATION

Keywords

     ELLIPSE
     MATRIX
     ROTATION
     TRANSFORMATION

Declarations

     IMPLICIT NONE

     DOUBLE PRECISION      SYMMAT ( 2, 2 )
     DOUBLE PRECISION      DIAG   ( 2, 2 )
     DOUBLE PRECISION      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   is 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(1,1)
                 SYMMAT(1,2)
                 SYMMAT(2,2)

              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
              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(2,1) is not used in this routine's
         computations, so the condition

            SYMMAT(1,2)  .NE.  SYMMAT(2,1)

         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 modeled 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.0D0    4.0D0  |
            |                  |
            |  4.0D0   -5.0D0  |
            +-                -+

         Then SYMMAT is similar to the diagonal matrix

            +-                -+
            |  3.0D0    0.0D0  |
            |                  |
            |  0.0D0   -7.0D0  |
            +-                -+

         so

            DIAG(1,1) =  3.D0
            DIAG(2,1) =  0.D0
            DIAG(1,2) =  0.D0
            DIAG(2,2) = -7.D0

         and ROTATE is

            +-                                   -+
            |   0.894427191          -0.447213595 |
            |                                     |
            |   0.447213595           0.894427191 |
            +-                                   -+

        which is an approximation to

            +-                                   -+
            |  0.4 * 5**(1/2)     -0.2 * 5**(1/2) |
            |                                     |
            |  0.2 * 5**(1/2)      0.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(1,1)  =  27.D0
            SYMMAT(2,1)  =   5.D0
            SYMMAT(1,2)  =   5.D0
            SYMMAT(2,2)  =   3.D0

            CALL DIAGS2 ( 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(1,1)   =  28.D0
            DIAG(2,1)   =  0.D0
            DIAG(1,2)   =  0.D0
            DIAG(2,2)   =  2.D0

          and

            ROTATE(1,1) =  0.980580676D0
            ROTATE(2,1) =  0.196116135D0
            ROTATE(1,2) = -0.196116135D0
            ROTATE(2,2) =  0.980580676D0

         The columns of ROTATE give the ellipse's axes, after scaling
         them by

                   1                            1
            ----------------     and     ---------------
              ____________                 ____________
            \/  DIAG(1,1)                \/  DIAG(2,2)

         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. Negative or zero
         eigenvalues will occur only as a result of round-off error;
         mathematically, the eigenvalues of the matrix SYMMAT are
         guaranteed to be positive, since they are the reciprocals of
         the squares of the lengths of the ellipse's semi-axes.

            DO I = 1, 2
               SMAJOR(I) = ROTATE(I,1)  /  DSQRT( DIAG(1,1) )
               SMINOR(I) = ROTATE(I,2)  /  DSQRT( DIAG(2,2) )
            END DO

Restrictions

     None.

Literature_References

     [1]  T. Apostol, "Calculus, Vol. II," chapter 5, "Eigenvalues of
          Operators Acting on Euclidean Spaces," John Wiley & Sons,
          1969.

Author_and_Institution

     N.J. Bachman       (JPL)
     J. Diaz del Rio    (ODC Space)
     W.L. Taber         (JPL)
     E.D. Wright        (JPL)

Version

    SPICELIB Version 1.3.0, 17-JUN-2021 (JDR)

        Added IMPLICIT NONE statement.

        Edited the header to comply with NAIF standard. Added
        "Error free." to $Exceptions section. Removed unnecessary
        $Revisions section.

    SPICELIB Version 1.2.0, 06-SEP-2005 (NJB)

        Updated to remove non-standard use of duplicate arguments
        in VHATG and SWAPD calls.

    SPICELIB Version 1.1.0, 24-JAN-2002 (EDW)

        Edited incorrect examples in the header. The example
        outputs did not correspond to the actual function
        of the routine.

    SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)

        Comment section for permuted index source lines was added
        following the header.

    SPICELIB Version 1.0.0, 04-NOV-1990 (NJB)