;PRO lecture_poisson2D_draft
;Aim: Solve Poisson equation in 2D
;-Laplace phi=rho

;device,decompose=0
;Setup the grid
Lx=12.0
Ly=10.0
h=0.4; grid resolution; try also h=0.2 and h=0.1
nx=long(Lx/h+1)
ny=long(Ly/h+1)
print,'nx=',nx,' ,ny=',ny
x=Lx*findgen(nx)/(nx-1)
y=Ly*findgen(ny)/(ny-1)
phi=fltarr(nx,ny)

;
;Some initialization
rho=fltarr(nx,ny)
xv=fltarr(nx,ny)
yv=fltarr(nx,ny)
for iy=0,ny-1 do xv[*,iy]=x
for ix=0,nx-1 do yv[ix,*]=y
;
;Sources rho(x,y)
x0=4.0 & y0=4.0
x1=8.0 & y1=4.0
l=1.0; How extended are the charges?
;Feel Free to add more or less electric charges
rho=5.0*exp(-((xv-x0)^2+(yv-y0)^2)/l) $
-5.0*exp(-((xv-x1)^2+(yv-y1)^2)/l)
phi_new=phi
contour,rho,x,y,/fill,nlevels=20,charsize=1.5,$
min_value=-10,max_value=10, xtitle='x',ytitle='y',$
title='Electric charge distribution rho(x,y)'



tol=1.0e-6; Tolerance level for error
w=1.0 ; for overrelextion: 1.0<w<2.0
      ;Try to find the optimal value for w
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Solve Laplace equation                            ;;
;Only slight modifications between Jacobi, Gauss-Seidel and
;Succesive over-relaxation methods are neccessary
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
@lecture_poisson_task1.pro
END