高斯投影正反算程序C++

高斯投影正反算是地图制图中常用的一种方法，下面是C++实现的代码：

#include <iostream>
#include <cmath>
using namespace std;

const double PI = 3.14159265358979323846;
const double a = 6378245.0; // 长半轴
const double b = 6356863.0188; // 短半轴
const double e2 = (a*a - b*b) / (a*a); // 第一偏心率的平方
const double e12 = (a*a - b*b) / (b*b); // 第二偏心率的平方
const double k0 = 1.0; // 中央经线的缩放因子
const double X0 = 500000.0; // 中央经线的X坐标偏移量
const double Y0 = 0.0; // 中央经线的Y坐标偏移量

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

double deg(double x) {
    return x * 180.0 / PI;
}

// 计算子午线弧长
double S(double B) {
    double sinB = sin(B);
    double sin2B = sinB * sinB;
    double sin4B = sin2B * sin2B;
    double sin6B = sin4B * sin2B;
    double sin8B = sin4B * sin4B;
    double A0 = a * (1 - e2);
    double A2 = 1.25 * e2 - 1.71875 * e12;
    double A4 = 1.3125 * e2 - 1.71875 * e12 * e2;
    double A6 = 1.57291666666667 * e12 - 1.59722222222222 * e12 * e2;
    double A8 = 1.45833333333333 * e12 * e2;
    double S = A0 * (B - sinB * cos(B) * (A2 + sin2B * (A4 + sin2B * (A6 + sin2B * A8))));
    return S;
}

// 计算高斯投影正算
void GaussForward(double B, double L, double &X, double &Y) {
    double L0 = rad(117); // 中央经线
    double L1 = rad(L); // 经度转换为弧度
    double B1 = rad(B); // 纬度转换为弧度
    double t = tan(B1);
    double eta2 = e12 * cos(B1) * cos(B1);
    double N = a / sqrt(1 - e2 * sin(B1) * sin(B1));
    double M = a * (1 - e2) / pow(1 - e2 * sin(B1) * sin(B1), 1.5);
    double l = L1 - L0;
    double l2 = l * l;
    double l4 = l2 * l2;
    double l6 = l4 * l2;
    double l8 = l6 * l2;
    X = X0 + k0 * N * (l + l3 / 6.0 * (1 - t*t + eta2) + l5 / 120.0 * (5 - 18 * t*t + t*t*t*t + 14 * eta2 - 58 * t*t*eta2));
    Y = Y0 + k0 * (S(B1) + N * tan(B1) * (l2 / 2.0 + l4 / 24.0 * (5 + 3 * t*t + eta2 - 9 * t*t*eta2) + l6 / 720.0 * (61 - 58 * t*t + t*t*t*t + 270 * eta2 - 330 * t*t*eta2) + l8 / 40320.0 * (1385 - 3111 * t*t + 543 * t*t*t*t - t*t*t*t*t*t))));
}

// 计算高斯投影反算
void GaussInverse(double X, double Y, double &B, double &L) {
    double L0 = rad(117); // 中央经线
    double X1 = X - X0;
    double Y1 = Y - Y0;
    double M0 = a * (1 - e2);
    double M2 = 3.0 / 2.0 * e2 * M0;
    double M4 = 5.0 / 4.0 * e2 * M2;
    double M6 = 7.0 / 6.0 * e2 * M4;
    double M8 = 9.0 / 8.0 * e2 * M6;
    double Bf = Y1 / (k0 * M0);
    double Bf1 = Bf;
    double Bf2 = Bf;
    while (true) {
        Bf1 = Bf2;
        Bf2 = (Y1 + M2 * sin(2 * Bf1) - M4 * sin(4 * Bf1) + M6 * sin(6 * Bf1) - M8 * sin(8 * Bf1)) / (k0 * M0);
        if (fabs(Bf2 - Bf1) < 1e-10) {
            break;
        }
    }
    double Nf = a / sqrt(1 - e2 * sin(Bf2) * sin(Bf2));
    double tf = tan(Bf2);
    double tf2 = tf * tf;
    double tf4 = tf2 * tf2;
    double tf6 = tf4 * tf2;
    double tf8 = tf4 * tf4;
    double eta2 = e12 * cos(Bf2) * cos(Bf2);
    double l = X1 / (k0 * Nf);
    double l3 = l * l * l;
    double l5 = l3 * l * l;
    double l7 = l5 * l * l;
    double l9 = l7 * l * l;
    B = Bf2 - tf2 / 2.0 * Nf / M0 * l*l + tf4 / 24.0 * Nf / M0 * (5 + 3 * tf2 + eta2 - 9 * tf2*eta2) * l3 - tf6 / 720.0 * Nf / M0 * (61 + 90 * tf2 + 45 * tf4) * l5 + tf8 / 40320.0 * Nf / M0 * (1385 + 3633 * tf2 + 4095 * tf4 + 1575 * tf6) * l7;
    L = L0 + 1.0 / cos(Bf2) * (l - 1.0 / 6.0 * cos(Bf2) * sin(Bf2) * (1 + tf2 + eta2) * l3 + 1.0 / 120.0 * cos(Bf2) * sin(Bf2) * sin(Bf2) * (5 - 18 * tf2 + tf4 + 14 * eta2 - 58 * tf2*eta2) * l5 - 1.0 / 5040.0 * cos(Bf2) * sin(Bf2) * sin(Bf2) * sin(Bf2) * (61 - 479 * tf2 + 179 * tf4 - tf6) * l9);
    B = deg(B);
    L = deg(L);
}

int main() {
    double B = 39.9075; // 纬度
    double L = 116.39723; // 经度
    double X, Y;
    GaussForward(B, L, X, Y);
    cout << "高斯投影正算结果：" << endl;
    cout << "X = " << X << endl;
    cout << "Y = " << Y << endl;
    GaussInverse(X, Y, B, L);
    cout << "高斯投影反算结果：" << endl;
    cout << "B = " << B << endl;
    cout << "L = " << L << endl;
    return 0;
}