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
  Dim Cs(2) As Double : Cs(1) = 0.1 : Cs(2) = 0.01
  Dim d As Double = 0.0003
  Dim dlt As Double = 0.003
  Dim dose As Double = 10
  Dim f(10, 100) As Double
  Dim g As Double = 0
  Dim h As Double = 0.0001
  Dim i, k 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(2) As Double
  Dim ms As Double = 0
  Dim nwe As Double
  Dim nz(2) As Double
  Dim p As Double = 1300 : Dim pp As Integer = 0
  Dim pc As Double
  Dim pd As Double
  Dim q As Double = 0.1
  Dim r(2) As Double : Dim di(2)
  Dim rt(2) As Double
  Dim s(10, 2) As Double : Dim s1t As Double : Dim s2t As Double = 0 : Dim s3t As Double = 0
  Dim t(2) As Double
  Dim V As Double = 250
  Dim vd As Double = 100
  Dim vz(2) As Double
  Dim X As Double = 0
  Dim xmax As Double = 60
  Dim Xso(2) As Double
  Dim Y(10, 2) As Double
  Dim yn(10, 2) As Double
  Dim file As StreamWriter
  file = My.Computer.FileSystem.OpenTextFileWriter(
   "C:\Users\Kevin C. Johnson\precip.txt", False)
  Dim w As Double = 70
  Dim dmean(2) As Double : dmean(1) = 0.001 : dmean(2) = 0.0001
  Dim rmean(2) As Double : rmean(1) = dmean(1) / 2 : rmean(2) = dmean(2) / 2
  For k = 1 To 2
   r(k) = rmean(k)
   If r(k) < dlt Then t(k) = r(k) Else t(k) = dlt
   vz(k) = 4 / 3 * PI * r(k) ^ 3
   mz(k) = vz(k) * p
  Next
  For i = 1 To m
   For k = 1 To 2
    Y(i, k) = 0
    Xso(k) = 0
   Next
  Next
  Xso(1) = dose
  Y(1, 1) = Xso(1)
  nz(1) = Y(1, 1) / mz(1)
Start:
  rt(1) = (3 * Y(1, 1) / (4 * PI * p * nz(1))) ^ (1 / 3)
  t(1) = rt(1)
  If t(1) > dlt Then t(1) = dlt
  If Y(1, 1) = 0 Then t(1) = 0.0000000001
  ms = s2t / V / Cs(2)
  If ms > 1.5 And pp = 0 Then precip(Xso, dose, Y, nz, mz, pp)
  If pp = 1 Then
   rt(2) = (3 * Y(1, 2) / (4 * PI * p * nz(2))) ^ (1 / 3)
   t(2) = rt(2)
   If t(2) > dlt Then t(2) = dlt
   If Y(1, 2) = 0 Then t(2) = 0.0000000001
  End If
  For i = 1 To m
   For k = 1 To 2
    s(i, k) = 0
   Next
  Next
  For i = 1 To m
   For k = 1 To 2
    s(i, k) = s(i, k) + Y(i, k)
   Next
  Next
  s1t = s(1, 1) + s(1, 2)
  s2t = s(2, 1) + s(2, 2)
  s3t = s(3, 1) + s(3, 2)
  pc = s3t * 1000 / (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.Write(s1t.ToString("0.000E+0"))
  If g < h Then file.Write(" ")
  If g < h Then file.Write(s2t.ToString("0.000E+0"))
  If g < h Then file.Write(" ")
  If g < h Then file.WriteLine(s3t.ToString("0.000E+0"))
  g = g + h
  If g > q - h / 2 Then g = 0
  For i = 1 To m
   For k = 1 To 2
    eq(d, Xso, Y, p, t, r, Cs, s2t, V, ka, ba, cl, vd, k12, k21, f, i, k, nwe)
    yn(i, k) = Y(i, k)
    kn1(i, k) = f(i, k)
    X = X + h / 2
    Y(i, k) = yn(i, k) + kn1(i, k) * h / 2
    If Y(1, k) < 1.0E-17 Then Y(1, k) = 0
    eq(d, Xso, Y, p, t, r, Cs, s2t, V, ka, ba, cl, vd, k12, k21, f, i, k, nwe)
    kn2(i, k) = f(i, k)
    Y(i, k) = yn(i, k) + kn2(i, k) * h / 2
    If Y(1, k) < 1.0E-17 Then Y(1, k) = 0
    eq(d, Xso, Y, p, t, r, Cs, s2t, V, ka, ba, cl, vd, k12, k21, f, i, k, nwe)
    kn3(i, k) = f(i, k)
    X = X + h / 2
    Y(i, k) = yn(i, k) + kn3(i, k) * h
    If Y(1, k) < 1.0E-17 Then Y(1, k) = 0
    eq(d, Xso, Y, p, t, r, Cs, s2t, V, ka, ba, cl, vd, k12, k21, f, i, k, nwe)
    kn4(i, k) = f(i, k)
    Y(i, k) = yn(i, k) + h * (kn1(i, k) + 2 * kn2(i, k) + 2 * kn3(i, k) + kn4(i, k)) / 6
    If Y(1, k) < 1.0E-17 Then Y(1, k) = 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 s2t 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 i As Integer, ByRef k As Integer, ByRef nwe As Double)
  nwe = 3 * d * Xso(k) ^ (1 / 3) * Y(1, k) ^ (2 / 3) / p / t(k) / r(k) * (Cs(k) - s2t / V)
  f(1, k) = -nwe
  f(2, k) = nwe - ka * Y(2, k)
  f(3, k) = ka * Y(2, k) * ba - (cl / vd + k12) * Y(3, k) + k21 * Y(4, k)
  f(4, k) = k12 * Y(3, k) - k21 * Y(4, k)
 End Sub
 Sub precip(ByRef Xso() As Double,
            ByRef dose As Double, ByRef Y(,) As Double, ByRef nz() As Double,
            ByRef mz() As Double, ByRef pp As Integer)
  Xso(2) = dose / 1000
  Y(1, 2) = Xso(2)
  nz(2) = Y(1, 2) / mz(2)
  pp = 1
 End Sub
End Module