xpos = 0
ypos = 0
xwalks = NULL
ywalks = NULL
deterministic.size = .25
arena.size = 20
angs = seq(0,2*pi+.001,length=200)
total.steps=0
walk.colors = NULL
direction = .5
show.direction = FALSE
show.average = FALSE
# =========================
startwalk = function(arena=20, nwalkers=100, detsize=0, show.direction=FALSE,average=FALSE, npaths=min(5,nwalkers)) {
  arena.size <<- arena
  deterministic.size <<- detsize
  show.average <<- average
  total.steps <<- 0
  xpos <<- rep(0,nwalkers-npaths) # altogether there will be nwalkers
  ypos <<- rep(0,nwalkers-npaths) # some shown as paths
  xwalks <<- list()
  ywalks <<- list()
  for (k in 1:npaths) {
    xwalks[[k]] <<- 0
    ywalks[[k]] <<- 0
  }
  direction <<- runif(1,0, 2*pi)
  show.direction <<- show.direction
  walk.colors <<- rainbow(npaths,alpha=.3)
  walk.colors[1] <<- rgb(0,0,0,1)
  showwalk(average=show.average)
}
# =========================
showwalk = function(average=FALSE,plot.size=NULL,dotsize=.5){
  if (average) {
    divisor = max(total.steps,1)
    arena.size <<- deterministic.size
    show.average <<- TRUE
  }
  else {
    show.average <<- FALSE
    divisor=1
  }
  if( show.average )
    label = paste( "Average Dist. per Step for n =", total.steps)
  else
    label = paste("Total Distance for n =", total.steps, "steps")
  if( is.null(plot.size) )
     lims = 2.5*c(-arena.size,arena.size)
  else
     lims = c(-plot.size, plot.size)
     
  plot( arena.size*cos(angs), arena.size*sin(angs), bty='n',
       xlim=lims, ylim=lims,
       type='l', lwd=3,col=rgb(1,0,0,.3), xaxt='n',yaxt='n',
       xlab="",ylab="",main=label)
  text( 0, arena.size, paste("Radius =",arena.size),
       pos=3, col=rgb(1,0,0,.3)) 
  points( xpos/divisor, ypos/divisor, pch=20, cex=dotsize)
  for (k in 1:length(xwalks) ) {
    lines( cumsum(xwalks[[k]])/divisor, cumsum(ywalks[[k]])/divisor, col=walk.colors[k])
    points( sum(xwalks[[k]])/divisor, sum(ywalks[[k]])/divisor, pch=20,cex=dotsize)
  }
  points(0,0,pch=20,col=rgb(1,0,0,.3), cex=4)
}
# =========================
walk = function(nsteps=1) {
  total.steps <<- total.steps + nsteps
  for (k in 1:length(xpos) ) {
    xpos[k] <<- xpos[k] + sum( rnorm(nsteps) + deterministic.size*cos(direction) )
    ypos[k] <<- ypos[k] + sum( rnorm(nsteps) + deterministic.size*sin(direction) )
  }
  for (k in 1:length(xwalks)) {
    xwalks[[k]] <<- c(xwalks[[k]] , rnorm(nsteps)+deterministic.size*cos(direction))
    ywalks[[k]] <<- c(ywalks[[k]] , rnorm(nsteps)+deterministic.size*sin(direction))
  }
  showwalk(average=show.average)
  if (show.direction) show.center()
}
# ==========================
show.center = function() {
  if (show.average) divisor = total.steps
  else divisor=1
  
  dist = total.steps*deterministic.size/divisor
  points( dist*cos(direction), dist*sin(direction), pch=20, cex=2,col=rgb(1,0,0,.8))
}
# ==========================
# show a 95% confidence interval around the first walker
show.ci = function() {
  ci.dist = 2.5*sqrt(total.steps)
  # Why 2.5 and not 1.96?
  # So that the probability of a random walker being in the circle
  # is 95%  
  if (show.average) divisor = total.steps
  else divisor = 1

  polygon( (sum(xwalks[[1]]) + ci.dist*cos(angs))/divisor,
         (sum(ywalks[[1]]) + ci.dist*sin(angs))/divisor,
        lwd=3,col=rgb(0,0,0,.2))
}
    
#
#===========================
# start the walkers in a quasi-random way
start.clump = function(){
library(gsl)
  q <- qrng_alloc(dim = 2)
  qrng_name(q)
  foo = qrng_get(q, length(xpos))
  xpos <<- 3*foo[,1] - 1.5
  ypos <<- 3*foo[,2] - 1.5
  foo = qrng_get(q, length(xwalks))
  for (k in 1:length(xwalks) ) {
    xwalks[[k]] <<- foo[k,1]
    ywalks[[k]] <<- foo[k,2]
  }

}