pro euvi_hydro_fit_disp,hydrofile,ct
;+ 
; Project     : STEREO
;
; Name        : EUVI_HYDRO_FIT_DISP 
;
; Category    : Hydrodynamic modeling
;
; Explanation : This routine produces a fit of a 3D model of an active region
;		by filling loops according to a hydrodynamic solution.
;
; Syntax      : IDL> euvi_hydro_fit,fieldfile,image_pair,...
;
; Inputs      : fieldfile = datafile containing (x,y,z)-coordinates of loops (e.g.,'loop_A.dat')
;		para	= structure containing image parameters
;		ct	=5   ;IDL colortable
;
; Outputs     : postscript file of plot *.ps
;
; History     : Oct 2007, Version 1 written by Markus J. Aschwanden
;               Jun 2008, Version 2 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;______________________DATA + MODEL IMAGE PLOT________________________________
restore,hydrofile
ipos    =strpos(hydrofile,'hydro.sav')
filename =strmid(hydrofile,0,ipos)
dateobs =strmid(hydrofile,ipos-9,8)
wave_   =['171','195','284']
i1      =fov(0)
j1      =fov(1)
i2      =fov(2)
j2      =fov(3)
nx      =(i2-i1+1)
ny      =(j2-j1+1)
ix	=i1+findgen(nx)
iy	=j1+findgen(ny)
ii	=where((ix ge i1+nk) and (ix le i1+nx-1-nk))
jj	=where((iy ge j1+nk) and (iy le j1+ny-1-nk))
ixk	=ix(ii)
iyk	=iy(jj)
image2	=image0-image1
perc	=fltarr(6)

loadct,ct
sc_	=['a','b']
sc__	=['A','B']
for io=0,3,3 do begin
for isc=0,1 do begin
 fignr	  ='hydro_fit'+sc_(isc)
 plotname=filename
 form    =1     ;landscape format
 char    =0.5   ;smaller character size
 unit    =isc     
 ref     =''     

 fig_open,io,form,char,fignr,plotname,unit
 for iwave=0,2 do begin
  iw	 =iwave+isc*3
  for ipanel=0,2 do begin
   if (ipanel eq 0) then image=reform(image0(nk:nx-1-nk,nk:ny-1-nk,iw))
   if (ipanel eq 1) then image=reform(image1(nk:nx-1-nk,nk:ny-1-nk,iw))
   if (ipanel eq 2) then image=reform(image2(nk:nx-1-nk,nk:ny-1-nk,iw))
   zmax	=max(image0(*,*,iw))
   x1_  =0.10+0.3*ipanel 
   x2_  =x1_+0.28 
   y1_	=0.7-0.25*iwave
   y2_	=y1_+0.23
   if (io eq 0) then begin
    nxw  =long(!d.x_vsize*(x2_-x1_)+0.5)
    nyw  =long(!d.y_vsize*(y2_-y1_)+0.5)
    ima  =congrid(image,nxw,nyw)
   endif
   if (io ne 0) then ima=image
   if (ipanel eq 0) then begin
    statistic,image,zav,zsig
    bw      =300.   ;100% black/white contrast
    z1      =(zav-zsig*(bw/100.))>0.
    z2      = zav+zsig*(bw/100.)
    zmax    =max(image)
    total0  =total(image)
   endif
   if (ipanel eq 1) then begin
    total1  =total(image)
    perc(iw)=100.*total1/total0
   endif
   inv     =-1.
   if (inv eq +1.) then zz=ima
   if (inv eq -1.) then zz=z2-ima

   !p.position=[x1_,y1_,x2_,y2_]
   !x.range=[i1+nk,i2-nk+0.1]
   !y.range=[j1+nk,j2-nk+0.1]
   !x.title=' '
   !y.title=' '
   !x.style=1
   !y.style=1
   plot,[0,0],[0,0]
   tv,bytscl(zz,min=z1,max=z2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal
   oplot,[i1,i2,i2,i1,i1],[j1,j1,j2,j2,j1]
   !noeras=1
   if (ipanel eq 0) then xyouts_norm,0.05,0.95,dateobs+': '+wave_(iwave)+$
     ' A, STEREO/'+sc__(isc),1,0,0 
   if (ipanel eq 1) then xyouts_norm,0.05,0.95,'Model: '+wave_(iwave)+$
	' A, '+string(perc(iw),'(I3)')+' %',1,0,0 
   !noeras=1
  endfor
 endfor
 if (isc eq 0) then pa=avg(perc(0:2))
 if (isc eq 1) then pa=avg(perc(3:5))

 !p.position=[0,0,1,1]
 !x.style=5
 !y.style=5
 plot,[0,0],[0,0]
 xyouts_norm,0.1,0.15,hydrofile
 xyouts_norm,0.1,0.13,'T!Dmin!N,T!Dmax!N = '+string(temin,'(f4.1)')+', '+string(temax,'(f4.1)')+' MK',1
 xyouts_norm,0.1,0.11,'EM!Dmin!N = '+string(em44_min,'(f8.4)')+' 10!U44!N cm!U-3!N'
 xyouts_norm,0.1,0.09,'w = '+string(w_km,'(i8)')+' km',1
 xyouts_norm,0.5,0.14,'Retrieved flux = '+string(pa,'(f5.1)')+' %',1.5
 xyouts_norm,0.5,0.11,'N!Dloop!N  = '+string(nloopmodel,'(i8)'),1.5

 if (io ge 1) then fig_close,io,fignr,ref
endfor
endfor

;________________________CORRELATION PLOT______________________________________
window,2,xsize=512,ysize=512
clearplot
loadct,5
range   =[1,max(image1)]
!x.range=range
!y.range=range
plot_oo,range,range
for isc=0,1 do begin
for iwave=0,2 do begin
 iw	=iwave+isc*3
 image  =reform(image0(*,*,iw))
 model  =reform(image1(*,*,iw))
 print,'STATISTIC MODEL/IMAGE'
 statistic,(model>1)/(image>1)
 col    =50*(1+iwave)
 oplot,image,model,psym=3,color=col
 print,wave_(iwave)+' A: flux retrieved = ',perc(iw),' %',format='(a,f6.1,a)'
endfor
endfor
print,'Average flux retrieved = ',avg(perc),' %',format='(a,f6.1,a)'
end