#include <math.h>

#include "cctk.h"
#include "cctk_Parameters.h"

void BndApplyRobin3Di(cGH *GH,
		      int  gdim,
		      int *sw, 
		      int *doBC, 
		      int *lsh, 
		      CCTK_REAL *var, 
		      CCTK_REAL *r, 
		      CCTK_REAL finf, 
		      CCTK_INT npow) 
{

  int xgp1,xgp2,xgp3,xgp4;
  int ygp1,ygp2,ygp3,ygp4;
  int zgp1,zgp2,zgp3,zgp4;

  int i,j,k;

  if ((sw[0]!=1) || (sw[1]!=1) || (sw[2]!=1))
    CCTK_WARN(0,"Robin BC only with stencil_size=1");

  /* Lower x-bound: x(2,:,:) --> x[xgp2] */
  if (doBC[0]==1) {
    for (k=0;k<lsh[2];k++) {
      for (j=0;j<lsh[1];j++) {
	
        xgp1 = CCTK_GFINDEX3D(GH,0,j,k);
	xgp2 = CCTK_GFINDEX3D(GH,1,j,k);
	xgp3 = CCTK_GFINDEX3D(GH,2,j,k);
	xgp4 = CCTK_GFINDEX3D(GH,3,j,k);

	var[xgp1] = finf + (
                pow(r[xgp4],npow)*( var[xgp4] - finf )
                -3.0*pow(r[xgp3],npow)*( var[xgp3] - finf )
                +3.0*pow(r[xgp2],npow)*( var[xgp2] - finf )
                )/pow(r[xgp1],npow);
      }
    }
  }

  /* Upper x-bound: x(nx,:,:) --> xgp[xgp0] */
  if (doBC[1]==1){
    for (k=0;k<lsh[2];k++) {
      for (j=0;j<lsh[1];j++) {	
	xgp1 = CCTK_GFINDEX3D(GH,lsh[0]-1,j,k);
        xgp2 = CCTK_GFINDEX3D(GH,lsh[0]-2,j,k);
        xgp3 = CCTK_GFINDEX3D(GH,lsh[0]-3,j,k);
	xgp4 = CCTK_GFINDEX3D(GH,lsh[0]-4,j,k);
	
	var[xgp1] = finf + (
             pow(r[xgp3],npow)*( var[xgp4] - finf )
             -3.0*pow(r[xgp3],npow)*( var[xgp3] - finf )
             +3.0*pow(r[xgp2],npow)*( var[xgp2] - finf )
             )/pow(r[xgp1],npow);
      }
    }
  }
  
  /* Lower y-bound */
  if (doBC[2] == 1) {
    for (k=0;k<lsh[2];k++) {
      for (i=0;i<lsh[0];i++) {
	ygp1 = CCTK_GFINDEX3D(GH,i,0,k);
	ygp2 = CCTK_GFINDEX3D(GH,i,1,k);
	ygp3 = CCTK_GFINDEX3D(GH,i,2,k);
	ygp4 = CCTK_GFINDEX3D(GH,i,3,k);
	
	var[ygp1] = finf + (
                pow(r[ygp4],npow)*( var[ygp4] - finf )
                -3.0*pow(r[ygp3],npow)*( var[ygp3] - finf )
                +3.0*pow(r[ygp2],npow)*( var[ygp2] - finf )
                )/pow(r[ygp1],npow);

      }
    }
  }

  /* Upper y bound */
  if (doBC[3] == 1) {
    for (k=0;k<lsh[2];k++) {
      for (i=0;i<lsh[0];i++) {
	ygp1 = CCTK_GFINDEX3D(GH,i,lsh[1]-1,k);
	ygp2 = CCTK_GFINDEX3D(GH,i,lsh[1]-2,k);
 	ygp3 = CCTK_GFINDEX3D(GH,i,lsh[1]-3,k);
	ygp4 = CCTK_GFINDEX3D(GH,i,lsh[1]-4,k);

   	var[ygp1] = finf + (
             pow(r[ygp3],npow)*( var[ygp4] - finf )
             -3.0*pow(r[ygp3],npow)*( var[ygp3] - finf )
             +3.0*pow(r[ygp2],npow)*( var[ygp2] - finf )
             )/pow(r[ygp1],npow);
	  }
    }
  }

  /* Lower z-bound */
  if (doBC[4]==1) {
    for (j=0;j<lsh[1];j++) {
      for (i=0;i<lsh[0];i++) {
	zgp1 = CCTK_GFINDEX3D(GH,i,j,1);
        zgp2 = CCTK_GFINDEX3D(GH,i,j,2);
        zgp3 = CCTK_GFINDEX3D(GH,i,j,3);
        zgp4 = CCTK_GFINDEX3D(GH,i,j,4);
    
	var[zgp1] = finf + (
                pow(r[zgp4],npow)*( var[zgp4] - finf )
                -3.0*pow(r[zgp3],npow)*( var[zgp3] - finf )
                +3.0*pow(r[zgp2],npow)*( var[zgp2] - finf )
                )/pow(r[zgp1],npow);
      }
    }
  }

  /* Upper z-bound */
  if (doBC[5] == 1) {
    for (j=0;j<lsh[1];j++) {
      for (i=0;i<lsh[0];i++) {
	zgp1 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-1);
	zgp2 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-2);
	zgp3 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-3);
	zgp4 = 	CCTK_GFINDEX3D(GH,i,j,lsh[2]-4);

     	var[zgp1] = finf + (
             pow(r[zgp3],npow)*( var[zgp4] - finf )
             -3.0*pow(r[zgp3],npow)*( var[zgp3] - finf )
             +3.0*pow(r[zgp2],npow)*( var[zgp2] - finf )
             )/pow(r[zgp1],npow);

      }
    }
  }
}