nub = k * (k + 1) / 2
pvv1 = pvv1 - nb(k) * nb(k) / NX(nub) Next k
For k = n To 1 Step -1 '在k循环内回代求解未知数 u = nb(k)
If k < n Then
For i = k + 1 To n
nub = (i - 1) * i / 2 + k u = u - NX(nub) * nb(i) Next i nb(k) = u End If
nub = k * (k + 1) / 2 nb(k) = nb(k) / NX(nub) Next k pvv2 = ll
For k = 1 To n '按[pvv]=[pll]+W*X求[[pvv] pvv2 = pvv2 + UX(k) * nb(k) Next k
For i = 1 To zds '求坐标平差值
cha = charact(i, k) '自定义函数，查点i是否已知点，如不是,用k返回i前面有几个已知点 If cha = \序号i的点不是已知点 h = i - k
d = 2 * (h - 1) + 1 '计算i点x坐标未知数在未知数点集中的序号 x(i) = x(i) + nb(d) / 100: y(i) = y(i) + nb(d + 1) / 100 Else End If Next i
ma = Sqr(pvv1 / (nl(cds) + ns(cds) - n - cds)) '求单位权中误差
Call inversion(NX(), nb()) '将约化后的法方程系数阵送入通用过程inversion()求逆 ReDim nb(n) '工作数组清零
For i = 1 To zds '按点号循环,求点位中误差
cha = charact(i, k) '自定义函数，查点i是否已知点，如不是,用k返回i前面有几个已知点 If cha = \不是已知点
n1 = 2 * (i - k - 1) + 1 '计算i点x坐标未知数在未知数点集中的序号 n2 = n1 * (n1 + 1) / 2 'i点x坐标未知数协因数在协因数阵中序号
13
n3 = (n1 + 1) * (n1 + 2) / 2 'i点y坐标未知数协因数在协因数阵中序号 nb(i) = ma * Sqr(NX(n2) + NX(n3)) End If Next i
CommonDialog1.Flags = 2 '设置对话框flags属性为 2 ,当已有同名文件时询问是否覆盖。 CommonDialog1.ShowSave '打开公共对话框中的\对话框
fname = CommonDialog1.FileName '选择不覆盖同名文件时，返回\另存\对话框。 If fname <> \
Set ts = fso.CreateTextFile(fname, True) 'True参数指出同名文件可被覆盖 c = Space(8) & \平差成果表\
c = c & Space(8) & \点 名\坐标\坐标\点位中误差\ End If ts.Write c
ts.Write blanklines For i = 1 To zds d1 = dm(i)
d1 = Format(d1, \
x1 = x(i): x1 = Format(x1, \ y1 = y(i): y1 = Format(y1, \ mx = nb(i): mx = Format(mx, \
d = Space(5) & d1 & Space(5) & x1 & Space(5) & y1 & Space(5) & mx & Chr(13) & Chr(10) ts.Write d c = c & d Next i
Text1.Text = c '将字符变量c显示在文本框text1中 Text1.Visible = True
MsgBox (\法方程解算完毕\End Select End Sub
Private Function radian(angle As Double) As Double '角度转换为弧度函数 Dim A As Integer, B As String, c As String, d As String, s As String s = str(angle) '将数字转化为字符串 A = InStr(s, \
If A < 1 Then '角度值是整数 B = s Else
14
B = Mid(s, 1, A)
c = Mid(s, A + 2, 2) d = Mid(s, A + 4, 2) End If
radian = pi * (Val(B) + Val(c) / 60 + Val(d) / 3600) / 180 End Function
Private Function azimuth(x1 As Double, y1 As Double, x2 As Double, y2 As Double) '反算坐标方位角 dx = x2 - x1: dy = y2 - y1 azimuth = Atn(dy / dx) If dx < 0 Then
azimuth = azimuth + pi Else
If dy < 0 Then azimuth = azimuth + 2 * pi End If
End Function
Private Function seqn(str As String) As Integer '由点名查计算序号函数 For i = 1 To zds
If str = dm(i) Then seqn = i '将查到的序号赋给函数名,返回调用处 Next i
End Function
Private Function sid(n As Integer, m As Integer) As Double '提取边长函数，n是测站号，m是方向号 ss = 0
If ns(n) - ns(n - 1) - 1 >= 0 Then '该测站有观测边，确定m方向的边长 For k = ns(n - 1) + 1 To ns(n)
If sb(k) = lb(m) Then ss = s(k) '依次将边的照准点名与m方向照准点名对比，相等则赋予ss。 Next k End If
If ns(n) - ns(n - 1) - 1 < 0 Or ss < 1 Then '如测站n未找到m方向的边，则转到lb(m)站去找 h = seqn(lb(m)) '由点名查计算编号 For k = ns(h - 1) + 1 To ns(h)
If sb(k) = dm(n) Then ss = s(k) '依次将k边的照准点名与n测站点名对比，相等则赋予ss。 Next k End If
sid = ss '将查到的边长赋给函数名,返回调用处 End Function
Private Function charact(n1 As Integer, n2 As Integer) As String '根据序号判断是否已知点,小于该序号的点中有几个已知点15
n2 = 0: charact = \
For k = 1 To yds '查那些点是已知点
If n1 = m(k) Then 'm(k)存的是已知点序号，此判断检查序号为i的点是否为已知点 charact = \是已知点 Else
If n1 > m(k) Then n2 = k 'n1不是已知点,用n2记下该序号前有几个已知点,m()数组已经排过序 End If Next k
End Function
Private Sub equation(B() As Double, p As Double, l As Double) '组法方程通用过程 n = 2 * (zds - yds)
For i = 1 To n '按未知数循环
If Abs(B(i)) > 0 Then '第i个未知数的系数不为0 UX(i) = UX(i) + p * B(i) * l '组法方程常数项 For j = 1 To i
If Abs(B(j)) > 0 Then
h = (i - 1) * i / 2 + j '计算xi、xj的互乘系数在一维按列上三角存储法方程系数阵中的序号 NX(h) = NX(h) + B(j) * B(i) * p '组法方程系数阵 End If Next j End If Next i End Sub
Private Sub order(m() As Integer, n As Integer) '排序通用过程 For i = 1 To n - 1 For j = i + 1 To n If m(j) < m(i) Then
A = m(j): m(j) = m(i): m(i) = A End If Next j Next i End Sub
Public Function angle(l As Double) As Double '弧度转角度通用过程 l = l * 180 / pi b1 = Int(l)
16