; NAME :	curlb1
; PURPOSE :
;	calculate curl B from (BX,BY) (planar geometry)
; CATEGORY :
; 	solar magnetic field computation package (MAGPACK2)
; CALLING SEQUENCE :
;	curlb1, mparam, rdata, faceon=faceon, smooth=smooth
; INPUTS :
;	mparam, rdata
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	faceon = assume the region is at disk center
;	smooth = apply 3x3 smoothing filter to (BX,BY)
; OUTPUTS :
;	JN is given in rdata
; COMMON BLOCKS : 	none
; SIDE EFFECTS : 	mparam will be modified
; RESTRICTIONS :	none
; PROCEDURE :
;	use spline fitting to obtain derivatives of BX and BY
; MODIFICATION HISTORY :
;	T.Sakurai, 1992
;
;------------------------------------------------------------------------
;
pro curlb1, mparam, rdata, faceon=faceon, smooth=smooth, keep_bp=keep_bp
;
	kindm = mparam.kindm

	if kindm lt 0 then pol2b, mparam, rdata

	kindmb = kindm mod 10
	kindm = ( kindm/100 ) * 100

	if kindmb eq 1 or kindmb eq 5 then begin
		print, '(proc.curlb1) error: magnetogram is not vector magnetogram'
		return
	endif

	if kindmb eq 2 or kindmb eq 6 then begin
		vector, mparam, rdata, faceon=faceon, keep_bp=keep_bp
		kindmb = mparam.kindm mod 10
	endif

	nx = mparam.nx
	ny = mparam.ny
	dx = mparam.dx
	dy = mparam.dy

	pi = 3.141592
	nx2 = 2*nx
	ny2 = 2*ny

	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 curlb1) error:kx,ky,kz=',kx,ky,kz
                return
        endif
	kj = kloc( mparam, /j)
	if kj eq -1 then begin
		kj = newk( mparam, rdata )
		mparam.nframe = mparam.nframe +1
	endif	
	
	if keyword_set(faceon) then begin
		ainv= [ 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 ] 
		mparam.status.jl = kj
		mparam.kindm = kindm + 10L + kindmb
	endif else begin
		coord_t = settu (mparam)
		amat = coord_t.atmat
		ainv = coord_t.atinv
		mparam.status.jn = kj
		mparam.kindm = kindm + 20L + kindmb
	endelse
	
;	bx = rdata(*,*,kx)*amat(0,0) + rdata(*,*,ky)*amat(1,0) $
;	   + rdata(*,*,kz)*amat(2,0)
;	by = rdata(*,*,kx)*amat(0,1) + rdata(*,*,ky)*amat(1,1) $
;	   + rdata(*,*,kz)*amat(2,1)

	bx = ainv(0,0)*rdata(*,*,kx) + ainv(0,1)*rdata(*,*,ky) $
	   + ainv(0,2)*rdata(*,*,kz)
	by = ainv(1,0)*rdata(*,*,kx) + ainv(1,1)*rdata(*,*,ky) $
	   + ainv(1,2)*rdata(*,*,kz)

	if keyword_set(smooth) then begin
		filter = [ [ 1.0, 1.0, 1.0 ], $
		           [ 1,0, 4.0, 1.0 ], $
		           [ 1.0, 1.0, 1.0 ] ]
		bx = convol(bx, filter)
		by = convol(by, filter)
	endif

	for j=0, ny-1 do begin
		splin1, by(*,j), dx, wk
		rdata(*,j,kj) = wk *amat(0,0)
	endfor

	for i=0, nx-1 do begin
		splin1, reform(bx(i, *)), dy, wk
		rdata(i,*,kj) = rdata(i,*,kj) -wk*amat(1,1)
	endfor

	for j=0, ny-1 do begin
		splin1, bx(*, j), dx, wk
		rdata(*,j,kj) = rdata(*,j,kj) -wk*amat(0,1)
	endfor

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