Index: /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/relphot/src/StarOps.c
===================================================================
--- /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/relphot/src/StarOps.c	(revision 24979)
+++ /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/relphot/src/StarOps.c	(revision 24980)
@@ -8,5 +8,5 @@
 
   int i, j;
-  
+
   Nmax = 0;
   for (i = 0; i < Ncatalog; i++) {
@@ -18,5 +18,5 @@
   ALLOCATE (list, double, MAX (1, Nmax));
   ALLOCATE (dlist, double, MAX (1, Nmax));
-}  
+}
 
 float getMrel (Catalog *catalog, int meas, int cat) {
@@ -26,6 +26,6 @@
 
   ave = catalog[cat].measure[meas].averef;
-  if (catalog[cat].average[ave].flags & STAR_BAD) return (NAN);  
- 
+  if (catalog[cat].average[ave].flags & STAR_BAD) return (NAN);
+
   value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M;
   return (value);
@@ -44,44 +44,44 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      if (catalog[i].average[j].flags & STAR_BAD) continue;
       m = catalog[i].average[j].measureOffset;
 
       N = 0;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	if (catalog[i].measure[m].dbFlags & MEAS_BAD) {
-	  Nbad ++;
-	  continue;
-	}
-	// XXX allow REF stars (no Image Entry) to be included in the calculation this
-	// should be optionally set, and should allow for REF stars to be downweighted by
-	// more than their reported errors.  how such info is carried is unclear...
-	if (getImageEntry (m, i) < 0) {
-	  Mcal = Mmos = Mgrid = 0;
-	} else {
-	  Mcal  = getMcal  (m, i);
-	  if (isnan(Mcal)) { 
-	    Ncal ++;
-	    continue;
-	  }
-	  Mmos  = getMmos  (m, i);
-	  if (isnan(Mmos)) {
-	    Nmos ++;
-	    continue;
-	  }
-	  Mgrid = getMgrid (m, i);
-	  if (isnan(Mgrid)) {
-	    Ngrid++;
-	    continue;
-	  }
-	}
-
-	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-	if (isnan(Msys)) {
-	  Nsys++;
-	  continue;
-	}
-	list[N] = Msys - Mcal - Mmos - Mgrid;
-	dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-	N++;
+        if (catalog[i].measure[m].dbFlags & MEAS_BAD) {
+          Nbad ++;
+          continue;
+        }
+        // XXX allow REF stars (no Image Entry) to be included in the calculation this
+        // should be optionally set, and should allow for REF stars to be downweighted by
+        // more than their reported errors.  how such info is carried is unclear...
+        if (getImageEntry (m, i) < 0) {
+          Mcal = Mmos = Mgrid = 0;
+        } else {
+          Mcal  = getMcal  (m, i);
+          if (isnan(Mcal)) {
+            Ncal ++;
+            continue;
+          }
+          Mmos  = getMmos  (m, i);
+          if (isnan(Mmos)) {
+            Nmos ++;
+            continue;
+          }
+          Mgrid = getMgrid (m, i);
+          if (isnan(Mgrid)) {
+            Ngrid++;
+            continue;
+          }
+        }
+
+        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+        if (isnan(Msys)) {
+          Nsys++;
+          continue;
+        }
+        list[N] = Msys - Mcal - Mmos - Mgrid;
+        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+        N++;
       }
 
@@ -89,9 +89,9 @@
 
       if (N <= STAR_TOOFEW) { /* too few measurements */
-	catalog[i].average[j].flags |= ID_STAR_FEW;
-	Nfew ++;
+        catalog[i].average[j].flags |= ID_STAR_FEW;
+        Nfew ++;
       } else {
-	catalog[i].average[j].flags &= ~ID_STAR_FEW;
-      }	
+        catalog[i].average[j].flags &= ~ID_STAR_FEW;
+      }
 
       liststats (list, dlist, N, &stats);
@@ -128,28 +128,28 @@
 
       /* skip stars already calibrated */
-      if (catalog[i].found[j]) continue;  
+      if (catalog[i].found[j]) continue;
 
       N = 0;
       m = catalog[i].average[j].measureOffset;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
-	// XXX allow REF stars (no Image Entry) to be included in the calculation this
-	// should be optionally set, and should allow for REF stars to be downweighted by
-	// more than their reported errors.  how such info is carried is unclear...
-	if (getImageEntry (m, i) < 0) {
-	  Mcal = Mmos = Mgrid = 0;
-	} else {
-	  Mcal  = getMcal  (m, i);
-	  if (isnan(Mcal)) continue;
-	  Mmos  = getMmos  (m, i);
-	  if (isnan(Mmos)) continue;
-	  Mgrid = getMgrid (m, i);
-	  if (isnan(Mgrid)) continue;
-	}
-
-	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-	list[N] = Msys - Mcal - Mmos - Mgrid;
-	dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-	N++;
+        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
+        // XXX allow REF stars (no Image Entry) to be included in the calculation this
+        // should be optionally set, and should allow for REF stars to be downweighted by
+        // more than their reported errors.  how such info is carried is unclear...
+        if (getImageEntry (m, i) < 0) {
+          Mcal = Mmos = Mgrid = 0;
+        } else {
+          Mcal  = getMcal  (m, i);
+          if (isnan(Mcal)) continue;
+          Mmos  = getMmos  (m, i);
+          if (isnan(Mmos)) continue;
+          Mgrid = getMgrid (m, i);
+          if (isnan(Mgrid)) continue;
+        }
+
+        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+        list[N] = Msys - Mcal - Mmos - Mgrid;
+        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+        N++;
       }
       if (N < 1) continue;
@@ -195,33 +195,40 @@
       for (Ns = 0; Ns < Nsecfilt; Ns++) {
 
-	code = GetPhotcodebyNsec (Ns);
-	Nc = code[0].code;
-	
-	N = 0;
-	m = catalog[i].average[j].measureOffset;
-	for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	  if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
-	  if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue;
-
-	  Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-	  if (isnan(Msys)) continue;
-
-	  // XXX this is a hack for the 2MASS search; better to save an average value?
-	  if (catalog[i].measure[m].psfQual < 0.85) continue;
-
-	  list[N] = Msys - catalog[i].measure[m].Mcal;
-	  dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-	  N++;
-	}
-	if (N < 1) continue;
-
-	liststats (list, dlist, N, &stats);
-
-	/* use sigma or error in dM for output? */
-	catalog[i].secfilt[Nsecfilt*j+Ns].M  = stats.mean;
-	catalog[i].secfilt[Nsecfilt*j+Ns].dM = MAX (stats.error, stats.sigma);
-	catalog[i].secfilt[Nsecfilt*j+Ns].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
-	catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N;
-	catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas;
+        code = GetPhotcodebyNsec (Ns);
+        Nc = code[0].code;
+
+        N = 0;
+        m = catalog[i].average[j].measureOffset;
+        for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+          if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
+          if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue;
+
+          Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+          if (isnan(Msys)) continue;
+
+          // XXX this is a hack for the 2MASS search; better to save an average value?
+          if (catalog[i].measure[m].psfQual < 0.85) continue;
+
+          list[N] = Msys - catalog[i].measure[m].Mcal;
+          dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+          N++;
+        }
+        if (N < 1) continue;
+
+        liststats (list, dlist, N, &stats);
+
+        /* use sigma or error in dM for output? */
+        // XXXJester: why not report *both* the RMS and the sigma, where one is the error-propagated standard
+        // error on the mean, and the other is the RMS of the values about the mean [with or without scaling
+        // by sqrt(N)]? If the two diverge by a lot, we have some un-understood source of systematic errors
+        // that could be QA-ed by looking for this divergence, but will otherwise only be found if people
+        // think actively about the expectations for what the errors should be, in effect doing this
+        // calculation by hand.
+        catalog[i].secfilt[Nsecfilt*j+Ns].M  = stats.mean;
+        // XXXJester: Neither of stats.error,stats.sigma must be NaN for the MAX() comparison to be meaningful!
+        catalog[i].secfilt[Nsecfilt*j+Ns].dM = MAX (stats.error, stats.sigma);
+        catalog[i].secfilt[Nsecfilt*j+Ns].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
+        catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N;
+        catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas;
       }
     }
@@ -247,12 +254,12 @@
       m = catalog[i].average[j].measureOffset;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
-	Mcal  = getMcal  (m, i);
-	if (isnan(Mcal)) continue;
-	Mmos  = getMmos  (m, i);
-	if (isnan(Mmos)) continue;
-	Mgrid = getMgrid (m, i);
-	if (isnan(Mgrid)) continue;
-	catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
+        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
+        Mcal  = getMcal  (m, i);
+        if (isnan(Mcal)) continue;
+        Mmos  = getMmos  (m, i);
+        if (isnan(Mmos)) continue;
+        Mgrid = getMgrid (m, i);
+        if (isnan(Mgrid)) continue;
+        catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
       }
     }
@@ -273,5 +280,5 @@
   /* find Xm median -> ChiSq lim must be > median */
   for (i = Ntot = 0; i < Ncatalog; i++) {
-    Ntot += catalog[i].Naverage; 
+    Ntot += catalog[i].Naverage;
   }
   ALLOCATE (xlist, double, Ntot);
@@ -290,5 +297,5 @@
     }
   }
-  
+
   initstats ("MEAN");
   liststats (xlist, dlist, Ntot, &stats);
@@ -306,8 +313,8 @@
       mark = (dM > MaxScatter) || (Xm == NAN_S_SHORT) || (Chisq > MaxChisq);
       if (mark) {
-	catalog[i].average[j].flags |= ID_STAR_POOR;
-	Ndel ++;
+        catalog[i].average[j].flags |= ID_STAR_POOR;
+        Ndel ++;
       } else {
-	catalog[i].average[j].flags &= ~ID_STAR_POOR;
+        catalog[i].average[j].flags &= ~ID_STAR_POOR;
       }
       Nave ++;
@@ -341,5 +348,5 @@
   ALLOCATE (ilist, int, Nmax);
   ALLOCATE (tlist, double, Nmax);
-  
+
   /* it makes no sense to mark 3-sigma outliers with <5 measurements */
   TOOFEW = MAX (5, STAR_TOOFEW);
@@ -351,8 +358,8 @@
 
       /* skip bad stars to prevent them from becoming good (on inner sample) */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      if (catalog[i].average[j].flags & STAR_BAD) continue;
 
       /* on final processing, skip stars already measured */
-      if (final && catalog[i].found[j]) continue;  
+      if (final && catalog[i].found[j]) continue;
 
       /* accumulate list of valid measurements */
@@ -360,16 +367,16 @@
       N = 0;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	/* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
-	Mcal  = getMcal  (m, i);
-	if (isnan(Mcal)) continue;
-	Mmos  = getMmos  (m, i);
-	if (isnan(Mmos)) continue;
-	Mgrid = getMgrid (m, i);
-	if (isnan(Mgrid)) continue;
-
-	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-	list[N] = Msys - Mcal - Mmos - Mgrid;
-	dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-	N++;
+        /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
+        Mcal  = getMcal  (m, i);
+        if (isnan(Mcal)) continue;
+        Mmos  = getMmos  (m, i);
+        if (isnan(Mmos)) continue;
+        Mgrid = getMgrid (m, i);
+        if (isnan(Mgrid)) continue;
+
+        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+        list[N] = Msys - Mcal - Mmos - Mgrid;
+        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+        N++;
       }
       if (N <= TOOFEW) continue;
@@ -380,8 +387,8 @@
       stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */
       for (k = m = 0; k < N; k++) {
-	if (fabs (list[k] - stats.median) < Ns*stats.sigma) {
-	  list[m] = list[k];
-	  m++;
-	}
+        if (fabs (list[k] - stats.median) < Ns*stats.sigma) {
+          list[m] = list[k];
+          m++;
+        }
       }
       initstats ("MEAN");
@@ -395,18 +402,18 @@
       N = 0;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	/* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
-	Mcal  = getMcal  (m, i);
-	if (isnan(Mcal)) continue;
-	Mmos  = getMmos  (m, i);
-	if (isnan(Mmos)) continue;
-	Mgrid = getMgrid (m, i);
-	if (isnan(Mgrid)) continue;
-
-	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-	list[N] = Msys - Mcal - Mmos - Mgrid;
-	dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-	ilist[N] = m;
-	N++;
-	Nave ++;
+        /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
+        Mcal  = getMcal  (m, i);
+        if (isnan(Mcal)) continue;
+        Mmos  = getMmos  (m, i);
+        if (isnan(Mmos)) continue;
+        Mgrid = getMgrid (m, i);
+        if (isnan(Mgrid)) continue;
+
+        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+        list[N] = Msys - Mcal - Mmos - Mgrid;
+        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+        ilist[N] = m;
+        N++;
+        Nave ++;
       }
       if (N < TOOFEW) continue;
@@ -414,8 +421,8 @@
       /* mark bad measures */
       for (k = 0; k < N; k++) {
-	if (fabs (list[k] - stats.median) > Ns*stats.sigma) {
-	  catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
-	  Ndel ++;
-	}
+        if (fabs (list[k] - stats.median) > Ns*stats.sigma) {
+          catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
+          Ndel ++;
+        }
       }
       IMAGE_BAD = image_bad;
@@ -437,5 +444,5 @@
   Ntot = 0;
   for (i = 0; i < Ncatalog; i++) {
-    Ntot += catalog[i].Naverage; 
+    Ntot += catalog[i].Naverage;
   }
 
@@ -448,18 +455,18 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      if (catalog[i].average[j].flags & STAR_BAD) continue;
       m = catalog[i].average[j].measureOffset;
 
       N = 0;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	Mcal = getMcal  (m, i);
-	if (isnan(Mcal)) continue;
-	Mmos = getMmos  (m, i);
-	if (isnan(Mmos)) continue;
-	Mgrid = getMgrid (m, i);
-	if (isnan(Mgrid)) continue;
-	N++;
-      }
-      
+        Mcal = getMcal  (m, i);
+        if (isnan(Mcal)) continue;
+        Mmos = getMmos  (m, i);
+        if (isnan(Mmos)) continue;
+        Mgrid = getMgrid (m, i);
+        if (isnan(Mgrid)) continue;
+        N++;
+      }
+
       list[n] = N;
       dlist[n] = 1;
@@ -482,5 +489,5 @@
   Ntot = 0;
   for (i = 0; i < Ncatalog; i++) {
-    Ntot += catalog[i].Naverage; 
+    Ntot += catalog[i].Naverage;
   }
 
@@ -493,5 +500,5 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      if (catalog[i].average[j].flags & STAR_BAD) continue;
 
       Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
@@ -518,5 +525,5 @@
   Ntot = 0;
   for (i = 0; i < Ncatalog; i++) {
-    Ntot += catalog[i].Naverage; 
+    Ntot += catalog[i].Naverage;
   }
 
@@ -529,5 +536,5 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      if (catalog[i].average[j].flags & STAR_BAD) continue;
 
       dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
@@ -559,5 +566,5 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
+      if (catalog[i].average[j].flags & STAR_BAD) continue;
       dMrel = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
       bin = dMrel / 0.00025;
@@ -613,5 +620,5 @@
   Graphdata graphdata;
 
-  N = 0; 
+  N = 0;
   for (i = 0; i < Ncatalog; i++) {
     N += catalog[i].Naverage;