pro euvi_hydro_fit,fieldfile,fov,para,ct,hydrofile,test,index_a,index_b,nmulti,nloopmax,em44_min,em44_max,w_km,spike
;+ 
; Project     : STEREO
;
; Name        : EUVI_HYDRO_FIT 
;
; 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
;-

common powell_para,s_,hs,flux_obs,te_,resp_,flux_model,temp_model,area_h,par_model,temin,temax,hmax,rad_pix,frac_max,n1_h

cdelt1	=para.cdelt1
crpix1	=para.crpix1
crpix2	=para.crpix2
rsun   	=para.rsun
rad_pix	=rsun/cdelt1
file_exist=findfile(fieldfile,count=nfiles)
if (nfiles eq 0) then print,'No previous fieldfile found';
print,'Datafile with loop coordinates found  LOOPFILE=',fieldfile 
ipos    =strpos(fieldfile,'field')
filename =strmid(fieldfile,0,ipos)
imagefile0=strmid(filename,0,ipos-2)+'A_image.sav'	;171 image
imagefile1=strmid(filename,0,ipos-2)+'B_image.sav'	;195 image
imagefile2=strmid(filename,0,ipos-2)+'C_image.sav'	;284 image
i1      =fov(0)
j1      =fov(1)
i2      =fov(2)
j2      =fov(3)
nx      =(i2-i1+1)
ny      =(j2-j1+1)
r_n     =1.0
im      =(i1+i2)/2.
jm      =(j1+j2)/2.
rm      =sqrt((im-crpix1)^2+(jm-crpix2)^2)
im_b    =im
jm_b    =jm
if (rm lt rad_pix) then euvi_stereo_rot,im,jm,para,r_n,im_b,jm_b,zm_b
i1_b    =im_b-(im-i1)
i2_b    =i1_b+(i2-i1)
;...............................................................................
image0	=fltarr(nx,ny,6)				 ;Input images
image1  =fltarr(nx,ny,6)				 ;Model images
restore,imagefile0
image0(*,*,0)=reform(image_pair(1,i1:i2,j1:j2)) 	;171/SC A
image0(*,*,3)=reform(image_pair(0,i1_b:i2_b,j1:j2))     ;171/SC B
;...............................................................................
restore,imagefile1
image0(*,*,1)=reform(image_pair(1,i1:i2,j1:j2)) 	;195/SC A
image0(*,*,4)=reform(image_pair(0,i1_b:i2_b,j1:j2))     ;195/SC B
;...............................................................................
restore,imagefile2
image0(*,*,2)=reform(image_pair(1,i1:i2,j1:j2)) 	;284/SC A
image0(*,*,5)=reform(image_pair(0,i1_b:i2_b,j1:j2))     ;284/SC B
;...............................................................................

;_______________________READ LOOP FIELD LINES___________________________
idat	=strpos(fieldfile,'.dat')
isav	=strpos(fieldfile,'.sav')
if (idat gt 0) then readcol,fieldfile,iloop_,ix_,iy_,iz_,ih_,iz_raw,dh_,nsm_,wave_,chi2_,q1_,q2_
if (isav gt 0) then restore,fieldfile ;-->iloop_,ix_,iy_,iz_
nloop  =long(max(iloop_)+1)
print,'Number of previously traced loops    NLOOP   =',nloop

;_______________________READ EUVI RESPONSE FUNCTION____________________	
wave_	=['171','195','284']
nte_	=100
te_     =3.*findgen(nte_)/float(nte_)    ;MK loop temperature
resp_	=fltarr(nte_,6)
sc_	=['a','b']
for iwave=0,2 do begin
for isc=0,1 do begin
 wave	=wave_(iwave)
;euvi_temp_resp,wave,sc_(isc),te_,resp_iwave,0  ;analytical
 euvi_temp_resp,wave,sc_(isc),te_,resp_iwave,1  ;EUVI response function 
 resp_(*,iwave+3*isc)=resp_iwave
endfor
endfor

;_______________________DEFINES POINT SPREAD FUNCTION____________________
w	=(w_km/725.)/cdelt1			;width in pixels
w_res	=2.15					;resolution in pixels 
w_obs	=sqrt(w^2+w_res^2)			;observed width
sig_gauss=w_obs/(2.*sqrt(2*alog(2.)))		;gaussian half width in pixels
sig_gauss2=sig_gauss*4.
nk	=long(3.*sig_gauss2)
nk	=(2*nk)/2-1
kernel	=fltarr(nk,nk)
for i=0,nk-1 do begin
 for j=0,nk-1 do begin
 kernel(i,j)=0.9*exp(-((i-nk/2)^2+(j-nk/2)^2)/(2.*sig_gauss^2))$
            +0.1*exp(-((i-nk/2)^2+(j-nk/2)^2)/(2.*sig_gauss2^2))
 endfor
endfor
kernel	=kernel/total(kernel)			;normalization

;_______________________EUVI POINT SPREAD FUNCTION_______________________
;psf	=euvi_psf(nk,index_a)	
;kernel  =psf/total(psf)				;normalization

;_______________________BACKGROUND_____________________________________
backgr	=fltarr(3,2)
peak	=fltarr(3,2)
for iwave=0,2 do begin
for isc=0,1 do begin
 iw	=iwave+3*isc
 peak(iwave,isc)  =max(image0(*,*,iw))
 backgr(iwave,isc)=min(image0(*,*,iw))
;backgr(iwave,isc)=median([reform(image0(0,*,iw)),reform(image0(nx-1,*,iw)),$
;                reform(image0(*,0,iw)),reform(image0(*,ny-1,iw))])
 image1(*,*,iw)=image1(*,*,iw)+backgr(iwave,isc)
endfor
endfor

;_______________________SPIKE FILTER___________________________________
if (spike gt 0.) then begin
 for iw=0,5 do image0(*,*,iw)=trace_unspike(image0(*,*,iw),sens=spike)
endif

;_______________________SORTING LOOPS BY BRIGHTNESS______________________
length=fltarr(nloop)
qflux  =fltarr(nloop)
im1=image0(*,*,0) &tot1=total(im1)
im2=image0(*,*,1) &tot2=total(im2)
im3=image0(*,*,2) &tot3=total(im3)
for iloop=0,nloop-1 do begin
 ind    =where(iloop_ eq iloop,np)
 if (np ge 3) then begin
  xp     =ix_(ind)
  yp     =iy_(ind)
; len	 =total(sqrt((xp(1:np-1)-xp(0:np-2))^2+(yp(1:np-1)-yp(0:np-2))^2))
; length(iloop)=len
  ix	=long(xp-i1+0.5)
  iy	=long(yp-j1+0.5)
  ind_ij=where((ix ge 0) and (ix le nx-1) and (iy ge 0) and (iy le ny-1))
  ij 	=ix(ind_ij)+iy(ind_ij)*nx
  q1    =total(im1(ij))/tot1
  q2    =total(im2(ij))/tot2
  q3    =total(im3(ij))/tot3
  qflux(iloop)=(q1+q2+q3)
 endif
endfor
;prio	=max(length)-length
 prio	=max(qflux)-qflux
isort   =sort(prio) 
;print,'Range of loop lengths = ',minmax(prio)
 print,'Loop flux fractions    = ',minmax(prio)

;_______________________ITERATIVE HYDRO LOOP FITTING___________________
ds	=1.0
para_model=fltarr(nloopmax,11)
temin   =0.01 ;MK
temax   =10.0 ;MK
inum	=0

for iorder=0,nloop-1 do begin
 iloop  =isort(iorder)

 ind    =where(iloop_ eq iloop,np)
 if (np lt 3) then goto,endloop
 xp     =ix_(ind)
 yp     =iy_(ind)
 zp     =iz_(ind)
 euvi_spline,para,xp,yp,zp,x,y,z,d,r,r_n
 n      =n_elements(x)
 euvi_stereo_rot,x,y,para,r_n,x_b,y_b,z_b

 s	=fltarr(n)
 for is=1,n-1 do s(is)=s(is-1)+$
        sqrt((x(is)-x(is-1))^2+(y(is)-y(is-1))^2+(z(is)-z(is-1))^2)
 ns	=long(max(s)/ds)
 s_	=ds*findgen(ns)
 if (ns le 3) then goto,endloop
 xs	=interpol(x,s,s_)
 ys	=interpol(y,s,s_)
 rs	=interpol(r,s,s_)
 hs	=rs-rad_pix			;height in pixels
 ix	=long(xs-i1+0.5)
 iy	=long(ys-j1+0.5)
 slen	=max(s_)-min(s_)
 hmax	=max(hs)
 ind_ij	=where((ix ge 0) and (ix le nx-1) and (iy ge 0) and (iy le ny-1))
 ija	=ix(ind_ij)+iy(ind_ij)*nx

 xs_b	=interpol(x_b,s,s_)
 ys_b	=interpol(y_b,s,s_)
 ix_b	=long(xs_b-i1_b+0.5)
 iy_b	=long(ys_b-j1  +0.5)
 ijb	=ix_b(ind_ij)+iy_b(ind_ij)*nx
 w_sm	=7<(ns-1)

 for imulti=0,nmulti-1 do begin

  flux_obs=fltarr(ns,6)
  for isc=0,1 do begin
   if (isc eq 0) then ij=ija
   if (isc eq 1) then ij=ijb
   for iwave=0,2 do begin 
    iw	=iwave+isc*3
    newimage=image0(*,*,iw)-image1(*,*,iw)
    flux_ob1 =fltarr(ns)
    flux_ob1(ind_ij)=newimage(ij)
    flux_obs(*,iw)=smooth(flux_ob1,w_sm)>0.
    if (max(flux_obs(*,iw)) le 0) then goto,endloop
   endfor
  endfor
  tem	=(2.0-0.5*float(imulti))>1.0  
  lam	=0.5     ;[x100 Mm] pressure scale height
  hchr	=2.0*(cdelt1*0.725)
  h0	=hchr 	 ;first footpoint height 
  h1	=hchr    ;second footpoint height
  coeff	=[tem,h0,lam,h1,lam]
  vector=[0.2,5.0,0.2,5.0,0.2]
  ncoeff=n_elements(coeff)
  vv    =fltarr(ncoeff,ncoeff)
  for i=0,ncoeff-1 do vv(i,i)=vector(i)
  ftol   =1.0e-4 
  powell,coeff,vv,ftol,dev,'euvi_hydro_fit_powell'
  te      =par_model(0)          
  h0      =par_model(1)          
  lam     =par_model(2)  
  h1      =par_model(3)          
  lam1    =par_model(4)  
  em44    =par_model(5) 
  if (em44 lt em44_min) or (em44 gt em44_max) or $
     (te lt temin) or (te gt temax) then goto,endpara
  para_model(inum,*)=[par_model(0:7),w,ds,iloop]

  for iwave=0,2 do begin
   for isc=0,1 do begin
    iw	=iwave+3*isc
    image_slice=fltarr(nx,ny)
    image_convol=fltarr(nx,ny)
    image_slice(ij)=image_slice(ij)+reform(flux_model(ind_ij,iwave))
    zmax=total(image_slice) 
    if (zmax gt 0) then begin
     image_convol=convol(image_slice,kernel)
     image1(*,*,iw)=image1(*,*,iw)+image_convol
    endif
   endfor
  endfor

;__________________________Test display________________________________
if (test ge 1) then begin
 window,1,xsize=1024,ysize=768
 clearplot
 loadct,5
 for iw=0,5 do begin
  isc  =iw/3
  iwave=iw mod 3
  x1_=0.05+0.16*iw
  x2_=x1_+0.14
  y1_=0.76
  y2_=0.96
  ima	=image0(*,*,iw)
  io	=0
  if (io eq 0) then begin
   nxw   =long(!d.x_vsize*(x2_-x1_)+0.5)
   nyw   =long(!d.y_vsize*(y2_-y1_)+0.5)
   z     =congrid(ima,nxw,nyw)
  endif
  if (io ne 0) then z=ima
  !p.position=[x1_,y1_,x2_,y2_]
  !x.range=[0,nx]
  !y.range=[0,ny]
  !x.style=1
  !y.style=1
  !x.title=''
  !y.title=''
  plot,[0,0],[0,0]
  z1=backgr(iwave,isc)
  z2=peak(iwave,isc)
  tv,bytscl(z,min=z1,max=z2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal
  if (isc eq 0) then oplot,ix,iy
  if (isc eq 1) then oplot,ix_b,iy_b
  !noeras=1

  y1_=0.55
  y2_=0.75
  ima	=image1(*,*,iw)
  io	=0
  if (io eq 0) then begin
   nxw   =long(!d.x_vsize*(x2_-x1_)+0.5)
   nyw   =long(!d.y_vsize*(y2_-y1_)+0.5)
   z     =congrid(ima,nxw,nyw)
  endif
  if (io ne 0) then z=ima
  !p.position=[x1_,y1_,x2_,y2_]
  plot,[0,0],[0,0]
  tv,bytscl(z,min=z1,max=z2),x1_,y1_,xsize=(x2_-x1_),ysize=(y2_-y1_),/normal
  if (isc eq 0) then oplot,ix,iy
  if (isc eq 1) then oplot,ix_b,iy_b

  y1_=0.30
  y2_=0.50
  col=50*(1+iwave) 
  !p.position=[x1_,y1_,x2_,y2_]
  !x.range=minmax(s)
  !y.range=[0,max(flux_obs(*,iw))*1.1]
  plot,s_,flux_obs(*,iw),psym=10
  oplot,s_,flux_model(*,iwave),color=col

  y1_=0.05
  y2_=0.25
  !p.position=[x1_,y1_,x2_,y2_]
  !y.title='Temperature, Flux'
  !y.range=[0,temax]
  plot,s_,temp_model,thick=3
  oplot,!x.range,[1,1]*(1.0+0.5*iwave),color=col
 endfor
 print,'#,T0,h0,lam,h1,lam1,EM = ',iloop,par_model(0:5),format='(a,i6,6f8.2)'
 if (test eq 1) then begin
  read,'continue? [0=stop ,1=bad fit, 2=good fit]',yes
  if (yes eq 0) then stop
 endif
endif
;__________________________________________________________________________
if ((iorder mod 10) eq 0) and (imulti eq 0)  then begin
  qf1	=total(image1(*,*,0))/total(image0(*,*,0))
  qf2	=total(image1(*,*,1))/total(image0(*,*,1))
  qf3	=total(image1(*,*,2))/total(image0(*,*,2))
  qf	=(qf1+qf2+qf3)/3.
  print,'N_field, N_model, L, f_mod/f_obs = ',iorder,inum,length(iloop),qf,format='(a,2i8,f9.1,f9.3)'
endif

  inum=inum+1
  if (inum ge nloopmax) then goto,complete
  endpara:
 endfor  ;for iloop=0,nloop-1 do begin 	
 endloop: 
endfor	;for imulti=0,nmulti-1 do begin
complete:
para_model=para_model(0:inum-1,*)
nloopmodel=inum

;________________________SAVE_______________________________________________
save,filename=hydrofile,image0,image1,para_model,fov,para,nmulti,nloopmodel,nk,em44_min,w_km,temin,temax
print,'Hydro parameters written in file =',hydrofile
end