Index: trunk/Ohana/src/addstar/src/findskycell.c
===================================================================
--- trunk/Ohana/src/addstar/src/findskycell.c	(revision 34843)
+++ trunk/Ohana/src/addstar/src/findskycell.c	(revision 34844)
@@ -15,5 +15,5 @@
 
 void usage (void) {
-  fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir)\n");
+  fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir) [-nx Nx] [-ny Ny]\n");
   fprintf (stderr, "USAGE: findcell -tree (tree) (datafile)\n");
   fprintf (stderr, "   (datafile) should contain a list of RA,DEC pairs\n");
@@ -24,4 +24,8 @@
 int apply_tree (char *treefile, char *datafile);
 
+float SCALE = 1.0;
+int NX_SUB = 1;
+int NY_SUB = 1;
+
 int main (int argc, char **argv) {
 
@@ -37,4 +41,24 @@
   if (get_argument (argc, argv, "-help")) usage ();
   if (get_argument (argc, argv, "-h")) usage ();
+
+  NX_SUB = NY_SUB = 1;
+  if ((N = get_argument (argc, argv, "-nx"))) {
+    remove_argument (N, &argc, argv);
+    NX_SUB = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-ny"))) {
+    remove_argument (N, &argc, argv);
+    NY_SUB = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  /* pixel scale (arcsec/pixel) */
+  SCALE = 1.0;
+  if ((N = get_argument (argc, argv, "-scale"))) {
+    remove_argument (N, &argc, argv);
+    SCALE = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
 
   /* extra error messages */
@@ -109,4 +133,7 @@
 
   tree.Nzone = 46;
+  tree.NX_SUB = NX_SUB;
+  tree.NY_SUB = NY_SUB;
+  tree.dPix = SCALE/3600.0;
 
   ALLOCATE (tree.Nband, int, tree.Nzone);
@@ -122,4 +149,8 @@
   ALLOCATE (tree.ra,   double *, tree.Nzone);
   ALLOCATE (tree.dec,  double *, tree.Nzone);
+  ALLOCATE (tree.Xo,   double *, tree.Nzone);
+  ALLOCATE (tree.Yo,   double *, tree.Nzone);
+  ALLOCATE (tree.dX,      int *, tree.Nzone);
+  ALLOCATE (tree.dY,      int *, tree.Nzone);
   ALLOCATE (tree.cell,    int *, tree.Nzone);
   ALLOCATE (tree.name,  char **, tree.Nzone);
@@ -152,4 +183,8 @@
     ALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
     ALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
+    ALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
+    ALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
+    ALLOCATE (tree.dX[zone],      int, tree.NBAND[zone]);
+    ALLOCATE (tree.dY[zone],      int, tree.NBAND[zone]);
     ALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
     ALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
@@ -157,4 +192,8 @@
       tree.ra[zone][band] = NAN;
       tree.dec[zone][band] = NAN;
+      tree.Xo[zone][band] = NAN;
+      tree.Yo[zone][band] = NAN;
+      tree.dX[zone][band] = -1;
+      tree.dY[zone][band] = -1;
       tree.cell[zone][band] = -1;
       ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH);
@@ -179,9 +218,17 @@
       REALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
       REALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
-      REALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
+      REALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
+      REALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
+      REALLOCATE (tree.dX[zone],      int, tree.NBAND[zone]);
+      REALLOCATE (tree.dY[zone],      int, tree.NBAND[zone]);
+      REALLOCATE (tree.cell[zone],    int, tree.NBAND[zone]);
       REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
       for (j = start; j < tree.NBAND[zone]; j++) {
 	tree.ra[zone][j] = NAN;
 	tree.dec[zone][j] = NAN;
+	tree.Xo[zone][band] = NAN;
+	tree.Yo[zone][band] = NAN;
+	tree.dX[zone][band] = -1;
+	tree.dY[zone][band] = -1;
 	tree.cell[zone][j] = -1;
 	ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH);
@@ -190,4 +237,9 @@
     tree.ra[zone][band] = ra;
     tree.dec[zone][band] = dec;
+    tree.Xo[zone][band] = x;
+    tree.Yo[zone][band] = y;
+    tree.dX[zone][band] = image[i].NX / NX_SUB;
+    tree.dY[zone][band] = image[i].NY / NY_SUB;
+    
     tree.cell[zone][band] = i;
 
@@ -325,7 +377,30 @@
     }
 
-    fprintf (stdout, "%10.6f %10.6f  %3d %3d  %s\n", ra, dec, zone, band, tree->name[zone][band]);
+    // I know the projection cell (band,zone), but I need to find the skycell within that projection cell.
+    // the proj cell is divided into Nx, Ny bits.
+    // (ra,dec) for (Ro,Do) -> (x,y).  given (Xo,Yo),(dX,dY) I can find ix,iy
+    // I currently track Ro,Do (tree->ra[zone][band], tree->dec[zone][band])
+    // I need to have the skycell center in pixels (Xo,Yo) and the scale dX,dY
+    
+    // convert R,D to X,Y with hard-wired projection and scale, orientation?
+    
+    double x = 0.0;
+    double y = 0.0;
+    BoundaryTreeProjection (&x, &y, ra, dec, tree, zone, band);
+
+    int xi = x / tree->dX[zone][band];
+    int yi = y / tree->dY[zone][band];
+    
+    char format[24], skycellname[128];
+    int Ndigit = (int)(log10(tree->NX_SUB*tree->NY_SUB)) + 1 ;
+    snprintf (format, 24, "%s.%%0%dd", tree->name[zone][band], Ndigit);
+
+    int N = xi + tree->NX_SUB * yi;
+    snprintf (skycellname, 128, format, N);
+
+    fprintf (stdout, "%10.6f %10.6f  %8.3f %8.3f  %3d %3d  %s\n", ra, dec, x, y, zone, band, skycellname);
   }
 
   exit (0);
 }
+
Index: trunk/Ohana/src/addstar/src/sky_tessalation.c
===================================================================
--- trunk/Ohana/src/addstar/src/sky_tessalation.c	(revision 34843)
+++ trunk/Ohana/src/addstar/src/sky_tessalation.c	(revision 34844)
@@ -275,5 +275,5 @@
       // convert the SkyRectangles to Images for output
       sky_subdivide_image (&image[j*NX_SUB*NY_SUB], &ring[j], NX_SUB, NY_SUB);
-      printf("%s %8.2f %8.2f\n", ring[j].name, ring[j].coords.crval1, ring[j].coords.crval2);
+      // printf("%s %8.2f %8.2f\n", ring[j].name, ring[j].coords.crval1, ring[j].coords.crval2);
     }
 
@@ -691,5 +691,5 @@
     // fprintf (stdout, "%f %f  %f  %f  %f %f  %f %f  %f %f %f\n", dec, decUpper, dRA, arg, root1, root2, theta1, theta2, r1, r2, k*r2);
   }
-  fprintf (stdout, "%f %f  %f x %f (%d)\n", dec, decUpper, dRA, theta, nRA);
+  // fprintf (stdout, "%f %f  %f x %f (%d)\n", dec, decUpper, dRA, theta, nRA);
 
   // I think we need to return the value of dec for the next ring, but I am not sure...
@@ -776,5 +776,5 @@
       ring[N].coords.crval2 = (j == 0) ? dec / d2r : -dec / d2r;
 
-      printf(" \t %d   %25.20f   %25.20f\n", i, ring[N].coords.crval2, ring[N].coords.crval1);
+      // printf(" \t %d   %25.20f   %25.20f\n", i, ring[N].coords.crval2, ring[N].coords.crval1);
 
       ring[N].coords.pc1_1 = +1.0 * X_PARITY;
Index: trunk/Ohana/src/addstar/test/boundarytree.sh
===================================================================
--- trunk/Ohana/src/addstar/test/boundarytree.sh	(revision 34844)
+++ trunk/Ohana/src/addstar/test/boundarytree.sh	(revision 34844)
@@ -0,0 +1,28 @@
+#!/usr/bin/env dvo
+
+# make a RINGS tess with the full skycells
+exec rm -r RINGS.V3
+exec skycells -mode RINGS -scale 0.25 -nx 10 -ny 10 -fix-ns -D CATDIR RINGS.V3 -overlap 60 60 -skyparity
+
+# make a RINGS tess with just projection cell
+exec rm -r RINGS.V3.proj
+exec skycells -mode RINGS -scale 2.50 -fix-ns -D CATDIR RINGS.V3.proj -overlap 60 60 -skyparity
+
+# build a boundary tree from RINGS.V3.proj
+exec findskycell -mktree RINGS.V3.tree.fits RINGS.V3.proj -nx 10 -ny 10 -scale 2.5
+
+delete -q r d
+concat 180.0 r; concat 0.1 d
+concat 180.1 r; concat 0.2 d
+concat 180.0 r; concat 0.2 d
+concat 180.5 r; concat 0.5 d
+write boundarytree.dat r d
+
+# try some test points
+exec findskycell -tree RINGS.V3.tree.fits boundarytree.dat
+
+catdir RINGS.V3.proj
+skyregion 0 360 -90 90
+for i 0 r[]
+  gimages -pix r[$i] d[$i]
+end
Index: trunk/Ohana/src/delstar/include/delstar.h
===================================================================
--- trunk/Ohana/src/delstar/include/delstar.h	(revision 34843)
+++ trunk/Ohana/src/delstar/include/delstar.h	(revision 34844)
@@ -27,4 +27,5 @@
 double ALPHA;
 int    VERBOSE;
+int    VERBOSE2;
 int    ORPHAN;
 int    MISSED;
@@ -38,4 +39,6 @@
 int          HOST_ID;
 char        *HOSTDIR;
+
+int SKIP_IMAGES;
 
 time_t    START;
Index: trunk/Ohana/src/delstar/src/args.c
===================================================================
--- trunk/Ohana/src/delstar/src/args.c	(revision 34843)
+++ trunk/Ohana/src/delstar/src/args.c	(revision 34844)
@@ -44,4 +44,11 @@
   }
 
+  VERBOSE2 = FALSE;
+  if ((N = get_argument (argc, argv, "-vv"))) {
+    VERBOSE = TRUE;
+    VERBOSE2 = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
   // region of interest
   UserPatch.Rmin = 0;
@@ -58,4 +65,10 @@
     remove_argument (N, &argc, argv);
     UserPatch.Dmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  SKIP_IMAGES = FALSE;
+  if ((N = get_argument (argc, argv, "-skip-images"))) {
+    SKIP_IMAGES = TRUE;
     remove_argument (N, &argc, argv);
   }
@@ -205,4 +218,11 @@
   }
 
+  VERBOSE2 = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    VERBOSE = TRUE;
+    VERBOSE2 = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
   /* specify portion of the sky */
   UserPatch.Rmin = 0;
Index: trunk/Ohana/src/delstar/src/delete_photcodes.c
===================================================================
--- trunk/Ohana/src/delstar/src/delete_photcodes.c	(revision 34843)
+++ trunk/Ohana/src/delstar/src/delete_photcodes.c	(revision 34844)
@@ -46,10 +46,10 @@
 
     // an error exit status here is a significant error
-    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "a")) {
+    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "a")) {
       fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename);
       exit (2);
     }
     if (!catalog.Naves_disk) {
-      if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
+      if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
       dvo_catalog_unlock (&catalog);
       dvo_catalog_free (&catalog);
@@ -58,5 +58,5 @@
 
     delete_photcodes_catalog (&catalog, photcodes, Nphotcodes);
-    dvo_catalog_save_complete (&catalog, VERBOSE);
+    dvo_catalog_save_complete (&catalog, VERBOSE2);
     dvo_catalog_unlock (&catalog);
     dvo_catalog_free (&catalog);
@@ -97,4 +97,5 @@
     char tmpline[1024];
     if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                    strcpy (command, tmpline); }
+    if (VERBOSE2)      { snprintf (tmpline, 1024, "%s -vv",             command);                    strcpy (command, tmpline); }
 
     fprintf (stderr, "command: %s\n", command);
@@ -287,4 +288,5 @@
   FREE(measure);
   FREE(average);
+  FREE(secfilt);
   catalog[0].measure = measureOut;
   catalog[0].average = averageOut;
@@ -299,13 +301,13 @@
   for (i = 0; i < NaveOut; i++) {
     NmeasureTotal += catalog[0].average[i].Nmeasure;
-    if (VERBOSE && !(NmeasureTotal <= catalog[0].Nmeasure)) {
+    if (VERBOSE2 && !(NmeasureTotal <= catalog[0].Nmeasure)) {
       fprintf (stderr, "too few measurements: %d %d %d\n", (int) i, NmeasureTotal, (int) catalog[0].Nmeasure);
     }
     measureOffsetOK &= (catalog[0].average[i].measureOffset < catalog[0].Nmeasure);
-    if (VERBOSE && !(catalog[0].average[i].measureOffset < catalog[0].Nmeasure)) {
+    if (VERBOSE2 && !(catalog[0].average[i].measureOffset < catalog[0].Nmeasure)) {
       fprintf (stderr, "offset too large: %d %d %d\n", (int) i, catalog[0].average[i].Nmeasure, (int) catalog[0].Nmeasure);
     }
     measureOffsetOK &= (catalog[0].average[i].measureOffset + catalog[0].average[i].Nmeasure <= catalog[0].Nmeasure);
-    if (VERBOSE && !(catalog[0].average[i].measureOffset + catalog[0].average[i].Nmeasure <= catalog[0].Nmeasure)) {
+    if (VERBOSE2 && !(catalog[0].average[i].measureOffset + catalog[0].average[i].Nmeasure <= catalog[0].Nmeasure)) {
       fprintf (stderr, "orrset + Nmeasure too large: %d + %d > %d %d\n", (int) i, catalog[0].average[i].measureOffset, catalog[0].average[i].Nmeasure, (int) catalog[0].Nmeasure);
     }
@@ -313,5 +315,5 @@
     for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
       averefOK &= (catalog[0].measure[m+j].averef == i);
-      if (VERBOSE && !(catalog[0].measure[m+j].averef == i)) {
+      if (VERBOSE2 && !(catalog[0].measure[m+j].averef == i)) {
 	fprintf (stderr, "averef broken: %d vs %d (measure %d)\n", (int) i, catalog[0].measure[m+j].averef, (int) (m+j));
       }
@@ -340,4 +342,5 @@
   FREE (measureSeqRaw);
   FREE (measureSeqOut);
+  FREE (measureRefOut);
   FREE (measureAveRaw);
   FREE (averageNmeas);
Index: trunk/Ohana/src/delstar/src/delstar.c
===================================================================
--- trunk/Ohana/src/delstar/src/delstar.c	(revision 34843)
+++ trunk/Ohana/src/delstar/src/delstar.c	(revision 34844)
@@ -10,11 +10,12 @@
   args (argc, argv);
 
-  set_db (&db);
-  status = dvo_image_lock (&db, ImageCat, 60.0, LCK_XCLD); /* XCLD */
-  if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
-  if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
-
-  status = dvo_image_load (&db, VERBOSE, FALSE);
-  if (!status) Shutdown ("can't read image catalog %s", db.filename);
+  if (!SKIP_IMAGES) {
+    set_db (&db);
+    status = dvo_image_lock (&db, ImageCat, 60.0, LCK_XCLD); /* XCLD */
+    if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
+    if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
+    status = dvo_image_load (&db, VERBOSE, FALSE);
+    if (!status) Shutdown ("can't read image catalog %s", db.filename);
+  }
 
   switch (MODE) {
@@ -29,5 +30,7 @@
       break;
     case MODE_PHOTCODES: {
-      if (!delete_image_photcodes (&db)) exit (1);
+      if (!SKIP_IMAGES) {
+	if (!delete_image_photcodes (&db)) exit (1);
+      }
 
       if (!delete_photcodes (0, NULL)) {
Index: trunk/Ohana/src/delstar/src/delstar_client.c
===================================================================
--- trunk/Ohana/src/delstar/src/delstar_client.c	(revision 34843)
+++ trunk/Ohana/src/delstar/src/delstar_client.c	(revision 34844)
@@ -32,6 +32,9 @@
       break;
     case MODE_PHOTCODES:
-      delete_photcodes (HOST_ID, HOSTDIR);
-      break;
+      if (!delete_photcodes (HOST_ID, HOSTDIR)) {
+	fprintf (stderr, "failure deleting measurements from %s\n", HOSTDIR);
+	exit (1);
+      }
+      exit (0);
     case MODE_ORPHAN:
       break;
Index: trunk/Ohana/src/libdvo/include/dvo.h
===================================================================
--- trunk/Ohana/src/libdvo/include/dvo.h	(revision 34843)
+++ trunk/Ohana/src/libdvo/include/dvo.h	(revision 34844)
@@ -322,8 +322,17 @@
   int *NBAND;
 
-  double   **ra;
-  double  **dec;
-  int    **cell;
-  char  ***name;
+  double   **ra; // RA of projection cell center
+  double  **dec; // DEC of projection cell center
+  int    **cell; // zone,band -> proj cell sequence
+  char  ***name; // projection cell name
+  
+  int NX_SUB;
+  int NY_SUB;
+  double dPix;
+
+  double **Xo;
+  double **Yo;
+  int **dX;
+  int **dY;
 } BoundaryTree;
 
@@ -496,4 +505,5 @@
 float PhotKronInst (Measure *measure);
 float PhotKronAve (PhotCode *code, Average *average, SecFilt *secfilt);
+float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt);
 
 float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt);
@@ -701,4 +711,5 @@
 int BoundaryTreeSave(char *filename, BoundaryTree *tree);
 BoundaryTree *BoundaryTreeLoad(char *filename);
+int BoundaryTreeProjection (double *x, double *y, double r, double d, BoundaryTree *tree, int zone, int band);
 
 void dvo_average_init (Average *average);
Index: trunk/Ohana/src/libdvo/src/BoundaryTree.c
===================================================================
--- trunk/Ohana/src/libdvo/src/BoundaryTree.c	(revision 34843)
+++ trunk/Ohana/src/libdvo/src/BoundaryTree.c	(revision 34844)
@@ -47,4 +47,8 @@
   gfits_scan (&header, "DEC_ORI", "%lf", 1, &tree->DEC_origin);
   gfits_scan (&header, "DEC_OFF", "%lf", 1, &tree->DEC_offset);
+
+  gfits_scan (&header, "NX_SUB", "%d", 1, &tree->NX_SUB);
+  gfits_scan (&header, "NY_SUB", "%d", 1, &tree->NY_SUB);
+  gfits_scan (&header, "PIXSCALE", "%lf", 1, &tree->dPix);
 
   ftable.header = &theader;
@@ -75,10 +79,18 @@
     ALLOCATE (tree->ra,   double *, tree->Nzone);
     ALLOCATE (tree->dec,  double *, tree->Nzone);
-    ALLOCATE (tree->cell, int *, tree->Nzone);
-    ALLOCATE (tree->name, char **, tree->Nzone);
+    ALLOCATE (tree->Xo,   double *, tree->Nzone);
+    ALLOCATE (tree->Yo,   double *, tree->Nzone);
+    ALLOCATE (tree->dX,      int *, tree->Nzone);
+    ALLOCATE (tree->dY,      int *, tree->Nzone);
+    ALLOCATE (tree->cell,    int *, tree->Nzone);
+    ALLOCATE (tree->name,  char **, tree->Nzone);
     for (i = 0; i < tree->Nzone; i++) {
       ALLOCATE (tree->ra[i],   double, tree->Nband[i]);
       ALLOCATE (tree->dec[i],  double, tree->Nband[i]);
-      ALLOCATE (tree->cell[i], int,    tree->Nband[i]);
+      ALLOCATE (tree->Xo[i],   double, tree->Nband[i]);
+      ALLOCATE (tree->Yo[i],   double, tree->Nband[i]);
+      ALLOCATE (tree->dX[i],      int, tree->Nband[i]);
+      ALLOCATE (tree->dY[i],      int, tree->Nband[i]);
+      ALLOCATE (tree->cell[i],    int, tree->Nband[i]);
       ALLOCATE (tree->name[i], char *, tree->Nband[i]);
       for (j = 0; j < tree->Nband[i]; j++) {
@@ -102,4 +114,8 @@
     GET_COLUMN_NEW(band,  "BAND",        int);
     GET_COLUMN_NEW(index, "INDEX",       int);
+    GET_COLUMN_NEW(Xo,    "X_CENT",      double);
+    GET_COLUMN_NEW(Yo,    "Y_CENT",      double);
+    GET_COLUMN_NEW(dX,    "X_GRID",      int);
+    GET_COLUMN_NEW(dY,    "Y_GRID",      int);
     GET_COLUMN_NEW(name,  "NAME",        char); // XXX how is this done?
     gfits_free_header (&theader);
@@ -113,4 +129,8 @@
       tree->ra[nz][nb] = R[i];
       tree->dec[nz][nb] = D[i];
+      tree->Xo[nz][nb] = Xo[i];
+      tree->Yo[nz][nb] = Yo[i];
+      tree->dX[nz][nb] = dX[i];
+      tree->dY[nz][nb] = dY[i];
       tree->cell[nz][nb] = i; // XXX ?
       memcpy(tree->name[nz][nb], &name[i*BOUNDARY_TREE_NAME_LENGTH], BOUNDARY_TREE_NAME_LENGTH);
@@ -121,4 +141,8 @@
     free (zone  );
     free (band  );
+    free (Xo    );
+    free (Yo    );
+    free (dX    );
+    free (dY    );
     free (index );
     free (name  );
@@ -165,4 +189,8 @@
   gfits_modify (&header, "DEC_ORI", "%lf", 1, tree->DEC_origin);
   gfits_modify (&header, "DEC_OFF", "%lf", 1, tree->DEC_offset);
+
+  gfits_modify (&header, "NX_SUB", "%d", 1, tree->NX_SUB);
+  gfits_modify (&header, "NY_SUB", "%d", 1, tree->NY_SUB);
+  gfits_modify (&header, "PIXSCALE", "%lf", 1, tree->dPix);
 
   gfits_fwrite_header  (f, &header);
@@ -223,4 +251,8 @@
     gfits_define_bintable_column (&theader, "J", "BAND", "band sequence number", "none", 1.0, 0.0);
     gfits_define_bintable_column (&theader, "J", "INDEX","cell index", "none", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "D", "X_CENT", "projection cell center pixel", "none", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "D", "Y_CENT", "projection cell center pixel", "none", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "J", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "J", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0);
     gfits_define_bintable_column (&theader, fmt, "NAME", "cell name", "none", 1.0, 0.0);
 
@@ -240,4 +272,8 @@
     int    *band          ; ALLOCATE (band,  int,    Ncell);
     int    *index         ; ALLOCATE (index, int,    Ncell);
+    double *Xo            ; ALLOCATE (Xo,    double, Ncell);
+    double *Yo            ; ALLOCATE (Yo,    double, Ncell);
+    int    *dX            ; ALLOCATE (dX,    int,    Ncell);
+    int    *dY            ; ALLOCATE (dY,    int,    Ncell);
     char   *name          ; ALLOCATE (name,  char,   Ncell*BOUNDARY_TREE_NAME_LENGTH);
 
@@ -251,4 +287,8 @@
 	R[i]     = tree->ra[nz][nb];
 	D[i]     = tree->dec[nz][nb];
+	Xo[i]    = tree->Xo[nz][nb];
+	Yo[i]    = tree->Yo[nz][nb];
+	dX[i]    = tree->dX[nz][nb];
+	dY[i]    = tree->dY[nz][nb];
 	zone[i]  = nz;
 	band[i]  = nb;
@@ -265,4 +305,8 @@
     gfits_set_bintable_column (&theader, &ftable, "BAND",  band,  Ncell);
     gfits_set_bintable_column (&theader, &ftable, "INDEX", index, Ncell);
+    gfits_set_bintable_column (&theader, &ftable, "X_CENT", Xo,   Ncell);
+    gfits_set_bintable_column (&theader, &ftable, "Y_CENT", Yo,   Ncell);
+    gfits_set_bintable_column (&theader, &ftable, "X_GRID", dX,   Ncell);
+    gfits_set_bintable_column (&theader, &ftable, "Y_GRID", dY,   Ncell);
     gfits_set_bintable_column (&theader, &ftable, "NAME",  name,  Ncell);
 
@@ -272,4 +316,8 @@
     free (band  );
     free (index );
+    free (Xo    );
+    free (Yo    );
+    free (dX    );
+    free (dY    );
     free (name  );
 
@@ -318,2 +366,58 @@
 }
 
+
+// projection = TAN
+// need Ro, Do, Xo, Yo, dPix
+
+int BoundaryTreeProjection (double *x, double *y, double r, double d, BoundaryTree *tree, int zone, int band) {
+
+    double Xo = tree->Xo[zone][band];
+    double Yo = tree->Yo[zone][band];
+    double Ro = tree->ra[zone][band];
+    double Do = tree->dec[zone][band];
+    double dPix = tree->dPix;
+
+    // this block only depends on Ro, Do
+
+    double sdp  = sin(RAD_DEG*Do);
+    double cdp  = cos(RAD_DEG*Do);
+    double salp = sin(RAD_DEG*(r - Ro));
+    double calp = cos(RAD_DEG*(r - Ro));
+    double sdel = sin(RAD_DEG*d);
+    double cdel = cos(RAD_DEG*d);
+    
+    double stht = sdel*sdp + cdel*cdp*calp;    /* sin(theta) */
+    double sphi = cdel*salp;                   /* = cos(theta)*sin(phi) */
+    double cphi = cdel*sdp*calp - sdel*cdp;    /* = cos(theta)*cos(phi) */
+
+    // defines the TAN projection (one of zenithal projections available, libdvo/src/coordops.c
+    // R = cot (theta) = cos(theta) / sin(theta)
+    double L, M;
+    if (stht == 0) {
+	double Rc = hypot(sphi, cphi);
+	L = 180.0 * sphi / Rc;
+	M = 180.0 * cphi / Rc;
+    } else {
+	L = +DEG_RAD * sphi / stht;
+	M = -DEG_RAD * cphi / stht;
+    }
+
+    // scale, rotation, parity:
+    // rotation == 0.0 (pc1_1 == pc2_2 == 1.0, pc1_2 = pc2_1 = 0.0)
+
+    // if there were rotation or parity:
+    // Ro = (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
+    // Xo = (coords[0].pc2_2*L - coords[0].pc1_2*M) / Ro;
+    // Yo = (coords[0].pc1_1*M - coords[0].pc2_1*L) / Ro;
+
+    double X = L;
+    double Y = M;
+
+    // scale is dPix
+
+    *x = Xo - X / dPix;
+    *y = Yo + Y / dPix;
+    
+    return TRUE;
+}
+
Index: trunk/Ohana/src/libdvo/src/dbExtractAverages.c
===================================================================
--- trunk/Ohana/src/libdvo/src/dbExtractAverages.c	(revision 34843)
+++ trunk/Ohana/src/libdvo/src/dbExtractAverages.c	(revision 34844)
@@ -233,4 +233,8 @@
 	  break;
 
+	case MAG_KRON_ERR:
+	  value.Flt = PhotKronAveErr (field->photcode, average, secfilt);
+	  break;
+
 	case MAG_20:
 	  value.Flt = PhotM20 (field->photcode, average, secfilt);
Index: trunk/Ohana/src/libdvo/src/dvo_catalog_split.c
===================================================================
--- trunk/Ohana/src/libdvo/src/dvo_catalog_split.c	(revision 34843)
+++ trunk/Ohana/src/libdvo/src/dvo_catalog_split.c	(revision 34844)
@@ -863,5 +863,5 @@
     first  = 0;                    // first row in memory to write
     start  = catalog[0].Nsecf_off; // first disk row to write
-    Nitems = catalog[0].Nsecf_mem;
+    Nitems = catalog[0].Naverage*Nsecfilt;
     Nrows  = Nitems - first;
 
Index: trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c
===================================================================
--- trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c	(revision 34843)
+++ trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c	(revision 34844)
@@ -572,4 +572,16 @@
 }
 
+float PhotKronAveErr (PhotCode *code, Average *average, SecFilt *secfilt) {
+
+  int Ns;
+  float dMkron;
+
+  if (code == NULL) return NAN;
+
+  Ns = photcodes[0].hashNsec[code[0].code];
+  dMkron = (Ns == -1) ? NAN : secfilt[Ns].dMkron;
+  return (dMkron);
+}
+
 float PhotAveFluxPSF (PhotCode *code, Average *average, SecFilt *secfilt) {
 
Index: trunk/Ohana/src/opihi/dvo/gimages.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/gimages.c	(revision 34843)
+++ trunk/Ohana/src/opihi/dvo/gimages.c	(revision 34844)
@@ -212,4 +212,6 @@
 	status = XY_to_RD (&Rout, &Dout, X, Y, &image[i].coords);
 	fprintf (stderr, "r,d = %f,%f\n", Rout, Dout);
+	status = RD_to_XY (&X, &Y, Ra, Dec, &image[i].coords);
+	fprintf (stderr, "x,y = %f,%f\n", X, Y);
     }
 # endif
Index: trunk/Ohana/src/opihi/dvo/gstar.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/gstar.c	(revision 34843)
+++ trunk/Ohana/src/opihi/dvo/gstar.c	(revision 34844)
@@ -1,3 +1,20 @@
 # include "dvoshell.h"
+
+typedef enum {
+  GSTAR_AVE_MAG,
+  GSTAR_AVE_MAG_ERR,
+  GSTAR_AVE_MAG_CHISQ,
+  GSTAR_FILTER_NAMES,
+  GSTAR_SECF_FLAGS,
+  GSTAR_AVE_AP_MAG,
+  GSTAR_AVE_MAG_20,
+  GSTAR_AVE_MAG_80,
+  GSTAR_AVE_KRON_MAG,
+  GSTAR_AVE_KRON_MAG_ERR,
+  GSTAR_STACK_FLUX_PSF,
+  GSTAR_STACK_FLUX_PSF_ERR,
+  GSTAR_STACK_FLUX_KRON,
+  GSTAR_STACK_FLUX_KRON_ERR,
+} GSTAR_SECF_CODES;
 
 void initPhotcodeSequence (int Nsecfilt);
@@ -13,5 +30,5 @@
   off_t i, Nstars, *N1;
   off_t j, k, m, N, Nlo, Nhi;
-  int Nsecfilt, NPTS, QUIET, FULL_OUTPUT, INST;
+  int Nsecfilt, NPTS, QUIET, FULL_OUTPUT, STACK_OUTPUT, INST;
   int found, GetMeasures;
   int SaveVectors;
@@ -55,4 +72,10 @@
   if ((N = get_argument (argc, argv, "-full"))) {
     FULL_OUTPUT = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  STACK_OUTPUT = FALSE;
+  if ((N = get_argument (argc, argv, "-stack"))) {
+    STACK_OUTPUT = TRUE;
     remove_argument (N, &argc, argv);
   }
@@ -216,34 +239,60 @@
       
 	/* filter names */
-	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 3);
+	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_FILTER_NAMES);
 	gprint (GP_LOG, "\n");
 
 	/* average mags */
-	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 0);
+	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG);
 	gprint (GP_LOG, "\n");
 
 	/* average mag errors */
-	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 1);
+	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_ERR);
 	gprint (GP_LOG, "\n");
 
 	/* average mag chisq */
-	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 2);
+	for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_CHISQ);
 	gprint (GP_LOG, "\n");
 
 	if (FULL_OUTPUT) {
 	  /* Map */
-	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 5);
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_AP_MAG);
+	  gprint (GP_LOG, "\n");
+
+	  /* Mkron */
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_KRON_MAG);
+	  gprint (GP_LOG, "\n");
+
+	  /* dMkron */
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_KRON_MAG_ERR);
 	  gprint (GP_LOG, "\n");
 
 	  /* M_20 */
-	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 6);
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_20);
 	  gprint (GP_LOG, "\n");
 
 	  /* M_80 */
-	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 7);
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_80);;
 	  gprint (GP_LOG, "\n");
 
 	  /* secfilt flags */
-	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, 4);
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_SECF_FLAGS);
+	  gprint (GP_LOG, "\n");
+
+	}
+	if (STACK_OUTPUT) {
+	  /* FluxPSF */
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_PSF);
+	  gprint (GP_LOG, "\n");
+
+	  /* dFluxPSF */
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_PSF_ERR);
+	  gprint (GP_LOG, "\n");
+
+	  /* FluxKron */
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON);
+	  gprint (GP_LOG, "\n");
+
+	  /* dFluxKron */
+	  for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON_ERR);
 	  gprint (GP_LOG, "\n");
 
@@ -279,4 +328,6 @@
 		gprint (GP_LOG, "%6.3f ", catalog.measure[m].Mcal);
 		gprint (GP_LOG, "%6.3f ", catalog.measure[m].Map);
+		gprint (GP_LOG, "%6.3f ", catalog.measure[m].Mkron);
+		gprint (GP_LOG, "%6.3f ", catalog.measure[m].dMkron);
 		gprint (GP_LOG, "%5.1f ", pow(10.0, 0.4*catalog.measure[m].dt));
 		gprint (GP_LOG, "%5.3f ", catalog.measure[m].airmass);
@@ -296,4 +347,10 @@
 		gprint (GP_LOG, "%.1f ", catalog.measure[m].extNsigma);
 		gprint (GP_LOG, "%.1f ", FromShortDegrees(catalog.measure[m].theta));
+	    }
+	    if (STACK_OUTPUT) {
+		gprint (GP_LOG, "%10.3e ", catalog.measure[m].FluxPSF);
+		gprint (GP_LOG, "%10.3e ", catalog.measure[m].dFluxPSF);
+		gprint (GP_LOG, "%10.3e ", catalog.measure[m].FluxKron);
+		gprint (GP_LOG, "%10.3e ", catalog.measure[m].dFluxKron);
 	    }
 	    gprint (GP_LOG, "\n");
@@ -351,4 +408,11 @@
 }
 
+void print_double_exp (double value) {
+  if (isnan(value)) 
+    gprint (GP_LOG, "  NaN    ");
+  else 
+    gprint (GP_LOG, "%10.3e ", value);
+}
+
 void print_short (double value, short int ival) {
   if (ival == NAN_S_SHORT) 
@@ -386,5 +450,4 @@
 }
 
-
 void printPhotcodeSequence (Average *average, SecFilt *secfilt, int entry, int type) {
 
@@ -395,5 +458,5 @@
 
   switch (type) {
-    case 0: /* average mags */
+    case GSTAR_AVE_MAG: /* average mags */
       if (seq == -1) {
 	print_double (NAN);
@@ -403,5 +466,5 @@
       break;
 
-    case 1: /* average mags errors */
+    case GSTAR_AVE_MAG_ERR: /* average mags errors */
       if (seq == -1) {
 	print_double (NAN);
@@ -411,5 +474,5 @@
       break;
 
-    case 2: /* average mag chisq */
+    case GSTAR_AVE_MAG_CHISQ: /* average mag chisq */
       if (seq == -1) {
 	print_short (NAN_S_SHORT, NAN_S_SHORT);
@@ -419,10 +482,10 @@
       break;
 
-    case 3: /* filter names */
+    case GSTAR_FILTER_NAMES: /* filter names */
       code = GetPhotcodebyNsec (seq);
       gprint (GP_LOG, "  %-6s ", code[0].name);
       break;
 
-    case 4: /* secfilt flags */
+    case GSTAR_SECF_FLAGS: /* secfilt flags */
       if (seq == -1) {
 	gprint (GP_LOG, "%08x ", 0);
@@ -432,5 +495,5 @@
       break;
 
-    case 5: /* average ap mags */
+    case GSTAR_AVE_AP_MAG: /* average ap mags */
       if (seq == -1) {
 	print_double (NAN);
@@ -440,5 +503,5 @@
       break;
 
-    case 6: /* average ap mags */
+    case GSTAR_AVE_MAG_20: /* average ap mags */
       if (seq == -1) {
 	print_double (NAN);
@@ -448,12 +511,60 @@
       break;
 
-    case 7: /* average ap mags */
-      if (seq == -1) {
-	print_double (NAN);
-      } else {
-	print_short (0.001*secfilt[seq].M_20, secfilt[seq].M_80);
-      }
-      break;
-
-  }
-}
+    case GSTAR_AVE_MAG_80: /* average ap mags */
+      if (seq == -1) {
+	print_double (NAN);
+      } else {
+	print_short (0.001*secfilt[seq].M_80, secfilt[seq].M_80);
+      }
+      break;
+
+    case GSTAR_AVE_KRON_MAG: /* average ap mags */
+      if (seq == -1) {
+	print_double (NAN);
+      } else {
+	print_double (secfilt[seq].Mkron);
+      }
+      break;
+
+    case GSTAR_AVE_KRON_MAG_ERR: /* average ap mags */
+      if (seq == -1) {
+	print_double (NAN);
+      } else {
+	print_double (secfilt[seq].dMkron);
+      }
+      break;
+
+    case GSTAR_STACK_FLUX_PSF: /* average ap mags */
+      if (seq == -1) {
+	print_double (NAN);
+      } else {
+	print_double_exp (secfilt[seq].FluxPSF);
+      }
+      break;
+
+    case GSTAR_STACK_FLUX_PSF_ERR: /* average ap mags */
+      if (seq == -1) {
+	print_double (NAN);
+      } else {
+	print_double_exp (secfilt[seq].dFluxPSF);
+      }
+      break;
+
+    case GSTAR_STACK_FLUX_KRON: /* average ap mags */
+      if (seq == -1) {
+	print_double (NAN);
+      } else {
+	print_double_exp (secfilt[seq].FluxKron);
+      }
+      break;
+
+    case GSTAR_STACK_FLUX_KRON_ERR: /* average ap mags */
+      if (seq == -1) {
+	print_double (NAN);
+      } else {
+	print_double_exp (secfilt[seq].dFluxKron);
+      }
+      break;
+
+  }
+}
Index: trunk/Ohana/src/opihi/pcontrol/HostOps.c
===================================================================
--- trunk/Ohana/src/opihi/pcontrol/HostOps.c	(revision 34843)
+++ trunk/Ohana/src/opihi/pcontrol/HostOps.c	(revision 34844)
@@ -232,5 +232,4 @@
   FreeIOBuffer (&host[0].comms_buffer);
   FREE (host[0].hostname);
-  FREE (host[0].job);
   FREE (host);
 }
Index: trunk/Ohana/src/relphot/Makefile
===================================================================
--- trunk/Ohana/src/relphot/Makefile	(revision 34843)
+++ trunk/Ohana/src/relphot/Makefile	(revision 34844)
@@ -49,4 +49,5 @@
 $(SRC)/setExclusions.$(ARCH).o 	 \
 $(SRC)/setMrelFinal.$(ARCH).o 	 \
+$(SRC)/setMrelCatalog.$(ARCH).o 	 \
 $(SRC)/BoundaryTreeOps.$(ARCH).o 	 \
 $(SRC)/write_coords.$(ARCH).o
@@ -78,4 +79,5 @@
 $(SRC)/setExclusions.$(ARCH).o 	 \
 $(SRC)/setMrelFinal.$(ARCH).o    \
+$(SRC)/setMrelCatalog.$(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 34843)
+++ trunk/Ohana/src/relphot/include/relphot.h	(revision 34844)
@@ -63,4 +63,33 @@
 
 typedef struct {
+  int Nfew;
+  int Ncode;
+  int Nsys;
+  int Nbad;
+  int Ncal;
+  int Nmos;
+  int Ngrid;
+
+  // NOTE: the following arrays are (possibly) pre-allocated and carried down to each
+  // thread.  The first 3 (list, dlist, wlist) are used in all relphot analyses; the
+  // others are only used on the final output steps.
+
+  double *list;		      // list of measure.mag values for a given star
+  double *dlist;	      // mag errors for a star
+  double *wlist;	      // weights to use for mean mags
+  double *aplist;	      // ap mags for a star
+  double *kronlist;	      // kron mags for a star
+  double *dkronlist;	      // kron mag errors 
+
+  double *psfqf_list;	      // psfqf for all filters
+  double *psfqfperf_list;     // psfqfperf for all filters
+  double *stargal_list;	      // stargal for all filters
+
+  StatType stats;
+  StatType apstats;
+  StatType kronstats;
+} SetMrelInfo;
+
+typedef struct {
   AverageTiny *average;	      // array of (minimal) average data
   MeasureTiny *measure;	      // array of (minimal) measure data 
@@ -113,4 +142,6 @@
 
 char        *BOUNDARY_TREE;
+
+int SET_MREL_VERSION;
 
 double MAG_LIM;
@@ -239,4 +270,6 @@
 Coords       *getCoords           PROTO((off_t meas, int cat));
 off_t         getImageEntry       PROTO((off_t meas, int cat));
+
+float         getMcal_alt         PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, float Xccd, float Yccd));
 float         getMcal             PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
 float         getMflat            PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
@@ -348,2 +381,10 @@
 int load_tree (char *treefile);
 int BoundaryTreePrimaryCell (char *primaryCellName, double ra, double dec);
+
+int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure);
+int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
+int setMrelAverageExposure (off_t meas, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found);
+int setMrelAverageStack (off_t meas, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell);
+int setGlobalObjStats (Average *average, Measure *measure);
+
+void SetMrelInfoInit (SetMrelInfo *results, int allocLists);
Index: trunk/Ohana/src/relphot/src/BoundaryTreeOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/BoundaryTreeOps.c	(revision 34843)
+++ trunk/Ohana/src/relphot/src/BoundaryTreeOps.c	(revision 34844)
@@ -21,5 +21,29 @@
   }
   
+# define USE_PROJECTION_CELL 0
+# if (USE_PROJECTION_CELL)
   snprintf (primaryCellName, DVO_MAX_PATH, "RINGS.V3.%s", tree->name[zone][band]);
+# else
+
+  // I have ra, dec, and the primary projection cell.  In order to choose the primary skycell,
+  // I just need to project to ra,dec to X,Y based on the center of the cell and then get the subdivision right.
+    
+  double x = 0.0;
+  double y = 0.0;
+  BoundaryTreeProjection (&x, &y, ra, dec, tree, zone, band);
+  
+  int xi = x / tree->dX[zone][band];
+  int yi = y / tree->dY[zone][band];
+  int N = xi + tree->NX_SUB * yi;
+  
+  // XXX short-circuit this for now (we should use this code if we make tree have more variable NX,NY values
+  // char format[24], skycellname[128];
+  // int Ndigit = (int)(log10(tree->NX_SUB*tree->NY_SUB)) + 1 ;
+  // snprintf (format, 24, "%s.%%0%dd", tree->name[zone][band], Ndigit);
+  // snprintf (skycellname, 128, format, N);
+
+  snprintf (primaryCellName, DVO_MAX_PATH, "RINGS.V3.%s.%03d", tree->name[zone][band], N);
+# endif
+
   return TRUE;
 }
Index: trunk/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/ImageOps.c	(revision 34843)
+++ trunk/Ohana/src/relphot/src/ImageOps.c	(revision 34844)
@@ -352,4 +352,27 @@
 
 // returns image.Mcal - ff(x,y)
+float getMcal_alt (off_t meas, int cat, FlatCorrectionTable *flatcorr, float Xccd, float Yccd) {
+
+  off_t i;
+  float value, offset;
+
+  i = MeasureToImage[cat][meas];
+  if (i == -1) return (NAN);
+
+  if (image[i].flags & IMAGE_BAD) return (NAN);  
+  value = image[i].Mcal;
+  offset = 0.0;
+
+  // to do this, I need to pass in the catalog and flatcorr pointers
+  int flat_id = image[i].photom_map_id;
+  if (flat_id) {
+    offset = FlatCorrectionOffset (flatcorr, flat_id, Xccd, Yccd);
+  }
+  value -= offset;
+
+  return (value);
+}
+
+// returns image.Mcal - ff(x,y)
 short getUbercalDist (off_t meas, int cat) {
 
@@ -400,7 +423,9 @@
   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
+  // 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 or RINGS.V3.skycell.1495.027
+  // (if we use projection or skycell as the primary)
 
   if (!strncmp(image[i].name, name, strlen(name))) return TRUE;
Index: trunk/Ohana/src/relphot/src/StarOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/StarOps.c	(revision 34843)
+++ trunk/Ohana/src/relphot/src/StarOps.c	(revision 34844)
@@ -2,29 +2,4 @@
 
 static int Nmax;
-
-typedef struct {
-  int Nfew;
-  int Ncode;
-  int Nsys;
-  int Nbad;
-  int Ncal;
-  int Nmos;
-  int Ngrid;
-
-  // NOTE: the following arrays are (possibly) pre-allocated and carried down to each
-  // thread.  The first 3 (list, dlist, wlist) are used in all relphot analyses; the
-  // others are only used on the final output steps.
-
-  double *list;		      // list of measure.mag values for a given star
-  double *dlist;	      // mag errors for a star
-  double *wlist;	      // weights to use for mean mags
-  double *aplist;	      // ap mags for a star
-  double *kronlist;	      // kron mags for a star
-  double *dkronlist;	      // kron mag errors 
-
-  double *psfqf_list;	      // psfqf for all filters
-  double *psfqfperf_list;     // psfqfperf for all filters
-  double *stargal_list;	      // stargal for all filters
-} SetMrelInfo;
 
 enum {THREAD_RUN, THREAD_DONE};
@@ -154,5 +129,15 @@
   for (i = 0; i < Ncatalog; i++) {
     // pass == -1 for anything other than the final pass
-    setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
+    switch (SET_MREL_VERSION) {
+      case 0:
+	setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
+	break;
+      case 1:
+	setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt);
+	break;
+      default:
+	fprintf (stderr, "invalid setMrel version %d (use 0 or 1)\n", SET_MREL_VERSION);
+	exit (2);
+    }
     SetMrelInfoAccum (&summary, &results);
   }
@@ -182,5 +167,15 @@
 
   for (i = 0; i < Ncatalog; i++) {
-    setMrel_catalog  (catalog, i, pass, flatcorr, &results, Nsecfilt); // XXX add arguments as needed for options
+    switch (SET_MREL_VERSION) {
+      case 0:
+	setMrel_catalog (catalog, i, pass, flatcorr, &results, Nsecfilt);
+	break;
+      case 1:
+	setMrel_catalog_alt (catalog, i, pass, flatcorr, &results, Nsecfilt);
+	break;
+      default:
+	fprintf (stderr, "invalid setMrel version %d (use 0 or 1)\n", SET_MREL_VERSION);
+	exit (2);
+    }
     SetMrelInfoAccum (&summary, &results);
   }
@@ -286,5 +281,15 @@
 
     // pass == -1 for anything other than the final pass
-    setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
+    switch (SET_MREL_VERSION) {
+      case 0:
+	setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
+	break;
+      case 1:
+	setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt);
+	break;
+      default:
+	fprintf (stderr, "invalid setMrel version %d (use 0 or 1)\n", SET_MREL_VERSION);
+	exit (2);
+    }
     SetMrelInfoAccum (&threadinfo->summary, &results);
   }
@@ -332,7 +337,4 @@
   liststats_setmode (&kronstats, STATMODE);
 
-  // StatType psfstats;
-  // liststats_setmode (&psfstats, "MEAN");
-
   double *list      = results->list;
   double *dlist     = results->dlist;
@@ -363,4 +365,5 @@
 
     if (isSetMrelFinal) {
+      // set the name of the primary skycell (this is used in a strcmp to match the skycells in stack detections)
       BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
     }
Index: trunk/Ohana/src/relphot/src/args.c
===================================================================
--- trunk/Ohana/src/relphot/src/args.c	(revision 34843)
+++ trunk/Ohana/src/relphot/src/args.c	(revision 34844)
@@ -64,4 +64,11 @@
   }
 
+  SET_MREL_VERSION = 1;
+  if ((N = get_argument (argc, argv, "-set-mrel-version"))) {
+    remove_argument (N, &argc, argv);
+    SET_MREL_VERSION = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   VERBOSE = VERBOSE2 = FALSE;
   if ((N = get_argument (argc, argv, "-v"))) {
Index: trunk/Ohana/src/relphot/src/setMrelCatalog.c
===================================================================
--- trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 34844)
+++ trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 34844)
@@ -0,0 +1,588 @@
+# include "relphot.h"
+
+# define SKIP_THIS_MEAS(REASON) {				\
+    measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;	\
+    if (measure) {					\
+      measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;		\
+    }								\
+    results->REASON ++;						\
+    continue; }
+
+# define SKIP_THIS_MEAS_STACK(REASON) {			\
+    measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;		\
+    continue; }
+
+int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) {
+
+  int Ns;
+  off_t k;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  for (k = 0; k < average[0].Nmeasure; k++) {
+    fprintf (stderr, "meas: %08x\n", measure[k].dbFlags);
+  }
+
+  for (Ns = 0; Ns < Nsecfilt; Ns++) {
+    fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
+  }
+  return (TRUE);
+}
+
+int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
+
+  off_t j;
+
+  liststats_setmode (&results->stats, STATMODE);
+  liststats_setmode (&results->apstats, STATMODE);
+  liststats_setmode (&results->kronstats, STATMODE);
+
+  SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
+
+  int isSetMrelFinal = (pass >= 0);
+
+  char *primaryCell = NULL;
+  if (isSetMrelFinal) {
+    ALLOCATE (primaryCell, char, DVO_MAX_PATH);
+  }
+
+  for (j = 0; j < catalog[Nc].Naverage; j++) {
+    off_t m = catalog[Nc].average[j].measureOffset;
+    setMrelAverageExposure (m, Nc, pass, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], &catalog[Nc].measureT[m], &catalog[Nc].found[Nsecfilt*j]);
+
+    // only apply Stack operation on setMrelFinal in first pass 
+    if (isSetMrelFinal && (pass == 0)) {
+      setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], primaryCell);
+      setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
+    }
+  }
+  if (primaryCell) free (primaryCell);
+  return (TRUE);
+}
+
+int setMrelAverageExposure (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) {
+
+  off_t k;
+  int N;
+  float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;
+
+  double *list      = results->list;
+  double *dlist     = results->dlist;
+  double *wlist     = results->wlist;
+  double *aplist    = results->aplist;
+  double *kronlist  = results->kronlist;
+  double *dkronlist = results->dkronlist;
+
+  StatType *stats = &results->stats;
+  StatType *apstats = &results->apstats;
+  StatType *kronstats = &results->kronstats;
+
+  int isSetMrelFinal = (pass >= 0);
+
+  // option for a test print
+  if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
+    fprintf (stderr, "test obj\n");
+    print_measure_set_alt (average, secfilt, measure);
+  }
+
+  int GoodPS1 = FALSE;
+  int Good2MASS = FALSE;
+  int Galaxy2MASS = FALSE;
+
+  int NextPS1 = 0;
+  int NpsfPS1 = 0;
+
+  int Ns;
+  for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+    int thisCode = photcodes[Ns][0].code;
+    int Nsec = GetPhotcodeNsec(thisCode);
+
+    /* calculate the average mag in this SEC photcode for a single star */
+
+    /* star/photcodes already calibrated */
+    if ( isSetMrelFinal && found[Nsec]) continue;  
+      
+    // skip bad stars
+    if (!isSetMrelFinal && (secfilt[Nsec].flags & STAR_BAD)) continue;
+
+    int Ncode = 0;
+    int Next = 0;
+    int haveSynth = FALSE;
+    int havePS1   = FALSE;
+
+    int forceSynth = FALSE;
+    int forceSynthEntry = -1;
+
+    int haveUbercal = FALSE;
+
+    int minUbercalDist = 1000;
+    
+    off_t meas = measureOffset;
+
+    N = 0;
+    for (k = 0; k < averageT[0].Nmeasure; k++, meas++) {
+
+      // skip measurements that do not match the current photcode
+      PhotCode *code = GetPhotcodebyCode (measureT[k].photcode);
+      if (!code) continue;
+      if (code->equiv != thisCode) { continue; }
+      Ncode ++;
+
+      if (measureT[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 
+
+      if (getImageEntry (meas, cat) < 0) {
+	// measurements without an image are either external reference photometry or
+	// data for which the associated image has not been loaded (probably because of
+	// overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
+	Mmos = Mgrid = 0;
+	Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value
+      } else {
+	Mcal  = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd);
+	if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
+	Mmos  = getMmos  (meas, cat);
+	if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
+	Mgrid = getMgrid (meas, cat);
+	if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
+      }
+
+      Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0]);
+      if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
+
+      int myUbercalDist = getUbercalDist(meas, cat);
+      minUbercalDist = MIN(minUbercalDist, myUbercalDist);
+
+      if (isSetMrelFinal) {
+	// special options for PS1 data
+	if ((measure[k].photcode >= 10000) && (measure[k].photcode <= 10500)) {
+	  // count the extended detections
+	  if (!isnan(measure[k].Map)) {
+	    float dMagAp = measure[k].M - measure[k].Map;
+	    float SigmaAp = hypot(0.1, 2.5*measure[k].dM);
+	    // XXX this is still quite ad hoc, but at least it:
+	    // (a) converges to 0.1 mag offset at the bright end
+	    // (b) converges to 0.5 mag offset at the faint end (dM = 0.2)
+	    if (dMagAp > SigmaAp) {
+	      Next ++;
+	      NextPS1 ++;
+	    } else {
+	      NpsfPS1 ++;
+	    }
+	  }
+	  havePS1 = TRUE;
+	}
+
+	// gpc1 stack data
+	if ((measure[k].photcode >= 11000) && (measure[k].photcode <= 11400)) {
+	  continue;
+	}
+
+	// count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
+	if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) {
+	  if (measure[k].photFlags & 0x00c00000) {
+	    Next ++;
+	    Galaxy2MASS = TRUE;
+	  }
+	  if (pass == 0) {
+	    if (measure[k].photFlags & 0x00000007) {
+	      Good2MASS = TRUE;
+	    } else {
+	      // detections without one of these bits should only be used in PASS_1
+	      SKIP_THIS_MEAS(Nbad);
+	    }
+	  }
+	}
+
+	// Blindly accepth the SYNTH mags if we are above saturation, otherwise, 
+	// ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars)
+	if ((measure[k].photcode >= 3001) && (measure[k].photcode <= 3005)) {
+	  // something of a hack: force object to use synth values if synth mags >>
+	  // saturation (3pi instrumental mags < -15)
+	  float MaxMagForceSynth = NAN;
+	  switch (measure[k].photcode) {
+	    case 3001:
+	      MaxMagForceSynth = 13.64;
+	      break;
+	    case 3002:
+	      MaxMagForceSynth = 13.76;
+	      break;
+	    case 3003:
+	      MaxMagForceSynth = 13.74;
+	      break;
+	    case 3004:
+	      MaxMagForceSynth = 12.94;
+	      break;
+	    case 3005:
+	      MaxMagForceSynth = 12.01;
+	      break;
+	  }
+	  if (measureT[k].M < MaxMagForceSynth) {
+	    forceSynth = TRUE;
+	    forceSynthEntry = N;
+	  } else {
+	    if (pass < 4) {
+	      SKIP_THIS_MEAS(Nbad);
+	    }
+	    haveSynth = TRUE;
+	  }
+	}
+
+	// Map (aplist) and Mkron (kronlist,dkronlist) are used to calculate mean mags per filter
+	float Map = PhotAper (&measure[k]);
+	aplist[N] = Map - Mcal - Mmos - Mgrid;
+
+	float Mkron = PhotKron (&measure[k]);
+	kronlist[N] = Mkron - Mcal - Mmos - Mgrid;
+	dkronlist[N] = measure[k].dMkron;
+      }
+
+      // dlist gives the error per measurement, wlist gives the weight
+      // we can modify the error and weight in a few ways:
+      // 1) MIN_ERROR guarantees a floor
+      // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
+      // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
+      // 4) some reference photcode of some kind can be specified as fixed and have a high weight
+      dlist[N] = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
+      wlist[N] = 1.0;
+      list[N]  = Msys - Mcal - Mmos - Mgrid;
+
+      // up-weight the ubercal values (or convergence can take a long time...)
+      if (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
+	haveUbercal = TRUE;
+	wlist[N] = 10.0;
+      }
+
+      // tie down reference photometry if the -refcode (code) option is selected
+      // eg, -refcode g_SDSS
+      // this probably makes no sense in the context of multifilter analysis
+      if (refPhotcode) {
+	if (code->code == refPhotcode->code) {
+	  wlist[N] = 100.0;
+	}
+      }
+      N++;
+    }
+
+    int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
+
+    // when performing the grid analysis, STAR_TOOFEW should be set to 1;
+    if (N < Nminmeas) { /* too few measurements */
+      // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, N, STAR_TOOFEW); 
+      secfilt[Nsec].flags |= ID_STAR_FEW;
+      if (Ncode == 0) {
+	results->Ncode ++;
+      } else {
+	results->Nfew ++;
+      }
+      continue;
+    } else {
+      secfilt[Nsec].flags &= ~ID_STAR_FEW;
+    }	
+
+    if (forceSynth) {
+      // use the single SYNTH value instead of the other mags here
+      myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error");
+      list[0]  = list[forceSynthEntry];
+      dlist[0] = dlist[forceSynthEntry];
+      wlist[0] = wlist[forceSynthEntry];
+      N = 1;
+    }
+    liststats (list, dlist, wlist, N, stats);
+
+    secfilt[Nsec].M  = stats->mean;
+    secfilt[Nsec].dM = stats->error;
+    secfilt[Nsec].Xm = (stats->Nmeas > 1) ? 100.0*log10(stats->chisq + 1e-4) : NAN_S_SHORT;
+
+    // when running -averages, we have no information about the images, so we cannot set this
+    if (minUbercalDist > -1) {
+      secfilt[Nsec].ubercalDist = minUbercalDist;
+    }
+
+    if (isSetMrelFinal) {
+      found[Nsec] = TRUE;
+
+      secfilt[Nsec].Mstdev = 1000.0*stats->sigma; // Mstdev is in millimags (not enough space for more precision)
+      secfilt[Nsec].Ncode = Ncode;
+      secfilt[Nsec].Nused = stats->Nmeas;
+
+      secfilt[Nsec].M_80 = 1000 * stats->Upper80;
+      secfilt[Nsec].M_20 = 1000 * stats->Lower20;
+
+      // NOTE : use the modified weight for apmags as well as psf mags
+      liststats (aplist, dlist, wlist, N, apstats);
+      secfilt[Nsec].Map  = apstats->mean; 
+
+      liststats (kronlist, dkronlist, wlist, N, kronstats);
+      secfilt[Nsec].Mkron  = kronstats->mean; 
+      secfilt[Nsec].dMkron = kronstats->error; 
+
+      // NOTE: for 2MASS measurements, Next should be 1, as should N
+      if ((Next > 0) && (Next > 0.5*N)) {
+	secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
+      }
+
+      switch (pass) {
+	case 0:
+	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_0;
+	  if (havePS1) GoodPS1 = TRUE;
+	  break;
+	case 1:
+	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_1;
+	  if (havePS1) GoodPS1 = TRUE;
+	  break;
+	case 2:
+	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_2;
+	  if (havePS1) GoodPS1 = TRUE;
+	  break;
+	case 3:
+	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_3;
+	  break;
+	case 4:
+	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_4;
+	  break;
+      }
+      if (haveSynth) {
+	secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
+      }	
+      if (havePS1) {
+	secfilt[Nsec].flags |= ID_SECF_HAS_PS1;
+      }	
+      if (haveUbercal) {
+	secfilt[Nsec].flags |= ID_SECF_USE_UBERCAL;
+      }	
+    }
+  }
+
+  if (isSetMrelFinal) {
+    if (pass == 0) {
+      DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
+      // we attempt to set a few flags here; reset those bits before trying:
+      average[0].flags &= ~flagBits;
+    }
+
+    if (NextPS1 && (NextPS1 > NpsfPS1)) {
+      average[0].flags |= ID_OBJ_EXT;
+    }
+    if (GoodPS1) {
+      average[0].flags |= ID_OBJ_GOOD;
+    }
+    if (Galaxy2MASS) {
+      average[0].flags |= ID_OBJ_EXT_ALT;
+    }
+    if (Good2MASS) {
+      average[0].flags |= ID_OBJ_GOOD_ALT;
+    }
+  }
+  return (TRUE);
+}
+
+// only apply Stack operation on setMrelFinal in first pass 
+int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell) {
+
+  off_t k, ID;
+
+  float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;
+
+  // set the name of the primary skycell (this is used in a strcmp to match the skycells in stack detections)
+  BoundaryTreePrimaryCell(primaryCell, average[0].R, average[0].D);
+
+  int NstackGood = 0;
+  int NstackSuspect = 0;
+  int haveStackObject = FALSE;
+
+  int Ns;
+  for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+    int thisCode = photcodes[Ns][0].code;
+    int Nsec = GetPhotcodeNsec(thisCode);
+
+    /* calculate the average mag in this SEC photcode for a single star */
+
+    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;
+    off_t stackCenterMeasureMin = -1;
+
+    float stackPrimaryOffsetMin = 1e9;
+    off_t stackPrimaryMeasureMin = -1;
+
+    int isBad = FALSE;
+    int isSuspect = FALSE;
+
+    off_t meas = measureOffset;
+    for (k = 0; k < average[0].Nmeasure; k++, meas++) {
+
+      // skip measurements that do not match the current photcode
+      PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
+      if (!code) continue;
+      if (code->equiv != thisCode) { continue; }
+
+      // clear this bit for all measurements
+      measure[k].dbFlags &= ~ID_MEAS_STACK_PRIMARY;
+
+      // only examine gpc1 stack data
+      // XXX this is absurdly hardwired (along with several photcode tests)
+      if (measure[k].photcode < 11000) continue;
+      if (measure[k].photcode > 11400) continue;
+
+      haveStack = TRUE;
+      haveStackObject = TRUE;
+      
+      if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); 
+
+      if (getImageEntry (meas, cat) < 0) continue;
+      // measurements without an image are either external reference photometry or
+      // data for which the associated image has not been loaded (probably because of
+      // overlaps).  we only want measurements associated with stack images in this loop
+
+      Mcal  = getMcal_alt  (meas, cat, flatcorr, measure[k].Xccd, measure[k].Yccd);
+      if (isnan(Mcal))  SKIP_THIS_MEAS_STACK(Ncal);
+      Mmos  = getMmos  (meas, cat);
+      if (isnan(Mmos))  SKIP_THIS_MEAS_STACK(Nmos);
+      Mgrid = getMgrid (meas, cat);
+      if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid);
+
+      Msys = PhotSys (&measure[k], &average[0], &secfilt[0]);
+      if (isnan(Msys)) SKIP_THIS_MEAS_STACK(Nsys);
+
+      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 (meas, cat, primaryCell)) {
+	float stackPrimaryOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID);
+	if (stackPrimaryOffset < stackPrimaryOffsetMin) {
+	  stackPrimaryOffsetMin = stackPrimaryOffset;
+	  // stackPrimaryIDmin = stackImageID;
+	  stackPrimaryMeasureMin = k;
+	}
+      }
+
+      // get the center distance for the generic case:
+      float stackCenterOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID);
+      if (stackCenterOffset < stackCenterOffsetMin) {
+	stackCenterOffsetMin = stackCenterOffset;
+	// stackCenterIDmin = stackImageID;
+	stackCenterMeasureMin = k;
+      }
+
+      // data quality assessment
+      isBad |= (measure[k].photFlags & code->photomBadMask);
+      isBad |= (measure[k].psfQF < 0.85);
+      isBad |= isnan(measure[k].psfQF);
+      isBad |= measure[k].dM > 0.2; // S/N < 5.0
+
+      isSuspect |= (measure[k].photFlags & code->photomPoorMask);
+      isSuspect |= (measure[k].psfQFperf < 0.85);
+    }
+
+    if (!haveStack) continue;
+
+    if (!isSuspect && !isBad) {
+      NstackGood ++;
+    }
+    if (isSuspect && !isBad) {
+      NstackSuspect ++;
+    }
+
+    // measurements which are bad will not have a valid stack entry and are skipped
+    k = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;
+    if (k < 0) continue;
+
+
+    // we are now populating stackDetectID not stack Image ID in secfilt
+    // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
+    ID = measure[k].extID; // for the stack, this is the stackDetectID
+
+    // get the zero point for the selected image
+    float zp = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (Mcal + Mmos + Mgrid);
+
+    // 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?)
+    secfilt[Nsec].FluxPSF   = zpFactor * measure[k].FluxPSF;  
+    secfilt[Nsec].dFluxPSF  = zpFactor * measure[k].dFluxPSF; 
+    secfilt[Nsec].FluxKron  = zpFactor * measure[k].FluxKron; 
+    secfilt[Nsec].dFluxKron = zpFactor * measure[k].dFluxKron;
+
+    secfilt[Nsec].stackDetectID = ID;
+
+    measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
+    if (stackPrimaryMeasureMin >= 0) {
+      measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
+    }
+
+    secfilt[Nsec].flags |= ID_SECF_HAS_STACK;
+  } // Nsecfilt loop
+
+  if (NstackGood >= 2) {
+    average[0].flags |= ID_OBJ_GOOD_STACK;
+  } else if (NstackGood + NstackSuspect >= 2) {
+    average[0].flags |= ID_OBJ_SUSPECT_STACK;
+  } else if (haveStackObject) {
+    average[0].flags |= ID_OBJ_BAD_STACK;
+  }
+  return (TRUE);
+}
+
+int setGlobalObjStats (Average *average, Measure *measure) {
+
+  off_t k;
+
+  float psfQFmax     = 0.0;
+  float psfQFperfmax = 0.0;
+  float stargalmax   = 0.0;
+
+  for (k = 0; k < average[0].Nmeasure; k++) {
+
+    if (isfinite(measure[k].psfQF)) {
+      psfQFmax     = MAX (psfQFmax,     measure[k].psfQF);
+    }
+    if (isfinite(measure[k].psfQFperf)) {
+      psfQFperfmax = MAX (psfQFperfmax, measure[k].psfQFperf);
+    }
+
+    // only count psfQF > 0.85 here
+    if (isfinite(measure[k].extNsigma) && (measure[k].psfQF > 0.85)) {
+      stargalmax   = MAX (stargalmax,   measure[k].extNsigma);
+    }
+  }
+  average[0].psfQF     = psfQFmax;
+  average[0].psfQFperf = psfQFperfmax;
+  average[0].stargal   = stargalmax;
+
+  return (TRUE);
+}
Index: trunk/Ohana/src/uniphot/src/initialize_setposangle.c
===================================================================
--- trunk/Ohana/src/uniphot/src/initialize_setposangle.c	(revision 34843)
+++ trunk/Ohana/src/uniphot/src/initialize_setposangle.c	(revision 34844)
@@ -6,4 +6,9 @@
     fprintf (stderr, "    -v : verbose mode\n");
     fprintf (stderr, "    -update : actually write results to detections tables\n");
+    fprintf (stderr, "    -parallel : run in parallel mode\n");
+    fprintf (stderr, "    -h     : this help list\n");
+    fprintf (stderr, "    -help  : this help list\n");
+    fprintf (stderr, "    --h    : this help list\n");
+    fprintf (stderr, "    --help : this help list\n");
     fprintf (stderr, "    Note that the dvo db can be specified by -D CATDIR (directory)\n");
     exit (2);
Index: trunk/Ohana/src/uniphot/src/update_catalog_setposangle.c
===================================================================
--- trunk/Ohana/src/uniphot/src/update_catalog_setposangle.c	(revision 34843)
+++ trunk/Ohana/src/uniphot/src/update_catalog_setposangle.c	(revision 34844)
@@ -25,9 +25,10 @@
       Mosaic *mosaic = getMosaicForImage(seq);
       Coords *coords = &image[seq].coords;
+      Coords *mosaicCoords = mosaic ? &mosaic->coords : NULL;
 
       double Xccd = catalog[0].measure[m].Xccd;
       double Yccd = catalog[0].measure[m].Yccd;
       
-      setposangle_local_astrometry (&posAngle, &pltScale, Xccd, Yccd, &mosaic->coords, coords);
+      setposangle_local_astrometry (&posAngle, &pltScale, Xccd, Yccd, mosaicCoords, coords);
 
       catalog[0].measure[m].posangle = ToShortDegrees(posAngle);
@@ -52,11 +53,26 @@
   // calculate the astrometry for the coordinate of interest
   XY_to_LM (&Lx, &Mx, x,       y,       coords);
-  XY_to_LM (&Po, &Qo, Lx,      Mx,      mosaic);
+  if (mosaic) {
+    XY_to_LM (&Po, &Qo, Lx,      Mx,      mosaic);
+  } else {
+    Po = Lx;
+    Qo = Mx;
+  }
 
   XY_to_LM (&Lx, &Mx, x + 1.0, y,       coords);
-  XY_to_LM (&Px, &Qx, Lx,      Mx,      mosaic);
+  if (mosaic) {
+    XY_to_LM (&Px, &Qx, Lx,      Mx,      mosaic);
+  } else {
+    Px = Lx;
+    Qx = Mx;
+  }
 
   XY_to_LM (&Lx, &Mx, x,       y + 1.0, coords);
-  XY_to_LM (&Py, &Qy, Lx,      Mx,      mosaic);
+  if (mosaic) {
+    XY_to_LM (&Py, &Qy, Lx,      Mx,      mosaic);
+  } else {
+    Py = Lx;
+    Qy = Mx;
+  }
 
   // XXX units for the resulting Tangent Plane coordinates??