3dで拡散

• こちらとこちらの続き
• 線維方向を考えたうえで3Dに拡張する
  • 6近傍への移動を考える
• 3Dで線維方向を含めた誘電率の値はこちらに
• 頂点への移動が少ない可能性がある

library(rgl)

Nx<-30
Ny<-30
Nz<-30
Nt<-100
U<-array(0,c(Nx,Ny,Nz))
U[11:20,11:20,11:20]<-1   #初期分布

kl<-0.15
kt<-0.02

theta<-pi/4
theta<-array(theta,c(Nx,Ny,Nz))
phi<-pi/4
phi<-array(phi,c(Nx,Ny,Nz))

kxx<- kl*cos(theta)^2*cos(phi)^2 + kt*(1-cos(theta)^2*cos(phi)^2)
kxy<- (kl - kt)*cos(theta)*sin(theta)*cos(phi)^2
kxz<-(kl - kt) *cos(theta) *cos(phi) *sin(phi)

kyx<-(kl - kt)*cos(theta) *sin(theta) *cos(phi)^2
kyy<-kl *sin(theta)^2 *cos(phi)^2 +kt*(1-sin(theta)^2 *cos(phi)^2)
kyz<-(kl -kt) *sin(theta) *cos(phi) *sin(phi)

kzx<-(kl - kt) *cos(theta) *cos(phi) *sin(phi)
kzy<-(kl -kt) *sin(theta) *cos(phi) *sin(phi)
kzz<-kl *sin(phi)^2 + kt *cos(phi)^2

kxx<-(kxx[-1,,]+kxx[-Nx,,])/2
kxy<-(kxy[,-1,]+kxy[,-Ny,])/2
kxz<-(kxz[,,-1]+kxz[,,-Nz])/2

kyx<-(kyx[-1,,]+kyx[-Nx,,])/2
kyy<-(kyy[,-1,]+kyy[,-Ny,])/2
kyz<-(kyz[,,-1]+kyz[,,-Nz])/2

kzx<-(kzx[-1,,]+kzx[-Nx,,])/2
kzy<-(kzy[,-1,]+kzy[,-Ny,])/2
kzz<-(kzz[,,-1]+kzz[,,-Nz])/2

for (t in 1:Nt){
	dUx<-U[-1,,]-U[-Nx,,]
	dUy<-U[,-1,]-U[,-Ny,]
	dUz<-U[,,-1]-U[,,-Nz]
	
	Ixy<- -dUy*kxy
	Iyx<- -dUx*kyx

	Iyz<- -dUz*kyz
	Izy<- -dUy*kzy

	Izx<- -dUx*kzx
	Ixz<- -dUz*kxz

	Ixypp<-(kxy>0)*(-dUy>0)*Ixy
	Ixypm<-(kxy<0)*(-dUy>0)*Ixy
	Ixymp<-(kxy<0)*(-dUy<0)*Ixy
	Ixymm<-(kxy>0)*(-dUy<0)*Ixy
	
	Iyxpp<-(kyx>0)*(-dUx>0)*Iyx
	Iyxpm<-(kyx<0)*(-dUx>0)*Iyx
	Iyxmp<-(kyx<0)*(-dUx<0)*Iyx
	Iyxmm<-(kyx>0)*(-dUx<0)*Iyx
	
	Iyzpp<-(kyz>0)*(-dUz>0)*Iyz
	Iyzpm<-(kyz<0)*(-dUz>0)*Iyz
	Iyzmp<-(kyz<0)*(-dUz<0)*Iyz
	Iyzmm<-(kyz>0)*(-dUz<0)*Iyz
	
	Izypp<-(kzy>0)*(-dUy>0)*Izy
	Izypm<-(kzy<0)*(-dUy>0)*Izy
	Izymp<-(kzy<0)*(-dUy<0)*Izy
	Izymm<-(kzy>0)*(-dUy<0)*Izy

	Izxpp<-(kzx>0)*(-dUx>0)*Izx
	Izxpm<-(kzx<0)*(-dUx>0)*Izx
	Izxmp<-(kzx<0)*(-dUx<0)*Izx
	Izxmm<-(kzx>0)*(-dUx<0)*Izx
	
	Ixzpp<-(kxz>0)*(-dUz>0)*Ixz
	Ixzpm<-(kxz<0)*(-dUz>0)*Ixz
	Ixzmp<-(kxz<0)*(-dUz<0)*Ixz
	Ixzmm<-(kxz>0)*(-dUz<0)*Ixz
	
	Ix<--dUx*kxx
	Ix[,-Ny,]<-Ix[,-Ny,]+Ixypp[-Nx,,]+Ixypm[-1,,]
	Ix[,-1,]<-Ix[,-1,]+Ixymp[-Nx,,]+Ixymm[-1,,]	
	Ix[,,-Nz]<-Ix[,,-Nz]+Ixzpp[-Nx,,]+Ixzpm[-1,,]
	Ix[,,-1]<-Ix[,,-1]+Ixzmp[-Nx,,]+Ixzmm[-1,,]	
	
	Iy<--dUy*kyy
	Iy[-Nx,,]<-Iy[-Nx,,]+Iyxpp[,-Ny,]+Iyxpm[,-1,]
	Iy[-1,,]<-Iy[-1,,]+Iyxmp[,-Ny,]+Iyxmm[,-1,]
	Iy[,,-Nz]<-Iy[,,-Nz]+Iyzpp[,-Ny,]+Iyzpm[,-1,]
	Iy[,,-1]<-Iy[,,-1]+Iyzmp[,-Ny,]+Iyzmm[,-1,]
	
	Iz<--dUz*kzz
	Iz[,-Ny,]<-Iz[,-Ny,]+Izypp[,,-Nz]+Izypm[,,-1]
	Iz[,-1,]<-Iz[,-1,]+Izymp[,,-Nz]+Izymm[,,-1]
	Iz[-Nx,,]<-Iz[-Nx,,]+Izxpp[,,-Nz]+Izxpm[,,-1]
	Iz[-1,,]<-Iz[-1,,]+Izxmp[,,-Nz]+Izxmm[,,-1]
	
	U[-1,,]<-U[-1,,]+Ix
	U[-Nx,,]<-U[-Nx,,]-Ix
	U[,-1,]<-U[,-1,]+Iy
	U[,-Ny,]<-U[,-Ny,]-Iy
	U[,,-1]<-U[,,-1]+Iz
	U[,,-Nz]<-U[,,-Nz]-Iz

	plot3d(slice.index(U,1),slice.index(U,2),slice.index(U,3),col=rainbow(10)[ifelse(U < 0.1,1,U*10)],alpha=ifelse(U < 0.1,0,1))
#	plot3d(slice.index(U,1),slice.index(U,2),slice.index(U,3),alpha=ifelse(U == 0,0,1))
#	image(U[,,2],col=topo.colors(100))
}