No Subject

[dwhitney at localhost dot localdomain]

#ifndef _LINALGEBRA_H_
#define _LINALGEBRA_H_

/*
 * LinAlgebra.h
 */

#include <math.h>
#include <string.h>


/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

template <class T1,class T2>
void CrossProduct3(const T1 *v1,const T2 *v2,T1 *r);

template <class T1,class T2>
inline T1 DotProduct3(const T1 *v1,const T2 *v2) {
    return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
}

template <class T1,class T2>
inline T1 DotProduct4(const T1 *v1,const T2 *v2) {
    return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]+v1[3]*v2[3];
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

template<class T> class Vector4t;

template<class T>
class Vector3t {
public:
    T vec[3];

    Vector3t() { }
    Vector3t(const T& t1,const T& t2,const T& t3) {
        vec[0]=t1; vec[1]=t2; vec[2]=t3;
    }
    void set(const T& t1,const T& t2, const T& t3) {
        vec[0]=t1; vec[1]=t2; vec[2]=t3;
    }
    Vector3t(const T *tp) {
        vec[0]=tp[0]; vec[1]=tp[1]; vec[2]=tp[2];
    }
    Vector3t(const Vector3t& v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2];
    }
    Vector3t(const Vector4t<T>& v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2];
    }
    Vector3t& operator=(const Vector3t& v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2];
        return *this;
    }
    Vector3t& operator=(const T *v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2];
        return *this;
    }
    Vector3t& operator=(const Vector4t<T>& v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2];
        return *this;
    }

    friend int operator!=(const Vector3t& v1,const Vector3t& v2) {
        return (v1[0]!=v2[0] || v1[1]!=v2[1] || v1[2]!=v2[2]);
    }
    friend int operator==(const Vector3t& v1,const Vector3t& v2) {
        return (v1[0]==v2[0] && v1[1]==v2[1] && v1[2]==v2[2]);
    }

    Vector3t<T>& operator+=(const Vector3t<T>& v) {
        vec[0] += v.vec[0]; vec[1] += v.vec[1]; vec[2] += v.vec[2];
        return *this;
    }
    Vector3t<T>& operator-=(const Vector3t<T>& v) {
        vec[0] -= v.vec[0]; vec[1] -= v.vec[1]; vec[2] -= v.vec[2];
        return *this;
    }
    Vector3t<T>& operator*=(const T& f) {
        vec[0] *= f;  vec[1] *= f; vec[2] *= f;
        return *this;
    }
    Vector3t<T>& operator/=(const T& f) {
        vec[0] /= f;  vec[1] /= f; vec[2] /= f;
        return *this;
    }

    friend Vector3t<T> operator-(const Vector3t<T>& v) {
        return Vector3t<T>(-v.vec[0], -v.vec[1], -v.vec[2]);
    }
    friend Vector3t<T> operator+(const Vector3t<T>& lhs,
                                 const Vector3t<T>& rhs) {
        return Vector3t<T>(lhs.vec[0] + rhs.vec[0],
                           lhs.vec[1] + rhs.vec[1],
                           lhs.vec[2] + rhs.vec[2]);
    }

    friend Vector3t<T> operator-(const Vector3t<T>& lhs,
                                 const Vector3t<T>& rhs) {
        return Vector3t<T>(lhs.vec[0] - rhs.vec[0],
                           lhs.vec[1] - rhs.vec[1],
                           lhs.vec[2] - rhs.vec[2]);
    }
    friend Vector3t<T> operator*(float lhs, const Vector3t<T>& rhs) {
        return Vector3t<T>(lhs * rhs.vec[0], lhs * rhs.vec[1], lhs * rhs.vec[2]);
    }
    friend Vector3t<T> operator*(const Vector3t<T>& lhs, float rhs) {
        return Vector3t<T>(lhs.vec[0] * rhs, lhs.vec[1] * rhs, lhs.vec[2] * rhs);
    }
    friend Vector3t<T> operator/(const Vector3t<T>& lhs, float rhs) {
        return Vector3t<T>(lhs.vec[0] / rhs, lhs.vec[1] / rhs, lhs.vec[2] / rhs);
    }
    friend Vector3t<T> operator*(double lhs, const Vector3t<T>& rhs) {
        return Vector3t<T>(lhs * rhs.vec[0], lhs * rhs.vec[1], lhs * rhs.vec[2]);
    }
    friend Vector3t<T> operator*(const Vector3t<T>& lhs, double rhs) {
        return Vector3t<T>(lhs.vec[0] * rhs, lhs.vec[1] * rhs, lhs.vec[2] * rhs);
    }
    friend Vector3t<T> operator/(const Vector3t<T>& lhs, double rhs) {
        return Vector3t<T>(lhs.vec[0] / rhs, lhs.vec[1] / rhs, lhs.vec[2] / rhs);
    }

    operator T *() { return (T *)vec; }
    operator const T *() const { return (const T *)vec; }

    T operator [](int i) const { return vec[i]; }
    T &operator [](int i) { return vec[i]; }

    T length() const {
        return (T)sqrt((double)(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]));
    }

    void scale(const T &t) { vec[0]*=t; vec[1]*=t; vec[2]*=t; }

    void normalize() { scale(1.0/length()); }
    Vector3t<T> normalized() const { return ((*this) * (1.0/length())); }

    T dot(const Vector3t<T> &t) const {
        return DotProduct3(vec,t);
    }
};


template <class T1,class T2>
Vector3t<T1> operator*(const Vector3t<T1>& v1, 
                       const Vector3t<T2>& v2)
{
    Vector3t<T1> r;
    CrossProduct3(v1.vec,v2.vec,r.vec);
    return r;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

template<class T>
class Vector4t {
public:
    T vec[4];

    Vector4t() { }
    Vector4t(const T& t1,const T& t2, const T& t3,const T& t4) {
        vec[0]=t1; vec[1]=t2; vec[2]=t3; vec[3]=t4;
    }
    void set(const T& t1,const T& t2, const T& t3,const T& t4) {
        vec[0]=t1; vec[1]=t2; vec[2]=t3; vec[3]=t4;
    }
    Vector4t(const T *tp) {
        vec[0]=tp[0]; vec[1]=tp[1]; vec[2]=tp[2]; vec[3]=tp[3];
    }
    Vector4t(const Vector4t& v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2]; vec[3]=v[3];
    }
    Vector4t(const Vector3t<T>& v,int i3=0) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2];
        vec[3] = i3;
    }
    Vector4t& operator=(const Vector4t& v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2]; vec[3]=v[3];
        return *this;
    }
    Vector4t& operator=(const Vector3t<T>& v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2]; vec[3]=v[3];
        vec[3] = 0;
        return *this;
    }
    Vector4t& operator=(const T *v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2]; vec[3]=v[3];
        return *this;
    }
    void set(const Vector3t<T>& v,int i3) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2];
        vec[3] = i3;
    }
    void load3(const Vector3t<T>& v) {
        vec[0]=v[0]; vec[1]=v[1]; vec[2]=v[2];
    }

    friend int operator!=(const Vector4t& v1,const Vector4t& v2) {
        return (v1[0]!=v2[0] || v1[1]!=v2[1] || v1[2]!=v2[2] ||
                v1[3]!=v2[3]);
    }
    friend int operator==(const Vector4t& v1,const Vector4t& v2) {
        return (v1[0]==v2[0] && v1[1]==v2[1] && v1[2]==v2[2] ||
                v1[3]!=v2[3]);
    }

    Vector4t<T>& operator-=(const Vector4t<T>& v) {
        vec[0] -= v.vec[0]; vec[1] -= v.vec[1]; vec[2] -= v.vec[2]; vec[3] -= v.vec[3];
        return *this;
    }
    Vector4t<T>& operator+=(const Vector4t<T>& v) {
        vec[0] += v.vec[0]; vec[1] += v.vec[1]; vec[2] += v.vec[2]; vec[3] += v.vec[3];
        return *this;
    }
    Vector4t<T>& operator*=(const T& f) {
        vec[0] *= f;  vec[1] *= f; vec[2] *= f; vec[3] *= f;
        return *this;
    }
    Vector4t<T>& operator/=(const T& f) {
        vec[0] /= f;  vec[1] /= f; vec[2] /= f; vec[3] /= f;
        return *this;
    }
    
    friend Vector4t<T> operator-(const Vector4t<T>& v) {
        return Vector4t<T>(-v.vec[0], -v.vec[1], -v.vec[2], -v.vec[3]);
    }
    friend Vector4t<T> operator+(const Vector4t<T>& lhs,
                                 const Vector4t<T>& rhs) {
        return Vector4t<T>(lhs.vec[0] + rhs.vec[0],
                           lhs.vec[1] + rhs.vec[1],
                           lhs.vec[2] + rhs.vec[2],
                           lhs.vec[3] + rhs.vec[3]);
    }
    friend Vector4t<T> operator-(const Vector4t<T>& lhs,
                                 const Vector4t<T>& rhs) {
        return Vector4t<T>(lhs.vec[0] - rhs.vec[0],
                           lhs.vec[1] - rhs.vec[1],
                           lhs.vec[2] - rhs.vec[2],
                           lhs.vec[3] - rhs.vec[3]);
    }

    friend Vector4t<T> operator*(float lhs, const Vector4t<T>& rhs) {
        return Vector4t<T>(lhs * rhs.vec[0], lhs * rhs.vec[1],
                           lhs * rhs.vec[2], lhs * rhs.vec[3]);
    }
    friend Vector4t<T> operator*(const Vector4t<T>& lhs, float rhs) {
        return Vector4t<T>(lhs.vec[0] * rhs, lhs.vec[1] * rhs,
                           lhs.vec[2] * rhs, lhs.vec[3] * rhs);
    }
    friend Vector4t<T> operator/(const Vector4t<T>& lhs, float rhs) {
        return Vector4t<T>(lhs.vec[0] / rhs, lhs.vec[1] / rhs,
                           lhs.vec[2] / rhs, lhs.vec[3] / rhs);
    }
    friend Vector4t<T> operator*(double lhs, const Vector4t<T>& rhs) {
        return Vector4t<T>(lhs * rhs.vec[0], lhs * rhs.vec[1],
                           lhs * rhs.vec[2], lhs * rhs.vec[3]);
    }
    friend Vector4t<T> operator*(const Vector4t<T>& lhs, double rhs) {
        return Vector4t<T>(lhs.vec[0] * rhs, lhs.vec[1] * rhs,
                           lhs.vec[2] * rhs, lhs.vec[3] * rhs);
    }
    friend Vector4t<T> operator/(const Vector4t<T>& lhs, double rhs) {
        return Vector4t<T>(lhs.vec[0] / rhs, lhs.vec[1] / rhs,
                           lhs.vec[2] / rhs, lhs.vec[3] / rhs);
    }

    operator T *() { return (T *)vec; }
    operator const T *() const { return (const T *)vec; }

    operator Vector3t<T>&() { return *(Vector3t<T> *)&vec; }

    T operator [](int i) const { return vec[i]; }
    T &operator [](int i) { return vec[i]; }

    T length() const { 
        return (T)sqrt((double)(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2] + vec[3]*vec[3]));
    }

    void scale(const T &t) { vec[0]*=t; vec[1]*=t; vec[2]*=t; vec[3]*=t; }

    void normalize() { scale(1.0/length()); }
    Vector4t<T> normalized() const { return ((*this) * (1.0/length())); }

    T dot(const Vector4t<T> &t) const {
        return DotProduct4(vec,t);
    }
    T dot3(const Vector4t<T> &t) const {
        return DotProduct3(vec,t);
    }

    void cross3(const Vector4t &v1,const Vector4t &v2) {
        CrossProduct3((const T *)v1,(const T *)v2,vec);
        vec[3] = 0;
    }
};


/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

template<class T>
class Matrix33t {
    static const int MATSIZE;
public:
    T mat[3][3];

    static const T identity[9];

    Matrix33t() {}
    Matrix33t(const T *tp) {         // ctor takes array of 9 T's
        memcpy(mat,tp,MATSIZE);
    }
    Matrix33t(const T& t1,const T& t2,const T& t3,
              const T& t4,const T& t5,const T& t6,
              const T& t7,const T& t8,const T& t9) {
        mat[0][0]=t1; mat[0][1]=t2; mat[0][2]=t3;
        mat[1][0]=t4; mat[1][1]=t5; mat[1][2]=t6;
        mat[2][0]=t7; mat[2][1]=t8; mat[2][2]=t9;
    }
    void set(const T& t1,const T& t2,const T& t3,
             const T& t4,const T& t5,const T& t6,
             const T& t7,const T& t8,const T& t9) {
        mat[0][0]=t1; mat[0][1]=t2; mat[0][2]=t3;
        mat[1][0]=t4; mat[1][1]=t5; mat[1][2]=t6;
        mat[2][0]=t7; mat[2][1]=t8; mat[2][2]=t9;
    }
    Matrix33t(const Vector3t<T>& v1,const Vector3t<T>& v2,const Vector3t<T>& v3) {
        mat[0][0]=v1[0]; mat[0][1]=v1[1]; mat[0][2]=v1[2];
        mat[1][0]=v2[0]; mat[1][1]=v2[1]; mat[1][2]=v2[2];
        mat[2][0]=v3[0]; mat[2][1]=v3[1]; mat[2][2]=v3[2];
    }
    Matrix33t(const Matrix33t& m) {  // copy ctor
        memcpy(mat,m.mat,MATSIZE);
    }

    Matrix33t<T>& operator=(const Matrix33t<T>& m) {
        memcpy(mat,m.mat,MATSIZE);
        return *this;
    }
    Matrix33t<T>& operator=(const T m[3][3]) {
        memcpy(mat,m,MATSIZE);
        return *this;
    }

    operator T *() const { return (T *)mat; }

    const Vector3t<T>& operator [](int i) const {
        return *(Vector3t<T> *)mat[i];
    }
    Vector3t<T>& operator [](int i) {
        return *(Vector3t<T> *)mat[i];
    }

    void setIdentity() { *this = identity; }

    Matrix33t<T> inverse() const {
        Matrix33t<T> m;
        Invert(*this,&m);
        return m;
    }

    Matrix33t<T> transpose() const {
        Matrix33t<T> m;
        Transpose(*this,&m);
        return m;
    }

    Matrix33t<T>& operator*=(const T& f) {
        mat[0][0] *= f;  mat[0][1] *= f; mat[0][2] *= f;
        mat[1][0] *= f;  mat[1][1] *= f; mat[1][2] *= f;
        mat[2][0] *= f;  mat[2][1] *= f; mat[2][2] *= f;
        return *this;
    }

    friend Matrix33t<T> operator+(const Matrix33t<T>& lhs,
                                  const Matrix33t<T>& rhs) {
        return Matrix33t<T>(lhs[0] + rhs[0],
                            lhs[1] + rhs[1],
                            lhs[2] + rhs[2]);
    }

    friend Matrix33t<T> operator*(double lhs,
                                  const Matrix33t<T>& rhs) {
        return Matrix33t<T>(lhs * rhs[0],
                            lhs * rhs[1],
                            lhs * rhs[2]);
    }

    friend T Determinant(const Matrix33t<T>&);
    friend void Invert_rational_terms(const Matrix33t<T>&,
                                      Matrix33t<T>*, T*);
    friend int Invert(const Matrix33t<T>&, Matrix33t<T>*);

    friend void Transpose(const Matrix33t<T>&, Matrix33t<T>*);
};


template <class T1, class T2>
Matrix33t<T1> operator*(const Matrix33t<T1>& m1, 
                        const Matrix33t<T2>& m2);

template<class T> const T Matrix33t<T>::identity[9] = {
    1, 0, 0,
    0, 1, 0,
    0, 0, 1
};


/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

template<class T>
class Matrix44t {
    static const int MATSIZE;
public:
    T mat[4][4];

    static const T identity[16];

    Matrix44t() {}
    Matrix44t(const T *tp) {         // ctor takes array of 16 T's
        memcpy(mat,tp,MATSIZE);
    }
    Matrix44t(const T& t1,const T& t2,const T& t3,const T& t4,
              const T& t5,const T& t6,const T& t7,const T& t8,
              const T& t9,const T& t10,const T& t11,const T& t12,
              const T& t13,const T& t14,const T& t15,const T& t16) {
        mat[0][0]=t1; mat[0][1]=t2; mat[0][2]=t3; mat[0][3]=t4;
        mat[1][0]=t5; mat[1][1]=t6; mat[1][2]=t7; mat[1][3]=t8;
        mat[2][0]=t9; mat[2][1]=t10; mat[2][2]=t11; mat[2][3]=t12;
        mat[3][0]=t13; mat[3][1]=t14; mat[3][2]=t15; mat[3][3]=t16;
    }
    void set(const T& t1,const T& t2,const T& t3,const T& t4,
             const T& t5,const T& t6,const T& t7,const T& t8,
             const T& t9,const T& t10,const T& t11,const T& t12,
             const T& t13,const T& t14,const T& t15,const T& t16) {
        mat[0][0]=t1; mat[0][1]=t2; mat[0][2]=t3; mat[0][3]=t4;
        mat[1][0]=t5; mat[1][1]=t6; mat[1][2]=t7; mat[1][3]=t8;
        mat[2][0]=t9; mat[2][1]=t10; mat[2][2]=t11; mat[2][3]=t12;
        mat[3][0]=t13; mat[3][1]=t14; mat[3][2]=t15; mat[3][3]=t16;
    }
    Matrix44t(const Matrix44t& m) {  // copy ctor
        memcpy(mat,m.mat,MATSIZE);
        for (int i=0; i<16; ++i)
            ((T *)mat)[i] = ((T *)m.mat)[i];
    }

    Matrix44t<T>& operator=(const Matrix44t<T>& m) {
        memcpy(mat,m.mat,MATSIZE);
        return *this;
    }
    Matrix44t<T>& operator=(const T m[4][4]) {
        memcpy(mat,m,MATSIZE);
        return *this;
    }

    void load33(const Matrix44t<T>& m) {
        mat[0][0]=m[0][0]; mat[0][1]=m[0][1]; mat[0][2]=m[0][2];
        mat[1][0]=m[1][0]; mat[1][1]=m[1][1]; mat[1][2]=m[1][2];
        mat[2][0]=m[2][0]; mat[2][1]=m[2][1]; mat[2][2]=m[2][2];
    }
    void load33(const Matrix33t<T>& m) {
        mat[0][0]=m[0][0]; mat[0][1]=m[0][1]; mat[0][2]=m[0][2];
        mat[1][0]=m[1][0]; mat[1][1]=m[1][1]; mat[1][2]=m[1][2];
        mat[2][0]=m[2][0]; mat[2][1]=m[2][1]; mat[2][2]=m[2][2];
    }

    operator T *() const { return (T *)mat; }

    const Vector4t<T>& operator [](int i) const {
        return *(Vector4t<T> *)mat[i];
    }
    Vector4t<T>& operator [](int i) {
        return *(Vector4t<T> *)mat[i];
    }

    void setIdentity() { *this = identity; }

    void restrictedInvert() { ::restrictedInvert(*this,this); }
    Matrix44t<T>  restrictedInverse() const {
        Matrix44t<T> m;
        ::restrictedInvert(*this,&m);
        return m;
    }

    Matrix44t<T> inverse() const {
        Matrix44t<T> m;
        Invert(*this,&m);
        return m;
    }

    Matrix44t<T> transpose() const {
        Matrix44t<T> m;
        Transpose(*this,&m);
        return m;
    }

    void translate(const T &x,const T &y,const T &z);
    void translate(const T *v) { translate(v[0],v[1],v[2]); }
    Matrix44t<T> translated(const T &x,const T &y,const T &z) const {
        Matrix44t<T> m(*this);
        m.translate(x,y,z);
        return m;
    }
    Matrix44t<T> translated(const T *v) const {
        return translated(v[0],v[1],v[2]);
    }

    void rotate(char axis,const T &angle) {
        *this = rotated(axis,angle);
    }
    Matrix44t<T> rotated(char axis,const T &angle) const;

    void scale(const T &x,const T &y,const T &z);
    void scale(const T &s) { scale(s,s,s); }
    Matrix44t<T> scaled(const T &x,const T &y,const T &z) const;
    Matrix44t<T> scaled(const T &s) const { return scaled(s,s,s); }

    friend int Invert(const Matrix44t<T>&, Matrix44t<T>*);
    friend void restrictedInvert(const Matrix44t<T>&, Matrix44t<T>*);

    friend void Transpose(const Matrix44t<T>&, Matrix44t<T>*);
};

 
template <class T1, class T2>
Vector4t<T2> operator*(const Matrix44t<T1>&, const Vector4t<T2>&);

template <class T1, class T2>
Vector4t<T2> operator*(const Vector4t<T2>&, const Matrix44t<T1>&); 

template <class T1, class T2>
Vector4t<T2> Transform(const Matrix44t<T1>&, const Vector3t<T2>&);

template <class T1, class T2>
Vector4t<T2> Transform(const Vector3t<T2>&, const Matrix44t<T1>&); 

template <class T1, class T2>
Matrix44t<T1> operator*(const Matrix44t<T1>& m1, 
                        const Matrix44t<T2>& m2);

template<class T> const T Matrix44t<T>::identity[16] = {
    1, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, 1, 0,
    0, 0, 0, 1
};


/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */


const float identmat9f[9] = {
    1.0, 0.0, 0.0,
    0.0, 1.0, 0.0,
    0.0, 0.0, 1.0
};

const double identmat9d[9] = {
    1.0, 0.0, 0.0,
    0.0, 1.0, 0.0,
    0.0, 0.0, 1.0
};

const float identmat16f[16] = {
    1.0, 0.0, 0.0, 0.0,
    0.0, 1.0, 0.0, 0.0,
    0.0, 0.0, 1.0, 0.0,
    0.0, 0.0, 0.0, 1.0
};

const double identmat16d[16] = {
    1.0, 0.0, 0.0, 0.0,
    0.0, 1.0, 0.0, 0.0,
    0.0, 0.0, 1.0, 0.0,
    0.0, 0.0, 0.0, 1.0
};

typedef Vector3t<float> Vector3f;
typedef Vector3t<double> Vector3d;

typedef Vector4t<float> Vector4f;
typedef Vector4t<double> Vector4d;

typedef Matrix33t<float> Matrix33f;
typedef Matrix33t<double> Matrix33d;

typedef Matrix44t<float> Matrix44f;
typedef Matrix44t<double> Matrix44d;

typedef Vector3f Point3f;
typedef Vector3d Point3d;

typedef Vector4f Point4f;
typedef Vector4d Point4d;

const Matrix33f IdentMat33f( Matrix33f::identity );
const Matrix33d IdentMat33d( Matrix33d::identity );
const Matrix44f IdentMat44f( Matrix44f::identity );
const Matrix44d IdentMat44d( Matrix44d::identity );

const Vector3f Xaxis3f ( 1.0, 0.0, 0.0);
const Vector3f Yaxis3f ( 0.0, 1.0, 0.0);
const Vector3f Zaxis3f ( 0.0, 0.0, 1.0);
const Vector3f Origin3f( 0.0, 0.0, 0.0);
const Vector3f Zero3f  ( 0.0, 0.0, 0.0);

const Vector4f Xaxis4f ( 1.0, 0.0, 0.0, 0.0);
const Vector4f Yaxis4f ( 0.0, 1.0, 0.0, 0.0);
const Vector4f Zaxis4f ( 0.0, 0.0, 1.0, 0.0);
const Vector4f Origin4f( 0.0, 0.0, 0.0, 1.0);
const Vector4f Zero4f  ( 0.0, 0.0, 0.0, 0.0);

const Vector3d Xaxis3d ( 1.0, 0.0, 0.0);
const Vector3d Yaxis3d ( 0.0, 1.0, 0.0);
const Vector3d Zaxis3d ( 0.0, 0.0, 1.0);
const Vector3d Origin3d( 0.0, 0.0, 0.0);
const Vector3d Zero3d  ( 0.0, 0.0, 0.0);

const Vector4d Xaxis4d ( 1.0, 0.0, 0.0, 0.0);
const Vector4d Yaxis4d ( 0.0, 1.0, 0.0, 0.0);
const Vector4d Zaxis4d ( 0.0, 0.0, 1.0, 0.0);
const Vector4d Origin4d( 0.0, 0.0, 0.0, 1.0);
const Vector4d Zero4d  ( 0.0, 0.0, 0.0, 0.0);

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */


/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
#endif