# include "data.h"

int integrate (int argc, char **argv) {
  
  int i, N, VERBOSE, isFloat;
  opihi_flt *Yf;
  opihi_int *Yi;
  double start, end, value, range, dvalue;
  Vector *vecx, *vecy;

  VERBOSE = FALSE;
  if ((N = get_argument (argc, argv, "-v"))) {
    VERBOSE = TRUE;
    remove_argument (N, &argc, argv);
  }

  if (argc != 5) {
    gprint (GP_ERR, "USAGE: integrate <x> <y> start end\n");
    return (FALSE);
  }

  if ((vecx = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((vecy = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);

  start = atof (argv[3]);
  end   = atof (argv[4]);

  isFloat = vecy[0].type == OPIHI_FLT;
  Yf = vecy[0].elements.Flt;
  Yi = vecy[0].elements.Int;

  value = 0;
  range = 0;

  if (vecx[0].type == OPIHI_FLT) {
    opihi_flt *X = vecx[0].elements.Flt;
    for (i = 0; i < vecx[0].Nelements-1; i++, X++, Yi++, Yf++) {
      if ((*X >= start) && (*X <= end)) {
	dvalue = (isFloat) ? *Yf : *Yi;
	value += dvalue * (X[1] - X[0]);
	range += (X[1] - X[0]);
      }
    }      
  } else {
    opihi_int *X = vecx[0].elements.Int;
    for (i = 0; i < vecx[0].Nelements-1; i++, X++, Yi++, Yf++) {
      if ((*X >= start) && (*X <= end)) {
	dvalue = (isFloat) ? *Yf : *Yi;
	value += dvalue * (X[1] - X[0]);
	range += (X[1] - X[0]);
      }
    }      
  }

  set_variable ("sum", value); 
  set_variable ("range", range);
  if (VERBOSE) gprint (GP_LOG, "sum: %f\n", value);
  return (TRUE);
}

