测量专业基本计算与编程(10)

w(i) = 0 For j = i To 4

IJ = j * (j - 1) / 2 + i c(IJ) = 0 Next j: Next i

For i = 1 To newnum c(1) = c(1) + 1 c(2) = 0

c(3) = c(3) + 1

c(4) = c(4) + oldx(i) c(5) = c(5) + oldy(i)

c(6) = c(6) + oldx(i) * oldx(i) + oldy(i) * oldy(i) c(7) = c(5) c(8) = -c(4) c(9) = 0 c(10) = c(6)

w(1) = w(1) + newx(i) w(2) = w(2) + newy(i)

w(3) = w(3) + oldx(i) * newx(i) + oldy(i) * newy(i) w(4) = w(4) + oldy(i) * newx(i) - oldx(i) * newy(i) Next i

'解法方程系数矩阵求逆 nx = 4

ReDim xx(1 To nx)

Dim H(1 To 2000) As Double, h1 As Double, h2 As Double For k = nx To 1 Step -1 h1 = c(1): i1 = 1 For i = 2 To nx

i2 = i1: i1 = i1 + i h2 = c(i2 + 1) H(i) = h2 / h1

If i <= k Then H(i) = -H(i) J1 = i2 + 2 For j = J1 To i1

c(j - i) = c(j) + h2 * H(j - i2) Next j Next i

i2 = i2 - 1: c(i1) = 1 / h1 For i = 2 To nx c(i2 + i) = H(i) Next i Next k

'求转换未知数 For i = 1 To nx

xx(i) = 0

For j = 1 To nx If j > i Then

KK0 = j * (j - 1) / 2 + i Else

KK0 = i * (i - 1) / 2 + j End If

xx(i) = xx(i) + c(KK0) * w(j) Next j Next i

m = Sqr(xx(3) * xx(3) + xx(4) * xx(4)) - 1 q = Atn(xx(4) / xx(3)): q = qdms(q) For i = 1 To oldnum

zhx(i) = xx(1) + xx(3) * oldx(i) + xx(4) * oldy(i) zhy(i) = xx(2) + xx(3) * oldy(i) - xx(4) * oldx(i) Next i

Call save_result

Open projectdir + \

Print #1, \

Print #1, \ 4.转换参数 Print #1, \Print #1, \轴平移参数= \Print #1, \轴平移参数= \Print #1, \旋转参数= \度.分秒)\Print #1, \尺度比参数= \

Print #1, \Close #1

Value = MsgBox(\坐标转换计算已完成!\坐标转换(\相似变换法）\Exit Sub Line1:

msg = \错误：\

Value = MsgBox(msg, 32, \错误信息窗口\End Sub

Private Sub computer3() ' 多项式逼近法

Dim KK As Integer, IJ As Integer

Dim a(15) As Double, l(1 To 8) As Double, m(1 To 8) As Double Dim xy1 As Double, xy2 As Double, xy3 As Double Dim ll0 As Double, ll1 As Double ReDim c(200), w(15) Dim w1(15) As Double

Dim px As Double, py As Double On Error GoTo Line1

\ Open projectdir + \oldnum = 0

Do While Not EOF(1) oldnum = oldnum + 1 Input #1, nno$, xy1, xy2

oldname(oldnum) = nno$: oldx(oldnum) = xy1: oldy(oldnum) = xy2 Loop Close #1

Open projectdir + \newnum = 0

Do While Not EOF(1)

newnum = newnum + 1 Input #1, nno$, xy1, xy2

newname(newnum) = nno$: newx(newnum) = xy1: newy(newnum) = xy2 Loop Close #1

'---根据新坐标调整旧坐标的顺序并确定公共点--- For J1 = 1 To newnum For j2 = 1 To oldnum

If J1 <> j2 And newname(J1) = oldname(j2) Then

nno$ = oldname(J1): oldname(J1) = oldname(j2): oldname(j2) = nno$ xy3 = oldx(J1): oldx(J1) = oldx(j2): oldx(j2) = xy3 xy3 = oldy(J1): oldy(J1) = oldy(j2): oldy(j2) = xy3 Exit For End If Next j2 Next J1

'--多项式逼近组成法方程-- nx = 6

For i = 1 To nx

w(i) = 0#: w1(i) = 0# For j = i To nx

IJ = j * (j - 1) / 2 + i c(IJ) = 0# Next j Next i

'--求公共点重心坐标-- px = 0: py = 0

For i = 1 To newnum px = px + oldx(i) py = py + oldy(i) Next i

px = px / newnum: py = py / newnum '--求误差方程系数并组成法方程--

For i = 1 To newnum

a(1) = 1#：a(2) = oldx(i) - px： a(3) = oldy(i) - py

a(4) = a(2) * a(2)：a(5) = a(3) * a(3)：a(6) = a(2) * a(3) ll0 = oldx(i) - newx(i) ll1 = oldy(i) - newy(i) For J1 = 1 To nx

w(J1) = w(J1) + a(J1) * ll0 w1(J1) = w1(J1) + a(J1) * ll1 For j2 = J1 To nx

IJ = j2 * (j2 - 1) / 2 + J1 c(IJ) = c(IJ) + a(J1) * a(j2) Next j2 Next J1 Next i

'-法方程系数矩阵求逆- ReDim xx(1 To nx)

Dim H(1 To 200) As Double, h1 As Double, h2 As Double For k = nx To 1 Step -1 h1 = c(1): i1 = 1 For i = 2 To nx

i2 = i1: i1 = i1 + i h2 = c(i2 + 1) H(i) = h2 / h1

If i <= k Then H(i) = -H(i) J1 = i2 + 2 For j = J1 To i1

c(j - i) = c(j) + h2 * H(j - i2) Next j Next i

i2 = i2 - 1: c(i1) = 1 / h1 For i = 2 To nx c(i2 + i) = H(i) Next i Next k

'-求法方程未知数 For i = 1 To nx xx(i) = 0

For j = 1 To nx If j > i Then

IJ = j * (j - 1) / 2 + i Else

IJ = i * (i - 1) / 2 + j End If

xx(i) = xx(i) - c(IJ) * w(j)

Next j Next i

'---------------------------------- For i = 1 To nx l(i) = xx(i) Next i

'--------------------------------- For i = 1 To nx xx(i) = 0

For j = 1 To nx If j > i Then

IJ = j * (j - 1) / 2 + i Else

IJ = i * (i - 1) / 2 + j End If

xx(i) = xx(i) - c(IJ) * w1(j) Next j Next i

For i = 1 To nx m(i) = xx(i) Next i

' 求转换坐标

For i = 1 To oldnum

a(1) = 1#: a(2) = oldx(i) - px：a(3) = oldy(i) - py

a(4) = a(2) * a(2)：a(5) = a(3) * a(3)：a(6) = a(2) * a(3) zhx(i) = oldx(i) + l(1) + l(2) * a(2) + l(3) * a(3) zhx(i) = zhx(i) + l(4) * a(4) + l(5) * a(5)

zhy(i) = oldy(i) + m(1) + m(2) * a(2) + m(3) * a(3) zhy(i) = zhy(i) + m(4) * a(4) + m(5) * a(5) Next i

Call save_result

Value = MsgBox(\坐标转换计算已完成!\坐标转换(\多项式拟合法）\Exit Sub Line1:

msg = \错误：\

Value = MsgBox(msg, 32, \错误信息窗口\End Sub

Private Sub save_result()

Open projectdir + \

Print #1, \ 新旧坐标转换计算成果 \Print #1, \

Print #1, \ 1.旧坐标系统坐标 \Print #1, \