[2018/08 started by Maruoka and Maeda]
MEEPに関すること分からないことをmaru@flex.phys.tohoku.ac.jpに送っていただければ答えます。丸岡の居室（842号室）に来てもらってもOKです。
*目次 [#t1296280]
#contents
*MEEP [#g2bb328a]
**MEEPとは [#fbfd5425]
時間領域差分法（FDTD法）による電磁場のシミュレーションを行うプログラムです。
Simpetusによってオープンソースで開発されており、フリーで使えます。
商用ソフトではLumerical(FDTD)やCOMSOL(FEM)、富士通のPoytingがありますが、かなり高価です。
例えば、富士通のPoyntingは[https://pr.fujitsu.com/jp/news/2005/10/3.html 350〜万円]するようです。
フリーのFDTDソフトとしてはMEEPがほぼ唯一の選択肢だと思います。
オープンソースであることやCUIベースであることから、使用者による機能の追加が簡単で自由度が高いのが特徴だと思います。
**PyMeepのインストールと実行 [#fc543fa9]
Linux環境でPythonを用いてMEEPが使えるようになることを目標とします。
MEEPはLinuxとMacでしか動作しないため、Windowsを使用する場合はWSL(Windows Subsystem for Linux)をインストールしてください。
WSLのインストール方法は[[こちら:http://www.atmarkit.co.jp/ait/articles/1608/08/news039.html]]。~
WSL上ではLinuxと同じ手順で環境が構築できます。~
***インストール [#p4cd7dd9]
まずはAnacondaを用いてPythonの実行環境を構築します。([[参考:https://qiita.com/t2y/items/2a3eb58103e85d8064b6]])~
[[Anacondaのサイト:https://www.anaconda.com/download/]]にてPython3.x系のインストーラーをダウンロードしてください。~
インストーラーはBourne Shell(sh)のシェルスクリプトで書かれているので以下のように実行してください。
$ bash [ダウンロードフォルダ]/Anaconda3-5.2.0-Linux-x86_64.sh
次にPyMeepをインストールします。([[参考:https://meep.readthedocs.io/en/latest/Installation/]])
以下のコマンドを実行してください。
$ conda create -n mp -c chogan -c conda-forge pymeep
これでPyMeepのインストールは完了です。
次のエラーが出るときは、~/anaconda3/binにパスを通してください。
bash: conda: command not found
tcshを使用している場合は次のコマンドを用いると、~/anaconda3/binにパスが通ります。
setenv PATH $PATH\:<path to anaconda3/bin>
bashを使用している場合は次のコマンドを用いてください。
export PATH=$PATH\:<path to anaconda3/bin>
相対パスではなく絶対パスを使用することに注意してください。
***実行 [#padb0b5f]
次のコマンドでPyMeepの仮想環境をアクティベートしてください。
$ source activate mp
"-bash: activate: No such file or directory" とエラーが出るときは
export PATH=~/anaconda3/bin:$PATH
を実行してください。
flexやtube（齋藤グループのワークステーションの名前）ではPyMeepの仮想環境をアクティベートしようとすると
#!/bin/sh
_CONDA_ROOT="/home/students/maru/anaconda3"
\. "$_CONDA_ROOT/etc/profile.d/conda.sh" || return $?
_conda_activate "$@"
というエラーが出ます。
これはログインシェルがtcshであるにも関わらず、shのシェルスクリプトをsourceしているのが原因です。
tcshからはtcshスクリプトしかsource出来ません。
bash
と打って、ログインシェルをbashに変更してからPyMeepの仮想環境をアクティベートしてください。
実行はhogehoge.pyの部分を任意のファイルに置き換えてください。
(mp) $ python hogehoge.py
サンプルソースはGithub等で入手できるので[[このソースコード:https://github.com/stevengj/meep/blob/master/python/examples/wvg-src.py]]で試してみてください。~
PyMeepの仮想環境を無効化するには以下のコマンドを実行してください。~
(mp) $ source deactivate
**MEEPの並列化 [#k1f4cae9]
***インストール([https://meep.readthedocs.io/en/latest/Installation/#mpi 参考]) [#v7b0b7df]
並列処理の可能なPyMeepをインストールするにはAnacondaのインストールを行ったうえで次のコマンドを実行してください。
conda create -n pmp -c conda-forge pymeep=*=mpi_mpich_*
!注意! 以前は以下のコマンドでのインストールが標準でしたが、変更されたようです(2019/05/23)。こちらはバージョン1.6から更新が止まっているので上のコマンドでインストールしてください。
conda create -n pmp -c chogan -c conda-forge pymeep-parallel
***実行 [#j3e26633]
$ source activate pmp
にて並列処理をサポートしたPyMeepの仮想環境をアクティベートしたのち、
$ mpirun -np 4 python <script_name>.py
を実行してください。
anacondaの一部のファイルがbashで書かれているため、tcshからはPyMeepの環境をアクティベートできません。
エラーの出ているbashで書かれたスクリプトをtcshに書き直すか、shellをbashに変更してからアクティベートしてください。[https://teratail.com/questions/144776 参考]
オプションの" -np 4 "は4プロセッサで並列処理することを指定しています。
並列に計算するプロセッサの数は適切な値に変更してください。
***別の並列化 [#e415ae6a]
独立なプログラムを並列に計算することで計算時間を削減することもできます。（Embarassingly parallelといいます）。
meep <file-name>.ctl &
とすると、バックグランドでプログラムを実行できます。
例えば、次のようなシェルスクリプトを使用しています。
for((i=0;i<90;i+=5))
do
meep incident-angle=${i} Main.ctl &
done
入射角を0から85度まで5度ずつ変えた18個のプログラムを同時に走らせています。
**MEEPで用いられるMaxwell方程式 [#t352b35e]
MEEPで用いられるMaxwell方程式は
++単位長さaをユーザーが設定する。
++真空の誘電率&mathjax{ε_0};が1。
++真空の透磁率&mathjax{μ_0};が1。
など通常のMaxwell方程式と異なります。
MEEPで用いられるMaxwell方程式は次のような式変形の最後の形です。
（[https://meep.readthedocs.io/en/latest/Introduction/ ドキュメント参考]）
通常の以下の形のMaxwell方程式を考えます。
#mathjax(\nabla\cdot\left(\varepsilon_0\varepsilon_r\bm{E}\right)=\rho)
#mathjax(\nabla\cdot\left(\mu_0 \bm{H}\right)= 0)
#mathjax(\nabla\times\bm{E}= - \frac{\partial}{\partial t}\left(\mu_0 \bm{H}\right))
#mathjax(\nabla\times\bm{H}=\bm{J} + \frac{\partial}{\partial t}\left(\varepsilon_0\varepsilon_r\bm{E}\right))
次のように変形します。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/maxwell_eq2.png
&mathjax{a};はシミュレーションで使用する単位長さです。
&mathjax{c};は光速です。
以下のように&mathjax{x',y',z',t',ω',E',H',σ'};を定めます。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/maxwell_eq3.png
このとき、Maxwell方程式は次のようになります。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/maxwell_eq4.png
この形のMaxwell方程式をMEEPで使っています。
このとき、真空の誘電率、透磁率はどちらも1です。
また、単位長さは&mathjax{a};[m]です。
MEEP中の単位時間&mathjax{a/c};に真空中の光は単位長さaだけ進みます。
この光の速度をSI単位系に直すと、&mathjax{a / (a/c) = c}; となり、矛盾していないことが確認できますね。
MEEPの単位系では光は単位時間&mathjax{c/a};に単位長さ&mathjax{a};だけ進むことから、周波数 &mathjax{f [c/a]}; の真空中の光の波長は &mathjax{1/f [a]}; です。
&mathjax{x',y',z',t',ω',σ'};はすべて無次元です。
金属が存在しない場合(&mathjax{σ=0};)は、&mathjax{a};の具体的な値を決めずにシミュレーションを行ったあとに、&mathjax{a};の値を決めてスケール変換することで&mathjax{a};の値に沿った正しいシミュレーション結果を得られます。
つまり、一回のシミュレーションで任意の単位長さaに対するシミュレーション結果が実質的に得られます。
一方、金属が存在する場合(σ!=0)では、MEEPのシミュレーションに用いる伝導率が単位長さaの値に依存するので、事前にaの値を定めておく必要があります。
電場E'と磁場H'のMEEP中での単位は等しくなっており、平面波の場合には電場と磁場の値が等しくなっています。
確かめてみましょう。平面波の場合、インピーダンスZ=\sqrt(μ/ε)を用いると磁場H'=\sqrt(μ)H ={\sqrt(μ)/Z}E=\sqrt(ε)E=E' となります。 H'=E'であることが確認できましたね。
角振動数ω'（MEEP単位系）で 比誘電率を &mathjax{ε=ε_{re} + i ε_{im}}; と定めたいときは
#mathjax( \epsilon_{r} = \epsilon_{re})
#mathjax( \sigma' = \omega' \frac{\epsilon_{im}}{ \epsilon_{re}})
とすればいいことも上の式から分かります。
**波源 [#n24a815a]
***TM波の作り方 [#j61c2744]
MEEPでは、電流を波源の代わりに用いています。
例えば、Ex偏光の点光源を設置するときには、代わりにJ = (Jx, 0, 0)exp(-iωt)という点双極子を計算しています。
入射角ΘのTMモードの平面波を作るにはxy平面に
J = J_0(exp(i(qsinΘ)x - ωt), 0, 0)
という電流を作ればよいです(ただし q = ω/c)。
この電流の放射する電磁場は
E_x = E_0 exp(iqsinΘx + κz)
E_z = -iqsinΘ/κ E_0 exp(iqsinΘx + κz)
H_y = iωε/κ E_0 exp(iqsinΘx + κz)
と表せて(MaierのPlasmonic Fundamentalsのページ30やShoufieさんの博論を参照)、
κ = \sqrt{(qsinΘ)^2 - (ω/c)^2}
= -i qcosΘ
κ = +i q cos Θ は波源に向かって収束する波になるので却下です。
κ = -iqcosΘを代入すると、
E_x = E_0 exp(i(qsinΘ)x + i(qcosΘ)z)
E_z = tanΘ E_0 exp(i(qsinΘ)x + i(qcosΘ)z)
H_y = - \sqrt(μ/ε) E_0 exp(i(qsinΘ)x + i(qcosΘ)z)
となります。これは、入射角ΘのTMモードになっていますね。
電場の振幅は\sqrt(1+(tanΘ)^2) = 1/cosΘ です。
波源を立ち上げた一番最初の部分が斜めに入射していない（xy平面に水平）なのが気になる場合は、点光源を波面に平行に並べても、斜め入射のTM波を作るという方法もあります。
***TE波の作り方 [#z315d617]
電流ではなく磁流を使う。
**左円偏光 [#u95702e4]
次のコードで左回りの円偏光になります。
sourcesクラスのamplitudeによって、y方向の電磁場の位相をx方向に比べてπ/2だけ進めています。
(set! sources
(list
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ex)
(center 0 0 8)
(size 50 50 0)
)
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ey)
(center 0 0 8)
(size 50 50 0)
(amplitude (exp (* 0+1i 1.570796327)))
)
)
)
**右円偏光 [#u95702e4]
次のコードで右回りの円偏光になります。
sourcesクラスのamplitudeによって、y方向の電磁場の位相をx方向に比べてπ/2だけ遅らせています。
(set! sources
(list
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ey)
(center 0 0 8)
(size 50 50 0)
)
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ex)
(center 0 0 8)
(size 50 50 0)
(amplitude (exp (* 0+1i 1.570796327)))
)
)
)
**45度直線偏光 [#u95702e4]
次のコードで45度方向((x,y)=(1,1)方向)に直線偏光した平面波になります。
x方向に直線偏光した平面波とy方向に直線偏光した平面波を重ねることで45度方向の直線偏光にしています。
(set! sources
(list
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ex)
(center 0 0 8)
(size 50 50 0)
)
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ey)
(center 0 0 8)
(size 50 50 0)
)
)
)
**斜め入射 [#e6b41572]
斜め入射のシミュレーションの方法は3つあります。
***全反射直前のTM波(ブロッホ周期境界条件) [#s8084013]
ブロッホの周期的境界条件を用いて電磁波の斜め入射をシミュレーションすることができます。
次の図のように画面に平行な方向にxz平面を取ります。
右方向にx軸、縦方向にz軸、画面垂直方向にy軸を取ります。
z<0には比誘電率2.25、z>0には比誘電率1.25の誘電体を敷き詰めます。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/config_oblique.PNG
誘電体界面に入射角48度で電磁場が入射したときのシミュレーションを考えます。
全反射が入射角arcsin(sqrt(1.25/2.25))=48.2度で起きるので入射角48度は全反射直前の角度です。
z軸に垂直な境界にはPMLを、x軸に垂直な境界にはブロッホの周期的境界条件を用います。
計算領域の大きさはx,y,z方向それぞれについて10,0,10です。
電磁場はy方向依存性を持たないので、y方向の計算領域の大きさは0となっています。
入射角をΘとすると、波数ベクトルkはk=k_0(sinΘ,0,cosΘ)です。
波源はx軸に平行な線として定義するので、位相のx依存性としてexp(i(x k_0 sinΘ))を持つ必要があります。
Schemeではsourceクラスのオプションであるamp-funcを設定することで位相のx座標依存性を設定できます。
また、計算する領域内（今回の場合は大きさ10×0×10の長方形）を単位胞として見た場合の基本並進ベクトルをRとすると、
x軸と垂直な計算領域の境界の点rについてE(r)exp(i k・R) = E(r+R)が成り立ちます。
磁場についても同様です。
点rと点r+Rはともに境界上の点になっています。E(r)は点rでの電場ベクトルです。
このような境界条件はブロッホの周期的境界条件を用いることで実現できます。
Schemeの場合はk-point変数を設定することでブロッホの周期的境界条件を設定することができます。
(set! geometry-lattice (make lattice (size 10 no-size 10)))
(set-param! resolution 10)
(define (pw-amp k x0) (lambda (x)
(exp (* 0+1i (vector3-dot k (vector3+ x x0))))))
(define-param theta-deg 48)
(define theta-rad (deg->rad theta-deg))
(define epsilon1 2.25)
(define epsilon2 1.25)
(define-param fcen 0.8) ; pulse center frequency
(define-param df 0.02) ; turn-on bandwidth
(define-param kdir (vector3 (sin theta-rad) 0 (cos theta-rad))) ; direction of k (length is irrelevant)
(define k (vector3-scale (* 2 pi fcen (sqrt epsilon1))
(unit-vector3 kdir))) ; k with correct length
(set! pml-layers (list (make pml (thickness 2)(direction Z))))
(set! k-point (vector3* (* fcen (sqrt epsilon1) ) (vector3 (sin theta-rad) 0 (cos theta-rad))))
(set! geometry
(list
(make block
(center 0 0 -2.5)
(size 10 infinity 5)
(material (make dielectric(epsilon epsilon1)))
)
(make block
(center 0 0 2.5)
(size 10 infinity 5)
(material (make dielectric(epsilon epsilon2)))
)
)
)
(set! sources
(list
(make source
(src (make continuous-src (frequency fcen) (fwidth df)))
(component Ez) (center 0 0 -3) (size 10 0 0)
(amp-func (pw-amp k (vector3 0 0 -3))))))
(define-param T 40) ; run time
(run-until T
(to-appended "ez" (at-every 0.1 output-efield-z))
(to-appended "ex" (at-every 0.1 output-efield-x))
)
Ezのプロット
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/main-ez.gif
Exのプロット
屈折角が90度近いので、透過光のExの大きさは小さくなっています。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/main-ex.gif
***光源を2つ並べる [#d178c899]
[https://github.com/stevengj/meep/blob/master/scheme/examples/pw-source.ctl コード]
入射角45度の場合に限ってのみ、2つの波原で斜め入射の平面波が作れます。
計算領域がxz平面上に乗っている場合を例として考えます。
境界条件は四方全てPMLとします。
左側の境界と下側の境界にそれぞれ線状の波原を設置します。
波原の位相のx,z依存性は波数ベクトルk=(k_0sinθ, 0, k_0cosθ)であることを考えるとそれぞれexp(ixk_0sinθ),exp(ixk_0cosθ)です。
Schemeではsourceクラスのオプションであるamp-funcを設定することで位相のx,z座標依存性を設定できます。
Ezのプロット (dfを0.1に変更しています）
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/ez.gif
光源が一つしかない場合のEzのプロット（dfを0.1に変更しています）
このプロットから一様な振幅の大きさを持つ斜め入射の電磁場を作るには光源が二つ必要なことが分かります。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/ez-single.gif
***光源を斜めにたくさん並べる [#t6b81813]
位相の揃った点光源を一直線に斜めに並べると、光源の並んだ方向と垂直な方向に進行する平面波が作れます。
直線状には点光源の強度がガウス分布になるように並べています。
光源の強度をすべて同じにすると干渉によって強度分布が不均一になります。
半値幅を調節することでビームの横幅を変更できます。
(define-param pi 3.14159265)
(define-param dpml 3.0)
(define-param len (+ 20 dpml))
(define-param wid (+ 20 dpml))
(set! geometry-lattice (make lattice (size len wid no-size)))
(define-param beam-waist 2.5) ; beam sigma (gaussian beam width)
(define-param rotation-angle (* (/ 22.5 360) 2 pi)) ; Degrees if you like them
(define-param source-points 60) ; should be a big number
(define-param source-size (* 4 beam-waist)) ; should be bigger than beam-waist
(define-param src_list (list ))
(do
((
r_0
(/ source-size -2)
(+ r_0 (/ source-size (- source-points 1)))
))
((> r_0 (/ source-size 2)))
;for r_0 = -source-size/2 : source- size/source-points : source-size/2
(set! src_list (append src_list
(list (make source
(src (make gaussian-src
(wavelength 1)
(width 3)
))
(amplitude (exp (- 0 (/ (* r_0 r_0) (* 2 beam-waist beam-waist)))))
(component Ez)
(center (* r_0 (sin rotation-angle)) (* r_0 (cos rotation-angle)))
))
))
)
(set! sources src_list)
(set! pml-layers (list (make pml
(thickness dpml)
)))
(set! resolution 10)
(use-output-directory)
(run-until (* 2 len)
(to-appended "ez" (at-every 0.1 output-efield-z))
)
plot of Ez
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/array.gif
[http://article.gmane.org/gmane.comp.science.electromagnetism.meep.general/3058 参考1]
[http://kenji-ishikawa.cocolog-nifty.com/plasma/2010/01/post-ea91.html 参考2]
**Otto geometry [#y9a02bff]
TM表面プラズモンを励起するには、入射光の波数の界面に平行な成分および周波数が、励起されるTM表面プラズモンの波長および周波数と一致する必要があります。
この励起状態を満たすためにOtto geometryという構造が使われることがあります。
import meep as mp
import math
import sys
resolution = 100
# 系の大きさ size_x,size_y,size_z
sx = 10
sy = 10
sz = 0
cell = mp.Vector3(sx, sy, 0)
# PML層の厚さ（MEEP単位）depth_pml
dpml = 1
pml_layers = [mp.PML(thickness=dpml)]
# 何秒（MEEP単位）おきにデータを取るか
dT = 0.2
sources = []
for i in range(100):
x_=-4 + 4 * i / 100.
y_=0 + 4 * i / 100.
r=math.sqrt(abs(x_+2)**2+abs(y_-2)**2)
coe=math.exp(-r*r)
sources.append(mp.Source(src=mp.ContinuousSource(wavelength=1),
center=mp.Vector3(x=x_, y=y_, z=0),
component=mp.Ex,
amplitude=coe
))
sources.append(mp.Source(src=mp.ContinuousSource(wavelength=1),
center=mp.Vector3(x=x_, y=y_, z=0),
component=mp.Ey,
amplitude=coe
))
material_=mp.Medium(
epsilon=1,
E_susceptibilities = [mp.DrudeSusceptibility(frequency=math.sqrt(3), gamma=0, sigma=1)],
)
geometry = [
mp.Block(size=mp.Vector3(x=20, y=5, z=mp.inf),
center=mp.Vector3(x=0, y=2.6, z=0),
material=mp.Medium(epsilon=4)
),
mp.Block(size=mp.Vector3(x=20, y=5, z=mp.inf),
center=mp.Vector3(x=0, y=-2.6, z=0),
material=material_
)
]
sim = mp.Simulation(cell_size=cell,
boundary_layers=pml_layers,
geometry=geometry,
sources=sources,
resolution=resolution
)
sim.run(
mp.to_appended("ey",
mp.at_every(dT, mp.in_volume(mp.Volume(center=mp.Vector3(0, 0, 0), size=mp.Vector3(x=sx, y=sy, z=0)), mp.output_efield_y))),
until=20
)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/otto_geometry2/ey.gif
**負の屈折率 [#t4f04921]
屈折率nはn=√mu√eps（mu,epsはそれぞれ比透磁率と比誘電率）と書けるので、muとepsがともに負のときは屈折率nが負になります。
このとき波の位相速度は通常の位相速度と進行方向が逆向きになります。
シミュレーションの際に、負の屈折率を持つ媒質をPMLと重ねないように注意してください。
PMLと負の屈折率を持つ媒質を重ねると、PML中で電磁波は吸収されずに増大します。
import meep as mp
import math
import sys
resolution = 20
# 系の大きさ size_x,size_y,size_z
sx = 20
sy = 15
sz = 0
cell = mp.Vector3(sx, sy, 0)
# PML層の厚さ（MEEP単位）depth_pml
dpml = 1
pml_layers = [mp.PML(thickness=dpml)]
# 何秒（MEEP単位）おきにデータを取るか
dT = 0.5
sources = []
for i in range(100):
x_=-4 + 4 * i / 100.
y_=0 + 4 * i / 100.
r=math.sqrt(abs(x_+2)**2+abs(y_-2)**2)
coe=math.exp(-r*r)
sources.append(mp.Source(src=mp.ContinuousSource(wavelength=1),
center=mp.Vector3(x=x_, y=y_, z=0),
component=mp.Ez,
amplitude=coe
))
material_=mp.Medium(
epsilon=1,
E_susceptibilities = [mp.DrudeSusceptibility(frequency=2, gamma=0.00145, sigma=1)],
mu=1,
H_susceptibilities = [mp.DrudeSusceptibility(frequency=2, gamma=0.00145, sigma=1)],
)
geometry = [
mp.Block(size=mp.Vector3(x=15, y=2, z=mp.inf),
center=mp.Vector3(x=0, y=-3, z=0),
material=material_
)
]
sim = mp.Simulation(cell_size=cell,
boundary_layers=pml_layers,
geometry=geometry,
sources=sources,
resolution=resolution
)
sim.run(
mp.to_appended("ez",
mp.at_every(dT, mp.in_volume(mp.Volume(center=mp.Vector3(0, 0, 0), size=mp.Vector3(x=sx, y=sy, z=0)), mp.output_efield_z))),
until=30
)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/negative_index/ez.gif
**反射率、吸収率、透過率を計算したい [#cee00bb6]
物体のあるときとないとき両方についてガウシアンソースを入射したときの、ポインティングベクトルを物体の前後で積分することで、反射率・透過率・吸収率が計算できます。
例えば透過率ならば透過光のポインティングベクトルの積分/入射光のポインティングベクトルの積分を求めれば良いです。
試していませんが、frequency solverを使って単一周波数のContinuous Sourceに対する応答を計算しても、反射率・透過率・吸収率が求まると思います。
ここでは一つ目の方法を紹介します。
***TM波の反射・透過率 [#e20bae0b]
誘電体にTM波が斜め入射したときの反射率・透過率を計算しましょう。
z>0の領域に誘電体（誘電率9）が存在し、入射面がxz平面の場合を考えます（z<0は真空）。
入射角Θのとき、入射光の波数ベクトルが \vec{k} = k_0(sinΘ, 0, cosΘ) で、
この波源の波数ベクトルがブロッホの周期境界条件の波数kになります。
電磁場の周期性がないz方向も周期境界条件が設定されてしまうのですが、
z方向にはPMLを設置すれば疑似的にz方向のセルの大きさが無限大になるので問題ありません。
したがって、このシミュレーションの場合にはブロッホの周期境界条件の波数ベクトルのうち、z成分の値にあまり意味はありません。
; sz : computational domain size in z direction
; theta-rad : incident angle in radian
(set-param! resolution 200)
(define-param sz 3)
(set! geometry-lattice (make lattice(size no-size no-size sz)))
(define-param no-obstacle? false)
(define-param theta-deg 0)
(define theta-rad (deg->rad theta-deg))
(define fcen 3)
(define df 0.01)
(define epsilon1 1)
(define epsilon2 9)
(define nfreq 1)
(set! dimensions (if (= theta-rad 0) 1 3))
(set! k-point (vector3* fcen (vector3 (sin theta-rad) 0 (cos theta-rad))))
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ex)
(center 0 0 -1))))
(if (not no-obstacle?)
(set! geometry
(list
(make block
(center 0 0 -0.5)
(size infinity infinity 1)
(material (make medium(index (sqrt epsilon1)))))
(make block
(center 0 0 0.5)
(size infinity infinity 1)
(material (make medium(index (sqrt epsilon2))))))))
(set! pml-layers (list (make pml (thickness 1))))
(define trans
(add-flux fcen df nfreq
(make flux-region
(center 0 0 0.5)))
)
(define refl
(add-flux fcen df nfreq
(make flux-region
(center 0 0 -0.5))))
(run-sources+ 50
(stop-when-fields-decayed 50 Ex (vector3 0 0 0) 1e-9))
(display-fluxes trans refl)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/TM/out.png
***TE波の反射・透過率 [#lde2976a]
; sz : computational domain size in z direction
; theta-rad : incident angle in radian
; fcen : pulse center frequency (500 nm)
; df : pulse width in frequenhcy
; nfreq : number of frequencies at which to compute flux
; unit length: 100 micro meter
(set-param! resolution 200)
(define-param sz 3)
(set! geometry-lattice (make lattice(size no-size no-size sz)))
(define-param no-obstacle? false)
(define-param theta-deg 0)
(define theta-rad (deg->rad theta-deg))
(define fcen 3)
(define df 0.01)
(define epsilonA 1)
(define epsilonB 9)
(define nfreq 1)
(set! k-point (vector3* fcen (vector3 (sin theta-rad) 0 (cos theta-rad))))
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ey)
(center 0 0 -1)
)))
(if (not no-obstacle?)
(set! geometry
(list
(make block
(center 0 0 -0.5)
(size infinity infinity 1)
(material (make medium(index (sqrt epsilonA))))
)
(make block
(center 0 0 0.5)
(size infinity infinity 1)
(material (make medium(index (sqrt epsilonB))))
)
)
)
)
(set! pml-layers (list (make pml (thickness 1))))
(define trans
(add-flux fcen df nfreq
(make flux-region
(center 0 0 0.5)
)
)
)
(define refl
(add-flux fcen df nfreq
(make flux-region
(center 0 0 -0.5)
)
)
)
(run-sources+ 50
(stop-when-fields-decayed 50 Ey (vector3 0 0 0) 1e-9)
;(at-beginning output-epsilon)
;(to-appended "ex" (at-every 1 output-efield-x))
)
(display-fluxes trans refl)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/TE/out.png
***誘電体とグラフェンの多層構造の反射・透過・吸収率 [#n9c84eaa]
アンドープのグラフェンを誘電体で挟むと特殊な条件で吸収率が2%から50%に上昇することが以下の論文で指摘されています。
(S. A. Nulli, M. S. Ukhtary, R. Saito, Appl. Phys. Lett. 112, 073101-1-4, (2018))。
(J. Blumberg, M. S. Ukhtary, R. Saito, Phys. Rev. Applied 10, 064015-1-14, (2018))。
MEEPで再現できます。
グラフェンを二種類の誘電体A, Bで挟んで例えばABABGBABAのような並びにします。
ここでGはグラフェンの層です。
それぞれの誘電体の厚みは垂直入射したときの光路長が1/4波長になるようにします。
誘電体Aの層の数を2nとします。
このとき、誘電体AとBの誘電率の比がε_A/ε_B={2/(Zσ)}^{1/n}のときグラフェン近傍の電場が増強されて光の吸収率が50%に達します。
ここで、Zは真空のインピーダンス、σはグラフェンの電気伝導率です。
シミュレーションではグラフェンの伝導率σはσ=πe^2/(2h)としています(eは電荷素量、hはプランク定数)。
平面波が誘電体とグラフェンの層に対して垂直入射するとき電磁場の値は層と平行な方向に対して依存性を持ちません。
従ってシミュレーションは一次元で行うことができます(マクスウェル方程式に含まれる微分のうち、層と平行な方向の微分を０に置換するとマクスウェル方程式に含まれる位置の変数が一つになります)。
; a : unit length in the program
; e : elementary charge
; hbar : reduced planck constant
; epsilon_vaccum : permittivity of vaccum
; sigma : the conductivity of graphene in the defenition of 2D material
; sigma_d : the conductivity of graphene in the defenition of MEEP
; d : the thickness of graphene in the SI unit
; dg : the thickness of graphene in this program unit
; fcen : pulse center frequency (500 nm)
; df : pulse width in frequency
; nfreq : number of frequencies at which to compute flux
; isEmpty : Whether graphene exist or not.If false, we calculate in vacuum
; sz : the size of computatinal domain in Z direction in this program unit
; Z_0 : the impedance of the vaccum
; s : the number of repetitions of the AB layers before mirroring.
; alpha : When the ratio of index of dielectric A and B is alpha, the absorbance is up to 50%
; da : thickness of dielectric A
; db : thickness of dielectric B
(define a 1e-9)
(define e (* 1.60217662 1e-19))
(define c 299792458)
(define hbar (* 1.05457180013 1e-34))
(define epsilon_vacuum (* 8.854187817 1e-12))
(define pi (* 4 (atan 1)))
(define d (* 0.3 1e-9))
(define dg 0.3)
(define sigma (/ (* e e) (* 4 hbar)))
(define sigma_d (/ (* sigma a) (* c epsilon_vacuum d)))
(define fcen 0.002)
(define df 0.0001)
(define nfreq 100)
(define s 2)
(define Z_0 376.730313461)
(define alpha (expt
(/ 2 (* Z_0 sigma))
(/ (* 2 s))))
(define epsilon_a (* alpha alpha))
(define epsilon_b (* 1 1))
(define da (/ (* 4 fcen (sqrt epsilon_a) )))
(define db (/ (* 4 fcen (sqrt epsilon_b) )))
(set-param! resolution 10)
(define sz (* 2 400))
(set! geometry-lattice (make lattice(size no-size no-size sz)))
(define-param isEmpty "DO NOT SET VALUE DIRECTLY")
(set! dimensions 1)
(if (string? isEmpty) ((print "Set the parameter \"isEmpty\" true or false from command line.\nDo not use the parameter isEmpty as it is.\n")(exit)))
(if (not isEmpty)
(set! geometry
(list
(make block
(center 0 0 0)
(size infinity infinity dg)
(material (make medium (epsilon 1)(D-conductivity sigma_d)))
)
(make block
(center 0 0 (+ (/ dg 2) (/ db 2)))
(size infinity infinity db)
(material (make dielectric(epsilon epsilon_b)))
)
(make block
(center 0 0 (+ (/ dg 2) db (/ da 2)))
(size infinity infinity da)
(material (make dielectric(epsilon epsilon_a)))
)
(make block
(center 0 0 (+ (/ dg 2) db da (/ db 2)))
(size infinity infinity db)
(material (make dielectric(epsilon epsilon_b)))
)
(make block
(center 0 0 (+ (/ dg 2) db da db (/ da 2)))
(size infinity infinity da)
(material (make dielectric(epsilon epsilon_a)))
)
(make block
(center 0 0 (- (+ (/ dg 2) (/ db 2))) )
(size infinity infinity db)
(material (make dielectric(epsilon epsilon_b)))
)
(make block
(center 0 0 (- (+ (/ dg 2) db (/ da 2))))
(size infinity infinity da)
(material (make dielectric(epsilon epsilon_a)))
)
(make block
(center 0 0 (- (+ (/ dg 2) db da (/ db 2))))
(size infinity infinity db)
(material (make dielectric(epsilon epsilon_b)))
)
(make block
(center 0 0 (- (+ (/ dg 2) db da db (/ da 2))))
(size infinity infinity da)
(material (make dielectric(epsilon epsilon_a)))
)
)
)
)
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ex)
(center 0 0 350)
)))
(set! pml-layers (list (make pml (thickness 1))))
(define trans
(add-flux fcen df nfreq
(make flux-region
(center 0 0 -330)
)
)
)
(define refl
(add-flux fcen df nfreq
(make flux-region
(center 0 0 330)
)
)
)
(run-sources+ 150000
(stop-when-fields-decayed 50 Ex (vector3 0 0 -330)
1e-3
)
(at-beginning output-epsilon)
(to-appended "ex" (at-every 10 output-efield-x))
)
(display-fluxes trans refl)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/dielectric_graphene/dielectric_graphene_n=2/out.png
電場分布は次のようになります（パルス幅や誘電率、入射光の波長を上のコードから変更しています。nは屈折率です）。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/dielectric_graphene/dielectric_graphene_n=2/multilayer_efield.png
**電場ベクトルをプロットしたい。 [#fee3ee23]
全反射における電場ベクトルの分布
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/evanescent/evanescent_afterEfield.gif
シミュレーションのためのコード
(define-param dx 1); cell size in x axis
(define-param dy 2); cell size along y axis
(define-param fp 1.5)
(define-param tpml 0.1); PML thickness
(define-param deg 48)
(define-param theta (* pi 0.88889)); incident angle of the plane wave
(set! theta (* deg (/ pi 180)))
(set! resolution 40);define the resolution for simulation domain
;define the lattice
(set! geometry-lattice (make lattice (size (+ dx (* 2 tpml)) (+ dy (* 2 tpml)) no-size)))
;define PML thickness
(set! pml-layers (list (make pml (thickness tpml))))
;define the blocks
(set! geometry (list
(make block (center 0 (* 0.25 dy))(size infinity (+ tpml (* 0.5 dy)) infinity)
(material (make dielectric (epsilon 1))))
(make block (center 0 (* -0.25 dy))(size infinity (+ tpml (* 0.5 dy)) infinity)
(material (make dielectric (epsilon 2.25))))
))
(set! pml-layers (list (make pml (thickness tpml))))
;(+ tpml (* 0.5 dy))
;define the wave vector
(define kx (* fp (sqrt 2.25) (sin theta)))
;give the amplitude function
(define (f_amp p) (exp (* 0+2i pi kx (vector3-x p))))
(set! k-point (vector3 kx 0 0))
;define the simulation domain to be complex field
(set! force-complex-fields? true)
;define the Gaussian beam
(set! sources
(list
(make source
(src (make continuous-src (frequency fp)))
(component Ex)
(component Ey)
(center 0 (* dy -0.50))
(size (+ dx (* 2 tpml)) 0 )
(amp-func f_amp))))
(set! pml-layers (list (make pml (thickness tpml) (direction Y))))
;extract the data to .png file
(run-until 40
(at-beginning output-epsilon)
;(at-end output-efield-y)
(to-appended "ey"(at-every 0.05 (in-volume (volume (center 0 0 0)(size dx dy 0))output-efield-y)))
(to-appended "ex"(at-every 0.05 (in-volume (volume (center 0 0 0)(size dx dy 0))output-efield-x)))
)
電場ベクトルをプロットするには、HDF5ファイルから直接プロットするコードを自分で書く必要があります。
pythonのmatplotlibのquiver関数を用いるといいと思います。
**電荷分布を見たい。 [#l0ec2ad8]
MEEPには電荷分布を計算する機能がないため、自作する必要があります。
ガウスの法則を差分化すればよいです。
前述の全反射直前のTM波(ブロッホ周期境界条件)のシミュレーションについて例を示します。
電荷分布(+電場ベクトル)のプロット
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/charge.gif
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/charge.py コード]
誘電体界面で誘電分極が起きていることと波源で電荷の移動があることが分かります。
**PMLをある方向にだけ設置する（Python） [#uf2f674c]
PMLをある方向にだけ設置するにはオプションdirectionを設定すればよいです。
例えば、Y方向だけに設置するには次のようにすればよいです。
pml_layers = [mp.PML(thickness=1.0, direction=mp.Y)]
PMLが設置されていない方向は固定端の全反射が起きます。
真ん中にwaveguideを置いて、Y方向にだけPMLを設置したときの電場は次のようになります。
X方向ではPMLがないため全反射が起きて定在波が生じています。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/pml-ydirection.gif
コードは次の通りです。
import meep as mp
cell = mp.Vector3(16,8,0)
geometry = [mp.Block(mp.Vector3(1e20,1,1e20),
center=mp.Vector3(0,0),
material=mp.Medium(epsilon=9))]
sources = [mp.Source(mp.ContinuousSource(frequency=0.1),
component=mp.Ez,
center=mp.Vector3(-7,0))]
pml_layers = [mp.PML(thickness=1.0, direction=mp.Y)]
resolution = 10
sim = mp.Simulation(cell_size=cell,
boundary_layers=pml_layers,
geometry=geometry,
sources=sources,
resolution=resolution)
dT = 0.2
sim.run(
mp.at_beginning(mp.output_epsilon),
mp.to_appended("ez",
mp.at_every(dT, mp.output_efield_z)),
until=100
)
**S字や螺旋の作り方 [#u5436192]
MEEPでは後に置いた物体が先に置いた物体を完全に上書きするため、S字や螺旋を小数のパーツで作ることは難しいです。
しかし、たくさんのディスクを配置することでS字や螺旋を形成することができます。
***アルキメデスの螺旋 [#fd807100]
import meep as mp
import math
resolution = 10
# 系の大きさ size_x,size_y,size_z
sx = 20
sy = 20
sz = 20
cell = mp.Vector3(sx, sy, sz)
# PML層の厚さ（MEEP単位）depth_pml
dpml = 3
pml_layers = [mp.PML(dpml)]
# 何秒（MEEP単位）おきにデータを取るか
dT = 0.1
sources = []
geometry_list=[]
for i in range(1000):
theta=i/1000.*2*math.pi*3
radius_=i/200.
geometry_list.append(mp.Cylinder(radius=0.5,height=0.35,center=mp.Vector3(x=radius_*math.cos(theta),y=radius_*math.sin(theta),z=0),material=mp.Medium(epsilon=12)))
geometry=geometry_list
sim = mp.Simulation(cell_size=cell,
boundary_layers=pml_layers,
geometry=geometry,
sources=sources,
resolution=resolution)
sim.run(
mp.at_beginning(mp.in_volume(mp.Volume(center=mp.Vector3(0, 0, 0), size=mp.Vector3(x=sx - 2 * dpml, y=sy - 2 * dpml, z=0)), mp.output_epsilon)),
until=1
)
z=0平面における誘電率のプロットは次の通りです（PMLは除いた範囲をプロットしています）。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/archimedes_spiral/after-S-drude-Au-eps-000000.00.png
**物性値のライブラリ [#ue28f25c]
***Schemeの場合 [#z7a9ca6a]
物性値のライブラリが公開されています。
[https://github.com/NanoComp/meep/blob/master/scheme/materials.scm materials-library.scm]をダウンロードしてください。
ライブラリを使用するには、プログラムの最初に
(include "/path/to/materials-library.scm")
を追加してライブラリをロードしてください。
"/path/to/materials-library.scm"はmaterials-library.scmを設置した場所のアドレスに書き換えてください。
丸岡がダウンロードしたmaterials.scmを使用する場合は次のパスを使用してください。
(include "/home/students/maru/meep/scheme/materials.scm")
ただし最新のmaterials.scmではない可能性があることに注意してください。
materials.scmはtube61, 62, 63でしか動きません（flex, tube60では動きません）。
詳しい使い方は https://meep.readthedocs.io/en/latest/Materials/ を参照してください。
materials.scmはMEEPの単位長さが1μmであるものとして設計されています。
異なる単位長さを使用したい場合、特別に指定する必要があります。
例えば、0.1µmを単位長さにしたいときは
(set! um-scale (* um-scale 0.1))
としてください。
***Pythonの場合 [#d7e5fefa]
物性値のライブラリmaterials.pyをロードするのに特別な操作は必要ありません。
materials.pyはMEEPの単位長さが1μmであるものとして設計されています。
異なる単位長さを使用したい場合、materials.pyを直接書き換える必要があります。
materials.pyのum_scaleが宣言されている行以降でum_scaleを目的の単位長さに書き換えてください。
例えば、0.1µmを単位長さにしたいときは
um_scale=um_scale*0.1
としてください。
materials.pyの場所は次のようにして調べてください。
(pmp) maru@tube60:~$ python
Python 3.6.6 | packaged by conda-forge | (default, Jul 26 2018, 09:53:17)
[GCC 4.8.2 20140120 (Red Hat 4.8.2-15)] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import meep
Using MPI version 3.1, 1 processes
>>> meep.__file__
'/home/students/maru/anaconda3/envs/pmp/lib/python3.6/site-packages/meep/__init__.py'
この場合、materials.pyのアドレスは/home/students/maru/anaconda3/envs/pmp/lib/python3.6/site-packages/meep/materials.pyです。
このmaterials.pyの書き換えは(pmp)の仮想環境でだけ有効であることに注意してください。
((mp)の仮想環境でmaterials.pyの単位長さを変える場合は(mp)の仮想環境のmaterials.pyを書き換えてください）
**ユーザー定義のドルーデ伝導 [#df7ff749]
ドルーデ伝導で表される次のような比誘電率を用いる方法を説明します。
（&mathjax{\varepsilon_\infty};は周波数が十分高いところでの比誘電率（定数）。&mathjax{f_p};はプラズマ周波数）
#mathjax(\varepsilon_r = \varepsilon_{\infty} -\frac{{f_p}^2}{f(f+iγ)})
MEEPのドルーデ伝導はMEEPの単位系の&mathjax{f_p'};, &mathjax{γ'};, &mathjax{f'};で表されているので、SI単位系の&mathjax{f_p};, &mathjax{γ};, &mathjax{f};との関係を求める必要があります。
#mathjax(\begin{align} \varepsilon_r &= \varepsilon_{\infty} -\frac{{f'_p}^2}{f'(f'+iγ')} \\ &=\varepsilon_{\infty} - \frac{{f'_p}^2}{(\frac{a}{c}f)(\frac{a}{c}f+iγ')} \\ &=\varepsilon_{\infty} - \frac{{f'_p}^2}{(\frac{1}{λ'})(\frac{1}{λ'}+iγ')} \\ &= \varepsilon_{\infty} - \frac{{(\frac{c}{a}f'_p)}^2}{f(f+i(\frac{c}{a}γ'))} \end{align})
この式から、SI単位系の&mathjax{f_p};, &mathjax{γ}; とMEEPの単位系の &mathjax{f_p'};, &mathjax{γ'};は次の関係にあります。
#mathjax(f_p'=\frac{a}{c}f_p)
#mathjax(γ'=\frac{a}{c}γ)
上の関係式から求めた&mathjax{f_p'};, &mathjax{γ'};と&mathjax{\varepsilon_\infty};について、
user_defined_material = mp.Medium(epsilon=ε_∞, E_susceptibilities=[mp.DrudeSusceptibility(frequency=f_p', gamma=γ', sigma=1)])
とすると、その&mathjax{f_p'};, &mathjax{γ'};と&mathjax{\varepsilon_\infty};で表されるドルーデ伝導を持つ物質が定義できます。
この物質を使用する際は、geometryのオプション materials = user_defined_material を設定する必要があります。
#mathjax(\begin{align} \varepsilon_r &= \varepsilon_\infty + i \frac{\sigma}{\omega\varepsilon_0} \\ &= \varepsilon_\infty + i \frac{\sigma_0}{\omega\varepsilon_0ω(1-iτω)} \\ &= \varepsilon_\infty - \frac{\sigma_0}{(2\pi)^2\varepsilon_0τf(\frac{i}{2\piτ}+f)} \end{align})
&mathjax{\sigma_0=f_p(2\pi)^2\varepsilon_0\tau};を代入すると、
#mathjax(\varepsilon_r = \varepsilon_\infty - \frac{f_p^2}{f(\frac{i}{2\piτ}+f)})
散乱時間τとγの関係は
#mathjax(γ=\frac{1}{2\piτ})
となる。よって散逸がない場合&mathjax{\gamma=0};とする。
伝導率と&mathjax{f_p'};, &mathjax{γ'};と&mathjax{\varepsilon_\infty};の関係は
#mathjax(\begin{align}\sigma &=-\left(\frac{i}{2\pi f\varepsilon_0}\right)^{-1} \frac{f_p^2}{f(iγ+f)} \\ &=\frac{2\pi\varepsilon_0f_p^2}{γ-if} \\ &=i\frac{2\pi\varepsilon_0 f_p^2}{f+γi} \\ &=i\frac{2\pi\varepsilon_0 {f'}_p^2 a^2}{f'ca+c^2γi}\end{align})
**HDF5の取り扱い [#ic8d7598]
MEEPのシミュレーション結果はHDF5ファイルの形で出力されます。
PythonではHDF5ファイルを編集するためのh5pyというライブラリが提供されています。
以下のコマンドでインストールできます。
apt-get install python-h5py
bash on windowsでは"pip install h5py"でもインストールできますが、次のエラーによりインポートできませんでした。
ImportError: No module named h5py
import h5py #以降のプログラムのpukiwikiでの掲載ではこの行は省略する
f = h5py.File("file-name.h5", "r") #rとすると読み込み専用で開く
print(f)
とすると、<HDF5 file "file-name" (mode r)>と表示され、読み込みができていることが確認できます。
print(list(f).keys())
とすると、 h5fileのkeyが表示されます。
例えば、['ex']と表示されたとします。
Ex = f['ex'][:]
とすると、キー'ex'に紐づいた配列をExに取り出せます。
このとき、numpyの配列に変換されてExに配列が渡されています。
numpyの配列に変換せずにf['ex'][i]のようにして、何回もアクセスすると処理に時間がかかるので注意が必要です。
print(Ex)
とすると、Exの配列の中身が確認できます。
要素が多すぎるときは・・・と表示されて省略されます。
HDF5ファイルの書き込みは[http://docs.h5py.org/en/stable/quick.html#appendix-creating-a-file こちら]。
**電磁場のh5ファイルをmp4に変換するためのスクリプト [#x4773146]
以下のスクリプトを.bash_profileに貼り付けて使用すると、h5ファイルをmp4ファイルに自動で変換できて便利です。
function h5mp4(){
name=${1:0:-3} # remove ".h5"
echo $name
mkdir $name
cp $name.h5 $name
cd $name
h5topng $name.h5 -Z -c dkbluered -R -t 0:$2
convert $name*png $name.gif
ffmpeg -i $name.gif -qscale:v 1 $name.mp4
}
例えば時間の区切り回数が1000回あるようなhogehoge.h5を動画にしたいときは
h5mp4 hogehoge.h5 999
とするとhogehoge.gifとhogehoge.mp4とhogehoge-<数字>.pngが自動で生成されます。
時間の区切り回数が知りたいときは
h5ls hogehoge.h5
としてください。
ez Dataset{100, 200, 1000/Inf}
と表示される数字のうち、1000が区切り回数です。
**印加電場と実際の電場の時間変化を折れ線グラフでプロット [#q9e259c0]
import h5py
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from PIL import Image
import moviepy.editor as edit
from mpl_toolkits.axes_grid1 import make_axes_locatable
plt.rcParams["font.size"] = 18
#----------input----------------
surfaceExFilename ="After-ex_+0.5.h5"
surfaceEyFilename = "After-ey_+0.5.h5"
surfacePxFilename = "Pre-ex_+0.5.h5"
surfacePyFilename = "Pre-ey_+0.5.h5"
savename = "current-line-withPreE-takingMiddleCurrent.png"
startTime = 0
endTime = 1000
waveLength = 6
timeResolution = 10
xid=16 #array[xid][yid][:]の電場をプロット
yid=46
#--------------------------------
def plotf(arr,label_):
time=np.array([])
amplitude=np.array([])
for i in range(startTime, endTime):
x=i * 1./(waveLength * timeResolution)
y=arr[xid][yid][i]
time=np.append(time,[x])
amplitude=np.append(amplitude,[y])
plt.plot(time,amplitude,label=label_)
plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left")
h5file_A_x = h5py.File(surfaceExFilename, "r")
h5file_A_y = h5py.File(surfaceEyFilename, "r")
h5file_P_x = h5py.File(surfacePxFilename, "r")
h5file_P_y = h5py.File(surfacePyFilename, "r")
Ax = h5file_A_x['ex'][:]
Ay = h5file_A_y['ey'][:]
Px = h5file_P_x['ex'][:]
Py = h5file_P_y['ey'][:]
**遠方界での指向性がみたい [#caf6cf04]
MEEPでは近傍界のフラックスデータから遠方界のフラックスデータを構成することができます。~
近傍界を基に遠方界のフラックスを得る利点は小さい系で計算させるためシミュレーションが速く終わることです。~
MEEPに遠方界を計算させるための条件
+観測したい方向に伝播する電磁波すべてを捕捉できるようなフラックス面を設置
+sourceはパルス波源 or [[continuous wave via frequency domain solver>https://meep.readthedocs.io/en/latest/Python_User_Interface/#frequency-domain-solver]]（フラックス面で周波数領域についてフーリエ変換を実行するため）
+遠方界の領域を円の縁で指定するなら円の半径は波長に比べて十分長くとる
+遠方界の領域を正方形の縁で指定するなら正方形の1辺の長さの半分が波長に比べて十分長くとる
アンテナは計算領域の中心に配置してください。
フラックス面はPMLより内側に設置してください。
パルス波源を使う場合は、電磁波がPMLの外に完全に発散するまでシミュレーションを続ける必要があることに注意してください。
以下に点光源を疑似的ダイポールアンテナと見立てた場合の遠方界のプロットを目的としたシミュレーションのコードを示します。~
[[該当ドキュメント>https://meep.readthedocs.io/en/latest/Python_Tutorials/Near_to_Far_Field_Spectra/]]~
[[シミュレーションのコード>http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/antenna-radiation.php]]
このシミュレーションでは円で遠方界を計算します。
計算領域はxy平面の二次元です。
Near2FarRegion()はFluxRegion()と全く同じ（関数名が違うだけ）
nearfield_boxに遠方場に必要なデータが入っているので，
get_farfieldで引数にnearfield_boxを引数にとり，遠方界の座標をVector3()で与えます。~
[[プロット用のコード>http://flex.phys.tohoku.ac.jp/~maeda/code/stable-antenna-radiation-plot.py]]~
実行コマンド群（上からシミュレーション、データの整形、プロット）
$ python antenna-radiation.py | tee source_Jz_farfields.out
$ grep farfield: source_Jz_farfields.out | cut -d , -f2- > source_Jz_farfields.dat
$ python stable-antenna-radiation-plot.py
実行からプロットまでを1つのプログラムで実現したものが
[[こちら>http://flex.phys.tohoku.ac.jp/~maeda/php/antenna-radiation.php]]
出力結果~
この結果は"antenna-radiation.py"のsrc_cmptにmp.Ex, mp.Ey, mp.Ezのどれかを入れれば
上から順にこれらの画像を作成することができます。~
それぞれxy平面をz軸方向から見ている図に相当しています。
ダイポールの振動方向と垂直な方向への放射が最も強いことが分かります。
放射の強さの最大値が1となるように規格化しています。
http://flex.phys.tohoku.ac.jp/~maeda/fig/Jx.png~
http://flex.phys.tohoku.ac.jp/~maeda/fig/Jy.png~
http://flex.phys.tohoku.ac.jp/~maeda/fig/Jz.png~
**バンド構造(分散関係)の計算 ( [https://mpb.readthedocs.io/en/latest/Python_Tutorial/ 参考] ) [#u9114db4]
まずはmpb(MIT Photonic Bands)を次のようにimportしてください。
from meep import mpb
それぞれのk点でいくつのバンドを求めるかをnum_bandsで指定してください。デフォルトは1です。
例えば8つのバンドを求めたいときは
num_bands = 8
としてください。
計算したいk点をk_pointsリストに入れてください。例えば(0,0,0)が計算したいときは次のようにしてください。
k_points=[mp.Vector3(0,0,0)]
k点の単位は&mathjax{2\pi/a};であることに注意してください(aは単位長さ)。
harminvを用いる場合のkも同様です。従って、SI単位系に直すには&mathjax{2\pi/a};を掛ける必要がある。
単位胞の大きさはgeometry_latticeによって次のように指定してください。
geometry_lattice = mp.Lattice(size=mp.Vector3(1, 1))
resolutionは2のべき乗だと計算が速くなります。
経験則として、resolution>=8でないと十分な精度が得られません。
バンド構造を計算するには次のようにしてください。
ms = mpb.ModeSolver(num_bands=num_bands,
k_points=k_points,
geometry=geometry,
geometry_lattice=geometry_lattice,
resolution=resolution)
ms.run()
TEモードとTMモードを分けたいときは、ms.run()の代わりにそれぞれms.run_te(), ms.run_tm()を使用してください。このとき、波はin-plane(k_z=0)です。
計算結果はnumpyの二次元配列 ms.all_freqs に格納されます。
i行目にはi番目(0-origin)のk点での計算結果が格納されます。
格納されている情報は次のようになっています。
ms.all_freqs[i] = [1つ目のバンドの周波数, 2つ目のバンドの周波数, 3つ目のバンドの周波数, ... ]
***誘電率のマッピング [#l15f56f7]
誘電率の画像化は次のようにしてください。
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(ms.get_epsilon(),cmap="binary")
ax.axis("off")
fig.savefig("epsilon.png",transparent=True)
単位胞を立方体(二次元の時は正方形)に取っていないと画像は歪んでしまいます。これは、誘電率データとプロットプログラムの単位ベクトルの取り方が異なるためです。
そのような場合には次のように変換した誘電率データ(converted_eps)をプロットしてください。
md = mpb.MPBData(rectify=True, periods=3, resolution=32)
eps = ms.get_epsilon()
converted_eps = md.convert(eps)
***六方(三角)格子の設定 [#ibf8e399]
[https://mpb.readthedocs.io/en/latest/Python_Data_Analysis_Tutorial/ 公式docs]のように、単位胞を
geometry_lattice = mp.Lattice(size=mp.Vector3(1, 1),
basis1=mp.Vector3(math.sqrt(3)/2, 0.5),
basis2=mp.Vector3(math.sqrt(3)/2, -0.5))
と設定してください。
basis1, basis2 が実空間の基本ベクトル。
k点を設定するときは単位が基本逆格子ベクトルになっていることに注意する（単に2π/aではない）。
***誘電体多層構造のバンド構造 [#sdf7ab3e]
[http://ab-initio.mit.edu/book/photonic-crystals-book.pdf Photonic Crystals: Molding the Flow of Light (second edition) ]の44ページに掲載されている構造をシミュレーションします。
異なる誘電体の2種類の層を交互に積層させた多層構造のバンド構造を考えます。
このような構造はレーリー卿によって1887年に最初に考えられました。
層に垂直に伝搬する（z方向とする）光のみを考えます。
GaAs（ヒ化ガリウム、誘電率=13）のバルクのバンド構造は次のようになります。
これはk=(ω/c)*sqrt(ε)のlight lineに相当します。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/multilayer_film/GaAs/dispersion.png
以下のコードをバンド構造の計算に用いました。
import meep as mp
from meep import mpb
from mpl_toolkits.mplot3d import Axes3D
import math
import cmath
import moviepy.editor as edit
from PIL import Image
import matplotlib.animation as animation
import sys
import matplotlib.pyplot as plt
import time
import h5py
import numpy as np
import matplotlib
import copy
import collections
matplotlib.use('Agg')
plt.rcParams["xtick.labelsize"] = 22 # デフォルトは12。tickのサイズ。
plt.rcParams["ytick.labelsize"] = 22 # デフォルトは12
plt.rcParams["axes.labelsize"] = 28 # デフォルトは12。軸のラベルのサイズ。
plt.rcParams['legend.fontsize'] = 20
plt.rcParams["figure.autolayout"] = True # ラベルが切れないように自動調整
plt.rcParams["axes.titlesize"] = 25 # 時刻を表示する以外には使用していない
plt.rcParams['font.family'] = 'sans-serif' # 'serif'がひげ飾りを意味するようで、いわゆる明朝系を表して、'sans-serif'はひげ飾りがないということでゴシック系を表しています。論文では、文章は明朝系、図はゴシック系という場合が多いのでここではゴシック系にしています。
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams["legend.fancybox"] = False # fancyboxは角の丸み
plt.rcParams["legend.edgecolor"] = 'black' # 枠線は黒色
plt.rcParams['xtick.major.pad'] = 8
plt.rcParams['ytick.major.pad'] = 8
num_bands = 14
resolution = 32
geometry_lattice = mp.Lattice(size=mp.Vector3(x=0, y=0, z=1))
geometry = [mp.Block(size=mp.Vector3(x=mp.inf, y=mp.inf, z=mp.inf),center=mp.Vector3(x=0,y=0,z=0), material=mp.Medium(epsilon=13))]
k_points = [
mp.Vector3(x=0,y=0,z=-0.5),
mp.Vector3(x=0,y=0,z=0),
mp.Vector3(x=0,y=0,z=0.5)
]
k_points = mp.interpolate(8, k_points)
ms = mpb.ModeSolver(
geometry=geometry,
geometry_lattice=geometry_lattice,
k_points=k_points,
resolution=resolution,
num_bands=num_bands
)
ms.run_te()
te_freqs = ms.all_freqs
te_gaps = ms.gap_list
fig=plt.figure()
ax=fig.add_subplot(111)
for i in range(len(te_gaps)):
if(te_gaps[i][0]<1):
continue
ax.axhspan(te_gaps[i][1],te_gaps[i][2],color="red",alpha=0.2)
ax.plot(te_freqs,color="red")
for i in range(len(te_freqs)):
for freq in te_freqs[i]:
ax.scatter(i,freq,color="red")
ax.set_xticks([0,9,18])
ax.set_xticklabels(["-0.5","0","0.5"])
ax.set_xlim([0,18])
ax.set_ylim([0,0.9])
ax.set_ylabel("frequency [c/a]")
ax.set_xlabel("wave vector [2π/a]")
fig.savefig("dispersion",transparent=True)
次にGaAs（ヒ化ガリウム、誘電率=13）と空気の層(誘電率=1)を交互に配置した多層構造を考えます。
単位胞の1辺の長さをaとしたとき、それぞれの層の厚さは0.5*aです。
バンド構造は次のようになります。
誘電率の異なる2種類の層を交互に並べるとバンドギャップがブリルアンゾーンの端（k=±π/a)できます。
ブリルアンゾーンの端は波長2a（2格子分）に対応し、単位胞の中心対称性から波の節がそれぞれの多層構造の真ん中にある2つの配置が存在します。
一般に低周波数のモードは高誘電率に局在したモード、高周波数のモードは低誘電率に局在したモードであるので、ブリルアンゾーンの端ではバンドギャップが開きます。
（詳しくは[http://ab-initio.mit.edu/book/photonic-crystals-book.pdf Photonic Crystals: Molding the Flow of Light (second edition) ]の47ページ)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/multilayer_film/GaAs_air/dispersion.png
バンド構造の計算は次のコードを用いました。
import meep as mp
from meep import mpb
from mpl_toolkits.mplot3d import Axes3D
import math
import cmath
import moviepy.editor as edit
from PIL import Image
import matplotlib.animation as animation
import sys
import matplotlib.pyplot as plt
import time
import h5py
import numpy as np
import matplotlib
import copy
import collections
matplotlib.use('Agg')
plt.rcParams["xtick.labelsize"] = 22 # デフォルトは12。tickのサイズ。
plt.rcParams["ytick.labelsize"] = 22 # デフォルトは12
plt.rcParams["axes.labelsize"] = 28 # デフォルトは12。軸のラベルのサイズ。
plt.rcParams['legend.fontsize'] = 20
plt.rcParams["figure.autolayout"] = True # ラベルが切れないように自動調整
plt.rcParams["axes.titlesize"] = 25 # 時刻を表示する以外には使用していない
plt.rcParams['font.family'] = 'sans-serif' # 'serif'がひげ飾りを意味するようで、いわゆる明朝系を表して、'sans-serif'はひげ飾りがないということでゴシック系を表しています。論文では、文章は明朝系、図はゴシック系という場合が多いのでここではゴシック系にしています。
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams["legend.fancybox"] = False # fancyboxは角の丸み
plt.rcParams["legend.edgecolor"] = 'black' # 枠線は黒色
plt.rcParams['xtick.major.pad'] = 8
plt.rcParams['ytick.major.pad'] = 8
num_bands = 8
resolution = 64
geometry_lattice = mp.Lattice(size=mp.Vector3(x=0, y=0, z=1))
geometry = [mp.Block(size=mp.Vector3(x=mp.inf, y=mp.inf, z=0.5), material=mp.Medium(epsilon=13))]
k_points = [
mp.Vector3(x=0,y=0,z=-0.5),
mp.Vector3(x=0,y=0,z=0),
mp.Vector3(x=0,y=0,z=0.5)
]
k_points = mp.interpolate(8, k_points)
ms = mpb.ModeSolver(
geometry=geometry,
geometry_lattice=geometry_lattice,
k_points=k_points,
resolution=resolution,
num_bands=num_bands
)
ms.run_te()
te_freqs = ms.all_freqs
te_gaps = ms.gap_list
fig=plt.figure()
ax=fig.add_subplot(111)
for i in range(len(te_gaps)):
if(te_gaps[i][0]<1):
continue
ax.axhspan(te_gaps[i][1],te_gaps[i][2],color="red",alpha=0.2)
ax.plot(te_freqs,color="red")
for i in range(len(te_freqs)):
for freq in te_freqs[i]:
ax.scatter(i,freq,color="red")
ax.set_xticks([0,9,18])
ax.set_xticklabels(["-0.5","0","0.5"])
ax.set_xlim([0,18])
ax.set_ylim([0,0.9])
ax.set_ylabel("frequency [c/a]")
ax.set_xlabel("wave vector [2π/a]")
fig.savefig("dispersion",transparent=True)
k=π/aにおける電束密度の分布は次のようになります。
バンド1は一番下のバンド、バンド2は下から二つ目のバンドです。
低周波数側のバンド1の方が高周波数側のバンド2より高誘電率部分にエネルギーが局在していることが分かります。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/multilayer_film/GaAs_air/dfield_+.png
誘電率の分布は次の図のようになっています（黒い部分がGaAs、白い部分が空気）。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/multilayer_film/GaAs_air/epsilon.png
プロットには次のコードを用いました。
import meep as mp
from meep import mpb
from mpl_toolkits.mplot3d import Axes3D
import math
import cmath
import moviepy.editor as edit
from PIL import Image
import matplotlib.animation as animation
import sys
import matplotlib.pyplot as plt
import time
import h5py
import numpy as np
import matplotlib
import copy
import collections
matplotlib.use('Agg')
plt.rcParams["xtick.labelsize"] = 22 # デフォルトは12。tickのサイズ。
plt.rcParams["ytick.labelsize"] = 22 # デフォルトは12
plt.rcParams["axes.labelsize"] = 28 # デフォルトは12。軸のラベルのサイズ。
plt.rcParams['legend.fontsize'] = 20
plt.rcParams["figure.autolayout"] = True # ラベルが切れないように自動調整
plt.rcParams["axes.titlesize"] = 25 # 時刻を表示する以外には使用していない
plt.rcParams['font.family'] = 'sans-serif' # 'serif'がひげ飾りを意味するようで、いわゆる明朝系を表して、'sans-serif'はひげ飾りがないということでゴシック系を表しています。論文では、文章は明朝系、図はゴシック系という場合が多いのでここではゴシック系にしています。
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams["legend.fancybox"] = False # fancyboxは角の丸み
plt.rcParams["legend.edgecolor"] = 'black' # 枠線は黒色
plt.rcParams['xtick.major.pad'] = 8
plt.rcParams['ytick.major.pad'] = 8
num_bands = 8
resolution = 64
geometry_lattice = mp.Lattice(size=mp.Vector3(x=1, y=0, z=0))
geometry = [mp.Block(size=mp.Vector3(x=0.5, y=mp.inf, z=mp.inf),
material=mp.Medium(epsilon=13))]
k_points = [
mp.Vector3(x=-0.5,y=0,z=0),
mp.Vector3(x=0,y=0,z=0),
mp.Vector3(x=+0.5,y=0,z=0)
]
position=["-","0","+"]
for t in range(3):
ms = mpb.ModeSolver(
geometry=geometry,
geometry_lattice=geometry_lattice,
k_points=[k_points[t]],
resolution=resolution,
num_bands=num_bands
)
def get_dfields(ms, band):
dfields.append(ms.get_dfield(band, bloch_phase=True))
dfields = []
ms.run_tm(mpb.output_at_kpoint(k_points[t], mpb.fix_dfield_phase, get_dfields))
md = mpb.MPBData(rectify=True, resolution=64, periods=3)
converted_dfields=[]
for f in dfields:
f=f[...,0,2] # f[i][j][k][l] = 位置(i,j)の第l成分の電場(l=0,1,2がそれぞれx,y,zに対応)。kに意味はない(kは1要素しかない)。z成分のみを取り出す。
converted_dfields.append(md.convert(f))
fig = plt.figure()
for i in range(2):
f=converted_dfields[i]
ax = fig.add_subplot(121+i)
if(np.sum(np.abs(np.real(f)))>np.sum(np.abs(np.imag(f)))):
ax.imshow(np.real(f).T, interpolation='spline36', cmap="seismic")
else:
ax.imshow(np.imag(f).T, interpolation='spline36', cmap="seismic")
ax.axis("off")
ax.set_title("band "+str(i+1)+" at "+position[t])
fig.savefig("dfield_"+position[t])
fig = plt.figure()
ax = fig.add_subplot(111)
md = mpb.MPBData(rectify=True, periods=3, resolution=64)
eps = ms.get_epsilon()
converted_eps = md.convert(eps)
converted_eps2 = np.empty((len(converted_eps),len(converted_eps)),float)
for i in range(len(converted_eps2)):
for j in range(len(converted_eps2[i])):
converted_eps2[i][j]=converted_eps[i]
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(converted_eps2.T,cmap="binary")
ax.axis("off")
fig.savefig("epsilon.png",transparent=True)
***Harminvを用いた多層構造のバンド構造の計算 [#xec1100a]
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/multilayer_film/GaAs_harminv/out.png
sim.run_k_pointsの引数及び返り値の波数の単位は[&mathjax{2\pi / a};]であることに注意する(a は MeepUnit)。
import math
import meep as mp
matplotlib.use('Agg')
fcen = 0.25
df = 1.5
resolution = 32
cell = mp.Vector3(x=1, y=0, z=0)
geometry = [
mp.Block(size=mp.Vector3(x=0.5, y=mp.inf, z=mp.inf),center=mp.Vector3(x=0,y=0,z=0), material=mp.Medium(epsilon=13))]
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), component=mp.Ez,
center=mp.Vector3(x=0.1234, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[],
sources=s,
resolution=resolution)
k_interp = 19
k_point_list = mp.interpolate(k_interp, [mp.Vector3(x=0), mp.Vector3(x=0.5)])
res = sim.run_k_points(300, k_point_list)
fig=plt.figure()
ax=fig.add_subplot(111)
for i in range(len(k_point_list)):
for j in range(len(res[i])):
ax.scatter([k_point_list[i].x],[res[i][j].real],color=plt.get_cmap("tab10")(j),marker="x")
ax.set_xlim(left=0,right=0.5)
ax.set_ylim(bottom=0,top=1)
ax.set_xlabel("wave vector [2π/a]")
ax.set_ylabel("frequency [2πc/a]")
fig.savefig("out",transparent=True)
***Harminvを用いた銀の表面プラズモンのバンド構造の計算 [#hb527d3a]
場所: /abinitio/maru/cpa/spl_bulk/TM/dispersion/main.py
厚さは200nm（十分に厚く無視できると仮定）。幅は無限。単位長さは100nm。物性値はZao。
プラズマ周波数&mathjax{\omega_p=2\pi\times(c/a)\times0.7255=1.37\times10^{16}\ \mathrm{s^{-1}}};
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/surface_plasmon_bluk_TM/dispersion.png
import math
import meep as mp
fcen = 0.25
df = 0.5
resolution = 1024
sx = 10
sy = 0
sz = 0
cell = mp.Vector3(x=sx, y=sy, z=sz)
plasma_freq=0.7255
Ag_susc = [mp.DrudeSusceptibility(frequency=plasma_freq,
gamma=0.001,
sigma=1)]
Ag = mp.Medium(epsilon=1.0, E_susceptibilities=Ag_susc)
geometry = [
mp.Block(size=mp.Vector3(x=sx/2, y=mp.inf, z=mp.inf),center=mp.Vector3(x=sx/4,y=0,z=0), material=Ag)
]
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), component=mp.Hz,
center=mp.Vector3(x=0, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[mp.PML(2, direction=mp.X)],
sources=s,
resolution=resolution
)
k_interp = 20
k_point_list = mp.interpolate(k_interp, [mp.Vector3(y=0), mp.Vector3(y=2)])
res = sim.run_k_points(1000, k_point_list)
f = open("dispersion.txt", mode="w")
a=100e-9 # unit length
c=3e8 # light velocity
for i in range(len(k_point_list)):
for j in range(len(res[i])):
k = k_point_list[i].y * (2*math.pi/a)
freq = res[i][j].real * (2*math.pi*c/a)
freq_imag = res[i][j] * (2*math.pi*c/a)
Q = -res[i][j].real/(res[i][j].imag*2)
f.write(str(k)+" "+str(freq)+" "+str(freq_imag)+" "+str(Q) + "\n")
f.close()
- confinement length
- propagation length
harminv から exp(-δt) の δ が周波数の虚数成分として求まるので、伝搬長が計算できます。
まず、分散関係を多項式近似してから、微分を行うことで群速度dω/dkを求めましょう。
多項式近似は numpy の [https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html polyfit] 関数を使うと一発で出来ます。
微分は同じく numpy の [https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyder.html polyder] 関数を使うとできます。
https://kx.lumerical.com/t/can-i-simulate-spp-propagation-length-of-imi-waveguide-using-fdtd/19186
***MIM( Metal-Insulator-Metal )構造中のプラズモンのバンド構造の計算 [#cc968278]
E.N. Economou, Phy. Rev. Vol.182, 539-554 (1969)のMIM構造を計算します。
MIM(metal-insulator-metal)構造とは絶縁体薄膜を金属で挟んだ構造の事です。
プラズモンの分散関係がlight-lineと交差するため、特殊な構造を使わなくてもプラズモンを励起できる特異な構造です。
下図の×がMEEPによって計算したバンド構造です。
紫色の部分は解析解から計算したバンド構造です。
絶縁体薄膜の両面に励起されたプラズモンのエバネッセント波が互いに重なって横波の定在波を作るため、バンドが３つに分裂しています。
バンドを周波数の低い方から順にバンド１、バンド２、バンド３と名付けます。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/MIM/dispersion.png
下図では、バンド１のk_x = 1における層に垂直な電場成分をプロットしています。
横軸が位置、縦軸が時間です。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/MIM/output_efield-ex-0.31778699322664195.png
下図では、バンド２のk_x = 1における層に垂直な電場成分をプロットしています。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/MIM/output_efield-ex-0.6327492762599471.png
電場分布は金属・誘電体界面に局在し、誘電体に対してバンド１では対称、バンド２では非対称な分布になっています。
誘電体に対して非対称なプラズモンを励起するには、波源を誘電体に対して非対称に設置しなければならないことに注意してください。
この点についてメーリングリストで[https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg00809.html 議論]されています。
例えば、波源を誘電体の中心に設置すると、バンド２は得られません。
バンド２を得るには波源を誘電体の中心からわずかにずらした部分に設置する必要があります。
バンド３は金属中の横波の電磁波に対する分散関係ω^2=ω_p^2+(ck)^2の真下にあります。
このため、バンド３は金属中の横波の電磁波分布に類似した、バンド１，２に比べて界面に局在していない（電磁場の強度が小さい）解です。
fcen = 0.5
df = 2
sx = 6
sy = 0
sz = 0
resolution = 256
cell = mp.Vector3(x=sx, y=sy, z=sz)
plasma_freq=0.7255
def fun(q,freq):
eps_metal = 1-(plasma_freq/freq)**2
eps_insulator = 1
a = eps_metal*cmath.sqrt(q**2-eps_insulator*(freq**2))
b = eps_insulator * cmath.sqrt(q**2-eps_metal*(freq**2))
return min(abs((a+b)+(a-b)*cmath.exp(-0.1*2*math.pi*cmath.sqrt(q**2-(freq**2)))), abs((a+b)-(a-b)*cmath.exp(-0.1*2*math.pi*cmath.sqrt(q**2-(freq**2)))))
# plasmon polariton の分散関係（解析解）
fig=plt.figure()
ax=fig.add_subplot(111)
x = np.empty((1000,1000), float)
y = np.empty((1000,1000), float)
z = np.empty((1000,1000),float)
for i in range(1000):
for j in range(1000):
x[i][j]=i/1000.
y[i][j]=j/1000.*1.3
z[i][j]=min(0.1,fun(x[i][j],y[i][j]))
ax.pcolormesh(x,y,z,alpha=0.2)
ax.set_xlim(left=0,right=1)
ax.set_ylim(bottom=0, top=1.2)
Ag_susc = [mp.DrudeSusceptibility(frequency=plasma_freq,
gamma=0,#0.00145,
sigma=1)]
Ag = mp.Medium(epsilon=1.0, E_susceptibilities=Ag_susc)
geometry = [
mp.Block(size=mp.Vector3(x=sx, y=mp.inf,
z=mp.inf),center=mp.Vector3(x=0,y=0,z=0), material=Ag),
mp.Block(size=mp.Vector3(x=0.1, y=mp.inf,
z=mp.inf),center=mp.Vector3(x=0,y=0,z=0), material=mp.Medium(epsilon=1)),]
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), component=mp.Hz,
center=mp.Vector3(x=0.5, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[mp.PML(1, direction=mp.X)],
sources=s,
resolution=resolution)
k_interp = 19
k_point_list = mp.interpolate(k_interp, [mp.Vector3(y=0), mp.Vector3(y=1)])
res = sim.run_k_points(300, k_point_list)
for i in range(len(k_point_list)):
for j in range(len(res[i])):
#ax.scatter([k_point_list[i].y],[res[i][j].real],color=plt.get_cmap("tab10")(j),marker="x")
ax.scatter([k_point_list[i].y],[res[i][j].real],color="black",marker="x")
ax.plot([0,1],
[plasma_freq/cmath.sqrt(2),plasma_freq/cmath.sqrt(2)],color="black",linestyle="--")
#ax.fill_between([0,1],[0,1],[1,1],facecolor="gray",alpha=0.5)
ax.set_xlim(left=0,right=1)
#ax.set_ylim(bottom=0,top=0.65)
ax.set_xlabel("wave vector [2π/a]")
ax.set_ylabel("frequency [2πc/a]")
fig.savefig("dispersion")
***IMI( Insulator-Metal-Insulator )構造中のプラズモンのバンド構造の計算 [#l7bab8fc]
E.N. Economou, Phy. Rev. Vol.182, 539-554 (1969)のIMI構造を計算します。
IMI(insulator-metal-insulator)構造とは金属薄膜を絶縁体で挟んだ構造の事です。
下図の×がMEEPによって計算したバンド構造です。
紫色の部分は解析解から計算したバンド構造です。
金属薄膜両面に励起されるプラズモンのエバネッセント波が重なって横波の定在波を形成するため、バンドは2つに分裂しています。
バンドを周波数の低い方から順にバンド１、バンド２と名付けます。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/IMI/dispersion.png
下図では、バンド１のk_x = 1における層に垂直な電場成分をプロットしています。
横軸が位置、縦軸が時間です。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/IMI/output_efield-ex-0.4214435875466942.png
下図では、バンド２のk_x = 1における層に垂直な電場成分をプロットしています。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/IMI/output_efield-ex-0.566896159247526.png
電場分布は金属・誘電体界面に局在し、誘電体に対してバンド１では非対称、バンド２では対称な分布になっています。
誘電体および金属中の横波の分散関係ω=ck,ω^2=(cq)^2+(ω_p)^2上にあるプロットは界面に局在しない電磁場に相当します。
fcen = 0.5
df = 1
sx = 6
sy = 0
sz = 0
resolution = 256
cell = mp.Vector3(x=sx, y=sy, z=sz)
plasma_freq=0.7255
def fun(q,freq):
eps_metal = 1-(plasma_freq/freq)**2
eps_insulator = 1
a = eps_insulator * cmath.sqrt(q**2-eps_metal*(freq**2))
b = eps_metal*cmath.sqrt(q**2-eps_insulator*(freq**2))
return min(abs((a+b)+(a-b)*cmath.exp(-0.1*2*math.pi*cmath.sqrt(q**2-eps_metal*(freq**2)))), abs((a+b)-(a-b)*cmath.exp(-0.1*2*math.pi*cmath.sqrt(q**2-eps_metal*(freq**2)))))
# plasmon polariton の分散関係（解析解）
fig=plt.figure()
ax=fig.add_subplot(111)
x = np.empty((1000,1000), float)
y = np.empty((1000,1000), float)
z = np.empty((1000,1000),float)
for i in range(1000):
for j in range(1000):
x[i][j]=i/1000.*2
y[i][j]=j/1000.*1.3
z[i][j]=min(0.1,fun(x[i][j],y[i][j]))
ax.pcolormesh(x,y,z,alpha=0.2)
ax.hlines(y=[plasma_freq,plasma_freq/math.sqrt(2)],xmin=0,xmax=2,linestyle="--")
ax.set_xlim(left=0,right=2)
ax.set_ylim(bottom=0, top=1.2)
Ag_susc = [mp.DrudeSusceptibility(frequency=plasma_freq,
gamma=0,#0.00145,
sigma=1)]
Ag = mp.Medium(epsilon=1.0, E_susceptibilities=Ag_susc)
geometry = [
mp.Block(size=mp.Vector3(x=0.1, y=mp.inf, z=mp.inf),center=mp.Vector3(x=0,y=0,z=0), material=Ag)]
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), component=mp.Hz,
center=mp.Vector3(x=0.1, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[mp.PML(1, direction=mp.X)],
sources=s,
resolution=resolution)
k_interp = 19
k_point_list = mp.interpolate(k_interp, [mp.Vector3(y=0), mp.Vector3(y=2)])
res = sim.run_k_points(500, k_point_list)
for i in range(len(k_point_list)):
for j in range(len(res[i])):
#ax.scatter([k_point_list[i].y],[res[i][j].real],color=plt.get_cmap("tab10")(j),marker="x")
ax.scatter([k_point_list[i].y],[res[i][j].real],color="black",marker="x")
ax.plot([0,1], [plasma_freq/cmath.sqrt(2),plasma_freq/cmath.sqrt(2)],color="black",linestyle="--")
#ax.fill_between([0,1],[0,1],[1,1],facecolor="gray",alpha=0.5)
ax.set_xlim(left=0,right=1)
#ax.set_ylim(bottom=0,top=0.65)
ax.set_xlabel("wave vector [2π/a]")
ax.set_ylabel("frequency [2πc/a]")
fig.savefig("dispersion", bbox_inches="tight")
***誘電体円柱の正方格子 [#m892c5e3]
[http://ab-initio.mit.edu/book/photonic-crystals-book.pdf Photonic Crystals: Molding the Flow of Light (second edition) ]の67ページに掲載されている構造をシミュレーションします。
アルミナ(誘電率8.9)の円柱が正方格子状に並んでいるフォトニック結晶のバンド構造を計算します。
円柱の半径rと正方格子の一辺の長さaがr=0.2aの関係にあるとします。
誘電率の分布は次のようになっています。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_of_dielectric_columns/epsilon.png
バンド構造は次のようになります。
青がTMモード、赤がTEモードです。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_of_dielectric_columns/dispersion.png
バンド構造の極大値、極小値（バンドギャップの大きさにおいて重要）はブリルアンゾーンの端、特に角で出現することが多いため、ブリルアンゾーンの端のみをプロットしています。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_of_dielectric_columns/real_and_reciprocal_space.png
Γ, X, M 点は上のように定めています。
バンド構造の計算には次のコードを用いました。
import meep as mp
from meep import mpb
from mpl_toolkits.mplot3d import Axes3D
import math
import cmath
import moviepy.editor as edit
from PIL import Image
import matplotlib.animation as animation
import sys
import matplotlib.pyplot as plt
import time
import h5py
import numpy as np
import matplotlib
import copy
import collections
matplotlib.use('Agg')
plt.rcParams["xtick.labelsize"] = 22 # デフォルトは12。tickのサイズ。
plt.rcParams["ytick.labelsize"] = 22 # デフォルトは12
plt.rcParams["axes.labelsize"] = 28 # デフォルトは12。軸のラベルのサイズ。
plt.rcParams['legend.fontsize'] = 20
plt.rcParams["figure.autolayout"] = True # ラベルが切れないように自動調整
plt.rcParams["axes.titlesize"] = 25 # 時刻を表示する以外には使用していない
plt.rcParams['font.family'] = 'sans-serif' # 'serif'がひげ飾りを意味するようで、いわゆる明朝系を表して、'sans-serif'はひげ飾りがないということでゴシック系を表しています。論文では、文章は明朝系、図はゴシック系という場合が多い のでここではゴシック系にしています。
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams["legend.fancybox"] = False # fancyboxは角の丸み
plt.rcParams["legend.edgecolor"] = 'black' # 枠線は黒色
plt.rcParams['xtick.major.pad'] = 8
plt.rcParams['ytick.major.pad'] = 8
num_bands = 8
resolution = 32
geometry_lattice = mp.Lattice(size=mp.Vector3(x=1, y=1))
geometry = [mp.Cylinder(radius=0.2, material=mp.Medium(epsilon=8.9))]
k_points = [
mp.Vector3(x=0,y=0), # Gamma
mp.Vector3(x=0.5,y=0), # M
mp.Vector3(x=0.5, y=0.5), # K
mp.Vector3(x=0,y=0), # Gamma
]
k_points = mp.interpolate(4, k_points)
ms = mpb.ModeSolver(
geometry=geometry,
geometry_lattice=geometry_lattice,
k_points=k_points,
resolution=resolution,
num_bands=num_bands
)
ms.run_tm(mpb.output_at_kpoint(mp.Vector3(x=0, y=0), mpb.fix_efield_phase, mpb.output_efield_z))
tm_freqs = ms.all_freqs
tm_gaps = ms.gap_list
ms.run_te()
te_freqs = ms.all_freqs
te_gaps = ms.gap_list
fig=plt.figure()
ax=fig.add_subplot(111)
for i in range(len(tm_gaps)):
if(tm_gaps[i][0]<1):
continue
ax.axhspan(tm_gaps[i][1],tm_gaps[i][2],color="blue",alpha=0.2)
for i in range(len(te_gaps)):
if(te_gaps[i][0]<1):
continue
ax.axhspan(te_gaps[i][1],te_gaps[i][2],color="red",alpha=0.2)
ax.plot(tm_freqs,color="blue")
ax.plot(te_freqs,color="red")
for i in range(len(tm_freqs)):
for freq in tm_freqs[i]:
ax.scatter(i,freq,color="blue")
for i in range(len(tm_freqs)):
for freq in te_freqs[i]:
ax.scatter(i,freq,color="red")
ax.set_xticks([0,5,10,15])
ax.set_xticklabels(["Γ","X","M","Γ"])
ax.set_xlim([0,15])
ax.set_ylim([0,1])
ax.set_ylabel("frequency [c/a]")
fig.savefig("dispersion",transparent=True)
fig = plt.figure()
ax = fig.add_subplot(111)
md = mpb.MPBData(rectify=True, periods=3, resolution=32)
ax.imshow(md.convert(ms.get_epsilon()),cmap="binary")
ax.axis("off")
fig.savefig("epsilon.png",transparent=True)
周波数0.4 [c/a]付近にTMモードのバンドギャップがあります。
一つ目のバンド(dielectric band)と二つ目のバンド(air band)におけるTMモードの電束密度の分布を次に示します。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_of_dielectric_columns/dfield_gamma.png
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_of_dielectric_columns/dfield_X.png
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_of_dielectric_columns/dfield_M.png
一般に、低周波数のバンドほど電場のエネルギーが高誘電率の場所に局在します。
TEモードでは、電束密度は円柱に対して垂直に走っており、必ず円柱外（低誘電率部分）を通る必要があるため高誘電率部分にのみ局在することができません。
したがって、TEモードではTMモードと比べてバンドギャップが小さくなります（詳しくは[http://ab-initio.mit.edu/book/photonic-crystals-book.pdf Photonic Crystals: Molding the Flow of Light (second edition) ]の71ページ）。
以下のコードを用いました。
mpb.ModeSolverにk_pointsを一要素ずつ渡さないと、電束密度の分布が壊れることに注意してください（原因不明）。
import meep as mp
from meep import mpb
from mpl_toolkits.mplot3d import Axes3D
import math
import cmath
import moviepy.editor as edit
from PIL import Image
import matplotlib.animation as animation
import sys
import matplotlib.pyplot as plt
import time
import h5py
import numpy as np
import matplotlib
import copy
import collections
matplotlib.use('Agg')
plt.rcParams["xtick.labelsize"] = 22 # デフォルトは12。tickのサイズ。
plt.rcParams["ytick.labelsize"] = 22 # デフォルトは12
plt.rcParams["axes.labelsize"] = 28 # デフォルトは12。軸のラベルのサイズ。
plt.rcParams['legend.fontsize'] = 20
plt.rcParams["figure.autolayout"] = True # ラベルが切れないように自動調整
plt.rcParams["axes.titlesize"] = 25 # 時刻を表示する以外には使用していない
plt.rcParams['font.family'] = 'sans-serif' # 'serif'がひげ飾りを意味するようで、いわゆる明朝系を表して、'sans-serif'はひげ飾りがないということでゴシック系を表しています。論文では、文章は明朝系、図はゴシック系という場合が多いのでここではゴシック系にしています。
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams["legend.fancybox"] = False # fancyboxは角の丸み
plt.rcParams["legend.edgecolor"] = 'black' # 枠線は黒色
plt.rcParams['xtick.major.pad'] = 8
plt.rcParams['ytick.major.pad'] = 8
num_bands = 8
resolution = 64
geometry_lattice = mp.Lattice(size=mp.Vector3(x=1, y=1))
geometry = [mp.Cylinder(radius=0.2, material=mp.Medium(epsilon=8.9))]
k_points = [
mp.Vector3(x=0.5, y=0.5), # M
mp.Vector3(x=0,y=0,z=0), # Gamma
mp.Vector3(x=0.5,y=0,z=0) # X
]
position=["M","Gamma","X"]
for t in range(3):
ms = mpb.ModeSolver(
geometry=geometry,
geometry_lattice=geometry_lattice,
k_points=[k_points[t]],
resolution=resolution,
num_bands=num_bands
)
def get_dfields(ms, band):
dfields.append(ms.get_dfield(band, bloch_phase=True))
dfields = []
ms.run_tm(mpb.output_at_kpoint(k_points[t], mpb.fix_dfield_phase, get_dfields))
md = mpb.MPBData(rectify=True, resolution=64, periods=3)
converted_dfields=[]
for f in dfields:
f=f[...,0,2] # f[i][j][k][l] = 位置(i,j)の第l成分の電場(l=0,1,2がそれぞれx,y,zに対応)。kに意味はない(kは1要素しかない)。z成分のみを取り出す。
converted_dfields.append(md.convert(f))
fig = plt.figure()
for i in range(2):
f=converted_dfields[i]
ax = fig.add_subplot(121+i)
if(np.sum(np.abs(np.real(f)))>np.sum(np.abs(np.imag(f)))):
ax.imshow(np.real(f).T, interpolation='spline36', cmap="seismic")
else:
ax.imshow(np.imag(f).T, interpolation='spline36', cmap="seismic")
ax.axis("off")
ax.set_title("band "+str(i+1)+" at "+position[t])
fig.savefig("dfield_"+position[t])
*** square lattice of dielectric veins [#ma1f1c0e]
[http://ab-initio.mit.edu/book/ Photonic Crystals: Molding the Flow of Light Page72]の図5, 6, 7の再現。
import meep as mp
from meep import mpb
import matplotlib.pyplot as plt
import numpy as np
import matplotlib
matplotlib.use('Agg')
plt.rcParams["xtick.labelsize"] = 22 # デフォルトは12。tickのサイズ。
plt.rcParams["ytick.labelsize"] = 22 # デフォルトは12
plt.rcParams["axes.labelsize"] = 28 # デフォルトは12。軸のラベルのサイズ。
plt.rcParams['legend.fontsize'] = 20
plt.rcParams["figure.autolayout"] = True # ラベルが切れないように自動調整
plt.rcParams["axes.titlesize"] = 25 # 時刻を表示する以外には使用していない
plt.rcParams['font.family'] = 'sans-serif' # 'serif'がひげ飾りを意味するようで、いわゆる明朝系を表して、'sans-serif'はひげ飾りがないということでゴシック系を表しています。論文では、文章は明朝系、図はゴシック系という場合が多いのでここではゴシック系にしています。
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams["legend.fancybox"] = False # fancyboxは角の丸み
plt.rcParams["legend.edgecolor"] = 'black' # 枠線は黒色
plt.rcParams['xtick.major.pad'] = 8
plt.rcParams['ytick.major.pad'] = 8
sx=1
sy=1
num_bands = 4
resolution = 200
geometry_lattice = mp.Lattice(size=mp.Vector3(x=sx, y=sy))
geometry = [
mp.Block(center=mp.Vector3(), size=mp.Vector3(x=sx*0.165, y=sy), material=mp.Medium(epsilon=8.9)),
mp.Block(center=mp.Vector3(), size=mp.Vector3(x=sx, y=sy*0.165), material=mp.Medium(epsilon=8.9))
]
k_points = [
mp.Vector3(x=0,y=0,z=0), # Gamma
mp.Vector3(x=0.5,y=0,z=0), # X
mp.Vector3(x=0.5, y=0.5), # M
mp.Vector3(x=0,y=0,z=0) # Gamma
]
k_points = mp.interpolate(4, k_points)
ms = mpb.ModeSolver(
geometry=geometry,
geometry_lattice=geometry_lattice,
k_points=k_points,
resolution=resolution,
num_bands=num_bands
)
ms.run_tm()
tm_freqs = ms.all_freqs
tm_gaps = ms.gap_list
ms.run_te()
te_freqs = ms.all_freqs
te_gaps = ms.gap_list
fig=plt.figure()
ax=fig.add_subplot(111)
#================================================
# plotting gap
#------------------------------------------------
for i in range(len(tm_gaps)):
if(tm_gaps[i][0]<1):
continue
ax.axhspan(tm_gaps[i][1],tm_gaps[i][2],color="blue",alpha=0.2)
for i in range(len(te_gaps)):
if(te_gaps[i][0]<1):
continue
ax.axhspan(te_gaps[i][1],te_gaps[i][2],color="red",alpha=0.2)
#================================================
ax.plot(tm_freqs,color="blue")
ax.plot(te_freqs,color="red")
for i in range(len(tm_freqs)):
for freq in tm_freqs[i]:
ax.scatter(i,freq,color="blue")
for i in range(len(tm_freqs)):
for freq in te_freqs[i]:
ax.scatter(i,freq,color="red")
ax.set_xticks([0,5,10,15])
ax.set_xticklabels(["Γ","X","M","Γ"])
ax.set_xlim([0,15])
ax.set_ylim([0,0.8])
ax.set_ylabel("frequency [c/a]")
fig.savefig("dispersion",transparent=True)
for t in range(2):
fields = []
def get_fields(ms, band): # bandはバンドのインデックス
if(t==0):
fields.append(ms.get_hfield(band, bloch_phase=True))
elif(t==1):
fields.append(ms.get_dfield(band, bloch_phase=True))
if(t==0):
ms.run_te(mpb.output_at_kpoint(mp.Vector3(x=0.5, y=0), mpb.fix_dfield_phase, get_fields))
elif(t==1):
ms.run_tm(mpb.output_at_kpoint(mp.Vector3(x=0.5, y=0), mpb.fix_dfield_phase, get_fields))
md = mpb.MPBData(rectify=True, resolution=200, periods=3)
converted_dfields=[]
for f in fields:
f=f[...,0,2] # f[i][j][k][l] = 位置(i,j)の第l成分の電場(l=0,1,2がそれぞれx,y,zに対応)。kに意味はない(kは1要素しかない)。z成分のみを取り出す。
converted_dfields.append(md.convert(f))
fig = plt.figure()
ax = fig.add_subplot(111)
for i in range(4):
f = converted_dfields[i]
if(np.sum(np.abs(np.real(f)))>np.sum(np.abs(np.imag(f)))):
ax.imshow(np.real(f).T, interpolation='spline36', cmap="RdBu")
else:
ax.imshow(np.imag(f).T, interpolation='spline36', cmap="RdBu")
ax.axis("off")
if(t==0):
fig.savefig("te_hfield"+str(i+1), transparent=True)
elif(t==1):
fig.savefig("tm_dfield"+str(i+1), transparent=True)
ax.clear()
ax.imshow(mpb.MPBData(periods=3, resolution=200).convert(ms.get_epsilon()).T ,cmap="binary") # x,y の定義がmatplotlibとMEEPで逆なので、転置(.T)する。
ax.axis("off")
fig.savefig("epsilon.png",transparent=True)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_veins/dispersion.png
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_veins/tm_dfield1.png
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_veins/tm_dfield2.png
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_veins/te_hfield1.png
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/square_lattice_veins/te_hfield2.png
***triangular lattice of air columns [#q58923ea]
[http://ab-initio.mit.edu/book/ Photonic Crystals: Molding the Flow of Light Page76]の図9, 10, 11の再現。
import meep as mp
from meep import mpb
from mpl_toolkits.mplot3d import Axes3D
import math
import cmath
import moviepy.editor as edit
from PIL import Image
import matplotlib.animation as animation
import sys
import matplotlib.pyplot as plt
import time
import h5py
import numpy as np
import matplotlib
import copy
import collections
matplotlib.use('Agg')
plt.rcParams["xtick.labelsize"] = 22 # デフォルトは12。tickのサイズ。
plt.rcParams["ytick.labelsize"] = 22 # デフォルトは12
plt.rcParams["axes.labelsize"] = 28 # デフォルトは12。軸のラベルのサイズ。
plt.rcParams['legend.fontsize'] = 20
plt.rcParams["figure.autolayout"] = True # ラベルが切れないように自動調整
plt.rcParams["axes.titlesize"] = 25 # 時刻を表示する以外には使用していない
plt.rcParams['font.family'] = 'sans-serif' # 'serif'がひげ飾りを意味するようで、いわゆる明朝系を表して、'sans-serif'はひげ飾りがないということでゴシック系を表しています。論文では、文章は明朝系、 図はゴシック系という場合が多いのでここではゴシック系にしています。
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams["legend.fancybox"] = False # fancyboxは角の丸み
plt.rcParams["legend.edgecolor"] = 'black' # 枠線は黒色
plt.rcParams['xtick.major.pad'] = 8
plt.rcParams['ytick.major.pad'] = 8
sx=1
sy=1
num_bands = 8
resolution = 200
geometry_lattice = mp.Lattice(size=mp.Vector3(x=sx, y=sy),
basis1=mp.Vector3(math.sqrt(3)/2, 0.5),
basis2=mp.Vector3(math.sqrt(3)/2, -0.5)
)
geometry = [
mp.Block(center=mp.Vector3(), size=mp.Vector3(x=sx, y=sy), material=mp.Medium(epsilon=13)),
mp.Cylinder(center=mp.Vector3(), radius=0.48, material=mp.Medium(epsilon=1))
]
k_points = [
mp.Vector3(), # Gamma
mp.Vector3(y=0.5), # X
mp.Vector3(x=-1/3, y=1/3), # M
mp.Vector3() # Gamma
]
k_points = mp.interpolate(4, k_points)
ms = mpb.ModeSolver(
geometry=geometry,
geometry_lattice=geometry_lattice,
k_points=k_points,
resolution=resolution,
num_bands=num_bands
)
ms.run_tm()
tm_freqs = ms.all_freqs
tm_gaps = ms.gap_list
ms.run_te()
te_freqs = ms.all_freqs
te_gaps = ms.gap_list
fig=plt.figure()
ax=fig.add_subplot(111)
#================================================
# plotting gap
#------------------------------------------------
for i in range(len(tm_gaps)):
if(tm_gaps[i][0]<1):
continue
ax.axhspan(tm_gaps[i][1],tm_gaps[i][2],color="blue",alpha=0.2)
for i in range(len(te_gaps)):
if(te_gaps[i][0]<1):
continue
ax.axhspan(te_gaps[i][1],te_gaps[i][2],color="red",alpha=0.2)
#================================================
ax.plot(tm_freqs,color="blue")
ax.plot(te_freqs,color="red")
for i in range(len(tm_freqs)):
for freq in tm_freqs[i]:
ax.scatter(i,freq,color="blue")
for i in range(len(tm_freqs)):
for freq in te_freqs[i]:
ax.scatter(i,freq,color="red")
ax.set_xticks([0,5,10,15])
ax.set_xticklabels(["Γ","M","K","Γ"])
ax.set_xlim([0,15])
ax.set_ylim([0,0.8])
ax.set_ylabel("frequency [c/a]")
fig.savefig("dispersion",transparent=True)
ax.clear()
ax.imshow(mpb.MPBData(periods=3, resolution=200, rectify=True).convert(ms.get_epsilon()).T ,cmap="binary") # x,y の定義がmatplotlibとMEEPで逆なので、転置(.T)する。
ax.axis("off")
fig.savefig("epsilon.png",transparent=True)
***状態密度（DOS） [#sa17a2a4]
https://www.mail-archive.com/mpb-discuss@ab-initio.mit.edu/msg00041.html
f=[0, df, 2df, 3df, ...]の各周波数についてDOSを求めて、プロットしたい。
周波数f[i]におけるDOSは次の式で近似できる。
#mathjax( \sum_j{ \frac{1}{\sqrt{\pi}} \exp\left(-\left(\frac{f_{resonant}[j]-f[i]}{df}\right)^2 \right) } );
**共振モードの計算 [#n9e267ed]
短パルスに対する応答を計算することで、調和振動する共振モードを見つけることができます。
短パルスを用いるのは、短パルスが周波数領域で広がっているため（さまざまな周波数の波を含むため）です。
***円柱型の誘電体の共振モード [#g98983fa]
MEEP公式の[https://meep.readthedocs.io/en/latest/Python_Tutorials/Basics/#modes-of-a-ring-resonator Tutorial]に従って、円柱型の誘電体の共振モードを計算します。
円柱はz軸方向に無限に伸びているとします。
電場がz方向に偏向したモードを計算します。
共振モードの周波数が分からないため、中心周波数1,周波数幅0.1のz方向に偏向したガウス波源を用いて、共振モードを励起します。
円柱の屈折率は3.4です。
円柱の周りは真空です。
ガウス波源が十分小さくなった後に残っている電磁場がHarminvによって、調和振動した共振モード群に分離されます。
次のrun(step_functions...,until_after_sources=time)関数は全ての波源がturn-offになるまでシミュレーションを行います。
run(step_functions..., until_after_sources=time)
また、各タイムステップ毎にstep_functions...を呼び出します。
until_after_sourcesに渡した数(time)だけ、全ての波源がturn-offになった後もシミュレーションを行います。
次のafter_sources関数はstep_functionの一種です。ラップしたstep_functionを波源がturn-offされているときのみ呼び出します。
after_sources(step_functions...)
次のHarminv関数は、点ptにおける電磁場成分c(例えばmp.Ez)を記録し、[fcen-df/2,fcen+df/2]の範囲で記録した電磁場を
f(t)=Σa_j exp(-i*omega_j*t+δt)
の形に分解します。
Harminv(c, pt, fcen, df)
Harminvは次のような標準出力を行います。
harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error
harminv0:, 0.11678193348699654, -0.0006938113872636076, 84.15971230133923, 0.0035411967397259462, -0.002461536200934776-0.002545763909111224i, 7.611818184260265e-06+0.0i
harminv0:, 0.14635420376193028, -0.00019625091272617196, 372.87521807895405, 0.0317495337301081, 0.014496676972289664+0.028246756430435668i, 2.4460004142673657e-08+0.0i
harminv0:, 0.17476005848070048, -4.962385886158849e-05, 1760.8471256552566, 0.008460040745912884, 0.0003936115065882616-0.008450879208957346i, 8.162174597685e-08+0.0i
frequency=omega[2πc], imag.freq=δ, Q = omega/(-2δ), |a_n| = abs(a_j), a_n = a_j, error = frequencyとimag.freqの大まかな誤差に相当します。 Q>>frequencyまたQ>>imag.freqのときはQの値は正確ではありません。
得られた共振周波数での実際の電磁場の分布を調べたいときは、fcenをその値に設定しdfを小さくしたうえでシミュレーションを実行してください。
TMモードの共振モードとして次の電場（Ez）分布が見つかりました。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/cylinder_resonant/fcen_0.1168564818622051_TM.gif
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/cylinder_resonant/fcen_0.14634615459699182_TM.gif
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/cylinder_resonant/fcen_0.17475492815676294_TM.gif
また、TEモードの共振モードとして次の磁場（Hz）分布が見つかりました。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/cylinder_resonant/fcen_0.14275359824538983_TE.gif
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/cylinder_resonant/fcen_0.17171322951171442_TE.gif
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/cylinder_resonant/fcen_0.19967038987644908_TE.gif
z方向には周期境界条件を用いているため、k_z=0なる共振モードのみが見つかっていることに注意してください(k_z≠0のものを見つけたいときはブロッホの周期境界条件を用いてください)。
TMモードのEzは
E_z = J_n(hr)cos(nΦ)exp(jωt-γz)
TEモードのHzは
H_z = J_n(hr)cos(nΦ)exp(jωt-γz)
と解析的に表せます。ここでγ=sqrt(h^2 - (ω/c)^2)かつhはJ_n(ha)=0の第p個目の解です([http://demonstrations.wolfram.com/CylindricalCavityResonator/ 参考])。
上のシミュレーションではp=1のもののみが見つかっています。
TMの共振モードの探索に使ったコードは次の通りです。
(TEの共振モードの探索には、下のコードのEzをHzに置き換えるだけでよい)
n = 3.4
w = 1
r = 1
pad = 4
dpml = 2
sxy = 2*(r+w+pad+dpml)
fcen = 0.15
df = 0.1
geometry=[
mp.Cylinder(radius=r+w, material=mp.Medium(index=n))
]
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), mp.Hz, mp.Vector3(r+0.1))]
sim = mp.Simulation(cell_size=mp.Vector3(sxy, sxy),
geometry=geometry,
sources=sources,
resolution=10,
boundary_layers=[mp.PML(dpml)])
harminv_instance=mp.Harminv(c=mp.Hz, pt=mp.Vector3(x=r+0.1,y=0), fcen=fcen, df=df)
sim.run(mp.after_sources(harminv_instance),
until_after_sources=300)
for mode in harminv_instance.modes:
fcen=mode.freq
df=0.05
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), mp.Hz, mp.Vector3(r+0.1))]
sim = mp.Simulation(cell_size=mp.Vector3(sxy, sxy),
geometry=geometry,
sources=sources,
resolution=10,
boundary_layers=[mp.PML(dpml)])
sim.run(until_after_sources=300)
hfields=[]
def get_hfields(sim):
hfields.append(sim.get_hfield_z())
sim.run(mp.at_every(1/fcen/20, get_hfields), until=1/fcen)
eps_data = sim.get_epsilon()
fig = plt.figure()
ax = fig.add_subplot(111)
ims=[]
ma=np.abs(hfields).max()
for i in range(len(hfields)):
im = ax.imshow(np.transpose(np.real(hfields[i])), interpolation='spline36', cmap='seismic', alpha=0.9, vmin=-ma, vmax=ma)
ax.axis('off')
ims.append([im])
ani = animation.ArtistAnimation(fig,ims,interval=100)
ani.save("fcen_"+str(fcen)+"_TE.gif",writer="imagemagick")
** frequency-domain solver ( 単一周波数の光源に対する応答を計算 ) [#n4f5e86f]
MEEPの[https://meep.readthedocs.io/en/latest/Python_Tutorials/Frequency_Domain_Solver/ Tutorial]に従って、frequency-domain solverを動かしてみます。
frequency-domain solverでは単一周波数の光源に対する応答を計算します。
FDTD法( time-domain solver )で単一周波数の光源に対する応答を計算しようと思っても、
波源を立ち上げる時に目的外の周波数の波が混ざってしまいます。
time-domain solverでは、この目的外の周波数の波が励起したモードが減衰するまで、（場合にはよっては長時間）シミュレーションを行う必要があります。
frequency-domain solverでは原理的に単一周波数の光源に対する応答しか計算しないので、目的外の周波数が励起したモードが混ざる可能性がありません。
frequency-domain solverの定式化については[http://ab-initio.mit.edu/~oskooi/papers/Oskooi10.pdf Computer Physics Communications 181 (2010) 687–702]の5.3. Frequency-domain solverを参照してください。
+波源 ContinuousSrc を使う。言い換えれば、GaussianSrc を使わない。
+ mp.Simulation() に force_complex_fields=True を渡すことで、電磁場に複素数を使う。
+シミュレーションの前に sim.init_sim() を呼び出して初期化を行う
+sim.solve_cm(tl, maxiters, L) によって、frequency-domain solverを呼び出す。tol は 行列を反復法によって解くときのエラーの許容値（デフォルトは 10^{-8} ）。maxiters は 反復の最大回数（デフォルトは 10^4 ）。L は反復法のメモリと計算量のトレードオフを決める値で、L が大きいほどメモリを使用する。L = 10 が推奨値。
Lorenz-Drude伝導は frequency-domain solver では使用できない。
**円筒座標におけるシミュレーション [#u3c69825]
[https://meep.readthedocs.io/en/latest/Python_Tutorials/Ring_Resonator_in_Cylindrical_Coordinates/ Tutorial]に従って、円筒座標においてシミュレーションを実行してみます。
mp.Simulation に dimensions = mp.CYLINDRICAL を設定すると、円筒座標系でシミュレーションが行われます。
系が回転対称性を持つ時、電磁場分布は exp(mΦ) の依存性を持ちます（一周すると元の位相に戻ってくる）。
mp.Simulation のオプション m で、この動径方向の依存性を設定できます（デフォルトは m = 3）。
cell は 直交座標系のときと同様に mp.Vector3(x, y, z) で定めますが、x, y, z はr, Φ, z に対応していることに注意してください。
平面の吸収層として設計された PML を曲げて使っているため、反射が直交座標系のときよりも強くなります。
cell の半径を大きくすると、PMLが平面に近づき反射は小さくなります。
**Q&A [#y2d911b6]
***一つのファイル内で複数回シミュレーションをすると計算結果がおかしくなる[#b9c73d94]
1つのファイル内で複数回シミュレーションをするときは、各シミュレーションを終了するたびに(reset-meep) を呼び出さなければなりません。
(reset-meep)を呼び出さない場合、正しい計算結果が得られません。
***（入射波＋反射波）のフラックスが入射波のみの時より、小さくなる [#ua9ad9fe]
フラックスは座標軸の正の方向を正に取っているためです（ポインティングベクトルで定義されている）。
正の方向に入射波が進行し、負の方向に反射波が進行するとき、入射波だけの時よりフラックスは小さくなります。
***dielectricよりmediumでindexを指定した方が計算が早い（計算結果は同じなのに） [#xe0c2da7]
理由は不明です。
***resolutionをsetしても変化しない [#g280b7dd]
resolutionはset! geometry-latticeの前に置かないとデフォルト値10になってしまうことがあります。
***変数vを使っていないのにunbounded variable vというエラーが出る [#ubd3a6e9]
(define lambda 5)とするとunbounded variable vとしてエラーが出る。
Schemeではlambdaを用いて関数を定義するからです（予約語になっています）。
***meep: Could not determine normal direction for given grid_volume. [#ta6b98ee]
fluxを3次元で設定していることが原因です。
fluxは2次元の面または1次元の点で設定しなければなりません。
***実行時にGUILE_WARN_DEPRECATEDという警告が出る [#ma3f1cad]
放っておいても問題ありませんが、.loginに以下を記述すると、非推奨(deprecated)な機能を使用しているというエラーを握りつぶせます。
setenv GUILE_WARN_DEPRECATED no
***unbound variable: とエラーが出て、しかも : 以降が空白になっている [#i0e856b1]
「unbound variable: 」とエラーが出てるにも関わらず、unbound variableの名前が表示されないときは、全角空白文字がプログラムのどこかに含まれています。
全角空白を削除すればよいです。
***コマンドラインから変数を渡しているのに無視される。 [#dbd0a076]
meep hoehoge.ctl var=1
とすると、var=1は無視されます。
meep var=1 hogehoge.ctl
の順にしてください。
***直方体を斜めに配置したい [#x6b589ab]
e1,e2,e3[vector3]により、方向を定めることができます。
例えば、xy平面内でx軸からΘだけ傾けて置くには、次のようにしてください。
(make block
(size a b c)
(e1 (vector3 (cos theta) (sin theta) 0))
(e2 (vector3 (- (sin theta)) (cos theta) 0)
(e3 (vector3 0 0 1))
)
***x,y,z方向のPMLについて、厚さをそれぞれ別々にしたい [#z61ce405]
X, Y, Z方向のPMLをそれぞれ厚さa, b, cにするには
(set! pml-layers (list (make pml (thickness a)(direction X))))
(set! pml-layers (list (make pml (thickness b)(direction Y))))
(set! pml-layers (list (make pml (thickness c)(direction Z))))
としてください。
***負の透過率、反射率、吸収率を示す [#j49576fc]
負の透過率、反射率、吸収率は電磁波がPMLの外に発散しきっていないのにシミュレーションを中断したことが原因のことが多いです。
特にstop-when-fields-decayed関数を使用しているときは注意してください。
(run-sources+ T
(stop-when-fields-decayed dt Ex (vector3 x y z) 1e-3))
とすると、T秒間シミュレーションした後に、(x,y,z)におけるdt秒間毎のExの最大値が、全ての時間を通しての最大値と比較して10の-3乗倍以下になるまでシミュレーションを続けてくれます。
しかしこの関数は、長時間共鳴するような構造の場合は、まだ、電磁波の発振が続いているのにシミュレーションを中断する場合があります。
電磁波が発散しきっているかどうかの判断は、h5topngなどで電磁波を画像化して直接確認するしかありません。
***伝導度の理論値を複数のローレンツの重ね合わせで表したい [#p43befa3]
下のリンクの方法で取り組んでいました。
理論値を複数のローレンツの重ね合わせで合わせるのはよくなさそうです。
https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg05147.html
***インターフェイスはPythonとSchemeどちらを使うべきか [#m1cc6868]
マニュアルのトップページにSchemeがobsoleteでPythonに変えろと書いてあったことと、メーリングリストで（バグはまだ多いですが）Pythonに機能が追加されていることもあって、先週からPythonに移行のために動作検証しているところです。
***PMLの厚さはどのくらいにすべきか[#l8406a8a]
電磁波がPMLに垂直入射するときはPMLの厚さを半波長にしてください。
以下、MEEP作成者のPMLに関するドキュメント( http://math.mit.edu/~stevenj/18.369/pml.pdf )からの引用。
With PML, however, the constant factor is very good to start with, so experience shows that a simple quadratic or cubic turn-on of the PML absorption usually produces negligible refections for a PML layer of only half a wavelength or thinner [1,21].
(Increasing the resolution also increases the resolution alos increase the effectiveness of the PML, because it approaches the exact wave equation.)
斜め入射のときはPMLの吸収率が減少するので、厚みを増やす必要があります。
反射率は入射角をXとすると、ある定数δを用いて
e^{-2*k*d*cos(X)*δ}
と書けます。
dはPMLの厚さ、kは波数。
詳細は[https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=477075 UPMLの論文]を参照してください。
***tubeとflexにおいてmaterials.pyが使えない[#m4142e79]
PyMeepのバージョンが1.5だと次のエラーによりmaterials.pyが実行できません。
tube60:~$ python
Python 3.6.7 | packaged by conda-forge | (default, Nov 21 2018, 03:09:43)
[GCC 7.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import meep as mp
>>> from meep.materials import Au
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'meep.materials'
PyMeepのバージョンは
conda list
から確認できます。
次のようにオプションをつけてPyMeepを1.5から1.6にアップデートするとmaterials.pyが実行できるようになります。
(mp) maru@flex:~$ conda update -c chogan -c conda-forge pymeep
Solving environment: done
## Package Plan ##
environment location: /home/students/maru/anaconda3/envs/mp
added / updated specs:
- pymeep
The following packages will be downloaded:
package | build
---------------------------|-----------------
libctl-4.1.4 | 1 97 KB chogan
libgdsii-0.2.dev | 0 800 KB chogan
python-3.6.6 | hd21baee_1003 29.0 MB conda-forge
openssl-1.0.2p | h14c3975_1002 3.1 MB conda-forge
mpb-1.7.0 | nomkl_1 74 KB chogan
pymeep-1.7.0 | py36_nomkl_1 3.3 MB chogan
openblas-0.3.5 | h9ac9557_1000 15.8 MB conda-forge
nomkl-3.0 | 0 48 KB
ca-certificates-2018.11.29 | ha4d7672_0 143 KB conda-forge
certifi-2018.11.29 | py36_1000 145 KB conda-forge
------------------------------------------------------------
Total: 52.5 MB
The following NEW packages will be INSTALLED:
libgdsii: 0.2.dev-0 chogan
nomkl: 3.0-0
The following packages will be UPDATED:
ca-certificates: 2018.03.07-0 --> 2018.11.29-ha4d7672_0 conda-forge
certifi: 2018.11.29-py36_0 --> 2018.11.29-py36_1000 conda-forge
libctl: 4.1.2-0 chogan --> 4.1.4-1 chogan
mpb: 1.6.2-4 chogan --> 1.7.0-nomkl_1 chogan [nomkl]
openblas: 0.2.20-8 conda-forge --> 0.3.5-h9ac9557_1000 conda-forge
pymeep: 1.5-py36_nomkl_2 chogan [nomkl] --> 1.7.0-py36_nomkl_1 chogan [nomkl]
The following packages will be DOWNGRADED:
openssl: 1.1.1a-h7b6447c_0 --> 1.0.2p-h14c3975_1002 conda-forge
python: 3.6.8-h0371630_0 --> 3.6.6-hd21baee_1003 conda-forge
Proceed ([y]/n)? y
Downloading and Extracting Packages
libctl-4.1.4 | 97 KB | ########################################################################################### | 100%
libgdsii-0.2.dev | 800 KB | ########################################################################################### | 100%
python-3.6.6 | 29.0 MB | ########################################################################################### | 100%
openssl-1.0.2p | 3.1 MB | ########################################################################################### | 100%
mpb-1.7.0 | 74 KB | ########################################################################################### | 100%
pymeep-1.7.0 | 3.3 MB | ########################################################################################### | 100%
openblas-0.3.5 | 15.8 MB | ########################################################################################### | 100%
nomkl-3.0 | 48 KB | ########################################################################################### | 100%
ca-certificates-2018 | 143 KB | ########################################################################################### | 100%
certifi-2018.11.29 | 145 KB | ########################################################################################### | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
pymeep-parallelにおいてmaterials.pyを使用したいときは、次のように1.7.0-py36_nomklh39e3cac_1にアップデートしてください。
(pmp) maru@tube60:~$ conda update -c chogan -c conda-forge pymeep-parallel
Solving environment: done
## Package Plan ##
environment location: /home/students/maru/anaconda3/envs/pmp
added / updated specs:
- pymeep-parallel
The following packages will be downloaded:
package | build
---------------------------|-----------------
libgcc-ng-7.3.0 | hdf63c60_0 6.1 MB conda-forge
pymeep-parallel-1.7.0 |py36_nomklh39e3cac_1 3.3 MB chogan
------------------------------------------------------------
Total: 9.4 MB
The following NEW packages will be INSTALLED:
libgdsii: 0.2.dev-0 chogan
nomkl: 3.0-0
The following packages will be UPDATED:
ca-certificates: 2018.8.24-ha4d7672_0 conda-forge --> 2018.11.29-ha4d7672_0 conda-forge
certifi: 2018.8.24-py36_1 conda-forge --> 2018.11.29-py36_1000 conda-forge
libctl: 4.1.2-0 chogan --> 4.1.4-1 chogan
libgcc-ng: 7.2.0-hdf63c60_3 conda-forge --> 7.3.0-hdf63c60_0 conda-forge
mpb: 1.6.2-4 chogan --> 1.7.0-nomkl_1 chogan [nomkl]
openssl: 1.0.2p-h470a237_0 conda-forge --> 1.0.2p-h14c3975_1002 conda-forge
pymeep-parallel: 1.5-py36_nomklh39e3cac_4 chogan [nomkl] --> 1.7.0-py36_nomklh39e3cac_1 chogan [nomkl]
Proceed ([y]/n)? y
Downloading and Extracting Packages
libgcc-ng-7.3.0 | 6.1 MB |
################################################################################################################################################################################################## | 100%
pymeep-parallel-1.7. | 3.3 MB |
##################################################################################################################################################### ############################################# | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
***error on line 497 of h5file.cpp: error opening HDF5 output file [#s9480bc6]
以下のようなエラーが出た場合は、計算を行っているディレクトリおよびファイルの書き込み・変更を許可してください。
Working in 3D dimensions.
Computational cell is 20 x 20 x 20 with resolution 10
block, center = (0,0,0)
size (1,3,1)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon diagonal = (1,1,1)
time for set_epsilon = 11.7491 s
-----------
creating output file "./hogehoge.h5"...
on time step 1 (time=0.05), 14.4708 s/step
HDF5-DIAG: Error detected in HDF5 (1.10.2) MPI-process 0:
#000: H5F.c line 445 in H5Fcreate(): unable to create file
major: File accessibilty
minor: Unable to open file
meep: error on line 497 of h5file.cpp: error opening HDF5 output file
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 1
***ImportError: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.14' not found (required by /home/students/maru/anaconda3/envs/mp/lib/python3.6/site-packages/meep/_meep.so) [#zcd4847a]
flexではGLIBCのバージョンが2.14よりも古いため次のエラーによりPyMeepを実行できません。
ImportError: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.14' not found (required by /home/students/maru/anaconda3/envs/mp/lib/python3.6/site-packages/meep/_meep.so)
tube70に移動してPyMeepを実行してください。
***Internal MPI error! Cannot read from remote process [#i054b206]
仙台高専からflexにsshしてmpiを使用すると次のエラーが出ます。
(pmp) maeda@tube60:/abinitio2/maeda/pymeep/Rectangle/drude/Omote$ mpirun -np 5 python after-Rec-drude-Au.py 45
Working in 3D dimensions.
Computational cell is 20 x 20 x 20 with resolution 10
block, center = (0,0,0)
size (1,3,1)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon diagonal = (1,1,1)
time for set_epsilon = 5.05642 s
lorentzian susceptibility: frequency=3.47141, gamma=2.01155
lorentzian susceptibility: frequency=2.39466, gamma=0.701702
lorentzian susceptibility: frequency=0.66944, gamma=0.278261
lorentzian susceptibility: frequency=0.33472, gamma=0.19438
drude susceptibility: frequency=1e-10, gamma=0.0427474
-----------
creating output file "./after-Rec-drude-Au-ex-surface-45deg.h5"...
creating output file "./after-Rec-drude-Au-ey-surface-45deg.h5"...
creating output file "./after-Rec-drude-Au-ez-inner-45deg.h5"...
creating output file "./after-Rec-drude-Au-ez-outer-45deg.h5"...
Fatal error in PMPI_Waitall: Other MPI error, error stack:
PMPI_Waitall(405)...............: MPI_Waitall(count=92, req_array=0x9adaad0, status_array=0xc9e3670) failed
MPIR_Waitall_impl(221)..........: fail failed
PMPIDI_CH3I_Progress(623).......: fail failed
pkt_RTS_handler(317)............: fail failed
do_cts(662).....................: fail failed
MPID_nem_lmt_dcp_start_recv(302): fail failed
dcp_recv(165)...................: Internal MPI error! Cannot read from remote process
Two workarounds have been identified for this issue:
1) Enable ptrace for non-root users with:
echo 0 | sudo tee /proc/sys/kernel/yama/ptrace_scope
2) Or, use:
I_MPI_SHM_LMT=shm
Fatal error in PMPI_Waitall: Other MPI error, error stack:
PMPI_Waitall(405)...............: MPI_Waitall(count=85, req_array=0xa328890, status_array=0xd21c7f0) failed
MPIR_Waitall_impl(221)..........: fail failed
PMPIDI_CH3I_Progress(658).......: fail failed
MPID_nem_handle_pkt(1439).......: fail failed
pkt_RTS_handler(317)............: fail failed
do_cts(662).....................: fail failed
MPID_nem_lmt_dcp_start_recv(302): fail failed
dcp_recv(165)...................: Internal MPI error! Cannot read from remote process
Two workarounds have been identified for this issue:
1) Enable ptrace for non-root users with:
echo 0 | sudo tee /proc/sys/kernel/yama/ptrace_scope
2) Or, use:
I_MPI_SHM_LMT=shm
Fatal error in PMPI_Waitall: Other MPI error, error stack:
PMPI_Waitall(405)...............: MPI_Waitall(count=85, req_array=0x9c2b8f0, status_array=0xcc31ed0) failed
MPIR_Waitall_impl(221)..........: fail failed
PMPIDI_CH3I_Progress(658).......: fail failed
MPID_nem_handle_pkt(1439).......: fail failed
pkt_RTS_handler(317)............: fail failed
do_cts(662).....................: fail failed
MPID_nem_lmt_dcp_start_recv(302): fail failed
dcp_recv(165)...................: Internal MPI error! Cannot read from remote process
Two workarounds have been identified for this issue:
1) Enable ptrace for non-root users with:
echo 0 | sudo tee /proc/sys/kernel/yama/ptrace_scope
2) Or, use:
I_MPI_SHM_LMT=shm
エラーメッセージで2つめのワークアラウンドとして提示されている
export I_MPI_SHM_LMT=shm
を実行するとmpiが実行できるようになります。
***静電場はどうシミュレーションすればいいか[#u810bd8a]
波源の周波数を十分低くします。
さらに波源の立ち上げを緩やかにします。
波源の立ち上げを緩やかにするにはContinuousSrcのwidthというパラメータを変更してください。
***ブロッホの周期的境界条件の下でシミュレーションした結果をh5topngすると"h5topng error: data can have at most two dimensions (try specifying a slice)"というエラーが出る [#h2f15931]
Schemeを使用している場合、ブロッホの周期的境界条件の下でシミュレーションした結果をh5topngすると
h5topng error: data can have at most two dimensions (try specifying a slice)
というエラーが出ます。
これはブロッホの周期的境界条件の下でシミュレーションすると計算結果に実部と虚部の両方が含まれるため、実部と虚部のどちらを出力するか選択する必要があるからです。
オプション"-d"によって出力するデータを指定してください。
出力できるデータはh5lsで調べることができます。
例えば
h5ls output.h5
ez.i Dataset {101, 101, 133/Inf}
ez.r Dataset {101, 101, 133/Inf}
となっているときは、電場のz成分の虚部ez.iと実部ez.rのどちらかを選択できます。
電場のz成分の実部ez.rを出力するときはh5opngにオプション"-d ez.r"を設定してください。
***共鳴周波数付近でシミュレーションするとノイズが入る [#f2efc046]
共鳴周波数付近では、波源の立ち上げの電磁場に対する応答のノイズがなかなかPMLの外に放出されません。
ノイズを防ぐには波源の立ち上げを緩やかにしてください。
例えばPyMeepを使用している場合はmp.ContinuousSourceのオプションwidthを設定してください。
***PMLと物体はどの程度離すべきか [#m1cd9b33]
PMLと物体は1/4波長、離してください。
PMLと物体が近接していると、PMLに対する電磁波の入射角が大きくなるためPMLの吸収率が減少してしまいます。
***なぜYee格子を使うのか [#ka49f596]
+Divergence free.
+Physical boundary condition are naturally satisfied.
+Elegant arrangement to approximate curl equations.
The proof of the divergence free nature of the Yee grid is derived by in Taflove's book on FDTD.By "satisfying natural boundary conditions" I mean that if you have a material interface meandering through your grid, you will not have to do anything special to simulate it accurately. If it were not for the Yee grid, you would observe crazy things like 800% reflection from a surface for some polarizations.The only remaining problem really is the staircasing of a curved surface due to the Cartesian grid. There are fixes for this like dielectric smoothing, unstructured grid, etc.
https://www.youtube.com/watch?v=hv5lIx4u8mY&t=566s
***fluxを取る面と波源を重ねるとfluxの値がおかしくなる。 [#b3ea0973]
fluxを取る面と波源を重ねないでください。
***ドルーデ伝導の媒質中で市松模様の定在場ができる。 [#nc8a54dc]
ドルーデ伝導の媒質中で市松模様の非物理的な共振モードが発生することがあります。
resolutionを上げると解決することが多いです。
メーリングリストに[https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg04502.html 議論]があります。
***電磁場の強度が減衰しないまま残った後に、無限大に発散する。 [#a65bf6fc]
PMLと分散性媒質を近すぎて、PMLと分散性媒質の電磁場が結合していることが原因のことがあります。
（グラフェンに光を斜めに入射してプラズモンを励起するシミュレーションの時におきました）
メーリングリストに[https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg04453.html 議論]があります。
***Harminv が 共鳴モードを見つけてくれない [#r5251f75]
[https://meep.readthedocs.io/en/latest/FAQ/#harminv-is-unable-to-find-the-resonant-modes-of-my-structure MEEPのFAQ]を日本語訳しただけです。
Harminv が 共鳴モードを見つけてくれないのは次の6つの理由が考えられます。
(1) シミュレーション時間が十分でなく、減衰率が小さいため、 Harminv に十分なデータを渡せていない。
(2) Harminvをafter_sourcesの中に入れていない。(sim.run_k_pointsから呼び出している場合は問題ないと思います)
(3) Harminvがデータを集める点がちょうど共鳴モードの定在波の節になっている。
(4) シミュレーションが不安定で電磁場が無限大に発散している。
(5) 減衰率が大きすぎる（50周期以内に減衰するモードをHarminvは無視する）。
バンドのある部分がごっそり途切れているときは、この点を疑ってください。ちょうど途切れている直前の点のQ値が50よりも十分大きいか確かめてください。Qが50台だと、Harminvが無視している可能性大です。散乱時間を長くするなどしてQ値を大きくして計算結果が変わらないか確認してください。
グラフェンナノリボンのプラズモンの分散関係を計算するときにこの点には嵌りました。
(6) PML が エバネッセント場と重なっていて、エバネッセント場を吸収し、共鳴モードが減衰してしまっている。
*** Q値 が 負 の モードが存在する。 [#vc2e9a60]
Q値が負ということは、指数関数的に強度が大きくなるモードであるということであり、非物理的です。
経験上、Harminvに渡すデータ量が足りないのが原因です。
until_after_sourceを大きくするとQ値が負のモードは消えることが多いです。
*** convert: UnableToOpenConfigureFile `delegates.xml' @ [#g9a5b2e3]
convert: UnableToOpenConfigureFile `delegates.xml' @
warning/configure.c/GetConfigureOptions/706.
convert: NoDecodeDelegateForThisImageFormat `PNG' @
error/constitute.c/ReadImage/501.
convert: NoImagesDefined `ey_in_graphene.gif' @
error/convert.c/ConvertImageCommand/3212.
login後にsource ~/.bash_profileをすると出現する。
loginし直すと消える。
*** エバネッセント波の発生している場所で、存在しないはずの電荷のプロットが現れる。 [#lcbd41a6]
エッジプラズモンの電荷分布をプロットする際、ナノリボンの端と接する真空側で存在しないはずの電荷が現れてしまいました。
電荷がナノリボンの端の数ピクセルに局在し、減衰長も数ピクセルという状況です。
たぶん、resolutionが足りないせいで、差分法の誤差が大きくなっているのだと思います。
エバネッセント光が指数関数的に減衰するため、1/r^kで減衰する電磁場に比べて誤差が現れやすいと感じています（全反射のシミュレーションの際も、界面のエバネッセント波の部分に存在しないはずの電荷が現れてしまったことがあります）。
低次元物質の端に局在する伝搬型プラズモンは特に注意が必要だと思います。
resolutionを上げることで改善できると思っています。調査中です。
*** paralell-meepにおいて、ディレクトリやファイルの作成・書き込みが複数回行われてしまう。 [#d2797f39]
ディレクトリやファイルの作成・書き込みは
if mp.am_master():
でラップをしないと、並列したプロセッサ数だけ同じioが行われてしまう。
***meep: invalid boundary absorbers for this grid_volume [#q3f4c827]
meep: invalid boundary absorbers for this grid_volume
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=1
PMLのサイズがセルのサイズより大きいと発生するようです（例えばx方向のcellの長さが3しかないのにPMLの厚さを2とすると、2*2>3なのでエラーが発生する）。
設定を間違えているはずなので修正しましょう。
PMLの厚さを減らすか、セルのサイズを大きくするといいと思います。
***エバネッセント波のPoynting vectorの積分が0にならない。 [#z5ac6365]
PoyntingVectorの積を取っている電場と磁場の時刻が同一であるかどうか確認してください。
時間領域をΔtで差分化しているとき、電場と磁場はそれぞれE(nΔt), H(nΔt+Δ/2)が計算されています（つまり同一時刻でない！！）。
E, Hの位相差がπ/2からπ/2+Δtになってしまい、
sin(ωt)cos(ωt+Δt) = {1/2}sin(2ωt+Δt)-{1/2}sinΔt
ですから、時間について積分を行っても{1/2}sinΔtが消えずに残ってしまいますね。
[https://meep.readthedocs.io/en/latest/Synchronizing_the_Magnetic_and_Electric_Fields/ 電場と磁場の時刻を同期]してください。
***負の屈折率の物質とPMLが重なっている部分で電磁場が発散する。 [#e9108eca]
PMLの実装上のバグ。
メーリングリストに[https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg06349.html 議論]あり。
***グラフェンの吸収率の計算値が入射角90付近で解析解から数十パーセントずれる。 [#u689c7c9]
周波数領域について離散フーリエ変換を実行し、poyhting vectorの積分を行う周波数の幅を小さくして、結果が変化しないか確認してみてください。
ガウス波源の中心周波数がf_cen、中心周波数における入射角がθ_cenだとします。
このとき、ガウス波に含まれる周波数fの波の入射角θは
θ = arcsin{(f_cen/f)sinθ_cen}
となり、sin関数は[0, π/2]で上に凸なので、θ_cenがπ/2に近づくほど、|θ-θ_cen|が大きくなってしまいます。
周波数幅をdfとすれば、θの取る範囲は
arcsin{(f_cen/(f+df))sinθ_cen} <
θ < arcsin{(f_cen/(f-df))sinθ_cen}
ですね。
***pymeepのmpbにmpiを使うと遅くなる [#gd81ea4c]
Pythonではmpiに対応していない[https://mpb.readthedocs.io/en/latest/Python_User_Interface/#mpb-with-mpi-parallelization 参考]。
Schemeは対応している。
*** On entry to ZGEBAL parameter number 3 had an illegal value [#dbb0b8ae]
発散した電磁場の値をharminvに渡すとこのエラーが発生する。
*** 'Simulation' object has no attribute 'run_k_point' [#j12afecf]
ノートPCの環境で発生しないが、tubeで実行すると発生する。原因不明。
***AttributeError: 'NoneType' object has no attribute 'get_epsilon' [#d662d316]
Must call ModeSolver.init_params or any run function before calling get_epsilon.
***Harminvの分散関係にω=(定数)の分散関係が現れる。 [#nf150ab8]
http://flex.phys.tohoku.ac.jp/~maru/drive-open/harminv190911.png
周波数を上げすぎ。(青線がlight-line)
**参考になるサイト [#f39b395a]
-[https://core.ac.uk/download/pdf/38256115.pdf ターゲット特性による超高強度レーザー生成高速電子の制御に関する研究(大阪大学リポジトリ]
-[http://kenji-ishikawa.cocolog-nifty.com/plasma/ プラズマナノテクノロジー(ブログ)]
-[http://www.alulab.org/group.htm Alu-sensei]
-[http://www.opt.ip.titech.ac.jp/research_j.html 梶川先生]
-[https://kx.lumerical.com/ lumerical knowledge exchange]
-[https://www.fujitsu.com/jp/solutions/business-technology/tc/fields/cae/poynting/index.html#optics 富士通Poyntingの解析事例]
*他の電磁界シミュレーションソフト [#qf1fc505]
株式会社EEMが無料で公開しているOpenFDTDとOpenMOMは[[OpenFDTD DLサイト:http://www.e-em.co.jp/OpenFDTD/index.html]]、[[OpenMOM DLサイト:http://www.e-em.co.jp/OpenMOM/index.html]]からダウンロードすることができます。
使い方はダウンロードサイトならびに研究室が所蔵している[[RFワールド No.39 2017年 08月号 :https://www.amazon.co.jp/RF%E3%83%AF%E3%83%BC%E3%83%AB%E3%83%89-No-39-2017%E5%B9%B4-08-%E6%9C%88%E5%8F%B7/dp/B07233YBBQ]]に記載されています。
書籍は丸岡が借りています（2018/6/12)。
=================まとめること
https://kb.lumerical.com/ref_sim_obj_planewave_edge.html
https://kb.lumerical.com/layout_analysis_diverging_simulations.html
![[PukiWiki]](image/c60.png "[PukiWiki]")
