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
+