二重スリット

二重スリットのシミュレーションをしてみましょう。 スリットの材質は完全導体とします。 Main.ctlというファイル名で次のようなコードを保存してください。
(set! geometry-lattice (make lattice(size 20 20 no-size)))
(set! pml-layers (list(make pml(thickness 1))))
(set! geometry
 (list
  (make block (center 0 6.75)(size 0.1 6.5 infinity)
   (material(make perfect-metal))
  )
  (make block (center 0 -6.75)(size 0.1 6.5 infinity)
   (material (make perfect-metal))
  )
  (make block(center 0 0)(size 0.1 3 infinity)
   (material (make perfect-metal))
  )
 )
)
(set! sources
 (list
  (make source
   (src (make continuous-src(frequency 0.5)))
   (component Ez)
   (center -9 0)
   (size 0 20)
  )
 )
)
(set! resolution 40)
(run-until 30
 (at-beginning output-epsilon)
 (to-appended "Ez" (at-every 0.2 output-efield-z))
)
コードの要素の意味は次の通りです(参考)。
(set! <変数> <値>)
<変数>に<値>を代入します。
geometry-lattice
計算領域のサイズを定めるグローバル変数。型はlatticeクラスです。
(make <クラス> (<プロパティ1> <値1>) (<プロパティ2> <値2>) ...)
<クラス>を作成し、<プロパティ1>に<値1>、<プロパティ2>に<値2>、...を設定します。
pml-layers
境界条件であるpml層の設定を定めるグローバル変数。型はpmlクラスのlistです。
geometry
誘電率、透磁率、伝導率などを定めた物体を配置するグローバル変数。型はgeometry-objectクラス(例えば後述のblock)のlistです。
(list <オブジェクト1> <オブジェクト2> ...)
オブジェクト1、オブジェクト2、...が格納されたリストをreturnします。
block
平行六面体のクラス。
(size <x軸方向の長さ> <y軸方向の長さ> <z軸方向の長さ>)
blockの各辺の長さを定めるプロパティ。
(center <x> <y> <z>)
オブジェクトの中心を座標(<x>,<y>,<z>)に定めるプロパティ。
material
オブジェクトの誘電率、透磁率、伝導率などを定めるクラス。
perfect-metal
materialの完全導体(抵抗率0)のプロパティ。
sources
電磁場の発信源を設定するグローバル変数。
continuous-src
$e^{ -i\omega t}$で発振する連続波を定めるプロパティ。
frequency
周波数を定めるsourceのプロパティ。MEEPは(光速)=1となる単位系を採用しているため、周波数=<値>のとき、波長は$<値>^{-1}$となります。
component
sourcesの成分を定めるプロパティ。例えばEzとすると、z方向に直線偏光した電場となります。
size
sourceの大きさを定めるプロパティ。
center
sourceの中心を定めるプロパティ。
resolution
単位長さ当たりのグリッド数を定めるプロパティ。
(run-until <数値> <関数>)
<数値>ステップだけシミュレーションを実行し、各ステップごとに<関数>を実行します。
(at-beginning <関数>)
シミュレーションの1ステップ目だけに<関数>を実行します。
(to-appended <ファイル名> <関数>)
シミュレーションの各ステップごとに<関数>を実行し、生成されたすべてのファイルを1つのファイルにまとめて、名前の末尾に<ファイル名>を追加後、保存します。
次のコマンドでシミュレーションが始まります。
meep Main.ctl
電磁波のシミュレーション結果Main-ez.h5と誘電率のデータMain-eps-000000.00.h5がMain.ctlと同じディレクトリに生成されます。 コマンドh5topngによって、HDF5ファイル(拡張子h5)をpngに変換できます。
h5topng Main-ez.h5 -S3 -Z -c dkbluered -t 0:149 -R -a yarg -A Main-eps-000000.00.h5 -m -1.4 -M 1.4
ez.h5がpngに変換されます。 オプションは次の通りです(参考)。
-S3
線形補間によって、出力する画像のピクセル数を入力の3倍にします。
-Z
カラースケールの中心(白)が0になるように設定します。
-c colormap
スケールの表示にcolormapを使用します。ここではmapをdkbluered(暗青色(負)-白(零)-暗赤色(正))にしています。
-t start:end
ここではh5ファイルはx、y、tの3つの次元を持っているので、二次元の画像にするには、何かの次元でスライスする必要があります。-t start:endとすると、データをt=start,start+1,...,endでスライスし、end-start+1個に分割します。
-R
-t start:endによって分割されたデータのスケールを揃えます。スケールの最大値、最小値はすべてのデータのうちの最大値と最小値に自動的に設定されます。最大値と最小値の絶対値が等しいとは限りません。-Rを設定しない場合、スケールの最小値と最大値が各スライスごとに独立に設定されます。
-a colormap -A file
インプットのデータセットと同じサイズのfileをcolormapのカラースケールで半透明にして重ねます。ここではcolormapとしてグレースケールのyargを用いて、誘電率のマップを電磁波のマップに重ねています。
-m min
スケールの最小値をminに設定します。データに含まれる値の最小値はオプション-vを設定して実行すると知ることができます。
-M max
スケールの最大値をmaxに設定します。データに含まれる値の最大値はオプション-vを設定して実行すると知ることができます。
convert *.png ez.gif
を実行すると、画像がgifに変換され、下のような動画になります。 gifアニメーションです。 400×30ピクセルのスケールバーが次のコマンドで生成できます(参考)。
echo "0 1" | h5fromtxt bar.h5
h5topng -X 200 -Y 30 -c dkbluered bar.h5
生成されたスケールバーは次の画像です。 スケールバーの画像です。 これは、0と1の2つしか数値を持たないHDF5ファイルを-X 200 -Y 30によって、x方向に200倍、y方向に30倍に引き伸ばしてdkblueredで着色することで生成しています。 経過時刻のgif上での表示は、ファイル名を経過時刻に編集してからPhotoScape Xの一括編集で"[name]sec"を挿入することで可能です。