module Cube
where

import Data.Maybe (fromJust)
import qualified Data.Vector as V (
  Vector,
  findIndex,
  map,
  minimum,
  singleton,
  snoc,
  unsafeIndex
  )
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, fv)
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)


tetrahedron :: Cube -> Int -> Tetrahedron

tetrahedron c 0 =
    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

tetrahedron c 1 =
    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

tetrahedron c 2 =
    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

tetrahedron c 3 =
    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

tetrahedron c 4 =
    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

tetrahedron c 5 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 6 =
    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
                       $ fv c
      vol = tetrahedra_volume c

tetrahedron c 7 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 8 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 9 =
    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
                       $ fv c
      vol = tetrahedra_volume c

tetrahedron c 10 =
    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
                       $ fv c

      vol = tetrahedra_volume c

tetrahedron c 11 =
    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
                       $ fv c
      vol = tetrahedra_volume c

tetrahedron c 12 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 13 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 14 =
    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
                        $ fv c
      vol = tetrahedra_volume c

tetrahedron c 15 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 16 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 17 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 18 =
    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
                        $ fv c
      vol = tetrahedra_volume c

tetrahedron c 19 =
    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
                        $ fv c
      vol = tetrahedra_volume c

tetrahedron c 20 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 21 =
    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 $ fv c
      vol = tetrahedra_volume c

tetrahedron c 22 =
    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
                       $ fv c
      vol = tetrahedra_volume c

tetrahedron c 23 =
    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
                       $ fv c
      vol = tetrahedra_volume c

-- Feels dirty, but whatever.
tetrahedron _ _ = error "asked for a nonexistent tetrahedron"


-- Only used in tests, so we don't need the added speed
-- of Data.Vector.
tetrahedra :: Cube -> [Tetrahedron]
tetrahedra c = [ tetrahedron c n | n <- [0..23] ]

front_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
front_left_top_tetrahedra  c =
  V.singleton (tetrahedron c 0) `V.snoc`
    (tetrahedron c 3) `V.snoc`
    (tetrahedron c 6) `V.snoc`
    (tetrahedron c 7) `V.snoc`
    (tetrahedron c 20) `V.snoc`
    (tetrahedron c 21)

front_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
front_left_down_tetrahedra  c =
  V.singleton (tetrahedron c 0) `V.snoc`
    (tetrahedron c 2) `V.snoc`
    (tetrahedron c 3) `V.snoc`
    (tetrahedron c 12) `V.snoc`
    (tetrahedron c 15) `V.snoc`
    (tetrahedron c 21)

front_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
front_right_top_tetrahedra  c =
  V.singleton (tetrahedron c 0) `V.snoc`
    (tetrahedron c 1) `V.snoc`
    (tetrahedron c 5) `V.snoc`
    (tetrahedron c 6) `V.snoc`
    (tetrahedron c 16) `V.snoc`
    (tetrahedron c 19)

front_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
front_right_down_tetrahedra  c =
  V.singleton (tetrahedron c 1) `V.snoc`
    (tetrahedron c 2) `V.snoc`
    (tetrahedron c 12) `V.snoc`
    (tetrahedron c 13) `V.snoc`
    (tetrahedron c 18) `V.snoc`
    (tetrahedron c 19)

back_left_top_tetrahedra :: Cube -> V.Vector Tetrahedron
back_left_top_tetrahedra  c =
  V.singleton (tetrahedron c 0) `V.snoc`
    (tetrahedron c 3) `V.snoc`
    (tetrahedron c 6) `V.snoc`
    (tetrahedron c 7) `V.snoc`
    (tetrahedron c 20) `V.snoc`
    (tetrahedron c 21)

back_left_down_tetrahedra :: Cube -> V.Vector Tetrahedron
back_left_down_tetrahedra  c =
  V.singleton (tetrahedron c 8) `V.snoc`
    (tetrahedron c 11) `V.snoc`
    (tetrahedron c 14) `V.snoc`
    (tetrahedron c 15) `V.snoc`
    (tetrahedron c 22) `V.snoc`
    (tetrahedron c 23)

back_right_top_tetrahedra :: Cube -> V.Vector Tetrahedron
back_right_top_tetrahedra  c =
  V.singleton (tetrahedron c 4) `V.snoc`
    (tetrahedron c 5) `V.snoc`
    (tetrahedron c 9) `V.snoc`
    (tetrahedron c 10) `V.snoc`
    (tetrahedron c 16) `V.snoc`
    (tetrahedron c 17)

back_right_down_tetrahedra :: Cube -> V.Vector Tetrahedron
back_right_down_tetrahedra  c =
  V.singleton (tetrahedron c 8) `V.snoc`
    (tetrahedron c 9) `V.snoc`
    (tetrahedron c 13) `V.snoc`
    (tetrahedron c 14) `V.snoc`
    (tetrahedron c 17) `V.snoc`
    (tetrahedron c 18)

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 =
  candidates `V.unsafeIndex` (fromJust lucky_idx)
  where
    front_half = in_front_half c p
    top_half = in_top_half c p
    left_half = in_left_half c p

    candidates =
      if front_half then

        if left_half then
          if top_half then
            front_left_top_tetrahedra c
          else
            front_left_down_tetrahedra c
        else
          if top_half then
            front_right_top_tetrahedra c
          else
            front_right_down_tetrahedra c

      else -- bottom half

        if left_half then
          if top_half then
            back_left_top_tetrahedra c
          else
            back_left_down_tetrahedra c
        else
          if top_half then
            back_right_top_tetrahedra c
          else
            back_right_down_tetrahedra c

    -- Use the dot product instead of 'distance' here to save a
    -- sqrt(). So, "distances" below really means "distances squared."
    distances = V.map ((dot p) . center) candidates
    shortest_distance = V.minimum distances
    lucky_idx = V.findIndex
                  (\t -> (center t) `dot` p == shortest_distance)
                  candidates