Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/Makefile	(revision 26462)
+++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/Makefile	(revision 26463)
@@ -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	\
Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian1d.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian1d.c	(revision 26463)
+++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian1d.c	(revision 26463)
@@ -0,0 +1,171 @@
+# include "astro.h"
+
+int hermitian1d (int argc, char **argv) {
+  
+  int i, order, Nvec, Poly;
+  int *inI;
+  float *out, inF;
+  float mean, sigma;
+  Vector *xvec, *yvec;
+
+  Poly = FALSE;
+  if ((N = get_argument (argc, argv, "-poly"))) {
+    Poly = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 6) {
+    gprint (GP_ERR, "USAGE: hermitian1d (x) (y) (mean) (sigma) (order)\n");
+    return (FALSE);
+  }
+
+  /* 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]);
+
+  out = (float *) yvec[0].elements.Flt;
+  inF = (float *) xvec[0].elements.Flt;
+  inI = (int *)   xvec[0].elements.Int;
+
+  // 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 *= exp (-0.25*x*x);
+    }
+  }
+  return (TRUE);
+}
+
+  // {\psi}_n(x) = \frac{1}{\sqrt{n! \, 2^n\sqrt{\pi}}}\, \mathrm{e}^{-x^2/2}H_n(x).\,\! 
+
+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);
+    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: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian2d.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian2d.c	(revision 26463)
+++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian2d.c	(revision 26463)
@@ -0,0 +1,68 @@
+# include "astro.h"
+
+int hermitian2d (int argc, char **argv) {
+  
+  int i, j, Nx, Ny, N;
+  float *in;
+  double Sig_x, Sig_y, Theta;
+  double root1, root2, R, A1, A2, A3;
+  double Sx, Sy, Sxy;
+  double x, y, r, f, Xo, Yo;
+  Buffer *buf;
+
+  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) {
+    gprint (GP_ERR, "USAGE: hermitian2d (buffer) (sigma) Xorder Yorder\n");
+    return (FALSE);
+  }
+
+  /* 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]);
+
+  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) {
+	value *= xf*yf;
+      }
+
+      *in += xf*yf;
+    }
+  }
+
+  return (TRUE);
+}
Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/init.c	(revision 26462)
+++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/init.c	(revision 26463)
@@ -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 **));
@@ -173,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"},