module Cube
where

import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
    
import Cardinal
import qualified Face (Face(Face, v0, v1, v2, v3))
import FunctionValues
import Point
import Tetrahedron hiding (c)
import ThreeDimensional

data Cube = Cube { h :: Double,
                   i :: Int,
                   j :: Int,
                   k :: Int,
                   fv :: FunctionValues }
            deriving (Eq)


instance Arbitrary Cube where
    arbitrary = do
      (Positive h') <- arbitrary :: Gen (Positive Double)
      i' <- choose (coordmin, coordmax)
      j' <- choose (coordmin, coordmax)
      k' <- choose (coordmin, coordmax)
      fv' <- arbitrary :: Gen FunctionValues
      return (Cube h' i' j' k' fv')
        where
          coordmin = -268435456 -- -(2^29 / 2)
          coordmax = 268435456  -- +(2^29 / 2)


instance Show Cube where
    show c =
        "Cube_" ++ subscript ++ "\n" ++
        " h: " ++ (show (h c)) ++ "\n" ++
        " Center: " ++ (show (center c)) ++ "\n" ++
        " xmin: " ++ (show (xmin c)) ++ "\n" ++
        " xmax: " ++ (show (xmax c)) ++ "\n" ++
        " ymin: " ++ (show (ymin c)) ++ "\n" ++
        " ymax: " ++ (show (ymax c)) ++ "\n" ++
        " zmin: " ++ (show (zmin c)) ++ "\n" ++
        " zmax: " ++ (show (zmax c)) ++ "\n" ++
        " fv: " ++ (show (Cube.fv c)) ++ "\n"
        where
          subscript =
              (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c))


-- | Returns an empty 'Cube'.
empty_cube :: Cube
empty_cube = Cube 0 0 0 0 empty_values


-- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
xmin :: Cube -> Double
xmin c = (2*i' - 1)*delta / 2
    where
      i' = fromIntegral (i c) :: Double
      delta = h c

-- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
xmax :: Cube -> Double
xmax c = (2*i' + 1)*delta / 2
    where
      i' = fromIntegral (i c) :: Double
      delta = h c

-- | The front boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
ymin :: Cube -> Double
ymin c = (2*j' - 1)*delta / 2
    where
      j' = fromIntegral (j c) :: Double
      delta = h c

-- | The back boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
ymax :: Cube -> Double
ymax c = (2*j' + 1)*delta / 2
    where
      j' = fromIntegral (j c) :: Double
      delta = h c

-- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
zmin :: Cube -> Double
zmin c = (2*k' - 1)*delta / 2
    where
      k' = fromIntegral (k c) :: Double
      delta = h c

-- | The top boundary of the cube. See Sorokina and Zeilfelder,
--   p. 76.
zmax :: Cube -> Double
zmax c = (2*k' + 1)*delta / 2
    where
      k' = fromIntegral (k c) :: Double
      delta = h c

instance ThreeDimensional Cube where
    -- | The center of Cube_ijk coincides with v_ijk at
    --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.    
    center c = (x, y, z)
           where
             delta = h c
             i' = fromIntegral (i c) :: Double
             j' = fromIntegral (j c) :: Double
             k' = fromIntegral (k c) :: Double
             x = delta * i'
             y = delta * j'
             z = delta * k'

    -- | It's easy to tell if a point is within a cube; just make sure
    --   that it falls on the proper side of each of the cube's faces.
    contains_point c (x, y, z)
        | x < (xmin c) = False
        | x > (xmax c) = False
        | y < (ymin c) = False
        | y > (ymax c) = False
        | z < (zmin c) = False
        | z > (zmax c) = False
        | otherwise = True



-- Face stuff.

-- | The top (in the direction of z) face of the cube.
top_face :: Cube -> Face.Face
top_face c = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2)*(h c)
      v0' = (center c) + (delta, -delta, delta)
      v1' = (center c) + (delta, delta, delta)
      v2' = (center c) + (-delta, delta, delta)
      v3' = (center c) + (-delta, -delta, delta)



-- | The back (in the direction of x) face of the cube.
back_face :: Cube -> Face.Face
back_face c = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2)*(h c)
      v0' = (center c) + (delta, -delta, -delta)
      v1' = (center c) + (delta, delta, -delta)
      v2' = (center c) + (delta, delta, delta)
      v3' = (center c) + (delta, -delta, delta)


-- The bottom face (in the direction of -z) of the cube.
down_face :: Cube -> Face.Face
down_face c = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2)*(h c)
      v0' = (center c) + (-delta, -delta, -delta)
      v1' = (center c) + (-delta, delta, -delta)
      v2' = (center c) + (delta, delta, -delta)
      v3' = (center c) + (delta, -delta, -delta)



-- | The front (in the direction of -x) face of the cube.
front_face :: Cube -> Face.Face
front_face c = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2)*(h c)
      v0' = (center c) + (-delta, -delta, delta)
      v1' = (center c) + (-delta, delta, delta)
      v2' = (center c) + (-delta, delta, -delta)
      v3' = (center c) + (-delta, -delta, -delta)

-- | The left (in the direction of -y) face of the cube.
left_face :: Cube -> Face.Face
left_face c = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2)*(h c)
      v0' = (center c) + (delta, -delta, delta)
      v1' = (center c) + (-delta, -delta, delta)
      v2' = (center c) + (-delta, -delta, -delta)
      v3' = (center c) + (delta, -delta, -delta)


-- | The right (in the direction of y) face of the cube.
right_face :: Cube -> Face.Face
right_face c = Face.Face v0' v1' v2' v3'
    where
      delta = (1/2)*(h c)
      v0' = (center c) + (-delta, delta, delta)
      v1' = (center c) + (delta, delta, delta)
      v2' = (center c) + (delta, delta, -delta)
      v3' = (center c) + (-delta, delta, -delta)


tetrahedron0 :: Cube -> Tetrahedron
tetrahedron0 c =
    Tetrahedron (Cube.fv c) v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (front_face c)
      v2' = Face.v0 (front_face c)
      v3' = Face.v1 (front_face c)

tetrahedron1 :: Cube -> Tetrahedron
tetrahedron1 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (front_face c)
      v2' = Face.v1 (front_face c)
      v3' = Face.v2 (front_face c)
      fv' = rotate ccwx (Cube.fv c)

tetrahedron2 :: Cube -> Tetrahedron
tetrahedron2 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (front_face c)
      v2' = Face.v2 (front_face c)
      v3' = Face.v3 (front_face c)
      fv' = rotate ccwx $ rotate ccwx $ Cube.fv c

tetrahedron3 :: Cube -> Tetrahedron
tetrahedron3 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (front_face c)
      v2' = Face.v3 (front_face c)
      v3' = Face.v0 (front_face c)
      fv' = rotate cwx (Cube.fv c)

tetrahedron4 :: Cube -> Tetrahedron
tetrahedron4 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (top_face c)
      v2' = Face.v0 (top_face c)
      v3' = Face.v1 (top_face c)
      fv' = rotate cwy (Cube.fv c)

tetrahedron5 :: Cube -> Tetrahedron
tetrahedron5 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (top_face c)
      v2' = Face.v1 (top_face c)
      v3' = Face.v2 (top_face c)
      fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)

tetrahedron6 :: Cube -> Tetrahedron
tetrahedron6 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (top_face c)
      v2' = Face.v2 (top_face c)
      v3' = Face.v3 (top_face c)
      fv' = rotate cwy $ rotate cwz $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)

tetrahedron7 :: Cube -> Tetrahedron
tetrahedron7 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (top_face c)
      v2' = Face.v3 (top_face c)
      v3' = Face.v0 (top_face c)
      fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)

tetrahedron8 :: Cube -> Tetrahedron
tetrahedron8 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (back_face c)
      v2' = Face.v0 (back_face c)
      v3' = Face.v1 (back_face c)
      fv' = rotate cwy $ rotate cwy $ (Tetrahedron.fv (tetrahedron0 c))

tetrahedron9 :: Cube -> Tetrahedron
tetrahedron9 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (back_face c)
      v2' = Face.v1 (back_face c)
      v3' = Face.v2 (back_face c)
      fv' = rotate cwy $ rotate cwy $ rotate cwx $ Tetrahedron.fv (tetrahedron0 c)

tetrahedron10 :: Cube -> Tetrahedron
tetrahedron10 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (back_face c)
      v2' = Face.v2 (back_face c)
      v3' = Face.v3 (back_face c)
      fv' = rotate cwy $ rotate cwy
                       $ rotate cwx
                       $ rotate cwx
                       $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron11 :: Cube -> Tetrahedron
tetrahedron11 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (back_face c)
      v2' = Face.v3 (back_face c)
      v3' = Face.v0 (back_face c)
      fv' = rotate cwy $ rotate cwy
                       $ rotate ccwx
                       $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron12 :: Cube -> Tetrahedron
tetrahedron12 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (down_face c)
      v2' = Face.v0 (down_face c)
      v3' = Face.v1 (down_face c)
      fv' = rotate ccwy (Tetrahedron.fv (tetrahedron0 c))


tetrahedron13 :: Cube -> Tetrahedron
tetrahedron13 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (down_face c)
      v2' = Face.v1 (down_face c)
      v3' = Face.v2 (down_face c)
      fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron14 :: Cube -> Tetrahedron
tetrahedron14 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (down_face c)
      v2' = Face.v2 (down_face c)
      v3' = Face.v3 (down_face c)
      fv' = rotate ccwy $ rotate ccwz
                        $ rotate ccwz
                        $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron15 :: Cube -> Tetrahedron
tetrahedron15 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (down_face c)
      v2' = Face.v3 (down_face c)
      v3' = Face.v0 (down_face c)
      fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron16 :: Cube -> Tetrahedron
tetrahedron16 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (right_face c)
      v2' = Face.v0 (right_face c)
      v3' = Face.v1 (right_face c)
      fv' = rotate ccwz (Tetrahedron.fv (tetrahedron0 c))


tetrahedron17 :: Cube -> Tetrahedron
tetrahedron17 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (right_face c)
      v2' = Face.v1 (right_face c)
      v3' = Face.v2 (right_face c)
      fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron18 :: Cube -> Tetrahedron
tetrahedron18 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (right_face c)
      v2' = Face.v2 (right_face c)
      v3' = Face.v3 (right_face c)
      fv' = rotate ccwz $ rotate cwy
                        $ rotate cwy
                        $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron19 :: Cube -> Tetrahedron
tetrahedron19 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (right_face c)
      v2' = Face.v3 (right_face c)
      v3' = Face.v0 (right_face c)
      fv' = rotate ccwz $ rotate ccwy
                        $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron20 :: Cube -> Tetrahedron
tetrahedron20 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (left_face c)
      v2' = Face.v0 (left_face c)
      v3' = Face.v1 (left_face c)
      fv' = rotate cwz (Tetrahedron.fv (tetrahedron0 c))


tetrahedron21 :: Cube -> Tetrahedron
tetrahedron21 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (left_face c)
      v2' = Face.v1 (left_face c)
      v3' = Face.v2 (left_face c)
      fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron22 :: Cube -> Tetrahedron
tetrahedron22 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (left_face c)
      v2' = Face.v2 (left_face c)
      v3' = Face.v3 (left_face c)
      fv' = rotate cwz $ rotate ccwy
                       $ rotate ccwy
                       $ Tetrahedron.fv (tetrahedron0 c)


tetrahedron23 :: Cube -> Tetrahedron
tetrahedron23 c =
    Tetrahedron fv' v0' v1' v2' v3'
    where
      v0' = center c
      v1' = center (left_face c)
      v2' = Face.v3 (left_face c)
      v3' = Face.v0 (left_face c)
      fv' = rotate cwz $ rotate cwy
                       $ Tetrahedron.fv (tetrahedron0 c)


tetrahedra :: Cube -> [Tetrahedron]
tetrahedra c =
    [tetrahedron0 c,
     tetrahedron1 c,
     tetrahedron2 c,
     tetrahedron3 c,
     tetrahedron4 c,
     tetrahedron5 c,
     tetrahedron6 c,
     tetrahedron7 c,
     tetrahedron8 c,
     tetrahedron9 c,
     tetrahedron10 c,
     tetrahedron11 c,
     tetrahedron12 c,
     tetrahedron13 c,
     tetrahedron14 c,
     tetrahedron15 c,
     tetrahedron16 c,
     tetrahedron17 c,
     tetrahedron18 c,
     tetrahedron19 c,
     tetrahedron20 c,
     tetrahedron21 c,
     tetrahedron22 c,
     tetrahedron23 c]


-- | Takes a 'Cube', and returns all Tetrahedra belonging to it that
--   contain the given 'Point'.
find_containing_tetrahedra :: Cube -> Point -> [Tetrahedron]
find_containing_tetrahedra c p =
    filter contains_our_point all_tetrahedra
    where
      contains_our_point = flip contains_point p
      all_tetrahedra = tetrahedra c