Index: /trunk/Ohana/src/addstar/Makefile
===================================================================
--- /trunk/Ohana/src/addstar/Makefile	(revision 26279)
+++ /trunk/Ohana/src/addstar/Makefile	(revision 26280)
@@ -51,4 +51,5 @@
 $(SRC)/find_matches_closest.$(ARCH).o \
 $(SRC)/find_matches_refstars.$(ARCH).o \
+$(SRC)/find_matches_closest_refstars.$(ARCH).o \
 $(SRC)/find_subset.$(ARCH).o \
 $(SRC)/getgsc.$(ARCH).o \
@@ -197,4 +198,5 @@
 $(SRC)/get2mass_full.$(ARCH).o \
 $(SRC)/find_matches_refstars.$(ARCH).o \
+$(SRC)/build_links.$(ARCH).o \
 $(SRC)/args_load2mass.$(ARCH).o \
 $(SRC)/replace_match.$(ARCH).o \
Index: /trunk/Ohana/src/addstar/include/addstar.h
===================================================================
--- /trunk/Ohana/src/addstar/include/addstar.h	(revision 26279)
+++ /trunk/Ohana/src/addstar/include/addstar.h	(revision 26280)
@@ -172,4 +172,5 @@
 int        find_matches_closest   PROTO((SkyRegion *region, Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions options));
 int        find_matches_refstars  PROTO((SkyRegion *region, Stars **stars, int Nstars, Catalog *catalog, AddstarClientOptions options));
+int        find_matches_closest_refstars  PROTO((SkyRegion *region, Stars **stars, int Nstars, Catalog *catalog, AddstarClientOptions options));
 
 Stars    **find_subset            PROTO((SkyRegion *region, Stars *stars, int Nstars, int *NSTARS));
Index: /trunk/Ohana/src/addstar/src/addstar.c
===================================================================
--- /trunk/Ohana/src/addstar/src/addstar.c	(revision 26279)
+++ /trunk/Ohana/src/addstar/src/addstar.c	(revision 26280)
@@ -138,5 +138,9 @@
       case M_REFLIST:
 	subset = find_subset (skylist[0].regions[i], stars, Nstars, &Nsubset);
-	Nmatch += find_matches_refstars (skylist[0].regions[i], subset, Nsubset, &catalog, options);
+	if (options.closest) {
+	  Nmatch += find_matches_closest_refstars (skylist[0].regions[i], subset, Nsubset, &catalog, options);
+	} else {
+	  Nmatch += find_matches_refstars (skylist[0].regions[i], subset, Nsubset, &catalog, options);
+	}
 	if (Nsubset) free (subset);
 	break;
Index: /trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c
===================================================================
--- /trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c	(revision 26280)
+++ /trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c	(revision 26280)
@@ -0,0 +1,359 @@
+# include "addstar.h"
+
+int find_matches_closest_refstars (SkyRegion *region, Stars **stars, int NstarsIn, Catalog *catalog, AddstarClientOptions options) {
+
+  int i, j, n, N, J, Jmin, status;
+  double RADIUS, RADIUS2, Rmin;
+  double *X1, *Y1, *X2, *Y2;
+  double dX, dY, dR;
+  int *N1, *N2, *next_meas;
+  int Nave, NAVE, Nmeas, NMEAS, Nstars, Nmatch;
+  unsigned int objID, catID;
+  Coords tcoords;
+  int Nsecfilt;
+
+  if ((NREFSTAR_GROUP != 1) && (NREFSTAR_GROUP != 3)) {
+    fprintf (stderr, "ERROR: NREFSTAR_GROUP NOT SET!\n");
+    exit (1);
+  }
+
+  /* photcode data -- should not have to modify secfilt / average */
+  Nsecfilt = GetPhotcodeNsecfilt ();
+
+  /** allocate local arrays (stars) **/
+  ALLOCATE (X1, double, NstarsIn);
+  ALLOCATE (Y1, double, NstarsIn);
+  ALLOCATE (N1, int,    NstarsIn);
+
+  /** allocate local arrays (catalog) **/
+  NAVE = Nave = catalog[0].Naverage;
+  ALLOCATE (X2, double, NAVE);
+  ALLOCATE (Y2, double, NAVE);
+  ALLOCATE (N2, int, NAVE);
+  ALLOCATE (catalog[0].found, int, NAVE);
+  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
+
+  /* internal counters */
+  Nmatch = 0;
+  NMEAS = Nmeas = catalog[0].Nmeasure;
+  
+  // current max obj ID for this catalog
+  objID = catalog[0].objID;
+  catID = catalog[0].catID;
+
+  /* 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.
+   */
+  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.crpix1 = 0.0;
+  tcoords.crpix2 = 0.0;
+  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+  tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
+  tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
+  tcoords.Npolyterms = 0;
+  strcpy (tcoords.ctype, "RA---ARC");
+  
+  /* build spatial index (RA sort) referencing input array sequence */
+  Nstars = 0;
+  for (i = 0; i < NstarsIn; i++) {
+    status = RD_to_XY (&X1[Nstars], &Y1[Nstars], stars[i][0].average.R, stars[i][0].average.D, &tcoords);
+    if (!status) continue;
+    N1[Nstars] = i;
+    Nstars ++;
+  }
+  if (Nstars < 1) {
+    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename);
+    free (catalog[0].found);
+    free (X1);
+    free (Y1);
+    free (N1);
+    free (X2);
+    free (Y2);
+    free (N2);
+    return (0);
+  }
+  if (Nstars > 1) sort_coords_index (X1, Y1, N1, Nstars);
+  
+  /* 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;
+    catalog[0].found[N2[i]] = -1;
+  }
+  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
+
+  /* set up pointers for linked list of measure */
+  if (catalog[0].sorted && (catalog[0].Nmeasure == catalog[0].Nmeas_disk)) {
+    // this version is only valid if we have done a full catalog load, and if the catalog
+    // is sorted while processed
+    next_meas = init_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
+  } else {
+    next_meas = build_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
+  }    
+
+  /* choose a radius for matches */
+  if (options.radius == 0) {
+    RADIUS = 2.0; /* hardwired default for refstars */
+  } else {
+    RADIUS = options.radius; /* provided by config */
+  }
+  RADIUS2 = RADIUS*RADIUS;
+
+  /** find matched stars **/
+  // XXX could use NREFSTAR_GROUP to do this match more quicky
+  for (i = j = 0; (i < Nstars) && (j < Nave); ) {
+    if (!finite(X1[i]) || !finite(Y1[i])) { 
+      i++; 
+      continue;
+    }
+    if (!finite(X2[j]) || !finite(Y2[j])) { 
+      j++; 
+      continue;
+    }
+    
+    /* negative dX: j is too large */
+    dX = X1[i] - X2[j];
+    if (dX <= -1.02*RADIUS) {
+      i++;
+      continue;
+    }
+    /* positive dX, i is too large */
+    if (dX >= 1.02*RADIUS) {
+      j++;
+      continue;
+    }
+
+    /* within match range; look for matches */
+    Jmin = -1;
+    Rmin = RADIUS2;
+    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;
+      if (dR < Rmin) {
+	Rmin = dR;
+	Jmin  = J;
+      }
+    }
+
+    /* no match, try next detection */ 
+    if (Jmin == -1) {
+      i++;
+      continue;
+    }
+
+    /* make sure there is space for next entry */
+    if (Nmeas >= NMEAS) {
+      NMEAS = Nmeas + 1000;
+      REALLOCATE (next_meas, int, NMEAS);
+      REALLOCATE (catalog[0].measure, Measure, NMEAS);
+    }
+
+    /*** a match is found, add to average, measure ***/
+    Nmatch ++;
+    n = N2[Jmin];
+    N = N1[i];
+
+    /** in replace mode, search for entry and replace values M, dM, R, D */
+    // XXX this fails for unsorted catalogs, right?
+    if (options.replace && replace_match (&catalog[0].average[n], catalog[0].measure, stars[N])) continue;
+
+    /* add to end of measurement list */
+    add_meas_link (&catalog[0].average[n], next_meas, Nmeas, NMEAS);
+	
+    /** add measurements for this star **/
+
+    // set the new measurements
+    catalog[0].measure[Nmeas]          = stars[N][0].measure;
+
+    /** *** dR,dD now in arcsec *** **/
+    catalog[0].measure[Nmeas].dR       = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
+    catalog[0].measure[Nmeas].dD       = 3600.0*(catalog[0].average[n].D - stars[N][0].average.D);
+    catalog[0].measure[Nmeas].dbFlags  = 0;
+    catalog[0].measure[Nmeas].averef   = n;
+    catalog[0].measure[Nmeas].objID    = catalog[0].average[n].objID;
+    catalog[0].measure[Nmeas].catID    = catalog[0].catID;
+
+    // rationalize dR:
+    if (catalog[0].measure[Nmeas].dR > +180.0*3600.0) {
+      // average on high end of boundary, move star up
+      stars[N][0].average.R += 360.0;
+      catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
+    }
+    if (catalog[0].measure[Nmeas].dR < -180.0*3600.0) {
+      // average on low end of boundary, move star down
+      stars[N][0].average.R -= 360.0;
+      catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
+    }
+
+    // the reference data may not have per-object time specified
+    catalog[0].measure[Nmeas].t        = (TIMEREF == 0) ? stars[N][0].measure.t      : TIMEREF; /** careful : time_t vs e_time **/
+    catalog[0].measure[Nmeas].t_msec   = (TIMEREF == 0) ? stars[N][0].measure.t_msec : TIMEREF; /** careful : time_t vs e_time **/
+
+    // we can choose to accept the proper-motion and parallax from the reference catalog
+    if (ACCEPT_MOTION) {
+      catalog[0].average[n].uR         = stars[N][0].average.uR;
+      catalog[0].average[n].uD         = stars[N][0].average.uD;
+      catalog[0].average[n].duR        = stars[N][0].average.duR;
+      catalog[0].average[n].duD        = stars[N][0].average.duD;
+      catalog[0].average[n].P          = stars[N][0].average.P;
+      catalog[0].average[n].dP         = stars[N][0].average.dP;
+    }
+
+    /** don't update average / secfilt values for REF photcodes **/
+
+    /*** handle multiple stars */
+    /* this catalog star matches more than one image star */
+    // XXX should this be an average flag?
+    if (catalog[0].found[n] > -1) {
+      catalog[0].measure[catalog[0].found[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
+      catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_OBJ;
+    } else {
+      catalog[0].found[n] = Nmeas;
+    }
+
+    stars[N][0].found = Nmeas;
+    catalog[0].average[n].Nmeasure ++;
+    Nmeas ++;
+
+    // make this optional:
+    // update_coords (&catalog[0].average[n], &catalog[0].measure[0], next_meas);
+    i++;
+  }
+
+  /* Incorporate unmatched refcat stars.  Skip this step if we want to
+     require matches combined with -replace, this lets us keep the
+     reference up-to-date with known stars only */
+
+  for (i = 0; (i < Nstars) && !options.only_match; i+=NREFSTAR_GROUP) {
+    /* make sure there is space for next entry */
+    if (Nmeas >= NMEAS - NREFSTAR_GROUP) {
+      NMEAS = Nmeas + 1000;
+      REALLOCATE (next_meas, int, NMEAS);
+      REALLOCATE (catalog[0].measure, Measure, NMEAS);
+    }
+    if (Nave >= NAVE) {
+      NAVE = Nave + 1000;
+      REALLOCATE (catalog[0].average, Average, NAVE);
+      REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
+    }
+
+    N = N1[i];
+    if (stars[N][0].found >= 0) continue;
+    if (!IN_REGION (stars[N][0].average.R, stars[N][0].average.D)) continue;
+
+    catalog[0].average[Nave].R         	   = stars[N][0].average.R;
+    catalog[0].average[Nave].D         	   = stars[N][0].average.D;
+
+    catalog[0].average[Nave].Nmeasure      = NREFSTAR_GROUP;
+    catalog[0].average[Nave].Nmissing      = 0;
+    catalog[0].average[Nave].Nextend       = 0;
+
+    catalog[0].average[Nave].measureOffset = Nmeas;
+    catalog[0].average[Nave].missingOffset = -1;
+    catalog[0].average[Nave].extendOffset  = -1;
+
+    if (ACCEPT_MOTION) {
+      catalog[0].average[Nave].dR    	   = stars[N][0].average.dR;
+      catalog[0].average[Nave].dD    	   = stars[N][0].average.dD;
+      catalog[0].average[Nave].uR    	   = stars[N][0].average.uR;
+      catalog[0].average[Nave].uD    	   = stars[N][0].average.uD;
+      catalog[0].average[Nave].duR   	   = stars[N][0].average.duR;
+      catalog[0].average[Nave].duD   	   = stars[N][0].average.duD;
+      catalog[0].average[Nave].P     	   = stars[N][0].average.P;
+      catalog[0].average[Nave].dP    	   = stars[N][0].average.dP;
+    } else {
+      catalog[0].average[Nave].dR    	   = 0;
+      catalog[0].average[Nave].dD    	   = 0;
+      catalog[0].average[Nave].uR    	   = 0;
+      catalog[0].average[Nave].uD    	   = 0;
+      catalog[0].average[Nave].duR   	   = 0;
+      catalog[0].average[Nave].duD   	   = 0;
+      catalog[0].average[Nave].P     	   = 0;
+      catalog[0].average[Nave].dP    	   = 0;
+      catalog[0].average[Nave].Xp    	   = 0;
+    }
+
+    catalog[0].average[Nave].Xp            = 0;
+    catalog[0].average[Nave].ChiSq    	   = 0.0;
+    catalog[0].average[Nave].Npos    	   = 0;
+
+    catalog[0].average[Nave].objID     	   = objID;
+    catalog[0].average[Nave].catID     	   = catID;
+    catalog[0].average[Nave].flags     	   = 0;
+
+    objID ++;
+
+    for (j = 0; j < Nsecfilt; j++) {
+      catalog[0].secfilt[Nave*Nsecfilt+j].M  	= NAN;
+      catalog[0].secfilt[Nave*Nsecfilt+j].dM 	= NAN;
+      catalog[0].secfilt[Nave*Nsecfilt+j].Xm 	= NAN_S_SHORT;
+      catalog[0].secfilt[Nave*Nsecfilt+j].M_20 	= NAN_S_SHORT;
+      catalog[0].secfilt[Nave*Nsecfilt+j].M_80 	= NAN_S_SHORT;
+      catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0;
+      catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0;
+    }
+
+    for (j = 0; j < NREFSTAR_GROUP; j++) {
+      N = N1[i + j];
+
+      catalog[0].measure[Nmeas]          = stars[N][0].measure;
+
+      catalog[0].measure[Nmeas].dR       = 0.0; // astrometric offset, not error
+      catalog[0].measure[Nmeas].dD       = 0.0; // astrometric offset, not error
+      catalog[0].measure[Nmeas].dbFlags  = 0;
+      catalog[0].measure[Nmeas].averef   = Nave;
+      catalog[0].measure[Nmeas].objID    = catalog[0].average[Nave].objID;
+      catalog[0].measure[Nmeas].catID    = catalog[0].catID;
+
+      catalog[0].measure[Nmeas].t        = (stars[N][0].measure.t == 0) ? TIMEREF : stars[N][0].measure.t;      /** careful : time_t vs e_time **/
+      catalog[0].measure[Nmeas].t_msec   = (stars[N][0].measure.t == 0) ?       0 : stars[N][0].measure.t_msec; /** careful : time_t vs e_time **/
+
+      stars[N][0].found = Nmeas;
+      next_meas[Nmeas] = -1;
+      Nmeas ++;
+    }
+    Nave ++;
+  }
+      
+  REALLOCATE (catalog[0].average, Average, Nave);
+  REALLOCATE (catalog[0].measure, Measure, Nmeas);
+
+  // XXX allow for unsorted output?
+  catalog[0].measure = sort_measure (catalog[0].average, Nave, catalog[0].measure, Nmeas, next_meas);
+
+  /* note stars which have been found in this catalog */
+  for (i = 0; i < Nstars; i++) {
+    if (stars[i][0].found > -1) {
+      stars[i][0].found = -2;
+    } else {
+      stars[i][0].found = -3;
+    }
+  }
+
+  catalog[0].objID    = objID; // new max value, save on catalog close
+  catalog[0].Naverage = Nave;
+  catalog[0].Nmeasure = Nmeas;
+  catalog[0].Nsecf_mem = Nave*Nsecfilt;
+  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmatch);
+
+  free (catalog[0].found);
+  free (X1);
+  free (Y1);
+  free (N1);
+  free (N2);
+  free (X2);
+  free (Y2);
+  free (next_meas);
+
+  return (Nmatch);
+}
Index: /trunk/Ohana/src/addstar/src/getusnob.c
===================================================================
--- /trunk/Ohana/src/addstar/src/getusnob.c	(revision 26279)
+++ /trunk/Ohana/src/addstar/src/getusnob.c	(revision 26280)
@@ -178,6 +178,5 @@
     free (buffer);
     Nitemsum += Nitems;
-    fprintf (stderr, "Nusno: %d, Nitems: %d, Nitemsum : %d, Nbytes stars: %d\n",
-	     Nusno, Nitems, Nitemsum, (unsigned int) (Nusno*sizeof(Stars)));
+    if (VERBOSE) fprintf (stderr, "Nusno: %d, Nitems: %d, Nitemsum : %d, Nbytes stars: %d\n", Nusno, Nitems, Nitemsum, (unsigned int) (Nusno*sizeof(Stars)));
     fclose (f);
   }