POV-Rayで描く「2重振り子」

# はじめに
これは [POV-Rayによる数学お絵かき入門 Advent Calendar 2017](https://qiita.com/advent-calendar/2017/povray4math) の15日目の記事です.

作例3回目となる今回は「2重振り子のアニメーション」を描きます.
最終的な出力は次です.

![](https://i.imgur.com/QBowZ42.gif)

![](https://i.imgur.com/lEzv6eX.gif)

# 方針
以下の方針でアニメーションを作成していきます.

1. Wolfram言語を用いて運動方程式の近似解を出力.
2. POV-Rayでテキストデータを読み込み, 連番画像を生成.
3. 連番画像をgifアニメーションに変換.

# 運動方程式を解く
次の図のような2重振り子を考えます.

![](https://i.imgur.com/5pi2pUg.png)

この振り子の運動の近似解をLagrangeの運動方程式から求めましょう.
各振り子の重心座標は次で与えられます.

```math
\begin{align}
	(x_1,y_1)
	&=\lambda_1(\sin\theta_1,-\cos\theta_1) \\
	(x_2,y_2)
	&=l_1(\sin\theta_1,-\cos\theta_1)+\lambda_2(\sin\theta_2,-\cos\theta_2)
\end{align}
```

$m_1, m_2$を各振り子の質量とすれば系の位置エネルギーは次で与えられます.

```math
\begin{align}
	U
	&=g(m_1y_1+m_2y_2) \\
	%&=-g(m_1(\lambda_1\cos\theta_1)+m_2(l_1\cos\theta_1\lambda_2\cos\theta_2))
\end{align}
```

$I_1, I_2$を各振り子の慣性モーメントとすれば系の運動エネルギーは次で与えられます.

```math
\begin{align}
	K
	=\frac{1}{2}m_1(\dot{x_1}^2+\dot{y_1}^2)
	+\frac{1}{2}m_2(\dot{x_2}^2+\dot{y_2}^2)
	+\frac{1}{2}I_1\dot{\theta_1}^2
	+\frac{1}{2}I_2\dot{\theta_2}^2
\end{align}
```

Lagrange運動方程式は$L=K-U$とおいて

```math
\begin{align}
\frac{d}{dt}\frac{\partial L}{\partial\dot{\theta_i}}-\frac{\partial L}{\partial\theta_i}=0
\end{align}
```

で与えられるのでした.
ここではWolfram言語を用いて運動方程式の近似解を求めます.
Wolfram言語を使えば, 上の微分方程式の面倒な展開をも任せて近似解を得ることが出来ます.
以下のWolfram Programming LabのTry It Yourselfから登録/課金なしで, Wolfram言語の実行が可能です.
https://lab.open.wolframcloud.com/app/

上計算と近似解を求めるためのWolfram言語のコードは次です.
ただし$m_1=m_2=1$, $l_1=l_2=1$, $\lambda_1=\lambda_2=1/2$, $I_1=I_2=m_1{l_1}^2/12$, $g=1$としました.

```mathematica
x1=\[Lambda]1{Sin[\[Theta]1[t]],-Cos[\[Theta]1[t]]};
x2=l1{Sin[\[Theta]1[t]],-Cos[\[Theta]1[t]]}+\[Lambda]2{Sin[\[Theta]2[t]],-Cos[\[Theta]2[t]]};
v1=D[x1,t];
v2=D[x2,t];
\[Omega]1=D[\[Theta]1[t],t];
\[Omega]2=D[\[Theta]2[t],t];
k=1/2 m1 v1.v1+1/2 m2 v2.v2+1/2 I1 \[Omega]1^2+1/2 I2 \[Omega]2^2;
u=m1 g x1[[2]]+m2 g x2[[2]];
l=k-u//Simplify;

l1=l2=1;
\[Lambda]1=\[Lambda]2=1/2;
m1=m2=1;
I1=(m1 l1^2)/12;
I2=(m2 l2^2)/12;
g=1;
eq1=D[D[l,\[Theta]1'[t]],t]-D[l,\[Theta]1[t]]//Simplify;
eq2=D[D[l,\[Theta]2'[t]],t]-D[l,\[Theta]2[t]]//Simplify;

s=NDSolve[{eq1==0,eq2==0,\[Theta]1[0]==\[Theta]2[0]==3,\[Theta]1'[0]==\[Theta]2'[0]==0},{\[Theta]1,\[Theta]2},{t,20}];

Pend1[L_List,\[CapitalTheta]_List]:={Sin[#],-Cos[#]}&/@\[CapitalTheta] L
Pend2[L_List,\[CapitalTheta]_List]:=Prepend[#,{0,0}]&@Table[Total@Pend1[L,\[CapitalTheta]][[1;;i]],{i,1,Length[L]}]
PendFig[L_List,\[CapitalTheta]_List]:=Graphics[Line[Pend2[L,\[CapitalTheta]]]]

Manipulate[
Show[ParametricPlot[Pend2[{l1,l2},Evaluate[{\[Theta]1[t],\[Theta]2[t]}/.s][[1]]],{t,0,T},PlotRange->{-2,2}],Evaluate[PendFig[{l1,l2},{\[Theta]1[T],\[Theta]2[T]}]/.s][[]]]
,{T,0.00001,20}]
```

実行結果を次に示します.
出力のスライダを動かせば振り子の動く様子が確認できます.

![](https://i.imgur.com/RJHxQqD.png)

これをPOV-Rayが読めるように次のように出力してファイルに保存します.
会員登録なしのWolfram Programming Labからはファイル出力が使えないため, クリップボードを経由することになります.

```mathematica

ExportString[List@Flatten@Table[Evaluate[{\[Theta]1[t],\[Theta]2[t]}/.s][[1]],{t,0,20,0.1}],"CSV"]
```

![](https://i.imgur.com/SuCUvsM.png)

結果のデータだけ欲しい方は以下からダウンロード出来ます.↓
https://gist.github.com/hyrodium/413417a5083711006db5acbcc5478aef

# POV-Rayのレンダリング
## カメラ等の設定
[第8回](https://qiita.com/Hyrodium/items/af91b1ddb8ea2c4359c2)を参考に, 次のようにカメラ等を設定します.

```
#version 3.7;
global_settings{assumed_gamma 1.0}

#macro SCS(lng,lat) <cos(radians(lat))*cos(radians(lng)),cos(radians(lat))*sin(radians(lng)),sin(radians(lat))> #end

#declare Lng=10;
#declare Lat=10;
#declare Pers=0.1;
#declare Zoom=0.45;
#declare LookAt=<0,0,0>;

#declare AspectRatio=image_width/image_height;
#declare Z=SCS(Lng,Lat);
#declare X=<-sin(radians(Lng)),cos(radians(Lng)),0>;
#declare Y=vcross(Z,X);
#declare Loc=LookAt+SCS(Lng,Lat)/(Zoom*Pers);
camera{
    perspective
    location Loc
    right -2*X*sqrt(AspectRatio)/Zoom
    up 2*Y/(sqrt(AspectRatio)*Zoom)
    direction Z/(Zoom*Pers)
    sky Y
    look_at LookAt
}
light_source{
    Loc
    color rgb<1,1,1>
}
background{rgb<1,1,1>}
```


## 外部ファイルの読み込み
[第13回](https://qiita.com/Hyrodium/items/ded2ae07269e47af83cf)を参考に外部ファイルを読み込みます.
時刻$t$を0.1刻みで0から20まで出力したので, 201組の$(\theta_1,\theta_2)$のデータを読み込みます.
ここでは読み込むついでに$(x_2,y_2)$の座標も配列`P2`に計算しています.

```Povray
#fopen MyFile "DoublePendulum.txt" read
#declare M=200;
#declare T1=array[M];
#declare T2=array[M];
#declare P2=array[M];

#declare i=0;
#while(i<M)
	#read (MyFile,t1)
	#read (MyFile,t2)
	#declare T1[i]=t1;
	#declare T2[i]=t2;
	#declare P2[i]=<0.5,sin(t1)+sin(t2),-cos(t1)-cos(t2)>;
	#declare i=i+1;
#end
```

## 振り子の描画
まず, 振り子1節分のobjectを次で描画します.

```Povray
#declare h=0.15;
#declare rod=union{
	box{<h/2,h,h>,<-h/2,-h,-1-h> pigment{rgb<0.4,0.4,0.4>}}
	cylinder{<h,0,0>,<-h,0,0>,0.075 pigment{rgb<0.05,0.05,0.05>}}
	cylinder{<h,0,-1>,<-h,0,-1>,0.075 pigment{rgb<0.05,0.05,0.05>}}
};
rod
```

![](https://i.imgur.com/CpHZGdy.png)

[第7回](https://qiita.com/Hyrodium/items/a30ff5e5304457ff102a)で解説した`rotate`, `translate`を使って次のようにします.

```Povray
#declare N=frame_number;
object{rod rotate<degrees(T1[N]),0,0>}
object{rod rotate<degrees(T2[N]),0,0> translate<0,sin(T1[N]),-cos(T1[N])> translate<0.3,0,0>}
```

ここで`flame_number`は[第9回](https://qiita.com/Hyrodium/items/089b2cf1466d0d040d61#frame_number)で解説した整数で, アニメーションのフレーム番号なのでした.
以上で次のような出力が得られます.

![](https://i.imgur.com/FqYbaWc.gif)

## 座標軸の描画
このままでは, 斜め上から振り子を見ていることが分かりにくいので, 座標軸の方向を表す矢印を描画します.
矢印を表すためのマクロを次のように書きます.

```Povray
#macro Arrow(a,b,r)
	union{
		difference{
			cylinder{a,b,r}
			cylinder{b-6*r*vnormalize(b-a),2*b-a,2*r}
		}
		cone{b-6*r*vnormalize(b-a),2*r,b,0}
	}
#end
```

これで`a`を始点, `b`を終点, `r`を半径とする矢印が描画できます.
これを用いて左下に3方向に矢印を配置します.

```Povray
union{
	sphere{<0,0,0>,0.05}
	object{Arrow(<0,0,0>,<1,0,0>,0.02)}
	object{Arrow(<0,0,0>,<0,1,0>,0.02)}
	object{Arrow(<0,0,0>,<0,0,1>,0.02)}
	translate <0,-2,-2>
}
```

これで冒頭に置いた1つ目のアニメーションが出来ました.

![](https://i.imgur.com/QBowZ42.gif)



## 軌跡の描画
次は配列`P2`を元に軌跡を表す曲線を描画します.
曲線を書くには[前回](https://qiita.com/Hyrodium/items/beae1f0562a259a3ec2b)説明した`sphere_sweep`を使います.
天下り的ですが, 以下のように書きます.

```Povray
#if(N>4)
	sphere_sweep {
	    cubic_spline
	    N+2,
	    #declare i=0;
	    #while(i<N+2)
		    P2[i],0.01
		    #declare i=i+1;
	    #end
	    pigment{rgb<1,0,0>}
	}
#end

sphere{<0.5,sin(T1[N])+sin(T2[N]),-cos(T1[N])-cos(T2[N])>,0.05 pigment{rgb<1,0,0>}}
```

以下に箇条書きで解説します.

* 離散点を滑らかに繋ぎたいので`cubic_spline`を使っています
* `cubic_spline`を使うには少なくとも4点必要なので`#if`で囲んでいます.
* `N`が小さい時は描画されない問題がありますが, 最後の`sphere`で隠れるのであまり問題になりません.

以上を合わせてレンダリングすれば以下のアニメーションが得られます.

![](https://i.imgur.com/lEzv6eX.gif)