#include <stdio.h>
#include <math.h>
#include "comphys.h"
#include "comphys.c"

/* Computational Physics -Posted solutions to Exercise 9.5 */
/* C Thaxton - Appalachian State University, 2007 */
/* Fits a sine series n = 3 to data95.dat */

int main()
{
  float x,y;
  float **A,**b;
  int M;
  int i,j,k,l;
  FILE *in,*out;
  double pi = 3.141592653589793;

  M = 4;

  A = matrix(1,M,1,M);
  b = matrix(1,M,1,1);

  for(j=1;j<=M;j++) {
    b[j][1] = 0.0;
    for(k=1;k<=M;k++) A[j][k] = 0.0;
  }

  if((in = fopen("data95.dat","r")) == NULL) {
    printf("Cannot find data95.dat\n");
    return(1);
  }

  l = 0;
  while(fscanf(in,"%f %f",&x,&y) != EOF) {
    for(k=1;k<=M;k++) {
      if(k == 1) b[k][1] += y;
      else b[k][1] += y*sin((k-1)*pi*x);
      for(j=1;j<=M;j++) {
        if(k == 1 && j == 1) A[k][j] += 1.0;
        else if(k == 1 && j != 1) A[k][j] += sin((j-1)*pi*x);
        else if(k != 1 && j == 1) A[k][j] += sin((k-1)*pi*x); 
        else A[k][j] += sin((j-1)*pi*x)*sin((k-1)*pi*x);  
      }
    }
    l++;
  }
  
  fclose(in);
  
  if((out = fopen("output.dat","w")) == NULL) {
    printf("Cannot find output.dat\n");
    exit(1);
  }
  fprintf(out,"\nThe matrix A is defined as:\n"); 
  for(i=1;i<=M;i++) {
    fprintf(out,"%f %f %f %f\n",A[i][1],A[i][2],A[i][3],A[i][4]);
  }

  gaussj(A,M,b,1);

  for(k=1;k<=M;k++) fprintf(out,"a[%d] = %f +/- %f\n",k,b[k][1],
                    sqrt(A[k][k])/(float)(l-M));
  fclose(out);
  printf("\nGo to file 'output.dat' to see results.\n");
  return(0); 
}