#include <math.h>
#include <R.h>
#include <Rmath.h>
#include "geom3.h"
#include "functable.h"

/*
	$Revision: 1.5 $	$Date: 2010/10/24 10:54:43 $

	pair correlation function of 3D point pattern
	(Epanechnikov kernel) 

	pcf3trans	  	translation correction

	pcf3isot		isotropic correction

*/

#define FOURPI (2.0 * M_2PI)


void
pcf3trans(p, n, b, pcf, delta)
     Point *p;
     int n;
     Box *b;
     Ftable *pcf;
     double delta;
{
  register int i, j, l, lmin, lmax;
  register double dx, dy, dz, dist;
  register double  vx, vy, vz, tval;
  Point *ip, *jp;
  double dt, vol, lambda, denom;
  double coef, twocoef, frac, invweight, kernel;

  double sphesfrac(), sphevol();

  /* compute denominator & initialise numerator*/
  vol = (b->x1 - b->x0) * (b->y1 - b->y0) * (b->z1 - b->z0);
  lambda = ((double) n )/ vol;
  denom = lambda * lambda;

  for(l = 0; l < pcf->n; l++) {
    (pcf->denom)[l] = denom;
    (pcf->num)[l]   = 0.0;
  }

  /* spacing of argument in result vector */
  dt = (pcf->t1 - pcf->t0)/(pcf->n - 1);

  /* compute numerator */
  for( i = 0; i < n; i++) {
    ip = p + i;
    for(j = i + 1; j < n; j++) {
      /* compute pairwise distance */
      jp = p + j;
      dx = jp->x - ip->x;
      dy = jp->y - ip->y;
      dz = jp->z - ip->z;
      dist = sqrt(dx * dx + dy * dy + dz * dz);
      lmin = ceil( ((dist - delta) - pcf->t0) / dt );
      lmax = floor( ((dist + delta) - pcf->t0) / dt );
      if(lmax >= 0 && lmin < pcf->n) {
	/* kernel centred at 'dist' has nonempty intersection 
	   with specified range of t values */
	/* compute intersection */
	if(lmin < 0)
	  lmin = 0;
	if(lmax >= pcf->n)
	  lmax = pcf->n - 1;
	/* compute (inverse) edge correction weight */
	vx = b->x1 - b->x0 - (dx > 0 ? dx : -dx);
	vy = b->y1 - b->y0 - (dy > 0 ? dy : -dy);
	vz = b->z1 - b->z0 - (dz > 0 ? dz : -dz);
	invweight = vx * vy * vz * FOURPI * dist * dist;
	if(invweight > 0.0) {
	  for(l = lmin; l < pcf->n; l++) {
	    tval = pcf->t0 + l * dt;
	    /* unnormalised Epanechnikov kernel with halfwidth delta */
	    frac = (dist - tval)/delta;
	    kernel = (1 - frac * frac);
	    if(kernel > 0) 	    
	      (pcf->num)[l] += kernel / invweight;
	  }
	}
      }
    }
  }

  /* constant factor in kernel */
  coef = 3.0/(4.0 * delta);
  /* multiplied by 2 because we only visited i < j pairs */
  twocoef = 2.0 * coef; 

  /* normalise kernel and compute ratio estimate */
  for(l = 0; l < pcf->n; l++) {
    (pcf->num)[l] *= twocoef;
    (pcf->f)[l] = ((pcf->denom)[l] > 0.0)?
      (pcf->num)[l] / (pcf->denom)[l]
      : 0.0;
  }

}

void
pcf3isot(p, n, b, pcf, delta)
     Point *p;
     int n;
     Box *b;
     Ftable *pcf;
     double delta;
{
  register int i, j, l, lmin, lmax;
  register double dx, dy, dz, dist;
  Point *ip, *jp;
  double dt, vol, denom, mass, tval;
  double coef, frac, kernel;

  double sphesfrac(), sphevol();
  Point vertex;
  Box   half;

  /* compute denominator & initialise numerator*/
  vol = (b->x1 - b->x0) * (b->y1 - b->y0) * (b->z1 - b->z0);
  denom = ((double) (n * n))/vol;

  for(l = 0; l < pcf->n; l++) {
    (pcf->denom)[l] = denom;
    (pcf->num)[l]   = 0.0;
  }

  /* spacing of argument in result vector */
  dt = (pcf->t1 - pcf->t0)/(pcf->n - 1);

  /* set up for volume correction */

  vertex.x = b->x0;
  vertex.y = b->y0;
  vertex.z = b->z0;
  half.x1  = b->x1;
  half.y1  = b->y1;
  half.z1  = b->z1;
  half.x0  = (b->x0 + b->x1)/2.0;
  half.y0  = (b->y0 + b->y1)/2.0;
  half.z0  = (b->z0 + b->z1)/2.0;

	/* compute numerator */
  for( i = 0; i < n; i++) {
    ip = p + i;
    for(j = i + 1; j < n; j++) {
      jp = p + j;
      dx = jp->x - ip->x;
      dy = jp->y - ip->y;
      dz = jp->z - ip->z;
      dist = sqrt(dx * dx + dy * dy + dz * dz);
      lmin = ceil( ((dist - delta) - pcf->t0) / dt );
      lmax = floor( ((dist + delta) - pcf->t0) / dt );
      if(lmax >= 0 && lmin < pcf->n) {
	/* kernel centred at 'dist' has nonempty intersection 
	   with specified range of t values */
	/* compute intersection */
	if(lmin < 0)
	  lmin = 0;
	if(lmax >= pcf->n)
	  lmax = pcf->n - 1;
	/* compute edge correction weight */
	mass = (1.0 / sphesfrac(ip, b, dist)) 
	  + (1.0 / sphesfrac(jp, b, dist)); 
	mass *= 
	  1.0 - 8.0 * sphevol(&vertex, &half, dist) / vol;
	if(mass > 0.0) {
	  mass /= FOURPI * dist * dist;
	  for(l = lmin; l < pcf->n; l++) {
	    tval = pcf->t0 + l * dt;
	    /* unnormalised Epanechnikov kernel with halfwidth delta */
	    frac = (dist - tval)/delta;
	    kernel = (1 - frac * frac);
	    if(kernel > 0) 	    
	      (pcf->num)[l] += kernel * mass;
	  }
	}
      }
    }
  }

  /* constant factor in kernel */
  coef = 3.0/(4.0 * delta);

  /* normalise kernel and compute ratio estimate */
  for(l = 0; l < pcf->n; l++) {
    (pcf->num)[l] *= coef;
    (pcf->f)[l] = ((pcf->denom)[l] > 0.0)?
      (pcf->num)[l] / (pcf->denom)[l]
      : 0.0;
  }
}