src/Cube.hs

   1 module Cube
   2 where
   3
   4 import Grid
   5 import Point
   6 import ThreeDimensional
   7
   8 class Gridded a where
   9     back :: a -> Cube
  10     down :: a -> Cube
  11     front :: a -> Cube
  12     left :: a -> Cube
  13     right :: a -> Cube
  14     top :: a -> Cube
  15
  16
  17 data Cube = Cube { grid :: Grid,
  18                    i :: Int,
  19                    j :: Int,
  20                    k :: Int,
  21                    datum :: Double }
  22             deriving (Eq)
  23
  24
  25 instance Show Cube where
  26     show c =
  27         "Cube_" ++ (show (i c)) ++ "," ++ (show (j c)) ++ "," ++ (show (k c)) ++
  28         " (Grid: " ++ (show (grid c)) ++ ")" ++
  29         " (Center: " ++ (show (center c)) ++ ")" ++
  30         " (xmin: " ++ (show (xmin c)) ++ ")" ++
  31         " (xmax: " ++ (show (xmax c)) ++ ")" ++
  32         " (ymin: " ++ (show (ymin c)) ++ ")" ++
  33         " (ymax: " ++ (show (ymax c)) ++ ")" ++
  34         " (zmin: " ++ (show (zmin c)) ++ ")" ++
  35         " (zmax: " ++ (show (zmax c)) ++ ")" ++
  36         " (datum: " ++ (show (datum c)) ++ ")\n\n"
  37
  38 empty_cube :: Cube
  39 empty_cube = Cube empty_grid 0 0 0 0
  40
  41 -- TODO: The paper considers 'i' to be the front/back direction,
  42 -- whereas I have it in the left/right direction.
  43 instance Gridded Cube where
  44     back c = cube_at (grid c) ((i c) + 1) (j c) (k c)
  45     down c = cube_at (grid c) (i c) (j c) ((k c) - 1)
  46     front c = cube_at (grid c) ((i c) - 1) (j c) (k c)
  47     left c = cube_at (grid c) (i c) ((j c) - 1) (k c)
  48     right c = cube_at (grid c) (i c) ((j c) + 1) (k c)
  49     top c = cube_at (grid c) (i c) (j c) ((k c) + 1)
  50
  51 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
  52 --   p. 76.
  53 xmin :: Cube -> Double
  54 xmin c = (2*i' - 1)*delta / 2
  55     where
  56       i' = fromIntegral (i c) :: Double
  57       delta = h (grid c)
  58
  59 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
  60 --   p. 76.
  61 xmax :: Cube -> Double
  62 xmax c = (2*i' + 1)*delta / 2
  63     where
  64       i' = fromIntegral (i c) :: Double
  65       delta = h (grid c)
  66
  67 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
  68 --   p. 76.
  69 ymin :: Cube -> Double
  70 ymin c = (2*j' - 1)*delta / 2
  71     where
  72       j' = fromIntegral (j c) :: Double
  73       delta = h (grid c)
  74
  75 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
  76 --   p. 76.
  77 ymax :: Cube -> Double
  78 ymax c = (2*j' + 1)*delta / 2
  79     where
  80       j' = fromIntegral (j c) :: Double
  81       delta = h (grid c)
  82
  83 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
  84 --   p. 76.
  85 zmin :: Cube -> Double
  86 zmin c = (2*k' - 1)*delta / 2
  87     where
  88       k' = fromIntegral (k c) :: Double
  89       delta = h (grid c)
  90
  91 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
  92 --   p. 76.
  93 zmax :: Cube -> Double
  94 zmax c = (2*k' + 1)*delta / 2
  95     where
  96       k' = fromIntegral (k c) :: Double
  97       delta = h (grid c)
  98
  99 instance ThreeDimensional Cube where
 100     -- | The center of Cube_ijk coincides with v_ijk at
 101     --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
 102     center c = (x, y, z)
 103            where
 104              delta = h (grid c)
 105              i' = fromIntegral (i c) :: Double
 106              j' = fromIntegral (j c) :: Double
 107              k' = fromIntegral (k c) :: Double
 108              x = (delta * i')
 109              y = (delta * j')
 110              z = (delta * k')
 111
 112     contains_point c p
 113         | (x_coord p) < (xmin c) = False
 114         | (x_coord p) > (xmax c) = False
 115         | (y_coord p) < (ymin c) = False
 116         | (y_coord p) > (ymax c) = False
 117         | (z_coord p) < (zmin c) = False
 118         | (z_coord p) > (zmax c) = False
 119         | otherwise = True
 120
 121
 122 instance Num Cube where
 123     (Cube g1 i1 j1 k1 d1) + (Cube _ i2 j2 k2 d2) =
 124         Cube g1 (i1 + i2) (j1 + j2) (k1 + k2) (d1 + d2)
 125
 126     (Cube g1 i1 j1 k1 d1) - (Cube _ i2 j2 k2 d2) =
 127         Cube g1 (i1 - i2) (j1 - j2) (k1 - k2) (d1 - d2)
 128
 129     (Cube g1 i1 j1 k1 d1) * (Cube _ i2 j2 k2 d2) =
 130         Cube g1 (i1 * i2) (j1 * j2) (k1 * k2) (d1 * d2)
 131
 132     abs (Cube g1 i1 j1 k1 d1) =
 133         Cube g1 (abs i1) (abs j1) (abs k1) (abs d1)
 134
 135     signum (Cube g1 i1 j1 k1 d1) =
 136         Cube g1 (signum i1) (signum j1) (signum k1) (signum d1)
 137
 138     fromInteger x = empty_cube { datum = (fromIntegral x) }
 139
 140 instance Fractional Cube where
 141     (Cube g1 i1 j1 k1 d1) / (Cube _ _ _ _ d2) =
 142         Cube g1 i1 j1 k1 (d1 / d2)
 143
 144     recip (Cube g1 i1 j1 k1 d1) =
 145         Cube g1 i1 j1 k1 (recip d1)
 146
 147     fromRational q = empty_cube { datum = fromRational q }
 148
 149 reverse_cube :: Grid -> Int -> Int -> Int -> Double -> Cube
 150 reverse_cube g k' j' i' datum' = Cube g i' j' k' datum'
 151
 152
 153 cube_at :: Grid -> Int -> Int -> Int -> Cube
 154 cube_at g i' j' k'
 155         | i' >= length (function_values g) = Cube g i' j' k' 0
 156         | i' < 0                          = Cube g i' j' k' 0
 157         | j' >= length ((function_values g) !! i') = Cube g i' j' k' 0
 158         | j' < 0                                  = Cube g i' j' k' 0
 159         | k' >= length (((function_values g) !! i') !! j') = Cube g i' j' k' 0
 160         | k' < 0                                          = Cube g i' j' k' 0
 161         | otherwise =
 162             Cube g i' j' k' ((((function_values g) !! i') !! j') !! k')
 163
 164 -- These next three functions basically form a 'for' loop, looping
 165 -- through the xs, ys, and zs in that order.
 166 cubes_from_values :: Grid -> Int -> Int -> [Double] -> [Cube]
 167 cubes_from_values g i' j' vals =
 168     zipWith (reverse_cube g i' j') [0..] vals
 169
 170 cubes_from_planes :: Grid -> Int -> [[Double]] -> [[Cube]]
 171 cubes_from_planes g i' planes =
 172     zipWith (cubes_from_values g i') [0..] planes
 173
 174 cubes :: Grid -> [[[Cube]]]
 175 cubes g = zipWith (cubes_from_planes g) [0..] (function_values g)