/******************************************************************
****Program to compute the integral of x^3sin(x)e^(-x) between the 
limits 0 to 1000. This program uses the simpson's 3/8th rule to compute
the integral and does it to an interval of 6000. Each interval has 2 bins
so the total number of bins can go upto 12000.
Functions: simpsonthree declared implicitly and used as a separate
piece of code
Author: Rakesh Nath
*******************************************************************/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<float.h>
#include<string.h>
#include<file.h>
#include<omp.h>
//#include<thread.h>
/*Precompiler definitions*/
#define N 6000/*Number of interval N=2xnumber of bins*/
#define ERR 1.0e-5//error tolerance 
/*formal declaration of the function */
double simpsonthree(double (*func)(double),double,double,int);
/*************************************************************
Function: f
Description: Function merely returns double value of the math function
declared within it.
Returns: double
Input: Double values which will be supplied by the main program
**************************************************************/

double f(double x)
	{	
		return pow(x,3)*sin(x)*exp(-x);
		//return sin(x);	
		//return exp(pow(x,2));
	}
/**************************************************************
Function Name:main()
Description: The main program does the integration using Simpson's 3/8th rule 
by calling the simpsonthree() function. The error values are calculated using
relative error and are saved to a file named error_values.dat. 
Returns: If successful 0
Input:none
***************************************************************/
int main()
{	FILE *fp;
	/*file pointer fp*/
	double *dSim= malloc((N+1)*sizeof(double));
	/*size is allocated based on the number of intervals, 
this will have the array of values of the result of the 
simpson's integration.*/
	double iLowerLimit,iUpperLimit;
	//upper and lower limits
	int i;
	//loop variable
	double dRes,dErr[N],dFinres;
	//results, errors and intermediate results

	iLowerLimit=0;
	iUpperLimit=1000;
/*The lower limit is set to 0 and the upper limit to a 1000*/

	fp=fopen("error_values.dat","w+");
	//open the file error_values.dat
	for(i=1;i<N;i+=1)
{	
		dSim[i]=simpsonthree(f,iLowerLimit,iUpperLimit,i);
//the first set of values computed with an interval i
		dFinres=simpsonthree(f,iLowerLimit,iUpperLimit,i-1);
//the preveious value with interval i-1
		dErr[i]=fabs((dSim[i]-dFinres)/dSim[i]);
//error calculation
		fprintf(fp,"%e\n",dErr[i]);
/*if the error is less than the 10^-5 then print, 
this is the convergence of the integral*/
		if(dErr[i]<=ERR)
		{	
		printf("at i=%d\t",i);
		printf("The value using Simpson's 3/8th rule is %e\n",dSim[i]);
		printf("the value of error is %e\n",dErr[i]);
		//break;
		}
	}

fclose(fp);
free(dSim);
	return 0;
}