pro magn_field,runfile,para
;+
; Project     : SOHO/MDI, SDO/HMI, STEREO
;
; Name        : MAGN_FIELD
;
; Category    : Magnetic data modeling
;
; Explanation : Calculating NF magnetic field lines FIELD_LINES[NS,5,NF] 
;		magnetic field coefficients [B1,X1,Y1,Z1,A1] 
;
; Syntax      : IDL>magn_field,runfile,para
;
; Inputs      : runfile   = input file *_coeff.dat
;               para      = input parameters
;                          [ds,ngrid,nmag,nsm,qloop,nseg,acc,amax,hmax,dpix,
;                           niter,qzone,eps,qv,thresh,q_mag,halt,iopt]
;
; Outputs     ; runfile   = output files *_fit.sav
;
; History     : 16-May-2011, Version 1 written by Markus J. Aschwanden
;               30-Jan-2012, Version 2 released to SSW
;
; Contact     : aschwanden@lmsal.com
;-


print,'________________MAGN_FIELD____________________________________'
ds      =para[0]
qloop   =para[4]
dpix	=para[9]
halt	=para[16]

;______________________READ DATA_____________________________________
savefile=runfile+'_fit.sav'
restore,savefile      ;-->bzmap,x,y,field_lines,q_scale,angle,cpu,iter
dpix	=x[1]-x[0]
;zmax	=max(field_lines(*,2,*)) 
;halt	=(zmax-1.+dpix)*1.5	
nz	=long(halt/ds)
;ns	=long(halt*4*!pi/ds)
nx	=n_elements(x)
ny	=n_elements(y)
dim	=size(field_lines)
if (dim(0) le 2) then nf=1
if (dim(0) eq 3) then nf=dim(3)
field_loops=field_lines

coeff_file=runfile+'_coeff.dat'
readcol,coeff_file,m1,x1,y1,z1,a1
coeff	=[[m1],[x1],[y1],[z1],[a1]]

;'______________________FIELD LINES_______________________________________'
nsmax	=1000
field_lines=fltarr(nsmax,5,nf)
for k=0,nf-1 do begin
 ind	=where(field_loops(*,2,k) ne 0,np)
 is	=long(qloop*(np-1))
 xm	=field_loops(is,0,k)
 ym	=field_loops(is,1,k)
 zm	=field_loops(is,2,k)
 magn_fieldline,xm,ym,zm,ds,halt,coeff,field,ns_,nm
 ns_	=ns_ < nsmax
 field_lines(0:ns_-1,0,k)=field(0:ns_-1,0)
 field_lines(0:ns_-1,1,k)=field(0:ns_-1,1)
 field_lines(0:ns_-1,2,k)=field(0:ns_-1,2)
 field_lines(0:ns_-1,3,k)=field(0:ns_-1,3)	;b
 field_lines(0:ns_-1,4,k)=field(0:ns_-1,4)	;alpha
endfor

;_______________________SAVE PARAMETERS________________________
save,filename=savefile,bzmap,x,y,field_lines,alpha,q_scale,angle,cpu,iter,dev_loop
print,'parameters saved in file = ',savefile
end