pro magn_fit,modelfile,runfile,para,meth
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO 
;
; Name        : MAGN_FIT
;
; Category    : Magnetic data modeling 
;
; Explanation : Fitting of twisted field lines (nonlinear force-free
;		field) NLFFF-model to observed loop projections.
;		The fitting constraints can be stereoscopic
;               loop coordinates (qv=0.0), photospheric 3D magnetic
;		field vector data (qv=1.0) or a combination of both.
;
; Syntax      : IDL>magn_fit,modelfile,runfile,para,meth
;
; Inputs      : modelfile = input filename *_cube.sav, *_loops.sav
;		runfile   = input/output files *_coeff.dat, *_fit.sav
;               para      = input parameters
;                          [ds,ngrid,nmag,nsm,qloop,nseg,acc,amax,hmax,dpix,
;                           niter,qzone,eps,qv,thresh,q_mag,halt,iopt]
;		meth	 = 'A' optimization method
;
; Outputs     ; runfile   = output files *_coeff.dat, *_fit.sav
;
; History     : 16-May-2011, Version 1 written by Markus J. Aschwanden
;               30-Jan-2012, Version 2 released to SSW
;		20-Feb-2012, Adjust da0=amax_lin (if iter ge 1)
;
; Contact     : aschwanden@lmsal.com
;-

print,'___________________MAGN_FIT_______________________________'
t1     =systime(0,/seconds)
ds      =para[0]
ngrid   =para[1]
nmag    =para[2]
nsm     =para[3]
qloop   =para[4]
nseg    =para[5]
acc     =para[6]
amax	=para[7]
hmax	=para[8]
dpix	=para[9]
niter   =para[10]
qzone	=para[11]
eps	=para[12] 
qv	=para[13]
thresh	=para[14]
iter	=0

;___________________input of 3D vector magnetograph data______________
cubefile=modelfile+'_cube.sav'
restore,cubefile		;-->x,y,z,bx,by,bz,b,a
if (qv eq 0.) then field_vec=fltarr(2,3,1)
if (qv gt 0.) then begin
 dim	=size(b)
 nx	=dim[1]
 ny	=dim[2]
 nz	=dim[3]
 ind	=where(b(*,*,0) ge thresh,nv)	&print,'Number of field vectors = ',nv
 field_vec=fltarr(2,3,nv)
 k	=0
 m	=0
 for i=0,nx-1 do begin
  for j=0,ny-1 do begin
   if (b(i,j,k) ge thresh) then begin
    cos_x=bx(i,j,k)/b(i,j,k)
    cos_y=by(i,j,k)/b(i,j,k)
    cos_z=bz(i,j,k)/b(i,j,k)
    field_vec(*,0,m)=x[i]+[0,ds*cos_x]
    field_vec(*,1,m)=y[j]+[0,ds*cos_y]
    field_vec(*,2,m)=z[0]+[0,ds*cos_z]
    m	=m+1
   endif
  endfor
 endfor
endif

;______________________READ DATA_____________________________________
savefile  =modelfile+'_loops.sav'
restore,savefile      ;-->bzmap,x,y,field_lines,q_scale
dim	=size(field_lines)
if (dim(0) le 2) then nf=1
if (dim(0) eq 3) then nf=dim(3)
nx	=n_elements(x)
ny	=n_elements(y)
polar	=1.
xrange	=x[nx-1]-x[0]
dpix	=x[1]-x[0]
qstep	=0.5

coeff_file=runfile+'_coeff.dat'
readcol,coeff_file,m1,x1,y1,z1,a1
nmag	=n_elements(m1)
a1	=a1*0.				;reset twist 
coeff	=[[m1],[x1],[y1],[z1],[a1]]

;_____________________START OF FITTING ITERATION______________________'
magn_misalign,field_lines,nf,nseg,ds,coeff,polar,dev_min,dev_loop,field_vec,qv
angle	=90.
if (niter eq 0) then goto,end_iter
for iter=0,niter-1 do begin

;______________HIERARCHICAL MAGNETIC ZONES (METHOD A)______________'
 n_zone =nmag
 dist0	=max(sqrt((x1-x1(0))^2+(y1-y1(0))^2+(z1-z1(0))^2))>(dpix*3.)
 rad	=(dist0*qzone^iter) > (dpix/2.)  
 i_zone  =lonarr(nmag)
 for iz=1,nmag do begin
  ind	=where(i_zone eq 0,nind)
  if (nind ge 1) then begin
   i0	=ind(0)
   x0	=x1[i0]
   y0	=y1[i0]
   z0	=z1[i0]
   dist=sqrt((x1-x0)^2+(y1-y0)^2+(z1-z0)^2)
   irad=where(dist le rad) 	
   i_zone(irad)=iz
  endif
 endfor
 n_zone=max(i_zone)

;______________________DISPLAY MAGNETIC ZONES______________________
 window,1,xsize=512,ysize=512
 clearplot
 if (iter eq 0) then loadct,5
 sym	=[1,2,4,5,6,7]
 plot,[0,0],[0,0],xrange=minmax(x),yrange=minmax(y)
 phi	=2*!pi*findgen(101)/100.
 for iz=1,n_zone do begin
  ind	=where(i_zone eq iz)
  i0	=ind(0)
  color	=50+205*float(iz)/float(n_zone)
  oplot,x1(ind),y1(ind),psym=sym(iz mod 5),color=color
  oplot,x1(i0)+rad*cos(phi),y1(i0)+rad*sin(phi),color=color
 endfor

 ;__________________________OPTIMIZATION_________________________
 a_min	=0.
 b_min	=0.
 for iz=1,n_zone do begin
  ind_zone=where(i_zone eq iz,nind)
  a0	=avg(coeff[ind_zone,4])	
  da0	=amax
  if (iter ge 1) then da0=amax_lin
  a_	=a0+da0*[-1.,0,1.]
  nitmax=long(alog10(acc)/alog10(qstep)) > 1
  for it=0,nitmax-1 do begin
   nv	=n_elements(a_)
   b_	=fltarr(nv)
   for iv=0,nv-1 do begin 		;vary parameter
    coeff2=coeff
    coeff2(ind_zone,4)=a_(iv) 
    magn_misalign,field_lines,nf,nseg,ds,coeff2,polar,dev,dev_loop,field_vec,qv
    b_(iv)=dev
   endfor            
   if (it eq 0) then begin &a_all=a_         &b_all=b_         &endif
   if (it ge 1) then begin &a_all=[a_all,a_] &b_all=[b_all,b_] &endif
   b_min =min(b_all,imin)
   a_min =a_all(imin)
   minimum_iterate,a_all,b_all,a1_,a2_,it,qstep
   a_	=[a1_,a2_]			;new values 
  endfor	;for it=0,nitmax-1
  coeff(ind_zone,4)=a_min
  NEXT_ZONE:
  print,'ITER=',iter+1,'/',niter,'  IZ=',iz,'/',n_zone,a_min,b_min,$
    format='(4(a,i3),2f10.3)'
  if (iter eq 0) and (iz eq 1) then amax_lin=abs(a_min)	;linear force-free alpha
 endfor  ;for iz=0,n_zone-1 
 magn_misalign,field_lines,nf,nseg,ds,coeff,polar,dev_min,dev_loop,field_vec,qv
 print,'MISALIGNMENT ANGLE ------------------> ',dev_min,format='(a,f12.3)'
 if (dev_min ge (angle-eps)) and (rad le dpix) then goto,end_iter
 angle=dev_min
endfor	;for iter=0,niter-1 
iter	=niter-1
END_ITER:
angle	=angle < dev_min
t2     =systime(0,/seconds)
cpu    =t2-t1
print,'FORWARD FITTING : CPU time = ',t2-t1

;_______________________SAVE COEFFICIENTS_____________________
openw,2,coeff_file
for is=0,nmag-1 do begin
 m1     =coeff(is,0)
 x1     =coeff(is,1)
 y1     =coeff(is,2)
 z1     =coeff(is,3)
 a1     =coeff(is,4)
 printf,2,m1,x1,y1,z1,a1
endfor
close,2

;_______________________SAVE PARAMETERS________________________
SAVE_DATA:
savefile2=runfile+'_fit.sav'
save,filename=savefile2,bzmap,x,y,field_lines,alpha,q_scale,angle,cpu,iter,para,dev_loop
help,dev_loop
statistic,dev_loop
print,'parameters saved in file = ',savefile2
end