/* relax.c
Simple-minded solution of Laplace's equation using relaxation method, with
overrelaxation w. The program is simple rather than elegant, and has no
error checking at all.  It should work with nearly any C compiler, and
probably even with most C++ compilers.

Output is written into a file called TEMP.$$$, with no formatting
except for a carriage return at the end of each row.  To write
output to console, just replace the "fprintf(wfile,..." statements
with "printf(...".   Some additional formatting would be nice if the
output is to be read by humans, rather than spreadsheets or a graphics program.

Last modified 10/06/09 by E**3.
*/

#include <stdio.h>
#include <math.h>
#define VSIZE 20
#define HSIZE 20

float phi[VSIZE][HSIZE];    /* Could be double if more precision needed */
main()
{
double vtop,vbot,vleft,vright,maxdiff,newval,w,tol;
int i,j,last_j,iter;
FILE *wfile;

printf("Enter V(top): ");
scanf("%lf",&vtop);
printf("Enter V(bottom): ");
scanf("%lf",&vbot);
printf("Enter V(left): ");
scanf("%lf",&vleft);
printf("Enter V(right): ");
scanf("%lf",&vright);
printf("Enter tolerance: ");
scanf("%lf",&tol);
printf("Enter overrelaxation (1-2): ");
scanf("%lf",&w);

/* Set the entire square array to zero initially */
for (i=0; i<VSIZE; i++)
  for (j=0; j<HSIZE; j++)
    phi[i][j]=0.0;

/* Now enter the boundary conditions. */

for (j=0;j<HSIZE;j++)
  {
  phi[0][j]=vtop;
  phi[VSIZE-1][j]=vbot;
  }
for (i=0;i<VSIZE;i++) 
  {
  phi[i][0]=vleft;
  phi[i][HSIZE-1]=vright;   // For an angled boundary, need to modify this.
  }

/*
Do relaxation to a desired level of convergence set by "tol."  To get
started, set maxdiff to a value larger than tol.
*/
maxdiff = 2*tol;

for (iter=1; maxdiff>tol; iter++) {

  /* Do an iteration.  Be careful not to overwrite the boundaries! */
  maxdiff=0;
  for (i=1; i<VSIZE-1; i++) {
    last_j = HSIZE-1;     // Modify to achieve an angled boundary. */
    for (j=1; j<last_j; j++) {
      newval = (phi[i][j+1] + phi[i+1][j] + phi[i][j-1] + phi[i-1][j]) * w/4
		+ (1.0-w)*phi[i][j];
      if (fabs(newval-phi[i][j]) > maxdiff)  // Adjust maximum difference if needed
		maxdiff = fabs(newval-phi[i][j]);
      phi[i][j] = newval;
      }
    }
}

/*Print the results */
printf("Maximum difference is %le\n",maxdiff);
printf("Number of iterations is %d\n",iter);
printf("Printing output to file \"temp.$$$\"\n");
wfile = fopen("temp.$$$","w");
for (i=0; i<VSIZE; i++)
  {
  for (j=0; j<HSIZE; j++) {
    fprintf (wfile,"%6.4f,",phi[i][j]);
    }
  fprintf(wfile,"\n",i);
    }
fclose(wfile);

printf("Done!  Type 'Enter' to quit: ");
getchar();  getchar();

return(0);

}