/*
 *  ==========================================================================
 *  prog_trapezoid.c : Approximate integral of a function via Trapezoidal Rule
 *
 *  Copyright (C) 1994-96 by Mark Austin and David Mazzoni.
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software on any computer system,
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 *
 *  Written By: M. Austin                                         January 1994  
 *  ==========================================================================
 */

#include <stdio.h>

double Trapezoid( double (*Func)( double ), float, float, int );    
double MySquare( double );                                              
double MyCubic( double );

int main( void ) {
double dIntegral;

    dIntegral = Trapezoid( MySquare, 0.0, 10.0, 100);
    printf("*** INTEGRAL OF (X*X)   from 0 to 10 by 0.1 = %10.4f\n", dIntegral);

    dIntegral = Trapezoid( MyCubic, 0.0, 10.0, 100);
    printf("*** INTEGRAL OF (X*X*X) from 0 to 10 by 0.1 = %10.4f\n", dIntegral);
}

/*
 *  =====================================
 *  Compute integral via Trapezoidal-Rule 
 *  =====================================
 */

double
Trapezoid( double (*Func)( double ), float fMinValue,
           float  fMaxValue,         int   iNoIntervals ) {
double dSum = 0.0;
double dIncrement, dX;

      dIncrement = ( fMaxValue - fMinValue )/iNoIntervals;
      for(dX = fMinValue; dX <= fMaxValue - dIncrement; dX += dIncrement)
          dSum += (dIncrement/2.0)*( Func(dX) + Func(dX+dIncrement));

      return ( dSum );
}

/*
 *  ========================================== 
 *  User-Defined Quadratic and Cubic Functions
 *  ========================================== 
 */

double MySquare( double dX ) {
   return ( dX*dX );
}

double MyCubic( double dX ) {
   return ( dX*dX*dX );
}