#!/usr/bin/perl

# Copyright (C) 2010 Thorsten Kukuk
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# in Version 2 as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
# MA  02110-1301, USA.

=head1 NAME

    geo-distance - Calculate distance and bearing between two coordinates

=head1 SYNOPSIS

    geo-distance [options] latitude1 longitude1 latitude2 longitude2

=head1 DESCRIPTION

    geo-distance calculates the distance between two points and
    the bearing.

    Acceptable input formats are (WGS84 Datum):

       Decimal       DegDec (decimal degrees)
                     49.4664 11.06903
                     "N 49.4664 E 11.06903"
                     "N 49.4664" "E 11.06903"

       DDD MM.MMM    MinDec (decimal minutes)
                     "N 49° 27.984" "E 011° 04.142"
                     "N 49° 27.984 E 011° 04.142"
                      N49.27.984 E011.04.142
                     "49 27.984" "011 04.142"
                      49.27.984 011.04.142

       DDD MM SS.SSS DMS (degrees, minutes, seconds)
                     "N 49° 27' 59.0\"" "E 11° 4' 08.5\""
                      "49 27 59.0"      "11 4 08.5"

       Instead of S and W a minus sign (-) can be used as prefix.


=head1 OPTIONS

    -d|--distance       Print only the distance value with unit.
    -f|--formula f      Formula used for calculation. See GEO::Coords.
    -r|--radius r       Earth Radius for calculation. Default is 6371.64
    -u|--unit unit      Use "unit" for the distance, default is km.
                        Possible values: km, m, mi, in, ft, yd
    --version           Print version number and exit
    --usage             Print short usage help
    -h|-?|--help	Help

=cut

use strict;
use warnings;
use Pod::Usage;
use GEO::Coords;
use Math::Trig;

my $Version = '(geo-tools) 1.23';

#
# process command line arguments
#
use Getopt::Long;
my $help = 0;
my $man = 0;
my $version = 0;
my $usage = 0;
my $unit;
my $dist_only = 0;
my $formula;
my $radius;

my $geo = new GEO::Coords;

GetOptions('d|distance' => \$dist_only,
	   'u|unit=s' => \$unit,
	   'r|radius=s' => \$radius,
	   'f|formula=s' => \$formula,
	   'version' => \$version,
	   'man' => \$man,
	   'usage' => \$usage,
	   'help|h|?' => \$help) or pod2usage(2);
pod2usage(0) if $help;
pod2usage(-exitstatus => 0, -verbose => 2) if $man;
pod2usage(-exitstatus => 0, -verbose => 0) if $usage;

if ($version) {
  print "geo-distance $Version\n";
  exit(0);
}

if ($unit) {
  chomp ($unit);
  $unit = $geo->unit($unit);
} else {
  $unit = $geo->unit;
}

if ($formula) {
  chomp ($formula);
  $geo->formula($formula);
}

$geo->earth_radius($radius) if ($radius);

if ($#ARGV >= 0) {
  parse_input(@ARGV);
} elsif (($#ARGV < 0) ||($#ARGV == 0 && $ARGV[0] eq "-")) {
  # read from tty
  while (<>) {
    my @Fields = SplitInputLine($_);
    parse_input(@Fields);
  }
}

exit 0;


sub parse_input
{
  my @Fields = @_;
  my $cnt = 0;
  my ($lat1, $lon1, $lat2, $lon2);

  foreach (@Fields) {
    my $input = $_;

    chomp ($input);
    next if (length($input) == 0);

    $cnt++;

    if ($cnt == 1) {
      ($lat1, $lon1) = GEO::Coords::parse2($input);

      if ($lat1 && $lon1) {
	$cnt++;
      } else {
	$lat1 = GEO::Coords::parse($input);
	if (!$lat1) {
	  print "Parse Error: $input\n";
	  exit(1);
	}
      }
    } elsif ($cnt == 2) {
      $lon1 = GEO::Coords::parse($input);
    } elsif ($cnt == 3) {
      ($lat2, $lon2) = GEO::Coords::parse2($input);

      if ($lat2 && $lon2) {
	$cnt++;
      } else {
	$lat2 = GEO::Coords::parse($input);
	if (!$lat2) {
	  print "Parse Error: $input\n";
	  exit(1);
	}
      }
    } elsif ($cnt == 4) {
      $lon2 = GEO::Coords::parse($input);

    }
    if ($cnt == 4) {
      my  $dist = $geo->distance($lat1, $lon1, $lat2, $lon2);

      if ($dist_only) {
	printf "%Lf%s\n", $dist, $unit;
      } else {
	printf "Formula: %s with earth radius: %s km\n",
	  $geo->formula, $geo->earth_radius;
	printf "%s -> %s\n", GEO::Coords::sprintf_mindec($lat1, $lon1),
	  GEO::Coords::sprintf_mindec($lat2, $lon2);
	printf "Distance: %Lf%s\n", $dist, $unit;
	my ($alpha1,$alpha2) = $geo->direction($lat1, $lon1, $lat2, $lon2);
	my ($alpha3,$alpha4) = $geo->direction($lat2, $lon2, $lat1, $lon1);
	if ($alpha2 && $alpha4) {
	  printf "1 -> 2: %3.2Lf° (Initial bearing/forward azimuth)\n",
	    $alpha1;
	  printf "1 -> 2: %3.2Lf° (Final bearing)\n",
	    $alpha2;
	  printf "2 -> 1: %3.2Lf° (Initial bearing/forward azimuth)\n",
	    $alpha3;
	  printf "2 -> 1: %3.2Lf° (Final bearing)\n",
	    $alpha4;
	} else {
	  printf "1 -> 2: %3.2Lf°\n", $alpha1;
	  printf "2 -> 1: %3.2Lf°\n", $alpha3;
	}
      }
     $cnt = 0;
    }
  }
  if ($cnt != 0) {
    printf STDERR "Incomplete coordinates\n";
    exit(1);
  }
}

sub SplitInputLine
{
   my $Line = $_[0];
   my @Fields = ();
   chomp($Line);
   if ($Line =~ /^\s*\#/) {
      return @Fields;
   }
   else {
      $Line =~ s/\s*\#.*//;
   }
   @Fields = split(/\s/, $Line);
   return @Fields;
}
