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

• 昨日の続き
• 初期状態が波打った地面だと…
• 山が増幅される
• 山の頂点がへこみがちなのは、折れ線近似のせい？？

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)
ys[1,] <- sin(x*pi*10)
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")