subroutine  integ_1(imode,a,b,dx0,acc,si)
C
C This routine calculates the definite (or indefinite integral) from a to 
C  b (or infinity) of a function called fun_1(x).
C Parameters:
C     
C  imode [in]:     if 0, calculate the definite integral \int_a^b. if not 0,
C                  calculate the indefinite integral \int_a^\infty.
C  a     [in]:     lower limit of the integration
C  b     [in/out]: upper limit of the integration (input) if
C                  imode=0. The maximum value of x (output) if imode!=0.
C  dx0   [in]:     a characteristic scale over which fun_1 changes
C                  appreciably. It will affect how fast the integral is
C                  calculated, but not the final success. if you don't
C                  know, set dx0 to a or 1. Should be positive.
C  acc   [in]:     relative accuracy of the integral
C  si    [out]:    value of the integral.
C
C *******************************************************
C
C  First, initialize variables:
C   dx is the current interval over which we integrate.
C   x1,x2 are current limits of the interval. x1 initially is a.
C   f1,f2 are function values at the limits
C
      dx = dx0
      si = 0.0
      x1 = a
      f1 = fun_1(x1)
      x2 = x1
C
C  Start the loop over all sub-intervals. if we reached b (for a definite
C  integral, then stop.
C
 10   if(imode.eq.0 .and. x2.ge.b) goto 100
C
C  get the upper limits of the current interval
C
      x2 = x1 + dx
C
C  make sure x2<=b if imode=0
C
      if (imode.eq.0 .and. x2.gt.b) x2 = b
      f2 = fun_1(x2)
C
C  calculate the integral between x1 and x2 using part_1
C
      call  part_1(x1,f1,x2,f2,acc,0.0,sn,iter)
C
C  if part_1 made only a few iterations (<3), the function is too
C  smooth over the interval dx. Then increase the size of the next
C  interval by 50%.
C
      if(iter .lt. 3) dx = dx*1.5
C
C  if part_1 made less than 10 iterations, it was successful.
C  Then add the integral from x1 to x2 (sn) to the total si.
C
      if(iter .lt. 10) then
         si = si + sn
C
C  If we are doing indefinite integral, check that the contribution
C  of sn to so is small so that we can stop. To be on the safe side,
C  I adopt 10 times more stringent tolerance limit than the one used to
C  compute sn.
C
         if(imode .ne. 0) then
            b = x2
            if(abs(sn) .le. 0.1*acc*abs(si)) goto 100
         endif
         x1 = x2
         f1 = f2
      else
C
C  If part_1 did 10 iterations, it did not finish. Thus, the interval is
C  too large and the function changes on a smaller scale. Decrease the
C  interval size by half and redo the integral.
C
         dx = 0.5*dx
         x2 = x1
C
C  If something is wrong, we can get a very small interval. If it
C  get smaller than 10^{-5} of the original one, stop the integration
C  and report the error.
C
         if(dx .lt. 1.0e-5*dx0) then
            write(6,*) 'dx is too small:',dx,x1,x2
            stop
         endif
      endif
C
C  Continue the loop over intervals
C
      goto 10
C
C  We are done!!! celebrate!
C
 100  continue
      
      return
      end
C
C
C
      subroutine part_1(x1,f1,x2,f2,acc,so,sn,iter)
C
C This routine calculates the definite integral from x1 to x2
C of a function called fun_1(x).
C Parameters:
C     
C  x1,x2 [in]:  lower and upper limits of the integration
C  f1,f2 [in]:  values of the function at the limits. Should be
C               supplied at the input
C  acc   [in]:  relative accuracy of the integral
C  so    [in]:  value of the integral over other intervals. The
C               accuracy acc is satisfied for the sum (so+sn), and
C               if so is larger than sn, sn is computed with smaller  
C               accuracy. Set so to zero if you need sn to precisely
C               acc accuracy.
C  sn    [out]: value of the computed integral
C  iter  [out]: number of iterations performed.
C
C *******************************************************
C
C  First, initialize variables:
C   n:  number of subdivisions of the interval, alway a power of 2
C   dx: width of the interval
C   h: current size of the subdivision, h*n = dx
C   sp: value of the integral with n/2 subdivisions
C   sf: sum of all function values, S^O_N+S^E_N in the notes.
C   ds: current tolerance, abs(sn/sp-1). Initially very large.
C
      n = 1
      iter = 0
      dx = x2-x1
      h = dx
      sn = 0.5*(f1+f2)*h
      sf = 0.0
      ds = 1.0 + acc
C
C  Loop over iterations. Stop if the accuracy is achieved (ds<acc) or
C  a number of iterations is 10.
C
 10   continue
C
C  Assign sn to sp to save previous value
C
      sp = sn
C
C  Update the iteration counter
C
      iter = iter + 1
C
C  Compute the sum minus S^O_{2N}. Use sf=S^O_N+S^E_N as S^E_{2N}
C
      sn = f2 + f1 + 2.0*sf
C
C  double n and halve h. m is n-1, the last odd point
C
      h = h*0.5
      n = n+n
      m = n-1
C
C  Compute s=S^O_{2N} as the sum over odd points.
C
      s = 0.0
      do i=1,m,2
         x = x1 + h*i      
         s = s + fun_1(x)
      enddo	
C
C  Now sn is the total sum
C
      sn = sn + 4.0*s
C
C  Update sf. Now sf = S^O_{2N} + S^E_{2N}, ready for the next iteration.
C
      sf = sf + s
C
C  Compute the integral as the sum times h divided by 3.
C
      sn = sn*h/3.0
C
C  Compute the tolerance as the relative difference between so+sp and so+sn.
C
      if(iter.lt.10 .and. abs(sp-sn).gt.acc*abs(so+sn)) goto 10
C
C  We come here when we are done. 
C
 100  continue    

      return
      end