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

import matplotlib
import matplotlib.pyplot as plt
from scipy import signal
import numpy as np

import ahkab
from ahkab import circuit, time_functions


# シンボル列を引き延ばす
def extend_symbol_to_simulation_sampling_rate(symbols,t_symbol,t_simulation_sampling):
    # シンボル数
    size_of_symbols=symbols.size
    # 最後のシンボルが終わる時間
    end_time=size_of_symbols*t_symbol
    # 計算タイミング
    t=np.arange(0,end_time,t_simulation_sampling)
    # 各シンボルの開始（終了）時間
    t_periods_of_symbols=np.append(0,np.arange(1,size_of_symbols,1)*t_symbol)
    # 結果データ領域確保
    ex_symbols=np.zeros_like(t)
    # 先頭データは入れておく
    ex_symbols[0]=symbols[0]
    # 各シンボルについて、その期間のインデックスを取得し、データを入れ込む
    for i in range(0,size_of_symbols):
        ex_symbols[np.where(t>t_periods_of_symbols[i])]=symbols[i]
    return ex_symbols,end_time

t_samp=0.5e-6
t_sym=100e-6
freq_carrier=459441
amplitude_carrier=5 # in [Vp]
symbols_in_base=np.array([0,1,0,1,0,1,0])
t_org=np.arange(0,t_sym*symbols_in_base.size,t_samp)
waveform_carrier=amplitude_carrier*np.cos(2*np.pi*freq_carrier*t_org)
[symbols,t_end]=extend_symbol_to_simulation_sampling_rate(symbols_in_base,t_sym,t_samp)
waveform_gen=waveform_carrier*symbols
vsource_pwl=ahkab.time_functions.pwl(t_org,waveform_gen,repeat=False,repeat_time=t_end,td=0)

R1_val=100e0
L1_val=600e-6
C1_val=200e-12
 
R2_val=100e3
L2_val=3e-3
C2_val=40e-12
 
Rg_val=10e6

k_val=0.00074

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))

print ('Resonance frequency=',end='')
print(1/(2*np.pi*np.sqrt(C2_val*L2_val)))
print ('Calculated Q=',end='')
print(R2_val*np.sqrt(C2_val/L2_val))

op_analysis = ahkab.new_op()
tran_analysis = ahkab.new_tran(tstart=0, tstop=t_end, tstep=t_samp, x0=None)

mycircuit = circuit.Circuit(title="mutual coupling circuit")
gnd = mycircuit.get_ground_node()
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_resistor("R2", n1="n4", n2="n5", value=R2_val)
mycircuit.add_inductor("L2", n1="n4", n2="n5", value=L2_val)
mycircuit.add_capacitor("C2", n1="n4", n2="n5", value=C2_val)
mycircuit.add_resistor("Rg", n1="n5", n2=gnd, value=Rg_val)
mycircuit.add_inductor_coupling("M1","L1","L2",k_val)
mycircuit.add_vsource("V1", n1="n1", n2=gnd, dc_value=0, ac_value=0, function=vsource_pwl)

#print(mycircuit)
r = ahkab.run(mycircuit, an_list=[op_analysis, tran_analysis])

t=r['tran']['T']
tI_R1=(r['tran']['Vn3']-r['tran']['Vn2'])/R1_val
tV_R2=r['tran']['Vn5']-r['tran']['Vn4']

matplotlib.use('TkAgg') 

fig = plt.figure()
ax1=plt.subplot(2,1,1)
ax1.grid(True)
ax1.plot(t, np.real(tI_R1), '-')
ax1.set_ylabel('(Vn3-Vn2)/R1 [A]')

ax2=plt.subplot(2,1,2)
ax2.grid(True)
ax2.plot(t, np.real(tV_R2), '-')
ax2.set_ylabel('(Vn5-Vn4) [V]')

ax2.set_xlabel('time [s]')

fig.tight_layout()
plt.show()