#!/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-projection - Calculate new coordinate from distance and bearing

=head1 SYNOPSIS

geo-projection [options] latitude longitude bearing distance[unit]

=head1 DESCRIPTION

geo-projection projects a waypoint by given start coordinates,
bearing and distance.

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.
The distance can have an optional suffix for the unit:
km, m, mi, in, ft and yd. km is the default.

=head1 OPTIONS

  -f|--formula f      Formula used for calculation. See GEO::Coords.
  -r|--radius r       Earth Radius for calculation. Default is 6371.64
  --degdec            Print only decimal format
  --degdec1           Print only decimal format with prefix
  --mindec            Print in DDD° MM.MMM format
  --mindec1           Shorter variant
  --mindec2           Print in DDD.MM.MMM format
  --utm               Print in UTM format
  --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;

#
# process command line arguments
#
use Getopt::Long;
my $help = 0;
my $man = 0;
my $version = 0;
my $usage = 0;
my $formula;
my $radius;
my $print_degdec = 0;
my $print_degdec1 = 0;
my $print_mindec = 0;
my $print_mindec1 = 0;
my $print_mindec2 = 0;
my $print_dms = 0;
my $print_utm = 0;

my $geo = new GEO::Coords;

GetOptions('r|radius=s' => \$radius,
	   'f|formula=s' => \$formula,
	   'degdec' => \$print_degdec,
           'degdec1' => \$print_degdec1,
           'mindec' => \$print_mindec,
           'mindec1' => \$print_mindec1,
           'mindec2' => \$print_mindec2,
           'dms' => \$print_dms,
           'utm' => \$print_utm,
	   '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 (geo-tools) 1.23\n";
  exit(0);
}

$print_mindec = !($print_degdec || $print_degdec1 || $print_mindec ||
                  $print_mindec1 ||$print_mindec2 || $print_dms ||
                  $print_utm);

#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, $bearing, $distance, $unit);

  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) {
      $bearing = $input;
    } elsif ($cnt == 4) {
      $input =~ /([\d\.]+)(\w*)/;
      $distance = $1;
      $unit = $2;
    } else {
      printf STDERR "Incomplete data\n";
      exit(1);
    }

    if ($cnt == 4) {
      $geo->unit($unit) if ($unit);

      my ($lat2,$lon2) = $geo->destination($lat1, $lon1, $bearing, $distance);

      printf "Formula: %s with earth radius: %s km\n",
	$geo->formula, $geo->earth_radius;

      GEO::Coords::print_decdeg($lat2, $lon2) if ($print_degdec);
      GEO::Coords::print_decdeg($lat2, $lon2, 1) if ($print_degdec1);
      GEO::Coords::print_mindec($lat2, $lon2, 0) if ($print_mindec);
      GEO::Coords::print_mindec($lat2, $lon2, 1) if ($print_mindec1);
      GEO::Coords::print_mindec($lat2, $lon2, 2) if ($print_mindec2);
      GEO::Coords::print_dms($lat2, $lon2, 0) if ($print_dms);
      GEO::Coords::print_utm($lat2, $lon2) if ($print_utm);

      $cnt = 0;
    }
  }
  if ($cnt != 0) {
    printf STDERR "Incomplete data\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;
}
