Index: /trunk/Ohana/src/opihi/cmd.data/test/vgauss.sh
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/test/vgauss.sh	(revision 10069)
+++ /trunk/Ohana/src/opihi/cmd.data/test/vgauss.sh	(revision 10069)
@@ -0,0 +1,99 @@
+
+list tests
+ test1
+ test2
+ test3
+end
+
+macro test1
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+ set y = 5 * exp(-0.5*x^2/3^2)
+
+ $C0 = 0.5
+ $C1 = 2
+ $C2 = 10
+ $C3 = 1
+ set dy = sqrt(y)
+
+ vgauss -q x y dy yfit
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 5.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.01)
+   $PASS = 0
+ end
+end
+
+# noise of 0.1
+macro test2
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+ set y = 1000 * exp(-0.5*x^2/3^2)
+ set dy = (rnd(y) - 0.5)/0.5
+ set y = y + dy
+
+ $C0 = 0.5
+ $C1 = 2
+ $C2 = 900
+ $C3 = 1
+
+ vgauss -q x y dy yfit
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 1000.0) > 0.1)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.1)
+   $PASS = 0
+ end
+end
+
+# poisson-distributed noise
+macro test3
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+ set y = 1000 * exp(-0.5*x^2/3^2)
+
+ gaussdev dY y[] 0.0 1.0
+ set dy = dY * sqrt(y)
+ set y = y + dy
+
+ $C0 = 0.5
+ $C1 = 2
+ $C2 = 900
+ $C3 = 1
+
+ vgauss -q x y dy yfit
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 1000.0) > 1)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.2)
+   $PASS = 0
+ end
+end
Index: /trunk/Ohana/src/opihi/cmd.data/test/vmaxwell.sh
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/test/vmaxwell.sh	(revision 10069)
+++ /trunk/Ohana/src/opihi/cmd.data/test/vmaxwell.sh	(revision 10069)
@@ -0,0 +1,103 @@
+
+list tests
+# test1
+ test2
+# test3
+end
+
+macro test1
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+ set y = 5 * exp(-0.5*x^2/3^2)
+
+ $C0 = 0.5
+ $C1 = 2
+ $C2 = 10
+ $C3 = 1
+ set dy = sqrt(y)
+
+ vgauss -q x y dy yfit
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 5.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.01)
+   $PASS = 0
+ end
+end
+
+# noise of 0.1
+macro test2
+ $PASS = 1
+ break -auto off
+
+ create x 0 1000 1
+ set y = 1000 * (x-200)^2 * exp(-0.5*(x-300)^2/30^2)
+ set dy = (rnd(y) - 0.5)/0.5
+ set y = y + dy
+
+ lim x y; clear; box; plot -x 1 -c black x y
+
+ $C0 = 250
+ $C1 = 25
+ $C2 = 900
+ $C3 = 1
+ $C4 = 190
+
+ vmaxwell x y dy yfit
+ plot x yfit -c red
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 1000.0) > 0.1)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.1)
+   $PASS = 0
+ end
+end
+
+# poisson-distributed noise
+macro test3
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+ set y = 1000 * exp(-0.5*x^2/3^2)
+
+ gaussdev dY y[] 0.0 1.0
+ set dy = dY * sqrt(y)
+ set y = y + dy
+
+ $C0 = 0.5
+ $C1 = 2
+ $C2 = 900
+ $C3 = 1
+
+ vgauss -q x y dy yfit
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 1000.0) > 1)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.2)
+   $PASS = 0
+ end
+end
Index: /trunk/Ohana/src/opihi/cmd.data/vmaxwell.c
===================================================================
--- /trunk/Ohana/src/opihi/cmd.data/vmaxwell.c	(revision 10069)
+++ /trunk/Ohana/src/opihi/cmd.data/vmaxwell.c	(revision 10069)
@@ -0,0 +1,106 @@
+# include "data.h"
+
+/* local private functions */
+float fmaxwellOD (float, float *, int, float *);
+
+# define GET_VAR(V,A) \
+  c = get_variable (A); \
+  if (c == NULL) { \
+    gprint (GP_ERR, "missing fit parameter A\n"); \
+    return (FALSE); \
+  } \
+  V = atof (c);
+
+int vmaxwell (int argc, char **argv) {
+
+  int i, N, Npts, Npar, Quiet;
+  float par[5], *v1, *v2, *dy, **covar;
+  float chisq, ochisq, dchisq;
+  Vector *xvec, *yvec, *svec, *ovec;
+  char *c, name[16];
+
+  Quiet = FALSE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-quiet"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: vmaxwell <x> <y> <dy> (out)\n");
+    gprint (GP_ERR, " uses guesses: C0 (mean), C1 (sigma), C2 (norm), C3 (sky), C4 (ref)\n");
+    return (FALSE);
+  }
+  
+  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((svec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((ovec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  Npts = xvec[0].Nelements;
+  ALLOCATE (dy, float, Npts);
+  REALLOCATE (ovec[0].elements, float, Npts);
+
+  GET_VAR (par[0], "C0");
+  GET_VAR (par[1], "C1");
+  GET_VAR (par[2], "C2");
+  GET_VAR (par[3], "C3");
+  GET_VAR (par[4], "C4");
+  Npar = 5;
+  /* careful of variable renomalization */
+
+  v1 = svec[0].elements;
+  v2 = dy;
+  for (i = 0; i < Npts; i++, v1++, v2++) *v2 = 1.0 / (*v1 * *v1);
+  
+  ochisq = mrqinit (xvec[0].elements, yvec[0].elements, dy, Npts, par, Npar, fmaxwellOD, !Quiet);
+  dchisq = ochisq + 2*Npts;
+
+  for (i = 0; (i < 20) && ((dchisq > 0.1*(Npts - Npar)) || (dchisq <= 0.0)); i++) {
+    chisq = mrqmin (xvec[0].elements, yvec[0].elements, dy, Npts, par, Npar, fmaxwellOD, !Quiet);
+    dchisq = ochisq - chisq;
+    ochisq = chisq;
+    if (!Quiet) gprint (GP_ERR, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar);
+  }  
+  if (!Quiet) gprint (GP_ERR, "%d iterations\n", i); 
+
+  for (i = 0; i < Npts; i++) {
+    ovec[0].elements[i] = fmaxwellOD (xvec[0].elements[i], par, Npar, dy);
+  }
+  ovec[0].Nelements = Npts;
+  /* set output *before* variable renomalization */
+
+  covar = mrqcovar (Npar);
+  for (i = 0; i < Npar; i++) {
+    sprintf (name, "C%d", i);
+    set_variable (name, par[i]);
+    if (!Quiet) gprint (GP_ERR, "%d  %f  %f\n", i, par[i], sqrt(covar[i][i]));
+    sprintf (name, "dC%d", i);
+    set_variable (name, sqrt(covar[i][i]));
+  }
+
+  mrqfree (Npar);
+  return (TRUE);
+}
+
+/* pars: x_o, -0.5/sigma^2, I, back, ref */
+float fmaxwellOD (float x, float *par, int Npar, float *dpar) {
+
+  float z, r, f;
+
+  z = (x - par[0])/par[1];
+  r = SQ(x - par[4])*exp (-0.5*SQ(z));
+  f = par[2]*r + par[3];
+
+  dpar[0] = par[2]*r*z/par[1];
+  dpar[1] = par[2]*r*z*z/par[1];
+  dpar[2] = r;
+  dpar[3] = 1;
+  dpar[4] = -par[2]*(x - par[4])*exp(-0.5*SQ(z));
+  
+  return (f);
+
+}