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]);
}
}