; NAME :	cff1n
; 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, Bl, bxp, byp
; INPUTS :
;	Bl   -- longitudinal magnetic field
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	position  :  position on the solar disk, [ ew, ns ]
;			 disk center => [0,0],  limb => +-1.
;	dx    : scale of 1 pixel in unit of 10000km
; 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
;	K.Ichimoto, 92/10/17   modify cff1.pro
;
;----------------------------------------------------------------------
;
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 cff1n, Bl, bxp, byp, position = position
;
	s=size(Bl)
	nx = s(1)
	ny = s(2)

	nx2 = 2*nx
	ny2 = 2*ny
	alx = nx2/(2.0*!pi)
	aly = ny2/(2.0*!pi)

	bl0 = total(Bl)/(nx*ny)

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

	amat=fltarr(3,3)
	if keyword_set(position) then begin
		amat = matrix_tn( position )
	endif else begin
		amat(*,0)= [ 1.0, 0.0, 0.0 ]
		amat(*,1)= [ 0.0, 1.0, 0.0 ]
		amat(*,2)= [ 0.0, 0.0, 1.0 ]
	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.cff1n) end'
	return
end