Index: trunk/Ohana/src/relphot/Makefile
===================================================================
--- trunk/Ohana/src/relphot/Makefile	(revision 34088)
+++ trunk/Ohana/src/relphot/Makefile	(revision 34260)
@@ -49,4 +49,5 @@
 $(SRC)/setExclusions.$(ARCH).o 	 \
 $(SRC)/setMrelFinal.$(ARCH).o 	 \
+$(SRC)/BoundaryTreeOps.$(ARCH).o 	 \
 $(SRC)/write_coords.$(ARCH).o
 
@@ -77,4 +78,5 @@
 $(SRC)/setExclusions.$(ARCH).o 	 \
 $(SRC)/setMrelFinal.$(ARCH).o    \
+$(SRC)/BoundaryTreeOps.$(ARCH).o 	 \
 $(SRC)/write_coords.$(ARCH).o
 
Index: trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- trunk/Ohana/src/relphot/include/relphot.h	(revision 34088)
+++ trunk/Ohana/src/relphot/include/relphot.h	(revision 34260)
@@ -111,4 +111,6 @@
 char        *BCATALOG;
 ModeType     MODE;
+
+char        *BOUNDARY_TREE;
 
 double MAG_LIM;
@@ -242,4 +244,5 @@
 float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
 short         getUbercalDist      PROTO((off_t meas, int cat));
+float         getCenterOffset     PROTO((off_t meas, int cat, Measure *measure, unsigned int *myID));
 Image        *getimage            PROTO((off_t N));
 Image        *getimages           PROTO((off_t *N, off_t **LineNumber));
@@ -340,2 +343,6 @@
 int client_logger_message (char *format,...);
 
+int MatchImageName (off_t meas, int cat, char *name);
+
+int load_tree (char *treefile);
+int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec);
Index: trunk/Ohana/src/relphot/src/BoundaryTreeOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/BoundaryTreeOps.c	(revision 34260)
+++ trunk/Ohana/src/relphot/src/BoundaryTreeOps.c	(revision 34260)
@@ -0,0 +1,36 @@
+# include "relphot.h"
+
+// XXX for the moment, only load one boundary tree at a time
+// XXX in fact, only allow RINGS.V3...
+
+static BoundaryTree *tree = NULL;
+
+int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec) {
+
+  int zone, band;
+
+  if (!primaryCellName) return FALSE;
+
+  primaryCellName[0] = 0;
+
+  if (!tree) return FALSE;
+
+  if (!BoundaryTreeCellCoords (tree, &zone, &band, ra, dec)) {
+    fprintf (stderr, "mismatch!\n");
+    return FALSE;
+  }
+  
+  snprintf (primaryCellName, DVO_MAX_PATH, "RINGS.V3.%s", tree->name[zone][band]);
+  return TRUE;
+}
+
+int load_tree (char *treefile) {
+
+  tree = BoundaryTreeLoad (treefile);
+  if (!tree) {
+    fprintf (stderr, "failed to load boundary tree %s\n", treefile);
+    exit (2);
+  }
+
+  return TRUE;
+}
Index: trunk/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/ImageOps.c	(revision 34088)
+++ trunk/Ohana/src/relphot/src/ImageOps.c	(revision 34260)
@@ -365,4 +365,45 @@
   distance = image[i].ubercalDist; // was dummy3 in structure
   return (distance);
+}
+
+// returns image.Mcal - ff(x,y)
+float getCenterOffset (off_t meas, int cat, Measure *measure, unsigned int *myID) {
+
+  off_t i;
+  float distance;
+
+  if (!MeasureToImage) return -1;
+
+  i = MeasureToImage[cat][meas];
+  if (i == -1) return (1000);
+
+  float Xcenter = 0.5*image[i].NX;
+  float Ycenter = 0.5*image[i].NY;
+
+  *myID = image[i].imageID;
+
+  distance = hypot (measure[0].Xccd - Xcenter, measure[0].Yccd - Ycenter);
+  return (distance);
+}
+
+// returns image.Mcal - ff(x,y)
+int MatchImageName (off_t meas, int cat, char *name) {
+
+  off_t i;
+
+  if (!name) return FALSE;
+  if (!name[0]) return FALSE;
+
+  if (!MeasureToImage) return FALSE;
+
+  i = MeasureToImage[cat][meas];
+  if (i == -1) return FALSE;
+
+  // this is a bit crude: stack image names are of the form 
+  // RINGS.V3.skycell.1495.027.sky.191211.stk.988232.cmf 
+  // the primaryCell has a name of the form RINGS.V3.skycell.1495
+
+  if (!strncmp(image[i].name, name, strlen(name))) return TRUE;
+  return FALSE;
 }
 
Index: trunk/Ohana/src/relphot/src/StarOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/StarOps.c	(revision 34088)
+++ trunk/Ohana/src/relphot/src/StarOps.c	(revision 34260)
@@ -15,5 +15,6 @@
   double *wlist;
   double *aplist;
-  double *daplist;
+  double *kronlist;
+  double *dkronlist;
 } SetMrelInfo;
 
@@ -164,5 +165,6 @@
   SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist
   ALLOCATE (results.aplist, double, Nmax);
-  ALLOCATE (results.daplist, double, Nmax);
+  ALLOCATE (results.kronlist, double, Nmax);
+  ALLOCATE (results.dkronlist, double, Nmax);
 
   for (i = 0; i < Ncatalog; i++) {
@@ -174,5 +176,6 @@
   SetMrelInfoFree (&results);
   free (results.aplist);
-  free (results.daplist);
+  free (results.kronlist);
+  free (results.dkronlist);
   return (TRUE);
 }
@@ -303,21 +306,26 @@
 int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
 
-  off_t j, k, m;
+  off_t j, k, m, ID;
   int N;
   float Msys, Mcal, Mmos, Mgrid;
 
-  StatType stats, apstats;
+  StatType stats, apstats, kronstats;
   liststats_setmode (&stats, STATMODE);
   liststats_setmode (&apstats, STATMODE);
-
-  double *list    = results->list;
-  double *dlist   = results->dlist;
-  double *wlist   = results->wlist;
-  double *aplist  = results->aplist;
-  double *daplist = results->daplist;
+  liststats_setmode (&kronstats, STATMODE);
+
+  double *list      = results->list;
+  double *dlist     = results->dlist;
+  double *wlist     = results->wlist;
+  double *aplist    = results->aplist;
+  double *kronlist  = results->kronlist;
+  double *dkronlist = results->dkronlist;
 
   SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
 
   int isSetMrelFinal = (pass >= 0);
+
+  char *primaryCell = NULL;
+  ALLOCATE (primaryCell, char, DVO_MAX_PATH);
 
   for (j = 0; j < catalog[Nc].Naverage; j++) {
@@ -325,8 +333,10 @@
 
     // option for a test print
-    if (FALSE && (catalog[Nc].average[j].objID == 0x46a4) && (catalog[Nc].average[j].catID == 0xf40e)) {
+    if (FALSE && (catalog[Nc].average[j].objID == 0x7146) && (catalog[Nc].average[j].catID == 0x49d8)) {
       fprintf (stderr, "test obj\n");
       print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure);
     }
+
+    BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
 
     int GoodPS1 = FALSE;
@@ -355,4 +365,15 @@
       int haveSynth = FALSE;
       int haveStack = FALSE;
+
+      // need to find the measurement closest to the center of its skycell, as well as the
+      // closest for the subset of primary projection cells
+
+      float stackCenterOffsetMin = 1e9;
+      int stackCenterIDmin = -1;
+      off_t stackCenterMeasureMin = -1;
+
+      float stackPrimaryOffsetMin = 1e9;
+      int stackPrimaryIDmin = -1;
+      off_t stackPrimaryMeasureMin = -1;
 
       int forceSynth = FALSE;
@@ -401,4 +422,8 @@
 	  float Map = PhotAper (&catalog[Nc].measure[m]);
 	  aplist[N] = Map - Mcal - Mmos - Mgrid;
+
+	  float Mkron = PhotKron (&catalog[Nc].measure[m]);
+	  kronlist[N] = Mkron - Mcal - Mmos - Mgrid;
+	  dkronlist[N] = catalog[Nc].measure[m].dMkron;
 
 	  // special options for PS1 data
@@ -416,6 +441,27 @@
 	  // gpc1 stack data
 	  if ((catalog[Nc].measure[m].photcode >= 11000) && (catalog[Nc].measure[m].photcode <= 11400)) {
-	    if (pass < 2) continue;
+	    // if (pass < 2) continue;
 	    haveStack = TRUE;
+
+	    unsigned int stackImageID;
+
+	    // which stack image should we use for the mean value?
+	    // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell
+	    if (MatchImageName (m, Nc, primaryCell)) {
+	      float stackPrimaryOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID);
+	      if (stackPrimaryOffset < stackPrimaryOffsetMin) {
+		stackPrimaryOffsetMin = stackPrimaryOffset;
+		stackPrimaryIDmin = stackImageID;
+		stackPrimaryMeasureMin = m;
+	      }
+	    }
+
+	    // get the center distance for the generic case:
+	    float stackCenterOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID);
+	    if (stackCenterOffset < stackCenterOffsetMin) {
+	      stackCenterOffsetMin = stackCenterOffset;
+	      stackCenterIDmin = stackImageID;
+	      stackCenterMeasureMin = m;
+	    }
 	  }
 
@@ -549,7 +595,58 @@
 
 	// NOTE : use the modified weight for apmags as well as psf mags
-	liststats (aplist, daplist, wlist, N, &apstats);
-
+	liststats (aplist, dlist, wlist, N, &apstats);
 	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map  = apstats.mean; 
+
+	liststats (kronlist, dkronlist, wlist, N, &kronstats);
+	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mkron  = kronstats.mean; 
+	catalog[Nc].secfilt[Nsecfilt*j+Nsec].dMkron = kronstats.error; 
+
+	if (haveStack) {
+	  m  = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;
+	  ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
+
+	  // get the zero point for the selected image
+	  float zp = Mcal + Mmos + Mgrid + PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]);
+
+	  // flux_cgs : erg sec^1 cm^-2 Hz^-1
+	  // mag_inst : -2.5 log (cts/sec)
+	  // mag_inst : -2.5 log (flux_inst)
+	  // flux_inst = ten(-0.4*mag_inst)
+
+	  // mag_AB = -2.5 log (flux_cgs) - 48.6 (~by definition) [~Vega flux in V-band]
+	  // flux_cgs = ten(-0.4*(mag_AB + 48.6))
+
+	  // flux_AB : ten(-0.4*mag_AB)
+
+	  // flux_cgs = ten(-0.4*48.6) * flux_AB
+	  // flux_AB  = ten(+0.4*48.6) * flux_cgs
+	    
+	  // flux_Jy : flux_cgs * 10^23
+
+	  // flux_AB = ten(+0.4*48.6) * ten(-23) * flux_Jy
+
+	  // mag_AB = mag_inst + ZP
+
+	  // flux_inst = ten(-0.4*(mag_AB - ZP)) = ten(0.4*ZP) * flux_AB
+
+	  // flux_AB = flux_inst * ten(-0.4*ZP)
+
+	  // flux_inst * ten(-0.4*ZP) = ten(+0.4*48.6 - 23) * flux_Jy
+
+	  // flux_inst = flux_Jy * ten(0.4*ZP + 0.4*48.6 - 23)
+	  // flux_inst = flux_Jy * ten(0.4*ZP - 3.56)
+	  // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56)
+
+	  // zpFactor to go from instrumental flux to Janskies
+	  float zpFactor = pow(10.0, -0.4*zp + 3.56);
+
+	  // need to put in AB mag factor to get to Janskies (or uJy?)
+	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxPSF   = zpFactor * catalog[Nc].measure[m].FluxPSF;  
+	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxPSF  = zpFactor * catalog[Nc].measure[m].dFluxPSF; 
+	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxKron  = zpFactor * catalog[Nc].measure[m].FluxKron; 
+	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron;
+
+	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackID   = ID;
+	}
 
 	// NOTE: for 2MASS measurements, Next should be 1, as should N
@@ -607,4 +704,5 @@
     }
   }
+  if (primaryCell) free (primaryCell);
   return (TRUE);
 }
Index: trunk/Ohana/src/relphot/src/args.c
===================================================================
--- trunk/Ohana/src/relphot/src/args.c	(revision 34088)
+++ trunk/Ohana/src/relphot/src/args.c	(revision 34260)
@@ -197,4 +197,10 @@
   }
 
+  if ((N = get_argument (argc, argv, "-boundary-tree"))) {
+    remove_argument (N, &argc, argv);
+    load_tree (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   SHOW_PARAMS = FALSE;
   if ((N = get_argument (argc, argv, "-params"))) {