      PROGRAM AREA
************************************************************************
* Program to approximate the integral of a function over the interval  *
* [A,B] using the trapezoidal method. This approximation is calculated *
* by the function subprogram DEFINT; the integrand, the interval of    *
* integration, and the # of subintervals are passed as arguments to    *
* DEFINT.  Identifiers used are:                                       *
*     A, B   : the endpoints of the interval of integration            *
*     SIN    : the integrand (library function)                        *
*     NSUBS  : the number of subintervals used                         *
*                                                                      *
* Input:  A, B, and NSUBS                                              *
* Output: Approximation to integral of F on [A, B]                     *
************************************************************************

      REAL A, B, DEFINT
      INTEGER NSUBS
      INTRINSIC SIN

      PRINT *, 'ENTER THE INTERVAL ENDPOINTS AND THE # OF SUBINTERVALS'
      READ *, A, B, NSUBS

      PRINT 10, NSUBS, DEFINT(SIN, A, B, NSUBS)
 10   FORMAT (1X, 'APPROXIMATE VALUE USING', I4, 
     +            ' SUBINTERVALS IS', F10.5)

      END


** DEFINT **************************************************************
* Function to calculate the trapezoidal approximation of the integral  *
* of the function F over the interval [A,B] using N subintervals.      *
* Local variables used are:                                            *
*     I    : counter                                                   *
*     DELX : the length of the subintervals                            *
*     X    : a point of subdivision                                    *
*     Y    : the value of the function  at X                           *   
* Accepts: Function F, endpoints A and B, and number N of subintervals *
* Returns: Approximate value of integral of F over [A, B]              * 
************************************************************************

      FUNCTION DEFINT(F, A, B, N)

      INTEGER N, I
      REAL DELX, X, DEFINT, F, A, B

* Calculate subinterval length 
* and initialize the approximating sum and X

      DELX = (B - A) / REAL(N)
      X = A 
      DEFINT = 0.0

* Now calculate the approximating sum

      DO 10 I = 1, N - 1
         X = X + DELX
         Y = F(X)
         DEFINT = DEFINT + Y
 10   CONTINUE
      DEFINT = DELX * ((F(A) + F(B)) / 2.0 + DEFINT)

      END

**POLY******************************************************************
*                         The integrand                                *
************************************************************************

      FUNCTION POLY(X)

      REAL X, POLY
      POLY = X ** 2 + 1

      END
