Projection: 3-degree Gauss-Kruger CM 108E

For Projected Coordinate System OSGB 1936 / British National Grid

Parameters:
Ellipsoid  Airy 1830  a = 6377563.396 m  1/f = 299.32496
then e'^2 = 0.00671534 and e^2 = 0.00667054

Latitude of natural origin (LatO) = 49°00'00"N = 0.85521133 rad
Longitude of natural origin (LonO) = 2°00'00"W = -0.03490659 rad
Scale factor (ko) = 0.9996013
False Eastings (FE) = 400000.00 m
False Northings (FN) = -100000.00 m

Forward calculation for: 
Latitude = 50°30'00.00"N = 0.88139127 rad
Longitude = 00°30'00.00"E = 0.00872665 rad
Constants of the projection:

n = 0.00167322
B = 6366914.609
h1 = 0.0008347452
h2 = 0.0000007554
h3 = 1.18487E-09
h4 = 2.40864E-12
QO = 0.9787671618
?O0 = 0.8518980373
?O1 = 0.0008273732
?O2 = -0.0000001986
?O3 = -1.0918E-09
?O4 = 1.2218E-12
Mo = 5429228.602
Q = 1.0191767215
? = 0.8781064142
?0 = 0.0278629616
?0 = 0.8785743280
?1 = -0.0000086229
?1 = 0.0008215669
?2 = -0.0000000786
?2 = -0.0000002768
?3 = 1.05551E-10
?3 = -1.01855E-09
?4 = 3.97791E-13
?4 = 1.67447E-12
? = 0.0278542603
? = 0.8793956171
Then
Easting	E = 577274.99 metres
Northing	N =  69740.50 metres

Reverse calculation for same easting and northing first gives:

h1' = 0.0008347455
h2' = 0.0000000586
h3' = 1.65563E-10
h4' = 2.13692E-13
Then
?' = 0.87939562
?' = 0.0278542603
?1' = 0.0008213109
?1' = -0.0000086953
?2' = -0.0000000217
?2' = -0.0000000061
?3' = -1.41881E-10
?3' = 1.486E-11
?4' = 1.49609E-13
?4' = 3.50657E-14
?0' = 0.8785743280
?0' = 0.0278629616
?' = 0.8781064142
Q' = 1.0191767215
Q" 1st iteration	= 1.0243166838
Q" 2nd iteration	= 1.0243306667
Q" 3rd iteration	= 1.0243307046
Q" 4th iteration	= 1.0243307047

Then
Latitude	(Lat) = 50°30'00.000"N
Longitude (Lon) =	00°30'00.000"E

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 the pole (within 2" or 50m), the tangent function in the equation for Qo is unstable and may fail. Mo may instead be calculated as:
Mo = a[(1 – e^2/4 – 3e^4/64 – 5e^6/256 –....)LatO – (3e^2/8 + 3e^4/32 + 45e^6/1024+....)sin2LatO 
+ (15e^4/256 + 45e^6/1024 +.....)sin4LatO – (35e^6/3072 + ....)sin6LatO  + .....]
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 ?')