vector3d.h

Go to the documentation of this file.

/***************************************************************************
                          vector3d.h  -  description
                             -------------------
    begin                : Mon Mar 3 2003
    copyright            : (C) 2003-2006 by Ben Swerts
    email                : bswerts@users.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; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/





#ifndef VECTOR3D_H
#define VECTOR3D_H


// C++ header files
#include <cmath>
#include <iostream>

// Xbrabo header files
#include <point3d.h>
#include <quaternion.h>

template<class T> class Vector3D
{
  public:
    Vector3D();                         // default constructor
    Vector3D(const T x, const T y, const T z);    // constructs a vector (xyz)
    Vector3D(const T x1, const T y1, const T z1, const T x2, const T y2, const T z2);     // constructs a vector from point1 to point2
    Vector3D(const Point3D<T> point1, const Point3D<T> point2);       // constructs a vector from point1 to point2
    Vector3D(const Vector3D& v);        // copy constructor
    ~Vector3D();                        // destructor

    // public member functions for retrieving data
    bool isUnit() const;                // returns true of it's a unit vector
    bool isZero() const;                // returns true if the vector is zero
    T x() const;                        // returns the x-value
    T y() const;                        // returns the y-value
    T z() const;                        // returns the z-value
    T length() const;                   // returns the length of the vector

    // public member functions for changing data
    void setValues(const T x, const T y, const T z);        // sets the vector
    void normalize();                   // normalizes the vector
    void setLength(const T newLength);  // sets a new length for the vector
    void changeLength(const T amount);  // changes the length of the vector by amount
    void invert();                      // inverts the direction of the vector
    void add(const Vector3D v);         // adds the vector v
    void rotate(const Vector3D axis, const T angle);        // rotate the vector around the axis by an angle
    void setTorsion(const T angle, const Vector3D refBond, const Vector3D centralBond);   // sets the torsion angle

    // public member functions for data generation
    Vector3D cross(const Vector3D v) const;       // calculates the cross product of 2 vectors
    T dot(const Vector3D v) const;      // calculates the dot product of 2 vectors
    T angle(const Vector3D v) const;    // calculates the angle between 2 vectors
    T torsion(const Vector3D v1, const Vector3D v2) const;  // calculates the torsion angle with another vector

  private:
    // private member variables
    T xVect, yVect, zVect;              
};


template <class T> Vector3D<T>::Vector3D()
{
  setValues(0.0, 0.0, 0.0);
}

template <class T> Vector3D<T>::Vector3D(const T x, const T y, const T z)
{
  setValues(x, y, z);
}

template <class T> Vector3D<T>::Vector3D(const T x1, const T y1, const T z1, const T x2, const T y2, const T z2)
{
  setValues(x2 - x1, y2 - y1, z2 - z1);
}

template <class T> Vector3D<T>::Vector3D(const Point3D<T> point1, const Point3D<T> point2)
{
  setValues(point2.x() - point1.x(), point2.y() - point1.y(), point2.z() - point1.z());
}

template <class T> Vector3D<T>::Vector3D(const Vector3D& v)
{
  xVect = v.xVect;
  yVect = v.yVect;
  zVect = v.zVect;
}

template <class T> Vector3D<T>::~Vector3D()
{

}

template <class T> bool Vector3D<T>::isUnit() const
{
  if(fabs(length() - 1.0) < Point3D<T>::TOLERANCE)
    return true;
  else
    return false;
}

template <class T> bool Vector3D<T>::isZero() const
{
  if(fabs(length()) < Point3D<T>::TOLERANCE)
    return true;
  else
    return false;
}

template <class T> T Vector3D<T>::x() const
{
  return xVect;
}

template <class T> T Vector3D<T>::y() const
{
  return yVect;
}
template <class T> T Vector3D<T>::z() const
{
  return zVect;
}

template <class T> T Vector3D<T>::length() const
{
  return sqrt(xVect*xVect + yVect*yVect + zVect*zVect);
}

template <class T> void Vector3D<T>::setValues(const T x, const T y, const T z)
{
  xVect = x;
  yVect = y;
  zVect = z;
}

template <class T> void Vector3D<T>::normalize()
{
  T len = length();
  if(len < Point3D<T>::TOLERANCE)
  {
    xVect = 0.0;
    yVect = 0.0;
    zVect = 0.0;
  }
  else
  {
    xVect = xVect/len;
    yVect = yVect/len;
    zVect = zVect/len;
  }
}

template <class T> void Vector3D<T>::setLength(const T newLength)
{
  // determine the scaling factor
  T scalingFactor = newLength/length();
  // scale the values
  xVect *= scalingFactor;
  yVect *= scalingFactor;
  zVect *= scalingFactor;

}

template <class T> void Vector3D<T>::changeLength(const T amount)
{
  // determine the scaling factor
  T oldLength = length();
  if(oldLength < Point3D<T>::TOLERANCE)
    return;
  T scalingFactor = (oldLength + amount)/oldLength;
  // scale the values
  xVect *= scalingFactor;
  yVect *= scalingFactor;
  zVect *= scalingFactor;

}

template <class T> void Vector3D<T>::invert()
{
  xVect = -xVect;
  yVect = -yVect;
  zVect = -zVect;
}

template <class T> void Vector3D<T>::add(const Vector3D v)
{
  xVect += v.xVect;
  yVect += v.yVect;
  zVect += v.zVect;
}

template <class T> void Vector3D<T>::rotate(const Vector3D axis, const T angle)
{
  Quaternion<T> rotation(axis, angle);
  Quaternion<T> inverse = rotation;
  inverse.inverse();
  Quaternion<T> thisVector(0.0, xVect, yVect, zVect);
  Quaternion<T> result = rotation*thisVector*inverse;
  result.normalize();
  setValues(result.x(), result.y(), result.z());
}

template <class T> Vector3D<T> Vector3D<T>::cross(const Vector3D v) const
{
  T resultX = yVect * v.zVect - v.yVect * zVect;
  T resultY = zVect * v.xVect - v.zVect * xVect;
  T resultZ = xVect * v.yVect - v.xVect * yVect;

  return Vector3D(resultX, resultY, resultZ);
}

template <class T> T Vector3D<T>::dot(const Vector3D v) const
{
  return xVect*v.xVect + yVect * v.yVect + zVect * v.zVect;
}

template <class T> T Vector3D<T>::angle(const Vector3D v) const
{
  Vector3D<T> v1 = v;
  Vector3D<T> v2 = *this;
  v1.normalize();
  v2.normalize();
  return acos(v1.dot(v2))*Point3D<T>::RADTODEG;
}

template <class T> T Vector3D<T>::torsion(const Vector3D v1, const Vector3D v2) const
{
  Vector3D<T> va = *this; // the first vector
  Vector3D<T> vb = v2; // the central vector
  Vector3D<T> vc = v1; // the second vector

  va.normalize();
  vb.normalize();
  vc.normalize();

  T argx = - va.dot(vc) + va.dot(vb) * vb.dot(vc);
  T argy = va.dot(vb.cross(vc));

  T acosAngle = acos(argx/sqrt(argx*argx + argy*argy))*Point3D<T>::RADTODEG;
//#ifdef Q_OS_WIN32
//  T asinAngle = asin(argy/sqrt(hypot(argx, argy)))*RADTODEG; //=> doesn't work too well when using templates
//#else
  T asinAngle = asin(argy/sqrt(argx*argx + argy*argy))*Point3D<T>::RADTODEG;
//#endif

#ifdef Q_OS_WIN32
  if(acosAngle == 90.0)
    return 90.0; // because the 'hypot' trick doesn't work when changing to template code
  else if(acosAngle > 90.0)
#else      
  if(acosAngle >= 90.0) // 1st and 4th
#endif
    return -asinAngle;
  else if(asinAngle < 0.0) // 3rd
    return asinAngle + 180.0;
  else
    return asinAngle - 180.0; // 2nd

  return 0.0; // omits warnings
}

template <class T> void Vector3D<T>::setTorsion(const T angle, const Vector3D refBond, const Vector3D centralBond)
{
  // new implementation: - get the current torsion angle
  //                     - calculate the difference with the requested torsion angle
  //                     - call rotate() to do the rotation around the central bond

  /*
  T torsionAngle = angle * DEGTORAD;
  Vector3D va = *this;
  Vector3D vb = centralBond;
  Vector3D vc = centralBond;
  vc.invert();
  Vector3D vd = refBond;

  T orgLength = va.length();
  va.normalize();
  vb.normalize();
  vc.normalize();
  vd.normalize();

  Vector3D v1 = va.cross(vb);
  Vector3D v2 = vc.cross(vd);

  while(torsionAngle > 180.0f)
    torsionAngle -= 360.0f;
  while(torsionAngle < -180.0f)
    torsionAngle += 360.0f;

  if(torsionAngle <= 90.0f && torsionAngle >= -90.0f) // 1st or 4th quadrant
  {
    v1.rotate(vb, torsionAngle);
  }
  */

  std::cerr << "Vector3D::setTorsion not implemented yet. DO NOT USE" << std::endl;
  exit(1);
}

#endif

---