Gaussian projection positive and inverse code

 

// Gaussian projection positive and Inverse Calculation // 6-degree bandwidth 54-year Beijing coordinate system // Gaussian projection by latitude and longitude (unit: dd) returns the geocoordinate (with a number, unit: metres) void gaussprojcal (double longpolling, double latitude, double * X, double * Y) {int projno = 0; int zonewide; /// bandwidth double longitude1, latitude1, longitude0, latitude0, x0, y0, xval, yval; Double A, F, E2, EE, NN, T, C,, m, IPI; IPI = 0.0174532925199433; // 3.1415926535898/180.0; zonewide = 6; // 6 degree bandwidth a = 6378245.0; F = 1.0/298.3; // Beijing coordinate system parameters for 54 years // A = 6378140.0; F = 1/298. 257; // Xi'an coordinate system parameter projno = (INT) (longpolling/zonewide); longitude0 = projno * zonewide + zonewide/2; longitude0 = longitude0 * IPI; latitude0 = 0; longitude1 = longitude * IPI; // converts longitude to radian latitude1 = latitude * IPI; // converts latitude to radian e2 = 2 * f-f * F; EE = e2 * (1.0-e2); nn = A/SQRT (1.0-e2 * sin (latitude1) * sin (latitude1); t = tan (latitude1) * Tan (latitude1); C = EE * Cos (latitude1) * Cos (latitude1); A = (longitude1-longitude0) * Cos (latitude1 ); M = A * (1-e2/4-3 * e2 * E2/64-5 * e2 * e2 * E2/256) * latitude1-(3 * E2/8 + 3 * e2 * E2/32 + 45 x e2 * e2 * E2/1024) * sin (2 * latitude1) + (15 * e2 * E2/256 + 45 * e2 * e2 * E2/1024) * sin (4 * latitude1)-(35 * e2 * e2 * E2/3072) * sin (6 * l atitude1); xval = nn * (a + (1-T + C) * A * a * A/6 + (5-18 * t + T x T + 72 * C-58 * ee) * A * a * A/120 ); yval = m + nn * Tan (latitude1) * (A * A/2 + (5-T + 9 * C + 4 * C) * A * a * A/24 + (61-58 * t + T * t + 600 * C-330 * ee) * A * a * A/720); X0 = 000000l * (projno + 1) + 500000l; Y0 = 0; xval = xval + x0; yval = yval + y0; * x = xval; * Y = yval;} // Gaussian projection returns the latitude and longitude (unit: dd) by the coordinate (unit: metres) void gaussprojinvcal (Double X, Double Y, double * longdistance, double * latitude) string 9 {int projno; int zonewide; // bandwidth double longitude1, latitude1, longitude0, latitude0, x0, y0, xval, yval; double E1, E2, F, A, EE, NN, T, C, M, D, R, U, Fai, IPI; IPI = 0.0174532925199433; // 3.1415926535898/180.0; A = 6378245.0; F = 1.0/298.3; // Beijing coordinate system parameter for 54 years /// A = 6378140.0; F = 1/298. 257; // Xi'an coordinate system parameter zonewide = 6 in 80 years; // projno = (INT) (X/000000l) in 6-degree bandwidth ); // find the belt number longitude0 = (ProjNo-1) * zonewide + zonewide/2; longitude0 = longitude0 * IPI; // The Central longitude line X0 = projno * 000000l + 500000l; Y0 = 0; xval = X-X0; yval = Y-Y0; // in-band geocoordinate e2 = 2 * f-f * F; e1 = (1.0-sqrt (1-e2 )) /(1.0 + SQRT (1-e2); EE = e2/(1-e2); M = yval; U = m/(A * (1-E2/4-3 * e2 * E2/64-5 * e2 * e2 * E2/256 )); fai = u + (3 * E1/2-27 * E1 * E1 * E1/32) * sin (2 * U) + (21 * E1 * E1/16-55 * E1 * E1 * E1 * E1/32) * sin (4 * U) + (151 * E1 * E1 * E1/96) * sin (6 * u) + (1097 * E1 * E1 * E1 * E1 * E1/512) * sin (8 * U); C = EE * Cos (FAI); t = tan (FAI) * Tan (FAI ); nn = A/SQRT (1.0-e2 * sin (FAI); string 1 R = A * (1-e2)/SQRT (1-e2 * sin (FAI) * sin (FAI) * (1-e2 * sin (FAI) * (1-e2 * sin (FAI ))); D = xval/nn; // calculate the longitude (longlongitude) latitude (latitude) longitude1 = longitude0 + (D-(1 + 2 * t + C) * D/6 + (5-2 * C + 28 * T-3 * C + 8 * EE + 24 * T * t) * D/120)/cos (FAI); latitude1 = Fai-(NN * Tan (FAI)/R) * (D * D/2-(5 + 3 * t + 10 * C-4 * C * C-9 * ee) * D/24 + (61 + 90 * t + 298 * C + 45 * T * T-256 * ee-3 * C) * D/720); // convert to degree dd * longpolling = longitude1/IPI; * latitude = latitude1/IPI ;} // NN (curvature radius), which is represented by N in measurement. // M indicates the meridian arc length, and X indicates the string 2 in measurement. // Fai indicates the base point latitude, obtained by the Inverse Calculation Formula of the Meridian arc length. BF is used in measurement to represent the curvature radius of the string 4 // R as the base point, and NF is used in measurement.