主要算法修改自网络上一个VB 脚本实现版。已经通过Google 地图 API提供的方法进行了验证。具体计算的经纬度存在一定的误差,可以容忍。
/**
* ** a, b, c, alpha, e, e2, W, V As Double'a 长轴半径
'b 短轴
'c 极曲率半径
'alpha 扁率
'e 第一偏心率
'e2 第二偏心率
'W 第一基本纬度函数
'V 第二基本纬度函数
** B1, L1, B2, L2 As Double
'B1 点1的纬度
'L1 点1的经度
'B2 点1的纬度
'L2 点2的经度
** S As Double '''''大地线长度
** A1, A2 As Double
'A1 点1到点2的方位角
'A2 点2到点1的方位角
*/
private final static double pi=Math.PI;
private final static double M_RAD = 1.74532925199433E-02;
private double a,b,c,alpha,e,e2,w,V;
private double B1, L1, B2, L2;
private double s;
private double A1,A2;
private final double R = 6371229; //地球的半径
/**
* 椭球体 长半轴 a(米) 短半轴b(米)
Krassovsky (北京54采用)6378245 6356863.0188
IAG 75(西安80采用) 6378140 6356755.2882
WGS 84 6378137 6356752.3142
*/
a = 6378140;//6378245;//6378137
b = 6356755.2882;
c = Math.pow(a,2)/b;
alpha = (a-b)/a;
e = Math.sqrt(Math.pow(a,2)-Math.pow(b,2))/a;
e2 = Math.sqrt(Math.pow(a,2)-Math.pow(b,2))/b;
}
private void computerw(){
w = Math.sqrt(1 - Math.pow(e,2)*(Math.pow(Math.sin(B1),2)));