Latitude/longitude conversion

# Given a latitude and longitude as parameters, this function

# creates a returns a Point (from the graphics library). It is

# only intended to work for points in Arkansas (though it should

# do reasonably well for any point in the central US). It uses

# the UTM coordinate system.

def latlong_to_xy(lat, long):

    a =  6377397               # Uses Bessel 1841 ellipsoid

    ecc_sqr =  0.006674372

    k0 = 0.9996



    lat = math.radians(lat)

    long = math.radians(long)

    long_origin = -1.623156    # This changes for other regions of the world



    ecc_sqr2 = ecc_sqr ** 2

    ecc_sqr3 = ecc_sqr ** 3

    ecc_sqr_prime = (ecc_sqr)/(1-ecc_sqr)



    N = a / math.sqrt(1 - ecc_sqr * math.sin(lat)**2)

    T = math.tan(lat)**2

    C = ecc_sqr_prime * math.cos(lat)**2

    A = math.cos(lat) * (long - long_origin)



    M = a*((1   - ecc_sqr/4     - 3*ecc_sqr2/64 - 5*ecc_sqr3/256)*lat 

                - (3*ecc_sqr/8  + 3*ecc_sqr2/32 + 45*ecc_sqr3/1024)*math.sin(2*lat)

                                    + (15*ecc_sqr2/256 + 45*ecc_sqr3/1024)*math.sin(4*lat) 

                                    - (35*ecc_sqr3/3072)*math.sin(6*lat))

    

    easting = (k0*N*(A+(1-T+C)*A**3/6

                    + (5-18*T+T*T+72*C-58*ecc_sqr_prime)*A**5/120)

                    + 500000.0)



    northing = (k0*(M+N*math.tan(lat)*(A*A/2+(5-T+9*C+4*C*C)*A**4/24

                 + (61-58*T+T*T+600*C-330*ecc_sqr_prime)*A**6/720)))

    if lat < 0:

        northing += 10000000.0 # 10000000 meter offset for southern hemisphere

    return Point((easting - 357398.24) * 0.001, 381.62 - (northing - 3656759.44) * 0.001)