# Loading our favorite python libraries

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

from __future__ import division

plt.style.use('fivethirtyeight')

from IPython.html.widgets import interact
from IPython.display import clear_output, display, HTML

def STreactor(C_A0 = 10.2, C_B0 = 0, C_C0 = 0, cat_load = 0.2):
    ''' Batch STR reactor model
	initial conditions (concentrations of A, B and C)
	should be given in [mol*L-1]
    
    catalyst load is given in %
	'''
	
    C_A0 = C_A0 # mol*L-1
    C_B0 = C_B0 # mol*L-1
    C_C0 = C_C0 # mol*L-1
	
	
    C_D0 = 0     # mol*L-1
    C_H0 = 19.2 # mol*L-1

    mtot = 200 # g of A in a semibatch reactor

    # Molecular weights
    MW = [84.118, 86.132, 88.148, 160, 1.008]  # g mol-1

    dens = 0.86 # g.mL-1 average denstiy of A
    # Note:  Assumes similar densities of all the compounds    
    # Reaction volume
    VL = mtot/MW[0]/C_A0 # in L


    To = 348 # K - Initial temperature (75ºC)


    # Catalyst
    CatalystMass = 0.3829 # in [g]
    CatAvalaible = cat_load/100 # Percent 

    CatMW = 106.42

    nCat = (CatalystMass*CatAvalaible)/CatMW


    # System reaction paramiters
    # Matrix with stoichiometry coefficients 
    # i reactions, j components (RxS)
    # j components in this order:  1    2   3   4   5     
    #                               A   B   C   D   H
    alfa = np.array([
            [-1, 1, 0, 0, -1],
            [0, -1, 1, 0, -1],
            [-1, 0, 1, 0, -2],
            [-2, 0, 0, 1, -2],
            ], dtype=float)

    Nreactions, Ncomponents = alfa.shape


    def f(y,t):
        ''' y is a vector with the concentration
		t is the indpendt variable (time)
		'''

        C = y[:] # mol*L-1
        
		# If you don't understand this
		# you can uncomment the following line and
		# run this part of the code line by line
		
		#import pdb; pdb.set_trace()
        # k=ko*exp(-Ea/(R*T));
		
		# Kinetic parameters

        kp1 = 99627.28  # L.molCat-1.s-1
        kp2 = 1884.84  # L.molCat-1.s-1
        kp3 = 915.49  # L.molCat-1.s-1
        kp4 = 479.29  # L.molCat-1.s-1
        K_A = 17.50  # L.mol-1
        K_B = 1.97   # L.mol-1
        K_H = 0.04  # non-dimensional
        Kest = 9.13  # L.mol-1

        #(1xR) en mol.molCat-1.s-1
        r1 = kp1*C[0] / (1+K_A*C[0]+np.sqrt(K_H)+K_B*C[1]+Kest*(C[2]+C[3]))**2          
        r2 = kp2*C[1] / (1+K_A*C[0]+np.sqrt(K_H)+K_B*C[1]+Kest*(C[2]+C[3]))**2 
        r3 = kp3*C[0] / (1+K_A*C[0]+np.sqrt(K_H)+K_B*C[1]+Kest*(C[2]+C[3]))**2
        r4 = kp4*C[0] / (1+K_A*C[0]+np.sqrt(K_H)+K_B*C[1]+Kest*(C[2]+C[3]))**2


        # D.Eqs

        dCAdt = nCat/VL*(-r1-r3-2*r4)
        dCBdt = nCat/VL*(r1-r2)
        dCCdt = nCat/VL*(r2+r3)
        dCDdt = nCat/VL*r4
        dCHdt = 0


        return [dCAdt, dCBdt, dCCdt, dCDdt, dCHdt]


    Ci0 =  [C_A0, C_B0, C_C0, C_D0, C_H0]  # mol*L-1
    t  = np.linspace(0, 400.*60, 100000)   # time grid

    # solving the DEs
    soln = odeint(f, Ci0, t)

    C_A0 = soln[:, 0] # mol*L-1
    C_B0 = soln[:, 1] # mol*L-1
    C_C0 = soln[:, 2] # mol*L-1
    C_D0   = soln[:, 3] # mol*L-1
    C_H20  = soln[:, 4] # mol*L-1

    t = t/60 # conversion to min

    
    # Plotting the results

    fig, ax1 = plt.subplots(figsize=(8,6))

    ax1.plot(t, C_A0, t, C_B0, t, C_C0)
    ax1.set_xlabel('time (min)')

    ax1.set_ylabel('$C_{i}$ (mol/L)')

    ax1.legend(['$C_A$','$C_B$','$C_C$'])

    plt.title('Heterogeneous STR: A -> B -> C')
    plt.suptitle('(CAChemE.org)')
    
    ax1.set_xlim([0,400])
    ax1.set_ylim([0,12])
    
	#ax2 = ax1.twinx()
    #ax2.plot(t, C_D0, color='y', ylim=(0,0.20) )
    #ax2.set_ylabel('Concentration [mol/L] of D')
    
    plt.show()


interact(STreactor,
         C_A0=(0.1,20,0.1), C_B0=(0.1,20,0.1),
         C_C0=(0.1,20,0.1), cat_load=(0.01,2.01,0.01) )