Movement within a circle
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:
On Dec 15, 2008, at 5:20 PM, Juliane Struve wrote:
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 .
That is a bit vague for an audience that is numerically oriented.
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)
You have initialized this object with 10 rows, but why would you ? specify the distance to the walls as 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]
See below. That looks all wrong. should just be [i] not [day=i]
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
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ^^^^ ? ? ? ? ? ? ? ? ? ? ? ? ? ?here you need to learn to use <- for ? assignment
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])
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
#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
You might want to look at the worked example here: http://bm2.genes.nig.ac.jp/RGM2/R_current/library/animation/man/brown...
?Dr. Juliane Struve Environmental Scientist 10, Lynwood Crescent Sunningdale SL5 0BL 01344 620811
______________________________________________ R-h... at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
______________________________________________ R-h... at r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.