pro aia_mag_fieldline,xm,ym,zm,ds,hmin,hmax,coeff,curr,a_mis,field,ns,nm
;+
; 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
;-

;____________________FIELD LINE HEIGHT LOOP_________________________
nfmax   =long(1./ds)
f   	=fltarr(nfmax,4)
f(0,0)	=xm
f(0,1)	=ym
f(0,2)	=zm
polar	=1.
for is=1,nfmax-1 do begin
 x1	=f(is-1,0)
 y1	=f(is-1,1)
 z1	=f(is-1,2)
 aia_mag_vector,x1,y1,z1,ds,hmin,hmax,coeff,curr,a_mis,polar,x2,y2,z2,b_tot,endv
 r2	=sqrt(x2^2+y2^2+z2^2)
 if (r2 le (1.+hmin)) or (r2 ge (1.+hmax)) or (is ge nfmax-1) or (endv eq 1) then $
  goto,first_footpoint
 f(is,0)=x2
 f(is,1)=y2
 f(is,2)=z2
 f(is,3)=b_tot
endfor
FIRST_FOOTPOINT:
nf1	=is < (nfmax-1)
nm	=nf1-1
if (nf1 ge 2) then begin
 f[nf1-1,3]=f[nf1-2,3]
 for j=0,3 do f[0:nf1-1,j]=reverse(f[0:nf1-1,j])
endif
polar	=-1.
for is=nf1,nfmax-1 do begin
 x1	=f(is-1,0)
 y1	=f(is-1,1)
 z1	=f(is-1,2)
 aia_mag_vector,x1,y1,z1,ds,hmin,hmax,coeff,curr,a_mis,polar,x2,y2,z2,b_tot,endv
 r2	=sqrt(x2^2+y2^2+z2^2)
 if (r2 lt (1.+hmin)) or (r2 gt (1.+hmax)) or (is ge nfmax-1) or (endv eq 1) then $
  goto,end_fieldline
 f(is,0)=x2
 f(is,1)=y2
 f(is,2)=z2
 f(is,3)=b_tot
endfor
END_FIELDLINE:
ns	=is < (nfmax-1)
field   =fltarr(ns,4)
for j=0,3 do field[*,j]=f[0:ns-1,j]
end