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 (1 d, 22 h)
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