方向を自由に与える

拡散
  • こちらのつづき
  • kの値を各位置について与える
    • まず \theta を各位置で与え、その \theta から k を計算する
Nx<-30
Ny<-30
Nt<-100
U<-tempU<-matrix(0,Nx,Ny)
U[11:20,11:20]<-1   #初期分布

#これがl方向とt方向の拡散の定数
kl<-0.20
kt<-0.05

#線維走向の傾き
theta<-matrix(0,Nx,Ny)
for(x in 1:Nx){
for(y in 1:Ny){
phi<-pi/2
theta[x,y] <- atan2(y-(Ny+1)/2,x-(Nx+1)/2)+phi
#theta<-pi/4
}}

kxx<-kl*cos(theta)^2+kt*sin(theta)^2
kyy<-kl*sin(theta)^2+kt*cos(theta)^2
kxy<-(kl-kt)*sin(theta)*cos(theta)
kyx<-(kl-kt)*sin(theta)*cos(theta)

kxx<-(kxx[-1,]+kxx[-Nx,])/2
kxy<-(kxy[,-1]+kxy[,-Ny])/2
kyx<-(kyx[-1,]+kyx[-Nx,])/2
kyy<-(kyy[,-1]+kyy[,-Ny])/2

for (t in 1:Nt){
	dUx<-U[-1,]-U[-Nx,]
	dUy<-U[,-1]-U[,-Ny]

	Ixy<- -dUy*kxy
	Iyx<- -dUx*kyx
	
	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

	
	Ix<--dUx*kxx
	Ix[,-Ny]<-Ix[,-Ny]+Ixypp[-Nx,]+Ixypm[-1,]
	Ix[,-1]<-Ix[,-1]+Ixymp[-Nx,]+Ixymm[-1,]
	
	Iy<--dUy*kyy
	Iy[-Nx,]<-Iy[-Nx,]+Iyxpp[,-Ny]+Iyxpm[,-1]
	Iy[-1,]<-Iy[-1,]+Iyxmp[,-Ny]+Iyxmm[,-1]
		
	U[-1,]<-U[-1,]+Ix
	U[-Nx,]<-U[-Nx,]-Ix
	U[,-1]<-U[,-1]+Iy
	U[,-Ny]<-U[,-Ny]-Iy

	
#	persp(U,col=heat.colors(Nx*Ny),theta=30,phi=30,zlim=c(0,1))
	image(U,col=topo.colors(100))
#	print(sum(U))
}