[meep] オープンソースのFDTDシミュレータを試してみよう

無料で本格的な光導波路のFDTDシミュレーションを始めよう

FDTD（Finite-Difference Time-Domain）法は光学デバイス設計において最も広く利用されているシミュレーション手法の一つです。
しかし商用ソフトウェアは強力な一方で、費用面のハードルが高く、気軽に試すのが難しいと感じるユーザーも多いはずです。

そこで今回は、オープンソースでありながら本格的な電磁界解析が可能な「Meep」をご紹介し、さらに実際に動作する導波路の透過スペクトル計算スクリプトを解説します。

「まずは無料でFDTDを試したい」「光導波路の基礎計算をやってみたい」
そんな方に向けた入門記事です。

1. Meepとは？

**Meep（MIT Electromagnetic Equation Propagation）**は、MITにより開発されたオープンソースのFDTD電磁界シミュレータです(https://meep.readthedocs.io/en/latest/Installation/)。
Python から操作できる PyMEEP も研究用途・教育用途で多く利用されています。

特徴
無料で利用可能（オープンソース）
商用FDTDのライセンス費用に手を出しにくいユーザーにとって大きなメリットです。

2D/3D FDTD をサポート

光導波路、レンズ、フォトニック結晶など幅広い構造解析に対応。

Python から柔軟に制御

Numpy・Matplotlib と併用して、
「入力 → 計算 → 可視化 → 解析」
のループを Python だけで完結できます。

導波路モード励振、PML、フラックス計算など高機能

商用ソフトと同様のツール（モードソース、PML、パワーモニタなど）が揃っており、波動素子の解析が十分実用的に行えます。

2. 今回実施するシミュレーション

シリコン導波路（幅 0.5 µm）を通過する光の透過スペクトルを計算します。

入力側：固有モード励振（EigenModeSource）

出力側：フラックスモニタで透過パワーを取得

波長 1.5–1.6 µm の範囲をスキャン

計算結果をグラフ表示

「有償FDTDでよく行う透過スペクトル解析」を PyMEEP だけで再現するイメージです。

3.スクリプト全体

こちらが今回作成したスクリプトになります。
とはいっても生成AIがコードを生成してくれる時代なので、やったことと言えば要件定義だけですが。。。

import meep as mp
import numpy as np
import matplotlib.pyplot as plt

def run_waveguide_transmission_spectrum(
    wl_min=1.5, wl_max=1.6, nfreq=200,
    filename="transmission_spectrum.csv"
):
    # ---------- 基本パラメータ ----------
    resolution = 50           # [pixels/um]

    # 波長レンジ → 周波数レンジ
    f_min = 1 / wl_max
    f_max = 1 / wl_min
    fcen = 0.5 * (f_min + f_max)
    df = f_max - f_min

    # ---------- 導波路パラメータ ----------
    n_core = 3.45
    n_clad = 1.0
    wg_width = 0.5

    # セルサイズ
    dpml = 1.0
    sx = 16.0
    sy = 8.0
    cell_size = mp.Vector3(sx, sy, 0)
    pml_layers = [mp.PML(dpml, direction=mp.X)]

    # ---------- 構造 ----------
    core = mp.Medium(index=n_core)
    clad = mp.Medium(index=n_clad)

    geometry = [
        mp.Block(
            size=mp.Vector3(mp.inf, wg_width, mp.inf),
            center=mp.Vector3(0, 0),
            material=core
        )
    ]

    # ---------- ソース ----------
    src_x = -0.5 * sx + dpml + 1.0
    sources = [
        mp.EigenModeSource(
            src=mp.GaussianSource(frequency=fcen, fwidth=df),
            center=mp.Vector3(src_x, 0),
            size=mp.Vector3(0, sy - 2 * dpml, 0),
            direction=mp.X,
            eig_band=1,
            eig_parity=mp.ODD_Z
        )
    ]

    # ---------- シミュレーション ----------
    sim = mp.Simulation(
        cell_size=cell_size,
        boundary_layers=pml_layers,
        geometry=geometry,
        sources=sources,
        default_material=clad,
        resolution=resolution,
        dimensions=2
    )

    # ---------- フラックスモニタ ----------
    flux_height = wg_width + 2.0
    flux_x_in = src_x + 0.5
    flux_x_out = 0.5 * sx - dpml - 0.5

    flux_in = sim.add_flux(
        fcen, df, nfreq,
        mp.FluxRegion(center=mp.Vector3(flux_x_in, 0),
                      size=mp.Vector3(0, flux_height, 0))
    )
    flux_out = sim.add_flux(
        fcen, df, nfreq,
        mp.FluxRegion(center=mp.Vector3(flux_x_out, 0),
                      size=mp.Vector3(0, flux_height, 0))
    )

    # ---------- 実行 ----------
    sim.run(until_after_sources=mp.stop_when_fields_decayed(
        50, mp.Ez, mp.Vector3(0, 0), 1e-7
    ))

    # ---------- スペクトル計算 ----------
    freqs = np.array(mp.get_flux_freqs(flux_in))
    wavelengths = 1 / freqs

    Pin = np.array(mp.get_fluxes(flux_in))
    Pout = np.array(mp.get_fluxes(flux_out))

    T = np.where(Pin != 0, Pout / Pin, 0.0)

    # 並べ替え（λ昇順）
    idx = np.argsort(wavelengths)
    wavelengths = wavelengths[idx]
    T = T[idx]
    Pin = Pin[idx]
    Pout = Pout[idx]

    # ---------- CSV 保存 ----------
    #data = np.column_stack([wavelengths, T, Pin, Pout])
    #header = "wavelength_um,transmission,Pin,Pout"
    #np.savetxt(filename, data, delimiter=",", header=header, comments="")
    #print(f"Saved CSV: {filename}")

    # ---------- グラフ表示 ----------
    plt.figure(figsize=(7,4))
    plt.plot(wavelengths, T)
    plt.xlabel("Wavelength (µm)")
    plt.ylabel("Transmission")
    plt.title("Transmission Spectrum")
    plt.grid(True)
    plt.tight_layout()
    plt.show()

    return wavelengths, T


if __name__ == "__main__":
    run_waveguide_transmission_spectrum()

今回はwindowsのwsl2を入れてシミュレーションを実施しました。

4. スクリプトの解説

4-1. 基本パラメータ設定

resolution = 50  # pixels/um
wl_min=1.5; wl_max=1.6
f_min = 1 / wl_max
f_max = 1 / wl_min

FDTD は周波数領域で計算するため、波長レンジを周波数に直変換します。

resolution=50 は 1 µm あたり 50 ピクセルの空間離散化。

商用FDTDと同じく、波長スキャンは周波数帯域を与えて一度の計算で済ませてます。

4-2. 導波路構造の定義

n_core = 3.45  # Si
n_clad = 1.0   # Air
wg_width = 0.5
geometry = [
    mp.Block(
        size=mp.Vector3(mp.inf, wg_width, mp.inf),
        material=core
    )
]

コア層はシリコン（n=3.45）

2D-FDTD なので高さ方向 z は mp.inf で無限扱い

導波路幅 0.5 µm の典型的なSOI導波路

4-3. モード励振（EigenModeSource）

mp.EigenModeSource(
    src=mp.GaussianSource(frequency=fcen, fwidth=df),
    eig_band=1,
    eig_parity=mp.ODD_Z
)

PyMEEPの強力な機能で、

TE/TMモードを指定して励振

商用ソフトの「Mode Source」「Eigenmode Source」と完全に同じ機能

4-4. フラックス計算

flux_in  = sim.add_flux(...)
flux_out = sim.add_flux(...)
Pin = np.array(mp.get_fluxes(flux_in))
Pout = np.array(mp.get_fluxes(flux_out))
T = Pout / Pin

Meep の add_flux を使って、周波数ごとのパワーを同時に取得

商用FDTDにおける「パワーモニタ」「フラックス計算」と同等

4-5. 計算終了条件

sim.run(until_after_sources=mp.stop_when_fields_decayed(
    50, mp.Ez, mp.Vector3(0, 0), 1e-7
))

「場が十分減衰したら自動停止」という Meep の便利機能
→ 不要な計算を省いて高速化できます。

4-6. スペクトルのプロット

plt.plot(wavelengths, T)

結果は次のような透過スペクトルになります。

5. コード全体のポイントまとめ

無料でここまで本格的なFDTD解析が可能

Pythonだけで完結する再現性の高い研究ワークフロー

有償FDTDとほとんど同じ機能で、そのまま商用ソフトへスキル移行可能

導波路、フォトニック結晶、共振器などにも応用可能

特に光導波路系の研究者・学生には大変おすすめです。

生成AIにコードを書かせて色々と試してみると良いでしょう

6. まとめ

本記事では、オープンソースの FDTD シミュレータ PyMEEP を紹介し、光導波路の透過スペクトルを計算する具体例を示しました。

有償FDTDは高機能だがハードルが高い

PyMEEPは無料で導入・学習が容易

モードソース、PML、フラックス計算など基本機能は十分

Pythonと組み合わせることで柔軟な解析が可能

まずは PyMEEP で FDTD の基礎を学び、
必要に応じて商用ソフトへスムーズにステップアップする、という使い方がおすすめです。

※本記事は筆者個人の見解であり、所属組織の公式見解を示すものではありません。

問い合わせ

光学シミュレーションソフトの導入や技術相談、
設計解析委託をお考えの方はサイバネットシステムにお問合せください。

光学ソリューションサイトについては以下の公式サイトを参照：
👉 光学ソリューションサイト(サイバネット)

光学分野のエンジニアリングサービスについては以下の公式サイトを参照：
👉 光学エンジニアリングサービス(サイバネット)