# include "dvo.h"
# define BYTES_STAR 23
# define BLOCK 1000
# define DNSTARS 1000

int cmatch (int argc, char **argv) {
  
  FILE *f;
  Catalog catalog1, catalog2;
  int i, Nbytes, nbytes, Nitems, nitems;
  char *tbuffer, filename[128];
  int nstar, NSTARS;
  double R, D, M, radius;
  char catdir[256], gscdir[256];
  Vector *rvec, *dvec, *mvec, *drvec, *ddvec, *dmvec;

  VarConfig ("CATDIR", "%s", catdir);
  VarConfig ("GSCDIR", "%s", gscdir);
      
  if (argc != 9) {
    fprintf (stderr, "USAGE: cmatch file radius (RA) (DEC) (Mag) (dRA) (dDEC) (dMag)\n");
    return (FALSE);
  }

  radius = atof (argv[2]);

  if ((rvec  = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
  if ((dvec  = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
  if ((mvec  = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
  if ((drvec = SelectVector (argv[6], ANYVECTOR, TRUE)) == NULL) return (FALSE);
  if ((ddvec = SelectVector (argv[7], ANYVECTOR, TRUE)) == NULL) return (FALSE);
  if ((dmvec = SelectVector (argv[8], ANYVECTOR, TRUE)) == NULL) return (FALSE);

  /* load data from the photometry database file */
  sprintf (filename, "%s/%s", catdir, argv[1]);
  if (!fits_read_header (filename, &catalog1.header)) {
    fprintf (stderr, "no stars in %s\n", filename);
    return (FALSE);
  }
  f = fopen (filename, "r");
  if (f == (FILE *) NULL) {
    fprintf (stderr, "no stars in %s\n", filename);
    return (FALSE);
  }
  fseek (f, catalog1.header.size, SEEK_SET); 
  
  /** find number of stars, measurements **/
  catalog1.Naverage = 0;
  fits_scan (&catalog1.header, "NSTARS", "%d", 1, &catalog1.Naverage);
  if (catalog1.Naverage == 0) {
    fprintf (stderr, "no stars in catalog %s\n", filename);
    fclose (f);
    return (TRUE);
  }
  ALLOCATE (catalog1.average, Average, catalog1.Naverage);
  Nitems = catalog1.Naverage;
  nitems = Fread (catalog1.average, sizeof(Average), Nitems, f, "average");
  if (nitems != Nitems) {
    fprintf (stderr, "ERROR: failed to read data from catalog file %s (1)\n", filename);
    fclose (f);
    free (catalog1.average);
    return (FALSE);
  }
  fclose (f);
  fprintf (stderr, "read %d stars from phot catalog file %s\n", catalog1.Naverage, filename);

  /* load data from the GSC file */
  sprintf (filename, "%s/%s", gscdir, argv[1]);
  f = fopen (filename, "r");
  if (f == (FILE *) NULL) {
    fprintf (stderr, "no stars in %s\n", filename);
    free (catalog1.average);
    return (FALSE);
    /* return (FALSE); */
  }
  
  nstar = 0;
  NSTARS = DNSTARS;
  ALLOCATE (tbuffer, char, (BLOCK*BYTES_STAR));
  ALLOCATE (catalog2.average, Average, NSTARS);
  Nbytes = BLOCK*BYTES_STAR;
  while ((nbytes = fread (tbuffer, 1, Nbytes, f)) > 0) {
    for (i = 0; i < nbytes / BYTES_STAR; i++) {
      dparse (&R, 1, &tbuffer[i*BYTES_STAR]);
      dparse (&D, 2, &tbuffer[i*BYTES_STAR]);
      dparse (&M, 3, &tbuffer[i*BYTES_STAR]);
      catalog2.average[nstar].R = R;
      catalog2.average[nstar].D = D;
      catalog2.average[nstar].M = M * 1000.0;
      nstar++;
      if (nstar == NSTARS - 1) {
	NSTARS += DNSTARS;
	REALLOCATE (catalog2.average, Average, NSTARS);
      }
    }
  }
  free (tbuffer);
  REALLOCATE (catalog2.average, Average, MAX (nstar, 1));
  catalog2.Naverage = nstar;
  fclose (f);

  /* sort data in order of RA */
  sortave (catalog1.average, catalog1.Naverage);
  sortave (catalog2.average, catalog2.Naverage);

  /* data has been loaded, use gcompare algorithm to match */
  compare (&catalog1, &catalog2, rvec, dvec, mvec, drvec, ddvec, dmvec, radius);

  free (catalog1.average);
  free (catalog2.average);

  return (TRUE);

}
  
