Index: /branches/eam_branches/ipp-20150616/Ohana/src/addstar/Makefile
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/Makefile	(revision 38482)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/Makefile	(revision 38483)
@@ -24,4 +24,5 @@
 sedstar      : $(BIN)/sedstar.$(ARCH)
 load2mass    : $(BIN)/load2mass.$(ARCH)
+loadgalphot  : $(BIN)/loadgalphot.$(ARCH)
 loadstarpar  : $(BIN)/loadstarpar.$(ARCH)
 loadstarpar_client : $(BIN)/loadstarpar_client.$(ARCH)
@@ -36,7 +37,7 @@
 mkcmf        : $(BIN)/mkcmf.$(ARCH)
 
-all: addstar addstar_client sedstar load2mass loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf loadwise loadtycho loadsupercos dumpskycells findskycell
-
-INSTALL = addstar addstar_client sedstar load2mass loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf loadwise loadtycho loadsupercos dumpskycells findskycell
+all: addstar addstar_client sedstar load2mass loadgalphot loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf loadwise loadtycho loadsupercos dumpskycells findskycell
+
+INSTALL = addstar addstar_client sedstar load2mass loadgalphot loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf loadwise loadtycho loadsupercos dumpskycells findskycell
 
 # I need to fix the client/server version of addstar now that I have dropped Stars
@@ -244,4 +245,23 @@
 $(SRC)/psps_ids.$(ARCH).o
 
+LOAD-GALPHOT = \
+$(SRC)/loadgalphot.$(ARCH).o \
+$(SRC)/ConfigInit.$(ARCH).o \
+$(SRC)/SetSignals.$(ARCH).o \
+$(SRC)/Shutdown.$(ARCH).o \
+$(SRC)/args_loadgalphot.$(ARCH).o \
+$(SRC)/find_matches_galphot.$(ARCH).o \
+$(SRC)/loadgalphot_catalog.$(ARCH).o \
+$(SRC)/loadgalphot_make_subset.$(ARCH).o \
+$(SRC)/loadgalphot_readstars.$(ARCH).o \
+$(SRC)/loadgalphot_fit2d.$(ARCH).o \
+$(SRC)/loadgalphot_table.$(ARCH).o \
+$(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-STARPAR = \
 $(SRC)/loadstarpar.$(ARCH).o \
@@ -404,4 +424,5 @@
 $(FINDSKYCELL)	       : $(INC)/addstar.h
 $(LOAD-2MASS)  	       : $(INC)/addstar.h $(INC)/2mass.h
+$(LOAD-GALPHOT)	       : $(INC)/addstar.h $(INC)/loadgalphot.h
 $(LOAD-STARPAR)	       : $(INC)/addstar.h $(INC)/loadstarpar.h
 $(LOAD-STARPAR-CLIENT) : $(INC)/addstar.h $(INC)/loadstarpar.h
@@ -421,4 +442,5 @@
 $(BIN)/loadwise.$(ARCH)       : $(LOAD-WISE)
 $(BIN)/loadtycho.$(ARCH)      : $(LOAD-TYCHO)
+$(BIN)/loadgalphot.$(ARCH)    : $(LOAD-GALPHOT)
 $(BIN)/loadstarpar.$(ARCH)    : $(LOAD-STARPAR)
 $(BIN)/loadstarpar_client.$(ARCH) : $(LOAD-STARPAR-CLIENT)
Index: /branches/eam_branches/ipp-20150616/Ohana/src/addstar/include/loadgalphot.h
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/include/loadgalphot.h	(revision 38483)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/include/loadgalphot.h	(revision 38483)
@@ -0,0 +1,49 @@
+
+typedef struct {
+  double **s;
+  double **b;
+  double **c;
+  double *Cii;
+  int order;
+  int nterm;
+  int wterm;
+  int mterm;
+  float mean;
+  float sigma;
+  double c00;
+  double c10;
+  double c20;
+  double c01;
+  double c11;
+  double c02;
+} Fit2D;
+
+typedef struct {
+  double R, D;
+  GalPhot galphot;
+  int flag; // in a subset?
+  int found; // assigned to an object?
+} GalPhot_Stars;
+
+AddstarClientOptions args_loadgalphot (int argc, char **argv, AddstarClientOptions options);
+
+int loadgalphot_table (SkyList *skylistInput, HostTable *hosts, char *filename, AddstarClientOptions *options);
+
+GalPhot_Stars *loadgalphot_readstars (char *filename, int *nstars);
+
+int loadgalphot_sortStars (GalPhot_Stars *stars, int Nstars);
+
+GalPhot_Stars *loadgalphot_make_subset (GalPhot_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset);
+
+int loadgalphot_catalog (GalPhot_Stars *stars, int Nstars, SkyRegion *region, char *filename, AddstarClientOptions *options);
+
+int find_matches_galphot (SkyRegion *region, GalPhot_Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions *options);
+
+int FitChisqMinimum (Fit2D *fit, GalPhot *galphot, 
+		     float *chisq, float *flux, float *fluxErr, int Npts, 
+		     float MajorMin, float MajorMax, float MajorDel, 
+		     float MinorMin, float MinorMax, float MinorDel);
+
+Fit2D *fit2d_init (int order);
+void fit2d_free (Fit2D *fit);
+int fit2d (Fit2D *fit, float *xval, float *yval, float *zval, float *zfit, char *mask, int Npts);
Index: /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/args_loadgalphot.c
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/args_loadgalphot.c	(revision 38483)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/args_loadgalphot.c	(revision 38483)
@@ -0,0 +1,116 @@
+# include "addstar.h"
+# include "loadgalphot.h"
+
+static void help (void);
+
+AddstarClientOptions args_loadgalphot (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);
+  }
+
+  /* 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 loadgalphot */
+  options.photcode = 0;
+  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: loadgalphot [options] (fitsfile) [..more files]\n");
+    exit (2);
+  }
+  return (options);
+}
+
+static void help () {
+
+  fprintf (stderr, "USAGE: loadgalphot [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/eam_branches/ipp-20150616/Ohana/src/addstar/src/find_matches_galphot.c
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/find_matches_galphot.c	(revision 38483)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/find_matches_galphot.c	(revision 38483)
@@ -0,0 +1,218 @@
+# include "addstar.h"
+# include "loadgalphot.h"
+
+// we find matching stars by (R,D).  we add the galphot to the array of galphots with appropriate links.  
+int find_matches_galphot (SkyRegion *region, GalPhot_Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions *options) {
+  
+  off_t i, j, n, N, J, Jmin;
+  double RADIUS, RADIUS2, Rmin;
+  double *X1, *Y1, *X2, *Y2;
+  double dX, dY, dR;
+  off_t *N1, *N2;
+  off_t Nave, NAVE, Ngalphot, NGALPHOT, Nmatch;
+
+  /** allocate local arrays (stars) **/
+  ALLOCATE (X1, double, Nstars);
+  ALLOCATE (Y1, double, Nstars);
+  ALLOCATE (N1, off_t,  Nstars);
+
+  /** allocate local arrays (catalog) **/
+  NAVE = Nave = catalog[0].Naverage;
+  ALLOCATE (X2, double, NAVE);
+  ALLOCATE (Y2, double, NAVE);
+  ALLOCATE (N2, off_t,  NAVE);
+
+  /* internal counters */
+  Nmatch = 0;
+  NGALPHOT = Ngalphot = catalog[0].Ngalphot;
+
+  // current max obj ID for this catalog
+  unsigned int objID = catalog[0].objID;
+  unsigned int catID = catalog[0].catID;
+
+  int Nsecfilt = catalog[0].Nsecfilt;
+
+  /* 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.
+   */
+  Coords tcoords;
+  InitCoords (&tcoords, "DEC--ARC");
+  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
+  tcoords.crval2 = (region[0].Dmax < 90) ? 0.5*(region[0].Dmin + region[0].Dmax) : 90.0;
+  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+
+  /* build spatial index (RA sort) referencing input array sequence */
+  for (i = 0; i < Nstars; i++) {
+    RD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords);
+    N1[i] = i;
+  }
+  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;
+  }
+  sort_coords_index (X2, Y2, N2, Nave);
+
+  /* choose a radius for matches */
+  RADIUS = options->radius; /* provided by config */
+  RADIUS2 = RADIUS*RADIUS;
+
+  /** find matched stars **/
+  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;
+    }
+
+    if (stars[N1[i]].found) {
+      /* this star has already been assigned to an object in this or another catalog */
+      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 */
+    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;
+    }
+
+    /*** a match is found, add to galphot with links to average ***/
+    Nmatch ++;
+    n = N2[Jmin];
+    N = N1[i];
+
+    /* make sure there is space for next entry */
+    if (Ngalphot >= NGALPHOT) {
+      NGALPHOT = Ngalphot + 1000;
+      REALLOCATE (catalog[0].galphot, GalPhot, NGALPHOT);
+    }
+
+    /** add galphot for this star **/
+
+    // set the new galphotments
+    catalog[0].galphot[Ngalphot]          = stars[N].galphot;
+    catalog[0].galphot[Ngalphot].averef   = n;
+    catalog[0].galphot[Ngalphot].objID    = catalog[0].average[n].objID;
+    catalog[0].galphot[Ngalphot].catID    = catID;
+
+    // NOTE: include R,D in galphot?
+    // XXX if I add R,D to galphot, I should rationalize R,D to the same boundary?
+
+    stars[N].found = TRUE;
+    catalog[0].average[n].Ngalphot ++;
+    Ngalphot ++;
+    i++;
+  }
+
+# if (1)
+  /** incorporate unmatched image stars? **/
+  for (i = 0; (i < Nstars) && !options->only_match; i++) {
+
+    // skip already matched stars
+    if (stars[i].found) continue;
+    if (!IN_REGION (stars[i].R, stars[i].D)) continue;
+
+    /* make sure there is space for next entry */
+    if (Ngalphot >= NGALPHOT) {
+      NGALPHOT = Ngalphot + 1000;
+      REALLOCATE (catalog[0].galphot, GalPhot, NGALPHOT);
+    }
+    if (Nave >= NAVE) {
+      NAVE = Nave + 1000;
+      REALLOCATE (catalog[0].average, Average, NAVE);
+      REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*Nsecfilt);
+    }
+
+    dvo_average_init (&catalog[0].average[Nave]);
+    catalog[0].average[Nave].R         	   = stars[i].R;
+    catalog[0].average[Nave].D         	   = stars[i].D;
+
+    catalog[0].average[Nave].Ngalphot  	   = 1;
+    catalog[0].average[Nave].galphotOffset = Ngalphot;
+    catalog[0].average[Nave].objID     	   = objID;
+    catalog[0].average[Nave].catID     	   = catID;
+
+    if (PSPS_ID) {
+      catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D);
+    }
+
+    objID ++;
+
+    for (j = 0; j < Nsecfilt; j++) {
+      dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j], SECFILT_RESET_ALL);
+    }
+
+    catalog[0].galphot[Ngalphot]        = stars[i].galphot;
+    catalog[0].galphot[Ngalphot].averef = Nave;
+    catalog[0].galphot[Ngalphot].objID  = catalog[0].average[Nave].objID;
+    catalog[0].galphot[Ngalphot].catID  = catID;
+
+    stars[i].found = TRUE;
+    Ngalphot ++;
+    Nave ++;
+  }
+# endif      
+
+  REALLOCATE (catalog[0].average, Average, Nave);
+  REALLOCATE (catalog[0].galphot, GalPhot, Ngalphot);
+ 
+  catalog[0].sorted = FALSE;
+
+  /* check if the catalog has changed?  if no change, no need to write */
+  catalog[0].objID     = objID; // new max value, save on catalog close
+  catalog[0].Naverage  = Nave;
+  catalog[0].Ngalphot  = Ngalphot;
+  catalog[0].Nsecfilt_mem = Nave*Nsecfilt;
+  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Ngalphot: %d "OFF_T_FMT" "OFF_T_FMT" ("OFF_T_FMT" matches)\n",  Nstars, Nave, Ngalphot, Nmatch);
+
+  free (X1);
+  free (Y1);
+  free (N1);
+  free (N2);
+  free (X2);
+  free (Y2);
+
+  return (Nmatch);
+}
Index: /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot.c
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot.c	(revision 38483)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot.c	(revision 38483)
@@ -0,0 +1,30 @@
+# include "addstar.h"
+# include "loadgalphot.h"
+
+/* This is the DVO program to upload fullforce summary cmf data for galaxies into a DVO database.
+
+   USAGE: loadgalphot -D CATDIR (catdir) (file.cmf)
+*/
+
+int main (int argc, char **argv) {
+
+  int i;
+  AddstarClientOptions options;
+
+  // need to construct these options with args_loadgalphot...
+  options = ConfigInit (&argc, argv);
+  options = args_loadgalphot (argc, argv, options);
+
+  // load the full sky description table (dvodb must exist)
+  SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE);
+  SkyTableSetFilenames (sky, CATDIR, "cpt");
+  
+  // generate the subset matching the user-selected region
+  SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch);
+
+  for (i = 1; i < argc; i++) {
+    fprintf (stderr, "loading %s\n", argv[i]);
+    loadgalphot_table (skylist, NULL, argv[i], &options);
+  }
+  exit (0);
+}  
Index: /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_catalog.c
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_catalog.c	(revision 38483)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_catalog.c	(revision 38483)
@@ -0,0 +1,28 @@
+# include "addstar.h"
+# include "loadgalphot.h"
+
+int loadgalphot_catalog (GalPhot_Stars *stars, int Nstars, SkyRegion *region, char *filename, AddstarClientOptions *options) {
+
+  Catalog catalog;
+
+  // now we have all of the loaded stars in this catalog
+  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_GALPHOT;
+  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_galphot (region, stars, Nstars, &catalog, options);
+    
+  dvo_catalog_save (&catalog, VERBOSE);
+  dvo_catalog_unlock (&catalog);
+  dvo_catalog_free (&catalog);
+
+  return (TRUE);
+}
Index: /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_fit2d.c
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_fit2d.c	(revision 38483)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_fit2d.c	(revision 38483)
@@ -0,0 +1,217 @@
+# include "addstar.h"
+# include "loadgalphot.h"
+
+Fit2D *fit2d_init (int order) {
+
+  int i;
+
+  Fit2D *fit;
+  ALLOCATE (fit, Fit2D, 1);
+
+  // allocate static arrays
+  fit->nterm = order + 1;
+  fit->wterm = fit->nterm*(fit->nterm + 1)/2;
+  fit->mterm = 2*order + 1;
+
+  /* allocate the summation matrices */
+  ALLOCATE (fit->s, double *, fit->mterm);
+  ALLOCATE (fit->b, double *, fit->wterm);
+  ALLOCATE (fit->c, double *, fit->wterm);
+  for (i = 0; i < fit->wterm; i++) {
+    ALLOCATE (fit->c[i], double, fit->wterm);
+    ALLOCATE (fit->b[i], double, 1);
+  }
+  for (i = 0; i < fit->mterm; i++) {
+    ALLOCATE (fit->s[i], double, fit->mterm);
+  }
+  return fit;
+}
+
+void fit2d_free (Fit2D *fit) {
+
+  int i;
+
+  for (i = 0; i < fit->wterm; i++) {
+    free (fit->c[i]);
+    free (fit->b[i]);
+  }
+  free (fit->b);
+  free (fit->c);
+
+  for (i = 0; i < fit->mterm; i++) {
+    free (fit->s[i]);
+  }
+  free (fit->s);
+  free (fit);
+}
+
+// unweighted 2D 2nd order fit (6 free parameters)
+int fit2d (Fit2D *fit, float *xval, float *yval, float *zval, float *zfit, char *mask, int Npts) {
+  
+  int N, i, j, nx, ny, K, k, n, Npt;
+
+  double mean, sigma, maxsigma;
+
+  float *x, *y, *z, *zf;
+
+  float ClipNSigma = 0;
+  int ClipNiter  = 1;
+
+  for (N = 0; N < ClipNiter; N++) {
+    /* init registers for current pass */
+    for (i = 0; i < fit->wterm; i++) {
+      memset (fit->c[i], 0, fit->wterm*sizeof(double));
+      memset (fit->b[i], 0, sizeof(double));
+    }
+    for (i = 0; i < fit->mterm; i++) {
+      memset (fit->s[i], 0, fit->mterm*sizeof(double));
+    }
+
+    // pointers which loop over datapoints
+    x = xval;
+    y = yval;
+    z = zval;
+
+    /* add up the x,y values */
+    for (i = 0; i < Npts; i++, x++, y++) {
+      if (mask[i]) continue;
+      if (!finite(*x) || !finite(*y) || !finite(*z)) continue;
+      double Y = 1.0;
+      for (ny = 0; ny < fit->mterm; ny++) {
+	double X = Y;
+	for (nx = 0; nx < fit->mterm - ny; nx++) {
+	  fit->s[nx][ny] += X;
+	  X = X * (*x);
+	}
+	Y = Y * (*y);
+      }
+    }
+
+    // pointers which loop over datapoints
+    x = xval;
+    y = yval;
+    z = zval;
+
+    /* add up the z values */
+    for (i = 0; i < Npts; i++, x++, y++, z++) {
+      if (mask[i]) continue;
+      if (!finite(*x) || !finite(*y) || !finite(*z)) continue;
+      double Y = *z;
+      for (j = 0, ny = 0; ny < fit->nterm; ny++) {
+	double X = Y;
+	for (nx = 0; nx < fit->nterm - ny; nx++, j++) {
+	  fit->b[j][0] += X;
+	  X = X * (*x);
+	}
+	Y = Y * (*y);
+      }
+    }
+
+    /* re-sort mterm x mterm matrix to wterm matrix */
+    for (k = j = 0; j < fit->nterm; j++) {
+      for (i = 0; i < fit->nterm - j; i++, k++) {
+	for (K = ny = 0; ny < fit->nterm; ny++) {
+	  for (nx = 0; nx < fit->nterm - ny; nx++, K++) {
+	    fit->c[K][k] = fit->s[nx+i][ny+j];
+	  }
+	}
+      }
+    }
+
+    // invert the c,b matrix equation
+    dgaussjordan (fit->c, fit->b, fit->wterm, 1);
+
+    /* the b[][0] terms are in the following order:
+       y^0 x^0, y^0 x^1, ... y^0 x^N
+       y^1 x^0, y^1 x^1, ... y^1 x^N
+       ...
+       y^N x^0, y^N x^1, ... y^N x^N
+    */
+
+    /* generate fitted values */
+    x  = xval;
+    y  = yval;
+    zf = zfit;
+    for (n = 0; n < Npts; n++, x++, y++, zf++) {
+      if (!finite(*x) || !finite(*y) || !finite(*z)) continue;
+      *zf = 0;
+      double X = 1;
+      for (i = ny = 0; ny < fit->nterm; ny++) {
+	double Y = X;
+	for (nx = 0; nx < fit->nterm - ny; nx++, i++) {
+	  *zf += fit->b[i][0]*Y;
+	  Y = Y * (*y);
+	}
+	X = X * (*x);
+      }
+    }
+
+    /* measure fit residual scatter */
+    x  = xval;
+    y  = yval;
+    z  = zval;
+    zf = zfit;
+    float dZ  = 0.0;
+    float dZ2 = 0.0;
+    for (i = Npt = 0; i < Npts; i++, x++, y++, z++, zf++) {
+      if (mask[i]) continue;
+      if (!finite(*x) || !finite(*y) || !finite(*z)) continue;
+      float dZi = (*z - *zf);
+      dZ  += dZi;
+      dZ2 += SQ(dZi);
+      Npt ++;
+    }
+    mean  = dZ / Npt;
+    sigma = sqrt (fabs(dZ2/Npt - SQ(mean)));
+    maxsigma = ClipNSigma * sigma;
+
+    // if (VERBOSE) gprint (GP_ERR, "mean: %g, sigma: %g, maxsigma: %g\n", mean, sigma, maxsigma);
+
+    /* mask outlier points */
+    x  = xval;
+    y  = yval;
+    z  = zval;
+    zf = zfit;
+    int Nmask = 0;
+    for (i = 0; ClipNSigma && (i < Npts); i++, x++, y++, z++, zf++) {
+      float dZi = (*z - *zf);
+      if (fabs(dZi) > maxsigma) {
+	mask[i] = TRUE;
+	Nmask ++;
+      } else {
+	mask[i] = FALSE;
+      }	
+    }
+    // if (VERBOSE) gprint (GP_ERR, "pass: %d, Nmask: %d\n", N, Nmask);
+  }
+    
+  // save the fit
+  for (i = 0; i < fit->wterm; i++) {
+    fit->Cii[i] = fit->b[i][0];
+  }
+  switch (fit->order) {
+    case 0:
+      fit->c00 = fit->Cii[0];
+      break;
+    case 1:
+      fit->c00 = fit->Cii[0];
+      fit->c10 = fit->Cii[1];
+      fit->c01 = fit->Cii[2];
+      break;
+    case 2:
+      fit->c00 = fit->Cii[0];
+      fit->c10 = fit->Cii[1];
+      fit->c20 = fit->Cii[2];
+      fit->c01 = fit->Cii[3];
+      fit->c11 = fit->Cii[4];
+      fit->c02 = fit->Cii[5];
+      break;
+    default:
+      myAbort("invalid order");
+  }
+
+  fit->mean  = mean;
+  fit->sigma = sigma;
+
+  return (TRUE);
+}
Index: /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_make_subset.c
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_make_subset.c	(revision 38483)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_make_subset.c	(revision 38483)
@@ -0,0 +1,44 @@
+# include "addstar.h"
+# include "loadgalphot.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)
+
+GalPhot_Stars *loadgalphot_make_subset (GalPhot_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset) {
+
+  int i;
+
+  GalPhot_Stars *subset = NULL;
+
+  // collect array of (Stars *) stars in a new output catalog
+  int Nsubset = 0;
+  int NSUBSET = 3000;
+  ALLOCATE (subset, GalPhot_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, GalPhot_Stars, NSUBSET, Nsubset, 10000);
+  }
+
+  *nsubset = Nsubset;
+  return subset;
+}
Index: /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_readstars.c
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_readstars.c	(revision 38483)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_readstars.c	(revision 38483)
@@ -0,0 +1,275 @@
+# include "addstar.h"
+# include "loadgalphot.h"
+
+# define GET_COLUMN(OUT,NAME,TYPE) \
+  TYPE *OUT = gfits_get_bintable_column_data (&header, &ftable, NAME, type, &NrowNew, &Ncol); \
+  myAssert (!strcmp(type, #TYPE), "wrong column type"); \
+  if (firstCol) { Nrow = NrowNew; firstCol = FALSE; } \
+  else { myAssert (Nrow == NrowNew, "table length mismatch"); }
+
+static float *Xpt = NULL;
+static float *Ypt = NULL;
+static float *Rpt = NULL;
+static float *Rsr = NULL;
+static char *mask = NULL;
+
+static float *chisqFit = NULL;
+
+GalPhot_Stars *loadgalphot_readstars (char *filename, int *nstars) {
+
+  // read in the full FITS files ('cause I don't have a partial read option)
+  FILE *f = fopen (filename, "r");
+  if (f == NULL) Shutdown ("can't read stellar parameter file: %s", filename);
+
+  int i, Ncol, Nsample;
+  off_t Nrow = 0;
+  off_t NrowNew;
+
+  Header PHU;
+  Header header;
+  FTable ftable;
+  
+  // load in PHU for astrometry
+  if (!gfits_fread_header (f, &PHU)) {
+    if (VERBOSE) fprintf (stderr, "can't read image subset header\n");
+    fclose (f);
+    return NULL;
+  }
+
+  if (!gfits_find_Xheader (f, &header, "SkyChip.xgal")) {
+    if (VERBOSE) fprintf (stderr, "can't read galaxy photometry header\n");
+    gfits_free_header (&header);
+    gfits_free_header (&PHU);
+    fclose (f);
+    return NULL;
+  }
+
+  ftable.header = &header;
+
+  if (!gfits_fread_ftable_data (f, &ftable, FALSE)) {
+    if (VERBOSE) fprintf (stderr, "can't read galaxy photometry data\n");
+    gfits_free_header (&header);
+    gfits_free_header (&PHU);
+    fclose (f);
+    return (NULL);
+  }
+
+  char type[16];
+  int firstCol = TRUE;
+
+  // GET_COLUMN (ID,           "IPP_IDET",     int);
+  // GET_COLUMN (MODEL_TYPE,   "MODEL_TYPE",   int);
+  GET_COLUMN (X_FIT,        "X_FIT", 	    float);
+  GET_COLUMN (Y_FIT,        "Y_FIT", 	    float);
+  // GET_COLUMN (NPIX,         "NPIX",  	    float);
+  GET_COLUMN (GAL_FLUX,     "GAL_FLUX",     float); Nsample = Ncol;
+  GET_COLUMN (GAL_FLUX_ERR, "GAL_FLUX_ERR", float); myAssert (Ncol == Nsample, "invalid table");
+  GET_COLUMN (GAL_CHISQ,    "GAL_CHISQ",    float); myAssert (Ncol == Nsample, "invalid table");
+  GET_COLUMN (FR_MAJOR_MIN, "FR_MAJOR_MIN", float);
+  GET_COLUMN (FR_MAJOR_MAX, "FR_MAJOR_MAX", float);
+  GET_COLUMN (FR_MAJOR_DEL, "FR_MAJOR_DEL", float);
+  GET_COLUMN (FR_MINOR_MIN, "FR_MINOR_MIN", float);
+  GET_COLUMN (FR_MINOR_MAX, "FR_MINOR_MAX", float);
+  GET_COLUMN (FR_MINOR_DEL, "FR_MINOR_DEL", float);
+
+  // free the memory associated with the FITS files
+  gfits_free_header (&header);
+  gfits_free_table (&ftable);
+
+  GalPhot_Stars *stars = NULL;
+  ALLOCATE (stars, GalPhot_Stars, Nrow);
+
+  Coords coords;
+  GetCoords (&coords, &PHU);
+
+  // we reuse these arrays : allocate just once:
+  ALLOCATE (Xpt, float, Nsample);
+  ALLOCATE (Ypt, float, Nsample);
+  ALLOCATE (Rpt, float, Nsample);
+  ALLOCATE (Rsr, float, Nsample);
+  ALLOCATE (mask, char, Nsample);
+  ALLOCATE (chisqFit, float, Nsample);
+
+  Fit2D *fit = fit2d_init (2);
+
+  for (i = 0; i < Nrow; i++) {
+
+    double R, D;
+    XY_to_RD (&R, &D, X_FIT[i], Y_FIT[i], &coords);
+
+    stars[i].R = R;
+    stars[i].D = D;
+    stars[i].flag  = FALSE;
+    stars[i].found = FALSE;
+
+    // stars[i].galphot.R = R;
+    // stars[i].galphot.D = D;
+    stars[i].galphot.Xfit = X_FIT[i];
+    stars[i].galphot.Yfit = Y_FIT[i];
+
+    // I need to match up with the entries in xfit to get:
+    // theta, theta_err, index, 
+    // note that FR_MAJOR, etc are the fractions of the stack fit value
+    stars[i].galphot.theta = NAN;
+    stars[i].galphot.thetaErr = NAN;
+
+    // I have a grid of measurements with (Flux, dFlux, Chisq) at each point
+    // I need to find the minimum position (interpolated) in this 2D space
+
+    FitChisqMinimum (fit, &stars[i].galphot, &GAL_CHISQ[Nsample*i], &GAL_FLUX[Nsample*i], &GAL_FLUX_ERR[Nsample*i], Nsample, FR_MAJOR_MIN[i], FR_MAJOR_MAX[i], FR_MAJOR_DEL[i], FR_MINOR_MIN[i], FR_MINOR_MAX[i], FR_MINOR_DEL[i]);
+  }
+
+  free (Xpt);
+  free (Ypt);
+  free (Rpt);
+  free (Rsr);
+  free (mask);
+  free (chisqFit);
+
+  fit2d_free (fit);
+
+  *nstars = Nrow;
+  return (stars);
+}
+
+int loadgalphot_sortStars (GalPhot_Stars *stars, int Nstars) {
+
+# define SWAPFUNC(A,B){ GalPhot_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;
+}
+
+int FitChisqMinimum (Fit2D *fit, GalPhot *galphot, 
+		     float *chisq, float *flux, float *fluxErr, int Npts, 
+		     float MajorMin, float MajorMax, float MajorDel, 
+		     float MinorMin, float MinorMax, float MinorDel) {
+
+  int i;
+
+  // validate the square dimensions??
+
+  int Nx = (MajorMax - MajorMin) / MajorDel;
+  int Ny = (MinorMax - MinorMin) / MinorDel;
+  myAssert (Nx*Ny == Npts, "inconsistent grid");
+  
+  // here are the steps:
+  // find the minimum chisq grid point (set X,Y values)
+  // choose the nearest 9 points
+  // fit 2D parabola (6 free parameters)
+  // use this to outlier-clip remaining points
+  // fit 2D parabola to the rest
+
+  memset (mask, 0, Npts*sizeof(char));
+
+  float chisqMin = FLT_MAX;
+  int iMin = -1;
+  int Nvalid = 0;
+  for (i = 0; i < Npts; i++) {
+    if (!isfinite(chisq[i])) {
+      mask[i] = TRUE;
+      continue;
+    }
+    int iY = i / Npts;
+    int iX = i % Npts;
+    Xpt[i] = MajorMin + MajorDel*iX;
+    Ypt[i] = MinorMin + MinorDel*iY;
+    if (chisq[i] < chisqMin) {
+      chisqMin = chisq[i];
+      iMin = i;
+      Nvalid ++;
+    }
+  }
+  if (!Nvalid) return FALSE; // we cannot do anything if there is no finite chisq
+  
+  if (Nvalid < 6) {
+    galphot->majorAxis = Xpt[iMin];
+    galphot->minorAxis = Ypt[iMin];
+    galphot->majorAxisErr = MajorDel;
+    galphot->minorAxisErr = MinorDel;
+    return TRUE;
+  }
+
+  // find the list of distances from the min point
+  for (i = 0; i < Npts; i++) {
+    Rpt[i] = hypot(Xpt[i] - Xpt[iMin], Ypt[i] - Ypt[iMin]);
+    // only include unmasked points
+    if (mask[i]) {
+      Rsr[i] = FLT_MAX;
+    } else {
+      Rsr[i] = Rpt[i];
+    }
+  }
+
+  // sorts the distances in ascending order
+  fsort (Rsr, Npts);
+  int Nmin = MIN(Nvalid, 9);
+
+  float Rmax = Rsr[Nmin-1];
+  for (i = 0; i < Npts; i++) {
+    if (Rpt[i] > Rmax) mask[i] = TRUE;
+  }
+
+  // fit the inner 9 points
+  fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts);
+
+  // exclude any extreme outliers
+  float maxSigma = 5.0*fit->sigma;
+
+  memset (mask, 0, Npts*sizeof(char));
+  for (i = 0; i < Npts; i++) {
+    if (!isfinite(chisq[i])) { mask[i] = TRUE; continue; }
+    float dZ = chisq[i] - chisqFit[i];
+    if (fabs(dZ) > maxSigma) { mask[i] = TRUE; continue; }
+  }
+  
+  // re-fit all except the most extreme set
+  fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts);
+
+  // get X,Y for the chisq min from the fit:
+  float Det = 1.0 / (4*fit->c20*fit->c02 - SQ(fit->c11));
+  float Xmin = (fit->c01*fit->c11 - 2.0*fit->c02*fit->c10) * Det;
+  float Ymin = (fit->c10*fit->c11 - 2.0*fit->c20*fit->c01) * Det;
+
+  // chisqMin @ (Xmin,Ymin)
+  chisqMin = 
+    fit->c00 + 
+    fit->c10*Xmin + fit->c20*Xmin*Xmin +
+    fit->c01*Ymin + fit->c02*Ymin*Ymin +
+    fit->c11*Xmin*Ymin;
+
+  float chisqOff = chisqMin + 1.0;
+
+  float A, B, C, Xoff, Yoff;
+
+  // Xoff @ chisqMin + 1.0
+  A = fit->c20;
+  B = fit->c10 + fit->c11*Ymin;
+  C = fit->c00 + fit->c01*Ymin + fit->c02*Ymin*Ymin - chisqOff;
+  Xoff = (-B + sqrt(B*B - 4*A*C)) / (2.0 * A);
+
+  // Xoff @ chisqMin + 1.0
+  A = fit->c02;
+  B = fit->c01 + fit->c11*Xmin;
+  C = fit->c00 + fit->c10*Xmin + fit->c20*Xmin*Xmin - chisqOff;
+  Yoff = (-B + sqrt(B*B - 4*A*C)) / (2.0 * A);
+
+  galphot->majorAxis = Xmin;
+  galphot->minorAxis = Ymin;
+  galphot->majorAxisErr = fabs(Xoff - Xmin);
+  galphot->minorAxisErr = fabs(Yoff - Ymin);
+
+  // use bilinear interpolation to choose Flux @ (Xmin,Ymin)?
+  galphot->mag    = -2.5*log10(flux[iMin]); // correct for exptime?
+  galphot->magErr = -2.5*log10(fluxErr[iMin]); // correct for exptime?
+  galphot->chisq  = chisqMin;
+
+  
+  return TRUE;
+}
+
Index: /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_table.c
===================================================================
--- /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_table.c	(revision 38483)
+++ /branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_table.c	(revision 38483)
@@ -0,0 +1,49 @@
+# include "addstar.h"
+# include "loadgalphot.h"
+
+int loadgalphot_table (SkyList *skylistInput, HostTable *hosts, char *filename, AddstarClientOptions *options) {
+  
+  int i, Nstars;
+
+  GalPhot_Stars *stars = loadgalphot_readstars (filename, &Nstars);
+
+  // sort the stars by RA
+  loadgalphot_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;
+    GalPhot_Stars *subset = loadgalphot_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: 
+    // loadgalphot_save_remote (subset, Nsubset, hosts, region, skylist[0].filename[0], options);
+    loadgalphot_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;
+}
+