====== 計算プログラム ======

Microsoft ExcelのVBAのプログラム




<code vb prof1d.vba>
Option Explicit
Private Const nix As Integer = 31
Private Const ndx As Integer = 8
Private Const nout As Integer = 20

Private h(nix) As Double
Private z(nix) As Double
Private hs(nix) As Double
Private b(nix) As Double
Private x(nix) As Double
Private us(nix) As Double
Private sib(nix, ndx) As Double
Private sib0(nix, ndx) As Double
Private a(nix) As Double
Private u(nix) As Double
Private si(nix) As Double
Private tsm(nix) As Double
Private dd(ndx) As Double
Private d(ndx) As Double
Private p0(nix, ndx) As Double
Private p(nix, ndx) As Double
Private w0(ndx) As Double
Private uci(ndx) As Double
Private tsci0(ndx) As Double
Private tsi(nix, ndx) As Double
Private wc(nix) As Double
Private d10(nix) As Double
Private d50(nix) As Double
Private d90(nix) As Double
Private dm(nix) As Double
Private dzdt(nix) As Double
Private ct(nix) As Double
Private tscm(nix) As Double
Private qb(nix) As Double
Private qbi(nix, ndx) As Double
Private qs(nix) As Double
Private qsi(nix, ndx) As Double
Private c(nix, ndx) As Double
Private cb(nix, ndx) As Double
Private dm0(nix) As Double
Private z0(nix) As Double

Private sn(nix) As Double
Private dx(nix) As Double


Public Sub pro1d()

    'パラメータ
    Const g As Double = 9.8
    Const s As Double = 1.65
    Const ram As Double = 0.4
    Const delta As Double = 0.1
      
    Const bs0 As Double = 0.143
    Const snu As Double = 0.000001
    Const rs As Double = 1.436
    Const sk As Double = 0.008
    Const als As Double = 0.14
    
    
    Dim i As Integer, j As Integer, k As Integer
    Dim ni As Integer ', dx As Double
    Dim ax As Double
    Dim sechour As Double
    Dim tuk As Double
    Dim etime As Double
    Dim dtmax As Double
    Dim al23 As Double
    Dim nd As Integer
    Dim ttime As Double
    Dim qq As Double
    Dim h00  As Double
    Dim d10_0 As Double, d50_0 As Double, d90_0 As Double
    Dim tt As Double
    Dim ucm As Double
    Dim aa As Double, bb As Double, cc As Double
    Dim fh As Double, fdh As Double, dh As Double
    Dim xi As Double
    Dim tsci As Double
    Dim usci As Double, xxi As Double
    Dim bs As Double
    Dim ome As Double
    Dim bet As Double
    Dim alfx As Double
    Dim amax As Double
    Dim dqb As Double
    Dim ah As Double
    Dim bh1 As Double, bh2 As Double, bh3 As Double
    Dim bh As Double, abh As Double
    Dim tmax As Double
    Dim dt As Double
    Dim dz As Double
    Dim dqbi As Double, dpi As Double
    
    Dim ncount As Long, mcount As Long, mcol As Integer
    
    Dim ws_out As Worksheet
    Dim ws_zd As Worksheet
    Dim ws_h As Worksheet
    Dim ws_z As Worksheet
    Dim ws_dm As Worksheet, ws_d10 As Worksheet, ws_d50 As Worksheet, ws_d90 As Worksheet
    
    Set ws_out = Sheets("計算結果")
    Set ws_zd = Sheets("河床変動量")
    Set ws_h = Sheets("水位")
    Set ws_z = Sheets("河床標高")
    Set ws_dm = Sheets("平均粒径")
    Set ws_d10 = Sheets("d10粒径")
    Set ws_d50 = Sheets("d50粒径")
    Set ws_d90 = Sheets("d90粒径")
    
    
    'Call param
        
    '粒径分布は対数正規分布(8区分)
    For i = 1 To ndx
        dd(1) = 0.1
        dd(2) = 0.2
        dd(3) = 0.5
        dd(4) = 0.7
        dd(5) = 1#
        dd(6) = 2#
        dd(7) = 4#
        dd(8) = 10#
    Next i

    '初期値
    Dim ws1 As Worksheet
    Set ws1 = Sheets("初期値")
    
    qq = ws1.Cells(3, "A").Text     '500#
    h00 = ws1.Cells(3, "B").Text    '3
    'sn = ws1.Cells(3, "C").Text     '0.025
    tuk = ws1.Cells(3, "D").Text    '5#
    etime = ws1.Cells(3, "E").Text  '30#
    
    d10_0 = ws1.Cells(7, "A").Text  '0.2
    d50_0 = ws1.Cells(7, "B").Text  '0.5
    d90_0 = ws1.Cells(7, "C").Text  '1#

    'Call datainput
    
    '河道断面数と初期河床標高の読み込み
    Dim ws0 As Worksheet
    Set ws0 = Sheets("河道データ")
    
    ni = ws0.Cells(1, "B").Text     '31
    'dx = ws0.Cells(1, "C").Text     '200
    
    For i = 1 To ni
       'read(17,*) b(i),z(i)
       b(i) = ws0.Cells(i + 2, "E").Text
       z(i) = ws0.Cells(i + 2, "D").Text
       dx(i) = ws0.Cells(i + 2, "C").Text
       sn(i) = ws0.Cells(i + 2, "F").Text
       z0(i) = z(i)
       d10(i) = d10_0
       d50(i) = d50_0
       d90(i) = d90_0
    Next i

    '計算時間
    ax = 1# / delta
    sechour = 3600#
    tuk = sechour * tuk
    etime = etime * sechour
    dtmax = tuk
    
    al23 = Log10(23#)

    nd = ndx - 1
    For k = 2 To nd
        d(k) = (dd(k) + dd(k - 1)) * 0.5 / 1000#
    Next k
    For i = 1 To ni
        x(i) = (i - 1) * dx(i)
    Next i

    '
    ' ------ initial grain size distribution ------
    '
    'Call initgrain
    
    For i = 1 To ni
        p0(i, 1) = 0#
        p0(i, nd) = 1#
        For k = 2 To nd - 1
            p0(i, k) = ppx(dd(k), d10(i), d50(i), d90(i))
        Next k
    Next i
    For i = 1 To ni
        For k = 2 To nd
            sib(i, k) = p0(i, k) - p0(i, k - 1)
            sib0(i, k) = sib(i, k)
        Next k
        For k = 1 To nd
            p(i, k) = p0(i, k)
        Next k
    Next i

    '
    ' ------ fall velocity -----
    '
    'Call fallvel
    
    For k = 2 To nd
        If d(k) > 0.001 Then
            w0(k) = 32.8 * Sqr(d(k) * 100#) * 0.01
        Else
            w0(k) = Sqr(2# / 3# * s * g * d(k) + _
                    (6# * snu / d(k)) ^ 2) - 6# * snu / d(k)
        End If
    Next k

    '
    ' ------ critical shear velocity ------
    '
    'Call critcalshear
    
    For k = 2 To nd
        Call usc(d(k), uci(k), s, snu, g)
        tsci0(k) = uci(k) ^ 2 / (s * g * d(k))
    Next k

    '
    '------ initial mean diameter ------
    '
    'Call meandia
    
    For i = 1 To ni
        dm0(i) = 0#
        For k = 1 To nd
            dm0(i) = dm0(i) + d(k) * sib(i, k)
        Next k
    Next i

    ttime = 0#
    tt = tuk

    ' ***********************************
    ' 計算開始
    ' ***********************************

    ncount = 0  '出力用のカウンタ
    mcol = 3    '出力用のカウンタ
    
g100: Debug.Print "tt = "; tt
        
    '
    '------  dm ------
    '
    'Call meandia
    
    For i = 1 To ni
        dm(i) = 0#
        For k = 1 To nd
            dm(i) = dm(i) + d(k) * sib(i, k)
        Next k
    Next i

    '
    ' ------  tau*cm ------
    '
    'Call taucm
    For i = 1 To ni
        Call usc(dm(i), ucm, s, snu, g)
        tscm(i) = ucm ^ 2 / (s * g * dm(i))
    Next i

    '
    '------ h ------
    '
    'Call height
    
    h(1) = h00
    hs(1) = h(1) - z(1)
    For i = 2 To ni
        cc = -z(i - 1) + z(i) - hs(i - 1) - _
             qq ^ 2 / (2# * g * b(i - 1) ^ 2 * hs(i - 1) ^ 2) _
             - dx(i - 1) * sn(i - 1) ^ 2 * qq ^ 2 / _
             - (2# * b(i - 1) ^ 2 * hs(i - 1) ^ (10# / 3#))
        hs(i) = hs(i - 1)
g120:   aa = qq ^ 2 / (2# * g * b(i) ^ 2)
        bb = -dx(i) * sn(i) ^ 2 * qq ^ 2 / (2# * b(i) ^ 2)
        fh = hs(i) + aa / hs(i) ^ 2 + bb / hs(i) ^ (10# / 3#) + cc
        fdh = 1# - 2# * aa / hs(i) ^ 3 - 10# / 3# * bb / hs(i) ^ (13# / 3#)
        dh = -fh / fdh
        'Debug.Print cc, hs(i), dh
        If Abs(dh) < 0.0001 Then GoTo g135:
        hs(i) = hs(i) + dh
        GoTo g120:
g135:   h(i) = z(i) + hs(i)
    Next i

    '
    ' ------ a,u,si ------
    '
    'Call ausi
    
    For i = 1 To ni
        a(i) = b(i) * hs(i)
        u(i) = qq / a(i)
        si(i) = (sn(i) * qq / b(i)) ^ 2 / hs(i) ^ 3.3333
        us(i) = Sqr(g * hs(i) * si(i))
        tsm(i) = us(i) ^ 2 / (s * g * dm(i))
    Next i
    
    '
    ' ------ qb (bed load) and qs (suspended load) ------
    '
    'Call bedload
    
    For i = 1 To ni
        qb(i) = 0#
        qs(i) = 0#
        For k = 2 To nd
            tsi(i, k) = us(i) ^ 2 / (s * g * d(k))
            xi = (al23 / Log10(21# * d(k) / dm(i) + 2#)) ^ 2
            tsci = xi * tscm(i)
            If tsi(i, k) <= tsci Then
                qbi(i, k) = 0#
                qsi(i, k) = 0#
            Else
                usci = Sqr(tsci * s * g * d(k))
                qbi(i, k) = sib(i, k) * 17# * Sqr(s * g * d(k) ^ 3) * _
                            tsi(i, k) ^ 1.5 * (1# - tsci / tsi(i, k)) * (1# - usci / us(i))
                If qbi(i, k) < 1E-20 Then qbi(i, k) = 0#
                qb(i) = qb(i) + qbi(i, k)

                xxi = tsci / tsci0(k)
                bs = bs0 * xxi
                Call omega(bs, tsi(i, k), ome)
                qsi(i, k) = sib(i, k) * sk * (als * ome * s * g * d(k) / (rs * us(i)) - w0(k))
                If qsi(i, k) < 1E-20 Then qsi(i, k) = 0#
                qs(i) = qs(i) + qsi(i, k)
            End If
        Next k
    Next i
    
    '
    ' ------ cal. of c ------
    '
    'Call calc
    
    For k = 2 To nd
        bet = 15# * w0(k) / us(ni)
        alfx = alf(bet)
        c(ni, k) = qsi(ni, k) / (w0(k) * alfx)
        cb(ni, k) = c(ni, k) * alfx
        For i = ni - 1 To 1 Step -1
            bet = 15# * w0(k) / us(i)
            alfx = alf(bet)
            c(i, k) = (qq * c(i + 1, k) + dx(i) * b(i) * qsi(i, k)) / _
                      (qq + dx(i) * b(i) * w0(k) * alfx)
            If c(i, k) <= 0.0000000001 Then
                c(i, k) = 0#
            End If
            cb(i, k) = c(i, k) * alfx
        Next i
    Next k

    '
    ' ------ wc, ct------
    '
    'Call wcct
    
    For i = 1 To ni
        ct(i) = 0#
        wc(i) = 0#
        For k = 2 To nd
            ct(i) = ct(i) + c(i, k)
            wc(i) = wc(i) + w0(k) * cb(i, k)
        Next k
    Next i

    '
    ' ------ cal. of dt ------
    '
    'Call caldt
    
    amax = 0#
    For i = 1 To ni - 1
        dqb = (qb(i + 1) * b(i + 1) - qb(i) * b(i)) / (b(i) * dx(i))
        dzdt(i) = (dqb + wc(i) - qs(i)) / (1# - ram)
        For k = 2 To nd
            If (sib(i, k) = 0# Or sib(i + 1, k) = 0#) Then GoTo g130:
            ah = ax / (1# - ram) * qbi(i, k) / sib(i, k)
            bh1 = (qbi(i, k) / sib(i, k) - qbi(i + 1, k) / sib(i + 1, k)) / dx(i)
            bh2 = qbi(i, k) / (b(i) * sib(i, k) * dx(i)) * (b(i) - b(i + 1))
            bh3 = qsi(i, k) / sib(i, k)
            bh = ax / (1# - ram) * (bh1 + bh2 + bh3)
            If dzdt(i) > 0# Then bh = bh + ax * dzdt(i)
            abh = ah / dx(i) + bh * 0.5
            amax = max(amax, abh)
g130:   Next k
    Next i

    If amax > 0# Then
        tmax = 1# / amax
        dt = min(tmax, dtmax)
    Else
        dt = dtmax
    End If
      
    If tt + dt > tuk Then dt = tuk - tt

    '
    '-------- output ------
    '
    'Call dataoutput
    
    mcount = 0  '出力用のカウンタ
    
    If tt >= tuk Then
        Debug.Print "tt=", tt, "dt=", dt, "tuk=", tuk
        tt = 0#
        '
        ' ------  p ------
        '
        For i = 1 To ni
            For k = 2 To nd - 1
                p(i, k) = p(i, k - 1) + sib(i, k)
            Next k
        Next i

        ' ------ d10,d50,d90 ------
        '
        For i = 1 To ni
            For k = 2 To nd
                If (p(i, k) >= 0.1 And p(i, k - 1) <= 0.1) Then
                    d10(i) = 10# ^ ( _
                             Log10(dd(k - 1)) + _
                             (0.1 - p(i, k - 1)) / (p(i, k) - p(i, k - 1)) * _
                             (Log10(dd(k)) - Log10(dd(k - 1))))
                End If
                If (p(i, k) >= 0.5 And p(i, k - 1) <= 0.5) Then
                    d50(i) = 10# ^ ( _
                             Log10(dd(k - 1)) + _
                             (0.5 - p(i, k - 1)) / (p(i, k) - p(i, k - 1)) * _
                             (Log10(dd(k)) - Log10(dd(k - 1))))
                End If
                If (p(i, k) >= 0.9 And p(i, k - 1) <= 0.9) Then
                    d90(i) = 10# ^ ( _
                             Log10(dd(k - 1)) + _
                             (0.9 - p(i, k - 1)) / (p(i, k) - p(i, k - 1)) * _
                             (Log10(dd(k)) - Log10(dd(k - 1))))
                End If
            Next k
        Next i
       
        Debug.Print "time=", ttime / sechour
        
        '計算結果のシート出力
        ncount = ncount + 1
        ws_out.Cells(ncount, "A") = "time:"
        ws_out.Cells(ncount, "B") = Format(ttime / sechour, "######.###0")
        
        '計算結果の出力
        ncount = ncount + 1
        ws_out.Cells(ncount, "C") = "断面i"
        ws_out.Cells(ncount, "D") = "断面距離(m)"
        ws_out.Cells(ncount, "E") = "河床標高(m)"
        ws_out.Cells(ncount, "F") = "河川幅(m)"
        ws_out.Cells(ncount, "G") = "水位(m)"
        ws_out.Cells(ncount, "H") = "平均粒径(mm)"
        ws_out.Cells(ncount, "I") = "d10粒径(mm)"
        ws_out.Cells(ncount, "J") = "d50粒径(mm)"
        ws_out.Cells(ncount, "K") = "d90粒径(mm)"
        ws_out.Cells(ncount, "L") = "初期河床からの変動量(m)"
        
        '河床変動量のみ出力
        mcount = mcount + 1
        ws_zd.Cells(mcount, "A") = "断面i"
        ws_zd.Cells(mcount, "B") = "断面距離(m)"
        ws_zd.Cells(mcount, "C") = "河床標高(m)"
        '河床標高のみ出力
        ws_z.Cells(mcount, "A") = "断面i"
        ws_z.Cells(mcount, "B") = "断面距離(m)"
        ws_z.Cells(mcount, "C") = "河床標高(m)"
        '水位のみ出力
        ws_h.Cells(mcount, "A") = "断面i"
        ws_h.Cells(mcount, "B") = "断面距離(m)"
        ws_h.Cells(mcount, "C") = "河床標高(m)"
        '平均粒径のみ出力
        ws_dm.Cells(mcount, "A") = "断面i"
        ws_dm.Cells(mcount, "B") = "断面距離(m)"
        'd10粒径のみ出力
        ws_d10.Cells(mcount, "A") = "断面i"
        ws_d10.Cells(mcount, "B") = "断面距離(m)"
        'd50粒径のみ出力
        ws_d50.Cells(mcount, "A") = "断面i"
        ws_d50.Cells(mcount, "B") = "断面距離(m)"
        'd90粒径のみ出力
        ws_d90.Cells(mcount, "A") = "断面i"
        ws_d90.Cells(mcount, "B") = "断面距離(m)"
        mcol = mcol + 1
        
        For j = 1 To ni
            ncount = ncount + 1
            mcount = mcount + 1
            ws_out.Cells(ncount, "C") = j
            ws_out.Cells(ncount, "D") = Format(x(j), "####.##0")
            ws_out.Cells(ncount, "E") = Format(z(j), "####.##0")
            ws_out.Cells(ncount, "F") = Format(b(j), "####.##0")
            ws_out.Cells(ncount, "G") = Format(h(j), "####.##0")
            ws_out.Cells(ncount, "H") = Format(dm(j) * 1000#, "###.##0")
            ws_out.Cells(ncount, "I") = Format(d10(j), "####.##0")
            ws_out.Cells(ncount, "J") = Format(d50(j), "####.##0")
            ws_out.Cells(ncount, "K") = Format(d90(j), "####.##0")
            ws_out.Cells(ncount, "L") = Format(z(j) - z0(j), "####.##0")
            
            '河床変動量のみ出力
            ws_zd.Cells(mcount, "A") = j
            ws_zd.Cells(mcount, "B") = Format(x(j), "####.##0")
            ws_zd.Cells(mcount, "C") = Format(z0(j), "####.##0")
            ws_zd.Cells(1, mcol) = "time:" & Format(ttime / sechour, "######.##0")
            ws_zd.Cells(mcount, mcol) = Format(z(j) - z0(j), "####.##0")
            
            '河床標高のみ出力
            ws_z.Cells(mcount, "A") = j
            ws_z.Cells(mcount, "B") = Format(x(j), "####.##0")
            ws_z.Cells(mcount, "C") = Format(z0(j), "####.##0")
            ws_z.Cells(1, mcol) = "time:" & Format(ttime / sechour, "######.###0")
            ws_z.Cells(mcount, mcol) = Format(z(j), "####.##0")
            
            '水位のみ出力
            ws_h.Cells(mcount, "A") = j
            ws_h.Cells(mcount, "B") = Format(x(j), "####.##0")
            ws_h.Cells(mcount, "C") = Format(z0(j), "####.##0")
            ws_h.Cells(1, mcol) = "time:" & Format(ttime / sechour, "######.###0")
            ws_h.Cells(mcount, mcol) = Format(h(j), "####.##0")
        
            '平均粒径のみ出力
            ws_dm.Cells(mcount, "A") = j
            ws_dm.Cells(mcount, "B") = Format(x(j), "####.##0")
            ws_dm.Cells(1, mcol - 1) = "time:" & Format(ttime / sechour, "######.###0")
            ws_dm.Cells(mcount, mcol - 1) = Format(dm(j) * 1000#, "####.##0")
            'd10粒径のみ出力
            ws_d10.Cells(mcount, "A") = j
            ws_d10.Cells(mcount, "B") = Format(x(j), "####.##0")
            ws_d10.Cells(1, mcol - 1) = "time:" & Format(ttime / sechour, "######.###0")
            ws_d10.Cells(mcount, mcol - 1) = Format(d10(j), "####.##0")
            'd50粒径のみ出力
            ws_d50.Cells(mcount, "A") = j
            ws_d50.Cells(mcount, "B") = Format(x(j), "####.##0")
            ws_d50.Cells(1, mcol - 1) = "time:" & Format(ttime / sechour, "######.###0")
            ws_d50.Cells(mcount, mcol - 1) = Format(d50(j), "####.##0")
            'd90粒径のみ出力
            ws_d90.Cells(mcount, "A") = j
            ws_d90.Cells(mcount, "B") = Format(x(j), "####.##0")
            ws_d90.Cells(1, mcol - 1) = "time:" & Format(ttime / sechour, "######.###0")
            ws_d90.Cells(mcount, mcol - 1) = Format(d90(j), "####.##0")
        
        Next j
    
    End If

    '
    ' ------ cal. of z, sib ------
    '
    'Call calzsib
    
    For i = 1 To ni - 1
        dz = dt * dzdt(i)
        z(i) = z(i) + dz
        For k = 2 To nd
            If dz > 0# Then
                xi = sib(i, k)
            Else
                xi = sib0(i, k)
            End If
            dqbi = (qbi(i + 1, k) * b(i + 1) - qbi(i, k) * b(i)) _
                   / (b(i) * dx(i))
            dpi = ax * (dt / (1# - ram) * _
                  (dqbi + w0(k) * cb(i, k) - qsi(i, k)) - xi * dz)
            sib(i, k) = sib(i, k) + dpi
            If sib(i, k) < 0.0000000001 Then sib(i, k) = 0#
        Next k
    Next i


    ttime = ttime + dt
    tt = tt + dt
    If ttime > etime Then GoTo g999:


    GoTo g100:

g999: Close

    ' ***********************************
    ' 計算終了
    ' ***********************************

End Sub

'
' ***********************************
'
Private Sub usc(d As Double, uc As Double, s As Double, snu As Double, g As Double)
    Dim rst As Double
    
    rst = Sqr(s * g) * d ^ (3# / 2#) / snu
    If rst <= 2.14 Then uc = 0.14 * s * g * d
    If rst > 2.14 Then uc = (0.1235 * s * g) ^ (25# / 32#) *_
                             snu ^ (7# / 16#) * d ^ (11# / 32#)
    If rst > 54.2 Then uc = 0.034 * s * g * d
    If rst > 162.7 Then uc = (0.01505 * s * g) ^ (25# / 22#) * _
                              snu ^ (-3# / 11#) * d ^ (31# / 22#)
    If rst > 671# Then uc = 0.05 * s * g * d
    
    uc = Sqr(uc)

End Sub

'
' ***********************************
'
Private Sub omega(bs As Double, tsi As Double, ome As Double)
    Dim pai As Double
    Dim a1 As Double, a2 As Double, a3 As Double
    Dim ad As Double
    Dim t As Double, zx As Double, px As Double
    Dim er1 As Double, er As Double
    Dim xxx As Double
    
    a1 = 0.4361836
    a2 = -0.1201676
    a3 = 0.937298
    pai = 3.141592
    
    If tsi <= 0.0000001 Then
        ome = 0#
        Exit Sub
    End If
    
    ad = bs / tsi - 2#
    If ad >= 0# Then
        xxx = ad * Sqr(2#)
    End If
    If ad < 0# Then
        xxx = -ad * Sqr(2#)
    End If
    
    t = 1# / (1# + 0.33627 * xxx)
    zx = 1# / Sqr(2# * pai) * Exp(-xxx ^ 2# / 2#)
    px = 1# - zx * (a1 * t + a2 * t ^ 2 + a3 * t ^ 3)
    
    er1 = 2# - 2# * px
    If ad >= 0# Then
        er = er1 / 2#
    End If
    If ad < 0# Then
        er = (2# - er1) / 2#
    End If
    
    If (bs = 0# Or er = 0#) Then
        ome = 0#
    Else
        ome = tsi / bs / (2# * Sqr(pai)) * Exp(-ad ^ 2) / er + tsi * 2# / bs - 1#
    End If
    
    'if(ome.lt.1e-10) ome=0.

End Sub

'
' ***********************************
'
Private Function alf(bet As Double)
    'if bet > 20. then
    '   alf = bet
    'else
       alf = bet / (1# - Exp(-bet))
    'end if

End Function

'
' ***********************************
'
Private Function ppx(dd As Double, d10 As Double, d50 As Double, d90 As Double)
    Dim s1 As Double, s2 As Double
    Dim x As Double
      
    s1 = 1.282 / (Log10(d90) - Log10(d50))
    s2 = -1.282 / (Log10(d10) - Log10(d50))

    If dd >= d50 Then
        x = (Log10(dd) - Log10(d50)) * s1
    Else
        x = (Log10(dd) - Log10(d50)) * s2
    End If
    
    ppx = pc(x)

End Function

'
' ***********************************
'
Private Function tfunc(x As Double)
    If x >= 0# Then
        tfunc = 1# / (1# + 0.33267 * x)
    Else
        tfunc = 1# / (1# - 0.33267 * x)
    End If

End Function

'
' ***********************************
'
Private Function pc(x As Double)
    Dim pai As Double
    Dim a1 As Double, a2 As Double, a3 As Double
    Dim tx As Double
    
    pai = 3.141592
    
    a1 = 0.4361836
    a2 = -0.1201676
    a3 = 0.937298
    tx = tfunc(x)
      
    If x >= 0# Then
       pc = 1# - 1# / Sqr(2# * pai) * _
       Exp(-x ^ 2 / 2#) * (a1 * tx + a2 * tx ^ 2 + a3 * tx ^ 3)
    Else
       pc = 1# / Sqr(2# * pai) * _
       Exp(-x ^ 2 / 2#) * (a1 * tx + a2 * tx ^ 2 + a3 * tx ^ 3)
    End If
    
End Function

'
' ***********************************
'
Private Function max(a As Double, b As Double)
    max = a
    If max < b Then max = b
End Function

'
' ***********************************
'
Private Function min(a As Double, b As Double)
    min = a
    If min > b Then min = b
End Function

Private Function Log10(ByVal Val As Double)
    Log10 = Log(Val) / Log(10)
End Function

</code>