(2014/10/16作成,2016/02/26追記)
ここではグラフ描画ができるgnuplot(wgnuplot)を用いて光渦を描く方法について説明します。
このページでは数式を表示するためにMathJaxを利用しています。javascriptがONでありかつこの(MathJax公式ページ)リストに載っているWebブラウザーを利用しているときのみに数式が正しく表示されます。ご了承ください。
ここでは、光渦を例にとって、gnuplotを用いた3次元グラフ描画法について述べます。具体的には、ラゲールガウスビームのビーム断面内強度分布と位相分布を描きますが、それ以外のグラフを描きたい場合でもこのページは有用であると考えます。
このページで描けるグラフは3種類です(図1)。一つ目はメッシュで表現した3次元グラフです。二つ目は、その3次元グラフの高さに応じて色付けを行ったグラフです。三つ目は、二つ目のグラフを真上からみたグラフ(密度プロット)です。
図1 このページを読んで描けるグラフ
gnuplotで数式を描く流れは以下のとおりです。
特に2番の「gnuplotに入力する式の形を考える」が重要です。すなわち、1番で求めた描きたいグラフからgnuplotに入力する式の形の変換に慣れれば、様々な数式をグラフとして描けるようになります。
gnuplotでグラフを描くには、gnuplotに入力する関数を考えなくてはなりません。このページでは、例として、光渦の強度分布と位相分布を3次元プロットで描きます。
まず、光渦の中でもラゲールガウスビーム(LGビーム)と呼ばれるラゲールガウスモードで記述するビームの強度分布と位相分布を描くことにします。ラゲールガウスビームには\(l, p\)の二つの量子化された指数があり、その電場振幅\(u^\mathrm{LG}_{l,p}\)は\(kz-\omega t\)の位相を取る系(時刻が進むに従い\(+z\)の向きにビームが伝播する)において
\begin{multline} u^\mathrm{LG}_{l,p}(r,\phi ,z) = A\left (\frac{\sqrt{2}r}{w} \right )^{|l|} L^{|l|}_{p}\left ( \frac{2r^2}{w^2} \right )\frac{w_0}{w} \\ \times \exp \left [ -\frac{r^2}{w^2} \right ] \exp \left [ i\left ( k\frac{r^2}{2R} +l\phi -(2p+|l|+1)\psi \right ) \right ] \tag{1} \end{multline}と表されます[1]。ここで、\(r, \phi, z\)は円筒座標系の座標で、\(A\)は振幅、\(k\)は波数を示します。そして、\(w_0, w, R, \phi\)(ビームウエスト半径、ビーム半径、波面の曲率半径、グイ位相)は
\begin{align} w_0&=\left ( \frac{2z_0}{k}\right )^{1/2} \tag{2}\\ w &=w(z)= \left ( \frac{z^2+z_0^2}{kz_0} \right )^{1/2}=w_0\left ( 1+\frac{z^2}{z_0^2} \right )^{1/2} \tag{3}\\ R &=R(z)= \frac{z^2+z_0^2}{z} \quad(z\neq 0) \tag{4}\\ \psi &= \mathrm{arctan} \left ( \frac{z}{z_0} \right ) \tag{5} \end{align}であり、\(z_0 \)はレイリー長と呼ばれています。さらに、\(L^{k}_n(\xi)\)(\(k\):整数,\(n\):0以上の整数)はラゲール陪多項式と呼ばれ、低次の\(n = 0,1,2\)においては、
\begin{align} L^{|l|}_0(\xi) &= 1 \tag{6}\\ L^{|l|}_1(\xi) &= -\xi + 1+|l| \tag{7}\\ L^{|l|}_2(\xi) &= \frac{1}{2}\xi^2 - (2+|l|)\xi + \frac{(2+|l|)(1+|l|)}{2} \tag{8} \end{align}と記述されます[2]。
式(1)から式(8)までの情報があれば原理的にグラフを描くことはできますが、gnuplot上でそれを行うのは、骨が折れる作業です。そこで、以下のような条件を付けて式を簡単にします。
これらの条件により、電場振幅\(u^\mathrm{LG}_{l,p}\)は
\begin{equation} u^\mathrm{LG}_{l,p}(r,\phi ,z) = u^\mathrm{LG}_{1,0}(r,\phi ,z=0) = \sqrt{2}A\tilde r \exp \left ( -\tilde r^2 \right ) \exp \left ( i \phi \right ) \tag{9} \end{equation}まで、簡単になります。(もう少し条件を緩くした場合を発展編で取り扱う予定です。)
ビーム断面内強度\(I(r,\phi,z=0)\)は電場振幅の絶対値の自乗に比例し、ビーム断面内位相\(\Phi(r,\phi,z=0)\)は電場振幅の位相であることから、
\begin{align} I(\tilde r,\phi,z=0) &= \left | u^\mathrm{LG}_{1,0}(\tilde r,\phi ,z=0) \right |^2 = 2|A|^2 \tilde r^{2} \exp \left ( -2\tilde r^2 \right ) \tag{10}\\ \Phi(\tilde r,\phi,z=0) &= \phi \tag{11} \end{align}となります。(便宜上、ビーム断面内強度は\(I(r,\phi,z=0)\propto \)とすべきですが、\(I(r,\phi,z=0)= \)にしています。)式(1)と比べると非常に簡単になりました。
さて、式(10),式(11)は極座標で記述されていますが、これを直交座標に変換します。(gnuplotでも極座標プロットはできますが、便宜上式を直交座標で表します。)\(r = \sqrt{x^2+y^2}, \phi = \mathrm{arctan}(y/x)\)であるので、
\begin{align} I(\tilde x, \tilde y,z=0) &= 2|A|^2 (\tilde x^2+\tilde y^2) \exp \left ( -2(\tilde x^2+\tilde y^2) \right ) \tag{12}\\ \Phi(\tilde x, \tilde y,z=0) &= \mathrm{arctan}\left (\frac{\tilde y}{\tilde x} \right ) \tag{13} \end{align}と書き変えることができます。(ただし、\(\tilde x = x/w_0, \tilde y = y/w_0\))
また、強度分布の最大値を1に規格化することを考えると、\(|A|^2= \frac{w(z)}{w_0}\left ( \frac{e}{|l|} \right )^{|l|} = e\)とすることができます。
これで、gnuplotに打ち込む数式の準備は完成です。
gnuplotでグラフを描くには、gnuplotを起動し、ターミナル画面に命令を打つ必要があります。今回は、そのような方法ではなく、命令のひと固まりをテキストファイルに書き保存して、gnuplotに読み込ませる方法をとります。この方法でも、入力するべき命令は同じです。(すなわち、これから紹介するファイルの内容を1行1行gnuplotに手入力すれば全く同じことができます。)
まず、メモ帳やTeraPad、秀丸などのテキストエディタを用いてコード1の内容を作ってください。ファイル名は"code1.plt"とします。(拡張子がpltであれば名前は何でもよいです。)これは、図1の中のメッシュで表現した強度分布の3次元グラフを作る内容になっています。。
set xrange[-2:2] set yrange[-2:2] set xtics 1 set ytics 1 set isosamples 100,100 f(x,y) = 2 * (x**2 + y**2) * exp(1 - 2*(x**2 + y**2)) splot f(x,y)
これをgnuplotで読み込んでみましょう。まずは、gnuplotを起動します。gnuplotのツールバーに「開く」というボタンがあります(図2)。このボタンをクリックして、先ほど作成したpltファイル("code1.plt")を開いてください。そうすると、新たなウィンドウが開いて、グラフが表示されます(図3)。
図2 「開く」ボタンの位置
図3 表示されたグラフ
コード1の内容を詳しく説明します。コード1は主に、グラフの描画設定を行う部分(1-5行目)、描画する式を設定する部分(8行目)、グラフを描画する部分(9行目)の3つに分かれています。
まずグラフの描画設定を行う部分(1-5行目)について説明します。1,2行目の"set xrange[-2:2]","set yrange[-2:2]"はx,y軸の表示範囲をそれぞれ-2から2に指定する命令です。3,4行目の"set xtics 1","set ytics 1"はx,y軸の見出しがつく大きい目盛の間隔を1に指定しています。5行目の"set isosamples 100,100"はメッシュの細かさを指定しています。100,100という数字はx方向に100本、y方向に100本のメッシュを用いて3次元グラフを描画せよということを意味しています。
8行目では"f(x,y) = 2 * (x**2 + y**2) * exp(1 - 2*(x**2 + y**2))"と式(12)の内容を関数f(x,y)に代入しています。(数式について詳しく知りたい場合は、gnuplotの説明書[3](バージョン4.4の説明書には日本語版があります。)の第13章「式」を参考にしてください。そのほかのwebサイトでも有用な情報を得ることができます。)
最後の行では"splot f(x,y)"と、6行目で作った関数を用いてプロットする命令を出しています。(この時点で、グラフ描画が行われて、図3のグラフが現れます。)
図1の3次元グラフの高さに応じて色付けを行ったグラフを描くpltファイルの内容をコード2に示しました。コード1との違いは9行目に"with pm3d"が加わったのみです。
set xrange[-2:2] set yrange[-2:2] set xtics 1 set ytics 1 set isosamples 100,100 f(x,y) = 2 * (x**2 + y**2) * exp(1 - 2*(x**2 + y**2)) splot f(x,y) with pm3d
最後に、密度プロットを行うpltファイルの内容をコード3に示します。
[注意]
set xrange[-2:2] set yrange[-2:2] set xtics 1 set ytics 1 set isosamples 100,100 set size square set pm3d map f(x,y) = 2 * (x**2 + y**2) * exp(1 - 2*(x**2 + y**2)) splot f(x,y)
コード1との違いは、6,7行目です。6行目はグラフのサイズを指定する命令で、描画領域を正方形にするように指定しています。(これを指定しない場合では、横長なってグラフが表示されてしまいます。)7行目では、密度プロットのモードする命令を行っています。7行目でpm3dのモードに既になっているので9行目ではwith pm3dを打つ必要はありません。
[注意]コード3の7行目で"set pm3d map"と密度プロットモードにしているために、コード3を実行した後にコード1やコード2などを実行すると、表示されるグラフが全て密度プロットになってしまいます。これを回避するためには、一度密度プロットのモードが抜けなければなりません。そのためには、コード4の内容をgnuplotのターミナルに打ち込んでください。
密度プロットモードは実質的には、pm3dモードで描画した3次元グラフを真上(0,0)から見ることをしています。"unset pm3d"でpm3dモードを解除し、"unset view"でデフォルトの位置(60,30)から3次元グラフを見るようにしています。
最後に、位相分布の表示方法について説明します。コード1,2,3の8,9行目をコード5のように書き換えてください。
8行目のatan2(y,x)という関数ですが、数学的には、atan(y/x)と同じ意味です。しかし、コンピュータの世界では、atan2(y,x)とatan(y/x)は少し異なった振る舞いをします。興味のある方は、8行目をg(x,y) = atan(y/x)としてみてください。
さて、もう少し発展させてみます。ここでは、先ほどのメッシュ3次元グラフの軸に名前を付けて、強度分布、位相分布をそれぞれPNGファイルで保存する方法について説明します。
コード6の内容をpltファイルとして作成し保存してください。これをgnuplotで実行すると、強度分布、位相分布のメッシュ3次元グラフがそれぞれ"I.png","phi.png"として保存されます。
set xrange[-2:2] set yrange[-2:2] set xtics 1 set ytics 1 set isosamples 100,100 set xlabel "x/w_0" set ylabel "y/w_0" set zlabel "I" f(x,y) = 2 * (x**2 + y**2) * exp(1 - 2*(x**2 + y**2)) g(x,y) = atan2(y,x) set term png set output "I.png" splot f(x,y) title "" set zlabel "Phi" set output "phi.png" splot g(x,y) title ""
コード6の6-8行目と13-17行目以外はこれまでに説明をしているのでここでは説明しません。
6-8行目では、x,y,z軸に軸の名前を設定する命令をしています。set [軸名]label "[軸の名前]"と打つことで[軸名]軸の名前を[軸の名前]にすることができます。
13行目では出力先をPNGファイルに変更しています。実は、gnuplotは起動時のデフォルトとして、wxtというウィンドウへグラフが表示される設定("set term wxt"が命令された状態)になっています。"set term [出力先]"と打つことによって、デフォルトのwxtから出力先を変えて様々なファイル形式にグラフを出力することができます。他の出力先についてはgnuplotの教科書[3]の第83章を参考にしてみてください。出力先を変えてグラフ描画を行っているため、前とはちがって、グラフを表示するウィンドウ(wxt)は現れません。
14,17行目では出力ファイルのファイル名を指定しています。適宜お好みのファイル名に設定してください。
15,18行目のsplotは前のコードにもありましたが、表示する関数の後にtitle ""と書かれています。これは、凡例に書く文字を無しにするという意味です。もし、title ""がなければ、I.pngには、グラフの右上にf(x,y)という文字が出ます。しかし、これは3次元プロットにおいてはあまり意味をなさないので文字を消しています。もちろん、splot f(x,y) title "abc"と書くと、f(x,y)の代わりにabcが書かれます。
16行目のset zlabel "Phi"を書き忘れると、位相分布のz軸の名前が"I"のままになってしまいます。
メッシュで表現した3次元グラフ以外のグラフについても、コード2,3を参考に適宜コード6を修正すればpngファイルに強度・位相分布を出力できます。
作製中・・・
gnuplotの使い方を勉強するには、gnuplotに関する書籍を読む、インターネット上のgnuplotに関するwebページを読む、説明書を読むなどの方法があります。
インターネット上のgnuplotに関するwebページは数多ありますが、例えば、GNUPLOT not so Frequently Asked Questionsは物理系で必要となるグラフ作成の全てが網羅されていると言っても過言ではない程の詳細な解説が載っています。このサイトの中の「gnuplot入門」をまずは始めてみるとよいかもしれません。(注:GNUPLOT not so Frequently Asked Questionsの本家サイトは閉鎖されていますが、有志によりミラーリングがなされ、このリンク先もそもミラーサイトの一つです。)
