src/math/SquareMatrix3.hpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program 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.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program 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 program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */

/**
 * @file SquareMatrix3.hpp
 * @author Teng Lin
 * @date 10/11/2004
 * @version 1.0
 */
#ifndef MATH_SQUAREMATRIX3_HPP
#define  MATH_SQUAREMATRIX3_HPP

#include "Quaternion.hpp"
#include "SquareMatrix.hpp"
#include "Vector3.hpp"

namespace oopse {

    template<typename Real>
    class SquareMatrix3 : public SquareMatrix<Real, 3> {
        public:
            
            /** default constructor */
            SquareMatrix3() : SquareMatrix<Real, 3>() {
            }

            /** copy  constructor */
            SquareMatrix3(const SquareMatrix<Real, 3>& m)  : SquareMatrix<Real, 3>(m) {
            }

            SquareMatrix3( const Vector3<Real>& eulerAngles) {
                setupRotMat(eulerAngles);
            }
            
            SquareMatrix3(Real phi, Real theta, Real psi) {
                setupRotMat(phi, theta, psi);
            }

            SquareMatrix3(const Quaternion<Real>& q) {
                *this = q.toRotationMatrix3();
            }

            SquareMatrix3(Real w, Real x, Real y, Real z) {
                Quaternion<Real> q(w, x, y, z);
                *this = q.toRotationMatrix3();
            }
            
            /** copy assignment operator */
            SquareMatrix3<Real>& operator =(const SquareMatrix<Real, 3>& m) {
                if (this == &m)
                    return *this;
                 SquareMatrix<Real, 3>::operator=(m);
                 return *this;
            }

            /**
             * Sets this matrix to a rotation matrix by three euler angles
             * @ param euler
             */
            void setupRotMat(const Vector3<Real>& eulerAngles) {
                setupRotMat(eulerAngles[0], eulerAngles[1], eulerAngles[2]);
            }

            /**
             * Sets this matrix to a rotation matrix by three euler angles
             * @param phi
             * @param theta
             * @psi theta
             */
            void setupRotMat(Real phi, Real theta, Real psi) {
                Real sphi, stheta, spsi;
                Real cphi, ctheta, cpsi;

                sphi = sin(phi);
                stheta = sin(theta);
                spsi = sin(psi);
                cphi = cos(phi);
                ctheta = cos(theta);
                cpsi = cos(psi);

                data_[0][0] = cpsi * cphi - ctheta * sphi * spsi;
                data_[0][1] = cpsi * sphi + ctheta * cphi * spsi;
                data_[0][2] = spsi * stheta;
                
                data_[1][0] = -spsi * ctheta - ctheta * sphi * cpsi;
                data_[1][1] = -spsi * stheta + ctheta * cphi * cpsi;
                data_[1][2] = cpsi * stheta;

                data_[2][0] = stheta * sphi;
                data_[2][1] = -stheta * cphi;
                data_[2][2] = ctheta;
            }


            /**
             * Sets this matrix to a rotation matrix by quaternion
             * @param quat
            */
            void setupRotMat(const Quaternion<Real>& quat) {
                *this = quat.toRotationMatrix3();
            }

            /**
             * Sets this matrix to a rotation matrix by quaternion
             * @param w the first element 
             * @param x the second element
             * @param y the third element
             * @param z the fourth element
            */
            void setupRotMat(Real w, Real x, Real y, Real z) {
                Quaternion<Real> q(w, x, y, z);
                *this = q.toRotationMatrix3();
            }

            /**
             * Returns the quaternion from this rotation matrix
             * @return the quaternion from this rotation matrix
             * @exception invalid rotation matrix
            */            
            Quaternion<Real> toQuaternion() {
                Quaternion<Real> q;
                Real t, s;
                Real ad1, ad2, ad3;    
                t = data_[0][0] + data_[1][1] + data_[2][2] + 1.0;

                if( t > 0.0 ){

                    s = 0.5 / sqrt( t );
                    q[0] = 0.25 / s;
                    q[1] = (data_[1][2] - data_[2][1]) * s;
                    q[2] = (data_[2][0] - data_[0][2]) * s;
                    q[3] = (data_[0][1] - data_[1][0]) * s;
                } else {

                    ad1 = fabs( data_[0][0] );
                    ad2 = fabs( data_[1][1] );
                    ad3 = fabs( data_[2][2] );

                    if( ad1 >= ad2 && ad1 >= ad3 ){

                        s = 2.0 * sqrt( 1.0 + data_[0][0] - data_[1][1] - data_[2][2] );
                        q[0] = (data_[1][2] + data_[2][1]) / s;
                        q[1] = 0.5 / s;
                        q[2] = (data_[0][1] + data_[1][0]) / s;
                        q[3] = (data_[0][2] + data_[2][0]) / s;
                    } else if ( ad2 >= ad1 && ad2 >= ad3 ) {
                        s = sqrt( 1.0 + data_[1][1] - data_[0][0] - data_[2][2] ) * 2.0;
                        q[0] = (data_[0][2] + data_[2][0]) / s;
                        q[1] = (data_[0][1] + data_[1][0]) / s;
                        q[2] = 0.5 / s;
                        q[3] = (data_[1][2] + data_[2][1]) / s;
                    } else {

                        s = sqrt( 1.0 + data_[2][2] - data_[0][0] - data_[1][1] ) * 2.0;
                        q[0] = (data_[0][1] + data_[1][0]) / s;
                        q[1] = (data_[0][2] + data_[2][0]) / s;
                        q[2] = (data_[1][2] + data_[2][1]) / s;
                        q[3] = 0.5 / s;
                    }
                }             

                return q;
                
            }

            /**
             * Returns the euler angles from this rotation matrix
             * @return the euler angles in a vector 
             * @exception invalid rotation matrix
             * We use so-called "x-convention", which is the most common definition. 
             * In this convention, the rotation given by Euler angles (phi, theta, psi), where the first 
             * rotation is by an angle phi about the z-axis, the second is by an angle  
             * theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
             * z-axis (again). 
            */            
            Vector3<Real> toEulerAngles() {
                Vector<Real> myEuler;
                Real phi,theta,psi,eps;
                Real ctheta,stheta;
                
                // set the tolerance for Euler angles and rotation elements

                theta = acos(min(1.0,max(-1.0,data_[2][2])));
                ctheta = data_[2][2]; 
                stheta = sqrt(1.0 - ctheta * ctheta);

                // when sin(theta) is close to 0, we need to consider singularity
                // In this case, we can assign an arbitary value to phi (or psi), and then determine 
                // the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
                // in cases of singularity.  
                // we use atan2 instead of atan, since atan2 will give us -Pi to Pi. 
                // Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
                // change the sign of both of the parameters passed to atan2.

                if (fabs(stheta) <= oopse::epsilon){
                    psi = 0.0;
                    phi = atan2(-data_[1][0], data_[0][0]);  
                }
                // we only have one unique solution
                else{    
                    phi = atan2(data_[2][0], -data_[2][1]);
                    psi = atan2(data_[0][2], data_[1][2]);
                }

                //wrap phi and psi, make sure they are in the range from 0 to 2*Pi
                if (phi < 0)
                  phi += M_PI;

                if (psi < 0)
                  psi += M_PI;

                myEuler[0] = phi;
                myEuler[1] = theta;
                myEuler[2] = psi;

                return myEuler;
            }
            
            /** Returns the determinant of this matrix. */
            Real determinant() const {
                Real x,y,z;

                x = data_[0][0] * (data_[1][1] * data_[2][2] - data_[1][2] * data_[2][1]);
                y = data_[0][1] * (data_[1][2] * data_[2][0] - data_[1][0] * data_[2][2]);
                z = data_[0][2] * (data_[1][0] * data_[2][1] - data_[1][1] * data_[2][0]);

                return(x + y + z);
            }            
            
            /**
             * Sets the value of this matrix to  the inversion of itself. 
             * @note since simple algorithm can be applied to inverse the 3 by 3 matrix, we hide the 
             * implementation of inverse in SquareMatrix class
             */
            SquareMatrix3<Real>  inverse() {
                SquareMatrix3<Real> m;
                double det = determinant();
                if (fabs(det) <= oopse::epsilon) {
                //"The method was called on a matrix with |determinant| <= 1e-6.",
                //"This is a runtime or a programming error in your application.");
                }

                m(0, 0) = data_[1][1]*data_[2][2] - data_[1][2]*data_[2][1];
                m(1, 0) = data_[1][2]*data_[2][0] - data_[1][0]*data_[2][2];
                m(2, 0) = data_[1][0]*data_[2][1] - data_[1][1]*data_[2][0];
                m(0, 1) = data_[2][1]*data_[0][2] - data_[2][2]*data_[0][1];
                m(1, 1) = data_[2][2]*data_[0][0] - data_[2][0]*data_[0][2];
                m(2, 1) = data_[2][0]*data_[0][1] - data_[2][1]*data_[0][0];
                m(0, 2) = data_[0][1]*data_[1][2] - data_[0][2]*data_[1][1];
                m(1, 2) = data_[0][2]*data_[1][0] - data_[0][0]*data_[1][2];
                m(2, 2) = data_[0][0]*data_[1][1] - data_[0][1]*data_[1][0];

                m /= det;
                return m;
            }

            void diagonalize(SquareMatrix3<Real>& a, Vector3<Real>& w, SquareMatrix3<Real>& v) {
                int i,j,k,maxI;
                Real tmp, maxVal;
                Vector3<Real> v_maxI, v_k, v_j;

                // diagonalize using Jacobi
                jacobi(a, w, v);

                // if all the eigenvalues are the same, return identity matrix
                if (w[0] == w[1] && w[0] == w[2] ){
                      v = SquareMatrix3<Real>::identity();
                      return
                }

                // transpose temporarily, it makes it easier to sort the eigenvectors
                v = v.tanspose(); 
                
                // if two eigenvalues are the same, re-orthogonalize to optimally line
                // up the eigenvectors with the x, y, and z axes
                for (i = 0; i < 3; i++) {
                    if (w((i+1)%3) == w((i+2)%3)) {// two eigenvalues are the same
                    // find maximum element of the independant eigenvector
                    maxVal = fabs(v(i, 0));
                    maxI = 0;
                    for (j = 1; j < 3; j++) {
                        if (maxVal < (tmp = fabs(v(i, j)))){
                            maxVal = tmp;
                            maxI = j;
                        }
                    }
                    
                    // swap the eigenvector into its proper position
                    if (maxI != i) {
                        tmp = w(maxI);
                        w(maxI) = w(i);
                        w(i) = tmp;

                        v.swapRow(i, maxI);
                    }
                    // maximum element of eigenvector should be positive
                    if (v(maxI, maxI) < 0) {
                        v(maxI, 0) = -v(maxI, 0);
                        v(maxI, 1) = -v(maxI, 1);
                        v(maxI, 2) = -v(maxI, 2);
                    }

                    // re-orthogonalize the other two eigenvectors
                    j = (maxI+1)%3;
                    k = (maxI+2)%3;

                    v(j, 0) = 0.0; 
                    v(j, 1) = 0.0; 
                    v(j, 2) = 0.0;
                    v(j, j) = 1.0;

                    /** @todo */
                    v_maxI = v.getRow(maxI);
                    v_j = v.getRow(j);
                    v_k = cross(v_maxI, v_j);
                    v_k.normailze();
                    v_j = cross(v_k, v_maxI);
                    v.setRow(j, v_j);
                    v.setRow(k, v_k);


                    // transpose vectors back to columns
                    v = v.transpose();
                    return;
                    }
                }

                // the three eigenvalues are different, just sort the eigenvectors
                // to align them with the x, y, and z axes

                // find the vector with the largest x element, make that vector
                // the first vector
                maxVal = fabs(v(0, 0));
                maxI = 0;
                for (i = 1; i < 3; i++) {
                    if (maxVal < (tmp = fabs(v(i, 0)))) {
                        maxVal = tmp;
                        maxI = i;
                    }
                }

                // swap eigenvalue and eigenvector
                if (maxI != 0) {
                    tmp = w(maxI);
                    w(maxI) = w(0);
                    w(0) = tmp;
                    v.swapRow(maxI, 0);
                }
                // do the same for the y element
                if (fabs(v(1, 1)) < fabs(v(2, 1))) {
                    tmp = w(2);
                    w(2) = w(1);
                    w(1) = tmp;
                    v.swapRow(2, 1);
                }

                // ensure that the sign of the eigenvectors is correct
                for (i = 0; i < 2; i++) {
                    if (v(i, i) < 0) {
                        v(i, 0) = -v(i, 0);
                        v(i, 1) = -v(i, 1);
                        v(i, 2) = -v(i, 2);
                    }
                }

                // set sign of final eigenvector to ensure that determinant is positive
                if (determinant(v) < 0) {
                    v(2, 0) = -v(2, 0);
                    v(2, 1) = -v(2, 1);
                    v(2, 2) = -v(2, 2);
                }

                // transpose the eigenvectors back again
                v = v.transpose();
                return ;
            }
    };

    typedef SquareMatrix3<double> Mat3x3d;
    typedef SquareMatrix3<double> RotMat3x3d;

} //namespace oopse
#endif // MATH_SQUAREMATRIX_HPP
Revision:	101
Committed:	Mon Oct 18 23:13:23 2004 UTC (21 years ago) by tim
File size:	15645 byte(s)
Log Message:	more tests on math library
#	User	Rev	Content
1	tim	70	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.org
5			*
6			* This program is free software; you can redistribute it and/or
7			* modify it under the terms of the GNU Lesser General Public License
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
22			* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23			*
24			*/
25
26			/**
27			* @file SquareMatrix3.hpp
28			* @author Teng Lin
29			* @date 10/11/2004
30			* @version 1.0
31			*/
32	tim	99	#ifndef MATH_SQUAREMATRIX3_HPP
33			#define MATH_SQUAREMATRIX3_HPP
34	tim	70
35	tim	93	#include "Quaternion.hpp"
36	tim	70	#include "SquareMatrix.hpp"
37	tim	93	#include "Vector3.hpp"
38
39	tim	70	namespace oopse {
40
41			template<typename Real>
42			class SquareMatrix3 : public SquareMatrix<Real, 3> {
43			public:
44
45			/** default constructor */
46			SquareMatrix3() : SquareMatrix<Real, 3>() {
47			}
48
49			/** copy constructor */
50			SquareMatrix3(const SquareMatrix<Real, 3>& m) : SquareMatrix<Real, 3>(m) {
51			}
52
53	tim	93	SquareMatrix3( const Vector3<Real>& eulerAngles) {
54			setupRotMat(eulerAngles);
55			}
56
57			SquareMatrix3(Real phi, Real theta, Real psi) {
58			setupRotMat(phi, theta, psi);
59			}
60
61			SquareMatrix3(const Quaternion<Real>& q) {
62			*this = q.toRotationMatrix3();
63			}
64
65			SquareMatrix3(Real w, Real x, Real y, Real z) {
66			Quaternion<Real> q(w, x, y, z);
67			*this = q.toRotationMatrix3();
68			}
69
70	tim	70	/** copy assignment operator */
71			SquareMatrix3<Real>& operator =(const SquareMatrix<Real, 3>& m) {
72			if (this == &m)
73			return *this;
74			SquareMatrix<Real, 3>::operator=(m);
75	tim	101	return *this;
76	tim	70	}
77	tim	76
78			/**
79			* Sets this matrix to a rotation matrix by three euler angles
80			* @ param euler
81			*/
82	tim	93	void setupRotMat(const Vector3<Real>& eulerAngles) {
83			setupRotMat(eulerAngles[0], eulerAngles[1], eulerAngles[2]);
84			}
85	tim	76
86			/**
87			* Sets this matrix to a rotation matrix by three euler angles
88			* @param phi
89			* @param theta
90			* @psi theta
91			*/
92	tim	93	void setupRotMat(Real phi, Real theta, Real psi) {
93			Real sphi, stheta, spsi;
94			Real cphi, ctheta, cpsi;
95	tim	76
96	tim	93	sphi = sin(phi);
97			stheta = sin(theta);
98			spsi = sin(psi);
99			cphi = cos(phi);
100			ctheta = cos(theta);
101			cpsi = cos(psi);
102	tim	76
103	tim	93	data_[0][0] = cpsi * cphi - ctheta * sphi * spsi;
104			data_[0][1] = cpsi * sphi + ctheta * cphi * spsi;
105			data_[0][2] = spsi * stheta;
106
107			data_[1][0] = -spsi * ctheta - ctheta * sphi * cpsi;
108			data_[1][1] = -spsi * stheta + ctheta * cphi * cpsi;
109			data_[1][2] = cpsi * stheta;
110
111			data_[2][0] = stheta * sphi;
112			data_[2][1] = -stheta * cphi;
113			data_[2][2] = ctheta;
114			}
115
116
117	tim	76	/**
118			* Sets this matrix to a rotation matrix by quaternion
119			* @param quat
120			*/
121	tim	93	void setupRotMat(const Quaternion<Real>& quat) {
122			*this = quat.toRotationMatrix3();
123			}
124	tim	76
125			/**
126			* Sets this matrix to a rotation matrix by quaternion
127	tim	93	* @param w the first element
128			* @param x the second element
129			* @param y the third element
130	tim	101	* @param z the fourth element
131	tim	76	*/
132	tim	93	void setupRotMat(Real w, Real x, Real y, Real z) {
133			Quaternion<Real> q(w, x, y, z);
134			*this = q.toRotationMatrix3();
135			}
136	tim	76
137			/**
138			* Returns the quaternion from this rotation matrix
139			* @return the quaternion from this rotation matrix
140			* @exception invalid rotation matrix
141			*/
142	tim	93	Quaternion<Real> toQuaternion() {
143			Quaternion<Real> q;
144			Real t, s;
145			Real ad1, ad2, ad3;
146			t = data_[0][0] + data_[1][1] + data_[2][2] + 1.0;
147	tim	76
148	tim	93	if( t > 0.0 ){
149
150			s = 0.5 / sqrt( t );
151			q[0] = 0.25 / s;
152			q[1] = (data_[1][2] - data_[2][1]) * s;
153			q[2] = (data_[2][0] - data_[0][2]) * s;
154			q[3] = (data_[0][1] - data_[1][0]) * s;
155			} else {
156
157			ad1 = fabs( data_[0][0] );
158			ad2 = fabs( data_[1][1] );
159			ad3 = fabs( data_[2][2] );
160
161			if( ad1 >= ad2 && ad1 >= ad3 ){
162
163			s = 2.0 * sqrt( 1.0 + data_[0][0] - data_[1][1] - data_[2][2] );
164			q[0] = (data_[1][2] + data_[2][1]) / s;
165			q[1] = 0.5 / s;
166			q[2] = (data_[0][1] + data_[1][0]) / s;
167			q[3] = (data_[0][2] + data_[2][0]) / s;
168			} else if ( ad2 >= ad1 && ad2 >= ad3 ) {
169			s = sqrt( 1.0 + data_[1][1] - data_[0][0] - data_[2][2] ) * 2.0;
170			q[0] = (data_[0][2] + data_[2][0]) / s;
171			q[1] = (data_[0][1] + data_[1][0]) / s;
172			q[2] = 0.5 / s;
173			q[3] = (data_[1][2] + data_[2][1]) / s;
174			} else {
175
176			s = sqrt( 1.0 + data_[2][2] - data_[0][0] - data_[1][1] ) * 2.0;
177			q[0] = (data_[0][1] + data_[1][0]) / s;
178			q[1] = (data_[0][2] + data_[2][0]) / s;
179			q[2] = (data_[1][2] + data_[2][1]) / s;
180			q[3] = 0.5 / s;
181			}
182			}
183
184			return q;
185
186			}
187
188	tim	76	/**
189			* Returns the euler angles from this rotation matrix
190	tim	93	* @return the euler angles in a vector
191	tim	76	* @exception invalid rotation matrix
192	tim	93	* We use so-called "x-convention", which is the most common definition.
193			* In this convention, the rotation given by Euler angles (phi, theta, psi), where the first
194			* rotation is by an angle phi about the z-axis, the second is by an angle
195			* theta (0 <= theta <= 180)about the x-axis, and thethird is by an angle psi about the
196			* z-axis (again).
197	tim	76	*/
198	tim	93	Vector3<Real> toEulerAngles() {
199			Vector<Real> myEuler;
200			Real phi,theta,psi,eps;
201			Real ctheta,stheta;
202
203			// set the tolerance for Euler angles and rotation elements
204
205			theta = acos(min(1.0,max(-1.0,data_[2][2])));
206			ctheta = data_[2][2];
207			stheta = sqrt(1.0 - ctheta * ctheta);
208
209			// when sin(theta) is close to 0, we need to consider singularity
210			// In this case, we can assign an arbitary value to phi (or psi), and then determine
211			// the psi (or phi) or vice-versa. We'll assume that phi always gets the rotation, and psi is 0
212			// in cases of singularity.
213			// we use atan2 instead of atan, since atan2 will give us -Pi to Pi.
214			// Since 0 <= theta <= 180, sin(theta) will be always non-negative. Therefore, it never
215			// change the sign of both of the parameters passed to atan2.
216
217			if (fabs(stheta) <= oopse::epsilon){
218			psi = 0.0;
219			phi = atan2(-data_[1][0], data_[0][0]);
220			}
221			// we only have one unique solution
222			else{
223			phi = atan2(data_[2][0], -data_[2][1]);
224			psi = atan2(data_[0][2], data_[1][2]);
225			}
226
227			//wrap phi and psi, make sure they are in the range from 0 to 2*Pi
228			if (phi < 0)
229			phi += M_PI;
230
231			if (psi < 0)
232			psi += M_PI;
233
234			myEuler[0] = phi;
235			myEuler[1] = theta;
236			myEuler[2] = psi;
237
238			return myEuler;
239			}
240	tim	70
241	tim	101	/** Returns the determinant of this matrix. */
242			Real determinant() const {
243			Real x,y,z;
244
245			x = data_[0][0] * (data_[1][1] * data_[2][2] - data_[1][2] * data_[2][1]);
246			y = data_[0][1] * (data_[1][2] * data_[2][0] - data_[1][0] * data_[2][2]);
247			z = data_[0][2] * (data_[1][0] * data_[2][1] - data_[1][1] * data_[2][0]);
248
249			return(x + y + z);
250			}
251
252	tim	70	/**
253			* Sets the value of this matrix to the inversion of itself.
254			* @note since simple algorithm can be applied to inverse the 3 by 3 matrix, we hide the
255			* implementation of inverse in SquareMatrix class
256			*/
257	tim	101	SquareMatrix3<Real> inverse() {
258			SquareMatrix3<Real> m;
259			double det = determinant();
260			if (fabs(det) <= oopse::epsilon) {
261			//"The method was called on a matrix with \|determinant\| <= 1e-6.",
262			//"This is a runtime or a programming error in your application.");
263			}
264	tim	70
265	tim	101	m(0, 0) = data_[1][1]data_[2][2] - data_[1][2]data_[2][1];
266			m(1, 0) = data_[1][2]data_[2][0] - data_[1][0]data_[2][2];
267			m(2, 0) = data_[1][0]data_[2][1] - data_[1][1]data_[2][0];
268			m(0, 1) = data_[2][1]data_[0][2] - data_[2][2]data_[0][1];
269			m(1, 1) = data_[2][2]data_[0][0] - data_[2][0]data_[0][2];
270			m(2, 1) = data_[2][0]data_[0][1] - data_[2][1]data_[0][0];
271			m(0, 2) = data_[0][1]data_[1][2] - data_[0][2]data_[1][1];
272			m(1, 2) = data_[0][2]data_[1][0] - data_[0][0]data_[1][2];
273			m(2, 2) = data_[0][0]data_[1][1] - data_[0][1]data_[1][0];
274
275			m /= det;
276			return m;
277	tim	99	}
278	tim	76
279	tim	101	void diagonalize(SquareMatrix3<Real>& a, Vector3<Real>& w, SquareMatrix3<Real>& v) {
280			int i,j,k,maxI;
281			Real tmp, maxVal;
282			Vector3<Real> v_maxI, v_k, v_j;
283	tim	99
284	tim	101	// diagonalize using Jacobi
285			jacobi(a, w, v);
286
287			// if all the eigenvalues are the same, return identity matrix
288			if (w[0] == w[1] && w[0] == w[2] ){
289			v = SquareMatrix3<Real>::identity();
290			return
291			}
292
293			// transpose temporarily, it makes it easier to sort the eigenvectors
294			v = v.tanspose();
295
296			// if two eigenvalues are the same, re-orthogonalize to optimally line
297			// up the eigenvectors with the x, y, and z axes
298			for (i = 0; i < 3; i++) {
299			if (w((i+1)%3) == w((i+2)%3)) {// two eigenvalues are the same
300			// find maximum element of the independant eigenvector
301			maxVal = fabs(v(i, 0));
302			maxI = 0;
303			for (j = 1; j < 3; j++) {
304			if (maxVal < (tmp = fabs(v(i, j)))){
305			maxVal = tmp;
306			maxI = j;
307			}
308			}
309
310			// swap the eigenvector into its proper position
311			if (maxI != i) {
312			tmp = w(maxI);
313			w(maxI) = w(i);
314			w(i) = tmp;
315
316			v.swapRow(i, maxI);
317			}
318			// maximum element of eigenvector should be positive
319			if (v(maxI, maxI) < 0) {
320			v(maxI, 0) = -v(maxI, 0);
321			v(maxI, 1) = -v(maxI, 1);
322			v(maxI, 2) = -v(maxI, 2);
323			}
324
325			// re-orthogonalize the other two eigenvectors
326			j = (maxI+1)%3;
327			k = (maxI+2)%3;
328
329			v(j, 0) = 0.0;
330			v(j, 1) = 0.0;
331			v(j, 2) = 0.0;
332			v(j, j) = 1.0;
333
334			/** @todo */
335			v_maxI = v.getRow(maxI);
336			v_j = v.getRow(j);
337			v_k = cross(v_maxI, v_j);
338			v_k.normailze();
339			v_j = cross(v_k, v_maxI);
340			v.setRow(j, v_j);
341			v.setRow(k, v_k);
342
343
344			// transpose vectors back to columns
345			v = v.transpose();
346			return;
347			}
348			}
349
350			// the three eigenvalues are different, just sort the eigenvectors
351			// to align them with the x, y, and z axes
352
353			// find the vector with the largest x element, make that vector
354			// the first vector
355			maxVal = fabs(v(0, 0));
356			maxI = 0;
357			for (i = 1; i < 3; i++) {
358			if (maxVal < (tmp = fabs(v(i, 0)))) {
359			maxVal = tmp;
360			maxI = i;
361			}
362			}
363
364			// swap eigenvalue and eigenvector
365			if (maxI != 0) {
366			tmp = w(maxI);
367			w(maxI) = w(0);
368			w(0) = tmp;
369			v.swapRow(maxI, 0);
370			}
371			// do the same for the y element
372			if (fabs(v(1, 1)) < fabs(v(2, 1))) {
373			tmp = w(2);
374			w(2) = w(1);
375			w(1) = tmp;
376			v.swapRow(2, 1);
377			}
378
379			// ensure that the sign of the eigenvectors is correct
380			for (i = 0; i < 2; i++) {
381			if (v(i, i) < 0) {
382			v(i, 0) = -v(i, 0);
383			v(i, 1) = -v(i, 1);
384			v(i, 2) = -v(i, 2);
385			}
386			}
387
388			// set sign of final eigenvector to ensure that determinant is positive
389			if (determinant(v) < 0) {
390			v(2, 0) = -v(2, 0);
391			v(2, 1) = -v(2, 1);
392			v(2, 2) = -v(2, 2);
393			}
394
395			// transpose the eigenvectors back again
396			v = v.transpose();
397			return ;
398	tim	99	}
399	tim	70	};
400
401	tim	99	typedef SquareMatrix3<double> Mat3x3d;
402			typedef SquareMatrix3<double> RotMat3x3d;
403	tim	93
404			} //namespace oopse
405			#endif // MATH_SQUAREMATRIX_HPP