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

/* A safe getline function. */
int get_line(char** buf, 
	     size_t* bufsize, 
	     FILE* F)
{
  char c;
  size_t j;
  char* tmp;

  /* Only return -1 if first character is EOF. */
  c = getc(F);
  if (c == EOF) return -1;

  /* The current position in the buffer. */
  j = 0;

  while (1)
    {
      if (c == '\n') break;
      
      /* Need to grow the buffer. */
      if (j >= *bufsize-2)
        {
          tmp = (char*)malloc(*bufsize+1000);
          memcpy(tmp, *buf, *bufsize);
          free(*buf);
          *buf = tmp;
          *bufsize += 1000;
        }
      
      (*buf)[j++] = c;

      /* Get the next character. */
      c = getc(F);
      if (c == EOF) break;
    }

  /* Terminate the string. */
  (*buf)[j++] = '\0';
  
  return 0;
}



double* read_table(const char* fname, 
		   int* nrow, 
		   int* ncol) 
{
  FILE* F;
  size_t bufsize;
  char *buf, *tok;
  int xs, xincr, f, d, jj;
  double *x, *x1;

  /* Open the file. */
  F = fopen(fname, "r");
  if (F == 0) 
    {
      printf("Unable to open %s.\n", fname);
      exit(0);
    }

  /* Storage for get_line. */
  bufsize = 1000;
  buf = (char*)malloc(bufsize*sizeof(char));

  /* Size of data array. */
  xs = 1000;
  
  /* Increment size for data array. */
  xincr = 1000;

  /* Array of data values. */
  x = (double*)malloc(xs*sizeof(double));

  /* Current location in x. */
  jj = 0;

  /* Count the number of rows and columns. */
  *nrow = 0;
  *ncol = 0;

  /* Loop over the lines. */
  f = 1;
  while (1) {

    /* Break on end of file or blank line. */
    if (get_line(&buf, &bufsize, F) == -1) break;
    if (buf[0] == '#') continue;
    if (strpbrk(buf, "1234567890") == NULL) break;

    /* Begin splitting off tokens. */
    tok = strtok(buf, " \n\t");
    d = 0;
    ++(*nrow);

    while (1) 
      {
	/* Resize x if we are out of room. */
	if (jj >= xs) 
	  {
	    x1 = (double*)malloc((xs+xincr)*sizeof(double));
	    memcpy(x1, x, xs*sizeof(double));
	    xs += xincr;
	    free(x);
	    x = x1;
	  }
	
	/* Stick the new value on the end. */
	x[jj++] = (double)atof(tok);
	++d;
	
	/* Get the next token. */
	tok = strtok(NULL, " \n\t");
	if (tok == NULL) break;
      }
    
    if (f == 1) 
      {
	*ncol = d;
	f = 0;
      }
    else if (d != *ncol) 
      {
	printf("d=%i\n", d);
	printf("Inconsistent number of elements at line %i of %s\n", 
	       *nrow, fname);
	return NULL;
      }
  }
  
  /* Clean up and return. */
  free(buf);
  return x;
}


int main(int argc, char** argv)
{
  int nrow, ncol;
  double* Z;

  /* This is an example of how to read in a data file. */
  Z = read_table("test.dat", &nrow, &ncol);

  return 1;
}