CHAPTER 3

Vector and Matrix

 

3-dimensional vector and matrix are frequently used in crystallography, for example, atomic coordinates are three numbers of a vector; a rotation matrix has 3 ´ 3 elements.  Two template classes mac::vector3D<T> and mac::matrix3D<T> are constructed to handle operations of vector and matrix.  These two classes are also modified slightly in order to be shadowed into Python by SWIG.

 

3.1     Vectors mac::vector3D<T> and cpl.mac.vm.vector3DTYP

 

3.1.1 The template class mac::vector3D<T>

 

This class mac::vector3D<T> represents the most widely used data type in crystallography due to the 3-dimensional nature of this science.  Although historically this class was made after the non-template class vector3 (obsolete),  it is now the primary implementation of 3-dimensional vectors.

 

Besides the basic operations of vector, several geometric functions are also implemented in this class, such as length of a vector, distance and angle between two vectors, projection, orientation, area, and volume.  It must be warned that these geometric functions are under a common assumption, that is, orthogonal coordinate system (see Chapter 12 Unit Cell and Coordinate Systems).  These geometric functions will not perform correctly, if the vector involved is under a non-orthogonal system, e.g., a crystallographic coordinate system.  These functions offer no coordinate system checking other than this warning.

 

3.1.1.1 Constructors

vector3D(const T = 0, const T = 0, const T = 0);

The 3-dimensional vector we refer to here contains only three numerical values of the same type.  These vectors cannot contain strings or other non-numerical types.  The general constructor takes no more than three values of the same type, and constructs a vector.  The missing values will be filled by the default 0, where T can only be a numerical type.

 

Another constructor taking the length and orientation of the vector may be added in the future release.

 

Default constructor

vector3D(const T u = 0,

         const T v = 0,

         const T w = 0);

The default constructor mac::vector3D<T>() simply constructs a vector of three elements of 0.

 

Copy constructor

vector3D(const vector3D<T>& v);

The copy constructor returns an identical vector to v.

 

Constructor from a polar orientation

vector3D(const std::complex<T>& p,

         const T s = 1);

Constructs a vector along an orientation specified by polar angles p (see 3.1.1.14 for definition of polar angles).  The newly constructed vector has a length of s.

 

3.1.1.2 Destructor

~vector3D();

There is nothing to do in the destructor.

 

3.1.1.3 Conversion

operator vector3         () const;

operator vector3D<double>() const;

operator vector3D<int>   () const;

Operators converting a vector to mac::vector3D<double>, mac::vector3D<int>, and the obsolete vector3 are available in some releases.  When vector type mac::vector3D<T> is supplied where those vector types are required, auto-conversion will take place.  However, the conversion may be associated with loosing precision.

 

3.1.1.4 Element access

inline T  operator[](const int i) const;

inline T& operator[](const int i);

inline T  operator()(const int i) const;

inline T& operator()(const int i);

The operators [] and () are overloaded in two forms, const and non-const.  The const form provides a means to get the element value.  The non-const form provides an l-value.  The operator [] accepts 0-offset subscript (C style), and () accepts 1-offset subscript (Fortran style).  For example, one may write the following code to test these overloaded element accessing operators.  This program shall yield no message, if everything is all right.

 

#include <assert.h>

#include <mac/vector3D.h>

 

main()

{

   mac::vector3D<int> v;

// mac::vector3D<int> v(); Error when compiled, since v is considered as

//                         a pointer-to-function.

 

   v[0] = 11;             // Use the non-const form of operator [] as

   v[1] = -5;             // l-values.

   v(3) = 28;             // Use the non-const form of operator ().

   assert (v[0] == v(1)); // Use the const form of operators [] and ().

   assert (v[1] == v(2));

   assert (v[2] == v(3));

 

   const mac::vector3D<int> u(v); // Copy constructor.

   assert (u == v);

 

   mac::vector3D<double> w(M_PI / 2 * std::complex<double>(1, 1), 2);

   assert (w &= mac::vector3D<double>(0, 0, -2));

 

   return 0;

}

 

Listing 3.1.1.4.1 Testing element access by const and non-const forms of operators [] and ().

 

Explicit element access

T getX() const;

T getY() const;

T getZ() const;

Three functions are also available to get the elements explicitly.

 

3.1.1.5 Assignment

const vector3D<T> &operator=(const vector3D<T>&);

Assignment operator allows assignment such as u = v.

 

3.1.1.6 Addition and subtraction

vector3D<T> operator+ (const vector3D<T>&) const;

vector3D<T> operator- (const vector3D<T>&) const;

vector3D<T> operator+=(const vector3D<T>&);

vector3D<T> operator-=(const vector3D<T>&);

vector3D<T> operator+() const;

vector3D<T> operator-() const;

Operators +, -, +=, and -= are overloaded to perform the common addition and subtraction.  The binary operators + and ¨C obviously handle vector addition u + v and subtraction u ¨C v, respectively.  The unary operator + maintains the sign of a vector +v, and ¨C inverts the sign ¨Cv.  The operators += and -= perform the anticipated vector addition u += v and subtraction u -= v, respectively.

 

3.1.1.7 Scaling

Operators *, /, *=, and /= are overloaded to perform isotropic and anisotropic scaling of a vector.

 

Isotropic scaling

Isotropic scaling means that all three elements of a vector are scaled by a common factor.

 

vector3D<T> scale(const T) const;

Returns the vector scaled by the given scalar number.

 

vector3D<T> operator* (const T) const;

vector3D<T> operator/ (const T) const;

These two overloaded binary operators * and / take one argument of a value as the isotropic scale factor.  Thus, v * a and v / a scale the vector v by the factors a and 1/a, respectively.  In case of operator /, if a is 0, an assertion error will occur.

 

template <class T>

vector3D<T> operator*(const T& x,

                      const vector3D<T>& v);

An unbounded function allows x * v has the same effect as v * x.  Obviously, writing x / v shall and will have very different effect from v / x.  See Anisotropic scaling below for detail.

 

vector3D<T> operator*=(const T);

vector3D<T> operator/=(const T);

Operators *= and /= are overloaded to perform isotropic scaling as well.  v *= a and v /= a return the scaled vectors by factors a and 1/a, respectively.  Again, in case if operator /=, a zero a will cause an assertion error.

 

Anisotropic scaling

If the three elements of a vector are scaled by different factors, we say this scaling is anisotropic.

 

vector3D<T> scale(const vector3D<T>&) const;

Returns the vector scaled by given vector.

 

vector3D<T> operator/(const vector3D<T>&) const;

The binary operator * that connects two vectors is reserved for dot multiplication (See 3.1.1.8), but / can be used for anisotropic scaling.  u / v returns a vector (u.x/v.x, u.y/v.y, u.z/v.z).  If any of the three elements of v is 0, an assertion error will occur.

 

template <class T>

vector3D<T> operator/(const T& a,

                      const vector3D<T>& v);

An unbounded function overloads the operator / to connect a scalar value with a vector a / v, and returns a vector (a/v.x, a/v.y, a/v.z).  If any of the three elements of v is 0, an assertion error will occur.

 

vector3D<T> operator*=(const vector3D<T>&);

vector3D<T> operator/=(const vector3D<T>&);

Operators *= and /= are overloaded to perform anisotropic scaling.  u *= v returns a vector (u.x*v.x, u.y*v.y, u.z*v.z), and u /= v returns (u.x/v.x, u.y/v.y, u.z/v.z).  Also in the latter case, if any of the three elements of v is 0, an assertion error will occur.

 

The following program tests these operators:

 

#include <assert.h>

#include <mac/vector3D.h>

 

main()

{

   typedef mac::vector3D<int> V;

 

   // Isotropic scaling

   assert (V(3, 0, -2) * 4 == V(12, 0, -8));

   assert (4 * V(3, 0, -2) == V(12, 0, -8));

   assert (V(6, 0, -2) / 2 == V( 3, 0, -1));

 

   V u(6, 0, -2);

   u *= 3;

   assert (u == V(18, 0, -6));

   u /= 2;

   assert (u == V( 9, 0, -3));

 

   // Anisotropic scaling

   assert (V(9, 0, -2) / V(-3, 1, 2) == V(-3, 0, -1));

   assert (6           / V(-3, 1, 2) == V(-2, 6,  3));

 

   u *= V(-1, 1, 2);

   assert (u == V(-9, 0, -6));

   u /= V(3, 2, 1);

   assert (u == V(-3, 0, -6));

 

   return 0;

}

 

Listing 3.1.1.7.1 Testing isotropic and anisotropic scaling.

 

3.1.1.8 Dot (scalar) multiplication

T operator*(const vector3D<T>&) const;

T dot(const vector3D<T>&) const;

Operator * is overloaded to perform dot or scalar multiplication, so does the function T dot(const vector3D<T> &v) const.  Both u * v and u.dot(v) return the dot product of the vectors u and v.

 

static T dot(const vector3D<T>&,

             const vector3D<T>&);

template <class T>

T dot_prod(const vector3D<T>&,

           const vector3D<T>&);

A static function and an unbounded function are also available to calculate the dot product.

 

3.1.1.9 Cross (vector) multiplication

static vector3D<T> cross(const vector3D<T>&,

                         const vector3D<T>&);

vector3D<T> cross(const vector3D<T>&) const;

vector3D<T> operator->*(

                  const vector3D<T>&) const;

Similar to the dot multiplication, a member function cross, a static function cross and the overloaded operator ->* all perform cross or vector multiplication.  One may get vector product of two vectors u and v by writing u.cross(v), mac::vector3D<T>::cross(u, v), or simply u ->* v.

 

3.1.1.10 Pre-multiplication by a matrix

vector3D<T> operator*(

   const mac::matrix3D<T>&) const;

The operator * is also overloaded to perform pre-multiplication by a matrix such as v * r, v and r are a vector and a matrix mac::matrix3D<T> (see 3.2.1), respectively.  However, written in this way it is inconsistent with the mathematical convention.

 

template <class T>

vector3D<T> operator*(const matrix3D<T>& m,

                      const vector3D<T>& v);

An unbounded function lets you write r * v for pre-multiplication of a vector v by a matrix r.

 

vector3D<T> operator*=(const mac::matrix3D<T>&);

vector3D<T> operator/=(const mac::matrix3D<T>&);

The other two operators *= and /= are also overloaded to perform pre-multiplication.  v *= r returns the vector that is v pre-multiplied by matrix r, and v /= r returns the vector that is v pre-multiplied by the inversed matrix of r.

 

This program tests pre-multiplication:

 

#include <assert.h>

#include <mac/vector3D.h>

 

main()

{

   typedef mac::vector3D<int> V;

   typedef mac::matrix3D<int> M;

 

   assert (V(3, 0, -2) * M(0,0,1, 0,1,0, 1,0,0) == V(-2, 0, 3));

      // Inconsistent with the mathematical convention.

   assert (M(0,0,1, 0,1,0, 1,0,0) * V(3, 0, -2) == V(-2, 0, 3));

      // Mathematical convention.  Either way is the same.

 

   V u(6, 0, -2);

   u *= M(0,0,1, 0,1,0, 1,0,0);

   assert (u == V(-2, 0, 6));

   u /= M(0,0,1, 0,1,0, 1,0,0);  // Pre-multiplied by the inversed matrix.

   assert (u == V(6, 0, -2));

 

   return 0;

}

 

Listing 3.1.1.10.1 Testing pre-multiplication by a matrix.

 

3.1.1.11 Length

double vabs() const;

T vnorm() const;

The overloaded functions vabs and vnorm mimic abs, fabs, and norm.  The function T vnorm() const returns the square of the vector length, and the function double vabs() const returns the vector length.  The return type of vabs is double since a square-root must be taken.  Therefore, precision may be lost when calculating the length of a long double vector by vabs.  You should use vnorm instead to avoid precision loss.

 

static double vabs(const T, const T, const T);

static double vabs(const vector3D<T>&);

static T vnorm(const T, const T, const T);

static T vnorm(const vector3D<T>&);

Several static functions are available to calculate vector length and its square.  Caution must be made to guard from precision loss.

 

template <class T> double  abs(

   const mac::vector3D<T>&);

template <class T> double  abs(const T&,

                               const T&,

                               const T&);

template <class T> double norm(

   const mac::vector3D<T>&);

template <class T> double norm(const T&,

                               const T&,

                               const T&);

template <class T> double vabs(

   const mac::vector3D<T>&);

template <class T> T     vnorm(

   const mac::vector3D<T>&);

template <class T> double vabs(const T&,

                               const T&,

                               const T&);

template <class T> T     vnorm(const T&,

                               const T&,

                               const T&);

Several forms of unbounded functions vabs, vnorm, abs, and norm are also available.  See the following example for different way of getting the vector length.

 

#include <assert.h>

#include <mac/vector3D.h>

 

main()

{

   typedef mac::vector3D<int> V;

 

   // There are many ways to calculate the vector length.

   assert (V(3, 0, -4).vabs()    == 5);

   assert (V::vabs(3, 0, -4)     == 5);

   assert (V::vabs(V(3, 0, -4))  == 5);

   assert (vabs(3, 0, -4)        == 5);

   assert (vabs(V(3, 0, -4))     == 5);

   assert (abs(3, 0, -4)         == 5);

   assert (abs(V(3, 0, -4))      == 5);

 

   assert (V(3, 0, -4).vnorm()   == 25);

   assert (V::vnorm(3, 0, -4)    == 25);

   assert (V::vnorm(V(3, 0, -4)) == 25);

   assert (vnorm(3, 0, -4)       == 25);

   assert (vnorm(V(3, 0, -4))    == 25);

   assert (norm(3, 0, -4)        == 25);

   assert (norm(V(3, 0, -4))     == 25);

 

   return 0;

}

 

Listing 3.1.1.11.1 Testing various forms of vabs, vnorm, abs, and norm.

 

3.1.1.12 Distance

static double distance(const vector3D<T> &u,

                       const vector3D<T> &v);

double distance(const vector3D<T> &v) const;

double operator|(const vector3D<T> &v) const;

A member function distance and a static function distance calculate the distance between one vector u to another v.  The operator | is also overloaded to calculate the distance.  u.distance(v), mac::vector3D<T>::distance(u, v), and u | v all return the distance between vectors u and v.

 

3.1.1.13 Projection

If there are two vectors in 3-dimensional space, one may always find the projection of one vector on the other.  Several functions are available to calculate the length of the projection.  The sign of the result is determined by whether the projection is along the direction of the vector that is projected on.

 

double projectionOn(const vector3D<T>& v) const;

complex<double> projectionOn(

   const vector3D<T>&,

   const vector3D<T>&) const;

vector3D<double> projectionOn(

   const vector3D<T>&,

   const vector3D<T>&,

   const vector3D<T>&) const;

u.projectionOn(v) returns a double as the projection length of vector u on vector v.  u.projectionOn(v, w) returns a complex<double> that contains the projection lengths of vector u on vectors v and w, respectively.  Similarly, t.projectionOn(u, v, w) returns a vector3D<double> that contains the projection lengths of vector t on vectors u, v, and w, respectively.

 

double projectionOf(const vector3D<T>& v) const;

u.projectionOf(v) returns a double as the projection length of vector v on vector u.

 

static double projection(const vector3D<T>& u,

                         const vector3D<T>& v);

The static function below returns the projection length of vector u on vector v.

 

All these functions above may have precision loss when dealing with vector of long double.

 

static T projectionSq(const vector3D<T>&,

                      const vector3D<T>&);

If working with high precision vector such as vector3D<long double>, another static function projectionSq should be used to avoid precision loss.  This static function returns the squared projection length of vector u on vector v.  The return type maintains as T, therefore no precision loss occurs.  However, if T is an integer type, the result is also that integer type, which is usually not what the user would intend.  Do not use this function with any integer type.  Please also note that the sign of the result follows the same rule as the other functions above, indicating whether the direction of the projection is along the direction of the vector projected on.

 

If attempting to project on a 0-length vector that defines no orientation, an assertion error will occur at runtime.  The following program tests these functions:

 

#include <assert.h>

#include <mac/vector3D.h>

 

main()

{

   typedef mac::vector3D<double> U;

   typedef mac::vector3D<int>    V;

 

   const V O(0, 0, 0); // Origin

   const V X(1, 0, 0); // X-axis

   const V Y(0, 1, 0); // Y-axis

   const V Z(0, 0, 1); // Z-axis

 

   assert (O.projectionOn(V(-6, 0, 8)) == 0);

// assert (V(3, 0, -4).projectionOn(O) == -5);

      // Error.  Cannot project on origin.

 

   assert (V(3, 0, -4).projectionOn(V(-6, 0, 8)) == -5);

   assert (V(3, 0, -4).projectionOf(V(-6, 0, 8)) == -10);

 

   assert (V(3, 0, -4).projectionOn(X) ==  3);

   assert (V(3, 0, -4).projectionOn(Z) == -4);

   assert (real(V(3, 0, -4).projectionOn(X, Z)) ==  3);

   assert (imag(V(3, 0, -4).projectionOn(X, Z)) == -4);

   assert (V(3, 0, -4).projectionOn(X, Y, Z)(1) ==  3);

   assert (V(3, 0, -4).projectionOn(X, Y, Z)(2) ==  0);

   assert (V(3, 0, -4).projectionOn(X, Y, Z)(3) == -4);

 

   assert (V::projection(V(3, 0, -4), V(-6, 0, 8)) == -5);

   assert (V::projectionSq(V(3, 2, 4), V(1, 2, 3)) <=

           U::projectionSq(U(3, 2, 4), U(1, 2, 3)));

      // The squared projection length returned by projectionSq with T being

      // integer types may be smaller than the true value, ...

   assert (V::projectionSq(V(3, 2, 4), V(-1, -2, -3)) >=

           U::projectionSq(U(3, 2, 4), U(-1, -2, -3)));

      // ... or greater depending on the direction of the projection.

 

   return 0;

}

 

Listing 3.1.1.13.1 Testing functions that calculate projection length.

 

3.1.1.14 Orientation

An orientation in 3-dimensional space can be described by three angles between this orientation and the three orthogonal axes.  These three angles are called the orientation angles, and their cosine values are called the orientation cosine.

 

vector3D<double> vcos() const;

vector3D<double> vsin() const;

vector3D<double> vtan() const;

vector3D<double> orientationAngle() const;

vcos() and orientationAngle() return the orientation cosine and angles in radian in the form of vector3D<double>, respectively.  vsin() and vtan() are defined to calculate the sine and tangent values of the orientation angles.  They all return vector3D<double>.

 

template <class T>

vector3D<double> sin(const vector3D<T>&);

template <class T>

vector3D<double> cos(const vector3D<T>&);

template <class T>

vector3D<double> tan(const vector3D<T>&);

Unbounded functions sin(), cos(), and tan() are overloaded to calculate the orientation sine, cosine, and tangent of a vector, respectively.

 

Polar angles

The three values of orientation cosine are not independent.  They must satisfy cos2a + cos2b + cos2g = 1.  Polar angles y and j are often used in crystallography to define an orientation.  The angle y can be thought equivalent to the latitude but starting from the north pole as 0, and ending at the south pole as p.  The angle j is equivalent to the longitude ranging from 0 to 2p.  There are different conventions when matching the latitude-longitude system with the orthogonal axes.  This is the origin of the most frequent errors.  Caution must be taken.  MAC sets Y-axis at the north pole, that is, y = 0, and X-axis at j = 0.  Z-axis may then be determined by the right-hand coordinates system.  The positive direction of j is defined by the right-hand rule about the Y-axis, which is consistent with the geographic convention.  Under this convention, a vector (x, y, z) in orthogonal coordinate system and its orientation¡¯s polar angles y and j have this relationship:

 

x =  sinycosj,

y =  cosy,

z = -sinysinj.

 

complex<double> polarAngle() const;

This member function returns the polar angles y and j in radian as the real and imaginary parts of the complex under the above definition.

 

All the functions calculating orientation cosine, orientation angles, and polar angles may have precision loss, when dealing with long double.  There is no functions in the current release that maintains the precision.  Such functions may be added in the future release.  Obviously, a 0-length vector defines no orientation.  All the functions will cause assertion error, when finding the orientation of a 0-length vector.

 

Parallel and perpendicular

bool parallel     (const vector3D<T> &v) const;

bool antiparallel (const vector3D<T> &v) const;

bool perpendicular(const vector3D<T> &v) const;

bool operator||(const vector3D<T> &v) const;

bool operator&&(const vector3D<T> &v) const;

Three functions test whether the given vector v is parallel, antiparallel or perpendicular.  u.parallel(v), u.antiparallel(v), and u.perpendicular(v) return Booleans.  The operator || is also overloaded to test parallelism, e.g., u || v.  u && v returns true if these two vectors are anti-parallel.  Any attempts involving 0-length vector will cause assertion error.

 

The following program tests these functions:

 

#include <assert.h>

#include <mac/vector3D.h>

 

main()

{

   typedef mac::vector3D<double> V;

 

   const V O(0, 0, 0); // Origin

   const V X(1, 0, 0); // X-axis

   const V Y(0, 1, 0); // Y-axis

   const V Z(0, 0, 1); // Z-axis

 

// O.vcos();             // Error

// O.orientationAngle(); // Error

 

   assert (X.vcos() == X);

   assert (Y.vcos() == Y);

   assert (Z.vcos() == Z);

   assert (X.orientationAngle() == V(   0, M_PI, M_PI) / 2);

   assert (Y.orientationAngle() == V(M_PI,    0, M_PI) / 2);

   assert (Z.orientationAngle() == V(M_PI, M_PI,    0) / 2);

 

   const double ang = acos(1 / sqrt(3.)); // 54.73561032 degrees

   assert ((X - Y + Z).vsin() &= V(1,  1, 1) * sqrt(2 / 3.));

   assert ((X - Y + Z).vcos() &= V(1, -1, 1) / sqrt(3.));

   assert ((X - Y + Z).vtan() &= V(1, -1, 1) * sqrt(2.));

   assert ((X - Y + Z).orientationAngle() &= V(ang, M_PI - ang, ang));

   assert (sin(X - Y + Z) &= V(1,  1, 1) * sqrt(2 / 3.));

   assert (cos(X - Y + Z) &= V(1, -1, 1) / sqrt(3.));

   assert (tan(X - Y + Z) &= V(1, -1, 1) * sqrt(2.));

      // Operator &= tests the approximate equality.

 

   assert (real(V(3, -sqrt(3.), 0).polarAngle()) - 2 * M_PI / 3 < FLT_EPSILON);

      // Psi

   assert (imag(V(3, -sqrt(3.), 0).polarAngle()) == 0); // Phi

 

   assert ( V(1, 2, 3).parallel(V( 2,  4,  6)));

   assert ( V(1, 2, 3) || V( 3,  6,  9));

   assert ( V(1, 2, 3).antiparallel(V(-2, -4, -6)));

   assert (!V(1, 2, 3).parallel(V(-2, -4,  6)));

 

   assert (X.perpendicular(Y));

   assert (Y.perpendicular(Z));

   assert (Z.perpendicular(X));

 

// assert (X || O); // Error

// assert (O || X); // Error

 

   return 0;

}

 

Listing 3.1.1.14.1 Testing various functions that calculate the orientation of a vector.

 

3.1.1.15 Angle

double angle(const vector3D<T> &v) const;

static double angle(

   const vector3D<T>& u,

   const vector3D<T>& v,

   const vector3D<T>& w = mac::vector3D<T>());

double operator^(const vector3D<T> &v) const;

A member function and a static function calculate the angle in radian between two vectors.  The operator ^ returns the angle too, e.g., u ^ v.  The third, optional vector3D<T> w defines the sign of the angle.  If the projection of u ´ v on w has the same direction as w, the angle between u and v is positive.  This is consistent with the right-hand rule.  If the optional vector is omitted, the returned angle is always non-negative.

 

static double cosAngle(const vector3D<T>&,

                       const vector3D<T>&);

double cosAngle(const vector3D<T>&) const;

static double sinAngle(

   const vector3D<T>&,

   const vector3D<T>&,

   const vector3D<T>& = mac::vector3D<T>());

double sinAngle(const vector3D<T>&) const;

static double tanAngle(

   const vector3D<T>&,

   const vector3D<T>&,

   const vector3D<T>& = mac::vector3D<T>());

double tanAngle(const vector3D<T>&) const;

The other functions sinAngle, cosAngle, and tanAngle calculate the corresponding values.  Each has the static and non-static forms.  If the third, optional vector is available, the sign of returned value is defined.  See above for definition.  Cosine function does not need to define the sign, since it is symmetric.  All these function may have precision loss, and no effort is made to maintain high precision such as long double in the current release.

 

3.1.1.16 Area

double triangleVertices(

   const vector3D<T>&, const vector3D<T>&) const;

static double triangleVertices(

                        const vector3D<T>&,

                        const vector3D<T>&,

                        const vector3D<T>&);

A triangle in 3-dimensional space can be defined by three vectors at the vertices of the triangle.  The triangular area is returned these functions.

 

double triangleEdges(const vector3D<T> &v) const;

static double triangleEdges(

   const vector3D<T> &u, const vector3D<T> &v);

These functions return the triangular area defined by two vectors as two of its edges. 

 

3.1.1.17 Volume

Scalar triple product

T scalarTripleProduct(

   const vector3D<T> &u,

   const vector3D<T> &v) const;

static T scalarTripleProduct(const vector3D<T>&,

                             const vector3D<T>&,

                             const vector3D<T>&);

Scalar triple product of vectors u, v, w is defined by u ¡¤ (v ´ w).  Two functions calculate this value.

 

Volume of parallelepiped

T parallelepipedVolume(const vector3D<T>&,

                       const vector3D<T>&) const;

static T parallelepipedVolume(

   const vector3D<T>&,

   const vector3D<T>&,

   const vector3D<T>&);

A parallelepiped can be defined by three vectors.  Its volume can be calculated by these two functions.  A unit cell volume can be easily calculated by these functions (See Chapter 4).

 

3.1.1.18 Location

bool  inside(const vector3D<T> &u,

             const vector3D<T> &v) const;

bool outside(const vector3D<T> &u,

             const vector3D<T> &v) const;

These two functions return whether the vector is inside or outside the parallelepiped whose two body diagonal vertices are u and v, respectively.  If a vector is right on the surface of the parallelepiped defined by u and v, it is considered inside.

 

3.1.1.19 Comparison

Equality and inequality

bool operator==(const vector3D<T> &v) const;

bool operator!=(const vector3D<T> &v) const;

Two vectors are equal, if and only if all three elements of one vector are equal to those of the other, respectively.  They are unequal, if and only if they are not equal.  The operators == and != are overloaded to test equality and inequality.

 

bool operator&=(const vector3D<T> &v) const;

In addition, an approximate equality is defined as the distance between two vectors is less than a small value, e.g., 10 times of FLT_EPSILON, which is about 10-6.  The operator &= is used to test approximate equality.  This release lacks the ability to change the threshold of approximate equality, which is likely to be added in a later release.

 

Less than

static bool lessThan(const vector3D<T>&,

                     const vector3D<T>&);

class less {public: bool operator()(

   const vector3D<T>&,

   const vector3D<T>&) const;};

The most important comparison method is ¡®less than¡¯, which is often used by sorting and sorted containers.  A static function lessThan and a functional object or functor less are defined.  The earlier element of the vector has higher priority to determine the comparison result.  For two vectors u and v, the relationship of u.x and v.x determines that of u and v, unless u.x == v.x, then the relationship of u.y and v.y determines that of u and v, unless u.y == v.y, then the relationship of u.z and v.z determines that of u and v.  There is no other functor such as equal or greater defined yet.

 

bool operator> (const vector3D<T>&) const;

bool operator< (const vector3D<T>&) const;

bool operator>=(const vector3D<T>&) const;

bool operator<=(const vector3D<T>&) const;

The operator < is defined in the same sense as the static function lessThan().  With the operators == and < defined, the operators <=, >, and >= are then defined accordingly.

 

static bool shorterThan(const vector3D<T>&,

                        const vector3D<T>&);

bool shorterThan(const vector3D<T>&) const;

class shorter {public: bool operator()(

   const vector3D<T>&,

   const vector3D<T>&) const;};

Another set of comparison is based on the length of the vector.  For two vectors u and v, the relationship of u.vnorm and v.vnorm determines that of u and v, unless u.vnorm == v.vnorm, then the previous definition takes effect.

 

0-length

bool operator!() const;

The operator ! is defined to test 0-length, therefore !v is true if v is 0 in length, and !!v is true if not.

 

The following program tests vector comparison:

 

#include <assert.h>

#include <mac/vector3D.h>

 

main()

{

   typedef mac::vector3D<int> V;

   V O(0, 0, 0); // 0-length vector

   V X(1, 0, 0); // X-axis

 

   assert (V::lessThan(V(1, 2, 3),  V(2, 6, 3)));

   assert (            V(1, 2, 3) < V(2, 6, 3));

      // When the first vector is shorter than the second, ...

 

   assert (V::lessThan(V(1, 2, 3),  V(2, 1, 3)));

   assert (            V(1, 2, 3) < V(2, 1, 3));

      // or they are equal in length, the earlier element determines the

      // comparison.

 

   assert (V::lessThan(V(1, 2, 3),  V(2, 1, 1)));

   assert (          !(V(1, 2, 3) < V(2, 1, 1)));

      // When the first vector is longer than the second, lessThan and the

      // operator < behave differently.

 

   assert (!O);  // 0-length

   assert (!!X); // non-0-length

 

   return 0;

}

 

Listing 3.1.1.19.1 Testing vector comparison.

 

3.1.1.20 Output

friend ostream& operator<<(ostream&,

                           const vector3D<T>&);

The operator << is defined to print three elements separated by commas (,) to ostream.

 

void write(ofstream*) const;

void write(const char* =

           message::typeGeneral) const;

Two write functions are defined to work with the class message of TUI.  See message for detail.  It is not yet finalized whether MAC should depend on TUI or be isolated.  The second write function may not be available in some releases.

 

The following program tests the rest of the functionalities that are not included in the examples above.

 

#include <assert.h>

#include <mac/vector3D.h>

#include <iostream.h>

 

main()

{

   typedef mac::vector3D<double> U;

   typedef mac::vector3D<int>    V;

 

// const V O();        // Error, since the compiler confuses with a function

                       // declaration.

// const V O = V();    // OK.

   const V O;          // Origin, by default constructor

   const V X(1);       // X-axis, the missing two arguments are 0

   const V Y(0, 1);    // Y-axis, the missing argument is 0

   const V Z(0, 0, 1); // Z-axis, by the general constructor

   cout << "Origin " << O << endl;

   cout << "X-axis " << X << endl;

   cout << "Y-axis " << Y << endl;

   cout << "Z-axis " << Z << endl;

   cout << endl;

 

   U u(1.2, 3.4, -2.7);

   V v = u; // vector3D<double> is converted and assigned to vector3D<int>

   V w(u);  // vector3D<int> w is a copy of vector3D<double> u, but implicit

            // conversion will happen.

   cout << "A vector3D<double>                              " << u << endl;

   cout << "A vector3D<int> converted from vector3D<double> " << v << endl;

   cout << "A vector3D<int> copied    from vector3D<double> " << w << endl;

   cout << endl;

 

   cout << "(" << v << ") + (" << X << ") is " << v + X << endl;

   cout << "(" << v << ") - (" << Y << ") is " << v - Y << endl;

   cout << "+(" << v << ") is " << +v << endl;

   cout << "-(" << v << ") is " << -v << endl;

   cout << "v is " << v << endl;

   cout << "v += (" << X + Y + Z << ") is " << (v += (X + Y + Z)) << endl;

   cout << "v -= (" << X + Y + Z << ") is " << (v -= (X + Y + Z)) << endl;

   cout << endl;

 

   assert (X ->* Y == Z);

   cout << "(" << v << ") . (" << Y << ") is " << v  *  Y << endl;

   cout << "(" << v << ") X (" << Y << ") is " << v ->* Y << endl;

   cout << endl;

 

   cout << "From O to " << V(3, 4, 0) << " is " << (O | V(3, 4, 0)) << endl;

   cout << endl;

 

   cout << "Angle between " << X << " and " << V(2, 1, 0) << " is "

        << (X ^ V(2, 1, 0)) << endl;

   cout << "Sine    of the angle is " << V::sinAngle(X, V(2, 1, 0)) << endl;

   cout << "Cosine  of the angle is " << V::cosAngle(X, V(2, 1, 0)) << endl;

   cout << "Tangent of the angle is " << V::tanAngle(X, V(2, 1, 0)) << endl;

   cout << endl;

 

   const V v1(0, 0, 1);

   const V v2(3, 0, 1);

   const V v3(3, 4, 1);

   cout << "The area of triangle with three vertices (" << v1 << ")("

                                                        << v2 << ")("

                                                        << v3 << ") is "

        << V::triangleVertices(v1, v2, v3) << endl;

   cout << "The area of triangle with two edges (" << V(3) << ")("

                                                   << V(0, 4) << ") is "

        << V::triangleEdges(V(3), V(0, 4)) << endl;

   cout << "The volume of parallelepiped with three edges (" << X << ")("

                                                             << 2*Y << ")("

                                                             << 3*Z << ") is "

        << V::parallelepipedVolume(X, 2*Y, 3*Z) << endl;

 

   assert (O.inside(X, Y));

   assert (O.outside(-X, -X-Y-Z));

 

   return 0;

}

 

Origin 0,0,0

X-axis 1,0,0

Y-axis 0,1,0

Z-axis 0,0,1

 

A vector3D<double>                              1.2,3.4,-2.7

A vector3D<int> converted from vector3D<double> 1,3,-2

A vector3D<int> copied    from vector3D<double> 1,3,-2

 

(1,3,-2) + (1,0,0) is 2,3,-2

(1,3,-2) - (0,1,0) is 1,2,-2

+(1,3,-2) is 1,3,-2

-(1,3,-2) is -1,-3,2

v is 1,3,-2

v += (1,1,1) is 2,4,-1

v -= (1,1,1) is 1,3,-2

 

(1,3,-2) . (0,1,0) is 3

(1,3,-2) X (0,1,0) is 2,0,1

 

From O to 3,4,0 is 5

 

Angle between 1,0,0 and 2,1,0 is 0.463648

Sine    of the angle is 0.447214

Cosine  of the angle is 0.894427

Tangent of the angle is 0.5

 

The area of triangle with three vertices (0,0,1)(3,0,1)(3,4,1) is 6

The area of triangle with two edges (3,0,0)(0,4,0) is 6

The volume of parallelepiped with three edges (1,0,0)(0,2,0)(0,0,3) is 6

 

Listing 3.1.1.0.1 Other testing and the output result.

 

3.1.2 The shadow classes cpl.mac.vm.vector3DTYP

 

The template class  mac::vector3D<T> is slightly modified and then shadowed into Python.  A shadow class in Python is really a Python class.  It calls various C++ functions generated by SWIG.  See Appendix 1 for examples of SWIG.  The last three letters of the shadow class name TYP are all capitalized, and indicate the type T.  The available TYPs are DBL, FLT, LIN, INT, SIN, LUI, UIN, and SUI (Table 1.3.2.0.1).

 

The member functions obviously can be accessed by the operator . following a variable identifier v, u or w.  Static functions are prefixed by cpl.mac.vm.vector3DTYP_ or cpl.mac.vm.vector3DTYP. as specified below, which depends on whether they are overloaded and renamed while using SWIG.  In general, prefix cpl.mac.vm.vector3DTYP. always works.

 

3.1.2.2 Constructors

vector3DTYP(x = 0, y = 0, z = 0)

The constructor takes no more than three values and constructs a vector3DTYP.  The missing values will be filled by the default 0.

 

vector3DTYP(v)

The copy constructor constructs a new object same as the given one.

 

vector3DTYP(p, s = 1)

Constructs a vector of length s along an orientation specified by polar angles p, where p is a complex (y + jj).

 

No assignment is needed.  Python offers the operator =, but the new object is only a reference to the existing one.

 

3.1.2.3 Element access

v(i)

The operator () takes a 1-offset subscript, and returns the element.  No 0-offset operator is available, since the operator [] does not work in Python.

 

v.getX()

v.getY()

v.getZ()

Explicit access is available through functions getX(), getY(), and getZ().

 

3.1.2.4 Addition and subtraction

u + v

u ¨C v

u += v

u -= v

+v

-v

All the binary operators +, -, +=, -= and the unary operators +, - are exactly working as anticipated.

 

3.1.2.5 Scaling and multiplication

v *= s

v /= s

Operators *= and /= perform isotropic scaling, where s or 1/s is a scale factor.  However, operator *= can only accept integer.

 

v * s

v / s

Return vectors scaled by s and 1/s, respectively.  (v * s does not work right with Python 2.2.1 and earlier.)

 

v.scale(s)

v.scale(u)

Return vectors scaled by a value s isotropically and by a vector u anisotropically.  v remains unmodified.

 

v *= u

v / u

v /= u

Operators *=, / and /= perform anisotropic scaling.  v is anisotropically scaled by u.

 

v * u

v ** u

Operator * is reserved for dot (scalar) multiplication, and operator ** is for cross (vector) multiplication.

 

v.dot(u)

v.cross(u)

Calculate dot and cross products, respectively.

 

vector3DTYP_dot_mul(u, v)

vector3DTYP.dot_mul(u, v)

vector3DTYP_cross_mul(u, v)

vector3DTYP.cross_mul(u, v)

These static functions return dot and cross products of u and v, respectively.

 

u.scalarTripleProduct(v, w)

vector3DTYP.scalarTripleProduct(u, v, w)

Return scalar triple product of u, v, w.

 

v.pre_mul(m)

Returns a vector equal to v pre-multiplied by matrix m.  v remains unmodified.

 

v *= m

v /= m

v is pre-multiplied by m and the inverse matrix of m, respectively.

 

3.1.2.6 Length, distance and projection

abs(v)

The special function __abs__ is defined, therefore abs(v) returns the length of the vector v.

 

vector3DTYP_abs(v)

vector3DTYP.abs(v)

Return the length of vector v.

 

vector3DTYP_abs(x, y, z)

vector3DTYP.abs(x, y, z)

Return the length of vector (x, y, z).

 

Vector length calculated by the functions above is of type of float even for vector3DINT.  However, these functions below return squared length of a vector in type of TYP.

 

v.norm()

vector3DTYP.norm(v)

vector3DTYP_norm_xyz(x, y, z)

vector3DTYP.norm_xyz(x, y, z)

The squared length is calculated by these functions.

 

v.is_zero()

Tests whether the vector is 0-length.

 

v.__nonzero__()

The special function is defined for truth-value testing.  It returns true if the vector is not 0-length.

 

vector3DTYP_shorterThan(u, v)

vector3DTYP.shorterThan(u, v)

u << v

Return true if u is shorter than v in length.  If they are equal in length, return true if u < v is true.  See 3.1.2.10 Comparison for operator <.

 

u.distance(v)

vector3DTYP.distance(u, v)

Return the distance between two vectors u and v.

 

u.projectionOn(v)

u.projectionOf(v)

Return the projection length.  See 3.1.1.13 for difference.

 

vector3DTYP_projection(u, v)

vector3DTYP.projection(u, v)

vector3DTYP_projectionSq(u, v)

vector3DTYP.projectionSq(u, v)

These two static functions return the length and squared length of projection, respectively.  See 3.1.1.13 for detail.

 

x.projectionOn(u, v)

Returns a complex that holds the projection of vector x on vectors u and v.

 

x.projectionOn(u, v, w)

Returns a vector3DDBL that holds the projection of vector x on vectors u, v, and w.

 

3.1.2.7 Orientation and inter-vector angle

v.orientationAngle()

Returns a vector3DDBL that contains the three orientation angles in radian.

 

v.sin()

v.cos()

v.tan()

Return vector3DDBL containing sine, cosine, and tangent values of the orientation angles, respectively.

 

v.polarAngle()

v.polarAngle(kappa)

The function polarAngle() returns a complex y + ji, and polarAngle(kappa) returns a vector3DDBL that contains the polar rotation angles y, j, and kappa.  See 3.1.1.14 for more details.

 

u.parallel(v)

u.antiparallel(v)

u.perpendicular(v)

Three functions parallel(v), antiparallel(v), and perpendicular(v) test the conditions.

 

u | v

The operator | can also be used to test parallelism.

 

u & v

Returns true if two vectors are anti-parallel.

 

u.angle(v, w = cpl.mac.vm.vector3DDBL())

vector3DTYP.angle(u, v,

                  w = cpl.mac.vm.vector3DDBL())

Return the inter-vector angle in radian.  The optional vector w defines the sign of the angle.  See 3.1.1.15 for detail.

 

u ^ v

The operator ^ is also used to calculate the angle.

 

u.sinAngle(v)

u.cosAngle(v)

u.tanAngle(v)

Three functions sinAngle(v), cosAngle(v), and tanAngle(v) return the trigonometric function values of the inter-vector angle.

 

vector3DTYP.sinAngle(u, v, w =

                     cpl.mac.vm.vector3DDBL())

vector3DTYP.cosAngle(u, v)

vector3DTYP.tanAngle(u, v, w =

                     cpl.mac.vm.vector3DDBL())

Calculate the trigonometric function values of the inter-vector angle.  The optional vector w determines the sign of the results.  See 3.1.1.15 for detail.

 

3.1.2.8 Area and volume

u.triangleVertices(v, w)

vector3DTYP.triangleVertices(u, v, w)

Calculates the triangular area formed by vertices u, v, and w.

 

u.triangleEdges(v)

vector3DTYP.triangleEdges(u, v)

Calculates the triangular area formed by two edges u and v.

 

u.parallelepipedVolume(v, w)

vector3DTYP.parallelepipedVolume(u, v, w)

Calculate the volume the parallelepiped.

 

3.1.2.9 Location

u.inside(v, w)

u.outside(v, w)

Test the location of the vector u.  See 3.1.1.18 for definition.

 

3.1.2.10 Comparison

u == v

u != v

Equality and inequality are tested by operator == and !=, respectively.

 

u.ap_eq(v)

The function ap_eq(v) returns true if two vectors are approximately equal.

 

u < v

u <= v

u > v

u >= v

Operators <, <=, >, and >= test the other relationship.

 

vector3DTYP_lessThan(u, v)

vector3DTYP.lessThan(u, v)

Returns true if u < v is true.

 

u.__cmp__(v)

The special function __cmp__ is defined to return ¨C1, 0, or 0, if this vector u is less, equal to, or greater than the given vector v.

 

v.__hash__()

Returns a long integer for hash key.

 

3.1.2.11 Output

v.__str__()

The special function __str__ is defined to print vector3DTYP.

 

The following program tests cpl.mac.vm.vector3DDBL.  It only prints four vector3DDBLs without other message, if everything is OK.

 

#!/usr/bin/python

 

import math

import cpl

 

V = cpl.mac.vm.vector3DDBL

M = cpl.mac.vm.matrix3DDBL

 

O = V()        # Origin, by the default constructor

X = V(1)       # X-axis

Y = V(0, 1)    # Y-axis

Z = V(0, 0, 1) # Z-axis, by the general constructor

z = Z          # Assignment is not done by C++ operator =

assert repr(Z) == repr(z) # They are the same object

 

y = cpl.mac.vm.vector3DDBL(Y) # Copy constructor

assert repr(Y) != repr(y) # They are different objects

 

A = Y + Z # Addition

B = Z - X # Subtraction

C = X + Y

D = X + Y + Z

 

print 'Origin', O # Using the added function __str__

print 'X-axis', X

print 'Y-axis', Y

print 'Z-axis', z

 

assert X(1) == D(1)     # Element access by 1-offset subscript

assert Y(2) == D(2)     # Operator [] is not working in Python

assert Z(3) == D.getZ() # Explicit element access

# X(1) = 0              # SyntaxError: can't assign to function call

 

# assert V(1.2, 2.7, -3.4).nearest() == V(1, 3, -3)

# nearest() does not work, since it returns vector3D<int>.

 

B += X # Inplace addition and subtraction

B -= Z

assert B == O

 

E = cpl.mac.vm.vector3DDBL(D)

E *= 2 # Isotropic scaling

E /= 2 # Isotropic scaling

assert D == E

 

assert D / V(1, 2, 4) == V(1, 0.5, 0.25) # anisotropic scaling

 

assert D  *  V(1, 2, 3)  == 6 # Dot product

assert D.dot(V(1, 2, 3)) == 6 # Dot product

assert cpl.mac.vm.vector3DDBL_dot_mul(D, V(1, 2, 3)) == 6

 

assert D.cross(V(1, 2, 3)) == V(1, -2, 1) # Cross product

assert D ** V(1, 2, 3)     == V(1, -2, 1) # Cross product

assert cpl.mac.vm.vector3DDBL_cross_mul(D, V(1, 2, 3)) == V(1, -2, 1)

 

assert D.pre_mul(M()) == D # Pre-multiplication by matrix

 

assert abs(V(3, 4, 0)) == 5 # Vector length, by special function __abs__

assert cpl.mac.vm.vector3DDBL_abs(V(3, 4, 0))  == 5

assert cpl.mac.vm.vector3DDBL_abs(3, 4, 0) == 5

assert O.is_zero() # Testing 0-length vector

 

assert V(3, 4, 0).norm() == 25 # Squared length

assert cpl.mac.vm.vector3DDBL.norm(V(3, 4, 0)) == 25

assert cpl.mac.vm.vector3DDBL_norm_xyz(3, 4, 0)  == 25

 

assert O << X # O is shorter than X

assert Y << X # If equal in length, Y << X when Y < X

assert cpl.mac.vm.vector3DDBL_shorterThan(Z, Y)

assert cpl.mac.vm.vector3DDBL.shorterThan(Z, Y)

 

assert O.distance(V(3, 4, 0)) == 5 # Distance

assert cpl.mac.vm.vector3DDBL.distance(O, V(3, 4, 0)) == 5

 

assert D.projectionOn(X) == 1 # Projection of D on X

assert X.projectionOf(D) == 1 # Projection of D on X

assert cpl.mac.vm.vector3DDBL_projection(D, X)   == 1 # Projection of D on X

assert cpl.mac.vm.vector3DDBL.projection(D, X)   == 1 # Projection of D on X

assert cpl.mac.vm.vector3DDBL_projectionSq(D, X) == 1 # Squared projection

assert cpl.mac.vm.vector3DDBL.projectionSq(D, X) == 1 # Squared projection

 

assert X.cos() == V(1)           # Orientation cosine

assert X.sin() == V(0, 1, 1)     # Orientation sine

D = D * math.sqrt(2)

assert D.tan() == D              # Orientation tangent

 

D = V(1, 1, 1)

D = D.scale(math.atan(math.sqrt(2)))

assert D.orientationAngle() == D # Orientation angle

assert Y.polarAngle()  == 0      # Polar angle

assert Y.polarAngle(0) == O      # Polar rotation angles psi, phi, kappa

 

assert (D / D) | E               # Testing parallelism

assert (D / D).parallel(E)

assert (D / D).antiparallel(-E)  # Unary operator -

assert X.perpendicular(Y)        # perpendicular

 

assert math.sin(X ^ Y) == 1      # Angle between vectors

assert math.sin(X.angle(Y)) == 1

assert math.sin(cpl.mac.vm.vector3DDBL.angle(X, Y)) == 1

assert X.sinAngle(Y) == 1

assert X.cosAngle(Y) == 0

assert X.tanAngle(X) == 0

assert cpl.mac.vm.vector3DDBL.sinAngle(X, Y) == 1

assert cpl.mac.vm.vector3DDBL.cosAngle(X, Y) == 0

assert cpl.mac.vm.vector3DDBL.tanAngle(X, X) == 0

 

assert X.triangleVertices(Y, Z) == math.sqrt(3) / 2 # Triangle area

assert cpl.mac.vm.vector3DDBL.triangleVertices(X, Y, Z) == math.sqrt(3) / 2

assert X.triangleEdges(Y) == 0.5

assert cpl.mac.vm.vector3DDBL.triangleEdges(X, Y) == 0.5

 

assert X.scalarTripleProduct(Y,  Z) ==  1

assert X.scalarTripleProduct(Y, -Z) == -1

assert X.parallelepipedVolume(Y,  Z) == 1

assert X.parallelepipedVolume(Y, -Z) == 1

assert cpl.mac.vm.vector3DDBL.scalarTripleProduct(X, Y,  Z) ==  1

assert cpl.mac.vm.vector3DDBL.scalarTripleProduct(X, Y, -Z) == -1

assert cpl.mac.vm.vector3DDBL.parallelepipedVolume(X, Y,  Z) ==  1

assert cpl.mac.vm.vector3DDBL.parallelepipedVolume(X, Y, -Z) ==  1

 

assert O.inside(D, -D) # Testing location

assert X.outside(Y, Z)

 

ang = math.acos(1 / math.sqrt(3)) # 54.73561032 degree

assert((X - Y + Z).orientationAngle().ap_eq(V(ang, 3.1415926 - ang, ang)))

   # Approximately equal

assert (X - Y + Z).orientationAngle() != V(ang, 3.1415926 - ang, ang)

   # Not equal

 

assert V(1, 2, 3) < V(2, 6, 3) # Less than

assert V(1, 2, 3) < V(2, 1, 3)

assert V(1, 2, 3) < V(2, 1, 1)

assert cpl.mac.vm.vector3DDBL_lessThan(V(1, 2, 3), V(2, 6, 3))

assert cpl.mac.vm.vector3DDBL.lessThan(V(1, 2, 3), V(2, 1, 3))

assert cpl.mac.vm.vector3DDBL.lessThan(V(1, 2, 3), V(2, 1, 1))

 

assert V(1, 2, 3) <= V(2, 6, 3) # Less than or equal to

assert V(2, 6, 3) >  V(1, 2, 3) # Greater than

assert V(2, 6, 3) >= V(1, 2, 3) # Greater than or equal to

 

Listing 3.1.2.0.1 Testing the shadow class cpl.mac.vm.vector3DDBL.

 

3.1.3 Container of 3D vectors

 

Some of std::vector<mac::vector3D<T> > are instantiated, and can be accessed by using cpl.mac.vm.vectorV3DTYP, where available TYPs are DBL, FLT, and INT.

 

3.2     Matrices mac::matrix3D<T> and cpl.mac.vm.matrix3DTYP

 

3.2.1 The template class mac::matrix3D<T>

 

This template class mac::matrix3D<T> is the primary implementation of 3 ´ 3 matrix.  The original non-template class matrix is obsolete.

 

3.2.1.1 Constructors

matrix3D(const T a11 = 1,

         const T a12 = 0,

         const T a13 = 0,

         const T a21 = 0,

         const T a22 = 1,

         const T a23 = 0,

         const T a31 = 0,

         const T a32 = 0,

         const T a33 = 1);

The 3 ´ 3 matrix modeled by this class contains only numerical elements of a same type.  The general constructor takes no more than 9 numerical values to construct a matrix.

 

Default constructor

The general constructor is also the default constructor when all arguments are missing.  It constructs a unit matrix with 1 on the diagonal and 0 elsewhere.

 

Copy constructor

matrix3D(const matrix3D<T>&);

The copy constructor copies an existing matrix to construct an identical one.

 

3.2.1.2 Destructor

~matrix3D();

There is nothing to do in the destructor.

 

3.2.1.3 Conversion

operator matrix() const;

operator matrix3D<double>() const;

Conversion operator automatically converts this type to mac::matrix3D<double> is available.  The first conversion operator may not be available in all releases, since the class matrix is obsolete.

 

3.2.1.4 Element access

T  operator()(const int, const int) const;

T& operator()(const int, const int);

mac::vector3D<T> operator[](const int) const;

The const and non-const forms of the operator () and only the const form of operator [] are available to provide element access by 1- and 0-offset subscripts, respectively.  The const forms are used to get values.  The non-const form returns an l-value of an element in the matrix.

 

See the example in Listing 3.1.1.4.1 for the usage of these operators of mac::vector3D<T>.  The operator () of this class behaves identically except that it takes exactly two subscripts, such as r(1, 3).  The operator [] takes one subscript for the row, the first subscript, and returns a vector of the row.  Another operator [] of vector3D<T> could then be used to get the value of an element, for example, r[0][2].

 

mac::vector3D<T> row(const int) const;

Returns a row as a vector.  This function takes a 1-offset subscript.

 

3.2.1.5 Assignment

const matrix3D<T> &operator=(const matrix3D<T>&);

The assignment operator = is defined.

 

3.2.1.6 Addition and subtraction

matrix3D<T> operator+ (const matrix3D<T>&) const;

matrix3D<T> operator- (const matrix3D<T>&) const;

matrix3D<T> operator+=(const matrix3D<T>&);

matrix3D<T> operator-=(const matrix3D<T>&);

matrix3D<T> operator+() const;

matrix3D<T> operator-() const;

Binary operators +, -, +=, and -= perform matrix addition and subtraction as anticipated.  Unary operator + keeps and ¨C flips the sign of every element of the matrix.

 

3.2.1.7 Multiplication by a value (scaling)

matrix3D<T> operator*=(const T);

matrix3D<T> operator/=(const T);

matrix3D<T>      operator*(const T) const;

matrix3D<double> operator/(const T) const;

template <class T> mac::matrix3D<T> operator*(

   const T&,

   const mac::matrix3D<T>&);

The binary operators *, /, *=, and /= are overloaded to perform scaling by a value, for example, r * a, r / a, r *= a, and r /= a.  Also, unbounded function is defined to allow a * r to return the scaled matrix.  Among all these operators, only / returns matrix3D<double> to avoid unanticipated result from some integer division, the others keep the matrix type.  Therefore, the operator / may cause precision lose with long double.

 

3.2.1.8 Multiplication

matrix3D<T>      operator*(

   const matrix3D<T>&) const;

matrix3D<double> operator/(

   const matrix3D<T>&) const;

matrix3D<T> operator*=(const matrix3D<T>&);

matrix3D<T> operator/=(const matrix3D<T>&)

The binary operator * also performs matrix multiplication like r * s.  Its counterpart / is defined so that r / s means r is pre-multiplied by inverse of s, that is, ~s * r.  See 3.2.1.12 for inversion operator ~.  The operators *= are defined to perform pre-multiplication so that r *= s means s * r and after return r is modified to store the returned result; r /= s means ~s * r and r is then modified to store the result.

 

All these operators return the same type matrix3D<T>, except that / always returns matrix3D<double>, since some integer division may cause unexpected result during matrix inversion.  However, operator may result in precision lose in case of long double.  Therefore, as a rule of thumb, use operator / with matrix3D<int>, and operator /= with matrix3D<long double>.   Mixing operator / with matrix3D<long double> will cause precision lose, and mixing operator /= with matrix3D<int> may result in unanticipated integer division.

 

3.2.1.9 Pre-multiplication to a vector

mac::vector3D<T> operator*(

   const mac::vector3D<T>&) const;

The binary operator * can connect a matrix with a vector.  r * v returns a vector3D<T> that is the result of v pre-multiplied by matrix r.

 

3.2.1.10 Determinant

T det() const;

T mabs() const;

template <class T> T  abs(

   const mac::matrix3D<T>&);

template <class T> T mabs(

   const mac::matrix3D<T>&);

bool operator!() const;

Two member functions det() and mabs() are defined to calculate the determinant of the matrix.  Unbounded functions are also defined so that abs(r) and mabs(r) returns the determinant of r.  The operator ! is overloaded to return true if its determinant is 0, and false otherwise.

 

T sub(const int, const int) const;

This member function returns the determinant of the sub-matrix about the given 1-offset subscripts.

 

3.2.1.11 Transpose

matrix3D<T> transpose();

matrix3D<T> operator*() const;

A function transpose() and the unary operator * are defined to return the transpose of the matrix.  The difference is whether the original is transposed.  Both r.transpose() and *r return the transpose of matrix r, but the former transposes r itself, and the latter does not.

 

static matrix3D<T> transpose(const matrix3D<T>&);

The static function is defined to return the transpose as well.

 

bool     symmetric() const;

bool antisymmetric() const;

bool    orthogonal() const;

Return Booleans that test the type of matrix.

 

3.2.1.12 Inversion

matrix3D<T> inverse();

matrix3D<double> operator~() const;

static matrix3D<double> inverse(

   const matrix3D<T>&);

A function inverse() and the unary operator ~ are defined to return the inversion of the matrix.  The function also inverses the original matrix, and it does not change the matrix type.  The operator ~ does not inverse the original matrix, but it always returns matrix3D<double>, since some integer division may behave not as always anticipated.  The operator ~ may cause precision lose when used with long double.  Therefore, use the function inverse() with matrix3D<long double>, and use the operator ~ with matrix3D<int> or the static function.

 

3.2.1.13 Comparison

bool operator==(const matrix3D<T>&) const;

bool operator!=(const matrix3D<T>&) const;

bool operator> (const matrix3D<T>&) const;

bool operator< (const matrix3D<T>&) const;

bool operator>=(const matrix3D<T>&) const;

bool operator<=(const matrix3D<T>&) const;

The comparison operators == and < are defined first to check equality and less than relationship, respectively.  Two matrices are equal if and only if all elements from one matrix are equal to their equivalent ones in the other.  A matrix is less than the other, if its determinant is less, or the earlier elements of the matrices determine the relationship, if their determinants are equal.  The sequence of the elements is the same as that in the constructor (3.2.1.1).  Then the other comparison operators !=, >, <=, and >= are defined accordingly.

 

bool operator&=(const matrix3D<T>&) const;

Returns true if the given matrix is approximately equal to this one.  Approximate equality is defined by all rows in the matrices are approximately equal to their equivalent ones.  See 3.1.1.19 for detail.

 

3.2.1.14 Output

friend ostream &operator<<(ostream&,

                           const matrix3D<T>&);

The operator << is used to print the matrix to an ostream.

 

The following example tests this template class.  This program shall print two matrices only, if everything is fine.

 

#include <iostream.h>

#include <assert.h>

#include <mac/vector3D.h>

#include <mac/matrix3D.h>

 

main()

{

   mac::matrix3D<int> id;         // By default constructor

   mac::matrix3D<int> r(1, 2, 3,

                        4, 5, 6,

                        7, 8, 9); // By general constructor

 

   cout << "Identity matrix" << endl;

   cout << id << endl;    // Printing

   cout << "Another matrix" << endl;

   cout << r << endl;

 

   assert (r(1, 1) == r[0][0] and // Using the const forms of operators () and

           r(1, 2) == r[0][1] and // [].

           r(1, 3) == r[0][2] and

           r(2, 1) == r[1][0] and

           r(2, 2) == r[1][1] and

           r(2, 3) == r[1][2] and

           r(3, 1) == r[2][0] and

           r(3, 2) == r[2][1] and

           r(3, 3) == r[2][2]);

 

   mac::matrix3D<int> s = r; // Assignment

 

   for (int i = 1; i < 4; i++)

   for (int j = 1; j < 4; j++) r(i, j) = i + j;

      // Using the non-const form of the operator () for l-value.

 

   assert ((r + s) == (r += s)); // Addition and subtraction

   assert ((r - s) == (r -= s));

   assert (-(-r) == +r);         // Unary operators + and -.

   assert ((r * 2) == (r *= 2)); // Scaling

   assert ((r / 2) == (r /= 2));

   assert ((r * s) == (s *= r)); // Multiplication

   assert ((r / id) == (r /= id));

 

   mac::vector3D<int> v(1, 2, 3);

   assert (id * v == v); // Vector pre-multiplied by a matrix.

 

   assert (r.det()  == 0); // Determinant

   assert (r.mabs() == 0);

   assert (abs(r)   == 0);

   assert (!r); // Determinant is 0.

 

   r = s;

   assert (r == *(s.transpose())); // Transpose

   assert (r != s);

   assert (r == mac::matrix3D<int>::transpose(s));

 

   r = s = mac::matrix3D<int>(2,5,3,5,3,7,8,5,1);

   s.inverse(); // Inversion results in 0-matrix.

   assert (!s);

   s = r;

   assert (mac::matrix3D<double>(r) &= ~(~s)); // Approximately equal.

 

   return 0;

}

 

Listing 3.2.1.0.1 Testing the template class mac::matrix3D<T>.

 

3.2.2 The shadow classes cpl.mac.vm.matrix3DTYP

 

mac::matrix3D<T> is shadowed into Python and called cpl.mac.vm.matrix3DTYP.  The available TYPs are DBL and FLT.  In the following sections, r and s are two matrices.  The header file matrix3D.h is slightly modified.  The source file matrix3D.c++ contains the implementation of the class is too long to be merged into the header.  It is instead included in the headed by using #include at the end of the header file.  See Appendix 1 for detail.

 

3.2.3.1 Constructors

matrix3DTYP(r11 = 1, r12 = 0, r13 = 0,

            r21 = 0, r22 = 1, r23 = 0,

            r31 = 0, r32 = 0, r33 = 1)

The general constructor takes no more than 9 values to construct a matrix.  The default constructor returns a identity matrix.

 

matrix3DTYP(r)

The copy constructor copies one matrix to another new, independent one.

 

3.2.2.2 Element access

r(i, j)

Only the operator () is working in Python.  It takes two 1-offset subscripts, and returns the element value.  One cannot assign value to this operator.

 

r(i)

Returns the i-th row as a vector3DTYP, where i is a 1-offset subscript.

 

3.2.2.3 Addition and subtraction

r + s

r ¨C s

r += s

r -= s

+r

-r

All binary operators +, -, +=, -= and unary operators +, - are working as they should.

 

3.2.2.4 Multiplication by a value (scaling)

r *= a

r /= a

The in-place operators *= and /= are reserved for scaling by a value a.  However, operator *= can only accept an integer.

 

r * a

r / a

Return matrix scaled by a or 1/a.  (r * a does not work right with Python 2.2.1 and earlier.)

 

3.2.2.5 Matrix multiplication

r * s

r / s

The operator * is used for matrix multiplication.  The operator / is defined to perform pre-multiplication by the matrix inversion, that is, r / s returns ~s * r, where ~s is the inversed s.

 

r *= s

r /= s

r is pre-multiplied by s and its inverse matrix, respectively.  r is then replaced by the product.

 

3.2.2.6 Pre-multiplication to a vector

r * v

Returns a vector3DTYP equal to the vector v pre-multiplied by matrix r.

 

3.2.2.7 Determinant

abs(r)

r.det()

Return the determinant of the matrix r.

 

r.sub(i, j)

Returns the determinant of the sub-matrix about i and j.

 

r.is_zero()

Returns true if r has 0-determinant.

 

r.__nonzero__()

Returns true if the determinant of r is not 0.

 

3.2.2.8 Transpose

r.transpose()

Returns the transposed matrix of r, also transposes r in place.

 

matrix3DTYP_trans(r)

matrix3DTYP.trans(r)

Returns the transposed matrix of r.  r remains.

 

r.symmetric()

r.antisymmetric()

r.orthogonal()

Return Booleans that test the type of matrix r.

 

3.2.2.9 Inversion

r.inverse()

Returns the inversion of r.  r will also be inversed after return.

 

matrix3DTYP_inv(r)

matrix3DTYP.inv(r)

~r

Return the inversion of r.  r remains.

 

3.2.2.10 Comparison

r == s

r != s

r < s

r <= s

r > s

r >= s

All comparison operators are available.  See 3.2.1.13 for detail.

 

r.ap_eq(s)

Returns true if the matrices r and s are approximately equal to each other.

 

r.__cmp__(s)

Returns ¨C1, 0, or 0, if r is less, equal to, or greater than s, respectively.

 

r.__hash__()

Returns a long integer as the hash key.

 

3.2.2.11 Output

r.__str__()

The special function __str__ is defined to print a matrix.  It takes three rows to print a matrix, therefore mixing matrix with other messages may be difficult.

 

This example program tests the shadow class.  It prints two matrices only if everything is fine.

 

#!/usr/bin/python

 

import cpl

 

V = cpl.mac.vm.vector3DDBL

M = cpl.mac.vm.matrix3DDBL

 

I = M()         # Identity matrix, by the default constructor

r = M(1, 2, 3,\

      4, 5, 6,\

      7, 8, 9)  # By general constructor

print 'Identity matrix'

print I # Printing

print 'Another matrix'

print r

 

assert r(1, 1) == 1 # Element access

assert r(1, 2) == 2

assert r(1, 3) == 3

assert r(1) == V(1, 2, 3) # A row

 

R = r                         # Assignment

assert repr(R) == repr(r)     # They are the same object.

s = cpl.mac.vm.matrix3DDBL(r) # By copy constructor

assert repr(s) != repr(r)     # They are different objects.

 

s -= I             # In place subtraction

r += s             # In place addition

t = r - s          # Subtraction

t *= 2             # Scaling

assert t == r + I  # Addition

 

t = r * s          # Multiplication

t = t / r          # t pre-multiplied by inverse of r

assert t.ap_eq(s)  # Approximately equal

 

assert I * V(1, 2, 3) == V(1, 2, 3) # Vector pre-multiplied by matrix

 

assert abs(r) == r.det() # Two ways of calculating determinant

assert not r.is_zero()   # 0-determinant

 

t = r.transpose() # Both t and r after return are transpose of the original r.

assert t == r

t = cpl.mac.vm.matrix3DDBL_trans(r) # Only t is transpose of r.

t = cpl.mac.vm.matrix3DDBL.trans(r) # Only t is transpose of r.

assert t.transpose() == r           # t is transposed again.

 

assert I.symmetric()

assert not I.antisymmetric()

assert I.orthogonal()

 

t = r.inverse() # Both t and r after return are inverse of the original r.

assert t == r

t = cpl.mac.vm.matrix3DDBL_inv(r)   # Only t is inverse of r.

t = cpl.mac.vm.matrix3DDBL.inv(r)   # Only t is inverse of r.

assert r.ap_eq(~t)                  # Operator ~ returns inverse of t, but

                                    # t is not inversed again.

assert r.ap_eq(~~r)

 

Listing 3.2.2.0.1 Testing the shadow class cpl.mac.vm.matrix3DDBL.

 

3.2.3 Container of 3D Matrices

 

Some of std::vector<mac::matrix3D<T> > are instantiated, and can be accessed by using cpl.mac.vm.vectorM3DTYP.  The available TYPs are DBL, FLT, and INT.

 

3.3     Rotation Matrices

 

A rotation matrix is an orthogonal 3-dimentional matrix, that is, its inverse is its transpose.  A rotation matrix represents a pure rotation about the origin in 3-dimentional space.  A vector pre-multiplied by the rotation matrix results in another vector that is rotated about the origin from the original one.  Obviously, rotation matrices are very useful in crystallography.  There are a variety of ways to define a rotational operation, therefore rotation matrices should also be constructed by different ways to fit the needs.  All these rotation matrices discussed below are derived from mac::matrix3D<double>.

 

3.3.1 Polar Rotation Matrices mac::polarRotationMatrix and cpl.mac.vm.polarRotationMatrix

 

Polar angle defining an orientation in 3-dimensional space is described in 3.1.1.14.  Polar rotation refers to a rotational operation about the axis defined by a pair of polar angles y and j.  The rotation angle is k, whose positive direction is defined by the right-hand rule about its rotation axis.  Therefore, three angular values y, j and k can define the polar rotation matrix.

 

3.3.1.1 mac::polarRotationMatrix

Constructors

polarRotationMatrix(const double psi   = 0,

                    const double phi   = 0,

                    const double kappa = 0);

polarRotationMatrix(const complex<double> polar,

                    const double kappa);

polarRotationMatrix(

   const vector3D<double> angle);

polarRotationMatrix(const vector3D<double> axis,

                    const double kappa);

polarRotationMatrix(const matrix3D<double>&);

polarRotationMatrix(const int axis,

                    const double kappa);

The default constructor polarRotationMatrix() constructs an identity matrix that defines 0 rotation.  The first constructor, obviously, takes three doubles for y, j and k in radian to construct a polar rotation matrix.  The polar angles y and j can be merged into a complex<double> in the second constructor, and all three angles can be merged into a vector3D<double> in the third constructor.  The fourth constructor takes a vector3D<double> that defines an axis and a double as a rotation angle.  The fifth is like a copy constructor, which will allow, for example, conversion from a multiplication of two rotation matrices back to the same type of rotation matrix.  This is necessary since the result of multiplication has the type of matrix3D<double>.  The sixth is a special constructor that constructs rotation matrices about three principal axes X, Y, and Z only.  The first int can be a single letter ¡®x¡¯, ¡®y¡¯, ¡®z¡¯, ¡®X¡¯, ¡®Y¡¯, or ¡®Z¡¯.

 

void set(const double,

         const double,

         const double);

If the matrix already exists, this public function can be used to set the matrix by y, j and k.

 

void set(const double,

         const double,

         const double,

         const double,

         const double,

         const double,

         const double,

         const double,

         const double);

Another overloaded function set() takes all 9 elements of the matrix.  A matrix set by this function must be orthogonal so that it is still a rotation matrix.  Non-orthogonal matrix may result in error when back calculating the rotation angles (see below).

 

Destructor

~polarRotationMatrix();

There is nothing to do in destructor.

 

Back calculate the polar rotation angles from the matrix

If a matrix is the result of the multiplication of two rotation matrices, how large is the total rotation, and about which orientation?  One needs to back calculate the rotation angles from the matrix.

 

double kappa() const;

This function returns the rotation angle k in radian.

 

vector3D<double> polar() const;

This function returns y, j and k in radian as a vector.  The back calculation is non-trivial and out of the scope of this manual, I will leave it behind the scene.

 

These functions for back calculation of polar angles may not return the same values as the original ones used to construct the matrix.  The returned polar angles are only the equivalent ones.

 

Partial derivatives to the polar angles

static mac::matrix3D<double> derivative(

   const double psi,

   const double phi,

   const double kappa,

   const int ang);

static mac::matrix3D<double> derivative(

   const mac::vector3D<double>& v,

   const int ang);

Return the partial direvative of a polar rotation matrix at the polar angles psi, phi, and kappa to the ang-th polar angle, where ang is a 1-offset subscript from 1 to 3.  The second function collects the polar angles in a vector v.  The returned derivative is obviously a matrix, but not necessarily a rotation matrix.

 

static mac::vector3D<double> derivative(

   const double psi,

   const double phi,

   const double kappa,

   const int i, const int j);

static mac::vector3D<double> derivative(

   const mac::vector3D<double>&,

   const int, const int);

Let a vector (x1, x2, x3) = Ru, where R is the polar rotation matrix at psi, phi, kappa.  These static functions return a vector v so that xi/pj = v ¡¤ u, where p1 = psi, p2 = phi, p3 = kappa.  The integers i and j are 1-offset subscripts.

 

Testing

The following example tests this class.  See comments therein for more details.

 

#include <assert.h>

#include <mac/vector3D.h>

#include <mac/rotationMatrix.h>

 

main()

{

   typedef mac::polarRotationMatrix P;

 

   assert (P().kappa() == 0);

      // The default constructor constructs 0 rotation.

   assert (P(0, 0, 1) == P('y', 1));

      // Polar angle (0, 0) is Y-axis.

   assert (P(M_PI / 2 * complex<double>(1 , -1), 1) == P('z', 1)); // Z-axis.

   assert (P(mac::vector3D<double>(M_PI / 2, 0, 1)) ==

           P(mac::vector3D<double>(1), 1)); // X-axis.

   P p, q;

   p.set(0, 0, M_PI);

   q.set(-1, 0,  0,

          0, 1,  0,

          0, 0, -1);

   assert (p &= q);

   assert (p.polar() == mac::vector3D<double>(0, 0, M_PI));

 

   p.set(M_PI / 2, 0, M_PI);

   assert (P(p * q).polar() == M_PI * mac::vector3D<double>(0.5, 0.5, 1));

      // The copy constructor allows the product call member function polar().

   assert (P('y', M_PI) * P('x', M_PI) &= P('z', M_PI));

 

   return 0;

}

 

Listing 3.3.1.1.1 Testing the class mac::polarRotationMatrix.

 

3.3.1.2 cpl.mac.vm.polarRotationMatrix

Polar rotation matrix is shadowed into a Python class that inherits from its base class matrix3DDBL.  r is a polar rotation matrix in this section below.

 

polarRotationMatrix(psi = 0, phi = 0, kappa = 0)

The general constructor takes no more than three angular values in radian.  By default, it constructs a identity matrix.

 

polarRotationMatrix(polar, kappa = 0)

Constructs a polarRotationMatrix from a complex of polar orientation polar and a polar angle kappa.

 

polarRotationMatrix(angle)

Constructs a polarRotationMatrix from a vector3DDBL that contains polar angles y, j and k.

 

polarRotationMatrix(axis, kappa)

Constructs a polarRotationMatrix that represents a rotation of angle kappa around an axis specified by a vector3DDBL axis.

 

polarRotationMatrix(r)

The copy constructor converts matrix3DDBL to polarRotationMatrix, which would allow the result from an inherited function have the member functions of polarRotationMatrix.  Since a polarRotationMatrix is also a matrix3DDBL, this constructor truly functions as a copy constructor.

 

r.set(psi, phi, kappa)

The set function sets the polar angles again.

 

r.set(r11, r12, r13,

      r21, r22, r23,

      r31, r32, r33)

Sets the matrix by elements.  Use this function cautiously, since there is an internal assertion that the newly set matrix must be orthogonal.

 

r.kappa()

r.polar()

Return kappa and polar angles, respectively.

 

derivative(psi, phi, kappa, ang)

derivative(polar, ang)

Return the partial derivative of a polar rotation matrix at psi, phi, kappa to the ang-th polar angle.  The second function collects the polar angles in a vector3DDBL polar.

 

derivative(psi, phi, kappa, i, j)

derivative(polar, i, j)

Let a vector (x1, x2, x3) = Ru, where R is the polar rotation matrix at psi, phi, kappa.  These functions return a vector v so that xi/pj = v ¡¤ u, where p1 = psi, p2 = phi, p3 = kappa.  The integers i and j are 1-offset subscripts.  The second function collects the polar angles in a vector3DDBL polar.

 

#!/usr/bin/python

 

import math

import cpl

 

V = cpl.mac.vm.vector3DDBL

P = cpl.mac.vm.polarRotationMatrix

 

p = P()                     # The default constructor.

assert p.kappa() == 0

p.set(0, 0, math.pi)        # Change the matrix by set.

q = P(math.pi, 0, -math.pi) # The general constructor.

assert p.ap_eq(q)           # Using inherited functions.

assert p.polar() == V(0, 0, math.pi)

p.set(math.pi / 2, 0, math.pi)

assert cpl.mac.vm.polarRotationMatrix(p * q).polar() ==\

       V(math.pi / 2, math.pi / 2, math.pi)

   # The copy constructor allows the product call member function polar().

assert (p * q).orthogonal() # Using inherited functions.

 

Listing 3.3.1.2.1 Testing the shadow class cpl.mac.vm.polarRotationMatrix.

 

3.3.1.3 Container of polar rotation matrices

std::vector<mac::polarRotationMatrix> is instantiated, and can be accessed by using cpl.mac.vm.vectorPRM.

 

3.3.2 Polar-to-polar Rotation Matrices mac::polarToPolarRotationMatrix and cpl.mac.vm.polarToPolarRotationMatrix

 

Sometimes it may be desirable to rotate one orientation to another given orientation.  Both of these orientations can be defined as polar angles.  A polar-to-polar rotation matrix will do exactly this, but the rotation about these orientations is arbitrary.

 

3.3.2.1 mac::polarToPolarRotationMatrix

polarToPolarRotationMatrix(

   const complex<double> = 0,

   const complex<double> = 0);

polarToPolarRotationMatrix(

   const mac::vector3D<double>,

   const mac::vector3D<double>);

polarToPolarRotationMatrix(

   const complex<double>,

   const mac::vector3D<double>);

polarToPolarRotationMatrix(

   const mac::vector3D<double>,

   const complex<double>);

polarToPolarRotationMatrix(

   const mac::matrix3D<double>&);

The default constructor polarToPolarRotationMatrix() constructs an identity matrix that defines 0 rotation.  The general constructors can take two complex<double>s, two vector3D<double>s, or mix-and-match of these that define two orientations.  A copy constructor converts matrix3D<double> to polarToPolarRotationMatrix.

 

~polarToPolarRotationMatrix();

There is nothing to do in destructor, and no other public function.

 

Since the rotation about the two given orientations is not guaranteed, the total rotation carried out by this matrix may be greater than the angle between the starting and ending orientations.

 

#include <assert.h>

#include <mac/vector3D.h>

#include <mac/rotationMatrix.h>

 

main()

{

   typedef mac::polarToPolarRotationMatrix P;

   typedef mac::vector3D<double>           V;

 

   P o;

   assert (mac::polarRotationMatrix(o).kappa() == 0);

      // The default constructor constructs 0 rotation.

   const complex<double> x(M_PI / 2);

   const complex<double> y;

   const complex<double> z(M_PI / 2, -M_PI / 2);

   P p(x, y); // Constructed by two complexes.

   assert (p * V(1) &= V(0, 1));

 

   for (double u = -3; u < 3; u += 0.3)

   for (double v = -3; v < 3; v += 0.3)

   for (double w = -3; w < 3; w += 0.3)

   {

      V t(u, v, w);

      P q(t, x); // Constructed by a vector and a complex.

      assert (q * t &= abs(t) * V(1));

      assert (mac::polarRotationMatrix(q).kappa() - (t ^ V(1)) > -FLT_EPSILON);

         // Polar-to-polar rotation is no less than the angle between the two

         // polar orientations.

   }

 

   return 0;

}

 

Listing 3.3.2.1.1 Testing polar-to-polar rotation matrix.

 

3.3.2.2 cpl.mac.vm.polarToPolarRotationMatrix

polarToPolarRotationMatrix(p1 = 0, p2 = 0)

Constructs a polarToPolarRotationMatrix that rotates from the first polar orientation p1 to the second p2.  p1 and p2 can either be complex or vector3DDBL or mix-and-match.

 

polarToPolarRotationMatrix(r)

Copies a matrix3DDBL to a polarToPolarRotationMatrix.

 

#!/usr/bin/python

 

import math

import cpl

 

V = cpl.mac.vm.vector3DDBL

P = cpl.mac.vm.polarToPolarRotationMatrix

 

p = P() # The default constructor.

assert cpl.mac.vm.polarRotationMatrix(p).kappa() == 0

 

x = math.pi / 2

y = 0

z = math.pi / 2 * (1 - 1j)

p = P(x, y)

assert (p * V(1)).ap_eq(V(0, 1))

 

t = V(3.4, 0.46, -2.5)

s = V(0, 0, 1)

s = s * abs(t)

q = P(t.polarAngle(), z)

assert (q * t).ap_eq(s)

assert cpl.mac.vm.polarRotationMatrix(q).kappa() >= t.angle(V(1))

   # Polar-to-polar rotation is no less than the angle between the two

   # polar orientations.

 

Listing 3.3.2.2.1 Testing polar-to-polar rotation matrix as the shadow class in Python.

 

3.3.2.3 Container of polar-to-polar rotation matrices

std::vector<mac::polarToPolarRotationMatrix> is instantiated, and can be accessed by using cpl.mac.vm.vectorPPM.

 

3.3.3 Euler Rotation Matrices mac::EulerRotationMatrix and cpl.mac.vm.EulerRotationMatrix

 

Another commonly-used rotation matrix is Euler rotation matrix.  The rotation operation described by this matrix is in fact composed from three rotation operations about three principal axes, therefore, this matrix can be considered as the product of three polar rotation matrices.  Since these rotations are about the principal axes, this matrix is easier to compute.  However, it is not directly obvious how large is the total rotation angle if this matrix is not converted to polar rotation matrix.  If all three rotations are small, the total rotation is obviously small, but it is entirely possible that some large rotations may compose a small total rotation.  This behavior makes it a bad choice to refine in a least-squares minimization.  It can be shown that the total rotation angle is greater than or equal to the smallest absolute Euler angle, and less than or equal to the sum of all absolute Euler angles.

 

min(|q1|, |q2|, |q3|) £ k £ |q1| + |q2| + |q3|.

 

3.3.3.1 mac::EulerRotationMatrix

EulerRotationMatrix(const double t1 = 0,

                    const double t2 = 0,

                    const double t3 = 0);

EulerRotationMatrix(const mac::vector3D<double>);

EulerRotationMatrix(

   const mac::matrix3D<double>&);

The default constructor EulerRotationMatrix() constructs an identity matrix that defines 0 rotation.  The general constructors can take either three angles in radian or a vector3D<double> that consists of three angles.  A copy constructor converts a matrix3D<double> to an EulerRotationMatrix.

 

void set(const mac::vector3D<double>);

A function set() takes a vector3D<double> of three angles to set the matrix to another value.

 

mac::vector3D<double> Euler() const;

A public function Euler() back-calculates the Euler angles and returns a vector3D<double>.  The returned Euler angles are only a set of equivalent values, may not be the same as the original one used to construct the matrix.

 

Partial derivatives to the Euler angles

static mac::matrix3D<double> derivative(

   const double t1,

   const double t2,

   const double t3,

   const int i);

static mac::matrix3D<double> derivative(

   const mac::vector3D<double>&, const int);

Return the partial derivative of the Euler rotation matrix at Euler angles t1, t2, t3 to one of the Euler angles specified by i.  The second function collects the Euler angles in a vector.

 

static mac::vector3D<double> derivative(

   const double t1,

   const double t2,

   const double t3,

   const int i, const int j);

static mac::vector3D<double> derivative(

   const mac::vector3D<double>&,

   const int, const int);

Let a vector (x1, x2, x3) = Ru, where R is the Euler rotation matrix at t1, t2, t3.  These static functions return a vector v so that xi/tj = v ¡¤ u.  The integers i and j are 1-offset subscripts.

 

#include <assert.h>

#include <mac/vector3D.h>

#include <mac/rotationMatrix.h>

 

main()

{

   typedef mac::EulerRotationMatrix E;

   typedef mac::polarRotationMatrix P;

   typedef mac::vector3D<double>    V;

 

   E e;

   assert (P(e).kappa() == 0);

      // The default constructor constructs 0 rotation.

   assert (E(0.1, 0, 0.2) &= E(0.3, 0, 0));

   assert (E(0.1, 0.5, 0.25).Euler() &= V(0.1, 0.5, 0.25));

 

   for (double t1 = -M_PI; t1 < M_PI; t1 += 0.1)

   for (double t2 = -M_PI; t2 < M_PI; t2 += 0.1)

   for (double t3 = -M_PI; t3 < M_PI; t3 += 0.1)

   {

      E m(t1, t2, t3);

      const double kappa = P(m).kappa();

      assert (kappa <= fabs(t1) + fabs(t2) + fabs(t3));

      assert (kappa > min(min(fabs(t1), fabs(t2)), fabs(t3)) - FLT_EPSILON);

   }

 

   return 0;

}

 

Listing 3.3.3.1.1 Testing Euler rotation matrix.

 

3.3.3.2 cpl.mac.vm.EulerRotationMatrix

EulerRotationMatrix(t1 = 0, t2 = 0, t3 = 0)

EulerRotationMatrix(t)

Construct an EulerRotationMatrix from three Euler angles.  The second constructor collects the Euler angles to a vector3DDBL.

 

EulerRotationMatrix(r)

Copies a matrix3DDBL r to an EulerRotationMatrix.

 

r.set(t)

Sets the rotation matrix r by a new set of Euler angles in a vector3DDBL t.

 

r.Euler()

The function Euler() returns the Euler angles of the rotation matrix.

 

derivative(t1, t2, t3, i)

derivative(t, i)

Return the partial derivative of a Euler rotation matrix at t1, t2, t3 to the i-th Euler angle.  The second function collects the Euler angles to a vector3DDBL t.

 

derivative(t1, t2, t3, i, j)

derivative(t, i, j)

Let a vector (x1, x2, x3) = Ru, where R is the Euler rotation matrix at t1, t2, t3.  This function returns a vector v so that xi/tj = v ¡¤ u.  The integers i and j are 1-offset subscripts.  The second function collects the Euler angles to a vector3DDBL t.

 

#!/usr/bin/python

 

import math

import cpl

 

V = cpl.mac.vm.vector3DDBL

E = cpl.mac.vm.EulerRotationMatrix

P = cpl.mac.vm.polarRotationMatrix

 

e = E() # The default constructor.

assert P(e).kappa() == 0

 

assert E(0.1, 0, 0.2).ap_eq(E(0.3, 0, 0))

assert E(0.1, 0.5, 0.25).Euler().ap_eq(V(0.1, 0.5, 0.25))

 

for i in range(30):

   t1 = -math.pi + i * 0.2

   for j in range(30):

      t2 = -math.pi + j * 0.2

      for k in range(30):

         t3 = -math.pi + k * 0.2

         m = E(t1, t2, t3)

         kappa = P(m).kappa()

         assert kappa <= abs(t1) + abs(t2) + abs(t3)

         assert kappa > min(min(abs(t1), abs(t2)), abs(t3)) - 1e-6

 

Listing 3.3.3.2.1 Testing Euler rotation matrix in Python.

 

3.3.3.3 Container of Euler rotation matrices

std::vector<mac::EulerRotationMatrix> is instantiated, and can be accessed by using cpl.mac.vm.vectorERM.