/* Successive Over-Relaxation solver */

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

extern NPTS;

fpt
relax(double alpha, fpt *field, int *bound) {
  int i, j, k;
  fpt tmp, change, new, foo;
  change = 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) == -1) {
		tmp = F(field,i,j,k);
		new =
		  (((1.0-alpha)*
		    ((fpt) F(field, i, j, k))));
		foo = ((fpt) F(field,i,j-1,k));
		foo+= ((fpt) F(field,i-1,j,k));
		foo+= ((fpt) F(field,i+1,j,k));
		foo+=((fpt) F(field,i,j+1,k));
		foo+=((fpt) F(field,i,j,k-1));
		foo+=((fpt) F(field,i,j,k+1));
		foo/=6;
		foo*=alpha;
		new += foo;
		(F(field,i,j,k)) = new;
		tmp = (F(field,i,j,k)-tmp);
		change += (tmp*tmp);
	      }
	      else {
			F(field, i, j, k) = (fpt) F(bound, i, j, k);
		}
	    }
	  }
	}
	
	for(i=0; i< NPTS; ++i)
	  for(j=0;j < NPTS; ++j){
	    F(field, i, j, 0) = F(field, i, j, 1);
	    F(field, i, 0, j) = F(field, i, 1, j);
	    F(field, 0, i, j) = F(field, 1, i, j);
	    F(field, i, j, NPTS-1) = F(field, i, j, NPTS-2);
	    F(field, i, NPTS-1, j) = F(field, i, NPTS-2, j);
	    F(field, NPTS-1, i, j) = F(field, NPTS-2, i, j);
	  }
	return change;
}

fpt
charge(fpt *field, int *bound, int objnum) {
  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) == objnum){
		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;
}