pro aia_mag_loop,fileset,fov,ds,hmin,hmax,model,nseg,qs,nrad,test
;+
; 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,ds,hmax,model,nseg,test
;
; 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
;		model	 : 'p'=potential, 'n'=nonpotential
;               nseg     = number of segments (misalignment angle)
;		test	 = display initial-guess values 
;
; Outputs:      *_field.sav = savefile containing fieldlines
;
; History     : 12-May-2011, Version 1 written by Markus J. Aschwanden
;                6-Sep-2011, modifications
;
; Contact     : aschwanden@lmsal.com
;-

print,'_____________READ MAGNETIC FIELD MODEL_________________________'
coeff_file=file_search(fileset+'_coeff.dat')
readcol,coeff_file(0),fmag,xmag,ymag,zmag 
coeff	=[[fmag],[xmag],[ymag],[zmag]]
ncomp	=n_elements(fmag)
print,'NUMBER OF MAGNETIC COMPONENTS = ',ncomp

print,'_____________READ LOOP COORDINATES_____________________________'
loopfile  =fileset+'*_loop.dat'	
readcol,loopfile,iloop,x,y,z,s
nloop   =long(max(iloop)+1)
print,'NUMBER OF LOOPS = ',nloop
nfmax   =long(1./ds)
print,'MAXIMUM NUMBER OF POINTS PER LOOP = ',nfmax
dev_loop=fltarr(nloop)
alpha0	=0.
hcurr	=-ds		;minimum height for current loops
t1	=systime(0,/seconds)

if (model eq 'p') then begin
 print,'_____________POTENTIAL FIELD___________________________________'
 field_p =fltarr(nfmax,nloop,4)
 for il=0,nloop-1 do begin
  ind   =where(il eq iloop,ns)
  xl	=x[ind]		
  yl	=y[ind]	
  zl	=z[ind]		
  imid	=long(qs*(ns-1))
  xm	=xl(imid)
  ym	=yl(imid)
  zm	=zl(imid)
  aia_mag_fieldline, xm,ym,zm,ds,hcurr,hmax,coeff,curr0,alpha0,curr,nf,nm ;>curr
  aia_mag_fieldline, xm,ym,zm,ds,hmin,hmax,coeff,curr ,alpha0,field,nf,nm
  aia_mag_fieldangle,xl,yl,zl,ds,hmin,hmax,coeff,curr ,alpha0,nseg,dev,dev_
  for j=0,3 do field_p(0:nf-1,il,j)=field(*,j)
  dev_loop[il]=dev
 endfor
endif

if (model eq 'n') then begin
 print,'_____________NON-POTENTIAL FIELD_______________________________'
 alpha_=fltarr(nloop)
 xf1_  =fltarr(nloop)
 yf1_  =fltarr(nloop)
 curr_ =fltarr(nfmax,nloop,4)			  ;current field lines
 field_=fltarr(nfmax,nloop,4)			  ;nonpot field lines
 for il=0,nloop-1 do begin
  ind   =where(il eq iloop,ns)
  xl	=x[ind]					  ;loop coordinates
  yl	=y[ind]	
  zl	=z[ind]		
  imid	=long(qs*(ns-1))			  ;midpoint of loop
  xm	=xl(imid)
  ym	=yl(imid)
  zm	=zl(imid)
  aia_mag_fieldline,xm,ym,zm,ds,hmin,hmax,coeff,curr0,alpha0,field_p,nf,nm  
  im	=(imid > 1) < (ns-2)
  nm	=(nm   > 1) < (nf-2)
  xp	=field_p(*,0)
  yp	=field_p(*,1)
  zp	=field_p(*,2)
  devbest=90.
  if (nf le 2) or (nm le 2) then goto,end_loop
  v_loop=[xl(im+1)-xl(im-1),yl(im+1)-yl(im-1),zl(im+1)-zl(im-1)]  
  v_pot =[xp(nm+1)-xp(nm-1),yp(nm+1)-yp(nm-1),zp(nm+1)-zp(nm-1)]  
  vector_product,v_pot,v_loop,v_rad,v_rad_norm,a_mis 
  aia_mag_len,xl,yl,zl,ns,slen,s
  a_mis_deg=a_mis*(180./!pi)
  ntwistmax=2.

  rad_	=ds*(findgen(nrad)-nrad/2)
  xr_	=fltarr(nrad,3)
  dev_r =fltarr(nrad)
  field_test=fltarr(nfmax,nrad,4)		  ;nonpot field lines
  curr_test =fltarr(nfmax,nrad,4)		  ;current field lines
 for ir=0,nrad-1 do begin
   rad	=rad_(ir)
   tan_mu=tan(a_mis) 
   rmin	=abs(slen*tan_mu/(2.*!pi*ntwistmax))
   if (ir eq 0) then rmin=0.			  ;at least one case 
   if (rad ne 0.) then alpha1=(2./rad)*tan_mu/(1.+tan_mu^2)
   if (rad eq 0.) then alpha1=0.
   alpha_(ir)=alpha1
   xr	=xm+v_rad_norm(0)*rad		;v_rad points towards twist axis
   yr	=ym+v_rad_norm(1)*rad  
   zr	=zm+v_rad_norm(2)*rad	
   hr	=sqrt(xr^2+yr^2+zr^2)-1.
   if (hr ge hmin) and (hr le hmax) and (abs(rad) gt rmin) then begin
    xr_[ir,0]=xr
    xr_[ir,1]=yr
    xr_[ir,2]=zr
    aia_mag_fieldline, xr,yr,zr,ds,hcurr,hmax,coeff,curr0,alpha0,curr,nc,nm 
    aia_mag_fieldline, xm,ym,zm,ds,hmin,hmax,coeff,curr,alpha1,field,nf,nm  
    aia_mag_fieldangle,xl,yl,zl,ds,hmin,hmax,coeff,curr,alpha1,nseg,dev,dev_
    dev_r[ir]=dev
    for j=0,3 do begin
     field_test(*,ir,j)=0.
     curr_test(*,ir,j)=0.
     field_test(0:nf-1,ir,j)=field(*,j)
     curr_test(0:nc-1,ir,j)=curr(*,j)
    endfor
    if (abs(dev) lt devbest) then begin
     devbest=abs(dev)
     alpha=alpha1
     radbest=rad
     xf1=curr(0,0)
     yf1=curr(0,1) 
     for j=0,3 do begin
      field_(*,il,j)=0.			;delete previous residuals
      curr_(*,il,j)=0.
      field_(0:nf-1,il,j)=field(*,j)
      curr_(0:nc-1,il,j)=curr(*,j)
     endfor
    endif
    if (test ge 1) then print,'rad,alpha,a_mis,dev=',rad,alpha1,a_mis_deg,dev 
   endif
  endfor	;for ir=0,nrad-1 
  dev_loop[il]=devbest
  print,'LOOP,RAD,ALPHA,A_MIS,DEV=',il,radbest,alpha1,a_mis_deg,devbest,$
   format='(a,i4,4f10.4)'
  if (test ge 1) then aia_mag_loop_test,il,xl,yl,zl,nseg,hmax,imid_,xf1,yf1,$
   field_test,curr_test,xr_,rad_,dev_r,field_,curr_,test
  END_LOOP:
 endfor		;for il=0,nloop-1 
endif		;if model eq 'n'
t2      =systime(0,/seconds)
cpu     =t2-t1
print,'Computation time = ',cpu,' s'

print,'_____________MEDIAN MISALIGNMENT ANGLE_________________________'
dev     =median(abs(dev_loop))
print,'Misalignment median = ',dev,' deg',format='(a,f6.2,a)'

print,'_____________SAVE FIELDLINES___________________________________'
code	='h'	;helical flux tube fitting
fieldfile=fileset+'_field_'+model+code+'.sav'
if (model eq 'p') then save,filename=fieldfile,field_p,dev_loop,cpu
if (model eq 'n') then save,filename=fieldfile,field_ ,dev_loop,cpu
print,'Field lines saved in file : ',fieldfile
end