src/openbabel/vector3.cpp

/**********************************************************************
vector3.cpp - Handle 3D coordinates.
 
Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
 
This file is part of the Open Babel project.
For more information, see <http://openbabel.sourceforge.net/>
 
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.
***********************************************************************/

#include <math.h>

#include "mol.hpp"
#include "vector3.hpp"

using namespace std;

namespace OpenBabel
{

/*! \class vector3
   \brief Represents a vector in the 3-dimensional real space.

The vector3 class was designed to simplify operations with doubleing
point coordinates. To this end many of the common operations have been
overloaded for simplicity. Vector addition, subtraction, scalar
multiplication, dot product, cross product, magnitude and a number of
other utility functions are built in to the vector class. For a full
description of the class member functions please consult the header
file vector3.h. The following code demonstrates several of the
functions of the vector class:
\code
vector3 v1,v2,v3;
v1 = VX;
v2 = VY;
v3 = cross(v1,v2);
v3 *= 2.5;
v3.normalize();
\endcode
*/

/*! This (slow) method allows to access the elements of the
  vector as if it were an array of doubles. If the index is > 2,
  then a warning is printed, and the program is terminated via
  exit(-1). Otherwise, if i is 0, 1 or 2, then a reference to x,
  y or z is returned, respectively.
  
  \warning This method is primarily designed to facilitate the
  integration ('Open Babelization') of code that uses arrays of
  doubles rather than the vector class. Due to the error checks
  the method is of course very slow and should therefore be
  avoided in production code.
*/
double& vector3::operator[] ( unsigned int i)
{
    if (i > 2)
    {
        cerr << "ERROR in OpenBabel::vector3::operator[]" << endl
        << "The method has been called with an illegal index i=" << i << "." << endl
        << "Please contact the author of the offending program immediately." << endl;
        exit(-1);
    }
    if (i == 0)
        return _vx;
    if (i == 1)
        return _vy;
    return _vz;
}

/*! replaces *this with a random unit vector, which is (supposed
  to be) uniformly distributed over the unit sphere. Uses the
  random number generator obRand, or uses the system number
  generator with a time seed if obRand == NULL.
     
  @param obRandP random number generator to use, or 0L, if the
  system random number generator (with time seed) should be used
*/
void vector3::randomUnitVector(OBRandom *obRandP)
{
    OBRandom *ptr;
    if (!obRandP)
    {
        ptr = new OBRandom(true);
        ptr->TimeSeed();
    }
    else
        ptr = obRandP;

    // obtain a random vector with 0.001 <= length^2 <= 1.0, normalize
    // the vector to obtain a random vector of length 1.0.
    double l;
    do
    {
        this->Set(ptr->NextFloat()-0.5, ptr->NextFloat()-0.5, ptr->NextFloat()-0.5);
        l = length_2();
    }
    while ( (l > 1.0) || (l < 1e-4) );
    this->normalize();

    if (!obRandP)
        delete ptr;
}

OBAPI ostream& operator<< ( ostream& co, const vector3& v )
{
    co << "< " << v._vx << ", " << v._vy << ", " << v._vz << " >" ;
    return co ;
}

OBAPI int operator== ( const vector3& v1, const vector3& v2 )
{
    if ( ( v1._vx == v2._vx ) &&
            ( v1._vy == v2._vy ) &&
            ( v1._vz == v2._vz ) )
        return ( true ) ;
    else
        return ( false ) ;
}

OBAPI int operator!= ( const vector3& v1, const vector3& v2 )
{
    if ( ( v1._vx != v2._vx ) ||
            ( v1._vy != v2._vy ) ||
            ( v1._vz != v2._vz ) )
        return ( true ) ;
    else
        return ( false ) ;
}

/*! This method checks if the current vector has length() ==
  0.0.  If so, *this remains unchanged. Otherwise, *this is
  scaled by 1.0/length().

  \warning If length() is very close to zero, but not == 0.0,
  this method may behave in unexpected ways and return almost
  random results; details may depend on your particular doubleing
  point implementation. The use of this method is therefore
  highly discouraged, unless you are certain that length() is in
  a reasonable range, away from 0.0 (Stefan Kebekus)

  \deprecated This method will probably replaced by a safer
  algorithm in the future.

  \todo Replace this method with a more fool-proof version.

  @returns a reference to *this
*/
vector3& vector3 :: normalize ()
{
    double l = length ();

    if (IsNearZero(l))
        return(*this);

    _vx = _vx / l ;
    _vy = _vy / l ;
    _vz = _vz / l ;

    return(*this);
}

OBAPI double dot ( const vector3& v1, const vector3& v2 )
{
    return v1._vx*v2._vx + v1._vy*v2._vy + v1._vz*v2._vz ;
}

OBAPI vector3 cross ( const vector3& v1, const vector3& v2 )
{
    vector3 vv ;

    vv._vx =   v1._vy*v2._vz - v1._vz*v2._vy ;
    vv._vy = - v1._vx*v2._vz + v1._vz*v2._vx ;
    vv._vz =   v1._vx*v2._vy - v1._vy*v2._vx ;

    return ( vv ) ;
}


/*! This method calculates the angle between two vectors
     
  \warning If length() of any of the two vectors is == 0.0,
  this method will divide by zero. If the product of the
  length() of the two vectors is very close to 0.0, but not ==
  0.0, this method may behave in unexpected ways and return
  almost random results; details may depend on your particular
  doubleing point implementation. The use of this method is
  therefore highly discouraged, unless you are certain that the
  length()es are in a reasonable range, away from 0.0 (Stefan
  Kebekus)

  \deprecated This method will probably replaced by a safer
  algorithm in the future.

  \todo Replace this method with a more fool-proof version.

  @returns the angle in degrees (0-360)
*/
OBAPI double vectorAngle ( const vector3& v1, const vector3& v2 )
{
    double mag;
    double dp;

    mag = v1.length() * v2.length();
    dp = dot(v1,v2)/mag;

    if (dp < -0.999999)
        dp = -0.9999999;

    if (dp > 0.9999999)
        dp = 0.9999999;

    if (dp > 1.0)
        dp = 1.0;

    return((RAD_TO_DEG * acos(dp)));
}

OBAPI double CalcTorsionAngle(const vector3 &a, const vector3 &b,
                        const vector3 &c, const vector3 &d)
{
    double torsion;
    vector3 b1,b2,b3,c1,c2,c3;

    b1 = a - b;
    b2 = b - c;
    b3 = c - d;

    c1 = cross(b1,b2);
    c2 = cross(b2,b3);
    c3 = cross(c1,c2);

    if (c1.length() * c2.length() < 0.001)
        torsion = 0.0;
    else
    {
        torsion = vectorAngle(c1,c2);
        if (dot(b2,c3) > 0.0)
            torsion *= -1.0;
    }

    return(torsion);
}

/*! This method checks if the current vector *this is zero
  (i.e. if all entries == 0.0). If so, a warning message is
  printed, and the whole program is aborted with exit(0).
  Otherwise, a vector of length one is generated, which is
  orthogonal to *this, and stored in v. The resulting vector is
  not random.

  \warning If the entries of the *this (in particular the
  z-component) are very close to zero, but not == 0.0, this
  method may behave in unexpected ways and return almost random
  results; details may depend on your particular floating point
  implementation. The use of this method is therefore highly
  discouraged, unless you are certain that all components of
  *this are in a reasonable range, away from 0.0 (Stefan
  Kebekus)

  \deprecated This method will probably replaced by a safer
  algorithm in the future.

  \todo Replace this method with a more fool-proof version that
  does not call exit()

  @param res a reference to a vector where the result will be
  stored
*/
void vector3::createOrthoVector(vector3 &res) const
{
    vector3 cO;

    if ( ( IsNearZero(this->x())) && (IsNearZero(this->y())) )
    {
        if ( IsNearZero(this->z()) )
        {
            cerr << "makeorthovec zero vector" << endl;
            exit(0);
        }
        cO.SetX(1.0);
    }
    else
    {
        cO.SetZ(1.0);
    }
    res= cross(cO,*this);
    res.normalize();
}

const vector3 VZero ( 0.0, 0.0, 0.0 ) ;
const vector3 VX    ( 1.0, 0.0, 0.0 ) ;
const vector3 VY    ( 0.0, 1.0, 0.0 ) ;
const vector3 VZ    ( 0.0, 0.0, 1.0 ) ;

/* Calculate the distance of point a to the plane determined by b,c,d */
double Point2Plane(vector3 a, vector3 b, vector3 c, vector3 d)
{
  double angle =0;
  double dist_ab =0;
  vector3 v_ba = a-b;
  vector3 v_normal = cross(c-b, d-b).normalize();
  angle = vectorAngle(v_normal, v_ba);
  dist_ab = v_ba.length();
  return fabs(dist_ab * cos(DEG_TO_RAD * angle));
}

} // namespace OpenBabel

//! \file vector3.cpp
//! \brief Handle 3D coordinates.
Revision:	2440
Committed:	Wed Nov 16 19:42:11 2005 UTC (18 years, 8 months ago) by tim
File size:	9072 byte(s)
Log Message:	adding openbabel
#	User	Rev	Content
1	tim	2440	/**********************************************************************
2			vector3.cpp - Handle 3D coordinates.
3
4			Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
5			Some portions Copyright (C) 2001-2005 by Geoffrey R. Hutchison
6
7			This file is part of the Open Babel project.
8			For more information, see <http://openbabel.sourceforge.net/>
9
10			This program is free software; you can redistribute it and/or modify
11			it under the terms of the GNU General Public License as published by
12			the Free Software Foundation version 2 of the License.
13
14			This program is distributed in the hope that it will be useful,
15			but WITHOUT ANY WARRANTY; without even the implied warranty of
16			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17			GNU General Public License for more details.
18			***********************************************************************/
19
20			#include <math.h>
21
22			#include "mol.hpp"
23			#include "vector3.hpp"
24
25			using namespace std;
26
27			namespace OpenBabel
28			{
29
30			/*! \class vector3
31			\brief Represents a vector in the 3-dimensional real space.
32
33			The vector3 class was designed to simplify operations with doubleing
34			point coordinates. To this end many of the common operations have been
35			overloaded for simplicity. Vector addition, subtraction, scalar
36			multiplication, dot product, cross product, magnitude and a number of
37			other utility functions are built in to the vector class. For a full
38			description of the class member functions please consult the header
39			file vector3.h. The following code demonstrates several of the
40			functions of the vector class:
41			\code
42			vector3 v1,v2,v3;
43			v1 = VX;
44			v2 = VY;
45			v3 = cross(v1,v2);
46			v3 *= 2.5;
47			v3.normalize();
48			\endcode
49			*/
50
51			/*! This (slow) method allows to access the elements of the
52			vector as if it were an array of doubles. If the index is > 2,
53			then a warning is printed, and the program is terminated via
54			exit(-1). Otherwise, if i is 0, 1 or 2, then a reference to x,
55			y or z is returned, respectively.
56
57			\warning This method is primarily designed to facilitate the
58			integration ('Open Babelization') of code that uses arrays of
59			doubles rather than the vector class. Due to the error checks
60			the method is of course very slow and should therefore be
61			avoided in production code.
62			*/
63			double& vector3::operator[] ( unsigned int i)
64			{
65			if (i > 2)
66			{
67			cerr << "ERROR in OpenBabel::vector3::operator[]" << endl
68			<< "The method has been called with an illegal index i=" << i << "." << endl
69			<< "Please contact the author of the offending program immediately." << endl;
70			exit(-1);
71			}
72			if (i == 0)
73			return _vx;
74			if (i == 1)
75			return _vy;
76			return _vz;
77			}
78
79			/! replaces this with a random unit vector, which is (supposed
80			to be) uniformly distributed over the unit sphere. Uses the
81			random number generator obRand, or uses the system number
82			generator with a time seed if obRand == NULL.
83
84			@param obRandP random number generator to use, or 0L, if the
85			system random number generator (with time seed) should be used
86			*/
87			void vector3::randomUnitVector(OBRandom *obRandP)
88			{
89			OBRandom *ptr;
90			if (!obRandP)
91			{
92			ptr = new OBRandom(true);
93			ptr->TimeSeed();
94			}
95			else
96			ptr = obRandP;
97
98			// obtain a random vector with 0.001 <= length^2 <= 1.0, normalize
99			// the vector to obtain a random vector of length 1.0.
100			double l;
101			do
102			{
103			this->Set(ptr->NextFloat()-0.5, ptr->NextFloat()-0.5, ptr->NextFloat()-0.5);
104			l = length_2();
105			}
106			while ( (l > 1.0) \|\| (l < 1e-4) );
107			this->normalize();
108
109			if (!obRandP)
110			delete ptr;
111			}
112
113			OBAPI ostream& operator<< ( ostream& co, const vector3& v )
114			{
115			co << "< " << v._vx << ", " << v._vy << ", " << v._vz << " >" ;
116			return co ;
117			}
118
119			OBAPI int operator== ( const vector3& v1, const vector3& v2 )
120			{
121			if ( ( v1._vx == v2._vx ) &&
122			( v1._vy == v2._vy ) &&
123			( v1._vz == v2._vz ) )
124			return ( true ) ;
125			else
126			return ( false ) ;
127			}
128
129			OBAPI int operator!= ( const vector3& v1, const vector3& v2 )
130			{
131			if ( ( v1._vx != v2._vx ) \|\|
132			( v1._vy != v2._vy ) \|\|
133			( v1._vz != v2._vz ) )
134			return ( true ) ;
135			else
136			return ( false ) ;
137			}
138
139			/*! This method checks if the current vector has length() ==
140			0.0. If so, this remains unchanged. Otherwise, this is
141			scaled by 1.0/length().
142
143			\warning If length() is very close to zero, but not == 0.0,
144			this method may behave in unexpected ways and return almost
145			random results; details may depend on your particular doubleing
146			point implementation. The use of this method is therefore
147			highly discouraged, unless you are certain that length() is in
148			a reasonable range, away from 0.0 (Stefan Kebekus)
149
150			\deprecated This method will probably replaced by a safer
151			algorithm in the future.
152
153			\todo Replace this method with a more fool-proof version.
154
155			@returns a reference to *this
156			*/
157			vector3& vector3 :: normalize ()
158			{
159			double l = length ();
160
161			if (IsNearZero(l))
162			return(*this);
163
164			_vx = _vx / l ;
165			_vy = _vy / l ;
166			_vz = _vz / l ;
167
168			return(*this);
169			}
170
171			OBAPI double dot ( const vector3& v1, const vector3& v2 )
172			{
173			return v1._vxv2._vx + v1._vyv2._vy + v1._vz*v2._vz ;
174			}
175
176			OBAPI vector3 cross ( const vector3& v1, const vector3& v2 )
177			{
178			vector3 vv ;
179
180			vv._vx = v1._vyv2._vz - v1._vzv2._vy ;
181			vv._vy = - v1._vxv2._vz + v1._vzv2._vx ;
182			vv._vz = v1._vxv2._vy - v1._vyv2._vx ;
183
184			return ( vv ) ;
185			}
186
187
188			/*! This method calculates the angle between two vectors
189
190			\warning If length() of any of the two vectors is == 0.0,
191			this method will divide by zero. If the product of the
192			length() of the two vectors is very close to 0.0, but not ==
193			0.0, this method may behave in unexpected ways and return
194			almost random results; details may depend on your particular
195			doubleing point implementation. The use of this method is
196			therefore highly discouraged, unless you are certain that the
197			length()es are in a reasonable range, away from 0.0 (Stefan
198			Kebekus)
199
200			\deprecated This method will probably replaced by a safer
201			algorithm in the future.
202
203			\todo Replace this method with a more fool-proof version.
204
205			@returns the angle in degrees (0-360)
206			*/
207			OBAPI double vectorAngle ( const vector3& v1, const vector3& v2 )
208			{
209			double mag;
210			double dp;
211
212			mag = v1.length() * v2.length();
213			dp = dot(v1,v2)/mag;
214
215			if (dp < -0.999999)
216			dp = -0.9999999;
217
218			if (dp > 0.9999999)
219			dp = 0.9999999;
220
221			if (dp > 1.0)
222			dp = 1.0;
223
224			return((RAD_TO_DEG * acos(dp)));
225			}
226
227			OBAPI double CalcTorsionAngle(const vector3 &a, const vector3 &b,
228			const vector3 &c, const vector3 &d)
229			{
230			double torsion;
231			vector3 b1,b2,b3,c1,c2,c3;
232
233			b1 = a - b;
234			b2 = b - c;
235			b3 = c - d;
236
237			c1 = cross(b1,b2);
238			c2 = cross(b2,b3);
239			c3 = cross(c1,c2);
240
241			if (c1.length() * c2.length() < 0.001)
242			torsion = 0.0;
243			else
244			{
245			torsion = vectorAngle(c1,c2);
246			if (dot(b2,c3) > 0.0)
247			torsion *= -1.0;
248			}
249
250			return(torsion);
251			}
252
253			/! This method checks if the current vector this is zero
254			(i.e. if all entries == 0.0). If so, a warning message is
255			printed, and the whole program is aborted with exit(0).
256			Otherwise, a vector of length one is generated, which is
257			orthogonal to *this, and stored in v. The resulting vector is
258			not random.
259
260			\warning If the entries of the *this (in particular the
261			z-component) are very close to zero, but not == 0.0, this
262			method may behave in unexpected ways and return almost random
263			results; details may depend on your particular floating point
264			implementation. The use of this method is therefore highly
265			discouraged, unless you are certain that all components of
266			*this are in a reasonable range, away from 0.0 (Stefan
267			Kebekus)
268
269			\deprecated This method will probably replaced by a safer
270			algorithm in the future.
271
272			\todo Replace this method with a more fool-proof version that
273			does not call exit()
274
275			@param res a reference to a vector where the result will be
276			stored
277			*/
278			void vector3::createOrthoVector(vector3 &res) const
279			{
280			vector3 cO;
281
282			if ( ( IsNearZero(this->x())) && (IsNearZero(this->y())) )
283			{
284			if ( IsNearZero(this->z()) )
285			{
286			cerr << "makeorthovec zero vector" << endl;
287			exit(0);
288			}
289			cO.SetX(1.0);
290			}
291			else
292			{
293			cO.SetZ(1.0);
294			}
295			res= cross(cO,*this);
296			res.normalize();
297			}
298
299			const vector3 VZero ( 0.0, 0.0, 0.0 ) ;
300			const vector3 VX ( 1.0, 0.0, 0.0 ) ;
301			const vector3 VY ( 0.0, 1.0, 0.0 ) ;
302			const vector3 VZ ( 0.0, 0.0, 1.0 ) ;
303
304			/* Calculate the distance of point a to the plane determined by b,c,d */
305			double Point2Plane(vector3 a, vector3 b, vector3 c, vector3 d)
306			{
307			double angle =0;
308			double dist_ab =0;
309			vector3 v_ba = a-b;
310			vector3 v_normal = cross(c-b, d-b).normalize();
311			angle = vectorAngle(v_normal, v_ba);
312			dist_ab = v_ba.length();
313			return fabs(dist_ab * cos(DEG_TO_RAD * angle));
314			}
315
316			} // namespace OpenBabel
317
318			//! \file vector3.cpp
319			//! \brief Handle 3D coordinates.