{-
   Haskell Dynamics Engine, Copyright (C) 2007 Ruben Henner Zilibowitz
   All rights reserved. Email: rubenz@cse.unsw.edu.au Web: www.cse.unsw.edu.au/~rubenz
   
   This is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License that is include with this
   package in the file LICENSE-GPL.TXT.
-}

-- Quaternions
-- Code by Ruben Henner Zilibowitz
-- Created on 20/3/07

module Quaternions where

data (RealFloat a) => Quaternion a = Q !a !a !a !a   deriving (Eq, Read, Show)

instance RealFloat a => Num (Quaternion a) where
   fromInteger a = Q (fromInteger a) 0 0 0
   (Q a b c d) + (Q a' b' c' d') = Q (a+a') (b+b')
                                     (c+c') (d+d')
   (Q a b c d) - (Q a' b' c' d') = Q (a-a') (b-b')
                                     (c-c') (d-d')
   (Q a b c d) * (Q a' b' c' d') = Q (a*a'-b*b'-c*c'-d*d') (a*b'+b*a'+c*d'-d*c') (a*c'+c*a'+d*b'-b*d') (a*d'+d*a'+b*c'-c*b')
   abs q = Q (magnitude q) 0 0 0
   signum (Q 0 0 0 0) = 0
   signum q@(Q a b c d) = let r = magnitude q in Q (a/r) (b/r) (c/r) (d/r)

instance  (RealFloat a) => Fractional (Quaternion a)  where
   q / r = q * (inverse r)
   fromRational a = Q (fromRational a) 0 0 0

magnitude (Q a b c d) = sqrt (a*a + b*b + c*c + d*d)

conjugate (Q a b c d) = Q a (-b) (-c) (-d)

inverse (Q a b c d) = Q (a*r) (-b*r) (-c*r) (-d*r)   where r = recip (a*a + b*b + c*c + d*d)

normalise q@(Q a b c d) = Q (a*r) (b*r) (c*r) (d*r)   where r = recip (magnitude q)