Pythonで電磁場解析やってみた！
2021.12.20
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で動かしてみる
今回、解析する内容はコチラ。
| 物体形状 | 球と平面基盤 |
| 媒質 | 真空、シリカ、金および銀 |
このあたりは、自分で使うものに後から変更すれば良いと思います。
そして、環境の紹介です。
| OS | Windows10 |
| python | 3.7 |
早速動かしてみましょう。
「Anaconda（アナコンダ）」から、Python（ファイル名）で動きます。
ですが、肝心のFDTD関連の3つのファイルに、エラーが起きます。
バックスラッシュのあとのタグ削除
メモ帳で開くと、FDTD の3つのファイルだけ、行間が1行おきになっています。
Spyder（スパイダー）で開くと、エディタから以下のメッセージが表示されます。
このまま「自動的に修正」で進み run します。
すると、「unexpected intend」と言われる箇所が続出しました。
（左ウィンドがコード。右ウィンドがコンソール）
バックスラッシュの次行の数式前スペースが「unexpected」とのことです。
こちらの自動修正エラーですが、「newline」や「コード変換」などで試行錯誤した結果、以下のコードで解決しました。
行間を削除するコード
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | # fdtd.pyファイルから文字列を読み込みリスト化後、複数種の改行コードを全削除して、\nコードでjoin with open('fdtd.py', 'r', encoding='utf-8-sig') as a_file: txt = a_file.read().splitlines() txta=[s for s in txt if s!=""] join_txta="\n".join(txta) #Nfdtd.pyとして保存 with open('Nfdtd.py', 'wb') as a_file: a_file.write(join_txta.encode('utf-8-sig')) # 同様にpreprocess.pyとrunfdtd.pyを、それぞれNpreprocess.pyとrunfdtd.pyというファイルで保存 with open('preprocess.py', 'r', encoding='utf-8-sig') as b_file: txt = b_file.read().splitlines() txtb=[s for s in txt if s!=""] join_txtb="\n".join(txtb) with open('Npreprocess.py', 'wb') as b_file: b_file.write(join_txtb.encode('utf-8-sig')) with open('runfdtd.py', 'r', encoding='utf-8-sig') as c_file: txt = c_file.read().splitlines() txtc=[s for s in txt if s!=""] join_txtc="\n".join(txtc) with open('Nrunfdtd.py', 'wb') as c_file: 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 の行列」で、「空白区切り」なので、下記のコードに。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | import numpy as np import matplotlib.pyplot as plt x = np. arange(0,96,1) #ｘ軸 y = np. arange(0,97,1) #y軸 x, y = np.meshgrid(x,y) #xとyのmesh作成 z = np.loadtxt(“Ex_y040_004.txt') #読みたいデータをnumpyで読み込み plt.pcolormesh(x, y, z, cmap='hsv') #等高線図 pp=plt.colorbar(orientation="vertical") #カラーバー表示 pp.set_label("Intensity", fontname="Arial",fontsize=18) #カラーバーラベル plt.xlabel('x', fontsize=12) plt.ylabel('y', fontsize=12) plt.show() |
結果
さいごに
参考となる書籍があるので、改行コードの問題さえ解いてしまえば、「MEEP for windows」より中身が分かりやすく、作り変えやすいと思います。
特に、書籍「プラズモニクス―基礎と応用」と合わせると、実際にどんな理論が、どんな数式で計算されているのかが分かります。
特に、日本人の電磁場学習者に、とても良いと思いました。
良ければ、そちらも試して、理解を深めていただければと思います！
また、今回ノートPCを使ったのですが、計算しながらでも他の作業は軽快でした！
こちらの記事もオススメ！
書いた人はこんな人
- 「好きを仕事にするエンジニア集団」の（株）ライトコードです！
ライトコードは、福岡、東京、大阪の3拠点で事業展開するIT企業です。
現在は、国内を代表する大手IT企業を取引先にもち、ITシステムの受託事業が中心。
いずれも直取引で、月間PV数1億を超えるWebサービスのシステム開発・運営、インフラの構築・運用に携わっています。
システム開発依頼・お見積もり大歓迎！
また、現在「WEBエンジニア」「モバイルエンジニア」「営業」「WEBデザイナー」「WEBディレクター」を積極採用中です！
インターンや新卒採用も行っております。
以下よりご応募をお待ちしております！
https://rightcode.co.jp/recruit
ライトコードの日常12月 1, 2023ライトコードクエスト〜東京オフィス歴史編〜- ITエンタメ10月 13, 2023Netflixの成功はレコメンドエンジン？
- ライトコードの日常8月 30, 2023退職者の最終出社日に密着してみた！
- ITエンタメ8月 3, 2023世界初の量産型ポータブルコンピュータを開発したのに倒産!?アダム・オズボーン
関連記事
採用情報
\ あの有名サービスに参画!? /
バックエンドエンジニア
\ クリエイティブの最前線 /
フロントエンドエンジニア
\ 世界を変える…！ /
Androidエンジニア
\ みんなが使うアプリを創る /
iOSエンジニア
