; NAME :	cff1
; PURPOSE :
;	calculate current-free field based on the longitudinal
;	field (BL) in rdata (planar geometry)
; CATEGORY :
;	solar magnetic field computation package (MAGPACK2)
; CALLING SEQUENCE :
;	cff1, mparam, rdata, bxp, byp
; INPUTS :
;	mparam, rdata
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	faceon: assume that given magnetogram is for BN, regardless
;	of the actual position on the solar disk
; OUTPUTS :
;	bxp, byp : current-free transverse fields
; COMMON BLOCKS : 	none
; SIDE EFFECTS : 	mparam will be modified
; RESTRICTIONS :	none
; PROCEDURE :
;	use FFT to calculate current-free fields
; MODIFICATION HISTORY :
;	T.Sakurai, 1992
;
;----------------------------------------------------------------------
;
pro setbtp, kwx,kwy, alx,aly, cx,cy, amat
;	calculate complex coefficients (cx,cy) to be multiplied
;	on the Fourier transform of BL
;	kwx,kwy = integer indices
;	alx = x-width /2*pi
;	amat = matrix of coordinate transform between data and tangential
;	coordinats
;	
	ci = complex (0.0, 1.0)
	akwx=kwx/alx
	akwy=kwy/aly
	akx=akwx*amat(0,0)
	aky=akwx*amat(0,1)+akwy*amat(1,1)
	akz=sqrt(akx*akx+aky*aky)
	cc=ci*akz-akwx*amat(0,2)-akwy*amat(1,2)
	cx=(akwx+cc*amat(0,2))/(cc*amat(2,2))
	cy=(akwy+cc*amat(1,2))/(cc*amat(2,2))
	return
	end
;
pro cff1, mparam, rdata, bxp, byp, faceon = faceon
;
	nx = mparam.nx
	ny = mparam.ny
	dx = mparam.dx
	dy = mparam.dy

	pi = 3.141592
	nx2 = 2*nx
	ny2 = 2*ny
	alx = dx*nx2/(2.0*pi)
	aly = dy*ny2/(2.0*pi)

	kx = kloc( mparam, /x)
	ky = kloc( mparam, /y)
	kz = kloc( mparam, /z)

        if kx eq -1 or ky eq -1 or kz eq -1 then begin
                print,'(proc cff1  ) error:kx,ky,kz=',kx,ky,kz
                return
        endif

	bl0 = total(rdata(*,*,kz))/(nx*ny)

	rtemp = fltarr(nx2,ny2) + bl0
	rtemp( 0:(nx-1), 0:(ny-1) ) = rdata(*,*,kz)
	cftbz = fft( rtemp, -1)
	cx = cftbz
	cy = cx

	if keyword_set(faceon) then $
		amat= [ 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 ] $
	else begin
		coord_t = settu(mparam)
		amat = coord_t.atmat
	endelse

	bn0 = bl0/amat(2,2)
	bx0 = amat(0,2)*bn0
	by0 = amat(1,2)*bn0

	for j=1, ny2 do begin
	for i=1, nx2 do begin
		if i eq 1 and j eq 1 then goto, endloop
		kwx = i-1
		kwy = j-1
		if kwx gt nx then kwx = kwx -nx2
		if kwy gt ny then kwy = kwy -ny2
		setbtp, kwx,kwy, alx,aly, cx1, cy1, amat
		cx(i-1,j-1) = cx1
		cy(i-1,j-1) = cy1
	endloop:
	endfor
	endfor

	kwy = ny
	for i=1, nx2 do begin
		kwx = i-1
		if kwx gt nx then kwx = kwx -nx2
		setbtp, kwx, kwy, alx,aly, cx1, cy1, amat
		setbtp, kwx,-kwy, alx,aly, cx2, cy2, amat
		cx(i-1, kwy) = 0.5*(cx1+cx2)
		cy(i-1, kwy) = 0.5*(cy1+cy2)
	endfor

	kwx = nx
	for j=1, ny2 do begin
		kwy = j-1
		if kwy gt ny then kwy = kwy -ny2
		setbtp,  kwx,kwy, alx,aly, cx1, cy1, amat
		setbtp, -kwx,kwy, alx,aly, cx2, cy2, amat
		cx(kwx, j-1) = 0.5*(cx1+cx2)
		cy(kwx, j-1) = 0.5*(cy1+cy2)
	endfor

	kwx = nx
	kwy = ny
	setbtp,  kwx, kwy, alx,aly, cx1, cy1, amat
	setbtp, -kwx, kwy, alx,aly, cx2, cy2, amat
	setbtp,  kwx,-kwy, alx,aly, cx3, cy3, amat
	setbtp, -kwx,-kwy, alx,aly, cx4, cy4, amat
	cx(kwx, kwy) = 0.25*(cx1+cx2+cx3+cx4)
	cy(kwx, kwy) = 0.25*(cy1+cy2+cy3+cy4)

	cx(0,0) = complex( 0.0, 0.0 )
	cy(0,0) = complex( 0.0, 0.0 )

	cx = cx*cftbz
	cy = cy*cftbz

	cftbz = fft( cx, 1)
	bxp = float( cftbz( 0:(nx-1), 0:(ny-1) ) ) +bx0

	cftbz = fft( cy, 1)
	byp = float( cftbz( 0:(nx-1), 0:(ny-1) ) ) +by0

	print, '(proc.cff1) end'
	return
end