/*  Solves Laplaces equation in a square Cartesian plane */
/* bounded by conducting walls. Compile with;             */                   /* gcc -O4 -o Laplace 2dlaplace.c -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM   -lICE -lXext -lX11 -lm */
/* The output reguires GNU plotutils-2.2 or above                   */
/* J. R. Schmidt 1999                                               */
/* Voltage illustrated with coloration and height; plots z=V(x,y)   */
/* Entire lattice is massivley vectorized to eliminate cache misses */
/* and to make stride considerations irrelevent.                    */
/* Point V(x,y) is translated to V[x+N*y]                           */
/* keep this in mind when you impose B.C.s                          */

#include<stdio.h>
#include<strings.h>
#include<stdlib.h>
#include<math.h>
#include<plot.h>

#define N 61       /* should be odd , this is grid size   */
/* creates a N X N grid, with middle being x-y axes       */
/* rotates figure by two Euler angles for different       */
/* camera angles                                          */
 
#define PI 3.1415927

double V[N*N],Vnew[N*N],X[N*N],Y[N*N],Xn[N*N],Yn[N*N];
double theta,phi,tmpx,tmpy,tmpz,dx,dy,dt,x,y,dimensionx, dimensiony,t,V0;
double ct,st,sp,cp;
double axes[4][2];
int mod(int A, int B);

int n,m,k,l;
int pl_handle,timestep,num,COLOR,displ;
FILE *fptr;

main (int argc, char *argv[])
{
  V0=6.0;  /* sets maximum voltage scale */
  dx = 0.2;
  dy = 0.2;
  dt = 0.05;
  dimensionx=(double)N*dx;
  dimensiony=(double)N*dy;
  /* Do some error checking to avoid segmentation faults */
 /* I hate segmentation faults.                          */

  if(argc < 5 || argc > 5){
    printf("usage: Laplace theta phi timesteps display\n");
    printf("Display=0 for X, 1 for GIF\n");
    exit(1);}

  theta=(double)atof(argv[1]);  /* azimuthal rotation angle */
  phi=(double)atof(argv[2]);    /* xy-plane rotation angle  */
  num=atoi(argv[3]);            /* number of time steps in relaxation */
  displ=atoi(argv[4]);
  cp=cos(phi);
  sp=sin(phi);
  ct=cos(theta);
  st=sin(theta);

  if(displ==0)
  pl_handle = pl_newpl ("X", stdin, stdout, stderr);
  else{
      fptr=fopen("laplace.gif","w+");
      pl_parampl ("BITMAPSIZE", "500x500");
      pl_parampl ("BG_COLOR", "white");
      pl_parampl ("TRANSPARENT_COLOR", "white");
      pl_handle = pl_newpl ("GIF", stdin, fptr, stderr);}
  

     
  pl_selectpl (pl_handle);
  if (pl_openpl () < 0)         /* open Plotter */
    {
      fprintf (stderr, "Couldn't open Plotter\n");
      return 1;
    }
  /* set scales for the size of the plotting window */
  pl_fspace (-dimensionx, -dimensiony,dimensionx, dimensiony);
  pl_flinewidth (0.01);

  /* this is why N must be odd */
 axes[0][0]=-(double)(N+1)*dx/2; /* points to draw x,y axes */
  axes[0][1]=0.0;
  axes[1][0]=(double)(N+1)*dx/2;
  axes[1][1]=0.0;
  axes[2][1]=-(double)(N+1)*dy/2;
  axes[2][0]=0.0;
  axes[3][1]=(double)(N+1)*dy/2;
  axes[3][0]=0.0;

  for(n=0;n<N*N;n++){
    X[n]=(double)(mod(n,N))*dx;
    Y[n]=(double)(mod((n-X[n])/N,N))*dy;
    /* now apply the offset */
    X[n]=X[n]-(double)(N-1)*dx/2.0;
    Y[n]=Y[n]-(double)(N-1)*dy/2.0;
    V[n]=0.0;
    Vnew[n]=0.0;
}
  /* this section sets up the boundary point voltages */
  for(n=0;n<N;n++){
    V[n]=0.0;  
    V[n+N*(N-1)]=V0; 
    Vnew[n]=0.0;       /* grounded edge , these are not updated*/
    Vnew[n+N*(N-1)]=V0; }




  /* begin the timesteps */
 for(timestep=1;timestep<=num;timestep++){
/* now update the lattice */
 for(k=1;k<N-1;k++){
   for(l=1;l<N-1;l++){
     /* update only interior points , edges are fixed */
     n=N*k+l; /* l*dx is the x-coord, k*dy is the y-coord */
     Vnew[n]=0.25*(V[n+1]+V[n-1]+V[n+N]+V[n-N]);
}}


  /* reset the lattice for the next sweep */
for(n=0;n<N*N;n++)
V[n]=Vnew[n];
 

 }/* end timesteps */ 
 /* graphics subsection */
 
/* now rotate the plot using phi and theta, to vantage point of choice*/
 /* the voltage at x,y becomes z-value of surface; (x,y, V(x,y))      */
  for(n=0;n<N*N;n++){
tmpx=cp*X[n]+sp*Y[n];
tmpy=-ct*sp*X[n]+ct*cp*Y[n]+st*V[n];
Xn[n]=tmpx;
Yn[n]=tmpy; }

  /* rotate axis lines too */
  for(n=0;n<4;n++){
tmpx=cp*axes[n][0]+sp*axes[n][1];
tmpy=-ct*sp*axes[n][0]+ct*cp*axes[n][1];
axes[n][0]=tmpx;
 axes[n][1]=tmpy;}

  /* we now create the graphics object */
pl_erase();
pl_pencolorname ("black");
pl_filltype(1);

 for(k=0;k<N-1;k++){
   for(l=0;l<N-1;l++){
n=N*N-2*N+l-N*k;
COLOR=(int)(65535-65535.0*V[n]/V0);
pl_fillcolor(COLOR,COLOR,COLOR);  /* grayscale coloration scheme */
/* pl_pencolor(COLOR,COLOR,COLOR);  */
pl_fmove(Xn[n],Yn[n]);             /* begin path at base point   */
pl_fline(Xn[n],Yn[n],Xn[n+1],Yn[n+1]);
pl_fline(Xn[n+1],Yn[n+1],Xn[n+1+N],Yn[n+1+N]);
pl_fline(Xn[n+1+N],Yn[n+1+N],Xn[n+N],Yn[n+N]);
pl_fline(Xn[n+N],Yn[n+N],Xn[n],Yn[n]);
pl_endpath();                     /* end closed 4-point path for filling */
;}}

  /* draw the axes in red */
 pl_pencolorname ("red");
pl_fline(axes[0][0],axes[0][1],axes[1][0],axes[1][1]);
pl_fline(axes[2][0],axes[2][1],axes[3][0],axes[3][1]);
  /* put in some labels */
pl_fmove(axes[1][0],axes[1][1]);
pl_ffontname("HersheySerif");
pl_ffontsize(0.6);
pl_label("x");
pl_fmove(axes[3][0],axes[3][1]);
pl_label("y"); 
pl_fmove(-0.75*dimensionx,0.5*dimensionx);
pl_label("\\gr\\sp2\\epV(x,y)=0");
/* finished; delete plotters */
        
  pl_closepl ();
  pl_selectpl (0);
  pl_deletepl (pl_handle);
  if(displ !=0) fclose(fptr);

}


 int mod(int A, int B){
  if(A>=0)
return(A-B*(A/B));
  else{
       do{A=A+B;}while(A<0);
       return(A);}
  }