Index: /branches/czw_branch/20160809/Ohana/src/addstar/Makefile
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/addstar/Makefile	(revision 39709)
+++ /branches/czw_branch/20160809/Ohana/src/addstar/Makefile	(revision 39710)
@@ -25,4 +25,5 @@
 load2mass    : $(BIN)/load2mass.$(ARCH)
 loadgalphot  : $(BIN)/loadgalphot.$(ARCH)
+loadgaia     : $(BIN)/loadgaia.$(ARCH)
 loadstarpar  : $(BIN)/loadstarpar.$(ARCH)
 loadstarpar_client : $(BIN)/loadstarpar_client.$(ARCH)
@@ -42,5 +43,5 @@
 # programs in 'SERVER' use the client-server concept and are out of date
 
-INSTALL = addstar addstar_client sedstar loadgalphot loadstarpar loadstarpar_client setobjflags setobjflags_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell load2mass loadwise loadtycho loadbsc loadsupercos 
+INSTALL = addstar addstar_client sedstar loadgalphot loadgaia loadstarpar loadstarpar_client setobjflags setobjflags_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell load2mass loadwise loadtycho loadbsc loadsupercos 
 SERVER  = addstarc addstard addstart
 
@@ -267,8 +268,20 @@
 $(SRC)/psps_ids.$(ARCH).o
 
-# $(SRC)/SkyRegionUtils.$(ARCH).o 
-# $(SRC)/loadstarpar_io.$(ARCH).o 
-# $(SRC)/loadgalphot_remote_hosts.$(ARCH).o 
-# $(SRC)/loadgalphot_save_remote.$(ARCH).o 
+LOAD-GAIA = \
+$(SRC)/loadgaia.$(ARCH).o \
+$(SRC)/ConfigInit.$(ARCH).o \
+$(SRC)/SetSignals.$(ARCH).o \
+$(SRC)/Shutdown.$(ARCH).o \
+$(SRC)/args_loadgaia.$(ARCH).o \
+$(SRC)/find_matches_gaia.$(ARCH).o \
+$(SRC)/loadgaia_catalog.$(ARCH).o \
+$(SRC)/loadgaia_make_subset.$(ARCH).o \
+$(SRC)/loadgaia_readstars.$(ARCH).o \
+$(SRC)/loadgaia_table.$(ARCH).o \
+$(SRC)/resort_catalog.$(ARCH).o \
+$(SRC)/build_links.$(ARCH).o \
+$(SRC)/strhash.$(ARCH).o \
+$(SRC)/sortIDs.$(ARCH).o \
+$(SRC)/psps_ids.$(ARCH).o
 
 LOAD-STARPAR = \
@@ -493,4 +506,5 @@
 $(LOAD-2MASS)  	       : $(INC)/addstar.h $(INC)/2mass.h
 $(LOAD-GALPHOT)	       : $(INC)/addstar.h $(INC)/loadgalphot.h
+$(LOAD-GAIA)	       : $(INC)/addstar.h $(INC)/gaia.h
 $(LOAD-STARPAR)	       : $(INC)/addstar.h $(INC)/loadstarpar.h
 $(LOAD-STARPAR-CLIENT) : $(INC)/addstar.h $(INC)/loadstarpar.h
@@ -514,4 +528,5 @@
 $(BIN)/loadbsc.$(ARCH)        : $(LOAD-BSC)
 $(BIN)/loadgalphot.$(ARCH)    : $(LOAD-GALPHOT)
+$(BIN)/loadgaia.$(ARCH)       : $(LOAD-GAIA)
 $(BIN)/loadstarpar.$(ARCH)    : $(LOAD-STARPAR)
 $(BIN)/loadstarpar_client.$(ARCH) : $(LOAD-STARPAR-CLIENT)
Index: /branches/czw_branch/20160809/Ohana/src/addstar/include/gaia.h
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/addstar/include/gaia.h	(revision 39710)
+++ /branches/czw_branch/20160809/Ohana/src/addstar/include/gaia.h	(revision 39710)
@@ -0,0 +1,44 @@
+
+typedef struct {
+  double R, D;
+  Measure measure;
+  int flag; // in a subset?
+  int found; // assigned to an object?
+} Gaia_Stars;
+
+int   HOST_ID;
+char *HOSTDIR;
+char *CPT_FILE;
+char *INPUT;
+
+AddstarClientOptions args_loadgaia (int *argc, char **argv, AddstarClientOptions options);
+AddstarClientOptions args_loadgaia_client (int *argc, char **argv, AddstarClientOptions options);
+
+int loadgaia_table (SkyList *skylistInput, HostTable *hosts, char *filename, AddstarClientOptions *options);
+
+Gaia_Stars *loadgaia_make_subset (Gaia_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset);
+
+int loadgaia_save_remote (Gaia_Stars *stars, int Nstars, HostTable *hosts, SkyRegion *region, char *fullname, AddstarClientOptions *options);
+
+int save_remote_host (HostInfo *host);
+
+int init_remote_hosts (void);
+void free_remote_hosts (void);
+int find_empty_slot (void);
+int harvest_all (void);
+int harvest_host (void);
+
+int loadgaia_catalog (Gaia_Stars *stars, int Nstars, SkyRegion *region, char *filename, AddstarClientOptions *options);
+
+int galactic_to_celestial (double *R, double *D, double l, double b);
+
+int find_matches_gaia (SkyRegion *region, Gaia_Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions *options);
+
+int loadgaia_save_stars (char *filename, Gaia_Stars *stars, int Nstars);
+Gaia_Stars *loadgaia_load_stars (char *filename, int *nstars);
+
+Gaia_Stars *loadgaia_readstars (char *filename, int *nstars);
+
+int loadgaia_sortStars (Gaia_Stars *stars, int Nstars);
+
+int loadgaia_tmpdir (void);
Index: /branches/czw_branch/20160809/Ohana/src/addstar/src/args_loadgaia.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/addstar/src/args_loadgaia.c	(revision 39710)
+++ /branches/czw_branch/20160809/Ohana/src/addstar/src/args_loadgaia.c	(revision 39710)
@@ -0,0 +1,127 @@
+# include "addstar.h"
+# include "gaia.h"
+
+static void help (void);
+
+AddstarClientOptions args_loadgaia (int *argc, char **argv, AddstarClientOptions options) {
+  
+  int N;
+
+  /* check for help request */
+  if (get_argument (*argc, argv, "-help") ||
+      get_argument (*argc, argv, "-h")) {
+    help ();
+  }
+
+  /* basic mode: image, list, refcat */
+  options.mode = ADDSTAR_MODE_REFCAT;
+
+  /* we do not allow a subset to be extracted -- all or nothing, babe */
+  UserPatch.Rmin = 0;
+  UserPatch.Rmax= 360;
+  UserPatch.Dmin = -90;
+  UserPatch.Dmax = +90;
+
+  // XXX for the moment, make this selection manual.  it needs to be automatic 
+  // based on the state of the SkyTable
+  HOST_ID = 0;
+  PARALLEL = FALSE;
+  if ((N = get_argument (*argc, argv, "-parallel"))) {
+    PARALLEL = TRUE;
+    remove_argument (N, argc, argv);
+  }
+  // this is a test mode : rather than launching the remote jobs and waiting for completion,
+  // relphot will simply list the remote command and wait for the user to signal completion
+  PARALLEL_MANUAL = FALSE;
+  if ((N = get_argument (*argc, argv, "-parallel-manual"))) {
+    PARALLEL = TRUE; // -parallel-manual implies -parallel
+    PARALLEL_MANUAL = TRUE;
+    remove_argument (N, argc, argv);
+  }
+  // this is a test mode : rather than launching the relphot_client jobs remotely, they are 
+  // run in serial via 'system'
+  PARALLEL_SERIAL = FALSE;
+  if ((N = get_argument (*argc, argv, "-parallel-serial"))) {
+    if (PARALLEL_MANUAL) {
+      fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n");
+      exit (1);
+    }
+    PARALLEL = TRUE; // -parallel-serial implies -parallel
+    PARALLEL_SERIAL = TRUE;
+    remove_argument (N, argc, argv);
+  }
+
+  /* only add to existing regions */
+  options.existing_regions = FALSE;
+  if ((N = get_argument (*argc, argv, "-existing-regions"))) {
+    options.existing_regions = TRUE;
+    remove_argument (N, argc, argv);
+  }
+  /* only add to existing objects */
+  options.only_match = FALSE;
+  if ((N = get_argument (*argc, argv, "-only-match"))) {
+    options.only_match = TRUE;
+    remove_argument (N, argc, argv);
+  }
+  /* replace measurement, don't duplicate (ref/cat only) */
+  options.replace = FALSE;
+  if ((N = get_argument (*argc, argv, "-replace"))) {
+    options.replace = TRUE;
+    remove_argument (N, argc, argv);
+  }
+  /* override any header PHOTCODE values */
+  options.photcode = 0;
+  if ((N = get_argument (*argc, argv, "-p"))) {
+    remove_argument (N, argc, argv);
+    options.photcode = GetPhotcodeCodebyName (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+  if ((N = get_argument (*argc, argv, "-photcode"))) {
+    remove_argument (N, argc, argv);
+    options.photcode = GetPhotcodeCodebyName (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
+  /* extra error messages */
+  VERBOSE = FALSE;
+  if ((N = get_argument (*argc, argv, "-v"))) {
+    VERBOSE = TRUE;
+    remove_argument (N, argc, argv);
+  }
+
+  /* other addstar options which cannot be used in loadgaia */
+  options.timeref = 0; 
+  options.mosaic = FALSE;
+  options.skip_missed = FALSE;
+  options.closest = FALSE;
+  options.nosort = FALSE;
+  options.update = FALSE;
+  options.only_images = FALSE;
+  options.calibrate = FALSE;
+  options.quality_airmass = FALSE;
+  ACCEPT_ASTROM = FALSE;
+  FORCE_READ = FALSE;
+  TEXTMODE = FALSE;
+  SUBPIX = FALSE;
+  DUMP = NULL;
+
+  if (*argc < 2) {
+    fprintf (stderr, "USAGE: loadgaia [options] (fitsfile) [..more files]\n");
+    exit (2);
+  }
+  return (options);
+}
+
+static void help () {
+
+  fprintf (stderr, "USAGE: loadgaia [options] (file) [..more files]\n");
+  fprintf (stderr, "  add data from stellar parameter file to DVO\n\n");
+
+  fprintf (stderr, "  optional flags:\n");
+  fprintf (stderr, "  -only-match           	  : only add measurements to existing objects\n");
+  fprintf (stderr, "  -replace              	  : replace time/photcode measurements (no duplication)\n");
+  fprintf (stderr, "  -v                    	  : verbose mode\n");
+  fprintf (stderr, "  -help                 	  : this list\n");
+  fprintf (stderr, "  -h                    	  : this list\n\n");
+  exit (2);
+}
Index: /branches/czw_branch/20160809/Ohana/src/addstar/src/find_matches_gaia.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/addstar/src/find_matches_gaia.c	(revision 39710)
+++ /branches/czw_branch/20160809/Ohana/src/addstar/src/find_matches_gaia.c	(revision 39710)
@@ -0,0 +1,318 @@
+# include "addstar.h"
+# include "gaia.h"
+
+int find_matches_gaia (SkyRegion *region, Gaia_Stars *stars, int NstarsIn, Catalog *tgtcat, AddstarClientOptions *options) {
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+  int Nsec     = GetPhotcodeNsec (options->photcode);
+
+  /** allocate local arrays (stars) **/
+  ALLOCATE_PTR (X1, double, NstarsIn);
+  ALLOCATE_PTR (Y1, double, NstarsIn);
+  ALLOCATE_PTR (N1, off_t,  NstarsIn);
+
+  /** allocate local arrays (tgtcat) **/
+  off_t NAVE = tgtcat[0].Naverage;
+  off_t Nave = tgtcat[0].Naverage;
+  ALLOCATE_PTR (X2, double, NAVE);
+  ALLOCATE_PTR (Y2, double, NAVE);
+  ALLOCATE_PTR (N2, off_t,  NAVE);
+
+  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
+
+  /* internal counters */
+  off_t Nmatch = 0;
+  off_t NMEAS = tgtcat[0].Nmeasure;
+  off_t Nmeas = tgtcat[0].Nmeasure;
+
+  off_t *next_meas = NULL;
+
+  // current max obj ID for this tgtcat
+  unsigned int objID = tgtcat[0].objID;
+  unsigned int catID = tgtcat[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
+   * substantially far from the projection pole. We use the center of the region (tgtcat)
+   * for crval1,2.
+   */
+
+  Coords tcoords;
+  InitCoords (&tcoords, "DEC--ARC");
+  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
+  if ((region[0].Dmax < 90) && (region[0].Dmin > -90)) {
+    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
+  } else {
+    tcoords.crval2 = (region[0].Dmax >= 90) ? 90.0 : -90.0;
+  }
+  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+
+  /* build spatial index (RA sort) referencing input array sequence */
+  off_t Nstars = 0;
+  for (off_t i = 0; i < NstarsIn; i++) {
+    int status = RD_to_XY (&X1[Nstars], &Y1[Nstars], stars[i].R, stars[i].D, &tcoords);
+    if (!status) continue;
+    N1[Nstars] = i;
+    Nstars ++;
+  }
+  if (Nstars < 1) {
+    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", tgtcat[0].filename);
+    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 (off_t i = 0; i < Nave; i++) {
+    RD_to_XY (&X2[i], &Y2[i], tgtcat[0].average[i].R, tgtcat[0].average[i].D, &tcoords);
+    N2[i] = i;
+  }
+  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);
+
+  /* set up pointers for linked list of measure */
+  if (tgtcat[0].sorted && (tgtcat[0].Nmeasure == tgtcat[0].Nmeasure_disk)) {
+    // this version is only valid if we have done a full tgtcat load, and if the tgtcat
+    // is sorted while processed
+    next_meas = init_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas);
+  } else {
+    next_meas = build_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas);
+  }    
+
+  /* choose a radius for matches */
+  double RADIUS = (options->radius == 0) ? 2.0 : options->radius; /* provided by config */
+  double RADIUS2 = RADIUS*RADIUS;
+
+  /****************** find matched stars ********************/
+
+  for (off_t i = 0, 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 */
+    double 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;
+    }
+
+    // skip this star if already assigned to an object
+    if (stars[N1[i]].found) {
+        i++;
+        continue;
+    }
+
+    /* this block will match a given detection to the closest object within range of that detection.
+       XXX note that this matches ALL detections within range of the single object to that same object 
+       this is bad, but I cannot just go in linear order (ie, mark off each object as they are
+       used).  I should make a list of all Nobj * Ndet pairs in range and choose the matches
+       based on their separations.  UGH
+     */
+    
+    /* within match range; look for matches */
+    off_t Jmin = -1;
+    double Rmin = RADIUS2;
+    for (off_t J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) {
+      /* find closest match for this detection */
+      dX = X1[i] - X2[J];
+      double dY = Y1[i] - Y2[J];
+      double 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;
+    }
+
+    /*** a match is found, add to average, measure ***/
+    Nmatch ++;
+    off_t n = N2[Jmin];
+    off_t N = N1[i];
+
+    /* make sure there is space for next entry */
+    if (Nmeas >= NMEAS - 1) {
+      NMEAS = Nmeas + 1000;
+      REALLOCATE (next_meas, off_t, NMEAS);
+      REALLOCATE (tgtcat[0].measure, Measure, NMEAS);
+    }
+
+    /** add measurement for this star **/
+
+    /* add to end of measurement list */
+    add_meas_link (&tgtcat[0].average[n], next_meas, Nmeas, NMEAS);
+
+    // set the new measurements
+    tgtcat[0].measure[Nmeas]          = stars[N].measure;
+
+    // measure now carries R,D (not dR,dD) 
+    tgtcat[0].measure[Nmeas].dbFlags  = 0;
+    tgtcat[0].measure[Nmeas].averef   = n;
+    tgtcat[0].measure[Nmeas].objID    = tgtcat[0].average[n].objID;
+    tgtcat[0].measure[Nmeas].catID    = tgtcat[0].catID;
+    
+    float dRoff = dvoOffsetR(&tgtcat[0].measure[Nmeas], &tgtcat[0].average[n]);
+
+    // rationalize dR
+    if (dRoff > +180.0*3600.0) {
+      // average on high end of boundary, move star up
+      tgtcat[0].measure[Nmeas].R += 360.0;
+      dRoff -= 360.0*3600.0;
+    }
+    if (dRoff < -180.0*3600.0) {
+      // average on low end of boundary, move star down
+      tgtcat[0].measure[Nmeas].R -= 360.0;
+      dRoff += 360.0*3600.0;
+    }
+    if (fabs(dRoff) > 10*RADIUS) {
+      // take declination into account and check again.
+      double cosD = cos(RAD_DEG*tgtcat[0].average[n].D);
+      if (fabs(dRoff*cosD) > 10*RADIUS) {
+	fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n", 
+		 tgtcat[0].average[n].R, tgtcat[0].average[n].D, 
+		 stars[N].R, stars[N].D,
+		 X1[i], X2[Jmin], 
+		 Y1[i], Y2[Jmin]);
+	// XXX abort on this? -- this is a bad failure...
+      }
+    }
+
+    /* Nm is updated, but not written out in -update mode (for existing entries)
+       Nm is recalculated in build_meas_links if loaded table is not sorted */
+    tgtcat[0].average[n].Nmeasure ++;
+    Nmeas ++;
+
+    /* if we choose to flag close encounters, see find_matches.c */
+    /* if we choose to calculate RA,DEC averages, see update_coords.c */
+
+    stars[N].found = TRUE;
+    i++;
+  }
+
+  /*************** add unmatched stars *************************/
+
+  /** incorporate unmatched image stars, if this star is in field of this tgtcat **/
+  /* these new entries are all written out in UPDATE mode */ 
+  for (off_t i = 0; (i < Nstars) && !options->only_match; i++) {
+    /* make sure there is space for next entry */
+    if (Nmeas >= NMEAS - 1) {
+      NMEAS = Nmeas + 1000;
+      REALLOCATE (next_meas, off_t, NMEAS);
+      REALLOCATE (tgtcat[0].measure, Measure, NMEAS);
+    }
+    if (Nave >= NAVE) {
+      NAVE = Nave + 1000;
+      REALLOCATE (tgtcat[0].average, Average, NAVE);
+      if (tgtcat[0].secfilt) {
+	// we only update the secfilt table if it has been allocated for output
+	REALLOCATE (tgtcat[0].secfilt, SecFilt, NAVE*tgtcat[0].Nsecfilt);
+      }
+    }
+
+    if (stars[i].found) continue;
+    if (!IN_REGION (stars[i].R, stars[i].D)) continue;
+
+    dvo_average_init (&tgtcat[0].average[Nave]);
+    tgtcat[0].average[Nave].R         	   = stars[i].R;
+    tgtcat[0].average[Nave].D         	   = stars[i].D;
+
+    tgtcat[0].average[Nave].Nmeasure  	   = 1;
+    tgtcat[0].average[Nave].measureOffset  = Nmeas;
+    tgtcat[0].average[Nave].objID     	   = objID;
+    tgtcat[0].average[Nave].catID     	   = catID;
+
+    if (PSPS_ID) {
+        tgtcat[0].average[Nave].extID = CreatePSPSObjectID(tgtcat[0].average[Nave].R, tgtcat[0].average[Nave].D);
+    }
+
+    objID ++;
+
+    // we only update the secfilt table if it has been allocated for output
+    for (int j = 0; tgtcat[0].secfilt && (j < Nsecfilt); j++) {
+      dvo_secfilt_init (&tgtcat[0].secfilt[Nave*Nsecfilt+j], SECFILT_RESET_ALL);
+    }
+
+    // supply the measurements from this detection
+    dvo_measure_init (&tgtcat[0].measure[Nmeas]);
+    tgtcat[0].measure[Nmeas] = stars[i].measure;
+
+    // the following measure elements cannot be set until here:
+    tgtcat[0].measure[Nmeas].dbFlags  = 0;
+    tgtcat[0].measure[Nmeas].averef   = Nave;
+    tgtcat[0].measure[Nmeas].objID    = tgtcat[0].average[Nave].objID;
+    tgtcat[0].measure[Nmeas].catID    = tgtcat[0].catID;
+
+    /* set the average magnitude if not already set and the photcode.equiv is not 0 */
+    /* in UPDATE mode, this value is not saved; use relphot to recalculate */
+    if (Nsec > -1) { 
+      tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF);
+    }
+
+    Nmeas ++;
+
+    // update the next_meas pointer for this entry (last one for this star is -1)
+    next_meas[Nmeas-1] = -1;
+    stars[i].found = TRUE;
+    Nave ++;
+  }
+
+  REALLOCATE (tgtcat[0].average, Average, Nave);
+  REALLOCATE (tgtcat[0].measure, Measure, Nmeas);
+ 
+  if (options->nosort) {
+    tgtcat[0].sorted = FALSE;
+  } else {
+    tgtcat[0].sorted = TRUE;
+    tgtcat[0].measure = sort_measure (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas, next_meas);
+  }
+
+  /* check if the tgtcat has changed?  if no change, no need to write */
+  tgtcat[0].objID    = objID; // new max value, save on tgtcat close
+  tgtcat[0].Naverage = Nave;
+  tgtcat[0].Nmeasure = Nmeas;
+  tgtcat[0].Nsecfilt_mem = tgtcat[0].secfilt ? Nave*Nsecfilt : 0;
+  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas, Nmatch);
+
+  free (X1);
+  free (Y1);
+  free (N1);
+  free (X2);
+  free (Y2);
+  free (N2);
+  free (next_meas);
+
+  return (Nmatch);
+}
+
+/* 
+   notes:
+   
+   for finding if a tgtcat star is in an image or an image star is in the tgtcat:
+   
+   tgtcats have boundaries defined by RA and DEC, but they may curve in projection
+   images have boundaries which are lines in pixels coords, but curve in RA and DEC
+   
+   tgtcat[0].found[Ncat] but stars[Nstar].found
+   
+*/
Index: /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia.c	(revision 39710)
+++ /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia.c	(revision 39710)
@@ -0,0 +1,36 @@
+# include "addstar.h"
+# include "gaia.h"
+
+/* This is the DVO program to upload gaia detections from D. Finkbeiner into a DVO database.
+   It is modeled on the loadtycho program but the source files are FITS tables.
+
+   USAGE: loadgaia -D CATDIR (catdir) (gaiafile) [...more files]
+
+*/
+
+int main (int argc, char **argv) {
+
+  int i;
+  SkyTable *sky;
+  SkyList *skylist = NULL;
+  AddstarClientOptions options;
+
+  // need to construct these options with args_loadtycho...
+  SetSignals ();
+  options = ConfigInit (&argc, argv);
+  options = args_loadgaia (&argc, argv, options);
+
+  // load the full sky description table:
+  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
+  SkyTableSetFilenames (sky, CATDIR, "cpt");
+  
+  // generate the subset matching the user-selected region
+  skylist = SkyListByPatch (sky, -1, &UserPatch);
+
+  for (i = 1; i < argc; i++) {
+    fprintf (stderr, "loading %s\n", argv[i]);
+    loadgaia_table (skylist, NULL, argv[i], &options);
+  }
+  exit (0);
+}  
+
Index: /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_catalog.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_catalog.c	(revision 39710)
+++ /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_catalog.c	(revision 39710)
@@ -0,0 +1,32 @@
+# include "addstar.h"
+# include "gaia.h"
+
+int loadgaia_catalog (Gaia_Stars *stars, int Nstars, SkyRegion *region, char *filename, AddstarClientOptions *options) {
+
+  Catalog catalog;
+
+  // now we have all of the loaded stars in this catalog
+  dvo_catalog_init (&catalog, TRUE);
+  catalog.filename = filename;
+  catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
+  catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
+  catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_SECFILT | DVO_LOAD_MEASURE;
+  catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
+    
+  // an error exit status here is a significant error
+  if (!dvo_catalog_open (&catalog, region, VERBOSE, "w")) {
+    fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename);
+    exit (2);
+  }
+
+  find_matches_gaia (region, stars, Nstars, &catalog, options);
+    
+  SetProtect (TRUE);
+  if (!dvo_catalog_save (&catalog, VERBOSE)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog.filename); exit (1); }
+  if (!dvo_catalog_unlock (&catalog)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog.filename); exit (1); }
+  SetProtect (FALSE);
+
+  dvo_catalog_free (&catalog);
+
+  return (TRUE);
+}
Index: /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_make_subset.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_make_subset.c	(revision 39710)
+++ /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_make_subset.c	(revision 39710)
@@ -0,0 +1,44 @@
+# include "addstar.h"
+# include "gaia.h"
+
+// assign stars in the given region to the subset (NOTE: stars are sorted by RA, start is
+// first entry in stars array in this region)
+
+Gaia_Stars *loadgaia_make_subset (Gaia_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset) {
+
+  int i;
+
+  Gaia_Stars *subset = NULL;
+
+  // collect array of (Stars *) stars in a new output catalog
+  int Nsubset = 0;
+  int NSUBSET = 3000;
+  ALLOCATE (subset, Gaia_Stars, NSUBSET);
+
+  // find the rest of the stars in this output region
+  for (i = start; i < Nstars; i++) {
+    if (stars[i].flag) continue;
+
+    // check if in skyregion
+    if (stars[i].R <  region[0].Rmin) continue;
+    if (stars[i].R >= region[0].Rmax) break;
+    if (stars[i].D <  region[0].Dmin) continue;
+    if (stars[i].D >= region[0].Dmax) continue;
+	  
+    // check if in UserPatch (a GLOBAL)
+    if (stars[i].R < UserPatch.Rmin) continue;
+    if (stars[i].R > UserPatch.Rmax) break;
+    if (stars[i].D < UserPatch.Dmin) continue;
+    if (stars[i].D > UserPatch.Dmax) continue;
+	  
+    subset[Nsubset] = stars[i];
+    Nsubset ++;
+
+    stars[i].flag = TRUE;
+
+    CHECK_REALLOCATE (subset, Gaia_Stars, NSUBSET, Nsubset, 10000);
+  }
+
+  *nsubset = Nsubset;
+  return subset;
+}
Index: /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_readstars.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_readstars.c	(revision 39710)
+++ /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_readstars.c	(revision 39710)
@@ -0,0 +1,136 @@
+# include "addstar.h"
+# include "gaia.h"
+
+# define GET_COLUMN(OUT,NAME,TYPE) \
+  TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
+  myAssert (!strcmp(type, #TYPE), "wrong column type");
+
+Gaia_Stars *loadgaia_readstars (char *filename, int *nstars) {
+
+  // read in the full FITS files
+  FILE *f = fopen (filename, "r");
+  if (f == NULL) Shutdown ("can't read gaia file: %s", filename);
+
+  int i, Ncol;
+  off_t Nrow;
+
+  Header header;
+  Matrix matrix;
+  Header theader;
+  FTable ftable;
+  
+  // load in PHU segment (ignore)
+  if (!gfits_fread_header (f, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read image subset header\n");
+    fclose (f);
+    return NULL;
+  }
+  if (!gfits_fread_matrix (f, &matrix, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n");
+    gfits_free_header (&header);
+    fclose (f);
+    return NULL;
+  }
+
+  ftable.header = &theader;
+
+  // load data for this header 
+  if (!gfits_load_header (f, &theader)) {
+    fclose (f);
+    return NULL;
+  }
+  if (!gfits_fread_ftable_data (f, &ftable, FALSE)) {
+    fclose (f);
+    return (NULL);
+  }
+
+  char type[16];
+
+  GET_COLUMN (RA,      "RA",              	 double);
+  GET_COLUMN (DEC,     "DEC",             	 double);
+  GET_COLUMN (dRA,     "RA_ERROR",        	 float);
+  GET_COLUMN (dDEC,    "DEC_ERROR",       	 float);
+  GET_COLUMN (gMag,    "PHOT_G_MEAN_MAG", 	 float);
+  GET_COLUMN (dgFlux,  "PHOT_G_MEAN_FLUX_ERROR", float);
+  GET_COLUMN (Nobs,    "PHOT_G_N_OBS",           float);
+
+  // free the memory associated with the FITS files
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  gfits_free_header (&theader);
+  gfits_free_table (&ftable);
+
+  fclose (f);
+
+  int NstarsIn = matrix.Naxis[2];
+
+  double Rmin = +360.0;
+  double Rmax = -360.0;
+  double Dmin = +360.0;
+  double Dmax = -360.0;
+
+  int Nstars = 0;
+  int NSTARS = 0.1*NstarsIn;
+
+  Gaia_Stars *stars = NULL;
+  ALLOCATE (stars, Gaia_Stars, NSTARS);
+
+  for (i = 0; i < NstarsIn; i++) {
+
+    Rmin = MIN (Rmin, RA[i]);
+    Rmax = MAX (Rmax, RA[i]);
+    Dmin = MIN (Dmin, DEC[i]);
+    Dmax = MAX (Dmax, DEC[i]);
+
+    float flux = pow(10.0, -0.4*(gMag[i] - 25.524770));
+
+    stars[Nstars].R = RA[i];
+    stars[Nstars].D = DEC[i];
+    stars[Nstars].flag  = FALSE;
+    stars[Nstars].found = FALSE;
+
+    dvo_measure_init (&stars[Nstars].measure);
+
+    stars[Nstars].measure.R = RA[i];
+    stars[Nstars].measure.D = DEC[i];
+    stars[Nstars].measure.dXccd = (int)(0x7fff*MAX(0.000,MIN(0.999, ((log10(dRA[i]) + 3.0)/6.0)))); 
+    stars[Nstars].measure.dYccd = (int)(0x7fff*MAX(0.000,MIN(0.999, ((log10(dDEC[i]) + 3.0)/6.0)))); 
+    // dXccd,dYccd range from 0 (10^-3) mas to 0x7fff (10^3 mas)
+    stars[Nstars].measure.M = gMag[i];
+    stars[Nstars].measure.dM = dgFlux[i] / flux;
+
+    stars[Nstars].measure.FluxPSF = flux;
+    stars[Nstars].measure.dFluxPSF = dgFlux[i];
+
+    Nstars ++;
+
+    CHECK_REALLOCATE (stars, Gaia_Stars, NSTARS, Nstars, 10000);
+  }
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+
+  free (RA);
+  free (DEC);
+  free (dRA);
+  free (dDEC);
+  free (gMag);
+  free (dgFlux);
+  free (Nobs);
+
+  *nstars = Nstars;
+  return (stars);
+}
+
+int loadgaia_sortStars (Gaia_Stars *stars, int Nstars) {
+
+# define SWAPFUNC(A,B){ Gaia_Stars temp = stars[A]; stars[A] = stars[B]; stars[B] = temp; }
+# define COMPARE(A,B)(stars[A].R < stars[B].R)
+
+  OHANA_SORT (Nstars, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+  
+  return TRUE;
+}
+
Index: /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_table.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_table.c	(revision 39710)
+++ /branches/czw_branch/20160809/Ohana/src/addstar/src/loadgaia_table.c	(revision 39710)
@@ -0,0 +1,51 @@
+# include "addstar.h"
+# include "gaia.h"
+
+int loadgaia_table (SkyList *skylistInput, HostTable *hosts, char *filename, AddstarClientOptions *options) {
+  OHANA_UNUSED_PARAM(hosts);
+  
+  int i, Nstars;
+
+  Gaia_Stars *stars = loadgaia_readstars (filename, &Nstars);
+  if (!stars) return FALSE;
+
+  // sort the stars by RA
+  loadgaia_sortStars (stars, Nstars);
+
+  // scan through the stars, loading the containing catalogs
+  // skip through table for unsaved stars
+  for (i = 0; i < Nstars; i++) {
+    if (stars[i].flag) continue;
+
+    // scan forward until we read the UserPatch
+    if (stars[i].R < UserPatch.Rmin) continue;
+    if (stars[i].R > UserPatch.Rmax) break;
+    if (stars[i].D < UserPatch.Dmin) continue;
+    if (stars[i].D > UserPatch.Dmax) continue;
+
+    // identify the relevant catalog
+    SkyList *skylist = SkyRegionByPoint_List (skylistInput, -1, stars[i].R, stars[i].D);
+    if (skylist[0].Nregions == 0) {
+      SkyListFree (skylist);
+      continue;
+    }
+    SkyRegion *region = skylist[0].regions[0];
+
+    // select stars matching this region
+    int Nsubset;
+    Gaia_Stars *subset = loadgaia_make_subset (stars, Nstars, i, region, &Nsubset);
+
+    // In parallel mode, write out the subset to a disk file.  Block until a remote host
+    // is available.  In serial mode, just match against the appropriate region and save
+    // NOTE: disable parallel mode for now: 
+    // loadgaia_save_remote (subset, Nsubset, hosts, region, skylist[0].filename[0], options);
+    loadgaia_catalog (subset, Nsubset, region, skylist[0].filename[0], options);
+  }
+
+  // wait for last remote clients to finish
+  // NOTE: disable parallel mode for now: 
+  // harvest_all ();
+
+  return TRUE;
+}
+