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;

      void randomRotation(OBRandom &rnd);

      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