'***********************************************************************************
'    本程序用于單一附（閉）合導線嚴密平差計算，采用按角度條件平差法。以方向觀測值中
'誤差的先驗值作為單位權中誤差。計算結果可求得各待定點的坐標平差值及其點位精度Mx，My，
'及M，并計算出各待定點誤差橢圓元素E，F，Z
'
'參考文獻： 郭久訓．《控制網平差程序設計》　北京：原子能出版社，2004．8
'
'平差數據來源：潘正風等．《數字測圖原理與方法》　武漢大學出版社　　186頁表6-5
'
'等級：所用平差數據為首級圖根導線數據（精度很不高）。本程序中寫出了方位角和導線全長相對閉合差的判
'別，但考慮到程序的通用性，將這些限差判斷當作了注釋處理，而不實際運行。
'
'程序不足：沒有導線網的圖形表達。
'***********************************************************************************

Private i%, j%, n%, tc#, tb#, B_x!, B_y!, e1!, e2!, m!, m0#, z#, aa As Boolean, bb As Boolean, cc As Boolean ',dd As Boolean
Private Naa#(2, 2), Naa逆#(2, 2), W#(2), K#(2), qq#(2), fx#(2), fy#(2)
Private A#(), Q#(), V#(), C#(), mx#(), my#(), mk#(), e#(), f#(), zz#()
Private Po() As Point





'文件格式說明：
'文件格式詳見文件"平差數據.txt"
'
Private Sub 打開文件_Click()                  '打開文件
    Dim ff$, temp$, A_name$, A_x!, A_y!, A_l#, A_s!, B_name$, B_l#, C_name$, C_x!, C_y!, D_name$, D_x!, D_y!
    Form1.Cls  '清屏
    CommonDialog1.DialogTitle = "打開數據文件"
    CommonDialog1.FileName = ""
    CommonDialog1.ShowOpen
    '出錯處理
    On Error GoTo FileErr
        ff = CommonDialog1.FileName           'ff是文件路徑名
        Open ff For Input As #1               '以順序文件方式打開文件，使用input
        
        Line Input #1, temp                   '讀取文件中的說明語句
        Line Input #1, temp                   '讀取文件中的說明語句
        Line Input #1, temp                   '讀取文件中的說明語句
        Input #1, n                           '讀取n，n為 （測站數-1）
        ReDim Po(n) As Point                  '定義Po(n)，其中Po(0)存A點數據，Po(n)存B點數據，Po(1)到Po(n-1)存n-1個未知點數據。
        Line Input #1, temp                   '讀取文件中的說明語句
        Line Input #1, temp                   '讀取文件中的說明語句
        Input #1, m, e1, e2                   '讀取先驗方向觀測值中誤差m，測距儀固定誤差e1，比例誤差e2
        Line Input #1, temp                   '讀取文件中的說明語句
        Line Input #1, temp                   '讀取文件中的說明語句
        Input #1, C_name, C_x, C_y            '讀取已知點C
        Input #1, A_name, A_x, A_y, A_l, A_s  '讀取已知點A
        
        For i = 1 To n - 1                    '讀取n-1個未知點
            Input #1, Po(i).name, Po(i).l, Po(i).s
        Next i
        
        Input #1, B_name, B_x, B_y, B_l       '讀取已知點B
        Input #1, D_name, D_x, D_y            '讀取已知點D
                      
        Close #1
        
        tc = ZBiaoFSuan(C_x, C_y, A_x, A_y)    '坐標反算,求點C到A的坐標方位角，并記作tc，單位是度
        tb = ZBiaoFSuan(B_x, B_y, D_x, D_y)    '坐標反算,求點B到D的坐標方位角，并記作tb，單位是度
        
        Po(0).name = A_name: Po(0).x = A_x: Po(0).y = A_y: Po(0).l = A_l: Po(0).s = A_s  'Po(0)存A點數據
        Po(n).name = B_name: Po(n).x = B_x: Po(n).y = B_y: Po(n).l = B_l                 'Po(n)存B點數據
        For i = 0 To n
            Po(i).l = deg(Po(i).l)             '將觀測方向左角的單位度分秒化作度
        Next i
        
        ReDim Q#(2 * n), A#(2, 2 * n), V#(2 * n), C#(1, 2 * n - 2), mx#(n - 1), my#(n - 1), _
        mk#(n - 1), e#(n - 1), f#(n - 1), zz#(n - 1)        '變量重新定義
        bb = True
        MsgBox "文件已成功打開", , "提示"
        顯示平差數據 ff
        Exit Sub
FileErr:
        MsgBox "您的文件未打開或打開的文件格式有誤！注意：請重新運行本程序！！", , "提示"
End Sub

'
'
'
Private Sub 開始計算_Click()
    
    If bb = True Then
        推算方位角和坐標         '求近似坐標
        求條件式
        'If dd = True Then        '若平差數據沒有超限，則進行下面的計算
        組法方程式
        求逆求K
        求改正數和平差值
        推算方位角和坐標         '求平差后坐標
        精度評定
        'Else
            'MsgBox "限差超限！", , "提示"
            'Exit Sub
        'End If
    Else
        MsgBox "數據文件未打開！！", , "提示"
        Exit Sub
    End If
    aa = True
    MsgBox "計算完畢！", , "提示"
    顯示平差結果
End Sub

'
'
'
Private Sub 保存_Click()
    Dim ff$
    
    If aa = True Then
        CommonDialog1.DialogTitle = "保存平差結果"
        CommonDialog1.Filter = "(*.txt)|*.txt"
        CommonDialog1.FileName = "平差結果.txt"
        CommonDialog1.ShowSave
        ff = CommonDialog1.FileName
        If ff = "" Then
            MsgBox "文件名不能為空", , "警告"
            Exit Sub
        End If
          
        Open ff For Output As #2                    '以順序文件方式保存文件，使用output
            
            Print #2, "平差結果："
            Print #2, "---------------------------------------------------------------------------------------------"
            Print #2, " 單位權中誤差 m0=" & Format(m0, "####.##") & " 秒"
            Print #2, "---------------------------------------------------------------------------------------------"
            Print #2, "點名", "x(m)", "y(m)", "點位誤差(cm)", "E(cm)", "F(cm)", "Z(度分秒)"
            For i = 1 To n - 1
                Print #2, " " & Po(i).name, Format(Po(i).x, "########.####"), Format(Po(i).y, "########.####"), "     " & Format(mk(i), "#####.##"), _
                        "      " & Format(e(i), "#####.##"), "      " & Format(f(i), "#####.##"), "      " & Format(zz(i), "###.#####")
            Next i
            Print #2, "---------------------------------------------------------------------------------------------"
        Close #2
          
        cc = True
        MsgBox "保存完畢！", , "提示"
    Else
        MsgBox "沒有需要保存平差結果！", , "提示"
        Exit Sub
    End If
End Sub

'
'
'
Private Sub 退出_Click()
    If aa = True And cc = False Then
            提示.Show                         '若平差結果未保存則進行提示
    Else
            End                               '若沒有平差結果或平差結果已保存則退出
    End If
End Sub

'
'
'
Private Sub 推算方位角和坐標()
    Po(0).t = tc + Po(0).l - 180
    If Po(0).t > 360 Then Po(0).t = Po(0).t - 360
    If Po(0).t < 0 Then Po(0).t = Po(0).t + 360
    
    For i = 1 To n
        '求方位角
        Po(i).t = Po(i - 1).t + Po(i).l - 180
        If Po(i).t > 360 Then Po(i).t = Po(i).t - 360
        If Po(i).t < 0 Then Po(i).t = Po(i).t + 360
        '求坐標
        Po(i).x = Po(i - 1).x + Po(i - 1).s * Cos(Po(i - 1).t / p0)   't要以度為單位
        Po(i).y = Po(i - 1).y + Po(i - 1).s * Sin(Po(i - 1).t / p0)
    Next i
End Sub

'
'角度改正數的單位是秒，邊長改正數的單位是厘米
'方向觀測值中誤差為單位權
'方向觀測值中誤差的先驗值為 m(秒)
'角度觀測值的權為 0.5
'距離觀測值中誤差的先驗值為 ms=e1+e2*s*0.0001(厘米)
'
Private Sub 求條件式()
    Dim sums!
    
    '求A
    For i = 0 To n       '求角度改正數系數
        A(0, n + i) = 1
        A(1, n + i) = (Po(i).y - Po(n).y) / 2062.65
        A(2, n + i) = (Po(n).x - Po(i).x) / 2062.65
        Q(n + i) = 2     '求角度觀測值的權倒數
    Next i
    For i = 0 To n - 1   '求邊長改正數系數
        A(1, i) = Cos(Po(i).t / p0)
        A(2, i) = Sin(Po(i).t / p0)
        Q(i) = (e1 + e2 * Po(i).s * 0.0001) ^ 2 / m ^ 2 '求距離觀測值的權倒數
    Next i
    '求W
    W(0) = (Po(n).t - tb) * 3600  '單位是s
    W(1) = (Po(n).x - B_x) * 100  '單位是cm
    W(2) = (Po(n).y - B_y) * 100  '單位是cm
    For i = 0 To n - 1
        sums = sums + Po(i).s
    Next i
   
    '限差判斷
    '這是首級圖根導線的精度要求，其中 n+1 為測站數
    'If Abs(W(0)) > 40 * Sqr(n + 1) Or Sqr((W(1) / 100) ^ 2 + (W(2) / 100) ^ 2) / sums > 0.00025 Then '方位角閉合差和導線全長相對閉合差
        'Exit Sub
    'End If
    'dd = True
End Sub

'
'
'
Private Sub 組法方程式()
    Dim g%
    For i = 0 To 2
        For j = i To 2        '由于Naa是對稱的，只求其上三角即可。
            For g = 0 To 2 * n
                Naa(i, j) = Naa(i, j) + A(i, g) * A(j, g) * Q(g)
            Next g
        Next j
    Next i
End Sub

'
'
'對于3*3的矩陣來說，利用公式　Naa逆=Naa*/|Naa| 　去求Naa的逆即可，其中Naa*是伴隨矩陣。
Private Sub 求逆求K()
    Dim ee#, ii#, jj#, kk#, det#       'det是|Naa|的值
    ee = 1
    For i = 0 To 2
        ee = ee * Naa(i, i)
    Next i
    ee = ee + 2 * Naa(0, 1) * Naa(1, 2) * Naa(0, 2)
    ii = Naa(0, 1) ^ 2
    jj = Naa(0, 2) ^ 2
    kk = Naa(1, 2) ^ 2
    det = ee - (ii * Naa(2, 2) + jj * Naa(1, 1) + kk * Naa(0, 0))
    '求伴隨矩陣Naa*，用Naa逆存儲
    Naa逆(0, 0) = Naa(1, 1) * Naa(2, 2) - kk
    Naa逆(0, 1) = Naa(0, 2) * Naa(1, 2) - Naa(0, 1) * Naa(2, 2)
    Naa逆(1, 0) = Naa逆(0, 1)
    Naa逆(1, 1) = Naa(0, 0) * Naa(2, 2) - jj
    Naa逆(0, 2) = Naa(0, 1) * Naa(1, 2) - Naa(0, 2) * Naa(1, 1)
    Naa逆(2, 0) = Naa逆(0, 2)
    Naa逆(2, 2) = Naa(0, 0) * Naa(1, 1) - ii
    Naa逆(1, 2) = Naa(0, 2) * Naa(0, 1) - Naa(0, 0) * Naa(1, 2)
    Naa逆(2, 1) = Naa逆(1, 2)
    '求Naa逆
    For i = 0 To 2
        For j = 0 To 2
            Naa逆(i, j) = Naa逆(i, j) / det
        Next j
    Next i
    '求K
    For i = 0 To 2
        For j = 0 To 2
            K(i) = K(i) - Naa逆(i, j) * W(j)
        Next j
    Next i
End Sub

'
'
'
'
Private Sub 求改正數和平差值()
    Dim pvv#
    pvv = 0
    For i = 0 To 2 * n
        For j = 0 To 2
            V(i) = V(i) + K(j) * A(j, i) * Q(i)     '求改正數V
        Next j
        pvv = pvv + V(i) * V(i) / Q(i)              '求 pvv
    Next i
    m0 = Sqr(pvv / 3)                               '求單位權中誤差 m0
    For i = 0 To n - 1
        Po(i).s = Po(i).s + V(i) / 100              '求導線邊長的平差值，其中邊長改正數的單位是厘米，邊長的單位是米
    Next i
    For i = 0 To n
        Po(i).l = Po(i).l + V(n + i) / 3600         '求導線左角的平差值，其中角度改正數的單位是秒，角度的單位是度
    Next i
End Sub

'
'列出待定點坐標的權函數式，據此進行精度評定。
'
Private Sub 精度評定()
    Dim uu#, vv#, ww#, ss#, r#, kk!
    For j = 1 To n - 1
        uu = 0: vv = 0: ww = 0
        For i = 0 To j - 1
            C(0, i) = Cos(Po(i).t / p0)
            C(1, i) = Sin(Po(i).t / p0)
            C(0, n + i) = (Po(i).y - Po(j).y) / 2062.65
            C(1, n + i) = (Po(j).x - Po(i).x) / 2062.65
            uu = uu + C(0, i) ^ 2 * Q(i) + C(0, n + i) ^ 2 * Q(n + i)
            vv = vv + C(1, i) ^ 2 * Q(i) + C(1, n + i) ^ 2 * Q(n + i)
            ww = ww + C(0, i) * C(1, i) * Q(i) + C(0, n + i) * C(1, n + i) * Q(n + i)
        Next i
        For kk = 0 To 2
            ss = 0: r = 0
            For i = 0 To j - 1
                ss = ss + A(kk, i) * C(0, i) * Q(i) + A(kk, n + i) * C(0, n + i) * Q(n + i)
                r = r + A(kk, i) * C(1, i) * Q(i) + A(kk, i + n) * C(1, n + i) * Q(n + i)
            Next i
            fx(kk) = ss: fy(kk) = r
        Next kk
        For kk = 0 To 2
            ss = 0: r = 0
            For i = 0 To 2
                ss = ss - Naa逆(kk, i) * fx(i)
                r = r - Naa逆(kk, i) * fy(i)
            Next i
            K(kk) = ss: qq(kk) = r
        Next kk
        For i = 0 To 2
            uu = uu + fx(i) * K(i)
            vv = vv + fy(i) * qq(i)
            ww = ww + fx(i) * qq(i)
        Next i
        mx(j) = m0 * Sqr(Abs(uu))
        my(j) = m0 * Sqr(Abs(vv))
        mk(j) = Sqr((mx(j)) ^ 2 + (my(j)) ^ 2)        '每個點的點位中誤差，單位cm
        r = uu + vv
        ss = Sqr((uu - vv) ^ 2 + 4 * ww ^ 2)
        e(j) = m0 * Sqr((r + ss) / 2)                 '點位誤差橢圓元素E，單位cm
        f(j) = m0 * Sqr(Abs(r - ss) / 2)              '點位誤差橢圓元素F，單位cm
        z = Atn(2 * ww / (uu - vv))
        z = z / 2 * p0
        If z < 0 Then z = z + 90
        If ww < 0 Then z = z + 180
        zz(j) = dms(z)                                '點位誤差橢圓元素Z
    Next j
End Sub

'
'
'
Private Function 顯示平差數據(ByVal fname$)
    Dim linesfromFile$, Nextline$
    
    Open fname For Input As #1
        Do Until EOF(1)
            Line Input #1, Nextline
            linesfromFile = linesfromFile + Nextline + Chr(13)
        Loop
    Close #1
    Print linesfromFile
End Function

'
'
'
Private Sub 顯示平差結果()
    Print "平差結果："
    Print "---------------------------------------------------------------------------------------------"
    Print " 單位權中誤差 m0=" & Format(m0, "####.##") & " 秒"
    Print "---------------------------------------------------------------------------------------------"
    Print "點名", "x(m)", "y(m)", "點位誤差(cm)", "E(cm)", "F(cm)", "Z(度分秒)"
    For i = 1 To n - 1
        Print " " & Po(i).name, Format(Po(i).x, "########.####"), Format(Po(i).y, "########.####"), "  " & Format(mk(i), "#####.##"), _
                  " " & Format(e(i), "#####.##"), Format(f(i), "#####.##"), Format(zz(i), "###.#####")
    Next i
    Print "---------------------------------------------------------------------------------------------"
    '
    '
    '
    '可以將平差所得的B點坐標和B點的實際坐標輸出進行比較，以驗證平差結果是否正確
    'Print Format(Po(n).x, "########.######"), Format(Po(n).y, "########.######")
    'Print B_x, B_y
End Sub