高斯投影坐标正反算VB程序文件.doc
高斯投影坐标正反算学院：班级：学号：姓名：课程名称: 指导教师：实验目的：1.了解高斯投影坐标正反算的根本思想；2.学会编写高斯正反算程序，加深了解。实验原理：高斯投影正算公式中应满足的三个条件：1. 中央子午线投影后为直线；2. 中央子午线投影后长度不变；3. 投影具有正形性质，即正形投影条件。高斯投影反算公式中应满足的三个条件：1. x坐标轴投影成中央子午线，是投影的对称轴；2. x轴上的长度投影保持不变；3. 正形投影条件，即高斯面上的角度投影到椭球面上后角度没有变形，仍然相等。操作工具：计算机中的VB6.0代码：Dim a As Double, b As Double, x As Double, y As Double, y_#Dim l_ As Double, b_ As Double, a0#, a2#, a4#, a6#, a8#, m2#, m4#, m6#, m8#, m0#, l0#, e#, e1#Dim deg1 As Double, min1 As Double, sec1 As Double, deg2 As Double, min2 As Double, sec2 As DoublePrivate Sub Command1_Click()Dim x_ As Double, t#, eta#, N#, W#, k1#, k2#, ik1%, ik2%, dh%deg1 = Val(Text1.Text)min1 = Val(Text2.Text)sec1 = Val(Text3.Text)deg2 = Val(Text4.Text)min2 = Val(Text5.Text)sec2 = Val(Text6.Text)l_ = (deg1 * 3600 + min1 * 60 + sec1) / 206265b_ = (deg2 * 3600 + min2 * 60 + sec2) / 206265dh = Val(Text9.Text)k1 = (l_ * 180 / 3.14159 + 3) / 6)k2 = (l_ * 180 / 3.14159 / 3)ik1 = Round(k1, 0)ik2 = Round(k2, 0)If dh = 6 Thenl0 = 6 * ik1 - 3ElseIf dh = 3 Thenl0 = 3 * ik2ElseMsgBox "error", 48, "error": Exit SubEnd IfEnd Ifl = l_ - l0 * 3.14159 / 180e = Sqr(a * a - b * b) / am0 = a * (1 - e * e)m2 = e * e * m0 * 3 / 2m4 = e * e * m2 * 5 / 4m6 = m4 * e * e * 7 / 6m8 = e * e * m6 * 9 / 8a0 = m0 + m2 / 2 + m4 * 3 / 8 + m6 * 5 / 16 + m8 * 35 / 128a2 = m2 / 2 + m4 / 2 + m6 * 15 / 32 + m8 * 7 / 16a4 = m4 / 8 + m6 * 3 / 16 + m8 * 7 / 32a6 = m6 / 32 + m8 / 16a8 = m8 / 128x_ = a0 * b_ - a2 * Sin(2 * b_) / 2 + a4 * Sin(4 * b_) / 4 - Sin(6 * b_) * a6 / 6 + Sin(8 * b_) * a8 / 8t = Tan(b_)e1 = Sqr(a * a - b * b) / (b * b)eta = Sqr(e1 * e1 * Cos(b) * Cos(b)W = Sqr(1 - e * e * Sin(b_) * Sin(b_)N = a / Wx = x_ + N * Sin(b_) * Cos(b_) * l * l / 2 + N * Sin(b_) * Cos(b_) 3 * (5 - t * t + 9 * eta * eta + 4 * eta 4) * l 4 / 24 + N * Sin(b_) * Cos(b_) 5 * (61 - 58 * t * t + t 4) * l 6 / 720y = N * Cos(b_) * l + N * Cos(b_) 3 * (1 - t * t + eta * eta) * l * l * l / 6 + N * Cos(b_) 5 * (5 - 18 * t * t + t 4 + 14 * eta * eta - 58 * eta * eta * t * t) * l 5 / 120Text7 = xIf dh = 6 Theny_ = y + 500000 + 1000000 * ik1ElseIf dh = 3 Theny_ = y + 500000 + 1000000 * ik2ElseMsgBox "error", 48, "error": Exit SubEnd IfEnd IfText8 = y_End SubPrivate Sub Command2_Click()Dim bf#, j%, Wf#, Vf#, Nf#, Mf#, c#, tf#, etaf#, dh%, ik%x = Val(Text7.Text)y_ = Val(Text8.Text)e = Sqr(a * a - b * b) / (a * a) m0 = a * (1 - e * e) m2 = e * e * m0 * 3 / 2 m4 = e * e * m2 * 5 / 4 m6 = m4 * e * e * 7 / 6 m8 = e * e * m6 * 9 / 8 a0 = m0 + m2 / 2 + m4 * 3 / 8 + m6 * 5 / 16 + m8 * 35 / 128 a2 = m2 / 2 + m4 / 2 + m6 * 15 / 32 + m8 * 7 / 16 a4 = m4 / 8 + m6 * 3 / 16 + m8 * 7 / 32 a6 = m6 / 32 + m8 / 16 a8 = m8 / 128 bf = x / a0 For j = 1 To 10 bf = (x + a2 / 2 * Sin(2 * bf) - a4 * Sin(4 * bf) / 4 + a6 * Sin(6 * bf) / 6 - a8 * Sin(8 * bf) / 8) / a0 Next j e1 = Sqr(a * a - b * b) / b Vf = Sqr(1 + e1 * e1 * Cos(bf) * Cos(bf) Wf = Sqr(1 - e * e * Sin(bf) * Sin(bf) Nf = a / Wf c = a * a / b Mf = c / Vf 3 tf = Tan(bf) e1 = Sqr(a * a - b * b) / (b * b) etaf = Sqr(e1 * e1 * Cos(bf) * Cos(bf) ik = Int(y_ / 1000000) y = y_ - ik * 1000000 - 500000 b_ = bf - tf * y * y / (2 * Mf * Nf) + tf * (5 + 3 * tf * tf + etaf * etaf - 9 * etaf * etaf * tf * tf) * y * y * y * y / (24 * Mf * Nf * Nf * Nf) + tf * (61 + 90 * tf * tf + 45 * tf * tf * tf * tf) * y * y * y * y * y * y / (720 * Mf * Nf * Nf * Nf * Nf * Nf) l = y / (Nf * Cos(bf) - (1 + 2 * tf * tf + etaf * etaf) * y * y * y / (6 * Nf * Nf * Nf * Cos(bf) + (5 + 28 * tf * tf + 24 * tf * tf * tf * tf + 6 * etaf * etaf + 8 * etaf * etaf * tf * tf) * y * y * y * y * y / (120 * Nf * Nf * Nf * Nf * Nf * Cos(bf)dh = Val(Text9.Text)If dh = 6 Thenl0 = 6 * ik - 3ElseIf dh = 3 Thenl0 = 3 * ikElseMsgBox "error", 48, "error": Exit SubEnd IfEnd Ifl_ = l0 * 3.14159 / 180 + lsec1 = l_ * 206265 deg1 = Int(sec1 / 3600) min1 = Int(sec1 - deg1 * 3600) / 60) sec1 = sec1 - deg1 * 3600 - min1 * 60 sec2 = b_ * 206265 deg2 = Int(sec2 / 3600) min2 = Int(sec2 - deg2 * 3600) / 60) sec2 = sec2 - deg2 * 3600 - min2 * 60 Text1 = deg1 Text2 = min1 Text3 = Round(sec1, 5) Text4 = deg2 Text5 = min2 Text6 = Round(sec2, 5)End SubPrivate Sub Command3_Click()EndEnd SubPrivate Sub Command4_Click()Text1.Text = ""Text2.Text = ""Text3.Text = ""Text4.Text = ""Text5.Text = ""Text6.Text = ""Text7.Text = ""Text8.Text = ""Text9.Text = ""Text10.Text = ""End SubPrivate Sub Option1_Click()a = 6378245b = 6356863.01877305End SubPrivate Sub Option2_Click()a = 6378140b = 6356755.28815753End SubPrivate Sub Option3_Click()a = 6378137b = 6356752.3142End SubPrivate Sub Option4_Click()a = 6378137b = 6356752.3141End SubPrivate Sub Option5_Click()a = Val(InputBox("a", "plsase input")b = Val(InputBox("b", "please input")End Sub界面截图如下：使用方法：1. 如下所示，高斯投影坐标正算时，在经度纬度相应的文本框里输入经纬度，再在度号对应的文本框里输入是六度带还是三度带，最后按坐标正算按钮即可，答案会显示在X，Y相应的文本框里；2. 如果要进展坐标反算，那么输入X，Y与度号，最后按坐标反算按钮即可得到所需的经纬度；3. 注意选择需要的椭球。- 12 - / 12
