import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import sympy as sp

# Constants
L = 1  # inductance
R = 1.5  # resistance
t_final = 2.2  # final time
omega = 50 * 2 * np.pi  # 50 Hz angular frequency
ph0 = 0  # initial phase
V_peak = 120  # peak voltage

# Symbolic variables and functions
t = sp.symbols('t', real=True)
i = sp.Function('i')(t)
V_t = V_peak * sp.sin(omega * t + ph0)

# ODE: di/dt = (V(t) - R*i(t)) / L
ode = sp.Eq(i.diff(t), (V_t - R * i) / L)
sol = sp.dsolve(ode, ics={i.subs(t, 0): 0})

# Extract solution as a lambda function for plotting
sol_func = sp.lambdify(t, sol.rhs, 'numpy')

# Calculate the settling time (5 percent settling time)
circuit_pole = -R / L
t_settling = -3/circuit_pole  # 5% settling time, approx
print(t_settling)

# Plot the transient solution
time_points = np.linspace(0, t_final, 1500)
i_values = sol_func(time_points)

plt.figure(figsize=(10, 6))
plt.plot(time_points, i_values, label="Transient", color='blue')
plt.axhline(0, color='black', linewidth=1)
plt.axvline(t_settling, color='green', linestyle='--', label="Settling Time (5%)")
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Current [A]')
plt.legend(loc="best")
plt.title('Transient Response of the Linear Inductor Circuit')
plt.show()

# No resistor case (R = 0)
# In this case, the ODE becomes: di/dt = V(t) / L
i_noR = sp.integrate(V_t / L, (t, 0, t))  # No resistor case
i_noR_simplified = sp.simplify(i_noR)

# Display symbolic result for no-resistor case
print(f"No resistor current i(t) = {i_noR_simplified}")

# Plot the no resistor solution alongside the transient solution
i_noR_func = sp.lambdify(t, i_noR_simplified, 'numpy')
i_noR_values = i_noR_func(time_points)

plt.figure(figsize=(10, 6))
plt.plot(time_points, i_values, label="Transient", color='blue')
plt.plot(time_points, i_noR_values, label="No Resistor", color='red', linestyle='--')
plt.axhline(0, color='black', linewidth=1)
plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Current [A]')
plt.legend(loc="best")
plt.title('Comparison of Transient and No Resistor Response')
plt.show()

# Steady-state sinusoidal regime (Z = R + jωL)
Z = R + L * omega * 1j
amplI_ss = abs(1 / Z) * V_peak  # Current amplitude
phaseI_ss = np.angle(1 / Z)  # Current phase

# Steady-state current
I_ss = amplI_ss * np.sin(omega * time_points + phaseI_ss + ph0)

# Plot the transient and steady-state solutions
plt.figure(figsize=(10, 6))
plt.plot(time_points, i_values, label="Transient", color='blue')
plt.plot(time_points, I_ss, label="Steady State", color='red', linestyle='--')
plt.axhline(0, color='black', linewidth=1)
plt.axvline(t_settling, color='green', linestyle='--', label="Settling Time (5%)")
for value in [-2 * amplI_ss, -amplI_ss, amplI_ss, 2 * amplI_ss]:
    plt.axhline(value, color='gray', linestyle=':', linewidth=1)

plt.grid(True)
plt.xlabel('Time [s]')
plt.ylabel('Current [A]')
plt.legend(loc="best")
plt.title('Comparison of Transient and Steady-State Response')
plt.show()

# Display the amplitude and phase of steady-state current
print(f"Steady-state current amplitude: {amplI_ss}")
print(f"Steady-state current phase (degrees): {np.degrees(phaseI_ss)}")