// ---------------------------------------------------------------------------- // Numerical diagonalization of 3x3 matrcies // Copyright (C) 2006 Joachim Kopp // ---------------------------------------------------------------------------- // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either // version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU // Lesser General Public License for more details. // // You should have received a copy of the GNU Lesser General Public // License along with this library; if not, write to the Free Software // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA // ---------------------------------------------------------------------------- /* * Changes for IMOD: combined include, moved declaration up * \$Id\$ */ #include #include #include "dsyev3.h" // Constants #define M_SQRT3 1.73205080756887729352744634151 // sqrt(3) // Macros #define SQR(x) ((x)*(x)) // x^2 // ---------------------------------------------------------------------------- int dsyevc3(double A[3][3], double w[3]) // ---------------------------------------------------------------------------- // Calculates the eigenvalues of a symmetric 3x3 matrix A using Cardano's // analytical algorithm. // Only the diagonal and upper triangular parts of A are accessed. The access // is read-only. // ---------------------------------------------------------------------------- // Parameters: // A: The symmetric input matrix // w: Storage buffer for eigenvalues // ---------------------------------------------------------------------------- // Return value: // 0: Success // -1: Error // ---------------------------------------------------------------------------- { double m, c1, c0; // Determine coefficients of characteristic poynomial. We write // | a d f | // A = | d* b e | // | f* e* c | double de = A[0][1] * A[1][2]; // d * e double dd = SQR(A[0][1]); // d^2 double ee = SQR(A[1][2]); // e^2 double ff = SQR(A[0][2]); // f^2 double p, sqrt_p, q, c, s, phi; m = A[0][0] + A[1][1] + A[2][2]; c1 = (A[0][0]*A[1][1] + A[0][0]*A[2][2] + A[1][1]*A[2][2]) // a*b + a*c + b*c - d^2 - e^2 - f^2 - (dd + ee + ff); c0 = A[2][2]*dd + A[0][0]*ee + A[1][1]*ff - A[0][0]*A[1][1]*A[2][2] - 2.0 * A[0][2]*de; // c*d^2 + a*e^2 + b*f^2 - a*b*c - 2*f*d*e) p = SQR(m) - 3.0*c1; q = m*(p - (3.0/2.0)*c1) - (27.0/2.0)*c0; sqrt_p = sqrt(fabs(p)); phi = 27.0 * ( 0.25*SQR(c1)*(p - c1) + c0*(q + 27.0/4.0*c0)); phi = (1.0/3.0) * atan2(sqrt(fabs(phi)), q); c = sqrt_p*cos(phi); s = (1.0/M_SQRT3)*sqrt_p*sin(phi); w[1] = (1.0/3.0)*(m - c); w[2] = w[1] + s; w[0] = w[1] + c; w[1] -= s; return 0; }