Index: /tags/ipp-ps1-20210510/Ohana/src/kapa2/src/Center.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/kapa2/src/Center.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/kapa2/src/Center.c	(revision 41695)
@@ -17,6 +17,7 @@
 
   // enforce integer center here?
-  image[0].picture.Xc = X;
-  image[0].picture.Yc = Y;
+  // NAN value means retain existing center
+  if (isfinite(X)) image[0].picture.Xc = X;
+  if (isfinite(Y)) image[0].picture.Yc = Y;
   if ((zoom != 0) && (zoom != -1)) {
     image[0].picture.expand = zoom;
Index: /tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_get_column.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_get_column.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_get_column.c	(revision 41695)
@@ -310,5 +310,5 @@
    C - complex float
    M - complex double
-   P - var lenght array descrpt (64 bit)
+   P - var length array descrpt (64 bit)
    
    all can be preceeded by integer N which specified number
@@ -341,5 +341,4 @@
 
   if (!gfits_table_format (format, type, Nval, &Nbytes)) return (FALSE);
-  *Nval = 1;
 
   return (TRUE);
@@ -379,16 +378,25 @@
   gfits_scan (header, "NAXIS2",  OFF_T_FMT, 1,  &Ny);
 
-  /* scan columns to find insert point */
-  Nstart = 0;
-  for (i = 1; i < N; i++) {
-    snprintf (field, 256, "TFORM%d", i);
-    gfits_scan (header, field, "%s", 1, format);
-    gfits_table_format (format, tmp, &Nv, &Nb);
-    Nstart += Nv*Nb;
+  // find the starting byte for this column
+  // FITS ASCII table is supposed to have TBCOLn to specify the starting point
+  // but if it is missing, we can try to find by counting
+  snprintf (field, 256, "TBCOL%d", N);
+  int status = gfits_scan (header, field, "%d", 1, &Nstart);
+  if (!status) {
+    /* scan columns to find insert point */
+    Nstart = 0;
+    for (i = 1; i < N; i++) {
+      snprintf (field, 256, "TFORM%d", i);
+      gfits_scan (header, field, "%s", 1, format);
+      gfits_table_format (format, tmp, &Nv, &Nb);
+      Nstart += Nv*Nb;
+    }
+  } else {
+    Nstart --;
   }
 
   /* allocate temporary line, init pointers */
-  ALLOCATE (line, char, Nval+1);
-  bzero (line, Nval+1);
+  ALLOCATE (line, char, Nval*Nbytes+1);
+  bzero (line, Nval*Nbytes+1);
   Pin  = table[0].buffer + Nstart;
 
Index: /tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_set_column.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_set_column.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_set_column.c	(revision 41695)
@@ -478,22 +478,22 @@
   /* print line with appropriate formatting */
   ALLOCATE (array, char, Nbytes*Nval*Nrow);
-  ALLOCATE (line, char, Nval+1);
+  ALLOCATE (line, char, Nbytes+1);
   Pin = data;
   Pout = array;
   if (!strcmp (type, "char")) {
-    for (i = 0; i < Nval*Nrow; i++, Pin++, Pout++) {
+    for (i = 0; i < Nbytes*Nval*Nrow; i++, Pin++, Pout++) {
       *Pout = *Pin;
     }
   }
   if (!strcmp (type, "int")) {
-    for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nval) {
-      snprintf (line, Nval + 1, cformat, *(int *)Pin);
-      memcpy (Pout, line, Nval);
+    for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nbytes) {
+      snprintf (line, Nbytes + 1, cformat, *(int *)Pin);
+      memcpy (Pout, line, Nbytes);
     }
   }
   if (!strcmp (type, "float")) {
-    for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nval) {
-      snprintf (line, Nval + 1, cformat, *(float *)Pin);
-      memcpy (Pout, line, Nval);
+    for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nbytes) {
+      snprintf (line, Nbytes + 1, cformat, *(float *)Pin);
+      memcpy (Pout, line, Nbytes);
     }
   }
Index: /tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_table_format.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_table_format.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/libfits/table/F_table_format.c	(revision 41695)
@@ -5,5 +5,5 @@
 // get the format of a table column
 // Nval : number of joined columns
-// Nbytes : width of column
+// Nbytes : width of column, determined by the type (e.g., float = 4 bytes)
 int gfits_bintable_format (char *format, char *type, int *Nval, int *Nbytes) {
 
@@ -90,4 +90,7 @@
 
 /***********************/
+// get the format of a table column
+// Nval : number of fields (1 for all types except string)
+// Nbytes : width of FITS table column, ie number of bytes in ASCII table
 int gfits_table_format (char *format, char *type, int *Nval, int *Nbytes) {
 
@@ -109,9 +112,9 @@
   int isLong = FALSE;
   char Type = 'x';
-  if (Fchar == 'D') { *Nbytes = 1;  strcpy (type, "double"); Type = 'e'; *Nval = Nv; isLong = TRUE; }
-  if (Fchar == 'E') { *Nbytes = 1;  strcpy (type, "float");  Type = 'e'; *Nval = Nv; }
-  if (Fchar == 'F') { *Nbytes = 1;  strcpy (type, "float");  Type = 'f'; *Nval = Nv; }
-  if (Fchar == 'I') { *Nbytes = 1;  strcpy (type, "int");    Type = 'd'; *Nval = Nv; }
-  if (Fchar == 'A') { *Nbytes = 1;  strcpy (type, "char");   Type = 's'; *Nval = Nv; }
+  if (Fchar == 'D') { *Nbytes = Nv;  strcpy (type, "double"); Type = 'e'; *Nval =  1; isLong = TRUE; }
+  if (Fchar == 'E') { *Nbytes = Nv;  strcpy (type, "float");  Type = 'e'; *Nval =  1; }
+  if (Fchar == 'F') { *Nbytes = Nv;  strcpy (type, "float");  Type = 'f'; *Nval =  1; }
+  if (Fchar == 'I') { *Nbytes = Nv;  strcpy (type, "int");    Type = 'd'; *Nval =  1; }
+  if (Fchar == 'A') { *Nbytes =  1;  strcpy (type, "char");   Type = 's'; *Nval = Nv; }
   if (!*Nbytes) { return (FALSE); }
   
@@ -133,6 +136,6 @@
 /*
   valid TABLE column formats:
-  FN.N  - floating point
-  INN   - integer
+  FN.N  - floating point taking NN characters in the table
+  INN   - integer taking NN characters in the table
   ANN   - NN char string
   
Index: /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/Makefile	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/Makefile	(revision 41695)
@@ -42,4 +42,5 @@
 $(SRC)/fixcols.$(ARCH).o	   \
 $(SRC)/fiximage.$(ARCH).o	   \
+$(SRC)/forcedphot.$(ARCH).o	   \
 $(SRC)/gauss.$(ARCH).o		   \
 $(SRC)/getvel.$(ARCH).o	           \
@@ -90,7 +91,8 @@
 $(SRC)/imfit-qgauss-psf.$(ARCH).o	   \
 $(SRC)/imfit-sgauss.$(ARCH).o	   \
-$(SRC)/imfit-sgauss-psf.$(ARCH).o	   \
+$(SRC)/imfit-sgauss-psf.$(ARCH).o  \
 $(SRC)/imfit-qfgauss.$(ARCH).o	   \
-$(SRC)/imfit-qrgauss.$(ARCH).o	   
+$(SRC)/imfit-qrgauss.$(ARCH).o	   \
+$(SRC)/imfit-trail.$(ARCH).o
 
 # dependancy rules for include files ########################
Index: /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/forcedphot.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/forcedphot.c	(revision 41695)
+++ /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/forcedphot.c	(revision 41695)
@@ -0,0 +1,146 @@
+# include "astro.h"
+
+float Nsigma = 4.0;
+float sigma = 4.0 / 2.35; // FWHM = 4 pixels, 1 arcsec on PS1 GPC1
+float calc_forcedphot (Buffer *buf, float x, float y);
+float norm = 2.0;
+
+int forcedphot (int argc, char **argv) {
+
+  int N = 0;
+  Buffer *buf = NULL;
+  Vector *output = NULL;
+  Vector *xvec = NULL;
+  Vector *yvec = NULL;
+
+  int VERBOSE = TRUE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    VERBOSE = FALSE;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-quiet"))) {
+    VERBOSE = FALSE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-nsigma"))) {
+    remove_argument (N, &argc, argv);
+    Nsigma  = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  
+  if ((N = get_argument (argc, argv, "-sigma"))) {
+    remove_argument (N, &argc, argv);
+    sigma  = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  
+  int FLUXNORMAL = TRUE;
+  if ((N = get_argument (argc, argv, "-unit-normal"))) {
+    FLUXNORMAL = FALSE;
+    remove_argument (N, &argc, argv);
+  }
+
+  // XXX this should be the output name of a variable for 
+  char *outname = NULL;
+  if ((N = get_argument (argc, argv, "-output"))) {
+    remove_argument (N, &argc, argv);
+    outname = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  } else {
+    outname = strcreate ("fflux");
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: psfqf (buffer) x y [-sigma sigma] [-nsigma nsigma] [-output outname]\n");
+    gprint (GP_ERR, "  measures forced-photometry using defined PSF for supplied coordinate(s)\n");
+    gprint (GP_ERR, "  supplied coordinates may be scalar or vector(s)\n");
+    gprint (GP_ERR, "  output values are saved in $fflux or vector fflux, or supplied name\n");
+    goto escape;
+  }
+
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) goto escape;
+
+  if (FLUXNORMAL) {
+    norm = 2.0;
+  } else {
+    norm = 1.0 / (2.0 * M_PI * SQ(sigma));
+  }
+
+  // Scalar coordinate input:
+  double Xin = NAN;
+  double Yin = NAN;
+  if (SelectScalar (argv[2], &Xin)) {
+    if (!SelectScalar (argv[3], &Yin)) {
+      gprint (GP_ERR, "syntax error: mixed vector and scalars for (x y)\n");
+      goto escape;
+    }
+    float flux = calc_forcedphot (buf, Xin, Yin);
+    if (VERBOSE) gprint (GP_LOG, "forced flux at %10.6f %10.6f = %f\n", Xin, Yin, flux);
+    set_variable (outname, flux);
+    free (outname);
+    return TRUE;
+  }
+
+  // Vector coordinate input:
+  if ((xvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if ((yvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) goto escape;
+  if (xvec[0].Nelements != yvec[0].Nelements) {
+    fprintf (stderr, "mis-matched vector lengths\n");
+    goto escape;
+  }
+  REQUIRE_VECTOR_FLT (xvec, FALSE); 
+  REQUIRE_VECTOR_FLT (yvec, FALSE); 
+  
+  if ((output = SelectVector (outname, ANYBUFFER, TRUE)) == NULL) goto escape;
+  
+  if (output) ResetVector (output, OPIHI_FLT, xvec->Nelements);
+
+  for (int i = 0; i < xvec->Nelements; i++) {
+    Xin = xvec->elements.Flt[i];
+    Yin = yvec->elements.Flt[i];
+    float flux = calc_forcedphot (buf, Xin, Yin);
+    output->elements.Flt[i] = flux;
+  }
+  if (VERBOSE) gprint (GP_LOG, "forced flux calculated for %d points (output is %s)\n", output->Nelements, output->name);
+  
+  free (outname);
+  return (TRUE);
+
+ escape:
+  free (outname);
+  return FALSE;
+}
+
+// XXX things to consider:
+// if source is totally off image, flux should be nan, not 0.0
+// if any pixel is nan, skip or nan the source?
+float calc_forcedphot (Buffer *buf, float x, float y) {
+
+  int Nx = buf->matrix.Naxis[0];
+  int Ny = buf->matrix.Naxis[1];
+
+  int Xmin = x - Nsigma*sigma;
+  int Xmax = x + Nsigma*sigma;
+  int Ymin = y - Nsigma*sigma;
+  int Ymax = y + Nsigma*sigma;
+
+  float *bvalue = (float *)buf->matrix.buffer;
+  float f1 = -0.5 / SQ(sigma);
+
+  float sum = 0.0;
+  for (int iy = Ymin; iy <= Ymax; iy++) {
+    if (iy < 0) continue;
+    if (iy >= Ny) continue;
+    for (int ix = Xmin; ix <= Xmax; ix++) {
+      if (ix < 0) continue;
+      if (ix >= Nx) continue;
+      float r = SQ(ix - x) + SQ(iy - y);
+      
+      float gvalue = bvalue[ix + iy*Nx]*exp(f1*r);
+      sum += gvalue;
+    }
+  }
+  sum *= norm;
+  return sum;
+}
Index: /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/imfit-trail.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/imfit-trail.c	(revision 41695)
+++ /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/imfit-trail.c	(revision 41695)
@@ -0,0 +1,133 @@
+# include "imfit.h"
+
+opihi_flt trailTD (opihi_flt, opihi_flt, opihi_flt *, int, opihi_flt *);
+void  trailCL ();
+
+# define PAR_X      0
+# define PAR_Y      1
+# define PAR_THETA  2
+# define PAR_SIGMA  3
+# define PAR_LENGTH 4
+# define PAR_I0     5
+# define PAR_SKY    6
+
+void trail_setup (char *name) {
+
+  if (strcmp(name, "trail")) return;
+
+  fitfunc = trailTD;
+  imfit_cleanup = trailCL;
+  Npar = 7;
+  Nfpar = 0;
+
+  /* allocate free and fixed parameters */
+  ALLOCATE (par, opihi_flt, MAX (Npar, 1));
+  bzero (par, Npar*sizeof(opihi_flt));
+  ALLOCATE (fpar, opihi_flt, MAX (Nfpar, 1));
+  bzero (fpar, Nfpar*sizeof(opihi_flt));
+
+  par[PAR_X      ] = get_variable_default ("Xg",      0.0);
+  par[PAR_Y      ] = get_variable_default ("Yg",      0.0);
+  par[PAR_THETA  ] = get_variable_default ("Tg",      0.0);
+  par[PAR_SIGMA  ] = get_variable_default ("Wg",      2.0);
+  par[PAR_LENGTH ] = get_variable_default ("Lg",     10.0);
+  par[PAR_I0     ] = get_variable_default ("Zpk", 10000.0);
+  par[PAR_SKY    ] = get_variable_default ("Sg",      0.0);
+  sky = &par[PAR_SKY];
+}
+
+void trailCL () {
+  set_variable ("Xg",  par[PAR_X     ]);
+  set_variable ("Yg",  par[PAR_Y     ]);
+  set_variable ("Wg",  par[PAR_SIGMA ]);
+  set_variable ("Tg",  par[PAR_THETA ]);
+  set_variable ("Lg",  par[PAR_LENGTH]);
+  set_variable ("Zpk", par[PAR_I0    ]);
+  set_variable ("Sg",  par[PAR_SKY   ]);
+}
+
+/* trailed gaussian -- x, y, Sigma, Theta, Length, I, sky */
+opihi_flt trailTD (opihi_flt x, opihi_flt y, opihi_flt *par, int Npar, opihi_flt *dpar) {
+  OHANA_UNUSED_PARAM(Npar);
+
+  opihi_flt X = x - par[PAR_X];
+  opihi_flt Y = y - par[PAR_Y];
+  
+  opihi_flt S2 = 2.0 * SQ(par[PAR_SIGMA]);
+
+  opihi_flt ST = sin(RAD_DEG*par[PAR_THETA]);
+  opihi_flt CT = cos(RAD_DEG*par[PAR_THETA]);
+
+  opihi_flt Zp = (X*CT + Y*ST + 0.5*par[PAR_LENGTH]) / sqrt(S2);
+  opihi_flt Zm = (X*CT + Y*ST - 0.5*par[PAR_LENGTH]) / sqrt(S2);
+
+  opihi_flt Ep = erf(Zp);
+  opihi_flt Em = erf(Zm);
+
+  opihi_flt Rxy = Y*CT - X*ST;
+  opihi_flt Gxy = exp(-Rxy*Rxy/S2);
+
+  opihi_flt Pxy = Gxy * (Ep - Em);
+  opihi_flt f = Pxy * par[PAR_I0] + par[PAR_SKY];
+
+  if (dpar != NULL) {
+    dpar[PAR_SKY]    = 1.0;
+    dpar[PAR_I0]     = Pxy;
+
+    float dGdR = -2.0 * Rxy * Gxy / S2; // -R Gxy / (2 Sigma^2)
+
+    // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo
+    float dRdX = +ST;
+    float dRdY = -CT;
+    float dRdT = -Y*ST - X*CT;
+
+    float dGdX = dGdR * dRdX;
+    float dGdY = dGdR * dRdY;
+    float dGdT = dGdR * dRdT;
+
+    // are these signs correct? I think so: (dR/dXo = -dR/dX); dRdX below is actually dR/dXo
+    float dZpdX = -CT / sqrt(S2);
+    float dZmdX = dZpdX; // float dZmdX = -CT / sqrt(S2); dZmdX = dZpdX
+
+    float dZpdY = -ST / sqrt(S2); // float dZmdY = -ST / sqrt(S2); dZmdY = dZpdY
+    float dZmdY = dZpdY;
+
+    float dZpdL = +0.5 / sqrt(S2);
+    float dZmdL = -0.5 / sqrt(S2);
+
+    float dZpdT = (-X*ST + Y*CT) / sqrt(S2);
+    float dZmdT = dZpdT; // dZpdT = dZmdT
+
+    float dEdZp = exp (-Zp*Zp) * M_2_SQRTPI;
+    float dEdZm = exp (-Zm*Zm) * M_2_SQRTPI;
+
+    float dEpdX = dEdZp * dZpdX;
+    float dEmdX = dEdZm * dZmdX;
+
+    float dEpdY = dEdZp * dZpdY;
+    float dEmdY = dEdZm * dZmdY;
+
+    float dEpdL = dEdZp * dZpdL;
+    float dEmdL = dEdZm * dZmdL;
+
+    float dEpdT = dEdZp * dZpdT;
+    float dEmdT = dEdZm * dZmdT;
+
+    float dPdX = dGdX * (Ep - Em) + Gxy * (dEpdX - dEmdX);
+    float dPdY = dGdY * (Ep - Em) + Gxy * (dEpdY - dEmdY);
+
+    // dGdL is 0.0 because dRdL is 0.0
+    float dPdL = Gxy * (dEpdL - dEmdL);
+
+    float dPdT = dGdT * (Ep - Em) + Gxy * (dEpdT - dEmdT);
+
+    dpar[PAR_X]      = par[PAR_I0] * dPdX;
+    dpar[PAR_Y]      = par[PAR_I0] * dPdY;
+
+    dpar[PAR_LENGTH] = par[PAR_I0] * dPdL;
+    dpar[PAR_THETA]  = par[PAR_I0] * dPdT;
+    dpar[PAR_SIGMA]  = 0;	// we don't actually allow this to vary, so we do not need to calculate it
+  }
+
+  return (f);
+}
Index: /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/imfit.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/imfit.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/imfit.c	(revision 41695)
@@ -65,4 +65,5 @@
     qfgauss_setup (argv[N]);
     qrgauss_setup (argv[N]);
+    trail_setup (argv[N]);
     if (fitfunc == NULL) {
       gprint (GP_ERR, "unknown function %s\n", argv[N]);
@@ -124,8 +125,12 @@
   dchisq = ochisq;
   chisq  = ochisq;
-  for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
+
+//for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
+
+  for (i = 0; (i < 25); i++) {
     chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
     dchisq = ochisq - chisq;
     ochisq = chisq;
+    fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
   }  
   set_int_variable ("Niter",  i);
Index: /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/init.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.astro/init.c	(revision 41695)
@@ -27,4 +27,5 @@
 int fixcols                 PROTO((int, char **));
 int fixrows                 PROTO((int, char **));
+int forcedphot              PROTO((int, char **));
 int gauss                   PROTO((int, char **));
 int gaussfit                PROTO((int, char **));
@@ -100,4 +101,5 @@
   {1, "fixcols",     fixcols,      "fix bad columns by comparing with others"},
   {1, "fixrows",     fixrows,      "fix bad rows by comparing with others"},
+  {1, "forcedphot",  forcedphot,   "forced photometry on a star or set of stars, assuming gaussian profile"},
   {1, "gauss",       gauss,        "get statistics on a star, assuming gaussian profile"},
   {1, "getvel",      getvel,       "rotcurve to velocities"},
Index: /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.data/read_vectors.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.data/read_vectors.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.data/read_vectors.c	(revision 41695)
@@ -706,4 +706,7 @@
     char name[80];
       
+    // for both BINARY and ASCII tables, Nval is the number of fields of the given type
+    // e.g., there may be 3 related floating point fields fields or 21 connected ascii characters.
+
     Nval = 0;
     if (Binary) {
Index: /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.data/wd.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.data/wd.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/opihi/cmd.data/wd.c	(revision 41695)
@@ -1,3 +1,4 @@
 # include "data.h"
+FILE *gfits_open_and_extend (char *filename);
 
 int wd (int argc, char **argv) {
@@ -98,47 +99,14 @@
   // it updates NEXTEND and set EXTEND to TRUE, and modifies the PHU header (neither should happen)
   // if those keywords do not exist...
-  if (Extend) {
-    Header Xhead;
-    FILE *f;
-    off_t nbytes;
-    int status, Nextend;
-
-    /* assume failure means non-existent file */
-    if (!gfits_read_header (argv[2], &Xhead)) {
-
-      gfits_init_header (&Xhead);
-      Xhead.bitpix = 16;
-      Xhead.extend = TRUE;
-      gfits_create_header (&Xhead);
-
-      gfits_modify (&Xhead, "NEXTEND", "%d", 1, 0);
-      f = fopen (argv[2], "w");
-      fclose (f);
-    }
-
-    gfits_modify_alt (&Xhead, "EXTEND", "%t", 1, TRUE);
-
-    Nextend = 0;
-    gfits_scan (&Xhead, "NEXTEND", "%d", 1, &Nextend);
-    Nextend ++;
-    gfits_modify (&Xhead, "NEXTEND", "%d", 1, Nextend);
-
-    /* write the main header to the start of the file */
-    f = fopen (argv[2], "r+");
-    if (f == NULL) {
-      gprint (GP_ERR, "failed to write file\n");
+  if (Extend && !CompressMode) {
+    int status = TRUE;
+
+    FILE *f = gfits_open_and_extend (argv[2]);
+    if (!f) {
+      gprint (GP_ERR, "failed to extend file %s\n", argv[2]);
       status = FALSE;
       goto done1;
     }
-    
-    /* position to begining of file to write header */
-    fseeko (f, 0LL, SEEK_SET);
-    nbytes = fwrite (Xhead.buffer, 1, Xhead.datasize, f);
-    if (nbytes != Xhead.datasize) {
-      gprint (GP_ERR, "ERROR: failed writing data to image header\n");
-      status = FALSE;
-      goto done1;
-    }
-    
+
     /* fix up header */
     {
@@ -153,20 +121,23 @@
     /* position to end of file to write new extend */
     fseeko (f, 0LL, SEEK_END);
-    nbytes = fwrite (temp_header.buffer, 1, temp_header.datasize, f);
-    fclose (f);
+    off_t nbytes = fwrite (temp_header.buffer, 1, temp_header.datasize, f);
     if (nbytes != temp_header.datasize) {
+      fclose (f);
       gprint (GP_ERR, "failed to write file\n");
       status = FALSE;
       goto done1;
     }
+
     /* write the matrix buffer (automatically goes to end of file */
-    if (!gfits_write_matrix (argv[2], &temp_matrix)) {
+    if (!gfits_fwrite_matrix (f, &temp_matrix)) {
+      fclose (f);
       gprint (GP_ERR, "failed to write file\n");
       status = FALSE;
       goto done1;
     }
+    fclose (f);
     status = TRUE;
+
   done1:
-    gfits_free_header (&Xhead);
     gfits_free_header (&temp_header);
     gfits_free_matrix (&temp_matrix);
@@ -180,19 +151,28 @@
     Matrix myMatrix;
 
-    FILE *f = fopen (argv[2], "w");
-    if (!f) {
-      fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", argv[2]);
-      FREE (CompressMode);
-      return FALSE;
-    }
-    
-    gfits_init_header (&myHeader);
-    myHeader.extend = TRUE;
-    gfits_create_header (&myHeader);
-    gfits_create_matrix (&myHeader, &myMatrix);
-    gfits_fwrite_header  (f, &myHeader);
-    gfits_fwrite_matrix  (f, &myMatrix);
-    gfits_free_header (&myHeader);
-    gfits_free_matrix (&myMatrix);
+    FILE *f;
+
+    if (Extend) {
+      // keeps an existing file
+      f = gfits_open_and_extend (argv[2]);
+      fseeko (f, 0LL, SEEK_END);
+    } else {
+      // replaces an existing file
+      f = fopen (argv[2], "w");
+      if (!f) {
+	fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", argv[2]);
+	FREE (CompressMode);
+	return FALSE;
+      }
+   
+      gfits_init_header (&myHeader);
+      myHeader.extend = TRUE;
+      gfits_create_header (&myHeader);
+      gfits_create_matrix (&myHeader, &myMatrix);
+      gfits_fwrite_header  (f, &myHeader);
+      gfits_fwrite_matrix  (f, &myMatrix);
+      gfits_free_header (&myHeader);
+      gfits_free_matrix (&myMatrix);
+    }
 
     FTable ftable;
@@ -236,2 +216,84 @@
   return (TRUE);
 }
+
+// prepare the file to add an extension:
+// read the PHU, update NEXTEND, write back header
+// return the file pointer open and ready for write
+FILE *gfits_open_and_extend (char *filename) {
+
+  Header Xhead;
+  FILE *f = NULL;
+  struct stat filestat;
+
+  // does the file already exist?
+  int status = stat (filename, &filestat);
+  if (status == 0) { /* file exists, are permissions OK? */
+    f = fopen (filename, "r+");
+    if (f == NULL) {
+      gprint (GP_ERR, "failed to open file %s for write\n", filename);
+      return NULL;
+    }
+
+    status = gfits_fread_header (f, &Xhead);
+    if (!status) {
+      gprint (GP_ERR, "failed to read header for existing file %s\n", filename);
+      return NULL;
+    }
+  } else {
+    f = fopen (filename, "w");
+    if (f == NULL) {
+      gprint (GP_ERR, "failed to create file %s for write\n", filename);
+      return NULL;
+    }
+
+    // create an empty header
+    gfits_init_header (&Xhead);
+    Xhead.bitpix = 16;
+    Xhead.extend = TRUE;
+    gfits_create_header (&Xhead);
+
+    gfits_modify (&Xhead, "NEXTEND", "%d", 1, 0);
+  }
+
+  gfits_modify_alt (&Xhead, "EXTEND", "%t", 1, TRUE);
+
+  int Nextend = 0;
+  gfits_scan (&Xhead, "NEXTEND", "%d", 1, &Nextend);
+  Nextend ++;
+  gfits_modify (&Xhead, "NEXTEND", "%d", 1, Nextend);
+
+  /* position to begining of file to write header */
+  fseeko (f, 0LL, SEEK_SET);
+  off_t nbytes = fwrite (Xhead.buffer, 1, Xhead.datasize, f);
+  if (nbytes != Xhead.datasize) {
+    gfits_free_header (&Xhead);
+    gprint (GP_ERR, "ERROR: failed writing data to image header\n");
+    return NULL;
+  }
+  gfits_free_header (&Xhead);
+    
+  return (f);
+}
+
+
+///    /* fix up header */
+///    {
+///      static char simple[] = "XTENSION= 'IMAGE  '            / Image extension";
+///      int Ns, No;
+///      Ns = strlen (simple);
+///      No = 80 - Ns;
+///      strncpy (temp_header.buffer, simple, Ns);
+///      memset (&temp_header.buffer[Ns], ' ', No);
+///    }
+///
+///    /* position to end of file to write new extend */
+///    fseeko (f, 0LL, SEEK_END);
+///    nbytes = fwrite (temp_header.buffer, 1, temp_header.datasize, f);
+///    fclose (f);
+///    if (nbytes != temp_header.datasize) {
+///      gprint (GP_ERR, "failed to write file\n");
+///      status = FALSE;
+///      goto done1;
+///    }
+///
+///  }
Index: /tags/ipp-ps1-20210510/Ohana/src/opihi/include/imfit.h
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/opihi/include/imfit.h	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/opihi/include/imfit.h	(revision 41695)
@@ -19,2 +19,4 @@
 void qgauss_psf_setup (char *name);
 void sgauss_psf_setup (char *name);
+
+void trail_setup (char *name);
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/include/relphot.h
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/include/relphot.h	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/include/relphot.h	(revision 41695)
@@ -163,12 +163,16 @@
 
 typedef enum {
-  STATS_NONE, 
-  STATS_MEAN, 
-  STATS_MEDIAN, 
-  STATS_WT_MEAN, 
-  STATS_INNER_MEAN, 
-  STATS_INNER_WTMEAN, 
-  STATS_CHI_INNER_MEAN, 
-  STATS_CHI_INNER_WTMEAN
+  // these modes are primary and mutually exclusive:
+  STATS_NONE               = 0x0000, 
+  STATS_MEAN               = 0x0001, 
+  STATS_MEDIAN             = 0x0002, 
+  STATS_WT_MEAN            = 0x0003, 
+  STATS_INNER_MEAN         = 0x0004, 
+  STATS_INNER_WTMEAN       = 0x0005, 
+  STATS_CHI_INNER_MEAN     = 0x0006, 
+  STATS_CHI_INNER_WTMEAN   = 0x0007, 
+  // these modes are additional options
+  STATS_PRIMARY            = 0x0007, // use this to mask out the optional bits
+  STATS_VARSTATS           = 0x0010, 
 } ListStatsMode;
 
@@ -185,4 +189,8 @@
   double Lower20;
   double Lower10;
+  double Upper90Nsig;
+  double Upper80Nsig;
+  double Lower20Nsig;
+  double Lower10Nsig;
   double total;
   int    Nmeas;
@@ -419,4 +427,6 @@
 double IMFIT_SYS_SIGMA_LIM;
 double CLOUD_TOLERANCE;
+int    VARIABILITY_STATS;
+int    USE_OLS_FOR_AVERAGES;
 
 int    PARALLEL;
@@ -874,3 +884,9 @@
 void ResetAverageActivePhotcodes (SecFilt *secfilt);
 
-
+void   rationalize_zeropoints (int Niter);
+double get_median_zpt_images (short photcode);
+void   set_median_zpt_images (short photcode, double zpt);
+double get_median_zpt_mosaics (short photcode);
+void   set_median_zpt_mosaics (short photcode, double zpt);
+double get_median_zpt_tgroups (short photcode);
+void   set_median_zpt_tgroups (short photcode, double zpt);
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/src/GridOps.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/src/GridOps.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/src/GridOps.c	(revision 41695)
@@ -274,4 +274,6 @@
     if (!GridCorr[code]) continue;
 
+    // float GridSum = 0.0;
+    // int   GridCnt =   0;
     for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
       for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
@@ -292,6 +294,10 @@
 	GridCorr[code]->dMgrid[ix][iy] = sqrt(r*(Mgrid2 - Mgrid*Mgrid)); // sample stdev
 	// fprintf (stderr, "grid code %d, %d x %d : %f +/- %f : %d\n", code, ix, iy, GridCorr[code]-> Mgrid[ix][iy], GridCorr[code]->dMgrid[ix][iy], GridCorr[code]->nMgrid[ix][iy]);
+	// GridSum += Mgrid;
+	// GridCnt ++;
       }
     }
+    // float GridAve = GridSum / GridCnt;
+    // fprintf (stderr, "grid average: %f\n", GridAve);
   }
   return;
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/src/ImageOps.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/src/ImageOps.c	(revision 41695)
@@ -951,4 +951,104 @@
 }
 
+// find the median zero point and subtract it (for each active average photcode)
+void rationalize_zeropoints (int Niter) {
+
+  if (VERBOSE) fprintf (stderr, "rationalize zero points\n");
+
+  // if we are calculating zero points for tgroups,
+  // rationalize based on the tgroups
+  if (TGROUP_ZEROPT) {
+    for (int ic = 0; ic < Nphotcodes; ic++) {
+      double zpt = get_median_zpt_tgroups (photcodes[ic][0].code);
+      fprintf (stderr, "rationalize zero points by tgroup for %s (%d) : %f\n", photcodes[ic][0].name, photcodes[ic][0].code, zpt);
+      if (isnan(zpt)) continue;
+      set_median_zpt_tgroups (photcodes[ic][0].code, zpt);
+      set_median_zpt_mosaics (photcodes[ic][0].code, zpt);
+      set_median_zpt_images  (photcodes[ic][0].code, zpt);
+    }
+    return;
+  }
+
+  // if we are calculating zero points for mosaics, but not tgroups,
+  // rationalize based on the mosaics
+  if (MOSAIC_ZEROPT) {
+    for (int ic = 0; ic < Nphotcodes; ic++) {
+      double zpt = get_median_zpt_mosaics (photcodes[ic][0].code);
+      fprintf (stderr, "rationalize zero points by mosaic for %s (%d) : %f\n", photcodes[ic][0].name, photcodes[ic][0].code, zpt);
+      if (isnan(zpt)) continue;
+      set_median_zpt_tgroups (photcodes[ic][0].code, zpt);
+      set_median_zpt_mosaics (photcodes[ic][0].code, zpt);
+      set_median_zpt_images  (photcodes[ic][0].code, zpt);
+    }
+    return;
+  }    
+
+  // otherwise, rationalize based on the images
+  for (int ic = 0; ic < Nphotcodes; ic++) {
+    double zpt = get_median_zpt_images (photcodes[ic][0].code);
+    fprintf (stderr, "rationalize zero points by image for %s (%d) : %f\n", photcodes[ic][0].name, photcodes[ic][0].code, zpt);
+    if (isnan(zpt)) continue;
+    set_median_zpt_tgroups (photcodes[ic][0].code, zpt);
+    set_median_zpt_mosaics (photcodes[ic][0].code, zpt);
+    set_median_zpt_images  (photcodes[ic][0].code, zpt);
+  }
+}
+
+double get_median_zpt_images (short photcode) {
+
+  double *mlist;
+
+  ALLOCATE (mlist, double, Nimage);
+
+  int N = 0;
+  for (int i = 0; i < Nimage; i++) {
+    int mycode = GetPhotcodeEquivCodebyCode (image[i].photcode);
+    if (mycode != photcode) continue;
+    if (!(image[i].flags & ID_IMAGE_IMAGE_PHOTCAL)) continue;
+    mlist[N] = image[i].McalPSF;
+    N++;
+  }
+
+  if (N == 0) {
+    free (mlist);
+    return NAN;
+  }
+
+  StatType stats;
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (mlist, NULL, NULL, N, &stats);
+  free (mlist);
+
+  fprintf (stderr, "rationalize by image using %d pts\n", N);
+  return stats.median; 
+}
+
+void set_median_zpt_images (short photcode, double zpt) {
+
+  if (!isfinite(zpt)) return; // do not break the zero points
+
+  for (int i = 0; i < Nimage; i++) {
+    int mycode = GetPhotcodeEquivCodebyCode (image[i].photcode);
+    if (mycode != photcode) continue;
+    // fprintf (stderr, "IMAGE %d zpt %f -> ", i, image[i].McalPSF);
+
+    int applyOffset = TRUE;
+    TGroup *mygrp = getTGroupForImage (i);
+    if (mygrp && (mygrp->flags & ID_IMAGE_TGROUP_PHOTCAL)) applyOffset = FALSE;
+
+    Mosaic *mymos = getMosaicForImage (i);
+    if (mymos && (mymos->flags & ID_IMAGE_MOSAIC_PHOTCAL)) applyOffset = FALSE;
+
+    if (applyOffset) {
+      image[i].McalPSF -= zpt;
+      image[i].McalAPER -= zpt;
+    }
+    // fprintf (stderr, "%f (%d)\n", image[i].McalPSF, applyOffset);
+  }
+
+  return; 
+}
+
 static int setMcal_init_done = FALSE;
 void plot_setMcal (double *list, int Npts) {
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/src/MosaicOps.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/src/MosaicOps.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/src/MosaicOps.c	(revision 41695)
@@ -1574,4 +1574,71 @@
 }
 
+TGroup *getTGroupForMosaic (int mos) {
+  if (mos >= Nmosaic) return NULL;
+  if (mos < 0) return NULL;
+
+  if (!mosaic[mos].inTGroup) return NULL;
+
+  int imageIdx = MosaicToImage[mos][0];
+
+  TGroup *mygrp = getTGroupForImage (imageIdx);
+  return mygrp;
+}
+
+double get_median_zpt_mosaics (short photcode) {
+
+  double *mlist;
+
+  if (!MOSAIC_ZEROPT) return NAN;
+
+  ALLOCATE (mlist, double, Nmosaic);
+
+  int N = 0;
+  for (int i = 0; i < Nmosaic; i++) {
+    if (mosaic[i].photcode != photcode) continue;
+    if (!(mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL)) continue;
+    mlist[N] = mosaic[i].McalPSF;
+    N++;
+  }
+
+  if (N == 0) {
+    free (mlist);
+    return NAN;
+  }
+
+  StatType stats;
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (mlist, NULL, NULL, N, &stats);
+  free (mlist);
+
+  fprintf (stderr, "rationalize by mosaic using %d pts\n", N);
+  return stats.median; 
+}
+
+void set_median_zpt_mosaics (short photcode, double zpt) {
+
+  if (!MOSAIC_ZEROPT) return;
+  if (!isfinite(zpt)) return; // do not break the zero points
+
+  for (int i = 0; i < Nmosaic; i++) {
+    if (mosaic[i].photcode != photcode) continue;
+    // fprintf (stderr, "MOSAIC %d zpt %f -> ", i, mosaic[i].McalPSF);
+
+    int applyOffset = TRUE;
+    TGroup *mygrp = getTGroupForMosaic (i);
+    if (mygrp && (mygrp->flags & ID_IMAGE_TGROUP_PHOTCAL)) applyOffset = FALSE;
+    if (!(mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL)) applyOffset = FALSE;
+
+    if (applyOffset) {
+      mosaic[i].McalPSF -= zpt;
+      mosaic[i].McalAPER -= zpt;
+    }
+    // fprintf (stderr, "%f (%d)\n", mosaic[i].McalPSF, applyOffset);
+  }
+
+  return; 
+}
+
 void plot_mosaic_fields (Catalog *catalog) {
 
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/src/TGroupOps.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/src/TGroupOps.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/src/TGroupOps.c	(revision 41695)
@@ -1335,5 +1335,5 @@
       }
       // too few exposures (configure)
-      if (tgroup[j].Nmosaic < 4) {
+      if ((tgroup[j].Nmosaic < 4) && (tgroup[j].Nimage < 4)) {
 	mark = TRUE;
 	NfewExp ++;
@@ -1365,4 +1365,55 @@
 }
 
+double get_median_zpt_tgroups (short photcode) {
+
+  double *mlist;
+
+  if (!TGROUP_ZEROPT) return NAN;
+
+  ALLOCATE (mlist, double, NtgroupTimes);
+
+  int N = 0;
+  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].photcode != photcode) continue;
+      if (!(tgroup[j].flags & ID_IMAGE_TGROUP_PHOTCAL)) continue;
+      mlist[N] = tgroup[j].McalPSF;
+      N++;
+    }
+  }
+
+  if (N == 0) {
+    free (mlist);
+    return NAN;
+  }
+
+  StatType stats;
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (mlist, NULL, NULL, N, &stats);
+  free (mlist);
+
+  fprintf (stderr, "rationalize by tgroup using %d pts\n", N);
+  return stats.median; 
+}
+
+void set_median_zpt_tgroups (short photcode, double zpt) {
+
+  if (!TGROUP_ZEROPT) return;
+  if (!isfinite(zpt)) return; // do not break the zero points
+
+  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].photcode != photcode) continue;
+      if (!(tgroup[j].flags & ID_IMAGE_TGROUP_PHOTCAL)) continue;
+      tgroup[j].McalPSF -= zpt;
+      tgroup[j].McalAPER -= zpt;
+    }
+  }
+  return; 
+}
+
 void plot_tgroup_fields (Catalog *catalog) {
 
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/src/args.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/src/args.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/src/args.c	(revision 41695)
@@ -434,4 +434,15 @@
     remove_argument (N, &argc, argv);
     KEEP_UBERCAL = FALSE;
+  }
+
+  VARIABILITY_STATS = FALSE;
+  if ((N = get_argument (argc, argv, "-varstats"))) {
+    remove_argument (N, &argc, argv);
+    VARIABILITY_STATS = TRUE;
+  }
+  USE_OLS_FOR_AVERAGES = FALSE;
+  if ((N = get_argument (argc, argv, "-use-ols-for-averages"))) {
+    remove_argument (N, &argc, argv);
+    USE_OLS_FOR_AVERAGES = TRUE;
   }
 
@@ -865,4 +876,15 @@
 
 
+  VARIABILITY_STATS = FALSE;
+  if ((N = get_argument (argc, argv, "-varstats"))) {
+    remove_argument (N, &argc, argv);
+    VARIABILITY_STATS = TRUE;
+  }
+  USE_OLS_FOR_AVERAGES = FALSE;
+  if ((N = get_argument (argc, argv, "-use-ols-for-averages"))) {
+    remove_argument (N, &argc, argv);
+    USE_OLS_FOR_AVERAGES = TRUE;
+  }
+
   MIN_ERROR = 0.001;
   if ((N = get_argument (argc, argv, "-minerror"))) {
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/src/liststats.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/src/liststats.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/src/liststats.c	(revision 41695)
@@ -3,5 +3,5 @@
 void liststats_setmode (StatType *stats, char *strmode) {
 
-  stats->statmode = -1;
+  stats->statmode = STATS_NONE;
   if (!strcmp (strmode, "MEAN"))             { stats->statmode = STATS_MEAN; return; }
   if (!strcmp (strmode, "MEDIAN"))           { stats->statmode = STATS_MEDIAN; return; }
@@ -17,17 +17,21 @@
 
 int liststats_init (StatType *stats) {
-  stats->median  = NAN;
-  stats->mean    = NAN;
-  stats->sigma   = NAN;
-  stats->error   = NAN;
-  stats->chisq   = NAN;
-  stats->min     = NAN;
-  stats->max     = NAN;
-  stats->Upper80 = NAN;
-  stats->Lower20 = NAN;
-  stats->Upper90 = NAN;
-  stats->Lower10 = NAN;
-  stats->total   = NAN;
-  stats->Nmeas = 0;
+  stats->median      = NAN;
+  stats->mean        = NAN;
+  stats->sigma       = NAN;
+  stats->error       = NAN;
+  stats->chisq       = NAN;
+  stats->min         = NAN;
+  stats->max         = NAN;
+  stats->Upper80     = NAN;
+  stats->Lower20     = NAN;
+  stats->Upper90     = NAN;
+  stats->Lower10     = NAN;
+  stats->Upper80Nsig = NAN;
+  stats->Lower20Nsig = NAN;
+  stats->Upper90Nsig = NAN;
+  stats->Lower10Nsig = NAN;
+  stats->total       = NAN;
+  stats->Nmeas       = 0;
   return TRUE;
 }
@@ -39,4 +43,5 @@
 
   myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");
+  ListStatsMode myMode = stats->statmode & STATS_PRIMARY; // exclude the option bits
 
   liststats_init (stats);
@@ -75,5 +80,5 @@
   stats[0].Lower10 = value[N10];
 
-  switch (stats->statmode) {
+  switch (myMode) {
     case STATS_MEDIAN:
       ks = 0;
@@ -100,4 +105,5 @@
       break;
     case STATS_NONE:
+    default:
       myAbort ("undefined stats");
   }    
@@ -105,5 +111,5 @@
   // for these two modes, I need a vector of the chi-square contribution
   // I'm actually just using chisq to get the correct sorting order
-  if ((stats->statmode == STATS_CHI_INNER_MEAN) || (stats->statmode == STATS_CHI_INNER_WTMEAN)) {
+  if ((myMode == STATS_CHI_INNER_MEAN) || (myMode == STATS_CHI_INNER_WTMEAN)) {
     ALLOCATE (chi, double, N);
     for (i = 0; i < N; i++) {
@@ -119,7 +125,7 @@
 
   int WeightedMean = FALSE;
-  WeightedMean |= (stats->statmode == STATS_WT_MEAN);
-  WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN);
-  WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN);
+  WeightedMean |= (myMode == STATS_WT_MEAN);
+  WeightedMean |= (myMode == STATS_INNER_WTMEAN);
+  WeightedMean |= (myMode == STATS_CHI_INNER_WTMEAN);
   if (!dvalue) WeightedMean = FALSE;  // warn the user?
 
@@ -278,5 +284,4 @@
 
   // make the initial guess based on the median (not weighted mean)
-  // ALLOCATE_PTR (values, double, Npoints);
   for (int i = 0; i < Npoints; i++) {
     dataset->values[i] = dataset->flxlist[i];
@@ -288,8 +293,4 @@
   // if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, NULL, Npoints)) return FALSE;
   
-  // XXX add to dataset elements?
-  // ALLOCATE_PTR (wtvals, double, Npoints);
-  // ALLOCATE_PTR (wtlist, double, Npoints);
-
   int converged = FALSE;
   for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) {
@@ -315,20 +316,11 @@
       
   // calculate the weight thresholds to mask the bad points:
-  // double Sum_W = 0.0;
   for (int i = 0; i < Npoints; i++) {
     dataset->wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
     dataset->wtlist[i] = dataset->wtvals[i];
-    // Sum_W += dataset->wtvals[i];
   }
   dsort (dataset->wtlist, Npoints);
   double WtMedian = (Npoints % 2) ? dataset->wtlist[midpt] : 0.5*(dataset->wtlist[midpt] + dataset->wtlist[midpt-1]);
-  // double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
   double WtThreshold = WEIGHT_THRESHOLD * WtMedian;
-
-  // generate the unmasked subset
-  // XXX add these to the dataset elements?
-  // ALLOCATE_PTR ( ykeep, double, Npoints);
-  // ALLOCATE_PTR (dykeep, double, Npoints);
-  // ALLOCATE_PTR (wtkeep, double, Npoints);
 
   // save unmasked points
@@ -343,5 +335,5 @@
       continue;
     }
-     dataset->ykeep[Nkeep] = dataset->flxlist[i];
+    dataset-> ykeep[Nkeep] = dataset->flxlist[i];
     dataset->dykeep[Nkeep] = dataset->errlist[i];
     dataset->wtkeep[Nkeep] = dataset->wgtlist[i]; // externally-supplied weight
@@ -359,10 +351,36 @@
   stats->sigma = sqrt (dSig / (Nkeep - 1));
 
-  // bootstrap resampling to generate the errorbars
-  // XXX add these to the dataset elements?
-  // ALLOCATE_PTR (ysample,  double, Nkeep);
-  // ALLOCATE_PTR (dysample, double, Nkeep);
-  // ALLOCATE_PTR (wtsample, double, Nkeep);
-  // ALLOCATE_PTR (bvalue,   double, NBOOTSTRAP); // vector to save the bootstrap values
+  // if percentile ranges are desired, calculate them
+  if (stats->statmode & STATS_VARSTATS) {
+    // save ykeep values in a vector to be sorted
+    for (int i = 0; i < Nkeep; i++) {
+      dataset->wtvals[i] = dataset->ykeep[i];
+    }
+    dsort (dataset->wtvals, Nkeep);
+
+    int N90 = MIN (Nkeep-1, 0.9*Nkeep);
+    int N80 = MIN (Nkeep-1, 0.8*Nkeep);
+    int N20 = MAX (0, 0.2*Nkeep);
+    int N10 = MAX (0, 0.1*Nkeep);
+
+    stats->Upper90 = dataset->wtvals[N90];
+    stats->Upper80 = dataset->wtvals[N80];
+    stats->Lower20 = dataset->wtvals[N20];
+    stats->Lower10 = dataset->wtvals[N10];
+    
+    int NkeepMid = 0.5 * Nkeep;
+    double ClipMedian = (Nkeep % 2) ? dataset->wtlist[NkeepMid] : 0.5*(dataset->wtlist[NkeepMid] + dataset->wtlist[NkeepMid-1]);
+
+    // save ykeep values in a vector to be sorted
+    for (int i = 0; i < Nkeep; i++) {
+      dataset->wtvals[i] = (dataset->ykeep[i] - ClipMedian) / dataset->dykeep[i];
+    }
+    dsort (dataset->wtvals, Nkeep);
+
+    stats->Upper90Nsig = dataset->wtvals[N90];
+    stats->Upper80Nsig = dataset->wtvals[N80];
+    stats->Lower20Nsig = dataset->wtvals[N20];
+    stats->Lower10Nsig = dataset->wtvals[N10];
+  }
 
   int Nboot = 0;
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/src/reload_catalogs.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/src/reload_catalogs.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/src/reload_catalogs.c	(revision 41695)
@@ -237,4 +237,5 @@
     if (USE_BASIC_CHECK)   {  strextend (&command, "-basic-image-search"); }
     // if (USE_ALL_IMAGES)    { strextend (&command, "-use-all-images"); }
+    if (VARIABILITY_STATS) { strextend (&command, "-varstats"); }
     if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }
 
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/src/relphot_images.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/src/relphot_images.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/src/relphot_images.c	(revision 41695)
@@ -96,4 +96,6 @@
     for (i = 0; i < NLOOP; i++) {
       SetZptIteration (i);
+
+      rationalize_zeropoints (i);
 
       setMrel  (catalog, Ncatalog); // threaded (calls setMrelCatalog)
Index: /tags/ipp-ps1-20210510/Ohana/src/relphot/src/setMrelCatalog.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/relphot/src/setMrelCatalog.c	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/relphot/src/setMrelCatalog.c	(revision 41695)
@@ -1,3 +1,13 @@
 # include "relphot.h"
+
+# define UPPER90_VAL MpsfStk
+# define UPPER80_VAL FpsfStk
+# define LOWER20_VAL MkronStk
+# define LOWER10_VAL FkronStk
+# define UPPER90_SIG MpsfWrp 
+# define UPPER80_SIG FpsfWrp 
+# define LOWER20_SIG MkronWrp
+# define LOWER10_SIG FkronWrp
+
 
 # if (0)
@@ -71,5 +81,13 @@
   SetMrelInfoReset (results); // reset the counters
 
-  useOLS = UseStandardOLS (ZPT_STARS);
+  // XXX if we are iterating on the image zero points, we should use OLS
+  // until we have performed a few iterations.  UseStandardOLS() returns
+  // true (use OLS) or false depending on the zero point mode (TGROUP, MOSAIC, IMAGE) and the
+  // number of iterations.
+
+  // On the final pass (isSetMrelFinal), we should use OLS or IRLS depending on
+  // the user choice (USE_OLS_FOR_AVERAGES)
+
+  useOLS = isSetMrelFinal ? USE_OLS_FOR_AVERAGES : UseStandardOLS (ZPT_STARS);
   for (j = 0; j < catalog[Nc].Naverage; j++) {
 
@@ -435,4 +453,7 @@
     StatType *psfstats = &results->psfstats;
 
+    // if isSetMrelFinal and varstats are requested, request percentile ranges
+    if (isSetMrelFinal && VARIABILITY_STATS) psfstats->statmode |= STATS_VARSTATS;
+
     int NrankingPSF = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats);
     if (NrankingPSF < Nminmeas) { 
@@ -469,4 +490,15 @@
 	secfilt[Nsec].Mmax     = psfstats->max;
 	secfilt[Nsec].Mmin     = psfstats->min;
+      }
+      if (VARIABILITY_STATS) {
+	// XXX NOTE: these are overloaded on warp and stack mags and fluxes
+	secfilt[Nsec].UPPER90_VAL = psfstats->Upper90;
+	secfilt[Nsec].UPPER80_VAL = psfstats->Upper80;
+	secfilt[Nsec].LOWER20_VAL = psfstats->Lower20;
+	secfilt[Nsec].LOWER10_VAL = psfstats->Lower10;
+	secfilt[Nsec].UPPER90_SIG = psfstats->Upper90Nsig;
+	secfilt[Nsec].UPPER80_SIG = psfstats->Upper80Nsig;
+	secfilt[Nsec].LOWER20_SIG = psfstats->Lower20Nsig;
+	secfilt[Nsec].LOWER10_SIG = psfstats->Lower10Nsig;
       }
       secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
Index: /tags/ipp-ps1-20210510/Ohana/src/tools/Makefile
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/tools/Makefile	(revision 41694)
+++ /tags/ipp-ps1-20210510/Ohana/src/tools/Makefile	(revision 41695)
@@ -61,4 +61,25 @@
 $(PROGRAMS) : % : $(BIN)/%.$(ARCH)
 
+$(BIN)/mpcorb_predict.$(ARCH) : $(SRC)/mpcorb_predict.$(ARCH).o
+	@if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi
+	$(CC) $(FULL_CFLAGS) -o $@ $^ $(FULL_LDFLAGS) -lsla
+	@echo "compiled $*"
+	@echo ""
+
+$(DESTBIN)/mpcorb_predict: $(BIN)/mpcorb_predict.$(ARCH)
+	@if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi
+	rm -f $(DESTBIN)/mpcorb_predict
+	cp $(BIN)/mpcorb_predict.$(ARCH) $(DESTBIN)/mpcorb_predict
+	@echo "installed mpcorb_predict"
+	@echo ""
+
+$(DESTLIB)/libsla.a:
+	echo "missing SLALIB libsla.a: install in $(DESTLIB)"
+
+$(DESTINC)/slalib.h:
+	echo "missing SLALIB slalib.h: install in $(DESTINC)"
+
+mpcorb_predict: $(DESTBIN)/mpcorb_predict $(DESTLIB)/libsla.a $(DESTINC)/slalib.h
+
 %.clean:
 	rm -f $(SRC)/$*.$(ARCH).o
Index: /tags/ipp-ps1-20210510/Ohana/src/tools/src/mpcorb_predict.c
===================================================================
--- /tags/ipp-ps1-20210510/Ohana/src/tools/src/mpcorb_predict.c	(revision 41695)
+++ /tags/ipp-ps1-20210510/Ohana/src/tools/src/mpcorb_predict.c	(revision 41695)
@@ -0,0 +1,915 @@
+# include <ohana.h>
+# include <gfitsio.h>
+# include <convert.h>
+# include <slalib.h>
+
+/* convert an MPCORB.DAT file to a useful format that can be loaded by mknight */
+
+/* The MPCORB file is NOT fixed bytes/row.  We need to handle fractional lines at the end
+ * of each read block.
+ */
+
+/* read in chunks of ~64MB */
+# define NBYTE 0x400000
+# define DEBUG 0
+# define VERBOSE 1
+# define TESTING 0
+
+typedef struct {
+  char   ID[8];
+  double epoch;
+  double mean_anomaly;
+  double perihelion;
+  double ascend_node;
+  double inclination;
+  double eccentricity;
+  double semimajor_a;
+  double mjdMin;
+  double mjdMax;
+  double RatMin;
+  double RatMax;
+  double DatMin;
+  double DatMax;
+} Planets;
+
+int Shutdown (char *format, ...);
+int mpcorb_parseline (char *line, Planets *planet, int Nmax);
+Planets *mpcorb_read_text (char *filename, int *nplanets);
+Planets *mpcorb_read_fits (char *filename, int *nplanets);
+void     mpcorb_save_fits (char *filename, Planets *planets, int Nplanets);
+
+void mpcorb_predict (Planets *planet, double mjdObs, double *Robs, double *Dobs, int FullCalc);
+
+void mpcorb_trange (int argc, char **argv);
+void mpcorb_moment (int argc, char **argv);
+
+void mpcorb_testing (int argc, char **argv);
+
+void mpcorb_predict_test ();
+void mpcorb_predict_test_calc (Planets *planet, double mjdObs, double *Robs, double *Dobs);
+
+# define PS1_LONGITUDE -156.2559041668965 /* degrees East */
+# define PS1_LATITUDE    20.7070994446013 /* degrees North */
+# define PS1_ALTITUDE  3067.7 /* meters */
+
+int main (int argc, char **argv) {
+
+  if (TESTING) mpcorb_testing (argc, argv);
+
+  if ((argc != 9) && (argc != 10)) goto usage;
+
+  if (!strcasecmp (argv[1], "trange")) mpcorb_trange (argc, argv);
+  if (!strcasecmp (argv[1], "moment")) mpcorb_moment (argc, argv);
+
+usage:
+  fprintf (stderr, "USAGE: %s trange (MPCORB.DAT) (Tmin) (Tmax) (Rmin) (Rmax) (Dmin) (Dmax) (range.fits)\n", argv[0]);
+  fprintf (stderr, "USAGE: %s moment (range.fits) (Tobs) (Rmin) (Rmax) (Dmin) (Dmax) (output)\n", argv[0]);
+
+  fprintf (stderr, " trange : generate a table of asteroid orbital elements which fall within the region at Tmin or Tmax (in MJD)\n");
+  fprintf (stderr, " moment : generate a table of asteroid positions which fall within the region at Tobs (in MJD)\n");
+
+  exit (2);
+}
+
+// generate a table of asteroid orbital elements which fall within the region at Tmin or Tmax (in MJD)
+void mpcorb_trange (int argc, char **argv) {
+
+  if (argc != 10) Shutdown ("USAGE: %s trange (MPCORB.DAT) (Tmin) (Tmax) (Rmin) (Rmax) (Dmin) (Dmax) (range.fits)\n", argv[0]);
+
+  char  *filename = argv[2];
+  double mjdMin   = atof (argv[3]);
+  double mjdMax   = atof (argv[4]);
+  double Rmin     = atof (argv[5]);
+  double Rmax     = atof (argv[6]);
+  double Dmin     = atof (argv[7]);
+  double Dmax     = atof (argv[8]);
+  char  *output   = argv[9];
+
+  int Nplanets = 0;
+  Planets *planets = mpcorb_read_text (filename, &Nplanets);
+  if (VERBOSE) fprintf (stderr, "loaded %d planets\n", Nplanets);
+
+  int NplanetsSave = 0;
+  int NPLANETSSAVE = 1000;
+  ALLOCATE_PTR (planetsSave, Planets, NPLANETSSAVE);
+
+  // transform all objects and identify those in the target region:
+  for (int i = 0; i < Nplanets; i++) {
+    if (i % 1000 == 0) fprintf (stderr, ".");
+
+    // is the object in the region for the start of the period?
+    double RatMin, DatMin;
+    mpcorb_predict (&planets[i], mjdMin, &RatMin, &DatMin, FALSE);
+    int inRangeMin = TRUE;
+    inRangeMin = inRangeMin && (RatMin > Rmin);
+    inRangeMin = inRangeMin && (RatMin < Rmax);
+    inRangeMin = inRangeMin && (DatMin > Dmin);
+    inRangeMin = inRangeMin && (DatMin < Dmax);
+    
+    // is the object in the region for the end of the period?
+    double RatMax, DatMax;
+    mpcorb_predict (&planets[i], mjdMax, &RatMax, &DatMax, FALSE);
+    int inRangeMax = TRUE;
+    inRangeMax = inRangeMax && (RatMax > Rmin);
+    inRangeMax = inRangeMax && (RatMax < Rmax);
+    inRangeMax = inRangeMax && (DatMax > Dmin);
+    inRangeMax = inRangeMax && (DatMax < Dmax);
+
+    // NOTE: if an object is moving so fast that it traverses the entire region, we will
+    // miss it.  I'm going to ignore this case for now.
+    if (!inRangeMin && !inRangeMax) continue;
+
+    planetsSave[NplanetsSave] = planets[i];
+    planetsSave[NplanetsSave].mjdMin = mjdMin;
+    planetsSave[NplanetsSave].RatMin = RatMin;
+    planetsSave[NplanetsSave].DatMin = DatMin;
+    planetsSave[NplanetsSave].mjdMax = mjdMax;
+    planetsSave[NplanetsSave].RatMax = RatMax;
+    planetsSave[NplanetsSave].DatMax = DatMax;
+    
+    NplanetsSave++;
+    CHECK_REALLOCATE (planetsSave, Planets, NPLANETSSAVE, NplanetsSave, 1000);
+  }
+  fprintf (stderr, "\n");
+
+  mpcorb_save_fits (output, planetsSave, NplanetsSave);
+
+  exit (0);
+}
+
+// generate a table of asteroid positions which fall within the region at Tobs (in MJD)
+void mpcorb_moment (int argc, char **argv) {
+
+  if (argc != 9) Shutdown ("USAGE: %s moment (range.fits) (Tobs) (Rmin) (Rmax) (Dmin) (Dmax) (output)\n", argv[0]);
+
+  char  *filename = argv[2];
+  double mjdObs   = atof (argv[3]);
+  double Rmin     = atof (argv[4]);
+  double Rmax     = atof (argv[5]);
+  double Dmin     = atof (argv[6]);
+  double Dmax     = atof (argv[7]);
+  char  *output   = argv[8];
+
+  FILE *fout = fopen (output, "w");
+  if (!fout) Shutdown ("ERROR: unable to open file for output %s\n", output);
+
+  int Nplanets = 0;
+  Planets *planets = mpcorb_read_fits (filename, &Nplanets);
+  if (VERBOSE) fprintf (stderr, "loaded %d planets\n", Nplanets);
+
+  // transform all objects and identify those in the target region:
+  for (int i = 0; i < Nplanets; i++) {
+    if (i % 1000 == 0) fprintf (stderr, ".");
+
+    // is the object in the region for the start of the period?
+    // first, use the fast, but inaccurate, calculation
+    double Robs, Dobs;
+    mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, FALSE);
+    if (Robs < Rmin) continue;
+    if (Robs > Rmax) continue;
+    if (Dobs < Dmin) continue;
+    if (Dobs > Dmax) continue;
+    
+    // if it is in the region, use the more accurate calculation
+    mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, TRUE);
+    fprintf (fout, "%8s %12.6f %12.6f\n", planets[i].ID, Robs, Dobs);
+  }
+  fclose (fout);
+  fprintf (stderr, "\n");
+
+  exit (0);
+}
+
+// STATUS is value expected for success
+# define CHECK_STATUS(STATUS,MSG,...) if (!(STATUS)) { Shutdown (MSG, __VA_ARGS__); }
+
+// write minor planet orbital elements (plus (R,D,T)[min,max]) to a FITS table
+void mpcorb_save_fits (char *filename, Planets *planets, int Nplanets) {
+
+  Header header;
+  Header theader;
+  Matrix matrix;
+  FTable ftable;
+
+  gfits_init_header (&header);
+  header.extend = TRUE;
+  gfits_create_header (&header);
+  gfits_create_matrix (&header, &matrix);
+
+  gfits_create_table_header (&theader, "BINTABLE", "MPCORB");
+
+  gfits_define_bintable_column (&theader, "8A", "ID",           "name",         "none",     1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "EPOCH",        "epoch",        "MJD",      1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "MEAN_ANOMALY", "mean_anomaly", "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "PERIHELION",   "perihelion",   "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "ASCEND_NODE",  "ascend_node",  "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "INCLINATION",  "inclination",  "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "ECCENTRICITY", "eccentricity", "unitless", 1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "SEMIMAJOR_A",  "semimajor_a",  "AU",       1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "MJD_MIN",      "mjdMin",       "MJD",      1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "MJD_MAX",      "mjdMax",       "MJD",      1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "RAT_MIN",      "RatMin",       "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "RAT_MAX",      "RatMax",       "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "DAT_MIN",      "DatMin",       "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "DAT_MAX",      "DatMax",       "degrees",  1.0, 0.0);
+
+  // generate the output array that carries the data
+  gfits_create_table (&theader, &ftable);
+
+  // create intermediate storage arrays
+  ALLOCATE_PTR (ID,                     char,   Nplanets*8);
+  ALLOCATE_PTR (epoch,                  double, Nplanets);
+  ALLOCATE_PTR (mean_anomaly,	        double, Nplanets);
+  ALLOCATE_PTR (perihelion,	        double, Nplanets);
+  ALLOCATE_PTR (ascend_node,	        double, Nplanets);
+  ALLOCATE_PTR (inclination,	        double, Nplanets);
+  ALLOCATE_PTR (eccentricity,	        double, Nplanets);
+  ALLOCATE_PTR (semimajor_a,	        double, Nplanets);
+  ALLOCATE_PTR (mjdMin, 	        double, Nplanets);
+  ALLOCATE_PTR (mjdMax, 	        double, Nplanets);
+  ALLOCATE_PTR (RatMin, 	        double, Nplanets);
+  ALLOCATE_PTR (RatMax, 	        double, Nplanets);
+  ALLOCATE_PTR (DatMin, 	        double, Nplanets);
+  ALLOCATE_PTR (DatMax, 	        double, Nplanets);
+
+  // set intermediate storage arrays
+  for (int i = 0; i < Nplanets; i++) {
+    epoch[i]           = planets[i].epoch;
+    mean_anomaly[i]    = planets[i].mean_anomaly;
+    perihelion[i]      = planets[i].perihelion;
+    ascend_node[i]     = planets[i].ascend_node;
+    inclination[i]     = planets[i].inclination;
+    eccentricity[i]    = planets[i].eccentricity;
+    semimajor_a[i]     = planets[i].semimajor_a;
+    mjdMin[i]          = planets[i].mjdMin;
+    mjdMax[i]          = planets[i].mjdMax;
+    RatMin[i]          = planets[i].RatMin;
+    RatMax[i]          = planets[i].RatMax;
+    DatMin[i]          = planets[i].DatMin;
+    DatMax[i]          = planets[i].DatMax;
+
+    memcpy(&ID[i*8], planets[i].ID, 8);
+  }
+
+  // add the columns to the output array
+  gfits_set_bintable_column (&theader, &ftable, "ID",           ID,           Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "EPOCH",        epoch,        Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "MEAN_ANOMALY", mean_anomaly, Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "PERIHELION",   perihelion,   Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "ASCEND_NODE",  ascend_node,  Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "INCLINATION",  inclination,  Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "ECCENTRICITY", eccentricity, Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "SEMIMAJOR_A",  semimajor_a,  Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "MJD_MIN",      mjdMin,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "MJD_MAX",      mjdMax,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "RAT_MIN",      RatMin,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "RAT_MAX",      RatMax,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "DAT_MIN",      DatMin,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "DAT_MAX",      DatMax,       Nplanets);
+
+  // free intermediate storage arrays
+  free (ID);
+  free (epoch);
+  free (mean_anomaly);
+  free (perihelion);
+  free (ascend_node);
+  free (inclination);
+  free (eccentricity);
+  free (semimajor_a);
+  free (mjdMin);
+  free (mjdMax);
+  free (RatMin);
+  free (RatMax);
+  free (DatMin);
+  free (DatMax);
+
+  // open the output file
+  FILE *f = fopen (filename, "w");
+  if (!f) Shutdown ("ERROR: cannot open mpcorb file for output %s\n", filename);
+
+  int status;
+  status = gfits_fwrite_header  (f, &header);
+  CHECK_STATUS (status, "ERROR: cannot write header for mpcorbs %s\n", filename);
+
+  status = gfits_fwrite_matrix  (f, &matrix);
+  CHECK_STATUS (status, "ERROR: cannot write matrix for mpcorbs %s\n", filename);
+
+  status = gfits_fwrite_Theader (f, &theader);
+  CHECK_STATUS (status, "ERROR: cannot write table header for mpcorbs %s\n", filename);
+
+  status = gfits_fwrite_table  (f, &ftable);
+  CHECK_STATUS (status, "ERROR: cannot write table data for mpcorbs %s\n", filename);
+
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  gfits_free_header (&theader);
+  gfits_free_table (&ftable);
+
+  int fd = fileno (f);
+
+  status = fflush (f);
+  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
+
+  status = fsync (fd);
+  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
+
+  status = fclose (f);
+  CHECK_STATUS (!status, "ERROR: problem closing mpcorbs file %s\n", filename);
+
+  return;
+}
+
+# define GET_COLUMN(OUT,NAME,TYPE) \
+  TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
+  myAssert (!strcmp(type, #TYPE), "wrong column type");
+
+// load planet data from a fits table
+Planets *mpcorb_read_fits (char *filename, int *nplanets) {
+
+  int Ncol;
+  off_t Nrow;
+  Header header;
+  Header theader;
+  Matrix matrix;
+  FTable ftable;
+
+  FILE *f = fopen (filename, "r");
+  if (!f) Shutdown ("ERROR: cannot open mpcorb fits file %s\n", filename);
+
+  // load in PHU segment (ignore)
+  if (!gfits_fread_header (f, &header)) Shutdown ("can't read mpcorb fits file header %s\n", filename);
+  if (!gfits_fread_matrix (f, &matrix, &header)) Shutdown ("can't read mpcorb fits file matrix %s\n", filename);
+  ftable.header = &theader;
+
+  // load data for this header 
+  if (!gfits_load_header (f, &theader)) Shutdown ("can't read mpcorb fits file table header %s\n", filename);
+  if (!gfits_fread_ftable_data (f, &ftable, FALSE)) Shutdown ("can't read mpcorb fits file table data %s\n", filename);
+
+  // we read the entire block of data, then extract the columns, then set the Planets structure values.
+  // I free the FITS table data after extracting the colums to avoid having 3 copies in memory.
+
+  char type[16]; // used by the functions called by GET_COLUMN
+
+  GET_COLUMN (ID,           "ID",            char);
+  GET_COLUMN (epoch,        "EPOCH",         double);
+  GET_COLUMN (mean_anomaly, "MEAN_ANOMALY",  double);
+  GET_COLUMN (perihelion,   "PERIHELION",    double);
+  GET_COLUMN (ascend_node,  "ASCEND_NODE",   double);
+  GET_COLUMN (inclination,  "INCLINATION",   double);
+  GET_COLUMN (eccentricity, "ECCENTRICITY",  double);
+  GET_COLUMN (semimajor_a,  "SEMIMAJOR_A",   double);
+  GET_COLUMN (mjdMin,       "MJD_MIN",       double);
+  GET_COLUMN (mjdMax,       "MJD_MAX",       double);
+  GET_COLUMN (RatMin,       "RAT_MIN",       double);
+  GET_COLUMN (RatMax,       "RAT_MAX",       double);
+  GET_COLUMN (DatMin,       "DAT_MIN",       double);
+  GET_COLUMN (DatMax,       "DAT_MAX",       double);
+
+  // free the memory associated with the FITS files
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  gfits_free_header (&theader);
+  gfits_free_table (&ftable);
+
+  ALLOCATE_PTR (planets, Planets, Nrow);
+  for (int i = 0; i < Nrow; i++) {
+    planets[i].epoch        = epoch[i];
+    planets[i].mean_anomaly = mean_anomaly[i];
+    planets[i].perihelion   = perihelion[i];
+    planets[i].ascend_node  = ascend_node[i];
+    planets[i].inclination  = inclination[i];
+    planets[i].eccentricity = eccentricity[i];
+    planets[i].semimajor_a  = semimajor_a[i];
+    planets[i].mjdMin       = mjdMin[i];
+    planets[i].mjdMax       = mjdMax[i];
+    planets[i].RatMin       = RatMin[i];
+    planets[i].RatMax       = RatMax[i];
+    planets[i].DatMin       = DatMin[i];
+    planets[i].DatMax       = DatMax[i];
+
+    memcpy(planets[i].ID, &ID[i*8], 8);
+    planets[i].ID[7] = 0;
+  }
+  fprintf (stderr, "loaded data for %lld planets\n", (long long) Nrow);
+
+  free (ID          );
+  free (epoch       );
+  free (mean_anomaly);
+  free (perihelion  );
+  free (ascend_node );
+  free (inclination );
+  free (eccentricity);
+  free (semimajor_a );
+  free (mjdMin      );
+  free (mjdMax      );
+  free (RatMin      );
+  free (RatMax      );
+  free (DatMin      );
+  free (DatMax      );
+
+  *nplanets = Nrow;
+  return planets;
+}
+
+void mpcorb_testing (int argc, char **argv) {
+
+  if (argc != 3) Shutdown ("USAGE: %s (MPCORB.DAT) (epoch)\n", argv[0]);
+
+  char  *filename = argv[1];
+  double mjdObs   = atof (argv[2]);
+
+  int Nplanets = 0;
+  Planets *planets = mpcorb_read_text (filename, &Nplanets);
+  if (VERBOSE) fprintf (stderr, "loaded %d planets\n", Nplanets);
+
+  if (TESTING == 1) { mpcorb_predict_test(planets, Nplanets); exit (0); }
+
+  if (TESTING == 2) {
+    double Robs, Dobs;
+    for (int i = 0; i < 10; i++) {
+      mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, TRUE);
+      char ra_string[64], de_string[64];
+      hms_format (ra_string, 64, Robs/15.0);
+      hms_format (de_string, 64, Dobs);
+      fprintf (stderr, "result: %12.6f %12.6f : %s %s\n", Robs, Dobs, ra_string, de_string);
+    }
+    exit (0);
+  }
+  fprintf (stderr, "ERROR: invalid value of TESTING : %d\n", TESTING);
+  exit (2);
+}
+    
+// Supply this test with the full list of minor planets from MPCORB.DAT.
+// It will compare the positions for the first 10 minor planets with the
+// reference data downloaded from the MPC for:
+// MJD = 59396.0 @ (Long,Lat,Alt) = (0.0, 0.0, 0.0)
+
+void mpcorb_predict_test (Planets *planets, int Nplanets) {
+  
+  int    Nref = 10;
+  double Rref[10] = { 49.14708333, 352.90125000, 251.95416667, 175.74666667, 266.25541667, 298.38000000,  52.84458333, 121.94916667, 192.29500000, 144.20583333 };
+  double Dref[10] = { 11.48944444,   8.68722222,  -4.19666667,   9.37111111, -16.97583333,  -8.30500000,  22.13472222,  21.50888889,  -1.52000000,  11.57138889 };
+
+  for (int i = 0; (i < Nref) && (i < Nplanets); i++) {
+
+    if (VERBOSE) fprintf (stderr, "%8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.7f %12.7f\n",
+	     planets[i].ID,
+	     planets[i].epoch,
+	     planets[i].mean_anomaly,
+	     planets[i].perihelion,
+	     planets[i].ascend_node,
+	     planets[i].inclination,
+	     planets[i].eccentricity,
+	     planets[i].semimajor_a);
+    
+    // generate predictions for (Long,Lat,Alt) = (0.0, 0.0, 0.0)
+    double Rpre, Dpre;
+    mpcorb_predict_test_calc (&planets[i], 59396.0, &Rpre, &Dpre);
+
+    fprintf (stderr, "result: %8s : %12.6f %12.6f : %6.2f %6.2f\n", planets[i].ID, Rpre, Dpre, 3600*(Rref[i]-Rpre)*cos(Dref[i]*RAD_DEG), 3600*(Dref[i]-Dpre));
+  }
+}
+
+// attempt to predict as a test for Long,Lat = 0,0, with all other corrections
+// this attempt uses slaPertel to calculate perturbed versions of the orbital elements
+// NOTE this test seems to have worked: the calculated positions of the first 10 asteroids reported
+// by MPC match to within 2 arcsec.  This is much larger than the Planetary abberation,
+// so slaPlante must already be accounting for that effect.
+
+void mpcorb_predict_test_calc (Planets *planet, double mjdObs, double *Robs, double *Dobs) {
+
+  int jstat;
+  
+  double dTT = slaDtt (mjdObs);
+  double mjdObsTT = mjdObs + dTT / 86400.0;
+
+  Planets tmpPlanet;
+
+  // correct the orbital elements to the perturbed version
+  slaPertel (2,
+	     planet->epoch,
+	     mjdObsTT,
+	     planet->epoch,
+	     planet->inclination*RAD_DEG,
+	     planet->ascend_node*RAD_DEG,
+	     planet->perihelion*RAD_DEG,
+	     planet->semimajor_a,
+	     planet->eccentricity,
+	     planet->mean_anomaly*RAD_DEG,
+	     &tmpPlanet.epoch,
+	     &tmpPlanet.inclination,  // returned as radians and used as radians below
+	     &tmpPlanet.ascend_node,  // returned as radians and used as radians below
+	     &tmpPlanet.perihelion,   // returned as radians and used as radians below
+	     &tmpPlanet.semimajor_a,
+	     &tmpPlanet.eccentricity,
+	     &tmpPlanet.mean_anomaly, // returned as radians and used as radians below
+	     &jstat);
+
+  // this function returns the topocentric apparent R,D.  I assume it
+  // has no atm corrections.
+
+  double Rtopo, Dtopo, dist;
+  slaPlante (mjdObsTT,
+	     0.0, // PS1_LONGITUDE*RAD_DEG,
+	     0.0, // PS1_LATITUDE*RAD_DEG,
+	     2,
+	     tmpPlanet.epoch,
+	     tmpPlanet.inclination,
+	     tmpPlanet.ascend_node,
+	     tmpPlanet.perihelion,
+	     tmpPlanet.semimajor_a,
+	     tmpPlanet.eccentricity,
+	     tmpPlanet.mean_anomaly,
+	     0.0, &Rtopo, &Dtopo, &dist, &jstat);
+
+  if (VERBOSE) fprintf (stderr, "Plante   : %12.6f %12.6f\n", Rtopo*DEG_RAD, Dtopo*DEG_RAD);
+
+  // NOTE: the atm effects seem to be large and incompatible with MPC outputs.
+  // Setting the pressure to 0.0 seems to avoid those corrections.
+
+  double Rap, Dap;
+  slaOap ("R",
+	  Rtopo,
+	  Dtopo,
+	  mjdObs,
+	  -0.4,
+	  0.0, // LONGITUDE
+	  0.0, // LATITUDE
+	  0.0, // ALTITUDE,
+	  0.0, // xp
+	  0.0, // yp
+	  273.15, // T(K)
+	  0.0, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
+	  0.0, // RH(%)
+	  0.7, // Wavelength (um)
+	  0.0, // lapse rate
+	  &Rap, &Dap);
+
+  if (VERBOSE) fprintf (stderr, "Oap   : %12.6f %12.6f : %6.2f %6.2f\n", Rap*DEG_RAD, Dap*DEG_RAD, 3600*(Rap-Rtopo)*DEG_RAD*cos(Dtopo), 3600*(Dap-Dtopo)*DEG_RAD);
+
+  double Rmean, Dmean;
+  slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
+
+  if (VERBOSE) fprintf (stderr, "Amp   : %12.6f %12.6f : %6.2f %6.2f\n", Rmean*DEG_RAD, Dmean*DEG_RAD, 3600*(Rmean-Rap)*DEG_RAD*cos(Dap), 3600*(Dmean-Dap)*DEG_RAD);
+
+  if (0) {
+    // slaAmp applies both precession/nutation as well as abberation
+    // but, I think abberation is not appropriate to solar system bodies.
+
+    // here, I instead calculate just the precession and nutation with slalib
+
+    double Rpre, Dpre;
+    Rpre = Rap;
+    Dpre = Dap;
+    slaPreces ("FK5", 2021.5, 2000.000000, &Rpre, &Dpre);
+    if (VERBOSE) fprintf (stderr, "Prec     : %12.6f %12.6f : %6.2f %6.2f\n", Rpre*DEG_RAD, Dpre*DEG_RAD, 3600*(Rpre-Rap)*DEG_RAD*cos(Dap), 3600*(Dpre-Dap)*DEG_RAD);
+  }
+
+  *Robs = Rmean * DEG_RAD;
+  *Dobs = Dmean * DEG_RAD;
+}
+
+// mjdObs is NOT in TT
+void mpcorb_predict (Planets *planet, double mjdObs, double *Robs, double *Dobs, int FullCalc) {
+
+  int jstat;
+  
+  double dTT = slaDtt (mjdObs);
+  double mjdObsTT = mjdObs + dTT / 86400.0;
+
+  Planets tmpPlanet;
+
+  // correct the orbital elements to the perturbed version
+  if (FullCalc) { slaPertel (2,
+	     planet->epoch,
+	     mjdObsTT,
+	     planet->epoch,
+	     planet->inclination*RAD_DEG,
+	     planet->ascend_node*RAD_DEG,
+	     planet->perihelion*RAD_DEG,
+	     planet->semimajor_a,
+	     planet->eccentricity,
+	     planet->mean_anomaly*RAD_DEG,
+	     &tmpPlanet.epoch,
+	     &tmpPlanet.inclination,  // returned as radians and used as radians below
+	     &tmpPlanet.ascend_node,  // returned as radians and used as radians below
+	     &tmpPlanet.perihelion,   // returned as radians and used as radians below
+	     &tmpPlanet.semimajor_a,
+	     &tmpPlanet.eccentricity,
+	     &tmpPlanet.mean_anomaly, // returned as radians and used as radians below
+	     &jstat);
+  } else {
+    tmpPlanet = *planet;
+  }
+
+  // this function returns the topocentric apparent R,D.  I assume it
+  // has no atm corrections.
+
+  double Rtopo, Dtopo, dist;
+  slaPlante (mjdObsTT,
+	     PS1_LONGITUDE*RAD_DEG,
+	     PS1_LATITUDE*RAD_DEG,
+	     2,
+	     planet->epoch, 
+	     planet->inclination*RAD_DEG,
+	     planet->ascend_node*RAD_DEG,
+	     planet->perihelion*RAD_DEG,
+	     planet->semimajor_a,
+	     planet->eccentricity,
+	     planet->mean_anomaly*RAD_DEG,
+	     0.0, &Rtopo, &Dtopo, &dist, &jstat);
+
+  // NOTE: the atm effects seem to be large and incompatible with MPC outputs.
+  // Setting the pressure to 0.0 seems to avoid those corrections.
+
+  double Rap, Dap;
+  if (FullCalc) {
+    slaOap ("R",
+	  Rtopo,
+	  Dtopo,
+	  mjdObs,
+	  -0.4,
+	  PS1_LONGITUDE*RAD_DEG,
+	  PS1_LATITUDE*RAD_DEG,
+	  PS1_ALTITUDE,
+	  0.0, // xp
+	  0.0, // yp
+	  273.15, // T(K)
+	  0.0, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
+	  0.0, // RH(%)
+	  0.7, // Wavelength (um)
+	  0.0065, // lapse rate
+	  &Rap, &Dap);
+  } else { Rap = Rtopo; Dap = Dtopo; }
+
+  double Rmean, Dmean;
+  if (FullCalc) {
+    slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
+  } else {Rmean = Rap; Dmean = Dap; }
+
+  *Robs = Rmean * DEG_RAD;
+  *Dobs = Dmean * DEG_RAD;
+}
+
+# define MPC_DOUBLE(FIELD,START,END) {FIELD = strtod(&line[START-1], NULL);}
+
+// part the data in a single line of the MPCORB.DAT file
+int mpcorb_parseline (char *line, Planets *planet, int Nmax) {
+
+  // fields I need:
+
+  // epoch       : 21 - 25
+  // inclination : 60 - 68
+  // longitude of ascending node (J2000) : 49 - 57
+  // argument of perihelion (J2000) : 38 - 46
+  // semimajor axis (AU) : 93 - 103
+  // eccentricity : 71 - 79
+  // Mean anomaly : 27 - 35
+  
+  // these fields are saved in the data file with whitespace between, so I can
+  // set the byte after the range to NULL and scanf the values.  or just use strtod as needed:
+
+  strncpy (planet->ID, line, 7); planet->ID[7] = 0;
+
+  MPC_DOUBLE (planet->mean_anomaly, 27, 35);
+  MPC_DOUBLE (planet->perihelion,   38, 46);
+  MPC_DOUBLE (planet->ascend_node,  49, 57);
+  MPC_DOUBLE (planet->inclination,  60, 68);
+  MPC_DOUBLE (planet->eccentricity, 71, 79);
+  MPC_DOUBLE (planet->semimajor_a,  93, 103);
+
+  // the epoch is stored in a compressed format:
+  int year = 1800 + 100*(line[20] - 'I') + 10*(line[21] - '0') + line[22] - '0';
+  char mchar = line[23];
+  int month = (mchar > 64) ? mchar - 'A' + 10 : mchar - '0';
+  char dchar = line[24];
+  int day = (dchar > 64) ? dchar - 'A' + 10 : dchar - '0';
+
+  double jd = day - 32075 + (int)(1461*(year + 4800 + (int)((month-14)/12))/4)
+    + (int)(367*(month - 2 - (int)((month - 14)/12)*12)/12)
+    - (int)(3*(int)((year + 4900 + (int)((month - 14)/12))/100)/4) - 0.5;
+
+  planet->epoch = jd - 2400000.5;;
+  return TRUE;
+}
+
+// read the full MPCORB file, storing desired information in Planets structure
+Planets *mpcorb_read_text (char *filename, int *nplanets) {
+
+  ALLOCATE_PTR (buffer, char, NBYTE);
+
+  // scan through the entire MPCORB file
+  FILE *f = fopen (filename, "r");
+  if (f == NULL) Shutdown ("can't open data file: %s", filename);
+
+  int Nbyte  = 0;  // number of bytes read on each pass
+  int Nextra = 0;  // number excess bytes from last partial row
+  int Nline  = 0;  // number of lines read so far in file
+
+  int Nplanets = 0;
+  int NPLANETS = 10000;
+  ALLOCATE_PTR (planets, Planets, NPLANETS);
+
+  while ((Nbyte = fread (&buffer[Nextra], 1, NBYTE-Nextra, f)) != 0) {
+    if (Nbyte == -1) Shutdown ("error reading from data file %s", filename);
+    if (DEBUG) fprintf (stderr, "read %d bytes", Nbyte);
+
+    Nbyte += Nextra;
+
+    if (VERBOSE) fprintf (stderr, "read .. ");
+
+    /* find bounds on first complete line */
+    char *p = buffer;
+    char *q = memchr (p, '\n', Nbyte);
+    if (q == NULL) Shutdown ("incomplete line at end of file\n");
+    int offset = p - buffer; // offset within this scan
+
+    // scan through entire buffer for star coords
+    while (1) {
+
+      Nline ++;
+      if ((Nline > 43) && (q - p > 5)) {
+	mpcorb_parseline (p, &planets[Nplanets], Nbyte - offset);
+	Nplanets ++;
+	CHECK_REALLOCATE (planets, Planets, NPLANETS, Nplanets, 10000);
+      }
+
+      /* start of the next line */
+      p = q + 1;
+      offset = p - buffer; // offset within this scan
+      if (offset == Nbyte) {
+	// last buffer read is a complete line
+	Nextra = 0;
+	break;
+      }
+
+      /* end of the next line */
+      q = memchr (p, '\n', Nbyte - offset);
+      if (q == NULL) {
+	// last, incomplete line in buffer
+	Nextra = Nbyte - offset;
+	break;
+      } 
+    }
+
+    if (VERBOSE) fprintf (stderr, "scan %d planets (%s).. ", Nplanets, planets[Nplanets-1].ID);
+
+    // at end, p points at the start of last, partial line
+    if (Nextra) memmove (buffer, p, Nextra);
+
+    // TESTING:
+    if (TESTING && (Nplanets > 50)) goto skip;
+
+    // XX temporary:
+    // if (Nplanets > 100000) goto skip;
+
+  }
+  if (VERBOSE) fprintf (stderr, "\n");
+  
+skip:
+
+  fclose (f);
+  free (buffer);
+
+  *nplanets = Nplanets;
+  return planets;
+}
+
+// print an error message and exit
+int Shutdown (char *format, ...) {  
+  va_list argp;
+  char formatplus[1024];
+  
+  strcpy (formatplus, format);
+  strcat (formatplus, "\n");
+
+  va_start (argp, format);
+  vfprintf (stderr, formatplus, argp);
+  va_end (argp);
+
+  fprintf (stderr, "ERROR: mpcorb_convert halted\n");
+  exit (1);
+}
+
+// The following test case provided by Serge can be perfectly reproduced, but only if I
+// use Dtt for the epoch, and this is clearly wrong.
+void mpcorb_predict_test_from_serge () {
+
+  double epochUT = 57435.5;
+
+  Planets testPlanet;
+  // testPlanet.epoch        = epochUT + slaDtt(epochUT) / 86400.0; // XXX this is what the epoch should be
+  testPlanet.epoch        = slaDtt(epochUT); // XXX this is the value that works with the test values
+  testPlanet.inclination  = 0.043306854729735354 * DEG_RAD;
+  testPlanet.ascend_node  = 4.728915154005972    * DEG_RAD;
+  testPlanet.perihelion   = 2.2756608097151414   * DEG_RAD;
+  testPlanet.mean_anomaly = 2.888163341284418    * DEG_RAD;
+  testPlanet.eccentricity = 0.09528339999999975;
+  testPlanet.semimajor_a  = 2.782049599999998;
+		
+  double Robs, Dobs;
+  mpcorb_predict (&testPlanet, epochUT, &Robs, &Dobs, TRUE);
+
+  fprintf (stderr, "%f vs %f : %f\n",  3.92105439135726600*DEG_RAD, Robs, 3600*(3.92105439135726600*DEG_RAD - Robs));
+  fprintf (stderr, "%f vs %f : %f\n", -0.33499850519808244*DEG_RAD, Dobs, 3600*(-0.33499850519808244*DEG_RAD - Dobs));
+  exit (0);
+}
+
+// The following test code uses unperturbed orbital elements and gets the
+// wrong answer by 20 - 60 arcsec for the top 10 minor planets.
+void mpcorb_predict_t0 (Planets *planet, double mjdObs, double *Robs, double *Dobs) {
+
+  int jstat;
+  
+  double dTT = slaDtt (mjdObs);
+  double mjdObsTT = mjdObs + dTT / 86400.0;
+
+  // this function returns the topocentric apparent R,D.  I assume it
+  // has no atm corrections.
+
+  double Rrad, Drad, dist;
+  slaPlante (mjdObsTT,
+	     0.0, // PS1_LONGITUDE*RAD_DEG,
+	     0.0, // PS1_LATITUDE*RAD_DEG,
+	     2,
+	     planet->epoch,
+	     planet->inclination*RAD_DEG,
+	     planet->ascend_node*RAD_DEG,
+	     planet->perihelion*RAD_DEG,
+	     planet->semimajor_a,
+	     planet->eccentricity,
+	     planet->mean_anomaly*RAD_DEG,
+	     0.0, &Rrad, &Drad, &dist, &jstat);
+
+  fprintf (stderr, "Plante   : %12.6f %12.6f\n", Rrad*DEG_RAD, Drad*DEG_RAD);
+
+  // I'm not certain if I should use realistic (STD) atm parameters or not.
+  // If I set them to zero, is the correction not applied?
+ 
+  double Rap, Dap;
+  slaOap ("R",
+	  Rrad,
+	  Drad,
+	  mjdObs,
+	  -0.4,
+	  0.0, // LONGITUDE
+	  0.0, // LATITUDE
+	  0.0, // ALTITUDE,
+	  0.0, // xp
+	  0.0, // yp
+	  273.15, // T(K)
+	  1013.25, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
+	  0.0, // RH(%)
+	  0.7, // Wavelength (um)
+	  0.0065, // lapse rate
+	  &Rap, &Dap);
+
+  // XXX the atm effects seem to be large and incompatible with MPC outputs
+  // fprintf (stderr, "Oap (STD): %12.6f %12.6f\n", Rap*DEG_RAD, Dap*DEG_RAD);
+
+  slaOap ("R",
+	  Rrad,
+	  Drad,
+	  mjdObs,
+	  -0.4,
+	  0.0, // LONGITUDE
+	  0.0, // LATITUDE
+	  0.0, // ALTITUDE,
+	  0.0, // xp
+	  0.0, // yp
+	  273.15, // T(K)
+	  0.0, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
+	  0.0, // RH(%)
+	  0.7, // Wavelength (um)
+	  0.0, // lapse rate
+	  &Rap, &Dap);
+
+  fprintf (stderr, "Oap (0.0): %12.6f %12.6f : %6.2f %6.2f\n", Rap*DEG_RAD, Dap*DEG_RAD, 3600*(Rap-Rrad)*DEG_RAD*cos(Drad), 3600*(Dap-Drad)*DEG_RAD);
+
+  double Rmean, Dmean;
+  slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
+
+  fprintf (stderr, "Amp1     : %12.6f %12.6f : %6.2f %6.2f\n", Rmean*DEG_RAD, Dmean*DEG_RAD, 3600*(Rmean-Rap)*DEG_RAD*cos(Dap), 3600*(Dmean-Dap)*DEG_RAD);
+
+  // apply Amp to output from Plante
+  double Ralt, Dalt;
+  slaAmp (Rrad, Drad, mjdObsTT, 2000.0, &Ralt, &Dalt);
+
+  fprintf (stderr, "Amp2     : %12.6f %12.6f : %6.2f %6.2f\n", Ralt*DEG_RAD, Dalt*DEG_RAD, 3600*(Ralt-Rrad)*DEG_RAD*cos(Drad), 3600*(Dalt-Drad)*DEG_RAD);
+
+
+  // slaAmp applies both precession/nutation as well as abberation
+  // but, I think abberation is not appropriate to solar system bodies.
+
+  // here, I instead calculate just the precession and nutation with slalib
+
+  double Rpre, Dpre;
+  Rpre = Rap;
+  Dpre = Dap;
+  slaPreces ("FK5", 2021.5, 2000.000000, &Rpre, &Dpre);
+  fprintf (stderr, "Prec     : %12.6f %12.6f : %6.2f %6.2f\n", Rpre*DEG_RAD, Dpre*DEG_RAD, 3600*(Rpre-Rap)*DEG_RAD*cos(Dap), 3600*(Dpre-Dap)*DEG_RAD);
+
+  *Robs = Rmean * DEG_RAD;
+  *Dobs = Dmean * DEG_RAD;
+}
+