C++：高斯坐标，大地坐标转经纬度

实现CGCS2000大地坐标系XY值转化为对应经纬度信息，
注意输入的XY值与给定的值相反，则参数X为已知的Y，参数Y为已知的X。 得出的结果为 [纬度，经度]，不要应用错误。
L0参数为中央子午线的经线值，应用方法前需要确定L0的值，否则得出的值会有很大的偏差。

void xylb(double l0, double x, double y, double *l, double *b)
{
    double bf,vf,nf,ynf,tf,yf2,hbf;
    double sa,sb,se2,sep2,mf;
    double w1,w2,w,w3,w4;
    double pi = 3.1415926;
    x = x/1000000.0;
    y = y - 500000.0;
    bf = 9.04353692458*x-0.00001007623*pow(x,2.0)-0.00074438304*pow(x,3.0)-0.00000463064*pow(x,4.0)+0.00000505846*pow(x,5.0)-0.00000016754*pow(x,6.0);
    hbf = bf * pi/ 180.0;
    //北京
//    sa = 6378245.0;
//    sb = 6356863.019;
    //西安
//    sa = 6378140;
//    sb = 6356755.28820;
    //高斯
    sa = 6378137.0; //长半轴 m
    sb = 6356752.31414; //短半轴 m

    se2 = 0.006693421623;
    sep2 = 0.006738525415;
    w1 = sin(hbf);
    w2 = 1.0 - se2 * pow(w1,(double)2);
    w = sqrt(w2);
    mf = sa*(1.0-se2)/pow(w,(double)3);
    w3 = cos(hbf);
    w4 = pow(sa,(double)2)*pow(w3,(double)2) + pow(sb,(double)2)*pow(w1,(double)2);
    nf = pow(sa,(double)2) / sqrt(w4);
    ynf = y/nf;
    vf = nf/mf;
    tf = tan(hbf);
    yf2 = sep2 * pow(w3, (double)2);
    *b = bf - 1.0/2.0 * vf * tf * (pow(ynf,(double)2)-1.0/12.0*(5.0+3.0*pow(tf,(double)2)+yf2-9.0*yf2*pow(tf,(double)2))*pow(ynf,(double)4))*180.0/pi;
    *l = 1.0/w3*ynf*(1.0-1.0/6.0*(1.0+2.0*pow(tf,(double)2)+yf2)*pow(ynf,(double)2)+1.0/120.0*(5.0+28.0*pow(tf,(double)2)+24.0*pow(tf,(double)2)+6.0*yf2+8.0*yf2*pow(tf,(double)2))*pow(ynf,(double)4))*180.0/pi;
    *l = l0 + *l;
}