平差核心代码

平差核心代码：

Const PI = 3.14159265358932

'求AB的坐标方位角，输入是两点坐标，输出的是弧度值

Public Function DirectAB(Xa#, Ya#, Xb#, Yb#) As Double

Dim detX#, detY#, tana#

detX = Xb - Xa

detY = Yb - Ya

If Abs(detX) < 0.000001 Then

If detY > 0 Then

DirectAB = PI / 2

Else

DirectAB = PI * 3 / 2

End If

Else

tana = detY / detX

DirectAB = Atn(tana)

If detX < 0 Then

DirectAB = PI + DirectAB

ElseIf detX > 0 And detY < 0 Then

DirectAB = PI * 2 + DirectAB

End If

End Function

'弧度化为度.分秒的形式：输入弧度值，输出度.分秒（各占两位）

Public Function HuToDo(ByVal Hu As Double) As Single

Dim du%, fen%, miao%

Hu = Hu * 180 / PI

du = Fix(Hu)

Hu = (Hu - du) * 60

fen = Fix(Hu)

Hu = (Hu - fen) * 60

miao = Fix(Hu + 0.5)

If miao = 60 Then

fen = fen + 1

miao = 0

End If

If fen = 60 Then

du = du + 1

fen = 0

End If

HuToDo = du + fen / 100 + miao / 10000

End Function

'将度.分秒形式化为弧度：输入为度.分秒形式，输出为弧度

Public Function DoToHu(ByVal DoFenMiao As Double) As Single

Dim du%, fen%, miao%, angle#

du = Fix(DoFenMiao)

DoFenMiao = (DoFenMiao - du) * 100

fen = Fix(DoFenMiao)

miao = (DoFenMiao - fen) * 100

angle = du + fen / 60 + miao / 3600

DoToHu = angle * PI / 180

End Function

'矩阵转置的通用过程

Public Sub MatrixTrans(A, c)

Dim i%, j%

Dim R1%, C1%

On Error Resume Next

C1 = UBound(A, 2) - LBound(A, 2) + 1

If Err Then

MsgBox "输入的矩阵维数不对！"

Exit Sub

End If

R1 = UBound(A, 1) - LBound(A, 1) + 1

ReDim c(1 To C1, 1 To R1)

For i = 1 To R1

For j = 1 To C1

c(j, i) = A(i, j)

Next j

Next i

End Sub

'矩阵相加的通用过程

Public Sub MatrixPlus(A, b, c)

Dim i%, j%

Dim R1%, C1%, R2%, C2%

On Error Resume Next

C1 = UBound(A, 2) - LBound(A, 2) + 1

If Err Then

MsgBox "第一个矩阵维数不对！"

Exit Sub

End If

On Error Resume Next

C2 = UBound(b, 2) - LBound(b, 2) + 1

If Err Then

MsgBox "第二个矩阵维数不对！"

Exit Sub

End If

R1 = UBound(A, 1) - LBound(A, 1) + 1

R2 = UBound(b, 1) - LBound(b, 1) + 1

If R1 <> R2 Or C1 <> C2 Then

MsgBox "输入的两个矩阵维数不等，不能相加！"

Exit Sub

End If

ReDim c(1 To m, 1 To n) As Double

For i = 1 To m

For j = 1 To n

c(i, j) = A(i, j) + b(i, j)

Next j

Next i

End Sub

'矩阵相减的通用过程

Public Sub MatrixMinus(A, b, c)

Dim i%, j%

Dim R1%, C1%, R2%, C2%

On Error Resume Next

C1 = UBound(A, 2) - LBound(A, 2) + 1

If Err Then

MsgBox "第一个矩阵维数不对！"

Exit Sub

End If

On Error Resume Next

C2 = UBound(b, 2) - LBound(b, 2) + 1

If Err Then

MsgBox "第二个矩阵维数不对！"

Exit Sub

End If

R1 = UBound(A, 1) - LBound(A, 1) + 1

R2 = UBound(b, 1) - LBound(b, 1) + 1

If R1 <> R2 Or C1 <> C2 Then

MsgBox "输入的两个矩阵维数不等，不能相减！"

Exit Sub

End If

ReDim c(1 To m, 1 To n) As Double

For i = 1 To m

For j = 1 To n

c(i, j) = A(i, j) - b(i, j)

Next j

Next i

End Sub

'矩阵相乘：输入矩阵或数Qa、Qb，自动识别它们的维数，并输出它们的乘积Qn

Public Sub Matrix_Multy(Qn, Qa, Qb)

Dim ia%, ib%, ic%

Dim ai%, bi%, ci%

Dim e1 As Boolean, e2 As Boolean, e3 As Boolean, e4 As Boolean, e5 As Boolean, e6 As Boolean, e7 As Boolean

On Error Resume Next '看Qa是不是一维数组

ic = UBound(Qa, 2) - LBound(Qa, 2)

If Err Then e1 = True

On Error Resume Next '看Qa是不是一维数组

ib = UBound(Qb, 2) - LBound(Qb, 2)

If Err Then e2 = True

If e1 = False And e2 = False Then '二维矩阵相乘

For ai = LBound(Qa, 1) To UBound(Qa, 1)

For bi = LBound(Qb, 2) To UBound(Qb, 2)

For ci = LBound(Qa, 2) To UBound(Qa, 2)

Qn(ai, bi) = Qn(ai, bi) + Qa(ai, ci) * Qb(ci, bi)

Next ci

Next bi

Next ai

ElseIf e1 = True And e2 = False Then

On Error Resume Next

ia = UBound(Qa) - LBound(Qa)

If Err Then e6 = True

If e6 Then '数乘以二维矩阵

For ai = LBound(Qb, 1) To UBound(Qb, 1)

For bi = LBound(Qb, 2) To UBound(Qb, 2)

Qn(ai, bi) = Qa * Qb(ai, bi)

Next bi

Next ai

Else '一维矩阵乘以二维矩阵

For ci = LBound(Qb, 2) To UBound(Qb, 2)

For ai = LBound(Qa, 1) To UBound(Qa, 1)

Qn(ci) = Qn(ci) + Qa(ai) * Qb(ai, ci)

Next ai

Next ci

End If

ElseIf e1 = False And e2 = True Then

On Error Resume Next

ic = UBound(Qb) - LBound(Qb)

If Err Then e7 = True

If e7 Then '二维矩阵乘以数

For ai = LBound(Qa, 1) To UBound(Qa, 1)

For bi = LBound(Qa, 2) To UBound(Qa, 2)

Qn(ai, bi) = Qa(ai, bi) * Qb

Next bi

Next ai

Else '二维矩阵乘以一维矩阵

For ai = LBound(Qa, 1) To UBound(Qa, 1)

For bi = LBound(Qa, 2) To UBound(Qa, 2)

Qn(ai) = Qn(ai) + Qa(ai, bi) * Qb(bi)

Next bi

Next ai

End If

Else

Dim errT As Integer

On Error Resume Next '结果是否是一个数

errT = UBound(Qn)

If Err Then e3 = True

If e3 Then '一维矩阵乘以一维矩阵得一个数

For ai = LBound(Qa, 1) To UBound(Qa, 1)

For bi = LBound(Qa, 2) To UBound(Qa, 2)

Qn = Qn + Qa(ai) * Qb(bi)

Next bi

Next ai

Exit Sub

End If

On Error Resume Next '是否是数乘一维矩阵

ia = UBound(Qa) - LBound(Qa)

If Err Then e4 = True

If e4 Then

For bi = LBound(Qa, 2) To UBound(Qa, 2)

Qn(bi) = Qa * Qb(bi)

Next bi

Exit Sub

End If

On Error Resume Next '是否是一维矩阵乘数

ib = UBound(Qb) - LBound(Qb)

If Err Then e5 = True

If e5 Then

For ai = LBound(Qa, 1) To UBound(Qa, 1)

Qn(ai) = Qa(ai) * Qb

Next ai

Exit Sub

End If

'一维矩阵相乘结果是二维矩阵

For ai = LBound(Qa, 1) To UBound(Qa, 1)

For bi = LBound(Qa, 2) To UBound(Qa, 2)

Qn(ai, bi) = Qa(ai) * Qb(bi)

Next bi

Next ai

End If

End Sub

'矩阵相乘的通用过程

Public Sub MatrixMulti(A, b, c)

Dim i%, j%, K%

Dim R1%, C1%, R2%, C2%

On Error Resume Next

C1 = UBound(A, 2) - LBound(A, 2) + 1

If Err Then

MsgBox "第一个矩阵维数不对！"

Exit Sub

End If

On Error Resume Next

C2 = UBound(b, 2) - LBound(b, 2) + 1

If Err Then

MsgBox "第二个矩阵维数不对！"

Exit Sub

End If

R1 = UBound(A, 1) - LBound(A, 1) + 1

R2 = UBound(b, 1) - LBound(b, 1) + 1

If C1 <> R2 Then

MsgBox "输入的两个矩阵大小不对，不能相乘！"

Exit Sub

End If

m = R1: s = C1: n = C2

ReDim c(1 To m, 1 To n) As Double

For i = 1 To m

For j = 1 To n

For K = 1 To s

c(i, j) = c(i, j) + A(i, K) * b(K, j)

Next K

Next j

Next i

End Sub

'列选主元法Guass约化求解线性方程组

Public Sub MajorInColGuass(A, b, X)

Dim Row%, Col%, n% '矩阵大小

Dim iStep%, iRow%, iCol% '循环变量

Dim L() As Double '各行的约化系数

'计算并检查矩阵的大小

Row = UBound(A, 1) - LBound(A, 1) + 1

Col = UBound(A, 2) - LBound(A, 2) + 1

If Row <> Col Then

MsgBox "方程组的系数矩阵有误！"

Exit Sub

End If

'准备约化过程的变量和数组

n = UBound(b) - LBound(b) + 1

If n <> Row Then

MsgBox "方程组的系数矩阵与常数项大小不符！"

Exit Sub

End If

ReDim L(2 To Row) As Double

Dim sumAX As Double, iPos%, temp#

'约化过程

For iStep = 1 To n - 1

'列选主元

iPos = 0

For iRow = iStep + 1 To n

If Abs(A(iRow, iStep)) > Abs(A(iStep, iStep)) Then

iPos = iRow

End If

Next iRow

If iPos > iStep Then '需要换主元

For iCol = iStep To n

temp = A(iStep, iCol)

A(iStep, iCol) = A(iPos, iCol)

A(iPos, iCol) = temp

Next iCol

temp = b(iStep)

b(iStep) = b(iPos)

b(iPos) = temp

End If

'约化过程

For iRow = iStep + 1 To n

L(iRow) = A(iRow, iStep) / A(iStep, iStep)

For iCol = iStep To n

A(iRow, iCol) = A(iRow, iCol) - L(iRow) * A(iStep, iCol)

Next iCol

b(iRow) = b(iRow) - L(iRow) * b(iStep)

Next iRow

ShowMatrix A

Next iStep

'回代过程

X(n) = b(n) / A(n, n)

For iRow = n - 1 To 1 Step -1

sumAX = 0

For iCol = n To iRow + 1 Step -1

sumAX = sumAX + A(iRow, iCol) * X(iCol)

Next iCol

X(iRow) = (b(iRow) - sumAX) / A(iRow, iRow)

Next iRow

End Sub

'Guass－Seidel迭代法求解线性方程组

Private Function Seidel(A, b, X, eps#) As Boolean

Dim i%, j%

Dim P#, Q#, s#, t#

Dim Row%, Col%, n%

Row = UBound(A, 1) - LBound(A, 1) + 1

Col = UBound(A, 2) - LBound(A, 2) + 1

n = UBound(b) - LBound(b) + 1

If n <> Row Then

MsgBox "方程组的系数矩阵与常数项大小不符！"

Exit Function

End If

For i = 1 To n

P = 0#

X(i) = 0#

For j = 1 To n

If i <> j Then P = P + Abs(A(i, j))

Next j

If P >= Abs(A(i, i)) Then

Seidel = False

Exit Function

End If

Next i

P = eps + 1#

While P >= eps

P = 0#

For i = 1 To n

t = X(i)

s = 0#

For j = 1 To n

If j <> i Then s = s + A(i, j) * X(j)

Next j

X(i) = (b(i) - s) / (A(i, i))

Q = Abs(X(i) - t) '/ (1# + Abs(x(i)))

If Q > P Then P = Q

Next i

Wend

Seidel = True

End Function

Public Sub ShowMatrix(tt)

Dim i%, j%, n%, m%

m = UBound(tt, 1) - LBound(tt, 1) + 1

n = UBound(tt, 2) - LBound(tt, 2) + 1

For i = 1 To m

For j = 1 To n

Debug.Print tt(i, j),

Next j

Debug.Print

Next i

End Sub

'通用的间接平差解算过程：输入系数矩阵A、权矩阵P、常数向量L和解向量X，求出X，并通过参数传出去

Public Sub InAdjust(A, P, L, X)

Dim a1%, a2%, p1%, p2%, L1%, x1% '输入矩阵或向量的大小

Dim At() As Double, AtP() As Double, Naa#(), W() As Double '几个中间矩阵

'计算并检查输入矩阵或向量的大小

On Error Resume Next

a1 = UBound(A, 1) - LBound(A, 1) + 1

If Err Then

MsgBox "系数矩阵A大小错误！"

Exit Sub

End If

On Error Resume Next

a2 = UBound(A, 2) - LBound(A, 2) + 1

If Err Then

MsgBox "系数矩阵A大小错误！"

Exit Sub

End If

On Error Resume Next

L1 = UBound(L) - LBound(L) + 1

If Err Then

MsgBox "常数向量L大小错误！"

Exit Sub

End If

On Error Resume Next

x1 = UBound(X) - LBound(X) + 1

If Err Then

MsgBox "解向量X大小错误！"

Exit Sub

End If

On Error Resume Next

p1 = UBound(P, 1) - LBound(P, 1) + 1

If Err Then

MsgBox "权矩阵P大小错误！"

Exit Sub

End If

On Error Resume Next

p2 = UBound(P, 2) - LBound(P, 2) + 1

If Err Then

MsgBox "权矩阵P大小错误！"

Exit Sub

End If

If p1 <> p2 Then

MsgBox "权矩阵P不是方阵！"

Exit Sub

End If

If p1 <> a1 Or p2 <> a1 Then

MsgBox "权矩阵P与系数矩阵A大小不符！"

Exit Sub

End If

If a2 <> x1 Then

MsgBox "系数矩阵A大小与解向量X大小不符！"

Exit Sub

End If

If a1 <> L1 Then

MsgBox "系数矩阵A大小与常数向量L大小不符！"

Exit Sub

End If

'定义中间矩阵的大小

ReDim At(1 To a2, 1 To a1), AtP(1 To a2, 1 To a1)

ReDim Naa(1 To a2, 1 To a2), W(1 To a2)

'组成法方程并计算

Debug.Print "The A matrix is:"

ShowMatrix A

MatrixTrans A, At '求A的转置矩阵

Debug.Print "The At matrix is:"

ShowMatrix At

Debug.Print "The P matrix is:"

ShowMatrix P

Matrix_Multy AtP, At, P '求AtP

Debug.Print "and The AtP matrix is:"

ShowMatrix AtP

Matrix_Multy Naa, AtP, A '法方程系数矩阵

Debug.Print "the Naa matrix is:"

ShowMatrix Naa

Debug.Print "the L matrix is:"

For x1 = LBound(L) To UBound(L)

Debug.Print L(x1)

Next x1

Matrix_Multy W, AtP, L '法方程常数向量

Debug.Print "the W matrix is:"

For x1 = LBound(W) To UBound(W)

Debug.Print W(x1)

Next x1

MajorInColGuass Naa, W, X

Debug.Print "the X matrix is:"

For x1 = LBound(X) To UBound(X)

Debug.Print X(x1)

Next x1

'Seidel Naa, W, x, 0.000001

End Sub

'通用的条件平差解算过程：输入系数矩阵A、权矩阵P、常数向量L和解向量X，求出X，并通过参数传出去

Public Sub CondiAdjust(b, P, W, V)

Dim b1%, b2%, p1%, p2%, w1%, v1% '输入矩阵或向量的大小

Dim Q#(), Bt#(), QBt#(), Nbb#(), K#(), i% '几个中间矩阵

'计算并检查输入矩阵或向量的大小

On Error Resume Next

b1 = UBound(b, 1) - LBound(b, 1) + 1

If Err Then

MsgBox "系数矩阵B大小错误！"

Exit Sub

End If

On Error Resume Next

b2 = UBound(b, 2) - LBound(b, 2) + 1

If Err Then

MsgBox "系数矩阵B大小错误！"

Exit Sub

End If

On Error Resume Next

w1 = UBound(W) - LBound(W) + 1

If Err Then

MsgBox "常数向量W大小错误！"

Exit Sub

End If

On Error Resume Next

v1 = UBound(V) - LBound(V) + 1

If Err Then

MsgBox "改正数向量V大小错误！"

Exit Sub

End If

On Error Resume Next

p1 = UBound(P, 1) - LBound(P, 1) + 1

If Err Then

MsgBox "权矩阵P大小错误！"

Exit Sub

End If

On Error Resume Next

p2 = UBound(P, 2) - LBound(P, 2) + 1

If Err Then

MsgBox "权矩阵P大小错误！"

Exit Sub

End If

If p1 <> p2 Then

MsgBox "权矩阵P不是方阵！"

Exit Sub

End If

If p1 <> b2 Then

MsgBox "权矩阵P与系数矩阵A大小不符！"

Exit Sub

End If

If b2 <> v1 Then

MsgBox "系数矩阵B大小与解向量V大小不符！"

Exit Sub

End If

If b1 <> w1 Then

MsgBox "系数矩阵B大小与常数向量W大小不符！"

Exit Sub

End If

'定义中间矩阵的大小

ReDim Bt(1 To b2, 1 To b1), QBt(1 To b2, 1 To b1)

ReDim Nbb(1 To b1, 1 To b1), K(1 To b1), Q(1 To p1, 1 To p2)

'组成法方程并计算

For i = 1 To p1 '求Q矩阵

Q(i, i) = 1 / P(i, i)

Next i

MatrixTrans b, Bt

Matrix_Multy QBt, Q, Bt

Matrix_Multy Nbb, b, QBt '法方程系数矩阵

ShowMatrix Nbb

MajorInColGuass Nbb, W, K '解法方程

'Seidel Nbb, W, K, 0.0000001

Matrix_Multy V, QBt, K '求改正数

End Sub

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。