# include "dvoshell.h"

/* This function uses the 'find_match' algorithm to select the objects of interest.
   Each entry in the match vectors (RA, DEC, RADIUS) yields a result value -- if no
   source matches, the resulting fields are all NAN.

   * choose the sky regions based on the provided RA,DEC points
   * loop over the catalogs
   * within a catalog, use the find_match code to find the matching coordinates
   * use dbExtractAverages to get the fields for the matched entry 
 */

int avmatch (int argc, char **argv) {
  
  off_t i, j, n, m, *index;
  int N, Ncat, Npts, NPTS, last, Nfields, Nsecfilt;
  int VERBOSE;
  char name[1024];
  void *Signal;
  float RADIUS;

  Catalog catalog;

  Vector **vec, *RAvec, *DECvec;
  dbField *fields;
  dbValue *values;
  SkyList *skylist;

  /* defaults */
  vec = NULL;
  fields = NULL;
  values = NULL;
  skylist = NULL;

  if ((N = get_argument (argc, argv, "-h"))) goto help;
  if ((N = get_argument (argc, argv, "--help"))) goto help;
  if (argc < 5) goto help;

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

  dvo_catalog_init (&catalog, TRUE);

  /* load photcode information */
  if (!InitPhotcodes ()) goto escape;
  Nsecfilt = GetPhotcodeNsecfilt ();

  // get vectors corresponding to coordinates of interest
  if ((RAvec  = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help;
  if ((DECvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help;
  RADIUS = atof (argv[3]);

  // strip off RA, DEC, RADIUS arguments
  remove_argument (1, &argc, argv);
  remove_argument (1, &argc, argv);
  remove_argument (1, &argc, argv);

  // parse the fields to be extracted and returned
  fields = dbCmdlineFields (argc, argv, DVO_TABLE_AVERAGE, &last, &Nfields);
  if (fields == NULL) goto help;
  if ((Nfields == 0) || (last != argc)) {
    dbFreeFields (fields, Nfields);
    dvo_catalog_free (&catalog);
    goto help;
  }

  /* load regions which contain all supplied RA,DEC coordinates */
  if ((skylist = SelectRegionsByCoordVectors (RAvec, DECvec)) == NULL) goto escape;

  /* create output storage vectors */
  NPTS = RAvec->Nelements;
  ALLOCATE (values, dbValue, Nfields);
  ALLOCATE (vec, Vector *, Nfields);
  for (i = 0; i < Nfields; i++) {
    if (ISNUM(fields[i].name[0])) {
      sprintf (name, "v_%s", fields[i].name);
    } else {
      sprintf (name, "%s", fields[i].name);
    }
    if ((vec[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) goto escape;
    ResetVector (vec[i], fields[i].type, NPTS);
    for (n = 0; n < vec[i][0].Nelements; n++) {
      if (vec[i][0].type == OPIHI_FLT) {
	vec[i][0].elements.Flt[n] = NAN;
      } else {
	vec[i][0].elements.Int[n] = 0; // or NAN_INT?
      }
    }
  }
  ALLOCATE (index, off_t, NPTS);

  // grab data from all selected sky regions
  Signal = signal (SIGINT, handle_interrupt);
  interrupt = FALSE;
  for (i = 0; (i < skylist[0].Nregions) && !interrupt; i++) {
    /* lock, load, unlock catalog */
    catalog.filename = skylist[0].filename[i];
    catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
    catalog.Nsecfilt = 0;

    if (VERBOSE) gprint (GP_ERR, "trying %s (%lld of %lld)\n", catalog.filename, (long long) i, (long long) skylist[0].Nregions);
      
    // an error exit status here is a significant error
    if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) {
      gprint (GP_ERR, "ERROR: failure to open catalog file %s\n", catalog.filename);
      exit (2);
    }
    dvo_catalog_unlock (&catalog);  /// we can unlock here since this is read-only (do not block other access)
    if (catalog.Naverage == 0) {
      dvo_catalog_free (&catalog);
      continue;
    }

    find_matches_by_vectors (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);

    for (j = 0; (j < NPTS) && !interrupt; j++) {
      Ncat = index[j];
      Npts = j;

      if (Ncat == -1) continue;
      if (Ncat == -2) continue;

      // XXX set a 'found' vector to double check we catch everything?
      m = catalog.average[Ncat].measureOffset;

      // reset counters for saved fields, extract fields
      dbExtractAveragesInit (); 
      for (n = 0; n < Nfields; n++) {
	values[n] = dbExtractAverages (&catalog.average[Ncat], &catalog.secfilt[Ncat*Nsecfilt], &catalog.measure[m], &fields[n]);
      }

      // XXX if we are allowed to return more rows than the supplied RA,DEC we will need to create an output RA,DEC
      for (n = 0; n < Nfields; n++) {
	if (vec[n][0].type == OPIHI_FLT) {
	  vec[n][0].elements.Flt[Npts] = values[n].Flt;
	} else {
	  vec[n][0].elements.Int[Npts] = values[n].Int;
	}
      }
    }
    dvo_catalog_free (&catalog);
  }
  signal (SIGINT, Signal);
  interrupt = FALSE;

  if (vec) free (vec);
  if (values) free (values);
  dbFreeFields (fields, Nfields);
  SkyListFree (skylist);
  return (TRUE);

 escape:
  if (vec) free (vec);
  if (values) free (values);
  dbFreeFields (fields, Nfields);
  SkyListFree (skylist);
  return (FALSE);

 help:
  gprint (GP_ERR, "USAGE: avmatch (RA) (DEC) (RADIUS) field[,field,field...]\n");

  if ((argc > N + 1) && !strcasecmp (argv[N+1], "fields")) {
    gprint (GP_ERR, "  RA : right ascension (J2000) [degrees]\n");
    gprint (GP_ERR, "  DEC : declination [degrees]\n");
    gprint (GP_ERR, "  GLON : galactic longitude [degrees]\n");
    gprint (GP_ERR, "  GLAT : galactic latitude [degrees]\n");
    gprint (GP_ERR, "  ELON : ecliptic longitude [degrees]\n");
    gprint (GP_ERR, "  ELAT : ecliptic latitude [degrees]\n");
    gprint (GP_ERR, "  dRA : ra scatter [degrees]\n");
    gprint (GP_ERR, "  dDEC : dec scatter [degrees]\n");
    gprint (GP_ERR, "  uRA : proper motion in ra [arcseconds]\n");
    gprint (GP_ERR, "  uDEC : proper motion in dec [arcseconds]\n");
    gprint (GP_ERR, "  duRA : proper motion error in ra [arcseconds]\n");
    gprint (GP_ERR, "  duDEC : proper motion error in dec [arcseconds]\n");
    gprint (GP_ERR, "  PAR : parallax\n");
    gprint (GP_ERR, "  dPAR : parallax error \n");

    gprint (GP_ERR, "  ChiSqPos : chi square of position fit \n");
    gprint (GP_ERR, "  ChiSqPM  : chi square of proper-motion fit \n");
    gprint (GP_ERR, "  ChiSqPar : chi square of parallax fit \n");

    gprint (GP_ERR, "  Tmean : mean epoch (reference for proper motion)\n");
    gprint (GP_ERR, "  Trange : range of times used for proper motion/parallax fit\n");

    gprint (GP_ERR, "  Nmeas : number of measurements\n");
    gprint (GP_ERR, "  Nmiss : number of non-detections\n");
    gprint (GP_ERR, "  Npos  : number of measurments used for astrometry\n");
    gprint (GP_ERR, "  Nastrom  : number of measurments used for astrometry (= Npos)\n");

    gprint (GP_ERR, "  flags     : object flags\n");
    gprint (GP_ERR, "  objflags  : object flags\n");
    gprint (GP_ERR, "  obj_flags : object flags\n");

    gprint (GP_ERR, "  objID : object ID (32 bit, unique in catalog)\n");
    gprint (GP_ERR, "  catID : catalog ID (32 bit)\n");
    gprint (GP_ERR, "  extID_hi : external ID (upper 32 of 64 bit) -- eg, PSPS ID\n");
    gprint (GP_ERR, "  extID_lo : external ID (lower 32 of 64 bit) -- eg, PSPS ID\n");

    gprint (GP_ERR, "  <photcode>:ave : average magnitude for <photcode>\n");
    gprint (GP_ERR, "  <photcode>:ref : reference magnitude system for <photcode>\n");
    gprint (GP_ERR, "  <photcode>:inst : first instrumental magnitude for <photcode>\n");
    gprint (GP_ERR, "  <photcode>:cat : first catalog magnitude for <photcode>\n");
    gprint (GP_ERR, "  <photcode>:sys : first system magnitude for <photcode>\n");
    gprint (GP_ERR, "  <photcode>:rel : first relative magnitude for <photcode>\n");
    gprint (GP_ERR, "  <photcode>:cal : first calibrated magnitude for <photcode> \n");
    gprint (GP_ERR, "  <photcode>:err : magnitude error for photcode\n");
    gprint (GP_ERR, "  <photcode>:chisq : raw chi-square of magnitude fit\n");
    gprint (GP_ERR, "  <photcode>:ncode : number of measurements matching photcode \n");
    gprint (GP_ERR, "  <photcode>:nphot : number of measurements used for average magnitude in this photcode\n");

    // gprint (GP_ERR, "  Xp : NOT VALID\n");
    // gprint (GP_ERR, "  type : dophot type (unused)\n");
    // gprint (GP_ERR, "  typefrac : dophot type fraction (unused)\n");
    return (FALSE);
  }
  gprint (GP_ERR, " avextract --help fields : for a complete listing of allowed fields\n");
  return (FALSE);
}
