pro magn_disp1,modelfile,runfile,para,io
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : MAGN_DISP1
;
; 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>magn_disp1,modelfile,runfile,para,io
;
; Inputs      : modelfile = input filename *_loops.sav
;               runfile   = input filename *_coeff.dat, *_fit.sav
;               para      = input parameters
;                          [ds,ngrid,nmag,nsm,qloop,nseg,acc,amax,hmax,dpix,
;                           niter,qzone,eps,qv,thresh,q_mag,halt,iopt]
;		io	  = ;0=X-window, 1=ps, 2=eps, 3=color ps
;
; Outputs     ; postscript file <runfile>.ps
;
; History     : 16-May-2011, Version 1 written by Markus J. Aschwanden
;               30-Jan-2012, Version 2 released to SSW
;
; Contact     : aschwanden@lmsal.com
;-

ds      =para[0]
ngrid   =para[1]
nmag    =para[2]
nsm     =para[3]
qloop   =para[4]
nseg    =para[5]
acc     =para[6]
amax    =para[7]
hmax    =para[8]
dpix    =para[9]
niter   =para[10]
qzone   =para[11]
eps     =para[12]
qv      =para[13]
thresh  =para[14]
q_mag	=para[15]
halt	=para[16]
iopt	=para[17]

;____________________READ DATA_____________________________________
savefile1=modelfile+'_loops.sav'	
restore,savefile1 ;-->bzmap,x,y,field_lines,alpha,q_scale,angle,cpu,merit
field_lines0=field_lines
alpha0	=alpha
alpha0_med=median(alpha0(*,0))
help,field_lines0

savefile2=modelfile+'_field.sav'		;-->*_field.sav
restore,savefile2 ;-->bzmap,x,y,field_lines,alpha,q_scale,angle,cpu,merit,dev_loop
savefile2=runfile+'_fit.sav'			;-->*_fit.sav
restore,savefile2 ;-->bzmap,x,y,field_lines,alpha,q_scale,angle,cpu,merit,dev_loop
help,field_lines

coeff_file=runfile+'_coeff.dat'
readcol,coeff_file,m1,x1,y1,z1,a1
imag    =n_elements(m1)
coeff	=[[m1],[x1],[y1],[z1]]
aia_mag_map,coeff,1.,x,y,bzmodel

;____________________DISPLAY________________________________________
form    =1      ;0=landscape, 1=portrait
char    =0.8    ;character size
fignr	=''
ref     =''
dlat    =5
unit	=0
fig_open,io,form,char,fignr,runfile,unit
loadct,0

;_____________________DISPLAY MAGNETOGRAM______________________________
nx	=n_elements(x)
ny	=n_elements(y)
x1_     =0.05
x2_     =0.75
y1_     =0.25
y2_     =0.75
!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='y-axis'
!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(bzmodel,nxw,nyw)		;bzmap
endif
if (io ne 0) then image=bzmodel 		;bzmap
plot,[0,0],[0,0]
c2      =max(abs(bzmodel))			;bzmap
c1      =-c2*2.
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

;______________________DISPLAY FIELD LINES_____________________________
loadct,5
red	=110
blue	=50
dim	=size(field_lines0)
if (dim[0] le 2) then nf0=dim[2]
if (dim[0] ge 3) then nf0=dim[3]
xf0	=fltarr(nf0)
yf0	=fltarr(nf0)
zf0	=fltarr(nf0)
nf0	=dim(3)
if (dim(0) eq 2) then nf0=1
len	=fltarr(nf0)
for k=0,nf0-1 do begin
 ind	=where(field_lines0(*,2,k) ne 0,ns) 
 if (ns ge 1) then begin 
  xl	=field_lines0(0:ns-1,0,k)
  yl	=field_lines0(0:ns-1,1,k)
  zl	=field_lines0(0:ns-1,2,k)
  is	=long(qloop*(ns-1))
  xf0[k]=xl[is]
  yf0[k]=yl[is]
  zf0[k]=zl[is]
  ds_len=sqrt((xl(is)-xl(is-1))^2+(yl(is)-yl(is-1))^2+(zl(is)-zl(is-1))^2)
  for is=1,ns-1 do len(k)=len(k)+ds_len
  oplot,xl,yl,color=blue,thick=3
 endif
endfor
oplot,xf0,yf0,psym=4

;_______________________DISPLAY FITTED LINES____________________________
dim	=size(field_lines)
nf	=dim(3)
if (dim(0) eq 2) then nf=1
istart	=lonarr(nf0)
iend	=lonarr(nf0)
 for k=0,nf-1 do begin
 ind	=where(field_lines(*,2,k) ne 0,ns) 
 if (ns ge 1) then begin 
  xl	=field_lines(0:ns-1,0,k)
  yl	=field_lines(0:ns-1,1,k)
  zl	=field_lines(0:ns-1,2,k)
  dist  =min((xl-xf0(k))^2+(yl-yf0(k))^2+(zl-zf0(k))^2,imid)
  nhalf	=long(len(k)/(2.*ds))
  istart(k)=((imid-nhalf)>0)<(ns-1)
  iend(k)  =((imid+nhalf)>0)<(ns-1)
  oplot,xl,yl,color=red,thick=1,linestyle=2
  oplot,xl(istart(k):iend(k)),yl(istart(k):iend(k)),color=red,thick=2
 endif
endfor
zmid	=median(field_lines(0:10,2,*))
zrange	=(max(x)-min(x))*0.55		;aspect ratio yrange/xrange=0.55

;________________________DISPLAY VERTICAL PROJECTION______________________
y1_     =0.75
y2_     =0.96
!p.position=[x1_,y1_,x2_,y2_]
ip	=strpos(runfile,'MDI')
np	=strlen(runfile)
!p.title=strmid(runfile,ip,np-ip+1)+' : div-free='+string(merit[0],'(f8.6)')+','+$
                   ' force-free='+string(merit[1],'(f8.6)')
!y.range=zmid+(zrange/2.)*[-1,1]
!y.title=' '
phi     =2.*!pi*findgen(361)/360.
plot,cos(phi),sin(phi)
coordsphere,1,0,0,-90,0,dlat
is      =0
for k=0,nf0-1 do begin
 ind	=where(field_lines0(*,2,k) ne 0,ns) 
 if (ns ge 1) then begin 
  xl	=field_lines0(0:ns-1,0,k)
  yl	=field_lines0(0:ns-1,1,k)
  zl	=field_lines0(0:ns-1,2,k)
  oplot,xl,zl,color=blue,thick=3
 endif
 ind	=where(field_lines(*,2,k) ne 0,ns) 
 if (ns ge 1) then begin 
  xl	=field_lines(0:ns-1,0,k)
  yl	=field_lines(0:ns-1,1,k)
  zl	=field_lines(0:ns-1,2,k)
  oplot,xl,zl,color=red,thick=1,linestyle=2
  oplot,xl(istart(k):iend(k)),zl(istart(k):iend(k)),color=red,thick=2
 endif
endfor
oplot,xf0,zf0,psym=4

;________________________DISPLAY Z-Y PROJECTION______________________
x1_     =0.75
x2_     =0.98
y1_     =0.25
y2_     =0.75
!p.position=[x1_,y1_,x2_,y2_]
!p.title=' '
!x.range=zmid+(zrange/2.)*[-1,1]
!y.range=[y[0],y[ny-1]]
phi     =2.*!pi*findgen(361)/360.
plot,cos(phi),sin(phi)
coordsphere,1,0,0,-90,0,dlat
is      =0
for k=0,nf0-1 do begin
 ind	=where(field_lines0(*,2,k) ne 0,ns) 
 if (ns ge 1) then begin 
  xl	=field_lines0(0:ns-1,0,k)
  yl	=field_lines0(0:ns-1,1,k)
  zl	=field_lines0(0:ns-1,2,k)
  oplot,zl,yl,color=blue,thick=3
 endif
 ind	=where(field_lines(*,2,k) ne 0,ns) 
 if (ns ge 1) then begin 
  xl	=field_lines(0:ns-1,0,k)
  yl	=field_lines(0:ns-1,1,k)
  zl	=field_lines(0:ns-1,2,k)
  oplot,zl,yl,color=red,thick=1,linestyle=2
  oplot,zl(istart(k):iend(k)),yl(istart(k):iend(k)),color=red,thick=2
 endif
endfor
oplot,zf0,yf0,psym=4

;_________________________DISPLAY PARAMETERS_____________________
!p.position=[0.77,0.75,1.00,0.98]
plot,[0,0],[0,0],xrange=[0,1],yrange=[0,1],xstyle=5,ystyle=5
xyouts_norm,0.,0.99,'nfield = '+string(nf,'(I5)'),char
xyouts_norm,0.,0.92,'dpix   = '+string(dpix,'(f10.6)'),char
xyouts_norm,0.,0.85,'ds     = '+string(ds,'(f10.6)'),char
xyouts_norm,0.,0.78,'ngrid  = '+string(ngrid,'(i5)'),char
xyouts_norm,0.,0.71,'thresh = '+string(thresh,'(f6.1)')+' G',char
xyouts_norm,0.,0.64,'imag/nmag='+string(imag,'(I3)')+'/'+string(nmag,'(I3)'),char
xyouts_norm,0.,0.57,'nsm = '+string(nsm,'(I5)')+$
                  ', iopt = '+string(iopt,'(I1)'),char
xyouts_norm,0.,0.50,'qloop  = '+string(qloop,'(f5.2)'),char
xyouts_norm,0.,0.43,'nseg   = '+string(nseg,'(I5)'),char
xyouts_norm,0.,0.36,'acc    = '+string(acc,'(f8.4)'),char
xyouts_norm,0.,0.29,'amax   = '+string(amax,'(f8.1)'),char
xyouts_norm,0.,0.22,'iter/niter='+string(iter+1,'(i3)')+'/'+string(niter,'(i3)'),char
xyouts_norm,0.,0.15,'qzone  = '+string(qzone,'(f5.3)'),char
xyouts_norm,0.,0.08,'qv     = '+string(qv,'(f6.1)'),char
xyouts_norm,0.,0.01,'cpu    = '+string(cpu,'(f8.1)')+' s',char

;______________DISPLAY MISALIGNMENT ANGLE HISTOGRAM_______________
xh1	=0
xh2	=90
dh	=2.
nh	=long((xh2-xh1)/dh)
xh	=dh*(findgen(nh)+0.5)
yh	=histogram(dev_loop,min=xh1,max=xh2,binsize=dh)
xmed	=median(dev_loop)

dx2	=1
nx2	=long((xh2-xh1)/dx2)
xg      =xh1+dx2*(findgen(nx2))
h       =histogram(dev_loop,min=xh1,max=xh2,bin=dx)
h0      =min(yh(1:nh-1))
result  =gaussfit(xh(0:nh-1),yh(0:nh-1)-h0,a,nterms=3)
peak    =a(0)
x0      =a(1)
s0      =a(2)
z       =(xg-x0)/s0
yg      =peak*exp(-z^2/2.)+h0

x1_     =0.05
x2_     =0.75
y1_     =0.05
y2_     =0.22
!p.position=[x1_,y1_,x2_,y2_]
!p.title=' '
!y.title='Number of loops'
!x.title='Misalignment angle !9a!3!Dmis!N (deg)'
!x.range=[0,90]
!y.range=[0,max(yh)*1.1]
!y.style=1
plot,xh,yh,psym=10
ymin	=0
histo_step,xh,yh,xpoly,ypoly,ymin
polyfill,xpoly,ypoly,spacing=0.1,orientation=45
oplot,angle*[1,1],!y.range,thick=2
oplot,xmed*[1,1],!y.range,thick=2,linestyle=2
xyouts,angle,0.9*!y.crange[1],' A',size=char
xyouts,xmed,0.9*!y.crange[1],' M',size=char
oplot,xg,yg,thick=3
oplot,x0*[1,1],[0,h0+peak]
xyouts,x0,h0+peak,'P',size=char

xyouts_norm,1.02,0.9,'peak (P) ='+string(x0,'(f4.1)')+'+'+$
                                  string(s0,'(f4.1)'),char*1.4
xyouts_norm,1.02,0.8,'median (M) ='+string(xmed,'(f5.1)'),char*1.4
xyouts_norm,1.02,0.7,'average (A) ='+string(angle,'(f5.1)'),char*1.4
print,'Gauss peak=',x0,s0
print,'Median=',xmed
print,'Average=',angle

fig_close,io,fignr,ref
end