结构示意图:
![](https://i-blog.csdnimg.cn/blog_migrate/0264d84d8fe9bbdba8f1528c9f4e9e2a.jpeg)
参考椭球类:
![ContractedBlock.gif](https://www.cnblogs.com/Images/OutliningIndicators/ContractedBlock.gif)
![ExpandedBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
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
![ContractedBlock.gif](https://www.cnblogs.com/Images/OutliningIndicators/ContractedBlock.gif)
坐标投影类和高斯坐标投影类:
1
abstract
class
MapProjection
2 {
3 protected ReferenceEllipsoid m_Ellipsoid;
4 }
2 {
3 protected ReferenceEllipsoid m_Ellipsoid;
4 }
![ContractedBlock.gif](https://www.cnblogs.com/Images/OutliningIndicators/ContractedBlock.gif)
![ExpandedBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
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
两个辅助类:
![ContractedBlock.gif](https://www.cnblogs.com/Images/OutliningIndicators/ContractedBlock.gif)
![ExpandedBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
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
![ContractedBlock.gif](https://www.cnblogs.com/Images/OutliningIndicators/ContractedBlock.gif)
![ExpandedBlockStart.gif](https://www.cnblogs.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
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 }