大地测量学基础算法实现

由于没有时间，目前只实现了大地主题解算、高斯投影坐标正反算和距离改正，但应该很容易扩充。

结构示意图：

参考椭球类：

  1 //参考椭球类，可选择克拉索夫斯基椭球体、1975年国际椭球体、WGS-84椭球体，也可以自定义椭球参数。
  2     class ReferenceEllipsoid
  3     {
  4         #region 椭球参数
  5         //克拉索夫斯基椭球体参数
  6         private const double a_1=6378245.0000000000;
  7         private const double b_1=6356863.0187730473;
  8         private const double c_1=6399698.9017827110;
  9         private const double r_alpha_1=298.3;
 10         private const double e2_1=0.006693421622966;
 11         private const double ei2_1=0.006738525414683;
 12         //1975年国际椭球体参数
 13         private const double a_2=6378140.0000000000;
 14         private const double b_2=6356755.2881575287;
 15         private const double c_2=6399596.6519880105;
 16         private const double r_alpha_2=298.257;
 17         private const double e2_2=0.006694384999588;
 18         private const double ei2_2=0.006739501819473;
 19         //WGS-84椭球体参数
 20         private const double a_3=6378137.0000000000;
 21         private const double b_3=6356752.3142;
 22         private const double c_3=6399593.6258;
 23         private const double r_alpha_3=298.257223563;
 24         private const double e2_3=0.0066943799013;
 25         private const double ei2_3=0.00673949674227;
 26         //自定义椭球体参数
 27         private double a_4 = 0;
 28         private double b_4 = 0;
 29         private double c_4 = 0;
 30         private double r_alpha_4 = 0;
 31         private double e2_4 = 0;
 32         private double ei2_4 = 0;
 33         #endregion
 34         #region 成员变量
 35         //标志椭球类型，1,2,3,4分别表示克拉索夫斯基椭球体、1975年国际椭球体、WGS-84椭球体、自定义椭球体
 36         private int m_type = 0;
 37         public double m_a 
 38         {
 39             get
 40             {
 41                 switch (m_type)
 42                 {
 43                     case 1:
 44                         return a_1;
 45                     case 2:
 46                         return a_2;
 47                     case 3:
 48                         return a_3;
 49                     default:
 50                         return a_4;
 51                 }
 52             }
 53         }
 54         public double m_b
 55         {
 56             get
 57             {
 58                 switch (m_type)
 59                 {
 60                     case 1:
 61                         return b_1;
 62                     case 2:
 63                         return b_2;
 64                     case 3:
 65                         return b_3;
 66                     default:
 67                         return b_4;
 68                 }
 69             }
 70         }
 71         public double m_c
 72         {
 73             get
 74             {
 75                 switch (m_type)
 76                 {
 77                     case 1:
 78                         return c_1;
 79                     case 2:
 80                         return c_2;
 81                     case 3:
 82                         return c_3;
 83                     default:
 84                         return c_4;
 85                 }
 86             }
 87         }
 88         public double m_r_alpha
 89         {
 90             get
 91             {
 92                 switch (m_type)
 93                 {
 94                     case 1:
 95                         return r_alpha_1;
 96                     case 2:
 97                         return r_alpha_2;
 98                     case 3:
 99                         return r_alpha_3;
                     default:
                         return r_alpha_4;
                 }
             }
         }
         public double m_e2
         {
             get
             {
                 switch (m_type)
                 {
                     case 1:
                         return e2_1;
                     case 2:
                         return e2_2;
                     case 3:
                         return e2_3;
                     default:
                         return e2_4;
                 }
             }
         }
         public double m_ei2
         {
             get
             {
                 switch (m_type)
                 {
                     case 1:
                         return ei2_1;
                     case 2:
                         return ei2_2;
                     case 3:
                         return ei2_3;
                     default:
                         return ei2_4;
                 }
             }
         }
         #endregion
         #region 公共函数
         public ReferenceEllipsoid(int type)
         {
             m_type = type;
         }
         public ReferenceEllipsoid(double a, double e2)
         {
             m_type = 4;
             a_4 = a;
             e2_4 = e2;
             b_4 = Math.Sqrt(a * a - a * a * e2);
             c_4 = a * a / b_4;
             r_alpha_4 = a / (a - b_4);
             ei2_4 = e2 / (1 - e2);
         }
         public void SetType(int type)
         {
             m_type = type;
         }
         public double Get_t(double B)
         {
             return Math.Tan(B);
         }
         public double Get_eta(double B)
         {
             return Math.Sqrt(m_ei2) * Math.Cos(B);
         }
         public double Get_W(double B)
         {
             return Math.Sqrt(1 - m_e2 * Math.Sin(B) * Math.Sin(B));
         }
         public double Get_V(double B)
         {
             return Math.Sqrt(1 + m_ei2 * Math.Cos(B) * Math.Cos(B));
         }
         //子午圈曲率半径
         public double Get_M(double B)
         {
             return m_a * (1 - m_e2) / Math.Pow(Get_W(B), 3.0);
         }
         //卯酉圈的曲率半径
         public double Get_N(double B)
         {
             return m_a / Get_W(B);
         }
         //任意法截弧的曲率半径
         public double Get_RA(double B, double A)
         {
             return Get_N(B) / (1 + m_ei2 * Math.Pow(Math.Cos(B) * Math.Cos(A), 2.0));
         }
         //平均曲率半径
         public double Get_R(double B)
         {
             return Math.Sqrt(Get_M(B) * Get_N(B));
         }
         //子午线弧长
         public double Get_X(double B)
         {
             double m0 = m_a * (1 - m_e2);
             double m2 = 3.0 * m_e2 * m0 / 2.0;
             double m4 = 5.0 * m_e2 * m2 / 4.0;
             double m6 = 7.0 * m_e2 * m4 / 6.0;
             double m8 = 9.0 * m_e2 * m6 / 8.0;
             double a0 = m0 + m2 / 2.0 + 3.0 * m4 / 8.0 + 5.0 * m6 / 16.0 + 35.0 * m8 / 128;
             double a2 = m2 / 2.0 + m4 / 2.0 + 15.0 * m6 / 32.0 + 7.0 * m8 / 16.0;
             double a4 = m4 / 8.0 + 3.0 * m6 / 16.0 + 7.0 * m8 / 32.0;
             double a6 = m6 / 32.0 + m8 / 16.0;
             double a8 = m8 / 128.0;
             return a0 * B - a2 * Math.Sin(2.0 * B) / 2.0 + a4 * Math.Sin(4.0 * B) / 4.0 -
                 a6 * Math.Sin(6.0 * B) / 6.0 + a8 * Math.Sin(8.0 * B) / 8.0;
         }
         //高斯平均引数法大地主题正算Gauss Average Coordination Direct Solution of Geodesic Problem
         public void GACDS_GP(double L1, double B1, double S, double A12, 
             ref double L2, ref double B2, ref double A21)
         {
             double delta_B0 = S * Math.Cos(A12) / Get_M(B1);
             double delta_L0 = S * Math.Sin(A12) / (Get_N(B1) * Math.Cos(B1));
             double delta_A0 = delta_L0 * Math.Sin(B1);
             double delta_B = delta_B0;
             double delta_L = delta_L0;
             double delta_A = delta_A0;
             do
             {
                 delta_B0 = delta_B;
                 delta_L0 = delta_L;
                 delta_A0 = delta_A;
                 double Bm = B1 + delta_B0 / 2.0;
                 double Lm = L1 + delta_L0 / 2.0;
                 double Am = A12 + delta_A0 / 2.0;
                 double Nm = Get_N(Bm);
                 double Mm = Get_M(Bm);
                 double tm = Get_t(Bm);
                 double eta2m = Math.Pow(Get_eta(Bm), 2.0);
                 delta_B = S * Math.Cos(Am) * (1 + S * S * (Math.Sin(Am) * Math.Sin(Am) * (2 + 3 * tm
                     * tm + 2 * eta2m) + 3 * Math.Cos(Am) * Math.Cos(Am) * eta2m * (tm* tm - 1 - eta2m 
                     - 4 * eta2m * tm * tm)) / (24 *Nm * Nm)) / Mm;
                 delta_L = S * Math.Sin(Am) * (1 + S * S * (tm *tm * Math.Sin(Am) * Math.Sin(Am) - Math.Cos(Am) 
                     * Math.Cos(Am) * (1 + eta2m - 9 * eta2m * tm * tm)) / (24 * Nm * Nm)) / (Nm * Math.Cos(Bm));
                 delta_A = S * Math.Sin(Am) * tm * (1 + S * S * (Math.Cos(Am) * Math.Cos(Am) * (2 + 7 * eta2m 
                     + 9 * eta2m * tm * tm + 5 * eta2m * eta2m) + Math.Sin(Am) * Math.Sin(Am) * (2 + tm * tm 
                     + 2 * eta2m)) / (24 * Nm * Nm)) /Nm;
             }
             while (Math.Abs(delta_B - delta_B0) >= PreciseControl.EPS ||
                 Math.Abs(delta_L - delta_L0) >= PreciseControl.EPS ||
                 Math.Abs(delta_A - delta_A0) >= PreciseControl.EPS);
             B2 = B1 + delta_B0;
             L2 = L1 + delta_L0;
             if (A12 > Math.PI)
             {
                 A21 = A12 + delta_A0 - Math.PI;
             }
             else
             {
                 A21 = A12 + delta_A0 + Math.PI;
             }
         }
 
         //高斯平均引数法大地主题反算的第一个版本，使用材料“06大地主题解算”上的方法，精度不高。
         private void past_GACIS_GP1(double L1, double B1, double L2, double B2,
             ref double S, ref double A12, ref double A21)
         {
             double Bm = (B1 + B2) / 2;
             double deta_B = B2 - B1;
             double deta_L = L2 - L1;
             double Nm = Get_N(Bm);
             double etam = Math.Pow(Get_eta(Bm), 2.0);
             double tm = Get_t(Bm);
             double Vm = Get_V(Bm);
             double r01 = Nm * Math.Cos(Bm);
             double r21 = Nm * Math.Cos(Bm) * (1 + Math.Pow(etam, 2.0) - Math.Pow(3 * etam * tm, 2.0)
                 + Math.Pow(etam, 4.0)) / (24 * Math.Pow(Get_V(Bm), 4));
             double r03 = -Nm * Math.Pow(Math.Cos(Bm), 3) * tm * Get_t(Bm) / 24;
             double s10 = Nm / (Vm * Get_V(Bm));
             double s12 = -Nm * Math.Cos(Bm) * Math.Cos(Bm) * (2 + 3 * tm * tm
                 + 2 * Math.Pow(etam, 2.0)) / (24 * Vm * Vm);
             double s30 = Nm * Math.Cos(Bm) * Math.Cos(Bm) * (2 + 3 * tm * tm
                 + 2 * Math.Pow(etam, 2.0)) / (24 * Vm * Vm);
             double t01 = tm * Math.Cos(Bm);
             double t21 = Math.Cos(Bm) * tm * (2 + 7 * Math.Pow(etam, 2.0) + Math.Pow(3 * etam * tm,2.0))
                 / (24 * Math.Pow(Vm, 4));
             double t03 = Math.Pow(Math.Cos(Bm), 3) * tm * (2 + tm * tm + 2
                 * Math.Pow(etam, 2.0)) / 24;
             double U = r01 * deta_L + r21 * deta_B * deta_B * deta_L + r03 * Math.Pow(deta_L, 3);
             double V = s10 * deta_B + s12 * deta_B * deta_L * deta_L + s30 * Math.Pow(deta_B, 3);
             double deta_A = t01 * deta_B + t21 * deta_B * deta_B * deta_L + t03 * Math.Pow(deta_L, 3);
             double Am = Math.Atan(U / V);
             double C = Math.Abs(V / U);
             double T = 0;
             if (Math.Abs(deta_B) >= Math.Abs(deta_L))
             {
                 T = Math.Atan(Math.Abs(U / V));
             }
             else
             {
                 T = Math.PI / 4 + Math.Atan(Math.Abs((1 - C) / (1 + C)));
             }
             if (deta_B > 0 && deta_L >= 0)
             {
                 Am = T;
             }
             else if (deta_B < 0 && deta_L >= 0)
             {
                 Am = Math.PI - T;
             }
             else if (deta_B <= 0 && deta_L < 0)
             {
                 Am = Math.PI + T;
             }
             else if (deta_B > 0 && deta_L < 0)
             {
                 Am = 2 * Math.PI - T;
             }
             else if (deta_B == 0 && deta_L > 0)
             {
                 Am = Math.PI / 2;
             }
             S = U / Math.Sin(Am);
             A12 = Am - deta_A / 2;
             A21 = Am + deta_A / 2;
             if (A21 >= -Math.PI && A21 < Math.PI)
             {
                 A21 = A21 + Math.PI;
             }
             else
             {
                 A21 = A21 - Math.PI;
             }
         }
 
         //高斯平均引数反算中由Am、Bm、S计算参数alpha
         private double GI_Alpha(double Am, double Bm, double S)
         {
             return Math.Pow(S / Get_N(Bm), 2.0) * (Math.Pow(Math.Sin(Am) * Get_t(Bm), 2.0) - Math.Pow(Math.Cos(Am), 2.0) * (1 + Math.Pow(Get_eta(Bm), 2.0)
                 - Math.Pow(3.0 * Get_t(Bm) * Get_eta(Bm), 2.0))) / 24.0;
         }
         //高斯平均引数反算中由Am、Bm、S计算参数beta
         private double GI_Beta(double Am, double Bm, double S)
         {
             return Math.Pow(S / Get_N(Bm), 2.0) * (Math.Pow(Math.Sin(Am), 2.0) * (2 + 3 * Math.Pow(Get_t(Bm), 2.0) + 2 * Math.Pow(Get_eta(Bm), 2.0)) + 3 *
                 Math.Pow(Math.Cos(Am) * Get_eta(Bm), 2.0) * (-1 + Math.Pow(Get_t(Bm), 2.0) - Math.Pow(Get_eta(Bm), 2.0) - Math.Pow(2.0 * Get_t(Bm) * 
                 Get_eta(Bm), 2.0))) / 24.0;
         }
         //高斯平均引数反算中由Am、Bm、S计算参数gamma
         private double GI_Gamma(double Am, double Bm, double S)
         {
             return Math.Pow(S / Get_N(Bm), 2.0) * (Math.Pow(Math.Cos(Am), 2.0) * (2 + 7 * Math.Pow(Get_eta(Bm), 2.0) + Math.Pow(3.0 * Get_t(Bm) * Get_eta(Bm), 2.0)
                 + 5 * Math.Pow(Get_eta(Bm), 4.0)) + Math.Pow(Math.Sin(Am), 2.0) * (2.0 + Math.Pow(Get_t(Bm), 2.0) + 2.0 * Math.Pow(Get_eta(Bm), 2.0))) / 24.0;
         }
         //高斯平均引数反算中计算参数u和v
         private void GI_uv(ref double u,ref double v,double del_L,double del_B,double Am,double Bm,double S)
         {
             u = del_L * Get_N(Bm) * Math.Cos(Bm) / (1 + GI_Alpha(Am, Bm, S));
             v = del_B * Get_N(Bm) / (Math.Pow(Get_V(Bm), 2.0) * (1 + GI_Beta(Am, Bm, S)));
         }
         //高斯平均引数法大地主题反算Gauss Average Coordination Invert Solution of Geodesic Problem
         public void GACIS_GP(double L1, double B1, double L2, double B2,
             ref double S, ref double A12, ref double A21)
         {
             double del_L = L2 - L1;
             double del_B = B2 - B1;
             double Bm = (B1 + B2) / 2;
             double u0 = del_L * Get_N(Bm) * Math.Cos(Bm);
             double v0 = del_B * Get_N(Bm) / Math.Pow(Get_V(Bm), 2.0);
             double u1 = u0;
             double v1 = v0;
             double Am = 0;
             do
             {
                 u0 = u1;
                 v0 = v1;
                 S = Math.Sqrt(u0 * u0 + v0 * v0);
                 Am = Math.Atan(u0 / v0);
                 GI_uv(ref u1, ref v1, del_L, del_B, Am, Bm, S);
             }
             while (Math.Abs(u1 - u0) > PreciseControl.EPS || Math.Abs(v1 - v0) > PreciseControl.EPS);
             double T = 0;
             if (Math.Abs(del_B) >= Math.Abs(del_L))
             {
                 T = Math.Atan(Math.Abs(u0 / v0));
             }
             else
             {
                 T = Math.PI / 4 + Math.Atan(Math.Abs((1 - Math.Abs(u0 / v0)) / (1 + Math.Abs(u0 / v0))));
             }
             if (del_B > 0 && del_L >= 0)
             {
                 Am = T;
             }
             else if (del_B < 0 && del_L >= 0)
             {
                 Am = Math.PI - T;
             }
             else if (del_B > 0 && del_L < 0)
             {
                 Am = Math.PI + T;
             }
             else if (del_B > 0 && del_L < 0)
             {
                 Am = 2 * Math.PI - T;
             }
             else if (del_B == 0 && del_L > 0)
             {
                 Am = Math.PI / 2;
             }
             double del_A = u0 * Get_t(Bm) * (1 + GI_Gamma(Am, Bm, S)) / Get_N(Bm);
             
             A12 = Am - del_A / 2;
             A21 = Am + del_A / 2;
             if (A21 >= -Math.PI && A21 < Math.PI)
             {
                 A21 = A21 + Math.PI;
             }
             else
             {
                 A21 = A21 - Math.PI;
             }
         }
         //白塞尔大地主题正算函数 Direct Solution of Bessel's Geodetic Problem
         public void DS_BGP(double L1, double B1, double S, double A12,
             ref double L2, ref double B2, ref double A21)
         {
             //将椭球面上的元素投影到球面上
             double u1 = Math.Atan(Math.Sqrt(1 - m_e2) * Math.Tan(B1));
             double A0 = Math.Asin(Math.Cos(u1) * Math.Sin(A12));
             double sigma1 = Math.Atan(Math.Tan(u1) / Math.Cos(A12));
             double k2 = m_ei2 * Math.Cos(A0) * Math.Cos(A0);
             double A = 1 + k2 / 4 - 3 * k2 * k2 / 64 - 5 * k2 * k2 * k2 / 256;
             double B = k2 / 4 - k2 * k2 / 16 + 15 * k2 * k2 * k2 / 512;
             double C = k2 * k2 / 128 - 3 * k2 * k2 * k2 / 512;
             double alpha = 1 / (m_b * A);
             double beta = B / A;
             double gamma = C / A;
             double sigma0 = alpha * S;
             double sigma = sigma0;
             do
             {
                 sigma0 = sigma;
                 sigma = alpha * S + beta * Math.Sin(sigma0) * Math.Cos(2 * sigma1 + sigma0)
                     + gamma * Math.Sin(2 * sigma0) * Math.Cos(4 * sigma1 + 2 * sigma0);
             }
             while (Math.Abs(sigma - sigma0) >= PreciseControl.EPS);
             //解算球面三角形
             A21 = Math.Atan(Math.Cos(u1) * Math.Sin(A12) / (Math.Cos(u1) * Math.Cos(sigma) 
                 * Math.Cos(A12) - Math.Sin(u1) * Math.Sin(sigma)));
             double u2 = Math.Asin(Math.Sin(u1) * Math.Cos(sigma) + Math.Cos(u1) 
                 * Math.Cos(A12) * Math.Sin(sigma));
             double lambda = Math.Atan(Math.Sin(A12) * Math.Sin(sigma) / (Math.Cos(u1) 
                 * Math.Cos(sigma) - Math.Sin(u1) * Math.Sin(sigma) * Math.Cos(A12)));
             //判定象限
             double Abs_lambda = 0;
             double Abs_A21 = 0;
             Operation_Angle.ToFirstQuadrant(lambda,ref Abs_lambda);
             Operation_Angle.ToFirstQuadrant(A21, ref Abs_A21);
             if (Math.Sin(A12) >= 0 && Math.Tan(lambda) >= 0)
             {
                 lambda = Abs_lambda;
             }
             else if (Math.Sin(A12) >= 0 && Math.Tan(lambda) < 0)
             {
                 lambda = Math.PI - Abs_lambda;
             }
             else if (Math.Sin(A12) < 0 && Math.Tan(lambda) >= 0)
             {
                 lambda = -Abs_lambda;
             }
             else if (Math.Sin(A12) < 0 && Math.Tan(lambda) < 0)
             {
                 lambda = Abs_lambda - Math.PI;
             }
             if (Math.Sin(A12) < 0 && Math.Tan(A21) >= 0)
             {
                 A21 = Abs_A21;
             }
             else if (Math.Sin(A12) < 0 && Math.Tan(A21) < 0)
             {
                 A21 = Math.PI - Abs_A21;
             }
             else if (Math.Sin(A12) >= 0 && Math.Tan(A21) >= 0)
             {
                 A21 = Math.PI + Abs_A21;
             }
             else if (Math.Sin(A12) >= 0 && Math.Tan(A21) < 0)
             {
                 A21 = 2 * Math.PI - Abs_A21;
             }
             //将球面元素换算到椭球面上
             B2 = Math.Atan(Math.Sqrt(1 + m_ei2) * Math.Tan(u2));
             double ki2 = m_e2 * Math.Cos(A0) * Math.Cos(A0);
             double alphai = (m_e2 / 2 + m_e2 * m_e2 / 8 + m_e2 * m_e2 * m_e2 / 16) 
                 - m_e2 * (1 + m_e2) * ki2 / 16 + 3 * m_e2 * ki2 * ki2 / 128;
             double betai = m_e2 * (1 + m_e2) * ki2 / 16 - m_e2 * ki2 * ki2 / 32;
             double gammai = m_e2 * ki2 * ki2 / 256;
             double l = lambda - Math.Sin(A0) * (alphai * sigma + betai * Math.Sin(sigma) 
                 * Math.Cos(2 * sigma1 + sigma) + gammai * Math.Sin(2 * sigma) * 
                 Math.Cos(4 * sigma1 + 2 * sigma));
             L2 = L1 + l;
         }
         //白塞尔大地主题反算函数 Inverse Solution of Bessel's Geodetic Problem
         public void IS_BGP(double L1, double B1, double L2, double B2,
             ref double S, ref double A12, ref double A21)
         {
             //将椭球面元素投影到球面上
             double u1 = Math.Atan(Math.Sqrt(1 - m_e2) * Math.Tan(B1));
             double u2 = Math.Atan(Math.Sqrt(1 - m_e2) * Math.Tan(B2));
             double l = L2 - L1;
             double a1 = Math.Sin(u1) * Math.Sin(u2);
             double a2 = Math.Cos(u1) * Math.Cos(u2);
             double b1 = Math.Cos(u1) * Math.Sin(u2);
             double b2 = Math.Sin(u1) * Math.Cos(u2);
             //迭代求lambda
             double lambda = l;
             double lambda0 = l;
             double sigma = 0;
             double sigma1 = 0;
             double A0 = 0;
             do
             {
                 lambda0 = lambda;
                 double p = Math.Cos(u2) * Math.Sin(lambda);
                 double q = b1 - b2 * Math.Cos(lambda);
                 A12 = Math.Atan(p / q);
                 double Abs_A12 = 0;
                 Operation_Angle.ToFirstQuadrant(A12, ref Abs_A12);
                 if (p >= 0 && q >= 0)
                 {
                     A12 = Abs_A12;
                 }
                 else if (p >= 0 && q < 0)
                 {
                     A12 = Math.PI - Abs_A12;
                 }
                 else if (p < 0 && q < 0)
                 {
                     A12 = Math.PI + Abs_A12;
                 }
                 else if (p < 0 && q >= 0)
                 {
                     A12 = 2 * Math.PI - Abs_A12;
                 }
                 double sins = Math.Cos(u2) * Math.Sin(lambda) * Math.Sin(A12) + (Math.Cos(u1) * Math.Sin(u2)
                     - Math.Sin(u1) * Math.Cos(u2) * Math.Cos(lambda)) * Math.Cos(A12);
                 double coss = Math.Sin(u1) * Math.Sin(u2) + Math.Cos(u1) * Math.Cos(u2) * Math.Cos(lambda);
                 sigma = Math.Atan(sins / coss);
                 double Abs_sigma = 0;
                 Operation_Angle.ToFirstQuadrant(sigma, ref Abs_sigma);
                 if (coss >= 0)
                 {
                     sigma = Abs_sigma;
                 }
                 else
                 {
                     sigma = Math.PI - Abs_sigma;
                 }
                 A0 = Math.Asin(Math.Cos(u1) * Math.Sin(A12));
                 sigma1 = Math.Atan(Math.Tan(u1) / Math.Cos(A12));
                 double ki2 = m_e2 * Math.Cos(A0) * Math.Cos(A0);
                 double alphai = (m_e2 / 2 + m_e2 * m_e2 / 8 + m_e2 * m_e2 * m_e2 / 16)
                     - m_e2 * (1 + m_e2) * ki2 / 16 + 3 * m_e2 * ki2 * ki2 / 128;
                 double betai = m_e2 * (1 + m_e2) * ki2 / 16 - m_e2 * ki2 * ki2 / 32;
                 double gammai = m_e2 * ki2 * ki2 / 256;
                 lambda = l + Math.Sin(A0) * (alphai * sigma + betai * Math.Sin(sigma)
                     * Math.Cos(2.0 * sigma1 + sigma) + gammai * Math.Sin(2.0 * sigma)
                     * Math.Cos(4.0 * sigma1 + 2.0 * sigma));
             }
             while (Math.Abs(lambda - lambda0) >= PreciseControl.EPS);
             //将球面元素换算到椭球面上
             double k2 = m_ei2 * Math.Cos(A0) * Math.Cos(A0);
             double A = 1 + k2 / 4 - 3 * k2 * k2 / 64 - 5 * k2 * k2 * k2 / 256;
             double B = k2 / 4 - k2 * k2 / 16 + 15 * k2 * k2 * k2 / 512;
             double C = k2 * k2 / 128 - 3 * k2 * k2 * k2 / 512;
             double alpha = 1 / (m_b * A);
             double beta = B / A;
             double gamma = C / A;
             S = (sigma - beta * Math.Sin(sigma) * Math.Cos(2 * sigma1 + sigma) - gamma
                 * Math.Sin(2 * sigma) * Math.Cos(4 * sigma1 + 2 * sigma)) / alpha;
             A21 = Math.Atan(Math.Cos(u1) * Math.Sin(lambda) / (Math.Cos(u1) * Math.Sin(u2)
                 * Math.Cos(lambda) - Math.Sin(u1) * Math.Cos(u2)));
             double Abs_A21 = 0;
             Operation_Angle.ToFirstQuadrant(A21, ref Abs_A21);
             if (Math.Sin(A12) < 0 && Math.Tan(A21) >= 0)
             {
                 A21 = Abs_A21;
             }
             else if (Math.Sin(A12) < 0 && Math.Tan(A21) < 0)
             {
                 A21 = Math.PI - Abs_A21;
             }
             else if (Math.Sin(A12) >= 0 && Math.Tan(A21) >= 0)
             {
                 A21 = Math.PI + Abs_A21;
             }
             else if (Math.Sin(A12) >= 0 && Math.Tan(A21) < 0)
             {
                 A21 = 2 * Math.PI - Abs_A21;
             }
         }
         //Convert Measured Value of Distance to Ellipsoid 将测量值归算至椭球面
         public double CMVTE_Distance(double D, double H1, double H2, double B1, double L1, double B2, double L2)
         {
             //此函数假设不知道两点连线的方位角，却知道两点的大地坐标，由大地坐标算出方位角。
             double A12 = 0;
             double S = 0;
             double A21 = 0;
             //由于距离一般较短，采用高斯平均引数法
             GACIS_GP(L1, B1, L2, B2, ref S, ref A12, ref A21);
             //计算起点曲率半径
             double RA = Get_N(B1) / (1 + m_ei2 * Math.Pow(Math.Cos(B1) * Math.Cos(A12), 2));
             //三个临时变量
             double temp1 = 1 - Math.Pow((H2 - H1) / D, 2);
             double temp2 = (1 + H1 / RA) * (1 + H2 / RA);
             double temp3 = Math.Pow(D, 3) / (24 * RA * RA);
             //距离换算公式
             return D * Math.Sqrt(temp1 / temp2) + temp3;
         }
         public double CMVTE_Distance(double D, double H1, double H2, double B1,double A12)
         {
             //此函数假设知道方位角。
 
             //计算起点曲率半径
             double RA = Get_N(B1) / (1 + m_ei2 * Math.Pow(Math.Cos(B1) * Math.Cos(A12), 2));
             //三个临时变量
             double temp1 = 1 - Math.Pow((H2 - H1) / D, 2);
             double temp2 = (1 + H1 / RA) * (1 + H2 / RA);
             double temp3 = Math.Pow(D, 3) / (24 * RA * RA);
             //距离换算公式
             return D * Math.Sqrt(temp1 / temp2) + temp3;
         }
         #endregion
 

坐标投影类和高斯坐标投影类：

 abstract class MapProjection
     {
         protected ReferenceEllipsoid m_Ellipsoid;
     }

  1 class GaussProjection : MapProjection
  2     {
  3         #region 成员变量
  4         private double m_D;                      //每带的宽度（弧度）
  5         private double m_Starting_L;             //起点经度
  6         #endregion
  7         #region 接口函数
  8         public GaussProjection(ReferenceEllipsoid RE, double D, double Starting_L)
  9         {
 10             m_Ellipsoid = RE;
 11             m_D = D;
 12             m_Starting_L = Starting_L;
 13         }
 14         //坐标正算
 15         public void Positive_Computation(double L, double B, ref double x, ref double y, ref int n)
 16         {
 17             //求带号
 18             n = Convert.ToInt32((L - m_Starting_L) / m_D);
 19             //求中央经度
 20             double L0 = m_Starting_L + (Convert.ToDouble(n) - 0.5) * m_D;
 21             //求点与中央子午线的经度差
 22             double l = L - L0;
 23             //辅助变量
 24             double m = Math.Cos(B) * l;
 25             double X = m_Ellipsoid.Get_X(B);
 26             double t = m_Ellipsoid.Get_t(B);
 27             double eta2 = Math.Pow(m_Ellipsoid.Get_eta(B), 2.0);
 28             double N = m_Ellipsoid.Get_N(B);
 29             //坐标正算公式
 30             x = X + N * t * ((0.5 + ((5 - t * t + 9 * eta2 + 4 * eta2 * eta2) / 24 + (61 - 58 * t * t 
 31                 + Math.Pow(t, 4)) * m * m / 720) * m * m) * m * m);
 32             y = N * ((1 + ((1 - t * t + eta2) / 6 + (5 - 18 * t * t + Math.Pow(t, 4) + 14 * eta2 - 58 
 33                 * eta2 * t * t) * m * m / 120) * m * m) * m);
 34         }
 35         //坐标反算
 36         public void Negative_Computation(int n, double x, double y, ref double L, ref double B)
 37         {
 38             //此函数采用迭代算法
 39 
 40             //辅助变量
 41             double a = m_Ellipsoid.m_a;
 42             double ei2 = m_Ellipsoid.m_ei2;
 43             double e2 = m_Ellipsoid.m_e2;
 44             double c = m_Ellipsoid.m_c;
 45             double beta0 = 1 - 3 * ei2 / 4 + 45 * ei2 * ei2 / 64 - 175 * Math.Pow(ei2, 3) / 256 
 46                 + 11025 * Math.Pow(ei2, 4) / 16384;
 47             double beta2 = beta0 - 1;
 48             double beta4 = 15 * ei2 * ei2 / 32 - 175 * Math.Pow(ei2, 3) / 384 + 3675 * Math.Pow(ei2, 4) 
 49                 / 8192;
 50             double beta6 = -35 * Math.Pow(ei2, 3) / 96 + 735 * Math.Pow(ei2, 4) / 2048;
 51             double beta8 = 315 * Math.Pow(ei2, 4) / 1024;
 52             //赋初值
 53             double B0 = x / (c * beta0);
 54             double a10 =a * Math.Cos(B0) / Math.Sqrt(1 - e2 * Math.Sin(B0) * Math.Sin(B0));
 55             double l0 = y / a10;
 56             double B1 = B0;
 57             double a11 = a10;
 58             double l1 = l0;
 59             do
 60             {
 61                 B0 = B1;
 62                 a10 = a11;
 63                 l0 = l1;
 64                 double N = m_Ellipsoid.Get_N(B0);
 65                 double t = m_Ellipsoid.Get_t(B0);
 66                 double eta2 = Math.Pow(m_Ellipsoid.Get_eta(B0), 2.0);
 67                 double a2 = N * Math.Sin(B0) * Math.Cos(B0) / 2;   
 68                 double a4 = N * Math.Sin(B0) * Math.Pow(Math.Cos(B0), 3) * (5 - t * t + 9 * eta2 + 4 
 69                     * eta2 * eta2) / 24;
 70                 double a6 = N * Math.Sin(B0) * Math.Pow(Math.Cos(B0), 5) * (61 - 58 * t * t + Math.Pow(t, 4))
 71                     / 720;
 72                 double FxB = (c * beta2 + (c * beta4 + (c * beta6 + c * beta8 * Math.Pow(Math.Cos(B0), 2)) 
 73                     * Math.Pow(Math.Cos(B0), 2)) * Math.Pow(Math.Cos(B0), 2)) * Math.Sin(B0) * Math.Cos(B0);
 74                 double FxBl = a2 * l0 * l0 + a4 * Math.Pow(l0, 4) + a6 * Math.Pow(l0, 6);
 75                 B1 = (x - FxB - FxBl) / (c * beta0);
 76                 a11 = a * Math.Cos(B1) / Math.Sqrt(1 - e2 * Math.Sin(B1) * Math.Sin(B1));
 77                 double N1 = m_Ellipsoid.Get_N(B1);
 78                 double t1 = m_Ellipsoid.Get_t(B1);
 79                 double eta21 = Math.Pow(m_Ellipsoid.Get_eta(B1), 2.0);
 80                 double a3 = N1 * Math.Pow(Math.Cos(B1), 3) * (1 - t1 * t1 + eta21) / 6;
 81                 double a5 = N1 * Math.Pow(Math.Cos(B1), 5) * (5 - 18 * t1 * t1 + Math.Pow(t1, 4) + 14 * 
 82                     eta21 - 58 * eta21 * t1 * t1) / 120;
 83                 double FyBl = a3 * Math.Pow(l0, 3) + a5 * Math.Pow(l0, 5);
 84                 l1 = (y - FyBl) / a11;
 85             }
 86             while (Math.Abs(B1 - B0) > PreciseControl.EPS || Math.Abs(l1 - l0) > PreciseControl.EPS);
 87             double L0 = m_Starting_L + (n - 0.5) * m_D;
 88             L = L0 + l0;
 89             B = B0;
 90         }
 91         //对于6度带，将坐标转换成国家统一坐标形式 National Coordinate System
 92         static public void To_NCS(int n, double x, double y, ref double NCS_x, ref double NCS_y)
 93         {
 94             NCS_x = x;
 95             NCS_y = n * 1000000 + y + 500000;
 96         }
 97         static public void From_NCS(double NCS_x, double NCS_y, ref int n, ref double x, ref double y)
 98         {
 99             x = NCS_x;
             n = (int)(NCS_y / 1000000);
             y = NCS_y - n * 1000000 - 500000;
         }
 
         //距离改化公式，将椭球面是的距离化算至平面距离
         public double Convert_Distance(double S, double L1, double B1, double L2, double B2)
         {
             int n = 0;
             double x1 = 0;
             double y1 = 0;
             double x2 = 0;
             double y2 = 0;
             Positive_Computation(L1, B1, ref x1, ref y1, ref n);
             Positive_Computation(L2, B2, ref x2, ref y2, ref n);
             double R1 = m_Ellipsoid.Get_R(B1);
             double R2 = m_Ellipsoid.Get_R(B2);
             double Rm = (R1 + R2) / 2.0;
             double ym = (y1 + y2) / 2.0;
             double deta_y = y1 - y2;
             return S * (1 + ym * ym / (2 * Rm * Rm) + deta_y * deta_y / (24 * Rm * Rm) + Math.Pow(ym / Rm, 4) / 24);
         }
         #endregion
 

两个辅助类：

 1 //角度运算类
 2     class Operation_Angle
 3     {
 4         //角度转弧度
 5         static public bool AngleToRadian(int degree, int minite, double second, ref double radian)
 6         {
 7             if (degree < 0 || degree >= 360)
 8             {
 9                 return false;
             }
             if (minite < 0 || minite >= 60)
             {
                 return false;
             }
             if (second < 0 || second >= 60)
             {
                 return false;
             }
             double temp = degree + minite / 60.0 + second / 3600.0;
             radian = Math.PI * temp / 180.0;
             return true;
         }
         static public bool  AngleToRadian(double degree, ref double radian)
         {
             if (degree < 0 || degree >= 360)
             {
                 return false;
             }
             radian = Math.PI * degree / 180.0;
             return true;
         }
         //弧度转角度
         static public bool RadianToAngle(double radian, ref int degree, ref int minite, ref double second)
         {
             if (radian < 0 || radian >= 2 * Math.PI)
             {
                 return false;
             }
             double temp = 180.0 * radian / Math.PI;
             degree = (int)temp;
             temp = (temp - (double)degree) * 60;
             minite = (int)temp;
             temp = (temp - (double)minite) * 60;
             second = temp;
             return true;
         }
         //将一个角度通过加减90度化到第一象限
         static public void ToFirstQuadrant(double radian, ref double f_radian)
         {
             f_radian = radian;
             while (f_radian >= Math.PI / 2)
             {
                 f_radian -= Math.PI / 2;
             }
             while (f_radian < 0)
             {
                 f_radian += Math.PI / 2;
             }
         }
         static public void ToFirstQuadrant(int degree, int minite, double second,
             ref int f_degree, ref int f_minite, ref double f_second)
         {
             f_degree = degree;
             f_minite = minite;
             f_second = second;
             while (f_degree >= 90)
             {
                 f_degree -= 90;
             }
             while (f_degree < 0)
             {
                 f_degree += 90;
             }
         }
 

 1 public class PreciseControl
 2     {
 3         //程序迭代精度控制
 4         public static double EPS
 5         {
 6             get
 7             {
 8                 return Math.Pow(10.0, -10.0);
 9             }
         }
     }