00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include <vector>
00026
00027 #include "WTensorFunctions.h"
00028
00029 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m )
00030 {
00031
00032
00033
00034 const double M_SQRT3 = 1.73205080756887729352744634151;
00035 double de = m( 1, 2 ) * m( 1, 0 );
00036 double dd = m( 1, 2 ) * m( 1, 2 );
00037 double ee = m( 1, 0 ) * m( 1, 0 );
00038 double ff = m( 2, 0 ) * m( 2, 0 );
00039 double m0 = m( 0, 0 ) + m( 1, 1 ) + m( 2, 2 );
00040 double c1 = m( 0, 0 ) * m( 1, 1 ) + m( 0, 0 ) * m( 2, 2 ) + m( 1, 1 ) * m( 2, 2 )
00041 - ( dd + ee + ff );
00042 double c0 = m( 2, 2 ) * dd + m( 0, 0 ) * ee + m( 1, 1 ) * ff - m( 0, 0 ) * m( 1, 1 ) * m( 2, 2 ) - 2. * m( 2, 0 ) * de;
00043
00044 double p, sqrt_p, q, c, s, phi;
00045 p = m0 * m0 - 3. * c1;
00046 q = m0 * ( p - ( 3. / 2. ) * c1 ) - ( 27. / 2. ) * c0;
00047 sqrt_p = sqrt( fabs( p ) );
00048
00049 phi = 27. * ( 0.25 * c1 * c1 * ( p - c1 ) + c0 * ( q + 27. / 4. * c0 ) );
00050 phi = ( 1. / 3. ) * atan2( sqrt( fabs( phi ) ), q );
00051
00052 c = sqrt_p * cos( phi );
00053 s = ( 1. / M_SQRT3 ) * sqrt_p * sin( phi );
00054
00055 std::vector< double > w( 3 );
00056
00057 w[ 2 ] = ( 1. / 3. ) * ( m0 - c );
00058 w[ 1 ] = w[ 2 ] + s;
00059 w[ 0 ] = w[ 2 ] + c;
00060 w[ 2 ] -= s;
00061 return w;
00062 }