I completed an initial test today, and found the results were not quite what I expected. I gave it 2 straight lines: the first true NS, the second true EW, both 1 degree long - since such distances can be calculated by hand, to compare results. The numbers given by your calculations were ~0.1 greater than expected. After much fiddling about, I realized that it was the radius of the earth that was used. According to
Wikipedia, the radius given is the equatorial radius, which is significantly larger than the polar radius. After numerous tests, I determined that a custom radius based on latitude would be more appropriate.
One way to do this would be to hardcode the radius based on the expected latitude. For the Galveston area, a latitude of about 29.5 would be acceptable. However, this would require adjusting the program based on what area was being measured. And what if the area being measured covered 60 degrees?
The other way to do it would be to calculate the radius on the spot based on the current latitude. So, first the program computes the average of the latitudes of the current set of points $lat0 & $lat1. Then, using the complex formula found at the above Wikipedia page, it computes a radius for that latitude, $r0. This is the radius that gets tossed into the great circle arc computation.
Well, it works. Here's a set of test data:
Code: Select all
>
-93.000000 17.000000
-93.000000 18.000000
>
-93.000000 27.000000
-93.000000 28.000000
>
-93.000000 37.000000
-93.000000 38.000000
>
-93.000000 47.000000
-93.000000 48.000000
>
And the results:
Code: Select all
Subtotal: 60.0894 nautical miles
Subtotal: 60.0646 nautical miles
Subtotal: 60.0328 nautical miles
Subtotal: 59.9979 nautical miles
Total: 240.1846 nautical miles
As you can see, the differences are somewhat significant, and they will tend to add up over an entire series. Granted, for the current application, where the latitudinal variation is very low, with values between 28.9788 and 29.7991, the difference ain't all that great.
Here's the current version of the program:
Code: Select all
#!/usr/bin/perl -w
use strict;
$/ = '>';
my $total = 0;
open(FILE,"<$ARGV[0]");
while(<FILE>) {
chomp();
next unless($_);
my $subtotal = 0;
my @list = split(/\n/);
for(0..$#list) {
next unless($list[$_]);
if ($list[$_+1]) {
my ($lon0, $lat0) = split(/\s+/, $list[$_]);
my ($lon1, $lat1) = split(/\s+/, $list[$_+1]);
$subtotal += calculate($lon1, (90-$lat1), $lon0, (90-$lat0));
}
}
printf("Subtotal: %.4f nautical miles\n",$subtotal) if $subtotal;
$total += $subtotal;
}
close(FILE);
printf("Total: %.4f nautical miles\n", $total);
sub calculate {
my ($lon0, $lat0, $lon1, $lat1) = @_;
my $pi = atan2(1, 1) * 4;
$lon0 *= ($pi/180);
$lat0 *= ($pi/180);
$lon1 *= ($pi/180);
$lat1 *= ($pi/180);
$lat0 = $pi/2 - $lat0;
$lat1 = $pi/2 - $lat1;
# Radius computation for geocentric latitude
my $r1 = 3443.917;
# equatorial radius in n-miles, use 3963.189 for s-miles or 6378.135 for km
my $r2 = 3432.37;
# polar radius in n-miles, use 3949.901 for s-miles or 6356.75 for km
my $avlat = ($lat0 + $lat1)/2;
my $r0 = ($r1 * $r2) / (sqrt($r1 * $r1 - (cos($avlat) * cos($avlat) * ($r1 * $r1 - $r2 * $r2))));
# Great circle arc distance calculation
my $x = cos($lat0) * cos($lat1) * cos($lon0 - $lon1) + sin($lat0) * sin($lat1);
return (atan2(sqrt(1 - $x**2), $x) * $r0);
}
