pro nlfff_disp_v1,dir,catfile,iflare,it,run
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : NLFFF_DISP
;
; Category    : display of magnetic data modeling
;
; Explanation : displays magnetogram and observed and fitted magnetic
;		field lines in top-view and side-view projection.
;
; Syntax      : IDL>nlfff_disp,dir,catfile,iflare,it,run,io
;
; Inputs      : dir,catfile
;		it	  = time step (1,...,nt)
;               run       = label to mark particular run
;
; Outputs     ; postscript file <imagefile>.ps
;
; History     : 16-May-2011, Version 1 written by Markus J. Aschwanden
;               30-Jan-2012, Version 2 released to SSW
;		 9-Oct-2016, display both screen (io=0) and postscript (io=3) 
;             ;  1-May-2017, hmin=input.hmin, hmax=input.hmax
;
; Contact     : aschwanden@lmsal.com
;-

print,'===================NLFFF_DISP====================='
readcol,dir+catfile,event_,dateobs_,noaa_,helpos_,instr_,$
  w1_,w2_,w3_,w4_,w5_,w6_,w7_,w8_,w9_,w10_,w11_,$
  format='(I,A,A,A,A,A,A,A,A,A,A,A,A,A,A,A)',/silent
ind     =where(event_ eq it,nind)
if (nind le 0) then stop,'NLFFF_EVENT : NO MATCHING EVENT FOUND'
i       =ind(0)                 ;selected entry
string0,3,it,it_nr
savefile =dateobs_(i)+'_'+instr_(i)+'_'+it_nr+'_'+run+'.sav'

;_______________________READ DATA________________________
print,'read savefile = ',dir+savefile
restore,dir+savefile    ;--> INPUT,PARA
instr	=input.instr
noaa	=input.noaa
helpos	=input.helpos
fov0    =input.fov0	;field-of-view [solar radii]
nmag_p  =input.nmag_p  	;number of magnetic charges
nmag_np =input.nmag_np  ;number of magnetic charges
nsm1    =input.nsm1    	;lowpass filter
nstruc  =input.nstruc   	;max number of structures 
nitmin	=input.nitmin
amis	=input.amis
prox_min=input.prox_min
lmin_wid=input.lmin_wid
rmin_wid=input.rmin_wid
qthresh1=input.qthresh1
qthresh2=input.qthresh2
ngap	=input.ngap
hmin	=input.hmin
hmax	=input.hmax

event	=para.event
dateobs	=para.dateobs
wave_	=para.wave_
filename=strmid(savefile,0,15)
x1_sun  =para.x1_sun
x2_sun  =para.x2_sun
y1_sun  =para.y1_sun
y2_sun  =para.y2_sun
dpix_euv=para.dpix_euv
dpix_mag=para.dpix_mag
qe_rebin=para.qe_rebin
qe_model=para.qe_model
iter	=para.iter
cpu	=para.cpu
angle	=para.angle
angle2	=para.angle2
merit	=para.merit
qiso_corr=para.qiso_corr
ndet	=para.ndet
nitmax	=para.nitmax

ds	=0.002		;step aplong field line
nsmax   =long(fov0/ds)  
nseg    =9
nwave	=n_elements(wave_)
ct	=1
blue	=50
purple	=75
red	=125
yellow	=200 
dim     =size(bzfull)
nx      =dim(1)
ny      =dim(2)
x       =x1_sun+dpix_mag*findgen(nx)
y       =y1_sun+dpix_mag*findgen(ny)
nf	=n_elements(ns_loop)

;____________________DISPLAY________________________________________
for io=0,3,3 do begin        ;0=screen, 3=color postscript file
form    =1      ;0=landscape, 1=portrait
char    =0.9    ;character size
if (io eq 0) then char=1.2
fignr	=''
ref     =''
dlat    =5
unit	=0
ipos	=strpos(savefile,'.sav')
plotname=dir+strmid(savefile,0,ipos)
fig_open,io,form,char,fignr,plotname,unit

;_____________________DISPLAY MAGNETOGRAM______________________________
loadct,ct,/silent
x1_     =0.1
x2_     =0.65
y1_     =0.02
y2_     =0.45
!p.title=' '
!p.position=[x1_,y1_,x2_,y2_]
!x.range=[x[0],x[nx-1]]
!y.range=[y[0],y[ny-1]]
!x.title=''
!y.title=''
!x.style=1
!y.style=1
if (io eq 0) then begin
 nxw   =long(!d.x_vsize*(x2_-x1_)+0.5)
 nyw   =long(!d.y_vsize*(y2_-y1_)+0.5)
 image =congrid(bzfull,nxw,nyw)			;magnetogram 
endif
if (io ne 0) then image=bzfull 			;bzmap
plot,[0,0],[0,0],charsize=0.8
statistic,bzfull,c0,csig
c1      =c0-5*csig
c2      =c0+1*csig
tv,bytscl(image,min=c1,max=c2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal
coordsphere,1,0,0,0,0,dlat
!noeras=1

;______________________SELECT DISPLAYED LOOPS__________________________
dev_loop=fltarr(nf)
for k=0,nf-1 do begin
 ifit   =k*nseg+findgen(nseg)
 dev_loop(k)=median(dev_deg(ifit))
endfor
nf_sel	=nf
isort	=sort(dev_loop)
k_sel	=isort(0:nf_sel-1)

;______________________DISPLAY FIELD LINES_______________________
loadct,5,/silent
len1    =fltarr(nf)
len2    =fltarr(nf)
for k_sel=0,nf_sel-1 do begin
 k	=isort(k_sel)
 np	=ns_field[k]
 xl	=field_lines(0:np-1,0,k)
 yl	=field_lines(0:np-1,1,k)
 oplot,xl,yl,color=red,thick=1
endfor

;______________________DISPLAY TRACED LOOPS_________________________
for k_sel=0,nf_sel-1 do begin
 k	=isort(k_sel)
 ns	=ns_loop[k]
 xl	=field_loop(0:ns-1,0,k)
 yl	=field_loop(0:ns-1,1,k)
 zl	=field_loop(0:ns-1,2,k)
 oplot,xl,yl,color=yellow,thick=1
endfor

;_____________________DISPLAY ROTATED MAGNETOGRAM_________________________
x0	=(min(x)+max(x))/2.
y0	=(min(y)+max(y))/2.
z0	=sqrt(1.-x0^2-y0^2)
yrot	=z0+(y-y0)
bzrot	=bzfull*0.
bmin	=min(bzfull)
bmax	=max(bzfull)
for j=ny-1,0,-1 do begin
 r	=sqrt(x^2+yrot(j)^2)
 ind	=where(r gt 1,nind)
 if (nind ge 1) then bzrot(ind,j)=bmin+(bmax-bmin)*0.4*exp((r(ind)-1.)/0.2)
endfor
for j=ny-1,0,-1 do begin
 z	=sqrt(1.-x^2-y(j)^2)
 iz	=long((z-yrot[0])/(y[1]-y[0])) < (ny-1) 
 for ix=0,nx-1 do if (iz(ix) ge 0) and (iz(ix) le nx-1) then $
  bzrot(ix,iz(ix))=bzfull(ix,j)
endfor

loadct,ct,/silent
y1_     =0.47
y2_     =0.92
!p.title=' '
!p.position=[x1_,y1_,x2_,y2_]
!x.range=[x[0],x[nx-1]]
!y.range=[yrot[0],yrot[ny-1]]
!x.title=''
!y.title=''
!x.style=1
!y.style=1
if (io eq 0) then begin
 nxw   =long(!d.x_vsize*(x2_-x1_)+0.5)
 nyw   =long(!d.y_vsize*(y2_-y1_)+0.5)
 image =congrid(bzrot,nxw,nyw)			;magnetogram 
endif
if (io ne 0) then image=bzrot 			;bzmap
plot,[0,0],[0,0],charsize=0.8
statistic,bzfull,c0,csig
tv,bytscl(image,min=c1,max=c2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal
coordsphere,1,0,0,90,0,dlat

;_______________________DISPLAY ROTATED FITTED LINES____________________________
loadct,5,/silent
xf0	=fltarr(nf)
yf0	=fltarr(nf)
zf0	=fltarr(nf)
i1_	=fltarr(nf)
i2_	=fltarr(nf)
zmin	=2.0
zmax	=0.0
len3	=len1>len2
for k_sel=0,nf_sel-1 do begin
 k	=isort(k_sel)
 np	=ns_field[k]
 xl	=field_lines(0:np-1,0,k)
 yl	=field_lines(0:np-1,1,k)
 zl	=field_lines(0:np-1,2,k)
 zmin	=zmin < min(zl)
 zmax	=zmax > max(zl)
 d0	=min(sqrt((xl-xf0[k])^2+(yl-yf0[k])^2),is)
 l1	=0.
 for i1=is-1,0,-1 do begin
  l1=l1+sqrt((xl(i1)-xl(i1+1))^2+(yl(i1)-yl(i1+1))^2)
  if (l1 ge len3[k]) then goto,field_start2
 endfor
FIELD_START2:
 i1	=i1>0
 i1_(k)=i1 
 l2	=0.
 for i2=is+1,np-1 do begin
  l2=l2+sqrt((xl(i2)-xl(i2-1))^2+(yl(i2)-yl(i2-1))^2)
  if (l2 ge len3[k]) then goto,field_end2
 endfor
FIELD_END2:
 i2	=i2<(np-1)
 i2_(k)=i2
 oplot,xl,zl,color=red,thick=th
;oplot,xl(i1:i2),zl(i1:i2),color=purple,thick=1
endfor

;______________________DISPLAY ROTATED LOOPS____________________________
for k_sel=0,nf_sel-1 do begin
 k	=isort(k_sel)
 ns	=ns_loop[k]
 xloop	=field_loop(0:ns-1,0,k)
 yloop	=field_loop(0:ns-1,1,k)
 zloop	=field_loop(0:ns-1,2,k)
 oplot,xloop,zloop,color=yellow,thick=2
 is	=0
 xf0[k]=xloop[is]
 yf0[k]=yloop[is]
 zf0[k]=zloop[is]
endfor

;______________DISPLAY MISALIGNMENT ANGLE HISTOGRAM_______________
 xh1	=0
 xh2	=45
 loop_per_bin=1.
 dh_bin	=(long((xh2-xh1)*loop_per_bin/nf) < 5.) >2. 	
 nh_bin	=long((xh2-xh1)/dh_bin)
 xh	=dh_bin*(findgen(nh_bin)+0.5)
 yh 	=histogram(dev_deg,min=xh1,max=xh2,binsize=dh_bin)
 xmed	=median(dev_deg)
 yh2 	=histogram(dev_deg2,min=xh1,max=xh2,binsize=dh_bin)
 xmed2	=median(dev_deg2)
 
 x1_     =0.74
 x2_     =0.98
 y1_     =0.75
 y2_     =0.92
 !p.position=[x1_,y1_,x2_,y2_]
 !p.title=' '
 !x.title=' '
 !x.title='Misalignment angle !9m!3 (deg)'
 !x.range=[xh1,xh2]
 !y.range=[0,max(yh2)*1.2]
 !y.style=1
 plot,xh,yh,psym=10,thick=1
 oplot,xh,yh2,psym=10,thick=3
 ymin	=0
 histo_step,xh,yh2,xpoly,ypoly,ymin
 polyfill,xpoly,ypoly,spacing=0.1,orientation=45
 oplot,xmed*[1,1],!y.range,linestyle=1
 oplot,xmed2*[1,1],!y.range,linestyle=2
 xyouts_norm,0.5,0.85,'!9m!3!D2!N='+string(xmed2,'(f4.1)')+'!U0!N',char*1.4
 xyouts_norm,0.5,0.70,'!9m!3!D3!N='+string(xmed,'(f4.1)')+'!U0!N',char*1.4

;_________________________DISPLAY PARAMETERS_____________________
!p.position=[0.,0.,1.,1.]
!p.title=''
!x.title=''
!y.title=''
!x.range=[0,1]
!y.range=[0,1]
!x.style=5
!y.style=5
plot,[0,0],[0,0]
xyouts, 0.1,0.94,filename+', time step='+strtrim(string(it,'(I4)'),2)+$
		 ', Catalog='+catfile,size=char*1.3

xyouts, 0.68,0.68,'INPUT:',size=char
xyouts, 0.70,0.66,'Instrument',size=char
xyouts, 0.82,0.66,'='+instr,size=char
xyouts, 0.70,0.64,'NOAA',size=char
xyouts, 0.82,0.64,'='+noaa+', '+helpos,size=char
xyouts, 0.70,0.62,'Wavelengths',size=char
wave_str=wave_[0] &for i=1,nwave-1 do wave_str=wave_str+','+wave_[i]
xyouts, 0.82,0.62,'='+wave_str+' A',size=char
xyouts, 0.70,0.60,'fov, hmax',size=char
xyouts, 0.82,0.60,'='+string(fov0,'(f5.3)')+', '+string(hmax,'(F4.2)'),size=char
xyouts, 0.70,0.58,'amis',size=char
xyouts, 0.82,0.58,'='+string(amis,'(F6.1)'),size=char
xyouts, 0.70,0.56,'nsm1',size=char
xyouts, 0.82,0.56,'='+string(nsm1,'(I4)'),size=char
xyouts, 0.70,0.54,'nmag_p,np',size=char
xyouts, 0.82,0.54,'='+string(nmag_p,'(I4)')+', '+string(nmag_np,'(I4)'),size=char
xyouts, 0.70,0.52,'nstruc,ngap',size=char
xyouts, 0.82,0.52,'='+string(nstruc,'(I4)')+', '+string(ngap,'(I4)'),size=char
xyouts, 0.70,0.50,'nitmin,nitmax',size=char
xyouts, 0.82,0.50,'='+string(nitmin,'(I4)')+', '+string(nitmax,'(I4)'),size=char
xyouts, 0.70,0.48,'prox_min',size=char
xyouts, 0.82,0.48,'='+string(prox_min,'(F4.1)'),size=char
xyouts, 0.70,0.46,'lmin,rmin',size=char
xyouts, 0.82,0.46,'='+string(lmin_wid,'(I4)')+', '+string(rmin_wid,'(I4)'),size=char
xyouts, 0.70,0.44,'qthresh1,2',size=char
xyouts, 0.82,0.44,'='+string(qthresh1,'(F5.2)')+', '+string(qthresh2,'(F5.2)'),size=char

xyouts, 0.68,0.38,'OUTPUT:',size=char
xyouts, 0.70,0.36,'[x1,x2]',size=char
xyouts, 0.82,0.36,'='+string(x1_sun,'(F7.4)')+', '+string(x2_sun,'(f7.4)'),size=char
xyouts, 0.70,0.34,'[y1,y2]',size=char
xyouts, 0.82,0.34,'='+string(y1_sun,'(F7.4)')+', '+string(y2_sun,'(f7.4)'),size=char
xyouts, 0.70,0.32,'dpix HMI',size=char
xyouts, 0.82,0.32,'='+string(dpix_mag,'(f6.4)')+$
                 ', '+string(dpix_mag*960.,'(f5.2)')+' arcsec',size=char
xyouts, 0.70,0.30,'dpix '+instr,size=char
xyouts, 0.82,0.30,'='+string(dpix_euv,'(f6.4)')+$,size=char
                 ', '+string(dpix_euv*960.,'(f5.2)')+' arcsec',size=char
xyouts, 0.70,0.28,'nloop,ndet',size=char
xyouts, 0.82,0.28,'='+string(nf,'(I4)')+' / '+string(ndet,'(I4)'),size=char
qloop	=float(nf)/float(ndet)
xyouts, 0.70,0.26,'nloop/ndet',size=char
xyouts, 0.82,0.26,'='+string(qloop,'(F6.3)'),size=char
xyouts, 0.70,0.24,'misalign',size=char
xyouts, 0.82,0.24,'='+string(xmed2,'(f5.1)')+', '+string(xmed,'(f5.1)')+' deg',size=char
xyouts, 0.70,0.22,'div-free',size=char
xyouts, 0.82,0.22,'='+string(merit(0),'(e8.1)'),size=char
xyouts, 0.70,0.20,'force-free',size=char
xyouts, 0.82,0.20,'='+string(merit(1),'(e8.1)'),size=char
xyouts, 0.70,0.18,'weight curr',size=char
xyouts, 0.82,0.18,'='+string(merit(2),'(e8.1)'),size=char
xyouts, 0.70,0.16,'qe_rebin',size=char
xyouts, 0.82,0.16,'='+string(qe_rebin,'(f6.3)'),size=char
xyouts, 0.70,0.14,'qe_model',size=char
xyouts, 0.82,0.14,'='+string(qe_model,'(f6.3)'),size=char
xyouts, 0.70,0.12,'qiso_corr',size=char
xyouts, 0.82,0.12,'='+string(qiso_corr,'(f6.3)'),size=char
xyouts, 0.70,0.10,'E_P',size=char
xyouts, 0.82,0.10,'='+string(e_p/1.e30,'(F7.1)')+' x 10!U30!N erg',size=char
xyouts, 0.70,0.08,'E_free',size=char
xyouts, 0.82,0.08,'='+string((e_np-e_p)/1.e30,'(F6.1)')+' x 10!U30!N erg',size=char
xyouts, 0.70,0.06,'E_NP/E_P',size=char
xyouts, 0.82,0.06,'='+string(e_np/e_p,'(f6.3)'),size=char
xyouts, 0.70,0.04,'Iterations',size=char
xyouts, 0.82,0.04,'='+string(iter,'(I3)'),size=char
xyouts, 0.70,0.02,'CPU',size=char
xyouts, 0.82,0.02,'='+string(cpu,'(f8.1)')+' s',size=char
endfor	;io=0 screen, io=3 postscript file

fig_close,io,fignr,ref
end