高斯克吕格与地理坐标相互转换算法（JS版本）

最近一段时间在研究高斯克吕格与地理坐标的互换算法，刚才的时候写了一个只能用于标准分带的算法，发现并不符合实际的一些地方坐标系的互换操作。经过研究最终写出了即可以应用于标准分带的和地方性的高斯克吕格与地理坐标系的算法，现在贴出来供大家参考，也希望大家批评指正。

代码如下：

/*
 * 投影变换算法。用于地方坐标系

 *MZ(f, 2)的MZ()方法为平方函数
 */
SuperMap.Web.Tool.ProjectionTransfor=function(TuoqiuCanshu, CentralMeridian,OriginLatitude,EastOffset,NorthOffset)
{

/* '说明: 用于初始化转换参数

'TuoqiuCanshu  枚举类型，提供北京54、西安80和WGS84三个椭球参数

'CentralMeridian 中央经线

OriginLatitude 原点纬度，如果是标准的分幅，则该参数是0

'EastOffset东偏移

NorthOffset 北偏移

*/

 ///'基本变量定义
 var a ;//'椭球体长半轴
 var b;// '椭球体短半轴
 var f; //'扁率
 var e;// '第一偏心率
 var e1; //'第二偏心率

 var FE ;//'东偏移
 var FN ;//'北偏移
 var L0 ;//'中央经度
 var W0;//'原点纬线
 var k0 ;//'比例因子
 

var PI = 3.14159265358979;

 /*
  *  Canshu
  *  Beijing54 = 0
     Xian80 = 1
     WGS84 = 2
  * 
 */
 /*'Krassovsky （北京54采用） 6378245 6356863.0188
 'IAG 75（西安80采用） 6378140 6356755.2882
 'WGS 84 6378137 6356752.3142
 */
 if(TuoqiuCanshu==0)//北京五四
 {
     a = 6378245;
     b = 6356863.0188;
 }
 else if(TuoqiuCanshu==1)// '西安八零
 {
     a = 6378140;
     b = 6356755.2882;
 }
 if(TuoqiuCanshu==2)//'WGS84
 {
     a = 6378137;
     b = 6356752.3142;
 }
 f = (a - b) / a;//扁率
 //e = Math.sqrt(1 - MZ((b / a) ,2));//'第一偏心率
 e = Math.sqrt(2*f - MZ(f ,2));//'第一偏心率
 
 //eq = Math.sqrt(MZ((a / b) , 2) - 1);//'第二偏心率
  e1 = e / Math.sqrt(1 - MZ(e , 2));//'第二偏心率

L0 = CentralMeridian;//中央经
W0= OriginLatitude;//原点纬线
 k0 = 1;//'比例因子
 FE = EastOffset;//东偏移
 FN = NorthOffset;//北偏移
 /*
  * 输入参数分别是：经度、纬度
  */
 this.JWgetGK=function( J,  W)
 {
 //'给出经纬度坐标，转换为高克投影坐标
 var resultP=new SuperMap.Web.Core.Point2D() ;
 var BR = (W - W0) * PI / 180;//纬度弧长
 var lo = (J - L0)*PI/180; //经差弧度
 
var N = a / Math.sqrt(1 - MZ((e * Math.sin(BR)) , 2)) //卯酉圈曲率半径

//求解参数s
var B0;
var B2;
var B4;
var B6;
var B8;
var C = MZ(a , 2)/ b;
B0 = 1 - 3 * MZ(e1 , 2) / 4 + 45 *MZ( e1 ,4) / 64 - 175 * MZ(e1 , 6) / 256 + 11025 * MZ(e1 , 8 )/ 16384;
B2 = B0 - 1
B4 = 15 / 32 * MZ(e1 , 4) - 175 / 384 * MZ(e1 , 6 )+ 3675 / 8192 *MZ( e1 , 8);
B6 = 0 - 35 / 96 *MZ( e1 , 6) + 735 / 2048 * MZ(e1 , 8);
B8 = 315 / 1024 * MZ(e1 , 8);
s = C * (B0 * BR + Math.sin(BR) * (B2 * Math.cos(BR) + B4 * MZ((Math.cos(BR)) , 3) + B6 * MZ((Math.cos(BR)) , 5 )+ B8 * MZ((Math.cos(BR)) , 7)))

var t = Math.tan(BR);
var g = e1 * Math.cos(BR);

 var XR= s + MZ(lo , 2) / 2 * N * Math.sin(BR) * Math.cos(BR) + MZ(lo , 4 )* N * Math.sin(BR) * MZ((Math.cos(BR)) , 3) / 24 * (5 -MZ( t , 2 )+ 9 * MZ(g , 2) + 4 *MZ( g , 4)) + MZ(lo , 6) * N * Math.sin(BR) * MZ((Math.cos(BR)) , 5) * (61 - 58 *MZ( t , 2) + MZ(t , 4)) / 720;
 var YR= lo * N * Math.cos(BR) + MZ(lo , 3 )* N / 6 *MZ( (Math.cos(BR)) , 3) * (1 - MZ(t , 2) + MZ(g , 2)) + MZ(lo , 5) * N / 120 * MZ((Math.cos(BR)) , 5) * (5 - 18 * MZ(t , 2) + MZ(t , 4) + 14 * MZ(g , 2) - 58 * MZ(g , 2) * MZ(t , 2));

resultP.x=YR+FE;
resultP.y=XR+FN;
 return resultP;
 }
 /*
  * 输入参数分别为：X、Y
  */
 this.GKgetJW=function( X,  Y)
{
 //'给出高克投影坐标，转换为经纬度坐标
var resultP=new SuperMap.Web.Core.Point2D();
var El1 = (1 - Math.sqrt(1 - MZ(e , 2))) / (1 + Math.sqrt(1 -MZ( e , 2)));

var Mf = (Y- FN)/ k0 ;//真实坐标值

var Q = Mf / (a * (1 - MZ(e , 2) / 4 - 3 * MZ(e , 4) / 64 - 5 *MZ( e , 6) / 256));//角度

Bf = Q + (3 * El1 / 2 - 27 *MZ( El1 , 3) / 32) * Math.sin(2 * Q) + (21 *MZ( El1 , 2) / 16 - 55 *MZ( El1 , 4 )/ 32) * Math.sin(4 * Q) + (151 *MZ( El1 , 3 )/ 96) * Math.sin(6 * Q) + 1097 / 512 * MZ(El1 , 4) * Math.sin(8 * Q);
Rf = a * (1 -MZ( e , 2)) / Math.sqrt(MZ((1 - MZ((e * Math.sin(Bf)) ,2)) , 3));
Nf = a / Math.sqrt(1 - MZ((e * Math.sin(Bf)) , 2));//卯酉圈曲率半径
Tf = MZ((Math.tan(Bf)) , 2);
D =(X - FE) / (k0 * Nf);

Cf =MZ( e1 , 2) * MZ((Math.cos(Bf)) , 2);

var B = Bf - Nf * Math.tan(Bf) / Rf * (MZ(D , 2) / 2 - (5 + 3 * Tf + 10 * Cf - 9 * Tf * Cf - 4 *MZ( Cf , 2) - 9 * MZ(e1 , 2)) *MZ( D , 4) / 24 + (61 + 90 * Tf + 45 * MZ(Tf , 2) - 256 * MZ(e1 , 2) - 3 * MZ(Cf , 2)) *MZ( D , 6) / 720);
var L = CentralMeridian*PI/180 + 1 / Math.cos(Bf) * (D - (1 + 2 * Tf + Cf) *MZ( D , 3) / 6 + (5 - 2 * Cf + 28 * Tf - 3 *MZ( Cf , 2) + 8 * MZ(e1 , 2) + 24 * MZ(Tf , 2)) * MZ(D , 5 )/ 120);

var Bangle = B * 180 / PI;
var Langle = L * 180 / PI;
 
 
 resultP.x = Langle;//RW * 180 / PI;
 resultP.y = Bangle + W0;//RJ * 180 / PI;
 return resultP;
}

};
SuperMap.Web.Tool.ProjectionTransfor.registerClass('SuperMap.Web.Tool.ProjectionTransfor');