Imports System.Math
Imports System.IO
Module Module1
 Sub Main()
  Dim auc As Double = 0
  Dim ba As Double = 1
  Dim cl As Double = 0.3
  Dim Cs As Double = 0.8
  Dim d As Double = 0.0003
  Dim dlt As Double = 0.003
  Dim dose As Double = 1
  Dim f(10, 100, 100) As Double
  Dim g As Double = 0
  Dim h As Double = 0.01
  Dim i, j, k As Integer
  Dim ka As Double = 0.03
  Dim kn1(10, 100, 100) As Double
  Dim kn2(10, 100, 100) As Double
  Dim kn3(10, 100, 100) As Double
  Dim kn4(10, 100, 100) As Double
  Dim k12 As Double = 0.0
  Dim k21 As Double = 0.0
  Dim m As Integer = 4
  Dim mz(100) As Double
  Dim n As Integer = 1
  Dim nwe As Double
  Dim nz(100) As Double
  Dim p As Double = 1300
  Dim pc As Double
  Dim pd As Double
  Dim q As Double = 15
  Dim r(100) As Double
  Dim rt(100, 100) As Double
  Dim s(10) As Double
  Dim t(100, 100) As Double
  Dim V As Double = 250
  Dim vd As Double = 100
  Dim vz(100) As Double
  Dim X As Double = 0
  Dim xmax As Double = 5760
  Dim Xso(100, 100) As Double
  Dim Y(10, 100, 100) As Double
  Dim yn(10, 100, 100) As Double
  Dim file As StreamWriter
  file = My.Computer.FileSystem.OpenTextFileWriter(
   "C:\Users\Kevin C. Johnson\multidose.txt", False)
  Dim di As Double = 240
  Dim o As Integer = xmax / di
  Dim dn As Integer = 1
  Dim w As Double = 70
  Dim sd As Double = 2
  If sd = 1 Then n = 1
  Dim dmean As Double = 0.001
  Dim rmean As Double = dmean / 2
  Dim rmin As Double = rmean / sd ^ 3
  Dim rmax As Double = rmean * sd ^ 3
  Dim rinc As Double
  If n > 1 Then rinc = Log(rmax / rmin) / (n - 1)
  Dim freq(100) As Double
  Dim tfreq As Double = 0
  For i = 1 To n
   If n > 1 Then r(i) = Exp(rinc * (i - 1) + Log(rmin)) Else r(i) = rmean
   If n > 1 Then
    freq(i) = 1 / Log(sd) / Sqrt(2 * PI) * Exp(-((Log(r(i)) - Log(rmean)) / Log(sd)) ^ 2 / 2)
   Else
    freq(i) = 1
   End If
   tfreq = tfreq + freq(i)
   vz(i) = 4 / 3 * PI * r(i) ^ 3
   mz(i) = vz(i) * p
  Next
  Dim s2(100) As Double
  Dim mdc As Double = 0
  Dim dy As Double = 0
  For i = 1 To m
   For j = 1 To n
    For k = 1 To o
     Y(i, j, k) = 0
     Xso(j, k) = 0
    Next
   Next
  Next
Start:
  For i = 1 To m
   s(i) = 0
  Next
  For i = 1 To o
   s2(i) = 0
  Next
  If mdc > di - h / 2 Then dn = dn + 1
  If mdc > di - h / 2 Then mdc = 0
  For j = 1 To n
   If mdc = 0 Then Xso(j, dn) = freq(j) / tfreq * dose
   If mdc = 0 Then Y(1, j, dn) = Xso(j, dn)
  Next
  mdc = mdc + h
  For i = 1 To m
   For j = 1 To n
    For k = 1 To o
     s(i) = s(i) + Y(i, j, k)
    Next
   Next
  Next
  For j = 1 To n
   For k = 1 To o
    s2(k) = s2(k) + Y(2, j, k)
    rt(j, k) = (3 * Y(1, j, k) / (4 * PI * p * nz(j))) ^ (1 / 3)
    t(j, k) = rt(j, k)
    If t(j, k) > dlt Then t(j, k) = dlt
    If Y(1, j, k) = 0 Then t(j, k) = 0.0000000001
   Next
  Next
  pc = s(3) * 1000 / (vd * w)
  auc = auc + pc * h
  dy = X / 1440
  If g < h Then file.Write(dy.ToString("0.00E+0"))
  If g < h Then file.Write(" ")
  If g < h Then file.WriteLine(pc.ToString("0.00E+0"))
  g = g + h
  If g > q - h / 2 Then g = 0
  For i = 1 To m
   For j = 1 To n
    For k = 1 To o
     eq(d, Xso, Y, p, t, r, Cs, V, ka, ba, cl, vd, k12, k21, f, j, k, s2, nwe)
     yn(i, j, k) = Y(i, j, k)
     kn1(i, j, k) = f(i, j, k)
     X = X + h / 2
     Y(i, j, k) = yn(i, j, k) + kn1(i, j, k) * h / 2
     If Y(1, j, k) < 1.0E-17 Then Y(1, j, k) = 0
     eq(d, Xso, Y, p, t, r, Cs, V, ka, ba, cl, vd, k12, k21, f, j, k, s2, nwe)
     kn2(i, j, k) = f(i, j, k)
     Y(i, j, k) = yn(i, j, k) + kn2(i, j, k) * h / 2
     If Y(1, j, k) < 1.0E-17 Then Y(1, j, k) = 0
     eq(d, Xso, Y, p, t, r, Cs, V, ka, ba, cl, vd, k12, k21, f, j, k, s2, nwe)
     kn3(i, j, k) = f(i, j, k)
     X = X + h / 2
     Y(i, j, k) = yn(i, j, k) + kn3(i, j, k) * h
     If Y(1, j, k) < 1.0E-17 Then Y(1, j, k) = 0
     eq(d, Xso, Y, p, t, r, Cs, V, ka, ba, cl, vd, k12, k21, f, j, k, s2, nwe)
     kn4(i, j, k) = f(i, j, k)
     Y(i, j, k) = yn(i, j, k) + h * (kn1(i, j, k) + 2 * kn2(i, j, k) + 2 * kn3(i, j, k) + kn4(i, j, k)) / 6
     If Y(1, j, k) < 1.0E-17 Then Y(1, j, k) = 0
     X = X - h
    Next
   Next
  Next
  X = X + h
  If X < xmax + h Then GoTo Start
  file.Close()
 End Sub
 Sub eq(ByRef d As Double, ByRef Xso(,) As Double, ByRef Y(,,) As Double, ByRef p As Double,
        ByRef t(,) As Double, ByRef r() As Double, ByRef Cs As Double, ByRef V As Double,
        ByRef ka As Double, ByRef ba As Double, ByRef cl As Double, ByRef vd As Double,
        ByRef k12 As Double, ByRef k21 As Double, ByRef f(,,) As Double, ByRef j As Integer,
        ByRef k As Integer, ByRef s2() As Double, ByRef nwe As Double)
  nwe = 3 * d * Xso(j, k) ^ (1 / 3) * Y(1, j, k) ^ (2 / 3) / p / t(j, k) / r(j) * (Cs - s2(k) / V)
  f(1, j, k) = -nwe
  f(2, j, k) = nwe - ka * Y(2, j, k)
  f(3, j, k) = ka * Y(2, j, k) * ba - (cl / vd + k12) * Y(3, j, k) + k21 * Y(4, j, k)
  f(4, j, k) = k12 * Y(3, j, k) - k21 * Y(4, j, k)
 End Sub
End Module