3dの拡散

• 修正したのはこちらへ
• 昨日の続き
  • 3dで拡散させて、値の大きさを色で表す
• こちらにZ軸を付け加える
  • ただし以下の書き方はナナメ方向には移動していない
• ドーナツはこの日のセミナーで

library(rgl)
Nx<-30
Ny<-30
Nz<-30

Nt<-30
R<-10
U<-array(0,c(Nx,Ny,Nz))
A<-array(0,c(Nx,Ny,Nz))
K1<-array(0.1,c(Nx-1,Ny,Nz))
K2<-array(0.1,c(Nx,Ny-1,Nz))
K3<-array(0.1,c(Nx,Ny,Nz-1))

#球を作る
for(x in 1:Nx){
for(y in 1:Nx){
for(z in 1:Nz){
	if((x-Nx/2)^2+(y-Ny/2)^2+(z-Nz/2)^2 <= R^2){A[x,y,z]<-1} 
}}}

#初期分布
U[1:ceiling(Nx/2),1:ceiling(Ny/2),1:ceiling(Nz/2)]<-1

U<-U*A
K1<-K1*A[-Nx,,]*A[-1,,]
K2<-K2*A[,-Ny,]*A[,-1,]
K3<-K3*A[,,-Nz]*A[,,-1]

#拡散
for (t in 1:Nt){
dUx<-U[2:Nx,,]-U[1:(Nx-1),,]
dUy<-U[,2:Ny,]-U[,1:(Ny-1),]
dUz<-U[,,2:Nz]-U[,,1:(Nz-1)]

U[1:(Nx-1),,]<-U[1:(Nx-1),,] +K1*dUx
U[2:Nx,,]<-U[2:Nx,,] -K1*dUx

U[,1:(Ny-1),]<-U[,1:(Ny-1),]+K2*dUy
U[,2:Ny,]<-U[,2:Ny,]-K2*dUy

U[,,1:(Nz-1)]<-U[,,1:(Nz-1)]+K3*dUz
U[,,2:Nz]<-U[,,2:Nz]-K3*dUz

plot3d(slice.index(U,1),slice.index(U,2),slice.index(U,3),col=rainbow(1000)[ifelse(U < 0.001, 1,U*1000)],alpha=ifelse(U < 0.001,0,1))
rgl.viewpoint(theta=300,phi=30)
}