Changeset 41387

Timestamp:

Jul 7, 2020, 9:30:05 AM (6 years ago)

Author:

eugene

Message:

updates to get tgroups working

Location:

branches/eam_branches/relphot.20200519

Files:

• include/relphot.h (modified) (8 diffs)
• src/GridOps.c (modified) (10 diffs)
• src/ImageOps.c (modified) (8 diffs)
• src/MosaicOps.c (modified) (5 diffs)
• src/StarOps.c (modified) (4 diffs)
• src/TGroupOps.c (modified) (22 diffs)
• src/args.c (modified) (3 diffs)
• src/liststats.c (modified) (2 diffs)
• src/reload_catalogs.c (modified) (3 diffs)
• src/relphot_images.c (modified) (3 diffs)
• src/setMrelCatalog.c (modified) (6 diffs)

branches/eam_branches/relphot.20200519/include/relphot.h

@@ -96 +96 @@
   float McalPSF;
   float McalAPER;
+  float dKlam;
   float dMcal;
   float dMsys;
@@ -110 +111 @@
   off_t *measure;
 
-  off_t Ncatalog;
-  off_t NCATALOG;
+  off_t Ncatalog; // XX not needed
+  off_t NCATALOG; // XX not needed
   off_t *catalog;
 } TGroup; 
@@ -332 +333 @@
 int    MOSAIC_ZEROPT;
 int    TGROUP_ZEROPT;
+int    TGROUP_FIT_AIRMASS;
 
 int    TEST_IMAGE1;
@@ -457 +459 @@
 void          freeImageBins       PROTO((int Ncatalog, int doImageList));
 void          freeMosaicBins      PROTO((int Ncatalog, int doMosaicList));
-void          freeTGroupBins      PROTO((int Ncatalog, int doTGroupList));
+void          freeTGroupBins      PROTO((int Ncatalog));
 
 void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
@@ -468 +470 @@
 float         getMgrid            PROTO((off_t meas, int cat));
 float         getMmos             PROTO((off_t meas, int cat));
-float         getMgrp             PROTO((off_t meas, int cat));
+float         getMgrp             PROTO((off_t meas, int cat, float airmass));
 float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source));
 short         getUbercalDist      PROTO((off_t meas, int cat));
@@ -485 +487 @@
 void          initMosaics         PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage));
 void          initTGroups         PROTO((Image *subset, off_t *LineNumber, off_t Nsubset));
+void          freeTGroups         PROTO((void));
 void          initTGroupBins      PROTO((Catalog *catalog, int Ncatalog));
 void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
@@ -493 +496 @@
 int           liststats_init      PROTO((StatType *stats));
 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));
 
 unsigned int *ReadTGroupFile      PROTO((FILE *f, int *nelem));
@@ -721 +725 @@
 void RepairWarpMeasuresOpenLogfile ();
 void RepairWarpMeasuresCloseLogfile ();
+
+int save_images_backup (FITS_DB *db);

branches/eam_branches/relphot.20200519/src/GridOps.c

@@ -289 +289 @@
 
   int **gotstar, **gridmeas;
-  int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys;
+  int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys, Ngrp;
   double **A, **B, *Mjx, *Wjx;
-  float Msys, Mcal, Mmos, Merr, Wsys;
+  float Msys, Mcal, Mmos, Mgrp, Merr, Wsys;
   double Mj, Wj;
   
@@ -305 +305 @@
   }
 
-  Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = 0;
+  Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = Ngrp = 0;
 
   ALLOCATE (gotstar, int *, Ncatalog);
@@ -411 +411 @@
         if (isnan(Mmos)) {
           Nmos ++;
+          continue;
+        }
+
+        // skip mosaics marked as BAD
+        Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+        if (isnan(Mgrp)) {
+          Ngrp ++;
           continue;
         }
@@ -526 +533 @@
 void setMgrid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
 
-  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys, Nfit;
+  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ngrp, Ncal, Nrel, Nsys, Nfit;
   double *list, *dlist;
-  float Msys, Mrel, Mcal, Mmos;
+  float Msys, Mrel, Mcal, Mmos, Mgrp;
 
   StatType stats;
@@ -544 +551 @@
   ALLOCATE (dlist, double, Nmax);
 
-  Nbad = Ncal = Nmos = Nrel = Nsys = Nfit = 0;
+  Nbad = Ncal = Nmos = Ngrp = Nrel = Nsys = Nfit = 0;
 
   for (i = 0; i < Ngrid; i++) {
@@ -568 +575 @@
         continue;
       }
+      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      if (isnan(Mgrp)) {
+        Ngrp ++;
+        continue;
+      }
       Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
       if (isnan(Mrel)) {
@@ -580 +592 @@
         continue;
       }
-      list[N] = Msys - Mrel - Mcal - Mmos;
+      list[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
       dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
       N++;
@@ -616 +628 @@
 
   int i, j, m, c, n, N, Narea;
-  float Msys, Mrel, Mcal, Mmos;
+  float Msys, Mrel, Mcal, Mmos, Mgrp;
   double *xlist, *Mlist, *dlist, *ylist;
   Graphdata graphdata;
@@ -649 +661 @@
       Mmos  = getMmos  (m, c);
       if (isnan(Mmos)) continue;
+      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      if (isnan(Mgrp)) continue;
       Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
       if (isnan(Mrel)) continue;
@@ -657 +671 @@
       xlist[N] = Xmeas[c][m];
       ylist[N] = Ymeas[c][m];
-      Mlist[N] = Msys - Mrel - Mcal - Mmos;
-      dlist[N] = Msys - Mrel - Mcal - Mmos - gridM[i];
+      Mlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
+      dlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp - gridM[i];
       N++;
     }

branches/eam_branches/relphot.20200519/src/ImageOps.c

@@ -586 +586 @@
 
   off_t i, j, m, c, n;
-  int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys;
+  int mark, bad, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys;
 
   StatType stats;
@@ -612 +612 @@
   StatDataSet *brightStars = StatDataSetAlloc (1, Nmax);
 
-  Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;
+  Nfew = Nbad = Nmos = Ngrp = Ngrid = Nrel = Nsys = 0;
 
   int Ncalibrated = 0;
@@ -631 +631 @@
 
     // only freeze the gpc1 chips
-    if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
+    // XXX : when do we FREEZE_IMAGES in a mix of GPC1, GPC2, etc?
+    // if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
+    if (FREEZE_IMAGES) continue;
     if (image[i].photcode == 0) continue; // skip the PHU images
 
@@ -666 +668 @@
           continue;
       }
+      float Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      if (isnan(Mgrp)) {
+          Ngrp ++;
+          continue;
+      }
       float Mgrid = getMgrid (m, c);
       if (isnan(Mgrid)) {
@@ -706 +713 @@
 
     skip:
-      psfStars->flxlist[Nref] = MsysPSF - MrelPSF - Mmos - Mgrid + Mflat;
+      psfStars->flxlist[Nref] = MsysPSF - MrelPSF - Mmos - Mgrp - Mgrid + Mflat;
       psfStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
       psfStars->wgtlist[Nref] = 1;
@@ -718 +725 @@
 
       if (isfinite(MrelKron) && isfinite(MsysKron)) {
-        kronStars->flxlist[Nkron] = MsysKron - MrelKron - Mmos - Mgrid + Mflat;
+        kronStars->flxlist[Nkron] = MsysKron - MrelKron - Mmos - Mgrp - Mgrid + Mflat;
         kronStars->errlist[Nkron] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
         kronStars->wgtlist[Nkron] = 1;
@@ -1023 +1030 @@
   off_t i, j, m, c, n, N;
   double *list, *dlist;
-  float Mcal, Mmos, Mgrid;
+  float Mcal, Mmos, Mgrp, Mgrid;
 
   StatType stats;
@@ -1050 +1057 @@
       Mmos  = getMmos  (m, c);
       if (isnan(Mmos)) continue;
+      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      if (isnan(Mgrp)) continue;
       Mgrid = getMgrid (m, c);
       if (isnan(Mgrid)) continue;

branches/eam_branches/relphot.20200519/src/MosaicOps.c

@@ -919 +919 @@
   int Nbad;
   int Ncal;
+  int Ngrp;
   int Ngrid;
   int Nrel;
@@ -1154 +1155 @@
     if (fout) {
       float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
+      // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass);
       float Mgrid    = getMgrid (m, c);
       float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
@@ -1175 +1177 @@
     if (isnan(Mcal)) {
       info->Ncal++;
+      continue;
+    }
+    float Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+    if (isnan(Mgrp)) {
+      info->Ngrp++;
       continue;
     }
@@ -1692 +1699 @@
   off_t i, j, m, c, n, N;
   double *list, *dlist;
-  float Mcal, Mgrid, Mrel;
+  float Mcal, Mgrp, Mgrid, Mrel;
   StatType stats;
 
@@ -1716 +1723 @@
       Mcal = getMcal  (m, c, MAG_CLASS_PSF);
       if (isnan(Mcal)) continue;
+      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
+      if (isnan(Mgrp)) continue;
       Mgrid = getMgrid (m, c);
       if (isnan(Mgrid)) continue;

branches/eam_branches/relphot.20200519/src/StarOps.c

@@ -491 +491 @@
         float Mmos  = getMmos  (m, i);
         if (isnan(Mmos)) continue;
-        float Mgrp  = getMgrp  (m, i);
+        float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass);
         if (isnan(Mgrp)) continue;
         float Mgrid = getMgrid (m, i);
@@ -507 +507 @@
 
         // set the output calibration
-        float Moff = Mmos + Mgrp + Mgrid;
+        // float Moff = Mmos + Mgrp + Mgrid;
+        float Moff = Mmos + Mgrid; // Mgrp has already been applied to the image.Mcal
         catalog[i].measure[m].McalPSF  = McalPSF + Moff;
         catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Moff : McalPSF + Moff;
@@ -730 +731 @@
           Mmos  = getMmos  (m, i);
           if (isnan(Mmos)) { Nmos ++; continue; }
-          Mgrp  = getMgrp  (m, i);
+          Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass);
           if (isnan(Mgrp)) { Ngrp ++; continue; }
           Mgrid = getMgrid (m, i);
@@ -736 +737 @@
 
           Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
-          list[N] = Msys - Mcal - Mmos - Mgrid;
+          list[N] = Msys - Mcal - Mmos - Mgrp - Mgrid;
           dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
           N++;

branches/eam_branches/relphot.20200519/src/TGroupOps.c

@@ -36 +36 @@
 // 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.
 void loadTGroups (char *filename) {
 
@@ -73 +74 @@
       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;
@@ -111 +113 @@
 
     // ignore non-GPC1 images
-    if (!isGPC1chip(subset[i].photcode)) {
-      Nsimple ++;
+    // if (!isGPC1chip(subset[i].photcode)) {
+    //   Nsimple ++;
+    //   continue;
+    // }
+
+    // ignore mosaic images (photcode == 0)
+    if (!subset[i].photcode) {
       continue;
     }
@@ -138 +145 @@
 
   fprintf (stderr, "matched %d images to %d tgroups, %d simple chips not matched to tgroups\n", (int) (Nsubset - Nsimple), (int) Ntgroup, (int) Nsimple);
+  return;
+}
+
+void freeTGroups (void) {
+
+  if (!TGROUP_ZEROPT) return;
+
+  free (ImageToTGroup);
+  free (tgroupTimes);
+
+  for (int i = 0; i < Ntgroup; i++) {
+    free (tgroup[i].image);
+    free (tgroup[i].measure);
+    free (tgroup[i].catalog);
+  }
+  free (tgroup);
   return;
 }
@@ -221 +244 @@
   // if start is not in this range, return -1
   if (start < tgroupTimes[Nlo]) return (-1);
-  if (start > tgroupTimes[Nhi]) return (-1);
+  if (start > tgroupTimes[Nhi] + 86399) return (-1); // 1-day range for tgroup for now
 
   while (Nhi - Nlo > 4) {
@@ -259 +282 @@
     for (j = 0; j < tgroup[i].Nimage; j++) {
       im = tgroup[i].image[j];
-      image[im].McalPSF    += tgroup[i].McalPSF;
-      image[im].McalAPER   += tgroup[i].McalAPER;
+      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;
@@ -267 +291 @@
       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;
   }      
 }
@@ -287 +314 @@
 }
 
-void freeTGroupBins (int Ncatalog, int doTGroupList) {
+void freeTGroupBins (int Ncatalog) {
 
   off_t i;
@@ -329 +356 @@
     }
   }
-  // fprintf (stderr, "matched %d detections to tgroups\n", Nmatch);
+  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);
+  }
   return (TRUE);
 }
@@ -380 +410 @@
 }
 
-float getMgrp (off_t meas, int cat) {
+float getMgrp (off_t meas, int cat, float airmass) {
 
   off_t i;
@@ -392 +422 @@
 
   // if (tgroup[i].flags & IMAGE_BAD) return (NAN);  
-  value = tgroup[i].McalPSF;
+  value = tgroup[i].McalPSF + tgroup[i].dKlam*(airmass - 1.0);
   return (value);
 }
@@ -419 +449 @@
   int PoorImages;
   double *psfMagList;
+  double *psfAirList;
   double *psfErrList;
   double *brightMagList;
@@ -456 +487 @@
   if (allocLists) {
     ALLOCATE (info->psfMagList, double, Nmax);
+    ALLOCATE (info->psfAirList, double, Nmax);
     ALLOCATE (info->psfErrList, double, Nmax);
     ALLOCATE (info->kronMagList, double, Nmax);
@@ -466 +498 @@
 void SetMgrpInfoFree (SetMgrpInfo *info) {
   free (info->psfMagList   );
+  free (info->psfAirList   );
   free (info->psfErrList   );
   free (info->kronMagList  );
@@ -522 +555 @@
   image = getimages (&N, NULL);
 
-  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
-
   if (PoorImages) {
     // XXX use bad stars and measurements for PoorImages? or not?
@@ -543 +574 @@
   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);
   }
   SetMgrpInfoFree (&info);
@@ -575 +607 @@
 
   double *psfMagList    = info->psfMagList;
+  double *psfAirList    = info->psfAirList;
   double *psfErrList    = info->psfErrList;
   double *kronMagList   = info->kronMagList;
@@ -593 +626 @@
 
   int testImage = FALSE;
-  // testImage |= (abs(myTGroup[0].start - 1324104046) < 10);
+  testImage |= (abs(myTGroup[0].start - 1245283200) < 10);
   // testImage |= (abs(myTGroup[0].start - 1324103823) < 10);
 
@@ -633 +666 @@
       int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
 
-      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad);
+      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, catalog[c].measureT[m].airmass, delta, isBad);
     }
 
@@ -697 +730 @@
 
     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;
@@ -725 +759 @@
   }
 
-  liststats (psfMagList, psfErrList, NULL, N, &stats);
+  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);
 
@@ -733 +772 @@
   myTGroup[0].McalChiSq  = stats.chisq;
   myTGroup[0].nFitPhotom = N;
+  myTGroup[0].dKlam      = dk;
 
   // drop this for tgroups?

branches/eam_branches/relphot.20200519/src/args.c

@@ -317 +317 @@
   }
 
+  FREEZE_IMAGES = FALSE;
+  if ((N = get_argument (argc, argv, "-imfreeze"))) {
+    remove_argument (N, &argc, argv);
+    FREEZE_IMAGES = TRUE;
+  }
+  FREEZE_MOSAICS = FALSE;
+  if ((N = get_argument (argc, argv, "-mosfreeze"))) {
+    remove_argument (N, &argc, argv);
+    FREEZE_MOSAICS = TRUE;
+  }
+  CALIBRATE_STACKS_AND_WARPS = FALSE;
+  if ((N = get_argument (argc, argv, "-only-stacks-and-warps"))) {
+    remove_argument (N, &argc, argv);
+    CALIBRATE_STACKS_AND_WARPS = TRUE;
+  }
+  USE_MCAL_PSF_FOR_STACK_APER = FALSE;
+  if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) {
+    remove_argument (N, &argc, argv);
+    USE_MCAL_PSF_FOR_STACK_APER = TRUE;
+  }
+
   // calibrate the photometry by night / time period
   TGROUP_ZEROPT = FALSE;
+  char *TGROUP_FILENAME = NULL;
   if ((N = get_argument (argc, argv, "-tgroup"))) {
     TGROUP_ZEROPT = TRUE;
     remove_argument (N, &argc, argv);
-    loadTGroups (argv[N]);
-    remove_argument (N, &argc, argv);
-  }
-
-  FREEZE_IMAGES = FALSE;
-  if ((N = get_argument (argc, argv, "-imfreeze"))) {
+    TGROUP_FILENAME = strcreate(argv[N]);
     remove_argument (N, &argc, argv);
     FREEZE_IMAGES = TRUE;
   }
-  FREEZE_MOSAICS = FALSE;
-  if ((N = get_argument (argc, argv, "-mosfreeze"))) {
-    remove_argument (N, &argc, argv);
-    FREEZE_MOSAICS = TRUE;
-  }
-  CALIBRATE_STACKS_AND_WARPS = FALSE;
-  if ((N = get_argument (argc, argv, "-only-stacks-and-warps"))) {
-    remove_argument (N, &argc, argv);
-    CALIBRATE_STACKS_AND_WARPS = TRUE;
-  }
-  USE_MCAL_PSF_FOR_STACK_APER = FALSE;
-  if ((N = get_argument (argc, argv, "-use-mcal-psf-for-stack-aper"))) {
-    remove_argument (N, &argc, argv);
-    USE_MCAL_PSF_FOR_STACK_APER = TRUE;
+  if ((N = get_argument (argc, argv, "-tgroups"))) {
+    if (TGROUP_ZEROPT) {
+      fprintf (stderr, "use only one of -tgroups and -tgroup (same meaning)\n");
+      exit (2);
+    }
+    TGROUP_ZEROPT = TRUE;
+    remove_argument (N, &argc, argv);
+    TGROUP_FILENAME = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+    FREEZE_IMAGES = TRUE;
+  }
+  TGROUP_FIT_AIRMASS = FALSE;
+  if ((N = get_argument (argc, argv, "-tgroup-fit-airmass"))) {
+    if (!TGROUP_ZEROPT) {
+      fprintf (stderr, "-tgroup-fit-airmass requires -tgroup (filename) option\n");
+      exit (2);
+    }
+    TGROUP_FIT_AIRMASS = TRUE;
+    remove_argument (N, &argc, argv);
   }
 
@@ -528 +550 @@
   if (argc != 1) relphot_usage ();
 
+  if (TGROUP_ZEROPT) {
+    loadTGroups (TGROUP_FILENAME);
+    free (TGROUP_FILENAME);
+  }
+
   return mode;
 }
@@ -551 +578 @@
   FreePhotcodeTable();
   
+  freeTGroups();
+
   ohana_memcheck (VERBOSE);
   ohana_memdump (VERBOSE);

branches/eam_branches/relphot.20200519/src/liststats.c

@@ -172 +172 @@
   for (i = ks; i < ke; i++) {
     M  = SQ (value[i] - Mo);
-    dM = SQ (dvalue[i]);
-    X2 += M / dM;
+    if (dvalue) {
+      dM = SQ (dvalue[i]);
+      X2 += M / dM;
+    }
     dS += M;
   }
-  X2 = X2 / (Nm - 1);
+  X2 = dvalue ? X2 / (Nm - 1) : NAN;
   dS = sqrt (dS / (Nm - 1));
 
@@ -187 +189 @@
   return (TRUE);
 }
+
+int liststats_fit1d (double *value, double *err, double *x, int Npts, StatType *stats, double *dk) {
+
+  double M = 0.0, MZ = 0.0, Z1 = 0.0, Z2 = 0.0, R = 0.0;
+
+  for (int i = 0; i < Npts; i++) {
+    double wt = 1.0 / (err[i]*err[i]);
+    double V = value[i];
+    double Z = x[i];
+    M  += V   * wt;
+    MZ += V*Z * wt;
+    Z1 += Z   * wt;
+    Z2 += Z*Z * wt;
+    R  += 1.0 * wt;
+  }
+
+  double DetInv = (R*Z2 - Z1*Z1);
+  double Det = (fabs(DetInv) < -1e6) ? 0.0 : 1.0 / DetInv; // do not allow unstable solutions
+
+  double zp = (M*Z2 - Z1*MZ)*Det;
+  *dk = (MZ*R - M*Z1)*Det;
+
+  stats->mean = zp;
+  stats->chisq = NAN;
+  return TRUE;
+}
+
 
 // These should probably be tunable:

branches/eam_branches/relphot.20200519/src/reload_catalogs.c

@@ -86 +86 @@
     initImageBins  (&catalog, 1, FALSE);
     initMosaicBins (&catalog, 1, FALSE);
+    initTGroupBins (&catalog, 1);
+
     initGridBins   (&catalog, 1); 
     TIMESTAMP(time3b);
@@ -91 +93 @@
     findImages (&catalog, 1, FALSE);
     findMosaics (&catalog, 1, FALSE);
+    findTGroups (&catalog, 1);
+
     TIMESTAMP(time4);
 
@@ -117 +121 @@
     TIMESTAMP(time6);
 
+    // XXX freeTgroupBins here: 
+
     freeImageBins (1, FALSE);
     freeMosaicBins (1, FALSE);
+    freeTGroupBins (1);
     freeGridBins (1);
     TIMESTAMP(time7);

branches/eam_branches/relphot.20200519/src/relphot_images.c

@@ -81 +81 @@
     MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
 
-    findTGroups (catalog, Ncatalog);  /* also sets Grid values */
+    findTGroups (catalog, Ncatalog);
     MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
 
@@ -158 +158 @@
       if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog);
       if ((i > 8) && (i % 8 == 5)) clean_mosaics ();
+      if ((i > 8) && (i % 8 == 5)) clean_tgroups ();
       if ((i > 8) && (i % 8 == 5)) clean_images ();
 
@@ -258 +259 @@
     freeImageBins (Ncatalog, TRUE);
     freeMosaicBins (Ncatalog, TRUE);
+    freeTGroupBins (Ncatalog);
     freeGridBins (Ncatalog);
 

branches/eam_branches/relphot.20200519/src/setMrelCatalog.c

@@ -221 +221 @@
 
     // ** Choose the calibration (depends on the mode : do I have an image reference or not?) 
-    float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0;
+    float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0, Mgrp = 0.0;
     off_t meas = measureOffset + k;
     if (getImageEntry (meas, cat) < 0) {
@@ -235 +235 @@
       Mmos  = getMmos  (meas, cat);
       if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
+      Mgrp  = getMgrp  (meas, cat, measureT[k].airmass);
+      if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
       Mgrid = getMgrid (meas, cat);
       if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
@@ -261 +263 @@
       if (CHECK_VALID_MAG(Map, dMap)) {
         int Nap = results->aperData[Nsec].Nlist;
-        results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
+        results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
         results->aperData[Nsec].errlist[Nap] = dMap;
         results->aperData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0;
@@ -281 +283 @@
     if (CHECK_VALID_MAG(Mkron, dMkron)) {
       int Nkron = results->kronData[Nsec].Nlist;
-      results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
+      results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
       results->kronData[Nsec].errlist[Nkron] = dMkron;
       results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
@@ -294 +296 @@
     if (CHECK_VALID_MAG(Mpsf, dMpsf)) {
       int Npsf = results->psfData[Nsec].Nlist;
-      results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
+      results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
       results->psfData[Nsec].errlist[Npsf] = dMpsf;
       results->psfData[Nsec].wgtlist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
@@ -711 +713 @@
       Mmos      = getMmos  (measSeq, cat);
       Mgrid     = getMgrid (measSeq, cat);
+      // XXX can Mmos and Mgrid exist for stacks?
     }