%matplotlib inline 

import numpy as np
import matplotlib.pyplot as plt 
from IPython.html.widgets import interact
from scipy.integrate import odeint

plt.style.use('bmh')

class Reactor(object):
    
    def __init__(self, h, Ta, adiabatico = False):
        
        #condiciones iniciales y constantes
        self.h = h # [m]
        
        self.T0 = 323. # [K] Tª inicial
        self.k0 = 3e3 # [min^-1]
        self.E = 30000. # [J/mol]
        self.R = 8.314 # [J·mol^-1·K^-1]
        self.Ca0 = 0.015 # [mol] inicial
        self.DH = -25000. # [J·mol^-1]
        
        if adiabatico == True:
            self.U = 0 # [J·(cm·K·min)^-1]
        else:
            self.U = 2.4 # [J·(cm·K·min)^-1]
        
        self.Ta = Ta # Cooling temperature [K]
        self.Cp = 125 # [J·mol^-1·K^-1]
     
    # La altura del reactor = al diametro interno
    # Volumen y área de reacción se pueden simplificar:
    
    def area(self):
        return (np.pi/4.*self.h**2) + (np.pi*self.h**2)
    
    def volumen(self):
        return np.pi/4.*self.h**3
   
    # Cinética
    def reac_veloc(self, T, na):
        
        # Arrhenius
        self.k = self.k0*np.exp(-1.*self.E/(self.R*T))
        
        # Orden 1
        self.r_a = -1.*self.k*self.concent(na)
        
        return self.r_a
    
    # Concentración moles/volumen
    def concent(self, na):
        V = self.volumen()
        return na/V
    
    # Balance de energía
    def dTdt(self, T, na):
        a = self.DH*self.reac_veloc(T,na)*self.volumen()
        b = self.U*self.area()*(self.Ta-T)
        c = self.Ca0*self.volumen()*self.Cp
        return (a+b)/c
    
    # Balance de matería
    def dNadt(self, T, na):
        return self.reac_veloc(T,na)*self.volumen()
    
    # Sistema de ecuaciones diferenciales
    def ode(self, y, t):
        
        # Condiciones iniciales
        T = y[0]
        na = y[1]
        
        # Ecuaciones diferenciales
        dTdt = self.dTdt(T,na)
        dnadt = self.dNadt(T,na)
        
        return [dTdt,dnadt]
        
        

def batch_reactor_plot(height, T_cooling, adiabatic):
    
    # Llamamos a la clase
    reac = Reactor(height,T_cooling, adiabatic)
    
    # Vector con espacio temporal a simular
    t = np.linspace(0,60,200)
    
    # Condiciones inciales
    y0 = [reac.T0, reac.Ca0*reac.volumen()]
    
    # Integración
    y_t = odeint(reac.ode, y0, t)
    
    
    # Representación de resultados    
    plt.plot(t,y_t[:,0])   
    plt.ylabel('Temperature / K')
    plt.xlabel('time / min')
    
    plt.title('Scale-up of a Bacth Reactor: \n'+
              'Reactor inner temperature during time')
    
    plt.ylim([300, 600])

    
    area_reactor =  reac.area()
    volume_reactor =  reac.volumen()
    
    area_to_volume = area_reactor/volume_reactor
    
    #Use 8 characters, give f decimal places
    a_reac_string = "%8.0f" % area_reactor
    volume_reactor_string = "%8.0f" % volume_reactor
    area2volume_string = "%8.5f" % area_to_volume
    
    # Adding these annotations to the graph
    
    plt.text(65, 323, #position (time, temperature)
             
             'Area to volume:' + area2volume_string + \
             ' $ \mathrm{cm^{-1}} $ \n'+ \
             
             'Surface area:' + a_reac_string + \
             ' $ \mathrm{cm^2} $ \n'+ \
             
             'Reacting volume:' + volume_reactor_string + \
             ' $ \mathrm{cm^3} $ \n \n'+ \
             
             '- CAChemE.org')
    
    #print(reac.h)
    
    

interact(batch_reactor_plot,
         height = (10, 200, 0.5),
         T_cooling = (300, 350, 0.5),
         adiabatic = False)