当前位置:文档之家› 高斯正反算(零误差)

高斯正反算(零误差)

//84的椭球final double a = 6378137;final double Alfa = 1.0 / 298.257223563;double centreL, x, y, b, e1, ee;double a0, a2, a4, a6, a8, Bf0;double[] Coeficient_a0 = new double[5];double sinBf, cosBf;double FBf, Bf1, dB, bf;double c, v, Nf, Mf, tf;double itaf, dietaB, dietaL;double B, L;double dmsB, dmsL, dmsCentreL1;double radlat, radlon, radl0, l;double sb, cb, t, ita, X, N;public String MakeProject(double L, double B, double CentreLon) //高斯正算{/*输入已知数据:经度\纬度\ 中央子午线*/dmsB = B;dmsL = L;dmsCentreL1 = CentreLon;radlat = DMSTORAD(dmsB);radlon = DMSTORAD(dmsL);radl0 = DMSTORAD(dmsCentreL1);l = radlon - radl0;b = a * (1 - Alfa);sb = Math.sin(radlat);cb = Math.cos(radlat);t = sb / cb;e1 = Math.sqrt((a / b) * (a / b) - 1);ee = Math.sqrt(1 - (b / a) * (b / a));ita = e1 * cb;a0a2a4a6a8(a, ee, Coeficient_a0);a0 = Coeficient_a0[0];a2 = Coeficient_a0[1];a6 = Coeficient_a0[3];a8 = Coeficient_a0[4];X = a0 * radlat - sb * cb * ((a2 - a4 + a6) + (2 * a4 - 16 * a6 / 3) * sb * sb + 16 * a6 * sb * sb * sb * sb / 3.0);//计算卯酉圈半径Nc = a * a / b;v = Math.sqrt(1 + e1 * e1 * cb * cb);N = c / v;//计算未知点的坐标x = X + N * sb * cb * l * l / 2 + N * sb * cb * cb * cb * (5 - t * t + 9 * ita * ita + 4 * ita * ita * ita * ita) * l * l * l * l / 24 + N * sb * cb * cb * cb * cb * cb * (61 - 58 * t * t + t * t * t * t) * l * l * l * l * l * l / 720;y = N * cb * l + N * cb * cb * cb * (1 - t * t + ita * ita) * l * l * l / 6 + N * cb * cb * cb * cb * cb * (5 - 18 * t * t + t * t * t * t + 14 * ita * ita - 58 * ita * ita * t * t) * l * l * l * l * l / 120;//输出未知点坐标return (y + 500000)+","+x;}public Point MoveProject(double CentreLon, double X1, double Y1){centreL = CentreLon * Math.PI / 180.0;x = X1;while(Y1>1000000)Y1=Y1-1000000;y = Y1 - 500000.0;b = a * (1 - Alfa);e1 = Math.sqrt((a / b) * (a / b) - 1);//第二偏心率ee = Math.sqrt(1 - (b / a) * (b / a));a0a2a4a6a8(a, ee, Coeficient_a0);a0 = Coeficient_a0[0];a2 = Coeficient_a0[1];a4 = Coeficient_a0[2];a8 = Coeficient_a0[4];Bf0 = x / a0;do{sinBf = Math.sin(Bf0); cosBf = Math.cos(Bf0);FBf = -sinBf * cosBf * ((a2 - a4 + a6) + (2.0 * a4 - 16.0 * a6 / 3.0) * sinBf * sinBf +(16.0 / 3.0) * a6 * (sinBf * sinBf * sinBf * sinBf));Bf1 = (x - FBf) / a0;dB = Bf1 - Bf0;Bf0 = Bf1;} while (Math.abs(dB * 180.0 / Math.PI * 3600.0) > 0.00001);bf = Bf1;c = a * a / b;v = Math.sqrt(1 + e1 * e1 * Math.cos(bf) * Math.cos(bf));Nf = c / v;Mf = c / (v * v * v);tf = Math.sin(bf) / Math.cos(bf);itaf = e1 * Math.cos(bf);dietaB = tf * Math.pow(y, 2) / (2 * Mf * Nf) - tf * (5 + 3 * Math.pow(tf, 2) + Math.pow(itaf, 2) - 9 * Math.pow(tf, 2) * Math.pow(itaf, 2)) * Math.pow(y, 4) / (24 * Mf * Math.pow(Nf, 3)) + (61 + 90 * Math.pow(tf, 2) + 45 * Math.pow(tf, 4)) * Math.pow(y, 6) / (720 * Mf * Math.pow(Nf, 5));dietaL = y / (Nf * Math.cos(bf) + (1 + 2 * tf * tf + itaf * itaf) * Math.cos(bf) * y * y / (6 * Nf)) + (5 + 44 * tf * tf + 32 * tf * tf * tf * tf - 2 * itaf * itaf - 16 * itaf * itaf * tf * tf) / (360 * Nf * Nf * Nf * Mf * Mf * Math.cos(bf));B = bf - dietaB;L = centreL + dietaL;dmsB = RADTODMS(B);dmsL = RADTODMS(L);return new Point(dmsL , dmsB);}public void a0a2a4a6a8(double a, double e, double[] Coeficient_a0){double m0, m2, m4, m6, m8;m0 = a * (1 - e * e);m2 = 3 * e * e * m0 / 2;m4 = 5 * e * e * m2 / 4;m6 = 7 * e * e * m4 / 6;m8 = 9 * e * e * m6 / 8;Coeficient_a0[0] = m0 + m2 / 2 + 3 * m4 / 8 + 5 * m6 / 16 + 35 * m8 / 128;Coeficient_a0[1] = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;Coeficient_a0[2] = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;Coeficient_a0[3] = m6 / 32 + m8 / 16;Coeficient_a0[4] = m8 / 128;}static int intSignOfDms, intSignOfRad;static double radAngle, dmsAngle;public static double DMSTORAD(double dmsAngle1){intSignOfDms = 1;if (dmsAngle1 < 0) intSignOfDms = -1;dmsAngle1 = Math.abs(dmsAngle1);radAngle = dmsAngle1 * Math.PI / 180.0;radAngle = radAngle * intSignOfDms;return radAngle;}public static double RADTODMS(double radAngle){intSignOfRad = 1;if (radAngle < 0) intSignOfRad = -1;radAngle = Math.abs(radAngle);dmsAngle = radAngle * 180 / Math.PI;dmsAngle = dmsAngle * intSignOfRad;return dmsAngle;}。

相关主题