/* --------------------------------------------------------------------- **
   Figure 6.28:  Integration by Simpson's rule.
   Figure 6.29:  Summing areas by Simpson's rule.
** --------------------------------------------------------------------- */

#include "tools.h"                                                  
#define N  100

double simpson( double lower, double upper );
/* ---------------------- This is the function that we will integrate. */
double f (double x) {  return x * x + x;  }                                         
    
void main()                                                               
{
    double lower, upper;    /* Lower and upper limits of integration. */           
    double swap;            /* For interchanging limits, if necessary. */
    double area;            /* Result of integration. */

    banner();                                                   
    puts( "\n Integration by Simpson's rule.\n" ); 

    printf( " The function we will integrate is:  f(x)=x*x+x \n"    
            " Please enter the lower and upper bounds for integration: " );                                         
    scanf ( "%lg%lg", &lower, &upper );

    if ( upper <  lower ) { swap=lower; lower=upper; upper=swap; }                                     
    printf( "\n Integrating f(x)=x*x+x from %g to %g\n", lower, upper );

    area = simpson( lower, upper );
    printf( " Area = %g \n", area);
    bye();                                       
}

/* --------------------------------- Integration by Simpson's Method  */
double 
simpson( double low, double high )
{
    double h;               /* Interval for integration. */          
    double odds, evens;     /* Sums for odd and even intervals. */ 
    double x;               /* Evaluate the function at this point. */
    int k;                  /* Loop counter. */                                                     
    double area;            /* Weighted sums of all intervals. */

    h = (high - low) / N;

    /* Compute the odd sum; k takes on the odd values 1,3,5,7....*/    
    /* Maximum k value is number of intervals - 1 */                    
    odds = 0.0;                 /* initial value of odd sum      */       
    for (k = 1; k < N; k += 2) {
        x = low + k * h;                                                        
        odds += f( x );                                                          
    }                         
 
    /* Compute the even sum; k takes on even values 2, 4,6,8.... */ 
    /* Maximum value for k is number of intervals - 2            */              
    evens = 0.0;                   /* initial value of even sum  */    
    for (k = 2; k < N; k += 2 ) {
        x = low + k * h;                                                       
        evens += f( x );                                           
    } 

    area = (high-low) / (3 * N) *
           (f( low ) + f( high ) + 4 * odds + 2 * evens);   
    return area;      
}