そんなことより豆まきしようぜ!!

循環器内科・心臓血管外科学を学んでいる。
8割方動脈硬化の話だ。MIだ。
そういえば生理学で、血管内の血流はずり応力による層流になるとかならんとか言ってた気がする…。
NKの拡散の話を聞いて、まあなんか簡単に書いてみた。
　
x,y,z空間で、最初はz=0（時間t=0）から始まる。
xy平面において、血管腔は半径R。
血管は剛体。形は変わらない。
血液を限りなく微小にし、質量1の質点的なものとする。
あるzでのxy平面上の血液は、自分の存在するマスの周囲8マスから、隣接する血液の速度に応じた力を受ける。
つまり、自分の隣の血液が1の速度で動いていたら、fric倍の力を受ける。これを8マス分合算。
速度と書いたけど、xy平面に垂直に移動しているものとする。
高校物理を思い出すと

とかあったので、加速度から時間tにおける速度、変位を出す。
ニョキッとしたやつが描きたい。

# パラメータ
n<- 31                 # 考える世界の広さ。
R<- 5                  # 血管半径。
fric<- 0.0001          # 摩擦係数的な。

A<- B<-  diag(0,n)
Astart<- diag(1,n)
Layout<- matrix(c(2,1,1,1,3,
                  1,2,1,3,1,
                  1,1,1,1,1,
                  1,1,1,1,1,
                  1,1,1,1,1),5,5,byrow=T)
for(i in 1:n){
  for(j in 1:n){
    if(((i-ceiling(n/2))^2+(j-ceiling(n/2))^2)<=R^2){
      B[i,j]<- 1       # 0でないところが血管腔。
    }
  }
}
V<- diag(0,n+2)           # 端の影響のため、付け足す。
accel<- diag(0,n)         # 加速度。
pos<- diag(0,n)           # z軸方向への変位。
V[2:(n+1),2:(n+1)]<-  B
C<- B
Z<- 100
for(t in 1:Z){
  V[2:(n+1),2:(n+1)]<- C
  for(a in 2:(n+1)){
    for(b in 2:(n+1)){
      accel[a-1,b-1]<- fric*(V[a-1,b+1]+  V[a,b+1]+V[a+1,b+1]+
                             V[a-1,b  ]+0*V[a,b  ]+V[a+1,b  ]+
                             V[a-1,b-1]+  V[a,b-1]+V[a+1,b-1])
    }
  }
  C<- C+accel*t           # Cは速度だな。
  pos<- pos+C*t+(accel*t^2)/2
  pos[B==0]<- 0
  #file.name=paste("test",t,".png",sep="")
  #png(file=file.name)
    layout(Layout)
    persp(pos,col=2,phi=20,theta=60,zlim=c(0,80000),
          xlab="x",main="position")
    image(accel,main="accel",col=topo.colors(1000))  
    image(C,main="velocity",col=topo.colors(1000))
  #dev.off()
}

左上がz軸上方から眺めた加速度、右上が同じ視点の速度、perspがシミュレーション上の血流。
発展させる予定は今のところなし。
NKがセミナーで言っていたように、xyを一辺に計算できたらもっと計算速度があがるのになぁ。
グラフィックもぎこちないぞ。層っていうくらいだから同じ速度のところを色分けするのもいいと思う。
スクリプト書いてからwikiったが、もっと真面目にモデル化しておいたほうがきれいな図が描けそう。
　
そんなことより豆まきしようぜ!!