/* Solution to the Laplace Equation by relaxation                 */
/* compile line " gcc -o Electrode Electrode.c                    */
/* Useage; "Laplace volts filename" will create a file called     */
/* "filename" of points ("volts" must be a floating num)          */
/* on the equipotential V(x,y)=volts                              */
/* to set up electrode voltages, edit "restore_bc() routine       */
/* at the very end of the program.                                */

#include<stdio.h>
#include<stdlib.h>
#define N  100


float V[N + 1][N + 1];
float V_0[N + 1][N + 1];
float v0, testval,volt;
int count, n, m, k, p, m1, p1, m2, p2, m3, p3, X, Y;
FILE *fptr;
void restore_bc(float (*A)[N+1]); /* routine to restore electrode voltages */



int main (int argc, char *argv[])
{
  
 float (*Vptr)[N + 1], (*V0ptr)[N + 1];


  /* Do some error checking to avoid segmentation faults */
 /* I hate segmentation faults.                          */

  if(argc < 3 || argc > 3){
    printf("usage: Laplace voltage filename\n");
    exit(1);}

volt=(float)atof(argv[1]);  
fptr=fopen(argv[2],"w+");
v0=10.0;

  /* Set up grounded bounding box.    */
  for (n = 0; n <= N; n++)
    {
      V_0[n][0] = 0.0;
      V_0[0][n] = 0.0;
      V_0[N][n] = 0.0;
      V_0[n][N] = 0.0;
      V[n][0] = 0.0;
      V[0][n] = 0.0;
      V[N][n] = 0.0;
      V[n][N] = 0.0;
    }
  /* Set up the zeroth estimate of interior voltages. */
  for (m = 1; m < N; m++)
    {
      for (p = 1; p < N; p++)
	{
	  V_0[m][p] = 1.0;
	  V[m][p] = 0.0;
	}
    }


  /* Set up the array pointer. */
  Vptr = &V[0][0];
  V0ptr = &V_0[0][0];

  /* set the electrode voltages for the first time */
restore_bc(V0ptr);

  /* Begin a sequence of 4000 updates for the interior.       */
  /* If you have all day, you can replace this with a do loop */
  /* that computes until arbitrary precision is attained, or  */
  /* the cows come home, whichever comes first.               */

  for (count = 0; count < 4000; count++)
    {


      for (m1 = 1; m1 < N; m1++)
	{
	  for (p1 = 1; p1 < N; p1++)
	    {
	      *(*(Vptr + m1) + p1) =
		0.25 * (*(*(V0ptr + m1 + 1) + p1) + *(*(V0ptr + m1 - 1) + p1) +
		*(*(V0ptr + m1) + p1 + 1) + *(*(V0ptr + m1) + p1 - 1));
	    }
	}


      for (m2 = 1; m2 < N; m2++)
	{
	  for (p2 = 1; p2 < N; p2++)
	    {
	      *(*(V0ptr + m2) + p2) = *(*(Vptr + m2) + p2);
	    }
	}
      /* restore the electrode voltages*/

 restore_bc(V0ptr);


    }
  
  /* Now we can dump out equipotentials between 0 < volt < 10 V  */
  /* within a tolerance of sqrt(0.005) volts.                    */
  
  fprintf(fptr,"# Points on the %f volt equipotential.\n",volt);
  for (m3 = 1; m3 < N; m3++)
    {
      for (p3 = 1; p3 < N; p3++)
	{
	  testval = *(*(Vptr + m3) + p3);
	  if ((testval - volt) * (testval - volt) < 0.005)
	    fprintf (fptr,"%d\t%d\n", m3, p3);

	}
    }

  fclose(fptr);
  return (0);
}

void restore_bc( float (*A)[N+1])
{

  /* this is where you specify the voltages V[n][m] on electrodes */
int n,m;

/* for example, a strip 40 units width, five deep, at end of rod*/
 for(n=1;n<30;n++){
   *(*(A+49)+n)= 10.0;
   *(*(A+50)+n)= 10.0;
   *(*(A+51)+n)= 10.0;
}

 for(n=30;n<70;n++){
*(*(A+n)+31)= 10.0;
*(*(A+n)+32)= 10.0;
*(*(A+n)+33)= 10.0;
*(*(A+n)+34)= 10.0;
 *(*(A+n)+35)= 10.0;}

}