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
  42 instance Gridded Cube where
  43     back c = cube_at (grid c) ((i c) + 1) (j c) (k c)
  44     down c = cube_at (grid c) (i c) (j c) ((k c) - 1)
  45     front c = cube_at (grid c) ((i c) - 1) (j c) (k c)
  46     left c = cube_at (grid c) (i c) ((j c) - 1) (k c)
  47     right c = cube_at (grid c) (i c) ((j c) + 1) (k c)
  48     top c = cube_at (grid c) (i c) (j c) ((k c) + 1)
  49
  50 -- | The left-side boundary of the cube. See Sorokina and Zeilfelder,
  51 --   p. 76.
  52 xmin :: Cube -> Double
  53 xmin c = (2*i' - 1)*delta / 2
  54     where
  55       i' = fromIntegral (i c) :: Double
  56       delta = h (grid c)
  57
  58 -- | The right-side boundary of the cube. See Sorokina and Zeilfelder,
  59 --   p. 76.
  60 xmax :: Cube -> Double
  61 xmax c = (2*i' + 1)*delta / 2
  62     where
  63       i' = fromIntegral (i c) :: Double
  64       delta = h (grid c)
  65
  66 -- | The front boundary of the cube. See Sorokina and Zeilfelder,
  67 --   p. 76.
  68 ymin :: Cube -> Double
  69 ymin c = (2*j' - 1)*delta / 2
  70     where
  71       j' = fromIntegral (j c) :: Double
  72       delta = h (grid c)
  73
  74 -- | The back boundary of the cube. See Sorokina and Zeilfelder,
  75 --   p. 76.
  76 ymax :: Cube -> Double
  77 ymax c = (2*j' + 1)*delta / 2
  78     where
  79       j' = fromIntegral (j c) :: Double
  80       delta = h (grid c)
  81
  82 -- | The bottom boundary of the cube. See Sorokina and Zeilfelder,
  83 --   p. 76.
  84 zmin :: Cube -> Double
  85 zmin c = (2*k' - 1)*delta / 2
  86     where
  87       k' = fromIntegral (k c) :: Double
  88       delta = h (grid c)
  89
  90 -- | The top boundary of the cube. See Sorokina and Zeilfelder,
  91 --   p. 76.
  92 zmax :: Cube -> Double
  93 zmax c = (2*k' + 1)*delta / 2
  94     where
  95       k' = fromIntegral (k c) :: Double
  96       delta = h (grid c)
  97
  98 instance ThreeDimensional Cube where
  99     -- | The center of Cube_ijk coincides with v_ijk at
 100     --   (ih, jh, kh). See Sorokina and Zeilfelder, p. 76.
 101     center c = (x, y, z)
 102            where
 103              delta = h (grid c)
 104              i' = fromIntegral (i c) :: Double
 105              j' = fromIntegral (j c) :: Double
 106              k' = fromIntegral (k c) :: Double
 107              x = delta * i'
 108              y = delta * j'
 109              z = delta * k'
 110
 111     contains_point c p
 112         | (x_coord p) < (xmin c) = False
 113         | (x_coord p) > (xmax c) = False
 114         | (y_coord p) < (ymin c) = False
 115         | (y_coord p) > (ymax c) = False
 116         | (z_coord p) < (zmin c) = False
 117         | (z_coord p) > (zmax c) = False
 118         | otherwise = True
 119
 120
 121 instance Num Cube where
 122     (Cube g1 i1 j1 k1 d1) + (Cube _ i2 j2 k2 d2) =
 123         Cube g1 (i1 + i2) (j1 + j2) (k1 + k2) (d1 + d2)
 124
 125     (Cube g1 i1 j1 k1 d1) - (Cube _ i2 j2 k2 d2) =
 126         Cube g1 (i1 - i2) (j1 - j2) (k1 - k2) (d1 - d2)
 127
 128     (Cube g1 i1 j1 k1 d1) * (Cube _ i2 j2 k2 d2) =
 129         Cube g1 (i1 * i2) (j1 * j2) (k1 * k2) (d1 * d2)
 130
 131     abs (Cube g1 i1 j1 k1 d1) =
 132         Cube g1 (abs i1) (abs j1) (abs k1) (abs d1)
 133
 134     signum (Cube g1 i1 j1 k1 d1) =
 135         Cube g1 (signum i1) (signum j1) (signum k1) (signum d1)
 136
 137     fromInteger x = empty_cube { datum = (fromIntegral x) }
 138
 139 instance Fractional Cube where
 140     (Cube g1 i1 j1 k1 d1) / (Cube _ _ _ _ d2) =
 141         Cube g1 i1 j1 k1 (d1 / d2)
 142
 143     recip (Cube g1 i1 j1 k1 d1) =
 144         Cube g1 i1 j1 k1 (recip d1)
 145
 146     fromRational q = empty_cube { datum = fromRational q }
 147
 148 -- | Constructs a cube, switching the x and z axes.
 149 reverse_cube :: Grid -> Int -> Int -> Int -> Double -> Cube
 150 reverse_cube g k' j' i' = Cube g i' j' k'
 151
 152
 153 -- | Return the cube corresponding to the grid point i,j,k. The list
 154 --   of cubes is stored as [z][y][x] but we'll be requesting it by
 155 --   [x][y][z] so we flip the indices in the last line.
 156 cube_at :: Grid -> Int -> Int -> Int -> Cube
 157 cube_at g i' j' k'
 158         | i' >= length (function_values g) = Cube g i' j' k' 0
 159         | i' < 0                          = Cube g i' j' k' 0
 160         | j' >= length ((function_values g) !! i') = Cube g i' j' k' 0
 161         | j' < 0                                  = Cube g i' j' k' 0
 162         | k' >= length (((function_values g) !! i') !! j') = Cube g i' j' k' 0
 163         | k' < 0                                          = Cube g i' j' k' 0
 164         | otherwise =
 165             (((cubes g) !! k') !! j') !! i'
 166
 167
 168 -- These next three functions basically form a 'for' loop, looping
 169 -- through the xs, ys, and zs in that order.
 170
 171 -- | The cubes_from_values function will return a function that takes
 172 --   a list of values and returns a list of cubes. It could just as
 173 --   well be written to take the values as a parameter; the omission
 174 --   of the last parameter is known as an eta reduce.
 175 cubes_from_values :: Grid -> Int -> Int -> ([Double] -> [Cube])
 176 cubes_from_values g i' j' =
 177     zipWith (reverse_cube g i' j') [0..]
 178
 179 -- | The cubes_from_planes function will return a function that takes
 180 --   a list of planes and returns a list of cubes. It could just as
 181 --   well be written to take the planes as a parameter; the omission
 182 --   of the last parameter is known as an eta reduce.
 183 cubes_from_planes :: Grid -> Int -> ([[Double]] -> [[Cube]])
 184 cubes_from_planes g i' =
 185     zipWith (cubes_from_values g i') [0..]
 186
 187 -- | Takes a grid as an argument, and returns a three-dimensional list
 188 --   of cubes centered on its grid points.
 189 cubes :: Grid -> [[[Cube]]]
 190 cubes g = zipWith (cubes_from_planes g) [0..] (function_values g)