坐标转换代码：我在网上找到下面代码，其中A1, A2, A3, A4代表什么？

Question

我在网上找到下面代码，其中A1,   A2,   A3,   A4代表什么？
坐标转换代码
double   a,   f,   e2,   e12;   //   基本椭球参数   
  double   A1,   A2,   A3,   A4;   //   用于计算X的椭球参数   
Const PI = 3.14159265358979
Const Y0 = 500000#
Dim A1#, A2#, A3#, A4#, a#, f#
'高斯正算求X
Public Function X(B#, L#, L0#, zbx%) As Double  // B纬度L经度 L0中央子午线
zbcs zbx
e2# = 1 - ((f - 1) / f) * ((f - 1) / f)
e12# = (f / (f - 1)) * (f / (f - 1)) - 1
B = D_R(B): L = D_R(L): L0 = D_R(L0)
Dim n, t, t2, m, m2, ng2
Dim sinB, cosB
X = A1 * B * 180# / 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 / Sqr(1 - e2 * sinB * sinB)
m = cosB * (L - L0)
m2 = m * m
ng2 = cosB * cosB * e2 / (1 - e2)
X = X + n * t * ((0.5 + ((5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24# + (61 - 58 * t2 + t2 * t2) * m2 / 720#) * m2) * m2)
End Function
'高斯正算求Y
Public Function Y(B#, L#, L0#, zbx%) As Double
zbcs zbx
e2# = 1 - ((f - 1) / f) * ((f - 1) / f)
e12# = (f / (f - 1)) * (f / (f - 1)) - 1
B = D_R(B): L = D_R(L): L0 = D_R(L0)
Dim n, t, t2, m, m2, ng2
Dim sinB, cosB
sinB = Sin(B)
cosB = Cos(B)
t = Tan(B)
t2 = t * t
n = a / Sqr(1 - e2 * sinB * sinB)
m = cosB * (L - L0)
m2 = m * m
ng2 = cosB * cosB * e2 / (1 - e2)
Y = n * m * (1 + m2 * ((1 - t2 + ng2) / 6# + m2 * (5 - 18 * t2 + t2 * t2 + 14 * ng2 - 58 * ng2 * t2) / 120#))
Y = Y + Y0
End Function
'高斯反算求B（纬度）
Public Function B(X#, Y#, L0#, zbx%) As Double
zbcs zbx
e2# = 1 - ((f - 1) / f) * ((f - 1) / f)
e12# = (f / (f - 1)) * (f / (f - 1)) - 1
L0 = D_R(L0)
Dim sinB, cosB, t, t2, n, ng2, V, yN
Dim preB0, B0
Dim eta
Y = Y - Y0
B0 = X / A1
Do While True
    preB0 = B0
    B0 = B0 * PI / 180#
    B0 = (X - (A2 * Sin(2 * B0) + A3 * Sin(4 * B0) + A4 * Sin(6 * B0))) / A1
    eta = Abs(B0 - preB0)
    If eta < 0.0000001 Then Exit Do
Loop
B0 = B0 * PI / 180#sinB = Sin(B0)
cosB = Cos(B0)
t = Tan(B0)
t2 = t * t
n = a / Sqr(1 - e2 * sinB * sinB)
ng2 = cosB * cosB * e2 / (1 - e2)
V = Sqr(1 + ng2)
yN = Y / n
B = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN / 12# + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360#) * V * V * t / 2
B = R_D(B)
End Function'高斯反算求L（经度）
Public Function L(X#, Y#, L0#, zbx%) As Double
zbcs zbx
e2# = 1 - ((f - 1) / f) * ((f - 1) / f)
e12# = (f / (f - 1)) * (f / (f - 1)) - 1
L0 = D_R(L0)
Dim sinB, cosB, t, t2, n, ng2, V, yN
Dim preB0, B0
Dim eta
Y = Y - Y0
B0 = X / A1
Do While True
    preB0 = B0
    B0 = B0 * PI / 180#
    B0 = (X - (A2 * Sin(2 * B0) + A3 * Sin(4 * B0) + A4 * Sin(6 * B0))) / A1
    eta = Abs(B0 - preB0)
    If eta < 0.0000001 Then Exit Do
Loop
B0 = B0 * PI / 180#sinB = Sin(B0)
cosB = Cos(B0)
t = Tan(B0)
t2 = t * t
n = a / Sqr(1 - e2 * sinB * sinB)
ng2 = cosB * cosB * e2 / (1 - e2)
V = Sqr(1 + ng2)
yN = Y / n
L = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6# + (5 + 28 * t2 + 24 * t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120#) / cosB
L = R_D(L)
End Function
'把弧度转化为DD.FFMM格式
Public Function R_D(R) As Double
R = R * 180 / PI
DD = Int(R)
FF = Int((R - DD) * 60#)
MM = ((R - DD) * 60# - FF) * 60#If Abs(MM - 60) < 0.0001 Then
MM = 0: FF = FF + 1
End If
If FF = 60 Then
FF = 0: DD = DD + 1
End IfR_D = DD + FF / 100 + MM / 10000
End Function
'把DD.FFMM转化为弧度
Public Function D_R(DFM) As Double
'DFM = Round(DFM, 4)
Dim DD%, FF%, MM#
DD = Int(DFM)
FF = Int((DFM - Int(DFM)) * 100 + 0.001)
MM = ((DFM * 100) - Int(Round(DFM * 100, 3))) * 100
D_R = (DD + FF / 60 + MM / 3600) * PI / 180End FunctionSub zbcs(DaiHao As Integer)
    Select Case DaiHao
      Case 1
      A1 = 111134.8611
      A2 = -16036.4803
      A3 = 16.8281
      A4 = -0.022
      a = 6378245#
      f = 298.3
      Case 2
      A1 = 111133.0047
      A2 = -16038.5282
      A3 = 16.8326
      A4 = -0.022
      a = 6378140     1975 year，IUUG，长轴 6378140      f = 298.257(扁率)
      Case Else
      A1 = 111134.8611
      A2 = -16036.4803
      A3 = 16.8281
      A4 = -0.022
      a = 6378245#
      f = 298.3
    End Select
End Sub
几种常见的椭球体参数值
        克拉索夫斯基椭球体        1975年国际椭球体        WGS-84椭球体
a
        6 378 245.000 000 000 0（m）        6 378 140.000 000 000 0（m）        6 378 137.000 000 000 0（m）
b        6 356 863.018 773 047 3（m）        6 356 755.288 157 528 7（m）        6 356 752.314 2（m）
c        6 399 698.901 782 711 0（m）        6 399 596.651 988 010 5（m）        6 399 593.625 8（m）
α        1／298.3        1／298.257        1/298.257 223 563
e2        0.006 693 421 622 966        0.006 694 384 999 588        0.006 694 379 901 3
e’2        0.006 738 525 414 683        0.006 739 501 819 473        0.006 739 496 742 27 
A1,   A2,   A3,   A4

Answer 1

坐标转换
我国规定1：1万、1：2.5万、1：5万、1：10万、1：25万、1：50万比例尺地形图，均采用高斯克吕格投影。1：2.5至1：50万比例尺地形图采用经差6度分带，1：1万比例尺地形图采用经差3度分带。
6度带是从0度子午线起，自西向东每隔经差6为一投影带，全球分为60带，各带的带号用自然序数1，2，3，…60表示。即以东经0-6为第1带，其中央经线为3E，东经6-12为第2带，其中央经线为9E，其余类推（图4-13）。
3度带，是从东经1度30分的经线开始，每隔3度为一带，全球划分为120个投影带。图4-13表示出6度带与3度带的中央经线与带号的关系。
在高斯克吕格投影上，规定以中央经线为X轴，赤道为Y轴，两轴的交点为坐标原点。
X坐标值在赤道以北为正，以南为负；Y坐标值在中央经线以东为正，以西为负。我国在北半球，X坐标皆为正值。Y坐标在中央经线以西为负值，运用起来很不方便。为了避免Y坐标出现负值，将各带的坐标纵轴西移500公里，即将所有Y值都加500公里。
如何计算当地的中央子午线？ 
当地中央子午线决定于当地的直角坐标系统，首先确定您的直角坐标系统是3度带还是6度带投影，然后再根据如下公式推算： 
6度带中央子午线计算公式：　　　　当地经度/6=N；中央子午线L=6 X N 
当没有除尽，N有余数时，　　　　　中央子午线L=6 X N - 3 
3度带中央子午线计算公式： 
当地经度/3=N；中央子午线L=3 X N
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; 
}