module Vectors where

infix 7 `scale`
infix 7 `mulVec`

infinite :: Scal 
infinite = 1.0e20

type Scal = Float --Int

data Vector 
    = Vector 
      {-# UNPACK #-}!Scal 
      {-# UNPACK #-}!Scal
      {-# UNPACK #-}!Scal deriving Eq

instance Show Vector where
  show (Vector x y z) = show (x,y,z)

data Matrix 
    = Matrix 
      {-# UNPACK #-}!Scal 
      {-# UNPACK #-}!Scal
      {-# UNPACK #-}!Scal 
      {-# UNPACK #-}!Scal
      {-# UNPACK #-}!Scal 
      {-# UNPACK #-}!Scal
      {-# UNPACK #-}!Scal 
      {-# UNPACK #-}!Scal
      {-# UNPACK #-}!Scal

norm = sqrt . sqnorm

sqnorm (Vector x y z) = x*x+y*y+z*z

normalized v = scale (1.0 / norm v) v

scale f (Vector x y z) = Vector (x*f) (y*f) (z*f)


cross (Vector vx vy vz) (Vector ux uy uz) = 
    Vector
    (vy * uz - vz * uy)
    (vz * ux - vx * uz)
    (vx * uy - vy * ux)

constVec c  = Vector c c c
liftV op (Vector vx vy vz) = Vector (op vx) (op vy) (op vz)
liftV2 op (Vector vx vy vz) (Vector ux uy uz) = Vector (vx`op`ux) (vy`op`uy) (vz`op`uz)

dotProd (Vector vx vy vz) (Vector ux uy uz) =
    vx*ux + vy*uy + vz*uz

(.*) = dotProd

rowProd = liftV2 (*)

instance Num Vector where
    (+) = liftV2 (+)
    (-) = liftV2 (-)
    (*) = cross
    fromInteger i = constVec $ fromInteger i
    signum = liftV signum
    negate = liftV negate
    abs = liftV abs

vec = Vector

rotateMatrix theta axis@(Vector a b c) = Matrix
-- rotates "th" radians around axis "ax"
-- !!! axis must be normalized !!!
					 
   (a*a*oct+ct)   (a*b*oct-c*st) (a*c*oct+b*st)
   (b*a*oct+c*st) (b*b*oct+ct)   (b*c*oct-a*st)
   (c*a*oct-b*st) (c*b*oct+a*st) (c*c*oct+ct)
      where ct = cos theta
 	    st = sin theta
 	    oct = 1.0 - ct


transpose (Matrix
	   m11 m12 m13
	   m21 m22 m23
	   m31 m32 m33) = 
	  (Matrix
	   m11 m21 m31
	   m21 m22 m32
	   m13 m23 m33)


mulVec (Matrix
	   m11 m12 m13
	   m21 m22 m23
	   m31 m32 m33) (Vector x y z) =
  Vector 
  (m11*x + m12*y + m13*z)
  (m21*x + m22*y + m23*z)
  (m31*x + m32*y + m33*z)