EPSG:9813

Note : these formulas have been transcribed from IGN Document NT/G 74. Users are encouraged to use that document rather than the text which follows as reference because limitations in the transcription will be avoided.

From the defining parameters the following constants for the map projection may be calculated:

B = {1+[e^2 cos^4(phiC)]/(1– e^2)}^0.5
phiS = asin[sin(phiC) / B]
R = a  kC {(1–e^2)^0.5 / [1–e^2 sin^2(phiC)]}
C = ln[tan(pi/4+phiS /2)] – B. ln{tan(pi/4+phiC /2) ([1 – e sin(phiC)]/[1+e sin(phiC)])^(e/2)}

Forward case: To compute (E,N) from a given (lat,lon)
L = B.(lon–lonC)
q = C + B . ln{tan(pi/4+lat/2) ([1–e sin(lat)] / [1+e sin(lat)])^(e/2)}
P = 2.atan[e ^q] – pi/2   where e is the base of natural logarithms
U = cos(P).cos(L).cos(phiS) + sin(P).sin(phiS) 
V = cos(P).cos(L).sin(phiS) - sin(P).cos(phiS)
W = cos(P).sin(L)
d = (U^2+V^2)^0.5
if d <> 0 then L' = 2.atan(V/(U+d)) and P' = atan(W/d)
if d = 0 then L' = 0 and P' = sign(W).pi/2
H = –L' + i.ln(tan(pi/4+P'/2))                 where i^2 = –1
G = (1 – cos(2.alphaC) + i.sin(2.alphaC))/12
E = FE + R . IMAGINARY(H + G.H^3)
N = FN + R . REAL(H + G.H^3)

Reverse case: To compute (lat, lon) from a given (E,N):
G = (1–cos(2.alphaC) + i.sin(2.alphaC))/12   where i^2 = –1
To solve for Latitude and Longitude, a re-iterative solution is required, where the first two elements are
H0 = (N–FN)/R + i.(E–FE)/R  ie k = 0
H1 = H0/(H0 + G.H0^3),    i.e. k = 1,
and in subsequent reiterations, k increments by 1
Hk+1 = (H0+2.G.Hk^3)/(3.G.Hk^2+1)
Re-iterate until ABSOLUTE(REAL([H0-Hk-G.Hk^3)])) < 1E-11

L' = –1.REAL(Hk)
P' = 2.atan{ e ^[IMAGINARY(Hk)]} – pi/2  where e is the base of natural logarithms.
U' = cos(P').cos(L').cos(phiS) + cos(P').sin(L').sin(phiS)   
V' = sin(P')
W' = cos(P').cos(L').sin(phiS) – cos(P').sin(L').cos(phiS)
d = (U'^2+ V'^2)^0.5
if d <> 0 then L = 2 atan[V'/( U'+d)] and P = atan(W'/d)
if d = 0   then L = 0 and P = SIGN(W') . pi/2
lon = lonC + (L/B)

q' = {ln[tan(pi/4+P/2)] – C}/B
The final solution for latitude requires a second re-iterative process, where the first element is
lat'(0) = 2.atan(e ^q') – pi/2    where e is the base of natural logarithms.
And the subsequent elements are
lat'(k) = 2.atan{({1+e.sin[lat(k-1)]} / {1–e.sin[lat(k-1)]})^(e/2).e ^q'} – pi/2  for K =1 ?
Iterate until ABSOLUTE(lat(k)-lat(k-1))  < 1E-11	
lat = lat(k)