Python（パイソン）で電磁場解析とは？

「電磁場解析」とは、マクスウェルの方程式を解くことにより、対象物と電磁場の相互作用を解析することです。

過去、マクスウェルの方程式から導出される偏微分方程式を、解析的に解くことを指していました。

現代では、COMSOL（コムソル）等の市販ソフトでも、最新光学が組み込まれた「電磁場解析」があります。

そして、論文で発表されるコードなども、最近は「Python（パイソン）」で書かれております。

そんな事情を踏まえて、Python でできる、電磁場解析を紹介したいと思います！

Pythonを使った光電磁場解析

今回は、2019年8月に発売された書籍「Pythonを使った光電磁場解析」を、取り上げたいと思います！

コードの説明が、日本語なのが嬉しいところです！

【参考書籍】
Pythonを使った光電磁場解析

「Pythonを使った光電磁場解析」のFDTD法プログラム

早速、「Pythonを使った光電磁場解析」を参考に進めていきたいと思います！

この教科書の良い点は、「電磁場解析の数式」と「Pythonコード」の両方が、一冊に、系統立てて載っているので、コードの理解がしやすいことです。

「Pythonを使った光電磁場解析」と「MEEP」の比較

少し脱線しますが、同じFDTD法によるシミュレーションプログラムに、「MEEP（ミープ）」というものがあります。

参考の数式とコードがある以外に、何が違うか比較してみましょう。

FDTD法のプログラムを見ると、MEEPにはない「runfdtd.py」、「fdtd.py」、「preprocess.py」の３つのファイルが存在しました。

このプログラムは、「PML終端（計算領域の壁で電磁波が反射しないように完全に吸収させる方法）」で、平面波を「TF/SF法（計算領域境界で入射波をうまく補正する方法）」で導入しています。

Pythonで動かしてみる

今回、解析する内容はコチラ。

このあたりは、自分で使うものに後から変更すれば良いと思います。

そして、環境の紹介です。

早速動かしてみましょう。

「Anaconda（アナコンダ）」から、Python（ファイル名）で動きます。

ですが、肝心のFDTD関連の3つのファイルに、エラーが起きます。

バックスラッシュのあとのタグ削除

メモ帳で開くと、FDTD の3つのファイルだけ、行間が1行おきになっています。

Spyder（スパイダー）で開くと、エディタから以下のメッセージが表示されます。

このまま「自動的に修正」で進み run します。

すると、「unexpected intend」と言われる箇所が続出しました。

（左ウィンドがコード。右ウィンドがコンソール）

バックスラッシュの次行の数式前スペースが「unexpected」とのことです。

こちらの自動修正エラーですが、「newline」や「コード変換」などで試行錯誤した結果、以下のコードで解決しました。

行間を削除するコード

1# fdtd.pyファイルから文字列を読み込みリスト化後、複数種の改行コードを全削除して、\nコードでjoin
2with open('fdtd.py', 'r', encoding='utf-8-sig') as a_file:
3    txt = a_file.read().splitlines()
4    txta=[s for s in txt if s!=""]
5    join_txta="\n".join(txta)
6
7#Nfdtd.pyとして保存
8with open('Nfdtd.py', 'wb') as a_file:
9    a_file.write(join_txta.encode('utf-8-sig'))
10
11# 同様にpreprocess.pyとrunfdtd.pyを、それぞれNpreprocess.pyとrunfdtd.pyというファイルで保存
12with open('preprocess.py', 'r', encoding='utf-8-sig') as b_file:
13    txt = b_file.read().splitlines()
14    txtb=[s for s in txt if s!=""]
15    join_txtb="\n".join(txtb)
16with open('Npreprocess.py', 'wb') as b_file:
17    b_file.write(join_txtb.encode('utf-8-sig'))
18
19with open('runfdtd.py', 'r', encoding='utf-8-sig') as c_file:
20    txt = c_file.read().splitlines()
21    txtc=[s for s in txt if s!=""]
22    join_txtc="\n".join(txtc)
23with open('Nrunfdtd.py', 'wb') as c_file:
24    c_file.write(join_txtc.encode('utf-8-sig'))

修正後

FDTD計算

いよいよ、計算していきたいと思います！

今回、計算に使うのは、

Intel(R) core(™) m5-6Y54 CPU@1.10 GHz 1.51 GHz

実装メモリ8 GBの、持ち運びPCです。

計算手順

１、新しい計算用フォルダを作成
２、先ほどのコードで作成された「N」から始まる3つのファイルを移動
３、ファイル名から「N」を抜く
４、計算用フォルダ内で、runfdtd.pyを、Anacondaで回す

計算する構造由来のパラメータは、runfdtd.py で設定します。

初期状態では、ガラス基板上に置かれた「半径25 nmの金のナノ粒子」に、真空中で「波長561 nmのパルス波」をいれた時の応答が計算されます。

光源は、双極子になっています。

もし、球以外のものを置きたい場合は、「preprocess.py」中の「def _sphere」の中を変えれば、望みの構造にできます。

メモ帳で確認

計算した構造体が、望みのものになっているかは、「eps」で始まる、誘電率の空間分布を示すファイルを、メモ帳などで開くと簡単に確認できます。

※確認する際は、「右端で折り返す」設定を外してください。

「epsx_y」を開くと、基板上（誘電率１）に、球（誘電率12）が載っている様子が見られます。

この、epsで始まるファイルは、計算途中でも開きます。

計算終了を待つ前に確認すると、時間が無駄にならずに済みます。

計算途中でも結果がテキストファイルで作られていく

計算用フォルダの中に、「field」と「_pycache_」という2つのフォルダができます。

計算中は、「fieldフォルダ」の中に電場・磁場計算結果のテキストファイルが、どんどん作られていきます。

いかにも、FDTD法らしいです。

待つこと4時間。

テキストの電場・磁場ファイルが出来上がりました。

テキストを空間分布へ

変数を全くいじらないダウンロード状態で、出来上がったファイルは、

「誘電率分布の、x、y、zファイル3種類」
「電場・磁場ファイル（各場で時間変化31ファイル）」
「ResponseとSpectraのファイル」

で、計133個。

中を開くと、どのファイルも数字が羅列しています。

次に、この数値の羅列を、空間分布として可視化します。

初学者であれば、Excel で開き、「条件付き書式」で表示するだけでカラーの電場分布が得られます。

Pythonで表示する場合のコード

さて、Pythonで結果を表示する場合です。

こちらのテキストデータは、初期設定では「縦97、横96 の行列」で、「空白区切り」なので、下記のコードに。

1import numpy as np
2import matplotlib.pyplot as plt
3
4x = np. arange(0,96,1) #ｘ軸
5y = np. arange(0,97,1) #y軸
6x, y = np.meshgrid(x,y) #xとyのmesh作成
7
8z = np.loadtxt(“Ex_y040_004.txt') #読みたいデータをnumpyで読み込み
9
10plt.pcolormesh(x, y, z, cmap='hsv') #等高線図
11pp=plt.colorbar(orientation="vertical") #カラーバー表示
12pp.set_label("Intensity", fontname="Arial",fontsize=18) #カラーバーラベル
13
14plt.xlabel('x', fontsize=12)
15plt.ylabel('y', fontsize=12)
16
17plt.show()

結果

さいごに

参考となる書籍があるので、改行コードの問題さえ解いてしまえば、「MEEP for windows」より中身が分かりやすく、作り変えやすいと思います。

特に、書籍「プラズモニクス―基礎と応用」と合わせると、実際にどんな理論が、どんな数式で計算されているのかが分かります。

特に、日本人の電磁場学習者に、とても良いと思いました。

良ければ、そちらも試して、理解を深めていただければと思います！

また、今回ノートPCを使ったのですが、計算しながらでも他の作業は軽快でした！

こちらの記事もオススメ！