package
{
	import flash.geom.Point // for return type

	/*
	excerpts from Ed Williams' Aviation Formulary:
	http://williams.best.vwh.net/avform.htm#Dist

	The great circle ^distance d between two points with coordinates {lat1,lon1} and
	{lat2,lon2} is given by: 

	d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2))

	A mathematically equivalent formula, which is less subject to rounding error for
	short ^distances is:

	d=2*asin(sqrt((sin((lat1-lat2)/2))^2 + cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))

	http://williams.best.vwh.net/avform.htm#Intermediate

	Here we find points (lat,lon) a given fraction of the ^distance (d) between them.
	Suppose the starting point is (lat1,lon1) and the final point (lat2,lon2) and we
	want the point a fraction f along the great circle route. f=0 is point 1. f=1 is
	point 2. The two points cannot be antipodal ( i.e. lat1+lat2=0 and abs(lon1-lon2)=pi)
	because then the route is undefined. The intermediate latitude and longitude is then
	given by: 

        A=sin((1-f)*d)/sin(d)
        B=sin(f*d)/sin(d)
        x = A*cos(lat1)*cos(lon1) +  B*cos(lat2)*cos(lon2)
        y = A*cos(lat1)*sin(lon1) +  B*cos(lat2)*sin(lon2)
        z = A*sin(lat1)           +  B*sin(lat2)
        lat=atan2(z,sqrt(x^2+y^2))
        lon=atan2(y,x)

	^distance means "arc length" (it is actually distance for R = 1 only)
	see eq. 5 at http://mathworld.wolfram.com/GreatCircle.html
	
	(thanks to Jeff Whitaker <jeffrey.s.whitaker@noaa.gov> for his python code with pointers)
	*/

	public class GreatCircle
	{
		private var lon1:Number, lat1:Number, lon2:Number, lat2:Number; // start and end points (in radians)
		private var d:Number; // arc length (in radians)

		public function GreatCircle (p_nLon1:Number, p_nLat1:Number, p_nLon2:Number, p_nLat2:Number)
		{
			lon1 = p_nLon1; lat1 = p_nLat1;
			lon2 = p_nLon2; lat2 = p_nLat2;
			var sinLat12:Number = Math.sin ((lat1 - lat2) / 2);
			var sinLon12:Number = Math.sin ((lon1 - lon2) / 2);
			d = 2 * Math.asin (Math.sqrt (
				sinLat12 * sinLat12 +
				sinLon12 * sinLon12 * Math.cos (lat1) * Math.cos (lat2)
			));
			if (d == Math.PI)
			{
				trace ('throw something here');
			}
		}

		public function getPointAt (f:Number)
		{
			var A:Number = Math.sin ((1-f)* d) / Math.sin (d);
			var B:Number = Math.sin     (f* d) / Math.sin (d);
			var x:Number = A * Math.cos (lat1) * Math.cos (lon1) + B * Math.cos (lat2) * Math.cos(lon2);
			var y:Number = A * Math.cos (lat1) * Math.sin (lon1) + B * Math.cos (lat2) * Math.sin(lon2);
			var z:Number = A * Math.sin (lat1)                   + B * Math.sin (lat2);
			return new Point (
				/* lon(f) */ Math.atan2 (y, x),
				/* lat(f) */ Math.atan2 (z, Math.sqrt (x * x + y * y))
			);
		}
	}
}