最近一段時間在研究高斯克呂格與地理座標的互換演算法，剛才的時候寫了一個只能用於標準分帶的演算法，發現並不符合實際的一些地方座標系的互換操作。經過研究最終寫出了即可以應用於標準分帶的和地方性的高斯克呂格與地理座標系的演算法，現在貼出來供大家參考，也希望大家批評指正。

程式碼如下：

/*
 * 投影變換演算法。用於地方座標系
 */
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');