/* 
   Name:           AdaptiveQuadraturSimpson  
   Author:         Manfred  Vogel
   Description:    Adaptives Quadratur-Verfahren mit Simpson-Regel
   Date:           27. Januar 2005
   Copyright: 
*/

//#include "stdafx.h"
#include "iostream.h"
#include <stdlib.h>
#include <math.h>

#define Pi 3.1415926535

// Integrationsintervall [a, b]
#define a0 0.0
#define b0 10.0

// # Teilintervalle für die Auswertung des Integrals
// mit der zusammengesetzten Simpson-Regel
#define n 4


// definiert die zu integrierende Funktion f
double f(double x)
{
	if( x==0.0) return 0.0;
	return sin(x*x*x) * exp(-x/exp(1.0));
	// return 0.0;
}

// implementiert die zusammengesetzte Simpson'sche Regel
double simpson(double a, double b)
{
	double h = (b-a)/n;
	double sum1=0.0, sum2=0.0;
	for (int i=1; i<=(n/2)-1; i++)
	{
		sum1 += f(a + 2*h*i);
	}
	for (int j=1; j<=n/2; j++)
	{
		sum2 += f(a + h*(2*j-1));
	}
	return h/3 * ( f(a) + f(b) + 2*sum1 + 4*sum2 );
	// return 0.0;	
}

double rek(double a, double b, double eps)
{
	double intTot=0.0, intLeft=0.0, intRight=0.0;
	intTot   = simpson( a      , b      );
	intLeft  = simpson( a      , (a+b)/2);
	intRight = simpson( (a+b)/2, b      );
	if( fabs( intTot - intLeft - intRight ) > eps )
	{       
		// System.out.println("a: " + a + "   b: " + b);
		return rek(a, (a+b)/2, eps/2.0) + rek((a+b)/2, b, eps/2.0);
	}
	else
	{
		// return intTot;
		return intLeft + intRight;	
	}
	//return 0.0;
}               

void main()
{
	double integral = rek(a0, b0, 0.00000001);
	
	cout.precision(14);
	cout << "\n\n" << endl;
	cout << " ** Adaptive Quadratur mit Simpson-Verfahren ** " << "\n\n";

	cout << "    Integral( f(x), x=a..b ) = ";
	cout << integral << endl;

	cout << "\n\n" << endl;

	system("Pause");
}