 # include "data.h"

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

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

  if (argc < 2) {
    fprintf (stderr, "USAGE: datafile (filename)\n");
    return (FALSE);
  }
  
  strcpy (filename, argv[1]);
  if (f != (FILE *) NULL) { fclose (f); }
  f = fopen (filename, "r");
  if (f == (FILE *) NULL) {
    fprintf (stderr, "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;
  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);
  }

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

  if (f == (FILE *) NULL) {
    fprintf (stderr, "no open file for read\n");
    return (FALSE);
  }
  fseek (f, 0, SEEK_SET);

  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) {
      fprintf (stderr, "USAGE: read name N name N ...\n");
      free (vec);
      free (col);
      return (FALSE);    
    }
    colstr = argv[2*i+2];
    for (j = 0; j < strlen (colstr); j++) {
      if (!isdigit(colstr[j])) {
	fprintf (stderr, "USAGE: read name N name N ...\n");
	free (vec);
	free (col);
	return (FALSE);    
      }
    }
    col[i] = atof (colstr);
  }
  
  NELEM = 1000;
  for (i = 0; i < Nvec; i++) {
    REALLOCATE (vec[i][0].elements, float, 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 (&value, col[i], c0);
	  vec[i][0].elements[N] = value;
	  if (!status) vec[i][0].elements[N] = 0.0/0.0;
	}
	if (status) N++;
      }
      c0 = c1 + 1;
      if (N == NELEM) {
	NELEM += 1000;
	for (i = 0; i < Nvec; i++) {
	  REALLOCATE (vec[i][0].elements, float, NELEM);
	}
      }
	
    }
  }
  for (i = 0; i < Nvec; i++) {
    REALLOCATE (vec[i][0].elements, float, MAX (N,1));
    vec[i][0].Nelements = N;
  }
  
  free (vec);
  free (col);
  free (buffer);
  return (TRUE);

}

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

  int i, j, k, Nbytes, Nread, Nextend, Ny, Binary;
  char type[16];
  FTable table;
  Header header;
  Vector **vec;

  if (argc < 2) {
    fprintf (stderr, "USAGE: read -fits extension name name ...\n");
    fprintf (stderr, "  ('extension' may be a name or number)\n");
    return (FALSE);
  }

  if (f == (FILE *) NULL) return (FALSE);
  fseek (f, 0, SEEK_SET);
  table.header = &header;

  /* load appropriate extension (if extname is a number, use count) */
  if (isdigit(extname[0])) {
    Nextend = atoi (extname);

    for (i = 0; i <= Nextend; i++) {
      if (!fits_load_header (f, &header)) return (FALSE);
      Nbytes = fits_matrix_size (&header);
      /* skip the prior data buffers */
      if (i < Nextend) { 
	fseek (f, Nbytes, SEEK_CUR);
	fits_free_header (&header);
	continue;
      }
      /* on the desired header, load buffer */
      ALLOCATE (table.buffer, char, Nbytes);
      Nread = fread (table.buffer, sizeof (char), Nbytes, f);
      if (Nread != Nbytes) {
	fits_free_table  (&table);
	fits_free_header (&header);
	return (FALSE);
      }
      table.size = Nbytes;
    }
  } else {
    if (!fits_fread_ftable (f, &table, extname)) return (FALSE);
  }

  /* identify table type (ascii / binary) */
  Binary = FALSE;
  fits_scan (&header, "XTENSION", "%s", 1, type);
  if (!strcmp (type, "BINTABLE")) {
    Binary = TRUE;
  } 
  Ny = header.Naxis[1];

  /* find columns which match requested vectors */
  for (i = 1; i < argc; i++) {
    void   *data;
    short  *Ps;
    int    *Pi;
    float  *Pf;
    double *Pd;
    int Nval;
    char name[80];
      
    Nval = 0;
    if (Binary) {
      if (!fits_get_bintable_column_type (&header, argv[i], type, &Nval)) return (FALSE);
      if (!strcmp (type, "char")) return (FALSE);
      if (!fits_get_bintable_column (&header, &table, argv[i], &data)) return (FALSE);
    } else {
      if (!fits_get_table_column_type (&header, argv[i], type)) return (FALSE);
      if (!strcmp (type, "char")) return (FALSE);
      if (!fits_get_table_column (&header, &table, argv[i], &data)) return (FALSE);
    }
    if (Nval == 0) return (FALSE);

    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) return (FALSE);    
      REALLOCATE (vec[j][0].elements, float, MAX (Ny,1));
      vec[j][0].Nelements = Ny;
    }

    if (!strcmp (type, "double")) {
      Pd = (double *) data;
      for (j = 0; j < Ny; j++) {
	for (k = 0; k < Nval; k++, Pd++) {
	  vec[k][0].elements[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[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[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[j] = *Ps;
	}
      }
    }
    free (data);
    free (vec);
  }
  fits_free_table (&table);
  fits_free_header (&header);
  return (TRUE);
}
