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 Tshirt 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 Tshirt 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;
}

RE: Lat/Long distance calculator by Fastolfe (Vicar) on Nov 07, 2000 at 20:59 UTC 
Also see Math::Trig's great_circle_distance function. From the documentation:
# To calculate the distance between London (51.3N 0.5W)
# and Tokyo (35.7N 139.8E) in kilometers:
use Math::Trig qw(great_circle_distance deg2rad);
# Notice the 90  latitude: phi zero is at the North Pole.
@L = (deg2rad(0.5), deg2rad(90  51.3));
@T = (deg2rad(139.8),deg2rad(90  35.7));
$km = great_circle_distance(@L, @T, 6378);
 [reply] [d/l] [select] 

Typical. I searched CPAN for latitude, longitude, Geo, and a few others, but not Math::Trig. Hardly an intuitive place to look.
Chris
email jcwren
 [reply] 

I agree. I only happened to remember it because I had examined Math::Trig for unrelated functions in the past and noticed it was in there.
 [reply] 
RE: Lat/Long distance calculator by chromatic (Archbishop) on Nov 08, 2000 at 10:02 UTC 
 [reply] 

 [reply] 
Back to Code Catacombs

