c      deadtime in an n-deep fifo buffer     6/14/96

c       code set for maximum of n=10.

      real c(11), aa(101), dead(10,101)

      amax = 10.
      aa(1) = 0.001 ! amin
      nstep = 100
      astep = exp(alog(amax/aa(1))/float(nstep))

      do i = 1,nstep+1
        if (i .gt. 1) aa(i) = aa(i-1)*astep
c        print 10, i, aa(i)
   10   format (' i, a ',i4,f9.5)
      enddo ! i

      do n = 1,6 ! loop on number of buffers
      
        do ia = 1,nstep+1 ! loop on a = t_0/t_1
          a = aa(ia)

c         calculate the coefficents c(i)
          do i = 0,n-1
            sum = 0.

            do j = 0,i
              product = 1.
              if (j .le. i-1) then

                do k = j,i-1
                  product = product*c(k+1)     
                enddo ! k

              endif ! j
              sum = sum + q(i,j,n,a)*product
            enddo ! j

            c(i+1) = q(i,i+1,n,a)/(1. - sum)
          enddo ! i
       
          sum = 0.
          do i = 0,n ! calculate the deadtime denominator

            product = 1.
            if (i .lt. n) then
              do j = i,n-1
                product = product*c(j+1)
              enddo ! j

            endif ! i

          sum = sum + product
          enddo ! i

        dead(n,ia) = 1./sum
c         set tiny negative values to zero
        if (dead(n,ia) .lt. 0.) dead(n,ia) = 0.  
        enddo ! ia
      enddo ! n          

c       print out

      do ia = 1, nstep+1
        write(1,100) aa(ia), (dead(n,ia),n=1,6)
  100   format (f9.5,',',e10.3,',',e10.3,',',e10.3,',',e10.3,','
     .          ,e10.3,',',e10.3,',')
      enddo ! ia

      stop
      end

      function q(i,j,n,a)

c         factorial(i) = (i-1)!  or conversely, i! = factorial(i+1)
      real factorial(11) / 1., 1., 2., 6., 24., 120., 720., 5040.,
     .                     40320., 362880., 3628800. /

      q = 0.
      if (i .ge. n) return  ! skip case of i = n

      if (j .eq. 0) then
        q = a**i*exp(-a)/factorial(i+1)
      elseif (j .gt. i+1) then
        q = 0.
      elseif (j .le. i+1 .and. i .lt. n) then   ! does {i,j} = {n-1,n} wrong
        q = a**(i-j+1)*exp(-a)/factorial(i-j+2)
      endif ! j

      if (j .eq. n .and. i .eq. n-1) q = (1. - exp(-a))/a ! fixup this case

      return
      end