// OS to lat/lon conversions
// Bill Chadwick, http://www.bdcc.co.uk/Gmaps/BdccGmapBits.htm
// modified for Trailwise with LandRanger funcs, speedups. version 0.2 Ross K 9/2007

//global constants
var deg2rad=Math.PI/180;
var rad2deg=180.0/Math.PI;

function WGS84toNGR(lat,lon) // return AA NNNN NNNN style (10m LR grid)
{
	var tempOGB=WGS84ToOGB(lat,lon,0);
	var tempEN=LLtoNE(tempOGB.lat,tempOGB.lon);
	var tempNGR=NE2NGR(tempEN.east,tempEN.north);
	return tempNGR.substr(0,7)+tempNGR.substr(8,5);
}

function NGRtoWGS84(landranger) // convert LR grid to WGS84, false on parse fail
{
	var tempEN=NGRtoNE(landranger);
	if (tempEN.east==0 && tempEN.north==0) {return false}
	tempOGB=NEtoLL(tempEN.east,tempEN.north);
	return OGBToWGS84(tempOGB.lat,tempOGB.lon,0);
}


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 Lat Lon to OS 1936 Lat Lon
function WGS84ToOGB(WGlat,WGlon, height)
{
	//convert to radians
var radWGlat=WGlat*deg2rad;
var radWGlon=WGlon*deg2rad;
	//values for WGS84(GRS80) to OSGB36(Airy) datum
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;
var sf=0.0000204894; //s*0.000001

	//convert to cartesian; lat,lon are in radians
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 OS 1936 Lat Lon to WGS84 Lat Lon
function OGBToWGS84 (OGlat,OGlon,height)
{
	//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;
var sf=-0.0000204894; //s*0.000001

	//convert to cartesian; lat,lon are in radians
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 e0=400000; //easting of false origin
var n0=-100000; //northing of false origin
var e2=0.0066705397616; //OSGB eccentricity squared
var lam0=-0.034906585039886591; //OSGB false east
var phi0=0.85521133347722145; //OSGB false north
var af0=6375020.48099; //(a) OSGB semi-major axis * (f0) OSGB scale factor on central meridian
var bf0=6353722.49049; //(b) semi-minor axis *f0

var phi=lat*deg2rad; //convert latitude to radians
var lam=lon*deg2rad; //convert longitude to radians

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

east=Math.round(east); //round to whole number of meters
north=Math.round(north); 
return new OGBNorthEast(east,north);
}
////helper
function Marc(bf0,n,phi0,phi)
{
var n2=n*n;
var n3=n*n*n;
return bf0*(((1+n+((5/4)*(n2))+((5/4)*(n3)))*(phi-phi0))
-(((3*n)+(3*(n2))+((21/8)*(n3)))*(Math.sin(phi-phi0))*(Math.cos(phi+phi0)))
+((((15/8)*(n2))+((15/8)*(n3)))*(Math.sin(2*(phi-phi0)))*(Math.cos(2*(phi+phi0))))
-(((35/24)*(n3))*(Math.sin(3*(phi-phi0)))*(Math.cos(3*(phi+phi0)))));
}


////convert northing and easting to letter and number grid system at 1m
function NE2NGR(east,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}

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} //pad leading zeros
while(nstr.length<5) {nstr="0"+nstr}
var ngr=String.fromCharCode(tmp+65) + String.fromCharCode(eX+65) + " " + estr + " " + nstr;
return ngr;
}


//// EN metres in, OSGB degrees out
function NEtoLL (east,north){
//var lat=0.8552113334772214; //OriginLat 49.0 deg,  in radians
var OriginLong=-2.0;    
var a=6375020.48099; //Airy Spheroid 6377563.396 * K0 0.9996012717 grid scale factor
var b=6353722.49049; //6356256.910 * 0.9996012717

	//compute interim values
var a2 = 40640886133041.97; //a squared
var n1=0.0016732202502415632; //(a-b)/(a+b)
var n2=0.000002799666005010; //n1 squared
var n3=4.684457854848326E-9; //n1 cubed
                                                                                
var e2=0.00667053976126558; //(a*a-b*b)/(a*a); //first eccentricity
var ex=0.00671533466817988; //(a*a-b*b)/(b*b); //second eccentricity

var OriginNorthings=5427063.815694233;

var easting=east-400000; //east-OriginX;

	//Evaluate M term: latitude of the northing on the centre meridian
var northing=north+100000+OriginNorthings; //north-OriginY+ON
var phid=northing/(b*(1.0+n1+5.0*(n2+n3)/4.0))-1.0;
var phid2=phid+1.0;                                
var nphid,dnphid; //temps
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;
}

var c=Math.cos(phid);
var s=Math.sin(phid);
var t=Math.tan(phid);
var t2=t*t;
var e2s2=e2*s*s;
var q2=easting*easting;
var nu2=(a2)/(1.0-e2s2);
var nu=Math.sqrt(nu2);
var nudivrho=a2*c*c/(b*b)-c*c+1.0;
var eta2=nudivrho-1;
var rho=nu/nudivrho;
var invnurho=((1.0-e2s2)*(1.0-e2s2))/(a2*(1.0-e2));

var lat=phid-t*q2*invnurho/2.0+(q2*q2*(t/(24*rho*nu2*nu)*(5+(3*t2)+eta2-(9*t2*eta2))));
var 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*rad2deg,(lon*rad2deg)+OriginLong);        
}


//// global array for NGRtoNE
var LRprefixes = new Array (
	new Array("SV","SW","SX","SY","SZ","TV","TW"),
	new Array("SQ","SR","SS","ST","SU","TQ","TR"),
	new Array("SL","SM","SN","SO","SP","TL","TM"),
	new Array("SF","SG","SH","SJ","SK","TF","TG"),
	new Array("SA","SB","SC","SD","SE","TA","TB"),
	new Array("NV","NW","NX","NY","NZ","OV","OW"),
	new Array("NQ","NR","NS","NT","NU","OQ","OR"),
	new Array("NL","NM","NN","NO","NP","OL","OM"),
	new Array("NF","NG","NH","NJ","NK","OF","OG"),
	new Array("NA","NB","NC","ND","NE","OA","OB"),
	new Array("HV","HW","HX","HY","HZ","JV","JW"),
	new Array("HQ","HR","HS","HT","HU","JQ","JR"),
	new Array("HL","HM","HN","HO","HP","JL","JM"));

//// parse Landranger string to eastings,northings; zeros if invalid
function NGRtoNE( landranger ) {
var tempLR=landranger.replace(/\s+/g,'');
tempLR=tempLR.toUpperCase();
var ok=false;
var northings=0;
var eastings=0;
for (var precision=5; precision>=1; precision--) {
	var pattern=new RegExp("^([A-Z]{2})\\s*(\\d{"+precision+"})\\s*(\\d{"+precision+"})$","i")
	var gridRef=tempLR.match(pattern);
	if (gridRef) {
		var gridSheet=gridRef[1];
		var gridEast=0;
		var gridNorth=0;
			//5x1 4x10 3x100 2x1000 1x10000 
		if (precision>0) {
			var mult=Math.pow(10,5-precision);
			gridEast=parseInt(gridRef[2],10)*mult;
			gridNorth=parseInt(gridRef[3],10)*mult;
		}
		var x,y;
		search: for(y=0; y<LRprefixes.length; y++) {
			for(x=0; x<LRprefixes[y].length; x++) 
				if (LRprefixes[y][x]==gridSheet) {
					eastings=(x*100000)+gridEast;
					northings=(y*100000)+gridNorth;
					ok=true;
					break search;
				}
		}
	}
}
return new OGBNorthEast(eastings,northings);
}