;+
; NAME: epot_fft.pro
;
; USAGE: pot_e = epot_fft(bz)
;
; PURPOSE: This function finds the volume integral of the potential
;     magnetic field squared, given normal field on the bottom of a
;     Cartesian box and assuming periodicity in the horizontal
;     directions.  With suitable normalization, this gives the
;     magnetic energy (e.g., 1/8\pi).
;
; METHOD: With periodic Bz, solve for potential field w/FFT's.
;
;     If B = -grad(S), w/B & S periodic in x & y and decaying
;     in z, and div(B) = 0 ==> laplacian(S) = 0, then
;
;      S(x,y,z) = sum_lm s_lm exp(i(k_l x + k_m y)- k_lm z) &
;     Bz(x,y,z) = sum_lm b_lm exp(i(k_l x + k_m y)- k_lm z) ,
;
;     with b_lm = + k_lm s_lm, and k_lm = sqrt(k_l^2 + k_m^2).
;
; NOTES: Error handling is non-existent.
;  
; HISTORY: Written 08 Jun 2005 by BTW.
;
;-

function epot_fft, bz, potential=potential

nx1 = n_elements(bz(*,0))
ny1 = n_elements(bz(0,*))

; Define arrays on z = 0 plane
;========================================
potspect = complexarr(nx1,ny1) 	; FFT spectrum of scalar potential


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

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

bzspect = fft(bz, -1) ; forward FFT of Bz

for i = 0,nx1-1 do begin
    for j = 0,ny1-1 do begin

	; k_n = sqrt(k_l^2 + k_m^2)
	;==========================
        zfreq = sqrt(xfreqs(i)^2 + yfreqs(j)^2) 

	; s_nlm = b_nlm / k_n
	;==================================
        if (zfreq ne 0) then potspect(i,j) = bzspect(i,j)/(2*!pi)/ zfreq 

    endfor
endfor

potential = real_part(fft(potspect,1))
emag =  .5*total(potential*bz)

return,emag

end