module1.bas

本程序所反应的是在测量中被广泛应用的全球定位系统（Global Position System）在空间中的形象以及基本的工作原理。GPS由三部分组成：空间部分、地面控制系统以及用户部分。

Attribute VB_Name = "Module1"
Public Const pi = 3.14159365358, radiu = 6370, orbit = 20200 + radiu
Public xe As Double, ye As Double, ze As Double, sx As Double, sy As Double
Public sita1 As Double, fai1 As Double, rou1 As Double, d1 As Double, zzj As Double
Public stltqishi As Single '控制卫星起始角
Public stltx(1 To 6, 1 To 4) As Double, stlty(1 To 6, 1 To 4) As Double, stltz(1 To 6, 1 To 4) As Double  '放24颗卫星的三维坐标
Public stltx1(1 To 6, 1 To 4) As Double, stlty1(1 To 6, 1 To 4) As Double, stltz1(1 To 6, 1 To 4) As Double  '放24颗卫星的三维坐标
Public stltx2(1 To 6, 1 To 4) As Double, stlty2(1 To 6, 1 To 4) As Double          ' 存放24颗卫星的2维坐标
Public stlti(1 To 4) As Double, stltj(1 To 4) As Double                '存放卫星的下标
Public eyex As Double, eyey As Double, eyez As Double
Public x3d(1 To 2) As Double, y3d(1 To 2) As Double, z3d(1 To 2) As Double
Public x2d(1 To 2) As Double, y2d(1 To 2) As Double
Public linx As Double, liny As Double, linz As Double
Public qingjiao As Double
Public ymin(-4000 To 4000) As Double, ymax(-4000 To 4000) As Double    '用于隐藏线的消除
Public maxx As Double, maxy As Double     '用于存放画地球时最大的sx和sy
Public minx As Double, miny As Double     '用于存放画地球时最小的sx和sy
Public qishi1(1 To 24)               '为了使卫星运行时不相撞,而使之开始的起始角不同
Public locationx As Double, locationy As Double, locationz As Double   '存放点击的地面点的三维坐标
Public pic1control As Integer   '控制是否点出了地面点
Public picx As Double, picy As Double         '存放地面点的屏幕2维坐标
Public Sub drawearth(sita As Double, fai As Double, rou As Double, d As Double, zzhuanjiao As Double)
 qingjiao = 10
 'For i = -4000 To 4000
     'ymin(i) = 5000
     'ymax(i) = -5000
 'Next i
 maxx = -4000
 maxy = -4000
 minx = 4000
 miny = 4000
 '画坐标轴
 
 '画不动的坐标轴
 
 x2d(1) = zbchangex(8000, 0, 0, sita, fai, rou, d)
 y2d(1) = zbchangey(8000, 0, 0, sita, fai, rou, d)
 Form1.Picture1.DrawStyle = 0
 For i = 0 To 50 Step 10
    Form1.Picture1.Circle (x2d(1), y2d(1)), i, vbWhite
 Next i
 Form1.Picture1.DrawStyle = 2
 Form1.Picture1.Line (0, 0)-(x2d(1), y2d(1)), vbWhite
 x2d(1) = zbchangex(0, 8000, 0, sita, fai, rou, d)
 y2d(1) = zbchangey(0, 8000, 0, sita, fai, rou, d)
 Form1.Picture1.DrawStyle = 0
 For i = 0 To 50 Step 10
    Form1.Picture1.Circle (x2d(1), y2d(1)), i, vbWhite
 Next i
 Form1.Picture1.DrawStyle = 2
 Form1.Picture1.Line (0, 0)-(x2d(1), y2d(1)), vbWhite
 x2d(1) = zbchangex(0, 0, 8000, sita, fai, rou, d)
 y2d(1) = zbchangey(0, 0, 8000, sita, fai, rou, d)
 Form1.Picture1.DrawStyle = 0
 For i = 0 To 50 Step 10
    Form1.Picture1.Circle (x2d(1), y2d(1)), i, vbWhite
 Next i
 Form1.Picture1.DrawStyle = 2
 Form1.Picture1.Line (0, 0)-(x2d(1), y2d(1)), vbWhite
 
 '画地固坐标系X
 Form1.Picture1.DrawStyle = 0
 linx = 8000
 liny = 0
 linz = 0
 x3d(1) = linx * Cos(zzhuanjiao * pi / 180) - liny * Sin(zzhuanjiao * pi / 180)
 y3d(1) = linx * Sin(zzhuanjiao * pi / 180) + liny * Cos(zzhuanjiao * pi / 180)
 z3d(1) = linz

 linx = x3d(1)
 liny = y3d(1)
 linz = z3d(1)
 x3d(1) = linx
 y3d(1) = liny * Cos(qingjiao * pi / 180) - linz * Sin(qingjiao * pi / 180)
 z3d(1) = liny * Sin(qingjiao * pi / 180) + linz * Cos(qingjiao * pi / 180)

 
     
     
 x2d(1) = zbchangex(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
 y2d(1) = zbchangey(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
 For i = 0 To 50 Step 10
    Form1.Picture1.Circle (x2d(1), y2d(1)), i, vbGreen
 Next i
 Form1.Picture1.Line (0, 0)-(x2d(1), y2d(1)), vbGreen
 
 '画地固坐标系Y
 linx = 0
 liny = 8000
 linz = 0
 x3d(1) = linx * Cos(zzhuanjiao * pi / 180) - liny * Sin(zzhuanjiao * pi / 180)
 y3d(1) = linx * Sin(zzhuanjiao * pi / 180) + liny * Cos(zzhuanjiao * pi / 180)
 z3d(1) = linz

 linx = x3d(1)
 liny = y3d(1)
 linz = z3d(1)
 x3d(1) = linx
 y3d(1) = liny * Cos(qingjiao * pi / 180) - linz * Sin(qingjiao * pi / 180)
 z3d(1) = liny * Sin(qingjiao * pi / 180) + linz * Cos(qingjiao * pi / 180)

 
     
     
 x2d(1) = zbchangex(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
 y2d(1) = zbchangey(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
 For i = 0 To 50 Step 10
    Form1.Picture1.Circle (x2d(1), y2d(1)), i, vbGreen
 Next i
 Form1.Picture1.Line (0, 0)-(x2d(1), y2d(1)), vbGreen

 '画地固坐标系Z
 linx = 0
 liny = 0
 linz = 8000
 x3d(1) = linx * Cos(zzhuanjiao * pi / 180) - liny * Sin(zzhuanjiao * pi / 180)
 y3d(1) = linx * Sin(zzhuanjiao * pi / 180) + liny * Cos(zzhuanjiao * pi / 180)
 z3d(1) = linz

 linx = x3d(1)
 liny = y3d(1)
 linz = z3d(1)
 x3d(1) = linx
 y3d(1) = liny * Cos(qingjiao * pi / 180) - linz * Sin(qingjiao * pi / 180)
 z3d(1) = liny * Sin(qingjiao * pi / 180) + linz * Cos(qingjiao * pi / 180)

 
     
     
 x2d(1) = zbchangex(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
 y2d(1) = zbchangey(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
 For i = 0 To 50 Step 10
    Form1.Picture1.Circle (x2d(1), y2d(1)), i, vbGreen
 Next i
 Form1.Picture1.Line (0, 0)-(x2d(1), y2d(1)), vbGreen

 '画纬度
 For i = 0 To 180 Step 180 / 10
     x3d(1) = radiu * Sin(i * pi / 180) * Cos(0)
     y3d(1) = radiu * Sin(i * pi / 180) * Sin(0)
     z3d(1) = radiu * Cos(i * pi / 180)
     '沿Z轴旋转
     linx = x3d(1)
     liny = y3d(1)
     linz = z3d(1)
     x3d(1) = linx * Cos(zzhuanjiao * pi / 180) - liny * Sin(zzhuanjiao * pi / 180)
     y3d(1) = linx * Sin(zzhuanjiao * pi / 180) + liny * Cos(zzhuanjiao * pi / 180)
     z3d(1) = linz
     '沿X轴旋转
     linx = x3d(1)
     liny = y3d(1)
     linz = z3d(1)
     x3d(1) = linx
     y3d(1) = liny * Cos(qingjiao * pi / 180) - linz * Sin(qingjiao * pi / 180)
     z3d(1) = liny * Sin(qingjiao * pi / 180) + linz * Cos(qingjiao * pi / 180)

     
     
     x2d(1) = zbchangex(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
     y2d(1) = zbchangey(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
     For j = 0 To 360 Step 2
         x3d(2) = radiu * Sin(i * pi / 180) * Cos(j * pi / 180)
         y3d(2) = radiu * Sin(i * pi / 180) * Sin(j * pi / 180)
         z3d(2) = radiu * Cos(i * pi / 180)
         
         linx = x3d(2)
         liny = y3d(2)
         linz = z3d(2)
         x3d(2) = linx * Cos(zzhuanjiao * pi / 180) - liny * Sin(zzhuanjiao * pi / 180)
         y3d(2) = linx * Sin(zzhuanjiao * pi / 180) + liny * Cos(zzhuanjiao * pi / 180)
         z3d(2) = linz
         
         linx = x3d(2)
         liny = y3d(2)
         linz = z3d(2)
         x3d(2) = linx
         y3d(2) = liny * Cos(qingjiao * pi / 180) - linz * Sin(qingjiao * pi / 180)
         z3d(2) = liny * Sin(qingjiao * pi / 180) + linz * Cos(qingjiao * pi / 180)
         x2d(2) = zbchangex(x3d(2), y3d(2), z3d(2), sita, fai, rou, d)
         y2d(2) = zbchangey(x3d(2), y3d(2), z3d(2), sita, fai, rou, d)
         'If ((x3d(2) - eyex) ^ 2 + (y3d(2) - eyey) ^ 2 + (z3d(2) - eyez) ^ 2) < ((eyex ^ 2 + eyey ^ 2 + eyez ^ 2) + radiu ^ 2) Then
           Form1.Picture1.Line (x2d(1), y2d(1))-(x2d(2), y2d(2)), QBColor(5)
         'End If
            x2d(1) = x2d(2)
            y2d(1) = y2d(2)
            If maxx < x2d(1) Then maxx = x2d(1)
            If maxy < y2d(1) Then maxy = y2d(1)
            If minx > x2d(1) Then minx = x2d(1)
            If miny > y2d(1) Then miny = y2d(1)

            'ymin(x2d(2)) = y2d(2)
         'ElseIf y2d(2) >= ymax(x2d(2)) Then
           ' Form1.Picture1.Line (x2d(1), y2d(1))-(x2d(2), y2d(2)), vbRed
           ' x2d(1) = x2d(2)
           ' y2d(1) = y2d(2)
            'ymax(x2d(2)) = y2d(2)
         'ElseIf y2d(2) > ymin(x2d(2)) And y2d(2) < ymax(x2d(2)) Then
            'x2d(1) = x2d(2)
            'y2d(1) = y2d(2)
        ' End If
     Next j
 Next i
 
 '画经度
 For j = 0 To (360 - 15) Step 30
     x3d(1) = radiu * Sin(0) * Cos(j * pi / 180)
     y3d(1) = radiu * Sin(0) * Sin(j * pi / 180)
     z3d(1) = radiu * Cos(0)
     
     linx = x3d(1)
     liny = y3d(1)
     linz = z3d(1)
     x3d(1) = linx * Cos(zzhuanjiao * pi / 180) - liny * Sin(zzhuanjiao * pi / 180)
     y3d(1) = linx * Sin(zzhuanjiao * pi / 180) + liny * Cos(zzhuanjiao * pi / 180)
     z3d(1) = linz
     
     linx = x3d(1)
     liny = y3d(1)
     linz = z3d(1)
     x3d(1) = linx
     y3d(1) = liny * Cos(qingjiao * pi / 180) - linz * Sin(qingjiao * pi / 180)
     z3d(1) = liny * Sin(qingjiao * pi / 180) + linz * Cos(qingjiao * pi / 180)
     x2d(1) = zbchangex(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
     y2d(1) = zbchangey(x3d(1), y3d(1), z3d(1), sita, fai, rou, d)
     For i = 0 To 180 Step 2
         x3d(2) = radiu * Sin(i * pi / 180) * Cos(j * pi / 180)
         y3d(2) = radiu * Sin(i * pi / 180) * Sin(j * pi / 180)
         z3d(2) = radiu * Cos(i * pi / 180)
         linx = x3d(2)
         liny = y3d(2)
         linz = z3d(2)
         x3d(2) = linx * Cos(zzhuanjiao * pi / 180) - liny * Sin(zzhuanjiao * pi / 180)
         y3d(2) = linx * Sin(zzhuanjiao * pi / 180) + liny * Cos(zzhuanjiao * pi / 180)
         z3d(2) = linz
         linx = x3d(2)
         liny = y3d(2)
         linz = z3d(2)
         x3d(2) = linx
         y3d(2) = liny * Cos(qingjiao * pi / 180) - linz * Sin(qingjiao * pi / 180)
         z3d(2) = liny * Sin(qingjiao * pi / 180) + linz * Cos(qingjiao * pi / 180)
         x2d(2) = zbchangex(x3d(2), y3d(2), z3d(2), sita, fai, rou, d)
         y2d(2) = zbchangey(x3d(2), y3d(2), z3d(2), sita, fai, rou, d)
         If j = 0 Or j = 180 Then
             Form1.Picture1.Line (x2d(1), y2d(1))-(x2d(2), y2d(2)), vbRed
         Else
             Form1.Picture1.Line (x2d(1), y2d(1))-(x2d(2), y2d(2)), QBColor(3)
         End If
         x2d(1) = x2d(2)
         y2d(1) = y2d(2)
         If maxx < x2d(1) Then maxx = x2d(1)
         If maxy < y2d(1) Then maxy = y2d(1)
         If minx > x2d(1) Then minx = x2d(1)
         If miny > y2d(1) Then miny = y2d(1)
         
     Next i
 Next j
 
 
 '画轨道及卫星

 For i = 60 To 360 Step 60
    x3d(1) = orbit * Sin(0) * Cos(i * pi / 180)
    y3d(1) = orbit * Sin(0) * Sin(i * pi / 180)
    z3d(1) = orbit * Cos(0)
    
    linx = x3d(1)
    liny = y3d(1)
    linz = z3d(1)
    x3d(1) = linx * Cos(-i * pi / 180) - liny * Sin(-i * pi / 180)
    y3d(1) = linx * Sin(-i * pi / 180) + liny * Cos(-i * pi / 180)
    z3d(1) = linz
    
    linx = x3d(1)
    liny = y3d(1)
    linz = z3d(1)
    x3d(1) = linx
    y3d(1) = liny * Cos(-35 * pi / 180) - linz * Sin(-35 * pi / 180)
    z3d(1) = liny * Sin(-35 * pi / 180) + linz * Cos(-35 * pi / 180)
    
    linx = x3d(1)
    liny = y3d(1)
    linz = z3d(1)
    x3d(1) = linx * Cos(i * pi / 180) - liny * Sin(i * pi / 180)
    y3d(1) = linx * Sin(i * pi / 180) + liny * Cos(i * pi / 180)