;+
; NAME : HR_CURRENT.pro (procedure)
; PURPOSE :
;	Calculate the vertical current from calibrated Huairou magnetic field data saved in a
;		FITS format file. Write it to a FITS file if set /FITSOUT
; CATEGORY :
;	idl/huairou
; CALLING SEQUENCE :
;	hr_current,fname,current,dir=dir,ext=ext,smth=smth,/fitsout,/fourier
; INPUTS :
;	fname: FITS file name containing the calibrated magnetic field data without extension
; OPTIONAL INPUT PARAMETERS :
;	none
; OUTPUTS :
;	current: calculated vertical current (Ampere/m^2)
; KEYWORD PARAMETERS :
;	dir=dir: directory where the file is located. default: current directory.
;	ext=ext: extension of the file. default: 'fit'
;	fourier: if set, use FFT method to calculate the current. default: difference
;	smth   : smooth the result by a square of SMOOTH by SMOOTH pixel if set
;	fitsout: if set, write the result in a FITS format file named CURRhhmm.fit, where
;			 hh and mm are hour and minute of the observation time in UT.
; COMMON BLOCKS :   none
; SIDE EFFECTS :
; RESTRICTIONS :
;	call hr_calib to calibrate the data before call this procedure to calculate
;	the vertical current
; PROCEDURE :
;	readfits, writefits, sxpar and sxaddpar are called
; MODIFICATION HISTORY :
;   Hui LI, write on 2001/02/26
;------------------------------------------------------------------------
;-
pro hr_current,fname,current,dir=dir,ext=ext,fourier=fourier,fitsout=fitsout,$
	smth=smth

on_error,2
!except=2
t0=systime(1)

;--------------Handle keywords and constants ---------------------------

xsize=5.23 & ysize=3.63  		; FOV in arc min
cail=1.e-4/(7.26e5*4.*!pi*1.e-7)

if keyword_set(ext) then ext='.'+ext else ext='.fit'
if keyword_set(dir) then begin
    if !version.os_family eq 'Windows' then begin
        if strmid(dir,strlen(dir)-1,1) eq '\' then dir=dir else dir=dir+'\'
    endif else begin
        if strmid(dir,strlen(dir)-1,1) eq '/' then dir=dir else dir=dir+'/'
    endelse
endif else begin
    dir=''
endelse
fname=dir+fname+ext
ff=findfile(fname)
!error=0
if ff(0) eq '' then begin
    !error=1
    if keyword_set(dir) then $
        print,'File '+fname+' not found ! Procedure stopped.' else $
        print,'File '+fname+' in current directory not found ! Procedure stopped.'
    return
endif

;-------------- read in data & get position parameters ---------------

array=readfits(fname,head,/sil)
bx=array(*,*,0)
by=array(*,*,1)
bz=array(*,*,2)
nx=sxpar(head,'NAXIS1')
ny=sxpar(head,'NAXIS2')
if sxpar(head,'NAXIS2') eq 4 then spot=array(*,*,3)
time=sxpar(head,'TIME-OBS')
dx=xsize*60./nx & dy=ysize*60/ny

if keyword_set(fourier) then begin
	bxx=complex(bx,fltarr(nx,ny))
	byy=complex(by,fltarr(nx,ny))

	bxx=fft(bxx,-1,/overwrite)
	byy=fft(byy,-1,/overwrite)

	bxi=imaginary(bxx) & bxr=float(bxx)
	byi=imaginary(byy) & byr=float(byy)

	jzr=fltarr(nx,ny) & jzi=jzr
	for m=0,nx-1 do begin
		for n=0,ny-1 do begin
			tmp1=sin(2*!pi*(m+1)/nx)/dx
			tmp2=sin(2*!pi*(n+1)/ny)/dy
			jzr(m,n)=-byi(m,n)*tmp1+bxi(m,n)*tmp2
			jzi(m,n)=+byr(m,n)*tmp1-bxr(m,n)*tmp2
		endfor
	endfor

	jzz=fft(complex(jzr,jzi),1)
	disp=1.e5*cail*float(jzz)
    disp(0,*)=0.0 & disp(nx-1,*)=0.0
    disp(*,0)=0.0 & disp(*,ny-1)=0.0
    method='Fourier'
endif else begin
    dxx=(shift(bx,0,-1)-shift(bx,0,1))/dy/2.0
    dyy=(shift(by,-1,0)-shift(by,1,0))/dx/2.0
    disp=(dyy-dxx)*cail*1.0e5
    disp(0,*)=0.0 & disp(nx-1,*)=0.0
    disp(*,0)=0.0 & disp(*,ny-1)=0.0
    method='Difference'
endelse

if keyword_set(smth) then current=smooth(disp,smth,/edge)*1.e-5 $
else current=disp*1.e-5

print,'MAX =',max(current,min=min),'         MIN =',min

;--------------Write FITS file ----------------------

if keyword_set(fitsout) then begin
	sxaddpar,head,'CURRENT',method,' Method for calculating of vertical current',/pdu
	file=dir+'curr'+strmid(time,0,2)+strmid(time,3,2)+ext
	writefits,file,disp,head
	print,'Result was save in FITS file: '+file
endif

;----------

print,'		 IDL take '+string(systime(1)-t0,format='(f6.2)')+' sec to finish the task.'


end