use Math::Trig 'asin';
use Math::Trig 'atan';
use Math::Trig 'acos';


sub station_distance {
	local($latitude_1, $longitude_1) = split('/', shift);
	local($latitude_2, $longitude_2) = split('/', shift);

	return(distance($latitude_1, $longitude_1, $latitude_2, $longitude_2));
}

#$PI =  3.1415926535897932384626433832795;
$PI =  3.14;

sub DEG2RAD {
	local($num) = @_;
	return($num * 3.1415926535897932384626433832795 / 180.0);	
}
sub RAD2DEG {
	local($num) = @_;
	return($num * 180.0 / 3.1415926535897932384626433832795);
}

sub SQR {
	local($num) = @_;
	return($num * $num);
}

sub distance {
	local($latitude_1, $longitude_1, $latitude_2, $longitude_2) = @_;

	if (($latitude_1 == $latitude_2) && ($longitude_1 == $longitude_2)) {
		return(0);
	}

	$cos_latitude1 = cos (DEG2RAD ($latitude_1));
	$sin_latitude1 = sin (DEG2RAD ($latitude_1));
	$cos_latitude2 = cos (DEG2RAD ($latitude_2));
	$sin_latitude2 = sin (DEG2RAD ($latitude_2));
 
	$sin_hdlatitude = sin (DEG2RAD (0.5 * ($latitude_1 - $latitude_2)));
	$sin_hdlongitude = sin (DEG2RAD (0.5 * ($longitude_1 - $longitude_2)));
	$sin_dlongitude = sin (DEG2RAD ($longitude_1 - $longitude_2));

	$sin_hdist_sq = SQR ($sin_hdlatitude) +
		$cos_latitude1 * $cos_latitude2 * SQR ($sin_hdlongitude);
	$distance_rad = 2.0 * asin (sqrt ($sin_hdist_sq));
	return(60.0 * RAD2DEG ($distance_rad));
}

sub course {
	local($latitude_1, $longitude_1, $latitude_2, $longitude_2) = @_;

	if (($latitude_1 == $latitude_2) && ($longitude_1 == $longitude_2)) {
		return(0);
	}

	$lat1 = DEG2RAD ($latitude_1);
	$lon1 = DEG2RAD ($longitude_1);
	$lat2 = DEG2RAD ($latitude_2);
	$lon2 = DEG2RAD ($longitude_2);
	$cos_latitude1 = cos (DEG2RAD ($lat1));
	$sin_latitude1 = sin (DEG2RAD ($lat1));
	$cos_latitude2 = cos (DEG2RAD ($lat2));
	$sin_latitude2 = sin (DEG2RAD ($lat2));
 
	$sin_hdlatitude = sin (DEG2RAD (0.5 * ($latitude_1 - $latitude_2)));
	$sin_hdlongitude = sin (DEG2RAD (0.5 * ($longitude_1 - $longitude_2)));
	$sin_dlongitude = sin (DEG2RAD ($longitude_1 - $longitude_2));

	$sin_hdist_sq = SQR ($sin_hdlatitude) +
		$cos_latitude1 * $cos_latitude2 * SQR ($sin_hdlongitude);
	$distance_rad = 2.0 * asin (sqrt ($sin_hdist_sq));
	
	#return(60.0 * RAD2DEG ($distance_rad));

	if (sin ($lon2 - $lon1) < 0) {
		$tc1 = acos((sin($lat2)-sin($lat1)*cos($distance_rad))/(sin($distance_rad)*cos($lat1))); 
	} else {
		$tc1 = 2*$PI-acos((sin($lat2)-sin($lat1)*cos($distance_rad))/(sin($distance_rad)*cos($lat1)));
	}
	
	return(360 - RAD2DEG($tc1));
}