Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/include/relastro.h
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/include/relastro.h	(revision 33566)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/include/relastro.h	(revision 33567)
@@ -417,5 +417,6 @@
 
 int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit);
-int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog);
+int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *catalog);
+int LimitDensityCatalog_RandomSample (Catalog *subcatalog, Catalog *catalog);
 
 int initializeConstraints();
Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c	(revision 33566)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c	(revision 33567)
@@ -114,5 +114,5 @@
   // limit the total number of stars in the catalog
   if (MaxDensityUse) {
-    LimitDensityCatalog (subcatalog, catalog);
+    LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
   }
 
@@ -123,6 +123,22 @@
 }
 
+// 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
+
+}
+
 /* this version does NOT use AverageTiny, MeasureTiny */ 
-int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
+int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *catalog) {
 
   Catalog tmpcatalog;
@@ -154,16 +170,14 @@
   // 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;
+  // we are going to select Nmax entries by choosing the brightest objects
+  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] = catalog[0].average[i].Nmeasure;
+  }
+  sort_by_Nmeasure (value, index, Naverage);
 
   // count the number of measurements this selection will yield
@@ -202,4 +216,6 @@
   }
 
+  free (index);
+  free (value);
   free (subcatalog[0].average);
   free (subcatalog[0].measureT);
@@ -217,2 +233,96 @@
 }
 
+/* this version does NOT use AverageTiny, MeasureTiny */ 
+int LimitDensityCatalog_RandomSample (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].average[ave].Nmeasure;
+  }
+
+  // allocate the output data 
+  ALLOCATE (tmpcatalog.average,  Average,     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.average[i] = subcatalog[0].average[ave];
+    tmpcatalog.average[i].measureOffset = Nmeasure;
+    for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
+      off_t offset = subcatalog[0].average[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);
+  }
+
+  free (subcatalog[0].average);
+  free (subcatalog[0].measureT);
+  free (subcatalog[0].secfilt);
+
+  subcatalog[0].average = tmpcatalog.average;
+  subcatalog[0].measureT = tmpcatalog.measureT;
+  subcatalog[0].secfilt = tmpcatalog.secfilt;
+  subcatalog[0].Naverage = Nmax;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt;
+
+  return (TRUE);
+}
+