/*
*           LEAST SQUARES FITTING OF EXPERIMENTAL DATA
*
* Peter Signell, 3/02
*/

#include"lib351.h"

void main(){

    // output initialization
	ofstream dataOutStr;
    outFileMaker("m359p1c.out", &dataOutStr);
	ifstream dataInStr;
	inputFileFinder("m359p1c.dat", &dataInStr);

    
    // program variable definitions
	const int wPage = 61;
	const double xScale = 6.;

	// write header
	dataOutStr<<"m359p1c.out: Least Squares Straight Line Fit"<< endl<<endl;

	//loop over runs
	bool goout = false;
	while (!goout){

		// data input
		cout<<"The data will be read from your file."<<endl;
		cout<<"It must be in the following order:"<<endl;
		cout<<"the first line has N, the number of data,"<<endl;
		cout<<"succeeding lines have x(n), y(n), dy(n)."<<endl;

		int N;

		char *imageLine;
		imageLine = new char [wPage];

		readANumber(&N, &dataInStr);

		double *x, *y, *deltaY;
		x = new double [N];
		y = new double [N];
		deltaY = new double [N];

		int setm;
		double sX, sXX, sY, sXY, s1, chi, xMax, xMin, yMax, yMin;
		double sig, a1, a2, x1, y1, residual;
		double da1, da2, rav, deltaX, deltaY1, II;
      	int m1, m2;

		readArray(x, N, &dataInStr);
		readArray(y, N, &dataInStr);
		readArray(deltaY, N, &dataInStr);

		dataOutStr<<"N = "<<N<<endl<<endl;
		dataOutStr<<"    X      Y        DY     RESIDUALS"<<endl;

		if(deltaY[0] > 0){
			sig = pow(deltaY[0], 2);
		}
		else{
			sig = 1;
		}
		s1 = 1./sig;
		sX = x[0]/sig;
		sXX = pow(x[0], 2)/sig;
		sY = y[0]/sig;
		sXY = x[0]*y[0]/sig;
		xMax = x[0];
		xMin = x[0];
  		yMax = y[0];
  		yMin = y[0];

		for(int i = 1; i < N; i++){
			if(deltaY[i] > 0){
				sig = pow(deltaY[i], 2);
			}
			else{
				sig = 1;
			}
			s1 += 1./sig;
			sX +=  x[i]/sig;
			sXX += pow(x[i], 2)/sig;
			sY += y[i]/sig;
			sXY += x[i]*y[i]/sig;

			if(x[i] < xMin){
				xMin = x[i];
			}
			if(x[i] > xMax){
				xMax = x[i];
			}
			if(y[i] < yMin){
				yMin = y[i];
			}
			if(y[i] > yMax){
				yMax = y[i];
			}
		}

		a1 = (sXX * sY - sX * sXY)/(s1 * sXX - sX * sX);
		a2 = (sX * sY - s1 * sXY)/(sX * sX - sXX * s1);

		chi = 0;
		for(int j = 0; j < N; j++){
			y1 = a1 + a2 * x[j];
			residual = y[j] - y1;
			if(deltaY[j] > 0){
				sig = pow(deltaY[j], 2);
			}
			else{
				sig = 1;
			}
			chi += residual*residual/sig;

  	   		dataOutStr<<setiosflags(ios::fixed)<<"  "<<setw(3)<<j+1<<"  "<<setw(8)<<setprecision(5)<<y[j]<<"  "<<setw(8)<<setprecision(5)<<deltaY[j]<<"  "<<setw(8)<<setprecision(5)<<residual<<endl;
		}
  	   	dataOutStr<<endl;

		// plot the graph
		da1 = sqrt(sXX/(s1 * sXX - sX * sX));
		da2 = sqrt(s1/(s1 * sXX - sX * sX));

		if(deltaY[0] == 0){
			rav = sqrt(chi)/N;
			da1 *= rav;
			da2 *= rav;
 		}

		deltaX = 1.2 * (xMax - xMin)/50;
		deltaY1 = 1.2 * (yMax - yMin)/(wPage - 1);

		II = int(a2 * deltaX/deltaY1 + 0.5);
		deltaX = II * deltaY1/a2;

		x1 = xMin - 5 * deltaX;

		dataOutStr<<"a1 = "<<a1<<" +/- "<<da1<<endl;
		dataOutStr<<"a2 = "<<a2<<" +/- "<<da2<<endl;
		dataOutStr<<"chisq = "<<chi<<endl;
	   	dataOutStr<<endl;
		//lineOutput(61, '_', &dataOutStr);

		//lineOutput(61, '_', &dataOutStr);
   		dataOutStr<<"\f"<<endl;
   		dataOutStr<<"                ";
		lineOutput(58, '_', &dataOutStr);
        dataOutStr<<endl;

		for(int m = 0; m < wPage; m++){
			imageLine[m] = ' ';
		}

		for(i = 0; i < 51; i++){
			imageLine[0] = '|';
         	imageLine[wPage -2] = '|';
			imageLine[wPage -1] = '\n';

	  	y1 = a1 + a2*x1;
	  	m1 = int((y1 - yMin)/deltaY1 + xScale + 0.5);

	  	if(m1 >= 1 && m1 < wPage){
				imageLine[m1] = '*';
                }
    		setm = 0;

			for(int j = 0; j < N; j++){
				if(fabs(x1 - x[j]) <= deltaX/2){
					m2 = int((y[j] - yMin)/deltaY1 + xScale + 0.5);
					imageLine[m2] = '0';
		  		setm = 1;
				}
			}
			dataOutStr<<setiosflags(ios::fixed)<<setw(7)<<setprecision(3)<<x1<<setw(7)<<setprecision(3)<<y1<<" ";
	 		for(int m = 0; m < wPage; m++){
				dataOutStr<<imageLine[m];
			}

	 		if(m1 >= 1 && m1 < wPage){
				imageLine[m1] = ' ';
	 		}
	 		if(setm == 1){
				imageLine[m2] = ' ';
	 		}
	 		x1 += deltaX;
		}
   		dataOutStr<<"                ";
		lineOutput(58, '~', &dataOutStr);
        dataOutStr<<endl;

		// another run?
		goout = query();
   }
    // close the output file
	fileClose("m359p1c.out", &dataOutStr);
};