Untitled
raw download clone
TEXT
views 50
,
size 5080 b
import numpy as np
from math import *
import matplotlib.pyplot as plt
 
#Гармонический сигнал
def harm_signal(t):
    return sin(2.0*pi*fc*t)
 
#Производная гармонического сигнала
def d_harm_signal(t):
    return 2.0*pi*fc*cos(2.0*pi*fc*t)
 
#Перевод частоты в циклическую
def f2w(f):
    return 2.0*pi*f
 
def Z1(f, C1):   
    return 2.0/(1j*f2w(f)*C1)
 
def Z2 (f, C2):
    return 1.0/(1j*f2w(f)*C2)
 
def Z3(f, L):
    return 1.0j*f2w(f)*L
 
#Постоянная распространения отдельной ячейки
def Gam(f, L, C1, C2):
    ZY = (Z2(f, C2)+Z3(f, L))/Z1(f, C1)
    return 2.0 * np.arcsinh(np.sqrt(ZY))
 
#Характеристическое сопротивление отдельной ячейки
def Zw(f, L, C1, C2):
    return np.sqrt((Z1(f, C1)**2*(Z2(f, C2)+Z3(f, L)))/(2*Z1(f, C1)+Z2(f, C2)+Z3(f, L)))
 
def d2V1():
    return (1.0/(L*C1)*(aV[1]-aV[0]+aU[0])+1.0/(Z0*K0*C1)*(A0*d_harm_signal(time[it])-dV[0]))
 
def d2VN():
    return (1.0/(L*C1)*(aV[Nc-1]-aV[Nc]-aU[Nc-1])+1.0/(Z0*KN*C1)*(AN*d_harm_signal(time[it])-dV[Nc]))
 
def d2Vs():
    return (1.0/(L*C1)*(aV[ic-1]-2.0*aV[ic]+aV[ic+1]+aU[ic]-aU[ic-1]))
 
def d2Us():
    return (1.0/(L*C2)*(aV[ic]-aV[ic+1]-aU[ic])-G/C2*dU[ic])
 
global fc, L, C1, C2, G, aU, dU, aV, dV, time, A0, AN, K0, KN
 
fc = float(input('Частота сигнала возбуждения ЛП '))
Tc = float(input('Временной интервал '))
#nt = int(input('Число временных отсчетов '))
 
fl = float(input('Нижняя граничная частота ЛП '))
fh = float(input('Верхняя граничная частота ЛП '))
#f0 = (fl + fh) * 0.5
Z0 = float(input('Характеристическое сопротивление одного звена ЛП на частоте сигнала '+str(fc)+' '))
Nc = int(input('Число ячеек в ЛП '))
 
L = (sqrt(Z0**2*f2w(fc)**2*(2*f2w(fh)**2-f2w(fl)**2-f2w(fc)**2)/
    ((f2w(fh)**2-f2w(fl)**2)**2*(f2w(fc)**2-f2w(fl)**2))))
C1 = 2.0 / L / (f2w(fh)**2 - f2w(fl)**2)
C2 = 1.0 / (f2w(fl)**2 * L)
G = 0
 
npp = 20
dt = 1/(fh*npp)
num = int(Tc / dt)
 
print('Параметры отдельной ячейки ЛП:')
print('C1 = {0: f}\nC2 = {1: f}\nL = {2: f}'.format(C1, C2, L))
 
freq = np.linspace(0.8*fl, fh*1.2, num)
 
Gama = Gam(freq, L, C1, C2)
Zw = Zw(freq, L, C1, C2)
dF = (Gam(freq+0.1, L, C1, C2).imag-Gam(freq-0.1, L, C1, C2).imag)/0.2
 
#######
A0 = 1 #Амплитуда сигнала слева
AN = 0 #Амплитуда сигнала справа
K0 = KN = 1 #Коэффициенты при нагрузочных сопротивлениях
 
aU = [0] * Nc     #Массив напряжений на емкости C2
dU = [0] * Nc     #Массив производных напряжений на емкости C2
aV = [0] * (Nc+1) #Массив напряжений на емкости C1
dV = [0] * (Nc+1) #Массив производных напряжений на емкости C1
 
Vinp = [0] * num  #Массив входных напряжений
Vout = [0] * num  #Массив выходных напряжений
time = [0] * num  #Массив временных отсчетов
 
#Решение уравнений возбуждения линии передачи во времени
for it in range(num):
    time[it] = dt * it
    dV[0] += d2V1()*dt
    for ic in range (Nc):
        dU[ic] = d2Us()*dt
        if ic == 0:
            continue
        dV[ic] = d2Vs()*dt
    dV[Nc] += d2VN()*dt
 
    for ic in range(Nc):
        aV[ic] += dV[ic]*dt
        aU[ic] += dU[ic]*dt
    aV[Nc] += dV[Nc]*dt
 
    Vinp[it] = aV[0]
    Vout[it] = aV[Nc]
    if it % 100 == 0:
        print('{0: 7.3f} {1: 7.3f} {2: 7.3f} '.format(time[it], Vinp[it], Vout[it]))
spectr_inp = np.fft.fft(Vinp)
spectr_out = np.fft.fft(Vout)
fft_freq = np.fft.fftfreq(num, Tc/num)
 
#######
plt.plot(freq, Gama.real, color='tab:blue', label=r'$\alpha(f)$')
plt.tick_params(axis='y', labelcolor='tab:blue')
plt.legend(loc='lower right')
plt.twinx()
plt.plot(freq, Gama.imag, color='tab:orange', label=r'$\phi(f)$')
plt.tick_params(axis='y', labelcolor='tab:orange')
plt.legend(loc='upper left')
plt.show()
 
plt.plot(freq, abs(Zw), label='$|Z_0|(f)$')
plt.plot(freq, Zw.real, label='$Re(Z_0)(f)$')
plt.plot(freq, Zw.imag, label='$Im(Z_0)(f)$')
plt.vlines(fc, 0, Z0, color='tab:olive', linestyles='dashdot', lw=1)
plt.hlines(Z0, freq[0], fc, color='tab:olive', linestyles='dashdot', lw=1)
plt.legend(loc='best')
plt.show()
 
plt.plot(freq, dF)
plt.show()
 
plt.plot(time, Vinp, time, Vout)
plt.show()
 
plt.plot(fft_freq[0:num//2], (np.hypot(spectr_inp.real, spectr_inp.imag)/num*2)[0:num//2], label='$V_{inp}$')
plt.plot(fft_freq[0:num//2], (np.hypot(spectr_out.real, spectr_out.imag)/num*2)[0:num//2], label='$V_{out}$')
plt.legend(loc='best')
plt.show()
close fullscreen
Login or Register to edit or fork this paste. It's free.