#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Using Cactus infrastructure */
#include "cctk.h" 

/* Using Cactus parameters */
#include "cctk_Parameters.h"

/* Using Cactus arguments lists */
#include "cctk_Arguments.h"

int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d);
int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d);

void WaveBinaryC(CCTK_ARGUMENTS)
{
  
  /* Declare variables in argument list */
  DECLARE_CCTK_ARGUMENTS
    
  /* Declare parameters */
  DECLARE_CCTK_PARAMETERS

  int i,j,k,d,index;
  
  /* the start/end points in local index space of the binary charge */
  int lowerloc[3];
  int upperloc[3];
  
  /* lower/upper physical coord. of one binary charge */
  CCTK_REAL lbin[3],ubin[3];
  CCTK_REAL minc[3],maxc[3];
  
  /*  lower/upper grid points, that we own (no ghostzones) */ 
  int mingp[3], maxgp[3];

  /* some flags */
  int onthisproc, xproc, sign;
  static int firstcall;
  
  CCTK_REAL xs,ys,zs,rad;
  CCTK_REAL zmin,zmax,charge_factor;

  /* Because binary_charge and binary_size are a steerable parameters now,
     charge_factor needs to be recomputed at every iteration. */
  charge_factor = 3.0*binary_charge/(4.0*3.1415*binary_size*binary_size*binary_size);

  /* Initialize the range arrays */
  for (d=0;d<3;d++)
  {
    mingp[d] = 0;
    maxgp[d] = cctk_lsh[d]-1;
  } 


  /* The physical minimum/maximum coordinates of our grid chunk;
     will use those to check if the charges are on our proc. */
  index   = CCTK_GFINDEX3D(cctkGH,0,0,0);
  minc[0] = x[index];
  minc[1] = y[index];
  minc[2] = z[index];
  index   = CCTK_GFINDEX3D(cctkGH,cctk_lsh[0]-1,cctk_lsh[1]-1,cctk_lsh[2]-1);
  maxc[0] = x[index];
  maxc[1] = y[index];
  maxc[2] = z[index];
        

  /*  we have two charges opposite, origin is the center
      sign mulitplies the charge xy-positions by +1 and -1 */
  for (sign=1; sign>=-1; sign-=2)
  {

    /* default flag: the charge is on this proc. */
    onthisproc = 0;

    /* calculate the center of a single binary source */
    CCTK_CoordRange(cctkGH, &zmin, &zmax, -1, "z", "cart3d");
    xs = sign * binary_radius * cos(binary_omega*cctk_time);
    ys = sign * binary_radius * sin(binary_omega*cctk_time);
    zs = (zmax-zmin)/2 + zmin;

    /* shortcuts for the physical boundaries of 
       a single charge */
    lbin[0] = xs - binary_size;
    lbin[1] = ys - binary_size;
    lbin[2] = zs - binary_size;

    ubin[0] = xs + binary_size;
    ubin[1] = ys + binary_size;
    ubin[2] = zs + binary_size;

   
    /* loop over all three dimensions */
    for (d=0;d<3;d++)
    {   
     
      /* check where our binaries are:
	 - the can be on our grid chunk completly
	 - we have a lower part of the source on our grid chunk
	 - we have a upper part
	 - the source covers a grid chunk completetly 
	 
	 Get the floor-index of the lower bound and 
	 ceiling-index of the upper bound.

	 Increase onthisproc : if equal to three, we have the source
      */
      
      xproc = 0;
      lowerloc[d]=mingp[d];
      upperloc[d]=maxgp[d];

      if ((lbin[d]<=minc[d]) && (ubin[d]>=maxc[d]))
      {    
	xproc = 1;
      } 
      else 
      {
	if ((minc[d]<=lbin[d]) && (lbin[d]<=maxc[d]))
	{   
	  lowerloc[d]=IndexFloorC(cctkGH, lbin[d],d);
	  xproc = 1;
	}
	if ((minc[d]<=ubin[d]) && (ubin[d]<=maxc[d]))
	{
	  upperloc[d]=IndexCeilC (cctkGH, ubin[d],d);
	  xproc = 1;
	}	
      }
      onthisproc+=xproc;
    }

    /* Debugging debugging */
    if (((firstcall==1)&&(CCTK_EQUALS(binary_verbose,"yes"))) ||
	(CCTK_EQUALS(binary_verbose,"debug")))
    {
      printf("Charge: %f \n",charge_factor);
      printf("On this proc: %d \n",onthisproc);
      printf("cctk_lsh: %d %d %d\n",cctk_lsh[0],cctk_lsh[0],cctk_lsh[0]);
      printf("Charge center: %f %f %f \n", xs,ys,zs);
      printf("Charge extension: \n");
      printf("     x-extension: %f -> %f \n",lbin[0],ubin[0]);
      printf("     y-extension: %f -> %f \n",lbin[1],ubin[1]);
      printf("     z-extension: %f -> %f \n",lbin[2],ubin[2]);
      printf("Charge local index range: \n");
      printf("     x-index   %d %d \n", lowerloc[0],upperloc[0]);
      printf("     y-index   %d %d \n", lowerloc[1],upperloc[1]);
      printf("     z-index   %d %d \n", lowerloc[2],upperloc[2]);
      index = CCTK_GFINDEX3D(cctkGH, lowerloc[0],lowerloc[1],lowerloc[2]);
      printf("Charge corner  LOW: [%f %f %f] \n",x[index],y[index],z[index]);
      index = CCTK_GFINDEX3D(cctkGH, upperloc[0],upperloc[1],upperloc[2]);
      printf("Charge corner  UP : [%f %f %f] \n",x[index],y[index],z[index]);
      printf(" MINC/MAXC  [%f %f %f] [%f %f %f] \n\n\n",
	     minc[0],minc[1],minc[2],
	     maxc[0],maxc[1],maxc[2]);
    }
      
    /* we have the source, no loop over the index boudary */
    if (onthisproc==3)
    {
      for (k=lowerloc[2];k<=upperloc[2];k++)
	for (j=lowerloc[1];j<=upperloc[1];j++)
	  for (i=lowerloc[0];i<=upperloc[0];i++)
	  {
	    index =  CCTK_GFINDEX3D(cctkGH,i,j,k);
	    rad = 
	      ((x[index]-xs)*(x[index]-xs)) +
	      ((y[index]-ys)*(y[index]-ys)) +
	      ((z[index]-zs)*(z[index]-zs));

	    if (rad<binary_size*binary_size) 
	      phi_n[index] = phi_n[index] + charge_factor;
	  }
    }

    /* end of the sign loop */
  }


  /* reset firstcall to zero */
  firstcall = 0;
  
  /*  Note, that we do not need to sync anything, since each grid 
      patch has filled out its ghostzones. */

  return;
}