Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile	(revision 40661)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile	(revision 40662)
@@ -182,4 +182,6 @@
 $(SRC)/vsmooth.$(ARCH).o	   \
 $(SRC)/vstats.$(ARCH).o		   \
+$(SRC)/virls.$(ARCH).o		   \
+$(SRC)/vwtmean.$(ARCH).o \
 $(SRC)/xsection.$(ARCH).o          \
 $(SRC)/vsh.$(ARCH).o	           \
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c	(revision 40661)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c	(revision 40662)
@@ -163,4 +163,7 @@
 int vzload           PROTO((int, char **));
 int vstats           PROTO((int, char **));
+int vstats           PROTO((int, char **));
+int virls            PROTO((int, char **));
+int vwtmean          PROTO((int, char **));
 int vroll            PROTO((int, char **));
 int vshift           PROTO((int, char **));
@@ -360,4 +363,6 @@
   {1, "vsmooth",      vsmooth,          "Gaussian smooth of a vector"},
   {1, "vstats",       vstats,           "statistics on a vector"},
+  {1, "vwtmean",      vwtmean,          "weighted mean of a vector"},
+  {1, "virls",        virls,            "IRLS mean of a vector"},
   {1, "vzload",       vzload,           "load vectors as overlay on image display (scaled points)"},
   {1, "vsh",          vsh,              "Vector Spherical Harmonics"},
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c	(revision 40661)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/medimage_commands.c	(revision 40662)
@@ -1,3 +1,8 @@
 # include "data.h"
+
+float weight_cauchy_square_flt (float x2);
+float irls_mean (float *val, float *wgt, int N);
+
+# define IRLS_TOLERANCE 1e-4
 
 int medimage_list (int argc, char **argv) {
@@ -15,15 +20,16 @@
 int medimage_add (int argc, char **argv) {
 
+  int N;
   Buffer *image;
-
-  Buffer *var = NULL;
-  if ((N = get_argument (argc, argv, "-var"))) {
-    remove_argument (N, &argc, argv);
-    if ((var = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  Buffer *sig = NULL;
+
+  if ((N = get_argument (argc, argv, "-sigma"))) {
+    remove_argument (N, &argc, argv);
+    if ((sig = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     remove_argument (N, &argc, argv);
   }
 
   if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-var var]\n");
+    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-sigma sigma]\n");
     gprint (GP_ERR, "       add the given image to the set of images to be medianed\n");
     gprint (GP_ERR, "       optionally supply variance image (for weighted calculations)\n");
@@ -33,5 +39,11 @@
   if ((image = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
 
-XXX: match dimensions wit var image
+  if (sig) {
+    if ((sig->matrix.Naxis[0] != image->matrix.Naxis[0]) ||
+	(sig->matrix.Naxis[1] != image->matrix.Naxis[1])) {
+      gprint (GP_ERR, "sigma buffer does not match image buffer dimensions\n");
+      return FALSE;
+    }
+  }
 
   MedImageType *median = FindMedImage (argv[1]);
@@ -56,34 +68,46 @@
   // new image should match existing medimage dimensions
 
+  // AddMedImage (median, image, sig);
   int Ninput = median->Ninput;
   median->Ninput ++;
+  REALLOCATE (median->flx, float *, median->Ninput);
   REALLOCATE (median->sig, float *, median->Ninput);
-  REALLOCATE (median->var, float *, median->Ninput);
-
-  ALLOCATE (median->sig[Ninput], float, median->Nx*median->Ny);
-  memcpy (median->sig[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
-
-  median->var[Ninput] = NULL;
-  if (var) {
-    ALLOCATE (median->var[Ninput], float, median->Nx*median->Ny);
-    memcpy (median->var[Ninput], var->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
-  }
-
-  return TRUE;
-}
-
+
+  ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny);
+  memcpy (median->flx[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+
+  median->sig[Ninput] = NULL;
+  if (sig) {
+    ALLOCATE (median->sig[Ninput], float, median->Nx*median->Ny);
+    memcpy (median->sig[Ninput], sig->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
+  }
+
+  return TRUE;
+}
+
+enum {CALC_MEDIAN, CALC_MEAN, CALC_IRLS, CALC_WTMEAN};
 int medimage_calc (int argc, char **argv) {
 
   int ix, iy, n, N;
   Buffer *output;
-
-  int CALC_MEAN = FALSE;
+  
+  int mode = CALC_MEDIAN;
   if ((N = get_argument (argc, argv, "-mean"))) {
-    CALC_MEAN = TRUE;
+    mode = CALC_MEAN;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-irls"))) {
+    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -irls, -wtmean\n"); return FALSE; }
+    mode = CALC_IRLS;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-wtmean"))) {
+    if (mode != CALC_MEDIAN) { gprint (GP_ERR, "supply only one of -mean, -irls, -wtmean\n"); return FALSE; }
+    mode = CALC_WTMEAN;
     remove_argument (N, &argc, argv);
   }
 
   if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean]\n");
+    gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean,-irls,-wtmean]\n");
     gprint (GP_ERR, "       calculate the median image for the median image set\n");
     return FALSE;
@@ -102,6 +126,6 @@
   int Ny = median->Ny;
 
-  float *value = NULL;
-  ALLOCATE (value, float, Ninput);
+  ALLOCATE_PTR (val, float, Ninput);
+  ALLOCATE_PTR (wgt, float, Ninput);
 
   gfits_free_matrix (&output->matrix);
@@ -117,67 +141,135 @@
       int Npix = ix + Nx*iy;
       for (n = 0; n < Ninput; n++) {
-	float v = median->sig[n][Npix];
+	float v = median->flx[n][Npix];
 	if (!isfinite(v)) continue;
-	value[N] = v;
+	val[N] = v;
+	wgt[N] = 1.0;
+	if (median->sig[n]) {
+	  float s = median->sig[n][Npix];
+	  if (!isfinite(s)) continue;
+	  if (fabs(s) < 2*FLT_MIN) s = 2*FLT_MIN;
+	  wgt[N] = 1.0 / SQ(s);
+	}
 	N++;
       }
       if (N == 0) continue;
 
-    XXXX: add in calc options:
-      weighted mean
-      irls
-
-      if (CALC_MEAN) {
-	float sum = 0.0;
-	for (n = 0; n < N; n++) {
-	  sum += value[n];
+      switch (mode) {
+	case CALC_MEDIAN:
+	  fsort (val, N);
+	  outvalue[Npix] = val[(int)(0.5*N)];
+	  break;
+	case CALC_MEAN: {
+	  float sum = 0.0;
+	  for (n = 0; n < N; n++) {
+	    sum += val[n];
+	  }
+	  outvalue[Npix] = sum / (float) N;
+	  break;
 	}
-	outvalue[Npix] = sum / (float) N;
-      } else {
-	fsort (value, N);
-	outvalue[Npix] = value[(int)(0.5*N)];
+	case CALC_WTMEAN: {
+	  float S1 = 0.0, S2 = 0.0;
+	  for (n = 0; n < N; n++) {
+	    S1 += wgt[n] * val[n];
+	    S2 += wgt[n];
+	  }
+	  outvalue[Npix] = S1 / S2;
+	  break;
+	}
+	case CALC_IRLS: 
+	  outvalue[Npix] = irls_mean (val, wgt, N);
       }
     }
   }
-
-  return TRUE;
-}
+  return TRUE;
+}
+
+float irls_mean (float *val, float *wgt, int N) {
+
+  // calculate weighted mean
+  float S1 = 0.0, S2 = 0.0;
+  for (int n = 0; n < N; n++) {
+    S1 += wgt[n] * val[n];
+    S2 += wgt[n];
+  }
+  float Value = S1 / S2;
+  
+  int converged = FALSE;
+  for (int i = 0; (i < 10) && !converged; i++) {
+    float ValueLast = Value;
+
+    float S1 = 0.0, S2 = 0.0;
+    
+    // calculate weight modification based on distances (squared).
+    // use modifier to calculate new weighted mean
+    for (int n = 0; n < N; n++) {
+      float dV = (val[n] - Value);
+      float d2 = SQ(dV) * wgt[n];
+      
+      float Mod = weight_cauchy_square_flt (d2);
+      S1 += Mod * wgt[n] * val[n];
+      S2 += Mod * wgt[n];
+    }
+    Value = S1 / S2;
+
+    float delta = fabs(Value - ValueLast);
+    if (delta < Value * IRLS_TOLERANCE) converged = TRUE;
+  }
+  return Value;
+}
+
+// exp(-(x^2/s^2)/2) = (1/2)
+//     -(x^2/s^2)/2  = ln(1/2)
+//      (x^2/s^2)/2  = ln(2)
+//      (x^2/s^2)    = 2ln(2)
+//      (x  /s)      = sqrt(2ln(2)) : half-width at half-max
+//       FWHM        = 2sqrt(2ln(2))
+
+// R2 = (X / 2.385)^2 = (X^2 / 2.385^2)
+
+# define CAUCY_FACTOR 1.0
+
+float weight_cauchy_square_flt (float x2) {
+  float r2 = x2 / CAUCY_FACTOR;
+  return (1.0 / (1.0 + r2));
+}
+
 
 /* 
-int medimage_save (int argc, char **argv) {
-
-  int N;
-
-  int APPEND = FALSE;
-  if ((N = get_argument (argc, argv, "-append"))) {
-    APPEND = TRUE;
-    remove_argument (N, &argc, argv);
-  }
-
-  if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
-    return FALSE;
-  }
-
-  if (!SaveMedImage(argv[2], argv[1], APPEND)) {
-    gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
-    return (FALSE);
-  }
-  return TRUE;
-}
-
-int medimage_load (int argc, char **argv) {
-
-  if (argc != 3) {
-    gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
-    return FALSE;
-  }
-
-  if (!LoadMedImage(argv[2], argv[1])) {
-    gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
-    return (FALSE);
-  }
-  return TRUE;
-}
+   int medimage_save (int argc, char **argv) {
+
+   int N;
+
+   int APPEND = FALSE;
+   if ((N = get_argument (argc, argv, "-append"))) {
+   APPEND = TRUE;
+   remove_argument (N, &argc, argv);
+   }
+
+   if (argc != 3) {
+   gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
+   return FALSE;
+   }
+
+   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
+   gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
+   return (FALSE);
+   }
+   return TRUE;
+   }
+
+   int medimage_load (int argc, char **argv) {
+
+   if (argc != 3) {
+   gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
+   return FALSE;
+   }
+
+   if (!LoadMedImage(argv[2], argv[1])) {
+   gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
+   return (FALSE);
+   }
+   return TRUE;
+   }
 */
 
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/spline.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/spline.c	(revision 40661)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/spline.c	(revision 40662)
@@ -7,4 +7,5 @@
 int spline_load (int argc, char **argv);
 int spline_save (int argc, char **argv);
+int spline_print (int argc, char **argv);
 int spline_delete (int argc, char **argv);
 int spline_rename (int argc, char **argv);
@@ -18,4 +19,5 @@
   {1, "load",       spline_load,       "write a spline to a FITS file"},
   {1, "save",       spline_save,       "read a spline from a FITS file"},
+  {1, "print",      spline_print,      "print a spline to stdout"},
   {1, "delete",     spline_delete,     "delete a spline"},
   {1, "rename",     spline_rename,     "rename a spline"},
@@ -35,4 +37,5 @@
   gprint (GP_ERR, "    spline load   (spline) (filename)            : load a spline from a FITS file\n");
   gprint (GP_ERR, "    spline save   (spline) (filename) [-append]  : save a spline in FITS format\n");
+  gprint (GP_ERR, "    spline print  (spline)                       : print a spline to stdout\n");
 
   return FALSE;
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/spline_commands.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/spline_commands.c	(revision 40661)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/spline_commands.c	(revision 40662)
@@ -93,4 +93,26 @@
 }
 
+int spline_print (int argc, char **argv) {
+
+  int i;
+
+  if (argc != 2) {
+    gprint (GP_ERR, "USAGE: spline print (name)\n");
+    return FALSE;
+  }
+
+  Spline *myspline = FindSpline (argv[1]);
+  if (!myspline) {
+    gprint (GP_ERR, "spline %s not found\n", argv[1]);
+    return (FALSE);
+  }
+
+  for (i = 0; i < myspline->Nknots; i++) {
+    gprint (GP_LOG, "%lf : %lf : %lf\n", myspline->xk[i], myspline->yk[i], myspline->y2[i]);
+  }    
+
+  return TRUE;
+}
+
 int spline_load (int argc, char **argv) {
 
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/virls.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/virls.c	(revision 40662)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/virls.c	(revision 40662)
@@ -0,0 +1,103 @@
+# include "data.h"
+
+double opihi_irls_mean_dbl (double *val, double *wgt, int N, int *Converged);
+double weight_cauchy_square_dbl (double x2);
+
+# define IRLS_TOLERANCE 1e-4
+
+int virls (int argc, char **argv) {
+  
+  if (argc != 3) {
+    gprint (GP_ERR, "USAGE: virls (values) (sigmas)\n");
+    gprint (GP_ERR, "NOTE: sigmas with abs values < %e will be truncated to that value\n", FLT_MIN);
+    return (FALSE);
+  }
+
+  Vector *value, *sigma;
+
+  if ((value = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((sigma = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  if (value->Nelements != sigma->Nelements) {
+    gprint (GP_ERR, "vectors are not the same length\n");
+    return FALSE;
+  }
+  if (value->type != OPIHI_FLT) {
+    gprint (GP_ERR, "%s is not a floating point vector\n", argv[1]);
+    return FALSE;
+  }
+  if (sigma->type != OPIHI_FLT) {
+    gprint (GP_ERR, "%s is not a floating point vector\n", argv[2]);
+    return FALSE;
+  }
+
+  ALLOCATE_PTR (val, opihi_flt, value->Nelements);
+  ALLOCATE_PTR (wgt, opihi_flt, value->Nelements);
+
+  int N = 0;
+  for (int i = 0; i < value->Nelements; i++) {
+    if (!isfinite(value->elements.Flt[i])) continue;
+
+    double s = sigma->elements.Flt[i];
+    if (!isfinite(s)) continue;
+
+    if (fabs(s) < FLT_MIN) s = FLT_MIN;
+
+    val[N] = value->elements.Flt[i];
+    wgt[N] = 1.0 / SQ(s);
+    N ++;
+  }
+
+  int converged = TRUE;
+  double Value = opihi_irls_mean_dbl (val, wgt, N, &converged);
+
+  set_variable ("WTMEAN",     Value);
+  set_int_variable ("NUSED", N);
+  set_int_variable ("CONVERGED", converged);
+  return (TRUE);
+}
+
+double opihi_irls_mean_dbl (double *val, double *wgt, int N, int *Converged) {
+
+  // calculate weighted mean
+  double S1 = 0.0, S2 = 0.0;
+  for (int n = 0; n < N; n++) {
+    S1 += wgt[n] * val[n];
+    S2 += wgt[n];
+  }
+  double Value = S1 / S2;
+  
+  int converged = FALSE;
+  for (int i = 0; (i < 10) && !converged; i++) {
+    double ValueLast = Value;
+
+    double S1 = 0.0, S2 = 0.0;
+
+    // calculate weight modification based on distances (squared).
+    // use modifier to calculate new weighted mean
+    for (int n = 0; n < N; n++) {
+      double dV = (val[n] - Value);
+      double d2 = SQ(dV) * wgt[n];
+      
+      double Mod = weight_cauchy_square_dbl (d2);
+      S1 += Mod * wgt[n] * val[n];
+      S2 += Mod * wgt[n];
+      fprintf (stderr, "%f %f : %f %f : %f\n", val[n], wgt[n], dV, sqrt(d2), Mod);
+    }
+    Value = S1 / S2;
+
+    double delta = fabs(Value - ValueLast);
+    if (delta < Value * IRLS_TOLERANCE) converged = TRUE;
+
+    // XXX if the answer is close to zero, we might not converge
+  }
+  *Converged = converged;
+  return Value;
+}
+
+# define CAUCY_FACTOR 1.0
+
+double weight_cauchy_square_dbl (double x2) {
+  double r2 = x2 / CAUCY_FACTOR;
+  return (1.0 / (1.0 + r2));
+}
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vwtmean.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vwtmean.c	(revision 40662)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vwtmean.c	(revision 40662)
@@ -0,0 +1,70 @@
+# include "data.h"
+
+double opihi_wt_mean_dbl (double *val, double *wgt, int N, double *Sigma);
+
+int vwtmean (int argc, char **argv) {
+  
+  if (argc != 3) {
+    gprint (GP_ERR, "USAGE: vwtmean (values) (sigmas)\n");
+    gprint (GP_ERR, "NOTE: sigmas with abs values < %e will be truncated to that value\n", FLT_MIN);
+    return (FALSE);
+  }
+
+  Vector *value, *sigma;
+
+  if ((value = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((sigma = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  if (value->Nelements != sigma->Nelements) {
+    gprint (GP_ERR, "vectors are not the same length\n");
+    return FALSE;
+  }
+  if (value->type != OPIHI_FLT) {
+    gprint (GP_ERR, "%s is not a floating point vector\n", argv[1]);
+    return FALSE;
+  }
+  if (sigma->type != OPIHI_FLT) {
+    gprint (GP_ERR, "%s is not a floating point vector\n", argv[2]);
+    return FALSE;
+  }
+
+  ALLOCATE_PTR (val, opihi_flt, value->Nelements);
+  ALLOCATE_PTR (wgt, opihi_flt, value->Nelements);
+
+  int N = 0;
+  for (int i = 0; i < value->Nelements; i++) {
+    if (!isfinite(value->elements.Flt[i])) continue;
+
+    double s = sigma->elements.Flt[i];
+    if (!isfinite(s)) continue;
+
+    if (fabs(s) < FLT_MIN) s = FLT_MIN;
+
+    val[N] = value->elements.Flt[i];
+    wgt[N] = 1.0 / SQ(s);
+    N ++;
+  }
+
+  double Sigma = NAN;
+  double Value = opihi_wt_mean_dbl (val, wgt, N, &Sigma);
+
+  set_variable ("WTMEAN",     Value);
+  set_variable ("WTSIGMA",    Sigma);
+  set_int_variable ("NUSED", N);
+  return (TRUE);
+}
+
+double opihi_wt_mean_dbl (double *val, double *wgt, int N, double *Sigma) {
+
+  // calculate weighted mean
+  double S1 = 0.0, S2 = 0.0;
+  for (int n = 0; n < N; n++) {
+    S1 += wgt[n] * val[n];
+    S2 += wgt[n];
+  }
+  double Value = S1 / S2;
+  *Sigma = 1.0 / S2;
+  
+  return Value;
+}
+
Index: /branches/eam_branches/ohana.20190329/src/opihi/include/data.h
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/include/data.h	(revision 40661)
+++ /branches/eam_branches/ohana.20190329/src/opihi/include/data.h	(revision 40662)
@@ -48,6 +48,6 @@
   int Nx;
   int Ny;
+  float **flx;
   float **sig;
-  float **var;
 } MedImageType;
 
Index: /branches/eam_branches/ohana.20190329/src/opihi/lib.data/MedImageOps.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/lib.data/MedImageOps.c	(revision 40661)
+++ /branches/eam_branches/ohana.20190329/src/opihi/lib.data/MedImageOps.c	(revision 40662)
@@ -32,6 +32,8 @@
   free (medimage[0].name);
   for (i = 0; i < medimage[0].Ninput; i++) {
-    free (medimage[0].sig[i]);
+    free (medimage[0].flx[i]);
+    FREE (medimage[0].sig[i]);
   }
+  free (medimage[0].flx);
   free (medimage[0].sig);
   free (medimage);
@@ -75,4 +77,5 @@
   medimage->Nx = Nx;
   medimage->Ny = Ny;
+  ALLOCATE (medimage->flx, float *, 1);
   ALLOCATE (medimage->sig, float *, 1);
 
Index: /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c	(revision 40661)
+++ /branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkobj.c	(revision 40662)
@@ -190,8 +190,19 @@
     if (psf) Npof -= Npsf;
 
+    // integral and fractional pixel offsets of PSF (none, by default)
+    int dxi = 0;
+    float dxf = 0.0, dxr = 1.0;
+    
     if (psf_trace) {
-      float dx = -spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
-      Npof += dx;
-    }
+      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
+      dxi = floor(dx);
+      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
+      dxr = 1 - dxf;
+      Npof += dxi;
+      if (iy % 200 == 0) { fprintf (stderr, "%d : %f : %d : %f %f\n", iy, dx, dxi, dxf, dxr); }
+    }
+
+    // if fractional offset is small, do not interpolate
+    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
 
     // flux = obj * PSF + sky
@@ -205,5 +216,14 @@
 	// only add in the flux if we are in range of the PSF
 	if ((n >= 0) && (n < NpsfFull)) {
-	  value += objVy * psfV[n];
+	  if (doInterp) {
+	    if ((n > 0) && (n < NpsfFull)) {
+	      value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
+	    }
+	    if (n == 0) {
+	      value += objVy*psfV[n]*dxr;
+	    }
+	  } else {
+	    value += objVy*psfV[n];
+	  }
 	}
       } else {