; PRO lecture_poisson_task1
;solves poisson-equation in 2D
;with Jacobi, Gauss-Seidel and SOR-method
dummy=' '& print,'' &print,'Press Enter to continue' & read,dummy
test=1000.0; initial value for error
it=0
while (test gt tol) and (it lt 10000) do begin
it=it+1
for ix=1,nx-2 do begin
for iy=1,ny-2 do begin
phi_new[ix,iy]=(1.0-w)*phi[ix,iy]+w* $
(0.25*(h^2*rho[ix,iy]+phi[ix+1,iy]+phi[ix,iy+1]+$
phi[ix-1,iy]+phi[ix,iy-1])) ;Jacobi method
;phi_new[ix-1,iy]+phi_new[ix,iy-1])) ;Gauss-Seidel method
endfor & endfor
test=total((phi_new-phi)^2);/nx/ny
print,'Step ',it, ' ,Error=',test
if it eq 1 then error_vec=[test] else $
error_vec=[error_vec,test]
phi=phi_new
endwhile
it_vec=findgen(N_ELEMENTS(error_vec))+1
contour,phi,x,y,/fill,nlevels=20,$
min_value=-10,max_value=10
;
;compute gradients for E=-grad phi
Ex=fltarr(nx,ny)
Ey=fltarr(nx,ny)
Ex=(shift(phi,1,0)-shift(phi,-1,0))/(2.0*h)
Ey=(shift(phi,0,1)-shift(phi,0,-1))/(2.0*h)
;vel,ex,ey,title='E-field'

plot,it_vec,error_vec,/xlog,/ylog,xtitle='Iteration steps',$
title='Error= Total(Phi_new-Phi)^2)',charsize=1.5
dummy=' '& print,' ' &print,'Press Enter to continue' & read,dummy
contour,phi,x,y,nlevels=20,charsize=1.5, $
xtitle='x',ytitle='y',title='Electric potential and electric field'
velovect,ex,ey,x,y,/overplot,length=2