#!/usr/bin/perl -w
# NAME: distance.pl
# AIM: Use perl trib function for distance London to Tokyo...
# from : http://perldoc.perl.org/Math/Trig.html
use strict;
use warnings;
use Math::Trig qw(great_circle_distance great_circle_direction deg2rad rad2deg);
require 'fg_wsg84.pl' or die "Unable to load fg_wsg84.pl ...\n";
 
# 1 knots = 1.85200 kph
my $K2KPH = 1.85200;

my $lonlon = -0.5;
my $lonlat = 51.3;
my $toklon = 139.8;
my $toklat = 35.7;
my $g_ias = 80; # Knots - c152-c172
my $g_speed = $g_ias * $K2KPH; # Kilometers/Hour
my $CP_EPSILON = 0.00001;

 # Notice the 90 - latitude: phi zero is at the North Pole
sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }

sub show_sg_distance_vs_est($$$$$$) {
    my ($lat1,$lat2,$lon1,$lon2,$estdist,$esthdg) = @_;
    my ($sg_az1,$sg_az2,$sg_dist,$res,$sg_km,$sg_ikm);
    my ($eddiff,$edpc,$ahdg,$ahdiff,$hdgpc,$adhdg,$distmsg,$hdgmsg);
    $res = fg_geo_inverse_wgs_84 ($lat1,$lon1,$lat2,$lon2,\$sg_az1,\$sg_az2,\$sg_dist);
    $sg_km = $sg_dist / 1000;
    $sg_ikm = int($sg_km + 0.5);
    $sg_az1 = int(($sg_az1 * 10) + 0.05) / 10;
    $sg_az2 = int(($sg_az2 * 10) + 0.05) / 10;
    $eddiff = abs($sg_km - $estdist);
    $edpc = ($eddiff / $sg_km) * 100;
    $edpc = int(($edpc * 10) + 0.05) / 10;
    if ($eddiff < $CP_EPSILON) {
        $edpc = "SAME";
    } elsif ($edpc < 1) {
        $edpc = 'lt 1';
    }
    $edpc = " $edpc" while (length($edpc) < 5);

    $adhdg = abs($sg_az1 - $esthdg);
    $ahdiff = ($adhdg / $sg_az1);
    $hdgpc = ($ahdiff / 360);
    $hdgpc = int(($hdgpc * 10) + 0.05) / 10;
    if ($ahdiff < $CP_EPSILON) {
        $hdgpc = "SAME";
    } elsif ($hdgpc < 1) {
        $hdgpc = 'lt 1';
    }
    $hdgpc = " $hdgpc" while (length($hdgpc) < 5);
    print "SG_Dist : $sg_ikm kilometers ($sg_km) ($edpc \% diff)\n";
    print "SG_Head : $sg_az1 degs (inverse $sg_az2) ($hdgpc \% diff)\n";
}

sub show_distance($$$$) {
    my ($lon1,$lat1,$lon2,$lat2) = @_;
    my @Pos1 = NESW( $lon1, $lat1);
    my @Pos2 = NESW( $lon2, $lat2);
    my $d_km = great_circle_distance(@Pos1, @Pos2, 6378); # About 9600 km.
    my $t_degs = rad2deg(great_circle_distance(@Pos1, @Pos2));
    my $hdg = rad2deg(great_circle_direction(@Pos1, @Pos2));
    $hdg = int(($hdg * 10) + 0.05) / 10;
    my $hrs = $d_km / $g_speed;
    my $eta = get_hhmmss($hrs);
    my $ikm = int($d_km + 0.5);
    #$degs = int($degs + 0.5);
    #$degs = int($degs + 0.5);
    $t_degs = sprintf("%0.6f",$t_degs);
    print "From (lon,lat) $lon1,$lon2 to $lon2,$lat2 is about -\n";
    print "Distance: $ikm kilometers ($d_km)\n";
    print "Heading : $hdg, for $t_degs degs\n";
    print "ETA     : $eta, at $g_ias Knots\n";
    # print "km $km / spd $g_speed = $hours\n";
    show_sg_distance_vs_est($lat1,$lat2,$lon1,$lon2,$d_km,$t_degs);
    # =================================================
}

sub get_hhmmss($) {
    my ($hours) = @_;
    my $hrs = int($hours);
    my $mins = ($hours * 60) - ($hrs * 60);
    my $min = int($mins);
    my $secs = ($mins * 60) - ($min * 60);
    my $days = 0;
    if ($hrs >= 24) {
        $days++;
        $hrs -= 24;
    }
    $min = ($min < 10) ? "0$min" : $min;
    $secs = int( $secs + 0.5 );
    $secs = ($secs < 10) ? "0$secs" : $secs;
    $hrs = ($hrs < 10) ? "0$hrs" : $hrs;
    my $stg = '';
    $stg .= "Days $days, " if ($days);
    $stg .= "$hrs:$min:$secs";
    return $stg;
}

my ($rkm,$heading,@L,@T,$km,$degs,$argc,$hours,$hms);
$argc = scalar @ARGV;
if ($argc == 4) {
    $lonlon = $ARGV[0];
    $lonlat = $ARGV[1];
    $toklon = $ARGV[2];
    $toklat = $ARGV[3];
    @L = NESW( $lonlon, $lonlat);
    @T = NESW( $toklon, $toklat);
    $km = great_circle_distance(@L, @T, 6378); # About 9600 km.
    $degs = rad2deg(great_circle_distance(@L, @T));
    $heading = rad2deg(great_circle_direction(@L, @T));
    $heading = int(($heading * 10) + 0.05) / 10;
    $hours = $km / $g_speed;
    $hms = get_hhmmss($hours);
    $rkm = int($km + 0.5);
    #$degs = int($degs + 0.5);
    #$degs = int($degs + 0.5);
    $degs = sprintf("%0.6f",$degs);
    print "From (lon,lat) $lonlon,$lonlat to $toklon,$toklat is about -\n";
    print "Distance: $rkm kilometers ($km)\n";
    print "Heading : $heading, for $degs degs\n";
    print "ETA     : $hms, at $g_ias Knots\n";
    # print "km $km / spd $g_speed = $hours\n";
} else {
    print "ERROR: Need 4 arguments, not $argc!\n" if ($argc);
    print "Given two pairs of longitude, latitude, and will calculate the distance, and heading...\n";
    show_distance( $lonlon, $lonlat, $toklon, $toklat );
}

# eof - distance.pl
