| Subject: | great_circle_distance is inaccurate when points are close together |
| Date: | Mon, 13 Aug 2012 18:59:31 +0200 |
| To: | <bug-Math-Complex [...] rt.cpan.org> |
| From: | Daniel Burr <dburr [...] topcon.com> |
The current implementation of <Math::Trig::great_circle_distance> gives
incorrect results when the points are very close together (similar to
the the problem with <Math::Trig::great_circle_direction> that I
reported in #62559). For example:
print great_circle_distance(deg2rad(149.1314), deg2rad(125.2828),
deg2rad(149.13141), deg2rad(125.2828), 6378137), "\n";
print great_circle_distance(deg2rad(149.1314), deg2rad(125.2828),
deg2rad(149.131401), deg2rad(125.2828), 6378137), "\n";
print great_circle_distance(deg2rad(149.1314), deg2rad(125.2828),
deg2rad(149.1314001), deg2rad(125.2828), 6378137), "\n";
Gives the following output:
0.90164423653156
0
0
The first value (90.2cm) is nearly correct, but the second two values
(0) are wrong. The cause is two-fold:
1. $theta0 - $theta1 is so close to zero that cos($theta0 - $theta1) gives 1
2. Even if (1) is accounted for (which can be done for example by using
the Chebyshev method) the final acos() will convert the "nearly zero"
sum to 0
This better implementation is based on the formula at
http://williams.best.vwh.net/avform.htm#Dist
sub great_circle_distance {
my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_;
$rho = 1 unless defined $rho; # Default to the unit sphere.
my $lat0 = pip2 - $phi0;
my $lat1 = pip2 - $phi1;
return $rho * 2 * asin(sqrt((sin(($lat0-$lat1)/2))**2 +
cos($lat0)*cos($lat1)*(sin(($theta1-$theta0)/2))**2));
}
Which gives the following (logically consistent) output for the example
above:
0.90871327210253
0.0908713255917341
0.0090871302470036
Confidentiality Notice: This message (including attachments) is a private communication solely for use of the intended recipient(s).
If you are not the intended recipient(s) or believe you received this message in error, notify the sender immediately and then delete this
message. Any other use, retention, dissemination or copying is prohibited and may be a violation of law, including the Electronic
Communication Privacy Act of 1986."