EPSG:9817

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 compute the Lambert Conic Near-Conformal the following formulae are used. First compute constants for the projection:

n = f / (2-f)
A = 1 / (6 rhoO nuO) 
A’ = a [ 1- n + 5 (n^2 - n^3 ) / 4 + 81 ( n^4 - n^5 ) / 64]*pi /180 
B’ = 3 a [ n - n^2 + 7 ( n^3 - n^4 ) / 8 + 55 n^5 / 64] / 2
C’ = 15 a [ n^2 -n^3 + 3 ( n^4 - n^5 ) / 4 ] / 16
D’ = 35 a [ n^3 - n^4 + 11 n^5 / 16 ] / 48
E’ = 315 a [ n^4 - n^5 ] / 512
r0 = ko nu0 / tan(lat0)
s0 = A’ latO - B’ sin(2 latO) + C’ sin(4 latO) - D’ sin(6 latO) + E’ sin(8 latO) where in the first term latO is in degrees, in the other terms latO is in radians.

Then for the computation of easting and northing from latitude and longitude:

s = A’ lat - B’ sin(2 lat) + C’ sin(4 lat) - D’ sin(6 lat) + E’ sin(8 lat) where in the first term latO is in degrees, in the other terms latO is in radians.
M = s - sO
M = ko ( m + A m^3)
r = rO - M
theta = (lon - lonO) sin(latO)

and
E = FE + r sin(theta)
N = FN + M + r sin(theta) tan(theta/2)

The reverse formulas for latitude and longitude from Easting and Northing are:

theta' = atan2 {(E –  FE) , [rO –  (N –  FN)]} (see GN7-2 implementation notes in preface for atan2 convention)
r' = +/- {(E –  FE)^2 + [rO –  (N –  FN)]}^2}^0.5, taking the sign of latO 
M' = rO – r'

If an exact solution is required, it is necessary to solve for m and lat using iteration of the two equations:
m'= 	m' – [M' – ko m' – ko A (m')^3] / [– ko – 3 ko A (m')^2]
using M' for m' in the first iteration. This will usually converge (to within 1mm) in a single iteration. Then
lat' = lat' +{m' + sO – [A' lat' (180/pi) – B' sin(2 lat')  + C' sin(4 lat')  –  D' sin(6lat') + E' sin(8 lat')]}/A' (pi/180)
first using lat' = latO + m'/A' (pi/180).

However the following non-iterative solution is accurate to better than 0.001" (3mm) within 5 degrees latitude of the projection origin and should suffice for most purposes:
m' = 	M' – [M' ko M' – ko A (M')^3] / [– ko – 3 ko A (M')^2]
lat' = latO + m'/A' (pi/180)
s' = A	' lat' –  B' sin(2 lat')  + C' sin(4 lat')  –  D' sin(6 lat') + E' sin(8 lat')
		where in the first term lat' is in degrees, in the other terms lat' is in radians.
Ds' = 	A'(180 / pi) – 2B' cos(2 lat')  + 4C' cos(4 lat')  –  6D' cos(6 lat') + 8E' cos(8 lat')
lat = lat' – [(m' + sO – s') / (–ds')] radians

Then after solution of lat using either method above
lon = lonO + theta' / sin(latO) where lonO and lon are in radians.