# 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];
    if (y0 == y1) {
      Xvalue = 0.5*(x0 + x1);
    } else {
      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];
    if (y0 == y1) {
      Xvalue = 0.5*(x0 + x1);
    } else {
      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);
}
