(* ::Package:: *)

(*initialisation.m*)

sd=10^6; tol=10^(-8)(*default = 10^(-8)*); nt=1; epsilon0=1.0;rhov=1.0*epsilon0;
boxsize=2.0;(*default =2.0*)


(*Provide N, the "resolution" of the system we wish to simulate*)
NStep=20;
iini=0;ilast=NStep;
jini=0;jlast=NStep;
kini=0;klast=NStep;

iprime=IntegerPart[(ilast-iini)/2.0];
jprime=IntegerPart[(jlast-jini)/2.0];
kprime=IntegerPart[(klast-kini)/2.0];


(*End of provide N, the "resolution" of the system we wish to simulate*)

(*The coordinates of the boundaries*)
x[iini]=-boxsize;x[ilast]=boxsize;
y[jini]=-boxsize;y[jlast]=boxsize;
z[kini]=-boxsize;z[klast]=boxsize;
Deltax=(x[ilast]-x[iini])/NStep;(*Step size*)
(*End of the coordinates of the boundaries*)


(*The boundary values of V at the boundaries*)
vx[iini]=0.0;vx[ilast]=0.0;
vy[jini]=0.0;vy[jlast]=0.0;
vz[kini]=0.0;vz[klast]=0.0;
(*End of the boundary values of V at the boundaries*)

(*Defining the values of coordinates of x[i],y[j],z[k]*)
For[j=jini,j<=jlast,j++,
	y[j]=y[jini]+(j-jini)*Deltax;
];(*end for j*)

For[i=iini,i<=ilast,i++,
	x[i]=x[iini]+(i-iini)*Deltax;
];(*end for j*)

For[k=kini,k<=klast,k++,
	z[k]=z[kini]+(k-kini)*Deltax;
];(*end for k*)

(*End of defining the values of coordinates of x[i],y[j],z[k]*)

(*Emptying all values in V[i,j,k]*)

For[k=kini,k<=klast,k++,
For[j=jini,j<=jlast,j++,
For[i=iini,i<=jlast,i++,
rho[i,j,k]=0;
V[i,j,k]=0;
Vold[i,j,k]=0;
];(*end for i*)
];(*end for j*)
];(*end for k*)
(*End of emptying all values in V[i,j,k]*)

rho[iprime,jprime,kprime]=rhov;


(*Overwritting V at the boundaries with the boundary conditions*)
For[k=kini,k<=klast,k++,
For[j=jini,j<=jlast,j++,
For[i=iini,i<=ilast,i++,

If[i==iini,V[i,j,k]=vx[iini]];(*end if i=iini*)
If[i==ilast,V[i,j,k]=vx[ilast]];(*end if i=ilast*)

If[j==jini,V[i,j,k]=vy[jini]];(*end if i=iini*)
If[j==jlast,V[i,j,k]=vy[jlast]];(*end if i=ilast*)

If[k==kini,V[i,j,k]=vz[kini]];(*end if i=iini*)
If[k==klast,V[i,j,k]=vz[klast]];(*end if i=ilast*)

];(*end for i*)
];(*end for j*)
];(*end for k*)
(*Overwritting V at the boundaries with the boundary conditions*)



tV=Table[
{x[i],y[j],z[k],V[i,j,k]},{i,iini,ilast},{j,jini,jlast},{k,kini,klast}
];


tV1=Flatten[tV,1];


tVz=Table[{x[i],y[j],V[i,j,0]},{i,iini,ilast},{j,jini,jlast}];

tVz1=Flatten[tVz,1];

Print["Print from initialisation.m"];
Print[ListContourPlot[tVz1,Mesh->None,ContourLabels->All]];


Deltax
