/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/* Vincenty Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2010             */
/*                                                                                                */
/* 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                                             */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

/**
* Calculates geodetic distance between two points specified by latitude/longitude using 
* Vincenty inverse formula for ellipsoids
*
* @param   {Number} lat1, lon1: first point in decimal degrees
* @param   {Number} lat2, lon2: second point in decimal degrees
* @returns (Number} distance in metres between points
*/
function distVincenty(lat1, lon1, lat2, lon2) {
    var a = 6378137, b = 6356752.3142, f = 1 / 298.257223563;  // WGS-84 ellipsoid params
    var L = (lon2 - lon1).toRad();
    var U1 = Math.atan((1 - f) * Math.tan(lat1.toRad()));
    var U2 = Math.atan((1 - f) * Math.tan(lat2.toRad()));
    var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
    var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);

    var lambda = L, lambdaP, iterLimit = 100;
    do {
        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 sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
        var cosSqAlpha = 1 - sinAlpha * sinAlpha;
        var cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
        if (isNaN(cos2SigmaM)) cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (§6)
        var C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
        lambdaP = lambda;
        lambda = L + (1 - C) * f * sinAlpha *
      (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
    } while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0);

    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);

    s = s.toFixed(3); // round to 1mm precision
    return s;
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */


// ellipse parameters
var e = { WGS84: { a: 6378137, b: 6356752.3142, f: 1 / 298.257223563 },
    Airy1830: { a: 6377563.396, b: 6356256.910, f: 1 / 299.3249646}
};

// helmert transform parameters
var h = { WGS84toOSGB36: { tx: -446.448, ty: 125.157, tz: -542.060,   // m
    rx: -0.1502, ry: -0.2470, rz: -0.8421,  // sec
    s: 20.4894
},                               // ppm
    OSGB36toWGS84: { tx: 446.448, ty: -125.157, tz: 542.060,
        rx: 0.1502, ry: 0.2470, rz: 0.8421,
        s: -20.4894}
    };


    function convertOSGB36toWGS84(p1) {
        var p2 = convert(p1, e.Airy1830, h.OSGB36toWGS84, e.WGS84);
        return p2;
    }


    function convertWGS84toOSGB36(p1) {
        var p2 = convert(p1, e.WGS84, h.WGS84toOSGB36, e.Airy1830);
        return p2;
    }


    function convert(p1, e1, t, e2) {
        // -- convert polar to cartesian coordinates (using ellipse 1)

        p1.lat = p1.lat.toRad(); p1.lon = p1.lon.toRad();

        var a = e1.a, b = e1.b;

        var sinPhi = Math.sin(p1.lat), cosPhi = Math.cos(p1.lat);
        var sinLambda = Math.sin(p1.lon), cosLambda = Math.cos(p1.lon);
        var H = p1.height;

        var eSq = (a * a - b * b) / (a * a);
        var nu = a / Math.sqrt(1 - eSq * sinPhi * sinPhi);

        var x1 = (nu + H) * cosPhi * cosLambda;
        var y1 = (nu + H) * cosPhi * sinLambda;
        var z1 = ((1 - eSq) * nu + H) * sinPhi;


        // -- apply helmert transform using appropriate params

        var tx = t.tx, ty = t.ty, tz = t.tz;
        var rx = t.rx / 3600 * Math.PI / 180;  // normalise seconds to radians
        var ry = t.ry / 3600 * Math.PI / 180;
        var rz = t.rz / 3600 * Math.PI / 180;
        var s1 = t.s / 1e6 + 1;              // normalise ppm to (s+1)

        // apply transform
        var x2 = tx + x1 * s1 - y1 * rz + z1 * ry;
        var y2 = ty + x1 * rz + y1 * s1 - z1 * rx;
        var z2 = tz - x1 * ry + y1 * rx + z1 * s1;


        // -- convert cartesian to polar coordinates (using ellipse 2)

        a = e2.a, b = e2.b;
        var precision = 4 / a;  // results accurate to around 4 metres

        eSq = (a * a - b * b) / (a * a);
        var p = Math.sqrt(x2 * x2 + y2 * y2);
        var phi = Math.atan2(z2, p * (1 - eSq)), phiP = 2 * Math.PI;
        while (Math.abs(phi - phiP) > precision) {
            nu = a / Math.sqrt(1 - eSq * Math.sin(phi) * Math.sin(phi));
            phiP = phi;
            phi = Math.atan2(z2 + eSq * nu * Math.sin(phi), p);
        }
        var lambda = Math.atan2(y2, x2);
        H = p / Math.cos(phi) - nu;

        return new LatLon(phi.toDeg(), lambda.toDeg(), H);
    }


    // ---- the following are duplicated from LatLong.html ---- //


    /*
    * construct a LatLon object: arguments in numeric degrees & metres
    *
    * note all LatLong methods expect & return numeric degrees (for lat/long & for bearings)
    */
    function LatLon(lat, lon, height) {
        if (arguments.length < 3) height = 0;
        this.lat = lat;
        this.lon = lon;
        this.height = height;
    }

    /*
    * represent point {lat, lon} in standard representation
    */
    LatLon.prototype.toString = function() {
        return this.lat.toLat() + ', ' + this.lon.toLon();
    }


    // extend String object with method for parsing degrees or lat/long values to numeric degrees
    //
    // this is very flexible on formats, allowing signed decimal degrees, 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).

    String.prototype.parseDeg = function() {
        if (!isNaN(this)) return Number(this);                 // signed decimal degrees without NSEW

        var degLL = this.replace(/^-/, '').replace(/[NSEW]/i, '');  // strip off any sign or compass dir'n
        var dms = degLL.split(/[^0-9.]+/);                     // split out separate d/m/s
        for (var i in dms) if (dms[i] == '') dms.splice(i, 1);    // remove empty elements (see note below)
        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:                                              // decimal or non-separated dddmmss
                if (/[NS]/i.test(this)) degLL = '0' + degLL;       // - normalise N/S to 3-digit degrees
                var deg = dms[0].slice(0, 3) / 1 + dms[0].slice(3, 5) / 60 + dms[0].slice(5) / 3600; break;
            default: return NaN;
        }
        if (/^-/.test(this) || /[WS]/i.test(this)) deg = -deg; // take '-', west and south as -ve
        return deg;
    }
    // note: whitespace at start/end will split() into empty elements (except in IE)


    // extend Number object with methods for converting degrees/radians

    Number.prototype.toRad = function() {  // convert degrees to radians
        return this * Math.PI / 180;
    }

    Number.prototype.toDeg = function() {  // convert radians to degrees (signed)
        return this * 180 / Math.PI;
    }

    // extend Number object with methods for presenting bearings & lat/longs

    Number.prototype.toDMS = function(dp) {  // convert numeric degrees to deg/min/sec
        if (arguments.length < 1) dp = 0;      // if no decimal places argument, round to int seconds
        var d = Math.abs(this);  // (unsigned result ready for appending compass dir'n)
        var deg = Math.floor(d);
        var min = Math.floor((d - deg) * 60);
        var sec = ((d - deg - min / 60) * 3600).toFixed(dp);
        // fix any nonsensical rounding-up
        if (sec == 60) { sec = (0).toFixed(dp); min++; }
        if (min == 60) { min = 0; deg++; }
        if (deg == 360) deg = 0;
        // 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';
    }

    Number.prototype.toLat = function(dp) {  // convert numeric degrees to deg/min/sec latitude
        return this.toDMS(dp).slice(1) + (this < 0 ? 'S' : 'N');  // knock off initial '0' for lat!
    }

    Number.prototype.toLon = function(dp) {  // convert numeric degrees to deg/min/sec longitude
        return this.toDMS(dp) + (this > 0 ? 'E' : 'W');
    }