Index: /branches/eam_branches/ipp-20150112/Ohana/src/opihi/lib.shell/check_stack.c
===================================================================
--- /branches/eam_branches/ipp-20150112/Ohana/src/opihi/lib.shell/check_stack.c	(revision 37856)
+++ /branches/eam_branches/ipp-20150112/Ohana/src/opihi/lib.shell/check_stack.c	(revision 37857)
@@ -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: /branches/eam_branches/ipp-20150112/Ohana/src/opihi/lib.shell/stack_math.c
===================================================================
--- /branches/eam_branches/ipp-20150112/Ohana/src/opihi/lib.shell/stack_math.c	(revision 37856)
+++ /branches/eam_branches/ipp-20150112/Ohana/src/opihi/lib.shell/stack_math.c	(revision 37857)
@@ -127,11 +127,14 @@
 int MMM_trinary (StackVar *OUT, StackVar *V1, StackVar *V2, StackVar *V3, char *op) {
 
-  int i, Nx, Ny;
+  int i, Npix;
   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];
-
+  Npix = 1;
+  for (i = 0; i < V1[0].buffer[0].matrix.Naxes; i++) {
+    if (V1[0].buffer[0].matrix.Naxis[i] == 0) break;
+    Npix *= V1[0].buffer[0].matrix.Naxis[i];
+  }
+  
   if (V1[0].type == ST_MATRIX_TMP) {  /** use V1 as temp buffer **/
     OUT[0].buffer = V1[0].buffer;
@@ -154,5 +157,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 +534,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 +629,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 +725,14 @@
 int MM_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {
 
-  int i, Nx, Ny;
+  int i, Npix;
   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];
-
+  Npix = 1;
+  for (i = 0; i < V1[0].buffer[0].matrix.Naxes; i++) {
+    if (V1[0].buffer[0].matrix.Naxis[i] == 0) break;
+    Npix *= V1[0].buffer[0].matrix.Naxis[i];
+  }
+  
   if (V1[0].type == ST_MATRIX_TMP) {  /** use V1 as temp buffer **/
     OUT[0].buffer = V1[0].buffer;
@@ -748,5 +754,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 +809,12 @@
 int MS_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {
 
-  int i, Nx, Ny;
+  int i, Npix;
   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];
+  Npix = 1;
+  for (i = 0; i < V1[0].buffer[0].matrix.Naxes; i++) {
+    if (V1[0].buffer[0].matrix.Naxis[i] == 0) break;
+    Npix *= V1[0].buffer[0].matrix.Naxis[i];
+  }
 
   /* if possible, use V1 as temp buffer, otherwise create new one */
@@ -825,5 +834,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 +841,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 +890,13 @@
 int SM_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {
 
-  int i, Nx, Ny;
+  int i, Npix;
   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];
-
+  Npix = 1;
+  for (i = 0; i < V2[0].buffer[0].matrix.Naxes; i++) {
+    if (V2[0].buffer[0].matrix.Naxis[i] == 0) break;
+    Npix *= V2[0].buffer[0].matrix.Naxis[i];
+  }
+  
   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 +914,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 +921,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;					\
       }							\
@@ -1259,10 +1271,13 @@
 int M_unary (StackVar *OUT, StackVar *V1, char *op) {
 
-  int i, j, Nx, Ny;
+  int i, j, Npix;
   float *out, *M1;
   
-  Nx = V1[0].buffer[0].matrix.Naxis[0];
-  Ny = V1[0].buffer[0].matrix.Naxis[1];
-
+  Npix = 1;
+  for (i = 0; i < V1[0].buffer[0].matrix.Naxes; i++) {
+    if (V1[0].buffer[0].matrix.Naxis[i] == 0) break;
+    Npix *= V1[0].buffer[0].matrix.Naxis[i];
+  }
+  
   if (V1[0].type == ST_MATRIX_TMP) {
     OUT[0].buffer = V1[0].buffer;
@@ -1277,46 +1292,48 @@
 
   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")) {
+    int Nx = V1[0].buffer[0].matrix.Naxis[0];
+    int Ny = V1[0].buffer[0].matrix.Naxis[1];
     for (j = 0; j < Ny; j++) {
       for (i = 0; i < Nx; i++, out++, M1++) {
@@ -1326,4 +1343,6 @@
   }
   if (!strcmp (op, "yramp")) {
+    int Nx = V1[0].buffer[0].matrix.Naxis[0];
+    int Ny = V1[0].buffer[0].matrix.Naxis[1];
     for (j = 0; j < Ny; j++) {
       for (i = 0; i < Nx; i++, out++, M1++) {