Skip to content

Movement within a circle

3 messages · Juliane Struve, David Winsemius, Andrew

Juliane Struve
# 
Dear list,
?
I am trying to program semi-random movement within a circle, with no particles leaving the circle. I would like them to bounce back when they come to close to the wall, but?I don't seem to be able to get this right. ?Would somebody?be able to?give me a hint ? This is my code so far, the particle starts at some point and moves towards the wall, but I don't get the "bouncing off" part right?.
?
Any help would be much appreciated.
?
Juliane 
?
days=10
circularspace=data.frame(day=c(0:days),xcoord=1,ycoord=1,xvelocity=1,yvelocity=1,xdistwall=0, ydistwall=0, wallxvel=0, wallyvel=0,stochasticxvel=0,stochasticyvel=0)
xmax=10
xmin=-10
ymax=10
ymin=-10
mindist=8
plot(xmin:xmax, ymin:ymax, type = "n")
circularspace
radius=10
timesteplength=1
weightfactor=1
for(i in 1:days)
{
#This is the stochastic component of the movement
circularspace$stochasticxvel[day=i+1]=runif(1,min=-1,max=1)
circularspace$stochasticyvel[day=i+1]=runif(1,min=-1,max=1)
#This is calculating the movment speed 
circularspace$xvelocity[day=i+1]=weightfactor*(circularspace$xvelocity[day=i])+ (1-weightfactor)*circularspace$stochasticxvel[day=i+1]
circularspace$yvelocity[day=i+1]=weightfactor*(circularspace$yvelocity[day=i])+ (1-weightfactor)*circularspace$stochasticyvel[day=i+1]
#This is calculating the new position
circularspace$xcoord[day=i+1]=circularspace$xcoord[day=i]+circularspace$xvelocity[day=i]/timesteplength
circularspace$ycoord[day=i+1]=circularspace$ycoord[day=i]+circularspace$yvelocity[day=i]/timesteplength
#Tis is checking the distance from the wall
circularspace$xdistwall[day=i+1]=xmax-abs(circularspace$xcoord[day=i])
circularspace$ydistwall[day=i+1]=ymax-abs(circularspace$ycoord[day=i])
#This is updating the movement if distance to wall too small
if (circularspace$xdistwall[day=i+1] <= mindist){
circularspace$wallxvel[day=i+1]= -1*(circularspace$xcoord[day=i+1])
circularspace$xvelocity[day=i+1]=weightfactor*circularspace$xvelocity[day=i]+ (1-weightfactor)*circularspace$stochasticxvel[day=i+1]+ circularspace$wallxvel[day=i+1]
circularspace$xcoord[day=i+1]=circularspace$xcoord[day=i]+circularspace$xvelocity[day=i]/timesteplength
}
if (circularspace$ydistwall[day=i+1] <= mindist){
circularspace$wallyvel[day=i+1]= -1*(circularspace$ycoord[day=i+1])
circularspace$yvelocity[day=i+1]=weightfactor*circularspace$yvelocity[day=i]+ (1-weightfactor)*circularspace$stochasticyvel[day=i+1]+ circularspace$wallyvel[day=i+1]
circularspace$ycoord[day=i+1]=circularspace$ycoord[day=i]+circularspace$yvelocity[day=i]/timesteplength
}
points(circularspace$xcoord[day=i+1],circularspace$ycoord[day=i+1])
symbols(0,0,circles=radius,inches=F,add=T)
symbols(0,0,circles=radius-2,inches=F,add=T)
}
circularspace
?Dr. Juliane Struve
Environmental Scientist
10, Lynwood Crescent
Sunningdale SL5 0BL
01344 620811
David Winsemius
#
On Dec 15, 2008, at 5:20 PM, Juliane Struve wrote:

            

      
      
      
    

        
That is a bit vague for an audience that is numerically oriented.

            

      
      
      
    

        
You have initialized this object with 10 rows, but why would you  
specify the distance to the walls as 0?

            

      
      
      
    

        
See below. That looks all wrong. should just be [i] not [day=i]

            

      
      
      
    

        
^^^^
                           here you need to learn to use <- for  
assignment

            

      
      
      
    

        
I would have expected to see code that checked whether the radial  
distance were greater than 10,
To do so would require calculating x^ + y^2. At the moment, you appear  
to have a square bounding box.

And why is the distance on day 5 calculated in terms of day 4?

And you need to look at the definitions of logical operators. "=" is  
not a logical operator. Above you were making unnecessary assignments,  
now you are conflating "=" , one of the assignment operators, with  
"==", the logical test for equality.

?"==" # or
?Comparison

            

      
      
      
    

        
You might want to look at the worked example here:

http://bm2.genes.nig.ac.jp/RGM2/R_current/library/animation/man/brownian.motion.html

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Andrew
    

    

      
      #
    

  


  

      
The following has a bug in it, but it does a reasonable job.  The bug
occurs in the lines that says

if(inner(xnew,xnew) > R^2)
   xnew <- x[n,] - eps ##error - sometimes x[n,] - eps is outside
circle

If it is changed to something like xnew <- proj(x[n,] + eps) where
proj is a projection function onto the circle, you should be ok.

HTH,

Andrew.

##=========================================================
N <- 1000 #steps
x <- matrix(0, nrow=N+1, ncol =2)
R <- 1 #radius of circle
delta <- 0.5 #step size

inner <- function(x,y) {
	if(length(x) != length(y)){
		print("Wrong length")
		return(0)
	}

	return(t(x) %*% (y))
}

for(n in 1:N){
	eps <- delta*runif(2)
	xnew <- x[n,] + eps
	if(inner(xnew,xnew) > R^2)
		xnew <- x[n,] - eps ##error - sometimes x[n,] - eps is outside
circle

	x[n+1, ]<- xnew

}

semicircle <- function(x, centre = c(0,0), radius=R){
	return(sqrt(radius^2 - x^2))
}
xvaries <- seq(-R,R,by=0.01)
plot(xvaries, semicircle(xvaries), type = 'l', xlim= c(-R, R), ylim =c
(-R,R))
lines(xvaries,- semicircle(xvaries))
points(x)

        
On Dec 16, 10:34?am, David Winsemius <dwinsem... at comcast.net> wrote: