#!/usr/bin/perl

# Calculate angle between 2 baselines given lat/long of 2 far
# locations an 1 center.

# Usage: bl.pl lat1 long1 lat2 long2 latc longc

use Math::Trig;

die "usage: bl.pl lat1 long1 lat2 long2.  (in deg)" if @ARGV != 6;

$PI = 3.14159265358979;

$lat1=$ARGV[0]*$PI/180.;
$long1=$ARGV[1]*$PI/180.;
$lat2=$ARGV[2]*$PI/180.;
$long2=$ARGV[3]*$PI/180.;
$latc=$ARGV[4]*$PI/180.;
$longc=$ARGV[5]*$PI/180.;

$rearth=6373;			# km

$x1 = $rearth * sin(0.5*$PI - $lat1) * cos($long1);
$x2 = $rearth * sin(0.5*$PI - $lat2) * cos($long2);
$xc = $rearth * sin(0.5*$PI - $latc) * cos($longc);

$y1 = $rearth * sin(0.5*$PI - $lat1) * sin($long1);
$y2 = $rearth * sin(0.5*$PI - $lat2) * sin($long2);
$yc = $rearth * sin(0.5*$PI - $latc) * sin($longc);

$z1 = $rearth * cos(0.5*$PI - $lat1);
$z2 = $rearth * cos(0.5*$PI - $lat2);
$zc = $rearth * cos(0.5*$PI - $latc);

$dx1=$x1-$xc;
$dy1=$y1-$yc;
$dz1=$z1-$zc;

$dx2=$x2-$xc;
$dy2=$y2-$yc;
$dz2=$z2-$zc;

$l1 = sqrt($dx1*$dx1 + $dy1*$dy1 + $dz1*$dz1);
$l2 = sqrt($dx2*$dx2 + $dy2*$dy2 + $dz2*$dz2);

$dot= $dx1*$dx2 + $dy1*$dy2 + $dz1*$dz2;
$dot /= $l1 * $l2;
$rad=acos($dot);
$deg=$rad/$PI*180.0;

print "bl1=$l1, bl2=$l2, angle=$rad ($deg deg)\n";