/** * Space earth Cartesian coordinates-> Earth coordinates */ Public geopoint xyz_blh (INT ellipse, point3d point) { Geopoint = new geopoint (); If (geopoint = NULL | point = NULL) { System. Out. println ("the object is blank! "); } Double A = 0, B = 0, e1 = 0, E2 = 0; // A is a long half axis, B is a short half axis, E1 is the first eccentric heart rate, E2 is the square of the first eccentric Heart Rate Double DETA = 0.0000001; Switch (ellipse) // select an elliptical body { Case 0: // WGS84 elliptical body { A = 6378137.0; // long Half Axis B = 6356752.3142; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = 0.006694384999588; Break; } Case 1: // Beijing 54 elliptical { A = 6378245.0; // long Half Axis B = 6356863.0187730473; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = 0.006693421622966; Break; } Case 2: // Xi'an 80 elliptical { A = 6378140.0; // long Half Axis B = 6356755.2881575287; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = 0.006694384999588; Break; } Default: break; } // Double W = math. SQRT (1.0-E1 * E1 * Math. Pow (math. Sin (point. getx (), 2 )); // Double N = A/W; Double X = point. getx (); Double Y = point. Gety (); Double Z = point. Getz (); Double M = math. SQRT (math. Pow (x, 2) + math. Pow (Y, 2 )); Double L = geopoint. getl (); // geopoint longitude Double B = geopoint. getb (); // Double H = geopoint. geth (); // large height L = math. atan (y/X ); If (L <0) { L + = math. Pi; } Double E2 _ = e2/(1-e2 ); Double C = A * Math. SQRT (1 + E2 _); Double CE2 = C * E2; Double K = 1 + E2 _; Double front = z/m; Double temp = front; Int COUNT = 0; // number of iterations Do { Front = temp; M = math. SQRT (math. Pow (x, 2) + math. Pow (Y, 2 )); Temp = z/m + CE2 * front/(m * Math. SQRT (K + math. Pow (front, 2 ))); } While (math. Abs (front-Temp) <DETA & count <100000); // whether it is within the error range B = math. atan (temp); // obtain the latitude. If (B <0) { B + = math. Pi; } Double W = math. SQRT (1-e1 * E1 * Math. Sin (B) * Math. Sin (B )); Double N = A/W; // N = (A * m-C * C)/(2 * B * z ); System. Out. println ("n =" + n ); H = m/Math. Cos (B)-N; // height // H = z/Math. Sin (B)-N * (1.0-E1 * E1 ); // H = x/(math. Cos (B) * Math. Cos (L)-N; // Convert to the angle value L = math. todegrees (L ); B = math. todegrees (B ); // H = math. todegrees (h ); Geopoint. SETl (L ); Geopoint. SETB (B ); Geopoint. Seth (h ); Return geopoint; } /** * Coordinate-> space ry Cartesian coordinates * @ Param ellipse elliptical body * @ Param geopoint coordinate before conversion * @ Return returns the Cartesian coordinates of the space. */ Public point3d blh_xyz (INT ellipse, geopoint) { Point3d point1 = new point3d (1, 1 ); If (geopoint = NULL | point1 = NULL) { System. Out. println ("the object is blank! "); } Double A = 0, B, e1 = 0, E2 = 0; // A is a long half axis, B is a short half axis, E1 is the first eccentric heart rate, E2 is the square of the first eccentric Heart Rate Switch (ellipse) // select an elliptical body { Case 0: // WGS84 elliptical body { A = 6378137.0; // long Half Axis B = 6356752.3142; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = 0.006694384999588; Break; } Case 1: // Beijing 54 elliptical { A = 6378245.0; // long Half Axis B = 6356863.0187730473; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = 0.006693421622966; Break; } Case 2: // Xi'an 80 elliptical { A = 6378140.0; // long Half Axis B = 6356755.2881575287; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = 0.006694384999588; Break; } Default: break; } Double B = geopoint. getb () * Math. PI/180; // converts to radians. Double L = geopoint. getl () * Math. PI/180; Double H = geopoint. geth (); // Calculate the radius of curvature of a uniform Ring Double N = A/Math. SQRT (1.0-E1 * E1 * Math. Pow (math. Sin (B), 2 )); // Calculate space Cartesian coordinates Double X = (n + H) * Math. Cos (B) * Math. Cos (L ); Double Y = (n + H) * Math. Cos (B) * Math. Sin (L ); Double Z = (N * (1-e1 * E1) + H) * Math. Sin (B ); Point1.x = X; Point1.y = y; Point1.z = z; Return point1; } /** * Gaussian-kerluge coordinate normal calculation (from the earth coordinate to the plane Cartesian coordinate) * @ Param ellipse elliptical body * @ Param zonewide bandwidth (3 or 6 degrees) * @ Param geopoint ry * @ Param point Cartesian coordinates */ Public void bl_xy (INT ellipse, int zonewide, geopoint, point) { If (geopoint = NULL | point = NULL) { System. Out. println ("the object is blank! "); } Double A = 0, B, e1 = 0, E2 = 0; // A is a long half axis, B is a short half axis, E1 is the first eccentric heart rate, E2 is the square of the first eccentric Heart Rate Switch (ellipse) // select an elliptical body { Case 0: // WGS84 elliptical body { A = 6378137.0; // long Half Axis B = 6356752.3142; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B; // E2e = 0.006694384999588; Break; } Case 1: // Beijing 54 elliptical { A = 6378245.0; // long Half Axis B = 6356863.0187730473; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B; // E2e = 0.006693421622966; Break; } Case 2: // Xi'an 80 elliptical { A = 6378140.0; // long Half Axis B = 6356755.2881575287; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B; // E2e = 0.006694384999588; Break; } Default: break; }
Double B = geopoint. getb (); Double L = geopoint. getl (); Double centerl = 0.0; // degree of projection with central longitude Int n = 0; // The belt number of the projection band (6 ° band: 1-60, 3 ° band: 1-120) Double L = 0; // The Difference Between the longitude of the region and the central longitude If (6 = zonewide) // if it is a 6-degree band { N = (INT) L/6; If (L % 6> 0) { N = n + 1; } Centerl = 6 * n-3; L = L-centerl; } If (3 = zonewide) // if it is a 3-degree band { N = (INT) (L-1.5)/3; If (L-1.5) % 3> 0) { N = n + 1; } Centerl = 3 * N; L = L-centerl; } // Convert to radians B = B * Math. PI/180; L = L * Math. PI/180; L = L * Math. PI/180; Double X; // the meridian arc length calculated from the equator // Calculate the coefficient of meridian Arc Length Double a0 = 1.0 + 3.0/4.0 * Math. Pow (E1, 2) + 45.0/64.0 * Math. Pow (E1, 4) + 350.0/512.0 * Math. Pow (E1, 6) + 11025.0/16384.0 * Math. Pow (E1, 8 ); Double a2 =-3.0/4 * Math. Pow (E1, 2) + 60.0/64 * Math. Pow (E1, 4) + 525.0/512 * Math. Pow (E1, 6) + 17640.0/16384.0 * Math. Pow (E1, 8)/2.0; Double A4 = 15.0/64.0 * Math. Pow (E1, 4) + 210.0/512.0 * Math. Pow (E1, 6) + 8820.0/16384.0 * Math. Pow (E1, 8)/4.0; Double A6 =-35.0/512.0 * Math. Pow (E1, 6) + 2520.0/16384.0 * Math. Pow (E1, 8)/6.0; Double A8 = 315.0/16384.0 * Math. Pow (E1, 8)/8.0; // Calculate the meridian arc length x X = A * (1.0-math. pow (E1, 2) * a0 * B + A2 * Math. sin (2 * B) + A4 * Math. sin (4 * B) + A6 * Math. Sin (6 * B) + A8 * Math. Sin (8 * B ); Double T = math. Tan (B ); Double Anke = e2 * Math. Cos (B ); Double N = A/Math. SQRT (1.0-math. Pow (E1, 2) * Math. Pow (math. Sin (B), 2 )); // Coordinate Calculation Double X = point. getx (); Double Y = point. Gety (); X = x + 1.0/2.0 * n * T * Math. Pow (math. Cos (B), 2) * Math. Pow (L, 2) + 1.0/24.0 * n * T * (5.0-math. pow (T, 2) + 9 * Math. pow (Anke, 2) + 4 * Math. pow (Anke, 4) * Math. pow (math. cos (B), 4) * Math. pow (L, 4) + 1.0/720 * n * T * (61.0-58 * Math. pow (T, 2) + math. pow (T, 4) + 270.0 * Math. pow (Anke, 2)-330.0 * Math. pow (Anke, 2) * Math. pow (T, 2) * Math. pow (math. cos (B), 6 );
Y = N * Math. cos (B) * L + 1.0/6 * n * (1.0-math. pow (T, 2) + math. pow (Anke, 2) * Math. pow (math. cos (B), 3) * Math. pow (L, 3) + 1.0/120 * n * (5.0-18 * T * t + math. pow (T, 4) + 14 * Math. pow (Anke, 2)-58 * Anke * T * t )* Math. Pow (math. Cos (B), 3) * Math. Pow (L, 5 ); Y + = 500000.0; // Add 500 km Point. setx (X ); Point. sety (y ); } /** * The Cartesian coordinate transformation of the Gaussian plane is the earth coordinate. * @ Param ellipse elliptical body * @ Param zonewide bandwidth * @ Param point plane coordinate * @ Param geopoint ry */ Public void xy_bl (INT ellipse, double centerl, point, geopoint) { If (geopoint = NULL | point = NULL) { System. Out. println ("the object is blank! "); } Double A = 0, B, e1 = 0, E2 = 0; // A is a long half axis, B is a short half axis, E1 is the first eccentric heart rate, E2 is the second eccentric Heart Rate Switch (ellipse) // select an elliptical body { Case 0: // WGS84 elliptical body { A = 6378137.0; // long Half Axis B = 6356752.3142; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B; // E2e = 0.006694384999588; Break; } Case 1: // Beijing 54 elliptical { A = 6378245.0; // long Half Axis B = 6356863.0187730473; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B; // E2e = 0.006693421622966; Break; } Case 2: // Xi'an 80 elliptical { A = 6378140.0; // long Half Axis B = 6356755.2881575287; // short half axis E1 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/; E2 = math. SQRT (math. Pow (A, 2)-math. Pow (B, 2)/B; // E2e = 0.006694384999588; Break; } Default: break; } Double X = point. getx (); Double Y = point. Gety ()-500000; Double BF = 0.0; // latitude of the base point Double b0 = 1.0 + 3.0/4.0 * E1 * E1 + 45.0/64.0 * Math. Pow (E1, 4) + 350.0/312.0 * Math. Pow (E1, 6) + 11025.0/16384.0 * Math. Pow (E1, 8) + 43659.0/65536.0 * Math. Pow (E1, 10 ); BF = x/(A * (1.0-math. Pow (E1, 2) * B0 ); Double TF = math. Tan (BF ); Double mf = A * (1.0-E1 * E1)/Math. SQRT (math. Pow (1.0-math. Pow (E1, 2), 3 )); Double NF = (A/math. SQRT (1.0-math. pow (E1, 2)/math. SQRT (1.0 + math. pow (e2 * Math. cos (BF), 2 )); Double Anke = e2 * Math. Cos (BF ); Double B = geopoint. getb (); Double L = geopoint. getl (); // Start Coordinate Calculation B = BF-TF/(2 * MF * NF * Math. cos (BF) * Math. pow (Y, 2) + TF/(24 * MF * Math. pow (NF, 3 ))* (5.0 + 3.0 * tF * TF + Anke * anke-9.0 * Anke * tF) * Math. Pow (Y, 4 )- 1.0/(720.0 * Math. pow (NF, 5) * Math. cos (BF) * (61.0 + 90.0 * tF * TF + 45 * Math. pow (TF, 4) * Math. pow (Y, 6 ); L = y/(NF * Math. cos (BF)-(1.0 + 2 * tF * TF + Anke * Anke) * Math. pow (Y, 3)/(6.0 * Math. pow (NF, 3) * Math. cos (BF )) + (5.0 + 28 * tF * TF + 24 * Math. pow (TF, 4) + 6 * Anke + 8 * Anke * tF) * Math. pow (Y, 5)/(120.0 * Math. pow (NF, 5) * Math. cos (BF ));
B = B * 180/Math. Pi; L = L * 180/Math. Pi; L = L + centerl; // central longitude + longitude difference Geopoint. SETB (B ); Geopoint. SETl (L ); } The possible reason for these codes is that, when the space Cartesian coordinates are converted to the Earth coordinates, the calculation of the Earth height is very different. |