; this program solves poission's equation:
; the laplacian of the 2-d scalar output 
; field gives the 2-d scalar input field
;
; Approach: for each possible value of (k_x, k_y), find output
; FFT spectrum by dividing input FFT spectrum by (k_x^2 + k_y^2)

function bpoisson, infield, inspect=inspect, outspect=outspect


dims = size(infield, /dim)

nx = dims(0)	; define dimensionality
ny = dims(1)

; define arrays
outfield = fltarr(nx,ny)  	; this will be output array
outspect = complexarr(nx,ny) 	; this spectrum w/FFT gives outfield


; define frequency arrays as IDL does -- for even and odd cases
if ((nx)mod(2) eq 0) then xfreqs = $
[0,(findgen(nx/2) + 1)/(nx),-reverse(findgen(nx/2-1) + 1)/(nx)] $
else xfreqs = $
[0,(findgen(nx/2.-.5) + 1)/(nx),-reverse(findgen(nx/2.-.5) + 1)/(nx)]


if ((ny)mod(2) eq 0) then yfreqs = $
[0,(findgen(ny/2) + 1)/(ny),-reverse(findgen(ny/2-1) + 1)/(ny)] $
else yfreqs = $
[0,(findgen(ny/2.-.5) + 1)/(ny),-reverse(findgen(ny/2.-.5) + 1)/(ny)]

inspect = fft(infield, -1)

for i = 0,nx-1 do begin
	for j = 0,ny-1 do begin

		; laplacian of 0-modes is zero, so start w/non-zero modes.
		if ((i ne 0) or (j ne 0)) then $
		outspect(i,j) = -(2*!pi)^(-2)/ $
				((xfreqs(i))^2+(yfreqs(j))^2) $
				*inspect(i,j)

	endfor
endfor

outfield = real_part(fft(outspect,1))

return, outfield

end