Imports System.Math
Imports System.IO
Module Module1
 Sub Main()
  Dim auc As Double = 0
  Dim ba As Double = 0.45
  Dim cl As Double = 4.0
  Dim Cs As Double = 0.01
  Dim d As Double = 0.0003
  Dim dlt As Double = 0.003
  Dim dose As Double = 0
  Dim f(10, 100, 1500) As Double
  Dim g As Double = 0
  Dim h As Double = 0.01
  Dim i, j, k As Integer
  Dim ka As Double = 0.07
  Dim kn1(10, 100, 1500), kn2(10, 100, 1500), kn3(10, 100, 1500), kn4(10, 100, 1500) As Double
  Dim k12 As Double = 0.028
  Dim k21 As Double = 0.0088
  Dim m As Integer = 4
  Dim mz(100) As Double
  Dim n As Integer = 1
  Dim nwe As Double
  Dim nz(100, 1500) As Double
  Dim p As Double = 1300
  Dim pc As Double
  Dim pd As Double
  Dim q As Double = 10
  Dim r(100) As Double
  Dim rt(100, 1500) As Double
  Dim s(10) As Double
  Dim t(100, 1500) As Double
  Dim V As Double = 250
  Dim vd As Double = 600
  Dim vz(100) As Double
  Dim X As Double = 0
  Dim xmax As Double = 1440
  Dim Xso(100, 1500) As Double
  Dim Y(10, 100, 1500) As Double
  Dim yn(10, 100, 1500) As Double
  Dim file As StreamWriter
  file = My.Computer.FileSystem.OpenTextFileWriter(
            "C:\Users\Kevin C. Johnson\cr.txt", False)
  Dim w As Double = 75.3
  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
  Dim cr As Double = 0
  Dim crc As Double = 0
  Dim rr As Double = 1.7 / 60
  Dim tdose As Double = 0
  Dim o As Integer = xmax
  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
  For i = 1 To m
   For j = 1 To n
    For k = 0 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 j = 1 To n
   If cr = 0 Then Xso(j, crc) = freq(j) / tfreq * rr
   If cr = 0 Then Y(1, j, crc) = Xso(j, crc)
   If cr = 0 Then nz(j, crc) = Y(1, j, crc) / mz(j)
  Next
  For k = 0 To o
   For i = 1 To n
    If nz(i, k) > 0 Then rt(i, k) = (3 * Y(1, i, k) / (4 * PI * p * nz(i, k))) ^ (1 / 3)
    If nz(i, k) > 0 Then t(i, k) = rt(i, k)
    If t(i, k) > dlt Then t(i, k) = dlt
    If Y(1, i, k) = 0 Then t(i, k) = 0.0000000001
   Next
  Next
  For i = 1 To m
   For j = 1 To n
    For k = 0 To o
     s(i) = s(i) + Y(i, j, k)
    Next
   Next
  Next
  If cr = 0 And X >= 1 - h / 2 Then tdose = tdose + rr
  pc = s(3) * 1000000 / (vd * w)
  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(pc.ToString("0.00E+0"))
  If g < h Then Console.WriteLine(X / xmax * 100)
  g = g + h
  If g > q - h / 2 Then g = 0
  cr = cr + h
  If cr > 1 - h / 2 Then crc = crc + 1
  If cr > 1 - h / 2 Then cr = 0
  For i = 1 To m
   For j = 1 To n
    For k = 0 To o
     eq(d, Xso, Y, p, t, r, Cs, s, V, ka, ba, cl, vd, k12, k21, f, j, k, 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, s, V, ka, ba, cl, vd, k12, k21, f, j, k, 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, s, V, ka, ba, cl, vd, k12, k21, f, j, k, 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, s, V, ka, ba, cl, vd, k12, k21, f, j, k, 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 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 k As Integer, 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 - s(2) / 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