#!/usr/bin/perl
# calculate maximum latitude on Great Circle route
# by Ian Kluft, October 2003
#
# uses formulas from the Aviation Formulary site
#    http://williams.best.vwh.net/avform.htm
use strict;
use Math::Trig;
use POSIX qw( fmod );

# function: convert between degrees and radians
#sub deg2rad { return $_[0] * pi / 180; }
#sub rad2deg { return $_[0] * 180 / pi; }

# origin point: Newark NJ (EWR) Airport
my $lat1 = 40.6924972;
my $lon1 = -74.1686606;

# destination point: Amsterdam (AMS) Schiphol Airport
my $lat2 = 52.3086;
my $lon2 = 4.76389;

# compute initial true course
my $tc1=fmod(
	atan2(sin(deg2rad($lon1-$lon2))
		*cos(deg2rad($lat2)),
		cos(deg2rad($lat1))*sin(deg2rad($lat2))
	   	-sin(deg2rad($lat1))
		*cos(deg2rad($lat2))
		*cos(deg2rad($lon1-$lon2))),
	2*pi);

# compute maximum latitude along Great Cirle route
my $latmx = rad2deg(acos(abs(sin($tc1)*cos(deg2rad($lat1)))));

# display result
print "maximum latitude along Great Circle: $latmx\n";