subr gen_corks, delta, itime, x, y
 ;assumes current grid in $gx, $gy already computed
 ;assumes delta is a mean flow field with 3 dimensions or
 ;a time dependent with 4, ignores dt in latter case
 ;computes positions up to itime assuming that flow displacements are in units of
 ;pixels per time step, if delta is time dependent, then itime must be less
 ;than the t dimension of delta
 ;x and y are returned in units of image pixels
 zeroifnotdefined, $worm_flag, $sparse_flag, $newrate
 ngx = dimen($gx,0)	ngy = dimen($gy,1)
 if num_dim(delta) eq 4 then {
 ;check if enough times in delta
   if dimen(delta,3) lt itime then { ty,'GEN_CORKS, itime too large'
 	x = 0	y = 0  return }
 }
 nt = itime
 ;now compute gx and gy to be just 1-D arrays for the x and y grid points
 gx = long($gx(*,0))		gy = long($gy(0,*))
 newtotal = $newrate*nt
 ty,'number of new corks = ',newtotal
 ncorks = ngx * ngy 
 if $sparse_flag le 1 then {
 x = fltarr(ncorks)		y = fltarr(ncorks)
 x(0) = $gx			y(0) = $gy
 } else {
 $sparse_flag = fix($sparse_flag)
 ncorks = fix((ngx+$sparse_flag-1)/$sparse_flag)*fix((ngy+$sparse_flag-1)/$sparse_flag)
 x = fltarr(ncorks)		y = fltarr(ncorks)
 d,x,y
 ic = 0
 for i=0,ngx-1,$sparse_flag do for j = 0,ngy-1,$sparse_flag do {
 x(ic) = $gx(i,j)
 y(ic) = $gy(i,j)
 ic += 1
 }
 }
 ind = lonarr(ncorks)	;for computed cork positions
 indlast = ind
 zero, indlast
 ;need some new random ones ?
 if newtotal gt 0 then {
 ty,'computing random new corks'
 iq = rfix( randomu(newtotal) * ($nx*$ny-1) )
 d, iq
 ty, 'max,min(iq)=', max(iq), min(iq)
 inq = indgen(lonarr($nx,$ny),0)
 newcork_x = float(inq(iq)
 inq = indgen(lonarr($nx,$ny),1)
 newcork_y = float(inq(iq)
 inq = 0
 indnew = lonarr(newtotal)
 indnewlast = lonarr(newtotal)
 zero, indnew, indnewlast
 }
 nx = $nx		ny = $ny
 nxm = nx -1		nym = ny -1
 t1 = !systime

 ;the mean flow case
 if num_dim(delta) eq 3 then {
 dx = delta(0,*,*)	dy = delta(1,*,*)
 for i = 0, nt-1 do {
 vx = bilinxy(dx,gx,gy,x,y)
 vy = bilinxy(dy,gx,gy,x,y)
 x = ((x + vx) > 0) < nxm		y = ((y +vy) > 0) < nym
 if newtotal gt 0 then {
 ;do the new ones separately
 new_count = (i+1) * $newrate - 1
 xq = newcork_x(0:new_count)
 yq = newcork_y(0:new_count)
 vx = bilinxy(dx,gx,gy,xq,yq)
 vy = bilinxy(dy,gx,gy,xq,yq)
 xq = ((xq + vx) > 0) < nxm
 newcork_x(0) = xq
 yq = ((yq +vy) > 0) < nym
 newcork_y(0) = yq
 }
 }

 ;the time dependent flow case
 } else {
 for i = 0, nt-1 do {
 vx = bilinxy(delta(0,*,*,i),gx,gy,x,y)
 vy = bilinxy(delta(1,*,*,i),gx,gy,x,y)
 x = ((x + vx) > 0) < nxm		y = ((y +vy) > 0) < nym

 if newtotal gt 0 then {
 ;do the new ones separately
 new_count = (i+1) * $newrate - 1
 xq = newcork_x(0:new_count)
 yq = newcork_y(0:new_count)
 vx = bilinxy(delta(0,*,*,i),gx,gy,xq,yq)
 vy = bilinxy(delta(1,*,*,i),gx,gy,xq,yq)
 xq = ((xq + vx) > 0) < nxm
 newcork_x(0) = xq
 yq = ((yq +vy) > 0) < nym
 newcork_y(0) = yq
 }
 }
 
 }
 ;now combine the originals and the new ones (if any)
 if newtotal gt 0 then {
 x = [x, newcork_x]
 y = [y, newcork_y]
 } 
 t2 = !systime
 ty,'time to compute and display cork positions =',t2-t1
 ty,'time per position =', (t2-t1)/nt

 endsubr
 ;===============================================================================