Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F74864242
cartesian.h
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Jul 30, 02:49
Size
15 KB
Mime Type
text/x-c
Expires
Thu, Aug 1, 02:49 (2 d)
Engine
blob
Format
Raw Data
Handle
19352519
Attached To
rLAMMPS lammps
cartesian.h
View Options
/****************************************************************************
* 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);
inline numtyp &operator[](unsigned i) {
switch(i) {
case X: return x;
case Y: return y;
case Z: return z;
}
return x;
}
inline numtyp operator[](unsigned i) const {
switch(i) {
case X: return x;
case Y: return y;
case Z: return z;
}
return x;
}
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);
inline void operator *= (const numtyp &two) {
x*=two; y*=two; z*=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
inline numtyp dist(const ThreeD<numtyp> &two) {
return (*this-two).norm();
}
/// Distance squared between two points
numtyp dist2(const ThreeD<numtyp> &two);
/// Converts \b *this to the unit vector
inline void normalize() {
numtyp temp=norm();
#ifdef NANCHECK
assert(temp!=0);
#endif
*this/=temp;
}
/// 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
Event Timeline
Log In to Comment