# include "dvo.h"
# define D_NSTARS 1000
# define BYTES_STAR 66
# define BLOCK 1000

int cmpread (int argc, char **argv) {
  
  int i, nstar, Nin, Nextra;
  int doneread, done, Nskip, Nbytes, nbytes, Ninstar, RAnDEC;
  char *c, *c2, *name, *file;
  float *R, *D, *V;
  double tR, tD, tmp, X, Y;
  FILE *f;
  char *buffer;
  int Nfield, Nra, Ndec;
  Vector *rvec, *dvec, *vec;
  Header header;
  Coords coords;

  if ((argc != 4) && (argc != 6)) {
    fprintf (stderr, "USAGE: cmpread name Nfield <filename>\n");
    fprintf (stderr, "USAGE: cmpread ra N dec N <filename>\n");
    return (FALSE);
  }

  if (argc == 4) {
    name = strcreate (argv[1]);
    Nfield = atof (argv[2]);
    file = strcreate (argv[3]);
    RAnDEC = FALSE;
  } else {
    Nra  = atof (argv[2]);
    Ndec = atof (argv[4]);
    file = strcreate (argv[5]);
    RAnDEC = TRUE;
  }    

  if (!fits_read_header (file, &header)) {
    fprintf (stderr, "ERROR: can't read header for %s\n", file);
    free (file);
    if (!RAnDEC) free (name);
    return (FALSE);
  }

  /* find expected number of stars */
  fits_scan (&header, "NSTARS", "%d", 1, &nstar);
  if (nstar == 0) {
    fprintf (stderr, "ERROR: can't get NSTARS from header\n");
    free (file);
    if (!RAnDEC) free (name);
    fits_free_header (&header);
    return (FALSE);
  }

  f = fopen (file, "r");
  if (f == NULL) {
    fprintf (stderr, "ERROR: can't read data from %s\n", file);
    free (file);
    if (!RAnDEC) free (name);
    fits_free_header (&header);
    return (FALSE);
  }
  fseek (f, header.size, SEEK_SET); 

  if (RAnDEC) {
    if ((rvec = SelectVector ("ra",  ANYVECTOR, TRUE)) == NULL) return (FALSE);
    if ((dvec = SelectVector ("dec",  ANYVECTOR, TRUE)) == NULL) return (FALSE);
    REALLOCATE (rvec[0].elements, float, nstar);
    REALLOCATE (dvec[0].elements, float, nstar);
    R = rvec[0].elements;
    D = dvec[0].elements;
    rvec[0].Nelements = dvec[0].Nelements = nstar;
    if (!GetCoords (&coords, &header)) {
      fprintf (stderr, "can't get WCS info from header\n");
      return (FALSE);
    }
  } else {
    if ((vec = SelectVector (name, ANYVECTOR, TRUE)) == NULL) return (FALSE);
    REALLOCATE (vec[0].elements, float, nstar);
    V = vec[0].elements;
    vec[0].Nelements = nstar;
  }

  /* load in stars by blocks of 1000 */
  Nin = 0;
  ALLOCATE (buffer, char, (BLOCK*BYTES_STAR) + 1);
  buffer[BLOCK*BYTES_STAR] = 0;
  Nextra = 0;
  doneread = FALSE;
  while (!doneread) {
    Nbytes = BYTES_STAR * BLOCK - Nextra;
    nbytes = fread (&buffer[Nextra], 1, Nbytes, f);
    if (nbytes == 0) {
      doneread = TRUE;
      continue;
    }
    nbytes += Nextra;
    /* check line-by-line integrity */
    c = buffer;
    done = FALSE;
    while ((c < buffer + nbytes) && (!done)) { 
      for (c2 = c; *c2 == '\n'; c2++);
      if (c2 > c) { /* extra return chars */
	memmove (c, c2, (int)(buffer + nbytes - c2));
	Nskip = c2 - c;
	nbytes -= Nskip;
	bzero (buffer + nbytes, Nskip);
      }
      c2 = strchr (c, '\n');
      if (c2 == (char *) NULL) {
	done = TRUE;	
	continue;
      }
      c2++;
      if ((c2 - c) != BYTES_STAR) { /* bad line, delete it */
	memmove (c, c2, (int)(buffer + nbytes - c2));
	Nskip = c2 - c;
	nbytes -= Nskip;
	bzero (buffer + nbytes, Nskip);
      } else {
	c = c2;
      }
    }
    Ninstar = nbytes / BYTES_STAR;
    Nextra = nbytes % BYTES_STAR;
    for (i = 0; i < Ninstar; i++, Nin++) {
      if (RAnDEC) {
	dparse (&X, Nra,  &buffer[i*BYTES_STAR]);
	dparse (&Y, Ndec, &buffer[i*BYTES_STAR]);
	XY_to_RD (&tR, &tD, X, Y, &coords);
	R[Nin] = tR; D[Nin] = tD;
      } else {
	dparse (&tmp, Nfield, &buffer[i*BYTES_STAR]);
	V[Nin] = tmp;
      }
    }
  }
  fclose (f);

  fits_free_header (&header);
  free (buffer);
  free (file);
  if (!RAnDEC) free (name);

  fprintf (stderr, "loaded %d objects\n", Nin);
  return (TRUE);
}

