dep -> govendor

vendor/github.com/go-gl/mathgl/mgl64/conv.go

@@ -0,0 +1,96 @@
+// This file is generated from mgl32/conv.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mgl64
+
+import (
+	"math"
+)
+
+// Converts 3-dimensional cartesian coordinates (x,y,z) to spherical
+// coordinates with radius r, inclination theta, and azimuth phi.
+//
+// All angles are in radians.
+func CartesianToSpherical(coord Vec3) (r, theta, phi float64) {
+	r = coord.Len()
+	theta = float64(math.Acos(float64(coord[2] / r)))
+	phi = float64(math.Atan2(float64(coord[1]), float64(coord[0])))
+
+	return
+}
+
+// Converts 3-dimensional cartesian coordinates (x,y,z) to cylindrical
+// coordinates with radial distance r, azimuth phi, and height z.
+//
+// All angles are in radians.
+func CartesianToCylindical(coord Vec3) (rho, phi, z float64) {
+	rho = float64(math.Hypot(float64(coord[0]), float64(coord[1])))
+
+	phi = float64(math.Atan2(float64(coord[1]), float64(coord[0])))
+
+	z = coord[2]
+
+	return
+}
+
+// Converts spherical coordinates with radius r, inclination theta,
+// and azimuth phi to cartesian coordinates (x,y,z).
+//
+// Angles are in radians.
+func SphericalToCartesian(r, theta, phi float64) Vec3 {
+	st, ct := math.Sincos(float64(theta))
+	sp, cp := math.Sincos(float64(phi))
+
+	return Vec3{r * float64(st*cp), r * float64(st*sp), r * float64(ct)}
+}
+
+// Converts spherical coordinates with radius r, inclination theta,
+// and azimuth phi to cylindrical coordinates with radial distance r,
+// azimuth phi, and height z.
+//
+// Angles are in radians
+func SphericalToCylindrical(r, theta, phi float64) (rho, phi2, z float64) {
+	s, c := math.Sincos(float64(theta))
+
+	rho = r * float64(s)
+	z = r * float64(c)
+	phi2 = phi
+
+	return
+}
+
+// Converts cylindrical coordinates with radial distance r,
+// azimuth phi, and height z to spherical coordinates with radius r,
+// inclination theta, and azimuth phi.
+//
+// Angles are in radians
+func CylindircalToSpherical(rho, phi, z float64) (r, theta, phi2 float64) {
+	r = float64(math.Hypot(float64(rho), float64(z)))
+	phi2 = phi
+	theta = float64(math.Atan2(float64(rho), float64(z)))
+
+	return
+}
+
+// Converts cylindrical coordinates with radial distance r,
+// azimuth phi, and height z to cartesian coordinates (x,y,z)
+//
+// Angles are in radians.
+func CylindricalToCartesian(rho, phi, z float64) Vec3 {
+	s, c := math.Sincos(float64(phi))
+
+	return Vec3{rho * float64(c), rho * float64(s), z}
+}
+
+// Converts degrees to radians
+func DegToRad(angle float64) float64 {
+	return angle * float64(math.Pi) / 180
+}
+
+// Converts radians to degrees
+func RadToDeg(angle float64) float64 {
+	return angle * 180 / float64(math.Pi)
+}

vendor/github.com/go-gl/mathgl/mgl64/doc.go

@@ -0,0 +1,24 @@
+// This file is generated from mgl32/doc.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+/*
+Package mgl[32|64] (an abbreviation of mathgl since the packages were split between 32 and 64-bit versions)
+is a pure Go math package specialized for 3D math, with inspiration from GLM. It provides statically-sized vectors and matrices with
+compile-time generated calculations for most basic math operations. It also provides several basic graphics utilities such as bezier curves and surfaces,
+generation of basic primitives like circles, easy creation of common matrices such as perspective or rotation, and common operations like converting
+to/from screen/OpenGL coordinates or Projecting/Unprojecting from an MVP matrix. Quaternions are also supported.
+
+The basic vectors and matrices are written with code generation, so looking directly at the source will probably be a bit confusing. I recommend looking at the Godoc
+instead, as all basic functions are documented.
+
+This package is written in Column Major Order to make it easier with OpenGL. This means for uniform blocks you can use the default ordering, and when you call
+pass-in functions you can leave the "transpose" argument as false.
+
+The package now contains variable sized vectors and matrices. Using these is discouraged. They exist for corner cases where you need "small" matrices that are still
+bigger than 4x4. An example may be a Jacobean used for inverse kinematics. Things like computer vision or general linear algebra are best left to packages
+more directly suited for that task -- OpenCV, BLAS, LAPACK, numpy, gonum (if you want to stay in Go), and so on.
+*/
+package mgl64

vendor/github.com/go-gl/mathgl/mgl64/matmn.go

@@ -0,0 +1,496 @@
+// This file is generated from mgl32/matmn.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mgl64
+
+import (
+	"math"
+)
+
+// An arbitrary mxn matrix backed by a slice of floats.
+//
+// This is emphatically not recommended for hardcore n-dimensional
+// linear algebra. For that purpose I recommend github.com/gonum/matrix or
+// well-tested C libraries such as BLAS or LAPACK.
+//
+// This is meant to complement future algorithms that may require matrices larger than
+// 4x4, but still relatively small (e.g. Jacobeans for inverse kinematics).
+//
+// It makes use of the same memory sync.Pool set that VecN does, with the same sizing rules.
+//
+// MatMN will always check if the receiver is nil on any method. Meaning MathMN(nil).Add(dst,m2)
+// should always work. Except for the Reshape function, the semantics of this is to "propogate" nils
+// forward, so if an invalid operation occurs in a long chain of matrix operations, the overall result will be nil.
+type MatMxN struct {
+	m, n int
+	dat  []float64
+}
+
+// Creates a matrix backed by a new slice of size m*n
+func NewMatrix(m, n int) (mat *MatMxN) {
+	if shouldPool {
+		return &MatMxN{m: m, n: n, dat: grabFromPool(m * n)}
+	} else {
+		return &MatMxN{m: m, n: n, dat: make([]float64, m*n)}
+	}
+}
+
+// Returns a matrix with data specified by the data in src
+//
+// For instance, to create a 3x3 MatMN from a Mat3
+//
+//    m1 := mgl32.Rotate3DX(3.14159)
+//    mat := mgl32.NewBackedMatrix(m1[:],3,3)
+//
+// will create an MN matrix matching the data in the original
+// rotation matrix. This matrix is NOT backed by the initial slice;
+// it's a copy of the data
+//
+// If m*n > cap(src), this function will panic.
+func NewMatrixFromData(src []float64, m, n int) *MatMxN {
+	var internal []float64
+	if shouldPool {
+		internal = grabFromPool(m * n)
+	} else {
+		internal = make([]float64, m*n)
+	}
+	copy(internal, src[:m*n])
+
+	return &MatMxN{m: m, n: n, dat: internal}
+}
+
+// Copies src into dst. This Reshapes dst
+// to the same size as src.
+//
+// If dst or src is nil, this is a no-op
+func CopyMatMN(dst, src *MatMxN) {
+	if dst == nil || src == nil {
+		return
+	}
+	dst.Reshape(src.m, src.n)
+	copy(dst.dat, src.dat)
+}
+
+// Stores the NxN identity matrix in dst, reallocating as necessary.
+func IdentN(dst *MatMxN, n int) *MatMxN {
+	dst = dst.Reshape(n, n)
+
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			if i == j {
+				dst.Set(i, j, 1)
+			} else {
+				dst.Set(i, j, 0)
+			}
+		}
+	}
+
+	return dst
+}
+
+// Creates an NxN diagonal matrix seeded by the diagonal vector
+// diag. Meaning: for all entries, where i==j, dst.At(i,j) = diag[i]. Otherwise
+// dst.At(i,j) = 0
+//
+// This reshapes dst to the correct size, returning/grabbing from the memory pool as necessary.
+func DiagN(dst *MatMxN, diag *VecN) *MatMxN {
+	dst = dst.Reshape(len(diag.vec), len(diag.vec))
+	n := len(diag.vec)
+
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			if i == j {
+				dst.Set(i, j, diag.vec[i])
+			} else {
+				dst.Set(i, j, 0)
+			}
+		}
+	}
+
+	return dst
+}
+
+// Reshapes the matrix to m by n and zeroes out all
+// elements.
+func (mat *MatMxN) Zero(m, n int) {
+	if mat == nil {
+		return
+	}
+
+	mat.Reshape(m, n)
+	for i := range mat.dat {
+		mat.dat[i] = 0
+	}
+}
+
+// Returns the underlying matrix slice to the memory pool
+func (mat *MatMxN) destroy() {
+	if mat == nil {
+		return
+	}
+
+	if shouldPool && mat.dat != nil {
+		returnToPool(mat.dat)
+	}
+	mat.m, mat.n = 0, 0
+	mat.dat = nil
+}
+
+// Reshapes the matrix to the desired dimensions.
+// If the overall size of the new matrix (m*n) is bigger
+// than the current size, the underlying slice will
+// be grown, sending the current slice to the memory pool
+// and grabbing a bigger one if necessary
+//
+// If the caller is a nil pointer, the return value will be a new
+// matrix, as if NewMatrix(m,n) had been called. Otherwise it's
+// simply the caller.
+func (mat *MatMxN) Reshape(m, n int) *MatMxN {
+	if mat == nil {
+		return NewMatrix(m, n)
+	}
+
+	if m*n <= cap(mat.dat) {
+		if mat.dat != nil {
+			mat.dat = mat.dat[:m*n]
+		} else {
+			mat.dat = []float64{}
+		}
+		mat.m, mat.n = m, n
+		return mat
+	}
+
+	if shouldPool && mat.dat != nil {
+		returnToPool(mat.dat)
+	}
+	(*mat) = (*NewMatrix(m, n))
+
+	return mat
+}
+
+// Infers an MxN matrix from a constant matrix from this package. For instance,
+// a Mat2x3 inferred with this function will work just like NewMatrixFromData(m[:],2,3)
+// where m is the Mat2x3. This uses a type switch.
+//
+// I personally recommend using NewMatrixFromData, because it avoids a potentially costly type switch.
+// However, this is also more robust and less error prone if you change the size of your matrix somewhere.
+//
+// If the value passed in is not recognized, it returns an InferMatrixError.
+func (mat *MatMxN) InferMatrix(m interface{}) (*MatMxN, error) {
+	switch raw := m.(type) {
+	case Mat2:
+		return NewMatrixFromData(raw[:], 2, 2), nil
+	case Mat2x3:
+		return NewMatrixFromData(raw[:], 2, 3), nil
+	case Mat2x4:
+		return NewMatrixFromData(raw[:], 2, 4), nil
+	case Mat3:
+		return NewMatrixFromData(raw[:], 3, 3), nil
+	case Mat3x2:
+		return NewMatrixFromData(raw[:], 3, 2), nil
+	case Mat3x4:
+		return NewMatrixFromData(raw[:], 3, 4), nil
+	case Mat4:
+		return NewMatrixFromData(raw[:], 4, 4), nil
+	case Mat4x2:
+		return NewMatrixFromData(raw[:], 4, 2), nil
+	case Mat4x3:
+		return NewMatrixFromData(raw[:], 4, 3), nil
+	default:
+		return nil, InferMatrixError{}
+	}
+}
+
+// Returns the trace of a square matrix (sum of all diagonal elements). If the matrix
+// is nil, or not square, the result will be NaN.
+func (mat *MatMxN) Trace() float64 {
+	if mat == nil || mat.m != mat.n {
+		return float64(math.NaN())
+	}
+
+	var out float64
+	for i := 0; i < mat.m; i++ {
+		out += mat.At(i, i)
+	}
+
+	return out
+}
+
+// Takes the transpose of mat and puts it in dst.
+//
+// If dst is not of the correct dimensions, it will be Reshaped,
+// if dst and mat are the same, a temporary matrix of the correct size will
+// be allocated; these resources will be released via the memory pool.
+//
+// This should be improved in the future.
+func (mat *MatMxN) Transpose(dst *MatMxN) (t *MatMxN) {
+	if mat == nil {
+		return nil
+	}
+
+	if dst == mat {
+		dst = NewMatrix(mat.n, mat.m)
+
+		// Copy data to correct matrix,
+		// delete temporary buffer,
+		// and set the return value to the
+		// correct one
+		defer func() {
+			copy(mat.dat, dst.dat)
+
+			mat.m, mat.n = mat.n, mat.m
+
+			dst.destroy()
+			t = mat
+		}()
+
+		return mat
+	} else {
+		dst = dst.Reshape(mat.n, mat.m)
+	}
+
+	for r := 0; r < mat.m; r++ {
+		for c := 0; c < mat.n; c++ {
+			dst.dat[r*dst.m+c] = mat.dat[c*mat.m+r]
+		}
+	}
+
+	return dst
+}
+
+// Returns the raw slice backing this matrix
+func (mat *MatMxN) Raw() []float64 {
+	if mat == nil {
+		return nil
+	}
+
+	return mat.dat
+}
+
+// Returns the number of rows in this matrix
+func (mat *MatMxN) NumRows() int {
+	return mat.m
+}
+
+// Returns the number of columns in this matrix
+func (mat *MatMxN) NumCols() int {
+	return mat.n
+}
+
+// Returns the number of rows and columns in this matrix
+// as a single operation
+func (mat *MatMxN) NumRowCols() (rows, cols int) {
+	return mat.m, mat.n
+}
+
+// Returns the element at the given row and column.
+// This is garbage in/garbage out and does no bounds
+// checking. If the computation happens to lead to an invalid
+// element, it will be returned; or it may panic.
+func (mat *MatMxN) At(row, col int) float64 {
+	return mat.dat[col*mat.m+row]
+}
+
+// Sets the element at the given row and column.
+// This is garbage in/garbage out and does no bounds
+// checking. If the computation happens to lead to an invalid
+// element, it will be set; or it may panic.
+func (mat *MatMxN) Set(row, col int, val float64) {
+	mat.dat[col*mat.m+row] = val
+}
+
+func (mat *MatMxN) Add(dst *MatMxN, addend *MatMxN) *MatMxN {
+	if mat == nil || addend == nil || mat.m != addend.m || mat.n != addend.n {
+		return nil
+	}
+
+	dst = dst.Reshape(mat.m, mat.n)
+
+	// No need to care about rows and columns
+	// since it's element-wise anyway
+	for i, el := range mat.dat {
+		dst.dat[i] = el + addend.dat[i]
+	}
+
+	return dst
+}
+
+func (mat *MatMxN) Sub(dst *MatMxN, subtrahend *MatMxN) *MatMxN {
+	if mat == nil || subtrahend == nil || mat.m != subtrahend.m || mat.n != subtrahend.n {
+		return nil
+	}
+
+	dst = dst.Reshape(mat.m, mat.n)
+
+	// No need to care about rows and columns
+	// since it's element-wise anyway
+	for i, el := range mat.dat {
+		dst.dat[i] = el - subtrahend.dat[i]
+	}
+
+	return dst
+}
+
+// Performs matrix multiplication on MxN matrix mat and NxO matrix mul, storing the result in dst.
+// This returns dst, or nil if the operation is not able to be performed.
+//
+// If mat == dst, or mul == dst a temporary matrix will be used.
+//
+// This uses the naive algorithm (though on smaller matrices,
+// this can actually be faster; about len(mat)+len(mul) < ~100)
+func (mat *MatMxN) MulMxN(dst *MatMxN, mul *MatMxN) *MatMxN {
+	if mat == nil || mul == nil || mat.n != mul.m {
+		return nil
+	}
+
+	if dst == mul {
+		mul = NewMatrix(mul.m, mul.n)
+		copy(mul.dat, dst.dat)
+
+		// If mat==dst==mul, we need to change
+		// mat too or we have a bug
+		if mat == dst {
+			mat = mul
+		}
+
+		defer mul.destroy()
+	} else if dst == mat {
+		mat = NewMatrix(mat.m, mat.n)
+		copy(mat.dat, dst.dat)
+
+		defer mat.destroy()
+	}
+
+	dst = dst.Reshape(mat.m, mul.n)
+	for r1 := 0; r1 < mat.m; r1++ {
+		for c2 := 0; c2 < mul.n; c2++ {
+
+			dst.dat[c2*mat.m+r1] = 0
+			for i := 0; i < mat.n; i++ {
+				dst.dat[c2*mat.m+r1] += mat.dat[i*mat.m+r1] * mul.dat[c2*mul.m+i]
+			}
+
+		}
+	}
+
+	return dst
+}
+
+// Performs a scalar multiplication between mat and some constant c,
+// storing the result in dst. Mat and dst can be equal. If dst is not the
+// correct size, a Reshape will occur.
+func (mat *MatMxN) Mul(dst *MatMxN, c float64) *MatMxN {
+	if mat == nil {
+		return nil
+	}
+
+	dst = dst.Reshape(mat.m, mat.n)
+
+	for i, el := range mat.dat {
+		dst.dat[i] = el * c
+	}
+
+	return dst
+}
+
+// Multiplies the matrix by a vector of size n. If mat or v is
+// nil, this returns nil. If the number of columns in mat does not match
+// the Size of v, this also returns nil.
+//
+// Dst will be resized if it's not big enough. If dst == v; a temporary
+// vector will be allocated and returned via the realloc callback when complete.
+func (mat *MatMxN) MulNx1(dst, v *VecN) *VecN {
+	if mat == nil || v == nil || mat.n != len(v.vec) {
+		return nil
+	}
+	if dst == v {
+		v = NewVecN(len(v.vec))
+		copy(v.vec, dst.vec)
+
+		defer v.destroy()
+	}
+
+	dst = dst.Resize(mat.m)
+
+	for r := 0; r < mat.m; r++ {
+		dst.vec[r] = 0
+
+		for c := 0; c < mat.n; c++ {
+			dst.vec[r] += mat.At(r, c) * v.vec[c]
+		}
+	}
+
+	return dst
+}
+
+func (mat *MatMxN) ApproxEqual(m2 *MatMxN) bool {
+	if mat == m2 {
+		return true
+	}
+	if mat.m != m2.m || mat.n != m2.n {
+		return false
+	}
+
+	for i, el := range mat.dat {
+		if !FloatEqual(el, m2.dat[i]) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (mat *MatMxN) ApproxEqualThreshold(m2 *MatMxN, epsilon float64) bool {
+	if mat == m2 {
+		return true
+	}
+	if mat.m != m2.m || mat.n != m2.n {
+		return false
+	}
+
+	for i, el := range mat.dat {
+		if !FloatEqualThreshold(el, m2.dat[i], epsilon) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (mat *MatMxN) ApproxEqualFunc(m2 *MatMxN, comp func(float64, float64) bool) bool {
+	if mat == m2 {
+		return true
+	}
+	if mat.m != m2.m || mat.n != m2.n {
+		return false
+	}
+
+	for i, el := range mat.dat {
+		if !comp(el, m2.dat[i]) {
+			return false
+		}
+	}
+
+	return true
+}
+
+type InferMatrixError struct{}
+
+func (me InferMatrixError) Error() string {
+	return "could not infer matrix. Make sure you're using a constant matrix such as Mat3 from within the same package (meaning: mgl32.MatMxN can't handle a mgl64.Mat2x3)."
+}
+
+type RectangularMatrixError struct{}
+
+func (mse RectangularMatrixError) Error() string {
+	return "the matrix was the wrong shape, needed a square matrix."
+}
+
+type NilMatrixError struct{}
+
+func (me NilMatrixError) Error() string {
+	return "the matrix is nil"
+}

vendor/github.com/go-gl/mathgl/mgl64/matstack/matstack.go

@@ -0,0 +1,61 @@
+// This file is generated from mgl32/matstack\matstack.go; DO NOT EDIT
+
+package matstack
+
+import (
+	"errors"
+
+	"github.com/go-gl/mathgl/mgl64"
+)
+
+// A MatStack is an OpenGL-style matrix stack,
+// usually used for things like scenegraphs. This allows you
+// to easily maintain matrix state per call level.
+type MatStack []mgl64.Mat4
+
+func NewMatStack() *MatStack {
+	return &MatStack{mgl64.Ident4()}
+}
+
+// Copies the top element and pushes it on the stack.
+func (ms *MatStack) Push() {
+	(*ms) = append(*ms, (*ms)[len(*ms)-1])
+}
+
+// Removes the first element of the matrix from the stack, if there is only one element left
+// there is an error.
+func (ms *MatStack) Pop() error {
+	if len(*ms) == 1 {
+		return errors.New("Cannot pop from mat stack, at minimum stack length of 1")
+	}
+	(*ms) = (*ms)[:len(*ms)-1]
+
+	return nil
+}
+
+// Right multiplies the current top of the matrix by the
+// argument.
+func (ms *MatStack) RightMul(m mgl64.Mat4) {
+	(*ms)[len(*ms)-1] = (*ms)[len(*ms)-1].Mul4(m)
+}
+
+// Left multiplies the current top of the matrix by the
+// argument.
+func (ms *MatStack) LeftMul(m mgl64.Mat4) {
+	(*ms)[len(*ms)-1] = m.Mul4((*ms)[len(*ms)-1])
+}
+
+// Returns the top element.
+func (ms *MatStack) Peek() mgl64.Mat4 {
+	return (*ms)[len(*ms)-1]
+}
+
+// Rewrites the top element of the stack with m
+func (ms *MatStack) Load(m mgl64.Mat4) {
+	(*ms)[len(*ms)-1] = m
+}
+
+// A shortcut for Load(mgl.Ident4())
+func (ms *MatStack) LoadIdent() {
+	(*ms)[len(*ms)-1] = mgl64.Ident4()
+}

vendor/github.com/go-gl/mathgl/mgl64/matstack/transformStack.go

@@ -0,0 +1,176 @@
+// This file is generated from mgl32/matstack\transformStack.go; DO NOT EDIT
+
+package matstack
+
+import (
+	"errors"
+	"fmt"
+
+	"github.com/go-gl/mathgl/mgl64"
+)
+
+// A transform stack is a linear fully-persistent data structure of matrix multiplications
+// Each push to a TransformStack multiplies the current top of the stack with thew new matrix
+// and appends it to the top. Each pop undoes the previous multiplication.
+//
+// This allows arbitrary unwinding of transformations, at the cost of a lot of memory. A notable feature
+// is the reseed and rebase, which allow invertible transformations to be rewritten as if a different transform
+// had been made in the middle.
+type TransformStack []mgl64.Mat4
+
+// Returns a matrix stack where the top element is the identity.
+func NewTransformStack() *TransformStack {
+	ms := make(TransformStack, 1)
+	ms[0] = mgl64.Ident4()
+
+	return &ms
+}
+
+// Multiplies the current top matrix by m, and pushes the result
+// on the stack.
+func (ms *TransformStack) Push(m mgl64.Mat4) {
+	prev := (*ms)[len(*ms)-1]
+	(*ms) = append(*ms, prev.Mul4(m))
+}
+
+// Pops the current matrix off the top of the stack and returns it.
+// If the matrix stack only has one element left, this will return an error.
+func (ms *TransformStack) Pop() (mgl64.Mat4, error) {
+	if len(*ms) == 1 {
+		return mgl64.Mat4{}, errors.New("attempt to pop last element of the stack; Matrix Stack must have at least one element")
+	}
+
+	retVal := (*ms)[len(*ms)-1]
+
+	(*ms) = (*ms)[:len(*ms)-1]
+
+	return retVal, nil
+}
+
+// Returns the value of the current top element of the stack, without
+// removing it.
+func (ms *TransformStack) Peek() mgl64.Mat4 {
+	return (*ms)[len(*ms)-1]
+}
+
+// Returns the size of the matrix stack. This value will never be less
+// than 1.
+func (ms *TransformStack) Len() int {
+	return len(*ms)
+}
+
+// This cuts down the matrix as if Pop had been called n times. If n would
+// bring the matrix down below 1 element, this does nothing and returns an error.
+func (ms *TransformStack) Unwind(n int) error {
+	if n > len(*ms)-1 {
+		return errors.New("Cannot unwind a matrix to below 1 value")
+	}
+
+	(*ms) = (*ms)[:len(*ms)-n]
+	return nil
+}
+
+// Copy will create a new "branch" of the current matrix stack,
+// the copy will contain all elements of the current stack in a new stack. Changes to
+// one will never affect the other.
+func (ms *TransformStack) Copy() *TransformStack {
+	v := append(TransformStack{}, (*ms)...)
+	return &v
+}
+
+// Reseed is tricky. It attempts to seed an arbitrary point in the matrix and replay all transformations
+// as if that point in the push had been the argument "change" instead of the original value.
+// The matrix stack does NOT keep track of arguments so this is done via consecutive inverses.
+// If the inverse of element i can be found, we can calculate the transformation that was given at point i+1.
+// This transformation can then be multiplied by the NEW matrix at point i to complete the "what if".
+// If no such inverse can be found at any given point along the rebase, it will be aborted, and the original
+// stack will NOT be visibly affected. The error returned will be of type NoInverseError.
+//
+// If n is out of bounds (n <= 0 || n >= len(*ms)), a generic error from the errors package will be returned.
+//
+// If you have the old transformations retained, it is recommended
+// that you use Unwind followed by Push(change) and then further calling Push for each transformation. Rebase is
+// imprecise by nature, and sometimes impossible. It's also expensive due to the inverse calculation at each point.
+func (ms *TransformStack) Reseed(n int, change mgl64.Mat4) error {
+	if n >= len(*ms) || n <= 0 {
+		return errors.New("Cannot rebase at the given point on the stack, it is out of bounds.")
+	}
+
+	return ms.reseed(n, change)
+}
+
+// Operates like reseed with no bounds checking; allows us to overwrite
+// the leading identity matrix with Rebase.
+func (ms *TransformStack) reseed(n int, change mgl64.Mat4) error {
+	backup := []mgl64.Mat4((*ms)[n:])
+	backup = append([]mgl64.Mat4{}, backup...) // copy into new slice
+
+	curr := (*ms)[n]
+	(*ms)[n] = (*ms)[n-1].Mul4(change)
+
+	for i := n + 1; i < len(*ms); i++ {
+		inv := curr.Inv()
+
+		blank := mgl64.Mat4{}
+		if inv == blank {
+			ms.undoRebase(n, backup)
+			return NoInverseError{Loc: i - 1, Mat: curr}
+		}
+
+		ghost := inv.Mul4((*ms)[i])
+
+		curr = (*ms)[i]
+		(*ms)[i] = (*ms)[i-1].Mul4(ghost)
+	}
+
+	return nil
+}
+
+func (ms *TransformStack) undoRebase(n int, prev []mgl64.Mat4) {
+	for i := n; i < len(*ms); i++ {
+		(*ms)[i] = prev[i-n]
+	}
+}
+
+// Rebase replays the current matrix stack as if the transformation that occurred at index "from"
+// in ms had instead started at the top of m.
+//
+// This returns a brand new stack containing all of m followed by all transformations
+// at from and after on ms as if they has been done on m instead.
+func Rebase(ms *TransformStack, from int, m *TransformStack) (*TransformStack, error) {
+	if from <= 0 || from >= len(*ms) {
+		return nil, errors.New("Cannot rebase, index out of range")
+	}
+
+	// Shift tmp so that the element immediately
+	// preceding our target is the "top" element of the list.
+	tmp := ms.Copy()
+	if from == 1 {
+		(*tmp) = append(*tmp, mgl64.Mat4{})
+	}
+	copy((*tmp)[1:], (*tmp)[from-1:])
+	if from-2 > 0 {
+		(*tmp) = (*tmp)[:len(*tmp)-(from-2)]
+	}
+
+	err := tmp.Reseed(1, m.Peek())
+	if err != nil {
+		return nil, err
+	}
+
+	(*tmp) = append(*m, (*tmp)[2:]...)
+
+	return tmp, nil
+}
+
+// A NoInverseError is returned on rebase when an inverse cannot be found along the chain,
+// due to a transformation projecting the matrix into a singularity. The values include the matrix
+// no inverse can be found for, and the location of that matrix.
+type NoInverseError struct {
+	Mat mgl64.Mat4
+	Loc int
+}
+
+func (nie NoInverseError) Error() string {
+	return fmt.Sprintf("cannot find inverse of matrix %v at location %d in matrix stack, aborting rebase/reseed", nie.Mat, nie.Loc)
+}

vendor/github.com/go-gl/mathgl/mgl64/mempool.go

@@ -0,0 +1,123 @@
+// This file is generated from mgl32/mempool.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mgl64
+
+import (
+	"sync"
+)
+
+var (
+	slicePools []*sync.Pool
+	listLock   sync.RWMutex
+)
+
+var shouldPool = true
+
+func DisableMemoryPooling() {
+	shouldPool = false
+}
+
+// Returns the given memory pool. If the pool doesn't exist, it will
+// create all pools up to element i. The number "i" corresponds to "p"
+// in most other comments. That is, it's Ceil(log_2(whatever)). So i=0
+// means you'll get the pool for slices of size 1, i=1 for size 2, i=2 for size 4,
+// and so on.
+//
+// This is concurrency safe and uses an RWMutex to protect the list expansion.
+func getPool(i int) *sync.Pool {
+	listLock.RLock()
+	if i >= len(slicePools) {
+
+		// Promote to a write lock because we now
+		// need to mutate the pool
+		listLock.RUnlock()
+		listLock.Lock()
+		defer listLock.Unlock()
+
+		for n := i - len(slicePools); n >= 0; n-- {
+			newFunc := genPoolNew(1 << uint(len(slicePools)))
+			slicePools = append(slicePools, &sync.Pool{New: newFunc})
+		}
+	} else {
+		defer listLock.RUnlock()
+	}
+
+	return slicePools[i]
+}
+
+func genPoolNew(i int) func() interface{} {
+	return func() interface{} {
+		return make([]float64, 0, i)
+	}
+}
+
+// Grabs a slice from the memory pool, such that its cap
+// is 2^p where p is Ceil(log_2(size)). It will be downsliced
+// such that the len is size.
+func grabFromPool(size int) []float64 {
+	pool, exact := binLog(size)
+
+	// Tried to grab something of size
+	// zero or less
+	if pool == -1 {
+		return nil
+	}
+
+	// If the log is not exact, we
+	// need to "overallocate" so we have
+	// log+1
+	if !exact {
+		pool++
+	}
+
+	slice := getPool(pool).Get().([]float64)
+	slice = slice[:size]
+	return slice
+}
+
+// Returns a slice to the appropriate pool. If the slice does not have a cap that's precisely
+// a power of 2, this will panic.
+func returnToPool(slice []float64) {
+	if cap(slice) == 0 {
+		return
+	}
+
+	pool, exact := binLog(cap(slice))
+
+	if !exact {
+		panic("attempt to pool slice with non-exact cap. If you're a user, please file an issue with github.com/go-gl/mathgl about this bug. This should never happen.")
+	}
+
+	getPool(pool).Put(slice)
+}
+
+// This returns the integer base 2 log of the value
+// and whether the log is exact or rounded down.
+//
+// This is only for positive integers.
+//
+// There are faster ways to do this, I'm open to suggestions. Most rely on knowing system endianness
+// which Go makes hard to do. I'm hesistant to use float conversions and the math package because of off-by-one errors.
+func binLog(val int) (int, bool) {
+	if val <= 0 {
+		return -1, false
+	}
+
+	exact := true
+	l := 0
+	for ; val > 1; val = val >> 1 {
+		// If the current lsb is 1 and the number
+		// is not equal to 1, this is not an exact
+		// log, but rather a rounding of it
+		if val&1 != 0 {
+			exact = false
+		}
+		l++
+	}
+
+	return l, exact
+}

vendor/github.com/go-gl/mathgl/mgl64/project.go

@@ -0,0 +1,100 @@
+// This file is generated from mgl32/project.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mgl64
+
+import (
+	"errors"
+	"math"
+)
+
+func Ortho(left, right, bottom, top, near, far float64) Mat4 {
+	rml, tmb, fmn := (right - left), (top - bottom), (far - near)
+
+	return Mat4{float64(2. / rml), 0, 0, 0, 0, float64(2. / tmb), 0, 0, 0, 0, float64(-2. / fmn), 0, float64(-(right + left) / rml), float64(-(top + bottom) / tmb), float64(-(far + near) / fmn), 1}
+}
+
+// Equivalent to Ortho with the near and far planes being -1 and 1, respectively
+func Ortho2D(left, right, bottom, top float64) Mat4 {
+	return Ortho(left, right, bottom, top, -1, 1)
+}
+
+func Perspective(fovy, aspect, near, far float64) Mat4 {
+	// fovy = (fovy * math.Pi) / 180.0 // convert from degrees to radians
+	nmf, f := near-far, float64(1./math.Tan(float64(fovy)/2.0))
+
+	return Mat4{float64(f / aspect), 0, 0, 0, 0, float64(f), 0, 0, 0, 0, float64((near + far) / nmf), -1, 0, 0, float64((2. * far * near) / nmf), 0}
+}
+
+func Frustum(left, right, bottom, top, near, far float64) Mat4 {
+	rml, tmb, fmn := (right - left), (top - bottom), (far - near)
+	A, B, C, D := (right+left)/rml, (top+bottom)/tmb, -(far+near)/fmn, -(2*far*near)/fmn
+
+	return Mat4{float64((2. * near) / rml), 0, 0, 0, 0, float64((2. * near) / tmb), 0, 0, float64(A), float64(B), float64(C), -1, 0, 0, float64(D), 0}
+}
+
+func LookAt(eyeX, eyeY, eyeZ, centerX, centerY, centerZ, upX, upY, upZ float64) Mat4 {
+	return LookAtV(Vec3{eyeX, eyeY, eyeZ}, Vec3{centerX, centerY, centerZ}, Vec3{upX, upY, upZ})
+}
+
+//  LookAtV generates a transform matrix from world space into the specific eye space
+func LookAtV(eye, center, up Vec3) Mat4 {
+	f := center.Sub(eye).Normalize()
+	s := f.Cross(up.Normalize()).Normalize()
+	u := s.Cross(f)
+
+	M := Mat4{
+		s[0], u[0], -f[0], 0,
+		s[1], u[1], -f[1], 0,
+		s[2], u[2], -f[2], 0,
+		0, 0, 0, 1,
+	}
+
+	return M.Mul4(Translate3D(float64(-eye[0]), float64(-eye[1]), float64(-eye[2])))
+}
+
+// Transform a set of coordinates from object space (in obj) to window coordinates (with depth)
+//
+// Window coordinates are continuous, not discrete (well, as continuous as an IEEE Floating Point can be), so you won't get exact pixel locations
+// without rounding or similar
+func Project(obj Vec3, modelview, projection Mat4, initialX, initialY, width, height int) (win Vec3) {
+	obj4 := obj.Vec4(1)
+
+	vpp := projection.Mul4(modelview).Mul4x1(obj4)
+	vpp = vpp.Mul(1 / vpp.W())
+	win[0] = float64(initialX) + (float64(width)*(vpp[0]+1))/2
+	win[1] = float64(initialY) + (float64(height)*(vpp[1]+1))/2
+	win[2] = (vpp[2] + 1) / 2
+
+	return win
+}
+
+// Transform a set of window coordinates to object space. If your MVP (projection.Mul(modelview) matrix is not invertible, this will return an error
+//
+// Note that the projection may not be perfect if you use strict pixel locations rather than the exact values given by Projectf.
+// (It's still unlikely to be perfect due to precision errors, but it will be closer)
+func UnProject(win Vec3, modelview, projection Mat4, initialX, initialY, width, height int) (obj Vec3, err error) {
+	inv := projection.Mul4(modelview).Inv()
+	var blank Mat4
+	if inv == blank {
+		return Vec3{}, errors.New("Could not find matrix inverse (projection times modelview is probably non-singular)")
+	}
+
+	obj4 := inv.Mul4x1(Vec4{
+		(2 * (win[0] - float64(initialX)) / float64(width)) - 1,
+		(2 * (win[1] - float64(initialY)) / float64(height)) - 1,
+		2*win[2] - 1,
+		1.0,
+	})
+	obj = obj4.Vec3()
+
+	//if obj4[3] > MinValue {}
+	obj[0] /= obj4[3]
+	obj[1] /= obj4[3]
+	obj[2] /= obj4[3]
+
+	return obj, nil
+}

vendor/github.com/go-gl/mathgl/mgl64/quat.go

@@ -0,0 +1,460 @@
+// This file is generated from mgl32/quat.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mgl64
+
+import (
+	"math"
+)
+
+// A rotation order is the order in which
+// rotations will be transformed for the purposes of AnglesToQuat
+type RotationOrder int
+
+const (
+	XYX RotationOrder = iota
+	XYZ
+	XZX
+	XZY
+	YXY
+	YXZ
+	YZY
+	YZX
+	ZYZ
+	ZYX
+	ZXZ
+	ZXY
+)
+
+// A Quaternion is an extension of the imaginary numbers; there's all sorts of
+// interesting theory behind it. In 3D graphics we mostly use it as a cheap way of
+// representing rotation since quaternions are cheaper to multiply by, and easier to
+// interpolate than matrices.
+//
+// A Quaternion has two parts: W, the so-called scalar component,
+// and "V", the vector component. The vector component is considered to
+// be the part in 3D space, while W (loosely interpreted) is its 4D coordinate.
+type Quat struct {
+	W float64
+	V Vec3
+}
+
+// The quaternion identity: W=1; V=(0,0,0).
+//
+// As with all identities, multiplying any quaternion by this will yield the same
+// quaternion you started with.
+func QuatIdent() Quat {
+	return Quat{1., Vec3{0, 0, 0}}
+}
+
+// Creates an angle from an axis and an angle relative to that axis.
+//
+// This is cheaper than HomogRotate3D.
+func QuatRotate(angle float64, axis Vec3) Quat {
+	// angle = (float32(math.Pi) * angle) / 180.0
+
+	c, s := float64(math.Cos(float64(angle/2))), float64(math.Sin(float64(angle/2)))
+
+	return Quat{c, axis.Mul(s)}
+}
+
+// A convenient alias for q.V[0]
+func (q Quat) X() float64 {
+	return q.V[0]
+}
+
+// A convenient alias for q.V[1]
+func (q Quat) Y() float64 {
+	return q.V[1]
+}
+
+// A convenient alias for q.V[2]
+func (q Quat) Z() float64 {
+	return q.V[2]
+}
+
+// Adds two quaternions. It's no more complicated than
+// adding their W and V components.
+func (q1 Quat) Add(q2 Quat) Quat {
+	return Quat{q1.W + q2.W, q1.V.Add(q2.V)}
+}
+
+// Subtracts two quaternions. It's no more complicated than
+// subtracting their W and V components.
+func (q1 Quat) Sub(q2 Quat) Quat {
+	return Quat{q1.W - q2.W, q1.V.Sub(q2.V)}
+}
+
+// Multiplies two quaternions. This can be seen as a rotation. Note that
+// Multiplication is NOT commutative, meaning q1.Mul(q2) does not necessarily
+// equal q2.Mul(q1).
+func (q1 Quat) Mul(q2 Quat) Quat {
+	return Quat{q1.W*q2.W - q1.V.Dot(q2.V), q1.V.Cross(q2.V).Add(q2.V.Mul(q1.W)).Add(q1.V.Mul(q2.W))}
+}
+
+// Scales every element of the quaternion by some constant factor.
+func (q1 Quat) Scale(c float64) Quat {
+	return Quat{q1.W * c, Vec3{q1.V[0] * c, q1.V[1] * c, q1.V[2] * c}}
+}
+
+// Returns the conjugate of a quaternion. Equivalent to
+// Quat{q1.W, q1.V.Mul(-1)}
+func (q1 Quat) Conjugate() Quat {
+	return Quat{q1.W, q1.V.Mul(-1)}
+}
+
+// Returns the Length of the quaternion, also known as its Norm. This is the same thing as
+// the Len of a Vec4
+func (q1 Quat) Len() float64 {
+	return float64(math.Sqrt(float64(q1.W*q1.W + q1.V[0]*q1.V[0] + q1.V[1]*q1.V[1] + q1.V[2]*q1.V[2])))
+}
+
+// Norm() is an alias for Len() since both are very common terms.
+func (q1 Quat) Norm() float64 {
+	return q1.Len()
+}
+
+// Normalizes the quaternion, returning its versor (unit quaternion).
+//
+// This is the same as normalizing it as a Vec4.
+func (q1 Quat) Normalize() Quat {
+	length := q1.Len()
+
+	if FloatEqual(1, length) {
+		return q1
+	}
+	if length == 0 {
+		return QuatIdent()
+	}
+	if length == InfPos {
+		length = MaxValue
+	}
+
+	return Quat{q1.W * 1 / length, q1.V.Mul(1 / length)}
+}
+
+// The inverse of a quaternion. The inverse is equivalent
+// to the conjugate divided by the square of the length.
+//
+// This method computes the square norm by directly adding the sum
+// of the squares of all terms instead of actually squaring q1.Len(),
+// both for performance and precision.
+func (q1 Quat) Inverse() Quat {
+	return q1.Conjugate().Scale(1 / q1.Dot(q1))
+}
+
+// Rotates a vector by the rotation this quaternion represents.
+// This will result in a 3D vector. Strictly speaking, this is
+// equivalent to q1.v.q* where the "."" is quaternion multiplication and v is interpreted
+// as a quaternion with W 0 and V v. In code:
+// q1.Mul(Quat{0,v}).Mul(q1.Conjugate()), and
+// then retrieving the imaginary (vector) part.
+//
+// In practice, we hand-compute this in the general case and simplify
+// to save a few operations.
+func (q1 Quat) Rotate(v Vec3) Vec3 {
+	cross := q1.V.Cross(v)
+	// v + 2q_w * (q_v x v) + 2q_v x (q_v x v)
+	return v.Add(cross.Mul(2 * q1.W)).Add(q1.V.Mul(2).Cross(cross))
+}
+
+// Returns the homogeneous 3D rotation matrix corresponding to the quaternion.
+func (q1 Quat) Mat4() Mat4 {
+	w, x, y, z := q1.W, q1.V[0], q1.V[1], q1.V[2]
+	return Mat4{
+		1 - 2*y*y - 2*z*z, 2*x*y + 2*w*z, 2*x*z - 2*w*y, 0,
+		2*x*y - 2*w*z, 1 - 2*x*x - 2*z*z, 2*y*z + 2*w*x, 0,
+		2*x*z + 2*w*y, 2*y*z - 2*w*x, 1 - 2*x*x - 2*y*y, 0,
+		0, 0, 0, 1,
+	}
+}
+
+// The dot product between two quaternions, equivalent to if this was a Vec4
+func (q1 Quat) Dot(q2 Quat) float64 {
+	return q1.W*q2.W + q1.V[0]*q2.V[0] + q1.V[1]*q2.V[1] + q1.V[2]*q2.V[2]
+}
+
+// Returns whether the quaternions are approximately equal, as if
+// FloatEqual was called on each matching element
+func (q1 Quat) ApproxEqual(q2 Quat) bool {
+	return FloatEqual(q1.W, q2.W) && q1.V.ApproxEqual(q2.V)
+}
+
+// Returns whether the quaternions are approximately equal with a given tolerence, as if
+// FloatEqualThreshold was called on each matching element with the given epsilon
+func (q1 Quat) ApproxEqualThreshold(q2 Quat, epsilon float64) bool {
+	return FloatEqualThreshold(q1.W, q2.W, epsilon) && q1.V.ApproxEqualThreshold(q2.V, epsilon)
+}
+
+// Returns whether the quaternions are approximately equal using the given comparison function, as if
+// the function had been called on each individual element
+func (q1 Quat) ApproxEqualFunc(q2 Quat, f func(float64, float64) bool) bool {
+	return f(q1.W, q2.W) && q1.V.ApproxFuncEqual(q2.V, f)
+}
+
+// Returns whether the quaternions represents the same orientation
+//
+// Different values can represent the same orientation (q == -q) because quaternions avoid singularities
+// and discontinuities involved with rotation in 3 dimensions by adding extra dimensions
+func (q1 Quat) OrientationEqual(q2 Quat) bool {
+	return q1.OrientationEqualThreshold(q2, Epsilon)
+}
+
+// Returns whether the quaternions represents the same orientation with a given tolerence
+func (q1 Quat) OrientationEqualThreshold(q2 Quat, epsilon float64) bool {
+	return Abs(q1.Normalize().Dot(q2.Normalize())) > 1-epsilon
+}
+
+// Slerp is *S*pherical *L*inear Int*erp*olation, a method of interpolating
+// between two quaternions. This always takes the straightest path on the sphere between
+// the two quaternions, and maintains constant velocity.
+//
+// However, it's expensive and QuatSlerp(q1,q2) is not the same as QuatSlerp(q2,q1)
+func QuatSlerp(q1, q2 Quat, amount float64) Quat {
+	q1, q2 = q1.Normalize(), q2.Normalize()
+	dot := q1.Dot(q2)
+
+	// If the inputs are too close for comfort, linearly interpolate and normalize the result.
+	if dot > 0.9995 {
+		return QuatNlerp(q1, q2, amount)
+	}
+
+	// This is here for precision errors, I'm perfectly aware that *technically* the dot is bound [-1,1], but since Acos will freak out if it's not (even if it's just a liiiiitle bit over due to normal error) we need to clamp it
+	dot = Clamp(dot, -1, 1)
+
+	theta := float64(math.Acos(float64(dot))) * amount
+	c, s := float64(math.Cos(float64(theta))), float64(math.Sin(float64(theta)))
+	rel := q2.Sub(q1.Scale(dot)).Normalize()
+
+	return q1.Scale(c).Add(rel.Scale(s))
+}
+
+// *L*inear Int*erp*olation between two Quaternions, cheap and simple.
+//
+// Not excessively useful, but uses can be found.
+func QuatLerp(q1, q2 Quat, amount float64) Quat {
+	return q1.Add(q2.Sub(q1).Scale(amount))
+}
+
+// *Normalized* *L*inear Int*erp*olation between two Quaternions. Cheaper than Slerp
+// and usually just as good. This is literally Lerp with Normalize() called on it.
+//
+// Unlike Slerp, constant velocity isn't maintained, but it's much faster and
+// Nlerp(q1,q2) and Nlerp(q2,q1) return the same path. You should probably
+// use this more often unless you're suffering from choppiness due to the
+// non-constant velocity problem.
+func QuatNlerp(q1, q2 Quat, amount float64) Quat {
+	return QuatLerp(q1, q2, amount).Normalize()
+}
+
+// Performs a rotation in the specified order. If the order is not
+// a valid RotationOrder, this function will panic
+//
+// The rotation "order" is more of an axis descriptor. For instance XZX would
+// tell the function to interpret angle1 as a rotation about the X axis, angle2 about
+// the Z axis, and angle3 about the X axis again.
+//
+// Based off the code for the Matlab function "angle2quat", though this implementation
+// only supports 3 single angles as opposed to multiple angles.
+func AnglesToQuat(angle1, angle2, angle3 float64, order RotationOrder) Quat {
+	var s [3]float64
+	var c [3]float64
+
+	s[0], c[0] = math.Sincos(float64(angle1 / 2))
+	s[1], c[1] = math.Sincos(float64(angle2 / 2))
+	s[2], c[2] = math.Sincos(float64(angle3 / 2))
+
+	ret := Quat{}
+	switch order {
+	case ZYX:
+		ret.W = float64(c[0]*c[1]*c[2] + s[0]*s[1]*s[2])
+		ret.V = Vec3{float64(c[0]*c[1]*s[2] - s[0]*s[1]*c[2]),
+			float64(c[0]*s[1]*c[2] + s[0]*c[1]*s[2]),
+			float64(s[0]*c[1]*c[2] - c[0]*s[1]*s[2]),
+		}
+	case ZYZ:
+		ret.W = float64(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float64(c[0]*s[1]*s[2] - s[0]*s[1]*c[2]),
+			float64(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+			float64(s[0]*c[1]*c[2] + c[0]*c[1]*s[2]),
+		}
+	case ZXY:
+		ret.W = float64(c[0]*c[1]*c[2] - s[0]*s[1]*s[2])
+		ret.V = Vec3{float64(c[0]*s[1]*c[2] - s[0]*c[1]*s[2]),
+			float64(c[0]*c[1]*s[2] + s[0]*s[1]*c[2]),
+			float64(c[0]*s[1]*s[2] + s[0]*c[1]*c[2]),
+		}
+	case ZXZ:
+		ret.W = float64(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float64(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+			float64(s[0]*s[1]*c[2] - c[0]*s[1]*s[2]),
+			float64(c[0]*c[1]*s[2] + s[0]*c[1]*c[2]),
+		}
+	case YXZ:
+		ret.W = float64(c[0]*c[1]*c[2] + s[0]*s[1]*s[2])
+		ret.V = Vec3{float64(c[0]*s[1]*c[2] + s[0]*c[1]*s[2]),
+			float64(s[0]*c[1]*c[2] - c[0]*s[1]*s[2]),
+			float64(c[0]*c[1]*s[2] - s[0]*s[1]*c[2]),
+		}
+	case YXY:
+		ret.W = float64(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float64(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+			float64(s[0]*c[1]*c[2] + c[0]*c[1]*s[2]),
+			float64(c[0]*s[1]*s[2] - s[0]*s[1]*c[2]),
+		}
+	case YZX:
+		ret.W = float64(c[0]*c[1]*c[2] - s[0]*s[1]*s[2])
+		ret.V = Vec3{float64(c[0]*c[1]*s[2] + s[0]*s[1]*c[2]),
+			float64(c[0]*s[1]*s[2] + s[0]*c[1]*c[2]),
+			float64(c[0]*s[1]*c[2] - s[0]*c[1]*s[2]),
+		}
+	case YZY:
+		ret.W = float64(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float64(s[0]*s[1]*c[2] - c[0]*s[1]*s[2]),
+			float64(c[0]*c[1]*s[2] + s[0]*c[1]*c[2]),
+			float64(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+		}
+	case XYZ:
+		ret.W = float64(c[0]*c[1]*c[2] - s[0]*s[1]*s[2])
+		ret.V = Vec3{float64(c[0]*s[1]*s[2] + s[0]*c[1]*c[2]),
+			float64(c[0]*s[1]*c[2] - s[0]*c[1]*s[2]),
+			float64(c[0]*c[1]*s[2] + s[0]*s[1]*c[2]),
+		}
+	case XYX:
+		ret.W = float64(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float64(c[0]*c[1]*s[2] + s[0]*c[1]*c[2]),
+			float64(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+			float64(s[0]*s[1]*c[2] - c[0]*s[1]*s[2]),
+		}
+	case XZY:
+		ret.W = float64(c[0]*c[1]*c[2] + s[0]*s[1]*s[2])
+		ret.V = Vec3{float64(s[0]*c[1]*c[2] - c[0]*s[1]*s[2]),
+			float64(c[0]*c[1]*s[2] - s[0]*s[1]*c[2]),
+			float64(c[0]*s[1]*c[2] + s[0]*c[1]*s[2]),
+		}
+	case XZX:
+		ret.W = float64(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float64(c[0]*c[1]*s[2] + s[0]*c[1]*c[2]),
+			float64(c[0]*s[1]*s[2] - s[0]*s[1]*c[2]),
+			float64(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+		}
+	default:
+		panic("Unsupported rotation order")
+	}
+	return ret
+}
+
+// Mat4ToQuat converts a pure rotation matrix into a quaternion
+func Mat4ToQuat(m Mat4) Quat {
+	// http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm
+
+	if tr := m[0] + m[5] + m[10]; tr > 0 {
+		s := float64(0.5 / math.Sqrt(float64(tr+1.0)))
+		return Quat{
+			0.25 / s,
+			Vec3{
+				(m[6] - m[9]) * s,
+				(m[8] - m[2]) * s,
+				(m[1] - m[4]) * s,
+			},
+		}
+	}
+
+	if (m[0] > m[5]) && (m[0] > m[10]) {
+		s := float64(2.0 * math.Sqrt(float64(1.0+m[0]-m[5]-m[10])))
+		return Quat{
+			(m[6] - m[9]) / s,
+			Vec3{
+				0.25 * s,
+				(m[4] + m[1]) / s,
+				(m[8] + m[2]) / s,
+			},
+		}
+	}
+
+	if m[5] > m[10] {
+		s := float64(2.0 * math.Sqrt(float64(1.0+m[5]-m[0]-m[10])))
+		return Quat{
+			(m[8] - m[2]) / s,
+			Vec3{
+				(m[4] + m[1]) / s,
+				0.25 * s,
+				(m[9] + m[6]) / s,
+			},
+		}
+
+	}
+
+	s := float64(2.0 * math.Sqrt(float64(1.0+m[10]-m[0]-m[5])))
+	return Quat{
+		(m[1] - m[4]) / s,
+		Vec3{
+			(m[8] + m[2]) / s,
+			(m[9] + m[6]) / s,
+			0.25 * s,
+		},
+	}
+}
+
+// QuatLookAtV creates a rotation from an eye vector to a center vector
+//
+// It assumes the front of the rotated object at Z- and up at Y+
+func QuatLookAtV(eye, center, up Vec3) Quat {
+	// http://www.opengl-tutorial.org/intermediate-tutorials/tutorial-17-quaternions/#I_need_an_equivalent_of_gluLookAt__How_do_I_orient_an_object_towards_a_point__
+	// https://bitbucket.org/sinbad/ogre/src/d2ef494c4a2f5d6e2f0f17d3bfb9fd936d5423bb/OgreMain/src/OgreCamera.cpp?at=default#cl-161
+
+	direction := center.Sub(eye).Normalize()
+
+	// Find the rotation between the front of the object (that we assume towards Z-,
+	// but this depends on your model) and the desired direction
+	rotDir := QuatBetweenVectors(Vec3{0, 0, -1}, direction)
+
+	// Recompute up so that it's perpendicular to the direction
+	// You can skip that part if you really want to force up
+	//right := direction.Cross(up)
+	//up = right.Cross(direction)
+
+	// Because of the 1rst rotation, the up is probably completely screwed up.
+	// Find the rotation between the "up" of the rotated object, and the desired up
+	upCur := rotDir.Rotate(Vec3{0, 1, 0})
+	rotUp := QuatBetweenVectors(upCur, up)
+
+	rotTarget := rotUp.Mul(rotDir) // remember, in reverse order.
+	return rotTarget.Inverse()     // camera rotation should be inversed!
+}
+
+// QuatBetweenVectors calculates the rotation between two vectors
+func QuatBetweenVectors(start, dest Vec3) Quat {
+	// http://www.opengl-tutorial.org/intermediate-tutorials/tutorial-17-quaternions/#I_need_an_equivalent_of_gluLookAt__How_do_I_orient_an_object_towards_a_point__
+	// https://github.com/g-truc/glm/blob/0.9.5/glm/gtx/quaternion.inl#L225
+	// https://bitbucket.org/sinbad/ogre/src/d2ef494c4a2f5d6e2f0f17d3bfb9fd936d5423bb/OgreMain/include/OgreVector3.h?at=default#cl-654
+
+	start = start.Normalize()
+	dest = dest.Normalize()
+	epsilon := float64(0.001)
+
+	cosTheta := start.Dot(dest)
+	if cosTheta < -1.0+epsilon {
+		// special case when vectors in opposite directions:
+		// there is no "ideal" rotation axis
+		// So guess one; any will do as long as it's perpendicular to start
+		axis := Vec3{1, 0, 0}.Cross(start)
+		if axis.Dot(axis) < epsilon {
+			// bad luck, they were parallel, try again!
+			axis = Vec3{0, 1, 0}.Cross(start)
+		}
+
+		return QuatRotate(math.Pi, axis.Normalize())
+	}
+
+	axis := start.Cross(dest)
+	s := float64(math.Sqrt(float64(1.0+cosTheta) * 2.0))
+
+	return Quat{
+		s * 0.5,
+		axis.Mul(1.0 / s),
+	}
+}

vendor/github.com/go-gl/mathgl/mgl64/shapes.go

@@ -0,0 +1,308 @@
+// This file is generated from mgl32/shapes.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mgl64
+
+import (
+	"fmt"
+	"math"
+)
+
+// Generates a circle centered at (0,0) with a given radius.
+// The radii are assumed to be in GL's coordinate sizing.
+//
+// Technically this draws an ellipse with two axes that match with the X and Y axes, the reason it has a radiusX and radiusY is because GL's coordinate system
+// is proportional to screen width and screen height. So if you have a non-square viewport, a single radius will appear
+// to "squash" the circle in one direction (usually the Y direction), so the X and Y radius allow for a circle to be made.
+// A good way to get the correct radii is with mathgl.ScreenToGLCoords(radius, radius, screenWidth, screenHeight) which will get you the correct
+// proportional GL coords.
+//
+// The numSlices argument specifies how many triangles you want your circle divided into, setting this
+// number to too low a value may cause problem (and too high will cause it to take a lot of memory and time to compute
+// without much gain in resolution).
+//
+// This uses discrete triangles, not a triangle fan
+func Circle(radiusX, radiusY float64, numSlices int) []Vec2 {
+	twoPi := float64(2.0 * math.Pi)
+
+	circlePoints := make([]Vec2, 0, numSlices*3)
+	center := Vec2{0.0, 0.0}
+	previous := Vec2{radiusX, 0.0}
+
+	for theta := twoPi / float64(numSlices); !FloatEqual(theta, twoPi); theta = Clamp(theta+twoPi/float64(numSlices), 0.0, twoPi) {
+		sin, cos := math.Sincos(float64(theta))
+		curr := Vec2{float64(cos) * radiusX, float64(sin) * radiusY}
+
+		circlePoints = append(circlePoints, center, previous, curr)
+		previous = curr
+	}
+
+	// Now add the final point at theta=2pi
+	circlePoints = append(circlePoints, center, previous, Vec2{radiusX, 0.0})
+	return circlePoints
+}
+
+// Generates a 2-triangle rectangle for use with GL_TRIANGLES. The width and height should use GL's proportions (that is, where a width of 1.0
+// is equivalent to half of the width of the render target); however, the y-coordinates grow downwards, not upwards. That is, it
+// assumes you want the origin of the rectangle with the top-left corner at (0.0,0.0).
+//
+// Keep in mind that GL's coordinate system is proportional, so width=height will not result in a square unless your viewport is square.
+// If you want to maintain proportionality regardless of screen size, use the results of w,h := ScreenToGLCoordsf(absoluteWidth, absoluteHeight, screenWidth, screenHeight);
+// w,h=w+1,h-1 in the call to this function. (The w+1,h-1 step maps the coordinates to start at 0.0 rather than -1.0)
+func Rect(width, height float64) []Vec2 {
+	return []Vec2{
+		{0.0, 0.0},
+		{0.0, -height},
+		{width, -height},
+
+		{0.0, 0.0},
+		{width, -height},
+		{width, 0.0},
+	}
+}
+
+func QuadraticBezierCurve2D(t float64, cPoint1, cPoint2, cPoint3 Vec2) Vec2 {
+	if t < 0.0 || t > 1.0 {
+		panic("Can't interpolate on bezier curve with t out of range [0.0,1.0]")
+	}
+
+	return cPoint1.Mul((1.0 - t) * (1.0 - t)).Add(cPoint2.Mul(2 * (1 - t) * t)).Add(cPoint3.Mul(t * t))
+}
+
+func QuadraticBezierCurve3D(t float64, cPoint1, cPoint2, cPoint3 Vec3) Vec3 {
+	if t < 0.0 || t > 1.0 {
+		panic("Can't interpolate on bezier curve with t out of range [0.0,1.0]")
+	}
+
+	return cPoint1.Mul((1.0 - t) * (1.0 - t)).Add(cPoint2.Mul(2 * (1 - t) * t)).Add(cPoint3.Mul(t * t))
+}
+
+func CubicBezierCurve2D(t float64, cPoint1, cPoint2, cPoint3, cPoint4 Vec2) Vec2 {
+	if t < 0.0 || t > 1.0 {
+		panic("Can't interpolate on bezier curve with t out of range [0.0,1.0]")
+	}
+
+	return cPoint1.Mul((1 - t) * (1 - t) * (1 - t)).Add(cPoint2.Mul(3 * (1 - t) * (1 - t) * t)).Add(cPoint3.Mul(3 * (1 - t) * t * t)).Add(cPoint4.Mul(t * t * t))
+}
+
+func CubicBezierCurve3D(t float64, cPoint1, cPoint2, cPoint3, cPoint4 Vec3) Vec3 {
+	if t < 0.0 || t > 1.0 {
+		panic("Can't interpolate on bezier curve with t out of range [0.0,1.0]")
+	}
+
+	return cPoint1.Mul((1 - t) * (1 - t) * (1 - t)).Add(cPoint2.Mul(3 * (1 - t) * (1 - t) * t)).Add(cPoint3.Mul(3 * (1 - t) * t * t)).Add(cPoint4.Mul(t * t * t))
+}
+
+// Returns the point at point t along an n-control point Bezier curve
+//
+// t must be in the range 0.0 and 1.0 or this function will panic. Consider [0.0,1.0] to be similar to a percentage,
+// 0.0 is first control point, and the point at 1.0 is the last control point. Any point in between is how far along the path you are between 0 and 1.
+//
+// This function is not sensitive to the coordinate system of the control points. It will correctly interpolate regardless of whether they're in screen coords,
+// gl coords, or something else entirely
+func BezierCurve2D(t float64, cPoints []Vec2) Vec2 {
+	if t < 0.0 || t > 1.0 {
+		panic("Input to bezier has t not in range [0,1]. If you think this is a precision error, use mathgl.Clamp[f|d] before calling this function")
+	}
+
+	n := len(cPoints) - 1
+	point := cPoints[0].Mul(float64(math.Pow(float64(1.0-t), float64(n))))
+
+	for i := 1; i <= n; i++ {
+		point = point.Add(cPoints[i].Mul(float64(float64(choose(n, i)) * math.Pow(float64(1-t), float64(n-i)) * math.Pow(float64(t), float64(i))))) // P += P_i * nCi * (1-t)^(n-i) * t^i
+	}
+
+	return point
+}
+
+// Same as the 2D version, except the line is in 3D space
+func BezierCurve3D(t float64, cPoints []Vec3) Vec3 {
+	if t < 0.0 || t > 1.0 {
+		panic("Input to bezier has t not in range [0,1]. If you think this is a precision error, use mathgl.Clamp[f|d] before calling this function")
+	}
+
+	n := len(cPoints) - 1
+	point := cPoints[0].Mul(float64(math.Pow(float64(1.0-t), float64(n))))
+
+	for i := 1; i <= n; i++ {
+		point = point.Add(cPoints[i].Mul(float64(float64(choose(n, i)) * math.Pow(float64(1-t), float64(n-i)) * math.Pow(float64(t), float64(i))))) // P += P_i * nCi * (1-t)^(n-i) * t^i
+	}
+
+	return point
+}
+
+// Generates a bezier curve with controlPoints cPoints. The numPoints argument
+// determines how many "samples" it makes along the line. For instance, a
+// call to this with numPoints 2 will have exactly two points: the start and end points
+// For any points above 2 it will divide it into numPoints-1 chunks (which means it will generate numPoints-2 vertices other than the beginning and end).
+// So for 3 points it will divide it in half, 4 points into thirds, and so on.
+//
+// This is likely to get rather expensive for anything over perhaps a cubic curve.
+func MakeBezierCurve2D(numPoints int, cPoints []Vec2) (line []Vec2) {
+	line = make([]Vec2, numPoints)
+	if numPoints == 0 {
+		return
+	} else if numPoints == 1 {
+		line[0] = cPoints[0]
+		return
+	} else if numPoints == 2 {
+		line[0] = cPoints[0]
+		line[1] = cPoints[len(cPoints)-1]
+		return
+	}
+
+	line[0] = cPoints[0]
+	for i := 1; i < numPoints-1; i++ {
+		line[i] = BezierCurve2D(Clamp(float64(i)/float64(numPoints-1), 0.0, 1.0), cPoints)
+	}
+	line[numPoints-1] = cPoints[len(cPoints)-1]
+
+	return
+}
+
+// Same as the 2D version, except with the line in 3D space
+func MakeBezierCurve3D(numPoints int, cPoints []Vec3) (line []Vec3) {
+	line = make([]Vec3, numPoints)
+	if numPoints == 0 {
+		return
+	} else if numPoints == 1 {
+		line[0] = cPoints[0]
+		return
+	} else if numPoints == 2 {
+		line[0] = cPoints[0]
+		line[1] = cPoints[len(cPoints)-1]
+		return
+	}
+
+	line[0] = cPoints[0]
+	for i := 1; i < numPoints-1; i++ {
+		line[i] = BezierCurve3D(Clamp(float64(i)/float64(numPoints-1), 0.0, 1.0), cPoints)
+	}
+	line[numPoints-1] = cPoints[len(cPoints)-1]
+
+	return
+}
+
+// Creates a 2-dimensional Bezier surface of arbitrary degree in 3D Space
+// Like the curve functions, if u or v are not in the range [0.0,1.0] the function will panic, use Clamp[f|d]
+// to ensure it is correct.
+//
+// The control point matrix must not be jagged, or this will end up panicking from an index out of bounds exception
+func BezierSurface(u, v float64, cPoints [][]Vec3) Vec3 {
+	if u < 0.0 || u > 1.0 || v < 1.0 || v > 1.0 {
+		panic("u or v not in range [0.0,1.0] in BezierSurface")
+	}
+
+	n := len(cPoints) - 1
+	m := len(cPoints[0]) - 1
+
+	point := cPoints[0][0].Mul(float64(math.Pow(float64(1.0-u), float64(n)) * math.Pow(float64(1.0-v), float64(m))))
+
+	for i := 0; i <= n; i++ {
+		for j := 0; j <= m; j++ {
+			if i == 0 && j == 0 {
+				continue
+			}
+
+			point = point.Add(cPoints[i][j].Mul(float64(float64(choose(n, i)) * math.Pow(float64(u), float64(i)) * math.Pow(float64(1.0-u), float64(n-i)) * float64(choose(m, j)) * math.Pow(float64(v), float64(j)) * math.Pow(float64(1.0-v), float64(m-j)))))
+		}
+	}
+
+	return point
+}
+
+// Does interpolation over a spline of several bezier curves. Each bezier curve must have a finite range,
+// though the spline may be disjoint. The bezier curves are not required to be in any particular order.
+//
+// If t is out of the range of all given curves, this function will panic
+func BezierSplineInterpolate2D(t float64, ranges [][2]float64, cPoints [][]Vec2) Vec2 {
+	if len(ranges) != len(cPoints) {
+		panic("Each bezier curve needs a range")
+	}
+
+	for i, curveRange := range ranges {
+		if t >= curveRange[0] && t <= curveRange[1] {
+			return BezierCurve2D((t-curveRange[0])/(curveRange[1]-curveRange[0]), cPoints[i])
+		}
+	}
+
+	panic("t is out of the range of all bezier curves in this spline")
+}
+
+// Does interpolation over a spline of several bezier curves. Each bezier curve must have a finite range,
+// though the spline may be disjoint. The bezier curves are not required to be in any particular order.
+//
+// If t is out of the range of all given curves, this function will panic
+func BezierSplineInterpolate3D(t float64, ranges [][2]float64, cPoints [][]Vec3) Vec3 {
+	if len(ranges) != len(cPoints) {
+		panic("Each bezier curve needs a range")
+	}
+
+	for i, curveRange := range ranges {
+		if t >= curveRange[0] && t <= curveRange[1] {
+			return BezierCurve3D((t-curveRange[0])/(curveRange[1]-curveRange[0]), cPoints[i])
+		}
+	}
+
+	panic("t is out of the range of all bezier curves in this spline")
+}
+
+// Reticulates ALL the Splines
+//
+// For the overly serious: the function is just for fun. It does nothing except prints a Maxis reference. Technically you could "reticulate splines"
+// by joining a bunch of splines together, but that ruins the joke.
+func ReticulateSplines(ranges [][][2]float64, cPoints [][][]Vec2, withLlamas bool) {
+	if !withLlamas {
+		fmt.Println("You can't reticulate splines without llamas, silly.")
+	} else {
+		fmt.Println("Actually, you can't even reticulate splines WITH llamas")
+	}
+}
+
+// Transform from pixel coordinates to GL coordinates.
+//
+// This assumes that your pixel coordinate system considers its origin to be in the top left corner (GL's is in the bottom left).
+// The coordinates x and y may be out of the range [0,screenWidth-1] and [0,screeneHeight-1].
+//
+// GL's coordinate system maps [screenWidth-1,0] to [1.0,1.0] and [0,screenHeight-1] to [-1.0,-1.0]. If x and y are out of the range, they'll still
+// be mapped correctly, just off the screen. (e.g. if y = 2*(screenHeight-1) you'll get -3.0 for yOut)
+//
+// This is similar to Unproject, except for 2D cases and much simpler (especially since an inverse may always be found)
+func ScreenToGLCoords(x, y int, screenWidth, screenHeight int) (xOut, yOut float64) {
+	xOut = 2.0*float64(x)/float64(screenWidth-1) - 1.0
+	yOut = -2.0*float64(y)/float64(screenHeight-1) + 1.0
+
+	return
+}
+
+// Transform from GL's proportional system to pixel coordinates.
+//
+// Assumes the pixel coordinate system has its origin in the top left corner. (GL's is in the bottom left)
+//
+// GL's coordinate system maps [screenWidth-1,0] to [1.0,1.0] and [0,screenHeight-1] to [-1.0,-1.0]. If x and y are out of the range, they'll still
+// be mapped correctly, just off the screen. (e.g. if y=-3.0, you'll get 2*(screenHeight-1) for yOut)
+//
+// This is similar to Project, except for 2D cases and much simpler
+func GLToScreenCoords(x, y float64, screenWidth, screenHeight int) (xOut, yOut int) {
+	xOut = int((x + 1.0) * float64(screenWidth-1) / 2.0)
+	yOut = int((1.0 - y) * float64(screenHeight-1) / 2.0)
+
+	return
+}
+
+func choose(n, k int) (result int) {
+	if k == 0 {
+		return 1
+	} else if n == 0 {
+		return 0
+	}
+	result = (n - (k - 1))
+	for i := 2; i <= k; i++ {
+		result *= (n - (k - i)) / i
+	}
+
+	return result
+}

vendor/github.com/go-gl/mathgl/mgl64/transform.go

@@ -0,0 +1,225 @@
+// This file is generated from mgl32/transform.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mgl64
+
+import "math"
+
+// Rotate2D returns a rotation Matrix about a angle in 2-D space. Specifically about the origin.
+// It is a 2x2 matrix, if you need a 3x3 for Homogeneous math (e.g. composition with a Translation matrix)
+// see HomogRotate2D
+func Rotate2D(angle float64) Mat2 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float64(math.Sin(float64(angle))), float64(math.Cos(float64(angle)))
+	return Mat2{cos, sin, -sin, cos}
+}
+
+// Rotate3DX returns a 3x3 (non-homogeneous) Matrix that rotates by angle about the X-axis
+//
+// Where c is cos(angle) and s is sin(angle)
+//    [1 0 0]
+//    [0 c -s]
+//    [0 s c ]
+func Rotate3DX(angle float64) Mat3 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float64(math.Sin(float64(angle))), float64(math.Cos(float64(angle)))
+
+	return Mat3{1, 0, 0, 0, cos, sin, 0, -sin, cos}
+}
+
+// Rotate3DY returns a 3x3 (non-homogeneous) Matrix that rotates by angle about the Y-axis
+//
+// Where c is cos(angle) and s is sin(angle)
+//    [c 0 s]
+//    [0 1 0]
+//    [s 0 c ]
+func Rotate3DY(angle float64) Mat3 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float64(math.Sin(float64(angle))), float64(math.Cos(float64(angle)))
+
+	return Mat3{cos, 0, -sin, 0, 1, 0, sin, 0, cos}
+}
+
+// Rotate3DZ returns a 3x3 (non-homogeneous) Matrix that rotates by angle about the Z-axis
+//
+// Where c is cos(angle) and s is sin(angle)
+//    [c -s 0]
+//    [s c 0]
+//    [0 0 1 ]
+func Rotate3DZ(angle float64) Mat3 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float64(math.Sin(float64(angle))), float64(math.Cos(float64(angle)))
+
+	return Mat3{cos, sin, 0, -sin, cos, 0, 0, 0, 1}
+}
+
+// Translate2D returns a homogeneous (3x3 for 2D-space) Translation matrix that moves a point by Tx units in the x-direction and Ty units in the y-direction
+//
+//    [[1, 0, Tx]]
+//    [[0, 1, Ty]]
+//    [[0, 0, 1 ]]
+func Translate2D(Tx, Ty float64) Mat3 {
+	return Mat3{1, 0, 0, 0, 1, 0, float64(Tx), float64(Ty), 1}
+}
+
+// Translate3D returns a homogeneous (4x4 for 3D-space) Translation matrix that moves a point by Tx units in the x-direction, Ty units in the y-direction,
+// and Tz units in the z-direction
+//
+//    [[1, 0, 0, Tx]]
+//    [[0, 1, 0, Ty]]
+//    [[0, 0, 1, Tz]]
+//    [[0, 0, 0, 1 ]]
+func Translate3D(Tx, Ty, Tz float64) Mat4 {
+	return Mat4{1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, float64(Tx), float64(Ty), float64(Tz), 1}
+}
+
+// Same as Rotate2D, except homogeneous (3x3 with the extra row/col being all zeroes with a one in the bottom right)
+func HomogRotate2D(angle float64) Mat3 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float64(math.Sin(float64(angle))), float64(math.Cos(float64(angle)))
+	return Mat3{cos, sin, 0, -sin, cos, 0, 0, 0, 1}
+}
+
+// Same as Rotate3DX, except homogeneous (4x4 with the extra row/col being all zeroes with a one in the bottom right)
+func HomogRotate3DX(angle float64) Mat4 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float64(math.Sin(float64(angle))), float64(math.Cos(float64(angle)))
+
+	return Mat4{1, 0, 0, 0, 0, cos, sin, 0, 0, -sin, cos, 0, 0, 0, 0, 1}
+}
+
+// Same as Rotate3DY, except homogeneous (4x4 with the extra row/col being all zeroes with a one in the bottom right)
+func HomogRotate3DY(angle float64) Mat4 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float64(math.Sin(float64(angle))), float64(math.Cos(float64(angle)))
+	return Mat4{cos, 0, -sin, 0, 0, 1, 0, 0, sin, 0, cos, 0, 0, 0, 0, 1}
+}
+
+// Same as Rotate3DZ, except homogeneous (4x4 with the extra row/col being all zeroes with a one in the bottom right)
+func HomogRotate3DZ(angle float64) Mat4 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float64(math.Sin(float64(angle))), float64(math.Cos(float64(angle)))
+	return Mat4{cos, sin, 0, 0, -sin, cos, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}
+}
+
+// Scale3D creates a homogeneous 3D scaling matrix
+// [[ scaleX, 0     , 0     , 0 ]]
+// [[ 0     , scaleY, 0     , 0 ]]
+// [[ 0     , 0     , scaleZ, 0 ]]
+// [[ 0     , 0     , 0     , 1 ]]
+func Scale3D(scaleX, scaleY, scaleZ float64) Mat4 {
+
+	return Mat4{float64(scaleX), 0, 0, 0, 0, float64(scaleY), 0, 0, 0, 0, float64(scaleZ), 0, 0, 0, 0, 1}
+}
+
+// Scale2D creates a homogeneous 2D scaling matrix
+// [[ scaleX, 0     , 0 ]]
+// [[ 0     , scaleY, 0 ]]
+// [[ 0     , 0     , 1 ]]
+func Scale2D(scaleX, scaleY float64) Mat3 {
+	return Mat3{float64(scaleX), 0, 0, 0, float64(scaleY), 0, 0, 0, 1}
+}
+
+// ShearX2D creates a homogeneous 2D shear matrix along the X-axis
+func ShearX2D(shear float64) Mat3 {
+	return Mat3{1, 0, 0, float64(shear), 1, 0, 0, 0, 1}
+}
+
+// ShearY2D creates a homogeneous 2D shear matrix along the Y-axis
+func ShearY2D(shear float64) Mat3 {
+	return Mat3{1, float64(shear), 0, 0, 1, 0, 0, 0, 1}
+}
+
+// ShearX3D creates a homogeneous 3D shear matrix along the X-axis
+func ShearX3D(shearY, shearZ float64) Mat4 {
+
+	return Mat4{1, float64(shearY), float64(shearZ), 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}
+}
+
+// ShearY3D creates a homogeneous 3D shear matrix along the Y-axis
+func ShearY3D(shearX, shearZ float64) Mat4 {
+	return Mat4{1, 0, 0, 0, float64(shearX), 1, float64(shearZ), 0, 0, 0, 1, 0, 0, 0, 0, 1}
+}
+
+// ShearZ3D creates a homogeneous 3D shear matrix along the Z-axis
+func ShearZ3D(shearX, shearY float64) Mat4 {
+	return Mat4{1, 0, 0, 0, 0, 1, 0, 0, float64(shearX), float64(shearY), 1, 0, 0, 0, 0, 1}
+}
+
+// HomogRotate3D creates a 3D rotation Matrix that rotates by (radian) angle about some arbitrary axis given by a Vector.
+// It produces a homogeneous matrix (4x4)
+//
+// Where c is cos(angle) and s is sin(angle), and x, y, and z are the first, second, and third elements of the axis vector (respectively):
+//
+//    [[ x^2(1-c)+c, xy(1-c)-zs, xz(1-c)+ys, 0 ]]
+//    [[ xy(1-c)+zs, y^2(1-c)+c, yz(1-c)-xs, 0 ]]
+//    [[ xz(1-c)-ys, yz(1-c)+xs, z^2(1-c)+c, 0 ]]
+//    [[ 0         , 0         , 0         , 1 ]]
+func HomogRotate3D(angle float64, axis Vec3) Mat4 {
+	x, y, z := axis[0], axis[1], axis[2]
+	s, c := float64(math.Sin(float64(angle))), float64(math.Cos(float64(angle)))
+	k := 1 - c
+
+	return Mat4{x*x*k + c, x*y*k + z*s, x*z*k - y*s, 0, x*y*k - z*s, y*y*k + c, y*z*k + x*s, 0, x*z*k + y*s, y*z*k - x*s, z*z*k + c, 0, 0, 0, 0, 1}
+}
+
+// Extracts the 3d scaling from a homogeneous matrix
+func Extract3DScale(m Mat4) (x, y, z float64) {
+	return float64(math.Sqrt(float64(m[0]*m[0] + m[1]*m[1] + m[2]*m[2]))),
+		float64(math.Sqrt(float64(m[4]*m[4] + m[5]*m[5] + m[6]*m[6]))),
+		float64(math.Sqrt(float64(m[8]*m[8] + m[9]*m[9] + m[10]*m[10])))
+}
+
+// Extracts the maximum scaling from a homogeneous matrix
+func ExtractMaxScale(m Mat4) float64 {
+	scaleX := float64(m[0]*m[0] + m[1]*m[1] + m[2]*m[2])
+	scaleY := float64(m[4]*m[4] + m[5]*m[5] + m[6]*m[6])
+	scaleZ := float64(m[8]*m[8] + m[9]*m[9] + m[10]*m[10])
+
+	return float64(math.Sqrt(math.Max(scaleX, math.Max(scaleY, scaleZ))))
+}
+
+// Calculates the Normal of the Matrix (aka the inverse transpose)
+func Mat4Normal(m Mat4) Mat3 {
+	n := m.Inv().Transpose()
+	return Mat3{n[0], n[1], n[2], n[4], n[5], n[6], n[8], n[9], n[10]}
+}
+
+// Multiplies a 3D vector by a transformation given by
+// the homogeneous 4D matrix m, applying any translation.
+// If this transformation is non-affine, it will project this
+// vector onto the plane w=1 before returning the result.
+//
+// This is similar to saying you're transforming and projecting a point.
+//
+// This is effectively equivalent to the GLSL code
+//     vec4 r = (m * vec4(v,1.));
+//     r = r/r.w;
+//     vec3 newV = r.xyz;
+func TransformCoordinate(v Vec3, m Mat4) Vec3 {
+	t := v.Vec4(1)
+	t = m.Mul4x1(t)
+	t = t.Mul(1 / t[3])
+
+	return t.Vec3()
+}
+
+// Multiplies a 3D vector by a transformation given by
+// the homogeneous 4D matrix m, NOT applying any translations.
+//
+// This is similar to saying you're applying a transformation
+// to a direction or normal. Rotation still applies (as does scaling),
+// but translating a direction or normal is meaningless.
+//
+// This is effectively equivalent to the GLSL code
+//    vec4 r = (m * vec4(v,0.));
+//    vec3 newV = r.xyz
+func TransformNormal(v Vec3, m Mat4) Vec3 {
+	t := v.Vec4(0)
+	t = m.Mul4x1(t)
+
+	return t.Vec3()
+}

vendor/github.com/go-gl/mathgl/mgl64/util.go

@@ -0,0 +1,144 @@
+// This file is generated from mgl32/util.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+//#go:generate go run codegen.go -template vector.tmpl -output vector.go
+//#go:generate go run codegen.go -template matrix.tmpl -output matrix.go
+//#go:generate go run codegen.go -mgl64
+
+package mgl64
+
+import (
+	"math"
+)
+
+// Epsilon is some tiny value that determines how precisely equal we want our floats to be
+// This is exported and left as a variable in case you want to change the default threshold for the
+// purposes of certain methods (e.g. Unproject uses the default epsilon when determining
+// if the determinant is "close enough" to zero to mean there's no inverse).
+//
+// This is, obviously, not mutex protected so be **absolutely sure** that no functions using Epsilon
+// are being executed when you change this.
+var Epsilon float64 = 1e-10
+
+// A direct copy of the math package's Abs. This is here for the mgl32
+// package, to prevent rampant type conversions during equality tests.
+func Abs(a float64) float64 {
+	if a < 0 {
+		return -a
+	} else if a == 0 {
+		return 0
+	}
+
+	return a
+}
+
+// FloatEqual is a safe utility function to compare floats.
+// It's Taken from http://floating-point-gui.de/errors/comparison/
+//
+// It is slightly altered to not call Abs when not needed.
+func FloatEqual(a, b float64) bool {
+	return FloatEqualThreshold(a, b, Epsilon)
+}
+
+// FloatEqualFunc is a utility closure that will generate a function that
+// always approximately compares floats like FloatEqualThreshold with a different
+// threshold.
+func FloatEqualFunc(epsilon float64) func(float64, float64) bool {
+	return func(a, b float64) bool {
+		return FloatEqualThreshold(a, b, epsilon)
+	}
+}
+
+var (
+	MinNormal = float64(1.1754943508222875e-38) // 1 / 2**(127 - 1)
+	MinValue  = float64(math.SmallestNonzeroFloat64)
+	MaxValue  = float64(math.MaxFloat64)
+
+	InfPos = float64(math.Inf(1))
+	InfNeg = float64(math.Inf(-1))
+	NaN    = float64(math.NaN())
+)
+
+// FloatEqualThreshold is a utility function to compare floats.
+// It's Taken from http://floating-point-gui.de/errors/comparison/
+//
+// It is slightly altered to not call Abs when not needed.
+//
+// This differs from FloatEqual in that it lets you pass in your comparison threshold, so that you can adjust the comparison value to your specific needs
+func FloatEqualThreshold(a, b, epsilon float64) bool {
+	if a == b { // Handles the case of inf or shortcuts the loop when no significant error has accumulated
+		return true
+	}
+
+	diff := Abs(a - b)
+	if a*b == 0 || diff < MinNormal { // If a or b are 0 or both are extremely close to it
+		return diff < epsilon*epsilon
+	}
+
+	// Else compare difference
+	return diff/(Abs(a)+Abs(b)) < epsilon
+}
+
+// Clamp takes in a value and two thresholds. If the value is smaller than the low
+// threshold, it returns the low threshold. If it's bigger than the high threshold
+// it returns the high threshold. Otherwise it returns the value.
+//
+// Useful to prevent some functions from freaking out because a value was
+// teeeeechnically out of range.
+func Clamp(a, low, high float64) float64 {
+	if a < low {
+		return low
+	} else if a > high {
+		return high
+	}
+
+	return a
+}
+
+// ClampFunc generates a closure that returns its parameter
+// clamped to the range [low,high].
+func ClampFunc(low, high float64) func(float64) float64 {
+	return func(a float64) float64 {
+		return Clamp(a, low, high)
+	}
+}
+
+/* The IsClamped functions use strict equality (meaning: not the FloatEqual function)
+there shouldn't be any major issues with this since clamp is often used to fix minor errors*/
+
+// Checks if a is clamped between low and high as if
+// Clamp(a, low, high) had been called.
+//
+// In most cases it's probably better to just call Clamp
+// without checking this since it's relatively cheap.
+func IsClamped(a, low, high float64) bool {
+	return a >= low && a <= high
+}
+
+// If a > b, then a will be set to the value of b.
+func SetMin(a, b *float64) {
+	if *b < *a {
+		*a = *b
+	}
+}
+
+// If a < b, then a will be set to the value of b.
+func SetMax(a, b *float64) {
+	if *a < *b {
+		*a = *b
+	}
+}
+
+// Round shortens a float32 value to a specified precision (number of digits after the decimal point)
+// with "round half up" tie-braking rule. Half-way values (23.5) are always rounded up (24).
+func Round(v float64, precision int) float64 {
+	p := float64(precision)
+	t := float64(v) * math.Pow(10, p)
+	if t > 0 {
+		return float64(math.Floor(t+0.5) / math.Pow(10, p))
+	}
+	return float64(math.Ceil(t-0.5) / math.Pow(10, p))
+}

vendor/github.com/go-gl/mathgl/mgl64/vecn.go

@@ -0,0 +1,352 @@
+// This file is generated from mgl32/vecn.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mgl64
+
+import (
+	"math"
+)
+
+// A vector of N elements backed by a slice
+//
+// As with MatMxN, this is not for hardcore linear algebra with large dimensions. Use github.com/gonum/matrix
+// or something like BLAS/LAPACK for that. This is for corner cases in 3D math where you require
+// something a little bigger that 4D, but still relatively small.
+//
+// This VecN uses several sync.Pool objects as a memory pool. The rule is that for any sized vector, the backing slice
+// has CAPACITY (not length) of 2^p where p is Ceil(log_2(N)) -- or in other words, rounding up the base-2
+// log of the size of the vector. E.G. a VecN of size 17 will have a backing slice of Cap 32.
+type VecN struct {
+	vec []float64
+}
+
+// Creates a new vector with a backing slice filled with the contents
+// of initial. It is NOT backed by initial, but rather a slice with cap
+// 2^p where p is Ceil(log_2(len(initial))), with the data from initial copied into
+// it.
+func NewVecNFromData(initial []float64) *VecN {
+	if initial == nil {
+		return &VecN{}
+	}
+	var internal []float64
+	if shouldPool {
+		internal = grabFromPool(len(initial))
+	} else {
+		internal = make([]float64, len(initial))
+	}
+	copy(internal, initial)
+	return &VecN{vec: internal}
+}
+
+// Creates a new vector with a backing slice of
+// 2^p where p = Ceil(log_2(n))
+func NewVecN(n int) *VecN {
+	if shouldPool {
+		return &VecN{vec: grabFromPool(n)}
+	} else {
+		return &VecN{vec: make([]float64, n)}
+	}
+}
+
+// Returns the raw slice backing the VecN
+//
+// This may be sent back to the memory pool at any time
+// and you aren't advised to rely on this value
+func (vn VecN) Raw() []float64 {
+	return vn.vec
+}
+
+// Gets the element at index i from the vector.
+// This does not bounds check, and will panic if i is
+// out of range.
+func (vn VecN) Get(i int) float64 {
+	return vn.vec[i]
+}
+
+func (vn *VecN) Set(i int, val float64) {
+	vn.vec[i] = val
+}
+
+// Sends the allocated memory through the callback if it exists
+func (vn *VecN) destroy() {
+	if vn == nil || vn.vec == nil {
+		return
+	}
+
+	if shouldPool {
+		returnToPool(vn.vec)
+	}
+	vn.vec = nil
+}
+
+// Resizes the underlying slice to the desired amount, reallocating or retrieving from the pool
+// if necessary. The values after a Resize cannot be expected to be related to the values before a Resize.
+//
+// If the caller is a nil pointer, this returns a value as if NewVecN(n) had been called,
+// otherwise it simply returns the caller.
+func (vn *VecN) Resize(n int) *VecN {
+	if vn == nil {
+		return NewVecN(n)
+	}
+
+	if n <= cap(vn.vec) {
+		if vn.vec != nil {
+			vn.vec = vn.vec[:n]
+		} else {
+			vn.vec = []float64{}
+		}
+		return vn
+	}
+
+	if shouldPool && vn.vec != nil {
+		returnToPool(vn.vec)
+	}
+	*vn = (*NewVecN(n))
+
+	return vn
+}
+
+// Sets the vector's backing slice to the given
+// new one.
+func (vn *VecN) SetBackingSlice(newSlice []float64) {
+	vn.vec = newSlice
+}
+
+// Return the len of the vector's underlying slice.
+// This is not titled Len because it conflicts the package's
+// convention of calling the Norm the Len.
+func (vn *VecN) Size() int {
+	return len(vn.vec)
+}
+
+// Returns the cap of the vector's underlying slice.
+func (vn *VecN) Cap() int {
+	return cap(vn.vec)
+}
+
+// Sets the vector's size to n and zeroes out the vector.
+// If n is bigger than the vector's size, it will realloc.
+func (vn *VecN) Zero(n int) {
+	vn.Resize(n)
+	for i := range vn.vec {
+		vn.vec[i] = 0
+	}
+}
+
+// Adds vn and addend, storing the result in dst.
+// If dst does not have sufficient size it will be resized
+// Dst may be one of the other arguments. If dst is nil, it will be allocated.
+// The value returned is dst, for easier method chaining
+//
+// If vn and addend are not the same size, this function will add min(vn.Size(), addend.Size())
+// elements.
+func (vn *VecN) Add(dst *VecN, subtrahend *VecN) *VecN {
+	if vn == nil || subtrahend == nil {
+		return nil
+	}
+	size := intMin(len(vn.vec), len(subtrahend.vec))
+	dst = dst.Resize(size)
+
+	for i := 0; i < size; i++ {
+		dst.vec[i] = vn.vec[i] + subtrahend.vec[i]
+	}
+
+	return dst
+}
+
+// Subtracts addend from vn, storing the result in dst.
+// If dst does not have sufficient size it will be resized
+// Dst may be one of the other arguments. If dst is nil, it will be allocated.
+// The value returned is dst, for easier method chaining
+//
+// If vn and addend are not the same size, this function will add min(vn.Size(), addend.Size())
+// elements.
+func (vn *VecN) Sub(dst *VecN, addend *VecN) *VecN {
+	if vn == nil || addend == nil {
+		return nil
+	}
+	size := intMin(len(vn.vec), len(addend.vec))
+	dst = dst.Resize(size)
+
+	for i := 0; i < size; i++ {
+		dst.vec[i] = vn.vec[i] - addend.vec[i]
+	}
+
+	return dst
+}
+
+// Takes the binary cross product of vn and other, and stores it in dst.
+// If either vn or other are not of size 3 this function will panic
+//
+// If dst is not of sufficient size, or is nil, a new slice is allocated.
+// Dst is permitted to be one of the other arguments
+func (vn *VecN) Cross(dst *VecN, other *VecN) *VecN {
+	if vn == nil || other == nil {
+		return nil
+	}
+	if len(vn.vec) != 3 || len(other.vec) != 3 {
+		panic("Cannot take binary cross product of non-3D elements (7D cross product not implemented)")
+	}
+
+	dst = dst.Resize(3)
+	dst.vec[0], dst.vec[1], dst.vec[2] = vn.vec[1]*other.vec[2]-vn.vec[2]*other.vec[1], vn.vec[2]*other.vec[0]-vn.vec[0]*other.vec[2], vn.vec[0]*other.vec[1]-vn.vec[1]*other.vec[0]
+
+	return dst
+}
+
+func intMin(a, b int) int {
+	if a < b {
+		return a
+	}
+
+	return b
+}
+
+// Computes the dot product of two VecNs, if
+// the two vectors are not of the same length -- this
+// will return NaN.
+func (vn *VecN) Dot(other *VecN) float64 {
+	if vn == nil || other == nil || len(vn.vec) != len(other.vec) {
+		return float64(math.NaN())
+	}
+
+	var result float64 = 0.0
+	for i, el := range vn.vec {
+		result += el * other.vec[i]
+	}
+
+	return result
+}
+
+// Computes the vector length (also called the Norm) of the
+// vector. Equivalent to math.Sqrt(vn.Dot(vn)) with the appropriate
+// type conversions.
+//
+// If vn is nil, this returns NaN
+func (vn *VecN) Len() float64 {
+	if vn == nil {
+		return float64(math.NaN())
+	}
+	if len(vn.vec) == 0 {
+		return 0
+	}
+
+	return float64(math.Sqrt(float64(vn.Dot(vn))))
+}
+
+// Normalizes the vector and stores the result in dst, which
+// will be returned. Dst will be appropraitely resized to the
+// size of vn.
+//
+// The destination can be vn itself and nothing will go wrong.
+//
+// This is equivalent to vn.Mul(dst, 1/vn.Len())
+func (vn *VecN) Normalize(dst *VecN) *VecN {
+	if vn == nil {
+		return nil
+	}
+
+	return vn.Mul(dst, 1/vn.Len())
+}
+
+// Multiplied the vector by some scalar value and stores the result in dst, which
+// will be returned. Dst will be appropraitely resized to the
+// size of vn.
+//
+// The destination can be vn itself and nothing will go wrong.
+func (vn *VecN) Mul(dst *VecN, c float64) *VecN {
+	if vn == nil {
+		return nil
+	}
+	dst = dst.Resize(len(vn.vec))
+
+	for i, el := range vn.vec {
+		dst.vec[i] = el * c
+	}
+
+	return dst
+}
+
+// Performs the vector outer product between vn and v2.
+// The outer product is like a "reverse" dot product. Where the dot product
+// aligns both vectors with the "sized" part facing "inward" (Vec3*Vec3=Mat1x3*Mat3x1=Mat1x1=Scalar).
+// The outer product multiplied them with it facing "outward"
+// (Vec3*Vec3=Mat3x1*Mat1x3=Mat3x3).
+//
+// The matrix dst will be Reshaped to the correct size, if vn or v2 are nil,
+// this returns nil.
+func (vn *VecN) OuterProd(dst *MatMxN, v2 *VecN) *MatMxN {
+	if vn == nil || v2 == nil {
+		return nil
+	}
+
+	dst = dst.Reshape(len(vn.vec), len(v2.vec))
+
+	for c, el1 := range v2.vec {
+		for r, el2 := range vn.vec {
+			dst.Set(r, c, el1*el2)
+		}
+	}
+
+	return dst
+}
+
+func (vn *VecN) ApproxEqual(vn2 *VecN) bool {
+	if vn == nil || vn2 == nil || len(vn.vec) != len(vn2.vec) {
+		return false
+	}
+
+	for i, el := range vn.vec {
+		if !FloatEqual(el, vn2.vec[i]) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (vn *VecN) ApproxEqualThreshold(vn2 *VecN, epsilon float64) bool {
+	if vn == nil || vn2 == nil || len(vn.vec) != len(vn2.vec) {
+		return false
+	}
+
+	for i, el := range vn.vec {
+		if !FloatEqualThreshold(el, vn2.vec[i], epsilon) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (vn *VecN) ApproxEqualFunc(vn2 *VecN, comp func(float64, float64) bool) bool {
+	if vn == nil || vn2 == nil || len(vn.vec) != len(vn2.vec) {
+		return false
+	}
+
+	for i, el := range vn.vec {
+		if !comp(el, vn2.vec[i]) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (vn *VecN) Vec2() Vec2 {
+	raw := vn.Raw()
+	return Vec2{raw[0], raw[1]}
+}
+
+func (vn *VecN) Vec3() Vec3 {
+	raw := vn.Raw()
+	return Vec3{raw[0], raw[1], raw[2]}
+}
+
+func (vn *VecN) Vec4() Vec4 {
+	raw := vn.Raw()
+	return Vec4{raw[0], raw[1], raw[2], raw[3]}
+}

vendor/github.com/go-gl/mathgl/mgl64/vector.go

@@ -0,0 +1,565 @@
+// This file is generated from mgl32/vector.go; DO NOT EDIT
+
+// Copyright 2014 The go-gl/mathgl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// This file is generated by codegen.go; DO NOT EDIT
+// Edit vector.tmpl and run "go generate" to make changes.
+
+package mgl64
+
+import (
+	"math"
+
+	"golang.org/x/image/math/f64"
+)
+
+type Vec2 f64.Vec2
+type Vec3 f64.Vec3
+type Vec4 f64.Vec4
+
+func (v Vec2) Vec3(z float64) Vec3 {
+	return Vec3{v[0], v[1], z}
+}
+
+func (v Vec2) Vec4(z, w float64) Vec4 {
+	return Vec4{v[0], v[1], z, w}
+}
+
+func (v Vec3) Vec2() Vec2 {
+	return Vec2{v[0], v[1]}
+}
+
+func (v Vec3) Vec4(w float64) Vec4 {
+	return Vec4{v[0], v[1], v[2], w}
+}
+
+func (v Vec4) Vec2() Vec2 {
+	return Vec2{v[0], v[1]}
+}
+
+func (v Vec4) Vec3() Vec3 {
+	return Vec3{v[0], v[1], v[2]}
+}
+
+// Elem extracts the elements of the vector for direct value assignment.
+func (v Vec2) Elem() (x, y float64) {
+	return v[0], v[1]
+}
+
+// Elem extracts the elements of the vector for direct value assignment.
+func (v Vec3) Elem() (x, y, z float64) {
+	return v[0], v[1], v[2]
+}
+
+// Elem extracts the elements of the vector for direct value assignment.
+func (v Vec4) Elem() (x, y, z, w float64) {
+	return v[0], v[1], v[2], v[3]
+}
+
+// The vector cross product is an operation only defined on 3D vectors. It is equivalent to
+// Vec3{v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1] - v1[1]*v2[0]}.
+// Another interpretation is that it's the vector whose magnitude is |v1||v2|sin(theta)
+// where theta is the angle between v1 and v2.
+//
+// The cross product is most often used for finding surface normals. The cross product of vectors
+// will generate a vector that is perpendicular to the plane they form.
+//
+// Technically, a generalized cross product exists as an "(N-1)ary" operation
+// (that is, the 4D cross product requires 3 4D vectors). But the binary
+// 3D (and 7D) cross product is the most important. It can be considered
+// the area of a parallelogram with sides v1 and v2.
+//
+// Like the dot product, the cross product is roughly a measure of directionality.
+// Two normalized perpendicular vectors will return a vector with a magnitude of
+// 1.0 or -1.0 and two parallel vectors will return a vector with magnitude 0.0.
+// The cross product is "anticommutative" meaning v1.Cross(v2) = -v2.Cross(v1),
+// this property can be useful to know when finding normals,
+// as taking the wrong cross product can lead to the opposite normal of the one you want.
+func (v1 Vec3) Cross(v2 Vec3) Vec3 {
+	return Vec3{v1[1]*v2[2] - v1[2]*v2[1], v1[2]*v2[0] - v1[0]*v2[2], v1[0]*v2[1] - v1[1]*v2[0]}
+}
+
+// Add performs element-wise addition between two vectors. It is equivalent to iterating
+// over every element of v1 and adding the corresponding element of v2 to it.
+func (v1 Vec2) Add(v2 Vec2) Vec2 {
+	return Vec2{v1[0] + v2[0], v1[1] + v2[1]}
+}
+
+// Sub performs element-wise subtraction between two vectors. It is equivalent to iterating
+// over every element of v1 and subtracting the corresponding element of v2 from it.
+func (v1 Vec2) Sub(v2 Vec2) Vec2 {
+	return Vec2{v1[0] - v2[0], v1[1] - v2[1]}
+}
+
+// Mul performs a scalar multiplication between the vector and some constant value
+// c. This is equivalent to iterating over every vector element and multiplying by c.
+func (v1 Vec2) Mul(c float64) Vec2 {
+	return Vec2{v1[0] * c, v1[1] * c}
+}
+
+// Dot returns the dot product of this vector with another. There are multiple ways
+// to describe this value. One is the multiplication of their lengths and cos(theta) where
+// theta is the angle between the vectors: v1.v2 = |v1||v2|cos(theta).
+//
+// The other (and what is actually done) is the sum of the element-wise multiplication of all
+// elements. So for instance, two Vec3s would yield v1.x * v2.x + v1.y * v2.y + v1.z * v2.z.
+//
+// This means that the dot product of a vector and itself is the square of its Len (within
+// the bounds of floating points error).
+//
+// The dot product is roughly a measure of how closely two vectors are to pointing in the same
+// direction. If both vectors are normalized, the value will be -1 for opposite pointing,
+// one for same pointing, and 0 for perpendicular vectors.
+func (v1 Vec2) Dot(v2 Vec2) float64 {
+	return v1[0]*v2[0] + v1[1]*v2[1]
+}
+
+// Len returns the vector's length. Note that this is NOT the dimension of
+// the vector (len(v)), but the mathematical length. This is equivalent to the square
+// root of the sum of the squares of all elements. E.G. for a Vec2 it's
+// math.Hypot(v[0], v[1]).
+func (v1 Vec2) Len() float64 {
+
+	return float64(math.Hypot(float64(v1[0]), float64(v1[1])))
+
+}
+
+// Normalize normalizes the vector. Normalization is (1/|v|)*v,
+// making this equivalent to v.Scale(1/v.Len()). If the len is 0.0,
+// this function will return an infinite value for all elements due
+// to how floating point division works in Go (n/0.0 = math.Inf(Sign(n))).
+//
+// Normalization makes a vector's Len become 1.0 (within the margin of floating point error),
+// while maintaining its directionality.
+//
+// (Can be seen here: http://play.golang.org/p/Aaj7SnbqIp )
+func (v1 Vec2) Normalize() Vec2 {
+	l := 1.0 / v1.Len()
+	return Vec2{v1[0] * l, v1[1] * l}
+}
+
+// ApproxEqual takes in a vector and does an element-wise
+// approximate float comparison as if FloatEqual had been used
+func (v1 Vec2) ApproxEqual(v2 Vec2) bool {
+	for i := range v1 {
+		if !FloatEqual(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxThresholdEq takes in a threshold for comparing two floats, and uses it to do an
+// element-wise comparison of the vector to another.
+func (v1 Vec2) ApproxEqualThreshold(v2 Vec2, threshold float64) bool {
+	for i := range v1 {
+		if !FloatEqualThreshold(v1[i], v2[i], threshold) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxFuncEq takes in a func that compares two floats, and uses it to do an element-wise
+// comparison of the vector to another. This is intended to be used with FloatEqualFunc
+func (v1 Vec2) ApproxFuncEqual(v2 Vec2, eq func(float64, float64) bool) bool {
+	for i := range v1 {
+		if !eq(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec2) X() float64 {
+	return v[0]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec2) Y() float64 {
+	return v[1]
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 2x2 matrix. E.G. a Vec2 * Vec2 = Mat2.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec2) OuterProd2(v2 Vec2) Mat2 {
+	return Mat2{v1[0] * v2[0], v1[1] * v2[0], v1[0] * v2[1], v1[1] * v2[1]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 2x3 matrix. E.G. a Vec2 * Vec3 = Mat2x3.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec2) OuterProd3(v2 Vec3) Mat2x3 {
+	return Mat2x3{v1[0] * v2[0], v1[1] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[0] * v2[2], v1[1] * v2[2]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 2x4 matrix. E.G. a Vec2 * Vec4 = Mat2x4.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec2) OuterProd4(v2 Vec4) Mat2x4 {
+	return Mat2x4{v1[0] * v2[0], v1[1] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[0] * v2[3], v1[1] * v2[3]}
+}
+
+// Add performs element-wise addition between two vectors. It is equivalent to iterating
+// over every element of v1 and adding the corresponding element of v2 to it.
+func (v1 Vec3) Add(v2 Vec3) Vec3 {
+	return Vec3{v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2]}
+}
+
+// Sub performs element-wise subtraction between two vectors. It is equivalent to iterating
+// over every element of v1 and subtracting the corresponding element of v2 from it.
+func (v1 Vec3) Sub(v2 Vec3) Vec3 {
+	return Vec3{v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2]}
+}
+
+// Mul performs a scalar multiplication between the vector and some constant value
+// c. This is equivalent to iterating over every vector element and multiplying by c.
+func (v1 Vec3) Mul(c float64) Vec3 {
+	return Vec3{v1[0] * c, v1[1] * c, v1[2] * c}
+}
+
+// Dot returns the dot product of this vector with another. There are multiple ways
+// to describe this value. One is the multiplication of their lengths and cos(theta) where
+// theta is the angle between the vectors: v1.v2 = |v1||v2|cos(theta).
+//
+// The other (and what is actually done) is the sum of the element-wise multiplication of all
+// elements. So for instance, two Vec3s would yield v1.x * v2.x + v1.y * v2.y + v1.z * v2.z.
+//
+// This means that the dot product of a vector and itself is the square of its Len (within
+// the bounds of floating points error).
+//
+// The dot product is roughly a measure of how closely two vectors are to pointing in the same
+// direction. If both vectors are normalized, the value will be -1 for opposite pointing,
+// one for same pointing, and 0 for perpendicular vectors.
+func (v1 Vec3) Dot(v2 Vec3) float64 {
+	return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]
+}
+
+// Len returns the vector's length. Note that this is NOT the dimension of
+// the vector (len(v)), but the mathematical length. This is equivalent to the square
+// root of the sum of the squares of all elements. E.G. for a Vec2 it's
+// math.Hypot(v[0], v[1]).
+func (v1 Vec3) Len() float64 {
+
+	return float64(math.Sqrt(float64(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2])))
+
+}
+
+// Normalize normalizes the vector. Normalization is (1/|v|)*v,
+// making this equivalent to v.Scale(1/v.Len()). If the len is 0.0,
+// this function will return an infinite value for all elements due
+// to how floating point division works in Go (n/0.0 = math.Inf(Sign(n))).
+//
+// Normalization makes a vector's Len become 1.0 (within the margin of floating point error),
+// while maintaining its directionality.
+//
+// (Can be seen here: http://play.golang.org/p/Aaj7SnbqIp )
+func (v1 Vec3) Normalize() Vec3 {
+	l := 1.0 / v1.Len()
+	return Vec3{v1[0] * l, v1[1] * l, v1[2] * l}
+}
+
+// ApproxEqual takes in a vector and does an element-wise
+// approximate float comparison as if FloatEqual had been used
+func (v1 Vec3) ApproxEqual(v2 Vec3) bool {
+	for i := range v1 {
+		if !FloatEqual(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxThresholdEq takes in a threshold for comparing two floats, and uses it to do an
+// element-wise comparison of the vector to another.
+func (v1 Vec3) ApproxEqualThreshold(v2 Vec3, threshold float64) bool {
+	for i := range v1 {
+		if !FloatEqualThreshold(v1[i], v2[i], threshold) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxFuncEq takes in a func that compares two floats, and uses it to do an element-wise
+// comparison of the vector to another. This is intended to be used with FloatEqualFunc
+func (v1 Vec3) ApproxFuncEqual(v2 Vec3, eq func(float64, float64) bool) bool {
+	for i := range v1 {
+		if !eq(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec3) X() float64 {
+	return v[0]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec3) Y() float64 {
+	return v[1]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec3) Z() float64 {
+	return v[2]
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 3x2 matrix. E.G. a Vec3 * Vec2 = Mat3x2.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec3) OuterProd2(v2 Vec2) Mat3x2 {
+	return Mat3x2{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 3x3 matrix. E.G. a Vec3 * Vec3 = Mat3.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec3) OuterProd3(v2 Vec3) Mat3 {
+	return Mat3{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[2] * v2[2]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 3x4 matrix. E.G. a Vec3 * Vec4 = Mat3x4.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec3) OuterProd4(v2 Vec4) Mat3x4 {
+	return Mat3x4{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[2] * v2[2], v1[0] * v2[3], v1[1] * v2[3], v1[2] * v2[3]}
+}
+
+// Add performs element-wise addition between two vectors. It is equivalent to iterating
+// over every element of v1 and adding the corresponding element of v2 to it.
+func (v1 Vec4) Add(v2 Vec4) Vec4 {
+	return Vec4{v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2], v1[3] + v2[3]}
+}
+
+// Sub performs element-wise subtraction between two vectors. It is equivalent to iterating
+// over every element of v1 and subtracting the corresponding element of v2 from it.
+func (v1 Vec4) Sub(v2 Vec4) Vec4 {
+	return Vec4{v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2], v1[3] - v2[3]}
+}
+
+// Mul performs a scalar multiplication between the vector and some constant value
+// c. This is equivalent to iterating over every vector element and multiplying by c.
+func (v1 Vec4) Mul(c float64) Vec4 {
+	return Vec4{v1[0] * c, v1[1] * c, v1[2] * c, v1[3] * c}
+}
+
+// Dot returns the dot product of this vector with another. There are multiple ways
+// to describe this value. One is the multiplication of their lengths and cos(theta) where
+// theta is the angle between the vectors: v1.v2 = |v1||v2|cos(theta).
+//
+// The other (and what is actually done) is the sum of the element-wise multiplication of all
+// elements. So for instance, two Vec3s would yield v1.x * v2.x + v1.y * v2.y + v1.z * v2.z.
+//
+// This means that the dot product of a vector and itself is the square of its Len (within
+// the bounds of floating points error).
+//
+// The dot product is roughly a measure of how closely two vectors are to pointing in the same
+// direction. If both vectors are normalized, the value will be -1 for opposite pointing,
+// one for same pointing, and 0 for perpendicular vectors.
+func (v1 Vec4) Dot(v2 Vec4) float64 {
+	return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] + v1[3]*v2[3]
+}
+
+// Len returns the vector's length. Note that this is NOT the dimension of
+// the vector (len(v)), but the mathematical length. This is equivalent to the square
+// root of the sum of the squares of all elements. E.G. for a Vec2 it's
+// math.Hypot(v[0], v[1]).
+func (v1 Vec4) Len() float64 {
+
+	return float64(math.Sqrt(float64(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2] + v1[3]*v1[3])))
+
+}
+
+// Normalize normalizes the vector. Normalization is (1/|v|)*v,
+// making this equivalent to v.Scale(1/v.Len()). If the len is 0.0,
+// this function will return an infinite value for all elements due
+// to how floating point division works in Go (n/0.0 = math.Inf(Sign(n))).
+//
+// Normalization makes a vector's Len become 1.0 (within the margin of floating point error),
+// while maintaining its directionality.
+//
+// (Can be seen here: http://play.golang.org/p/Aaj7SnbqIp )
+func (v1 Vec4) Normalize() Vec4 {
+	l := 1.0 / v1.Len()
+	return Vec4{v1[0] * l, v1[1] * l, v1[2] * l, v1[3] * l}
+}
+
+// ApproxEqual takes in a vector and does an element-wise
+// approximate float comparison as if FloatEqual had been used
+func (v1 Vec4) ApproxEqual(v2 Vec4) bool {
+	for i := range v1 {
+		if !FloatEqual(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxThresholdEq takes in a threshold for comparing two floats, and uses it to do an
+// element-wise comparison of the vector to another.
+func (v1 Vec4) ApproxEqualThreshold(v2 Vec4, threshold float64) bool {
+	for i := range v1 {
+		if !FloatEqualThreshold(v1[i], v2[i], threshold) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxFuncEq takes in a func that compares two floats, and uses it to do an element-wise
+// comparison of the vector to another. This is intended to be used with FloatEqualFunc
+func (v1 Vec4) ApproxFuncEqual(v2 Vec4, eq func(float64, float64) bool) bool {
+	for i := range v1 {
+		if !eq(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec4) X() float64 {
+	return v[0]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec4) Y() float64 {
+	return v[1]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec4) Z() float64 {
+	return v[2]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec4) W() float64 {
+	return v[3]
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 4x2 matrix. E.G. a Vec4 * Vec2 = Mat4x2.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec4) OuterProd2(v2 Vec2) Mat4x2 {
+	return Mat4x2{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[3] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[3] * v2[1]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 4x3 matrix. E.G. a Vec4 * Vec3 = Mat4x3.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec4) OuterProd3(v2 Vec3) Mat4x3 {
+	return Mat4x3{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[3] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[3] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[2] * v2[2], v1[3] * v2[2]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 4x4 matrix. E.G. a Vec4 * Vec4 = Mat4.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec4) OuterProd4(v2 Vec4) Mat4 {
+	return Mat4{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[3] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[3] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[2] * v2[2], v1[3] * v2[2], v1[0] * v2[3], v1[1] * v2[3], v1[2] * v2[3], v1[3] * v2[3]}
+}