
/*
 * LatLng object - methods summary
 *
 *   p = new LatLng('512839N', '0002741W')
 *   p = new LatLng(53.123, -1.987)
 *
 *   dist = LatLng.distHaversine(p1, p2)
 *   dist = LatLng.distCosineLaw(p1, p2)
 *   dist = LatLng.distVincenty(p1, p2)
 *
 *   brng = LatLng.bearing(p1, p2)
 *   dist = p1.distAlongVector(orig, dirn)
 *   p = LatLng.midPoint(p1, p2)
 *   p2 = p1.destPoint(initBrng, dist)
 *   brng = p.finalBrng(initBrng, dist)
 *
 *   dist = LatLng.distRhumb(p1, p2)
 *   brng = LatLng.brngRhumb(p1, p2)
 *   p2 = p1.destPointRhumb(brng, dist)
 *
 *   rad = LatLng.llToRad('51º28'39"N')
 *   latDms = p.latitude()
 *   lngDms = p.longitude()
 *   dms = LatLng.radToDegMinSec(0.1284563)
 *   dms = LatLng.radToBrng(0.1284563)
 *
 * properties:
 *   p.lat - latitude in radians (0=equator, pi/2=N.pole)
 *   p.lng - longitude in radians (0=Greenwich, E=+ve)
 *
 * © 2002-2005 Chris Veness, www.movable-type.co.uk
 */


/*
 * LatLng constructor:
 *
 *	arguments are in degrees: signed decimal or d-m-s + NSEW as per LatLng.llToRad()
 *
*/

function LatLng(degLat, degLng) 
{
  this.lat = LatLng.llToRad(degLat);
  this.lng = LatLng.llToRad(degLng);
}


/*
 * Calculate distance (in km) between two points specified by latitude/longitude with Haversine formula
 *
 * from: Haversine formula - R. W. Sinnott, "Virtues of the Haversine",
 *       Sky and Telescope, vol 68, no 2, 1984
 *       http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
 */
LatLng.distHaversine = function(p1, p2) 
{
  var R = 6371; // earth's mean radius in km
  var dLat  = p2.lat - p1.lat;
  var dLng = p2.lng - p1.lng;

  var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
          Math.cos(p1.lat) * Math.cos(p2.lat) * Math.sin(dLng/2) * Math.sin(dLng/2);
  var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
  var d = R * c;

  return d;
}


/*
 * Calculate distance (in km) between two points specified by latitude/longitude using law of cosines.
 */
LatLng.distCosineLaw = function(p1, p2) 
{
  var R = 6371; // earth's mean radius in km
  var d = Math.acos(Math.sin(p1.lat)*Math.sin(p2.lat) +
                    Math.cos(p1.lat)*Math.cos(p2.lat)*Math.cos(p2.lng-p1.lng)) * R;
  return d;
}


/*
 * calculate (initial) bearing (in radians clockwise) between two points
 *
 * from: Ed Williams' Aviation Formulary, http://williams.best.vwh.net/avform.htm#Crs
 */
LatLng.bearing = function(p1, p2) 
{
  var y = Math.sin(p2.lng-p1.lng) * Math.cos(p2.lat);
  var x = Math.cos(p1.lat)*Math.sin(p2.lat) -
          Math.sin(p1.lat)*Math.cos(p2.lat)*Math.cos(p2.lng-p1.lng);
  return Math.atan2(y, x);
}


/*
 * calculate distance of point along a given vector defined by origin point
 * and direction in radians (uses planar not spherical geometry, so only valid
 * for small distances).
 */
LatLng.prototype.distAlongVector = function(orig, dirn) 
{
  var dist = LatLng.distHaversine(this, orig);  // distance from orig to point
  var brng = LatLng.bearing(this, orig);        // bearing between orig and point
  return dist * Math.cos(brng-dirn);
}


/*
 * calculate midpoint of great circle line between p1 & p2.
 *   see http://mathforum.org/library/drmath/view/51822.html for derivation
 */
LatLng.midPoint = function(p1, p2) 
{
  var dLng = p2.lng - p1.lng;

  var Bx = Math.cos(p2.lat) * Math.cos(dLng);
  var By = Math.cos(p2.lat) * Math.sin(dLng);

  lat3 = Math.atan2(Math.sin(p1.lat)+Math.sin(p2.lat),
                    Math.sqrt((Math.cos(p1.lat)+Bx)*(Math.cos(p1.lat)+Bx) + By*By ) );
  lng3 = p1.lng + Math.atan2(By, Math.cos(p1.lat) + Bx);

  if (isNaN(lat3) || isNaN(lng3)) return null;
  return new LatLng(lat3*180/Math.PI, lng3*180/Math.PI);
}


/*
 * calculate destination point given start point, initial bearing and distance
 *   see http://williams.best.vwh.net/avform.htm#LL
 */
LatLng.prototype.destPoint = function(brng, dist, miles) 
{
	if (!miles) var R = 6371; // earth's mean radius in km
	else var R = 3963;// earth's mean radius in miles // 6371; // earth's mean radius in km
  var p1 = this, p2 = new LatLng(0,0), d = parseFloat(dist)/R;  // d = angular distance covered on earth's surface

	brng = LatLng.degToRad(brng);

  p2.lat = Math.asin( Math.sin(p1.lat)*Math.cos(d) + Math.cos(p1.lat)*Math.sin(d)*Math.cos(brng) );
  p2.lng = p1.lng + Math.atan2(Math.sin(brng)*Math.sin(d)*Math.cos(p1.lat), Math.cos(d)-Math.sin(p1.lat)*Math.sin(p2.lat));

	if (isNaN(p2.lat) || isNaN(p2.lng)) return null;
//  p2 = llDmsToDec(p2);alert(p2);
  return p2;
}


/*
 * calculate final bearing arriving at destination point given start point, initial bearing and distance
 */
LatLng.prototype.finalBrng = function(brng, dist) 
{
  var p1 = this, p2 = p1.destPoint(brng, dist);
  // get reverse bearing point 2 to point 1 & reverse it by adding 180º
  var h2 = (LatLng.bearing(p2, p1) + Math.PI) % (2*Math.PI);
  return h2;
}


/*
 * calculate distance, bearing, destination point on rhumb line
 *   see http://williams.best.vwh.net/avform.htm#Rhumb
 */
LatLng.distRhumb = function(p1, p2) 
{
  var R = 6371; // earth's mean radius in km
  var dLat = p2.lat-p1.lat, dLng = Math.abs(p2.lng-p1.lng);
  var dPhi = Math.log(Math.tan(p2.lat/2+Math.PI/4)/Math.tan(p1.lat/2+Math.PI/4));
  var q = dLat/dPhi;
  if (!isFinite(q)) q = Math.cos(p1.lat);
  // if dLng over 180° take shorter rhumb across 180° meridian:
  if (dLng > Math.PI) dLng = 2*Math.PI - dLng;
  var d = Math.sqrt(dLat*dLat + q*q*dLng*dLng); 
  return d * R;
}


LatLng.brngRhumb = function(p1, p2) 
{
  var dLng = p2.lng-p1.lng;
  var dPhi = Math.log(Math.tan(p2.lat/2+Math.PI/4)/Math.tan(p1.lat/2+Math.PI/4));
  if (Math.abs(dLng) > Math.PI) dLng = dLng>0 ? -(2*Math.PI-dLng) : (2*Math.PI+dLng);
  return Math.atan2(dLng, dPhi);
}


LatLng.prototype.destPointRhumb = function(brng, dist) 
{
  var R = 6371; // earth's mean radius in km
  var p1 = this, p2 = new LatLng(0,0);
  var d = parseFloat(dist)/R;  // d = angular distance covered on earth's surface

  brng = LatLng.degToRad(brng);

  p2.lat = p1.lat + d*Math.cos(brng);
  var dPhi = Math.log(Math.tan(p2.lat/2+Math.PI/4)/Math.tan(p1.lat/2+Math.PI/4));
  var q = (p2.lat-p1.lat)/dPhi;
  if (!isFinite(q)) q = Math.cos(p1.lat);
  var dLng = d*Math.sin(brng)/q;
  // check for some daft bugger going past the pole
  if (Math.abs(p2.lat) > Math.PI/2) p2.lat = p2.lat>0 ? Math.PI-p2.lat : -Math.PI-p2.lat;
  p2.lng = (p1.lng+dLng+Math.PI)%(2*Math.PI) - Math.PI;
 
  if (isNaN(p2.lat) || isNaN(p2.lng)) return null;
  return p2;
}


/*
 * convert lat/long in degrees to radians, for handling input values
 *
 *   this is very flexible on formats, allowing signed decimal degrees (numeric or text), or
 *   deg-min-sec suffixed by compass direction (NSEW). A variety of separators are accepted 
 *   (eg 3º 37' 09"W) or fixed-width format without separators (eg 0033709W). Seconds and minutes
 *   may be omitted. Minimal validation is done.
 *
 */
/*
LatLng.llToBrng = function(brng) 
{alert(brng);
	return brng;
}
 */


/*
 * convert lat/long in degrees to radians, for handling input values
 *
 *   this is very flexible on formats, allowing signed decimal degrees (numeric or text), or
 *   deg-min-sec suffixed by compass direction (NSEW). A variety of separators are accepted 
 *   (eg 3º 37' 09"W) or fixed-width format without separators (eg 0033709W). Seconds and minutes
 *   may be omitted. Minimal validation is done.
 */
LatLng.llToRad = function(brng) 
{//alert(brng);
  if (!isNaN(brng)) return brng * Math.PI / 180;  // signed decimal degrees without NSEW

  brng = brng.replace(/[\s]*$/,'');               // strip trailing whitespace
  var dir = brng.slice(-1).toUpperCase();         // compass dir'n
  if (!/[NSEW]/.test(dir)) return NaN;            // check for correct compass direction
  brng = brng.slice(0,-1);                        // and lose it off the end
//  var dms = brng.split(/[\s:,°º?\'?\"]/);         // check for separators indicating d/m/s
	var dms = brng;
  switch (dms.length) 
	{                           // convert to decimal degrees...
    case 3:                                       // interpret 3-part result as d/m/s
      var deg = dms[0]/1 + dms[1]/60 + dms[2]/3600; break;
    case 2:                                       // interpret 2-part result as d/m
      var deg = dms[0]/1 + dms[1]/60; break;
    case 1:                                       // non-separated format dddmmss
      if (/[NS]/.test(dir)) brng = '0' + brng;    // - normalise N/S to 3-digit degrees
      var deg = brng.slice(0,3)/1 + brng.slice(3,5)/60 + brng.slice(5)/3600; break;
    default: return NaN;
  }
  if (/[WS]/.test(dir)) deg = -deg;               // take west and south as -ve
  return deg * Math.PI / 180;                     // then convert to radians
}


/* 
 * convert degrees to radians - used for bearing, so 360º with no N/S/E/W suffix
 *   can accept d/m/s, d/m, or decimal degrees
 */
LatLng.degToRad = function(brng) 
{//alert(brng);
//  var dms = brng.split(/[\s:,º°\'\"??]/);          // check for separators indicating d/m/s (PB added semi-colon)
//	var dms = brng.split(/[\s:,°º?\'?\"]/);         // check for separators indicating d/m/s (PB added semi-colon)
	var dms = brng;
  switch (dms.length) 
	{                           // convert to decimal degrees...
    case 3:                                       // interpret 3-part result as d/m/s
      var deg = dms[0]/1 + dms[1]/60 + dms[2]/3600; break;
    case 2:                                       // interpret 2-part result as d/m
      var deg = dms[0]/1 + dms[1]/60; break;
    default: 
      var deg = parseFloat(brng); break;          // otherwise decimal degrees
  }
  return deg * Math.PI / 180;                     // then convert to radians
}


/*
 * convert latitude into degrees, minutes, seconds; eg 51º28'38"N
 */
LatLng.prototype.latitude = function() 
{
  return LatLng._dms(this.lat).slice(1) + (this.lat<0 ? 'S' : 'N');
}


/*
 * convert longitude into degrees, minutes, seconds; eg 000º27'41"W
 */
LatLng.prototype.longitude = function() 
{
  return LatLng._dms(this.lng) + (this.lng>0 ? 'E' : 'W');
}


/*
 * convert radians to (signed) degrees, minutes, seconds; eg -0.1rad = -000°05'44"
 */
LatLng.radToDegMinSec = function(rad) 
{
  return (rad<0?'-':'') + LatLng._dms(rad);
}


/*
 * convert radians to compass bearing - 0°-360° rather than +ve/-ve
 */
LatLng.radToBrng = function(rad) 
{
  return LatLng.radToDegMinSec((rad+2*Math.PI) % (2*Math.PI));
}


/*
 * convert radians to deg/min/sec, with no sign or compass dirn (internal use)
 */
LatLng._dms = function(rad) 
{
  var d = Math.abs(rad * 180 / Math.PI);
  d += 1/7200;  // add ½ second for rounding
  var deg = Math.floor(d);
  var min = Math.floor((d-deg)*60);
  var sec = Math.floor((d-deg-min/60)*3600);
  // add leading zeros if required
  if (deg<100) deg = '0' + deg; if (deg<10) deg = '0' + deg;
  if (min<10) min = '0' + min;
  if (sec<10) sec = '0' + sec;
  return deg + '\u00B0' + min + '\u2032' + sec + '\u2033';
}


/*
 * override toPrecision method with one which displays trailing zeros in place
 *   of exponential notation
 *
 * (for Haversine, use 4 sf to reflect reasonable indication of accuracy)
 */
Number.prototype.toPrecision = function(fig) 
{
  var scale = Math.ceil(Math.log(this)*Math.LOG10E);
  var mult = Math.pow(10, fig-scale);
  return Math.round(this*mult)/mult;
}


/*
 * it's good form to include a toString method...
 */
LatLng.prototype.toString = function() 
{
  return this.latitude() + ', ' + this.longitude();
}


