Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian1d.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian1d.c	(revision 26466)
+++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian1d.c	(revision 26467)
@@ -1,10 +1,11 @@
-# include "astro.h"
+# include "data.h"
 
 int hermitian1d (int argc, char **argv) {
   
-  int i, order, Nvec, Poly;
-  int *inI;
-  float *out, inF;
-  float mean, sigma;
+  int i, N, order, Nvec, Poly;
+  opihi_int *inI;
+  opihi_flt *out, *inF;
+  opihi_flt mean, sigma, x;
+  double nf, norm;
   Vector *xvec, *yvec;
 
@@ -15,8 +16,5 @@
   }
 
-  if (argc != 6) {
-    gprint (GP_ERR, "USAGE: hermitian1d (x) (y) (mean) (sigma) (order)\n");
-    return (FALSE);
-  }
+  if (argc != 6) goto usage;
 
   /* select input / output buffers */
@@ -27,11 +25,16 @@
 
   /* gaussian parameters */
-  mean = atof (argv[3]);
+  mean  = atof (argv[3]);
   sigma = atof (argv[4]);
   order = atoi (argv[5]);
+  if (order < 0) goto usage;
+  if (order > 10) goto usage;
 
-  out = (float *) yvec[0].elements.Flt;
-  inF = (float *) xvec[0].elements.Flt;
-  inI = (int *)   xvec[0].elements.Int;
+  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?
@@ -44,128 +47,15 @@
     *out = hermitian_polynomial (x, order);
     if (!Poly) {
-      *out *= exp (-0.25*x*x);
+      *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).\,\! 
-
-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;
-}
+// {\psi}_n(x) = \frac{1}{\sqrt{n! \, 2^n\sqrt{\pi}}}\, \mathrm{e}^{-x^2/2}H_n(x).\,\! 
Index: /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian2d.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian2d.c	(revision 26466)
+++ /branches/eam_branches/20091201/Ohana/src/opihi/cmd.data/hermitian2d.c	(revision 26467)
@@ -3,13 +3,13 @@
 int hermitian2d (int argc, char **argv) {
   
-  int i, j, Nx, Ny, N;
+  int i, j, N, Nx, Ny, Xorder, Yorder, Poly;
   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;
+  double Xo, Yo, sigma;
+  double x, y, xf, yf, value, norm, weight;
+  double nf, mf;
   Buffer *buf;
 
-  Poly = FALSE;
+  // optionally return the Hermitian Polynomial or the Hermitian Function
+  Poly = FALSE; 
   if ((N = get_argument (argc, argv, "-poly"))) {
     Poly = TRUE;
@@ -26,8 +26,5 @@
   }
 
-  if (argc > 5) {
-    gprint (GP_ERR, "USAGE: hermitian2d (buffer) (sigma) Xorder Yorder\n");
-    return (FALSE);
-  }
+  if (argc != 5) goto usage;
 
   /* select input / output buffers */
@@ -43,4 +40,12 @@
   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;
@@ -58,11 +63,17 @@
       value = xf*yf;
       if (!Poly) {
-	value *= xf*yf;
+	weight = norm*exp(-0.25*(x*x + y*y));
+	value *= weight;
       }
 
-      *in += xf*yf;
+      *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: /branches/eam_branches/20091201/Ohana/src/opihi/include/data.h
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/opihi/include/data.h	(revision 26466)
+++ /branches/eam_branches/20091201/Ohana/src/opihi/include/data.h	(revision 26467)
@@ -165,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: /branches/eam_branches/20091201/Ohana/src/opihi/lib.data/Makefile
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/opihi/lib.data/Makefile	(revision 26466)
+++ /branches/eam_branches/20091201/Ohana/src/opihi/lib.data/Makefile	(revision 26467)
@@ -28,4 +28,5 @@
 $(SDIR)/precess.$(ARCH).o		\
 $(SDIR)/starfuncs.$(ARCH).o		\
+$(SDIR)/hermitian.$(ARCH).o		\
 $(SDIR)/gaussian.$(ARCH).o		\
 $(SDIR)/graphtools.$(ARCH).o            \
Index: /branches/eam_branches/20091201/Ohana/src/opihi/lib.data/hermitian.c
===================================================================
--- /branches/eam_branches/20091201/Ohana/src/opihi/lib.data/hermitian.c	(revision 26467)
+++ /branches/eam_branches/20091201/Ohana/src/opihi/lib.data/hermitian.c	(revision 26467)
@@ -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;
+}