Gaussian projection positive and Inverse Calculation

// Gaussian projection positive and Inverse Calculation
///// Bandwidth of 6 degrees 54 years Beijing Coordinate System
// Gaussian projection inverse returns the geocoordinate (including band number, unit: metres) from the latitude and longitude (unit: dd)
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, A, M, IPI;
IPI = 0.0174532925199433; // 3.1415926535898/180.0;
Zonewide = 6; // bandwidth of 6 degrees
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 parameters for 80 years
Projno = (INT) (longpolling/zonewide );
Longitude0 = projno * zonewide + zonewide/2;
Longitude0 = longitude0 * IPI;
Latitude0 = 0;
Longitude1 = longitude * IPI; // converts the longitude to a radian.
Latitude1 = latitude * IPI; // converts a latitude to a 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 * 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 * 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 inverse latitude and longitude (unit: dd) by the coordinate (unit: metres)
Void gaussprojinvcal (Double X, Double Y, double * longpolling, double * latitude)

{
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 parameters for 54 years
/// A = 6378140.0; F = 1/298. 257; // Xi'an coordinate system parameters for 80 years
Zonewide = 6; // bandwidth of 6 degrees
Projno = (INT) (X/000000l); // search for a carry ID
Longitude0 = (ProjNo-1) * zonewide + zonewide/2;
Longitude0 = longitude0 * IPI; // central longitude
X0 = projno * 000000l + 500000l;
Y0 = 0;
Xval = X-X0; yval = Y-Y0; // in-band Earth Coordinate
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 ));
R = A * (1-e2)/SQRT (1-e2 * sin (FAI) * (1-e2 * sin (FAI )) * (1-e2 * sin
(FAI) * sin (FAI )));
D = xval/nn;
// Calculate the longitude (latitude)
Longitude1 = longitude0 + (D-(1 + 2 * t + C) * D/6 + (5-2 * C + 28 * T-3 * C + 8 * EE + 24 * T * t) * d
* 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 returns the radius of curvature in the ring, which is represented by N in measurement.

M is the meridian arc length, which is expressed by X in measurement.

FAI is the latitude of the base point, which is obtained by the Inverse Calculation Formula of the Meridian arc length. BF is used in the measurement.

R is the curvature radius of the base point, expressed by NF in measurement.