#include <math.h>
#include <stdio.h>
#include <assert.h>

/* useful constants: pi, 2*pi, pi/2 */
#define PI   3.14159265358979323846
#define PI2  6.28318530717958647692
#define PI_2 1.57079632679489661923
/* The limit of fourier series elements to calculate */
#define max 100
int nsigma;

/* Defines the function f(x) to be integrated
 * This is the only part that needs changing for different f(x) functions
 */
double function ( double x )
{
    double o;
    o = x*cos(x*x);
    return o;
}

/* Defines the function f(x)*cos(nx) in order to work out 'an'
 * This does not need to be changed for different f(x) functions as it calls up 
 * 'function' to define it
 */
double function_cos ( double x )
{
    double o;
    o = function(x)*cos(nsigma*x);
    return o;
}

/* This is the same as above but for the sin part of the equation 'bn'
 */
double function_sin ( double x )
{
    double o;
    o = function(x)*sin(nsigma*x);
    return o;
}

/* integrate the function f(x) numerically between limits x=a and x=b,
 * using the simpsons rule with n steps. Returns the calculated integral.
 */
double integrate ( double (*f)(double), double a, double b, int n )
{
	int x;
	double d,o;
	/* The first part of the equation, 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 multiplies the whole equation by 'd/3' */
	o = (d/3) * o;
 
	/* Returns the answer */
	return o;
}


/* print function f(x) for equally space values of x between -pi and pi,
 * and also an approximation to f(x) given by its first few Fourier
 * coefficients
 */
void fourier_print ( double (*f)(double),  /* function f(x) */
		     int nmax,             /* number of terms of Fourier series */
		     double *a, double *b, /* Fourier coefficients in arrays */
		     int ndiv )            /* number of equal parts to divide
					      interval -pi to pi into */
{	
	/* Set up the doubles and integers - variables */
	double x,fx;
	int n;
	
	/* Printing out the series in the required format, so it can be used
     * with the plotter.
     * Prints out the values of f(x)
     */
	for (x = -PI; x <= PI; x = x + (PI2 / ndiv))
	printf("%f %f\n",x,f(x));
	printf("&\n");
	
	/* This will work out the fourier series values to print out */
	for (x = -PI; x <= PI; x = x + (PI2 / ndiv))
	{
		fx = 0;
		for (n = 1; n <= nmax; n++)
			fx = fx + (a[0]/2) + ( a[n]*cos(n*x) ) + ( b[n]*sin(n*x) );
		
		printf("%f %f\n",x ,fx);
	}
	
	return;
}

int main ( int argc, char *argv[] )
{
	/* Set up the doubles and integers - variables */
    int nmax,ndiv,idiv;
	double a[max], b[max];
 
    /* Take in the divisions of integration */
    fprintf(stderr, "Please enter the amount of divisons for integration\n");
    scanf("%i", &idiv);
    fflush(stdin);
	
	/* Take in the 'nmax' value which is the elements of the fourier series */
	fprintf(stderr, "Please enter the maximum amount of terms for the fourier series\n");
	scanf("%i", &nmax);
    fflush(stdin);
    	
	/* If 'nmax' is greater than 'max' then abort with a exit error statement */
	if (nmax > max)
	{
		fprintf(stderr, "nmax is too big\n");
		exit(1);
	}
	
	/* Take in the 'ndiv' value which is the amount of divisions between the limits */
	fprintf(stderr, "Please enter the amount of divisions between -pi and +pi\n");
	scanf("%i", &ndiv);
    fflush(stdin);
	
	/* The value for 'a[0]' - from fourier series */
	a[0] = (1/PI)*integrate(function,-PI,PI,idiv);
	
	/* The value for 'a[n]' - the even function from fourier series */
	for(nsigma = 1; nsigma <= nmax; nsigma++)
		a[nsigma] = (1/PI)*integrate(function_cos,-PI,PI,idiv);
		
	/* The value for 'b[0]' (always = 0) from fourier series */
	b[0] = 0;
	
	/* The value for 'b[n]' - the odd function from fourier series */
	for(nsigma = 1; nsigma <= nmax; nsigma++)
		b[nsigma] = (1/PI)*integrate(function_sin,-PI,PI,idiv);
	
	/* call up the fuction 'fourier_print' in order to print out the data */
	fourier_print(function,nmax,a,b,ndiv);
    
    /* Exit the main function with '0' */
	return 0;
}