← Back to branch summary

~ubuntu-branches/debian/stretch/openbabel/stretch

« back to all changes in this revision

Viewing changes to src/math/matrix3x3.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-07-22 23:54:58 UTC
  • mfrom: (3.1.10 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080722235458-3o606czluviz4akx
Tags: 2.2.0-2
* Upload to unstable.
* debian/control: Updated descriptions.
* debian/patches/gauss_cube_format.patch: New patch, makes the 
  gaussian cube format available again.
* debian/rules (DEB_DH_MAKESHLIBS_ARGS_libopenbabel3): Removed.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Likewise.
* debian/libopenbabel3.install: Adjust formats directory to single 
  version hierarchy.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/math/matrix3x3.cpp
2
2
matrix3x3.cpp - Handle 3D rotation matrix.
3
3
 
4
4
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
5
 
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
6
 
 
 
5
Some portions Copyright (C) 2001-2006 by Geoffrey R. Hutchison
 
6
Some portions Copyright (C) 2006 by Benoit Jacob
 
7
 
7
8
This file is part of the Open Babel project.
8
9
For more information, see <http://openbabel.sourceforge.net/>
9
10
 
16
17
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17
18
GNU General Public License for more details.
18
19
***********************************************************************/
19
 
 
20
 
#include <math.h>
21
 
 
22
 
#include "mol.h"
23
 
#include "math/matrix3x3.h"
24
 
 
25
 
#ifndef true
26
 
#define true 1
27
 
#endif
28
 
 
29
 
#ifndef false
30
 
#define false 0
31
 
#endif
 
20
#include <openbabel/babelconfig.h>
 
21
 
 
22
#include <openbabel/math/matrix3x3.h>
 
23
#include <openbabel/obutil.h>
32
24
 
33
25
using namespace std;
34
26
 
35
27
namespace OpenBabel
36
28
{
37
29
 
38
 
/** \class matrix3x3
39
 
    \brief Represents a real 3x3 matrix.
40
 
 
41
 
 Rotating points in space can be performed by a vector-matrix
42
 
 multiplication. The matrix3x3 class is designed as a helper to the
43
 
 vector3 class for rotating points in space. The rotation matrix may be
44
 
 initialised by passing in the array of floating point values, by
45
 
 passing euler angles, or a rotation vector and angle of rotation about
46
 
 that vector. Once set, the matrix3x3 class can be used to rotate
47
 
 vectors by the overloaded multiplication operator. The following
48
 
 demonstrates the usage of the matrix3x3 class:
49
 
 
50
 
\code
51
 
  matrix3x3 mat;
52
 
  mat.SetupRotMat(0.0,180.0,0.0); //rotate theta by 180 degrees
53
 
  vector3 v = VX;
54
 
  v *= mat; //apply the rotation
55
 
\endcode
56
 
 
57
 
*/
 
30
  /** \class matrix3x3 matrix3x3.h <openbabel/math/matrix3x3.h>
 
31
      \brief Represents a real 3x3 matrix.
 
32
 
 
33
      Rotating points in space can be performed by a vector-matrix
 
34
      multiplication. The matrix3x3 class is designed as a helper to the
 
35
      vector3 class for rotating points in space. The rotation matrix may be
 
36
      initialised by passing in the array of floating point values, by
 
37
      passing euler angles, or a rotation vector and angle of rotation about
 
38
      that vector. Once set, the matrix3x3 class can be used to rotate
 
39
      vectors by the overloaded multiplication operator. The following
 
40
      demonstrates the usage of the matrix3x3 class:
 
41
 
 
42
      \code
 
43
      matrix3x3 mat;
 
44
      mat.SetupRotMat(0.0,180.0,0.0); //rotate theta by 180 degrees
 
45
      vector3 v = VX;
 
46
      v *= mat; //apply the rotation
 
47
      \endcode
 
48
 
 
49
  */
58
50
 
59
 
/*! the axis of the rotation will be uniformly distributed on
60
 
  the unit sphere, the angle will be uniformly distributed in
61
 
  the interval 0..360 degrees. */
62
 
void matrix3x3::randomRotation(OBRandom &rnd)
63
 
{
 
51
  /*! The axis of the rotation will be uniformly distributed on
 
52
    the unit sphere and the angle will be uniformly distributed in
 
53
    the interval 0..360 degrees. */
 
54
  void matrix3x3::randomRotation(OBRandom &rnd)
 
55
  {
64
56
    double rotAngle;
65
57
    vector3 v1;
66
58
 
67
59
    v1.randomUnitVector(&rnd);
68
60
    rotAngle = 360.0 * rnd.NextFloat();
69
61
    this->RotAboutAxisByAngle(v1,rotAngle);
70
 
}
 
62
  }
71
63
 
72
 
void matrix3x3::SetupRotMat(double phi,double theta,double psi)
73
 
{
 
64
  void matrix3x3::SetupRotMat(double phi,double theta,double psi)
 
65
  {
74
66
    double p  = phi * DEG_TO_RAD;
75
67
    double h  = theta * DEG_TO_RAD;
76
68
    double b  = psi * DEG_TO_RAD;
93
85
    ele[2][0] = cx*sy*cz+sx*sz;
94
86
    ele[2][1] = cx*sy*sz-sx*cz;
95
87
    ele[2][2] = cx*cy;
96
 
}
 
88
  }
97
89
 
98
 
/*! Replaces *this with a matrix that represents reflection on
99
 
  the plane through 0 which is given by the normal vector norm.
100
 
        
101
 
  \warning If the vector norm has length zero, this method will
102
 
  generate the 0-matrix. If the length of the axis is close to
103
 
  zero, but not == 0.0, this method may behave in unexpected
104
 
  ways and return almost random results; details may depend on
105
 
  your particular floating point implementation. The use of this
106
 
  method is therefore highly discouraged, unless you are certain
107
 
  that the length is in a reasonable range, away from 0.0
108
 
  (Stefan Kebekus)
109
 
        
110
 
  \deprecated This method will probably replaced by a safer
111
 
  algorithm in the future.
112
 
        
113
 
  \todo Replace this method with a more fool-proof version.
114
 
        
115
 
  @param norm specifies the normal to the plane
116
 
*/
117
 
void matrix3x3::PlaneReflection(const vector3 &norm)
118
 
{
 
90
  /*! Replaces *this with a matrix that represents reflection on
 
91
    the plane through 0 which is given by the normal vector norm.
 
92
        
 
93
    \warning If the vector norm has length zero, this method will
 
94
    generate the 0-matrix. If the length of the axis is close to
 
95
    zero, but not == 0.0, this method may behave in unexpected
 
96
    ways and return almost random results; details may depend on
 
97
    your particular floating point implementation. The use of this
 
98
    method is therefore highly discouraged, unless you are certain
 
99
    that the length is in a reasonable range, away from 0.0
 
100
    (Stefan Kebekus)
 
101
        
 
102
    \deprecated This method will probably replaced by a safer
 
103
    algorithm in the future.
 
104
        
 
105
    \todo Replace this method with a more fool-proof version.
 
106
        
 
107
    @param norm specifies the normal to the plane
 
108
  */
 
109
  void matrix3x3::PlaneReflection(const vector3 &norm)
 
110
  {
119
111
    //@@@ add a safety net
120
112
 
121
113
    vector3 normtmp = norm;
124
116
    SetColumn(0, vector3(1,0,0) - 2*normtmp.x()*normtmp);
125
117
    SetColumn(1, vector3(0,1,0) - 2*normtmp.y()*normtmp);
126
118
    SetColumn(2, vector3(0,0,1) - 2*normtmp.z()*normtmp);
127
 
}
128
 
 
129
 
#define x vtmp.x()
130
 
#define y vtmp.y()
131
 
#define z vtmp.z()
132
 
 
133
 
/*! Replaces *this with a matrix that represents rotation about the
134
 
  axis by a an angle. 
135
 
        
136
 
  \warning If the vector axis has length zero, this method will
137
 
  generate the 0-matrix. If the length of the axis is close to
138
 
  zero, but not == 0.0, this method may behave in unexpected ways
139
 
  and return almost random results; details may depend on your
140
 
  particular floating point implementation. The use of this method
141
 
  is therefore highly discouraged, unless you are certain that the
142
 
  length is in a reasonable range, away from 0.0 (Stefan
143
 
  Kebekus)
144
 
        
145
 
  \deprecated This method will probably replaced by a safer
146
 
  algorithm in the future.
147
 
        
148
 
  \todo Replace this method with a more fool-proof version.
149
 
        
150
 
  @param v specifies the axis of the rotation
151
 
  @param angle angle in degrees (0..360)
152
 
*/
153
 
void matrix3x3::RotAboutAxisByAngle(const vector3 &v,const double angle)
154
 
{
 
119
  }
 
120
 
 
121
  /*! Replaces *this with a matrix that represents rotation about the
 
122
    axis by a an angle. 
 
123
        
 
124
    \warning If the vector axis has length zero, this method will
 
125
    generate the 0-matrix. If the length of the axis is close to
 
126
    zero, but not == 0.0, this method may behave in unexpected ways
 
127
    and return almost random results; details may depend on your
 
128
    particular floating point implementation. The use of this method
 
129
    is therefore highly discouraged, unless you are certain that the
 
130
    length is in a reasonable range, away from 0.0 (Stefan
 
131
    Kebekus)
 
132
        
 
133
    \deprecated This method will probably replaced by a safer
 
134
    algorithm in the future.
 
135
        
 
136
    \todo Replace this method with a more fool-proof version.
 
137
        
 
138
    @param v specifies the axis of the rotation
 
139
    @param angle angle in degrees (0..360)
 
140
  */
 
141
  void matrix3x3::RotAboutAxisByAngle(const vector3 &v,const double angle)
 
142
  {
155
143
    double theta = angle*DEG_TO_RAD;
156
144
    double s = sin(theta);
157
145
    double c = cos(theta);
160
148
    vector3 vtmp = v;
161
149
    vtmp.normalize();
162
150
 
 
151
    double x = vtmp.x(),
 
152
           y = vtmp.y(),
 
153
           z = vtmp.z();
 
154
 
163
155
    ele[0][0] = t*x*x + c;
164
156
    ele[0][1] = t*x*y + s*z;
165
157
    ele[0][2] = t*x*z - s*y;
171
163
    ele[2][0] = t*x*z + s*y;
172
164
    ele[2][1] = t*y*z - s*x;
173
165
    ele[2][2] = t*z*z + c;
174
 
}
175
 
 
176
 
#undef x
177
 
#undef y
178
 
#undef z
179
 
 
180
 
void matrix3x3::SetColumn(int col, const vector3 &v) throw(OBError)
181
 
{
 
166
  }
 
167
 
 
168
  void matrix3x3::SetColumn(int col, const vector3 &v)
 
169
#ifdef OB_OLD_MATH_CHECKS
 
170
  throw(OBError)
 
171
#endif
 
172
  {
 
173
#ifdef OB_OLD_MATH_CHECKS
182
174
    if (col > 2)
183
 
    {
 
175
      {
184
176
        OBError er("matrix3x3::SetColumn(int col, const vector3 &v)",
185
177
                   "The method was called with col > 2.",
186
178
                   "This is a programming error in your application.");
187
179
        throw er;
188
 
    }
 
180
      }
 
181
#endif
189
182
 
190
183
    ele[0][col] = v.x();
191
184
    ele[1][col] = v.y();
192
185
    ele[2][col] = v.z();
193
 
}
 
186
  }
194
187
 
195
 
void matrix3x3::SetRow(int row, const vector3 &v) throw(OBError)
196
 
{
 
188
  void matrix3x3::SetRow(int row, const vector3 &v)
 
189
#ifdef OB_OLD_MATH_CHECKS
 
190
  throw(OBError)
 
191
#endif
 
192
  {
 
193
#ifdef OB_OLD_MATH_CHECKS
197
194
    if (row > 2)
198
 
    {
 
195
      {
199
196
        OBError er("matrix3x3::SetRow(int row, const vector3 &v)",
200
197
                   "The method was called with row > 2.",
201
198
                   "This is a programming error in your application.");
202
199
        throw er;
203
 
    }
 
200
      }
 
201
#endif
204
202
 
205
203
    ele[row][0] = v.x();
206
204
    ele[row][1] = v.y();
207
205
    ele[row][2] = v.z();
208
 
}
 
206
  }
209
207
 
210
 
vector3 matrix3x3::GetColumn(unsigned int col) const throw(OBError)
211
 
{
 
208
  vector3 matrix3x3::GetColumn(unsigned int col) const
 
209
#ifdef OB_OLD_MATH_CHECKS
 
210
  throw(OBError)
 
211
#endif
 
212
  {
 
213
#ifdef OB_OLD_MATH_CHECKS
212
214
    if (col > 2)
213
 
    {
 
215
      {
214
216
        OBError er("matrix3x3::GetColumn(unsigned int col) const",
215
217
                   "The method was called with col > 2.",
216
218
                   "This is a programming error in your application.");
217
219
        throw er;
218
 
    }
 
220
      }
 
221
#endif
219
222
 
220
223
    return vector3(ele[0][col], ele[1][col], ele[2][col]);
221
 
}
 
224
  }
222
225
 
223
 
vector3 matrix3x3::GetRow(unsigned int row) const throw(OBError)
224
 
{
 
226
  vector3 matrix3x3::GetRow(unsigned int row) const
 
227
#ifdef OB_OLD_MATH_CHECKS
 
228
  throw(OBError)
 
229
#endif
 
230
  {
 
231
#ifdef OB_OLD_MATH_CHECKS
225
232
    if (row > 2)
226
 
    {
 
233
      {
227
234
        OBError er("matrix3x3::GetRow(unsigned int row) const",
228
235
                   "The method was called with row > 2.",
229
236
                   "This is a programming error in your application.");
230
237
        throw er;
231
 
    }
 
238
      }
 
239
#endif
232
240
 
233
241
    return vector3(ele[row][0], ele[row][1], ele[row][2]);
234
 
}
 
242
  }
235
243
 
236
 
/*! calculates the product m*v of the matrix m and the column
237
 
  vector represented by v
238
 
*/
239
 
vector3 operator *(const matrix3x3 &m,const vector3 &v)
240
 
{
 
244
  /*! Calculates the product m*v of the matrix m and the column
 
245
    vector represented by v
 
246
  */
 
247
  vector3 operator *(const matrix3x3 &m,const vector3 &v)
 
248
  {
241
249
    vector3 vv;
242
250
 
243
 
    vv._vx = v._vx*m.ele[0][0] + v._vy*m.ele[0][1] + v._vz*m.ele[0][2];
244
 
    vv._vy = v._vx*m.ele[1][0] + v._vy*m.ele[1][1] + v._vz*m.ele[1][2];
245
 
    vv._vz = v._vx*m.ele[2][0] + v._vy*m.ele[2][1] + v._vz*m.ele[2][2];
 
251
    vv.x() = v.x()*m.ele[0][0] + v.y()*m.ele[0][1] + v.z()*m.ele[0][2];
 
252
    vv.y() = v.x()*m.ele[1][0] + v.y()*m.ele[1][1] + v.z()*m.ele[1][2];
 
253
    vv.z() = v.x()*m.ele[2][0] + v.y()*m.ele[2][1] + v.z()*m.ele[2][2];
246
254
 
247
255
    return(vv);
248
 
}
 
256
  }
249
257
 
250
 
matrix3x3 operator *(const matrix3x3 &A,const matrix3x3 &B)
251
 
{
 
258
  matrix3x3 operator *(const matrix3x3 &A,const matrix3x3 &B)
 
259
  {
252
260
    matrix3x3 result;
253
261
 
254
262
    result.ele[0][0] = A.ele[0][0]*B.ele[0][0] + A.ele[0][1]*B.ele[1][0] + A.ele[0][2]*B.ele[2][0];
264
272
    result.ele[2][2] = A.ele[2][0]*B.ele[0][2] + A.ele[2][1]*B.ele[1][2] + A.ele[2][2]*B.ele[2][2];
265
273
 
266
274
    return(result);
267
 
}
 
275
  }
268
276
 
269
 
/*! calculates the product m*(*this) of the matrix m and the
270
 
  column vector represented by *this
271
 
*/
272
 
vector3 &vector3::operator *= (const matrix3x3 &m)
273
 
{
 
277
  /*! calculates the product m*(*this) of the matrix m and the
 
278
    column vector represented by *this
 
279
  */
 
280
  vector3 &vector3::operator *= (const matrix3x3 &m)
 
281
  {
274
282
    vector3 vv;
275
283
 
276
 
    vv.SetX(_vx*m.Get(0,0) + _vy*m.Get(0,1) + _vz*m.Get(0,2));
277
 
    vv.SetY(_vx*m.Get(1,0) + _vy*m.Get(1,1) + _vz*m.Get(1,2));
278
 
    vv.SetZ(_vx*m.Get(2,0) + _vy*m.Get(2,1) + _vz*m.Get(2,2));
279
 
    _vx = vv.x();
280
 
    _vy = vv.y();
281
 
    _vz = vv.z();
 
284
    vv.SetX(x()*m.Get(0,0) + y()*m.Get(0,1) + z()*m.Get(0,2));
 
285
    vv.SetY(x()*m.Get(1,0) + y()*m.Get(1,1) + z()*m.Get(1,2));
 
286
    vv.SetZ(x()*m.Get(2,0) + y()*m.Get(2,1) + z()*m.Get(2,2));
 
287
    x() = vv.x();
 
288
    y() = vv.y();
 
289
    z() = vv.z();
282
290
 
283
291
    return(*this);
284
 
}
 
292
  }
285
293
 
286
 
/*! This method checks if the absolute value of the determinant is smaller than 1e-6. If
287
 
  so, nothing is done and an exception is thrown. Otherwise, the
288
 
  inverse matrix is calculated and returned. *this is not changed.
 
294
  /*! This method checks if the absolute value of the determinant is smaller than 1e-6. If
 
295
    so, nothing is done and an exception is thrown. Otherwise, the
 
296
    inverse matrix is calculated and returned. *this is not changed.
289
297
        
290
 
  \warning If the determinant is close to zero, but not == 0.0,
291
 
  this method may behave in unexpected ways and return almost
292
 
  random results; details may depend on your particular floating
293
 
  point implementation. The use of this method is therefore highly
294
 
  discouraged, unless you are certain that the determinant is in a
295
 
  reasonable range, away from 0.0 (Stefan Kebekus)
296
 
*/
297
 
matrix3x3 matrix3x3::inverse(void) const throw(OBError)
298
 
{
 
298
    \warning If the determinant is close to zero, but not == 0.0,
 
299
    this method may behave in unexpected ways and return almost
 
300
    random results; details may depend on your particular floating
 
301
    point implementation. The use of this method is therefore highly
 
302
    discouraged, unless you are certain that the determinant is in a
 
303
    reasonable range, away from 0.0 (Stefan Kebekus)
 
304
  */
 
305
  matrix3x3 matrix3x3::inverse(void) const
 
306
#ifdef OB_OLD_MATH_CHECKS
 
307
  throw(OBError)
 
308
#endif
 
309
  {
299
310
    double det = determinant();
 
311
 
 
312
#ifdef OB_OLD_MATH_CHECKS
300
313
    if (fabs(det) <= 1e-6)
301
 
    {
 
314
      {
302
315
        OBError er("matrix3x3::invert(void)",
303
316
                   "The method was called on a matrix with |determinant| <= 1e-6.",
304
317
                   "This is a runtime or a programming error in your application.");
305
318
        throw er;
306
 
    }
307
 
 
308
 
    matrix3x3 inverse;
309
 
    inverse.ele[0][0] = ele[1][1]*ele[2][2] - ele[1][2]*ele[2][1];
310
 
    inverse.ele[1][0] = ele[1][2]*ele[2][0] - ele[1][0]*ele[2][2];
311
 
    inverse.ele[2][0] = ele[1][0]*ele[2][1] - ele[1][1]*ele[2][0];
312
 
    inverse.ele[0][1] = ele[2][1]*ele[0][2] - ele[2][2]*ele[0][1];
313
 
    inverse.ele[1][1] = ele[2][2]*ele[0][0] - ele[2][0]*ele[0][2];
314
 
    inverse.ele[2][1] = ele[2][0]*ele[0][1] - ele[2][1]*ele[0][0];
315
 
    inverse.ele[0][2] = ele[0][1]*ele[1][2] - ele[0][2]*ele[1][1];
316
 
    inverse.ele[1][2] = ele[0][2]*ele[1][0] - ele[0][0]*ele[1][2];
317
 
    inverse.ele[2][2] = ele[0][0]*ele[1][1] - ele[0][1]*ele[1][0];
318
 
 
319
 
    inverse /= det;
320
 
 
321
 
    return(inverse);
322
 
}
323
 
 
324
 
/* This method returns the transpose of a matrix. The original
325
 
   matrix remains unchanged. */
326
 
matrix3x3 matrix3x3::transpose(void) const
327
 
{
328
 
    matrix3x3 transpose;
 
319
      }
 
320
#endif
 
321
 
 
322
    matrix3x3 returnValue;
 
323
    returnValue.ele[0][0] = ele[1][1]*ele[2][2] - ele[1][2]*ele[2][1];
 
324
    returnValue.ele[1][0] = ele[1][2]*ele[2][0] - ele[1][0]*ele[2][2];
 
325
    returnValue.ele[2][0] = ele[1][0]*ele[2][1] - ele[1][1]*ele[2][0];
 
326
    returnValue.ele[0][1] = ele[2][1]*ele[0][2] - ele[2][2]*ele[0][1];
 
327
    returnValue.ele[1][1] = ele[2][2]*ele[0][0] - ele[2][0]*ele[0][2];
 
328
    returnValue.ele[2][1] = ele[2][0]*ele[0][1] - ele[2][1]*ele[0][0];
 
329
    returnValue.ele[0][2] = ele[0][1]*ele[1][2] - ele[0][2]*ele[1][1];
 
330
    returnValue.ele[1][2] = ele[0][2]*ele[1][0] - ele[0][0]*ele[1][2];
 
331
    returnValue.ele[2][2] = ele[0][0]*ele[1][1] - ele[0][1]*ele[1][0];
 
332
 
 
333
    returnValue /= det;
 
334
 
 
335
    return(returnValue);
 
336
  }
 
337
 
 
338
  /* This method returns the transpose of a matrix. The original
 
339
     matrix remains unchanged. */
 
340
  matrix3x3 matrix3x3::transpose(void) const
 
341
  {
 
342
    matrix3x3 returnValue;
329
343
 
330
344
    for(unsigned int i=0; i<3; i++)
331
 
        for(unsigned int j=0; j<3; j++)
332
 
            transpose.ele[i][j] = ele[j][i];
333
 
 
334
 
    return(transpose);
335
 
}
336
 
 
337
 
double matrix3x3::determinant(void) const
338
 
{
339
 
    double x,y,z;
340
 
 
341
 
    x = ele[0][0] * (ele[1][1] * ele[2][2] - ele[1][2] * ele[2][1]);
342
 
    y = ele[0][1] * (ele[1][2] * ele[2][0] - ele[1][0] * ele[2][2]);
343
 
    z = ele[0][2] * (ele[1][0] * ele[2][1] - ele[1][1] * ele[2][0]);
344
 
 
345
 
    return(x + y + z);
346
 
}
347
 
 
348
 
/*! This method returns false if there are indices i,j such that
349
 
  fabs(*this[i][j]-*this[j][i]) > 1e-6. Otherwise, it returns
350
 
  true. */
351
 
bool matrix3x3::isSymmetric(void) const
352
 
{
353
 
    if (fabs(ele[0][1] - ele[1][0]) > 1e-6)
354
 
        return false;
355
 
    if (fabs(ele[0][2] - ele[2][0]) > 1e-6)
356
 
        return false;
357
 
    if (fabs(ele[1][2] - ele[2][1]) > 1e-6)
358
 
        return false;
359
 
    return true;
360
 
}
361
 
 
362
 
/*! This method returns false if there are indices i != j such
363
 
  that fabs(*this[i][j]) > 1e-6. Otherwise, it returns true. */
364
 
bool matrix3x3::isDiagonal(void) const
365
 
{
366
 
    if (fabs(ele[0][1]) > 1e-6)
367
 
        return false;
368
 
    if (fabs(ele[0][2]) > 1e-6)
369
 
        return false;
370
 
    if (fabs(ele[1][2]) > 1e-6)
371
 
        return false;
372
 
 
373
 
    if (fabs(ele[1][0]) > 1e-6)
374
 
        return false;
375
 
    if (fabs(ele[2][0]) > 1e-6)
376
 
        return false;
377
 
    if (fabs(ele[2][1]) > 1e-6)
378
 
        return false;
379
 
 
380
 
    return true;
381
 
}
382
 
 
383
 
/*! This method returns false if there are indices i != j such
384
 
  that fabs(*this[i][j]) > 1e-6, or if there is an index i such
385
 
  that fabs(*this[i][j]-1) > 1e-6. Otherwise, it returns
386
 
  true. */
387
 
bool matrix3x3::isUnitMatrix(void) const
388
 
{
389
 
    if (!isDiagonal())
390
 
        return false;
391
 
 
392
 
    if (fabs(ele[0][0]-1) > 1e-6)
393
 
        return false;
394
 
    if (fabs(ele[1][1]-1) > 1e-6)
395
 
        return false;
396
 
    if (fabs(ele[2][2]-1) > 1e-6)
397
 
        return false;
398
 
 
399
 
    return true;
400
 
}
401
 
 
402
 
/*! This method employs the static method matrix3x3::jacobi(...)
403
 
  to find the eigenvalues and eigenvectors of a symmetric
404
 
  matrix. On entry it is checked if the matrix really is
405
 
  symmetric: if isSymmetric() returns 'false', an OBError is
406
 
  thrown.
407
 
 
408
 
  \note The jacobi algorithm is should work great for all
409
 
  symmetric 3x3 matrices. If you need to find the eigenvectors
410
 
  of a non-symmetric matrix, you might want to resort to the
411
 
  sophisticated routines of LAPACK.
412
 
 
413
 
  @param eigenvals a reference to a vector3 where the
414
 
  eigenvalues will be stored. The eigenvalues are ordered so
415
 
  that eigenvals[0] <= eigenvals[1] <= eigenvals[2].
416
 
 
417
 
  @return an orthogonal matrix whose ith column is an
418
 
  eigenvector for the eigenvalue eigenvals[i]. Here 'orthogonal'
419
 
  means that all eigenvectors have length one and are mutually
420
 
  orthogonal. The ith eigenvector can thus be conveniently
421
 
  accessed by the GetColumn() method, as in the following
422
 
  example.
423
 
  \code
424
 
  // Calculate eigenvectors and -values
425
 
  vector3 eigenvals;
426
 
  matrix3x3 eigenmatrix = somematrix.findEigenvectorsIfSymmetric(eigenvals);
427
 
  
428
 
  // Print the 2nd eigenvector
429
 
  cout << eigenmatrix.GetColumn(1) << endl;
430
 
  \endcode
431
 
  With these conventions, a matrix is diagonalized in the following way:
432
 
  \code
433
 
  // Diagonalize the matrix
434
 
  matrix3x3 diagonalMatrix = eigenmatrix.inverse() * somematrix * eigenmatrix;
435
 
  \endcode
436
 
  
437
 
*/
438
 
matrix3x3 matrix3x3::findEigenvectorsIfSymmetric(vector3 &eigenvals) const throw(OBError)
439
 
{
 
345
      for(unsigned int j=0; j<3; j++)
 
346
        returnValue.ele[i][j] = ele[j][i];
 
347
 
 
348
    return(returnValue);
 
349
  }
 
350
 
 
351
  double matrix3x3::determinant(void) const
 
352
  {
 
353
    return( ele[0][0] * (ele[1][1] * ele[2][2] - ele[1][2] * ele[2][1])
 
354
          + ele[0][1] * (ele[1][2] * ele[2][0] - ele[1][0] * ele[2][2])
 
355
          + ele[0][2] * (ele[1][0] * ele[2][1] - ele[1][1] * ele[2][0]) );
 
356
  }
 
357
 
 
358
  /*! \return False if there are indices i,j such that
 
359
    fabs(*this[i][j]-*this[j][i]) > 1e-6. Otherwise, it returns
 
360
    true. */
 
361
  bool matrix3x3::isSymmetric(void) const
 
362
  {
 
363
    return( IsApprox( ele[0][1], ele[1][0], 1e-6 )
 
364
         && IsApprox( ele[0][2], ele[2][0], 1e-6 )
 
365
         && IsApprox( ele[1][2], ele[2][1], 1e-6 ) );
 
366
  }
 
367
 
 
368
  /*! This method returns true if and only if the matrix is
 
369
   * (approximately) a diagonal matrix. The precision used
 
370
   * by this function is 1e-6. */
 
371
  bool matrix3x3::isDiagonal(void) const
 
372
  {
 
373
    return( IsNegligible( ele[1][0], ele[0][0], 1e-6 )
 
374
         && IsNegligible( ele[2][0], ele[0][0], 1e-6 )
 
375
         && IsNegligible( ele[0][1], ele[1][1], 1e-6 )
 
376
         && IsNegligible( ele[2][1], ele[1][1], 1e-6 )
 
377
         && IsNegligible( ele[0][2], ele[2][2], 1e-6 )
 
378
         && IsNegligible( ele[1][2], ele[2][2], 1e-6 ) );
 
379
  }
 
380
 
 
381
  /*! This method returns true if and only if the matrix is
 
382
   * (approximately) equal to the identity matrix. The precision used
 
383
   * by this function is 1e-6. */
 
384
  bool matrix3x3::isUnitMatrix(void) const
 
385
  {
 
386
    return ( isDiagonal()
 
387
          && IsApprox( ele[0][0], 1.0, 1e-6 )
 
388
          && IsApprox( ele[1][1], 1.0, 1e-6 )
 
389
          && IsApprox( ele[2][2], 1.0, 1e-6 ) );
 
390
  }
 
391
 
 
392
  /*! This method employs the static method matrix3x3::jacobi(...)
 
393
    to find the eigenvalues and eigenvectors of a symmetric
 
394
    matrix. On entry it is checked if the matrix really is
 
395
    symmetric: if isSymmetric() returns 'false', an OBError is
 
396
    thrown.
 
397
 
 
398
    \note The jacobi algorithm is should work great for all
 
399
    symmetric 3x3 matrices. If you need to find the eigenvectors
 
400
    of a non-symmetric matrix, you might want to resort to the
 
401
    sophisticated routines of LAPACK.
 
402
 
 
403
    @param eigenvals a reference to a vector3 where the
 
404
    eigenvalues will be stored. The eigenvalues are ordered so
 
405
    that eigenvals[0] <= eigenvals[1] <= eigenvals[2].
 
406
 
 
407
    @return an orthogonal matrix whose ith column is an
 
408
    eigenvector for the eigenvalue eigenvals[i]. Here 'orthogonal'
 
409
    means that all eigenvectors have length one and are mutually
 
410
    orthogonal. The ith eigenvector can thus be conveniently
 
411
    accessed by the GetColumn() method, as in the following
 
412
    example.
 
413
    \code
 
414
    // Calculate eigenvectors and -values
 
415
    vector3 eigenvals;
 
416
    matrix3x3 eigenmatrix = somematrix.findEigenvectorsIfSymmetric(eigenvals);
 
417
  
 
418
    // Print the 2nd eigenvector
 
419
    cout << eigenmatrix.GetColumn(1) << endl;
 
420
    \endcode
 
421
    With these conventions, a matrix is diagonalized in the following way:
 
422
    \code
 
423
    // Diagonalize the matrix
 
424
    matrix3x3 diagonalMatrix = eigenmatrix.inverse() * somematrix * eigenmatrix;
 
425
    \endcode
 
426
  
 
427
  */
 
428
  matrix3x3 matrix3x3::findEigenvectorsIfSymmetric(vector3 &eigenvals) const
 
429
#ifdef OB_OLD_MATH_CHECKS
 
430
  throw(OBError)
 
431
#endif
 
432
  {
440
433
    matrix3x3 result;
441
434
 
 
435
#ifdef OB_OLD_MATH_CHECKS
442
436
    if (!isSymmetric())
443
 
    {
 
437
      {
444
438
        OBError er("matrix3x3::findEigenvectorsIfSymmetric(vector3 &eigenvals) const throw(OBError)",
445
439
                   "The method was called on a matrix that was not symmetric, i.e. where isSymetric() == false.",
446
440
                   "This is a runtime or a programming error in your application.");
447
441
        throw er;
448
 
    }
 
442
      }
 
443
#endif
449
444
 
450
445
    double d[3];
451
446
    matrix3x3 copyOfThis = *this;
454
449
    eigenvals.Set(d);
455
450
 
456
451
    return result;
457
 
}
458
 
 
459
 
matrix3x3 &matrix3x3::operator/=(const double &c)
460
 
{
461
 
    for (int row = 0;row < 3; row++)
462
 
        for (int col = 0;col < 3; col++)
463
 
            ele[row][col] /= c;
464
 
 
465
 
    return(*this);
466
 
}
467
 
 
468
 
#ifndef SQUARE
469
 
#define SQUARE(x) ((x)*(x))
470
 
#endif
471
 
 
472
 
void matrix3x3::FillOrth(double Alpha,double Beta, double Gamma,
473
 
                         double A, double B, double C)
474
 
{
 
452
  }
 
453
 
 
454
  static inline double SQUARE( double x ) { return x*x; }
 
455
 
 
456
  void matrix3x3::FillOrth(double Alpha,double Beta, double Gamma,
 
457
                           double A, double B, double C)
 
458
  {
475
459
    double V;
476
460
 
477
461
    Alpha *= DEG_TO_RAD;
484
468
 
485
469
    // since we'll ultimately divide by (a * b), we've factored those out here
486
470
    V = C * sqrt(1 - SQUARE(cos(Alpha)) - SQUARE(cos(Beta)) - SQUARE(cos(Gamma))
487
 
                 + 2 * cos(Alpha) * cos(Beta) * cos(Gamma));
 
471
                 + 2 * cos(Alpha) * cos(Beta) * cos(Gamma));
488
472
 
489
473
    ele[0][0] = A;
490
474
    ele[0][1] = B*cos(Gamma);
497
481
    ele[2][0] = 0.0;
498
482
    ele[2][1] = 0.0;
499
483
    ele[2][2] = V / (sin(Gamma)); // again, we factored out A * B when defining V
500
 
}
501
 
 
502
 
ostream& operator<< ( ostream& co, const matrix3x3& m )
503
 
 
504
 
{
505
 
    co << "[ "
506
 
    << m.ele[0][0]
507
 
    << ", "
508
 
    << m.ele[0][1]
509
 
    << ", "
510
 
    << m.ele[0][2]
511
 
    << " ]" << endl;
512
 
 
513
 
    co << "[ "
514
 
    << m.ele[1][0]
515
 
    << ", "
516
 
    << m.ele[1][1]
517
 
    << ", "
518
 
    << m.ele[1][2]
519
 
    << " ]" << endl;
520
 
 
521
 
    co << "[ "
522
 
    << m.ele[2][0]
523
 
    << ", "
524
 
    << m.ele[2][1]
525
 
    << ", "
526
 
    << m.ele[2][2]
527
 
    << " ]" << endl;
 
484
  }
 
485
 
 
486
  /** Print a text representation of the matrix in the standardized form:
 
487
      [ a, b, c ] <br>
 
488
      [ d, e, f ] <br>
 
489
      [ g, h, i ] <br>
 
490
      where the letters represent the appropriate entries in the matrix.
 
491
      Uses the standard output format for the individual entries, separated
 
492
      by ", " for each column, and [ ] indicating each row.
 
493
   **/
 
494
  ostream& operator<< ( ostream& co, const matrix3x3& m )
 
495
 
 
496
  {
 
497
    co << "[ "
 
498
       << m.ele[0][0]
 
499
       << ", "
 
500
       << m.ele[0][1]
 
501
       << ", "
 
502
       << m.ele[0][2]
 
503
       << " ]" << endl;
 
504
 
 
505
    co << "[ "
 
506
       << m.ele[1][0]
 
507
       << ", "
 
508
       << m.ele[1][1]
 
509
       << ", "
 
510
       << m.ele[1][2]
 
511
       << " ]" << endl;
 
512
 
 
513
    co << "[ "
 
514
       << m.ele[2][0]
 
515
       << ", "
 
516
       << m.ele[2][1]
 
517
       << ", "
 
518
       << m.ele[2][2]
 
519
       << " ]" << endl;
528
520
 
529
521
    return co ;
530
 
}
 
522
  }
531
523
 
532
 
/*! This static function computes the eigenvalues and
533
 
  eigenvectors of a SYMMETRIC nxn matrix. This method is used
534
 
  internally by OpenBabel, but may be useful as a general
535
 
  eigenvalue finder.
536
 
  
537
 
  The algorithm uses Jacobi transformations. It is described
538
 
  e.g. in Wilkinson, Reinsch "Handbook for automatic computation,
539
 
  Volume II: Linear Algebra", part II, contribution II/1. The
540
 
  implementation is also similar to the implementation in this
541
 
  book. This method is adequate to solve the eigenproblem for
542
 
  small matrices, of size perhaps up to 10x10. For bigger
543
 
  problems, you might want to resort to the sophisticated routines
544
 
  of LAPACK.
545
 
  
546
 
  \note If you plan to find the eigenvalues of a symmetric 3x3
547
 
  matrix, you will probably prefer to use the more convenient
548
 
  method findEigenvectorsIfSymmetric()
549
 
  
550
 
  @param n the size of the matrix that should be diagonalized
551
 
  
552
 
  @param a array of size n^2 which holds the symmetric matrix
553
 
  whose eigenvectors are to be computed. The convention is that
554
 
  the entry in row r and column c is addressed as a[n*r+c] where,
555
 
  of course, 0 <= r < n and 0 <= c < n. There is no check that the
556
 
  matrix is actually symmetric. If it is not, the behaviour of
557
 
  this function is undefined.  On return, the matrix is
558
 
  overwritten with junk.
559
 
  
560
 
  @param d pointer to a field of at least n doubles which will be
561
 
  overwritten. On return of this function, the entries d[0]..d[n-1]
562
 
  will contain the eigenvalues of the matrix.
563
 
  
564
 
  @param v an array of size n^2 where the eigenvectors will be
565
 
  stored. On return, the columns of this matrix will contain the
566
 
  eigenvectors. The eigenvectors are normalized and mutually
567
 
  orthogonal.
568
 
*/
569
 
void matrix3x3::jacobi(unsigned int n, double *a, double *d, double *v)
570
 
{
 
524
  /*! This static function computes the eigenvalues and
 
525
    eigenvectors of a SYMMETRIC nxn matrix. This method is used
 
526
    internally by OpenBabel, but may be useful as a general
 
527
    eigenvalue finder.
 
528
  
 
529
    The algorithm uses Jacobi transformations. It is described
 
530
    e.g. in Wilkinson, Reinsch "Handbook for automatic computation,
 
531
    Volume II: Linear Algebra", part II, contribution II/1. The
 
532
    implementation is also similar to the implementation in this
 
533
    book. This method is adequate to solve the eigenproblem for
 
534
    small matrices, of size perhaps up to 10x10. For bigger
 
535
    problems, you might want to resort to the sophisticated routines
 
536
    of LAPACK.
 
537
  
 
538
    \note If you plan to find the eigenvalues of a symmetric 3x3
 
539
    matrix, you will probably prefer to use the more convenient
 
540
    method findEigenvectorsIfSymmetric()
 
541
  
 
542
    @param n the size of the matrix that should be diagonalized
 
543
  
 
544
    @param a array of size n^2 which holds the symmetric matrix
 
545
    whose eigenvectors are to be computed. The convention is that
 
546
    the entry in row r and column c is addressed as a[n*r+c] where,
 
547
    of course, 0 <= r < n and 0 <= c < n. There is no check that the
 
548
    matrix is actually symmetric. If it is not, the behaviour of
 
549
    this function is undefined.  On return, the matrix is
 
550
    overwritten with junk.
 
551
  
 
552
    @param d pointer to a field of at least n doubles which will be
 
553
    overwritten. On return of this function, the entries d[0]..d[n-1]
 
554
    will contain the eigenvalues of the matrix.
 
555
  
 
556
    @param v an array of size n^2 where the eigenvectors will be
 
557
    stored. On return, the columns of this matrix will contain the
 
558
    eigenvectors. The eigenvectors are normalized and mutually
 
559
    orthogonal.
 
560
  */
 
561
  void matrix3x3::jacobi(unsigned int n, double *a, double *d, double *v)
 
562
  {
571
563
    double onorm, dnorm;
572
564
    double b, dma, q, t, c, s;
573
565
    double  atemp, vtemp, dtemp;
579
571
    // Set v to the identity matrix, set the vector d to contain the
580
572
    // diagonal elements of the matrix a
581
573
    for (j = 0; j < static_cast<int>(n); j++)
582
 
    {
 
574
      {
583
575
        for (i = 0; i < static_cast<int>(n); i++)
584
 
            v[n*i+j] = 0.0;
 
576
          v[n*i+j] = 0.0;
585
577
        v[n*j+j] = 1.0;
586
578
        d[j] = a[n*j+j];
587
 
    }
 
579
      }
588
580
 
589
581
    nrot = MAX_SWEEPS;
590
582
    for (l = 1; l <= nrot; l++)
591
 
    {
 
583
      {
592
584
        // Set dnorm to be the maximum norm of the diagonal elements, set
593
585
        // onorm to the maximum norm of the off-diagonal elements
594
586
        dnorm = 0.0;
595
587
        onorm = 0.0;
596
588
        for (j = 0; j < static_cast<int>(n); j++)
597
 
        {
 
589
          {
598
590
            dnorm += (double)fabs(d[j]);
599
591
            for (i = 0; i < j; i++)
600
 
                onorm += (double)fabs(a[n*i+j]);
601
 
        }
 
592
              onorm += (double)fabs(a[n*i+j]);
 
593
          }
602
594
        // Normal end point of this algorithm.
603
595
        if((onorm/dnorm) <= 1.0e-12)
604
 
            goto Exit_now;
 
596
          goto Exit_now;
605
597
 
606
598
        for (j = 1; j < static_cast<int>(n); j++)
607
 
        {
 
599
          {
608
600
            for (i = 0; i <= j - 1; i++)
609
 
            {
 
601
              {
610
602
 
611
603
                b = a[n*i+j];
612
604
                if(fabs(b) > 0.0)
613
 
                {
 
605
                  {
614
606
                    dma = d[j] - d[i];
615
607
                    if((fabs(dma) + fabs(b)) <=  fabs(dma))
616
 
                        t = b / dma;
 
608
                      t = b / dma;
617
609
                    else
618
 
                    {
 
610
                      {
619
611
                        q = 0.5 * dma / b;
620
612
                        t = 1.0/((double)fabs(q) + (double)sqrt(1.0+q*q));
621
613
                        if (q < 0.0)
622
 
                            t = -t;
623
 
                    }
 
614
                          t = -t;
 
615
                      }
624
616
 
625
617
                    c = 1.0/(double)sqrt(t*t + 1.0);
626
618
                    s = t * c;
627
619
                    a[n*i+j] = 0.0;
628
620
 
629
621
                    for (k = 0; k <= i-1; k++)
630
 
                    {
 
622
                      {
631
623
                        atemp = c * a[n*k+i] - s * a[n*k+j];
632
624
                        a[n*k+j] = s * a[n*k+i] + c * a[n*k+j];
633
625
                        a[n*k+i] = atemp;
634
 
                    }
 
626
                      }
635
627
 
636
628
                    for (k = i+1; k <= j-1; k++)
637
 
                    {
 
629
                      {
638
630
                        atemp = c * a[n*i+k] - s * a[n*k+j];
639
631
                        a[n*k+j] = s * a[n*i+k] + c * a[n*k+j];
640
632
                        a[n*i+k] = atemp;
641
 
                    }
 
633
                      }
642
634
 
643
635
                    for (k = j+1; k < static_cast<int>(n); k++)
644
 
                    {
 
636
                      {
645
637
                        atemp = c * a[n*i+k] - s * a[n*j+k];
646
638
                        a[n*j+k] = s * a[n*i+k] + c * a[n*j+k];
647
639
                        a[n*i+k] = atemp;
648
 
                    }
 
640
                      }
649
641
 
650
642
                    for (k = 0; k < static_cast<int>(n); k++)
651
 
                    {
 
643
                      {
652
644
                        vtemp = c * v[n*k+i] - s * v[n*k+j];
653
645
                        v[n*k+j] = s * v[n*k+i] + c * v[n*k+j];
654
646
                        v[n*k+i] = vtemp;
655
 
                    }
 
647
                      }
656
648
 
657
649
                    dtemp = c*c*d[i] + s*s*d[j] - 2.0*c*s*b;
658
650
                    d[j] = s*s*d[i] + c*c*d[j] +  2.0*c*s*b;
659
651
                    d[i] = dtemp;
660
 
                } /* end if */
661
 
            } /* end for i */
662
 
        } /* end for j */
663
 
    } /* end for l */
 
652
                  } /* end if */
 
653
              } /* end for i */
 
654
          } /* end for j */
 
655
      } /* end for l */
664
656
 
665
 
Exit_now:
 
657
  Exit_now:
666
658
 
667
659
    // Now sort the eigenvalues (and the eigenvectors) so that the
668
660
    // smallest eigenvalues come first.
669
661
    nrot = l;
670
662
 
671
663
    for (j = 0; j < static_cast<int>(n)-1; j++)
672
 
    {
 
664
      {
673
665
        k = j;
674
666
        dtemp = d[k];
675
667
        for (i = j+1; i < static_cast<int>(n); i++)
676
 
            if(d[i] < dtemp)
 
668
          if(d[i] < dtemp)
677
669
            {
678
 
                k = i;
679
 
                dtemp = d[k];
 
670
              k = i;
 
671
              dtemp = d[k];
680
672
            }
681
673
 
682
674
        if(k > j)
683
 
        {
 
675
          {
684
676
            d[k] = d[j];
685
677
            d[j] = dtemp;
686
678
            for (i = 0; i < static_cast<int>(n); i++)
687
 
            {
 
679
              {
688
680
                dtemp = v[n*i+k];
689
681
                v[n*i+k] = v[n*i+j];
690
682
                v[n*i+j] = dtemp;
691
 
            }
692
 
        }
693
 
    }
694
 
}
 
683
              }
 
684
          }
 
685
      }
 
686
  }
695
687
 
696
688
} // end namespace OpenBabel
697
689