经纬度转化平面坐标

#define MAXFIELD 64
#define GGALENGTH 12
#define RMCLENGTH 12
#define VTGLENGTH 12
#define PIP 3.1415926535897932 /* pi /
#define rad (PIP/180.0) / deg to rad /
#define deg (180.0/PIP) / rad to deg /
#define WGS84_a 6378137.0
#define WGS84_f 1.0/298.257223563
#define WGS84_ee 0.00669437999013
void GaussCal( double Latitude, double Longitude, double X, double Y)
{
double a, f, ee, et;
double Latitude0, Longitude0, Latitude1, Longitude1,X0,Y0;
double N, T, C, A, M;
int BeltNumber = 0;
int Beltwidth = 6;
BeltNumber = (int)(Longitude / Beltwidth);
Longitude0 = (BeltNumberBeltwidth + Beltwidth / 2)rad;
Latitude0 = 0rad;
Latitude1 = Latituderad;
Longitude1 = Longituderad;
X0 = 500000 + (BeltNumber + 1) * 1000000;
Y0 = 0;
a = WGS84_a;
f = WGS84_f;
ee = 2 * f - ff;
et = ee * (1.0 - ee);
N = a / sqrt(1.0-eesin(Latitude1)sin(Latitude1));
T = tan(Latitude1)tan(Latitude1);
C = etcos(Latitude1)cos(Latitude1);
A = (Longitude1 - Longitude0)cos(Latitude1);
M = a((1 - ee / 4 - 3 * eeee / 64 - 5 * eeeeee / 256)Latitude1 - (3 * ee / 8 + 3 * eeee / 32 + 45 * eeee*ee / 1024)*sin(2 * Latitude1)

• (15 * eeee / 256 + 45 * eeeeee / 1024)sin(4 * Latitude1) - (35 * eeeeee / 3072)sin(6 * Latitude1));
  X = X0 + N(A + (1 - T + C)AAA / 6 + (5 - 18 * T + TT + 72 * C - 58 * et)pow(A, 5) / 120);
  Y = Y0+M + Ntan(Latitude1)(AA / 2 + (5 - T + 9 * C + 4 * CC)AAA*A / 24
• (61 - 58 * T + T*T + 600 * C - 330 * et)*pow(A, 6) / 720);
  }