content: -- Copyright (c) 2000 Galois Connections, Inc.
  -- All rights reserved.  This software is distributed as
  -- free software under the license in the file "LICENSE",
  -- which is included in the distribution.
  
  module Geometry
      ( Coords
      , Ray
      , Point  -- abstract
      , Vector -- abstract
      , Matrix -- abstract
      , Color  -- abstract
      , Box(..)
      , Radian
      , matrix
      , coord
      , color
      , uncolor
      , xCoord , yCoord , zCoord
      , xComponent , yComponent , zComponent
      , point
      , vector
      , nearV
      , point_to_vector
      , vector_to_point
      , dot
      , cross
      , tangents
      , addVV
      , addPV
      , subVV
      , negV
      , subPP
      , norm
      , normalize
      , dist2
      , sq
      , distFrom0Sq
      , distFrom0
      , multSV
      , multMM
      , transposeM
      , multMV
      , multMP
      , multMQ
      , multMR
      , white
      , black
      , addCC
      , subCC
      , sumCC
      , multCC
      , multSC
      , nearC
      , offsetToPoint
      , epsilon
      , inf
      , nonZero
      , eqEps
      , near
      , clampf
      ) where
  
  
  type Coords = (Double,Double,Double)
  
  type Ray = (Point,Vector)    -- origin of ray, and unit vector giving direction
  
  data Point  = P !Double !Double !Double -- implicit extra arg of 1
      deriving (Show)
  data Vector = V !Double !Double !Double -- implicit extra arg of 0
      deriving (Show, Eq)
  data Matrix = M !Quad   !Quad   !Quad   !Quad
      deriving (Show)
  
  data Color  = C !Double !Double !Double
      deriving (Show, Eq)
  
  data Box = B !Double !Double !Double !Double !Double !Double
      deriving (Show)
  
  data Quad   = Q !Double !Double !Double !Double
      deriving (Show)
  
  type Radian = Double
  
  type Tup4 a = (a,a,a,a)
  
  --{-# INLINE matrix #-}
  matrix :: Tup4 (Tup4 Double) -> Matrix
  matrix ((m11, m12, m13, m14),
            (m21, m22, m23, m24),
            (m31, m32, m33, m34),
            (m41, m42, m43, m44))
    = M (Q m11 m12 m13 m14)
        (Q m21 m22 m23 m24)
        (Q m31 m32 m33 m34)
        (Q m41 m42 m43 m44)
  
  coord x y z = (x, y, z)
  
  color r g b = C r g b
  
  uncolor (C r g b) = (r,g,b)
  
  {-# INLINE xCoord #-}
  xCoord (P x y z) = x
  {-# INLINE yCoord #-}
  yCoord (P x y z) = y
  {-# INLINE zCoord #-}
  zCoord (P x y z) = z
  
  {-# INLINE xComponent #-}
  xComponent (V x y z) = x
  {-# INLINE yComponent #-}
  yComponent (V x y z) = y
  {-# INLINE zComponent #-}
  zComponent (V x y z) = z
  
  point :: Double -> Double -> Double -> Point
  point x y z = P x y z
  
  vector :: Double -> Double -> Double -> Vector
  vector x y z = V x y z
  
  nearV :: Vector -> Vector -> Bool
  nearV (V a b c) (V d e f) = a `near` d && b `near` e && c `near` f
  
  point_to_vector :: Point -> Vector
  point_to_vector (P x y z) = V x y z
  
  vector_to_point :: Vector -> Point
  vector_to_point (V x y z)  = P x y z
  
  {-# INLINE vector_to_quad #-}
  vector_to_quad :: Vector -> Quad
  vector_to_quad (V x y z) = Q x y z 0
  
  {-# INLINE point_to_quad #-}
  point_to_quad :: Point -> Quad
  point_to_quad (P x y z) = Q x y z 1
  
  {-# INLINE quad_to_point #-}
  quad_to_point :: Quad -> Point
  quad_to_point (Q x y z _) = P x y z
  
  {-# INLINE quad_to_vector #-}
  quad_to_vector :: Quad -> Vector
  quad_to_vector (Q x y z _) = V x y z
  
  --{-# INLINE dot #-}
  dot :: Vector -> Vector -> Double
  dot (V x1 y1 z1) (V x2 y2 z2) = x1 * x2 + y1 * y2 + z1 * z2
  
  cross :: Vector -> Vector -> Vector
  cross (V x1 y1 z1) (V x2 y2 z2)
    = V (y1 * z2 - z1 * y2) (z1 * x2 - x1 * z2) (x1 * y2 - y1 * x2)
  
  -- assumption: the input vector is a normal
  tangents :: Vector -> (Vector, Vector)
  tangents v@(V x y z)
    = (v1, v `cross` v1)
    where v1 | x == 0    = normalize (vector 0 z (-y))
  	   | otherwise = normalize (vector (-y) x 0)
  
  {-# INLINE dot4 #-}
  dot4 :: Quad -> Quad -> Double
  dot4 (Q x1 y1 z1 w1) (Q x2 y2 z2 w2) = x1 * x2 + y1 * y2 + z1 * z2 + w1 * w2
  
  addVV :: Vector -> Vector -> Vector
  addVV (V x1 y1 z1) (V x2 y2 z2)
      = V (x1 + x2) (y1 + y2) (z1 + z2)
  
  addPV :: Point -> Vector -> Point
  addPV (P x1 y1 z1) (V x2 y2 z2)
      = P (x1 + x2) (y1 + y2) (z1 + z2)
  
  subVV :: Vector -> Vector -> Vector
  subVV (V x1 y1 z1) (V x2 y2 z2)
      = V (x1 - x2) (y1 - y2) (z1 - z2)
  
  negV :: Vector -> Vector
  negV (V x1 y1 z1)
      = V (-x1) (-y1) (-z1)
  
  subPP :: Point -> Point -> Vector
  subPP (P x1 y1 z1) (P x2 y2 z2)
      = V (x1 - x2) (y1 - y2) (z1 - z2)
  
  --{-# INLINE norm #-}
  norm :: Vector -> Double
  norm (V x y z) = sqrt (sq x + sq y + sq z)
  
  --{-# INLINE normalize #-}
  -- normalize a vector to a unit vector
  normalize :: Vector -> Vector
  normalize v@(V x y z)
               | norm /= 0 = multSV (1/norm) v
  	     | otherwise = error "normalize empty!"
      where norm = sqrt (sq x + sq y + sq z)
  
  -- This does computes the distance *squared*
  dist2 :: Point -> Point -> Double
  dist2 us vs = sq x + sq y + sq z
      where
         (V x y z) = subPP us vs
  
  {-# INLINE sq #-}
  sq :: Double -> Double
  sq d = d * d
  
  {-# INLINE distFrom0Sq #-}
  distFrom0Sq :: Point -> Double  -- Distance of point from origin.
  distFrom0Sq (P x y z) = sq x + sq y + sq z
  
  {-# INLINE distFrom0 #-}
  distFrom0 :: Point -> Double  -- Distance of point from origin.
  distFrom0 p = sqrt (distFrom0Sq p)
  
  --{-# INLINE multSV #-}
  multSV :: Double -> Vector -> Vector
  multSV k (V x y z) = V (k*x) (k*y) (k*z)
  
  --{-# INLINE multMM #-}
  multMM :: Matrix -> Matrix -> Matrix
  multMM m1@(M q1 q2 q3 q4) m2
       = M (multMQ m2' q1)
           (multMQ m2' q2)
           (multMQ m2' q3)
           (multMQ m2' q4)
    where
       m2' = transposeM m2
  
  {-# INLINE transposeM #-}
  transposeM :: Matrix -> Matrix
  transposeM (M (Q e11  e12  e13  e14)
                (Q e21  e22  e23  e24)
                (Q e31  e32  e33  e34)
                (Q e41  e42  e43  e44)) = (M (Q e11  e21  e31  e41)
                                             (Q e12  e22  e32  e42)
                                             (Q e13  e23  e33  e43)
                                             (Q e14  e24  e34  e44))
  
  
  --multMM m1 m2 = [map (dot4 row) (transpose m2) | row <- m1]
  
  --{-# INLINE multMV #-}
  multMV :: Matrix -> Vector -> Vector
  multMV m v = quad_to_vector (multMQ m (vector_to_quad v))
  
  --{-# INLINE multMP #-}
  multMP :: Matrix -> Point -> Point
  multMP m p = quad_to_point (multMQ m (point_to_quad p))
  
  -- mat vec = map (dot4 vec) mat
  
  {-# INLINE multMQ #-}
  multMQ :: Matrix -> Quad -> Quad
  multMQ (M q1 q2 q3 q4) q
         = Q (dot4 q q1)
             (dot4 q q2)
             (dot4 q q3)
             (dot4 q q4)
  
  {-# INLINE multMR #-}
  multMR :: Matrix -> Ray -> Ray
  multMR m (r, v) = (multMP m r, multMV m v)
  
  white :: Color
  white = C 1 1 1
  black :: Color
  black = C 0 0 0
  
  addCC :: Color -> Color -> Color
  addCC (C a b c) (C d e f) = C (a+d) (b+e) (c+f)
  
  subCC :: Color -> Color -> Color
  subCC (C a b c) (C d e f) = C (a-d) (b-e) (c-f)
  
  sumCC :: [Color] -> Color
  sumCC = foldr addCC black
  
  multCC :: Color -> Color -> Color
  multCC (C a b c) (C d e f) = C (a*d) (b*e) (c*f)
  
  multSC :: Double -> Color -> Color
  multSC k       (C a b c) = C (a*k) (b*k) (c*k)
  
  nearC :: Color -> Color -> Bool
  nearC (C a b c) (C d e f) = a `near` d && b `near` e && c `near` f
  
  offsetToPoint :: Ray -> Double -> Point
  offsetToPoint (r,v) i = r `addPV` (i `multSV` v)
  
  --
  
  epsilon, inf :: Double      -- aproximate zero and infinity
  epsilon = 1.0e-10
  inf = 1.0e20
  
  nonZero :: Double -> Double         -- Use before a division. It makes definitions
  nonZero x | x > epsilon  = x        -- more complete and I bet the errors that get
            | x < -epsilon = x        -- introduced will be undetectable if epsilon
            | otherwise    = epsilon  -- is small enough
  
  
  eqEps x y = abs (x-y) < epsilon
  near = eqEps
  
  clampf :: Double -> Double
  clampf p | p < 0 = 0
           | p > 1 = 1
           | True  = p
  