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/brownian.motion.html
Dr. Juliane Struve
Environmental Scientist
10, Lynwood Crescent
Sunningdale SL5 0BL
01344 620811
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
