From 740aa86cb551d6388f5cf4a8f39568d52fac6ed7 Mon Sep 17 00:00:00 2001
From: Robert Snedegar <viking@gmail.com>
Date: Thu, 11 Feb 2021 23:42:55 +0000
Subject: [PATCH] s1: Update ChordAngle comments

---
 s1/chordangle.go | 88 +++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 79 insertions(+), 9 deletions(-)

diff --git a/s1/chordangle.go b/s1/chordangle.go
index 406c69e..77d7164 100644
--- a/s1/chordangle.go
+++ b/s1/chordangle.go
@@ -26,14 +26,84 @@ import (
 // to be calculated and compared. Otherwise it is simpler to use Angle.
 //
 // ChordAngle loses some accuracy as the angle approaches π radians.
-// Specifically, the representation of (π - x) radians has an error of about
-// (1e-15 / x), with a maximum error of about 2e-8 radians (about 13cm on the
-// Earth's surface). For comparison, for angles up to π/2 radians (10000km)
-// the worst-case representation error is about 2e-16 radians (1 nanonmeter),
-// which is about the same as Angle.
-//
-// ChordAngles are represented by the squared chord length, which can
-// range from 0 to 4. Positive infinity represents an infinite squared length.
+// There are several different ways to measure this error, including the
+// representational error (i.e., how accurately ChordAngle can represent
+// angles near π radians), the conversion error (i.e., how much precision is
+// lost when an Angle is converted to an ChordAngle), and the measurement
+// error (i.e., how accurate the ChordAngle(a, b) constructor is when the
+// points A and B are separated by angles close to π radians). All of these
+// errors differ by a small constant factor.
+//
+// For the measurement error (which is the largest of these errors and also
+// the most important in practice), let the angle between A and B be (π - x)
+// radians, i.e. A and B are within "x" radians of being antipodal. The
+// corresponding chord length is
+//
+//    r = 2 * sin((π - x) / 2) = 2 * cos(x / 2)
+//
+// For values of x not close to π the relative error in the squared chord
+// length is at most 4.5 * dblEpsilon (see MaxPointError below).
+// The relative error in "r" is thus at most 2.25 * dblEpsilon ~= 5e-16. To
+// convert this error into an equivalent angle, we have
+//
+//    |dr / dx| = sin(x / 2)
+//
+// and therefore
+//
+//    |dx| = dr / sin(x / 2)
+//         = 5e-16 * (2 * cos(x / 2)) / sin(x / 2)
+//         = 1e-15 / tan(x / 2)
+//
+// The maximum error is attained when
+//
+//    x  = |dx|
+//       = 1e-15 / tan(x / 2)
+//      ~= 1e-15 / (x / 2)
+//      ~= sqrt(2e-15)
+//
+// In summary, the measurement error for an angle (π - x) is at most
+//
+//    dx  = min(1e-15 / tan(x / 2), sqrt(2e-15))
+//      (~= min(2e-15 / x, sqrt(2e-15)) when x is small)
+//
+// On the Earth's surface (assuming a radius of 6371km), this corresponds to
+// the following worst-case measurement errors:
+//
+//     Accuracy:             Unless antipodal to within:
+//     ---------             ---------------------------
+//     6.4 nanometers        10,000 km (90 degrees)
+//     1 micrometer          81.2 kilometers
+//     1 millimeter          81.2 meters
+//     1 centimeter          8.12 meters
+//     28.5 centimeters      28.5 centimeters
+//
+// The representational and conversion errors referred to earlier are somewhat
+// smaller than this. For example, maximum distance between adjacent
+// representable ChordAngle values is only 13.5 cm rather than 28.5 cm. To
+// see this, observe that the closest representable value to r^2 = 4 is
+// r^2 =  4 * (1 - dblEpsilon / 2). Thus r = 2 * (1 - dblEpsilon / 4) and
+// the angle between these two representable values is
+//
+//    x  = 2 * acos(r / 2)
+//       = 2 * acos(1 - dblEpsilon / 4)
+//      ~= 2 * asin(sqrt(dblEpsilon / 2)
+//      ~= sqrt(2 * dblEpsilon)
+//      ~= 2.1e-8
+//
+// which is 13.5 cm on the Earth's surface.
+//
+// The worst case rounding error occurs when the value halfway between these
+// two representable values is rounded up to 4. This halfway value is
+// r^2 = (4 * (1 - dblEpsilon / 4)), thus r = 2 * (1 - dblEpsilon / 8) and
+// the worst case rounding error is
+//
+//    x  = 2 * acos(r / 2)
+//       = 2 * acos(1 - dblEpsilon / 8)
+//      ~= 2 * asin(sqrt(dblEpsilon / 4)
+//      ~= sqrt(dblEpsilon)
+//      ~= 1.5e-8
+//
+// which is 9.5 cm on the Earth's surface.
 type ChordAngle float64
 
 const (
@@ -225,7 +295,7 @@ func (c ChordAngle) Sin() float64 {
 // It is more efficient than Sin.
 func (c ChordAngle) Sin2() float64 {
 	// Let a be the (non-squared) chord length, and let A be the corresponding
-	// half-angle (a = 2*sin(A)).  The formula below can be derived from:
+	// half-angle (a = 2*sin(A)). The formula below can be derived from:
 	//   sin(2*A) = 2 * sin(A) * cos(A)
 	//   cos^2(A) = 1 - sin^2(A)
 	// This is much faster than converting to an angle and computing its sine.