Fix value_at for 2d slices.

diff --git a/src/FunctionValues.hs b/src/FunctionValues.hs
index 8400c80988bf30a43bacbfdecc0cd83abe0d0b30..74e212386a77b535b42c290b5f52180f63ca207a 100644 (file)
--- a/src/FunctionValues.hs
+++ b/src/FunctionValues.hs
@@ -135,6 +135,7 @@ empty_values :: FunctionValues
 empty_values =
     FunctionValues 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
+
 -- | The eval function is where the magic happens for the
 --   FunctionValues type. Given a 'Cardinal' direction and a
 --   'FunctionValues' object, eval will return the value of the
@@ -175,6 +176,7 @@ eval f (Difference x y) = (eval f x) - (eval f y)
 eval f (Product x y) = (eval f x) * (eval f y)
 eval f (Quotient x y) = (eval f x) / (eval f y)
 
+
 -- | Takes a three-dimensional list of 'Double' and a set of 3D
 --   coordinates (i,j,k), and returns the value at (i,j,k) in the
 --   supplied list. If there is no such value, we calculate one
@@ -207,25 +209,46 @@ value_at v3d i j k
   -- have been added where the indices are one-too-big. These are the
   -- "one index is bad" cases.
   | not (valid_i i) =
-      if (i == -1)
+      if (dim_i == 1)
       then
-        2*(value_at v3d 0 j k) - (value_at v3d 1 j k)
+        -- We're one-dimensional in our first coordinate, so just
+        -- return the data point that we do have. If we try to use
+        -- the formula from remark 7.3, we go into an infinite loop.
+        value_at v3d 0 j k
       else
-        2*(value_at v3d (i-1) j k) - (value_at v3d (i-2) j k)
+        if (i == -1)
+        then
+          2*(value_at v3d 0 j k) - (value_at v3d 1 j k)
+        else
+          2*(value_at v3d (i-1) j k) - (value_at v3d (i-2) j k)
 
   | not (valid_j j) =
-      if (j == -1)
+      if (dim_j == 1)
       then
-        2*(value_at v3d i 0 k) - (value_at v3d i 1 k)
+        -- We're one-dimensional in our second coordinate, so just
+        -- return the data point that we do have. If we try to use
+        -- the formula from remark 7.3, we go into an infinite loop.
+        value_at v3d i 0 k
       else
-        2*(value_at v3d i (j-1) k) - (value_at v3d i (j-2) k)
+        if (j == -1)
+        then
+          2*(value_at v3d i 0 k) - (value_at v3d i 1 k)
+        else
+          2*(value_at v3d i (j-1) k) - (value_at v3d i (j-2) k)
 
   | not (valid_k k) =
-      if (k == -1)
+      if (dim_k == 1)
       then
-        2*(value_at v3d i j 0) - (value_at v3d i j 1)
+        -- We're one-dimensional in our third coordinate, so just
+        -- return the data point that we do have. If we try to use
+        -- the formula from remark 7.3, we go into an infinite loop.
+        value_at v3d i j 0
       else
-        2*(value_at v3d i j (k-1)) - (value_at v3d i j (k-2))
+        if (k == -1)
+        then
+          2*(value_at v3d i j 0) - (value_at v3d i j 1)
+        else
+          2*(value_at v3d i j (k-1)) - (value_at v3d i j (k-2))
 
   | otherwise =
       let istr = show i
@@ -235,16 +258,16 @@ value_at v3d i j k
       in
         error $ "value_at called outside of domain: " ++ coordstr
   where
-    (xsize, ysize, zsize) = dims v3d
+    (dim_i, dim_j, dim_k) = dims v3d
 
     valid_i :: Int -> Bool
-    valid_i i' = (i' >= 0) && (i' < xsize)
+    valid_i i' = (i' >= 0) && (i' < dim_i)
 
     valid_j :: Int -> Bool
-    valid_j j' = (j' >= 0) && (j' < ysize)
+    valid_j j' = (j' >= 0) && (j' < dim_j)
 
     valid_k :: Int -> Bool
-    valid_k k' = (k' >= 0) && (k' < zsize)
+    valid_k k' = (k' >= 0) && (k' < dim_k)