Index: trunk/Ohana/src/opihi/cmd.data/Makefile
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 40382)
+++ trunk/Ohana/src/opihi/cmd.data/Makefile	(revision 40383)
@@ -166,4 +166,5 @@
 $(SRC)/vclip.$(ARCH).o		   \
 $(SRC)/vgauss.$(ARCH).o            \
+$(SRC)/vlorentz.$(ARCH).o          \
 $(SRC)/vellipse.$(ARCH).o          \
 $(SRC)/vmaxwell.$(ARCH).o          \
@@ -178,5 +179,5 @@
 $(SRC)/vstats.$(ARCH).o		   \
 $(SRC)/xsection.$(ARCH).o          \
-$(SRC)/vsh.$(ARCH).o	   \
+$(SRC)/vsh.$(ARCH).o	           \
 $(SRC)/vshfit.$(ARCH).o	   \
 $(SRC)/shterms.$(ARCH).o	   \
Index: trunk/Ohana/src/opihi/cmd.data/init.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/init.c	(revision 40382)
+++ trunk/Ohana/src/opihi/cmd.data/init.c	(revision 40383)
@@ -153,4 +153,5 @@
 int vgrid            PROTO((int, char **));
 int vgauss           PROTO((int, char **));
+int vlorentz         PROTO((int, char **));
 int vellipse         PROTO((int, char **));
 int vmaxwell         PROTO((int, char **));
@@ -338,4 +339,5 @@
   {1, "vtype",        vtype,            "return the vector type (FLT or INT)"},
   {1, "vgauss",       vgauss,           "fit a Gaussian to a vector"},
+  {1, "vlorentz",     vlorentz,         "fit a Lorentzian to a vector"},
   {1, "vellipse",     vellipse,         "fit a Ellipse to a vector pair"},
   {1, "vgrid",        vgrid,            "generate an image from a triplet of vectors"},
Index: trunk/Ohana/src/opihi/cmd.data/test/vlorentz.sh
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/test/vlorentz.sh	(revision 40383)
+++ trunk/Ohana/src/opihi/cmd.data/test/vlorentz.sh	(revision 40383)
@@ -0,0 +1,175 @@
+
+list tests
+ test1
+ test2
+ test3
+end
+
+macro test1
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+
+ $xo = 0.0
+ $go = 2.0
+ $Io = 10.0
+ $so = 0.0
+
+ set ytru = $Io * $go / ((x - $xo)^2 + $go^2 / 4.0) / 2 / $M_PI + $so
+ gaussdev ystd ytru[] 0.0 1.0
+ # set yobs = ytru + ystd * sqrt(abs(ytru))
+ set yobs = ytru + ystd * 0.05
+
+ $C0 = 0.5
+ $C1 = 3
+ $C2 = 20
+ $C3 = 0
+ # set dy = sqrt(abs(yobs))
+ set dy = 0.05 + zero(yobs)
+
+ vlorentz x yobs dy yfit
+ lim x yobs; clear; box; plot -x 1 -c black x yobs; plot -c red x 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
+
+macro test1a
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+
+ $xo = 0.0
+ $go = 2.0
+ $Io = 10.0
+ $so = 0.0
+
+ set ytru = $Io * $go / ((x - $xo)^2 + $go^2 / 4.0) / 2 / $M_PI + $so
+ gaussdev ystd ytru[] 0.0 1.0
+ # set yobs = ytru + ystd * sqrt(abs(ytru))
+ set yobs = ytru + ystd * 0.05
+ set dy = 0.05 + zero(yobs)
+
+ delete rC0 rC1 rC2 rC3 vC0 vC1 vC2 vC3
+
+ for i 0 500
+
+   $C0 = -1 + 2*rnd(0)
+   $C1 = 0.5 + 5*rnd(0)
+   $C2 = 2.0 + 30*rnd(0)
+   $C3 = -10 + 20*rnd(0)
+   
+   concat $C0 rC0
+   concat $C1 rC1
+   concat $C2 rC2
+   concat $C3 rC3
+
+   vlorentz -q x yobs dy yfit
+   # lim x yobs; clear; box; plot -x 1 -c black x yobs; plot -c red x yfit
+
+   concat $C0 vC0
+   concat $C1 vC1
+   concat $C2 vC2
+   concat $C3 vC3
+ end
+
+ 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
+ # lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit
+ # cursor
+
+ 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
+
+ $C0o = 0.5
+ $C1o = 2
+ $C2o = 900
+ $C3o = 1
+
+ create x -10 10 0.1
+ set y = $C3o + $C2o * exp(-0.5*(x - $C0o)^2/$C1o^2)
+
+ gaussdev dY y[] 0.0 1.0
+ set dy = dY * sqrt(y)
+ set y = y + dy
+
+ $C0 = $C0o + 2
+ $C1 = $C1o + 2
+ $C2 = $C2o + 50
+ $C3 = $C3o + 10
+
+ vgauss -q x y dy yfit
+ # lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit
+ # cursor
+
+ $dS = 3.0 / sqrt(y[])
+
+ if (abs($C0 - $C0o) > $dS)
+   $PASS = 0
+ end
+ if (abs($C1 - $C1o)/$C1o > $dS)
+   $PASS = 0
+ end
+ if (abs($C2 - $C2o)/$C2o > $dS)
+   $PASS = 0
+ end
+ if (abs($C3 - $C3o) > $dS)
+   $PASS = 0
+ end
+end
Index: trunk/Ohana/src/opihi/cmd.data/vlorentz.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/vlorentz.c	(revision 40383)
+++ trunk/Ohana/src/opihi/cmd.data/vlorentz.c	(revision 40383)
@@ -0,0 +1,146 @@
+# include "data.h"
+
+/* local private functions */
+opihi_flt florentzOD (opihi_flt, opihi_flt *, int, opihi_flt *);
+
+# 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); \
+  free (c);
+
+int vlorentz (int argc, char **argv) {
+
+  int i, N, Npts, Npar, Quiet;
+  opihi_flt par[4], *v1, *v2, *dy, **covar;
+  opihi_flt 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);
+  }
+
+//   Guess = FALSE;
+//   if ((N = get_argument (argc, argv, "-guess"))) {
+//     Guess = TRUE;
+//     remove_argument (N, &argc, argv);
+//   }
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: vlorentz <x> <y> <dy> (out)\n");
+    gprint (GP_ERR, "  <dy> may be the words 'con[stant]' or 'poi[sson]', in which case the error is constructed'\n");
+    gprint (GP_ERR, " uses guesses: C0 (mean), C1 (sigma), C2 (norm), C3 (sky)\n");
+    return (FALSE);
+  }
+  
+  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((ovec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  CastVector (xvec, OPIHI_FLT);
+  CastVector (yvec, OPIHI_FLT);
+
+  int Nsvec = strlen(argv[3]);
+
+  if ((Nsvec >= 3) && (!strncasecmp ("poisson", argv[3], Nsvec) || !strncasecmp ("constant", argv[3], Nsvec))) {
+    svec = SelectVector ("vlorentz_err", ANYVECTOR, TRUE);
+    if (svec == NULL) return (FALSE);
+    MatchVector (svec, xvec, OPIHI_FLT);
+    v1 = svec[0].elements.Flt;
+    v2 = yvec[0].elements.Flt;
+    if (!strncasecmp ("poisson", argv[3], Nsvec)) {
+      for (i = 0; i < svec[0].Nelements; i++) {
+	v1[i] = (isfinite(v2[i]) && (v2[i] > 0)) ? sqrt(v2[i]) : 1.0;
+      }
+    } else {
+      for (i = 0; i < svec[0].Nelements; i++) {
+	v1[i] = 1.0;
+      }
+    }
+  } else {
+    if ((svec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  }
+
+  CastVector (svec, OPIHI_FLT);
+
+  Npts = xvec[0].Nelements;
+  ResetVector (ovec, OPIHI_FLT, Npts);
+
+  ALLOCATE (dy, opihi_flt, Npts);
+
+  GET_VAR (par[0], "C0");
+  GET_VAR (par[1], "C1");
+  GET_VAR (par[2], "C2");
+  GET_VAR (par[3], "C3");
+  Npar = 4;
+
+  // mrqmin takes the inverse variance (do not generate NANs)
+  v1 = svec[0].elements.Flt;
+  v2 = dy;
+  for (i = 0; i < Npts; i++, v1++, v2++) {
+      *v2 = (*v1 == 0.0) ? 0.0 : 1.0 / (*v1 * *v1);
+  } 
+  
+  ochisq = mrqinit (xvec[0].elements.Flt, yvec[0].elements.Flt, dy, Npts, par, Npar, florentzOD, !Quiet);
+  dchisq = ochisq + 2*Npts;
+
+  for (i = 0; (i < 20) && ((dchisq > 0.1*(Npts - Npar)) || (dchisq <= 0.0)); i++) {
+    chisq = mrqmin (xvec[0].elements.Flt, yvec[0].elements.Flt, dy, Npts, par, Npar, florentzOD, !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.Flt[i] = florentzOD (xvec[0].elements.Flt[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]));
+  }
+
+  free (dy);
+  mrqfree (Npar);
+  return (TRUE);
+}
+
+/* pars: x_o, g, Io, sky */
+opihi_flt florentzOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
+  OHANA_UNUSED_PARAM(Npar);
+
+  opihi_flt z, f, L;
+
+  opihi_flt xo = par[0];
+  opihi_flt go = par[1];
+  opihi_flt Io = par[2];
+  opihi_flt So = par[3];
+
+  z = SQ(x - xo) + SQ(go)/4.0;
+  f = 1.0 / z;
+  L = 0.5*Io*go*f/M_PI + So;
+
+  dpar[0] = Io*go*SQ(f)*(x - xo) / M_PI;
+  dpar[1] = 0.5*Io*(f - 0.5*SQ(f)*SQ(go)) / M_PI;
+  dpar[2] = 0.5*go*f / M_PI;
+  dpar[3] = 1;
+  
+  return (L);
+}