EPSG:9837

Note: These formulas have been transcribed from EPSG Guidance Note #7-2. Users are encouraged to use that document rather than the text which follows as reference because limitations in the transcription will be avoided.

To convert latitude P, longitude L and ellipsoidal height h into topocentric coordinates U,V,W: 

U = (nu + h) cos(P) sin (L – Lo)
V = (nu + h) [sin(P) cos(Po) – cos(P) sin(Po) cos(L – Lo)] + e^2 (nuO sin(Po) – nu sin(P)) cos(Po) 
W = (nu + h) [sin(P) sin(Po) + cos(P) cos(Po) cos(L – Lo)] + e^2 (nuO sin(Po) – nu sin(P)) sin(Po) – (nuO – ho)

where Po, Lo and ho are the ellipsoidal coordinates of the topocentric origin

and nu is the radius of curvature in the prime vertical at latitude P = a /(1 – e^2 sin^2 P)^0.5
nuO is the radius of curvature in the prime vertical at latitude Po = a /(1 – e^2 sin^2 Po)^0.5
	e is the eccentricity of the ellipsoid here e^2 = (a^2 – b^2)/a^2 = 2f – f^2

The reverse formulae to convert topocentric coordinates (U, V, W) into latitude, longitude and ellipsoidal height (P, L, h) first draws on the reverse case of method 9836 to derive geocentric coordinates X, Y, Z and then on the reverse case of method 9602 to derive latitude, longitude and height.

First,
X = Xo – U sin(Lo) – V sin(Po) cos(Lo) + W cos(Po) cos(Lo) 
Y = Yo + U cos(Lo) – V sin(Po) sin(Lo) + W cos(Po) sin(Lo)
Z = Zo + V cos(Po) + W sin(Po)

where,
XO = (nuO + hO) cos(Po) cos(Lo) 
YO = (nuO + hO) cos(Po) sin(Lo)
ZO = [(1 – e2) nuO + hO] sin(Po)
Po, Lo, hO are the ellipsoidal coordinates of the topocentric origin,
nuO is the radius of curvature in the prime vertical at latitude Po = a /(1 – e^2 sin^2(Po))^0.5, and
e is the eccentricity of the ellipsoid where e^2 = (a^2 – b^2)/a^2 = 2f – f^2.

Then,
lat = atan2[(Z + ε b sin^3(q)) , (p – e^2 a cos^3(q))]
lon = atan2 (Y , X)
where
	ε = e^2 / (1 – e^2)
	b = a(1 – f)
	p = (X^2 + Y^2)^0.5
	q = atan2[(Z a) , (p b)]
	L is relative to the Greenwich prime meridian.

and
h = (p / cos(P)) – nu
where
nu is the radius of curvature in the prime vertical at latitude lat = a /(1 – e^2sin^2(P))^0.5
(see implementation notes in preface for atan2 convention)