高斯投影正反算c#

以下是C#中实现高斯投影正反算的示例代码，供参考：

using System;

public static class GaussProj
{
    private const double a = 6378137.0; // 长半轴
    private const double b = 6356752.3142; // 短半轴
    private const double f = 0.00335281066474748; // 扁率
    private const double e = 0.0818191908426215; // 第一偏心率
    private const double e2 = 0.00669437999019758; // 第二偏心率
    private const double m0 = 1.0; // 中央经线的缩比因子

    // 弧度转角度
    private static double RadToDeg(double rad)
    {
        return rad / Math.PI * 180.0;
    }

    // 角度转弧度
    private static double DegToRad(double deg)
    {
        return deg / 180.0 * Math.PI;
    }

    // 计算子午线弧长
    private static double MeridianDist(double B)
    {
        double A = a;
        double C = Math.Pow(a, 2) / b;
        double alpha = (A - b) / (A + b);
        double beta = (3.0 / 2.0) * alpha - (27.0 / 32.0) * Math.Pow(alpha, 3);
        double gamma = (21.0 / 16.0) * Math.Pow(alpha, 2) - (55.0 / 32.0) * Math.Pow(alpha, 4);
        double delta = (151.0 / 96.0) * Math.Pow(alpha, 3);
        double epsilon = (1097.0 / 512.0) * Math.Pow(alpha, 4);
        double B_ = B - DegToRad(1.0 / 60.0);
        double C1 = beta * Math.Sin(2 * B_) * Math.Cosh(2 * Math.Atan(Math.Sqrt((1 + e) * Math.Tan(B_)) - e * Math.Sqrt(1 - Math.Pow(e * Math.Sin(B_), 2))));
        double C2 = gamma * Math.Sin(4 * B_) * Math.Cosh(4 * Math.Atan(Math.Sqrt((1 + e) * Math.Tan(B_)) - e * Math.Sqrt(1 - Math.Pow(e * Math.Sin(B_), 2))));
        double C3 = delta * Math.Sin(6 * B_) * Math.Cosh(6 * Math.Atan(Math.Sqrt((1 + e) * Math.Tan(B_)) - e * Math.Sqrt(1 - Math.Pow(e * Math.Sin(B_), 2))));
        double C4 = epsilon * Math.Sin(8 * B_) * Math.Cosh(8 * Math.Atan(Math.Sqrt((1 + e) * Math.Tan(B_)) - e * Math.Sqrt(1 - Math.Pow(e * Math.Sin(B_), 2))));
        double M = C * (alpha * (B + C1 + C2 + C3 + C4) - Math.Sin(B));
        return M;
    }

    // 计算投影坐标
    public static (double X, double Y) GaussProjForward(double B, double L, double L0)
    {
        L = L - L0;
        double l = L;
        double A = a;
        double T = Math.Pow(Math.Tan(B), 2);
        double C = Math.Pow(e * Math.Cos(B), 2);
        double L2 = Math.Pow(L, 2);
        double N = A / Math.Sqrt(1.0 - (e2 * Math.Pow(Math.Sin(B), 2)));
        double M = MeridianDist(B);
        double X = m0 * N * (l + (1.0 / 2.0) * Math.Pow(l, 2) * N * Math.Cos(B) * Math.Sin(B) +
            (1.0 / 24.0) * Math.Pow(l, 4) * N * Math.Cos(B) * Math.Pow(Math.Sin(B), 3) *
            (5 - T + 9 * C + 4 * Math.Pow(C, 2)) +
            (1.0 / 720.0) * Math.Pow(l, 6) * N * Math.Cos(B) * Math.Pow(Math.Sin(B), 5) *
            (-61 + 58 * T + Math.Pow(T, 2) + 600 * C - 330 * f * f)) +
            500000.0;
        double Y = m0 * (M + N * Math.Tan(B) * (L2 / 2 + (1.0 / 6.0) * Math.Pow(L2, 2) * Math.Pow(Math.Cos(B), 2) *
            (1 + 2 * T + C) + (1.0 / 120.0) * Math.Pow(L2, 4) * Math.Pow(Math.Cos(B), 4) *
            (-1 + 2 * Math.Pow(T, 2) + Math.Pow(C, 2) - 2 * C - 7 * T))) +
            0.0;
        return (X, Y);
    }

    // 计算经纬度
    public static (double B, double L) GaussProjInverse(double X, double Y, double L0)
    {
        X = X - 500000.0;
        double A = a;
        double Mf = MeridianDist(0.0);
        double M = Y / m0 + Mf;
        double mu = M / (A * (1.0 - Math.Pow(e, 2) / 4.0 - 3 * Math.Pow(e, 4) / 64.0 - 5 * Math.Pow(e, 6) / 256.0));
        double e1 = (1.0 - Math.Sqrt(1.0 - Math.Pow(e, 2))) / (1.0 + Math.Sqrt(1.0 - Math.Pow(e, 2)));
        double J1 = (3 * e1 / 2 - 27 * Math.Pow(e1, 3) / 32.0);
        double J2 = (21 * Math.Pow(e1, 2) / 16 - 55 * Math.Pow(e1, 4) / 32.0);
        double J3 = (151 * Math.Pow(e1, 3) / 96.0);
        double J4 = (1097 * Math.Pow(e1, 4) / 512.0);
        double fp = mu + J1 * Math.Sin(2 * mu) + J2 * Math.Sin(4 * mu) + J3 * Math.Sin(6 * mu) + J4 * Math.Sin(8 * mu);
        double C1 = e2 * Math.Pow(Math.Cos(fp), 2);
        double T1 = Math.Pow(Math.Tan(fp), 2);
        double R1 = a * (1 - Math.Pow(e, 2)) / Math.Pow(1 - Math.Pow(e * Math.Sin(fp), 2), 1.5);
        double N1 = a / Math.Sqrt(1 - Math.Pow(e * Math.Sin(fp), 2));
        double D = X / (N1 * m0);
        double Q1 = N1 * Math.Tan(fp) / R1;
        double Q2 = D * D / 2;
        double Q3 = (5 + 3 * T1 + 10 * C1 - 4 * Math.Pow(C1, 2) - 9 * e2) * Math.Pow(D, 4) / 24;
        double Q4 = (61 + 90 * T1 + 298 * C1 + 45 * Math.Pow(T1, 2) - 252 * e2 - 3 * Math.Pow(C1, 2)) * Math.Pow(D, 6) / 720;
        double Bf = fp - Q1 * (Q2 - Q3 + Q4);
        double Lf = L0 + D - (1 + 2 * T1 + C1) * Math.Pow(D, 3) / 6 + (5 - 2 * C1 + 28 * T1 - 3 * Math.Pow(C1, 2) +
            8 * e2 + 24 * Math.Pow(T1, 2)) * Math.Pow(D, 5) / 120;
        return (RadToDeg(Bf), RadToDeg(Lf));
    }
}

使用示例：

// 正算示例
var x = 500000.0;
var y = 4000000.0;
var l0 = 117.0;
var (B, L) = GaussProj.GaussProjInverse(x, y, l0);
Console.WriteLine($"({x}, {y}) -> ({B}, {L})");

// 反算示例
var B = 31.0;
var L = 118.0;
var l0 = 117.0;
var (x, y) = GaussProj.GaussProjForward(DegToRad(B), DegToRad(L), DegToRad(l0));
Console.WriteLine($"({B}, {L}) -> ({x}, {y})");

注意事项：

- 本示例代码仅适用于中国国家标准的高斯投影坐标系，其他坐标系需要修改相关参数。
- 代码中涉及的参数值为标准值，实际使用中可能需要根据具体情况进行调整。