#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_fit.h>

/* Maximum name size + 1. */

#define MAX_NAME_SZ 256
#define MAX_LINE_SIZE 256

int file_lines (FILE *file);

int file_read (FILE *file, int n);

int main (int argc, char *argv[]) {

  int lines, i;
  FILE *input;
  char *name, line[MAX_LINE_SIZE];
  double **data, **bet;
  double *x, *y; // PA

  double a, b, cov00, cov01, cov11, chisq; // Fit variables
  int n, status;

  if (!(name = malloc(sizeof(char)*MAX_NAME_SZ))) {
    printf("No memory!\n");
    return 1;
  }

  if (argc == 1) {      // No extra arguments on the command line
    printf("Name of the input file: ");
    fgets(name, MAX_NAME_SZ, stdin);
  } else if (argc == 2) {
    strcpy(name, argv[1]);
  }
  else {
    printf ("Too many arguments! Exiting.\n");
    exit(1);
  }
  if ((strlen(name) > 0) && (name[strlen (name) - 1] == '\n'))
    name[strlen (name) - 1] = '\0';
  
  //printf ("Loading file: %s\n", name);

  if (!(input=fopen(name, "r")))
    exit(1);

  lines = file_lines (input);
  //  printf ("Lines in the file: %d\n", lines);

  if(!(data = malloc(2 * sizeof(double *))))
    exit(1);
  for (i=0; i<2; i++)
    if(!(data[i] = malloc((lines) * sizeof(double))))
      exit(1);
  
  for (i=0; i<lines; i++) {
    fgets(line, MAX_LINE_SIZE, input);
    line[strlen (line) - 1] = '\0';
    data[0][i] = atof(strtok(line, " \t"));
    data[1][i] = atof(strtok(NULL, " \t"));
  }

  // File loaded into data[][]. 


  if(!(bet = malloc(2 * sizeof(double *))))
    exit(1);
  for (i=0; i<=2; i++)
    if(!(bet[i] = malloc((lines) * sizeof(double))))
      exit(1);

  // PA: allocate x and y
  x = malloc(lines * sizeof(double));
  y = malloc(lines * sizeof(double));

  n = 0;
  for (i=0; i<lines; i++) {
    if (data[0][i] >= 0.05 && data[0][i] <= 0.3) {
      // PA: set x and y vectors
      x[n] = data[0][i];
      y[n] = 1/(data[1][i]*((1/data[0][i])-1));

      bet[0][i] = data[0][i];
      bet[1][i] = 1/(data[1][i]*((1/data[0][i])-1));     // y' = 1/(y*(x-1))
      bet[2][i] = 1/(bet[1][i]*0.05);              // weight for the fits, based on 5% estimated error
      printf ("%.8e\t%.8e\t%.8e\n", bet[0][i], bet[1][i], bet[2][i]);
      n++;
    }
  }

  // Linear fit to BET-transformed data from 0.05 < x < 0.3 (y' = a + bx')
#if 1
  {
    printf ("Performing fit with %d points...", n);
    status = gsl_fit_linear (x, 1, y, 1, n, &a, &b, &cov00, &cov01, &cov11, &chisq);
    if (status) printf ("not ");
    printf ("done! (status = %d)\n", status);
  }
#else
  printf ("Performing fit with %d points...", n);
  status = gsl_fit_linear (bet[0], 1, bet[1], 1, n, &a, &b, &cov00, &cov01, &cov11, &chisq);
  if (status) printf ("not ");
  printf ("done! (status = %d)\n", status);
#endif

  printf ("====== BET results ======\n");
  printf ("y = %lf + %lf x\n", a, b);
  printf ("[%g, %g\n %g, %g]\n", 
          cov00, cov01, cov01, cov11);
  printf ("Chi Square = %lf\n", chisq);
  printf ("=========================\n");
  
  fclose (input);
  free (name);
  free (data);
  free (bet);
  return 0;
  
}


int file_lines (FILE *file) {  // Counts the lines in a text file

  int i=0, c;

  while (!feof(file)) {
    c = fgetc(file);
    if (c == '\n')
      i++;
  }

  rewind(file);
  return i;
  
}