pro bgvect, vx, vy, bgarr=bgarr, overplot=overplot, help=help, $
            bgct=bgct, vct=vct, vcolor=vcolor, aspect=aspect, $ 
            arrwfac=arrwfac, arrwmod=arrwmod, xmod=xmod, ymod=ymod,$
            xtitle=xtitle, ytitle=ytitle, title=title, min=min, max=max, $
            missing=missing, $
            scale=scale, scolor=scolor, scstr=scstr, sci=sci, sct=sct, $
            charsize=charsize,thick=thick, SCLTH=SCLTH,  $
            xfrac=xfrac,yfrac=yfrac, xs=xs, ys=ys,vmax=vmax, $
            axcol=axcol, plot_charsize=plot_charsize, $
            position = position, abs_scl=abs_scl

if (keyword_set(help)) then begin
    print,'%'
    print,'% PRO BGVECT: VELOVECTs vector field on DISPLAYed background.'
    print,'%'
    print,'% USAGE: bgvect,vx,vy,bgarr=bgarr,arrwfac=2,arrwmod=3,vct=12,vcol=108'
    print,'%'
    print,'% PURPOSE: VELOVECTS reduced density vector field, over BGARR.' 
    print,'%  Can layer successive vectors by calls with OVERPLOT set.'
    print,'%  Can XYOUTS vector scale at position [XFRAC, YFRAC]'
    print,'%'
    print,'% REQUIRES:' 
    print,'%'
    print,'%  OTHER ROUTINES: reduce.pro and sig_dig.pro'
    print,'%  VX = array of x-comps. of vector field.'
    print,'%  VY = array of y-comps. of vector field.'
    print,'%  BGARR = DISPLAYed as background (REQUIRED if OVERPLOT not set)'
    print,'%'
    print,'% OPTIONAL ARGS:' 
    print,'%'
    print,'%  OVERPLOT = set this to suppress DISPLAY of bgarr'
    print,'%  BGCT = sets background color table, default 0'
    print,'%  X-/Y-/TITLE = axes title strings sent display if not /overplot, default ""'
    print,'%  XS/YS = 1-d axis arrays, passed to DISPLAY'
    print,'%'
    print,'%  ARRWMOD = REBINs to REDUCE vector density modulo arrwmod, default 4'
    print,'%  ARRWFAC = sets length for VELOVECT, default = 1.'
    print,'%'
    print,'%  MIN/MAX = passed to DISPLAY'
    print,'%  ASPECT = passed to DISPLAY, defaults to 1 unless set to 0'
    print,'%  MISSING = passed to VELOVECT'
    print,'%  '
    print,'%  VCT = sets vector color table, default 0'
    print,'%  VCOLOR = sets vector color, default 0'
    print,'% ' 
    print,'%  SCALE = set to show vector scale, at [XFRAC, YFRAC]'
    print,'%  VMAX = max. length for scaling vectors; useful with /OVER'
    print,'%  XFRAC = default .05 (left)'
    print,'%  YFRAC = default .9 (upper)'
    print,'%  SCT = if SCALE, sets color table, default 0'
    print,'%  SCOLOR = sets color for scale text, default 0'
    print,'%  SCSTR = if scale, string for units, default ""'
    print,'%  SCI = if scale, makes scale number in sci. notation'
    print,'%  CHARSIZE = passed to XYOUTS, default 1'
    print,'%  THICK = passed to velovect'
    print,'%  SCLTH = thickness passed to XYOUTS, default 2'
    print,'%  HELP = prints this documentation, then exits.'
    print,'%'
    print,'% WARNING: Error handling is primitive -- halts mid-procedure.'
    print,'%'
    print,'% HISTORY: Written Brian Welsch, Dec. 2003'
    print,'%    BTW: Fixed bug in scaling, 21 Apr. 2005'
    print,'%    BTW: Mucked w/XS,YS keywords, 07 Jun. 2005'
    print,'%    BTW: Only loads new ct if (ct ne 0), 23 Jan 2006.'
    print,'%    BTW: Fixed bug with scale placement, 25 Jan 2006.'
    print,'%    BTW: Added compatibility with nonzero !p.multi'
    print,'%    BTW: Last Modified 10 May 2006.'
    print,'%'
    return
endif


if (n_params() eq 0) then begin
    print,'%'
    print,'% PRO BGVECT: VELOVECTs vector field on DISPLAYed background.'
    print,'%'
    print,'% USAGE: bgvect,vx,vy,bgarr=bgarr,arrwfac=2,arrwmod=3,vct=12,vcol=108'
    print,'%'
    print,'% DOCUMENTATION: call with /HELP set to see more information.'
    print,'%'
    print,'% PURPOSE: VELOVECTS reduced density vector field, over BGARR. '
    print,'%  Can layer successive vectors by calls with OVERPLOT set.'
    print,'%  Can XYOUTS vector scale at position [XFRAC, YFRAC]'
    print,'%'
    print,'% REQUIRED: '
    print,'%  OTHER ROUTINES: reduce.pro and sig_dig.pro'
    print,'%  VX = array of x-comps. of vector field.'
    print,'%  VY = array of y-comps. of vector field.'
    print,'%  BGARR = DISPLAYed as background (REQUIRED if OVERPLOT not set)'
    print,'%'
    print,'% HISTORY: Written Brian Welsch, Dec. 2003'
    print,'%    BTW: Last Modified 10-May-2006.'
    print,'%'
    return
endif 


nx = n_elements(vx(*,0))
ny = n_elements(vx(0,*))

if not(keyword_set(arrwfac)) then arrwfac = 1. 
arrwfac_in = arrwfac

if not(keyword_set(arrwmod)) then begin
    if not(keyword_set(xmod)) then xmod = 4. 
    if not(keyword_set(ymod)) then ymod = 4. 
    arrwmod = sqrt(xmod^2 + ymod^2)
endif
if not(keyword_set(xmod)) then xmod = arrwmod
if not(keyword_set(ymod)) then ymod = arrwmod

if not(keyword_set(filename)) then filename = 'bgvect.ps'

if (n_elements(bgct) eq 0) then bgct = 0
if (n_elements(vct) eq 0) then vct = 0
if (n_elements(vcolor) eq 0) then vcolor = 0
if (n_elements(sct) eq 0) then sct = 0
if (n_elements(scolor) eq 0) then scolor = 0
if (n_elements(aspect) eq 0) then aspect = 1


; IF NOT OVERPLOT, DISPLAY BGARR
if not(keyword_set(overplot)) then begin
    dims = size(bgarr)
    if not(keyword_set(xs)) then xs = findgen(dims(1))
    if not(keyword_set(ys)) then ys = findgen(dims(2))
    if not(keyword_set( title)) then  title = ''
    if not(keyword_set(xtitle)) then xtitle = ''
    if not(keyword_set(ytitle)) then ytitle = ''

    if (!p.multi(1) ne 0) or (!p.multi(2) ne 0) then begin
 
       if (keyword_set(bgct)) then loadct,bgct,/silent

       tvscale, bgarr,keep_aspect_ratio=aspect,position=position

       if (keyword_set(bgct)) then loadct,0,/silent

       plot,[min(xs),max(xs)],[min(ys),max(ys)],xs=1,ys=1, $
         /nodata,/noerase, $
         position=temporary(position), $
         xtit=xtitle, ytit=ytitle, tit=title, $
         charsize = plot_charsize, color = axcol

       no_xs = temporary(xs)
       no_ys = temporary(ys)

    endif else begin
        if (keyword_set(bgct)) then loadct,bgct,/silent
        display,bgarr,xs,ys,min=min,max=max,$
          tit=title,xtit=xtitle,ytit=ytitle,aspect=aspect, $
          axcol=axcol;, charsize = plot_charsize
    endelse
    if (keyword_set(bgct)) then loadct,0,/silent
endif

; PLOT VECTORS
if (keyword_set(vct)) then loadct,vct,/silent

; REDUCE VECTOR DENSITY, IF CALLED FOR
if (arrwmod gt 1) then begin
    rvx = reduce(vx, xmod, ymod, sx=sx, sy=sy) ; sx,sy are "reduced" axis arrays
    rvy = reduce(vy, xmod, ymod, sx=sx, sy=sy)

    if (not(keyword_set(xs))) then begin
        xsmin = !x.crange(0)
        xsmax = !x.crange(1)
        dxs = xsmax - xsmin
        ysmin = !y.crange(0)
        ysmax = !y.crange(1)
        dys = ysmax - ysmin

        rnx = n_elements(rvx(*,0))
        rny = n_elements(rvx(0,*))

        rxfrac = rnx*xmod/float(nx)
        ryfrac = rny*ymod/float(ny)

        rxs = dxs*(findgen(rnx)+.5)/(rnx)*rxfrac + xsmin
        rys = dys*(findgen(rny)+.5)/(rny)*ryfrac + ysmin
    endif else begin
        rxs = xs(sx)
        rys = ys(sy)        
    endelse

    if keyword_set(abs_scl) then begin
        if not(keyword_set(missing)) then missing = 0.
        notmiss = where((rvx ne missing) and (rvy ne missing), n_notmiss)
        vmax = max(sqrt(rvx(notmiss)^2+rvy(notmiss)^2)) 
        arrwfac = vmax/arrwfac_in
    endif 

    velovect, rvx, rvy, rxs, rys, /overplot,$  
      length=arrwfac, color=vcolor, missing=missing, thick=thick
endif

fail_check0 = check_math() 
if ((fail_check0 ne 0) and (fail_check0 ne 32)) then begin
    print,"Check Math has detected a non-trivial error."
    print,fail_variable0
endif


; OR PLOT AS IS!
if (arrwmod eq 1) then begin
    rvx = vx 
    rvy = vy 
    sx = findgen(nx)
    sy = findgen(ny)
    rxs = sx 
    rys = sy 
    if keyword_set(abs_scl) then begin
        if not(keyword_set(missing)) then missing = 0.
        notmiss = where((rvx ne missing) and (rvy ne missing), n_notmiss)
        vmax = max(sqrt(rvx(notmiss)^2+rvy(notmiss)^2)) 
        arrwfac = vmax/arrwfac_in
    endif 
    velovect, rvx, rvy, /overplot,$  
      length=arrwfac, color=vcolor, missing=missing, thick=thickness
endif  

fail_check0_5 = check_math() 
if ((fail_check0_5 ne 0) and (fail_check0_5 ne 32)) then begin
    print,"Check Math has detected a non-trivial error."
    print,fail_variable0_5
endif

if not(keyword_set(vmax)) then begin
    if not(keyword_set(missing)) then $
      vmax = max(sqrt(rvx^2+rvy^2)) else begin 
        notmiss = where((rvx ne missing) and (rvy ne missing), n_notmiss)
        if (n_notmiss eq 0) then begin
            print,'All elements equal "missing" value. Returning.'
            return
        endif
        vmax = max(sqrt(rvx(notmiss)^2+rvy(notmiss)^2)) 
    endelse 
endif 


fail_check1 = check_math() 
if ((fail_check1 ne 0) and (fail_check1 ne 32)) then begin
    print,"Check Math has detected a non-trivial error."
    print,fail_variable1
endif

if (keyword_set(vct)) then loadct,0,/silent


; PRINT SCALE ON PLOT?
if keyword_set(scale) then begin

    if not(keyword_set(scstr)) then scstr = ''
    if not(keyword_set(thickness)) then thickness = 1
    
    if (keyword_set(xs)) then minx = min(xs) else minx = !x.crange(0)
    if (keyword_set(xs)) then maxx = max(xs) else maxx = !x.crange(1)
    deltax = maxx-minx

    if (keyword_set(ys)) then miny = min(ys) else miny = !y.crange(0)
    if (keyword_set(ys)) then maxy = max(ys) else maxy = !y.crange(1)
    deltay = maxy-miny

    nsx = n_elements(sx)
    nsy = n_elements(sy)
    dx = (sx(1)-sx(0))/nx*deltax
    dy = (sy(1)-sy(0))/ny*deltay
    
    if (keyword_set(sct)) then loadct, sct,/silent
    
    if not(keyword_set(xfrac)) then xfrac =.05
    if not(keyword_set(yfrac)) then yfrac =.9
    
    plots,$ ; horizontal bar
      [xfrac*deltax - .5*dx, xfrac*deltax + .5*dx]+minx, $
      [yfrac*deltay + .5*dy, yfrac*deltay + .5*dy]+miny, $
      thick=sclth, color=scolor
    plots, $ ; left vertical bar
      [xfrac*deltax - .5*dx, xfrac*deltax - .5*dx]+minx, $
      [yfrac*deltay        , yfrac*deltay +    dy]+miny, $
      thick=sclth, color=scolor
    plots, $ ; right vertical bar
      [xfrac*deltax + .5*dx, xfrac*deltax + .5*dx]+minx, $
      [yfrac*deltay        , yfrac*deltay +    dy]+miny, $
      thick=sclth, color=scolor
    
    if keyword_set(abs_scl) then len = arrwfac_in else len = vmax/arrwfac
    len = sig_dig(len,2,sci=sci, /plot) 
    xyouts, xfrac*deltax + 1.5*dx+minx, $
      yfrac*deltay + miny, $
      len+' '+scstr, color=scolor, charsize=charsize, charthick=sclth    

    if (keyword_set(sct)) then loadct, 0,/silent
endif 

; Error trap, ignore floating underflows
fail_check = check_math() 
if ((fail_check ne 0) and (fail_check ne 32)) then print,fail_variable

arrwfac = arrwfac_in

end