/*
  loccums.h

  C template for loccum.c

  data-to-data functions

  $Revision: 1.5 $ $Date: 2013/09/18 04:28:45 $

  macros: 

  FNAME    function name
  NULVAL   initial value (empty sum = 0, empty product = 1)
  INC(A,B) increment operation A += B or A *= B

*/

void FNAME(n, x, y, v,
	   nr, rmax, 
	   ans)
     /* inputs */
     int *n, *nr;
     double *x, *y, *v;
     double *rmax;
     /* output */
     double *ans;  /* matrix of column vectors of functions for each point */
{
  int N, Nr, Nans;
  double Rmax;

  int i, j, k, kmin, maxchunk, columnstart;
  double Rmax2, rstep, xi, yi;
  double dx, dy, dx2, d2, d, contrib;

  N    = *n;
  Nr   = *nr;
  Rmax = *rmax;

  if(N == 0) 
    return;

  rstep = Rmax/(Nr-1);
  Rmax2 = Rmax * Rmax;
  Nans  = Nr * N;

  /* initialise products to 1 */
  OUTERCHUNKLOOP(k, Nans, maxchunk, 8196) {
    R_CheckUserInterrupt();
    INNERCHUNKLOOP(k, Nans, maxchunk, 8196) {
      ans[k] = NULVAL;
    }
  }
   
  OUTERCHUNKLOOP(i, N, maxchunk, 8196) {
    R_CheckUserInterrupt();
    INNERCHUNKLOOP(i, N, maxchunk, 8196) {
      xi = x[i];
      yi = y[i];
      columnstart = Nr * i; /* start position for f_i(.) in 'ans' */
      /* 
	 process backward until |dx| > Rmax
      */
      if(i > 0) {
	for(j=i-1; j >= 0; j--) {
	  dx = x[j] - xi;
	  dx2 = dx * dx;
	  if(dx2 > Rmax2) 
	    break;
	  dy = y[j] - yi;
	  d2 = dx2 + dy * dy;
	  if(d2 <= Rmax2) {
	    d = sqrt(d2);
	    kmin = (int) ceil(d/rstep);
	    if(kmin < Nr) {
	      contrib = v[j];
	      for(k = kmin; k < Nr; k++) 
		INC(ans[columnstart + k] , contrib);
	    }
	  }
	}
      }
      /* 
	 process forward until |dx| > Rmax
      */
      if(i < N - 1) {
	for(j=i+1; j < N; j++) {
	  dx = x[j] - xi;
	  dx2 = dx * dx;
	  if(dx2 > Rmax2) 
	    break;
	  dy = y[j] - yi;
	  d2 = dx2 + dy * dy;
	  if(d2 <= Rmax2) {
	    d = sqrt(d2);
	    kmin = (int) ceil(d/rstep);
	    if(kmin < Nr) {
	      contrib = v[j];
	      for(k = kmin; k < Nr; k++) 
		INC(ans[columnstart + k] , contrib);
	    }
	  }
	}
      }
    }   
  }
}