Index: /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.data/threshold.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.data/threshold.c	(revision 33601)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/opihi/cmd.data/threshold.c	(revision 33601)
@@ -0,0 +1,88 @@
+# include "data.h"
+
+int threshold (int argc, char **argv) {
+  
+  int N, QUIET;
+  double value;
+  Vector *vecx, *vecy;
+
+  QUIET = FALSE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    QUIET = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: threshold <x> <y> (value)\n");
+    gprint (GP_ERR, "  find the x coordinate at which we pass the specified value\n");
+    gprint (GP_ERR, "  by default, y must be monotonically increasing\n");
+    return (FALSE);
+  }
+  
+  if ((vecx = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((vecy = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  
+  value = atof (argv[3]);
+  
+  // if (argc == 6) {
+  //   start = atof (argv[4]);
+  //   end   = atof (argv[5]);
+  // } else {
+  //   start = (vecx[0].type == OPIHI_FLT) ? vecx[0].elements.Flt[0] : vecx[0].elements.Int[0];
+  //   end   = (vecx[0].type == OPIHI_FLT) ? vecx[0].elements.Flt[vecx[0].Nelements - 1] : vecx[0].elements.Int[vecx[0].Nelements - 1];
+  // }
+
+  int isFltX = (vecx[0].type == OPIHI_FLT);
+  int isFltY = (vecy[0].type == OPIHI_FLT);
+
+  // use bisection to find the value
+  int Nlo, Nhi;
+  int Nelements = vecx[0].Nelements;
+
+  // find the last entry before start
+  Nlo = 0;
+  Nhi = Nelements - 1;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N];
+    if (testval < value) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N, Nelements - 1);
+    }
+  }
+  // v[Nlo] < value <= v[Nhi]
+  for (N = Nlo; N <= Nhi; N++) {
+    double testval = isFltY ? vecy[0].elements.Flt[N] : vecy[0].elements.Int[N];
+    if (testval > value) {
+      Nhi = N;
+      break;
+    }
+  }
+  // v[Nhi] is transition bin
+  
+  double x0, x1, y0, y1, Xvalue;
+  if (Nhi == 0) {
+    // interpolate to value:
+    y0 = isFltY ? vecy[0].elements.Flt[Nhi] : vecy[0].elements.Int[Nhi];
+    y1 = isFltY ? vecy[0].elements.Flt[Nhi+1]   : vecy[0].elements.Int[Nhi+1];
+    x0 = isFltX ? vecx[0].elements.Flt[Nhi] : vecy[0].elements.Int[Nhi];
+    x1 = isFltX ? vecx[0].elements.Flt[Nhi+1]   : vecy[0].elements.Int[Nhi+1];
+    Xvalue = (value - y0) * (x1 - x0) / (y1 - y0) + x0;
+  } else {
+    // interpolate to value:
+    y0 = isFltY ? vecy[0].elements.Flt[Nhi-1] : vecy[0].elements.Int[Nhi-1];
+    y1 = isFltY ? vecy[0].elements.Flt[Nhi]   : vecy[0].elements.Int[Nhi];
+    x0 = isFltX ? vecx[0].elements.Flt[Nhi-1] : vecy[0].elements.Int[Nhi-1];
+    x1 = isFltX ? vecx[0].elements.Flt[Nhi]   : vecy[0].elements.Int[Nhi];
+    Xvalue = (value - y0) * (x1 - x0) / (y1 - y0) + x0;
+  }
+
+  set_variable ("threshval", y1);
+  set_int_variable ("threshbin", Nhi);
+  set_variable ("threshold", Xvalue);
+
+  if (!QUIET) gprint (GP_LOG, "theshold %f (bin %d is %f)\n", Xvalue, Nhi, y1);
+
+  return (TRUE);
+}