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()

{