高斯正反算

VB下高斯坐标变换的实现
曾圣陈伟
高斯-克吕格(Gauss-Kruger)投影，是一种“等角横切圆柱投影”。

德国数学家、物理学家、天文学家高斯（Carl Friedrich Ｇauss，1777一 1855）于十九世纪二十年代拟定，后经德国大地测量学家克吕格（Johannes Kruger，1857～1928）于 1912年对投影公式加以补充，故名。

设想用一个圆柱横切于球面上投影带的中央经线，按照投影带中央经线投影为直线且长度不变和赤道投影为直线的条件，将中央经线两侧一定经差范围内的球面正形投影于圆柱面.
高斯-克吕格(Gauss-Kruger)投影是横轴墨卡托投影的变种，高斯-克吕格投影是“等角横切圆柱投影”，投影后中央经线保持长度不变，即比例系数为1。

高斯-克吕格投影的东伪偏移是500公里，投影北伪偏移为零。

高斯-克吕格投影正解公式：已知（B,L）值求解(X,Y)，（原点纬度 0，中央经度L0）
上面公式中东纬偏移FE = 500000米（本程序中中只设计加500000米的常数，如还要带号的话FE = 500000+ 带号 * 1000000米）；高斯-克吕格投影比例因子k0 = 1
高斯-克吕格投影反解公式：已知(X,Y)求解（B,L），（原点纬度 0，中央经度
L0）
a -- 椭球体长半轴
b -- 椭球体短半轴
f -- 扁率
e -- 第一偏心率
e’ -- 第二偏心率
N -- 卯酉圈曲率半径
R -- 子午圈曲率半径
B -- 纬度，L -- 经度，单位弧度(RAD)
X -- 纵直角坐标, Y -- 横直角坐标，单位米(M) 关于椭球体参数，我国常用的3个椭球体参数如下：
程序实现:
'高斯正算求X
Public Function X(ByVal B#, ByVal L#, ByVal L0#) As Double
Dim n#, T#, T2#, m#, m2#, ng2#
Dim S#, C#
X = A1 * B + A2 * Sin(2 * B) + A3 * Sin(4 * B) + A4 * Sin(6 * B)? '子午线弧长
S = Sin(B)
C = Cos(B)
T = Tan(B)
T2 = T * T
n = a / Sqr(1 - e12 * S * S)'卯酉圈曲率半径
m = C * (L - L0)
m2 = m * m
ng2 = C * C * e12 / (1 - e12)????
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(ByVal B#, ByVal L#, ByVal L0#) As Double
Dim n#, T#, T2#, m#, m2#, ng2#
Dim S#, C#
S = Sin(B)
C = Cos(B)
T = Tan(B)
T2 = T * T
n = a / Sqr(1 - e12 * S * S)
m = C * (L - L0)
m2 = m * m
ng2 = C * C * e12 / (1 - e12)
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(ByVal X#, ByVal Y#) As Double
Dim S#, C#, T#, T2#, n#, ng2#, V#, yN#
Dim preB0#, B0#
Dim eta#
Y = Y - Y0
B0 = X / A1
Do
preB0 = B0
B0 = (X - (A2 * Sin(2 * B0) + A3 * Sin(4 * B0) + A4 * Sin(6 * B0))) / A1 If Abs(B0 - preB0) < 0.000000001 Then Exit Do
Loop
S = Sin(B0)
C = Cos(B0)
T = Tan(B0)
T2 = T * T
n = a / Sqr(1 - e12 * S * S)
ng2 = C * C * e12 / (1 - e12)
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#
End Function
程序运行界面:
检测数据:
（1）、已知在北京坐标系下
中央子午线：117度，纬度B：28度32分14.5秒，经度L：116度54分12.3秒
正算求解北方向X、东方向Y值。

得:
X= 3158054.118米Y= 490547.265米
（2）、已知在北京坐标系下
中央子午线：114度，北方向X值：2834561.381米、东方向Y值：583726.735米
反算求解B、L值。

得:
B= 25度36分56.1秒L= 114度50分00.8秒
（3）、已知在北京坐标系下
换带前中央子午线：117度，换带后中央子午线：116度20分，北方向X值：3158038.731米，东方向Y值：490539.097米
换带计算中央子午线116度20分下新的X、Y值。

得
X= 3158167.489米Y= 555787.193米。