/*
*           FOURIER SYNTHESIS
*/

#include"lib351.h"

void main(){

	// data
	
	bool goout = false;

    	ofstream dataOutStr;

		const int nRows = 51;
		const int nCols = 65;
		const int nColsZ = 33;
		const double pi = 3.14159265;
		// profile analysis coefficients
		const double at[48] = { .1322, .0138, .0816,-.0360, .0122,
                               -.0379,-.0046,-.0147, .0089, .0036, 
                                .0084, .0050,-.0027,-.0015,-.0055, 
                                .0013,-.0042, .0050,-.0017, .0030,
                               -.0041, .0010,-.0047,-.0002,-.0019, 
                                .0020, .0001, .0017, .0002,-.0011,
                               -.0016,-.0021,-.0001,-.0008, .0023, 
                                .0002, .0015,-.0016, .0003,-.0016, 
                                .0001,-.0002, .0013, .0004, .0007,
			                   -.0003,-.0006,-.0009};
		const double bt[48] = {-.1534,-.0468, .0777,-.0440, .0593,
                               -.0544, .0387,-.0457, .0348,-.0322,
                                .0381,-.0229, .0306,-.0210, .0205,
                               -.0246, .0133,-.0202, .0127,-.0130,
                                .0142,-.0103, .0099,-.0121, .0066,
                               -.0117, .0075,-.0081, .0088,-.0081,
                                .0079,-.0090, .0065,-.0093, .0081,
                               -.0081, .0093,-.0074, .0089,-.0092,
                                .0076,-.0098, .0074,-.0088, .0091,
						       -.0077, .0084,-.0088};
		int nRow, nCol, nWave, nWaves, nType, mWave, mCol;
		double a,b,x,dx,f,fn;
		double ds[101];
		char tic;

    outFileMaker("m352p1c.out", &dataOutStr);
	
	while (!goout){
		
		// data input and restriction

		cout <<"Please supply values:"<< endl;
		cout <<"\nnumber of waves, N ="<< endl;
			cin  >>nWaves;
		cout <<"\nProblem Type (0,1,2,3) ="<< endl;
			cin  >>nType;
		// Type: 0 = square wave; 1 = triangle; 2 = spike; 3 = profile
		// profile is limited to 48 waves or less
		if ((nType == 3) && (nWaves > 48)) {nWaves = 48;}

        // write the header
		dataOutStr<<"\f"<<endl;
		dataOutStr<<"    FOURIER SYNTHESIS"<<endl;
        dataOutStr<<"number of waves = "<<nWaves<<""<<endl;
        dataOutStr<<"problem type = "<<nType<<""<<endl;

		//write horiz line before table
		dataOutStr<<"    X     F(X)   ";
		for (nCol = 1; nCol <= nCols; nCol++) {dataOutStr<<".";}
		dataOutStr<<' '<<endl;

		// initialize the average deviation
		for (nWave = 1; nWave <= 100; nWave++) {
			ds[nWave] = 0.0;
		}

		// loop on x values
		dx = 2*pi/(nRows - 1);
		for (nRow = 1; nRow <= nRows; nRow++){

			x = - pi + (nRow - 1)*dx;
			
			// get f, the actual value of the function at this x

			// one square wave cycle
			if (nType == 0) {
				f = 0;
				if ((nRow > 1) && (nRow < (nRows/2 + 1))) {
					f = -1;
				}
				if ((nRow > (nRows/2 + 1)) && (nRow < nRows)) {
					f = +1;
				}
			}
			
			// one saw-tooth cycle 
			if (nType == 1) {f = 1.0 - 2.0*fabsf(x)/pi;}
				
			// one spike
			if (nType == 2) {
				f = 0.0;
				if (nRow == (nRows/2 + 1)) {
					f = 1.0;
				}
			}

			// the 'profile' is exact at the given x-values if all waves are used
			if (nType == 3) {
				f = 0;
				for (mWave = 1; mWave <= 48; mWave++) {
					f = f + at[mWave-1]*cos(mWave*x) + bt[mWave-1]*sin(mWave*x);
				}
			}

			// get fn, the synthesis value of the function at this x
			fn = 0;
			for (int nWave = 1; nWave <= nWaves; nWave++) {
			
				// square wave; only odd sine waves contribute
				if (nType == 0) {
					a = 0;
					b = 0;
					if (nWave != 2*(nWave/2)) {
						b = 4.0/(pi*nWave);
					}
				}

				// saw-tooth; only odd cosine waves contribute
				if (nType == 1) {
					b = 0;
					a = 0;
					if (nWave != 2*(nWave/2)) {
						a = 8.0/pow((pi*nWave),2);
					}
				}
				
				// spike; only cosine waves contribute
				if (nType == 2) {
					a = 1.0/nWaves;
					b = 0.0;
				}
				
				// profile; get coefficients from analysis tables
				if (nType == 3) {
					a = at[nWave-1];
					b = bt[nWave-1];
				}

				fn = fn + a*cos(nWave*x) + b*sin(nWave*x);
				if (nWaves == 100) {
					ds[nWave]=ds[nWave]+fabs(f-fn)/51.0;
				}
			}
			dataOutStr <<setiosflags(ios::fixed)<<setw(7)<<setprecision(4)<<x<<' '<<setw(7)<<setprecision(4)<<fn<<"  ";
			mCol = int(16.0*(fn+2.06251));
			if (nType == 3) {mCol = int(30.0*fn) + 32;}
			for (int nCol = 1; nCol <= nCols; nCol++) {
				tic = ' ';
				if (nCol == 1) {tic = '.';}
				if (nCol == nCols) {tic = '.';}
				if ((nCol == nColsZ) && (nType != 3)) {tic = '|';}
				if (nCol == mCol) {tic = '*';}
				dataOutStr <<tic;
			}
			dataOutStr <<endl;
		}
		// write horiz line after table
		dataOutStr<<"                 ";
		for (nCol = 1; nCol <= nCols; nCol++) {dataOutStr<<".";}
		dataOutStr<<endl;

		// print the deviation graph
		if (nWaves == 100) {
			// header
			dataOutStr<<"\f"<<endl;
			dataOutStr<<"Average Deviations for N = 2,4,...,100 Waves of Type "<<nType<<endl<<endl;
			dataOutStr<<"      N    DEV"<<endl;
			// horizontal dotted line
			dataOutStr<<"                 ";
			for (nCol = 1; nCol <= nCols; nCol++) {dataOutStr<<".";}
			dataOutStr<<' '<<endl;
			// loop on number of waves
			for (nWave = 2; nWave <= 100; nWave = nWave + 2) {
				dataOutStr<<setiosflags(ios::fixed)<<setw(7)<<nWave<<' '<<setw(7)<<setprecision(4)<<ds[nWave-1]<<"  ";
				//dataOutStr<<setiosflags(ios::fixed)<<setw(7)<<setprecision(4)<<nWave<<' '<<setw(7)<<setprecision(4)<<ds[nWave-1]<<"  ";
				double dev = ds[nWave];
				mCol = int(60.0*ds[nWave]/ds[2]) + 1;
				for (nCol = 1; nCol <= nCols; nCol++) {
					tic = ' ';
					if (nCol == 1) {tic = '.';}
					if (nCol == nCols) {tic = '.';}
					if (nCol == mCol) {tic = '*';}
					dataOutStr <<tic;
				}
				dataOutStr<<endl;
			}
			// horizontal dotted line
			dataOutStr<<"                 ";
			for (nCol = 1; nCol <= nCols; nCol++) {dataOutStr<<".";}
			dataOutStr<<' '<<endl;
		}

		goout = query();
	}
};