This program is a text companion to OeufMayo's drawmap.pl. It is no longer necessary to place this in the same directory as monks.xml
Summary: crowflight.plx is a simple script to calculate airline distance and bearing between map locations. The utility functions are based on my Cartesian 3-Vectors Snippet. Additional functions support handling of geographic data
Syntax:
$ ./crowflight.plx Zaxo # prints list of Monks reverse sorted by dist
+ance
baku,18308,62.8
greenFox,18161,-71.4
pope,16309,-17.2
.
.
notsoevil,281,-91.4
extremely,280,-92.2
jeffa,276,-92.2
Zaxo,0, 0.0
$ ./crowflight OeufMayo Corion #prints distance between these two Mon
+ks
Corion is 483 km from OeufMayo, bearing 70.5.
$ ./crowflight.plx Thelonious # no such Monk
Brother Thelonious has left the map
$
TODO
Add bearing. Partly doneNicer formatting DoneMetric/English cmd line switch (tip'o the hat to jcwren) DoneTidy up cl parsing (dimmesdale)Your suggestion here.
Please hammer on this, particularly wrt numerical stability. I'd appreciate suggestions for handling the bearing calculation. It seems to work, but there are untested corner cases.
After Compline,
Zaxo
Update: TODO suggestion by jcwren, added error case to screenshots. 6/22 added command line options for monkmap xml file location (-m=path) and English units (-u=[mi.*|E.*|e.*]). Added first cut at bearing calculation. Changed format to comma-separated fields, no filler for the one-monk output. 6/25 Another Todo sugg. by dimmesdale 6/26 normalize() picked up off the floor and inserted. Thanks, pope!
Da Code:
#!/usr/bin/perl -w
# crowflight.plx
# Geodesic code to get great circle distance between monk locations
# by Zaxo
# This code is intended as a contribution to the drawmap project
# Modified code from drawmap.pl gets monk locations from xml
# by Briac 'OeufMayo' Pilpré
# 2001/06/15
use strict;
sub usage {
print "Usage: ./crowflight.plx [-m <file>] [-u e|mi] Monkname [Mon
+kname]\n";
exit 0;
}
use XML::Simple;
use Getopt::Mixed 'nextOption';
Getopt::Mixed::init('m=s u=s xml>m units>u');
# Fetch the command line parameters
my ($xml,$units);
while (my ($option, $value, $pretty) = nextOption()) {
$xml = $value if $option eq 'm';
$units = $value if $option eq 'u';
}
Getopt::Mixed::cleanup();
my $monk = shift or die usage();
my $monk2 = shift || '';
if ($xml && $xml =~ /^(http:\/\/|ftp:\/\/)/) {
$xml = '' &&
warn "Warning: Network Monk file grab not yet implemented. Try
+ing local default file.";
}
$xml = './monks.xml' if !$xml;
die "No Monk map file found, quitting.\n" unless -f $xml;
$units = 'km' if ! $units;
$units = ($units =~ /^(e|mi)/i)? 'mi.': 'km';
my %monks;
# Parse the monks coordinates XML file I fetched from jcwren's stats s
+ite.
# ( code to fetch & create the XML is available on request )
my $xs = new XML::Simple();
my $ref = $xs->XMLin("./monks.xml");
# Fill the monks hash with their respective locations
foreach (keys %{$ref->{monk}}){
push (@{$monks{$_}}, (
$ref->{monk}->{$_}->{location}->{latitude},
$ref->{monk}->{$_}->{location}->{longitude},
));
}
# End OeufMayo's code
die "Brother $monk has left the map\n" unless $monks{$monk};
die "Brother $monk2 has left the map\n" unless !$monk2 || $monks{$monk
+2};
# From my Cartesian 3-Vector Snippet
# These functions all take references to arrays over [0..2] without
# modifying the arguments. Arguments assume Cartesian(x,y,z) geometry
# This could be expanded to a module.
# Cartesian 3-vector dot product
sub dot {
$_[0]->[0]*$_[1]->[0]+$_[0]->[1]*$_[1]->[1]+$_[0]->[2]*$_[1]->[2];
}
# Cartesian 3-vector cross product
# returns reference to array over [0..2]
sub cross {
[$_[0]->[1]*$_[1]->[2]-$_[0]->[2]*$_[1]->[1],
$_[0]->[2]*$_[1]->[0]-$_[0]->[0]*$_[1]->[2],
$_[0]->[0]*$_[1]->[1]-$_[0]->[1]*$_[1]->[0]];
}
# length of a Cartesian 3-vector
# returns number
sub norm {
sqrt(dot($_[0],$_[0]));
}
# angle between two Cartesian 3-vectors
# returns number: angle in radians
sub arc {
atan2 norm(cross($_[0],$_[1])), dot($_[0],$_[1]);
}
# End Cartesian 3-Vector Snippet code
# new 3-Vector code
sub normalize {
my $len = norm($_[0]);
($len)? [map {$_/$len} @{$_[0]}]
: [0,0,0];
}
use constant PI => 4 * atan2 1, 1;
my $EARTH_RADIUS = ($units eq 'mi.')? 3960: 6370;
# Convert latitude and longitude in decimal degrees
# to polar coordinates in radians. North Pole is "latitude" 0,
# South pole is "latitude" PI
sub geo2polar {
[ PI/2 - PI*$_[0]->[0]/180, PI*$_[0]->[1]/180];
}
# Map polar coords into a unit length Cartesian 3-Vector
# x-axis is near West Africa, y-axis in Indian Ocean,
# z-axis at North Pole
sub xyz1 {
[sin($_[0]->[0])*cos($_[0]->[1]),
sin($_[0]->[0])*sin($_[0]->[1]),
cos($_[0]->[0])];
}
# Great circle distance between given locations
# takes two references to array [latitude,longitude]
# returns number: surface distance in km.
sub distance {
$EARTH_RADIUS *
arc(xyz1(geo2polar($_[0])),xyz1(geo2polar($_[1])));
}
# Bearing of one location from another
# fails at the poles ( special case them? TODO: find out convention fo
+r that)
# and perhaps elsewhere.
# takes two references to array [latitude,longitude]
# returns number: signed north referenced bearing
sub bearing {
my ($here,$there) = ( xyz1(geo2polar($_[0])),
xyz1(geo2polar($_[1])));
my $east = normalize( cross([0,0,1], $here));
my $north = normalize( cross( $here, $east ));
my $beeline = normalize(cross($here, cross($there, $here)));
180 * (atan2 dot($beeline,$east), dot($beeline, $north))/PI;
}
if ($monk2) {
my $range = int(0.5+distance($monks{$monk},$monks{$monk2}));
my $bearing = bearing($monks{$monk},$monks{$monk2});
printf( "%s is %d %s from %s, bearing %5.1f.\n",
$monk2, $range, $units, $monk, $bearing);
}
else {
print
sort {
${[split ',', $b]}[-2] <=> ${[split ',', $a]}[-2]
}
map
{ sprintf( "%s,%d,%4.1f\n",
$_,
int(0.5 + distance($monks{$monk}, $monks{$_})),
bearing($monks{$monk}, $monks{$_}));}
keys %monks;
}
