VERSION 5.00
Object = "{F9043C88-F6F2-101A-A3C9-08002B2F49FB}#1.2#0"; "Comdlg32.ocx"
Begin VB.Form frmvkform 
   AutoRedraw      =   -1  'True
   Caption         =   "產(chǎn)生垂向滲透率值"
   ClientHeight    =   6780
   ClientLeft      =   45
   ClientTop       =   270
   ClientWidth     =   6615
   LinkTopic       =   "Form1"
   ScaleHeight     =   6780
   ScaleWidth      =   6615
   StartUpPosition =   3  '窗口缺省
   Begin MSComDlg.CommonDialog cmndia 
      Left            =   3120
      Top             =   6120
      _ExtentX        =   688
      _ExtentY        =   688
      _Version        =   393216
   End
   Begin VB.CommandButton cmdsave 
      Caption         =   "保存(&S)"
      Height          =   372
      Left            =   1320
      TabIndex        =   10
      Top             =   6120
      Width           =   1332
   End
   Begin VB.CommandButton cmdexit 
      Caption         =   "退出(&E)"
      Height          =   372
      Left            =   4200
      TabIndex        =   9
      Top             =   6120
      Width           =   1332
   End
   Begin VB.PictureBox picout 
      AutoSize        =   -1  'True
      Height          =   2415
      Left            =   240
      ScaleHeight     =   2355
      ScaleWidth      =   5835
      TabIndex        =   7
      Top             =   3360
      Width           =   5892
   End
   Begin VB.CommandButton cmdcal 
      Caption         =   "計(jì) 算 各 層 滲 透 率 值"
      Height          =   372
      Left            =   1440
      TabIndex        =   6
      Top             =   2640
      Width           =   3372
   End
   Begin VB.TextBox txtnlayer 
      Height          =   372
      Left            =   3360
      TabIndex        =   5
      Top             =   1920
      Width           =   2172
   End
   Begin VB.TextBox txtvk 
      Height          =   372
      Left            =   3360
      TabIndex        =   3
      Top             =   1080
      Width           =   2172
   End
   Begin VB.TextBox txtavk 
      Height          =   372
      Left            =   3360
      TabIndex        =   1
      Top             =   360
      Width           =   2172
   End
   Begin VB.Label Label4 
      Caption         =   "點(diǎn)擊“計(jì)算各層滲透率值”按鈕，在以下框中查看結(jié)果"
      ForeColor       =   &H00C00000&
      Height          =   252
      Left            =   240
      TabIndex        =   8
      Top             =   3960
      Width           =   5052
   End
   Begin VB.Label label3 
      Caption         =   "請(qǐng)輸入垂向上的層數(shù)"
      ForeColor       =   &H000040C0&
      Height          =   252
      Left            =   840
      TabIndex        =   4
      Top             =   2040
      Width           =   1812
   End
   Begin VB.Label Label2 
      Caption         =   "請(qǐng)輸入滲透率變異系數(shù)"
      ForeColor       =   &H000040C0&
      Height          =   252
      Left            =   720
      TabIndex        =   2
      Top             =   1200
      Width           =   1812
   End
   Begin VB.Label Label1 
      Caption         =   "請(qǐng)輸入滲透率的平均值 (mdc)"
      ForeColor       =   &H000040C0&
      Height          =   252
      Left            =   480
      TabIndex        =   0
      Top             =   480
      Width           =   2412
   End
End
Attribute VB_Name = "frmvkform"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Dim avk As Single, vk As Single, nlayer As Integer
Dim perm() As Single, frach() As Single, phi() As Single
Dim hh() As Single
Const pi = 3.1415926

Private Sub cmdcal_Click()
    Dim i As Integer
    txtnlayer.Enabled = False
    cmdsave.Enabled = True
    avk = txtavk.Text
    vk = txtvk.Text
'    nlayer(層數(shù)) = txtnlayer.Text
    nlayer = 5
    Call calperm
    picout.Cls
    For i = 1 To nlayer
        picout.Print "第" & i & "層的滲透率為：" & Format(perm(i), "######.##") _
        & "   孔隙度為：" & Format(phi(i), "##.###")
    Next i
End Sub

Private Sub cmdexit_Click()
    Unload Me
End Sub

Private Sub calperm()
    ReDim frach(1 To nlayer)
    ReDim hh(1 To nlayer), perm(1 To nlayer), phi(1 To nlayer)
    Dim fy1 As Single, fy2 As Single
    Dim i As Integer, yy As Single, th As Single, ah As Single
    Dim slope As Single, cut As Single
    Dim k50 As Single, k841 As Single
    Dim aa As Double, bb As Double
    hh(1) = 3: hh(2) = 3: hh(3) = 3: hh(4) = 3: hh(5) = 3
    k50 = avk
    k841 = k50 - vk * k50
    fy1 = bidivide(0.5)
    fy2 = bidivide(0.841)
    slope = (log10(k50) - log10(k841)) / (fy1 - fy2)
    cut = log10(k50) - slope * fy1
    th = 0#
    For i = 1 To nlayer
      th = th + hh(i)
    Next i
    ah = 0#
    For i = 1 To nlayer - 1
        ah = ah + hh(i)
        frach(i) = ah / th
        If frach(i) < 0.5 Then frach(i) = 1.3 * frach(i)
        If frach(i) > 0.5 Then frach(i) = frach(i) / 1.3
        
        yy = bidivide(frach(i))
        perm(i) = slope * yy + cut
        perm(i) = pow(10, perm(i))
    Next i
    yy = 3.71
    perm(nlayer) = slope * yy + cut
    perm(nlayer) = pow(10, perm(nlayer))
    aa = 0.0551
    bb = 0.134587
    For i = 1 To nlayer
        phi(i) = aa * log10(perm(i)) + bb
    Next i
End Sub


Private Function bidivide(aa As Single)
    Dim xx As Single, err As Single, yy As Single
    Dim min As Single, max As Single, nn As Integer
    err = 0.00001
    min = -4.9
    max = 4.9
    Do While (1)
        xx = (max + min) / 2
        yy = erf(xx)
        If (Abs(yy - aa) < err) Then
            bidivide = xx
            Exit Function
        ElseIf (yy < aa) Then
            min = xx
        ElseIf (yy > aa) Then
            max = xx
        End If
    Loop
End Function

Private Function erf(z) As Single
    Dim a1 As Single, a2 As Single, a3 As Double, sum0 As Double
    Dim i As Integer, kk As Integer, mm As Integer, mm1 As Integer, mm2 As Integer
    Dim nn1 As Double
    Dim x1 As Single, x2 As Single, err As Single, zz1 As Single
    If (z < 0) Then
        zz1 = -z
    ElseIf (z > 0) Then
        zz1 = z
    Else
        erf = 0.5
        Exit Function
    End If
    err = 0.000001
    a1 = 1 / Sqr(2 * pi)
    
    a2 = 1 / Exp(zz1 * zz1 / 2)
    
    kk = 1
    '================================
    x1 = 0.5 + a1 * a2 * zz1
    
    sum0 = zz1
    
    Do While (1)
        a3 = pow(zz1, 2 * kk + 1)
        
        mm = 2 * kk + 1
        
        nn1 = 1
        
        For i = 1 To mm Step 2
        
            nn1 = nn1 * i
            
        Next i
        
        sum0 = sum0 + a3 / nn1
        
        x2 = 0.5 + a1 * a2 * sum0
        
        If (Abs(x1 - x2) < err) Then
        
            If (z < 0) Then
               erf = 1 - x2
            Else
                erf = x2
            End If
            '======
            Exit Function
        Else
            x1 = x2
            kk = kk + 1
        End If
    Loop
End Function

Private Function pow(xx, yy) As Double
    Dim zz As Single
    zz = yy * Log(xx)
    pow = Exp(zz)
End Function

Private Function log10(xx As Single)
    log10 = Log(xx) / Log(10)
End Function


Private Sub cmdsave_Click()
    Dim fileout As String, i As Integer
    cmndia.ShowSave
    cmndia.Filter = "數(shù)據(jù)文件(*.dat)|*.dat|文本文件(*.txt)|*.txt"
    cmndia.FilterIndex = 1
    fileout = cmndia.FileName
    On Error GoTo errhandler
    Open fileout For Output As #1
    For i = 1 To nlayer
'        Write #1, Format(perm(i), "######.##"), Format(phi(i), "##.###")
    Print #1, Format(perm(i), "0.#"), Format(phi(i), "######.###"), Format(hh(i), "######.##")
'    Write #1, perm(i), hh(i)
    Next i
    Close #1
    Exit Sub
errhandler:
    MsgBox "打開(kāi)文件時(shí)出錯(cuò)！", vbExclamation, "警告！"
    Exit Sub
End Sub

Private Sub Form_Load()
    cmdsave.Enabled = False
    txtnlayer.Visible = False
End Sub