RCWA(FMM)のフリーライブラリーS4の紹介 - 光学設計とその周辺、そしてたまに全く関係ないやつ

RCWA(Rigorous Coupled-Wave Analysis)またはFMM(Fourier Modal Method)とはマクスウェル方程式を厳密に解く電磁場解析の手法の1つで, 回折格子のような周期性のある構造のみ解析できるが,FDTDに比べ高速に解析できるのが特徴です.

有償のソフトウェアももちろんありますが,オープンソースの無料のライブラリもちらほら存在する. 特にS4(Stanford Stratified Structure Solver)というライブラリが有名なようで自分もちょっと触ろうと思ったが, 日本語の情報も無く, 思ってたより環境を用意するのが大変だったので忘備録も兼ねて記事にしてみます.

S4の概要
実行環境の用意
実際の解析

S4の概要

S4とは名前の通りもともとはStanford大学で開発されたオープンソースのRCWAのライブラリです. 本元のHPは以下のリンク先にあり,ここに各種情報があります. ライセンスはGPLのため, 商用利用でも社内で使う分にはめんどくさいことはないはずです.
web.stanford.edu
また開発者によるS4を紹介した論文もオープンアクセスとして以下のリンク先にありますが理論的なことはそちらが詳しい.
www.sciencedirect.com

実行環境の用意

環境は以下の通り.仮想環境は商用利用に適しているvirtualboxを選択しましたがvmwareでも動くことは一応確認しました.
os:ubuntu-20.04.1
仮想環境:VirtualBox on Windows(v6.1.14)

当初直接Windows上でトライしたがうまくいかず断念. またWSL1でもうまくいかなかった記憶がありますが, より互換性の高い最新のWSL2なら動くかもしれない.

さてVirtualBoxでubuntuを初めて起動した状態が以下のとおりです.

ここからターミナル(windowsでいうコンソール)を開いて以下のコマンドを順に実行する.
sudo rm /var/lib/apt/lists/lock
sudo rm /var/lib/dpkg/lock

sudo apt-get update
sudo apt-get upgrade

S4はもともとはLuaというプログラミング言語用に開発されたライブラリだがPython用のインターフェイスも用意されている. 少なくともこの記事を書いている2021年はPythonのほうが圧倒的に利用者が多いため今回もPythonで実行することを目指す. ということで以下のコマンドを実行しpython環境を用意していきます.

sudo apt install python3-venv
sudo apt install python3-pip
sudo apt-get install libblas-dev liblapack-dev
alias python=python3

以上で環境の用意は終わりです. 恐らく使用する環境によってもっとやることはるのかもしれないですが,ここまでは通常のLinux, Python環境の設定なので困ったことがあったらネットに答えは転がっていると思います.
最後にS4自体のインストールを行います.
sudo apt install git
git clone https://github.com/phoebe-p/S4.git

ディレクトリに移動し、最後に以下のコマンドを実行
make S4_pyext
setup.py

pythonを起動させ、以下のコードを実行してもエラーが起きなければ問題なくインストールができています.

import S4
S = S4.NewSimulation(Lattice=((1,0),(0,1)), NumBasis=100)

実際の解析

S4をダウンロードするとそのサブフォルダであるS4/examples/に0次元、1次元、2次元それぞれにRCWAの解析を利用した論文と対応するサンプルコードがある.先述した通りいずれもluaであるが、コード自体は単純であるためpythonに移植するのは簡単なはず.このサンプルを見ていけばS4の使い方は十分理解できる.

せっかくなので、このうちの１つJ. Pietarinen, et al., "Wideband four-level transmission gratings with flattened spectral efficiency," Opt. Express 14, 2583-2588 (2006)をS4とpythonで解析してみよう。こちらもオープンアクセスなので無料で閲覧できるのがありがたい。
Optica Publishing Group

解析対象となるのはこのうちのFig2(a)の青線である。

さてこの構造を解析するためのコードは以下のとおりです.ざっくり言えば最初のS4.Newでフーリエ係数の数や周期構造の向きなどを指定し, SetMaterialで使用する光学特性を, AddLayerとSetRegionRectangleで各層の光学特性と形状を指定していく.

import S4						

d=11
x1=0.25
x2= 0.5
x3=0.75
z1=0.722
z2=1.444
z3=2.166
angle=0

S=S4.New(Lattice=((d,0),(0,0)),NumBasis=21)	

#define substance and structure.SiO2の屈折率はとりあえず適当
S.SetMaterial('Vacuum',1)
S.SetMaterial('SiO2',1.445*1.445)	

S.AddLayer('Substrate',0,'SiO2')

S.AddLayer('Level1',z1,'Vacuum')
S.SetRegionRectangle('Level1','SiO2',(0.5*(1+x1)*d,0),0, (0.5*(1-x1)*d,0))

S.AddLayer('Level2',z2-z1,'Vacuum')
S.SetRegionRectangle('Level2','SiO2',(0.5*(1+x2)*d,0),0, (0.5*(1-x2)*d,0))

S.AddLayer('Level3',z3-z2,'Vacuum')
S.SetRegionRectangle('Level3','SiO2',(0.5*(1+x3)*d,0),0, (0.5*(1-x3)*d,0))

S.AddLayer('Air',0,'Vacuum')

for i in range(100,205,5):
    if i==100:
        wl=1.001
    else:
        wl=i/100

    S.SetFrequency(1/wl)
    S.SetExcitationPlanewave((angle,0),1,0,1)
    (inc,back)=S.GetPowerFlux('Substrate',0)
    P=S.GetPowerFluxByOrder('Air',0)
    T=P[3][0]/inc
    print(wl)
    print(T.real)

元の論文の画像とS4で得た結果は以下の通り.結果が一致していることがわかってもらえるだろう。

simulation結果左図が元論文(出典：J. Pietarinen, et al., "Wideband four-level transmission gratings with flattened spectral efficiency"(2006) Fig2(a)より)、右図がS4による解析結果

何点かS4特有の注意点をまとめるとまず光速度と真空の誘電率、透磁率は1に規格化されているため波長を指定するときはSetFrequency(1/波長値)とします. また長さのスケールも規格化されていることに注意.
この点詳しくはhttps://web.stanford.edu/group/fan/S4/units.htmlを参照.