Index: /branches/eam_branches/relphot.20210521/Makefile
===================================================================
--- /branches/eam_branches/relphot.20210521/Makefile	(revision 41603)
+++ /branches/eam_branches/relphot.20210521/Makefile	(revision 41604)
@@ -29,5 +29,5 @@
 $(SRC)/StatDataSetOps.$(ARCH).o	 \
 $(SRC)/ConfigInit.$(ARCH).o	 \
-$(SRC)/GridOpsNew.$(ARCH).o	 \
+$(SRC)/GridOps.$(ARCH).o	 \
 $(SRC)/ImageOps.$(ARCH).o	 \
 $(SRC)/ImageSubset.$(ARCH).o	 \
@@ -39,8 +39,4 @@
 $(SRC)/Shutdown.$(ARCH).o 	 \
 $(SRC)/StarOps.$(ARCH).o	 \
-$(SRC)/WarpImageMaps.$(ARCH).o	 \
-$(SRC)/StackImageMaps.$(ARCH).o	 \
-$(SRC)/RepairWarpMeasures.$(ARCH).o	 \
-$(SRC)/myIndex.$(ARCH).o	 \
 $(SRC)/psps_ids.$(ARCH).o	 \
 $(SRC)/args.$(ARCH).o		 \
@@ -92,5 +88,5 @@
 $(SRC)/StatDataSetOps.$(ARCH).o	 \
 $(SRC)/ConfigInit.$(ARCH).o	 \
-$(SRC)/GridOpsNew.$(ARCH).o	 \
+$(SRC)/GridOps.$(ARCH).o	 \
 $(SRC)/ImageOps.$(ARCH).o	 \
 $(SRC)/ImageSubset.$(ARCH).o	 \
@@ -102,8 +98,4 @@
 $(SRC)/Shutdown.$(ARCH).o 	 \
 $(SRC)/StarOps.$(ARCH).o	 \
-$(SRC)/WarpImageMaps.$(ARCH).o	 \
-$(SRC)/StackImageMaps.$(ARCH).o	 \
-$(SRC)/RepairWarpMeasures.$(ARCH).o	 \
-$(SRC)/myIndex.$(ARCH).o	 \
 $(SRC)/psps_ids.$(ARCH).o	 \
 $(SRC)/bcatalog.$(ARCH).o	 \
Index: /branches/eam_branches/relphot.20210521/include/relphot.h
===================================================================
--- /branches/eam_branches/relphot.20210521/include/relphot.h	(revision 41603)
+++ /branches/eam_branches/relphot.20210521/include/relphot.h	(revision 41604)
@@ -49,9 +49,9 @@
 
 typedef enum {
-    MODE_NONE   = 0,
-    MODE_LOAD   = 1,
-    MODE_UPDATE = 2,
+    MODE_NONE           = 0,
+    MODE_LOAD           = 1,
+    MODE_UPDATE         = 2,
     MODE_UPDATE_OBJECTS = 3,
-    MODE_SYNTH_PHOT = 4,
+    MODE_SYNTH_PHOT     = 4,
 } ModeType;
 
@@ -814,28 +814,9 @@
 void free_astrom_table ();
 
-int RepairWarpMeasures (Catalog *catalog);
-
-int FindWarpGroups (void);
-int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y);
-void FreeWarpGroups (void);
-
-myIndexType *myIndexInit ();
-int myIndexFree (myIndexType *myIndex);
-int myIndexUpdateLimits (myIndexType *myIndex, int value);
-int myIndexSetRange (myIndexType *myIndex);
-int myIndexSetEntry (myIndexType *myIndex, int value, int entry);
-int myIndexGetEntry (myIndexType *myIndex, int value);
-
 uint64_t CreatePSPSObjectID(double ra, double dec);
 uint64_t CreatePSPSStackDetectionID(int sourceID, int imageID, int detID);
 uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID);
 
-int MakeStackIndex (void);
-void FreeStackGroups (void);
-int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y);
 void sort_by_ra (double *R, double *D, int *I, int *S, int N);
-
-void RepairWarpMeasuresOpenLogfile ();
-void RepairWarpMeasuresCloseLogfile ();
 
 void dump_tgroups (Catalog *catalog, int Npass);
Index: /branches/eam_branches/relphot.20210521/src/ImageOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/ImageOps.c	(revision 41603)
+++ /branches/eam_branches/relphot.20210521/src/ImageOps.c	(revision 41604)
@@ -507,11 +507,4 @@
     offset = 0.0;
   }
-
-  // to do this, I need to pass in the catalog and flatcorr pointers
-  // int flat_id = image[i].photom_map_id;
-  // if (flat_id) {
-  //   offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd);
-  // }
-
   return (offset);
 }
Index: anches/eam_branches/relphot.20210521/src/RepairWarpMeasures.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/RepairWarpMeasures.c	(revision 41603)
+++ 	(revision )
@@ -1,181 +1,0 @@
-# include "relphot.h"
-
-FILE *logfile = NULL;
-
-void RepairWarpMeasuresOpenLogfile () {
-
-  char name[DVO_MAX_PATH];
-  snprintf (name, DVO_MAX_PATH, "%s/repair.warp.log", HOSTDIR);
-  logfile = fopen (name, "w");
-  if (!logfile) fprintf (stderr, "failed to open repair warp logfile %s\n", name);
-}
-
-void RepairWarpMeasuresCloseLogfile () {
-  if (!logfile) return;
-  fclose (logfile);
-}
-
-// XXX need to load_images first and generate some lookup tables
-int RepairWarpMeasures (Catalog *catalog) {
-
-  off_t Nimage;
-  Image *image = getimages (&Nimage, NULL);
-
-  myAssert (!catalog->Nmeasure || catalog->measure, "programming error");
-
-  int NfixChipID = 0, NfixStackID = 0, NfixWarpID = 0, NfixWarpImageID = 0, NmissWarp = 0, NmissStack = 0, NbadWarp = 0;
-
-  int onePercent = catalog->Nmeasure / 100;
-
-  for (off_t j = 0; j < catalog->Nmeasure; j++) {
-    if (j % onePercent == 0) fprintf (stderr, ".");
-
-    Measure *measure = &catalog->measure[j];
-
-    // get the associated average value:
-    int averef = measure->averef;
-    Average *average = &catalog->average[averef];
-
-    // double-check for consistency
-    myAssert (average->objID == measure->objID, "objID mismatch: %f %f : %d %d\n", average->R, average->D, average->objID, measure->objID);
-
-    // repair extID for non-warps:
-    if (isGPC1chip(measure->photcode))  {
-      double mjd = ohana_sec_to_mjd (measure->t);
-      int ccdnum = measure->photcode % 100;
-      uint64_t extID = CreatePSPSDetectionID (mjd, ccdnum, measure->detID);
-      if (extID != measure->extID) {
-	measure->extID = extID;
-	NfixChipID ++;
-      }
-      continue;
-    }
-
-    // repair extID for non-warps:
-    if (isGPC1stack(measure->photcode))  {
-      int im = getImageByID (measure->imageID);
-      if (im < 0) {
-	// if we don't have a matching image, we are probably assigned to an image on the wrong side of the sky
-	// reconstruct the imageID from the R,D,photcode & X,Y consistency
-	im = GetStackSeq (image, average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd);
-	if (im < 0) {
-	  // still failed, this is bad
-	  if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode);
-	  NmissStack ++;
-	  continue;
-	}
-	// the old imageID was wrong.  replace with the new one:
-	if (logfile) fprintf (logfile, "fix stack %10.6f %10.6f %5d %6.1f %6.1f 0x%08x %d -> %d\n", average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd, measure->detID, measure->imageID, image[im].imageID);
-	measure->imageID = image[im].imageID;
-      } else {
-	// we have a matching image, but check that it is correct by re-projecting the pixel position
-	double Xtst, Ytst;
-	RD_to_XY (&Xtst, &Ytst, average->R, average->D, &image[im].coords);
-	
-	// find the pixel offset
-	double dX = (Xtst - measure->Xccd);
-	double dY = (Ytst - measure->Yccd);
-	
-	// if dPos is small, we have the right image
-	double dPos = hypot(dX,dY);
-	if (dPos > 10.0) {
-	  // if dPos is large, we don't the right imageID, reconstruct the imageID from
-	  // the R,D,photcode & X,Y consistency
-	  im = GetStackSeq (image, average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd);
-	  if (im < 0) {
-	    // still failed, this is bad
-	    if (NmissStack < 10) fprintf (stderr, "missing stack exposure: %f %f : %d %d\n", measure->R, measure->D, measure->imageID, measure->photcode);
-	    NmissStack ++;
-	    continue;
-	  }
-	  // the old imageID was wrong.  replace with the new one:
-	  if (logfile) fprintf (logfile, "fix stack %10.6f %10.6f %5d %6.1f %6.1f 0x%08x %d -> %d\n", average->R, average->D, measure->photcode, measure->Xccd, measure->Yccd, measure->detID, measure->imageID, image[im].imageID);
-	  measure->imageID = image[im].imageID;
-	}
-      }
-      uint64_t extID = CreatePSPSStackDetectionID (35, image[im].externID, measure->detID);
-      if (extID != measure->extID) {
-	measure->extID = extID;
-	NfixStackID ++;
-      }
-      continue;
-    }
-
-    // we are only going to repair warp detections
-    if (!isGPC1warp(measure->photcode)) continue;
-
-    // warp coordinates to confirm warp
-    double X = measure->Xccd;
-    double Y = measure->Yccd;
-
-    // check if this detection is far from its correct location
-    double Rwrp = measure->R;
-    double Dwrp = measure->D;
-
-    double Rave = average->R;
-    double Dave = average->D;
-      
-    float csdec = cos(Dave * RAD_DEG);
-
-    // find the ra,dec displacement in arcsec:
-    double dR = 3600.0*(Rwrp - Rave)*csdec;
-    double dD = 3600.0*(Dwrp - Dave);
-
-    // skip detections which are within a small distance of the expected location
-    // NOTE: this actually works surprisingly well near the pole
-    double dPos = hypot(dR,dD);
-    if (dPos < 2.0) { 
-      int im = getImageByID (measure->imageID);
-      uint64_t extID = CreatePSPSStackDetectionID (34, image[im].externID, measure->detID);
-      if (extID != measure->extID) {
-	measure->extID = extID;
-	NfixWarpID ++;
-      }
-      continue;
-    }
-
-    // find the corrected warp ID:
-    int warpSeq = GetWarpSeq (image, measure->t, measure->photcode, Rave, Dave, X, Y);
-    if (warpSeq < 0) {
-      if (NmissWarp < 10) fprintf (stderr, "missing warp exposure: %f %f : %d %d\n", Rave, Dave, measure->t, measure->photcode);
-      NmissWarp ++;
-      continue;
-    }
-
-    // assert on coords being TAN?
-    Coords *imcoords = &image[warpSeq].coords;
-
-    double Rnew, Dnew;
-    XY_to_RD (&Rnew, &Dnew, X, Y, imcoords);
-
-    Rnew = ohana_normalize_angle_to_midpoint(Rnew, 180.0);
-    Dnew = ohana_normalize_angle_to_midpoint(Dnew,   0.0);
-
-    // find the ra,dec displacement in arcsec:
-    double dRnew = 3600.0*(Rnew - Rave)*csdec;
-    double dDnew = 3600.0*(Dnew - Dave);
-
-    // detections should now be repaired; fail if not
-    dPos = hypot(dRnew,dDnew);
-    if (dPos > 5.0) {
-      fprintf (stderr, "measurement still far from average location: %f %f vs %f %f : %d %d\n", Rave, Dave, Rnew, Dnew, measure->t, measure->photcode);
-      NbadWarp ++;
-    }
-
-    measure->imageID = image[warpSeq].imageID;
-    measure->R = Rnew;
-    measure->D = Dnew;
-
-    uint64_t extID = CreatePSPSStackDetectionID (34, image[warpSeq].externID, measure->detID);
-    if (extID != measure->extID) {
-      measure->extID = extID;
-      NfixWarpID ++;
-    }
-    NfixWarpImageID ++;
-  }
-
-  fprintf (stderr, "\n");
-  fprintf (stderr, "repaired %s : warp image ID %d : ext ID chip %d stack %d warp %d : missed warp %d stack %d : bad warp: %d\n", catalog->filename, NfixWarpImageID, NfixChipID, NfixStackID, NfixWarpID, NmissWarp, NmissStack, NbadWarp);
-
-  return (TRUE);
-}
Index: anches/eam_branches/relphot.20210521/src/StackImageMaps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/StackImageMaps.c	(revision 41603)
+++ 	(revision )
@@ -1,165 +1,0 @@
-# include "relphot.h"
-
-// generate an index of stack image ra,dec centers sorted indexes to the image IDs
-
-# define D_NSTACKS 1000
-# define NGROUPS 5
-
-typedef struct {
-  unsigned short photcode;
-  double *Rcenter;
-  double *Dcenter;
-  int    *imageID;
-  int    *imageSeq;
-  int Nstacks;
-  int NSTACKS;
-} StackGroupType;
-
-// one group for each photcode
-StackGroupType *stackgroup = NULL;
-
-int MakeStackIndex (void) {
-
-  off_t Nimage;
-  Image *image = getimages (&Nimage, NULL);
-
-  ALLOCATE (stackgroup, StackGroupType, NGROUPS);
-
-  for (int i = 0; i < NGROUPS; i++) {
-    stackgroup[i].photcode = 11000 + i*100;
-
-    stackgroup[i].Nstacks = 0;
-    stackgroup[i].NSTACKS = D_NSTACKS;
-    ALLOCATE (stackgroup[i].Rcenter,  double, stackgroup[i].NSTACKS);
-    ALLOCATE (stackgroup[i].Dcenter,  double, stackgroup[i].NSTACKS);
-    ALLOCATE (stackgroup[i].imageID,  int,    stackgroup[i].NSTACKS);
-    ALLOCATE (stackgroup[i].imageSeq, int,    stackgroup[i].NSTACKS);
-  }
-
-  // first generate a list of obstime values
-  for (off_t i = 0; i < Nimage; i++) {
-    if (!isGPC1stack(image[i].photcode)) continue;
-
-    int Ng = -1;
-    switch (image[i].photcode) {
-      case 11000: Ng = 0; break;
-      case 11100: Ng = 1; break;
-      case 11200: Ng = 2; break;
-      case 11300: Ng = 3; break;
-      case 11400: Ng = 4; break;
-      default: myAbort("impossible photcode");
-    }
-    myAssert (Ng < NGROUPS, "oops");
-
-    int N = stackgroup[Ng].Nstacks;
-    XY_to_RD (&stackgroup[Ng].Rcenter[N], &stackgroup[Ng].Dcenter[N], 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords);
-    stackgroup[Ng].Rcenter[N] = ohana_normalize_angle (stackgroup[Ng].Rcenter[N]);
-
-    stackgroup[Ng].imageSeq[N] = i;
-    stackgroup[Ng].imageID[N] = image[i].imageID;
-
-    stackgroup[Ng].Nstacks ++;
-
-    if (stackgroup[Ng].Nstacks >= stackgroup[Ng].NSTACKS) {
-      stackgroup[Ng].NSTACKS += D_NSTACKS;
-      REALLOCATE (stackgroup[Ng].Rcenter,  double, stackgroup[Ng].NSTACKS);
-      REALLOCATE (stackgroup[Ng].Dcenter,  double, stackgroup[Ng].NSTACKS);
-      REALLOCATE (stackgroup[Ng].imageID,  int,    stackgroup[Ng].NSTACKS);
-      REALLOCATE (stackgroup[Ng].imageSeq, int,    stackgroup[Ng].NSTACKS);
-    }
-  }
-
-  // sort the list:
-  for (int i = 0; i < NGROUPS; i++) { 
-    sort_by_ra (stackgroup[i].Rcenter, stackgroup[i].Dcenter, stackgroup[i].imageID, stackgroup[i].imageSeq, stackgroup[i].Nstacks);
-  }
-
-  return TRUE;
-}
-
-// find the stack which yields the given Rstk, Dstk for the given X,Y (and photcode)
-int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y) {
-
-  if (!stackgroup) return -1;
-
-  // we have the stack centers; find all stacks within 0.3 degrees of this point
-  Rstk = ohana_normalize_angle (Rstk);
-  double dD = 0.3;
-  double dR = dD / cos(RAD_DEG*Dstk);
-  double Rmin = Rstk - dR;
-  double Rmax = Rstk + dR;
-  double Dmin = Dstk - dD;
-  double Dmax = Dstk + dD;
-
-  int Ng = (photcode / 100) % 10;
-  int N = ohana_bisection_double (stackgroup[Ng].Rcenter, stackgroup[Ng].Nstacks, Rmin);
-
-  double dPosMin = NAN;
-  int    nPosMin = -1;
-
-  for (; N < stackgroup[Ng].Nstacks; N++) {
-    if (stackgroup[Ng].Dcenter[N] < Dmin) continue;
-    if (stackgroup[Ng].Dcenter[N] > Dmax) continue;
-    if (stackgroup[Ng].Rcenter[N] > Rmax) break;
-    
-    int im = stackgroup[Ng].imageSeq[N];
-
-    // project to image pixels
-    double Xtst, Ytst;
-    RD_to_XY (&Xtst, &Ytst, Rstk, Dstk, &image[im].coords);
-    
-    // find the pixel offset
-    double dX = (Xtst - X);
-    double dY = (Ytst - Y);
-    
-    // if dPos is small, we have the right image
-    double dPos = hypot(dX,dY);
-    if (dPos < 10.0) return im; // 10 pixels == 2.0 arcsec
-
-    // check for multiple possible matches??
-
-    // record the min pos so we can see how things fail
-    if (isnan(dPosMin)) {
-      dPosMin = dPos;
-      nPosMin = im;
-    } else {
-      if (dPos < dPosMin) {
-	dPosMin = dPos;
-	nPosMin = im;
-      }
-    }
-  }
-  
-  // we did not find the right image, return what we found (-1 or closest match)
-  return nPosMin;
-}
-
-void FreeStackGroups (void) {
-
-  if (!stackgroup) return;
-
-  for (int i = 0; i < NGROUPS; i++) {
-    free (stackgroup[i].Rcenter);
-    free (stackgroup[i].Dcenter);
-    free (stackgroup[i].imageID);
-    free (stackgroup[i].imageSeq);
-  }    
-  free (stackgroup);
-}
-
-void sort_by_ra (double *R, double *D, int *I, int *S, int N) {
-
-# define SWAPFUNC(A,B){ double dtmp; int itmp; 	\
-    dtmp = R[A]; R[A] = R[B]; R[B] = dtmp;	\
-    dtmp = D[A]; D[A] = D[B]; D[B] = dtmp;	\
-    itmp = I[A]; I[A] = I[B]; I[B] = itmp;	\
-    itmp = S[A]; S[A] = S[B]; S[B] = itmp;	\
-  }
-# define COMPARE(A,B)(R[A] < R[B])
-  
-  OHANA_SORT (N, COMPARE, SWAPFUNC);
-  
-# undef SWAPFUNC
-# undef COMPARE
-  
-}
Index: anches/eam_branches/relphot.20210521/src/WarpImageMaps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/WarpImageMaps.c	(revision 41603)
+++ 	(revision )
@@ -1,279 +1,0 @@
-# include "relphot.h"
-
-// we have a list of warp images, each with obstime & photcode.  I need to make a list of the obstime groups
-
-# define D_NOBSTIMES 10000
-
-typedef struct {
-  int obstime;
-  unsigned short photcode;
-  int Nwarps;
-  int NWARPS;
-  off_t *warps;
-  double *Rmin;
-  double *Rmax;
-  double *Dmin;
-  double *Dmax;
-  int    *onBoundary;
-} WarpGroup;
-
-WarpGroup *warpgroup = NULL;
-int Nwarpgroup = 0;
-int NWARPGROUP = 0;
-
-myIndexType *warpObstimeIndex = NULL;
-
-int FindWarpGroups (void) {
-
-  off_t Nimage;
-  Image *image = getimages (&Nimage, NULL);
-
-  int *obstimes = NULL;
-  int Nobstimes = 0;
-  int NOBSTIMES = D_NOBSTIMES;
-  ALLOCATE (obstimes, int, NOBSTIMES);
-
-  // first generate a list of obstime values
-  for (off_t i = 0; i < Nimage; i++) {
-    if (!isGPC1warp(image[i].photcode)) continue;
-
-    obstimes[Nobstimes] = image[i].tzero;
-    Nobstimes ++;
-
-    CHECK_REALLOCATE (obstimes, int, NOBSTIMES, Nobstimes, D_NOBSTIMES);
-  }
-  if (!Nobstimes) { 
-    free (obstimes);
-    return FALSE;
-  }
-
-  // sort the list:
-  isort (obstimes, Nobstimes);
-
-  // find the unique obstimes
-  int *uniqtimes = NULL;
-  int *uniqcount = NULL;
-  int Nuniqtimes = 0;
-  ALLOCATE (uniqtimes, int, Nobstimes);
-  ALLOCATE (uniqcount, int, Nobstimes);
-
-  // generate a uniq set of obstimes
-  for (int i = 0; i < Nobstimes; Nuniqtimes ++) {
-    uniqtimes[Nuniqtimes] = obstimes[i];
-    int Ndup = 0;
-    int lastValue = uniqtimes[Nuniqtimes];
-    while ((i < Nobstimes) && (obstimes[i] == lastValue)) {
-      i++;
-      Ndup ++;
-      uniqcount[Nuniqtimes] = Ndup;
-    }
-  }
-  REALLOCATE (uniqtimes, int, Nuniqtimes);
-  REALLOCATE (uniqcount, int, Nuniqtimes);
-  
-  // now I need to assign all warp images to a warpgroup
-  // I have 2 options to go from warp obstime to warpgroup:
-  // a) use a bisection lookup [may be slow]
-  // b) use an index on obstime [require ~5*3.14*1e7*4 bytes ~ 600MB for index]
-
-  // generate an index on obstime
-  warpObstimeIndex = myIndexInit ();
-  warpObstimeIndex->minID = obstimes[0];
-  warpObstimeIndex->maxID = obstimes[Nobstimes-1];
-
-  myIndexSetRange (warpObstimeIndex);
-
-  // generate the warp groups (create the index as we go)
-  Nwarpgroup = Nuniqtimes;
-  ALLOCATE (warpgroup, WarpGroup, Nwarpgroup);
-  
-  for (int i = 0; i < Nuniqtimes; i++) {
-    warpgroup[i].obstime = uniqtimes[i];
-    warpgroup[i].photcode = 0; // not yet set
-    warpgroup[i].Nwarps = 0;
-    warpgroup[i].NWARPS = 100;
-    ALLOCATE (warpgroup[i].warps, off_t, warpgroup[i].NWARPS);
-    myIndexSetEntry (warpObstimeIndex, uniqtimes[i], i);
-  }
-
-  // assign all warps to one of the warpgroups
-  for (off_t i = 0; i < Nimage; i++) {
-    if (!isGPC1warp(image[i].photcode)) continue;
-
-    int seq = myIndexGetEntry (warpObstimeIndex, image[i].tzero);
-    myAssert (warpgroup[seq].obstime == image[i].tzero, "oops");
-
-    if (!warpgroup[seq].photcode) {
-      warpgroup[seq].photcode = image[i].photcode;
-    } else {
-      myAssert (warpgroup[seq].photcode == image[i].photcode, "oops");
-    }
-
-    int N = warpgroup[seq].Nwarps;
-    warpgroup[seq].warps[N] = i;
-    warpgroup[seq].Nwarps ++;
-    
-    CHECK_REALLOCATE (warpgroup[seq].warps, off_t, warpgroup[seq].NWARPS, warpgroup[seq].Nwarps, 100);
-  }
-
-  // we now have the warps assigned to groups.  now we need to generate the grid data to find the warp assignments
-  for (int i = 0; i < Nwarpgroup; i++) {
-
-    myAssert (warpgroup[i].Nwarps == uniqcount[i], "failure");
-    ALLOCATE (warpgroup[i].Rmin, double, warpgroup[i].Nwarps);
-    ALLOCATE (warpgroup[i].Rmax, double, warpgroup[i].Nwarps);
-    ALLOCATE (warpgroup[i].Dmin, double, warpgroup[i].Nwarps);
-    ALLOCATE (warpgroup[i].Dmax, double, warpgroup[i].Nwarps);
-    ALLOCATE (warpgroup[i].onBoundary, int, warpgroup[i].Nwarps);
-
-    for (int j = 0; j < warpgroup[i].Nwarps; j++) {
-
-      off_t N = warpgroup[i].warps[j];
-      
-      int Nx = image[N].NX;
-      int Ny = image[N].NY;
-
-      warpgroup[i].onBoundary[j] = FALSE;
-
-      double R, D;
-      XY_to_RD (&R, &D, 0.0, 0.0, &image[N].coords);
-      R = ohana_normalize_angle_to_midpoint(R, 180.0);
-      warpgroup[i].Rmin[j] = R;
-      warpgroup[i].Rmax[j] = R;
-      warpgroup[i].Dmin[j] = D;
-      warpgroup[i].Dmax[j] = D;
-
-      // XXX need to worry about 0,360 boundary
-      XY_to_RD (&R, &D, Nx, 0.0, &image[N].coords);
-      R = ohana_normalize_angle_to_midpoint(R, 180.0);
-      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
-      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
-      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
-      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
-
-      // XXX need to worry about 0,360 boundary
-      XY_to_RD (&R, &D, 0.0, Ny, &image[N].coords);
-      R = ohana_normalize_angle_to_midpoint(R, 180.0);
-      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
-      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
-      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
-      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
-
-      // XXX need to worry about 0,360 boundary
-      XY_to_RD (&R, &D, Nx, Ny, &image[N].coords);
-      R = ohana_normalize_angle_to_midpoint(R, 180.0);
-      warpgroup[i].Rmin[j] = MIN (warpgroup[i].Rmin[j], R);
-      warpgroup[i].Rmax[j] = MAX (warpgroup[i].Rmax[j], R);
-      warpgroup[i].Dmin[j] = MIN (warpgroup[i].Dmin[j], D);
-      warpgroup[i].Dmax[j] = MAX (warpgroup[i].Dmax[j], D);
-
-      // bump Rmin,Rmax,Dmin,Dmax 10 arcsec worth of padding
-      double dR = 10.0/3600.0 / cos (RAD_DEG*warpgroup[i].Dmin[j]);
-      double dD = 10.0/3600.0;
-
-      // at north pole, force test of all nearby skycells
-      if (warpgroup[i].Dmax[j] > 89.8) {
-	warpgroup[i].Rmin[j] =   0.0;
-	warpgroup[i].Rmax[j] = 360.0;
-	warpgroup[i].Dmax[j] =  90.0;
-	warpgroup[i].Dmin[j] -=  dD;
-	continue;
-      } 
-
-      if (warpgroup[i].Rmax[j] - warpgroup[i].Rmin[j] > 270.0) {
-	// Rmin and Rmax are in the range 0 - 360.  For images at the 0,360 boundary,
-	// "Rmax" is the lower edge, and "Rmin" is the upper edge.  we need to flip them 
-	// and then break the 0-360 range:
-	double tmp = warpgroup[i].Rmin[j];
-	warpgroup[i].Rmin[j] = warpgroup[i].Rmax[j] - 360.0;
-	warpgroup[i].Rmax[j] = tmp;
-	warpgroup[i].onBoundary[j] = TRUE;
-      }
-
-      warpgroup[i].Rmin[j] -=  dR;
-      warpgroup[i].Rmax[j] +=  dR;
-      warpgroup[i].Dmin[j] -=  dD;
-      warpgroup[i].Dmax[j] +=  dD;
-    }
-  }
-
-  free (uniqtimes); 
-  free (uniqcount); 
-  free (obstimes);
-  return TRUE;
-}
-
-int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y) {
-
-  if (!warpObstimeIndex) return -1;
-
-  int seq = myIndexGetEntry (warpObstimeIndex, obstime);
-  myAssert (seq > -1, "ooops");
-  myAssert (warpgroup[seq].photcode == photcode, "oops");
-  
-  double dPosMin = NAN;
-  int nPosMin = -1;
-
-  // we now have the warp group, but which is the correct warp?
-  for (int i = 0; i < warpgroup[seq].Nwarps; i++) {
-    if (warpgroup[seq].onBoundary[i]) {
-      int inRange1 = (Rave >= warpgroup[seq].Rmin[i]) && (Rave <= warpgroup[seq].Rmax[i]);
-      int inRange2 = (Rave >= warpgroup[seq].Rmin[i] + 360.0) && (Rave <= warpgroup[seq].Rmax[i] + 360.0);
-      if (!inRange1 && !inRange2) continue;
-    } else {
-      if (Rave < warpgroup[seq].Rmin[i]) continue;
-      if (Rave > warpgroup[seq].Rmax[i]) continue;
-    }
-    if (Dave < warpgroup[seq].Dmin[i]) continue;
-    if (Dave > warpgroup[seq].Dmax[i]) continue;
-
-    // this is a possible image: check the coordinates:
-    int N = warpgroup[seq].warps[i];
-
-    // project Rave,Dave to image pixels
-    double Xtst, Ytst;
-    RD_to_XY (&Xtst, &Ytst, Rave, Dave, &image[N].coords);
-
-    // find the pixel offset
-    double dX = (Xtst - X);
-    double dY = (Ytst - Y);
-    
-    // skip detections which are within a small distance from the expected location
-    double dPos = hypot(dX,dY);
-    if (dPos < 20.0) return N; // 20.0 pixels = 5.0 arcsec
-
-    // check for multiple possible matches??
-
-    // record the min pos so we can see how things fail
-    if (isnan(dPosMin)) {
-      dPosMin = dPos;
-      nPosMin = i;
-    } else {
-      if (dPos < dPosMin) {
-	dPosMin = dPos;
-	nPosMin = i;
-      }
-    }
-  }
-  
-  // if we did not find any matched, give up
-  if (nPosMin < 0) return -1;
-
-  // return closest image
-  return warpgroup[seq].warps[nPosMin];
-}
-
-void FreeWarpGroups (void) {
-
-  for (int i = 0; i < Nwarpgroup; i++) {
-    free (warpgroup[i].warps);
-    free (warpgroup[i].Rmin);
-    free (warpgroup[i].Rmax);
-    free (warpgroup[i].Dmin);
-    free (warpgroup[i].Dmax);
-    free (warpgroup[i].onBoundary);
-  }    
-
-  free (warpgroup);
-  myIndexFree (warpObstimeIndex);
-}
Index: /branches/eam_branches/relphot.20210521/src/relphot.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/relphot.c	(revision 41603)
+++ /branches/eam_branches/relphot.20210521/src/relphot.c	(revision 41604)
@@ -34,5 +34,7 @@
   switch (mode) {
     case UPDATE_IMAGES:
+      // calculate zero points for images (may group my exposure [mosaic], night [tgroups]; may measure flat-correction)
       // IF CALLED WITH NLOOP == 0, DOES NOT LOAD bcatalog, just loads image table and generates ImageSubset.dat
+      // If called with -update, calls reload_catalogs() just like apply_offsets
       relphot_images (skylist);
       relphot_free (sky, skylist);
@@ -52,5 +54,6 @@
 
     case PARALLEL_REGIONS:
-      // run image updates in parallel across multiple remote machines
+      // calculate zero points for images (may group my exposure [mosaic], night [tgroups]; may measure flat-correction)
+      // equivalent to relphot_images, but run in parallel across multiple remote machines
       relphot_parallel_regions (sky);
       relphot_free (sky, skylist);
Index: /branches/eam_branches/relphot.20210521/src/relphot_images.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/relphot_images.c	(revision 41603)
+++ /branches/eam_branches/relphot.20210521/src/relphot_images.c	(revision 41604)
@@ -45,12 +45,5 @@
   if (!UPDATE) dvo_image_unlock (&db); 
 
-  // load the flat correction table (if defined)
-  // we no longer need to do this: the flatcorrection is now saved in measure.Mflat
-# if (0)
-  char flatcorrfile[256];
-  sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR);
-  FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE);
-# endif
-  FlatCorrectionTable *flatcorr = NULL;
+  FlatCorrectionTable *flatcorr = NULL; // no longer used here 
 
   initGridBins (); // allocates the empty array of corrections (elements are build below)
Index: /branches/eam_branches/relphot.20210521/src/relphot_objects.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/relphot_objects.c	(revision 41603)
+++ /branches/eam_branches/relphot.20210521/src/relphot_objects.c	(revision 41604)
@@ -43,12 +43,4 @@
   // load the ZP corrections here
   if (SYNTH_ZERO_POINTS) SynthZeroPointsLoad (SYNTH_ZERO_POINTS);
-
-  if (REPAIR_WARPS) { 
-    FindWarpGroups ();
-    MARKTIME("setup warp groups: %f sec\n", dtime);
-    MakeStackIndex ();
-    MARKTIME("setup stack index: %f sec\n", dtime);
-    RepairWarpMeasuresOpenLogfile ();
-  }
 
   // load data from each region file, only use bright stars
@@ -128,7 +120,4 @@
     findMosaics (&catalog, 1, FALSE);
 
-    // update the detection coordinates using the new image parameters
-    if (REPAIR_WARPS) RepairWarpMeasures (&catalog);
-
     initMrel (&catalog, 1);
     setMrelFinal (&catalog, NULL, TRUE);
@@ -153,13 +142,8 @@
     if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog.filename);
     
-    char history[128];
     struct timeval now;
     gettimeofday (&now, (void *) NULL);
     char *moddate = ohana_sec_to_date (now.tv_sec);
     gfits_modify (&catalog.header, "RELPHOT", "%s", 1, moddate);      
-    if (REPAIR_WARPS) {
-      snprintf (history, 128, "repair warp measure.imageID: %s", moddate);
-      gfits_modify_alt (&catalog.header, "HISTORY", "%S", 0, history); // adds a new entry
-    }
     free (moddate);
 
@@ -184,10 +168,4 @@
     freeMosaicBins (1, FALSE);
     freeGridBins ();
-  }
-
-  if (REPAIR_WARPS) {
-    FreeWarpGroups ();
-    FreeStackGroups ();
-    RepairWarpMeasuresCloseLogfile ();
   }
 
@@ -252,5 +230,4 @@
     if (RESET)         	   { strextend (&command, "-reset"); }
     if (RESET_ZEROPTS) 	   { strextend (&command, "-reset-zpts"); }
-    if (REPAIR_WARPS)      { strextend (&command, "-repair-warps"); }
     if (PRESERVE_PS1)      { strextend (&command, "-preserve-ps1"); }
     if (IS_DIFF_DB)        { strextend (&command, "-is-diff-db"); }
Index: /branches/eam_branches/relphot.20210521/test/relphot.grid.dvo
===================================================================
--- /branches/eam_branches/relphot.20210521/test/relphot.grid.dvo	(revision 41604)
+++ /branches/eam_branches/relphot.20210521/test/relphot.grid.dvo	(revision 41604)
@@ -0,0 +1,520 @@
+# -*-sh-*-
+
+$TAP_BREAK = 0
+$TAP_VERBOSE = 1
+input tap.dvo
+input ptolemy.dvo
+
+# set globals
+if (not($?PLOT)) set PLOT = 0
+
+$USE_CELL_OFFSETS = 1; # turn on non-zero cell offsets
+$SET_UC_RESID = 1; # write a non-zero resid value to the ubercal table (should not have an impact)
+$DO_2MASS = 1
+$DO_SYNTH = 1
+$DO_TYCHO = 1
+$NC_CLOUDS = 0.20
+$UC_CLOUDS = 0.02
+
+## full test set: generate a set of test files, insert into a catdir, check raw catalog, run relphot, check results
+macro go
+  if ($0 != 3)
+    echo "go (fileroot) (catdir)"
+    break
+  end
+
+  local i fileroot rootdir catdir
+
+  $fileroot = $1
+  $catdir = $2
+
+  dirname $1 -var rootdir
+  mkdir $rootdir
+
+  break -auto off
+  exec rm -rf $catdir
+  exec rm -rf $fileroot*
+  break -auto on
+
+  # define a set of photometric (ubercal) and non-photometric exposure times
+  init
+
+  # generate the common list of stars (coords & mags) for all exposures
+  mkstars 1000
+
+  catdir $catdir
+
+  tapPLAN {((mjd_uc[]*4 + mjd_nc[]*4) + (mjd_uc[] + mjd_nc[]) + (mjd_uc[] + mjd_nc[])) / 3}
+
+  # XXX first, run the analysis with only ubercal data (photometric data)
+  # we should be able to recover the nightly zero points fitting only tgroups
+
+  # generate the basic images and check they were correctly ingested by dvo
+  for i 0 mjd_uc[]
+    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 + dR_uc[$i]} {$DEC_CENTER + dD_uc[$i]}
+    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
+
+  # run relphot on the db and check that the images now match the expected values
+  echo relphot -grid -tgroup-fit-airmass -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 25
+  # 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
+
+  break
+
+  # reset the per-filter ZPT_OFF:$filter values
+  $ZPT_OFF:$FILTER = NAN
+
+  for i 0 mjd_uc[]
+    ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $FILTER relphot_uc
+  end
+  for i 0 mjd_nc[]
+    ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $FILTER relphot_nc
+  end
+  tapDONE
+end
+
+# set various global variables
+macro init
+
+  # generates ptolemy.rc, dvo.photcodes, $zpt_nominal:g,r $klam_nominal:g,r
+  mkptolemy.rc
+
+  $RA_CENTER = 10.0
+  $DEC_CENTER = 20.0
+  $PLATE_SCALE = 0.25
+  $NCHIP_X = 1
+  $NCHIP_Y = 1
+  $NCELL_X = 2
+  $NCELL_Y = 2
+  $CHIP_DX = 1000
+  $CHIP_DY = 1000
+  $DEFECT_FRAC = 0.00
+
+  # we have two sets of images: ubercal (photometric) and not-ubercaled
+
+  # ********* UBERCAL IMAGES ***********
+
+  # ubercal nights each have a distinct zero point and slope: 
+  vlist night_uc 55000.0 55010.0 55030.0 55060.0 55080.0 55090.0
+  vlist dz_nt_uc +0.0000 +0.0300 -0.0200 -0.0400 0.0100 0.0150
+# vlist dk_nt_uc +0.0000 +0.0005 -0.0005 -0.0002
+  vlist dk_nt_uc +0.0000 +0.0200 -0.0200 +0.0300 0.0005 0.0050
+
+  write tgroups.dat night_uc  
+
+  $EXPTIME = 10.0
+  $FILTER = g
+  $ZPT_OFF:$FILTER = NAN
+
+  delete -q mjd_uc seq_uc zpt_uc klam_uc secz_uc exptime_uc dR_uc dD_uc
+
+  # generate 5 exposures for each night
+  for i 0 16
+    set secz_tmp = rnd(night_uc) + 1.0; # airmass in range 1.0 - 2.0
+    set klam_tmp = $klam_nominal:$FILTER + dk_nt_uc
+    set exp_tmp  = zero(night_uc) + $EXPTIME
+    set zpt_tmp  = $zpt_nominal:$FILTER  + dz_nt_uc + 2.5*log(exp_tmp) + klam_tmp*(secz_tmp - 1.0)
+
+    set seq_tmp  = ramp(night_uc)
+    set mjd_tmp  = night_uc + 0.33 + 0.01*$i; # exposures are 14 minutes apart starting at 10pm HST
+    set dR_tmp   = zero(night_uc) + (0.25*($i % 4) - 0.5) * $PLATE_SCALE * $CHIP_DX / 3600.0 ; # offset in arcsec 
+    set dD_tmp   = zero(night_uc) + (0.25*int($i / 4) - 0.5) * $PLATE_SCALE * $CHIP_DY / 3600.0 ; # offset in arcsec 
+
+    concat klam_tmp klam_uc
+    concat secz_tmp secz_uc
+    concat  zpt_tmp  zpt_uc
+    concat  seq_tmp  seq_uc
+    concat  mjd_tmp  mjd_uc
+    concat   dR_tmp   dR_uc
+    concat   dD_tmp   dD_uc
+
+    concat  exp_tmp  exptime_uc
+  end    
+
+  # ********* NON-UBERCAL IMAGES ***********
+
+  # non-ubercal nights each have a random zero points and fixed slope (klam_nominal)
+  # non-ubercal nights must not be the same set as the ubercal nights above
+  vlist night_nc 55005.0 55015.0 55035.0 55065.0
+
+  delete -q mjd_nc seq_nc zpt_nc klam_nc secz_nc
+
+  # generate 5 exposures for each night
+  for i 0 5
+    set secz_tmp = rnd(night_nc)  + 1.0; # airmass in range 1.0 - 2.0
+    set klam_tmp = zero(night_nc) + $klam_nominal:$FILTER
+    set exp_tmp  = zero(night_nc) + $EXPTIME
+    set zpt_tmp  = $zpt_nominal:$FILTER + 2.5*log(exp_tmp) + klam_tmp*(secz_tmp - 1.0) - $NC_CLOUDS*(0.1 + rnd(secz_tmp))
+    # observed zero points for nc nights are 0.1 - 1.1 lower than the ubercal nights
+
+    set seq_tmp  = ramp(night_nc)
+    set mjd_tmp  = night_nc + 0.33 + 0.01*$i; # exposures are 14 minutes apart starting at 10pm HST
+
+    concat klam_tmp klam_nc
+    concat secz_tmp secz_nc
+    concat  zpt_tmp  zpt_nc
+    concat  seq_tmp  seq_nc
+    concat  mjd_tmp  mjd_nc
+
+    concat  exp_tmp  exptime_nc
+  end    
+end
+
+macro try.relphot
+  if ($0 != 3)
+    echo "go (fileroot) (catdir)"
+    break
+  end
+
+  local i fileroot rootdir catdir
+
+  $fileroot = $1
+  $catdir = $2
+
+  # run relphot on the db and check that the images now match the expected values
+  tapEXEC relphot -images g,r -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -statmode WT_MEAN -cloud-limit 0.5 -update -nloop 12
+
+  # reset the per-filter ZPT_OFF:$filter values
+  $ZPT_OFF:$FILTER = NAN
+
+  for i 0 mjd_uc[]
+    # ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $FILTER relphot_uc
+      ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $FILTER relphot_uc
+  end
+end
+
+macro ckexposure
+  if ($0 != 8)
+    echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)"
+    echo "  mode == raw or corr"
+    break
+  end
+
+  local CATDIR MJD_IMAGE ZPT_REAL ZPT_NOMINAL EXPTIME FILTER AIRMASS
+  
+  $CATDIR      = $1
+  $MJD_IMAGE   = $2
+  $ZPT_REAL    = $3
+  $EXPTIME     = $4
+  $AIRMASS     = $5
+  $FILTER      = $6
+  $FUNC_MODE   = $7
+
+  # XXX need a function to extract the nominal zpt for a given filter / photcode from the db
+  $ZPT_NOMINAL = $zpt_nominal:$FILTER
+  $KLAM_NOMINAL = $klam_nominal:$FILTER
+
+  $TIMEFORMAT = mjd
+  $TIMEREF = 0.0
+
+  catdir $CATDIR
+
+  # assume we still have stars_ra, stars_dec, stars_mag in hand
+  if ($PLOT)
+    dev -n 0
+    region $RA_CENTER $DEC_CENTER 0.2
+    images
+    pmeasure -all -m 15 20
+  end
+
+  # RA_MIN, etc are defined by mkstars
+  skyregion $RA_MIN $RA_MAX $DEC_MIN $DEC_MAX
+
+  ## calculate the average zero point offset
+  if (isnan($ZPT_OFF:$FILTER)) 
+    avextract ra dec $FILTER
+    match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2
+    reindex stars_tru_m = stars_mag using index1
+    set dm_tru = $FILTER - stars_tru_m
+    vstat -q dm_tru
+    $ZPT_OFF:$FILTER = $MEDIAN
+  end
+
+  mextract ra dec mag xccd yccd dbflags photflags time mag:sys mag:rel mag:inst mag:ave where (abs(time - $MJD_IMAGE) < 0.0001)
+ 
+  match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2
+
+  reindex stars_ra_m  = stars_ra  using index1
+  reindex stars_dec_m = stars_dec using index1
+  reindex stars_mag_m = stars_mag using index1
+  set dr = 3600*(ra - stars_ra_m)
+  set dd = 3600*(dec - stars_dec_m)
+  set dm = mag - stars_mag_m
+
+  set good_photom = (not(photflags & 0x02))
+  foreach field dm xccd yccd
+    subset $field\_good = $field if good_photom
+  end
+
+  if ($PLOT) 
+    dev -n 1
+    lim mag dm; clear; box; plot mag dm; plot -c red -pt 7 mag dm where (photflags & 0x00000002)
+    dev -n 2
+    lim xccd dm; clear; box; plot xccd dm; plot -c red -pt 7 xccd dm where (photflags & 0x00000002)
+    dev -n 3
+    lim ra dm; clear; box; plot ra dm; plot -c red -pt 7 ra dm where (photflags & 0x00000002)
+    echo "red points are known bad photometry"
+  end
+
+  $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0)
+
+  if ("$FUNC_MODE" == "raw") 
+    set dm_adjust = dm_good - $ZPT_NOMINAL + $ZPT_REAL_NORM
+
+    vstat -q dm_adjust
+    sprintf line "addstar raw ZP %8.2f $FILTER %6.3f %6.3f %6.3f (%7.4f): %7.4f %7.4f : %6.3f" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL {$ZPT_REAL_NORM - $ZPT_NOMINAL} $MEAN $SIGMA dm_adjust[]
+    tapOK {abs($MEAN) < 0.005} "RAW  : $line"
+    return
+  end
+
+  if (("$FUNC_MODE" == "relphot_uc") || ("$FUNC_MODE" == "relphot_nc"))
+
+    ## relphot does not set an absolute zero point, so we need to find that and subtract
+    set dm_adjust = dm_good - $ZPT_OFF:$FILTER
+
+    vstat -q dm_adjust
+    sprintf line "relphot     ZP %8.2f $FILTER %6.3f %6.3f %6.3f (%7.4f) : %7.4f %7.4f : %3d" $MJD_IMAGE $ZPT_REAL $ZPT_REAL_NORM $ZPT_NOMINAL {$ZPT_REAL_NORM - $ZPT_NOMINAL} $MEAN $SIGMA dm_adjust[]
+
+    if ("$FUNC_MODE" == "relphot_uc")
+      tapOK {abs($MEAN)  < 0.001} "MEAN : $line"
+      tapOK {abs($SIGMA) < 0.001} "SIGMA: $line"
+    end
+    if ("$FUNC_MODE" == "relphot_nc")
+      # XXX NOTE relaxed condition : bad photometry points are partly clipped but still cause a bias
+      tapOK {abs($MEAN)  < 0.020} "MEAN : $line"
+      tapOK {abs($SIGMA) < 0.020} "SIGMA: $line"
+    end
+
+    if ($PLOT)
+      cursor 1
+      if ("$KEY" == "Escape") break
+    end
+    return
+  end
+
+  echo "ERROR: unknown mode for ckexposure : $FUNC_MODE"
+end
+
+macro mkexposure
+ if ($0 != 10)
+   echo "mkexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (fileroot) (ra) (dec)"
+   break
+ end
+
+ local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR cmd
+
+ $CATDIR  = $1
+ $MJD     = $2
+ $ZPT     = $3
+ $EXPTIME = $4
+ $AIRMASS = $5 
+ $FILTER  = $6
+ $ROOT    = $7
+ $RAo     = $8
+ $DECo    = $9
+
+ $TIMEFORMAT = mjd
+ $TIMEREF = 0.0
+
+ # XX ctimes -ref $MJD -var datetime
+ # XX substr $datetime 0 10 date
+ # XX substr $datetime 11 8 time
+
+ for ix 0 $NCHIP_X
+  for iy 0 $NCHIP_Y
+    $dx  = $CHIP_DX * $PLATE_SCALE * ($ix - 0.5*$NCHIP_X + 0.5)
+    $dy  = $CHIP_DY * $PLATE_SCALE * ($iy - 0.5*$NCHIP_Y + 0.5)
+    $ra  = $RAo  - $dx / 3600.0 / dcos($DECo)
+    $dec = $DECo - $dy / 3600.0
+    # echo $ra $dec $dx $dy
+    mkinput $ROOT.$ix.$iy.in.txt $ra $dec $ZPT 
+    
+    # ra,dec is the center of this chip
+    local options
+    $options = -mjd $MJD 
+    $options = $options -radec $ra $dec 
+
+    #$options = $options -type PS1_V2 
+    #$options = $options -type PS1_V3 
+    #$options = $options -type PS1_V4 
+    $options = $options -type PS1_V5 
+
+    $options = $options -coords 
+    $options = $options -photcode GPC1.$FILTER.XY$ix\$iy 
+    $options = $options -no-noise 
+    $options = $options -flags
+    $options = $options -bad-psfqf-frac $DEFECT_FRAC
+    $options = $options -size $CHIP_DX $CHIP_DY 
+    $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} 
+    $options = $options -airmass $AIRMASS
+    $options = $options -exptime $EXPTIME
+
+    tapEXEC mkcmf  $ROOT.$ix.$iy.in.txt $ROOT.$ix.$iy.cmf $options
+
+    # the fake images have inconsistent ra,dec and airmass,sidtime values
+    tapEXEC addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass
+  end
+ end
+end
+
+# make a simple input file for mkcmf
+macro mkinput
+  if ($0 != 5)
+    echo "mkinput (filename) (ra) (dec) (zpt)"
+    break
+  end
+
+  # chip coordinate of the stars
+  # ra,dec is at the chip center, which corresponds to pixel (CHIP_DX,CHIP_DY)/2
+  set stars_X = (stars_ra  - $2) * 3600.0 * dcos($3) / $PLATE_SCALE + $CHIP_DX*0.5
+  set stars_Y = (stars_dec - $3) * 3600.0            / $PLATE_SCALE + $CHIP_DY*0.5
+  set stars_M = (stars_mag - $4)
+
+  set valid = (stars_X > 0) && (stars_X < $CHIP_DX) && (stars_Y > 0) && (stars_Y < $CHIP_DX)
+  subset stars_x  = stars_X if valid
+  subset stars_y  = stars_Y if valid
+  subset stars_m  = stars_M if valid
+  subset stars_r  = stars_ra if valid
+  subset stars_d  = stars_dec if valid
+  subset stars_mr = stars_mag if valid
+
+  # if we want to modify the observed magnitude based on the camera locations, add code here:
+
+  # cx,cy can have values of 0,1 (since NCELL_X, NCELL_Y = 2)
+  set stars_cx = int($NCELL_X * stars_x / $CHIP_DX)
+  set stars_cy = int($NCELL_Y * stars_y / $CHIP_DY)
+
+  # if I've done this right the flat-field offsets are
+  # (0,0) = 0.02, (0,1) = -0.02, (1,0) = 0.03, (1,1) = -0.03
+  set stars_cm = (stars_cx == 0) ? ((stars_cy == 0) ? zero(stars_cx) + 0.02 : zero(stars_cx) - 0.02) : ((stars_cy == 0) ? zero(stars_cx) + 0.03 : zero(stars_cx) - 0.03)
+  set stars_mi = stars_m + stars_cm
+
+  # create bad detections / detections to ignore or avoid, save to a file?
+
+  # create 3 types of things to ignore / avoid
+  # * badPhot (deviate the photometry)
+  # * badFlag (set a bad flag bit)
+  # * poorFlag (set a bad flag bit)
+  set badPhot  = rnd(stars_r) <  $DEFECT_FRAC
+  set badFlag  = rnd(stars_r) <  $DEFECT_FRAC
+  set poorFlag = rnd(stars_r) <  $DEFECT_FRAC
+
+  # create deviant photometry (1%)
+  # I'm setting the flag value 0x02 (EXTMODEL) if the object is deviated.  This is a simple way to pass the info
+  # note that outlier rejection only occurs for Nmeas > 5
+  set stars_mi = stars_mi + 0.2*badPhot
+  set stars_fl = zero(stars_mi) + 0x00000080*badFlag + 0x00000100*poorFlag + 0x00000002*badPhot
+
+  # bad flags that we sould respect:
+  # PM_SOURCE_MODE_FAIL             = 0x00000008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero)
+  # PM_SOURCE_MODE_SATSTAR          = 0x00000080, ///< Source model peak is above saturation
+  # PM_SOURCE_MODE_BADPSF           = 0x00000400, ///< Failed to get good estimate of object's PSF
+  # PM_SOURCE_MODE_DEFECT           = 0x00000800, ///< Source is thought to be a defect
+  # PM_SOURCE_MODE_SATURATED        = 0x00001000, ///< Source is thought to be saturated pixels (bleed trail)
+  # PM_SOURCE_MODE_CR_LIMIT         = 0x00002000, ///< Source has crNsigma above limit
+  # PM_SOURCE_MODE_MOMENTS_FAILURE  = 0x00008000, ///< could not measure the moments
+  # PM_SOURCE_MODE_SKY_FAILURE      = 0x00010000, ///< could not measure the local sky
+  # PM_SOURCE_MODE_SKYVAR_FAILURE   = 0x00020000, ///< could not measure the local sky variance
+  # PM_SOURCE_MODE_SIZE_SKIPPED     = 0x10000000, ///< size could not be determined
+
+  # poor flags that we should respect:
+  # PM_SOURCE_MODE_POOR             = 0x00000010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?)
+  # PM_SOURCE_MODE_PAIR             = 0x00000020, ///< Source fitted with a double psf
+  # PM_SOURCE_MODE_BLEND            = 0x00000100, ///< Source is a blend with other sources
+  # PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000, ///< moments not measured due to low S/N
+  # PM_SOURCE_MODE_BLEND_FIT        = 0x00400000, ///< source was fitted as a blend
+  # PM_SOURCE_MODE_ON_SPIKE         = 0x20000000, ///< peak lands on diffraction spike
+  # PM_SOURCE_MODE_ON_GHOST         = 0x40000000, ///< peak lands on ghost or glint
+  # PM_SOURCE_MODE_OFF_CHIP         = 0x80000000, ///< peak lands off edge of chip
+
+  write -f "%10.6f %10.6f %8.4f  %6.1f %6.1f %8.4f 0x%x" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi stars_fl
+end
+
+# generate a set of stars with raw RA, DEC, MAG values to use for all of the output cmf files
+macro mkstars
+  if ($0 != 2)
+    echo "USAGE: mkstars (Nstars)"
+    break
+  end
+
+  local RA_RANGE DEC_RANGE
+
+  # the images are oriented along N-S, E-W lines
+  # we make a dithered grid of images in each night
+
+  # size of region of interest in linear degrees
+  $RA_RANGE  = 2.5 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE / 3600
+  $DEC_RANGE = 2.5 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE / 3600 
+
+  $RA_MIN  =  $RA_CENTER - 0.5*$RA_RANGE/dcos($DEC_CENTER)
+  $RA_MAX  =  $RA_CENTER + 0.5*$RA_RANGE/dcos($DEC_CENTER)
+  $DEC_MIN = $DEC_CENTER - 0.5*$DEC_RANGE
+  $DEC_MAX = $DEC_CENTER + 0.5*$DEC_RANGE
+
+  create tmp 0 $1
+  set stars_ra  = $RA_CENTER  + $RA_RANGE  * (rnd(tmp) - 0.5) / dcos ($DEC_CENTER)
+  set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5)
+  set stars_mag = 10.0 + 10.0 * rnd(tmp) 
+  set stars_dmag = 0.02 + 0.1 * rnd(tmp)
+end
+
+macro plot.test.files
+  if ($0 != 2)
+    echo "USAGE: plot.test.files (N)"
+    break
+  end
+
+  clear -s
+  resize 2000 1200
+
+  section a0 0.00 0.50 0.25 0.50
+  section a1 0.25 0.50 0.25 0.50
+  section a2 0.50 0.50 0.25 0.50
+  section a3 0.75 0.50 0.25 0.50
+
+  section a4 0.00 0.00 0.25 0.50
+  section a5 0.25 0.00 0.25 0.50
+  section a6 0.50 0.00 0.25 0.50
+  section a7 0.75 0.00 0.25 0.50
+
+  label -fn helvetica 12
+  $BSTYLE = -tickpad 0.2 +xpad 0.2 +ypad 0.2 -ypad 5
+
+  for i 0 8
+    sprintf name tgrp.%05d.%02d.dat $i $1
+    data $name
+    exec ls -l --full-time $name
+    read Msys 4 Mrel 5 Mcal 7 Mmos 8 Mgrp 9 airmass 11 dMag 12
+  
+    section a$i
+    lim Msys dMag; box $BSTYLE; plot Msys dMag -pt cir -op 0.2 -c blue
+  end
+end
+
+macro plot.tcat.data
+  if ($0 != 2)
+    echo "USAGE: plot.tcat.data (Nlast)"
+    break
+  end
+
+  $BSTYLE = -tickpad 0.2 +xpad 0.2 +ypad 0.2 -ypad 5
+
+  for i 0 $1
+    sprintf name tcat.%02d.dat $i
+    data $name
+    dev -n $i
+    resize 900 500
+    clear;
+    box $BSTYLE
+    read Msys 4 Mrel 5 Mcal 7 Mmos 8 Mgrp 9 airmass 11 dMag 12;
+    lim Msys dMag; plot Msys dMag -pt cir -op 0.2 -c blue
+    echo $name
+  end
+end