Index: trunk/Ohana/src/libfits/include/gfitsio.h
===================================================================
--- trunk/Ohana/src/libfits/include/gfitsio.h	(revision 41175)
+++ trunk/Ohana/src/libfits/include/gfitsio.h	(revision 41176)
@@ -259,5 +259,5 @@
 int     gfits_read_table               PROTO((char *filename, FTable *ftable)); 
 int     gfits_set_bintable_column      PROTO((Header *header, FTable *table, char *label, void *data, off_t Nrow));
-int     gfits_set_bintable_column_reformat PROTO((Header *header, FTable *table, char *label, char *intype, void *data, off_t Nrow, char nativeOrder));
+int     gfits_set_bintable_column_reformat PROTO((Header *header, FTable *table, char *label, char *intype, void *data, off_t Nrow, int element, char nativeOrder));
 int     gfits_set_table_column         PROTO((Header *header, FTable *table, char *label, void *data, off_t Nrow));
 int     gfits_table_column             PROTO((FTable *ftable, char *field, char *mode,...)) OHANA_FORMAT(printf, 3, 4);
Index: trunk/Ohana/src/libfits/table/F_define_column.c
===================================================================
--- trunk/Ohana/src/libfits/table/F_define_column.c	(revision 41175)
+++ trunk/Ohana/src/libfits/table/F_define_column.c	(revision 41176)
@@ -56,4 +56,5 @@
   char type[64], field[64];
   
+  // this call supported multiple columns per named field
   if (!gfits_bintable_format (format, type, &Nval, &Nbytes)) return (FALSE);
   
Index: trunk/Ohana/src/libfits/table/F_set_column.c
===================================================================
--- trunk/Ohana/src/libfits/table/F_set_column.c	(revision 41175)
+++ trunk/Ohana/src/libfits/table/F_set_column.c	(revision 41176)
@@ -253,5 +253,5 @@
 /***********************/
 // convert the input data array (of the specified intype) to the desired table data type. swap unless nativeOrder is requested
-int gfits_set_bintable_column_reformat (Header *header, FTable *table, char *label, char *intype, void *data, off_t Nrow, char nativeOrder) {
+int gfits_set_bintable_column_reformat (Header *header, FTable *table, char *label, char *intype, void *data, off_t Nrow, int element, char nativeOrder) {
 
   off_t Nx, Ny;
@@ -288,4 +288,8 @@
 
   if (!gfits_bintable_format (format, outtype, &Nval, &NbytesOut)) return (FALSE);
+  if (element >= Nval) {
+    fprintf (stderr, "programming error: element >= Nval for field %s, format %s\n", label, format);
+    return FALSE;
+  }
   
   /* check existing table dimensions */
@@ -307,9 +311,13 @@
   }
 
+  // NOTE: we are inserting a single column into the output
+  // table.  If the output field is multi-value, the resulting
+  // array is inserted into the appropriate bytes in that output field
+
   /* make duplicate of data with correct type
      byte swap and Bzero/Bscale */
-  ALLOCATE (array, char, NbytesOut*Nval*Nrow);
+  ALLOCATE (array, char, NbytesOut*Nrow);
   Pin = data;
-  Pout = array;
+  Pout = array; // 
 
   // # define ASSIGN_DATA(ITYPE,INAME,ISIZE,OTYPE,ONAME)
@@ -321,5 +329,5 @@
   if (!strcmp (outtype, OUTNAME) && !strcmp (intype, INNAME)) {		\
     int NbytesIn = NBYTES_IN;						\
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { \
+    for (i = 0; i < Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) { \
       if (directCopy) {							\
 	*(OUTTYPE *)Pout = *(INTYPE *)Pin;				\
@@ -385,425 +393,4 @@
   SET_VALUES("double",   double, "double", double, SWAP_DBLE, 8);
 
-# if (0)
-  /** input == char **/
-  if (!strcmp (outtype, "char") && !strcmp (intype, "char")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) {
-      if (directCopy) { *(char *)Pout = *(char *)Pin; } else {
-	*(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
-      }
-    }
-  }
-  if (!strcmp (outtype, "byte") && !strcmp (intype, "char")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) {
-      if (directCopy) { *(char *)Pout = *(char *)Pin; } else {
-	*(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
-      }
-    }
-  }
-  if (!strcmp (outtype, "short") && !strcmp (intype, "char")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      if (directCopy) { *(short *)Pout = *(char *)Pin; } else {
-	*(short *)Pout = (*(char *)Pin - Bzero) / Bscale;
-      }
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_BYTE; }
-# endif
-    }  
-  }
-  if (!strcmp (outtype, "int") && !strcmp (intype, "char")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int *)Pout = (*(char *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "char")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int64_t *)Pout = (*(char *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "float") && !strcmp (intype, "char")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(float *)Pout = (*(char *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "double") && !strcmp (intype, "char")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(double *)Pout = (*(char *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-
-  /** input == byte **/
-  if (!strcmp (outtype, "char") && !strcmp (intype, "byte")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) {
-      *(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "byte") && !strcmp (intype, "byte")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) {
-      *(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "short") && !strcmp (intype, "byte")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(short *)Pout = (*(char *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_BYTE; }
-# endif
-    }  
-  }
-  if (!strcmp (outtype, "int") && !strcmp (intype, "byte")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int *)Pout = (*(char *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "byte")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int64_t *)Pout = (*(char *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "float") && !strcmp (intype, "byte")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(float *)Pout = (*(char *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "double") && !strcmp (intype, "byte")) {
-    int NbytesIn = 1;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(double *)Pout = (*(char *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-
-  /** input == short **/
-  if (!strcmp (outtype, "char") && !strcmp (intype, "short")) {
-    int NbytesIn = 2;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(short *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "byte") && !strcmp (intype, "short")) {
-    int NbytesIn = 2;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(short *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "short") && !strcmp (intype, "short")) {
-    int NbytesIn = 2;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(short *)Pout = (*(short *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_BYTE; }
-# endif
-    }  
-  }
-  if (!strcmp (outtype, "int") && !strcmp (intype, "short")) {
-    int NbytesIn = 2;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int *)Pout = (*(short *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "short")) {
-    int NbytesIn = 2;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int64_t *)Pout = (*(short *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "float") && !strcmp (intype, "short")) {
-    int NbytesIn = 2;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(float *)Pout = (*(short *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "double") && !strcmp (intype, "short")) {
-    int NbytesIn = 2;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(double *)Pout = (*(short *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-
-  /** input == int **/
-  if (!strcmp (outtype, "char") && !strcmp (intype, "int")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(int *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "byte") && !strcmp (intype, "int")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(int *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "short") && !strcmp (intype, "int")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(short *)Pout = (*(int *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_BYTE; }
-# endif
-    }  
-  }
-  if (!strcmp (outtype, "int") && !strcmp (intype, "int")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int *)Pout = (*(int *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "int")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int64_t *)Pout = (*(int *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "float") && !strcmp (intype, "int")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(float *)Pout = (*(int *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "double") && !strcmp (intype, "int")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(double *)Pout = (*(int *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-
-  /** input == int64_t **/
-  if (!strcmp (outtype, "char") && !strcmp (intype, "int64_t")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "byte") && !strcmp (intype, "int64_t")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "short") && !strcmp (intype, "int64_t")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(short *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_BYTE; }
-# endif
-    }  
-  }
-  if (!strcmp (outtype, "int") && !strcmp (intype, "int64_t")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "int64_t")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int64_t *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "float") && !strcmp (intype, "int64_t")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(float *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "double") && !strcmp (intype, "int64_t")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(double *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-
-  /** input == float **/
-  if (!strcmp (outtype, "char") && !strcmp (intype, "float")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(float *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "byte") && !strcmp (intype, "float")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(float *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "short") && !strcmp (intype, "float")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(short *)Pout = (*(float *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_BYTE; }
-# endif
-    }  
-  }
-  if (!strcmp (outtype, "int") && !strcmp (intype, "float")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int *)Pout = (*(float *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "float")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int64_t *)Pout = (*(float *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "float") && !strcmp (intype, "float")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(float *)Pout = (*(float *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "double") && !strcmp (intype, "float")) {
-    int NbytesIn = 4;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(double *)Pout = (*(float *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-
-  /** input == double **/
-  if (!strcmp (outtype, "char") && !strcmp (intype, "double")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(double *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "byte") && !strcmp (intype, "double")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(char *)Pout = (*(double *)Pin - Bzero) / Bscale;
-    }
-  }
-  if (!strcmp (outtype, "short") && !strcmp (intype, "double")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(short *)Pout = (*(double *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_BYTE; }
-# endif
-    }  
-  }
-  if (!strcmp (outtype, "int") && !strcmp (intype, "double")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int *)Pout = (*(double *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "double")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(int64_t *)Pout = (*(double *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "float") && !strcmp (intype, "double")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(float *)Pout = (*(double *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_WORD; }
-# endif
-    }
-  }
-  if (!strcmp (outtype, "double") && !strcmp (intype, "double")) {
-    int NbytesIn = 8;
-    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
-      *(double *)Pout = (*(double *)Pin - Bzero) / Bscale;
-# ifdef BYTE_SWAP
-      if (!nativeOrder) { SWAP_DBLE; }
-# endif
-    }
-  }
-# endif
-
   /* check array space */
   if (Nx*Ny < Nx*(Nrow - 1) + Nstart + Nval*NbytesOut) {
@@ -813,8 +400,8 @@
 
   /* insert bytes from array into appropriate section of buffer */
-  Pout = table[0].buffer + Nstart;
+  Pout = table[0].buffer + Nstart + element*NbytesOut;
   Pin  = array;
-  for (i = 0; i < Nrow; i++, Pout += Nx, Pin += Nval*NbytesOut) {
-    memcpy (Pout, Pin, Nval*NbytesOut);
+  for (i = 0; i < Nrow; i++, Pout += Nx, Pin += NbytesOut) {
+    memcpy (Pout, Pin, NbytesOut);
   }
 
Index: trunk/Ohana/src/libohana/src/isolate_elements.c
===================================================================
--- trunk/Ohana/src/libohana/src/isolate_elements.c	(revision 41175)
+++ trunk/Ohana/src/libohana/src/isolate_elements.c	(revision 41176)
@@ -186,9 +186,4 @@
   if (!strncmp (c, "+",  1)) return (TRUE);
   if (!strncmp (c, "-",  1)) return (TRUE);
-  if (!strncmp (c, "&",  1)) return (TRUE);
-  if (!strncmp (c, "|",  1)) return (TRUE);
-
-  if (!strncmp (c, "<",  1)) return (TRUE);
-  if (!strncmp (c, ">",  1)) return (TRUE);
 
   if (!strncmp (c, "(",  1)) return (TRUE);
@@ -203,4 +198,11 @@
   if (!strncmp (c, "<=", 2)) return (TRUE);
   if (!strncmp (c, ">=", 2)) return (TRUE);
+  if (!strncmp (c, "<~", 2)) return (TRUE);
+  if (!strncmp (c, "~>", 2)) return (TRUE);
+
+  if (!strncmp (c, "<",  1)) return (TRUE);
+  if (!strncmp (c, ">",  1)) return (TRUE);
+  if (!strncmp (c, "&",  1)) return (TRUE);
+  if (!strncmp (c, "|",  1)) return (TRUE);
 
   return (FALSE);
@@ -218,4 +220,6 @@
   if (!strncmp (c, "<=", 2)) return (TRUE);
   if (!strncmp (c, ">=", 2)) return (TRUE);
+  if (!strncmp (c, "<~", 2)) return (TRUE);
+  if (!strncmp (c, "~>", 2)) return (TRUE);
 
   return (FALSE);
Index: trunk/Ohana/src/opihi/cmd.basic/test/stackmath.sh
===================================================================
--- trunk/Ohana/src/opihi/cmd.basic/test/stackmath.sh	(revision 41176)
+++ trunk/Ohana/src/opihi/cmd.basic/test/stackmath.sh	(revision 41176)
@@ -0,0 +1,128 @@
+
+input tap.sh
+
+macro go
+ plan {14 + 14}
+ test1
+ test2
+ done
+end
+
+# test unary scalar ops
+macro test1
+ $x = 5;                 ok ($x == 5) "scalar equality"
+ $x = abs(-3);           ok ($x == 3) "scalar abs"
+ $x = int(5.23);         ok ($x == 5) "scalar int"
+ $x = floor(7.1234);     ok ($x == 7) "scalar floor"
+ $x = floor(7.9876);     ok ($x == 7) "scalar floor"
+ $x = ceil(7.1234);      ok ($x == 8) "scalar ceil"
+ $x = ceil(7.9876);      ok ($x == 8) "scalar ceil"
+		         
+ $x = exp(1.0);          ok ($x == 2.71828182846) "scalar exp"
+ $x = ten(3.0);          ok ($x == 1000) "scalar ten"
+ $x = log(10000);        ok ($x == 4) "scalar log"
+ $x = ln(2.71828182846); ok ($x == 1) "scalar ln"
+ $x = sqrt(9);           ok ($x == 3) "scalar sqrt"
+ $x = dsin(30);          ok ($x == 0.5) "scalar dsin"
+ $x = sin($M_PI/6);      ok ($x == 0.5) "scalar sin"
+end
+
+# test unary vector ops
+macro test2
+ create x -3 3
+
+ set y = x;                 ok (y[1] == -2) "vector equality"
+ set y = abs(x);            ok (y[1] ==  2) "vector abs"
+ 
+ set y = ten(x);            ok (y[1] == 0.01) "vector ten"
+       	 		    ok (y[5] == 100) "vector ten"
+ set z = log(y);            ok (z[5] == 2) "vector log"
+       	 		    ok (z[1] == -2) "vector ten"
+
+ set y = exp(x);            ok (y[4] == 2.71828182846) "vector exp"
+       	 		    ok (y[2] == 0.367879441171) "vector exp"
+ set z = ln(y);             ok (z[2] == -1) "vector ln"
+                            ok (z[4] == 1) "vector ln"
+ create x 0 26
+ set y = sqrt(x);           ok (y[4] == 2) "vector sqrt"
+       	 		    ok (y[9] == 3) "vector sqrt"
+
+ create x 0 360 30
+ set y = dsin(x);           ok (y[1] == 0.5) "vector dsin"
+                            ok (y[3] == 1.0) "vector dsin"
+end
+
+# here are all defined unary functions
+# SVM =
+# SVM abs    
+# SVM int    
+# SVM floor  
+#     round  
+# SVM ceil   
+#     rint   
+# SVM exp    
+# SVM ten    
+# SVM log    
+# SVM ln     
+# SVM sqrt   
+# SVM erf    
+# SVM sinh   
+# SVM cosh   
+# SVM asinh  
+# SVM acosh  
+# SVM sin    
+# SVM cos    
+# SVM tan    
+# SVM dsin   
+# SVM dcos   
+# SVM dtan   
+# SVM asin   
+# SVM acos   
+# SVM atan   
+# SVM dasin  
+# SVM dacos  
+# SVM datan  
+# SVM lgamma 
+# SVM rnd    
+# SVM drnd   
+# SVM lrnd   
+# SVM mrnd   
+# SVM not    
+# SVM negate
+# SVM isinf  
+# SVM isnan  
+#  VM ramp   
+#  VM xramp  
+#  VM yramp  
+#  VM zramp  
+#  VM zero   
+
+## tests are not reporting an error when using one of the unimplemented functions
+
+# here are all defined binary functions
+
+# V "^"
+# V max     "U"
+# V min     "D"
+# V atan2   "a"
+# V datan2  "d"
+# * <-      "l"
+# * ->      "r"
+# V @       
+# V /       
+# V *       
+# V %       
+# V +       
+# V -       
+# V &       
+# V |       
+# V <       
+# V >       
+# V ==      "E"
+# V !=      "N"
+# V <=      "L"
+# V >=      "G"
+# V >>      "U"
+# V <<      "D"
+# V &&      "A"
+# V ||      "O"
Index: trunk/Ohana/src/opihi/cmd.basic/test/tap.sh
===================================================================
--- trunk/Ohana/src/opihi/cmd.basic/test/tap.sh	(revision 41176)
+++ trunk/Ohana/src/opihi/cmd.basic/test/tap.sh	(revision 41176)
@@ -0,0 +1,61 @@
+
+macro plan
+ if ($0 != 2)
+  echo "USAGE: plan (Ntest)"
+  break
+ end
+
+ $tap_Nplan = $1
+ $tap_Npass = 0
+ $tap_Nfail = 0
+ $tap_Ntest = 0
+end
+
+macro done
+ if (not($?tap_Nplan))
+  echo "ERROR: tap tests not planned"
+  break
+ end
+
+ echo "planned $tap_Nplan tests, ran $tap_Ntest"
+ echo "passed $tap_Npass tests, failed $tap_Nfail tests"
+
+ if (($tap_Nfail == 0) && ($tap_Ntest == $tap_Nplan) && ($tap_Npass == $tap_Ntest))
+   echo "PASS: all tests run and passed"
+   return 1
+ end
+
+ if ($tap_Nfail > 0)
+   echo "FAIL: $tap_Nfail tests failed"
+ end
+
+ if ($tap_Ntest != $tap_Nplan)
+   echo "FAIL: not all tests run"
+ end
+
+ if ($tap_Ntest != $tap_Npass + $tap_Nfail)
+   echo "FAIL: inconsistent count: $tap_Ntest tests != $tap_Npass pass + $tap_Nfail fail"
+ end
+end
+
+macro ok
+ if ($0 != 3)
+  echo "USAGE (condition) (label)"
+  break
+ end
+ 
+ if (not($?tap_Nplan))
+  echo "ERROR: tap tests not planned"
+  break
+ end
+ 
+ if ($1)
+  echo "OK: $2"
+  $tap_Npass ++
+ else     
+  echo "not OK: $2"
+  $tap_Nfail ++
+ end
+ $tap_Ntest ++
+end
+
Index: trunk/Ohana/src/opihi/lib.data/SplineOps.c
===================================================================
--- trunk/Ohana/src/opihi/lib.data/SplineOps.c	(revision 41175)
+++ trunk/Ohana/src/opihi/lib.data/SplineOps.c	(revision 41176)
@@ -181,7 +181,7 @@
 
   // NOTE: if we want to compress the output table, use native byte order here (last element)
-  gfits_set_bintable_column_reformat (&theader, &ftable, "X_KNOT", "double", myspline->xk, myspline->Nknots, FALSE);
-  gfits_set_bintable_column_reformat (&theader, &ftable, "Y_KNOT", "double", myspline->yk, myspline->Nknots, FALSE);
-  gfits_set_bintable_column_reformat (&theader, &ftable, "DY2_DX", "double", myspline->y2, myspline->Nknots, FALSE);
+  gfits_set_bintable_column_reformat (&theader, &ftable, "X_KNOT", "double", myspline->xk, myspline->Nknots, 0, FALSE);
+  gfits_set_bintable_column_reformat (&theader, &ftable, "Y_KNOT", "double", myspline->yk, myspline->Nknots, 0, FALSE);
+  gfits_set_bintable_column_reformat (&theader, &ftable, "DY2_DX", "double", myspline->y2, myspline->Nknots, 0, FALSE);
 
   if (!append) {
Index: trunk/Ohana/src/opihi/lib.shell/VectorIO.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/VectorIO.c	(revision 41175)
+++ trunk/Ohana/src/opihi/lib.shell/VectorIO.c	(revision 41176)
@@ -7,17 +7,28 @@
   int j;
 
-  char *tformat = NULL;
-
   Header *theader = ftable->header;
   gfits_create_table_header (theader, "BINTABLE", extname);
 
-  ALLOCATE (tformat, char, 2*Nvec);
+  // allocate an array of strings to represent the format for each output field
+  // formats include single column formats (BIJKDE) and multi-column formats (e.g., 2I)
+  // we will have no more than Nvec fields (but we can have fewer)
+  int Nfield = 0;
+  ALLOCATE_PTR (tformat, char *, Nvec);
+  ALLOCATE_PTR (Nelement, int, Nvec);
+
   if (format) {
-    // the bintable format string can only define the byte-width of each field.  valid output columns are currently:
+    // the bintable format string can defines the byte-width of each field and number of elements (columns per field).
+    // valid output columns are currently:
     // B (char), I (16 bit short), J (32 bit int), E (32 bit float), D (64 bit double).
-    // the format string is just the sequence of types, eg: LIIJEED
-    // validate the format string
+    // the format string is just the sequence of types, eg: LIIJEED.
+    // it may have spaces or integer element counts:
+    // "2D 4I EEJ"
+
+    // *** validate the format string
+
+    // as I parse each elements, if it is a digit, I need parse that value
+
     char *ptr = format;
-    for (j = 0; j < Nvec; j++) {
+    for (j = 0; j < Nvec; ) {
       while (*ptr && OHANA_WHITESPACE (*ptr)) ptr++;
       if (*ptr == 0) {
@@ -25,11 +36,33 @@
 	goto escape;
       }
+
+      // is there a leading integer?
+      char *endptr;
+      Nelement[Nfield] = strtol (ptr, &endptr, 10);
+      if (endptr == ptr) {
+	Nelement[Nfield] = 1;
+      }
+      ptr = endptr; // this should now point at the letter that is the format type
       if ((*ptr != 'B') && (*ptr != 'I') && (*ptr != 'J') && (*ptr != 'K') && (*ptr != 'D') && (*ptr != 'E')) {
 	gprint (GP_ERR, "error in binary table format %s: invalid character %c\n", format, *ptr);
 	goto escape;
       }
-      tformat[2*j + 0] = *ptr;
-      tformat[2*j + 1] = 0; // a bit sleazy : use a 2xN string to store N 1-byte strings
+
+      int Nchar = snprintf (tformat[Nfield], 0, "%d%c", Nelement[Nfield], *ptr);
+      ALLOCATE (tformat[Nfield], char, Nchar + 1);
+      int Nout = snprintf (tformat[Nfield], Nchar + 1, "%d%c", Nelement[Nfield], *ptr);
+      myAssert (Nout <= Nchar, "oops");
+
+      // tformat[2*j + 0] = *ptr;
+      // tformat[2*j + 1] = 0; // a bit sleazy : use a 2xN string to store N 1-byte strings
+
+      j += Nelement[Nfield]; // advance past 
+      if (j > Nvec) {
+	gprint (GP_ERR, "error in binary table format %s (too few vectors for listed field)\n", format);
+	goto escape;
+      }
+
       ptr ++;
+      Nfield ++;
     }
     while (*ptr && OHANA_WHITESPACE (*ptr)) ptr++;
@@ -41,7 +74,10 @@
     for (j = 0; j < Nvec; j++) {
       // if the format is not defined, just use the native byte-widths
-      tformat[2*j + 0] = (vec[j][0].type == OPIHI_FLT) ? 'D' : 'K'; // this depends on opihi_int == int64_t for Int
-      tformat[2*j + 1] = 0;
-    }
+      ALLOCATE (tformat[j], char, 2);
+      tformat[j][0] = (vec[j][0].type == OPIHI_FLT) ? 'D' : 'K'; // this depends on opihi_int == int64_t for Int
+      tformat[j][1] = 0;
+      Nelement[j] = 1;
+    }
+    Nfield = Nvec;
   }
 
@@ -49,6 +85,15 @@
   // output table (because the data goes to the named column below).  need to enforce
   // this somehow
-  for (j = 0; j < Nvec; j++) {
-    gfits_define_bintable_column (theader, &tformat[2*j], vec[j][0].name, NULL, NULL, 1.0, 0.0);
+  int ivec = 0;
+  for (j = 0; j < Nfield; j++) {
+    // XX need to loop over fields, and skip the additional vectors that are part of a field
+    // this call supported multiple columns per named field
+    gfits_define_bintable_column (theader, tformat[j], vec[ivec][0].name, NULL, NULL, 1.0, 0.0);
+    ivec += Nelement[j];
+  }
+
+  // need to free the array 
+  for (j = 0; j < Nfield; j++) {
+    free (tformat[j]);
   }
   free (tformat);
@@ -57,11 +102,18 @@
   gfits_create_table (theader, ftable);
 
+  // I need to add each vector in order, but I need to
+  // track which field it corresponds to.
+
   // add the vectors to the output array
-  for (j = 0; j < Nvec; j++) {
-    if (vec[j][0].type == OPIHI_FLT) {
-      gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "double",  vec[j][0].elements.Flt, vec[j][0].Nelements, nativeOrder);
-    } else {
-//    gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "int",     vec[j][0].elements.Int, vec[j][0].Nelements, nativeOrder);
-      gfits_set_bintable_column_reformat (theader, ftable, vec[j][0].name, "int64_t", vec[j][0].elements.Int, vec[j][0].Nelements, nativeOrder);
+  for (ivec = 0, j = 0; j < Nfield; j++) {
+    // the first vector provides the name for the field
+    char *fieldname = vec[ivec][0].name;
+    for (int k = 0; k < Nelement[j]; k++, ivec++) {
+      Vector *thisvec = vec[ivec];
+      if (thisvec->type == OPIHI_FLT) {
+	gfits_set_bintable_column_reformat (theader, ftable, fieldname, "double",  thisvec->elements.Flt, thisvec->Nelements, k, nativeOrder);
+      } else {
+	gfits_set_bintable_column_reformat (theader, ftable, fieldname, "int64_t", thisvec->elements.Int, thisvec->Nelements, k, nativeOrder);
+      }
     }
   }
Index: trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c	(revision 41175)
+++ trunk/Ohana/src/opihi/lib.shell/convert_to_RPN.c	(revision 41176)
@@ -78,4 +78,8 @@
     if (!strcmp (argv[i], "atan2"))  { type = ST_BINARY; strcpy (argv[i], "a"); goto gotit; }
     if (!strcmp (argv[i], "datan2")) { type = ST_BINARY; strcpy (argv[i], "d"); goto gotit; }
+
+    if (!strcmp (argv[i], "<~"))     { type = ST_BINARY; strcpy (argv[i], "l"); goto gotit; }
+    if (!strcmp (argv[i], "~>"))     { type = ST_BINARY; strcpy (argv[i], "r"); goto gotit; }
+
     if (!strcmp (argv[i], ","))      { type = ST_COMMA; goto gotit; }
 
@@ -99,4 +103,8 @@
     if (!strcmp (argv[i], ">>"))     { type = ST_LOGIC; strcpy (argv[i], "U"); goto gotit; }
     if (!strcmp (argv[i], "<<"))     { type = ST_LOGIC; strcpy (argv[i], "D"); goto gotit; }
+
+    /* XXX I would like to change the syntax to allow << and >> to mean bitshifts
+       but that means breaking these older values which means MIN and MAX */
+    // for now, use <- and -> to mean bitshift 
 
     if (!strcmp (argv[i], "&&"))     { type = ST_AND; strcpy (argv[i], "A"); goto gotit; }
Index: trunk/Ohana/src/opihi/lib.shell/parse.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/parse.c	(revision 41175)
+++ trunk/Ohana/src/opihi/lib.shell/parse.c	(revision 41176)
@@ -66,5 +66,5 @@
 	/* dvomath returns a new string, or NULL, with the result of the expression */
 	val = dvomath (1, &V1, &size, 0);
-	if (val == NULL) goto error;
+	if (val == NULL) { print_error (); goto error; }
 	fval += atof(val);
 	// save the result
@@ -88,5 +88,5 @@
 	/* dvomath returns a new string, or NULL, with the result of the expression */
 	val = dvomath (1, &V1, &size, 0);
-	if (val == NULL) goto error;
+	if (val == NULL) { print_error (); goto error; }
 	fval -= atof(val);
 	// save the result
Index: trunk/Ohana/src/opihi/lib.shell/stack_math.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/stack_math.c	(revision 41175)
+++ trunk/Ohana/src/opihi/lib.shell/stack_math.c	(revision 41176)
@@ -34,5 +34,5 @@
     case '?': SSS_FUNC(M1 ? M2: M3);
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined as scalar trinary op!", op[0]);
       push_error (line);
       return (FALSE);
@@ -162,5 +162,5 @@
     case '?': VVV_FUNC(*M1 ? *M2: *M3);
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined as vector trinary op!", op[0]);
       push_error (line);
       return (FALSE);
@@ -227,5 +227,5 @@
     case '?': MMM_FUNC(*M1 ? *M2: *M3);
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined as matrix trinary op!", op[0]);
       push_error (line);
       return (FALSE);
@@ -357,6 +357,48 @@
     case 'A': VV_FUNC(ST_SCALAR_INT, (*M1 && *M2) ? 1 : 0);
     case 'O': VV_FUNC(ST_SCALAR_INT, (*M1 || *M2) ? 1 : 0);
+
+    // for the bitshift operators, we have to treat the INT and FLT values differently
+    // this makes the operator incompatible with the macros used above
+    case 'l': {
+      CopyVector (OUT[0].vector, V1[0].vector);
+      if ((V1->vector->type == OPIHI_FLT) || (V2->vector->type == OPIHI_FLT)) {
+	// bitshift is not valid with float valuess
+	for (i = 0; i < Nx; i++) {
+	  OUT[0].vector[0].elements.Flt[i] = NAN;
+	}
+	break;
+      }
+      // I could just do this for all types and bitshift regardless...
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M1++, M2++) {
+  	*out = *M1 << *M2;
+      }
+      break;
+    }
+
+    case 'r': {
+      CopyVector (OUT[0].vector, V1[0].vector);
+      if ((V1->vector->type == OPIHI_FLT) || (V2->vector->type == OPIHI_FLT)) {
+	// bitshift is not valid with float valuess
+	CopyVector (OUT[0].vector, V1[0].vector);
+	for (i = 0; i < Nx; i++) {
+	  OUT[0].vector[0].elements.Flt[i] = NAN;
+	}
+	break;
+      }
+      // I could just do this for all types and bitshift regardless...
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M1++, M2++) {
+  	*out = *M1 >> *M2;
+      }
+      break;
+    }
+
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (vector OP vector)!", op[0]);
       push_error (line);
       return (FALSE);
@@ -469,6 +511,45 @@
     case 'A': SV_FUNC(ST_SCALAR_INT, (M1 && *M2) ? 1 : 0);
     case 'O': SV_FUNC(ST_SCALAR_INT, (M1 || *M2) ? 1 : 0);
+
+    // for the bitshift operators, we have to treat the INT and FLT values differently
+    // this makes the operator incompatible with the macros used above
+    case 'l': {
+      CopyVector (OUT[0].vector, V2[0].vector);
+      if ((V1->type == ST_SCALAR_FLT) || (V2->vector->type == OPIHI_FLT)) {
+	// bitshift is not valid with float valuess
+	for (i = 0; i < Nx; i++) {
+	  OUT[0].vector[0].elements.Flt[i] = NAN;
+	}
+	break;
+      }
+      opihi_int  M1  =  V1[0].IntValue;					\
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M2++) {
+  	*out = M1 << *M2;
+      }
+      break;
+    }
+
+    case 'r': {
+      CopyVector (OUT[0].vector, V2[0].vector);
+      if ((V1->type == ST_SCALAR_FLT) || (V2->vector->type == OPIHI_FLT)) {
+	// bitshift is not valid with float valuess
+	for (i = 0; i < Nx; i++) {
+	  OUT[0].vector[0].elements.Flt[i] = NAN;
+	}
+	break;
+      }
+      opihi_int  M1  =  V1[0].IntValue;					\
+      opihi_int *M2  =  V2[0].vector[0].elements.Int;
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M2++) {
+  	*out = M1 >> *M2;
+      }
+      break;
+    }
+
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (scalar OP vector)!", op[0]);
       push_error (line);
       return (FALSE);
@@ -577,6 +658,44 @@
     case 'A': VS_FUNC(ST_SCALAR_INT, (*M1 && M2) ? 1 : 0);
     case 'O': VS_FUNC(ST_SCALAR_INT, (*M1 || M2) ? 1 : 0);
+
+    case 'l': {
+      CopyVector (OUT[0].vector, V1[0].vector);
+      if ((V1->vector->type == OPIHI_FLT) || (V2->type == ST_SCALAR_FLT)) {
+	// bitshift is not valid with float valuess
+	for (i = 0; i < Nx; i++) {
+	  OUT[0].vector[0].elements.Flt[i] = NAN;
+	}
+	break;
+      }
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;
+      opihi_int  M2  =  V2[0].IntValue;					\
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M1++) {
+  	*out = *M1 << M2;
+      }
+      break;
+    }
+
+    case 'r': {
+      CopyVector (OUT[0].vector, V1[0].vector);
+      if ((V1->vector->type == OPIHI_FLT) || (V2->type == ST_SCALAR_FLT)) {
+	// bitshift is not valid with float valuess
+	CopyVector (OUT[0].vector, V1[0].vector);
+	for (i = 0; i < Nx; i++) {
+	  OUT[0].vector[0].elements.Flt[i] = NAN;
+	}
+	break;
+      }
+      opihi_int *M1  =  V1[0].vector[0].elements.Int;
+      opihi_int  M2  =  V2[0].IntValue;					\
+      opihi_int *out = OUT[0].vector[0].elements.Int;
+      for (i = 0; i < Nx; i++, out++, M1++) {
+  	*out = *M1 >> M2;
+      }
+      break;
+    }
+
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (vector OP scalar)!", op[0]);
       push_error (line);
       return (FALSE);
@@ -670,6 +789,7 @@
     case 'A': MV_FUNC((*M1 && *M2) ? 1 : 0);
     case 'O': MV_FUNC((*M1 || *M2) ? 1 : 0);
+
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (matrix OP vector)!", op[0]);
       push_error (line);
       return (FALSE);
@@ -768,5 +888,5 @@
     case 'O': VM_FUNC((*M1 || *M2) ? 1 : 0);
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (vector OP matrix)!", op[0]);
       push_error (line);
       return (FALSE);
@@ -849,5 +969,5 @@
     case 'O': MM_FUNC((*M1 || *M2) ? 1 : 0);
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (matrix OP matrix)!", op[0]);
       push_error (line);
       return (FALSE);
@@ -935,5 +1055,5 @@
     case 'O': MS_FUNC((*M1 || M2) ? 1 : 0);
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (matrix OP scalar)!", op[0]);
       push_error (line);
       return (FALSE);
@@ -1012,5 +1132,5 @@
     case 'O': SM_FUNC((M1 || *M2) ? 1 : 0);
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (scalar OP matrix)!", op[0]);
       push_error (line);
       return (FALSE);
@@ -1094,6 +1214,37 @@
     case 'A': SS_FUNC(ST_SCALAR_INT, (M1 && M2) ? 1 : 0);
     case 'O': SS_FUNC(ST_SCALAR_INT, (M1 || M2) ? 1 : 0);
+
+    // for the bitshift operators, we have to treat the INT and FLT values differently
+    // this makes the operator incompatible with the macros used above
+    case 'l': {
+      if ((V1->type == ST_SCALAR_FLT) || (V2->type == ST_SCALAR_FLT)) {
+	// bitshift is not valid with float valuess
+	OUT[0].type = ST_SCALAR_FLT;
+	OUT[0].FltValue = NAN;
+	break;
+      }
+      opihi_int M1 = V1[0].IntValue;
+      opihi_int M2 = V2[0].IntValue;
+      OUT[0].type = ST_SCALAR_INT;
+      OUT[0].IntValue = M1 << M2;
+      break;
+    }
+
+    case 'r': {
+      if ((V1->type == ST_SCALAR_FLT) || (V2->type == ST_SCALAR_FLT)) {
+	// bitshift is not valid with float valuess
+	OUT[0].type = ST_SCALAR_FLT;
+	OUT[0].FltValue = NAN;
+	break;
+      }
+      opihi_int M1 = V1[0].IntValue;
+      opihi_int M2 = V2[0].IntValue;
+      OUT[0].type = ST_SCALAR_INT;
+      OUT[0].IntValue = M1 >> M2;
+      break;
+    }
+
     default:
-      snprintf (line, 512, "error: op %c not defined!", op[0]);
+      snprintf (line, 512, "error: op %c not defined for (scalar OP scalar)!", op[0]);
       push_error (line);
       return (FALSE);
@@ -1168,6 +1319,4 @@
 int S_unary (StackVar *OUT, StackVar *V1, char *op) {
 
-  char line[512]; // this is only used to report an error 
-  
 # define S_FUNC(OP,FTYPE) {						\
     if (V1->type == ST_SCALAR_FLT) {					\
@@ -1210,5 +1359,4 @@
   if (!strcmp (op, "asinh"))  S_FUNC(asinh (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "acosh"))  S_FUNC(acosh (M1), ST_SCALAR_FLT);
-  if (!strcmp (op, "lgamma")) S_FUNC(lgamma (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "sin"))    S_FUNC(sin (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "cos"))    S_FUNC(cos (M1), ST_SCALAR_FLT);
@@ -1223,4 +1371,5 @@
   if (!strcmp (op, "dacos"))  S_FUNC(acos (M1)*DEG_RAD, ST_SCALAR_FLT);
   if (!strcmp (op, "datan"))  S_FUNC(atan (M1)*DEG_RAD, ST_SCALAR_FLT);
+  if (!strcmp (op, "lgamma")) S_FUNC(lgamma (M1), ST_SCALAR_FLT);
   if (!strcmp (op, "rnd"))    S_FUNC(0.0*M1 + drand48(), ST_SCALAR_FLT);
   if (!strcmp (op, "drnd"))   S_FUNC(0.0*M1 + drand48(), ST_SCALAR_FLT);
@@ -1235,8 +1384,9 @@
 
   clear_stack (V1);
-  snprintf (line, 512, "error: op %s not defined!", op);
+
+  char line[512]; // this is only used to report an error 
+  snprintf (line, 512, "error: op %s not defined as unary scalar op!", op);
   push_error (line);
   return (FALSE);
-
 }
 
@@ -1295,5 +1445,4 @@
   if (!strcmp (op, "asinh"))  V_FUNC(asinh(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "acosh"))  V_FUNC(acosh(*M1), ST_SCALAR_FLT);
-  if (!strcmp (op, "lgamma")) V_FUNC(lgamma(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "sin"))    V_FUNC(sin(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "cos"))    V_FUNC(cos(*M1), ST_SCALAR_FLT);
@@ -1308,23 +1457,23 @@
   if (!strcmp (op, "dacos"))  V_FUNC(acos(*M1)*DEG_RAD, ST_SCALAR_FLT);
   if (!strcmp (op, "datan"))  V_FUNC(atan(*M1)*DEG_RAD, ST_SCALAR_FLT);
+  if (!strcmp (op, "lgamma")) V_FUNC(lgamma(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "rnd"))    V_FUNC(drand48(), ST_SCALAR_FLT);
   if (!strcmp (op, "drnd"))   V_FUNC(drand48(), ST_SCALAR_FLT);
   if (!strcmp (op, "lrnd"))   V_FUNC(lrand48(), ST_SCALAR_INT);
   if (!strcmp (op, "mrnd"))   V_FUNC(mrand48(), ST_SCALAR_INT);
-  if (!strcmp (op, "ramp"))   V_FUNC(i, ST_SCALAR_INT);
-  if (!strcmp (op, "zero"))   V_FUNC(0, ST_SCALAR_INT);
   if (!strcmp (op, "not"))    V_FUNC(!(*M1), ST_SCALAR_INT);
   if (!strcmp (op, "--"))     V_FUNC(-1*(*M1), ST_SCALAR_INT); // NOTE: opihi_int is signed
   if (!strcmp (op, "isinf"))  V_FUNC(!finite(*M1), ST_SCALAR_FLT);
   if (!strcmp (op, "isnan"))  V_FUNC(isnan((opihi_flt)(*M1)), ST_SCALAR_FLT);
+  if (!strcmp (op, "ramp"))   V_FUNC(i, ST_SCALAR_INT);
   if (!strcmp (op, "xramp"))  V_FUNC(i, ST_SCALAR_INT);
   if (!strcmp (op, "yramp"))  V_FUNC(0, ST_SCALAR_INT);
   if (!strcmp (op, "zramp"))  V_FUNC(0, ST_SCALAR_INT);
+  if (!strcmp (op, "zero"))   V_FUNC(0, ST_SCALAR_INT);
   /* xramp, yramp, zramp above only make sense for matrices. for vectors, xramp = ramp, yramp = zero */
 
 # undef V_FUNC
 
-escape:
-
+  // free the temp vector if needed
   if (V1[0].type == ST_VECTOR_TMP) {
     free (V1[0].vector[0].elements.Ptr);
@@ -1334,9 +1483,24 @@
 
   clear_stack (V1);
+
+  char line[512]; // this is only used to report an error 
+  snprintf (line, 512, "error: op %s not defined as unary vector op!", op);
+  push_error (line);
+  return (FALSE);
+
+escape:
+
+  if (V1[0].type == ST_VECTOR_TMP) {
+    free (V1[0].vector[0].elements.Ptr);
+    free (V1[0].vector);
+    V1[0].vector = NULL;
+  }  
+
+  clear_stack (V1);
   return (TRUE);
 
 }
 
-# define M_FUNC(OP) { for (i = 0; i < Npix; i++, out++, M1++) { *out = (OP); }}
+# define M_FUNC(OP) { for (i = 0; i < Npix; i++, out++, M1++) { *out = (OP); } goto escape; }
 
 int M_unary (StackVar *OUT, StackVar *V1, char *op) {
@@ -1360,5 +1524,5 @@
 // if (!strcmp (op, "rint"))  { for (i = 0; i < Npix; i++, out++, M1++) { *out = nearbyint (*M1); }}
   
-  if (!strcmp (op, "="))     { }
+  if (!strcmp (op, "="))      { goto escape; }
   if (!strcmp (op, "abs"))    M_FUNC(fabs(*M1));
   if (!strcmp (op, "int"))    M_FUNC((opihi_flt)(long long)(*M1));
@@ -1375,5 +1539,4 @@
   if (!strcmp (op, "asinh"))  M_FUNC(asinh(*M1));
   if (!strcmp (op, "acosh"))  M_FUNC(acosh(*M1));
-  if (!strcmp (op, "lgamma")) M_FUNC(lgamma(*M1));
   if (!strcmp (op, "sin"))    M_FUNC(sin(*M1));
   if (!strcmp (op, "cos"))    M_FUNC(cos(*M1));
@@ -1388,14 +1551,15 @@
   if (!strcmp (op, "dacos"))  M_FUNC(acos(*M1)*DEG_RAD);
   if (!strcmp (op, "datan"))  M_FUNC(atan(*M1)*DEG_RAD);
-  if (!strcmp (op, "not"))    M_FUNC(!(*M1));
-  if (!strcmp (op, "--"))     M_FUNC(-(*M1));
+  if (!strcmp (op, "lgamma")) M_FUNC(lgamma(*M1));
   if (!strcmp (op, "rnd"))    M_FUNC(drand48());
   if (!strcmp (op, "drnd"))   M_FUNC(drand48());
   if (!strcmp (op, "lrnd"))   M_FUNC(lrand48());
   if (!strcmp (op, "mrnd"))   M_FUNC(mrand48());
+  if (!strcmp (op, "not"))    M_FUNC(!(*M1));
+  if (!strcmp (op, "--"))     M_FUNC(-(*M1));
   if (!strcmp (op, "ramp"))   M_FUNC(i);
-  if (!strcmp (op, "zero"))   M_FUNC(0);
   if (!strcmp (op, "isinf"))  M_FUNC(!finite(*M1));
   if (!strcmp (op, "isnan"))  M_FUNC(isnan(*M1));
+  if (!strcmp (op, "zero"))   M_FUNC(0);
 
   /* xrm and yrm only make sense for 2D matrices. see special meaning for vectors */
@@ -1411,4 +1575,5 @@
       }
     }
+    goto escape; 
   }
   if (!strcmp (op, "yramp")) {
@@ -1423,4 +1588,5 @@
       }
     }
+    goto escape; 
   }
   if (!strcmp (op, "zramp")) {
@@ -1435,6 +1601,22 @@
       }
     }
-  }
-  
+    goto escape; 
+  }
+  
+  if (V1[0].type == ST_MATRIX_TMP) {
+    free (V1[0].buffer[0].header.buffer);
+    free (V1[0].buffer[0].matrix.buffer);
+    free (V1[0].buffer);
+  }
+
+  clear_stack (V1);
+
+  char line[512]; // this is only used to report an error 
+  snprintf (line, 512, "error: op %s not defined as unary matrix op!", op);
+  push_error (line);
+  return (FALSE);
+
+ escape:
+
   if (V1[0].type == ST_MATRIX_TMP) {
     free (V1[0].buffer[0].header.buffer);