--- | The matrix used in the tetrahedron volume calculation as given in
--- Lai & Schumaker, Definition 15.4, page 436.
-vol_matrix :: Tetrahedron -> Matrix Double
-vol_matrix t = (4><4)
- [1, 1, 1, 1,
- x1, x2, x3, x4,
- y1, y2, y3, y4,
- z1, z2, z3, z4 ]
- where
- (x1, y1, z1) = v0 t
- (x2, y2, z2) = v1 t
- (x3, y3, z3) = v2 t
- (x4, y4, z4) = v3 t
+-- | Compute the determinant of the 4x4 matrix,
+--
+-- [1]
+-- [x]
+-- [y]
+-- [z]
+--
+-- where [1] = [1, 1, 1, 1],
+-- [x] = [x1,x2,x3,x4],
+--
+-- et cetera.
+--
+-- The termX nonsense is an attempt to prevent Double overflow.
+-- which has been observed to happen with large coordinates.
+--
+det :: Point -> Point -> Point -> Point -> Double
+det p0 p1 p2 p3 =
+ term5 + term6
+ where
+ (x1, y1, z1) = p0
+ (x2, y2, z2) = p1
+ (x3, y3, z3) = p2
+ (x4, y4, z4) = p3
+ term1 = ((x2 - x4)*y1 - (x1 - x4)*y2 + (x1 - x2)*y4)*z3
+ term2 = ((x2 - x3)*y1 - (x1 - x3)*y2 + (x1 - x2)*y3)*z4
+ term3 = ((x3 - x4)*y2 - (x2 - x4)*y3 + (x2 - x3)*y4)*z1
+ term4 = ((x3 - x4)*y1 - (x1 - x4)*y3 + (x1 - x3)*y4)*z2
+ term5 = term1 - term2
+ term6 = term3 - term4
+