Python trig precision problem

[Cary West]

Hello all, having a small problem with a trig routine using python math 
module. This code is designed to convert geodetic coordinates to lambert 
conformal conic coordinates. It implements the formulas found at 
http://mathworld.wolfram.com/LambertConformalConicProjection.html .  The 
problem is that, at least on my machine, the precision is off to the tune of 
around 1 km, which is unacceptable. (using the calculator found at 
http://www.deh.gov.au/erin/tools/geo2lam-gda.html ) Anyone have any ideas? 
here is the offending code.

from math import *

def lambert(lat,lon,ref_lat,ref_lon,st_parallel_1,st_parallel_2):

    earth_radius=  6371.9986
    lon*=pi/180.0
    lat *=pi/180.0

    ref_lon*=pi/180.0
    ref_lat *=pi/180.0

    st_parallel_1 *=pi/180.0
    st_parallel_2 *=pi/180.0

    def sec(theta):
        return 1.0/cos(theta)

    def cot(theta):
        return 1.0/tan(theta)

    n = log( cos( st_parallel_1 ) * sec( st_parallel_2 ) ) / ( log( tan ( pi 
/ 4.0 + st_parallel_2 / 2.0 ) * cot( pi / 4.0 + st_parallel_1 / 2.0 ) ) )

    F = ( cos( st_parallel_1 ) * tan ( pi / 4.0 + st_parallel_1 / 2.0 ) ** 
n ) / n

    rho = F * ( cot ( pi / 4.0 + lat / 2.0 ) ** n )

    ref_rho = F * ( cot ( pi / 4.0 + ref_lat / 2.0 ) ** n )

    x = rho * sin( n * ( lon-ref_lon ) )
    y = ref_rho - ( rho * cos( n * ( lon-ref_lon ) ) )

    return earth_radius*x,earth_radius*y
#######################################

lat,lon=35,-90

print lambert(lat,lon,40,-97,33,45)