pro aia_mag_disp1,fileset,fov,hmax,ds,io,model,code
;+
; Project     : AIA/SDO
;
; Name        : AIA_MAG_LOOP
;
; Category    : SDO/AIA data modeling 
;
; Explanation : Fitting of twisted field lines to observed loop projections
;
; Syntax      : IDL>aia_mag_loop,fileset,fov,hmax,ntwist,np
;
; Inputs      : fileset	 = filename 
;	        fov[4]   = [x1,y1,x2,y2] field-of-view (in solar radii)
;		ds	 = step size 
;		hmax	 = height range of magnetic field model
;
; Outputs:      *_disp1.ps = plotfile
;
; History     : 12-May-2011, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;__________________READ MAGNETOGRAM___________________________________
aia_mag_fov,fileset,fov,image_mag,x,y,nx,ny,dpix
x1      =fov(0)
x2      =fov(2)
y1      =fov(1)
y2      =fov(3)
bgauss	 =max(abs(image_mag))	&print,'B[Gauss]=',bgauss
statistic,image_mag,zavg0,zsig0

;__________________READ MAGNETIC FIELD LINES__________________________
fieldfile=fileset+'_field_'+model+code+'.sav'
restore,fieldfile       ;-->field_,dev_loop
bx_diff=bx_np-bx_p
by_diff=by_np-by_p
bz_diff=bz_np-bz_p
dim=size(bz_diff)
ni	=dim(1)
nj	=dim(2)
bz_diff(*,nj-3:nj-1)=0.

;__________________DISPLAY____________________________________________
form    =1      ;0=landscape, 1=portrait
char    =0.8    ;character size
fignr	=''
plotname=fileset+'_disp1_'+model+code 	;filename
ref     =''     
unit    =0      ;window number
dlat	=5

fig_open,io,form,char,fignr,plotname,unit
loadct,0
for j=0,2 do begin
for i=0,2 do begin
 if (i eq 0) then begin
  if (j eq 0) then begin &ima=bz_p &!p.title='B!Dz!Upot!N' &endif
  if (j eq 1) then begin &ima=by_p &!p.title='B!Dy!Upot!N' &endif
  if (j eq 2) then begin &ima=bx_p &!p.title='B!Dx!Upot!N' &endif
 endif
 if (i eq 1) then begin
  if (j eq 0) then begin &ima=bz_diff &!p.title='B!Dz!UNP!N-B!Dz!Upot!N' &endif
  if (j eq 1) then begin &ima=by_diff &!p.title='B!Dy!UNP!N-B!Dy!Upot!N' &endif
  if (j eq 2) then begin &ima=bx_diff &!p.title='B!Dx!UNP!N-B!Dx!Upot!N' &endif
 endif
 if (i eq 2) then begin
  if (j eq 0) then begin &ima=image_mag &!p.title='B!Dz!N!Uobs!N' &endif
  if (j eq 1) then begin &ima=alpha &!p.title='alpha' &endif
  if (j eq 2) then begin &ima=jcurr &!p.title='j-current' &endif
 endif
 if (i ne 1) then dim=size(ima)
 ni	=dim(1)
 nj	=dim(2)
 ima	=rebin(ima(0:ni-1,0:nj-1),ni*2,nj*2)
 if (i ge 1) and (j ge 1) then ima=smooth(ima,3)
 x1_	=0.05+0.33*i
 x2_	=x1_+0.27
 y2_	=0.8-0.26*j
 y1_	=y2_-0.2
 !p.position=[x1_,y1_,x2_,y2_]  
 !x.range=[x1,x2]
 !y.range=[y1,y2]
 !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(ima,nxw,nyw)
 endif
 if (io ne 0) then image=ima
 plot,[0,0],[0,0]
 statistic,ima,zavg,zsig
 nsig=3.
 if (i ge 1) and (j ge 1) then begin
  nsig=0.1
  ind0	=where(ima ne 0)
  statistic,ima(ind0),zavg,zsig
 endif
 c1	=-nsig*zsig*2
 c2	=+nsig*zsig
 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

if (model eq 'p') then label1='Potential field'
if (model eq 'n') then label1='Non-potential field'
if (code  eq 'h') then label2='Helical loop code'
if (code  eq 'c') then label2='NLFFF cube code'
title=fileset+' - '+label1+' - '+label2
if (i eq 0) and (j eq 0) then xyouts_norm,0.,1.2,title,1.5

endfor
endfor

fig_close,io,fignr,ref
end