Changeset 41453 for trunk

Timestamp:

Dec 7, 2020, 12:27:41 PM (6 years ago)

Author:

eugene

Message:

working on adapting IRLS fitting with Priors for zero points

Location:

trunk/Ohana/src/relphot

Files:

• Makefile (modified) (4 diffs)
• include/relphot.h (modified) (7 diffs)
• src/ConfigInit.c (modified) (4 diffs)
• src/ImageOps.c (modified) (10 diffs)
• src/MosaicOps.c (modified) (2 diffs)
• src/StarOps.c (modified) (1 diff)
• src/TGroupOps.c (modified) (43 diffs)
• src/args.c (modified) (1 diff)
• src/fit1d_irls.c (added)
• src/initialize.c (modified) (1 diff)
• src/liststats.c (modified) (3 diffs)
• src/setMrelCatalog.c (modified) (2 diffs)
• src/test_fitirls.c (added)
• test/relphot.tgroups.dvo (modified) (1 diff)

trunk/Ohana/src/relphot/Makefile

@@ -21 +21 @@
 
 test_liststats: $(BIN)/test_liststats.$(ARCH)
+
+test_fitirls: $(BIN)/test_fitirls.$(ARCH)
 
 install: $(DESTBIN)/relphot $(DESTBIN)/relphot_client
@@ -48 +50 @@
 $(SRC)/global_stats.$(ARCH).o    \
 $(SRC)/initialize.$(ARCH).o      \
+$(SRC)/fit1d_irls.$(ARCH).o      \
 $(SRC)/liststats.$(ARCH).o       \
 $(SRC)/load_catalogs.$(ARCH).o   \
@@ -107 +110 @@
 $(SRC)/synthetic_mags.$(ARCH).o  \
 $(SRC)/plotstuff.$(ARCH).o       \
+$(SRC)/fit1d_irls.$(ARCH).o      \
 $(SRC)/liststats.$(ARCH).o       \
 $(SRC)/initialize.$(ARCH).o      \
@@ -137 +141 @@
 $(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)
+

trunk/Ohana/src/relphot/include/relphot.h

@@ -91 +91 @@
 
 typedef struct {
-  unsigned int start;
-  unsigned int stop;
   short photcode;
   float McalPSF;
@@ -110 +108 @@
   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 +168 @@
 
 } 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 +365 @@
 double IMAGE_SCATTER;
 double IMAGE_OFFSET;
+double NIGHT_SCATTER;
+double NIGHT_OFFSET;
 double STAR_SCATTER;
 double STAR_CHISQ;
@@ -497 +553 @@
 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 +566 @@
 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 +591 @@
 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));

trunk/Ohana/src/relphot/src/ConfigInit.c

@@ -1 +1 @@
 # 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 +26 @@
   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 +53 @@
   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 +98 @@
   return;
 }
+

trunk/Ohana/src/relphot/src/ImageOps.c

@@ -1 +1 @@
 # 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 +586 @@
   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 +603 @@
   }
 
-  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 +710 @@
 
     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 +717 @@
 
       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 +750 @@
 
     // 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 +766 @@
     }
 
-    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 +783 @@
     image[i].ubercalDist = minUbercalDist + 1;
   }
-
-  StatDataSetFree (brightStars, 1);
-  StatDataSetFree (kronStars, 1);
-  StatDataSetFree (psfStars, 1);
 
   fprintf (stderr, "%d images calibrated\n", Ncalibrated);
@@ -825 +792 @@
     MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
   }
+
+  FitDataSetFree (&brightStars);
+  FitDataSetFree (&kronStars);
+  FitDataSetFree (&psfStars);
+
   return;
 }
@@ -891 +863 @@
 
 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;

trunk/Ohana/src/relphot/src/MosaicOps.c

@@ -1 +1 @@
 # 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 +1280 @@
   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);
   }
 

trunk/Ohana/src/relphot/src/StarOps.c

@@ -228 +228 @@
 
 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);

trunk/Ohana/src/relphot/src/TGroupOps.c

@@ -1 +1 @@
 # 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 +41 @@
     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 +62 @@
 
   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 +158 @@
    */
 
-  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 +176 @@
     }
 
-    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 +191 @@
 
     // 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 +201 @@
   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 +210 @@
 
   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 +235 @@
   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 +283 @@
   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 +291 @@
 // 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 +307 @@
   // 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 +348 @@
   // 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 +380 @@
   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 +430 @@
   }
   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 +450 @@
   }
 
-  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 +465 @@
     // 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 +474 @@
   // 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 +490 @@
 float getMgrp (off_t meas, int cat, float airmass) {
 
-  off_t i;
   float value;
 
@@ -421 +495 @@
 
   // 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 +508 @@
 
   // 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 +602 @@
 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 +620 @@
 int setMgrp (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
 
-  off_t i, N, Nmax;
+  off_t N;
   Image *image;
 
@@ -564 +640 @@
   }
 
-  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 +654 @@
   // 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 +683 @@
 # 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 +820 @@
     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 +842 @@
     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 +852 @@
       }
     } 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 +922 @@
 
   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 +984 @@
              threadinfo[i].info.Nrel, 
              threadinfo[i].info.Nsys,
-             threadinfo[i].info.Nskip);
+             threadinfo[i].info.Nskip);
     SetMgrpInfoAccum (&summary, &threadinfo[i].info);
   }
@@ -924 +1027 @@
     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 +1038 @@
 }
 
+# if (0)
 StatType statsTGroupM (Catalog *catalog) {
   OHANA_UNUSED_PARAM(catalog);
@@ -943 +1050 @@
   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 +1179 @@
   return (stats);
 }
+# endif
 
 /* mark tgroup if: abs(Mcal - <Mcal>) too large, dMcal too large */
@@ -1080 +1192 @@
   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 +1211 @@
 
   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 +1255 @@
 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 +1265 @@
   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 +1276 @@
   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 +1310 @@
   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 +1330 @@
   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 +1343 @@
   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);

trunk/Ohana/src/relphot/src/args.c

@@ -579 +579 @@
   
   freeTGroups();
+  free_error();
 
   ohana_memcheck (VERBOSE);

trunk/Ohana/src/relphot/src/initialize.c

@@ -3 +3 @@
 RelphotMode initialize (int argc, char **argv) {
 
+  init_error();
   relphot_help (argc, argv);
   ConfigInit (&argc, argv);

trunk/Ohana/src/relphot/src/liststats.c

@@ -209 +209 @@
 
   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 +238 @@
 
 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 +445 @@
 // \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) 
 // R = \sum (weight_i^2)
+

trunk/Ohana/src/relphot/src/setMrelCatalog.c

@@ -65 +65 @@
   liststats_setmode (&results->kronstats, STATMODE);
 
-  SetMrelInfoReset (results); // reset the countesrs
+  SetMrelInfoReset (results); // reset the counters
 
   for (j = 0; j < catalog[Nc].Naverage; j++) {
@@ -94 +94 @@
 // 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

trunk/Ohana/src/relphot/test/relphot.tgroups.dvo

@@ -52 +52 @@
   # 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