ahkab(python)で回路シミュレーション(3)

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

import ahkab
from ahkab import circuit
mycircuit = circuit.Circuit(title="series resonant circuit")

gnd = mycircuit.get_ground_node()

R1_val=100e0
L1_val=600e-6
C1_val=200e-12

mycircuit.add_capacitor("C1", n1="n1", n2="n2", value=C1_val)
mycircuit.add_resistor("R1", n1="n2", n2="n3", value=R1_val)
mycircuit.add_inductor("L1", n1="n3", n2=gnd, value=L1_val)

mycircuit.add_vsource("V1", n1="n1", n2=gnd, dc_value=0, ac_value=-1)

print(mycircuit)

op_analysis = ahkab.new_op()
ac_analysis = ahkab.new_ac(start=300e3, stop=600e3, points=1001)

r = ahkab.run(mycircuit, an_list=[op_analysis, ac_analysis])

print ('Resonance frequency=',end='')
print(1/(2*np.pi*np.sqrt(C1_val*L1_val)))
print ('Calculated Q=',end='')
print((1/R1_val)*np.sqrt(L1_val/C1_val))

f=r['ac']['f']
I_R1=(r['ac']['Vn3']-r['ac']['Vn2'])/R1_val

I_R1_abs=np.abs(I_R1)
I_R1_max=np.max(I_R1_abs)
f_band=f[I_R1_abs>=(I_R1_max/np.sqrt(2))]
q_sim=(f[I_R1_abs==I_R1_max]/(f_band[-1]-f_band[0]))[0]
print ('Q from waveform=',end='')
print(q_sim)

matplotlib.use('TkAgg') 

fig = plt.figure()
ax1=plt.subplot(2,1,1)
ax1.grid(True)
ax1.plot(f, np.abs(I_R1), '-')
ax1.set_yscale('log')
ax1.set_ylabel('abs(I@R1) [I]')
ax1.set_title(mycircuit.title + " - AC Simulation")
ax2=plt.subplot(2,1,2)
ax2.grid(True)
ax2.plot(f, np.angle(I_R1), '-')
ax2.set_xlabel('frequency [Hz]')
ax2.set_ylabel('arg(I@R1) [rad]')
fig.tight_layout()
plt.show()