diff --git a/Quat.hs b/Quat.hs
deleted file mode 100644
--- a/Quat.hs
+++ /dev/null
@@ -1,64 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-{-# Language StandaloneDeriving #-}
-{-# Language DeriveDataTypeable #-}
-{-# Language DeriveFunctor #-}
-
-module Quat ( Quat(..)
-            , zipWithQuat
-            , inv
-            , norm
-            , normalize
-            , qmult
-            , qmult'
-            ) where
-
-import Data.Data ( Data )
-import Data.Typeable ( Typeable1 )
-
-data Quat a = Quat a a a a deriving (Show, Eq)
-
-deriving instance Typeable1 Quat
-deriving instance Data a => Data (Quat a)
-deriving instance Functor Quat
-
-zipWithQuat :: (a -> b -> c) -> Quat a -> Quat b -> Quat c
-zipWithQuat f (Quat p0 p1 p2 p3) (Quat q0 q1 q2 q3) = Quat (f p0 q0) (f p1 q1) (f p2 q2) (f p3 q3)
-
-instance (Num a, Ord a) => Num (Quat a) where
-  (+) = zipWithQuat (+)
-  (-) = zipWithQuat (-)
-  negate = fmap negate
-  (*) = qmult
-  abs = fmap abs
-  signum = fmap signum
-  fromInteger = error "fromInteger undefined for Quat"
-
--- | q_out = q_in^-1
-inv :: Num a => Quat a -> Quat a
-inv (Quat q0 q1 q2 q3) = Quat q0 (-q1) (-q2) (-q3)
-
--- | return ||q||
-norm :: Floating a => Quat a -> a
-norm (Quat q0 q1 q2 q3) = sqrt $ q0*q0 + q1*q1 + q2*q2 + q3*q3
-
--- | q /= ||q||
-normalize :: Floating a => Quat a -> Quat a
-normalize q = fmap (* normInv) q
-  where
-    normInv = 1/(norm q)
-
--- | quaternion multiply: qa * qb
-qmult :: (Num a, Ord a) => Quat a -> Quat a -> Quat a
-qmult (Quat p0 p1 p2 p3) (Quat q0 q1 q2 q3)
-  | r0 < 0 = negate qOut
-  | otherwise = qOut
-  where
-    qOut = Quat r0 r1 r2 r3
-    r0 = p0*q0 - p1*q1 - p2*q2 - p3*q3
-    r1 = p0*q1 + p1*q0 + p2*q3 - p3*q2
-    r2 = p0*q2 - p1*q3 + p2*q0 + p3*q1
-    r3 = p0*q3 + p1*q2 - p2*q1 + p3*q0
-
--- | quaternion multiply then normalize
-qmult' :: (Floating a, Ord a) => Quat a -> Quat a -> Quat a
-qmult' p q = normalize (qmult q p)
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,10 @@
+spatial-math
+===
+
+[![Build Status](https://secure.travis-ci.org/ghorn/spatial-math.png?branch=master)](http://travis-ci.org/ghorn/spatial-math)
+
+3d math including quaternions/euler angles/dcms and utility functions.
+
+This is a port of my 'mathlib' C library: https://github.com/ghorn/mathlib and may someday be merged with it.
+
+
diff --git a/SpatialMath.hs b/SpatialMath.hs
--- a/SpatialMath.hs
+++ b/SpatialMath.hs
@@ -1,13 +1,12 @@
 {-# OPTIONS_GHC -Wall #-}
 {-# Language StandaloneDeriving #-}
 {-# Language DeriveDataTypeable #-}
+{-# Language DeriveGeneric #-}
 {-# Language DeriveFunctor #-}
+{-# Language DeriveFoldable #-}
+{-# Language DeriveTraversable #-}
 
-module SpatialMath ( module Xyz
-                   , module Quat
-                   , Xyz(..)
-                   , Quat(..)
-                   , Euler(..)
+module SpatialMath ( Euler(..)
                    , rotateXyzAboutX
                    , rotateXyzAboutY
                    , rotateXyzAboutZ
@@ -25,56 +24,105 @@
                    , rotVecByQuatB2A
                    , rotVecByEuler
                    , rotVecByEulerB2A
+                     -- * re-exported from linear
+                   , M33
+                   , V3(..)
+                   , Quaternion(..)
                    ) where
 
-import qualified Xyz
-import qualified Quat
-import Xyz ( Xyz(..) )
-import Quat ( Quat(..) )
-
-import Numeric.LinearAlgebra
-import Foreign.Storable ( Storable )
 import Data.Data ( Data )
+import Data.Foldable ( Foldable )
+import Data.Traversable ( Traversable )
 import Data.Typeable ( Typeable1 )
+import GHC.Generics ( Generic, Generic1 )
+import Linear
 
-data Euler a = Euler a a a deriving (Eq, Show) -- {yaw, pitch, roll}
+normalize' :: Floating a => Quaternion a -> Quaternion a
+normalize' q = fmap (* normInv) q
+  where
+    normInv = 1/(norm q)
 
+--normalize' :: (Floating a, Epsilon a) => Quaternion a -> Quaternion a
+--normalize' = normalize
+
+-- | 3-2-1 Euler angle rotation sequence
+data Euler a = Euler { eYaw :: a
+                     , ePitch :: a
+                     , eRoll :: a
+                     } deriving (Eq, Show, Generic, Generic1, Functor, Foldable, Traversable, Ord)
+
 deriving instance Typeable1 Euler
 deriving instance Data a => Data (Euler a)
-deriving instance Functor Euler
 
-rotateXyzAboutX :: Floating a => Xyz a -> a -> Xyz a
-rotateXyzAboutX (Xyz ax ay az) rotAngle = Xyz bx by bz
+-- | Rotate a vector about the X axis
+--
+-- >>> rotateXyzAboutX (V3 0 1 0) (pi/2)
+-- V3 0.0 6.123233995736766e-17 1.0
+--
+-- >>> rotateXyzAboutX (V3 0 0 1) (pi/2)
+-- V3 0.0 (-1.0) 6.123233995736766e-17
+rotateXyzAboutX :: Floating a => V3 a -> a -> V3 a
+rotateXyzAboutX (V3 ax ay az) rotAngle = V3 bx by bz
   where
     cosTheta = cos rotAngle
     sinTheta = sin rotAngle
 
     bx =  ax
-    by =  ay*cosTheta + az*sinTheta
-    bz = -ay*sinTheta + az*cosTheta
+    by =  ay*cosTheta - az*sinTheta
+    bz =  ay*sinTheta + az*cosTheta
 
-rotateXyzAboutY :: Floating a => Xyz a -> a -> Xyz a
-rotateXyzAboutY (Xyz ax ay az) rotAngle = Xyz bx by bz
+-- | Rotate a vector about the Y axis
+--
+-- >>> rotateXyzAboutY (V3 0 0 1) (pi/2)
+-- V3 1.0 0.0 6.123233995736766e-17
+--
+-- >>> rotateXyzAboutY (V3 1 0 0) (pi/2)
+-- V3 6.123233995736766e-17 0.0 (-1.0)
+rotateXyzAboutY :: Floating a => V3 a -> a -> V3 a
+rotateXyzAboutY (V3 ax ay az) rotAngle = V3 bx by bz
   where
     cosTheta = cos rotAngle
     sinTheta = sin rotAngle
 
-    bx =  ax*cosTheta - az*sinTheta
+    bx =  ax*cosTheta + az*sinTheta
     by =  ay
-    bz =  ax*sinTheta + az*cosTheta
+    bz = -ax*sinTheta + az*cosTheta
 
-rotateXyzAboutZ :: Floating a => Xyz a -> a -> Xyz a
-rotateXyzAboutZ (Xyz ax ay az) rotAngle = Xyz bx by bz
+-- | Rotate a vector about the Z axis
+--
+-- >>> rotateXyzAboutZ (V3 1 0 0) (pi/2)
+-- V3 6.123233995736766e-17 1.0 0.0
+--
+-- >>> rotateXyzAboutZ (V3 0 1 0) (pi/2)
+-- V3 (-1.0) 6.123233995736766e-17 0.0
+--
+rotateXyzAboutZ :: Floating a => V3 a -> a -> V3 a
+rotateXyzAboutZ (V3 ax ay az) rotAngle = V3 bx by bz
   where
     cosTheta = cos rotAngle
     sinTheta = sin rotAngle
 
-    bx =  ax*cosTheta + ay*sinTheta
-    by = -ax*sinTheta + ay*cosTheta
+    bx =  ax*cosTheta - ay*sinTheta
+    by =  ax*sinTheta + ay*cosTheta
     bz =  az
 
-euler321OfQuat :: RealFloat a => Quat a -> Euler a
-euler321OfQuat (Quat q0 q1 q2 q3) = Euler yaw pitch roll
+
+-- | Convert quaternion to Euler angles
+--
+-- >>> euler321OfQuat (Quaternion 1.0 (V3 0.0 0.0 0.0))
+-- Euler {eYaw = 0.0, ePitch = -0.0, eRoll = 0.0}
+--
+-- >>> euler321OfQuat (Quaternion (sqrt(2)/2) (V3 (sqrt(2)/2) 0.0 0.0))
+-- Euler {eYaw = 0.0, ePitch = -0.0, eRoll = 1.5707963267948966}
+--
+-- >>> euler321OfQuat (Quaternion (sqrt(2)/2) (V3 0.0 (sqrt(2)/2) 0.0))
+-- Euler {eYaw = 0.0, ePitch = 1.5707963267948966, eRoll = 0.0}
+--
+-- >>> euler321OfQuat (Quaternion (sqrt(2)/2) (V3 0.0 0.0 (sqrt(2)/2)))
+-- Euler {eYaw = 1.5707963267948966, ePitch = -0.0, eRoll = 0.0}
+--
+euler321OfQuat :: RealFloat a => Quaternion a -> Euler a
+euler321OfQuat (Quaternion q0 (V3 q1 q2 q3)) = Euler yaw pitch roll
   where
     r11 = q0*q0 + q1*q1 - q2*q2 - q3*q3
     r12 = 2.0*(q1*q2 + q0*q3)
@@ -90,31 +138,67 @@
     pitch = asin mr13
     roll  = atan2 r23 r33
 
-quatOfDcm :: (Storable a, RealFloat a) => Matrix a -> Quat a
+-- | convert a DCM to a quaternion
+--
+-- >>> quatOfDcm $ V3 (V3 1 0 0) (V3 0 1 0) (V3 0 0 1)
+-- Quaternion 1.0 (V3 0.0 0.0 0.0)
+--
+-- >>> quatOfDcm $ V3 (V3 0 1 0) (V3 (-1) 0 0) (V3 0 0 1)
+-- Quaternion 0.7071067811865476 (V3 0.0 0.0 0.7071067811865475)
+--
+-- >>> let s = sqrt(2)/2 in quatOfDcm $ V3 (V3 s s 0) (V3 (-s) s 0) (V3 0 0 1)
+-- Quaternion 0.9238795325112867 (V3 0.0 0.0 0.3826834323650898)
+--
+quatOfDcm :: RealFloat a => M33 a -> Quaternion a
 quatOfDcm = quatOfEuler321 . euler321OfDcm
 
-quatOfDcmB2A :: (Storable a, RealFloat a) => Matrix a -> Quat a
-quatOfDcmB2A = Quat.inv . quatOfDcm
+quatOfDcmB2A :: (Conjugate a, RealFloat a) => M33 a -> Quaternion a
+quatOfDcmB2A = conjugate . quatOfDcm
 
-euler321OfDcm :: (RealFloat a, Storable a) => Matrix a -> Euler a
-euler321OfDcm r = Euler yaw pitch roll
+-- | Convert DCM to euler angles
+--
+-- >>> euler321OfDcm $ V3 (V3 1 0 0) (V3 0 1 0) (V3 0 0 1)
+-- Euler {eYaw = 0.0, ePitch = -0.0, eRoll = 0.0}
+--
+-- >>> euler321OfDcm $ V3 (V3 0 1 0) (V3 (-1) 0 0) (V3 0 0 1)
+-- Euler {eYaw = 1.5707963267948966, ePitch = -0.0, eRoll = 0.0}
+--
+-- >>> let s = sqrt(2)/2 in euler321OfDcm $ V3 (V3 s s 0) (V3 (-s) s 0) (V3 0 0 1)
+-- Euler {eYaw = 0.7853981633974483, ePitch = -0.0, eRoll = 0.0}
+--
+euler321OfDcm :: RealFloat a => M33 a -> Euler a
+euler321OfDcm
+  (V3
+   (V3 r11 r12 r13)
+   (V3   _   _ r23)
+   (V3   _   _ r33)) = Euler yaw pitch roll
   where
-    r11 = r @@> (0,0)
-    r12 = r @@> (0,1)
-    mr13' = -(r @@> (0,2))
+    mr13' = -r13
     mr13 -- nan protect
       | mr13' >  1 =  1
       | mr13' < -1 = -1
       | otherwise = mr13'
-    r23 = r @@> (1,2)
-    r33 = r @@> (2,2)
-  
+
     yaw   = atan2 r12 r11
     pitch = asin mr13
     roll  = atan2 r23 r33
 
-quatOfEuler321 :: (Floating a, Ord a) => Euler a -> Quat a
-quatOfEuler321 (Euler yaw pitch roll) = Quat.normalize q
+-- | Convert Euler angles to quaternion
+--
+-- >>> quatOfEuler321 (Euler 0 0 0)
+-- Quaternion 1.0 (V3 0.0 0.0 0.0)
+--
+-- >>> quatOfEuler321 (Euler (pi/2) 0 0)
+-- Quaternion 0.7071067811865476 (V3 0.0 0.0 0.7071067811865475)
+--
+-- >>> quatOfEuler321 (Euler 0 (pi/2) 0)
+-- Quaternion 0.7071067811865476 (V3 0.0 0.7071067811865475 0.0)
+--
+-- >>> quatOfEuler321 (Euler 0 0 (pi/2))
+-- Quaternion 0.7071067811865476 (V3 0.7071067811865475 0.0 0.0)
+--
+quatOfEuler321 :: (Floating a, Ord a) => Euler a -> Quaternion a
+quatOfEuler321 (Euler yaw pitch roll) = normalize' q
   where
     sr2 = sin $ 0.5*roll
     cr2 = cos $ 0.5*roll
@@ -127,57 +211,78 @@
     q2 = cr2*sp2*cy2 + sr2*cp2*sy2
     q3 = cr2*cp2*sy2 - sr2*sp2*cy2
 
-    q' = Quat q0 q1 q2 q3
-    
+    q' = Quaternion q0 (V3 q1 q2 q3)
+
     q
-      | q0 < 0 = negate q'
+      | q0 < 0 = Quaternion (-q0) (V3 (-q1) (-q2) (-q3))
       | otherwise = q'
 
-dcmOfQuat :: (Num a, Element a) => Quat a -> Matrix a
-dcmOfQuat (Quat q0 q1 q2 q3) = fromLists [ [r0, r1, r2]
-                                         , [r3, r4, r5]
-                                         , [r6, r7, r8]
-                                         ]
+-- | convert a quaternion to a DCM
+--
+-- >>> dcmOfQuat $ Quaternion 1.0 (V3 0.0 0.0 0.0)
+-- V3 (V3 1.0 0.0 0.0) (V3 0.0 1.0 0.0) (V3 0.0 0.0 1.0)
+--
+-- >>> let s = sqrt(2)/2 in dcmOfQuat $ Quaternion s (V3 0.0 0.0 s)
+-- V3 (V3 0.0 1.0000000000000002 0.0) (V3 (-1.0000000000000002) 0.0 0.0) (V3 0.0 0.0 1.0000000000000002)
+--
+-- >>> dcmOfQuat $ Quaternion 0.9238795325112867 (V3 0.0 0.0 0.3826834323650898)
+-- V3 (V3 0.7071067811865475 0.7071067811865476 0.0) (V3 (-0.7071067811865476) 0.7071067811865475 0.0) (V3 0.0 0.0 1.0)
+--
+dcmOfQuat :: Num a => Quaternion a -> M33 a
+dcmOfQuat (Quaternion q0 (V3 q1 q2 q3)) = V3 (V3 r0 r1 r2)
+                                             (V3 r3 r4 r5)
+                                             (V3 r6 r7 r8)
   where
     -- 1st column
     r0 = q0*q0 + q1*q1 - q2*q2 - q3*q3
     r3 = 2*(q1*q2 - q0*q3)
     r6 = 2*(q1*q3 + q0*q2)
-  
+
     -- 2nd column
     r1 = 2*(q1*q2 + q0*q3)
     r4 = q0*q0 - q1*q1 + q2*q2 - q3*q3
     r7 = 2*(q2*q3 - q0*q1)
-  
+
     -- 3rd column
     r2 = 2*(q1*q3 - q0*q2)
     r5 = 2*(q2*q3 + q0*q1)
     r8 = q0*q0 - q1*q1 - q2*q2 + q3*q3
 
-dcmOfEuler321 :: (Floating a, Element a, Ord a) => Euler a -> Matrix a
+-- | Convert DCM to euler angles
+--
+-- >>> dcmOfEuler321 $ Euler {eYaw = 0.0, ePitch = 0, eRoll = 0}
+-- V3 (V3 1.0 0.0 0.0) (V3 0.0 1.0 0.0) (V3 0.0 0.0 1.0)
+--
+-- >>> dcmOfEuler321 $ Euler {eYaw = pi/2, ePitch = 0, eRoll = 0}
+-- V3 (V3 2.220446049250313e-16 1.0 0.0) (V3 (-1.0) 2.220446049250313e-16 0.0) (V3 0.0 0.0 1.0)
+--
+-- >>> dcmOfEuler321 $ Euler {eYaw = pi/4, ePitch = 0, eRoll = 0}
+-- V3 (V3 0.7071067811865475 0.7071067811865476 0.0) (V3 (-0.7071067811865476) 0.7071067811865475 0.0) (V3 0.0 0.0 1.0)
+--
+dcmOfEuler321 :: (Floating a, Ord a) => Euler a -> M33 a
 dcmOfEuler321 = dcmOfQuat . quatOfEuler321
 
-dcmOfQuatB2A :: (Num a, Element a) => Quat a -> Matrix a
-dcmOfQuatB2A = dcmOfQuat . Quat.inv
+dcmOfQuatB2A :: (Conjugate a, RealFloat a) => Quaternion a -> M33 a
+dcmOfQuatB2A = dcmOfQuat . conjugate
 
 -- | vec_b = R_a2b * vec_a
-rotVecByDcm :: (Num a, Storable a) => Matrix a -> Xyz a -> Xyz a
-rotVecByDcm dcm vec = Xyz.mult3x3ByXyz dcm vec
+rotVecByDcm :: Num a => M33 a -> V3 a -> V3 a
+rotVecByDcm dcm vec = dcm !* vec
 
 -- | vec_a = R_a2b^T * vec_b
-rotVecByDcmB2A :: (Num a, Storable a) => Matrix a -> Xyz a -> Xyz a
-rotVecByDcmB2A dcm vec = Xyz.mult3x3TransposeByXyz dcm vec
+rotVecByDcmB2A :: Num a => M33 a -> V3 a -> V3 a
+rotVecByDcmB2A dcm vec = vec *! dcm
 
 -- | vec_b = q_a2b * vec_a * q_a2b^(-1)
 --   vec_b = R(q_a2b) * vec_a
-rotVecByQuat :: (Num a, Element a) => Quat a -> Xyz a -> Xyz a
+rotVecByQuat :: Num a => Quaternion a -> V3 a -> V3 a
 rotVecByQuat q = rotVecByDcm (dcmOfQuat q)
 
-rotVecByQuatB2A :: (Num a, Element a) => Quat a -> Xyz a -> Xyz a
+rotVecByQuatB2A :: Num a => Quaternion a -> V3 a -> V3 a
 rotVecByQuatB2A q = rotVecByDcmB2A (dcmOfQuat q)
 
-rotVecByEuler :: (Floating a, Element a, Ord a) => Euler a -> Xyz a -> Xyz a
+rotVecByEuler :: (Floating a, Ord a) => Euler a -> V3 a -> V3 a
 rotVecByEuler = rotVecByDcm . dcmOfEuler321
 
-rotVecByEulerB2A :: (Floating a, Element a, Ord a) => Euler a -> Xyz a -> Xyz a
+rotVecByEulerB2A :: (Floating a, Ord a) => Euler a -> V3 a -> V3 a
 rotVecByEulerB2A = rotVecByDcmB2A . dcmOfEuler321
diff --git a/Xyz.hs b/Xyz.hs
deleted file mode 100644
--- a/Xyz.hs
+++ /dev/null
@@ -1,128 +0,0 @@
-{-# OPTIONS_GHC -Wall #-}
-{-# Language StandaloneDeriving #-}
-{-# Language DeriveDataTypeable #-}
-
-module Xyz ( Xyz(..)
-           , zipWithXyz
-           , cross
-           , dot
-           , normSquared
-           , norm
-           , distance
-           , scale
-           , normalizeTo
-           , normalize
-           , mult3x3ByXyz
-           , mult3x3TransposeByXyz
-           ) where
-
-import Numeric.LinearAlgebra ( (@@>), Matrix )
-import Foreign.Storable ( Storable )
-import Data.Data ( Data )
-import Data.Typeable ( Typeable1 )
-import System.Random ( Random(..) )
-
-data Xyz a = Xyz a a a deriving (Show, Eq)
-
-deriving instance Typeable1 Xyz
-deriving instance Data a => Data (Xyz a)
-
-instance Functor Xyz where
-  fmap f (Xyz x y z) = Xyz (f x) (f y) (f z)
-
-instance Random a => Random (Xyz a) where
-  random g0 = (Xyz x y z, gz)
-    where
-      (x,gx) = random g0
-      (y,gy) = random gx
-      (z,gz) = random gy
-  randomR (Xyz x0 y0 z0, Xyz x1 y1 z1) g0 = (Xyz x y z, gz)
-    where
-      (x,gx) = randomR (x0,x1) g0
-      (y,gy) = randomR (y0,y1) gx
-      (z,gz) = randomR (z0,z1) gy
-
-zipWithXyz :: (a -> b -> c) -> Xyz a -> Xyz b -> Xyz c
-zipWithXyz f (Xyz x0 y0 z0) (Xyz x1 y1 z1) = Xyz (f x0 x1) (f y0 y1) (f z0 z1)
-
-instance (Num a) => Num (Xyz a) where
-  (+) = zipWithXyz (+)
-  (-) = zipWithXyz (-)
-  negate = fmap negate
-  (*) = zipWithXyz (*)
-  abs = fmap abs
-  signum = fmap signum
-  fromInteger k = fmap fromInteger (Xyz k k k)
-
-instance (Fractional a) => Fractional (Xyz a) where
-  fromRational r = fmap fromRational (Xyz r r r)
-  (/) = zipWithXyz (/)
-
-instance (Floating a) => Floating (Xyz a) where
-  pi = Xyz pi pi pi
-  exp   = fmap exp
-  log   = fmap log
-  sin   = fmap sin
-  cos   = fmap cos
-  tan   = fmap tan
-  asin  = fmap asin
-  acos  = fmap acos
-  atan  = fmap atan
-  sinh  = fmap sinh
-  cosh  = fmap cosh
-  tanh  = fmap tanh
-  asinh = fmap asinh
-  acosh = fmap acosh
-  atanh = fmap atanh
-
--- | c = a (cross) b
-cross :: Num a => Xyz a -> Xyz a -> Xyz a
-cross (Xyz ax ay az) (Xyz bx by bz) = Xyz cx cy cz
-  where
-    cx =   ay*bz - az*by
-    cy = - ax*bz + az*bx
-    cz =   ax*by - ay*bx
-
--- | c = a (dot) b
-dot :: Num a => Xyz a -> Xyz a -> a
-dot (Xyz ax ay az) (Xyz bx by bz) = ax*bx + ay*by + az*bz;
-
--- | c = vec (dot) vec
-normSquared :: Num a => Xyz a -> a
-normSquared x = dot x x
-
--- | norm(x)
-norm :: Floating a => Xyz a -> a
-norm x = sqrt $ dot x x
-
--- | norm(a - b)
-distance :: Floating a => Xyz a -> Xyz a -> a
-distance a b = norm $ a - b
-
--- | vec_out = vec_in*scale_factor
-scale :: Num a => a -> Xyz a -> Xyz a
-scale k = fmap (k *)
-
--- | vec_out = scale (new_norm/norm(vec_in)) vec_in
-normalizeTo :: Floating a => a -> Xyz a -> Xyz a -> Xyz a
-normalizeTo newNorm vec = scale (newNorm/(norm(vec) + 1e-12))
-
--- | vec_out = vec_in/norm(vec_in)
-normalize :: Floating a => Xyz a -> Xyz a -> Xyz a
-normalize = normalizeTo 1
-
--- | v_out = M*v
-mult3x3ByXyz :: (Num a, Storable a) => Matrix a -> Xyz a -> Xyz a
-mult3x3ByXyz mat (Xyz x y z) = Xyz x' y' z'
-  where
-    x' = (mat @@> (0,0))*x + (mat @@> (0,1))*y +  (mat @@> (0,2))*z
-    y' = (mat @@> (1,0))*x + (mat @@> (1,1))*y +  (mat @@> (1,2))*z
-    z' = (mat @@> (2,0))*x + (mat @@> (2,1))*y +  (mat @@> (2,2))*z
-
--- // v_out = M^T*v
-mult3x3TransposeByXyz :: (Num a, Storable a) => Matrix a -> Xyz a -> Xyz a
-mult3x3TransposeByXyz mat (Xyz x y z) = Xyz x' y' z'
-  where
-    x' = (mat @@> (0,0))*x + (mat @@> (1,0))*y +  (mat @@> (2,0))*z
-    y' = (mat @@> (0,1))*x + (mat @@> (1,1))*y +  (mat @@> (2,1))*z
-    z' = (mat @@> (0,2))*x + (mat @@> (1,2))*y +  (mat @@> (2,2))*z
diff --git a/changelog.txt b/changelog.txt
new file mode 100644
--- /dev/null
+++ b/changelog.txt
@@ -0,0 +1,4 @@
+0.2.0
+- convert to using `linear` V3, M33, Quaternion types
+- doctests
+- fix long unknown bug: rotateXyzAbout{X,Y,Z} was rotating opposite direction from intended
diff --git a/spatial-math.cabal b/spatial-math.cabal
--- a/spatial-math.cabal
+++ b/spatial-math.cabal
@@ -1,25 +1,34 @@
 name:                spatial-math
-version:             0.1.7
+version:             0.2.0
 synopsis:            3d math including quaternions/euler angles/dcms and utility functions
-description:         This is a port of my 'mathlib' C library: https://github.com/ghorn/mathlib
+description:         This is a port of my 'mathlib' C library: `https://github.com/ghorn/mathlib`
 license:             BSD3
 license-file:        LICENSE
 author:              Greg Horn
 maintainer:          gregmainland@gmail.com
--- copyright:           
+copyright:           Copyright (c) 2012, Greg Horn
 category:            Math
 build-type:          Simple
-cabal-version:       >=1.8
+cabal-version:       >=1.10
 
+extra-source-files:  README.md
+                     changelog.txt
+
 library
-  exposed-modules:     SpatialMath,
-                       Quat,
-                       Xyz
-  -- other-modules:       
+  exposed-modules:     SpatialMath
   build-depends:       base >= 4 && < 5,
-                       hmatrix >= 0.14 && < 0.15,
-                       random
+                       linear >= 1.3.1
+  default-language:    Haskell2010
 
 source-repository head
   type:     git
   location: git://github.com/ghorn/spatial-math.git
+
+test-suite doctests
+  type:                exitcode-stdio-1.0
+  main-is:             doctests.hs
+  build-depends:       base >= 4 && < 5,
+                       doctest >= 0.8
+  default-language:    Haskell2010
+  ghc-options:         -threaded
+  hs-source-dirs:      tests
diff --git a/tests/doctests.hs b/tests/doctests.hs
new file mode 100644
--- /dev/null
+++ b/tests/doctests.hs
@@ -0,0 +1,8 @@
+{-# OPTIONS_GHC -Wall #-}
+
+module Main where
+
+import Test.DocTest
+
+main :: IO ()
+main = doctest ["SpatialMath.hs"]