import matplotlib.pyplot as plt
import numpy as np

H = 0.079
D = 0.090
c = 4180
rho = 1000
k_tc = 0.2
L = 3e-3
T_room = 20

V = H*np.pi*(D/2)**2
C = c*rho*V
A = np.pi*D*H + 2*np.pi*(D/2)**2
G = (k_tc/L)*A

dt = 0.1
t_start = 0
t_stop = 200
N_sim = int((t_stop - t_start)/dt) + 1

P_on = 1777
P_off = 0

t_array = np.zeros(N_sim)
T_array = np.zeros(N_sim)
T_room_array = np.zeros(N_sim) + T_room
P_array = np.zeros(N_sim)

T_min = 0
T_max = 100

T_k = T_init = 20.0

for k in range(0, N_sim):

    t_k = k*dt

    if (0 <= t_k <= t_start): P_k = 0
    else: P_k = P_on
 
    dT_dt_k = (1/C)*(P_k + G*(T_room - T_k))
    T_kp1 = T_k + dt*dT_dt_k
    T_kp1 = np.clip(T_kp1, T_min, T_max)     

    t_array[k] = t_k
    T_array[k] = T_k
    P_array[k] = P_k

    T_k = T_kp1

T_100 = 100 # [oC]
indexes_less_than_100 = np.argwhere(T_array < T_100)
t_100 = t_array[np.amax(indexes_less_than_100)]
print('t_100 [s] =', f'{t_100:.2f}')

plt.close('all')
plt.figure(num=1, figsize=(12, 9))

plt.subplot(2, 1, 1)
plt.plot(t_array, T_array,'b', label='T')
plt.plot(t_array, T_room_array, 'g', label='T_room')
plt.legend()
plt.ylabel('[deg C]')
plt.grid()

plt.subplot(2, 1, 2)
plt.plot(t_array, P_array, 'r', label='P')
plt.legend()
plt.grid()
plt.xlabel('t [s]')
plt.ylabel('[W]')

# plt.savefig('prog_sim_kettle_detect_p.pdf')
plt.show()