Index: trunk/Ohana/src/relphot/Makefile
===================================================================
--- trunk/Ohana/src/relphot/Makefile	(revision 41451)
+++ trunk/Ohana/src/relphot/Makefile	(revision 41453)
@@ -21,4 +21,6 @@
 
 test_liststats: $(BIN)/test_liststats.$(ARCH)
+
+test_fitirls: $(BIN)/test_fitirls.$(ARCH)
 
 install: $(DESTBIN)/relphot $(DESTBIN)/relphot_client
@@ -48,4 +50,5 @@
 $(SRC)/global_stats.$(ARCH).o	 \
 $(SRC)/initialize.$(ARCH).o	 \
+$(SRC)/fit1d_irls.$(ARCH).o	 \
 $(SRC)/liststats.$(ARCH).o	 \
 $(SRC)/load_catalogs.$(ARCH).o	 \
@@ -107,4 +110,5 @@
 $(SRC)/synthetic_mags.$(ARCH).o	 \
 $(SRC)/plotstuff.$(ARCH).o	 \
+$(SRC)/fit1d_irls.$(ARCH).o	 \
 $(SRC)/liststats.$(ARCH).o	 \
 $(SRC)/initialize.$(ARCH).o	 \
@@ -137,2 +141,9 @@
 $(BIN)/test_liststats.$(ARCH): $(TEST_LISTSTATS)
 
+TEST_FITIRLS = \
+$(SRC)/fit1d_irls.$(ARCH).o	 \
+$(SRC)/test_fitirls.$(ARCH).o
+
+$(TEST_FITIRLS): $(INC)/relphot.h
+$(BIN)/test_fitirls.$(ARCH): $(TEST_FITIRLS)
+
Index: trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- trunk/Ohana/src/relphot/include/relphot.h	(revision 41451)
+++ trunk/Ohana/src/relphot/include/relphot.h	(revision 41453)
@@ -91,6 +91,4 @@
 
 typedef struct {
-  unsigned int start;
-  unsigned int stop;
   short photcode;
   float McalPSF;
@@ -110,9 +108,16 @@
   off_t NMEASURE;
   off_t *measure;
-
-  off_t Ncatalog; // XX not needed
-  off_t NCATALOG; // XX not needed
   off_t *catalog;
+
+  void *parent;
 } TGroup; 
+
+// we have an array of TGroup times, each pointing to N sets of data values
+typedef struct {
+  unsigned int start;
+  unsigned int stop;
+  TGroup *byCode; // each of these contains the collection of images for the time and photcode
+  int      nCode;
+} TGTimes;
 
 typedef enum {
@@ -163,4 +168,53 @@
 
 } StatDataSet;
+
+typedef struct {
+  double *xVector;	      // complete list of values in independent variable (optional)
+  double *yVector;	      // complete list of values available
+  double *dyVector;	      // complete list of errors available
+} FitDataType;
+  
+// this structure carries the data and pre-allocated arrays for
+// 1D fitting with arbitrary order for data with external weights,
+// optional priors, irls iterations, and bootstrap resampling
+typedef struct {
+  int     order;	      // order of fit (e.g., y = C0 + C1*x has order = 1)
+  int     nterm;	      // number of fitted parameters (order + 1)
+  int     mterm;	      // number of summations needed for OLS (2*order + 1)
+
+  int     Nlist;	      // total number of measurements in list (may be more than the number to use)
+  int     Nbootstrap;	      // number of bootstrap iterations
+
+  // input data & parameters: 
+  FitDataType *alldata;	      // full data set
+
+  double *bPriorValue; 	      // prior value
+  double *bPriorSigma;	      // prior sigma
+
+  // internal / temporary arrays:
+  FitDataType *keepdata;      // unclipped data for bootstrap analysis
+  FitDataType *sample;	      // sample for a bootstrap iteration
+
+  double *yOffVector;	      // difference between yFit and yObs
+
+  double *tmpVector;	      // internal copy of values to use in, e.g., calculation of median
+  double *wtIRLS; 	      // IRLS weights (distance from model value weighted by standard error)
+
+  double **cArray;
+  double  *sumVector;	      // pre-allocated to mterm (2*order + 1)
+
+  double **bArray;
+  double **bSaveArray;
+  double **bBootArray;	      // values generated by the bootstrap resampling 
+
+  // output results:
+  double *bSigma;	      // 1 sigma range of parameters
+
+  double  min;		      // min of bSigma[0] values
+  double  max;		      // max of bSigma[0] values
+  double  sigma;	      // (sample) standard deviation of bSigma[0] values
+  double  chisq;	      // chisq of fit (unmasked values only)
+  int     Nmeas;	      // number of unmasked values used in fit
+} FitDataSet;
 
 typedef struct {
@@ -311,4 +365,6 @@
 double IMAGE_SCATTER;
 double IMAGE_OFFSET;
+double NIGHT_SCATTER;
+double NIGHT_OFFSET;
 double STAR_SCATTER;
 double STAR_CHISQ;
@@ -497,4 +553,12 @@
 int           liststats_irls      PROTO((StatDataSet *dataset, int Npoints, StatType *stats));
 int           liststats_fit1d     PROTO((double *value, double *err, double *x, int Npts, StatType *stats, double *dk));
+double        weight_cauchy       PROTO((double x));
+double        VectorFractionInterpolate PROTO((double *values, float fraction, int Npoints));
+
+// fit1d_irls:
+int           fit1d_irls          PROTO((FitDataSet *dataset, int Npoints));
+void          FitDataSetFree      PROTO((FitDataSet *dataset));
+void          FitDataSetAlloc     PROTO((FitDataSet *dataset, int Nmax, int order, int Nbootstrap));
+void          FitDataSetAddPriors PROTO((FitDataSet *dataset));
 
 unsigned int *ReadTGroupFile      PROTO((FILE *f, int *nelem));
@@ -502,5 +566,5 @@
 void          initTGroupsMcal     PROTO((void));
 TGroup       *getTGroupForImage   PROTO((off_t im));
-off_t         findTGroup          PROTO((unsigned int start, int photcode));
+TGroup       *findTGroup          PROTO((unsigned int start, int photcode));
 
 Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile));
@@ -527,4 +591,5 @@
 void          plot_star_coords    PROTO((Catalog *catalog, int Ncatalog));
 void          plot_stars          PROTO((Catalog *catalog, int Ncatalog));
+void          plot_setMcal        PROTO((double *list, int Npts));
 
 void          plot_list_add       PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N));
Index: trunk/Ohana/src/relphot/src/ConfigInit.c
===================================================================
--- trunk/Ohana/src/relphot/src/ConfigInit.c	(revision 41451)
+++ trunk/Ohana/src/relphot/src/ConfigInit.c	(revision 41453)
@@ -1,3 +1,9 @@
 # include "relphot.h"
+
+// do not use with %s
+# define DefConfig(NAME, FMT, DEF, VAR) {				\
+  char *status = ScanConfig (config, NAME, FMT, 0, &VAR);		\
+  if (status == NULL) { VAR = DEF; }					\
+  }
 
 void ConfigInit (int *argc, char **argv) {
@@ -20,7 +26,7 @@
   GetConfig (config, "SIGMA_LIM",              "%lf", 0, &SIGMA_LIM);
 
-  if (!ScanConfig (config, "RELPHOT_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01;
 
   GetConfig (config, "STAR_SCATTER",           "%lf", 0, &STAR_SCATTER);
+
   GetConfig (config, "IMAGE_SCATTER",          "%lf", 0, &IMAGE_SCATTER);
   GetConfig (config, "IMAGE_OFFSET",           "%lf", 0, &IMAGE_OFFSET);
@@ -47,7 +53,9 @@
   snprintf (ImageCat, DVO_MAX_PATH, "%s/Images.dat", CATDIR);
 
-  if (!ScanConfig (config, "SKY_DEPTH",         "%d",  0, &SKY_DEPTH)) {
-    SKY_DEPTH = 2;
-  }
+  DefConfig ("RELPHOT_IMFIT_SYS_SIGMA_LIM", "%lf", 0.01, IMFIT_SYS_SIGMA_LIM);
+  DefConfig ("NIGHT_SCATTER",               "%lf", 0.05, NIGHT_SCATTER);
+  DefConfig ("NIGHT_OFFSET",                "%lf", 0.50, NIGHT_OFFSET);
+  DefConfig ("SKY_DEPTH",                   "%d",     2, SKY_DEPTH);
+
   if (!ScanConfig (config, "SKY_TABLE",         "%s",  0, SKY_TABLE)) {
     SKY_TABLE[0] = 0;
@@ -90,2 +98,3 @@
   return;
 }
+
Index: trunk/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/ImageOps.c	(revision 41451)
+++ trunk/Ohana/src/relphot/src/ImageOps.c	(revision 41453)
@@ -1,6 +1,4 @@
 # include "relphot.h"
 # define BASIC_STATS 0
-
-void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
 
 // the MeasureToImage, ImageToCatalog, and ImageToMeasure arrays are a substantial part of the memory footprint.
@@ -588,7 +586,4 @@
   int mark, bad, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys;
 
-  StatType stats;
-  liststats_setmode (&stats, "INNER_WTMEAN");
-
   // FREEZE_IMAGES only applies to mosaic data (eg, gpc1)
 
@@ -608,7 +603,9 @@
   }
 
-  StatDataSet *kronStars   = StatDataSetAlloc (1, Nmax);
-  StatDataSet *psfStars    = StatDataSetAlloc (1, Nmax);
-  StatDataSet *brightStars = StatDataSetAlloc (1, Nmax);
+  // we are making a 0-order fit and not doing bootstrap analysis:
+  FitDataSet psfStars, kronStars, brightStars;
+  FitDataSetAlloc (&psfStars,    Nmax, 0, 0);
+  FitDataSetAlloc (&kronStars,   Nmax, 0, 0);
+  FitDataSetAlloc (&brightStars, Nmax, 0, 0);
 
   Nfew = Nbad = Nmos = Ngrp = Ngrid = Nrel = Nsys = 0;
@@ -713,11 +710,6 @@
 
     skip:
-      psfStars->flxlist[Nref] = MsysPSF - MrelPSF - Mmos - Mgrp - Mgrid + Mflat;
-      psfStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-      psfStars->wgtlist[Nref] = 1;
-      psfStars->msklist[Nref] = 0;
-      if (fabs(psfStars->flxlist[Nref]) > 0.03) {
-	// fprintf (stderr, "deviant\n");
-      }
+      psfStars.alldata-> yVector[Nref] = MsysPSF - MrelPSF - Mmos - Mgrp - Mgrid + Mflat;
+      psfStars.alldata->dyVector[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
 
       float MrelKron = getMrel  (catalog, m, c, MAG_CLASS_KRON, MAG_SRC_CHP);
@@ -725,20 +717,16 @@
 
       if (isfinite(MrelKron) && isfinite(MsysKron)) {
-	kronStars->flxlist[Nkron] = MsysKron - MrelKron - Mmos - Mgrp - Mgrid + Mflat;
-	kronStars->errlist[Nkron] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-	kronStars->wgtlist[Nkron] = 1;
-	kronStars->msklist[Nkron] = 0;
+	kronStars.alldata-> yVector[Nkron] = MsysKron - MrelKron - Mmos - Mgrp - Mgrid + Mflat;
+	kronStars.alldata->dyVector[Nkron] = psfStars.alldata->dyVector[Nref];
 	Nkron ++;
       }
 
       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->flxlist[Nref]);
+	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]);
       }
 
       if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
-	brightStars->flxlist[Nbright] = psfStars->flxlist[Nref];
-	brightStars->errlist[Nbright] = psfStars->errlist[Nref];
-	brightStars->wgtlist[Nbright] = 1;
-	brightStars->msklist[Nbright] = 0;
+	brightStars.alldata-> yVector[Nbright] = psfStars.alldata-> yVector[Nref];
+	brightStars.alldata->dyVector[Nbright] = psfStars.alldata->dyVector[Nref];
 	Nbright ++;
       }
@@ -762,28 +750,15 @@
 
     // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
-# if (BASIC_STATS) 
-    liststats (psfStars->flxlist, psfStars->errlist, NULL, Nref, &stats);
-# else
-    liststats_irls (psfStars, Nref, &stats);
-# endif
-    image[i].McalPSF    = stats.mean;
-    image[i].dMcal      = stats.error;
-    image[i].nFitPhotom = Nref;
-    image[i].McalChiSq  = stats.chisq;
+    fit1d_irls (&psfStars, Nref);
+    image[i].McalPSF    = psfStars.bSaveArray[0][0];
+    image[i].dMcal      = psfStars.bSigma[0];
+    image[i].nFitPhotom = psfStars.Nmeas;
+    image[i].McalChiSq  = psfStars.chisq;
     Ncalibrated ++;
 
     // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
-# if (BASIC_STATS) 
-    liststats (kronStars->flxlist, kronStars->errlist, NULL, Nkron, &stats);
-# else
-    liststats_irls (kronStars, Nref, &stats);
-# endif
-    image[i].McalAPER   = stats.mean;
-
-    if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
-      for (j = 0; j < Nref; j++) {
-	fprintf (stderr, "%1d, %8d : %6.3f  %6.3f  %6.3f  %d\n", (int) i, (int) image[i].imageID, kronStars->flxlist[j], kronStars->errlist[j], kronStars->wgtlist[j], kronStars->msklist[j]);
-      }
-    }
+    // XXX: apply airmass from above and fit only zp?
+    fit1d_irls (&kronStars, Nkron);
+    image[i].McalAPER   = kronStars.bSaveArray[0][0];
 
     if ((image[i].imageID == TEST_IMAGE1) || (image[i].imageID == TEST_IMAGE2)) {
@@ -791,16 +766,12 @@
     }
 
-    if (!mark && VERBOSE_IMAGE) {
-      fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].McalAPER, image[i].dMcal);
-    }
-
     if (PLOTSTUFF) {
       fprintf (stderr, "Mcal for : %s : %7.4f %7.4f\n", image[i].name, image[i].McalAPER, image[i].dMcal);
-      plot_setMcal (psfStars->flxlist, Nref, &stats, CLOUD_TOLERANCE);
+      plot_setMcal (psfStars.alldata-> yVector, Nref);
     }
 
     // bright end scatter
-    liststats (brightStars->flxlist, brightStars->errlist, NULL, Nbright, &stats);
-    image[i].dMagSys = stats.sigma;
+    // fit1d_irls (&brightStars, Nbright);
+    // image[i].dMagSys = brightStars.sigma;
 
     if (image[i].McalPSF < -CLOUD_TOLERANCE) {
@@ -812,8 +783,4 @@
     image[i].ubercalDist = minUbercalDist + 1;
   }
-
-  StatDataSetFree (brightStars, 1);
-  StatDataSetFree (kronStars, 1);
-  StatDataSetFree (psfStars, 1);
 
   fprintf (stderr, "%d images calibrated\n", Ncalibrated);
@@ -825,4 +792,9 @@
     MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
   }
+
+  FitDataSetFree (&brightStars);
+  FitDataSetFree (&kronStars);
+  FitDataSetFree (&psfStars);
+
   return;
 }
@@ -891,7 +863,5 @@
 
 static int setMcal_init_done = FALSE;
-void plot_setMcal (double *list, int Npts, StatType *stats, float clouds) {
-  OHANA_UNUSED_PARAM(stats);
-  OHANA_UNUSED_PARAM(clouds);
+void plot_setMcal (double *list, int Npts) {
 
   off_t i;
Index: trunk/Ohana/src/relphot/src/MosaicOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 41451)
+++ trunk/Ohana/src/relphot/src/MosaicOps.c	(revision 41453)
@@ -1,4 +1,3 @@
 # include "relphot.h"
-void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
 off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start);
 
@@ -1281,5 +1280,5 @@
   if (PLOTSTUFF) {
     fprintf (stderr, "Mmos: %6.3f %6.3f +/- %6.3f %5d %5d | %s\n", stats.mean, stats.median, stats.sigma, stats.Nmeas, N, image[MosaicToImage[Nmos][0]].name);
-    plot_setMcal (psfMagList, N, &stats, CLOUD_TOLERANCE);
+    plot_setMcal (psfMagList, N);
   }
 
Index: trunk/Ohana/src/relphot/src/StarOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/StarOps.c	(revision 41451)
+++ trunk/Ohana/src/relphot/src/StarOps.c	(revision 41453)
@@ -228,7 +228,9 @@
 
 void SetMrelInfoFree (SetMrelInfo *results) {
+
   StatDataSetFree (results->psfData,  results->Nsecfilt);
   StatDataSetFree (results->aperData, results->Nsecfilt);
   StatDataSetFree (results->kronData, results->Nsecfilt);
+
   FREE (results->psfqf_list);
   FREE (results->psfqfperf_list);
Index: trunk/Ohana/src/relphot/src/TGroupOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/TGroupOps.c	(revision 41451)
+++ trunk/Ohana/src/relphot/src/TGroupOps.c	(revision 41453)
@@ -1,40 +1,34 @@
 # include "relphot.h"
 
-// XXX is this needed?
-void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
-
-// see discussion in ImagesOps.c re: IDX_T
-
-// array of tgroup definition structures
-static off_t  Ntgroup;
-static TGroup *tgroup;
-
-// tgroup index for given image : ImageToTGroup[ImageIndex] = TGroupIndex (ImageIndex : 0 < Nimage)
-static off_t   *ImageToTGroup;
+// tgroupTimes carries the times of the tgroups (initially, just the photometric nights)
+static     int  NtgroupTimes = 0;
+static TGTimes **tgroupTimes = NULL;
+
+// TGroup pointers for each image
+static TGroup **ImageToTGroup = NULL;
 
 // elsewhere, we have loaded a set of catalogs with measures (catalog[cat].measure[meas])
 
-// each tgroup has N_onTGroup[TGroupIndex] measurements
-// static off_t    *N_onTGroup;   // actual number of measurements on tgroup	 
-// static off_t    *N_ONTGROUP;   // allocated number of measurements on tgroup   
-
-// relationships between the measure,catalog set and the tgroups:
-static off_t   **MeasureToTGroup; // TGroup index from measure,catalog   : MeasureToTGroup[cat][meas] = TGroupIndex 
-// static off_t   **TGroupToCatalog; // catalog for given measure on tgroup : TGroupCatalog[TGroupIndex][i] = cat (i : 0 < NonTGroup[TGroupIndex])
-// static off_t   **TGroupToMeasure; // measure for given measure on tgroup : TGroupMeasure[TGroupIndex][i] = cat (i : 0 < NonTGroup[TGroupIndex])
-
-// tgroupTimes carries the times of the photometric nights
-static          int NtgroupTimes;
-static unsigned int *tgroupTimes;
+// TGroup pointers for each measure,catalog set
+static TGroup  ***MeasureToTGroup = NULL;
 
 void sort_times (unsigned int *T, int N);
 
-/* find tgroups (unique time periods) (NOTE : we do NOT require matching photcodes...)
-   we need to have a tgroup for each photcode available for a given photometric night
- */
+void sort_tgtimes (TGTimes **T, int N) {
+
+# define SWAPFUNC(A,B){ TGTimes *tmp;		\
+    tmp = T[A]; T[A] = T[B]; T[B] = tmp;	\
+  }
+# define COMPARE(A,B)(T[A][0].start < T[B][0].start)
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
 
 // load the photometric nights from a file : these will be used to define the tgroups 
-// XXX call this in args.c where TGROUP_ZEROPT is set
-// XXX cannot call this in args.c : the Nphotcodes have not yet been loaded.
+// this is called in args.c after TGROUP_ZEROPT is set
 void loadTGroups (char *filename) {
 
@@ -47,19 +41,20 @@
     exit (4);
   }
-  tgroupTimes = ReadTGroupFile (f, &NtgroupTimes);
-  if (!tgroupTimes) {
+  // ReadTGroupFile returns a list of the time values (int)
+  unsigned int *tgroupTimesRaw = ReadTGroupFile (f, &NtgroupTimes);
+  if (!tgroupTimesRaw) {
     fprintf (stderr, "error reading tgroup file %s\n", filename);
     exit (4);
   }
+  fprintf (stderr, "loaded %d tgroups from file %s\n", NtgroupTimes, filename);
 
   // sort the times for quick bisection
-  sort_times (tgroupTimes, NtgroupTimes);
-
-  // we define a tgroup for each photcode being considered
-  // below, if there are no images in this tgroup with a given
-  // photcode, it is left empty (no images)
-  Ntgroup = NtgroupTimes*Nphotcodes;
-
-  ALLOCATE (tgroup, TGroup, Ntgroup);
+  sort_times (tgroupTimesRaw, NtgroupTimes);
+
+  // we first define the collection of times (TGTimes) then for each time, we define a tgroup for
+  // each photcode being considered below, if there are no images in this tgroup with a
+  // given photcode, it is left empty (no images)
+
+  ALLOCATE (tgroupTimes, TGTimes *, NtgroupTimes);
 
   // the start of an MJD day is 14:00 HST, so the end time can be the start time + 23h 59m 59s (86399s)
@@ -67,33 +62,91 @@
 
   for (int i = 0; i < NtgroupTimes; i++) {
+    ALLOCATE (tgroupTimes[i], TGTimes, 1);
+    tgroupTimes[i][0].start = tgroupTimesRaw[i];
+    tgroupTimes[i][0].stop  = tgroupTimesRaw[i] + 86399;
+
+    ALLOCATE (tgroupTimes[i][0].byCode, TGroup, Nphotcodes);
+    tgroupTimes[i][0].nCode = Nphotcodes;
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+
     for (int j = 0; j < Nphotcodes; j++) {
-      /* a new tgroup, define ranges */
-      int Nt = i*Nphotcodes + j;
-      tgroup[Nt].start     = tgroupTimes[i];
-      tgroup[Nt].stop      = tgroupTimes[i] + 86399;
-      tgroup[Nt].McalPSF   = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
-      tgroup[Nt].McalAPER  = 0.0; // note : tgroup stores only offsets relative to the original image values
-      tgroup[Nt].dKlam     = 0.0; 
-      tgroup[Nt].dMcal     = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
-      tgroup[Nt].dMsys     = 0.0;
-      tgroup[Nt].McalChiSq = 0.0; // NAN or 0.0?
-      tgroup[Nt].flags     = 0;
-      tgroup[Nt].photcode  = photcodes[j][0].code;
+      tgroup[j].parent    = (void *) tgroupTimes[i];
+      tgroup[j].McalPSF   = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
+      tgroup[j].McalAPER  = 0.0; // note : tgroup stores only offsets relative to the original image values
+      tgroup[j].dKlam     = 0.0; 
+      tgroup[j].dMcal     = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
+      tgroup[j].dMsys     = 0.0;
+      tgroup[j].McalChiSq = 0.0; // NAN or 0.0?
+      tgroup[j].flags     = ID_IMAGE_PHOTOM_UBERCAL; // we start by treating the suspected photometric nights as photometric
+      tgroup[j].photcode  = photcodes[j][0].code;
     
-      tgroup[Nt].NIMAGE    = 1000;
-      tgroup[Nt].Nimage    = 0;
-      ALLOCATE (tgroup[Nt].image, off_t, tgroup[Nt].NIMAGE);
-      tgroup[Nt].image[0]  = -1;
-
-      tgroup[Nt].Nmeasure  = 0;
-      tgroup[Nt].NMEASURE  = 1000;
-      ALLOCATE (tgroup[Nt].measure, off_t, tgroup[Nt].NMEASURE);
-      ALLOCATE (tgroup[Nt].catalog, off_t, tgroup[Nt].NMEASURE);
-    }
-  }
+      tgroup[j].NIMAGE    = 1000;
+      tgroup[j].Nimage    = 0;
+      ALLOCATE (tgroup[j].image, off_t, tgroup[j].NIMAGE);
+      tgroup[j].image[0]  = -1;
+
+      tgroup[j].Nmeasure  = 0;
+      tgroup[j].NMEASURE  = 1000;
+      ALLOCATE (tgroup[j].measure, off_t, tgroup[j].NMEASURE);
+      ALLOCATE (tgroup[j].catalog, off_t, tgroup[j].NMEASURE);
+    }
+  }
+
+  free (tgroupTimesRaw);
   return;
 }
 
-// assign the images to the tgroups
+// create a new TGroup based on the supplied time and photcode
+void extendTGroups (unsigned int tzero) {
+
+  if (!TGROUP_ZEROPT) return;
+  myAssert (tgroupTimes, "oops, tgroupTimes not allocated");
+
+  // index of the new entry
+  int Ntgt = NtgroupTimes;
+
+  NtgroupTimes ++;
+  REALLOCATE (tgroupTimes, TGTimes *, NtgroupTimes);
+  
+  double mjd = ohana_sec_to_mjd (tzero);
+  double mjdInt = (int) mjd;
+  time_t start = ohana_mjd_to_sec (mjdInt);
+
+  ALLOCATE (tgroupTimes[Ntgt], TGTimes, 1);
+  tgroupTimes[Ntgt][0].start = start;
+  tgroupTimes[Ntgt][0].stop  = start + 86399;
+
+  ALLOCATE (tgroupTimes[Ntgt][0].byCode, TGroup, Nphotcodes);
+  tgroupTimes[Ntgt][0].nCode = Nphotcodes;
+  TGroup *tgroup = tgroupTimes[Ntgt][0].byCode;
+
+  for (int j = 0; j < Nphotcodes; j++) {
+    tgroup[j].parent    = (void *) tgroupTimes[Ntgt];
+    tgroup[j].McalPSF   = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
+    tgroup[j].McalAPER  = 0.0; // note : tgroup stores only offsets relative to the original image values
+    tgroup[j].dKlam     = 0.0; 
+    tgroup[j].dMcal     = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
+    tgroup[j].dMsys     = 0.0;
+    tgroup[j].McalChiSq = 0.0; // NAN or 0.0?
+    tgroup[j].flags     = 0;   // tgroups NOT in the original file are not suspected as photometric
+    tgroup[j].photcode  = photcodes[j][0].code;
+    tgroup[j].NIMAGE    = 1000;
+    tgroup[j].Nimage    = 0;   // the new TGroup is not suspected to be photometric
+    ALLOCATE (tgroup[j].image, off_t, tgroup[j].NIMAGE);
+    tgroup[j].image[0]  = -1;
+    tgroup[j].Nmeasure  = 0;
+    tgroup[j].NMEASURE  = 1000;
+    ALLOCATE (tgroup[j].measure, off_t, tgroup[j].NMEASURE);
+    ALLOCATE (tgroup[j].catalog, off_t, tgroup[j].NMEASURE);
+  }
+
+  // sort the times for quick bisection
+  sort_tgtimes (tgroupTimes, NtgroupTimes);
+  return;
+}
+
+// Assign the images to the tgroups.  The initial tgroups table contains only suspected
+// photometric nights.  As we find images which are not in the current tgroups list, add
+// new non-photometric entries.
 void initTGroups (Image *subset, off_t Nsubset) {
 
@@ -105,10 +158,10 @@
    */
 
-  ALLOCATE (ImageToTGroup, off_t, Nsubset); // tgroup to which image belongs
+  ALLOCATE (ImageToTGroup, TGroup *, Nsubset); // tgroup to which image belongs
 
   // assign each image to a tgroup
   int Nsimple = 0;
   for (off_t i = 0; i < Nsubset; i++) {
-    ImageToTGroup[i] = -1;
+    ImageToTGroup[i] = NULL;
 
     // ignore non-GPC1 images
@@ -123,13 +176,14 @@
     }
 
-    off_t j = findTGroup(subset[i].tzero, subset[i].photcode);
-    if (j == -1) {
-      // images which do not belong to a tgroup should be ignored in this analysis
-      // XXX mark the images and the measurements somehow
-      continue;
-    }
-
-    // add reference from image to tgroup
-    ImageToTGroup[i] = j;
+    TGroup *myGroup = findTGroup(subset[i].tzero, subset[i].photcode);
+    if (!myGroup) {
+      // generate a new TGroup which is not photometric
+      extendTGroups (subset[i].tzero);
+      myGroup = findTGroup(subset[i].tzero, subset[i].photcode);
+      myAssert (myGroup, "oops, we just extended to include this");
+    }
+
+    // add pointer to tgroup for this image
+    ImageToTGroup[i] = myGroup;
 
     // note the array 'subset' is registered in ImageOps.c:images by initImages
@@ -137,9 +191,9 @@
 
     // add image to tgroup image list
-    tgroup[j].image[tgroup[j].Nimage] = i; // reference to entry in 'subset'
-    tgroup[j].Nimage ++;
-    if (tgroup[j].Nimage == tgroup[j].NIMAGE) {
-      tgroup[j].NIMAGE += 1000;
-      REALLOCATE (tgroup[j].image, off_t, tgroup[j].NIMAGE);
+    myGroup[0].image[myGroup[0].Nimage] = i; // reference to entry in 'subset'
+    myGroup[0].Nimage ++;
+    if (myGroup[0].Nimage == myGroup[0].NIMAGE) {
+      myGroup[0].NIMAGE += 1000;
+      REALLOCATE (myGroup[0].image, off_t, myGroup[0].NIMAGE);
     }
   }
@@ -147,5 +201,5 @@
   initTGroupsMcal ();
 
-  fprintf (stderr, "matched %d images to %d tgroups, %d simple chips not matched to tgroups\n", (int) (Nsubset - Nsimple), (int) Ntgroup, (int) Nsimple);
+  fprintf (stderr, "matched %d images to %d tgroups, %d simple chips not matched to tgroups\n", (int) (Nsubset - Nsimple), (int) NtgroupTimes, (int) Nsimple);
   return;
 }
@@ -156,12 +210,18 @@
 
   free (ImageToTGroup);
+
+  for (int i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      free (tgroup[j].image);
+      free (tgroup[j].measure);
+      free (tgroup[j].catalog);
+    }
+    free (tgroup);
+    free (tgroupTimes[i]);
+  }
   free (tgroupTimes);
 
-  for (int i = 0; i < Ntgroup; i++) {
-    free (tgroup[i].image);
-    free (tgroup[i].measure);
-    free (tgroup[i].catalog);
-  }
-  free (tgroup);
+  // free (ImageToGroup); why not?
   return;
 }
@@ -175,38 +235,43 @@
   if (!TGROUP_ZEROPT) return;
 
+  // we use liststats to calculate the median Mcal for the input images
+  // in each tgroup as a starting point.
   StatType stats;
   liststats_setmode (&stats, "MEDIAN");
 
   // init the tgroup array values
-  for (int i = 0; i < Ntgroup; i++) {
-    tgroup[i].McalPSF   = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
-    tgroup[i].McalAPER  = 0.0; // note : tgroup stores only offsets relative to the original image values
-    tgroup[i].dMcal     = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
-    tgroup[i].dMsys     = 0.0;
-    tgroup[i].McalChiSq = 0.0;// NAN or 0.0?
-
-    if (tgroup[i].Nimage == 0) continue; // no images, ignore the tgroup
-
-    // calculate the median of the Mcal values for the images in this tgroup:
-
-    ALLOCATE_PTR (McalPSF,  double, tgroup[i].Nimage);
-    ALLOCATE_PTR (McalAPER, double, tgroup[i].Nimage);
-    for (int j = 0; j < tgroup[i].Nimage; j++) {
-      int Nsub = tgroup[i].image[j];
-      McalPSF[j] = images[Nsub].McalPSF;
-      McalAPER[j] = images[Nsub].McalAPER;
-
-      // set these back to 0.0 to have no future effect
-      images[Nsub].McalPSF = 0.0;
-      images[Nsub].McalAPER = 0.0;
-    }
+  for (int i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      tgroup[j].McalPSF   = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
+      tgroup[j].McalAPER  = 0.0; // note : tgroup stores only offsets relative to the original image values
+      tgroup[j].dMcal     = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
+      tgroup[j].dMsys     = 0.0;
+      tgroup[j].McalChiSq = 0.0;// NAN or 0.0?
+
+      if (tgroup[j].Nimage == 0) continue; // no images, ignore the tgroup
+
+      // calculate the median of the Mcal values for the images in this tgroup:
+
+      ALLOCATE_PTR (McalPSF,  double, tgroup[j].Nimage);
+      ALLOCATE_PTR (McalAPER, double, tgroup[j].Nimage);
+      for (int im = 0; im < tgroup[j].Nimage; im++) {
+	int Nsub = tgroup[j].image[im];
+	McalPSF[im] = images[Nsub].McalPSF;
+	McalAPER[im] = images[Nsub].McalAPER;
+
+	// set these back to 0.0 to have no future effect
+	images[Nsub].McalPSF = 0.0;
+	images[Nsub].McalAPER = 0.0;
+      }
     
-    liststats (McalPSF, NULL, NULL, tgroup[i].Nimage, &stats);
-    tgroup[i].McalPSF = stats.median;
-    FREE (McalPSF);
+      liststats (McalPSF, NULL, NULL, tgroup[j].Nimage, &stats);
+      tgroup[j].McalPSF = stats.median;
+      FREE (McalPSF);
   
-    liststats (McalAPER, NULL, NULL, tgroup[i].Nimage, &stats);
-    tgroup[i].McalAPER = stats.median;
-    FREE (McalAPER);
+      liststats (McalAPER, NULL, NULL, tgroup[j].Nimage, &stats);
+      tgroup[j].McalAPER = stats.median;
+      FREE (McalAPER);
+    }
   }
   return;
@@ -218,9 +283,7 @@
   if (!ImageToTGroup) return NULL;
 
-  off_t m = ImageToTGroup[im];
-  if (m < 0) return NULL;
-  if (m >= Ntgroup) return NULL;
-
-  return (&tgroup[m]);
+  // test if im > Nimages / Nsubset?
+  TGroup *myGroup = ImageToTGroup[im];
+  return (myGroup);
 }
 
@@ -228,14 +291,15 @@
 // tgroupTimes is a sorted, unique list of times
 // assume a 1day range for now
-off_t findTGroup (unsigned int start, int photcode) {
-
-  // first find the matching photcode
+TGroup *findTGroup (unsigned int start, int photcode) {
+
+  // first find the matching photcode (our list of active secondary photcodes
+  // does not necessarily match the list of all secondary photcodes)
   int Ns = -1;
+  int ecode = GetPhotcodeEquivCodebyCode (photcode);
   for (int i = 0; (Ns < 0) && (i < Nphotcodes); i++) {
-    int ecode = GetPhotcodeEquivCodebyCode (photcode);
     if (photcodes[i][0].code != ecode) continue;
     Ns = i;
   }
-  if (Ns < 0) return -1;
+  if (Ns < 0) return NULL;
 
   off_t Nlo, Nhi, N;
@@ -243,33 +307,35 @@
   // find the last tgroup before start
   Nlo = 0; // first valid tgroupTimes value
-  Nhi = Ntgroup - 1; // last valid tgroupTimes value
-
-  // if start is not in this range, return -1
-  if (start < tgroupTimes[Nlo]) return (-1);
-  if (start > tgroupTimes[Nhi] + 86399) return (-1); // 1-day range for tgroup for now
+  Nhi = NtgroupTimes - 1; // last valid tgroupTimes value
+
+  // if start is not in this range, return NULL
+  if (start < tgroupTimes[Nlo][0].start) return NULL;
+  if (start > tgroupTimes[Nhi][0].start + 86399) return NULL; // 1-day range for tgroup for now
 
   while (Nhi - Nlo > 4) {
     N = 0.5*(Nlo + Nhi);
-    if (tgroupTimes[N] < start) {
+    if (tgroupTimes[N][0].start < start) {
       Nlo = MAX(N, 0);
     } else {
-      Nhi = MIN(N, Ntgroup - 1);
+      Nhi = MIN(N, NtgroupTimes - 1);
     }
   }
   // we now have : tgroupTimes[Nlo] < start <= tgroupTimes[Nhi]
 
-  // find a matched tgroup starting from Nlo, or return -1
+  // find a matched tgroup starting from Nlo, or return NULL
   // tgroupTimes is a lower bound, the upper bound is tgroupTimes + 86399 (1d)
   for (N = Nlo; N <= Nhi; N++) { 
-    if (start < tgroupTimes[N]) continue;
-    if (start > tgroupTimes[N] + 86399) continue; // XXX since they are sorted, probably do not need this check
-    return (N*Nphotcodes + Ns); // we have the tgroupTime and photcodes, so we can construct the tgroup index
-  }
-  return (-1);
+    if (start < tgroupTimes[N][0].start) continue;
+    if (start > tgroupTimes[N][0].start + 86399) continue; // XXX since they are sorted, probably do not need this check
+    // N is the tgroupTimes entry we want.
+    // double-check that tgroupTimes[N][0].byCode[Ns].photcode = ecode?
+    return &tgroupTimes[N][0].byCode[Ns];
+  }
+  return (NULL);
 }
 
 void setMcalFromTGroups () {
 
-  off_t i, j, im, Nimage;
+  off_t Nimage;
   Image *image;
 
@@ -282,22 +348,26 @@
   // copy the tgroup results to the images.  set the tgroup Mcal to 0.0 since we have moved its
   // impact to the images
-  for (i = 0; i < Ntgroup; i++) {
-    for (j = 0; j < tgroup[i].Nimage; j++) {
-      im = tgroup[i].image[j];
-      double Mgrp = (TGROUP_FIT_AIRMASS) ? tgroup[i].McalPSF + tgroup[i].dKlam*(image[im].secz - 1.0) : tgroup[i].McalPSF;
-      image[im].McalPSF    += Mgrp;
-      image[im].McalAPER   += Mgrp;
-      image[im].dMcal       = tgroup[i].dMcal;
-      image[im].McalChiSq   = tgroup[i].McalChiSq;
-      image[im].dMagSys     = tgroup[i].dMsys;
-      image[im].nFitPhotom  = tgroup[i].nFitPhotom;
-      image[im].flags 	   |= (tgroup[i].flags & ID_IMAGE_PHOTOM_FEW);
-      image[im].flags 	   |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR);
-
-      fprintf (stderr, "TG to IMAGE: %f -> %f (%d)\n", tgroup[i].McalPSF, image[im].McalPSF, tgroup[i].nFitPhotom);
-    }
-    tgroup[i].McalPSF  = 0.0;
-    tgroup[i].McalAPER = 0.0;
-    tgroup[i].dKlam    = 0.0;
+  for (off_t k = 0; k < NtgroupTimes; k++) {
+    TGroup *tgroup = tgroupTimes[k][0].byCode;
+    for (off_t i = 0; i < tgroupTimes[k][0].nCode; i++) {
+      for (off_t j = 0; j < tgroup[i].Nimage; j++) {
+	off_t im = tgroup[i].image[j];
+	double Mgrp = (TGROUP_FIT_AIRMASS) ? tgroup[i].McalPSF + tgroup[i].dKlam*(image[im].secz - 1.0) : tgroup[i].McalPSF;
+	image[im].McalPSF    += Mgrp;
+	image[im].McalAPER   += Mgrp;
+	image[im].dMcal       = tgroup[i].dMcal;
+	image[im].McalChiSq   = tgroup[i].McalChiSq;
+	image[im].dMagSys     = tgroup[i].dMsys;
+	image[im].nFitPhotom  = tgroup[i].nFitPhotom;
+	image[im].flags 	   |= (tgroup[i].flags & ID_IMAGE_PHOTOM_FEW);
+	image[im].flags 	   |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR);
+	image[im].flags 	   |= (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL);
+	
+	fprintf (stderr, "TG to IMAGE: %f +/- %f -> %f (%d)\n", tgroup[i].McalPSF, tgroup[i].dMcal, image[im].McalPSF, tgroup[i].nFitPhotom);
+      }
+      tgroup[i].McalPSF  = 0.0;
+      tgroup[i].McalAPER = 0.0;
+      tgroup[i].dKlam    = 0.0;
+    }
   }      
 }
@@ -310,8 +380,8 @@
   if (!TGROUP_ZEROPT) return;
 
-  ALLOCATE (MeasureToTGroup, off_t *, Ncatalog);
+  ALLOCATE (MeasureToTGroup, TGroup **, Ncatalog);
   for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (MeasureToTGroup[i], off_t, MAX (catalog[i].Nmeasure, 1));
-    for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToTGroup[i][j] = -1;
+    ALLOCATE (MeasureToTGroup[i], TGroup *, MAX (catalog[i].Nmeasure, 1));
+    for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToTGroup[i][j] = NULL;
   }
 }
@@ -360,6 +430,9 @@
   }
   fprintf (stderr, "Matched %d detections to tgroups\n", Nmatch);
-  for (i = 0; i < Ntgroup; i++) {
-    fprintf (stderr, "  TGROUP time %d, Nmeasure: %d, Nimage: %d\n", tgroup[i].start, (int) tgroup[i].Nmeasure, (int) tgroup[i].Nimage);
+  for (i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      fprintf (stderr, "  TGROUP time %.0f photcode %d, Nmeasure: %d, Nimage: %d\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, (int) tgroup[j].Nmeasure, (int) tgroup[j].Nimage);
+    }
   }
   return (TRUE);
@@ -377,6 +450,7 @@
   }
 
-  off_t grpID = ImageToTGroup[idx];
-  if (grpID < 0) {
+  TGroup *myGroup = ImageToTGroup[idx];
+  TGTimes *parent = (TGTimes *) myGroup->parent;
+  if (!myGroup) {
     // Image *image = getimage(idx);
     // fprintf (stderr, "unmatched image %s\n", image[0].name);
@@ -391,5 +465,6 @@
     // unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1);
     unsigned int imageStart = image[0].tzero;
-    if ((imageStart < tgroup[grpID].start) || (imageStart > tgroup[grpID].stop)) {
+    // XXX to do this check I need a pointer to the group master
+    if ((imageStart < parent->start) || (imageStart > parent->stop)) {
       fprintf (stderr, "error in image to tgroup match\n");
       abort();
@@ -399,14 +474,14 @@
   // this measurement is on one of my tgroups, mark it as mine.
   catalog[cat].measureT[meas].myDet = TRUE;
-  MeasureToTGroup[cat][meas] = grpID;
-
-  tgroup[grpID].catalog[tgroup[grpID].Nmeasure] = cat;
-  tgroup[grpID].measure[tgroup[grpID].Nmeasure] = meas;
-  tgroup[grpID].Nmeasure ++;
-
-  if (tgroup[grpID].Nmeasure == tgroup[grpID].NMEASURE) {
-    tgroup[grpID].NMEASURE += 1000;
-    REALLOCATE (tgroup[grpID].catalog, off_t, tgroup[grpID].NMEASURE);
-    REALLOCATE (tgroup[grpID].measure, off_t, tgroup[grpID].NMEASURE);
+  MeasureToTGroup[cat][meas] = myGroup;
+
+  myGroup->catalog[myGroup->Nmeasure] = cat;
+  myGroup->measure[myGroup->Nmeasure] = meas;
+  myGroup->Nmeasure ++;
+
+  if (myGroup->Nmeasure == myGroup->NMEASURE) {
+    myGroup->NMEASURE += 1000;
+    REALLOCATE (myGroup->catalog, off_t, myGroup->NMEASURE);
+    REALLOCATE (myGroup->measure, off_t, myGroup->NMEASURE);
   }
   return;
@@ -415,5 +490,4 @@
 float getMgrp (off_t meas, int cat, float airmass) {
 
-  off_t i;
   float value;
 
@@ -421,9 +495,9 @@
 
   // unassigned measurements belong to simple chips
-  i = MeasureToTGroup[cat][meas];
-  if (i == -1) return (0.0);
-
-  // if (tgroup[i].flags & IMAGE_BAD) return (NAN);  
-  value = tgroup[i].McalPSF + tgroup[i].dKlam*(airmass - 1.0);
+  TGroup *myGroup = MeasureToTGroup[cat][meas];
+  if (!myGroup) return (0.0);
+
+  // if (myGroup->flags & IMAGE_BAD) return (NAN);  
+  value = myGroup->McalPSF + myGroup->dKlam*(airmass - 1.0);
   return (value);
 }
@@ -434,10 +508,12 @@
 
   // unassigned measurements belong to simple chips
-  int i = MeasureToTGroup[cat][meas];
-  if (i == -1) return (0);
-
-  return (tgroup[i].flags);
-}
-
+  TGroup *myGroup = MeasureToTGroup[cat][meas];
+  if (!myGroup) return (0);
+
+  return (myGroup->flags);
+}
+
+// XXX remove the lists and replace with StatDataSet structs?
+// XXX move the StatDataSetAlloc up?
 typedef struct {
   int Nfew;
@@ -526,10 +602,10 @@
 static int nextTGroup = 0;
 
-// we have an array of tgroups (tgroup, Ntgroup).  we need to hand out tgroups one at a time to
+// we have an array of (tgroupTimes, NtgroupTimes).  we need to hand out tgroupTimes one at a time to
 // the worker threads as they need
 off_t getNextTGroupForThread () {
 
   pthread_mutex_lock (&setMgrp_mutex);
-  if (nextTGroup >= Ntgroup) {
+  if (nextTGroup >= NtgroupTimes) {
     pthread_mutex_unlock (&setMgrp_mutex);
     return (-1);
@@ -544,5 +620,5 @@
 int setMgrp (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
 
-  off_t i, N, Nmax;
+  off_t N;
   Image *image;
 
@@ -564,7 +640,10 @@
   }
 
-  Nmax = 0;
-  for (i = 0; i < Ntgroup; i++) {
-    Nmax = MAX (Nmax, tgroup[i].Nmeasure);
+  off_t Nmax = 0;
+  for (off_t i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      Nmax = MAX (Nmax, tgroup[j].Nmeasure);
+    }
   }
 
@@ -575,7 +654,10 @@
   // if (PLOTSTUFF) PLOTDELAY = 0.0;
 
-  for (i = 0; i < Ntgroup; i++) {
-    setMgrp_tgroup (&tgroup[i], i, image, catalog, &info, flatcorr);
-    fprintf (stderr, "TGROUP Mcal: %f, dK: %f, limiting negative clouds to %f\n", tgroup[i].McalPSF, tgroup[i].dKlam, CLOUD_TOLERANCE);
+  for (off_t i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      setMgrp_tgroup (&tgroup[j], i, image, catalog, &info, flatcorr);
+      fprintf (stderr, "TGROUP time %.0f photcode %d Mcal: %f, dK: %f, dMcal: %f, limiting negative clouds to %f\n", ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, CLOUD_TOLERANCE);
+    }
   }
   SetMgrpInfoFree (&info);
@@ -601,34 +683,40 @@
 # define NIGHT_GOOD_FRACTION 0.1
 
-// 'tgroup' is a pointer to the current tgroup of interest (Nmos)
+// 'tgroup' is a pointer to the current tgroup of interest (entry Ngrp)
 int setMgrp_tgroup (TGroup *myTGroup, off_t Ngrp, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr) {
 
   off_t j, NimageReal;
 
-  StatType stats;
-  liststats_setmode (&stats, "INNER_WTMEAN");
-
-  double *psfMagList    = info->psfMagList;
-  double *psfAirList    = info->psfAirList;
-  double *psfErrList    = info->psfErrList;
-  double *kronMagList   = info->kronMagList;
-  double *kronErrList   = info->kronErrList;
-  double *brightMagList = info->brightMagList;
-  double *brightErrList = info->brightErrList;
+  /* on PoorImages run, skip good images */
+  if (info->PoorImages) {
+    int bad = myTGroup->flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
+    if (!bad) return TRUE;
+  }      
+
+  int fitOrder = TGROUP_FIT_AIRMASS ? 1 : 0;
+
+  // allocate the arrays used by the fit1d_irls analysis below.
+  // the last 0 means we are NOT using bootstrap analysis here
+
+  FitDataSet psfStars, kronStars, brightStars;
+  FitDataSetAlloc (&psfStars,    myTGroup->Nmeasure, fitOrder, 0); 
+  FitDataSetAlloc (&kronStars,   myTGroup->Nmeasure, fitOrder, 0); 
+  FitDataSetAlloc (&brightStars, myTGroup->Nmeasure, fitOrder, 0); 
+
+  // add prior constraint on airmass:
+  if (1) {
+    FitDataSetAddPriors (&psfStars);
+    psfStars.bPriorValue[1] = 0.0;  // note that we are fitting relative to the nominal slope
+    psfStars.bPriorSigma[1] = 0.04; // XXX this prior sigma needs to be user-configured
+  }
 
   // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used
   getimages (&NimageReal, NULL);
 
-  /* on PoorImages run, skip good images */
-  if (info->PoorImages) {
-    int bad = myTGroup[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
-    if (!bad) return TRUE;
-  }      
-
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
   int testImage = FALSE;
-  testImage |= (abs(myTGroup[0].start - 1245283200) < 10);
-  // testImage |= (abs(myTGroup[0].start - 1324103823) < 10);
+  TGTimes *tgroup = (TGTimes *) myTGroup->parent;
+  testImage |= (abs(tgroup->start - 1245283200) < 10);
 
   FILE *fout = NULL;
@@ -732,12 +820,16 @@
     assert (Nbright >= 0);
 
-    psfMagList[N]  = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
-    psfAirList[N]  = (catalog[c].measureT[m].airmass - 1.0);
-    psfErrList[N]  = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-    kronMagList[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
-    kronErrList[N] = psfErrList[N];
+    psfStars.alldata-> yVector[N] = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
+    psfStars.alldata-> xVector[N] = (catalog[c].measureT[m].airmass - 1.0);
+    psfStars.alldata->dyVector[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+
+    kronStars.alldata-> yVector[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
+    kronStars.alldata-> xVector[N] = psfStars.alldata-> xVector[N];
+    kronStars.alldata->dyVector[N] = psfStars.alldata->dyVector[N];
+
     if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
-      brightMagList[Nbright] = psfMagList[N];
-      brightErrList[Nbright] = psfErrList[N];
+      brightStars.alldata-> yVector[Nbright] = psfStars.alldata-> yVector[N];
+      brightStars.alldata-> xVector[Nbright] = psfStars.alldata-> xVector[N];
+      brightStars.alldata->dyVector[Nbright] = psfStars.alldata->dyVector[N];
       Nbright ++;
     }
@@ -750,7 +842,9 @@
     int mark = (N < NIGHT_TOOFEW) || (N < NIGHT_GOOD_FRACTION*myTGroup->Nmeasure);
     if (mark) {
-      // XXX is this right (effect of subset?) : image[tgroup[Ngrp].image[0]].name
-      if (VERBOSE2) { fprintf (stderr, "marked tgroup %s ("OFF_T_FMT"), (%d < %d) || (%d < %f*"OFF_T_FMT")\n", image[tgroup[Ngrp].image[0]].name, Ngrp, N, NIGHT_TOOFEW, N, NIGHT_GOOD_FRACTION, myTGroup->Nmeasure); }
-      myTGroup[0].flags |= ID_IMAGE_PHOTOM_FEW;
+      if (VERBOSE2) {
+	TGTimes *mygroup = (TGTimes *) myTGroup->parent;
+	fprintf (stderr, "marked tgroup %f,%d, (%d < %d) || (%d < %f*"OFF_T_FMT")\n", ohana_sec_to_mjd(mygroup->start), myTGroup->photcode, N, NIGHT_TOOFEW, N, NIGHT_GOOD_FRACTION, myTGroup->Nmeasure);
+      }
+      myTGroup->flags |= ID_IMAGE_PHOTOM_FEW;
       info->Nfew ++;
       if (testImage) {
@@ -758,48 +852,54 @@
       }
     } else {
-      myTGroup[0].flags &= ~ID_IMAGE_PHOTOM_FEW;
-    }
-  }
-
-  double dk = 0.0;
-  if (TGROUP_FIT_AIRMASS) {
-    liststats_fit1d (psfMagList, psfErrList, psfAirList, N, &stats, &dk);
-  } else {
-    liststats (psfMagList, psfErrList, NULL, N, &stats);
-  }
-  if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mgrp: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
+      myTGroup->flags &= ~ID_IMAGE_PHOTOM_FEW;
+    }
+  }
+
+  // do anything special with identified good nights?
+  // (myTGroup->flags & ID_IMAGE_PHOTOM_UBERCAL)
+
+  fit1d_irls (&psfStars, N);
+
+  // if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mgrp: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
 
   // for now, I have no reason to measure these separately for camera-level images
-  myTGroup[0].McalPSF    = stats.mean;
-  myTGroup[0].dMcal      = stats.error;
-  myTGroup[0].McalChiSq  = stats.chisq;
-  myTGroup[0].nFitPhotom = N;
-  myTGroup[0].dKlam      = dk;
-
-  // drop this for tgroups?
-  liststats (kronMagList, kronErrList, NULL, N, &stats);
-  myTGroup[0].McalAPER   = stats.mean;
+  myTGroup->McalPSF    = psfStars.bSaveArray[0][0];
+  myTGroup->dMcal      = psfStars.bSigma[0];
+  myTGroup->McalChiSq  = psfStars.chisq;
+  myTGroup->nFitPhotom = psfStars.Nmeas;
+  myTGroup->dKlam      = psfStars.bSaveArray[1][0];
+
+  // XXX drop this for tgroups?
+  // fit1d_irls (&kronStars, N);    // does it make sense to calculate a zero point based kron mags?
+  // myTGroup->McalAPER   = stats.mean;
 
   if (testImage) {
-    fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Ngrp, myTGroup[0].start, stats.mean, stats.error, myTGroup[0].nFitPhotom);
+    TGTimes *parent = (TGTimes *) myTGroup->parent;
+    fprintf (stderr, "test night %.0f aper: %f %f %d ... ", ohana_sec_to_mjd(parent->start), myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom);
   }
 
   if (PLOTSTUFF) {
-    fprintf (stderr, "Mgrp: %6.3f %6.3f +/- %6.3f %5d %5d | %s\n", stats.mean, stats.median, stats.sigma, stats.Nmeas, N, image[tgroup[Ngrp].image[0]].name);
-    plot_setMcal (psfMagList, N, &stats, CLOUD_TOLERANCE);
+    fprintf (stderr, "Mgrp: %6.3f +/- %6.3f %5d of %5d | %.0f\n", myTGroup->McalPSF, myTGroup->dMcal, myTGroup->nFitPhotom, N, ohana_sec_to_mjd(tgroup->start));
+    plot_setMcal (psfStars.alldata->yVector, N);
   }
 
   // bright end scatter
-  liststats (brightMagList, brightErrList, NULL, Nbright, &stats);
-  myTGroup[0].dMsys = stats.sigma;
+  // XXX this does not make sense: I need to apply the airmass slope calculated above first
+  // redo this by calculating the corrected bright stars mags 
+  // fit1d_irls (&brightStars, N);  // does it make sense to calculate a separate bright star zero point?
+  // myTGroup->dMsys = stats.sigma;
 
   // keep this??
-  if (myTGroup[0].McalPSF < -CLOUD_TOLERANCE) {
-    myTGroup[0].McalPSF = 0.0;
+  if (myTGroup->McalPSF < -CLOUD_TOLERANCE) {
+    myTGroup->McalPSF = 0.0;
   }
 
   if (testImage) {
-    fprintf (stderr, "%f %f  :  %f\n", myTGroup[0].McalPSF, myTGroup[0].dMsys, myTGroup[0].McalChiSq);
-  }
+    fprintf (stderr, "%f %f  :  %f\n", myTGroup->McalPSF, myTGroup->dMsys, myTGroup->McalChiSq);
+  }
+
+  FitDataSetFree (&brightStars);
+  FitDataSetFree (&kronStars);
+  FitDataSetFree (&psfStars);
 
   return TRUE;
@@ -822,6 +922,9 @@
 
   off_t Nmax = 0;
-  for (i = 0; i < Ntgroup; i++) {
-    Nmax = MAX (Nmax, tgroup[i].Nmeasure);
+  for (off_t i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      Nmax = MAX (Nmax, tgroup[j].Nmeasure);
+    }
   }
 
@@ -881,5 +984,5 @@
 	     threadinfo[i].info.Nrel, 
 	     threadinfo[i].info.Nsys,
-      	     threadinfo[i].info.Nskip);
+	     threadinfo[i].info.Nskip);
     SetMgrpInfoAccum (&summary, &threadinfo[i].info);
   }
@@ -924,6 +1027,9 @@
     Image *image = threadinfo->image;
 
-    setMgrp_tgroup (&tgroup[i], i, image, catalog, &results, flatcorr);
-    SetMgrpInfoAccum (&threadinfo->info, &results);
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      setMgrp_tgroup (&tgroup[j], i, image, catalog, &results, flatcorr);
+      SetMgrpInfoAccum (&threadinfo->info, &results);
+    }
   }
 
@@ -932,4 +1038,5 @@
 }
 
+# if (0)
 StatType statsTGroupM (Catalog *catalog) {
   OHANA_UNUSED_PARAM(catalog);
@@ -943,13 +1050,17 @@
   if (!TGROUP_ZEROPT) return (stats);
 
+  // XXX do this by time & photcode?
   ALLOCATE (list, double, Ntgroup);
   ALLOCATE (dlist, double, Ntgroup);
 
   n = 0;
-  for (i = 0; i < Ntgroup; i++) {
-    if (tgroup[i].flags & IMAGE_BAD) continue;
-    list[n] = tgroup[i].McalPSF;
-    dlist[n] = 1;
-    n++;
+  for (int i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      if (tgroup[j].flags & IMAGE_BAD) continue;
+      list[n] = tgroup[j].McalPSF;
+      dlist[n] = 1;
+      n++;
+    }
   }
 
@@ -1068,4 +1179,5 @@
   return (stats);
 }
+# endif
 
 /* mark tgroup if: abs(Mcal - <Mcal>) too large, dMcal too large */
@@ -1080,14 +1192,17 @@
   if (VERBOSE) fprintf (stderr, "marking poor tgroups\n");
 
-  ALLOCATE (mlist, double, Ntgroup);
-  ALLOCATE (slist, double, Ntgroup);
-  ALLOCATE (dlist, double, Ntgroup);
-
-  for (i = N = 0; i < Ntgroup; i++) {
-    if (tgroup[i].flags & IMAGE_BAD) continue;
-    mlist[N] = tgroup[i].McalPSF;
-    slist[N] = tgroup[i].dMcal;
-    dlist[N] = 1;
-    N++;
+  ALLOCATE (mlist, double, NtgroupTimes*Nphotcodes);
+  ALLOCATE (slist, double, NtgroupTimes*Nphotcodes);
+  ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
+
+  for (i = N = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      if (tgroup[j].flags & IMAGE_BAD) continue;
+      mlist[N] = tgroup[j].McalPSF;
+      slist[N] = tgroup[j].dMcal;
+      dlist[N] = 1;
+      N++;
+    }
   }
 
@@ -1096,32 +1211,36 @@
 
   liststats (mlist, dlist, NULL, N, &stats);
-  MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
+  MaxOffset = MAX (NIGHT_OFFSET, 2*stats.sigma);
   MedOffset = stats.median;
 
   liststats (slist, dlist, NULL, N, &stats);
-  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
+  MaxScatter = MAX (NIGHT_SCATTER, 2*stats.median);
   fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
   
   Nmark = Nscatter = Noffset = 0;
-  for (i = 0; i < Ntgroup; i++) {
-    // if we are keeping ubercal sacrosanct, then we should not be allowed to break them...
-    if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
-
-    if (tgroup[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue;
-
-    mark = FALSE;
-    if (tgroup[i].dMcal > MaxScatter) {
-      mark = TRUE;
-      Nscatter ++;
-    }
-    if (fabs(tgroup[i].McalPSF - MedOffset) > MaxOffset) {
-      mark = TRUE;
-      Noffset ++;
-    }
-    if (mark) { 
-      Nmark ++;
-      tgroup[i].flags |= ID_IMAGE_PHOTOM_POOR;
-    } else {
-      tgroup[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
+  for (i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      // if we are keeping ubercal sacrosanct, then we should not be allowed to break them...
+      if (KEEP_UBERCAL && (tgroup[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
+      
+      if (tgroup[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue;
+      
+      mark = FALSE;
+      if (tgroup[i].dMcal > MaxScatter) {
+	mark = TRUE;
+	Nscatter ++;
+      }
+      if (fabs(tgroup[i].McalPSF - MedOffset) > MaxOffset) {
+	mark = TRUE;
+	Noffset ++;
+      }
+      if (mark) { 
+	Nmark ++;
+	tgroup[i].flags |= ID_IMAGE_PHOTOM_POOR;
+      } else {
+	tgroup[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
+      }
+      // XXX : unset the ubercal flag if we mark the night as bad?
     }
   }
@@ -1136,5 +1255,5 @@
 void plot_tgroup_fields (Catalog *catalog) {
 
-  off_t i, j, m, c, N, Nimage;
+  off_t m, c, Nimage;
   double *xlist, *ylist;
   char string[64];
@@ -1146,7 +1265,10 @@
   getimages (&Nimage, NULL);
 
-  N = 0;
-  for (i = 0; i < Ntgroup; i++) {
-    N = MAX (N, tgroup[i].Nmeasure);
+  off_t N = 0;
+  for (off_t i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      N = MAX (N, tgroup[j].Nmeasure);
+    }
   }
 
@@ -1154,20 +1276,22 @@
   ALLOCATE (ylist, double, N);
 
-  for (i = 0; i < Ntgroup; i++) {
+  for (off_t k = 0; k < NtgroupTimes; k++) {
     N = 0;
-    for (j = 0; j < tgroup[i].Nmeasure; j++) {
+    TGroup *tgroup = tgroupTimes[k][0].byCode;
+    for (off_t i = 0; i < tgroupTimes[i][0].nCode; i++) {
+      for (off_t j = 0; j < tgroup[i].Nmeasure; j++) {
+	
+	m = tgroup[i].measure[j];
+	c = tgroup[i].catalog[j];
       
-      m = tgroup[i].measure[j];
-      c = tgroup[i].catalog[j];
-      
-      if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
-
-      // ave = catalog[c].measureT[m].averef;
-      xlist[N] = catalog[c].measureT[m].R;
-      ylist[N] = catalog[c].measureT[m].D;
-      N++;
-    }
-  
-    sprintf (string, "TGroup "OFF_T_FMT,  i);
+	if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
+
+	// ave = catalog[c].measureT[m].averef;
+	xlist[N] = catalog[c].measureT[m].R;
+	ylist[N] = catalog[c].measureT[m].D;
+	N++;
+      }
+    }
+    sprintf (string, "TGroup "OFF_T_FMT,  k);
     plot_defaults (&graphdata);
     plot_list (&graphdata, xlist, ylist, N, string, NULL);
@@ -1186,12 +1310,17 @@
   if (!TGROUP_ZEROPT) return;
 
-  ALLOCATE (xlist, double, Ntgroup);
-  ALLOCATE (dlist, double, Ntgroup);
-  ALLOCATE (Mlist, double, Ntgroup);
-
-  for (i = 0; i < Ntgroup; i++) {
-    Mlist[i] = tgroup[i].McalPSF;
-    dlist[i] = tgroup[i].dMcal;
-    xlist[i] = tgroup[i].start;
+  ALLOCATE (xlist, double, NtgroupTimes*Nphotcodes);
+  ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
+  ALLOCATE (Mlist, double, NtgroupTimes*Nphotcodes);
+
+  int Npts = 0;
+  for (i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      Mlist[Npts] = tgroup[j].McalPSF;
+      dlist[Npts] = tgroup[j].dMcal;
+      xlist[Npts] = tgroupTimes[i][0].start;
+      Npts ++;
+    }
   }
 
@@ -1201,9 +1330,9 @@
   graphdata.ymin = PlotdMmin;
   graphdata.ymax = PlotdMmax;
-  plot_list (&graphdata, xlist, Mlist, Ntgroup, "airmass vs Mcal", "%s.airmass.png", OUTROOT);
+  plot_list (&graphdata, xlist, Mlist, Npts, "airmass vs Mcal", "%s.airmass.png", OUTROOT);
   plot_defaults (&graphdata);
   graphdata.size = 1.5;
   graphdata.ptype = 7;
-  plot_list (&graphdata, Mlist, dlist, Ntgroup, "Mcal vs dMcal", "%s.MdM.png", OUTROOT);
+  plot_list (&graphdata, Mlist, dlist, Npts, "Mcal vs dMcal", "%s.MdM.png", OUTROOT);
 
 # define NBIN 200
@@ -1214,8 +1343,12 @@
   for (i = 0; i < NBIN; i++) xlist[i] = 0.00005*i;
   bzero (Mlist, NBIN*sizeof(double));
-  for (i = 0; i < Ntgroup; i++) {
-    bin = tgroup[i].dMcal / 0.00005;
-    bin = MAX (0, MIN (NBIN - 1, bin));
-    Mlist[bin] += 1.0;
+
+  for (i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      bin = tgroup[j].dMcal / 0.00005;
+      bin = MAX (0, MIN (NBIN - 1, bin));
+      Mlist[bin] += 1.0;
+    }
   }
   plot_defaults (&graphdata);
Index: trunk/Ohana/src/relphot/src/args.c
===================================================================
--- trunk/Ohana/src/relphot/src/args.c	(revision 41451)
+++ trunk/Ohana/src/relphot/src/args.c	(revision 41453)
@@ -579,4 +579,5 @@
   
   freeTGroups();
+  free_error();
 
   ohana_memcheck (VERBOSE);
Index: trunk/Ohana/src/relphot/src/fit1d_irls.c
===================================================================
--- trunk/Ohana/src/relphot/src/fit1d_irls.c	(revision 41453)
+++ trunk/Ohana/src/relphot/src/fit1d_irls.c	(revision 41453)
@@ -0,0 +1,369 @@
+# include "relphot.h"
+
+// These should probably be tunable:
+# define MAX_ITERATIONS 10
+# define FIT_TOLERANCE 1e-4
+# define FLT_TOLERANCE 1e-6
+# define WEIGHT_THRESHOLD 0.3
+
+# define IS_OLS TRUE
+
+void fit1d_copy_results (double **srcArray, double **tgtArray, int entry, int nterm);
+void get_yfit (FitDataSet *dataset, int Npoints);
+int fit1d_least_squares (FitDataSet *dataset, FitDataType *data, int isOLS, int Npoints);
+
+// 1D fit with arbitrary order.
+// NOTE: the number of points considered for the fit (Npoints)
+// may be different from the number of points stored in the array
+int fit1d_irls (FitDataSet *dataset, int Npoints) {
+
+  // initialize various things:
+  dataset->Nmeas = 0;
+  dataset->min = NAN;
+  dataset->max = NAN;
+  dataset->sigma = NAN;
+  dataset->chisq = NAN;
+  for (int i = 0; i < dataset->nterm; i++) {
+    dataset->bSigma[i] = 0;
+  }
+
+  // initial fit with OLS, results in dataset->bArray
+  if (!fit1d_least_squares (dataset, dataset->alldata, IS_OLS, Npoints)) return FALSE;
+  
+  // apply fit to data and save in yFitVector & yOffVector
+  get_yfit (dataset, Npoints);
+
+  int converged = FALSE;
+  for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) {
+    
+    for (int i = 0; i < Npoints; i++) {
+      // we are only including the formal error, not the weight in the definition of wt[]
+      dataset->wtIRLS[i] = weight_cauchy (dataset->yOffVector[i] / dataset->alldata->dyVector[i]);
+    }
+    
+    // save this solution to entry 0
+    fit1d_copy_results (dataset->bArray, dataset->bSaveArray, 0, dataset->nterm);
+
+    // recalculate the least squares, but now apply the modified weight above (arg 3 = isOLS)
+    if (!fit1d_least_squares (dataset, dataset->alldata, !IS_OLS, Npoints)) {
+      // restore the saved solution
+      fit1d_copy_results (dataset->bSaveArray, dataset->bArray, 0, dataset->nterm);
+      break;
+    }
+    
+    // apply fit to data and save in yFitVector & yOffVector (used above to calculate wtIRLS
+    get_yfit (dataset, Npoints);
+
+    // has the fit changed in a meaningful way?
+    converged = TRUE;
+    for (int i = 0; i < dataset->nterm; i++) {
+      if ((fabs(dataset->bArray[i][0] - dataset->bSaveArray[i][0]) > FIT_TOLERANCE * fabs(dataset->bArray[i][0])) && 
+	  (fabs(dataset->bArray[i][0] - dataset->bSaveArray[i][0]) > FLT_TOLERANCE))
+	converged = FALSE;
+    }
+  }
+
+  // save this solution to entry 0
+  fit1d_copy_results (dataset->bArray, dataset->bSaveArray, 0, dataset->nterm);
+      
+  // calculate the weight thresholds to mask the bad points:
+  for (int i = 0; i < Npoints; i++) {
+    dataset->wtIRLS[i] = weight_cauchy (dataset->yOffVector[i] / dataset->alldata->dyVector[i]);
+    dataset->tmpVector[i] = dataset->wtIRLS[i];
+  }
+  int midpt = 0.5 * Npoints;
+  dsort (dataset->tmpVector, Npoints);
+  double WtMedian = (Npoints % 2) ? dataset->tmpVector[midpt] : 0.5*(dataset->tmpVector[midpt] + dataset->tmpVector[midpt-1]);
+  double WtThreshold = WEIGHT_THRESHOLD * WtMedian;
+
+  // points with too low weight are ignored in the bootstrap analysis below
+  int Nkeep = 0;
+  double dChi = 0.0;
+  double dSig = 0.0;
+  for (int i = 0; i < Npoints; i++) {
+    if ((dataset->wtIRLS[i] < WtThreshold) || !isfinite(dataset->alldata->yVector[i])) {
+      continue; // skip the masked points
+    }
+    // save values with sufficient weight for error analysis below (either basic OLS or bootstrap)
+    dataset->keepdata->yVector[Nkeep]  = dataset->alldata->yVector[i];
+    dataset->keepdata->dyVector[Nkeep] = dataset->alldata->dyVector[i];
+    if (dataset->alldata->xVector) {
+      dataset->keepdata->xVector[Nkeep] = dataset->alldata->xVector[i]; // only defined if order > 0
+    }
+    Nkeep ++;
+    
+    // record min, max, chisq, sigma, Nmeas for the unmasked points
+    dataset->min = isfinite(dataset->min) ? MIN(dataset->alldata->yVector[i], dataset->min) : dataset->alldata->yVector[i];
+    dataset->max = isfinite(dataset->max) ? MAX(dataset->alldata->yVector[i], dataset->max) : dataset->alldata->yVector[i];
+    // fprintf (stderr, "%d : %f %f %f\n", i, dataset->alldata->yVector[i], dataset->min, dataset->max);
+    double dValue2 = SQ(dataset->yOffVector[i]);
+    dChi += dValue2 / SQ (dataset->alldata->dyVector[i]);
+    dSig += dValue2;
+  }
+  dataset->Nmeas = Nkeep;
+  dataset->chisq = dChi / (Nkeep - 1);
+  dataset->sigma = sqrt (dSig / (Nkeep - 1));
+
+  if (dataset->Nbootstrap) {
+    int Nboot = 0;
+    for (int iboot = 0; iboot < dataset->Nbootstrap; iboot++) {
+    
+      // resample
+      for (int i = 0; i < Nkeep; i++) {
+	// I need to draw Npoints random entries from 'points' with replacement:
+	int N = Nkeep * drand48();
+	dataset->sample->yVector[i]  = dataset->keepdata-> yVector[N];
+	dataset->sample->dyVector[i] = dataset->keepdata->dyVector[N];
+	if (dataset->keepdata->xVector) {
+	  dataset->sample->xVector[i] = dataset->keepdata->xVector[N]; // only defined if order > 0
+	}
+      }
+
+      if (!fit1d_least_squares (dataset, dataset->sample, IS_OLS, Nkeep)) continue;
+      // fprintf (stderr, "%d = %f : ", iboot, dataset->bArray[0][0]);
+
+      // save this solution in entry Nboot
+      fit1d_copy_results (dataset->bArray, dataset->bBootArray, Nboot, dataset->nterm);
+      // fprintf (stderr, "%d = %f : ", Nboot, dataset->bBootArray[0][Nboot]);
+      Nboot ++;
+    }
+
+    // loop over the nterm values
+    for (int i = 0; i < dataset->nterm; i++) {
+      // copy the bBootArray values to a tmp array
+      for (int j = 0; j < Nboot; j++) {
+	dataset->tmpVector[j] = dataset->bBootArray[i][j];
+      }
+      dsort (dataset->tmpVector, Nboot);
+      double Slo = VectorFractionInterpolate (dataset->tmpVector, 0.158655, Nboot);
+      double Shi = VectorFractionInterpolate (dataset->tmpVector, 0.841345, Nboot);
+      dataset->bSigma[i] = (Shi - Slo) / 2.0;
+    }
+  }
+
+  // calculate the formal fit error from the covariance matrix.  if we do not do
+  // bootstrap, this is the reported parameter error.  Even if we do use bootstrap,
+  // bootstrap with few values can sometimes yield an excessively-optimistic result for
+  // the error.  Do not let the reported error be smaller than the formal error
+  if (fit1d_least_squares (dataset, dataset->keepdata, IS_OLS, Nkeep)) {
+    for (int i = 0; i < dataset->nterm; i++) {
+      double errvalue = sqrt(dataset->cArray[i][i]);
+      dataset->bSigma[i] = MAX (dataset->bSigma[i], errvalue);
+    }
+  } 
+  return TRUE;
+}
+
+// full least-squares fit function including weights
+// results are returned to dataset->bArray
+int fit1d_least_squares (FitDataSet *dataset, FitDataType *data, int isOLS, int Npoints) { 
+
+  // these are preallocated so repeated evaluations are not slowed down by alloc
+  double *s  = dataset->sumVector;
+  double **b = dataset->bArray;
+  double **c = dataset->cArray;
+
+  int nterm  = dataset->nterm; // number of parameters being fitted, e.g., y = C0 + C1*x -> 1st order, 2 terms
+  int mterm  = dataset->mterm; // number of x^n terms calculated (e.g., for 3rd order polynomial, we need x^0 - x^6, i.e., 7 x^n terms)
+
+  double *x  = data->xVector;
+  double *y  = data->yVector;
+  double *dy = data->dyVector;
+
+  // XXX : the wtIRLS vector has the same sequence as the data vectors
+  double *wt = isOLS ? NULL : dataset->wtIRLS; // use IRLS weight if requested
+
+  // initialize summation registers for current pass
+  memset (s, 0, mterm*sizeof(double));
+  for (int i = 0; i < nterm; i++) {
+    memset (c[i], 0, nterm*sizeof(double));
+    memset (b[i], 0, sizeof(double));
+  }
+
+  // perform linear fit
+  for (int i = 0; i < Npoints; i++, y++, dy++) {
+    if (!finite(*y)) goto next; // skip invalid (NAN) data
+    if (x && !finite(*x)) goto next; // if we are fitting x, skip invalid (NAN) data
+
+    // XXX in the opihi version, I combine weights; in the relphot version I have a second weight vector to apply
+    double dY = wt ? *wt / SQ(*dy) : 1.0 / SQ(*dy);
+
+    double  X =  1*dY;
+    double  Y = *y*dY;
+
+    for (int j = 0; j < nterm; j++) {
+      s[j] += X;
+      b[j][0] += Y;
+      if (x) {
+	// x is NULL is 0-order fit
+	X = X * (*x);
+	Y = Y * (*x);
+      }
+    }
+    for (int j = nterm; j < mterm; j++) {
+      s[j] += X;
+      if (x) { X = X * (*x); }
+    }
+  next:
+    // these elements are optional and NULL if not used:
+    if (x) { x++; }
+    if (wt) { wt++; }
+  }
+  for (int i = 0; i < nterm; i++) {
+    for (int j = 0; j < nterm; j++) {
+      c[i][j] = s[i + j];
+      if (dataset->bPriorValue && (i == j)) {
+	if (!isnan(dataset->bPriorValue[i])) {
+	  c[i][i] += 1.0 / SQ(dataset->bPriorSigma[i]);
+	  b[i][0] += dataset->bPriorValue[i] / SQ(dataset->bPriorSigma[i]);
+	}
+      }
+    }
+  }
+
+  if (!dgaussjordan (c, b, nterm, 1)) {
+    gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
+    return FALSE;
+  }
+  return TRUE;
+}
+
+void get_yfit (FitDataSet *dataset, int Npoints) {
+
+  double *x = dataset->alldata->xVector;
+  double *y = dataset->alldata->yVector;
+
+  // XXX: keep this or not?
+  // double *yfit = dataset->yFitVector;
+  double *yoff = dataset->yOffVector;
+
+  int nterm = dataset->nterm;
+  double **b = dataset->bArray;
+
+  /* generate fitted values */
+  for (int i = 0; i < Npoints; i++, y++, yoff++) {
+    if (x && !finite(*x)) continue;
+    double yfit = 0;
+    double X = 1;
+    for (int j = 0; j < nterm; j++) {
+      yfit += b[j][0]*X;
+      if (x) { X = X * (*x); }
+    }
+    *yoff = *y - yfit;
+    if (x) { x ++; }
+  }
+  return;
+}
+
+void fit1d_copy_results (double **srcArray, double **tgtArray, int entry, int nterm) {
+  for (int i = 0; i < nterm; i++) {
+    tgtArray[i][entry] = srcArray[i][0];
+  }
+  // fprintf (stderr, "%d = %f -> %f: ", entry, srcArray[0][0], tgtArray[0][entry]);
+}
+
+// allocate the vectors for this FitDataType entry (if order == 0, xVector = NULL)
+FitDataType *FitDataTypeAlloc (int order, int N) {
+
+  ALLOCATE_PTR (data, FitDataType, 1);
+
+  ALLOCATE (data->yVector, double, N);
+  ALLOCATE (data->dyVector, double, N);
+  if (order) {
+    ALLOCATE (data->xVector, double, N);
+  } else {
+    data->xVector = NULL;
+  }
+  return data;
+}
+
+// allocate the vectors for this FitDataType entry (if order == 0, xVector = NULL)
+void FitDataTypeFree (FitDataType *data) {
+  if (!data) return;
+  FREE (data->yVector);
+  FREE (data->dyVector);
+  FREE (data->xVector);
+  FREE (data);
+}
+
+void FitDataSetAlloc (FitDataSet *dataset, int Nmax, int order, int Nbootstrap) {
+
+  dataset->order = order;
+  dataset->nterm = order + 1;
+  dataset->mterm = 2*order + 1;
+  dataset->Nbootstrap = Nbootstrap;
+
+  // arrays for actual data values
+  dataset->alldata  = FitDataTypeAlloc (order, Nmax);
+
+  // subset arrays to calculate formal errors
+  dataset->keepdata = FitDataTypeAlloc (order, Nmax);
+
+  dataset->sample   = NULL;  
+  if (Nbootstrap) {
+    // arrays for bootstrap samples
+    dataset->sample   = FitDataTypeAlloc (order, Nmax);
+  } 
+
+  // allocate internal vectors and arrays 
+  ALLOCATE (dataset->bSigma,     double,   dataset->nterm);
+  ALLOCATE (dataset->sumVector,  double,   dataset->mterm);
+  ALLOCATE (dataset->cArray,     double *, dataset->nterm);
+  ALLOCATE (dataset->bArray,     double *, dataset->nterm);
+  ALLOCATE (dataset->bSaveArray, double *, dataset->nterm);
+  ALLOCATE (dataset->bBootArray, double *, dataset->nterm);
+  for (int i = 0; i < dataset->nterm; i++) {
+    ALLOCATE (dataset->cArray[i],      double, dataset->nterm);
+    ALLOCATE (dataset->bArray[i],      double, 1);
+    ALLOCATE (dataset->bSaveArray[i],  double, 1);
+    ALLOCATE (dataset->bBootArray[i],  double, Nbootstrap);
+  }
+
+  ALLOCATE (dataset->yOffVector,  double, Nmax);
+  ALLOCATE (dataset->tmpVector,   double, MAX(Nmax, Nbootstrap));
+  ALLOCATE (dataset->wtIRLS,      double, Nmax);
+
+  // set to NULL; user can activate
+  dataset->bPriorValue = NULL;
+  dataset->bPriorSigma = NULL;
+}
+
+void FitDataSetFree (FitDataSet *dataset) {
+
+  FREE (dataset->bPriorValue);
+  FREE (dataset->bPriorSigma);
+  FREE (dataset->wtIRLS);
+  FREE (dataset->tmpVector);
+  FREE (dataset->yOffVector);
+
+  FitDataTypeFree (dataset->alldata);
+  FitDataTypeFree (dataset->keepdata);
+  if (dataset->Nbootstrap) {
+    FitDataTypeFree (dataset->sample);
+  }
+
+  // allocate internal vectors and arrays 
+  for (int i = 0; i < dataset->nterm; i++) {
+    FREE (dataset->cArray[i]);
+    FREE (dataset->bArray[i]);
+    FREE (dataset->bSaveArray[i]);
+    FREE (dataset->bBootArray[i]);
+  }
+  FREE (dataset->sumVector);
+  FREE (dataset->cArray);
+  FREE (dataset->bArray);
+  FREE (dataset->bSigma);
+  FREE (dataset->bSaveArray);
+  FREE (dataset->bBootArray);
+}  
+
+void FitDataSetAddPriors (FitDataSet *dataset) {
+
+  ALLOCATE (dataset->bPriorValue, double, dataset->nterm);
+  ALLOCATE (dataset->bPriorSigma, double, dataset->nterm);
+
+  for (int i = 0; i < dataset->nterm; i++) {
+    dataset->bPriorValue[i] = NAN;
+    dataset->bPriorSigma[i] = NAN;
+  }
+}
Index: trunk/Ohana/src/relphot/src/initialize.c
===================================================================
--- trunk/Ohana/src/relphot/src/initialize.c	(revision 41451)
+++ trunk/Ohana/src/relphot/src/initialize.c	(revision 41453)
@@ -3,4 +3,5 @@
 RelphotMode initialize (int argc, char **argv) {
 
+  init_error();
   relphot_help (argc, argv);
   ConfigInit (&argc, argv);
Index: trunk/Ohana/src/relphot/src/liststats.c
===================================================================
--- trunk/Ohana/src/relphot/src/liststats.c	(revision 41451)
+++ trunk/Ohana/src/relphot/src/liststats.c	(revision 41453)
@@ -209,11 +209,25 @@
 
   double zp = (M*Z2 - Z1*MZ)*Det;
-  *dk = (MZ*R - M*Z1)*Det;
-
+  double dK = (MZ*R - M*Z1)*Det;
+
+  /* find sigma and chisq */
+  double Mo = 0.0, dM = 0.0;
+  double X2 = 0.0, dS = 0.0;
+  for (int i = 0; i < Npts; i++) {
+    Mo = zp + dK*x[i];
+    M  = SQ (value[i] - Mo);
+    dM = SQ (err[i]);
+    X2 += M / dM;
+    dS += M;
+  }
+  X2 = X2 / (Npts - 1);
+  dS = sqrt (dS / (Npts - 1));
+
+  *dk = dK;
   stats->mean = zp;
-  stats->chisq = NAN;
+  stats->chisq = X2;
+  stats->error = dS;
   return TRUE;
 }
-
 
 // These should probably be tunable:
@@ -224,6 +238,4 @@
 
 int fit_least_squares (double *fit, double *err, double *y, double *dy, double *wgt, double *wt, int Npts);
-double VectorFractionInterpolate (double *values, float fraction, int Npts);
-double weight_cauchy (double x);
 
 // this is a zero-order fit (constant value only)
@@ -433,2 +445,3 @@
 // \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) 
 // R = \sum (weight_i^2)
+
Index: trunk/Ohana/src/relphot/src/setMrelCatalog.c
===================================================================
--- trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 41451)
+++ trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 41453)
@@ -65,5 +65,5 @@
   liststats_setmode (&results->kronstats, STATMODE);
 
-  SetMrelInfoReset (results); // reset the countesrs
+  SetMrelInfoReset (results); // reset the counters
 
   for (j = 0; j < catalog[Nc].Naverage; j++) {
@@ -94,5 +94,5 @@
 // eg, -refcode g_SDSS
 // this probably makes no sense in the context of multifilter analysis
-// XXX probably need to use the photocde table to assign reference mag weights.
+// XXX probably need to use the photcode table to assign reference mag weights.
 
 // dlist gives the error per measurement, wlist gives the weight
Index: trunk/Ohana/src/relphot/src/test_fitirls.c
===================================================================
--- trunk/Ohana/src/relphot/src/test_fitirls.c	(revision 41453)
+++ trunk/Ohana/src/relphot/src/test_fitirls.c	(revision 41453)
@@ -0,0 +1,175 @@
+# include "relphot.h"
+
+int mkdata (FitDataSet *dataset, int fitSlope, double valueTru, double slopeTru, int Npoints, int Nbad, int GaussianOutliers);
+
+int main (int argc, char **argv) {
+  
+  int Ntests     =  100;
+  int Npoints    =  100;
+  int Noutliers  =   10;
+  int Nbootstrap =  100;
+
+  int N;
+  if ((N = get_argument (argc, argv, "-Ntests"))) {
+    remove_argument (N, &argc, argv);
+    Ntests = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Npoints"))) {
+    remove_argument (N, &argc, argv);
+    Npoints = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Noutliers"))) {
+    remove_argument (N, &argc, argv);
+    Noutliers = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Nbootstrap"))) {
+    remove_argument (N, &argc, argv);
+    Nbootstrap = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  int GaussianOutliers = TRUE;
+  if ((N = get_argument (argc, argv, "-flat-outliers"))) {
+    remove_argument (N, &argc, argv);
+    GaussianOutliers = FALSE;
+  }
+
+  int FitSlope = FALSE;
+  if ((N = get_argument (argc, argv, "-fit-slope"))) {
+    remove_argument (N, &argc, argv);
+    FitSlope = TRUE;
+  }
+  
+  double ValuePrior = NAN;
+  double ValueSigma = NAN;
+  if ((N = get_argument (argc, argv, "-value-prior"))) {
+    remove_argument (N, &argc, argv);
+    ValuePrior = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ValueSigma = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  double SlopePrior = NAN;
+  double SlopeSigma = NAN;
+  if ((N = get_argument (argc, argv, "-slope-prior"))) {
+    if (!FitSlope) { fprintf (stderr, "-slope-prior only valid if -fit-slope is selected\n"); exit (2); }
+    remove_argument (N, &argc, argv);
+    SlopePrior = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    SlopeSigma = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 1) {
+    fprintf (stderr, "USAGE: %s\n", argv[0]);
+    fprintf (stderr, " some options: -Ntests [3000], -Npoints [100], -Noutliers [10], -flat-outliers\n");
+    fprintf (stderr, " more options: -fit-slope, -value-prior (value) (sigma), -slope-prior (slope) (sigma)\n");
+    exit (2);
+  }
+
+  // plan_tests (14);
+  // diag ("relastro fitpm tests");
+  
+  // init the random seed
+  { 
+    struct timeval now;
+    gettimeofday (&now, NULL);
+    long A = now.tv_sec + now.tv_usec * 1000000;
+    srand48(A);
+    // srand48(1);
+  }
+
+  ohana_gaussdev_init ();
+
+  FitDataSet dataset;
+  FitDataSetAlloc (&dataset, Npoints + Noutliers, FitSlope, Nbootstrap);
+
+  if (isfinite(ValuePrior) || isfinite(SlopePrior)) {
+    FitDataSetAddPriors (&dataset);
+    if (isfinite(ValuePrior)) {
+      dataset.bPriorValue[0] = ValuePrior;
+      dataset.bPriorSigma[0] = ValueSigma;
+    }
+    if (isfinite(SlopePrior)) {
+      dataset.bPriorValue[1] = SlopePrior;
+      dataset.bPriorSigma[1] = SlopeSigma;
+    }
+  }
+
+  for (int i = 0; i < Ntests; i++) {
+
+    double valueTru = 25.0;
+    double slopeTru = 0.05;
+
+    // generate a single fake star with Npoints real points and Noutliers bad points
+    mkdata (&dataset, FitSlope, valueTru, slopeTru, Npoints, Noutliers, GaussianOutliers);
+
+    fit1d_irls (&dataset, Npoints + Noutliers);
+
+    if (FitSlope) {
+      fprintf (stdout, "%8.3f %8.5f : %6.3f %6.3f : %8.5f %8.3f %6.3f : %6.3f %6.3f : %6.3f : %4d %4d %4d\n", valueTru, dataset.bSaveArray[0][0], dataset.bSigma[0], slopeTru, dataset.bSaveArray[1][0], dataset.bSigma[1], dataset.sigma, dataset.min, dataset.max, dataset.chisq, Npoints, Noutliers, dataset.Nmeas);
+    } else {
+      fprintf (stdout, "%8.3f %8.5f : %6.3f %6.3f : %6.3f %6.3f : %6.3f : %4d %4d %4d\n", valueTru, dataset.bSaveArray[0][0], dataset.bSigma[0], dataset.sigma, dataset.min, dataset.max, dataset.chisq, Npoints, Noutliers, dataset.Nmeas);
+    }
+  }
+
+  ohana_gaussdev_free();
+
+  FitDataSetFree (&dataset);
+
+  ohana_memcheck (TRUE);
+  ohana_memdump (TRUE);
+
+  exit (0);
+}
+
+// we are using IRLS 1D fitting to solve for zero points and airmass slope
+// I will generate data samples with a nominal zero point and an airmass slope
+
+// valueObs = valueTru + slopeTru*trend
+
+// choose a random value for 'trend' between MIN_TREND and MAX_TREND
+// set the raw value to valueTru + slopeTru*trend
+// vary by Gaussian deviate using sigma of VAL_ERROR
+
+# define VAL_ERROR 0.03
+# define MIN_TREND 0.0
+# define MAX_TREND 1.0
+int mkdata (FitDataSet *dataset, int FitSlope, double valueTru, double slopeTru, int Npoints, int Nbad, int GaussianOutliers) {
+  
+  for (int i = 0; i < Npoints; i++) {
+    double dF = ohana_gaussdev_rnd(0.0, VAL_ERROR);
+
+    double trend = MIN_TREND + (MAX_TREND - MIN_TREND)*drand48();
+
+    double valueRaw = FitSlope ? valueTru + slopeTru*trend : valueTru;
+    double valueObs = valueRaw + dF;
+
+    dataset->alldata->yVector[i] = valueObs;
+    dataset->alldata->dyVector[i] = VAL_ERROR;
+    if (FitSlope) { 
+      dataset->alldata->xVector[i] = trend;
+    }
+  }
+
+  for (int i = 0; i < Nbad; i++) {
+    double dF = GaussianOutliers ? ohana_gaussdev_rnd(10.0*VAL_ERROR, 10.0*VAL_ERROR) : 10.0*VAL_ERROR*drand48();
+
+    double trend = MIN_TREND + (MAX_TREND - MIN_TREND)*drand48();
+
+    double valueRaw = FitSlope ? valueTru + slopeTru*trend : valueTru;
+    double valueObs = valueRaw + dF;
+
+    dataset->alldata->yVector[i + Npoints] = valueObs;
+    dataset->alldata->dyVector[i + Npoints] = VAL_ERROR;
+    if (FitSlope) { 
+      dataset->alldata->xVector[i + Npoints] = trend;
+    }
+  }
+
+  return TRUE;
+}
Index: trunk/Ohana/src/relphot/test/relphot.tgroups.dvo
===================================================================
--- trunk/Ohana/src/relphot/test/relphot.tgroups.dvo	(revision 41451)
+++ trunk/Ohana/src/relphot/test/relphot.tgroups.dvo	(revision 41453)
@@ -52,11 +52,17 @@
   # generate the basic images and check they were correctly ingested by dvo
   for i 0 mjd_uc[]
-    sprintf filerootout "%s.%02d" $fileroot $i
+    sprintf filerootout "%s.uc.%02d" $fileroot $i
     mkexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] secz_uc[$i] $FILTER $filerootout $RA_CENTER $DEC_CENTER
     ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] secz_uc[$i] $FILTER raw
   end
 
+  for i 0 mjd_nc[]
+    sprintf filerootout "%s.nc.%02d" $fileroot $i
+    mkexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] secz_nc[$i] $FILTER $filerootout $RA_CENTER $DEC_CENTER
+    ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] secz_nc[$i] $FILTER raw
+  end
+
   break
-  
+
   # run relphot on the db and check that the images now match the expected values
   tapEXEC relphot -tgroup-fit-airmass -tgroups tgroups.dat -images g -v -region $RA_MIN $RA_MAX $DEC_MIN $DEC_MAX -D CATDIR $catdir -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12