pro vpd,vx1,vy1,step,thr,im
;
; Displays the flow-map VX1,VY1 like arrows at each STEP-th point.
; THR is a maximum threshold for VX1,VY1. IM is a name of array
; containing the underlying image. Black background, white axes.
;
vx=vx1>(-thr)
vx=vx<thr
vy=vy1>(-thr)
vy=vy<thr
;
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px(1)-px(0)+1
sy=py(1)-py(0)+1
erase
;set_plot,'ps'
;device,filename='d:\data\palma\map.ps',/color
;loadct,3
tv,congrid(bytscl(im),sx,sy),px(0),py(0)
;tv,congrid(bytscl(im),sx,sy),0,0,xsize=12.7,ysize=12.7,/centimeters
;
mag=sqrt(vx^2+vy^2)
s=size(mag)
lmax=max(mag)
print,lmax,mean(mag)
vect=fltarr(s(1),s(2))
mul=1.5*step/lmax		; arrow length...
k=fix(step)
x0=0
x1=0
y0=0
y1=0
;!p.position=[0,0,0.7155,1]
;plot,[x0,x1],[y0,y1],/nodata,xrange=[0,s(1)],yrange=[0,s(2)],title='vel',$
;    xstyle=1,ystyle=1,/noerase,/normal,xthick=5,ythick=5,charthick=5
plot,[x0,x1],/nodata,xrange=[0,s(1)],yrange=[0,s(2)],xstyle=1,ystyle=1,$
/noerase
for i=k,s(1)-k,k do begin
  for j=k,s(1)-k,k do begin
    x0=i
    x1=i+mul*vx(i,j)*2
    angle=22.5*!dtor
    a1=abs(atan(vx(i,j)/vy(i,j)))
    st=mul*mag(i,j)*sin(angle)
    ct=mul*mag(i,j)*cos(angle)
    lte=sqrt(st^2+ct^2)/2.
      x2=x1-sin(angle+a1)*lte*vx(i,j)/abs(vx(i,j))
      x3=x1-sin(a1-angle)*lte*vx(i,j)/abs(vx(i,j))
    y0=j
    y1=j+mul*vy(i,j)*2
      y2=y1-cos(a1+angle)*lte*vy(i,j)/abs(vy(i,j))
      y3=y1-cos(a1-angle)*lte*vy(i,j)/abs(vy(i,j))
;      plots,[x0,x1,x2,x1,x3],[y0,y1,y2,y1,y3],color=235,thick=4
     plots,[x0,x1,x2,x3,x1],[y0,y1,y2,y3,y1],color=0,thick=1
    endfor
endfor
;div=divergence(vx,vy)
;lh=max(div)/2.
;l1=lh/5.
;contour,div,xstyle=5,ystyle=5,/noerase,$
;       nlevels=10,color=235,/normal,c_thick=1
;device,/close
;set_plot,'x'
end