现在用VB.net的人也很少了，程序都是本科的时候写的，属于自学成才，很怀念当年熬夜写代码debug的自己，没人看也无所谓，算是给自己留个纪念。
1. 窗口界面
这个就个人风格很强烈了，我喜欢简洁的。
2.Form中的代码
Imports System.IO
Public Class Form1
'*************************************************读取.o文件************************************************
Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
OpenFileDialog1.InitialDirectory = "C:\Users\dell\Desktop\2019工作学习\GPS课程设计"
OpenFileDialog1.Filter = "所有文件(*.*)|*.*"
OpenFileDialog1.FilterIndex = 1
OpenFileDialog1.RestoreDirectory = True
If OpenFileDialog1.ShowDialog = DialogResult.OK Then
TextFileName1.Text = OpenFileDialog1.FileName
End If
If TextFileName1.Text <> "" Then
Dim fs As FileStream
fs = New FileStream(TextFileName1.Text, FileMode.OpenOrCreate, FileAccess.Read)
Dim sr As StreamReader
sr = New StreamReader(fs)
sum1 = sr.ReadToEnd
sr.Close()
fs.Close()
End If
Dim sum11(), sum12(), sum13(), sum14(), sum15() As String
sum11 = Split(Trim(sum1), vbLf)
Dim num, num2 As Integer
Dim i, j, k As Integer
'读取卫星的近似位置X0,Y0,Z0
num = UBound(sum11)
sum12 = Split(Trim(sum11(8)), " ")
X0 = CDbl(CDec(Trim(sum12(0))))
Y0 = CDbl(CDec(Trim(sum12(1))))
Z0 = CDbl(CDec(Trim(sum12(2))))
'读取时间
sum13 = Split(Trim(sum11(74)), " ")
' 读取年月日时分秒
Dim prn As String
nian = CInt(sum13(0))
yue = CInt(sum13(1))
ri = CInt(sum13(2))
ho = CInt(sum13(3))
mi = CInt(sum13(4))
se = CInt(sum13(5))
'计算观测时刻瞬时历元 周内秒
'计算这天是星期几
Dim yc As Integer
Dim mc As Integer
Dim fc As Integer
If IsRunnian(nian) Then
mc = Val(Mid("512503514624", yue, 1))
Else
mc = Val(Mid("622503514624", yue, 1))
End If
yc = ((nian Mod 100) \ 4 + (nian Mod 100)) Mod 7
fc = (yc + mc + ri) Mod 7
t = fc * 24 * 3600 + ho * 3600 + mi * 60 + se
'读取序列号PRN（为方便区分省略字母），本次读取10个
prn = sum13(6)
For i = 0 To 9
sum14 = Split(Trim(prn), "G")
prnn(i) = sum14(1 + i)
Next
'读取伪距row
For j = 0 To 9
sum15 = Split(Trim(sum11(76 + j * 2)), " ")
row(j) = CDbl(Replace(sum15(0), " ", ""))
Next
End Sub
'*************************************************读取.n文件************************************************
Private Sub Button3_Click(sender As Object, e As EventArgs) Handles Button3.Click
OpenFileDialog1.InitialDirectory = "C:\Users\dell\Desktop\2019工作学习\GPS课程设计"
OpenFileDialog1.Filter = "所有文件(*.*)|*.*"
OpenFileDialog1.FilterIndex = 1
OpenFileDialog1.RestoreDirectory = True
If OpenFileDialog1.ShowDialog = DialogResult.OK Then
TextFileName2.Text = OpenFileDialog1.FileName
End If
If TextFileName2.Text <> "" Then
Dim fs As FileStream
fs = New FileStream(TextFileName2.Text, FileMode.OpenOrCreate, FileAccess.Read)
Dim sr As StreamReader
sr = New StreamReader(fs)
sum2 = sr.ReadToEnd
sr.Close()
fs.Close()
End If
Dim sum21(), sum22(), sum23(), sum24(), sum25(), sum26, sum27, sum28, sum29, sum30, sum31, sum32 As String
sum21 = Split(Trim(sum2), vbLf)
Dim num, num2 As Integer
Dim i, j, k As Integer
'读取 ALPHA ,BETA
sum22 = Split(Trim(sum21(3)), " ")
sum23 = Split(Trim(sum21(4)), " ")
For i = 0 To 3
ALPHA(i) = CDec(sum22(i)) '将科学计数法转换成十进制
BETA(i) = CDec(sum23(i))
Next
'********************************读取导航电文等参数*******************************
'a0(9), a1(9), a2(9), IODE(9), Crs(9), dn(9), M0(9), Cuc(9), ee(9), Cus(9), sqrtA(9), TOE(9), Cic(9), omig(9), Cis(9), i0(9), Crc(9), w(9), omigp(9), II(9)
'将从O文件中读取的序列号与N文件匹配，然后读取对应的参数
For i = 0 To 9
For j = 8 To 87
sum26 = Trim(sum21(j))
sum24 = Split(Trim(sum21(j)), " ")
num = Len(sum26)
If sum24(0) = prnn(i) Then
year(i) = sum24(1)
mon(i) = sum24(2)
day(i) = sum24(3)
shi(i) = sum24(4)
fen(i) = sum24(5)
sum25 = Split(Trim(sum24(7)), ".")
miao(i) = sum25(0)
a0(i) = CDec(Trim(Mid(sum26, 23, 19)))
a1(i) = CDec(Trim(Mid(sum26, 42, 19)))
a2(i) = CDec(Trim(Mid(sum26, 61, 19)))
sum27 = sum21(j + 1)
IODE(i) = CDec(Trim(Mid(sum27, 4, 19)))
Crs(i) = CDec(Trim(Mid(sum27, 23, 19)))
dn(i) = CDec(Trim(Mid(sum27, 42, 19)))
M0(i) = CDec(Trim(Mid(sum27, 61, 19)))
sum28 = sum21(j + 2)
Cuc(i) = CDec(Trim(Mid(sum28, 4, 19)))
ee(i) = CDec(Trim(Mid(sum28, 23, 19)))
Cus(i) = CDec(Trim(Mid(sum28, 42, 19)))
sqrtA(i) = CDec(Trim(Mid(sum28, 61, 19)))
sum29 = sum21(j + 3)
TOE(i) = CDec(Trim(Mid(sum29, 4, 19)))
Cic(i) = CDec(Trim(Mid(sum29, 23, 19)))
omig(i) = CDec(Trim(Mid(sum29, 42, 19)))
Cis(i) = CDec(Trim(Mid(sum29, 61, 19)))
sum30 = sum21(j + 4)
TOE(i) = CDec(Trim(Mid(sum30, 4, 19)))
Cic(i) = CDec(Trim(Mid(sum30, 23, 19)))
omig(i) = CDec(Trim(Mid(sum30, 42, 19)))
Cis(i) = CDec(Trim(Mid(sum30, 61, 19)))
sum31 = sum21(j + 5)
i0(i) = CDec(Trim(Mid(sum31, 4, 19)))
Crc(i) = CDec(Trim(Mid(sum31, 23, 19)))
w(i) = CDec(Trim(Mid(sum31, 42, 19))) / 10000
omigp(i) = CDec(Trim(Mid(sum31, 61, 19)))
sum32 = sum21(j + 6)
II(i) = CDec(Trim(Mid(sum32, 4, 19)))
End If
Next
Next
End Sub
Private Sub Button2_Click(sender As Object, e As EventArgs) Handles Button2.Click
'************************************************计算卫星位置***********************************************
Dim i, j, q As Integer
Dim jE(9) As Double '偏近点角E
Dim a, n0, n, tk, M, E1, E2, de As Double
Dim f, u1, du, dr, di, u, r, i2, L As Double
Dim Xs, Ys, Zs As Double '接收机坐标
Dim Xs1(9), Ys1(9), Zs1(9) As Double '卫星在协议地球坐标系位置
Dim x1, y1, z1 As Double '卫星在轨道平面坐标系位置
Dim X0, Y0, Z0 As Double
X0 = -2148744.2864 '卫星坐标XYZ近似值
Y0 = 4426639.1875
Z0 = 4044655.3434
For i = 0 To 9 '10个卫星 10次循环
'根据广播星历计算卫星无摄运动参数
a = Math.Pow(sqrtA(i), 2)
n0 = Math.Sqrt(GM / Math.Pow(a, 3))
n = n0 + dn(i)
tk = t - TOE(i)
M = M0(i) + n * tk
'迭代计算偏近点角E
E1 = M
Do
E2 = M + ee(i) * Math.Sin(E1)
de = E2 - E1
E1 = E2
Loop While (Math.Abs(de) < 0.00000001) '取绝对值
jE(i) = E2 '赋值偏近点角E
f = Math.Atan2((Math.Sqrt(1 - Math.Pow(ee(i), 2)) * Math.Sin(E2)), (Math.Cos(E2) - ee(i)))
If ((Math.Sqrt(1 - Math.Pow(ee(i), 2)) * Math.Sin(E2)) < 0) Then
f = 2 * Math.PI - Math.Abs(f)
End If
u1 = w(i) + f
'计算摄动改正项
du = Cuc(i) * Math.Cos(2 * u1) + Cus(i) * Math.Sin(2 * u1)
dr = Crc(i) * Math.Cos(2 * u1) + Crs(i) * Math.Sin(2 * u1)
di = Cic(i) * Math.Cos(2 * u1) + Cis(i) * Math.Sin(2 * u1)
u = u1 + du
r = a * (1 - ee(i) * Math.Cos(E2)) + dr
i2 = i0(i) + di + II(i) * tk
'计算卫星在轨道平面坐标系中的位置
x1 = r * Math.Cos(u)
y1 = r * Math.Sin(u)
'计算卫星在协议地球坐标系中的位置
L = omig(i) + omigp(i) * tk - we * t
Xs1(i) = x1 * Math.Cos(L) - y1 * Math.Cos(i2) * Math.Sin(L)
Ys1(i) = x1 * Math.Sin(L) + y1 * Math.Cos(i2) * Math.Cos(L)
Zs1(i) = y1 * Math.Sin(i2)
Next
'************************************************计算电离层延迟改正***********************************************
Dim idl(9) As Double
Dim b = 0.000000005
Dim s = 14 * 3600.0
Dim el(9) As Double '卫星在测站P处的截止高度角 , 单位度
For i = 0 To 9
el(i) = Math.Atan2(Math.Abs(Zs1(i) - Z0), Math.Sqrt(Math.Pow(Xs1(i) - X0, 2) + Math.Pow(Ys1(i) - Y0, 2)))
el(i) = el(i) * 180.0 / Math.PI '弧度转度
Next
'测站的地心经纬度
Dim lp, faip As Double
lp = 0
faip = 0
Dim c = 6399593.6258
Dim EE1, EE2 As Double
EE1 = 0.00669437999013
EE2 = 0.00673949674227
lp = Math.Acos(X0 / Math.Sqrt(Math.Pow(X0, 2) + Math.Pow(Y0, 2))) * 180.0 / Math.PI
'迭代计算faip
Dim t0, P, K, t2, t1, dt As Double
t0 = Z0 / Math.Sqrt(Math.Pow(X0, 2) + Math.Pow(Y0, 2))
P = (c * EE1) / Math.Sqrt(Math.Pow(X0, 2) + Math.Pow(Y0, 2))
K = 1 + E2
t1 = t0
Do
t2 = t0 + (P * t1) / Math.Sqrt(K + Math.Pow(t1, 2))
dt = t2 - t1
t1 = t2
Loop While (dt < 0.0000001)
faip = (Math.Atan(t2) * 180.0) / Math.PI
Dim alf(9) As Double
For i = 0 To 9
alf(i) = Math.Atan(Math.Abs(Ys1(i)) / Math.Abs(Xs1(i)))
If (Xs1(i) < 0 And Ys1(i) > 0) Then
alf(i) = Math.PI - alf(i)
ElseIf (Xs1(i) < 0 And Ys1(i) < 0) Then
alf(i) = Math.PI + alf(i)
ElseIf (Xs1(i) > 0 And Ys1(i) < 0) Then
alf(i) = 2 * Math.PI - alf(i)
End If
Next
Dim UT As Double
Dim seita(9), lp1(9), faip1(9), tcc(9), faim(9), AMP(9), PER(9), Fqx(9), secZ(9) As Double
UT = ho '观测时刻的世界时 单位小时
' seita 测站与穿刺点对应的地心夹角 单位度
' lp1 穿刺点地心经纬度
' tcc 观测瞬间穿刺点的地方时 单位秒
' faim 穿刺点的地磁纬度 单位度
' AMP 振幅 PER 周期 fqx 倾斜校正因子 secZ 天顶距
'***************************************************计算电离层延迟 Idl**********************************************
For i = 0 To 9
seita(i) = 445.0 / (el(i) + 20.0) - 4
lp1(i) = lp + seita(i) * Math.Sin(alf(i)) / Math.Cos(faip * Math.PI / 180.0)
faip1(i) = faip + seita(i) * Math.Cos(alf(i))
tcc(i) = (UT + lp1(i) / 15.0) * 3600.0
faim(i) = faip1(i) + (90 - 79.93) * Math.Cos((faip1(i) - 288.04) * Math.PI / 180.0)
faim(i) = faim(i) * Math.PI / 180.0
AMP(i) = ALPHA(0) + ALPHA(1) * faim(i) + ALPHA(2) * Math.Pow(faim(i), 2) + ALPHA(3) * Math.Pow(faim(i), 3)
PER(i) = BETA(0) + BETA(1) * faim(i) + BETA(2) * Math.Pow(faim(i), 2) + BETA(3) * Math.Pow(faim(i), 3)
Fqx(i) = 1.0 + 16.0 * Math.Pow(0.53 - el(i) * Math.PI / 180.0, 3)
idl(i) = Fqx(i) * (b + AMP(i) * Math.Cos(2 * Math.PI * (tcc(i) - s) / PER(i)))
secZ(i) = 1 + 2 * Math.Pow((96 - el(i)) / 90, 3)
idl(i) = idl(i) * secZ(i)
idl(i) = idl(i) * c '光速
Next
'************************************************计算卫星钟差改正VtS***********************************************
Dim VtS(9) As Double
For i = 0 To 9
VtS(i) = a0(i) + a1(i) * (t - TOE(i)) + a2(i) * (t - TOE(i)) * (t - TOE(i))
Next
'************************************************* 平差计算************************************************
Dim sigema, mx, my, mz, num As Double
Dim temp(0, 0) As Double
Dim row0(9), LL(9), MM(9), NN(9) As Double
Dim X01, Y01, Z01 As Double
Dim Bb(9, 3), xx(3, 0), lf(9, 0), V(9, 0), PP(9, 9) As Double
Dim BT(3, 9), BTPB(3, 3), BTPL(3, 0), BTP(3, 9), BX(9, 0), VT(0, 9), VTp(0, 9), qx(9, 9), BTB(3, 3) As Double
num = 10 '平差计算迭代次数
sigema = 0 ' 单位权中误差
mx = 0
My = 0
mz = 0 ' 坐标中误差
For q = 0 To 9
For i = 0 To 9
row0(i) = Math.Sqrt(Math.Pow(Xs1(i) - X0, 2) + Math.Pow(Ys1(i) - Y0, 2) + Math.Pow(Zs1(i) - Z0, 2))
LL(i) = (Xs1(i) - X0) / row0(i)
MM(i) = (Ys1(i) - Y0) / row0(i)
NN(i) = (Zs1(i) - Z0) / row0(i)
Bb(i, 0) = -LL(i)
Bb(i, 1) = -MM(i)
Bb(i, 2) = -NN(i)
Bb(i, 3) = 1
lf(i, 0) = (row(i) - row0(i) + c * VtS(i) - idl(i)) / 100000
Next
'定权
For i = 0 To 9
For j = 0 To 9
If (i <> j) Then
PP(i, j) = 0
ElseIf (jE(i) < 60) Then
PP(i, j) = Math.Pow(Math.Sin(el(i) * Math.PI / 180.0), 2)
Else
PP(i, j) = 1
End If
Next
Next
'矩阵运算
BT = ZZ(Bb, BT) 'BB的转置
BTPB = Mmul(Mmul(BT, PP, BTP), Bb, BTPB)
BTPL = Mmul(BTP, lf, BTPL)
xx = Mmul(QN(BTPB), BTPL, xx)
X01 = X0 + xx(0, 0)
Y01 = Y0 + xx(1, 0)
Z01 = Z0 + xx(2, 0)
Next
'接收机坐标 赋值
Xs = X01
Ys = Y01
Zs = Z01
'精度评定
V = QH(Mmul(Bb, xx, BX), FZ(lf))
VT = ZZ(V, VT)
temp = Mmul(Mmul(VT, PP, VTp), V, temp)
sigema = Math.Sqrt(temp(0, 0) / (9 - 4)) '单位权中误差
'协因数矩阵
qx = QN(Mmul(BT, Bb, BTB))
mx = sigema * Math.Sqrt(qx(0, 0))
my = sigema * Math.Sqrt(qx(1, 1))
mz = sigema * Math.Sqrt(qx(2, 2))
'输出
TextX.Text = CStr(Xs)
TextY.Text = CStr(Ys)
TextZ.Text = CStr(Zs)
TextX1.Text = CStr(mx)
TextY1.Text = CStr(my)
TextZ1.Text = CStr(mz)
wucha.Text = CStr(sigema)
End Sub
Private Sub Form1_Load(sender As Object, e As EventArgs) Handles MyBase.Load
End Sub
End Class
Module中的代码
模块就是用来写全局变量或者函数的，很必要。
Module Module1
'*********************************************声明部分*************************************************
Public GM = 398600500000000.0 '万有引力常数G与地球总质量M的乘积
Public we = 0.0000729211567 ' 地球自转角速度
Public sum1 As String
Public sum2 As String
Public sj As String
Public t As Double '观测时刻瞬时历元 周内秒
Public nian, yue, ri, ho, mi, se As Integer
Public prn As String
Public prnn(9) As String ' 序列号
Public row(9) As Double ' 伪距
Public ALPHA(3) As Double
Public BETA(3) As Double
Public year(9), mon(9), day(9), shi(9), fen(9), miao(9) As Double
Public X0 As Double '卫星坐标XYZ近似值
Public Y0 As Double
Public Z0 As Double
Public a0(9), a1(9), a2(9), IODE(9), Crs(9), dn(9), M0(9), Cuc(9), ee(9), Cus(9), sqrtA(9), TOE(9), Cic(9), omig(9), Cis(9), i0(9), Crc(9), w(9), omigp(9), II(9) As Double
'******************************************函数部分**************************************************
'将string转换为矩阵
Public Function Str2Mit(str As String)
Dim mat(,) As Double, matRows() As String, matCells() As String
Dim i, j, k As Integer
matRows = str.Trim().Split(vbCrLf)
matCells = matRows(0).Trim().Split(New Char() {" ", ",", "，"})
ReDim mat(matRows.Length - 1, matCells.Length - 1)
For i = 0 To matRows.Length - 1
For j = 0 To matCells.Length - 1
matCells = matRows(i).Trim().Split(New Char() {" ", ",", "，"})
mat(i, j) = matCells(j)
Next
Next
Return mat
End Function
'矩阵相乘函数
Function Mmul(mtxA(,) As Double, mtxB(,) As Double, ByRef mtxC(,) As Double)
Dim m As Integer
Dim n As Integer
Dim l As Integer
Dim i As Integer, j As Integer, K As Integer
m = UBound(mtxA, 1) - LBound(mtxA, 1)
n = UBound(mtxA, 2) - LBound(mtxA, 2)
l = UBound(mtxB, 2) - LBound(mtxB, 2)
For i = 0 To m
For j = 0 To l
mtxC(i, j) = 0#
For K = 0 To n
mtxC(i, j) = mtxC(i, j) + mtxA(i, K) * mtxB(K, j)
Next K
Next j
Next i
Return mtxC
End Function
'矩阵转置函数
Function ZZ(mtxA(,) As Double, mtxB(,) As Double)
Dim m As Integer
Dim n As Integer
Dim i As Integer, j As Integer
m = UBound(mtxA, 1) - LBound(mtxA, 1) '原来的行数
n = UBound(mtxA, 2) - LBound(mtxA, 2) '原来的列数
For i = 0 To m
For j = 0 To n
mtxB(j, i) = mtxA(i, j)
Next j
Next i
Return mtxB
End Function
'矩阵求逆函数
Function QN(Mtx(,) As Double)
Dim N As Double
N = UBound(Mtx, 1) - LBound(Mtx, 1) '矩阵的行数/列数
Dim row(0 To N) As Integer, col(0 To N) As Integer
Dim d As Double, p As Double
For k = 0 To N
d = 0.0#
For i = k To N
For j = k To N
p = Math.Abs(Mtx(i, j))
If (p > d) Then
d = p
row(k) = i
col(k) = j '记录变换的位置
End If
Next j
Next i
If (d + 1.0# = 1.0#) Then '找到主元，并存储到行列号
End If
If (row(k) <> k) Then
For j = 0 To N
p = Mtx(k, j)
Mtx(k, j) = Mtx(row(k), j)
Mtx(row(k), j) = p '行变换
Next j
End If
If (col(k) <> k) Then
For i = 0 To N
p = Mtx(i, k)
Mtx(i, k) = Mtx(i, col(k))
Mtx(i, col(k)) = p '列变换
Next i
End If
Mtx(k, k) = 1.0# / Mtx(k, k)
For j = 0 To N
If (j <> k) Then
Mtx(k, j) = Mtx(k, j) * Mtx(k, k)
End If
Next j
For i = 0 To N
If (i <> k) Then
For j = 0 To N
If (j <> k) Then
Mtx(i, j) = Mtx(i, j) - Mtx(i, k) * Mtx(k, j)
End If
Next j
End If
Next i
For i = 0 To N
If (i <> k) Then
Mtx(i, k) = -Mtx(i, k) * Mtx(k, k)
End If
Next i
Next k
For k = N To 0 Step -1
If (col(k) <> k) Then
For j = 0 To N
p = Mtx(k, j)
Mtx(k, j) = Mtx(col(k), j)
Mtx(col(k), j) = p
Next j
End If
If (row(k) <> k) Then
For i = 0 To N
p = Mtx(i, k)
Mtx(i, k) = Mtx(i, row(k))
Mtx(i, row(k)) = p
Next i
End If
Next k
Return Mtx
End Function
'求一个矩阵的负值
Function FZ(Mtx(,) As Double)
Dim m, n As Integer
Dim C(,) As Double
Dim I, J As Integer
m = UBound(Mtx, 1) - LBound(Mtx, 1)
n = UBound(Mtx, 2) - LBound(Mtx, 2)
ReDim C(m, n)
For I = 0 To m
For J = 0 To n
C(I, J) = -Mtx(I, J)
Next
Next
Return C
End Function
'矩阵相加
Function QH(mtxA(,) As Double, mtxB(,) As Double)
Dim m, n As Double
Dim SUM(,) As Double
Dim I, J As Integer
m = UBound(mtxA, 1) - LBound(mtxB, 1)
n = UBound(mtxA, 2) - LBound(mtxB, 2)
ReDim SUM(m, n)
For I = 0 To m
For J = 0 To n
SUM(I, J) = mtxA(I, J) + mtxB(I, J)
Next
Next
Return sum
End Function
Public Function IsRunnian(ByVal n As Integer) As Boolean
Dim res As Boolean
res = False
If n Mod 400 = 0 Then
res = True
End If
If (n Mod 100 <> 0) And (n Mod 4 = 0) Then
res = True
End If
IsRunnian = res
End Function
End Module
一看就知道是测绘人！希望学弟学妹们就算看到了也别直接抄，最起码可以换换变量名，哈哈。
可以的话还是建议用C#，VB写界面方便，但为了个人编程能力的提高还是多用一些主流的语言，如C++，java，Python
扫一扫
4479
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
