src/Cube.hs

   1 module Cube
   2 where
   3
   4 import Data.List ( (\\) )
   5 import Test.QuickCheck (Arbitrary(..), Gen, Positive(..), choose)
   6
   7 import Cardinal
   8 import qualified Face (Face(Face, v0, v1, v2, v3))
   9 import FunctionValues
  10 import Point
  11 import Tetrahedron hiding (c)
  12 import ThreeDimensional
  13
  14 data Cube = Cube { h :: Double,
  15                    i :: Int,
  16                    j :: Int,
  17                    k :: Int,
  18                    fv :: FunctionValues,
  19                    tetrahedra_volume :: Double }
  20             deriving (Eq)
  21
  22
  23 instance Arbitrary Cube where
  24     arbitrary = do
  25       (Positive h') <- arbitrary :: Gen (Positive Double)
  26       i' <- choose (coordmin, coordmax)
  27       j' <- choose (coordmin, coordmax)
  28       k' <- choose (coordmin, coordmax)
  29       fv' <- arbitrary :: Gen FunctionValues
  30       (Positive tet_vol) <- arbitrary :: Gen (Positive Double)
  31       return (Cube h' i' j' k' fv' tet_vol)
  32         where
  33           coordmin = -268435456 -- -(2^29 / 2)
  34           coordmax = 268435456  -- +(2^29 / 2)
  35
  36
  37 instance Show Cube where
  38     show c =
  39         "Cube_" ++ subscript ++ "\n" ++
  40         " h: " ++ (show (h c)) ++ "\n" ++
  41         " Center: " ++ (show (center c)) ++ "\n" ++
  42         " xmin: " ++ (show (xmin c)) ++ "\n" ++
  43         " xmax: " ++ (show (xmax c)) ++ "\n" ++
  44         " ymin: " ++ (show (ymin c)) ++ "\n" ++
  45         " ymax: " ++ (show (ymax c)) ++ "\n" ++
  46         " zmin: " ++ (show (zmin c)) ++ "\n" ++
  47         " zmax: " ++ (show (zmax c)) ++ "\n" ++
  48         " fv: " ++ (show (Cube.fv c)) ++ "\n"
  49         where
  50           subscript =
  51               (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c))
  52
  53
  54 -- | Returns an empty 'Cube'.
  55 empty_cube :: Cube
  56 empty_cube = Cube 0 0 0 0 empty_values 0
  57
  58
  59 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
  60 --   p. 76.
  61 xmin :: Cube -> Double
  62 xmin c = (2*i' - 1)*delta / 2
  63     where
  64       i' = fromIntegral (i c) :: Double
  65       delta = h c
  66
  67 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
  68 --   p. 76.
  69 xmax :: Cube -> Double
  70 xmax c = (2*i' + 1)*delta / 2
  71     where
  72       i' = fromIntegral (i c) :: Double
  73       delta = h c
  74
  75 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
  76 --   p. 76.
  77 ymin :: Cube -> Double
  78 ymin c = (2*j' - 1)*delta / 2
  79     where
  80       j' = fromIntegral (j c) :: Double
  81       delta = h c
  82
  83 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
  84 --   p. 76.
  85 ymax :: Cube -> Double
  86 ymax c = (2*j' + 1)*delta / 2
  87     where
  88       j' = fromIntegral (j c) :: Double
  89       delta = h c
  90
  91 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
  92 --   p. 76.
  93 zmin :: Cube -> Double
  94 zmin c = (2*k' - 1)*delta / 2
  95     where
  96       k' = fromIntegral (k c) :: Double
  97       delta = h c
  98
  99 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
 100 --   p. 76.
 101 zmax :: Cube -> Double
 102 zmax c = (2*k' + 1)*delta / 2
 103     where
 104       k' = fromIntegral (k c) :: Double
 105       delta = h c
 106
 107 instance ThreeDimensional Cube where
 108     -- | The center of Cube_ijk coincides with v_ijk at
 109     --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
 110     center c = (x, y, z)
 111            where
 112              delta = h c
 113              i' = fromIntegral (i c) :: Double
 114              j' = fromIntegral (j c) :: Double
 115              k' = fromIntegral (k c) :: Double
 116              x = delta * i'
 117              y = delta * j'
 118              z = delta * k'
 119
 120     -- | It's easy to tell if a point is within a cube; just make sure
 121     --   that it falls on the proper side of each of the cube's faces.
 122     contains_point c (x, y, z)
 123         | x < (xmin c) = False
 124         | x > (xmax c) = False
 125         | y < (ymin c) = False
 126         | y > (ymax c) = False
 127         | z < (zmin c) = False
 128         | z > (zmax c) = False
 129         | otherwise = True
 130
 131
 132
 133 -- Face stuff.
 134
 135 -- | The top (in the direction of z) face of the cube.
 136 top_face :: Cube -> Face.Face
 137 top_face c = Face.Face v0' v1' v2' v3'
 138     where
 139       delta = (1/2)*(h c)
 140       v0' = (center c) + (delta, -delta, delta)
 141       v1' = (center c) + (delta, delta, delta)
 142       v2' = (center c) + (-delta, delta, delta)
 143       v3' = (center c) + (-delta, -delta, delta)
 144
 145
 146
 147 -- | The back (in the direction of x) face of the cube.
 148 back_face :: Cube -> Face.Face
 149 back_face c = Face.Face v0' v1' v2' v3'
 150     where
 151       delta = (1/2)*(h c)
 152       v0' = (center c) + (delta, -delta, -delta)
 153       v1' = (center c) + (delta, delta, -delta)
 154       v2' = (center c) + (delta, delta, delta)
 155       v3' = (center c) + (delta, -delta, delta)
 156
 157
 158 -- The bottom face (in the direction of -z) of the cube.
 159 down_face :: Cube -> Face.Face
 160 down_face c = Face.Face v0' v1' v2' v3'
 161     where
 162       delta = (1/2)*(h c)
 163       v0' = (center c) + (-delta, -delta, -delta)
 164       v1' = (center c) + (-delta, delta, -delta)
 165       v2' = (center c) + (delta, delta, -delta)
 166       v3' = (center c) + (delta, -delta, -delta)
 167
 168
 169
 170 -- | The front (in the direction of -x) face of the cube.
 171 front_face :: Cube -> Face.Face
 172 front_face c = Face.Face v0' v1' v2' v3'
 173     where
 174       delta = (1/2)*(h c)
 175       v0' = (center c) + (-delta, -delta, delta)
 176       v1' = (center c) + (-delta, delta, delta)
 177       v2' = (center c) + (-delta, delta, -delta)
 178       v3' = (center c) + (-delta, -delta, -delta)
 179
 180 -- | The left (in the direction of -y) face of the cube.
 181 left_face :: Cube -> Face.Face
 182 left_face c = Face.Face v0' v1' v2' v3'
 183     where
 184       delta = (1/2)*(h c)
 185       v0' = (center c) + (delta, -delta, delta)
 186       v1' = (center c) + (-delta, -delta, delta)
 187       v2' = (center c) + (-delta, -delta, -delta)
 188       v3' = (center c) + (delta, -delta, -delta)
 189
 190
 191 -- | The right (in the direction of y) face of the cube.
 192 right_face :: Cube -> Face.Face
 193 right_face c = Face.Face v0' v1' v2' v3'
 194     where
 195       delta = (1/2)*(h c)
 196       v0' = (center c) + (-delta, delta, delta)
 197       v1' = (center c) + (delta, delta, delta)
 198       v2' = (center c) + (delta, delta, -delta)
 199       v3' = (center c) + (-delta, delta, -delta)
 200
 201
 202 tetrahedron0 :: Cube -> Tetrahedron
 203 tetrahedron0 c =
 204     Tetrahedron (Cube.fv c) v0' v1' v2' v3' vol
 205     where
 206       v0' = center c
 207       v1' = center (front_face c)
 208       v2' = Face.v0 (front_face c)
 209       v3' = Face.v1 (front_face c)
 210       vol = tetrahedra_volume c
 211
 212 tetrahedron1 :: Cube -> Tetrahedron
 213 tetrahedron1 c =
 214     Tetrahedron fv' v0' v1' v2' v3' vol
 215     where
 216       v0' = center c
 217       v1' = center (front_face c)
 218       v2' = Face.v1 (front_face c)
 219       v3' = Face.v2 (front_face c)
 220       fv' = rotate ccwx (Cube.fv c)
 221       vol = tetrahedra_volume c
 222
 223 tetrahedron2 :: Cube -> Tetrahedron
 224 tetrahedron2 c =
 225     Tetrahedron fv' v0' v1' v2' v3' vol
 226     where
 227       v0' = center c
 228       v1' = center (front_face c)
 229       v2' = Face.v2 (front_face c)
 230       v3' = Face.v3 (front_face c)
 231       fv' = rotate ccwx $ rotate ccwx $ Cube.fv c
 232       vol = tetrahedra_volume c
 233
 234 tetrahedron3 :: Cube -> Tetrahedron
 235 tetrahedron3 c =
 236     Tetrahedron fv' v0' v1' v2' v3' vol
 237     where
 238       v0' = center c
 239       v1' = center (front_face c)
 240       v2' = Face.v3 (front_face c)
 241       v3' = Face.v0 (front_face c)
 242       fv' = rotate cwx (Cube.fv c)
 243       vol = tetrahedra_volume c
 244
 245 tetrahedron4 :: Cube -> Tetrahedron
 246 tetrahedron4 c =
 247     Tetrahedron fv' v0' v1' v2' v3' vol
 248     where
 249       v0' = center c
 250       v1' = center (top_face c)
 251       v2' = Face.v0 (top_face c)
 252       v3' = Face.v1 (top_face c)
 253       fv' = rotate cwy (Cube.fv c)
 254       vol = tetrahedra_volume c
 255
 256 tetrahedron5 :: Cube -> Tetrahedron
 257 tetrahedron5 c =
 258     Tetrahedron fv' v0' v1' v2' v3' vol
 259     where
 260       v0' = center c
 261       v1' = center (top_face c)
 262       v2' = Face.v1 (top_face c)
 263       v3' = Face.v2 (top_face c)
 264       fv' = rotate cwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
 265       vol = tetrahedra_volume c
 266
 267 tetrahedron6 :: Cube -> Tetrahedron
 268 tetrahedron6 c =
 269     Tetrahedron fv' v0' v1' v2' v3' vol
 270     where
 271       v0' = center c
 272       v1' = center (top_face c)
 273       v2' = Face.v2 (top_face c)
 274       v3' = Face.v3 (top_face c)
 275       fv' = rotate cwy $ rotate cwz
 276                        $ rotate cwz
 277                        $ Tetrahedron.fv (tetrahedron0 c)
 278       vol = tetrahedra_volume c
 279
 280 tetrahedron7 :: Cube -> Tetrahedron
 281 tetrahedron7 c =
 282     Tetrahedron fv' v0' v1' v2' v3' vol
 283     where
 284       v0' = center c
 285       v1' = center (top_face c)
 286       v2' = Face.v3 (top_face c)
 287       v3' = Face.v0 (top_face c)
 288       fv' = rotate cwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
 289       vol = tetrahedra_volume c
 290
 291 tetrahedron8 :: Cube -> Tetrahedron
 292 tetrahedron8 c =
 293     Tetrahedron fv' v0' v1' v2' v3' vol
 294     where
 295       v0' = center c
 296       v1' = center (back_face c)
 297       v2' = Face.v0 (back_face c)
 298       v3' = Face.v1 (back_face c)
 299       fv' = rotate cwy $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
 300       vol = tetrahedra_volume c
 301
 302 tetrahedron9 :: Cube -> Tetrahedron
 303 tetrahedron9 c =
 304     Tetrahedron fv' v0' v1' v2' v3' vol
 305     where
 306       v0' = center c
 307       v1' = center (back_face c)
 308       v2' = Face.v1 (back_face c)
 309       v3' = Face.v2 (back_face c)
 310       fv' = rotate cwy $ rotate cwy
 311                        $ rotate cwx
 312                        $ Tetrahedron.fv (tetrahedron0 c)
 313       vol = tetrahedra_volume c
 314
 315 tetrahedron10 :: Cube -> Tetrahedron
 316 tetrahedron10 c =
 317     Tetrahedron fv' v0' v1' v2' v3' vol
 318     where
 319       v0' = center c
 320       v1' = center (back_face c)
 321       v2' = Face.v2 (back_face c)
 322       v3' = Face.v3 (back_face c)
 323       fv' = rotate cwy $ rotate cwy
 324                        $ rotate cwx
 325                        $ rotate cwx
 326                        $ Tetrahedron.fv (tetrahedron0 c)
 327
 328       vol = tetrahedra_volume c
 329
 330 tetrahedron11 :: Cube -> Tetrahedron
 331 tetrahedron11 c =
 332     Tetrahedron fv' v0' v1' v2' v3' vol
 333     where
 334       v0' = center c
 335       v1' = center (back_face c)
 336       v2' = Face.v3 (back_face c)
 337       v3' = Face.v0 (back_face c)
 338       fv' = rotate cwy $ rotate cwy
 339                        $ rotate ccwx
 340                        $ Tetrahedron.fv (tetrahedron0 c)
 341       vol = tetrahedra_volume c
 342
 343
 344 tetrahedron12 :: Cube -> Tetrahedron
 345 tetrahedron12 c =
 346     Tetrahedron fv' v0' v1' v2' v3' vol
 347     where
 348       v0' = center c
 349       v1' = center (down_face c)
 350       v2' = Face.v0 (down_face c)
 351       v3' = Face.v1 (down_face c)
 352       fv' = rotate ccwy (Tetrahedron.fv (tetrahedron0 c))
 353       vol = tetrahedra_volume c
 354
 355
 356 tetrahedron13 :: Cube -> Tetrahedron
 357 tetrahedron13 c =
 358     Tetrahedron fv' v0' v1' v2' v3' vol
 359     where
 360       v0' = center c
 361       v1' = center (down_face c)
 362       v2' = Face.v1 (down_face c)
 363       v3' = Face.v2 (down_face c)
 364       fv' = rotate ccwy $ rotate ccwz $ Tetrahedron.fv (tetrahedron0 c)
 365       vol = tetrahedra_volume c
 366
 367
 368 tetrahedron14 :: Cube -> Tetrahedron
 369 tetrahedron14 c =
 370     Tetrahedron fv' v0' v1' v2' v3' vol
 371     where
 372       v0' = center c
 373       v1' = center (down_face c)
 374       v2' = Face.v2 (down_face c)
 375       v3' = Face.v3 (down_face c)
 376       fv' = rotate ccwy $ rotate ccwz
 377                         $ rotate ccwz
 378                         $ Tetrahedron.fv (tetrahedron0 c)
 379       vol = tetrahedra_volume c
 380
 381
 382 tetrahedron15 :: Cube -> Tetrahedron
 383 tetrahedron15 c =
 384     Tetrahedron fv' v0' v1' v2' v3' vol
 385     where
 386       v0' = center c
 387       v1' = center (down_face c)
 388       v2' = Face.v3 (down_face c)
 389       v3' = Face.v0 (down_face c)
 390       fv' = rotate ccwy $ rotate cwz $ Tetrahedron.fv (tetrahedron0 c)
 391       vol = tetrahedra_volume c
 392
 393
 394 tetrahedron16 :: Cube -> Tetrahedron
 395 tetrahedron16 c =
 396     Tetrahedron fv' v0' v1' v2' v3' vol
 397     where
 398       v0' = center c
 399       v1' = center (right_face c)
 400       v2' = Face.v0 (right_face c)
 401       v3' = Face.v1 (right_face c)
 402       fv' = rotate ccwz (Tetrahedron.fv (tetrahedron0 c))
 403       vol = tetrahedra_volume c
 404
 405
 406 tetrahedron17 :: Cube -> Tetrahedron
 407 tetrahedron17 c =
 408     Tetrahedron fv' v0' v1' v2' v3' vol
 409     where
 410       v0' = center c
 411       v1' = center (right_face c)
 412       v2' = Face.v1 (right_face c)
 413       v3' = Face.v2 (right_face c)
 414       fv' = rotate ccwz $ rotate cwy $ Tetrahedron.fv (tetrahedron0 c)
 415       vol = tetrahedra_volume c
 416
 417
 418 tetrahedron18 :: Cube -> Tetrahedron
 419 tetrahedron18 c =
 420     Tetrahedron fv' v0' v1' v2' v3' vol
 421     where
 422       v0' = center c
 423       v1' = center (right_face c)
 424       v2' = Face.v2 (right_face c)
 425       v3' = Face.v3 (right_face c)
 426       fv' = rotate ccwz $ rotate cwy
 427                         $ rotate cwy
 428                         $ Tetrahedron.fv (tetrahedron0 c)
 429       vol = tetrahedra_volume c
 430
 431
 432 tetrahedron19 :: Cube -> Tetrahedron
 433 tetrahedron19 c =
 434     Tetrahedron fv' v0' v1' v2' v3' vol
 435     where
 436       v0' = center c
 437       v1' = center (right_face c)
 438       v2' = Face.v3 (right_face c)
 439       v3' = Face.v0 (right_face c)
 440       fv' = rotate ccwz $ rotate ccwy
 441                         $ Tetrahedron.fv (tetrahedron0 c)
 442       vol = tetrahedra_volume c
 443
 444
 445 tetrahedron20 :: Cube -> Tetrahedron
 446 tetrahedron20 c =
 447     Tetrahedron fv' v0' v1' v2' v3' vol
 448     where
 449       v0' = center c
 450       v1' = center (left_face c)
 451       v2' = Face.v0 (left_face c)
 452       v3' = Face.v1 (left_face c)
 453       fv' = rotate cwz (Tetrahedron.fv (tetrahedron0 c))
 454       vol = tetrahedra_volume c
 455
 456
 457 tetrahedron21 :: Cube -> Tetrahedron
 458 tetrahedron21 c =
 459     Tetrahedron fv' v0' v1' v2' v3' vol
 460     where
 461       v0' = center c
 462       v1' = center (left_face c)
 463       v2' = Face.v1 (left_face c)
 464       v3' = Face.v2 (left_face c)
 465       fv' = rotate cwz $ rotate ccwy $ Tetrahedron.fv (tetrahedron0 c)
 466       vol = tetrahedra_volume c
 467
 468
 469 tetrahedron22 :: Cube -> Tetrahedron
 470 tetrahedron22 c =
 471     Tetrahedron fv' v0' v1' v2' v3' vol
 472     where
 473       v0' = center c
 474       v1' = center (left_face c)
 475       v2' = Face.v2 (left_face c)
 476       v3' = Face.v3 (left_face c)
 477       fv' = rotate cwz $ rotate ccwy
 478                        $ rotate ccwy
 479                        $ Tetrahedron.fv (tetrahedron0 c)
 480       vol = tetrahedra_volume c
 481
 482
 483 tetrahedron23 :: Cube -> Tetrahedron
 484 tetrahedron23 c =
 485     Tetrahedron fv' v0' v1' v2' v3' vol
 486     where
 487       v0' = center c
 488       v1' = center (left_face c)
 489       v2' = Face.v3 (left_face c)
 490       v3' = Face.v0 (left_face c)
 491       fv' = rotate cwz $ rotate cwy
 492                        $ Tetrahedron.fv (tetrahedron0 c)
 493       vol = tetrahedra_volume c
 494
 495
 496 tetrahedra :: Cube -> [Tetrahedron]
 497 tetrahedra c =
 498     [tetrahedron0 c,
 499      tetrahedron1 c,
 500      tetrahedron2 c,
 501      tetrahedron3 c,
 502      tetrahedron4 c,
 503      tetrahedron5 c,
 504      tetrahedron6 c,
 505      tetrahedron7 c,
 506      tetrahedron8 c,
 507      tetrahedron9 c,
 508      tetrahedron10 c,
 509      tetrahedron11 c,
 510      tetrahedron12 c,
 511      tetrahedron13 c,
 512      tetrahedron14 c,
 513      tetrahedron15 c,
 514      tetrahedron16 c,
 515      tetrahedron17 c,
 516      tetrahedron18 c,
 517      tetrahedron19 c,
 518      tetrahedron20 c,
 519      tetrahedron21 c,
 520      tetrahedron22 c,
 521      tetrahedron23 c]
 522
 523 -- | All completely contained in the front half of the cube.
 524 front_half_tetrahedra :: Cube -> [Tetrahedron]
 525 front_half_tetrahedra c =
 526   [tetrahedron0 c,
 527    tetrahedron1 c,
 528    tetrahedron2 c,
 529    tetrahedron3 c,
 530    tetrahedron6 c,
 531    tetrahedron12 c,
 532    tetrahedron19 c,
 533    tetrahedron21 c]
 534
 535
 536 -- | All tetrahedra completely contained in the top half of the cube.
 537 top_half_tetrahedra :: Cube -> [Tetrahedron]
 538 top_half_tetrahedra c =
 539   [tetrahedron4 c,
 540    tetrahedron5 c,
 541    tetrahedron6 c,
 542    tetrahedron7 c,
 543    tetrahedron0 c,
 544    tetrahedron10 c,
 545    tetrahedron16 c,
 546    tetrahedron20 c]
 547
 548
 549 -- | All tetrahedra completely contained in the back half of the cube.
 550 back_half_tetrahedra :: Cube -> [Tetrahedron]
 551 back_half_tetrahedra c =
 552   [tetrahedron8 c,
 553    tetrahedron9 c,
 554    tetrahedron10 c,
 555    tetrahedron11 c,
 556    tetrahedron4 c,
 557    tetrahedron14 c,
 558    tetrahedron17 c,
 559    tetrahedron23 c]
 560
 561
 562 -- | All tetrahedra completely contained in the down half of the cube.
 563 down_half_tetrahedra :: Cube -> [Tetrahedron]
 564 down_half_tetrahedra c =
 565   [tetrahedron12 c,
 566    tetrahedron13 c,
 567    tetrahedron14 c,
 568    tetrahedron15 c,
 569    tetrahedron2 c,
 570    tetrahedron8 c,
 571    tetrahedron18 c,
 572    tetrahedron22 c]
 573
 574
 575 -- | All tetrahedra completely contained in the right half of the cube.
 576 right_half_tetrahedra :: Cube -> [Tetrahedron]
 577 right_half_tetrahedra c =
 578   [tetrahedron16 c,
 579    tetrahedron17 c,
 580    tetrahedron18 c,
 581    tetrahedron19 c,
 582    tetrahedron1 c,
 583    tetrahedron5 c,
 584    tetrahedron9 c,
 585    tetrahedron13 c]
 586
 587
 588 -- | All tetrahedra completely contained in the left half of the cube.
 589 left_half_tetrahedra :: Cube -> [Tetrahedron]
 590 left_half_tetrahedra c =
 591   [tetrahedron20 c,
 592    tetrahedron21 c,
 593    tetrahedron22 c,
 594    tetrahedron23 c,
 595    tetrahedron3 c,
 596    tetrahedron7 c,
 597    tetrahedron11 c,
 598    tetrahedron15 c]
 599
 600
 601 in_top_half :: Cube -> Point -> Bool
 602 in_top_half c (_,_,z) =
 603   distance_from_top <= distance_from_bottom
 604   where
 605     distance_from_top = abs $ (zmax c) - z
 606     distance_from_bottom = abs $ (zmin c) - z
 607
 608 in_front_half :: Cube -> Point -> Bool
 609 in_front_half c (x,_,_) =
 610     distance_from_front <= distance_from_back
 611   where
 612     distance_from_front = abs $ (xmin c) - x
 613     distance_from_back = abs $ (xmax c) - x
 614
 615
 616 in_left_half :: Cube -> Point -> Bool
 617 in_left_half c (_,y,_) =
 618     distance_from_left <= distance_from_right
 619   where
 620     distance_from_left = abs $ (ymin c) - y
 621     distance_from_right = abs $ (ymax c) - y
 622
 623
 624 -- | Takes a 'Cube', and returns the Tetrahedra belonging to it that
 625 --   contain the given 'Point'. This should be faster than checking
 626 --   every tetrahedron individually, since we determine which half
 627 --   (hemisphere?) of the cube the point lies in three times: once in
 628 --   each dimension. This allows us to eliminate non-candidates
 629 --   quickly.
 630 --
 631 --   This can throw an exception, but the use of 'head' might
 632 --   save us some unnecessary computations.
 633 --
 634 find_containing_tetrahedron :: Cube -> Point -> Tetrahedron
 635 find_containing_tetrahedron c p =
 636   head containing_tetrahedra
 637   where
 638     candidates = tetrahedra c
 639     non_candidates_x =
 640         if (in_front_half c p) then
 641           back_half_tetrahedra c
 642         else
 643           front_half_tetrahedra c
 644
 645     candidates_x = candidates \\ non_candidates_x
 646
 647     non_candidates_y =
 648       if (in_left_half c p) then
 649         right_half_tetrahedra c
 650       else
 651         left_half_tetrahedra c
 652
 653     candidates_xy = candidates_x \\ non_candidates_y
 654
 655     non_candidates_z =
 656       if (in_top_half c p) then
 657         down_half_tetrahedra c
 658       else
 659         top_half_tetrahedra c
 660
 661     candidates_xyz = candidates_xy \\ non_candidates_z
 662
 663     contains_our_point = flip contains_point p
 664     containing_tetrahedra = filter contains_our_point candidates_xyz