2013-12-25

らせん

R らせん

昨日はDNA二重らせんを描いてみた
３次元空間におけるらせんについて少し考えてみる
結果としてこんな絵ができる

DNAは２重らせんを作った上で、さらに高次構造を作っている。ヒストンとのタンパク複合体を作ったり、とその現実界的折りたたみ構造についてわかっていることは多いが、さて一般に、ひも状のものをくるくると畳むにはどうすればよいのか、という点に照らして考えてみる
そもそも、(２重)らせんの「らせん」が「まっすぐに伸ばせは長いひも」をくるくる巻きにすることで、「線の太さ」はくるくる直径分に長くなるけれど、長く伸びたひもっぽい長さを短くする、という仕組みになっている
これを一般化しよう
まず普通のらせんを「１次のらせん」と呼ぶことにして、「高次化」への準備をする
「１次のらせん」を円だけで表現する
- 昨日の記事でも書いたように、１次のらせんは、xy平面に円を描きながら、z方向に等速で移動したときに描かれる軌跡である
- ここには、「円」と「直線」という２種類の動きが出てくる
- 「高次化したらせん」では、らせん軌道を進みながら、その軌道上の点を中心として円を描かせたい
- そのようなとき「直線」を進みながら「周回円」を描くという「１次のらせん」ではなく、「円的な部分を持った曲線」を進みながら「周回円」を描きたい
- したがって「１次のらせん」で登場する「直線」を「円」的なものとして表すことにしよう
- 直線は半径が無限大の円弧であると考えれば、直線は円として表現できる
- このようにすることで、登場するすべての動きを「円」としてしまう算段がついた
「１次のらせん」とトーラスの関係
- 「１次のらせん」を円弧に沿って進みながら、その円弧の垂直断面に円を描くと、それは、円の組合せであって、いわゆるトーラス上の軌道になる
「いわゆるトーラス」上の軌道
- トーラス上の軌道は
  - $x=R1\cos(\theta1 t) + R2\cos(\theta2 t) * \cos(\theta1 t)$
  - $y=R1\sin(\theta1 t) + R2\cos(\theta2 t) * \sin(\theta1 t)$
  - $z=R2\sin(\theta2 t)$
- と表すことができる。ここで $R1,R2$ は第１、第２の円の半径であり、 $\theta1,\theta2$ は第１、第２の円の周回の速度を決めるパラメタである
「いわゆるトーラス」上の軌道を表現の見通しを良くする
- 第１の円
  - 第１の円は、原点を中心に、ある１点をスタート点とした円である
  - そのような円の上の点を $\mathbf{x}_1$ とし、原点を $\mathbf{x}_0$ とする
  - このとき、この点は、ベクトル $\mathbf{x}_1-\mathbf{x}_0$ に垂直な進行方向ベクトル $\mathbf{v}_1$ がある
  - ここで、２つのベクトル $\mathbf{x}_1-\mathbf{x}_0$ と $\mathbf{v}_1$ とに垂直なベクトルを、この２ベクトルの外積ベクトルとしてとりそれを $\mathbf{z}_1$ と置くことにしよう
- 第２の円
  - 第２の円は $\mathbf{x}_1$ を中心とし、 $\mathbf{x}_1-\mathbf{x}_0$ と $\mathbf{z}_1$ とが張る平面上の円である
  - このような第２の円について、半径 $R2$ 角 $\theta2$ だけ回った点は、xy平面上でのそのような点を、 $(1,0,0),(0,1,0),(0,0,1)$ という座標軸を $(\mathbf{x}_1-\mathbf{x}_0)_{st},(\mathbf{z}_1)_{st},(\mathbf{v}_)_{st}$ という座標軸(ただし $v_{st}$ はベクトル $v$ を単位ベクトル化したもの)に回転したものと考えればよい
- 結局、上記の処理は、ある回転運動aをしている点を中心に回転bしているときの点を、存在位置座標と、そこでの回転aの法線ベクトルと回転方向ベクトルとから、その外積ベクトルを介して回転bによる座標の変化分を加える、という読み方にしたものである
- これは、第１の円、その上の第２の円、さらにその上の第３の円…と積み重ねることが容易な表現になっている
以下はこれをRで実装したもの

周期が $2\pi$ の整数倍だと、軌道は「元に戻る」ので、曲線が次のように閉じる

library(rgl)
# まずは全体の原点を決める
ori <- c(0,0,0)
# 第１、第２、...の円の半径を与える。だんだん短くするのが見やすい絵ができる
# 逆に周期の係数は第１、第２、…の順にだんだん早くするのがよい
Rs <- 1.3^(10:1)
Ts <- 2^(1:10)
# 円の数
n.steps <- length(Rs)
# 軌道の長さを決めるパラメタTs * max.tの角度まで回る
# この値が2pi　で一周
max.t <- 4
# 定めた角に等間隔で打点する、その点の数
t.n <- 10000
t <- seq(from=0,to=max.t,length=t.n)
# 行列を与えて、各行のベクトルを単位ベクトル化する関数
unit.vector <- function(X){
	v <- sqrt(apply(X^2,1,sum))
	X/v
}
# ３次元ベクトル間の外積ベクトル
outer.prod <- function(v1,v2){
	c(v1[2]*v2[3]-v1[3]*v2[2],v1[3]*v2[1]-v1[1]*v2[3],v1[1]*v2[2]-v1[2]*v2[3])
}
# ２つの行列を引数に取り、対応する行をベクトルとしてその外積ベクトルを返す
outer.prod.V <- function(V1,V2){
	cbind(V1[,2]*V2[,3]-V1[,3]*V2[,2],V1[,3]*V2[,1]-V1[,1]*V2[,3],V1[,1]*V2[,2]-V1[,2]*V2[,1])
}
# ここから、座標計算
# Xは座標ベクトルを第０，１，２、…の円で積み重ねて格納、Vはそれぞれ第i円での速度ベクトルを格納
X <- list()
V <- list()
# X[[1]],V[[1]]は第０円(これは原点で、静止)
X[[1]] <- V[[1]] <- cbind(rep(0,t.n),rep(0,t.n),rep(0,t.n))
# X[[2]],V[[2]]は第１円の座標と速度ベクトル。xy平面上に取る
X[[2]] <- cbind(Rs[1]*cos(Ts[1]*t),Rs[1]*sin(Ts[1]*t),rep(0,t.n))
V[[2]] <- cbind(-Rs[1]*sin(Ts[1]*t),Rs[1]*cos(Ts[1]*t),rep(0,t.n))
# 第２円からそれ以降を順次計算する
for(i in 2:n.steps){
# 計算するのは第i円で、それはX[[i+1],V[[i+1]]に格納する
# 一つ前の段階の第i-1円に関する、中心と座標との差ベクトル
	tmp.x <- X[[i]]-X[[i-1]]
# その単位ベクトル化
	tmp.x <- unit.vector(tmp.x)
# 第i円の中心の速度ベクトルの単位ベクトル化
	tmp.v <- unit.vector(V[[i]])
# 外積ベクトルを取り出し
	tmp.z <- outer.prod.V(tmp.x,tmp.v)
# xy平面上の原点を中心とした円だったらどのような座標に相当するか、そこでの速度ベクトルはいくつかを計算し
	tmp.X <- cbind(Rs[i]*cos(Ts[i]*t),Rs[i]*sin(Ts[i]*t),rep(0,t.n))
	V[[i+1]] <- cbind(-Rs[i]*sin(Ts[i]*t),Rs[i]*cos(Ts[i]*t),rep(0,t.n))
# 座標系に合わせて回転する
	for(j in 1:t.n){
		tmp.X[j,] <- cbind(tmp.x[j,],tmp.z[j,],tmp.v[j,]) %*% tmp.X[j,]
		V[[i+1]][j,] <- cbind(tmp.x[j,],tmp.z[j,],tmp.v[j,]) %*% V[[i+1]][j,]
	}
# 第i円まで積み重ねた座標は、一段階前の円までを積み重ねた座標に加算する
	X[[i+1]] <- X[[i]] + tmp.X
}
# ３次元プロットするにあたり、座標軸を見せず、また３軸のスケールがそろうように、表示枠を作る
plot3d.minmax <- function(X){
	min.X <- min(X)
	max.X <- max(X)
	tmp <- rbind(X,rep(min.X,length(X[1,])),rep(max.X,length(X[1,])))
	plot3d(tmp,col="white",axes=FALSE)
}
# 各点の球の大きさは、全体の枠と最近点間距離とで適当に選ぶので、そのためにいくつかの値を計算する
tmp.d <- dist(X[[n.steps+1]])
max.X <- max(X[[n.steps+1]])
min.X <- min(X[[n.steps+1]])
dd <- max.X-min.X
# 全部の円を積み重ねた座標をプロットしよう
plot3d.minmax(X[[n.steps+1]])
# 各点に球を打点する
spheres3d(X[[n.steps+1]],col=2,radius=max(min(tmp.d)*0.5,dd*0.001))
# 隣接点間に線分を引く
s <- rep(1:t.n,each=2)
s <- s[2:(length(s)-1)]
segments3d(X[[n.steps+1]][s,])