/* Actually do the math */

#include "main.h"
#include "xgraphics.h"
#include "sor.h"


fpt
relax(double alpha, fpt field[NPTS][NPTS][NPTS], int bound[NPTS][NPTS][NPTS]) {
	int i, j, k;
	fpt tmp, change;

	fprintf(stderr, "Slice ");
	fflush(stderr);
	for (i = 1; i < NPTS-1; ++i) {
	  fprintf(stderr, "%d, ", i);
	  fflush(stderr);
	  for (j = 1; j < NPTS-1; ++j) {
	    for (k = 1; k < NPTS-1; ++k) {
	      if(F(bound,i,j,k) == -1) {
		/*		if(F(field, i,j,k) > .01)
				fprintf(stderr, "F(%d, %d, %d) = %f\n", i, j, k, F(field,i,j,k));
				*/
		tmp = F(field,i,j,k);
		F(field,i,j,k) = (((1.0-alpha)*((fpt) F(field, i, j, k))) +
				  alpha*
				  (
				   (
				    ((fpt) F(field,i-1,j,k))+
				    ((fpt) F(field,i,j-1,k))+
				    ((fpt) F(field,i+1,j,k))+
				    ((fpt) F(field,i,j+1,k))+
				    ((fpt) F(field,i,j,k-1))+
				    ((fpt) F(field,i,j,k+1))
				    )/
				   6.0)
				  );
		tmp = (F(field,i,j,k)-tmp);
		change += (tmp*tmp);
	      }
	      else
		F(field, i, j, k) = (fpt) F(bound, i, j, k);
	    }
	  }
	}
      return change;
}

fpt
charge(fpt field[NPTS][NPTS][NPTS], int bound[NPTS][NPTS][NPTS]) {
  int i,j,k;
  fpt chg=0.0;

    for (i = 1; i < NPTS-1; ++i)
         for (j = 1; j < NPTS-1; ++j)
	    for (k = 1; k < NPTS-1; ++k) {
	      if(F(bound,i,j,k) >=0){
		chg+= (F(field,i,j,k) - F(field,i,j,k+1));
		chg+= (F(field,i,j,k) - F(field,i,j,k-1));
		chg+= (F(field,i,j,k) - F(field,i+1,j,k));
		chg+= (F(field,i,j,k) - F(field,i-1,j,k));
		chg+= (F(field,i,j,k) - F(field,i,j-1,k));
		chg+= (F(field,i,j,k) - F(field,i,j+1,k));
	      }
	    }
    return chg;
}