EPSG:9807

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.

For the calculation of easting and northing from latitude and longitude, first calculate constants for the projection:

n = f / (2-f)
B = [a/(1+n)] (1 + n^2/4 + n^4/64)
		
h1 = n/2 – (2/3)n^2 + (5/16)n^3 + (41/180)n^4
h2 = (13/48)n^2 – (3/5)n^3 + (557/1440)n^4
h3 = (61/240)n^3 – (103/140)n^4
h4 = (49561/161280)n^4

Then the meridional arc distance from equator to the projection origin (Mo) is computed from:
If LatO = 0 then Mo = 0
else if LatO ≡ 90°N ≡ π/2 radians
   Mo = B (π/2)
else if LatO ≡ 90°S ≡ -π/2 radians
   Mo = B (-π/2)	
else
 Qo = asinh(tan LatO) – [e atanh(e sin LatO)]
 βo = atan(sinh Qo)
 ξO0 = asin (sin βo)
Note: The previous two steps are taken from the generic calculation flow given below for latitude Lat, but here for LatO may be simplified to ξO0 = βo = atan(sinh Qo).
 ξO1 = h1 sin(2ξOo)
 ξO2 = h2 sin(4ξOo) 
 ξO3 = h3 sin(6ξOo) 
 ξO4 = h4 sin(8ξOo) 
 ξO = ξO0+ ξO1+ ξO2+ ξO3+ ξO4
 Mo = B ξO
end

Note: if the projection grid origin is very close to but not exactly at the pole (within 2" or 50m), in the equation above for Qo the tangent function is unstable and may fail. Rather than using Qo as above, Mo may instead be calculated from:

Mo = a[(1 – e^2/4 – 3e^4/64 – 5e^6/256 –....)LatO – (3e^2/8 + 3e^4/32 + 45e^6/1024+....)sin(2.LatO) 
+ (15e^4/256 + 45e^6/1024 +.....)sin(4.LatO) – (35e^6/3072 + ....)sin(6.LatO)  + .....]
with LatO in radians.

Then
Q = asinh(tan Lat) – [e atanh(e sin Lat)]
β = atan(sinh Q)
η0 = atanh [cos β sin(Lon – LonO)]
ξ0 = asin (sin β  cosh η0)	
ξ1 = h1 sin(2ξ0) cosh(2η0)
η1 = h1 cos(2ξ0) sinh(2η0)
ξ2 = h2 sin(4ξ0) cosh(4η0)
η2 = h2 cos(4ξ0) sinh(4η0)
ξ3 = h3 sin(6ξ0) cosh(6η0)
η3 = h3 cos(6ξ0) sinh(6η0)
ξ4 = h4 sin(8ξ0) cosh(8η0)
η4 = h4 cos(8ξ0) sinh(8η0)
ξ = ξ0 + ξ1 + ξ2 + ξ3 + ξ4
η = η0 + η1 + η2 + η3 + η4
and
Easting, E =  FE + ko B η
Northing, N =  FN + ko (B ξ – Mo)

For the reverse formulas to convert Easting and Northing projected coordinates to latitude and longitude first calculate constants of the projection where n is as for the forward conversion, as are B and Mo:
h1' = n/2 – (2/3)n^2 + (37/96)n^3 – (1/360)n^4
h2' = (1/48)n^2 + (1/15)n^3 – (437/1440)n^4
h3' = (17/480)n^3 – (37/840)n^4
h4' = (4397/161280)n^4

Then 
η' = (E –  FE) / (B ko)				
ξ' = [(N – FN) + ko Mo] / (B ko)
ξ1' = h1' sin(2ξ') cosh(2η')
η1' = h1' cos(2ξ') sinh(2η')
ξ2' = h2' sin(4ξ') cosh(4η')
η2' = h2' cos(4ξ') sinh(4η')
ξ3' = h3' sin(6ξ') cosh(6η')
η3' = h3' cos(6ξ') sinh(6η')
ξ4' = h4' sin(8ξ') cosh(8η')
η4' = h4' cos(8ξ') sinh(8η')
ξ0' = ξ' – (ξ1' + ξ2' + ξ3' + ξ4')
η0' = η' – (η1' + η2' + η3' + η4')
						
β' = asin(sin ξ0' / cosh η0')
Q' = asinh(tan β')
Q" = Q' + [e atanh(e tanh Q')] = Q' + [e atanh(e tanh Q")] which should be iterated until the change in Q" is insignificant. Then
Lat = atan(sinh Q")
Lon = LonO + asin(tanh(η0') / cos β')