地心地固坐标系（ECEF）和WGS-84坐标系（WGS84）互转

//经纬度WGS-84坐标系（WGS84）转 地心地固坐标系（ECEF）
myPoint WGS84toECEF(double latitude, double longitude, double height)
{
	double X;
	double Y;
	double Z;
	double a = 6378137;
	double b = 6356752.314245;
	double E = (a * a - b * b) / (a * a);
	double COSLAT = cos(latitude * PI / 180);
	double SINLAT = sin(latitude * PI / 180);
	double COSLONG = cos(longitude * PI / 180);
	double SINLONG = sin(longitude * PI / 180);
	double N = a / (sqrt(1 - E * SINLAT * SINLAT));
	double NH = N + height;
	X = NH * COSLAT * COSLONG;
	Y = NH * COSLAT * SINLONG;
	Z = (b * b * N / (a * a) + height) * SINLAT;
	return myPoint(X, Y, Z);
}

//地心地固坐标系（ECEF）转 经纬度WGS-84坐标系（WGS84）
void ECEFtoWGS84(myPoint pt, double& lon, double& lat, double& altitude)
{
	double x = pt.X(), y = pt.Y(), z = pt.Z();
	double a, b, c, d;
	double p, q;
	double N;
	a = 6378137.0;
	b = 6356752.31424518;
	c = sqrt(((a * a) - (b * b)) / (a * a));
	d =	sqrt(((a * a) - (b * b)) / (b * b));
	p = sqrt((x * x) + (y * y));
	q = atan2((z * a), (p * b));
	lon = atan2(y, x);
	lat = atan2((z + (d * d) * b * pow(sin(q), 3)), (p - (c * c) * a * pow(cos(q), 3)));
	N = a / sqrt(1 - ((c * c) * pow(sin(lat), 2)));
	altitude = (p / cos(lat)) - N;
	lon = lon * 180.0 / PI;
	lat = lat * 180.0 / PI;
}