#include <stdio.h>
#include <math.h>
double a,b,o;
int n;

/* integrate the function f(x) numerically between limits x=a and x=b,
 * using the trapezium rule with n steps. Returns the calculated integral.
 */
double integrate ( double (*f)(double), double a, double b, int n )
{
	int x;
	double d;
	/* The first part of the equation, it halfs the function of
	 * a added to the function of b
	 */
	o = 0.5*(f(a) + f(b));
	/* This part works out d which is the strip seperation */
	d = (b - a) / n;
	/* This part adds on the the sum part of the equation */
	for (x = 1; x < n; x++)
	{
		o = o + f(a + (x*d));
	}
	/* This part adds on Simpsons part of the sum */
	for (x = 1; x <= n; x++)
	{
		o = o + (2 * f(0.5*( (a + (x*d)) + (a + ((x-1)*d) )) ));
	}
	/* This part just the whole brackets by d */
	o = (d/3) * o;
	/* This returns the answer */
	return o;
}

/* return value of x*sin(x*x) given x
 */
double formula(double x)
{
	/* This will work out the function, sin is included in
	 * the maths library
	 */
	x = x*sin(x*x);
	return x;
}

int main ( int argc, char *argv[] )
{
	/* Print the welcome message into standard error and take in
	 * the values for a,b and n. Where a and b are double format
	 */
	fprintf(stderr, "Welcome to Integration by the trapezium rule\n");
	fprintf(stderr, "Please enter the value for a :-\n");
	scanf("%lf", &a);
	fflush(stdin);
	fprintf(stderr, "Please enter the value for b :-\n");
	scanf("%lf", &b);
	fflush(stdin);
	fprintf(stderr, "Please enter the value for n :-\n");
	scanf("%i", &n);
	fflush(stdin);
	/* Run the integral program */
	o = integrate(formula,a,b,n);
	/* Print out the values a,b and n to standard error for\
	 * a check
	 */
	fprintf(stderr, "a - %f, b - %f, n - %i\n", a, b, n);
	/* Print out the answer */
	printf("%f\n",o);
    system("PAUSE");
	return 0;
}