[MOAB-dev] r6022 - MOAB/trunk/tools/mbcslam

[iulian at mcs.anl.gov]

Author: iulian
Date: 2013-03-03 22:14:27 -0600 (Sun, 03 Mar 2013)
New Revision: 6022

Modified:
   MOAB/trunk/tools/mbcslam/CslamUtils.cpp
   MOAB/trunk/tools/mbcslam/CslamUtils.hpp
   MOAB/trunk/tools/mbcslam/intx_on_sphere_test.cpp
   MOAB/trunk/tools/mbcslam/spherical_area_test.cpp
Log:
use l'Huiller formula to compute the area of triangles on a sphere
http://williams.best.vwh.net/avform.htm
compared to Girard's formula, this is numerically better for small triangles


Modified: MOAB/trunk/tools/mbcslam/CslamUtils.cpp
===================================================================
--- MOAB/trunk/tools/mbcslam/CslamUtils.cpp	2013-03-01 19:45:57 UTC (rev 6021)
+++ MOAB/trunk/tools/mbcslam/CslamUtils.cpp	2013-03-04 04:14:27 UTC (rev 6022)
@@ -671,7 +671,66 @@
   return Radius*Radius *correction;
 
 }
+/*
+ *  l'Huiller's formula for spherical triangle
+ *  http://williams.best.vwh.net/avform.htm
+ *  a, b, c are arc measures in radians, too
+ *  A, B, C are angles on the sphere, for which we already have formula
+ *               c
+ *         A -------B
+ *          \       |
+ *           \      |
+ *            \b    |a
+ *             \    |
+ *              \   |
+ *               \  |
+ *                \C|
+ *                 \|
+ *
+ *  (The angle at B is not necessarily a right angle)
+ *
+ *    sin(a)  sin(b)   sin(c)
+ *    ----- = ------ = ------
+ *    sin(A)  sin(B)   sin(C)
+ *
+ * In terms of the sides (this is excess, as before, but numerically stable)
+ *
+ *  E = 4*atan(sqrt(tan(s/2)*tan((s-a)/2)*tan((s-b)/2)*tan((s-c)/2)))
+ */
+double area_spherical_triangle_lHuiller(double * ptA, double * ptB, double * ptC, double Radius)
+{
+  // now, a is angle BOC, O is origin
+  CartVect vA(ptA), vB(ptB), vC(ptC);
+  double a = angle(vB, vC);
+  double b = angle(vC, vA);
+  double c = angle(vA, vB);
+  int sign =1;
+  if ((vA*vB)%vC < 0)
+    sign = -1;
+  double s = (a+b+c)/2;
+  double E = 4*atan(sqrt(tan(s/2)*tan((s-a)/2)*tan((s-b)/2)*tan((s-c)/2)));
+  return sign * E * Radius*Radius;
+}
 
+
+double area_spherical_polygon_lHuiller (double * A, int N, double Radius)