Index: trunk/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.astro/Makefile	(revision 38062)
@@ -29,4 +29,5 @@
 $(SRC)/czplot.$(ARCH).o	   \
 $(SRC)/cdensify.$(ARCH).o	   \
+$(SRC)/cdhistogram.$(ARCH).o	   \
 $(SRC)/drizzle.$(ARCH).o	   \
 $(SRC)/flux.$(ARCH).o		   \
Index: trunk/Ohana/src/opihi/cmd.astro/cdensify.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.astro/cdensify.c	(revision 38062)
@@ -111,8 +111,4 @@
     switch (PSFTYPE) {
       case IS_DOT:
-	RD_to_XY (&x, &y, rn, *d, &graphmode.coords);
-	Xb = (x - Xmin) / dX;
-	Yb = (y - Ymin) / dY;
-
 	RD_to_XY (&x, &y, rn, *d, &newcoords);
 	Xb = (int) x;
Index: trunk/Ohana/src/opihi/cmd.astro/cdhistogram.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/cdhistogram.c	(revision 38062)
+++ trunk/Ohana/src/opihi/cmd.astro/cdhistogram.c	(revision 38062)
@@ -0,0 +1,135 @@
+# include "data.h"
+
+# define CHECKVAL(ARG) if (!isfinite(ARG)) { gprint (GP_ERR, "illegal value for %s: %f\n", #ARG, ARG); return (FALSE); }
+
+int cdhistogram (int argc, char **argv) {
+
+  int i, Nz, N, Xpix, Ypix;
+  double Xmin, Xmax, dX, Ymin, Ymax, dY;
+  float *val;
+  Buffer *bf;
+  Vector *vr, *vd, *vz, *range;
+  opihi_flt *r, *d, *z, x, y;
+  int kapa;
+  Graphdata graphmode;
+
+  range = NULL;
+  if ((N = get_argument (argc, argv, "-range"))) {
+    remove_argument (N, &argc, argv);
+    if ((range = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  }
+
+  if (!style_args (&graphmode, &argc, argv, &kapa)) return FALSE;
+  double Rmin = graphmode.coords.crval1 - 182.0;
+  double Rmax = graphmode.coords.crval1 + 182.0;
+
+  float dZ = NAN;
+  if ((N = get_argument (argc, argv, "-delta"))) {
+    remove_argument (N, &argc, argv);
+    dZ = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 7) {
+    gprint (GP_ERR, "USAGE: cdhistogram buffer R D value (min) (max) [-delta dval]\n");
+    gprint (GP_ERR, " output buffer is 3D\n");
+    return (FALSE);
+  }
+  
+  if ((bf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((vr = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vd = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vz = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  if (vr[0].Nelements != vd[0].Nelements) return (FALSE);
+  if (vr[0].Nelements != vz[0].Nelements) return (FALSE);
+
+  float Zmin = atof(argv[5]);
+  float Zmax = atof(argv[6]);
+  if (!isfinite(dZ)) {
+    Nz = 100;
+    dZ = (Zmax - Zmin) / (Nz - 1);
+  } else {
+    Nz = (Zmax - Zmin) / dZ + 1;
+  }
+  if (dZ < 0) {
+    gprint (GP_ERR, "invalid value for delta: %f\n", dZ);
+    return (FALSE);
+  }
+  if (Nz > 1000) {
+    gprint (GP_ERR, "warning: delta of %f will result in %d histogram bins\n", dZ, Nz);
+    return (FALSE);
+  }
+
+  if (range) {
+    ResetVector (range, OPIHI_FLT, Nz);
+    for (i = 0; i < range[0].Nelements; i++) {
+      range[0].elements.Flt[i] = Zmin + i*dZ;
+    }
+  }
+
+  REQUIRE_VECTOR_FLT (vr, FALSE); 
+  REQUIRE_VECTOR_FLT (vd, FALSE); 
+  REQUIRE_VECTOR_FLT (vz, FALSE); 
+
+  KapaGetImageRange (kapa, &Xmin, &Xmax, &Ymax, &Ymin, &Xpix, &Ypix);
+  Xmax = graphmode.xmax;
+  Xmin = graphmode.xmin;
+  Ymax = graphmode.ymax;
+  Ymin = graphmode.ymin;
+  dX = (Xmax - Xmin) / (Xpix - 1);
+  dY = (Ymax - Ymin) / (Ypix - 1);
+
+  CHECKVAL(Xmin);
+  CHECKVAL(Xmax);
+  CHECKVAL(dX);
+
+  CHECKVAL(Ymin);
+  CHECKVAL(Ymax);
+  CHECKVAL(dY);
+
+  int Nx = (Xmax - Xmin) / dX + 1;
+  int Ny = (Ymax - Ymin) / dY + 1;
+  
+  Coords newcoords = graphmode.coords;
+  newcoords.cdelt1 *= dX;
+  newcoords.cdelt2 *= dY;
+  newcoords.crpix1 = (newcoords.crpix1 - Xmin) / dX;
+  newcoords.crpix2 = (newcoords.crpix2 - Ymin) / dY;
+
+  gfits_free_matrix (&bf[0].matrix);
+  gfits_free_header (&bf[0].header);
+  CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0);
+  strcpy (bf[0].file, "(empty)");
+  PutCoords (&newcoords, &bf[0].header);
+  
+  // generate the PSF in a local tangent plane
+  Coords coords;
+  InitCoords (&coords, "DEC--TAN");
+
+  r = vr[0].elements.Flt;
+  d = vd[0].elements.Flt;
+  z = vz[0].elements.Flt;
+
+  val = (float *)bf[0].matrix.buffer;
+
+  for (i = 0; i < vr[0].Nelements; i++, r++, d++, z++) {
+    double rn = ohana_normalize_angle (*r);
+    while (rn < Rmin) rn += 360.0;
+    while (rn > Rmax) rn -= 360.0;
+    RD_to_XY (&x, &y, rn, *d, &newcoords);
+    int Xb = x;
+    int Yb = y;
+    int Zb = (*z - Zmin) / dZ;
+
+    if (Xb >= Nx) continue;
+    if (Yb >= Ny) continue;
+    if (Zb >= Nz) continue;
+    if (Xb < 0) continue;
+    if (Yb < 0) continue;
+    if (Zb < 0) continue;
+    val[Xb + Yb*Nx + Zb*Nx*Ny] ++;
+  }
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.astro/init.c	(revision 38062)
@@ -14,4 +14,5 @@
 int czcplot                 PROTO((int, char **));
 int cdensify                PROTO((int, char **));
+int cdhistogram             PROTO((int, char **));
 int drizzle                 PROTO((int, char **));
 int flux                    PROTO((int, char **));
@@ -77,5 +78,6 @@
   {1, "czplot",      czplot,       "plot scaled vectors in sky coordinates"},
   {1, "czcplot",     czcplot,      "plot color-scaled vectors in sky coordinates"},
-  {1, "cdensify",    cdensify,      "vectors to density history on projection"},
+  {1, "cdensify",    cdensify,     "vectors to density history on projection"},
+  {1, "cdhistogram", cdhistogram,  "vectors to 3D histogram on projection"},
   {1, "drizzle",     drizzle,      "transform image to image"},
   {1, "flux",        flux,         "flux in a convex contour"},
Index: trunk/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 38062)
@@ -61,4 +61,5 @@
 $(SRC)/ungridify.$(ARCH).o     \
 $(SRC)/histogram.$(ARCH).o	\
+$(SRC)/tdhistogram.$(ARCH).o	\
 $(SRC)/hermitian1d.$(ARCH).o	\
 $(SRC)/hermitian2d.$(ARCH).o	\
@@ -89,4 +90,5 @@
 $(SRC)/medacc.$(ARCH).o	\
 $(SRC)/mget.$(ARCH).o		\
+$(SRC)/mget3d.$(ARCH).o		\
 $(SRC)/minterpolate.$(ARCH).o	\
 $(SRC)/medimage.$(ARCH).o	\
@@ -130,4 +132,5 @@
 $(SRC)/imspline_apply.$(ARCH).o	\
 $(SRC)/imspline_construct.$(ARCH).o \
+$(SRC)/squash3d.$(ARCH).o	   \
 $(SRC)/imstats.$(ARCH).o	   \
 $(SRC)/style.$(ARCH).o		   \
Index: trunk/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/init.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.data/init.c	(revision 38062)
@@ -51,4 +51,5 @@
 int ungridify        PROTO((int, char **));
 int histogram        PROTO((int, char **));
+int tdhistogram      PROTO((int, char **));
 int hermitian1d      PROTO((int, char **));
 int hermitian2d      PROTO((int, char **));
@@ -80,4 +81,5 @@
 int medacc           PROTO((int, char **));
 int mget             PROTO((int, char **));
+int mget3d           PROTO((int, char **));
 int minterp          PROTO((int, char **));
 int medimage_command PROTO((int, char **));
@@ -117,4 +119,5 @@
 int imspline_apply   PROTO((int, char **));
 int imspline_construct PROTO((int, char **));
+int squash3d         PROTO((int, char **));
 int stats            PROTO((int, char **));
 int imstats          PROTO((int, char **));
@@ -217,4 +220,5 @@
   {1, "header",       header,           "print image header"},
   {1, "histogram",    histogram,        "generate histogram from vector"},
+  {1, "tdhistogram",  tdhistogram,      "generate 2D histogram image from vector set"},
   {1, "hermitian1d",  hermitian1d,      "generate 1-D Hermitian Polynomial"},
   {1, "hermitian2d",  hermitian2d,      "generate 2-D Hermitian Polynomial"},
@@ -243,4 +247,5 @@
   {1, "medacc",       medacc,           "accumulate vector values in another vector"},
   {1, "mget",         mget,             "extract a vector from an image"},
+  {1, "mget3d",       mget3d,           "extract a vector from a 3D image"},
   {1, "imget",        mget,             "extract a vector from an image"},
   {1, "minterp",      minterp,          "interpolate image pixels"},
@@ -289,4 +294,5 @@
   {1, "imspline.apply", imspline_apply, "apply spline fit to generate an image"},
   {1, "imspline.const", imspline_construct, "create spline 2nd deriv. terms"},
+  {1, "squash3d",     squash3d,         "squash 3d buffer to 2d"},
   {1, "stats",        imstats,          "statistics on a portion of an image"},
   {1, "style",        style,            "set the style for graph plots"},
Index: trunk/Ohana/src/opihi/cmd.data/interpolate.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/interpolate.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.data/interpolate.c	(revision 38062)
@@ -1,4 +1,5 @@
 # include "data.h"
 
+// XXX use 'threshold' to interpolate to a value
 int interpolate (int argc, char **argv) {
 
@@ -11,6 +12,7 @@
     gprint (GP_ERR, "USAGE: interpolate Xi Yi Xo Yo\n");
     gprint (GP_ERR, "  Xi Yi - sorted reference vectors\n");
-    gprint (GP_ERR, "  Xo    - output positions\n");
-    gprint (GP_ERR, "  Yo    - output values\n");
+    gprint (GP_ERR, "  Xo    - output positions (vector)\n");
+    gprint (GP_ERR, "  Yo    - output values (vector)\n");
+    gprint (GP_ERR, "  (use 'threshold' to interpolate to a value)\n");
     return (FALSE);
   }
@@ -18,4 +20,6 @@
   if ((xin  = SelectVector (argv[1],  OLDVECTOR, TRUE)) == NULL) return (FALSE);
   if ((yin  = SelectVector (argv[2],  OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  // target positions are a vector
   if ((xout = SelectVector (argv[3],  OLDVECTOR, TRUE)) == NULL) return (FALSE);
   if ((yout = SelectVector (argv[4],  ANYVECTOR, TRUE)) == NULL) return (FALSE);
Index: trunk/Ohana/src/opihi/cmd.data/mcreate.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/mcreate.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.data/mcreate.c	(revision 38062)
@@ -3,20 +3,32 @@
 int mcreate (int argc, char **argv) {
   
-  int Nx, Ny;
+  int N;
   Buffer *buf;
 
+  int Nz = 0;
+  if ((N = get_argument (argc, argv, "-nz"))) {
+    remove_argument (N, &argc, argv);
+    Nz = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   if (argc != 4) {
-    gprint (GP_ERR, "USAGE: mcreate <buffer> Nx Ny\n");
+    gprint (GP_ERR, "USAGE: mcreate <buffer> Nx Ny [-nz Nz]\n");
     return (FALSE);
   }
 
   if ((buf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
-  Nx = atof (argv[2]);
-  Ny = atof (argv[3]);
+  int Nx = atof (argv[2]);
+  int Ny = atof (argv[3]);
 
   /* I should encapsulate this in a create_default_buffer */
   gfits_free_matrix (&buf[0].matrix);
   gfits_free_header (&buf[0].header);
-  CreateBuffer (buf, Nx, Ny, -32, 1.0, 0.0);
+
+  if (Nz) {
+    CreateBuffer3D (buf, Nx, Ny, Nz, -32, 1.0, 0.0);
+  } else {
+    CreateBuffer (buf, Nx, Ny, -32, 1.0, 0.0);
+  }
   return (TRUE);
 }
Index: trunk/Ohana/src/opihi/cmd.data/mget3d.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/mget3d.c	(revision 38062)
+++ trunk/Ohana/src/opihi/cmd.data/mget3d.c	(revision 38062)
@@ -0,0 +1,46 @@
+# include "data.h"
+
+int mget3d (int argc, char **argv) {
+  
+  int i;
+  Buffer *buf;
+  Vector *vec;
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: mget <buffer> <vector> x y\n");
+    return (FALSE);
+  }
+
+  int x = atoi(argv[3]);
+  int y = atoi(argv[4]);
+
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if (buf[0].matrix.Naxes < 3) {
+    gprint (GP_ERR, "buffer is not 3D\n");
+    return FALSE;
+  }
+
+  int Nx = buf[0].matrix.Naxis[0];
+  int Ny = buf[0].matrix.Naxis[1];
+  int Nz = buf[0].matrix.Naxis[2];
+
+  int invalid = FALSE;
+  invalid = invalid || (x < 0);
+  invalid = invalid || (x >= Nx);
+  invalid = invalid || (y < 0);
+  invalid = invalid || (y >= Ny);
+  if (invalid) {
+    gprint (GP_ERR, "selection (%d,%d) out of range\n", x, y);
+    return (FALSE);
+  }
+
+  if ((vec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  ResetVector (vec, OPIHI_FLT, Nz);
+  float *in  = (float *) buf[0].matrix.buffer + x + y*Nx;
+  opihi_flt *out = vec[0].elements.Flt;
+  for (i = 0; i < Nz; i++, in += Nx*Ny, out++) {
+    *out = *in;
+  }
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.data/rd.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/rd.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.data/rd.c	(revision 38062)
@@ -4,11 +4,10 @@
 int rd (int argc, char **argv) {
   
-  int i, N, status, plane, Nplane, extend, Nextend, Nskip, JustHead, blank;
-  int ccdsel, done, Nword, IsCompressed;
-  char region[512], *ccdid, *filename;
-  FILE *f;
+  int i, N, Nskip, blank;
+  int done, Nword;
+  char region[512];
   Buffer *buf;
 
-  JustHead = FALSE;
+  int JustHead = FALSE;
   if ((N = get_argument (argc, argv, "-head"))) {
     remove_argument (N, &argc, argv);
@@ -16,5 +15,5 @@
   }
 
-  plane = 1;
+  int plane = -1;
   if ((N = get_argument (argc, argv, "-plane"))) {
     remove_argument (N, &argc, argv);
@@ -23,6 +22,6 @@
   }
 
-  extend = FALSE;
-  Nextend = -1;
+  int extend = FALSE;
+  int Nextend = -1;
   if ((N = get_argument (argc, argv, "-x"))) {
     remove_argument (N, &argc, argv);
@@ -32,6 +31,6 @@
   }
 
-  ccdsel = FALSE;
-  ccdid = (char *) NULL;
+  int ccdsel = FALSE;
+  char *ccdid = NULL;
   if ((N = get_argument (argc, argv, "-n"))) {
     remove_argument (N, &argc, argv);
@@ -52,5 +51,5 @@
 
   /* test if file exists */
-  f = fopen (argv[2], "r");
+  FILE *f = fopen (argv[2], "r");
   if (f == (FILE *) NULL) {
     gprint (GP_ERR, "file %s not found\n", argv[2]);
@@ -67,10 +66,10 @@
 
   /* save file name */
-  filename = filebasename (argv[2]);
+  char *filename = filebasename (argv[2]);
   strcpy (buf[0].file, filename);
   free (filename);
 
-  status = FALSE;
-  IsCompressed = FALSE;
+  int status = FALSE;
+  int IsCompressed = FALSE;
 
   /*** advance to the correct FITS extension ***/
@@ -168,15 +167,24 @@
 
   /* check for valid plane */
-  Nplane = buf[0].header.Naxis[2];
-  if (Nplane == 0) Nplane = 1;
-  if (plane > Nplane) {
-    gprint (GP_ERR, "-plane is too large: %d total planes\n", Nplane);
-    DeleteBuffer (buf);
-    fclose (f);
-    return (FALSE);
+  int Nz = buf[0].header.Naxis[2];
+  if (plane >= 0) {
+    // we are requesting a specific plane (-1 : all data)
+    int tooFar = Nz ? (plane >= Nz) : (plane > Nz);
+    if (tooFar) {
+      gprint (GP_ERR, "-plane is too large: %d total planes\n", Nz);
+      DeleteBuffer (buf);
+      fclose (f);
+      return (FALSE);
+    }
   }
 
   /* load matrix data */
   if (IsCompressed) {
+    if (plane > -1) {
+      gprint (GP_ERR, "-plane incompatible with compressed image\n");
+      DeleteBuffer (buf);
+      fclose (f);
+      return (FALSE);
+    }
     FTable ftable;
     Header theader;
@@ -191,6 +199,15 @@
     // XXX this currently does not work for a cube (we get a cube back, not a specific plane)
   } else {
-    sprintf (region, "-1 -1 -1 -1 %d %d", (plane - 1), plane);
-    status = gfits_fread_matrix_segment (f, &buf[0].matrix, &buf[0].header, region);
+    if (plane > -1) {
+      // read a single plane into a 2D image
+      sprintf (region, "-1 -1 -1 -1 %d %d", (plane - 1), plane);
+      status = gfits_fread_matrix_segment (f, &buf[0].matrix, &buf[0].header, region);
+      buf[0].header.Naxis[2] = 0;
+      buf[0].header.Naxes = 2;
+      gfits_modify (&buf[0].header, "NAXIS", "%d", 1, 2);
+      gfits_modify (&buf[0].header, "NAXIS3", "%d", 1, 0);
+    } else {
+      status = gfits_fread_matrix (f, &buf[0].matrix, &buf[0].header);
+    }
   }
   fclose (f);
@@ -200,12 +217,4 @@
     DeleteBuffer (buf);
     return (FALSE);
-  }
-
-  /* adjust buffer to represent 2D data */
-  if (Nplane > 1) {
-    buf[0].header.Naxis[2] = 0;
-    buf[0].header.Naxes = 2;
-    gfits_modify (&buf[0].header, "NAXIS", "%d", 1, 2);
-    gfits_delete (&buf[0].header, "NAXIS3", 1);
   }
 
@@ -215,4 +224,6 @@
     buf[0].header.Naxis[1] = 1;
     buf[0].matrix.Naxis[1] = 1;
+    gfits_modify (&buf[0].header, "NAXIS", "%d", 1, 2);
+    gfits_modify (&buf[0].header, "NAXIS2", "%d", 1, 1);
   }    
 
Index: trunk/Ohana/src/opihi/cmd.data/read_vectors.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/read_vectors.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.data/read_vectors.c	(revision 38062)
@@ -372,4 +372,10 @@
   }
 
+  int VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = TRUE;
+  }
+
   int start = 0;
   int Nrows = -1; // -1 : read entire table
@@ -399,5 +405,6 @@
   // }
 
-  if (argc < 2) ESCAPE ("USAGE: read -fits extension [-extnum] [-keyword key] name name ...\n");
+  if ((argc < 2) && !getSizes) ESCAPE ("USAGE: read -fits extension [-extnum] [-keyword key] name name ...\n");
+  if ((argc != 1) && getSizes) ESCAPE ("USAGE: read -sizes -fits extension [-extnum] (does not read data values)\n");
 
   if (f == NULL) ESCAPE ("file not found\n");
@@ -436,5 +443,9 @@
     if (start >= header.Naxis[1]) ESCAPE ("invalid range: start >= Ny (%d)\n", header.Naxis[1]);
     if (Nrows < 0) ESCAPE ("invalid range: Nrows < 0\n");
-    if (start + Nrows > header.Naxis[1]) ESCAPE ("invalid range: start + Nrows > Ny (%d)\n", header.Naxis[1]);
+
+    // just a warning:
+    if (start + Nrows > header.Naxis[1]) {
+      if (VERBOSE) gprint (GP_ERR, "NOTE: reading last block will return only %d rows\n", header.Naxis[1] - start);
+    }
 
     // Ny = 100, start = 0, Nrows = -1 -> Nrows => 100
@@ -472,4 +483,11 @@
       }
 
+      if (getSizes) {
+	read_table_sizes (&header);
+	if (CCDKeyword != NULL) free (CCDKeyword); 
+	gfits_free_header (&header); 
+	return TRUE;
+      }
+
       if (Nrows == -1) {
 	Nrows = header.Naxis[1] - start;
@@ -478,5 +496,9 @@
       if (start >= header.Naxis[1]) ESCAPE ("invalid range: start >= Ny (%d)\n", header.Naxis[1]);
       if (Nrows < 0) ESCAPE ("invalid range: Nrows < 0\n");
-      if (start + Nrows > header.Naxis[1]) ESCAPE ("invalid range: start + Nrows > Ny (%d)\n", header.Naxis[1]);
+
+      // just a warning:
+      if (start + Nrows > header.Naxis[1]) {
+	if (VERBOSE) gprint (GP_ERR, "NOTE: reading last block will return only %d rows\n", header.Naxis[1] - start);
+      }
 
       if (!gfits_fread_ftable_range (f, padIfShort, &table, start, Nrows)) ESCAPE ("error reading table for extension %d\n", Nextend);
@@ -603,5 +625,5 @@
 
   set_int_variable ("table:Nx", header->Naxis[0]);
-  set_int_variable ("table:Nx", header->Naxis[0]);
+  set_int_variable ("table:Ny", header->Naxis[1]);
   set_int_variable ("table:Nfields", Nfields);
 
Index: trunk/Ohana/src/opihi/cmd.data/squash3d.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/squash3d.c	(revision 38062)
+++ trunk/Ohana/src/opihi/cmd.data/squash3d.c	(revision 38062)
@@ -0,0 +1,113 @@
+# include "data.h"
+
+enum {SQUASH_NONE, SQUASH_X, SQUASH_Y, SQUASH_Z};
+
+int squash3d (int argc, char **argv) {
+  
+  int ox, oy, oz;
+  Buffer *src;
+  Buffer *tgt;
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: squash3d <3d> <2d> (dir)\n");
+    gprint (GP_ERR, " dir: x -> squash in x-dir, y -> squash in y-dir, z -> squash in z-dir\n");
+    return (FALSE);
+  }
+
+  int dir = SQUASH_NONE;
+  if (!strcasecmp (argv[3], "x")) dir = SQUASH_X;
+  if (!strcasecmp (argv[3], "y")) dir = SQUASH_Y;
+  if (!strcasecmp (argv[3], "z")) dir = SQUASH_Z;
+  if (!dir) {
+    gprint (GP_ERR, "invalid direction %s\n", argv[3]);
+    return FALSE;
+  }
+
+  if ((src = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if (src[0].matrix.Naxes < 3) {
+    gprint (GP_ERR, "buffer is not 3D\n");
+    return FALSE;
+  }
+
+  float *iBuf  = (float *) src[0].matrix.buffer;
+
+  if ((tgt = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  gfits_free_matrix (&tgt[0].matrix);
+  gfits_free_header (&tgt[0].header);
+
+  switch (dir) {
+    case SQUASH_X:
+      {     
+	int iNx = src[0].matrix.Naxis[0];
+	int iNy = src[0].matrix.Naxis[1];
+	int iNz = src[0].matrix.Naxis[2];
+
+	// output is Nz,Ny
+	int oNx = iNz;
+	int oNy = iNy;
+	CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0);
+	float *oBuf  = (float *) tgt[0].matrix.buffer;
+
+	for (ox = 0; ox < oNx; ox ++) {
+	  for (oy = 0; oy < oNy; oy ++) {
+	    float val = 0.0;
+	    for (oz = 0; oz < iNx; oz ++) { // src x is tgt z
+	      val += iBuf[oz + oy*iNx + ox*iNx*iNy];
+	    }
+	    oBuf[ox + oy*oNx] = val;
+	  }
+	}
+      }
+      break;
+      
+    case SQUASH_Y:
+      {     
+	int iNx = src[0].matrix.Naxis[0];
+	int iNy = src[0].matrix.Naxis[1];
+	int iNz = src[0].matrix.Naxis[2];
+
+	// output is Nx,Nz
+	int oNx = iNx;
+	int oNy = iNz;
+	CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0);
+	float *oBuf  = (float *) tgt[0].matrix.buffer;
+
+	for (ox = 0; ox < oNx; ox ++) {
+	  for (oy = 0; oy < oNy; oy ++) { // src z is tgt y
+	    float val = 0.0;
+	    for (oz = 0; oz < iNy; oz ++) {
+	      val += iBuf[ox + oz*iNx + oy*iNx*iNy];
+	    }
+	    oBuf[ox + oy*oNx] = val;
+	  }
+	}
+      }
+      break;
+      
+    case SQUASH_Z:
+      {     
+	int iNx = src[0].matrix.Naxis[0];
+	int iNy = src[0].matrix.Naxis[1];
+	int iNz = src[0].matrix.Naxis[2];
+
+	// output is Nx,Ny
+	int oNx = iNx;
+	int oNy = iNy;
+	CreateBuffer (tgt, oNx, oNy, -32, 0.0, 1.0);
+	float *oBuf  = (float *) tgt[0].matrix.buffer;
+
+	for (ox = 0; ox < oNx; ox ++) {
+	  for (oy = 0; oy < oNy; oy ++) {
+	    float val = 0.0;
+	    for (oz = 0; oz < iNz; oz ++) {
+	      val += iBuf[ox + oy*iNx + oz*iNx*iNy];
+	    }
+	    oBuf[ox + oy*oNx] = val;
+	  }
+	}
+      }
+      break;
+  }
+
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.data/tdhistogram.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/tdhistogram.c	(revision 38062)
+++ trunk/Ohana/src/opihi/cmd.data/tdhistogram.c	(revision 38062)
@@ -0,0 +1,176 @@
+# include "data.h"
+
+# define CHECKVAL(ARG) if (!isfinite(ARG)) { gprint (GP_ERR, "illegal value for %s: %f\n", #ARG, ARG); return (FALSE); }
+
+int tdhistogram (int argc, char **argv) {
+
+  int i, Nx, Ny, Nz, N;
+  float *val;
+  Buffer *bf;
+  Vector *vx, *vy, *vz, *range;
+  opihi_flt *x, *y, *z;
+
+  int reuse = FALSE;
+  if ((N = get_argument (argc, argv, "-reuse"))) {
+    remove_argument (N, &argc, argv);
+    reuse = TRUE;
+  }
+
+  range = NULL;
+  if ((N = get_argument (argc, argv, "-range"))) {
+    remove_argument (N, &argc, argv);
+    if ((range = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  }
+
+  double dx = NAN;
+  if ((N = get_argument (argc, argv, "-dx"))) {
+    remove_argument (N, &argc, argv);
+    dx = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  double dy = NAN;
+  if ((N = get_argument (argc, argv, "-dy"))) {
+    remove_argument (N, &argc, argv);
+    dy = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  double dz = NAN;
+  if ((N = get_argument (argc, argv, "-dz"))) {
+    remove_argument (N, &argc, argv);
+    dz = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  int valid = (!reuse && (argc == 11)) || (reuse && (argc == 5));
+  if (!valid) {
+    gprint (GP_ERR, "USAGE: tdhistogram buffer x y z (Xmin) (Xmax) (Ymin) (Ymax) (Zmin) (Zmax) [-dx dx] [-dy dy] [-dz dz]\n");
+    gprint (GP_ERR, "   OR: tdhistogram buffer x y z -reuse\n");
+    gprint (GP_ERR, " output buffer is 3D\n");
+    return (FALSE);
+  }
+  
+  if ((bf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((vx = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vy = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vz = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+
+  if (vx[0].Nelements != vy[0].Nelements) return (FALSE);
+  if (vx[0].Nelements != vz[0].Nelements) return (FALSE);
+
+  double Xmin, Xmax, Ymin, Ymax, Zmin, Zmax;
+
+  if (!reuse) {
+    Xmin = atof(argv[5]);
+    Xmax = atof(argv[6]);
+    if (!isfinite(dx)) {
+      Nx = 100;
+      dx = (Xmax - Xmin) / (Nx - 1);
+    } else {
+      Nx = (Xmax - Xmin) / dx + 1;
+    }
+    if (dx < 0) {
+      gprint (GP_ERR, "invalid value for delta: %f\n", dx);
+      return (FALSE);
+    }
+    Ymin = atof(argv[7]);
+    Ymax = atof(argv[8]);
+    if (!isfinite(dy)) {
+      Ny = 100;
+      dy = (Ymax - Ymin) / (Ny - 1);
+    } else {
+      Ny = (Ymax - Ymin) / dy + 1;
+    }
+    if (dy < 0) {
+      gprint (GP_ERR, "invalid value for delta: %f\n", dy);
+      return (FALSE);
+    }
+    Zmin = atof(argv[9]);
+    Zmax = atof(argv[10]);
+    if (!isfinite(dz)) {
+      Nz = 100;
+      dz = (Zmax - Zmin) / (Nz - 1);
+    } else {
+      Nz = (Zmax - Zmin) / dz + 1;
+    }
+    if (dz < 0) {
+      gprint (GP_ERR, "invalid value for delta: %f\n", dz);
+      return (FALSE);
+    }
+
+    if (Nz > 1000) {
+      gprint (GP_ERR, "warning: delta of %f will result in %d histogram bins\n", dz, Nz);
+      return (FALSE);
+    }
+  } else {
+    gfits_scan (&bf[0].header, "XMIN", "%lf", 1, &Xmin);
+    gfits_scan (&bf[0].header, "XMAX", "%lf", 1, &Xmax);
+    gfits_scan (&bf[0].header, "XDEL", "%lf", 1, &dx);
+    gfits_scan (&bf[0].header, "YMIN", "%lf", 1, &Ymin);
+    gfits_scan (&bf[0].header, "YMAX", "%lf", 1, &Ymax);
+    gfits_scan (&bf[0].header, "YDEL", "%lf", 1, &dy);
+    gfits_scan (&bf[0].header, "ZMIN", "%lf", 1, &Zmin);
+    gfits_scan (&bf[0].header, "ZMAX", "%lf", 1, &Zmax);
+    gfits_scan (&bf[0].header, "ZDEL", "%lf", 1, &dz);
+    Nx = bf[0].header.Naxis[0];
+    Ny = bf[0].header.Naxis[1];
+    Nz = bf[0].header.Naxis[2];
+  }
+
+  REQUIRE_VECTOR_FLT (vx, FALSE); 
+  REQUIRE_VECTOR_FLT (vy, FALSE); 
+  REQUIRE_VECTOR_FLT (vz, FALSE); 
+
+  CHECKVAL(Xmin);
+  CHECKVAL(Xmax);
+  CHECKVAL(dx);
+
+  CHECKVAL(Ymin);
+  CHECKVAL(Ymax);
+  CHECKVAL(dy);
+
+  if (range) {
+    ResetVector (range, OPIHI_FLT, Nz);
+    for (i = 0; i < range[0].Nelements; i++) {
+      range[0].elements.Flt[i] = Zmin + i*dz;
+    }
+  }
+
+  if (!reuse) {
+    gfits_free_matrix (&bf[0].matrix);
+    gfits_free_header (&bf[0].header);
+    CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0);
+    strcpy (bf[0].file, "(empty)");
+
+    gfits_modify (&bf[0].header, "XMIN", "%lf", 1, Xmin);
+    gfits_modify (&bf[0].header, "XMAX", "%lf", 1, Xmax);
+    gfits_modify (&bf[0].header, "XDEL", "%lf", 1, dx);
+    gfits_modify (&bf[0].header, "YMIN", "%lf", 1, Ymin);
+    gfits_modify (&bf[0].header, "YMAX", "%lf", 1, Ymax);
+    gfits_modify (&bf[0].header, "YDEL", "%lf", 1, dy);
+    gfits_modify (&bf[0].header, "ZMIN", "%lf", 1, Zmin);
+    gfits_modify (&bf[0].header, "ZMAX", "%lf", 1, Zmax);
+    gfits_modify (&bf[0].header, "ZDEL", "%lf", 1, dz);
+  }
+  
+  x = vx[0].elements.Flt;
+  y = vy[0].elements.Flt;
+  z = vz[0].elements.Flt;
+
+  val = (float *) bf[0].matrix.buffer;
+
+  for (i = 0; i < vx[0].Nelements; i++, x++, y++, z++) {
+    int Xb = (*x - Xmin) / dx;
+    int Yb = (*y - Ymin) / dy;
+    int Zb = (*z - Zmin) / dz;
+
+    if (Xb >= Nx) continue;
+    if (Yb >= Ny) continue;
+    if (Zb >= Nz) continue;
+    if (Xb < 0) continue;
+    if (Yb < 0) continue;
+    if (Zb < 0) continue;
+    val[Xb + Yb*Nx + Zb*Nx*Ny] ++;
+  }
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.data/tv.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/tv.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.data/tv.c	(revision 38062)
@@ -28,4 +28,15 @@
   }
 
+  int plane = 0;
+  if ((N = get_argument (argc, argv, "-plane"))) {
+    remove_argument (N, &argc, argv);
+    plane = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if (plane < 0) {
+    gprint (GP_ERR, " ERROR: -plane (plane) : cannot be negative\n");
+    return (FALSE);
+  }
+
   /* shell exits on pipe close, FIX */
   if ((N = get_argument (argc, argv, "-kill"))) {
@@ -55,7 +66,16 @@
   GetCoords (&coords, &buf[0].header);
   
-  image.data1d = (float *) buf[0].matrix.buffer;
   image.Nx = buf[0].matrix.Naxis[0];
   image.Ny = buf[0].matrix.Naxis[1];
+
+  int tooBig = buf[0].matrix.Naxis[2] ? (plane >= buf[0].matrix.Naxis[2]) : plane > 0;
+  if (tooBig) {
+    gprint (GP_ERR, " ERROR: -plane (plane) : out of bounds (%d vs %d)\n", plane, buf[0].matrix.Naxis[2]);
+    return (FALSE);
+  }
+  int Npix2D = image.Nx * image.Ny;
+
+  float *imdata = (float *) buf[0].matrix.buffer;
+  image.data1d = &imdata[plane*Npix2D];
 
   // send only the root of the file, not the full path
Index: trunk/Ohana/src/opihi/cmd.data/wd.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/wd.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.data/wd.c	(revision 38062)
@@ -81,14 +81,9 @@
   memcpy (temp_header.buffer, buf[0].header.buffer, temp_header.datasize);
 
-  if (temp_header.Naxes) {
-    // the inBlank value probably does not matter: temp_matrix is float, so nan is used
-    gfits_convert_format (&temp_header, &temp_matrix, outBitpix, outScale, outZero, 0xffff, outUnsign);
-  } else {
-    gfits_modify (&temp_header, "BITPIX", "%d", 1, outBitpix);
-    gfits_modify (&temp_header, "BSCALE", "%lf", 1, outScale);
-    gfits_modify (&temp_header, "BZERO",  "%lf", 1, outZero);
-    gfits_modify_alt (&temp_header, "UNSIGN", "%t", 1, outUnsign);
-  }
+  gfits_convert_format (&temp_header, &temp_matrix, outBitpix, outScale, outZero, 0xffff, outUnsign);
 
+  // Extend puts the output matrix in the first available non-PHU slot (ie, the last one)
+  // 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;
Index: trunk/Ohana/src/opihi/cmd.data/zplot.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/zplot.c	(revision 38002)
+++ trunk/Ohana/src/opihi/cmd.data/zplot.c	(revision 38062)
@@ -4,11 +4,34 @@
   
   char *outname = NULL;
-  int i, kapa, valid, size;
+  int i, N, kapa, valid, size;
   opihi_flt *out;
   double min, range;
   Graphdata graphmode;
-  Vector *xvec, *yvec, *zvec, Zvec;
+  Vector *xvec, *yvec, *zvec, *dxmvec, *dxpvec, *dymvec, *dypvec, Zvec;
 
   if (!style_args (&graphmode, &argc, argv, &kapa)) return (FALSE);
+
+  /* decide on error bars */
+  dxmvec = dxpvec = dymvec = dypvec = NULL;
+  if ((N = get_argument (argc, argv, "-dx"))) {
+    remove_argument (N, &argc, argv);
+    if ((dxmvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "+dx"))) {
+    remove_argument (N, &argc, argv);
+    if ((dxpvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-dy"))) {
+    remove_argument (N, &argc, argv);
+    if ((dymvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "+dy"))) {
+    remove_argument (N, &argc, argv);
+    if ((dypvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  }
 
   valid  = (argc == 6);
@@ -40,4 +63,13 @@
   }
 
+  /* set errorbar mode (these are NOT sticky) */
+  graphmode.etype = 0;
+  if ((dymvec != NULL) && (dypvec == NULL)) dypvec = dymvec;
+  if ((dypvec != NULL) && (dymvec == NULL)) dymvec = dypvec;
+  if ((dypvec != NULL) || (dymvec != NULL)) graphmode.etype |= 0x01;
+  if ((dxmvec != NULL) && (dxpvec == NULL)) dxpvec = dxmvec;
+  if ((dxpvec != NULL) && (dxmvec == NULL)) dxmvec = dxpvec;
+  if ((dxpvec != NULL) || (dxmvec != NULL)) graphmode.etype |= 0x02;
+  
   /* find vectors */
   if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
@@ -56,7 +88,13 @@
     return (FALSE);
   }
+  if (dypvec && (dypvec->Nelements != xvec->Nelements)) goto mismatch;
+  if (dymvec && (dymvec->Nelements != xvec->Nelements)) goto mismatch;
+  if (dxpvec && (dxpvec->Nelements != xvec->Nelements)) goto mismatch;
+  if (dxmvec && (dxmvec->Nelements != xvec->Nelements)) goto mismatch;
+
   SetVector (&Zvec, OPIHI_FLT, zvec[0].Nelements);
   out = Zvec.elements.Flt;
  
+  // note actual size is 3.3x plot -sz sizes. (DrawObjects.c:399)
   if (zvec[0].type == OPIHI_FLT) {
     opihi_flt *in = zvec[0].elements.Flt;
@@ -81,13 +119,31 @@
   graphmode.style = 2; /* plot points */
   graphmode.size = -1; /* point size determined by Zvec */
-  graphmode.etype = 0; /* no errorbars */
   PlotVectorTriplet (kapa, xvec, yvec, &Zvec, mask, &graphmode);
+  if (graphmode.etype & 0x01) {
+    PlotVectorSingle (kapa, dymvec, mask, "dym");
+    PlotVectorSingle (kapa, dypvec, mask, "dyp");
+  }
+  if (graphmode.etype & 0x02) {
+    PlotVectorSingle (kapa, dxmvec, mask, "dxm");
+    PlotVectorSingle (kapa, dxpvec, mask, "dxp");
+  }
 
   free (Zvec.elements.Ptr);
   if (mask) free (mask);
-  DeleteNamedVector (outname);
+
+  if (outname) {
+    DeleteNamedVector (outname);
+    free (outname);
+  }
 
   return (TRUE);
 
+mismatch:
+  gprint (GP_ERR, "error and data vector lengths are mismatched\n");
+  if (outname) {
+    DeleteNamedVector (outname);
+    free (outname);
+  }
+  return (FALSE);
 }
 
Index: trunk/Ohana/src/opihi/doc/BufferMath.txt
===================================================================
--- trunk/Ohana/src/opihi/doc/BufferMath.txt	(revision 38062)
+++ trunk/Ohana/src/opihi/doc/BufferMath.txt	(revision 38062)
@@ -0,0 +1,46 @@
+
+Opihi has the main data types: scalar, vector, buffer.  these are all
+handled somewhat different, with different internal structures.  
+
+In most places, the buffer is essentially an image or a 2D matrix.
+However, it is not really restricted as such.  I recently (2015.01.16)
+updated dvomath operations to handle 3D images as well as 2D images.  
+
+Here are some of the places where the difference matters:
+
+* stack_math.c & vector OP matrix and matrix OP vector: these
+  operations are explicitly only valid for 2D matrices (and order
+  determines if the vector is applied to a row (vector first) or a
+  column (vector last).
+
+* stack_math.c & xramp(a) or yramp(a) : these only are applied to the
+  1st 2D plane of an N-D matrix
+
+* check_stack only tests Nx,Ny,Nz.  it probably ignores higher
+  dimensions, but probably does not disallow them.
+
+* concepts like A[x][y] only work for 2D buffers
+
+* lots of the functions which operate on buffers assume the buffer is
+  an image.  Most of the functions in cmd.astro would not be
+  well-defined for a higher-dimension object.  these should probably
+  at least test for 2D and exit / warn if the selected buffer is not
+  2D
+
+* for some of the functions in cmd.data, a higher-dimension object
+  could make sense, but is currently not allowed.  for example, tv A
+  could be extended to display a given slice.  
+
+* rd and wd should really be extended to allow read/write to/from a 3D
+  file.  actually, rd already can read a plane and wd used to be able
+  to write to a plane.  
+
+* some of the subset-like operations should be extended:
+
+  * mget, mset -- i've added mget3d, but a better abstraction would be
+    good
+   
+  * extract is a problem anyway
+
+  * dimendown and dimenup have similar issues
+
Index: trunk/Ohana/src/opihi/dvo/avextract.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/avextract.c	(revision 38002)
+++ trunk/Ohana/src/opihi/dvo/avextract.c	(revision 38062)
@@ -221,4 +221,5 @@
 
       for (n = 0; n < Nfields; n++) {
+	// we are passing in the *first* measure, but average->Nmeasure gives the count
 	values[n] = dbExtractAverages (average, secfilt, measure, lensobj, starpar, &fields[n]);
       }
Index: trunk/Ohana/src/opihi/include/dvomath.h
===================================================================
--- trunk/Ohana/src/opihi/include/dvomath.h	(revision 38002)
+++ trunk/Ohana/src/opihi/include/dvomath.h	(revision 38062)
@@ -187,4 +187,5 @@
 int           ListBuffersToList     PROTO((char *name));
 int           CreateBuffer          PROTO((Buffer *buf, int Nx, int Ny, int bitpix, float bzero, float bscale));
+int           CreateBuffer3D        PROTO((Buffer *buf, int Nx, int Ny, int Nz, int bitpix, float bzero, float bscale));
 int           ResetBuffer           PROTO((Buffer *buf, int Nx, int Ny, int bitpix, float bzero, float bscale));
 Buffer       *SelectBuffer          PROTO((char *name, int mode, int verbose));
Index: trunk/Ohana/src/opihi/lib.shell/BufferOps.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/BufferOps.c	(revision 38002)
+++ trunk/Ohana/src/opihi/lib.shell/BufferOps.c	(revision 38062)
@@ -120,4 +120,31 @@
 }
   
+// buffer is 1D array referenced in the order:
+// buffer[x + Nx*y + Nx*Ny*z + ...]
+int CreateBuffer3D (Buffer *buf, int Nx, int Ny, int Nz, int bitpix, float bzero, float bscale) {
+
+  /* store the default output values */
+  gfits_init_header (&buf[0].header);
+
+  /* assign the necessary internal values */
+  buf[0].header.bitpix   = -32;
+  buf[0].header.Naxes = 3;
+  buf[0].header.Naxis[0] = Nx;
+  buf[0].header.Naxis[1] = Ny;
+  buf[0].header.Naxis[2] = Nz;
+
+  buf[0].bitpix = bitpix;
+  buf[0].bzero  = bzero;
+  buf[0].bscale = bscale;
+  
+  /* make some test of the validity of the values */
+
+  /* create the appropriate header and matrix */
+  gfits_create_header (&buf[0].header);
+  gfits_create_matrix (&buf[0].header, &buf[0].matrix);
+
+  return (TRUE);
+}
+  
 /* copy data from in to out - new memory space */
 int CopyNamedBuffer (char *out, char *in) {
Index: trunk/Ohana/src/opihi/lib.shell/check_stack.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/check_stack.c	(revision 38002)
+++ trunk/Ohana/src/opihi/lib.shell/check_stack.c	(revision 38062)
@@ -6,8 +6,8 @@
 int check_stack (StackVar *stack, int Nstack, int validsize) {
 
-  int i, Nx, Ny, Nv, size;
+  int i, Nx, Ny, Nz, Nv, size;
   char *c1, *c2;
 
-  Nv = Nx = Ny = -1;
+  Nv = Nx = Ny = Nz = -1;
 
   for (i = 0; i < Nstack; i++) {
@@ -46,6 +46,9 @@
 	  Nx = stack[i].buffer[0].matrix.Naxis[0];
 	  Ny = stack[i].buffer[0].matrix.Naxis[1];
+	  Nz = stack[i].buffer[0].matrix.Naxis[2];
 	} 
-	if ((Nx != stack[i].buffer[0].matrix.Naxis[0]) && (Ny != stack[i].buffer[0].matrix.Naxis[1])) {
+	if ((Nx != stack[i].buffer[0].matrix.Naxis[0]) && 
+	    (Ny != stack[i].buffer[0].matrix.Naxis[1]) && 
+	    (Nz != stack[i].buffer[0].matrix.Naxis[2])) {
 	  push_error ("dimensions don't match");
 	  return (-1);
Index: trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c	(revision 38002)
+++ trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c	(revision 38062)
@@ -60,4 +60,5 @@
     if (!strcmp (argv[i], "xramp"))  { type = ST_UNARY; goto gotit; }
     if (!strcmp (argv[i], "yramp"))  { type = ST_UNARY; goto gotit; }
+    if (!strcmp (argv[i], "zramp"))  { type = ST_UNARY; goto gotit; }
     if (!strcmp (argv[i], "ramp"))   { type = ST_UNARY; goto gotit; }
     if (!strcmp (argv[i], "zero"))   { type = ST_UNARY; goto gotit; }
Index: trunk/Ohana/src/opihi/lib.shell/evaluate_stack.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/evaluate_stack.c	(revision 38002)
+++ trunk/Ohana/src/opihi/lib.shell/evaluate_stack.c	(revision 38062)
@@ -121,5 +121,5 @@
       got_three_op:
 	if (!status) {
-	  snprintf (line, 512, "syntax error: invalid operand for binary operation: %s or %s or %s\n", stack[i-1].name, stack[i-2].name, stack[i-3].name);
+	  snprintf (line, 512, "syntax error: invalid operand for trinary operation: %s or %s or %s\n", stack[i-1].name, stack[i-2].name, stack[i-3].name);
 	  push_error (line);
 	  clear_stack (&tmp_stack);
Index: trunk/Ohana/src/opihi/lib.shell/stack_math.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/stack_math.c	(revision 38002)
+++ trunk/Ohana/src/opihi/lib.shell/stack_math.c	(revision 38062)
@@ -80,4 +80,70 @@
       break;								\
     }									\
+    if ((V1->vector->type == OPIHI_FLT) && (V2->vector->type == OPIHI_FLT) && (V3->vector->type == OPIHI_INT)) { \
+      CopyVector (OUT[0].vector, V1[0].vector);				\
+      opihi_flt *M1  =  V1[0].vector[0].elements.Flt;			\
+      opihi_flt *M2  =  V2[0].vector[0].elements.Flt;			\
+      opihi_int *M3  =  V3[0].vector[0].elements.Int;			\
+      opihi_flt *out = OUT[0].vector[0].elements.Flt;			\
+      for (i = 0; i < Nx; i++, out++, M1++, M2++, M3++) {		\
+  	*out = OP;							\
+      }									\
+      break;								\
+    }									\
+    if ((V1->vector->type == OPIHI_FLT) && (V2->vector->type == OPIHI_INT) && (V3->vector->type == OPIHI_FLT)) { \
+      CopyVector (OUT[0].vector, V1[0].vector);				\
+      opihi_flt *M1  =  V1[0].vector[0].elements.Flt;			\
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;			\
+      opihi_flt *M3  =  V3[0].vector[0].elements.Flt;			\
+      opihi_flt *out = OUT[0].vector[0].elements.Flt;			\
+      for (i = 0; i < Nx; i++, out++, M1++, M2++, M3++) {		\
+  	*out = OP;							\
+      }									\
+      break;								\
+    }									\
+    if ((V1->vector->type == OPIHI_FLT) && (V2->vector->type == OPIHI_INT) && (V3->vector->type == OPIHI_INT)) { \
+      CopyVector (OUT[0].vector, V2[0].vector);				\
+      opihi_flt *M1  =  V1[0].vector[0].elements.Flt;			\
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;			\
+      opihi_int *M3  =  V3[0].vector[0].elements.Int;			\
+      opihi_int *out = OUT[0].vector[0].elements.Int;			\
+      for (i = 0; i < Nx; i++, out++, M1++, M2++, M3++) {		\
+  	*out = OP;							\
+      }									\
+      break;								\
+    }									\
+    if ((V1->vector->type == OPIHI_INT) && (V2->vector->type == OPIHI_FLT) && (V3->vector->type == OPIHI_FLT)) { \
+      CopyVector (OUT[0].vector, V2[0].vector);				\
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;			\
+      opihi_flt *M2  =  V2[0].vector[0].elements.Flt;			\
+      opihi_flt *M3  =  V3[0].vector[0].elements.Flt;			\
+      opihi_flt *out = OUT[0].vector[0].elements.Flt;			\
+      for (i = 0; i < Nx; i++, out++, M1++, M2++, M3++) {		\
+  	*out = OP;							\
+      }									\
+      break;								\
+    }									\
+    if ((V1->vector->type == OPIHI_INT) && (V2->vector->type == OPIHI_FLT) && (V3->vector->type == OPIHI_INT)) { \
+      CopyVector (OUT[0].vector, V2[0].vector);				\
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;			\
+      opihi_flt *M2  =  V2[0].vector[0].elements.Flt;			\
+      opihi_int *M3  =  V3[0].vector[0].elements.Int;			\
+      opihi_flt *out = OUT[0].vector[0].elements.Flt;			\
+      for (i = 0; i < Nx; i++, out++, M1++, M2++, M3++) {		\
+  	*out = OP;							\
+      }									\
+      break;								\
+    }									\
+    if ((V1->vector->type == OPIHI_INT) && (V2->vector->type == OPIHI_INT) && (V3->vector->type == OPIHI_FLT)) { \
+      CopyVector (OUT[0].vector, V3[0].vector);				\
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;			\
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;			\
+      opihi_flt *M3  =  V3[0].vector[0].elements.Flt;			\
+      opihi_flt *out = OUT[0].vector[0].elements.Flt;			\
+      for (i = 0; i < Nx; i++, out++, M1++, M2++, M3++) {		\
+  	*out = OP;							\
+      }									\
+      break;								\
+    }									\
     if ((V1->vector->type == OPIHI_INT) && (V2->vector->type == OPIHI_INT) && (V3->vector->type == OPIHI_INT)) { \
       CopyVector (OUT[0].vector, V1[0].vector);				\
@@ -127,11 +193,10 @@
 int MMM_trinary (StackVar *OUT, StackVar *V1, StackVar *V2, StackVar *V3, char *op) {
 
-  int i, Nx, Ny;
+  int i;
   float *out, *M1, *M2, *M3;
   char line[512]; // this is only used to report an error 
   
-  Nx = V1[0].buffer[0].matrix.Naxis[0];
-  Ny = V1[0].buffer[0].matrix.Naxis[1];
-
+  int Npix = gfits_npix_matrix (&V1[0].buffer[0].matrix);
+  
   if (V1[0].type == ST_MATRIX_TMP) {  /** use V1 as temp buffer **/
     OUT[0].buffer = V1[0].buffer;
@@ -154,5 +219,5 @@
 
 # define MMM_FUNC(OP)					\
-  for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++, M3++) {	\
+  for (i = 0; i < Npix; i++, out++, M1++, M2++, M3++) {	\
     *out = OP;						\
   }							\
@@ -531,5 +596,5 @@
 }
 
-// the vector is applied to each column
+// the vector is applied to each column (currently only valid for 2D matrix)
 int MV_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {
 
@@ -626,5 +691,5 @@
 }
 
-// the vector is applied to each row
+// the vector is applied to each row (currently only valid for 2D matrix)
 int VM_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {
 
@@ -722,11 +787,10 @@
 int MM_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {
 
-  int i, Nx, Ny;
+  int i;
   float *out, *M1, *M2;
   char line[512]; // this is only used to report an error 
   
-  Nx = V1[0].buffer[0].matrix.Naxis[0];
-  Ny = V1[0].buffer[0].matrix.Naxis[1];
-
+  int Npix = gfits_npix_matrix (&V1[0].buffer[0].matrix);
+  
   if (V1[0].type == ST_MATRIX_TMP) {  /** use V1 as temp buffer **/
     OUT[0].buffer = V1[0].buffer;
@@ -748,5 +812,5 @@
 
 # define MM_FUNC(OP)					\
-  for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) {	\
+  for (i = 0; i < Npix; i++, out++, M1++, M2++) {	\
     *out = OP;						\
   }							\
@@ -803,9 +867,8 @@
 int MS_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {
 
-  int i, Nx, Ny;
+  int i;
   char line[512]; // this is only used to report an error 
   
-  Nx = V1[0].buffer[0].matrix.Naxis[0];
-  Ny = V1[0].buffer[0].matrix.Naxis[1];
+  int Npix = gfits_npix_matrix (&V1[0].buffer[0].matrix);
 
   /* if possible, use V1 as temp buffer, otherwise create new one */
@@ -825,5 +888,5 @@
     if (V2->type == ST_SCALAR_FLT)  {				\
       opihi_flt M2 = V2[0].FltValue;			\
-      for (i = 0; i < Nx*Ny; i++, out++, M1++) {	\
+      for (i = 0; i < Npix; i++, out++, M1++) {	\
 	*out = OP;					\
       }							\
@@ -832,5 +895,5 @@
     if (V2->type == ST_SCALAR_INT)  {				\
       opihi_int M2 = V2[0].IntValue;			\
-      for (i = 0; i < Nx*Ny; i++, out++, M1++) {	\
+      for (i = 0; i < Npix; i++, out++, M1++) {	\
 	*out = OP;					\
       }							\
@@ -881,10 +944,9 @@
 int SM_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {
 
-  int i, Nx, Ny;
+  int i;
   char line[512]; // this is only used to report an error 
   
-  Nx = V2[0].buffer[0].matrix.Naxis[0];
-  Ny = V2[0].buffer[0].matrix.Naxis[1];
-
+  int Npix = gfits_npix_matrix (&V2[0].buffer[0].matrix);
+  
   if (V2[0].type == ST_MATRIX_TMP) {  /* V2[0] is NOT temporary, we can't use it for storage */
     OUT[0].buffer = V2[0].buffer;
@@ -902,5 +964,5 @@
     if (V1->type == ST_SCALAR_FLT)  {				\
       opihi_flt M1 = V1[0].FltValue;			\
-      for (i = 0; i < Nx*Ny; i++, out++, M2++) {	\
+      for (i = 0; i < Npix; i++, out++, M2++) {	\
 	*out = OP;					\
       }							\
@@ -909,5 +971,5 @@
     if (V1->type == ST_SCALAR_INT)  {				\
       opihi_int M1 = V1[0].IntValue;			\
-      for (i = 0; i < Nx*Ny; i++, out++, M2++) {	\
+      for (i = 0; i < Npix; i++, out++, M2++) {	\
 	*out = OP;					\
       }							\
@@ -1240,5 +1302,6 @@
   if (!strcmp (op, "xramp"))  V_FUNC(i, ST_SCALAR_INT);
   if (!strcmp (op, "yramp"))  V_FUNC(0, ST_SCALAR_INT);
-  /* xramp and yramp above only make sense for matrices. for vectors, xramp = ramp, yramp = zero */
+  if (!strcmp (op, "zramp"))  V_FUNC(0, ST_SCALAR_INT);
+  /* xramp, yramp, zramp above only make sense for matrices. for vectors, xramp = ramp, yramp = zero */
 
 # undef V_FUNC
@@ -1259,10 +1322,9 @@
 int M_unary (StackVar *OUT, StackVar *V1, char *op) {
 
-  int i, j, Nx, Ny;
+  int i, j, k;
   float *out, *M1;
   
-  Nx = V1[0].buffer[0].matrix.Naxis[0];
-  Ny = V1[0].buffer[0].matrix.Naxis[1];
-
+  int Npix = gfits_npix_matrix (&V1[0].buffer[0].matrix);
+  
   if (V1[0].type == ST_MATRIX_TMP) {
     OUT[0].buffer = V1[0].buffer;
@@ -1277,56 +1339,78 @@
 
   if (!strcmp (op, "="))     { }
-  if (!strcmp (op, "abs"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = fabs(*M1);         }}
-  if (!strcmp (op, "int"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = (opihi_flt)(long long)(*M1); }}
-
-  if (!strcmp (op, "floor")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = floor (*M1); }}
-  if (!strcmp (op, "ceil"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = ceil (*M1); }}
-  // if (!strcmp (op, "rint"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = nearbyint (*M1); }}
-
-  if (!strcmp (op, "exp"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = exp(*M1);          }}
-  if (!strcmp (op, "ten"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = pow(10.0,*M1);     }}
-  if (!strcmp (op, "log"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = log10(*M1);        }}
-  if (!strcmp (op, "ln"))    { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = log(*M1);          }}
-  if (!strcmp (op, "sqrt"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = sqrt(*M1);         }}
-  if (!strcmp (op, "erf"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = erf(*M1);          }}
-
-  if (!strcmp (op, "sinh"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = sinh(*M1);         }}
-  if (!strcmp (op, "cosh"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = cosh(*M1);         }}
-  if (!strcmp (op, "asinh")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = asinh(*M1);        }}
-  if (!strcmp (op, "acosh")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = acosh(*M1);        }}
-  if (!strcmp (op, "lgamma")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = lgamma(*M1);      }}
-
-  if (!strcmp (op, "sin"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = sin(*M1);          }}
-  if (!strcmp (op, "cos"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = cos(*M1);          }}
-  if (!strcmp (op, "tan"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = tan(*M1);          }}
-  if (!strcmp (op, "dsin"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = sin(*M1*RAD_DEG);  }}
-  if (!strcmp (op, "dcos"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = cos(*M1*RAD_DEG);  }}
-  if (!strcmp (op, "dtan"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = tan(*M1*RAD_DEG);  }}
-  if (!strcmp (op, "asin"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = asin(*M1);         }}
-  if (!strcmp (op, "acos"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = acos(*M1);         }}
-  if (!strcmp (op, "atan"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = atan(*M1);         }}
-  if (!strcmp (op, "dasin")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = asin(*M1)*DEG_RAD; }}
-  if (!strcmp (op, "dacos")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = acos(*M1)*DEG_RAD; }}
-  if (!strcmp (op, "datan")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = atan(*M1)*DEG_RAD; }}
-  if (!strcmp (op, "not"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = !(*M1);            }}
-  if (!strcmp (op, "--"))    { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = -(*M1);            }}
-  if (!strcmp (op, "rnd"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = drand48();         }}
-  if (!strcmp (op, "ramp"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = i;                 }}
-  if (!strcmp (op, "zero"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = 0;                 }}
-  if (!strcmp (op, "isinf")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = !finite(*M1);      }}
-  if (!strcmp (op, "isnan")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = isnan(*M1);        }}
-
-  /* xrm and yrm only make sense in for matrices. see special meaning for vectors */
+  if (!strcmp (op, "abs"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = fabs(*M1);         }}
+  if (!strcmp (op, "int"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = (opihi_flt)(long long)(*M1); }}
+
+  if (!strcmp (op, "floor")) { for (i = 0; i < Npix; i++, out++, M1++) { *out = floor (*M1); }}
+  if (!strcmp (op, "ceil"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = ceil (*M1); }}
+  // if (!strcmp (op, "rint"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = nearbyint (*M1); }}
+
+  if (!strcmp (op, "exp"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = exp(*M1);          }}
+  if (!strcmp (op, "ten"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = pow(10.0,*M1);     }}
+  if (!strcmp (op, "log"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = log10(*M1);        }}
+  if (!strcmp (op, "ln"))    { for (i = 0; i < Npix; i++, out++, M1++) { *out = log(*M1);          }}
+  if (!strcmp (op, "sqrt"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = sqrt(*M1);         }}
+  if (!strcmp (op, "erf"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = erf(*M1);          }}
+
+  if (!strcmp (op, "sinh"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = sinh(*M1);         }}
+  if (!strcmp (op, "cosh"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = cosh(*M1);         }}
+  if (!strcmp (op, "asinh")) { for (i = 0; i < Npix; i++, out++, M1++) { *out = asinh(*M1);        }}
+  if (!strcmp (op, "acosh")) { for (i = 0; i < Npix; i++, out++, M1++) { *out = acosh(*M1);        }}
+  if (!strcmp (op, "lgamma")) { for (i = 0; i < Npix; i++, out++, M1++) { *out = lgamma(*M1);      }}
+
+  if (!strcmp (op, "sin"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = sin(*M1);          }}
+  if (!strcmp (op, "cos"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = cos(*M1);          }}
+  if (!strcmp (op, "tan"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = tan(*M1);          }}
+  if (!strcmp (op, "dsin"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = sin(*M1*RAD_DEG);  }}
+  if (!strcmp (op, "dcos"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = cos(*M1*RAD_DEG);  }}
+  if (!strcmp (op, "dtan"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = tan(*M1*RAD_DEG);  }}
+  if (!strcmp (op, "asin"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = asin(*M1);         }}
+  if (!strcmp (op, "acos"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = acos(*M1);         }}
+  if (!strcmp (op, "atan"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = atan(*M1);         }}
+  if (!strcmp (op, "dasin")) { for (i = 0; i < Npix; i++, out++, M1++) { *out = asin(*M1)*DEG_RAD; }}
+  if (!strcmp (op, "dacos")) { for (i = 0; i < Npix; i++, out++, M1++) { *out = acos(*M1)*DEG_RAD; }}
+  if (!strcmp (op, "datan")) { for (i = 0; i < Npix; i++, out++, M1++) { *out = atan(*M1)*DEG_RAD; }}
+  if (!strcmp (op, "not"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = !(*M1);            }}
+  if (!strcmp (op, "--"))    { for (i = 0; i < Npix; i++, out++, M1++) { *out = -(*M1);            }}
+  if (!strcmp (op, "rnd"))   { for (i = 0; i < Npix; i++, out++, M1++) { *out = drand48();         }}
+  if (!strcmp (op, "ramp"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = i;                 }}
+  if (!strcmp (op, "zero"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = 0;                 }}
+  if (!strcmp (op, "isinf")) { for (i = 0; i < Npix; i++, out++, M1++) { *out = !finite(*M1);      }}
+  if (!strcmp (op, "isnan")) { for (i = 0; i < Npix; i++, out++, M1++) { *out = isnan(*M1);        }}
+
+  /* xrm and yrm only make sense for 2D matrices. see special meaning for vectors */
   if (!strcmp (op, "xramp")) {
-    for (j = 0; j < Ny; j++) {
-      for (i = 0; i < Nx; i++, out++, M1++) {
-	*out = i;
+    int Nx = V1[0].buffer[0].matrix.Naxis[0];
+    int Ny = V1[0].buffer[0].matrix.Naxis[1];
+    int Nz = MAX (1, V1[0].buffer[0].matrix.Naxis[2]);
+    for (k = 0; k < Nz; k++) {
+      for (j = 0; j < Ny; j++) {
+	for (i = 0; i < Nx; i++, out++, M1++) {
+	  *out = i;
+	}
       }
     }
   }
   if (!strcmp (op, "yramp")) {
-    for (j = 0; j < Ny; j++) {
-      for (i = 0; i < Nx; i++, out++, M1++) {
-	*out = j;
+    int Nx = V1[0].buffer[0].matrix.Naxis[0];
+    int Ny = V1[0].buffer[0].matrix.Naxis[1];
+    int Nz = MAX (1, V1[0].buffer[0].matrix.Naxis[2]);
+    for (k = 0; k < Nz; k++) {
+      for (j = 0; j < Ny; j++) {
+	for (i = 0; i < Nx; i++, out++, M1++) {
+	  *out = j;
+	}
+      }
+    }
+  }
+  if (!strcmp (op, "zramp")) {
+    int Nx = V1[0].buffer[0].matrix.Naxis[0];
+    int Ny = V1[0].buffer[0].matrix.Naxis[1];
+    int Nz = MAX (1, V1[0].buffer[0].matrix.Naxis[2]);
+    for (k = 0; k < Nz; k++) {
+      for (j = 0; j < Ny; j++) {
+	for (i = 0; i < Nx; i++, out++, M1++) {
+	  *out = k;
+	}
       }
     }
Index: trunk/Ohana/src/opihi/test/ellipse.sh
===================================================================
--- trunk/Ohana/src/opihi/test/ellipse.sh	(revision 38062)
+++ trunk/Ohana/src/opihi/test/ellipse.sh	(revision 38062)
@@ -0,0 +1,122 @@
+
+macro ShapeToAxes
+  if ($0 != 7)
+    echo "USAGE: ShapeToAxes (sxx) (sxy) (syy) (major) (minor) (theta)"
+    echo " theta is returned in degrees"
+    break
+  end
+
+  # I need the concept of a local vector...
+  # I should be able to test if these are vectors or scalars
+  set _sxx = $1
+  set _sxy = $2
+  set _syy = $3
+
+  set f1 = _syy^-2 + _sxx^-2
+  set f2 = _syy^-2 - _sxx^-2
+  set f3 = sqrt(f2^2 + 4*_sxy^2)
+
+  set _minor = sqrt (2.0 / (f1 + f3))
+
+  # this returns theta in degrees
+  # @ == atan2 -- I should really replace this..
+  set _theta = -0.5 * (+2.0*_sxy @ f2)
+  set aratio2 = (f1 - f3) / (f1 + f3)
+
+  # I can test here if aratio2 is too large/small/nan
+  set _major = sqrt (2.0 / (f1 - f3))
+
+  set $4 = _major
+  set $5 = _minor
+  set $6 = _theta
+end
+
+macro AxesToShapeAlt
+  if ($0 != 7)
+    echo "USAGE: ShapeToAxes (major) (minor) (theta) (sxx) (sxy) (syy)"
+    echo " theta is supplied in degrees"
+    break
+  end
+
+  set _major = $1
+  set _minor = $2
+  set _theta = -1*$3
+
+  set tc2 = dcos(_theta)^2
+  set ts2 = dsin(_theta)^2
+
+  set f1 = tc2*_major^-2 + ts2*_minor^-2
+  set f2 = ts2*_major^-2 + tc2*_minor^-2
+  set f3 = _minor^-2 - _major^-2
+
+  set $4 = +1.0 / sqrt(f1);
+  set $6 = +1.0 / sqrt(f2);
+  set $5 = 0.5*f3*dsin(2*_theta);
+end
+
+macro AxesToShape
+  if ($0 != 7)
+    echo "USAGE: ShapeToAxes (major) (minor) (theta) (sxx) (sxy) (syy)"
+    echo " theta is supplied in degrees"
+    break
+  end
+
+  set _major = $1
+  set _minor = $2
+  set _theta = $3
+
+  set f1 = _minor^-2 + _major^-2
+  set f2 = _minor^-2 - _major^-2
+
+  set sxr = 0.5*f1 - 0.5*f2*dcos(2*_theta);
+  set syr = 0.5*f1 + 0.5*f2*dcos(2*_theta);
+
+  set $4 = +1.0 / sqrt(sxr);
+  set $6 = +1.0 / sqrt(syr);
+  set $5 = -0.5*f2*dsin(2*_theta);
+end
+
+macro test1
+
+  # create a series of ellipses with theta spinning and AR = 2.0
+  $DX = 201
+
+  create theta 0 360 10
+  set major = 20.0 + zero(theta)
+  set minor = 10.0 + zero(theta)
+
+  # AxesToShape major minor theta sxx sxy syy
+  AxesToShape major minor theta sxx sxy syy
+
+  mcreate mosaic {6*($DX + 2)} {6*($DX + 2)} 
+  mcreate base $DX $DX
+  set x = xramp(base) - int(0.5*$DX)
+  set y = yramp(base) - int(0.5*$DX)
+  
+  for i 0 theta[]
+    set r = 0.5*(x/sxx[$i])^2 + 0.5*(y/syy[$i])^2 + sxy[$i]*x*y
+    set object = exp(-r)
+    $ix = $i % 6
+    $iy = int($i / 6)
+    extract object mosaic 0 0 $DX $DX {$ix * ($DX + 1)} {$iy * ($DX + 1)}  {6*($DX + 2)} {6*($DX + 2)} 
+  end    
+
+  tv mosaic -0.01 0.5
+  lim -image
+  clear; box
+
+  for i 0 theta[]
+    $ix = $i % 6
+    $iy = int($i / 6)
+    
+    $dX = $ix * ($DX + 1) + int(0.5*$DX)
+    $dY = $iy * ($DX + 1) + int(0.5*$DX)
+
+    create T 0 360 0.1
+    set Xr = major[$i]*dcos(T)
+    set Yr = minor[$i]*dsin(T)
+    set Xo = Xr*dcos(theta[$i]) - Yr*dsin(theta[$i]) + $dX
+    set Yo = Yr*dcos(theta[$i]) + Xr*dsin(theta[$i]) + $dY
+    plot Xo Yo -x 0 -c red
+  end    
+end