平差数据处理.bas

一个附和导线的严密计算平差程序

Attribute VB_Name = "modJiSuan"

Public Type Shurushuju              '输入数据记录类型
  ICount As Integer
  PiontName As String
  Guancejiao As Double
  Bianchang As Double
  x As Double
  y As Double
End Type
Public Type ZBJieGuo               '坐标值和坐标精度类型
  x As Double                      '平差后的x坐标
  y As Double                      '平差后的y坐标
  m As Single                      '平差后的坐标精度
  mx As Single                     '平差后的x坐标方向精度
  my As Single                     '平差后的y坐标方向精度
  E As Single                      '点误差椭圆长轴方向长度
  F As Single                      '点误差椭圆短轴方向长度
  Q As Double                      '点误差椭圆长轴方位角
  a As Double
End Type
Public Type WZDJTY           '求未知点间椭圆元素、方位角和边长误差类型
  E As Single
  F As Single
  Q As Double
  s As Double
  a As Double
  Ms As Single
  ma As Single
  MsS As Long
End Type
Public InputDat(100) As Shurushuju
Public ZBJG(100) As ZBJieGuo
Public Wzdty(100) As WZDJTY            '求点间椭圆元素、方位角和边长误差及相对误差
Public mm As Single, ma As Single, mb As Single  '观测方向中误差和边长丈量中误差系数
Public m0 As Single                 '计算单位权中误差
Public aa(4) As Double                    '起始方位角和末边方位角
Public a0(100) As Double                 '近似方位角
Public StrFilename As String              '打开文件名及路径
Public xx(4) As Double, yy(4) As Double                  '始末点的坐标
Public PointHao(4) As Integer                       '点号
Public PointName(4) As String         '点名
Public Const RuoC = 206.265          '若C常数
Public IntCount As Integer          '记录的个数
Public Ps(100) As Single        '观测边长的权
Public Pb(100) As Single        '观测角的权
Public Wa As Single             '方位角闭合差
Public Wx As Single             'X坐标闭合差
Public Wy As Single             'Y坐标闭合差
Public X0(100) As Double, Y0(100) As Double     '近似坐标
Public a(100, 100) As Single               '条件方程系数
Public Q(100, 100) As Single                'Q阵
Public n(10, 10) As Single                'N阵
Public nnn()  As Double                   'N的逆矩阵
Public QLL(100, 100) As Single            '观测值的权阵
Public Vaxy(100, 100) As Single           '点位改正数系数
Public QFF(100, 100) As Single            '点位协因数阵
Public Qas(100, 100) As Single            '方位角和边长中误差
Public K(10) As Single                    'K阵
Public v(100) As Single                   'V为观测角和观测边长改正数
Public W(10) As Single                    'W为条件方程的闭合差
Public Va(100) As Single                  'Va为方位角改正数
Public Vs(100) As Single                  '边长改正数
Public s(100) As Double                   '平差后的边长
Public Vx(100) As Single, Vy(100) As Single    '纵横坐标改正数
Public i As Integer, j As Integer, m As Integer    '循环变量
Public Const PI = 3.141592654             'PI为圆周率
Public Function ZuobiaoFangwei(ByVal x1 As Double, ByVal y1 As Double, ByVal x2 As Double, ByVal y2 As Double) As Double  '由坐标反求方位角，方位角为弧度
  Dim Dx As Double, Dy As Double
  Dim T As Double
  Dx = x2 - x1
  Dy = y2 - y1
  If Abs(Dx) < Abs(Dy) Then
    T = Atn(Abs(Dx) / Abs(Dy))
    T = 0.5 * PI - T
  End If
  If Abs(Dx) > Abs(Dy) Then
    T = Atn(Abs(Dy) / Abs(Dx))
  End If
  If Dx >= 0 And Dy > 0 Then
 T = Atn(Abs(Dx) / Abs(Dy))
    T = 0.5 * PI - T
    End If
  If Dx >= 0 And Dy < 0 Then
    T = 2 * PI - T
  End If
  If Dx < 0 And Dy < 0 Then
    T = PI + T
  End If
  If Dx < 0 And Dy > 0 Then
    T = PI - T
  End If
  ZuobiaoFangwei = T
End Function

Public Function DFMDu(ByVal x As Double) As Double      '度分秒化为度
Dim x1 As Integer, x2 As Integer
Dim x3 As Double
  x1 = Int(x)
  x = (x - x1) * 100
  x2 = Int(x)
  x = (x - x2) * 100
  x3 = x1 + (x2 * 60 + x) / 3600
  DFMDu = x3
End Function

Public Function DuHudu(ByVal x As Double) As Double      '度化为弧度
Dim Hu As Double
  Hu = x * PI / 180
  DuHudu = Hu
End Function

Public Function HuduDu(ByVal x As Double) As Double      '弧度化为度
  HuduDu = x * 180 / PI
End Function

Public Function DuDFM(ByVal x As Double) As Double       '度化为度分秒
Dim x1 As Integer, x2 As Integer
  x1 = Int(x)
  x = x - x1
  x = x * 3600
  x2 = Int(x / 60)
  DuDFM = x1 + x2 / 100 + (x - x2 * 60) / 10000
End Function

Public Function DuM(ByVal x As Double) As Double
 DuM = x * 3600
End Function
Public Sub jz_con(a() As Double)                            '矩阵求逆
Dim b() As Double, temp As Double, tt As Double, c() As Double
Dim F As Boolean, jszb() As Double
Dim n As Integer, cir As Integer, T As Integer
Dim i%, j%, K%
n = UBound(a)
ReDim b(n, n), c(n, 2 * n + 1)
' 定义P
For i = 0 To n
      For j = 0 To n
       c(i, j) = a(i, j)
       c(i, n + j + 1) = 0
     Next j
        c(i, n + i + 1) = 1
Next i
For i = 0 To n
    For j = 0 To 2 * n + 1
     ' Text2.Text = Text2.Text & c(i, j) & "     "
      Next j
     'Text2.Text = Text2.Text + vbCrLf
     Next i
 For i = 0 To n
       K = i
      Do While c(K, i) = 0 And K < n And Not F
           K = K + 1
         Loop
       If c(K, i) = 0 Then
           K = K + 1
           End If
      If K <= n Then
         F = True
         Else
         F = False
        Text1.Text = Text1.Text + "没有逆矩阵"
      '  Picture1.Print "没有逆矩阵"
       ' text2.Text = "no"
         End If               '判0
      If K <= n And F Then
        For j = i To 2 * n + 1
        temp = c(i, j)
        c(i, j) = c(K, j)
        c(K, j) = temp
        Next j
        For cir = i To n
          If Abs(c(cir, i)) >= 0.000001 Then
          temp = c(cir, i)
       For j = i To 2 * n + 1
       c(cir, j) = c(cir, j) / temp   '变为 1
        Next j
        End If
        Next cir
       End If
        
    If F And i < n Then
       For j = i + 1 To n
       temp = c(j, i)
       If Abs(temp) >= 0.000001 Then
         For cir = i To 2 * n + 1
          c(j, cir) = c(j, cir) - c(i, cir)
           Next cir
           End If
           Next j
           End If
      Next i                 '上三角距阵
     If F Then
        For i = 0 To n - 1
            If i < n Then
          For j = i + 1 To n
             temp = c(i, j)
            For cir = j To 2 * n + 1
               c(i, cir) = c(i, cir) - temp * c(j, cir)
             Next cir
             Next j
             End If
             Next i
      End If
   For i = 0 To n
      For j = 0 To n
         a(i, j) = c(i, n + 1 + j)
         If Abs(a(i, j)) <= 0.000001 Then
         a(i, j) = 0
         End If
         Next j
         Next i
    For i = 0 To n
    For j = 0 To 2 * n + 1
'      Text2.Text = Text2.Text & c(i, j) & "    "
      Next j
'     Text2.Text = Text2.Text + vbCrLf
     Next i
End Sub

Public Sub QiuDianweiXS(Vxy() As Single)              '求点位改正数系数
 For j = 1 To IntCount * 2
   For i = 1 To IntCount * 2 - 1
     Vxy(j, i) = 0
   Next i
 Next j
 For j = 3 To (IntCount - 1) * 2 Step 2
   m = j \ 2 + 1
   For i = 1 To m - 1
    Vxy(j, i) = -(Y0(m) - Y0(i)) / RuoC
    Vxy(j, i + IntCount) = Cos(DuHudu(a0(i)))
    Vxy(j + 1, i) = (X0(m) - X0(i)) / RuoC
    Vxy(j + 1, i + IntCount) = Sin(DuHudu(a0(i)))
   Next i
 Next j
End Sub


Public Sub QiuQLLzhen(QL() As Single)                  '观测值权阵
  Dim i As Integer, j As Integer, m As Integer
  Dim Nk(100, 100) As Single
  Dim Nkk(100, 100) As Single
  
  For i = 1 To IntCount * 2 - 1
    For j = 1 To 3
      Nk(i, j) = 0
     For m = 1 To IntCount * 2 - 1
      Nk(i, j) = Nk(i, j) + Q(i, m) * a(j, m)
     Next m
    Next j
  Next i
  For i = 1 To IntCount * 2 - 1
    For j = 1 To 3
     Nkk(i, j) = 0
      For m = 1 To 3
        Nkk(i, j) = Nkk(i, j) + Nk(i, m) * nnn(m - 1, j - 1)
      Next m
    Next j
  Next i
  For i = 1 To IntCount * 2 - 1
    For j = 1 To IntCount * 2 - 1
      Nk(i, j) = 0
      For m = 1 To 3
        Nk(i, j) = Nk(i, j) + Nkk(i, m) * a(m, j)
      Next m
    Next j
  Next i
  For i = 1 To IntCount * 2 - 1
    For j = 1 To IntCount * 2 - 1
     Nkk(i, j) = 0
       For m = 1 To IntCount * 2 - 1
         Nkk(i, j) = Nkk(i, j) + Nk(i, m) * Q(m, j)
       Next m
    Next j
  Next i
  For i = 1 To IntCount * 2 - 1
    For j = 1 To IntCount * 2 - 1
      QL(i, j) = Q(i, j) - Nkk(i, j)
    Next j
  Next i
End Sub

Public Sub QiuDianweiQuan(Qf() As Single)
                                           '求点位改正数协因素阵
Dim i As Integer
Dim j As Integer
Dim m As Integer
Dim Nk(100, 100) As Single
  For i = 1 To IntCount * 2
    For j = 1 To IntCount * 2 - 1
     Nk(i, j) = 0
     For m = 1 To IntCount * 2 - 1
     Nk(i, j) = Nk(i, j) + Vaxy(i, m) * QLL(m, j)
     Next m
    Next j
  Next i
  For i = 1 To IntCount * 2
   For j = 1 To IntCount * 2
    Qf(i, j) = 0
    For m = 1 To IntCount * 2 - 1
     Qf(i, j) = Qf(i, j) + Nk(i, m) * Vaxy(j, m)
    Next m
   Next j
  Next i
End Sub

Public Sub QiuDianweiJD(ZBJD() As ZBJieGuo)           '求点位精度
 For i = 1 To IntCount
   With ZBJD(i)
     .mx = m0 * Sqr(QFF(i * 2 - 1, i * 2 - 1))
     .my = m0 * Sqr(QFF(i * 2, i * 2))
     .m = Sqr(.mx * .mx + .my * .my)
   End With
 Next i
End Sub

Public Function TuoyuanZhuidazhiFWJ(Qxx As Single, Qyy As Single, Qxy As Single) As Double
                                              '求误差椭圆长半轴方位角
  TuoyuanZhuidazhiFWJ = DuDFM(HuduDu(ZuobiaoFangwei(Qyy, 0, Qxx, 2 * Qxy) / 2))
End Function

Public Sub QiuDianTuoyuan(DianTY() As ZBJieGuo)
                                              '求点位误差椭圆元素
Dim d As Single
 For i = 2 To IntCount - 1
   d = Tan(DuHudu(DFMDu(TuoyuanZhuidazhiFWJ(QFF(i * 2 - 1, i * 2 - 1), QFF(i * 2, i * 2), QFF(i * 2 - 1, i * 2)))))
   With DianTY(i)
        .E = m0 * Sqr(QFF(i * 2 - 1, i * 2 - 1) + QFF(i * 2 - 1, i * 2) * d)
        .F = m0 * Sqr(QFF(i * 2, i * 2) - QFF(i * 2 - 1, i * 2) * d)
        .Q = TuoyuanZhuidazhiFWJ(QFF(i * 2 - 1, i * 2 - 1), QFF(i * 2, i * 2), QFF(i * 2 - 1, i * 2))
      End With
 Next i
End Sub

Public Sub QiuWZDJTY(DianjianTY() As WZDJTY)
                                              '求未知点间误差椭圆元素的过程
 Dim Qdxx As Single
 Dim Qdyy As Single
 Dim Qdxy As Single
 Dim d As Single
  For i = 1 To IntCount - 1
    Qdxx = QFF(i * 2 - 1, i * 2 - 1) + QFF((i + 1) * 2 - 1, (i + 1) * 2 - 1) - 2 * QFF(i * 2 - 1, (i + 1) * 2 - 1)
    Qdyy = QFF(i * 2, i * 2) + QFF((i + 1) * 2, (i + 1) * 2) - 2 * QFF(i * 2, (i + 1) * 2)
    Qdxy = QFF(i * 2 - 1, i * 2) - QFF((i + 1) * 2 - 1, i * 2) - QFF(i * 2 - 1, (i + 1) * 2) + QFF((i + 1) * 2 - 1, (i + 1) * 2)
    d = Tan(DuHudu(DFMDu(TuoyuanZhuidazhiFWJ(Qdxx, Qdyy, Qdxy))))
    With DianjianTY(i)
        .E = m0 * Sqr(Qdxx + Qdxy * d)
        .F = m0 * Sqr(Qdyy - Qdxy * d)
        .Q = TuoyuanZhuidazhiFWJ(Qdxx, Qdyy, Qdxy)
   End With
  Next i
End Sub

Public Sub QiuBianHeFangweiJiaoJD(Qasf() As Single)           '求方位角和边长改正数Q阵
 Dim a() As Single               '方位角和边长改正数系数矩阵
 Dim Nk() As Single
 ReDim a(IntCount * 2 - 2, IntCount * 2 - 1)
 For i = 1 To IntCount * 2 - 2
  For j = 1 To IntCount * 2 - 1
    If i < IntCount Then
      If i >= j Then
        a(i, j) = 1
        Else
        a(i, j) = 0
       End If
       Else
       If i = (j - 1) Then
         a(i, j) = 1
        Else
         a(i, j) = 0
       End If
    End If
  Next j
 Next i
 
 ReDim Nk(IntCount * 2 - 1, IntCount * 2 - 1)
  For i = 1 To IntCount * 2 - 2
    For j = 1 To IntCount * 2 - 1
     Nk(i, j) = 0
     For m = 1 To IntCount * 2 - 1
     Nk(i, j) = Nk(i, j) + a(i, m) * QLL(m, j)
     Next m
    Next j
  Next i
  For i = 1 To IntCount * 2 - 2
   For j = 1 To IntCount * 2 - 2
    Qasf(i, j) = 0
    For m = 1 To IntCount * 2 - 1
    Qasf(i, j) = Qasf(i, j) + Nk(i, m) * a(j, m)
    Next m
   Next j
  Next i
End Sub

Public Sub QiuMsMa()
 For i = 1 To IntCount - 1
  With Wzdty(i)
       .s = s(i)
       .a = DuDFM(ZBJG(i).a)
       .ma = m0 * Sqr(Qas(i, i))
       .Ms = m0 * Sqr(Qas((i + IntCount - 1), (i + IntCount - 1)))
       .MsS = Int(s(i) * 1000 / .Ms)
  End With
 Next i
End Sub