#!/usr/bin/perl

######################################################################
#
# distaz:
#
# Calculate the arc distance, azimuth and back azimuth for two pairs
# of geographic/geodetic coordinates.  The reported azimuth and back
# azimuth are with respect to the first coordinate pair.
#
# An enhanced version of an old "delaz" program, cobbled together by
# Chad Trabant, IRIS Data Management Center
######################################################################


use strict;
use Math::Trig;

my $version = "2006.341";

my $lat1  = undef;
my $lon1  = undef;
my $lat2  = undef;
my $lon2  = undef;

if ( $#ARGV == 3 ) {
  if ( $ARGV[0] <= 90.0 && $ARGV[0] >= -90.0 ) {
    $lat1 = $ARGV[0];
  } else {
    die "Latitude 1 specified incorrectly: $ARGV[0]\n";
  }
  
  if ( $ARGV[1] <= 180.0 && $ARGV[1] >= -180.0 ) {
    $lon1 = $ARGV[1];
  } else {
    die "Longitude 1 specified incorrectly: $ARGV[1]\n";
  }

  if ( $ARGV[2] <= 90.0 && $ARGV[2] >= -90.0 ) {
    $lat2 = $ARGV[2];
  } else {
    die "Latitude 2 specified incorrectly: $ARGV[2]\n";
  }
  
  if ( $ARGV[3] <= 180.0 && $ARGV[3] >= -180.0 ) {
    $lon2 = $ARGV[3];
  } else {
    die "Longitude 2 specified incorrectly: $ARGV[3]\n";
  }
  
  my ($delta, $azimuth, $backazimuth) = &distaz ($lat1, $lon1, $lat2, $lon2);
  
  # 111.32 km/deg and 69.17 miles/deg
  my $mdist = $delta * 111.32;
  my $edist = $delta * 69.17;
  
  printf "Distance: %g (deg), %g (km), %g (miles)\n", $delta, $mdist, $edist;
  printf "Azimuth: %g (deg), Backazimuth: %g (deg)\n", $azimuth, $backazimuth;
}
elsif ( $#ARGV == 0 ) {
  my $delta = $ARGV[0];
  
  # 111.32 km/deg and 69.17 miles/deg
  my $mdist = $delta * 111.32;
  my $edist = $delta * 69.17;
  
  printf "Distance: %g (deg), %g (km), %g (miles)\n", $delta, $mdist, $edist;
}
else {
  print "distaz ($version)\n";
  print "Compute arc distance, azimuth and back azimuth for Earth coordinates.\n";
  print "The arc distance is given as degrees, kilometers and miles.\n\n";
  print "Usage: distaz [lat1 lon1 lat2 lon2] | [delta]\n\n";
  print "For two coordinate pairs distance, azimuth and back azimuth will be printed.\n";
  print "For arc distance in degrees (delta) the distance will be printed.\n";
}

# Calculate delta (angular distance) and azimuth between two
# geographic/geodetic coordinates.  The returned azimuth and
# back-azimuth are with respect to the first coordinate pair.
# Latitudes are converted to geocentric latitudes using the WGS84
# spheriod to correct for ellipticity.
sub distaz {
  my ($lat1, $lon1, $lat2, $lon2) = @_;
  
  # Major and minor axies for WGS84 spheriod
  my $semimajor = 6378137.0;
  my $semiminor = 6356752.3142;
  
  my $ratio2 = (($semiminor * $semiminor) / ($semimajor * $semimajor));
  
  my $pirad  = acos(-1.0) / 180.0;
  my $halfpi = acos(-1.0) / 2;
  
  # Convert latitude to geocentric coordinates
  my $nlat1 = atan ($ratio2 * tan ($lat1 * $pirad));
  my $nlat2 = atan ($ratio2 * tan ($lat2 * $pirad));
  
  # Great circle calculation for delta and azimuth
  my $gamma = ( $lon2 - $lon1) * $pirad;
  my $a = ($halfpi - $nlat2);
  my $b = ($halfpi - $nlat1);
  
  my $sita;
  my $bsita;
  
  if ( $a == 0.0 ) {
    $sita = 1.0;
  }
  elsif ( $nlat2 == 0.0 ) {
    $sita = 0.0;
  }
  else {
    $sita = sin($b) / tan($a);
  }
  
  if ( $b == 0.0 ) {
    $bsita = 1.0;
  }
  elsif ( $nlat1 == 0.0 ) {
    $bsita = 0.0;
  }
  else {
    $bsita = sin($a) / tan($b);
  }
  
  my $delta = acos (cos($a) * cos($b) + sin($a) * sin($b) * cos($gamma)) / $pirad;
  
  my $azimuth = atan2 (sin($gamma), $sita - cos($gamma) * cos($b)) / $pirad;
  $azimuth = 0.0 if ( abs($azimuth) < 0.000001 );
  $azimuth += 360 if ( $azimuth < 0.0 );
  
  my $backazimuth = atan2 (-sin($gamma), $bsita - cos($gamma) * cos($a)) / $pirad;
  $backazimuth = 0.0 if ( abs($backazimuth) < 0.000001 );
  $backazimuth += 360 if ( $backazimuth < 0.0 );
  
  return ($delta, $azimuth, $backazimuth);
}