EPSG:1041

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.

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

A = a(1 - e^2)^0.5 / [1 - e^2 sin^2.(latC)]
B = {1 + [e^2 * cos^4(latC) / (1 - e^2)]}^0.5
gammao = asin[sin(latC) / B]
to = tan(pi/4 + gammao/2).[(1 + e sin(latC)) / (1 - e sin(latC))]^(e.B/2) / [tan(pi/4 + latC/2)]^B
n = sin(latp)
ro = kp.A / tan(latp)

To derive the projected Southing and Westing coordinates of a point with geographical coordinates (lat, lon) the formulas for the Krovak are:

U = 2(atan{to.tan^B(lat/2 + pi/4) / [(1 + e sin(lat)) / (1 - e sin(lat))]^[e.B/2]} - pi/4)
V = B(lonO - lon) where lonO and lon must both be referenced to the same prime meridian.
T = asin[cos(alphaC).sin(U) + sin(alphaC).cos(U). cos(V)]
D = asin[cos(U).sin(V)/cos(T)]
theta = n.D
r = ro.tan^n(pi/4 + latp/2) / tan^n(T/2 + pi/4)
Xp = r.cos(theta)
Yp = r.sin(theta)
Southing = Xp + FN
Westing = Yp + FE

Then Easting = -Westing
     Northing = -Southing

Note also that the formula for D is satisfactory for the normal use of the projection within the pseudo-longitude range on the conformal sphere of ±90 degrees from the central line of the projection. Should there be a need to exceed this range (which is not necessary for application in Czechia and Slovakia) then for the calculation of D use:
sin(D1) = cos(U) * sin(V) / cos(T)
cos(D1) = {[cos(alphaC)*sin(T) – sin(U)] / [sin(alphaC)*cos(T)]}
D =  atan2(sin(D1) , cos(D1))  (see implementation notes in preface for atan2 convention)

The reverse formulas to derive the latitude and longitude of a point from its Easting and Northing values are:
Southing = -Northing
Westing = -Easting
Xp' = Southing – FN
Yp' = Westing – FE
r' = [(Yp')^2 + (Xp')^2]^(1/2)  
theta' = atan2[Yp' , Xp'] (see GN7-2 implementation notes in preface for atan2 convention)
D' = theta' / sin(latp)
T' = 2{atan[((ro / r')^(1/n)).tan(pi/4 + latp/2)] - pi/4}
U' = asin[cos(alphaC).sin(T') - sin(alphaC).cos(T').cos(D')]
V' = asin(cos(T').sin(D') / cos(U'))

Then latitude lat is found by iteration using U' as the value for lat(j-1) in the first iteration:
lat(j) = 2*(atan{to^(-1/B) tan^(1/B).(U'/2 + pi/4).[(1 + e sin(lat(j-1)) / (1 - e sin(lat(j-1))]^(e/2)} - pi/4) 

Then
lon = lonO - V' / B where lon is referenced to the same prime meridian as lonO.