はじめに

前回の記事(https://qiita.com/Takeshi_Baba/items/be840825c2c344e12d0e)
では光変調器の位相シフタ部分のDC特性の解析事例について紹介しました。

しかしながら、進行波型の光変調器の場合、電極部分の特性インピーダンスや損失、高周波の屈折率も光変調器の周波数応答を決める重要なパラメータになります。
そこで今回は前回の事例の記事の構造をAnsys HFSS上に設定してSパラメータ解析を実施する事例を紹介します。
＊今回は下図の伝送線路をAnsys HFSS上に設定し、PyAEDTを用いてSパラメータ解析を行う方法をご紹介

解析の流れ

解析の流れは以下の通りです
①HFSSのプロジェクトファイルをGUIベースで作成する(変数化したいパラメータも設定しておく。例えば電極の長さやシグナル、グランド電極の幅など)
②作成したプロジェクトファイルを呼び出して、python上で追加設定するスクリプトなどを追加して実行
③出力された結果の表示と追加解析の実行

となります。

解析の実際

HFSSのプロジェクトファイル
こちらが事前に用意した伝送線路のHFSSのプロジェクトファイルです。
左下の変数欄にあるように電極の長さ(L)や、シグナル(WS)、グランドの幅(WG)、シグナルーグランド間のギャップ(gap)などのパラメータを定義しておきます。

PyAEDTのスクリプト
以下がPyAEDTのスクリプトです。
このスクリプトでは、信号線幅やギャップ、電極長といった幾何パラメータを変更しながら自動で解析を行うことを想定して、HFSS による解析から S パラメータ取得、さらに特性インピーダンスや屈折率、損失までを一括に取得できるように設定しています。

# metal_pad.py
# HFSS (PyAEDT) を用いて S パラ→Z0, n_rf を抽出し、NPZ/CSVに保存する完全版（周波数単位の混乱対策込み）
# ★方法①: Sweep一覧の差分から “今作ったSweep” を特定して、そのSweepからデータ取得する版
# ★追加: RF損失（alpha, dB/mm, dB/cm）を算出・保存・プロット

from pyaedt import Desktop, Hfss
import numpy as np
import argparse, os, sys
import pandas as pd
import matplotlib.pyplot as plt


def _normalize_freq_to_hz(freq_raw: np.ndarray) -> np.ndarray:
    """
    PyAEDT/HFSS API から返る primary_sweep_values が、設計単位（例: GHz）で返ってくる場合があるため Hz に正規化する。
    """
    freq_raw = np.asarray(freq_raw, dtype=float)
    if freq_raw.size == 0:
        return freq_raw

    if np.nanmax(freq_raw) < 1e6:
        print("[warn] primary sweep values look NOT in Hz (likely GHz). Converting *1e9 to Hz.")
        return freq_raw * 1e9

    return freq_raw


def _get_sweep_names(hfss: Hfss, setup_name: str) -> list:
    """Setup配下の sweep 名一覧を（文字列で）取得する。"""
    try:
        sweeps = hfss.get_sweeps(setup_name)
        names = [getattr(sw, "name", str(sw)) for sw in sweeps]
        return names
    except Exception as e:
        print(f"[debug] get_sweeps({setup_name}) failed:", e)
        return []


def _pick_newest_created_sweep(before: list, after: list) -> str:
    """
    before/after の差分から「今回増えた sweep」を特定する。
    """
    set_before = set(before)
    added = [x for x in after if x not in set_before]

    if len(added) == 1:
        return added[0]

    if len(added) >= 2:
        print("[warn] multiple sweeps added in this run:", added)
        return after[-1] if after else ""

    print("[warn] no sweep name difference detected. Falling back to the last sweep in 'after'.")
    return after[-1] if after else ""


def _compute_rf_loss_from_gamma(gamma: np.ndarray):
    """
    gamma = alpha + j*beta [1/m] から損失を算出。
    - alpha_np_per_m: [Np/m]
    - amp_db_per_m : 振幅減衰[dB/m] = 8.686 * alpha
    - pwr_db_per_m : 電力減衰[dB/m] = 17.372 * alpha
    - db_per_mm / db_per_cm も返す
    """
    gamma = np.asarray(gamma, dtype=np.complex128)
    alpha = np.real(gamma)  # [Np/m]
    amp_db_per_m = 8.686 * alpha
    pwr_db_per_m = 17.372 * alpha

    amp_db_per_mm = amp_db_per_m / 1e3
    amp_db_per_cm = amp_db_per_m / 1e2
    pwr_db_per_mm = pwr_db_per_m / 1e3
    pwr_db_per_cm = pwr_db_per_m / 1e2

    return alpha, amp_db_per_m, pwr_db_per_m, amp_db_per_mm, pwr_db_per_mm, amp_db_per_cm, pwr_db_per_cm


def run_hfss_and_get_impedance(freqs_target_Hz,
                               design_params,
                               cpw_length_m,
                               aedt_path,
                               design_name="HFSSDesign1",
                               version="2024.2",
                               non_graphical=False,
                               new_desktop=True,
                               close_projects=False,
                               close_desktop=False,
                               use_arccosh=False):
    """
    HFSS を開いて S→Z0 と n_rf を取得する。周波数は必ず Hz に正規化して返す。
    追加: use_arccosh=True の場合は gamma と RF損失も返す（損失はこの場合のみ物理的に意味がある）
    """
    c0 = 3.0e8
    Z_ref = 50.0
    d = float(cpw_length_m)

    setup_name = "Setup1"

    # --- 起動モード ---
    dkt = Desktop(version=version, new_desktop=new_desktop, non_graphical=non_graphical)
    hfss = Hfss(project=aedt_path,
                design=design_name,
                version=version,
                non_graphical=non_graphical,
                new_desktop=False)

    # --- 単位ログ ---
    print(f"[metal_pad] d (for n_rf) = {d:.6e} m (= {d*1e3:.6f} mm)")
    for k, v in design_params.items():
        print(f"[metal_pad] HFSS param {k} = {v} mm (sent as 'mm')")

    # --- 設計変数（mm 前提） ---
    for param, value in design_params.items():
        hfss[param] = f"{value}mm"

    # --- Setup 準備 ---
    if setup_name not in hfss.existing_analysis_setups:
        setup = hfss.create_setup(setup_name)
        setup.props["Frequency"] = "100GHz"
        setup.props["MaximumPasses"] = 15
        setup.update()

    # --- Sweep作成前の一覧 ---
    sweeps_before = _get_sweep_names(hfss, setup_name)
    print("[debug] sweeps BEFORE create:", sweeps_before)

    start_f = float(np.min(freqs_target_Hz))
    stop_f  = float(np.max(freqs_target_Hz))
    npts    = int(len(freqs_target_Hz))
    if npts < 2 or start_f <= 0 or stop_f <= start_f:
        raise ValueError("周波数設定が不正です（点数>=2、start>0、stop>start を満たすこと）")

    print(f"[debug] request sweep (Hz): start={start_f:.6e}, stop={stop_f:.6e}, points={npts}")

    hfss.create_linear_count_sweep(
        setup=setup_name,
        units="Hz",
        start_frequency=start_f,
        stop_frequency=stop_f,
        num_of_freq_points=npts,
        name="Sweep1"
    )

    # --- 解析実行 ---
    hfss.analyze_setup(setup_name)

    # --- 解析後の sweep 一覧 ---
    sweeps_after = _get_sweep_names(hfss, setup_name)
    print("[debug] sweeps AFTER solve:", sweeps_after)

    newest_sweep = _pick_newest_created_sweep(sweeps_before, sweeps_after)
    if newest_sweep == "":
        raise RuntimeError("Could not determine the created sweep name. sweeps_after is empty.")

    setup_sweep_name = f"{setup_name}:{newest_sweep}"
    print("[debug] using setup_sweep_name =", setup_sweep_name)

    # --- 結果取得 ---
    sdata = hfss.post.get_solution_data(
        expressions=["S11", "S12", "S21", "S22"],
        setup_sweep_name=setup_sweep_name
    )

    # primary sweep -> Hz 正規化
    freq_raw = np.array(sdata.primary_sweep_values, dtype=float)
    if freq_raw.size > 0:
        print(f"[debug] freq_raw min/max = {np.min(freq_raw):.6g} / {np.max(freq_raw):.6g}  (len={len(freq_raw)})")
    else:
        print("[debug] freq_raw is empty")

    freq_Hz = _normalize_freq_to_hz(freq_raw)
    if freq_Hz.size > 0:
        print(f"[debug] freq_Hz  min/max = {np.min(freq_Hz):.6e} / {np.max(freq_Hz):.6e}  (len={len(freq_Hz)})")

    # S
    S11 = np.array(sdata.data_real("S11")) + 1j * np.array(sdata.data_imag("S11"))
    S12 = np.array(sdata.data_real("S12")) + 1j * np.array(sdata.data_imag("S12"))
    S21 = np.array(sdata.data_real("S21")) + 1j * np.array(sdata.data_imag("S21"))
    S22 = np.array(sdata.data_real("S22")) + 1j * np.array(sdata.data_imag("S22"))

    # Z0
    Z0 = Z_ref * np.sqrt(((1 + S11)**2 - S21**2) / ((1 - S11)**2 - S21**2))

    # n_rf & gamma
    k0 = 2 * np.pi * freq_Hz / c0
    cosh_arg = (1 - S11**2 + S21**2) / (2 * S21)

    gamma = None
    rf_loss = None

    if use_arccosh:
        gamma = np.arccosh(cosh_arg) / d  # [1/m]
        n_rf = gamma / k0                 # 複素 n_rf
        rf_loss = _compute_rf_loss_from_gamma(gamma)

        # 確定ログ（損失）
        alpha, amp_db_m, pwr_db_m, amp_db_mm, pwr_db_mm, amp_db_cm, pwr_db_cm = rf_loss
        print(f"[debug] alpha (Np/m) min/max = {np.min(alpha):.6e} / {np.max(alpha):.6e}")
        print(f"[debug] loss (amp dB/mm) min/max = {np.min(amp_db_mm):.6e} / {np.max(amp_db_mm):.6e}")
        print(f"[debug] loss (pwr dB/mm) min/max = {np.min(pwr_db_mm):.6e} / {np.max(pwr_db_mm):.6e}")

    else:
        cos_arg = np.real(cosh_arg)
        cos_arg = np.clip(cos_arg, -1.0, 1.0)
        n_rf = (1.0 / (k0 * d)) * np.arccos(cos_arg)  # 実数
        print("[info] use_arccosh=False -> gamma/loss not available (phase-only approximation).")

    # 健全性
    if len(freq_Hz) == 0:
        raise RuntimeError("HFSSの結果が取得できませんでした（freq が空）")
    if np.any(~np.isfinite(freq_Hz)):
        raise RuntimeError("freq に非有限値が含まれています")

    hfss.release_desktop(close_projects=close_projects, close_desktop=close_desktop)
    return freq_Hz, Z0, n_rf, gamma, rf_loss


# ===== CLI / Spyder 対応 =====
if __name__ == "__main__":

    def _running_in_spyder():
        return ("spyder_kernels" in sys.modules) or ("SPYDER_PID" in os.environ)

    parser = argparse.ArgumentParser(description="Run HFSS via PyAEDT and extract Z0, n_rf, and RF loss.")
    parser.add_argument("--aedt", type=str,
                        default=r"C:\Users\babatake\Desktop\test_CPW2.aedt",
                        help="Path to the AEDT project (.aedt)")
    parser.add_argument("--design", type=str, default="HFSSDesign1", help="Design name")
    parser.add_argument("--version", type=str, default="2024.2", help="AEDT version")

    parser.add_argument("--non-graphical", action="store_true", help="Headless mode")
    parser.add_argument("--attach", action="store_true",
                        help="Attach to existing AEDT instead of launching new (recommended for Spyder)")
    parser.add_argument("--close-desktop", action="store_true",
                        help="Close AEDT at the end (default: keep open)")

    parser.add_argument("--start", type=float, default=1e9, help="Start freq [Hz]")
    parser.add_argument("--stop", type=float, default=1e11, help="Stop freq [Hz]")
    parser.add_argument("--points", type=int, default=21, help="Number of points")

    parser.add_argument("--length_mm", type=float, default=0.1, help="CPW length d [mm] (for n_rf)")
    parser.add_argument("--WG_mm", type=float, default=0.1, help="WG width [mm]")
    parser.add_argument("--WS_mm", type=float, default=0.01, help="Signal width [mm]")
    parser.add_argument("--gap_mm", type=float, default=0.01, help="Gap [mm]")
    parser.add_argument("--L_mm", type=float, default=0.1, help="HFSS param L [mm]")

    parser.add_argument("--save", type=str, default="hfss_cpw_results.npz", help="Output NPZ path")
    parser.add_argument("--csv", type=str, default="hfss_cpw_results.csv", help="Output CSV path")
    parser.add_argument("--plot", action="store_true", help="Plot results")

    parser.add_argument("--use-arccosh", action="store_true",
                        help="Use arccosh to compute complex n_rf (required for RF loss extraction)")

    args = parser.parse_args()
    args.plot = True  # ★要望どおり残す
    args.use_arccosh = True

    print(f"[debug] args.start/stop/points = {args.start} / {args.stop} / {args.points}")

    in_spyder = _running_in_spyder()
    new_desktop = not (args.attach or in_spyder)
    non_graphical = args.non_graphical and (not in_spyder)
    close_desktop = args.close_desktop and (not in_spyder)

    freqs_target_Hz = np.linspace(args.start, args.stop, args.points)
    design_params = {"WG": args.WG_mm, "WS": args.WS_mm, "gap": args.gap_mm, "L": args.L_mm}
    d_m = args.length_mm * 1e-3

    print(f"AEDT: {args.aedt}")
    print(f"Design: {args.design} | Version: {args.version}")
    print(f"Mode: {'attach(existing GUI)' if not new_desktop else 'launch new'} | "
          f"{'GUI' if not non_graphical else 'headless'} | "
          f"{'keep GUI' if not close_desktop else 'close GUI'}")
    print(f"Freq request: {args.start:.3e} Hz → {args.stop:.3e} Hz | points={args.points}")

    print("Running sweep ...")
    freq_Hz, Z0, n_rf, gamma, rf_loss = run_hfss_and_get_impedance(
        freqs_target_Hz=freqs_target_Hz,
        design_params=design_params,
        cpw_length_m=d_m,
        aedt_path=args.aedt,
        design_name=args.design,
        version=args.version,
        non_graphical=non_graphical,
        new_desktop=new_desktop,
        close_projects=False,
        close_desktop=close_desktop,
        use_arccosh=args.use_arccosh
    )
    print("Done.")

    # NPZ保存（gamma/lossも入れる）
    np.savez(args.save,
             freq_Hz=freq_Hz,
             Z0=Z0,
             n_rf=n_rf,
             gamma=gamma if gamma is not None else np.array([]),
             design_params=design_params,
             d_m=d_m,
             use_arccosh=args.use_arccosh)
    print(f"[NPZ] saved: {os.path.abspath(args.save)}")

    # CSV保存
    freq_Hz = np.asarray(freq_Hz, dtype=np.float64)

    out = {
        "freq_Hz":   freq_Hz,
        "Z0_real":   np.asarray(np.real(Z0), dtype=np.float64),
        "Z0_imag":   np.asarray(np.imag(Z0), dtype=np.float64),
        "n_rf_real": np.asarray(np.real(n_rf), dtype=np.float64),
        "n_rf_imag": np.asarray(np.imag(n_rf), dtype=np.float64),
    }

    # RF損失（use_arccosh=True の時のみ追加）
    if rf_loss is not None:
        alpha, amp_db_m, pwr_db_m, amp_db_mm, pwr_db_mm, amp_db_cm, pwr_db_cm = rf_loss
        out.update({
            "alpha_Np_per_m": np.asarray(alpha, dtype=np.float64),
            "loss_amp_dB_per_m": np.asarray(amp_db_m, dtype=np.float64),
            "loss_pwr_dB_per_m": np.asarray(pwr_db_m, dtype=np.float64),
            "loss_amp_dB_per_mm": np.asarray(amp_db_mm, dtype=np.float64),
            "loss_pwr_dB_per_mm": np.asarray(pwr_db_mm, dtype=np.float64),
            "loss_amp_dB_per_cm": np.asarray(amp_db_cm, dtype=np.float64),
            "loss_pwr_dB_per_cm": np.asarray(pwr_db_cm, dtype=np.float64),
        })

    df = pd.DataFrame(out)
    df.to_csv(args.csv, index=False, float_format="%.12e")
    print(f"[CSV] saved: {os.path.abspath(args.csv)}")

    # プロット（表示用に GHz）
    if args.plot:
        f_GHz = freq_Hz / 1e9

        # (1) |Z0|
        plt.figure(figsize=(8, 5))
        plt.plot(f_GHz, np.abs(Z0), label="|Z0| [Ohm]")
        plt.xscale("log")
        plt.xlabel("Frequency [GHz]")
        plt.ylabel("|Z0| [Ohm]")
        plt.grid(True, which="both")
        plt.legend()
        plt.tight_layout()
        plt.show()

        # (2) Re/Im(Z0)
        plt.figure(figsize=(8, 5))
        plt.plot(f_GHz, np.real(Z0), label="Re(Z0) [Ohm]")
        plt.plot(f_GHz, np.imag(Z0), label="Im(Z0) [Ohm]")
        plt.xscale("log")
        plt.xlabel("Frequency [GHz]")
        plt.ylabel("Z0 [Ohm]")
        plt.grid(True, which="both")
        plt.legend()
        plt.tight_layout()
        plt.show()

        # (3) n_rf
        plt.figure(figsize=(8, 5))
        if args.use_arccosh:
            plt.plot(f_GHz, np.real(n_rf), label="Re(n_rf) [-]")
            plt.plot(f_GHz, np.imag(n_rf), label="Im(n_rf) [-]")
        else:
            plt.plot(f_GHz, np.real(n_rf), label="n_rf (real) [-]")
        plt.xscale("log")
        plt.xlabel("Frequency [GHz]")
        plt.ylabel("n_rf [-]")
        plt.grid(True, which="both")
        plt.legend()
        plt.tight_layout()
        plt.show()

        # (4) RF損失（dB/mm, dB/cm）
        if rf_loss is not None:
            _, _, _, amp_db_mm, pwr_db_mm, amp_db_cm, pwr_db_cm = rf_loss

            plt.figure(figsize=(8, 5))
            plt.plot(f_GHz, amp_db_mm, label="Loss (amp) [dB/mm]")
            plt.plot(f_GHz, pwr_db_mm, label="Loss (power) [dB/mm]")
            plt.xscale("log")
            plt.xlabel("Frequency [GHz]")
            plt.ylabel("Loss [dB/mm]")
            plt.grid(True, which="both")
            plt.legend()
            plt.tight_layout()
            plt.show()

            plt.figure(figsize=(8, 5))
            plt.plot(f_GHz, amp_db_cm, label="Loss (amp) [dB/cm]")
            plt.plot(f_GHz, pwr_db_cm, label="Loss (power) [dB/cm]")
            plt.xscale("log")
            plt.xlabel("Frequency [GHz]")
            plt.ylabel("Loss [dB/cm]")
            plt.grid(True, which="both")
            plt.legend()
            plt.tight_layout()
            plt.show()
        else:
            print("[info] RF loss plot skipped (set --use-arccosh to enable complex gamma/loss).")

解析後に取得した S パラメータ（S11、S21 など）からは、進行波型伝送線路の理論式に基づいて特性インピーダンス(Z0)を求めてます。を算出しています。基準インピーダンスは 50 Ω としており、一般的な二端子ネットワークの関係式をそのまま用いています。また、RF 有効屈折率nrfは線路長dと自由空間波数k0から求めています。計算結果は NumPy の NPZ 形式および CSV 形式で保存されるため、後段で Lumerical や Python を用いた回路解析、光学解析、あるいは電気・光学連成シミュレーションにそのまま利用することができます。

解析結果

こちらが解析結果になります。

これらの結果は、上記のスクリプトを実施すると自動でcsvで保存されるようになっていますので、それらの解析結果を利用して、変調器の光F特などを求めるのに有効です。

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

問い合わせ

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

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

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