#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
#define max 100

/* The f(x) function */
double fout ( double x )
{
	/* Sets up the square wave as stated */
	if ( -PI <= x && x <= -PI_2)
		return -1;
	if ( -PI_2 < x && x <= PI_2)
		return 1;
	if ( PI_2 < x && x <= PI)
		return -1;
	return 0;
}

/* 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 double and integers - variables */
	double x,fx;
	int n;
	
	/* Printing out the series in the required format */
	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 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[] )
{
	int nmax,ndiv,n;
	double a[max], b[max];
	
	/* Take in the nmax value */
	fprintf(stderr, "Please enter the maximum amount of terms for the fourier series\n");
	scanf("%d", &nmax);
	
	/* If nmax is greater than max then abort */
	if (nmax > max)
	{
		fprintf(stderr, "nmax is too big\n");
		exit(1);
	}
	
	/* Take in the ndiv value */
	fprintf(stderr, "Please enter the amount of divisions between -pi and +pi\n");
	scanf("%d", &ndiv);
	
	/* The value for a[0] */
	a[0] = 0;
	
	/* The value for a[n] - the even function */
	for(n = 1; n <= nmax; n++)
		a[n] = (4/(n*PI)) * sin((n*PI)/2);
		
	/* The value for b[0] (always = 0) */
	b[0] = 0;
	
	/* The value for b[n] - the odd function */
	for(n = 1; n <= nmax; n++)
		b[n] = 0;
	
	/* call up the fuction fourier_print */
	fourier_print(fout,nmax,a,b,ndiv);
	system("PAUSE");
	return 0;
}