Index: trunk/Ohana/src/kapa2/doc/communication.txt
===================================================================
--- trunk/Ohana/src/kapa2/doc/communication.txt	(revision 26891)
+++ trunk/Ohana/src/kapa2/doc/communication.txt	(revision 26891)
@@ -0,0 +1,119 @@
+
+These are the allowed commands from the pipe, and their required additional messages:
+
+NCUR : end cursor interaction (only valid if already in active cursor mode)
+     NONE
+QUIT : quit kapa
+     NONE
+CURS : turn on active cursor
+     NONE
+PSIT : create PS file
+  KiiSendMessage (fd, "%s", filename);
+PNGF : create PNG file
+  KiiSendMessage (fd, "%s", filename);
+PPMF : create PPM file
+  KiiSendMessage (fd, "%s", filename);
+DBOX : 
+  KiiSendMessage (fd, "%12.6g %12.6g %12.6g %12.6g", graphdata[0].xmin, graphdata[0].xmax, graphdata[0].ymin, graphdata[0].ymax);
+  KiiSendMessage (fd, "%lf %d", graphdata[0].lweight, graphdata[0].color);
+  KiiSendMessage (fd, "%s %s %s", graphdata[0].axis, graphdata[0].labels, graphdata[0].ticks);
+PLOT : 
+  KiiSendMessage (fd, "%8d %8d %d %d %d %d %d %f %f", Npts, data[0].style, data[0].ptype, data[0].ltype, data[0].etype, data[0].ebar, data[0].color, data[0].lweight, data[0].size);
+  KiiSendMessage (fd, "%g %g %g %g", data[0].xmin, data[0].xmax, data[0].ymin, data[0].ymax);
+LABL : 
+  KiiSendMessage (fd, "%6d", mode);
+  KiiSendData (fd, string, strlen(string));
+PTXT : 
+  KiiSendMessage (fd, "%f %f %f", x, y, angle);
+  KiiSendData (fd, string, strlen(string));
+RSIZ : 
+  KiiSendMessage (fd, "%d %d", Nx, Ny); 
+CHAN : 
+  KiiSendMessage (fd, "%1d", channel);
+CMAP : 
+  KiiSendMessage (fd, "%s", colormap);
+CSVE : 
+  KiiSendMessage (fd, "FILE: %d %s", n, file);
+DSEC : 
+  KiiSendMessage (fd, "%s %6.3f %6.3f %6.3f %6.3f", section[0].name, section[0].x, section[0].y, section[0].dx, section[0].dy);
+ERSC : 
+     NONE
+ERSI : 
+     NONE
+ERSO : 
+  KiiSendCommand (fd, 16, "OVERLAY %7d ", n);
+ERSP : 
+     NONE
+ERSS : 
+     NONE
+FONT : 
+  KiiSendCommand (fd, 16, "%s", name);
+  KiiSendCommand (fd, 16, "%d", size);
+GIMC : 
+  KiiScanMessage (fd, "%f %f %f %f", &coords[0].pc1_1, &coords[0].pc2_2, &coords[0].pc1_2, &coords[0].pc2_1); 
+  KiiScanMessage (fd, "%s", coords[0].ctype);
+  KiiScanMessage (fd, "%lf %lf %f %f %f %f", &coords[0].crval1, &coords[0].crval2, &coords[0].crpix1, &coords[0].crpix2, &coords[0].cdelt1, &coords[0].cdelt2);
+GIMD : 
+  KiiScanMessage (fd, "%lf %lf %s %s", &data[0].zero, &data[0].range, data[0].name, data[0].file);
+GIMR : 
+  KiiScanMessage (fd, "%lf %lf %lf %lf", Xmin, Xmax, Ymin, Ymax);
+GLIM : 
+  KiiScanMessage (fd, "%f %f", dx, dy); 
+GSTY : 
+  /* send kapa the plot details */
+  KiiScanMessage (fd, "%d %d %d %d %d %d %lf %lf", &data[0].style, &data[0].ptype, &data[0].ltype, &data[0].etype, &data[0].ebar, &data[0].color, &data[0].lweight, &data[0].size);
+  KiiScanMessage (fd, "%lf %lf %lf %lf", &data[0].xmin, &data[0].xmax, &data[0].ymin, &data[0].ymax);
+  KiiScanMessage (fd, "%f %f %f %f", &data[0].coords.pc1_1, &data[0].coords.pc2_2, &data[0].coords.pc1_2, &data[0].coords.pc2_1);
+  KiiScanMessage (fd, "%d %d %s", &data[0].flipeast, &data[0].flipnorth, data[0].coords.ctype);
+  KiiScanMessage (fd, "%lf %lf %f %f %f %f", &data[0].coords.crval1, &data[0].coords.crval2, &data[0].coords.crpix1, &data[0].coords.crpix2, &data[0].coords.cdelt1, &data[0].coords.cdelt2);
+LOAD : (eac write follewed by 'DONE')
+  KiiSendMessage (fd, "%d %d %d %d", overnum, Noverlay, Ntext, Ntextdata);
+  write (fd, buffer, Noverlay*sizeof(KiiOverlayBase)); 
+  write (fd, textdata, Ntextdata);
+LSEC : 
+  KiiSendMessage (fd, "%s", name);
+MOVE : 
+  KiiSendMessage (fd, "%d %d", x, y); 
+MSEC : 
+  KiiSendMessage (fd, "%s %s", name, direction);
+READ : 
+  KiiScanMessage (fd, "%d", &Ncolors);
+  KiiSendMessage (fd, "%8d %8d", image[0].Nx, image[0].Ny);
+  KiiSendMessage (fd, "%f %f %s %s", data[0].zero, data[0].range, data[0].name, data[0].file);
+  KiiSendMessage (fd, "%f %f %d ", min, max, size);
+  KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval1, coords[0].crpix1, coords[0].cdelt1, coords[0].pc1_1, coords[0].pc1_2);
+  KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval2, coords[0].crpix2, coords[0].cdelt2, coords[0].pc2_1, coords[0].pc2_2);
+  KiiSendMessage (fd, "%s", coords[0].ctype);
+  write (fd, image[0].data1d, size);
+  (note that 1D, 2D, and 1D_8bit have the same messaging)
+SAVE : 
+  KiiSendMessage (fd, "FILE: %d %s", n, file);
+SIMC : 
+  KiiSendMessage (fd, "%g %g %g %g", coords[0].pc1_1, coords[0].pc2_2, coords[0].pc1_2, coords[0].pc2_1);
+  KiiSendMessage (fd, "%s", coords[0].ctype);
+  KiiSendMessage (fd, "%g %g %g %g %g %g", coords[0].crval1, coords[0].crval2, coords[0].crpix1, coords[0].crpix2, coords[0].cdelt1, coords[0].cdelt2);
+SIMD : 
+  KiiSendMessage (fd, "%g %g %s %s", data[0].zero, data[0].range, data[0].name, data[0].file);
+SLIM : 
+  KiiSendMessage (fd, "%g %g %g %g ", data[0].xmin, data[0].xmax, data[0].ymin, data[0].ymax);
+SSEC : 
+  KiiSendMessage (fd, "%s", name);
+SSTY : 
+  KiiSendMessage (fd, "%8d %d %d %d %d %d %f %f", data[0].style, data[0].ptype, data[0].ltype, data[0].etype, data[0].ebar, data[0].color, data[0].lweight, data[0].size);
+  KiiSendMessage (fd, "%g %g %g %g", data[0].xmin, data[0].xmax, data[0].ymin, data[0].ymax);
+  KiiSendMessage (fd, "%g %g %g %g", data[0].coords.pc1_1, data[0].coords.pc2_2, data[0].coords.pc1_2, data[0].coords.pc2_1);
+  KiiSendMessage (fd, "%d %d %s", data[0].flipeast, data[0].flipnorth, data[0].coords.ctype);
+  KiiSendMessage (fd, "%g %g %g %g %g %g", data[0].coords.crval1, data[0].coords.crval2, data[0].coords.crpix1, data[0].coords.crpix2, data[0].coords.cdelt1, data[0].coords.cdelt2);
+TOOL : 
+  KiiSendMessage (fd, "%d", location); 
+JPEG : create JPEG file
+  KiiSendMessage (fd, "%s", filename);
+CENT : 
+  KiiSendMessage (fd, "%8.3f %8.3f %8d ", x, y, zoom);
+PARI : 
+  KiiSendMessage (fd, "%2d %2d ", xflip, yflip);
+PLOB :
+  KiiSendMessage (fd, "%s %d %d %d", type, Npts, Nbytes, swap); 
+  write (fd, values, Nbytes);
+
+NPIX : ?? deprecated?
Index: trunk/Ohana/src/kapa2/include/constants.h
===================================================================
--- trunk/Ohana/src/kapa2/include/constants.h	(revision 26884)
+++ trunk/Ohana/src/kapa2/include/constants.h	(revision 26891)
@@ -11,8 +11,8 @@
 
 # define NCHANNELS 3
-# define NPIXELS_DYNAMIC 4600
+# define NPIXELS_DYNAMIC 128
 
 // XXX for the moment, this is set to match the values in SetColorScale3D_CC
-# define NPIXELS_STATIC 4600
+# define NPIXELS_STATIC 128
 
 # define PAD1  3
Index: trunk/Ohana/src/kapa2/src/ColorHistogram.c
===================================================================
--- trunk/Ohana/src/kapa2/src/ColorHistogram.c	(revision 26884)
+++ trunk/Ohana/src/kapa2/src/ColorHistogram.c	(revision 26891)
@@ -174,5 +174,5 @@
     node->pixel = current[0];
     current[0] ++;
-    assert (current[0] <= Npixels);
+    assert (current[0] < Npixels);
     return TRUE;
   }
Index: trunk/Ohana/src/kapa2/src/EraseOverlay.c
===================================================================
--- trunk/Ohana/src/kapa2/src/EraseOverlay.c	(revision 26884)
+++ trunk/Ohana/src/kapa2/src/EraseOverlay.c	(revision 26891)
@@ -8,10 +8,11 @@
   KapaImageWidget *image;
 
+  // We have to accept the incoming message even if we cannot perform the action
+  KiiScanCommand (sock, 16, "%*s %d", &N);
+
   graphic = GetGraphic();
   section = GetActiveSection();
   image = section->image;
   if (image == NULL) return (TRUE);
-
-  KiiScanCommand (sock, 16, "%*s %d", &N);
 
   if (N > NOVERLAYS) {
Index: trunk/Ohana/src/kapa2/src/LoadPicture.c
===================================================================
--- trunk/Ohana/src/kapa2/src/LoadPicture.c	(revision 26884)
+++ trunk/Ohana/src/kapa2/src/LoadPicture.c	(revision 26891)
@@ -73,4 +73,5 @@
     if (status == 0) {  /* No more pipe */
       fprintf (stderr, "error: pipe closed\n");
+      fcntl (sock, F_SETFL, !O_NONBLOCK);  
       return (FALSE);
     }
@@ -92,6 +93,4 @@
   SetColorScale (graphic, image);
 
-  SetColorScale (graphic, image);
-
   if (!USE_XWINDOW) return (TRUE);
 
Index: trunk/Ohana/src/kapa2/src/PNGit.c
===================================================================
--- trunk/Ohana/src/kapa2/src/PNGit.c	(revision 26884)
+++ trunk/Ohana/src/kapa2/src/PNGit.c	(revision 26891)
@@ -16,9 +16,5 @@
 
   /* expect a line telling the number of bytes and a filename */
-  status = read (sock, filename, 16);
-  filename[16] = 0;
-  sscanf (filename, "%*s %d", &Nbytes);
-  status = read (sock, filename, Nbytes);
-  filename[status] = 0; /* make the string easy to parse */
+  KiiScanMessage (sock, "%s", filename);
 
   f = fopen (filename, "w");
Index: trunk/Ohana/src/kapa2/src/PPMit.c
===================================================================
--- trunk/Ohana/src/kapa2/src/PPMit.c	(revision 26884)
+++ trunk/Ohana/src/kapa2/src/PPMit.c	(revision 26891)
@@ -13,9 +13,5 @@
 
   /* expect a line telling the number of bytes and a filename */
-  status = read (sock, filename, 16);
-  filename[16] = 0;
-  sscanf (filename, "%*s %d", &Nbytes);
-  status = read (sock, filename, Nbytes);
-  filename[status] = 0; /* make the string easy to parse */
+  KiiScanMessage (sock, "%s", filename);
 
   f = fopen (filename, "w");
Index: trunk/Ohana/src/kapa2/src/Resize.c
===================================================================
--- trunk/Ohana/src/kapa2/src/Resize.c	(revision 26884)
+++ trunk/Ohana/src/kapa2/src/Resize.c	(revision 26891)
@@ -11,4 +11,5 @@
   graphic = GetGraphic();
 
+  // must scan the message before possible return
   KiiScanMessage (sock, "%d %d", &NX, &NY);
 
Index: trunk/Ohana/src/kapa2/src/SaveOverlay.c
===================================================================
--- trunk/Ohana/src/kapa2/src/SaveOverlay.c	(revision 26884)
+++ trunk/Ohana/src/kapa2/src/SaveOverlay.c	(revision 26891)
@@ -9,10 +9,10 @@
   KapaImageWidget *image;
 
+  KiiScanMessage (sock, "%*s %d %s", &N, filename);
+  
   section = GetActiveSection();
   image   = section->image;
   if (image == NULL) return (TRUE);
 
-  KiiScanMessage (sock, "%*s %d %s", &N, filename);
-  
   f = fopen (filename, "w");
   if (f == NULL) {
Index: trunk/Ohana/src/kapa2/src/SetColorScale.c
===================================================================
--- trunk/Ohana/src/kapa2/src/SetColorScale.c	(revision 26884)
+++ trunk/Ohana/src/kapa2/src/SetColorScale.c	(revision 26891)
@@ -34,12 +34,13 @@
   float slope;
   float start;
-  unsigned short MaxValue;
+  unsigned short MaxValue, NANValue;
   Matrix *matrix;
 
   // define the color transform parameters
-  MaxValue = graphic[0].Npixels - 1;
+  NANValue = graphic[0].Npixels - 1;
+  MaxValue = graphic[0].Npixels - 2;
   if (image[0].image[0].range != 0.0) {
-    slope = graphic[0].Npixels / image[0].image[0].range;
-    start = graphic[0].Npixels * image[0].image[0].zero / image[0].image[0].range;
+    slope = (graphic[0].Npixels - 1) / image[0].image[0].range;
+    start = slope * image[0].image[0].zero;
   } else {
     slope = 1.0;
@@ -63,7 +64,17 @@
   // convert pixel data values to pixel index values (0 - Npixel)
   for (i = 0; i < nPixels; i++, iData++, oData++) {
+    if (isnan(*iData) || isinf(*iData)) {
+      *oData = NANValue;
+      continue;
+    }
     value = *iData * slope - start;
-    if (value < 0) value = 0;
-    if (value > MaxValue) value = MaxValue;
+    if (value < 0) {
+      *oData = 0;
+      continue;
+    }
+    if (value > MaxValue) {
+      *oData = MaxValue;
+      continue;
+    }
     *oData = value;
   }
@@ -169,4 +180,13 @@
   if (DX != image[0].channel[KAPA_CHANNEL_BLUE].matrix.Naxis[0]) return FALSE;
   if (DY != image[0].channel[KAPA_CHANNEL_BLUE].matrix.Naxis[1]) return FALSE;
+
+  // XXX a bit bogus: we are going to store a special color at Npixels-1 for 
+  // marking the NAN pixels.  
+  { 
+      // XXX set pure green for now
+      graphic[0].cmap[graphic[0].Npixels-1].red = 0;
+      graphic[0].cmap[graphic[0].Npixels-1].green = 0xffff;
+      graphic[0].cmap[graphic[0].Npixels-1].blue = 0;
+  }  
 
   // create the top-level cube
@@ -214,5 +234,5 @@
   // convert histogram 3D values to cmap colors
   Npixels = 0;
-  CCNodeSetColorMap (cube, graphic[0].cmap, graphic[0].Npixels, &Npixels);
+  CCNodeSetColorMap (cube, graphic[0].cmap, graphic[0].Npixels - 1, &Npixels);
 
   // store the colors
Index: trunk/Ohana/src/kapa2/src/SetColormap.c
===================================================================
--- trunk/Ohana/src/kapa2/src/SetColormap.c	(revision 26884)
+++ trunk/Ohana/src/kapa2/src/SetColormap.c	(revision 26891)
@@ -10,4 +10,9 @@
     VAR = tmp; \
   } }   
+
+// the image arrives from the client program as an array of floats, stored in image->matrix.
+// whenever we load a new image or change tv channels, we need to map the currently active
+// image from float to pixel index (image->pixmap).  When we actually display (or change the
+// display of) an image, we simply remap the lit pixels (in Remap32.c, etc).
 
 int SetColormap (char *name) {
@@ -26,9 +31,18 @@
   }
 
+  // XXX a bit bogus: we are going to store a special color at Npixels-1 for 
+  // marking the NAN pixels.  
+  { 
+      // XXX set pure green for now
+      graphic[0].cmap[graphic[0].Npixels-1].red = 0;
+      graphic[0].cmap[graphic[0].Npixels-1].green = 0xffff;
+      graphic[0].cmap[graphic[0].Npixels-1].blue = 0;
+  }  
+
   // very simple color model: evenly spaced cube 
   if (!strcasecmp (name, "ruffcolor")) {
       graphic[0].nRed  = pow (graphic[0].Npixels, 0.333);
       graphic[0].nBlue = pow (graphic[0].Npixels, 0.333);
-      graphic[0].nGreen = graphic[0].Npixels / (graphic[0].nRed * graphic[0].nBlue);
+      graphic[0].nGreen = (graphic[0].Npixels - 1) / (graphic[0].nRed * graphic[0].nBlue);
 
       // red,green,blue are values in range 0x0000 to 0xffff
@@ -48,4 +62,5 @@
 	  }
       }
+      fprintf (stderr, "ruff Npix: %d vs %d\n", i, graphic[0].Npixels);
 
       // all other modes are 1D: set the flag:
@@ -58,9 +73,9 @@
 
   // red,green,blue are values in range 0x0000 to 0xffff
-  scale = 0xffff / (graphic[0].Npixels - 1);
+  scale = 0xffff / (graphic[0].Npixels - 2);
 
   /* greyscale */
   if ((!strcasecmp (name, "grayscale")) || (!strcasecmp (name, "greyscale"))) {
-    for (i = 0; i < graphic[0].Npixels; i++) {  
+    for (i = 0; i < graphic[0].Npixels - 1; i++) {  
       SETVALUE (graphic[0].cmap[i].red,   0xffff - i*scale, 0, 0xffff);
       SETVALUE (graphic[0].cmap[i].green, 0xffff - i*scale, 0, 0xffff);
@@ -72,5 +87,5 @@
   /* -grayscale */
   if ((!strcasecmp (name, "-grayscale")) || (!strcasecmp (name, "-greyscale"))) {
-    for (i = 0; i < graphic[0].Npixels; i++) {  
+    for (i = 0; i < graphic[0].Npixels - 1; i++) {  
       SETVALUE (graphic[0].cmap[i].red,   i*scale, 0, 0xffff);
       SETVALUE (graphic[0].cmap[i].green, i*scale, 0, 0xffff);
@@ -102,5 +117,5 @@
       graphic[0].cmap[i].flags = DoRed | DoGreen | DoBlue;
     }
-    for (i = (int)(0.75*graphic[0].Npixels); i < graphic[0].Npixels; i++) {  
+    for (i = (int)(0.75*graphic[0].Npixels); i < graphic[0].Npixels - 1; i++) {  
       graphic[0].cmap[i].red   = 0xffff;
       graphic[0].cmap[i].green = 0xffff;
@@ -134,5 +149,5 @@
     }
     blueRef  = 0.75*graphic[0].Npixels*scale*4.0;
-    for (i = (int)(0.75*graphic[0].Npixels); i < graphic[0].Npixels; i++) {  
+    for (i = (int)(0.75*graphic[0].Npixels); i < graphic[0].Npixels - 1; i++) {  
       graphic[0].cmap[i].red   = 0xffff;
       graphic[0].cmap[i].green = 0xffff;
Index: trunk/Ohana/src/libkapa/src/KiiOverlay.c
===================================================================
--- trunk/Ohana/src/libkapa/src/KiiOverlay.c	(revision 26884)
+++ trunk/Ohana/src/libkapa/src/KiiOverlay.c	(revision 26891)
@@ -92,10 +92,10 @@
 
   write (fd, textdata, Ntextdata);
-  KiiWaitAnswer (fd, "DONE");
+  KiiWaitAnswer (fd, "DONE"); // this 'DONE' notes the end of the textdata buffer
 
   free (buffer);
   free (textdata);
 
-  KiiWaitAnswer (fd, "DONE");
+  KiiWaitAnswer (fd, "DONE"); // this 'DONE' notes the end of the command
   return (TRUE);
 }
Index: trunk/Ohana/src/opihi/cmd.astro/mkgauss.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/mkgauss.c	(revision 26884)
+++ trunk/Ohana/src/opihi/cmd.astro/mkgauss.c	(revision 26891)
@@ -11,5 +11,6 @@
   Buffer *buf;
 
-  Xo = Yo = 0;
+  // this should be Nx/2, Ny/2 if not set
+  Xo = Yo = NAN;
   if ((N = get_argument (argc, argv, "-c"))) {
     remove_argument (N, &argc, argv);
@@ -30,4 +31,7 @@
   Ny = buf[0].header.Naxis[1];
   
+  if (isnan(Xo)) Xo = Nx / 2;
+  if (isnan(Yo)) Yo = Ny / 2;
+
   /* gaussian parameters */
   Sig_x = atof (argv[2]);
Index: trunk/Ohana/src/opihi/cmd.astro/profile.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/profile.c	(revision 26884)
+++ trunk/Ohana/src/opihi/cmd.astro/profile.c	(revision 26891)
@@ -31,12 +31,12 @@
   ResetVector (yvec, OPIHI_FLT, (int)SQ(2*N+1));
 
-  bzero (yvec[0].elements.Flt, (int)SQ(2*N+1)*sizeof(opihi_flt)+1);
+  bzero (yvec[0].elements.Flt, (int)SQ(2*N+1)*sizeof(opihi_flt));
   V = (float *)(buf[0].matrix.buffer); 
   Npt = 0;
   Nx = buf[0].matrix.Naxis[0];
-  for (i = sx - N; i < sx + N; i++) {
-    for (j = sy - N; j < sy + N; j++, Npt++) {
+  for (i = sx - N; i <= sx + N; i++) {
+    for (j = sy - N; j <= sy + N; j++, Npt++) {
       yvec[0].elements.Flt[Npt] = V[i + j*Nx];
-      xvec[0].elements.Flt[Npt] = hypot (i - sx, j - sy);
+      xvec[0].elements.Flt[Npt] = hypot (i + 0.5 - sx, j + 0.5 - sy);
     }
   }
Index: trunk/Ohana/src/opihi/cmd.astro/star.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/star.c	(revision 26884)
+++ trunk/Ohana/src/opihi/cmd.astro/star.c	(revision 26891)
@@ -24,4 +24,5 @@
   if ((argc != 4) && (argc != 5)) {
     gprint (GP_ERR, "USAGE: star (buffer) x y [dx] [-border N] [-sat cnts]\n");
+    gprint (GP_ERR, " dx is the aperture diameter, but is adjusted up to the next odd number\n");
     return (FALSE);
   }
Index: trunk/Ohana/src/opihi/cmd.basic/Makefile
===================================================================
--- trunk/Ohana/src/opihi/cmd.basic/Makefile	(revision 26884)
+++ trunk/Ohana/src/opihi/cmd.basic/Makefile	(revision 26891)
@@ -39,4 +39,5 @@
 $(SRC)/mkdir.$(ARCH).o     \
 $(SRC)/module.$(ARCH).o     \
+$(SRC)/nop.$(ARCH).o        \
 $(SRC)/output.$(ARCH).o     \
 $(SRC)/quit.$(ARCH).o	     \
Index: trunk/Ohana/src/opihi/cmd.basic/init.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.basic/init.c	(revision 26884)
+++ trunk/Ohana/src/opihi/cmd.basic/init.c	(revision 26891)
@@ -25,4 +25,5 @@
 int mkdir_opihi     PROTO((int, char **));
 int module          PROTO((int, char **));
+int nop             PROTO((int, char **));
 int output          PROTO((int, char **));
 int pwd             PROTO((int, char **));
@@ -66,4 +67,5 @@
   {1, "?",             list_help,          "list commands *"},
   {1, "??",            list_vars,          "list variables *"},
+  {1, "#",             nop,                "a NOP function"},
   {1, "local",         local,              "define local variables"},
   {1, "macro",         macro,              "deal with the macros *"}, 
@@ -71,4 +73,5 @@
   {1, "mkdir",         mkdir_opihi,        "built-in mkdir command"},
   {1, "module",        module,             "load script file from the modules directories"},
+  {1, "nop",           nop,                "a NOP function"},
   {1, "output",        output,             "redirect output to file"},
   {1, "pwd",           pwd,                "print current working directory"},
Index: trunk/Ohana/src/opihi/cmd.basic/nop.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.basic/nop.c	(revision 26891)
+++ trunk/Ohana/src/opihi/cmd.basic/nop.c	(revision 26891)
@@ -0,0 +1,5 @@
+# include "basic.h"
+
+int nop (int argc, char **argv) {
+  return (TRUE);
+}
Index: trunk/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 26884)
+++ trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 26891)
@@ -53,4 +53,6 @@
 $(SRC)/ungridify.$(ARCH).o     \
 $(SRC)/histogram.$(ARCH).o	\
+$(SRC)/hermitian1d.$(ARCH).o	\
+$(SRC)/hermitian2d.$(ARCH).o	\
 $(SRC)/imcut.$(ARCH).o	 	\
 $(SRC)/imhist.$(ARCH).o	\
@@ -69,4 +71,5 @@
 $(SRC)/load.$(ARCH).o		\
 $(SRC)/lookup.$(ARCH).o	\
+$(SRC)/matrix.$(ARCH).o	\
 $(SRC)/mkrgb.$(ARCH).o	\
 $(SRC)/mcreate.$(ARCH).o	\
Index: trunk/Ohana/src/opihi/cmd.data/hermitian1d.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/hermitian1d.c	(revision 26891)
+++ trunk/Ohana/src/opihi/cmd.data/hermitian1d.c	(revision 26891)
@@ -0,0 +1,61 @@
+# include "data.h"
+
+int hermitian1d (int argc, char **argv) {
+  
+  int i, N, order, Nvec, Poly;
+  opihi_int *inI;
+  opihi_flt *out, *inF;
+  opihi_flt mean, sigma, x;
+  double nf, norm;
+  Vector *xvec, *yvec;
+
+  Poly = FALSE;
+  if ((N = get_argument (argc, argv, "-poly"))) {
+    Poly = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 6) goto usage;
+
+  /* select input / output buffers */
+  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  Nvec = xvec[0].Nelements;
+  ResetVector (yvec, OPIHI_FLT, Nvec);
+
+  /* gaussian parameters */
+  mean  = atof (argv[3]);
+  sigma = atof (argv[4]);
+  order = atoi (argv[5]);
+  if (order < 0) goto usage;
+  if (order > 10) goto usage;
+
+  out = yvec[0].elements.Flt;
+  inF = xvec[0].elements.Flt;
+  inI = xvec[0].elements.Int;
+
+  nf = exp(lgamma(order + 1));
+  norm = 1.0 / sqrt(nf*sqrt(2*M_PI)*sigma);
+
+  // a little sub-optimal : split this up with macros?
+  for (i = 0; i < Nvec; i++, out++) {
+    if (xvec[0].type == OPIHI_FLT) {
+      x = (inF[i] - mean) / sigma;
+    } else {
+      x = (inI[i] - mean) / sigma;
+    }
+    *out = hermitian_polynomial (x, order);
+    if (!Poly) {
+      *out *= norm*exp (-0.25*x*x);
+    }
+  }
+  return (TRUE);
+
+usage:
+  gprint (GP_ERR, "USAGE: hermitian1d (x) (y) (mean) (sigma) (order)\n");
+  gprint (GP_ERR, "  Note : only orders up to 10 are supported\n");
+  return (FALSE);
+
+}
+
+// {\psi}_n(x) = \frac{1}{\sqrt{n! \, 2^n\sqrt{\pi}}}\, \mathrm{e}^{-x^2/2}H_n(x).\,\! 
Index: trunk/Ohana/src/opihi/cmd.data/hermitian2d.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/hermitian2d.c	(revision 26891)
+++ trunk/Ohana/src/opihi/cmd.data/hermitian2d.c	(revision 26891)
@@ -0,0 +1,79 @@
+# include "astro.h"
+
+int hermitian2d (int argc, char **argv) {
+  
+  int i, j, N, Nx, Ny, Xorder, Yorder, Poly;
+  float *in;
+  double Xo, Yo, sigma;
+  double x, y, xf, yf, value, norm, weight;
+  double nf, mf;
+  Buffer *buf;
+
+  // optionally return the Hermitian Polynomial or the Hermitian Function
+  Poly = FALSE; 
+  if ((N = get_argument (argc, argv, "-poly"))) {
+    Poly = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  Xo = Yo = NAN;
+  if ((N = get_argument (argc, argv, "-c"))) {
+    remove_argument (N, &argc, argv);
+    Xo = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    Yo = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 5) goto usage;
+
+  /* select input / output buffers */
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  Nx = buf[0].header.Naxis[0];
+  Ny = buf[0].header.Naxis[1];
+  
+  if (isnan(Xo)) Xo = Nx / 2;
+  if (isnan(Yo)) Yo = Ny / 2;
+
+  /* gaussian parameters */
+  sigma = atof (argv[2]);
+  Xorder = atof (argv[3]);
+  Yorder = atof (argv[4]);
+  if (Xorder < 0) goto usage;
+  if (Yorder < 0) goto usage;
+  if (Xorder > 10) goto usage;
+  if (Yorder > 10) goto usage;
+
+  nf = exp(lgamma(Xorder + 1));
+  mf = exp(lgamma(Yorder + 1));
+  norm = 1.0 / sqrt(nf*mf*2*M_PI) / sigma;
+
+  in = (float *) buf[0].matrix.buffer;
+  for (j = 0; j < Ny; j++) {
+    for (i = 0; i < Nx; i++, in++) {
+
+      x = (i - Xo) / sigma;
+      y = (j - Yo) / sigma;
+
+      // not sure what a 2D hermitian looks like
+      // is it H_i(x) * H_j(y)?
+
+      xf = hermitian_polynomial (x, Xorder);
+      yf = hermitian_polynomial (y, Yorder);
+      value = xf*yf;
+      if (!Poly) {
+	weight = norm*exp(-0.25*(x*x + y*y));
+	value *= weight;
+      }
+
+      *in += value;
+    }
+  }
+
+  return (TRUE);
+
+usage:
+  gprint (GP_ERR, "USAGE: hermitian2d (buffer) (sigma) Xorder Yorder\n");
+  gprint (GP_ERR, "  Note : only orders up to 10 are supported\n");
+  return (FALSE);
+}
Index: trunk/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/init.c	(revision 26884)
+++ trunk/Ohana/src/opihi/cmd.data/init.c	(revision 26891)
@@ -42,4 +42,6 @@
 int ungridify        PROTO((int, char **));
 int histogram        PROTO((int, char **));
+int hermitian1d      PROTO((int, char **));
+int hermitian2d      PROTO((int, char **));
 int imcut            PROTO((int, char **));
 int imhist           PROTO((int, char **));
@@ -58,4 +60,5 @@
 int load             PROTO((int, char **));
 int lookup           PROTO((int, char **));
+int matrix           PROTO((int, char **));
 int mkrgb            PROTO((int, char **));
 int mcreate          PROTO((int, char **));
@@ -172,4 +175,6 @@
   {1, "header",       header,           "print image header"},
   {1, "histogram",    histogram,        "generate histogram from vector"},
+  {1, "hermitian1d",  hermitian1d,      "generate 1-D Hermitian Polynomial"},
+  {1, "hermitian2d",  hermitian2d,      "generate 2-D Hermitian Polynomial"},
   {1, "imbin",        rebin,            "rebin image data by factor of N"},
   {1, "imclip",       imclip,           "clip values in an image to be within a range"},
@@ -195,4 +200,5 @@
   {1, "minterp",      minterp,          "interpolate image pixels"},
   {1, "iminterp",     minterp,          "interpolate image pixels"},
+  {1, "matrix",       matrix,           "matrix math operations"},
   {1, "mkrgb",        mkrgb,            "convert 3 images to rgb jpeg (use Kapa for better control)"},
   {1, "mset",         mset,             "insert a vector in an image"},
Index: trunk/Ohana/src/opihi/cmd.data/matrix.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/matrix.c	(revision 26891)
+++ trunk/Ohana/src/opihi/cmd.data/matrix.c	(revision 26891)
@@ -0,0 +1,106 @@
+# include "data.h"
+
+int matrix (int argc, char **argv) {
+  
+  int ix, iy, ik, Nx, Ny, Nk;
+  float *Va, *Vb, *Vc, sum;
+  Buffer *bufA, *bufB, *bufC;
+
+  if ((argc != 5) && (argc != 7)) {
+    gprint (GP_ERR, "USAGE: matrix set A = B * C : matrix multiplication\n");
+    gprint (GP_ERR, "USAGE: matrix transpose A to B\n");
+    return (FALSE);
+  }
+
+  if (!strcasecmp (argv[1], "set")) {
+
+    if ((argc != 7) || strcasecmp (argv[3], "=") || strcasecmp (argv[5], "*")) {
+      gprint (GP_ERR, "USAGE: matrix set A = B * C : matrix multiplication\n");
+      return (FALSE);
+    }
+    if ((bufB = SelectBuffer (argv[4], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+    if ((bufC = SelectBuffer (argv[6], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+
+    if (bufB[0].matrix.Naxis[0] != bufC[0].matrix.Naxis[1]) {
+      gprint (GP_ERR, "size mis-match in matrices: (%d x %d) * (%d x %d)\n", 
+	      bufB[0].matrix.Naxis[0], bufB[0].matrix.Naxis[1], 
+	      bufC[0].matrix.Naxis[0], bufC[0].matrix.Naxis[1]);
+      return (FALSE);
+    }
+    if ((bufA = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+
+    Nx = bufC[0].matrix.Naxis[0];
+    Ny = bufB[0].matrix.Naxis[1];
+    Nk = bufB[0].matrix.Naxis[0];
+
+    gfits_free_matrix (&bufA[0].matrix);
+    gfits_free_header (&bufA[0].header);
+
+    bufA[0].bitpix = bufB[0].bitpix;
+    bufA[0].unsign = bufB[0].unsign;
+    bufA[0].bscale = bufB[0].bscale;
+    bufA[0].bzero  = bufB[0].bzero;
+    gfits_copy_header (&bufB[0].header, &bufA[0].header);
+    gfits_modify (&bufA[0].header, "NAXIS1", "%d", 1, Nx);
+    gfits_modify (&bufA[0].header, "NAXIS2", "%d", 1, Ny);
+    bufA[0].header.Naxis[0] = Nx;
+    bufA[0].header.Naxis[1] = Ny;
+    gfits_create_matrix (&bufA[0].header, &bufA[0].matrix);
+
+    Va = (float *)bufA[0].matrix.buffer;
+    Vb = (float *)bufB[0].matrix.buffer;
+    Vc = (float *)bufC[0].matrix.buffer;
+
+    for (iy = 0; iy < Ny; iy++) {
+      for (ix = 0; ix < Nx; ix++) {
+	sum = 0.0;
+	for (ik = 0; ik < Nk; ik++) {
+	  sum += Vb[iy*Nk + ik] * Vc[ik*Nx + ix];
+	}
+	Va[iy*Nx + ix] = sum;
+      }
+    }
+    return (TRUE);
+  }
+
+  if (!strcasecmp (argv[1], "transpose")) {
+    if ((argc != 5) || strcasecmp (argv[3], "to")) {
+      gprint (GP_ERR, "USAGE: matrix transpose A to B\n");
+      return (FALSE);
+    }
+    if ((bufA = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+    if ((bufB = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+
+    Nx = bufA[0].matrix.Naxis[0];
+    Ny = bufA[0].matrix.Naxis[1];
+
+    gfits_free_matrix (&bufB[0].matrix);
+    gfits_free_header (&bufB[0].header);
+
+    bufB[0].bitpix = bufA[0].bitpix;
+    bufB[0].unsign = bufA[0].unsign;
+    bufB[0].bscale = bufA[0].bscale;
+    bufB[0].bzero  = bufA[0].bzero;
+    gfits_copy_header (&bufA[0].header, &bufB[0].header);
+    gfits_modify (&bufB[0].header, "NAXIS1", "%d", 1, Ny);
+    gfits_modify (&bufB[0].header, "NAXIS2", "%d", 1, Nx);
+    bufB[0].header.Naxis[0] = Ny;
+    bufB[0].header.Naxis[1] = Nx;
+    gfits_create_matrix (&bufB[0].header, &bufB[0].matrix);
+
+    Va = (float *)bufA[0].matrix.buffer;
+    Vb = (float *)bufB[0].matrix.buffer;
+
+    for (iy = 0; iy < Ny; iy++) {
+      for (ix = 0; ix < Nx; ix++) {
+	Vb[ix*Ny + iy] = Va[iy*Nx + ix];
+      }
+    }
+
+    return (TRUE);
+  }
+
+    gprint (GP_ERR, "USAGE: matrix set A = B * C : matrix and/or vector multiplication\n");
+    gprint (GP_ERR, "USAGE: matrix transpose A to B\n");
+    return (FALSE);
+}
Index: trunk/Ohana/src/opihi/cmd.data/svd.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/svd.c	(revision 26884)
+++ trunk/Ohana/src/opihi/cmd.data/svd.c	(revision 26891)
@@ -57,4 +57,56 @@
   /* use a bcopy instead? */
 
+  // try C.R. Bond's version -- requires matrices in the form A[row][col] not A[row*Ncol + col]
+  if (1) { 
+      int j;
+      double **a, **u, **v, *q;
+      ALLOCATE (a, double *, Ny);
+      ALLOCATE (u, double *, Ny);
+      ALLOCATE (v, double *, Ny);
+      for (i = 0; i < Ny; i++) {
+	  ALLOCATE (a[i], double, Nx);
+	  ALLOCATE (u[i], double, Nx);
+	  ALLOCATE (v[i], double, Nx);
+      }	  
+      ALLOCATE (q, double, Nx);
+
+      for (j = 0; j < Ny; j++) {
+	  for (i = 0; i < Nx; i++) {
+	      a[j][i] = A[j*Nx + i];
+	      u[j][i] = 0;
+	      v[j][i] = 0;
+	  }
+      }
+      for (i = 0; i < Nx; i++) {
+	  q[i] = 0;
+      }
+
+      status = svdcmp_bond_raw (Ny, Nx, 1, 1, FLT_EPSILON, 1e-6, a, q, u, v);
+      fprintf (stderr, "status: %d\n", status);
+
+      // copy u q v back to U W V:
+      for (j = 0; j < Ny; j++) {
+	  for (i = 0; i < Nx; i++) {
+	      U[j*Nx + i] = u[j][i];
+	      V[j*Nx + i] = v[j][i];
+	  }
+      }
+      for (i = 0; i < Nx; i++) {
+	  W[i] = q[i];
+      }
+
+      for (j = 0; j < Ny; j++) {
+	  free(a[j]);
+	  free(u[j]);
+	  free(v[j]);
+      }
+      free (a);
+      free (u);
+      free (v);
+      free (q);
+
+      return TRUE;
+  }
+
   status = svdcmp (U, W, V, Nx, Ny);
   if (!status) {
Index: trunk/Ohana/src/opihi/cmd.data/test/matrix.sh
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/test/matrix.sh	(revision 26891)
+++ trunk/Ohana/src/opihi/cmd.data/test/matrix.sh	(revision 26891)
@@ -0,0 +1,78 @@
+
+list tests
+ test1
+ test2
+end
+
+# a very simple diagonal matrix equation
+macro test1
+ $PASS = 1
+ break -auto off
+
+ delete -q A B C
+ mcreate B 2 2
+ mcreate C 2 2
+
+ B[0][0] = 2
+ B[1][0] = -1
+ B[0][1] = -1
+ B[1][1] = 4
+
+ C[0][0] = 2
+ C[1][0] = 1
+ C[0][1] = 1
+ C[1][1] = 1
+
+ matrix set A = B * C
+
+ # A[0][0] = B[0][0]*C[0][0] + B[1][0]*C[0][1] =  4 - 1 = 3
+ # A[1][0] = B[0][0]*C[1][0] + B[1][0]*C[1][1] =  2 - 1 = 1
+ # A[0][1] = B[0][1]*C[0][0] + B[1][1]*C[0][1] = -2 + 4 = 2
+ # A[1][1] = B[0][1]*C[1][0] + B[1][1]*C[1][1] = -1 + 4 = 3
+
+ if (abs(A[0][0] - 3.0) > 0.0001)
+  $PASS = 0
+ end
+ if (abs(A[0][1] - 2.0) > 0.0001)
+  $PASS = 0
+ end
+ if (abs(A[1][0] - 1.0) > 0.0001)
+  $PASS = 0
+ end
+ if (abs(A[1][1] - 3.0) > 0.0001)
+  $PASS = 0
+ end
+end
+
+# a very simple off-diagonal matrix equation
+macro test2
+ $PASS = 1
+ break -auto off
+
+ delete -q A B
+
+ mcreate A 2 2
+
+ A[0][0] =  2
+ A[1][0] = +1
+ A[0][1] = -1
+ A[1][1] =  4
+
+ matrix transpose A to B
+ if (not($STATUS))
+   $PASS = 0
+ end
+
+ if (abs(B[0][0] - 2.0) > 0.0001)
+  $PASS = 0
+ end
+ if (abs(B[0][1] - 1.0) > 0.0001)
+  $PASS = 0
+ end
+ if (abs(B[1][0] + 1.0) > 0.0001)
+  $PASS = 0
+ end
+ if (abs(B[1][1] - 4.0) > 0.0001)
+  $PASS = 0
+ end
+end
Index: trunk/Ohana/src/opihi/cmd.data/test/svd.sh
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/test/svd.sh	(revision 26891)
+++ trunk/Ohana/src/opihi/cmd.data/test/svd.sh	(revision 26891)
@@ -0,0 +1,91 @@
+
+list tests
+ test1
+ test2
+end
+
+# a very simple diagonal matrix equation
+macro test1
+ $PASS = 1
+ break -auto off
+
+ delete -q A U w Vt W A1
+
+ mcreate A 2 2
+
+ A[0][0] = 1.0
+ A[1][0] = 0.5
+
+ A[0][1] = 0.5
+ A[1][1] = 2.0
+
+ # svd give A = U w Vt, but this function returns V, not Vt:
+ svd A = U w V
+ if (not($STATUS))
+   $PASS = 0
+ end
+
+ # can we recreate A?
+ mcreate W w[] w[]
+ for i 0 w[]
+   W[$i][$i] = w[$i]
+ end
+
+ matrix transpose V to Vt
+ matrix set t1 = W * Vt
+ matrix set A1 = U * t1
+ set dA = A1 - A
+
+ stats -q dA
+ if ($MIN < -0.0001)
+   $PASS = 0
+ end
+ if ($MAX >  0.0001)
+   $PASS = 0
+ end
+end
+
+# a very simple diagonal matrix equation
+macro test2
+ $PASS = 1
+ break -auto off
+
+ delete -q A U w Vt W A1
+
+ mcreate A 2 3
+
+ A[0][0] = 1.0
+ A[1][0] = 0.0
+
+ A[0][1] = 0.0
+ A[1][1] = 2.0
+
+ A[0][2] = 0.0
+ A[1][2] = 1.0
+
+ # svd give A = U w Vt, but this function returns V, not Vt:
+ svd A = U w V
+ if (not($STATUS))
+   $PASS = 0
+ end
+
+ # can we recreate A?
+ mcreate W w[] w[]
+ for i 0 w[]
+   W[$i][$i] = w[$i]
+ end
+
+ matrix transpose V to Vt
+ matrix set t1 = W * Vt
+ matrix set A1 = U * t1
+ set dA = A1 - A
+
+ stats -q dA
+ if ($MIN < -0.0001)
+   $PASS = 0
+ end
+ if ($MAX >  0.0001)
+   $PASS = 0
+ end
+end
+
Index: trunk/Ohana/src/opihi/doc/svdcmp_eispack.f
===================================================================
--- trunk/Ohana/src/opihi/doc/svdcmp_eispack.f	(revision 26891)
+++ trunk/Ohana/src/opihi/doc/svdcmp_eispack.f	(revision 26891)
@@ -0,0 +1,358 @@
+      subroutine svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1)
+c
+      integer i,j,k,l,m,n,ii,i1,kk,k1,ll,l1,mn,nm,its,ierr
+      double precision a(nm,n),w(n),u(nm,n),v(nm,n),rv1(n)
+      double precision c,f,g,h,s,x,y,z,tst1,tst2,scale,pythag
+      logical matu,matv
+c
+c     this subroutine is a translation of the algol procedure svd,
+c     num. math. 14, 403-420(1970) by golub and reinsch.
+c     handbook for auto. comp., vol ii-linear algebra, 134-151(1971).
+c
+c     this subroutine determines the singular value decomposition
+c          t
+c     a=usv  of a real m by n rectangular matrix.  householder
+c     bidiagonalization and a variant of the qr algorithm are used.
+c
+c     on input
+c
+c        nm must be set to the row dimension of two-dimensional
+c          array parameters as declared in the calling program
+c          dimension statement.  note that nm must be at least
+c          as large as the maximum of m and n.
+c
+c        m is the number of rows of a (and u).
+c
+c        n is the number of columns of a (and u) and the order of v.
+c
+c        a contains the rectangular input matrix to be decomposed.
+c
+c        matu should be set to .true. if the u matrix in the
+c          decomposition is desired, and to .false. otherwise.
+c
+c        matv should be set to .true. if the v matrix in the
+c          decomposition is desired, and to .false. otherwise.
+c
+c     on output
+c
+c        a is unaltered (unless overwritten by u or v).
+c
+c        w contains the n (non-negative) singular values of a (the
+c          diagonal elements of s).  they are unordered.  if an
+c          error exit is made, the singular values should be correct
+c          for indices ierr+1,ierr+2,...,n.
+c
+c        u contains the matrix u (orthogonal column vectors) of the
+c          decomposition if matu has been set to .true.  otherwise
+c          u is used as a temporary array.  u may coincide with a.
+c          if an error exit is made, the columns of u corresponding
+c          to indices of correct singular values should be correct.
+c
+c        v contains the matrix v (orthogonal) of the decomposition if
+c          matv has been set to .true.  otherwise v is not referenced.
+c          v may also coincide with a if u is not needed.  if an error
+c          exit is made, the columns of v corresponding to indices of
+c          correct singular values should be correct.
+c
+c        ierr is set to
+c          zero       for normal return,
+c          k          if the k-th singular value has not been
+c                     determined after 30 iterations.
+c
+c        rv1 is a temporary storage array.
+c
+c     calls pythag for  dsqrt(a*a + b*b) .
+c
+c     questions and comments should be directed to burton s. garbow,
+c     mathematics and computer science div, argonne national laboratory
+c
+c     this version dated august 1983.
+c
+c     ------------------------------------------------------------------
+c
+      ierr = 0
+c
+      do 100 i = 1, m
+c
+         do 100 j = 1, n
+            u(i,j) = a(i,j)
+  100 continue
+c     .......... householder reduction to bidiagonal form ..........
+      g = 0.0d0
+      scale = 0.0d0
+      x = 0.0d0
+c
+      do 300 i = 1, n
+         l = i + 1
+         rv1(i) = scale * g
+         g = 0.0d0
+         s = 0.0d0
+         scale = 0.0d0
+         if (i .gt. m) go to 210
+c
+         do 120 k = i, m
+  120    scale = scale + dabs(u(k,i))
+c
+         if (scale .eq. 0.0d0) go to 210
+c
+         do 130 k = i, m
+            u(k,i) = u(k,i) / scale
+            s = s + u(k,i)**2
+  130    continue
+c
+         f = u(i,i)
+         g = -dsign(dsqrt(s),f)
+         h = f * g - s
+         u(i,i) = f - g
+         if (i .eq. n) go to 190
+c
+         do 150 j = l, n
+            s = 0.0d0
+c
+            do 140 k = i, m
+  140       s = s + u(k,i) * u(k,j)
+c
+            f = s / h
+c
+            do 150 k = i, m
+               u(k,j) = u(k,j) + f * u(k,i)
+  150    continue
+c
+  190    do 200 k = i, m
+  200    u(k,i) = scale * u(k,i)
+c
+  210    w(i) = scale * g
+         g = 0.0d0
+         s = 0.0d0
+         scale = 0.0d0
+         if (i .gt. m .or. i .eq. n) go to 290
+c
+         do 220 k = l, n
+  220    scale = scale + dabs(u(i,k))
+c
+         if (scale .eq. 0.0d0) go to 290
+c
+         do 230 k = l, n
+            u(i,k) = u(i,k) / scale
+            s = s + u(i,k)**2
+  230    continue
+c
+         f = u(i,l)
+         g = -dsign(dsqrt(s),f)
+         h = f * g - s
+         u(i,l) = f - g
+c
+         do 240 k = l, n
+  240    rv1(k) = u(i,k) / h
+c
+         if (i .eq. m) go to 270
+c
+         do 260 j = l, m
+            s = 0.0d0
+c
+            do 250 k = l, n
+  250       s = s + u(j,k) * u(i,k)
+c
+            do 260 k = l, n
+               u(j,k) = u(j,k) + s * rv1(k)
+  260    continue
+c
+  270    do 280 k = l, n
+  280    u(i,k) = scale * u(i,k)
+c
+  290    x = dmax1(x,dabs(w(i))+dabs(rv1(i)))
+  300 continue
+c     .......... accumulation of right-hand transformations ..........
+      if (.not. matv) go to 410
+c     .......... for i=n step -1 until 1 do -- ..........
+      do 400 ii = 1, n
+         i = n + 1 - ii
+         if (i .eq. n) go to 390
+         if (g .eq. 0.0d0) go to 360
+c
+         do 320 j = l, n
+c     .......... double division avoids possible underflow ..........
+  320    v(j,i) = (u(i,j) / u(i,l)) / g
+c
+         do 350 j = l, n
+            s = 0.0d0
+c
+            do 340 k = l, n
+  340       s = s + u(i,k) * v(k,j)
+c
+            do 350 k = l, n
+               v(k,j) = v(k,j) + s * v(k,i)
+  350    continue
+c
+  360    do 380 j = l, n
+            v(i,j) = 0.0d0
+            v(j,i) = 0.0d0
+  380    continue
+c
+  390    v(i,i) = 1.0d0
+         g = rv1(i)
+         l = i
+  400 continue
+c     .......... accumulation of left-hand transformations ..........
+  410 if (.not. matu) go to 510
+c     ..........for i=min(m,n) step -1 until 1 do -- ..........
+      mn = n
+      if (m .lt. n) mn = m
+c
+      do 500 ii = 1, mn
+         i = mn + 1 - ii
+         l = i + 1
+         g = w(i)
+         if (i .eq. n) go to 430
+c
+         do 420 j = l, n
+  420    u(i,j) = 0.0d0
+c
+  430    if (g .eq. 0.0d0) go to 475
+         if (i .eq. mn) go to 460
+c
+         do 450 j = l, n
+            s = 0.0d0
+c
+            do 440 k = l, m
+  440       s = s + u(k,i) * u(k,j)
+c     .......... double division avoids possible underflow ..........
+            f = (s / u(i,i)) / g
+c
+            do 450 k = i, m
+               u(k,j) = u(k,j) + f * u(k,i)
+  450    continue
+c
+  460    do 470 j = i, m
+  470    u(j,i) = u(j,i) / g
+c
+         go to 490
+c
+  475    do 480 j = i, m
+  480    u(j,i) = 0.0d0
+c
+  490    u(i,i) = u(i,i) + 1.0d0
+  500 continue
+c     .......... diagonalization of the bidiagonal form ..........
+  510 tst1 = x
+c     .......... for k=n step -1 until 1 do -- ..........
+      do 700 kk = 1, n
+         k1 = n - kk
+         k = k1 + 1
+         its = 0
+c     .......... test for splitting.
+c                for l=k step -1 until 1 do -- ..........
+  520    do 530 ll = 1, k
+            l1 = k - ll
+            l = l1 + 1
+            tst2 = tst1 + dabs(rv1(l))
+            if (tst2 .eq. tst1) go to 565
+c     .......... rv1(1) is always zero, so there is no exit
+c                through the bottom of the loop ..........
+            tst2 = tst1 + dabs(w(l1))
+            if (tst2 .eq. tst1) go to 540
+  530    continue
+c     .......... cancellation of rv1(l) if l greater than 1 ..........
+  540    c = 0.0d0
+         s = 1.0d0
+c
+         do 560 i = l, k
+            f = s * rv1(i)
+            rv1(i) = c * rv1(i)
+            tst2 = tst1 + dabs(f)
+            if (tst2 .eq. tst1) go to 565
+            g = w(i)
+            h = pythag(f,g)
+            w(i) = h
+            c = g / h
+            s = -f / h
+            if (.not. matu) go to 560
+c
+            do 550 j = 1, m
+               y = u(j,l1)
+               z = u(j,i)
+               u(j,l1) = y * c + z * s
+               u(j,i) = -y * s + z * c
+  550       continue
+c
+  560    continue
+c     .......... test for convergence ..........
+  565    z = w(k)
+         if (l .eq. k) go to 650
+c     .......... shift from bottom 2 by 2 minor ..........
+         if (its .eq. 30) go to 1000
+         its = its + 1
+         x = w(l)
+         y = w(k1)
+         g = rv1(k1)
+         h = rv1(k)
+         f = 0.5d0 * (((g + z) / h) * ((g - z) / y) + y / h - h / y)
+         g = pythag(f,1.0d0)
+         f = x - (z / x) * z + (h / x) * (y / (f + dsign(g,f)) - h)
+c     .......... next qr transformation ..........
+         c = 1.0d0
+         s = 1.0d0
+c
+         do 600 i1 = l, k1
+            i = i1 + 1
+            g = rv1(i)
+            y = w(i)
+            h = s * g
+            g = c * g
+            z = pythag(f,h)
+            rv1(i1) = z
+            c = f / z
+            s = h / z
+            f = x * c + g * s
+            g = -x * s + g * c
+            h = y * s
+            y = y * c
+            if (.not. matv) go to 575
+c
+            do 570 j = 1, n
+               x = v(j,i1)
+               z = v(j,i)
+               v(j,i1) = x * c + z * s
+               v(j,i) = -x * s + z * c
+  570       continue
+c
+  575       z = pythag(f,h)
+            w(i1) = z
+c     .......... rotation can be arbitrary if z is zero ..........
+            if (z .eq. 0.0d0) go to 580
+            c = f / z
+            s = h / z
+  580       f = c * g + s * y
+            x = -s * g + c * y
+            if (.not. matu) go to 600
+c
+            do 590 j = 1, m
+               y = u(j,i1)
+               z = u(j,i)
+               u(j,i1) = y * c + z * s
+               u(j,i) = -y * s + z * c
+  590       continue
+c
+  600    continue
+c
+         rv1(l) = 0.0d0
+         rv1(k) = f
+         w(k) = x
+         go to 520
+c     .......... convergence ..........
+  650    if (z .ge. 0.0d0) go to 700
+c     .......... w(k) is made non-negative ..........
+         w(k) = -z
+         if (.not. matv) go to 700
+c
+         do 690 j = 1, n
+  690    v(j,k) = -v(j,k)
+c
+  700 continue
+c
+      go to 1001
+c     .......... set error -- no convergence to a
+c                singular value after 30 iterations ..........
+ 1000 ierr = k
+ 1001 return
+      end
Index: trunk/Ohana/src/opihi/include/data.h
===================================================================
--- trunk/Ohana/src/opihi/include/data.h	(revision 26884)
+++ trunk/Ohana/src/opihi/include/data.h	(revision 26891)
@@ -89,4 +89,7 @@
 int svdcmp (float *a, opihi_flt *w, float *v, int Nx, int Ny);
 
+/* in svdcmp_bond_raw.c */
+int svdcmp_bond_raw(int m, int n, int withu, int withv, double eps, double tol, double **a, double *q, double **u, double **v);
+
 /* mrqmin.c */
 opihi_flt mrqcof (opihi_flt *x, opihi_flt *y, opihi_flt *dy, int Npts, 
@@ -162,3 +165,17 @@
 void FreeBooks (void);
 
+/* hermitian functions */
+double hermitian_polynomial (double x, int order);
+double hermitian_00(double x);
+double hermitian_01(double x);
+double hermitian_02(double x);
+double hermitian_03(double x);
+double hermitian_04(double x);
+double hermitian_05(double x);
+double hermitian_06(double x);
+double hermitian_07(double x);
+double hermitian_08(double x);
+double hermitian_09(double x);
+double hermitian_10(double x);
+
 # endif
Index: trunk/Ohana/src/opihi/include/pclient.h
===================================================================
--- trunk/Ohana/src/opihi/include/pclient.h	(revision 26884)
+++ trunk/Ohana/src/opihi/include/pclient.h	(revision 26891)
@@ -52,4 +52,6 @@
 void FreePclient (void);
 void gotsignal (int signum);
+void pipe_signal (int signum);
+void pipe_signal_clear(void);
 
 # define DTIME(A,B) ((A.tv_sec - B.tv_sec) + 1e-6*(A.tv_usec - B.tv_usec))
Index: trunk/Ohana/src/opihi/lib.data/Makefile
===================================================================
--- trunk/Ohana/src/opihi/lib.data/Makefile	(revision 26884)
+++ trunk/Ohana/src/opihi/lib.data/Makefile	(revision 26891)
@@ -20,4 +20,5 @@
 $(SDIR)/fft.$(ARCH).o			\
 $(SDIR)/svdcmp.$(ARCH).o		\
+$(SDIR)/svdcmp_bond_raw.$(ARCH).o		\
 $(SDIR)/convert.$(ARCH).o		\
 $(SDIR)/bracket.$(ARCH).o		\
@@ -27,4 +28,5 @@
 $(SDIR)/precess.$(ARCH).o		\
 $(SDIR)/starfuncs.$(ARCH).o		\
+$(SDIR)/hermitian.$(ARCH).o		\
 $(SDIR)/gaussian.$(ARCH).o		\
 $(SDIR)/graphtools.$(ARCH).o            \
Index: trunk/Ohana/src/opihi/lib.data/hermitian.c
===================================================================
--- trunk/Ohana/src/opihi/lib.data/hermitian.c	(revision 26891)
+++ trunk/Ohana/src/opihi/lib.data/hermitian.c	(revision 26891)
@@ -0,0 +1,120 @@
+# include "data.h"
+
+double hermitian_polynomial (double x, int order) {
+  double value;
+    switch (order) {
+      case 0:
+	value = hermitian_00(x);
+	break;
+      case 1:
+	value = hermitian_01(x);
+	break;
+      case 2:
+	value = hermitian_02(x);
+	break;
+      case 3:
+	value = hermitian_03(x);
+	break;
+      case 4:
+	value = hermitian_04(x);
+	break;
+      case 5:
+	value = hermitian_05(x);
+	break;
+      case 6:
+	value = hermitian_06(x);
+	break;
+      case 7:
+	value = hermitian_07(x);
+	break;
+      case 8:
+	value = hermitian_08(x);
+	break;
+      case 9:
+	value = hermitian_09(x);
+	break;
+      case 10:
+	value = hermitian_10(x);
+	break;
+      default:
+	value = NAN;
+	break;
+    }
+    return value;
+}
+
+double hermitian_00(double x) {
+    double value;
+    // H_0(x) = 1
+    value = 1;
+    return value;
+}
+double hermitian_01(double x) {
+    double value;
+    // H_1(x) = x
+    value = x;
+    return value;
+}
+double hermitian_02(double x) {
+    double value, x2;
+    // H_2(x) = x^2-1
+    x2 = x*x;
+    value = x2 - 1.0;
+    return value;
+}
+double hermitian_03(double x) {
+    double value, x2;
+    // H_3(x) = x^3-3x
+    x2 = x*x;
+    value = x*(x2 - 3.0);
+    return value;
+}
+double hermitian_04(double x) {
+    double value, x2;
+    // H_4(x) = x^4-6x^2+3
+    x2 = x*x;
+    value = (x2 - 6.0)*x2 + 3.0;
+    return value;
+}
+double hermitian_05(double x) {
+    double value, x2;
+    // H_5(x) = x^5-10x^3+15x
+    x2 = x*x;
+    value = ((x2 - 10.0)*x2 + 15.0)*x;
+    return value;
+}
+double hermitian_06(double x) {
+    double value, x2;
+    // H_6(x) = x^6-15x^4+45x^2-15
+    x2 = x*x;
+    value = (((x2 - 15.0)*x2 + 45.0)*x2) - 15.0;
+    return value;
+}
+double hermitian_07(double x) {
+    double value, x2;
+    // H_7(x) = x^7-21x^5+105x^3-105x
+    x2 = x*x;
+    value = (((x2 - 21.0)*x2+105.0)*x2 - 105.0)*x;
+    return value;
+}
+double hermitian_08(double x) {
+    double value, x2;
+    // H_8(x) = x^8-28x^6+210x^4-420x^2+105
+    x2 = x*x;
+    value = ((((x2 - 28.0)*x2 + 210.0)*x2 - 420.0)*x2 + 105.0);
+    return value;
+}
+double hermitian_09(double x) {
+    double value, x2;
+    // H_9(x) = x^9-36x^7+378x^5-1260x^3+945x
+    x2 = x*x;
+    value = ((((x2 - 36.0)*x2 + 378.0)*x2 - 1260.0)*x2 + 945.0)*x;
+    return value;
+}
+double hermitian_10(double x) {
+    double value, x2;
+    // H_{10}(x) = x^{10}-45x^8+630x^6-3150x^4+4725x^2-945 
+    x2 = x*x;
+    value = (((((x2 - 45.0)*x2 + 630.0)*x2 - 3150.0)*x2 + 4725.0)*x2 - 945.0);
+    return value;
+}
Index: trunk/Ohana/src/opihi/lib.data/svdcmp_bond_raw.c
===================================================================
--- trunk/Ohana/src/opihi/lib.data/svdcmp_bond_raw.c	(revision 26891)
+++ trunk/Ohana/src/opihi/lib.data/svdcmp_bond_raw.c	(revision 26891)
@@ -0,0 +1,244 @@
+/*  svd.c -- Singular value decomposition. Translated to 'C' from the
+ *           original Algol code in "Handbook for Automatic Computation,
+ *           vol. II, Linear Algebra", Springer-Verlag.
+ *
+ *  (C) 2000, C. Bond. All rights reserved.
+ *
+ *  This is almost an exact translation from the original, except that
+ *  an iteration counter is added to prevent stalls. This corresponds
+ *  to similar changes in other translations.
+ *
+ *  Returns an error code = 0, if no errors and 'k' if a failure to
+ *  converge at the 'kth' singular value.
+ * 
+ */
+#include <malloc.h> /* for array allocation */
+#include <math.h>    /* for 'fabs'           */
+
+int svdcmp_bond_raw(int m, int n, int withu, int withv, double eps, double tol, double **a, double *q, double **u, double **v)
+{
+	int i,j,k,l,l1,iter,retval;
+	double c,f,g,h,s,x,y,z;
+	double *e;
+
+	e = (double *)calloc(n,sizeof(double));
+	retval = 0;
+
+/* Copy 'a' to 'u' */    
+	for (i=0;i<m;i++) {
+		for (j=0;j<n;j++)
+			u[i][j] = a[i][j];
+	}
+/* Householder's reduction to bidiagonal form. */
+	g = x = 0.0;    
+	for (i=0;i<n;i++) {
+		e[i] = g;
+		s = 0.0;
+		l = i+1;
+		for (j=i;j<m;j++)
+			s += (u[j][i]*u[j][i]);
+		if (s < tol)
+			g = 0.0;
+		else {
+			f = u[i][i];
+			g = (f < 0) ? sqrt(s) : -sqrt(s);
+			h = f * g - s;
+			u[i][i] = f - g;
+			for (j=l;j<n;j++) {
+				s = 0.0;
+				for (k=i;k<m;k++)
+					s += (u[k][i] * u[k][j]);
+				f = s / h;
+				for (k=i;k<m;k++)
+					u[k][j] += (f * u[k][i]);
+			} /* end j */
+		} /* end s */
+		q[i] = g;
+		s = 0.0;
+		for (j=l;j<n;j++)
+			s += (u[i][j] * u[i][j]);
+		if (s < tol)
+			g = 0.0;
+		else {
+			f = u[i][i+1];
+			g = (f < 0) ? sqrt(s) : -sqrt(s);
+			h = f * g - s;
+			u[i][i+1] = f - g;
+			for (j=l;j<n;j++) 
+				e[j] = u[i][j]/h;
+			for (j=l;j<m;j++) {
+				s = 0.0;
+				for (k=l;k<n;k++) 
+					s += (u[j][k] * u[i][k]);
+				for (k=l;k<n;k++)
+					u[j][k] += (s * e[k]);
+			} /* end j */
+		} /* end s */
+		y = fabs(q[i]) + fabs(e[i]);                         
+		if (y > x)
+			x = y;
+	} /* end i */
+
+/* accumulation of right-hand transformations */
+	if (withv) {
+		for (i=n-1;i>=0;i--) {
+			if (g != 0.0) {
+				h = u[i][i+1] * g;
+				for (j=l;j<n;j++)
+					v[j][i] = u[i][j]/h;
+				for (j=l;j<n;j++) {
+					s = 0.0;
+					for (k=l;k<n;k++) 
+						s += (u[i][k] * v[k][j]);
+					for (k=l;k<n;k++)
+						v[k][j] += (s * v[k][i]);
+
+				} /* end j */
+			} /* end g */
+			for (j=l;j<n;j++)
+				v[i][j] = v[j][i] = 0.0;
+			v[i][i] = 1.0;
+			g = e[i];
+			l = i;
+		} /* end i */
+ 
+	} /* end withv, parens added for clarity */
+
+/* accumulation of left-hand transformations */
+	if (withu) {
+		for (i=n;i<m;i++) {
+			for (j=n;j<m;j++)
+				u[i][j] = 0.0;
+			u[i][i] = 1.0;
+		}
+	}
+	if (withu) {
+		for (i=n-1;i>=0;i--) {
+			l = i + 1;
+			g = q[i];
+			for (j=l;j<m;j++)  /* upper limit was 'n' */
+				u[i][j] = 0.0;
+			if (g != 0.0) {
+				h = u[i][i] * g;
+				for (j=l;j<m;j++) { /* upper limit was 'n' */
+					s = 0.0;
+					for (k=l;k<m;k++)
+						s += (u[k][i] * u[k][j]);
+					f = s / h;
+					for (k=i;k<m;k++) 
+						u[k][j] += (f * u[k][i]);
+				} /* end j */
+				for (j=i;j<m;j++) 
+					u[j][i] /= g;
+			} /* end g */
+			else {
+				for (j=i;j<m;j++)
+					u[j][i] = 0.0;
+			}
+			u[i][i] += 1.0;
+		} /* end i*/
+	} /* end withu, parens added for clarity */
+
+/* diagonalization of the bidiagonal form */
+	eps *= x;
+	for (k=n-1;k>=0;k--) {
+		iter = 0;
+test_f_splitting:
+		for (l=k;l>=0;l--) {
+			if (fabs(e[l]) <= eps) goto test_f_convergence;
+			if (fabs(q[l-1]) <= eps) goto cancellation;
+		} /* end l */
+
+/* cancellation of e[l] if l > 0 */
+cancellation:
+		c = 0.0;
+		s = 1.0;
+		l1 = l - 1;
+		for (i=l;i<=k;i++) {
+			f = s * e[i];
+			e[i] *= c;
+			if (fabs(f) <= eps) goto test_f_convergence;
+			g = q[i];
+			h = q[i] = sqrt(f*f + g*g);
+			c = g / h;
+			s = -f / h;
+			if (withu) {
+				for (j=0;j<m;j++) {
+					y = u[j][l1];
+					z = u[j][i];
+					u[j][l1] = y * c + z * s;
+					u[j][i] = -y * s + z * c;
+				} /* end j */
+			} /* end withu, parens added for clarity */
+		} /* end i */
+test_f_convergence:
+		z = q[k];
+		if (l == k) goto convergence;
+
+/* shift from bottom 2x2 minor */
+		iter++;
+		if (iter > 30) {
+			retval = k;
+			break;
+		}
+		x = q[l];
+		y = q[k-1];
+		g = e[k-1];
+		h = e[k];
+		f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2*h*y);
+		g = sqrt(f*f + 1.0);
+		f = ((x-z)*(x+z) + h*(y/((f<0)?(f-g):(f+g))-h))/x;
+/* next QR transformation */
+		c = s = 1.0;
+		for (i=l+1;i<=k;i++) {
+			g = e[i];
+			y = q[i];
+			h = s * g;
+			g *= c;
+			e[i-1] = z = sqrt(f*f+h*h);
+			c = f / z;
+			s = h / z;
+			f = x * c + g * s;
+			g = -x * s + g * c;
+			h = y * s;
+			y *= c;
+			if (withv) {
+				for (j=0;j<n;j++) {
+					x = v[j][i-1];
+					z = v[j][i];
+					v[j][i-1] = x * c + z * s;
+					v[j][i] = -x * s + z * c;
+				} /* end j */
+			} /* end withv, parens added for clarity */
+			q[i-1] = z = sqrt(f*f + h*h);
+			c = f/z;
+			s = h/z;
+			f = c * g + s * y;
+			x = -s * g + c * y;
+			if (withu) {
+				for (j=0;j<m;j++) {
+					y = u[j][i-1];
+					z = u[j][i];
+					u[j][i-1] = y * c + z * s;
+					u[j][i] = -y * s + z * c;
+				} /* end j */
+			} /* end withu, parens added for clarity */
+		} /* end i */
+		e[l] = 0.0;
+		e[k] = f;
+		q[k] = x;
+		goto test_f_splitting;
+convergence:
+		if (z < 0.0) {
+/* q[k] is made non-negative */
+			q[k] = - z;
+			if (withv) {
+				for (j=0;j<n;j++)
+					v[j][k] = -v[j][k];
+			} /* end withv, parens added for clarity */
+		} /* end z */
+	} /* end k */
+	
+	free(e);
+	return retval;
+}
Index: trunk/Ohana/src/opihi/lib.shell/stack_math.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/stack_math.c	(revision 26884)
+++ trunk/Ohana/src/opihi/lib.shell/stack_math.c	(revision 26891)
@@ -628,5 +628,5 @@
 
 # define MS_FUNC(OP) {					\
-    if (V2->type == 's')  {				\
+    if (V2->type == 'S')  {				\
       opihi_flt M2 = V2[0].FltValue;			\
       for (i = 0; i < Nx*Ny; i++, out++, M1++) {	\
@@ -635,5 +635,5 @@
       break;						\
     }							\
-    if (V2->type == 'S')  {				\
+    if (V2->type == 's')  {				\
       opihi_int M2 = V2[0].IntValue;			\
       for (i = 0; i < Nx*Ny; i++, out++, M1++) {	\
@@ -704,5 +704,5 @@
 
 # define SM_FUNC(OP) {					\
-    if (V1->type == 's')  {				\
+    if (V1->type == 'S')  {				\
       opihi_flt M1 = V1[0].FltValue;			\
       for (i = 0; i < Nx*Ny; i++, out++, M2++) {	\
@@ -711,5 +711,5 @@
       break;						\
     }							\
-    if (V1->type == 'S')  {				\
+    if (V1->type == 's')  {				\
       opihi_int M1 = V1[0].IntValue;			\
       for (i = 0; i < Nx*Ny; i++, out++, M2++) {	\
Index: trunk/Ohana/src/opihi/pclient/ChildOps.c
===================================================================
--- trunk/Ohana/src/opihi/pclient/ChildOps.c	(revision 26884)
+++ trunk/Ohana/src/opihi/pclient/ChildOps.c	(revision 26891)
@@ -1,6 +1,6 @@
 # include "pclient.h"
-#include <sys/ioctl.h>
-#include <sys/types.h>
-#include <unistd.h>
+# include <sys/ioctl.h>
+# include <sys/types.h>
+# include <unistd.h>
 // #include <stropts.h>
 
@@ -67,4 +67,5 @@
     case 0:   /* pipe is closed */
       /** change child state? **/
+      pipe_signal_clear();
       break;
     default:  /* data in pipe */
@@ -81,4 +82,5 @@
     case 0:   /* pipe is closed */
       /** change child state? **/
+      pipe_signal_clear();
       break;
     default:  /* data in pipe */
Index: trunk/Ohana/src/opihi/pclient/pclient.c.in
===================================================================
--- trunk/Ohana/src/opihi/pclient/pclient.c.in	(revision 26884)
+++ trunk/Ohana/src/opihi/pclient/pclient.c.in	(revision 26891)
@@ -46,5 +46,15 @@
 
   /* ignore the history file.  to change this, see, eg, mana.c */
-  signal (SIGPIPE, gotsignal);
+
+  /* We do not want pclient to exit just because the connection to the 
+     child process is broken.  However, if the connection to the parent 
+     process is broken, we should exit.  To handle this, we define 
+     pipe_signal() which tracks the number of times it has been called.
+     Every time an I/O operation on the child file descriptor raises this signal, 
+     we need to manage the crashed child and then clear the count.  Signals raised by 
+     I/O to the parent do not result in clearing the count.  If the count passes a maximum,
+     we exit pclient */
+     
+  signal (SIGPIPE, pipe_signal);
   signal (SIGTSTP, gotsignal);
   signal (SIGTTIN, gotsignal);
@@ -67,4 +77,19 @@
 }
 
+static int Npipe = 0;
+
+void pipe_signal (int signum) {
+  Npipe ++;
+  if (Npipe > 100) {
+    exit (5);
+  }
+  gprint (GP_ERR, "closed pipe : did child process exit unexpectedly?\n");
+  return;
+}
+
+void pipe_signal_clear(void) {
+  Npipe = 0;
+}
+
 void gotsignal (int signum) {
   gprint (GP_ERR, "got signal : %d\n", signum);