program micros
      implicit none
      save

c..educational version of an s-process network. 
c..one hundred isotopes are included and the initial abundance of the first, 
c..nominally fe56, is taken equal to one. these abudances are evolved
c..until the user chosen ending time. 
c..
c..here is a little theory on what this program does:
c.. 
c..let ratng be the reaction rate for (n,g) reactions.
c..then the ordinary differentail equations describing the change
c..in the abundances y are:
c..dydt(1) = -ratng(1) * y(1)
c..dydt(i) =  signng(i-1) * y(i-1) - ratng(i) * y(i)   i=2,n-1
c..dydt(n) =  signng(n-1) * y(n-1)
c..
c..for the first-order accurate euler method, 
c..each abundance is updated over a timestep h as 
c..ynew(i) = y(i) + dy(i) 
c..where the change dy(i) is obtained from solving the linear system
c.. (I/h - J) = f(i)
c..
c..with only (n,g) reactions, the jacobian matrix J has the special form
c..
c..| -ratng(1)                                                 |
c..|  ratng(1)  -ratng(2)                                      |
c..|             ratng(2)    -ratng(3)                         |
c..|                                     ...                   |
c..|                                           ratng(n-1)  0.0 |
c..
c..this system of linear equations can be easily solved by hand:
c.. dy(1) = -y(1)*ratng(1)/(dt + ratng(1)  
c.. dy(i) = (-y(i)*ratng(i) + (y(i-1) + dy(i-1))*ratng(i-1))/(1/h + ratng(i))
c.. dy(n) = y(n-1)*ratng(n-1)*h
c..
c..hence the very succint updating scheme implemented below.



c..declare
      integer          i,j,jmax,nmax
      parameter        (jmax=10000, nmax=100)
      double precision y(nmax),dy(nmax),ratng(nmax),tend,dt,dtinv,t,xx


c..time step controllers
c..chimin is smallest abundance that can affect the time step
c..delchi is maximum fractional change allowed 
c..fdtn is the maximum enlargement factor for new time step

      double precision chimin,delchi,fdtn
      parameter        (chimin = 1.0d-3, 
     1                  delchi = 0.05d0, 
     2                  fdtn   = 2.0d0)


c..popular format statements
 01   format (1x,a,1pe11.2,/,1x,a,a)
 02   format (i5,1pe11.2)



c..initialize, this effectively says all the cross sections are equal.
c..as an exercise in bottlenecks, try setting a few not equal to one.
 10   do i=1,nmax
       y(i)     = 0.0d0 
       ratng(i) = 1.0d0
      enddo
      y(1) = 1.0d0


c..get the ending time. note the ending time may be interpreted as the 
c..exposure tau, the integral of the neutron density, velocity, and time.
      write (6,*) ' give ending time (less than zero to stop) =>'
      read (5,*) tend
      if (tend.le.0.) stop ' normal termination'


c..initial time, ending time and timestep
      t    = 0.0d0
      dt   = 1.0d-6*tend


c..start the evolution
      do j=1,jmax


c..increment the total time
       t = t + dt


c..here is the very succint update scheme
       dtinv = 1.0d0/dt
       dy(1) = -y(1)*ratng(1)/(dtinv + ratng(1))
       y(1)  = y(1) + dy(1) 
       do i=2,nmax-1
        dy(i) = (-y(i)*ratng(i) + y(i-1)*ratng(i-1))/(dtinv + ratng(i))
        y(i)  = y(i) + dy(i)
       enddo
       dy(nmax) = y(nmax-1)*ratng(nmax-1)*dt
       y(nmax)  = y(nmax) + dy(nmax)


c..this is the normal exit point
       if ( t-tend .ge. 0.0d0  .or.  j .eq. jmax ) then
        write (6,01) 'time=',t,'   i', '   abundance'
        write (6,02) (i,y(i), i=1,nmax)
        goto 10
       end if


c..get a new timestep
       xx = fdtn*dt
       do i = 1,nmax
        if (y(i) .gt. chimin) then
         xx = min(xx,abs(y(i)/dy(i))*delchi*dt)
        endif
       enddo
       dt = xx

c..if the step can overshoot the stop point, cut it down
       if ( t+dt-tend .gt. 0.0d0) dt = tend - t


c..back for another update with this new timestep
      enddo 

      return
      end