#include <stdio.h>
#include <math.h>

/* Computational Physics - posted solution to Exercise 9.1 */
/* C Thaxton - Appalachian State University, 2007 */
/* Determines least square fit y = a + bx to data in data61.dat */

int main()
{
  float x,y;
  float S,Sx,Sxx,Sy,Sxy,D;
  float a,b,siga,sigb;
  int k;
  FILE *in;

  if((in = fopen("data91.dat","r")) == NULL) {
    printf("\nSilly Rabbit, data61.dat does not exist!!\n");
    return(1);
  }

  S = Sx = Sxx = Sy = Sxy = 0.0;

  /* Note: data61.dat does not contain errors for each point, so we
     simply assume that the sigma for each point = 1 */
  /* Read in the data and form the sums */

  k = 0;
  while(fscanf(in,"%f %f",&x,&y) != EOF) {
    S += 1.0;
    Sx += x;
    Sy += y;
    Sxx += x*x;
    Sxy += x*y;
    k++;
  }

  /* Calculate the determinant */

  D = S*Sxx - Sx*Sx;

  a = (Sxx*Sy - Sx*Sxy)/D;
  b = (S*Sxy - Sx*Sy)/D;

  siga = sqrt(Sxx/D)/(float)(k-2);
  sigb = sqrt(S/D)/(float)(k-2);

  printf("\na = %f +/- %f",a,siga);
  printf("\nb = %f +/- %f\n",b,sigb);

  return(1);
}