Index: trunk/Ohana/src/opihi/dvo/Makefile
===================================================================
--- trunk/Ohana/src/opihi/dvo/Makefile	(revision 40522)
+++ trunk/Ohana/src/opihi/dvo/Makefile	(revision 40523)
@@ -46,4 +46,5 @@
 $(SRC)/avextract.$(ARCH).o	  	\
 $(SRC)/avmatch.$(ARCH).o	  	\
+$(SRC)/avselect.$(ARCH).o	  	\
 $(SRC)/avperiodogram.$(ARCH).o	  	\
 $(SRC)/avperiodomatch.$(ARCH).o	  	\
Index: trunk/Ohana/src/opihi/dvo/avmatch.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/avmatch.c	(revision 40522)
+++ trunk/Ohana/src/opihi/dvo/avmatch.c	(revision 40523)
@@ -209,5 +209,5 @@
     }
 
-    find_matches_by_vectors (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);
+    find_matches_by_vectors_closest (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);
 
     for (j = 0; (j < NPTS) && !interrupt; j++) {
Index: trunk/Ohana/src/opihi/dvo/avperiodomatch.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/avperiodomatch.c	(revision 40522)
+++ trunk/Ohana/src/opihi/dvo/avperiodomatch.c	(revision 40523)
@@ -264,5 +264,5 @@
     dvo_catalog_unlock (&catalog);
 
-    find_matches_by_vectors (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, idxValue);
+    find_matches_by_vectors_closest (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, idxValue);
 
     // Scan the catalog for objects which match the WHERE clause
Index: trunk/Ohana/src/opihi/dvo/avselect.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/avselect.c	(revision 40523)
+++ trunk/Ohana/src/opihi/dvo/avselect.c	(revision 40523)
@@ -0,0 +1,405 @@
+# 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 avselect (int argc, char **argv) {
+  
+  off_t i, j, n, m;
+  int N, Ncat, Npts, NPTS, last, Nfields, Nsecfilt, Ninvec;
+  int VERBOSE;
+  char name[1024];
+  float RADIUS;
+
+  Catalog catalog;
+
+  Vector **vec, **invec, *RAvec, *DECvec, *IDXvec, *RADvec;
+  dbField *fields;
+  dbValue *values;
+  SkyList *skylist;
+
+  /* defaults */
+  vec = NULL;
+  invec = NULL;
+  fields = NULL;
+  values = NULL;
+  skylist = NULL;
+  Ninvec = 0;
+
+  if ((N = get_argument (argc, argv, "-h"))) goto help;
+  if ((N = get_argument (argc, argv, "--help"))) goto help;
+
+  VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = TRUE;
+  }
+
+  IDXvec = NULL;
+  if ((N = get_argument (argc, argv, "-index"))) {
+    remove_argument (N, &argc, argv);
+    IDXvec = SelectVector (argv[N], ANYVECTOR, TRUE);
+    if (IDXvec == NULL) goto help;
+    remove_argument (N, &argc, argv);
+  }
+
+  RADvec = NULL;
+  if ((N = get_argument (argc, argv, "-radius"))) {
+    remove_argument (N, &argc, argv);
+    RADvec = SelectVector (argv[N], ANYVECTOR, TRUE);
+    if (RADvec == NULL) goto help;
+    remove_argument (N, &argc, argv);
+  }
+
+  int PARALLEL = FALSE;
+  if ((N = get_argument (argc, argv, "-parallel"))) {
+    remove_argument (N, &argc, argv);
+    PARALLEL = TRUE;
+  }
+
+  // dump results directly to fits file (esp for parallel dvo)
+  char *CoordsFile = NULL;
+  if ((N = get_argument (argc, argv, "-coords"))) {
+    remove_argument (N, &argc, argv);
+    CoordsFile = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (!CoordsFile && (argc < 5)) goto help;
+  if ( CoordsFile && (argc < 3)) goto help;
+
+  dvo_catalog_init (&catalog, TRUE);
+
+  /* load photcode information */
+  if (!InitPhotcodes ()) goto escape;
+  Nsecfilt = GetPhotcodeNsecfilt ();
+
+  // init locally static variables (time refs)
+  dbExtractAveragesInit (); 
+
+  RAvec  = NULL;
+  DECvec = NULL;
+  // get vectors corresponding to coordinates of interest
+  if (CoordsFile) {
+    // read RAvec, DECvec from coords file (1st 2 fields?)
+    Ninvec = 0;
+    invec = ReadVectorTableFITS (CoordsFile, "COORDS", &Ninvec);
+    RAvec = invec[0];
+    DECvec = invec[1];
+  } else {
+    if ((RAvec  = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help;
+    if ((DECvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help;
+    // strip off RA & DEC
+    remove_argument (1, &argc, argv);
+    remove_argument (1, &argc, argv);
+  }
+
+  /* load regions which contain all supplied RA,DEC coordinates */
+  if ((skylist = SelectRegionsByCoordVectors (RAvec, DECvec)) == NULL) goto escape;
+
+  // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results
+  if (PARALLEL && !HOST_ID) {
+
+    // We need to copy the args to a temp array and modify them so that we send the
+    // correct set to the remote client.  The args list looks like this:
+    // if (!CoordsFile) : avmatch (RADIUS) field, ... [we removed RA & DEC above]
+    // if ( CoordsFile) : avmatch (RADIUS) field, ... [because we stripped off the -coords filename elements]
+
+    // allocate the temp array and copy all but (RA) (DEC)
+    int targc = 0;
+    char **targv = NULL;
+    ALLOCATE (targv, char *, argc + 2);
+    for (i = 0; i < argc; i++) {
+      targv[targc] = strcreate (argv[i]);
+      targc ++;
+    }
+
+    // if not specified, create the coords.fits input file
+    // NOTE: RAvec, DECvec were set above
+    if (!CoordsFile) {
+      ALLOCATE (vec, Vector *, 2);
+      vec[0] = RAvec;
+      vec[1] = DECvec;
+
+      CoordsFile = abspath("coords.fits", 1024);
+      int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
+      if (!status) goto escape;
+    }
+
+    // add the coords file to the args list
+    targv[targc+0] = strcreate ("-coords");
+    targv[targc+1] = CoordsFile; // this gets freed with targv
+    targc += 2;
+    
+    // if needed, add the index vector to the args list
+    if (IDXvec) {
+      REALLOCATE (targv, char *, targc + 2);
+      targv[targc+0] = strcreate ("-index");
+      targv[targc+1] = strcreate (IDXvec[0].name);
+      targc += 2;
+    }      
+    if (RADvec) {
+      REALLOCATE (targv, char *, targc + 2);
+      targv[targc+0] = strcreate ("-radius");
+      targv[targc+1] = strcreate (RADvec[0].name);
+      targc += 2;
+    }      
+
+    // I need to pass the RA & DEC vectors to the remote clients...
+    int status = HostTableParallelOps (skylist, targc, targv, RESULT_FILE, TRUE, 0, VERBOSE);
+    if (vec) free (vec);
+    
+    // free up targv
+    for (i = 0; i < targc; i++) {
+      free (targv[i]);
+    }
+    free (targv);
+
+    return status;
+  }
+
+  RADIUS = atof (argv[1]);
+  remove_argument (1, &argc, argv);
+
+  // parse the fields to be extracted and returned
+  // XXX with this block here, we do not get a check on the syntax until after the clients have launched
+  // (see mmatch.c for a fix)
+  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;
+  }
+
+  // check the requested fields
+  int needMeasure = dbFieldNeedMeasure (fields, Nfields);
+  int needLensobj = dbFieldNeedLensobj (fields, Nfields);
+  int needStarpar = dbFieldNeedStarpar (fields, Nfields, TRUE);
+
+  /* create output storage vectors */
+  Npts = 0;
+  NPTS = 1000;
+  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);
+  }
+  if (IDXvec) {
+    ResetVector (IDXvec, OPIHI_INT, NPTS);
+  }
+  if (RADvec) {
+    ResetVector (RADvec, OPIHI_FLT, NPTS);
+  }
+
+  // grab data from all selected sky regions
+  struct sigaction *old_sigaction = SetInterrupt();
+  for (i = 0; (i < skylist[0].Nregions) && !interrupt; i++) {
+
+    // does this host ID match the desired location for the table?
+    if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue;
+
+    /* lock, load, unlock catalog */
+    char hostfile[1024];
+    snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name);
+
+    dvo_catalog_init (&catalog, TRUE);
+    catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i];
+    catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_SECFILT;
+    catalog.catflags |= needMeasure ? DVO_LOAD_MEASURE : DVO_SKIP_MEASURE;
+    catalog.catflags |= needLensobj ? DVO_LOAD_LENSOBJ : DVO_SKIP_LENSOBJ;
+    catalog.catflags |= needStarpar ? DVO_LOAD_STARPAR : DVO_SKIP_STARPAR;
+    catalog.Nsecfilt = 0;
+
+    if (VERBOSE) gprint (GP_ERR, "trying %s ("OFF_T_FMT" of "OFF_T_FMT")\n", catalog.filename,  i,  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 all matches within the radius. index returns the elements of catalog which are matches
+    off_t Nresult;
+    AvselectResult *result = find_matches_by_vectors_allmatch (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, &Nresult);
+
+    for (j = 0; (j < Nresult) && !interrupt; j++) {
+      Ncat = result[j].Ncat;
+      Average *average = &catalog.average[Ncat];
+
+      m = average->measureOffset;
+      Measure *measure = needMeasure ? &catalog.measure[m] : NULL;
+
+      m = average->lensobjOffset;
+      Lensobj *lensobj = needLensobj ? &catalog.lensobj[m] : NULL;
+
+      m = average->starparOffset;
+      StarPar *starpar = needStarpar ? &catalog.starpar[m] : NULL;
+
+      m = Ncat*Nsecfilt;
+      SecFilt *secfilt = &catalog.secfilt[m];
+
+      // reset counters for saved fields, extract fields
+      dbExtractAveragesInitAve (); 
+      for (n = 0; n < Nfields; n++) {
+	values[n] = dbExtractAverages (average, secfilt, measure, lensobj, starpar, NULL, &fields[n]);
+      }
+
+      // set resulting values
+      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;
+	}
+      }
+      // set (optional) IDXvec and RADvec
+      if (IDXvec) {
+	IDXvec[0].elements.Int[Npts] = result[j].Nseq;;
+      }
+      if (RADvec) {
+	RADvec[0].elements.Flt[Npts] = result[j].Roff;;
+      }
+
+      // extend length of output vectors
+      Npts++;
+      if (Npts >= NPTS) {
+	NPTS += 2000;
+	for (n = 0; n < Nfields; n++) {
+	  REALLOCATE (vec[n][0].elements.Flt, opihi_flt, NPTS);
+	}
+	if (IDXvec) {
+	  REALLOCATE (IDXvec[0].elements.Int, opihi_int, NPTS);
+	}
+	if (RADvec) {
+	  REALLOCATE (RADvec[0].elements.Flt, opihi_flt, NPTS);
+	}
+      }
+    }
+    dvo_catalog_free (&catalog);
+  }
+  ClearInterrupt (old_sigaction);
+
+  for (i = 0; i < Nfields; i++) {
+    vec[i][0].Nelements = Npts;
+  }
+  if (IDXvec) {
+    IDXvec[0].Nelements = Npts;
+  }
+  if (RADvec) {
+    RADvec[0].Nelements = Npts;
+  }
+
+  // write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere)
+  // only write the fields which were in a valid catalog
+  if (RESULT_FILE) {
+    if (IDXvec) {
+      // extend the array by one to hold index array
+      Nfields ++;
+      REALLOCATE (vec, Vector *, Nfields);
+      vec[Nfields-1] = IDXvec;
+    }
+    if (RADvec) {
+      // extend the array by one to hold index array
+      Nfields ++;
+      REALLOCATE (vec, Vector *, Nfields);
+      vec[Nfields-1] = RADvec;
+    }
+    int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nfields, FALSE, FALSE, NULL, 0);
+    if (!status) goto escape;
+  }
+
+  if (vec) free (vec);
+  if (values) free (values);
+  if (invec) FreeVectorArray (invec, Ninvec);
+  dbFreeFields (fields, Nfields);
+  SkyListFree (skylist);
+  return (TRUE);
+
+ escape:
+  if (vec) free (vec);
+  if (values) free (values);
+  if (invec) FreeVectorArray (invec, Ninvec);
+  dbFreeFields (fields, Nfields);
+  SkyListFree (skylist);
+  return (FALSE);
+
+ help:
+  gprint (GP_ERR, "USAGE: avselect (RA) (DEC) (RADIUS) field[,field,field...]\n");
+  gprint (GP_ERR, "   OR: avselect -coords (filename.fits) (RADIUS) field[,field,field...]\n");
+  gprint (GP_ERR, "   RADIUS is in arcseconds\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, "  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);
+}
Index: trunk/Ohana/src/opihi/dvo/find_matches.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/find_matches.c	(revision 40522)
+++ trunk/Ohana/src/opihi/dvo/find_matches.c	(revision 40523)
@@ -5,5 +5,5 @@
 // attempt to match every RA,DEC entry with an entry in the catalog.  the result is stored in 'index'
 // index >= 0 : valid match; index == -2 : no valid match; index == -1 : not contained by this catalog
-int find_matches_by_vectors (SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *index) {
+int find_matches_by_vectors_closest (SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *index) {
 
   off_t i, j, J, Jmin, Npoints, Nmatch;
@@ -148,2 +148,141 @@
   return (TRUE);
 }
+
+// attempt to match every RA,DEC entry with an entry in the catalog.  the result is stored in 'result'
+AvselectResult *find_matches_by_vectors_allmatch (SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *nresult) {
+
+  off_t i, j, J;
+  double dX, dY, dR, Doff, Dmin, Dmax, Rmin, Rmax;
+  int status;
+  double RADIUS2;
+  Coords tcoords;
+
+  // XXX dvoconvert does not correctly maintain 'sorted' 
+  // assert(catalog[0].sorted);
+
+  off_t Npoints = RAvec->Nelements;
+
+  /** allocate local arrays (points) **/
+  ALLOCATE_PTR (X1, double, Npoints);
+  ALLOCATE_PTR (Y1, double, Npoints);
+  ALLOCATE_PTR (N1, off_t,  Npoints);
+
+  ALLOCATE_PTR (inCatalog, int,  Npoints);
+
+  /** allocate local arrays (catalog) **/
+  off_t Nave = catalog[0].Naverage;
+  ALLOCATE_PTR (X2, double, Nave);
+  ALLOCATE_PTR (Y2, double, Nave);
+  ALLOCATE_PTR (N2, off_t,  Nave);
+
+  off_t Nresult = 0;
+  off_t NRESULT = 1000;
+  ALLOCATE_PTR (result, AvselectResult, NRESULT);
+
+  /* project onto rectilinear grid with 1 arcsec pixels. the choice of ARC projection has
+   * the advantage that every point in R,D has a mapping to a unique X,Y.  However, note
+   * that not all possible X,Y points map back to R,D and the local plate scale changes
+   * far from the projection pole. We use the center of the region (catalog) for crval1,2.
+   */
+  InitCoords (&tcoords, "DEC--ARC");
+  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
+  if (region[0].Dmax < 90) {
+    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
+  } else {
+    tcoords.crval2 = 90.0;
+  }
+  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+
+  // this region includes a boundary layer of size RADIUS
+  if (fabs(region[0].Dmin) < fabs(region[0].Dmax)) {
+    Doff = RAD_DEG*region[0].Dmax;
+  } else {
+    Doff = RAD_DEG*region[0].Dmin;
+  }    
+  if (Doff < 80) {
+    Rmin = region[0].Rmin - RADIUS / 3600.0 / cos(Doff);
+    Rmax = region[0].Rmax + RADIUS / 3600.0 / cos(Doff);
+  } else {
+    Rmin = 0.0;
+    Rmax = 360.0;
+  }
+  Dmin = region[0].Dmin - RADIUS / 3600.0;
+  Dmax = region[0].Dmax + RADIUS / 3600.0;
+
+  // identify the entries contained by this catalog
+  for (i = 0; i < Npoints; i++) {
+    inCatalog[i] = FALSE;
+    if (RAvec->elements.Flt[i] < Rmin) continue;
+    if (RAvec->elements.Flt[i] > Rmax) continue;
+    if (DECvec->elements.Flt[i] < Dmin) continue;
+    if (DECvec->elements.Flt[i] > Dmax) continue;
+    inCatalog[i] = TRUE;
+  }
+
+  /* build spatial index (RA sort) referencing input array sequence */
+  for (i = 0; i < Npoints; i++) {
+    // we need to have a finite number for the sort below; inCatalog == 0 entries are skipped
+    // in the matching process
+    X1[i] = Y1[i] = 0.0; 
+    N1[i] = i;
+    if (!inCatalog[i]) continue;
+    status = RD_to_XY (&X1[i], &Y1[i], RAvec->elements.Flt[i], DECvec->elements.Flt[i], &tcoords);
+    assert (status);
+  }
+  if (Npoints > 1) sort_coords_index (X1, Y1, N1, Npoints);
+
+  /* build spatial index (RA sort) */
+  for (i = 0; i < Nave; i++) {
+    RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
+    N2[i] = i;
+  }
+  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
+
+  /* choose a radius for matches */
+  RADIUS2 = RADIUS*RADIUS;
+
+# define NEXTi { i++; continue; }
+# define NEXTj { j++; continue; }
+                                                                                                                                                                                 
+  /** find matched stars **/
+  for (i = j = 0; (i < Npoints) && (j < Nave); ) {
+    if (!inCatalog[N1[i]]) NEXTi;
+    if (!finite(X1[i]) || !finite(Y1[i])) NEXTi;
+    if (!finite(X2[j]) || !finite(Y2[j])) NEXTj;
+    
+    /* negative dX: j is too large */
+    dX = X1[i] - X2[j];
+    if (dX <= -1.02*RADIUS) NEXTi;
+
+    /* positive dX, i is too large */
+    if (dX >= 1.02*RADIUS) NEXTj;
+
+    /* within match range; look for matches */
+    for (J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) {
+      /* find closest match for this detection */
+      dX = X1[i] - X2[J];
+      dY = Y1[i] - Y2[J];
+      dR = dX*dX + dY*dY;
+      if (dR > RADIUS2) continue;
+
+      /*** a match is found, set the result for this entry ***/
+      result[Nresult].Ncat = N2[J];
+      result[Nresult].Nseq = N1[i];
+      result[Nresult].Roff = sqrt(dR);
+
+      Nresult ++;
+      CHECK_REALLOCATE (result, AvselectResult, NRESULT, Nresult, 1000);
+    }
+    NEXTi;
+  }
+
+  free (X1);
+  free (Y1);
+  free (N1);
+  free (N2);
+  free (X2);
+  free (Y2);
+
+  *nresult = Nresult;
+  return (result);
+}
Index: trunk/Ohana/src/opihi/dvo/init.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/init.c	(revision 40522)
+++ trunk/Ohana/src/opihi/dvo/init.c	(revision 40523)
@@ -3,4 +3,5 @@
 int avextract       PROTO((int, char **));
 int avmatch         PROTO((int, char **));
+int avselect        PROTO((int, char **));
 int avperiodogram   PROTO((int, char **));
 int avperiodomatch  PROTO((int, char **));
@@ -69,4 +70,5 @@
   {1, "avextract",   avextract,    "extract average data values"},
   {1, "avmatch",     avmatch,      "extract average data values matched to RA,DEC points"},
+  {1, "avselect",    avselect,     "extract average data values within range of a list of RA,DEC points"},
   {1, "avperiodogram",  avperiodogram, "perform periodogram on objects based on restrictions"},
   {1, "avperiodomatch", avperiodomatch, "perform periodogram on objects based on ra,dec list"},
Index: trunk/Ohana/src/opihi/dvo/mmatch.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/mmatch.c	(revision 40522)
+++ trunk/Ohana/src/opihi/dvo/mmatch.c	(revision 40523)
@@ -84,8 +84,7 @@
 
   // parse the fields to be extracted and returned
-  int first = 3;
-  if (CoordsFile) {
-    first = 1;
-  }
+  // this is a syntax check that we can perform before reading the input vectors
+  // NOTE: This block is missing in avmatch / avselect
+  int first = CoordsFile ? 1 : 3;
   fields = dbCmdlineFields (argc-first, &argv[first], DVO_TABLE_MEASURE, &last, &Nfields);
   if (fields == NULL) goto help;
@@ -166,4 +165,5 @@
     }      
 
+    // NOTE: image metadata is only needed to mmatch (not avmatch) since measures are tied to images, not averages
     if (loadImages) {
       Image *image;
@@ -245,4 +245,5 @@
     char hostfile[1024];
     snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name);
+
     dvo_catalog_init (&catalog, TRUE);
     catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i];
@@ -265,5 +266,5 @@
 
     // returns the matches to AVERAGE objects; what do we do for each MEASURE?
-    find_matches_by_vectors (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);
+    find_matches_by_vectors_closest (skylist[0].regions[i], &catalog, RAvec, DECvec, RADIUS, index);
 
     for (j = 0; (j < Nelem) && !interrupt; j++) {
@@ -331,4 +332,5 @@
   ClearInterrupt (old_sigaction);
 
+  // XXXX ???? this seems odd
   for (n = 0; n < Nfields; n++) {
     ResetVector (vec[n], fields[n].type, Npts);
Index: trunk/Ohana/src/opihi/include/dvoshell.h
===================================================================
--- trunk/Ohana/src/opihi/include/dvoshell.h	(revision 40522)
+++ trunk/Ohana/src/opihi/include/dvoshell.h	(revision 40523)
@@ -52,4 +52,10 @@
 } PeriodogramResult;
 
+typedef struct {
+  off_t Ncat;
+  off_t Nseq;
+  float Roff;
+} AvselectResult;
+
 /*** Periodogram functions (avperiodogram and avperiodomatch) ***/
 void PeriodogramResultFree (PeriodogramResult *result);
@@ -98,5 +104,6 @@
 SkyList      *SelectRegionsByCoordVectors PROTO((Vector *RA, Vector *DEC));
 
-int           find_matches_by_vectors PROTO((SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *index));
+int             find_matches_by_vectors_closest PROTO((SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *index));
+AvselectResult *find_matches_by_vectors_allmatch PROTO((SkyRegion *region, Catalog *catalog, Vector *RAvec, Vector *DECvec, float RADIUS, off_t *Nresult));
 
 int           SetImageSelection     PROTO((int mode, SkyRegionSelection *selection));