Here's something in python. Note it doesn't use the Z coordinate:
R = 6360000 # earth radius
(x, y, z) = xy # does not use `z`
r = math.sqrt(x**2 + y**2)
long =180 * math.atan2(y,x)/math.pi
lat = 180 * math.acos(r/R)/math.pi
return (lat, long)
To test on a known point, I've take the location of the statue of Eric Morecambe in Morecambe, which is approximately at sea level:
eric_google = (54.072878,-2.8680394)
The online converter returns the following OSGB and Cartesian coordinates for these decimal lat-long coordinates:
eric_osgb = (343296.36, 464452.65)
eric_xyz = (3745951.536, -187666.857, 5141507.615)
Test script is:
print "Eric according to google ", eric_google
print "Eric in cartesian ", eric_xyz
print "Eric converted ",cc_to_ll(eric_xyz)
The code above returns:
Eric according to google (54.072878, -2.8680394)
Eric in cartesian (3745951.536, -187666.857, 5141507.615)
Eric converted (53.862520806716375, -2.868039396573805)
The converted longitude is exact to 6dp, but latitude is off by a bit. This is because the code assumes a spherical earth, and a point on the surface of the sphere. Hence its dependent on the radius,
R, and the ellipsoidal shape used in the EPSG 4326 coordinate system.
I don't think you can get 10m precision without a Z coordinate in the cartesian system. The Z axis is everywhere parallel to the earth's axis, and in the UK that is at an angle of about 37 degrees, so a difference in Z coordinate of 100m will result in a shift of about 30m.
I think a better solution might be obtained by working out the intersection of the Z axis-aligned point with the EPSG 4326 ellipsoid rather than a sphere.