00240 {
00241   double vol, l;
00242   double err;
00243 
00244 
00245   Vec3d v1 (p1, p2);
00246   Vec3d v2 (p1, p3);
00247   Vec3d v3 (p1, p4);
00248 
00249   vol = -Determinant (v1, v2, v3) / 6;
00250 
00251   double l1 = v1.Length();
00252   double l2 = v2.Length();
00253   double l3 = v3.Length();
00254   double l4 = Dist (p2, p3);
00255   double l5 = Dist (p2, p4);
00256   double l6 = Dist (p3, p4);
00257 
00258   l = l1 + l2 + l3 + l4 + l5 + l6;
00259 
00260   if (vol <= 1e-24 * l * l * l)
00261     { 
00262       return 1e24;
00263     }
00264 
00265   err = (l*l*l) / (1832.82 * vol);    // 6^4 * sqrt(2)
00266   
00267   if (h > 0)
00268     err += l / h + 
00269       h * (1 / l1 + 1/l2 + 1/l3 + 1/l4 + 1/l5 + 1/l6) - 12;
00270 
00271   return pow (err, teterrpow);
00272 }