pro magn_disp3,modelfile,runfile,para,itype,inr,io,error

;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : MAGN_DISP3
;
; Category    : display of magnetic data modeling
;
; Explanation : displays photospheric maps of magnetic field components
;		Bx(x,y), By(x,y), Bz(x,y) and force-free parameter alpha(x,y)
;
; Syntax      : IDL>magn_disp3,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
;		error	  = display error bars (variation of alpha)
;		regress	  = display linear regression fit
;
; Outputs     ; postscript file <runfile>_alpha.ps
;
; History     : 16-May-2011, Version 1 written by Markus J. Aschwanden
;               30-Jan-2012, Version 2 released to SSW
;
; Contact     : aschwanden@lmsal.com
;-

;___________________INPUT PARAMETERS_______________________________
ds      =para[0]
thresh  =para[1]
nmag    =para[2]
nsm     =para[3]
qloop   =para[4]
nseg    =para[5]
acc     =para[6]
amax    =para[7]
scaling =1.0		;0.4	empirical scaling factor for alpha theoretical

;____________________READ DATA_____________________________________
savefile1=modelfile+'_loops.sav'
restore,savefile1 ;-->bzmap,x,y,field_lines,alpha,q_scale,angle,cpu,merit
field_lines0=field_lines
a_loop	=alpha[*,0]*scaling
s_loop	=alpha[*,1]

savefile2=runfile+'_fit.sav'
restore,savefile2 ;-->bzmap,x,y,field_lines,alpha,q_scale,angle,cpu,merit
dim	=size(alpha)
nf	=dim(1)
 a_fit	=alpha[*,0]		;cube-calculated alpha
 s_fit	=alpha[*,1]
;a_fit	=alpha[*,2]		;theoretical alpha
;s_fit	=alpha[*,3]
y2	=max([abs(a_fit),abs(a_loop)])
y2	=(y2 > amax*0.1) 

;____________________READ DATA CUBES________________________________
savefile=modelfile+'_cube.sav'
restore,savefile 	;-->x,y,z,bx,by,bz,b,a
nx      =n_elements(x)
ny      =n_elements(y)
nz      =n_elements(z)
map_model=reform(a(1:nx-2,1:ny-2,1))*scaling

;____________________READ DATA CUBES________________________________
savefile=runfile+'_cube.sav'
restore,savefile	;-->x,y,z,bx,by,bz,b,a
nx      =n_elements(x)
ny      =n_elements(y)
nz      =n_elements(z)
map_fit =reform( a(1:nx-2,1:ny-2,1))

;____________________DISPLAY________________________________________
form    =1      ;0=landscape, 1=portrait
char    =1.0    ;character size
fignr	=''
ref     =''
dlat    =5
unit	=0
sym	=[4,1]
plotname=runfile+'_alpha'
fig_open,io,form,char,fignr,plotname,unit
loadct,0
for iplot=0,1 do begin
 if (iplot eq 0) then begin &xx=a_loop    &yy=a_fit   &endif
 if (iplot eq 1) then begin &xx=map_model &yy=map_fit &endif
 x1_     =0.15
 x2_     =0.85
 y1_     =0.55-0.5*iplot
 y2_     =0.95-0.5*iplot
 !p.title=runfile
 !p.position=[x1_,y1_,x2_,y2_]
 !x.title='Force-free alpha (loop)'
 !y.title='Force-free alpha (fit)'
 !x.range=[-1,1]*y2
 !y.range=[-1,1]*y2
 !x.style=0
 !y.style=0
 plot,xx,yy,psym=sym[iplot],symsize=1 
 n	=n_elements(xx)<n_elements(yy)
 if (iplot eq 0) and (error eq 1) then $
  for j=0,n-1 do oplot,a_loop[j]*[1,1],a_fit[j]+[-1,1]*s_fit[j]
 xyouts_norm,0.1,0.9,'N='+string(n,'(i4)'),char*2
 oplot,!x.crange,!x.crange
 oplot,!x.crange,[0,0],linestyle=1
 oplot,[0,0],!y.crange,linestyle=1
 !noeras=1

 qav	=0.
 qsig   =0.
 if (itype eq 0) or (abs(inr) ge 7) then begin
  text  ='!9a!3!Dfit!N!3/!9a!3!Dloop!N='
  q	=yy/xx
  statistic,q,qav,qsig	
  statistic,q	
  xfit=[-1,1]*!x.crange[1]
  yfit=xfit*qav
  oplot,xfit,yfit,linestyle=2
 endif
 if (itype eq 1) and (abs(inr) le 6) then begin
  text  ='!9a!3!Dfit!N!3='
  q	=yy
  statistic,q,qav,qsig	
  statistic,q	
 endif
 xyouts_norm,0.1,0.8,text+string(qav,'(f6.2)')+'+'+string(qsig,'(f6.2)'),char*2
endfor
fig_close,io,fignr,ref
end