pro magn_disp4,modelfile,runfile,para,io
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : MAGN_DISP4
;
; Category    : display of magnetic data modeling
;
; Explanation : displays observed and fitted force-free parameters alpha
;
; Syntax      : IDL>magn_disp4,modelfile,runfile,para,io
;
; Inputs      : modelfile = model input filename *_cube.sav
;               runfile   = fit   input filename *_cube.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
;-

amax	=para[7]

;______________________READ DATA_____________________________________
cubefile=modelfile+'_cube.sav'
restore,cubefile        ;help,x,y,z,bx,by,bz,b,a
nx      =n_elements(x)
ny      =n_elements(y)
nz      =n_elements(z)
bx0	=reform(bx(1:nx-2,1:ny-2,1))
by0	=reform(by(1:nx-2,1:ny-2,1))
bz0	=reform(bz(1:nx-2,1:ny-2,1))
amap0	=reform( a(1:nx-2,1:ny-2,1))
xx	=x[1:nx-2]
yy	=y[1:ny-2]

savefile=runfile+'_fit.sav'
restore,savefile      ;-->bzmap,x,y,field_lines,q_scale,angle,cpu,merit

cubefile=runfile+'_cube.sav'
restore,cubefile        ;help,x,y,z,bx,by,bz,b,a
bx1	=reform(bx(1:nx-2,1:ny-2,1))
by1	=reform(by(1:nx-2,1:ny-2,1))
bz1	=reform(bz(1:nx-2,1:ny-2,1))
amap1	=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
plotname=runfile+'_maps'
fig_open,io,form,char,fignr,plotname,unit
label	=['B!Dx!N','B!Dy!N','B!Dz!N','!9a!3']
for j=0,3 do begin
 for i=0,1 do begin
  x1_     =0.1+0.4*i
  x2_     =x1_+0.32
  y2_     =0.96-0.23*j
  y1_     =y2_-0.20
  !p.position=[x1_,y1_,x2_,y2_]
  !p.title=' '
  if (j eq 0) and (i eq 0) then !p.title=modelfile
  if (j eq 0) and (i eq 1) then !p.title=runfile
  !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 (i eq 0) and (j eq 0) then image=bx0
  if (i eq 1) and (j eq 0) then image=bx1
  if (i eq 0) and (j eq 1) then image=by0
  if (i eq 1) and (j eq 1) then image=by1
  if (i eq 0) and (j eq 2) then image=bz0
  if (i eq 1) and (j eq 2) then image=bz1
  if (i eq 0) and (j eq 3) then image=amap0
  if (i eq 1) and (j eq 3) then image=amap1
  if (i eq 0) then image0=image
  if (i eq 1) then image1=image
  if (i eq 1) then q=total(abs(image1))/total(abs(image0))
  if (i eq 0) then begin
   c2	=max(abs(image))
   level1 =c2*(-1.+findgen(11)/float(10))
   level2 =c2*findgen(11)/float(10)
  endif
  plot,[0,0],[0,0]
  contour,image,xx,yy,level=level1,c_linestyle=2
  contour,image,xx,yy,level=level2,/overplot
  coordsphere,1,0,0,0,0,10
  if (i eq 0) then xyouts_norm,0.1,0.8,label(j),char*2
  if (i eq 1) then xyouts_norm,0.1,0.9,'|!9a!3| ratio: '+string(q,'(f6.2)'),char
  if (j eq 3) then begin
   if (i eq 0) then dyy=(yy(ny-3)-yy(ny/2-1))/max(image(*,ny/2))
   oplot,xx,dyy*image(*,ny/2),thick=3
  endif
  !noeras=1
 endfor
endfor
fig_close,io,fignr,ref
end