import matplotlib.pyplot as plt
import numpy as np

F = 0.2
V = 1

cA_min = 0
cA_max = np.inf

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

t_array = np.zeros(N_sim)
cA_array = np.zeros(N_sim)
wA_array = np.zeros(N_sim)

cA_k = cA_init = 0

for k in range(0, N_sim):

    t_k = k*dt
    
    # State lim:
    cA_k = np.clip(cA_k, cA_min, cA_max)

    if (0 <= t_k <= 10): wA_k = 0
    else: wA_k = 0.8
    
    dcA_dt_k = (1/V)*(wA_k - cA_k*F)
    cA_kp1 = cA_k + dt*dcA_dt_k
    
    t_array[k] = t_k
    cA_array[k] = cA_k
    wA_array[k] = wA_k

    cA_k = cA_kp1

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

plt.subplot(2, 1, 1)
plt.plot(t_array, cA_array, 'b', label='cA')
plt.legend()
plt.xlabel('t [min]')
plt.ylabel('[kg/m3]')
plt.grid()

plt.subplot(2, 1, 2)
plt.plot(t_array, wA_array, 'r', label='wA')
plt.legend()
plt.grid()
plt.xlabel('t [min]')
plt.ylabel('[kg/min]')

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