function OGBLatLng(lat, lon)
{
	this.lat = lat;
	this.lon = lon;
}

function WGS84LatLng(lat, lon)
{
	this.lat = lat;
	this.lon = lon;
}

function OGBNorthEast(east, north)
{
	this.north = north;
	this.east = east;
}

function OGBRect(bottomLeft, topRight)
{
	this.bl = bottomLeft;
	this.tr = topRight;
}



//convert WGS84 Latitude and Longitude to Ordnance Survey 1936 Latitude Longitude

function WGS84ToOGB(WGlat, WGlon, height)
{
var deg2rad = Math.PI / 180;
var rad2deg = 180.0 / Math.PI;

//first off convert to radians
var radWGlat = WGlat * deg2rad;
var radWGlon = WGlon * deg2rad;
//these are the values for WGS84(GRS80) to OSGB36(Airy)
var a = 6378137; // WGS84_AXIS
var e = 0.00669438037928458; // WGS84_ECCENTRIC
var h = height; // height above datum (from $GPGGA sentence)
var a2 = 6377563.396; // OSGB_AXIS
var e2 = 0.0066705397616; // OSGB_ECCENTRIC 
var xp = -446.448;
var yp = 125.157;
var zp = -542.06;
var xr = -0.1502;
var yr = -0.247;
var zr = -0.8421;
var s = 20.4894;

// convert to cartesian; lat, lon are in radians
var sf = s * 0.000001;
var v = a / (Math.sqrt(1 - (e * (Math.sin(radWGlat) * Math.sin(radWGlat)))));
var x = (v + h) * Math.cos(radWGlat) * Math.cos(radWGlon);
var y = (v + h) * Math.cos(radWGlat) * Math.sin(radWGlon);
var z = ((1 - e) * v + h) * Math.sin(radWGlat);

// 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
var newLon = Math.atan(hy / hx);
var p = Math.sqrt((hx * hx) + (hy * hy));
var newLat = Math.atan(hz / (p * (1 - e2)));
v = a2 / (Math.sqrt(1 - e2 * (Math.sin(newLat) * Math.sin(newLat))));
var errvalue = 1.0;
var lat0 = 0;
while (errvalue > 0.001)
{
lat0 = Math.atan((hz + e2 * v * Math.sin(newLat)) / p);
errvalue = Math.abs(lat0 - newLat);
newLat = lat0;
}

//convert back to degrees
newLat = newLat * rad2deg;
newLon = newLon * rad2deg;

return new OGBLatLng(newLat, newLon);

}


//convert Ordnance Survey 1936 Latitude Longitude to WGS84 Latitude and Longitude

function OGBToWGS84(OGlat, OGlon, height)
{
var deg2rad = Math.PI / 180;
var rad2deg = 180.0 / Math.PI;

//first off convert to radians
var radOGlat = OGlat * deg2rad;
var radOGlon = OGlon * deg2rad;
//these are the values for WGS84(GRS80) to OSGB36(Airy)

var a2 = 6378137; // WGS84_AXIS
var e2 = 0.00669438037928458; // WGS84_ECCENTRIC

var h = height; // height above datum (from $GPGGA sentence)

var a = 6377563.396; // OSGB_AXIS
var e = 0.0066705397616; // OSGB_ECCENTRIC 

var xp = 446.448;
var yp = -125.157;
var zp = 542.06;
var xr = 0.1502;
var yr = 0.247;
var zr = 0.8421;

var s = -20.4894;

// convert to cartesian; lat, lon are in radians
var sf = s * 0.000001;
var v = a / (Math.sqrt(1 - (e * (Math.sin(radOGlat) * Math.sin(radOGlat)))));
var x = (v + h) * Math.cos(radOGlat) * Math.cos(radOGlon);
var y = (v + h) * Math.cos(radOGlat) * Math.sin(radOGlon);
var z = ((1 - e) * v + h) * Math.sin(radOGlat);

// 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
var newLon = Math.atan(hy / hx);
var p = Math.sqrt((hx * hx) + (hy * hy));
var newLat = Math.atan(hz / (p * (1 - e2)));
v = a2 / (Math.sqrt(1 - e2 * (Math.sin(newLat) * Math.sin(newLat))));
var errvalue = 1.0;
var lat0 = 0;
while (errvalue > 0.001)
{
lat0 = Math.atan((hz + e2 * v * Math.sin(newLat)) / p);
errvalue = Math.abs(lat0 - newLat);
newLat = lat0;
}

//convert back to degrees
newLat = newLat * rad2deg;
newLon = newLon * rad2deg;

return new WGS84LatLng(newLat, newLon);

}


//converts lat and lon (OSGB36) to OS northings and eastings
function LLtoNE(lat, lon)
{
var deg2rad = Math.PI / 180;
var rad2deg = 180.0 / Math.PI;

var phi = lat * deg2rad; // convert latitude to radians
var lam = lon * deg2rad; // convert longitude to radians
var a = 6377563.396; // OSGB semi-major axis
var b = 6356256.91; // OSGB semi-minor axis
var e0 = 400000; // easting of false origin
var n0 = -100000; // northing of false origin
var f0 = 0.9996012717; // OSGB scale factor on central meridian
var e2 = 0.0066705397616; // OSGB eccentricity squared
var lam0 = -0.034906585039886591; // OSGB false east
var phi0 = 0.85521133347722145; // OSGB false north
var af0 = a * f0;
var bf0 = b * f0;

// easting
var slat2 = Math.sin(phi) * Math.sin(phi);
var nu = af0 / (Math.sqrt(1 - (e2 * (slat2))));
var rho = (nu * (1 - e2)) / (1 - (e2 * slat2));
var eta2 = (nu / rho) - 1;
var p = lam - lam0;
var IV = nu * Math.cos(phi);
var clat3 = Math.pow(Math.cos(phi), 3);
var tlat2 = Math.tan(phi) * Math.tan(phi);
var V = (nu / 6) * clat3 * ((nu / rho) - tlat2);
var clat5 = Math.pow(Math.cos(phi), 5);
var tlat4 = Math.pow(Math.tan(phi), 4);
var VI = (nu / 120) * clat5 * ((5 - (18 * tlat2)) + tlat4 + (14 * eta2) - (58 * tlat2 * eta2));
var east = e0 + (p * IV) + (Math.pow(p, 3) * V) + (Math.pow(p, 5) * VI);

// northing
var n = (af0 - bf0) / (af0 + bf0);
var M = Marc(bf0, n, phi0, phi);
var I = M + (n0);
var II = (nu / 2) * Math.sin(phi) * Math.cos(phi);
var III = ((nu / 24) * Math.sin(phi) * Math.pow(Math.cos(phi), 3)) * (5 - Math.pow(Math.tan(phi), 2) + (9 * eta2));
var IIIA = ((nu / 720) * Math.sin(phi) * clat5) * (61 - (58 * tlat2) + tlat4);
var north = I + ((p * p) * II) + (Math.pow(p, 4) * III) + (Math.pow(p, 6) * IIIA);

// make whole number values
east = Math.round(east); // round to whole number of meters
north = Math.round(north); 

return new OGBNorthEast(east, north);
}


//convert northing and easting to letter and number grid system
function NE2NGR( east,  north)
{
east = Math.round(east);
north = Math.round(north);
var eX = east / 500000;
var nX = north / 500000;
var tmp = Math.floor(eX) - 5.0 * Math.floor(nX) + 17.0; 
nX = 5 * (nX - Math.floor(nX));
eX = 20 - 5.0 * Math.floor(nX) + Math.floor(5.0 * (eX - Math.floor(eX)));
if (eX > 7.5) eX = eX + 1; // I is not used
if (tmp > 7.5) tmp = tmp + 1; // I is not used

var eing = east - (Math.floor(east / 100000)*100000);
var ning = north - (Math.floor(north / 100000)*100000);
var estr = eing.toString();
var nstr = ning.toString();
while(estr.length < 5)
	estr = "0" + estr;
while(nstr.length < 5)
	nstr = "0" + nstr;

var ngr = String.fromCharCode(tmp + 65) + 
          String.fromCharCode(eX + 65) + 
          " " + estr + " " + nstr;

return ngr;
}

//helper
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)))));
}

function NEtoLL (east, north){
 
//metres in, degrees out

var K0 = 0.9996012717; // grid scale factor on central meridean
var OriginLat = 49.0;
var OriginLong = -2.0;    
var OriginX = 400000; // 400 kM
var OriginY = -100000; // 100 kM
var a = 6377563.396; // Airy Spheroid
var b = 6356256.910;

var 	e2;
var 	ex;
var 	n1;
var 	n2;
var 	n3;
var 	OriginNorthings;


// compute interim values
a = a * K0;
b = b * K0;
                
n1 = (a - b)/(a + b);
n2 = n1 * n1;
n3 = n2 * n1;   
                                                                                
lat = OriginLat * Math.PI / 180.0; // to radians                                                    
                

e2 = (a*a - b*b)/(a*a);  // first eccentricity
ex = (a*a - b*b)/(b*b);  // second eccentricity


OriginNorthings = b*lat + b*(n1*(1.0 + 5.0*n1*(1.0+n1)/4.0)*lat         
					- 3.0*n1*(1.0+n1*(1.0+7.0*n1/8.0))*Math.sin(lat)*Math.cos(lat)
					+ (15.0*n1*(n1+n2)/8.0)*Math.sin(2.0*lat)*Math.cos(2.0*lat)
					- (35.0*n3/24.0)*Math.sin(3.0*lat)*Math.cos(3.0*lat) );
	
var OriginLat = 49.0;
var OriginLong = -2.0;    
var OriginX = 400000; // 400 kM
var OriginY = -100000; // 100 kM


var lat;    // what we calculate
var lon;

var northing = north - OriginY; 
var easting = east - OriginX;

var nu, phid, phid2, t2, t, q2, c, s, nphid, dnphid; // temps
var nu2, nudivrho, invnurho, rho, eta2;


/* Evaluate M term: latitude of the northing on the centre meridian */  

northing += OriginNorthings;  
           
phid = northing / (b*(1.0 + n1 + 5.0*(n2+n3)/4.0)) - 1.0;
phid2 = phid + 1.0;                                

while (Math.abs(phid2 - phid) > 1E-6)
   {
   phid = phid2;
   nphid = b*phid + b*(n1*(1.0 + 5.0*n1*(1.0+n1)/4.0)*phid
           - 3.0*n1*(1.0+n1*(1.0+7.0*n1/8.0))*Math.sin(phid)*Math.cos(phid)
           + (15.0*n1*(n1+n2)/8.0)*Math.sin(2.0*phid)*Math.cos(2.0*phid)
           - (35.0*n3/24.0)*Math.sin(3.0*phid)*Math.cos(3.0*phid) );
           
   dnphid = b*((1.0+n1+5.0*(n2+n3)/4.0)-3.0*(n1+n2+7.0*n3/8.0)*Math.cos(2.0*phid)
            +(15.0*(n2+n3)/4.0)*Math.cos(4*phid)-(35.0*n3/8.0)*Math.cos(6.0*phid));
           
   phid2 = phid - (nphid - northing)/dnphid;
   }

c  = Math.cos(phid);
s  = Math.sin(phid);
t  = Math.tan(phid);
t2 = t*t;
q2 = easting*easting;


nu2 = (a*a)/(1.0 - e2*s*s);
nu = Math.sqrt(nu2);

nudivrho = a*a*c*c/(b*b) - c*c + 1.0;

eta2 = nudivrho - 1;

rho = nu / nudivrho;

invnurho = ((1.0-e2*s*s)*(1.0-e2*s*s))/(a*a*(1.0-e2));

lat = phid - t*q2*invnurho/2.0 + (q2*q2*(t/(24*rho*nu2*nu)*(5 + (3*t2) + eta2 -(9*t2*eta2))));

lon = (easting/(c*nu))
        - (easting*q2*((nudivrho+2.0*t2)/(6.0*nu2))/(c*nu))
		+ (q2*q2*easting*(5 + (28*t2) + (24*t2*t2))/(120*nu2*nu2*nu*c));
		
return new OGBLatLng(lat * 180.0 / Math.PI,(lon * 180.0 / Math.PI) + OriginLong);        

}

//return rect in OSGB N/E coords that encloses the given wgs84 rect
function enclosingOsgbRect(WGleft,WGbottom,WGtop,WGright){

	var blOGB = WGS84ToOGB(WGbottom,WGleft,0);
	var trOGB = WGS84ToOGB(WGtop,WGright,0);
	var brOGB = WGS84ToOGB(WGbottom,WGright,0);
	var tlOGB = WGS84ToOGB(WGtop,WGleft,0);
	
	var blEN = LLtoNE(blOGB.lat,blOGB.lon);
	var trEN = LLtoNE(trOGB.lat,trOGB.lon);
	var brEN = LLtoNE(brOGB.lat,brOGB.lon);
	var tlEN = LLtoNE(tlOGB.lat,tlOGB.lon);
	
	var e = Math.min(blEN.east,tlEN.east); 
	var w = Math.max(brEN.east,trEN.east); 
	var s = Math.min(blEN.north,brEN.north); 
	var n = Math.max(trEN.north,tlEN.north); 
	
	return new OGBRect(new OGBNorthEast(e,s),new OGBNorthEast(w,n));
}

function OsgbConvergence(lat,lon)
{

var K0 = 0.9996012717; // grid scale factor on central meridean
var a = 6377563.396; // Airy Spheroid
var b = 6356256.910;
var oDeg = -2.0;//longitude origin
a = a * K0;
b = b * K0;

/* Compute convergence, OS Geodetic Information Paper No 1 */
var phi = lat * Math.PI / 180.0;
var c = Math.cos(phi);
var s = Math.sin(phi);
var t = Math.tan(phi);
var t2 = t*t;
var p = (lon - oDeg) * Math.PI / 180.0;
var eta2 = a*a*c*c/(b*b) - c*c;

return (
( p*s ) +
( p*p*p * s*c*c * (1.0 + (3.0*eta2) + (2.0*eta2*eta2)) / 3.0 ) +
( p*p*p*p*p * s*c*c*c*c * (2.0 - t2) / 15.0 )
)
* 180.0 / Math.PI;

}