キルヒホッフの法則 — 電流則・電圧則による回路解析

     R1         R2
v1 --/\/\/-- v2 --/\/\/-- v3
|            |            |
Is1↑        R3           R4
|            |            |
GND ---------+------------ GND

    R1        R2
+--/\/\/--+--/\/\/--+
|         |         |
Vs ↑     R3        R4
|         |         |
+---------+---------+

import numpy as np
import matplotlib.pyplot as plt

# --- 3節点回路の節点解析 ---
# 回路構成:
#   節点1 -- R1=2Ω -- 節点2 -- R2=4Ω -- 節点3
#   節点1に電流源 Is=5A が流入
#   節点2 -- R3=6Ω -- GND
#   節点3 -- R4=3Ω -- GND

# 抵抗値
R1 = 2.0   # Ω
R2 = 4.0   # Ω
R3 = 6.0   # Ω
R4 = 3.0   # Ω

# コンダクタンス
G1 = 1 / R1
G2 = 1 / R2
G3 = 1 / R3
G4 = 1 / R4

# コンダクタンス行列の構成
# 対角要素: 節点iに接続された全コンダクタンスの和
# 非対角要素: -（節点iとjを結ぶコンダクタンス）
G = np.array([
    [G1,       -G1,        0      ],
    [-G1,      G1+G2+G3,  -G2     ],
    [0,        -G2,        G2+G4  ]
])

# 電流源ベクトル（節点への流入を正）
Is = np.array([5.0, 0.0, 0.0])

print("=== コンダクタンス行列 G ===")
print(G)
print(f"\nG は対称行列: {np.allclose(G, G.T)}")
print(f"G の条件数: {np.linalg.cond(G):.2f}")

# 連立方程式を解く: Gv = Is
v = np.linalg.solve(G, Is)

print(f"\n=== 節点電圧 ===")
for i, vi in enumerate(v):
    print(f"  v{i+1} = {vi:.4f} V")

# 各枝の電流を計算
I_R1 = (v[0] - v[1]) / R1
I_R2 = (v[1] - v[2]) / R2
I_R3 = v[1] / R3
I_R4 = v[2] / R4

print(f"\n=== 枝電流 ===")
print(f"  I_R1 (節点1→2): {I_R1:.4f} A")
print(f"  I_R2 (節点2→3): {I_R2:.4f} A")
print(f"  I_R3 (節点2→GND): {I_R3:.4f} A")
print(f"  I_R4 (節点3→GND): {I_R4:.4f} A")

# KCL検証
print(f"\n=== KCL検証 ===")
kcl_node1 = Is[0] - I_R1
kcl_node2 = I_R1 - I_R2 - I_R3
kcl_node3 = I_R2 - I_R4
print(f"  節点1: Is - I_R1 = {kcl_node1:.6f} A (= 0?)")
print(f"  節点2: I_R1 - I_R2 - I_R3 = {kcl_node2:.6f} A (= 0?)")
print(f"  節点3: I_R2 - I_R4 = {kcl_node3:.6f} A (= 0?)")

# 消費電力
P_total = Is[0] * v[0]
P_R1 = I_R1**2 * R1
P_R2 = I_R2**2 * R2
P_R3 = I_R3**2 * R3
P_R4 = I_R4**2 * R4

print(f"\n=== 電力収支 ===")
print(f"  電流源供給電力: {P_total:.4f} W")
print(f"  R1消費電力: {P_R1:.4f} W")
print(f"  R2消費電力: {P_R2:.4f} W")
print(f"  R3消費電力: {P_R3:.4f} W")
print(f"  R4消費電力: {P_R4:.4f} W")
print(f"  全抵抗消費電力合計: {P_R1+P_R2+P_R3+P_R4:.4f} W")
print(f"  電力収支誤差: {abs(P_total - (P_R1+P_R2+P_R3+P_R4)):.6f} W")

import numpy as np
import matplotlib.pyplot as plt

class CircuitSolver:
    """節点電圧法による回路解析ソルバー"""

    def __init__(self, n_nodes):
        """
        n_nodes: 節点数（GNDを除く）
        """
        self.n = n_nodes
        self.G = np.zeros((n_nodes, n_nodes))  # コンダクタンス行列
        self.Is = np.zeros(n_nodes)             # 電流源ベクトル

    def add_resistor(self, node_i, node_j, R):
        """
        抵抗を追加
        node_i, node_j: 接続する節点（0=GND, 1〜n）
        R: 抵抗値 [Ω]
        """
        G = 1.0 / R
        # 0はGNDノードなので、1-indexedを0-indexedに変換
        i = node_i - 1
        j = node_j - 1

        if node_i > 0 and node_j > 0:
            self.G[i, i] += G
            self.G[j, j] += G
            self.G[i, j] -= G
            self.G[j, i] -= G
        elif node_i > 0:  # node_j == 0 (GND)
            self.G[i, i] += G
        elif node_j > 0:  # node_i == 0 (GND)
            self.G[j, j] += G

    def add_current_source(self, node_into, current):
        """
        電流源を追加（node_into に流入する電流）
        node_into: 電流が流入する節点（1〜n）
        current: 電流値 [A]
        """
        if node_into > 0:
            self.Is[node_into - 1] += current

    def solve(self):
        """連立方程式を解いて節点電圧を返す"""
        self.v = np.linalg.solve(self.G, self.Is)
        return self.v

    def get_branch_current(self, node_i, node_j, R):
        """枝電流を計算（node_i → node_j）"""
        vi = self.v[node_i - 1] if node_i > 0 else 0.0
        vj = self.v[node_j - 1] if node_j > 0 else 0.0
        return (vi - vj) / R

# テスト: ホイートストンブリッジ回路
#
#          R1=100Ω
#   1 ----/\/\/---- 2
#   |               |
#   |    R5=50Ω     |
#   |  3---/\/\/    |
#   |  |        |   |
#   | R3=200Ω  R4=300Ω
#   |  |        |   |
#   +--+--------+---+
#       GND
#   Vs=10V → 節点1に10Vの電圧源
#   （電圧源を扱うために、節点1の電圧を10Vに固定するアプローチ）

# 簡略化: 電圧源の代わりに、節点1の電圧をVs=10Vと既知とし、
# 残りの節点について解く
# 節点: 2, 3（節点1は電圧既知）

# ホイートストンブリッジ
R1 = 100.0   # 節点1→節点2
R2 = 150.0   # 節点1→節点3
R3 = 200.0   # 節点2→GND
R4 = 300.0   # 節点3→GND
R5 = 50.0    # 節点2→節点3（ブリッジ抵抗）
Vs = 10.0    # 電圧源

# 節点1の電圧は Vs = 10V（既知）
# 節点2、節点3について KCL を立てる
# 節点2: (Vs-v2)/R1 + (v3-v2)/R5 - v2/R3 = 0
# 節点3: (Vs-v3)/R2 + (v2-v3)/R5 - v3/R4 = 0

G1 = 1/R1; G2 = 1/R2; G3 = 1/R3; G4 = 1/R4; G5 = 1/R5

# 行列形式: Gv = b
G_matrix = np.array([
    [G1 + G3 + G5,  -G5],
    [-G5,            G2 + G4 + G5]
])

b = np.array([G1 * Vs, G2 * Vs])

v = np.linalg.solve(G_matrix, b)
v2, v3 = v[0], v[1]

print("=== ホイートストンブリッジ回路 ===")
print(f"Vs = {Vs} V")
print(f"R1={R1}Ω, R2={R2}Ω, R3={R3}Ω, R4={R4}Ω, R5={R5}Ω")
print(f"\n節点電圧:")
print(f"  v1 = {Vs:.4f} V (電圧源)")
print(f"  v2 = {v2:.4f} V")
print(f"  v3 = {v3:.4f} V")

# ブリッジ電圧（検出電圧）
V_bridge = v2 - v3
print(f"\nブリッジ電圧 (v2 - v3) = {V_bridge:.4f} V")

# 平衡条件の確認: R1*R4 == R2*R3 のとき V_bridge = 0
print(f"\n平衡条件: R1*R4 = {R1*R4:.0f}, R2*R3 = {R2*R3:.0f}")
if abs(R1*R4 - R2*R3) < 1e-10:
    print("→ ブリッジは平衡状態")
else:
    print("→ ブリッジは非平衡状態")

# R3を変化させたときのブリッジ電圧
R3_range = np.linspace(50, 500, 100)
V_bridge_range = []

for R3_var in R3_range:
    G3_var = 1/R3_var
    G_var = np.array([
        [G1 + G3_var + G5,  -G5],
        [-G5,                G2 + G4 + G5]
    ])
    b_var = np.array([G1 * Vs, G2 * Vs])
    v_var = np.linalg.solve(G_var, b_var)
    V_bridge_range.append(v_var[0] - v_var[1])

V_bridge_range = np.array(V_bridge_range)

# 平衡点
R3_balance = R1 * R4 / R2
print(f"\n平衡条件の R3 = R1*R4/R2 = {R3_balance:.1f} Ω")

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(R3_range, V_bridge_range, 'b-', linewidth=2)
ax.axhline(y=0, color='gray', linestyle='-', linewidth=0.5)
ax.axvline(x=R3_balance, color='r', linestyle='--', alpha=0.7,
           label=f'$R_3$ = {R3_balance:.1f} Ω (balanced)')
ax.set_xlabel('$R_3$ [Ω]')
ax.set_ylabel('Bridge Voltage $V_{bridge}$ [V]')
ax.set_title('Wheatstone Bridge: Bridge Voltage vs $R_3$')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('wheatstone_bridge.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# --- 2メッシュ回路の網目電流法 ---
#
#     R1=10Ω      R2=20Ω
# +---/\/\/---+---/\/\/---+
# |           |           |
# Vs=12V↑    R3=30Ω      R4=40Ω
# |           |           |
# +-----------+-----------+
#
# メッシュ1（左ループ、時計回り）: i1
# メッシュ2（右ループ、時計回り）: i2

R1 = 10.0
R2 = 20.0
R3 = 30.0
R4 = 40.0
Vs = 12.0

# 抵抗行列の構成
# メッシュ1のKVL: Vs - R1*i1 - R3*(i1-i2) = 0
#   → (R1+R3)*i1 - R3*i2 = Vs
# メッシュ2のKVL: -R2*i2 - R4*i2 - R3*(i2-i1) = 0
#   → -R3*i1 + (R2+R3+R4)*i2 = 0

R_matrix = np.array([
    [R1 + R3,    -R3],
    [-R3,         R2 + R3 + R4]
])

V_source = np.array([Vs, 0])

# 網目電流を求解
i_mesh = np.linalg.solve(R_matrix, V_source)

print("=== 2メッシュ回路の解析 ===")
print(f"R1={R1}Ω, R2={R2}Ω, R3={R3}Ω, R4={R4}Ω, Vs={Vs}V")
print(f"\n抵抗行列 R:")
print(R_matrix)
print(f"\nR は対称行列: {np.allclose(R_matrix, R_matrix.T)}")

print(f"\n=== 網目電流 ===")
print(f"  i1 = {i_mesh[0]:.4f} A")
print(f"  i2 = {i_mesh[1]:.4f} A")

# 各枝の電流
I_R1 = i_mesh[0]
I_R2 = i_mesh[1]
I_R3 = i_mesh[0] - i_mesh[1]
I_R4 = i_mesh[1]

print(f"\n=== 枝電流 ===")
print(f"  I_R1 = {I_R1:.4f} A")
print(f"  I_R2 = {I_R2:.4f} A")
print(f"  I_R3 = {I_R3:.4f} A (i1 - i2)")
print(f"  I_R4 = {I_R4:.4f} A")

# 各素子の電圧降下
V_R1 = I_R1 * R1
V_R2 = I_R2 * R2
V_R3 = I_R3 * R3
V_R4 = I_R4 * R4

print(f"\n=== 電圧降下 ===")
print(f"  V_R1 = {V_R1:.4f} V")
print(f"  V_R2 = {V_R2:.4f} V")
print(f"  V_R3 = {V_R3:.4f} V")
print(f"  V_R4 = {V_R4:.4f} V")

# KVL検証
kvl_mesh1 = Vs - V_R1 - V_R3
kvl_mesh2 = -V_R2 - V_R4 + V_R3
print(f"\n=== KVL検証 ===")
print(f"  メッシュ1: Vs - V_R1 - V_R3 = {kvl_mesh1:.6f} V (= 0?)")
print(f"  メッシュ2: V_R3 - V_R2 - V_R4 = {kvl_mesh2:.6f} V (= 0?)")

# --- 節点電圧法で同じ回路を解いて検証 ---
# 節点: A（R1とR2の接続点）, B（Vsの+端子）
# GND = Vsの-端子 = R3,R4の下端
# 節点Bの電圧 = Vs = 12V（既知）
# 節点Aについて KCL: (VB-VA)/R1 + (0-VA)/R3... ではなく
# 節点Aの KCL: (Vs-vA)/R1 - vA/R3 - vA/(R2+R4) ... ではなく

# 正しい節点設定:
# 節点1 = R1とR3とR2の接続点（上辺中央）
# GND = 下辺
# Vs の + 端子 = 節点2（左上）
# R2 + R4 の先 = 右下 = GND

# 再整理:
# 節点A（上辺中央）: R1を通じてVs+端子に接続、R2を通じてGNDに接続(R2+R4直列)
# いや、より正確には:
# 上辺: [Vs+]--R1--[A]--R2--[B]
# 下辺: [Vs-]=[GND]
# A--R3--GND, B--R4--GND

# 節点A: (Vs-vA)/R1 = vA/R3 + (vA-vB)/R2  (不適切)
# ⇒ 正しくは 2節点(A, B):

G1n = 1/R1; G2n = 1/R2; G3n = 1/R3; G4n = 1/R4
Gn = np.array([
    [G1n + G2n + G3n, -G2n],
    [-G2n,             G2n + G4n]
])
bn = np.array([G1n * Vs, 0])
vn = np.linalg.solve(Gn, bn)

print(f"\n=== 節点電圧法による検証 ===")
print(f"  vA = {vn[0]:.4f} V")
print(f"  vB = {vn[1]:.4f} V")

# 節点電圧法から計算した枝電流
I_R1_node = (Vs - vn[0]) / R1
I_R3_node = vn[0] / R3
I_R2_node = (vn[0] - vn[1]) / R2
I_R4_node = vn[1] / R4

print(f"\n枝電流（節点電圧法）:")
print(f"  I_R1 = {I_R1_node:.4f} A")
print(f"  I_R2 = {I_R2_node:.4f} A")
print(f"  I_R3 = {I_R3_node:.4f} A")
print(f"  I_R4 = {I_R4_node:.4f} A")

# 両手法の一致確認
print(f"\n=== 両手法の一致確認 ===")
print(f"  I_R1 差: {abs(I_R1 - I_R1_node):.2e} A")
print(f"  I_R3 差: {abs(I_R3 - I_R3_node):.2e} A")

import numpy as np
import matplotlib.pyplot as plt

# --- 重ね合わせの定理の検証 ---
#
# 回路: 2つの電源を持つ回路
#     R1=10Ω      R2=20Ω
# +---/\/\/---A---/\/\/---+
# |                       |
# Vs=12V↑      R3=15Ω   ↑Is=0.5A
# |              |        |
# +------+------+---------+
#        GND
#
# 節点A から R3 を通じてGNDに接続

R1 = 10.0
R2 = 20.0
R3 = 15.0
Vs = 12.0
Is = 0.5

# === 全電源同時作用 ===
# 節点A の KCL: (Vs-vA)/R1 + Is - vA/R3 - vA/R2 = 0
# ⇒ vA*(1/R1 + 1/R2 + 1/R3) = Vs/R1 + Is

G_total = 1/R1 + 1/R2 + 1/R3
b_total = Vs/R1 + Is
vA_total = b_total / G_total

print("=== 全電源同時作用 ===")
print(f"  vA = {vA_total:.4f} V")
print(f"  I_R1 = {(Vs - vA_total)/R1:.4f} A")
print(f"  I_R2 = {vA_total/R2:.4f} A")
print(f"  I_R3 = {vA_total/R3:.4f} A")

# === Vs のみ作用（Is を開放）===
b_vs_only = Vs/R1
vA_vs = b_vs_only / G_total

print(f"\n=== Vs のみ作用 (Is = 0) ===")
print(f"  vA(Vs) = {vA_vs:.4f} V")

# === Is のみ作用（Vs を短絡）===
# Vs を短絡 → R1 はGNDに直接接続
G_is_only = 1/R1 + 1/R2 + 1/R3  # 同じG行列（Vsの影響はbのみ）
b_is_only = Is
vA_is = b_is_only / G_is_only

print(f"\n=== Is のみ作用 (Vs = 0) ===")
print(f"  vA(Is) = {vA_is:.4f} V")

# === 重ね合わせの検証 ===
vA_superposition = vA_vs + vA_is
print(f"\n=== 重ね合わせの検証 ===")
print(f"  vA(Vs) + vA(Is) = {vA_vs:.4f} + {vA_is:.4f} = {vA_superposition:.4f} V")
print(f"  vA(全電源) = {vA_total:.4f} V")
print(f"  誤差: {abs(vA_total - vA_superposition):.2e} V")

# 電力が重ね合わせに従わないことの確認
P_R3_total = (vA_total)**2 / R3
P_R3_vs = (vA_vs)**2 / R3
P_R3_is = (vA_is)**2 / R3

print(f"\n=== 電力の重ね合わせ（成り立たない！）===")
print(f"  P_R3(全電源) = {P_R3_total:.4f} W")
print(f"  P_R3(Vs) + P_R3(Is) = {P_R3_vs:.4f} + {P_R3_is:.4f} = {P_R3_vs + P_R3_is:.4f} W")
print(f"  差異: {abs(P_R3_total - (P_R3_vs + P_R3_is)):.4f} W")
print(f"  → 電力は重ね合わせに従わない！")

# 可視化: 各電源の寄与
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 節点電圧の分解
labels = ['$V_s$ only', '$I_s$ only', 'Total']
voltages = [vA_vs, vA_is, vA_total]
colors_bar = ['tab:blue', 'tab:orange', 'tab:green']

bars = axes[0].bar(labels, voltages, color=colors_bar, alpha=0.8, edgecolor='black')
axes[0].set_ylabel('Node Voltage $v_A$ [V]')
axes[0].set_title('Superposition of Node Voltage')
axes[0].grid(True, alpha=0.3, axis='y')

# 数値ラベル
for bar, v in zip(bars, voltages):
    axes[0].text(bar.get_x() + bar.get_width()/2, v + 0.1, f'{v:.2f} V',
                ha='center', fontsize=11, fontweight='bold')

# 積み上げ棒グラフ（重ね合わせの図解）
axes[0].bar(['Superposition'], [vA_vs], color='tab:blue', alpha=0.5, label='$V_s$ contribution')
axes[0].bar(['Superposition'], [vA_is], bottom=[vA_vs], color='tab:orange', alpha=0.5,
            label='$I_s$ contribution')
axes[0].legend()

# 電力の比較（重ね合わせが成り立たないことを示す）
labels_p = ['$V_s$ only', '$I_s$ only', 'Sum', 'Actual']
powers = [P_R3_vs, P_R3_is, P_R3_vs + P_R3_is, P_R3_total]
colors_p = ['tab:blue', 'tab:orange', 'tab:red', 'tab:green']

bars_p = axes[1].bar(labels_p, powers, color=colors_p, alpha=0.8, edgecolor='black')
axes[1].set_ylabel('Power in $R_3$ [W]')
axes[1].set_title('Power Does NOT Obey Superposition')
axes[1].grid(True, alpha=0.3, axis='y')

for bar, p in zip(bars_p, powers):
    axes[1].text(bar.get_x() + bar.get_width()/2, p + 0.02, f'{p:.3f} W',
                ha='center', fontsize=10, fontweight='bold')

plt.tight_layout()
plt.savefig('superposition_verification.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# --- RLC直列回路のインピーダンス解析 ---
# KVL: Vs = V_R + V_L + V_C = I*(R + jωL + 1/(jωC))

R = 100.0     # Ω
L = 10e-3     # H (10 mH)
C = 1e-6      # F (1 μF)
Vs = 10.0     # V (振幅)

# 共振周波数
f0 = 1 / (2 * np.pi * np.sqrt(L * C))
omega0 = 2 * np.pi * f0
print(f"=== RLC直列回路 ===")
print(f"R = {R} Ω, L = {L*1e3} mH, C = {C*1e6} μF")
print(f"共振周波数 f0 = {f0:.1f} Hz")
print(f"共振角周波数 ω0 = {omega0:.1f} rad/s")

# Q値
Q = omega0 * L / R
print(f"Q値 = {Q:.2f}")

# 周波数応答
f_range = np.logspace(2, 5, 1000)  # 100 Hz 〜 100 kHz
omega_range = 2 * np.pi * f_range

# 複素インピーダンス（KVL）
Z = R + 1j * omega_range * L + 1 / (1j * omega_range * C)

# 電流フェーザ（KVLから）
I_phasor = Vs / Z

# 各素子の電圧フェーザ
V_R = I_phasor * R
V_L = I_phasor * 1j * omega_range * L
V_C = I_phasor / (1j * omega_range * C)

# KVL検証: Vs = V_R + V_L + V_C
kvl_check = np.max(np.abs(Vs - (V_R + V_L + V_C)))
print(f"\nKVL検証（最大誤差）: {kvl_check:.2e} V")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 電流の大きさ（周波数応答）
axes[0, 0].semilogx(f_range, np.abs(I_phasor) * 1e3, 'b-', linewidth=2)
axes[0, 0].axvline(x=f0, color='r', linestyle='--', alpha=0.7, label=f'$f_0$ = {f0:.0f} Hz')
axes[0, 0].set_xlabel('Frequency [Hz]')
axes[0, 0].set_ylabel('|I| [mA]')
axes[0, 0].set_title('Current Magnitude (Frequency Response)')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3, which='both')

# インピーダンスの大きさと位相
axes[0, 1].semilogx(f_range, np.abs(Z), 'g-', linewidth=2, label='|Z|')
axes[0, 1].set_xlabel('Frequency [Hz]')
axes[0, 1].set_ylabel('|Z| [Ω]')
axes[0, 1].set_title('Impedance Magnitude')
axes[0, 1].axvline(x=f0, color='r', linestyle='--', alpha=0.7)
axes[0, 1].axhline(y=R, color='gray', linestyle=':', alpha=0.7, label=f'R = {R} Ω')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3, which='both')

# 各素子の電圧振幅
axes[1, 0].semilogx(f_range, np.abs(V_R), 'b-', linewidth=1.5, label='$|V_R|$')
axes[1, 0].semilogx(f_range, np.abs(V_L), 'r-', linewidth=1.5, label='$|V_L|$')
axes[1, 0].semilogx(f_range, np.abs(V_C), 'g-', linewidth=1.5, label='$|V_C|$')
axes[1, 0].axhline(y=Vs, color='gray', linestyle=':', alpha=0.7, label=f'$V_s$ = {Vs} V')
axes[1, 0].axvline(x=f0, color='gray', linestyle='--', alpha=0.5)
axes[1, 0].set_xlabel('Frequency [Hz]')
axes[1, 0].set_ylabel('Voltage [V]')
axes[1, 0].set_title('Voltage Across Each Element')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3, which='both')
axes[1, 0].set_yscale('log')

# 共振時のフェーザ図
omega_res = omega0
I_res = Vs / (R + 1j * omega_res * L + 1 / (1j * omega_res * C))
VR_res = I_res * R
VL_res = I_res * 1j * omega_res * L
VC_res = I_res / (1j * omega_res * C)

ax_phasor = axes[1, 1]
origin = [0, 0]

# 電圧フェーザを描画
phasors = {
    '$V_R$': (VR_res, 'blue'),
    '$V_L$': (VL_res, 'red'),
    '$V_C$': (VC_res, 'green'),
    '$V_s$': (Vs + 0j, 'black')
}

for label, (V, color) in phasors.items():
    ax_phasor.annotate('', xy=(V.real, V.imag), xytext=(0, 0),
                        arrowprops=dict(arrowstyle='->', color=color, lw=2))
    ax_phasor.text(V.real * 1.1, V.imag * 1.1 + 0.3, label, color=color,
                   fontsize=11, fontweight='bold')

ax_phasor.set_xlim(-12, 12)
ax_phasor.set_ylim(-12, 12)
ax_phasor.set_xlabel('Real [V]')
ax_phasor.set_ylabel('Imaginary [V]')
ax_phasor.set_title(f'Phasor Diagram at Resonance (f = {f0:.0f} Hz)')
ax_phasor.set_aspect('equal')
ax_phasor.grid(True, alpha=0.3)
ax_phasor.axhline(y=0, color='gray', linewidth=0.5)
ax_phasor.axvline(x=0, color='gray', linewidth=0.5)

plt.tight_layout()
plt.savefig('rlc_circuit_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

# 共振時の数値
print(f"\n=== 共振時 (f = {f0:.0f} Hz) ===")
print(f"  |I| = {np.abs(I_res)*1e3:.2f} mA (= Vs/R = {Vs/R*1e3:.2f} mA)")
print(f"  |V_R| = {np.abs(VR_res):.4f} V (= Vs = {Vs} V)")
print(f"  |V_L| = {np.abs(VL_res):.4f} V (= Q*Vs = {Q*Vs:.4f} V)")
print(f"  |V_C| = {np.abs(VC_res):.4f} V (= Q*Vs = {Q*Vs:.4f} V)")
print(f"  V_L + V_C = {np.abs(VL_res + VC_res):.2e} V (互いに打ち消し)")

import numpy as np
import matplotlib.pyplot as plt

def build_rc_ladder(n_stages, R, C, omega):
    """
    n段RCラダー回路のコンダクタンス行列を生成
    各段: 直列抵抗R + 並列コンデンサC（GNDへ）
    """
    n_nodes = n_stages
    G = np.zeros((n_nodes, n_nodes), dtype=complex)
    Is = np.zeros(n_nodes, dtype=complex)

    Yr = 1.0 / R           # 抵抗のアドミタンス
    Yc = 1j * omega * C     # コンデンサのアドミタンス

    for i in range(n_nodes):
        # 自己アドミタンス
        G[i, i] += Yc  # GNDへのコンデンサ
        if i > 0:
            G[i, i] += Yr     # 左の抵抗
            G[i, i-1] -= Yr
            G[i-1, i] -= Yr
        if i == 0:
            G[i, i] += Yr     # 入力側の抵抗
        if i < n_nodes - 1:
            G[i, i] += Yr     # 右の抵抗

    # 入力: 電圧源 Vs = 1V → 節点0に Vs/R の電流を注入
    Is[0] = 1.0 / R  # Vs = 1V

    return G, Is

# パラメータ
n_stages = 10
R = 1e3      # 1 kΩ
C = 10e-9    # 10 nF

# 周波数応答を計算
f_range = np.logspace(1, 6, 500)
transfer_function = np.zeros(len(f_range), dtype=complex)

for idx, f in enumerate(f_range):
    omega = 2 * np.pi * f
    G, Is = build_rc_ladder(n_stages, R, C, omega)
    v = np.linalg.solve(G, Is)
    transfer_function[idx] = v[-1]  # 最終段の電圧

# ボード線図
magnitude_dB = 20 * np.log10(np.abs(transfer_function))
phase_deg = np.degrees(np.angle(transfer_function))

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# ゲイン
axes[0].semilogx(f_range, magnitude_dB, 'b-', linewidth=2)
axes[0].set_ylabel('Gain [dB]')
axes[0].set_title(f'{n_stages}-Stage RC Ladder Network (R={R/1e3:.0f}kΩ, C={C*1e9:.0f}nF)')
axes[0].grid(True, alpha=0.3, which='both')
axes[0].axhline(y=-3, color='r', linestyle='--', alpha=0.5, label='-3 dB')

# -3dBカットオフ周波数の推定
idx_3dB = np.argmin(np.abs(magnitude_dB - (-3)))
f_3dB = f_range[idx_3dB]
axes[0].axvline(x=f_3dB, color='r', linestyle=':', alpha=0.5)
axes[0].text(f_3dB * 1.2, -3 + 2, f'$f_{{-3dB}}$ = {f_3dB:.0f} Hz', color='r', fontsize=10)
axes[0].legend()

# 理想的な1段RCの周波数応答（比較用）
f_1stage = 1 / (2 * np.pi * R * C)
H_1stage = 1 / (1 + 1j * f_range / f_1stage)
axes[0].semilogx(f_range, 20 * np.log10(np.abs(H_1stage)), 'g--', linewidth=1.5,
                  label='1-stage RC', alpha=0.7)
axes[0].legend()

# 位相
axes[1].semilogx(f_range, phase_deg, 'r-', linewidth=2)
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Phase [degrees]')
axes[1].grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('rc_ladder_bode.png', dpi=150, bbox_inches='tight')
plt.show()

# 行列サイズとスケーラビリティ
print(f"\n=== 行列のスケーラビリティ ===")
for n in [5, 10, 20, 50, 100]:
    import time
    omega = 2 * np.pi * 1000
    G_test, Is_test = build_rc_ladder(n, R, C, omega)
    t_start = time.perf_counter()
    v_test = np.linalg.solve(G_test, Is_test)
    t_end = time.perf_counter()
    print(f"  {n:3d}段: 行列サイズ {n}x{n}, 計算時間 {(t_end-t_start)*1e6:.1f} μs")

print(f"\n=== カットオフ周波数 ===")
print(f"  1段RC: f_c = {f_1stage:.0f} Hz")
print(f"  {n_stages}段ラダー: f_-3dB ≈ {f_3dB:.0f} Hz")