• トップ
  • ブログ一覧
  • Pythonで電磁場解析やってみた!
  • Pythonで電磁場解析やってみた!

    広告メディア事業部広告メディア事業部
    2020.03.09

    IT技術

    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」の3つのファイルが存在しました。

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

    Pythonで動かしてみる

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

    物体形状球と平面基盤
    媒質真空、シリカ、金および銀

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

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

    OSWindows10
    python3.7

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

    「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です。

    計算手順

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

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

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

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

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

    メモ帳で確認

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

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

    「epsx_y」を開くと、基板上(誘電率1)に、球(誘電率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) #x軸
    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を使ったのですが、計算しながらでも他の作業は軽快でした!

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

    featureImg2020.07.17ライトコード的「やってみた!」シリーズ「やってみた!」を集めました!(株)ライトコードが今まで作ってきた「やってみた!」記事を集めてみました!※作成日が新し...
    featureImg2020.07.30Python 特集実装編※最新記事順Responder + Firestore でモダンかつサーバーレスなブログシステムを作ってみた!P...

    広告メディア事業部

    広告メディア事業部

    おすすめ記事