EPSG:9812

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.

The following constants for the projection may be calculated :

B = {1 + [esq * cos^4(latc) / (1 - esq )]}^0.5
A = a * B * kc *(1 - esq )^0.5 / ( 1 - esq * sin^2(latc))
to = tan(pi/4 - latc/2) / ((1 - e*sin(latc)) / (1 + e*sin(latc)))^(e/2)
D = B (1 - esq)^0.5  / (cos(latc) * ( 1 - esq*sin^2(latc))^0.5)
if D < 1 to avoid problems with computation of F make D^2  = 1 
F = D + (D^2 - 1)^0.5  * SIGN(latc)
H = F*(to)^B
G = (F - 1/F) / 2
gammao = asin(sin(alphac) / D)
lonO = lonc - (asin(G*tan(gammao))) / B

Forward case: To compute (E,N) from a given (lat,lon) :

t = tan(pi/4 - lat/2) / ((1 - e sin (lat)) / (1 + e sin (lat)))^(e/2)
Q = H / t^B
S = (Q - 1 / Q) / 2
T = (Q + 1 / Q) / 2
V = sin(B (lon - lonO))
U = (- V cos(gammao) + S sin(gammao)) / T
v = A ln((1 - U) / (1 + U)) / 2 B
u = A atan2((S cos(gammao) + V sin(gammao)) , cos(B (lon - lonO))) / B  (see GN7-2 implementation notes in preface for atan2 convention)

The rectified skew co-ordinates are then derived from:
E = v cos(gammac) + u sin(gammac) + FE
N = u cos(gammac) - v sin(gammac) + FN

Reverse case: Compute (lat,lon)  from a given (E,N)  :

v’ = (E - FE) cos(gammac) - (N - FN) sin(gammac)
u’ = (N - FN) cos(gammac) + (E - FE) sin(gammac)

Q’ = e^- (B v ‘/ A)  where e is the base of natural logarithms.
S' = (Q’ - 1 / Q’) / 2
T’ = (Q’ + 1 / Q’) / 2
V’ = sin (B u’ / A)
U’ = (V’ cos(gammac) + S’ sin(gammac)) / T’
t’ = (H / ((1 + U’) / (1 - U’))^0.5)^(1 / B)

chi = pi / 2 - 2 atan(t’)

lat = chi + sin(2chi).( e^2 / 2 + 5*e^4 / 24 + e^6 / 12 + 13*e^8 / 360) +  sin(4*chi).( 7*e^4 /48 + 29*e^6 / 240 + 811*e8 / 11520) +  sin(6chi).( 7*e^6 / 120 + 81*e8 / 1120) +  sin(8chi).(4279 e^8 / 161280)

lon = lonO - atan2 ((S’ cos(gammao) - V’ sin(gammao)) , cos(B*u’ / A)) / B