经纬度BL换算到高斯平面直角坐标XY（高斯投影正算）的源码及算法

一、
经纬度BL换算到高斯平面直角坐标XY（高斯投影正算）的源码及算法

在   gis   的帖子 "用excel完成gps坐标转换的简易方法   "   基础上,
我给出对应的vb程序段,我在   evb   开发的gps定位功能中，用它实现坐标换算(具体的参数请参考gis   的帖子)。
感觉速度比较快，效果比较好。所以帖上来，希望与名位交流：
=====================================
'经纬度bl换算到高斯平面直角坐标xy（高斯投影正算）

 

private     function    bl2xy( byref    a2    as     double ,    byref    f2    as     double ,    byref    e2    as     double ,   _ 
byref    s2    as     double ,    byref    t2    as     double )    as     boolean  
' a2   输入中央子午线，以度.分形式输入，如115度30分则输入115.30;   起算数据l0 
' f2   以度小数形式输入经度值,   l 
' e2   以度小数形式输入纬度值，b 
' s2   计算结果，横坐标y 
' t2   计算结果，纵坐标x 
' 投影带号计算　n=[l/6]+1   如：测得经度103.xxxx，故n=[103.x/6]+1=17+1=18 
' 中央经线经度　l0   =   n*6-3   =   [l/6]*6+3 

dim    b2    as     double  
' dim   g2   as   double 
dim    h2    as     double  
dim    i2    as     double  
dim    j2    as     double  
dim    k2    as     double  
dim    l2    as     double  
dim    m2    as     double  
dim    n2    as     double  
dim    o2    as     double  
dim    p2    as     double  
dim    q2    as     double  
dim    r2    as     double  

b2    =     int (a2)    +    ( int (a2    *     100 )    -     int (a2)    *     100 )    /     60     +    (a2    *     10000     -     int (a2    *     100 )    *     100 )    /     3600  
' 把l0化成度(a2) 
' g2   =   f2   -   b2   '   l   -l0 
' h2   =   g2   /   57.2957795130823   '化作弧度 
h2    =    (f2    -    b2)    /     57.2957795130823     ' 将经差的单位化为弧度 
i2    =     tan (e2    /     57.2957795130823 )    ' tan   (b) 
j2    =     cos (e2    /     57.2957795130823 )    '    cos   (b) 
k2    =     0.006738525415     *    j2    *    j2 
l2    =    i2    *    i2 
m2    =     1     +    k2 
n2    =     6399698.9018     /     sqr (m2) 
o2    =    h2    *    h2    *    j2    *    j2 
p2    =    i2    *    j2 
q2    =    p2    *    p2 
r2    =    ( 32005.78006     +    q2    *    ( 133.92133     +    q2    *     0.7031 )) 
s2    =    ((((l2    -     18 )    *    l2    -    ( 58     *    l2    -     14 )    *    k2    +     5 )    *    o2    /     20     +    m2    -    l2)    *    o2    /     6     +     1 )    *    n2    *    (h2    *    j2) 
s2    =    s2    +     18500000     ' 在计算的基础上加上了“带号”（18）和“东移”（500km） 
' 计算结果,横坐标y 
t2    =     6367558.49686     *    e2    /     57.29577951308     -    p2    *    j2    *    r2    +    ((((l2    -     58 )    *    l2    +     61 )    *    _ 
o2    /     30     +    ( 4     *    k2    +     5 )    *    m2    -    l2)    *    o2    /     12     +     1 )    *    n2    *    i2    *    o2    /     2  
' 计算结果，纵坐标x 
' msgbox   "pts2= "   &   s2   &   "   pt   t2= "   &   t2 

bl2xy    =     true  
end     function

 二、

 

//    GaussBL2xy.cpp   :   Defines   the   entry   point   for   the   console   application. 
//  

#include    " stdafx.h  "  
#include    " math.h  "  
#include    " CoorTrans.h  "  
#include    < iostream >  

using     namespace    std; 

void    main( int    argc,    char *    argv[]) 
{ 
   double    MyL0    =     108 ;    // 中央子午线 
   double    MyB    =     33.44556666 ;    // 33   du   44   fen   55.6666   miao 
   double    MyL    =     109.22334444 ;    // 3度带，109   d   22   m   33.4444   s 
  
  PrjPoint_Krasovsky   MyPrj; 
  MyPrj.SetL0(MyL0); 
  MyPrj.SetBL(MyB,   MyL); 
   double    OutMyX;                        // 结果应该大致是:3736714.783,   627497.303 
   double    OutMyY; 
  OutMyX    =    MyPrj.x;            // 正算结果:   北坐标x 
  OutMyY    =    MyPrj.y;            // 结果：东坐标y 
  
   // 反算 /// / 
   double    InputMyX    =     3736714.783 ;      // 如果是独立计算，应该给出中央子午线L0 
   double    InputMyY    =     627497.303 ; 
  MyPrj.Setxy(InputMyX,   InputMyY); 
  MyPrj.GetBL( & MyPrj.B,    & MyPrj.L);    // 把计算出的BL的弧度值换算为dms形式 
   double    OutputMyB; 
   double    OutputMyL; 
  OutputMyB    =    MyPrj.B;      // 反算结果：B 
  OutputMyL    =    MyPrj.L;      // 反算结果：L 

   // 分析表明，此程序的结果和Coord4.2的转换结果是一样的，只差到毫米级 
         // 原程序有几个问题，1.Pi的值不对。2.SetBL中多了两行错误代码 
} 


double    Dms2Rad( double    Dms) 
{ 
   double    Degree,   Miniute; 
   double    Second; 
   int    Sign; 
   double    Rad; 
   if (Dms    >   =     0 ) 
    Sign    =     1 ; 
   else  
    Sign    =     - 1 ; 
  Dms    =    fabs(Dms); 
  Degree    =    floor(Dms); 
  Miniute    =    floor(fmod(Dms    *     100.0 ,    100.0 )); 
  Second    =    fmod(Dms    *     10000.0 ,    100.0 ); 
  Rad    =    Sign    *    (Degree    +    Miniute    /     60.0     +    Second    /     3600.0 )    *    PI    /     180.0 ; 
   return    Rad; 
} 
double    Rad2Dms( double    Rad) 
{ 
   double    Degree,   Miniute; 
   double    Second; 
   int    Sign; 
   double    Dms; 
   if (Rad    >   =     0 ) 
    Sign    =     1 ; 
   else  
    Sign    =     - 1 ; 
  Rad    =    fabs(Rad    *     180.0     /    PI); 
  Degree    =    floor(Rad); 
  Miniute    =    floor(fmod(Rad    *     60.0 ,    60.0 )); 
  Second    =    fmod(Rad    *     3600.0 ,    60.0 ); 
  Dms    =    Sign    *    (Degree    +    Miniute    /     100.0     +    Second    /     10000.0 ); 
   return    Dms; 
} 
///  
//    Definition   of   PrjPoint 
///  
bool    PrjPoint::BL2xy() 
{ 
   double    X,   N,   t,   t2,   m,   m2,   ng2; 
   double    sinB,   cosB; 
  X    =    A1    *    B    *     180.0     /    PI    +    A2    *    sin( 2     *    B)    +    A3    *    sin( 4     *    B)    +    A4    *    sin( 6     *    B); 
  sinB    =    sin(B); 
  cosB    =    cos(B); 
  t    =    tan(B); 
  t2    =    t    *    t; 
  N    =    a    /    sqrt( 1     -    e2    *    sinB    *    sinB); 
  m    =    cosB    *    (L    -    L0); 
  m2    =    m    *    m; 
  ng2    =    cosB    *    cosB    *    e2    /    ( 1     -    e2); 
   // x,y的计算公式见孔祥元等主编武汉大学出版社2002年出版的《控制测量学》的第72页 
   // 书的的括号有问题，   (   和   [   应该交换 
  x    =    X    +    N    *    t    *    (( 0.5     +    (( 5     -    t2    +     9     *    ng2    +     4     *    ng2    *    ng2)    /     24.0     +    ( 61     -      
     58     *    t2    +    t2    *    t2)    *    m2    /     720.0 )    *    m2)    *    m2); 
  y    =    N    *    m    *    (    1     +    m2    *    (   ( 1     -    t2    +    ng2)    /     6.0     +    m2    *    (    5     -     18     *    t2    +    t2    *    t2 
     +     14     *    ng2    -     58     *    ng2    *    t2   )    /     120.0 )); 
  y    +=     500000 ; 
     return     true ; 
} 
bool    PrjPoint::xy2BL() 
{ 
   double    sinB,   cosB,   t,   t2,   N   ,ng2,   V,   yN; 
   double    preB0,   B0; 
   double    eta; 
  y    -=     500000 ; 
  B0    =    x    /    A1; 
   do  
  { 
    preB0    =    B0; 
    B0    =    B0    *    PI    /     180.0 ; 
    B0    =    (x    -    (A2    *    sin( 2     *    B0)    +    A3    *    sin( 4     *    B0)    +    A4    *    sin( 6     *    B0)))    /    A1; 
    eta    =    fabs(B0    -    preB0); 
  } while (eta    >     0.000000001 ); 
  B0    =    B0    *    PI    /     180.0 ; 
  B    =    Rad2Dms(B0); 
  sinB    =    sin(B0); 
  cosB    =    cos(B0); 
  t    =    tan(B0); 
  t2    =    t    *    t; 
  N    =    a    /    sqrt( 1     -    e2    *    sinB    *    sinB); 
  ng2    =    cosB    *    cosB    *    e2    /    ( 1     -    e2); 
  V    =    sqrt( 1     +    ng2); 
  yN    =    y    /    N; 
B    =    B0    -    (yN    *    yN    -    ( 5     +     3     *    t2    +    ng2    -     9     *    ng2    *    t2)    *    yN    *    yN    *    yN    *    yN    /  
     12.0     +    ( 61     +     90     *    t2    +     45     *    t2    *    t2)    *    yN    *    yN    *    yN    *    yN    *    yN    *    yN    /     360.0 ) 
     *    V    *    V    *    t    /     2 ; 
  L    =    L0    +    (yN    -    ( 1     +     2     *    t2    +    ng2)    *    yN    *    yN    *    yN    /     6.0     +    ( 5     +     28     *    t2    +     24  
     *    t2    *    t2    +     6     *    ng2    +     8     *    ng2    *    t2)    *    yN    *    yN    *    yN    *    yN    *    yN    /     120.0 )    /    cosB; 
     return     true ; 
} 
bool    PrjPoint::SetL0( double    dL0) 
{ 
  L0    =    Dms2Rad(dL0); 
     return     true ; 
} 
bool    PrjPoint::SetBL( double    dB,    double    dL) 
{ 
  B    =    Dms2Rad(dB); 
  L    =    Dms2Rad(dL); 
         // B   =   dB;          // 我靠，I   wana   say   fuck 
   // L   =   dL;    // del   it   ! 
  BL2xy(); 
     return     true ; 
} 
bool    PrjPoint::GetBL( double     * dB,    double     * dL) 
{ 
   * dB    =    Rad2Dms(B); 
   * dL    =    Rad2Dms(L); 
   return     true ; 
} 
bool    PrjPoint::Setxy( double    dx,    double    dy) 
{ 
  x    =    dx; 
  y    =    dy; 
  xy2BL(); 
   return     true ; 
} 
bool    PrjPoint::Getxy( double     * dx,    double     * dy) 
{ 
   * dx    =    x; 
   * dy    =    y; 
     return     true ; 
} 
///  
//    Definition   of   PrjPoint_IUGG1975 
///  
PrjPoint_IUGG1975::PrjPoint_IUGG1975()    // 在类外定义构造成员函数，要加上类名和域限定符   :: 
{ 
  a    =     6378140 ; 
  f    =     298.257 ; 
  e2    =     1     -    ((f    -     1 )    /    f)    *    ((f    -     1 )    /    f); 
  e12    =    (f    /    (f    -     1 ))    *    (f    /    (f    -     1 ))    -     1 ; 
  A1    =     111133.0047 ;      // 这几个A是什么意思？ 
  A2    =     - 16038.5282 ; 
  A3    =     16.8326 ; 
  A4    =     - 0.0220 ; 
} 
PrjPoint_IUGG1975:: ~ PrjPoint_IUGG1975()        // 析构函数，释放构造函数占用的内存 
{ 
} 
///  
//    Definition   of   PrjPoint_Krasovsky 
///  
PrjPoint_Krasovsky::PrjPoint_Krasovsky() 
{ 
  a    =     6378245 ; 
  f    =     298.3 ; 
  e2    =     1     -    ((f    -     1 )    /    f)    *    ((f    -     1 )    /    f); 
  e12    =    (f    /    (f    -     1 ))    *    (f    /    (f    -     1 ))    -     1 ; 
  A1    =     111134.8611 ; 
  A2    =     - 16036.4803 ; 
  A3    =     16.8281 ; 
  A4    =     - 0.0220 ; 
} 
PrjPoint_Krasovsky:: ~ PrjPoint_Krasovsky() 
{ 
} 

  

// CoorTrans.h   文件 
#ifndef   _COORTRANS_H_INCLUDED 
#define    _COORTRANS_H_INCLUDED 

#include    " stdafx.h  "  
#include    " math.h  "  

const     double    PI    =     3.141592653589793 ; 

double    Dms2Rad( double    Dms); 
double    Rad2Dms( double    Rad); 

class    PrjPoint 
{ 
public : 
   double    L0;    //    中央子午线经度 
   double    B,   L;    //    大地坐标 
   double    x,   y;    //    高斯投影平面坐标 
public : 
   bool    BL2xy(); 
   bool    xy2BL(); 

protected : 
   double    a,   f,   e2,   e12;    //    基本椭球参数 
   double    A1,   A2,   A3,   A4;    //    用于计算X的椭球参数 

public : 
   bool    SetL0( double    dL0); 
   bool    SetBL( double    dB,    double    dL); 
   bool    GetBL( double     * dB,    double     * dL); 
   bool    Setxy( double    dx,    double    dy); 
   bool    Getxy( double     * dx,    double     * dy); 
}; 

class    PrjPoint_Krasovsky   :    virtual     public    PrjPoint    // 类的继承，声明基类是   PrjPoint，虚基类 
{ 
public : 
  PrjPoint_Krasovsky();    // 定义构造函数，用来初始化。(函数名和类名相同) 
   ~ PrjPoint_Krasovsky(); 
}; 

class    PrjPoint_IUGG1975   :    virtual     public    PrjPoint 
{ 
public : 
  PrjPoint_IUGG1975(); 
   ~ PrjPoint_IUGG1975(); 
}; 

#endif    /*   ndef   _COORTRANS_H_INCLUDED   */