diff --git a/gglm/constants.go b/gglm/constants.go
new file mode 100755
index 0000000..94a0439
--- /dev/null
+++ b/gglm/constants.go
@@ -0,0 +1,7 @@
+package gglm
+
+const (
+	Pi      float32 = 3.14159265359
+	Deg2Rad float32 = Pi / 180
+	Rad2Deg float32 = 180 / Pi
+)
diff --git a/gglm/quat.go b/gglm/quat.go
new file mode 100755
index 0000000..3ce4ea3
--- /dev/null
+++ b/gglm/quat.go
@@ -0,0 +1,63 @@
+package gglm
+
+import (
+	"fmt"
+)
+
+var _ Swizzle4 = &Quat{}
+var _ fmt.Stringer = &Quat{}
+
+type Quat struct {
+	Vec4
+}
+
+//Euler takes rotations in radians and produces a rotation that
+//rotates around the z-axis, y-axis and lastly x-axis.
+func NewQuatEuler(v *Vec3) *Quat {
+
+	//Some other common terminology: x=roll, y=pitch, z=yaw
+	sinX, cosX := Sincos32(v.Data[0] * 0.5)
+	sinY, cosY := Sincos32(v.Data[1] * 0.5)
+	sinZ, cosZ := Sincos32(v.Data[2] * 0.5)
+
+	//This produces a z->y->x multiply order, but its written as XYZ.
+	//This is due to XYZ meaning independent rotation matrices, so Z is applied
+	//first, then Y matrix and lastly X.
+	//See this for more info: https://github.com/godotengine/godot/issues/6816#issuecomment-254592170
+	//
+	//Note: On most conversion tools putting the multiply order (e.g. ZYX for us) is required.
+	return &Quat{
+		Vec4: Vec4{
+			Data: [4]float32{
+				sinX*cosY*cosZ - cosX*sinY*sinZ,
+				cosX*sinY*cosZ + sinX*cosY*sinZ,
+				cosX*cosY*sinZ - sinX*sinY*cosZ,
+				cosX*cosY*cosZ + sinX*sinY*sinZ,
+			},
+		},
+	}
+}
+
+//AngleAxis rotates rotRad radians around the *normalized* vector rotAxisNorm
+func NewQuatAngleAxis(rotRad float32, rotAxisNorm *Vec3) *Quat {
+
+	s, c := Sincos32(rotRad * 0.5)
+	return &Quat{
+		Vec4: Vec4{
+			Data: [4]float32{
+				rotAxisNorm.Data[0] * s,
+				rotAxisNorm.Data[1] * s,
+				rotAxisNorm.Data[2] * s,
+				c,
+			},
+		},
+	}
+}
+
+func NewQuatId() *Quat {
+	return &Quat{
+		Vec4: Vec4{
+			Data: [4]float32{0, 0, 0, 1},
+		},
+	}
+}
diff --git a/gglm/scalar.go b/gglm/scalar.go
index 69c3e02..0a7bb25 100755
--- a/gglm/scalar.go
+++ b/gglm/scalar.go
@@ -14,3 +14,16 @@ func EqF32(f1, f2 float32) bool {
 func EqF32Epsilon(f1, f2, eps float32) bool {
 	return math.Abs(float64(f1-f2)) <= float64(eps)
 }
+
+func Sin32(x float32) float32 {
+	return float32(math.Sin(float64(x)))
+}
+
+func Cos32(x float32) float32 {
+	return float32(math.Cos(float64(x)))
+}
+
+func Sincos32(x float32) (sinx, cosx float32) {
+	a, b := math.Sincos(float64(x))
+	return float32(a), float32(b)
+}
diff --git a/gglm/transform.go b/gglm/transform.go
new file mode 100755
index 0000000..a2c4e79
--- /dev/null
+++ b/gglm/transform.go
@@ -0,0 +1,67 @@
+package gglm
+
+//TMat is the 4x4 transformation matrix
+// type TMat struct {
+// 	Mat4
+// }
+
+func GenTranslationMat(v *Vec3) *Mat4 {
+	return &Mat4{
+		Data: [16]float32{
+			1, 0, 0, v.Data[0],
+			0, 1, 0, v.Data[1],
+			0, 0, 1, v.Data[2],
+			0, 0, 0, 1,
+		},
+	}
+}
+
+func GenScaleMat(v *Vec3) *Mat4 {
+	return &Mat4{
+		Data: [16]float32{
+			v.Data[0], 0, 0, 0,
+			0, v.Data[1], 0, 0,
+			0, 0, v.Data[2], 0,
+			0, 0, 0, 1,
+		},
+	}
+}
+
+func GenRotationMat(q *Quat) *Mat4 {
+
+	xy := q.Data[0] * q.Data[1]
+	xz := q.Data[0] * q.Data[2]
+	xw := q.Data[0] * q.Data[3]
+
+	yz := q.Data[1] * q.Data[2]
+	yw := q.Data[1] * q.Data[3]
+
+	zw := q.Data[2] * q.Data[3]
+
+	xx := q.Data[0] * q.Data[0]
+	yy := q.Data[1] * q.Data[1]
+	zz := q.Data[2] * q.Data[2]
+	ww := q.Data[3] * q.Data[3]
+
+	return &Mat4{
+		Data: [16]float32{
+			xx + yy - zz - ww, 2 * (yz - xw), 2 * (xz + yw), 0,
+			2 * (xw + yz), xx - yy + zz - ww, 2 * (zw - xy), 0,
+			2 * (yw - xz), 2 * (xy + zw), xx - yy - zz + ww, 0,
+			0, 0, 0, 1,
+		},
+	}
+}
+
+// func NewTMatId() *TMat {
+// 	return &TMat{
+// 		Mat4: Mat4{
+// 			Data: [16]float32{
+// 				1, 0, 0, 0,
+// 				0, 1, 0, 0,
+// 				0, 0, 1, 0,
+// 				0, 0, 0, 1,
+// 			},
+// 		},
+// 	}
+// }
diff --git a/gglm/vec3.go b/gglm/vec3.go
index 16cbe8d..f3a27db 100755
--- a/gglm/vec3.go
+++ b/gglm/vec3.go
@@ -105,11 +105,24 @@ func (v *Vec3) Set(x, y, z float32) {
 	v.Data[2] = z
 }
 
-func (v *Vec3) Normalize() {
+//Normalize normalizes this vector and returns it (doesn't copy)
+func (v *Vec3) Normalize() *Vec3 {
 	mag := float32(math.Sqrt(float64(v.X()*v.X() + v.Y()*v.Y() + v.Z()*v.Z())))
 	v.Data[0] /= mag
 	v.Data[1] /= mag
 	v.Data[2] /= mag
+
+	return v
+}
+
+func (v *Vec3) AsRad() *Vec3 {
+	return &Vec3{
+		Data: [3]float32{
+			v.Data[0] * Deg2Rad,
+			v.Data[1] * Deg2Rad,
+			v.Data[2] * Deg2Rad,
+		},
+	}
 }
 
 //AddVec3 v3 = v1 + v2
diff --git a/main.go b/main.go
index 275abe6..3c2125f 100755
--- a/main.go
+++ b/main.go
@@ -118,12 +118,22 @@ func main() {
 	}
 
 	vec3A := gglm.Vec3{Data: [3]float32{1, 2, 3}}
-	lol := gglm.MulMat3Vec3(&mat3A, &vec3A)
-	println(lol.String())
+	mm3v3 := gglm.MulMat3Vec3(&mat3A, &vec3A)
+	println(mm3v3.String())
 
 	//ReflectVec2
 	vec2B := &gglm.Vec2{Data: [2]float32{4, 5}}
 	normA := &gglm.Vec2{Data: [2]float32{0, 1}}
 	rVec2A := gglm.ReflectVec2(vec2B, normA)
 	println(rVec2A.String())
+
+	//Quaternion
+	vRot := &gglm.Vec3{Data: [3]float32{60, 30, 20}}
+	q := gglm.NewQuatEuler(vRot.AsRad())
+	println("\n" + vRot.AsRad().String())
+	println(q.String(), "\n", q.Mag())
+
+	q = gglm.NewQuatAngleAxis(60*gglm.Deg2Rad, vRot.Normalize())
+	println("\n" + vRot.Normalize().String())
+	println(q.String())
 }