module Cube
where

import Data.List ( (\\) )
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,
                   tetrahedra_volume :: Double }
            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
      (Positive tet_vol) <- arbitrary :: Gen (Positive Double)
      return (Cube h' i' j' k' fv' tet_vol)
        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 0


-- | 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' vol
    where
      v0' = center c
      v1' = center (front_face c)
      v2' = Face.v0 (front_face c)
      v3' = Face.v1 (front_face c)
      vol = tetrahedra_volume c

tetrahedron1 :: Cube -> Tetrahedron
tetrahedron1 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c

tetrahedron2 :: Cube -> Tetrahedron
tetrahedron2 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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
      vol = tetrahedra_volume c

tetrahedron3 :: Cube -> Tetrahedron
tetrahedron3 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c

tetrahedron4 :: Cube -> Tetrahedron
tetrahedron4 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c

tetrahedron5 :: Cube -> Tetrahedron
tetrahedron5 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c

tetrahedron6 :: Cube -> Tetrahedron
tetrahedron6 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c

tetrahedron7 :: Cube -> Tetrahedron
tetrahedron7 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c

tetrahedron8 :: Cube -> Tetrahedron
tetrahedron8 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c

tetrahedron9 :: Cube -> Tetrahedron
tetrahedron9 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c

tetrahedron10 :: Cube -> Tetrahedron
tetrahedron10 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)

      vol = tetrahedra_volume c

tetrahedron11 :: Cube -> Tetrahedron
tetrahedron11 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c


tetrahedron12 :: Cube -> Tetrahedron
tetrahedron12 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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))
      vol = tetrahedra_volume c


tetrahedron13 :: Cube -> Tetrahedron
tetrahedron13 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c


tetrahedron14 :: Cube -> Tetrahedron
tetrahedron14 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c


tetrahedron15 :: Cube -> Tetrahedron
tetrahedron15 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c


tetrahedron16 :: Cube -> Tetrahedron
tetrahedron16 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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))
      vol = tetrahedra_volume c


tetrahedron17 :: Cube -> Tetrahedron
tetrahedron17 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c


tetrahedron18 :: Cube -> Tetrahedron
tetrahedron18 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c


tetrahedron19 :: Cube -> Tetrahedron
tetrahedron19 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c


tetrahedron20 :: Cube -> Tetrahedron
tetrahedron20 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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))
      vol = tetrahedra_volume c


tetrahedron21 :: Cube -> Tetrahedron
tetrahedron21 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c


tetrahedron22 :: Cube -> Tetrahedron
tetrahedron22 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume c


tetrahedron23 :: Cube -> Tetrahedron
tetrahedron23 c =
    Tetrahedron fv' v0' v1' v2' v3' vol
    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)
      vol = tetrahedra_volume 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]

-- | All completely contained in the front half of the cube.
front_half_tetrahedra :: Cube -> [Tetrahedron]
front_half_tetrahedra c =
  [tetrahedron0 c,
   tetrahedron1 c,
   tetrahedron2 c,
   tetrahedron3 c,
   tetrahedron6 c,
   tetrahedron12 c,
   tetrahedron19 c,
   tetrahedron21 c]


-- | All tetrahedra completely contained in the top half of the cube.
top_half_tetrahedra :: Cube -> [Tetrahedron]
top_half_tetrahedra c =
  [tetrahedron4 c,
   tetrahedron5 c,
   tetrahedron6 c,
   tetrahedron7 c,
   tetrahedron0 c,
   tetrahedron10 c,
   tetrahedron16 c,
   tetrahedron20 c]


-- | All tetrahedra completely contained in the back half of the cube.
back_half_tetrahedra :: Cube -> [Tetrahedron]
back_half_tetrahedra c =
  [tetrahedron8 c,
   tetrahedron9 c,
   tetrahedron10 c,
   tetrahedron11 c,
   tetrahedron4 c,
   tetrahedron14 c,
   tetrahedron17 c,
   tetrahedron23 c]


-- | All tetrahedra completely contained in the down half of the cube.
down_half_tetrahedra :: Cube -> [Tetrahedron]
down_half_tetrahedra c =
  [tetrahedron12 c,
   tetrahedron13 c,
   tetrahedron14 c,
   tetrahedron15 c,
   tetrahedron2 c,
   tetrahedron8 c,
   tetrahedron18 c,
   tetrahedron22 c]


-- | All tetrahedra completely contained in the right half of the cube.
right_half_tetrahedra :: Cube -> [Tetrahedron]
right_half_tetrahedra c =
  [tetrahedron16 c,
   tetrahedron17 c,
   tetrahedron18 c,
   tetrahedron19 c,
   tetrahedron1 c,
   tetrahedron5 c,
   tetrahedron9 c,
   tetrahedron13 c]


-- | All tetrahedra completely contained in the left half of the cube.
left_half_tetrahedra :: Cube -> [Tetrahedron]
left_half_tetrahedra c =
  [tetrahedron20 c,
   tetrahedron21 c,
   tetrahedron22 c,
   tetrahedron23 c,
   tetrahedron3 c,
   tetrahedron7 c,
   tetrahedron11 c,
   tetrahedron15 c]


in_top_half :: Cube -> Point -> Bool
in_top_half c (_,_,z) =
  distance_from_top <= distance_from_bottom
  where
    distance_from_top = abs $ (zmax c) - z
    distance_from_bottom = abs $ (zmin c) - z

in_front_half :: Cube -> Point -> Bool
in_front_half c (x,_,_) =
    distance_from_front <= distance_from_back
  where
    distance_from_front = abs $ (xmin c) - x
    distance_from_back = abs $ (xmax c) - x


in_left_half :: Cube -> Point -> Bool
in_left_half c (_,y,_) =
    distance_from_left <= distance_from_right
  where
    distance_from_left = abs $ (ymin c) - y
    distance_from_right = abs $ (ymax c) - y


-- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
--   contain the given 'Point'. This should be faster than checking
--   every tetrahedron individually, since we determine which half
--   (hemisphere?) of the cube the point lies in three times: once in
--   each dimension. This allows us to eliminate non-candidates
--   quickly.
--
--   This can throw an exception, but the use of 'head' might
--   save us some unnecessary computations.
--
find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
find_containing_tetrahedron c p =
  head containing_tetrahedra
  where
    candidates = tetrahedra c
    non_candidates_x =
        if (in_front_half c p) then
          back_half_tetrahedra c
        else
          front_half_tetrahedra c

    candidates_x = candidates \\ non_candidates_x

    non_candidates_y =
      if (in_left_half c p) then
        right_half_tetrahedra c
      else
        left_half_tetrahedra c

    candidates_xy = candidates_x \\ non_candidates_y

    non_candidates_z =
      if (in_top_half c p) then
        down_half_tetrahedra c
      else
        top_half_tetrahedra c

    candidates_xyz = candidates_xy \\ non_candidates_z

    contains_our_point = flip contains_point p
    containing_tetrahedra = filter contains_our_point candidates_xyz