; Computes Coulomb-gauge vector potential, A_p, consistent 
; with given normal field array B_z, by Green's function.   
;===========================================================

pro vectpot_gf, bz, ax=ax, ay=ay, dx=dx, $
                    reconbz=reconbz, plot=plot

nparms=N_Params()
if (nparms eq 0) then begin
   Message,'Usage:   vectpot_gf, bz, ax=ax, ay=ay, dx=dx, /plot',/info
   Message,'',/info
   Message,'Purpose: Uses Greens Function method to find Coulomb-gauge ',/info
   Message,'         vector potential, A, consistent with the potential ', $
     /info
   Message,'         magnetic field matching the input Bz.',/info
   Message,'',/info
   Message,'Input:   bz = (nx,ny) fltarr of z-comp of B [usu. in gauss]',/info
   Message,'Input:   dx = pixel size [usu. in cm], default = 1',/info
   Message,'',/info
   Message,'Output:  ax = (nx,ny) fltarr of x-comp of A (units of [B][dx])', $
     /info
   Message,'Output:  ay = (nx,ny) fltarr of y-comp of A (units of [B][dx])', $
     /info
   Message,'Output:  reconbz = curl of A, for comparison with Bz [G]',/info
   Message,'',/info
   Message,'Keyword: plot = set to generate four surface plots, of ',/info
   Message,'                Bz, reconbz (curl A), Ax, & Ay',/info
   Message,'',/info
   Message,'History: Written, B. Welsch, 04 Jan. 2004',/info
   Message,'         Commented, B. Welsch, 10 Jan 2007',/info
   return
endif

if (not(keyword_set(dx))) then dx = 1.

nx = n_elements(bz(*,0))
ny = n_elements(bz(0,*))

ax = fltarr(nx,ny)
ay = ax

x1d = findgen(nx)*dx
y1d = findgen(ny)*dx

x2d = x1d#replicate(1,ny) 
y2d = replicate(1,nx)#y1d 

for i=0.,float(nx*ny)-1. do begin

    xi = (i)mod(nx)*dx
    yi = floor(i/float(nx))*dx
        
    rx =  x2d - xi  
    ry =  y2d - yi
    r2 =  rx^2 + ry^2
    nzr = where(r2 ne 0, n_nzr)
        
    ax(i) =  total(bz(nzr)*ry(nzr)/r2(nzr))*dx^2/(2*!pi) 
    ay(i) = -total(bz(nzr)*rx(nzr)/r2(nzr))*dx^2/(2*!pi)

    ;if ((i)mod(nx) eq 0) then print,i

endfor

daydx = fltarr(nx,ny)
daxdy = fltarr(nx,ny)
for i=0.,ny-1. do daydx(*,i) = deriv(x1d,ay(*,i))
for i=0.,nx-1. do daxdy(i,*) = deriv(y1d,ax(i,*))
reconbz = daydx - daxdy


; CHECK accuracy, if called for
;===============================
if (keyword_set(plot)) then begin
    window,0
    surface,bz,tit='B!dz!n',chars=2
    window,1
    surface,reconbz,tit='Reconstructed B!dz!n',chars=2
    window,2
    surface,ax,tit='A!dx!n',chars=2
    window,3
    surface,ay,tit='A!dy!n',chars=2
endif


end