I am having problems with the following code converting to lat long.
The following code is always off by 50 metres and I can't see anything wrong with it:
For example if you enter easting and northing 517835.99 and 428747.39 it should translate to latlong 53.7418 and -0.21484 but instead it gives 53.74159 and -0.21312.
// Define the constants for the BNG coordinate system
var a = 6377563.396; // Semi-major axis of the Airy 1830 ellipsoid
var b = 6356256.909; // Semi-minor axis of the Airy 1830 ellipsoid
var e0 = 400000; // False easting
var n0 = -100000; // False northing
var F0 = 0.9996012717; // Central meridian scale factor
var phi0 = 49 * (Math.PI / 180); // Central meridian (latitude of true origin) in radians
var lambda0 = -2 * (Math.PI / 180); // Central meridian (longitude of true origin) in radians
// Define the easting and northing coordinates
var easting = 517835.99; // Replace with your easting value
var northing = 428747.39; // Replace with your northing value
// Calculate the transformation variables
var e2 = (Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(a, 2); // Eccentricity squared
var n = (a - b) / (a + b);
var n2 = Math.pow(n, 2);
var n3 = Math.pow(n, 3);
var phi = phi0;
var M = 0;
do {
phi = (northing - n0 - M) / (a * F0) + phi;
var Ma = (1 + n + (5 / 4) * n2 + (5 / 4) * n3) * (phi - phi0);
var Mb = (3 * n + 3 * n * n + (21 / 8) * n3) * Math.sin(phi - phi0) * Math.cos(phi + phi0);
var Mc = ((15 / 8) * n2 + (15 / 8) * n3) * Math.sin(2 * (phi - phi0)) * Math.cos(2 * (phi + phi0));
var Md = (35 / 24) * n3 * Math.sin(3 * (phi - phi0)) * Math.cos(3 * (phi + phi0));
M = b * F0 * (Ma - Mb + Mc - Md);
} while (northing - n0 - M >= 0.00001);
var cosPhi = Math.cos(phi);
var sinPhi = Math.sin(phi);
var tanPhi = Math.tan(phi);
var secPhi = 1 / cosPhi;
var nu = a * F0 / Math.sqrt(1 - e2 * Math.pow(sinPhi, 2));
var rho = a * F0 * (1 - e2) / Math.pow(1 - e2 * Math.pow(sinPhi, 2), 1.5);
var eta2 = nu / rho - 1;
var VII = tanPhi / (2 * rho * nu);
var VIII = tanPhi / (24 * rho * Math.pow(nu, 3)) * (5 + 3 * Math.pow(tanPhi, 2) + eta2 - 9 * Math.pow(tanPhi, 2) * eta2);
var IX = tanPhi / (720 * rho * Math.pow(nu, 5)) * (61 + 90 * Math.pow(tanPhi, 2) + 45 * Math.pow(tanPhi, 4));
var X = secPhi / nu;
var XI = secPhi / (6 * Math.pow(nu, 3)) * (nu / rho + 2 * Math.pow(tanPhi, 2));
var XII = secPhi / (120 * Math.pow(nu, 5)) * (5 + 28 * Math.pow(tanPhi, 2) + 24 * Math.pow(tanPhi, 4));
var XIIA = secPhi / (5040 * Math.pow(nu, 7)) * (61 + 662 * Math.pow(tanPhi, 2) + 1320 * Math.pow(tanPhi, 4) + 720 * Math.pow(tanPhi, 6));
var dE = easting - e0;
var lat = phi - VII * Math.pow(dE, 2) + VIII * Math.pow(dE, 4) - IX * Math.pow(dE, 6);
var lon = lambda0 + X * dE - XI * Math.pow(dE, 3) + XII * Math.pow(dE, 5) - XIIA * Math.pow(dE, 7);
// Convert latitude and longitude to degrees
var latitude = lat * (180 / Math.PI);
var longitude = lon * (180 / Math.PI);
// The single point is considered as the centroid
var centroid = [latitude, longitude];
// Output the results
console.log("Latitude: " + latitude);
console.log("Longitude: " + longitude);
// Output the centroid coordinates
console.log("Centroid: " + centroid); Solution works for kingston upon hull:
I have taken the above code and because I needed it in raw maths form or jsmaths. I have simply pushed the equated values. I have done a lot of research to no avail and this gives me required values up to 1.6 metres from actual values.
public class Main {
public static void main(String[] args) {
// Define the constants for the BNG coordinate system
double a = 6377563.396; // Semi-major axis of the Airy 1830 ellipsoid
double b = 6356256.909; // Semi-minor axis of the Airy 1830 ellipsoid
double e0 = 400000; // False easting
double n0 = -100000; // False northing
double F0 = 0.9996012717; // Central meridian scale factor
double phi0 = 49 * (Math.PI / 180); // Central meridian (latitude of true origin) in radians
double lambda0 = -2 * (Math.PI / 180); // Central meridian (longitude of true origin) in radians
// Define the easting and northing coordinates
double easting = 510124.08; // Replace with your easting value
double northing = 436343.05; // Replace with your northing value
// Calculate the transformation variables
double e2 = (Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(a, 2); // Eccentricity squared
double n = (a - b) / (a + b);
double n2 = Math.pow(n, 2);
double n3 = Math.pow(n, 3);
double phi = phi0;
double M = 0;
do {
phi = (northing - n0 - M) / (a * F0) + phi;
double Ma = (1 + n + (5 / 4.0) * n2 + (5 / 4.0) * n3) * (phi - phi0);
double Mb = (3 * n + 3 * n * n + (21 / 8.0) * n3) * Math.sin(phi - phi0) * Math.cos(phi + phi0);
double Mc = ((15 / 8.0) * n2 + (15 / 8.0) * n3) * Math.sin(2 * (phi - phi0)) * Math.cos(2 * (phi + phi0));
double Md = (35 / 24.0) * n3 * Math.sin(3 * (phi - phi0)) * Math.cos(3 * (phi + phi0));
M = b * F0 * (Ma - Mb + Mc - Md);
} while (Math.abs(northing - n0 - M) >= 0.00001); // Use absolute value for convergence check
double cosPhi = Math.cos(phi);
double sinPhi = Math.sin(phi);
double tanPhi = Math.tan(phi);
double secPhi = 1 / cosPhi;
double nu = a * F0 / Math.sqrt(1 - e2 * Math.pow(sinPhi, 2));
double rho = a * F0 * (1 - e2) / Math.pow(1 - e2 * Math.pow(sinPhi, 2), 1.5);
double eta2 = nu / rho - 1;
double VII = tanPhi / (2 * rho * nu);
double VIII = tanPhi / (24 * rho * Math.pow(nu, 3)) * (5 + 3 * Math.pow(tanPhi, 2) + eta2 - 9 * Math.pow(tanPhi, 2) * eta2);
double IX = tanPhi / (720 * rho * Math.pow(nu, 5)) * (61 + 90 * Math.pow(tanPhi, 2) + 45 * Math.pow(tanPhi, 4));
double X = secPhi / nu;
double XI = secPhi / (6 * Math.pow(nu, 3)) * (nu / rho + 2 * Math.pow(tanPhi, 2));
double XII = secPhi / (120 * Math.pow(nu, 5)) * (5 + 28 * Math.pow(tanPhi, 2) + 24 * Math.pow(tanPhi, 4));
double XIIA = secPhi / (5040 * Math.pow(nu, 7)) * (61 + 662 * Math.pow(tanPhi, 2) + 1320 * Math.pow(tanPhi, 4) + 720 * Math.pow(tanPhi, 6));
double dE = easting - e0;
double lat = phi - VII * Math.pow(dE, 2) + VIII * Math.pow(dE, 4) - IX * Math.pow(dE, 6);
double lon = lambda0 + X * dE - XI * Math.pow(dE, 3) + XII * Math.pow(dE, 5) - XIIA * Math.pow(dE, 7);
// Convert latitude and longitude to degrees
double latitude = lat * (180 / Math.PI);
double longitude = lon * (180 / Math.PI);
// Add the corrections 1.6 metres is down from 24 metres inaccuracy
latitude += 0.000239504;
longitude -= 0.001692888;
// The single point is considered as the centroid
double[] centroid = { latitude, longitude };
// Output the results
System.out.println("Latitude: " + latitude);
System.out.println("Longitude: " + longitude);
System.out.println("Centroid: " + java.util.Arrays.toString(centroid));
}
}