対数らせん、動標構、渦、接面

  • 対数らせんに沿って動標構を描く
  • 対数らせんをぐるりと回せば渦
  • 対数らせんの曲線上の動標構と渦とみたときの接面は少しずれたりすることを絵で描く

# 渦を作る
# 渦の軸を螺旋にする
t <- seq(from=-10,to=3,length=1000)*3
# 回転半径の拡大係数
kr <- 0.3
# 渦の進行方向の拡大係数
kz <- 0.2
# 回転の周波数を調整する係数
theta <- 3

x <- exp(kr*t) * cos(theta*t)
y <- exp(kr*t) * sin(theta*t)
z <- exp(kz*t)

xyz <- cbind(x,y,z)
xyz. <- xyz
xyz. <- rbind(xyz.,rep(min(xyz),3))
xyz. <- rbind(xyz.,rep(max(xyz),3))
library(rgl)
plot3d(xyz.)

phis <- seq(from=0,to=2*pi,length=51)
phis <- phis[-1]
# 1本の対数らせんをぐるりと回して渦にする
xyzall <- matrix(0,0,3)
for(i in 1:length(phis)){
	rot <- matrix(c(cos(phis[i]),sin(phis[i]),-sin(phis[i]),cos(phis[i])),2,2)
	tmp <- xyz
	tmp[,1:2] <- t(rot %*% t(tmp[,1:2]))
	xyzall <- rbind(xyzall,tmp)
}
xyzall. <- xyzall
xyzall. <- rbind(xyzall.,rep(min(xyzall),3))
xyzall. <- rbind(xyzall.,rep(max(xyzall),3))
library(rgl)
plot3d(xyzall.)

# 1本の対数らせんをもう一度作る
krA <- 0.35
kzA <- 0.25
thetaA <- 3.5
X <- exp(krA*t) * cos(thetaA*t)
Y <- exp(krA*t) * sin(thetaA*t)
Z <- exp(kzA*t)

XYZ <- cbind(X,Y,Z)
XYZ. <- XYZ./sqrt(apply(XYZ.^2,1,sum))

# 微分して速度ベクトル化
X. <- krA * X - theta * Y
Y. <- krA * Y + theta * X
Z. <- kzA * Z

XYZ. <- cbind(X.,Y.,Z.)

# さらに微分して加速度ベクトル化
tmpX.. <- krA * X. - theta * Y.
tmpY.. <- krA * Y. + theta * X.
tmpZ.. <- kzA * Z.

# 加速度ベクトルのうち、速度ベクトルと直交する成分を取り出す
XYZ.. <- cbind(tmpX..,tmpY..,tmpZ..)
range(apply(XYZ.*XYZ..,1,sum))

i.p <- apply(XYZ.*XYZ..,1,sum)
XYZ.. <- XYZ.. - i.p * XYZ.
range(apply(XYZ.*XYZ..,1,sum))
# 動標構
E1 <- E2 <- E3 <- XYZ
for(i in 1:length(XYZ[,1])){
	tmp.e2 <- XYZ.[i,]
	tmp.e3 <- XYZ..[i,]
	tmp.e1 <- -solve(rbind(tmp.e2[1:2],tmp.e3[1:2])) %*% c(tmp.e2[3],tmp.e3[3])
	tmp.e1 <- c(tmp.e1,1)
	E1[i,] <- tmp.e1/sqrt(sum(tmp.e1^2))
	E2[i,] <- tmp.e2/sqrt(sum(tmp.e2^2))
	E3[i,] <- tmp.e3/sqrt(sum(tmp.e3^2))
}
range(apply(E1*E2,1,sum))
range(apply(E2*E3,1,sum))
range(apply(E3*E1,1,sum))


E1 <- E2 <- E3 <- XYZ
for(i in 1:length(XYZ[,1])){
	tmp.e2 <- XYZ.[i,]
	tmp.e3 <- XYZ..[i,]
	tmp.e1 <- -solve(rbind(tmp.e2[1:2],tmp.e3[1:2])) %*% c(tmp.e2[3],tmp.e3[3])
	tmp.e1 <- c(tmp.e1,1)
	E1[i,] <- tmp.e1/sqrt(sum(tmp.e1^2))
	E2[i,] <- tmp.e2/sqrt(sum(tmp.e2^2))
	E3[i,] <- tmp.e3/sqrt(sum(tmp.e3^2))
}
range(apply(E1*E2,1,sum))
range(apply(E2*E3,1,sum))
range(apply(E3*E1,1,sum))
# 渦という曲面の法線
Zt <-  krA * exp(krA*t)
Xt <- - kzA * exp(kzA*t) * cos(thetaA*t)
Yt <- - kzA * exp(kzA*t) * sin(thetaA*t)
XYZt <- cbind(Xt,Yt,Zt)

# 3dplotの尺度の統一のために2点加える
XYZ.st <- XYZ
XYZ.st <- rbind(XYZ,rep(min(XYZ),3))
XYZ.st <- rbind(XYZ.st,rep(max(XYZ),3))

# 動標構築を描く
plot3d(XYZ.st)

ss <- 1:length(t)
for(i in 1:length(ss)){
	s <- ss[i]
#lines3d(rbind(XYZ[s,],XYZ[s,] + XYZ.[s,]*10))
#lines3d(rbind(XYZ[s,],XYZ[s,] + XYZ..[s,]*10))

lines3d(rbind(XYZ[s,],XYZ[s,] + E1[s,]*10),col=2)
lines3d(rbind(XYZ[s,],XYZ[s,] + E2[s,]*10),col=3)
lines3d(rbind(XYZ[s,],XYZ[s,] + E3[s,]*10),col=4)
}
# 接面を描く
plot3d(XYZ.st)

ss <- sample(1:length(t),50)
for(i in 1:length(ss)){
	s <- ss[i]
rgl.planes(a=E1[s,],d= - sum(XYZ[s,]*E1[s,]),alpha=0.5)

}
# 渦の接面を描く
plot3d(XYZ.st)

s <- 150
lines3d(rbind(XYZ[s,],XYZ[s,] + E1[s,]*10),col=2)
lines3d(rbind(XYZ[s,],XYZ[s,] + E2[s,]*10),col=3)
lines3d(rbind(XYZ[s,],XYZ[s,] + E3[s,]*10),col=4)

rgl.planes(a=XYZt[s,],d= - sum(XYZ[s,]*XYZt[s,]),alpha=0.5)

plot(apply(XYZt*E1,1,sum))
# 速度ベクトルは渦の接面の上
plot(apply(XYZt*E2,1,sum))
plot(apply(XYZt*E3,1,sum))