pro fieldlines,b1,b2,b3,factor,angx=angx,angz=angz,ps=ps
;
; traces field lines for a specified magnetic field
; configuration.  Allows for a non-uniform grid spacing. 
; Uses the bi-linear & tri-linear interpolation routines
; lin2interp and lin3interp, as well as the routine
; differential2.
;
if(n_params() lt 4) then begin
  print,'fieldlines,bx,by,bz,factor (nfl=number of field lines)'
  print,'angx=angx,angz=angz,ps=ps'
  return
endif
;
if (n_elements(ps) ne 0) then begin 
  set_plot,'ps'
  device,filename='out.ps',bits=8,/inches,xoffset=0.7,yoffset=2.7,$
  xsize=7.,ysize=7.
endif
if (n_elements(angx) eq 0) then begin
  angx=30.
endif
if (n_elements(angz) eq 0) then begin
  angz=30.
endif
;
common shared, bx,by,bz,x,y,z,x2,y2,z2,sign
bx=b1 & by=b2 & bz=b3
;
; begin calculation
;
nx=n_elements(bx(*,0,0))
ny=n_elements(bx(0,*,0))
nz=n_elements(bx(0,0,*))
x=findgen(nx)/float(nx-1) & y=findgen(ny)/float(ny-1) 
z=findgen(nz)/float(nz-1)
magb=sqrt(bx^2+by^2+bz^2)
;
; define grid 
;
x2=0.5d+0*x+0.5d+0
y2=0.5d+0*y+0.5d+0
z2=0.5d+0*z
magfac=2.*float(nx/ny)
x2=magfac*0.5d+0*x
y2=magfac*0.5d+0*y
z2=magfac*0.5d+0*z
xmin=min(x2) & xmax=max(x2) & ymin=min(y2) & ymax=max(y2) & zmin=min(z2)
zmax=max(z2)
;
; define grid spacing
;
dx=dblarr(nx-1) & dy=dblarr(ny-1) & dz=dblarr(nz-1)
for i=0,nx-2 do begin
  dx(i)=x(i+1)-x(i)
endfor
for j=0,ny-2 do begin
  dy(j)=y(j+1)-y(j)
endfor
for k=0,nz-2 do begin
  dz(k)=z(k+1)-z(k)
endfor
;
; open square window (otherwise box is distorted)
;
if (n_elements(ps) eq 0) then begin
    window,xsize=512,ysize=512
endif
;
; set 3D coordinate system
;
scale3,xrange=[0,1],yrange=[0,1],zrange=[0,1],$
ax=angx,az=angz
;
; specify color table (rainbow + white)
;
if (n_elements(ps) ne 0) then begin
    loadct,0
endif else begin
    loadct,0
endelse
;
if (n_elements(ps) ne 0) then begin 
    white=0 
    flcol=0
endif else begin
    white=255
    flcol=255
endelse
;
; draw image of bz in boundary
;
print,'regridding and drawing image..... '
xorig=[0,nx-1,0,nx-1]
yorig=[0,0,ny-1,ny-1]
s=size(bz(*,*,0))
xp=(xorig*!x.s(1)+!x.s(0))/(1.*(nx-1))
yp=(yorig*!y.s(1)+!y.s(0))/(1.*(ny-1))
p=[[xp],[yp],[fltarr(4)],[replicate(1,4)]]#!P.T
;
u=p(*,0)/p(*,3)*!d.x_vsize
v=p(*,1)/p(*,3)*!d.y_vsize
u0=min(u) & v0=min(v)
su=max(u)-u0+1 & sv=max(v)-v0+1
;
if (!d.flags and 1) eq 1 then begin     ;Scalable pixels (PostScript)?
    fact = 50               ;Yes, shrink it
    miss = 255              ;Missing values are white
    c_color=[0,0]           ;Contour in only one
endif else begin
    fact = 1                ;one pixel/output coordinate
    miss = 0                ;missing is black
    c_color=[150,200,250]
endelse
;
if (!d.flags and 512) ne 0 then $  ;White background?
    miss = 255 else miss = 0       ;Get polynomial coeff for warp 
if !d.n_colors gt 2 then top = !d.n_colors -1 else top = 255
;
u2=(u-u0)/fact
v2=(v-v0)/fact
m=[[u2(0),v2(0),0,0,1,0],[0,0,u2(0),v2(0),0,1],[u2(1),v2(1),0,0,1,0],$
[0,0,u2(1),v2(1),0,1],[u2(2),v2(2),0,0,1,0],[0,0,u2(2),v2(2),0,1]]
rhs=[xorig(0),yorig(0),xorig(1),yorig(1),xorig(2),yorig(2)]
;
ludc,m,index,/double
sol=lusol(m,index,rhs,/double)
;
kx=fltarr(2,2)
ky=fltarr(2,2)
kx(0,1)=sol(0)
kx(1,0)=sol(1)
ky(0,1)=sol(2)
ky(1,0)=sol(3)
kx(0,0)=sol(4)
ky(0,0)=sol(5)
;
interp = 2
;
; rebin boundary values on a uniform grid
;
image=dblarr(nx,ny)
;
xu=indgen(nx)/(1.0*(nx-1))
yu=indgen(ny)/(1.0*(ny-1))
;
for i=0,nx-1 do begin
  for j=0,ny-1 do begin
    image(i,j)=lin2interp(bz(*,*,0),x2,y2,xu(i),yu(j))
  endfor
endfor
;
a=poly_2d(bytscl(image,top=top),kx,ky,interp,su/fact,sv/fact,$
  missing=miss) ;Warp it
tv,a,u0,v0,xsize=su,ysize=sv,/device
;
; draw axes
;
plots,[xmin,xmax],[ymin,ymin],[zmin,zmin],/t3d,/data,color=white,$
thick=2.0 ; x-axis
plots,[xmin,xmin],[ymin,ymax],[zmin,zmin],/t3d,/data,color=white,$
thick=2.0 ; y-axis
plots,[xmin,xmin],[ymin,ymin],[zmin,zmax],/t3d,/data,color=white,$
thick=2.0 ; y-axis
;
; outline box by drawing edges (dashed)
;
plots,[xmax,xmax],[ymin,ymax],[zmin,zmin],/t3d,/data,linestyle=2,color=white,$
thick=2.0
plots,[xmax,xmin],[ymax,ymax],[zmin,zmin],/t3d,/data,linestyle=2,color=white,$
thick=2.0
;
; label axes
;
xyouts,xmax+.05,ymin,z=zmin,'x',/t3d,/data,size=3,color=white,$
charthick=2.0
xyouts,xmin,ymax+.05,z=zmin,'y',/t3d,/data,size=3,color=white,$
charthick=2.0
xyouts,xmin,ymin,z=zmax+0.05,'z',/t3d,/data,size=3,color=white,$
charthick=2.0
;
print,'drawing fieldlines...'
;
; field line tracing - loop over starting coords
;
bzmax=max(reform(bz(*,*,0)))/factor
step=6
for ii=0,nx-1,step do begin
  iy=where( ( abs(bz(ii,*,0)) gt bzmax ),icounty)
for jj=0,ny-1,step do begin
  ix=where( (reform(abs(bz(*,jj,0))) gt bzmax ),icountx)
for i=0,icountx-1 do begin
  for j=0,icounty-1 do begin
    x0=x2(ix(i)) & y0=y2(iy(j)) & z0=zmin
       ; h is stepsize, by default 1/8 of minimum gridsize
    h=min(dx)/8.0 & s=0.0 & coords=[x0,y0,z0]
       ; determine whether field is up or down
    sign=float(0)
    bznearby=lin3interp(bz,x2,y2,z2,x0,y0,z0)
       ; always integrate into the box
    if (bznearby gt 0.0) then sign=1.0 else sign=-1.0 
       ; integrate until field line leaves box
    while (coords(0) ge xmin) and (coords(0) le xmax) and $
        (coords(1) ge ymin) and (coords(1) le ymax) and $
        (coords(2) ge zmin) and (coords(2) le zmax) do begin
         ; RHS of ODEs determined by linear interpolation
      dydx=differential2(s,coords)
      result=rk4(coords,dydx,s,h,'differential2')
         ; plot increment in field line
      plots,[coords(0),result(0)],[coords(1),result(1)],$
        [coords(2),result(2)],/t3d,/data,thick=1,color=flcol
      coords=result
    endwhile
  endfor
endfor
endfor
print,'completed loop ',ii
endfor
;
if (n_elements(ps) ne 0) then begin
  device,/close
  set_plot,'x
endif
;
return
end