EPSG:9809

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.

Given the geodetic origin of the projection at the tangent point (lat0, lon0), the parameters defining the conformal sphere are:

R= sqrt( rho0 * nu0)
n= {1 + [e^2 * cos^4(lat0) / (1 - e^2)]}^0.5
c=  [(n+sin(lat0)) (1-sin(chi00))]/[(n-sin(lat0)) (1+sin(chi00))]

where:
sin(chi00) = (w1-1)/(w1+1)
w1 = (S1.(S2)^e)^n
S1 = (1+sin(lat0))/(1-sin(lat0))
S2 = (1-e sin(lat0))/(1+e sin(lat0))

The conformal latitude of the origin (chi0) is then computed from:

chi0 = asin[(w2-1)/(w2+1)]

where w2 = c w1
 
For any point with geodetic latitude (lat) the equivalent conformal latitude (chi) is then computed from: 
chi = asin[(w-1)/(w+1)]

where w = c (Sa (Sb)^e)^n
Sa = (1+sin(lat))/(1-sin(lat))
Sb = (1-e.sin(lat))/(1+e.sin(lat))
 
Then
dLambda = Lambda - Lamda0 = n(lon-lon0) where Lambda is the conformal longitude at geodetic longitude (lon)
B = [1+sin(chi) sin(chi0) + cos(chi) cos(chi0) cos(dLambda)]

and
E = FE + 2 R k0 cos(chi) sin(Lambda-Lambda0) / B
N = FN + 2 R k0 [sin(chi) cos(chi0) - cos(chi) sin(chi0) cos(dLambda)] / B

The reverse formulae to compute the geodetic coordinates from the grid coordinates involves computing the conformal values, then the isometric latitude and finally the geodetic values.

The parameters of the conformal sphere and conformal latitude at the origin are computed as above. Then for any point with Stereographic grid coordinates (E,N):

chi = chi0 + 2 atan[{(N-FN)-(E-FE) tan (j/2)} / (2 R k0)]

where 
g = 2 R k0 tan(pi/4 - chi0/2)
h = 4 R k0 tan(chi0) + g
i = atan2[(E-FE) , {h+(N-FN)}]
j = atan2[(E-FE) , (g-(N-FN)] - i
(see GN7-2 implementation notes for atan2 convention)

Geodetic longitude lon = (Lambda-Lambda0 ) / n +  Lambda0

Then isometric latitude psi = 0.5 ln [(1+ sin(chi)) / { c (1-  sin(chi))}] / n

The first approximation for geodetic latitude is found from:

lat1 = 2 atan(e^psi)  - pi/2  

where e=base of natural logarithms
            psi = isometric latitude at lat1
                   = ln[{tan(lat1 / 2 + pi/4}  {(1-e sin(lat1))/(1+e sin(lat1))}^(e/2)]
 
Then iterate: 
psi1 = ln[{tan(lati/2 + pi/4}  {(1-e sin(lati))/(1+e sin(lati))}^(e/2)]
and
lat(i+1) = lati - ( psii - psi ) cos(lati) (1 -e^2 sin^2(lati)) / (1 - e^2)

until the change in lat is sufficiently small.

For Oblique Stereographic projections centred on points in the southern hemisphere,  the signs of E, N, lon0 and lon must be reversed to be used in the equations and lat as a southerly latitude will be negative anyway.

If the projection is the equatorial case, lat0 and chi0 will be zero degrees and the formulas may be simplified as a result, but the above formulas remain valid.

For the polar case, lat0 and chi0 will be 90 degrees and the formulas become indeterminate. See polar steroegraphic method for formulas for the polar case.

An alternative approach is given by Snyder, where, instead of defining a single conformal sphere at the origin point, the conformal latitude at each point on the ellipsoid is computed.  The conformal longitude is then always equivalent to the geodetic longitude.  This approach is a valid alternative to the above, but gives slightly different results away from the origin point. It is therefore considered by EPSG to be a different coordinate operation method to that described above.