MEEP (公開)
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
|
ログイン
]
開始行:
[2018/08 started by Maruoka and Maeda]
MEEPに関すること分からないことをmaru@flex.phys.tohoku.ac....
*目次 [#t1296280]
#contents
*MEEP [#g2bb328a]
**MEEPとは [#fbfd5425]
時間領域差分法(FDTD法)による電磁場のシミュレーションを...
Simpetusによってオープンソースで開発されており、フリーで...
商用ソフトではLumerical(FDTD)やCOMSOL(FEM)、富士通のPoyti...
例えば、富士通のPoyntingは[https://pr.fujitsu.com/jp/news...
フリーのFDTDソフトとしてはMEEPがほぼ唯一の選択肢だと思い...
オープンソースであることやCUIベースであることから、使用者...
**PyMeepのインストールと実行 [#fc543fa9]
Linux環境でPythonを用いてMEEPが使えるようになることを目標...
MEEPはLinuxとMacでしか動作しないため、Windowsを使用する場...
WSLのインストール方法は[[こちら:http://www.atmarkit.co.jp...
WSL上ではLinuxと同じ手順で環境が構築できます。~
***インストール [#p4cd7dd9]
まずはAnacondaを用いてPythonの実行環境を構築します。([[参...
[[Anacondaのサイト:https://www.anaconda.com/download/]]に...
インストーラーはBourne Shell(sh)のシェルスクリプトで書か...
$ bash [ダウンロードフォルダ]/Anaconda3-5.2.0-Linux-x86_...
次にPyMeepをインストールします。([[参考:https://meep.read...
以下のコマンドを実行してください。
$ conda create -n mp -c chogan -c conda-forge pymeep
これでPyMeepのインストールは完了です。
次のエラーが出るときは、~/anaconda3/binにパスを通してくだ...
bash: conda: command not found
tcshを使用している場合は次のコマンドを用いると、~/anacond...
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(齋藤グループのワークステーションの名前)ではP...
#!/bin/sh
_CONDA_ROOT="/home/students/maru/anaconda3"
\. "$_CONDA_ROOT/etc/profile.d/conda.sh" || return $?
_conda_activate "$@"
というエラーが出ます。
これはログインシェルがtcshであるにも関わらず、shのシェル...
tcshからはtcshスクリプトしかsource出来ません。
bash
と打って、ログインシェルをbashに変更してからPyMeepの仮想...
実行はhogehoge.pyの部分を任意のファイルに置き換えてくださ...
(mp) $ python hogehoge.py
サンプルソースはGithub等で入手できるので[[このソースコー...
PyMeepの仮想環境を無効化するには以下のコマンドを実行して...
(mp) $ source deactivate
**MEEPの並列化 [#k1f4cae9]
***インストール([https://meep.readthedocs.io/en/latest/In...
並列処理の可能なPyMeepをインストールするにはAnacondaのイ...
conda create -n pmp -c conda-forge pymeep=*=mpi_mpich_*
!注意! 以前は以下のコマンドでのインストールが標準でしたが...
conda create -n pmp -c chogan -c conda-forge pymeep-para...
***実行 [#j3e26633]
$ source activate pmp
にて並列処理をサポートしたPyMeepの仮想環境をアクティベー...
$ mpirun -np 4 python <script_name>.py
を実行してください。
anacondaの一部のファイルがbashで書かれているため、tcshか...
エラーの出ているbashで書かれたスクリプトをtcshに書き直す...
オプションの" -np 4 "は4プロセッサで並列処理することを指...
並列に計算するプロセッサの数は適切な値に変更してください。
***別の並列化 [#e415ae6a]
独立なプログラムを並列に計算することで計算時間を削減する...
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(\begin{align}&\nabla\cdot\left(\varepsilon_0\var...
次のように変形します。
#mathjax(\begin{align}&(a\nabla)\cdot \left(1+i\frac{(\fr...
&mathjax{a};はシミュレーションで使用する単位長さです。
&mathjax{c};は光速です。
以下のように&mathjax{x',y',z',t',ω',E',H',σ'};を定めます。
#mathjax(\begin{align}&x'=\frac{x}{a},y'=\frac{y}{a},z'=\...
このとき、Maxwell方程式は次のようになります。
#mathjax(\begin{align}&\nabla'\cdot\left(1+i\frac{\sigma'...
この形のMaxwell方程式をMEEPで使っています。
このとき、真空の誘電率、透磁率はどちらも1です。
また、単位長さは&mathjax{a};[m]です。
MEEP中の単位時間&mathjax{a/c};に真空中の光は単位長さaだけ...
この光の速度をSI単位系に直すと、&mathjax{a / (a/c) = c}; ...
MEEPの単位系では光は単位時間&mathjax{c/a};に単位長さ&math...
&mathjax{x',y',z',t',ω',σ'};はすべて無次元です。
金属が存在しない場合(&mathjax{σ=0};)は、&mathjax{a};の具...
つまり、一回のシミュレーションで任意の単位長さaに対するシ...
一方、金属が存在する場合(σ!=0)では、MEEPのシミュレーショ...
電場E'と磁場H'のMEEP中での単位は等しくなっており、平面波...
確かめてみましょう。平面波の場合、インピーダンス&mathjax{...
角振動数ω'(MEEP単位系)で 比誘電率を &mathjax{ε=ε_{re} +...
#mathjax( \epsilon_{r} = \epsilon_{re})
#mathjax( \sigma' = \omega' \frac{\epsilon_{im}}{ \epsilo...
とすればいいことも上の式から分かります。
**波源 [#n24a815a]
***TM波の作り方 [#j61c2744]
MEEPでは、電流を波源の代わりに用いています。
例えば、Ex偏光の点光源を設置するときには、代わりにJ = (Jx...
入射角Θの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...
***TE波の作り方 [#z315d617]
電流ではなく磁流を使う。
**左円偏光 [#u95702e4]
次のコードで左回りの円偏光になります。
sourcesクラスのamplitudeによって、y方向の電磁場の位相をx...
(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...
(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方向に直線偏光した平面波を重...
(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/c...
誘電体界面に入射角48度で電磁場が入射したときのシミュレー...
全反射が入射角arcsin(sqrt(1.25/2.25))=48.2度で起きるので...
z軸に垂直な境界にはPMLを、x軸に垂直な境界にはブロッホの周...
計算領域の大きさはx,y,z方向それぞれについて10,0,10です。
電磁場はy方向依存性を持たないので、y方向の計算領域の大き...
入射角をΘとすると、波数ベクトルkはk=k_0(sinΘ,0,cosΘ)です。
波源はx軸に平行な線として定義するので、位相のx依存性とし...
Schemeではsourceクラスのオプションであるamp-funcを設定す...
また、計算する領域内(今回の場合は大きさ10×0×10の長方形)...
x軸と垂直な計算領域の境界の点rについてE(r)exp(i k・R) = E...
磁場についても同様です。
点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 (vect...
(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...
(define k (vector3-scale (* 2 pi fcen (sqrt epsilon1))
(unit-vector3 kdir))) ; k with c...
(set! pml-layers (list (make pml (thickness 2)(direction...
(set! k-point (vector3* (* fcen (sqrt epsilon1) ) (vecto...
(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) (fwi...
(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...
(to-appended "ex" (at-every 0.1 output-efield...
)
Ezのプロット
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
Exのプロット
屈折角が90度近いので、透過光のExの大きさは小さくなってい...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
***光源を2つ並べる [#d178c899]
[https://github.com/stevengj/meep/blob/master/scheme/exam...
入射角45度の場合に限ってのみ、2つの波原で斜め入射の平面波...
計算領域がxz平面上に乗っている場合を例として考えます。
境界条件は四方全てPMLとします。
左側の境界と下側の境界にそれぞれ線状の波原を設置します。
波原の位相のx,z依存性は波数ベクトルk=(k_0sinθ, 0, k_0cosθ...
Schemeではsourceクラスのオプションであるamp-funcを設定す...
Ezのプロット (dfを0.1に変更しています)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
光源が一つしかない場合のEzのプロット(dfを0.1に変更してい...
このプロットから一様な振幅の大きさを持つ斜め入射の電磁場...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
***光源を斜めにたくさん並べる [#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-si...
(define-param beam-waist 2.5) ; beam sigma (gaussian bea...
(define-param rotation-angle (* (/ 22.5 360) 2 pi)) ; De...
(define-param source-points 60) ; should be a big number
(define-param source-size (* 4 beam-waist)) ; should be ...
(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 = -sour...
(set! src_list (append src_list
(list (make source
(src (make gaussian-src
(wavelength 1)
(width 3)
))
(amplitude (exp (- 0 (/...
(component Ez)
(center (* r_0 (sin rot...
))
))
)
(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...
)
plot of Ez
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
[http://article.gmane.org/gmane.comp.science.electromagne...
[http://kenji-ishikawa.cocolog-nifty.com/plasma/2010/01/p...
**Otto geometry [#y9a02bff]
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(wav...
center=mp.Vector3(x=x_, y=y_...
component=mp.Ex,
amplitude=coe
))
sources.append(mp.Source(src=mp.ContinuousSource(wav...
center=mp.Vector3(x=x_, y=y...
component=mp.Ey,
amplitude=coe
))
material_=mp.Medium(
epsilon=1,
E_susceptibilities = [mp.DrudeSusceptibility(frequen...
)
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.Volum...
until=20
)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
**負の屈折率 [#t4f04921]
屈折率nはn=√mu√eps(mu,epsはそれぞれ比透磁率と比誘電率)...
このとき波の位相速度は通常の位相速度と進行方向が逆向きに...
シミュレーションの際に、負の屈折率を持つ媒質を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(wav...
center=mp.Vector3(x=x_, y=y_...
component=mp.Ez,
amplitude=coe
))
material_=mp.Medium(
epsilon=1,
E_susceptibilities = [mp.DrudeSusceptibility(frequen...
mu=1,
H_susceptibilities = [mp.DrudeSusceptibility(frequen...
)
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.Volum...
until=30
)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/n...
**反射率、吸収率、透過率を計算したい [#cee00bb6]
物体のあるときとないとき両方についてガウシアンソースを入...
例えば透過率ならば透過光のポインティングベクトルの積分/入...
試していませんが、frequency solverを使って単一周波数のCon...
ここでは一つ目の方法を紹介します。
***TM波の反射・透過率 [#e20bae0b]
誘電体にTM波が斜め入射したときの反射率・透過率を計算しま...
z>0の領域に誘電体(誘電率9)が存在し、入射面がxz平面の場...
入射角Θのとき、入射光の波数ベクトルが \vec{k} = k_0(sinΘ,...
この波源の波数ベクトルがブロッホの周期境界条件の波数kにな...
電磁場の周期性がないz方向も周期境界条件が設定されてしまう...
z方向にはPMLを設置すれば疑似的に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-siz...
(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 ...
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen)...
(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...
(display-fluxes trans refl)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/TM/out.png
***TE波の反射・透過率 [#lde2976a]
; sz : computational domain size in z direct...
; theta-rad : incident angle in radian
; fcen : pulse center frequency (500 nm)
; df : pu...
; nfreq : nu...
; unit length: 100 micro meter
(set-param! resolution 200)
(define-param sz 3)
(set! geometry-lattice (make lattice(size no-size no-siz...
(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 ...
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen)...
(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...
;(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
***誘電体とグラフェンの多層構造の反射・透過・吸収率 [#n9c...
アンドープのグラフェンを誘電体で挟むと特殊な条件で吸収率...
(S. A. Nulli, M. S. Ukhtary, R. Saito, Appl. Phys. Lett. ...
(J. Blumberg, M. S. Ukhtary, R. Saito, Phys. Rev. Applied...
MEEPで再現できます。
グラフェンを二種類の誘電体A, Bで挟んで例えばABABGBABAのよ...
ここでGはグラフェンの層です。
それぞれの誘電体の厚みは垂直入射したときの光路長が1/4波長...
誘電体Aの層の数を2nとします。
このとき、誘電体AとBの誘電率の比がε_A/ε_B={2/(Zσ)}^{1/n}...
ここで、Zは真空のインピーダンス、σはグラフェンの電気伝導...
シミュレーションではグラフェンの伝導率σはσ=πe^2/(2h)とし...
平面波が誘電体とグラフェンの層に対して垂直入射するとき電...
従ってシミュレーションは一次元で行うことができます(マクス...
; a : unit length in the program
; e : elementary charge
; hbar : reduced planck constant
; epsilon_vaccum : permittivity of vaccum
; sigma : the conductivity of graphene in t...
; sigma_d : the conductivity of graphene in t...
; d ...
; dg : the thickness of graphene in this...
; fcen ...
; df ...
; nfreq ...
; isEmpty : Whether graphene exist or not.If ...
; sz : the size of computatinal domain i...
; Z_0 : the impedance of the vaccum
; s : the number of repetitions of the ...
; alpha : When the ratio of index of dielec...
; 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-siz...
(define-param isEmpty "DO NOT SET VALUE DIRECTLY")
(set! dimensions 1)
(if (string? isEmpty) ((print "Set the parameter \"isEm...
(if (not isEmpty)
(set! geometry
(list
(make block
(center 0 0 0)
(size infinity infinity dg)
(material (make medium (epsilon 1)(D-conductivity s...
)
(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)...
(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...
電場分布は次のようになります(パルス幅や誘電率、入射光の...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/dielectric...
**電場ベクトルをプロットしたい。 [#fee3ee23]
全反射における電場ベクトルの分布
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/e...
シミュレーションのためのコード
(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 t...
(set! theta (* deg (/ pi 180)))
(set! resolution 40);define the resolution for simulatio...
;define the lattice
(set! geometry-lattice (make lattice (size (+ dx (* 2 tp...
;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 (+ tpm...
(material (make dielectric (epsilon 1))))
(make block (center 0 (* -0.25 dy))(size infinity (+ tp...
(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) (direc...
;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 (cent...
(to-appended "ex"(at-every 0.05 (in-volume (volume (cent...
)
電場ベクトルをプロットするには、HDF5ファイルから直接プロ...
pythonのmatplotlibのquiver関数を用いるといいと思います。
**電荷分布を見たい。 [#l0ec2ad8]
MEEPには電荷分布を計算する機能がないため、自作する必要が...
ガウスの法則を差分化すればよいです。
前述の全反射直前のTM波(ブロッホ周期境界条件)のシミュレー...
電荷分布(+電場ベクトル)のプロット
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/...
誘電体界面で誘電分極が起きていることと波源で電荷の移動が...
**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/p...
コードは次の通りです。
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字や螺旋を形成...
***アルキメデスの螺旋 [#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...
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.Vec...
until=1
)
z=0平面における誘電率のプロットは次の通りです(PMLは除い...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/a...
**物性値のライブラリ [#ue28f25c]
***Schemeの場合 [#z7a9ca6a]
物性値のライブラリが公開されています。
[https://github.com/NanoComp/meep/blob/master/scheme/mate...
ライブラリを使用するには、プログラムの最初に
(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, tube...
詳しい使い方は https://meep.readthedocs.io/en/latest/Mate...
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 2...
[GCC 4.8.2 20140120 (Red Hat 4.8.2-15)] on linux
Type "help", "copyright", "credits" or "license" for mor...
>>> import meep
Using MPI version 3.1, 1 processes
>>> meep.__file__
'/home/students/maru/anaconda3/envs/pmp/lib/python3.6/si...
この場合、materials.pyのアドレスは/home/students/maru/ana...
このmaterials.pyの書き換えは(pmp)の仮想環境でだけ有効であ...
((mp)の仮想環境でmaterials.pyの単位長さを変える場合は(mp)...
**ユーザー定義のドルーデ伝導 [#df7ff749]
ドルーデ伝導で表される次のような比誘電率を用いる方法を示...
(&mathjax{\varepsilon_\infty};は周波数が十分高いところで...
#mathjax(\varepsilon_r = \varepsilon_{\infty} -\frac{{f_p...
MEEPのドルーデ伝導はMEEPの単位系の&mathjax{f_p'};, &mathj...
#mathjax(\begin{align} \varepsilon_r &= \varepsilon_{\inf...
この式から、SI単位系の&mathjax{f_p};, &mathjax{γ}; とMEEP...
#mathjax(f_p'=\frac{a}{c}f_p)
#mathjax(γ'=\frac{a}{c}γ)
上の関係式から求めた&mathjax{f_p'};, &mathjax{γ'};と&math...
user_defined_material = mp.Medium(epsilon=ε_∞, E_suscept...
とすると、その&mathjax{f_p'};, &mathjax{γ'};と&mathjax{\v...
この物質を使用する際は、geometryのオプション materials = ...
#mathjax(\begin{align} \varepsilon_r &= \varepsilon_\inft...
&mathjax{\sigma_0=f_p(2\pi)^2\varepsilon_0\tau};を代入す...
#mathjax(\varepsilon_r = \varepsilon_\infty - \frac{f_p^2...
散乱時間τとγの関係は
#mathjax(γ=\frac{1}{2\piτ})
となる。よって散逸がない場合&mathjax{\gamma=0};とする。
伝導率と&mathjax{f_p'};, &mathjax{γ'};と&mathjax{\varepsi...
#mathjax(\begin{align}\sigma &=-\left(\frac{i}{2\pi f\var...
**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/q...
**電磁場のh5ファイルをmp4に変換するためのスクリプト [#x47...
以下のスクリプトを.bash_profileに貼り付けて使用すると、h5...
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が区切り回数です。
**印加電場と実際の電場の時間変化を折れ線グラフでプロット ...
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 do...
+遠方界の領域を円の縁で指定するなら円の半径は波長に比べて...
+遠方界の領域を正方形の縁で指定するなら正方形の1辺の長さ...
アンテナは計算領域の中心に配置してください。
フラックス面はPMLより内側に設置してください。
パルス波源を使う場合は、電磁波がPMLの外に完全に発散するま...
以下に点光源を疑似的ダイポールアンテナと見立てた場合の遠...
[[該当ドキュメント>https://meep.readthedocs.io/en/latest/...
[[シミュレーションのコード>http://flex.phys.tohoku.ac.jp/...
このシミュレーションでは円で遠方界を計算します。
計算領域はxy平面の二次元です。
Near2FarRegion()はFluxRegion()と全く同じ(関数名が違うだ...
nearfield_boxに遠方場に必要なデータが入っているので,
get_farfieldで引数にnearfield_boxを引数にとり,遠方界の座...
[[プロット用のコード>http://flex.phys.tohoku.ac.jp/~maeda...
実行コマンド群(上からシミュレーション、データの整形、プ...
$ python antenna-radiation.py | tee source_Jz_farfields....
$ grep farfield: source_Jz_farfields.out | cut -d , -f2-...
$ python stable-antenna-radiation-plot.py
実行からプロットまでを1つのプログラムで実現したものが
[[こちら>http://flex.phys.tohoku.ac.jp/~maeda/php/antenna...
出力結果~
この結果は"antenna-radiation.py"のsrc_cmptにmp.Ex, mp.Ey,...
上から順にこれらの画像を作成することができます。~
それぞれ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.i...
まずはmpb(MIT Photonic Bands)を次のようにimportしてくださ...
from meep import mpb
それぞれのk点でいくつのバンドを求めるかをnum_bandsで指定...
例えば8つのバンドを求めたいときは
num_bands = 8
としてください。
計算したいk点をk_pointsリストに入れてください。例えば(0,0...
k_points=[mp.Vector3(0,0,0)]
k点の単位は&mathjax{2\pi/a};であることに注意してください(...
harminvを用いる場合のkも同様です。従って、SI単位系に直す...
単位胞の大きさは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()の代わりにそ...
計算結果はnumpyの二次元配列 ms.all_freqs に格納されます。
i行目にはi番目(0-origin)のk点での計算結果が格納されます。
格納されている情報は次のようになっています。
ms.all_freqs[i] = [1つ目のバンドの周波数, 2つ目のバンドの...
***誘電率のマッピング [#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)
単位胞を立方体(二次元の時は正方形)に取っていないと画像は...
そのような場合には次のように変換した誘電率データ(converte...
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_Analysi...
geometry_lattice = mp.Lattice(size=mp.Vector3(1, 1),
basis1=mp.Vector3(math.sqrt...
basis2=mp.Vector3(math.sqrt...
と設定してください。
basis1, basis2 が実空間の基本ベクトル。
k点を設定するときは単位が基本逆格子ベクトルになっているこ...
***誘電体多層構造のバンド構造 [#sdf7ab3e]
[http://ab-initio.mit.edu/book/photonic-crystals-book.pdf...
異なる誘電体の2種類の層を交互に積層させた多層構造のバンド...
このような構造はレーリー卿によって1887年に最初に考えられ...
層に垂直に伝搬する(z方向とする)光のみを考えます。
GaAs(ヒ化ガリウム、誘電率=13)のバルクのバンド構造は次の...
これはk=(ω/c)*sqrt(ε)のlight lineに相当します。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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, ...
geometry = [mp.Block(size=mp.Vector3(x=mp.inf, y=mp.inf,...
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",al...
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種類の層を交互に並べるとバンドギャップがブ...
ブリルアンゾーンの端は波長2a(2格子分)に対応し、単位胞の...
一般に低周波数のモードは高誘電率に局在したモード、高周波...
(詳しくは[http://ab-initio.mit.edu/book/photonic-crystal...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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, ...
geometry = [mp.Block(size=mp.Vector3(x=mp.inf, y=mp.inf,...
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",al...
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/m...
誘電率の分布は次の図のようになっています(黒い部分がGaAs...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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, ...
geometry = [mp.Block(size=mp.Vector3(x=0.5, y=mp.inf, z=...
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=Tr...
dfields = []
ms.run_tm(mpb.output_at_kpoint(k_points[t], mpb.fix_d...
md = mpb.MPBData(rectify=True, resolution=64, periods...
converted_dfields=[]
for f in dfields:
f=f[...,0,2] # f[i][j][k][l] = 位置(i,j)の第l成分...
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.im...
ax.imshow(np.real(f).T, interpolation='spline...
else:
ax.imshow(np.imag(f).T, interpolation='spline...
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(conver...
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/m...
sim.run_k_pointsの引数及び返り値の波数の単位は[&mathjax{2...
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),c...
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), c...
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)...
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],c...
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を用いた銀の表面プラズモンのバンド構造の計算 [#...
場所: /abinitio/maru/cpa/spl_bulk/TM/dispersion/main.py
厚さは200nm(十分に厚く無視できると仮定)。幅は無限。単位...
プラズマ周波数&mathjax{\omega_p=2\pi\times(c/a)\times0.72...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
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),...
]
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), c...
center=mp.Vector3(x=0, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[mp.PML(2, direction=...
sources=s,
resolution=resolution
)
k_interp = 20
k_point_list = mp.interpolate(k_interp, [mp.Vector3(y=0)...
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)+"...
f.close()
- confinement length
- propagation length
harminv から exp(-δt) の δ が周波数の虚数成分として求まる...
まず、分散関係を多項式近似してから、微分を行うことで群速...
多項式近似は numpy の [https://docs.scipy.org/doc/numpy/r...
微分は同じく numpy の [https://docs.scipy.org/doc/numpy/r...
https://kx.lumerical.com/t/can-i-simulate-spp-propagation...
***MIM( Metal-Insulator-Metal )構造中のプラズモンのバンド...
E.N. Economou, Phy. Rev. Vol.182, 539-554 (1969)のMIM構造...
MIM(metal-insulator-metal)構造とは絶縁体薄膜を金属で挟ん...
プラズモンの分散関係がlight-lineと交差するため、特殊な構...
下図の×がMEEPによって計算したバンド構造です。
紫色の部分は解析解から計算したバンド構造です。
絶縁体薄膜の両面に励起されたプラズモンのエバネッセント波...
バンドを周波数の低い方から順にバンド1、バンド2、バンド...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/M...
下図では、バンド1のk_x = 1における層に垂直な電場成分をプ...
横軸が位置、縦軸が時間です。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/M...
下図では、バンド2のk_x = 1における層に垂直な電場成分をプ...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/M...
電場分布は金属・誘電体界面に局在し、誘電体に対してバンド...
誘電体に対して非対称なプラズモンを励起するには、波源を誘...
この点についてメーリングリストで[https://www.mail-archive...
例えば、波源を誘電体の中心に設置すると、バンド2は得られ...
バンド2を得るには波源を誘電体の中心からわずかにずらした...
バンド3は金属中の横波の電磁波に対する分散関係ω^2=ω_p^2+(...
このため、バンド3は金属中の横波の電磁波分布に類似した、...
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**...
return min(abs((a+b)+(a-b)*cmath.exp(-0.1*2*math.pi*c...
# 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.Me...
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), c...
center=mp.Vector3(x=0.5, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[mp.PML(1, direction=...
sources=s,
resolution=resolution)
k_interp = 19
k_point_list = mp.interpolate(k_interp, [mp.Vector3(y=0)...
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],...
ax.scatter([k_point_list[i].y],[res[i][j].real],c...
ax.plot([0,1],
[plasma_freq/cmath.sqrt(2),plasma_freq/cmath.sqrt(2)],co...
#ax.fill_between([0,1],[0,1],[1,1],facecolor="gray",alph...
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 )構造中のプラズモンのバ...
E.N. Economou, Phy. Rev. Vol.182, 539-554 (1969)のIMI構造...
IMI(insulator-metal-insulator)構造とは金属薄膜を絶縁体で...
下図の×がMEEPによって計算したバンド構造です。
紫色の部分は解析解から計算したバンド構造です。
金属薄膜両面に励起されるプラズモンのエバネッセント波が重...
バンドを周波数の低い方から順にバンド1、バンド2と名付け...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/I...
下図では、バンド1のk_x = 1における層に垂直な電場成分をプ...
横軸が位置、縦軸が時間です。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/I...
下図では、バンド2のk_x = 1における層に垂直な電場成分をプ...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/I...
電場分布は金属・誘電体界面に局在し、誘電体に対してバンド...
誘電体および金属中の横波の分散関係ω=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**...
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*c...
# 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=...
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),c...
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), c...
center=mp.Vector3(x=0.1, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[mp.PML(1, direction=...
sources=s,
resolution=resolution)
k_interp = 19
k_point_list = mp.interpolate(k_interp, [mp.Vector3(y=0)...
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],...
ax.scatter([k_point_list[i].y],[res[i][j].real],c...
ax.plot([0,1], [plasma_freq/cmath.sqrt(2),plasma_freq/cm...
#ax.fill_between([0,1],[0,1],[1,1],facecolor="gray",alph...
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...
アルミナ(誘電率8.9)の円柱が正方格子状に並んでいるフォトニ...
円柱の半径rと正方格子の一辺の長さaがr=0.2aの関係にあると...
誘電率の分布は次のようになっています。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
バンド構造は次のようになります。
青がTMモード、赤がTEモードです。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
バンド構造の極大値、極小値(バンドギャップの大きさにおい...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
Γ, 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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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(e...
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...
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",a...
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",al...
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)...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
一般に、低周波数のバンドほど電場のエネルギーが高誘電率の...
TEモードでは、電束密度は円柱に対して垂直に走っており、必...
したがって、TEモードではTMモードと比べてバンドギャップが...
以下のコードを用いました。
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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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(e...
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=Tr...
dfields = []
ms.run_tm(mpb.output_at_kpoint(k_points[t], mpb.fix_d...
md = mpb.MPBData(rectify=True, resolution=64, periods...
converted_dfields=[]
for f in dfields:
f=f[...,0,2] # f[i][j][k][l] = 位置(i,j)の第l成分...
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.im...
ax.imshow(np.real(f).T, interpolation='spline...
else:
ax.imshow(np.imag(f).T, interpolation='spline...
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: Moldin...
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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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...
mp.Block(center=mp.Vector3(), size=mp.Vector3(x=sx, ...
]
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",...
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",a...
#================================================
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_phas...
elif(t==1):
fields.append(ms.get_dfield(band, bloch_phas...
if(t==0):
ms.run_te(mpb.output_at_kpoint(mp.Vector3(x=0.5,...
elif(t==1):
ms.run_tm(mpb.output_at_kpoint(mp.Vector3(x=0.5,...
md = mpb.MPBData(rectify=True, resolution=200, perio...
converted_dfields=[]
for f in fields:
f=f[...,0,2] # f[i][j][k][l] = 位置(i,j)の第l成...
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.i...
ax.imshow(np.real(f).T, interpolation='splin...
else:
ax.imshow(np.imag(f).T, interpolation='splin...
ax.axis("off")
if(t==0):
fig.savefig("te_hfield"+str(i+1), transparen...
elif(t==1):
fig.savefig("tm_dfield"+str(i+1), transparen...
ax.clear()
ax.imshow(mpb.MPBData(periods=3, resolution=200).convert...
ax.axis("off")
fig.savefig("epsilon.png",transparent=True)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
***triangular lattice of air columns [#q58923ea]
[http://ab-initio.mit.edu/book/ Photonic Crystals: Moldin...
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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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...
basis2=mp.Vector3(math.sqrt...
)
geometry = [
mp.Block(center=mp.Vector3(), size=mp.Vector3(x=sx, y...
mp.Cylinder(center=mp.Vector3(), radius=0.48, materia...
]
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",a...
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",al...
#================================================
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...
ax.axis("off")
fig.savefig("epsilon.png",transparent=True)
***状態密度(DOS) [#sa17a2a4]
https://www.mail-archive.com/mpb-discuss@ab-initio.mit.ed...
f=[0, df, 2df, 3df, ...]の各周波数についてDOSを求めて、プ...
周波数f[i]におけるDOSは次の式で近似できる。
#mathjax( \sum_j{ \frac{1}{\sqrt{\pi}} \exp\left(-\left(\...
**共振モードの計算 [#n9e267ed]
短パルスに対する応答を計算することで、調和振動する共振モ...
短パルスを用いるのは、短パルスが周波数領域で広がっている...
***円柱型の誘電体の共振モード [#g98983fa]
MEEP公式の[https://meep.readthedocs.io/en/latest/Python_T...
円柱はz軸方向に無限に伸びているとします。
電場がz方向に偏向したモードを計算します。
共振モードの周波数が分からないため、中心周波数1,周波数幅0...
円柱の屈折率は3.4です。
円柱の周りは真空です。
ガウス波源が十分小さくなった後に残っている電磁場がHarminv...
次のrun(step_functions...,until_after_sources=time)関数は...
run(step_functions..., until_after_sources=time)
また、各タイムステップ毎にstep_functions...を呼び出します。
until_after_sourcesに渡した数(time)だけ、全ての波源がturn...
次のafter_sources関数はstep_functionの一種です。ラップし...
after_sources(step_functions...)
次のHarminv関数は、点ptにおける電磁場成分c(例えばmp.Ez)を...
f(t)=Σa_j exp(-i*omega_j*t+δt)
の形に分解します。
Harminv(c, pt, fcen, df)
Harminvは次のような標準出力を行います。
harminv0:, frequency, imag. freq., Q, |amp|, amplitude, ...
harminv0:, 0.11678193348699654, -0.0006938113872636076, ...
harminv0:, 0.14635420376193028, -0.00019625091272617196,...
harminv0:, 0.17476005848070048, -4.962385886158849e-05, ...
frequency=omega[2πc], imag.freq=δ, Q = omega/(-2δ), |a_n|...
得られた共振周波数での実際の電磁場の分布を調べたいときは...
TMモードの共振モードとして次の電場(Ez)分布が見つかりま...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
また、TEモードの共振モードとして次の磁場(Hz)分布が見つ...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
z方向には周期境界条件を用いているため、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...
上のシミュレーションでは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),...
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...
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=d...
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/...
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])),...
ax.axis('off')
ims.append([im])
ani = animation.ArtistAnimation(fig,ims,interval=100)
ani.save("fcen_"+str(fcen)+"_TE.gif",writer="imagemag...
** frequency-domain solver ( 単一周波数の光源に対する応答...
MEEPの[https://meep.readthedocs.io/en/latest/Python_Tutor...
frequency-domain solverでは単一周波数の光源に対する応答を...
FDTD法( time-domain solver )で単一周波数の光源に対する応...
波源を立ち上げる時に目的外の周波数の波が混ざってしまいま...
time-domain solverでは、この目的外の周波数の波が励起した...
frequency-domain solverでは原理的に単一周波数の光源に対す...
frequency-domain solverの定式化については[http://ab-initi...
+波源 ContinuousSrc を使う。言い換えれば、GaussianSrc を...
+ mp.Simulation() に force_complex_fields=True を渡すこと...
+シミュレーションの前に sim.init_sim() を呼び出して初期化...
+sim.solve_cm(tl, maxiters, L) によって、frequency-domain...
Lorenz-Drude伝導は frequency-domain solver では使用できな...
**円筒座標におけるシミュレーション [#u3c69825]
[https://meep.readthedocs.io/en/latest/Python_Tutorials/R...
mp.Simulation に dimensions = mp.CYLINDRICAL を設定すると...
系が回転対称性を持つ時、電磁場分布は exp(mΦ) の依存性を持...
mp.Simulation のオプション m で、この動径方向の依存性を設...
cell は 直交座標系のときと同様に mp.Vector3(x, y, z) で定...
平面の吸収層として設計された PML を曲げて使っているため、...
cell の半径を大きくすると、PMLが平面に近づき反射は小さく...
**Q&A [#y2d911b6]
***一つのファイル内で複数回シミュレーションをすると計算結...
1つのファイル内で複数回シミュレーションをするときは、各シ...
(reset-meep)を呼び出さない場合、正しい計算結果が得られま...
***(入射波+反射波)のフラックスが入射波のみの時より、小...
フラックスは座標軸の正の方向を正に取っているためです(ポ...
正の方向に入射波が進行し、負の方向に反射波が進行するとき...
***dielectricよりmediumでindexを指定した方が計算が早い(...
理由は不明です。
***resolutionをsetしても変化しない [#g280b7dd]
resolutionはset! geometry-latticeの前に置かないとデフォル...
***変数vを使っていないのにunbounded variable vというエラ...
(define lambda 5)とするとunbounded variable vとしてエラー...
Schemeではlambdaを用いて関数を定義するからです(予約語に...
***meep: Could not determine normal direction for given g...
fluxを3次元で設定していることが原因です。
fluxは2次元の面または1次元の点で設定しなければなりません。
***実行時にGUILE_WARN_DEPRECATEDという警告が出る [#ma3f1c...
放っておいても問題ありませんが、.loginに以下を記述すると...
setenv GUILE_WARN_DEPRECATED no
***unbound variable: とエラーが出て、しかも : 以降が空白...
「unbound variable: 」とエラーが出てるにも関わらず、unbou...
全角空白を削除すればよいです。
***コマンドラインから変数を渡しているのに無視される。 [#d...
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について、厚さをそれぞれ別々にしたい [#z...
X, Y, Z方向のPMLをそれぞれ厚さa, b, cにするには
(set! pml-layers (list (make pml (thickness a)(direction...
(set! pml-layers (list (make pml (thickness b)(direction...
(set! pml-layers (list (make pml (thickness c)(direction...
としてください。
***負の透過率、反射率、吸収率を示す [#j49576fc]
負の透過率、反射率、吸収率は電磁波がPMLの外に発散しきって...
特にstop-when-fields-decayed関数を使用しているときは注意...
(run-sources+ T
(stop-when-fields-decayed dt Ex (vector3 x...
とすると、T秒間シミュレーションした後に、(x,y,z)におけるd...
しかしこの関数は、長時間共鳴するような構造の場合は、まだ...
電磁波が発散しきっているかどうかの判断は、h5topngなどで電...
***伝導度の理論値を複数のローレンツの重ね合わせで表したい...
下のリンクの方法で取り組んでいました。
理論値を複数のローレンツの重ね合わせで合わせるのはよくな...
https://www.mail-archive.com/meep-discuss@ab-initio.mit.e...
***インターフェイスはPythonとSchemeどちらを使うべきか [#m...
マニュアルのトップページにSchemeがobsoleteでPythonに変え...
***PMLの厚さはどのくらいにすべきか[#l8406a8a]
電磁波がPMLに垂直入射するときはPMLの厚さを半波長にしてく...
以下、MEEP作成者のPMLに関するドキュメント( http://math.mi...
With PML, however, the constant factor is very good to s...
(Increasing the resolution also increases the resolution...
斜め入射のときはPMLの吸収率が減少するので、厚みを増やす必...
反射率は入射角をXとすると、ある定数δを用いて
e^{-2*k*d*cos(X)*δ}
と書けます。
dはPMLの厚さ、kは波数。
詳細は[https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumb...
***tubeとflexにおいてmaterials.pyが使えない[#m4142e79]
PyMeepのバージョンが1.5だと次のエラーによりmaterials.pyが...
tube60:~$ python
Python 3.6.7 | packaged by conda-forge | (default, Nov 2...
[GCC 7.3.0] on linux
Type "help", "copyright", "credits" or "license" for mor...
>>> 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にアップデ...
(mp) maru@flex:~$ conda update -c chogan -c conda-forge ...
Solving environment: done
## Package Plan ##
environment location: /home/students/maru/anaconda3/env...
added / updated specs:
- pymeep
The following packages will be downloaded:
package | build
---------------------------|-----------------
libctl-4.1.4 | 1 ...
libgdsii-0.2.dev | 0 ...
python-3.6.6 | hd21baee_1003 ...
openssl-1.0.2p | h14c3975_1002 ...
mpb-1.7.0 | nomkl_1 ...
pymeep-1.7.0 | py36_nomkl_1 ...
openblas-0.3.5 | h9ac9557_1000 ...
nomkl-3.0 | 0 ...
ca-certificates-2018.11.29 | ha4d7672_0 ...
certifi-2018.11.29 | py36_1000 ...
-----------------------------------------------------...
Total: ...
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 --> 20...
certifi: 2018.11.29-py36_0 --> 20...
libctl: 4.1.2-0 chogan --> 4....
mpb: 1.6.2-4 chogan --> 1....
openblas: 0.2.20-8 conda-forge --> 0....
pymeep: 1.5-py36_nomkl_2 chogan [nomkl...
The following packages will be DOWNGRADED:
openssl: 1.1.1a-h7b6447c_0 --> 1....
python: 3.6.8-h0371630_0 --> 3....
Proceed ([y]/n)? y
Downloading and Extracting Packages
libctl-4.1.4 | 97 KB | #####################...
libgdsii-0.2.dev | 800 KB | #####################...
python-3.6.6 | 29.0 MB | #####################...
openssl-1.0.2p | 3.1 MB | #####################...
mpb-1.7.0 | 74 KB | #####################...
pymeep-1.7.0 | 3.3 MB | #####################...
openblas-0.3.5 | 15.8 MB | #####################...
nomkl-3.0 | 48 KB | #####################...
ca-certificates-2018 | 143 KB | #####################...
certifi-2018.11.29 | 145 KB | #####################...
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
pymeep-parallelにおいてmaterials.pyを使用したいときは、次...
(pmp) maru@tube60:~$ conda update -c chogan -c conda-for...
Solving environment: done
## Package Plan ##
environment location: /home/students/maru/anaconda3/env...
added / updated specs:
- pymeep-parallel
The following packages will be downloaded:
package | build
---------------------------|-----------------
libgcc-ng-7.3.0 | hdf63c60_0 ...
pymeep-parallel-1.7.0 |py36_nomklh39e3cac_1 ...
-----------------------------------------------------...
Total: ...
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...
certifi: 2018.8.24-py36_1 conda-forge...
libctl: 4.1.2-0 chogan ...
libgcc-ng: 7.2.0-hdf63c60_3 conda-forge...
mpb: 1.6.2-4 chogan ...
openssl: 1.0.2p-h470a237_0 conda-forge...
pymeep-parallel: 1.5-py36_nomklh39e3cac_4 chogan ...
Proceed ([y]/n)? y
Downloading and Extracting Packages
libgcc-ng-7.3.0 | 6.1 MB |
########################################################...
pymeep-parallel-1.7. | 3.3 MB |
########################################################...
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
***error on line 497 of h5file.cpp: error opening HDF5 ou...
以下のようなエラーが出た場合は、計算を行っているディレク...
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 f...
major: File accessibilty
minor: Unable to open file
meep: error on line 497 of h5file.cpp: error opening HDF...
application called MPI_Abort(MPI_COMM_WORLD, 1) - proces...
***ImportError: /lib/x86_64-linux-gnu/libc.so.6: version ...
flexではGLIBCのバージョンが2.14よりも古いため次のエラーに...
ImportError: /lib/x86_64-linux-gnu/libc.so.6: version `G...
tube70に移動してPyMeepを実行してください。
***Internal MPI error! Cannot read from remote process [...
仙台高専からflexにsshしてmpiを使用すると次のエラーが出ま...
(pmp) maeda@tube60:/abinitio2/maeda/pymeep/Rectangle/dru...
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.01...
lorentzian susceptibility: frequency=2.39466, gamma=0.70...
lorentzian susceptibility: frequency=0.66944, gamma=0.27...
lorentzian susceptibility: frequency=0.33472, gamma=0.19...
drude susceptibility: frequency=1e-10, gamma=0.0427474
-----------
creating output file "./after-Rec-drude-Au-ex-surface-45...
creating output file "./after-Rec-drude-Au-ey-surface-45...
creating output file "./after-Rec-drude-Au-ez-inner-45de...
creating output file "./after-Rec-drude-Au-ez-outer-45de...
Fatal error in PMPI_Waitall: Other MPI error, error stack:
PMPI_Waitall(405)...............: MPI_Waitall(count=92, ...
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! C...
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, ...
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! C...
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, ...
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! C...
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とい...
***ブロッホの周期的境界条件の下でシミュレーションした結果...
Schemeを使用している場合、ブロッホの周期的境界条件の下で...
h5topng error: data can have at most two dimensions (try...
というエラーが出ます。
これはブロッホの周期的境界条件の下でシミュレーションする...
オプション"-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にオプション"-...
***共鳴周波数付近でシミュレーションするとノイズが入る [#f...
共鳴周波数付近では、波源の立ち上げの電磁場に対する応答の...
ノイズを防ぐには波源の立ち上げを緩やかにしてください。
例えばPyMeepを使用している場合はmp.ContinuousSourceのオプ...
***PMLと物体はどの程度離すべきか [#m1cd9b33]
PMLと物体は1/4波長、離してください。
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 i...
https://www.youtube.com/watch?v=hv5lIx4u8mY&t=566s
***fluxを取る面と波源を重ねるとfluxの値がおかしくなる。 [...
fluxを取る面と波源を重ねないでください。
***ドルーデ伝導の媒質中で市松模様の定在場ができる。 [#nc8...
ドルーデ伝導の媒質中で市松模様の非物理的な共振モードが発...
resolutionを上げると解決することが多いです。
メーリングリストに[https://www.mail-archive.com/meep-disc...
***電磁場の強度が減衰しないまま残った後に、無限大に発散す...
PMLと分散性媒質を近すぎて、PMLと分散性媒質の電磁場が結合...
(グラフェンに光を斜めに入射してプラズモンを励起するシミ...
メーリングリストに[https://www.mail-archive.com/meep-disc...
***Harminv が 共鳴モードを見つけてくれない [#r5251f75]
[https://meep.readthedocs.io/en/latest/FAQ/#harminv-is-un...
Harminv が 共鳴モードを見つけてくれないのは次の6つの理由...
(1) シミュレーション時間が十分でなく、減衰率が小さいため...
(2) Harminvをafter_sourcesの中に入れていない。(sim.run_k_...
(3) Harminvがデータを集める点がちょうど共鳴モードの定在波...
(4) シミュレーションが不安定で電磁場が無限大に発散してい...
(5) 減衰率が大きすぎる(50周期以内に減衰するモードをHarmi...
バンドのある部分がごっそり途切れているときは、この点を疑...
グラフェンナノリボンのプラズモンの分散関係を計算するとき...
(6) PML が エバネッセント場と重なっていて、エバネッセント...
*** Q値 が 負 の モードが存在する。 [#vc2e9a60]
Q値が負ということは、指数関数的に強度が大きくなるモードで...
経験上、Harminvに渡すデータ量が足りないのが原因です。
until_after_sourceを大きくするとQ値が負のモードは消えるこ...
*** convert: UnableToOpenConfigureFile `delegates.xml' @ ...
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し直すと消える。
*** エバネッセント波の発生している場所で、存在しないはず...
エッジプラズモンの電荷分布をプロットする際、ナノリボンの...
電荷がナノリボンの端の数ピクセルに局在し、減衰長も数ピク...
たぶん、resolutionが足りないせいで、差分法の誤差が大きく...
エバネッセント光が指数関数的に減衰するため、1/r^kで減衰す...
低次元物質の端に局在する伝搬型プラズモンは特に注意が必要...
resolutionを上げることで改善できると思っています。調査中...
*** paralell-meepにおいて、ディレクトリやファイルの作成・...
ディレクトリやファイルの作成・書き込みは
if mp.am_master():
でラップをしないと、並列したプロセッサ数だけ同じioが行わ...
***meep: invalid boundary absorbers for this grid_volume ...
meep: invalid boundary absorbers for this grid_volume
application called MPI_Abort(MPI_COMM_WORLD, 1) - proces...
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=1
PMLのサイズがセルのサイズより大きいと発生するようです(例...
設定を間違えているはずなので修正しましょう。
PMLの厚さを減らすか、セルのサイズを大きくするといいと思い...
***エバネッセント波のPoynting vectorの積分が0にならない。...
PoyntingVectorの積を取っている電場と磁場の時刻が同一であ...
時間領域をΔtで差分化しているとき、電場と磁場はそれぞれE(n...
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_...
***負の屈折率の物質とPMLが重なっている部分で電磁場が発散...
PMLの実装上のバグ。
メーリングリストに[https://www.mail-archive.com/meep-disc...
***グラフェンの吸収率の計算値が入射角90付近で解析解から数...
周波数領域について離散フーリエ変換を実行し、poyhting vect...
ガウス波源の中心周波数がf_cen、中心周波数における入射角が...
このとき、ガウス波に含まれる周波数fの波の入射角θは
θ = arcsin{(f_cen/f)sinθ_cen}
となり、sin関数は[0, π/2]で上に凸なので、θ_cenがπ/2に近づ...
周波数幅を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/...
Schemeは対応している。
*** On entry to ZGEBAL parameter number 3 had an illegal...
発散した電磁場の値をharminvに渡すとこのエラーが発生する。
*** 'Simulation' object has no attribute 'run_k_point' [#...
ノートPCの環境で発生しないが、tubeで実行すると発生する。...
***AttributeError: 'NoneType' object has no attribute 'ge...
Must call ModeSolver.init_params or any run function befo...
***Harminvの分散関係にω=(定数)の分散関係が現れる。 [#nf15...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/harminv190...
周波数を上げすぎ。(青線が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-technolog...
*他の電磁界シミュレーションソフト [#qf1fc505]
株式会社EEMが無料で公開しているOpenFDTDとOpenMOMは[[OpenF...
使い方はダウンロードサイトならびに研究室が所蔵している[[R...
書籍は丸岡が借りています(2018/6/12)。
=================まとめること
https://kb.lumerical.com/ref_sim_obj_planewave_edge.html
https://kb.lumerical.com/layout_analysis_diverging_simula...
終了行:
[2018/08 started by Maruoka and Maeda]
MEEPに関すること分からないことをmaru@flex.phys.tohoku.ac....
*目次 [#t1296280]
#contents
*MEEP [#g2bb328a]
**MEEPとは [#fbfd5425]
時間領域差分法(FDTD法)による電磁場のシミュレーションを...
Simpetusによってオープンソースで開発されており、フリーで...
商用ソフトではLumerical(FDTD)やCOMSOL(FEM)、富士通のPoyti...
例えば、富士通のPoyntingは[https://pr.fujitsu.com/jp/news...
フリーのFDTDソフトとしてはMEEPがほぼ唯一の選択肢だと思い...
オープンソースであることやCUIベースであることから、使用者...
**PyMeepのインストールと実行 [#fc543fa9]
Linux環境でPythonを用いてMEEPが使えるようになることを目標...
MEEPはLinuxとMacでしか動作しないため、Windowsを使用する場...
WSLのインストール方法は[[こちら:http://www.atmarkit.co.jp...
WSL上ではLinuxと同じ手順で環境が構築できます。~
***インストール [#p4cd7dd9]
まずはAnacondaを用いてPythonの実行環境を構築します。([[参...
[[Anacondaのサイト:https://www.anaconda.com/download/]]に...
インストーラーはBourne Shell(sh)のシェルスクリプトで書か...
$ bash [ダウンロードフォルダ]/Anaconda3-5.2.0-Linux-x86_...
次にPyMeepをインストールします。([[参考:https://meep.read...
以下のコマンドを実行してください。
$ conda create -n mp -c chogan -c conda-forge pymeep
これでPyMeepのインストールは完了です。
次のエラーが出るときは、~/anaconda3/binにパスを通してくだ...
bash: conda: command not found
tcshを使用している場合は次のコマンドを用いると、~/anacond...
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(齋藤グループのワークステーションの名前)ではP...
#!/bin/sh
_CONDA_ROOT="/home/students/maru/anaconda3"
\. "$_CONDA_ROOT/etc/profile.d/conda.sh" || return $?
_conda_activate "$@"
というエラーが出ます。
これはログインシェルがtcshであるにも関わらず、shのシェル...
tcshからはtcshスクリプトしかsource出来ません。
bash
と打って、ログインシェルをbashに変更してからPyMeepの仮想...
実行はhogehoge.pyの部分を任意のファイルに置き換えてくださ...
(mp) $ python hogehoge.py
サンプルソースはGithub等で入手できるので[[このソースコー...
PyMeepの仮想環境を無効化するには以下のコマンドを実行して...
(mp) $ source deactivate
**MEEPの並列化 [#k1f4cae9]
***インストール([https://meep.readthedocs.io/en/latest/In...
並列処理の可能なPyMeepをインストールするにはAnacondaのイ...
conda create -n pmp -c conda-forge pymeep=*=mpi_mpich_*
!注意! 以前は以下のコマンドでのインストールが標準でしたが...
conda create -n pmp -c chogan -c conda-forge pymeep-para...
***実行 [#j3e26633]
$ source activate pmp
にて並列処理をサポートしたPyMeepの仮想環境をアクティベー...
$ mpirun -np 4 python <script_name>.py
を実行してください。
anacondaの一部のファイルがbashで書かれているため、tcshか...
エラーの出ているbashで書かれたスクリプトをtcshに書き直す...
オプションの" -np 4 "は4プロセッサで並列処理することを指...
並列に計算するプロセッサの数は適切な値に変更してください。
***別の並列化 [#e415ae6a]
独立なプログラムを並列に計算することで計算時間を削減する...
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(\begin{align}&\nabla\cdot\left(\varepsilon_0\var...
次のように変形します。
#mathjax(\begin{align}&(a\nabla)\cdot \left(1+i\frac{(\fr...
&mathjax{a};はシミュレーションで使用する単位長さです。
&mathjax{c};は光速です。
以下のように&mathjax{x',y',z',t',ω',E',H',σ'};を定めます。
#mathjax(\begin{align}&x'=\frac{x}{a},y'=\frac{y}{a},z'=\...
このとき、Maxwell方程式は次のようになります。
#mathjax(\begin{align}&\nabla'\cdot\left(1+i\frac{\sigma'...
この形のMaxwell方程式をMEEPで使っています。
このとき、真空の誘電率、透磁率はどちらも1です。
また、単位長さは&mathjax{a};[m]です。
MEEP中の単位時間&mathjax{a/c};に真空中の光は単位長さaだけ...
この光の速度をSI単位系に直すと、&mathjax{a / (a/c) = c}; ...
MEEPの単位系では光は単位時間&mathjax{c/a};に単位長さ&math...
&mathjax{x',y',z',t',ω',σ'};はすべて無次元です。
金属が存在しない場合(&mathjax{σ=0};)は、&mathjax{a};の具...
つまり、一回のシミュレーションで任意の単位長さaに対するシ...
一方、金属が存在する場合(σ!=0)では、MEEPのシミュレーショ...
電場E'と磁場H'のMEEP中での単位は等しくなっており、平面波...
確かめてみましょう。平面波の場合、インピーダンス&mathjax{...
角振動数ω'(MEEP単位系)で 比誘電率を &mathjax{ε=ε_{re} +...
#mathjax( \epsilon_{r} = \epsilon_{re})
#mathjax( \sigma' = \omega' \frac{\epsilon_{im}}{ \epsilo...
とすればいいことも上の式から分かります。
**波源 [#n24a815a]
***TM波の作り方 [#j61c2744]
MEEPでは、電流を波源の代わりに用いています。
例えば、Ex偏光の点光源を設置するときには、代わりにJ = (Jx...
入射角Θの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...
***TE波の作り方 [#z315d617]
電流ではなく磁流を使う。
**左円偏光 [#u95702e4]
次のコードで左回りの円偏光になります。
sourcesクラスのamplitudeによって、y方向の電磁場の位相をx...
(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...
(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方向に直線偏光した平面波を重...
(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/c...
誘電体界面に入射角48度で電磁場が入射したときのシミュレー...
全反射が入射角arcsin(sqrt(1.25/2.25))=48.2度で起きるので...
z軸に垂直な境界にはPMLを、x軸に垂直な境界にはブロッホの周...
計算領域の大きさはx,y,z方向それぞれについて10,0,10です。
電磁場はy方向依存性を持たないので、y方向の計算領域の大き...
入射角をΘとすると、波数ベクトルkはk=k_0(sinΘ,0,cosΘ)です。
波源はx軸に平行な線として定義するので、位相のx依存性とし...
Schemeではsourceクラスのオプションであるamp-funcを設定す...
また、計算する領域内(今回の場合は大きさ10×0×10の長方形)...
x軸と垂直な計算領域の境界の点rについてE(r)exp(i k・R) = E...
磁場についても同様です。
点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 (vect...
(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...
(define k (vector3-scale (* 2 pi fcen (sqrt epsilon1))
(unit-vector3 kdir))) ; k with c...
(set! pml-layers (list (make pml (thickness 2)(direction...
(set! k-point (vector3* (* fcen (sqrt epsilon1) ) (vecto...
(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) (fwi...
(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...
(to-appended "ex" (at-every 0.1 output-efield...
)
Ezのプロット
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
Exのプロット
屈折角が90度近いので、透過光のExの大きさは小さくなってい...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
***光源を2つ並べる [#d178c899]
[https://github.com/stevengj/meep/blob/master/scheme/exam...
入射角45度の場合に限ってのみ、2つの波原で斜め入射の平面波...
計算領域がxz平面上に乗っている場合を例として考えます。
境界条件は四方全てPMLとします。
左側の境界と下側の境界にそれぞれ線状の波原を設置します。
波原の位相のx,z依存性は波数ベクトルk=(k_0sinθ, 0, k_0cosθ...
Schemeではsourceクラスのオプションであるamp-funcを設定す...
Ezのプロット (dfを0.1に変更しています)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
光源が一つしかない場合のEzのプロット(dfを0.1に変更してい...
このプロットから一様な振幅の大きさを持つ斜め入射の電磁場...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
***光源を斜めにたくさん並べる [#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-si...
(define-param beam-waist 2.5) ; beam sigma (gaussian bea...
(define-param rotation-angle (* (/ 22.5 360) 2 pi)) ; De...
(define-param source-points 60) ; should be a big number
(define-param source-size (* 4 beam-waist)) ; should be ...
(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 = -sour...
(set! src_list (append src_list
(list (make source
(src (make gaussian-src
(wavelength 1)
(width 3)
))
(amplitude (exp (- 0 (/...
(component Ez)
(center (* r_0 (sin rot...
))
))
)
(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...
)
plot of Ez
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
[http://article.gmane.org/gmane.comp.science.electromagne...
[http://kenji-ishikawa.cocolog-nifty.com/plasma/2010/01/p...
**Otto geometry [#y9a02bff]
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(wav...
center=mp.Vector3(x=x_, y=y_...
component=mp.Ex,
amplitude=coe
))
sources.append(mp.Source(src=mp.ContinuousSource(wav...
center=mp.Vector3(x=x_, y=y...
component=mp.Ey,
amplitude=coe
))
material_=mp.Medium(
epsilon=1,
E_susceptibilities = [mp.DrudeSusceptibility(frequen...
)
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.Volum...
until=20
)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
**負の屈折率 [#t4f04921]
屈折率nはn=√mu√eps(mu,epsはそれぞれ比透磁率と比誘電率)...
このとき波の位相速度は通常の位相速度と進行方向が逆向きに...
シミュレーションの際に、負の屈折率を持つ媒質を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(wav...
center=mp.Vector3(x=x_, y=y_...
component=mp.Ez,
amplitude=coe
))
material_=mp.Medium(
epsilon=1,
E_susceptibilities = [mp.DrudeSusceptibility(frequen...
mu=1,
H_susceptibilities = [mp.DrudeSusceptibility(frequen...
)
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.Volum...
until=30
)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/n...
**反射率、吸収率、透過率を計算したい [#cee00bb6]
物体のあるときとないとき両方についてガウシアンソースを入...
例えば透過率ならば透過光のポインティングベクトルの積分/入...
試していませんが、frequency solverを使って単一周波数のCon...
ここでは一つ目の方法を紹介します。
***TM波の反射・透過率 [#e20bae0b]
誘電体にTM波が斜め入射したときの反射率・透過率を計算しま...
z>0の領域に誘電体(誘電率9)が存在し、入射面がxz平面の場...
入射角Θのとき、入射光の波数ベクトルが \vec{k} = k_0(sinΘ,...
この波源の波数ベクトルがブロッホの周期境界条件の波数kにな...
電磁場の周期性がないz方向も周期境界条件が設定されてしまう...
z方向にはPMLを設置すれば疑似的に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-siz...
(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 ...
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen)...
(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...
(display-fluxes trans refl)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/TM/out.png
***TE波の反射・透過率 [#lde2976a]
; sz : computational domain size in z direct...
; theta-rad : incident angle in radian
; fcen : pulse center frequency (500 nm)
; df : pu...
; nfreq : nu...
; unit length: 100 micro meter
(set-param! resolution 200)
(define-param sz 3)
(set! geometry-lattice (make lattice(size no-size no-siz...
(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 ...
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen)...
(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...
;(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
***誘電体とグラフェンの多層構造の反射・透過・吸収率 [#n9c...
アンドープのグラフェンを誘電体で挟むと特殊な条件で吸収率...
(S. A. Nulli, M. S. Ukhtary, R. Saito, Appl. Phys. Lett. ...
(J. Blumberg, M. S. Ukhtary, R. Saito, Phys. Rev. Applied...
MEEPで再現できます。
グラフェンを二種類の誘電体A, Bで挟んで例えばABABGBABAのよ...
ここでGはグラフェンの層です。
それぞれの誘電体の厚みは垂直入射したときの光路長が1/4波長...
誘電体Aの層の数を2nとします。
このとき、誘電体AとBの誘電率の比がε_A/ε_B={2/(Zσ)}^{1/n}...
ここで、Zは真空のインピーダンス、σはグラフェンの電気伝導...
シミュレーションではグラフェンの伝導率σはσ=πe^2/(2h)とし...
平面波が誘電体とグラフェンの層に対して垂直入射するとき電...
従ってシミュレーションは一次元で行うことができます(マクス...
; a : unit length in the program
; e : elementary charge
; hbar : reduced planck constant
; epsilon_vaccum : permittivity of vaccum
; sigma : the conductivity of graphene in t...
; sigma_d : the conductivity of graphene in t...
; d ...
; dg : the thickness of graphene in this...
; fcen ...
; df ...
; nfreq ...
; isEmpty : Whether graphene exist or not.If ...
; sz : the size of computatinal domain i...
; Z_0 : the impedance of the vaccum
; s : the number of repetitions of the ...
; alpha : When the ratio of index of dielec...
; 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-siz...
(define-param isEmpty "DO NOT SET VALUE DIRECTLY")
(set! dimensions 1)
(if (string? isEmpty) ((print "Set the parameter \"isEm...
(if (not isEmpty)
(set! geometry
(list
(make block
(center 0 0 0)
(size infinity infinity dg)
(material (make medium (epsilon 1)(D-conductivity s...
)
(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)...
(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...
電場分布は次のようになります(パルス幅や誘電率、入射光の...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/dielectric...
**電場ベクトルをプロットしたい。 [#fee3ee23]
全反射における電場ベクトルの分布
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/e...
シミュレーションのためのコード
(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 t...
(set! theta (* deg (/ pi 180)))
(set! resolution 40);define the resolution for simulatio...
;define the lattice
(set! geometry-lattice (make lattice (size (+ dx (* 2 tp...
;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 (+ tpm...
(material (make dielectric (epsilon 1))))
(make block (center 0 (* -0.25 dy))(size infinity (+ tp...
(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) (direc...
;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 (cent...
(to-appended "ex"(at-every 0.05 (in-volume (volume (cent...
)
電場ベクトルをプロットするには、HDF5ファイルから直接プロ...
pythonのmatplotlibのquiver関数を用いるといいと思います。
**電荷分布を見たい。 [#l0ec2ad8]
MEEPには電荷分布を計算する機能がないため、自作する必要が...
ガウスの法則を差分化すればよいです。
前述の全反射直前のTM波(ブロッホ周期境界条件)のシミュレー...
電荷分布(+電場ベクトル)のプロット
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/o...
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/...
誘電体界面で誘電分極が起きていることと波源で電荷の移動が...
**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/p...
コードは次の通りです。
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字や螺旋を形成...
***アルキメデスの螺旋 [#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...
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.Vec...
until=1
)
z=0平面における誘電率のプロットは次の通りです(PMLは除い...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/a...
**物性値のライブラリ [#ue28f25c]
***Schemeの場合 [#z7a9ca6a]
物性値のライブラリが公開されています。
[https://github.com/NanoComp/meep/blob/master/scheme/mate...
ライブラリを使用するには、プログラムの最初に
(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, tube...
詳しい使い方は https://meep.readthedocs.io/en/latest/Mate...
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 2...
[GCC 4.8.2 20140120 (Red Hat 4.8.2-15)] on linux
Type "help", "copyright", "credits" or "license" for mor...
>>> import meep
Using MPI version 3.1, 1 processes
>>> meep.__file__
'/home/students/maru/anaconda3/envs/pmp/lib/python3.6/si...
この場合、materials.pyのアドレスは/home/students/maru/ana...
このmaterials.pyの書き換えは(pmp)の仮想環境でだけ有効であ...
((mp)の仮想環境でmaterials.pyの単位長さを変える場合は(mp)...
**ユーザー定義のドルーデ伝導 [#df7ff749]
ドルーデ伝導で表される次のような比誘電率を用いる方法を示...
(&mathjax{\varepsilon_\infty};は周波数が十分高いところで...
#mathjax(\varepsilon_r = \varepsilon_{\infty} -\frac{{f_p...
MEEPのドルーデ伝導はMEEPの単位系の&mathjax{f_p'};, &mathj...
#mathjax(\begin{align} \varepsilon_r &= \varepsilon_{\inf...
この式から、SI単位系の&mathjax{f_p};, &mathjax{γ}; とMEEP...
#mathjax(f_p'=\frac{a}{c}f_p)
#mathjax(γ'=\frac{a}{c}γ)
上の関係式から求めた&mathjax{f_p'};, &mathjax{γ'};と&math...
user_defined_material = mp.Medium(epsilon=ε_∞, E_suscept...
とすると、その&mathjax{f_p'};, &mathjax{γ'};と&mathjax{\v...
この物質を使用する際は、geometryのオプション materials = ...
#mathjax(\begin{align} \varepsilon_r &= \varepsilon_\inft...
&mathjax{\sigma_0=f_p(2\pi)^2\varepsilon_0\tau};を代入す...
#mathjax(\varepsilon_r = \varepsilon_\infty - \frac{f_p^2...
散乱時間τとγの関係は
#mathjax(γ=\frac{1}{2\piτ})
となる。よって散逸がない場合&mathjax{\gamma=0};とする。
伝導率と&mathjax{f_p'};, &mathjax{γ'};と&mathjax{\varepsi...
#mathjax(\begin{align}\sigma &=-\left(\frac{i}{2\pi f\var...
**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/q...
**電磁場のh5ファイルをmp4に変換するためのスクリプト [#x47...
以下のスクリプトを.bash_profileに貼り付けて使用すると、h5...
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が区切り回数です。
**印加電場と実際の電場の時間変化を折れ線グラフでプロット ...
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 do...
+遠方界の領域を円の縁で指定するなら円の半径は波長に比べて...
+遠方界の領域を正方形の縁で指定するなら正方形の1辺の長さ...
アンテナは計算領域の中心に配置してください。
フラックス面はPMLより内側に設置してください。
パルス波源を使う場合は、電磁波がPMLの外に完全に発散するま...
以下に点光源を疑似的ダイポールアンテナと見立てた場合の遠...
[[該当ドキュメント>https://meep.readthedocs.io/en/latest/...
[[シミュレーションのコード>http://flex.phys.tohoku.ac.jp/...
このシミュレーションでは円で遠方界を計算します。
計算領域はxy平面の二次元です。
Near2FarRegion()はFluxRegion()と全く同じ(関数名が違うだ...
nearfield_boxに遠方場に必要なデータが入っているので,
get_farfieldで引数にnearfield_boxを引数にとり,遠方界の座...
[[プロット用のコード>http://flex.phys.tohoku.ac.jp/~maeda...
実行コマンド群(上からシミュレーション、データの整形、プ...
$ python antenna-radiation.py | tee source_Jz_farfields....
$ grep farfield: source_Jz_farfields.out | cut -d , -f2-...
$ python stable-antenna-radiation-plot.py
実行からプロットまでを1つのプログラムで実現したものが
[[こちら>http://flex.phys.tohoku.ac.jp/~maeda/php/antenna...
出力結果~
この結果は"antenna-radiation.py"のsrc_cmptにmp.Ex, mp.Ey,...
上から順にこれらの画像を作成することができます。~
それぞれ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.i...
まずはmpb(MIT Photonic Bands)を次のようにimportしてくださ...
from meep import mpb
それぞれのk点でいくつのバンドを求めるかをnum_bandsで指定...
例えば8つのバンドを求めたいときは
num_bands = 8
としてください。
計算したいk点をk_pointsリストに入れてください。例えば(0,0...
k_points=[mp.Vector3(0,0,0)]
k点の単位は&mathjax{2\pi/a};であることに注意してください(...
harminvを用いる場合のkも同様です。従って、SI単位系に直す...
単位胞の大きさは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()の代わりにそ...
計算結果はnumpyの二次元配列 ms.all_freqs に格納されます。
i行目にはi番目(0-origin)のk点での計算結果が格納されます。
格納されている情報は次のようになっています。
ms.all_freqs[i] = [1つ目のバンドの周波数, 2つ目のバンドの...
***誘電率のマッピング [#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)
単位胞を立方体(二次元の時は正方形)に取っていないと画像は...
そのような場合には次のように変換した誘電率データ(converte...
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_Analysi...
geometry_lattice = mp.Lattice(size=mp.Vector3(1, 1),
basis1=mp.Vector3(math.sqrt...
basis2=mp.Vector3(math.sqrt...
と設定してください。
basis1, basis2 が実空間の基本ベクトル。
k点を設定するときは単位が基本逆格子ベクトルになっているこ...
***誘電体多層構造のバンド構造 [#sdf7ab3e]
[http://ab-initio.mit.edu/book/photonic-crystals-book.pdf...
異なる誘電体の2種類の層を交互に積層させた多層構造のバンド...
このような構造はレーリー卿によって1887年に最初に考えられ...
層に垂直に伝搬する(z方向とする)光のみを考えます。
GaAs(ヒ化ガリウム、誘電率=13)のバルクのバンド構造は次の...
これはk=(ω/c)*sqrt(ε)のlight lineに相当します。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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, ...
geometry = [mp.Block(size=mp.Vector3(x=mp.inf, y=mp.inf,...
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",al...
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種類の層を交互に並べるとバンドギャップがブ...
ブリルアンゾーンの端は波長2a(2格子分)に対応し、単位胞の...
一般に低周波数のモードは高誘電率に局在したモード、高周波...
(詳しくは[http://ab-initio.mit.edu/book/photonic-crystal...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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, ...
geometry = [mp.Block(size=mp.Vector3(x=mp.inf, y=mp.inf,...
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",al...
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/m...
誘電率の分布は次の図のようになっています(黒い部分がGaAs...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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, ...
geometry = [mp.Block(size=mp.Vector3(x=0.5, y=mp.inf, z=...
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=Tr...
dfields = []
ms.run_tm(mpb.output_at_kpoint(k_points[t], mpb.fix_d...
md = mpb.MPBData(rectify=True, resolution=64, periods...
converted_dfields=[]
for f in dfields:
f=f[...,0,2] # f[i][j][k][l] = 位置(i,j)の第l成分...
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.im...
ax.imshow(np.real(f).T, interpolation='spline...
else:
ax.imshow(np.imag(f).T, interpolation='spline...
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(conver...
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/m...
sim.run_k_pointsの引数及び返り値の波数の単位は[&mathjax{2...
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),c...
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), c...
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)...
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],c...
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を用いた銀の表面プラズモンのバンド構造の計算 [#...
場所: /abinitio/maru/cpa/spl_bulk/TM/dispersion/main.py
厚さは200nm(十分に厚く無視できると仮定)。幅は無限。単位...
プラズマ周波数&mathjax{\omega_p=2\pi\times(c/a)\times0.72...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
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),...
]
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), c...
center=mp.Vector3(x=0, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[mp.PML(2, direction=...
sources=s,
resolution=resolution
)
k_interp = 20
k_point_list = mp.interpolate(k_interp, [mp.Vector3(y=0)...
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)+"...
f.close()
- confinement length
- propagation length
harminv から exp(-δt) の δ が周波数の虚数成分として求まる...
まず、分散関係を多項式近似してから、微分を行うことで群速...
多項式近似は numpy の [https://docs.scipy.org/doc/numpy/r...
微分は同じく numpy の [https://docs.scipy.org/doc/numpy/r...
https://kx.lumerical.com/t/can-i-simulate-spp-propagation...
***MIM( Metal-Insulator-Metal )構造中のプラズモンのバンド...
E.N. Economou, Phy. Rev. Vol.182, 539-554 (1969)のMIM構造...
MIM(metal-insulator-metal)構造とは絶縁体薄膜を金属で挟ん...
プラズモンの分散関係がlight-lineと交差するため、特殊な構...
下図の×がMEEPによって計算したバンド構造です。
紫色の部分は解析解から計算したバンド構造です。
絶縁体薄膜の両面に励起されたプラズモンのエバネッセント波...
バンドを周波数の低い方から順にバンド1、バンド2、バンド...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/M...
下図では、バンド1のk_x = 1における層に垂直な電場成分をプ...
横軸が位置、縦軸が時間です。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/M...
下図では、バンド2のk_x = 1における層に垂直な電場成分をプ...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/M...
電場分布は金属・誘電体界面に局在し、誘電体に対してバンド...
誘電体に対して非対称なプラズモンを励起するには、波源を誘...
この点についてメーリングリストで[https://www.mail-archive...
例えば、波源を誘電体の中心に設置すると、バンド2は得られ...
バンド2を得るには波源を誘電体の中心からわずかにずらした...
バンド3は金属中の横波の電磁波に対する分散関係ω^2=ω_p^2+(...
このため、バンド3は金属中の横波の電磁波分布に類似した、...
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**...
return min(abs((a+b)+(a-b)*cmath.exp(-0.1*2*math.pi*c...
# 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.Me...
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), c...
center=mp.Vector3(x=0.5, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[mp.PML(1, direction=...
sources=s,
resolution=resolution)
k_interp = 19
k_point_list = mp.interpolate(k_interp, [mp.Vector3(y=0)...
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],...
ax.scatter([k_point_list[i].y],[res[i][j].real],c...
ax.plot([0,1],
[plasma_freq/cmath.sqrt(2),plasma_freq/cmath.sqrt(2)],co...
#ax.fill_between([0,1],[0,1],[1,1],facecolor="gray",alph...
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 )構造中のプラズモンのバ...
E.N. Economou, Phy. Rev. Vol.182, 539-554 (1969)のIMI構造...
IMI(insulator-metal-insulator)構造とは金属薄膜を絶縁体で...
下図の×がMEEPによって計算したバンド構造です。
紫色の部分は解析解から計算したバンド構造です。
金属薄膜両面に励起されるプラズモンのエバネッセント波が重...
バンドを周波数の低い方から順にバンド1、バンド2と名付け...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/I...
下図では、バンド1のk_x = 1における層に垂直な電場成分をプ...
横軸が位置、縦軸が時間です。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/I...
下図では、バンド2のk_x = 1における層に垂直な電場成分をプ...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/I...
電場分布は金属・誘電体界面に局在し、誘電体に対してバンド...
誘電体および金属中の横波の分散関係ω=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**...
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*c...
# 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=...
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),c...
s = [mp.Source(src=mp.GaussianSource(fcen, fwidth=df), c...
center=mp.Vector3(x=0.1, y=0, z=0))]
sim = mp.Simulation(cell_size=cell,
geometry=geometry,
boundary_layers=[mp.PML(1, direction=...
sources=s,
resolution=resolution)
k_interp = 19
k_point_list = mp.interpolate(k_interp, [mp.Vector3(y=0)...
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],...
ax.scatter([k_point_list[i].y],[res[i][j].real],c...
ax.plot([0,1], [plasma_freq/cmath.sqrt(2),plasma_freq/cm...
#ax.fill_between([0,1],[0,1],[1,1],facecolor="gray",alph...
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...
アルミナ(誘電率8.9)の円柱が正方格子状に並んでいるフォトニ...
円柱の半径rと正方格子の一辺の長さaがr=0.2aの関係にあると...
誘電率の分布は次のようになっています。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
バンド構造は次のようになります。
青がTMモード、赤がTEモードです。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
バンド構造の極大値、極小値(バンドギャップの大きさにおい...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
Γ, 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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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(e...
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...
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",a...
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",al...
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)...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
一般に、低周波数のバンドほど電場のエネルギーが高誘電率の...
TEモードでは、電束密度は円柱に対して垂直に走っており、必...
したがって、TEモードではTMモードと比べてバンドギャップが...
以下のコードを用いました。
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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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(e...
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=Tr...
dfields = []
ms.run_tm(mpb.output_at_kpoint(k_points[t], mpb.fix_d...
md = mpb.MPBData(rectify=True, resolution=64, periods...
converted_dfields=[]
for f in dfields:
f=f[...,0,2] # f[i][j][k][l] = 位置(i,j)の第l成分...
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.im...
ax.imshow(np.real(f).T, interpolation='spline...
else:
ax.imshow(np.imag(f).T, interpolation='spline...
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: Moldin...
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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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...
mp.Block(center=mp.Vector3(), size=mp.Vector3(x=sx, ...
]
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",...
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",a...
#================================================
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_phas...
elif(t==1):
fields.append(ms.get_dfield(band, bloch_phas...
if(t==0):
ms.run_te(mpb.output_at_kpoint(mp.Vector3(x=0.5,...
elif(t==1):
ms.run_tm(mpb.output_at_kpoint(mp.Vector3(x=0.5,...
md = mpb.MPBData(rectify=True, resolution=200, perio...
converted_dfields=[]
for f in fields:
f=f[...,0,2] # f[i][j][k][l] = 位置(i,j)の第l成...
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.i...
ax.imshow(np.real(f).T, interpolation='splin...
else:
ax.imshow(np.imag(f).T, interpolation='splin...
ax.axis("off")
if(t==0):
fig.savefig("te_hfield"+str(i+1), transparen...
elif(t==1):
fig.savefig("tm_dfield"+str(i+1), transparen...
ax.clear()
ax.imshow(mpb.MPBData(periods=3, resolution=200).convert...
ax.axis("off")
fig.savefig("epsilon.png",transparent=True)
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/s...
***triangular lattice of air columns [#q58923ea]
[http://ab-initio.mit.edu/book/ Photonic Crystals: Moldin...
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。...
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'が...
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['xtick.direction'] = 'in' # x軸の目盛線が...
plt.rcParams['ytick.direction'] = 'in' # y軸の目盛線が...
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...
basis2=mp.Vector3(math.sqrt...
)
geometry = [
mp.Block(center=mp.Vector3(), size=mp.Vector3(x=sx, y...
mp.Cylinder(center=mp.Vector3(), radius=0.48, materia...
]
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",a...
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",al...
#================================================
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...
ax.axis("off")
fig.savefig("epsilon.png",transparent=True)
***状態密度(DOS) [#sa17a2a4]
https://www.mail-archive.com/mpb-discuss@ab-initio.mit.ed...
f=[0, df, 2df, 3df, ...]の各周波数についてDOSを求めて、プ...
周波数f[i]におけるDOSは次の式で近似できる。
#mathjax( \sum_j{ \frac{1}{\sqrt{\pi}} \exp\left(-\left(\...
**共振モードの計算 [#n9e267ed]
短パルスに対する応答を計算することで、調和振動する共振モ...
短パルスを用いるのは、短パルスが周波数領域で広がっている...
***円柱型の誘電体の共振モード [#g98983fa]
MEEP公式の[https://meep.readthedocs.io/en/latest/Python_T...
円柱はz軸方向に無限に伸びているとします。
電場がz方向に偏向したモードを計算します。
共振モードの周波数が分からないため、中心周波数1,周波数幅0...
円柱の屈折率は3.4です。
円柱の周りは真空です。
ガウス波源が十分小さくなった後に残っている電磁場がHarminv...
次のrun(step_functions...,until_after_sources=time)関数は...
run(step_functions..., until_after_sources=time)
また、各タイムステップ毎にstep_functions...を呼び出します。
until_after_sourcesに渡した数(time)だけ、全ての波源がturn...
次のafter_sources関数はstep_functionの一種です。ラップし...
after_sources(step_functions...)
次のHarminv関数は、点ptにおける電磁場成分c(例えばmp.Ez)を...
f(t)=Σa_j exp(-i*omega_j*t+δt)
の形に分解します。
Harminv(c, pt, fcen, df)
Harminvは次のような標準出力を行います。
harminv0:, frequency, imag. freq., Q, |amp|, amplitude, ...
harminv0:, 0.11678193348699654, -0.0006938113872636076, ...
harminv0:, 0.14635420376193028, -0.00019625091272617196,...
harminv0:, 0.17476005848070048, -4.962385886158849e-05, ...
frequency=omega[2πc], imag.freq=δ, Q = omega/(-2δ), |a_n|...
得られた共振周波数での実際の電磁場の分布を調べたいときは...
TMモードの共振モードとして次の電場(Ez)分布が見つかりま...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
また、TEモードの共振モードとして次の磁場(Hz)分布が見つ...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/c...
z方向には周期境界条件を用いているため、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...
上のシミュレーションでは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),...
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...
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=d...
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/...
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])),...
ax.axis('off')
ims.append([im])
ani = animation.ArtistAnimation(fig,ims,interval=100)
ani.save("fcen_"+str(fcen)+"_TE.gif",writer="imagemag...
** frequency-domain solver ( 単一周波数の光源に対する応答...
MEEPの[https://meep.readthedocs.io/en/latest/Python_Tutor...
frequency-domain solverでは単一周波数の光源に対する応答を...
FDTD法( time-domain solver )で単一周波数の光源に対する応...
波源を立ち上げる時に目的外の周波数の波が混ざってしまいま...
time-domain solverでは、この目的外の周波数の波が励起した...
frequency-domain solverでは原理的に単一周波数の光源に対す...
frequency-domain solverの定式化については[http://ab-initi...
+波源 ContinuousSrc を使う。言い換えれば、GaussianSrc を...
+ mp.Simulation() に force_complex_fields=True を渡すこと...
+シミュレーションの前に sim.init_sim() を呼び出して初期化...
+sim.solve_cm(tl, maxiters, L) によって、frequency-domain...
Lorenz-Drude伝導は frequency-domain solver では使用できな...
**円筒座標におけるシミュレーション [#u3c69825]
[https://meep.readthedocs.io/en/latest/Python_Tutorials/R...
mp.Simulation に dimensions = mp.CYLINDRICAL を設定すると...
系が回転対称性を持つ時、電磁場分布は exp(mΦ) の依存性を持...
mp.Simulation のオプション m で、この動径方向の依存性を設...
cell は 直交座標系のときと同様に mp.Vector3(x, y, z) で定...
平面の吸収層として設計された PML を曲げて使っているため、...
cell の半径を大きくすると、PMLが平面に近づき反射は小さく...
**Q&A [#y2d911b6]
***一つのファイル内で複数回シミュレーションをすると計算結...
1つのファイル内で複数回シミュレーションをするときは、各シ...
(reset-meep)を呼び出さない場合、正しい計算結果が得られま...
***(入射波+反射波)のフラックスが入射波のみの時より、小...
フラックスは座標軸の正の方向を正に取っているためです(ポ...
正の方向に入射波が進行し、負の方向に反射波が進行するとき...
***dielectricよりmediumでindexを指定した方が計算が早い(...
理由は不明です。
***resolutionをsetしても変化しない [#g280b7dd]
resolutionはset! geometry-latticeの前に置かないとデフォル...
***変数vを使っていないのにunbounded variable vというエラ...
(define lambda 5)とするとunbounded variable vとしてエラー...
Schemeではlambdaを用いて関数を定義するからです(予約語に...
***meep: Could not determine normal direction for given g...
fluxを3次元で設定していることが原因です。
fluxは2次元の面または1次元の点で設定しなければなりません。
***実行時にGUILE_WARN_DEPRECATEDという警告が出る [#ma3f1c...
放っておいても問題ありませんが、.loginに以下を記述すると...
setenv GUILE_WARN_DEPRECATED no
***unbound variable: とエラーが出て、しかも : 以降が空白...
「unbound variable: 」とエラーが出てるにも関わらず、unbou...
全角空白を削除すればよいです。
***コマンドラインから変数を渡しているのに無視される。 [#d...
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について、厚さをそれぞれ別々にしたい [#z...
X, Y, Z方向のPMLをそれぞれ厚さa, b, cにするには
(set! pml-layers (list (make pml (thickness a)(direction...
(set! pml-layers (list (make pml (thickness b)(direction...
(set! pml-layers (list (make pml (thickness c)(direction...
としてください。
***負の透過率、反射率、吸収率を示す [#j49576fc]
負の透過率、反射率、吸収率は電磁波がPMLの外に発散しきって...
特にstop-when-fields-decayed関数を使用しているときは注意...
(run-sources+ T
(stop-when-fields-decayed dt Ex (vector3 x...
とすると、T秒間シミュレーションした後に、(x,y,z)におけるd...
しかしこの関数は、長時間共鳴するような構造の場合は、まだ...
電磁波が発散しきっているかどうかの判断は、h5topngなどで電...
***伝導度の理論値を複数のローレンツの重ね合わせで表したい...
下のリンクの方法で取り組んでいました。
理論値を複数のローレンツの重ね合わせで合わせるのはよくな...
https://www.mail-archive.com/meep-discuss@ab-initio.mit.e...
***インターフェイスはPythonとSchemeどちらを使うべきか [#m...
マニュアルのトップページにSchemeがobsoleteでPythonに変え...
***PMLの厚さはどのくらいにすべきか[#l8406a8a]
電磁波がPMLに垂直入射するときはPMLの厚さを半波長にしてく...
以下、MEEP作成者のPMLに関するドキュメント( http://math.mi...
With PML, however, the constant factor is very good to s...
(Increasing the resolution also increases the resolution...
斜め入射のときはPMLの吸収率が減少するので、厚みを増やす必...
反射率は入射角をXとすると、ある定数δを用いて
e^{-2*k*d*cos(X)*δ}
と書けます。
dはPMLの厚さ、kは波数。
詳細は[https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumb...
***tubeとflexにおいてmaterials.pyが使えない[#m4142e79]
PyMeepのバージョンが1.5だと次のエラーによりmaterials.pyが...
tube60:~$ python
Python 3.6.7 | packaged by conda-forge | (default, Nov 2...
[GCC 7.3.0] on linux
Type "help", "copyright", "credits" or "license" for mor...
>>> 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にアップデ...
(mp) maru@flex:~$ conda update -c chogan -c conda-forge ...
Solving environment: done
## Package Plan ##
environment location: /home/students/maru/anaconda3/env...
added / updated specs:
- pymeep
The following packages will be downloaded:
package | build
---------------------------|-----------------
libctl-4.1.4 | 1 ...
libgdsii-0.2.dev | 0 ...
python-3.6.6 | hd21baee_1003 ...
openssl-1.0.2p | h14c3975_1002 ...
mpb-1.7.0 | nomkl_1 ...
pymeep-1.7.0 | py36_nomkl_1 ...
openblas-0.3.5 | h9ac9557_1000 ...
nomkl-3.0 | 0 ...
ca-certificates-2018.11.29 | ha4d7672_0 ...
certifi-2018.11.29 | py36_1000 ...
-----------------------------------------------------...
Total: ...
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 --> 20...
certifi: 2018.11.29-py36_0 --> 20...
libctl: 4.1.2-0 chogan --> 4....
mpb: 1.6.2-4 chogan --> 1....
openblas: 0.2.20-8 conda-forge --> 0....
pymeep: 1.5-py36_nomkl_2 chogan [nomkl...
The following packages will be DOWNGRADED:
openssl: 1.1.1a-h7b6447c_0 --> 1....
python: 3.6.8-h0371630_0 --> 3....
Proceed ([y]/n)? y
Downloading and Extracting Packages
libctl-4.1.4 | 97 KB | #####################...
libgdsii-0.2.dev | 800 KB | #####################...
python-3.6.6 | 29.0 MB | #####################...
openssl-1.0.2p | 3.1 MB | #####################...
mpb-1.7.0 | 74 KB | #####################...
pymeep-1.7.0 | 3.3 MB | #####################...
openblas-0.3.5 | 15.8 MB | #####################...
nomkl-3.0 | 48 KB | #####################...
ca-certificates-2018 | 143 KB | #####################...
certifi-2018.11.29 | 145 KB | #####################...
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
pymeep-parallelにおいてmaterials.pyを使用したいときは、次...
(pmp) maru@tube60:~$ conda update -c chogan -c conda-for...
Solving environment: done
## Package Plan ##
environment location: /home/students/maru/anaconda3/env...
added / updated specs:
- pymeep-parallel
The following packages will be downloaded:
package | build
---------------------------|-----------------
libgcc-ng-7.3.0 | hdf63c60_0 ...
pymeep-parallel-1.7.0 |py36_nomklh39e3cac_1 ...
-----------------------------------------------------...
Total: ...
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...
certifi: 2018.8.24-py36_1 conda-forge...
libctl: 4.1.2-0 chogan ...
libgcc-ng: 7.2.0-hdf63c60_3 conda-forge...
mpb: 1.6.2-4 chogan ...
openssl: 1.0.2p-h470a237_0 conda-forge...
pymeep-parallel: 1.5-py36_nomklh39e3cac_4 chogan ...
Proceed ([y]/n)? y
Downloading and Extracting Packages
libgcc-ng-7.3.0 | 6.1 MB |
########################################################...
pymeep-parallel-1.7. | 3.3 MB |
########################################################...
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
***error on line 497 of h5file.cpp: error opening HDF5 ou...
以下のようなエラーが出た場合は、計算を行っているディレク...
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 f...
major: File accessibilty
minor: Unable to open file
meep: error on line 497 of h5file.cpp: error opening HDF...
application called MPI_Abort(MPI_COMM_WORLD, 1) - proces...
***ImportError: /lib/x86_64-linux-gnu/libc.so.6: version ...
flexではGLIBCのバージョンが2.14よりも古いため次のエラーに...
ImportError: /lib/x86_64-linux-gnu/libc.so.6: version `G...
tube70に移動してPyMeepを実行してください。
***Internal MPI error! Cannot read from remote process [...
仙台高専からflexにsshしてmpiを使用すると次のエラーが出ま...
(pmp) maeda@tube60:/abinitio2/maeda/pymeep/Rectangle/dru...
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.01...
lorentzian susceptibility: frequency=2.39466, gamma=0.70...
lorentzian susceptibility: frequency=0.66944, gamma=0.27...
lorentzian susceptibility: frequency=0.33472, gamma=0.19...
drude susceptibility: frequency=1e-10, gamma=0.0427474
-----------
creating output file "./after-Rec-drude-Au-ex-surface-45...
creating output file "./after-Rec-drude-Au-ey-surface-45...
creating output file "./after-Rec-drude-Au-ez-inner-45de...
creating output file "./after-Rec-drude-Au-ez-outer-45de...
Fatal error in PMPI_Waitall: Other MPI error, error stack:
PMPI_Waitall(405)...............: MPI_Waitall(count=92, ...
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! C...
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, ...
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! C...
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, ...
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! C...
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とい...
***ブロッホの周期的境界条件の下でシミュレーションした結果...
Schemeを使用している場合、ブロッホの周期的境界条件の下で...
h5topng error: data can have at most two dimensions (try...
というエラーが出ます。
これはブロッホの周期的境界条件の下でシミュレーションする...
オプション"-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にオプション"-...
***共鳴周波数付近でシミュレーションするとノイズが入る [#f...
共鳴周波数付近では、波源の立ち上げの電磁場に対する応答の...
ノイズを防ぐには波源の立ち上げを緩やかにしてください。
例えばPyMeepを使用している場合はmp.ContinuousSourceのオプ...
***PMLと物体はどの程度離すべきか [#m1cd9b33]
PMLと物体は1/4波長、離してください。
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 i...
https://www.youtube.com/watch?v=hv5lIx4u8mY&t=566s
***fluxを取る面と波源を重ねるとfluxの値がおかしくなる。 [...
fluxを取る面と波源を重ねないでください。
***ドルーデ伝導の媒質中で市松模様の定在場ができる。 [#nc8...
ドルーデ伝導の媒質中で市松模様の非物理的な共振モードが発...
resolutionを上げると解決することが多いです。
メーリングリストに[https://www.mail-archive.com/meep-disc...
***電磁場の強度が減衰しないまま残った後に、無限大に発散す...
PMLと分散性媒質を近すぎて、PMLと分散性媒質の電磁場が結合...
(グラフェンに光を斜めに入射してプラズモンを励起するシミ...
メーリングリストに[https://www.mail-archive.com/meep-disc...
***Harminv が 共鳴モードを見つけてくれない [#r5251f75]
[https://meep.readthedocs.io/en/latest/FAQ/#harminv-is-un...
Harminv が 共鳴モードを見つけてくれないのは次の6つの理由...
(1) シミュレーション時間が十分でなく、減衰率が小さいため...
(2) Harminvをafter_sourcesの中に入れていない。(sim.run_k_...
(3) Harminvがデータを集める点がちょうど共鳴モードの定在波...
(4) シミュレーションが不安定で電磁場が無限大に発散してい...
(5) 減衰率が大きすぎる(50周期以内に減衰するモードをHarmi...
バンドのある部分がごっそり途切れているときは、この点を疑...
グラフェンナノリボンのプラズモンの分散関係を計算するとき...
(6) PML が エバネッセント場と重なっていて、エバネッセント...
*** Q値 が 負 の モードが存在する。 [#vc2e9a60]
Q値が負ということは、指数関数的に強度が大きくなるモードで...
経験上、Harminvに渡すデータ量が足りないのが原因です。
until_after_sourceを大きくするとQ値が負のモードは消えるこ...
*** convert: UnableToOpenConfigureFile `delegates.xml' @ ...
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し直すと消える。
*** エバネッセント波の発生している場所で、存在しないはず...
エッジプラズモンの電荷分布をプロットする際、ナノリボンの...
電荷がナノリボンの端の数ピクセルに局在し、減衰長も数ピク...
たぶん、resolutionが足りないせいで、差分法の誤差が大きく...
エバネッセント光が指数関数的に減衰するため、1/r^kで減衰す...
低次元物質の端に局在する伝搬型プラズモンは特に注意が必要...
resolutionを上げることで改善できると思っています。調査中...
*** paralell-meepにおいて、ディレクトリやファイルの作成・...
ディレクトリやファイルの作成・書き込みは
if mp.am_master():
でラップをしないと、並列したプロセッサ数だけ同じioが行わ...
***meep: invalid boundary absorbers for this grid_volume ...
meep: invalid boundary absorbers for this grid_volume
application called MPI_Abort(MPI_COMM_WORLD, 1) - proces...
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=1
PMLのサイズがセルのサイズより大きいと発生するようです(例...
設定を間違えているはずなので修正しましょう。
PMLの厚さを減らすか、セルのサイズを大きくするといいと思い...
***エバネッセント波のPoynting vectorの積分が0にならない。...
PoyntingVectorの積を取っている電場と磁場の時刻が同一であ...
時間領域をΔtで差分化しているとき、電場と磁場はそれぞれE(n...
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_...
***負の屈折率の物質とPMLが重なっている部分で電磁場が発散...
PMLの実装上のバグ。
メーリングリストに[https://www.mail-archive.com/meep-disc...
***グラフェンの吸収率の計算値が入射角90付近で解析解から数...
周波数領域について離散フーリエ変換を実行し、poyhting vect...
ガウス波源の中心周波数がf_cen、中心周波数における入射角が...
このとき、ガウス波に含まれる周波数fの波の入射角θは
θ = arcsin{(f_cen/f)sinθ_cen}
となり、sin関数は[0, π/2]で上に凸なので、θ_cenがπ/2に近づ...
周波数幅を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/...
Schemeは対応している。
*** On entry to ZGEBAL parameter number 3 had an illegal...
発散した電磁場の値をharminvに渡すとこのエラーが発生する。
*** 'Simulation' object has no attribute 'run_k_point' [#...
ノートPCの環境で発生しないが、tubeで実行すると発生する。...
***AttributeError: 'NoneType' object has no attribute 'ge...
Must call ModeSolver.init_params or any run function befo...
***Harminvの分散関係にω=(定数)の分散関係が現れる。 [#nf15...
http://flex.phys.tohoku.ac.jp/~maru/drive-open/harminv190...
周波数を上げすぎ。(青線が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-technolog...
*他の電磁界シミュレーションソフト [#qf1fc505]
株式会社EEMが無料で公開しているOpenFDTDとOpenMOMは[[OpenF...
使い方はダウンロードサイトならびに研究室が所蔵している[[R...
書籍は丸岡が借りています(2018/6/12)。
=================まとめること
https://kb.lumerical.com/ref_sim_obj_planewave_edge.html
https://kb.lumerical.com/layout_analysis_diverging_simula...
ページ名: