{-

 Haskell Quaternion Module
 Author : Uğur Yörük , 040060244
 
 -}

module Number.Quaternion (
  Quaternion(..),
  Others(..),
  addQuats,
  mulQuats,
  subQuats,
  getLength,
  multiquat,
  divQuats,
  defQfromVec,
  neg,
  getEuclidianLength,
  conjugate,
  absQ,
  normalize,
  toAngleAxis,
  toMat3,
  toMat4
) where

data Quaternion = Q { w :: Float
                     ,x :: Float
                     ,y :: Float
                     ,z :: Float } deriving (Eq,Show,Read,Ord)

-- Simple data structures used by Quaternion
data Others = Vec3 Float Float Float | 
	      Mat3 Float Float Float Float Float Float Float Float Float | 
	      Mat4 Float Float Float Float Float Float Float Float Float Float Float Float Float Float Float Float 
		deriving (Show)

-- Define Quaternion from Vec3
defQfromVec :: Float -> Others -> Quaternion
defQfromVec a (Vec3 x y z) = Q a x y z

addQuats :: Quaternion -> Quaternion -> Quaternion
addQuats (Q a1 b1 c1 d1) (Q a2 b2 c2 d2) = Q (a1+a2) (b1+b2) (c1+c2) (d1+d2)

mulQuats :: Quaternion -> Quaternion -> Quaternion
mulQuats (Q a1 b1 c1 d1) (Q a2 b2 c2 d2) = Q (a1*a2) (b1*b2) (c1*c2) (d1*d2)

subQuats :: Quaternion -> Quaternion -> Quaternion
subQuats (Q a1 b1 c1 d1) (Q a2 b2 c2 d2) = Q (a1-a2) (b1-b2) (c1-c2) (d1-d2)

divQuats :: Quaternion -> Quaternion -> Quaternion
divQuats (Q a1 b1 c1 d1) (Q a2 b2 c2 d2) = Q (a1/a2) (b1/b2) (c1/c2) (d1/d2)

neg :: Quaternion -> Quaternion
neg (Q a b c d) = Q (-a) (-b) (-c) (-d)

absQ :: Quaternion -> Float
absQ (Q a b c d)
  | a < 0 = absQ (Q (-a) (b) (c) (d))
  | b < 0 = absQ (Q (a) (-b) (c) (d))
  | c < 0 = absQ (Q (a) (b) (-c) (d))
  | d < 0 = absQ (Q (a) (b) (c) (-d))
  | otherwise = sqrt (a*a + b*b + c*c + d*d)

instance Num Quaternion where
  (+) = addQuats
  (*) = mulQuats
  (-) = subQuats
  negate = neg
  abs = undefined
  fromInteger = undefined
  signum = undefined

instance Fractional Quaternion where
  (/) = divQuats
  fromRational = undefined
  

getEuclidianLength :: Quaternion -> Float
getEuclidianLength (Q a b c d) = sqrt (a*a + b*b + c*c + d*d)

getLength :: Quaternion -> Float
getLength (Q a b c d ) = (a*a + b*b + c*c + d*d)

multiquat :: Quaternion -> Float -> Quaternion
multiquat (Q a b c d) e = Q (a*e) (b*e) (c*e) (d*e)

conjugate :: Quaternion -> Quaternion
conjugate (Q a b c d) = Q (a) (-b) (-c) (-d)

normalize :: Quaternion -> Quaternion
normalize (Q a b c d) = Q (a*l) (b*l) (c*l) (d*l)
  where l = 1.0 / absQ (Q a b c d)


toAngleAxis :: Quaternion -> (Float,Others)
toAngleAxis (Q a b c d) = ((2.0*a),(Vec3 x y z))
  where Q a1 b1 c1 d1 = normalize (Q a b c d)
	w = max (min a1 1.0) (-1.0)
	w1 = acos w
	s = sin w1
	x = b1 / s
	y = c1 / s
	z = d1 / s

toMat3 :: Quaternion -> Others
toMat3 (Q a b c d) = Mat3 (1.0-bb-cc) (ab-cd) (ac+bd) 
			  (ab+cd) (1.0-aa-cc) (bc-ad) 
			  (ac-bd) (bc+ad) (1.0-aa-bb)
  where aa = a*a*2.0
	bb = b*b*2.0
	cc = c*c*2.0
	ab = a*b*2.0
	cd = c*d*2.0
	ac = a*c*2.0
	bd = b*d*2.0
	ad = a*d*2.0
	bc = b*c*2.0

toMat4 :: Quaternion -> Others
toMat4 (Q a b c d) = Mat4 (1.0-bb-cc) (ab-cd) (ac+bd) (0.0) 
			  (ab+cd) (1.0-aa-cc) (bc-ad) (0.0) 
			  (ac-bd) (bc+ad) (1.0-aa-bb) (0.0) 
			  (0.0) (0.0) (0.0) (1.0)
  where aa = a*a*2.0
	bb = b*b*2.0
	cc = c*c*2.0
	ab = a*b*2.0
	cd = c*d*2.0
	ac = a*c*2.0
	bd = b*d*2.0
	ad = a*d*2.0
	bc = b*c*2.0
