Use a slightly more stable (I think?) formula for the {x,y,z}{min,max}.

[spline3.git] / src / FunctionValues.hs

diff --git a/src/FunctionValues.hs b/src/FunctionValues.hs
index e9da25ffc235171c89bd6b878162d681554d5e0a..8400c80988bf30a43bacbfdecc0cd83abe0d0b30 100644 (file)
--- a/src/FunctionValues.hs
+++ b/src/FunctionValues.hs
@@ -177,8 +177,11 @@ 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 choose a nearby
---   point and use its value.
+--   supplied list. If there is no such value, we calculate one
+--   according to Sorokina and Zeilfelder, remark 7.3, p. 99.
+--
+--   We specifically do not consider values more than one unit away
+--   from our grid.
 --
 --   Examples:
 --
@@ -186,26 +189,64 @@ eval f (Quotient x y) = (eval f x) / (eval f y)
 --   1.0
 --
 --   >>> value_at Examples.trilinear (-1) 0 0
---   1.0
+--   0.0
 --
 --   >>> value_at Examples.trilinear 0 0 4
 --   1.0
 --
 --   >>> value_at Examples.trilinear 1 3 0
---   4.0
+--   5.0
 --
 value_at :: Values3D -> Int -> Int -> Int -> Double
 value_at v3d i j k
-         | i < 0 = value_at v3d 0 j k
-         | j < 0 = value_at v3d i 0 k
-         | k < 0 = value_at v3d i j 0
-         | xsize <= i = value_at v3d (xsize - 1) j k
-         | ysize <= j = value_at v3d i (ysize - 1) k
-         | zsize <= k = value_at v3d i j (zsize - 1)
-         | otherwise = idx v3d i j k
+  -- Put the most common case first!
+  | (valid_i i) && (valid_j j) && (valid_k k) =
+      idx v3d i j k
+
+  -- The next three are from the first line in (7.3). Analogous cases
+  -- 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)
+      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)
+      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)
+      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
+          jstr = show j
+          kstr = show k
+          coordstr = "(" ++ istr ++ "," ++ jstr ++ "," ++ kstr ++ ")"
+      in
+        error $ "value_at called outside of domain: " ++ coordstr
   where
     (xsize, ysize, zsize) = dims v3d
 
+    valid_i :: Int -> Bool
+    valid_i i' = (i' >= 0) && (i' < xsize)
+
+    valid_j :: Int -> Bool
+    valid_j j' = (j' >= 0) && (j' < ysize)
+
+    valid_k :: Int -> Bool
+    valid_k k' = (k' >= 0) && (k' < zsize)
+
+
 
 -- | Given a three-dimensional list of 'Double' and a set of 3D
 --   coordinates (i,j,k), constructs and returns the 'FunctionValues'