pro poissond_qd, f, u, ds=ds, verbose=verbose, eps=eps, maxits=maxits
 
; PURPOSE: This code was shamelessly ripped off from Dietmar
; Krauss-Varban & George Fisher; Dietmar's (FORTRAN) version was
; inspired by the SOR routine from numerical recipes 2nd version, but
; written completely from scratch, and George modified it.
;
; It solves Poisson's equation with source f for u, with Dirichlet
; boundary conditions on the outer two zones, which are kept at zero.
;
; Hence, Poisson's equation isn't solved on the outer two zones.
;
if (n_params() eq 0) then begin
    message,' NAME: POISSOND_QD.PRO',/info
    message,'',/info
    message,' PURPOSE: This quick & dirty Poisson solver finds 2D ',/info
    message,'    scalar potential "u" whose horizontal Laplacian equals',/info
    message,'    2D source array "f".  If input is double precision, ',/info
    message,'    then output is, too.',/info
    message,'',/info
    message,' USAGE: poissond_qd, f, u, ds=725., /verbose, eps=1e-5',/info
    message,'',/info
    message,' REQUIRED: f = 2D-array , either float or double.', /info
    message,'',/info
    message,' OPTIONAL KEYWORDS: DS = pixel size; default = 1. (float)',/info
    message,'    VERBOSE = set to print diagnostic info.',/info
    message,'    EPS = tolerance on integrated error as fraction of',/info
    message,'          integrated abs. source function; default = 1e-5',/info
    message,'    MAXITS = max. # of iterations; default = 10000L ',/info
    message,'',/info
    message,' ERROR HANDLING: Primitive: no input checks, errors cause halt.',/info
    message,'',/info
    message,' HISTORY: Last modified 15-Jan-2008.'/info
    return
endif 

dims = size(f)
imax = dims(1)	; define dimensionality
jmax = dims(2)

; Define arrays, either double or single precision, based upon input
;===================================================================
if (n_elements(u) eq 0) then $
  if (dims(3) eq 5) then u = dblarr(imax,jmax) else u = fltarr(imax,jmax) 

if not(keyword_set(ds)) then begin
    ds = 1.
    message,'Assuming unit pixel scale.',/info
endif

if not(keyword_set(eps)) then eps = 1.d-5
if not(keyword_set(maxits)) then maxits = 10000L
maxits = long(maxits)

anormf = total(abs(f))

omega=0.9                 ; this is unstable if over-relaxed, so under-relax

if (keyword_set(verbose)) then print,'Iter. #, sum |resid|, omega:'

continue = 1
n = 0L
while (continue) do begin
    anorm=0.d0

    resid = ( u(3:imax-2,2:jmax-3)/ds^2 + $
              u(1:imax-4,2:jmax-3)/ds^2 + $
              u(2:imax-3,3:jmax-2)/ds^2 + $
              u(2:imax-3,1:jmax-4)/ds^2 - $
              4.*u(2:imax-3,2:jmax-3)/ds^2 - $
              f(2:imax-3,2:jmax-3))

    anorm = total(abs(resid))

    u(2:imax-3,2:jmax-3) = u(2:imax-3,2:jmax-3) + ds^2*omega*resid/4.

    if (keyword_set(verbose) and ((n)mod(1000) eq 0)) then begin 
        print,n,anorm
        wait,0.5
    endif

    n = n+1L
    if ((n ge MAXITS) or (anorm lt EPS*anormf)) then continue = 0
    
endwhile

if (keyword_set(verbose) and (anorm lt EPS*anormf)) then print,n,anorm 

if (keyword_set(verbose) and (n ge MAXITS)) then $
  message,"MAXITS exceeded in poissond, returning anyway",/info



end