经纬度转换屏面坐标

package com.hhwy.webgis.map.util;
import java.util.*;
public class ProjectTranslate {

 
 
 
 
 /*-------------------计算坐标反算因数-----------------*/
    //fdao输入--扁率倒数
    //a输入--椭球半径，单位米
    //K0,K1,K2,K3为反算因数
    private void ComputeBLFactor(double a, double fdao,  double k[])
    {
        double e;
        double f;
        double b;

        double A1;
        double B1;
        double C1;
        double D1;
        double E1;
        double F1;

        double fa;
        double fb;
        double fc;
        double fd;
        double fe;
        double fg;

        double fc1, fc2, fc3, fc4, fb1, fb2, fb3, fb4, fd1, fd2, fd3, fd4;

        f = 1 / fdao;
        b = a - a * f;
        e = Math.sqrt(a * a - b * b) / a;

        A1 = 1 + Math.pow(e, 2) * 3 / 4 + Math.pow(e, 4) * 45 / 64 + Math.pow(e, 6) * 175 / 256 + Math.pow(e, 8) * 11025 / 16384 + Math.pow(e, 10) * 43659 / 65536;

        B1 = Math.pow(e, 2) * 3 / 4 + Math.pow(e, 4) * 15 / 16 + Math.pow(e, 6) * 525 / 512 + Math.pow(e, 8) * 2205 / 2048 + Math.pow(e, 10) * 72765 / 65536;

        C1 = Math.pow(e, 4) * 15 / 64 + Math.pow(e, 6) * 105 / 256 + Math.pow(e, 8) * 2205 / 4096 + Math.pow(e, 10) * 10395 / 16384;

        D1 = Math.pow(e, 6) * 35 / 512 + Math.pow(e, 8) * 315 / 2048 + Math.pow(e, 10) * 31185 / 13072;

        E1 = Math.pow(e, 8) * 315 / 16384 + Math.pow(e, 10) * 3465 / 65536;

        F1 = Math.pow(e, 10) * 693 / 13072;

        fa = A1 * a * (1 - e * e);
        fb = -B1 * a * (1 - e * e) / 2;
        fc = C1 * a * (1 - e * e) / 4;
        fd = -D1 * a * (1 - e * e) / 6;
        fe = E1 * a * (1 - e * e) / 8;
        fg = -F1 * a * (1 - e * e) / 10;

        fb1 = -fb / fa; fc1 = -fc / fa; fd1 = -fd / fa;

        fb2 = fb1 + fb1 * fc1 - 3.0 / 2.0 * fb1 * Math.pow(fb1, 2) - 2 * fc1 * fb1;
        fc2 = fc1 + fb1 * fb1;
        fd2 = fd1 + fb1 * fc1 + 1.0 / 2.0 * fb1 * Math.pow(fb1, 2) + 2 * fc1 * fb1;

        fb3 = fb1 + fb1 * fc2 - 3.0 / 2.0 * fb1 * Math.pow(fb2, 2) - 2 * fc1 * fb2;
        fc3 = fc1 + fb1 * fb2;
        fd3 = fd1 + fb1 * fc2 + 1.0 / 2.0 * fb1 * Math.pow(fb2, 2) + 2 * fc1 * fb2;

        fb4 = fb1 + fb1 * fc3 - 3.0 / 2.0 * fb1 * Math.pow(fb3, 2) - 2 * fc1 * fb3;
        fc4 = fc1 + fb1 * fb3;
        fd4 = fd1 + fb1 * fc3 + 1.0 / 2.0 * fb1 * Math.pow(fb3, 2) + 2 * fc1 * fb3;

        k[0] = 1.0 / fa;
        k[1] = (fb4 * 2 + fc4 * 4 + fd4 * 6);
        k[2] = (fc4 * 8 + fd4 * 32);
        k[3] = fd4 * 32;
    }
   
   
    /*-------------------计算坐标正算因数-----------------*/
    //fdao输入--扁率倒数
    //a输入--椭球半径，单位米
    //C0,C1,C2,C3为反算因数
    private void ComputeFactor(double a, double fdao,  double C[])
    {
        double e;
        double f;
        double b;

        double A11;
        double B11;
        double C11;
        double D11;
        double E11;
        double F11;

        double fa;
        double fb;
        double fc;
        double fd;
        double fe;
        double fg;

        f = 1 / fdao;
        b = a - a * f;
        e = Math.sqrt(a * a - b * b) / a;
        A11 = 1 + Math.pow(e, 2) * 3 / 4 + Math.pow(e, 4) * 45 / 64 + Math.pow(e, 6) * 175 / 256 + Math.pow(e, 8) * 11025 / 16384 + Math.pow(e, 10) * 43659 / 65536;

        B11 = Math.pow(e, 2) * 3 / 4 + Math.pow(e, 4) * 15 / 16 + Math.pow(e, 6) * 525 / 512 + Math.pow(e, 8) * 2205 / 2048 + Math.pow(e, 10) * 72765 / 65536;

        C11 = Math.pow(e, 4) * 15 / 64 + Math.pow(e, 6) * 105 / 256 + Math.pow(e, 8) * 2205 / 4096 + Math.pow(e, 10) * 10395 / 16384;

        D11 = Math.pow(e, 6) * 35 / 512 + Math.pow(e, 8) * 315 / 2048 + Math.pow(e, 10) * 31185 / 13072;

        E11 = Math.pow(e, 8) * 315 / 16384 + Math.pow(e, 10) * 3465 / 65536;

        F11 = Math.pow(e, 10) * 693 / 13072;

        fa = A11 * a * (1 - e * e);
        fb = -B11 * a * (1 - e * e) / 2;
        fc = C11 * a * (1 - e * e) / 4;
        fd = -D11 * a * (1 - e * e) / 6;
        fe = E11 * a * (1 - e * e) / 8;
        fg = -F11 * a * (1 - e * e) / 10;

        C[0] = fa;
        C[1] = -(fb * 2 + fc * 4 + fd * 6);
        C[2] = (fc * 8 + fd * 32);
        C[3] = -fd * 32;
        /**********************************************
        54 coordinate system computation result have difference from
        the known value;but I don't find fault.
        known: CO 6367558.49686
               C1 32005.79642
               C2 133.86115
               C3  0.7031
        compute:
               CO 6367558.4968746
               C1 32005.780305529
               C2 133.9203503
               c3 0.7041
    *****************************************************************/
    }
    /*-----------------------------------------------------------*/
    //lon输入--经度，lat输入--纬度，单位度
    //x输出--平面坐标x值，y输出--平面坐标y值，单位米
    //L0输入--中央子午线，单位度
    //fdao输入--扁率倒数
    //a输入--椭球半径，单位米
    public void GuassBLtoXY(double lon,double lat,  double coordXY[], float L0, double fdao, float a)
    {
        double t, Ita2, N, m0, l;
        double Temp1, Temp2, Temp3, Temp4, Temp5, Temp6, Temp7, Temp8;
        double Pi = Math.PI;
        double f;
        double b;
        double e2;
        double e12;
        double p2 = 3600.0 * 180.0 / Pi;
        double P0 = Pi / 180.0;
        /*  54 Coordonation     KeLaSuoFuSiJi */
        // double C0= 6367558.49686;
        // double C1=32005.79642;
        // double C2=133.86115;
        // double C3=0.7031;
        double C0, C1, C2, C3;

        C0 = C1 = C2 = C3 = 0;
        double []c = new double[4] ;
       
        ComputeFactor(a, fdao,c);
        C0=c[0] ;
        C1=c[1] ;
        C2=c[2] ;
        C3=c[3] ;
        double temp1, temp2, temp3, temp4, temp5;

        f = 1 / fdao;
        b = a - a * f;
        e2 = (a * a - b * b) / (a * a);
        e12 = e2 / (1 - e2);
        /* GausTransform */
        l = (lon - L0) * 3600;
        t = Math.tan(lat * P0);
        temp1 = t * t;

        Ita2 = e12 * Math.cos(lat * P0) * Math.cos(lat * P0);
        temp2 = Ita2 * Ita2;

        N = a / Math.sqrt(1 - e2 * Math.sin(lat * P0) * Math.sin(lat * P0));
        m0 = l * Math.cos(lat * P0) / p2;
        temp3 = m0 * m0;
        temp4 = temp3 * temp3;

        Temp1 = N * m0;
        Temp2 = C0 * lat * P0;

        temp5 = Math.sin(lat * P0) * Math.sin(lat * P0);
        Temp3 = Math.cos(lat * P0) * Math.sin(lat * P0) * (C1 + C2 * temp5
                + C3 * temp5 * temp5);
        Temp4 = 1.0 / 2.0 * N * t * temp3;
        Temp5 = 1 / 24.0 * (5.0 - (temp1 * temp1) + 9 * Ita2 + 4 * (temp2 * temp2)) * N * t * temp4;
        Temp6 = 1 / 720 * (61 - 58 * temp1 + (temp1 * temp1)) * N * t * (temp3 * temp4);
        Temp7 = 1 / 6.0 * (1 - temp1 + Ita2) * N * (m0 * temp3);
        Temp8 = 1 / 120.0 * (5 - 18 * temp1 + (temp1 * temp1) + 14 * temp2 - 58 * Ita2 * temp1) * N * (m0 * temp4);

        coordXY[0] = Temp2 - Temp3 + Temp4 + Temp5 + Temp6;
        coordXY[1] = Temp1 + Temp7 + Temp8 + 500000;

    }
    /*-----------------------------------------------------------*/
    //lon输入--经度，lat输入--纬度，单位度
    //x输出--平面坐标x值，y输出--平面坐标y值，单位米
    //L0输入--中央子午线，单位度
    //fdao输入--扁率倒数
    //a输入--椭球半径，单位米

    public void GuassXYtoBL( double[]coordBL , double x, double y, double L0, double fdao, double a)
    {
        double K0, K1, K2, K3;

        double f, b;

        double e2, e12;
        double Pi = Math.PI;
        K0 = K1 = K2 = K3 = 0;
        double[] k=new double[4] ;
        ComputeBLFactor(a, fdao, k);
        K0=k[0];
        K1=k[1];
        K2=k[2];
        K3=k[3];
        y -= 500000.00;
        f = 1 / fdao;
        b = a - a * f;
        e2 = (a * a - b * b) / (a * a);
        e12 = e2 / (1 - e2);

        double p0 = 180.0 / Pi;
        double p2 = 3600.0 * 180.0 / Pi;
        double P0 = Pi / 180.0;

        double B0f = K0 * x * p0;
        double SinB0f = Math.sin(B0f * P0);
        double SinB0f2 = Math.sin(B0f * P0) * Math.sin(B0f * P0);
        double Temp1 = p0 * (Math.cos(B0f * P0)) * SinB0f * (K1 - K2 * SinB0f2 + K3 * SinB0f2 * SinB0f2);
        double Bf = B0f + Temp1;

        double t = Math.tan(Bf * P0);
        double t2 = t * t;
        double Ita2 = e12 * Math.cos(Bf * P0) * Math.cos(Bf * P0);
        double V2 = 1.0 + Ita2;
        double N = a / Math.sqrt(1.0 - e2 * Math.sin(Bf * P0) * Math.sin(Bf * P0));
        double yN = y / N;
        double yN2 = yN * yN;

        double Temp2 = p2 * (-0.5) * V2 * t * yN2;
        double Temp3 = p2 * 1.0 / 24.0 * (5.0 + 3.0 * t2 + Ita2 - 9.0 * Ita2 * t2) * V2 * t * yN2 * yN2;
        double Temp4 = p2 * (-1.0) / 720.0 * (61.0 + 90.0 * t2 + 45.0 * t2 * t2) * V2 * yN2 * yN2 * yN2;

        double Temp5 = 1.0 / Math.cos(Bf * P0) * yN * p2;
        double Temp6 = p2 * (-1.0) / 6.0 * (1.0 + 2.0 * t2 + Ita2) * (1.0 / Math.cos(Bf * P0)) * yN * yN2;
        double Temp7 = p2 * 1.0 / 120.0 * (5.0 + 28.0 * t2 + 24.0 * t2 * t2 + 6.0 * Ita2 + 8.0 * Ita2 * t2)
                 * (1.0 / Math.cos(Bf)) * yN2 * yN2 * yN;

        coordBL[0] = Bf + (Temp2 + Temp3 + Temp4) / 3600.0;
        coordBL[1] = L0 + (Temp5 + Temp6 + Temp7) / 3600.0;

    }

    public static void main(String args[])
    {
     ProjectTranslate pt = new ProjectTranslate() ;
     double[] coordXY= new double[2] ;
     
     pt.GuassBLtoXY(117.1, 39, coordXY, 117, 298.3,6378245) ;
     System.out.println("====coordx="+coordXY[0]+"====coordy="+coordXY[1]);
     
    }
   
   
}