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

下面是一个简单的 C++ 程序，用于将北东坐标系转换为经纬度坐标系：

#include <iostream>
#include <cmath>

using namespace std;

const double pi = 3.14159265358979323846;
const double a = 6378137;
const double b = 6356752.3142;
const double e2 = 1 - pow(b / a, 2);

double rad(double d) {
    return d * pi / 180.0;
}

void NEtoLL(double ne_latitude, double ne_longitude, double &ll_latitude, double &ll_longitude) {
    double dx = 25000000;
    double dy = 0;
    double dz = 0;
    double L0 = rad(117);
    double sinBf = ne_latitude / (sqrt(pow(ne_latitude, 2) + pow(ne_longitude, 2)));
    double cosBf = ne_longitude / (sqrt(pow(ne_latitude, 2) + pow(ne_longitude, 2)));
    double Nf = a / sqrt(1 - e2 * pow(sinBf, 2));
    double Mf = a * (1 - e2) / pow(sqrt(1 - e2 * pow(sinBf, 2)), 3);
    double tf = tan(rad(1) * (ne_latitude / 3600 + 45));
    double Nf_2 = Nf * pow(cos(rad(ne_latitude / 3600)), 2);
    double x = dx + (Nf + dz) * cos(rad(ne_latitude / 3600)) * cos(rad(ne_longitude / 3600) - L0);
    double y = dy + (Nf + dz) * cos(rad(ne_latitude / 3600)) * sin(rad(ne_longitude / 3600) - L0);
    double z = dz + Nf_2 / (1 + e2 * pow(cos(rad(ne_latitude / 3600)), 2)) * tf * sin(rad(ne_latitude / 3600));
    double p = sqrt(pow(x, 2) + pow(y, 2));
    double theta = atan(z * a / (p * b));
    double sintheta = sin(theta);
    double costheta = cos(theta);
    ll_latitude = atan((z + pow(e2 * b * sintheta, 3)) / (p - e2 * a * costheta * costheta * costheta));
    ll_longitude = L0 + atan(y / x);
}

该程序使用的是WGS84椭球体参数，可根据需要修改。使用方法如下：

int main() {
    double ne_latitude = 4151811.16;
    double ne_longitude = 1320241.12;
    double ll_latitude, ll_longitude;
    NEtoLL(ne_latitude, ne_longitude, ll_latitude, ll_longitude);
    cout << "经度：" << ll_longitude * 180 / pi << endl;
    cout << "纬度：" << ll_latitude * 180 / pi << endl;
    return 0;
}

其中，`ne_latitude` 和 `ne_longitude` 分别代表北东坐标系下的纬度和经度，`ll_latitude` 和 `ll_longitude` 分别代表经纬度坐标系下的纬度和经度。运行程序后，会输出经度和纬度的值。