Lambert-93 to WGS84 converter

lam93toLatLon.py

def lam93toLatLon(x, y):
    constants = {
        'GRS80E': 0.081819191042816,
        'LONG_0': 3,
        'XS': 700000,
        'YS': 12655612.0499,
        'n': 0.7256077650532670,
        'C': 11754255.4261
    }

    delX = x - constants['XS']
    delY = y - constants['YS']
    gamma = math.atan(-delX / delY)
    R = math.sqrt(delX * delX + delY * delY)
    latiso = math.log(constants['C'] / R) / constants['n']
    sinPhiit0 = math.tanh(latiso + constants['GRS80E'] * math.atanh(constants['GRS80E'] * math.sin(1)))
    sinPhiit1 = math.tanh(latiso + constants['GRS80E'] * math.atanh(constants['GRS80E'] * sinPhiit0))
    sinPhiit2 = math.tanh(latiso + constants['GRS80E'] * math.atanh(constants['GRS80E'] * sinPhiit1))
    sinPhiit3 = math.tanh(latiso + constants['GRS80E'] * math.atanh(constants['GRS80E'] * sinPhiit2))
    sinPhiit4 = math.tanh(latiso + constants['GRS80E'] * math.atanh(constants['GRS80E'] * sinPhiit3))
    sinPhiit5 = math.tanh(latiso + constants['GRS80E'] * math.atanh(constants['GRS80E'] * sinPhiit4))
    sinPhiit6 = math.tanh(latiso + constants['GRS80E'] * math.atanh(constants['GRS80E'] * sinPhiit5))

    longRad = math.asin(sinPhiit6)
    latRad = gamma / constants['n'] + constants['LONG_0'] / 180 * math.pi

    longitude = latRad / math.pi * 180
    latitude = longRad / math.pi * 180

    return (latitude, longitude)