      PROGRAM integ
*     compute integral of f(x) from x = a to x = b
*     using rectangular approximation
      IMPLICIT NONE
      REAL*8 a, b, delta, dx, sum
      INTEGER k, n, kbhit, KBGETCH
      LOGICAL Ltmp1_
      CALL initial(a,b,sum,dx,delta,n)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL sumf(a,b,sum,dx,delta)
         CALL output(sum,dx,n)
         CALL double(dx,delta,n)
         Ltmp1_ = kbhit() .EQ. 0
      END DO
      k = KBGETCH(0)
      END
C
      SUBROUTINE initial(a,b,sum,dx,delta,n)
      IMPLICIT NONE
      REAL*8 a, b, sum, dx, delta, f
      INTEGER n
      WRITE(*,'(A,$)') 'lower limit a = '
      READ(*,*) a
      WRITE(*,'(A,$)') 'upper limit b = '
      READ(*,*) b
      n = 2                                 ! initial number of intervals
      dx = (b - a)/n
      delta = dx                            ! delta = dx only for n = 2
      sum = f(a)                            ! -> (f(a) + f(b))/2 for trapezoidal rule
      WRITE(*,*)
      WRITE(*,*) '       n rectangular approximation'
      WRITE(*,*)
      END
C
      SUBROUTINE sumf(a,b,sum,dx,delta)
      IMPLICIT NONE
      REAL*8 x, a, b, sum, dx, delta, f
      x = a + dx                            ! dx is integration interval
      DO WHILE(x .LT. b)
         sum = sum + f(x)
         x = x + delta
      END DO
      END
C
      SUBROUTINE double(dx,delta,n)
*     double number of intervals
      IMPLICIT NONE
      REAL*8 dx, delta
      INTEGER n
      delta = dx
      n = 2*n
      dx = 0.5D0*dx                         ! dx does not equal delta for n > 2
      END
C
      SUBROUTINE output(sum,dx,n)
      IMPLICIT NONE
      REAL*8 sum, dx
      INTEGER n
      WRITE(*, '(I9,4X,F12.7)') n, sum*dx
      END
C
      FUNCTION f(x)
      IMPLICIT NONE
      REAL*8 f, x
      f = cos(x)
      END
