"""大地坐标转墨卡托投影坐标"""
import math

#######################################################
# 计算参数
wd0 = math.pi / 180 * 24  # 基准纬度
wd = math.pi / 180 * (23 + 9 / 60)  # 经度
jd = math.pi / 180 * (119 + 14 / 60 + 23 / 3600)  # 纬度
m = 500  # 比例尺分母*100图上1cm代表实地多少米
########################################################
# 椭球参数
a = 6378137  # 椭圆的长半轴
b = 6356752.31414  # 椭圆的短半轴,准确应为6,356,752.3141403558478521068615295
e1 = math.sqrt(a ** 2 - b ** 2) / a  # 椭圆的第一偏心率0.0818191910428#
e2 = math.sqrt(a ** 2 - b ** 2) / b  # 椭圆的第二偏心率
f = a / (a - b)  # 扁率 f=1/298.257222101


def get_NB0(wd0):
    """计算卯酉圈曲率半径, wd0:基准维度"""
    _1 = a * a / b
    _2 = math.sqrt(1 + (e2 ** 2) * (math.cos(wd0) ** 2))
    _3 = _1 / _2
    return _3


NB0 = get_NB0(wd0)
k = NB0 * math.cos(wd0)  # 纬度圈半径

print('e1', e1)
print('e2', e2)
print('f', f)
print('NB0', NB0)
print('k  ', k)
######################################################
x = k * jd / m
t1 = math.tan(math.pi / 4 + wd / 2)
t2 = ((1 - e1 * math.sin(wd)) / (1 + e1 * math.sin(wd))) ** (e1 / 2)
y = k / m * math.log(t1 * t2, math.e)
print(t1, t2, math.log(math.e))
print('y', x)
print('y', y)
