高斯-克里格投影反算Python代码，由x,y到经纬度

高斯-克里格投影反算Python代码，由x,y,中央子午线,东加常数到经纬度BL,由豆包AI帮我写出，经过我的调试与验证，结果正确。

import math

def gauss_krueger_inverse(x, y, central_meridian, false_easting=500000):
    """ 已知条件：x为纵坐标,y为横坐标，central_meridian为中央子午线经度，false_easting为东加常数 """
    """ 高斯-克吕格投影反算：将平面直角坐标(x,y)转换为经纬度(B,L)坐标 """    
    
    y = y - false_easting  # 去除东加常数
    
    # WGS84椭球体参数
    a = 6378137.0  # 长半轴
    f = 1/298.257223563  # 扁率
    e2 = 2*f - f*f  # 第一偏心率平方
    e4 = e2 * e2
    e6 = e4 * e2
    
    # 计算底点纬度Bf
    X = x
    mu = X / (a * (1 - e2/4 - 3*e4/64 - 5*e6/256))
    
    e1 = (1 - math.sqrt(1 - e2)) / (1 + math.sqrt(1 - e2))
    e1_squared = e1 * e1  # e1的平方
    e1_cubed = e1 * e1 * e1  # e1的立方
    e1_quad = e1_squared * e1_squared  # e1的四次方
    
    Bf = mu + (3*e1/2 - 27*e1_cubed/32) * math.sin(2*mu) + \
         (21*e1_squared/16 - 55*e1_quad/32) * math.sin(4*mu) + \
         (151*e1_cubed/96) * math.sin(6*mu) + \
         (1097*e1_quad/512) * math.sin(8*mu)
    
    # 计算辅助参数
    Nf = a / math.sqrt(1 - e2 * math.sin(Bf)**2)
    Mf = a * (1 - e2) / (1 - e2 * math.sin(Bf)**2)**(3/2)
    tf = math.tan(Bf)
    n2f = e2 * math.cos(Bf)**2 / (1 - e2)
    
    # 计算经纬度
    dy = y
    dy2 = dy * dy
    dy4 = dy2 * dy2
    dy6 = dy4 * dy2
    
    B = Bf - (tf / (2*Mf*Nf)) * dy2 + \
        (tf / (24*Mf*Nf**3)) * (5 + 3*tf**2 + n2f - 9*tf**2*n2f) * dy4 - \
        (tf / (720*Mf*Nf**5)) * (61 + 90*tf**2 + 45*tf**4) * dy6
    
    L = central_meridian * math.pi/180 + \
        (1 / (Nf * math.cos(Bf))) * dy - \
        (1 / (6*Nf**3 * math.cos(Bf))) * (1 + 2*tf**2 + n2f) * dy**3 + \
        (1 / (120*Nf**5 * math.cos(Bf))) * (5 + 28*tf**2 + 24*tf**4 + 6*n2f + 8*tf**2*n2f) * dy**5
    
    # 返回 (纬度B, 经度L) 单位为度
    return (B * 180/math.pi, L * 180/math.pi)