def WGS84toOSGB36(lat, lon):
<EOL>lat_1 = lat*pi/<NUM_LIT><EOL>lon_1 = lon*pi/<NUM_LIT><EOL>a_1, b_1 = <NUM_LIT>, <NUM_LIT><EOL>e2_1 = <NUM_LIT:1> - (b_1*b_1)/(a_1*a_1)  <EOL>nu_1 = a_1/sqrt(<NUM_LIT:1>-e2_1*sin(lat_1)**<NUM_LIT:2>)<EOL>H = <NUM_LIT:0>  <EOL>x_1 = (nu_1 + H)*cos(lat_1)*cos(lon_1)<EOL>y_1 = (nu_1 + H)*cos(lat_1)*sin(lon_1)<EOL>z_1 = ((<NUM_LIT:1>-e2_1)*nu_1 + H)*sin(lat_1)<EOL>s = <NUM_LIT>*<NUM_LIT:10>**-<NUM_LIT:6>  <EOL>tx, ty, tz = -<NUM_LIT>, <NUM_LIT>, -<NUM_LIT><EOL>rxs, rys, rzs = -<NUM_LIT>, -<NUM_LIT>, -<NUM_LIT><EOL>rx, ry, rz = rxs*pi/(<NUM_LIT>*<NUM_LIT>), rys*pi/(<NUM_LIT>*<NUM_LIT>), rzs*pi/(<NUM_LIT>*<NUM_LIT>)<EOL>x_2 = tx + (<NUM_LIT:1>+s)*x_1 + (-rz)*y_1 + (ry)*z_1<EOL>y_2 = ty + (rz)*x_1 + (<NUM_LIT:1>+s)*y_1 + (-rx)*z_1<EOL>z_2 = tz + (-ry)*x_1 + (rx)*y_1 + (<NUM_LIT:1>+s)*z_1<EOL>a, b = <NUM_LIT>, <NUM_LIT><EOL>e2 = <NUM_LIT:1> - (b*b)/(a*a)  <EOL>p = sqrt(x_2**<NUM_LIT:2> + y_2**<NUM_LIT:2>)<EOL>lat = atan2(z_2, (p*(<NUM_LIT:1>-e2)))  <EOL>latold = <NUM_LIT:2>*pi<EOL>while abs(lat - latold) > <NUM_LIT:10>**-<NUM_LIT:16>:<EOL><INDENT>lat, latold = latold, lat<EOL>nu = a/sqrt(<NUM_LIT:1>-e2*sin(latold)**<NUM_LIT:2>)<EOL>lat = atan2(z_2+e2*nu*sin(latold), p)<EOL><DEDENT>lon = atan2(y_2, x_2)<EOL>H = p/cos(lat) - nu<EOL>F0 = <NUM_LIT>  <EOL>lat0 = <NUM_LIT>*pi/<NUM_LIT>  <EOL>lon0 = -<NUM_LIT:2>*pi/<NUM_LIT>  <EOL>N0, E0 = -<NUM_LIT>, <NUM_LIT>  <EOL>n = (a-b)/(a+b)<EOL>rho = a*F0*(<NUM_LIT:1>-e2)*(<NUM_LIT:1>-e2*sin(lat)**<NUM_LIT:2>)**(-<NUM_LIT>)<EOL>eta2 = nu*F0/rho-<NUM_LIT:1><EOL>M1 = (<NUM_LIT:1> + n + (<NUM_LIT:5>/<NUM_LIT:4>)*n**<NUM_LIT:2> + (<NUM_LIT:5>/<NUM_LIT:4>)*n**<NUM_LIT:3>) * (lat-lat0)<EOL>M2 = (<NUM_LIT:3>*n + <NUM_LIT:3>*n**<NUM_LIT:2> + (<NUM_LIT>/<NUM_LIT:8>)*n**<NUM_LIT:3>) * sin(lat-lat0) * cos(lat+lat0)<EOL>M3 = ((<NUM_LIT:15>/<NUM_LIT:8>)*n**<NUM_LIT:2> + (<NUM_LIT:15>/<NUM_LIT:8>)*n**<NUM_LIT:3>) * sin(<NUM_LIT:2>*(lat-lat0)) * cos(<NUM_LIT:2>*(lat+lat0))<EOL>M4 = (<NUM_LIT>/<NUM_LIT>)*n**<NUM_LIT:3> * sin(<NUM_LIT:3>*(lat-lat0)) * cos(<NUM_LIT:3>*(lat+lat0))<EOL>M = b * F0 * (M1 - M2 + M3 - M4)<EOL>I = M + N0<EOL>II = nu*F0*sin(lat)*cos(lat)/<NUM_LIT:2><EOL>III = nu*F0*sin(lat)*cos(lat)**<NUM_LIT:3>*(<NUM_LIT:5> - tan(lat)**<NUM_LIT:2> + <NUM_LIT:9>*eta2)/<NUM_LIT><EOL>IIIA = nu*F0*sin(lat)*cos(lat)**<NUM_LIT:5>*(<NUM_LIT> - <NUM_LIT>*tan(lat)**<NUM_LIT:2> + tan(lat)**<NUM_LIT:4>)/<NUM_LIT><EOL>IV = nu*F0*cos(lat)<EOL>V = nu*F0*cos(lat)**<NUM_LIT:3>*(nu/rho - tan(lat)**<NUM_LIT:2>)/<NUM_LIT:6><EOL>VI = nu*F0*cos(lat)**<NUM_LIT:5>*(<NUM_LIT:5> - <NUM_LIT>*tan(lat)**<NUM_LIT:2> + tan(lat)**<NUM_LIT:4> + <NUM_LIT>*eta2 - <NUM_LIT>*eta2*tan(lat)**<NUM_LIT:2>)/<NUM_LIT><EOL>N = I + II*(lon-lon0)**<NUM_LIT:2> + III*(lon-lon0)**<NUM_LIT:4> + IIIA*(lon-lon0)**<NUM_LIT:6><EOL>E = E0 + IV*(lon-lon0) + V*(lon-lon0)**<NUM_LIT:3> + VI*(lon-lon0)**<NUM_LIT:5><EOL>return E, N<EOL>
Accept latitude and longitude as used in GPS.
    Return OSGB grid coordinates: eastings and northings.

    Usage:
    >>> from latlon_to_bng import WGS84toOSGB36
    >>> WGS84toOSGB36(51.4778, -0.0014)
    (538890.1053365842, 177320.49650700082)
    >>> WGS84toOSGB36(53.50713, -2.71766)
    (352500.19520169357, 401400.01483428996)