由于没有时间，目前只实现了大地主题解算、高斯投影坐标正反算和距离改正，但应该很容易扩充。
结构示意图：
参考椭球类：
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_e2106 {107 get108 {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_ei2123 {124 get125 {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 #endregion140 #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 Problem212 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 do222 {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 * tm234 * 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 else252 {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 * tm274 + 2 * Math.Pow(etam, 2.0)) / (24 * Vm * Vm);275 double s30 = Nm * Math.Cos(Bm) * Math.Cos(Bm) * (2 + 3 * tm * tm276 + 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 + 2281 * 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 else293 {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 else324 {325 A21 = A21 - Math.PI;326 }327 }328 329 //高斯平均引数反算中由Am、Bm、S计算参数alpha330 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计算参数beta336 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计算参数gamma343 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和v349 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 Problem355 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 do367 {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 else381 {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 else413 {414 A21 = A21 - Math.PI;415 }416 }417 //白塞尔大地主题正算函数 Direct Solution of Bessel's Geodetic Problem418 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 do435 {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 Problem498 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 //迭代求lambda510 double lambda = l;511 double lambda0 = l;512 double sigma = 0;513 double sigma1 = 0;514 double A0 = 0;515 do516 {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 else550 {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) - gamma574 * 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 #endregion628
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 #endregion122
两个辅助类：
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 }
转载于:https://www.cnblogs.com/waynecraig/archive/2009/09/26/1574306.html
相关资源：【椭球大地测量学】MATLAB实现贝塞尔大地问题正反解计算编程（含流程图）
