pro aia_mag_vector,x1,y1,z1,ds,hmin,hmax,coeff,curr,alpha,polar,x2,y2,z2,b_tot,endv
;+
; Project     : AIA/SDO
;
; Name        : AIA_LOOP_FIELDLINE
;
; Category    : Data analysis   
;
; Explanation : Magnetic field line extrapolation
;               from midpoint position [xm,ym,zm] to both footpoints
;               [xf1,yf1,0] and [xf2,yf2,0]
;
; Syntax      : IDL>aia_loop_fieldline,xm,ym,zm,np,hmax,coeff,curr,field
;
; Inputs      : xm       = x-coordinate of loop midpoint (solar radii)
;               ym       = y-coordinate of loop midpoint 
;               zm       = z-coordinate of loop midpoint 
;               ds       = 3D loop increment 
;               hmax     = maximum height of (open) field line
;               coeff    = coefficients of magnetic potential field model
;
; Outputs     : field[ns,3] = coordinates of full loop field line
;
; History     : 1-Jun-2011, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

;_____________________POTENTIAL FIELD INPUT_________ __________________
dim     =size(coeff)
ncomp   =dim(1)
magn    =coeff(*,0)
xmag    =coeff(*,1)
ymag    =coeff(*,2)
zmag    =coeff(*,3)                     ;negative height
rmag    =sqrt(xmag^2+ymag^2+zmag^2)
hmag    =rmag-1.
pol     =magn/abs(magn)                 ;magnetic polarity

;....................potential field...................................
bx   	=0.
by   	=0.
bz   	=0.
for ic=0,ncomp-1 do begin
 d_mag=sqrt((x1-xmag(ic))^2+(y1-ymag(ic))^2+(z1-zmag(ic))^2)
 cos_x  =(x1-xmag(ic))/d_mag 
 cos_y  =(y1-ymag(ic))/d_mag
 cos_z  =(z1-zmag(ic))/d_mag
 b0_p   =magn(ic)*(hmag(ic)/d_mag)^2
 bx     =bx+b0_p*cos_x
 by     =by+b0_p*cos_y
 bz     =bz+b0_p*cos_z
endfor
b_pot	=[bx,by,bz]
b_p     =sqrt(total(b_pot^2))
b_tot   =b_p

;....................non-potential twisted field line..................
eps	=1.e-3
x2	=x1
y2	=y1
z2	=z1
endv	=0
if (alpha ne 0.) then begin
 dim	=size(curr)
 ncurr	=dim(1)
 if (ncurr le 3) then begin &endv=1 &goto,end_proc &endif
 x1_	=[x1,y1,z1]
 dist_perp_curve,x1_,curr(*,0),curr(*,1),curr(*,2),rc_,dr_curr,is_curr
 if (is_curr le 0.) or (is_curr ge (ncurr-2.)) then endv=1
 tan_mu	=0.
 if (alpha lt 0.) then dr_curr=-dr_curr
 c	=alpha*dr_curr				;alpha<0 for r<0 
 if (c ne 0.) then tan_mu =(1./c)+sqrt((1./c^2)-1.)  ;misalign angle mu(alpha)
 alpha2=(2./dr_curr)*tan_mu/(1.+tan_mu^2)
 if abs(alpha-alpha2) ge eps then print,alpha,alpha2
 b_para =b_p 	                      		;parallel B component 
 b_azi  =b_p*tan_mu                		;azimuthal B component 
 cos_para=b_pot/b_p	 			;cosine of potential field
 v_perp =reform(rc_-x1_)	    	    	;radial vector
 vector_product,v_perp,b_pot,v_azi,cos_azi,angle  
 b_np	=b_para*cos_para+b_azi*cos_azi		;non-potential vector
 if (dr_curr eq 0) then b_np=b_pot		;radial distance zero=potential
 bx     =b_np[0]
 by     =b_np[1]
 bz     =b_np[2]
 b_tot	=0.
 b_tot2	=total(b_np(0:2)^2)
 if (b_tot2 gt 0) then b_tot=sqrt(b_tot2)
endif
;......................incremental step along field line............
if (b_tot ne 0.) then begin
 x2	=x1 + ds*(bx/b_tot)*polar 
 y2	=y1 + ds*(by/b_tot)*polar 
 z2	=z1 + ds*(bz/b_tot)*polar 
 r2	=sqrt(x2^2+y2^2+z2^2)
 if (r2 le (1.+hmin)) then endv=1
endif
end_proc:
if (b_tot eq 0) then endv=1
end