当前位置:文档之家› 高斯投影坐标正反算VB程序

高斯投影坐标正反算VB程序

高斯投影坐标正反算VB程序
高斯投影坐标正反算VB程序

高斯投影坐标正反算

V B程序

文件编码(GHTU-UITID-GGBKT-POIU-WUUI-8968)

高斯投影坐标正反算

学院:

班级:

学号:

姓名:

课程名称:

指导老师:

实验目的:

1.了解高斯投影坐标正反算的基本思想;

2.学会编写高斯正反算程序,加深了解。

实验原理:

高斯投影正算公式中应满足的三个条件:

1. 中央子午线投影后为直线;

2. 中央子午线投影后长度不变;

3. 投影具有正形性质,即正形投影条件。

高斯投影反算公式中应满足的三个条件:

1. x坐标轴投影成中央子午线,是投影的对称轴;

2. x轴上的长度投影保持不变;

3. 正形投影条件,即高斯面上的角度投影到椭球面上后角度没有

变形,仍然相等。

操作工具:

计算机中的

代码:

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 Double

Private Sub Command1_Click()

Dim x_ As Double, t#, eta#, N#, W#, k1#, k2#, ik1%, ik2%, dh% deg1 = Val

min1 = Val

sec1 = Val

deg2 = Val

min2 = Val

sec2 = Val

l_ = (deg1 * 3600 + min1 * 60 + sec1) / 206265

b_ = (deg2 * 3600 + min2 * 60 + sec2) / 206265

dh = Val

k1 = ((l_ * 180 / + 3) / 6)

k2 = (l_ * 180 / / 3)

ik1 = Round(k1, 0)

ik2 = Round(k2, 0)

If dh = 6 Then

l0 = 6 * ik1 - 3

Else

If dh = 3 Then

l0 = 3 * ik2

Else

MsgBox "error", 48, "error": Exit Sub

End If

End If

l = l_ - l0 * / 180

e = Sqr(a * a - b * b) / 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

x_ = a0 * b_ - a2 * Sin(2 * b_) / 2 + a4 * Sin(4 * b_) / 4 - Sin(6 * b_) * a6 / 6 + Sin(8 * b_) * a8 / 8

t = 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 / W

x = 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 / 720

y = 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 / 120

Text7 = x

If dh = 6 Then

y_ = y + 500000 + 1000000 * ik1

Else

If dh = 3 Then

y_ = y + 500000 + 1000000 * ik2

Else

MsgBox "error", 48, "error": Exit Sub

End If

End If

Text8 = y_

End Sub

Private Sub Command2_Click()

Dim bf#, j%, Wf#, Vf#, Nf#, Mf#, c#, tf#, etaf#, dh%, ik%

x = Val

y_ = Val

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

If dh = 6 Then

l0 = 6 * ik - 3

Else

If dh = 3 Then

l0 = 3 * ik

Else

MsgBox "error", 48, "error": Exit Sub End If

End If

l_ = l0 * / 180 + l

sec1 = 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 Sub

Private Sub Command3_Click() End

End Sub

Private Sub Command4_Click() = ""

= ""

= ""

= ""

= ""

= ""

= ""

= ""

= ""

= ""

End Sub

Private Sub Option1_Click() a = 6378245

b =

End Sub

Private Sub Option2_Click()

a = 6378140

b = 6356755.

End Sub

Private Sub Option3_Click()

a = 6378137

b =

End Sub

Private Sub Option4_Click()

a = 6378137

b =

End Sub

Private Sub Option5_Click()

a = Val(InputBox("a", "plsase input"))

b = Val(InputBox("b", "please input")) End Sub

界面截图如下:

使用方法:

1.如下所示,高斯投影坐标正算时,在经度纬度相应的文本框里输入经

纬度,再在度号对应的文本框里输入是六度带还是三度带,最后按坐标正算按钮即可,答案会显示在X,Y相应的文本框里;

2.如果要进行坐标反算,则输入X,Y与度号,最后按坐标反算按钮即

可得到所需的大地经纬度;

3.注意选择需要的椭球。

相关主题
文本预览
相关文档 最新文档