pro fieldlines,b1,b2,b3,alpha,sfact,bmin=bmin,istep=istep,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 5) then begin
  print,'fieldlines,bx,by,bz,alpha,sfact,bmin=bmin,istep=istep,' 
  print,'angx=angx,angz=angz,ps=ps'
  return
endif
;
!p.background=255
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
if (n_elements(bmin) eq 0) then begin
  bmin=100.
endif
if (n_elements(istep) eq 0) then begin
  istep=2
endif
;
common shared, bx,by,bz,x,y,z,x2,y2,z2,sign
bx=b1 & by=b2 & bz=b3
;
; begin calculation
;
!p.background=255
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=800,ysize=800
;    window,xsize=400,ysize=400
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
;
set_plot,'z'
!p.background=255
loadct,0
erase
device,set_resolution=[800,800]
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 
    miss = 255 else miss = 255       ;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)]
rhs=[xorig(0),yorig(0),xorig(1)+1,yorig(1),xorig(2),yorig(2)+1]
;
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
;
; outline box by drawing edges (dashed)
;
plots,[xmax,xmax],[ymin,ymax],[zmin,zmin],/t3d,/data,linestyle=0,color=0,$
thick=3.0
plots,[xmax,xmin],[ymax,ymax],[zmin,zmin],/t3d,/data,linestyle=0,color=0,$
thick=3.0
;
print,'drawing fieldlines...'
;
; field line tracing - loop over starting coords
;
bzmax=max(reform(bz(*,*,0)))
alphamax=max(alpha)
alphamin=min(alpha)
avgalpha=total(alpha)/n_elements(alpha)
stdalpha=stddev(alpha)/sfact
bzmin=min(reform(bz(*,*,0)))
nbins=10
;istep=2
;   istep=6
istepmin=2
if (bzmax le 0.) then begin
  print,'bzmax lt 0.'
  return
endif
dbz=(bzmax-bzmin)/nbins
bin=fltarr(nbins+1)
bin(0)=bzmin
for i=1,nbins do begin 
  bin(i)=bin(i-1)+dbz
endfor
for n=1,nbins do begin
  bzb=bytscl(reform(abs(bz(*,*,0))),min=0.,max=bzmax,top=254)
;  alphab=bytscl(alpha,min=avgalpha-stdalpha,max=avgalpha+stdalpha,top=254)
  alphab=bytscl(abs(alpha),min=0.,max=avgalpha+stdalpha,top=254)
  index=where( ( (reform(bz(*,*,0)) ge bin(n-1)) and $
    (reform(bz(*,*,0)) lt bin(n)) and (abs(reform(bz(*,*,0))) gt bmin) ),jcount)
  iy=fix(index/ny)
  ix=index-iy*ny
  icount=n_elements(ix)
  for i=0,icount-1,istep do begin
    x0=x2(ix(i)) & y0=y2(iy(i)) & z0=zmin
      ; h is stepsize, by default 1/8 of minimum gridsize
    h=min(dx)/28.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
;    flcol=bzb(ix(i),iy(i))
    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')
      ixx=fix(nx*(coords(0)+result(0))/2.)
      iyy=fix(ny*(coords(1)+result(1))/2.)
      izz=fix(nz*(coords(2)+result(2))/2.)
      if(ixx ge nx) then ixx=nx-1
      if(iyy ge ny) then iyy=ny-1
      if(izz ge nz) then izz=nz-1
      flcol=alphab(ixx,iyy,izz)
      ; plot increment in field line
      plots,[coords(0),result(0)],[coords(1),result(1)],$
        [coords(2),result(2)],/t3d,/data,thick=4.0,color=flcol
      coords=result
    endwhile
  endfor
  istep=istep-1
  if(istep lt istepmin) then istep=istepmin
  print,'finished bin ',n
endfor
;
;  color bar
;
;  note: run cbar.jou first
;
image=tvrd()
set_plot,'x'
tv,image,xsize=su,ysize=sv,/device
cbox=0
psx0=0.35
psy0=0.
psx1=0.
psy1=0.
cbarchar=1.75
;cbar,vmin=0.,vmax=bzmax,cmin=1,cmax=255,/horizontal,$
;  position=[0.56-psx0,0.87-psy0,0.85-psx1,0.93-psy1],$
;  color=cbox,charsize=cbarchar,charthick=2.5
if (n_elements(ps) ne 0) then begin
  device,/close
  set_plot,'x
endif
;
;  annotations
;
loadct,0
;xyouts,0.07,0.89,'!3B!dz!n(0)',charsize=1.8,color=0,charthick=2.5
;
; draw axes
;
plots,[xmin,xmax],[ymin,ymin],[zmin,zmin],/t3d,/data,color=0,$
thick=3.0 ; x-axis
plots,[xmin,xmin],[ymin,ymax],[zmin,zmin],/t3d,/data,color=0,$
thick=3.0 ; y-axis
plots,[xmin,xmin],[ymin,ymin],[zmin,zmax],/t3d,/data,color=0,$
thick=3.0 ; y-axis
;
; label axes
;
charthick=1.75
;xyouts,xmax+.05,ymin,z=zmin,'X',/t3d,/data,size=4,color=0,$
charthick=charthick
;xyouts,xmin,ymax+.05,z=zmin,'Y',/t3d,/data,size=4,color=0,$
charthick=charthick
;xyouts,xmin,ymin,z=zmax+0.05,'Z',/t3d,/data,size=4,color=0,$
charthick=charthick
;
return
end