Index: /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/include/dvomerge.h
===================================================================
--- /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/include/dvomerge.h	(revision 39219)
+++ /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/include/dvomerge.h	(revision 39220)
@@ -164,4 +164,9 @@
 StarPar   *sort_starpar     	  PROTO((Average *average, off_t Naverage, StarPar *starpar, off_t Nstarpar, off_t *next));
 
+off_t 	  *build_galphot_links    PROTO((Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot));
+off_t 	  *init_galphot_links     PROTO((Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot));
+int   	   add_galp_link     	  PROTO((Average *average, off_t *next, off_t Ngalphot, off_t NGALPHOT));
+GalPhot   *sort_galphot     	  PROTO((Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot, off_t *next));
+
 off_t 	  *init_missing_links     PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing));
 int   	   add_miss_link     	  PROTO((Average *average, off_t *next, off_t Nmissing));
Index: /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/LoadCatalog.c
===================================================================
--- /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/LoadCatalog.c	(revision 39219)
+++ /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/LoadCatalog.c	(revision 39220)
@@ -1,3 +1,4 @@
 # include "dvomerge.h"
+# define SKIP_GALPHOT 0
 
 int LoadCatalog (Catalog *catalog, SkyRegion *region, char *filename, char *mode, int Nsecfilt) {
@@ -10,6 +11,10 @@
   // always load all of the data (if any exists)
   // XXXX TEMP HACK : skip GALPHOT
-  // XXXX catalog[0].catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_MISSING | DVO_LOAD_SECFILT | DVO_LOAD_LENSING | DVO_LOAD_LENSOBJ | DVO_LOAD_STARPAR | DVO_LOAD_GALPHOT;
+
+# if (SKIP_GALPHOT)  
   catalog[0].catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_MISSING | DVO_LOAD_SECFILT | DVO_LOAD_LENSING | DVO_LOAD_LENSOBJ | DVO_LOAD_STARPAR;
+# else
+  catalog[0].catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_MISSING | DVO_LOAD_SECFILT | DVO_LOAD_LENSING | DVO_LOAD_LENSOBJ | DVO_LOAD_STARPAR | DVO_LOAD_GALPHOT;
+# endif
   
   catalog[0].catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
Index: /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/build_links.c
===================================================================
--- /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/build_links.c	(revision 39219)
+++ /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/build_links.c	(revision 39220)
@@ -24,4 +24,6 @@
 */
 
+/*** Measure ****************************************************************************************/
+
 /* build the initial links assuming the table is sorted, 
    not partial, and has a correct set of average[].measureOffset,Nmeasure values */
@@ -34,6 +36,10 @@
 
   ALLOCATE (next_meas, off_t, Nmeasure);
-  for (i = 0; i < Naverage; i++, N++) {
+  for (i = 0; i < Naverage; i++) {
+    if (!average[i].Nmeasure) continue;
+    off_t m = average[i].measureOffset;
+    myAssert (measure[m].averef == i, "not sorted");
     for (j = 0; j < average[i].Nmeasure - 1; j++, N++) {
+      myAssert (measure[m+j+1].averef == i, "not sorted");
       next_meas[N] = N + 1;
       if (N >= Nmeasure) {
@@ -50,4 +56,5 @@
       abort ();
     }
+    N++;
   }
   return (next_meas);
@@ -122,5 +129,6 @@
   }
 
-  if (m == -1) {
+  // if Nmeasure is 0, m may have been mis-set; add to the end
+  if ((average[0].Nmeasure == 0) || (m == -1)) {
     average[0].measureOffset = Nmeasure;
   } else {
@@ -136,20 +144,26 @@
 Measure *sort_measure (Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next_meas) {
 
-  off_t i, k, n, N;
+  off_t i, k, n, np, N;
   Measure *tmpmeasure;
 
   /* fix order of Measure (memory intensive, but fast) */
+  np = -1; // previous entry in case of errors
   N = 0; 
   ALLOCATE (tmpmeasure, Measure, Nmeasure);
   for (i = 0; i < Naverage; i++) {
+    if (!average[i].Nmeasure) continue;
     n = average[i].measureOffset;
     average[i].measureOffset = N;
     int myObjID = average[i].objID;
     for (k = 0; k < average[i].Nmeasure; k++, N++) {
-      if (n == -1) abort();
+      if (n == -1) {
+	fprintf (stderr, "entry after %d has a problem\n", (int) np);
+	abort();
+      }
       tmpmeasure[N] = measure[n]; 
       myAssert (measure[n].averef == i, "error in averef?");
       myAssert (measure[n].objID == myObjID, "error in objID?");
       tmpmeasure[N].averef = i;
+      np = n;
       n = next_meas[n];
     }
@@ -159,5 +173,5 @@
 }
 
-/*******************************************************************************************/
+/*** Missing ****************************************************************************************/
 
 /* build the initial links assuming the table is sorted */
@@ -229,5 +243,5 @@
 }
 
-/*******************************************************************************************/
+/*** Lensing ****************************************************************************************/
 
 /* build the initial links assuming the table is sorted, 
@@ -352,12 +366,4 @@
   Lensing *tmplensing;
 
-  /*
-  for (i = 0; i < Naverage; i++) {
-    if (average[i].Nlensing != 4) {
-      fprintf (stderr, "check %d %d %d\n", (int) i, (int) average[i].Nlensing, (int) average[i].lensingOffset);
-    }
-  }
-  */
-
   /* fix order of Lensing (memory intensive, but fast) */
   np = -1;
@@ -368,4 +374,5 @@
     n = average[i].lensingOffset;
     average[i].lensingOffset = N;
+    int myObjID = average[i].objID;
     for (k = 0; k < average[i].Nlensing; k++, N++) {
       if (n == -1) {
@@ -374,5 +381,6 @@
       }
       tmplensing[N] = lensing[n]; 
-      if (lensing[n].averef != i) abort();
+      myAssert (lensing[n].averef == i, "error in averef");
+      myAssert (lensing[n].objID == myObjID, "error in objID?");
       tmplensing[N].averef = i;
       np = n;
@@ -384,5 +392,5 @@
 }
 
-/********  StarPar ********************************************************************************/
+/*** StarPar ********************************************************************************/
 
 /* build the initial links assuming the table is sorted, 
@@ -398,5 +406,8 @@
   for (i = 0; i < Naverage; i++) {
     if (!average[i].Nstarpar) continue;
+    off_t m = average[i].starparOffset;
+    myAssert (starpar[m].averef == i, "not sorted");
     for (j = 0; j < average[i].Nstarpar - 1; j++, N++) {
+      myAssert (starpar[m+j+1].averef == i, "not sorted");
       next_star[N] = N + 1;
       if (N >= Nstarpar) {
@@ -486,5 +497,6 @@
   }
 
-  if (m == -1) {
+  // if Nmeasure is 0, m may have been mis-set; add to the end
+  if ((average[0].Nstarpar == 0) || (m == -1)) {
     average[0].starparOffset = Nstarpar;
   } else {
@@ -503,14 +515,6 @@
   StarPar *tmpstarpar;
 
-  /*
-  for (i = 0; i < Naverage; i++) {
-    if (average[i].Nstarpar != 4) {
-      fprintf (stderr, "check %d %d %d\n", (int) i, (int) average[i].Nstarpar, (int) average[i].starparOffset);
-    }
-  }
-  */
-
   /* fix order of StarPar (memory intensive, but fast) */
-  np = -1;
+  np = -1; // previous entry in case of errors
   N = 0; 
   ALLOCATE (tmpstarpar, StarPar, Nstarpar);
@@ -519,4 +523,5 @@
     n = average[i].starparOffset;
     average[i].starparOffset = N;
+    int myObjID = average[i].objID;
     for (k = 0; k < average[i].Nstarpar; k++, N++) {
       if (n == -1) {
@@ -525,5 +530,6 @@
       }
       tmpstarpar[N] = starpar[n]; 
-      if (starpar[n].averef != i) abort();
+      myAssert (starpar[n].averef == i, "error in averef?");
+      myAssert (starpar[n].objID == myObjID, "error in objID?");
       tmpstarpar[N].averef = i;
       np = n;
@@ -535,3 +541,151 @@
 }
 
-
+/*** GalPhot ****************************************************************************************/
+
+/* build the initial links assuming the table is sorted, 
+   not partial, and has a correct set of average[].galphotOffset,Ngalphot values */
+off_t *init_galphot_links (Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot) {
+
+  off_t i, j, N;
+  off_t *next_galp;
+
+  N = 0;
+
+  ALLOCATE (next_galp, off_t, Ngalphot);
+  for (i = 0; i < Naverage; i++) {
+    if (!average[i].Ngalphot) continue;
+    off_t m = average[i].galphotOffset;
+    myAssert (galphot[m].averef == i, "not sorted");
+    for (j = 0; j < average[i].Ngalphot - 1; j++, N++) {
+      myAssert (galphot[m+j+1].averef == i, "not sorted");
+      next_galp[N] = N + 1;
+      if (N >= Ngalphot) {
+	fprintf (stderr, "WARNING: N out of bounds (1)\n");
+      }
+    }
+    next_galp[N] = -1;
+    if (N >= Ngalphot) {
+      fprintf (stderr, "WARNING: N out of bounds (2)\n");
+    }
+
+    if (N >= Ngalphot) {
+      fprintf (stderr, "overflow in init_galphot_links\n");
+      abort ();
+    }
+    N++;
+  }
+  return (next_galp);
+}
+
+/* construct galphot links which are valid FOR THIS LOAD
+ * - if we have a full load, we will get links which can
+ *   be used by other programs (eg, relphot, etc)
+ * - if we have a partial load, the links are only valid
+ *   for that partial load
+ */ 
+
+off_t *build_galphot_links (Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot) {
+
+  off_t i, m, k, Nm, averef;
+  off_t *next_galp;
+
+  ALLOCATE (next_galp, off_t, Ngalphot);
+
+  /* reset the Nm, offset values for average */
+  for (i = 0; i < Naverage; i++) {
+    average[i].galphotOffset = -1;
+    average[i].Ngalphot     =  0;
+  }
+
+  for (Nm = 0; Nm < Ngalphot; Nm++) {
+    averef = galphot[Nm].averef;
+    m = average[averef].galphotOffset;  
+    next_galp[Nm] = -1;
+
+    if (m == -1) { /* no links yet for source */
+      average[averef].galphotOffset = Nm;
+      average[averef].Ngalphot     = 1;
+      continue;
+    }
+
+    for (k = 0; next_galp[m] != -1; k++) {
+      m = next_galp[m];
+      if (m >= Ngalphot) {
+	fprintf (stderr, "WARNING: m out of bounds (1)\n");
+      }
+    }
+
+    average[averef].Ngalphot = k + 2;
+    next_galp[m] = Nm;
+    if (m >= Ngalphot) {
+      fprintf (stderr, "WARNING: m out of bounds (2)\n");
+    }
+  }
+  return (next_galp);
+}
+
+/* average[].galphotOffset, average[].Ngalphot are valid within an addstar run */
+int add_galp_link (Average *average, off_t *next_galp, off_t Ngalphot, off_t NGALPHOT) {
+
+  off_t k, m;
+
+  /* if we have trouble, check validity of next_galp[m] : m < Ngalphot */
+  m = average[0].galphotOffset;  
+
+  for (k = 0; k < average[0].Ngalphot - 1; k++)  {
+    m = next_galp[m];
+    if (m >= NGALPHOT) {
+      fprintf (stderr, "WARNING: m out of bounds (3)\n");
+    }
+  }
+
+  /* set up references */
+  next_galp[Ngalphot] = -1;
+  if (Ngalphot >= NGALPHOT) {
+    fprintf (stderr, "WARNING: Ngalphot out of bounds (1)\n");
+  }
+
+  // if Ngalphot is 0, m may have been mis-set; add to the end
+  if ((average[0].Ngalphot == 0) || (m == -1)) {
+    average[0].galphotOffset = Ngalphot;
+  } else {
+    next_galp[m] = Ngalphot;
+    if (m >= NGALPHOT) {
+      fprintf (stderr, "WARNING: m out of bounds (4)\n");
+    }
+  }
+
+  return (TRUE);
+}
+
+GalPhot *sort_galphot (Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot, off_t *next_galp) {
+
+  off_t i, k, n, np, N;
+  GalPhot *tmpgalphot;
+
+  /* fix order of GalPhot (memory intensive, but fast) */
+  np = -1;
+  N = 0; 
+  ALLOCATE (tmpgalphot, GalPhot, Ngalphot);
+  for (i = 0; i < Naverage; i++) {
+    if (!average[i].Ngalphot) continue;
+    n = average[i].galphotOffset;
+    average[i].galphotOffset = N;
+    int myObjID = average[i].objID;
+    for (k = 0; k < average[i].Ngalphot; k++, N++) {
+      if (n == -1) {
+	fprintf (stderr, "entry after %d has a problem\n", (int) np);
+	abort();
+      }
+      tmpgalphot[N] = galphot[n]; 
+      myAssert (galphot[n].averef == i, "error in averef?");
+      myAssert (galphot[n].objID == myObjID, "error in objID?");
+      tmpgalphot[N].averef = i;
+      np = n;
+      n = next_galp[n];
+    }
+  }
+  free (galphot);
+  return (tmpgalphot);
+}
+
Index: /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c
===================================================================
--- /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c	(revision 39219)
+++ /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c	(revision 39220)
@@ -234,4 +234,7 @@
     catalog[0].lensing[i].imageID = newID;
   }
+
+  // NOTE galphot also carries an imageID, but it is actually the EXTERNAL image ID reference.  Do NOT update it.
+
   return TRUE;
 }
Index: /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/dvoverify_utils.c
===================================================================
--- /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/dvoverify_utils.c	(revision 39219)
+++ /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/dvoverify_utils.c	(revision 39220)
@@ -350,6 +350,13 @@
       m = catalog.average[i].lensobjOffset;
       for (j = 0; j < catalog.average[i].Nlensobj; j++) {
+	if ((catalog.lensobj[m+j].objID == 0xffffffff) && (catalog.lensobj[m+j].catID == 0xffffffff)) continue;
 	objIDsOK &= (catalog.average[i].objID == catalog.lensobj[m+j].objID);
 	catIDsOK &= (catalog.average[i].catID == catalog.lensobj[m+j].catID);
+	if (!(catalog.average[i].objID == catalog.lensobj[m+j].objID)) {
+	  fprintf (stderr, "!");
+	}
+	if (!(catalog.average[i].catID == catalog.lensobj[m+j].catID)) {
+	  fprintf (stderr, "?");
+	}
       }
     }
Index: /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_new.c
===================================================================
--- /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_new.c	(revision 39219)
+++ /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_new.c	(revision 39220)
@@ -11,5 +11,5 @@
   
   off_t i, j, offset;
-  off_t NAVERAGE, NMEASURE, NLENSING, NSTARPAR, Naverage, Nmeasure, Nlensing, Nstarpar, NsecfiltIn, NsecfiltOut, Nm;
+  off_t NAVERAGE, NMEASURE, NLENSING, NSTARPAR, NGALPHOT, Naverage, Nmeasure, Nlensing, Nstarpar, Ngalphot, NsecfiltIn, NsecfiltOut, Nm;
 
   Naverage = output[0].Naverage;
@@ -17,4 +17,5 @@
   Nlensing = output[0].Nlensing;
   Nstarpar = output[0].Nstarpar;
+  Ngalphot = output[0].Ngalphot;
 
   NsecfiltOut = output[0].Nsecfilt;
@@ -25,4 +26,5 @@
   NLENSING = Nlensing + 1000;
   NSTARPAR = Nstarpar + 1000;
+  NGALPHOT = Ngalphot + 1000;
 
   REALLOCATE (output[0].average, Average, NAVERAGE);
@@ -31,4 +33,5 @@
   REALLOCATE (output[0].lensing, Lensing, NLENSING);
   REALLOCATE (output[0].starpar, StarPar, NSTARPAR);
+  REALLOCATE (output[0].galphot, GalPhot, NGALPHOT);
 
   // add all of these entries to the existing catalog
@@ -100,4 +103,20 @@
       output[0].average[Naverage].Nstarpar = Nm;
 
+      Nm = 0;
+      for (j = 0; j < input[0].average[i].Ngalphot; j++) {
+	  offset = input[0].average[i].galphotOffset + j;
+
+	  output[0].galphot[Ngalphot] = input[0].galphot[offset];
+	  output[0].galphot[Ngalphot].averef = Naverage;
+
+	  Ngalphot ++;
+	  Nm ++;
+	  if (Ngalphot == NGALPHOT) {
+	      NGALPHOT += 1000;
+	      REALLOCATE (output[0].galphot, GalPhot, NGALPHOT);
+	  }
+      }
+      output[0].average[Naverage].Ngalphot = Nm;
+
       Naverage ++;
       if (Naverage == NAVERAGE) {
@@ -111,8 +130,10 @@
   REALLOCATE (output[0].lensing, Lensing, MAX (Nlensing, 1));
   REALLOCATE (output[0].starpar, StarPar, MAX (Nstarpar, 1));
+  REALLOCATE (output[0].galphot, GalPhot, MAX (Ngalphot, 1));
   REALLOCATE (output[0].secfilt, SecFilt, NsecfiltOut*MAX (Naverage, 1));
   output[0].Naverage = Naverage;
   output[0].Nlensing = Nlensing;
   output[0].Nstarpar = Nstarpar;
+  output[0].Ngalphot = Ngalphot;
   output[0].Nsecfilt_mem = Naverage * NsecfiltOut;
 
@@ -124,9 +145,11 @@
   
   if (VERBOSE) {
-      fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures, "OFF_T_FMT" lensing) for catalog\n", 
+      fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures, "OFF_T_FMT" lensing, "OFF_T_FMT" starpar, "OFF_T_FMT" galphot) for catalog\n", 
 	        i, 
 	        output[0].Naverage, 
 	        output[0].Nmeasure, 
-	        output[0].Nlensing);
+	        output[0].Nlensing, 
+	        output[0].Nstarpar, 
+	        output[0].Ngalphot);
   }
   return (TRUE);
Index: /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_old.c
===================================================================
--- /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_old.c	(revision 39219)
+++ /branches/eam_branches/ipp-20151113/Ohana/src/dvomerge/src/merge_catalogs_old.c	(revision 39220)
@@ -14,6 +14,6 @@
   double *X1, *Y1, *X2, *Y2;
   double dX, dY, dR;
-  off_t *N1, *N2, *next_meas, *next_lens, *next_star;
-  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS, Nstar, NSTAR;
+  off_t *N1, *N2, *next_meas, *next_lens, *next_star, *next_galp;
+  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS, Nstar, NSTAR, Ngalp, NGALP;
   int NsecfiltIn;
   int NsecfiltOut;
@@ -54,4 +54,5 @@
   NLENS = Nlens = output[0].Nlensing;
   NSTAR = Nstar = output[0].Nstarpar;
+  NGALP = Ngalp = output[0].Ngalphot;
 
   // current max obj ID for this catalog
@@ -118,8 +119,10 @@
     next_lens = init_lensing_links (output[0].average, Nave, output[0].lensing, Nlens);
     next_star = init_starpar_links (output[0].average, Nave, output[0].starpar, Nstar);
+    next_galp = init_galphot_links (output[0].average, Nave, output[0].galphot, Ngalp);
   } else {
     next_meas = build_measure_links (output[0].average, Nave, output[0].measure, Nmeas);
     next_lens = build_lensing_links (output[0].average, Nave, output[0].lensing, Nlens);
     next_star = build_starpar_links (output[0].average, Nave, output[0].starpar, Nstar);
+    next_galp = build_galphot_links (output[0].average, Nave, output[0].galphot, Ngalp);
   }    
 
@@ -193,4 +196,9 @@
       REALLOCATE (output[0].starpar, StarPar, NSTAR);
     }
+    if (Ngalp + input[0].average[N].Ngalphot >= NGALP) {
+      NGALP = Ngalp + input[0].average[N].Ngalphot + 1000;
+      REALLOCATE (next_galp, off_t, NGALP);
+      REALLOCATE (output[0].galphot, GalPhot, NGALP);
+    }
 
     // 4) average properties from the input and the output db need to be properly merged.
@@ -214,5 +222,5 @@
       if (REPLACE_BY_PHOTCODE) {
 	// index to first measure for this object
-	// XXX this does not support lensing measurements
+	// XXX this does not support lensing, starpar, or galphot measurements
 	if (replace_match (&output[0].average[n], output[0].measure, next_meas, &input[0].average[N], &input[0].measure[offset])) {
 	  input[0].found_t[N] = output[0].average[n].measureOffset;  
@@ -307,4 +315,22 @@
     }
 
+    // if galphot measurements exist, add them too
+    if (output[0].galphot) {
+      for (Nin = 0; Nin < input[0].average[N].Ngalphot; Nin++) {
+	/* add to end of galphot list */
+	add_galp_link (&output[0].average[n], next_galp, Ngalp, NGALP);
+	
+	// set the new galphot
+	off_t galpoff = input[0].average[N].galphotOffset + Nin;
+	output[0].galphot[Ngalp] = input[0].galphot[galpoff];
+
+	output[0].galphot[Ngalp].averef   = n;
+	output[0].galphot[Ngalp].objID    = output[0].average[n].objID;
+	output[0].galphot[Ngalp].catID    = output[0].catID;
+	output[0].average[n].Ngalphot ++;
+	Ngalp ++;
+      }
+    }
+
     // update the average properties to reflect the incoming entries:
     // if the original value is NAN but the input value is not, accept the input:
@@ -357,4 +383,9 @@
       REALLOCATE (output[0].starpar, StarPar, NSTAR);
     }
+    if (Ngalp + input[0].average[N].Ngalphot >= NGALP) {
+      NGALP = Ngalp + input[0].average[N].Ngalphot + 1000;
+      REALLOCATE (next_galp, off_t, NGALP);
+      REALLOCATE (output[0].galphot, GalPhot, NGALP);
+    }
     if (Nave >= NAVE) {
       NAVE = Nave + 1000;
@@ -480,4 +511,30 @@
     }
 
+    /** add galphot for this input average object **/
+    if (output[0].galphot) {
+      output[0].average[Nave].galphotOffset  = Ngalp;
+      for (Nin = 0; Nin < input[0].average[N].Ngalphot; Nin ++) {
+	// supply the galphot values from this detection
+	off_t galpoff = input[0].average[N].galphotOffset + Nin;
+	output[0].galphot[Ngalp]           = input[0].galphot[galpoff];
+
+	// the following galphot elements cannot be set until here:
+	output[0].galphot[Ngalp].averef   = Nave;
+	output[0].galphot[Ngalp].objID    = output[0].average[Nave].objID;
+	output[0].galphot[Ngalp].catID    = output[0].catID;
+
+	// as we add galphot, update Ngalphot to match
+	output[0].average[Nave].Ngalphot ++;
+
+	/* we set next[Ngalp] to -1 here, and update correctly below */
+	next_galp[Ngalp] = -1;
+	Ngalp ++;
+      }
+      int Ngroup = input[0].average[N].Ngalphot;
+      for (j = 0; j < Ngroup - 1; j++) {
+	next_galp[Ngalp - Ngroup + j] = Ngalp - Ngroup + j + 1;
+      }
+    }
+
     Nave ++;
   }
@@ -489,4 +546,5 @@
   REALLOCATE (output[0].lensing, Lensing, Nlens);
   REALLOCATE (output[0].starpar, StarPar, Nstar);
+  REALLOCATE (output[0].galphot, GalPhot, Ngalp);
  
 # define NOSORT 0
@@ -498,4 +556,5 @@
     output[0].lensing = sort_lensing (output[0].average, Nave, output[0].lensing, Nlens, next_lens);
     output[0].starpar = sort_starpar (output[0].average, Nave, output[0].starpar, Nstar, next_star);
+    output[0].galphot = sort_galphot (output[0].average, Nave, output[0].galphot, Ngalp, next_galp);
   }
 
@@ -515,10 +574,12 @@
   output[0].Nlensing = Nlens;
   output[0].Nstarpar = Nstar;
+  output[0].Ngalphot = Ngalp;
   output[0].Nsecfilt_mem = Nave*NsecfiltOut;
-  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nlens: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nlens, Nmatch);
+  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nlens, Ngalp: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nlens, Ngalp, Nmatch);
 
   free (next_meas);
   free (next_lens);
   free (next_star);
+  free (next_galp);
 
   free (X2);