#!/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 [Monkname]\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. Trying 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 site.
# ( 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{$monk2};

# 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 for 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;
}