matrix3x3.h

Go to the documentation of this file.

/**********************************************************************
matrix3x3.cpp - Handle 3D Rotation matrix.

Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison
Some portions Copyright (C) 2006 by Benoit Jacob

This file is part of the Open Babel project.
For more information, see <http://openbabel.org/>

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation version 2 of the License.

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 General Public License for more details.
***********************************************************************/

#ifndef OB_MATRIX3x3_H
#define OB_MATRIX3x3_H

#include <ostream>

#include <openbabel/math/vector3.h> // includes rand.h, which includes <math.h>
#include <openbabel/oberror.h>

#ifndef RAD_TO_DEG
#define RAD_TO_DEG (180.0/M_PI)
#endif

#ifndef DEG_TO_RAD
#define DEG_TO_RAD (M_PI/180.0)
#endif

namespace OpenBabel
{
 class OBRandom; // class introduction in rand.h

 // class introduction in matrix3x3.cpp
 class OBAPI matrix3x3
 {

 double ele[3][3];

 public:
 matrix3x3(void)
 {
 // Loops are typically unrolled and/or vectorized
 for (unsigned int i = 0; i < 3; ++i)
 for (unsigned int j = 0; j < 3; ++j)
 ele[i][j] = 0.0;
 }

 matrix3x3(double s)
 {
 // Loops are typically unrolled and/or vectorized
 for (unsigned int i = 0; i < 3; ++i)
 for (unsigned int j = 0; j < 3; ++j)
 ele[i][j] = 0.0;

 for (unsigned int i = 0; i < 3; ++i)
 ele[i][i] = s;
 }

 matrix3x3(vector3 row1,vector3 row2,vector3 row3)
 {
 ele[0][0] = row1.x();
 ele[0][1] = row1.y();
 ele[0][2] = row1.z();
 ele[1][0] = row2.x();
 ele[1][1] = row2.y();
 ele[1][2] = row2.z();
 ele[2][0] = row3.x();
 ele[2][1] = row3.y();
 ele[2][2] = row3.z();
 }


 matrix3x3(double d[3][3])
 {
 // Loops are typically unrolled and/or vectorized
 for (unsigned int i = 0; i < 3; ++i)
 for (unsigned int j = 0; j < 3; ++j)
 ele[i][j] = d[i][j];

 // We could also potentially use memcpy here
 }

 ~matrix3x3() {}


 void GetArray(double *m)
 {
 for (unsigned int i = 0; i < 3; ++i)
 for (unsigned int j = 0; j < 3; ++j)
 m[3*i+j] = ele[i][j];
 }

 const double & operator() (int row, int column ) const
 {
 return ele[row][column];
 }

 double & operator() (int row, int column )
 {
 return ele[row][column];
 }

 matrix3x3 inverse(void) const
#ifdef OB_OLD_MATH_CHECKS
 throw(OBError)
#endif
 ;

 matrix3x3 transpose(void) const;

 double determinant() const;

 bool isSymmetric(void) const;


 bool isOrthogonal(void) const
 {
 return (*this * transpose()).isUnitMatrix();
 };

 bool isDiagonal(void) const;

 bool isUnitMatrix(void) const;


 double Get(int row,int column) const
 {
#ifdef OB_OLD_MATH_CHECKS
 if (row >= 0 && row <= 2 && column >= 0 && column <= 2)
 return(ele[row][column]);
 else
 return 0.0f;
#else
 return(ele[row][column]);
#endif
 }


 void Set(int row,int column, double v)
 {
#ifdef OB_OLD_MATH_CHECKS
 if (row >= 0 && row <= 2 && column >= 0 && column <= 2)
 ele[row][column]= v;
#else
 ele[row][column]= v;
#endif
 }


 void SetColumn(int column, const vector3 &v)
#ifdef OB_OLD_MATH_CHECKS
 throw(OBError)
#endif
 ;


 void SetRow(int row, const vector3 &v)
#ifdef OB_OLD_MATH_CHECKS
 throw(OBError)
#endif
 ;


 vector3 GetColumn(unsigned int col) const
#ifdef OB_OLD_MATH_CHECKS
 throw(OBError)
#endif
 ;


 vector3 GetRow(unsigned int row) const
#ifdef OB_OLD_MATH_CHECKS
 throw(OBError)
#endif
 ;

 matrix3x3 &operator*=(const double &c)
 {
 for( int i = 0; i < 3; i++ )
 for( int j = 0; j < 3; j++ )
 ele[i][j] *= c;
 return *this;
 }

 matrix3x3 &operator/=(const double &c)
 {
 return( (*this) *= ( 1.0 / c ) );
 }

 void SetupRotMat(double x, double y, double z);

 void PlaneReflection(const vector3 &norm);

 void RotAboutAxisByAngle(const vector3 &axis, const double angle);

 void FillOrth(double alpha, double beta, double gamma,
 double a, double b, double c);

 matrix3x3 findEigenvectorsIfSymmetric(vector3 &eigenvals) const
#ifdef OB_OLD_MATH_CHECKS
 throw(OBError)
#endif
 ;

 friend OBAPI vector3 operator *(const matrix3x3 &,const vector3 &);

 friend OBAPI matrix3x3 operator *(const matrix3x3 &,const matrix3x3 &);

 friend OBAPI std::ostream& operator<< ( std::ostream&, const matrix3x3 & ) ;

 static void jacobi(unsigned int n, double *a, double *d, double *v);
 };

#ifndef SWIG
 OBAPI vector3 center_coords(double*,int);
#endif
}

#endif // OB_MATRIX3x3_H