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

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

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

gnd = mycircuit.get_ground_node()

R1_val=100e3
L1_val=3e-3
C1_val=40e-12

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

mycircuit.add_isource("I1", 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(R1_val*np.sqrt(C1_val/L1_val))

f=r['ac']['f']
V_n1=r['ac']['Vn1']

V_n1_abs=np.abs(V_n1)
V_n1_max=np.max(V_n1_abs)
f_band=f[V_n1_abs>=(V_n1_max/np.sqrt(2))]
q_sim=(f[V_n1_abs==V_n1_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(V_n1), '-')
ax1.set_yscale('log')
ax1.set_ylabel('abs(V(n1)) [V]')
ax1.set_title(mycircuit.title + " - AC Simulation")
ax2=plt.subplot(2,1,2)
ax2.grid(True)
ax2.plot(f, np.angle(V_n1), '-')
ax2.set_xlabel('frequency [Hz]')
ax2.set_ylabel('arg(V(n1)) [rad]')
fig.tight_layout()
plt.show()