pro euvi_mdi_field,fileset,fov,ds,hmin,hmax,model,thresh
;+
; 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
;
; 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_____________________________'
nfmax   =long(1./ds)
print,'MAXIMUM NUMBER OF POINTS PER LOOP = ',nfmax
a0_mis	=0.
hcurr	=-ds		;minimum height for current loops

 print,'_____________POTENTIAL FIELD___________________________________'
 nx	=50
 ny	=50
 ncmax	=nx*ny
 x1	=fov(0)
 x2	=fov(2)
 x	=x1+(x2-x1)*findgen(nx)/float(nx-1)
 y1	=fov(1)
 y2	=fov(3)
 y	=y1+(y2-y1)*findgen(ny)/float(ny-1)
 field_p =fltarr(nfmax,ncmax,4)
 ic	=0
 for i=0,nx-1 do begin
 for j=0,ny-1 do begin
  xm	=x(i)
  ym	=y(j)
  rm	=1.
  zm	=sqrt(rm^2-xm^2-ym^2)
  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
  bgauss=max(field(*,3))
  if (bgauss ge thresh) then begin
   for k=0,3 do field_p(0:nf-1,ic,k)=field(*,k)
   print,'fieldline #i,j,ic,b=',i,j,ic,bgauss
   ic	=ic+1
  endif
  if (ic eq ncmax) then goto,end_field 
 endfor
 endfor
 end_field:
 nc	=ic
 field_	=fltarr(nfmax,nc,4)
 for ic=0,nc-1 do for k=0,3 do field_(*,ic,k)=field_p(*,ic,k)

print,'_____________SAVE FIELDLINES___________________________________'
code	='h'	;helical flux tube fitting
fieldfile=fileset+'_loop_'+model+code+'.sav'
save,filename=fieldfile,field_
print,'Field lines saved in file : ',fieldfile
end