トーラス上を歩かせる

• 歩き方が少しおおざっぱだけれど、トーラス上を歩く様子を２次元に表示させる
• そのためにトーラスドーナツを上下に分けて開いて片方を上面、もう片方を下面として描くこととする

cs<-matrix(c(2,0,-2,0),ncol=2,byrow=TRUE)

r1<-1
r2<-2



Niter<-100000
xy<-matrix(0,Niter,2)
rt<-matrix(0,Niter,2)
donut<-rep(0,Niter)

initr<-runif(1)+1
initt<-runif(1)*2*pi
xy[1,]<-c(initr*cos(initt),initr*sin(initt))+cs[1,]
rt[1,]<-c(initr,initt)
donut[1]<-1
dt<-0.1
direction<-runif(1)*2*pi
library(MCMCpack)
prev<-rdirichlet(1,rep(1,2))
prev<-prev/sqrt(sum(prev^2))*dt
lon<-rep(0,Niter)
lon[1]<-sum(prev^2)
for(i in 2:Niter){
	current<-prev+rdirichlet(1,rep(1,2))*dt*0.1
	#current<-prev
	current<-current/sqrt(sum(current^2))*dt
	tmpxy<-xy[i-1,]+current
	lon[i]<-sum(current^2)
	prev<-current
	tmpr<-sqrt(sum((tmpxy-cs[donut[i-1],])^2))
	if(tmpr>r2){
		difr<-tmpr-r2
		tmpxy[1]<--tmpxy[1]
		donut[i]<-abs(donut[i-1]-2)+1
		tmpxy2<-tmpxy-cs[donut[i],]
		tmpxy2<-(tmpr-2*difr)/tmpr * tmpxy2
		xy[i,]<-tmpxy2+cs[donut[i],]
		prev<-c(prev[1],-prev[2])
	}else if(tmpr<r1){
		difr<-r1-tmpr
		tmpxy[1]<--tmpxy[1]
		donut[i]<-abs(donut[i-1]-2)+1
		tmpxy2<-tmpxy-cs[donut[i],]
		tmpxy2<-(tmpr+2*difr)/tmpr * tmpxy2
		xy[i,]<-tmpxy2+cs[donut[i],]
		prev<-c(prev[1],-prev[2])
	}else{
		xy[i,]<-tmpxy
		donut[i]<-donut[i-1]
	}
}

plot(xy,cex=0.5,col=rainbow(1000),xlim=c(-5,5),ylim=c(-5,5),pch=16)
tt<-seq(from=0,to=1,length.out=100)*2*pi
par(new=TRUE)
en1x<-r1*cos(tt)+cs[1,1]
en1y<-r1*sin(tt)+cs[1,2]
en2x<-r2*cos(tt)+cs[1,1]
en2y<-r2*sin(tt)+cs[1,2]
en3x<-r1*cos(tt)+cs[2,1]
en3y<-r1*sin(tt)+cs[2,2]
en4x<-r2*cos(tt)+cs[2,1]
en4y<-r2*sin(tt)+cs[2,2]
enX<-c(en1x,en2x,en3x,en4x)
enY<-c(en1y,en2y,en3y,en4y)
plot(enX,enY,cex=0.1,xlim=c(-5,5),ylim=c(-5,5))