• トップ
  • ブログ一覧
  • 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...

    ライトコードでは、エンジニアを積極採用中!

    ライトコードでは、エンジニアを積極採用しています!社長と一杯しながらお話しする機会もご用意しております。そのほかカジュアル面談等もございますので、くわしくは採用情報をご確認ください。

    採用情報へ

    広告メディア事業部
    広告メディア事業部
    Show more...

    おすすめ記事

    エンジニア大募集中!

    ライトコードでは、エンジニアを積極採用中です。

    特に、WEBエンジニアとモバイルエンジニアは是非ご応募お待ちしております!

    また、フリーランスエンジニア様も大募集中です。

    background