/* 
 * Convert degrees to radians
 */
rad = function(x) {return x*Math.PI/180;}

/*
 * 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
 */
distHaversine = function(p1, p2) {
  var R = 6371; // earth's mean radius in km
  var dLat  = rad(p2.y - p1.y);
  var dLong = rad(p2.x - p1.x);

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

  return d.toFixed(3);
}


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

/*
 * Calculate geodesic distance (in km) between two points specified by latitude/longitude with Vincenty formula
 *
 * from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the 
 *       Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975
 *       http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
 */
distVincenty = function(p1, p2) {
  var a = 6378137, b = 6356752.3142,  f = 1/298.257223563;
  var L = rad(p2.x - p1.x);
  var U1 = Math.atan((1-f) * Math.tan(rad(p1.y)));
  var U2 = Math.atan((1-f) * Math.tan(rad(p2.y)));
  var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
  var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
  var lambda = L, lambdaP = 2*Math.PI;
  var iterLimit = 20;
  while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0) {
    var sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
    var sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + 
      (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
    if (sinSigma==0) return 0;  // co-incident points
    var cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
    var sigma = Math.atan2(sinSigma, cosSigma);
    var alpha = Math.asin(cosU1 * cosU2 * sinLambda / sinSigma);
    var cosSqAlpha = Math.cos(alpha) * Math.cos(alpha);
    var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
    var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
    lambdaP = lambda;
    lambda = L + (1-C) * f * Math.sin(alpha) *
      (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
  }
  if (iterLimit==0) return NaN  // formula failed to converge
  var uSq = cosSqAlpha * (a*a - b*b) / (b*b);
  var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
  var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
    B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
  var s = b*A*(sigma-deltaSigma);
  var d = s.toFixed(3)/1000; // round to 1mm precision
  return d;
}


bearing = function(p1,p2) {
  var p3 = new GPoint(p2.x,p1.y);
  var dy = distHaversine(p3,p2); 
  if (p1.y > p2.y) {dy = -dy;}
  var dx = distHaversine(p1,p3);
  if (p1.x > p2.x) {dx = -dx;}
  var b = Math.atan2(dx,dy) * 180/Math.PI;
  if (b<0) {b=b+360;}
      b = b.toFixed(1);
  return b;
}