#%%
#Libraries
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams.update({'font.size': 13, 'lines.linewidth': 2.5})
from matplotlib.widgets import Slider, Button

#%%
#Explicit equations
UA=16500;
FA0=0.0376;
FI0=0;
CpA=163;
CpB=83;
CT0=18.8;
mc=0.111;
Cpcoolant=34.5;
CpC=71;
To=1035;
CpI=32.68;
Tao=995.15
Ea=284.52
A=8.198*10**14
def ODEfun(Yfuncvec, V, UA,FA0,FI0,CpA,CpB,CT0,mc,Cpcoolant,CpC,To,CpI,Tao,Ea,A):
    Ta= Yfuncvec[0]
    T= Yfuncvec[1]
    X= Yfuncvec[2]
        #Explicit Equation Inline
    FT0=FI0+FA0;
    vo=(FT0/CT0)*1000;
    Ca0=(FA0/vo)*1000;
    delCp=CpB+CpC-CpA;
    thetaI=FI0/FA0;
    sumcp=CpA+thetaI*CpI;
    deltaH = 80770 + delCp * (T - 298); 
    k = A*np.exp(-Ea/(0.008314*T)); 
    ra = 0 - k*(Ca0 * (1 - X) * (To / T)) / (1 + X); 
     # Differential equations
    dTadV = -UA * (T - Ta) / mc / Cpcoolant; 
    dTdV = (ra * deltaH - (UA * (T - Ta))) / (FA0 * (sumcp + X * delCp)); 
    dXdV = 0 - (ra / FA0); 
    return np.array([dTadV,dTdV,dXdV])

Vspan = np.linspace(0, 0.001, 1000) # Range for the independent variable
y0 = np.array([Tao,1035,0]) # Initial values for the dependent variables

#%%
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
fig.suptitle("""LEP-12-2d:Production of Acetic Anhydride (Counter-Current)""", fontweight='bold', x = 0.2,y=0.99)
plt.subplots_adjust(left  = 0.48)
fig.subplots_adjust(wspace=0.25,hspace=0.3)
sol = odeint(ODEfun, y0, Vspan, (UA,FA0,FI0,CpA,CpB,CT0,mc,Cpcoolant,CpC,To,CpI,Tao,Ea,A))
Ta =sol[:, 0]
T = sol[:, 1]
X = sol[:, 2]
FT0=FI0+FA0;
vo=(FT0/CT0)*1000;
Ca0=(FA0/vo)*1000;
delCp=CpB+CpC-CpA;
thetaI=FI0/FA0;
sumcp=CpA+thetaI*CpI;
deltaH = 80770 + delCp * (T - 298); 
k = A*np.exp(-Ea/(0.008314*T)); 
ra = 0 - k*(Ca0 * (1 - X) * (To / T)) / (1 + X);
rate = 0 - ra;
Qg=ra * deltaH;
Qr=UA*(T-Ta);
p1,p2= ax1.plot(Vspan,T,Vspan,Ta)
ax1.legend(['T','$T_a$'], loc='upper left')
ax1.set_xlabel(r'$Volume  {(m^3)}$', fontsize='medium')
ax1.set_ylabel('Temperature (K)', fontsize='medium')
ax1.set_ylim(900,1300)
ax1.set_xlim(0,0.001)
ax1.grid()
ax1.ticklabel_format(style='sci',scilimits=(0,0),axis='x')

p3 = ax2.plot(Vspan,X)[0]
ax2.legend(['X'], loc='upper left')
ax2.set_ylim(0,1)
ax2.set_xlim(0,0.001)
ax2.grid()
ax2.set_xlabel(r'$Volume  {(m^3)}$', fontsize='medium')
ax2.set_ylabel('Conversion', fontsize='medium')
ax2.ticklabel_format(style='sci',scilimits=(0,0),axis='x')

p5,p6 = ax3.plot(Vspan,Qg,Vspan,Qr)
ax3.legend(['$Q_g$','$Q_r$'], loc='lower right')
ax3.set_ylim(-10*10**6,2*10**6)
ax3.set_xlim(0,0.001)
ax3.grid()
ax3.set_xlabel(r'$Volume  {(m^3)}$', fontsize='medium')
ax3.set_ylabel(r'$Q {(J/m^3.s)}$', fontsize='medium')
ax3.ticklabel_format(style='sci',scilimits=(3,4),axis='y')
ax3.ticklabel_format(style='sci',scilimits=(0,0),axis='x')

p4 = ax4.plot(Vspan, rate)[0]
ax4.legend(['$-r_A$'], loc='upper right')
ax4.set_ylim(0,100)
ax4.set_xlim(0,0.001)
ax4.grid()
ax4.set_xlabel(r'$Volume  {(m^3)}$', fontsize='medium')
ax4.set_ylabel('$Rate {(mol/m^3.s)}$', fontsize='medium')
ax4.ticklabel_format(style='sci',scilimits=(0,0),axis='x')

#ax1.axis('off')
ax1.text(-0.0025,450,'Note: While we used the expression k=$k_1$*exp(E/R*(1/$T_1$ - 1/$T_2$)) \n         in the textbook, in python we have to use k=A*exp(-E/RT) \n          in order to explore all the variables.',wrap = True, fontsize=13,
        bbox=dict(facecolor='none', edgecolor='red', pad=10))
ax1.text(-0.0025,580,'Differential Equations'
         '\n\n'
         r'$\dfrac{dT_a}{dV} = \dfrac{-Ua*(T-T_a)}{m_c*C_{P_{coolant}}} $'
         '\n'
         r'$\dfrac{dT}{dV} = \dfrac{r_{A}\Delta H_{Rx}-Ua*(T-T_a)}{F_{A0}.\left(\sum_{i}\theta_iC_{P_i}+X \Delta C_P\right)}$'
         '\n'
         r'$\dfrac{dX}{dV} = \dfrac{-r_{A}}{F_{A0}}$'
                  '\n \n'
                  
         'Explicit Equations'
                  '\n\n'  
         r'$A = 8.198*10^{14}\thinspace s^{-1}$'
         '\n'
         r'$F_{A0}=0.0376$'
         '\n'
         r'$F_{T0}=F_{A0}+F_{I0}$' 
         '\n'
          r'$v_{0}=\dfrac{F_{T0}}{C_{T0}}$' 
         '\n'
         r'$C_{A0}=\dfrac{F_{A0}}{v_{0}}$' 
         '\n'        
         r'$\theta_I = \dfrac{F_{I0}}{F_{A0}}}$'
                  '\n'
         r'$\sum_{i}\theta_iC_{pi} = C_{P_A} +\theta_IC_{P_I} $'
         '\n'
         r'$\Delta H_{Rx}=80770+ \Delta C_P* (T-298) $'
         '\n'
         r'$\Delta C_P= C_{P_B}+C_{P_C}-C_{P_A}$'
         '\n'
         r'$k = A*exp\left(\dfrac{-E}{0.00831*T}\right)$'
         '\n'
         r'$r_A = -\dfrac{kC_{A0}(1-X)}{1+X} \dfrac{T_0}{T}$'
                  '\n'
         r'$rate = -r_A$'
                 '\n'
         r'$Q_g = r_{A}.\Delta H_{Rx}$'
                  '\n'
         r'$Q_r = Ua(T-T_a)$', ha='left', wrap = True, fontsize=12,
        bbox=dict(facecolor='none', edgecolor='black', pad=10.0), fontweight='bold')
    
ax1.text(-0.0025, 290,'Instructions'
         '\n'
         'Guess a value of $T_a$ at V=0 and see if it matches $T_a$=1250 K at V=0.001 ${m^3}$'
         '\n'
        'So, Whenever you change any parameter,you also need to vary $T_{a0}$ to match '
        '\n'
        '$T_a$=1250 K at V=0.001 ${m^3}$', ha='left', wrap = True, fontsize=12,
        bbox=dict(facecolor='none', edgecolor='black', pad=10.0))

axcolor = 'black'
ax_FA0 = plt.axes([0.24, 0.84, 0.15, 0.015], facecolor=axcolor)
ax_FI0 = plt.axes([0.24, 0.80, 0.15, 0.015], facecolor=axcolor)
ax_CpA = plt.axes([0.24, 0.76, 0.15, 0.015], facecolor=axcolor)
ax_CpB = plt.axes([0.24, 0.72, 0.15, 0.015], facecolor=axcolor)
ax_CT0 = plt.axes([0.24, 0.68, 0.15, 0.015], facecolor=axcolor)
ax_CpC = plt.axes([0.24, 0.64, 0.15, 0.015], facecolor=axcolor)
ax_mc = plt.axes([0.24, 0.60, 0.15, 0.015], facecolor=axcolor)
ax_UA = plt.axes([0.24, 0.56, 0.15, 0.015], facecolor=axcolor)
ax_Tao = plt.axes([0.24, 0.52, 0.15, 0.015], facecolor=axcolor)
ax_Ea = plt.axes([0.24, 0.48, 0.15, 0.015], facecolor=axcolor)

sFA0 = Slider(ax_FA0, r'F$_{A0}$($\frac{mol}{s}$)', 0.01, 1, valinit=0.0376,valfmt='%1.3f')
sFI0 = Slider(ax_FI0, r'F$_{I0}$ ($\frac{mol}{s}$)', 0, 10, valinit=0,valfmt='%1.1f')
sCpA= Slider(ax_CpA, r'C$_{P_A}$ ($\frac{J}{mol.k}$)', 50, 200, valinit=163,valfmt='%1.1f')
sCpB = Slider(ax_CpB, r'C$_{P_B}$ ($\frac{J}{mol.k}$)',10, 200, valinit=83,valfmt='%1.2f')
sCT0 = Slider(ax_CT0,r'C$_{T0}$ ($\frac{mol}{m^3}$)', 10, 40, valinit= 18.8,valfmt='%1.1f')
sCpC = Slider(ax_CpC, r'C$_{P_C}$ ($\frac{J}{mol. K}$)', 1, 200, valinit=71,valfmt='%1.1f')
smc = Slider(ax_mc, r'm$_{C}$ ($\frac{mol}{s}$)', 0.01, 0.2, valinit=0.111,valfmt='%1.3f')
sUA = Slider(ax_UA, r'Ua ($\frac{J}{mol^3.s.K}$)', 1000, 30000, valinit=16500,valfmt='%1.0f')
sTao = Slider(ax_Tao,'$T_{ao}$ (K)', 500, 1200, valinit=995.15,valfmt='%1.1f')
sEa= Slider(ax_Ea, r'$E$ ($\frac{kJ}{mol}$)', 250, 320, valinit=284.5,valfmt='%1.1f')

def update_plot2(val):
    FA0 = sFA0.val
    FI0 =sFI0.val
    CpA =sCpA.val
    CpB = sCpB.val
    CT0 =sCT0.val
    CpC = sCpC.val
    mc = smc.val
    UA = sUA.val
    Tao = sTao.val
    Ea= sEa.val
    y0 = np.array([Tao,1035,0])
    sol = odeint(ODEfun, y0, Vspan, (UA,FA0,FI0,CpA,CpB,CT0,mc,Cpcoolant,CpC,To,CpI,Tao,Ea,A))
    Ta = sol[:, 0]
    T = sol[:, 1]
    X = sol[:, 2]
    FT0=FI0+FA0;
    vo=(FT0/CT0)*1000;
    Ca0=(FA0/vo)*1000;
    delCp=CpB+CpC-CpA;
    deltaH = 80770 + delCp * (T - 298); 
    k = A*np.exp(-Ea/(0.008314*T)); 
    ra = 0 - k*(Ca0 * (1 - X) * (To / T)) / (1 + X);
    rate = 0 - ra;   
    Qg=ra * deltaH;
    Qr=UA*(T-Ta);
    p1.set_ydata(T)
    p2.set_ydata(Ta)
    p3.set_ydata(X)
    p4.set_ydata(rate)
    p5.set_ydata(Qg)
    p6.set_ydata(Qr)
    fig.canvas.draw_idle()

sFA0.on_changed(update_plot2)
sFI0.on_changed(update_plot2)
sCpA.on_changed(update_plot2)
sCpB.on_changed(update_plot2)
sCT0.on_changed(update_plot2)
sCpC.on_changed(update_plot2)
smc.on_changed(update_plot2)
sUA.on_changed(update_plot2)
sTao.on_changed(update_plot2)
sEa.on_changed(update_plot2)

resetax = plt.axes([0.26, 0.88, 0.09, 0.04])
button = Button(resetax, 'Reset variables', color='cornflowerblue', hovercolor='0.975')

def reset(event):
    sFA0.reset()
    sFI0.reset()
    sCpA.reset()
    sCpB.reset()
    sCT0.reset()
    sCpC.reset()
    smc.reset()
    sUA.reset()
    sTao.reset()
    sEa.reset()
button.on_clicked(reset)