java代码实现CGCS2000大地坐标系经纬度转换为对应XY值

通过经纬度转换为XY值

输入的参数为经度，纬度

得到的结果是[X，Y]

代码内部longitude0代表*子午线的值，需要根据实际内容配置，也可以自己提到参数中

​​​​​​​public static double[] GaussToBLToGauss(double longitude, double latitude)//参数 经度，纬度
{
    double[] output = new double[2];
    double longitude1,latitude1, longitude0, X0,Y0, xval,yval;
    //NN曲率半径，测量学里面用N表示
    //M为子午线弧长，测量学里用大X表示
    //fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示
    //R为底点所对的曲率半径，测量学里用Nf表示
    double a,f, e2,ee, NN, T,C,A, M, iPI;
    iPI = 0.0174532925199433; //3.1415926535898/180.0;
    a=6378137.0; f=1/298.257222101; //CGCS2000坐标系参数
    //a=6378137.0; f=1/298.2572236; //wgs84坐标系参数
    longitude0 = 117.0;//*子午线 根据实际进行配置
    longitude0 = longitude0 * iPI ;//*子午线转换为弧度
    longitude1 = longitude * iPI ; //经度转换为弧度
    latitude1 = latitude * iPI ; //纬度转换为弧度
    e2=2*f-f*f;
    ee=e2*(1.0-e2);
    NN=a/Math.sqrt(1.0-e2*Math.sin(latitude1)*Math.sin(latitude1));
    T=Math.tan(latitude1)*Math.tan(latitude1);
    C=ee*Math.cos(latitude1)*Math.cos(latitude1);
    A=(longitude1-longitude0)*Math.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)*Math.sin(2*latitude1)
            +(15*e2*e2/256+45*e2*e2*e2/1024)*Math.sin(4*latitude1)-(35*e2*e2*e2/3072)*Math.sin(6*latitude1));
    xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120);
    yval = M+NN*Math.tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
            +(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);
    X0 = 500000L;
    Y0 = 0;
    xval = xval+X0; yval = yval+Y0;

    //转换为投影
    output[0] = xval;
    output[1] = yval;
    return output;
}