func Slerp(a, b *T, t float64) T {
	d := math.Acos(a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + a[3]*b[3])
	ooSinD := 1 / math.Sin(d)
	t1 := math.Sin(d*(1-t)) * ooSinD
	t2 := math.Sin(d*t) * ooSinD
	q := T{
		a[0]*t1 + b[0]*t2,
		a[1]*t1 + b[1]*t2,
		a[2]*t1 + b[2]*t2,
		a[3]*t1 + b[3]*t2,
	}
	return q.Normalized()
}
Slerp returns the spherical linear interpolation quaternion between a and b at t (0,1).
See http://en.wikipedia.org/wiki/Slerp
package quaternion
import (
	"fmt"
	"math"
	"github.com/ungerik/go3d/float64/vec3"
	"github.com/ungerik/go3d/float64/vec4"
)
var (
	Zero = T{}
	Ident = T{0, 0, 0, 1}
)
type T [4]float64
func FromAxisAngle(axis *vec3.T, angle float64) T {
	angle *= 0.5
	sin := math.Sin(angle)
	q := T{axis[0] * sin, axis[1] * sin, axis[2] * sin, math.Cos(angle)}
	return q.Normalized()
}
func FromXAxisAngle(angle float64) T {
	angle *= 0.5
	return T{math.Sin(angle), 0, 0, math.Cos(angle)}
}
func FromYAxisAngle(angle float64) T {
	angle *= 0.5
	return T{0, math.Sin(angle), 0, math.Cos(angle)}
}
func FromZAxisAngle(angle float64) T {
	angle *= 0.5
	return T{0, 0, math.Sin(angle), math.Cos(angle)}
}
func FromEulerAngles(yHead, xPitch, zRoll float64) T {
	qy := FromYAxisAngle(yHead)
	qx := FromXAxisAngle(xPitch)
	qz := FromZAxisAngle(zRoll)
	return Mul3(&qy, &qx, &qz)
}
func (quat *T) ToEulerAngles() (yHead, xPitch, zRoll float64) {
	z := quat.RotatedVec3(&vec3.T{0, 0, 1})
	yHead = math.Atan2(z[0], z[2])
	xPitch = -math.Atan2(z[1], math.Sqrt(z[0]*z[0]+z[2]*z[2]))
	quatNew := FromEulerAngles(yHead, xPitch, 0)
	x2 := quatNew.RotatedVec3(&vec3.T{1, 0, 0})
	x := quat.RotatedVec3(&vec3.T{1, 0, 0})
	r := vec3.Cross(&x, &x2)
	sin := vec3.Dot(&r, &z)
	cos := vec3.Dot(&x, &x2)
	zRoll = -math.Atan2(sin, cos)
	return
}
func FromVec4(v *vec4.T) T {
	return T(*v)
}
func (quat *T) Vec4() vec4.T {
	return vec4.T(*quat)
}
func Parse(s string) (r T, err error) {
	_, err = fmt.Sscan(s, &r[0], &r[1], &r[2], &r[3])
	return r, err
}
func (quat *T) String() string {
	return fmt.Sprint(quat[0], quat[1], quat[2], quat[3])
}
func (quat *T) AxisAngle() (axis vec3.T, angle float64) {
	cos := quat[3]
	sin := math.Sqrt(1 - cos*cos)
	angle = math.Acos(cos) * 2
	var ooSin float64
	if math.Abs(sin) < 0.0005 {
		ooSin = 1
	} else {
		ooSin = 1 / sin
	}
	axis[0] = quat[0] * ooSin
	axis[1] = quat[1] * ooSin
	axis[2] = quat[2] * ooSin
	return axis, angle
}
func (quat *T) Norm() float64 {
	return quat[0]*quat[0] + quat[1]*quat[1] + quat[2]*quat[2] + quat[3]*quat[3]
}
func (quat *T) Normalize() *T {
	norm := quat.Norm()
	if norm != 1 && norm != 0 {
		ool := 1 / math.Sqrt(norm)
		quat[0] *= ool
		quat[1] *= ool
		quat[2] *= ool
		quat[3] *= ool
	}
	return quat
}
func (quat *T) Normalized() T {
	norm := quat.Norm()
	if norm != 1 && norm != 0 {
		ool := 1 / math.Sqrt(norm)
		return T{
			quat[0] * ool,
			quat[1] * ool,
			quat[2] * ool,
			quat[3] * ool,
		}
	} else {
		return *quat
	}
}
func (quat *T) Negate() *T {
	quat[0] = -quat[0]
	quat[1] = -quat[1]
	quat[2] = -quat[2]
	quat[3] = -quat[3]
	return quat
}
func (quat *T) Negated() T {
	return T{-quat[0], -quat[1], -quat[2], -quat[3]}
}
func (quat *T) Invert() *T {
	quat[0] = -quat[0]
	quat[1] = -quat[1]
	quat[2] = -quat[2]
	return quat
}
func (quat *T) Inverted() T {
	return T{-quat[0], -quat[1], -quat[2], quat[3]}
}
func (quat *T) SetShortestRotation(other *T) *T {
	if !IsShortestRotation(quat, other) {
		quat.Negate()
	}
	return quat
}
func IsShortestRotation(a, b *T) bool {
	return Dot(a, b) >= 0
}
func (quat *T) IsUnitQuat(tolerance float64) bool {
	norm := quat.Norm()
	return norm >= (1.0-tolerance) && norm <= (1.0+tolerance)
}
func (quat *T) RotateVec3(v *vec3.T) {
	u := vec3.T{quat[0], quat[1], quat[2]}
	s := quat[3]
	vt1 := u.Scaled(2 * vec3.Dot(&u, v))
	vt2 := v.Scaled(s*s - vec3.Dot(&u, &u))
	vt3 := vec3.Cross(&u, v)
	vt3 = vt3.Scaled(2 * s)
	v[0] = vt1[0] + vt2[0] + vt3[0]
	v[1] = vt1[1] + vt2[1] + vt3[1]
	v[2] = vt1[2] + vt2[2] + vt3[2]
}
func (quat *T) RotatedVec3(v *vec3.T) vec3.T {
	u := vec3.T{quat[0], quat[1], quat[2]}
	s := quat[3]
	vt1 := u.Scaled(2 * vec3.Dot(&u, v))
	vt2 := v.Scaled(s*s - vec3.Dot(&u, &u))
	vt3 := vec3.Cross(&u, v)
	vt3 = vt3.Scaled(2 * s)
	return vec3.T{vt1[0] + vt2[0] + vt3[0], vt1[1] + vt2[1] + vt3[1], vt1[2] + vt2[2] + vt3[2]}
}
func Dot(a, b *T) float64 {
	return a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + a[3]*b[3]
}
func Mul(a, b *T) T {
	q := T{
		a[3]*b[0] + a[0]*b[3] + a[1]*b[2] - a[2]*b[1],
		a[3]*b[1] + a[1]*b[3] + a[2]*b[0] - a[0]*b[2],
		a[3]*b[2] + a[2]*b[3] + a[0]*b[1] - a[1]*b[0],
		a[3]*b[3] - a[0]*b[0] - a[1]*b[1] - a[2]*b[2],
	}
	return q.Normalized()
}
func Mul3(a, b, c *T) T {
	q := Mul(a, b)
	return Mul(&q, c)
}
func Mul4(a, b, c, d *T) T {
	q := Mul(a, b)
	q = Mul(&q, c)
	return Mul(&q, d)
}