Index: trunk/Ohana/src/getstar/include/getstar.h
===================================================================
--- trunk/Ohana/src/getstar/include/getstar.h	(revision 20184)
+++ trunk/Ohana/src/getstar/include/getstar.h	(revision 20266)
@@ -27,4 +27,8 @@
 int  SKY_DEPTH;
 
+int MaxDensityUse;
+float MaxDensityValue;
+int MinMagUse;
+float MinMagValue;
 int  MagLimitUse;
 float MagLimitValue;
Index: trunk/Ohana/src/getstar/src/args.c
===================================================================
--- trunk/Ohana/src/getstar/src/args.c	(revision 20184)
+++ trunk/Ohana/src/getstar/src/args.c	(revision 20266)
@@ -44,4 +44,24 @@
     remove_argument (N, &argc, argv);
     MagLimitValue = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  MaxDensityUse = FALSE;
+  if ((N = get_argument (argc, argv, "-max-density"))) {
+    if (MagLimitUse) {
+      fprintf (stderr, "-maglim and -max-density are mutually exclusive; ignoring -maglim");
+      // MagLimitValue will be overridded by a density-based limit
+    }
+    MaxDensityUse = TRUE;
+    remove_argument (N, &argc, argv);
+    MaxDensityValue = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  MinMagUse = FALSE;
+  if ((N = get_argument (argc, argv, "-minmag"))) {
+    MinMagUse = TRUE;
+    remove_argument (N, &argc, argv);
+    MinMagValue = atof (argv[N]);
     remove_argument (N, &argc, argv);
   }
Index: trunk/Ohana/src/getstar/src/select_by_region.c
===================================================================
--- trunk/Ohana/src/getstar/src/select_by_region.c	(revision 20184)
+++ trunk/Ohana/src/getstar/src/select_by_region.c	(revision 20266)
@@ -6,5 +6,5 @@
   int i, j, n, Nm, offset, m, Nsecfilt, code, Nsec;
   int Nave, NAVE, Nmeas, NMEAS;
-  double R, D;
+  double R, D, AREA;
   float mag;
 
@@ -37,4 +37,67 @@
   NMEAS = output[0].Nmeasure + 1000;
   REALLOCATE (output[0].measure, Measure, NMEAS);
+
+  // if MaxDensityUse is set, we need to determine the MagLimitValue for this catalog by
+  // examining the magnitude distribution.  generate a magnitude histogram
+  if (MaxDensityUse) {
+    int bin, Nmax, Nbins, Nsum, *Nmag;
+    float MagMin, MagMax, dMag;
+    double Rmin, Rmax, Dmin, Dmax;
+
+    MagLimitUse = TRUE; // we will use the MagLimit implied by the mag distribution below
+
+    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 (VERBOSE) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
+
+    AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));  
+    assert (AREA > 0);
+
+    Nmax = MaxDensityValue * AREA;
+
+    dMag = 0.1;
+    MagMin = 0.0;
+    MagMax = 32.0;
+    Nbins = (MagMax - MagMin) / dMag;
+    ALLOCATE (Nmag, int, Nbins);
+    bzero (Nmag, Nbins * sizeof(int));
+
+    for (i = 0; i < catalog[0].Naverage; i++) {
+      n = catalog[0].measure[i].averef;
+      mag = NAN;
+      if (Nsec != -1) {
+	mag = catalog[0].secfilt[i*Nsecfilt + Nsec].M;
+      } else {
+	offset = catalog[0].average[i].measureOffset;
+	for (m = 0; m < catalog[0].average[i].Nmeasure; m++) {
+	  if (catalog[0].measure[offset + m].photcode == code) {
+	    mag = PhotRel (&catalog[0].measure[offset + m], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
+	    break;
+	  }
+	}
+      }
+      if (isnan(mag)) continue;
+
+      bin = (mag - MagMin) / dMag;
+      if (bin < 0) continue;
+      if (bin >= Nbins) continue;
+
+      Nmag[bin] ++;
+    }
+      
+    if (!MinMagUse) MinMagValue = MagMin;
+    bin = (MinMagValue - MagMin) / dMag;
+    bin = MIN (MAX (0, bin), Nbins);
+
+    Nsum = 0;
+    for (i = bin; (i < Nbins) && (Nsum < Nmax); i++) {
+      Nsum += Nmag[i];
+    }
+    MagLimitValue = i*dMag + MagMin;
+    if (VERBOSE) fprintf (stderr, "using %d stars in mag range %f - %f\n", Nsum, MinMagValue, MagLimitValue);
+  }
 
   for (i = 0; i < catalog[0].Naverage; i++) {