;+
; NAME : HR_CURRENT.pro (procedure)
; PURPOSE :
;	Calculate the vertical current and current helicity (hc) from calibrated Huairou
;	magnetic field data saved in a FITS format file. Write the results to a FITS file
;	if set /FITSOUT. When calculating current helicity, only compute for points with
;	abs(Bl)>20 G, or Bt>100G, or current>0.001A/m^2.
; CATEGORY :
;	idl/huairou
; CALLING SEQUENCE :
;	hr_current,fname,current,helicity,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)
;	helicity: computed current helicity (Gauss^2/m)
; 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 CUHEhhmm.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
;   Hui LI, 2001/06/12 add current helicity computation
;------------------------------------------------------------------------
;-
pro hr_current,fname,current,helicity,dir=dir,ext=ext,fourier=fourier,$
	fitsout=fitsout,smth=smth

on_error,2
!except=2
t0=systime(1)
xsize=5.23 & ysize=3.63  		; FOV in arc min
rsun=6.95990e8

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

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)
bt=sqrt(bx^2.0+by^2.0)
nx=sxpar(head,'NAXIS1')
ny=sxpar(head,'NAXIS2')
if sxpar(head,'NAXIS2') eq 4 then spot=array(*,*,3)
time=sxpar(head,'TIME-OBS')
if sxpar(head,'PROJECT') eq 2 then begin
	dx=1 & dy=1
	cail=1/(4.*!pi*1.e-3)/(0.1*!dtor*rsun)
endif else begin
	cail=1/(7.25e5*4.*!pi*1.e-3)
	dx=xsize*60./nx & dy=ysize*60/ny
endelse

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=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
    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) $
else current=disp

helicity=4*!pi/1000*bz*current
ii=where(abs(bz) le 20 or bt le 100 or abs(current) le 0.001, ncount)
if ncount ne 0 then helicity(ii)=0.0

print,'MAX Current =',max(current,min=min),'         MIN =',min
print,'MAX Helicity=',max(helicity,min=min),'         MIN =',min

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

if keyword_set(fitsout) then begin
	tmp=findgen(nx,ny,2)
	tmp(*,*,0)=current
	tmp(*,*,1)=helicity
	sxaddpar,head,'CURR/HEL',method,$
		' Computing method of vertical current & helicity',/pdu
	file=dir+'CUHE'+strmid(time,0,2)+strmid(time,3,2)+ext
	writefits,file,tmp,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