• 作成例。Rを使い始めて11日目
n.t <- 6 # No. times
k <- 3  # 虫の種類数
n.init <- rep(1,k)  # 初期の虫の数
Loc.hx <- list()  # 位置の履歴格納
Stat.hx <- list()  # 虫の状態の格納
  S<-list()
  L<-list()
  S[[1]]<-list()
  L[[1]]<-list()
  #初期値の設定(t=1(初期)での種類ごとの値)
  for(i in 1:k){
  S[[1]][[i]]<-rep(1,5)
  L[[1]][[i]]<-matrix(0,5,2)
 }

##############################################

# 時間経過
for(t in 2:n.t){ 
  S[[t]]<-list()
  L[[t]]<-list()

##############################################
   #まず数を増やす
   for(i in 1:k){  
        
      #Sから1のものを選ぶ(S=1であれば必ず分裂するとする) (何番目に1があるか)
   f<-which(S[[t-1]][[i]]==1)
      #増やすべきものがfの個数あるので、SとLの要素数を増やす
   S[[t]][[i]]<-c(S[[t-1]][[i]],rep(1,length(f)))
   
      #増えた分の座標dL
      #親と子の座標は同じ
      dL<-L[[t-1]][[i]][f,]
   L[[t]][[i]]<-rbind(L[[t-1]][[i]],dL)
                   }

################################################
   #次に、酔歩させる
  for(i in 1:k){ 
       
   #歩いた分の座標
      N<-length(S[[t]][[i]])
      a<-rnorm(N)
      b<-rnorm(N)
   dzahyo<-cbind(a,b)
      #動いた後の座標
   L[[t]][[i]]<-L[[t]][[i]]+dzahyo
                   }
##################################################
#インターラクションする他種との距離を計算し、ステータスを決める
     #他種との組み合わせ
    for(i in 1:(k-1)){
    for(j in (i+1):k){
       
     #i種個体とj種個体の距離を計算して行列にする(#選ばれた二種との距離を要素とした行列をつくる)
                   
                   Dist<-matrix(0,length(S[[t]][[i]]),length(S[[t]][[j]]))
                   for(p in 1:length(S[[t]][[i]])){
                   for(q in 1:length(S[[t]][[j]])){
                                  #距離の計算関数
                                  a<-L[[t]][[i]][p,]
                                  b<-L[[t]][[j]][q,]
                                  d<-sqrt(sum((a-b)^2))
                                  #閾値lim以内ならT(1)と判定
                                  lim<-0.5
                                  if(d<lim){d<-1
                                  }else{d<-0
                                  }
                                  #個体間の閾値判定を行列に
                                  Dist[p,q]<-d
                   }}
                   s.ci<-which(apply(Dist,1,sum)>0)
                   s.cj<-which(apply(Dist,2,sum)>0)
                   #Tが一つでもあればSを1から0に変える
            S[[t]][[i]][s.ci]<-S[[t]][[i]][s.ci]-1
                   S[[t]][[j]][s.cj]<-S[[t]][[j]][s.cj]-1
                   }}

###############################################################################
#0の個体を無効に
                   for(i in 1:k){
          alive<-which(S[[t]][[i]]==1)
                   S[[t]][[i]]<-S[[t]][[i]][alive]
                   L[[t]][[i]]<-L[[t]][[i]][alive,]
                   }

    Stat.hx<- c(Stat.hx,  list(S)) # Stat.hxにステータスを格納
    Loc.hx <- c(Loc.hx,  list(L))
}

# オブジェクト L を簡単アニメーション化
oopt = ani.options(nmax = 20, ani.width = 600, ani.height = 500, 
interval = 0.2,outdir="./")
install.packages("animation")
library(animation)
ani.start()
#boot.iid()
for(i in 1:length(S)){
	for(j in 1:length(S[[i]])){
		if(j > 1){
			points(L[[i]][[j]],col=j,pch=pch.type[j])
		}else{
			plot(L[[i]][[j]],col=j,pch=20,xlim=xlim,ylim=
ylim)
			#par(ask = FALSE)
		}
		
	}
	#par(ask=TRUE)
}


ani.stop()
ani.options(oopt)