请写出高斯投影正反算代码，用c#语言

高斯投影正算代码：//高斯投影正算
//输入：经纬度坐标（B,L）
//输出：大地坐标（X,Y）
#include<math.h>
#define a 6378137.0 //椭球长半轴
#define b 6356752.3142 //椭球短半轴
#define pi 3.1415926
double B,L; //经纬度
double X,Y; //大地坐标

//以下为正算部分
double ee,N,T,C,A,M;
//椭球第一偏心率的平方
ee=pow(a*a-b*b,0.5)/a;
//卯酉圈曲率半径
N=a/pow(1-ee*ee*sin(B)*sin(B),0.5);
T=tan(B)*tan(B);
C=ee*ee*cos(B)*cos(B);
A=(L*pi)/180.0;
M=a*((1-ee*ee/4-3*ee*ee*ee*ee/64-5*ee*ee*ee*ee*ee*ee/256)*B-(3*ee*ee/8+3*ee*ee*ee*ee/32+45*ee*ee*ee*ee*ee*ee/1024)*sin(2*B)+(15*ee*ee*ee*ee/256+45*ee*ee*ee*ee*ee*ee/1024)*sin(4*B)-(35*ee*ee*ee*ee*ee*ee/3072)*sin(6*B));
X=N*cos(B)*A+M;
Y=N*sin(B)*A;

相关问题

高斯投影正反算代码 c++代码

高斯投影是一种常用的地理坐标系统转换方法，通过投影方式将地球上的经纬度坐标转换为平面直角坐标系。其正反算过程可以通过C代码实现。

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

#include <math.h>

// 定义一些常数
#define PI      3.1415926535897932384626433832795
#define L0      108     // 中央经线
#define A       6378137.0   // 长半轴
#define E2      0.006694380023
#define K0      1.0     // 比例因子
#define N       0.001679220394
#define a0      32140.404
#define a2      -135.3302
#define a4      0.7092
#define a6      0.0040
#define X_ORIGIN    500000.0
#define Y_ORIGIN    0.0

// 高斯正算函数
void GaussForward(double B, double L, double *X, double *Y) {
  double L_rad = L * PI / 180.0;
  double B_rad = B * PI / 180.0;
  double L0_rad = L0 * PI / 180.0;

  double cosB = cos(B_rad);
  double cos2B = cosB * cosB;
  double cos4B = cos2B * cos2B;

  double T = tan(B_rad);
  double T2 = T * T;
  double T4 = T2 * T2;
  double T6 = T2 * T4;

  double N_ = A / sqrt(1 - E2 * sin(B_rad) * sin(B_rad));
  double n = N_ / A;

  double deltaL = L_rad - L0_rad;

  double X_ = a0 * B_rad + a2 * sin(2.0 * B_rad) + a4 * sin(4.0 * B_rad) + a6 * sin(6.0 * B_rad);
  double X1 = K0 * N_ * deltaL * cosB;
  double X2 = K0 * N_ * deltaL * deltaL * deltaL * cosB * cosB * cosB / 6.0 * (1.0 - T2 + n + 5.0 * T4 * (1.0 - 18.0 * T2 + T4));
  double X3 = K0 * N_ * deltaL * deltaL * deltaL * deltaL * deltaL * cosB * cosB * cosB * cosB * cosB / 120.0 * (5.0 - 18.0 * T2 + T4 + 14.0 * n - 58.0 * T2 * n + 13.0 * n * n + 4.0 * n * n * n - 64.0 * T2 * n * n - 24.0 * T2 * T2);
  *X = X_ORIGIN + X_ + X1 + X2 + X3;

  double Y_ = N_ * tan(B_rad) * cos(deltaL) / 2.0;
  double Y1 = N_ * tan(B_rad) * cos(deltaL) * cos(deltaL) * cos(deltaL) / 24.0 * (5.0 - T2 + 9.0 * n + 4.0 * n * n);
  double Y2 = N_ * tan(B_rad) * cos(deltaL) * cos(deltaL) * cos(deltaL) * cos(deltaL) * cos(deltaL) / 720.0 * (61.0 - 58.0 * T2 + T4 + 270.0 * n - 330.0 * T2 * n + 445.0 * n * n + 324.0 * n * n * n - 680.0 * T2 * n * n + 88.0 * n * n * n * n - 600.0 * T2 * n * n * n - 192.0 * T2 * T2 - 48.0 * T2 * n * n * n * n);
  *Y = Y_ORIGIN + Y_ + Y1 + Y2;
}

// 测试函数
int main() {
  double B = 39.912345;
  double L = 116.56789;
  double X = 0.0;
  double Y = 0.0;

  GaussForward(B, L, &X, &Y);
  printf("经度：%f，纬度：%f 转换后的X坐标：%f，Y坐标：%f\n", B, L, X, Y);

  return 0;
}

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

#include <math.h>

// 定义一些常数
#define PI      3.1415926535897932384626433832795
#define L0      108     // 中央经线
#define A       6378137.0   // 长半轴
#define E2      0.006694380023
#define K0      1.0     // 比例因子
#define N       0.001679220394
#define a0      32140.404
#define a2      -135.3302
#define a4      0.7092
#define a6      0.0040
#define X_ORIGIN    500000.0
#define Y_ORIGIN    0.0

// 高斯反算函数
void GaussBackward(double X, double Y, double *B, double *L) {
  double X_ = X - X_ORIGIN;
  double Y_ = Y - Y_ORIGIN;

  double Bf = a0 * X_ + a2 * sin(2.0 * X_) + a4 * sin(4.0 * X_) + a6 * sin(6.0 * X_);
  double Ef = Bf - Y_ * Y_ * Y_ / (2.0 * 6378137.0 * 6378137.0 * N);
  double Tf = tan(Ef);
  double Nf = A / sqrt(1 - E2 * sin(Ef) * sin(Ef));
  double Q = X_ / Nf;

  double B_rad = Ef - Nf * tan(Ef) * tan(Ef) * (Q * Q) / (2.0 * 6378137.0 * 6378137.0);
  double Bf_rad = Bf * PI / 180.0;
  double L_rad = L0 * PI / 180.0 + Q - tan(Ef) * tan(Ef) * (Q * Q * Q) / (6.0 * 6378137.0 * 6378137.0 * Nf);
  double Lf_rad = L_rad * 180.0 / PI;

  *B = B_rad * 180.0 / PI;
  *L = Lf_rad;
}

// 测试函数
int main() {
  double X = 405000.123;
  double Y = 3456789.456;
  double B = 0.0;
  double L = 0.0;

  GaussBackward(X, Y, &B, &L);
  printf("X坐标：%f，Y坐标：%f 转换后的经度：%f，纬度：%f\n", X, Y, B, L);

  return 0;
}

这是一个简单的高斯投影正反算的C代码示例。使用时，可以将经纬度坐标传入GaussForward函数进行正算，或将平面直角坐标传入GaussBackward函数进行反算，即可得到相应的转换结果。请注意，具体的高斯投影参数值可能因地区而异，需要根据实际情况进行调整。

高斯投影正反算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})");

注意事项：

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