Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c	(revision 33609)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c	(revision 33610)
@@ -153,4 +153,27 @@
 }
 
+int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
+
+  off_t k;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  off_t m = average[0].measureOffset;
+
+  for (k = 0; k < average[0].Nmeasure; k++, m++) {
+    fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
+  }
+
+  int Ns;
+  for (Ns = 0; Ns < Nsecfilt; Ns++) {
+    fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
+  }
+  return 1;
+}
+
+# define MARK_SKIP_MEAS \
+  catalog[i].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \
+  catalog[i].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
+
 // setMrel and setMrelOutput are extremely similar, but have slightly different implications:
 // * setMrel uses the internal Tiny structures only
@@ -190,4 +213,9 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
 
+      if (FALSE && (catalog[i].average[j].objID == 0x46a4) && (catalog[i].average[j].catID == 0xf40e)) {
+	fprintf (stderr, "test obj\n");
+	print_measure_set (&catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt], catalog[i].measure);
+      }
+
       int minUbercalDist = 1000;
 
@@ -205,4 +233,5 @@
 	int haveSynth = FALSE;
 	int haveUbercal = FALSE;
+	int isRefPhot = FALSE;
 
 	N = 0;
@@ -221,15 +250,30 @@
 	    // these detetions have no image (eg, ref values such as 2MASS)
 	    Mcal = Mmos = Mgrid = 0;
+	    isRefPhot = TRUE;
 	  } else {
 	    Mcal  = getMcal  (m, i, flatcorr, catalog);
-	    if (isnan(Mcal)) continue;
+	    if (isnan(Mcal)) {
+	      MARK_SKIP_MEAS;
+	      continue;
+	    }
 	    Mmos  = getMmos  (m, i);
-	    if (isnan(Mmos)) continue;
+	    if (isnan(Mmos)) {
+	      MARK_SKIP_MEAS;
+	      continue;
+	    }
 	    Mgrid = getMgrid (m, i);
-	    if (isnan(Mgrid)) continue;
+	    if (isnan(Mgrid)) {
+	      MARK_SKIP_MEAS;
+	      continue;
+	    }
 	  }
 
 	  Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
-	  if (isnan(Msys)) { continue; }
+	  if (isnan(Msys)) { 
+	    if (!isRefPhot) {
+	      MARK_SKIP_MEAS;
+	    }
+	    continue; 
+	  }
 	  list[N] = Msys - Mcal - Mmos - Mgrid;
 
@@ -255,4 +299,5 @@
 	    if ((pass == 0) && !(catalog[i].measure[m].photFlags & 0x00000007)) {
 	      // detections without one of these bits should only be used in PASS_1
+	      MARK_SKIP_MEAS;
 	      continue;
 	    }
@@ -260,5 +305,8 @@
 	  // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars)
 	  if ((catalog[i].measure[m].photcode >= 3001) && (catalog[i].measure[m].photcode <= 3005)) {
-	    if (pass < 4) continue;
+	    if (pass < 4) {
+	      MARK_SKIP_MEAS;
+	      continue;
+	    }
 	    haveSynth = TRUE;
 	  }
Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/bcatalog.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/bcatalog.c	(revision 33609)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/bcatalog.c	(revision 33610)
@@ -1,3 +1,5 @@
 # include "relphot.h"
+
+int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog);
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
@@ -212,5 +214,5 @@
   // limit the total number of stars in the catalog
   if (MaxDensityUse) {
-    LimitDensityCatalog (subcatalog, catalog);
+    LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
   }
 
@@ -224,5 +226,21 @@
 }
 
-int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
+// sort by decreasing Nmeasure (X)
+void sort_by_Nmeasure (int *X, off_t *Y, off_t N) {
+
+# define SWAPFUNC(A,B){ int tmpI; off_t tmpT;	\
+  tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \
+  tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \
+}
+# define COMPARE(A,B)(X[A] > X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) {
 
   Catalog tmpcatalog;
@@ -232,8 +250,8 @@
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
-  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
-  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
-  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
-  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
+  gfits_scan (&oldcatalog[0].header, "RA0",  "%lf", 1, &Rmin);
+  gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin);
+  gfits_scan (&oldcatalog[0].header, "RA1",  "%lf", 1, &Rmax);
+  gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax);
 
   if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
@@ -255,15 +273,13 @@
 
   // we are going to select Nmax entries by generating a random-sorted index list
-  off_t *index, tmp, i, j, ave;
+  int *value;
+  off_t *index, i, j, ave;
   ALLOCATE (index, off_t, Naverage);
+  ALLOCATE (value, int, Naverage);
   for (i = 0; i < Naverage; i++) {
     index[i] = i;
-  }
-  for (i = 0; i < Naverage; i++) {
-    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
-    tmp = index[j];
-    index[j] = index[i];
-    index[i] = tmp;
-  }
+    value[i] = subcatalog[0].averageT[i].Nmeasure;
+  }
+  sort_by_Nmeasure (value, index, Naverage);
 
   // count the number of measurements this selection will yield
@@ -298,4 +314,101 @@
 
   if (VERBOSE) {
+    char *basename = filebasename (oldcatalog[0].filename);
+    fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", 
+	     Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
+    free (basename);
+  }
+
+  free (index);
+  free (value);
+  free (subcatalog[0].averageT);
+  free (subcatalog[0].measureT);
+  free (subcatalog[0].secfilt);
+
+  subcatalog[0].averageT = tmpcatalog.averageT;
+  subcatalog[0].measureT = tmpcatalog.measureT;
+  subcatalog[0].secfilt = tmpcatalog.secfilt;
+  subcatalog[0].Naverage = Nmax;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;
+
+  return (TRUE);
+}
+
+int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
+
+  Catalog tmpcatalog;
+
+  double Rmin, Rmax, Dmin, Dmax;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
+  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
+  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
+  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
+
+  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
+
+  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
+  assert (AREA > 0);
+
+  off_t Nmax = MaxDensityValue * AREA;
+  if (subcatalog[0].Naverage <= Nmax) {
+    if (VERBOSE) {
+      fprintf (stderr, "subcatalog has less than the max density\n");
+    }
+    return (TRUE);
+  }
+
+  off_t Naverage = subcatalog[0].Naverage;
+
+  // select a random subset of Nmax stars from subcatalog using Fisher-Yates
+
+  // we are going to select Nmax entries by generating a random-sorted index list
+  off_t *index, tmp, i, j, ave;
+  ALLOCATE (index, off_t, Naverage);
+  for (i = 0; i < Naverage; i++) {
+    index[i] = i;
+  }
+  for (i = 0; i < Naverage; i++) {
+    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
+    tmp = index[j];
+    index[j] = index[i];
+    index[i] = tmp;
+  }
+
+  // count the number of measurements this selection will yield
+  off_t NMEASURE = 0;
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    NMEASURE += subcatalog[0].averageT[ave].Nmeasure;
+  }
+
+  // allocate the output data 
+  ALLOCATE (tmpcatalog.averageT, AverageTiny, Nmax);
+  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
+  ALLOCATE (tmpcatalog.secfilt,  SecFilt, Nmax * Nsecfilt);
+
+  off_t Nmeasure = 0;
+
+  // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links)
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    tmpcatalog.averageT[i] = subcatalog[0].averageT[ave];
+    tmpcatalog.averageT[i].measureOffset = Nmeasure;
+    for (j = 0; j < tmpcatalog.averageT[i].Nmeasure; j++) {
+      off_t offset = subcatalog[0].averageT[ave].measureOffset + j;
+      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
+      tmpcatalog.measureT[Nmeasure].averef = i;
+      Nmeasure ++;
+    }
+    for (j = 0; j < Nsecfilt; j++) {
+      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
+    }
+  }
+
+  if (VERBOSE) {
     fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", 
 	     Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  catalog[0].filename);
Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/relphot.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/relphot.c	(revision 33609)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/relphot.c	(revision 33610)
@@ -230,20 +230,21 @@
       MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
     }
+
+    /* at this point, we have correct cal coeffs in the image/mosaic structures */
+    for (i = 0; i < Ncatalog; i++) {
+      free_tiny_values (&catalog[i]);
+      dvo_catalog_free (&catalog[i]);
+    }
+    freeImageBins (Ncatalog, TRUE);
+    freeMosaicBins (Ncatalog, TRUE);
+    freeGridBins (Ncatalog);
+
   } // (NLOOP > 0) : this loop determines the offsets per chip 
 
-  // only change the real database if -update is requested
+  reload_images (&db);
+
+  // only change the real database files if -update is requested
   if (!UPDATE) exit (0);
   
-  reload_images (&db);
-
-  /* at this point, we have correct cal coeffs in the image/mosaic structures */
-  for (i = 0; i < Ncatalog; i++) {
-    free_tiny_values (&catalog[i]);
-    dvo_catalog_free (&catalog[i]);
-  }
-  freeImageBins (Ncatalog, TRUE);
-  freeMosaicBins (Ncatalog, TRUE);
-  freeGridBins (Ncatalog);
-
   /* load catalog data from region files, update Mrel include all data */
   reload_catalogs (skylist, flatcorr, 0, NULL);
Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c	(revision 33609)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c	(revision 33610)
@@ -123,4 +123,25 @@
 // non-outliers with Imag range
 
+// int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
+// 
+//   off_t k;
+// 
+//   int Nsecfilt = GetPhotcodeNsecfilt ();
+// 
+//   off_t m = average[0].measureOffset;
+// 
+//   for (k = 0; k < average[0].Nmeasure; k++, m++) {
+//     fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
+//   }
+// 
+//   int Ns;
+//   for (Ns = 0; Ns < Nsecfilt; Ns++) {
+//     fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
+//   }
+//   return 1;
+// }
+
+int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
+
 void skip_measurements (Catalog *catalog, int pass, FlatCorrectionTable *flatcorr) {
 
@@ -137,9 +158,18 @@
   
   /* allow measures marked as outliers (POOR) and off image region (AREA) */
-  if (pass >= 2) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
+  if (pass >= 2) {
+    MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
+  } else {
+    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
+  }
 
   /* mark measures which should be ignored on this pass */
   for (i = 0; i < catalog[0].Naverage; i++) {
     Ntot += catalog[0].average[i].Nmeasure;
+
+    if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) {
+      fprintf (stderr, "test obj\n");
+      print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure);
+    }
 
     // mark measurements for each secfilt separately
@@ -173,15 +203,21 @@
 	// clear SKIP for all measures at first 
 	catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-
-	// skip measurements from BAD images and mosaics 
-	// do NOT skip measurements without a matching image 
-	if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
-	if (isnan(getMmos (m, 0))) goto skip;
-
-	// PASS 4 : skip measurements by inst mag limit 
+	catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
+
+	// skip measurements from BAD images and mosaics (not REF mags)
+	// do NOT skip measurements without a matching image (REF mags)
+	off_t Nim = getImageEntry (m, 0);
+	if (Nim > -1) {
+	  if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
+	  if (isnan(getMmos (m, 0))) goto skip;
+	}
+
+	// PASS 4 : skip measurements by inst mag limit (not REF mags)
 	if ((pass < 4) && ImagSelect) {
-	  mag = PhotInst (&catalog[0].measure[m]);
-	  if (mag < ImagMin) goto skip;
-	  if (mag > ImagMax) goto skip;
+	  if (Nim > -1) {
+	    mag = PhotInst (&catalog[0].measure[m]);
+	    if (mag < ImagMin) goto skip;
+	    if (mag > ImagMax) goto skip;
+	  }
 	}