Index: /branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- /branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/ImageOps.c	(revision 42704)
+++ /branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/ImageOps.c	(revision 42705)
@@ -531,6 +531,9 @@
 
   // until the analysis has converged a bit, do not use the IRLS analysis
-  // default is MaxIterations = 10
-  if (UseStandardOLS(ZPT_IMAGES)) {
+  // default is MaxIterations = 10.
+
+  // However, for stack and warp calibration, we are not iteratively solving for the
+  // zero points.  In this case, use IRLS in a single pass
+  if (!CALIBRATE_STACKS_AND_WARPS && UseStandardOLS(ZPT_IMAGES)) {
     brightStars.MaxIterations = 0;
     kronStars.MaxIterations = 0;
@@ -672,5 +675,5 @@
 
       if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
-	fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f : %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysKron, MrelKron, Mmos, Mgrid, Mflat, kronStars.alldata->yVector[Nkron]);
+	fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f %6.3f %6.3f | %5.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, MsysKron, MrelKron, Mmos, Mgrid, Mflat, Moff, psfStars.alldata->dyVector[Nref]);
       }
 
@@ -702,17 +705,5 @@
     if (mark) continue;
 
-    // use liststats to find the 20-pct, median, 80-pct points
-    StatType stats;
-    liststats_setmode (&stats, "MEDIAN");
-    liststats (psfStars.alldata->yVector, NULL, NULL, Nref, &stats);
-    double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
-
-    // soften the individual errors with 10% of the scatter above
-    for (j = 0; j < Nref; j++) {
-      double newSigma = hypot(psfStars.alldata->dyVector[j], 0.1*altSigma);
-      psfStars.alldata->dyVector[j] = newSigma;
-    }
-
-    // soften the errors based on the scatter
+    // soften the errors based on the core scatter
     FitDataSetSoften (&psfStars, Nref);
 
Index: /branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/bcatalog.c
===================================================================
--- /branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/bcatalog.c	(revision 42704)
+++ /branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/bcatalog.c	(revision 42705)
@@ -1,3 +1,23 @@
 # include "relphot.h"
+
+# define BAD_COUNT(FIELD) \
+  badCount.FIELD ++; \
+  if (doTest1 && (catalog[0].measure[offset].imageID == TEST_IMAGE1)) { badTest1.FIELD ++; } \
+  if (doTest2 && (catalog[0].measure[offset].imageID == TEST_IMAGE2)) { badTest2.FIELD ++; }
+
+typedef struct {
+  int Ncode;
+  int Ntime;
+  int Ndophot;
+  int Nmag;
+  int Nsigma;
+  int Nimag;
+  int Nfew;
+  int Ngalaxy;
+  int Npsfqf;
+  int Nnan;
+  int Nbad;
+  int Npoor;
+} MyBadCountType;
 
 int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog);
@@ -9,5 +29,6 @@
   off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
   float mag;
-  int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
+
+  // int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
 
   int Nsecfilt = GetPhotcodeNsecfilt ();
@@ -26,5 +47,15 @@
   int NmMin = catalog[0].Nmeasure;
 
-  Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
+  MyBadCountType badCount = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
+  MyBadCountType badTest1 = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
+  MyBadCountType badTest2 = {0,0,0, 0,0,0, 0,0,0, 0,0,0};
+
+  int doTest1 = (TEST_IMAGE1 >= 0);
+  int doTest2 = (TEST_IMAGE2 >= 0);
+
+  FILE *ftest1 = (doTest1) ? fopen ("test1.meas.dat", "a") : NULL;
+  FILE *ftest2 = (doTest2) ? fopen ("test2.meas.dat", "a") : NULL;
+
+  // Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
 
   // copy the following fields to the subcatalog:
@@ -57,4 +88,17 @@
 
       offset = catalog[0].average[i].measureOffset + j;
+
+      if (doTest1) {
+	Measure *m = &catalog[0].measure[offset];
+	if (m->imageID == TEST_IMAGE1) {
+	  fprintf (ftest1, "%5d %5d %5d | %5d %12d %5.3f %5.3f %5.3f %5.3f %4.3f %4.3f 0x%08x\n", m->detID, m->objID, m->catID, m->photcode, m->t, m->M, m->Map, m->Mkron, m->dM, m->psfQF, m->psfQFperf, m->photFlags);
+	}
+      }
+      if (doTest2) {
+	Measure *m = &catalog[0].measure[offset];
+	if (m->imageID == TEST_IMAGE2) {
+	  fprintf (ftest2, "%5d %5d %5d | %5d %12d %5.3f %5.3f %5.3f %5.3f %4.3f %4.3f 0x%08x\n", m->detID, m->objID, m->catID, m->photcode, m->t, m->M, m->Map, m->Mkron, m->dM, m->psfQF, m->psfQFperf, m->photFlags);
+	}
+      }
 
       /* select measurements by photcode */
@@ -69,5 +113,5 @@
       }
       if (!found) {
-	Ncode ++; 
+	BAD_COUNT(Ncode);
 	continue; 
       }
@@ -75,31 +119,31 @@
       /* select measurements by time */
       if (TimeSelect) {
-	if (catalog[0].measure[offset].t < TSTART) { Ntime ++; continue; }
-	if (catalog[0].measure[offset].t > TSTOP)  { Ntime ++; continue; }
+	if (catalog[0].measure[offset].t < TSTART) { BAD_COUNT(Ntime); continue; }
+	if (catalog[0].measure[offset].t > TSTOP)  { BAD_COUNT(Ntime); continue; }
       }
 
       /* select measurements by quality */
-      if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { Ndophot ++; continue; }
+      if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { BAD_COUNT(Ndophot); continue; }
 
       float Maper = USE_APER_FOR_STARGAL ? catalog[0].measure[offset].Map : catalog[0].measure[offset].Mkron;
 
       // skip garbage measurements
-      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF     < PSF_QF_MIN)) { Npsfqf ++; continue; }
-      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < PSF_QF_MIN)) { Npsfqf ++; continue; }
-      if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
-      if (isnan(Maper)) { Nnan ++; continue; }
+      if (isnan(catalog[0].measure[offset].psfQF)     || (catalog[0].measure[offset].psfQF     < PSF_QF_MIN)) { BAD_COUNT(Npsfqf); continue; }
+      if (isnan(catalog[0].measure[offset].psfQFperf) || (catalog[0].measure[offset].psfQFperf < PSF_QF_MIN)) { BAD_COUNT(Npsfqf); continue; }
+      if (isnan(catalog[0].measure[offset].M)) { BAD_COUNT(Nnan); continue; }
+      if (isnan(Maper)) { BAD_COUNT(Nnan); continue; }
 
       // require 0x01 in photFlags (fitted with a PSF) -- add to the test data
-      if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
+      if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { BAD_COUNT(Nbad); continue; }
 
       // very loose cut on PSF - APER mag (Map or Mkron) -- this lightly rejects galaxies & other oddities
       float Mkp = catalog[0].measure[offset].M - Maper;
-      if (fabs(Mkp) > 1.0) { Nbad ++; continue; }
+      if (fabs(Mkp) > 1.0) { BAD_COUNT(Nbad); continue; }
 
       if (catalog[0].measure[offset].photFlags & code->photomBadMask) { 
-	Nbad++; 
+	BAD_COUNT(Nbad); 
 	continue;
       }
-      if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { Npoor++; continue;}
+      if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { BAD_COUNT(Npoor); continue;}
 
       // weak test for galaxies
@@ -112,6 +156,6 @@
       /* select measurements by mag limit */
       mag = PhotCat (&catalog[0].measure[offset], MAG_CLASS_PSF);
-      if (mag > MAG_LIM) { Nmag ++; continue; } 
-      if (mag < MAG_MIN) { Nmag ++; continue; } 
+      if (mag > MAG_LIM) { BAD_COUNT(Nmag); continue; } 
+      if (mag < MAG_MIN) { BAD_COUNT(Nmag); continue; } 
       // note that these two cuts are applied to all photcodes. these are also applied to
       // the nominal apparent mags : if the expected sensitivity is very different between
@@ -120,11 +164,11 @@
 
       /* select measurements by measurement error */
-      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) { Nsigma ++; continue; }
+      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) { BAD_COUNT(Nsigma); continue; }
 
       /* select measurements by mag limit */
       if (ImagSelect) {
 	mag = PhotInst (&catalog[0].measure[offset], MAG_CLASS_PSF);
-	if (mag < ImagMin) { Nimag ++; continue; }
-	if (mag > ImagMax) { Nimag ++; continue; }
+	if (mag < ImagMin) { BAD_COUNT(Nimag); continue; }
+	if (mag > ImagMax) { BAD_COUNT(Nimag); continue; }
       }
       // note that instrumental mag cuts are reasonable for individual exposures using the same
@@ -155,5 +199,5 @@
     if (nEXT >= nPSF) {
       Nmeasure -= Nm;
-      Ngalaxy ++;
+      BAD_COUNT(Ngalaxy);
       continue; 
     }
@@ -165,5 +209,5 @@
     if (Nm <= STAR_TOOFEW) { /* enough measurements total? */
       Nmeasure -= Nm;
-      Nfew ++;
+      BAD_COUNT(Nfew);
       continue; 
     }
@@ -175,5 +219,5 @@
     if (!anySecfiltGood) {
       Nmeasure -= Nm;
-      Nfew ++;
+      BAD_COUNT(Nfew);
       continue; 
     }
@@ -212,5 +256,14 @@
 	     subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  catalog[0].Naverage, catalog[0].Nmeasure, catalog[0].filename, NmMin, NmMax);
     fprintf (stderr, "rejections: %d stars have too few measures:\n   %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
-	     Nfew, Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
+	     badCount.Nfew, badCount.Ncode, badCount.Ntime, badCount.Ndophot, badCount.Nmag, badCount.Nsigma, badCount.Nimag, badCount.Npsfqf, badCount.Nnan, badCount.Ngalaxy, badCount.Nbad, badCount.Npoor);
+  }
+
+  if (doTest1) {
+    fprintf (stderr, "TEST1 rejections: %d stars have too few measures: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
+	     badTest1.Nfew, badTest1.Ncode, badTest1.Ntime, badTest1.Ndophot, badTest1.Nmag, badTest1.Nsigma, badTest1.Nimag, badTest1.Npsfqf, badTest1.Nnan, badTest1.Ngalaxy, badTest1.Nbad, badTest1.Npoor);
+  }
+  if (doTest2) {
+    fprintf (stderr, "TEST2 rejections: %d stars have too few measures: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
+	     badTest2.Nfew, badTest2.Ncode, badTest2.Ntime, badTest2.Ndophot, badTest2.Nmag, badTest2.Nsigma, badTest2.Nimag, badTest2.Npsfqf, badTest2.Nnan, badTest2.Ngalaxy, badTest2.Nbad, badTest2.Npoor);
   }
 
@@ -220,5 +273,8 @@
   }
 
-  free (Nvalid)
+  free (Nvalid);
+
+  if (ftest1) { fclose (ftest1); }
+  if (ftest2) { fclose (ftest2); }
 
   return (TRUE);
Index: /branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/fit1d_irls.c
===================================================================
--- /branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/fit1d_irls.c	(revision 42704)
+++ /branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/fit1d_irls.c	(revision 42705)
@@ -374,5 +374,6 @@
   StatType stats;
   liststats_setmode (&stats, "MEDIAN");
-  liststats (dataset->alldata->yVector, NULL, NULL, Nvalues, &stats);
+  liststats (dataset->alldata->yVector, dataset->alldata->dyVector, NULL, Nvalues, &stats);
+  // liststats sorts yVector and dyVector together
 
   double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
@@ -380,5 +381,5 @@
   // soften the individual errors with 10% of the scatter above
   for (int j = 0; j < Nvalues; j++) {
-    double newSigma = hypot(dataset->alldata->dyVector[j], 0.1*altSigma);
+    double newSigma = hypot(dataset->alldata->dyVector[j], 0.5*altSigma);
     dataset->alldata->dyVector[j] = newSigma;
   }