2012-07-18

いびつな表面でのランダムな成長

MIKU

気相・液相・固相は物質の３相
気相から固相に移って、固体が大きくなるときの話
空間があって、時間がある
空間には状態が定められ、状態の境界が定まる
あるルールがあるときに境界の形状がどのように時間変化するかをモデル化する話
n次元空間の境界はn-1次元多様体
境界の「平等な」形は「(超)平面」？
平等である場合「(超)平面」には限らない？
曲面であって、平等ならば、定曲率の曲面(球面)
3次元球の曲率〜n次元球の曲率→曲率って？
平面は曲率一定でかつ、「対称」な面
曲率一定であるときには時間変化は平等
曲率一定でないけれども、曲率が滑らかに変化するときには、変化も滑らか、ただし、一定ではない
曲率の変化が離散的なときには…
曲率の変化(一次微分)ではなくて、より高次な微分の一定でない・離散的なときというのは、影響がある？？？
シミュレーションするとしたら
- 境界を取り扱う
- ランダムな過程を確率分布として扱う
- Rの扱いをマスターして、上記、実施してみる
２次元空間で底に水分子がぶつかって底面が成長すると仮定する
底面は折れ線とする

Cross2Segments <- function(s1, s2){
	a1 <- (s1[2,2]-s1[1,2])/(s1[2,1]-s1[1,1])
	b1 <- s1[1,2] - s1[1,1]*a1
	a2 <- (s2[2,2]-s2[1,2])/(s2[2,1]-s2[1,1])
	b2 <- s2[1,2] - s2[1,1]*a2
	
	x0 <- (b2-b1)/(a1-a2)
	inside <- FALSE
	y0 <- a2 * x0 + b2
	dist <- sqrt((x0-s2[1,1])^2 + (y0-s2[1,2])^2) 
	if(((x0-s1[1,1])*(x0-s1[2,1]) < 0) & ((x0-s2[1,1])*(x0-s2[2,1]) < 0)){
		inside <- TRUE
	}
	return(list(inside = inside, dist =dist))
}
s1 <- matrix(c(0,0,1,0),byrow=TRUE,2,2)
s2 <- matrix(c(-1,1,0,0.3),byrow=TRUE,2,2)

Cross2Segments(s1,s2)

n.iter <- 10

n.pt <- 100
x <- seq(from = 0 ,to = 1, length=n.pt)
ys <- matrix(0, n.iter,n.pt)
n <- 1000
L <- 1
D <- 0.1
for(i in 2:n.iter){
	r <- matrix(c(runif(n),runif(n,min = min(ys[i-1,]),max=max(ys[i-1,])+L)),ncol=2)
	theta <- runif(n)*2*pi
	hit<-rep(0,n.pt-1)
	for(j in 1:n){
		small <- max(which(x < r[j,1]))
		if(small <= n.pt){
			big <- small + 1
			x.small <- x[small]
			x.big <- x[big]
			height.small <- ys[i-1,small]
			height.big <- ys[i-1,big]
			
			height <- height.small + (height.big-height.small) * (r[j,1]-x.small)/(x.big-x.small)
			if(height < r[j,2]){
				# above the surface
				tmphit <- 0
				tmpdist <- Inf
				for(k in 1:(n.pt-1)){
					s1 <- matrix(c(x[k],ys[i-1,k],x[k+1],ys[i-1,k+1]),byrow=TRUE,2,2)
					s2 <- matrix(c(r[j,],r[j,]+L*c(cos(theta[j]),sin(theta[j]))),byrow=TRUE,2,2)
					
					tmp <- Cross2Segments(s1,s2)
					if(tmp$inside){
						if(tmp$dist < L){
							if(tmp$dist < tmpdist){
								tmpdist <- tmp$dist
								tmphit<-k
							}
							#hit[j] <- hit[j] + 1
						}
					}
				}
				hit[tmphit] <- hit[tmphit]+1
			}
		}
	}
	tmpys <- rep(0,length(ys[i-1,]))
	for(j in 1:length(hit)){
		tmpys[j] <- tmpys[j] + hit[j]/2 * D
		tmpys[j+1] <- tmpys[j+1] + hit[j]/2 *D
	}
	ys[i,] <- ys[i-1,] + tmpys
}

matplot(t(ys),type="l")