TES の fixed-tone マイクロ波読み出しを Python で直感的に理解する -- 共振器 + I/Q + MW-Mux のミニシミュレータ入門

0. はじめに

TES（Transition-Edge Sensor）のマイクロ波読み出し（いわゆる MW-Mux / uMux）では、

• 数 GHz 帯の超伝導共振器をたくさん並べて
• それぞれに少しずつ異なる共振周波数 $f_0$ を割り当て
• 各共振器に対応した 固定周波数トーン（fixed-tone） を当てっぱなしにして
• X 線が TES に当たったとき、TES ↔ コイル ↔ SQUID ↔ 共振器の結合を通じて生じる
  共振周波数のズレ を
• マイクロ波の 透過振幅 / 位相（I/Q）の変化 として読み出す

という方式が使われています。……と文章で読んでも、

「結局、何がどう変わっていて、それをどうやって読んでいるのか？」

がイメージしにくいことが多いと思います。

そこで本記事では、Python + matplotlib だけ を使って、

• 単一共振器に対する fixed-tone readout
• そのときの I/Q（複素 $S_{21}$）の動き
• 複数共振器を並べた MW-Mux 的シミュレーション

を行うミニコードを題材に、

「TES の fixed-tone マイクロ波読み出しの物理的イメージ」

を、数式と図を通して直感的に理解することを目指します。

S21の数学的な意味については下記の記事を参考にしてください。

0.1 ラジオとの違い

ラジオは、空間に存在する電磁波を「受ける（受信する）」ことで情報を得ます。
つまり基本は 受動的（passive） で、受信機は「入ってきた波」を増幅・復調しているだけです。

一方で、TES のマイクロ波読み出し（MW-Mux / uMux）は発想が逆で、

こちらから狙った周波数のマイクロ波トーンを能動的（active）に注入し，
回路を通って戻ってきた波の “振幅と位相” の変化を診断する

という 計測器（プローブ） としての読み出しです。

• 入れる：既知の連続波（fixed-tone, probe tone）
• 見る：戻ってくる波の 複素比（透過なら $S_{21}$，反射なら $S_{11}$）
• 変わる理由：TES イベントで共振器の 共振周波数 $f_\mathrm{res}$ がわずかに動く
• その結果：固定周波数 $f_\mathrm{probe}$ で見た $S_{21}(t)$ の 位相・振幅$I/Q$ が時間的に動く

要するに MW-Mux は、

「共振器という“周波数センサー”を、マイクロ波トーンで常時プロービングして、
共振周波数のズレを I/Q の動きに写像して読む」

という方式です。

この “能動的に入れて、返りを診断する” という点が、ラジオ（受動受信）と本質的に異なります。

1. TES のマイクロ波読み出しの「最小モデル」

1.1 何が起きているのかを 1 行で言うと

X 線が TES に入射すると、

1. TES が温まり、抵抗やインダクタンスが変化する
2. TES に結合した超伝導共振器の 共振周波数 $f_\mathrm{res}$ が一時的にずれる
3. 共振周波数のズレにより、固定周波数トーンで見ている $S_{21}$ の振幅・位相が変化する
4. その時間波形から、エネルギーや到来時刻を復元する

という流れになっています。

本記事のシミュレータでは、TES・RF-SQUID・マイクロ波線路などの回路詳細は追わず、
「共振周波数が時間的にずれる単純な共振器」 として振る舞いだけを取り出します。

1.2 ノッチ型共振器の周波数応答

実際の uMux 回路では、通過帯（$|S_{21}|\approx 1$）の中に
「共振周波数でノッチ（減衰）が落ちる」タイプの応答が典型的です。

補足説明: S パラメータと S_21 とは？

マイクロ波回路では、電圧や電流そのものではなく、

• 「ポート 1 からどのくらい入ってきて」
• 「ポート 2 からどのくらい出ていくか」

を表す 散乱パラメータ (S パラメータ) をよく使います。
2 ポート回路だと

• $S_{11}$ : ポート 1 に入れた信号のうち、ポート 1 に反射して戻ってくる成分
• $S_{21}$ : ポート 1 に入れた信号のうち、ポート 2 から通過して出てくる成分

というイメージです（$S_{12}, S_{22}$ も同様）。

$S_{21}$ は 複素数 になっていて、

• 大きさ $|S_{21}|$ ：「どのくらいの振幅で通過しているか」（減衰の量）
• 位相 $\arg S_{21}$：通過する過程でどれだけ位相が回転したか

を表しています。

ネットワークアナライザ（VNA）が測っているのはまさにこの $S_{21}$ で、

• 実部 $\Re(S_{21})$ → I 成分
• 虚部 $\Im(S_{21})$ → Q 成分

として読み出したものが、実験でいう「I/Q データ」です。

---

なぜ S_21(f) と書いてよいのか？

ここで $S_{21}$ は 周波数 $f$ の関数 $S_{21}(f)$（＝一価なマップ） だとみなしています。

• 回路が 線形・時間不変 (LTI:Linear Time-Invariant) で、
• 入力に 1 つの周波数 $f$ の正弦波（$e^{i 2\pi f t}$）だけを入れる

と、出力も同じ周波数 $f$ の正弦波で出てきて、その 振幅と位相の比 が
「周波数 $f$ における $S_{21}(f)$」として一意に定まります。
つまり、
正弦波を入れたら、同じ周波数の正弦波で返してくれる “お行儀の良い系”
ということ。

実際の TES–uMux 系では、TES がゆっくり変化することで共振周波数 $f_\mathrm{res}(t)$ が動きますが、

• マイクロ波の周期（数100ピコ秒）に比べて
• TES パルスの時間スケール（ミリ秒）が圧倒的に長い

ため、「その瞬間ごとにはほぼ定常な回路」とみなして

各瞬間の回路に対する $S_{21}(f)$ が定義されていて、
TES の揺らぎはそのパラメータ（主に $f_\mathrm{res}$）をゆっくり変えているだけ

という 準定常近似 を使うことができます。

この記事では、この近似のもとで
共振器を通過するマイクロ波の「透過係数」 として $S_{21}(f)$ を見ており、
「TES が揺らいだときに $S_{21}$ がどう変わるか？」という視点で話を進めています。

この記事のコードでは、その一番単純なモデルとして

S_{21}(f) = 1 - \frac{1}{1 + i x},
\qquad
x = \frac{f - f_\mathrm{res}}{\gamma},

を使います。ここで

\gamma = \frac{f_\mathrm{res}}{2Q}

は半値半幅（HWHM）で、Q 値が大きいほど共振はシャープになります。

このとき、$S_{21}$ を実部・虚部・振幅に分解すると

\frac{1}{1+i x}
= \frac{1 - i x}{1 + x^2}

なので

\begin{aligned}
S_{21}(f)
&= 1 - \frac{1 - i x}{1 + x^2} \\
&= \left(1 - \frac{1}{1 + x^2}\right)
 + i \frac{x}{1 + x^2}.
\end{aligned}

したがって、$S_{21}$ を実部・虚部・振幅、はそれぞれ、

\boxed{
\begin{aligned}
I(x) &= \Re(S_{21}) = \frac{x^2}{1 + x^2}, \\
Q(x) &= \Im(S_{21}) = \frac{x}{1 + x^2}, \\
|S_{21}(x)| &= \sqrt{I(x)^2 + Q(x)^2}.
\end{aligned}
}

と書けます。

• $f$ が共振から十分離れている ($|x|\gg 1$) とき
  $\displaystyle \frac{1}{1 + i x} \approx 0$ なので $S_{21} \approx 1$、ほぼ全通過。
• 逆に $f = f_\mathrm{res}（x=0)$では
  $\displaystyle S_{21} = 1 - 1 = 0$ となり、理想的には完全なノッチになります
  （実機では損失や結合条件のため 0 までは落ちません）。

ここまでが「静的な周波数応答」の話です。

1.3 X 線イベントによる共振周波数の時間変化

TES に X 線が入射した瞬間を $t = t_\mathrm{event}$ として、共振周波数が

• それ以前：
  $\displaystyle f_\mathrm{res}(t) = f_0$
• それ以降：
  $\displaystyle f_\mathrm{res}(t) = f_0 + \Delta f_\mathrm{max}\exp!\left(-\frac{t - t_\mathrm{event}}{\tau}\right)$

のように動く、と単純化しておきます：

f_\mathrm{res}(t) =
\begin{cases}
f_0, & t < t_\mathrm{event},\\[4pt]
f_0 + \Delta f_\mathrm{max}
 \exp\!\left(-\dfrac{t - t_\mathrm{event}}{\tau}\right),
 & t \ge t_\mathrm{event}.
\end{cases}

ここで

• $f_0$：基本の共振周波数（例：5 GHz）
• $\Delta f_\mathrm{max}$：X 線イベントでの最大シフト（例：10 MHz）
• $\tau$：緩和時間（例：数〜数十 ms）

となります。

大ざっぱな物理イメージとしては

• $\Delta f_\mathrm{max}$：TES パルスの「高さ」 → エネルギーに比例
• $\tau$：TES の熱が抜ける時間 → パルス幅

と対応づけて考えることができます。

補足説明: この S21 モデルの元になっている基礎方程式は？

もう少し「原理寄り」に見ると、ここで使っている

S_{21}(f) = 1 - \frac{1}{1 + i x},\quad
x = \frac{f - f_\mathrm{res}}{\gamma}

という形は、1 自由度の RLC 共振器に
「入力ポート」「出力ポート」を弱く結合したときに現れる

• 1 極 (single pole) の周波数応答
• ノッチ型の結合

を、できるだけ簡単な形に潰したものになっています。

元の微分方程式レベルでは、例えば

L\frac{\mathrm{d}^2 q}{\mathrm{d}t^2}
+ R\frac{\mathrm{d} q}{\mathrm{d}t}
+ \frac{1}{C} q
= V_\mathrm{drive}(t)

のような RLC 共振回路の運動方程式から出発できます。
ここで $q$ はコンデンサに溜まる電荷、$V_\mathrm{drive}$ は外から加える駆動電圧です。

これをサイン波駆動

V_\mathrm{drive}(t) = \Re\{ \tilde{V}\, e^{i\omega t} \}

とみなしてフーリエ領域に移すと、

Z(\omega)
= R + i\left(\omega L - \frac{1}{\omega C}\right)

というインピーダンスになり、

• 共振角周波数 $\omega_0 = 1/\sqrt{LC}$
• Q 値 $Q = \omega_0 L / R$

といったおなじみの関係が出てきます。

この RLC をマイクロ波ラインに結合した 2 ポート回路として解くと、

\frac{1}{1 + i x},\quad
x = \frac{f - f_\mathrm{res}}{f_\mathrm{res}/(2Q)}

のような「1 極ローレンツ型応答」が現れます。
さらに、ポートへの結合のさせ方（直列結合か並列結合か、結合の強さなど）によって

• ピーク型
• ノッチ型

のどちらの $S_{21}$ も作ることができます。

この記事では、uMux でよく使う ノッチ型のケースを最も単純化した形として

S_{21}(f) = 1 - \frac{1}{1 + i x}

を採用し、「TES が変えるのは $f_\mathrm{res}$ 側だけ」という
縮約された有効モデルとして扱っています。
厳密な結合定数や $Q_\mathrm{c}, Q_\mathrm{i}$ の議論は脇に置きつつも、

• 「共振周波数が動けば $S_{21}$ が動く」
• 「Q が大きいほど、少しの周波数シフトで大きく $S_{21}$ が変わる」

という本質だけを抽出した形だと思っていただければ十分です。

2. fixed-tone microwave readout の考え方

2.1 Fixed-tone とは？

理想的には、共振器の周波数を高速に掃引して

「その瞬間の共振周波数 $f_\mathrm{res}(t)$ を直接測る」

という方法も考えられます。

しかし実際の uMux では、

• 数百〜数千チャンネルの TES を同時に読む必要があり
• DAC・ADC の本数、FPGA 資源、帯域などの制約が厳しい

ため、実用的な方式として

各共振器ごとに 1 本の 固定周波数トーン（probe tone） を当てっぱなしにし、
そのトーンの振幅・位相の時間変化だけを観測する

という fixed-tone readout が使われます。

2.2 なぜ fixed-tone で読めるのか？（数式で確認）

固定トーンの周波数を $f_\mathrm{probe}$ とすると、時間 $t$ における「ずれ」は

\Delta f(t) = f_\mathrm{probe} - f_\mathrm{res}(t)

で与えられます。

これを HWHM で規格化した無次元量

x(t) = \frac{f_\mathrm{probe} - f_\mathrm{res}(t)}{\gamma}

を通じて、$S_{21}(t) = S_{21}(x(t))$ が 時間とともに動く複素数 になります。

つまり fixed-tone readout は

「周波数シフト $f_\mathrm{res}(t)$ の時間変化を
単一周波数で見た複素数 $S_{21}(t)$ の軌跡に写像している」

と理解できます。

実際、"I/Q" を測っている、という意味は、

I(t) = \Re(S_{21}(t)), \qquad Q(t) = \Im(S_{21}(t))

であり、これを適切な方向に射影することで TES パルス波形に変換したり、
さらに最適フィルタをかけてエネルギーを推定したりします。

補足説明 : なぜ「周波数を掃引せず fixed-tone で十分なのか？」

直感的には、

「共振周波数が動くなら、その都度 VNA のように周波数を掃引して
$f_\mathrm{res}(t)$ を直接測ればいいのでは？」

と思いたくなりますが、実際の uMux では 数百〜数千チャンネル を
1 本のマイクロ波ラインで読む必要があります。

• それぞれの共振器ごとに掃引していたら時間が足りない
• DAC / ADC / FPGA のリソースや帯域がすぐに限界を迎える

といった事情から、

「各共振器ごとに 1 本の固定トーンだけを当てっぱなしにして、
そのトーンの I/Q の時間変化だけを読む」

という方式が合理的になります。

しかも、周波数シフトが十分小さいときには、
$S_{21}$ は「ほぼ線形なセンサー」として振る舞ってくれます。

固定トーンの周波数を $f_\mathrm{probe}$ とすると、

S_{21}(t)
= S_{21}\bigl(f_\mathrm{probe} - \Delta f(t)\bigr)
\simeq S_{21}(f_\mathrm{probe})
 - \left.\frac{\mathrm{d}S_{21}}{\mathrm{d}f}\right|_{f_\mathrm{probe}}
   \Delta f(t).

ここで

• $\Delta f(t) = f_\mathrm{probe} - f_\mathrm{res}(t)$
• $\dfrac{\mathrm{d}S_{21}}{\mathrm{d}f}$ は、固定トーンの周波数での「勾配」

です。$\Delta f(t)$ が小さいなら、

• $I(t) = \Re[S_{21}(t)]$
• $Q(t) = \Im[S_{21}(t)]$

も $\Delta f(t)$ にほぼ比例して変化します。
つまり fixed-tone readout は、

「共振曲線のある一点の傾きをうまく利用して、
周波数シフトを I/Q の線形変化として読む方式」

とみなすことができます。

どの周波数をプローブに選ぶか（共振のどの辺りにトーンを置くか）によって、

• 勾配の大きさ（感度）
• 応答が主に I に乗るか、Q に乗るか
• 線形性がどこまで保たれるか

が変わってくるので、実機ではこの「プローブ周波数の取り方」も
重要な設計パラメータになります。

X線のように早い信号の場合は難しいですが、電波のようにゆっくりとパワーを検出する場合は、tone tracking という手法で、パワーが小さいところに tone の probe を動かす方式も研究が進められています。

• SLAC microresonator RF (SMuRF) electronics: A tone-tracking readout system for superconducting microwave resonator arrays, Cyndia Yu et al., Rev. Sci. Instrum. 94, 014712 (2023) link

3. シミュレータの全体構成

この記事の Python コードは、ざっくり

1. 共振器の複素伝達関数 resonator_s21_complex
   （ノッチ型 $S_{21}$ の最小モデル）
2. 単一共振器の時間発展 + I/Q 可視化
   （fixed-tone readout の I/Q ベクトルの動きを見る）
3. 複数共振器（MW-Mux 風）の時間発展可視化
   （周波数軸に共振器を並べて multi-tone で読む）

という 3 層構造になっています。

コード全文は記事末尾の <details> に置き、ここではポイントごとに抜粋して説明します。

4. ノッチ型共振器の複素 S_21 モデル

def resonator_s21_complex(f, f_res, Q):
    """
    単純化された 1ポール共振器の複素伝達関数 S21(f) を返す (ノッチ型)。
    """
    gamma = f_res / (2.0 * Q)  # HWHM [Hz]
    x = (f - f_res) / gamma
    S21 = 1.0 - 1.0 / (1.0 + 1j * x)
    return S21

4.1 I–Q 平面では円弧を描く

先ほどの式

I(x) = \frac{x^2}{1 + x^2}, \qquad
Q(x) = \frac{x}{1 + x^2}

から、$x$ を消去して I–Q の関係式を作ってみます。

まず $\displaystyle I = \frac{x^2}{1 + x^2} \Rightarrow
x^2 = \frac{I}{1-I}$。

これを $Q^2$ に代入すると

Q^2
= \frac{x^2}{(1 + x^2)^2}
= \frac{I}{1-I} \left(1-I\right)^2
= I(1-I).

両辺を整理すると

Q^2 + (I - 1/2)^2 = (1/2)^2.

つまり

I–Q 平面上の軌道は、
中心 ((I,Q)=(0.5, 0))、半径 (0.5) の円

になります。

X 線イベントの前後で $x(t)$ が変化することで、
I–Q 平面上を「円弧として」ぐるっと動き、時間とともに元の点へ戻っていく、
という幾何学的なイメージがここから理解できます。

（コード中の「I–Q trajectory at fixed tone」という図がまさにこの円弧です）

5. 単一共振器の fixed-tone + I/Q シミュレーション

5.1 時間発展部分：simulate_single_resonator

def simulate_single_resonator(
    f0=5.0e9,
    Q=300.0,
    df_max=10.0e6,
    tau_ms=10.0,
    t_total_ms=50.0,
    dt_us=20.0,
    t_event_ms=5.0,
    probe_offset_MHz=0.0,
):
    ...
    dt_s = dt_us * 1.0e-6
    t_total_s = t_total_ms * 1.0e-3
    t_event_s = t_event_ms * 1.0e-3
    tau_s = tau_ms * 1.0e-3

    t_s = np.arange(0.0, t_total_s, dt_s)
    t_ms = t_s * 1.0e3

    # 共振周波数の時間変化 f_res(t)
    f_res_Hz = np.empty_like(t_s)
    for i, t in enumerate(t_s):
        if t < t_event_s:
            f_res_Hz[i] = f0
        else:
            f_res_Hz[i] = f0 + df_max * np.exp(-(t - t_event_s) / tau_s)

    # fixed-tone
    f_probe_Hz = f0 + probe_offset_MHz * 1.0e6

    # S21(t) を計算
    S21_t = resonator_s21_complex(f_probe_Hz, f_res_Hz, Q)
    ...

• ここでは 1 本の共振器 + 1 本のトーン の世界だけを考えています。
• df_max と tau_ms が、TES パルスの大きさと時間幅の役割を担っています。
• probe_offset_MHz を変えることで、「共振中心からどれだけ detune した周波数で読むか」を実験できます。

5.2 可視化：周波数シフト・I/Q・|S21|・I–Q

plot_single_resonator では、

1. (左上) 時間 vs $\Delta f_\mathrm{res}(t)$
2. (右上) 時間 vs $|S_{21}(t)|$（線形スケール）
3. (左下) I–Q 平面上の軌跡
4. (右下) 静的な共振曲線（周波数 vs $|S_{21}(f)|$）と probe の位置

をまとめて表示しています。

single_resonator_linear.png

次に、dB表示をして、dBの感覚を養う図をプロットしてみましょう。

補足説明: dB 表示って何をしているのか

$S_{21}$ の振幅 $|S_{21}|$ は 0〜1 の間の 線形な比 ですが、現場ではよく

S_{21~\mathrm{dB}} = 20\log_{10}\bigl(|S_{21}|\bigr)

という dB (デシベル) 表示 を使います。

• $|S_{21}| = 1$ なら $0~\mathrm{dB}$
  → 「損失ゼロ（理想的には）」という意味
• $|S_{21}| = 0.5$ なら $20\log_{10}(0.5) \approx -6~\mathrm{dB}$
  → 振幅が半分（パワーは 1/4）
• $|S_{21}| = 1/\sqrt{2}$ なら $-3~\mathrm{dB}$
  → パワーがちょうど半分になる点（いわゆる「-3 dB 点」）

となります。

なぜ 20 なのか？というと、

• パワー比は $10\log_{10}(P/P_0)$ [dB]
• 振幅は「電圧」や「電流」で、パワーは振幅の 2 乗に比例

なので、

10\log_{10}\left(\frac{P}{P_0}\right)
= 10\log_{10}\left(\frac{|V|^2}{|V_0|^2}\right)
= 20\log_{10}\left(\frac{|V|}{|V_0|}\right)

となり、振幅を dB 表示するときは 20 がつきます。

dB 表示のメリットは、

• 小さな変化（-0.1 dB など）と大きな変化（-30 dB など）を
  1 本のスケールにまとめて見やすくできる
• 共振の「深さ」や「幅」を、足し算・引き算で扱える

といった点です。
このあと出てくる dB プロットも、線形スケールでは見えにくい細かい差を
人間の目で捉えやすくするためのものだと思ってください。

single_resonator_db.png

さらに plot_single_resonator_time_series では、

• 時間 vs $I(t)$
• 時間 vs $Q(t)$
• 時間 vs $|S_{21}(t)|$（左軸） + 時間 vs $\Delta f_\mathrm{res}(t)$（右軸）

を 3 段で描きます。

３段目は、
左軸（青）: 実際に readout が見ている量 $|S_{21}(t)|$
右軸（オレンジ）：その原因となっている共振周波数の変化 $\Delta f_\mathrm{res}(t)$

を表示し、

• 「原因：周波数シフト」
• 「結果：振幅変化」

の因果関係が一目で分かるように工夫してみました。

6. 複数共振器で MW-Mux をイメージする

6.1 周波数軸上に共振器を並べる

idx_offset = np.arange(n_res) - (n_res - 1) / 2.0
f0_list = f0_center + idx_offset * spacing_MHz * 1.0e6

ここでは、例として n_res=4 の共振器を

• 中心周波数 f0_center をはさんで
• spacing_MHz 間隔で

等間隔に並べています。

イメージとしては、周波数軸に

... 4.925 GHz, 4.975 GHz, 5.025 GHz, 5.075 GHz ...

のように TES チャンネルが並んでいる状況です。

6.2 各共振器に 1 本ずつトーンを当てる

f_probe_k = f0_k
S21_k = resonator_s21_complex(f_probe_k, f_res_k, Q)

各共振器は「自分の共振周波数」でロックしており、
そのトーンの I/Q だけを時間的に見ています。

6.3 1 本だけが X 線ヒットを受ける

if k == hit_index and t >= t_event_s:
    f_res_k[i] = f0_k + df_max * np.exp(-(t - t_event_s) / tau_s)
else:
    f_res_k[i] = f0_k

• hit_index で指定した共振器だけが周波数シフトを経験し、
• 他の共振器はずっとフラットのまま

になるようにしています。

その結果、

• plot_multi_resonators の time-domain プロットでは
  「ある 1 本のトレースだけがパルス状に変化」 し、
• 残りのトレースはほぼ水平線になります。

これは

「周波数軸に並んだ TES チャンネルのうち、どれに X 線が当たったかを
単一のマイクロ波線路で見分ける」

という MW-Mux の基本原理そのものです。

multi_resonator_time_linear.png

6.4 静的な共振曲線をまとめて描く

同じ関数 plot_multi_resonators の後半では、

• 横軸を $f - f_\mathrm{center}$ [MHz]
• 縦軸を $|S_{21}(f)|$ または $20\log_{10}|S_{21}(f)|$

として、複数のノッチ共振を一括でプロットしています。

• 各共振器に対応する縦線（probe tone）の位置も表示しているので、
• 「周波数ドメインでのチャンネル配置」を視覚的に確認できます。

multi_resonator_static_linear.png

dB表示をすると、下記になる。

multi_resonator_static_db.png

7. このシミュレータから学べること

このコードの重要ポイントを整理すると：

1. TES が温まる → 共振器の共振周波数が少しずれる
2. fixed-tone では、その周波数の動きを
   単一周波数で見た $|S_{21}|$ と I/Q の変化として読む
3. I–Q 平面で見ると、円弧をなぞるベクトルの動きとして可視化できる
4. 複数共振器を並べれば、周波数軸に TES が並んだマトリクス読み出しになる

という 物理のストーリー全体を、手元の PC で「動く図」として体感できることです。

さらに発展的なトピックとしては、

• ノイズ項を足して S/N と $Q$、$\Delta f_\mathrm{max}$ の関係を見る
• I–Q 軌跡をある方向に射影して パルス波形を作り、
  optimal filter を適用する
• 同時に 2 個のイベントを入れて パイルアップの挙動を見る

なども、このコードを少し拡張するだけで試すことができます。

8. 実行方法と「いじって遊ぶ」ためのヒント

8.1 実行方法

ファイル名を fixed_tone_mw_readout_sim.py とすると、

python fixed_tone_mw_readout_sim.py

で実行できます。

main() では

• demo_single_resonator()
• demo_multi_resonators()

を順に呼んでいるので、

1. 単一共振器の fixed-tone + I/Q の図
2. 複数共振器（MW-Mux 風）の図

の順に PNG ファイルがカレントディレクトリに保存されます。

Google Colab で使う場合は、main() の定義はそのままにしておきつつ、

# 単一共振器だけ見たいとき
demo_single_resonator()

# すべてのデモを一気に走らせたいとき
main()

のように、下のセルから明示的に呼ぶ形にすると混乱しにくいと思います。

8.2 触ってみると面白いパラメータ

単一共振器側（demo_single_resonator() 内）では特に：

• probe_offset_MHz

  • 0 → 5 → 10 と変えて I–Q 円のどの部分をなぞるかを比較

• Q

  • 100, 300, 1000 などに変え、「同じ $\Delta f$ でも応答の大きさが変わる」様子を確認

• df_max

  • 10 MHz → 5 MHz → 1 MHz と小さくして、波形が読み取りづらくなる様子を見る

• tau_ms

  • 1 → 10 → 30 ms に変えて、パルス幅の変化と実験データの時間応答の対応を考える

複数共振器側では：

• n_res と spacing_MHz を変えて
  「周波数軸にチャンネルを詰める」とどう見えるか試す
• hit_index を変え、どのトレースがどの共振器に対応しているか確認する

といった遊び方がおすすめです。

9. コード全文（参考）

コード全文を見る

#!/usr/bin/env python3
"""
Fixed-tone microwave readout simulator with:
  (1) Single resonator + complex S21 (I/Q) + vector plot
  (2) Multiple resonators (MW-Mux-like) demonstration

TES / MW-Mux を意識して、S21 は「ピーク型」ではなく
共振でノッチ (減衰) が出る「ノッチ型 (notch-type)」としてモデル化している。

このスクリプトの目的:
- 「S21 がノッチになる共振器」を数式とプロットでイメージする
- fixed-tone readout で I/Q がどう動くかを可視化する
- Q 値や共振周波数などの基礎パラメータを print で確認しながら学ぶ

実行すると、現在のディレクトリに PNG ファイルとして図が保存される:
- single_resonator_linear.png
- single_resonator_db.png
- multi_resonator_time_linear.png
- multi_resonator_static_linear.png
- multi_resonator_static_db.png
"""

import numpy as np
import matplotlib.pyplot as plt


# ============================================================
#  描画用のユーティリティ
# ============================================================

def apply_light_grid(ax):
    """
    学生さん向けに「薄めのグリッド」を統一的に適用する関数。
    grid(True) だけだと線が濃くて見づらいことがあるので、
    alpha, linewidth, linestyle を調整している。
    """
    ax.grid(True, alpha=0.3, linewidth=0.5, linestyle="--")


# ============================================================
#  共通ユーティリティ: 共振器の基礎パラメータを表示する関数
# ============================================================

def resonator_s21_complex(f, f_res, Q):
    """
    単純化された 1ポール共振器の複素伝達関数 S21(f) を返す (ノッチ型)。

    ここでは TES の MW-Mux でよく出てくる
    「通過帯 (|S21|≈1) の中に、共振周波数でノッチが落ちる」
    形を意識して、以下のような簡略モデルを採用する:

        S21(f) = 1 - 1 / (1 + j * x)

      ここで
        x     = (f - f_res) / gamma
        gamma = f_res / (2 Q)  (Half Width at Half Maximum: HWHM)

    ・f が f_res から十分離れているとき:
        x が大きくなり 1 / (1 + j x) ≈ 0 なので
        S21 ≈ 1  (ほぼ全通過)

    ・f = f_res の近傍:
        x = 0 なので 1 / (1 + j * 0) = 1
        よって S21 = 1 - 1 = 0 となり、ノッチ (深い落ち込み) になる。

    実際の MW-Mux では結合の強さや内部損失などでもう少し複雑になるが、
    ここでは「ノッチ型である」という直感をつかむことを目的にしている。
    """
    # ★ f はスカラーでも配列でもよいように計算している
    gamma = f_res / (2.0 * Q)  # HWHM [Hz]
    x = (f - f_res) / gamma
    S21 = 1.0 - 1.0 / (1.0 + 1j * x)
    return S21


def print_resonator_basic_params(f_res, Q, label="Resonator"):
    """
    共振器の基礎パラメータを print で表示する。

    引数:
        f_res : 共振周波数 [Hz]
        Q     : Q値 (loaded Q を想定)
        label : どの共振器か識別するラベル（例: "Res 0"）
    """
    # Half Width at Half Maximum (HWHM)
    gamma = f_res / (2.0 * Q)     # [Hz]
    # Full Width at Half Maximum (FWHM) ~ f_res / Q
    fwhm = f_res / Q              # [Hz]

    # ノッチ共振器の伝達特性を使って、
    # - 共振から十分離れたところ (f_res + 10*gamma)
    # - 共振ちょうど (f_res)
    # での透過率 |S21| を計算してみる
    S21_far = resonator_s21_complex(f_res + 10.0 * gamma, f_res, Q)
    S21_at_res = resonator_s21_complex(f_res, f_res, Q)

    amp_far = np.abs(S21_far)
    amp_res = np.abs(S21_at_res)

    # log10(0) にならないように微小量 1e-15 を足している
    amp_far_db = 20.0 * np.log10(amp_far + 1e-15)
    amp_res_db = 20.0 * np.log10(amp_res + 1e-15)

    print("========================================================")
    print(f"[{label}] 基本パラメータ")
    print(f"  共振周波数 f0 = {f_res:.3e} Hz  ({f_res*1e-9:.3f} GHz)")
    print(f"  Q 値 Q      = {Q:.3e}")
    print(f"  HWHM gamma  = {gamma:.3e} Hz")
    print(f"  FWHM ≈ f0/Q = {fwhm:.3e} Hz  ({fwhm*1e-6:.3f} MHz)")
    print("  共振から十分離れた周波数での透過率 (ほぼ全通過):")
    print(f"    |S21| ≈ {amp_far:.3f}   ({amp_far_db:.2f} dB)")
    print("  共振周波数ちょうどでの透過率 (ノッチの深さ):")
    print(f"    |S21| ≈ {amp_res:.3e}   ({amp_res_db:.2f} dB)")
    print("    -> この簡略モデルでは理想的にノッチが 0 まで落ちる (実機ではここまで深くない)。")
    print("========================================================\n")


# ============================================================
#  (1) 単一共振器: fixed-tone + I/Q
# ============================================================

def simulate_single_resonator(
    f0=5.0e9,
    Q=300.0,           # ★ デフォルト Q を 300 に変更
    df_max=10.0e6,
    tau_ms=10.0,
    t_total_ms=50.0,
    dt_us=20.0,
    t_event_ms=5.0,
    probe_offset_MHz=0.0,
):
    """
    単一共振器の fixed-tone readout を時間領域でシミュレートする。

    - 共振周波数 f_res(t) が、ある時刻 t_event_ms 以降で
      指数関数的にシフトして元に戻る簡単なモデル。
    - 実際には TES への X 線入射に対応しているとイメージすると良い。

    戻り値:
        t_ms      : 時間 [ms]
        f_res_Hz  : 各時刻での共振周波数 [Hz]
        S21_t     : 各時刻で固定トーンに対する S21(t)
        f_probe_Hz: 固定トーンの周波数 [Hz]
    """
    # --- 時間軸の設定 ---
    dt_s = dt_us * 1.0e-6            # サンプリング間隔 [s]
    t_total_s = t_total_ms * 1.0e-3  # 総観測時間 [s]
    t_event_s = t_event_ms * 1.0e-3  # 事象発生時刻 [s]
    tau_s = tau_ms * 1.0e-3          # 緩和時定数 [s]

    # 0 から t_total_s まで dt_s 刻みの時間配列
    # ★ arange は「終点を含まない」ので、少し余裕を持たせている
    t_s = np.arange(0.0, t_total_s, dt_s)
    t_ms = t_s * 1.0e3              # プロットしやすいよう [ms] も作る

    # --- 共振周波数 f_res(t) の時間変化 ---
    f_res_Hz = np.empty_like(t_s)
    for i, t in enumerate(t_s):
        if t < t_event_s:
            # 事象発生前は共振周波数 f0 のまま
            f_res_Hz[i] = f0
        else:
            # 事象発生後は df_max だけシフトして、指数関数的に元に戻る
            f_res_Hz[i] = f0 + df_max * np.exp(-(t - t_event_s) / tau_s)

    # --- fixed-tone の周波数 ---
    # 実際の readout では、共振周波数の近くに一定周波数のトーンを打ちっぱなしにする。
    f_probe_Hz = f0 + probe_offset_MHz * 1.0e6

    # 各時刻での S21(t)（固定トーン f_probe_Hz に対する伝達関数）
    # ★ 引数の順番に注意: (f_probe, f_res(t), Q)
    S21_t = resonator_s21_complex(f_probe_Hz, f_res_Hz, Q)

    return t_ms, f_res_Hz, S21_t, f_probe_Hz


def plot_single_resonator(t_ms, f_res_Hz, S21_t, f_probe_Hz, f0, Q):
    """
    単一共振器の結果をプロットする。

    Figure 1 (リニア表示):
      1) Δf_res(t) [MHz]
      2) |S21(t)| vs Time (linear)
      3) I-Q 平面での軌跡
      4) 静的な共振曲線 |S21(f)| (notch, linear)

    Figure 2 (dB 表示):
      1) S21(t) [dB] vs Time
      2) 静的な共振曲線 S21(f) [dB]

    ★ 各 Figure は PNG ファイルとしても保存する。
    """
    # 共振周波数の変化量 (f_res - f0) を [MHz] 単位で計算
    df_MHz = (f_res_Hz - f0) * 1.0e-6

    # S21 の振幅（透過率）と dB 表示
    amp_t = np.abs(S21_t)                      # |S21| (linear, 透過率)
    amp_db_t = 20.0 * np.log10(amp_t + 1e-15)  # [dB], log(0) 回避のため微小量を足す

    # I/Q 成分を取り出す (I=実部, Q=虚部)
    I_t = np.real(S21_t)
    Q_t = np.imag(S21_t)

    # -------------------------------
    # Figure 1: リニア表示まとめ
    # -------------------------------
    fig1 = plt.figure(figsize=(10, 8))

    # 1) Δf_res(t)
    ax1 = fig1.add_subplot(2, 2, 1)
    ax1.plot(t_ms, df_MHz)
    ax1.set_xlabel("Time [ms]")
    ax1.set_ylabel("Δf_res [MHz]")
    ax1.set_title("Resonance frequency shift vs time")
    apply_light_grid(ax1)

    # 2) |S21(t)| (linear)
    ax2 = fig1.add_subplot(2, 2, 2)
    ax2.plot(t_ms, amp_t)
    ax2.set_xlabel("Time [ms]")
    ax2.set_ylabel("Transmission |S21| (linear)")
    ax2.set_title("Fixed-tone microwave readout (amplitude, linear)")
    apply_light_grid(ax2)

    # 3) I-Q trajectory
    ax3 = fig1.add_subplot(2, 2, 3)
    ax3.plot(I_t, Q_t, "-")
    ax3.set_xlabel("I = Re(S21)")
    ax3.set_ylabel("Q = Im(S21)")
    ax3.set_title("I-Q trajectory at fixed tone")
    apply_light_grid(ax3)
    ax3.set_aspect("equal", "box")

    # 開始点・終了点・最大シフト付近をマーカーで示す
    ax3.plot(I_t[0], Q_t[0], "go", label="start")
    ax3.plot(I_t[-1], Q_t[-1], "ro", label="end")
    idx_event = np.argmax(df_MHz)  # 周波数シフト最大付近
    ax3.plot(I_t[idx_event], Q_t[idx_event], "bo", label="around max shift")
    ax3.legend(loc="best")

    # 4) 静的な共振曲線 |S21(f)| (linear)
    ax4 = fig1.add_subplot(2, 2, 4)
    span_MHz = 30.0
    f_MHz = np.linspace(-span_MHz, span_MHz, 1000)
    f_Hz = f0 + f_MHz * 1.0e6
    S21_static = resonator_s21_complex(f_Hz, f0, Q)
    amp_static = np.abs(S21_static)

    ax4.plot(f_MHz, amp_static, label="|S21(f)| (linear)")
    ax4.axvline((f_probe_Hz - f0) * 1.0e-6, linestyle="--", label="probe tone")
    ax4.set_xlabel("f - f0 [MHz]")
    ax4.set_ylabel("Transmission |S21| (linear)")
    ax4.set_title("Static resonance (notch) and probe frequency (linear)")
    apply_light_grid(ax4)
    ax4.legend(loc="best", fontsize=8)

    plt.tight_layout()

    # ★ 図を保存（カレントディレクトリに PNG として）
    fig1.savefig("single_resonator_linear.png", dpi=150, bbox_inches="tight")
    print("Saved: single_resonator_linear.png")

    plt.show()

    # -------------------------------
    # Figure 2: dB 表示まとめ
    # -------------------------------
    fig2, (ax5, ax6) = plt.subplots(1, 2, figsize=(10, 4))

    # 5) S21(t) [dB] vs time
    ax5.plot(t_ms, amp_db_t)
    ax5.set_xlabel("Time [ms]")
    ax5.set_ylabel("S21 [dB]")
    ax5.set_title("Fixed-tone microwave readout (amplitude, dB)")
    apply_light_grid(ax5)

    # 6) 静的な共振曲線 S21(f) [dB]
    amp_static_db = 20.0 * np.log10(amp_static + 1e-15)
    ax6.plot(f_MHz, amp_static_db)
    ax6.axvline((f_probe_Hz - f0) * 1.0e-6, linestyle="--", label="probe tone")
    ax6.set_xlabel("f - f0 [MHz]")
    ax6.set_ylabel("S21 [dB]")
    ax6.set_title("Static resonance (notch) in dB")
    apply_light_grid(ax6)
    ax6.legend(loc="best", fontsize=8)

    plt.tight_layout()

    # ★ 図を保存
    fig2.savefig("single_resonator_db.png", dpi=150, bbox_inches="tight")
    print("Saved: single_resonator_db.png")

    plt.show()


# ============================================================
#  (2) 複数共振器 (MW-Mux 的なデモ)
# ============================================================

def simulate_multi_resonators(
    n_res=4,
    f0_center=5.0e9,
    spacing_MHz=50.0,
    Q=300.0,           # ★ デフォルト Q を 300 に変更
    df_max=10.0e6,
    tau_ms=10.0,
    t_total_ms=50.0,
    dt_us=20.0,
    t_event_ms=5.0,
    hit_index=1,
):
    """
    n_res 本の共振器を周波数方向に等間隔に並べ、
    そのうち 1 本だけが X 線ヒットで周波数シフトを経験するモデル。

    - f0_center を中心として、spacing_MHz [MHz] 間隔で共振器を並べる。
    - 各共振器ごとにプローブトーンは「その共振器の共振周波数」に固定。
    - hit_index 番目の共振器だけが df_max だけ周波数シフトする。

    戻り値:
        t_ms     : 時間 [ms]
        f0_list  : 各共振器の中心共振周波数 [Hz] の配列 (長さ n_res)
        S21_multi: shape = (n_res, n_t) の複素配列
                   S21_multi[k, i] が「k 番目の共振器の i 番目の時刻での S21」
    """
    # 共振器のインデックス (0, 1, 2, ..., n_res-1) を
    # 中心からの相対オフセット (-..., 0, +...) に変換
    idx_offset = np.arange(n_res) - (n_res - 1) / 2.0
    f0_list = f0_center + idx_offset * spacing_MHz * 1.0e6  # 各共振器の f0

    # --- 時間軸 ---
    dt_s = dt_us * 1.0e-6
    t_total_s = t_total_ms * 1.0e-3
    t_event_s = t_event_ms * 1.0e-3
    tau_s = tau_ms * 1.0e-3

    t_s = np.arange(0.0, t_total_s, dt_s)
    t_ms = t_s * 1.0e3
    n_t = len(t_s)

    # S21_multi[k, i] : k 番目の共振器の i 番時刻での S21
    S21_multi = np.zeros((n_res, n_t), dtype=complex)

    # 各共振器ごとに、時間変化する f_res_k(t) と S21_k(t) を計算
    for k in range(n_res):
        f0_k = f0_list[k]
        f_probe_k = f0_k  # ここでは probe トーンを共振周波数に固定（オフセット 0）

        f_res_k = np.empty_like(t_s)
        for i, t in enumerate(t_s):
            if k == hit_index and t >= t_event_s:
                # この共振器だけ事象によって df_max シフトする
                f_res_k[i] = f0_k + df_max * np.exp(-(t - t_event_s) / tau_s)
            else:
                # 他の共振器、または事象前は f0_k のまま
                f_res_k[i] = f0_k

        # ★ ここも引数の順序: (f_probe_k, f_res_k(t), Q)
        S21_k = resonator_s21_complex(f_probe_k, f_res_k, Q)
        S21_multi[k, :] = S21_k

    return t_ms, f0_list, S21_multi


def plot_multi_resonators(t_ms, f0_center, f0_list, S21_multi, Q, spacing_MHz=50.0):
    """
    複数共振器の |S21_k(t)| と静的共振曲線をプロット。

    ・Figure 1: time-domain |S21_k(t)| (linear)
    ・Figure 2: static |S21(f)| (linear)
    ・Figure 3: static S21(f) [dB]

    ★ 各 Figure は PNG ファイルとしても保存する。
    """
    n_res, n_t = S21_multi.shape

    # -------------------------------
    # Figure 1: time-domain amplitude (linear)
    # -------------------------------
    fig1, ax1 = plt.subplots(figsize=(7, 4))
    for k in range(n_res):
        amp_k = np.abs(S21_multi[k, :])
        label = f"Resonator {k} (f0 = {f0_list[k]*1e-9:.3f} GHz)"
        ax1.plot(t_ms, amp_k, label=label)
    ax1.set_xlabel("Time [ms]")
    ax1.set_ylabel("|S21| at each probe tone (linear)")
    ax1.set_title("MW-Mux-like: amplitude at each tone vs time (notch-type, linear)")
    apply_light_grid(ax1)
    ax1.legend(loc="best", fontsize=8)

    plt.tight_layout()
    fig1.savefig("multi_resonator_time_linear.png", dpi=150, bbox_inches="tight")
    print("Saved: multi_resonator_time_linear.png")
    plt.show()

    # -------------------------------
    # Figure 2: static resonance curves (linear)
    # -------------------------------
    fig2, ax2 = plt.subplots(figsize=(7, 4))

    span_total_MHz = spacing_MHz * (n_res + 1)
    f_MHz = np.linspace(-span_total_MHz, span_total_MHz, 2000)
    f_Hz = f0_center + f_MHz * 1.0e6

    for k in range(n_res):
        f0_k = f0_list[k]
        S_static_k = resonator_s21_complex(f_Hz, f0_k, Q)
        amp_k = np.abs(S_static_k)
        label = f"Res {k} (f0={f0_k*1e-9:.3f} GHz)"
        ax2.plot(f_MHz, amp_k, label=label)

        # 各共振器の位置に縦線を引いて分かりやすくする
        offset_k_MHz = (f0_k - f0_center) * 1.0e-6
        ax2.axvline(offset_k_MHz, linestyle="--", alpha=0.3)

    ax2.set_xlabel("f - f_center [MHz]")
    ax2.set_ylabel("Transmission |S21| (linear)")
    ax2.set_title("Static resonance curves (notch-type, linear)")
    apply_light_grid(ax2)
    ax2.legend(loc="best", fontsize=8)

    plt.tight_layout()
    fig2.savefig("multi_resonator_static_linear.png", dpi=150, bbox_inches="tight")
    print("Saved: multi_resonator_static_linear.png")
    plt.show()

    # -------------------------------
    # Figure 3: static resonance curves (dB)
    # -------------------------------
    fig3, ax3 = plt.subplots(figsize=(7, 4))

    for k in range(n_res):
        f0_k = f0_list[k]
        S_static_k = resonator_s21_complex(f_Hz, f0_k, Q)
        amp_k = np.abs(S_static_k)
        amp_k_db = 20.0 * np.log10(amp_k + 1e-15)
        label = f"Res {k} (f0={f0_k*1e-9:.3f} GHz)"
        ax3.plot(f_MHz, amp_k_db, label=label)

        offset_k_MHz = (f0_k - f0_center) * 1.0e-6
        ax3.axvline(offset_k_MHz, linestyle="--", alpha=0.3)

    ax3.set_xlabel("f - f_center [MHz]")
    ax3.set_ylabel("S21 [dB]")
    ax3.set_title("Static resonance curves (notch-type, dB)")
    apply_light_grid(ax3)
    ax3.legend(loc="best", fontsize=8)

    plt.tight_layout()
    fig3.savefig("multi_resonator_static_db.png", dpi=150, bbox_inches="tight")
    print("Saved: multi_resonator_static_db.png")
    plt.show()


def plot_single_resonator_time_series(t_ms, S21_t, f_res_Hz=None, f0=None):
    """
    時間を横軸にして、
      - I(t) = Re(S21(t))
      - Q(t) = Im(S21(t))
      - |S21(t)| (linear)
    を 3 段のサブプロットで表示する。

    第3パネルでは、
      - 左軸 : |S21(t)| (linear)
      - 右軸 : Δf_res(t) [MHz]（オプション）
    とし、色と凡例で混乱しないようにする。
    """
    # --- I, Q, 振幅 ---
    I_t = np.real(S21_t)
    Q_t = np.imag(S21_t)
    amp_t = np.abs(S21_t)

    fig, axes = plt.subplots(3, 1, figsize=(8, 8), sharex=True)

    # ===============================
    # 1) I(t)
    # ===============================
    ax_I = axes[0]
    ax_I.plot(t_ms, I_t, color="tab:blue")
    ax_I.set_ylabel("I(t) = Re(S21)")
    ax_I.set_title("Time-domain view: I(t), Q(t), |S21(t)|")
    apply_light_grid(ax_I)

    # ===============================
    # 2) Q(t)
    # ===============================
    ax_Q = axes[1]
    ax_Q.plot(t_ms, Q_t, color="tab:green")
    ax_Q.set_ylabel("Q(t) = Im(S21)")
    apply_light_grid(ax_Q)

    # ===============================
    # 3) |S21(t)| + Δf(t)
    # ===============================
    ax_A = axes[2]

    # 左軸: |S21|
    line_amp, = ax_A.plot(
        t_ms, amp_t,
        color="tab:blue",
        label="|S21(t)| (linear)"
    )
    ax_A.set_ylabel("|S21(t)| (linear)", color="tab:blue")
    ax_A.tick_params(axis="y", labelcolor="tab:blue")
    ax_A.set_xlabel("Time [ms]")
    apply_light_grid(ax_A)

    lines = [line_amp]
    labels = [line_amp.get_label()]

    # 右軸: Δf_res(t)
    if f_res_Hz is not None and f0 is not None:
        df_MHz = (f_res_Hz - f0) * 1.0e-6

        ax_df = ax_A.twinx()
        line_df, = ax_df.plot(
            t_ms, df_MHz,
            linestyle="--",
            color="tab:orange",
            label="Δf_res(t) [MHz]"
        )
        ax_df.set_ylabel("Δf_res(t) [MHz]", color="tab:orange")
        ax_df.tick_params(axis="y", labelcolor="tab:orange")

        # 凡例用に追加
        lines.append(line_df)
        labels.append(line_df.get_label())

    # 凡例（左右軸をまたぐので手動指定）
    ax_A.legend(lines, labels, loc="best", fontsize=9)

    plt.tight_layout()

    # 保存
    fig.savefig(
        "single_resonator_time_I_Q_amp_df.png",
        dpi=150,
        bbox_inches="tight"
    )
    print("Saved: single_resonator_time_I_Q_amp_df.png")

    plt.show()


# ============================================================
#  main: デモ用ラッパ関数
# ============================================================

def demo_single_resonator():
    """単一共振器 + I/Q デモ (ノッチ型 S21)"""
    # ここで設定しているのは「5 GHz, Q=300 の共振器」をイメージ
    f0 = 5.0e9       # 共振周波数 [Hz]
    Q = 300.0        # ★ Q 値 (比較的低めにして共振を太く見せる)
    df_max = 10.0e6  # 事象による最大周波数シフト [Hz]
    tau_ms = 10.0    # 緩和時定数 [ms]
    t_total_ms = 50.0
    dt_us = 20.0
    t_event_ms = 5.0
    probe_offset_MHz = 0.0  # probe トーンを f0 ど真ん中に置く

    # --- 共振器の基本パラメータを print しておく ---
    print_resonator_basic_params(f0, Q, label="Single resonator")

    # シミュレーション実行
    t_ms, f_res_Hz, S21_t, f_probe_Hz = simulate_single_resonator(
        f0=f0,
        Q=Q,
        df_max=df_max,
        tau_ms=tau_ms,
        t_total_ms=t_total_ms,
        dt_us=dt_us,
        t_event_ms=t_event_ms,
        probe_offset_MHz=probe_offset_MHz,
    )

    # プロット
    plot_single_resonator(t_ms, f_res_Hz, S21_t, f_probe_Hz, f0, Q)

    # ★ 追加: 時間 vs I, Q, |S21| の図
    plot_single_resonator_time_series(t_ms, S21_t, f_res_Hz=f_res_Hz, f0=f0)


def demo_multi_resonators():
    """複数共振器 (MW-Mux 的, ノッチ型 S21) デモ"""
    n_res = 4
    f0_center = 5.0e9
    spacing_MHz = 50.0
    Q = 300.0       # ★ こちらも Q=300
    df_max = 10.0e6
    tau_ms = 10.0
    t_total_ms = 50.0
    dt_us = 20.0
    t_event_ms = 5.0
    hit_index = 1  # このインデックスの共振器だけ周波数シフトさせる

    # シミュレーション実行
    t_ms, f0_list, S21_multi = simulate_multi_resonators(
        n_res=n_res,
        f0_center=f0_center,
        spacing_MHz=spacing_MHz,
        Q=Q,
        df_max=df_max,
        tau_ms=tau_ms,
        t_total_ms=t_total_ms,
        dt_us=dt_us,
        t_event_ms=t_event_ms,
        hit_index=hit_index,
    )

    # --- 各共振器の基本パラメータを print ---
    print("=== Multi-resonator configuration ===")
    for k, f0_k in enumerate(f0_list):
        print_resonator_basic_params(f0_k, Q, label=f"Resonator {k}")

    # プロット
    plot_multi_resonators(
        t_ms=t_ms,
        f0_center=f0_center,
        f0_list=f0_list,
        S21_multi=S21_multi,
        Q=Q,
        spacing_MHz=spacing_MHz,
    )


def main():
    # ★ 学生向けには main() の中でどんな処理が呼ばれているかも重要
    #   - 単一共振器のデモ
    #   - 複数共振器 (MW-Mux 的) のデモ
    demo_single_resonator()
    demo_multi_resonators()


if __name__ == "__main__":
    main()

関連記事