PRO qfactor, Bx, By, Bz, xreg=xreg, yreg=yreg, zreg=zreg, factor=factor, $
                         twistFlag=twistFlag, csFlag=csFlag, $
                         nbridges=nbridges,  step=step, odir=odir, fstr=fstr
              
;+
; PURPOSE:
;   Calculate the squashing factor Q at the photosphere or a cross section
;   or a box volume, given a 3D magnetic field with uniform grids
; 
;   For details see:
;     Liu et al. (2016, ApJ)  
;
; INPUTS               
;   Bx, By, Bz:      3D magnetic field 
;   
; OPTIONAL INPUTS  
;   xreg, yreg, zreg:  pixel coordinates of the field region, in arrays of two elements
;                      default is to include all available pixels		               
;   	             --- If (xreg[0] ne xreg[1]) AND (yreg[0] ne yreg[1]) AND (zreg[0] ne zreg[1]) AND NOT csFlag
;   	                    caculate Q in a box volume
;   	                    invoke vFlag
;  		             --- If (xreg[0] eq xreg[1]) OR (yreg[0] eq yreg[1]) or (zreg[0] eq zreg[1])
;  		                    calculate Q in a cross section perpendicular to X or Y or Z axis
;  		                    invoke cFlag 
;  		             --- If zreg=[0, 0], calculate Q at the photosphere
;  		                    invoke z0Flag 
;  		             --- If csFlag is set, see below
;  		                    invoke cFlag 
;      
;   factor:     to bloat up the original resolution; i.e., grid spacing = 1/factor; default is 4 
;   
;   nbridges:   number of processors to engage; default is 8
;
;   twistFlag:  to calulate twist number Tw; see Liu et al. (2016, ApJ); default is 0
;   
;   csFlag:     to calulate Q in a cross section defined by three points; default is 0
;               point0=[xreg[0],yreg[0],zreg[0]]
;               point1=[xreg[1],yreg[1],zreg[0]]
;               point2=[xreg[0],yreg[0],zreg[1]]
;				   
;		step:       step size in tracing field lines; default=1/factor
;		
;		odir:       directory to save the results
;		
;		fstr:       filename to save the results 
; 
; OUTPUTS:
;   q/qcs/q3d:  squashing factor   
;   slogq:      sign(Bz) x log_10(Q), calulated with all the field lines
;   slogq_orig: only include field lines with both footpoints on the photoshere
;   twist:      twist number = \int \mu_0 J \cdot B /4\pi B^2 ds
;   rboundary:  nature of the ends of field lines
;               1 - end at zmin
;               2 - end at zmax
;               3 - end at ymin
;               4 - end at ymax
;               5 - end at xmin
;               6 - end at xmax
;               0 - others
;   
; PROCEDURES USED:
; ------IDL procedures used; those included in SolarSoft are not specified
;   check_slash()
;   doppler_color
;   doppler_color_mix            
;   sign2d()
;   rev_true_image()
;      
; ------FORTRAN procedures used
;   qfactor.f90
;   trace_bline.f90
;      if compiled by ifort:  ifort -o qfactor qfactor.f90 -parallel -openmp -mcmodel=large
;      
;      
; MODIFICATION HISTORY:
;   Developed by R. Liu and J. Chen @ USTC
;   
;   June 30, 2014, R. Liu @ USTC, IDL edition 
;   July 2,  2014 R. Liu, introduce nchunks > nbridges to utilize idle child processes
;   April 21, 2015 R. Liu and J. Chen, deal with field lines pass through the boudnary other than bottom
;   April 29, 2015 R. Liu and J. Chen, further optimization on shared memory
;   Arpil 30, 2015 R. Liu, correct the size of qmap
;   Apr 27, 2015 R. Liu, qcs
;   June 15, 2015 J. Chen, Fortran Edition, modify foot points
;   July 8, 2015 J. Chen, q3d
;	  Otc 29, 2015 J. Chen, deal with field lines touching the cut plane: use the plane quasi-perp to field line;
;	  
;	  Nov 1, 2015 J. Chen,
;		(1) combine qcs and qfactor on one procedure;  
;		(2) the cross section can be parallel to photosphere;
;		(3) set tmp_dir to avoid possible conflict by many users run same code on a PC;
;		(4) further optimization on thread task allocation method;
;		(5) deal with marginal Q.
;	  
;	  Nov 4, 2015 J. Chen,
;		(1) when field lines touching the cut plane, qcs[i, j]= min([q_cut, q_vertical]), 
;		    q_cut caculated at the cut plane, q_vert caculated at the plane which quasi vertical of field line.
;		(2) introduce rboundary3d (byte array, for save disk)
;		    rboundary3d=byte(rsboundary3d+7*reboundary3d)
;		So if rboundary3d[i, j, k] eq 8B, the definition of q3d[i, j, k] is based on the bottom surface.
;		
;		Jun 22, 2016
;		add map of field line length - J.C.
;
;
;   This software is provided without any warranty. Permission to use,
;   copy, modify. Distributing modified or unmodified copies is granted,
;   provided this disclaimer and information are included unchanged.
;-

sbx=size(Bx)
sby=size(By)
sbz=size(Bz)
if sbx[0] ne 3 or sby[0] ne 3 or sbz[0] ne 3 then message,'Bx, By and Bz must be 3D arrays!'
if sbx[1] ne sby[1] or sbx[1] ne sbz[1] or $
   sbx[2] ne sby[2] or sbx[2] ne sbz[2] or $
   sbx[3] ne sby[3] or sbx[3] ne sbz[3] then message,'Bx, By and Bz must have the same dimensions!' 
nx=sbz[1] & ny=sbz[2] & nz=sbz[3]
b2dz=Bz[*,*,0]

if keyword_set(xreg) then xreg=xreg else xreg=[0, nx-1]
if keyword_set(yreg) then yreg=yreg else yreg=[0, ny-1]
if keyword_set(zreg) then zreg=zreg else zreg=[0, 0]
if n_elements(xreg) ne 2 or n_elements(yreg) ne 2 or n_elements(zreg) ne 2 then $
   message,'xreg, yreg and zreg must be 2-element arrays!'

if keyword_set(factor) then factor=long(factor) else factor=4L
factor_str='f'+string(factor,format='(i02)')
if keyword_set(step) then step=step else step=1.0D/factor

qx=abs(xreg[1]-xreg[0])*factor+1
qy=abs(yreg[1]-yreg[0])*factor+1
qz=abs(zreg[1]-zreg[0])*factor+1

; calculate Q in a cross section parallel to the photosphere
zFlag = xreg[0] ne xreg[1] AND yreg[0] ne yreg[1] AND zreg[0] eq zreg[1] 
if zFlag then begin 
	q1=qx & q2=qy
endif
if zFlag and zreg[0] eq 0 then z0Flag=1B else z0Flag=0B

; calculate Q in a cross section perpendicular to X- or Y-axis 
xFlag = xreg[0] eq xreg[1] AND yreg[0] ne yreg[1] AND zreg[0] ne zreg[1] 
if xFlag then begin 
	q1=qy & q2=qz
endif

yFlag = xreg[0] ne xreg[1] AND yreg[0] eq yreg[1] AND zreg[0] ne zreg[1]
if yFlag then begin 
	q1=qx & q2=qz
endif

; calculate Q in a cross section perpendicular to the photosphere
if keyword_set(csFlag) then csFlag=1B else csFlag=0B
cFlag = xFlag OR yFlag OR zFlag OR csFlag AND NOT z0Flag
 

if csFlag then begin
	point0=[xreg[0],yreg[0],zreg[0]]
	point1=[xreg[1],yreg[1],zreg[0]]
	point2=[xreg[0],yreg[0],zreg[1]]
	if 	((total(point0 eq point1) eq 3 ) or (total(point0 eq point2) eq 3 ) or (total(point1 eq point2) eq 3 )) then $
      message,'Something is wrong with the cross section .......'
	q1=sqrt(total(Double(point1-point0)*(point1-point0)))*factor+1
	q2=sqrt(total(Double(point2-point0)*(point2-point0)))*factor+1
endif else begin
;	point* not use here but need input to head.txt
	point0=fltarr(3) & point1=fltarr(3) & point2=fltarr(3)	
endelse


if cFlag then begin
	if (q1 eq 1 OR q2 eq 1) then message,'Something is wrong with the cross section .......'
end

; caculate Q in a box volume 
vFlag = xreg[0] ne xreg[1] AND yreg[0] ne yreg[1] AND zreg[0] ne zreg[1] AND NOT csFlag
if vFlag then begin
	q1=qx & q2=qy
	if qx eq 1 OR qy eq 1 OR qz eq 1 then message,'Something is wrong with the box volume .......'
endif

if vflag then head_str='q3d_' else head_str='qcs_'
  cut_str='' 
  if xFlag then cut_str='_x'+string(xreg[0],format='(I0)')
  if yFlag then cut_str='_y'+string(yreg[0],format='(I0)')
  if zFlag then cut_str='_z'+string(zreg[0],format='(I0)')
if keyword_set(fstr) then fstr=fstr+'_' + head_str + factor_str + cut_str $
else fstr = head_str + factor_str + cut_str

if keyword_set(odir) then odir = check_slash(check_slash(odir)+'qfactor') $
                     else odir= check_slash(check_slash(curdir())+'qfactor')
if not dir_exist(odir) then file_mkdir,odir
print,'Results to be saved in ', odir+fstr+'.sav' 

tmp_dir= check_slash(odir+'tmp')
if not dir_exist(tmp_dir) then file_mkdir, tmp_dir

if keyword_set(nbridges) then nbridges=nbridges else nbridges=8
max_threads=!CPU.HW_NCPU
nbridges=nbridges < max_threads  

if keyword_set(twistFlag) then twistFlag=1B else twistFlag=0B

;###############################################################

dummy=file_search(tmp_dir+'*.txt',count=nf)
if nf ne 0 then spawn,'rm '+tmp_dir+'*.txt'
dummy=file_search(tmp_dir+'*.bin',count=nf)
if nf ne 0 then spawn,'rm '+tmp_dir+'*.bin'

cur_device=!D.name
SET_PLOT, 'Z' 
loadct,0
DEVICE, SET_RESOLUTION=[nx,ny]
tv,bytscl(b2dz,min=-1000,max=1000,/nan)

if csflag then begin
	plots,[xreg[0],xreg[1]],[yreg[0],yreg[1]],/dev
endif else begin
	plots,[xreg[0],xreg[1],xreg[1],xreg[0],xreg[0]],[yreg[0],yreg[0],yreg[1],yreg[1],yreg[0]],/dev
endelse
im=tvread()
write_png,odir+fstr+'_bz.png',im

;,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
;  transmit data to Fortran

get_lun,unit
openw, unit, tmp_dir+'head.txt'
printf, unit, nx, ny ,nz , nbridges, factor
printf, unit, float(xreg), float(yreg), float(zreg), long(twistFlag), step
printf, unit, qx, qy, qz, q1, q2
printf, unit, long(csflag), float(point0), float(point1), float(point2)
close, unit

openw, unit, tmp_dir+'b3d.bin'
writeu, unit, Bx, By, Bz
close, unit

;,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
; calculating in Fortran
cdir=curdir()
CD, tmp_dir
tnow=systime(1)
spawn, 'qfactor'  ; if not known by the system, specify the path
tend=systime(1)
tcalc=tend-tnow
if (tcalc ge 3600) then begin
	print,string(tcalc/3600.0,format='(f0.2)'),' hours elapsed for calculating qfactor'
endif else begin 
	if (tcalc ge 60) then begin 
		print,string(tcalc/60.0,format='(f0.2)'),' minutes elapsed for calculating qfactor' 
	endif else begin 
		print,string(tcalc,format='(f0.2)'),' seconds elapsed for calculating qfactor'
	endelse
endelse

; ################################### retrieving results ###################################################### 
; Q at the photosphere ----------------------------------------------------------------------------------------------
IF z0Flag THEN BEGIN
  ; get data From Fortran
	q=fltarr(q1,q2)
	rboundary=lonarr(q1,q2)
  length=fltarr(q1,q2)
	openr,unit, tmp_dir+'qfactor0.bin'
	readu, unit, q, rboundary, length
	if twistFlag then begin
		twist=fltarr(qx,qy)
		readu, unit, twist
	endif
	close, unit
	free_lun, unit

  ; save results and slogQ map 
	sign_congrid_b2dz=sign2d(congrid(b2dz[xreg[0]:xreg[1],yreg[0]:yreg[1]], qx, qy, /interp, /minus_one))
	
  ; only deal with field lines with both footpoints anchored at the bottom
	ss=where(rboundary ne 1,complement=nss)
	qtmp=q
	if(ss[0] ne -1) then qtmp[ss]=!values.F_NAN

	
  DEVICE, SET_RESOLUTION=[q1,q2]
  loadct, 0
  tv,bytscl(length,min=0,max=sqrt(nx^2.0+ny^2.0+nz^2.0),/nan)  
  im=tvread(/truecolor)
  write_png,odir+fstr+'_length.png',rev_true_image(im)
	
	doppler_color 
  slogq_orig=sign_congrid_b2dz*alog10(qtmp>1)
	tv,bytscl(slogq_orig,min=-5,max=5,/nan)
	im=tvread(/truecolor)
	write_png,odir+fstr+'_slogq_orig.png',rev_true_image(im)

  ; include all field lines
	doppler_color_mix
	slogq=sign_congrid_b2dz*alog10(q>1.)
	slogq_mix=slogq
	if(ss[0] ne -1) then slogq_mix[ss]=bytscl(slogq[ss],min=-5,max=5,/nan,top=127)+128B  ; green-white-yellow
	if(nss[0] ne -1) then slogq_mix[nss]=bytscl(slogq[nss],min=-5,max=5,/nan,top=127)     ; blue-white-red
	tv,slogq_mix
	im=tvread(/truecolor)
	write_png,odir+fstr+'_slogq_mix.png',rev_true_image(im)
  
  ; twist map
	if twistFlag then begin
		doppler_color
	  tv,bytscl(twist,min=-2,max=2,/nan)
	  im=tvread(/truecolor)  
	  write_png,odir+fstr+'_twist'+'.png',rev_true_image(im)
	endif
	
  if twistFlag then save,filename=odir+fstr+'.sav', slogq, slogq_orig, q, length, twist, rboundary, xreg, yreg, zreg, factor $ 
               else save,filename=odir+fstr+'.sav', slogq, slogq_orig, q, length, rboundary, xreg, yreg, zreg, factor
ENDIF 

; Q in the 3D box volume ----------------------------------------------------------------------------------------------
IF vflag THEN BEGIN
  ;get data From Fortran
	q3d=fltarr(qx,qy,qz)
	rboundary3d_tmp=lonarr(qx,qy,qz)	
	openr,unit, tmp_dir+'q3d.bin'
	readu, unit, q3d, rboundary3d_tmp
	close, unit
		
	rboundary3d=byte(rboundary3d_tmp)
	dummy=(temporary(rboundary3d_tmp))[0,0,0]
  ;	rboundary3d=rsboundary3d+7*reboundary3d  (It has been caculated in Fortran)
  ;	So if rboundary3d[i, j, k] eq 8B, the definition of q3d[i, j, k]  is based on the bottom surface

	if twistFlag then begin
		openr,unit, tmp_dir+'twist3d.bin'	
		twist3d=fltarr(qx,qy,qz)
		readu, unit, twist3d
		close, unit
	endif

  ; save resutls 
	if twistFlag then save,filename=odir+fstr+'.sav', q3d, rboundary3d, twist3d, xreg, yreg, zreg, factor $
			         else save,filename=odir+fstr+'.sav', q3d, rboundary3d, xreg, yreg, zreg, factor
			         
	sign_congrid_b2dz=sign2d(congrid(b2dz[xreg[0]:xreg[1],yreg[0]:yreg[1]], qx, qy, /interp, /minus_one))
  
  ; only deal with field lines with both footpoints anchored at the bottom  
  qtmp=q3d[*,*,0]
  rb=rboundary3d[*,*,0]
  ss=where(rb ne 8,complement=nss)
  if(ss[0] ne -1) then qtmp[ss]=!values.F_NAN

  doppler_color 
  slogq_orig=sign_congrid_b2dz*alog10(qtmp>1)
  
  DEVICE, SET_RESOLUTION=[qx,qy]  
  tv,bytscl(slogq_orig,min=-5,max=5,/nan)
  im=tvread(/truecolor)
  write_png,odir+fstr+'_z0_slogq_orig.png',rev_true_image(im)

  ; include all field lines
  doppler_color_mix
  qtmp=q3d[*,*,0]
  slogq=sign_congrid_b2dz*alog10(qtmp>1)
  slogq_mix=slogq
  if(ss[0] ne -1) then slogq_mix[ss]=bytscl(slogq[ss],min=-5,max=5,/nan,top=127)+128B  ; green-white-yellow
  if(nss[0] ne -1) then slogq_mix[nss]=bytscl(slogq[nss],min=-5,max=5,/nan,top=127)     ; blue-white-red
  tv,slogq_mix
  im=tvread(/truecolor)
  write_png,odir+fstr+'_z0_slogq_mix.png',rev_true_image(im)		         
ENDIF 

; Q at the cross section ----------------------------------------------------------------------------------------------
IF cFlag THEN BEGIN

	qcs=fltarr(q1,q2)
	rsboundary=lonarr(q1,q2)		
	reboundary=lonarr(q1,q2)
	length=fltarr(q1,q2)
  openr,unit, tmp_dir+'qcs.bin'
  readu, unit, qcs, rsboundary, reboundary, length			
	if twistFlag then begin
		twist=fltarr(q1,q2)
		readu, unit, twist
	endif
	close, unit

	
	DEVICE, SET_RESOLUTION=[q1,q2]
	loadct,0
	tv,bytscl(length,min=0,max=sqrt(nx^2.0+ny^2.0+nz^2.0),/nan)
  im=tvread(/truecolor)
  write_png,odir+fstr+'_length.png',rev_true_image(im)
	
	logq=alog10(qcs>10.) 
	tv,bytscl(logq,max=5,/nan)
	im=tvread(/truecolor)
	write_png,odir+fstr+'_logq.png',rev_true_image(im)

	qcs_orig=qcs
	ss=where(rsboundary ne 1,complement=nss)
	if(ss[0] ne -1) then qcs_orig[ss]=!values.F_NAN
	logq=alog10(qcs_orig>10.) 
	tv,bytscl(logq,max=5,/nan)
	im=tvread(/truecolor)
	write_png,odir+fstr+'_logq_orig.png',rev_true_image(im)

	if twistFlag then begin
		doppler_color
		tv,bytscl(twist,min=-1.5,max=1.5,/nan)
		im=tvread(/truecolor)  
		write_png,odir+fstr+'_twist.png',rev_true_image(im)
	endif	   

	if twistFlag then save,filename=odir+fstr+'.sav', qcs, qcs_orig, length, twist, rsboundary, reboundary, xreg, yreg, zreg, factor $
	else save,filename=odir+fstr+'.sav', qcs, qcs_orig, length, rsboundary, reboundary, xreg, yreg, zreg, factor
	
ENDIF
	
; hourse keeping
cd, cdir
spawn, 'rm -r '+tmp_dir
dummy=temporary(odir)
dummy=temporary(fstr)
set_plot, cur_device
free_lun, unit
END