结构示意图：

参考椭球类：

Code
  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;
100                     default:
101                         return r_alpha_4;
102                 }
103             }
104         }
105         public double m_e2
106         {
107             get
108             {
109                 switch (m_type)
110                 {
111                     case 1:
112                         return e2_1;
113                     case 2:
114                         return e2_2;
115                     case 3:
116                         return e2_3;
117                     default:
118                         return e2_4;
119                 }
120             }
121         }
122         public double m_ei2
123         {
124             get
125             {
126                 switch (m_type)
127                 {
128                     case 1:
129                         return ei2_1;
130                     case 2:
131                         return ei2_2;
132                     case 3:
133                         return ei2_3;
134                     default:
135                         return ei2_4;
136                 }
137             }
138         }
139         #endregion
140         #region 公共函数
141         public ReferenceEllipsoid(int type)
142         {
143             m_type = type;
144         }
145         public ReferenceEllipsoid(double a, double e2)
146         {
147             m_type = 4;
148             a_4 = a;
149             e2_4 = e2;
150             b_4 = Math.Sqrt(a * a - a * a * e2);
151             c_4 = a * a / b_4;
152             r_alpha_4 = a / (a - b_4);
153             ei2_4 = e2 / (1 - e2);
154         }
155         public void SetType(int type)
156         {
157             m_type = type;
158         }
159         public double Get_t(double B)
160         {
161             return Math.Tan(B);
162         }
163         public double Get_eta(double B)
164         {
165             return Math.Sqrt(m_ei2) * Math.Cos(B);
166         }
167         public double Get_W(double B)
168         {
169             return Math.Sqrt(1 - m_e2 * Math.Sin(B) * Math.Sin(B));
170         }
171         public double Get_V(double B)
172         {
173             return Math.Sqrt(1 + m_ei2 * Math.Cos(B) * Math.Cos(B));
174         }
175         //子午圈曲率半径
176         public double Get_M(double B)
177         {
178             return m_a * (1 - m_e2) / Math.Pow(Get_W(B), 3.0);
179         }
180         //卯酉圈的曲率半径
181         public double Get_N(double B)
182         {
183             return m_a / Get_W(B);
184         }
185         //任意法截弧的曲率半径
186         public double Get_RA(double B, double A)
187         {
188             return Get_N(B) / (1 + m_ei2 * Math.Pow(Math.Cos(B) * Math.Cos(A), 2.0));
189         }
190         //平均曲率半径
191         public double Get_R(double B)
192         {
193             return Math.Sqrt(Get_M(B) * Get_N(B));
194         }
195         //子午线弧长
196         public double Get_X(double B)
197         {
198             double m0 = m_a * (1 - m_e2);
199             double m2 = 3.0 * m_e2 * m0 / 2.0;
200             double m4 = 5.0 * m_e2 * m2 / 4.0;
201             double m6 = 7.0 * m_e2 * m4 / 6.0;
202             double m8 = 9.0 * m_e2 * m6 / 8.0;
203             double a0 = m0 + m2 / 2.0 + 3.0 * m4 / 8.0 + 5.0 * m6 / 16.0 + 35.0 * m8 / 128;
204             double a2 = m2 / 2.0 + m4 / 2.0 + 15.0 * m6 / 32.0 + 7.0 * m8 / 16.0;
205             double a4 = m4 / 8.0 + 3.0 * m6 / 16.0 + 7.0 * m8 / 32.0;
206             double a6 = m6 / 32.0 + m8 / 16.0;
207             double a8 = m8 / 128.0;
208             return a0 * B - a2 * Math.Sin(2.0 * B) / 2.0 + a4 * Math.Sin(4.0 * B) / 4.0 -
209                 a6 * Math.Sin(6.0 * B) / 6.0 + a8 * Math.Sin(8.0 * B) / 8.0;
210         }
211         //高斯平均引数法大地主题正算Gauss Average Coordination Direct Solution of Geodesic Problem
212         public void GACDS_GP(double L1, double B1, double S, double A12, 
213             ref double L2, ref double B2, ref double A21)
214         {
215             double delta_B0 = S * Math.Cos(A12) / Get_M(B1);
216             double delta_L0 = S * Math.Sin(A12) / (Get_N(B1) * Math.Cos(B1));
217             double delta_A0 = delta_L0 * Math.Sin(B1);
218             double delta_B = delta_B0;
219             double delta_L = delta_L0;
220             double delta_A = delta_A0;
221             do
222             {
223                 delta_B0 = delta_B;
224                 delta_L0 = delta_L;
225                 delta_A0 = delta_A;
226                 double Bm = B1 + delta_B0 / 2.0;
227                 double Lm = L1 + delta_L0 / 2.0;
228                 double Am = A12 + delta_A0 / 2.0;
229                 double Nm = Get_N(Bm);
230                 double Mm = Get_M(Bm);
231                 double tm = Get_t(Bm);
232                 double eta2m = Math.Pow(Get_eta(Bm), 2.0);
233                 delta_B = S * Math.Cos(Am) * (1 + S * S * (Math.Sin(Am) * Math.Sin(Am) * (2 + 3 * tm
234                     * tm + 2 * eta2m) + 3 * Math.Cos(Am) * Math.Cos(Am) * eta2m * (tm* tm - 1 - eta2m 
235                     - 4 * eta2m * tm * tm)) / (24 *Nm * Nm)) / Mm;
236                 delta_L = S * Math.Sin(Am) * (1 + S * S * (tm *tm * Math.Sin(Am) * Math.Sin(Am) - Math.Cos(Am) 
237                     * Math.Cos(Am) * (1 + eta2m - 9 * eta2m * tm * tm)) / (24 * Nm * Nm)) / (Nm * Math.Cos(Bm));
238                 delta_A = S * Math.Sin(Am) * tm * (1 + S * S * (Math.Cos(Am) * Math.Cos(Am) * (2 + 7 * eta2m 
239                     + 9 * eta2m * tm * tm + 5 * eta2m * eta2m) + Math.Sin(Am) * Math.Sin(Am) * (2 + tm * tm 
240                     + 2 * eta2m)) / (24 * Nm * Nm)) /Nm;
241             }
242             while (Math.Abs(delta_B - delta_B0) >= PreciseControl.EPS ||
243                 Math.Abs(delta_L - delta_L0) >= PreciseControl.EPS ||
244                 Math.Abs(delta_A - delta_A0) >= PreciseControl.EPS);
245             B2 = B1 + delta_B0;
246             L2 = L1 + delta_L0;
247             if (A12 > Math.PI)
248             {
249                 A21 = A12 + delta_A0 - Math.PI;
250             }
251             else
252             {
253                 A21 = A12 + delta_A0 + Math.PI;
254             }
255         }
256 
257         //高斯平均引数法大地主题反算的第一个版本，使用材料“06大地主题解算”上的方法，精度不高。
258         private void past_GACIS_GP1(double L1, double B1, double L2, double B2,
259             ref double S, ref double A12, ref double A21)
260         {
261             double Bm = (B1 + B2) / 2;
262             double deta_B = B2 - B1;
263             double deta_L = L2 - L1;
264             double Nm = Get_N(Bm);
265             double etam = Math.Pow(Get_eta(Bm), 2.0);
266             double tm = Get_t(Bm);
267             double Vm = Get_V(Bm);
268             double r01 = Nm * Math.Cos(Bm);
269             double r21 = Nm * Math.Cos(Bm) * (1 + Math.Pow(etam, 2.0) - Math.Pow(3 * etam * tm, 2.0)
270                 + Math.Pow(etam, 4.0)) / (24 * Math.Pow(Get_V(Bm), 4));
271             double r03 = -Nm * Math.Pow(Math.Cos(Bm), 3) * tm * Get_t(Bm) / 24;
272             double s10 = Nm / (Vm * Get_V(Bm));
273             double s12 = -Nm * Math.Cos(Bm) * Math.Cos(Bm) * (2 + 3 * tm * tm
274                 + 2 * Math.Pow(etam, 2.0)) / (24 * Vm * Vm);
275             double s30 = Nm * Math.Cos(Bm) * Math.Cos(Bm) * (2 + 3 * tm * tm
276                 + 2 * Math.Pow(etam, 2.0)) / (24 * Vm * Vm);
277             double t01 = tm * Math.Cos(Bm);
278             double t21 = Math.Cos(Bm) * tm * (2 + 7 * Math.Pow(etam, 2.0) + Math.Pow(3 * etam * tm,2.0))
279                 / (24 * Math.Pow(Vm, 4));
280             double t03 = Math.Pow(Math.Cos(Bm), 3) * tm * (2 + tm * tm + 2
281                 * Math.Pow(etam, 2.0)) / 24;
282             double U = r01 * deta_L + r21 * deta_B * deta_B * deta_L + r03 * Math.Pow(deta_L, 3);
283             double V = s10 * deta_B + s12 * deta_B * deta_L * deta_L + s30 * Math.Pow(deta_B, 3);
284             double deta_A = t01 * deta_B + t21 * deta_B * deta_B * deta_L + t03 * Math.Pow(deta_L, 3);
285             double Am = Math.Atan(U / V);
286             double C = Math.Abs(V / U);
287             double T = 0;
288             if (Math.Abs(deta_B) >= Math.Abs(deta_L))
289             {
290                 T = Math.Atan(Math.Abs(U / V));
291             }
292             else
293             {
294                 T = Math.PI / 4 + Math.Atan(Math.Abs((1 - C) / (1 + C)));
295             }
296             if (deta_B > 0 && deta_L >= 0)
297             {
298                 Am = T;
299             }
300             else if (deta_B < 0 && deta_L >= 0)
301             {
302                 Am = Math.PI - T;
303             }
304             else if (deta_B <= 0 && deta_L < 0)
305             {
306                 Am = Math.PI + T;
307             }
308             else if (deta_B > 0 && deta_L < 0)
309             {
310                 Am = 2 * Math.PI - T;
311             }
312             else if (deta_B == 0 && deta_L > 0)
313             {
314                 Am = Math.PI / 2;
315             }
316             S = U / Math.Sin(Am);
317             A12 = Am - deta_A / 2;
318             A21 = Am + deta_A / 2;
319             if (A21 >= -Math.PI && A21 < Math.PI)
320             {
321                 A21 = A21 + Math.PI;
322             }
323             else
324             {
325                 A21 = A21 - Math.PI;
326             }
327         }
328 
329         //高斯平均引数反算中由Am、Bm、S计算参数alpha
330         private double GI_Alpha(double Am, double Bm, double S)
331         {
332             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)
333                 - Math.Pow(3.0 * Get_t(Bm) * Get_eta(Bm), 2.0))) / 24.0;
334         }
335         //高斯平均引数反算中由Am、Bm、S计算参数beta
336         private double GI_Beta(double Am, double Bm, double S)
337         {
338             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 *
339                 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) * 
340                 Get_eta(Bm), 2.0))) / 24.0;
341         }
342         //高斯平均引数反算中由Am、Bm、S计算参数gamma
343         private double GI_Gamma(double Am, double Bm, double S)
344         {
345             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)
346                 + 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;
347         }
348         //高斯平均引数反算中计算参数u和v
349         private void GI_uv(ref double u,ref double v,double del_L,double del_B,double Am,double Bm,double S)
350         {
351             u = del_L * Get_N(Bm) * Math.Cos(Bm) / (1 + GI_Alpha(Am, Bm, S));
352             v = del_B * Get_N(Bm) / (Math.Pow(Get_V(Bm), 2.0) * (1 + GI_Beta(Am, Bm, S)));
353         }
354         //高斯平均引数法大地主题反算Gauss Average Coordination Invert Solution of Geodesic Problem
355         public void GACIS_GP(double L1, double B1, double L2, double B2,
356             ref double S, ref double A12, ref double A21)
357         {
358             double del_L = L2 - L1;
359             double del_B = B2 - B1;
360             double Bm = (B1 + B2) / 2;
361             double u0 = del_L * Get_N(Bm) * Math.Cos(Bm);
362             double v0 = del_B * Get_N(Bm) / Math.Pow(Get_V(Bm), 2.0);
363             double u1 = u0;
364             double v1 = v0;
365             double Am = 0;
366             do
367             {
368                 u0 = u1;
369                 v0 = v1;
370                 S = Math.Sqrt(u0 * u0 + v0 * v0);
371                 Am = Math.Atan(u0 / v0);
372                 GI_uv(ref u1, ref v1, del_L, del_B, Am, Bm, S);
373             }
374             while (Math.Abs(u1 - u0) > PreciseControl.EPS || Math.Abs(v1 - v0) > PreciseControl.EPS);
375             double T = 0;
376             if (Math.Abs(del_B) >= Math.Abs(del_L))
377             {
378                 T = Math.Atan(Math.Abs(u0 / v0));
379             }
380             else
381             {
382                 T = Math.PI / 4 + Math.Atan(Math.Abs((1 - Math.Abs(u0 / v0)) / (1 + Math.Abs(u0 / v0))));
383             }
384             if (del_B > 0 && del_L >= 0)
385             {
386                 Am = T;
387             }
388             else if (del_B < 0 && del_L >= 0)
389             {
390                 Am = Math.PI - T;
391             }
392             else if (del_B > 0 && del_L < 0)
393             {
394                 Am = Math.PI + T;
395             }
396             else if (del_B > 0 && del_L < 0)
397             {
398                 Am = 2 * Math.PI - T;
399             }
400             else if (del_B == 0 && del_L > 0)
401             {
402                 Am = Math.PI / 2;
403             }
404             double del_A = u0 * Get_t(Bm) * (1 + GI_Gamma(Am, Bm, S)) / Get_N(Bm);
405             
406             A12 = Am - del_A / 2;
407             A21 = Am + del_A / 2;
408             if (A21 >= -Math.PI && A21 < Math.PI)
409             {
410                 A21 = A21 + Math.PI;
411             }
412             else
413             {
414                 A21 = A21 - Math.PI;
415             }
416         }
417         //白塞尔大地主题正算函数 Direct Solution of Bessel's Geodetic Problem
418         public void DS_BGP(double L1, double B1, double S, double A12,
419             ref double L2, ref double B2, ref double A21)
420         {
421             //将椭球面上的元素投影到球面上
422             double u1 = Math.Atan(Math.Sqrt(1 - m_e2) * Math.Tan(B1));
423             double A0 = Math.Asin(Math.Cos(u1) * Math.Sin(A12));
424             double sigma1 = Math.Atan(Math.Tan(u1) / Math.Cos(A12));
425             double k2 = m_ei2 * Math.Cos(A0) * Math.Cos(A0);
426             double A = 1 + k2 / 4 - 3 * k2 * k2 / 64 - 5 * k2 * k2 * k2 / 256;
427             double B = k2 / 4 - k2 * k2 / 16 + 15 * k2 * k2 * k2 / 512;
428             double C = k2 * k2 / 128 - 3 * k2 * k2 * k2 / 512;
429             double alpha = 1 / (m_b * A);
430             double beta = B / A;
431             double gamma = C / A;
432             double sigma0 = alpha * S;
433             double sigma = sigma0;
434             do
435             {
436                 sigma0 = sigma;
437                 sigma = alpha * S + beta * Math.Sin(sigma0) * Math.Cos(2 * sigma1 + sigma0)
438                     + gamma * Math.Sin(2 * sigma0) * Math.Cos(4 * sigma1 + 2 * sigma0);
439             }
440             while (Math.Abs(sigma - sigma0) >= PreciseControl.EPS);
441             //解算球面三角形
442             A21 = Math.Atan(Math.Cos(u1) * Math.Sin(A12) / (Math.Cos(u1) * Math.Cos(sigma) 
443                 * Math.Cos(A12) - Math.Sin(u1) * Math.Sin(sigma)));
444             double u2 = Math.Asin(Math.Sin(u1) * Math.Cos(sigma) + Math.Cos(u1) 
445                 * Math.Cos(A12) * Math.Sin(sigma));
446             double lambda = Math.Atan(Math.Sin(A12) * Math.Sin(sigma) / (Math.Cos(u1) 
447                 * Math.Cos(sigma) - Math.Sin(u1) * Math.Sin(sigma) * Math.Cos(A12)));
448             //判定象限
449             double Abs_lambda = 0;
450             double Abs_A21 = 0;
451             Operation_Angle.ToFirstQuadrant(lambda,ref Abs_lambda);
452             Operation_Angle.ToFirstQuadrant(A21, ref Abs_A21);
453             if (Math.Sin(A12) >= 0 && Math.Tan(lambda) >= 0)
454             {
455                 lambda = Abs_lambda;
456             }
457             else if (Math.Sin(A12) >= 0 && Math.Tan(lambda) < 0)
458             {
459                 lambda = Math.PI - Abs_lambda;
460             }
461             else if (Math.Sin(A12) < 0 && Math.Tan(lambda) >= 0)
462             {
463                 lambda = -Abs_lambda;
464             }
465             else if (Math.Sin(A12) < 0 && Math.Tan(lambda) < 0)
466             {
467                 lambda = Abs_lambda - Math.PI;
468             }
469             if (Math.Sin(A12) < 0 && Math.Tan(A21) >= 0)
470             {
471                 A21 = Abs_A21;
472             }
473             else if (Math.Sin(A12) < 0 && Math.Tan(A21) < 0)
474             {
475                 A21 = Math.PI - Abs_A21;
476             }
477             else if (Math.Sin(A12) >= 0 && Math.Tan(A21) >= 0)
478             {
479                 A21 = Math.PI + Abs_A21;
480             }
481             else if (Math.Sin(A12) >= 0 && Math.Tan(A21) < 0)
482             {
483                 A21 = 2 * Math.PI - Abs_A21;
484             }
485             //将球面元素换算到椭球面上
486             B2 = Math.Atan(Math.Sqrt(1 + m_ei2) * Math.Tan(u2));
487             double ki2 = m_e2 * Math.Cos(A0) * Math.Cos(A0);
488             double alphai = (m_e2 / 2 + m_e2 * m_e2 / 8 + m_e2 * m_e2 * m_e2 / 16) 
489                 - m_e2 * (1 + m_e2) * ki2 / 16 + 3 * m_e2 * ki2 * ki2 / 128;
490             double betai = m_e2 * (1 + m_e2) * ki2 / 16 - m_e2 * ki2 * ki2 / 32;
491             double gammai = m_e2 * ki2 * ki2 / 256;
492             double l = lambda - Math.Sin(A0) * (alphai * sigma + betai * Math.Sin(sigma) 
493                 * Math.Cos(2 * sigma1 + sigma) + gammai * Math.Sin(2 * sigma) * 
494                 Math.Cos(4 * sigma1 + 2 * sigma));
495             L2 = L1 + l;
496         }
497         //白塞尔大地主题反算函数 Inverse Solution of Bessel's Geodetic Problem
498         public void IS_BGP(double L1, double B1, double L2, double B2,
499             ref double S, ref double A12, ref double A21)
500         {
501             //将椭球面元素投影到球面上
502             double u1 = Math.Atan(Math.Sqrt(1 - m_e2) * Math.Tan(B1));
503             double u2 = Math.Atan(Math.Sqrt(1 - m_e2) * Math.Tan(B2));
504             double l = L2 - L1;
505             double a1 = Math.Sin(u1) * Math.Sin(u2);
506             double a2 = Math.Cos(u1) * Math.Cos(u2);
507             double b1 = Math.Cos(u1) * Math.Sin(u2);
508             double b2 = Math.Sin(u1) * Math.Cos(u2);
509             //迭代求lambda
510             double lambda = l;
511             double lambda0 = l;
512             double sigma = 0;
513             double sigma1 = 0;
514             double A0 = 0;
515             do
516             {
517                 lambda0 = lambda;
518                 double p = Math.Cos(u2) * Math.Sin(lambda);
519                 double q = b1 - b2 * Math.Cos(lambda);
520                 A12 = Math.Atan(p / q);
521                 double Abs_A12 = 0;
522                 Operation_Angle.ToFirstQuadrant(A12, ref Abs_A12);
523                 if (p >= 0 && q >= 0)
524                 {
525                     A12 = Abs_A12;
526                 }
527                 else if (p >= 0 && q < 0)
528                 {
529                     A12 = Math.PI - Abs_A12;
530                 }
531                 else if (p < 0 && q < 0)
532                 {
533                     A12 = Math.PI + Abs_A12;
534                 }
535                 else if (p < 0 && q >= 0)
536                 {
537                     A12 = 2 * Math.PI - Abs_A12;
538                 }
539                 double sins = Math.Cos(u2) * Math.Sin(lambda) * Math.Sin(A12) + (Math.Cos(u1) * Math.Sin(u2)
540                     - Math.Sin(u1) * Math.Cos(u2) * Math.Cos(lambda)) * Math.Cos(A12);
541                 double coss = Math.Sin(u1) * Math.Sin(u2) + Math.Cos(u1) * Math.Cos(u2) * Math.Cos(lambda);
542                 sigma = Math.Atan(sins / coss);
543                 double Abs_sigma = 0;
544                 Operation_Angle.ToFirstQuadrant(sigma, ref Abs_sigma);
545                 if (coss >= 0)
546                 {
547                     sigma = Abs_sigma;
548                 }
549                 else
550                 {
551                     sigma = Math.PI - Abs_sigma;
552                 }
553                 A0 = Math.Asin(Math.Cos(u1) * Math.Sin(A12));
554                 sigma1 = Math.Atan(Math.Tan(u1) / Math.Cos(A12));
555                 double ki2 = m_e2 * Math.Cos(A0) * Math.Cos(A0);
556                 double alphai = (m_e2 / 2 + m_e2 * m_e2 / 8 + m_e2 * m_e2 * m_e2 / 16)
557                     - m_e2 * (1 + m_e2) * ki2 / 16 + 3 * m_e2 * ki2 * ki2 / 128;
558                 double betai = m_e2 * (1 + m_e2) * ki2 / 16 - m_e2 * ki2 * ki2 / 32;
559                 double gammai = m_e2 * ki2 * ki2 / 256;
560                 lambda = l + Math.Sin(A0) * (alphai * sigma + betai * Math.Sin(sigma)
561                     * Math.Cos(2.0 * sigma1 + sigma) + gammai * Math.Sin(2.0 * sigma)
562                     * Math.Cos(4.0 * sigma1 + 2.0 * sigma));
563             }
564             while (Math.Abs(lambda - lambda0) >= PreciseControl.EPS);
565             //将球面元素换算到椭球面上
566             double k2 = m_ei2 * Math.Cos(A0) * Math.Cos(A0);
567             double A = 1 + k2 / 4 - 3 * k2 * k2 / 64 - 5 * k2 * k2 * k2 / 256;
568             double B = k2 / 4 - k2 * k2 / 16 + 15 * k2 * k2 * k2 / 512;
569             double C = k2 * k2 / 128 - 3 * k2 * k2 * k2 / 512;
570             double alpha = 1 / (m_b * A);
571             double beta = B / A;
572             double gamma = C / A;
573             S = (sigma - beta * Math.Sin(sigma) * Math.Cos(2 * sigma1 + sigma) - gamma
574                 * Math.Sin(2 * sigma) * Math.Cos(4 * sigma1 + 2 * sigma)) / alpha;
575             A21 = Math.Atan(Math.Cos(u1) * Math.Sin(lambda) / (Math.Cos(u1) * Math.Sin(u2)
576                 * Math.Cos(lambda) - Math.Sin(u1) * Math.Cos(u2)));
577             double Abs_A21 = 0;
578             Operation_Angle.ToFirstQuadrant(A21, ref Abs_A21);
579             if (Math.Sin(A12) < 0 && Math.Tan(A21) >= 0)
580             {
581                 A21 = Abs_A21;
582             }
583             else if (Math.Sin(A12) < 0 && Math.Tan(A21) < 0)
584             {
585                 A21 = Math.PI - Abs_A21;
586             }
587             else if (Math.Sin(A12) >= 0 && Math.Tan(A21) >= 0)
588             {
589                 A21 = Math.PI + Abs_A21;
590             }
591             else if (Math.Sin(A12) >= 0 && Math.Tan(A21) < 0)
592             {
593                 A21 = 2 * Math.PI - Abs_A21;
594             }
595         }
596         //Convert Measured Value of Distance to Ellipsoid 将测量值归算至椭球面
597         public double CMVTE_Distance(double D, double H1, double H2, double B1, double L1, double B2, double L2)
598         {
599             //此函数假设不知道两点连线的方位角，却知道两点的大地坐标，由大地坐标算出方位角。
600             double A12 = 0;
601             double S = 0;
602             double A21 = 0;
603             //由于距离一般较短，采用高斯平均引数法
604             GACIS_GP(L1, B1, L2, B2, ref S, ref A12, ref A21);
605             //计算起点曲率半径
606             double RA = Get_N(B1) / (1 + m_ei2 * Math.Pow(Math.Cos(B1) * Math.Cos(A12), 2));
607             //三个临时变量
608             double temp1 = 1 - Math.Pow((H2 - H1) / D, 2);
609             double temp2 = (1 + H1 / RA) * (1 + H2 / RA);
610             double temp3 = Math.Pow(D, 3) / (24 * RA * RA);
611             //距离换算公式
612             return D * Math.Sqrt(temp1 / temp2) + temp3;
613         }
614         public double CMVTE_Distance(double D, double H1, double H2, double B1,double A12)
615         {
616             //此函数假设知道方位角。
617 
618             //计算起点曲率半径
619             double RA = Get_N(B1) / (1 + m_ei2 * Math.Pow(Math.Cos(B1) * Math.Cos(A12), 2));
620             //三个临时变量
621             double temp1 = 1 - Math.Pow((H2 - H1) / D, 2);
622             double temp2 = (1 + H1 / RA) * (1 + H2 / RA);
623             double temp3 = Math.Pow(D, 3) / (24 * RA * RA);
624             //距离换算公式
625             return D * Math.Sqrt(temp1 / temp2) + temp3;
626         }
627         #endregion
628 

Code

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

1 abstract class MapProjection
2     {
3         protected ReferenceEllipsoid m_Ellipsoid;
4     }

Code
  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;
100             n = (int)(NCS_y / 1000000);
101             y = NCS_y - n * 1000000 - 500000;
102         }
103 
104         //距离改化公式，将椭球面是的距离化算至平面距离
105         public double Convert_Distance(double S, double L1, double B1, double L2, double B2)
106         {
107             int n = 0;
108             double x1 = 0;
109             double y1 = 0;
110             double x2 = 0;
111             double y2 = 0;
112             Positive_Computation(L1, B1, ref x1, ref y1, ref n);
113             Positive_Computation(L2, B2, ref x2, ref y2, ref n);
114             double R1 = m_Ellipsoid.Get_R(B1);
115             double R2 = m_Ellipsoid.Get_R(B2);
116             double Rm = (R1 + R2) / 2.0;
117             double ym = (y1 + y2) / 2.0;
118             double deta_y = y1 - y2;
119             return S * (1 + ym * ym / (2 * Rm * Rm) + deta_y * deta_y / (24 * Rm * Rm) + Math.Pow(ym / Rm, 4) / 24);
120         }
121         #endregion
122 

两个辅助类：

Code
 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;
10             }
11             if (minite < 0 || minite >= 60)
12             {
13                 return false;
14             }
15             if (second < 0 || second >= 60)
16             {
17                 return false;
18             }
19             double temp = degree + minite / 60.0 + second / 3600.0;
20             radian = Math.PI * temp / 180.0;
21             return true;
22         }
23         static public bool  AngleToRadian(double degree, ref double radian)
24         {
25             if (degree < 0 || degree >= 360)
26             {
27                 return false;
28             }
29             radian = Math.PI * degree / 180.0;
30             return true;
31         }
32         //弧度转角度
33         static public bool RadianToAngle(double radian, ref int degree, ref int minite, ref double second)
34         {
35             if (radian < 0 || radian >= 2 * Math.PI)
36             {
37                 return false;
38             }
39             double temp = 180.0 * radian / Math.PI;
40             degree = (int)temp;
41             temp = (temp - (double)degree) * 60;
42             minite = (int)temp;
43             temp = (temp - (double)minite) * 60;
44             second = temp;
45             return true;
46         }
47         //将一个角度通过加减90度化到第一象限
48         static public void ToFirstQuadrant(double radian, ref double f_radian)
49         {
50             f_radian = radian;
51             while (f_radian >= Math.PI / 2)
52             {
53                 f_radian -= Math.PI / 2;
54             }
55             while (f_radian < 0)
56             {
57                 f_radian += Math.PI / 2;
58             }
59         }
60         static public void ToFirstQuadrant(int degree, int minite, double second,
61             ref int f_degree, ref int f_minite, ref double f_second)
62         {
63             f_degree = degree;
64             f_minite = minite;
65             f_second = second;
66             while (f_degree >= 90)
67             {
68                 f_degree -= 90;
69             }
70             while (f_degree < 0)
71             {
72                 f_degree += 90;
73             }
74         }
75 

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