;PRO demo_laplace
;Show smoothing for small and large scale
;perturbations
;Aim: Solve Poisson equation in 2D
;-Laplace phi=rho

;device,decompose=0
;Setup the grid
Lx=10.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)
nx=40 &ny=nx
h=Lx/(nx-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=5.0 & y0=5.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=0.0*exp(-((xv-x0)^2+(yv-y0)^2)/l) $
-0.0*exp(-((xv-x1)^2+(yv-y1)^2)/l)
phi_new=phi
;
k=1
print,'Input wave number of disturbance:'
read, k
hu=1.0*sin(k*!pi*xv/Lx)*sin(k*!pi*yv/Ly)
phi=phi+hu
;
method=1
print,' '
print,'Which method to use?'
print,'1=Jacobi, 2=Gauss-Seidel, 3=SOR'
read,method
tol=1.0e-6; Tolerance level for error
w=1.0 ; for overrelextion: 1.0<w<2.0
if method eq 3 then begin
print,'Input SOR-relaxation parameter 0<w<2'
read,w
endif




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Solve Laplace equation                            ;;
;Only slight modifications between Jacobi, Gauss-Seidel and
;Succesive over-relaxation methods;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
test=1000.0; initial value for error
it=long(0)
while (1 eq 1) and (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
if (method eq 1) then $
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
;
if (method eq 2) or (method eq 3) then $
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_new[ix-1,iy]+phi_new[ix,iy-1])) ;Gauss-Seidel/SOR method
endfor & endfor
surface,phi,charsize=2.0,zrange=[-1,1],xtitle='x',$
ytitle='y',ztitle='Phi(x,y)',$
title='After'+strcompress(string(it))+' Iterations'
wait,0.2
;
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

;

;
dummy=' '
print,'Press enter to continue'
read,dummy
plot,it_vec,error_vec,/xlog,/ylog,xtitle='Iteration steps',$
title='Error= Total(Phi_new-Phi)^2)',charsize=1.5

END