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_____________________________'
fileset0=strmid(fileset,0,17)			;CHANGE OCT 15, 2011
loopfile  =fileset0+'_a_coord.dat'		;Earth/MDI coordinates
print,'read : ',loopfile
readcol,loopfile,iloop,x,y,z 			;CHANGE OCT 15, 2011
nloop   =long(max(iloop)+1)
print,'NUMBER OF LOOPS = ',nloop
nfmax   =long(1./ds)
print,'MAXIMUM NUMBER OF POINTS PER LOOP = ',nfmax
a0_mis	=0.
hcurr	=-ds		;minimum height for current loops
dev_loop=fltarr(nloop)
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,a0_mis,curr,nf,nm ;>curr
  aia_mag_fieldline, xm,ym,zm,ds,hmin,hmax,coeff,curr ,a0_mis,field,nf,nm
  aia_mag_fieldangle,xl,yl,zl,ds,hmin,hmax,coeff,curr ,a0_mis,nseg,dev,dev_
  for j=0,3 do field_p(0:nf-1,il,j)=field(*,j)
  dev_loop[il]=dev
  print,'fieldline #',il,' ns=',ns,' nf=',nf,' dev=',dev
 endfor
endif

if (model eq 'n') then begin
 print,'_____________NON-POTENTIAL FIELD_______________________________'
 alpha_loop=fltarr(nloop)
 rad_loop=fltarr(nloop)
 xf1_loop=fltarr(nloop)
 yf1_loop=fltarr(nloop)
 len_loop=fltarr(nloop)
 amis_loop=fltarr(nloop)
 ntwist_loop=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,a0_mis,field_p,nf,nm  
  im	=(imid > 0) < (ns-2)
  nm	=(nm   > 0) < (nf-2)
  xp	=field_p(*,0)
  yp	=field_p(*,1)
  zp	=field_p(*,2)
  dev_best=90.
  if (nf le 2) or (nm le 2) then goto,end_loop
  v_loop=[xl(im+1)-xl(im),yl(im+1)-yl(im),zl(im+1)-zl(im)]  
  v_pot =[xp(nm+1)-xp(nm),yp(nm+1)-yp(nm),zp(nm+1)-zp(nm)]  
  vector_product,v_loop,v_pot,v_rad,v_rad_norm,aa_mis 
  aia_mag_len,xl,yl,zl,ns,slen,s
  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)
   if (rad ge 0) then a_mis=+aa_mis		  ;positive direction r>0
   if (rad lt 0) then a_mis=-aa_mis		  ;negative direction r<0
   a_mis_deg=a_mis*(180./!pi)
   tan_mu=tan(a_mis) 
   rmin	=(slen*abs(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.
   xr	=xm+v_rad_norm(0)*rad	
   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,a0_mis,curr,nc,nm 
    aia_mag_fieldline, xm,ym,zm,ds,hmin,hmax,coeff,curr,a_mis,field,nf,nm  
    aia_mag_fieldangle,xl,yl,zl,ds,hmin,hmax,coeff,curr,a_mis,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 dev_best) and (nc ge 3) and (nf ge 3) then begin
     dev_best=abs(dev)
     alpha_best=alpha1
     rad_best=rad
     xf1_best=curr(0,0)
     yf1_best=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]=dev_best
  rad_loop[il]=rad_best
  alpha_loop[il]=alpha_best
  xf1_loop[il]=xf1_best
  yf1_loop[il]=yf1_best
  len_loop[il]=slen
  amis_loop[il]=aa_mis
  ntwist_loop[il]=(slen*tan(aa_mis))/(2.*!pi*rad_best)
  print,'LOOP,RAD,ALPHA,A_MIS,DEV=',il,rad_best,alpha_best,a_mis_deg,dev_best,$ 
    format='(a,i4,4f10.4)'
  if (test ge 1) then begin
   a_mis_s=fltarr(ns)
   for is=0,ns-2 do begin 
    x1	=xl(is)
    y1	=yl(is)
    z1	=zl(is)
    c0	=0.	;curr
    a0	=0.	;a0_mis
    polar=1
    aia_mag_vector,x1,y1,z1,ds,hmin,hmax,coeff,c0,a0,polar,x2,y2,z2,b_tot,endv
    v_loop=[xl(is+1)-xl(is),yl(is+1)-yl(is),zl(is+1)-zl(is)]  
    v_pot =[x2-x1,y2-y1,z2-z1]
    vector_product,v_loop,v_pot,v_rad,v_rad_norm,aa_mis 
    a_mis_s(is)=aa_mis
   endfor
   aia_mag_loop_test,il,xl,yl,zl,nseg,hmax,imid_,xf1_best,yf1_best,$
     field_test,curr_test,xr_,rad_,dev_r,field_,curr_,fov,dev_best,$
     a_mis_s,test
  endif
  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+'_loop_'+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,rad_loop,$
 xf1_loop,yf1_loop,len_loop,amis_loop,ntwist_loop,alpha_loop,cpu
print,'Field lines saved in file : ',fieldfile
end