Meepの使い方（基礎編）

まずは簡単なサンプルプログラムを使ってMeepの使い方、シミュレーション方法の基礎を学んでいきます。英語に抵抗がない方はMeepの公式チュートリアルに従って進めることを推奨します。公式が一番わかりやすくMeepの機能を説明してくれています。この記事では公式チュートリアルを参考にしつつ日本語で説明していきます。プログラミング自体が初心者の方を想定しています。

基礎編：単一モード光導波路のシミュレーション

このガラスに光を入射すると、光はガラスの中を伝わっていきます。このとき光が入射面に対して少し傾いた角度で入射した場合はどうでしょうか。途中でガラスの内壁に光線がぶつかります。ガラスと空気で屈折率が異なるため、この境界面で屈折や反射といった現象が起きます。ただしこのときの入射角が臨界角を下回っていれば、全反射が起きて光はすべてガラス内を伝搬していきます（実際には全反射面でエバネッセント場が発生し、それがガラス外部への染み出し光となりますが今回は深入りしません）。光はそのままガラスの内壁で全反射を繰り返しながら進んでいくため、光はガラス内の領域に閉じ込められた状態で伝搬していきます。従ってこの棒状のガラスは光導波路としての役割を果たします。

それではこの例と同じ状況をMeepを使ってシミュレーションしていきます。まずはシミュレーションの概要を示しておきます。

横16, 縦8 の大きさの2次元平面の領域に対して、幅1 の光導波路を横に渡し、左端の地点から光を伝搬させるシミュレーションを行います。ここでの16や8といった数値の単位は一旦考えずに、相対的な大きさの比率を表すものとして捉えます。これからこの図と同じ状況をコードで再現していくので、この図を頭に浮かべながら読み進めてください。

まずはDebianとJupyter notebookを起動します。起動方法はWindows 10 にMeepをインストールする方法 の項目８にある通りです。インストール作業を終えてそのままの状態であれば、Jupyter notebookを起動しmeep_pythonのフォルダまで進むとテストとして作成したnotebookであるUntitled.ipynbのみが見える状態になるかと思います。

このnotebookを使って最初のMeepシミュレーションを行っていきましょう。まずは混乱を防ぐためにnotebookのファイル名を"single_waveguide.ipynb"など分かりやすいものに変更しておきます。"Waveguide"は導波路という意味です。今後も新しいファイルを作ったら適切な名前に変更する癖をつけておきましょう。

名前を変更したnotebookをクリックして開きます。一つ目のセルはインストール作業時のテストとしてMeepのバージョンを確認しただけのコードであるため、削除してもそのまま放置しても構いません。ここではそのままにして次のセルにシミュレーションのコードを書いていきます。早速2つ目のセルに

import meep as mp

と入力します。

これは「meep というライブラリを読み込んで（インポートして） mp という名前で使用する」といった意味のコードです。ライブラリとは機能や関数の集合体のようなものです。簡単には、これからMeepの機能を使いたいので、まず最初にMeepを呼び出したといったイメージです。

続いてこれは任意ですが、コードの意味を忘れないためにPythonのコメント機能を使ってコードの説明を書いておきます。半角"#"のあと「meepライブラリのインポート」などとコメントを入れておきます。"#"以降に書いたものはコードとして認識されないため、好きなように補足説明などを入れておくことができます。

続いてシミュレーションの本体に入っていきます。まずはシミュレーションを行う空間領域の大きさを決めるコードを書きます。改行した上で、

cell = mp.Vector3(16,8,0)

と書きます。

これは、geometry という変数を作成して、それに [...] で囲まれたオブジェクトを代入しています。Pythonにおいて[...] はリストと呼ばれ、配列を表します。

リストの中身を見ていきます。リストの中にはmp.Block(...) なる配列要素が一つあり、 (...) の中にさらにコードが並んでいます。mp.Block(...) はMeepのライブラリを使用してBlock（直方体）の構造（今回の光導波路）を作るコードです。これはPythonにおける関数の形になっていて、 (...) の中に関数の引数（変数）を書いていくことでBlockのサイズや性質を決めていきます。一つ目の引数 mp.Vector3(mp.inf,1,mp.inf)は先ほども出てきた3次元ベクトルを使ってBlockのサイズを定義するものです。mp.infは無限大を表すので、x と z 軸方向が無限大で y 軸方向が１というサイズです。x 軸方向が無限大と書いても、cell で領域サイズを (16,8,0) と決めたので自動的に横幅は16となります。z 軸方向は今回2次元平面上のシミュレーションであるため数値は何でもよいです。

続くcenter=mp.Vector3()はBlockの中心座標を定義するものです。Vector3()の括弧の中に何も書かない場合はデフォルトの(0, 0, 0)が自動的に入るため、Blockの中心座標はシミュレーション領域の原点となります。上で示したシミュレーションの概要図の通り、原点は領域の中心にあります。ちなみにPython関数のデフォルトの引数は関数ごとに定義されており、「引数を書かない場合デフォルトの引数が自動で入る」というのはPythonが持つ特性です。

3つ目の引数material=mp.Medium(epsilon=12)はBlockの材質を決めるものです。Meepのライブラリを使用してMediumなる材質を定義しています。Mediumは直訳すると光を伝搬させる「媒質」です。より具体的には「誘電体」というものをシミュレートしています。誘電体には「誘電率 (epsilon)」、「透磁率 (mu)」、「屈折率 (index)」などを設定することができ、この設定によって光の伝搬する際の挙動が変わります（設定項目について詳しくはhttps://meep.readthedocs.io/en/latest/Python_User_Interface/#medium ）。今回はepsilon=12と書き、誘電率が12と定義したことになります。透磁率は書かない場合自動的に１となります。ちなみに屈折率は誘電率と透磁率を使って

と表されるため、屈折率が約3.46と同じ意味になります。

次は光源（光の発生源）の性質や位置を定義するコードに移ります。

sources = [mp.Source(mp.ContinuousSource(frequency=0.15),

component=mp.Ez,

center=mp.Vector3(-7,0))]

と書きます。

geometry のときと同様で、リストの中に関数が一つあり、そのリストを引数sourcesに代入しています。関数mp.Source(...)は光源を定義するものです。関数の引数を一つずつ見ていきます。

mp.ContinuousSource(...)は光が連続波であることを意味し、単一周波数の波を連続的に発生させます。(...)内のfrequency=0.15によって、その連続波の周波数を0.15に設定しています。Meepでは周波数は2πcの単位で表現されるため、周波数の逆数をとれば波長に直すことができ、波長で表せば約6.67です。

component=mp.Ezによって光の偏光を決めています。ここでは、mp.Ezによって電場が z 軸方向に振動している電磁波、すなわちTE波を定義しています（偏光について詳しくはWikipediaページを参照）。center=mp.Vector3(-7,0)によって光源の位置を決めています。(x, y, z) = (-7, 0, 0) の位置、すなわち領域の左端としています。

今回は行っていませんが、4つ目の引数 size=mp.Vector3(...) を加えることで光源のサイズを指定することもできます。指定しない場合はデフォルトのサイズ (0, 0, 0) となり、点になります。

ここまでのコードで、シミュレーションの主要な部分である光導波路と光源の設定ができました。続いて、シミュレーションの境界条件を設定します。

pml_layers = [mp.PML(1.0)]

と書きます。

これも変数にリストを代入する形になっています。リスト内のmp.PML(1.0)は厚み1.0のPerfectly matched layers (PML) なるものを定義しています。PMLとは、光を完全に吸収するシミュレーション上の仮想の物質です。シミュレーションの概要図にあるようにcellの内壁に厚み1.0のPMLがぐるっと張り巡らされるため、cellの端に到達した光は完全に吸収されます。cellの端で反射するなどしてシミュレーション結果に不必要な影響を及ぼすと困るので、このような境界条件にしています。

次は、シミュレーションにおける計算の解像度を設定します。

resolution = 10

と書きます。

resolution（解像度）という変数を定義し、10を代入しました。この数値は、空間と時間をどのくらい細かく分割して計算を行うかに影響します。この場合、１の長さを10個の区間に分けて電場や磁場を計算していきます。resolution = 20とすれば20個の区間に分けます。下の図のように、自然界では連続的な波を有限の区間に分割してシミュレートするため、resolutionが大きいほど計算の精度も高くなると予測できます。

ただしresolutionを増やすほど、計算に掛かる時間も増していきます。空間と併せて時間の分割数（計算のステップ数）もresolutionの影響を受けるため、例えばresolutionを2倍にすると、x 方向に2倍、y 方向に2倍、時間の分割数も2倍に増えるため、8倍も計算量が増えることになります（必要なメモリの容量も増えます。PCのメモリ容量を超えるとPCがフリーズするため特に注意が必要です）。そのため精度と計算時間のバランスを常に考える必要がありますが、典型的にはresolutionは10～20で行うことが多いです。10未満にすると解像度が低過ぎて本来と異なる結果が得られる場合もあるため避けるべきですが、20より大きくしても結果がほとんど変わらないことが多いためです。

次は、シミュレーションを実行するための準備として、シミュレーションオブジェクト（変数）を作成します。

sim = mp.Simulation(cell_size=cell,

boundary_layers=pml_layers,

geometry=geometry,

sources=sources,

resolution=resolution)

と書きます。

少々長いですが、やっていることは単純で、関数mp.Simulation(...)の引数としてこれまで定義してきたcellやsourcesなどの変数をすべて代入しています。シミュレーションの情報をひとまとめにしているイメージです。これとセットになるコードが次のもので、

sim.run(until=200)

と書きます。

run(...) によってsim に格納されたシミュレーションを実行するコードです。run(...)の引数によってシミュレーションを実行する時間の長さや何を出力するかを指定します。今回はシンプルに until=200 とだけあります。時刻200までシミュレーションを実行せよという意味です。

これでシミュレーションのコーディングが完了しました。早速実行してみましょう。notebook上でキーボードのshift + Enter キーを同時押しして実行します。数秒～十数秒で以下のような実行結果が出力されます。

この出力は、前半部分でシミュレーションの条件が示され、後半部分でシミュレーションが無事完了したことが示されています。"Computational cell is 16 * 8 * 0 with ..." などのようにコーディングした通りの設計になっていることが確認できます。

次はシミュレーションの具体的な結果を画像として出力してみます。構造のプロットと、電場分布（光の分布）のプロットを2段階に分けて行います。

まずは結果を処理したりプロットするために必要なPythonライブラリを2つインポートします。notebookの次のセルをクリックして、

import numpy as np

import matplotlib.pyplot as plt

と書きます。

numpy はPythonでデータ処理や数値計算を行う際によく使用されるライブラリです。Matplotlib は Pythonで図を描画する際によく使用されるライブラリです。Matplotlib は巨大なライブラリであるため、全体をインポートせずに今回のデータプロットに必要なサブパッケージである pyplot のみをインポートしています。

続いて以下の一連のコードを書きます。

eps_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Dielectric)

plt.figure()

plt.imshow(eps_data.transpose(), interpolation='spline36', cmap='binary')

plt.axis('off')

plt.show()

キーボードのshift + Enter キーを同時押しして実行すると以下のような画像がnotebookに表示されます。

このコードによって、シミュレーションを行った誘電体の構造をプロットしました。画像の黒い棒が光導波路です。

コードの意味を説明します。一行目のeps_data = ...ではまずシミュレーション結果から「シミュレーション領域全体の誘電率の値」を配列として取得しています。sim という変数にシミュレーション結果が格納されているので、sim.get_array(...) という関数によって(...) の中で指定した要素（配列）を sim から取り出してきます。(...) の中を見ると、center=mp.Vector3() によって取得する配列の中心座標は (0,0,0)、size=cell によって取得する配列のサイズを cell（すなわちシミュレーションの領域全体）と指定し、最後の component=mp.Dielectric によって誘電体構造を表す配列を取得せよ指定しています。今回のシミュレーションでは誘電体構造に誘電率(epsilon=12)という性質を与えたので、シミュレーション領域全体の誘電率の値が配列として取得されます。これによって eps_data という変数に以下のような配列（行列）が代入されたとイメージしてください。

行列の要素数が 16*10, 8*10 である理由は、cell = (16, 8, 0) にresolution = 10 が掛かるためです。

続く4行のコードはこの配列をプロットするためのものです。

plt.figure() によってまずfigureオブジェクトを作り、図を描画する土台を作ります。

plt.imshow(eps_data.transpose(), ...)は引数で指定した2次元配列をカラーマップとして描画する関数です。引数に eps_data とあるため、上記の配列が描画されます。.transpose() とくっついていますが、これは配列を転置する関数で、Meepで作成した構造を配列に変換するとき反転してしまう仕様を補正するためのものであるため深い意味はありません。残りの引数は描画時の見た目を決めるオプションで、interpolation = 'spline36' によって画像を拡大しても滑らかに表示できるようになり、cmap = 'binary' によって2値のカラーマップになります。この場合 [1 , 12] の誘電率の値を [0, 1] に変換するということです。

plt.axis('off')は、デフォルトではプロットに縦軸と横軸が自動的に入るため、それを消すためのコードです。

最後にplt.show() によってプロットを画面に表示します。

シミュレーションで設定した通りの構造を図示できたので、次は電場の分布をプロットして光がどう伝搬したかを確認しましょう。以下の一連のコードをnotebookの次のセルに書き、実行します。

ez_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Ez)

plt.figure()

plt.imshow(eps_data.transpose(), interpolation='spline36', cmap='binary')

plt.imshow(ez_data.transpose(), interpolation='spline36', cmap='RdBu', alpha=0.9)

plt.axis('off')

plt.show()

実行結果を見ると、赤と青のカラーマップが確認できます。これが光を表していて、z 軸方向に振動する電場 Ez の分布です。赤が負の値を表し、青が正の値を表します（色が濃いほど絶対値が大きく、白は0です）。画面に対して垂直な z 軸方向に振動している波を表しています。光導波路の部分に電場が集中していることがこの図から読み取れます。

コードの内容は、一つ前のセルと見比べると分かりやすいはずです。ez_data という変数に電場を指す mp.Ez の配列を代入して、2つ目のplt.imshow(...) でそれを描画しています。cmap='RdBu' によって赤と青のカラーマップになっています（matplotlib のカラーマップコードを調べれば色を変更することも可能です）。alpha=0.9 はアルファ値、すなわち色の不透明度です。わずかに透明で後ろが見えるようにしている理由は、一つ上のplt.imshow(...) 関数にあります。1つ目のimshow() によって光導波路構造をプロットして、2つ目のimshow() でその上に重ねるようにして電場をプロットしました。

ちなみにこの電場分布は、光源から光が発生してから t=200 後の様子です。シミュレーション実行時に sim.run(until=200) と指定したためこの時刻の結果が得られています。試しにsim.run(until=25) などのように短い時刻に変更してシミュレーションを再び実行してみると、光が導波路の中間あたりまで来ている時刻の様子を見ることができます。この時刻の単位については別の記事で解説します。

以上でMeepの使い方の基礎編は終了です。Meepにおける基本的なプログラミングの書き方と、その結果を図でプロットする方法を紹介しました。