import matplotlib.pyplot as plt
import numpy as np

# ----------------------------------
# Function for Euler integration:
    
def fun_euler(T_k, P_k, params, T_limits, dt):
    C, G, T_room = params
    T_min, T_max = T_limits
    T_k = np.clip(T_k, T_min, T_max) 
    dT_dt_k = (1/C)*(P_k + G*(T_room - T_k))
    T_kp1 = T_k + dt*dT_dt_k
    return T_kp1
# ----------------------------------

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

# ----------------------------------
params = np.array([C, G, T_room])
# ----------------------------------

dt = 0.1
t_start = 10
t_stop = 400
N_sim = int((t_stop - t_start)/dt) + 1

P_on = 700
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_limits = np.array([T_min, T_max])
# ----------------------------------

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
 
    # ----------------------------------
    T_kp1 = fun_euler(T_k, P_k, params, T_limits, dt)
    # ----------------------------------
    
    t_array[k] = t_k
    T_array[k] = T_k
    P_array[k] = P_k

    T_k = T_kp1

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

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

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

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