原文：http://www.cnblogs.com/SuperMagic/archive/2007/12/04/982817.html
/* 功能說明： 將絕對高斯坐標(y,x)轉換成絕對的地理坐標(wd,jd)。 */
// double y; 輸入參數: 高斯坐標的橫坐標，以米為單位
// double x; 輸入參數: 高斯坐標的縱坐標，以米為單位
// short DH; 輸入參數: 帶號，表示上述高斯坐標是哪個帶的
// double *L; 輸出參數: 指向經度坐標的指針，其中經度坐標以秒為單位
// double *B; 輸出參數: 指向緯度坐標的指針，其中緯度坐標以秒為單位
void GaussToGeo(double y, double x, short DH, double *L, double *B, double LP)
{
double l0; // 經差
double tf; // tf = tg(Bf0),注意要將Bf轉換成以弧度為單位
double nf ; // n = y * sqrt( 1 + etf ** 2) / c, 其中etf = e'**2 * cos(Bf0) ** 2
double t_l0; // l0，經差，以度為單位
double t_B0; // B0，緯度，以度為單位
double Bf0; // Bf0
double etf; // etf,其中etf = e'**2 * cos(Bf0) ** 2
double X_3 ;
double PI=3.14159265358979;
double b_e2=0.0067385254147;
double b_c=6399698.90178271;
X_3 = x / 1000000.00 - 3 ; // 以兆米（1000000）為單位
// 對於克拉索夫斯基橢球，計算Bf0
Bf0 = 27.11115372595 + 9.02468257083 * X_3 - 0.00579740442 * pow(X_3,2)
- 0.00043532572 * pow(X_3,3) + 0.00004857285 * pow(X_3,4)
+ 0.00000215727 * pow(X_3,5) - 0.00000019399 * pow(X_3,6) ;
tf = tan(Bf0*PI/180); // tf = tg(Bf),注意這里將Bf轉換成以弧度為單位
etf = b_e2 * pow(cos(Bf0*PI/180),2); // etf = e'**2 * cos(Bf) ** 2
nf = y * sqrt( 1 + etf ) / b_c; // n = y * sqrt( 1 + etf ** 2) / c
// 計算緯度，注意這里計算出來的結果是以度為單位的
t_B0 = Bf0 - (1.0+etf) * tf / PI * (90.0 * pow(nf,2)
- 7.5 * (5.0 + 3 * pow(tf,2) + etf - 9 * etf * pow(tf,2)) * pow(nf,4)
+ 0.25 * (61 + 90 * pow(tf,2) + 45 * pow(tf,4)) * pow(nf,6)) ;
// 計算經差，注意這里計算出來的結果是以度為單位的
t_l0 = (180 * nf - 30 * ( 1 + 2 * pow(tf,2) + etf ) * pow(nf,3)
+ 1.5 * (5 + 28 * pow(tf,2) + 24 * pow(tf,4)) * pow(nf,5))
/ ( PI * cos(Bf0*PI/180) ) ;
l0 = (t_l0 * 3600.0); // 將經差轉成秒
if (LP == -1000)
{
*L = (double)((DH * 6 - 3) * 3600.0 + l0); // 根據帶號計算出以秒為單位的絕對經度，返回指針
}
else
{
*L = LP * 3600.0 + l0; // 根據帶號計算出以秒為單位的絕對經度，返回指針
}
//----------------------------------
*B = (double)(t_B0 * 3600.0) ; // 將緯差轉成秒，並返回指針
}
/* 功能說明： （1）將地理坐標(wd,jd)轉換成絕對的高斯坐標(y,x)
（2）本函數支持基於六度帶（或三度帶）、克拉索夫斯基橢球進行轉換 */
/* 適用范圍： 本函數適用於將地球東半球中北半球（即東經0度到東經180度，北緯0度至90度）范圍
內所有地理坐標到高斯坐標的轉換 */
/* 使用說明： 調用本函數后返回的結果應在滿足精度的條件下進行四舍五入 */
// double jd; 輸入參數: 地理坐標的經度，以秒為單位
// double wd; 輸入參數: 地理坐標的緯度，以秒為單位
// short DH; 輸入參數: 三度帶或六度帶的帶號
/* 六度帶(三度帶)的帶號是這樣得到的：從東經0度到東經180度自西向東按每6度(3度)順序編號
(編號從1開始)，這個順序編號就稱為六度帶(三度帶)的帶號。因此，六度帶的帶號的范圍是1-30，
三度帶的帶號的范圍是1-60。
如果一個點在圖號為TH的圖幅中，那麽該點所處的六度帶的帶號就可以這樣得到：將該圖號的
第3、4位組成的字符串先轉換成數字，再減去30。例如某點在圖幅06490701中，該點所在的帶號就
是49-30，即19。
如果調用本函數去進行一般的從地理坐標到基於六度帶高斯坐標的變換（非鄰帶轉換），則參
數DH的選取按前一段的方法去確定。
如果調用本函數去進行基於六度帶鄰帶轉換，則參數DH的選取先按上述方法去確定，然后看是
往前一個帶還是后一個帶進行鄰帶轉換再確定是加1還是減1。 */
void GeoToGauss(double jd, double wd, short DH, short DH_width, double *y, double *x, double LP)
{
double t; // t=tgB
double L; // 中央經線的經度
double l0; // 經差
double jd_hd,wd_hd; // 將jd、wd轉換成以弧度為單位
double et2; // et2 = (e' ** 2) * (cosB ** 2)
double N; // N = C / sqrt(1 + et2)
double X; // 克拉索夫斯基橢球中子午弧長
double m; // m = cosB * PI/180 * l0
double tsin,tcos; // sinB,cosB
double PI=3.14159265358979;
double b_e2=0.0067385254147;
double b_c=6399698.90178271;
jd_hd = jd / 3600.0 * PI / 180.0 ; // 將以秒為單位的經度轉換成弧度
wd_hd = wd / 3600.0 * PI / 180.0 ; // 將以秒為單位的緯度轉換成弧度
// 如果不設中央經線（缺省參數: -1000），則計算中央經線，
// 否則，使用傳入的中央經線，不再使用帶號和帶寬參數
//L = (DH - 0.5) * DH_width ; // 計算中央經線的經度
if (LP == -1000)
{
L = (DH - 0.5) * DH_width ; // 計算中央經線的經度
}
else
{
L = LP ;
}
l0 = jd / 3600.0 - L ; // 計算經差
tsin = sin(wd_hd); // 計算sinB
tcos = cos(wd_hd); // 計算cosB
// 計算克拉索夫斯基橢球中子午弧長X
X = 111134.8611 / 3600.0 * wd - (32005.7799 * tsin + 133.9238 * pow(tsin,3)
+ 0.6976 * pow(tsin,5) + 0.0039 * pow(tsin,7) ) * tcos;
et2 = b_e2 * pow(tcos,2) ; // et2 = (e' ** 2) * (cosB ** 2)
N = b_c / sqrt( 1 + et2 ) ; // N = C / sqrt(1 + et2)
t = tan(wd_hd); // t=tgB
m = PI/180 * l0 * tcos; // m = cosB * PI/180 * l0
*x = X + N * t * ( 0.5 * pow(m,2)
+ (5.0 - pow(t,2) + 9.0 * et2 + 4 * pow(et2,2)) * pow(m,4)/24.0
+ (61.0 - 58.0 * pow(t,2) + pow(t,4)) * pow(m,6) / 720.0 ) ;
*y = N * ( m + ( 1.0 - pow(t,2) + et2 ) * pow(m,3) / 6.0
+ ( 5.0 - 18.0 * pow(t,2) + pow(t,4) + 14.0 * et2
- 58.0 * et2 * pow(t,2) ) * pow(m,5) / 120.0 );
}