diff --git a/Data/Glome/Texture.hs b/Data/Glome/Texture.hs
new file mode 100644
--- /dev/null
+++ b/Data/Glome/Texture.hs
@@ -0,0 +1,126 @@
+{-
+module SolidTexture (square_wave, triangle_wave, sine_wave, 
+                     stripe, noise, turbulence
+                    ) where -}
+
+module Data.Glome.Texture where
+import Data.Glome.Vec
+import Data.Array.IArray
+
+-- INTERPOLATION FUNCTIONS --
+square_wave :: Flt -> Flt
+square_wave x =
+ let offset = x - (fromIntegral (floor x))
+ in if offset < 0.5 then 0 else 1
+
+triangle_wave :: Flt -> Flt
+triangle_wave x =
+ let offset = x - (fromIntegral (floor x))
+ in if offset < 0.5 
+    then (offset*2)
+    else (2-(offset*2))
+
+sine_wave :: Flt -> Flt
+sine_wave x = (sin (x*2*pi))*0.5 + 0.5
+
+
+lump_wave :: Flt -> Flt
+lump_wave x = 1 - x*x*x
+
+-- SCALAR TEXTURE FUNCTIONS --
+
+-- These are simple solid texture functions that take a 
+-- point as argument and return a number 0 < n < 1
+
+stripe :: Vec -> (Flt -> Flt) -> (Vec -> Flt)
+stripe axis interp =
+ let len = vlen axis 
+ in
+  (\pos -> let offset = vdot pos axis 
+           in interp offset)
+
+
+-- PERLIN NOISE --
+
+-- (-6 t^5 + 15 t^4 - 10t^3 +1)
+-- "realistic ray tracing 2nd edition" inconsistent 
+-- on whether it should be t^5 or t^6,
+-- but t^5 works and t^6 doesn't.
+omega :: Flt -> Flt
+omega t_ = 
+ let t     = fabs t_
+     tsqr  = t*t
+     tcube = tsqr*t
+ in (-6)*tcube*tsqr + 15*tcube*t - 10*tcube + 1
+
+-- questionably random
+phi :: Array Int Int
+phi = listArray (0,11) [3,0,2,7,4,1,5,11,8,10,9,6]
+
+grad :: Array Int Vec
+grad = listArray (0,11) 
+         $ filter (\x -> let l = vlen x in l < 1.5 && l > 1.1) 
+                  [Vec x y z | x <- [(-1),0,1],
+                               y <- [(-1),0,1],
+                               z <- [(-1),0,1]] 
+
+gamma :: Int -> Int -> Int -> Vec
+gamma i j k =
+ let a = phi!(mod (iabs k) 12)
+     b = phi!(mod (iabs (j+a)) 12)
+     c = phi!(mod (iabs (i+b)) 12)
+ in grad!c
+
+knot :: Int -> Int -> Int -> Vec -> Flt
+knot i j k v =
+ let Vec x y z = v
+ in (omega x) * (omega y) * (omega z) * (vdot (gamma i j k) v)
+
+intGamma :: Int -> Int -> Int
+intGamma i j =
+ let a = phi!(mod (iabs j) 16)
+     b = phi!(mod (iabs (i+a)) 16)
+ in b
+
+turbulence :: Vec -> Int -> Flt
+turbulence p 1 = fabs(noise(p))
+turbulence p n =
+ let newp = vscale p 0.5
+     t = fabs (noise p)
+ in t + (0.5 * (turbulence newp (n-1)))
+
+noise :: Vec -> Flt 
+noise (Vec x y z) =
+ let i = floor x
+     j = floor y
+     k = floor z
+     u = x-(fromIntegral i)
+     v = y-(fromIntegral j)
+     w = z-(fromIntegral k)
+ in knot i j k             (Vec u v w) +
+    knot (i+1) j k         (Vec (u-1) v w) +
+    knot i (j+1) k         (Vec u (v-1) w) +
+    knot i j (k+1)         (Vec u v (w-1)) +
+    knot (i+1) (j+1) k     (Vec (u-1) (v-1) w) +
+    knot (i+1) j (k+1)     (Vec (u-1) v (w-1)) +
+    knot i (j+1) (k+1)     (Vec u (v-1) (w-1)) +
+    knot (i+1) (j+1) (k+1) (Vec (u-1) (v-1) (w-1))
+
+perlin :: Vec -> Flt
+perlin v =
+ let p = ((noise v)+1)*0.5
+ in if p > 1 
+    then error $ "perlin noise error, 1 < " ++ (show p)
+    else if p < 0 
+         then error $ "perlin noise error, 0 > " ++ (show p)
+         else p
+
+--untested
+perlin_turb :: Vec -> Int -> Flt
+perlin_turb v l =
+ let p = turbulence v l
+ in if p > 1 
+    then error $ "perlin turbulence error, 1 < " ++ (show p)
+    else if p < 0 
+         then error $ "perlin turbulence error, 0 > " ++ (show p)
+         else p
diff --git a/Data/Glome/Vec.hs b/Data/Glome/Vec.hs
new file mode 100644
--- /dev/null
+++ b/Data/Glome/Vec.hs
@@ -0,0 +1,738 @@
+{-# OPTIONS_GHC -fexcess-precision #-}
+{-# OPTIONS_GHC -funbox-strict-fields #-}
+{-# LANGUAGE BangPatterns #-}
+
+module Data.Glome.Vec where
+
+-- | Performance is pretty similar with Floats or Doubles.
+-- Todo: make separate Float and Double instances of this library.
+type Flt = Double
+
+-- maybe this is defined somewhere?
+infinity :: Flt
+--infinity = 1.0 / 0.0
+infinity = 1000000.0
+
+-- | Convert from degrees to native angle format (radians).
+deg :: Flt -> Flt
+deg !x = (x*3.1415926535897)/180
+
+-- | Convert from radians to native format (noop).
+rad :: Flt -> Flt
+rad !x = x
+
+-- | Convert from rotations to native format.  (rot 1 == deg 360)
+rot :: Flt -> Flt
+rot !x = x*3.1415926535897*2
+
+-- | Trig with degrees instead of radians.
+dcos :: Flt -> Flt
+dcos d = cos $ deg d
+
+-- | Force a value to be within a range.  Usage: clamp min x max
+clamp :: Flt -> Flt -> Flt -> Flt
+clamp !min !x !max
+ | x < min = min
+ | x > max = max
+ | otherwise = x
+
+-- | Tuning parameter.
+delta = 0.0001 :: Flt
+
+-- | Non-polymorphic fmin; this speeds
+-- things up in ocaml, not sure about haskell.
+fmin :: Flt -> Flt -> Flt
+fmin !a !b = if a > b then b else a
+
+-- | Non-polymorphic fmax.
+fmax :: Flt -> Flt -> Flt
+fmax !a !b = if a > b then a else b
+
+-- | Non-polymorphic min of 3 values.
+fmin3 :: Flt -> Flt -> Flt -> Flt
+fmin3 !a !b !c = if a > b 
+                 then if b > c 
+                      then c
+                      else b
+                 else if a > c
+                      then c
+                      else a
+
+-- | Non-polymorphic max of 3 values.
+fmax3 :: Flt -> Flt -> Flt -> Flt
+fmax3 !a !b !c = if a > b
+                 then if a > c
+                      then a
+                      else c
+                 else if b > c
+                      then b
+                      else c
+
+-- | Min of 4 values.
+fmin4 :: Flt -> Flt -> Flt -> Flt -> Flt
+fmin4 !a !b !c !d = fmin (fmin a b) (fmin c d)
+
+-- | Max of 4 values.
+fmax4 :: Flt -> Flt -> Flt -> Flt -> Flt
+fmax4 !a !b !c !d = fmax (fmax a b) (fmax c d)
+
+-- | Non-polymorphic absolute value.
+fabs :: Flt -> Flt
+fabs !a = 
+ if a < 0 then (-a) else a
+
+-- | Non-polymorphic integer absolute value.
+iabs :: Int -> Int
+iabs !a =
+ if a < 0 then (-a) else a
+
+-- | Force user to use fabs or iabs, for performance reasons.  Not sure if
+-- this really helps, though.
+abs a = error "use non-polymorphic version, fabs"
+
+-- | Approximate equality for Flt.  True if a and b are "almost" equal.
+-- The (abs $ a-b) test doesn't work if
+-- a and b are large.
+about_equal :: Flt -> Flt -> Bool
+about_equal !a !b =
+ if a > 1 
+ then
+  fabs (1 - (a/b)) < (delta*10) 
+ else
+  (fabs $ a - b) < (delta*10)
+
+-- | 3d type represented as a record of unboxed floats.
+data Vec = Vec !Flt !Flt !Flt deriving Show
+
+-- | A Ray is made up of an origin and direction Vec.
+data Ray = Ray {origin, dir :: !Vec} deriving Show
+--data Plane = Plane {norm :: !Vec, offset :: !Flt} deriving Show
+
+-- | Vec constructor.
+vec :: Flt -> Flt -> Flt -> Vec
+vec !x !y !z = (Vec x y z)
+
+-- | Zero Vec.
+vzero :: Vec
+vzero = Vec 0.0 0.0 0.0
+
+-- | For when we need a unit vector, but we 
+-- don't care where it points.
+vunit :: Vec
+vunit = vx
+
+-- | Unit X vector.
+vx :: Vec
+vx  = Vec 1 0 0
+
+-- | Unit y vector.
+vy :: Vec
+vy  = Vec 0 1 0
+
+-- | Unit z vector.
+vz :: Vec
+vz  = Vec 0 0 1
+
+-- | Negative x vector.
+nvx :: Vec
+nvx = Vec (-1) 0 0
+
+-- | Negative y vector.
+nvy :: Vec
+nvy = Vec 0 (-1) 0
+
+-- | Negative z vector.
+nvz :: Vec
+nvz = Vec 0 0 (-1)
+
+-- Extract x coordinate.
+x :: Vec -> Flt
+x (Vec x_ _ _) = x_
+
+-- Extract y coordinate.
+y :: Vec -> Flt
+y (Vec _ y_ _) = y_
+
+-- Extract z coordinate.
+z :: Vec -> Flt
+z (Vec _ _ z_) = z_
+
+-- | Access the Vec as if it were an array indexed from 0..2.
+-- Note: this actually accounts for a noticeable amount of cpu 
+-- time in the Glome ray tracer.
+va :: Vec -> Int -> Flt
+va !(Vec x y z) !n = 
+ case n of
+  0 -> x
+  1 -> y
+  2 -> z
+
+-- | Create a new Vec with the Nth field overwritten by new value.
+-- I could have used record update syntax.
+vset :: Vec -> Int -> Flt -> Vec
+vset !(Vec x y z) !i !f =
+ case i of
+  0 -> Vec f y z
+  1 -> Vec x f z
+  2 -> Vec x y f
+
+-- | Dot product of 2 vectors.  We use this all the time.  Dot product of 2
+-- normal vectors is the cosine of the angle between them.
+vdot :: Vec -> Vec -> Flt
+vdot !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
+ (x1*x2)+(y1*y2)+(z1*z2)
+
+-- | Cross product of 2 vectors.  Produces a vector perpendicular 
+-- to the given vectors.  We use this for things like making the forward,
+-- up, and right camera vectors orthogonal.  If the input vectors are
+-- normalized, the output vector will be as well.
+vcross :: Vec -> Vec -> Vec
+vcross !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
+ Vec 
+  ((y1 * z2) - (z1 * y2))
+  ((z1 * x2) - (x1 * z2))
+  ((x1 * y2) - (y1 * x2))
+
+-- | Apply a unary Flt operator to each field of the Vec.
+vmap :: (Flt -> Flt) -> Vec -> Vec
+vmap f !v1 = 
+ Vec (f (x v1)) (f (y v1)) (f (z v1))
+
+-- | Apply a binary Flt operator to pairs of fields from 2 Vecs.
+vmap2 :: (Flt -> Flt -> Flt) -> Vec -> Vec -> Vec
+vmap2 f !v1 !v2 =
+ Vec (f (x v1) (x v2)) 
+     (f (y v1) (y v2)) 
+     (f (z v1) (z v2))
+
+-- | Reverse the direction of a Vec.
+vinvert :: Vec -> Vec
+vinvert !(Vec x1 y1 z1) =
+ Vec (-x1) (-y1) (-z1)
+
+-- | Get the length of a Vec squared.  We use this to avoid a slow sqrt. 
+vlensqr :: Vec -> Flt
+vlensqr !v1 = vdot v1 v1
+
+-- | Get the length of a Vec.  This is expensive because sqrt is slow.
+vlen :: Vec -> Flt
+vlen !v1 = sqrt (vdot v1 v1)
+
+-- | Add 2 vectors.
+vadd :: Vec -> Vec -> Vec
+vadd !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
+ Vec (x1 + x2)
+     (y1 + y2)
+     (z1 + z2)
+
+-- | Add 3 vectors.
+vadd3 :: Vec -> Vec -> Vec -> Vec
+vadd3 !(Vec x1 y1 z1) !(Vec x2 y2 z2) !(Vec x3 y3 z3) =
+    Vec (x1 + x2 + x3)
+        (y1 + y2 + y3)
+        (z1 + z2 + z3)
+
+-- | Subtract vectors.  "vsub b a" is the vector from a to b.
+vsub :: Vec -> Vec -> Vec
+vsub !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
+ Vec (x1 - x2)
+     (y1 - y2)
+     (z1 - z2)
+
+-- | Multiply corresponding fields.  Rarely useful.
+vmul :: Vec -> Vec -> Vec
+vmul !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
+ Vec (x1 * x2)
+     (y1 * y2)
+     (z1 * z2)
+
+-- | Add a value to all the fields of a Vec.  Useful, for instance, to get
+-- one corner of the bounding box around a sphere.
+vinc :: Vec -> Flt -> Vec
+vinc !(Vec x y z) !n =
+ Vec (x + n)
+     (y + n)
+     (z + n)
+
+-- | Subtract a value from all fields of a Vec.
+vdec :: Vec -> Flt -> Vec
+vdec !(Vec x y z) !n =
+ Vec (x - n)
+     (y - n)
+     (z - n)
+
+-- | Get the maximum of all corresponding fields between 2 Vecs.
+vmax :: Vec -> Vec -> Vec
+vmax !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
+ Vec (fmax x1 x2)
+     (fmax y1 y2)
+     (fmax z1 z2)
+
+-- | Get the minimum of all corresponding fields between 2 Vecs.
+vmin :: Vec -> Vec -> Vec
+vmin !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
+ Vec (fmin x1 x2)
+     (fmin y1 y2)
+     (fmin z1 z2)
+
+-- | Return the largest axis.  Often used with "va".
+vmaxaxis :: Vec -> Int
+vmaxaxis !(Vec x y z) =
+ if (x > y) 
+ then if (x > z) 
+      then 0
+      else 2
+ else if (y > z) 
+      then 1
+      else 2
+
+-- | Scale a Vec by some value.
+vscale :: Vec -> Flt -> Vec
+vscale !(Vec x y z) !fac =
+ Vec (x * fac)
+     (y * fac)
+     (z * fac)
+
+-- | Take the first Vec, and add to it the second Vec scaled by some amount.
+-- This is used quite a lot in Glome.
+vscaleadd :: Vec -> Vec -> Flt -> Vec
+vscaleadd !(Vec x1 y1 z1) !(Vec x2 y2 z2) fac =
+ Vec (x1 + (x2 * fac))
+     (y1 + (y2 * fac))
+     (z1 + (z2 * fac))
+            
+-- | Make the length of a Vec just a little shorter.
+vnudge :: Vec -> Vec
+vnudge x = vscale x (1-delta)
+
+-- | Normalize a vector.  Division is expensive, so we compute the reciprocol 
+-- of the length and multiply by that.  The sqrt is also expensive.
+vnorm :: Vec -> Vec
+vnorm !(Vec x1 y1 z1) = 
+ let !invlen = 1.0 / (sqrt ((x1*x1)+(y1*y1)+(z1*z1))) in
+ Vec (x1*invlen) (y1*invlen) (z1*invlen)
+
+-- | Throw an exception if a vector hasn't been normalized.
+assert_norm :: Vec -> Vec
+assert_norm v =
+ let l = vdot v v
+ in if l > (1+delta) 
+    then error $ "vector too long" ++ (show v)
+    else if l < (1-delta)
+         then error $ "vector too short: " ++ (show v)
+         else v
+
+-- | Get the victor bisecting two other vectors (which ought to be the same
+-- length).
+bisect :: Vec -> Vec -> Vec
+bisect !v1 !v2 = vnorm (vadd v1 v2)
+
+-- | Distance between 2 vectors.
+vdist :: Vec -> Vec -> Flt
+vdist v1 v2 = 
+ let d = vsub v2 v1 in vlen d
+
+-- | Reflect a vector "v" off of a surface with normal "norm".
+reflect :: Vec -> Vec -> Vec
+reflect !v !norm =
+  -- vadd v $ vscale norm $ (-2) * (vdot v norm)
+  vscaleadd v norm $ (-2) * (vdot v norm)
+
+-- | Reciprocol of all fields of a Vec.
+vrcp :: Vec -> Vec
+vrcp !(Vec x y z) =
+ Vec (1/x) (1/y) (1/z)
+
+-- | Test Vecs for approximate equality
+veq :: Vec -> Vec -> Bool
+veq !(Vec ax ay az) !(Vec bx by bz) =
+ (about_equal ax bx) && (about_equal ay by) && (about_equal az bz)
+
+-- | Test Vecs for matching sign on all fields.  Returns false if any value is
+-- zero.  Used by packet tracing.
+veqsign :: Vec -> Vec -> Bool
+veqsign !(Vec ax ay az) !(Vec bx by bz) =
+ ax*bx > 0 && ay*by > 0 && az*bz > 0
+
+-- | Translate a ray's origin in ray's direction by d amount.
+ray_move :: Ray -> Flt -> Ray
+ray_move !(Ray orig dir) !d =
+ (Ray (vscaleadd orig dir d) dir)
+
+-- | Find a pair of orthogonal vectors to the one given.
+orth :: Vec -> (Vec,Vec)
+orth v1 =
+ if about_equal (vdot v1 v1) 1
+ then
+  let x = (Vec 1 0 0)
+      y = (Vec 0 1 0)
+      dvx = vdot v1 x
+      v2 = if dvx < 0.8 && dvx > (-0.8) -- don't want to cross with a
+           then vnorm $ vcross v1 x     -- vector that's too similar
+           else vnorm $ vcross v1 y
+      v3 = vcross v1 v2
+  in (v2,v3)
+ else error $ "orth: unnormalized vector" ++ (show v1)
+
+-- | Intersect a ray with a plane 
+-- defined by a point "p" and a normal "norm".
+-- (Ray does not need to be normalized.)
+plane_int :: Ray -> Vec -> Vec -> Vec
+plane_int !(Ray orig dir) !p !norm =
+ let newo = vsub orig p
+     dist = -(vdot norm newo) / (vdot norm dir)
+ in vscaleadd orig dir dist
+
+-- | Find the distance along a ray until it intersects with a plane defined
+-- by a point "p" and normal "norm".
+plane_int_dist :: Ray -> Vec -> Vec -> Flt
+plane_int_dist !(Ray orig dir) !p !norm =
+ let newo = vsub orig p
+ in -(vdot norm newo) / (vdot norm dir)
+
+-- find intersection with plane
+-- from graphics gems -- an efficient ray-polygon intersection
+-- it seems that the ray need not be normalized
+-- let plane_intersect ray (n,d) =
+--  let t = -.((d +. (vdot n ray.origin)) /. (vdot n ray.dir))
+--  in vadd ray.origin (vscale ray.dir t)
+
+
+-- TRANSFORMATIONS --
+
+-- | 3x4 Transformation matrix.  These are described in most graphics texts.
+data Matrix = Matrix !Flt !Flt !Flt !Flt  
+                     !Flt !Flt !Flt !Flt  
+                     !Flt !Flt !Flt !Flt deriving Show
+
+-- | A transformation.  Inverting a matrix is expensive, so we keep a forward
+-- transformation matrix and a reverse transformation matrix.
+-- Note: This can be made a little faster if the matricies are non-strict.
+data Xfm = Xfm Matrix Matrix deriving Show
+
+-- | Identity matrix.  Transforming a vector by this matrix does nothing.
+ident_matrix :: Matrix
+ident_matrix = (Matrix 1 0 0 0  0 1 0 0  0 0 1 0)
+
+-- | Identity transformation.
+ident_xfm :: Xfm
+ident_xfm = Xfm ident_matrix ident_matrix
+
+-- | Multiply two matricies.  This is unrolled for efficiency, and it's also
+-- a little bit easier (in my opinion) to see what's going on.
+mat_mult :: Matrix -> Matrix -> Matrix
+mat_mult (Matrix a00 a01 a02 a03  a10 a11 a12 a13  a20 a21 a22 a23)
+         (Matrix b00 b01 b02 b03  b10 b11 b12 b13  b20 b21 b22 b23) =
+ Matrix
+   (a00*b00 + a01*b10 + a02*b20)
+   (a00*b01 + a01*b11 + a02*b21)
+   (a00*b02 + a01*b12 + a02*b22)
+   (a00*b03 + a01*b13 + a02*b23 + a03)
+
+   (a10*b00 + a11*b10 + a12*b20)
+   (a10*b01 + a11*b11 + a12*b21)
+   (a10*b02 + a11*b12 + a12*b22)
+   (a10*b03 + a11*b13 + a12*b23 + a13)
+
+   (a20*b00 + a21*b10 + a22*b20)
+   (a20*b01 + a21*b11 + a22*b21)
+   (a20*b02 + a21*b12 + a22*b22)
+   (a20*b03 + a21*b13 + a22*b23 + a23)
+
+-- | Multiply two tranformations.  This just multiplies the forward and 
+-- reverse transformations.
+xfm_mult :: Xfm -> Xfm -> Xfm
+xfm_mult (Xfm a inva) (Xfm b invb) =
+ Xfm (mat_mult a b) (mat_mult invb inva)
+
+-- TRANSFORM UTILITY FUNCTIONS --
+
+-- | There is a seemingly-magical property of transformation matricies, that
+-- we can combine the effects of any number of transformations into a single
+-- transformation just by multiplying them together in reverse order.  For 
+-- instance, we could move a point, then rotate it about the origin by some 
+-- angle around some vector, then move it again, and this can all be done by 
+-- a single transformation.
+-- This function combines transformations in this way, though it reverses the
+-- list first so the transformations take effect in their expected order.
+compose :: [Xfm] -> Xfm
+compose xfms = check_xfm $ foldr xfm_mult ident_xfm (reverse xfms)
+
+-- | Make sure a transformation is valid.  Multipy the forward and reverse
+-- matrix and verify that the result is the identity matrix.
+check_xfm :: Xfm -> Xfm
+check_xfm (Xfm m i) = 
+ let (Matrix m00 m01 m02 m03  
+             m10 m11 m12 m13  
+             m20 m21 m22 m23) = mat_mult m i
+     ae = about_equal
+ in
+  if ae m00 1 && ae m01 0 && ae m02 0 && ae m03 0 &&
+     ae m10 0 && ae m11 1 && ae m12 0 && ae m13 0 &&
+     ae m20 0 && ae m21 0 && ae m22 1 && ae m23 0
+  then (Xfm m i)
+  else error $ "corrupt matrix " ++ (show (Xfm m i)) ++ "\n" ++ (show (mat_mult m i)) 
+
+-- | Complex transformations: Rotate point (or vector) "pt" about ray by 
+-- angle c.  The angle is in radians,
+-- but using the angle conversion routines "deg", "rad" and "rot" is 
+-- recommended.
+vrotate :: Vec -> Ray -> Flt -> Vec
+vrotate pt (Ray orig axis_) angle =
+ let axis = assert_norm axis_
+     transform = compose [ translate (vinvert orig)
+                         , rotate axis angle
+                         , translate orig
+                         ]
+     new_pt = xfm_point transform pt
+ in if about_equal (vlen (vsub orig pt)) (vlen (vsub orig new_pt))
+    then new_pt
+    else error $ "something is wrong with vrotate" ++ 
+                 (show $ vlen (vsub orig pt)) ++ " " ++ 
+                 (show $ vlen (vsub orig new_pt))
+
+
+-- TRANSFORM APPLICATION --
+-- these need to be fast
+
+-- | Transform a point.  The point is treated as (x y z 1).
+xfm_point :: Xfm -> Vec -> Vec
+xfm_point !(Xfm (Matrix m00 m01 m02 m03  
+                        m10 m11 m12 m13  
+                        m20 m21 m22 m23) inv) 
+          !(Vec x y z) =
+ Vec (m00*x + m01*y + m02*z + m03)
+     (m10*x + m11*y + m12*z + m13)
+     (m20*x + m21*y + m22*z + m23)
+
+-- | Inverse transform a point.
+invxfm_point :: Xfm -> Vec -> Vec
+invxfm_point !(Xfm fwd (Matrix i00 i01 i02 i03  
+                               i10 i11 i12 i13  
+                               i20 i21 i22 i23)) 
+             !(Vec x y z) =
+  Vec (i00*x + i01*y + i02*z + i03)
+      (i10*x + i11*y + i12*z + i13)
+      (i20*x + i21*y + i22*z + i23)
+
+-- | Transform a vector.  The vector is treated as (x y z 0).
+xfm_vec :: Xfm -> Vec -> Vec
+xfm_vec !(Xfm (Matrix m00 m01 m02 m03  
+                      m10 m11 m12 m13  
+                      m20 m21 m22 m23) inv) 
+        !(Vec x y z) =
+ Vec (m00*x + m01*y + m02*z)
+     (m10*x + m11*y + m12*z)
+     (m20*x + m21*y + m22*z)
+
+-- | Inverse transform a vector.
+invxfm_vec :: Xfm -> Vec -> Vec
+invxfm_vec !(Xfm fwd (Matrix i00 i01 i02 i03  
+                             i10 i11 i12 i13  
+                             i20 i21 i22 i23)) 
+           !(Vec x y z) =
+  Vec (i00*x + i01*y + i02*z)
+      (i10*x + i11*y + i12*z)
+      (i20*x + i21*y + i22*z)
+
+-- | Inverse transform a normal.  This one is tricky: we need to transform 
+-- by the inverse transpose.
+invxfm_norm :: Xfm -> Vec -> Vec
+invxfm_norm !(Xfm fwd (Matrix i00 i01 i02 i03  
+                              i10 i11 i12 i13  
+                              i20 i21 i22 i23)) 
+            !(Vec x y z) =
+ Vec (i00*x + i10*y + i20*z)
+     (i01*x + i11*y + i21*z)
+     (i02*x + i12*y + i22*z)
+
+-- | Transform a Ray.
+xfm_ray :: Xfm -> Ray -> Ray
+xfm_ray !xfm !(Ray orig dir) =
+ Ray (xfm_point xfm orig) (vnorm (xfm_vec xfm dir))
+
+-- | Inverse transform a Ray.
+invxfm_ray :: Xfm -> Ray -> Ray
+invxfm_ray !xfm !(Ray orig dir) =
+ Ray (invxfm_point xfm orig) (vnorm (invxfm_vec xfm dir))
+
+-- BASIC TRANSFORMS --
+-- | Basic transforms: move by some displacement vector.
+translate :: Vec -> Xfm
+translate (Vec x y z) =
+ check_xfm $ Xfm (Matrix 1 0 0   x   0 1 0   y   0 0 1   z) 
+                 (Matrix 1 0 0 (-x)  0 1 0 (-y)  0 0 1 (-z))
+
+-- | Basic transforms: stretch along the three axes, by the amount
+-- in the given vector.  (If x==y==z, then it's uniform scaling.)
+scale :: Vec -> Xfm
+scale (Vec x y z) =
+ check_xfm $ Xfm (Matrix   x  0 0 0  0   y  0 0  0 0   z  0)
+                (Matrix (1/x) 0 0 0  0 (1/y) 0 0  0 0 (1/z) 0)
+
+-- | Basic transforms: rotate about a given axis by some angle.
+rotate :: Vec -> Flt -> Xfm
+rotate v@(Vec x y z) angle =
+ if not $ (vlen v) `about_equal` 1
+ then error $ "please use a normalized vector for rotation: " ++ (show (vlen v))
+ else 
+  let s = sin angle
+      c = cos angle 
+
+      m00 = ((x*x)+((1-(x*x))*c)) 
+      m01 = (((x*y)*(1-c))-(z*s)) 
+      m02 = ((x*z*(1-c))+(y*s))
+
+      m10 = (((x*y)*(1-c))+(z*s))
+      m11 = ((y*y)+((1-(y*y))*c))
+      m12 = ((y*z*(1-c))-(x*s))
+
+      m20 = ((x*z*(1-c))-(y*s))
+      m21 = ((y*z*(1-c))+(x*s))
+      m22 = ((z*z)+((1-(z*z))*c))
+  in
+  check_xfm $ Xfm (Matrix m00 m01 m02 0  m10 m11 m12 0  m20 m21 m22 0)
+                  (Matrix m00 m10 m20 0  m01 m11 m21 0  m02 m12 m22 0)
+
+-- | Basic transforms: Convert coordinate system from canonical xyz 
+-- coordinates to uvw coordinates.
+xyz_to_uvw :: Vec -> Vec -> Vec -> Xfm
+xyz_to_uvw u v w =
+ let Vec ux uy uz = u
+     Vec vx vy vz = v
+     Vec wx wy wz = w
+ in if (vdot u u) `about_equal` 1
+    then
+     if (vdot v v) `about_equal` 1
+     then
+      if (vdot w w) `about_equal` 1
+      then 
+       if ((vdot u v) `about_equal` 0) && 
+          ((vdot u w) `about_equal` 0) && 
+          ((vdot v w) `about_equal` 0)
+       then
+        check_xfm $ Xfm (Matrix ux vx wx 0  uy vy wy 0  uz vz wz 0)
+                        (Matrix ux uy uz 0  vx vy vz 0  wx wy wz 0)
+       else error  "vectors aren't orthogonal"
+      else error $ "unnormalized w " ++ (show w)
+     else error $ "unnormalized v " ++ (show v)
+    else error $ "unnormalized u " ++ (show u)
+
+-- | Basic transforms: Convert from uvw coordinates back to normal xyz 
+-- coordinates.
+uvw_to_xyz :: Vec -> Vec -> Vec -> Xfm
+uvw_to_xyz (Vec ux uy uz) (Vec vx vy vz) (Vec wx wy wz) =
+ check_xfm $ Xfm (Matrix ux uy uz 0  vx vy vz 0  wx wy wz 0)
+                 (Matrix ux vx wx 0  uy vy wy 0  uz vz wz 0)
+
+
+
+-- TRIANGLE UTILITY FUNCTIONS --
+
+-- | Given a side, angle, and side of a triangle, produce the length of the 
+-- opposite side.
+sas2s :: Flt -> Flt -> Flt -> Flt
+sas2s s1 a s2 =
+  sqrt (((s1 * s1) + (s2 * s2)) - ((2 * s1 * s2 * (dcos a))))
+
+
+
+-- BOUNDING BOXES --
+-- | Axis-aligned Bounding Box (AABB), defined by opposite corners.  P1 is the
+-- min values, p2 has the max values.
+data Bbox = Bbox {p1 :: !Vec, p2 :: !Vec} deriving Show
+
+-- | A near-far pair of distances.  Basically just a tuple.
+data Interval = Interval !Flt !Flt deriving Show -- used instead of a tuple
+
+-- | Bounding box that encloses two bounding boxes.
+bbjoin :: Bbox -> Bbox -> Bbox
+bbjoin (Bbox p1a p2a) (Bbox p1b p2b) =
+ (Bbox (vmin p1a p1b) (vmax p2a p2b))
+
+-- | Find the overlap of two bounding boxes.
+bboverlap :: Bbox -> Bbox -> Bbox
+bboverlap (Bbox p1a p2a) (Bbox p1b p2b) =
+ (Bbox (vmax p1a p1b) (vmin p2a p2b))
+
+-- | Test if a Vec is inside the bounding box.
+bbinside :: Bbox -> Vec -> Bool
+bbinside (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) (Vec x y z) =
+ p1x <= x && x <= p2x && p1y <= y && y <= p2y && p1z <= z && z <= p2z
+
+-- | Split a bounding box into two, given an axis and offset.  Throw exception
+-- if the offset isn't inside the bounding box.
+bbsplit :: Bbox -> Int -> Flt -> (Bbox,Bbox)
+bbsplit (Bbox p1 p2) axis offset =
+ if (offset < (va p1 axis)) || (offset > (va p2 axis))
+ then error "degenerate bounding box split"
+ else ((Bbox p1 (vset p2 axis offset)),
+       (Bbox (vset p1 axis offset) p2))
+
+-- | Generate a minimum bounding box that encloses a list of points.
+bbpts :: [Vec] -> Bbox
+bbpts [] = empty_bbox
+bbpts ((Vec x y z):[]) =
+ Bbox (Vec (x-delta) (y-delta) (z-delta)) 
+      (Vec (x+delta) (y+delta) (z+delta))
+
+bbpts ((Vec x y z):pts) =
+ let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bbpts pts
+     minx = fmin (x-delta) p1x
+     miny = fmin (y-delta) p1y
+     minz = fmin (z-delta) p1z
+     maxx = fmax (x+delta) p2x
+     maxy = fmax (y+delta) p2y
+     maxz = fmax (z+delta) p2z in
+ Bbox (Vec minx miny minz) (Vec maxx maxy maxz)
+
+-- | Surface area of a bounding box.  Useful for cost heuristics when attempting
+-- to build optimal bounding box heirarchies.  Undefined for degenerate bounding
+-- boxes.
+bbsa :: Bbox -> Flt
+bbsa (Bbox p1 p2) =
+ let Vec dx dy dz = vsub p2 p1 
+ in dx*dy + dx*dz + dy*dz
+
+-- | Volume of a bounding box.  Undefined for degenerate bounding boxes.
+bbvol :: Bbox -> Flt
+bbvol (Bbox p1 p2) =
+ let (Vec dx dy dz) = vsub p2 p1
+ in dx*dy*dz
+
+-- | Degenerate bounding box that contains an empty volume.
+empty_bbox :: Bbox
+empty_bbox = 
+ Bbox (Vec infinity infinity infinity) 
+      (Vec (-infinity) (-infinity) (-infinity))
+
+-- | "Infinite" bounding box.
+everything_bbox :: Bbox
+everything_bbox =
+ Bbox (Vec (-infinity) (-infinity) (-infinity))
+      (Vec infinity infinity infinity)
+
+-- | Find a ray's entrance and exit from a bounding 
+-- box.  If last entrance is before the first exit,
+-- we hit.  Otherwise, we miss. (It's up to the 
+-- caller to figure that out.)
+
+bbclip :: Ray -> Bbox -> Interval
+bbclip (Ray (Vec ox oy oz) (Vec dx dy dz)) 
+       (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) =
+ let dxrcp = 1/dx
+     dyrcp = 1/dy
+     dzrcp = 1/dz
+     Interval inx outx = if dx > 0 
+                         then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp)
+                         else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp)
+     Interval iny outy = if dy > 0
+                         then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp)
+                         else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp)
+     Interval inz outz = if dz > 0
+                         then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp)
+                         else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp)
+ in
+   Interval (fmax3 inx iny inz) (fmin3 outx outy outz)
+
+
diff --git a/Data/GlomeTexture.hs b/Data/GlomeTexture.hs
deleted file mode 100644
--- a/Data/GlomeTexture.hs
+++ /dev/null
@@ -1,126 +0,0 @@
-{-
-module SolidTexture (square_wave, triangle_wave, sine_wave, 
-                     stripe, noise, turbulence
-                    ) where -}
-
-module Data.GlomeTexture where
-import Data.GlomeVec
-import Data.Array.IArray
-
--- INTERPOLATION FUNCTIONS --
-square_wave :: Flt -> Flt
-square_wave x =
- let offset = x - (fromIntegral (floor x))
- in if offset < 0.5 then 0 else 1
-
-triangle_wave :: Flt -> Flt
-triangle_wave x =
- let offset = x - (fromIntegral (floor x))
- in if offset < 0.5 
-    then (offset*2)
-    else (2-(offset*2))
-
-sine_wave :: Flt -> Flt
-sine_wave x = (sin (x*2*pi))*0.5 + 0.5
-
-
-lump_wave :: Flt -> Flt
-lump_wave x = 1 - x*x*x
-
--- SCALAR TEXTURE FUNCTIONS --
-
--- These are simple solid texture functions that take a 
--- point as argument and return a number 0 < n < 1
-
-stripe :: Vec -> (Flt -> Flt) -> (Vec -> Flt)
-stripe axis interp =
- let len = vlen axis 
- in
-  (\pos -> let offset = vdot pos axis 
-           in interp offset)
-
-
--- PERLIN NOISE --
-
--- (-6 t^5 + 15 t^4 - 10t^3 +1)
--- "realistic ray tracing 2nd edition" inconsistent 
--- on whether it should be t^5 or t^6,
--- but t^5 works and t^6 doesn't.
-omega :: Flt -> Flt
-omega t_ = 
- let t     = fabs t_
-     tsqr  = t*t
-     tcube = tsqr*t
- in (-6)*tcube*tsqr + 15*tcube*t - 10*tcube + 1
-
--- questionably random
-phi :: Array Int Int
-phi = listArray (0,11) [3,0,2,7,4,1,5,11,8,10,9,6]
-
-grad :: Array Int Vec
-grad = listArray (0,11) 
-         $ filter (\x -> let l = vlen x in l < 1.5 && l > 1.1) 
-                  [Vec x y z | x <- [(-1),0,1],
-                               y <- [(-1),0,1],
-                               z <- [(-1),0,1]] 
-
-gamma :: Int -> Int -> Int -> Vec
-gamma i j k =
- let a = phi!(mod (iabs k) 12)
-     b = phi!(mod (iabs (j+a)) 12)
-     c = phi!(mod (iabs (i+b)) 12)
- in grad!c
-
-knot :: Int -> Int -> Int -> Vec -> Flt
-knot i j k v =
- let Vec x y z = v
- in (omega x) * (omega y) * (omega z) * (vdot (gamma i j k) v)
-
-intGamma :: Int -> Int -> Int
-intGamma i j =
- let a = phi!(mod (iabs j) 16)
-     b = phi!(mod (iabs (i+a)) 16)
- in b
-
-turbulence :: Vec -> Int -> Flt
-turbulence p 1 = fabs(noise(p))
-turbulence p n =
- let newp = vscale p 0.5
-     t = fabs (noise p)
- in t + (0.5 * (turbulence newp (n-1)))
-
-noise :: Vec -> Flt 
-noise (Vec x y z) =
- let i = floor x
-     j = floor y
-     k = floor z
-     u = x-(fromIntegral i)
-     v = y-(fromIntegral j)
-     w = z-(fromIntegral k)
- in knot i j k             (Vec u v w) +
-    knot (i+1) j k         (Vec (u-1) v w) +
-    knot i (j+1) k         (Vec u (v-1) w) +
-    knot i j (k+1)         (Vec u v (w-1)) +
-    knot (i+1) (j+1) k     (Vec (u-1) (v-1) w) +
-    knot (i+1) j (k+1)     (Vec (u-1) v (w-1)) +
-    knot i (j+1) (k+1)     (Vec u (v-1) (w-1)) +
-    knot (i+1) (j+1) (k+1) (Vec (u-1) (v-1) (w-1))
-
-perlin :: Vec -> Flt
-perlin v =
- let p = ((noise v)+1)*0.5
- in if p > 1 
-    then error $ "perlin noise error, 1 < " ++ (show p)
-    else if p < 0 
-         then error $ "perlin noise error, 0 > " ++ (show p)
-         else p
-
---untested
-perlin_turb :: Vec -> Int -> Flt
-perlin_turb v l =
- let p = turbulence v l
- in if p > 1 
-    then error $ "perlin turbulence error, 1 < " ++ (show p)
-    else if p < 0 
-         then error $ "perlin turbulence error, 0 > " ++ (show p)
-         else p
diff --git a/Data/GlomeVec.hs b/Data/GlomeVec.hs
deleted file mode 100644
--- a/Data/GlomeVec.hs
+++ /dev/null
@@ -1,738 +0,0 @@
-{-# OPTIONS_GHC -fexcess-precision #-}
-{-# OPTIONS_GHC -funbox-strict-fields #-}
-{-# LANGUAGE BangPatterns #-}
-
-module Data.GlomeVec where
-
--- | Performance is pretty similar with Floats or Doubles.
--- Todo: make separate Float and Double instances of this library.
-type Flt = Double
-
--- maybe this is defined somewhere?
-infinity :: Flt
---infinity = 1.0 / 0.0
-infinity = 1000000.0
-
--- | Convert from degrees to native angle format (radians).
-deg :: Flt -> Flt
-deg !x = (x*3.1415926535897)/180
-
--- | Convert from radians to native format (noop).
-rad :: Flt -> Flt
-rad !x = x
-
--- | Convert from rotations to native format.  (rot 1 == deg 360)
-rot :: Flt -> Flt
-rot !x = x*3.1415926535897*2
-
--- | Trig with degrees instead of radians.
-dcos :: Flt -> Flt
-dcos d = cos $ deg d
-
--- | Force a value to be within a range.  Usage: clamp min x max
-clamp :: Flt -> Flt -> Flt -> Flt
-clamp !min !x !max
- | x < min = min
- | x > max = max
- | otherwise = x
-
--- | Tuning parameter.
-delta = 0.0001 :: Flt
-
--- | Non-polymorphic fmin; this speeds
--- things up in ocaml, not sure about haskell.
-fmin :: Flt -> Flt -> Flt
-fmin !a !b = if a > b then b else a
-
--- | Non-polymorphic fmax.
-fmax :: Flt -> Flt -> Flt
-fmax !a !b = if a > b then a else b
-
--- | Non-polymorphic min of 3 values.
-fmin3 :: Flt -> Flt -> Flt -> Flt
-fmin3 !a !b !c = if a > b 
-                 then if b > c 
-                      then c
-                      else b
-                 else if a > c
-                      then c
-                      else a
-
--- | Non-polymorphic max of 3 values.
-fmax3 :: Flt -> Flt -> Flt -> Flt
-fmax3 !a !b !c = if a > b
-                 then if a > c
-                      then a
-                      else c
-                 else if b > c
-                      then b
-                      else c
-
--- | Min of 4 values.
-fmin4 :: Flt -> Flt -> Flt -> Flt -> Flt
-fmin4 !a !b !c !d = fmin (fmin a b) (fmin c d)
-
--- | Max of 4 values.
-fmax4 :: Flt -> Flt -> Flt -> Flt -> Flt
-fmax4 !a !b !c !d = fmax (fmax a b) (fmax c d)
-
--- | Non-polymorphic absolute value.
-fabs :: Flt -> Flt
-fabs !a = 
- if a < 0 then (-a) else a
-
--- | Non-polymorphic integer absolute value.
-iabs :: Int -> Int
-iabs !a =
- if a < 0 then (-a) else a
-
--- | Force user to use fabs or iabs, for performance reasons.  Not sure if
--- this really helps, though.
-abs a = error "use non-polymorphic version, fabs"
-
--- | Approximate equality for Flt.  True if a and b are "almost" equal.
--- The (abs $ a-b) test doesn't work if
--- a and b are large.
-about_equal :: Flt -> Flt -> Bool
-about_equal !a !b =
- if a > 1 
- then
-  fabs (1 - (a/b)) < (delta*10) 
- else
-  (fabs $ a - b) < (delta*10)
-
--- | 3d type represented as a record of unboxed floats.
-data Vec = Vec !Flt !Flt !Flt deriving Show
-
--- | A Ray is made up of an origin and direction Vec.
-data Ray = Ray {origin, dir :: !Vec} deriving Show
---data Plane = Plane {norm :: !Vec, offset :: !Flt} deriving Show
-
--- | Vec constructor.
-vec :: Flt -> Flt -> Flt -> Vec
-vec !x !y !z = (Vec x y z)
-
--- | Zero Vec.
-vzero :: Vec
-vzero = Vec 0.0 0.0 0.0
-
--- | For when we need a unit vector, but we 
--- don't care where it points.
-vunit :: Vec
-vunit = vx
-
--- | Unit X vector.
-vx :: Vec
-vx  = Vec 1 0 0
-
--- | Unit y vector.
-vy :: Vec
-vy  = Vec 0 1 0
-
--- | Unit z vector.
-vz :: Vec
-vz  = Vec 0 0 1
-
--- | Negative x vector.
-nvx :: Vec
-nvx = Vec (-1) 0 0
-
--- | Negative y vector.
-nvy :: Vec
-nvy = Vec 0 (-1) 0
-
--- | Negative z vector.
-nvz :: Vec
-nvz = Vec 0 0 (-1)
-
--- Extract x coordinate.
-x :: Vec -> Flt
-x (Vec x_ _ _) = x_
-
--- Extract y coordinate.
-y :: Vec -> Flt
-y (Vec _ y_ _) = y_
-
--- Extract z coordinate.
-z :: Vec -> Flt
-z (Vec _ _ z_) = z_
-
--- | Access the Vec as if it were an array indexed from 0..2.
--- Note: this actually accounts for a noticeable amount of cpu 
--- time in the Glome ray tracer.
-va :: Vec -> Int -> Flt
-va !(Vec x y z) !n = 
- case n of
-  0 -> x
-  1 -> y
-  2 -> z
-
--- | Create a new Vec with the Nth field overwritten by new value.
--- I could have used record update syntax.
-vset :: Vec -> Int -> Flt -> Vec
-vset !(Vec x y z) !i !f =
- case i of
-  0 -> Vec f y z
-  1 -> Vec x f z
-  2 -> Vec x y f
-
--- | Dot product of 2 vectors.  We use this all the time.  Dot product of 2
--- normal vectors is the cosine of the angle between them.
-vdot :: Vec -> Vec -> Flt
-vdot !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
- (x1*x2)+(y1*y2)+(z1*z2)
-
--- | Cross product of 2 vectors.  Produces a vector perpendicular 
--- to the given vectors.  We use this for things like making the forward,
--- up, and right camera vectors orthogonal.  If the input vectors are
--- normalized, the output vector will be as well.
-vcross :: Vec -> Vec -> Vec
-vcross !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
- Vec 
-  ((y1 * z2) - (z1 * y2))
-  ((z1 * x2) - (x1 * z2))
-  ((x1 * y2) - (y1 * x2))
-
--- | Apply a unary Flt operator to each field of the Vec.
-vmap :: (Flt -> Flt) -> Vec -> Vec
-vmap f !v1 = 
- Vec (f (x v1)) (f (y v1)) (f (z v1))
-
--- | Apply a binary Flt operator to pairs of fields from 2 Vecs.
-vmap2 :: (Flt -> Flt -> Flt) -> Vec -> Vec -> Vec
-vmap2 f !v1 !v2 =
- Vec (f (x v1) (x v2)) 
-     (f (y v1) (y v2)) 
-     (f (z v1) (z v2))
-
--- | Reverse the direction of a Vec.
-vinvert :: Vec -> Vec
-vinvert !(Vec x1 y1 z1) =
- Vec (-x1) (-y1) (-z1)
-
--- | Get the length of a Vec squared.  We use this to avoid a slow sqrt. 
-vlensqr :: Vec -> Flt
-vlensqr !v1 = vdot v1 v1
-
--- | Get the length of a Vec.  This is expensive because sqrt is slow.
-vlen :: Vec -> Flt
-vlen !v1 = sqrt (vdot v1 v1)
-
--- | Add 2 vectors.
-vadd :: Vec -> Vec -> Vec
-vadd !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
- Vec (x1 + x2)
-     (y1 + y2)
-     (z1 + z2)
-
--- | Add 3 vectors.
-vadd3 :: Vec -> Vec -> Vec -> Vec
-vadd3 !(Vec x1 y1 z1) !(Vec x2 y2 z2) !(Vec x3 y3 z3) =
-    Vec (x1 + x2 + x3)
-        (y1 + y2 + y3)
-        (z1 + z2 + z3)
-
--- | Subtract vectors.  "vsub b a" is the vector from a to b.
-vsub :: Vec -> Vec -> Vec
-vsub !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
- Vec (x1 - x2)
-     (y1 - y2)
-     (z1 - z2)
-
--- | Multiply corresponding fields.  Rarely useful.
-vmul :: Vec -> Vec -> Vec
-vmul !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
- Vec (x1 * x2)
-     (y1 * y2)
-     (z1 * z2)
-
--- | Add a value to all the fields of a Vec.  Useful, for instance, to get
--- one corner of the bounding box around a sphere.
-vinc :: Vec -> Flt -> Vec
-vinc !(Vec x y z) !n =
- Vec (x + n)
-     (y + n)
-     (z + n)
-
--- | Subtract a value from all fields of a Vec.
-vdec :: Vec -> Flt -> Vec
-vdec !(Vec x y z) !n =
- Vec (x - n)
-     (y - n)
-     (z - n)
-
--- | Get the maximum of all corresponding fields between 2 Vecs.
-vmax :: Vec -> Vec -> Vec
-vmax !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
- Vec (fmax x1 x2)
-     (fmax y1 y2)
-     (fmax z1 z2)
-
--- | Get the minimum of all corresponding fields between 2 Vecs.
-vmin :: Vec -> Vec -> Vec
-vmin !(Vec x1 y1 z1) !(Vec x2 y2 z2) =
- Vec (fmin x1 x2)
-     (fmin y1 y2)
-     (fmin z1 z2)
-
--- | Return the largest axis.  Often used with "va".
-vmaxaxis :: Vec -> Int
-vmaxaxis !(Vec x y z) =
- if (x > y) 
- then if (x > z) 
-      then 0
-      else 2
- else if (y > z) 
-      then 1
-      else 2
-
--- | Scale a Vec by some value.
-vscale :: Vec -> Flt -> Vec
-vscale !(Vec x y z) !fac =
- Vec (x * fac)
-     (y * fac)
-     (z * fac)
-
--- | Take the first Vec, and add to it the second Vec scaled by some amount.
--- This is used quite a lot in Glome.
-vscaleadd :: Vec -> Vec -> Flt -> Vec
-vscaleadd !(Vec x1 y1 z1) !(Vec x2 y2 z2) fac =
- Vec (x1 + (x2 * fac))
-     (y1 + (y2 * fac))
-     (z1 + (z2 * fac))
-            
--- | Make the length of a Vec just a little shorter.
-vnudge :: Vec -> Vec
-vnudge x = vscale x (1-delta)
-
--- | Normalize a vector.  Division is expensive, so we compute the reciprocol 
--- of the length and multiply by that.  The sqrt is also expensive.
-vnorm :: Vec -> Vec
-vnorm !(Vec x1 y1 z1) = 
- let !invlen = 1.0 / (sqrt ((x1*x1)+(y1*y1)+(z1*z1))) in
- Vec (x1*invlen) (y1*invlen) (z1*invlen)
-
--- | Throw an exception if a vector hasn't been normalized.
-assert_norm :: Vec -> Vec
-assert_norm v =
- let l = vdot v v
- in if l > (1+delta) 
-    then error $ "vector too long" ++ (show v)
-    else if l < (1-delta)
-         then error $ "vector too short: " ++ (show v)
-         else v
-
--- | Get the victor bisecting two other vectors (which ought to be the same
--- length).
-bisect :: Vec -> Vec -> Vec
-bisect !v1 !v2 = vnorm (vadd v1 v2)
-
--- | Distance between 2 vectors.
-vdist :: Vec -> Vec -> Flt
-vdist v1 v2 = 
- let d = vsub v2 v1 in vlen d
-
--- | Reflect a vector "v" off of a surface with normal "norm".
-reflect :: Vec -> Vec -> Vec
-reflect !v !norm =
-  -- vadd v $ vscale norm $ (-2) * (vdot v norm)
-  vscaleadd v norm $ (-2) * (vdot v norm)
-
--- | Reciprocol of all fields of a Vec.
-vrcp :: Vec -> Vec
-vrcp !(Vec x y z) =
- Vec (1/x) (1/y) (1/z)
-
--- | Test Vecs for approximate equality
-veq :: Vec -> Vec -> Bool
-veq !(Vec ax ay az) !(Vec bx by bz) =
- (about_equal ax bx) && (about_equal ay by) && (about_equal az bz)
-
--- | Test Vecs for matching sign on all fields.  Returns false if any value is
--- zero.  Used by packet tracing.
-veqsign :: Vec -> Vec -> Bool
-veqsign !(Vec ax ay az) !(Vec bx by bz) =
- ax*bx > 0 && ay*by > 0 && az*bz > 0
-
--- | Translate a ray's origin in ray's direction by d amount.
-ray_move :: Ray -> Flt -> Ray
-ray_move !(Ray orig dir) !d =
- (Ray (vscaleadd orig dir d) dir)
-
--- | Find a pair of orthogonal vectors to the one given.
-orth :: Vec -> (Vec,Vec)
-orth v1 =
- if about_equal (vdot v1 v1) 1
- then
-  let x = (Vec 1 0 0)
-      y = (Vec 0 1 0)
-      dvx = vdot v1 x
-      v2 = if dvx < 0.8 && dvx > (-0.8) -- don't want to cross with a
-           then vnorm $ vcross v1 x     -- vector that's too similar
-           else vnorm $ vcross v1 y
-      v3 = vcross v1 v2
-  in (v2,v3)
- else error $ "orth: unnormalized vector" ++ (show v1)
-
--- | Intersect a ray with a plane 
--- defined by a point "p" and a normal "norm".
--- (Ray does not need to be normalized.)
-plane_int :: Ray -> Vec -> Vec -> Vec
-plane_int !(Ray orig dir) !p !norm =
- let newo = vsub orig p
-     dist = -(vdot norm newo) / (vdot norm dir)
- in vscaleadd orig dir dist
-
--- | Find the distance along a ray until it intersects with a plane defined
--- by a point "p" and normal "norm".
-plane_int_dist :: Ray -> Vec -> Vec -> Flt
-plane_int_dist !(Ray orig dir) !p !norm =
- let newo = vsub orig p
- in -(vdot norm newo) / (vdot norm dir)
-
--- find intersection with plane
--- from graphics gems -- an efficient ray-polygon intersection
--- it seems that the ray need not be normalized
--- let plane_intersect ray (n,d) =
---  let t = -.((d +. (vdot n ray.origin)) /. (vdot n ray.dir))
---  in vadd ray.origin (vscale ray.dir t)
-
-
--- TRANSFORMATIONS --
-
--- | 3x4 Transformation matrix.  These are described in most graphics texts.
-data Matrix = Matrix !Flt !Flt !Flt !Flt  
-                     !Flt !Flt !Flt !Flt  
-                     !Flt !Flt !Flt !Flt deriving Show
-
--- | A transformation.  Inverting a matrix is expensive, so we keep a forward
--- transformation matrix and a reverse transformation matrix.
--- Note: This can be made a little faster if the matricies are non-strict.
-data Xfm = Xfm Matrix Matrix deriving Show
-
--- | Identity matrix.  Transforming a vector by this matrix does nothing.
-ident_matrix :: Matrix
-ident_matrix = (Matrix 1 0 0 0  0 1 0 0  0 0 1 0)
-
--- | Identity transformation.
-ident_xfm :: Xfm
-ident_xfm = Xfm ident_matrix ident_matrix
-
--- | Multiply two matricies.  This is unrolled for efficiency, and it's also
--- a little bit easier (in my opinion) to see what's going on.
-mat_mult :: Matrix -> Matrix -> Matrix
-mat_mult (Matrix a00 a01 a02 a03  a10 a11 a12 a13  a20 a21 a22 a23)
-         (Matrix b00 b01 b02 b03  b10 b11 b12 b13  b20 b21 b22 b23) =
- Matrix
-   (a00*b00 + a01*b10 + a02*b20)
-   (a00*b01 + a01*b11 + a02*b21)
-   (a00*b02 + a01*b12 + a02*b22)
-   (a00*b03 + a01*b13 + a02*b23 + a03)
-
-   (a10*b00 + a11*b10 + a12*b20)
-   (a10*b01 + a11*b11 + a12*b21)
-   (a10*b02 + a11*b12 + a12*b22)
-   (a10*b03 + a11*b13 + a12*b23 + a13)
-
-   (a20*b00 + a21*b10 + a22*b20)
-   (a20*b01 + a21*b11 + a22*b21)
-   (a20*b02 + a21*b12 + a22*b22)
-   (a20*b03 + a21*b13 + a22*b23 + a23)
-
--- | Multiply two tranformations.  This just multiplies the forward and 
--- reverse transformations.
-xfm_mult :: Xfm -> Xfm -> Xfm
-xfm_mult (Xfm a inva) (Xfm b invb) =
- Xfm (mat_mult a b) (mat_mult invb inva)
-
--- TRANSFORM UTILITY FUNCTIONS --
-
--- | There is a seemingly-magical property of transformation matricies, that
--- we can combine the effects of any number of transformations into a single
--- transformation just by multiplying them together in reverse order.  For 
--- instance, we could move a point, then rotate it about the origin by some 
--- angle around some vector, then move it again, and this can all be done by 
--- a single transformation.
--- This function combines transformations in this way, though it reverses the
--- list first so the transformations take effect in their expected order.
-compose :: [Xfm] -> Xfm
-compose xfms = check_xfm $ foldr xfm_mult ident_xfm (reverse xfms)
-
--- | Make sure a transformation is valid.  Multipy the forward and reverse
--- matrix and verify that the result is the identity matrix.
-check_xfm :: Xfm -> Xfm
-check_xfm (Xfm m i) = 
- let (Matrix m00 m01 m02 m03  
-             m10 m11 m12 m13  
-             m20 m21 m22 m23) = mat_mult m i
-     ae = about_equal
- in
-  if ae m00 1 && ae m01 0 && ae m02 0 && ae m03 0 &&
-     ae m10 0 && ae m11 1 && ae m12 0 && ae m13 0 &&
-     ae m20 0 && ae m21 0 && ae m22 1 && ae m23 0
-  then (Xfm m i)
-  else error $ "corrupt matrix " ++ (show (Xfm m i)) ++ "\n" ++ (show (mat_mult m i)) 
-
--- | Complex transformations: Rotate point (or vector) "pt" about ray by 
--- angle c.  The angle is in radians,
--- but using the angle conversion routines "deg", "rad" and "rot" is 
--- recommended.
-vrotate :: Vec -> Ray -> Flt -> Vec
-vrotate pt (Ray orig axis_) angle =
- let axis = assert_norm axis_
-     transform = compose [ translate (vinvert orig)
-                         , rotate axis angle
-                         , translate orig
-                         ]
-     new_pt = xfm_point transform pt
- in if about_equal (vlen (vsub orig pt)) (vlen (vsub orig new_pt))
-    then new_pt
-    else error $ "something is wrong with vrotate" ++ 
-                 (show $ vlen (vsub orig pt)) ++ " " ++ 
-                 (show $ vlen (vsub orig new_pt))
-
-
--- TRANSFORM APPLICATION --
--- these need to be fast
-
--- | Transform a point.  The point is treated as (x y z 1).
-xfm_point :: Xfm -> Vec -> Vec
-xfm_point !(Xfm (Matrix m00 m01 m02 m03  
-                        m10 m11 m12 m13  
-                        m20 m21 m22 m23) inv) 
-          !(Vec x y z) =
- Vec (m00*x + m01*y + m02*z + m03)
-     (m10*x + m11*y + m12*z + m13)
-     (m20*x + m21*y + m22*z + m23)
-
--- | Inverse transform a point.
-invxfm_point :: Xfm -> Vec -> Vec
-invxfm_point !(Xfm fwd (Matrix i00 i01 i02 i03  
-                               i10 i11 i12 i13  
-                               i20 i21 i22 i23)) 
-             !(Vec x y z) =
-  Vec (i00*x + i01*y + i02*z + i03)
-      (i10*x + i11*y + i12*z + i13)
-      (i20*x + i21*y + i22*z + i23)
-
--- | Transform a vector.  The vector is treated as (x y z 0).
-xfm_vec :: Xfm -> Vec -> Vec
-xfm_vec !(Xfm (Matrix m00 m01 m02 m03  
-                      m10 m11 m12 m13  
-                      m20 m21 m22 m23) inv) 
-        !(Vec x y z) =
- Vec (m00*x + m01*y + m02*z)
-     (m10*x + m11*y + m12*z)
-     (m20*x + m21*y + m22*z)
-
--- | Inverse transform a vector.
-invxfm_vec :: Xfm -> Vec -> Vec
-invxfm_vec !(Xfm fwd (Matrix i00 i01 i02 i03  
-                             i10 i11 i12 i13  
-                             i20 i21 i22 i23)) 
-           !(Vec x y z) =
-  Vec (i00*x + i01*y + i02*z)
-      (i10*x + i11*y + i12*z)
-      (i20*x + i21*y + i22*z)
-
--- | Inverse transform a normal.  This one is tricky: we need to transform 
--- by the inverse transpose.
-invxfm_norm :: Xfm -> Vec -> Vec
-invxfm_norm !(Xfm fwd (Matrix i00 i01 i02 i03  
-                              i10 i11 i12 i13  
-                              i20 i21 i22 i23)) 
-            !(Vec x y z) =
- Vec (i00*x + i10*y + i20*z)
-     (i01*x + i11*y + i21*z)
-     (i02*x + i12*y + i22*z)
-
--- | Transform a Ray.
-xfm_ray :: Xfm -> Ray -> Ray
-xfm_ray !xfm !(Ray orig dir) =
- Ray (xfm_point xfm orig) (vnorm (xfm_vec xfm dir))
-
--- | Inverse transform a Ray.
-invxfm_ray :: Xfm -> Ray -> Ray
-invxfm_ray !xfm !(Ray orig dir) =
- Ray (invxfm_point xfm orig) (vnorm (invxfm_vec xfm dir))
-
--- BASIC TRANSFORMS --
--- | Basic transforms: move by some displacement vector.
-translate :: Vec -> Xfm
-translate (Vec x y z) =
- check_xfm $ Xfm (Matrix 1 0 0   x   0 1 0   y   0 0 1   z) 
-                 (Matrix 1 0 0 (-x)  0 1 0 (-y)  0 0 1 (-z))
-
--- | Basic transforms: stretch along the three axes, by the amount
--- in the given vector.  (If x==y==z, then it's uniform scaling.)
-scale :: Vec -> Xfm
-scale (Vec x y z) =
- check_xfm $ Xfm (Matrix   x  0 0 0  0   y  0 0  0 0   z  0)
-                (Matrix (1/x) 0 0 0  0 (1/y) 0 0  0 0 (1/z) 0)
-
--- | Basic transforms: rotate about a given axis by some angle.
-rotate :: Vec -> Flt -> Xfm
-rotate v@(Vec x y z) angle =
- if not $ (vlen v) `about_equal` 1
- then error $ "please use a normalized vector for rotation: " ++ (show (vlen v))
- else 
-  let s = sin angle
-      c = cos angle 
-
-      m00 = ((x*x)+((1-(x*x))*c)) 
-      m01 = (((x*y)*(1-c))-(z*s)) 
-      m02 = ((x*z*(1-c))+(y*s))
-
-      m10 = (((x*y)*(1-c))+(z*s))
-      m11 = ((y*y)+((1-(y*y))*c))
-      m12 = ((y*z*(1-c))-(x*s))
-
-      m20 = ((x*z*(1-c))-(y*s))
-      m21 = ((y*z*(1-c))+(x*s))
-      m22 = ((z*z)+((1-(z*z))*c))
-  in
-  check_xfm $ Xfm (Matrix m00 m01 m02 0  m10 m11 m12 0  m20 m21 m22 0)
-                  (Matrix m00 m10 m20 0  m01 m11 m21 0  m02 m12 m22 0)
-
--- | Basic transforms: Convert coordinate system from canonical xyz 
--- coordinates to uvw coordinates.
-xyz_to_uvw :: Vec -> Vec -> Vec -> Xfm
-xyz_to_uvw u v w =
- let Vec ux uy uz = u
-     Vec vx vy vz = v
-     Vec wx wy wz = w
- in if (vdot u u) `about_equal` 1
-    then
-     if (vdot v v) `about_equal` 1
-     then
-      if (vdot w w) `about_equal` 1
-      then 
-       if ((vdot u v) `about_equal` 0) && 
-          ((vdot u w) `about_equal` 0) && 
-          ((vdot v w) `about_equal` 0)
-       then
-        check_xfm $ Xfm (Matrix ux vx wx 0  uy vy wy 0  uz vz wz 0)
-                        (Matrix ux uy uz 0  vx vy vz 0  wx wy wz 0)
-       else error  "vectors aren't orthogonal"
-      else error $ "unnormalized w " ++ (show w)
-     else error $ "unnormalized v " ++ (show v)
-    else error $ "unnormalized u " ++ (show u)
-
--- | Basic transforms: Convert from uvw coordinates back to normal xyz 
--- coordinates.
-uvw_to_xyz :: Vec -> Vec -> Vec -> Xfm
-uvw_to_xyz (Vec ux uy uz) (Vec vx vy vz) (Vec wx wy wz) =
- check_xfm $ Xfm (Matrix ux uy uz 0  vx vy vz 0  wx wy wz 0)
-                 (Matrix ux vx wx 0  uy vy wy 0  uz vz wz 0)
-
-
-
--- TRIANGLE UTILITY FUNCTIONS --
-
--- | Given a side, angle, and side of a triangle, produce the length of the 
--- opposite side.
-sas2s :: Flt -> Flt -> Flt -> Flt
-sas2s s1 a s2 =
-  sqrt (((s1 * s1) + (s2 * s2)) - ((2 * s1 * s2 * (dcos a))))
-
-
-
--- BOUNDING BOXES --
--- | Axis-aligned Bounding Box (AABB), defined by opposite corners.  P1 is the
--- min values, p2 has the max values.
-data Bbox = Bbox {p1 :: !Vec, p2 :: !Vec} deriving Show
-
--- | A near-far pair of distances.  Basically just a tuple.
-data Interval = Interval !Flt !Flt deriving Show -- used instead of a tuple
-
--- | Bounding box that encloses two bounding boxes.
-bbjoin :: Bbox -> Bbox -> Bbox
-bbjoin (Bbox p1a p2a) (Bbox p1b p2b) =
- (Bbox (vmin p1a p1b) (vmax p2a p2b))
-
--- | Find the overlap of two bounding boxes.
-bboverlap :: Bbox -> Bbox -> Bbox
-bboverlap (Bbox p1a p2a) (Bbox p1b p2b) =
- (Bbox (vmax p1a p1b) (vmin p2a p2b))
-
--- | Test if a Vec is inside the bounding box.
-bbinside :: Bbox -> Vec -> Bool
-bbinside (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) (Vec x y z) =
- p1x <= x && x <= p2x && p1y <= y && y <= p2y && p1z <= z && z <= p2z
-
--- | Split a bounding box into two, given an axis and offset.  Throw exception
--- if the offset isn't inside the bounding box.
-bbsplit :: Bbox -> Int -> Flt -> (Bbox,Bbox)
-bbsplit (Bbox p1 p2) axis offset =
- if (offset < (va p1 axis)) || (offset > (va p2 axis))
- then error "degenerate bounding box split"
- else ((Bbox p1 (vset p2 axis offset)),
-       (Bbox (vset p1 axis offset) p2))
-
--- | Generate a minimum bounding box that encloses a list of points.
-bbpts :: [Vec] -> Bbox
-bbpts [] = empty_bbox
-bbpts ((Vec x y z):[]) =
- Bbox (Vec (x-delta) (y-delta) (z-delta)) 
-      (Vec (x+delta) (y+delta) (z+delta))
-
-bbpts ((Vec x y z):pts) =
- let (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) = bbpts pts
-     minx = fmin (x-delta) p1x
-     miny = fmin (y-delta) p1y
-     minz = fmin (z-delta) p1z
-     maxx = fmax (x+delta) p2x
-     maxy = fmax (y+delta) p2y
-     maxz = fmax (z+delta) p2z in
- Bbox (Vec minx miny minz) (Vec maxx maxy maxz)
-
--- | Surface area of a bounding box.  Useful for cost heuristics when attempting
--- to build optimal bounding box heirarchies.  Undefined for degenerate bounding
--- boxes.
-bbsa :: Bbox -> Flt
-bbsa (Bbox p1 p2) =
- let Vec dx dy dz = vsub p2 p1 
- in dx*dy + dx*dz + dy*dz
-
--- | Volume of a bounding box.  Undefined for degenerate bounding boxes.
-bbvol :: Bbox -> Flt
-bbvol (Bbox p1 p2) =
- let (Vec dx dy dz) = vsub p2 p1
- in dx*dy*dz
-
--- | Degenerate bounding box that contains an empty volume.
-empty_bbox :: Bbox
-empty_bbox = 
- Bbox (Vec infinity infinity infinity) 
-      (Vec (-infinity) (-infinity) (-infinity))
-
--- | "Infinite" bounding box.
-everything_bbox :: Bbox
-everything_bbox =
- Bbox (Vec (-infinity) (-infinity) (-infinity))
-      (Vec infinity infinity infinity)
-
--- | Find a ray's entrance and exit from a bounding 
--- box.  If last entrance is before the first exit,
--- we hit.  Otherwise, we miss. (It's up to the 
--- caller to figure that out.)
-
-bbclip :: Ray -> Bbox -> Interval
-bbclip (Ray (Vec ox oy oz) (Vec dx dy dz)) 
-       (Bbox (Vec p1x p1y p1z) (Vec p2x p2y p2z)) =
- let dxrcp = 1/dx
-     dyrcp = 1/dy
-     dzrcp = 1/dz
-     Interval inx outx = if dx > 0 
-                         then Interval ((p1x-ox)*dxrcp) ((p2x-ox)*dxrcp)
-                         else Interval ((p2x-ox)*dxrcp) ((p1x-ox)*dxrcp)
-     Interval iny outy = if dy > 0
-                         then Interval ((p1y-oy)*dyrcp) ((p2y-oy)*dyrcp)
-                         else Interval ((p2y-oy)*dyrcp) ((p1y-oy)*dyrcp)
-     Interval inz outz = if dz > 0
-                         then Interval ((p1z-oz)*dzrcp) ((p2z-oz)*dzrcp)
-                         else Interval ((p2z-oz)*dzrcp) ((p1z-oz)*dzrcp)
- in
-   Interval (fmax3 inx iny inz) (fmin3 outx outy outz)
-
-
diff --git a/GlomeVec.cabal b/GlomeVec.cabal
--- a/GlomeVec.cabal
+++ b/GlomeVec.cabal
@@ -1,5 +1,5 @@
 Name:                GlomeVec
-Version:             0.1.1
+Version:             0.1.2
 Synopsis:            Simple 3D vector library
 Description:         A simple library for dealing with 3D vectors, suitable for graphics projects.  A small texture library with Perlin noise is included as well.
 License:             GPL
@@ -14,7 +14,7 @@
 Cabal-Version: >= 1.2
 extra-source-files:
   README.txt
-
 library
-  exposed-modules:   Data.GlomeVec, Data.GlomeTexture
+  exposed-modules:   Data.Glome.Vec
+                     Data.Glome.Texture
   Build-Depends:     base >= 3 && < 4, array