$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 {
=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);
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);
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";