Category: | Miscellaneous |
Author/Contact Info | jcwren@jcwren.com |
Description: | This code takes a latitude/longtitude pair and returns the distance and heading between them. As a practical (hah!) demonstration, it calculates how far your PM T-shirt has to travel to get from vroom to you (assuming you change the location, otherwise it's vroom to jcwren) |
#!/usr/local/bin/perl -w use strict; use POSIX qw(atan acos floor); # # Use www.mapblast.com, or your coordinates from the MonkMap page if +you're already there. # Latitude and longtitude are expressed in decimal notation. North a +nd east values are # positive, South and west are negative. # { my @location_of_me = qw(34.172500 -84.001667); # latitude=34.17 +2500, longitude=-84.001667 of jcwren my @location_of_vroom = qw(42.763333 -86.110556); # latitude=42.76 +3333, longitude=-86.110556 of vroom print sprintf ("\nMy T-shirt has to travel %0.f miles, on a heading + of %.1f degrees, to get from vroom to me\n\n", range_and_bearing (@location_of_vroom, @location_of_ +me)); } # # Returns range (in miles), and angle (relative to $latitude1/$longti +tude1) # # Uses basic great circle route calculation. Output compares favorab +ly with www.mapblast.com # www.mapquest.com, and Delorme Street Atlas. Minor deviations since + the earth is not a true # sphere, but bulges slighty at the equator. Accurate enough for com +mon distance calculations, # do not use to calculate fuel usage for cruise missiles targetting ( +0.2% error seems to be # typical worst case, though) # sub range_and_bearing { my ($latitude1, $longitude1, $latitude2, $longitude2) = @_; my $pi = 3.1415926535897; my $radToDeg = 180.0 / $pi; my $degToRad = $pi / 180.0; my $earthRadius = 3958.9; # Use 3958.9=miles, 6371.0=K +m; my $distance = 0; my $azimuth = 0; my $beta = 0; my $cosBeta = 0; my $cosAzimuth = 0; $latitude1 = $latitude1 * $degToRad; $longitude1 = $longitude1 * $degToRad; $latitude2 = $latitude2 * $degToRad; $longitude2 = $longitude2 * $degToRad; if (abs ($latitude1) < 90.0) { $cosBeta = (sin ($latitude1) * sin ($latitude2)) + ((cos ($latit +ude1) * cos ($latitude2)) * cos ($longitude2 - $longitude1)); if ($cosBeta >= 1.0) { return (0.0, 0.0); } # # Antipodes (return miles, 0 degrees) # if ($cosBeta <= -1.0) { return (floor ($earthRadius * $pi * 100.0) / 100.0, 0.0); } $beta = acos ($cosBeta); $distance = $beta * $earthRadius; $cosAzimuth = (sin ($latitude2) - sin ($latitude1) * cos ($beta) +) / (cos ($latitude1) * sin ($beta)); if ($cosAzimuth >= 1.0) { $azimuth = 0.0; } elsif ($cosAzimuth <= -1.0) { $azimuth = 180.0; } else { $azimuth = acos ($cosAzimuth) * $radToDeg; } if (sin ($longitude2 - $longitude1) < 0.0) { $azimuth = 360.0 - $azimuth; } return (floor ($distance * 100.0) / 100.0, floor ($azimuth * 100 +.0) / 100.0); } # # If P1 Is North Or South Pole, Then Azimuth Is Undefined # if (sgn ($latitude1) == sgn ($latitude2)) { $distance = $earthRadius * ($pi / 2 - abs ($latitude2)); } else { $distance = $earthRadius * ($pi / 2 + abs ($latitude2)); } return (floor ($distance * 100.0) / 100.0, 0.0); } # # Why is there no intrinsic sign function in Perl? # sub sgn { return $_[0] == 0 ? 0 : $_[0] < 0 ? -1 : 1; } |
|
---|
Replies are listed 'Best First'. | |
---|---|
RE: Lat/Long distance calculator
by Fastolfe (Vicar) on Nov 07, 2000 at 20:59 UTC | |
by jcwren (Prior) on Nov 07, 2000 at 21:00 UTC | |
by Fastolfe (Vicar) on Nov 07, 2000 at 21:02 UTC | |
RE: Lat/Long distance calculator
by chromatic (Archbishop) on Nov 08, 2000 at 10:02 UTC | |
by extremely (Priest) on Nov 08, 2000 at 13:04 UTC |
Back to
Code Catacombs