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...@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...@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...@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.
______________________________________________
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.
