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.0
  Dim Cs As Double = 0.00626
  Dim d As Double = 0.000404
  Dim dlt As Double = 0.003
  Dim dose As Double = 5
  Dim f(10, 100) As Double
  Dim g As Double = 0
  Dim h As Double = 0.0001
  Dim i, j As Integer
  Dim ka As Double = 0.0
  Dim kn1(10, 100) As Double
  Dim kn2(10, 100) As Double
  Dim kn3(10, 100) As Double
  Dim kn4(10, 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 = 19
  Dim nwe As Double
  Dim nz(100) As Double
  Dim p As Double = 1260
  Dim pc As Double
  Dim pd As Double
  Dim q As Double = 0.1
  Dim r(100) As Double
  Dim rt(100) As Double
  Dim s(10) As Double
  Dim t(100) As Double
  Dim V As Double = 900
  Dim vd As Double = 1000
  Dim vz(100) As Double
  Dim X As Double = 0
  Dim xmax As Double = 60
  Dim Xso(100) As Double
  Dim Y(10, 100) As Double
  Dim yn(10, 100) As Double
  Dim realps As StreamReader
  realps = My.Computer.FileSystem.OpenTextFileReader(
   "C:\Users\Kevin C. Johnson\realpsj.txt")
  Dim file As StreamWriter
  file = My.Computer.FileSystem.OpenTextFileWriter(
   "C:\Users\Kevin C. Johnson\jinnoj.txt", False)
  For i = 1 To n
   r(i) = realps.ReadLine() / 20000
   Xso(i) = realps.ReadLine() / 99.974 * dose
   vz(i) = 4 / 3 * PI * r(i) ^ 3
   mz(i) = vz(i) * p
   Y(1, i) = Xso(i)
   Y(2, i) = 0
   Y(3, i) = 0
   Y(4, i) = 0
   nz(i) = Y(1, i) / mz(i)
  Next
Start:
  For i = 1 To n
   rt(i) = (3 * Y(1, i) / (4 * PI * p * nz(i))) ^ (1 / 3)
   t(i) = rt(i)
   If t(i) > dlt Then t(i) = dlt
   If Y(1, i) = 0 Then t(i) = 0.0000000001
  Next
  If X > 0.1 Then h = 0.001
  If X > 1 Then h = 0.01
  For i = 1 To m
   s(i) = 0
  Next
  For i = 1 To m
   For j = 1 To n
    s(i) = s(i) + Y(i, j)
   Next
  Next
  pd = s(2) / dose * 100
  pc = s(3) * 1000 / vd
  auc = auc + pc * h
  If g < h Then file.Write(X.ToString("0.00E+0"))
  If g < h Then file.Write(" ")
  If g < h Then file.WriteLine(pd.ToString("0.0000E+0"))
  g = g + h
  If g > q - h / 2 Then g = 0
  For i = 1 To m
   For j = 1 To n
    eq(d, Xso, Y, p, t, r, Cs, s, V, ka, ba, cl, vd, k12, k21, f, j, nwe)
    yn(i, j) = Y(i, j)
    kn1(i, j) = f(i, j)
    X = X + h / 2
    Y(i, j) = yn(i, j) + kn1(i, j) * h / 2
    If Y(1, j) < 1.0E-17 Then Y(1, j) = 0
    eq(d, Xso, Y, p, t, r, Cs, s, V, ka, ba, cl, vd, k12, k21, f, j, nwe)
    kn2(i, j) = f(i, j)
    Y(i, j) = yn(i, j) + kn2(i, j) * h / 2
    If Y(1, j) < 1.0E-17 Then Y(1, j) = 0
    eq(d, Xso, Y, p, t, r, Cs, s, V, ka, ba, cl, vd, k12, k21, f, j, nwe)
    kn3(i, j) = f(i, j)
    X = X + h / 2
    Y(i, j) = yn(i, j) + kn3(i, j) * h
    If Y(1, j) < 1.0E-17 Then Y(1, j) = 0
    eq(d, Xso, Y, p, t, r, Cs, s, V, ka, ba, cl, vd, k12, k21, f, j, nwe)
    kn4(i, j) = f(i, j)
    Y(i, j) = yn(i, j) + h * (kn1(i, j) + 2 * kn2(i, j) + 2 * kn3(i, j) + kn4(i, j)) / 6
    If Y(1, j) < 1.0E-17 Then Y(1, j) = 0
    X = X - h
   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 s() 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 nwe As Double)
  nwe = 3 * d * Xso(j) ^ (1 / 3) * Y(1, j) ^ (2 / 3) / p / t(j) / r(j) * (Cs - s(2) / V)
  f(1, j) = -nwe
  f(2, j) = nwe - ka * Y(2, j)
  f(3, j) = ka * Y(2, j) * ba - (cl / vd + k12) * Y(3, j) + k21 * Y(4, j)
  f(4, j) = k12 * Y(3, j) - k21 * Y(4, j)
 End Sub
End Module