# include "data.h"

int medacc (int argc, char **argv) {
  
  int i, j, Nbins, Nvalues, N, N0, N1;
  double start, end, delta, k0, k1, fn;
  float *V, *K, *V1, *K1, *O, *tmpvec, *tmpkey;
  Vector *val, *key, *out;

  if ((argc != 6) && (argc != 7)) {
    gprint (GP_ERR, "USAGE: medacc <value> <vector> <key> start end [delta]\n");
    return (FALSE);
  }

  if ((val = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((key = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if (val[0].Nelements != key[0].Nelements) {
    gprint (GP_ERR, "key and value don't match\n");
    return (FALSE);
  }
  if ((out = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);

  start = atof (argv[4]);
  end   = atof (argv[5]);
  if (argc == 7) 
    delta = atof (argv[6]);
  else 
    delta = 1;

  delta = fabs (delta);
  if (end - start < 0) {
    delta = -1.0 * delta;
  }
  Nbins = 1 + (int)((end - start) / delta);

  out[0].Nelements = Nbins;
  REALLOCATE (out[0].elements, float, out[0].Nelements);
  bzero (out[0].elements, sizeof(float)*out[0].Nelements);

  /* copy vec and key to temp vectors */
  ALLOCATE (tmpvec, float, val[0].Nelements);
  ALLOCATE (tmpkey, float, val[0].Nelements);

  V = val[0].elements;
  K = key[0].elements;
  V1 = tmpvec;
  K1 = tmpkey;
  Nvalues = val[0].Nelements;
  for (i = 0; i < Nvalues; i++, V++, K++, V1++, K1++) {
    *V1 = *V;
    *K1 = *K;
  }      

  /* sort vec and key by key */
  fsortpair (tmpkey, tmpvec, Nvalues);

  O = out[0].elements;
  /* find the start and end key for each range */
  N0 = 0;
  N1 = 0;
  for (i = 0; i < Nbins; i++) {
    k0 = i*delta + start;
    k1 = (i+1)*delta + start;
    for (j = N1; (j < Nvalues) && (tmpkey[j] < k0); j++);
    N0 = j;
    for (j = N0; (j < Nvalues) && (tmpkey[j] < k1); j++);
    N1 = j;
    N = N1 - N0;
    fsort (&tmpvec[N0], N);
    fn = O[i] = 0;
    for (j = N0 + 0.25*N; j < N0 + 0.75*N; j++) {
      O[i] += tmpvec[j];
      fn += 1.0;
    }
    if (fn > 0) O[i] /= fn;
  }      
  
  free (tmpvec);
  free (tmpkey);

  return (TRUE);
}

