Index: /branches/eam_branches/ipp-20140904/Ohana/src/libautocode/def/AstromOffsetMap_Disk.d
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libautocode/def/AstromOffsetMap_Disk.d	(revision 37387)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libautocode/def/AstromOffsetMap_Disk.d	(revision 37387)
@@ -0,0 +1,26 @@
+STRUCT  AstromOffsetMap_Disk_6x6
+EXTNAME ASTROM_OFFSET_MAP_DISK_6x6
+TYPE    BINTABLE
+SIZE    312
+
+# we have one row per correction image, with the max correction dimensions hard-wired in the structure
+
+# ID is the ID of this correction image, image ID is the image being corrected
+FIELD   ID,        ID,             unsigned int
+FIELD   imageID,   IMAGE_ID,       unsigned int
+
+# 6x6 is the max size; Nx,Ny give the actual size of the array
+FIELD   Nx,	   NX,             int
+FIELD   Ny,	   NY,             int
+
+# for an image of size (NxBig,NyBig) and a map of size (Nx,Ny), the relationship between
+# real and map pixels is given by
+
+# Ix = ix * (Nx / NxBig)
+# dX = (Nx / NxBig)
+FIELD   dX,	   DX_SCALE,       float
+FIELD   dY,	   DY_SCALE,       float
+
+FIELD   dXv,       DX_VALUE,       float[6][6]
+FIELD   dYv,       DY_VALUE,       float[6][6]
+
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapIO.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapIO.c	(revision 37387)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapIO.c	(revision 37387)
@@ -0,0 +1,187 @@
+# include "dvo.h"
+
+/* AstromOffsetMap is the map of Astrometry Offsets for non-polynomial astrometric
+ * corrections for each chip.  
+ */
+
+AstromOffsetTable *AstromOffsetMapLoad (char *filename, int VERBOSE) {
+
+  off_t Nmap;
+
+  Header header;
+  Matrix matrix;
+  Header theader;
+  FTable ftable;
+
+  FILE *f = fopen (filename, "r");
+  if (f == NULL) {
+    if (VERBOSE) fprintf (stderr, "can't find Astrom Offset Map file %s\n", filename);
+    return (NULL);
+  }
+
+  /* load in table data */
+  ftable.header = &theader;
+  if (!gfits_fread_header (f, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map header\n");
+    fclose (f);
+    return (NULL);
+  }
+  if (!gfits_fread_matrix (f, &matrix, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map matrix\n");
+    gfits_free_header (&header);
+    fclose (f);
+    return (NULL);
+  }
+
+  // for now, we only have one flavor (6x6)
+
+  if (!gfits_fread_ftable (f, &ftable, "ASTROM_OFFSET_MAP_DISK_6x6")) {
+    if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map table\n");
+    gfits_free_header (&header);
+    gfits_free_matrix (&matrix);
+    fclose (f);
+    return (NULL);
+  }
+  AstromOffsetMap_Disk_6x6 *map_disk = gfits_table_get_AstromOffsetMap_Disk_6x6 (&ftable, &Nmap, NULL);
+  if (!map_disk) {
+    fprintf (stderr, "ERROR: failed to read Astrom Offset Map\n");
+    exit (2);
+  }
+
+  // convert the disk table format to internal format:
+  AstromOffsetTable *table = AstromOffsetMapToTable (map_disk, Nmap);
+  free (map_disk);
+
+  gfits_free_header (&theader);
+  gfits_free_matrix (&matrix);
+  gfits_free_table  (&ftable);
+
+  return (table);
+}
+
+int AstromOffsetMapSave (AstromOffsetTable *table, char *filename) {
+
+  off_t Nmap;
+
+  Header header;
+  Matrix matrix;
+  Header theader;
+  FTable ftable;
+
+  FILE *f;
+
+  /* make phu header (no matrix needed) */
+  gfits_init_header (&header);
+  header.extend = TRUE;
+  gfits_create_header (&header);
+  gfits_create_matrix (&header, &matrix);
+
+  AstromOffsetMap_Disk_6x6 *map_disk = AstromOffsetTableToMap (table, &Nmap); 
+
+  ftable.header = &theader;
+  gfits_table_set_AstromOffsetMap_Disk_6x6 (&ftable, map_disk, Nmap);
+
+  f = fopen (filename, "w");
+  if (f == (FILE *) NULL) { 
+    fprintf (stderr, "cannot open %s for output\n", filename);
+    return (FALSE);
+  }
+  
+  gfits_fwrite_header  (f, &header);
+  gfits_fwrite_matrix  (f, &matrix);
+  gfits_fwrite_Theader (f, &theader);
+  gfits_fwrite_table  (f, &ftable);
+  fclose (f);
+
+  return (TRUE);
+}
+
+AstromOffsetTable *AstromOffsetMapToTable(AstromOffsetMap_Disk_6x6 *map_disk, off_t Nmap) {
+
+  int i, j, k;
+
+  AstromOffsetTable *table = NULL;
+  ALLOCATE (table, AstromOffsetTable, 1);
+
+  table->Nmap = Nmap;
+  ALLOCATE (table->map, AstromOffsetMap, Nmap);
+
+  // find the max value of ID
+  int MaxID = 0;
+  for (i = 0; i < Nmap; i++) {
+    MaxID = MAX(map_disk[i].ID, MaxID);
+  }
+  MaxID ++; // we want the outer bound, not the last value
+
+  // generate the index and init values to -1
+  ALLOCATE (table->IDtoSeq, int, MaxID);
+  for (i = 0; i < MaxID; i++) {
+    table->IDtoSeq[i] = -1;
+  }
+  table->MaxID = MaxID;
+
+  // assign the ID values
+  for (i = 0; i < Nmap; i++) {
+    int ID = map_disk[i].ID;
+    myAssert (table->IDtoSeq[ID] == -1, "oops, duplicate map IDs");
+    table->IDtoSeq[ID] = i;
+  }
+
+  // assign the map values (this allocates just the area needed for each image, not the
+  // full 6x6, saving some memory while doing the analysis)
+  for (i = 0; i < Nmap; i++) {
+    table->map[i].Nx 	  = map_disk[i].Nx;
+    table->map[i].Ny 	  = map_disk[i].Ny;
+    table->map[i].ID 	  = map_disk[i].ID;
+    table->map[i].imageID = map_disk[i].imageID;
+    
+    // STORE THESE VALUES?
+    table->map[i].dY 	  = map_disk[i].dX;
+    table->map[i].dX 	  = map_disk[i].dY;
+
+    ALLOCATE (table->map[i].dXv, float *, map_disk[i].Nx);
+    ALLOCATE (table->map[i].dYv, float *, map_disk[i].Nx);
+
+    for (j = 0; j < map_disk[i].Nx; j++) {
+      ALLOCATE (table->map[i].dXv[j], float, map_disk[i].Ny);
+      ALLOCATE (table->map[i].dYv[j], float, map_disk[i].Ny);
+
+      for (k = 0; k < map_disk[i].Ny; k++) {
+	table->map[i].dXv[j][k] = map_disk[i].dXv[j][k];
+	table->map[i].dYv[j][k] = map_disk[i].dYv[j][k];
+      }
+    }
+  }
+
+  return table;
+}
+
+AstromOffsetMap_Disk_6x6 *AstromOffsetTableToMap(AstromOffsetTable *table, off_t *Nmap) {
+
+  int i, j, k;
+
+  AstromOffsetMap_Disk_6x6 *map_disk = NULL;
+  ALLOCATE (map_disk, AstromOffsetMap_Disk_6x6, Nmap);
+
+  // assign the map values (this allocates just the area needed for each image, not the
+  // full 6x6, saving some memory while doing the analysis)
+  for (i = 0; i < table->Nmap; i++) {
+    map_disk[i].Nx 	 = table->map[i].Nx;
+    map_disk[i].Ny 	 = table->map[i].Ny;
+    map_disk[i].ID 	 = table->map[i].ID;
+    map_disk[i].imageID  = table->map[i].imageID;
+    
+    map_disk[i].dX 	 = table->map[i].dX;
+    map_disk[i].dY 	 = table->map[i].dY;
+
+    for (j = 0; j < map_disk[i].Nx; j++) {
+      for (k = 0; k < map_disk[i].Ny; k++) {
+	map_disk[i].dXv[j][k] = table->map[i].dXv[j][k];
+	map_disk[i].dYv[j][k] = table->map[i].dYv[j][k];
+      }
+    }
+  }
+  *Nmap = table->Nmap;
+  return map_disk;
+}
+
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c	(revision 37387)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c	(revision 37387)
@@ -0,0 +1,328 @@
+# include "dvo.h"
+
+/* AstromOffsetMap functions:
+ */
+
+float AstromOffsetMapValue (AstromOffsetMap *map, float x, float y, int xdir) {
+
+  // given x,y find the offset (x-direction) at that location from the map
+
+  float **V = xdir ? map->dXv : map->dYv;
+
+  // if there is no spatial information, the value at the point is the value in the map
+  if ((map->Nx == 1) && (map->Ny == 1)) return V[0][0];
+
+  // if there is no spatial information in 1D, the interpolation is 1D
+  if (map->Nx == 1) {
+    float ymo = y * map->dY - 0.5;
+    int   ymi = floor(ymo);
+    ymi = MAX(0,MIN(map->Ny - 1, ymi)); // force range of ymi to be 0,Ny-1
+
+    float ymf = ymo - ymi;
+
+    float value = ymf * V[0][ymi+1] + (1.0 - ymf) * V[0][ymi];
+    return value;
+  }
+
+  // if there is no spatial information in 1D, the interpolation is 1D
+  if (map->Ny == 1) {
+    float xmo = x * map->dX - 0.5;
+    int   xmi = floor(xmo);
+    xmi = MAX(0,MIN(map->Nx - 1, xmi)); // force range of ymi to be 0,Ny-1
+
+    float xmf = xmo - xmi;
+
+    float value = xmf * V[xmi+1][0] + (1.0 - xmf) * V[xmi][0];
+    return value;
+  }
+
+  // x & y are in Big Image coordinates, convert to fractional map coordinates 
+
+  float xmo = x * map->dX - 0.5;
+  int   xmi = floor(xmo);
+  xmi = MAX(0,MIN(map->Nx - 1, xmi)); // force range of ymi to be 0,Ny-1
+  float xmf = xmo - xmi;
+
+  float ymo = y * map->dY - 0.5;
+  int   ymi = floor(ymo);
+  ymi = MAX(0,MIN(map->Ny - 1, ymi)); // force range of ymi to be 0,Ny-1
+  float ymf = ymo - ymi;
+
+  float V00 = V[xmi+0][ymi+0];
+  float V01 = V[xmi+0][ymi+1];
+  float V10 = V[xmi+1][ymi+0];
+  float V11 = V[xmi+1][ymi+1];
+
+  float Vx0 = V10*xmf + V00*(1.0 - xmf);
+  float Vx1 = V11*xmf + V01*(1.0 - xmf);
+
+  float value = ymf * Vx1 + (1.0 - ymf) * Vx0;
+  return value;
+}
+
+// given (x,y),value vector sets, choose the map that minimizes the difference between 
+// the values and bilinear interpolation of the map
+int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) {
+
+  int i, ix, iy, jx, jy;
+
+  // choose to map direction:
+  float **V = xdir ? map->dXv : map->dYv;
+
+  // special cases:
+  if ((map->Nx == 1) && (map->Ny == 1)) {
+
+    // no spatial information, just return (clipped?) median or mean or something
+    VStatsType stats;
+    stats.statmode = VSTATS_INNER_MEAN;
+
+    vstats_getstats_f (f, NULL, NULL, Npts, &stats);
+    V[0][0] = stats.mean;
+    return TRUE;
+  }
+
+  // special cases:
+  if (map->Nx == 1) {
+    fprintf (stderr, "1 x N not yet coded\n");
+    exit (1);
+    return TRUE;
+  }
+
+  if (map->Ny == 1) {
+    fprintf (stderr, "N x 1 not yet coded\n");
+    exit (1);
+    return TRUE;
+  }
+
+  // this bit is copied nearly intact from psImageMapFit.c
+
+  // set up the redirection table so we can use sA[-1][-1], etc
+  float SAm[3][3], *SAv[3], **sA;
+
+  for (i = 0; i < 3; i++) {
+    SAv[i] = SAm[i] + 1;
+  }
+  sA = SAv + 1;
+
+  int Nx = map->Nx;
+  int Ny = map->Ny;
+
+  float **A, **B;
+  ALLOCATE (A, float *, Nx*Ny);
+  ALLOCATE (B, float *, Nx*Ny);
+  for (i = 0; i < Nx*Ny; i++) {
+    ALLOCATE (A[i], float, Nx*Ny);
+    ALLOCATE (B[i], float, 1);
+  }    
+
+  // we are looping over the Nx,Ny image map elements;
+  // the matrix equation contains Nx*Ny rows and columns
+  // for (int n = 1; n < Nx - 1; n++) {
+  // for (int m = 1; m < Ny - 1; m++) {
+
+  // float Total = 0.0;
+  for (ix = 0; ix < Nx; ix++) {
+    for (iy = 0; iy < Ny; iy++) {
+      // define & init summing variables
+      float rx_rx_ry_ry = 0;
+      float rx_rx_dy_ry = 0;
+      float dx_rx_ry_ry = 0;
+      float dx_rx_dy_ry = 0;
+      float fi_rx_ry    = 0;
+      float rx_rx_py_py = 0;
+      float rx_rx_qy_py = 0;
+      float dx_rx_py_py = 0;
+      float dx_rx_qy_py = 0;
+      float fi_rx_py    = 0;
+      float px_px_ry_ry = 0;
+      float px_px_dy_ry = 0;
+      float qx_px_ry_ry = 0;
+      float qx_px_dy_ry = 0;
+      float fi_px_ry    = 0;
+      float px_px_py_py = 0;
+      float px_px_qy_py = 0;
+      float qx_px_py_py = 0;
+      float qx_px_qy_py = 0;
+      float fi_px_py    = 0;
+
+      // generate the sums for the fitting matrix element I,J
+      // I = n + nX*m
+      // J = (n + jn) + nX*(m + jm)
+      for (i = 0; i < Npts; i++) {
+
+	// data value & weight for this point
+	if (!isfinite(f[i])) continue;
+
+	// if (mask && (mask[i] & maskValue)) continue;
+
+	// base coordinate offset for this point (x,y) relative to this map element (n,m)
+	// float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5);
+	// float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
+
+	float dx = x[i] * map->dX - ix - 0.5;
+	float dy = y[i] * map->dY - iy - 0.5;
+
+	// XXX do I need to do something different at the edge or not?  I want to allow
+	// x[i] to extend beyond the grid
+
+# if (0)
+	// edge cases to include:
+	bool edgeX = false;
+	edgeX |= ((n == 1) && (dx < -1.0));
+	edgeX |= ((n == Nx - 2) && (dx > +1.0));
+
+	bool edgeY = false;
+	edgeY |= ((m == 1) && (dy < -1.0));
+	edgeY |= ((m == Ny - 2) && (dy > +1.0));
+
+	// skip points outside of 2x2 grid centered on n,m:
+	if (!edgeX && (fabs(dx) > 1.0)) continue;
+	if (!edgeY && (fabs(dy) > 1.0)) continue;
+# endif
+
+	if (fabs(dx) > 1.0) continue;
+	if (fabs(dy) > 1.0) continue;
+
+	// related offset values
+	float rx = 1.0 - dx;
+	float ry = 1.0 - dy;
+	float px = 1.0 + dx;
+	float py = 1.0 + dy;
+	float qx = 0.0 - dx;
+	float qy = 0.0 - dy;
+
+	// sum the appropriate elements for the different quadrants
+	int Qx = (dx >= 0) ? 1 : 0;
+	int Qy = (dy >= 0) ? 1 : 0;
+
+	// XXX why change the selection at the edges?
+	// if (n ==      0) Qx = 1;
+	// if (n == Nx - 1) Qx = 0;
+	// if (m ==      0) Qy = 1;
+	// if (m == Ny - 1) Qy = 0;
+
+	float fi = f[i];
+	assert (isfinite(fi));
+	assert (isfinite(rx));
+	assert (isfinite(ry));
+
+	// points at offset 1,1
+	if ((Qx == 1) && (Qy == 1)) {
+	  rx_rx_ry_ry += rx*rx*ry*ry;
+	  rx_rx_dy_ry += rx*rx*dy*ry;
+	  dx_rx_ry_ry += dx*rx*ry*ry;
+	  dx_rx_dy_ry += dx*rx*dy*ry;
+	  fi_rx_ry    += fi*rx*ry;
+	}
+	// points at offset 1,0
+	if ((Qx == 1) && (Qy == 0)) {
+	  rx_rx_py_py += rx*rx*py*py;
+	  rx_rx_qy_py += rx*rx*qy*py;
+	  dx_rx_py_py += dx*rx*py*py;
+	  dx_rx_qy_py += dx*rx*qy*py;
+	  fi_rx_py    += fi*rx*py;
+	}
+	// points at offset 0,1
+	if ((Qx == 0) && (Qy == 1)) {
+	  px_px_ry_ry += px*px*ry*ry;
+	  px_px_dy_ry += px*px*dy*ry;
+	  qx_px_ry_ry += qx*px*ry*ry;
+	  qx_px_dy_ry += qx*px*dy*ry;
+	  fi_px_ry    += fi*px*ry;
+	}
+	// points at offset 0,0
+	if ((Qx == 0) && (Qy == 0)) {
+	  px_px_py_py += px*px*py*py;
+	  px_px_qy_py += px*px*qy*py;
+	  qx_px_py_py += qx*px*py*py;
+	  qx_px_qy_py += qx*px*qy*py;
+	  fi_px_py    += fi*px*py;
+	}
+      }
+
+      // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm),
+      // jn,jm = -1 to +1. Convert the sums above into the correct coefficients
+      sA[-1][-1] = qx_px_qy_py;
+      sA[-1][ 0] = qx_px_ry_ry + qx_px_py_py;
+      sA[-1][+1] = qx_px_dy_ry;
+      sA[ 0][-1] = rx_rx_qy_py + px_px_qy_py;
+      sA[ 0][ 0] = rx_rx_ry_ry + px_px_ry_ry + rx_rx_py_py + px_px_py_py;
+      sA[ 0][+1] = rx_rx_dy_ry + px_px_dy_ry;
+      sA[+1][-1] = dx_rx_qy_py;
+      sA[+1][ 0] = dx_rx_ry_ry + dx_rx_py_py;
+      sA[+1][+1] = dx_rx_dy_ry;
+
+      // I[ 0][ 0] = index for this n,m element:
+      int I = ix + Nx * iy;
+      B[I][0] = fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py;
+
+      // insert these values into their corresponding locations in A, B
+      // float Sum = 0.0;
+      for (jx = -1; jx <= +1; jx++) {
+	if (ix + jx <   0) continue;
+	if (ix + jx >= Nx) continue;
+	for (jy = -1; jy <= +1; jy++) {
+	  if (iy + jy <   0) continue;
+	  if (iy + jy >= Ny) continue;
+	  int J = (ix + jx) + Nx * (iy + jy);
+	  A[J][I] = sA[jx][jy];
+	  // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]);
+	  // Sum += sA[jn][jm];
+	}
+      }
+      // fprintf (stderr, "B %d (%d %d) : %f  :  %f\n", I, n, m, B->data.F32[I], Sum);
+      // Total += Sum;
+    }
+  }
+  // fprintf (stderr, "Total: %f\n", Total);
+
+  // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0
+  int *Empty = NULL;
+  ALLOCATE (Empty, int, Nx*Ny);
+
+  for (i = 0; i < Nx*Ny; i++) {
+    Empty[i] = 0;
+    if (A[i][i] == 0.0) {
+      Empty[i] = 1;
+      for (j = 0; j < Nx*Ny; j++) {
+	A[i][j] = 0.0;
+	A[j][i] = 0.0;
+      }
+      A[i][i] = 1.0;
+      B[i][0] = 0.0;
+    }
+  }
+
+  if (!fgaussjordan(A, B, Nx*Ny, 1)) {
+    fprintf (stderr, "FAIL\n");
+    exit (1);
+  }
+
+  // set bad values to NaN
+  for (i = 0; i < Nx*Ny; i++) {
+    if (Empty[i]) {
+      B[i][0] = NAN;
+      A[i][i] = 0;
+    }
+  }
+
+  for (ix = 0; ix < Nx; ix++) {
+    for (iy = 0; iy < Ny; iy++) {
+      int I = ix + Nx * iy;
+      V[ix][iy] = B[I][0];
+      // error[ix][iy] = sqrt(A[I][I]);
+    }
+  }
+
+  // XXX free someone?
+  for (i = 0; i < Nx*Ny; i++) {
+    free (A[i]);
+    free (B[i]);
+  }    
+  free (A);
+  free (B);
+  free (Empty);
+
+  return TRUE;
+}
+
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/vstats.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/vstats.c	(revision 37387)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/vstats.c	(revision 37387)
@@ -0,0 +1,396 @@
+# include <ohana.h>
+
+int vstats_setmode (VStatsType *stats, char *mode) {
+
+  stats->statmode = VSTATS_NONE;
+  if (!strcmp (mode, "MEAN"))             { stats->statmode = VSTATS_MEAN; return TRUE; }
+  if (!strcmp (mode, "MEDIAN"))           { stats->statmode = VSTATS_MEDIAN; return TRUE; }
+  if (!strcmp (mode, "WT_MEAN"))          { stats->statmode = VSTATS_WT_MEAN; return TRUE; }
+  if (!strcmp (mode, "INNER_MEAN"))       { stats->statmode = VSTATS_INNER_MEAN; return TRUE; }
+  if (!strcmp (mode, "INNER_WTMEAN"))     { stats->statmode = VSTATS_INNER_WTMEAN; return TRUE; }
+  if (!strcmp (mode, "CHI_INNER_MEAN"))   { stats->statmode = VSTATS_CHI_INNER_MEAN; return TRUE; }
+  if (!strcmp (mode, "CHI_INNER_WTMEAN")) { stats->statmode = VSTATS_CHI_INNER_WTMEAN; return TRUE; }
+
+  push_error ("ERROR: invalid stats mode");
+  return FALSE;
+}
+
+// NOTE: this function resorts the input vectors.  if this is not what you want, you need
+// to either pass a vector copy or save the sequence
+int vstats_getstats (double *value, double *dvalue, double *weight, int N, VStatsType *stats) {
+  
+  int i, ks, ke;
+  double Mo, dMo, M, dM, Nm, X2, dS, R, W, *chi;
+
+  myAssert (stats->statmode != VSTATS_NONE, "programming error, liststats mode not set");
+
+  ke = ks = dMo = 0;
+
+  stats[0].Nmeas   = N;
+  stats[0].mean    = NAN;
+  stats[0].median  = NAN;
+  stats[0].sigma   = NAN;
+  stats[0].error   = NAN;
+  stats[0].chisq   = NAN;
+  stats[0].min     = NAN;
+  stats[0].max     = NAN;
+  stats[0].Upper80 = NAN;
+  stats[0].Lower20 = NAN;
+  stats[0].total   = NAN;
+
+  // take care of some ridiculous special cases
+  if (N < 1) return FALSE;
+
+  if (N == 1) {
+    stats[0].mean    = value[0];
+    stats[0].median  = value[0];
+    stats[0].sigma   = NAN;
+    stats[0].error   = dvalue ? dvalue[0] : NAN;
+    stats[0].chisq   = NAN;
+    stats[0].min     = value[0];
+    stats[0].max     = value[0];
+    stats[0].Upper80 = value[0];
+    stats[0].Lower20 = value[0];
+    stats[0].total   = value[0];
+    return TRUE;
+  }
+
+  // for less than 5 entries, this is equivalent to max,min
+  int N80 = MIN (N-1, 0.8*N);
+  int N20 = MAX (  0, 0.2*N);
+
+  if (weight) {
+    dsortthree (value, dvalue, weight, N);
+  } else {
+    if (dvalue) {
+      dsortpair (value, dvalue, N);
+    } else {
+      dsort (value, N);
+    }
+  }
+
+  // these values do not depend on the errors or weighting scheme
+  if (N % 2) {
+    stats[0].median  = value[(int)(N/2)];
+  } else {
+    stats[0].median  = 0.5*(value[(int)(N/2)] + value[(int)(N/2) - 1]);
+  }
+  stats[0].min     = value[0];
+  stats[0].max     = value[N-1];
+  stats[0].Upper80 = value[N80];
+  stats[0].Lower20 = value[N20];
+
+  switch (stats->statmode) {
+    case VSTATS_MEDIAN:
+      ks = 0;
+      ke = N;
+      Mo = stats[0].median;
+      Nm = N;
+      goto chisq;
+      break;
+    case VSTATS_MEAN:
+    case VSTATS_WT_MEAN:
+      ks = 0;
+      ke = N;
+      break;
+    case VSTATS_INNER_MEAN:
+    case VSTATS_INNER_WTMEAN:
+    case VSTATS_CHI_INNER_MEAN:
+    case VSTATS_CHI_INNER_WTMEAN:
+      ks = 0.25*N + 0.50;
+      ke = 0.75*N + 0.25;
+      if (N <= 3) {
+	ks = 0;
+	ke = N;
+      }
+      break;
+    case VSTATS_NONE:
+      myAbort ("undefined stats");
+  }    
+
+  // 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 == VSTATS_CHI_INNER_MEAN) || (stats->statmode == VSTATS_CHI_INNER_WTMEAN)) {
+    ALLOCATE (chi, double, N);
+    for (i = 0; i < N; i++) {
+      chi[i] = (value[i] - stats[0].median) / dvalue[i];
+    }
+    if (weight) {
+      dsortfour (chi, value, dvalue, weight, N);
+    } else {
+      dsortthree (chi, value, dvalue, N);
+    }
+    free (chi);
+  }
+
+  int WeightedMean = FALSE;
+  WeightedMean |= (stats->statmode == VSTATS_WT_MEAN);
+  WeightedMean |= (stats->statmode == VSTATS_INNER_WTMEAN);
+  WeightedMean |= (stats->statmode == VSTATS_CHI_INNER_WTMEAN);
+
+  /* calculate stats based on the desired weighting scheme */
+  M = dM = Nm = W = R = 0;
+  if (weight) {
+    // the weight value is multiplied by whichever nominal weighting scheme is provided 
+    // thus the user should set weight to 1.0 for nominal weight, and 100 for heavy weight (or so)
+    // and 0.01 for under-weight
+    if (WeightedMean) {
+      for (i = ks; i < ke; i++) {
+	M  += value[i] * weight[i] / SQ(dvalue[i]);
+	W  +=            weight[i] / SQ(dvalue[i]);
+	dM += SQ (weight[i] / dvalue[i]);
+	R  += weight[i] / SQ(dvalue[i]);
+	Nm += 1.0;  
+      }	
+      Mo  = M / W;
+      dMo = sqrt (dM) / R;
+    } else {
+      for (i = ks; i < ke; i++) {
+	M  += value[i] * weight[i];
+	W  +=            weight[i];
+	dM += SQ (weight[i] * dvalue[i]);
+	R  += weight[i];
+	Nm += 1.0;  
+      }	
+      Mo  = M / W;
+      dMo = sqrt (dM) / R;
+    }
+  } else {
+    // NULL weight vector is supplied, revert to standard form (above reverts if weight[i] == 1)
+    if (WeightedMean) {
+      // weighted by inverse-variance
+      for (i = ks; i < ke; i++) {
+	M   += value[i] / SQ (dvalue[i]);
+	dM  += 1.0 / SQ (dvalue[i]);
+	Nm  += 1.0;  
+      }	
+      Mo = M / dM;
+      dMo = sqrt (1.0 / dM);
+    } else {
+      // pure un-weighted
+      for (i = ks; i < ke; i++) {
+	M   += value[i];
+	dM  += SQ (dvalue[i]);
+	Nm  += 1.0;  
+      }	
+      Mo = M / Nm;
+      dMo = sqrt (dM) / Nm;
+    }
+  }
+
+ chisq:
+  /* find sigma and chisq */
+  X2 = dS = 0;
+  for (i = ks; i < ke; i++) {
+    M  = SQ (value[i] - Mo);
+    dM = SQ (dvalue[i]);
+    X2 += M / dM;
+    dS += M;
+  }
+  X2 = X2 / (Nm - 1);
+  dS = sqrt (dS / (Nm - 1));
+
+  stats[0].mean  = Mo;
+  stats[0].Nmeas = Nm;
+  stats[0].chisq = X2;
+  stats[0].sigma = dS;
+  stats[0].error = dMo;
+
+  return (TRUE);
+}
+
+// NOTE: this function resorts the input vectors.  if this is not what you want, you need
+// to either pass a vector copy or save the sequence
+int vstats_getstats_f (float *value, float *dvalue, float *weight, int N, VStatsType *stats) {
+  
+  int i, ks, ke;
+  float Mo, dMo, M, dM, Nm, X2, dS, R, W, *chi;
+
+  myAssert (stats->statmode != VSTATS_NONE, "programming error, liststats mode not set");
+
+  ke = ks = dMo = 0;
+
+  stats[0].Nmeas   = N;
+  stats[0].mean    = NAN;
+  stats[0].median  = NAN;
+  stats[0].sigma   = NAN;
+  stats[0].error   = NAN;
+  stats[0].chisq   = NAN;
+  stats[0].min     = NAN;
+  stats[0].max     = NAN;
+  stats[0].Upper80 = NAN;
+  stats[0].Lower20 = NAN;
+  stats[0].total   = NAN;
+
+  // take care of some ridiculous special cases
+  if (N < 1) return FALSE;
+
+  if (N == 1) {
+    stats[0].mean    = value[0];
+    stats[0].median  = value[0];
+    stats[0].sigma   = NAN;
+    stats[0].error   = dvalue ? dvalue[0] : NAN;
+    stats[0].chisq   = NAN;
+    stats[0].min     = value[0];
+    stats[0].max     = value[0];
+    stats[0].Upper80 = value[0];
+    stats[0].Lower20 = value[0];
+    stats[0].total   = value[0];
+    return TRUE;
+  }
+
+  // for less than 5 entries, this is equivalent to max,min
+  int N80 = MIN (N-1, 0.8*N);
+  int N20 = MAX (  0, 0.2*N);
+
+  if (weight) {
+    fsortthree (value, dvalue, weight, N);
+  } else {
+    if (dvalue) {
+      fsortpair (value, dvalue, N);
+    } else {
+      fsort (value, N);
+    }
+  }
+
+  // these values do not depend on the errors or weighting scheme
+  if (N % 2) {
+    stats[0].median  = value[(int)(N/2)];
+  } else {
+    stats[0].median  = 0.5*(value[(int)(N/2)] + value[(int)(N/2) - 1]);
+  }
+  stats[0].min     = value[0];
+  stats[0].max     = value[N-1];
+  stats[0].Upper80 = value[N80];
+  stats[0].Lower20 = value[N20];
+
+  switch (stats->statmode) {
+    case VSTATS_MEDIAN:
+      ks = 0;
+      ke = N;
+      Mo = stats[0].median;
+      Nm = N;
+      goto chisq;
+      break;
+    case VSTATS_MEAN:
+    case VSTATS_WT_MEAN:
+      ks = 0;
+      ke = N;
+      break;
+    case VSTATS_INNER_MEAN:
+    case VSTATS_INNER_WTMEAN:
+    case VSTATS_CHI_INNER_MEAN:
+    case VSTATS_CHI_INNER_WTMEAN:
+      ks = 0.25*N + 0.50;
+      ke = 0.75*N + 0.25;
+      if (N <= 3) {
+	ks = 0;
+	ke = N;
+      }
+      break;
+    case VSTATS_NONE:
+      myAbort ("undefined stats");
+  }    
+
+  // 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 == VSTATS_CHI_INNER_MEAN) || (stats->statmode == VSTATS_CHI_INNER_WTMEAN)) {
+    ALLOCATE (chi, float, N);
+    for (i = 0; i < N; i++) {
+      chi[i] = (value[i] - stats[0].median) / dvalue[i];
+    }
+    if (weight) {
+      fsortfour (chi, value, dvalue, weight, N);
+    } else {
+      fsortthree (chi, value, dvalue, N);
+    }
+    free (chi);
+  }
+
+  int WeightedMean = FALSE;
+  WeightedMean |= (stats->statmode == VSTATS_WT_MEAN);
+  WeightedMean |= (stats->statmode == VSTATS_INNER_WTMEAN);
+  WeightedMean |= (stats->statmode == VSTATS_CHI_INNER_WTMEAN);
+
+  /* calculate stats based on the desired weighting scheme */
+  M = dM = Nm = W = R = 0;
+  if (weight) {
+    // the weight value is multiplied by whichever nominal weighting scheme is provided 
+    // thus the user should set weight to 1.0 for nominal weight, and 100 for heavy weight (or so)
+    // and 0.01 for under-weight
+    if (WeightedMean) {
+      for (i = ks; i < ke; i++) {
+	M  += value[i] * weight[i] / SQ(dvalue[i]);
+	W  +=            weight[i] / SQ(dvalue[i]);
+	dM += SQ (weight[i] / dvalue[i]);
+	R  += weight[i] / SQ(dvalue[i]);
+	Nm += 1.0;  
+      }	
+      Mo  = M / W;
+      dMo = sqrt (dM) / R;
+    } else {
+      for (i = ks; i < ke; i++) {
+	M  += value[i] * weight[i];
+	W  +=            weight[i];
+	dM += SQ (weight[i] * dvalue[i]);
+	R  += weight[i];
+	Nm += 1.0;  
+      }	
+      Mo  = M / W;
+      dMo = sqrt (dM) / R;
+    }
+  } else {
+    // NULL weight vector is supplied, revert to standard form (above reverts if weight[i] == 1)
+    if (WeightedMean) {
+      // weighted by inverse-variance
+      for (i = ks; i < ke; i++) {
+	M   += value[i] / SQ (dvalue[i]);
+	dM  += 1.0 / SQ (dvalue[i]);
+	Nm  += 1.0;  
+      }	
+      Mo = M / dM;
+      dMo = sqrt (1.0 / dM);
+    } else {
+      // pure un-weighted
+      for (i = ks; i < ke; i++) {
+	M   += value[i];
+	dM  += SQ (dvalue[i]);
+	Nm  += 1.0;  
+      }	
+      Mo = M / Nm;
+      dMo = sqrt (dM) / Nm;
+    }
+  }
+
+ chisq:
+  /* find sigma and chisq */
+  X2 = dS = 0;
+  for (i = ks; i < ke; i++) {
+    M  = SQ (value[i] - Mo);
+    dM = SQ (dvalue[i]);
+    X2 += M / dM;
+    dS += M;
+  }
+  X2 = X2 / (Nm - 1);
+  dS = sqrt (dS / (Nm - 1));
+
+  stats[0].mean  = Mo;
+  stats[0].Nmeas = Nm;
+  stats[0].chisq = X2;
+  stats[0].sigma = dS;
+  stats[0].error = dMo;
+
+  return (TRUE);
+}
+
+// we could define the weight to be the only scale factor:
+// \mu      = \sum (value_i * weight_i) / \sum (weight_i)
+// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) 
+// R = \sum (weight_i^2)
+
+// or, we could define the weight to be a scale factor times the inverse error:
+// \mu      = \sum (value_i * weight_i / sigma_i) / \sum (weight_i)
+// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) 
+// R = \sum (weight_i^2)
+
+//