/*  
    FITLINE_SIMPLE.cxx

    Fit the straight line to y = ax + b, to find the best-fit values for a and b. The errors for a and b are also calculated. The algorithm is the simplest linear fitting program, without considering the errors in x and y.

    The algorithm is from Press et al. (1992), "Numerical Recipes", Pg 661-665.

    Make input file to be: input.dat, which contains x and y values.

    Compile: g++ -o fitlines fitline_simple.cxx
*/

#include <iostream.h>
#include <fstream.h>
#include <math.h>

main()
{
  int N;
  double x, y;

  /* Calculate the parameters */
  double S, Sx, Sy;
  double Sxx, Sxy, delta;
  double a, b;
  double siga, sigb;
  double chi2, sigma;

  ifstream in("input.dat");
  
  S = 0.0;
  Sx = 0.0;
  Sy = 0.0;
  Sxx = 0.0;
  Sxy = 0.0;
  
  int count = 0;

  while(in>>x>>y)
    {
      Sx += x;
      Sy += y;
      Sxx += x * x;
      Sxy += x * y;
      S += 1;
      count++;
    }

  N = count;
  delta = S * Sxx - Sx * Sx;
 
  a = (Sxx * Sy - Sx * Sxy) / delta;
  b = (S * Sxy - Sx * Sy) / delta;

  siga = sqrt(Sxx/delta);
  sigb = sqrt(S/delta);

  in.close();

  /* Calculate chi2 */
  ifstream ina("input.dat");

  chi2 = 0.0;
  for(int i=1; i<=N; i++) {    
    ina>>x>>y;

    chi2 += (y - a - b * x) * (y - a - b * x);
  }

  siga = siga * sqrt(chi2/(float(N-2)));
  sigb = sigb * sqrt(chi2/(float(N-2)));
  
  sigma = sqrt(chi2/(float(N-1)));

  ina.close();

  /* Output results */
  cout<<"  ***** Fitting Data to Straight Line (Simplest Case) *****\n\n";
  cout<<"Number of data = "<<N<<"\n";
  cout<<"Slope          = "<<b<<" +- "<<sigb<<"\n";
  cout<<"Zero-point     = "<<a<<" +- "<<siga<<"\n";
  cout<<"Chi2           = "<<chi2<<"\n";
  cout<<"Sigma          = "<<sigma<<"\n";
  
}