small refactoring of routine lat_edge_correction

* some reordering of operations
* improved documentation

src/core/area.c

@@ -214,50 +214,69 @@ static inline double XYZtoLon(double a[3]) {
   return atan2(a[1] , a[0]);
 }
 
-static double
-lat_edge_correction(double base_vec[3], double a[3], double b[3]) {
+static double lat_edge_correction(
+  double base_vec[3], double a[3], double b[3]) {
+
+  //-------------------------------------------------------------------
+  // compute the area that is enclosed between a great circle edges and
+  // a lat circle edge using that go through the points a and b
+  //-------------------------------------------------------------------
 
   double const tol = 1e-8;
 
   YAC_ASSERT(
-    fabs(a[2]-b[2]) <= tol, "ERROR: latitude of both corners is not identical")
+    fabs(a[2]-b[2]) <= tol,
+    "ERROR(lat_edge_correction_info): "
+    "latitude of both corners is not identical")
 
   double h = fabs(a[2]);
 
-  // if we are at the equator or at a pole
-  if (h < tol || fabs(1.0 - h) < tol)
-    return 0.0;
+  // if we are at the equator or at a pole -> area is zero
+  if (h < tol || fabs(1.0 - h) < tol) return 0.0;
 
+  // compute the norm vector of the plane going through a and b
+  // (if the angle between a and b is to small to compute a norm vector, then
+  //  the area is negligible)
+  double norm_ab[3];
+  if (compute_norm_vector(a, b, norm_ab)) return 0.0;
+
+  // compute the area of a triangle consisting of the lat edge and
+  // the reference pole
   double lat_area = fabs((1.0 - h) * get_angle(XYZtoLon(a), XYZtoLon(b)));
 
+  // compute the same area, but use a gc edge instead of a lat edge
   double pole[3] = {0, 0, (a[2] >= 0.0)?1.0:-1.0};
   double gc_area = tri_area(a, b, pole);
 
+  // the difference between the two triangle areas is the area enclosed
+  // between the points a and b using a lat and a gc edge
   double correction = MAX(lat_area - gc_area, 0.0);
 
+  //-------------------------------------------------------------------
+  // now we have to decide whether this correction needs to be added or
+  // subtracted
+  //-------------------------------------------------------------------
+
+  // compute the middle point of the lat edge
   double middle_lat[3] = {a[0]+b[0], a[1]+b[1], a[2]};
   double scale = sqrt(middle_lat[0]*middle_lat[0]+middle_lat[1]*middle_lat[1]);
-
   YAC_ASSERT(fabs(scale) >= 1e-18, "internal error")
-
   scale = sqrt(1.0 - a[2]*a[2]) / scale;
-
   middle_lat[0] *= scale;
   middle_lat[1] *= scale;
 
-  double norm_ab[3];
-
-  // if the angle between a and b is to small to compute a norm vector
-  if (compute_norm_vector(a, b, norm_ab)) return 0.0;
-
+  // compute the cosine between norm vector and the base vector
+  // --> their sign indicates the ordering of vertices
   double scalar_base = scalar_product(norm_ab, base_vec);
-  double scalar_middle_lat = scalar_product(norm_ab, middle_lat);
 
   // if the base point is on the same plane as a and b
   if (fabs(scalar_base) < 1e-11) {
 
+    // compute the norm vector of the lat edge middle point and the associated
+    // pole
+    // (if the angle between lat edge middle point and the pols is to small to
+    //  compute a norm vector, then the area is negligible)
     double norm_middle[3];
-
     if (compute_norm_vector(middle_lat, pole, norm_middle)) return 0.0;
 
     double scalar_a = scalar_product(norm_middle, a);
@@ -269,6 +288,9 @@ lat_edge_correction(double base_vec[3], double a[3], double b[3]) {
 
   } else {
 
+    // compute the cosine between the edge norm vector and the lat edge middle
+    // point
+    double scalar_middle_lat = scalar_product(norm_ab, middle_lat);
     if (scalar_middle_lat < 0)
       return correction;
     else