高斯投影正反算C语言程序代码,高斯投影正反算c代码

以下是高斯投影正反算的C语言程序代码：

正算：

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

#define PI 3.14159265358979323846
#define a 6378245.0
#define f 1.0 / 298.3
#define b (1.0 - f) * a
#define ee (a*a - b*b) / (a*a)

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

void Gauss(double lat, double lon, double *x, double *y)
{
    int ProjNo = 0;
    int ZoneWide = 6; //6度带宽
    double longitude0, X0, Y0, xval, yval;
    double a0, a2, a4, a6, a8, a10;
    double yita2, M, N, t, p;
    double l, _l, C, i, ii, iii, III, IV, V, VI;

    lat = rad(lat);
    lon = rad(lon);

    if (lon >= 75 && lon < 135)
    {
        longitude0 = 120.0;
        ProjNo = (int)((lon - longitude0) / ZoneWide) + 1;
    }
    else if (lon >= 135 && lon < 195)
    {
        longitude0 = 150.0;
        ProjNo = (int)((lon - longitude0) / ZoneWide) + 31;
    }
    else if (lon >= 195 && lon < 255)
    {
        longitude0 = 210.0;
        ProjNo = (int)((lon - longitude0) / ZoneWide) + 61;
    }
    else if (lon >= 255 && lon < 315)
    {
        longitude0 = 240.0;
        ProjNo = (int)((lon - longitude0) / ZoneWide) + 91;
    }
    else if (lon >= 315 || lon < 75)
    {
        longitude0 = 0.0;
        ProjNo = (int)((lon - longitude0) / ZoneWide) + 121;
        if (ProjNo > 180)
        {
            ProjNo = ProjNo - 180;
        }
    }

    a0 = a * (1 - ee) * (1 - ee*ee/4 - 3*ee*ee*ee/64 - 5*ee*ee*ee*ee/256);
    a2 = 3.0 / 2.0 * ee - 27.0 / 32.0 * ee*ee*ee + 269.0 / 512.0 * ee*ee*ee*ee - 6607.0 / 245760.0 * ee*ee*ee*ee*ee;
    a4 = 21.0 / 16.0 * ee*ee - 55.0 / 32.0 * ee*ee*ee + 699.0 / 512.0 * ee*ee*ee*ee - 8183.0 / 204840.0 * ee*ee*ee*ee*ee;
    a6 = 151.0 / 96.0 * ee*ee*ee - 417.0 / 128.0 * ee*ee*ee*ee + 52121.0 / 155520.0 * ee*ee*ee*ee*ee;
    a8 = 1097.0 / 512.0 * ee*ee*ee*ee - 1117.0 / 384.0 * ee*ee*ee*ee*ee;
    a10 = 8011.0 / 2560.0 * ee*ee*ee*ee*ee;

    yita2 = ee * cos(lat) * cos(lat);
    M = a * ((1 - ee / 4 - 3 * ee * ee / 64 - 5 * ee * ee * ee / 256) * lat - (3 * ee / 8 + 3 * ee * ee / 32 + 45 * ee * ee * ee / 1024) * sin(2 * lat) + (15 * ee * ee / 256 + 45 * ee * ee * ee / 1024) * sin(4 * lat) - 35 * ee * ee * ee / 3072 * sin(6 * lat));
    N = a / sqrt(1 - ee * sin(lat) * sin(lat));

    t = tan(lat) * tan(lat);
    p = lon - longitude0;

    l = p * cos(lat);
    _l = p * cos(lat) * cos(lat) * cos(lat) / 6 * (1 - tan(lat) * tan(lat) + yita2);
    C = yita2 * yita2 * cos(lat) * cos(lat) * (1 - t) / 24;
    i = p * cos(lat) * cos(lat) * cos(lat) * cos(lat) * (5 - t + 9 * yita2 + 4 * yita2 * yita2) / 120;
    ii = p * cos(lat) * cos(lat) * cos(lat) * cos(lat) * cos(lat) * (61 - 58 * t + t * t + 600 * yita2 - 330 * ee * ee) / 720;
    iii = p * cos(lat) * cos(lat) * cos(lat) * cos(lat) * cos(lat) * cos(lat) * (1385 - 3111 * t + 543 * t * t - t * t * t) / 40320;

    X0 = 1000000L * ProjNo + 500000L;
    Y0 = 0;

    xval = X0 + M + N * l * l * (0.5 + l * l * (a2 + a4 + a6 + a8 + a10));
    yval = Y0 + N * (l + _l + C + i + ii + iii);

    *x = xval;
    *y = yval;
}

int main()
{
    double lat, lon, x, y;

    lat = 31.2304;
    lon = 121.4737;

    Gauss(lat, lon, &x, &y);

    printf("x: %lf, y: %lf\n", x, y);

    return 0;
}

反算：

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

#define PI 3.14159265358979323846
#define a 6378245.0
#define f 1.0 / 298.3
#define b (1.0 - f) * a
#define ee (a*a - b*b) / (a*a)

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

void Gauss_Reverse(double x, double y, double *lat, double *lon)
{
    int ProjNo = 0;
    int ZoneWide = 6; //6度带宽
    double longitude0, X0, Y0, xval, yval;
    double e1, e2, f1, f2, a1, a2, j1, j2, m;

    X0 = 1000000L;
    Y0 = 0;
    x -= X0;
    y -= Y0;

    if (y >= 0)
    {
        ProjNo = (int)(y / 1000000L);
        longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
        yval = y - ProjNo * 1000000L - 500000L;
    }
    else
    {
        ProjNo = (int)((y + 1000000L) / 1000000L);
        longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
        yval = y + 1000000L - ProjNo * 1000000L - 500000L;
    }

    xval = x / (a * 1.0);

    e1 = (1 - sqrt(1 - ee)) / (1 + sqrt(1 - ee));
    e2 = e1 * e1;
    f1 = 1 - sqrt(1 - ee);
    f2 = f1 * f1;
    a1 = 1 + 3 * e1 / 2 + 45 * e1 * e1 / 16 + 175 * e1 * e1 * e1 / 64 + 11025 * e1 * e1 * e1 * e1 / 256;
    a2 = 1 - e1 - 3 * e1 * e1 / 4 - 45 * e1 * e1 * e1 / 64 - 175 * e1 * e1 * e1 * e1 / 256;
    j1 = 3 * e1 / 2 - 27 * e1 * e1 * e1 / 32 + 269 * e1 * e1 * e1 * e1 / 512;
    j2 = 21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 / 32;
    m = b * a1 * (lat - j1 * sin(2 * lat) + j2 * sin(4 * lat) - 5 * e1 * sin(6 * lat)) / a;

    *lat = (m + (1 + 2 * t + c * c) * sin(lat) * cos(lat) / 2 + (5 - 2 * c * c + 28 * t - 3 * c * c * c * c + 8 * e1 * e1 + 24 * t * t) * sin(2 * lat) * cos(2 * lat) / 24 + (61 - 90 * t + 45 * t * t + 46 * e1 * e1 - 252 * t * e1 * e1 - 3 * c * c) * sin(4 * lat) * cos(4 * lat) / 720) / (1 - e2);

    *lon = xval / cos(*lat) / a2 + longitude0;
    *lat = *lat * 180 / PI;
    *lon = *lon * 180 / PI;
}

int main()
{
    double x, y, lat, lon;

    x = 344931.80;
    y = 3453545.02;

    Gauss_Reverse(x, y, &lat, &lon);

    printf("lat: %lf, lon: %lf\n", lat, lon);

    return 0;
}