gaussbl2xy.cpp

本程序是常用地图投影系列小程序之一

//   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() 
{ 
}