平面控制网间接平差.frm

测量间接平差vb源代码,测量工作者的福音

'平差计算
Private Sub mnuAdjust_Click()
    Screen.MousePointer = 13
    
    '列立误差方程===================================================================================
    Dim strR As String, S3$             '用于显示内容的字符串
    Dim i%, j%, k%                      '循环变量
    Dim ii%, jj%, kk%                   '辅助记数变量
    Dim dX#, dY#, ss#, cosa#, sinA#     '列边长误差方程所需中间变量
    Dim Z0#, Z1#, a0#                   '列方向误差方程所需中间变量
    Dim ai#, bi#, P3(300) As Single     '一些辅助变量和数组
    Dim k1%, k2%, pPp(300) As Integer   '一些辅助变量和数组
    
    If Ne > 0 Then '按边长计算方程系数
        strR = "         边长误差方程：" & vbCrLf
        For i = 1 To Ne
            dX = X0(eE(i)) - X0(bE(i)):        dY = Y0(eE(i)) - Y0(bE(i))
            ss = Sqr(dX ^ 2 + dY ^ 2):         cosa = dX / ss: sinA = dY / ss
            '系数阵及常数项
            Pa(i, 1) = 2 * (bE(i) - nn) - 1:   Pa(i, 2) = -cosa:                  Pa(i, 3) = Pa(i, 1) + 1
            Pa(i, 4) = -sinA:                  Pa(i, 5) = 2 * (eE(i) - nn) - 1:   Pa(i, 6) = cosa
            Pa(i, 7) = Pa(i, 5) + 1:           Pa(i, 8) = sinA:                   Pa(i, 9) = 100 * (ss - s(i))
            '权
            qL(i) = mM ^ 2 + (ss * pP * 0.0001) ^ 2
            
            S3 = "V(" & i & ")=" & Format(Pa(i, 2), "0.000") & "X(" & Pname(bE(i))
            S3 = S3 & ")+" & Format(Pa(i, 4), "0.000") & "Y(" & Pname(bE(i)) & ")+"
            S3 = S3 & Format(Pa(i, 6), "0.000") & "X(" & Pname(eE(i)) & ")+"
            S3 = S3 & Format(Pa(i, 8), "0.000") & "Y(" & Pname(eE(i)) & ")" & "+" & Format(Pa(i, 9), "0.00") & " , "
            For j = 1 To Len(S3)
                If Mid(S3, j, 1) = "+" And Mid(S3, j + 1, 1) = "-" Then S3 = Left$(S3, j - 1) & Right(S3, Len(S3) - j)
            Next j
            strR = strR & S3 & "权为：" & Format(qL(i), "0.000000") & vbCrLf & vbCrLf
        Next i
    End If
    
    If Nd > 0 Then '按方向计算方程系数
        For i = 1 To tn
            jj = 5
            
            Z0 = 0              '统计虚拟方程的常数项
            For j = Si(i) To Si(i) + Ni(i) - 1
                a0 = DirectAB(X0(eD(j)), Y0(eD(j)), X0(bD(j)), Y0(bD(j)))
                Z1 = a0 - Dir(j)
                If Z1 < 0 Then Z1 = Z1 + 2 * PI
                Z0 = Z0 + Z1
            Next j
            Z0 = Z0 / Ni(i)
            
            k2 = 1
            P3(i) = 0
            strR = strR & vbCrLf & "         方向误差方程：" & vbCrLf
            For j = Si(i) To Si(i) + Ni(i) - 1
                dX = X0(eD(j)) - X0(bD(j)):                      dY = Y0(eD(j)) - Y0(bD(j))
                a0 = DirectAB(X0(eD(j)), Y0(eD(j)), X0(bD(j)), Y0(bD(j)))
                ss = Sqr((X0(eD(j)) - X0(bD(j))) ^ 2 + (Y0(eD(j)) - Y0(bD(j))) ^ 2)
                ai = dY / (ss ^ 2) * RU / 100
                bi = -dX / (ss ^ 2) * RU / 100
                
                ii = Ne + j
                Pa(ii, 1) = 2 * (bD(j) - nn) - 1:                Pa(ii, 2) = ai
                Pa(ii, 3) = Pa(ii, 1) + 1:                       Pa(ii, 4) = bi
                Pa(ii, 5) = 2 * (eD(j) - nn) - 1:                Pa(ii, 6) = -ai
                Pa(ii, 7) = Pa(ii, 5) + 1:                       Pa(ii, 8) = -bi
                qL(ii) = mD ^ 2
                
                ss = Dir(j) + Z0:                                If ss >= 2 * PI Then ss = ss - 2 * PI
                Pa(ii, 9) = (a0 - ss) * RU
                
                Pa3(i, jj) = Pa(ii, 5):                          Pa3(i, jj + 1) = Pa(ii, 6) / qL(ii)
                Pa3(i, jj + 2) = Pa(ii, 7):                      Pa3(i, jj + 3) = Pa(ii, 8) / qL(ii)
                Pa3(i, 2) = Pa3(i, 2) + Pa(ii, 2) / qL(ii):      Pa3(i, 4) = Pa3(i, 4) + Pa(ii, 4) / qL(ii)
                P3(i) = P3(i) + Pa(ii, 9) / qL(ii)
                
                jj = jj + 4
                
                S3 = vbCrLf & "V(" & ii & ")=-" & Format(Z0, "0.000") & "+" & Format(Pa(ii, 2), "0.000") & "X(" & Pname(bD(j)) & ")+"
                S3 = S3 & Format(Pa(ii, 4), "0.000") & "Y(" & Pname(bD(j)) & ")+" & Format(Pa(ii, 6), "0.000") & "X(" & Pname(eD(j)) & ")+"
                S3 = S3 & Format(Pa(ii, 8), "0.000") & "Y(" & Pname(eD(j)) & ")" & "+" & Format(Pa(ii, 9), "0.000")
                S3 = S3 & " , 权为：" & Format(qL(ii), "0.000000") & vbCrLf
                For k = 1 To Len(S3)
                    If Mid(S3, k, 1) = "+" And Mid(S3, k + 1, 1) = "-" Then S3 = Left$(S3, k - 1) & Right(S3, Len(S3) - k)
                Next k
                strR = strR & S3
                pPp(k2) = eD(j)
                k2 = k2 + 1
            Next j
            
            '列虚拟误差方程
            k2 = k2 - 1
            Pa3(i, 1) = Pa(ii, 1):            Pa3(i, 3) = Pa(ii, 3):            qLS(i) = -Ni(i) / (mD ^ 2)
            '显示虚拟误差方程
            strR = strR & "     虚拟误差方程(测站号" & Pname(i) & ")：" & vbCrLf
            S3 = "V(" & Pname(i) & ")=" & Format(Pa3(i, 2), "0.000") & "X(" & Pname(i) & ")+" & Format(Pa3(i, 4), "0.000") & "Y(" & Pname(i) & ")+"
            For k = 6 To jj Step 4
                S3 = S3 & Format(Pa3(i, k), "0.000") & "X(" & Pname(pPp((k - 2) \ 4)) & ")+"
                S3 = S3 & Format(Pa3(i, 4), "0.000") & "Y(" & Pname(pPp((k - 2) \ 4)) & ")+"
            Next k
            S3 = S3 & Format(P3(i), "0.000")
            For k = 1 To Len(S3)
                If Mid(S3, k, 1) = "+" And Mid(S3, k + 1, 1) = "-" Then S3 = Left$(S3, k - 1) & Right(S3, Len(S3) - k)
            Next k
            strR = strR & S3 & " , 权为：" & Format(qL(i), "0.000000") & vbCrLf
        Next i
    End If
  
    Open App.Path & "\误差方程.txt" For Output As #1
        Print #1, strR
    Close #1
    txtShow.Text = txtShow.Text & vbCrLf & strR & vbCrLf
        
    '解算法方程：由系数矩阵Pa、Pa3，常数向量qL、qLS，求解＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
    Dim m%, n%, st%                     '误差方程总数，未知数总数，定向角参数个数
    Dim Pvv#, sigma#(100)               '用于统计单位权中误差的误差累计量
    Dim detS, detDir                    '边长改正数、方向改正数
    Dim t#(800, 800), R#(100, 100), c#(100, 100)      '临时的数组，辅助进行平差解算
    
    '形成法方程
    m = Ne + Nd         '误差方程总数
    For i = 1 To m
        For j = 1 To 4
            jj = Int(Pa(i, 2 * j - 1))
            If jj > 0 Then
                W(jj) = W(jj) + Pa(i, 2 * j) * Pa(i, 9) / qL(i)
                For k = 1 To 4
                    kk = Int(Pa(i, 2 * k - 1))
                    If kk > 0 And jj <= kk Then c(jj, kk) = c(jj, kk) + Pa(i, 2 * k) * Pa(i, 2 * j) / qL(i)
                Next k
            End If
        Next j
    Next i
    
    '和误差方程组成法方程
    For i = 1 To tn
        If Ni(i) <> 0 Then
            For j = 1 To 2 * (Ni(i) + 1)
                jj = Int(Pa3(i, 2 * j - 1))
                If jj > 0 Then
                    For k = 1 To 2 * (Ni(i) + 1)
                        kk = Int(Pa3(i, 2 * k - 1))
                        If kk >= 0 And jj <= kk Then c(jj, kk) = c(jj, kk) + Pa3(i, 2 * k) * Pa3(i, 2 * j) / qLS(i)
                    Next k
                End If
            Next j
        End If
    Next i
    st = 0              '统计定向角参数个数
    For i = 1 To tn
        If Ni(i) <> 0 Then st = st + 1
    Next i

    '求解法方程
    n = 2 * un              '未知数的个数：un个未知点，每个点x,y两个坐标
    For j = 1 To n
        t(1, j) = c(1, j)
    Next j
    For i = 2 To n
        For j = i To n
            t(i, j) = c(i, j)
            For k = 1 To i - 1
                t(i, j) = t(i, j) - t(k, i) * t(k, j) / t(k, k)
            Next k
        Next j
    Next i
    For i = 1 To n
        R(i, i) = 1 / t(i, i)
    Next i
    For i = 1 To n
        For j = i + 1 To n
            R(i, j) = 0
            For k = i To j - 1
                R(i, j) = R(i, j) - R(i, k) * t(k, j) / t(j, j)
            Next k
        Next j
    Next i
    For i = 1 To n
        For j = i To n
            For k = j To n
                Q(i, j) = t(k, k) * R(i, k) * R(j, k) + Q(i, j)
            Next k
            Q(j, i) = Q(i, j)
        Next j
    Next i

    '求改正数
    For i = 1 To n
        For j = 1 To n
            sigma(i) = sigma(i) - Q(i, j) * W(j)
        Next j
    Next i
    '改化坐标
    For i = nn + 1 To tn
        X(i) = X0(i) + sigma((i - nn) * 2 - 1) / 100
        Y(i) = Y0(i) + sigma((i - nn) * 2) / 100
    Next i
    
    '求单位权中误差
    For i = 1 To n
        Pvv = Pvv + sigma(i) * W(i)
    Next i
    For i = 1 To m
        Pvv = Pvv + Pa(i, 9) * Pa(i, 9) / qL(i)
    Next i
    uW0 = Sqr(Pvv / (m - n - st))
     
    '显示计算结果
    txtShow.Text = txtShow.Text & vbCrLf
    Open App.Path & "\平差结果.txt" For Output As #1
        Print #1, "    坐标计算结果(m)："
        Print #1, "点号    原X坐标       X坐标平差值      Vx        原Y坐标         Y坐标平差值       Vy"
        txtShow.Text = txtShow.Text & "    坐标计算结果(m)：" & vbCrLf
        txtShow = txtShow & "点号    原X坐标       X坐标平差值     Vx        原Y坐标        Y坐标平差值      Vy" & vbCrLf
        For i = nn + 1 To tn
            '给Y坐标加上带号和500Km
            Y(i) = Y(i) + N500 * 100000
            Y0(i) = Y0(i) + N500 * 100000
            Print #1, Pname(i), Format(X(i), "0.0000"), Format(Y(i), "0.0000")
            txtShow = txtShow & Pname(i) & "   " & Format(X0(i), "0.0000") & "   " & Format(X(i), "0.0000") & "   " & Format(X0(i) - X(i), "0.0000")
            txtShow = txtShow & "   " & Format(Y0(i), "0.0000") & "   " & Format(Y(i), "0.0000") & "   " & Format(Y0(i) - Y(i), "0.0000") & vbCrLf
        Next i
        Print #1, "σ0=" & Format(uW0, "0.00"); "厘米"
        txtShow.Text = txtShow.Text & "σ0=" & Format(uW0, "0.00") & "厘米" & vbCrLf
    Close #1
    
    Screen.MousePointer = 0
End Sub

'由方向值求夹角值：夹角的三个点，角顶点在中间，返回夹角值；参数
Public Function GetBeta(i1%, i2%, i3%, k%) As Double
    Dim k1%, k2%, dir1#, dir2#, i%
    k = 0: k1 = 0: k2 = 0
    For i = Si(i2) To Si(i2) + Ni(i2) - 1
        If eD(i) = i1 Then
            k1 = 1
            dir1 = Dir(i)
        End If
        If eD(i) = i3 Then
            k2 = 1
            dir2 = Dir(i)
        End If
    Next i
    If k1 = 1 And k2 = 1 Then
        k = 1
        GetBeta = dir2 - dir1
        If GetBeta < 0 Then GetBeta = GetBeta + 2 * PI
        If GetBeta > PI Then GetBeta = 2 * PI - GetBeta
    End If
End Function