UTM投影正轉算法實作

    在Gissky上發現了戴勤奮老師關于坐标轉換的程式以及公式，程式的界面和轉換精度都值得稱贊，隻是可惜沒有提供實作的源碼，最近由于工作需要正好需要實作這部分功能，通過對戴勤奮老師提供的公式以及OGP（歐洲石油勘探組織）最新的《Coordinate Conversions and Transformations including Formulas》文檔進行了研究，利用程式将UTM投影正轉實作，現釋出出來供大家共同學習，也希望大家指出程式中的不足。

LLtoUTM( double  a, double  f, const   double  Lat,  const   double  Long,

           double  LongOrigin, double  FN, double   & UTMNorthing, double   & UTMEasting)

... {

      //e表示WGS84第一偏心率,eSquare表示e的平方,

     double eSquare =(2*f-f*f) ;

     double k0 = 0.9996;

     double e2Square;

     double V, T, C, A, M;

      //確定longtitude位于-180.00----179.9之間

     double LongTemp = (Long+180)-int((Long+180)/360)*360-180; 

     double LatRad = Lat*PI/180;

     double LongRad = LongTemp*PI/180;

     double LongOriginRad;

     LongOriginRad = LongOrigin * PI/180;

     e2Square = (eSquare)/(1-eSquare);

     V = a/sqrt(1-eSquare*sin(LatRad)*sin(LatRad));

     T = tan(LatRad)*tan(LatRad);

     C = e2Square*cos(LatRad)*cos(LatRad);

     A = cos(LatRad)*(LongRad-LongOriginRad);

     M = a*((1-eSquare/4-3*eSquare*eSquare/64-5*eSquare*eSquare*eSquare/256)*LatRad

         -(3*eSquare/8+3*eSquare*eSquare/32+45*eSquare*eSquare*eSquare/1024)*sin(2*LatRad)

         +(15*eSquare*eSquare/256+45*eSquare*eSquare*eSquare/1024)*sin(4*LatRad)

         -(35*eSquare*eSquare*eSquare/3072)*sin(6*LatRad));

     UTMEasting = (double)(k0*V*(A+(1-T+C)*A*A*A/6

                        + (5-18*T+T*T+72*C-58*e2Square)*A*A*A*A*A/120)+ 500000.0);

     UTMNorthing = (double)(k0*(M+V*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24

                        +(61-58*T+T*T+600*C-330*e2Square)*A*A*A*A*A*A/720)));

     //南半球緯度起點為10000000.0m

     UTMNorthing=UTMNorthing+FN;

}

本程式實作的公式請參考《Coordinate Conversions and Transformations including Formulas》文檔第35頁，如果你需要此文檔，請和我聯系！也歡迎大家提出程式的不足。 [轉載請注明出處，謝謝！]