#include <math.h>

void fit_line(
         float x[], 
         float y[], 
         int ndata, 
         float sig[], 
         float &a, 
         float &b, 
         float &siga, 
         float &sigb, 
         float &chi2, 
         float &q)
/*
 * From Numerical Recipes in C 
 * Given a set of data points x[1..ndata],y[1..ndata] with individual standard deviations
sig[1..ndata], fit them to a straight line y = a + bx by minimizing χ2. Returned are
a,b and their respective probable uncertainties siga and sigb, the chi-square chi2, and the
goodness-of-fit probability q (that the fit would have chi2 this large or larger). 
*/

{
//  float gammq(float a, float x);
  int i;
  float wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;
  b=0.0;
  a=0.0;

  ss=0.0;
  for (i=0;i<ndata;i++) { 
    wt=1.0/pow(sig[i],2);
    ss += wt;
    sx += x[i]*wt;
    sy += y[i]*wt;
  }

  sxoss=sx/ss;
  for (i=0;i<ndata;i++) {
    t=(x[i]-sxoss)/sig[i];
    st2 += t*t;
    b += t*y[i]/sig[i];
  }

  b /= st2; 
  a=(sy-sx*b)/ss;
  siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
  sigb=sqrt(1.0/st2);
  chi2=0.0; 
  q=1.0;
  
  for (i=0;i<ndata;i++) chi2 += pow(((y[i]-a-b*x[i])/sig[i]),2);
  //if (ndata>2) q=gammq(0.5*(ndata-2),0.5*chi2);

}