Gnuplotで等値線付きアニメーションを描画する
これはなに
おはこんばんにちはてれじょんです。もう8月60日ですね。もう数日で8月が終わって10月になるのが信じられません。みなさんはこの夏休み何をしましたか？私は実家で皿洗いと風呂掃除をやらされてやらせて頂いて、免許をとって睡眠して食事していたら終わっていました。なんで？休み前半の進捗は結構いい感じでしたが、後半の進捗を自動車学校に破壊されたのが非常に遺憾です。まあ進捗が微妙になったのは自動車学校のせいではなく自分のせいですが……
さて本題ですが、今日はGnuplotを使ってこういうgifアニメを作っていきます。gnuplotのことは記事にするつもりはあんまりなかったのですが、結構複雑になってきて忘れそうなので記事にしておきます。
ちなみに上のグラフは地学部気象班のtwitterアカウントでも投稿した浅水方程式の数値解です。境界は自由端としてt=0で中央に正弦波を与えています。
作るぜ
作っていきます。作りたいのは山々ですが、まずインプットするデータの体裁を整えていきます。
x-y平面上の格子点におけるzの値が与えられ、それらが時間変化するような場合を考えましょう。格子点の座標は(x_i,y_j)(0<=i<=l,0<=j<=m)とします。このとき、インプットするデータは
####x,y,z####
x_0 y_0 z(x_0,y_,t_0)
x_1 y_0 z(x_1,y_0,t_0)
...
x_l y_0 z(x_l,y_0,t_0)
x_0 y_1 z(x_0,y_1,t_0)
x_1 y_1 z(x_1,y_1,t_0)
...
x_l y_m z(x_l,y_m,t_0)
x_0 y_0 z(x_0,y_,t_1)
x_1 y_0 z(x_1,y_0,t_1)
...
x_l y_0 z(x_l,y_0,t_1)
x_0 y_1 z(x_0,y_1,t_1)
x_1 y_1 z(x_1,y_1,t_1)
...
x_l y_m z(x_l,y_m,t_1)
x_0 y_0 z(x_0,y_,t_n)
x_1 y_0 z(x_1,y_0,t_n)
...
x_l y_0 z(x_l,y_0,t_n)
x_0 y_1 z(x_0,y_1,t_n)
x_1 y_1 z(x_1,y_1,t_n)
...
x_l y_m z(x_l,y_m,t_n)
#########
のように、時間ステップが次のものに移るたびに2行空行を入れておきます。これだけでOKです。
gnuplotに移り、
gnuplot>set term gif animate
gnuplot>set output "nami.gif"
で出力をgifアニメにして、適当な名前をつけます。必要に応じて
gnuplot>set xrange[0:1]
gnuplot>set yrange[0:1]
などで範囲を指定します。
で、次に等値線の設定をします。インプットするデータは離散的なデータなので、平面のデータを作らなければなりません。それを代行してくれるのがdgrid3dです。
gnuplot>set dgrid3d 50,50,4
をするとグリッドデータを作ってくれます。便利ですね。
最初2つの引数はそれぞれx,y方向の格子数、3つめの引数は格子点の距離からの重みをいくらにするかを指定する引数です。このとき、格子点からの距離の4乗の逆数で重みがかかります。
また等値線間隔を
gnuplot>set cntrparam levels incremental -0.1,0.02,0.1
で指定します。このとき等値線の最小値は-0.1,間隔は0.02,最大値が0.1になります。この他にも
・set cntrparam levels autoでは等値線の本数が自動で設定されます
・set cntrparam levels ～では等値線の本数が～本になります
最後に
gnuplot>set contour
で等値線を出力するようセットします。これで等値線の設定がオワリました。
今度はアニメーションの設定です。次のような中身のテキストファイルを(適当な名前).pltで保存します。
##
if(exist("n")==0 || n<0)n=n0
splot "output3.d" index n using 1:2:3 w l
n=n+dn
if(n<n1)reread
undefine n
##
簡単なループなのでノリで理解できると思います。
これを実行すると、splot "output3.d" index n using 1:2:3 w l
のところでインプットデータの上からnつ目のブロックが読み込まれ、t_nでのグラフが出力されます。
n0が初期値でn1が読み込むブロックの総数なので
gnuplot>n0=1
gnuplot>n1=200
gnuplot>dn=1
とか適当に設定しておきます。
最後に、
gnuplot>load "(名前).plt"
を実行します。すると現在のディレクトリにファイルが出力されるという算段です。
まとめ
####
if(exist("n")==0 || n<0)n=n0
splot "output3.d" index n using 1:2:3 w l
n=n+dn
if(n<n1)reread
undefine n
###
をplt形式で現在のディレクトリに保存、んで以下の図のように実行していくと出力されます。便利だね。
参考サイト
以下のサイトを参考にさせていただきました。ありがとうございました。
