曲面_立体图m2.bas

<VB数理统计实用算法>书中的算法源程序

Attribute VB_Name = "modDraw"
'立体图方法
Public X(101) As Double, Y(101) As Double
Public H(101, 101) As Double, H1(101, 101) As Double
Public PX(6) As Double, PY(6) As Double, QX(6) As Double, QY(6) As Double
Public DPX(6) As Double, DPY(6) As Double, DXY(6) As Double, DS(6) As Double
Public DX As Double, DY As Double, WX As Double, WY As Double
Public TXX As Double, TXY As Double, TYX As Double, TYY As Double, TYZ As Double
Public DOT As Double, DOTC As Double
Public KPOINT As Integer, JS As Integer, JE As Integer, DJ As Integer, KS As Integer
Public IE As Integer, DI As Integer, NX As Integer, NY As Integer
Public J As Integer, I As Integer, JJ As Integer, II As Integer, K As Integer
Public CST As Double, CEN As Double, CIN As Double, CS0 As Double, CE0 As Double
Public XIC As Double, YJC As Double, XP As Double, YP As Double
Public COLB As Integer, COLM As Integer, COLG As Integer, COLC As Integer
Public Z0 As Double, VE As Double, HH1 As Double, HH2 As Double
Public NFIL As String, strData As String
Public Const PI = 3.14159 / 180

'读数据
Public Sub FDATA1()
    Dim intI As Integer, intJ As Integer
    Dim HMAX As Double, HMIN As Double
    strFileName = NFIL                  '文件名
    intFileNumber = FreeFile            '取得空闲的文件号码
    Open strFileName For Input As intFileNumber
    Input #intFileNumber, strData       '读列数
    intCol = Val(strData)               '取得列数
    If intCol >= 2 Then
        For intI = 2 To intCol          '空转，读*****
            Input #intFileNumber, strData
        Next intI
    End If
    Input #intFileNumber, strData       '读行数
    intRow = Val(strData)               '取得行数
    If intCol >= 2 Then
        For intI = 2 To intCol          '空转，读*****
            Input #intFileNumber, strData
        Next intI
    End If
    Input #intFileNumber, strData       '读总行数
    intRowAll = Val(strData)            '取得总行数
    If intCol >= 2 Then
        For intI = 2 To intCol          '空转，读*****
            Input #intFileNumber, strData
        Next intI
    End If
    blnTitle = False: blnRowLabel = False: blnColLabel = False
'优先考虑图题
    If intRowAll > intRow + 3 Then blnTitle = True      '有图题
'其次考虑行标
    If intRowAll > 2 * intRow + 3 Then
        blnRowLabel = True                              '有行标
        ReDim strRowLabel(1 To intRow)                  '重新定义行标数组
    End If
'最后考虑列标
    If intRowAll > 2 * intRow + 4 Then
        blnColLabel = True                              '有列标
        ReDim strColLabel(1 To intCol)                  '重新定义列标数组
    End If
    If blnTitle Then
        Input #intFileNumber, strData   '读图形标题
        strLabelName = strData          '保存图题
        If intCol >= 2 Then
            For intI = 2 To intCol      '空转，读*****号
                Input #intFileNumber, strData
            Next intI
        End If
    End If
    If blnRowLabel Then
        For intI = 1 To intRow
            Input #intFileNumber, strData               '读行标题
            If intCol >= 2 Then
                For intJ = 2 To intCol                  '空转，读*****号
                    Input #intFileNumber, strData
                Next intJ
            End If
        Next intI
    End If
    If blnColLabel Then
        For intI = 1 To intCol                          '读列标题
            Input #intFileNumber, strData
        Next intI
    End If
    For intI = 1 To intRow
        For intJ = 1 To intCol
            Input #intFileNumber, strData               '读图形数据
            H(intJ, intI) = Val(strData)
            H1(intJ, intI) = H(intJ, intI)
        Next intJ
    Next intI
    Close
    frmPicture.lblTitle = strLabelName
    NX = intCol: NY = intRow
'求极大值和极小值
    HMAX = H(1, 1): HMIN = H(1, 1)
    For I = 1 To NX
        For J = 1 To NY
            If HMAX < H(I, J) Then HMAX = H(I, J)
            If HMIN > H(I, J) Then HMIN = H(I, J)
        Next J
    Next I
    frmPicture.TextST.Text = HMIN                       '起始线高度
    frmPicture.TextEN.Text = HMAX                       '终止线高度
    frmPicture.TextIV.Text = (HMAX - HMIN) / 10         '线与线间隔
End Sub

'投影
Public Sub PROJCT()
'(X,Y,Z)投影到曲面(PX,PY) 和基面 (QX,QY)
    QX(K) = TXX * X(II) + TXY * Y(JJ)
    QY(K) = TYX * X(II) + TYY * Y(JJ)
    If H(II, JJ) > Z0 Then HH = (H(II, JJ) - Z0) * VE Else HH = 0
    PX(K) = QX(K)
    PY(K) = QY(K) + TYZ * HH
End Sub

'对每个网格绘等值线
Public Sub DRWC()
    Dim XP As Double, YP As Double
    Dim LXP As Double, LYP As Double
    XP = XIC
    YP = YJC
    If HH1 > Z0 Then HH2 = (HH1 - Z0) * VE Else HH2 = 0
    LXP = TXX * XP + TXY * YP
    LYP = TYX * XP + TYY * YP + TYZ * HH2
    If KPOINT = 1 Then GoTo ENDPOINT
'开始点
    frmPicture.pic.PSet (LXP, -LYP)
    KPOINT = 1
    Exit Sub
'结束点
ENDPOINT:
    frmPicture.pic.Line -(LXP, -LYP), QBColor(COLC)
    KPOINT = 0
End Sub

'对每个网格寻找等值点
Public Sub LCROSS()
    Dim LJ As Integer, LI As Integer, LJ1 As Integer, LI1 As Integer
    Dim LII As Integer, LJJ As Integer, ICS As Integer, ICE As Integer
    Dim IC As Integer
    Dim FBOT As Double, FTOP As Double, FLBP As Double, FLTP As Double
    Dim CON As Double, FC00 As Double, FC01 As Double, FC10 As Double
    Dim FC11 As Double
    LJ = J
    LJ1 = J + DJ
    LI = II - DI
    LI1 = II
'1。数据的最大值和最小值
    FBOT = H(LI, LJ): FTOP = H(LI, LJ)
    For LII = LI To LI1 Step DI
        For LJJ = LJ To LJ1 Step DJ
            If H(LII, LJJ) < FBOT Then FBOT = H(LII, LJJ)
            If H(LII, LJJ) > FTOP Then FTOP = H(LII, LJJ)
        Next
    Next
'2。通过网格（I，J）的等值线
    If FBOT < CST Then FLBT = CST Else FLBT = FBOT
    ICS = Int(FLBT / CIN + CS0)
    If FTOP > CEN Then FLTP = CEN Else FLTP = FTOP
    ICE = Int(FLTP / CIN + CE0)
    If ICS > ICE Then Exit Sub
'3。对网格（I，J）寻找等值点
'   从ICS到ICE的等值线
    For IC = ICS To ICE
        CON = CST + CIN * IC
        KPOINT = 0
        FC00 = H(LI, LJ) - CON: FC10 = H(LI1, LJ) - CON
        FC01 = H(LI, LJ1) - CON: FC11 = H(LI1, LJ1) - CON
'(3-1)。 检查H(LI,LJ) 和 H(LI+1,LJ)
        If FC10 * FC00 > 0 Then GoTo C1
        If H(LI1, LJ) = H(LI, LJ) Then GoTo C1
        XIC = -DI * DX * FC00 / (H(LI1, LJ) - H(LI, LJ)) + X(LI)
        YJC = Y(LJ)
        HH1 = DI * (XIC - X(LI)) * (H(LI1, LJ) - H(LI, LJ)) / DX + H(LI, LJ)
        DRWC
'(3-2)。 检查H(I+1,J) 和 H(I+1,J+1)
C1:
        If FC11 * FC10 > 0 Then GoTo C2
        If H(LI1, LJ1) = H(LI1, LJ) Then GoTo C2
        XIC = X(LI1)
        YJC = -DJ * DY * FC10 / (H(LI1, LJ1) - H(LI1, LJ)) + Y(LJ)
        HH1 = DJ * (YJC - Y(LJ)) * (H(LI1, LJ1) - H(LI1, LJ)) / DY + H(LI1, LJ)
        DRWC
'(3-3)。 检查H(I+1,J+1) 和 H(I,J+1)
C2:
        If FC01 * FC11 > 0 Then GoTo C3
        If H(LI1, LJ1) = H(LI, LJ1) Then GoTo C3
        XIC = -DI * DX * FC01 / (H(LI1, LJ1) - H(LI, LJ1)) + X(LI)
        YJC = Y(LJ1)
        HH1 = DI * (XIC - X(LI)) * (H(LI1, LJ1) - H(LI, LJ1)) / DX + H(LI, LJ1)
        DRWC
'(3-4)。 检查 H(I,J+1) 和 H(I,J)
C3:
        If FC01 * FC00 > 0 Then GoTo C4
        If H(LI, LJ1) = H(LI, LJ) Then GoTo C4
        XIC = X(LI)
        YJC = -DJ * DY * FC00 / (H(LI, LJ1) - H(LI, LJ)) + Y(LJ)
        HH1 = DJ * (YJC - Y(LJ)) * (H(LI, LJ1) - H(LI, LJ)) / DY + H(LI, LJ)
        DRWC
C4:
    Next
End Sub

Public Sub BLOCK1()
    For J = JS To JE - DJ Step DJ
        II = IE
        K = 2: JJ = J: PROJCT
        K = 3: JJ = J + DJ: PROJCT
        For II = KS To IE Step DI
            PX(1) = PX(2): PX(4) = PX(3)
            PY(1) = PY(2): PY(4) = PY(3)
            K = 2: JJ = J: PROJCT
            K = 3: JJ = J + DJ: PROJCT
            If frmPicture.checkXY.Value = 1 Then BPNT2
            If frmPicture.CheckZ.Value = 1 Then LCROSS
        Next II
    Next J
End Sub

Public Sub BLOCK2()
    For J = JS To JE - DJ Step DJ
        II = KS
        K = 2: JJ = J:  PROJCT
        K = 3: JJ = J + DJ:  PROJCT
        For II = KS + DI To IE Step DI
            PX(1) = PX(2): PX(4) = PX(3)
            PY(1) = PY(2): PY(4) = PY(3)
            K = 2: JJ = J:  PROJCT
            K = 3: JJ = J + DJ:  PROJCT
            If frmPicture.CheckZ.Value = 1 Then BPNT2
            If frmPicture.checkXY.Value = 1 Then LCROSS
        Next II
    Next J
End Sub

'画一个块的曲面
Public Sub BPNT1()
    Dim COLP1 As Integer, COLP2 As Integer, COLP As Integer
    If J = JS Then COLP1 = COLW Else COLP1 = COLG
    If II = IE Then COLP2 = COLW Else COLP2 = COLG
    If J = JE - DJ Then COLP3 = COLW Else COLP3 = COLG
    If II = KS + DI Then COLP4 = COLW Else COLP4 = COLG
    frmPicture.pic.Line (PX(1), -PY(1))-(PX(2), -PY(2)), QBColor(COLP1)
    frmPicture.pic.Line -(PX(3), -PY(3)), QBColor(COLP2)
    frmPicture.pic.Line -(PX(4), -PY(4)), QBColor(COLP3)
    frmPicture.pic.Line -(PX(1), -PY(1)), QBColor(COLP4)
End Sub

'画一个块的曲面
Public Sub BPNT2()
    Dim SS As Double, PP As Double, RR As Double
    Dim K1 As Integer, K2 As Integer, K3 As Integer
    Dim COLP1 As Integer, COLP2 As Integer, COLP As Integer
    Dim PXC As Double, PYC As Double
    PX(5) = PX(1): PY(5) = PY(1)
    PX(6) = PX(3): PY(6) = PY(3)
    For K = 1 To 5
        DPX(K) = PX(K + 1) - PX(K): DPY(K) = PY(K + 1) - PY(K)
        DXY(K) = Sqr(DPX(K) ^ 2 + DPY(K) ^ 2)
    Next
    K2 = 4
BP2:
    K1 = K2 - 1: K3 = K2 + 1
    PP = (DXY(K1) + DXY(K2) + DXY(K3)) / 2
    SS = (DPX(K1) * DPY(K2) - DPX(K2) * DPY(K1)) / 2
    RR = Abs(SS) / PP
    If Abs(SS / PP) <= DOT Then GoTo BP1
    frmPicture.pic.FillStyle = 0
    frmPicture.pic.Line (PX(K3), -PY(K3))-(PX(K1), -PY(K1)), QBColor(COLM)
    frmPicture.pic.Line -(PX(K2), -PY(K2)), QBColor(COLM)
    frmPicture.pic.Line -(PX(K3), -PY(K3)), QBColor(COLM)
    PXC = (PX(K1) + PX(K2) + PX(K3)) / 3
    PYC = (PY(K1) + PY(K2) + PY(K3)) / 3
    frmPicture.pic.PSet (PXC, -PYC)
BP1:
    COLP1 = COLG: COLP2 = COLG
    If J = JS And K2 = 2 Then COLP1 = COLW
    If J = JE - DJ And K2 = 4 Then COLP1 = COLW
    If II = KS + DI And K2 = 4 Then COLP2 = COLW
    If II = IE And K2 = 2 Then COLP2 = COLW
    frmPicture.pic.Line (PX(K3), -PY(K3))-(PX(K1), -PY(K1)), QBColor(0)
    frmPicture.pic.Line -(PX(K2), -PY(K2)), QBColor(COLP1)
    frmPicture.pic.Line -(PX(K3), -PY(K3)), QBColor(COLP2)
    If K2 = 4 Then K2 = 2: GoTo BP2
End Sub

Public Sub Surface()
    Dim HMAX As Double, HMIN As Double, D As Double
    Dim XO As Double, YO As Double
    Dim I As Integer, J As Integer
    Dim RAL As Double, RGM As Double, ALPHA As Double, GAMMA As Double
    Dim WX1 As Double, WX2 As Double, WY1 As Double, WY2 As Double
    COLB = 1: COLM = 4: COLG = 1: COLW = 6: COLC = 4
    ALPHA = Val(frmPicture.Text2.Text)          '旋转角
    GAMMA = Val(frmPicture.Text3.Text)          '视角
    Z0 = Val(frmPicture.TextZ0.Text)            '基面高度
    DSCALE = Val(frmPicture.Text4.Text)         '整体伸缩
    VE = Val(frmPicture.TextVE.Text)            '垂直伸缩
    CST = Val(frmPicture.TextST.Text)           '起始等值线
    CEN = Val(frmPicture.TextEN.Text)           '终止等值线
    CIN = Val(frmPicture.TextIV.Text)           '线与线间距
    DX = Val(frmPicture.txtX.Text)      'DX为X方向数据点间隔，即网格的X边长
    DY = Val(frmPicture.txtY.Text)      'DY为Y方向数据点间隔，即网格的Y边长
    For I = 1 To intRow
        For J = 1 To intCol
            H(J, I) = H1(J, I)
        Next J
    Next I
    If frmPicture.checkC.Value Then
'数据列倒转
        For I = 1 To intRow
            For J = 1 To intCol \ 2
                D = H(intCol - J + 1, I)
                H(intCol - J + 1, I) = H(J, I)
                H(J, I) = D
            Next J
        Next I
    End If
    If frmPicture.CheckR.Value Then
'数据行倒转
        For I = 1 To intRow \ 2
            For J = 1 To intCol
                D = H(J, intRow - I + 1)
                H(J, intRow - I + 1) = H(J, I)
                H(J, I) = D
            Next J
        Next I
    End If
    WX = DX * (NX - 1): WY = DY * (NY - 1)
    XO = -WX / 2: YO = -WY / 2
    For I = 0 To NX + 1: X(I) = XO + DX * (I - 1): Next
    For J = 0 To NY + 1: Y(J) = YO + DY * (J - 1): Next
    RAL = ALPHA * PI: RGM = GAMMA * PI
    TXX = Cos(RAL): TXY = -Sin(RAL)
    TYX = Sin(RAL) * Sin(RGM): TYY = Cos(RAL) * Sin(RGM): TYZ = Cos(RGM)
    If Cos(RAL) >= 0 Then JS = NY: JE = 1: DJ = -1 Else JS = 1: JE = NY: DJ = 1
    If Sin(RAL) >= 0 Then KS = NX: IE = 1: DI = -1 Else KS = 1: IE = NX: DI = 1
    If WX > WY Then WW = WX Else WW = WY
    WX1 = -1.6 * WW / DSCALE: WY1 = -1.25 * WW / DSCALE
    WX2 = -WX1: WY2 = 0.75 * WW / DSCALE
    frmPicture.pic.Scale (WX1, WY1)-(WX2, WY2)          '自定义坐标系
    DOT = WW / (200 * DSCALE): DOTC = 5 * DOT
    If CIN = 0 Then
        MsgBox "您必须先输入数据文件，然后才能绘图！", , "立体图"
    End
    End If
    CS0 = -CST / CIN + 0.9999
    CE0 = -CST / CIN + 0.0001
'基面(Z=Z0)的框架
    For K = 1 To 4
        If (K - 1) * (K - 4) = 0 Then II = 1 Else II = NX
        If (K - 1) * (K - 2) = 0 Then JJ = 1 Else JJ = NY
        PROJCT
    Next K
    frmPicture.pic.PSet (QX(4), -QY(4))
'曲面图和等值线图
    frmPicture.pic.Cls
    If frmPicture.Option1.Value = True Then BLOCK1
    If frmPicture.Option2.Value = True Then BLOCK2
'框架
    II = KS: JJ = JE: K = 1:  PROJCT
    II = IE: JJ = JE: K = 2:  PROJCT
    II = IE: JJ = JS: K = 3:  PROJCT
    frmPicture.pic.Line (QX(1), -QY(1))-(QX(2), -QY(2)), QBColor(7)
    frmPicture.pic.Line (QX(2), -QY(2))-(QX(3), -QY(3)), QBColor(7)
    frmPicture.pic.Line (PX(1), -PY(1))-(QX(1), -QY(1)), QBColor(7)
    frmPicture.pic.Line (PX(2), -PY(2))-(QX(2), -QY(2)), QBColor(7)
    frmPicture.pic.Line (PX(3), -PY(3))-(QX(3), -QY(3)), QBColor(7)
End Sub