diff --git a/gl3n/linalg.d b/gl3n/linalg.d
deleted file mode 100644
index 450c82f..0000000
--- a/gl3n/linalg.d
+++ /dev/null
@@ -1,2596 +0,0 @@
-/**
-gl3n.linalg
-
-Special thanks to:
-$(UL
-  $(LI Tomasz Stachowiak (h3r3tic): allowed me to use parts of $(LINK2 https://bitbucket.org/h3r3tic/boxen/src/default/src/xf/omg, omg).)
-  $(LI Jakob Øvrum (jA_cOp): improved the code a lot!)
-  $(LI Florian Boesch (___doc__): helps me to understand opengl/complex maths better, see: $(LINK http://codeflow.org/).)
-  $(LI #D on freenode: answered general questions about D.)
-)
-
-Authors: David Herberth
-License: MIT
-
-Note: All methods marked with pure are weakly pure since, they all access an instance member.
-All static methods are strongly pure.
-*/
-
-
-module gl3n.linalg;
-
-private {
-    import std.math : isNaN, isInfinity;
-    import std.conv : to;
-    import std.traits : isIntegral, isFloatingPoint, isStaticArray, isDynamicArray, isImplicitlyConvertible, isArray;
-    import std.string : format, rightJustify;
-    import std.array : join;
-    import std.algorithm : max, min, reduce;
-    import gl3n.math : clamp, PI, sqrt, sin, cos, acos, tan, asin, atan2, almost_equal;
-    import gl3n.util : is_vector, is_matrix, is_quaternion, TupleRange;
-}
-
-version(NoReciprocalMul) {
-    private enum rmul = false;
-} else {
-    private enum rmul = true;
-}
-
-/// Base template for all vector-types.
-/// Params:
-/// type = all values get stored as this type
-/// dimension = specifies the dimension of the vector, can be 1, 2, 3 or 4
-/// Examples:
-/// ---
-/// alias Vector!(int, 3) vec3i;
-/// alias Vector!(float, 4) vec4;
-/// alias Vector!(real, 2) vec2r;
-/// ---
-struct Vector(type, int dimension_) {
-    static assert(dimension > 0, "0 dimensional vectors don't exist.");
-
-    alias type vt; /// Holds the internal type of the vector.
-    static const int dimension = dimension_; ///Holds the dimension of the vector.
-
-    vt[dimension] vector; /// Holds all coordinates, length conforms dimension.
-
-    /// Returns a pointer to the coordinates.
-    @property auto value_ptr() const { return vector.ptr; }
-
-    /// Returns the current vector formatted as string, useful for printing the vector.
-    @property string as_string() {
-        return format("%s", vector);
-    }
-    alias as_string toString; /// ditto
-
-    @safe pure nothrow:
-    ///
-    private @property ref inout(vt) get_(char coord)() inout {
-        return vector[coord_to_index!coord];
-    }
-
-    alias get_!'x' x; /// static properties to access the values.
-    alias x u; /// ditto
-    alias x s; /// ditto
-    alias x r; /// ditto
-    static if(dimension >= 2) {
-        alias get_!'y' y; /// ditto
-        alias y v; /// ditto
-        alias y t; /// ditto
-        alias y g; /// ditto
-    }
-    static if(dimension >= 3) {
-        alias get_!'z' z; /// ditto
-        alias z b; /// ditto
-        alias z p; /// ditto
-    }
-    static if(dimension >= 4) {
-        alias get_!'w' w; /// ditto
-        alias w a; /// ditto
-        alias w q; /// ditto
-    }
-
-    static if(dimension == 2) {
-        enum Vector e1 = Vector(1.to!vt, 0.to!vt); /// canonical basis for Euclidian space
-        enum Vector e2 = Vector(0.to!vt, 1.to!vt); /// ditto
-    } else static if(dimension == 3) {
-        enum Vector e1 = Vector(1.to!vt, 0.to!vt, 0.to!vt); /// canonical basis for Euclidian space
-        enum Vector e2 = Vector(0.to!vt, 1.to!vt, 0.to!vt); /// ditto
-        enum Vector e3 = Vector(0.to!vt, 0.to!vt, 1.to!vt); /// ditto
-    } else static if(dimension == 4) {
-        enum Vector e1 = Vector(1.to!vt, 0.to!vt, 0.to!vt, 0.to!vt); /// canonical basis for Euclidian space
-        enum Vector e2 = Vector(0.to!vt, 1.to!vt, 0.to!vt, 0.to!vt); /// ditto
-        enum Vector e3 = Vector(0.to!vt, 0.to!vt, 1.to!vt, 0.to!vt); /// ditto
-        enum Vector e4 = Vector(0.to!vt, 0.to!vt, 0.to!vt, 1.to!vt); /// ditto
-    }
-
-    unittest {
-        assert(vec2.e1.vector == [1.0, 0.0]);
-        assert(vec2.e2.vector == [0.0, 1.0]);
-
-        assert(vec3.e1.vector == [1.0, 0.0, 0.0]);
-        assert(vec3.e2.vector == [0.0, 1.0, 0.0]);
-        assert(vec3.e3.vector == [0.0, 0.0, 1.0]);
-
-        assert(vec4.e1.vector == [1.0, 0.0, 0.0, 0.0]);
-        assert(vec4.e2.vector == [0.0, 1.0, 0.0, 0.0]);
-        assert(vec4.e3.vector == [0.0, 0.0, 1.0, 0.0]);
-        assert(vec4.e4.vector == [0.0, 0.0, 0.0, 1.0]);
-    }
-
-    static void isCompatibleVectorImpl(int d)(Vector!(vt, d) vec) if(d <= dimension) {
-    }
-
-    template isCompatibleVector(T) {
-        enum isCompatibleVector = is(typeof(isCompatibleVectorImpl(T.init)));
-    }
-
-    static void isCompatibleMatrixImpl(int r, int c)(Matrix!(vt, r, c) m) {
-    }
-
-    template isCompatibleMatrix(T) {
-        enum isCompatibleMatrix = is(typeof(isCompatibleMatrixImpl(T.init)));
-    }
-
-    private void construct(int i, T, Tail...)(T head, Tail tail) {
-        static if(i >= dimension) {
-            static assert(false, "Too many arguments passed to constructor");
-        } else static if(is(T : vt)) {
-            vector[i] = head;
-            construct!(i + 1)(tail);
-        } else static if(isDynamicArray!T) {
-            static assert((Tail.length == 0) && (i == 0), "dynamic array can not be passed together with other arguments");
-            vector[] = head[];
-        } else static if(isStaticArray!T) {
-            vector[i .. i + T.length] = head[];
-            construct!(i + T.length)(tail);
-        } else static if(isCompatibleVector!T) {
-            vector[i .. i + T.dimension] = head.vector[];
-            construct!(i + T.dimension)(tail);
-        } else {
-            static assert(false, "Vector constructor argument must be of type " ~ vt.stringof ~ " or Vector, not " ~ T.stringof);
-        }
-    }
-
-    private void construct(int i)() { // terminate
-        static assert(i == dimension, "Not enough arguments passed to constructor");
-    }
-
-    /// Constructs the vector.
-    /// If a single value is passed the vector, the vector will be cleared with this value.
-    /// If a vector with a higher dimension is passed the vector will hold the first values up to its dimension.
-    /// If mixed types are passed they will be joined together (allowed types: vector, static array, $(I vt)).
-    /// Examples:
-    /// ---
-    /// vec4 v4 = vec4(1.0f, vec2(2.0f, 3.0f), 4.0f);
-    /// vec3 v3 = vec3(v4); // v3 = vec3(1.0f, 2.0f, 3.0f);
-    /// vec2 v2 = v3.xy; // swizzling returns a static array.
-    /// vec3 v3_2 = vec3(1.0f); // vec3 v3_2 = vec3(1.0f, 1.0f, 1.0f);
-    /// ---
-    this(Args...)(Args args) {
-        construct!(0)(args);
-    }
-
-    /// ditto
-    this(T)(T vec) if(is_vector!T && is(T.vt : vt) && (T.dimension >= dimension)) {
-        foreach(i; TupleRange!(0, dimension)) {
-            vector[i] = vec.vector[i];
-        }
-    }
-
-    /// ditto
-    this()(vt value) {
-        clear(value);
-    }
-
-    /// Returns true if all values are not nan and finite, otherwise false.
-    @property bool isFinite() const {
-        static if(isIntegral!type) {
-            return true;
-        }
-        else {
-            foreach(v; vector) {
-                if(isNaN(v) || isInfinity(v)) {
-                    return false;
-                }
-            }
-            return true;
-        }
-    }
-    deprecated("Use isFinite instead of ok") alias ok = isFinite;
-
-    /// Sets all values of the vector to value.
-    void clear(vt value) {
-        foreach(i; TupleRange!(0, dimension)) {
-            vector[i] = value;
-        }
-    }
-
-    unittest {
-        vec3 vec_clear;
-        assert(!vec_clear.isFinite);
-        vec_clear.clear(1.0f);
-        assert(vec_clear.isFinite);
-        assert(vec_clear.vector == [1.0f, 1.0f, 1.0f]);
-        assert(vec_clear.vector == vec3(1.0f).vector);
-        vec_clear.clear(float.infinity);
-        assert(!vec_clear.isFinite);
-        vec_clear.clear(float.nan);
-        assert(!vec_clear.isFinite);
-        vec_clear.clear(1.0f);
-        assert(vec_clear.isFinite);
-
-        vec4 b = vec4(1.0f, vec_clear);
-        assert(b.isFinite);
-        assert(b.vector == [1.0f, 1.0f, 1.0f, 1.0f]);
-        assert(b.vector == vec4(1.0f).vector);
-
-        vec2 v2_1 = vec2(vec2(0.0f, 1.0f));
-        assert(v2_1.vector == [0.0f, 1.0f]);
-
-        vec2 v2_2 = vec2(1.0f, 1.0f);
-        assert(v2_2.vector == [1.0f, 1.0f]);
-
-        vec3 v3 = vec3(v2_1, 2.0f);
-        assert(v3.vector == [0.0f, 1.0f, 2.0f]);
-
-        vec4 v4_1 = vec4(1.0f, vec2(2.0f, 3.0f), 4.0f);
-        assert(v4_1.vector == [1.0f, 2.0f, 3.0f, 4.0f]);
-        assert(vec3(v4_1).vector == [1.0f, 2.0f, 3.0f]);
-        assert(vec2(vec3(v4_1)).vector == [1.0f, 2.0f]);
-        assert(vec2(vec3(v4_1)).vector == vec2(v4_1).vector);
-        assert(v4_1.vector == vec4([1.0f, 2.0f, 3.0f, 4.0f]).vector);
-
-        vec4 v4_2 = vec4(vec2(1.0f, 2.0f), vec2(3.0f, 4.0f));
-        assert(v4_2.vector == [1.0f, 2.0f, 3.0f, 4.0f]);
-        assert(vec3(v4_2).vector == [1.0f, 2.0f, 3.0f]);
-        assert(vec2(vec3(v4_2)).vector == [1.0f, 2.0f]);
-        assert(vec2(vec3(v4_2)).vector == vec2(v4_2).vector);
-        assert(v4_2.vector == vec4([1.0f, 2.0f, 3.0f, 4.0f]).vector);
-
-        float[2] f2 = [1.0f, 2.0f];
-        float[3] f3 = [1.0f, 2.0f, 3.0f];
-        float[4] f4 = [1.0f, 2.0f, 3.0f, 4.0f];
-        assert(vec2(1.0f, 2.0f).vector == vec2(f2).vector);
-        assert(vec3(1.0f, 2.0f, 3.0f).vector == vec3(f3).vector);
-        assert(vec3(1.0f, 2.0f, 3.0f).vector == vec3(f2, 3.0f).vector);
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f).vector == vec4(f4).vector);
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f).vector == vec4(f3, 4.0f).vector);
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f).vector == vec4(f2, 3.0f, 4.0f).vector);
-        // useful for: "vec4 v4 = […]" or "vec4 v4 = other_vector.rgba"
-
-        assert(vec3(vec3i(1, 2, 3)) == vec3(1.0, 2.0, 3.0));
-        assert(vec3d(vec3(1.0, 2.0, 3.0)) == vec3d(1.0, 2.0, 3.0));
-
-        static assert(!__traits(compiles, vec3(0.0f, 0.0f)));
-        static assert(!__traits(compiles, vec4(0.0f, 0.0f, 0.0f)));
-        static assert(!__traits(compiles, vec4(0.0f, vec2(0.0f, 0.0f))));
-        static assert(!__traits(compiles, vec4(vec3(0.0f, 0.0f, 0.0f))));
-    }
-
-    template coord_to_index(char c) {
-        static if((c == 'x') || (c == 'r') || (c == 'u') || (c == 's')) {
-            enum coord_to_index = 0;
-        } else static if((c == 'y') || (c == 'g') || (c == 'v') || (c == 't')) {
-            enum coord_to_index = 1;
-        } else static if((c == 'z') || (c == 'b') || (c == 'p')) {
-            static assert(dimension >= 3, "the " ~ c ~ " property is only available on vectors with a third dimension.");
-            enum coord_to_index = 2;
-        } else static if((c == 'w') || (c == 'a') || (c == 'q')) {
-            static assert(dimension >= 4, "the " ~ c ~ " property is only available on vectors with a fourth dimension.");
-            enum coord_to_index = 3;
-        } else {
-            static assert(false, "accepted coordinates are x, s, r, u, y, g, t, v, z, p, b, w, q and a not " ~ c ~ ".");
-        }
-    }
-
-    static if(dimension == 2) { void set(vt x, vt y) { vector[0] = x; vector[1] = y; } }
-    static if(dimension == 3) { void set(vt x, vt y, vt z) { vector[0] = x; vector[1] = y; vector[2] = z; } }
-    static if(dimension == 4) { void set(vt x, vt y, vt z, vt w) { vector[0] = x; vector[1] = y; vector[2] = z; vector[3] = w; } }
-
-    /// Updates the vector with the values from other.
-    void update(Vector!(vt, dimension) other) {
-        vector = other.vector;
-    }
-
-    unittest {
-        vec2 v2 = vec2(1.0f, 2.0f);
-        assert(v2.x == 1.0f);
-        assert(v2.y == 2.0f);
-        v2.x = 3.0f;
-        v2.x += 1;
-        v2.x -= 1;
-        assert(v2.vector == [3.0f, 2.0f]);
-        v2.y = 4.0f;
-        v2.y += 1;
-        v2.y -= 1;
-        assert(v2.vector == [3.0f, 4.0f]);
-        assert((v2.x == 3.0f) && (v2.x == v2.u) && (v2.x == v2.s) && (v2.x == v2.r));
-        assert(v2.y == 4.0f);
-        assert((v2.y == 4.0f) && (v2.y == v2.v) && (v2.y == v2.t) && (v2.y == v2.g));
-        v2.set(0.0f, 1.0f);
-        assert(v2.vector == [0.0f, 1.0f]);
-        v2.update(vec2(3.0f, 4.0f));
-        assert(v2.vector == [3.0f, 4.0f]);
-
-        vec3 v3 = vec3(1.0f, 2.0f, 3.0f);
-        assert(v3.x == 1.0f);
-        assert(v3.y == 2.0f);
-        assert(v3.z == 3.0f);
-        v3.x = 3.0f;
-        v3.x += 1;
-        v3.x -= 1;
-        assert(v3.vector == [3.0f, 2.0f, 3.0f]);
-        v3.y = 4.0f;
-        v3.y += 1;
-        v3.y -= 1;
-        assert(v3.vector == [3.0f, 4.0f, 3.0f]);
-        v3.z = 5.0f;
-        v3.z += 1;
-        v3.z -= 1;
-        assert(v3.vector == [3.0f, 4.0f, 5.0f]);
-        assert((v3.x == 3.0f) && (v3.x == v3.s) && (v3.x == v3.r));
-        assert((v3.y == 4.0f) && (v3.y == v3.t) && (v3.y == v3.g));
-        assert((v3.z == 5.0f) && (v3.z == v3.p) && (v3.z == v3.b));
-        v3.set(0.0f, 1.0f, 2.0f);
-        assert(v3.vector == [0.0f, 1.0f, 2.0f]);
-        v3.update(vec3(3.0f, 4.0f, 5.0f));
-        assert(v3.vector == [3.0f, 4.0f, 5.0f]);
-
-        vec4 v4 = vec4(1.0f, 2.0f, vec2(3.0f, 4.0f));
-        assert(v4.x == 1.0f);
-        assert(v4.y == 2.0f);
-        assert(v4.z == 3.0f);
-        assert(v4.w == 4.0f);
-        v4.x = 3.0f;
-        v4.x += 1;
-        v4.x -= 1;
-        assert(v4.vector == [3.0f, 2.0f, 3.0f, 4.0f]);
-        v4.y = 4.0f;
-        v4.y += 1;
-        v4.y -= 1;
-        assert(v4.vector == [3.0f, 4.0f, 3.0f, 4.0f]);
-        v4.z = 5.0f;
-        v4.z += 1;
-        v4.z -= 1;
-        assert(v4.vector == [3.0f, 4.0f, 5.0f, 4.0f]);
-        v4.w = 6.0f;
-        v4.w += 1;
-        v4.w -= 1;
-        assert(v4.vector == [3.0f, 4.0f, 5.0f, 6.0f]);
-        assert((v4.x == 3.0f) && (v4.x == v4.s) && (v4.x == v4.r));
-        assert((v4.y == 4.0f) && (v4.y == v4.t) && (v4.y == v4.g));
-        assert((v4.z == 5.0f) && (v4.z == v4.p) && (v4.z == v4.b));
-        assert((v4.w == 6.0f) && (v4.w == v4.q) && (v4.w == v4.a));
-        v4.set(0.0f, 1.0f, 2.0f, 3.0f);
-        assert(v4.vector == [0.0f, 1.0f, 2.0f, 3.0f]);
-        v4.update(vec4(3.0f, 4.0f, 5.0f, 6.0f));
-        assert(v4.vector == [3.0f, 4.0f, 5.0f, 6.0f]);
-    }
-
-    private void dispatchImpl(int i, string s, int size)(ref vt[size] result) const {
-        static if(s.length > 0) {
-            result[i] = vector[coord_to_index!(s[0])];
-            dispatchImpl!(i + 1, s[1..$])(result);
-        }
-    }
-
-    /// Implements dynamic swizzling.
-    /// Returns: a Vector
-    @property Vector!(vt, s.length) opDispatch(string s)() const {
-        vt[s.length] ret;
-        dispatchImpl!(0, s)(ret);
-        Vector!(vt, s.length) ret_vec;
-        ret_vec.vector = ret;
-        return ret_vec;
-    }
-
-    unittest {
-        vec2 v2 = vec2(1.0f, 2.0f);
-        assert(v2.xytsy == [1.0f, 2.0f, 2.0f, 1.0f, 2.0f]);
-
-        assert(vec3(1.0f, 2.0f, 3.0f).xybzyr == [1.0f, 2.0f, 3.0f, 3.0f, 2.0f, 1.0f]);
-        assert(vec4(v2, 3.0f, 4.0f).xyzwrgbastpq == [1.0f, 2.0f, 3.0f, 4.0f,
-                                                     1.0f, 2.0f, 3.0f, 4.0f,
-                                                     1.0f, 2.0f, 3.0f, 4.0f]);
-        assert(vec4(v2, 3.0f, 4.0f).wgyzax == [4.0f, 2.0f, 2.0f, 3.0f, 4.0f, 1.0f]);
-        assert(vec4(v2.xyst).vector == [1.0f, 2.0f, 1.0f, 2.0f]);
-    }
-
-    /// Returns the squared magnitude of the vector.
-    @property real magnitude_squared() const {
-        real temp = 0;
-
-        foreach(index; TupleRange!(0, dimension)) {
-            temp += vector[index]^^2;
-        }
-
-        return temp;
-    }
-
-    /// Returns the magnitude of the vector.
-    @property real magnitude() const {
-        return sqrt(magnitude_squared);
-    }
-
-    alias magnitude_squared length_squared; /// ditto
-    alias magnitude length; /// ditto
-
-    /// Normalizes the vector.
-    void normalize() {
-        real len = length;
-
-        if(len != 0) {
-            foreach(index; TupleRange!(0, dimension)) {
-                vector[index] = cast(type)(vector[index]/len);
-            }
-        }
-    }
-
-    /// Returns a normalized copy of the current vector.
-    @property Vector normalized() const {
-        Vector ret;
-        ret.update(this);
-        ret.normalize();
-        return ret;
-    }
-
-    Vector opUnary(string op : "-")() const {
-        Vector ret;
-
-        foreach(index; TupleRange!(0, dimension)) {
-            ret.vector[index] = -vector[index];
-        }
-
-        return ret;
-    }
-
-    unittest {
-        assert(vec2(1.0f, 1.0f) == -vec2(-1.0f, -1.0f));
-        assert(vec2(-1.0f, 1.0f) == -vec2(1.0f, -1.0f));
-
-        assert(-vec3(1.0f, 1.0f, 1.0f) == vec3(-1.0f, -1.0f, -1.0f));
-        assert(-vec3(-1.0f, 1.0f, -1.0f) == vec3(1.0f, -1.0f, 1.0f));
-
-        assert(vec4(1.0f, 1.0f, 1.0f, 1.0f) == -vec4(-1.0f, -1.0f, -1.0f, -1.0f));
-        assert(vec4(-1.0f, 1.0f, -1.0f, 1.0f) == -vec4(1.0f, -1.0f, 1.0f, -1.0f));
-    }
-
-    // let the math begin!
-    Vector opBinary(string op : "*")(vt r) const {
-        Vector ret;
-
-        foreach(index; TupleRange!(0, dimension)) {
-            ret.vector[index] = vector[index] * r;
-        }
-
-        return ret;
-    }
-
-    Vector opBinary(string op : "/")(vt r) const {
-        Vector ret;
-
-        foreach(index; TupleRange!(0, dimension)) {
-            ret.vector[index] = vector[index] / r;
-        }
-
-        return ret;
-    }
-
-    Vector opBinary(string op)(Vector r) const if((op == "+") || (op == "-")) {
-        Vector ret;
-
-        foreach(index; TupleRange!(0, dimension)) {
-            ret.vector[index] = mixin("vector[index]" ~ op ~ "r.vector[index]");
-        }
-
-        return ret;
-    }
-
-    vt opBinary(string op : "*")(Vector r) const {
-        return dot(this, r);
-    }
-
-    // vector * matrix (for matrix * vector -> struct Matrix)
-    Vector!(vt, T.cols) opBinary(string op : "*", T)(T inp) const if(isCompatibleMatrix!T && (T.rows == dimension)) {
-        Vector!(vt, T.cols) ret;
-        ret.clear(0);
-
-        foreach(c; TupleRange!(0, T.cols)) {
-            foreach(r; TupleRange!(0, T.rows)) {
-                ret.vector[c] += vector[r] * inp.matrix[r][c];
-            }
-        }
-
-        return ret;
-    }
-
-    auto opBinaryRight(string op, T)(T inp) const if(!is_vector!T && !is_matrix!T && !is_quaternion!T) {
-        return this.opBinary!(op)(inp);
-    }
-
-    unittest {
-        vec2 v2 = vec2(1.0f, 3.0f);
-        auto v2times2 = 2 * v2;
-        assert((v2*2.5f).vector == [2.5f, 7.5f]);
-        assert((v2+vec2(3.0f, 1.0f)).vector == [4.0f, 4.0f]);
-        assert((v2-vec2(1.0f, 3.0f)).vector == [0.0f, 0.0f]);
-        assert((v2*vec2(2.0f, 2.0f)) == 8.0f);
-
-        vec3 v3 = vec3(1.0f, 3.0f, 5.0f);
-        assert((v3*2.5f).vector == [2.5f, 7.5f, 12.5f]);
-        assert((v3+vec3(3.0f, 1.0f, -1.0f)).vector == [4.0f, 4.0f, 4.0f]);
-        assert((v3-vec3(1.0f, 3.0f, 5.0f)).vector == [0.0f, 0.0f, 0.0f]);
-        assert((v3*vec3(2.0f, 2.0f, 2.0f)) == 18.0f);
-
-        vec4 v4 = vec4(1.0f, 3.0f, 5.0f, 7.0f);
-        assert((v4*2.5f).vector == [2.5f, 7.5f, 12.5f, 17.5]);
-        assert((v4+vec4(3.0f, 1.0f, -1.0f, -3.0f)).vector == [4.0f, 4.0f, 4.0f, 4.0f]);
-        assert((v4-vec4(1.0f, 3.0f, 5.0f, 7.0f)).vector == [0.0f, 0.0f, 0.0f, 0.0f]);
-        assert((v4*vec4(2.0f, 2.0f, 2.0f, 2.0f)) == 32.0f);
-
-        mat2 m2 = mat2(1.0f, 2.0f, 3.0f, 4.0f);
-        vec2 v2_2 = vec2(2.0f, 2.0f);
-        assert((v2_2*m2).vector == [8.0f, 12.0f]);
-
-        mat3 m3 = mat3(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f);
-        vec3 v3_2 = vec3(2.0f, 2.0f, 2.0f);
-        assert((v3_2*m3).vector == [24.0f, 30.0f, 36.0f]);
-    }
-
-    void opOpAssign(string op : "*")(vt r) {
-        foreach(index; TupleRange!(0, dimension)) {
-            vector[index] *= r;
-        }
-    }
-
-    void opOpAssign(string op : "/")(vt r) {
-        foreach(index; TupleRange!(0, dimension)) {
-            vector[index] /= r;
-        }
-    }
-
-    void opOpAssign(string op)(Vector r) if((op == "+") || (op == "-")) {
-        foreach(index; TupleRange!(0, dimension)) {
-            mixin("vector[index]" ~ op ~ "= r.vector[index];");
-        }
-    }
-
-    unittest {
-        vec2 v2 = vec2(1.0f, 3.0f);
-        v2 *= 2.5f;
-        assert(v2.vector == [2.5f, 7.5f]);
-        v2 -= vec2(2.5f, 7.5f);
-        assert(v2.vector == [0.0f, 0.0f]);
-        v2 += vec2(1.0f, 3.0f);
-        assert(v2.vector == [1.0f, 3.0f]);
-        assert(almost_equal(v2.length, sqrt(10.0f)));
-        assert(v2.length_squared == 10.0f);
-        assert((v2.magnitude == v2.length) && (v2.magnitude_squared == v2.length_squared));
-        v2 /= 2.0f;
-        assert(v2.vector == [0.5f, 1.5f]);
-        assert(almost_equal(v2.normalized, vec2(1.0f/sqrt(10.0f), 3.0f/sqrt(10.0f))));
-
-        vec3 v3 = vec3(1.0f, 3.0f, 5.0f);
-        v3 *= 2.5f;
-        assert(v3.vector == [2.5f, 7.5f, 12.5f]);
-        v3 -= vec3(2.5f, 7.5f, 12.5f);
-        assert(v3.vector == [0.0f, 0.0f, 0.0f]);
-        v3 += vec3(1.0f, 3.0f, 5.0f);
-        assert(v3.vector == [1.0f, 3.0f, 5.0f]);
-        assert(almost_equal(v3.length, sqrt(35.0f)));
-        assert(v3.length_squared == 35.0f);
-        assert((v3.magnitude == v3.length) && (v3.magnitude_squared == v3.length_squared));
-        v3 /= 2.0f;
-        assert(v3.vector == [0.5f, 1.5f, 2.5f]);
-        assert(almost_equal(v3.normalized, vec3(1.0f/sqrt(35.0f), 3.0f/sqrt(35.0f), 5.0f/sqrt(35.0f))));
-
-        vec4 v4 = vec4(1.0f, 3.0f, 5.0f, 7.0f);
-        v4 *= 2.5f;
-        assert(v4.vector == [2.5f, 7.5f, 12.5f, 17.5]);
-        v4 -= vec4(2.5f, 7.5f, 12.5f, 17.5f);
-        assert(v4.vector == [0.0f, 0.0f, 0.0f, 0.0f]);
-        v4 += vec4(1.0f, 3.0f, 5.0f, 7.0f);
-        assert(v4.vector == [1.0f, 3.0f, 5.0f, 7.0f]);
-        assert(almost_equal(v4.length, sqrt(84.0f)));
-        assert(v4.length_squared == 84.0f);
-        assert((v4.magnitude == v4.length) && (v4.magnitude_squared == v4.length_squared));
-        v4 /= 2.0f;
-        assert(v4.vector == [0.5f, 1.5f, 2.5f, 3.5f]);
-        assert(almost_equal(v4.normalized, vec4(1.0f/sqrt(84.0f), 3.0f/sqrt(84.0f), 5.0f/sqrt(84.0f), 7.0f/sqrt(84.0f))));
-    }
-
-    int opCmp(ref const Vector vec) const {
-        foreach(i, a; vector) {
-            if(a < vec.vector[i]) {
-                return -1;
-            } else if(a > vec.vector[i]) {
-                return 1;
-            }
-        }
-
-        // Vectors are the same
-        return 0;
-    }
-
-    bool opEquals(T)(const T vec) const if(!isArray!T && T.dimension == dimension) {
-        return vector == vec.vector;
-    }
-
-    bool opEquals(T)(const(T)[] array) const if(!isArray!T && !is_vector!T) {
-        if(array.length != dimension) {
-            return false;
-        }
-
-        foreach(index; TupleRange!(0, dimension)) {
-            if(vector[index] != array[index]) {
-                return false;
-            }
-        }
-
-        return true;
-    }
-
-    bool opCast(T : bool)() const {
-        return isFinite;
-    }
-
-    unittest {
-        assert(vec2(1.0f, 2.0f) == vec2(1.0f, 2.0f));
-        assert(vec2(1.0f, 2.0f) != vec2(1.0f, 1.0f));
-        assert(vec2(1.0f, 2.0f) == vec2d(1.0, 2.0));
-        assert(vec2(1.0f, 2.0f) != vec2d(1.0, 1.0));
-        assert(vec2(1.0f, 2.0f) == vec2(1.0f, 2.0f).vector);
-        assert(vec2(1.0f, 2.0f) != vec2(1.0f, 1.0f).vector);
-        assert(vec2(1.0f, 2.0f) == vec2d(1.0, 2.0).vector);
-        assert(vec2(1.0f, 2.0f) != vec2d(1.0, 1.0).vector);
-
-        assert(vec3(1.0f, 2.0f, 3.0f) == vec3(1.0f, 2.0f, 3.0f));
-        assert(vec3(1.0f, 2.0f, 3.0f) != vec3(1.0f, 2.0f, 2.0f));
-        assert(vec3(1.0f, 2.0f, 3.0f) == vec3d(1.0, 2.0, 3.0));
-        assert(vec3(1.0f, 2.0f, 3.0f) != vec3d(1.0, 2.0, 2.0));
-        assert(vec3(1.0f, 2.0f, 3.0f) == vec3(1.0f, 2.0f, 3.0f).vector);
-        assert(vec3(1.0f, 2.0f, 3.0f) != vec3(1.0f, 2.0f, 2.0f).vector);
-        assert(vec3(1.0f, 2.0f, 3.0f) == vec3d(1.0, 2.0, 3.0).vector);
-        assert(vec3(1.0f, 2.0f, 3.0f) != vec3d(1.0, 2.0, 2.0).vector);
-
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4(1.0f, 2.0f, 3.0f, 4.0f));
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4(1.0f, 2.0f, 3.0f, 3.0f));
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4d(1.0, 2.0, 3.0, 4.0));
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4d(1.0, 2.0, 3.0, 3.0));
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4(1.0f, 2.0f, 3.0f, 4.0f).vector);
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4(1.0f, 2.0f, 3.0f, 3.0f).vector);
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4d(1.0, 2.0, 3.0, 4.0).vector);
-        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4d(1.0, 2.0, 3.0, 3.0).vector);
-
-        assert(!(vec4(float.nan)));
-        if(vec4(1.0f)) { }
-        else { assert(false); }
-    }
-
-}
-
-/// Calculates the product between two vectors.
-T.vt dot(T)(const T veca, const T vecb) @safe pure nothrow if(is_vector!T) {
-    T.vt temp = 0;
-
-    foreach(index; TupleRange!(0, T.dimension)) {
-        temp += veca.vector[index] * vecb.vector[index];
-    }
-
-    return temp;
-}
-
-/// Calculates the cross product of two 3-dimensional vectors.
-T cross(T)(const T veca, const T vecb) @safe pure nothrow if(is_vector!T && (T.dimension == 3)) {
-   return T(veca.y * vecb.z - vecb.y * veca.z,
-            veca.z * vecb.x - vecb.z * veca.x,
-            veca.x * vecb.y - vecb.x * veca.y);
-}
-
-/// Calculates the distance between two vectors.
-T.vt distance(T)(const T veca, const T vecb) @safe pure nothrow if(is_vector!T) {
-    return (veca - vecb).length;
-}
-
-unittest {
-    // dot is already tested in Vector.opBinary, so no need for testing with more vectors
-    vec3 v1 = vec3(1.0f, 2.0f, -3.0f);
-    vec3 v2 = vec3(1.0f, 3.0f, 2.0f);
-
-    assert(dot(v1, v2) == 1.0f);
-    assert(dot(v1, v2) == (v1 * v2));
-    assert(dot(v1, v2) == dot(v2, v1));
-    assert((v1 * v2) == (v1 * v2));
-
-    assert(cross(v1, v2).vector == [13.0f, -5.0f, 1.0f]);
-    assert(cross(v2, v1).vector == [-13.0f, 5.0f, -1.0f]);
-
-    assert(distance(vec2(0.0f, 0.0f), vec2(0.0f, 10.0f)) == 10.0);
-}
-
-/// reflect a vector using a surface normal
-T reflect(T)(const T vec, const T norm) @safe pure nothrow if(is_vector!T) {
-    return (2 * (vec * norm) * norm) - vec;
-}
-
-unittest
-{
-    assert(vec2(1,1).reflect(vec2(0,1)) == vec2(-1,1));
-    assert(vec2(-1,1).reflect(vec2(0,1)) == vec2(1,1));
-    assert(vec2(2,1).reflect(vec2(0,1)) == vec2(-2,1));
-
-    assert(vec3(1,1,1).reflect(vec3(0,1,0)) == vec3(-1,1,-1));
-}
-
-/// Pre-defined vector types, the number represents the dimension and the last letter the type (none = float, d = double, i = int).
-alias Vector!(float, 2) vec2;
-alias Vector!(float, 3) vec3; /// ditto
-alias Vector!(float, 4) vec4; /// ditto
-
-alias Vector!(double, 2) vec2d; /// ditto
-alias Vector!(double, 3) vec3d; /// ditto
-alias Vector!(double, 4) vec4d; /// ditto
-
-alias Vector!(int, 2) vec2i; /// ditto
-alias Vector!(int, 3) vec3i; /// ditto
-alias Vector!(int, 4) vec4i; /// ditto
-
-/*alias Vector!(ubyte, 2) vec2ub;
-alias Vector!(ubyte, 3) vec3ub;
-alias Vector!(ubyte, 4) vec4ub;*/
-
-
-/// Base template for all matrix-types.
-/// Params:
-///  type = all values get stored as this type
-///  rows_ = rows of the matrix
-///  cols_ = columns of the matrix
-/// Examples:
-/// ---
-/// alias Matrix!(float, 4, 4) mat4;
-/// alias Matrix!(double, 3, 4) mat34d;
-/// alias Matrix!(real, 2, 2) mat2r;
-/// ---
-struct Matrix(type, int rows_, int cols_) if((rows_ > 0) && (cols_ > 0)) {
-    alias type mt; /// Holds the internal type of the matrix;
-    static const int rows = rows_; /// Holds the number of rows;
-    static const int cols = cols_; /// Holds the number of columns;
-
-    /// Holds the matrix $(RED row-major) in memory.
-    mt[cols][rows] matrix; // In C it would be mt[rows][cols], D does it like this: (mt[foo])[bar]
-    alias matrix this;
-
-    unittest {
-        mat2 m2 = mat2(0.0f, 1.0f, 2.0f, 3.0f);
-        assert(m2[0][0] == 0.0f);
-        assert(m2[0][1] == 1.0f);
-        assert(m2[1][0] == 2.0f);
-        assert(m2[1][1] == 3.0f);
-        m2[0..1] = [2.0f, 2.0f];
-        assert(m2 == [[2.0f, 2.0f], [2.0f, 3.0f]]);
-
-        mat3 m3 = mat3(0.0f, 0.1f, 0.2f, 1.0f, 1.1f, 1.2f, 2.0f, 2.1f, 2.2f);
-        assert(m3[0][1] == 0.1f);
-        assert(m3[2][0] == 2.0f);
-        assert(m3[1][2] == 1.2f);
-        m3[0][0..$] = 0.0f;
-        assert(m3 == [[0.0f, 0.0f, 0.0f],
-                      [1.0f, 1.1f, 1.2f],
-                      [2.0f, 2.1f, 2.2f]]);
-
-        mat4 m4 = mat4(0.0f, 0.1f, 0.2f, 0.3f,
-                       1.0f, 1.1f, 1.2f, 1.3f,
-                       2.0f, 2.1f, 2.2f, 2.3f,
-                       3.0f, 3.1f, 3.2f, 3.3f);
-       assert(m4[0][3] == 0.3f);
-       assert(m4[1][1] == 1.1f);
-       assert(m4[2][0] == 2.0f);
-       assert(m4[3][2] == 3.2f);
-       m4[2][1..3] = [1.0f, 2.0f];
-       assert(m4 == [[0.0f, 0.1f, 0.2f, 0.3f],
-                     [1.0f, 1.1f, 1.2f, 1.3f],
-                     [2.0f, 1.0f, 2.0f, 2.3f],
-                     [3.0f, 3.1f, 3.2f, 3.3f]]);
-
-    }
-
-    /// Returns the pointer to the stored values as OpenGL requires it.
-    /// Note this will return a pointer to a $(RED row-major) matrix,
-    /// $(RED this means you've to set the transpose argument to GL_TRUE when passing it to OpenGL).
-    /// Examples:
-    /// ---
-    /// // 3rd argument = GL_TRUE
-    /// glUniformMatrix4fv(programs.main.model, 1, GL_TRUE, mat4.translation(-0.5f, -0.5f, 1.0f).value_ptr);
-    /// ---
-    @property auto value_ptr() const { return matrix[0].ptr; }
-
-    /// Returns the current matrix formatted as flat string.
-    @property string as_string() {
-        return format("%s", matrix);
-    }
-    alias as_string toString; /// ditto
-
-    /// Returns the current matrix as pretty formatted string.
-    @property string as_pretty_string() {
-        string fmtr = "%s";
-
-        size_t rjust = max(format(fmtr, reduce!(max)(matrix[])).length,
-                           format(fmtr, reduce!(min)(matrix[])).length) - 1;
-
-        string[] outer_parts;
-        foreach(mt[] row; matrix) {
-            string[] inner_parts;
-            foreach(mt col; row) {
-                inner_parts ~= rightJustify(format(fmtr, col), rjust);
-            }
-            outer_parts ~= " [" ~ join(inner_parts, ", ") ~ "]";
-        }
-
-        return "[" ~ join(outer_parts, "\n")[1..$] ~ "]";
-    }
-    alias as_pretty_string toPrettyString; /// ditto
-
-    @safe pure nothrow:
-    static void isCompatibleMatrixImpl(int r, int c)(Matrix!(mt, r, c) m) {
-    }
-
-    template isCompatibleMatrix(T) {
-        enum isCompatibleMatrix = is(typeof(isCompatibleMatrixImpl(T.init)));
-    }
-
-    static void isCompatibleVectorImpl(int d)(Vector!(mt, d) vec) {
-    }
-
-    template isCompatibleVector(T) {
-        enum isCompatibleVector = is(typeof(isCompatibleVectorImpl(T.init)));
-    }
-
-    private void construct(int i, T, Tail...)(T head, Tail tail) {
-        static if(i >= rows*cols) {
-            static assert(false, "Too many arguments passed to constructor");
-        } else static if(is(T : mt)) {
-            matrix[i / cols][i % cols] = head;
-            construct!(i + 1)(tail);
-        } else static if(is(T == Vector!(mt, cols))) {
-            static if(i % cols == 0) {
-                matrix[i / cols] = head.vector;
-                construct!(i + T.dimension)(tail);
-            } else {
-                static assert(false, "Can't convert Vector into the matrix. Maybe it doesn't align to the columns correctly or dimension doesn't fit");
-            }
-        } else static if(isDynamicArray!T) {
-            foreach(j; 0..cols*rows)
-                matrix[j / cols][j % cols] = head[j];
-        } else {
-            static assert(false, "Matrix constructor argument must be of type " ~ mt.stringof ~ " or Vector, not " ~ T.stringof);
-        }
-    }
-
-    private void construct(int i)() { // terminate
-        static assert(i == rows*cols, "Not enough arguments passed to constructor");
-    }
-
-    /// Constructs the matrix:
-    /// If a single value is passed, the matrix will be cleared with this value (each column in each row will contain this value).
-    /// If a matrix with more rows and columns is passed, the matrix will be the upper left nxm matrix.
-    /// If a matrix with less rows and columns is passed, the passed matrix will be stored in the upper left of an identity matrix.
-    /// It's also allowed to pass vectors and scalars at a time, but the vectors dimension must match the number of columns and align correctly.
-    /// Examples:
-    /// ---
-    /// mat2 m2 = mat2(0.0f); // mat2 m2 = mat2(0.0f, 0.0f, 0.0f, 0.0f);
-    /// mat3 m3 = mat3(m2); // mat3 m3 = mat3(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f);
-    /// mat3 m3_2 = mat3(vec3(1.0f, 2.0f, 3.0f), 4.0f, 5.0f, 6.0f, vec3(7.0f, 8.0f, 9.0f));
-    /// mat4 m4 = mat4.identity; // just an identity matrix
-    /// mat3 m3_3 = mat3(m4); // mat3 m3_3 = mat3.identity
-    /// ---
-    this(Args...)(Args args) {
-        construct!(0)(args);
-    }
-
-    /// ditto
-    this(T)(T mat) if(is_matrix!T && (T.cols >= cols) && (T.rows >= rows)) {
-        foreach(r; TupleRange!(0, rows)) {
-            foreach(c; TupleRange!(0, cols)) {
-                matrix[r][c] = mat.matrix[r][c];
-            }
-        }
-    }
-
-    /// ditto
-    this(T)(T mat) if(is_matrix!T && (T.cols < cols) && (T.rows < rows)) {
-        make_identity();
-
-        foreach(r; TupleRange!(0, T.rows)) {
-            foreach(c; TupleRange!(0, T.cols)) {
-                matrix[r][c] = mat.matrix[r][c];
-            }
-        }
-    }
-
-    /// ditto
-    this()(mt value) {
-        clear(value);
-    }
-
-    /// Returns true if all values are not nan and finite, otherwise false.
-    @property bool isFinite() const {
-        static if(isIntegral!type) {
-            return true;
-        }
-        else {
-            foreach(row; matrix) {
-                foreach(col; row) {
-                    if(isNaN(col) || isInfinity(col)) {
-                        return false;
-                    }
-                }
-            }
-            return true;
-        }
-
-    }
-    deprecated("Use isFinite instead of ok") alias ok = isFinite;
-
-    /// Sets all values of the matrix to value (each column in each row will contain this value).
-    void clear(mt value) {
-        foreach(r; TupleRange!(0, rows)) {
-            foreach(c; TupleRange!(0, cols)) {
-                matrix[r][c] = value;
-            }
-        }
-    }
-
-    unittest {
-        mat2 m2 = mat2(1.0f, 1.0f, vec2(2.0f, 2.0f));
-        assert(m2.matrix == [[1.0f, 1.0f], [2.0f, 2.0f]]);
-        m2.clear(3.0f);
-        assert(m2.matrix == [[3.0f, 3.0f], [3.0f, 3.0f]]);
-        assert(m2.isFinite);
-        m2.clear(float.nan);
-        assert(!m2.isFinite);
-        m2.clear(float.infinity);
-        assert(!m2.isFinite);
-        m2.clear(0.0f);
-        assert(m2.isFinite);
-
-        mat3 m3 = mat3(1.0f);
-        assert(m3.matrix == [[1.0f, 1.0f, 1.0f],
-                             [1.0f, 1.0f, 1.0f],
-                             [1.0f, 1.0f, 1.0f]]);
-
-        mat4 m4 = mat4(vec4(1.0f, 1.0f, 1.0f, 1.0f),
-                            2.0f, 2.0f, 2.0f, 2.0f,
-                            3.0f, 3.0f, 3.0f, 3.0f,
-                       vec4(4.0f, 4.0f, 4.0f, 4.0f));
-        assert(m4.matrix == [[1.0f, 1.0f, 1.0f, 1.0f],
-                             [2.0f, 2.0f, 2.0f, 2.0f],
-                             [3.0f, 3.0f, 3.0f, 3.0f],
-                             [4.0f, 4.0f, 4.0f, 4.0f]]);
-        assert(mat3(m4).matrix == [[1.0f, 1.0f, 1.0f],
-                                   [2.0f, 2.0f, 2.0f],
-                                   [3.0f, 3.0f, 3.0f]]);
-        assert(mat2(mat3(m4)).matrix == [[1.0f, 1.0f], [2.0f, 2.0f]]);
-        assert(mat2(m4).matrix == mat2(mat3(m4)).matrix);
-        assert(mat4(mat3(m4)).matrix == [[1.0f, 1.0f, 1.0f, 0.0f],
-                                         [2.0f, 2.0f, 2.0f, 0.0f],
-                                         [3.0f, 3.0f, 3.0f, 0.0f],
-                                         [0.0f, 0.0f, 0.0f, 1.0f]]);
-
-        Matrix!(float, 2, 3) mt1 = Matrix!(float, 2, 3)(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f);
-        Matrix!(float, 3, 2) mt2 = Matrix!(float, 3, 2)(6.0f, -1.0f, 3.0f, 2.0f, 0.0f, -3.0f);
-
-        assert(mt1.matrix == [[1.0f, 2.0f, 3.0f], [4.0f, 5.0f, 6.0f]]);
-        assert(mt2.matrix == [[6.0f, -1.0f], [3.0f, 2.0f], [0.0f, -3.0f]]);
-
-        static assert(!__traits(compiles, mat2(1, 2, 1)));
-        static assert(!__traits(compiles, mat3(1, 2, 3, 1, 2, 3, 1, 2)));
-        static assert(!__traits(compiles, mat4(1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3)));
-
-        auto m5 = mat2([0.0f,1,2,3]);
-        assert(m5.matrix == [[0.0f, 1.0f], [2.0f, 3.0f]]);
-
-        auto m6 = Matrix!(int, 2, 3)([0,1,2,3,4,5]);
-        assert(m6.matrix == [[0, 1, 2], [3, 4, 5]]);
-    }
-
-    static if(rows == cols) {
-        /// Makes the current matrix an identity matrix.
-        void make_identity() {
-            clear(0);
-            foreach(r; TupleRange!(0, rows)) {
-                matrix[r][r] = 1;
-            }
-        }
-
-        /// Returns a identity matrix.
-        static @property Matrix identity() {
-            Matrix ret;
-            ret.clear(0);
-
-            foreach(r; TupleRange!(0, rows)) {
-                ret.matrix[r][r] = 1;
-            }
-
-            return ret;
-        }
-
-        /// Transposes the current matrix;
-        void transpose() {
-            matrix = transposed().matrix;
-        }
-
-        unittest {
-            mat2 m2 = mat2(1.0f);
-            m2.transpose();
-            assert(m2.matrix == mat2(1.0f).matrix);
-            m2.make_identity();
-            assert(m2.matrix == [[1.0f, 0.0f],
-                                 [0.0f, 1.0f]]);
-            m2.transpose();
-            assert(m2.matrix == [[1.0f, 0.0f],
-                                 [0.0f, 1.0f]]);
-            assert(m2.matrix == m2.identity.matrix);
-
-            mat3 m3 = mat3(1.1f, 1.2f, 1.3f,
-                           2.1f, 2.2f, 2.3f,
-                           3.1f, 3.2f, 3.3f);
-            m3.transpose();
-            assert(m3.matrix == [[1.1f, 2.1f, 3.1f],
-                                 [1.2f, 2.2f, 3.2f],
-                                 [1.3f, 2.3f, 3.3f]]);
-
-            mat4 m4 = mat4(2.0f);
-            m4.transpose();
-            assert(m4.matrix == mat4(2.0f).matrix);
-            m4.make_identity();
-            assert(m4.matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
-                                 [0.0f, 1.0f, 0.0f, 0.0f],
-                                 [0.0f, 0.0f, 1.0f, 0.0f],
-                                 [0.0f, 0.0f, 0.0f, 1.0f]]);
-            assert(m4.matrix == m4.identity.matrix);
-        }
-
-    }
-
-    /// Returns a transposed copy of the matrix.
-    @property Matrix!(mt, cols, rows) transposed() const {
-        typeof(return) ret;
-
-        foreach(r; TupleRange!(0, rows)) {
-            foreach(c; TupleRange!(0, cols)) {
-                ret.matrix[c][r] = matrix[r][c];
-            }
-        }
-
-        return ret;
-    }
-
-    // transposed already tested in last unittest
-
-
-    static if((rows == 2) && (cols == 2)) {
-        @property mt det() const {
-            return (matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]);
-        }
-
-        private Matrix invert(ref Matrix mat) const {
-            static if(isFloatingPoint!mt && rmul) {
-                mt d = 1 / det;
-
-                mat.matrix = [[matrix[1][1]*d, -matrix[0][1]*d],
-                              [-matrix[1][0]*d, matrix[0][0]*d]];
-            } else {
-                mt d = det;
-
-                mat.matrix = [[matrix[1][1]/d, -matrix[0][1]/d],
-                              [-matrix[1][0]/d, matrix[0][0]/d]];
-            }
-
-            return mat;
-        }
-
-        static Matrix scaling(mt x, mt y) {
-            Matrix ret = Matrix.identity;
-
-            ret.matrix[0][0] = x;
-            ret.matrix[1][1] = y;
-
-            return ret;
-        }
-
-        Matrix scale(mt x, mt y) {
-            this = Matrix.scaling(x, y) * this;
-            return this;
-        }
-
-        unittest {
-            assert(mat2.scaling(3, 3).matrix == mat2.identity.scale(3, 3).matrix);
-            assert(mat2.scaling(3, 3).matrix == [[3.0f, 0.0f], [0.0f, 3.0f]]);
-        }
-
-    } else static if((rows == 3) && (cols == 3)) {
-        @property mt det() const {
-            return (matrix[0][0] * matrix[1][1] * matrix[2][2]
-                  + matrix[0][1] * matrix[1][2] * matrix[2][0]
-                  + matrix[0][2] * matrix[1][0] * matrix[2][1]
-                  - matrix[0][2] * matrix[1][1] * matrix[2][0]
-                  - matrix[0][1] * matrix[1][0] * matrix[2][2]
-                  - matrix[0][0] * matrix[1][2] * matrix[2][1]);
-        }
-
-        private Matrix invert(ref Matrix mat) const {
-            static if(isFloatingPoint!mt && rmul) {
-                mt d = 1 / det;
-                enum op = "*";
-            } else {
-                mt d = det;
-                enum op = "/";
-            }
-
-            mixin(`
-            mat.matrix = [[(matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])`~op~`d,
-                           (matrix[0][2] * matrix[2][1] - matrix[0][1] * matrix[2][2])`~op~`d,
-                           (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])`~op~`d],
-                          [(matrix[1][2] * matrix[2][0] - matrix[1][0] * matrix[2][2])`~op~`d,
-                           (matrix[0][0] * matrix[2][2] - matrix[0][2] * matrix[2][0])`~op~`d,
-                           (matrix[0][2] * matrix[1][0] - matrix[0][0] * matrix[1][2])`~op~`d],
-                          [(matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0])`~op~`d,
-                           (matrix[0][1] * matrix[2][0] - matrix[0][0] * matrix[2][1])`~op~`d,
-                           (matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0])`~op~`d]];
-            `);
-
-            return mat;
-        }
-    } else static if((rows == 4) && (cols == 4)) {
-        /// Returns the determinant of the current matrix (2x2, 3x3 and 4x4 matrices).
-        @property mt det() const {
-            return (matrix[0][3] * matrix[1][2] * matrix[2][1] * matrix[3][0] - matrix[0][2] * matrix[1][3] * matrix[2][1] * matrix[3][0]
-                  - matrix[0][3] * matrix[1][1] * matrix[2][2] * matrix[3][0] + matrix[0][1] * matrix[1][3] * matrix[2][2] * matrix[3][0]
-                  + matrix[0][2] * matrix[1][1] * matrix[2][3] * matrix[3][0] - matrix[0][1] * matrix[1][2] * matrix[2][3] * matrix[3][0]
-                  - matrix[0][3] * matrix[1][2] * matrix[2][0] * matrix[3][1] + matrix[0][2] * matrix[1][3] * matrix[2][0] * matrix[3][1]
-                  + matrix[0][3] * matrix[1][0] * matrix[2][2] * matrix[3][1] - matrix[0][0] * matrix[1][3] * matrix[2][2] * matrix[3][1]
-                  - matrix[0][2] * matrix[1][0] * matrix[2][3] * matrix[3][1] + matrix[0][0] * matrix[1][2] * matrix[2][3] * matrix[3][1]
-                  + matrix[0][3] * matrix[1][1] * matrix[2][0] * matrix[3][2] - matrix[0][1] * matrix[1][3] * matrix[2][0] * matrix[3][2]
-                  - matrix[0][3] * matrix[1][0] * matrix[2][1] * matrix[3][2] + matrix[0][0] * matrix[1][3] * matrix[2][1] * matrix[3][2]
-                  + matrix[0][1] * matrix[1][0] * matrix[2][3] * matrix[3][2] - matrix[0][0] * matrix[1][1] * matrix[2][3] * matrix[3][2]
-                  - matrix[0][2] * matrix[1][1] * matrix[2][0] * matrix[3][3] + matrix[0][1] * matrix[1][2] * matrix[2][0] * matrix[3][3]
-                  + matrix[0][2] * matrix[1][0] * matrix[2][1] * matrix[3][3] - matrix[0][0] * matrix[1][2] * matrix[2][1] * matrix[3][3]
-                  - matrix[0][1] * matrix[1][0] * matrix[2][2] * matrix[3][3] + matrix[0][0] * matrix[1][1] * matrix[2][2] * matrix[3][3]);
-        }
-
-        private Matrix invert(ref Matrix mat) const {
-            static if(isFloatingPoint!mt && rmul) {
-                mt d = 1 / det;
-                enum op = "*";
-            } else {
-                mt d = det;
-                enum op = "/";
-            }
-
-            mixin(`
-            mat.matrix = [[(matrix[1][1] * matrix[2][2] * matrix[3][3] + matrix[1][2] * matrix[2][3] * matrix[3][1] + matrix[1][3] * matrix[2][1] * matrix[3][2]
-                          - matrix[1][1] * matrix[2][3] * matrix[3][2] - matrix[1][2] * matrix[2][1] * matrix[3][3] - matrix[1][3] * matrix[2][2] * matrix[3][1])`~op~`d,
-                           (matrix[0][1] * matrix[2][3] * matrix[3][2] + matrix[0][2] * matrix[2][1] * matrix[3][3] + matrix[0][3] * matrix[2][2] * matrix[3][1]
-                          - matrix[0][1] * matrix[2][2] * matrix[3][3] - matrix[0][2] * matrix[2][3] * matrix[3][1] - matrix[0][3] * matrix[2][1] * matrix[3][2])`~op~`d,
-                           (matrix[0][1] * matrix[1][2] * matrix[3][3] + matrix[0][2] * matrix[1][3] * matrix[3][1] + matrix[0][3] * matrix[1][1] * matrix[3][2]
-                          - matrix[0][1] * matrix[1][3] * matrix[3][2] - matrix[0][2] * matrix[1][1] * matrix[3][3] - matrix[0][3] * matrix[1][2] * matrix[3][1])`~op~`d,
-                           (matrix[0][1] * matrix[1][3] * matrix[2][2] + matrix[0][2] * matrix[1][1] * matrix[2][3] + matrix[0][3] * matrix[1][2] * matrix[2][1]
-                          - matrix[0][1] * matrix[1][2] * matrix[2][3] - matrix[0][2] * matrix[1][3] * matrix[2][1] - matrix[0][3] * matrix[1][1] * matrix[2][2])`~op~`d],
-                          [(matrix[1][0] * matrix[2][3] * matrix[3][2] + matrix[1][2] * matrix[2][0] * matrix[3][3] + matrix[1][3] * matrix[2][2] * matrix[3][0]
-                          - matrix[1][0] * matrix[2][2] * matrix[3][3] - matrix[1][2] * matrix[2][3] * matrix[3][0] - matrix[1][3] * matrix[2][0] * matrix[3][2])`~op~`d,
-                           (matrix[0][0] * matrix[2][2] * matrix[3][3] + matrix[0][2] * matrix[2][3] * matrix[3][0] + matrix[0][3] * matrix[2][0] * matrix[3][2]
-                          - matrix[0][0] * matrix[2][3] * matrix[3][2] - matrix[0][2] * matrix[2][0] * matrix[3][3] - matrix[0][3] * matrix[2][2] * matrix[3][0])`~op~`d,
-                           (matrix[0][0] * matrix[1][3] * matrix[3][2] + matrix[0][2] * matrix[1][0] * matrix[3][3] + matrix[0][3] * matrix[1][2] * matrix[3][0]
-                          - matrix[0][0] * matrix[1][2] * matrix[3][3] - matrix[0][2] * matrix[1][3] * matrix[3][0] - matrix[0][3] * matrix[1][0] * matrix[3][2])`~op~`d,
-                           (matrix[0][0] * matrix[1][2] * matrix[2][3] + matrix[0][2] * matrix[1][3] * matrix[2][0] + matrix[0][3] * matrix[1][0] * matrix[2][2]
-                          - matrix[0][0] * matrix[1][3] * matrix[2][2] - matrix[0][2] * matrix[1][0] * matrix[2][3] - matrix[0][3] * matrix[1][2] * matrix[2][0])`~op~`d],
-                          [(matrix[1][0] * matrix[2][1] * matrix[3][3] + matrix[1][1] * matrix[2][3] * matrix[3][0] + matrix[1][3] * matrix[2][0] * matrix[3][1]
-                          - matrix[1][0] * matrix[2][3] * matrix[3][1] - matrix[1][1] * matrix[2][0] * matrix[3][3] - matrix[1][3] * matrix[2][1] * matrix[3][0])`~op~`d,
-                           (matrix[0][0] * matrix[2][3] * matrix[3][1] + matrix[0][1] * matrix[2][0] * matrix[3][3] + matrix[0][3] * matrix[2][1] * matrix[3][0]
-                          - matrix[0][0] * matrix[2][1] * matrix[3][3] - matrix[0][1] * matrix[2][3] * matrix[3][0] - matrix[0][3] * matrix[2][0] * matrix[3][1])`~op~`d,
-                           (matrix[0][0] * matrix[1][1] * matrix[3][3] + matrix[0][1] * matrix[1][3] * matrix[3][0] + matrix[0][3] * matrix[1][0] * matrix[3][1]
-                          - matrix[0][0] * matrix[1][3] * matrix[3][1] - matrix[0][1] * matrix[1][0] * matrix[3][3] - matrix[0][3] * matrix[1][1] * matrix[3][0])`~op~`d,
-                           (matrix[0][0] * matrix[1][3] * matrix[2][1] + matrix[0][1] * matrix[1][0] * matrix[2][3] + matrix[0][3] * matrix[1][1] * matrix[2][0]
-                          - matrix[0][0] * matrix[1][1] * matrix[2][3] - matrix[0][1] * matrix[1][3] * matrix[2][0] - matrix[0][3] * matrix[1][0] * matrix[2][1])`~op~`d],
-                          [(matrix[1][0] * matrix[2][2] * matrix[3][1] + matrix[1][1] * matrix[2][0] * matrix[3][2] + matrix[1][2] * matrix[2][1] * matrix[3][0]
-                          - matrix[1][0] * matrix[2][1] * matrix[3][2] - matrix[1][1] * matrix[2][2] * matrix[3][0] - matrix[1][2] * matrix[2][0] * matrix[3][1])`~op~`d,
-                           (matrix[0][0] * matrix[2][1] * matrix[3][2] + matrix[0][1] * matrix[2][2] * matrix[3][0] + matrix[0][2] * matrix[2][0] * matrix[3][1]
-                          - matrix[0][0] * matrix[2][2] * matrix[3][1] - matrix[0][1] * matrix[2][0] * matrix[3][2] - matrix[0][2] * matrix[2][1] * matrix[3][0])`~op~`d,
-                           (matrix[0][0] * matrix[1][2] * matrix[3][1] + matrix[0][1] * matrix[1][0] * matrix[3][2] + matrix[0][2] * matrix[1][1] * matrix[3][0]
-                          - matrix[0][0] * matrix[1][1] * matrix[3][2] - matrix[0][1] * matrix[1][2] * matrix[3][0] - matrix[0][2] * matrix[1][0] * matrix[3][1])`~op~`d,
-                           (matrix[0][0] * matrix[1][1] * matrix[2][2] + matrix[0][1] * matrix[1][2] * matrix[2][0] + matrix[0][2] * matrix[1][0] * matrix[2][1]
-                          - matrix[0][0] * matrix[1][2] * matrix[2][1] - matrix[0][1] * matrix[1][0] * matrix[2][2] - matrix[0][2] * matrix[1][1] * matrix[2][0])`~op~`d]];
-            `);
-
-            return mat;
-        }
-
-        // some static fun ...
-        // (1) glprogramming.com/red/appendixf.html - ortographic is broken!
-        // (2) http://fly.cc.fer.hr/~unreal/theredbook/appendixg.html
-        // (3) http://en.wikipedia.org/wiki/Orthographic_projection_(geometry)
-
-        static if(isFloatingPoint!mt) {
-            static private mt[6] cperspective(mt width, mt height, mt fov, mt near, mt far)
-                in { assert(height != 0); }
-                body {
-                    mt aspect = width/height;
-                    mt top = near * tan(fov*(PI/360.0));
-                    mt bottom = -top;
-                    mt right = top * aspect;
-                    mt left = -right;
-
-                    return [left, right, bottom, top, near, far];
-                }
-
-            /// Returns a perspective matrix (4x4 and floating-point matrices only).
-            static Matrix perspective(mt width, mt height, mt fov, mt near, mt far) {
-                mt[6] cdata = cperspective(width, height, fov, near, far);
-                return perspective(cdata[0], cdata[1], cdata[2], cdata[3], cdata[4], cdata[5]);
-            }
-
-            /// ditto
-            static Matrix perspective(mt left, mt right, mt bottom, mt top, mt near, mt far)
-                in {
-                    assert(right-left != 0);
-                    assert(top-bottom != 0);
-                    assert(far-near != 0);
-                }
-                body {
-                    Matrix ret;
-                    ret.clear(0);
-
-                    ret.matrix[0][0] = (2*near)/(right-left);
-                    ret.matrix[0][2] = (right+left)/(right-left);
-                    ret.matrix[1][1] = (2*near)/(top-bottom);
-                    ret.matrix[1][2] = (top+bottom)/(top-bottom);
-                    ret.matrix[2][2] = -(far+near)/(far-near);
-                    ret.matrix[2][3] = -(2*far*near)/(far-near);
-                    ret.matrix[3][2] = -1;
-
-                    return ret;
-                }
-
-            /// Returns an inverse perspective matrix (4x4 and floating-point matrices only).
-            static Matrix perspective_inverse(mt width, mt height, mt fov, mt near, mt far) {
-                mt[6] cdata = cperspective(width, height, fov, near, far);
-                return perspective_inverse(cdata[0], cdata[1], cdata[2], cdata[3], cdata[4], cdata[5]);
-            }
-
-            /// ditto
-            static Matrix perspective_inverse(mt left, mt right, mt bottom, mt top, mt near, mt far)
-                in {
-                    assert(near != 0);
-                    assert(far != 0);
-                }
-                body {
-                    Matrix ret;
-                    ret.clear(0);
-
-                    ret.matrix[0][0] = (right-left)/(2*near);
-                    ret.matrix[0][3] = (right+left)/(2*near);
-                    ret.matrix[1][1] = (top-bottom)/(2*near);
-                    ret.matrix[1][3] = (top+bottom)/(2*near);
-                    ret.matrix[2][3] = -1;
-                    ret.matrix[3][2] = -(far-near)/(2*far*near);
-                    ret.matrix[3][3] = (far+near)/(2*far*near);
-
-                    return ret;
-                }
-
-            // (2) and (3) say this one is correct
-            /// Returns an orthographic matrix (4x4 and floating-point matrices only).
-            static Matrix orthographic(mt left, mt right, mt bottom, mt top, mt near, mt far)
-                in {
-                    assert(right-left != 0);
-                    assert(top-bottom != 0);
-                    assert(far-near != 0);
-                }
-                body {
-                    Matrix ret;
-                    ret.clear(0);
-
-                    ret.matrix[0][0] = 2/(right-left);
-                    ret.matrix[0][3] = -(right+left)/(right-left);
-                    ret.matrix[1][1] = 2/(top-bottom);
-                    ret.matrix[1][3] = -(top+bottom)/(top-bottom);
-                    ret.matrix[2][2] = -2/(far-near);
-                    ret.matrix[2][3] = -(far+near)/(far-near);
-                    ret.matrix[3][3] = 1;
-
-                    return ret;
-                }
-
-            // (1) and (2) say this one is correct
-            /// Returns an inverse ortographic matrix (4x4 and floating-point matrices only).
-            static Matrix orthographic_inverse(mt left, mt right, mt bottom, mt top, mt near, mt far) {
-                Matrix ret;
-                ret.clear(0);
-
-                ret.matrix[0][0] = (right-left)/2;
-                ret.matrix[0][3] = (right+left)/2;
-                ret.matrix[1][1] = (top-bottom)/2;
-                ret.matrix[1][3] = (top+bottom)/2;
-                ret.matrix[2][2] = (far-near)/-2;
-                ret.matrix[2][3] = (far+near)/2;
-                ret.matrix[3][3] = 1;
-
-                return ret;
-            }
-
-            /// Returns a look at matrix (4x4 and floating-point matrices only).
-            static Matrix look_at(Vector!(mt, 3) eye, Vector!(mt, 3) target, Vector!(mt, 3) up) {
-                alias Vector!(mt, 3) vec3mt;
-                vec3mt look_dir = (target - eye).normalized;
-                vec3mt up_dir = up.normalized;
-
-                vec3mt right_dir = cross(look_dir, up_dir).normalized;
-                vec3mt perp_up_dir = cross(right_dir, look_dir);
-
-                Matrix ret = Matrix.identity;
-                ret.matrix[0][0..3] = right_dir.vector[];
-                ret.matrix[1][0..3] = perp_up_dir.vector[];
-                ret.matrix[2][0..3] = (-look_dir).vector[];
-
-                ret.matrix[0][3] = -dot(eye, right_dir);
-                ret.matrix[1][3] = -dot(eye, perp_up_dir);
-                ret.matrix[2][3] = dot(eye, look_dir);
-
-                return ret;
-            }
-
-            unittest {
-                mt[6] cp = cperspective(600f, 900f, 60f, 1f, 100f);
-                assert(cp[4] == 1.0f);
-                assert(cp[5] == 100.0f);
-                assert(cp[0] == -cp[1]);
-                assert((cp[0] < -0.38489f) && (cp[0] > -0.38491f));
-                assert(cp[2] == -cp[3]);
-                assert((cp[2] < -0.577349f) && (cp[2] > -0.577351f));
-
-                assert(mat4.perspective(600f, 900f, 60.0, 1.0, 100.0) == mat4.perspective(cp[0], cp[1], cp[2], cp[3], cp[4], cp[5]));
-                float[4][4] m4p = mat4.perspective(600f, 900f, 60.0, 1.0, 100.0).matrix;
-                assert((m4p[0][0] < 2.598077f) && (m4p[0][0] > 2.598075f));
-                assert(m4p[0][2] == 0.0f);
-                assert((m4p[1][1] < 1.732052) && (m4p[1][1] > 1.732050));
-                assert(m4p[1][2] == 0.0f);
-                assert((m4p[2][2] < -1.020201) && (m4p[2][2] > -1.020203));
-                assert((m4p[2][3] < -2.020201) && (m4p[2][3] > -2.020203));
-                assert((m4p[3][2] < -0.9f) && (m4p[3][2] > -1.1f));
-
-                float[4][4] m4pi = mat4.perspective_inverse(600f, 900f, 60.0, 1.0, 100.0).matrix;
-                assert((m4pi[0][0] < 0.384901) && (m4pi[0][0] > 0.384899));
-                assert(m4pi[0][3] == 0.0f);
-                assert((m4pi[1][1] < 0.577351) && (m4pi[1][1] > 0.577349));
-                assert(m4pi[1][3] == 0.0f);
-                assert(m4pi[2][3] == -1.0f);
-                assert((m4pi[3][2] < -0.494999) && (m4pi[3][2] > -0.495001));
-                assert((m4pi[3][3] < 0.505001) && (m4pi[3][3] > 0.504999));
-
-                // maybe the next tests should be improved
-                float[4][4] m4o = mat4.orthographic(-1.0f, 1.0f, -1.0f, 1.0f, -1.0f, 1.0f).matrix;
-                assert(m4o == [[1.0f, 0.0f, 0.0f, 0.0f],
-                               [0.0f, 1.0f, 0.0f, 0.0f],
-                               [0.0f, 0.0f, -1.0f, 0.0f],
-                               [0.0f, 0.0f, 0.0f, 1.0f]]);
-
-                float[4][4] m4oi = mat4.orthographic_inverse(-1.0f, 1.0f, -1.0f, 1.0f, -1.0f, 1.0f).matrix;
-                assert(m4oi == [[1.0f, 0.0f, 0.0f, 0.0f],
-                                [0.0f, 1.0f, 0.0f, 0.0f],
-                                [0.0f, 0.0f, -1.0f, 0.0f],
-                                [0.0f, 0.0f, 0.0f, 1.0f]]);
-
-                //TODO: look_at tests
-            }
-        }
-    }
-
-    static if((rows == cols) && (rows >= 3) && (rows <= 4)) {
-        /// Returns a translation matrix (3x3 and 4x4 matrices).
-        static Matrix translation(mt x, mt y, mt z) {
-            Matrix ret = Matrix.identity;
-
-            ret.matrix[0][cols-1] = x;
-            ret.matrix[1][cols-1] = y;
-            ret.matrix[2][cols-1] = z;
-
-            return ret;
-        }
-
-        /// ditto
-        static Matrix translation(Vector!(mt, 3) v) {
-            Matrix ret = Matrix.identity;
-            
-            ret.matrix[0][cols-1] = v.x;
-            ret.matrix[1][cols-1] = v.y;
-            ret.matrix[2][cols-1] = v.z;
-            
-            return ret;
-        }
-
-        /// Applys a translation on the current matrix and returns $(I this) (3x3 and 4x4 matrices).
-        Matrix translate(mt x, mt y, mt z) {
-            this = Matrix.translation(x, y, z) * this;
-            return this;
-        }
-
-        /// ditto
-        Matrix translate(Vector!(mt, 3) v) {
-            this = Matrix.translation(v) * this;
-            return this;
-        }
-
-        /// Returns a scaling matrix (3x3 and 4x4 matrices);
-        static Matrix scaling(mt x, mt y, mt z) {
-            Matrix ret = Matrix.identity;
-
-            ret.matrix[0][0] = x;
-            ret.matrix[1][1] = y;
-            ret.matrix[2][2] = z;
-
-            return ret;
-        }
-
-        /// Applys a scale to the current matrix and returns $(I this) (3x3 and 4x4 matrices).
-        Matrix scale(mt x, mt y, mt z) {
-            this = Matrix.scaling(x, y, z) * this;
-            return this;
-        }
-
-        unittest {
-            mat3 m3 = mat3.identity;
-            assert(m3.translate(1.0f, 2.0f, 3.0f).matrix == mat3.translation(1.0f, 2.0f, 3.0f).matrix);
-            assert(mat3.translation(1.0f, 2.0f, 3.0f).matrix == [[1.0f, 0.0f, 1.0f],
-                                                                 [0.0f, 1.0f, 2.0f],
-                                                                 [0.0f, 0.0f, 3.0f]]);
-            assert(mat3.identity.translate(0.0f, 1.0f, 2.0f).matrix == mat3.translation(0.0f, 1.0f, 2.0f).matrix);
-
-            mat3 m31 = mat3.identity;
-            assert(m31.translate(vec3(1.0f, 2.0f, 3.0f)).matrix == mat3.translation(vec3(1.0f, 2.0f, 3.0f)).matrix);
-            assert(mat3.translation(vec3(1.0f, 2.0f, 3.0f)).matrix == [[1.0f, 0.0f, 1.0f],
-                    [0.0f, 1.0f, 2.0f],
-                    [0.0f, 0.0f, 3.0f]]);
-            assert(mat3.identity.translate(vec3(0.0f, 1.0f, 2.0f)).matrix == mat3.translation(vec3(0.0f, 1.0f, 2.0f)).matrix);
-
-            assert(m3.scaling(0.0f, 1.0f, 2.0f).matrix == mat3.scaling(0.0f, 1.0f, 2.0f).matrix);
-            assert(mat3.scaling(0.0f, 1.0f, 2.0f).matrix == [[0.0f, 0.0f, 0.0f],
-                                                             [0.0f, 1.0f, 0.0f],
-                                                             [0.0f, 0.0f, 2.0f]]);
-            assert(mat3.identity.scale(0.0f, 1.0f, 2.0f).matrix == mat3.scaling(0.0f, 1.0f, 2.0f).matrix);
-
-            // same tests for 4x4
-
-            mat4 m4 = mat4(1.0f);
-            assert(m4.translation(1.0f, 2.0f, 3.0f).matrix == mat4.translation(1.0f, 2.0f, 3.0f).matrix);
-            assert(mat4.translation(1.0f, 2.0f, 3.0f).matrix == [[1.0f, 0.0f, 0.0f, 1.0f],
-                                                                 [0.0f, 1.0f, 0.0f, 2.0f],
-                                                                 [0.0f, 0.0f, 1.0f, 3.0f],
-                                                                 [0.0f, 0.0f, 0.0f, 1.0f]]);
-            assert(mat4.identity.translate(0.0f, 1.0f, 2.0f).matrix == mat4.translation(0.0f, 1.0f, 2.0f).matrix);
-
-            assert(m4.scaling(0.0f, 1.0f, 2.0f).matrix == mat4.scaling(0.0f, 1.0f, 2.0f).matrix);
-            assert(mat4.scaling(0.0f, 1.0f, 2.0f).matrix == [[0.0f, 0.0f, 0.0f, 0.0f],
-                                                             [0.0f, 1.0f, 0.0f, 0.0f],
-                                                             [0.0f, 0.0f, 2.0f, 0.0f],
-                                                             [0.0f, 0.0f, 0.0f, 1.0f]]);
-            assert(mat4.identity.scale(0.0f, 1.0f, 2.0f).matrix == mat4.scaling(0.0f, 1.0f, 2.0f).matrix);
-        }
-    }
-
-
-    static if((rows == cols) && (rows >= 3)) {
-        static if(isFloatingPoint!mt) {
-            /// Returns an identity matrix with an applied rotate_axis around an arbitrary axis (nxn matrices, n >= 3).
-            static Matrix rotation(real alpha, Vector!(mt, 3) axis) {
-                Matrix mult = Matrix.identity;
-
-                if(axis.length != 1) {
-                    axis.normalize();
-                }
-
-                real cosa = cos(alpha);
-                real sina = sin(alpha);
-
-                Vector!(mt, 3) temp = (1 - cosa)*axis;
-
-                mult.matrix[0][0] = to!mt(cosa + temp.x * axis.x);
-                mult.matrix[0][1] = to!mt(       temp.x * axis.y + sina * axis.z);
-                mult.matrix[0][2] = to!mt(       temp.x * axis.z - sina * axis.y);
-                mult.matrix[1][0] = to!mt(       temp.y * axis.x - sina * axis.z);
-                mult.matrix[1][1] = to!mt(cosa + temp.y * axis.y);
-                mult.matrix[1][2] = to!mt(       temp.y * axis.z + sina * axis.x);
-                mult.matrix[2][0] = to!mt(       temp.z * axis.x + sina * axis.y);
-                mult.matrix[2][1] = to!mt(       temp.z * axis.y - sina * axis.x);
-                mult.matrix[2][2] = to!mt(cosa + temp.z * axis.z);
-
-                return mult;
-            }
-
-            /// ditto
-            static Matrix rotation(real alpha, mt x, mt y, mt z) {
-                return Matrix.rotation(alpha, Vector!(mt, 3)(x, y, z));
-            }
-
-            /// Returns an identity matrix with an applied rotation around the x-axis (nxn matrices, n >= 3).
-            static Matrix xrotation(real alpha) {
-                Matrix mult = Matrix.identity;
-
-                mt cosamt = to!mt(cos(alpha));
-                mt sinamt = to!mt(sin(alpha));
-
-                mult.matrix[1][1] = cosamt;
-                mult.matrix[1][2] = -sinamt;
-                mult.matrix[2][1] = sinamt;
-                mult.matrix[2][2] = cosamt;
-
-                return mult;
-            }
-
-            /// Returns an identity matrix with an applied rotation around the y-axis (nxn matrices, n >= 3).
-            static Matrix yrotation(real alpha) {
-                Matrix mult = Matrix.identity;
-
-                mt cosamt = to!mt(cos(alpha));
-                mt sinamt = to!mt(sin(alpha));
-
-                mult.matrix[0][0] = cosamt;
-                mult.matrix[0][2] = sinamt;
-                mult.matrix[2][0] = -sinamt;
-                mult.matrix[2][2] = cosamt;
-
-                return mult;
-            }
-
-            /// Returns an identity matrix with an applied rotation around the z-axis (nxn matrices, n >= 3).
-            static Matrix zrotation(real alpha) {
-                Matrix mult = Matrix.identity;
-
-                mt cosamt = to!mt(cos(alpha));
-                mt sinamt = to!mt(sin(alpha));
-
-                mult.matrix[0][0] = cosamt;
-                mult.matrix[0][1] = -sinamt;
-                mult.matrix[1][0] = sinamt;
-                mult.matrix[1][1] = cosamt;
-
-                return mult;
-            }
-
-            Matrix rotate(real alpha, Vector!(mt, 3) axis) {
-                this = rotation(alpha, axis) * this;
-                return this;
-            }
-
-            /// Rotates the current matrix around the x-axis and returns $(I this) (nxn matrices, n >= 3).
-            Matrix rotatex(real alpha) {
-                this = xrotation(alpha) * this;
-                return this;
-            }
-
-            /// Rotates the current matrix around the y-axis and returns $(I this) (nxn matrices, n >= 3).
-            Matrix rotatey(real alpha) {
-                this = yrotation(alpha) * this;
-                return this;
-            }
-
-            /// Rotates the current matrix around the z-axis and returns $(I this) (nxn matrices, n >= 3).
-            Matrix rotatez(real alpha) {
-                this = zrotation(alpha) * this;
-                return this;
-            }
-
-            unittest {
-                assert(mat4.xrotation(0).matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
-                                                    [0.0f, 1.0f, -0.0f, 0.0f],
-                                                    [0.0f, 0.0f, 1.0f, 0.0f],
-                                                    [0.0f, 0.0f, 0.0f, 1.0f]]);
-                assert(mat4.yrotation(0).matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
-                                                    [0.0f, 1.0f, 0.0f, 0.0f],
-                                                    [0.0f, 0.0f, 1.0f, 0.0f],
-                                                    [0.0f, 0.0f, 0.0f, 1.0f]]);
-                assert(mat4.zrotation(0).matrix == [[1.0f, -0.0f, 0.0f, 0.0f],
-                                                    [0.0f, 1.0f, 0.0f, 0.0f],
-                                                    [0.0f, 0.0f, 1.0f, 0.0f],
-                                                    [0.0f, 0.0f, 0.0f, 1.0f]]);
-                mat4 xro = mat4.identity;
-                xro.rotatex(0);
-                assert(mat4.xrotation(0).matrix == xro.matrix);
-                assert(xro.matrix == mat4.identity.rotatex(0).matrix);
-                assert(xro.matrix == mat4.rotation(0, vec3(1.0f, 0.0f, 0.0f)).matrix);
-                mat4 yro = mat4.identity;
-                yro.rotatey(0);
-                assert(mat4.yrotation(0).matrix == yro.matrix);
-                assert(yro.matrix == mat4.identity.rotatey(0).matrix);
-                assert(yro.matrix == mat4.rotation(0, vec3(0.0f, 1.0f, 0.0f)).matrix);
-                mat4 zro = mat4.identity;
-                xro.rotatez(0);
-                assert(mat4.zrotation(0).matrix == zro.matrix);
-                assert(zro.matrix == mat4.identity.rotatez(0).matrix);
-                assert(zro.matrix == mat4.rotation(0, vec3(0.0f, 0.0f, 1.0f)).matrix);
-            }
-        } // isFloatingPoint
-
-
-        /// Sets the translation of the matrix (nxn matrices, n >= 3).
-        void set_translation(mt[] values...) // intended to be a property
-            in { assert(values.length >= (rows-1)); }
-            body {
-                foreach(r; TupleRange!(0, rows-1)) {
-                    matrix[r][rows-1] = values[r];
-                }
-            }
-
-        /// Copyies the translation from mat to the current matrix (nxn matrices, n >= 3).
-        void set_translation(Matrix mat) {
-            foreach(r; TupleRange!(0, rows-1)) {
-                matrix[r][rows-1] = mat.matrix[r][rows-1];
-            }
-        }
-
-        /// Returns an identity matrix with the current translation applied (nxn matrices, n >= 3)..
-        Matrix get_translation() {
-            Matrix ret = Matrix.identity;
-
-            foreach(r; TupleRange!(0, rows-1)) {
-                ret.matrix[r][rows-1] = matrix[r][rows-1];
-            }
-
-            return ret;
-        }
-
-        unittest {
-            mat3 m3 = mat3(0.0f, 1.0f, 2.0f,
-                           3.0f, 4.0f, 5.0f,
-                           6.0f, 7.0f, 1.0f);
-            assert(m3.get_translation().matrix == [[1.0f, 0.0f, 2.0f], [0.0f, 1.0f, 5.0f], [0.0f, 0.0f, 1.0f]]);
-            m3.set_translation(mat3.identity);
-            assert(mat3.identity.matrix == m3.get_translation().matrix);
-            m3.set_translation([2.0f, 5.0f]);
-            assert(m3.get_translation().matrix == [[1.0f, 0.0f, 2.0f], [0.0f, 1.0f, 5.0f], [0.0f, 0.0f, 1.0f]]);
-            assert(mat3.identity.matrix == mat3.identity.get_translation().matrix);
-
-            mat4 m4 = mat4(0.0f, 1.0f, 2.0f, 3.0f,
-                           4.0f, 5.0f, 6.0f, 7.0f,
-                           8.0f, 9.0f, 10.0f, 11.0f,
-                           12.0f, 13.0f, 14.0f, 1.0f);
-            assert(m4.get_translation().matrix == [[1.0f, 0.0f, 0.0f, 3.0f],
-                                       [0.0f, 1.0f, 0.0f, 7.0f],
-                                       [0.0f, 0.0f, 1.0f, 11.0f],
-                                       [0.0f, 0.0f, 0.0f, 1.0f]]);
-            m4.set_translation(mat4.identity);
-            assert(mat4.identity.matrix == m4.get_translation().matrix);
-            m4.set_translation([3.0f, 7.0f, 11.0f]);
-            assert(m4.get_translation().matrix == [[1.0f, 0.0f, 0.0f, 3.0f],
-                                       [0.0f, 1.0f, 0.0f, 7.0f],
-                                       [0.0f, 0.0f, 1.0f, 11.0f],
-                                       [0.0f, 0.0f, 0.0f, 1.0f]]);
-            assert(mat4.identity.matrix == mat4.identity.get_translation().matrix);
-        }
-
-        /// Sets the scale of the matrix (nxn matrices, n >= 3).
-        void set_scale(mt[] values...)
-            in { assert(values.length >= (rows-1)); }
-            body {
-                foreach(r; TupleRange!(0, rows-1)) {
-                    matrix[r][r] = values[r];
-                }
-            }
-
-        /// Copyies the scale from mat to the current matrix (nxn matrices, n >= 3).
-        void set_scale(Matrix mat) {
-            foreach(r; TupleRange!(0, rows-1)) {
-                matrix[r][r] = mat.matrix[r][r];
-            }
-        }
-
-        /// Returns an identity matrix with the current scale applied (nxn matrices, n >= 3).
-        Matrix get_scale() {
-            Matrix ret = Matrix.identity;
-
-            foreach(r; TupleRange!(0, rows-1)) {
-                ret.matrix[r][r] = matrix[r][r];
-            }
-
-            return ret;
-        }
-
-        unittest {
-            mat3 m3 = mat3(0.0f, 1.0f, 2.0f,
-                           3.0f, 4.0f, 5.0f,
-                           6.0f, 7.0f, 1.0f);
-            assert(m3.get_scale().matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 4.0f, 0.0f], [0.0f, 0.0f, 1.0f]]);
-            m3.set_scale(mat3.identity);
-            assert(mat3.identity.matrix == m3.get_scale().matrix);
-            m3.set_scale([0.0f, 4.0f]);
-            assert(m3.get_scale().matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 4.0f, 0.0f], [0.0f, 0.0f, 1.0f]]);
-            assert(mat3.identity.matrix == mat3.identity.get_scale().matrix);
-
-            mat4 m4 = mat4(0.0f, 1.0f, 2.0f, 3.0f,
-                           4.0f, 5.0f, 6.0f, 7.0f,
-                           8.0f, 9.0f, 10.0f, 11.0f,
-                           12.0f, 13.0f, 14.0f, 1.0f);
-            assert(m4.get_scale().matrix == [[0.0f, 0.0f, 0.0f, 0.0f],
-                                       [0.0f, 5.0f, 0.0f, 0.0f],
-                                       [0.0f, 0.0f, 10.0f, 0.0f],
-                                       [0.0f, 0.0f, 0.0f, 1.0f]]);
-            m4.set_scale(mat4.identity);
-            assert(mat4.identity.matrix == m4.get_scale().matrix);
-            m4.set_scale([0.0f, 5.0f, 10.0f]);
-            assert(m4.get_scale().matrix == [[0.0f, 0.0f, 0.0f, 0.0f],
-                                       [0.0f, 5.0f, 0.0f, 0.0f],
-                                       [0.0f, 0.0f, 10.0f, 0.0f],
-                                       [0.0f, 0.0f, 0.0f, 1.0f]]);
-            assert(mat4.identity.matrix == mat4.identity.get_scale().matrix);
-        }
-
-        /// Copies rot into the upper left corner, the translation (nxn matrices, n >= 3).
-        void set_rotation(Matrix!(mt, 3, 3) rot) {
-            foreach(r; TupleRange!(0, 3)) {
-                foreach(c; TupleRange!(0, 3)) {
-                    matrix[r][c] = rot[r][c];
-                }
-            }
-        }
-
-        /// Returns an identity matrix with the current rotation applied (nxn matrices, n >= 3).
-        Matrix!(mt, 3, 3) get_rotation() {
-            Matrix!(mt, 3, 3) ret = Matrix!(mt, 3, 3).identity;
-
-            foreach(r; TupleRange!(0, 3)) {
-                foreach(c; TupleRange!(0, 3)) {
-                    ret.matrix[r][c] = matrix[r][c];
-                }
-            }
-
-            return ret;
-        }
-
-        unittest {
-            mat3 m3 = mat3(0.0f, 1.0f, 2.0f,
-                           3.0f, 4.0f, 5.0f,
-                           6.0f, 7.0f, 1.0f);
-            assert(m3.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [3.0f, 4.0f, 5.0f], [6.0f, 7.0f, 1.0f]]);
-            m3.set_rotation(mat3.identity);
-            assert(mat3.identity.matrix == m3.get_rotation().matrix);
-            m3.set_rotation(mat3(0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 1.0f));
-            assert(m3.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [3.0f, 4.0f, 5.0f], [6.0f, 7.0f, 1.0f]]);
-            assert(mat3.identity.matrix == mat3.identity.get_rotation().matrix);
-
-            mat4 m4 = mat4(0.0f, 1.0f, 2.0f, 3.0f,
-                           4.0f, 5.0f, 6.0f, 7.0f,
-                           8.0f, 9.0f, 10.0f, 11.0f,
-                           12.0f, 13.0f, 14.0f, 1.0f);
-            assert(m4.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [4.0f, 5.0f, 6.0f], [8.0f, 9.0f, 10.0f]]);
-            m4.set_rotation(mat3.identity);
-            assert(mat3.identity.matrix == m4.get_rotation().matrix);
-            m4.set_rotation(mat3(0.0f, 1.0f, 2.0f, 4.0f, 5.0f, 6.0f, 8.0f, 9.0f, 10.0f));
-            assert(m4.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [4.0f, 5.0f, 6.0f], [8.0f, 9.0f, 10.0f]]);
-            assert(mat3.identity.matrix == mat4.identity.get_rotation().matrix);
-        }
-
-    }
-
-    static if((rows == cols) && (rows >= 2) && (rows <= 4)) {
-        /// Returns an inverted copy of the current matrix (nxn matrices, 2 >= n <= 4).
-        @property Matrix inverse() const {
-            Matrix mat;
-            invert(mat);
-            return mat;
-        }
-
-        /// Inverts the current matrix (nxn matrices, 2 >= n <= 4).
-        void invert() {
-            // workaround Issue #11238
-            // uses a temporary instead of invert(this)
-            Matrix temp;
-            invert(temp);
-            this.matrix = temp.matrix;
-        }
-    }
-
-    unittest {
-        mat2 m2 = mat2(1.0f, 2.0f, vec2(3.0f, 4.0f));
-        assert(m2.det == -2.0f);
-        assert(m2.inverse.matrix == [[-2.0f, 1.0f], [1.5f, -0.5f]]);
-
-        mat3 m3 = mat3(1.0f, -2.0f, 3.0f,
-                       7.0f, -1.0f, 0.0f,
-                       3.0f, 2.0f, -4.0f);
-        assert(m3.det == -1.0f);
-        assert(m3.inverse.matrix == [[-4.0f, 2.0f, -3.0f],
-                                     [-28.0f, 13.0f, -21.0f],
-                                     [-17.0f, 8.0f, -13.0f]]);
-
-        mat4 m4 = mat4(1.0f, 2.0f, 3.0f, 4.0f,
-                       -2.0f, 1.0f, 5.0f, -2.0f,
-                       2.0f, -1.0f, 7.0f, 1.0f,
-                       3.0f, -3.0f, 2.0f, 0.0f);
-        assert(m4.det == -8.0f);
-        assert(m4.inverse.matrix == [[6.875f, 7.875f, -11.75f, 11.125f],
-                                     [6.625f, 7.625f, -11.25f, 10.375f],
-                                     [-0.375f, -0.375f, 0.75f, -0.625f],
-                                     [-4.5f, -5.5f, 8.0f, -7.5f]]);
-    }
-
-    private void mms(mt inp, ref Matrix mat) const { // mat * scalar
-        for(int r = 0; r < rows; r++) {
-            for(int c = 0; c < cols; c++) {
-                mat.matrix[r][c] = matrix[r][c] * inp;
-            }
-        }
-    }
-
-    private void masm(string op)(Matrix inp, ref Matrix mat) const { // mat + or - mat
-        foreach(r; TupleRange!(0, rows)) {
-            foreach(c; TupleRange!(0, cols)) {
-                mat.matrix[r][c] = mixin("inp.matrix[r][c]" ~ op ~ "matrix[r][c]");
-            }
-        }
-    }
-
-    Matrix!(mt, rows, T.cols) opBinary(string op : "*", T)(T inp) const if(isCompatibleMatrix!T && (T.rows == cols)) {
-        Matrix!(mt, rows, T.cols) ret;
-
-        foreach(r; TupleRange!(0, rows)) {
-            foreach(c; TupleRange!(0, T.cols)) {
-                ret.matrix[r][c] = 0;
-
-                foreach(c2; TupleRange!(0, cols)) {
-                    ret.matrix[r][c] += matrix[r][c2] * inp.matrix[c2][c];
-                }
-            }
-        }
-
-        return ret;
-    }
-
-    Vector!(mt, rows) opBinary(string op : "*", T : Vector!(mt, cols))(T inp) const {
-        Vector!(mt, rows) ret;
-        ret.clear(0);
-
-        foreach(c; TupleRange!(0, cols)) {
-            foreach(r; TupleRange!(0, rows)) {
-                ret.vector[r] += matrix[r][c] * inp.vector[c];
-            }
-        }
-
-        return ret;
-    }
-
-    Matrix opBinary(string op : "*")(mt inp) const {
-        Matrix ret;
-        mms(inp, ret);
-        return ret;
-    }
-
-    Matrix opBinaryRight(string op : "*")(mt inp) const {
-        return this.opBinary!(op)(inp);
-    }
-
-    Matrix opBinary(string op)(Matrix inp) const if((op == "+") || (op == "-")) {
-        Matrix ret;
-        masm!(op)(inp, ret);
-        return ret;
-    }
-
-    void opOpAssign(string op : "*")(mt inp) {
-        mms(inp, this);
-    }
-
-    void opOpAssign(string op)(Matrix inp) if((op == "+") || (op == "-")) {
-        masm!(op)(inp, this);
-    }
-
-    void opOpAssign(string op)(Matrix inp) if(op == "*") {
-        this = this * inp;
-    }
-
-    unittest {
-        mat2 m2 = mat2(1.0f, 2.0f, 3.0f, 4.0f);
-        vec2 v2 = vec2(2.0f, 2.0f);
-        assert((m2*2).matrix == [[2.0f, 4.0f], [6.0f, 8.0f]]);
-        assert((2*m2).matrix == (m2*2).matrix);
-        m2 *= 2;
-        assert(m2.matrix == [[2.0f, 4.0f], [6.0f, 8.0f]]);
-        assert((m2*v2).vector == [12.0f, 28.0f]);
-        assert((v2*m2).vector == [16.0f, 24.0f]);
-        assert((m2*m2).matrix == [[28.0f, 40.0f], [60.0f, 88.0f]]);
-        assert((m2-m2).matrix == [[0.0f, 0.0f], [0.0f, 0.0f]]);
-        assert((m2+m2).matrix == [[4.0f, 8.0f], [12.0f, 16.0f]]);
-        m2 += m2;
-        assert(m2.matrix == [[4.0f, 8.0f], [12.0f, 16.0f]]);
-        m2 -= m2;
-        assert(m2.matrix == [[0.0f, 0.0f], [0.0f, 0.0f]]);
-
-        mat3 m3 = mat3(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f);
-        vec3 v3 = vec3(2.0f, 2.0f, 2.0f);
-        assert((m3*2).matrix == [[2.0f, 4.0f, 6.0f], [8.0f, 10.0f, 12.0f], [14.0f, 16.0f, 18.0f]]);
-        assert((2*m3).matrix == (m3*2).matrix);
-        m3 *= 2;
-        assert(m3.matrix == [[2.0f, 4.0f, 6.0f], [8.0f, 10.0f, 12.0f], [14.0f, 16.0f, 18.0f]]);
-        assert((m3*v3).vector == [24.0f, 60.0f, 96.0f]);
-        assert((v3*m3).vector == [48.0f, 60.0f, 72.0f]);
-        assert((m3*m3).matrix == [[120.0f, 144.0f, 168.0f], [264.0f, 324.0f, 384.0f], [408.0f, 504.0f, 600.0f]]);
-        assert((m3-m3).matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f]]);
-        assert((m3+m3).matrix == [[4.0f, 8.0f, 12.0f], [16.0f, 20.0f, 24.0f], [28.0f, 32.0f, 36.0f]]);
-        m3 += m3;
-        assert(m3.matrix == [[4.0f, 8.0f, 12.0f], [16.0f, 20.0f, 24.0f], [28.0f, 32.0f, 36.0f]]);
-        m3 -= m3;
-        assert(m3.matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f]]);
-
-        // test opOpAssign for matrix multiplication
-        auto m4 = mat4.translation(0,1,2);
-        m4 *= mat4.translation(0,-1,2);
-        assert(m4 == mat4.translation(0,0,4));
-
-        //TODO: tests for mat4, mat34
-    }
-
-    // opEqual => "alias matrix this;"
-
-    bool opCast(T : bool)() const {
-        return isFinite;
-    }
-
-    unittest {
-        assert(mat2(1.0f, 2.0f, 1.0f, 1.0f) == mat2(1.0f, 2.0f, 1.0f, 1.0f));
-        assert(mat2(1.0f, 2.0f, 1.0f, 1.0f) != mat2(1.0f, 1.0f, 1.0f, 1.0f));
-
-        assert(mat3(1.0f) == mat3(1.0f));
-        assert(mat3(1.0f) != mat3(2.0f));
-
-        assert(mat4(1.0f) == mat4(1.0f));
-        assert(mat4(1.0f) != mat4(2.0f));
-
-        assert(!(mat4(float.nan)));
-        if(mat4(1.0f)) { }
-        else { assert(false); }
-    }
-
-}
-
-/// Pre-defined matrix types, the first number represents the number of rows
-/// and the second the number of columns, if there's just one it's a nxn matrix.
-/// All of these matrices are floating-point matrices.
-alias Matrix!(float, 2, 2) mat2;
-alias Matrix!(float, 3, 3) mat3;
-alias Matrix!(float, 3, 4) mat34;
-alias Matrix!(float, 4, 4) mat4;
-
-private unittest {
-    Matrix!(float,  1, 1) A = 1;
-    Matrix!(double, 1, 1) B = 1;
-    Matrix!(real,   1, 1) C = 1;
-    Matrix!(int,    1, 1) D = 1;
-    Matrix!(float,  5, 1) E = 1;
-    Matrix!(double, 5, 1) F = 1;
-    Matrix!(real,   5, 1) G = 1;
-    Matrix!(int,    5, 1) H = 1;
-    Matrix!(float,  1, 5) I = 1;
-    Matrix!(double, 1, 5) J = 1;
-    Matrix!(real,   1, 5) K = 1;
-    Matrix!(int,    1, 5) L = 1;
-}
-
-/// Base template for all quaternion-types.
-/// Params:
-///  type = all values get stored as this type
-struct Quaternion(type) {
-    alias type qt; /// Holds the internal type of the quaternion.
-
-    qt[4] quaternion; /// Holds the w, x, y and z coordinates.
-
-    /// Returns a pointer to the quaternion in memory, it starts with the w coordinate.
-    @property auto value_ptr() const { return quaternion.ptr; }
-
-    /// Returns the current vector formatted as string, useful for printing the quaternion.
-    @property string as_string() {
-        return format("%s", quaternion);
-    }
-    alias as_string toString;
-
-    @safe pure nothrow:
-    private @property qt get_(char coord)() const {
-        return quaternion[coord_to_index!coord];
-    }
-    private @property void set_(char coord)(qt value) {
-        quaternion[coord_to_index!coord] = value;
-    }
-
-    alias get_!'w' w; /// static properties to access the values.
-    alias set_!'w' w;
-    alias get_!'x' x; /// ditto
-    alias set_!'x' x;
-    alias get_!'y' y; /// ditto
-    alias set_!'y' y;
-    alias get_!'z' z; /// ditto
-    alias set_!'z' z;
-
-    /// Constructs the quaternion.
-    /// Takes a 4-dimensional vector, where vector.x = the quaternions w coordinate,
-    /// or a w coordinate of type $(I qt) and a 3-dimensional vector representing the imaginary part,
-    /// or 4 values of type $(I qt).
-    this(qt w_, qt x_, qt y_, qt z_) {
-        w = w_;
-        x = x_;
-        y = y_;
-        z = z_;
-    }
-
-    /// ditto
-    this(qt w_, Vector!(qt, 3) vec) {
-        w = w_;
-        quaternion[1..4] = vec.vector[];
-    }
-
-    /// ditto
-    this(Vector!(qt, 4) vec) {
-        quaternion[] = vec.vector[];
-    }
-
-    /// Returns true if all values are not nan and finite, otherwise false.
-    @property bool isFinite() const {
-        foreach(q; quaternion) {
-            if(isNaN(q) || isInfinity(q)) {
-                return false;
-            }
-        }
-        return true;
-    }
-    deprecated("Use isFinite instead of ok") alias ok = isFinite;
-
-    unittest {
-        quat q1 = quat(0.0f, 0.0f, 0.0f, 1.0f);
-        assert(q1.quaternion == [0.0f, 0.0f, 0.0f, 1.0f]);
-        assert(q1.quaternion == quat(0.0f, 0.0f, 0.0f, 1.0f).quaternion);
-        assert(q1.quaternion == quat(0.0f, vec3(0.0f, 0.0f, 1.0f)).quaternion);
-        assert(q1.quaternion == quat(vec4(0.0f, 0.0f, 0.0f, 1.0f)).quaternion);
-
-        assert(q1.isFinite);
-        q1.x = float.infinity;
-        assert(!q1.isFinite);
-        q1.x = float.nan;
-        assert(!q1.isFinite);
-        q1.x = 0.0f;
-        assert(q1.isFinite);
-    }
-
-    template coord_to_index(char c) {
-        static if(c == 'w') {
-            enum coord_to_index = 0;
-        } else static if(c == 'x') {
-            enum coord_to_index = 1;
-        } else static if(c == 'y') {
-            enum coord_to_index = 2;
-        } else static if(c == 'z') {
-            enum coord_to_index = 3;
-        } else {
-            static assert(false, "accepted coordinates are x, y, z and w not " ~ c ~ ".");
-        }
-    }
-
-    /// Returns the squared magnitude of the quaternion.
-    @property real magnitude_squared() const {
-        return to!real(w^^2 + x^^2 + y^^2 + z^^2);
-    }
-
-    /// Returns the magnitude of the quaternion.
-    @property real magnitude() const {
-        return sqrt(magnitude_squared);
-    }
-
-    /// Returns an identity quaternion (w=1, x=0, y=0, z=0).
-    static @property Quaternion identity() {
-        return Quaternion(1, 0, 0, 0);
-    }
-
-    /// Makes the current quaternion an identity quaternion.
-    void make_identity() {
-        w = 1;
-        x = 0;
-        y = 0;
-        z = 0;
-    }
-
-    /// Inverts the quaternion.
-    void invert() {
-        x = -x;
-        y = -y;
-        z = -z;
-    }
-    alias invert conjugate; /// ditto
-
-    /// Returns an inverted copy of the current quaternion.
-    @property Quaternion inverse() const {
-        return Quaternion(w, -x, -y, -z);
-    }
-    alias inverse conjugated; /// ditto
-
-    unittest {
-        quat q1 = quat(1.0f, 1.0f, 1.0f, 1.0f);
-
-        assert(q1.magnitude == 2.0f);
-        assert(q1.magnitude_squared == 4.0f);
-        assert(q1.magnitude == quat(0.0f, 0.0f, 2.0f, 0.0f).magnitude);
-
-        quat q2 = quat.identity;
-        assert(q2.quaternion == [1.0f, 0.0f, 0.0f, 0.0f]);
-        assert(q2.x == 0.0f);
-        assert(q2.y == 0.0f);
-        assert(q2.z == 0.0f);
-        assert(q2.w == 1.0f);
-
-        assert(q1.inverse.quaternion == [1.0f, -1.0f, -1.0f, -1.0f]);
-        q1.invert();
-        assert(q1.quaternion == [1.0f, -1.0f, -1.0f, -1.0f]);
-
-        q1.make_identity();
-        assert(q1.quaternion == q2.quaternion);
-
-    }
-
-    /// Creates a quaternion from a 3x3 matrix.
-    /// Params:
-    ///  matrix = 3x3 matrix (rotation)
-    /// Returns: A quaternion representing the rotation (3x3 matrix)
-    static Quaternion from_matrix(Matrix!(qt, 3, 3) matrix) {
-        Quaternion ret;
-
-        auto mat = matrix.matrix;
-        qt trace = mat[0][0] + mat[1][1] + mat[2][2];
-
-        if(trace > 0) {
-            real s = 0.5 / sqrt(trace + 1.0f);
-
-            ret.w = to!qt(0.25 / s);
-            ret.x = to!qt((mat[2][1] - mat[1][2]) * s);
-            ret.y = to!qt((mat[0][2] - mat[2][0]) * s);
-            ret.z = to!qt((mat[1][0] - mat[0][1]) * s);
-        } else if((mat[0][0] > mat[1][1]) && (mat[0][0] > mat[2][2])) {
-            real s = 2.0 * sqrt(1.0 + mat[0][0] - mat[1][1] - mat[2][2]);
-
-            ret.w = to!qt((mat[2][1] - mat[1][2]) / s);
-            ret.x = to!qt(0.25f * s);
-            ret.y = to!qt((mat[0][1] + mat[1][0]) / s);
-            ret.z = to!qt((mat[0][2] + mat[2][0]) / s);
-        } else if(mat[1][1] > mat[2][2]) {
-            real s = 2.0 * sqrt(1 + mat[1][1] - mat[0][0] - mat[2][2]);
-
-            ret.w = to!qt((mat[0][2] - mat[2][0]) / s);
-            ret.x = to!qt((mat[0][1] + mat[1][0]) / s);
-            ret.y = to!qt(0.25f * s);
-            ret.z = to!qt((mat[1][2] + mat[2][1]) / s);
-        } else {
-            real s = 2.0 * sqrt(1 + mat[2][2] - mat[0][0] - mat[1][1]);
-
-            ret.w = to!qt((mat[1][0] - mat[0][1]) / s);
-            ret.x = to!qt((mat[0][2] + mat[2][0]) / s);
-            ret.y = to!qt((mat[1][2] + mat[2][1]) / s);
-            ret.z = to!qt(0.25f * s);
-        }
-
-        return ret;
-    }
-
-    /// Returns the quaternion as matrix.
-    /// Params:
-    ///  rows = number of rows of the resulting matrix (min 3)
-    ///  cols = number of columns of the resulting matrix (min 3)
-    Matrix!(qt, rows, cols) to_matrix(int rows, int cols)() const if((rows >= 3) && (cols >= 3)) {
-        static if((rows == 3) && (cols == 3)) {
-            Matrix!(qt, rows, cols) ret;
-        } else {
-            Matrix!(qt, rows, cols) ret = Matrix!(qt, rows, cols).identity;
-        }
-
-        qt xx = x^^2;
-        qt xy = x * y;
-        qt xz = x * z;
-        qt xw = x * w;
-        qt yy = y^^2;
-        qt yz = y * z;
-        qt yw = y * w;
-        qt zz = z^^2;
-        qt zw = z * w;
-
-        ret.matrix[0][0] = 1 - 2 * (yy + zz);
-        ret.matrix[0][1] = 2 * (xy - zw);
-        ret.matrix[0][2] = 2 * (xz + yw);
-
-        ret.matrix[1][0] = 2 * (xy + zw);
-        ret.matrix[1][1] = 1 - 2 * (xx + zz);
-        ret.matrix[1][2] = 2 * (yz - xw);
-
-        ret.matrix[2][0] = 2 * (xz - yw);
-        ret.matrix[2][1] = 2 * (yz + xw);
-        ret.matrix[2][2] = 1 - 2 * (xx + yy);
-
-        return ret;
-    }
-
-    unittest {
-        quat q1 = quat(4.0f, 1.0f, 2.0f, 3.0f);
-
-        assert(q1.to_matrix!(3, 3).matrix == [[-25.0f, -20.0f, 22.0f], [28.0f, -19.0f, 4.0f], [-10.0f, 20.0f, -9.0f]]);
-        assert(q1.to_matrix!(4, 4).matrix == [[-25.0f, -20.0f, 22.0f, 0.0f],
-                                              [28.0f, -19.0f, 4.0f, 0.0f],
-                                              [-10.0f, 20.0f, -9.0f, 0.0f],
-                                              [0.0f, 0.0f, 0.0f, 1.0f]]);
-        assert(quat.identity.to_matrix!(3, 3).matrix == Matrix!(qt, 3, 3).identity.matrix);
-        assert(q1.quaternion == quat.from_matrix(q1.to_matrix!(3, 3)).quaternion);
-
-        assert(quat(1.0f, 0.0f, 0.0f, 0.0f).quaternion == quat.from_matrix(mat3.identity).quaternion);
-
-        quat q2 = quat.from_matrix(mat3(1.0f, 3.0f, 2.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f));
-        assert(q2.x == 0.0f);
-        assert(almost_equal(q2.y, 0.7071067f));
-        assert(almost_equal(q2.z, -1.060660));
-        assert(almost_equal(q2.w, 0.7071067f));
-    }
-
-    /// Normalizes the current quaternion.
-    void normalize() {
-        qt m = to!qt(magnitude);
-
-        if(m != 0) {
-            w = w / m;
-            x = x / m;
-            y = y / m;
-            z = z / m;
-        }
-    }
-
-    /// Returns a normalized copy of the current quaternion.
-    Quaternion normalized() const {
-        Quaternion ret;
-        qt m = to!qt(magnitude);
-
-        if(m != 0) {
-            ret.w = w / m;
-            ret.x = x / m;
-            ret.y = y / m;
-            ret.z = z / m;
-        } else {
-            ret = Quaternion(w, x, y, z);
-        }
-
-        return ret;
-    }
-
-    unittest {
-        quat q1 = quat(1.0f, 2.0f, 3.0f, 4.0f);
-        quat q2 = quat(1.0f, 2.0f, 3.0f, 4.0f);
-
-        q1.normalize();
-        assert(q1.quaternion == q2.normalized.quaternion);
-        //assert(q1.quaternion == q1.normalized.quaternion);
-        assert(almost_equal(q1.magnitude, 1.0));
-    }
-
-    /// Returns the yaw.
-    @property real yaw() const {
-        return atan2(to!real(2.0*(w*z + x*y)), to!real(1.0 - 2.0*(y*y + z*z)));
-    }
-
-    /// Returns the pitch.
-    @property real pitch() const {
-        return asin(to!real(2.0*(w*y - z*x)));
-    }
-
-    /// Returns the roll.
-    @property real roll() const {
-        return atan2(to!real(2.0*(w*x + y*z)), to!real(1.0 - 2.0*(x*x + y*y)));
-    }
-
-    unittest {
-        quat q1 = quat.identity;
-        assert(q1.pitch == 0.0f);
-        assert(q1.yaw == 0.0f);
-        assert(q1.roll == 0.0f);
-
-        quat q2 = quat(1.0f, 1.0f, 1.0f, 1.0f);
-        assert(almost_equal(q2.yaw, q2.roll));
-        //assert(almost_equal(q2.yaw, 1.570796f));
-        assert(q2.pitch == 0.0f);
-
-        quat q3 = quat(0.1f, 1.9f, 2.1f, 1.3f);
-        //assert(almost_equal(q3.yaw, 2.4382f));
-        assert(isNaN(q3.pitch));
-        //assert(almost_equal(q3.roll, 1.67719f));
-    }
-
-    /// Returns a quaternion with applied rotation around the x-axis.
-    static Quaternion xrotation(real alpha) {
-        Quaternion ret;
-
-        alpha /= 2;
-        ret.w = to!qt(cos(alpha));
-        ret.x = to!qt(sin(alpha));
-        ret.y = 0;
-        ret.z = 0;
-
-        return ret;
-    }
-
-    /// Returns a quaternion with applied rotation around the y-axis.
-    static Quaternion yrotation(real alpha) {
-        Quaternion ret;
-
-        alpha /= 2;
-        ret.w = to!qt(cos(alpha));
-        ret.x = 0;
-        ret.y = to!qt(sin(alpha));
-        ret.z = 0;
-
-        return ret;
-    }
-
-    /// Returns a quaternion with applied rotation around the z-axis.
-    static Quaternion zrotation(real alpha) {
-        Quaternion ret;
-
-        alpha /= 2;
-        ret.w = to!qt(cos(alpha));
-        ret.x = 0;
-        ret.y = 0;
-        ret.z = to!qt(sin(alpha));
-
-        return ret;
-    }
-
-    /// Returns a quaternion with applied rotation around an axis.
-    static Quaternion axis_rotation(real alpha, Vector!(qt, 3) axis) {
-        if(alpha == 0) {
-            return Quaternion.identity;
-        }
-        Quaternion ret;
-
-        alpha /= 2;
-        qt sinaqt = to!qt(sin(alpha));
-
-        ret.w = to!qt(cos(alpha));
-        ret.x = axis.x * sinaqt;
-        ret.y = axis.y * sinaqt;
-        ret.z = axis.z * sinaqt;
-
-        return ret;
-    }
-
-    /// Creates a quaternion from an euler rotation.
-    static Quaternion euler_rotation(real roll, real pitch, real yaw) {
-        Quaternion ret;
-
-        auto cr = cos(roll / 2.0);
-        auto cp = cos(pitch / 2.0);
-        auto cy = cos(yaw / 2.0);
-        auto sr = sin(roll / 2.0);
-        auto sp = sin(pitch / 2.0);
-        auto sy = sin(yaw / 2.0);
-
-        ret.w = cr * cp * cy + sr * sp * sy;
-        ret.x = sr * cp * cy - cr * sp * sy;
-        ret.y = cr * sp * cy + sr * cp * sy;
-        ret.z = cr * cp * sy - sr * sp * cy;
-
-        return ret;
-    }
-
-    unittest {
-        enum startPitch = 0.1;
-        enum startYaw = -0.2;
-        enum startRoll = 0.6;
-        
-        auto q = quat.euler_rotation(startRoll,startPitch,startYaw);
-        
-        assert(almost_equal(q.pitch,startPitch));
-        assert(almost_equal(q.yaw,startYaw));
-        assert(almost_equal(q.roll,startRoll));
-    }
-
-    /// Rotates the current quaternion around the x-axis and returns $(I this).
-    Quaternion rotatex(real alpha) {
-        this = xrotation(alpha) * this;
-        return this;
-    }
-
-    /// Rotates the current quaternion around the y-axis and returns $(I this).
-    Quaternion rotatey(real alpha) {
-        this = yrotation(alpha) * this;
-        return this;
-    }
-
-    /// Rotates the current quaternion around the z-axis and returns $(I this).
-    Quaternion rotatez(real alpha) {
-        this = zrotation(alpha) * this;
-        return this;
-    }
-
-    /// Rotates the current quaternion around an axis and returns $(I this).
-    Quaternion rotate_axis(real alpha, Vector!(qt, 3) axis) {
-        this = axis_rotation(alpha, axis) * this;
-        return this;
-    }
-
-    /// Applies an euler rotation to the current quaternion and returns $(I this).
-    Quaternion rotate_euler(real heading, real attitude, real bank) {
-        this = euler_rotation(heading, attitude, bank) * this;
-        return this;
-    }
-
-    unittest {
-        assert(quat.xrotation(PI).quaternion[1..4] == [1.0f, 0.0f, 0.0f]);
-        assert(quat.yrotation(PI).quaternion[1..4] == [0.0f, 1.0f, 0.0f]);
-        assert(quat.zrotation(PI).quaternion[1..4] == [0.0f, 0.0f, 1.0f]);
-        assert((quat.xrotation(PI).w == quat.yrotation(PI).w) && (quat.yrotation(PI).w == quat.zrotation(PI).w));
-        //assert(quat.rotatex(PI).w == to!(quat.qt)(cos(PI)));
-        assert(quat.xrotation(PI).quaternion == quat.identity.rotatex(PI).quaternion);
-        assert(quat.yrotation(PI).quaternion == quat.identity.rotatey(PI).quaternion);
-        assert(quat.zrotation(PI).quaternion == quat.identity.rotatez(PI).quaternion);
-
-        assert(quat.axis_rotation(PI, vec3(1.0f, 1.0f, 1.0f)).quaternion[1..4] == [1.0f, 1.0f, 1.0f]);
-        assert(quat.axis_rotation(PI, vec3(1.0f, 1.0f, 1.0f)).w == quat.xrotation(PI).w);
-        assert(quat.axis_rotation(PI, vec3(1.0f, 1.0f, 1.0f)).quaternion ==
-               quat.identity.rotate_axis(PI, vec3(1.0f, 1.0f, 1.0f)).quaternion);
-
-        quat q1 = quat.euler_rotation(PI, PI, PI);
-        assert((q1.x > -1e-16) && (q1.x < 1e-16));
-        assert((q1.y > -1e-16) && (q1.y < 1e-16));
-        assert((q1.z > -1e-16) && (q1.z < 1e-16));
-        //assert(q1.w == -1.0f);
-        assert(quat.euler_rotation(PI, PI, PI).quaternion == quat.identity.rotate_euler(PI, PI, PI).quaternion);
-    }
-
-    Quaternion opBinary(string op : "*")(Quaternion inp) const {
-        Quaternion ret;
-
-        ret.w = -x * inp.x - y * inp.y - z * inp.z + w * inp.w;
-        ret.x = x * inp.w + y * inp.z - z * inp.y + w * inp.x;
-        ret.y = -x * inp.z + y * inp.w + z * inp.x + w * inp.y;
-        ret.z = x * inp.y - y * inp.x + z * inp.w + w * inp.z;
-
-        return ret;
-    }
-
-    auto opBinaryRight(string op, T)(T inp) const if(!is_quaternion!T) {
-        return this.opBinary!(op)(inp);
-    }
-
-    Quaternion opBinary(string op)(Quaternion inp) const  if((op == "+") || (op == "-")) {
-        Quaternion ret;
-
-        mixin("ret.w = w" ~ op ~ "inp.w;");
-        mixin("ret.x = x" ~ op ~ "inp.x;");
-        mixin("ret.y = y" ~ op ~ "inp.y;");
-        mixin("ret.z = z" ~ op ~ "inp.z;");
-
-        return ret;
-    }
-
-    Vector!(qt, 3) opBinary(string op : "*")(Vector!(qt, 3) inp) const {
-        Vector!(qt, 3) ret;
-
-        qt ww = w^^2;
-        qt w2 = w * 2;
-        qt wx2 = w2 * x;
-        qt wy2 = w2 * y;
-        qt wz2 = w2 * z;
-        qt xx = x^^2;
-        qt x2 = x * 2;
-        qt xy2 = x2 * y;
-        qt xz2 = x2 * z;
-        qt yy = y^^2;
-        qt yz2 = 2 * y * z;
-        qt zz = z * z;
-
-        ret.vector =  [ww * inp.x + wy2 * inp.z - wz2 * inp.y + xx * inp.x +
-                       xy2 * inp.y + xz2 * inp.z - zz * inp.x - yy * inp.x,
-                       xy2 * inp.x + yy * inp.y + yz2 * inp.z + wz2 * inp.x -
-                       zz * inp.y + ww * inp.y - wx2 * inp.z - xx * inp.y,
-                       xz2 * inp.x + yz2 * inp.y + zz * inp.z - wy2 * inp.x -
-                       yy * inp.z + wx2 * inp.y - xx * inp.z + ww * inp.z];
-
-       return ret;
-    }
-
-    Quaternion opBinary(string op : "*")(qt inp) const {
-        return Quaternion(w*inp, x*inp, y*inp, z*inp);
-    }
-
-    void opOpAssign(string op : "*")(Quaternion inp) {
-        qt w2 = -x * inp.x - y * inp.y - z * inp.z + w * inp.w;
-        qt x2 = x * inp.w + y * inp.z - z * inp.y + w * inp.x;
-        qt y2 = -x * inp.z + y * inp.w + z * inp.x + w * inp.y;
-        qt z2 = x * inp.y - y * inp.x + z * inp.w + w * inp.z;
-        w = w2; x = x2; y = y2; z = z2;
-    }
-
-    void opOpAssign(string op)(Quaternion inp) if((op == "+") || (op == "-")) {
-        mixin("w = w" ~ op ~ "inp.w;");
-        mixin("x = x" ~ op ~ "inp.x;");
-        mixin("y = y" ~ op ~ "inp.y;");
-        mixin("z = z" ~ op ~ "inp.z;");
-    }
-
-    void opOpAssign(string op : "*")(qt inp) {
-        quaternion[0] *= inp;
-        quaternion[1] *= inp;
-        quaternion[2] *= inp;
-        quaternion[3] *= inp;
-    }
-
-    unittest {
-        quat q1 = quat.identity;
-        quat q2 = quat(3.0f, 0.0f, 1.0f, 2.0f);
-        quat q3 = quat(3.4f, 0.1f, 1.2f, 2.3f);
-
-        assert((q1 * q1).quaternion == q1.quaternion);
-        assert((q1 * q2).quaternion == q2.quaternion);
-        assert((q2 * q1).quaternion == q2.quaternion);
-        quat q4 = q3 * q2;
-        assert((q2 * q3).quaternion != q4.quaternion);
-        q3 *= q2;
-        assert(q4.quaternion == q3.quaternion);
-        assert(almost_equal(q4.x, 0.4f));
-        assert(almost_equal(q4.y, 6.8f));
-        assert(almost_equal(q4.z, 13.8f));
-        assert(almost_equal(q4.w, 4.4f));
-
-        quat q5 = quat(1.0f, 2.0f, 3.0f, 4.0f);
-        quat q6 = quat(3.0f, 1.0f, 6.0f, 2.0f);
-
-        assert((q5 - q6).quaternion == [-2.0f, 1.0f, -3.0f, 2.0f]);
-        assert((q5 + q6).quaternion == [4.0f, 3.0f, 9.0f, 6.0f]);
-        assert((q6 - q5).quaternion == [2.0f, -1.0f, 3.0f, -2.0f]);
-        assert((q6 + q5).quaternion == [4.0f, 3.0f, 9.0f, 6.0f]);
-        q5 += q6;
-        assert(q5.quaternion == [4.0f, 3.0f, 9.0f, 6.0f]);
-        q6 -= q6;
-        assert(q6.quaternion == [0.0f, 0.0f, 0.0f, 0.0f]);
-
-        quat q7 = quat(2.0f, 2.0f, 2.0f, 2.0f);
-        assert((q7 * 2).quaternion == [4.0f, 4.0f, 4.0f, 4.0f]);
-        assert((2 * q7).quaternion == (q7 * 2).quaternion);
-        q7 *= 2;
-        assert(q7.quaternion == [4.0f, 4.0f, 4.0f, 4.0f]);
-
-        vec3 v1 = vec3(1.0f, 2.0f, 3.0f);
-        assert((q1 * v1).vector == v1.vector);
-        assert((v1 * q1).vector == (q1 * v1).vector);
-        assert((q2 * v1).vector == [-2.0f, 36.0f, 38.0f]);
-    }
-
-    int opCmp(ref const Quaternion qua) const {
-        foreach(i, a; quaternion) {
-            if(a < qua.quaternion[i]) {
-                return -1;
-            } else if(a > qua.quaternion[i]) {
-                return 1;
-            }
-        }
-
-        // Quaternions are the same
-        return 0;
-    }
-
-    bool opEquals(const Quaternion qu) const {
-        return quaternion == qu.quaternion;
-    }
-
-    bool opCast(T : bool)() const  {
-        return isFinite;
-    }
-
-    unittest {
-        assert(quat(1.0f, 2.0f, 3.0f, 4.0f) == quat(1.0f, 2.0f, 3.0f, 4.0f));
-        assert(quat(1.0f, 2.0f, 3.0f, 4.0f) != quat(1.0f, 2.0f, 3.0f, 3.0f));
-
-        assert(!(quat(float.nan, float.nan, float.nan, float.nan)));
-        if(quat(1.0f, 1.0f, 1.0f, 1.0f)) { }
-        else { assert(false); }
-    }
-
-}
-
-/// Pre-defined quaternion of type float.
-alias Quaternion!(float) quat;
diff --git a/gl3n/linalg/matrix.d b/gl3n/linalg/matrix.d
new file mode 100644
index 0000000..0939fd5
--- /dev/null
+++ b/gl3n/linalg/matrix.d
@@ -0,0 +1,1271 @@
+/**
+gl3n.linalg.matrix
+
+Special thanks to:
+$(UL
+  $(LI Tomasz Stachowiak (h3r3tic): allowed me to use parts of $(LINK2 https://bitbucket.org/h3r3tic/boxen/src/default/src/xf/omg, omg).)
+  $(LI Jakob Øvrum (jA_cOp): improved the code a lot!)
+  $(LI Florian Boesch (___doc__): helps me to understand opengl/complex maths better, see: $(LINK http://codeflow.org/).)
+  $(LI #D on freenode: answered general questions about D.)
+)
+
+Authors: David Herberth, Stephan Dilly
+License: MIT
+
+Note: All methods marked with pure are weakly pure since, they all access an instance member.
+All static methods are strongly pure.
+*/
+
+module gl3n.linalg.matrix;
+
+import gl3n.linalg.vec;
+import gl3n.math : PI, sin, cos, tan;
+
+import std.math : isNaN, isInfinity;
+import std.traits : isIntegral;
+import std.string : format, rightJustify;
+import std.array : join;
+import std.algorithm : min, reduce;
+import std.traits : isFloatingPoint;
+
+version(NoReciprocalMul) {
+    private enum rmul = false;
+} else {
+    private enum rmul = true;
+}
+
+/// Base template for all matrix-types.
+/// Params:
+///  type = all values get stored as this type
+///  rows_ = rows of the matrix
+///  cols_ = columns of the matrix
+/// Examples:
+/// ---
+/// alias Matrix!(float, 4, 4) mat4;
+/// alias Matrix!(double, 3, 4) mat34d;
+/// alias Matrix!(real, 2, 2) mat2r;
+/// ---
+struct Matrix(type, int rows_, int cols_) if((rows_ > 0) && (cols_ > 0)) {
+    alias type mt; /// Holds the internal type of the matrix;
+    static const int rows = rows_; /// Holds the number of rows;
+    static const int cols = cols_; /// Holds the number of columns;
+    
+    /// Holds the matrix $(RED row-major) in memory.
+    mt[cols][rows] matrix; // In C it would be mt[rows][cols], D does it like this: (mt[foo])[bar]
+    alias matrix this;
+    
+    unittest {
+        mat2 m2 = mat2(0.0f, 1.0f, 2.0f, 3.0f);
+        assert(m2[0][0] == 0.0f);
+        assert(m2[0][1] == 1.0f);
+        assert(m2[1][0] == 2.0f);
+        assert(m2[1][1] == 3.0f);
+        m2[0..1] = [2.0f, 2.0f];
+        assert(m2 == [[2.0f, 2.0f], [2.0f, 3.0f]]);
+        
+        mat3 m3 = mat3(0.0f, 0.1f, 0.2f, 1.0f, 1.1f, 1.2f, 2.0f, 2.1f, 2.2f);
+        assert(m3[0][1] == 0.1f);
+        assert(m3[2][0] == 2.0f);
+        assert(m3[1][2] == 1.2f);
+        m3[0][0..$] = 0.0f;
+        assert(m3 == [[0.0f, 0.0f, 0.0f],
+                [1.0f, 1.1f, 1.2f],
+                [2.0f, 2.1f, 2.2f]]);
+        
+        mat4 m4 = mat4(0.0f, 0.1f, 0.2f, 0.3f,
+            1.0f, 1.1f, 1.2f, 1.3f,
+            2.0f, 2.1f, 2.2f, 2.3f,
+            3.0f, 3.1f, 3.2f, 3.3f);
+        assert(m4[0][3] == 0.3f);
+        assert(m4[1][1] == 1.1f);
+        assert(m4[2][0] == 2.0f);
+        assert(m4[3][2] == 3.2f);
+        m4[2][1..3] = [1.0f, 2.0f];
+        assert(m4 == [[0.0f, 0.1f, 0.2f, 0.3f],
+                [1.0f, 1.1f, 1.2f, 1.3f],
+                [2.0f, 1.0f, 2.0f, 2.3f],
+                [3.0f, 3.1f, 3.2f, 3.3f]]);
+        
+    }
+    
+    /// Returns the pointer to the stored values as OpenGL requires it.
+    /// Note this will return a pointer to a $(RED row-major) matrix,
+    /// $(RED this means you've to set the transpose argument to GL_TRUE when passing it to OpenGL).
+    /// Examples:
+    /// ---
+    /// // 3rd argument = GL_TRUE
+    /// glUniformMatrix4fv(programs.main.model, 1, GL_TRUE, mat4.translation(-0.5f, -0.5f, 1.0f).value_ptr);
+    /// ---
+    @property auto value_ptr() const { return matrix[0].ptr; }
+    
+    /// Returns the current matrix formatted as flat string.
+    @property string as_string() {
+        return format("%s", matrix);
+    }
+    alias as_string toString; /// ditto
+    
+    /// Returns the current matrix as pretty formatted string.
+    @property string as_pretty_string() {
+        string fmtr = "%s";
+        
+        size_t rjust = max(format(fmtr, reduce!(max)(matrix[])).length,
+            format(fmtr, reduce!(min)(matrix[])).length) - 1;
+        
+        string[] outer_parts;
+        foreach(mt[] row; matrix) {
+            string[] inner_parts;
+            foreach(mt col; row) {
+                inner_parts ~= rightJustify(format(fmtr, col), rjust);
+            }
+            outer_parts ~= " [" ~ join(inner_parts, ", ") ~ "]";
+        }
+        
+        return "[" ~ join(outer_parts, "\n")[1..$] ~ "]";
+    }
+    alias as_pretty_string toPrettyString; /// ditto
+    
+    @safe pure nothrow:
+    static void isCompatibleMatrixImpl(int r, int c)(Matrix!(mt, r, c) m) {
+    }
+    
+    template isCompatibleMatrix(T) {
+        enum isCompatibleMatrix = is(typeof(isCompatibleMatrixImpl(T.init)));
+    }
+    
+    static void isCompatibleVectorImpl(int d)(Vector!(mt, d) vec) {
+    }
+    
+    template isCompatibleVector(T) {
+        enum isCompatibleVector = is(typeof(isCompatibleVectorImpl(T.init)));
+    }
+    
+    private void construct(int i, T, Tail...)(T head, Tail tail) {
+        static if(i >= rows*cols) {
+            static assert(false, "Too many arguments passed to constructor");
+        } else static if(is(T : mt)) {
+            matrix[i / cols][i % cols] = head;
+            construct!(i + 1)(tail);
+        } else static if(is(T == Vector!(mt, cols))) {
+            static if(i % cols == 0) {
+                matrix[i / cols] = head.vector;
+                construct!(i + T.dimension)(tail);
+            } else {
+                static assert(false, "Can't convert Vector into the matrix. Maybe it doesn't align to the columns correctly or dimension doesn't fit");
+            }
+        } else static if(isDynamicArray!T) {
+            foreach(j; 0..cols*rows)
+                matrix[j / cols][j % cols] = head[j];
+        } else {
+            static assert(false, "Matrix constructor argument must be of type " ~ mt.stringof ~ " or Vector, not " ~ T.stringof);
+        }
+    }
+    
+    private void construct(int i)() { // terminate
+        static assert(i == rows*cols, "Not enough arguments passed to constructor");
+    }
+    
+    /// Constructs the matrix:
+    /// If a single value is passed, the matrix will be cleared with this value (each column in each row will contain this value).
+    /// If a matrix with more rows and columns is passed, the matrix will be the upper left nxm matrix.
+    /// If a matrix with less rows and columns is passed, the passed matrix will be stored in the upper left of an identity matrix.
+    /// It's also allowed to pass vectors and scalars at a time, but the vectors dimension must match the number of columns and align correctly.
+    /// Examples:
+    /// ---
+    /// mat2 m2 = mat2(0.0f); // mat2 m2 = mat2(0.0f, 0.0f, 0.0f, 0.0f);
+    /// mat3 m3 = mat3(m2); // mat3 m3 = mat3(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f);
+    /// mat3 m3_2 = mat3(vec3(1.0f, 2.0f, 3.0f), 4.0f, 5.0f, 6.0f, vec3(7.0f, 8.0f, 9.0f));
+    /// mat4 m4 = mat4.identity; // just an identity matrix
+    /// mat3 m3_3 = mat3(m4); // mat3 m3_3 = mat3.identity
+    /// ---
+    this(Args...)(Args args) {
+        construct!(0)(args);
+    }
+    
+    /// ditto
+    this(T)(T mat) if(is_matrix!T && (T.cols >= cols) && (T.rows >= rows)) {
+        foreach(r; TupleRange!(0, rows)) {
+            foreach(c; TupleRange!(0, cols)) {
+                matrix[r][c] = mat.matrix[r][c];
+            }
+        }
+    }
+    
+    /// ditto
+    this(T)(T mat) if(is_matrix!T && (T.cols < cols) && (T.rows < rows)) {
+        make_identity();
+        
+        foreach(r; TupleRange!(0, T.rows)) {
+            foreach(c; TupleRange!(0, T.cols)) {
+                matrix[r][c] = mat.matrix[r][c];
+            }
+        }
+    }
+    
+    /// ditto
+    this()(mt value) {
+        clear(value);
+    }
+    
+    /// Returns true if all values are not nan and finite, otherwise false.
+    @property bool isFinite() const {
+        static if(isIntegral!type) {
+            return true;
+        }
+        else {
+            foreach(row; matrix) {
+                foreach(col; row) {
+                    if(isNaN(col) || isInfinity(col)) {
+                        return false;
+                    }
+                }
+            }
+            return true;
+        }
+        
+    }
+    deprecated("Use isFinite instead of ok") alias ok = isFinite;
+    
+    /// Sets all values of the matrix to value (each column in each row will contain this value).
+    void clear(mt value) {
+        foreach(r; TupleRange!(0, rows)) {
+            foreach(c; TupleRange!(0, cols)) {
+                matrix[r][c] = value;
+            }
+        }
+    }
+    
+    unittest {
+        mat2 m2 = mat2(1.0f, 1.0f, vec2(2.0f, 2.0f));
+        assert(m2.matrix == [[1.0f, 1.0f], [2.0f, 2.0f]]);
+        m2.clear(3.0f);
+        assert(m2.matrix == [[3.0f, 3.0f], [3.0f, 3.0f]]);
+        assert(m2.isFinite);
+        m2.clear(float.nan);
+        assert(!m2.isFinite);
+        m2.clear(float.infinity);
+        assert(!m2.isFinite);
+        m2.clear(0.0f);
+        assert(m2.isFinite);
+        
+        mat3 m3 = mat3(1.0f);
+        assert(m3.matrix == [[1.0f, 1.0f, 1.0f],
+                [1.0f, 1.0f, 1.0f],
+                [1.0f, 1.0f, 1.0f]]);
+        
+        mat4 m4 = mat4(vec4(1.0f, 1.0f, 1.0f, 1.0f),
+            2.0f, 2.0f, 2.0f, 2.0f,
+            3.0f, 3.0f, 3.0f, 3.0f,
+            vec4(4.0f, 4.0f, 4.0f, 4.0f));
+        assert(m4.matrix == [[1.0f, 1.0f, 1.0f, 1.0f],
+                [2.0f, 2.0f, 2.0f, 2.0f],
+                [3.0f, 3.0f, 3.0f, 3.0f],
+                [4.0f, 4.0f, 4.0f, 4.0f]]);
+        assert(mat3(m4).matrix == [[1.0f, 1.0f, 1.0f],
+                [2.0f, 2.0f, 2.0f],
+                [3.0f, 3.0f, 3.0f]]);
+        assert(mat2(mat3(m4)).matrix == [[1.0f, 1.0f], [2.0f, 2.0f]]);
+        assert(mat2(m4).matrix == mat2(mat3(m4)).matrix);
+        assert(mat4(mat3(m4)).matrix == [[1.0f, 1.0f, 1.0f, 0.0f],
+                [2.0f, 2.0f, 2.0f, 0.0f],
+                [3.0f, 3.0f, 3.0f, 0.0f],
+                [0.0f, 0.0f, 0.0f, 1.0f]]);
+        
+        Matrix!(float, 2, 3) mt1 = Matrix!(float, 2, 3)(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f);
+        Matrix!(float, 3, 2) mt2 = Matrix!(float, 3, 2)(6.0f, -1.0f, 3.0f, 2.0f, 0.0f, -3.0f);
+        
+        assert(mt1.matrix == [[1.0f, 2.0f, 3.0f], [4.0f, 5.0f, 6.0f]]);
+        assert(mt2.matrix == [[6.0f, -1.0f], [3.0f, 2.0f], [0.0f, -3.0f]]);
+        
+        static assert(!__traits(compiles, mat2(1, 2, 1)));
+        static assert(!__traits(compiles, mat3(1, 2, 3, 1, 2, 3, 1, 2)));
+        static assert(!__traits(compiles, mat4(1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3)));
+        
+        auto m5 = mat2([0.0f,1,2,3]);
+        assert(m5.matrix == [[0.0f, 1.0f], [2.0f, 3.0f]]);
+        
+        auto m6 = Matrix!(int, 2, 3)([0,1,2,3,4,5]);
+        assert(m6.matrix == [[0, 1, 2], [3, 4, 5]]);
+    }
+    
+    static if(rows == cols) {
+        /// Makes the current matrix an identity matrix.
+        void make_identity() {
+            clear(0);
+            foreach(r; TupleRange!(0, rows)) {
+                matrix[r][r] = 1;
+            }
+        }
+        
+        /// Returns a identity matrix.
+        static @property Matrix identity() {
+            Matrix ret;
+            ret.clear(0);
+            
+            foreach(r; TupleRange!(0, rows)) {
+                ret.matrix[r][r] = 1;
+            }
+            
+            return ret;
+        }
+        
+        /// Transposes the current matrix;
+        void transpose() {
+            matrix = transposed().matrix;
+        }
+        
+        unittest {
+            mat2 m2 = mat2(1.0f);
+            m2.transpose();
+            assert(m2.matrix == mat2(1.0f).matrix);
+            m2.make_identity();
+            assert(m2.matrix == [[1.0f, 0.0f],
+                    [0.0f, 1.0f]]);
+            m2.transpose();
+            assert(m2.matrix == [[1.0f, 0.0f],
+                    [0.0f, 1.0f]]);
+            assert(m2.matrix == m2.identity.matrix);
+            
+            mat3 m3 = mat3(1.1f, 1.2f, 1.3f,
+                2.1f, 2.2f, 2.3f,
+                3.1f, 3.2f, 3.3f);
+            m3.transpose();
+            assert(m3.matrix == [[1.1f, 2.1f, 3.1f],
+                    [1.2f, 2.2f, 3.2f],
+                    [1.3f, 2.3f, 3.3f]]);
+            
+            mat4 m4 = mat4(2.0f);
+            m4.transpose();
+            assert(m4.matrix == mat4(2.0f).matrix);
+            m4.make_identity();
+            assert(m4.matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
+                    [0.0f, 1.0f, 0.0f, 0.0f],
+                    [0.0f, 0.0f, 1.0f, 0.0f],
+                    [0.0f, 0.0f, 0.0f, 1.0f]]);
+            assert(m4.matrix == m4.identity.matrix);
+        }
+        
+    }
+    
+    /// Returns a transposed copy of the matrix.
+    @property Matrix!(mt, cols, rows) transposed() const {
+        typeof(return) ret;
+        
+        foreach(r; TupleRange!(0, rows)) {
+            foreach(c; TupleRange!(0, cols)) {
+                ret.matrix[c][r] = matrix[r][c];
+            }
+        }
+        
+        return ret;
+    }
+    
+    // transposed already tested in last unittest
+    
+    
+    static if((rows == 2) && (cols == 2)) {
+        @property mt det() const {
+            return (matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]);
+        }
+        
+        private Matrix invert(ref Matrix mat) const {
+            static if(isFloatingPoint!mt && rmul) {
+                mt d = 1 / det;
+                
+                mat.matrix = [[matrix[1][1]*d, -matrix[0][1]*d],
+                    [-matrix[1][0]*d, matrix[0][0]*d]];
+            } else {
+                mt d = det;
+                
+                mat.matrix = [[matrix[1][1]/d, -matrix[0][1]/d],
+                    [-matrix[1][0]/d, matrix[0][0]/d]];
+            }
+            
+            return mat;
+        }
+        
+        static Matrix scaling(mt x, mt y) {
+            Matrix ret = Matrix.identity;
+            
+            ret.matrix[0][0] = x;
+            ret.matrix[1][1] = y;
+            
+            return ret;
+        }
+        
+        Matrix scale(mt x, mt y) {
+            this = Matrix.scaling(x, y) * this;
+            return this;
+        }
+        
+        unittest {
+            assert(mat2.scaling(3, 3).matrix == mat2.identity.scale(3, 3).matrix);
+            assert(mat2.scaling(3, 3).matrix == [[3.0f, 0.0f], [0.0f, 3.0f]]);
+        }
+        
+    } else static if((rows == 3) && (cols == 3)) {
+        @property mt det() const {
+            return (matrix[0][0] * matrix[1][1] * matrix[2][2]
+                + matrix[0][1] * matrix[1][2] * matrix[2][0]
+                + matrix[0][2] * matrix[1][0] * matrix[2][1]
+                - matrix[0][2] * matrix[1][1] * matrix[2][0]
+                - matrix[0][1] * matrix[1][0] * matrix[2][2]
+                - matrix[0][0] * matrix[1][2] * matrix[2][1]);
+        }
+        
+        private Matrix invert(ref Matrix mat) const {
+            static if(isFloatingPoint!mt && rmul) {
+                mt d = 1 / det;
+                enum op = "*";
+            } else {
+                mt d = det;
+                enum op = "/";
+            }
+            
+            mixin(`
+            mat.matrix = [[(matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])`~op~`d,
+                           (matrix[0][2] * matrix[2][1] - matrix[0][1] * matrix[2][2])`~op~`d,
+                           (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])`~op~`d],
+                          [(matrix[1][2] * matrix[2][0] - matrix[1][0] * matrix[2][2])`~op~`d,
+                           (matrix[0][0] * matrix[2][2] - matrix[0][2] * matrix[2][0])`~op~`d,
+                           (matrix[0][2] * matrix[1][0] - matrix[0][0] * matrix[1][2])`~op~`d],
+                          [(matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0])`~op~`d,
+                           (matrix[0][1] * matrix[2][0] - matrix[0][0] * matrix[2][1])`~op~`d,
+                           (matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0])`~op~`d]];
+            `);
+            
+            return mat;
+        }
+    } else static if((rows == 4) && (cols == 4)) {
+        /// Returns the determinant of the current matrix (2x2, 3x3 and 4x4 matrices).
+        @property mt det() const {
+            return (matrix[0][3] * matrix[1][2] * matrix[2][1] * matrix[3][0] - matrix[0][2] * matrix[1][3] * matrix[2][1] * matrix[3][0]
+                - matrix[0][3] * matrix[1][1] * matrix[2][2] * matrix[3][0] + matrix[0][1] * matrix[1][3] * matrix[2][2] * matrix[3][0]
+                + matrix[0][2] * matrix[1][1] * matrix[2][3] * matrix[3][0] - matrix[0][1] * matrix[1][2] * matrix[2][3] * matrix[3][0]
+                - matrix[0][3] * matrix[1][2] * matrix[2][0] * matrix[3][1] + matrix[0][2] * matrix[1][3] * matrix[2][0] * matrix[3][1]
+                + matrix[0][3] * matrix[1][0] * matrix[2][2] * matrix[3][1] - matrix[0][0] * matrix[1][3] * matrix[2][2] * matrix[3][1]
+                - matrix[0][2] * matrix[1][0] * matrix[2][3] * matrix[3][1] + matrix[0][0] * matrix[1][2] * matrix[2][3] * matrix[3][1]
+                + matrix[0][3] * matrix[1][1] * matrix[2][0] * matrix[3][2] - matrix[0][1] * matrix[1][3] * matrix[2][0] * matrix[3][2]
+                - matrix[0][3] * matrix[1][0] * matrix[2][1] * matrix[3][2] + matrix[0][0] * matrix[1][3] * matrix[2][1] * matrix[3][2]
+                + matrix[0][1] * matrix[1][0] * matrix[2][3] * matrix[3][2] - matrix[0][0] * matrix[1][1] * matrix[2][3] * matrix[3][2]
+                - matrix[0][2] * matrix[1][1] * matrix[2][0] * matrix[3][3] + matrix[0][1] * matrix[1][2] * matrix[2][0] * matrix[3][3]
+                + matrix[0][2] * matrix[1][0] * matrix[2][1] * matrix[3][3] - matrix[0][0] * matrix[1][2] * matrix[2][1] * matrix[3][3]
+                - matrix[0][1] * matrix[1][0] * matrix[2][2] * matrix[3][3] + matrix[0][0] * matrix[1][1] * matrix[2][2] * matrix[3][3]);
+        }
+        
+        private Matrix invert(ref Matrix mat) const {
+            static if(isFloatingPoint!mt && rmul) {
+                mt d = 1 / det;
+                enum op = "*";
+            } else {
+                mt d = det;
+                enum op = "/";
+            }
+            
+            mixin(`
+            mat.matrix = [[(matrix[1][1] * matrix[2][2] * matrix[3][3] + matrix[1][2] * matrix[2][3] * matrix[3][1] + matrix[1][3] * matrix[2][1] * matrix[3][2]
+                          - matrix[1][1] * matrix[2][3] * matrix[3][2] - matrix[1][2] * matrix[2][1] * matrix[3][3] - matrix[1][3] * matrix[2][2] * matrix[3][1])`~op~`d,
+                           (matrix[0][1] * matrix[2][3] * matrix[3][2] + matrix[0][2] * matrix[2][1] * matrix[3][3] + matrix[0][3] * matrix[2][2] * matrix[3][1]
+                          - matrix[0][1] * matrix[2][2] * matrix[3][3] - matrix[0][2] * matrix[2][3] * matrix[3][1] - matrix[0][3] * matrix[2][1] * matrix[3][2])`~op~`d,
+                           (matrix[0][1] * matrix[1][2] * matrix[3][3] + matrix[0][2] * matrix[1][3] * matrix[3][1] + matrix[0][3] * matrix[1][1] * matrix[3][2]
+                          - matrix[0][1] * matrix[1][3] * matrix[3][2] - matrix[0][2] * matrix[1][1] * matrix[3][3] - matrix[0][3] * matrix[1][2] * matrix[3][1])`~op~`d,
+                           (matrix[0][1] * matrix[1][3] * matrix[2][2] + matrix[0][2] * matrix[1][1] * matrix[2][3] + matrix[0][3] * matrix[1][2] * matrix[2][1]
+                          - matrix[0][1] * matrix[1][2] * matrix[2][3] - matrix[0][2] * matrix[1][3] * matrix[2][1] - matrix[0][3] * matrix[1][1] * matrix[2][2])`~op~`d],
+                          [(matrix[1][0] * matrix[2][3] * matrix[3][2] + matrix[1][2] * matrix[2][0] * matrix[3][3] + matrix[1][3] * matrix[2][2] * matrix[3][0]
+                          - matrix[1][0] * matrix[2][2] * matrix[3][3] - matrix[1][2] * matrix[2][3] * matrix[3][0] - matrix[1][3] * matrix[2][0] * matrix[3][2])`~op~`d,
+                           (matrix[0][0] * matrix[2][2] * matrix[3][3] + matrix[0][2] * matrix[2][3] * matrix[3][0] + matrix[0][3] * matrix[2][0] * matrix[3][2]
+                          - matrix[0][0] * matrix[2][3] * matrix[3][2] - matrix[0][2] * matrix[2][0] * matrix[3][3] - matrix[0][3] * matrix[2][2] * matrix[3][0])`~op~`d,
+                           (matrix[0][0] * matrix[1][3] * matrix[3][2] + matrix[0][2] * matrix[1][0] * matrix[3][3] + matrix[0][3] * matrix[1][2] * matrix[3][0]
+                          - matrix[0][0] * matrix[1][2] * matrix[3][3] - matrix[0][2] * matrix[1][3] * matrix[3][0] - matrix[0][3] * matrix[1][0] * matrix[3][2])`~op~`d,
+                           (matrix[0][0] * matrix[1][2] * matrix[2][3] + matrix[0][2] * matrix[1][3] * matrix[2][0] + matrix[0][3] * matrix[1][0] * matrix[2][2]
+                          - matrix[0][0] * matrix[1][3] * matrix[2][2] - matrix[0][2] * matrix[1][0] * matrix[2][3] - matrix[0][3] * matrix[1][2] * matrix[2][0])`~op~`d],
+                          [(matrix[1][0] * matrix[2][1] * matrix[3][3] + matrix[1][1] * matrix[2][3] * matrix[3][0] + matrix[1][3] * matrix[2][0] * matrix[3][1]
+                          - matrix[1][0] * matrix[2][3] * matrix[3][1] - matrix[1][1] * matrix[2][0] * matrix[3][3] - matrix[1][3] * matrix[2][1] * matrix[3][0])`~op~`d,
+                           (matrix[0][0] * matrix[2][3] * matrix[3][1] + matrix[0][1] * matrix[2][0] * matrix[3][3] + matrix[0][3] * matrix[2][1] * matrix[3][0]
+                          - matrix[0][0] * matrix[2][1] * matrix[3][3] - matrix[0][1] * matrix[2][3] * matrix[3][0] - matrix[0][3] * matrix[2][0] * matrix[3][1])`~op~`d,
+                           (matrix[0][0] * matrix[1][1] * matrix[3][3] + matrix[0][1] * matrix[1][3] * matrix[3][0] + matrix[0][3] * matrix[1][0] * matrix[3][1]
+                          - matrix[0][0] * matrix[1][3] * matrix[3][1] - matrix[0][1] * matrix[1][0] * matrix[3][3] - matrix[0][3] * matrix[1][1] * matrix[3][0])`~op~`d,
+                           (matrix[0][0] * matrix[1][3] * matrix[2][1] + matrix[0][1] * matrix[1][0] * matrix[2][3] + matrix[0][3] * matrix[1][1] * matrix[2][0]
+                          - matrix[0][0] * matrix[1][1] * matrix[2][3] - matrix[0][1] * matrix[1][3] * matrix[2][0] - matrix[0][3] * matrix[1][0] * matrix[2][1])`~op~`d],
+                          [(matrix[1][0] * matrix[2][2] * matrix[3][1] + matrix[1][1] * matrix[2][0] * matrix[3][2] + matrix[1][2] * matrix[2][1] * matrix[3][0]
+                          - matrix[1][0] * matrix[2][1] * matrix[3][2] - matrix[1][1] * matrix[2][2] * matrix[3][0] - matrix[1][2] * matrix[2][0] * matrix[3][1])`~op~`d,
+                           (matrix[0][0] * matrix[2][1] * matrix[3][2] + matrix[0][1] * matrix[2][2] * matrix[3][0] + matrix[0][2] * matrix[2][0] * matrix[3][1]
+                          - matrix[0][0] * matrix[2][2] * matrix[3][1] - matrix[0][1] * matrix[2][0] * matrix[3][2] - matrix[0][2] * matrix[2][1] * matrix[3][0])`~op~`d,
+                           (matrix[0][0] * matrix[1][2] * matrix[3][1] + matrix[0][1] * matrix[1][0] * matrix[3][2] + matrix[0][2] * matrix[1][1] * matrix[3][0]
+                          - matrix[0][0] * matrix[1][1] * matrix[3][2] - matrix[0][1] * matrix[1][2] * matrix[3][0] - matrix[0][2] * matrix[1][0] * matrix[3][1])`~op~`d,
+                           (matrix[0][0] * matrix[1][1] * matrix[2][2] + matrix[0][1] * matrix[1][2] * matrix[2][0] + matrix[0][2] * matrix[1][0] * matrix[2][1]
+                          - matrix[0][0] * matrix[1][2] * matrix[2][1] - matrix[0][1] * matrix[1][0] * matrix[2][2] - matrix[0][2] * matrix[1][1] * matrix[2][0])`~op~`d]];
+            `);
+            
+            return mat;
+        }
+        
+        // some static fun ...
+        // (1) glprogramming.com/red/appendixf.html - ortographic is broken!
+        // (2) http://fly.cc.fer.hr/~unreal/theredbook/appendixg.html
+        // (3) http://en.wikipedia.org/wiki/Orthographic_projection_(geometry)
+        
+        static if(isFloatingPoint!mt) {
+            static private mt[6] cperspective(mt width, mt height, mt fov, mt near, mt far)
+            in { assert(height != 0); }
+            body {
+                mt aspect = width/height;
+                mt top = near * tan(fov*(PI/360.0));
+                mt bottom = -top;
+                mt right = top * aspect;
+                mt left = -right;
+                
+                return [left, right, bottom, top, near, far];
+            }
+            
+            /// Returns a perspective matrix (4x4 and floating-point matrices only).
+            static Matrix perspective(mt width, mt height, mt fov, mt near, mt far) {
+                mt[6] cdata = cperspective(width, height, fov, near, far);
+                return perspective(cdata[0], cdata[1], cdata[2], cdata[3], cdata[4], cdata[5]);
+            }
+            
+            /// ditto
+            static Matrix perspective(mt left, mt right, mt bottom, mt top, mt near, mt far)
+            in {
+                assert(right-left != 0);
+                assert(top-bottom != 0);
+                assert(far-near != 0);
+            }
+            body {
+                Matrix ret;
+                ret.clear(0);
+                
+                ret.matrix[0][0] = (2*near)/(right-left);
+                ret.matrix[0][2] = (right+left)/(right-left);
+                ret.matrix[1][1] = (2*near)/(top-bottom);
+                ret.matrix[1][2] = (top+bottom)/(top-bottom);
+                ret.matrix[2][2] = -(far+near)/(far-near);
+                ret.matrix[2][3] = -(2*far*near)/(far-near);
+                ret.matrix[3][2] = -1;
+                
+                return ret;
+            }
+            
+            /// Returns an inverse perspective matrix (4x4 and floating-point matrices only).
+            static Matrix perspective_inverse(mt width, mt height, mt fov, mt near, mt far) {
+                mt[6] cdata = cperspective(width, height, fov, near, far);
+                return perspective_inverse(cdata[0], cdata[1], cdata[2], cdata[3], cdata[4], cdata[5]);
+            }
+            
+            /// ditto
+            static Matrix perspective_inverse(mt left, mt right, mt bottom, mt top, mt near, mt far)
+            in {
+                assert(near != 0);
+                assert(far != 0);
+            }
+            body {
+                Matrix ret;
+                ret.clear(0);
+                
+                ret.matrix[0][0] = (right-left)/(2*near);
+                ret.matrix[0][3] = (right+left)/(2*near);
+                ret.matrix[1][1] = (top-bottom)/(2*near);
+                ret.matrix[1][3] = (top+bottom)/(2*near);
+                ret.matrix[2][3] = -1;
+                ret.matrix[3][2] = -(far-near)/(2*far*near);
+                ret.matrix[3][3] = (far+near)/(2*far*near);
+                
+                return ret;
+            }
+            
+            // (2) and (3) say this one is correct
+            /// Returns an orthographic matrix (4x4 and floating-point matrices only).
+            static Matrix orthographic(mt left, mt right, mt bottom, mt top, mt near, mt far)
+            in {
+                assert(right-left != 0);
+                assert(top-bottom != 0);
+                assert(far-near != 0);
+            }
+            body {
+                Matrix ret;
+                ret.clear(0);
+                
+                ret.matrix[0][0] = 2/(right-left);
+                ret.matrix[0][3] = -(right+left)/(right-left);
+                ret.matrix[1][1] = 2/(top-bottom);
+                ret.matrix[1][3] = -(top+bottom)/(top-bottom);
+                ret.matrix[2][2] = -2/(far-near);
+                ret.matrix[2][3] = -(far+near)/(far-near);
+                ret.matrix[3][3] = 1;
+                
+                return ret;
+            }
+            
+            // (1) and (2) say this one is correct
+            /// Returns an inverse ortographic matrix (4x4 and floating-point matrices only).
+            static Matrix orthographic_inverse(mt left, mt right, mt bottom, mt top, mt near, mt far) {
+                Matrix ret;
+                ret.clear(0);
+                
+                ret.matrix[0][0] = (right-left)/2;
+                ret.matrix[0][3] = (right+left)/2;
+                ret.matrix[1][1] = (top-bottom)/2;
+                ret.matrix[1][3] = (top+bottom)/2;
+                ret.matrix[2][2] = (far-near)/-2;
+                ret.matrix[2][3] = (far+near)/2;
+                ret.matrix[3][3] = 1;
+                
+                return ret;
+            }
+            
+            /// Returns a look at matrix (4x4 and floating-point matrices only).
+            static Matrix look_at(Vector!(mt, 3) eye, Vector!(mt, 3) target, Vector!(mt, 3) up) {
+                alias Vector!(mt, 3) vec3mt;
+                vec3mt look_dir = (target - eye).normalized;
+                vec3mt up_dir = up.normalized;
+                
+                vec3mt right_dir = cross(look_dir, up_dir).normalized;
+                vec3mt perp_up_dir = cross(right_dir, look_dir);
+                
+                Matrix ret = Matrix.identity;
+                ret.matrix[0][0..3] = right_dir.vector[];
+                ret.matrix[1][0..3] = perp_up_dir.vector[];
+                ret.matrix[2][0..3] = (-look_dir).vector[];
+                
+                ret.matrix[0][3] = -dot(eye, right_dir);
+                ret.matrix[1][3] = -dot(eye, perp_up_dir);
+                ret.matrix[2][3] = dot(eye, look_dir);
+                
+                return ret;
+            }
+            
+            unittest {
+                mt[6] cp = cperspective(600f, 900f, 60f, 1f, 100f);
+                assert(cp[4] == 1.0f);
+                assert(cp[5] == 100.0f);
+                assert(cp[0] == -cp[1]);
+                assert((cp[0] < -0.38489f) && (cp[0] > -0.38491f));
+                assert(cp[2] == -cp[3]);
+                assert((cp[2] < -0.577349f) && (cp[2] > -0.577351f));
+                
+                assert(mat4.perspective(600f, 900f, 60.0, 1.0, 100.0) == mat4.perspective(cp[0], cp[1], cp[2], cp[3], cp[4], cp[5]));
+                float[4][4] m4p = mat4.perspective(600f, 900f, 60.0, 1.0, 100.0).matrix;
+                assert((m4p[0][0] < 2.598077f) && (m4p[0][0] > 2.598075f));
+                assert(m4p[0][2] == 0.0f);
+                assert((m4p[1][1] < 1.732052) && (m4p[1][1] > 1.732050));
+                assert(m4p[1][2] == 0.0f);
+                assert((m4p[2][2] < -1.020201) && (m4p[2][2] > -1.020203));
+                assert((m4p[2][3] < -2.020201) && (m4p[2][3] > -2.020203));
+                assert((m4p[3][2] < -0.9f) && (m4p[3][2] > -1.1f));
+                
+                float[4][4] m4pi = mat4.perspective_inverse(600f, 900f, 60.0, 1.0, 100.0).matrix;
+                assert((m4pi[0][0] < 0.384901) && (m4pi[0][0] > 0.384899));
+                assert(m4pi[0][3] == 0.0f);
+                assert((m4pi[1][1] < 0.577351) && (m4pi[1][1] > 0.577349));
+                assert(m4pi[1][3] == 0.0f);
+                assert(m4pi[2][3] == -1.0f);
+                assert((m4pi[3][2] < -0.494999) && (m4pi[3][2] > -0.495001));
+                assert((m4pi[3][3] < 0.505001) && (m4pi[3][3] > 0.504999));
+                
+                // maybe the next tests should be improved
+                float[4][4] m4o = mat4.orthographic(-1.0f, 1.0f, -1.0f, 1.0f, -1.0f, 1.0f).matrix;
+                assert(m4o == [[1.0f, 0.0f, 0.0f, 0.0f],
+                        [0.0f, 1.0f, 0.0f, 0.0f],
+                        [0.0f, 0.0f, -1.0f, 0.0f],
+                        [0.0f, 0.0f, 0.0f, 1.0f]]);
+                
+                float[4][4] m4oi = mat4.orthographic_inverse(-1.0f, 1.0f, -1.0f, 1.0f, -1.0f, 1.0f).matrix;
+                assert(m4oi == [[1.0f, 0.0f, 0.0f, 0.0f],
+                        [0.0f, 1.0f, 0.0f, 0.0f],
+                        [0.0f, 0.0f, -1.0f, 0.0f],
+                        [0.0f, 0.0f, 0.0f, 1.0f]]);
+                
+                //TODO: look_at tests
+            }
+        }
+    }
+    
+    static if((rows == cols) && (rows >= 3) && (rows <= 4)) {
+        /// Returns a translation matrix (3x3 and 4x4 matrices).
+        static Matrix translation(mt x, mt y, mt z) {
+            Matrix ret = Matrix.identity;
+            
+            ret.matrix[0][cols-1] = x;
+            ret.matrix[1][cols-1] = y;
+            ret.matrix[2][cols-1] = z;
+            
+            return ret;
+        }
+        
+        /// ditto
+        static Matrix translation(Vector!(mt, 3) v) {
+            Matrix ret = Matrix.identity;
+            
+            ret.matrix[0][cols-1] = v.x;
+            ret.matrix[1][cols-1] = v.y;
+            ret.matrix[2][cols-1] = v.z;
+            
+            return ret;
+        }
+        
+        /// Applys a translation on the current matrix and returns $(I this) (3x3 and 4x4 matrices).
+        Matrix translate(mt x, mt y, mt z) {
+            this = Matrix.translation(x, y, z) * this;
+            return this;
+        }
+        
+        /// ditto
+        Matrix translate(Vector!(mt, 3) v) {
+            this = Matrix.translation(v) * this;
+            return this;
+        }
+        
+        /// Returns a scaling matrix (3x3 and 4x4 matrices);
+        static Matrix scaling(mt x, mt y, mt z) {
+            Matrix ret = Matrix.identity;
+            
+            ret.matrix[0][0] = x;
+            ret.matrix[1][1] = y;
+            ret.matrix[2][2] = z;
+            
+            return ret;
+        }
+        
+        /// Applys a scale to the current matrix and returns $(I this) (3x3 and 4x4 matrices).
+        Matrix scale(mt x, mt y, mt z) {
+            this = Matrix.scaling(x, y, z) * this;
+            return this;
+        }
+        
+        unittest {
+            mat3 m3 = mat3.identity;
+            assert(m3.translate(1.0f, 2.0f, 3.0f).matrix == mat3.translation(1.0f, 2.0f, 3.0f).matrix);
+            assert(mat3.translation(1.0f, 2.0f, 3.0f).matrix == [[1.0f, 0.0f, 1.0f],
+                    [0.0f, 1.0f, 2.0f],
+                    [0.0f, 0.0f, 3.0f]]);
+            assert(mat3.identity.translate(0.0f, 1.0f, 2.0f).matrix == mat3.translation(0.0f, 1.0f, 2.0f).matrix);
+            
+            mat3 m31 = mat3.identity;
+            assert(m31.translate(vec3(1.0f, 2.0f, 3.0f)).matrix == mat3.translation(vec3(1.0f, 2.0f, 3.0f)).matrix);
+            assert(mat3.translation(vec3(1.0f, 2.0f, 3.0f)).matrix == [[1.0f, 0.0f, 1.0f],
+                    [0.0f, 1.0f, 2.0f],
+                    [0.0f, 0.0f, 3.0f]]);
+            assert(mat3.identity.translate(vec3(0.0f, 1.0f, 2.0f)).matrix == mat3.translation(vec3(0.0f, 1.0f, 2.0f)).matrix);
+            
+            assert(m3.scaling(0.0f, 1.0f, 2.0f).matrix == mat3.scaling(0.0f, 1.0f, 2.0f).matrix);
+            assert(mat3.scaling(0.0f, 1.0f, 2.0f).matrix == [[0.0f, 0.0f, 0.0f],
+                    [0.0f, 1.0f, 0.0f],
+                    [0.0f, 0.0f, 2.0f]]);
+            assert(mat3.identity.scale(0.0f, 1.0f, 2.0f).matrix == mat3.scaling(0.0f, 1.0f, 2.0f).matrix);
+            
+            // same tests for 4x4
+            
+            mat4 m4 = mat4(1.0f);
+            assert(m4.translation(1.0f, 2.0f, 3.0f).matrix == mat4.translation(1.0f, 2.0f, 3.0f).matrix);
+            assert(mat4.translation(1.0f, 2.0f, 3.0f).matrix == [[1.0f, 0.0f, 0.0f, 1.0f],
+                    [0.0f, 1.0f, 0.0f, 2.0f],
+                    [0.0f, 0.0f, 1.0f, 3.0f],
+                    [0.0f, 0.0f, 0.0f, 1.0f]]);
+            assert(mat4.identity.translate(0.0f, 1.0f, 2.0f).matrix == mat4.translation(0.0f, 1.0f, 2.0f).matrix);
+            
+            assert(m4.scaling(0.0f, 1.0f, 2.0f).matrix == mat4.scaling(0.0f, 1.0f, 2.0f).matrix);
+            assert(mat4.scaling(0.0f, 1.0f, 2.0f).matrix == [[0.0f, 0.0f, 0.0f, 0.0f],
+                    [0.0f, 1.0f, 0.0f, 0.0f],
+                    [0.0f, 0.0f, 2.0f, 0.0f],
+                    [0.0f, 0.0f, 0.0f, 1.0f]]);
+            assert(mat4.identity.scale(0.0f, 1.0f, 2.0f).matrix == mat4.scaling(0.0f, 1.0f, 2.0f).matrix);
+        }
+    }
+    
+    
+    static if((rows == cols) && (rows >= 3)) {
+        static if(isFloatingPoint!mt) {
+            /// Returns an identity matrix with an applied rotate_axis around an arbitrary axis (nxn matrices, n >= 3).
+            static Matrix rotation(real alpha, Vector!(mt, 3) axis) {
+                Matrix mult = Matrix.identity;
+                
+                if(axis.length != 1) {
+                    axis.normalize();
+                }
+                
+                real cosa = cos(alpha);
+                real sina = sin(alpha);
+                
+                Vector!(mt, 3) temp = (1 - cosa)*axis;
+                
+                mult.matrix[0][0] = to!mt(cosa + temp.x * axis.x);
+                mult.matrix[0][1] = to!mt(       temp.x * axis.y + sina * axis.z);
+                mult.matrix[0][2] = to!mt(       temp.x * axis.z - sina * axis.y);
+                mult.matrix[1][0] = to!mt(       temp.y * axis.x - sina * axis.z);
+                mult.matrix[1][1] = to!mt(cosa + temp.y * axis.y);
+                mult.matrix[1][2] = to!mt(       temp.y * axis.z + sina * axis.x);
+                mult.matrix[2][0] = to!mt(       temp.z * axis.x + sina * axis.y);
+                mult.matrix[2][1] = to!mt(       temp.z * axis.y - sina * axis.x);
+                mult.matrix[2][2] = to!mt(cosa + temp.z * axis.z);
+                
+                return mult;
+            }
+            
+            /// ditto
+            static Matrix rotation(real alpha, mt x, mt y, mt z) {
+                return Matrix.rotation(alpha, Vector!(mt, 3)(x, y, z));
+            }
+            
+            /// Returns an identity matrix with an applied rotation around the x-axis (nxn matrices, n >= 3).
+            static Matrix xrotation(real alpha) {
+                Matrix mult = Matrix.identity;
+                
+                mt cosamt = to!mt(cos(alpha));
+                mt sinamt = to!mt(sin(alpha));
+                
+                mult.matrix[1][1] = cosamt;
+                mult.matrix[1][2] = -sinamt;
+                mult.matrix[2][1] = sinamt;
+                mult.matrix[2][2] = cosamt;
+                
+                return mult;
+            }
+            
+            /// Returns an identity matrix with an applied rotation around the y-axis (nxn matrices, n >= 3).
+            static Matrix yrotation(real alpha) {
+                Matrix mult = Matrix.identity;
+                
+                mt cosamt = to!mt(cos(alpha));
+                mt sinamt = to!mt(sin(alpha));
+                
+                mult.matrix[0][0] = cosamt;
+                mult.matrix[0][2] = sinamt;
+                mult.matrix[2][0] = -sinamt;
+                mult.matrix[2][2] = cosamt;
+                
+                return mult;
+            }
+            
+            /// Returns an identity matrix with an applied rotation around the z-axis (nxn matrices, n >= 3).
+            static Matrix zrotation(real alpha) {
+                Matrix mult = Matrix.identity;
+                
+                mt cosamt = to!mt(cos(alpha));
+                mt sinamt = to!mt(sin(alpha));
+                
+                mult.matrix[0][0] = cosamt;
+                mult.matrix[0][1] = -sinamt;
+                mult.matrix[1][0] = sinamt;
+                mult.matrix[1][1] = cosamt;
+                
+                return mult;
+            }
+            
+            Matrix rotate(real alpha, Vector!(mt, 3) axis) {
+                this = rotation(alpha, axis) * this;
+                return this;
+            }
+            
+            /// Rotates the current matrix around the x-axis and returns $(I this) (nxn matrices, n >= 3).
+            Matrix rotatex(real alpha) {
+                this = xrotation(alpha) * this;
+                return this;
+            }
+            
+            /// Rotates the current matrix around the y-axis and returns $(I this) (nxn matrices, n >= 3).
+            Matrix rotatey(real alpha) {
+                this = yrotation(alpha) * this;
+                return this;
+            }
+            
+            /// Rotates the current matrix around the z-axis and returns $(I this) (nxn matrices, n >= 3).
+            Matrix rotatez(real alpha) {
+                this = zrotation(alpha) * this;
+                return this;
+            }
+            
+            unittest {
+                assert(mat4.xrotation(0).matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
+                        [0.0f, 1.0f, -0.0f, 0.0f],
+                        [0.0f, 0.0f, 1.0f, 0.0f],
+                        [0.0f, 0.0f, 0.0f, 1.0f]]);
+                assert(mat4.yrotation(0).matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
+                        [0.0f, 1.0f, 0.0f, 0.0f],
+                        [0.0f, 0.0f, 1.0f, 0.0f],
+                        [0.0f, 0.0f, 0.0f, 1.0f]]);
+                assert(mat4.zrotation(0).matrix == [[1.0f, -0.0f, 0.0f, 0.0f],
+                        [0.0f, 1.0f, 0.0f, 0.0f],
+                        [0.0f, 0.0f, 1.0f, 0.0f],
+                        [0.0f, 0.0f, 0.0f, 1.0f]]);
+                mat4 xro = mat4.identity;
+                xro.rotatex(0);
+                assert(mat4.xrotation(0).matrix == xro.matrix);
+                assert(xro.matrix == mat4.identity.rotatex(0).matrix);
+                assert(xro.matrix == mat4.rotation(0, vec3(1.0f, 0.0f, 0.0f)).matrix);
+                mat4 yro = mat4.identity;
+                yro.rotatey(0);
+                assert(mat4.yrotation(0).matrix == yro.matrix);
+                assert(yro.matrix == mat4.identity.rotatey(0).matrix);
+                assert(yro.matrix == mat4.rotation(0, vec3(0.0f, 1.0f, 0.0f)).matrix);
+                mat4 zro = mat4.identity;
+                xro.rotatez(0);
+                assert(mat4.zrotation(0).matrix == zro.matrix);
+                assert(zro.matrix == mat4.identity.rotatez(0).matrix);
+                assert(zro.matrix == mat4.rotation(0, vec3(0.0f, 0.0f, 1.0f)).matrix);
+            }
+        } // isFloatingPoint
+        
+        
+        /// Sets the translation of the matrix (nxn matrices, n >= 3).
+        void set_translation(mt[] values...) // intended to be a property
+        in { assert(values.length >= (rows-1)); }
+        body {
+            foreach(r; TupleRange!(0, rows-1)) {
+                matrix[r][rows-1] = values[r];
+            }
+        }
+        
+        /// Copyies the translation from mat to the current matrix (nxn matrices, n >= 3).
+        void set_translation(Matrix mat) {
+            foreach(r; TupleRange!(0, rows-1)) {
+                matrix[r][rows-1] = mat.matrix[r][rows-1];
+            }
+        }
+        
+        /// Returns an identity matrix with the current translation applied (nxn matrices, n >= 3)..
+        Matrix get_translation() {
+            Matrix ret = Matrix.identity;
+            
+            foreach(r; TupleRange!(0, rows-1)) {
+                ret.matrix[r][rows-1] = matrix[r][rows-1];
+            }
+            
+            return ret;
+        }
+        
+        unittest {
+            mat3 m3 = mat3(0.0f, 1.0f, 2.0f,
+                3.0f, 4.0f, 5.0f,
+                6.0f, 7.0f, 1.0f);
+            assert(m3.get_translation().matrix == [[1.0f, 0.0f, 2.0f], [0.0f, 1.0f, 5.0f], [0.0f, 0.0f, 1.0f]]);
+            m3.set_translation(mat3.identity);
+            assert(mat3.identity.matrix == m3.get_translation().matrix);
+            m3.set_translation([2.0f, 5.0f]);
+            assert(m3.get_translation().matrix == [[1.0f, 0.0f, 2.0f], [0.0f, 1.0f, 5.0f], [0.0f, 0.0f, 1.0f]]);
+            assert(mat3.identity.matrix == mat3.identity.get_translation().matrix);
+            
+            mat4 m4 = mat4(0.0f, 1.0f, 2.0f, 3.0f,
+                4.0f, 5.0f, 6.0f, 7.0f,
+                8.0f, 9.0f, 10.0f, 11.0f,
+                12.0f, 13.0f, 14.0f, 1.0f);
+            assert(m4.get_translation().matrix == [[1.0f, 0.0f, 0.0f, 3.0f],
+                    [0.0f, 1.0f, 0.0f, 7.0f],
+                    [0.0f, 0.0f, 1.0f, 11.0f],
+                    [0.0f, 0.0f, 0.0f, 1.0f]]);
+            m4.set_translation(mat4.identity);
+            assert(mat4.identity.matrix == m4.get_translation().matrix);
+            m4.set_translation([3.0f, 7.0f, 11.0f]);
+            assert(m4.get_translation().matrix == [[1.0f, 0.0f, 0.0f, 3.0f],
+                    [0.0f, 1.0f, 0.0f, 7.0f],
+                    [0.0f, 0.0f, 1.0f, 11.0f],
+                    [0.0f, 0.0f, 0.0f, 1.0f]]);
+            assert(mat4.identity.matrix == mat4.identity.get_translation().matrix);
+        }
+        
+        /// Sets the scale of the matrix (nxn matrices, n >= 3).
+        void set_scale(mt[] values...)
+        in { assert(values.length >= (rows-1)); }
+        body {
+            foreach(r; TupleRange!(0, rows-1)) {
+                matrix[r][r] = values[r];
+            }
+        }
+        
+        /// Copyies the scale from mat to the current matrix (nxn matrices, n >= 3).
+        void set_scale(Matrix mat) {
+            foreach(r; TupleRange!(0, rows-1)) {
+                matrix[r][r] = mat.matrix[r][r];
+            }
+        }
+        
+        /// Returns an identity matrix with the current scale applied (nxn matrices, n >= 3).
+        Matrix get_scale() {
+            Matrix ret = Matrix.identity;
+            
+            foreach(r; TupleRange!(0, rows-1)) {
+                ret.matrix[r][r] = matrix[r][r];
+            }
+            
+            return ret;
+        }
+        
+        unittest {
+            mat3 m3 = mat3(0.0f, 1.0f, 2.0f,
+                3.0f, 4.0f, 5.0f,
+                6.0f, 7.0f, 1.0f);
+            assert(m3.get_scale().matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 4.0f, 0.0f], [0.0f, 0.0f, 1.0f]]);
+            m3.set_scale(mat3.identity);
+            assert(mat3.identity.matrix == m3.get_scale().matrix);
+            m3.set_scale([0.0f, 4.0f]);
+            assert(m3.get_scale().matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 4.0f, 0.0f], [0.0f, 0.0f, 1.0f]]);
+            assert(mat3.identity.matrix == mat3.identity.get_scale().matrix);
+            
+            mat4 m4 = mat4(0.0f, 1.0f, 2.0f, 3.0f,
+                4.0f, 5.0f, 6.0f, 7.0f,
+                8.0f, 9.0f, 10.0f, 11.0f,
+                12.0f, 13.0f, 14.0f, 1.0f);
+            assert(m4.get_scale().matrix == [[0.0f, 0.0f, 0.0f, 0.0f],
+                    [0.0f, 5.0f, 0.0f, 0.0f],
+                    [0.0f, 0.0f, 10.0f, 0.0f],
+                    [0.0f, 0.0f, 0.0f, 1.0f]]);
+            m4.set_scale(mat4.identity);
+            assert(mat4.identity.matrix == m4.get_scale().matrix);
+            m4.set_scale([0.0f, 5.0f, 10.0f]);
+            assert(m4.get_scale().matrix == [[0.0f, 0.0f, 0.0f, 0.0f],
+                    [0.0f, 5.0f, 0.0f, 0.0f],
+                    [0.0f, 0.0f, 10.0f, 0.0f],
+                    [0.0f, 0.0f, 0.0f, 1.0f]]);
+            assert(mat4.identity.matrix == mat4.identity.get_scale().matrix);
+        }
+        
+        /// Copies rot into the upper left corner, the translation (nxn matrices, n >= 3).
+        void set_rotation(Matrix!(mt, 3, 3) rot) {
+            foreach(r; TupleRange!(0, 3)) {
+                foreach(c; TupleRange!(0, 3)) {
+                    matrix[r][c] = rot[r][c];
+                }
+            }
+        }
+        
+        /// Returns an identity matrix with the current rotation applied (nxn matrices, n >= 3).
+        Matrix!(mt, 3, 3) get_rotation() {
+            Matrix!(mt, 3, 3) ret = Matrix!(mt, 3, 3).identity;
+            
+            foreach(r; TupleRange!(0, 3)) {
+                foreach(c; TupleRange!(0, 3)) {
+                    ret.matrix[r][c] = matrix[r][c];
+                }
+            }
+            
+            return ret;
+        }
+        
+        unittest {
+            mat3 m3 = mat3(0.0f, 1.0f, 2.0f,
+                3.0f, 4.0f, 5.0f,
+                6.0f, 7.0f, 1.0f);
+            assert(m3.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [3.0f, 4.0f, 5.0f], [6.0f, 7.0f, 1.0f]]);
+            m3.set_rotation(mat3.identity);
+            assert(mat3.identity.matrix == m3.get_rotation().matrix);
+            m3.set_rotation(mat3(0.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 1.0f));
+            assert(m3.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [3.0f, 4.0f, 5.0f], [6.0f, 7.0f, 1.0f]]);
+            assert(mat3.identity.matrix == mat3.identity.get_rotation().matrix);
+            
+            mat4 m4 = mat4(0.0f, 1.0f, 2.0f, 3.0f,
+                4.0f, 5.0f, 6.0f, 7.0f,
+                8.0f, 9.0f, 10.0f, 11.0f,
+                12.0f, 13.0f, 14.0f, 1.0f);
+            assert(m4.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [4.0f, 5.0f, 6.0f], [8.0f, 9.0f, 10.0f]]);
+            m4.set_rotation(mat3.identity);
+            assert(mat3.identity.matrix == m4.get_rotation().matrix);
+            m4.set_rotation(mat3(0.0f, 1.0f, 2.0f, 4.0f, 5.0f, 6.0f, 8.0f, 9.0f, 10.0f));
+            assert(m4.get_rotation().matrix == [[0.0f, 1.0f, 2.0f], [4.0f, 5.0f, 6.0f], [8.0f, 9.0f, 10.0f]]);
+            assert(mat3.identity.matrix == mat4.identity.get_rotation().matrix);
+        }
+        
+    }
+    
+    static if((rows == cols) && (rows >= 2) && (rows <= 4)) {
+        /// Returns an inverted copy of the current matrix (nxn matrices, 2 >= n <= 4).
+        @property Matrix inverse() const {
+            Matrix mat;
+            invert(mat);
+            return mat;
+        }
+        
+        /// Inverts the current matrix (nxn matrices, 2 >= n <= 4).
+        void invert() {
+            // workaround Issue #11238
+            // uses a temporary instead of invert(this)
+            Matrix temp;
+            invert(temp);
+            this.matrix = temp.matrix;
+        }
+    }
+    
+    unittest {
+        mat2 m2 = mat2(1.0f, 2.0f, vec2(3.0f, 4.0f));
+        assert(m2.det == -2.0f);
+        assert(m2.inverse.matrix == [[-2.0f, 1.0f], [1.5f, -0.5f]]);
+        
+        mat3 m3 = mat3(1.0f, -2.0f, 3.0f,
+            7.0f, -1.0f, 0.0f,
+            3.0f, 2.0f, -4.0f);
+        assert(m3.det == -1.0f);
+        assert(m3.inverse.matrix == [[-4.0f, 2.0f, -3.0f],
+                [-28.0f, 13.0f, -21.0f],
+                [-17.0f, 8.0f, -13.0f]]);
+        
+        mat4 m4 = mat4(1.0f, 2.0f, 3.0f, 4.0f,
+            -2.0f, 1.0f, 5.0f, -2.0f,
+            2.0f, -1.0f, 7.0f, 1.0f,
+            3.0f, -3.0f, 2.0f, 0.0f);
+        assert(m4.det == -8.0f);
+        assert(m4.inverse.matrix == [[6.875f, 7.875f, -11.75f, 11.125f],
+                [6.625f, 7.625f, -11.25f, 10.375f],
+                [-0.375f, -0.375f, 0.75f, -0.625f],
+                [-4.5f, -5.5f, 8.0f, -7.5f]]);
+    }
+    
+    private void mms(mt inp, ref Matrix mat) const { // mat * scalar
+        for(int r = 0; r < rows; r++) {
+            for(int c = 0; c < cols; c++) {
+                mat.matrix[r][c] = matrix[r][c] * inp;
+            }
+        }
+    }
+    
+    private void masm(string op)(Matrix inp, ref Matrix mat) const { // mat + or - mat
+        foreach(r; TupleRange!(0, rows)) {
+            foreach(c; TupleRange!(0, cols)) {
+                mat.matrix[r][c] = mixin("inp.matrix[r][c]" ~ op ~ "matrix[r][c]");
+            }
+        }
+    }
+    
+    Matrix!(mt, rows, T.cols) opBinary(string op : "*", T)(T inp) const if(isCompatibleMatrix!T && (T.rows == cols)) {
+        Matrix!(mt, rows, T.cols) ret;
+        
+        foreach(r; TupleRange!(0, rows)) {
+            foreach(c; TupleRange!(0, T.cols)) {
+                ret.matrix[r][c] = 0;
+                
+                foreach(c2; TupleRange!(0, cols)) {
+                    ret.matrix[r][c] += matrix[r][c2] * inp.matrix[c2][c];
+                }
+            }
+        }
+        
+        return ret;
+    }
+    
+    Vector!(mt, rows) opBinary(string op : "*", T : Vector!(mt, cols))(T inp) const {
+        Vector!(mt, rows) ret;
+        ret.clear(0);
+        
+        foreach(c; TupleRange!(0, cols)) {
+            foreach(r; TupleRange!(0, rows)) {
+                ret.vector[r] += matrix[r][c] * inp.vector[c];
+            }
+        }
+        
+        return ret;
+    }
+
+    /+Vector!(mt, rows) opBinaryRight(string op : "*", T : Vector!(mt, cols))(T inp) const {
+        return opBinary!(op,T)(inp);
+    }+/
+    
+    Matrix opBinary(string op : "*")(mt inp) const {
+        Matrix ret;
+        mms(inp, ret);
+        return ret;
+    }
+    
+    Matrix opBinaryRight(string op : "*")(mt inp) const {
+        return this.opBinary!(op)(inp);
+    }
+    
+    Matrix opBinary(string op)(Matrix inp) const if((op == "+") || (op == "-")) {
+        Matrix ret;
+        masm!(op)(inp, ret);
+        return ret;
+    }
+    
+    void opOpAssign(string op : "*")(mt inp) {
+        mms(inp, this);
+    }
+    
+    void opOpAssign(string op)(Matrix inp) if((op == "+") || (op == "-")) {
+        masm!(op)(inp, this);
+    }
+    
+    void opOpAssign(string op)(Matrix inp) if(op == "*") {
+        this = this * inp;
+    }
+    
+    unittest {
+        mat2 m2 = mat2(1.0f, 2.0f, 3.0f, 4.0f);
+        vec2 v2 = vec2(2.0f, 2.0f);
+        assert((m2*2).matrix == [[2.0f, 4.0f], [6.0f, 8.0f]]);
+        assert((2*m2).matrix == (m2*2).matrix);
+        m2 *= 2;
+        assert(m2.matrix == [[2.0f, 4.0f], [6.0f, 8.0f]]);
+        assert((m2*v2).vector == [12.0f, 28.0f]);
+        assert((v2*m2).vector == [16.0f, 24.0f]);
+        assert((m2*m2).matrix == [[28.0f, 40.0f], [60.0f, 88.0f]]);
+        assert((m2-m2).matrix == [[0.0f, 0.0f], [0.0f, 0.0f]]);
+        assert((m2+m2).matrix == [[4.0f, 8.0f], [12.0f, 16.0f]]);
+        m2 += m2;
+        assert(m2.matrix == [[4.0f, 8.0f], [12.0f, 16.0f]]);
+        m2 -= m2;
+        assert(m2.matrix == [[0.0f, 0.0f], [0.0f, 0.0f]]);
+        
+        mat3 m3 = mat3(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f);
+        vec3 v3 = vec3(2.0f, 2.0f, 2.0f);
+        assert((m3*2).matrix == [[2.0f, 4.0f, 6.0f], [8.0f, 10.0f, 12.0f], [14.0f, 16.0f, 18.0f]]);
+        assert((2*m3).matrix == (m3*2).matrix);
+        m3 *= 2;
+        assert(m3.matrix == [[2.0f, 4.0f, 6.0f], [8.0f, 10.0f, 12.0f], [14.0f, 16.0f, 18.0f]]);
+        assert((m3*v3).vector == [24.0f, 60.0f, 96.0f]);
+        assert((v3*m3).vector == [48.0f, 60.0f, 72.0f]);
+        assert((m3*m3).matrix == [[120.0f, 144.0f, 168.0f], [264.0f, 324.0f, 384.0f], [408.0f, 504.0f, 600.0f]]);
+        assert((m3-m3).matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f]]);
+        assert((m3+m3).matrix == [[4.0f, 8.0f, 12.0f], [16.0f, 20.0f, 24.0f], [28.0f, 32.0f, 36.0f]]);
+        m3 += m3;
+        assert(m3.matrix == [[4.0f, 8.0f, 12.0f], [16.0f, 20.0f, 24.0f], [28.0f, 32.0f, 36.0f]]);
+        m3 -= m3;
+        assert(m3.matrix == [[0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f], [0.0f, 0.0f, 0.0f]]);
+        
+        // test opOpAssign for matrix multiplication
+        auto m4 = mat4.translation(0,1,2);
+        m4 *= mat4.translation(0,-1,2);
+        assert(m4 == mat4.translation(0,0,4));
+        
+        //TODO: tests for mat4, mat34
+    }
+    
+    // opEqual => "alias matrix this;"
+    
+    bool opCast(T : bool)() const {
+        return isFinite;
+    }
+    
+    unittest {
+        assert(mat2(1.0f, 2.0f, 1.0f, 1.0f) == mat2(1.0f, 2.0f, 1.0f, 1.0f));
+        assert(mat2(1.0f, 2.0f, 1.0f, 1.0f) != mat2(1.0f, 1.0f, 1.0f, 1.0f));
+        
+        assert(mat3(1.0f) == mat3(1.0f));
+        assert(mat3(1.0f) != mat3(2.0f));
+        
+        assert(mat4(1.0f) == mat4(1.0f));
+        assert(mat4(1.0f) != mat4(2.0f));
+        
+        assert(!(mat4(float.nan)));
+        if(mat4(1.0f)) { }
+        else { assert(false); }
+    }
+    
+}
+
+/// Pre-defined matrix types, the first number represents the number of rows
+/// and the second the number of columns, if there's just one it's a nxn matrix.
+/// All of these matrices are floating-point matrices.
+alias Matrix!(float, 2, 2) mat2;
+alias Matrix!(float, 3, 3) mat3;
+alias Matrix!(float, 3, 4) mat34;
+alias Matrix!(float, 4, 4) mat4;
+
+private unittest {
+    Matrix!(float,  1, 1) A = 1;
+    Matrix!(double, 1, 1) B = 1;
+    Matrix!(real,   1, 1) C = 1;
+    Matrix!(int,    1, 1) D = 1;
+    Matrix!(float,  5, 1) E = 1;
+    Matrix!(double, 5, 1) F = 1;
+    Matrix!(real,   5, 1) G = 1;
+    Matrix!(int,    5, 1) H = 1;
+    Matrix!(float,  1, 5) I = 1;
+    Matrix!(double, 1, 5) J = 1;
+    Matrix!(real,   1, 5) K = 1;
+    Matrix!(int,    1, 5) L = 1;
+}
diff --git a/gl3n/linalg/package.d b/gl3n/linalg/package.d
new file mode 100644
index 0000000..e5472b3
--- /dev/null
+++ b/gl3n/linalg/package.d
@@ -0,0 +1,23 @@
+/**
+gl3n.linalg.matrix
+
+Special thanks to:
+$(UL
+  $(LI Tomasz Stachowiak (h3r3tic): allowed me to use parts of $(LINK2 https://bitbucket.org/h3r3tic/boxen/src/default/src/xf/omg, omg).)
+  $(LI Jakob Øvrum (jA_cOp): improved the code a lot!)
+  $(LI Florian Boesch (___doc__): helps me to understand opengl/complex maths better, see: $(LINK http://codeflow.org/).)
+  $(LI #D on freenode: answered general questions about D.)
+)
+
+Authors: David Herberth, Stephan Dilly
+License: MIT
+
+Note: All methods marked with pure are weakly pure since, they all access an instance member.
+All static methods are strongly pure.
+*/
+ 
+module gl3n.linalg;
+
+public import gl3n.linalg.vec;
+public import gl3n.linalg.matrix;
+public import gl3n.linalg.quaternion;
\ No newline at end of file
diff --git a/gl3n/linalg/quaternion.d b/gl3n/linalg/quaternion.d
new file mode 100644
index 0000000..01198a4
--- /dev/null
+++ b/gl3n/linalg/quaternion.d
@@ -0,0 +1,645 @@
+/**
+gl3n.linalg.quaternion
+
+Special thanks to:
+$(UL
+  $(LI Tomasz Stachowiak (h3r3tic): allowed me to use parts of $(LINK2 https://bitbucket.org/h3r3tic/boxen/src/default/src/xf/omg, omg).)
+  $(LI Jakob Øvrum (jA_cOp): improved the code a lot!)
+  $(LI Florian Boesch (___doc__): helps me to understand opengl/complex maths better, see: $(LINK http://codeflow.org/).)
+  $(LI #D on freenode: answered general questions about D.)
+)
+
+Authors: David Herberth, Stephan Dilly
+License: MIT
+
+Note: All methods marked with pure are weakly pure since, they all access an instance member.
+All static methods are strongly pure.
+*/
+
+module gl3n.linalg.quaternion;
+
+import gl3n.linalg.vec;
+import gl3n.linalg.matrix;
+import gl3n.math : atan2, asin, almost_equal;
+
+/// Base template for all quaternion-types.
+/// Params:
+///  type = all values get stored as this type
+struct Quaternion(type) {
+    alias type qt; /// Holds the internal type of the quaternion.
+    
+    qt[4] quaternion; /// Holds the w, x, y and z coordinates.
+    
+    /// Returns a pointer to the quaternion in memory, it starts with the w coordinate.
+    @property auto value_ptr() const { return quaternion.ptr; }
+    
+    /// Returns the current vector formatted as string, useful for printing the quaternion.
+    @property string as_string() {
+        return format("%s", quaternion);
+    }
+    alias as_string toString;
+    
+    @safe pure nothrow:
+    private @property qt get_(char coord)() const {
+        return quaternion[coord_to_index!coord];
+    }
+    private @property void set_(char coord)(qt value) {
+        quaternion[coord_to_index!coord] = value;
+    }
+    
+    alias get_!'w' w; /// static properties to access the values.
+    alias set_!'w' w;
+    alias get_!'x' x; /// ditto
+    alias set_!'x' x;
+    alias get_!'y' y; /// ditto
+    alias set_!'y' y;
+    alias get_!'z' z; /// ditto
+    alias set_!'z' z;
+    
+    /// Constructs the quaternion.
+    /// Takes a 4-dimensional vector, where vector.x = the quaternions w coordinate,
+    /// or a w coordinate of type $(I qt) and a 3-dimensional vector representing the imaginary part,
+    /// or 4 values of type $(I qt).
+    this(qt w_, qt x_, qt y_, qt z_) {
+        w = w_;
+        x = x_;
+        y = y_;
+        z = z_;
+    }
+    
+    /// ditto
+    this(qt w_, Vector!(qt, 3) vec) {
+        w = w_;
+        quaternion[1..4] = vec.vector[];
+    }
+    
+    /// ditto
+    this(Vector!(qt, 4) vec) {
+        quaternion[] = vec.vector[];
+    }
+    
+    /// Returns true if all values are not nan and finite, otherwise false.
+    @property bool isFinite() const {
+        foreach(q; quaternion) {
+            if(isNaN(q) || isInfinity(q)) {
+                return false;
+            }
+        }
+        return true;
+    }
+    deprecated("Use isFinite instead of ok") alias ok = isFinite;
+    
+    unittest {
+        quat q1 = quat(0.0f, 0.0f, 0.0f, 1.0f);
+        assert(q1.quaternion == [0.0f, 0.0f, 0.0f, 1.0f]);
+        assert(q1.quaternion == quat(0.0f, 0.0f, 0.0f, 1.0f).quaternion);
+        assert(q1.quaternion == quat(0.0f, vec3(0.0f, 0.0f, 1.0f)).quaternion);
+        assert(q1.quaternion == quat(vec4(0.0f, 0.0f, 0.0f, 1.0f)).quaternion);
+        
+        assert(q1.isFinite);
+        q1.x = float.infinity;
+        assert(!q1.isFinite);
+        q1.x = float.nan;
+        assert(!q1.isFinite);
+        q1.x = 0.0f;
+        assert(q1.isFinite);
+    }
+    
+    template coord_to_index(char c) {
+        static if(c == 'w') {
+            enum coord_to_index = 0;
+        } else static if(c == 'x') {
+            enum coord_to_index = 1;
+        } else static if(c == 'y') {
+            enum coord_to_index = 2;
+        } else static if(c == 'z') {
+            enum coord_to_index = 3;
+        } else {
+            static assert(false, "accepted coordinates are x, y, z and w not " ~ c ~ ".");
+        }
+    }
+    
+    /// Returns the squared magnitude of the quaternion.
+    @property real magnitude_squared() const {
+        return to!real(w^^2 + x^^2 + y^^2 + z^^2);
+    }
+    
+    /// Returns the magnitude of the quaternion.
+    @property real magnitude() const {
+        return sqrt(magnitude_squared);
+    }
+    
+    /// Returns an identity quaternion (w=1, x=0, y=0, z=0).
+    static @property Quaternion identity() {
+        return Quaternion(1, 0, 0, 0);
+    }
+    
+    /// Makes the current quaternion an identity quaternion.
+    void make_identity() {
+        w = 1;
+        x = 0;
+        y = 0;
+        z = 0;
+    }
+    
+    /// Inverts the quaternion.
+    void invert() {
+        x = -x;
+        y = -y;
+        z = -z;
+    }
+    alias invert conjugate; /// ditto
+    
+    /// Returns an inverted copy of the current quaternion.
+    @property Quaternion inverse() const {
+        return Quaternion(w, -x, -y, -z);
+    }
+    alias inverse conjugated; /// ditto
+    
+    unittest {
+        quat q1 = quat(1.0f, 1.0f, 1.0f, 1.0f);
+        
+        assert(q1.magnitude == 2.0f);
+        assert(q1.magnitude_squared == 4.0f);
+        assert(q1.magnitude == quat(0.0f, 0.0f, 2.0f, 0.0f).magnitude);
+        
+        quat q2 = quat.identity;
+        assert(q2.quaternion == [1.0f, 0.0f, 0.0f, 0.0f]);
+        assert(q2.x == 0.0f);
+        assert(q2.y == 0.0f);
+        assert(q2.z == 0.0f);
+        assert(q2.w == 1.0f);
+        
+        assert(q1.inverse.quaternion == [1.0f, -1.0f, -1.0f, -1.0f]);
+        q1.invert();
+        assert(q1.quaternion == [1.0f, -1.0f, -1.0f, -1.0f]);
+        
+        q1.make_identity();
+        assert(q1.quaternion == q2.quaternion);
+        
+    }
+    
+    /// Creates a quaternion from a 3x3 matrix.
+    /// Params:
+    ///  matrix = 3x3 matrix (rotation)
+    /// Returns: A quaternion representing the rotation (3x3 matrix)
+    static Quaternion from_matrix(Matrix!(qt, 3, 3) matrix) {
+        Quaternion ret;
+        
+        auto mat = matrix.matrix;
+        qt trace = mat[0][0] + mat[1][1] + mat[2][2];
+        
+        if(trace > 0) {
+            real s = 0.5 / sqrt(trace + 1.0f);
+            
+            ret.w = to!qt(0.25 / s);
+            ret.x = to!qt((mat[2][1] - mat[1][2]) * s);
+            ret.y = to!qt((mat[0][2] - mat[2][0]) * s);
+            ret.z = to!qt((mat[1][0] - mat[0][1]) * s);
+        } else if((mat[0][0] > mat[1][1]) && (mat[0][0] > mat[2][2])) {
+            real s = 2.0 * sqrt(1.0 + mat[0][0] - mat[1][1] - mat[2][2]);
+            
+            ret.w = to!qt((mat[2][1] - mat[1][2]) / s);
+            ret.x = to!qt(0.25f * s);
+            ret.y = to!qt((mat[0][1] + mat[1][0]) / s);
+            ret.z = to!qt((mat[0][2] + mat[2][0]) / s);
+        } else if(mat[1][1] > mat[2][2]) {
+            real s = 2.0 * sqrt(1 + mat[1][1] - mat[0][0] - mat[2][2]);
+            
+            ret.w = to!qt((mat[0][2] - mat[2][0]) / s);
+            ret.x = to!qt((mat[0][1] + mat[1][0]) / s);
+            ret.y = to!qt(0.25f * s);
+            ret.z = to!qt((mat[1][2] + mat[2][1]) / s);
+        } else {
+            real s = 2.0 * sqrt(1 + mat[2][2] - mat[0][0] - mat[1][1]);
+            
+            ret.w = to!qt((mat[1][0] - mat[0][1]) / s);
+            ret.x = to!qt((mat[0][2] + mat[2][0]) / s);
+            ret.y = to!qt((mat[1][2] + mat[2][1]) / s);
+            ret.z = to!qt(0.25f * s);
+        }
+        
+        return ret;
+    }
+    
+    /// Returns the quaternion as matrix.
+    /// Params:
+    ///  rows = number of rows of the resulting matrix (min 3)
+    ///  cols = number of columns of the resulting matrix (min 3)
+    Matrix!(qt, rows, cols) to_matrix(int rows, int cols)() const if((rows >= 3) && (cols >= 3)) {
+        static if((rows == 3) && (cols == 3)) {
+            Matrix!(qt, rows, cols) ret;
+        } else {
+            Matrix!(qt, rows, cols) ret = Matrix!(qt, rows, cols).identity;
+        }
+        
+        qt xx = x^^2;
+        qt xy = x * y;
+        qt xz = x * z;
+        qt xw = x * w;
+        qt yy = y^^2;
+        qt yz = y * z;
+        qt yw = y * w;
+        qt zz = z^^2;
+        qt zw = z * w;
+        
+        ret.matrix[0][0] = 1 - 2 * (yy + zz);
+        ret.matrix[0][1] = 2 * (xy - zw);
+        ret.matrix[0][2] = 2 * (xz + yw);
+        
+        ret.matrix[1][0] = 2 * (xy + zw);
+        ret.matrix[1][1] = 1 - 2 * (xx + zz);
+        ret.matrix[1][2] = 2 * (yz - xw);
+        
+        ret.matrix[2][0] = 2 * (xz - yw);
+        ret.matrix[2][1] = 2 * (yz + xw);
+        ret.matrix[2][2] = 1 - 2 * (xx + yy);
+        
+        return ret;
+    }
+    
+    unittest {
+        quat q1 = quat(4.0f, 1.0f, 2.0f, 3.0f);
+        
+        assert(q1.to_matrix!(3, 3).matrix == [[-25.0f, -20.0f, 22.0f], [28.0f, -19.0f, 4.0f], [-10.0f, 20.0f, -9.0f]]);
+        assert(q1.to_matrix!(4, 4).matrix == [[-25.0f, -20.0f, 22.0f, 0.0f],
+                [28.0f, -19.0f, 4.0f, 0.0f],
+                [-10.0f, 20.0f, -9.0f, 0.0f],
+                [0.0f, 0.0f, 0.0f, 1.0f]]);
+        assert(quat.identity.to_matrix!(3, 3).matrix == Matrix!(qt, 3, 3).identity.matrix);
+        assert(q1.quaternion == quat.from_matrix(q1.to_matrix!(3, 3)).quaternion);
+        
+        assert(quat(1.0f, 0.0f, 0.0f, 0.0f).quaternion == quat.from_matrix(mat3.identity).quaternion);
+        
+        quat q2 = quat.from_matrix(mat3(1.0f, 3.0f, 2.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f));
+        assert(q2.x == 0.0f);
+        assert(almost_equal(q2.y, 0.7071067f));
+        assert(almost_equal(q2.z, -1.060660));
+        assert(almost_equal(q2.w, 0.7071067f));
+    }
+    
+    /// Normalizes the current quaternion.
+    void normalize() {
+        qt m = to!qt(magnitude);
+        
+        if(m != 0) {
+            w = w / m;
+            x = x / m;
+            y = y / m;
+            z = z / m;
+        }
+    }
+    
+    /// Returns a normalized copy of the current quaternion.
+    Quaternion normalized() const {
+        Quaternion ret;
+        qt m = to!qt(magnitude);
+        
+        if(m != 0) {
+            ret.w = w / m;
+            ret.x = x / m;
+            ret.y = y / m;
+            ret.z = z / m;
+        } else {
+            ret = Quaternion(w, x, y, z);
+        }
+        
+        return ret;
+    }
+    
+    unittest {
+        quat q1 = quat(1.0f, 2.0f, 3.0f, 4.0f);
+        quat q2 = quat(1.0f, 2.0f, 3.0f, 4.0f);
+        
+        q1.normalize();
+        assert(q1.quaternion == q2.normalized.quaternion);
+        //assert(q1.quaternion == q1.normalized.quaternion);
+        assert(almost_equal(q1.magnitude, 1.0));
+    }
+    
+    /// Returns the yaw.
+    @property real yaw() const {
+        return atan2(to!real(2.0*(w*z + x*y)), to!real(1.0 - 2.0*(y*y + z*z)));
+    }
+    
+    /// Returns the pitch.
+    @property real pitch() const {
+        return asin(to!real(2.0*(w*y - z*x)));
+    }
+    
+    /// Returns the roll.
+    @property real roll() const {
+        return atan2(to!real(2.0*(w*x + y*z)), to!real(1.0 - 2.0*(x*x + y*y)));
+    }
+    
+    unittest {
+        quat q1 = quat.identity;
+        assert(q1.pitch == 0.0f);
+        assert(q1.yaw == 0.0f);
+        assert(q1.roll == 0.0f);
+        
+        quat q2 = quat(1.0f, 1.0f, 1.0f, 1.0f);
+        assert(almost_equal(q2.yaw, q2.roll));
+        //assert(almost_equal(q2.yaw, 1.570796f));
+        assert(q2.pitch == 0.0f);
+        
+        quat q3 = quat(0.1f, 1.9f, 2.1f, 1.3f);
+        //assert(almost_equal(q3.yaw, 2.4382f));
+        assert(isNaN(q3.pitch));
+        //assert(almost_equal(q3.roll, 1.67719f));
+    }
+    
+    /// Returns a quaternion with applied rotation around the x-axis.
+    static Quaternion xrotation(real alpha) {
+        Quaternion ret;
+        
+        alpha /= 2;
+        ret.w = to!qt(cos(alpha));
+        ret.x = to!qt(sin(alpha));
+        ret.y = 0;
+        ret.z = 0;
+        
+        return ret;
+    }
+    
+    /// Returns a quaternion with applied rotation around the y-axis.
+    static Quaternion yrotation(real alpha) {
+        Quaternion ret;
+        
+        alpha /= 2;
+        ret.w = to!qt(cos(alpha));
+        ret.x = 0;
+        ret.y = to!qt(sin(alpha));
+        ret.z = 0;
+        
+        return ret;
+    }
+    
+    /// Returns a quaternion with applied rotation around the z-axis.
+    static Quaternion zrotation(real alpha) {
+        Quaternion ret;
+        
+        alpha /= 2;
+        ret.w = to!qt(cos(alpha));
+        ret.x = 0;
+        ret.y = 0;
+        ret.z = to!qt(sin(alpha));
+        
+        return ret;
+    }
+    
+    /// Returns a quaternion with applied rotation around an axis.
+    static Quaternion axis_rotation(real alpha, Vector!(qt, 3) axis) {
+        if(alpha == 0) {
+            return Quaternion.identity;
+        }
+        Quaternion ret;
+        
+        alpha /= 2;
+        qt sinaqt = to!qt(sin(alpha));
+        
+        ret.w = to!qt(cos(alpha));
+        ret.x = axis.x * sinaqt;
+        ret.y = axis.y * sinaqt;
+        ret.z = axis.z * sinaqt;
+        
+        return ret;
+    }
+    
+    /// Creates a quaternion from an euler rotation.
+    static Quaternion euler_rotation(real roll, real pitch, real yaw) {
+        Quaternion ret;
+        
+        auto cr = cos(roll / 2.0);
+        auto cp = cos(pitch / 2.0);
+        auto cy = cos(yaw / 2.0);
+        auto sr = sin(roll / 2.0);
+        auto sp = sin(pitch / 2.0);
+        auto sy = sin(yaw / 2.0);
+        
+        ret.w = cr * cp * cy + sr * sp * sy;
+        ret.x = sr * cp * cy - cr * sp * sy;
+        ret.y = cr * sp * cy + sr * cp * sy;
+        ret.z = cr * cp * sy - sr * sp * cy;
+        
+        return ret;
+    }
+    
+    unittest {
+        enum startPitch = 0.1;
+        enum startYaw = -0.2;
+        enum startRoll = 0.6;
+        
+        auto q = quat.euler_rotation(startRoll,startPitch,startYaw);
+        
+        assert(almost_equal(q.pitch,startPitch));
+        assert(almost_equal(q.yaw,startYaw));
+        assert(almost_equal(q.roll,startRoll));
+    }
+    
+    /// Rotates the current quaternion around the x-axis and returns $(I this).
+    Quaternion rotatex(real alpha) {
+        this = xrotation(alpha) * this;
+        return this;
+    }
+    
+    /// Rotates the current quaternion around the y-axis and returns $(I this).
+    Quaternion rotatey(real alpha) {
+        this = yrotation(alpha) * this;
+        return this;
+    }
+    
+    /// Rotates the current quaternion around the z-axis and returns $(I this).
+    Quaternion rotatez(real alpha) {
+        this = zrotation(alpha) * this;
+        return this;
+    }
+    
+    /// Rotates the current quaternion around an axis and returns $(I this).
+    Quaternion rotate_axis(real alpha, Vector!(qt, 3) axis) {
+        this = axis_rotation(alpha, axis) * this;
+        return this;
+    }
+    
+    /// Applies an euler rotation to the current quaternion and returns $(I this).
+    Quaternion rotate_euler(real heading, real attitude, real bank) {
+        this = euler_rotation(heading, attitude, bank) * this;
+        return this;
+    }
+    
+    unittest {
+        assert(quat.xrotation(PI).quaternion[1..4] == [1.0f, 0.0f, 0.0f]);
+        assert(quat.yrotation(PI).quaternion[1..4] == [0.0f, 1.0f, 0.0f]);
+        assert(quat.zrotation(PI).quaternion[1..4] == [0.0f, 0.0f, 1.0f]);
+        assert((quat.xrotation(PI).w == quat.yrotation(PI).w) && (quat.yrotation(PI).w == quat.zrotation(PI).w));
+        //assert(quat.rotatex(PI).w == to!(quat.qt)(cos(PI)));
+        assert(quat.xrotation(PI).quaternion == quat.identity.rotatex(PI).quaternion);
+        assert(quat.yrotation(PI).quaternion == quat.identity.rotatey(PI).quaternion);
+        assert(quat.zrotation(PI).quaternion == quat.identity.rotatez(PI).quaternion);
+        
+        assert(quat.axis_rotation(PI, vec3(1.0f, 1.0f, 1.0f)).quaternion[1..4] == [1.0f, 1.0f, 1.0f]);
+        assert(quat.axis_rotation(PI, vec3(1.0f, 1.0f, 1.0f)).w == quat.xrotation(PI).w);
+        assert(quat.axis_rotation(PI, vec3(1.0f, 1.0f, 1.0f)).quaternion ==
+            quat.identity.rotate_axis(PI, vec3(1.0f, 1.0f, 1.0f)).quaternion);
+        
+        quat q1 = quat.euler_rotation(PI, PI, PI);
+        assert((q1.x > -1e-16) && (q1.x < 1e-16));
+        assert((q1.y > -1e-16) && (q1.y < 1e-16));
+        assert((q1.z > -1e-16) && (q1.z < 1e-16));
+        //assert(q1.w == -1.0f);
+        assert(quat.euler_rotation(PI, PI, PI).quaternion == quat.identity.rotate_euler(PI, PI, PI).quaternion);
+    }
+    
+    Quaternion opBinary(string op : "*")(Quaternion inp) const {
+        Quaternion ret;
+        
+        ret.w = -x * inp.x - y * inp.y - z * inp.z + w * inp.w;
+        ret.x = x * inp.w + y * inp.z - z * inp.y + w * inp.x;
+        ret.y = -x * inp.z + y * inp.w + z * inp.x + w * inp.y;
+        ret.z = x * inp.y - y * inp.x + z * inp.w + w * inp.z;
+        
+        return ret;
+    }
+    
+    auto opBinaryRight(string op, T)(T inp) const if(!is_quaternion!T) {
+        return this.opBinary!(op)(inp);
+    }
+    
+    Quaternion opBinary(string op)(Quaternion inp) const  if((op == "+") || (op == "-")) {
+        Quaternion ret;
+        
+        mixin("ret.w = w" ~ op ~ "inp.w;");
+        mixin("ret.x = x" ~ op ~ "inp.x;");
+        mixin("ret.y = y" ~ op ~ "inp.y;");
+        mixin("ret.z = z" ~ op ~ "inp.z;");
+        
+        return ret;
+    }
+    
+    Vector!(qt, 3) opBinary(string op : "*")(Vector!(qt, 3) inp) const {
+        Vector!(qt, 3) ret;
+        
+        qt ww = w^^2;
+        qt w2 = w * 2;
+        qt wx2 = w2 * x;
+        qt wy2 = w2 * y;
+        qt wz2 = w2 * z;
+        qt xx = x^^2;
+        qt x2 = x * 2;
+        qt xy2 = x2 * y;
+        qt xz2 = x2 * z;
+        qt yy = y^^2;
+        qt yz2 = 2 * y * z;
+        qt zz = z * z;
+        
+        ret.vector =  [ww * inp.x + wy2 * inp.z - wz2 * inp.y + xx * inp.x +
+            xy2 * inp.y + xz2 * inp.z - zz * inp.x - yy * inp.x,
+            xy2 * inp.x + yy * inp.y + yz2 * inp.z + wz2 * inp.x -
+            zz * inp.y + ww * inp.y - wx2 * inp.z - xx * inp.y,
+            xz2 * inp.x + yz2 * inp.y + zz * inp.z - wy2 * inp.x -
+            yy * inp.z + wx2 * inp.y - xx * inp.z + ww * inp.z];
+        
+        return ret;
+    }
+    
+    Quaternion opBinary(string op : "*")(qt inp) const {
+        return Quaternion(w*inp, x*inp, y*inp, z*inp);
+    }
+    
+    void opOpAssign(string op : "*")(Quaternion inp) {
+        qt w2 = -x * inp.x - y * inp.y - z * inp.z + w * inp.w;
+        qt x2 = x * inp.w + y * inp.z - z * inp.y + w * inp.x;
+        qt y2 = -x * inp.z + y * inp.w + z * inp.x + w * inp.y;
+        qt z2 = x * inp.y - y * inp.x + z * inp.w + w * inp.z;
+        w = w2; x = x2; y = y2; z = z2;
+    }
+    
+    void opOpAssign(string op)(Quaternion inp) if((op == "+") || (op == "-")) {
+        mixin("w = w" ~ op ~ "inp.w;");
+        mixin("x = x" ~ op ~ "inp.x;");
+        mixin("y = y" ~ op ~ "inp.y;");
+        mixin("z = z" ~ op ~ "inp.z;");
+    }
+    
+    void opOpAssign(string op : "*")(qt inp) {
+        quaternion[0] *= inp;
+        quaternion[1] *= inp;
+        quaternion[2] *= inp;
+        quaternion[3] *= inp;
+    }
+    
+    unittest {
+        quat q1 = quat.identity;
+        quat q2 = quat(3.0f, 0.0f, 1.0f, 2.0f);
+        quat q3 = quat(3.4f, 0.1f, 1.2f, 2.3f);
+        
+        assert((q1 * q1).quaternion == q1.quaternion);
+        assert((q1 * q2).quaternion == q2.quaternion);
+        assert((q2 * q1).quaternion == q2.quaternion);
+        quat q4 = q3 * q2;
+        assert((q2 * q3).quaternion != q4.quaternion);
+        q3 *= q2;
+        assert(q4.quaternion == q3.quaternion);
+        assert(almost_equal(q4.x, 0.4f));
+        assert(almost_equal(q4.y, 6.8f));
+        assert(almost_equal(q4.z, 13.8f));
+        assert(almost_equal(q4.w, 4.4f));
+        
+        quat q5 = quat(1.0f, 2.0f, 3.0f, 4.0f);
+        quat q6 = quat(3.0f, 1.0f, 6.0f, 2.0f);
+        
+        assert((q5 - q6).quaternion == [-2.0f, 1.0f, -3.0f, 2.0f]);
+        assert((q5 + q6).quaternion == [4.0f, 3.0f, 9.0f, 6.0f]);
+        assert((q6 - q5).quaternion == [2.0f, -1.0f, 3.0f, -2.0f]);
+        assert((q6 + q5).quaternion == [4.0f, 3.0f, 9.0f, 6.0f]);
+        q5 += q6;
+        assert(q5.quaternion == [4.0f, 3.0f, 9.0f, 6.0f]);
+        q6 -= q6;
+        assert(q6.quaternion == [0.0f, 0.0f, 0.0f, 0.0f]);
+        
+        quat q7 = quat(2.0f, 2.0f, 2.0f, 2.0f);
+        assert((q7 * 2).quaternion == [4.0f, 4.0f, 4.0f, 4.0f]);
+        assert((2 * q7).quaternion == (q7 * 2).quaternion);
+        q7 *= 2;
+        assert(q7.quaternion == [4.0f, 4.0f, 4.0f, 4.0f]);
+        
+        vec3 v1 = vec3(1.0f, 2.0f, 3.0f);
+        assert((q1 * v1).vector == v1.vector);
+        assert((v1 * q1).vector == (q1 * v1).vector);
+        assert((q2 * v1).vector == [-2.0f, 36.0f, 38.0f]);
+    }
+    
+    int opCmp(ref const Quaternion qua) const {
+        foreach(i, a; quaternion) {
+            if(a < qua.quaternion[i]) {
+                return -1;
+            } else if(a > qua.quaternion[i]) {
+                return 1;
+            }
+        }
+        
+        // Quaternions are the same
+        return 0;
+    }
+    
+    bool opEquals(const Quaternion qu) const {
+        return quaternion == qu.quaternion;
+    }
+    
+    bool opCast(T : bool)() const  {
+        return isFinite;
+    }
+    
+    unittest {
+        assert(quat(1.0f, 2.0f, 3.0f, 4.0f) == quat(1.0f, 2.0f, 3.0f, 4.0f));
+        assert(quat(1.0f, 2.0f, 3.0f, 4.0f) != quat(1.0f, 2.0f, 3.0f, 3.0f));
+        
+        assert(!(quat(float.nan, float.nan, float.nan, float.nan)));
+        if(quat(1.0f, 1.0f, 1.0f, 1.0f)) { }
+        else { assert(false); }
+    }
+    
+}
+
+/// Pre-defined quaternion of type float.
+alias Quaternion!(float) quat;
diff --git a/gl3n/linalg/vec.d b/gl3n/linalg/vec.d
new file mode 100644
index 0000000..a7eee21
--- /dev/null
+++ b/gl3n/linalg/vec.d
@@ -0,0 +1,734 @@
+/**
+gl3n.linalg.vec
+
+Special thanks to:
+$(UL
+  $(LI Tomasz Stachowiak (h3r3tic): allowed me to use parts of $(LINK2 https://bitbucket.org/h3r3tic/boxen/src/default/src/xf/omg, omg).)
+  $(LI Jakob Øvrum (jA_cOp): improved the code a lot!)
+  $(LI Florian Boesch (___doc__): helps me to understand opengl/complex maths better, see: $(LINK http://codeflow.org/).)
+  $(LI #D on freenode: answered general questions about D.)
+)
+
+Authors: David Herberth, Stephan Dilly
+License: MIT
+
+Note: All methods marked with pure are weakly pure since, they all access an instance member.
+All static methods are strongly pure.
+*/
+
+module gl3n.linalg.vec;
+
+import gl3n.linalg.matrix;
+import gl3n.util : is_vector, is_matrix, is_quaternion, TupleRange;
+import gl3n.math : sqrt, max, almost_equal;
+
+import std.math : isNaN, isInfinity;
+import std.string : format;
+import std.conv : to;
+import std.traits : isArray, isDynamicArray, isStaticArray, isIntegral;
+
+/// Base template for all vector-types.
+/// Params:
+/// type = all values get stored as this type
+/// dimension = specifies the dimension of the vector, can be 1, 2, 3 or 4
+/// Examples:
+/// ---
+/// alias Vector!(int, 3) vec3i;
+/// alias Vector!(float, 4) vec4;
+/// alias Vector!(real, 2) vec2r;
+/// ---
+struct Vector(type, int dimension_) {
+    static assert(dimension > 0, "0 dimensional vectors don't exist.");
+    
+    alias type vt; /// Holds the internal type of the vector.
+    static const int dimension = dimension_; ///Holds the dimension of the vector.
+    
+    vt[dimension] vector; /// Holds all coordinates, length conforms dimension.
+    
+    /// Returns a pointer to the coordinates.
+    @property auto value_ptr() const { return vector.ptr; }
+    
+    /// Returns the current vector formatted as string, useful for printing the vector.
+    @property string as_string() {
+        return format("%s", vector);
+    }
+    alias as_string toString; /// ditto
+    
+    @safe pure nothrow:
+    ///
+    private @property ref inout(vt) get_(char coord)() inout {
+        return vector[coord_to_index!coord];
+    }
+    
+    alias get_!'x' x; /// static properties to access the values.
+    alias x u; /// ditto
+    alias x s; /// ditto
+    alias x r; /// ditto
+    static if(dimension >= 2) {
+        alias get_!'y' y; /// ditto
+        alias y v; /// ditto
+        alias y t; /// ditto
+        alias y g; /// ditto
+    }
+    static if(dimension >= 3) {
+        alias get_!'z' z; /// ditto
+        alias z b; /// ditto
+        alias z p; /// ditto
+    }
+    static if(dimension >= 4) {
+        alias get_!'w' w; /// ditto
+        alias w a; /// ditto
+        alias w q; /// ditto
+    }
+    
+    static if(dimension == 2) {
+        enum Vector e1 = Vector(1.to!vt, 0.to!vt); /// canonical basis for Euclidian space
+        enum Vector e2 = Vector(0.to!vt, 1.to!vt); /// ditto
+    } else static if(dimension == 3) {
+        enum Vector e1 = Vector(1.to!vt, 0.to!vt, 0.to!vt); /// canonical basis for Euclidian space
+        enum Vector e2 = Vector(0.to!vt, 1.to!vt, 0.to!vt); /// ditto
+        enum Vector e3 = Vector(0.to!vt, 0.to!vt, 1.to!vt); /// ditto
+    } else static if(dimension == 4) {
+        enum Vector e1 = Vector(1.to!vt, 0.to!vt, 0.to!vt, 0.to!vt); /// canonical basis for Euclidian space
+        enum Vector e2 = Vector(0.to!vt, 1.to!vt, 0.to!vt, 0.to!vt); /// ditto
+        enum Vector e3 = Vector(0.to!vt, 0.to!vt, 1.to!vt, 0.to!vt); /// ditto
+        enum Vector e4 = Vector(0.to!vt, 0.to!vt, 0.to!vt, 1.to!vt); /// ditto
+    }
+    
+    unittest {
+        assert(vec2.e1.vector == [1.0, 0.0]);
+        assert(vec2.e2.vector == [0.0, 1.0]);
+        
+        assert(vec3.e1.vector == [1.0, 0.0, 0.0]);
+        assert(vec3.e2.vector == [0.0, 1.0, 0.0]);
+        assert(vec3.e3.vector == [0.0, 0.0, 1.0]);
+        
+        assert(vec4.e1.vector == [1.0, 0.0, 0.0, 0.0]);
+        assert(vec4.e2.vector == [0.0, 1.0, 0.0, 0.0]);
+        assert(vec4.e3.vector == [0.0, 0.0, 1.0, 0.0]);
+        assert(vec4.e4.vector == [0.0, 0.0, 0.0, 1.0]);
+    }
+    
+    static void isCompatibleVectorImpl(int d)(Vector!(vt, d) vec) if(d <= dimension) {
+    }
+    
+    template isCompatibleVector(T) {
+        enum isCompatibleVector = is(typeof(isCompatibleVectorImpl(T.init)));
+    }
+    
+    static void isCompatibleMatrixImpl(int r, int c)(Matrix!(vt, r, c) m) {
+    }
+    
+    template isCompatibleMatrix(T) {
+        enum isCompatibleMatrix = is(typeof(isCompatibleMatrixImpl(T.init)));
+    }
+    
+    private void construct(int i, T, Tail...)(T head, Tail tail) {
+        static if(i >= dimension) {
+            static assert(false, "Too many arguments passed to constructor");
+        } else static if(is(T : vt)) {
+            vector[i] = head;
+            construct!(i + 1)(tail);
+        } else static if(isDynamicArray!T) {
+            static assert((Tail.length == 0) && (i == 0), "dynamic array can not be passed together with other arguments");
+            vector[] = head[];
+        } else static if(isStaticArray!T) {
+            vector[i .. i + T.length] = head[];
+            construct!(i + T.length)(tail);
+        } else static if(isCompatibleVector!T) {
+            vector[i .. i + T.dimension] = head.vector[];
+            construct!(i + T.dimension)(tail);
+        } else {
+            static assert(false, "Vector constructor argument must be of type " ~ vt.stringof ~ " or Vector, not " ~ T.stringof);
+        }
+    }
+    
+    private void construct(int i)() { // terminate
+        static assert(i == dimension, "Not enough arguments passed to constructor");
+    }
+    
+    /// Constructs the vector.
+    /// If a single value is passed the vector, the vector will be cleared with this value.
+    /// If a vector with a higher dimension is passed the vector will hold the first values up to its dimension.
+    /// If mixed types are passed they will be joined together (allowed types: vector, static array, $(I vt)).
+    /// Examples:
+    /// ---
+    /// vec4 v4 = vec4(1.0f, vec2(2.0f, 3.0f), 4.0f);
+    /// vec3 v3 = vec3(v4); // v3 = vec3(1.0f, 2.0f, 3.0f);
+    /// vec2 v2 = v3.xy; // swizzling returns a static array.
+    /// vec3 v3_2 = vec3(1.0f); // vec3 v3_2 = vec3(1.0f, 1.0f, 1.0f);
+    /// ---
+    this(Args...)(Args args) {
+        construct!(0)(args);
+    }
+    
+    /// ditto
+    this(T)(T vec) if(is_vector!T && is(T.vt : vt) && (T.dimension >= dimension)) {
+        foreach(i; TupleRange!(0, dimension)) {
+            vector[i] = vec.vector[i];
+        }
+    }
+    
+    /// ditto
+    this()(vt value) {
+        clear(value);
+    }
+    
+    /// Returns true if all values are not nan and finite, otherwise false.
+    @property bool isFinite() const {
+        static if(isIntegral!type) {
+            return true;
+        }
+        else {
+            foreach(v; vector) {
+                if(isNaN(v) || isInfinity(v)) {
+                    return false;
+                }
+            }
+            return true;
+        }
+    }
+    deprecated("Use isFinite instead of ok") alias ok = isFinite;
+    
+    /// Sets all values of the vector to value.
+    void clear(vt value) {
+        foreach(i; TupleRange!(0, dimension)) {
+            vector[i] = value;
+        }
+    }
+    
+    unittest {
+        vec3 vec_clear;
+        assert(!vec_clear.isFinite);
+        vec_clear.clear(1.0f);
+        assert(vec_clear.isFinite);
+        assert(vec_clear.vector == [1.0f, 1.0f, 1.0f]);
+        assert(vec_clear.vector == vec3(1.0f).vector);
+        vec_clear.clear(float.infinity);
+        assert(!vec_clear.isFinite);
+        vec_clear.clear(float.nan);
+        assert(!vec_clear.isFinite);
+        vec_clear.clear(1.0f);
+        assert(vec_clear.isFinite);
+        
+        vec4 b = vec4(1.0f, vec_clear);
+        assert(b.isFinite);
+        assert(b.vector == [1.0f, 1.0f, 1.0f, 1.0f]);
+        assert(b.vector == vec4(1.0f).vector);
+        
+        vec2 v2_1 = vec2(vec2(0.0f, 1.0f));
+        assert(v2_1.vector == [0.0f, 1.0f]);
+        
+        vec2 v2_2 = vec2(1.0f, 1.0f);
+        assert(v2_2.vector == [1.0f, 1.0f]);
+        
+        vec3 v3 = vec3(v2_1, 2.0f);
+        assert(v3.vector == [0.0f, 1.0f, 2.0f]);
+        
+        vec4 v4_1 = vec4(1.0f, vec2(2.0f, 3.0f), 4.0f);
+        assert(v4_1.vector == [1.0f, 2.0f, 3.0f, 4.0f]);
+        assert(vec3(v4_1).vector == [1.0f, 2.0f, 3.0f]);
+        assert(vec2(vec3(v4_1)).vector == [1.0f, 2.0f]);
+        assert(vec2(vec3(v4_1)).vector == vec2(v4_1).vector);
+        assert(v4_1.vector == vec4([1.0f, 2.0f, 3.0f, 4.0f]).vector);
+        
+        vec4 v4_2 = vec4(vec2(1.0f, 2.0f), vec2(3.0f, 4.0f));
+        assert(v4_2.vector == [1.0f, 2.0f, 3.0f, 4.0f]);
+        assert(vec3(v4_2).vector == [1.0f, 2.0f, 3.0f]);
+        assert(vec2(vec3(v4_2)).vector == [1.0f, 2.0f]);
+        assert(vec2(vec3(v4_2)).vector == vec2(v4_2).vector);
+        assert(v4_2.vector == vec4([1.0f, 2.0f, 3.0f, 4.0f]).vector);
+        
+        float[2] f2 = [1.0f, 2.0f];
+        float[3] f3 = [1.0f, 2.0f, 3.0f];
+        float[4] f4 = [1.0f, 2.0f, 3.0f, 4.0f];
+        assert(vec2(1.0f, 2.0f).vector == vec2(f2).vector);
+        assert(vec3(1.0f, 2.0f, 3.0f).vector == vec3(f3).vector);
+        assert(vec3(1.0f, 2.0f, 3.0f).vector == vec3(f2, 3.0f).vector);
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f).vector == vec4(f4).vector);
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f).vector == vec4(f3, 4.0f).vector);
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f).vector == vec4(f2, 3.0f, 4.0f).vector);
+        // useful for: "vec4 v4 = […]" or "vec4 v4 = other_vector.rgba"
+        
+        assert(vec3(vec3i(1, 2, 3)) == vec3(1.0, 2.0, 3.0));
+        assert(vec3d(vec3(1.0, 2.0, 3.0)) == vec3d(1.0, 2.0, 3.0));
+        
+        static assert(!__traits(compiles, vec3(0.0f, 0.0f)));
+        static assert(!__traits(compiles, vec4(0.0f, 0.0f, 0.0f)));
+        static assert(!__traits(compiles, vec4(0.0f, vec2(0.0f, 0.0f))));
+        static assert(!__traits(compiles, vec4(vec3(0.0f, 0.0f, 0.0f))));
+    }
+    
+    template coord_to_index(char c) {
+        static if((c == 'x') || (c == 'r') || (c == 'u') || (c == 's')) {
+            enum coord_to_index = 0;
+        } else static if((c == 'y') || (c == 'g') || (c == 'v') || (c == 't')) {
+            enum coord_to_index = 1;
+        } else static if((c == 'z') || (c == 'b') || (c == 'p')) {
+            static assert(dimension >= 3, "the " ~ c ~ " property is only available on vectors with a third dimension.");
+            enum coord_to_index = 2;
+        } else static if((c == 'w') || (c == 'a') || (c == 'q')) {
+            static assert(dimension >= 4, "the " ~ c ~ " property is only available on vectors with a fourth dimension.");
+            enum coord_to_index = 3;
+        } else {
+            static assert(false, "accepted coordinates are x, s, r, u, y, g, t, v, z, p, b, w, q and a not " ~ c ~ ".");
+        }
+    }
+    
+    static if(dimension == 2) { void set(vt x, vt y) { vector[0] = x; vector[1] = y; } }
+    static if(dimension == 3) { void set(vt x, vt y, vt z) { vector[0] = x; vector[1] = y; vector[2] = z; } }
+    static if(dimension == 4) { void set(vt x, vt y, vt z, vt w) { vector[0] = x; vector[1] = y; vector[2] = z; vector[3] = w; } }
+    
+    /// Updates the vector with the values from other.
+    void update(Vector!(vt, dimension) other) {
+        vector = other.vector;
+    }
+    
+    unittest {
+        vec2 v2 = vec2(1.0f, 2.0f);
+        assert(v2.x == 1.0f);
+        assert(v2.y == 2.0f);
+        v2.x = 3.0f;
+        v2.x += 1;
+        v2.x -= 1;
+        assert(v2.vector == [3.0f, 2.0f]);
+        v2.y = 4.0f;
+        v2.y += 1;
+        v2.y -= 1;
+        assert(v2.vector == [3.0f, 4.0f]);
+        assert((v2.x == 3.0f) && (v2.x == v2.u) && (v2.x == v2.s) && (v2.x == v2.r));
+        assert(v2.y == 4.0f);
+        assert((v2.y == 4.0f) && (v2.y == v2.v) && (v2.y == v2.t) && (v2.y == v2.g));
+        v2.set(0.0f, 1.0f);
+        assert(v2.vector == [0.0f, 1.0f]);
+        v2.update(vec2(3.0f, 4.0f));
+        assert(v2.vector == [3.0f, 4.0f]);
+        
+        vec3 v3 = vec3(1.0f, 2.0f, 3.0f);
+        assert(v3.x == 1.0f);
+        assert(v3.y == 2.0f);
+        assert(v3.z == 3.0f);
+        v3.x = 3.0f;
+        v3.x += 1;
+        v3.x -= 1;
+        assert(v3.vector == [3.0f, 2.0f, 3.0f]);
+        v3.y = 4.0f;
+        v3.y += 1;
+        v3.y -= 1;
+        assert(v3.vector == [3.0f, 4.0f, 3.0f]);
+        v3.z = 5.0f;
+        v3.z += 1;
+        v3.z -= 1;
+        assert(v3.vector == [3.0f, 4.0f, 5.0f]);
+        assert((v3.x == 3.0f) && (v3.x == v3.s) && (v3.x == v3.r));
+        assert((v3.y == 4.0f) && (v3.y == v3.t) && (v3.y == v3.g));
+        assert((v3.z == 5.0f) && (v3.z == v3.p) && (v3.z == v3.b));
+        v3.set(0.0f, 1.0f, 2.0f);
+        assert(v3.vector == [0.0f, 1.0f, 2.0f]);
+        v3.update(vec3(3.0f, 4.0f, 5.0f));
+        assert(v3.vector == [3.0f, 4.0f, 5.0f]);
+        
+        vec4 v4 = vec4(1.0f, 2.0f, vec2(3.0f, 4.0f));
+        assert(v4.x == 1.0f);
+        assert(v4.y == 2.0f);
+        assert(v4.z == 3.0f);
+        assert(v4.w == 4.0f);
+        v4.x = 3.0f;
+        v4.x += 1;
+        v4.x -= 1;
+        assert(v4.vector == [3.0f, 2.0f, 3.0f, 4.0f]);
+        v4.y = 4.0f;
+        v4.y += 1;
+        v4.y -= 1;
+        assert(v4.vector == [3.0f, 4.0f, 3.0f, 4.0f]);
+        v4.z = 5.0f;
+        v4.z += 1;
+        v4.z -= 1;
+        assert(v4.vector == [3.0f, 4.0f, 5.0f, 4.0f]);
+        v4.w = 6.0f;
+        v4.w += 1;
+        v4.w -= 1;
+        assert(v4.vector == [3.0f, 4.0f, 5.0f, 6.0f]);
+        assert((v4.x == 3.0f) && (v4.x == v4.s) && (v4.x == v4.r));
+        assert((v4.y == 4.0f) && (v4.y == v4.t) && (v4.y == v4.g));
+        assert((v4.z == 5.0f) && (v4.z == v4.p) && (v4.z == v4.b));
+        assert((v4.w == 6.0f) && (v4.w == v4.q) && (v4.w == v4.a));
+        v4.set(0.0f, 1.0f, 2.0f, 3.0f);
+        assert(v4.vector == [0.0f, 1.0f, 2.0f, 3.0f]);
+        v4.update(vec4(3.0f, 4.0f, 5.0f, 6.0f));
+        assert(v4.vector == [3.0f, 4.0f, 5.0f, 6.0f]);
+    }
+    
+    private void dispatchImpl(int i, string s, int size)(ref vt[size] result) const {
+        static if(s.length > 0) {
+            result[i] = vector[coord_to_index!(s[0])];
+            dispatchImpl!(i + 1, s[1..$])(result);
+        }
+    }
+    
+    /// Implements dynamic swizzling.
+    /// Returns: a Vector
+    @property Vector!(vt, s.length) opDispatch(string s)() const {
+        vt[s.length] ret;
+        dispatchImpl!(0, s)(ret);
+        Vector!(vt, s.length) ret_vec;
+        ret_vec.vector = ret;
+        return ret_vec;
+    }
+    
+    unittest {
+        vec2 v2 = vec2(1.0f, 2.0f);
+        assert(v2.xytsy == [1.0f, 2.0f, 2.0f, 1.0f, 2.0f]);
+        
+        assert(vec3(1.0f, 2.0f, 3.0f).xybzyr == [1.0f, 2.0f, 3.0f, 3.0f, 2.0f, 1.0f]);
+        assert(vec4(v2, 3.0f, 4.0f).xyzwrgbastpq == [1.0f, 2.0f, 3.0f, 4.0f,
+                1.0f, 2.0f, 3.0f, 4.0f,
+                1.0f, 2.0f, 3.0f, 4.0f]);
+        assert(vec4(v2, 3.0f, 4.0f).wgyzax == [4.0f, 2.0f, 2.0f, 3.0f, 4.0f, 1.0f]);
+        assert(vec4(v2.xyst).vector == [1.0f, 2.0f, 1.0f, 2.0f]);
+    }
+    
+    /// Returns the squared magnitude of the vector.
+    @property real magnitude_squared() const {
+        real temp = 0;
+        
+        foreach(index; TupleRange!(0, dimension)) {
+            temp += vector[index]^^2;
+        }
+        
+        return temp;
+    }
+    
+    /// Returns the magnitude of the vector.
+    @property real magnitude() const {
+        return sqrt(magnitude_squared);
+    }
+    
+    alias magnitude_squared length_squared; /// ditto
+    alias magnitude length; /// ditto
+    
+    /// Normalizes the vector.
+    void normalize() {
+        real len = length;
+        
+        if(len != 0) {
+            foreach(index; TupleRange!(0, dimension)) {
+                vector[index] = cast(type)(vector[index]/len);
+            }
+        }
+    }
+    
+    /// Returns a normalized copy of the current vector.
+    @property Vector normalized() const {
+        Vector ret;
+        ret.update(this);
+        ret.normalize();
+        return ret;
+    }
+    
+    Vector opUnary(string op : "-")() const {
+        Vector ret;
+        
+        foreach(index; TupleRange!(0, dimension)) {
+            ret.vector[index] = -vector[index];
+        }
+        
+        return ret;
+    }
+    
+    unittest {
+        assert(vec2(1.0f, 1.0f) == -vec2(-1.0f, -1.0f));
+        assert(vec2(-1.0f, 1.0f) == -vec2(1.0f, -1.0f));
+        
+        assert(-vec3(1.0f, 1.0f, 1.0f) == vec3(-1.0f, -1.0f, -1.0f));
+        assert(-vec3(-1.0f, 1.0f, -1.0f) == vec3(1.0f, -1.0f, 1.0f));
+        
+        assert(vec4(1.0f, 1.0f, 1.0f, 1.0f) == -vec4(-1.0f, -1.0f, -1.0f, -1.0f));
+        assert(vec4(-1.0f, 1.0f, -1.0f, 1.0f) == -vec4(1.0f, -1.0f, 1.0f, -1.0f));
+    }
+    
+    // let the math begin!
+    Vector opBinary(string op : "*")(vt r) const {
+        Vector ret;
+        
+        foreach(index; TupleRange!(0, dimension)) {
+            ret.vector[index] = vector[index] * r;
+        }
+        
+        return ret;
+    }
+    
+    Vector opBinary(string op : "/")(vt r) const {
+        Vector ret;
+        
+        foreach(index; TupleRange!(0, dimension)) {
+            ret.vector[index] = vector[index] / r;
+        }
+        
+        return ret;
+    }
+    
+    Vector opBinary(string op)(Vector r) const if((op == "+") || (op == "-")) {
+        Vector ret;
+        
+        foreach(index; TupleRange!(0, dimension)) {
+            ret.vector[index] = mixin("vector[index]" ~ op ~ "r.vector[index]");
+        }
+        
+        return ret;
+    }
+    
+    vt opBinary(string op : "*")(Vector r) const {
+        return dot(this, r);
+    }
+    
+    // vector * matrix (for matrix * vector -> struct Matrix)
+    Vector!(vt, T.cols) opBinary(string op : "*", T)(T inp) const if(isCompatibleMatrix!T && (T.rows == dimension)) {
+        Vector!(vt, T.cols) ret;
+        ret.clear(0);
+        
+        foreach(c; TupleRange!(0, T.cols)) {
+            foreach(r; TupleRange!(0, T.rows)) {
+                ret.vector[c] += vector[r] * inp.matrix[r][c];
+            }
+        }
+        
+        return ret;
+    }
+    
+    auto opBinaryRight(string op, T)(T inp) const if(!is_vector!T && !is_matrix!T && !is_quaternion!T) {
+        return this.opBinary!(op)(inp);
+    }
+    
+    unittest {
+        import gl3n.linalg.matrix : mat2, mat3;
+
+        vec2 v2 = vec2(1.0f, 3.0f);
+        auto v2times2 = 2 * v2;
+        assert((v2*2.5f).vector == [2.5f, 7.5f]);
+        assert((v2+vec2(3.0f, 1.0f)).vector == [4.0f, 4.0f]);
+        assert((v2-vec2(1.0f, 3.0f)).vector == [0.0f, 0.0f]);
+        assert((v2*vec2(2.0f, 2.0f)) == 8.0f);
+        
+        vec3 v3 = vec3(1.0f, 3.0f, 5.0f);
+        assert((v3*2.5f).vector == [2.5f, 7.5f, 12.5f]);
+        assert((v3+vec3(3.0f, 1.0f, -1.0f)).vector == [4.0f, 4.0f, 4.0f]);
+        assert((v3-vec3(1.0f, 3.0f, 5.0f)).vector == [0.0f, 0.0f, 0.0f]);
+        assert((v3*vec3(2.0f, 2.0f, 2.0f)) == 18.0f);
+        
+        vec4 v4 = vec4(1.0f, 3.0f, 5.0f, 7.0f);
+        assert((v4*2.5f).vector == [2.5f, 7.5f, 12.5f, 17.5]);
+        assert((v4+vec4(3.0f, 1.0f, -1.0f, -3.0f)).vector == [4.0f, 4.0f, 4.0f, 4.0f]);
+        assert((v4-vec4(1.0f, 3.0f, 5.0f, 7.0f)).vector == [0.0f, 0.0f, 0.0f, 0.0f]);
+        assert((v4*vec4(2.0f, 2.0f, 2.0f, 2.0f)) == 32.0f);
+        
+        mat2 m2 = mat2(1.0f, 2.0f, 3.0f, 4.0f);
+        vec2 v2_2 = vec2(2.0f, 2.0f);
+        assert((v2_2*m2).vector == [8.0f, 12.0f]);
+        
+        mat3 m3 = mat3(1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f);
+        vec3 v3_2 = vec3(2.0f, 2.0f, 2.0f);
+        assert((v3_2*m3).vector == [24.0f, 30.0f, 36.0f]);
+    }
+    
+    void opOpAssign(string op : "*")(vt r) {
+        foreach(index; TupleRange!(0, dimension)) {
+            vector[index] *= r;
+        }
+    }
+    
+    void opOpAssign(string op : "/")(vt r) {
+        foreach(index; TupleRange!(0, dimension)) {
+            vector[index] /= r;
+        }
+    }
+    
+    void opOpAssign(string op)(Vector r) if((op == "+") || (op == "-")) {
+        foreach(index; TupleRange!(0, dimension)) {
+            mixin("vector[index]" ~ op ~ "= r.vector[index];");
+        }
+    }
+    
+    unittest {
+        vec2 v2 = vec2(1.0f, 3.0f);
+        v2 *= 2.5f;
+        assert(v2.vector == [2.5f, 7.5f]);
+        v2 -= vec2(2.5f, 7.5f);
+        assert(v2.vector == [0.0f, 0.0f]);
+        v2 += vec2(1.0f, 3.0f);
+        assert(v2.vector == [1.0f, 3.0f]);
+        assert(almost_equal(v2.length, sqrt(10.0f)));
+        assert(v2.length_squared == 10.0f);
+        assert((v2.magnitude == v2.length) && (v2.magnitude_squared == v2.length_squared));
+        v2 /= 2.0f;
+        assert(v2.vector == [0.5f, 1.5f]);
+        assert(almost_equal(v2.normalized, vec2(1.0f/sqrt(10.0f), 3.0f/sqrt(10.0f))));
+        
+        vec3 v3 = vec3(1.0f, 3.0f, 5.0f);
+        v3 *= 2.5f;
+        assert(v3.vector == [2.5f, 7.5f, 12.5f]);
+        v3 -= vec3(2.5f, 7.5f, 12.5f);
+        assert(v3.vector == [0.0f, 0.0f, 0.0f]);
+        v3 += vec3(1.0f, 3.0f, 5.0f);
+        assert(v3.vector == [1.0f, 3.0f, 5.0f]);
+        assert(almost_equal(v3.length, sqrt(35.0f)));
+        assert(v3.length_squared == 35.0f);
+        assert((v3.magnitude == v3.length) && (v3.magnitude_squared == v3.length_squared));
+        v3 /= 2.0f;
+        assert(v3.vector == [0.5f, 1.5f, 2.5f]);
+        assert(almost_equal(v3.normalized, vec3(1.0f/sqrt(35.0f), 3.0f/sqrt(35.0f), 5.0f/sqrt(35.0f))));
+        
+        vec4 v4 = vec4(1.0f, 3.0f, 5.0f, 7.0f);
+        v4 *= 2.5f;
+        assert(v4.vector == [2.5f, 7.5f, 12.5f, 17.5]);
+        v4 -= vec4(2.5f, 7.5f, 12.5f, 17.5f);
+        assert(v4.vector == [0.0f, 0.0f, 0.0f, 0.0f]);
+        v4 += vec4(1.0f, 3.0f, 5.0f, 7.0f);
+        assert(v4.vector == [1.0f, 3.0f, 5.0f, 7.0f]);
+        assert(almost_equal(v4.length, sqrt(84.0f)));
+        assert(v4.length_squared == 84.0f);
+        assert((v4.magnitude == v4.length) && (v4.magnitude_squared == v4.length_squared));
+        v4 /= 2.0f;
+        assert(v4.vector == [0.5f, 1.5f, 2.5f, 3.5f]);
+        assert(almost_equal(v4.normalized, vec4(1.0f/sqrt(84.0f), 3.0f/sqrt(84.0f), 5.0f/sqrt(84.0f), 7.0f/sqrt(84.0f))));
+    }
+    
+    int opCmp(ref const Vector vec) const {
+        foreach(i, a; vector) {
+            if(a < vec.vector[i]) {
+                return -1;
+            } else if(a > vec.vector[i]) {
+                return 1;
+            }
+        }
+        
+        // Vectors are the same
+        return 0;
+    }
+    
+    bool opEquals(T)(const T vec) const if(!isArray!T && T.dimension == dimension) {
+        return vector == vec.vector;
+    }
+    
+    bool opEquals(T)(const(T)[] array) const if(!isArray!T && !is_vector!T) {
+        if(array.length != dimension) {
+            return false;
+        }
+        
+        foreach(index; TupleRange!(0, dimension)) {
+            if(vector[index] != array[index]) {
+                return false;
+            }
+        }
+        
+        return true;
+    }
+    
+    bool opCast(T : bool)() const {
+        return isFinite;
+    }
+    
+    unittest {
+        assert(vec2(1.0f, 2.0f) == vec2(1.0f, 2.0f));
+        assert(vec2(1.0f, 2.0f) != vec2(1.0f, 1.0f));
+        assert(vec2(1.0f, 2.0f) == vec2d(1.0, 2.0));
+        assert(vec2(1.0f, 2.0f) != vec2d(1.0, 1.0));
+        assert(vec2(1.0f, 2.0f) == vec2(1.0f, 2.0f).vector);
+        assert(vec2(1.0f, 2.0f) != vec2(1.0f, 1.0f).vector);
+        assert(vec2(1.0f, 2.0f) == vec2d(1.0, 2.0).vector);
+        assert(vec2(1.0f, 2.0f) != vec2d(1.0, 1.0).vector);
+        
+        assert(vec3(1.0f, 2.0f, 3.0f) == vec3(1.0f, 2.0f, 3.0f));
+        assert(vec3(1.0f, 2.0f, 3.0f) != vec3(1.0f, 2.0f, 2.0f));
+        assert(vec3(1.0f, 2.0f, 3.0f) == vec3d(1.0, 2.0, 3.0));
+        assert(vec3(1.0f, 2.0f, 3.0f) != vec3d(1.0, 2.0, 2.0));
+        assert(vec3(1.0f, 2.0f, 3.0f) == vec3(1.0f, 2.0f, 3.0f).vector);
+        assert(vec3(1.0f, 2.0f, 3.0f) != vec3(1.0f, 2.0f, 2.0f).vector);
+        assert(vec3(1.0f, 2.0f, 3.0f) == vec3d(1.0, 2.0, 3.0).vector);
+        assert(vec3(1.0f, 2.0f, 3.0f) != vec3d(1.0, 2.0, 2.0).vector);
+        
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4(1.0f, 2.0f, 3.0f, 4.0f));
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4(1.0f, 2.0f, 3.0f, 3.0f));
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4d(1.0, 2.0, 3.0, 4.0));
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4d(1.0, 2.0, 3.0, 3.0));
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4(1.0f, 2.0f, 3.0f, 4.0f).vector);
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4(1.0f, 2.0f, 3.0f, 3.0f).vector);
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) == vec4d(1.0, 2.0, 3.0, 4.0).vector);
+        assert(vec4(1.0f, 2.0f, 3.0f, 4.0f) != vec4d(1.0, 2.0, 3.0, 3.0).vector);
+        
+        assert(!(vec4(float.nan)));
+        if(vec4(1.0f)) { }
+        else { assert(false); }
+    }
+    
+}
+
+/// Calculates the product between two vectors.
+T.vt dot(T)(const T veca, const T vecb) @safe pure nothrow if(is_vector!T) {
+    T.vt temp = 0;
+    
+    foreach(index; TupleRange!(0, T.dimension)) {
+        temp += veca.vector[index] * vecb.vector[index];
+    }
+    
+    return temp;
+}
+
+/// Calculates the cross product of two 3-dimensional vectors.
+T cross(T)(const T veca, const T vecb) @safe pure nothrow if(is_vector!T && (T.dimension == 3)) {
+    return T(veca.y * vecb.z - vecb.y * veca.z,
+        veca.z * vecb.x - vecb.z * veca.x,
+        veca.x * vecb.y - vecb.x * veca.y);
+}
+
+/// Calculates the distance between two vectors.
+T.vt distance(T)(const T veca, const T vecb) @safe pure nothrow if(is_vector!T) {
+    return (veca - vecb).length;
+}
+
+unittest {
+    // dot is already tested in Vector.opBinary, so no need for testing with more vectors
+    vec3 v1 = vec3(1.0f, 2.0f, -3.0f);
+    vec3 v2 = vec3(1.0f, 3.0f, 2.0f);
+    
+    assert(dot(v1, v2) == 1.0f);
+    assert(dot(v1, v2) == (v1 * v2));
+    assert(dot(v1, v2) == dot(v2, v1));
+    assert((v1 * v2) == (v1 * v2));
+    
+    assert(cross(v1, v2).vector == [13.0f, -5.0f, 1.0f]);
+    assert(cross(v2, v1).vector == [-13.0f, 5.0f, -1.0f]);
+    
+    assert(distance(vec2(0.0f, 0.0f), vec2(0.0f, 10.0f)) == 10.0);
+}
+
+/// reflect a vector using a surface normal
+T reflect(T)(const T vec, const T norm) @safe pure nothrow if(is_vector!T) {
+    return (2 * (vec * norm) * norm) - vec;
+}
+
+unittest
+{
+    assert(vec2(1,1).reflect(vec2(0,1)) == vec2(-1,1));
+    assert(vec2(-1,1).reflect(vec2(0,1)) == vec2(1,1));
+    assert(vec2(2,1).reflect(vec2(0,1)) == vec2(-2,1));
+    
+    assert(vec3(1,1,1).reflect(vec3(0,1,0)) == vec3(-1,1,-1));
+}
+
+/// Pre-defined vector types, the number represents the dimension and the last letter the type (none = float, d = double, i = int).
+alias Vector!(float, 2) vec2;
+alias Vector!(float, 3) vec3; /// ditto
+alias Vector!(float, 4) vec4; /// ditto
+
+alias Vector!(double, 2) vec2d; /// ditto
+alias Vector!(double, 3) vec3d; /// ditto
+alias Vector!(double, 4) vec4d; /// ditto
+
+alias Vector!(int, 2) vec2i; /// ditto
+alias Vector!(int, 3) vec3i; /// ditto
+alias Vector!(int, 4) vec4i; /// ditto
+
+/*alias Vector!(ubyte, 2) vec2ub;
+alias Vector!(ubyte, 3) vec3ub;
+alias Vector!(ubyte, 4) vec4ub;*/
diff --git a/libgl3n.a b/libgl3n.a
new file mode 100644
index 0000000..0731d89
Binary files /dev/null and b/libgl3n.a differ