經緯度坐标與高斯坐标的轉換代碼

// 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) ;     //  将緯差轉成秒，并傳回指針

}

// double jd;         輸入參數: 地理坐标的經度，以秒為機關

// double wd;         輸入參數: 地理坐标的緯度，以秒為機關

// short  DH;      輸入參數: 三度帶或六度帶的帶号

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

}