! React, version 1a.
! John Wise, initiated 2001/01/12, last updated 2001/01/12.
! Based on "Decay" in "Computational Physics" by N. Giordano.
program React1
   library "IrGraphics*"
   option nolet
   dim t(0), na(0), nb(0)
   call Init(nt, dt, t0, na0, nb0, ra, rb)
   mat redim t(0:nt), na(0:nt), nb(0:nt)
   t(0) = t0
   na(0) = na0
   nb(0) = nb0
   call Calc(nt, dt, t(), na(), nb(), ra, rb)
   id$ = "JHW: React1. Dashed line : N-a sample (ra:2). Solid Line : N-b sample (rb:5)."
   call Graph2A(id$, "Time", "Number", t(), na(), 2, t(), nb(), 1)
   print " t = "; t(nt)
   print " na = "; na(nt)
   print " nb = "; nb(nt)
   print "Done"
end
sub Init(nt, dt, t0, na0, nb0, ra, rb)
   nt = 50           ! number of time steps
   dt = 0.01          ! time step
   t0 = 0             ! initial time
   na0 = 1             ! initial value
   nb0 = 0             ! initial value
   ra = 2               ! decay rate
   rb = 5               ! decay rate
end sub
sub Calc(nt, dt, t(), na(), nb(), ra, rb)
   for i = 1 to nt
      nd1 = - ra*na(i-1)
      nd2 = ra*na(i-1) - rb*nb(i-1)
      t(i) = t(i-1) + dt
      na(i) = na(i-1) + nd1*dt
      nb(i) = nb(i-1) + nd2*dt
   next i
end sub