from math import radians, cos, degrees
from math import radians, degrees, sin, cos, atan2, sqrt


def convert_to_cartesian(lat, lon, reference_point):
    """将经纬度转换为基于参考点的直角坐标，考虑地球椭球模型"""
    # 地球椭球参数（WGS84）
    a = 6378137.0  # 长半轴，单位：米
    f = 1 / 298.257223563  # 扁率
    e2 = 2 * f - f ** 2  # 第一偏心率平方

    # 提取参考点坐标
    ref_lat, ref_lon = reference_point

    # 转换成弧度
    lat_rad = radians(lat)
    lon_rad = radians(lon)
    ref_lat_rad = radians(ref_lat)
    ref_lon_rad = radians(ref_lon)

    # 计算曲率半径
    N = a / sqrt(1 - e2 * sin(ref_lat_rad) ** 2)  # 参考点处的卯酉圈曲率半径

    # 计算基于参考点的平面直角坐标
    delta_lon = lon_rad - ref_lon_rad
    X = (N + 0) * cos(ref_lat_rad) * delta_lon
    Y = (a * (1 - e2)) / (1 - e2 * sin(ref_lat_rad) ** 2) * (lat_rad - ref_lat_rad)

    return X, Y


def convert_to_geodetic(x, y, reference_point):
    """将直角坐标转换为经纬度，考虑地球椭球模型"""
    # 地球椭球参数（WGS84）
    a = 6378137.0  # 长半轴，单位：米
    f = 1 / 298.257223563  # 扁率
    e2 = 2 * f - f ** 2  # 第一偏心率平方

    # 提取参考点坐标
    ref_lat, ref_lon = reference_point

    # 转换成弧度
    ref_lat_rad = radians(ref_lat)
    ref_lon_rad = radians(ref_lon)

    # 计算曲率半径
    N = a / sqrt(1 - e2 * sin(ref_lat_rad) ** 2)  # 参考点处的卯酉圈曲率半径

    # 计算纬度
    lat_rad = y * (1 - e2 * sin(ref_lat_rad) ** 2) / (a * (1 - e2)) + ref_lat_rad

    # 计算经度
    if cos(ref_lat_rad) == 0:
        lon_rad = 0
    else:
        lon_rad = x / ((N + 0) * cos(ref_lat_rad)) + ref_lon_rad

    # 转换回角度
    lat = degrees(lat_rad)
    lon = degrees(lon_rad)

    return lat, lon