北东地坐标系和经纬度坐标系相互转换的c++程序

北东地坐标系（简称“平面坐标系”）与经纬度坐标系（简称“地理坐标系”）之间的转换需要用到大地测量学中的一些公式，包括高斯投影、大地坐标系和空间直角坐标系之间的转换等。下面是一个基于 WGS84 参考椭球体的平面坐标系和地理坐标系之间相互转换的 C++ 程序：

#include <iostream>
#include <cmath>

using namespace std;

const double PI = 3.14159265358979323846;
const double a = 6378137.0;
const double f = 1 / 298.257223563; // WGS84椭球体参数

// 计算椭球体的第一偏心率
double e1(double f) {
    return sqrt(2 * f - f * f);
}

// 计算椭球体的第二偏心率
double e2(double f) {
    return sqrt(f * (2 - f) / (1 - f) / (1 - f));
}

// 经纬度坐标系转平面坐标系
void geo2plane(double B, double L, double B0, double L0, double &N, double &E) {
    double e12 = e1(f) * e1(f);
    double e22 = e2(f) * e2(f);
    double L_L0 = L - L0;
    double sinB = sin(B);
    double cosB = cos(B);
    double tanB = tan(B);
    double N1 = a / sqrt(1 - e12 * sinB * sinB);
    double t2 = tanB * tanB;
    double N2 = N1 / sqrt(1 + e22 * cosB * cosB);
    double M = a * ((1 - e12 / 4 - 3 * e12 * e12 / 64 - 5 * e12 * e12 * e12 / 256) * B
             - (3 * e12 / 8 + 3 * e12 * e12 / 32 + 45 * e12 * e12 * e12 / 1024) * sin(2 * B)
             + (15 * e12 * e12 / 256 + 45 * e12 * e12 * e12 / 1024) * sin(4 * B)
             - (35 * e12 * e12 * e12 / 3072) * sin(6 * B));
    N = M + N2 * tanB * (L_L0 * L_L0 / 2 + (5 - t2 + 9 * e22 * cosB * cosB + 4 * e22 * e22 * cosB * cosB * cosB * cosB) * L_L0 * L_L0 * L_L0 * L_L0 / 24
         + (61 - 58 * t2 + t2 * t2 + 600 * e22 * cosB * cosB + 324 * e22 * e22 * cosB * cosB * cosB * cosB) * L_L0 * L_L0 * L_L0 * L_L0 * L_L0 * L_L0 / 720);
    E = N2 * L_L0 + N2 * tanB * (1 + t2 + e22 * cosB * cosB) * L_L0 * L_L0 * L_L0 / 6
        + N2 * tanB * (5 - 18 * t2 + t2 * t2 + 72 * e22 * cosB * cosB - 58 * e22 * e22 * cosB * cosB * cosB * cosB) * L_L0 * L_L0 * L_L0 * L_L0 * L_L0 / 120;
}

// 平面坐标系转经纬度坐标系
void plane2geo(double N, double E, double B0, double L0, double &B, double &L) {
    double e12 = e1(f) * e1(f);
    double e22 = e2(f) * e2(f);
    double n = (a - sqrt(a * a - N * N)) / (a * f);
    double sinB = (E / a) / (1 - n + n * n / 4 - 3 * n * n * n / 64 + 5 * n * n * n * n / 256);
    double cosB = (N / a) / (1 + n - 3 * n * n / 8 + 3 * n * n * n / 32 - 45 * n * n * n * n / 1024);
    double tanB = sinB / cosB;
    double L_L0 = atan(tanB / sqrt(1 - e12 * sinB * sinB / (1 - e12)));
    double L1 = L_L0 + L0;
    double M = a * ((1 - e12 / 4 - 3 * e12 * e12 / 64 - 5 * e12 * e12 * e12 / 256) * sinB
             - (3 * e12 / 8 + 3 * e12 * e12 / 32 + 45 * e12 * e12 * e12 / 1024) * sin(2 * sinB)
             + (15 * e12 * e12 / 256 + 45 * e12 * e12 * e12 / 1024) * sin(4 * sinB)
             - (35 * e12 * e12 * e12 / 3072) * sin(6 * sinB));
    double N1 = a / sqrt(1 - e12 * sinB * sinB);
    double t2 = tanB * tanB;
    double N2 = N1 / sqrt(1 + e22 * cosB * cosB);
    double B1 = sinB - N2 * tanB * (M - N * N / N2) / a;
    int count = 0;
    while (fabs(B1 - sinB) > 1e-10 && count < 10) {
        count++;
        sinB = B1;
        cosB = sqrt(1 - e12 * sinB * sinB / (1 - e12));
        tanB = sinB / cosB;
        M = a * ((1 - e12 / 4 - 3 * e12 * e12 / 64 - 5 * e12 * e12 * e12 / 256) * sinB
                - (3 * e12 / 8 + 3 * e12 * e12 / 32 + 45 * e12 * e12 * e12 / 1024) * sin(2 * sinB)
                + (15 * e12 * e12 / 256 + 45 * e12 * e12 * e12 / 1024) * sin(4 * sinB)
                - (35 * e12 * e12 * e12 / 3072) * sin(6 * sinB));
        N1 = a / sqrt(1 - e12 * sinB * sinB);
        t2 = tanB * tanB;
        N2 = N1 / sqrt(1 + e22 * cosB * cosB);
        B1 = sinB - N2 * tanB * (M - N * N / N2) / a;
    }
    B = B1 * 180 / PI;
    L = L1 * 180 / PI;
}

int main() {
    double B = 31.2383; // 经度
    double L = 121.5015; // 纬度
    double B0 = 31.0; // 经线中央经度
    double L0 = 120.0; // 纬线中央纬度
    double N, E;
    geo2plane(B, L, B0, L0, N, E);
    cout << "平面坐标系(N, E)：" << N << ", " << E << endl;
    double B1, L1;
    plane2geo(N, E, B0, L0, B1, L1);
    cout << "经纬度坐标系(B, L)：" << B1 << ", " << L1 << endl;
    return 0;
}

其中 `geo2plane` 函数用于将经纬度坐标系下的点 `(B, L)` 转换为平面坐标系下的点 `(N, E)`，`plane2geo` 函数用于将平面坐标系下的点 `(N, E)` 转换为经纬度坐标系下的点 `(B, L)`。函数中输入的参数单位为度，输出的参数单位为米（平面坐标系）或度（经纬度坐标系）。