- 履歴一覧
- 履歴 を表示
- MEEP (Open) は削除されています。
[2018/08 started by Maruoka and Maeda]
#contents
*MEEP [#g2bb328a]
**PythonでMEEP [#fc543fa9]
Linux環境でPythonを用いてMEEPを使えるようすることを目標とします。~
WindowsではMEEPがLinuxかMacでしか動作しないのでLinuxのWSL(Windows Subsystem for Linux)のインストールをしてください。
やり方は[[こちら:http://www.atmarkit.co.jp/ait/articles/1608/08/news039.html]]。~
WindowsでLinuxが使えるようになったらLinuxのと同じ手順で環境が構築できます。~
MacではAnacondaをコマンドラインインストーラーを使ってインストールすればLinuxの時と同様にできると思います。(未確認)
***インストール [#p4cd7dd9]
まずはAnacondaを用いてPythonの環境を構築します。([[参考:https://qiita.com/t2y/items/2a3eb58103e85d8064b6]])~
[[Anacondaのサイト:https://www.anaconda.com/download/]]にてPython3.x系のインストーラーをダウンロードしてください。~
インストーラーは.shのファイルなので以下のように実行してください。
$ bash [ダウンロードフォルダ]/Anaconda3-5.2.0-Linux-x86_64.sh
次にPyMeepのインストールをします。([[参考:https://meep.readthedocs.io/en/latest/Installation/]])~
以下のコマンドを実行してください。
$ conda create -n mp -c chogan -c conda-forge pymeep
これで必要なセットアップは終わりました。
***実行 [#padb0b5f]
以下のコマンドでアクティベートしてください。
$ source activate mp
"-bash: activate: No such file or directory" とエラーが出たときは
export PATH=~/anaconda3/bin:$PATH
を実行してください。
flexやtubeでは
#!/bin/sh
_CONDA_ROOT="/home/students/maru/anaconda3"
\. "$_CONDA_ROOT/etc/profile.d/conda.sh" || return $?
_conda_activate "$@"
というエラーが出ると思います。これはログイン先のシェルが、tcshなのに、shのスクリプトをsourceしてしまっているのが原因です。
tcshからはtcshスクリプトしかsource出来ません。
bash
と打って、ログインシェルをbashに変更してから実行してみてください。
実行はhogehoge.pyの部分を任意のファイルに置き換えてください。
(mp) $ python hogehoge.py
サンプルソースはGithub等で入手できるので[[このソースコード:https://github.com/stevengj/meep/blob/master/python/examples/wvg-src.py]]で試してみてください。~
無効かするには以下のコマンドを実行してください。~
(mp) $ source deactivate
**MEEPの並列化 [#k1f4cae9]
***インストール([https://meep.readthedocs.io/en/latest/Installation/#mpi 参考]) [#v7b0b7df]
並列処理をサポートしたpymeepをインストールするには上記の作業をすべて行ったうえで次のコマンドを実行してください。
$ conda create -n pmp -c chogan -c conda-forge pymeep-parallel
***実行 [#j3e26633]
$ source activate pmp
にて並列処理をサポートしたpymeepをアクティベートしたのち、
$ mpirun -np 4 python <script_name>.py
を実行してください。-np 4とすると4プロセッサで並列処理されます。数字は適当に変更してください。
***別の並列化 [#e415ae6a]
meep <file-name>.ctl &
とすると、バックグランドでプログラムが実行される。可変パラメータを外部から渡すことで、多数のプログラムを同時に動かせる。例えば、次のようなシェルスクリプトを著者は使用している。入射角を0から85度まで5度ずつ変えたプログラムを同時に走らせている。
for((i=0;i<90;i+=5))
do
meep incident-angle=${i} Main.ctl &
done
**HDF5の編集 [#ic8d7598]
-python向けにh5pyというライブラリが提供されている。以下でインストール。
apt-get install python-h5py
丸岡のbash on windowsではpip install h5pyでは、インストールはできても実行できなかった（ImportError: No module named h5py）。
**左円偏光 [#u95702e4]
次のコードで左回りの円偏光になります。
amplitudeによって、ソースの位相を変化することができます。
(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]
次のコードで右回りの円偏光になります。
amplitudeによって、ソースの位相を変化することができます。
(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]
x方向に直線偏光した平面波とy方向に直線偏光した平面波を
中心(0,0,8)、大きさ50×50×0の平面から、xy平面について45度方向((x,y)=(1,1)方向)に直線偏光した平面波を出力するコード。
(set! sources
(list
(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)))
)
(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)))
)
)
)
**斜め入射 [#e6b41572]
斜め入射のシミュレーションの方法は3つあります。
***全反射直前のTM波(ブロッホ周期境界条件) [#s8084013]
(define (pw-amp k x0) (lambda (x)
(exp (* 0+1i (vector3-dot k (vector3+ x x0))))))
(define-param theta-deg 48)
(define theta-rad (deg->rad theta-deg))
(define epsilon1 2.25)
(define epsilon2 1.25)
(define-param fcen 0.8) ; pulse center frequency
(define-param df 0.02) ; turn-on bandwidth
(define-param kdir (vector3 (sin theta-rad) 0 (cos theta-rad))) ; direction of k (length is irrelevant)
(define k (vector3-scale (* 2 pi fcen (sqrt epsilon1))
(unit-vector3 kdir))) ; k with correct length
(set! pml-layers (list (make pml (thickness 2)(direction Z))))
(set! k-point (vector3* (* fcen (sqrt epsilon1) ) (vector3 (sin theta-rad) 0 (cos theta-rad))))
(set! geometry
(list
(make block
(center 0 0 -2.5)
(size 10 infinity 5)
(material (make dielectric(epsilon epsilon1)))
)
(make block
(center 0 0 2.5)
(size 10 infinity 5)
(material (make dielectric(epsilon epsilon2)))
)
)
)
(set! sources
(list
(make source
(src (make continuous-src (frequency fcen) (fwidth df)))
(component Ez) (center 0 0 -3) (size 10 0 0)
(amp-func (pw-amp k (vector3 0 0 -3))))))
(define-param T 40) ; run time
(run-until T
(to-appended "ez" (at-every 0.1 output-efield-z))
(to-appended "ex" (at-every 0.1 output-efield-x))
)
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/main-ez.gif Ezのプロット]
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/main-ex.gif Exのプロット]
***光源を2つ並べる。 [#d178c899]
[https://github.com/stevengj/meep/blob/master/scheme/examples/pw-source.ctl コード]
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/pw-source-ez.gif Ezのプロット]
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/pw-source-single-ez.gif 光源が一つしかない場合のEzのプロット]
***光源を斜めにたくさん並べる。 [#t6b81813]
[http://kenji-ishikawa.cocolog-nifty.com/plasma/2010/01/post-ea91.html 動画とコード]
**電荷分布を見たい [#l0ec2ad8]
MEEPのbuilt-in関数には電荷分布を求める機能がないので、自作する必要がある。
ガウスの法則を差分化すればよい。
前述の全反射直前のTM波(ブロッホ周期境界条件)のシミュレーションについて例を示す。
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/charge.gif 電荷分布(+電場ベクトル)のプロット]
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/charge.py コード]
**物性値のライブラリ [#ue28f25c]
物性値のライブラリが公開されています。[https://github.com/stevengj/meep/blob/master/scheme/examples/material-dispersion.ctl materials-library.scm]をダウンロードしてください。
ライブラリを使用するには、
(include "/path/to/materials-library.scm")
と、プログラムの最初に追記してください。pathは自分がmaterials-library.scmをダウンロードした箇所に書き換えてください。
**hdf5ファイルの取り扱い [#de839d03]
pythonによるhdf5ファイルの取り扱いを説明する。
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に取り出せる。注意として、f['ex'][i]のようにして、何回もアクセスしないこと。このようなアクセスは非常に時間がかかる。
print(Ex) #print(h5file['ex'].value)も同様の動作をする。
とすると、Exの配列の中身が表示される。要素が多すぎるときは・・・と表示されて省略される。
Exが二次元配列でEx[i][j]と参照するとしたとき、iを第一インデックス、jを第二インデックスと呼ぶことにする。
第一インデックスの長さ(iの動く長さは)
print(len(Ex))
で分かる。第二インデックスの長さは
print(len(Ex[0]))
で分かる。
***印加電場と実際の電場の時間変化を折れ線グラフでプロット [#q9e259c0]
import h5py
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from PIL import Image
import moviepy.editor as edit
from mpl_toolkits.axes_grid1 import make_axes_locatable
plt.rcParams["font.size"] = 18
#----------input----------------
surfaceExFilename ="After-ex_+0.5.h5"
surfaceEyFilename = "After-ey_+0.5.h5"
surfacePxFilename = "Pre-ex_+0.5.h5"
surfacePyFilename = "Pre-ey_+0.5.h5"
savename = "current-line-withPreE-takingMiddleCurrent.png"
startTime = 0
endTime = 1000
waveLength = 6
timeResolution = 10
xid=16 #array[xid][yid][:]の電場をプロット
yid=46
#--------------------------------
def plotf(arr,label_):
time=np.array([])
amplitude=np.array([])
for i in range(startTime, endTime):
x=i * 1./(waveLength * timeResolution)
y=arr[xid][yid][i]
time=np.append(time,[x])
amplitude=np.append(amplitude,[y])
plt.plot(time,amplitude,label=label_)
plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left")
h5file_A_x = h5py.File(surfaceExFilename, "r")
h5file_A_y = h5py.File(surfaceEyFilename, "r")
h5file_P_x = h5py.File(surfacePxFilename, "r")
h5file_P_y = h5py.File(surfacePyFilename, "r")
Ax = h5file_A_x['ex'][:]
Ay = h5file_A_y['ey'][:]
Px = h5file_P_x['ex'][:]
Py = h5file_P_y['ey'][:]
***楕円率をヒストグラムで表したい [#g4d319c5]
Ellは2次元配列でxyに応じて楕円率が入っていると想定~
ヒストグラムの正規化はググると「normed=True」が出てくるがうまくいかなかった
[https://stackoverflow.com/questions/3866520/plotting-histograms-whose-bar-heights-sum-to-1-in-matplotlib/16399202 stackoverflow]が役に立った~
[https://jb102.blogspot.com/2017/10/22-histogram.html 日本語解説ブログ]~
ecのキーワードでヒストグラムの縦線の色を指定できる([https://stackoverflow.com/questions/42542252/cannot-get-histogram-to-show-separated-bins-with-vertical-lines 参考])
Ell_ = Ell.reshape(-1,)#plt.histは1次元配列を対象とするのでEllを1次元化
weights_ = np.ones(len(Ell_))/float(len(Ell_))#ヒストグラムを正規化するための処理
fig = plt.figure()
ax = fig.add_subplot(111)
#binsはヒストグラムの縦棒の数
ax.hist(Ell_, weights=weights_, bins=21, range=(-1,0,1.0),ec='black')
ax.set_xlim(-1.0,1.0)
ax.set_ylim(0.0,0.5)
fig.savefig("histgram.png")
array.reshape(-1)の挙動
>>> a = np.linspace(1,8,num=8)
>>>
>>> a
array([1., 2., 3., 4., 5., 6., 7., 8.])
>>> b = a.reshape((2,4))
>>> b
array([[1., 2., 3., 4.],
[5., 6., 7., 8.]])
>>> c = a.reshape((2,2,2))
>>> c
array([[[1., 2.],
[3., 4.]],
[[5., 6.],
[7., 8.]]])
>>> d = c.reshape(-1,)
>>> d
array([1., 2., 3., 4., 5., 6., 7., 8.])
代わりにweightsオプションを指定すればこの問題は解決できる. これはdataの1つの値の重みを指定するパラメータで, デフォルトでは実質的に
weights = np.ones(len(data))
を指定していることになる. いま総和を1にしたいのだから, ひとつのデータの重みは1/len(data)である.
ndarrayの場合は配列 / scalarが可能。
>>> x = [1,2,3]
>>> x
[1, 2, 3]
>>> x /2
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: unsupported operand type(s) for /: 'list' and 'int'
>>> y = np.array([1,2,3])
>>> y / 2
array([0.5, 1. , 1.5])
楕円率のヒストグラムの作成の際の
hist()のキーワードに
range=(-1,1)
を入れなないと空気を読んで，最大値と最小値の間をbinsの数だけ分割するようです
ax.set_xlim(-1., 0, 1., 0)
だと，-1から0までにしかプロットされないので
ax.set_xlim(-1.,1.)
に修正
**Q&A [#y2d911b6]
***一つのファイル内で複数回シミュレーションを呼び出したときと、一つひとつ呼び出したときの、結果が整合しない [#b9c73d94]
-1つのファイル内で複数回シミュレーションをするときは、実行を終了するたびに(reset-meep) を呼び出す。
***（入射波＋反射波）のフラックスが入射波のみの時より、小さくなる [#ua9ad9fe]
--フラックスは座標軸の正の方向を正に取っているため（ポインティングベクトルで定義されている）。正の方向に入射波が進行し、負の方向に反射波が進行する時、入射波だけの時よりフラックスは小さくなる。
***dielectricよりmediumでindexを指定した方が計算が早い（計算結果は同じなのに） [#xe0c2da7]
***resolutionをsetしても変化しない [#g280b7dd]
-resolutionはset! geometry-latticeの前に置かないとデフォルト値10になってしまうことがある。
***(meep-time)でMEEP系の単位における現在時刻を出力 [#d35030a2]
***変数vを使っていないのにunbounded variable vというエラーが出る。 [#ubd3a6e9]
-(define lambda 5)とするとunbounded variable vとしてエラーが出る。schemeではlambdaを用いて関数を定義するからだと思う（予約語になっている）
***meep: Could not determine normal direction for given grid_volume. [#ta6b98ee]
-fluxのサイズの少なくとも一つを0にしないとエラーが出る。
***実行時にGUILE_WARN_DEPRECATEDという警告が出る。 [#ma3f1cad]
-.loginに以下を記述すると、非推奨(deprecated)な機能を使用しているというエラーを握りつぶせる。
setenv GUILE_WARN_DEPRECATED no
***unbound variable: とエラーが出るが、:以降が空白になっている。 [#i0e856b1]
「unbound variable: 」と、エラーが出てるにもかかわらず、unbound variableの名前が表示されないときは、全角空白文字がどこかに含まれている。全角空白を削除すればよい。
***コマンドラインから変数を渡しているのに無視される [#dbd0a076]
meep hoehoge.ctl var=1
とすると、var=1は無視される。
meep var=1 hogehoge.ctl
とする。
***直方体を斜めに配置したい。 [#x6b589ab]
e1,e2,e3[vector3]により、方向を定められる。
例えば、xy平面内でx軸からΘだけ傾けて置きたいならば、
(make block
(size a b c)
(e1 (vector3 (cos theta) (sin theta) 0))
(e2 (vector3 (- (sin theta)) (cos theta) 0)
(e3 (vector3 0 0 1))
)
***x,y,z方向のpmlについて、厚さをそれぞれ別々にしたい [#z61ce405]
X, Y, Z方向をそれぞれ厚さa, b, cにしたいなら、
(set! pml-layers (list (make pml (thickness a)(direction X))))
(set! pml-layers (list (make pml (thickness b)(direction Y))))
(set! pml-layers (list (make pml (thickness c)(direction Z))))
とする。
***負の透過率、反射率、吸収率を示す [#j49576fc]
負の透過率、反射率、吸収率は電磁波が発散しきっていないのにシミュレーションを中断したことが原因のことが丸岡の場合多い。
(run-sources+ T
(stop-when-fields-decayed dt Ex (vector3 x y z) 1e-3))
とすると、T秒間シミュレーションをとりあえず回した後に、(x,y,z)におけるdt秒間のExの最大値が、全ての時間を通しての最大値と比較して1e-3倍以下になるまで回してくれる。
しかしこの関数は落とし穴があって、非常に長時間共振するような構造の場合は、まだ、電磁波の発振が続いているのにシミュレーションを中断する場合がある。
電磁波が発散しきっているかどうかの判断は、h5topngなどで電磁波を画像化してを直接確認するしかない。
***伝導度の理論値を複数のローレンツの重ね合わせで表したい [#p43befa3]
下のリンクの方法で取り組んでいました。理論値を複数のローレンツの重ね合わせで合わせるのはよくなさそうです。
https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg05147.html
***インターフェイスはpython, schemeどちらを使うべきか [#m1cc6868]
マニュアルのトップページにschemeがobsoleteでpythonに変えろと書いてあったことと、メーリングリストで（バグはまだ多いですが）pythonに機能が追加されていることもあって、先週からpythonに移行のために動作検証しているところです。
***PMLの厚さはどのくらいにすべきか？ [#l8406a8a]
PMLに垂直に入射するときは半波長。
以下、MEEP作成者のPMLに関するドキュメン(http://math.mit.edu/~stevenj/18.369/pml.pdf)からの引用
With PML, however, the constant factor is very good to start with, so experience shows that a simple quadratic or cubic turn-on of the PML absorption usually produces negligible refections for a PML layer of only half a wavelength or thinner [1,21].
(Increasing the resolution also increases the resolution alos increase the effectiveness of the PML, because it approaches the exact wave equation.)
斜め入射の時はPMLの吸収率が落ちるので、厚みを増す必要がある。
反射率は入射角をXとすると、ある定数δを用いて
e^{-2*k*d*cos(X)*δ}
と書ける。dはPMLの厚さ。kは波数。(詳細はhttps://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=477075を参照してください)。
***tube,flexにおいてmaterials.pyが使えない。[#m4142e79]
高専ではmaterials.pyが使用できるが、tube,flexでは使用できていない(2019/01/20:使用できるようになった)。次のようなエラーが出る。
tube60:~$ python
Python 3.6.7 | packaged by conda-forge | (default, Nov 21 2018, 03:09:43)
[GCC 7.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import meep as mp
>>> from meep.materials import Au
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'meep.materials'
高専の環境では~/anaconda3/envs/mp/lib/python3.6/site-packages/meep/にmaterials.pyが存在するが、flex,tubeには存在せずこのことが原因の可能性がある。
flex,tubeの~/anaconda3/pkgs/pymeep-1.5-py36_nomkl_2が高専では~/anaconda3/pkgs/pymeep-1.6-py36_nomkl_2になっており、バージョンを1.5から1.6にすることで使用できるようになる可能性がある。
ネット上からmaterias.pyをダウンロードしてホームディレクトリに置いたのちに実行すると、実行はできるものの次のようなワーニングが出てくる。
>>> from materials import Au
/home/guest/maeda/.conda/envs/mp/lib/python3.6/site-packages/meep/geom.py:175: RuntimeWarning: Epsilon < 1 may require adjusting the Courant parameter. See the 'Numerical Stability' entry under the 'Materials'
section of the documentation
warnings.warn(eps_warning, RuntimeWarning)
高専ではワーニングなしに実行できる。
オプションをつけてpymeepをupdateすると実行できるようになった。
(mp) maru@flex:~$ conda update -c chogan -c conda-forge pymeep
Solving environment: done
## Package Plan ##
environment location: /home/students/maru/anaconda3/envs/mp
added / updated specs:
- pymeep
The following packages will be downloaded:
package | build
---------------------------|-----------------
libctl-4.1.4 | 1 97 KB chogan
libgdsii-0.2.dev | 0 800 KB chogan
python-3.6.6 | hd21baee_1003 29.0 MB conda-forge
openssl-1.0.2p | h14c3975_1002 3.1 MB conda-forge
mpb-1.7.0 | nomkl_1 74 KB chogan
pymeep-1.7.0 | py36_nomkl_1 3.3 MB chogan
openblas-0.3.5 | h9ac9557_1000 15.8 MB conda-forge
nomkl-3.0 | 0 48 KB
ca-certificates-2018.11.29 | ha4d7672_0 143 KB conda-forge
certifi-2018.11.29 | py36_1000 145 KB conda-forge
------------------------------------------------------------
Total: 52.5 MB
The following NEW packages will be INSTALLED:
libgdsii: 0.2.dev-0 chogan
nomkl: 3.0-0
The following packages will be UPDATED:
ca-certificates: 2018.03.07-0 --> 2018.11.29-ha4d7672_0 conda-forge
certifi: 2018.11.29-py36_0 --> 2018.11.29-py36_1000 conda-forge
libctl: 4.1.2-0 chogan --> 4.1.4-1 chogan
mpb: 1.6.2-4 chogan --> 1.7.0-nomkl_1 chogan [nomkl]
openblas: 0.2.20-8 conda-forge --> 0.3.5-h9ac9557_1000 conda-forge
pymeep: 1.5-py36_nomkl_2 chogan [nomkl] --> 1.7.0-py36_nomkl_1 chogan [nomkl]
The following packages will be DOWNGRADED:
openssl: 1.1.1a-h7b6447c_0 --> 1.0.2p-h14c3975_1002 conda-forge
python: 3.6.8-h0371630_0 --> 3.6.6-hd21baee_1003 conda-forge
Proceed ([y]/n)? y
Downloading and Extracting Packages
libctl-4.1.4 | 97 KB | ########################################################################################### | 100%
libgdsii-0.2.dev | 800 KB | ########################################################################################### | 100%
python-3.6.6 | 29.0 MB | ########################################################################################### | 100%
openssl-1.0.2p | 3.1 MB | ########################################################################################### | 100%
mpb-1.7.0 | 74 KB | ########################################################################################### | 100%
pymeep-1.7.0 | 3.3 MB | ########################################################################################### | 100%
openblas-0.3.5 | 15.8 MB | ########################################################################################### | 100%
nomkl-3.0 | 48 KB | ########################################################################################### | 100%
ca-certificates-2018 | 143 KB | ########################################################################################### | 100%
certifi-2018.11.29 | 145 KB | ########################################################################################### | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
***ImportError: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.14' not found (required by /home/students/maru/anaconda3/envs/mp/lib/python3.6/site-packages/meep/_meep.so) [#zcd4847a]
flexのGLIBCのバージョンが古いので実行できていない。t70に移動して実行するとよい。
**参考になるサイト [#f39b395a]
-[https://core.ac.uk/download/pdf/38256115.pdf ターゲット特性による超高強度レーザー生成高速電子の制御に関する研究(大阪大学リポジトリ]
-[http://kenji-ishikawa.cocolog-nifty.com/plasma/ プラズマナノテクノロジー(ブログ)]
*コードの解析 [#c4799b48]
git clone https://github.com/stevengj/meep
./autogen.sh
make
とするとコードの完全体が得られる。
[https://github.com/stevengj/libctl/ libctl]というschemeの自作ライブラリが読み込まれているため読解が困難になっている。
**meep.scm.in [#o9d56706]
-meep.scm.inのdefine群が読み込まれている。
-ctlファイルで使用する関数の定義が書かれている。
Autotools で処理されるべきファイルに .in という拡張子が使われます。 つまり、その段階ではまだ Scheme のコードのテンプレートであり、処理されることで Scheme のコードとして妥当なものになります。ビルドプロセスによりますが、通常は configure が済んだ段階では処理済みのもの (この場合は meep.scm) が出来ていると思います。 もちろん make が済んだ時点なら確実に出来ているでしょう。
***クラスの定義 [#qa52258e]
(define-class A B (define-property C D E))とすると、Aという名前のクラスが定義される。Bはオーバーライドするクラスで、BのプロパティをAは引き継ぐ。Bをno-parentにすると、オーバーライドは起きない。
Aに対してCというプロパティが設定され、その初期値はDとなる。またCの型はEとなる。
以下にmeep.scm.inに含まれるクラスの定義の例を示す。
(define-class susceptibility no-parent
(define-property sigma-offdiag (vector3 0 0 0) 'vector3)
(define-property sigma-diag no-default 'vector3))
(define-class lorentzian-susceptibility susceptibility
(define-property frequency no-default 'number)
**meep.hpp [#ybf06d67]
- シミュレーション本体(C++)のすべてのクラスの前方宣言が含まれている。
**ユーザー定義の物質を使用する方法 [#yd375859]
[https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg04155.html]
-epsilon_{\omega=\infty}は、structureクラスのset_epsilon関数を使用する（多分）
例えば次のような感じ。(test/stress_tensor.cpp)
double two_waveguides(const vec &p) {
if ((fabs(p.x()) >= 0.5*d) && (fabs(p.x()) <= 0.5*d+sw) && (p.y() <= 0.5*sw) && (p.y() >= -0.5*sw))
return 11.9;
else
return 1.0;
}
int main(){
structure s(gv, two_waveguides, pml(dpml,X)+pml(dpml,Y), S);
s.set_epsilon(two_waveguides);
}
*他の電磁界シミュレーションソフト [#qf1fc505]
株式会社EEMが無料で公開している、OpenFDTDとOpenMOMは[[OpenFDTD DLサイト:http://www.e-em.co.jp/OpenFDTD/index.html]]、[[OpenMOM DLサイト:http://www.e-em.co.jp/OpenMOM/index.html]]からダウンロードすることができます。使い方については、DLサイトならびに、研究室が所蔵している[[RFワールド No.39 2017年 08月号 :https://www.amazon.co.jp/RF%E3%83%AF%E3%83%BC%E3%83%AB%E3%83%89-No-39-2017%E5%B9%B4-08-%E6%9C%88%E5%8F%B7/dp/B07233YBBQ]]に記載されており、書籍の方は現在（2018/6/12)丸岡が借りています。
![[PukiWiki]](image/c60.png "[PukiWiki]")
