# include "data.h"

FILE *f = (FILE *) NULL;
char filename[2048];

int datafile (int argc, char **argv) {

  if (argc < 2) {
    gprint (GP_ERR, "USAGE: datafile (filename)\n");
    return (FALSE);
  }
  
  strcpy (filename, argv[1]);
  if (f != (FILE *) NULL) { fclose (f); }
  f = fopen (filename, "r");
  if (f == (FILE *) NULL) {
    gprint (GP_ERR, "failed to open file %s\n", argv[1]);
    return (FALSE);
  }
  return (TRUE);
}

int read_vectors (int argc, char **argv) {
  
    int i, j, Nskip, Nvec, *col, done, status, IsCSV, VERBOSE;
  int Nbytes, nbytes, Nstart, NELEM, N, nread;
  char *colstr, *c0, *c1, *buffer, *extname;
  double value;
  Vector **vec;

  /* auto-sense table type */
  if ((N = get_argument (argc, argv, "-fits"))) {
    remove_argument (N, &argc, argv);
    extname = strcreate (argv[N]);
    if (extname == (char *) NULL) return (FALSE);
    remove_argument (N, &argc, argv);
    status = read_table_vectors (argc, argv, extname);
    free (extname);
    return (status);
  }

  Nskip = 0;
  if ((N = get_argument (argc, argv, "-skip"))) {
    remove_argument (N, &argc, argv);
    Nskip = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

  IsCSV = FALSE;
  if ((N = get_argument (argc, argv, "-csv"))) {
    remove_argument (N, &argc, argv);
    IsCSV = TRUE;
  }

  VERBOSE = FALSE;
  if ((N = get_argument (argc, argv, "-v"))) {
    remove_argument (N, &argc, argv);
    VERBOSE = TRUE;
  }

  if ((argc < 3) || !(argc % 2)) {
    gprint (GP_ERR, "USAGE: read name N name N ...\n");
    return (FALSE);
  }
  /* read name N name N  */

  if (f == (FILE *) NULL) {
    gprint (GP_ERR, "no open file for read\n");
    return (FALSE);
  }
  fseeko (f, 0LL, SEEK_SET);

//  if (IsCSV) {
//    status = read_vectors_csv (argc, argv, Nskip, f);
//    return status;
//  }

  Nvec = (argc - 1) / 2;
  ALLOCATE (vec, Vector *, Nvec);
  ALLOCATE (col, int, Nvec);

  for (i = 0; i < Nvec; i++) {
    if ((vec[i] = SelectVector (argv[2*i + 1], ANYVECTOR, TRUE)) == NULL) {
      gprint (GP_ERR, "USAGE: read name N name N ...\n");
      free (vec);
      free (col);
      return (FALSE);    
    }
    // XXX we could allow flags (eg, N.i) to specify INT vs FLT types vectors...
    colstr = argv[2*i+2];
    for (j = 0; j < strlen (colstr); j++) {
      if (isdigit(colstr[j])) {
	col[i] = atof (colstr);
	continue;
      }
      // allow 'excel' columns names of the form A-Z, AA-AZ, BA-BZ, .. ZA-ZZ
      if (strlen(colstr) >= 3) goto bad_colname;
      if (colstr[0] < 'A') goto bad_colname;
      if (colstr[0] > 'Z') goto bad_colname;
      if (colstr[1] && colstr[1] < 'A') goto bad_colname;
      if (colstr[1] && colstr[1] > 'Z') goto bad_colname;

      col[i] = colstr[0] - 'A' + 1;
      if (colstr[1]) {
	col[i] *= 26;
	col[i] += colstr[1] - 'A' + 1;
      }
      continue;

    bad_colname:
      gprint (GP_ERR, "USAGE: read name N name N ...\n");
      free (vec);
      free (col);
      return (FALSE);    
    }
  }

  // currently, all read vectors are forced to be type FLT
  NELEM = 1000;
  for (i = 0; i < Nvec; i++) {
    ResetVector (vec[i], OPIHI_FLT, NELEM);
  }
  
  ALLOCATE (buffer, char, 0x10001);
  bzero (buffer, 0x10001);
  for (i = 0; i < Nskip; i++) {
    scan_line (f, buffer);
  }

  Nstart = 0;
  N = 0;
  done = FALSE;
  while (!done) {
    Nbytes = 0x10000 - Nstart;
    bzero (&buffer[Nstart], Nbytes);
    nread = fread (&buffer[Nstart], 1, Nbytes, f);
    if (ferror (f)) {
      perror ("error reading data file");
      break;
    }
    if (nread == 0) break;
    nbytes = nread + Nstart;
    
    status = TRUE;
    c0 = buffer; 
    while (status) {
      c1 = strchr (c0, '\n');
      if (c1 == (char *) NULL) {
	Nstart = strlen (c0);
	memmove (buffer, c0, Nstart);
	status = FALSE;
      } else {
	*c1 = 0;
      }      
      if ((*c0 != '#') && (*c0 != '!')) {
	for (i = 0; (i < Nvec) && status; i++) {
	  if (IsCSV) {
	    status = dparse_csv (&value, col[i], c0);
	  } else {
	    status = dparse (&value, col[i], c0);
	  }
	  if (status) {
	      vec[i][0].elements.Flt[N] = value;
	  } else {
	      vec[i][0].elements.Flt[N] = NAN;
	      if (VERBOSE) {
		  if (IsCSV) {
		      gprint (GP_ERR, "suspect field: %d (%s) in %s\n", col[i], argv[2*i+2], c0);
		  } else {
		      gprint (GP_ERR, "suspect field: %d in %s\n", col[i], c0);
		  }
	      }
	  }
	}
	if (status) {
	    N++;
	} else {
	    if (VERBOSE && FALSE) {
		char temp[32];
		strncpy (temp, c0, 32);
		temp[31] = 0;
		gprint (GP_ERR, "skip line %s\n\n", temp);
	    }
	}
      }
      c0 = c1 + 1;
      if (N == NELEM) {
	NELEM += 1000;
	for (i = 0; i < Nvec; i++) {
	  REALLOCATE (vec[i][0].elements.Flt, opihi_flt, NELEM);
	}
      }
    }
  }
  for (i = 0; i < Nvec; i++) {
    REALLOCATE (vec[i][0].elements.Flt, opihi_flt, MAX (N,1));
    vec[i][0].Nelements = N;
  }
  
  free (vec);
  free (col);
  free (buffer);
  return (TRUE);

}

# define ESCAPE(MSG) { \
  gprint (GP_ERR, "%s\n", MSG); \
  if (CCDKeyword != NULL) free (CCDKeyword); \
  gfits_free_table  (&table); \
  gfits_free_header (&header); \
  return (FALSE);  }

int read_table_vectors (int argc, char **argv, char *extname) {

  off_t Nbytes;
  int i, j, k, N, Nextend, Ny, Binary, vecType, padIfShort;
  char type[16], ID[80], *CCDKeyword;
  FTable table;
  Header header;
  Vector **vec;

  table.buffer = NULL;
  header.buffer = NULL;

  CCDKeyword = NULL;
  if ((N = get_argument (argc, argv, "-keyword"))) {
    remove_argument (N, &argc, argv);
    CCDKeyword = strcreate (argv[N]);
    remove_argument (N, &argc, argv);
  }

  padIfShort = FALSE;
  if ((N = get_argument (argc, argv, "-pad-if-short"))) {
    remove_argument (N, &argc, argv);
    padIfShort = TRUE;
  }

  Nextend = -1;
  if ((N = get_argument (argc, argv, "-extnum"))) {
    remove_argument (N, &argc, argv);
    Nextend = atoi (extname);
  }

  if (argc < 2) ESCAPE ("USAGE: read -fits extension [-extnum] [-keyword key] name name ...");

  if (f == NULL) ESCAPE ("file not found");
  fseeko (f, 0LL, SEEK_SET);
  table.header = &header;

  /* load appropriate extension (if extname is a number, use count) */
  if (Nextend > -1) {
    // first extension is PHU, cannot be a table. 
    // Nextend counts from 0 for first extension
    if (!gfits_load_header (f, &header)) ESCAPE ("error reading primary header for file");
    Nbytes = gfits_data_size (&header);
    fseeko (f, Nbytes, SEEK_CUR);
    gfits_free_header (&header);

    for (i = 0; i < Nextend; i++) {
      if (!gfits_load_header (f, &header)) ESCAPE ("extension not found");
      Nbytes = gfits_data_size (&header);
      /* skip the prior data buffers */
      fseeko (f, Nbytes, SEEK_CUR);
      gfits_free_header (&header);
    }
    if (!gfits_load_header (f, &header)) ESCAPE ("error reading header for extension");
    if (!gfits_fread_ftable_data (f, &table, padIfShort)) ESCAPE ("error reading table for extension");

  } else {
    if (CCDKeyword == NULL) {
      CCDKeyword = get_variable ("CCDKEYWORD");
    }
    if (CCDKeyword == NULL) {
      CCDKeyword = strcreate ("EXTNAME");
    }

    while (1) {
      if (!gfits_load_header (f, &header)) {
	gprint (GP_ERR, "extension %s not found in file\n", extname);
	if (CCDKeyword != NULL) free (CCDKeyword); 
	gfits_free_table  (&table); 
	gfits_free_header (&header); 
	return (TRUE);  
      }

      Nbytes = gfits_data_size (&header);

      if (!gfits_scan (&header, CCDKeyword, "%s", 1, ID)) {
	fseeko (f, Nbytes, SEEK_CUR);
	gfits_free_header (&header);
	continue;
      }
      if (strcmp (ID, extname)) {
	fseeko (f, Nbytes, SEEK_CUR);
	gfits_free_header (&header);
	continue;
      }
      if (!gfits_fread_ftable_data (f, &table, padIfShort)) ESCAPE ("error reading table for extension");
      break;
    }
  }

  /* identify table type (ascii / binary) */
  Binary = FALSE;
  gfits_scan (&header, "XTENSION", "%s", 1, type);
  if (strcmp (type, "BINTABLE") && strcmp (type, "TABLE")) {
    ESCAPE ("specified extension is not a table\n");
  }
  Binary = !strcmp (type, "BINTABLE");
  Ny = header.Naxis[1];

  /* find columns which match requested vectors */
  for (i = 1; i < argc; i++) {
    void   *data;
    char   *Pc;
    short  *Ps;
    int    *Pi;
    float  *Pf;
    double *Pd;
    int Nval;
    char name[80];
      
    Nval = 0;
    if (Binary) {
      if (!gfits_get_bintable_column_type (&header, argv[i], type, &Nval)) ESCAPE ("requested field not found");
      if (!gfits_get_bintable_column (&header, &table, argv[i], &data)) ESCAPE ("error reading data from specified field");
    } else {
      if (!gfits_get_table_column_type (&header, argv[i], type)) ESCAPE ("requested field not found");
      if (!gfits_get_table_column (&header, &table, argv[i], &data)) ESCAPE ("error reading data from specified field");
    }
    if (Nval == 0) ESCAPE ("no data for field in table");
    
    vecType = OPIHI_INT;
    if (!strcmp (type, "double") || !strcmp (type, "float")) {
      vecType = OPIHI_FLT;
    }
	
    // define the multifield vector names
    ALLOCATE (vec, Vector *, Nval);
    for (j = 0; j < Nval; j++) {
      if (Nval == 1) 
	sprintf (name, "%s", argv[i]);
      else
	sprintf (name, "%s:%d", argv[i], j);
      if ((vec[j] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) ESCAPE ("bad vector name");
      ResetVector (vec[j], vecType, MAX (Ny,1));
    }

    if (!strcmp (type, "double")) {
      Pd = (double *) data;
      for (j = 0; j < Ny; j++) {
	for (k = 0; k < Nval; k++, Pd++) {
	  vec[k][0].elements.Flt[j] = *Pd;
	}
      }
    }
    if (!strcmp (type, "float")) {
      Pf = (float *) data;
      for (j = 0; j < Ny; j++) {
	for (k = 0; k < Nval; k++, Pf++) {
	  vec[k][0].elements.Flt[j] = *Pf;
	}
      }
    }
    if (!strcmp (type, "int")) {
      Pi = (int *) data;
      for (j = 0; j < Ny; j++) {
	for (k = 0; k < Nval; k++, Pi++) {
	  vec[k][0].elements.Int[j] = *Pi;
	}
      }
    }
    if (!strcmp (type, "short")) {
      Ps = (short *) data;
      for (j = 0; j < Ny; j++) {
	for (k = 0; k < Nval; k++, Ps++) {
	  vec[k][0].elements.Int[j] = *Ps;
	}
      }
    }
    if (!strcmp (type, "char")) {
      Pc = (char *) data;
      for (j = 0; j < Ny; j++) {
	for (k = 0; k < Nval; k++, Pc++) {
	  vec[k][0].elements.Int[j] = *Pc;
	}
      }
    }
    free (data);
    free (vec);
  }
  if (CCDKeyword != NULL) free (CCDKeyword);
  gfits_free_table (&table);
  gfits_free_header (&header);
  return (TRUE);
}

# if (0)
int read_vectors_csv (int argc, char **argv, int Nskip, FILE *f) {
  
  int i, j, Nvec, *col, done, status;
  int Nbytes, nbytes, Nstart, NELEM, N, nread;
  char *colstr, *c0, *c1, *buffer, *extname;
  double value;
  Vector **vec;

  Nvec = (argc - 1) / 2;
  ALLOCATE (vec, Vector *, Nvec);
  ALLOCATE (col, int, Nvec);

  for (i = 0; i < Nvec; i++) {
    if ((vec[i] = SelectVector (argv[2*i + 1], ANYVECTOR, TRUE)) == NULL) {
      gprint (GP_ERR, "USAGE: read name N name N ...\n");
      free (vec);
      free (col);
      return (FALSE);    
    }
    // XXX we could allow flags (eg, N.i) to specify INT vs FLT types vectors...
    colstr = argv[2*i+2];
    for (j = 0; j < strlen (colstr); j++) {
      if (!isdigit(colstr[j])) {
	gprint (GP_ERR, "USAGE: read name N name N ...\n");
	free (vec);
	free (col);
	return (FALSE);    
      }
    }
    col[i] = atof (colstr);
  }

  // currently, all read vectors are forced to be type FLT
  NELEM = 1000;
  for (i = 0; i < Nvec; i++) {
    ResetVector (vec[i], OPIHI_FLT, NELEM);
  }
  
  ALLOCATE (buffer, char, 0x10001);
  bzero (buffer, 0x10001);
  for (i = 0; i < Nskip; i++) {
    scan_line (f, buffer);
  }

  Nstart = 0;
  N = 0;
  done = FALSE;
  while (!done) {
    Nbytes = 0x10000 - Nstart;
    bzero (&buffer[Nstart], Nbytes);
    nread = fread (&buffer[Nstart], 1, Nbytes, f);
    if (ferror (f)) {
      perror ("error reading data file");
      break;
    }
    if (nread == 0) break;
    nbytes = nread + Nstart;
    
    status = TRUE;
    c0 = buffer; 
    while (status) {
      c1 = strchr (c0, '\n');
      if (c1 == (char *) NULL) {
	Nstart = strlen (c0);
	memmove (buffer, c0, Nstart);
	status = FALSE;
      } else {
	*c1 = 0;
      }      
      if ((*c0 != '#') && (*c0 != '!')) {
	for (i = 0; (i < Nvec) && status; i++) {
	  status = dparse_csv (&value, col[i], c0);
	  vec[i][0].elements.Flt[N] = value;
	  if (!status) vec[i][0].elements.Flt[N] = NAN;
	}
	if (status) N++;
      }
      c0 = c1 + 1;
      if (N == NELEM) {
	NELEM += 1000;
	for (i = 0; i < Nvec; i++) {
	  REALLOCATE (vec[i][0].elements.Flt, opihi_flt, NELEM);
	}
      }
	
    }
  }
  for (i = 0; i < Nvec; i++) {
    REALLOCATE (vec[i][0].elements.Flt, opihi_flt, MAX (N,1));
    vec[i][0].Nelements = N;
  }
  
  free (vec);
  free (col);
  free (buffer);
  return (TRUE);

}

# endif
