= Great Circle Algorithms =

!OpenLayers is a wonderful !JavaScript framework for map making, but its functions assume a !FlatLand view of the Earth. This is fine for simply displaying maps, but if you're trying to interact with your map by measuring distances or angles over large areas you're not likely to get accurate results.

The Earth is essentially a sphere, and so it's more appropriate to measure distances and angles using spherical trigonometry. First, note that the shortest distance between two points on a sphere is the route traced out by a plane passing through the center of the sphere and the two points where it intersects the sphere's surface. The arc described by the intersection of any plane that passes through the center of a sphere and the surface of the sphere is known as a great circle. All the meridians (lines of longitude) are great circles as is the Equator. For detailed discussion of great circles and their application see the [http://en.wikipedia.org/wiki/Great_circle/ Great Circle] page at !WikiPedia.

The !OpenLayers.Geometry.Point class can easily be extended to provide great circle functionality for calculating distances and angles between points as well as computing waypoints along a great circle route (i.e. finding the coordinates of a point that is a specified distance and bearing from an origin point.) The mathematics for the following functions comes from Ed Williams' excellent [http://williams.best.vwh.net/avform.htm Aviation Formulary] site.

If you want to download some working !JavaScript code that implements these algorithms in !OpenLayers it's available at [http://206.168.217.244/gallery/lib/Point.js Point.js]. This is used in an [http://206.168.217.244/gallery/LineAnimator.html OpenLayers demonstration of animating a line].

== Great Circle Distance ==

In this discussion (and in the actual functions), we assume that coordinates are Latitude and Longitude as signed decimal degrees with negative longitude in the western hemisphere and negative latitude in the southern hemisphere. Coordinates are also converted to radians to do the calculations as this simplifies the spherical trigonometry a great deal. 
{{{
Constants:
Radius of the Earth in miles
EARTH_RADIUS = 3958.75

Ratio of a circle's circumference to its diameter
PI = 3.1415926535897932384626433832795

Factor to convert decimal degrees to radians
DEG2RAD =  0.01745329252

Factor to convert decimal degrees to radians
RAD2DEG = 57.29577951308

Great Circle Distance:
Inputs are 2 coordinates x1,y1 and x2,y2

  Convert to radians
  x1 = x1 * DEG2RAD
  y1 = y1 * DEG2RAD
  x2 = x2 * DEG2RAD
  y2 = y2 * DEG2RAD
  
  a = sin(( y2-y1 ) / 2.0 )^2
  b = sin(( x2-x1 ) / 2.0 )^2
  c = sqrt( a + cos( y2 ) * cos( y1 ) * b )

  distance = 2 * asin( c ) * EARTH_RADIUS
}}}

== Great Circle Bearing ==

This algorithm calculates the angle that point 2 bears to point 1 along a great circle route. Note that as you move along that route from point 1 to point 2 the bearing is always changing.

{{{
Inputs are 2 coordinates x1,y1 and x2,y2

  Convert to radians
  x1 = x1 * DEG2RAD
  y1 = y1 * DEG2RAD
  x2 = x2 * DEG2RAD
  y2 = y2 * DEG2RAD
  
  a = cos(y2) * sin(x2 - x1)
  b = cos(y1) * sin(y2) - sin(y1) * cos(y2) * cos(x2 - x1)
  adjust = 0
    
  if((a == 0) && (b == 0)) {
      bearing = 0
  } else if( b == 0) {
      if( a < 0)  
          bearing = 3 * PI / 2
      else
          bearing = PI / 2
  } else if( b < 0) 
      adjust = PI
  else {
      if( a < 0) 
          adjust = 2 * PI
      else
          adjust = 0
  }
  bearing = (atan(a/b) + adjust) * RAD2DEG
}}}

== Great Circle Waypoints ==

This algorithm returns a point that is a specified distance and bearing from the input point.
Inputs are 2 coordinates x1,y1, bearing (an angle) and distance (units match EARTH_RADIUS). It returns coordinates x2, y2.
{{{
  Convert to radians
  x1 = x1 * DEG2RAD
  y1 = y1 * DEG2RAD
  bearing = bearing * DEG2RAD

  // Convert arc distance to radians
  c = distance / EARTH_RADIUS

  y2 = asin( sin(y1) * cos(c) + cos(y1) * sin(c) * cos(bearing)) * RAD2DEG


  a = sin(c) * sin(bearing)
  b = cos(y1) * cos(c) - sin(y1) * sin(c) * cos(bearing)

  if( b == 0 ) 
      x2 = x1
  else
      x2 = x1 + atan(a/b) * RAD2DEG

  return x2,y2
}}}