大地测量学基础算法实现（转）
大地测量学基础算法实现（转）
原文http://www.cnblogs.com/waynecraig/archive/2009/09/26/1574306.html
由于没有时间，目前只实现了大地主题解算、高斯投影坐标正反算和距离改正，但应该很容易扩充。
结构示意图：
参考椭球类：
View Code
//参考椭球类，可选择克拉索夫斯基椭球体、1975年国际椭球体、WGS-84椭球体，也可以自定义椭球参数。 class ReferenceEllipsoid { #region 椭球参数 //克拉索夫斯基椭球体参数 private const double a_1=6378245.0000000000; private const double b_1=6356863.0187730473; private const double c_1=6399698.9017827110; private const double r_alpha_1=298.3; private const double e2_1=0.006693421622966; private const double ei2_1=0.006738525414683; //1975年国际椭球体参数 private const double a_2=6378140.0000000000; private const double b_2=6356755.2881575287; private const double c_2=6399596.6519880105; private const double r_alpha_2=298.257; private const double e2_2=0.006694384999588; private const double ei2_2=0.006739501819473; //WGS-84椭球体参数 private const double a_3=6378137.0000000000; private const double b_3=6356752.3142; private const double c_3=6399593.6258; private const double r_alpha_3=298.257223563; private const double e2_3=0.0066943799013; private const double ei2_3=0.00673949674227; //自定义椭球体参数 private double a_4 = 0; private double b_4 = 0; private double c_4 = 0; private double r_alpha_4 = 0; private double e2_4 = 0; private double ei2_4 = 0; #endregion #region 成员变量 //标志椭球类型，1,2,3,4分别表示克拉索夫斯基椭球体、1975年国际椭球体、WGS-84椭球体、自定义椭球体 private int m_type = 0; public double m_a { get { switch (m_type) { case 1: return a_1; case 2: return a_2; case 3: return a_3; default: return a_4; } } } public double m_b { get { switch (m_type) { case 1: return b_1; case 2: return b_2; case 3: return b_3; default: return b_4; } } } public double m_c { get { switch (m_type) { case 1: return c_1; case 2: return c_2; case 3: return c_3; default: return c_4; } } } public double m_r_alpha { get { switch (m_type) { case 1: return r_alpha_1; case 2: return r_alpha_2; case 3: 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; }
View Code
class GaussProjection : MapProjection { #region 成员变量 private double m_D; //每带的宽度（弧度） private double m_Starting_L; //起点经度 #endregion #region 接口函数 public GaussProjection(ReferenceEllipsoid RE, double D, double Starting_L) { m_Ellipsoid = RE; m_D = D; m_Starting_L = Starting_L; } //坐标正算 public void Positive_Computation(double L, double B, ref double x, ref double y, ref int n) { //求带号 n = Convert.ToInt32((L - m_Starting_L) / m_D); //求中央经度 double L0 = m_Starting_L + (Convert.ToDouble(n) - 0.5) * m_D; //求点与中央子午线的经度差 double l = L - L0; //辅助变量 double m = Math.Cos(B) * l; double X = m_Ellipsoid.Get_X(B); double t = m_Ellipsoid.Get_t(B); double eta2 = Math.Pow(m_Ellipsoid.Get_eta(B), 2.0); double N = m_Ellipsoid.Get_N(B); //坐标正算公式 x = X + N * t * ((0.5 + ((5 - t * t + 9 * eta2 + 4 * eta2 * eta2) / 24 + (61 - 58 * t * t + Math.Pow(t, 4)) * m * m / 720) * m * m) * m * m); y = N * ((1 + ((1 - t * t + eta2) / 6 + (5 - 18 * t * t + Math.Pow(t, 4) + 14 * eta2 - 58 * eta2 * t * t) * m * m / 120) * m * m) * m); } //坐标反算 public void Negative_Computation(int n, double x, double y, ref double L, ref double B) { //此函数采用迭代算法 //辅助变量 double a = m_Ellipsoid.m_a; double ei2 = m_Ellipsoid.m_ei2; double e2 = m_Ellipsoid.m_e2; double c = m_Ellipsoid.m_c; double beta0 = 1 - 3 * ei2 / 4 + 45 * ei2 * ei2 / 64 - 175 * Math.Pow(ei2, 3) / 256 + 11025 * Math.Pow(ei2, 4) / 16384; double beta2 = beta0 - 1; double beta4 = 15 * ei2 * ei2 / 32 - 175 * Math.Pow(ei2, 3) / 384 + 3675 * Math.Pow(ei2, 4) / 8192; double beta6 = -35 * Math.Pow(ei2, 3) / 96 + 735 * Math.Pow(ei2, 4) / 2048; double beta8 = 315 * Math.Pow(ei2, 4) / 1024; //赋初值 double B0 = x / (c * beta0); double a10 =a * Math.Cos(B0) / Math.Sqrt(1 - e2 * Math.Sin(B0) * Math.Sin(B0)); double l0 = y / a10; double B1 = B0; double a11 = a10; double l1 = l0; do { B0 = B1; a10 = a11; l0 = l1; double N = m_Ellipsoid.Get_N(B0); double t = m_Ellipsoid.Get_t(B0); double eta2 = Math.Pow(m_Ellipsoid.Get_eta(B0), 2.0); double a2 = N * Math.Sin(B0) * Math.Cos(B0) / 2; double a4 = N * Math.Sin(B0) * Math.Pow(Math.Cos(B0), 3) * (5 - t * t + 9 * eta2 + 4 * eta2 * eta2) / 24; double a6 = N * Math.Sin(B0) * Math.Pow(Math.Cos(B0), 5) * (61 - 58 * t * t + Math.Pow(t, 4)) / 720; double FxB = (c * beta2 + (c * beta4 + (c * beta6 + c * beta8 * Math.Pow(Math.Cos(B0), 2)) * Math.Pow(Math.Cos(B0), 2)) * Math.Pow(Math.Cos(B0), 2)) * Math.Sin(B0) * Math.Cos(B0); double FxBl = a2 * l0 * l0 + a4 * Math.Pow(l0, 4) + a6 * Math.Pow(l0, 6); B1 = (x - FxB - FxBl) / (c * beta0); a11 = a * Math.Cos(B1) / Math.Sqrt(1 - e2 * Math.Sin(B1) * Math.Sin(B1)); double N1 = m_Ellipsoid.Get_N(B1); double t1 = m_Ellipsoid.Get_t(B1); double eta21 = Math.Pow(m_Ellipsoid.Get_eta(B1), 2.0); double a3 = N1 * Math.Pow(Math.Cos(B1), 3) * (1 - t1 * t1 + eta21) / 6; double a5 = N1 * Math.Pow(Math.Cos(B1), 5) * (5 - 18 * t1 * t1 + Math.Pow(t1, 4) + 14 * eta21 - 58 * eta21 * t1 * t1) / 120; double FyBl = a3 * Math.Pow(l0, 3) + a5 * Math.Pow(l0, 5); l1 = (y - FyBl) / a11; } while (Math.Abs(B1 - B0) > PreciseControl.EPS || Math.Abs(l1 - l0) > PreciseControl.EPS); double L0 = m_Starting_L + (n - 0.5) * m_D; L = L0 + l0; B = B0; } //对于6度带，将坐标转换成国家统一坐标形式 National Coordinate System static public void To_NCS(int n, double x, double y, ref double NCS_x, ref double NCS_y) { NCS_x = x; NCS_y = n * 1000000 + y + 500000; } static public void From_NCS(double NCS_x, double NCS_y, ref int n, ref double x, ref double y) { 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
两个辅助类：
View Code
//角度运算类 class Operation_Angle { //角度转弧度 static public bool AngleToRadian(int degree, int minite, double second, ref double radian) { if (degree < 0 || degree >= 360) { 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; } }
View Code
public class PreciseControl { //程序迭代精度控制 public static double EPS { get { return Math.Pow(10.0, -10.0); } } }
浙公网安备 33010602011771号