git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@641 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
466
tools/pymol_asphere/src/cartesian.h
Executable file
466
tools/pymol_asphere/src/cartesian.h
Executable file
@ -0,0 +1,466 @@
|
||||
/****************************************************************************
|
||||
* cartesian.h
|
||||
* Shared stuff for dealing with Cartesian coordinates
|
||||
* Also contains quaternion structures and operations
|
||||
*
|
||||
* Cartesian point of doubles: cPt
|
||||
* Cartesian point of intergers: iPt
|
||||
* Vector of doubles: vectorPt
|
||||
* Color of doubles: colorPt
|
||||
* Quaternion: Quaternion
|
||||
*
|
||||
* Can be accessed as cPt.x, cPt[X], colorPt[GREEN], iPt[I], etc.
|
||||
*
|
||||
*
|
||||
* W. Michael Brown
|
||||
****************************************************************************/
|
||||
|
||||
/*! \file */
|
||||
|
||||
#ifndef CARTESIAN_H
|
||||
#define CARTESIAN_H
|
||||
|
||||
#include "miscm.h"
|
||||
#include "m_constants.h"
|
||||
#include "spherical.h"
|
||||
#include <assert.h>
|
||||
#include <math.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
using namespace std;
|
||||
|
||||
enum { X, ///<0
|
||||
Y, ///<1
|
||||
Z ///<2
|
||||
};
|
||||
enum { I, ///<0
|
||||
J, ///<1
|
||||
K, ///<2
|
||||
W ///<3
|
||||
};
|
||||
enum { RED, ///<0
|
||||
GREEN, ///<1
|
||||
BLUE ///<2
|
||||
};
|
||||
|
||||
// Other coordinates
|
||||
template<class numtyp> class Ball;
|
||||
|
||||
// Template friend declarations
|
||||
template<class numtyp> class TwoD;
|
||||
template<class numtyp>
|
||||
ostream & operator<< (ostream &out, const TwoD<numtyp> &t);
|
||||
template<class numtyp>
|
||||
istream & operator>> (istream &in, TwoD<numtyp> &t);
|
||||
|
||||
template<class numtyp> class ThreeD;
|
||||
template<class numtyp>
|
||||
ostream & operator<< (ostream &out, const ThreeD<numtyp> &t);
|
||||
template<class numtyp>
|
||||
istream & operator>> (istream &in, ThreeD<numtyp> &t);
|
||||
template<class numtyp>
|
||||
ThreeD<numtyp> operator+ (const numtyp, const ThreeD<numtyp> &two);
|
||||
template<class numtyp>
|
||||
ThreeD<numtyp> operator- (const numtyp, const ThreeD<numtyp> &two);
|
||||
template<class numtyp>
|
||||
ThreeD<numtyp> operator* (const numtyp, const ThreeD<numtyp> &two);
|
||||
template<class numtyp>
|
||||
ThreeD<numtyp> operator/ (const numtyp, const ThreeD<numtyp> &two);
|
||||
|
||||
/// Two dimensional vector
|
||||
/** The elements can be accessed directly .x or .y
|
||||
* or by using the operator [] ( [X], [Y] or [I], [J] )
|
||||
*
|
||||
* The following operators are currently overloaded:
|
||||
* +,-,*
|
||||
*
|
||||
* operators *,/ returns a vector with each element multiplied
|
||||
* times the corresponding element in the other vector (Matlab .* )
|
||||
*
|
||||
* the member function dot can be used to compute dot products
|
||||
*
|
||||
* Input and output are overloaded for element I/O of the form "x y"
|
||||
* <<, >>
|
||||
*
|
||||
* \sa cartesian.h for typedefs and defines*/
|
||||
template<class numtyp>
|
||||
class TwoD {
|
||||
public:
|
||||
/// Empty construct. Not necessarily initialized to [0 0]
|
||||
TwoD();
|
||||
/// Assign both x and y the value
|
||||
TwoD(numtyp x);
|
||||
/// Assignment Constructor
|
||||
TwoD(numtyp cx, numtyp cy);
|
||||
|
||||
/// Type conversion
|
||||
TwoD(const TwoD<float> &two);
|
||||
/// Ball projection (onto xy-plane)
|
||||
TwoD(const Ball<float> &ball);
|
||||
|
||||
numtyp x; ///< First element
|
||||
numtyp y; ///< Last element
|
||||
|
||||
numtyp &operator[](unsigned i);
|
||||
|
||||
friend ostream & operator<< <>(ostream &out, const TwoD<numtyp> &t);
|
||||
friend istream & operator>> <>(istream &in, TwoD<numtyp> &t);
|
||||
|
||||
TwoD<numtyp> operator + (const TwoD<numtyp> &two) const;
|
||||
void operator += (const TwoD<numtyp> &two);
|
||||
TwoD<numtyp> operator + (const numtyp two) const;
|
||||
TwoD<numtyp> operator - (const TwoD<numtyp> &two) const;
|
||||
TwoD<numtyp> operator - (const numtyp two) const;
|
||||
TwoD<numtyp> operator * (const numtyp two) const;
|
||||
TwoD<numtyp> operator * (const TwoD<numtyp> &two) const;
|
||||
void operator /= (const numtyp two);
|
||||
|
||||
bool operator != (const TwoD<numtyp> &two) const;
|
||||
|
||||
/// Dot Product
|
||||
numtyp dot(const TwoD<numtyp> &two) const;
|
||||
|
||||
/// Distance between two points
|
||||
numtyp dist(const TwoD<numtyp> &two) const;
|
||||
/// Distance squared between two points
|
||||
numtyp dist2(const TwoD<numtyp> &two) const;
|
||||
/// Returns one of two normals to a line represented by vector
|
||||
TwoD<numtyp> normal();
|
||||
|
||||
/// Move coordinates into array
|
||||
void to_array(numtyp *array);
|
||||
/// Set coordinates from array
|
||||
void from_array(numtyp *array);
|
||||
|
||||
// -------------- Weird functions that help with coord templating
|
||||
unsigned dimensionality();
|
||||
/// Returns the index of *this (for unsigned) in a square 3D array
|
||||
/** \param s_size s_size[0]=1Dsize **/
|
||||
numtyp array_index(vector<unsigned> &s_size);
|
||||
/// Increment a 2D index from min to max (same as nested for)
|
||||
/** Returns false when increment is complete **/
|
||||
bool increment_index(TwoD &minp,TwoD &maxp);
|
||||
/// Return false if x or y is not within the inclusive range
|
||||
bool check_bounds(numtyp min,numtyp max);
|
||||
private:
|
||||
};
|
||||
|
||||
///\var typedef TwoD<double> c2DPt
|
||||
/// Two dimensional vector of doubles
|
||||
typedef TwoD<double> c2DPt;
|
||||
|
||||
|
||||
/// Three dimensional vector
|
||||
/** The elements can be accessed directly .x or .y or .z
|
||||
* or by using the operator [] ( [X], [Y], [Z] or [I], [J], [K] or
|
||||
* [RED], [GREEN], [BLUE] )
|
||||
*
|
||||
* The following operators are currently overloaded:
|
||||
* +,-,*,/,+=,-=,*=,/=,==,!=
|
||||
*
|
||||
* operators *,/,*=,/= returns a vector with each element multiplied
|
||||
* times the corresponding element in the other vector (Matlab .* )
|
||||
* or with each element multiplied by a scalar
|
||||
*
|
||||
* the member function dot can be used to compute dot products
|
||||
|
||||
*
|
||||
* Input and output are overloaded for element I/O of the form "x y z"
|
||||
* <<, >>
|
||||
*
|
||||
* \sa cartesian.h for typedefs and defines*/
|
||||
template<class numtyp>
|
||||
class ThreeD {
|
||||
public:
|
||||
/// Assignment Constructor
|
||||
ThreeD(numtyp cx, numtyp cy, numtyp cz);
|
||||
/// Assign all the value
|
||||
ThreeD(numtyp in);
|
||||
/// Type Conversion
|
||||
ThreeD(const ThreeD<unsigned>&);
|
||||
/// Type Conversion
|
||||
ThreeD(const ThreeD<int>&);
|
||||
/// Type Conversion
|
||||
ThreeD(const ThreeD<float>&);
|
||||
/// Type Conversion
|
||||
ThreeD(const ThreeD<double>&);
|
||||
/// TwoD with (z-coordinate set to zero)
|
||||
ThreeD(const TwoD<float>&);
|
||||
/// Spherical Conversion
|
||||
ThreeD(const Ball<double>&);
|
||||
/// Spherical Conversion
|
||||
ThreeD(const Ball<float>&);
|
||||
/// Empty construct (Not necessarily initialized to zero)
|
||||
ThreeD();
|
||||
|
||||
numtyp x; ///< First Element
|
||||
numtyp y; ///< Second Element
|
||||
numtyp z; ///< Last Element
|
||||
|
||||
friend ostream & operator<< <>(ostream &out, const ThreeD<numtyp> &t);
|
||||
friend istream & operator>> <>(istream &in, ThreeD<numtyp> &t);
|
||||
friend ThreeD<numtyp> operator+ <>(const numtyp, const ThreeD<numtyp> &two);
|
||||
friend ThreeD<numtyp> operator- <>(const numtyp, const ThreeD<numtyp> &two);
|
||||
friend ThreeD<numtyp> operator* <>(const numtyp, const ThreeD<numtyp> &two);
|
||||
friend ThreeD<numtyp> operator/ <>(const numtyp, const ThreeD<numtyp> &two);
|
||||
|
||||
numtyp &operator[](unsigned i);
|
||||
numtyp operator[](unsigned i) const;
|
||||
|
||||
bool operator == (const ThreeD<numtyp> &two) const;
|
||||
bool operator != (const ThreeD<numtyp> &two) const;
|
||||
ThreeD<numtyp> operator + (const ThreeD<numtyp> &two) const;
|
||||
ThreeD<numtyp> operator + (const numtyp &two) const;
|
||||
ThreeD<numtyp> operator - (const ThreeD<numtyp> &two) const;
|
||||
ThreeD<numtyp> operator - (const numtyp &two) const;
|
||||
ThreeD<numtyp> operator * (const numtyp &two) const;
|
||||
ThreeD<numtyp> operator * (const ThreeD<numtyp> &two) const;
|
||||
ThreeD<numtyp> operator / (const numtyp &two) const;
|
||||
ThreeD<numtyp> operator / (const ThreeD<numtyp> &two) const;
|
||||
void operator = (const ThreeD &two);
|
||||
void operator += (const numtyp &two);
|
||||
void operator += (const ThreeD &two);
|
||||
void operator -= (const numtyp &two);
|
||||
void operator -= (const ThreeD &two);
|
||||
void operator *= (const numtyp &two);
|
||||
void operator /= (const numtyp &two);
|
||||
|
||||
/// Move coordinates into array
|
||||
void to_array(numtyp *array);
|
||||
/// Set coordinates from array
|
||||
void from_array(numtyp *array);
|
||||
|
||||
/// Returns the dot product of *this and two
|
||||
numtyp dot(const ThreeD<numtyp> &two) const;
|
||||
|
||||
/// Returns the cross product of \b *this and \b two
|
||||
/** The input vectors do not need to be normalized, however, the output
|
||||
* vector will not be */
|
||||
ThreeD<numtyp> cross(const ThreeD<numtyp> &two) const;
|
||||
/// Returns the angle between \b *this and \b two
|
||||
numtyp angle(const ThreeD<numtyp> &two) const;
|
||||
/// Returns an arbitrary vector that is perpendicular to the input
|
||||
/** The output vector is not normalized */
|
||||
ThreeD<numtyp> perpendicular();
|
||||
/// Rotate a vector in xz-plane by t radians
|
||||
ThreeD<numtyp> rotatey(double t) const;
|
||||
/// Rotate a vector in xy-plane by t radians
|
||||
ThreeD<numtyp> rotatez(double t) const;
|
||||
/// Magnitude of vector
|
||||
numtyp norm() const;
|
||||
/// Squared norm of a vector
|
||||
numtyp norm2() const;
|
||||
/// Magnitude of vector
|
||||
numtyp hypot() const;
|
||||
/// Distance between two points
|
||||
numtyp dist(const ThreeD<numtyp> &two);
|
||||
/// Distance squared between two points
|
||||
numtyp dist2(const ThreeD<numtyp> &two);
|
||||
/// Converts \b *this to the unit vector
|
||||
void normalize();
|
||||
/// Return the unit vector of \b *this
|
||||
ThreeD<numtyp> unit();
|
||||
|
||||
/// Returns the projection of the point or vector onto Z-plane
|
||||
c2DPt projz();
|
||||
|
||||
/// Set this to be the max of one and two for each dimension
|
||||
void max(ThreeD &one, ThreeD &two);
|
||||
/// Set this to be the min of one and two for each dimension
|
||||
void min(ThreeD &one, ThreeD &two);
|
||||
|
||||
// -------------- Weird functions that help with coord templating
|
||||
/// Returns 3
|
||||
unsigned dimensionality();
|
||||
/// Returns the index of *this (for unsigned) in a square 3D array
|
||||
/** \param s_size s_size[0]=1Dsize, and s_size[1]=1D size*1Dsize **/
|
||||
numtyp array_index(vector<unsigned> &s_size);
|
||||
/// Increment a 3D index from min to max (same as nested for)
|
||||
/** \note This is currently only implemented for unsigned numbers
|
||||
* Returns false when increment is complete **/
|
||||
bool increment_index(ThreeD &minp,ThreeD &maxp);
|
||||
/// Return false if x,y, or z is not within the inclusive range
|
||||
bool check_bounds(numtyp min,numtyp max);
|
||||
private:
|
||||
};
|
||||
|
||||
///\var typedef ThreeD<int> iPt;
|
||||
/// Three dimensional vector of integers
|
||||
typedef ThreeD<int> iPt;
|
||||
///\var typedef ThreeD<unsigned> uPt;
|
||||
/// Three dimensional vector of unsigned
|
||||
typedef ThreeD<unsigned> uPt;
|
||||
///\var typedef ThreeD<double> cPt;
|
||||
/// Three dimensional vector of doubles
|
||||
typedef ThreeD<double> cPt;
|
||||
///\var typedef ThreeD<double> vectorPt;
|
||||
/// Three dimensional vector of doubles
|
||||
typedef ThreeD<double> vectorPt;
|
||||
///\var typedef ThreeD<double> colorPt;
|
||||
/// Three dimensional vector of doubles
|
||||
typedef ThreeD<double> colorPt;
|
||||
|
||||
///\def ORIGIN
|
||||
/// Point at origin
|
||||
#define ORIGIN cPt(0.0,0.0,0.0)
|
||||
///\def XAXIS
|
||||
/// Unit vector for x-axis
|
||||
#define XAXIS vectorPt(1.0,0.0,0.0)
|
||||
///\def YAXIS
|
||||
/// Unit vector for y-axis
|
||||
#define YAXIS vectorPt(0.0,1.0,0.0)
|
||||
///\def ZAXIS
|
||||
/// Unit vector for z-axis
|
||||
#define ZAXIS vectorPt(0.0,0.0,1.0)
|
||||
|
||||
class RotMat;
|
||||
|
||||
/// Euler Rotation
|
||||
class EulerRot {
|
||||
public:
|
||||
EulerRot();
|
||||
EulerRot(double theta,double psi,double phi);
|
||||
|
||||
/// Rotate the rotation axis by a given rotation matrix
|
||||
void rotate_axis(const RotMat &rotmat);
|
||||
|
||||
void operator= (const EulerRot &two);
|
||||
friend ostream & operator<<(ostream &out, const EulerRot &rot);
|
||||
|
||||
double theta;
|
||||
double psi;
|
||||
double phi;
|
||||
};
|
||||
|
||||
///\def EU_NOROTATE
|
||||
/// EulerRot with no rotation
|
||||
#define EU_NOROTATE EulerRot(0.0,0.0,0.0)
|
||||
|
||||
//---------------------------Quaternion Stuff -------------------------
|
||||
|
||||
/// Class for handling Quaternion vectors
|
||||
/** The elements can be accessed directly .w .i .j or .k
|
||||
* or by using the operator[] ( [W], [X], [Y], [Z] or [W], [I], [J], [K] )
|
||||
*
|
||||
* The following operators are currently overloaded:
|
||||
* +=,*,*=,=,==
|
||||
*
|
||||
* Overloaded * concatenates two successive rotations \n
|
||||
* \e It \e is \e not \e communative \n
|
||||
* For \e join=q1*q2, \e join is equivalent to performing rotation \e q1
|
||||
* \b followed by \e q2.
|
||||
*
|
||||
* Input and output are overloaded for element I/O of the form "w i j k"
|
||||
* <<, >>
|
||||
*
|
||||
* \sa cartesian.h for typedefs and defines*/
|
||||
class Quaternion {
|
||||
public:
|
||||
Quaternion();
|
||||
~Quaternion();
|
||||
|
||||
double w; ///<Real Component
|
||||
double i; ///<Imaginary Component
|
||||
double j; ///<Imaginary Component
|
||||
double k; ///<Imaginary Component
|
||||
|
||||
double &operator[](unsigned index);
|
||||
|
||||
/// Copy Constructor
|
||||
Quaternion(const Quaternion &two);
|
||||
/// Assignment constructor
|
||||
Quaternion(double inw, double ini, double inj, double ink);
|
||||
/// Axis-Angle Rotation (\b v must be a \b UNIT vector)
|
||||
Quaternion(const vectorPt &v, double angle);
|
||||
/// Rotation from vector 1 to vector 2 (v1 and v2 need not be normalized)
|
||||
/** Angle is between \b v1 and \b v2 and axis is calculated as the cross
|
||||
* product of \b v1 and \b v2 */
|
||||
Quaternion(const vectorPt &v1, const vectorPt &v2);
|
||||
/// Spherical rotation
|
||||
Quaternion(double longitude, double latitude);
|
||||
/// From Euler angles
|
||||
Quaternion(const EulerRot &erot);
|
||||
|
||||
void operator += (const Quaternion &two);
|
||||
Quaternion operator * (const Quaternion &two) const;
|
||||
void operator*= (const Quaternion &two);
|
||||
Quaternion operator * (const double two) const;
|
||||
void operator*= (const double two);
|
||||
void operator= (const Quaternion &two);
|
||||
bool operator== (const Quaternion &two) const;
|
||||
|
||||
/// Returns the conjugate
|
||||
Quaternion conj() const;
|
||||
/// Returns the norm
|
||||
double norm();
|
||||
/// Form \b *this into the unit vector
|
||||
void normalize();
|
||||
/// Returns the unit vector of \b *this
|
||||
Quaternion unit();
|
||||
|
||||
friend ostream & operator<<(ostream &out, const Quaternion &q);
|
||||
friend istream & operator>>(istream &in, Quaternion &q);
|
||||
};
|
||||
|
||||
/// Rotation matrix (about origin)
|
||||
/** Rotation of ThreeD can be applied via overloaded * operator **/
|
||||
class RotMat {
|
||||
public:
|
||||
/// Unspecified matrix!
|
||||
RotMat();
|
||||
/// Initialize with quaternion
|
||||
RotMat(const Quaternion &q);
|
||||
/// Set based on quaternion
|
||||
void set(const Quaternion &q);
|
||||
/// Assert that the rotation is proper
|
||||
void proper();
|
||||
|
||||
cPt operator*(const cPt &in) const;
|
||||
|
||||
double x_x,x_y,x_z,y_x,y_y,y_z,z_x,z_y,z_z;
|
||||
};
|
||||
|
||||
///\def NOROTATE
|
||||
/// Quaternion with no rotation
|
||||
#define NOROTATE Quaternion(1.0,0.0,0.0,0.0)
|
||||
|
||||
/// Global geometry functions
|
||||
namespace c {
|
||||
/// Returns point on a line closest to an outside point
|
||||
/** The line is described by a directional vector and a point on the line
|
||||
*\param v Vector describing the line direction
|
||||
*\param v_point Point on the line
|
||||
*\param point The Outside point */
|
||||
cPt point_to_line(const vectorPt &v, const cPt &v_point,
|
||||
const cPt &point);
|
||||
|
||||
/// Return the closest distance between two line segments input as end points
|
||||
/**\param l1_1 End point of segment 1
|
||||
*\param l1_2 End point of segment 1
|
||||
*\param l2_1 End point of segment 2
|
||||
*\param l2_2 End point of segment 2 */
|
||||
double closest_approach(const cPt &l1_1, const cPt &l1_2,
|
||||
const cPt &l2_1, const cPt &l2_2);
|
||||
/// Calculates the points where two line segments are closest
|
||||
/**\param l1_1 End point of segment 1
|
||||
*\param l1_2 End point of segment 1
|
||||
*\param l2_1 End point of segment 2
|
||||
*\param l2_2 End point of segment 2
|
||||
*\param close_l1 Calculated closest point on line segment 1
|
||||
*\param close_l2 Calculated closest point on line segment 2 */
|
||||
void closest_approach_points(const cPt &l1_1, const cPt &l1_2,
|
||||
const cPt &l2_1, const cPt &l2_2,
|
||||
cPt &close_l1, cPt &close_l2);
|
||||
|
||||
/// Returns true if two line segments intersect
|
||||
bool intersect(const c2DPt &line1_start, const c2DPt &line1_end,
|
||||
const c2DPt &line2_start, const c2DPt &line2_end);
|
||||
|
||||
/// Liang-Barsky Intersection between line segment and line
|
||||
bool sline_intersect(const c2DPt &line1_start, const c2DPt &line1_end,
|
||||
const c2DPt &line2normal, const c2DPt &line2_point);
|
||||
/// Average position
|
||||
cPt mean(vector<cPt> &vec);
|
||||
}
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user