WGS84坐标转BJ54坐标

int main(int argc, char* argv[])
{
 //WGS84坐标
 double dLongitude = 117;
 double dLatitude =40;

 double EARTH_WGS84_A = 6378137.0000;
 double EARTH_WGS84_E2 = 0.00669437999013;
 double EARTH_WGS84_FLATTENING = 298.257223563;

 double EARTH_BJ54_A = 6378245.0;
 double EARTH_BJ54_E2 = 0.00667;
 double EARTH_BJ54_FLATTENING = 298.3;  

 double OMIGA = 206264.8062;
 double DELTA_X = -28.3;
 double DELTA_Y = 144.9;
 double DELTA_Z = 77.5;

 double sinRenda = sin(dLongitude);
    double cosRenda = cos(dLongitude);
    double sinFi = sin(dLatitude);
    double cosFi = cos(dLatitude);
    double sinFi2 = sinFi*sinFi;
    double eSinFi = sqrt(1-EARTH_WGS84_E2*sinFi2);
    double M = EARTH_WGS84_A*(1-EARTH_WGS84_E2)/(eSinFi*eSinFi*eSinFi);
    double N = EARTH_WGS84_A/eSinFi;
    double deltaE2 = 2*(1-1/EARTH_BJ54_FLATTENING)*(1/EARTH_WGS84_FLATTENING-1/EARTH_BJ54_FLATTENING);
    double deltaA = EARTH_BJ54_A - EARTH_WGS84_A;
    double deltaRenda = OMIGA*(cosRenda*DELTA_Y - sinRenda*DELTA_X)/(N*sinRenda);
    double deltaFi = OMIGA*((EARTH_WGS84_A*deltaE2 + EARTH_WGS84_E2*deltaA)*sinFi*cosFi + EARTH_WGS84_A*EARTH_WGS84_E2*deltaE2*sinFi*sinFi*sinFi*cosFi-sinFi*cosRenda*DELTA_X - sinFi*sinRenda*DELTA_X + cosFi*DELTA_Z)/M;

 //BJ54坐标
    double dLatitude54 = dLatitude + deltaFi/3600;
    double dLongitude54 = dLongitude + deltaRenda/3600;

 return 0;
}