正轴等角双标准纬线圆锥投影算法C++代码

以下是一个简单的 C++实现，实现了正轴等角双标准纬线圆锥投影算法：

#include <iostream>
#include <cmath>

using namespace std;

const double PI = 3.14159265358979323846;
const double R = 6371000; // 地球半径

// 经纬度转弧度
double deg2rad(double deg) {
    return deg * PI / 180;
}

// 弧度转经纬度
double rad2deg(double rad) {
    return rad * 180 / PI;
}

// 计算子午线弧长
double meridian_arc(double lat) {
    double a = 6378137; // 赤道半径
    double b = 6356752.3142; // 极半径
    double e = sqrt(1 - pow(b / a, 2));
    double e2 = pow(e, 2);
    double m1 = a * (1 - e2);
    double m2 = a * (1 - e2 / 4);
    double m3 = a * (1 - 3 * e2 / 4);
    double m4 = a * (1 - e2 * 0.75);
    double m5 = a * (1 - e2 * 15 / 16);
    double m6 = a * (1 - e2 * 35 / 36);
    double m7 = a * (1 - e2 * 315 / 324);
    double m8 = a * (1 - e2 * 693 / 784);
    double m9 = a * (1 - e2 * 639 / 680);
    double m10 = a * (1 - e2 * 175 / 192);
    double m11 = a * (1 - e2 * 229 / 256);
    double m12 = a * (1 - e2 * 100 / 81);
    double m13 = a * (1 - e2 * 361 / 400);
    double m14 = a * (1 - e2 * 99 / 100);
    double m15 = a * (1 - e2 * 51 / 64);
    double m16 = a * (1 - e2 * 533 / 576);
    double m17 = a * (1 - e2 * 49561 / 51840);
    double m18 = a * (1 - e2 * 175 / 168);
    double m19 = a * (1 - e2 * 81 / 80);
    double m20 = a * (1 - e2 * 1397 / 1280);
    double m21 = a * (1 - e2 * 601 / 560);
    double m22 = a * (1 - e2 * 34729 / 30240);
    double m23 = a * (1 - e2 * 193 / 180);
    double m24 = a * (1 - e2 * 779 / 720);
    double m25 = a * (1 - e2 * 4397 / 4032);
    double m26 = a * (1 - e2 * 4583 / 4320);
    double m27 = a * (1 - e2 * 126217 / 115200);
    double m28 = a * (1 - e2 * 63 / 64);
    double m29 = a * (1 - e2 * 364463 / 349920);
    double m30 = a * (1 - e2 * 4351 / 4032);
    double m31 = a * (1 - e2 * 79903 / 74880);
    double m32 = a * (1 - e2 * 9186421 / 8775360);

    double phi = deg2rad(lat);
    double m = m1 * phi
             - m2 * sin(2 * phi) / 2
             + m3 * sin(4 * phi) / 4
             - m4 * sin(6 * phi) / 6
             + m5 * sin(8 * phi) / 8
             - m6 * sin(10 * phi) / 10
             + m7 * sin(12 * phi) / 12
             - m8 * sin(14 * phi) / 14
             + m9 * sin(16 * phi) / 16
             - m10 * sin(18 * phi) / 18
             + m11 * sin(20 * phi) / 20
             - m12 * sin(22 * phi) / 22
             + m13 * sin(24 * phi) / 24
             - m14 * sin(26 * phi) / 26
             + m15 * sin(28 * phi) / 28
             - m16 * sin(30 * phi) / 30
             + m17 * sin(32 * phi) / 32
             - m18 * sin(34 * phi) / 34
             + m19 * sin(36 * phi) / 36
             - m20 * sin(38 * phi) / 38
             + m21 * sin(40 * phi) / 40
             - m22 * sin(42 * phi) / 42
             + m23 * sin(44 * phi) / 44
             - m24 * sin(46 * phi) / 46
             + m25 * sin(48 * phi) / 48
             - m26 * sin(50 * phi) / 50
             + m27 * sin(52 * phi) / 52
             - m28 * sin(54 * phi) / 54
             + m29 * sin(56 * phi) / 56
             - m30 * sin(58 * phi) / 58
             + m31 * sin(60 * phi) / 60
             - m32 * sin(62 * phi) / 62;
    return m;
}

// 计算投影后的坐标
void project(double lat, double lon, double std_lat, double r_lat, double &x, double &y) {
    double m0 = meridian_arc(std_lat); // 第一标准纬线的子午线弧长
    double m1 = meridian_arc(r_lat); // 第二标准纬线的子午线弧长
    double n = (log(cos(std_lat)) - log(cos(r_lat))) / (log(tan(PI / 4 + r_lat / 2)) - log(tan(PI / 4 + std_lat / 2))); // 投影系数
    double F = cos(std_lat) * pow(tan(PI / 4 + r_lat / 2), n) / n; // F 值
    double p = R * F / pow(tan(PI / 4 + deg2rad(lat) / 2), n); // 投影半径
    double theta = n * (deg2rad(lon) - deg2rad(r_lat)); // 投影角度
    x = p * sin(theta);
    y = R * F - p * cos(theta);
}

int main() {
    double lat = 31.2304; // 纬度
    double lon = 121.4737; // 经度
    double std_lat = 0; // 第一标准纬线
    double r_lat = 30; // 第二标准纬线
    double x, y;
    project(lat, lon, std_lat, r_lat, x, y);
    cout << "x: " << x << ", y: " << y << endl;
    return 0;
}

在此代码中，我们使用了地球半径 R = 6371000 米，以及第一标准纬线 std_lat = 0 度和第二标准纬线 r_lat = 30 度。在实际应用中，这些值可以根据需要进行修改。