pro aia_mag_disp0,fileset,fov,hmax,ds,io,model,code
;+
; Project     : AIA/SDO
;
; Name        : AIA_LOOP_MAGFIT
;
; Category    : Data analysis   
;
; Explanation : Fitting of potential magnetic field model parameterized with
;		unipolar charges to observed loop projections
;
; Syntax      : IDL>aia_loop_magfit,fileset,fov,rsun
;
; Inputs      : fileset	 = beginning of filename 
;	        fov[4]   = [x1,y1,x2,y2] field-of-view in solar radii
;
; Outputs     : magnfile = x,y coordinates of magnetic field lines 
;
; 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]
aspect	 =(x2-x1)/(y2-y1)

;__________________READ MAGNETIC FIELD MODEL__________________________
coeff_file=file_search(fileset+'_coeff.dat')
readcol,coeff_file(0),fmag,xmag,ymag,zmag
ncomp   =n_elements(fmag)

;__________________READ MAGNETIC FIELD LINES__________________________
fieldfile=fileset+'_loop_'+model+code+'.sav'
restore,fieldfile	;-->field_,dev_loop
print,'Read field lines from file : ',fieldfile
if (model eq 'p') then field_=field_p

;____________________READ LOOPFILE_________________________________
loopfile  =fileset+'_loop.dat'		
readcol,loopfile,iloop,x,y,z,s
nloop     =long(max(iloop)+1)
print,nloop,' loop structures'
z1	=min(z)
z2	=max(z)

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

fig_open,io,form,char,fignr,plotname,unit
loadct,0

;_____________________DISPLAY MAGNETOGRAM______________________________
x1_	=0.5*(1.-0.75*aspect)
x2_	=0.5*(1.+0.75*aspect)
y1_	=0.03
y2_	=0.60
!p.title=' '
!p.position=[x1_,y1_,x2_,y2_]  
!x.range=[x1,x2]
!y.range=[y1,y2]
!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(image_mag,nxw,nyw)
endif
if (io ne 0) then image=image_mag
plot,[0,0],[0,0]
statistic,image_mag,zavg,zsig
c1	=zavg-6*zsig
c2	=zavg+2*zsig
tv,bytscl(image,min=c1,max=c2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal
coordsphere,1,0,0,0,0,dlat
xyouts_norm,1.02,0.,'CPU='+string(cpu,'(f7.1)')+' s',char,90,128
!noeras=1

loadct,5
blue    =50     	
red	=110		
orange	=160
yellow	=200

for il=0,nloop-1 do begin
 ind    =where(il eq iloop,ns)
 ind2	=where(field_(*,il,0) ne 0,n2)
 xl	=x[ind]
 yl	=y[ind]
 zl	=z[ind]
 oplot,xl,yl,thick=1,color=blue
 if (n2 ge 1) then begin
  xf	=field_[ind2,il,0]
  yf	=field_[ind2,il,1]
  zf	=field_[ind2,il,2]
  oplot,xl,yl,thick=5,color=blue
; oplot,xf,yf,color=red,linestyle=2
  d1	=min((xf-xl(0))^2+(yf-yl(0))^2+(zf-zl(0))^2,i1)
  d2	=min((xf-xl(ns-1))^2+(yf-yl(ns-1))^2+(zf-zl(ns-1))^2,i2)
  if (i2 lt i1) then begin &i0=i1 &i1=i2 &i2=i0 &endif
  oplot,xf(i1:i2),yf(i1:i2),color=red,thick=3
 endif
endfor
;oplot,curr(*,0),curr(*,1),color=yellow,thick=10,linestyle=2

for j=0,ncomp-1 do begin
 if (j le 10) then begin
; oplot,xmag(j)*[1,1],ymag(j)*[1,1],psym=1,symsize=1,color=orange
; xyouts,xmag(j)+0.001,ymag(j),string(j,'(I2)'),size=char*0.7,color=orange
 endif
endfor

;________________________DISPLAY VERTICAL PROJECTION______________________
y1_	=0.63
y2_	=0.83
!p.position=[x1_,y1_,x2_,y2_]  
!p.title=' '
!y.range=(z1+z2)/2.+0.2*(x2-x1)*[-1,1]	
!y.title='z-axis'
phi	=2.*!pi*findgen(361)/360.
plot,cos(phi),sin(phi)
coordsphere,1,0,0,-90,0,dlat
is	=0
for il=0,nloop-1 do begin
 ind    =where(il eq iloop,ns)
 ind2	=where(field_(*,il,0) ne 0,n2)
 xl	=x[ind]
 yl	=y[ind]
 zl	=z[ind]
 oplot,xl,zl,thick=1,color=blue
 if (n2 ge 1) then begin
  xf	=field_[ind2,il,0]
  yf	=field_[ind2,il,1]
  zf	=field_[ind2,il,2]
  oplot,xl,zl,thick=5,color=blue
; oplot,xf,zf,color=red,linestyle=2
  d1	=min((xf-xl(0))^2+(yf-yl(0))^2+(zf-zl(0))^2,i1)
  d2	=min((xf-xl(ns-1))^2+(yf-yl(ns-1))^2+(zf-zl(ns-1))^2,i2)
  if (i2 lt i1) then begin &i0=i1 &i1=i2 &i2=i0 &endif
  oplot,xf(i1:i2),zf(i1:i2),color=red,thick=3
 endif
endfor

;_________________________DISPLAY MISALIGNMENT HISTOGRAM____________________
dh	=2.
nh	=long(90/dh)
h	=histogram(dev_loop,min=0,max=90,binsize=dh)
xh	=dh*(0.5+findgen(nh))
dev	=median(abs(dev_loop))

y1_	=0.86
y2_	=0.96
!p.position=[x1_,y1_,x2_,y2_]  
if (model eq 'p') then label1='Potential field'
if (model eq 'n') then label1='Non-potential constant !9a!3 '
if (code  eq 'h') then label2='Helical loop code'
if (code  eq 'c') then label2='NLFFF cube code'
!p.title=fileset+' - '+label1+' - '+label2 
!x.range=[0,90]
!y.range=[0,max(h)+1]
plot, xh,h,psym=10
oplot,dev*[1,1],!y.range,thick=2,linestyle=2
histo_step,xh,h,xpoly,ypoly
polyfill,xpoly,ypoly,spacing=0.1,orientation=45

xyouts_norm,0.7,0.85,'N!Dmag!N='+string(ncomp,'(I5)'),char
xyouts_norm,0.7,0.70,'N!Dloop!N='+string(nloop,'(I5)'),char
xyouts_norm,0.7,0.55,'h!Dmax!N='+string(hmax,'(f5.2)'),char
xyouts_norm,0.7,0.40,'ds ='+string(ds,'(f8.4)'),char
xyouts_norm,0.7,0.25,'!9a!3!Dmis!N='+string(dev,'(f5.2)')+' deg',char
fig_close,io,fignr,ref
end