;cork routines
 ;some routines to handle the 2/23/96 flowmaps
 ;=================================================================
 ;use more advanced versions of pixelmaskarea and xdraw here, eventually
 ;put in libraries and update browser
 ;============================================================================
subr pixelmaskarea(pw, window, in, nx, ny)
 ;pw is the pen width to use for the drawn mask
 ;returns indices of area masked
 pencolor,'red'
 ;ty,'to use: put cursor on start position, click quickly once, and move'
 ;ty,'mouse with button UP (this is easier for detailed work).'
 ;ty,'Click again to stop.'
 ;ty,'Do NOT hold button down or it will just give you a dot.'
 xdraw, pw, window
 ;drawn, read back and figure mask
 xq=tvread(window)
 nx = dimen(xq,0)	ny = dimen(xq,1)
 in=sieve(xq eq iq)
 ;now remark the area in green so it doesn't interfere with next area
 xq(in)=!green_pixel
 tv,xq
 endsubr
 ;============================================================================
subr show_corks, delta, ntformean
 ;also does colored worms, using current color table, if $worm_flag is set
 ;assumes current grid in $gx, $gy already computed
 ;assumes delta is a mean flow field with 3 dimensions or
 ;a time dependent with 4, if nt not present, uses dimension
 narg = !narg
 zeroifnotdefined, $worm_flag, $sparse_flag, $newrate, $display_scale, #xwin
 im = 0
 if $display_scale eq 0 then $display_scale = 1.0
 ngx = dimen($gx,0)	ngy = dimen($gy,1)
 if num_dim(delta) eq 4 then {
   nt = dimen(delta,3)
   if narg ge 2 then { if ntformean le nt then nt = ntformean else {
 	  ty,'warning, specified nt larger than time dimension of map'
	  ty,'only the first',nt,' will be shown' }
	  }
 } else { nt = ntformean }
 ;now compute gx and gy to be just 1-D arrays for the x and y grid points
 gx = $display_scale*long($gx(*,0))
 gy = $display_scale*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) = $display_scale*$gx	y(0) = $display_scale*$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)
   y(ic) = gy(j)
   ic += 1
   }
 }
 ty, 'ncorks =', ncorks
 ind = lonarr(ncorks)	;for computed cork positions
 indlast = ind
 zero, indlast
 nx = fix($nx*$display_scale)
 ny = fix($ny*$display_scale)
 ty,'nx, ny, nx*ny =', nx,ny,nx*ny
 xwin,#xwin,nx,ny
 nxm = nx -1		nym = ny -1
 if $movie_flag then {
   fname = ''
   $gif_color = xtvrct(0)
   read,'enter gif file template', fname
   ;;read,'enter animated gif file template', fname
   ;;gif3dopenw, fname, nx, ny, $gif_color
 }
 
 brim = zero(bytarr(nx, ny))+byte(!black_pixel)
 ;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
   ;show the randoms
   iq = fix(newcork_x) + nx * fix(newcork_y)
   brim(iq) = !white_pixel
   tv, reverse(brim,1), 0, 0, #xwin
   mess=''
   read,'enter return to go on', mess
   brim(iq) = !black_pixel
 }
 t1 = !systime

 ;show the initial corks
 ind(0) = fix(x) + nx * fix(y)
 brim(ind) = !white_pixel
 switch, ind, indlast
 tv, reverse(brim,1), 0, 0, #xwin
 
 if nt le 0 then return

 ;the mean flow case
 if num_dim(delta) eq 3 then {
   dx = $display_scale*delta(0,*,*)	dy = $display_scale*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
     ind(0) = fix(x) + nx * fix(y)
     if $worm_flag eq 0 then brim(indlast) = !black_pixel
     brim(ind) = !white_pixel
     switch, ind, indlast
     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
       iq = fix(xq) + nx * fix(yq)
       indnew(0) = iq
       if $worm_flag eq 0 then brim(indnewlast(0:new_count)) = !black_pixel
       brim(iq) = !white_pixel
       switch, indnew, indnewlast
     }
     tv, reverse(brim,1), 0, 0, #xwin
   }

 ;the time dependent flow case
 } else {
   ;set up a pixel color range
   ncol = !scalemax - !scalemin + 1
   pixv = fltarr(ncol)
   pixv = indgen(pixv) + !scalemin
   pixf = polate(indgen(pixv), pixv, (indgen(fltarr(nt)))*ncol/float(nt) )
   for i = 0, nt-1 do {
     vx = $display_scale*bilinxy(delta(0,*,*,i),gx,gy,x,y)
     vy = $display_scale*bilinxy(delta(1,*,*,i),gx,gy,x,y)
     x = ((x + vx) > 0) < nxm		y = ((y +vy) > 0) < nym
     ind(0) = fix(x) + nx * fix(y)
     if $worm_flag eq 0 then {
       brim(indlast) = !black_pixel
       brim(ind) = !white_pixel
      } else {
       brim(ind) = pixf(i)
      }
     switch, ind, indlast

     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 = $display_scale*bilinxy(delta(0,*,*,i),gx,gy,xq,yq)
       vy = $display_scale*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
       iq = fix(xq) + nx * fix(yq)
       indnew(0) = iq
       if $worm_flag eq 0 then {
	 brim(indnewlast(0:new_count)) = !black_pixel
	 brim(iq) = !white_pixel
       } else {
	 brim(ind) = pixf(i)
       }
       switch, indnew, indnewlast
     }

     tv, reverse(brim,1), 0, 0, #xwin
     if $movie_flag then {
       gifwrite, reverse(brim,1), fns(fname, im), $gif_color
       ;;gif3dappend, reverse(brim,1)
       im = im + 1
     }
   }
 
 }
 ;;if $movie_flag then gif3dclose
 t2 = !systime
 ty,'time to compute and display cork positions =',t2-t1
 ty,'time per position =', (t2-t1)/nt
 endsubr
 ;===============================================================================
subr show_seeded_corks, seedx, seedy, delta, ntformean
 ;a variation that starts with specified cork positions of 2 flavors
 ;in seed1 and seed2
 ;assumes current grid in $gx, $gy already computed
 ;assumes delta is a mean flow field with 3 dimensions or
 ;a time dependent with 4, if nt not present, uses dimension
 narg = !narg
 zeroifnotdefined, $worm_flag, $sparse_flag, $newrate
 ngx = dimen($gx,0)	ngy = dimen($gy,1)
 if num_dim(delta) eq 4 then {
 nt = dimen(delta,3)
 if narg ge 4 then { if ntformean le nt then nt = ntformean else {
 	ty,'warning, specified nt larger than time dimension of map'
	ty,'only the first',nt,' will be shown' }
	}
 } else { nt = ntformean }
 ;now compute gx and gy to be just 1-D arrays for the x and y grid points
 gx = $display_scale*long($gx(*,0))
 gy = $display_scale*long($gy(0,*))
 newtotal = $newrate*nt
 ty,'number of new corks = ',newtotal
 
 
 ncorks = num_elem(seedx) 
 x = float(seedx)		y = float(seedy)


 ind = long(seedx)	;for computed cork positions
 d, ind
 indlast = ind
 zero, indlast
 nx = fix($nx*$display_scale)
 ny = fix($ny*$display_scale)
 ty,'nx, ny, nx*ny =', nx,ny,nx*ny
 xwin,0,nx,ny
 nxm = nx -1		nym = ny -1
 brim = zero(bytarr(nx, ny))+byte(!black_pixel)
 ;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
 ;show the randoms
 iq = fix(newcork_x) + nx * fix(newcork_y)
 brim(iq) = !white_pixel
 tv, brim, 0, 0, #xwin
 mess=''
 read,'enter return to go on', mess
 brim(iq) = !black_pixel
 }
 t1 = !systime

 ;the mean flow case
 if num_dim(delta) eq 3 then {
 dx = $display_scale*delta(0,*,*)	dy = $display_scale*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
 ind(0) = fix(x) + nx * fix(y)
 if $worm_flag eq 0 then brim(indlast) = !black_pixel
 brim(ind(*,0)) = !red_pixel
 brim(ind(*,1)) = !green_pixel
 switch, ind, indlast
 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
 iq = fix(xq) + nx * fix(yq)
 indnew(0) = iq
 if $worm_flag eq 0 then brim(indnewlast(0:new_count)) = !black_pixel
 brim(iq) = !white_pixel
 switch, indnew, indnewlast
 }
 tv, brim, 0, 0, #xwin
 }

 ;the time dependent flow case
 } else {
 for i = 0, nt-1 do {
 vx = $display_scale*bilinxy(delta(0,*,*,i),gx,gy,x,y)
 vy = $display_scale*bilinxy(delta(1,*,*,i),gx,gy,x,y)
 x = ((x + vx) > 0) < nxm		y = ((y +vy) > 0) < nym
 ind(0) = fix(x) + nx * fix(y)
 if $worm_flag eq 0 then brim(indlast) = !black_pixel
 
 brim(ind(*,0)) = !red_pixel
 brim(ind(*,1)) = !green_pixel
 switch, ind, indlast

 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 = $display_scale*bilinxy(delta(0,*,*,i),gx,gy,xq,yq)
 vy = $display_scale*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
 iq = fix(xq) + nx * fix(yq)
 indnew(0) = iq
 if $worm_flag eq 0 then brim(indnewlast(0:new_count)) = !black_pixel
 brim(iq) = !white_pixel
 switch, indnew, indnewlast
 }
 
 tv, brim, 0, 0, #xwin
 }
 
 }
 t2 = !systime
 ty,'time to compute and display cork positions =',t2-t1
 ty,'time per position =', (t2-t1)/nt

 endsubr
 ;===============================================================================
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
 ;===============================================================================
subr gen_seeded_corks, seedx, seedy, delta, itime, x, y, newcork_x, newcork_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
 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 = num_elem(seedx) 
 x = float(seedx)		y = float(seedy)

 ind = long(seedx)	;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
 } else { newcork_x = 0		newcork_y = 0 }
 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 news (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
 ;===============================================================================
subr getmeanv(window, gind, in)
 ;assumes $gx and $gy defined and that a flow map is displayed using
 ;$display_scale
 ngx = dimen($gx,0)	ngy = dimen($gy,1)
 ;now compute gx and gy to be just 1-D arrays for the x and y grid points
 gx = $display_scale*long($gx(*,0))
 gy = $display_scale*long($gy(0,*))
 pw = 5
 pixelmaskarea, pw, window, in, nx, ny
 ty,'nx, ny = ', nx,ny
 ;now find the gx and gy's contained in that area
 ix = in%nx
 iy = in/nx
 igxs = sieve(gx ge min(ix) and gx le max(ix))
 igys = sieve(gy ge min(iy) and gy le max(iy))
 ;that narrows the cases to the rectangle
 ngxs = num_elem(igxs)
 ngys = num_elem(igys)
 gin = zero(word($gx)
 for j=0,ngys-1 do {
  jg = igys(j)
  igy = fix(gy(jg) )
  for i=0,ngxs-1 do {
  ig = igxs(i)
  igx = fix(gx(ig)) )
  gv = igx + nx * igy
  if min(abs(gv - in)) eq 0 then { gin(ig,jg) = 1 }
   }}
  gind = sieve(gin eq 1)
 endsubr
 ;===============================================================================
;define some parameters used
$worm_flag = 0   $movie_flag = 0  $display_scale = 1
;to use - 	1. run browser and load the flowmap
;		2. dataread, del, 'file'
;		3. @/net/shimmer/usr/people/shine/ana/tools/corkstuff2.ana
;		4. show_corks, del