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

T = 5              # cas
dt = 0.0001         # delka kroku
Nt = int(T/dt)     # pocet casovych intervalu

# hodnoty soucastek z prik. 3.6, viz str. 60
C = 0.001
R = 5
R_Z = 15
L = 2
U_A = 1
f = 1
I_0 = 0.000001
n = 1
U_T = 0.0257


U_out = np.zeros(Nt+1) # pole, do ktereho budeme ukladat vypoctene hodnoty
U_out_d = np.zeros(Nt+1)   # derivace U_out
t = np.linspace(0, T, Nt+1) # casova promenna
U_in = U_A * np.sin(2*m.pi*f*t) # vstupni signal
U_out[0] = 0  # poc. podminka 1
U_out_d[0] = 0 # poc. podminka 2

# logaritmus s rozsirenym def. oborem - v pripade nestability kolem nuly
def log_ext(x):
        if x < 10**(-30):
                return -69
        else:
                return np.log(x)


for i in range(0, Nt):
        """metoda konecnych diferenci, viz vztah (3.49) na str. 58"""
        z = (U_out[i])/(I_0*R_Z) + (C/I_0)*U_out_d[i] + 1
        s_1 = -((2*n*dt*U_T)/(L*C)) * log_ext(z)
        s_2 = -(dt/(C*R_Z) - 1 + (dt*R)/L)*U_out_d[i]
        s_3 = -( (dt*R)/(L*C*R_Z) + dt/(L*C) )*U_out[i]
        s_4 = (dt/(L*C))*abs(U_in[i])
        U_out_d[i+1] = s_1 + s_2 + s_3 + s_4
        U_out[i+1] = U_out_d[i]*dt + U_out[i]

# graficke zobrazeni
plt.plot(t, U_in, label=r'$U_{in}$') # vstupni zdroj
plt.plot(t, U_out, label=r'$U_{out}$') # vystupni napeti
plt.ylabel('u [V]', rotation=0)
plt.xlabel('čas [s]')
plt.legend()
plt.show()