c语言高斯投影正反算

高斯投影正反算是指将经纬度坐标转换为高斯平面坐标，或将高斯平面坐标转换为经纬度坐标的过程。下面是C语言实现高斯投影正反算的示例代码：

1.高斯投影正算

#include <stdio.h>
#include <math.h>

#define PI 3.14159265358979323846
#define a 6378137.0
#define b 6356752.3142
#define f (a-b)/a
#define e sqrt(2*f-f*f)
#define ee e*e
#define M_PI 3.14159265358979323846

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

void Gauss(double B, double L, double &X, double &Y)
{
    double L0 = rad(120);
    double l = L - L0;
    double t = tan(B);
    double n2 = ee * cos(B) * cos(B);
    double l2 = l * l;
    double l4 = l2 * l2;
    double l6 = l4 * l2;
    double N = a / sqrt(1 - ee * sin(B) * sin(B));
    X = 500000 + (1 + l2 / 2 + (5 - t * t + 9 * n2 + 4 * n2 * n2) * l4 / 24 + (61 - 58 * t * t + t * t * t * t) * l6 / 720) * N * cos(B) * l;    Y = (1 + (1 - t * t + n2) * l2 / 6 + (5 - 18 * t * t + t * t * t * t + 14 * n2 - 58 * n2 * t * t) * l4 / 120) * N * cos(B) * l * l * l / 6 + N * rad(B);
}

int main()
{
    double B = rad(31.177427);
    double L = rad(121.527210);
    double X, Y;
    Gauss(B, L, X, Y);
    printf("X: %lf, Y: %lf\n", X, Y);
    return 0;
}

2.高斯投影反算

#include <stdio.h>
#include <math.h>

#define PI 3.14159265358979323846
#define a 6378137.0
#define b 6356752.3142
#define f (a-b)/a
#define e sqrt(2*f-f*f)
#define ee e*e
#define M_PI 3.14159265358979323846

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

void Gauss_Reverse(double X, double Y, double &B, double &L)
{
    double L0 = rad(120);
    double x = X - 500000;
    double y = Y;
    double e1 = (1 - sqrt(1 - ee)) / (1 + sqrt(1 - ee));
    double M = y;
    double mu = M / (a * (1 - ee / 4 - 3 * ee * ee / 64 - 5 * ee * ee * ee / 256));
    double e2 = (e * e) / (1 - e * e);
    double C1 = e2 * cos(mu) * cos(mu);
    double T1 = tan(mu) * tan(mu);
    double N1 = a / sqrt(1 - ee * sin(mu) * sin(mu));
    double R1 = a * (1 - ee) / pow(1 - ee * sin(mu) * sin(mu), 1.5);
    double D = x / N1;
    B = mu - (N1 * tan(mu) / R1) * (D * D / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * e2) * pow(D, 4) / 24 + (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * e2 - 3 * C1 * C1) * pow(D, 6) / 720);
    L = L0 + (D - (1 + 2 * T1 + C1) * pow(D, 3) / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * e2 + 24 * T1 * T1) * pow(D, 5) / 120) / cos(mu);
    B = B * 180 / PI;
    L = L * 180 / PI;
}

int main()
{
    double X = 345704.0;
    double Y = 345704.0;
    double B, L;
    Gauss_Reverse(X, Y, B, L);
    printf("B: %lf, L: %lf\n", B, L);
    return 0;
}