pro aia_mag_disp2,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
;		np	 = number of loop points	
;
; Outputs:      *_field.sav = savefile containing fieldlines
;                             xfield_,yfield_,zfield_,dev_
;
; History     : 12-May-2011, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;__________________READ LOOP COORDINATES______________________________
loopfile  =fileset+'*_loop.dat'
readcol,loopfile,iloop,x,y,z,s
nloop   =long(max(iloop)+1)
print,'NUMBER OF LOOPS = ',nloop

;__________________READ LOOP PARAMETERS______________________________
fieldfile=fileset+'_loop_'+model+code+'.sav'
restore,fieldfile       ;-->field_,dev_loop,rad_loop,ntwist_loop,alpha_loop

;__________________READ FIELD LINES___________________________________
fieldfile=fileset+'_field_'+model+code+'.sav'
restore,fieldfile       ;-->field_,dev_loop,alpha_cgs,jcurr_cgs

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

fig_open,io,form,char,fignr,plotname,unit
loadct,0
xtitle	=['Curvature radius','Number ot twist N!Dtwist!N',$
          'Nonlinear force-free !9a!3','Current density j']
unit_	=['Mm','',' 10!U-11!N cm','Mx cm!U-2!N s!U-1!N']
for j=0,3 do begin 
for i=0,1 do begin
 x1_	=0.1+0.45*i 
 x2_	=x1_+0.3 
 y2_	=0.95-0.225*j 
 y1_	=y2_-0.175
 !p.position=[x1_,y1_,x2_,y2_]  
 if (j eq 0) then y=abs(rad_loop)*696.
 if (j eq 1) then y=abs(ntwist_loop)
 if (j eq 2) then y=abs(alpha_cgs)*1.0e11	
 if (j eq 3) then y=abs(jcurr_cgs)
 ymed	=median(y)
 y1	=min(y) &if (i eq 1) then y1=y1>0.1
 y2	=max(y)
 if (i eq 0) then y2=ymed*4
 nh	=15

 if (i eq 0) then begin
  bin	=(y2-y1)/float(nh)
  x_ 	=y1+(y2-y1)*(findgen(nh)+0.5)/float(nh)
  y_	=histogram(y,min=y1,max=y2,bin=bin)
 endif
 if (i eq 1) then begin
  ystepmin=1.
  histo_log,y,y1,y2,nh,x_,dx,y_,yerr,xterp,ystep,ystepmin
 endif
 histo_step,x_,y_,xh,yh		&if (i eq 1) then yh=yh > 1.e-2

 !p.title=''
 if (j eq 0) then !p.title=fileset
 !x.title=xtitle(j)+' ['+unit_(j)+']'
 !y.title='Number of loops'
 !x.range=[y1,y2]
 !y.range=minmax(yh) &if (i eq 1) then !y.range=minmax(yh)
 !x.style=0
 !y.style=0

 if (i eq 0) then plot,xh,yh,psym=10,thick=3
 if (i eq 1) then plot_oo,xh,yh,psym=10,thick=3
 ymed	=median(y)
 print,xtitle(j),' median=',ymed
 if (i eq 0) then $
  xyouts_norm,0.4-0.07*j,0.8,'Median = '+string(ymed,'(f7.2)')+' '+unit_[j],char
 oplot,ymed*[1,1],!y.range,linestyle=2
 !noeras=1
endfor
endfor
fig_close,io,fignr,ref
end