import numpy as np
import matplotlib.pyplot as plt

# Заданные параметры (их можно изменять)
R1 = 10
R2 = 40
R3 = 10
C_value = 200e-6  # Ф
w = 1000  # угловая частота для мгновенных значений, рад/с

# Входное напряжение (амплитуда и фаза в радианах)
Uent_rms = 5.0
Uent_phase = -0.2444  # радиан

# Перевод в комплексную форму амплитуды (Effective Value) -> Комплексная амплитуда
# Важно: входная величина Uent_rms уже среднеквадратичная, так что комплексное значение:
Uent = Uent_rms * np.exp(1j * Uent_phase)

# Импеданс конденсатора
Z_C = 1/(1j * w * C_value)

# Расчет эквивалентного входного импеданса параллельной части
ZII = (R2*(R3+Z_C)) / (R2 + R3 + Z_C)
Zent = R1 + ZII

# Входной ток
ient = Uent / Zent

# Напряжение на параллельной части
UII = ient * ZII

# Ток через R1
i1 = (Uent - UII) / R1

# Токи через элементы параллельной ветви
i2 = UII / R2
i3 = UII / (R3 + Z_C)

# Выходное напряжение
Uexit = i3 * (R3 + Z_C)

# Отношение амплитуд и фазовый сдвиг
abs_Uent = np.abs(Uent)
abs_Uexit = np.abs(Uexit)
rel = abs_Uexit / abs_Uent
relang = np.angle(Uexit) - np.angle(Uent)

# Вывод промежуточных и конечных результатов
print("Исходные данные:")
print(f"R1 = {R1} Ом, R2 = {R2} Ом, R3 = {R3} Ом, C = {C_value} Ф")
print(f"w = {w} рад/с")
print(f"U_вх(rms) = {Uent_rms} В, фаза(U_вх) = {Uent_phase} рад")
print("\nРассчитанные импедансы и токи:")
print(f"Z_C = {Z_C} Ом")
print(f"Z_II = {ZII} Ом")
print(f"Z_вх = Z_ent = {Zent} Ом")
print(f"i_вх = i_ent = {ient} А")
print(f"UII = {UII} В")
print(f"i_1 = {i1} А")
print(f"i_2 = {i2} А")
print(f"i_3 = {i3} А")
print("\nВыходное напряжение:")
print(f"U_вых = {Uexit} В")
print(f"|U_вых| = {abs_Uexit} В, фаза(U_вых) = {np.angle(Uexit)} рад")
print("\nОтношение амплитуд:")
print(f"|U_вых|/|U_вх| = {rel}")
print(f"Сдвиг фаз между выходным и входным напряжениями: {relang} рад")

# Мгновенные значения для сравнения
T = 2 * np.pi / w
tt = np.linspace(0, 2*T, 1000)
# Переводим rms величины обратно в амплитудные для мгновенных значений
# u(t) = U_rms * sqrt(2) * sin(...)
uent_time = Uent_rms*np.sqrt(2)*np.sin(w*tt + np.angle(Uent))
uexit_time = abs_Uexit*np.sqrt(2)*np.sin(w*tt + np.angle(Uexit))

# Построение временных диаграмм
plt.figure(figsize=(10, 6))
plt.plot(tt, uent_time, label="Входное напряжение $u_{вх}(t)$")
plt.plot(tt, uexit_time, label="Выходное напряжение $u_{вых}(t)$")
plt.xlabel("Время $t$, с")
plt.ylabel("Напряжение, В")
plt.title("Временные диаграммы входного и выходного напряжений")
plt.legend()
plt.grid(True)
plt.savefig('images/vrem_diagrams.png', dpi=300)
plt.close()

# Передаточная функция как функция s
def W(s):
    ZC = 1/(s*C_value)
    ZII = (R2*(R3+ZC))/(R2+R3+ZC)
    Zent = R1 + ZII
    return ZII/Zent

# Частотные характеристики
www = np.linspace(0.1, 5000, 500)  # частоты от 0.1 до 5000 рад/с
Ww = np.array([W(1j*om) for om in www])
Aw = np.abs(Ww)
Fw = np.angle(Ww)

# АЧХ
plt.figure(figsize=(10, 6))
plt.plot(www, Aw, label="АЧХ")
plt.xlabel("Частота, рад/с")
plt.ylabel("Амплитуда")
plt.title("Амплитудно-частотная характеристика")
plt.grid(True)
plt.legend()
plt.savefig('images/achx.png', dpi=300)
plt.close()

# ФЧХ
plt.figure(figsize=(10, 6))
plt.plot(www, Fw, label="ФЧХ")
plt.xlabel("Частота, рад/с")
plt.ylabel("Фаза, рад")
plt.title("Фазово-частотная характеристика")
plt.grid(True)
plt.legend()
plt.savefig('images/fchx.png', dpi=300)
plt.close()

# Вычисление специальных точек для годографа
w0 = 0.001     # Приближение к нулевой частоте
W0 = W(1j*w0)     
W1000 = W(1j*w) # ω=1000 рад/с
Winf = W(1j*1e6)   # Приближение для очень большой частоты (ω→∞)

print("\nСпециальные точки передаточной функции:")
print(f"W(j0) ≈ W(j{w0}) = {W0}")
print(f"W(j1000) = {W1000}")
print(f"W(j∞) ≈ W(j1e6) = {Winf}")

# Вычисление специальных точек для годографа
w0 = 0.001     # Приближение к нулевой частоте
W0 = W(1j*w0)     
W1000 = W(1j*w) # ω=1000 рад/с
Winf = W(1j*1e6)   # Приближение для очень большой частоты (ω→∞)

print("\nСпециальные точки передаточной функции:")
print(f"W(j0) ≈ W(j{w0}) = {W0}")
print(f"W(j1000) = {W1000}")
print(f"W(j∞) ≈ W(j1e6) = {Winf}")

A0, F0 = np.abs(W0), np.angle(W0)
A1000, F1000 = np.abs(W1000), np.angle(W1000)
Ainf, Finf = np.abs(Winf), np.angle(Winf)

# Годограф
# Построение полярного графика
plt.figure(figsize=(8,8))
ax = plt.subplot(111, projection='polar')

# Построение годографа
angles = np.angle(Ww)
magnitudes = np.abs(Ww)
ax.plot(angles, magnitudes, label="Годограф $W(j\omega)$")

# Отметка специальных точек
ax.plot(np.angle(W0), np.abs(W0), 'ro', label=f'ω=0 рад/с')
ax.plot(np.angle(W1000), np.abs(W1000), 'go', label=f'ω={w} рад/с')
# Для ω → ∞ можно использовать предел или выбрать очень большое ω
# Например, ω = 1e6
W_large = W(1j*1e6)
ax.plot(np.angle(W_large), np.abs(W_large), 'bo', label='ω→∞')

# Настройка заголовка и легенды
ax.set_title("Годограф передаточной функции $W(j\omega)$", va='bottom')
ax.legend(loc='upper right')

# Сохранение графика
plt.savefig('images/godograf_polar.png', dpi=300)
plt.close()