OpenWalnut
1.4.0
|
00001 //--------------------------------------------------------------------------- 00002 // 00003 // Project: OpenWalnut ( http://www.openwalnut.org ) 00004 // 00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS 00006 // For more information see http://www.openwalnut.org/copying 00007 // 00008 // This file is part of OpenWalnut. 00009 // 00010 // OpenWalnut is free software: you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as published by 00012 // the Free Software Foundation, either version 3 of the License, or 00013 // (at your option) any later version. 00014 // 00015 // OpenWalnut is distributed in the hope that it will be useful, 00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00018 // GNU Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public License 00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>. 00022 // 00023 //--------------------------------------------------------------------------- 00024 00025 #include <vector> 00026 00027 #include "WTensorFunctions.h" 00028 00029 std::ostream& operator<<( std::ostream& os, const RealEigenSystem& sys ) 00030 { 00031 os << sys[0].first << ", " << sys[0].second << std::endl; 00032 os << sys[1].first << ", " << sys[1].second << std::endl; 00033 os << sys[2].first << ", " << sys[2].second << std::endl; 00034 return os; 00035 } 00036 00037 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m ) 00038 { 00039 // this is copied from the gpu glyph shader 00040 // src/graphicsEngine/shaders/tensorTools.fs 00041 // originally implemented by Mario Hlawitschka 00042 const double M_SQRT3 = 1.73205080756887729352744634151; 00043 double de = m( 1, 2 ) * m( 1, 0 ); 00044 double dd = m( 1, 2 ) * m( 1, 2 ); 00045 double ee = m( 1, 0 ) * m( 1, 0 ); 00046 double ff = m( 2, 0 ) * m( 2, 0 ); 00047 double m0 = m( 0, 0 ) + m( 1, 1 ) + m( 2, 2 ); 00048 double c1 = m( 0, 0 ) * m( 1, 1 ) + m( 0, 0 ) * m( 2, 2 ) + m( 1, 1 ) * m( 2, 2 ) 00049 - ( dd + ee + ff ); 00050 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; 00051 00052 double p, sqrt_p, q, c, s, phi; 00053 p = m0 * m0 - 3. * c1; 00054 q = m0 * ( p - ( 3. / 2. ) * c1 ) - ( 27. / 2. ) * c0; 00055 sqrt_p = sqrt( fabs( p ) ); 00056 00057 phi = 27. * ( 0.25 * c1 * c1 * ( p - c1 ) + c0 * ( q + 27. / 4. * c0 ) ); 00058 phi = ( 1. / 3. ) * atan2( sqrt( fabs( phi ) ), q ); 00059 00060 c = sqrt_p * cos( phi ); 00061 s = ( 1. / M_SQRT3 ) * sqrt_p * sin( phi ); 00062 00063 std::vector< double > w( 3 ); 00064 // fix: swapped w[ 2 ] and w[ 1 ] 00065 w[ 2 ] = ( 1. / 3. ) * ( m0 - c ); 00066 w[ 1 ] = w[ 2 ] + s; 00067 w[ 0 ] = w[ 2 ] + c; 00068 w[ 2 ] -= s; 00069 return w; 00070 }