/* diags2.f -- translated by f2c (version 19980913).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Table of constant values */

static integer c__4 = 4;
static doublereal c_b6 = 1.;
static integer c__2 = 2;

/* $Procedure  DIAGS2   ( Diagonalize symmetric 2x2 matrix ) */
/* Subroutine */ int diags2_(doublereal *symmat, doublereal *diag, doublereal 
	*rotate)
{
    /* Initialized data */

    static doublereal ident[4]	/* was [2][2] */ = { 1.,0.,0.,1. };

    /* System generated locals */
    doublereal d__1, d__2, d__3;

    /* Local variables */
    doublereal a, b, c__, root1[2], root2[2], scale;
    extern /* Subroutine */ int chkin_(char *, ftnlen), vhatg_(doublereal *, 
	    integer *, doublereal *), moved_(doublereal *, integer *, 
	    doublereal *), rquad_(doublereal *, doublereal *, doublereal *, 
	    doublereal *, doublereal *), swapd_(doublereal *, doublereal *);
    doublereal eigvec[2];
    extern /* Subroutine */ int chkout_(char *, ftnlen);
    extern logical return_(void);

/* $ 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 */

/* $ Declarations */
/* $ 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(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. */

/* $ Files */

/*     None. */

/* $ Exceptions */

/*     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. */

/* $ 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.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]  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 */

/* -    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) */

/* -& */
/* $ Index_Entries */

/*     diagonalize symmetric 2x2_matrix */

/* -& */

/*     SPICELIB functions */

/*     Local variables */

/*     Saved variables */

/*     Initial values */

/*     Standard SPICE error handling. */

    if (return_()) {
	return 0;
    } else {
	chkin_("DIAGS2", (ftnlen)6);
    }

/*     We check for the case of a diagonal input matrix, since */
/*     eigenvector determination is simplified by ruling out this */
/*     case. */
    if (symmat[2] == 0.) {
	moved_(ident, &c__4, rotate);
	moved_(symmat, &c__4, diag);

/*        Explicity zero out the (2,1) entry of DIAG, since DIAG is */
/*        guaranteed to be diagonal. */

	diag[1] = 0.;
	chkout_("DIAGS2", (ftnlen)6);
	return 0;
    }

/*     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. */

/* Computing MAX */
    d__1 = abs(symmat[0]), d__2 = abs(symmat[2]), d__1 = max(d__1,d__2), d__2 
	    = abs(symmat[3]);
    scale = max(d__1,d__2);
    a = symmat[0] / scale;
    b = symmat[2] / scale;
    c__ = symmat[3] / 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. */

    d__1 = -(a + c__);
/* Computing 2nd power */
    d__3 = b;
    d__2 = a * c__ - d__3 * d__3;
    rquad_(&c_b6, &d__1, &d__2, root1, root2);

/*     ROOT1 is the root corresponding to the positive discriminant term; */
/*     this is guaranteed by RQUAD. */

    diag[0] = root1[0];
    diag[1] = 0.;
    diag[2] = 0.;
    diag[3] = 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 1 or 2. */

/*        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.  The */
/*        rigorous justification of this choice would literally take */
/*        pages of explanation, and we are not going to go through it */
/*        here.  In most cases, either choice is satisfactory, and in */
/*        the case where cancellation is a problem, our choice is */
/*        preferable. */

/*        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 = '2' instead */
/*        of x = '1'. */

    if ((d__1 = diag[0] - a, abs(d__1)) >= (d__2 = diag[0] - c__, abs(d__2))) 
	    {

/*        In this case, the second eigenvector component EIGVEC(2) */
/*        should be larger than |B|; we explain why in detail below. */
/*        We use the MAX function below to guard against reversal of the */
/*        inequality due to round-off error. */

	eigvec[0] = b;
/* Computing MAX */
	d__1 = diag[0] - a, d__2 = abs(b);
	eigvec[1] = max(d__1,d__2);

/*        Recall that DIAG(1,1) is an eigenvalue of the scaled version */
/*        of SYMMAT */

/*           +-      -+ */
/*           | A    B | */
/*           |        |. */
/*           | B    C | */
/*           +-      -+ */

/*        DIAG(1,1) is the positive-discriminant root of this matrix's */
/*        characteristic equation.  EIGVEC's y-component */

/*           DIAG(1,1) - 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(1,1) - A  =  -----   +     \  /  ----------    +    B */
/*                               2            \/       4 */

/*        and */
/*                                                 ______________________ */
/*                                                /         2 */
/*                             A - C             / ( A - C )           2 */
/*           DIAG(1,1) - 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(1,1) is the eigenvalue for the SECOND eigenvector, we */
/*        must swap the eigenvalues. */

/*        Unitize the eigenvector. */

	vhatg_(eigvec, &c__2, eigvec);
	rotate[0] = eigvec[1];
	rotate[1] = -eigvec[0];
	rotate[2] = eigvec[0];
	rotate[3] = eigvec[1];
	swapd_(diag, &diag[3]);
    } else {
/* Computing MAX */
	d__1 = diag[0] - c__, d__2 = abs(b);
	eigvec[0] = max(d__1,d__2);
	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_(eigvec, &c__2, eigvec);
	rotate[0] = eigvec[0];
	rotate[1] = eigvec[1];
	rotate[2] = -eigvec[1];
	rotate[3] = eigvec[0];
    }

/*     We must scale the eigenvalues. */

    diag[0] *= scale;
    diag[3] *= scale;
    chkout_("DIAGS2", (ftnlen)6);
    return 0;
} /* diags2_ */