From f27803a7c59e049a2fc70327db9e01f55c71df4e Mon Sep 17 00:00:00 2001 From: Peter John Acklam Date: Thu, 1 Sep 2022 20:42:35 +0200 Subject: [PATCH 1/7] Math-Complex: improve numerical accuracy - Compute the great circle distance with a formula that has better numerical properties. The formula used now is accurate for all distances. - Add a few tests to verify. - This fixes CPAN RT #78938 --- dist/Math-Complex/lib/Math/Trig.pm | 19 +++++++++++++------ dist/Math-Complex/t/Trig.t | 12 ++++++++++-- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/dist/Math-Complex/lib/Math/Trig.pm b/dist/Math-Complex/lib/Math/Trig.pm index 50daec57e5c6..5c0fcc239a6b 100644 --- a/dist/Math-Complex/lib/Math/Trig.pm +++ b/dist/Math-Complex/lib/Math/Trig.pm @@ -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; - - return $rho * - acos_real( cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) + - sin( $lat0 ) * sin( $lat1 ) ); + my $dphi = $phi1 - $phi0; + my $dtheta = $theta1 - $theta0; + + # 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 { diff --git a/dist/Math-Complex/t/Trig.t b/dist/Math-Complex/t/Trig.t index a9a12556b659..8e2e5d164577 100644 --- a/dist/Math-Complex/t/Trig.t +++ b/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"; From ffa942153e82a1a9067b2123cc8275d7bea643ed Mon Sep 17 00:00:00 2001 From: Peter John Acklam Date: Sat, 5 Nov 2022 11:21:14 +0100 Subject: [PATCH 2/7] Improve documentation for great_circle_distance() --- dist/Math-Complex/lib/Math/Trig.pm | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/dist/Math-Complex/lib/Math/Trig.pm b/dist/Math-Complex/lib/Math/Trig.pm index 5c0fcc239a6b..8b8279173462 100644 --- a/dist/Math-Complex/lib/Math/Trig.pm +++ b/dist/Math-Complex/lib/Math/Trig.pm @@ -545,26 +545,20 @@ points. =head2 great_circle_distance -You can compute spherical distances, called B, -by importing the great_circle_distance() function: +Returns the great circle distance between two points on a sphere. - use Math::Trig 'great_circle_distance'; + $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]); - $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]); - -The I 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 +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), therefore the distance defaults to radians. -If you think geographically the I are longitudes: zero at the -Greenwhich meridian, eastward positive, westward negative -- and the -I are latitudes: zero at the North Pole, northward positive, -southward negative. B: this formula thinks in mathematics, not -geographically: the I 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 (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); From e05216441959c091f660fd1b7608a4a2915d245f Mon Sep 17 00:00:00 2001 From: Peter John Acklam Date: Sat, 5 Nov 2022 19:00:09 +0100 Subject: [PATCH 3/7] Improve documentation for great_circle_distance() --- dist/Math-Complex/lib/Math/Trig.pm | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/dist/Math-Complex/lib/Math/Trig.pm b/dist/Math-Complex/lib/Math/Trig.pm index 8b8279173462..844571976405 100644 --- a/dist/Math-Complex/lib/Math/Trig.pm +++ b/dist/Math-Complex/lib/Math/Trig.pm @@ -551,8 +551,7 @@ Returns the great circle distance between two points on a sphere. 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), therefore the distance -defaults to radians. +is optional. It defaults to 1 (the unit sphere). 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 From 32d1d67d1899dffe48916e440143d0b4f61a3d4e Mon Sep 17 00:00:00 2001 From: Peter John Acklam Date: Thu, 1 Sep 2022 20:42:35 +0200 Subject: [PATCH 4/7] Math-Complex: improve numerical accuracy - Compute the great circle distance with a formula that has better numerical properties. The formula used now is accurate for all distances. - Add a few tests to verify. - This fixes CPAN RT #78938 --- dist/Math-Complex/lib/Math/Trig.pm | 19 +++++++++++++------ dist/Math-Complex/t/Trig.t | 12 ++++++++++-- 2 files changed, 23 insertions(+), 8 deletions(-) diff --git a/dist/Math-Complex/lib/Math/Trig.pm b/dist/Math-Complex/lib/Math/Trig.pm index 50daec57e5c6..5c0fcc239a6b 100644 --- a/dist/Math-Complex/lib/Math/Trig.pm +++ b/dist/Math-Complex/lib/Math/Trig.pm @@ -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; - - return $rho * - acos_real( cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) + - sin( $lat0 ) * sin( $lat1 ) ); + my $dphi = $phi1 - $phi0; + my $dtheta = $theta1 - $theta0; + + # 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 { diff --git a/dist/Math-Complex/t/Trig.t b/dist/Math-Complex/t/Trig.t index a9a12556b659..8e2e5d164577 100644 --- a/dist/Math-Complex/t/Trig.t +++ b/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"; From e05f34517b96a8bb2dc50232b8a52d9546df67c2 Mon Sep 17 00:00:00 2001 From: Peter John Acklam Date: Sat, 5 Nov 2022 11:21:14 +0100 Subject: [PATCH 5/7] Improve documentation for great_circle_distance() --- dist/Math-Complex/lib/Math/Trig.pm | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/dist/Math-Complex/lib/Math/Trig.pm b/dist/Math-Complex/lib/Math/Trig.pm index 5c0fcc239a6b..8b8279173462 100644 --- a/dist/Math-Complex/lib/Math/Trig.pm +++ b/dist/Math-Complex/lib/Math/Trig.pm @@ -545,26 +545,20 @@ points. =head2 great_circle_distance -You can compute spherical distances, called B, -by importing the great_circle_distance() function: +Returns the great circle distance between two points on a sphere. - use Math::Trig 'great_circle_distance'; + $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]); - $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]); - -The I 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 +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), therefore the distance defaults to radians. -If you think geographically the I are longitudes: zero at the -Greenwhich meridian, eastward positive, westward negative -- and the -I are latitudes: zero at the North Pole, northward positive, -southward negative. B: this formula thinks in mathematics, not -geographically: the I 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 (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); From 1d08d11309bc16d2c75ff8d872d7f19cb23998e9 Mon Sep 17 00:00:00 2001 From: Peter John Acklam Date: Sat, 5 Nov 2022 19:00:09 +0100 Subject: [PATCH 6/7] Improve documentation for great_circle_distance() --- dist/Math-Complex/lib/Math/Trig.pm | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/dist/Math-Complex/lib/Math/Trig.pm b/dist/Math-Complex/lib/Math/Trig.pm index 8b8279173462..844571976405 100644 --- a/dist/Math-Complex/lib/Math/Trig.pm +++ b/dist/Math-Complex/lib/Math/Trig.pm @@ -551,8 +551,7 @@ Returns the great circle distance between two points on a sphere. 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), therefore the distance -defaults to radians. +is optional. It defaults to 1 (the unit sphere). 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 From dd1bc8be4160fbb3ead8523867a9ad29f6d6f7b3 Mon Sep 17 00:00:00 2001 From: Peter John Acklam Date: Sun, 6 Nov 2022 18:14:44 +0100 Subject: [PATCH 7/7] Increase version numbers --- dist/Math-Complex/lib/Math/Complex.pm | 2 +- dist/Math-Complex/lib/Math/Trig.pm | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/dist/Math-Complex/lib/Math/Complex.pm b/dist/Math-Complex/lib/Math/Complex.pm index 498b0dcf33d5..f20ce31c6867 100644 --- a/dist/Math-Complex/lib/Math/Complex.pm +++ b/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; diff --git a/dist/Math-Complex/lib/Math/Trig.pm b/dist/Math-Complex/lib/Math/Trig.pm index 844571976405..32e98531459d 100644 --- a/dist/Math-Complex/lib/Math/Trig.pm +++ b/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