#delimit;
capture clear;
clear matrix;
set mem 80m;
set matsize 800;
capture log close;
set more off;

*********This do file converts easting/northings for longitude/latitude;
*********The code is the stata version of the javascript code from: http://www.movable-type.co.uk/scripts/latlong-gridref.html;
*********The Ordnance Survey Spreasheet giving these calculations is located here http://www.ordnancesurvey.co.uk/oswebsite/gps/osnetfreeservices/furtherinfo/spreadsheet.html;
*********Elaine Kelly - July 2012;
*********Note: the eastings variable must be called "eastings", and the northings variable must be called "northings", replace "dataset.dta" with your dataset;

use "dataset.dta", clear;

  gen a = 6377563.396;
  gen b = 6356256.910;             
  *Airy 1830 major & minor semi-axes;
  gen f0 = 0.9996012717;                          
  *NatGrid scale factor on central meridian;
  gen lat0 = 49*_pi/180;
  
  gen lon0 = -2*_pi/180;   
  sum lon0;
  *NatGrid true origin;
  gen N0 = -100000;
  gen e0 = 400000;                     
  *northing & easting of true origin, metres;
  gen e2 = 1 - (b*b)/(a*a);                          
  *eccentricity squared;
  gen n = (a-b)/(a+b);
  gen n2 = n*n;
  gen n3 = n*n*n;

  gen lat = lat0;
  gen M=0;
  gen Ma = .;
  gen Mb = .;
  gen Mc = .;
  gen Md =.;
  
  *********Remember that your northings variable needs to be called northings, or replace "northings" with your variable in the cap confirm line below;
  
  local i = 1;
  while `i' <2{;
  
    replace lat = (northings-N0-M)/(a*f0) + lat;

    replace Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (lat-lat0);
    replace Mb = (3*n + 3*n*n + (21/8)*n3) * sin(lat-lat0) * cos(lat+lat0);
    replace Mc = ((15/8)*n2 + (15/8)*n3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0));
    replace Md = (35/24)*n3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0));
    replace M = b * f0 * (Ma - Mb + Mc - Md);               
	*meridional arc;
	*The sequence below instructs STATA to keep going until the distance is sufficiently small;
	cap confirm (northings-N0-M >= 0.00001);
	 if !_rc {;
	 };
	 else{;
	 local i = `i' +1;
	 };
};

sum lon0;

  gen cosLat = cos(lat);
  gen sinLat = sin(lat);
  gen nu = a*f0/((1-e2*sinLat*sinLat)^0.5);              
  * transverse radius of curvature;
  gen rho = a*f0*(1-e2)/((1-e2*sinLat*sinLat)^1.5);  
* meridional radius of curvature;
  gen eta2 = (nu/rho)-1;

  sum lon0;
  gen tanLat = tan(lat);
  gen tan2lat = tanLat*tanLat;
  gen tan4lat = tan2lat*tan2lat;
  gen tan6lat = tan4lat*tan2lat;
  gen secLat = 1/cosLat;
  gen double nu3 = nu*nu*nu;
  gen double nu5 = nu3*nu*nu;
  *note nu7 in the original calculation is not included as there are too many digits for STATA (replaced by nu5*nu*nu);
  format nu* %40.0g;  
  gen double VII = tanLat/(2*rho*nu);
  gen double VIII = tanLat/(24*rho*nu3)*(5+3*tan2lat+eta2-9*tan2lat*eta2);
  gen double IX = (tanLat/(720*rho*nu5)*(61+90*tan2lat+45*tan4lat));
  gen double X = secLat/nu;
  gen double XI = secLat/(6*nu3)*(nu/rho+2*tan2lat);
  gen double XII = secLat/(120*nu5)*(5+28*tan2lat+24*tan4lat);
  gen double XIIA = (secLat/(5040*nu5*nu*nu)*(61+662*tan2lat+1320*tan4lat+720*tan6lat));
  format V* X* %40.0g;
  
  *********Remember that your eastings variable needs to be called eastings, or replace "eastings" with your variable in the line below;
	
  gen DE = (eastings-e0);
  gen de2 = DE*DE;
  gen de3 = de2*DE;
  gen de4 = de2*de2;
  gen de5 = de3*DE*DE;
  gen de6 = de4*de2;
  gen de7 = de5*de2;
  replace lat = lat - (VII*de2) + VIII*de4 - (IX*de6);
  gen lon = lon0 + (X*DE) - (XI*de3) + (XII*de5) - (XIIA*de7);
  
  ****Now need to convert to radians;
gen long_degrees = ( 180/_pi )*lon;
label variable long_degrees "Longitude in Degrees";
gen lat_degrees = ( 180/_pi )*lat;
label variable lat_degrees "Latitide in Degrees";