经纬度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   */