Convert WGS84 latitude/longitude coordinates into OSGB36 easting/northing coordinates (as used in OS maps and gov.uk web sites)

wgs84LatLon2Osgb36EastNorth.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// Convert Latitude and Longitude to Easting and Northing coordinates
// Distilled from Visual Basic code published by the Ordnance Survey.

// Constants
#define DEG2RAD         (M_PI / 180)
#define OSGB36_A        6377563.396         // OSGB36 semi-major & semi-minor axis
#define OSGB36_B        6356256.909
#define WGS84_A         6378137.0           // WGS84 semi-major & semi-minor axis
#define WGS84_B         6356752.314245
#define F0              0.9996012717        // Central meridian scale

// OSGB origin (lat/lon and easting/northing)
#define PHI0            (49.0 * DEG2RAD)
#define LAM0            (-2.0 * DEG2RAD)
#define E0              400000.0
#define N0              -100000.0

// Helmert transform constants (WGS84 to OSGB36; reverse signs for inverse transform)
#define HELMERT_DX      -446.448
#define HELMERT_DY      125.157
#define HELMERT_DZ      -542.060
#define HELMERT_S       20.4894e-06
#define HELMERT_XROT    (-0.1502 / 3600 * DEG2RAD)
#define HELMERT_YROT    (-0.2470 / 3600 * DEG2RAD)
#define HELMERT_ZROT    (-0.8421 / 3600 * DEG2RAD)

void wgs84LatLon2Osgb36EastNorth(double wgs84_latitude, double wgs84_longitude, double *easting, double *northing) {

    // Step 1: Convert WGS84 Lat/Lon to geocentric XYZ coordinates
    // (We don't have a value for height; let's just guess 50 m)
    double phi = wgs84_latitude * DEG2RAD;
    double lam = wgs84_longitude * DEG2RAD;
    double h = 50.0;

    double sinPhi = sin(phi);
    double sinPhi2 = sinPhi * sinPhi;
    double cosPhi = cos(phi);
    double sinLam = sin(lam);
    double cosLam = cos(lam);

    double wgs84_e2 = (WGS84_A * WGS84_A - WGS84_B * WGS84_B) / (WGS84_A * WGS84_A);
    double wgs84_v = WGS84_A / sqrt(1 - (wgs84_e2 * sinPhi2));
    double x = (wgs84_v + h) * cosPhi * cosLam;
    double y = (wgs84_v + h) * cosPhi * sinLam;
    double z = ((wgs84_v * (1 - wgs84_e2)) + h) * sinPhi;

    // Step 2: Perform a Helmert transformation to convert these
    // geocentric coordinates into the system used by OSGB36
    double helmertX = x + x * HELMERT_S - y * HELMERT_ZROT+ z * HELMERT_YROT + HELMERT_DX;
    double helmertY = x * HELMERT_ZROT + y + y * HELMERT_S - z * HELMERT_XROT + HELMERT_DY;
    double helmertZ = -x * HELMERT_YROT + y * HELMERT_XROT + z + z * HELMERT_S + HELMERT_DZ;

    x = helmertX;
    y = helmertY;
    z = helmertZ;

    // Step 3: Convert these geocentric coordinates back into Lat/Long
    // values in the OSGB36 coordinate system
    double osgb36_e2 = ((OSGB36_A * OSGB36_A) - (OSGB36_B * OSGB36_B)) / (OSGB36_A * OSGB36_A);
    double rootXY2 = hypot(x, y);
    double phi1, sinPhi1, osgb36_v;
    double phi2 = atan(z / (rootXY2 * (1 - osgb36_e2)));
    
    do {
        phi1 = phi2;
        sinPhi1 = sin(phi1);
        osgb36_v = OSGB36_A / sqrt(1 - osgb36_e2 * sinPhi1 * sinPhi1);
        phi2 = atan((z + (osgb36_e2 * osgb36_v * sinPhi1)) / rootXY2);
    } while (fabs(phi1 - phi2) > 1e-09);
    phi = phi2;

    if (x >= 0) {
        lam = atan(y / x);
    }
    else if (x < 0 && y >= 0) {
        lam = atan(y / x) + M_PI;
    }
    else if (x < 0 && y < 0) {
        lam = atan(y / x) - M_PI;
    }

    // Step 4: Map these Lat/Lon values to the Ordnance Survey grid
    sinPhi = sin(phi);
    sinPhi2 = sinPhi * sinPhi;
    cosPhi = cos(phi);
    sinLam = sin(lam);
    cosLam = cos(lam);

    double tanPhi = tan(phi);
    double cosPhi3 = cosPhi * cosPhi * cosPhi;
    double cosPhi5 = cosPhi3 * cosPhi * cosPhi;
    double tanPhi2 = tanPhi * tanPhi;
    double tanPhi4 = tanPhi2 * tanPhi2;

    double af0 = OSGB36_A * F0;
    double bf0 = OSGB36_B * F0;

    double n = (af0 - bf0) / (af0 + bf0);
    double nu = af0 / sqrt(1 - osgb36_e2 * sinPhi2);
    double rho = (nu * (1 - osgb36_e2)) / (1 - osgb36_e2 * sinPhi2);
    double eta2 = (nu / rho) - 1;
    double p = lam - LAM0;
    double n2 = n * n;
    double n3 = n2 * n;
    double m = bf0 * (((1 + n + 5 * n2 / 4 + 5 * n3 / 4) * (phi - PHI0)) - ((3 * n + 3 * n2 + (21 * n3 / 8))
           * sin(phi - PHI0) * cos(phi + PHI0)) + ((15 * n2 / 8 + 15 * n3 / 8) * sin(2 * (phi - PHI0))
           * cos(2 * (phi + PHI0))) - (35 * n3 * sin(3 * (phi - PHI0)) * cos(3 * (phi + PHI0)) / 24));
    
    double i = m + N0;
    double ii = nu * sinPhi * cosPhi / 2;
    double iii = (nu * sinPhi * cosPhi3 / 24) * (5 - tanPhi2 + 9 * eta2);
    double iiia = (nu * sinPhi * cosPhi5 / 720) * (61 - 58 * tanPhi2 + tanPhi4);
    double iv = nu * cosPhi;
    double v = nu * cosPhi3 * (nu / rho - tanPhi2) / 6;
    double vi = nu * cosPhi5 * (5 - 18 * tanPhi2 + tanPhi4 + 14 * eta2 - 58 * tanPhi2 * eta2) / 120;
    
    double p2 = p * p;
    double p3 = p2 * p;
    double p4 = p2 * p2;
    double p5 = p3 * p2;
    double p6 = p3 * p3;

    *easting = E0 + p * iv + p3 * v + p5 * vi;
    *northing = i + p2 * ii + p4 * iii + p6 * iiia;
}


int main(int argc, char *argv[]) {
    double latitude, longitude, easting, northing;
    if (argc != 3) {
        return puts("Please provide latitude and longitude as command line arguments");
    }
    latitude = atof(argv[1]);
    longitude = atof(argv[2]);
    wgs84LatLon2Osgb36EastNorth(latitude, longitude, &easting, &northing);
    printf("Easting: %lf\nNorthing: %lf\n", easting, northing);
    return 0;
}