#!/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-intersection - Calculate intersection of two lines

=head1 SYNOPSIS

geo-distance [options] lat1 long1 lat2 long2 lat3 long3 lat4 long4

=head1 DESCRIPTION

geo-distance calculates the intersection of two lines given
by four waypoints.

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"

   UTM ZONE X Y
                 "32U 649911 5481363"

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


=head1 OPTIONS

  --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
  -l|--lat        Output Latitude only
  -L|--long       Output Longitude only
  --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 $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 $printall = 1;
my $lat_only = 0;
my $long_only = 0;

my $geo = new GEO::Coords;

GetOptions('d|degdec' => \$print_degdec,
           'degdec1' => \$print_degdec1,
           'mindec' => \$print_mindec,
           'mindec1' => \$print_mindec1,
           'm|mindec2' => \$print_mindec2,
           'dms' => \$print_dms,
           'utm' => \$print_utm,
           'l|lat' => \$lat_only,
           'L|long' => \$long_only,
	   '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-intersection (geo-tools) 1.23\n";
  exit(0);
}

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

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, $lat3, $lon3, $lat4, $lon4);

  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);
    } elsif ($cnt == 5) {
      ($lat3, $lon3) = GEO::Coords::parse2($input);

      if ($lat3 && $lon3) {
	$cnt++;
      } else {
	$lat3 = GEO::Coords::parse($input);
	if (!$lat3) {
	  print "Parse Error: $input\n";
	  exit(1);
	}
      }
    } elsif ($cnt == 6) {
      $lon3 = GEO::Coords::parse($input);
    } elsif ($cnt == 7) {
      ($lat4, $lon4) = GEO::Coords::parse2($input);

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

    if ($cnt == 8) {
      my ($lat, $lon) = $geo->intersection ($lat1, $lon1, $lat2, $lon2, $lat3, $lon3, $lat4, $lon4);

      if ($lat && $lon) {
	  print_coordinates ($lat, $lon);
	} else {
	  print "Lines do not cross\n";
	}
      $cnt = 0;
    }
  }

  if ($cnt != 0) {
    printf STDERR "Incomplete coordinates\n";
    exit(1);
  }
}

sub print_coordinates
{
  my $dec1 = $_[0];
  my $dec2 = $_[1];

  GEO::Coords::print_decdeg($dec1, $dec2, 0,
                            $lat_only, $long_only) if ($printall || $print_degdec);
  GEO::Coords::print_decdeg($dec1, $dec2, 1,
                            $lat_only, $long_only) if ($printall || $print_degdec1);
  GEO::Coords::print_mindec($dec1, $dec2, 0,
                            $lat_only, $long_only) if ($printall || $print_mindec);
  GEO::Coords::print_mindec($dec1, $dec2, 1,
                            $lat_only, $long_only) if ($printall || $print_mindec1);
  GEO::Coords::print_mindec($dec1, $dec2, 2,
                            $lat_only, $long_only) if ($printall || $print_mindec2);
  GEO::Coords::print_dms($dec1, $dec2, 0) if ($printall || $print_dms);
  GEO::Coords::print_utm($dec1, $dec2) if ($printall || $print_utm);
}


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