Math-Complex: improve numerical accuracy

dist/Math-Complex/lib/Math/Complex.pm

@@ -10,7 +10,7 @@ package Math::Complex;
 { use 5.006; }
 use strict;
 
-our $VERSION = 1.59_03;
+our $VERSION = 1.59_04;
 
 use Config;
 

dist/Math-Complex/lib/Math/Trig.pm

@@ -15,7 +15,7 @@ require Exporter;
 
 our @ISA = qw(Exporter);
 
-our $VERSION = 1.23_01;
+our $VERSION = 1.59_04;
 
 my @angcnv = qw(rad2deg rad2grad
 		deg2rad deg2grad
@@ -154,12 +154,19 @@ sub great_circle_distance {
 
     $rho = 1 unless defined $rho; # Default to the unit sphere.
 
-    my $lat0 = pip2 - $phi0;
-    my $lat1 = pip2 - $phi1;
+    my $dphi   = $phi1 - $phi0;
+    my $dtheta = $theta1 - $theta0;
 
-    return $rho *
-	acos_real( cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
-		   sin( $lat0 ) * sin( $lat1 ) );
+    # A formula that is accurate for all distances is the following special
+    # case of the Vincenty formula for an ellipsoid with equal major and minor
+    # axes.  See
+    # https://en.wikipedia.org/wiki/Great-circle_distance#Computational_formulas
+
+    my $c1 = sin($phi1) * sin($dtheta);
+    my $c2 = sin($phi1) * cos($dtheta);
+    my $c3 = sin($phi0) * cos($phi1) - cos($phi0) * $c2;
+    my $c4 = cos($phi0) * cos($phi1) + sin($phi0) * $c2;
+    return $rho * atan2(sqrt($c1 * $c1 + $c3 * $c3), $c4);
 }
 
 sub great_circle_direction {
@@ -538,26 +545,19 @@ points.
 
 =head2 great_circle_distance
 
-You can compute spherical distances, called B<great circle distances>,
-by importing the great_circle_distance() function:
-
-  use Math::Trig 'great_circle_distance';
+Returns the great circle distance between two points on a sphere.
 
-  $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
+    $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
 
-The I<great circle distance> is the shortest distance between two
-points on a sphere.  The distance is in C<$rho> units.  The C<$rho> is
-optional, it defaults to 1 (the unit sphere), therefore the distance
-defaults to radians.
+Where ($theta0, $phi0) and ($theta1, $phi1) are the spherical coordinates of
+the two points, respectively. The distance is in C<$rho> units. The C<$rho>
+is optional. It defaults to 1 (the unit sphere).
 
-If you think geographically the I<theta> are longitudes: zero at the
-Greenwhich meridian, eastward positive, westward negative -- and the
-I<phi> are latitudes: zero at the North Pole, northward positive,
-southward negative.  B<NOTE>: this formula thinks in mathematics, not
-geographically: the I<phi> zero is at the North Pole, not at the
-Equator on the west coast of Africa (Bay of Guinea).  You need to
-subtract your geographical coordinates from I<pi/2> (also known as 90
-degrees).
+If you are using geographic coordinates, latitude and longitude, you need to
+adjust for the fact that latitude is zero at the equator increasing towards
+the north and decreasing towards the south. Assuming ($lat0, $lon0) and
+($lat1, $lon1) are the geographic coordinates in radians of the two points,
+the distance can be computed with
 
   $distance = great_circle_distance($lon0, pi/2 - $lat0,
                                     $lon1, pi/2 - $lat1, $rho);

dist/Math-Complex/t/Trig.t

@@ -10,7 +10,7 @@
 
 use strict;
 use warnings;
-use Test::More tests => 153;
+use Test::More tests => 157;
 
 use Math::Trig 1.18;
 use Math::Trig 1.18 qw(:pi Inf);
@@ -363,7 +363,15 @@ print "# great_circle_distance with small angles\n";
 for my $e (qw(1e-2 1e-3 1e-4 1e-5)) {
     # Can't assume == 0 because of floating point fuzz,
     # but let's hope for at least < $e.
-    cmp_ok(great_circle_distance(0, $e, 0, $e), '<', $e);
+    cmp_ok(great_circle_distance(0, $e, 0, $e), '<', $e,
+           "great_circle_distance(0, $e, 0, $e) < $e");
+}
+
+for my $e (qw(1e-5 1e-6 1e-7 1e-8)) {
+    # Verify that the distance is positive for points close together. A poor
+    # algorithm is likely to give a distance of zero in some of these cases.
+    cmp_ok(great_circle_distance(2, 2, 2, 2+$e), '>', 0,
+           "great_circle_distance(2, 2, 2, " . (2+$e) . ") > 0");
 }
 
 print "# asin_real, acos_real\n";