经纬度坐标系转北东地坐标系c++程序

以下是C++程序，可以实现经纬度坐标系转北东地坐标系：

#include <iostream>
#include <cmath>

const double PI = 3.14159265358979323846;
const double a = 6378137.0; // 地球长半轴
const double b = 6356752.3142; // 地球短半轴
const double f = 1 / 298.257223563; // 地球扁率

double toRadians(double degree) {
    return degree * PI / 180.0;
}

double toDegrees(double radians) {
    return radians * 180.0 / PI;
}

double calculateM(double latitude) {
    double e2 = 2 * f - f * f;
    return a * (1 - e2) / pow(1 - e2 * sin(toRadians(latitude)) * sin(toRadians(latitude)), 1.5);
}

double calculateN(double latitude) {
    double e2 = 2 * f - f * f;
    return a / sqrt(1 - e2 * sin(toRadians(latitude)) * sin(toRadians(latitude)));
}

double calculateX(double longitude, double latitude) {
    double longitude0 = toRadians(120.0);
    double m0 = calculateM(31.0);
    double n0 = calculateN(31.0);
    double m = calculateM(latitude);
    double n = calculateN(latitude);
    double t = tan(toRadians(latitude));
    double l = toRadians(longitude) - longitude0;
    double cos_latitude = cos(toRadians(latitude));
    double sin_latitude = sin(toRadians(latitude));
    double x = m * l * cos_latitude + m * pow(l, 3) * pow(cos_latitude, 3) * (1 / 6.0 + pow(t, 2) + pow(l, 2) / 20.0) + m * pow(l, 5) * pow(cos_latitude, 5) * (1 / 120.0 - pow(t, 2) / 24.0) + m * pow(l, 7) * pow(cos_latitude, 7) * (1 / 5040.0 + pow(t, 2) / 240.0);
    x += 500000.0;
    return x;
}

double calculateY(double longitude, double latitude) {
    double longitude0 = toRadians(120.0);
    double m0 = calculateM(31.0);
    double n0 = calculateN(31.0);
    double m = calculateM(latitude);
    double n = calculateN(latitude);
    double t = tan(toRadians(latitude));
    double l = toRadians(longitude) - longitude0;
    double cos_latitude = cos(toRadians(latitude));
    double sin_latitude = sin(toRadians(latitude));
    double y = n * (l * sin_latitude + pow(l, 3) * pow(cos_latitude, 2) * t / 3.0 + pow(l, 5) * pow(cos_latitude, 4) * (2 * pow(t, 2) + pow(l, 2)) / 15.0 + pow(l, 7) * pow(cos_latitude, 6) * (17 * pow(t, 3) + 16 * pow(t, 5)) / 315.0);
    double y0 = n0 * (m / n + sin_latitude * t * (pow(l, 2) / 2.0 + pow(t, 2) * pow(l, 4) / 24.0 + 5 * pow(t, 4) * pow(l, 6) / 624.0));
    y += y0;
    return y;
}

int main() {
    double longitude = 121.5; // 经度
    double latitude = 31.0; // 纬度
    double x = calculateX(longitude, latitude);
    double y = calculateY(longitude, latitude);
    std::cout << "x: " << x << std::endl;
    std::cout << "y: " << y << std::endl;
    return 0;
}

此程序以120度为中央经线，31度为基准纬线计算，也可以根据实际情况进行更改。