pro nlfff_tracing,dir,savefile,test,datadir    
;+
; Project     : SDO/HMI, SDO/AIA, IRIS 
;
; Name        : NLFFF_TRACING
;
; Category    : data analysis and automated feature recognition
;
; Explanation : Reads AIA or IRIS images, enables automated loop tracing,
;               and writes loop coordinates into temporary file 
;
; Syntax      : IDL>nlfff_tracing,dir,savefile,test
;
; Inputs      : dir,savefile,test
;
; Outputs     ; saves loop coordinates in temporary file *loop*.dat
;
; History     : 21-Sep-2012, Version 1 written by Markus J. Aschwanden
;		26-Nov-2014, include IRIS SJI images 
;		26-Nov-2014, correct macropixel (nbin-->nxx) in backgr subimage
;		26-Nov-2014: correct x=(i1-crpix1+ii+crval1/cdelt1)*dx_euv	
;		08-Oct-2015: replace ds with dss=sqrt(dxx^2+dyy^2) in hori/vert 
;		08-Oct-2015: nbin_image=para[32], rebinning EUV image 
;		24-Nov-2015: eliminate thresh, thresh2
;               21-Jan-2020, add DATADIR (data directory)
;
; Contact     : aschwanden@lmsal.com
;-

print,'__________________NLFFF TRACING_____________________________
restore,dir+savefile    ;--> PARA
instr	=input.instr
nmag0   =input.nmag_p   ;number of magnetic charges
fov0    =input.fov0	;field-of-view
nsm1	=input.nsm1	;lowpass filter
nstruc	=input.nstruc	;maximum number of analyzed structures 
prox_min=input.prox_min 	;magnetic separation 
lmin_wid=input.lmin_wid	;min loop length / loop width
rmin_wid=input.rmin_wid  ;min curvature radius / loop width
qthresh1=input.qthresh1  ;noise threshold / median(flux)
qthresh2=input.qthresh2  ;filter threshold / median(filter flux)
wave_	=para.wave_
x1_sun  =para.x1_sun
x2_sun  =para.x2_sun
y1_sun  =para.y1_sun
y2_sun  =para.y2_sun

filename=strmid(savefile,0,15)
eps	=0.1
ds	=0.002		;spatial resolution along loop
nsmax	=long(fov0/ds)	;max number of loop points 
ngap	=0
dn_sat	=2.^32
if (instr eq 'AIA') then dn_sat=2.^14		;saturation limit

yellow	=200
red	=125
blue	=75
char	=1
ipos    =strpos(savefile,'.sav')
autofile=datadir+strmid(savefile,0,ipos)+'_auto.tmp'
loopfile=datadir+strmid(savefile,0,ipos)+'_loop.tmp'

;________________WAVELENGTH LOOP_____________________________
nr	=0
nwave	=n_elements(wave_)
nloopw	=long(nstruc/nwave)
openw,1,loopfile
close,1
iwave_disp=0
for iwave=0,nwave-1 do begin
 wave	=wave_[iwave]
 if (wave eq '171') then iwave_disp=iwave
 if (wave eq '193') then iwave_disp=iwave
 imagefile=datadir+filename+'_'+instr+'_'+wave+'.fits'
 fovfile  =datadir+filename+'_'+instr+'_'+wave+'_fov.fits'

;print,'__________READ EUV IMAGE__________________________________
 search =findfile(imagefile,count=file_found)
 if (file_found le 0) then begin
  print,'FILE NOT FOUND : ',imagefile
  goto,end_wave
 endif
 data	=readfits(imagefile,h,/silent)
 dateobs=sxpar(h,'DATE_OBS')
 naxis1	=sxpar(h,'NAXIS1')
 naxis2	=sxpar(h,'NAXIS2')
 cdelt1 =sxpar(h,'CDELT1')
 cdelt2	=sxpar(h,'CDELT2')
 crpix1	=sxpar(h,'CRPIX1')
 crpix2	=sxpar(h,'CRPIX2')
 crval1	=sxpar(h,'CRVAL1')
 crval2	=sxpar(h,'CRVAL2')
 eph	=get_sun(dateobs)
 rsun 	=eph(1)			;rsun=index.rsun_obs in arcseconds 
 rpix   =rsun/cdelt1		;solar radius in number of pixels
 dx_euv	=1./rpix		;pixel size in units of solar radius
 i0_sun =long(crpix1-crval1/cdelt1)
 j0_sun =long(crpix2-crval2/cdelt2)
 nx     =long(1+(x2_sun-x1_sun)/dx_euv)
 ny     =long(1+(y2_sun-y1_sun)/dx_euv)
 image1	=fltarr(nx,ny)
 print,'Dimension of loop image = ',nx,ny
 i1	=long(i0_sun+x1_sun/dx_euv)
 i2	=long(i0_sun+x2_sun/dx_euv)
 j1	=long(j0_sun+y1_sun/dx_euv)
 j2	=long(j0_sun+y2_sun/dx_euv)
 for j=j1,j2 do begin
  if (j ge 0) and (j le naxis2-1) and (j-j1 le ny-1) then begin
   ii1	=(i1 > 0) < (naxis1-1)
   ii2	=(i2 > 0) < (naxis1-1)
   nnx	=ii2-ii1+1
   ii3	=(ii1-i1)>0
   ii4	=ii3+nnx-1
   if (ii4 gt ii3) then image1(ii3:ii4-1,j-j1)=float(data(ii1:ii2-1,j))
  endif
 endfor
 for j=1,ny-1 do begin
  max1=max(image1(*,j-1))
  max2=max(image1(*,j))
  if (max1 gt 0) and (max2 eq 0) then image1(*,j)=image1(*,j-1)
 endfor
 for j=ny-2,0,-1 do begin
  max1=max(image1(*,j+1)) 
  max2=max(image1(*,j))
  if (max1 gt 0) and (max2 eq 0) then image1(*,j)=image1(*,j+1)
 endfor

;_________________BYPASS-FILTERED IMAGE____________________________
 nsm2	=nsm1+2
 if (nsm1 le 2) then image2=image1-smooth(image1,nsm2)
 if (nsm1 ge 3) then image2=smooth(image1,nsm1)-smooth(image1,nsm2)
 writefits,fovfile,image1 

;_________________DISPLAY EUV IMAGE________________________________
if (test ge 1) then begin
 zoom   =long(1024/nx) >1
 if (nx gt 1024) then zoom=0.5
 nxx	=long(nx*zoom)
 nyy	=long(ny*zoom)
 mx	=long(nxx/zoom)
 my	=long(nyy/zoom)
 unit	=0
 window,unit,xsize=nxx,ysize=nyy
 !p.position=[0,0,1,1]
 !x.range=[x1_sun,x2_sun]		;FOV coordinates
 !y.range=[y1_sun,y2_sun]
 !x.style=1
 !y.style=1
 loadct,0,/silent
 plot,[0,0] 
 i3	=long(nx*0.25)			;avoid edge effects (IRIS)
 i4	=long(nx*0.75)
 statistic,image2(i3:i4,i3:i4),zavg,zsig
 c1	=(zavg-0.5*zsig)
 c2	=(zavg+0.5*zsig)
 image3=rebin(image2(0:mx-1,0:my-1),nxx,nyy)
 if (iwave eq iwave_disp) then begin &image4=image3 &c4=c1 &c5=c2 &endif
 tv,bytscl(image3,min=c1,max=c2)
 loadct,5,/silent
endif

;_______________________MAGNETIC SOURCES___________________________
 mm1	=coeff(*,0)
 xm1	=coeff(*,1)
 ym1	=coeff(*,2)
 zm1	=coeff(*,3)
 am1	=coeff(*,4)
 nmag	=n_elements(mm1)
 rm1	=sqrt(xm1^2+ym1^2+zm1^2)
 dm1	=1.-rm1

;_______________________AUTOMATED LOOP TRACING_____________________
 lmin	 =lmin_wid*(nsm1+1)		;minimum length [pixels]
 rmin	 =rmin_wid*(nsm1+1) 		;empirical min curvature radius [pixels]
 para2	 =[nsm1,rmin,lmin,nstruc,nloopw,ngap,qthresh1,qthresh2]
 ipos    =strpos(savefile,'.sav')
 looptracing_auto4,image1,image2,autofile,para2,output,test

;_______________________COORDINATE TRANSFORMATION__________________
 file	=findfile(autofile,count=nfiles)
 if (nfiles le 0) then goto,end_wave
 readcol,autofile,iloop,ix,iy,flux,is,/silent
 il1	=long(min(iloop))
 il2	=long(max(iloop))
 nloop	=il2-il1+1
 for il=il1,il2 do begin
  ind	=where(iloop eq il,nind)
  if (nind eq 0) then goto,skip_loop
  ii	=ix(ind)
  jj	=iy(ind)
  x	=x1_sun+ii*dx_euv
  y	=y1_sun+jj*dx_euv
  f	=flux(ind)
  s	=is(ind)/rpix
;______________ELIMINATION OF LOOPS REMOTE FROM NEXT MAGNETIC SOURCE_______
  dx1    =(x(0)-xm1)
  dx2    =(x(nind-1)-xm1)
  dy1    =(y(0)-ym1)
  dy2    =(y(nind-1)-ym1)
  prox1  =min(sqrt(dx1^2+dy1^2)/dm1)
  prox2  =min(sqrt(dx2^2+dy2^2)/dm1)
  prox	 =prox1 < prox2         ;r1=sqrt(x1^2+y1^2+z1^2), d1=1.-r1
  if (prox gt prox_min) then goto,skip_loop		;eliminate loop
;______________ELIMINATION OF SATURATED LOOPS_____________________
  fmax	 =max(image1(ii,jj))
  if (fmax ge 0.5*dn_sat) then goto,skip_loop		;eliminate saturated fluxes
;______________WRITE LOOP COORDINATES INTO FILE____________________
  openw,1,loopfile,/append
  for ii=0,nind-1 do printf,1,nr,x(ii),y(ii),f(ii),s(ii),wave,$
     format='(I6,2F12.6,F12.2,F10.6,2I6)'
  close,1
  if (test ge 1) then oplot,x,y,color=red,thick=3     ;"good curvi-linear structure"
  nr	=nr+1
SKIP_LOOP:
 endfor
 wait,3
 END_WAVE:
endfor			;wavelength loop
print,'Number of all loops        = ',nr
if (nr eq 0) then stop,'No loops found --> lower minimum length limit LMIN_WID'

;________________READ LOOP FILE____________________________
readcol,loopfile,iloop,xl,yl,fl,sl,wl,/silent
nf	=long(max(iloop)+1)

;________________MAXIMUM LENGTH OF FIELDLINES______________
ns_detmax=0
ns_loop_det  =lonarr(nf)
wave_loop_det=lonarr(nf)
for k=0,nf-1 do begin
 ind	=where(iloop eq k,ns)
 wave_loop_det[k]=wl(ind(0))
 ns_loop_det[k]=ns
 ns_detmax=ns_detmax > ns
endfor

;________________DISPLAY EUV IMAGE (IWAVE_DISP)______________
if (test ge 1) then begin
 loadct,0,/silent
 tv,bytscl(image4,min=c4,max=c5)
 loadct,5,/silent
endif

;________________FIELD LINE ARRAY___________________________
loadct,5,/silent
field_loop_det	=fltarr(nsmax,3,nf)
icount	=0
slope	=fltarr(nf)
for k=0,nf-1 do begin
 ind	=where(iloop eq k,ns)
 ns	=ns < nsmax
 if (ns ge 1) then begin
  xs	=xl(ind[0:ns-1])
  ys	=yl(ind[0:ns-1])
  fs	=fl(ind[0:ns-1])
  field_loop_det(0:ns-1,0,icount)=xs 		;x-coordinate loop
  field_loop_det(0:ns-1,1,icount)=ys		;y-coordinate loop
  field_loop_det(0:ns-1,2,icount)=fs		;flux of loop
  c1	=linfit(xs,ys)
  c2	=linfit(ys,xs)
  slope(k)=abs(c1(1)) < abs(c2(1))	;slope in pixel ratio 
  ns_loop_det[icount]=ns				
  wave_loop_det[icount]=wl(ind(0))
  icount=icount+1
 endif
endfor

;.................Display loops...................................
for k=0,nf-1 do begin
 ind    =where(iloop eq k,ns)
 if (ns ge 1) then begin
  xs    =xl(ind[0:ns-1])
  ys    =yl(ind[0:ns-1])
  zf    =fl(ind[0:ns-1])
  if (test ge 1) then begin
   oplot,xs,ys,color=red,thick=3
   if (slope(k) lt eps) then oplot,xs,ys,color=blue,thick=3
  endif
 endif
endfor
print,'Number of loops            NF    = ',nf 

;______________DISPLAY NLFFF SOURCES______________________________
for im=0,nmag-1 do begin
 mm     =coeff(im,0)
 xm     =coeff(im,1)
 ym     =coeff(im,2)
 zm     =coeff(im,3)
 rm	=sqrt(xm^2+ym^2+zm^2)
 dm 	=1.-rm
 phi    =2.*!pi*findgen(101)/100.
 size	=char*3
 if (test ge 1) then begin
  oplot,[xm,xm],[ym,ym],psym=1,color=yellow 
  oplot,xm+dm*cos(phi),ym+dm*sin(phi),thick=1,color=orange
  d	=min(abs(findgen(nmag)-im))
  if (d eq 0) then oplot,xm+dm*cos(phi),ym+dm*sin(phi),thick=3,color=yellow
 endif
endfor

;________________UPDATE LOOP FILE____________________________
l1	=min(iloop)
l2	=max(iloop)
openw,1,loopfile
for k=l1,l2 do begin
 ind   =where(iloop eq k,ns)
 if (ns ge 1) and (slope(k) ge eps) then begin
  for i=0,ns-1 do begin
   nr=iloop[ind[i]]
   x_i=xl[ind[i]]
   y_i=yl[ind[i]]
   f_i=fl[ind[i]]		 ;contains flux
   s_i=sl[ind[i]]
   w_i=wl[ind[i]]
   printf,1,nr,x_i,y_i,f_i,s_i,w_i,format='(I6,2F12.6,F12.2,F10.6,I6)'
  endfor
 endif
endfor
close,1

;________________SAVE PARAMETERS______________________________
para.ndet = nf
save,filename=dir+savefile,input,para,bzmap,bzfull,bzmodel,coeff,$
     field_loop_det,ns_loop_det,wave_loop_det
print,'parameters saved in file = ',savefile
end