(* ::Package:: *)

(*update.m*)

(*In this subroutine, Jacobi relaxation will be executed to update the potential fields until 
convergence requirement (determined by tol) is achieved. The output would in the List of Table 
named "tV02" containing the data of the format {x[i],y[j],V[i,j]}*)


While[sd>tol,
tV={};tt={};
tVz={};ttz={};

For [k=kini+1,k<=klast-1,k++,(*end for k*)
For [j=jini+1,j<=jlast-1,j++,(*end for j*)
For [i=iini+1,i<=ilast-1,i++,(*end for i*)

Vold[i,j,k]=V[i,j,k];
V[i,j,k]=(1/6.0)*(V[i-1,j,k]+V[i+1,j,k]+V[i,j+1,k]+V[i,j-1,k]+V[i,j,k+1]+V[i,j,k-1])+ Deltax^2*rho[i,j,k]/(6.0*epsilon0);
DeltaV[i,j,k]=Abs[V[i,j,k]-Vold[i,j,k]];
](*end for i*)
](*end for j*)
](*end for k*)

(*Creating a list the form of {x[i],y[j],V[i,j]} for contourplot *)

For[k=kini,k<=klast,k++,
	For[j=jini,j<=jlast,j++,
		For[i=iini,i<=ilast,i++,
		tt={x[i],y[j],z[k],V[i,j,k]};
		tV=Append[tV,tt];

		ttz={x[i],y[j],V[i,j,kprime]};
		tVz=Append[tVz,ttz];
		];
	];
];
(*tV2=Table[{x[i],y[j],V[i,j]},{i,iini,ilast},{j,jini,jlast}];*)


Print["Print from update.m"];
Print["nt, sd = ",nt, " ",sd];
Print[ListContourPlot[tVz,Mesh->None]];
Print[" "];

(*Test convergence*)
sd=
Sqrt[
Sum[
DeltaV[i,j,k]^2,{i,iini+1,ilast-1},{j,jini+1,jlast-1},{k,kini+1,klast-1}
]
]/
(
(ilast-iini-1)*(jlast-jini-1)*(klast-kini-1)
);


If[sd<tol,
Print["Convergence achieved, sd= ",sd]
];(*end if test tol*)
(*End of test convergence*)

nt=nt+1;

];(*end while*)



