/* Coordinate conversion functions */
 
var rad2deg = 57.29577951;
var deg2rad = 0.017453292;
 
// Constants for WGS84
var WGS84SemiMajor = 6378137.0;
var WGS84Eccentricity = 0.00669438037928458;
 
// Airy Spheroid constants for BNG
var OSGBSemiMajor = 6377563.396;            // OSGB Semi-major axis
var OSGBSemiMinor = 6356256.91;             // OSGB semi-minor axis
var OSGBEastingOrigin = 400000;             // OSGB easting of false origin
var OSGBNorthingOrigin = -100000;           // OSGB northing of false origin
var OSGBScaleFactor = 0.9996012717;         // OSGB scale factor on central meridian
var OSGBEccentricity = 0.0066705397616;     // OSGB eccentricity squared
var OSGBFalseEast = -0.034906585039886591;  // OSGB false east
var OSGBFalseNorth = 0.85521133347722145;   // OSGB false north
// Modified Airy Speriod constants for ING
var OSISemiMajor = 6377340.189;                         // OSI semi-major
var OSISemiMinor = 6356034.447;                         // OSI semi-minor
var OSIEastingOrigin = 200000;                          // OSI easting of false origin
var OSINorthingOrigin = 250000;                         // OSI northing of false origin
var OSIScaleFactor = 1.000035;                          // OSI scale factor on central meridian
var OSIEccentricity = 0.00667054015;                    // OSI eccentricity squared
var OSIFalseEast = -0.13962634015954636615389526147909; // OSI false east
var OSIFalseNorth = 0.93375114981696632365417456114141; // OSI false north
 
// Function to convert BNG coordinates into WGS84
function ConvertBNGtoWGS84(easting, northing)
{
    // Convert to lat/long using the Airy Spheroid
    var geoLL = NEtoLL(Number(easting), Number(northing),
                       OSGBSemiMajor, OSGBSemiMinor,
                       OSGBEastingOrigin, OSGBNorthingOrigin,
                       OSGBScaleFactor, OSGBEccentricity,
                       OSGBFalseEast, OSGBFalseNorth);
                       
    var h = 10; // Dummy height
    
    // Apply the dataum shift (using 3 parameter Helmert transformation to match the ESRI projection parameters)
    var geo = transform(geoLL.latitude, geoLL.longitude,
                        OSGBSemiMajor, OSGBEccentricity, h, 
                        WGS84SemiMajor, WGS84Eccentricity,
                        375.0, -111.0, 431.0, 0.0, 0.0, 0.0, 0.0);
                        
    // Covert from radians to degrees and return
    geoLL.latitude = geo.latitude * rad2deg;
    geoLL.longitude = geo.longitude * rad2deg;
    
    return geoLL;                
}
 
// Function to convert ING coordinates into WGS84
function ConvertINGtoWGS84(easting, northing)
{
    // Convert to lat/long using the Modified Airy Spheroid
    var geoLL = NEtoLL(Number(easting), Number(northing),
                       OSISemiMajor, OSISemiMinor,
                       OSIEastingOrigin, OSINorthingOrigin,
                       OSIScaleFactor, OSIEccentricity,
                       OSIFalseEast, OSIFalseNorth);
                       
    var h = 10; // Dummy height
    
    // Apply the dataum shift (using 3 parameter Helmert transformation to match the ESRI projection parameters)
    var geo = transform(geoLL.latitude, geoLL.longitude,
                        OSISemiMajor, OSIEccentricity, h, 
                        WGS84SemiMajor, WGS84Eccentricity,
                        506.0, -122.0, 611.0, 0.0, 0.0, 0.0, 0.0);
                        
    // Covert from radians to degrees and return
    geoLL.latitude = geo.latitude * rad2deg;
    geoLL.longitude = geo.longitude * rad2deg;
    
    return geoLL;                
}
 
// Helper function to convert Northings and Eastings to Lat/Long
// Input easting and northing are in metres
// Output latitude and longitude are in radians
function NEtoLL(easting, northing, a, b, e0, n0, f0, e2, lam0, phi0)
{
    var af0 = a * f0;
    var bf0 = b * f0;
    var n = (af0 - bf0) / (af0 + bf0);
    var Et = easting - e0;
    var phid = InitialLat(northing, n0, af0, phi0, n, bf0);
    var nu = af0 / (Math.sqrt(1 - (e2 * (Math.sin(phid) * Math.sin(phid)))));
    var rho = (nu * (1 - e2)) / (1 - (e2 * (Math.sin(phid)) * (Math.sin(phid))));
    var eta2 = (nu / rho) - 1;
    var tlat2 = Math.tan(phid) * Math.tan(phid);
    var tlat4 = Math.pow(Math.tan(phid), 4);
    var tlat6 = Math.pow(Math.tan(phid), 6);
    var clatm1 = Math.pow(Math.cos(phid), -1);
    var VII = Math.tan(phid) / (2 * rho * nu);
    var VIII = (Math.tan(phid) / (24 * rho * (nu * nu * nu))) * (5 + (3 * tlat2) + eta2 - (9 * eta2 * tlat2));
    var IX = ((Math.tan(phid)) / (720 * rho * Math.pow(nu, 5))) * (61 + (90 * tlat2) + (45 * Math.pow(Math.tan(phid), 4) ));
    var phip = (phid - ((Et * Et) * VII) + (Math.pow(Et, 4) * VIII) - (Math.pow(Et, 6) * IX)); 
    var X = Math.pow(Math.cos(phid), -1) / nu;
    var XI = (clatm1 / (6 * (nu * nu * nu))) * ((nu / rho) + (2 * (tlat2)));
    var XII = (clatm1 / (120 * Math.pow(nu, 5))) * (5 + (28 * tlat2) + (24 * tlat4));
    var XIIA = clatm1 / (5040 * Math.pow(nu, 7)) * (61 + (662 * tlat2) + (1320 * tlat4) + (720 * tlat6));
    var lambdap = (lam0 + (Et * X) - ((Et * Et * Et) * XI) + (Math.pow(Et, 5) * XII) - (Math.pow(Et, 7) * XIIA));
    
    var geo = { latitude: phip, longitude: lambdap };
    return geo;
}
 
// Function to apply the datum transformation using 7 parameter Helmert transformation
function transform(lat, lon, a, e, h, a2, e2, xp, yp, zp, xr, yr, zr, s)
{
    // Convert to cartesian; lat, lon are radians
    var sf = s * 0.000001;
    var v = a / (Math.sqrt(1 - (e *(Math.sin(lat) * Math.sin(lat)))));
    var x = (v + h) * Math.cos(lat) * Math.cos(lon);
    var y = (v + h) * Math.cos(lat) * Math.sin(lon);
    var z = ((1 - e) * v + h) * Math.sin(lat);
    
    // Transform cartesian
    var xrot = (xr / 3600) * deg2rad;
    var yrot = (yr / 3600) * deg2rad;
    var zrot = (zr / 3600) * deg2rad;
    var hx = x + (x * sf) - (y * zrot) + (z * yrot) + xp;
    var hy = (x * zrot) + y + (y * sf) - (z * xrot) + yp;
    var hz = (-1 * x * yrot) + (y * xrot) + z + (z * sf) + zp;
    
    // Convert back to lat, lon
    lon = Math.atan(hy / hx);
    var p = Math.sqrt((hx * hx) + (hy * hy));
    lat = Math.atan(hz / (p * (1 - e2)));
    v = a2 / (Math.sqrt(1 - e2 * (Math.sin(lat) * Math.sin(lat))));
    var errvalue = 1.0;
    var lat0 = 0;
    while (errvalue > 0.001)
    {
      lat0 = Math.atan((hz + e2 * v * Math.sin(lat)) / p);  
      errvalue = Math.abs(lat0 - lat);
      lat = lat0;
    }
    h = p / Math.cos(lat) - v;
    var geo = { latitude: lat, longitude: lon, height: h };  // object to hold lat and lon
    return geo;
}    
 
// Helper function to compute the initial value for latitide
function InitialLat(north, n0, af0, phi0, n, bf0)
{
  var phi1 = ((north - n0) / af0) + phi0;
  var M = Marc(bf0, n, phi0, phi1);
  var phi2 = ((north - n0 - M) / af0) + phi1;
  var ind = 0;
  while ((Math.abs(north - n0 - M) > 0.00001) && (ind < 20))  // max 20 iterations in case of error
  {  
            ind = ind + 1;
            phi2 = ((north - n0 - M) / af0) + phi1;
    M = Marc(bf0, n, phi0, phi2);
    phi1 = phi2;
  }
  
  return phi2;  
}
 
// Helper function to compute the meridional arc
function Marc(bf0, n, phi0, phi)
{
  return bf0 * (((1 + n + ((5 / 4) * (n * n)) + ((5 / 4) * (n * n * n))) * (phi - phi0)) -
                (((3 * n) + (3 * (n * n)) + ((21 / 8) * (n * n * n))) * (Math.sin(phi - phi0)) * (Math.cos(phi + phi0))) +
                ((((15 / 8) * (n * n)) + ((15 / 8) * (n * n * n))) * (Math.sin(2 * (phi - phi0))) * (Math.cos(2 * (phi + phi0)))) -
                (((35 / 24) * (n * n * n)) * (Math.sin(3 * (phi - phi0))) * (Math.cos(3 * (phi + phi0)))));
}