[Python×Ansys Lumerical×Ansys HFSS]進行波型LiNbO3光変調器の光周波数応答を求めてみよう

問い合わせフォームのご連絡

この記事の内容についてのご質問や、
Ansys（Lumerical / HFSS など）の使い方相談はもちろん、
フリーソフト・OSSを使ったシミュレーションや解析の技術的なご相談も受け付けています。
「この設定で合ってる？」「この解析、どう進めるのがよい？」
「商用ソフトとフリーソフト、どう使い分けるべき？」
といったレベル感でも問題ありません。
気軽にこちらの問い合わせフォームからご連絡ください👇
https://form.cybernet.co.jp/optical/contact/qiita_form/

はじめに

以前の記事では、Ansys Lumericalを用いて、位相シフタの群屈折率、Ansys HFSSでは進行波型電極の特性インピーダンス、屈折率を求めるところまでご紹介しました。
https://qiita.com/Takeshi_Baba/items/be840825c2c344e12d0e
https://qiita.com/Takeshi_Baba/items/07599ffb43020824bd31

次に今回は、Ansys Lumerical と Ansys HFSS を用いて事前に実行したシミュレーション結果を読み込み、電気光学（EO）変調器の周波数応答を解析的に評価するためのポストプロセス専用スクリプトをご紹介します。

光F特を求める解析式は、以前の記事の内容をもとにしています。
https://qiita.com/Takeshi_Baba/items/35e4a6f6bd920a3edad1

スクリプト

# -*- coding: utf-8 -*-
"""
EO frequency response (post-processing only)
- Assumes Lumerical and HFSS were run manually beforehand
- Reads their exported results and computes analytic EO response

Author: babatake
"""

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


# ============================================================
#  PATHS (edit as needed)
# ============================================================
PROJECT_DIR = r"D:\LN modulator"

# Lumerical output (from LN_main7.py)
LUM_CSV = os.path.join(PROJECT_DIR, "lumerical_loss_vpil.csv")

# HFSS output (from LN_CPW2.py)
HFSS_NPZ = os.path.join(PROJECT_DIR, "hfss_cpw_results.npz")
HFSS_CSV = os.path.join(PROJECT_DIR, "hfss_cpw_results.csv")  # optional fallback

# Output for EO response
OUT_EO_CSV = os.path.join(PROJECT_DIR, "eo_frequency_response.csv")


# ============================================================
#  USER SETTINGS
# ============================================================
# Use this if you don't have ng saved as a file.
NG_DEFAULT = 2.2　　#これは群屈折率の入力部分です。

# Modulator length used in EO response (override > HFSS(d_m) > default)
# - If L_M_OVERRIDE is not None, it will override HFSS-exported length (d_m) and L_M_DEFAULT.
# - If L_M_OVERRIDE is None, it will use HFSS NPZ d_m if available, otherwise L_M_DEFAULT.
L_M_DEFAULT  = 0     # [m] fallback
L_M_OVERRIDE = 1e-3     # [m] e.g. 0.1e-3 for 100 um, or None to auto

# Source/Load impedance
ZG = 50.0
ZL = 50.0
C0 = 3.0e8


# ============================================================
#  Helpers
# ============================================================
def _to_1d(a, dtype=float):
    a = np.asarray(a)
    a = np.ravel(a)
    return a.astype(dtype, copy=False)

def _safe_complex_from_npz(npz, key, fallback_key_real=None, fallback_key_imag=None):
    """
    Prefer complex array stored in NPZ. If not found, try reconstructing from real/imag keys.
    """
    if key in npz.files:
        arr = np.asarray(npz[key])
        return arr.astype(np.complex128) if not np.iscomplexobj(arr) else arr.astype(np.complex128)
    if (fallback_key_real in npz.files) and (fallback_key_imag in npz.files):
        re = npz[fallback_key_real]
        im = npz[fallback_key_imag]
        return np.asarray(re, dtype=float) + 1j * np.asarray(im, dtype=float)
    raise KeyError(f"Could not find complex array '{key}' in NPZ.")

def load_lumerical_csv(path=LUM_CSV):
    """
    Reads Lumerical CSV produced by LN_main7.py.
    Expected columns: Volt[V], L_pi[m], (optional) loss[dB/cm], VpiL_step[V*cm], VpiL_prod[V*cm]
    """
    if not os.path.isfile(path):
        raise FileNotFoundError(f"Lumerical CSV not found: {path}")

    df = pd.read_csv(path)
    if "Volt[V]" not in df.columns:
        raise ValueError("Lumerical CSV missing column: Volt[V]")

    volt = _to_1d(df["Volt[V]"].values, float)

    # Optional (not strictly needed for EO frequency response)
    lpi  = df["L_pi[m]"].values if "L_pi[m]" in df.columns else None
    loss = df["loss[dB/cm]"].values if "loss[dB/cm]" in df.columns else None

    return {
        "df": df,
        "Volt": volt,
        "L_pi_m": None if lpi is None else _to_1d(lpi, float),
        "loss_dB_per_cm": None if loss is None else _to_1d(loss, float),
    }

def load_hfss_results(npz_path=HFSS_NPZ, csv_path=HFSS_CSV):
    """
    Prefer NPZ because it can preserve complex arrays.
    If NPZ is absent, fall back to CSV (will reconstruct complex from real/imag columns).
    Returns:
      freq_Hz, Z0(complex), n_rf(complex), L_m(from NPZ d_m if any), source
    """
    if os.path.isfile(npz_path):
        npz = np.load(npz_path, allow_pickle=True)

        if "freq_Hz" not in npz.files:
            raise ValueError(f"HFSS NPZ missing key: freq_Hz ({npz_path})")
        if "Z0" not in npz.files:
            raise ValueError(f"HFSS NPZ missing key: Z0 ({npz_path})")
        if "n_rf" not in npz.files:
            raise ValueError(f"HFSS NPZ missing key: n_rf ({npz_path})")

        freq_Hz = np.asarray(npz["freq_Hz"], dtype=float)

        Z0 = np.asarray(npz["Z0"])
        n_rf = np.asarray(npz["n_rf"])

        # Ensure complex dtype
        Z0 = Z0.astype(np.complex128) if not np.iscomplexobj(Z0) else Z0.astype(np.complex128)
        n_rf = n_rf.astype(np.complex128) if not np.iscomplexobj(n_rf) else n_rf.astype(np.complex128)

        # Electrode length used in HFSS export (if present)
        L_m = None
        if "d_m" in npz.files:
            try:
                L_m = float(npz["d_m"])
            except Exception:
                L_m = None

        return {"freq_Hz": freq_Hz, "Z0": Z0, "n_rf": n_rf, "L_m": L_m, "source": "npz"}

    # --- CSV fallback ---
    if not os.path.isfile(csv_path):
        raise FileNotFoundError(
            "HFSS NPZ not found and CSV not found either:\n"
            f"  NPZ: {npz_path}\n"
            f"  CSV: {csv_path}"
        )

    df = pd.read_csv(csv_path)
    req = ["freq_Hz", "Z0_real", "Z0_imag", "n_rf_real", "n_rf_imag"]
    for c in req:
        if c not in df.columns:
            raise ValueError(f"HFSS CSV missing column: {c}")

    freq_Hz = _to_1d(df["freq_Hz"].values, float)
    Z0 = _to_1d(df["Z0_real"].values, float) + 1j * _to_1d(df["Z0_imag"].values, float)
    n_rf = _to_1d(df["n_rf_real"].values, float) + 1j * _to_1d(df["n_rf_imag"].values, float)

    return {"freq_Hz": freq_Hz, "Z0": Z0.astype(np.complex128), "n_rf": n_rf.astype(np.complex128), "L_m": None, "source": "csv"}

def compute_eo_response(freqs_Hz, Z0, n_rf, ng_used, L_m):
    """
    Analytic EO response used in your script (m_total = m1*m2).

    Definitions used (matching your implementation):
      omega     = 2*pi*f
      lambda_rf = c0/f
      beta_m    = 2*pi*Re(n_rf)/lambda_rf
      beta_o    = 2*pi*ng/lambda_rf
      alpha_m   = Im(n_rf)     (kept as-is to match your current formula)
      gamma_rf  = (omega/c0) * n_rf

      Z_in      = Z0 * (ZL + Z0*tanh(gamma_rf*L)) / (Z0 + ZL*tanh(gamma_rf*L))

      u_plus    = j*(beta_m - beta_o)*L + alpha_m*L
      u_minus   = j*(-beta_m - beta_o)*L - alpha_m*L

      F(u)      = (1 - exp(u))/u   (with small-u stabilization)

      m1        = | Z_in / (Z_in + ZG) |
      m2        = | ((ZL+ZG)*F(u_plus) + (ZL - Z0)*F(u_minus))
                   / ((ZL+ZG)*exp(gamma_rf*L) + (ZL - Z0)*exp(-gamma_rf*L)) |
      m_total   = m1*m2
    """
    freqs_Hz = np.asarray(freqs_Hz, dtype=float)
    Z0 = np.asarray(Z0, dtype=np.complex128)
    n_rf = np.asarray(n_rf, dtype=np.complex128)

    if freqs_Hz.size == 0:
        raise ValueError("freqs_Hz is empty.")
    if np.any(~np.isfinite(freqs_Hz)) or np.any(freqs_Hz <= 0):
        raise ValueError("freqs_Hz contains non-finite or non-positive values.")
    if not np.isfinite(L_m) or L_m <= 0:
        raise ValueError(f"L_m must be positive finite. Got {L_m}")

    omega = 2 * np.pi * freqs_Hz
    lambda_rf = C0 / freqs_Hz

    beta_m = 2 * np.pi * np.real(n_rf) / lambda_rf
    beta_o = 2 * np.pi * float(ng_used) / lambda_rf

    # NOTE: kept to match your existing implementation
    alpha_m = np.imag(n_rf)

    gamma_rf = (omega / C0) * n_rf  # complex propagation constant [1/m]

    # Input impedance of a transmission line with load ZL
    tanh_term = np.tanh(gamma_rf * L_m)
    Z_in = Z0 * (ZL + Z0 * tanh_term) / (Z0 + ZL * tanh_term)

    u_plus  = 1j * (beta_m - beta_o) * L_m + alpha_m * L_m
    u_minus = 1j * (-beta_m - beta_o) * L_m - alpha_m * L_m

    def _safe_F(u):
        eps = 1e-15
        u2 = np.where(np.abs(u) < eps, eps + 0j, u)
        return (1 - np.exp(u2)) / u2

    F_plus  = _safe_F(u_plus)
    F_minus = _safe_F(u_minus)

    m_1 = np.abs(Z_in / (Z_in + ZG))

    m_2_1 = (ZL + ZG) * F_plus + (ZL - Z0) * F_minus
    m_2_2 = (ZL + ZG) * np.exp(gamma_rf * L_m) + (ZL - Z0) * np.exp(-gamma_rf * L_m)
    m_2 = np.abs(m_2_1 / m_2_2)

    m_total = m_1 * m_2

    m_total_dB = 20 * np.log10(np.maximum(m_total, 1e-300))
    m_total_norm = m_total / np.nanmax(m_total)
    m_total_norm_dB = 20 * np.log10(np.maximum(m_total_norm, 1e-300))

    return {
        "m_total": m_total,
        "m_total_dB": m_total_dB,
        "m_total_norm": m_total_norm,
        "m_total_norm_dB": m_total_norm_dB,
        "m_1": m_1,
        "m_2": m_2,
    }

def estimate_3db_bandwidth(freqs_Hz, m_total_norm):
    """
    Finds the first frequency where 20log10(norm) <= -3 dB.
    If not reached, returns None.
    """
    freqs_Hz = np.asarray(freqs_Hz, dtype=float)
    y_dB = 20 * np.log10(np.maximum(np.asarray(m_total_norm, dtype=float), 1e-300))

    # Ensure frequency ascending
    idx_sort = np.argsort(freqs_Hz)
    f = freqs_Hz[idx_sort]
    y = y_dB[idx_sort]

    target = -3.0
    below = np.where(y <= target)[0]
    if len(below) == 0:
        return None

    k = below[0]
    if k == 0:
        return f[0]

    # Linear interpolation in dB
    f1, f2 = f[k-1], f[k]
    y1, y2 = y[k-1], y[k]
    if y2 == y1:
        return f2
    frac = (target - y1) / (y2 - y1)
    return f1 + frac * (f2 - f1)


# ============================================================
#  MAIN
# ============================================================
if __name__ == "__main__":

    # 1) Load Lumerical result (manual run)
    lum = load_lumerical_csv(LUM_CSV)

    # Decide ng used in EO response.
    # Manual-run assumption => use NG_DEFAULT (replace with your own selection logic if needed).
    ng_used = float(NG_DEFAULT)

    # 2) Load HFSS result (manual run)
    hf = load_hfss_results(HFSS_NPZ, HFSS_CSV)

    freqs_Hz = hf["freq_Hz"]
    Z0 = hf["Z0"]
    n_rf = hf["n_rf"]

    # 3) Decide modulator length L_m (override > HFSS(d_m) > default)
    if L_M_OVERRIDE is not None:
        L_m = float(L_M_OVERRIDE)
        L_src = "override"
    elif hf["L_m"] is not None:
        L_m = float(hf["L_m"])
        L_src = "hfss_npz(d_m)"
    else:
        L_m = float(L_M_DEFAULT)
        L_src = "default"

    if not np.isfinite(L_m) or L_m <= 0:
        raise ValueError(f"Invalid L_m: {L_m} (source={L_src})")

    print(f"[INFO] Loaded Lumerical: {LUM_CSV}")
    print(f"[INFO] Loaded HFSS ({hf['source']}): {HFSS_NPZ if hf['source']=='npz' else HFSS_CSV}")
    print(f"[INFO] Using ng = {ng_used}")
    print(f"[INFO] Using L_m = {L_m:.6e} m (source={L_src})")

    # 4) Compute EO response analytically
    eo = compute_eo_response(freqs_Hz=freqs_Hz, Z0=Z0, n_rf=n_rf, ng_used=ng_used, L_m=L_m)

    m_total = eo["m_total"]
    m_total_dB = eo["m_total_dB"]
    m_total_norm = eo["m_total_norm"]
    m_total_norm_dB = eo["m_total_norm_dB"]

    # 5) Estimate -3 dB bandwidth
    f_3db = estimate_3db_bandwidth(freqs_Hz, m_total_norm)
    if f_3db is None:
        print("[RESULT] -3 dB bandwidth: NOT reached within simulated range.")
    else:
        print(f"[RESULT] -3 dB bandwidth: {f_3db/1e9:.3f} GHz")

    # 6) Save CSV (also log L_m used)
    df_out = pd.DataFrame({
        "freq_Hz": np.asarray(freqs_Hz, dtype=float),
        "m_total": np.asarray(m_total, dtype=float),
        "m_total_dB": np.asarray(m_total_dB, dtype=float),
        "m_total_norm": np.asarray(m_total_norm, dtype=float),
        "m_total_norm_dB": np.asarray(m_total_norm_dB, dtype=float),
        "L_m_used": np.full_like(np.asarray(freqs_Hz, dtype=float), L_m, dtype=float),
        "ng_used": np.full_like(np.asarray(freqs_Hz, dtype=float), ng_used, dtype=float),
    })
    df_out.to_csv(OUT_EO_CSV, index=False, float_format="%.12e")
    print(f"[SAVE] EO response CSV: {OUT_EO_CSV}")

    # 7) Plot
    f_GHz = np.asarray(freqs_Hz, dtype=float) / 1e9

    plt.figure(figsize=(10, 6))
    plt.plot(f_GHz, m_total_dB, marker="o", linestyle="-", label="20log10(m_total)")
    plt.xscale("log")
    plt.xlabel("Frequency (GHz)")
    plt.ylabel("Response (dB)")
    plt.title("Electro-Optic Frequency Response")
    plt.grid(True, which="both", linestyle="--")
    if f_3db is not None:
        plt.axvline(f_3db / 1e9, linestyle="--", label=f"-3 dB ~ {f_3db/1e9:.2f} GHz")
    plt.legend()
    plt.show()

    plt.figure(figsize=(10, 6))
    plt.plot(f_GHz, m_total_norm_dB, marker="o", linestyle="-", label="Normalized 20log10(m_total)")
    plt.xscale("log")
    plt.xlabel("Frequency (GHz)")
    plt.ylabel("Normalized Response (dB)")
    plt.title("Electro-Optic Frequency Response (Normalized)")
    plt.grid(True, which="both", linestyle="--")
    plt.axhline(-3.0, linestyle="--", label="-3 dB")
    if f_3db is not None:
        plt.axvline(f_3db / 1e9, linestyle="--", label=f"-3 dB ~ {f_3db/1e9:.2f} GHz")
    plt.legend()
    plt.show()

EO 周波数応答の計算では、進行波型変調器の解析モデルに基づき、電極損失、速度不整合、インピーダンス不整合の影響をすべて含んだ形で変調効率を求めています。RF の伝搬定数は HFSS から得られた複素屈折率を用いて定義され、光波側は群屈折率 ng を用いた近似で扱っています。変調応答は、入力インピーダンスによる係数と、電極上の位相・減衰効果を含む係数の積として定義されており、これらを周波数ごとに評価することで全帯域の EO 応答が得られます。

解析結果

以下に解析結果を示します。今回は位相シフタ長さ1 mmで計算してみました。

計算された EO 応答は、絶対値、dB 表示、最大値で正規化した値、および正規化後の dB 表示として整理されます。さらに、正規化応答が −3 dB に到達する最初の周波数を探索することで、−3 dB 帯域幅を自動的に推定します。
CSV形式でも保存するようになってます。

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

問い合わせ

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

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

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