#%%
#Libraries
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import matplotlib 
matplotlib.rcParams.update({'font.size': 13, 'lines.linewidth': 2.5})
from matplotlib.widgets import Slider, Button
#%%
#Explicit equations
alpha = 0.0162
k = 5
CA0 = 0.2
FA0 = 2.5
epsilon=-0.5

def ODEfun(Yfuncvec, W, alpha,k, CA0, FA0, epsilon):
    
    X= Yfuncvec[0]
    p= Yfuncvec[1]
    CA =CA0*(1-X)*p/(1+epsilon*X)
    rA=-k*CA*CA
    # Differential equations
    dXdW = -rA/FA0
    dpdW =-alpha/2/p*(1+epsilon*X) 
    return np.array([dXdW, dpdW])

Wspan = np.linspace(0, 100, 10000) # Range for the independent variable
y0 = np.array([0, 1]) # Initial values for the dependent variables
#%%

fig, ax = plt.subplots()
plt.subplots_adjust( left  = 0.4)
p1, p2 = plt.plot(Wspan, odeint(ODEfun, y0, Wspan, (alpha,k, CA0, FA0, epsilon)))
plt.legend(['X', 'p'], loc='best')
plt.title('Conversion and pressure drop Profile')
plt.ylim(0,1)
plt.xlim(0,100)
ax.grid()
ax.set_xlabel('W (kg)', fontsize='medium')
ax.set_ylabel('X, p', fontsize='medium')
# Slider Code
fig.suptitle('Example: PBR Reactor with Pressure drop \n  Note: If the graph becomes unstable, it means pressure drop is so high that flow can not happen beyond a particular catalyst weight, try changing k values',fontsize=12, fontweight='bold', x = 0.5, y= 0.98)
# Slider Code

ax.text(-50, 0.1,'Differential Equations'
         '\n\n'
         r'$\dfrac{d{X}}{dW} = \dfrac{-r_A^\prime}{F_{A0}}$'
                  '\n\n'
         r'$\dfrac{dp}{dW} = -\dfrac{\alpha(1+\epsilon X)}{2p}$'
                  '\n\n'
         'Explicit Equations' '\n\n'
         r'$\epsilon = -0.5$' '\n\n'
         r'$C_A = \dfrac{C_{A0}(1-X)p}{1+\epsilon X}$' '\n\n'
         r'$r_A^\prime = -kC_A^2$' '\n\n'
        , ha='left', wrap = True, fontsize=12,
        bbox=dict(facecolor='none', edgecolor='black', pad=10), fontweight='bold')

#%%
axcolor = 'black'
ax_alpha = plt.axes([0.1, 0.75, 0.2, 0.02], facecolor=axcolor)
ax_k = plt.axes([0.1, 0.7, 0.2, 0.02], facecolor=axcolor)
ax_CA0 = plt.axes([0.1, 0.65, 0.2, 0.02], facecolor=axcolor)
ax_FA0 = plt.axes([0.1, 0.6, 0.2, 0.02], facecolor=axcolor)

salpha = Slider(ax_alpha, r'$\alpha$ (kg$^{-1}$)', 0.000001, 0.017, valinit=.0162)
sk = Slider(ax_k, r'k ($\frac{dm^6}{mol.s.kg}$)', 2, 20, valinit=5)
sCA0= Slider(ax_CA0, r'$C_{A0}$ ($\frac{mol}{dm^3}$)', 0.2, 1, valinit=0.2)
sFA0 = Slider(ax_FA0, r'$F_{A0}$ ($\frac{mol}{s}$)', 0.5, 20, valinit=2.5)

def update_plot(val):
    alpha = salpha.val
    k =sk.val
    CA0 =sCA0.val
    FA0 = sFA0.val
    x = odeint(ODEfun, y0, Wspan, (alpha,k, CA0, FA0, epsilon))
    p1.set_ydata(x[:,0])
    p2.set_ydata(x[:,1])
    fig.canvas.draw_idle()

salpha.on_changed(update_plot)
sk.on_changed(update_plot)
sCA0.on_changed(update_plot)
sFA0.on_changed(update_plot)

resetax = plt.axes([0.15, 0.8, 0.09, 0.05])
button = Button(resetax, 'Reset variables', color='cornflowerblue', hovercolor='0.975')


def reset(event):
    salpha.reset()
    sk.reset()
    sCA0.reset()
    sFA0.reset()
   
button.on_clicked(reset)

plt.show()