#!/usr/bin/env python
# coding: utf-8

# # Tutorial 3 - Part 2
# _applying our knowledge to solve Chem Reactor problems_
# 
# We now know ODEs
# 
# It is time everything together with other cool python modules to solve engineering problems.
# 
# If you don't know how to solve these problems it's ok, start small and make your way up.
# 
# The easiest example we could do would be to analyze the chemical reaction kinetics in the simplest case:
# 
# # Problem 1:
# 
# *A -> B *
# 
# We assume it is a first order reaction, and we use the design equation from the Batch reactor:
# 
# So that means that $-r_A = k C_A$.
# 
# Let's assume $k = 10^{-3} 1/s$ and $C_A(t=0) = 10$ moles
# 
# We can expect the half-time of A to be around 10000 seconds, which is ~3 hours.
# 
# Ok, let's see that through Python.
# 
# We want to track the concentration, which we will call $y_A$.
# 
# $dy_A/dt = - k C_A = -k y_A$
# 
# An important tool we need to use here is the ODEINT module:
# 
# `from scipy.integrate import odeint`
# 
# In this module if you have an equation $\frac{dy}{dt} = f(y,t)$ you need:
# 
# - To define at least the Python function f which implements the mathematical function $f$. It must take a vector y and a time t, and return the new vector f(y,t).
# 
# 
# Some additional commands that you may need in your script:
# 
# ```python
# #lets import the modules # plotting everything inline
# # %matplotlib inline 
# import matplotlib.pyplot as plt # plotting modules
# plt.style.use('presentation') # just have in your script for prettier plotting
# import numpy as np # our matlab-like module
# from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
# ```
# 
# To plot your data use these commands:
# 
# ```python
# plt.plot(times, y_A_exact,'k--', label='analytical')
# plt.xlabel('$t [s]$')
# plt.ylabel('$C_A [moles]$')
# plt.legend()
# plt.savefig('odeint.pdf')
# plt.show()
# 
# ```
# 

# In[ ]:


#let's import the modules # plotting everything inline
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt # plotting modules
#plt.style.use('presentation') # just have in your script for prettier plotting
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself

#reaction order
def f(y, t):
    return -k*y

k = 0.001 # L/(mol*sec)
time_start = 0 #s
time_finish = 10000 #s
N_points = 50

times = np.linspace(time_start,time_finish,N_points)
y_A0 = 10. # moles

y_A_calc = odeint(f, y_A0, times)
plt.plot(times, y_A_calc, 'ro-', label='integrated')


y_A_exact = y_A0*np.exp(-k*times)
plt.plot(times, y_A_exact,'k--', label='analytical')


plt.xlabel('$t [s]$')
plt.ylabel('$C_A [moles]$')
plt.legend()
plt.savefig('odeint.pdf')
plt.show()


# Now let's do the same calculation as case number 1 (of this problem), but this time we will use the conversion factor $X_A$ as our variable:
# 
# $X_A = \dfrac{N_{A0} - N_A}{N_{A0}} = \dfrac{C_{A0} - C_A}{C_{A0}}$ if the volume $V=const$
# 
# So the design equation looks like:
# 
# $\dfrac{d X_A}{dt} = (-r_A/N_{A0}) V $
# 

# In[ ]:


get_ipython().run_line_magic('reset', '')

get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt # plotting modules
plt.style.use('ggplot') # just have in your script for prettier plotting
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself

plt.clf()
k=0.001
y_A0 = 10.
time_start = 0 #s
time_finish = 10000 #s
N_points = 50

times = np.linspace(time_start,time_finish,N_points)

y_A_exact = y_A0*np.exp(-k*times)

x_A = ((y_A0 - y_A_exact)/y_A0)
# print(y_A_exact)
# print(x_A)
plt.plot(times, x_A, 'ro-')

plt.xlabel('$t [s]$')
plt.ylabel('$X_A$')
# plt.savefig('conversion factor.pdf')
plt.show()


# ## Problem 2:
# 
# *2A ->B*
# 
# We assume that it is a first order reaction and we use the design equation from the Batch reactor:
# 
# So that means that $-r_A = k [C_A]^2$.
# 
# Let's assume $k = 10^{-3} 1/s$ and $C_A(t=0) = 10$ moles
# 
# So we can expect the half-time of A to be around 10000 seconds which is ~3 hours.
# 
# Ok let's see that through Python.
# 
# We want to track the concentration which we will call $y_A$.
# 
# $dy_A/dt = - k C_A = -k y^2_A$
# 
# The exact solution here would be
# 
# $ 1/y_A = kt + C$
# $y_A(0) = 10$ => $C = 1/10$
# 
# 
# $1/y_A = kt + 1/y_A_0$
# 
# So now it is easier to track $1/y_A$ rather than the function itself
# 

# In[ ]:


# %matplotlib inline 
get_ipython().run_line_magic('reset', '')
import matplotlib.pyplot as plt # plotting modules
#plt.style.use('presentation') # just have in your script for prettier plotting
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself

#reaction order
def f(y, t):
    return -k*y**2

k = 0.001 # 1/sec
time_start = 0 #s
time_finish = 10 #s
N_points = 50

times = np.linspace(time_start,time_finish,N_points)
y_A0 = 10. # moles

y_A_calc = odeint(f, y_A0, times)
plt.plot(times, 1./y_A_calc, 'ro-', label='integrated')


inverse_y_A_exact = (1/y_A0 + k*times)
plt.plot(times, inverse_y_A_exact,'k--', label='analytical')


plt.xlabel('$t [s]$')
plt.ylabel('$1/C_A [moles]$')
plt.legend()
# plt.savefig('odeint.pdf')
plt.show()


# ## Problem 3: CSTR 
# 
# Having had solved ODEs, now we are going to solve "Simple" algebraic equations.
# 
# Given a continuously stirred tank reactor with a volume of $66 m^3$ where the reaction $A \to B$ occurs, at a rate of $−r_A= k C^2_A$ ($k=3 L/mol/h$), with an entering molar flow of $F_{A0} = 5 mol/h$ and a volumetric flowrate of $\upsilon = 10 L/h$, what is the exit concentration of A?
# 
# From a mole balance we know that at steady state $0=F_{A0}−F_A+ V \cdot r_A$. That equation simply states the sum of the molar flow of A in in minus the molar flow of A out plus the molar rate A is generated is equal to zero at steady state. This is directly the equation we need to solve. We need the following relationship:

# In[ ]:


#Copyright (C) 2013 by John Kitchin

from scipy.optimize import fsolve

Fa0 = 5.0
v0 = 10.

V = 66000.0  # reactor volume L^3
k = 3.0      # rate constant L/mol/h

def func(Ca):
    "Mole balance for a CSTR. Solve this equation for func(Ca)=0"
    Fa = v0 * Ca     # exit molar flow of A
    ra = -k * Ca**2  # rate of reaction of A L/mol/h
    return Fa0 - Fa + V * ra

# CA guess that 90 % is reacted away
CA_guess = 0.1 * Fa0 / v0
CA_sol, = fsolve(func, CA_guess)

print('The exit concentration is {0} mol/L'.format(CA_sol))


# In[ ]:


# Now let's solve the same problem using Python's 
# built-in functionality through a Python package odeint

# Step 1:
# Import modules into Python

import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
plt.style.use('fivethirtyeight')
# integration of ODEs so you don't have 
# to write your finite difference yourself
from scipy.integrate import odeint 

# Step 2:
# define the required parameters and functions
def f(y, t, k):
    dydt = -k*t
    return dydt

# define initial values, arrays, intervals
t_0 = 0.
t_f = 10.
npoints = 10
t_array = np.linspace(t_0, t_f, npoints + 1)
y_0 = 1.
k = 2.
h = t_array[1] - t_array[0]
y_num = np.zeros_like(t_array)

# let's also define a numerical solution
t_array_analyt = np.linspace(t_0, t_f, 1000 + 1)
y_analyt = -k*t_array_analyt**2/2. + y_0

# here the syntax is
# module odeint(right hand side function,
# initial value of y,
# array for the time values
# additional parameters like k in a tuple (k,))
# if we needed two parameters like k1, k2
# then args=(k1, k2)

y_num = odeint(f, y_0, t_array, args=(k,))


plt.xlabel('$t, s$')
plt.ylabel('$y$')
plt.plot(t_array_analyt, y_analyt, 'k--', lw=3,label='analytical')
plt.plot(t_array, y_num, 'go',alpha=0.6, label='numerical')
plt.legend()


# # Convergence
# 
# Let's compare a set number of intervals vs. 1001 intervals from the analytical, at a specific time. 
# For the sake of simplicity, I chose the last point. And because typically, error increases as we increase the number of steps, since every step is dependant on the last.
# 
# 
#     When we have 5 intervals the difference in $y$ is ~4.13922407e-08
# 
#     When we have 10 intervals the difference in $y$ is ~2.98024077e-08
# 
#     When we have 20 intervals the difference in $y$ is ~7.4505806e-09
# 
#     When we have 30 intervals the difference in $y$ is ~3.31137073e-09
# 
#     When we have 100 intervals the difference in $y$ is ~2.98030045e-10
# 
# 
# In these cases, we have an overestimate, which gets more accurate as we increase the number of intervals.
# 
# Although it is harder to see visually, most approximations are fairly accurate.

# In[ ]:


#https://en.wikipedia.org/wiki/Euler_method
#https://rosettacode.org/wiki/Euler_method

import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
plt.style.use('fivethirtyeight')
from scipy.integrate import odeint 

def f(y, t, k):
    dydt = -k*t
    return dydt

t_0 = 0.
t_f = 10

npoints0 = 5
t_array0 = np.linspace(t_0, t_f, npoints0 + 1)

npoints1 = 10
t_array1 = np.linspace(t_0, t_f, npoints1 + 1)

npoints2 = 20
t_array2 = np.linspace(t_0, t_f, npoints2 + 1)

npoints3 = 30
t_array3 = np.linspace(t_0, t_f, npoints3 + 1)

npoints4 = 100
t_array4 = np.linspace(t_0, t_f, npoints4 + 1)

y_0 = 1.
k = 2.
h = t_array1[1] - t_array1[0]

t_array_analyt = np.linspace(t_0, t_f, 1000 + 1)
y_analyt = -k*t_array_analyt**2/2. + y_0

y_num0 = np.zeros_like(t_array0)
y_num0 = odeint(f, y_0, t_array0, args=(k,))

y_num1 = np.zeros_like(t_array1)
y_num1 = odeint(f, y_0, t_array1, args=(k,))

y_num2 = np.zeros_like(t_array2)
y_num2 = odeint(f, y_0, t_array2, args=(k,))

y_num3 = np.zeros_like(t_array3)
y_num3 = odeint(f, y_0, t_array3, args=(k,))

y_num4 = np.zeros_like(t_array4)
y_num4 = odeint(f, y_0, t_array4, args=(k,))


plt.xlabel('$t, s$')
plt.ylabel('$y$')
plt.plot(t_array_analyt, y_analyt, 'k--', lw=3,label='analytical')

plt.plot(t_array0, y_num0, 'go',alpha=0.4, label='numerical0')
plt.plot(t_array1, y_num1, 'bo',alpha=0.4, label='numerical1')
plt.plot(t_array2, y_num2, 'co',alpha=0.4, label='numerical2')
plt.plot(t_array3, y_num3, 'mo',alpha=0.4, label='numerical3')
plt.plot(t_array4, y_num4, 'yo',alpha=0.4, label='numerical3')

plt.legend()

print(y_num0[5] - y_analyt[1000])
print(y_num1[10] - y_analyt[1000])
print(y_num2[20] - y_analyt[1000])
print(y_num3[30] - y_analyt[1000])
print(y_num4[100] - y_analyt[1000])


# # Problem 4
# 
# This problem should be familiar, but we will throw in a twist.
# 
# When $A -> B$ as simple reaction with the rate $-r_A = k C_A$ and the solution for A could be found in this way:
# $\dfrac{d C_A}{dt} = -k C_A(t)$
# $C_{A0} = 10$ mole/L

# In[ ]:


import matplotlib.pyplot as plt # plotting modules
get_ipython().run_line_magic('matplotlib', 'inline')
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself

#step 1 definition of the function:

def f(C_A, t):
    dcadt = - k*C_A
    return dcadt

k = 1
C_A0 = 10.    
time_start = 0.
time_finish = 10000. 
N_points = 100
time_array = np.linspace(0, 10., N_points)

C_num = odeint(f,C_A0, time_array)

plt.plot(time_array, C_num, 'r-')
plt.plot(time_array, C_num, 'go', alpha=0.8, fillstyle='none') # alpha - 80% transparency
# fillstyle - none - means no filling of the circles

plt.xlabel('t, s')
plt.ylabel('$C_A(t)$, mole/L')
plt.show()


# 
# Now the same thing for B
# $ -(r_B) = -(-r_A) = - k C_A$
# so the design equation will have an additional equation:
# 
# $\dfrac{d C_B}{dt} = k C_A(t)$
# $C_{B0} = 0$ mol/L
# 

# In[ ]:


import matplotlib.pyplot as plt # plotting modules
get_ipython().run_line_magic('matplotlib', 'inline')
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself

#step 1 definition of the function:

def f(C_list, t):
    C_A = C_list[0]
    C_B = C_list[1]
    dcadt = - k*C_A
    dcbdt = k*C_A
    return [dcadt, dcbdt]

k = 2
C_A0 = 10.    
C_B0 = 0.
time_start = 0.
time_finish = 10. 
N_points = 100
time_array = np.linspace(0, 10., N_points)

C0_list = [1, 0.]
C_num_list = odeint(f, C0_list, time_array)

print(C_num_list.shape)

C_A_num = C_num_list[:,0]
C_B_num = C_num_list[:,1]

plt.plot(time_array, C_A_num, 'r--', fillstyle='none', label='$A$')
plt.plot(time_array, C_B_num, 'g--', fillstyle='none', label='$B$') 

plt.xlabel('t, s')
plt.ylabel('$C(t)$, mole/L')
plt.legend()
plt.show()


# To see precisely what the values were from the integration, we can use 
# 
# ```python
# print(C_num_list[:,0])
# ```

# In[ ]:


C_num_list[:,0]


# # Problem 5
# This one is more complicated:
# 
# 
# A system of ODES
# 
# $ A + B \to C$, $k_1 = 1$ L/(mol s)
# 
# $ B + C \to D$, $k_1 = 1.5$ L/(mol s)
# 
# with initial concentrations 
# $C_{A0} = 1$ mol/L
# $C_{B0} = 1$ mol/L
# $C_{C0} = 0$ mol/L
# $C_{D0} = 0$ mol/L
# 
# 1. Calculate concentration of each component
# 2. Calculate the selectivity parameters $S = \dfrac{C}{(C+D)}$, S(t=0) = 1.0
# 
# r1 = k1*C_A*C_B
# r2 = k2*C_B*C_C
# 
# dcAdt = -r1
# dcBdt = -r1 - r2
# dcCdt = r1 - r2
# dcDdt = r2
#     
#     

# In[ ]:


import matplotlib.pyplot as plt # plotting modules
get_ipython().run_line_magic('matplotlib', 'inline')
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself

#step 1 definition of the function:

def f(C_list, t):

    k1 = 1.
    k2 = 1.5
    C_A = C_list[0]
    C_B = C_list[1]
    C_C = C_list[2]
    
    r1 = k1*C_A*C_B
    r2 = k2*C_B*C_C

    dcAdt = -r1
    dcBdt = -r1 - r2
    dcCdt = r1 - r2
    dcDdt = r2
    
    return [dcAdt,dcBdt,dcCdt,dcDdt]

k1 = 1.
k2 = 1.5

C0_list = [1.0, 1.0, 0.0, 0.0]

time_start = 0.
time_finish = 3. 
N_points = 100
time_array = np.linspace(0, 3., N_points)


C_num_list = odeint(f, C0_list, time_array)

print(C_num_list.shape)

C_A_num = C_num_list[:,0]
C_B_num = C_num_list[:,1]
C_C_num = C_num_list[:,2]
C_D_num = C_num_list[:,3]

plt.plot(time_array, C_A_num, 'r-', label='$A$')
plt.plot(time_array, C_B_num, 'g-', label='$B$') 
plt.plot(time_array, C_C_num, 'm--', label='$C$') 
plt.plot(time_array, C_D_num, 'k--', label='$D$') 

plt.xlabel('t, s')
plt.ylabel('$C(t)$, mole/L')
plt.legend()
plt.show()


# 
# # Additional problem : 
# 
# $$\require{mhchem}$$
# Normal butane, C<sub>4</sub>H<sub>10</sub>, is to be isomerized to isobutane in a plug-flow reactor, because isobutane is worth more. This elementary reversible reaction is to be carried out adiabatically in the liquid phase under high pressure using trace amounts of a liquid catalyst which gives a specific reaction rate of 31.1/hr at 360K. The feed enters at 330K. 
# 
# Additional information: 
# * $\Delta H^\circ_{Rx} = -6900$ J/mol $n$-butane
# * Activation energy $E_a = 65.7$ kJ/mol 
# * $K_C = 3.03$ at $60^\circ$ C  
# * $C_{P,n\text{-Butane}} = 141$ J/mol/K 
# * $C_{P,i\text{-Butane}} = 141$ J/mol/K 
# * $C_{P,n\text{-Pentane}} = 161$ J/mol/K
# * $C_{A0} = 9.3$ mol/dm3 = 9.3 kmol/m3 feed concentration of n-butane
# 
# a) Calculate the PFR volume necessary to process 100,000 gal/day (163 kmol/h) at 70% conversion of a mixture of 90 mol% n-butane and 10 mol% i-pentane, which is considered an inert.
# 
# b) Plot and analyze $X, X_e, T$, and $-r_A$ down the length of the reactor.
# 
# 
# You may use this code in your program.
# You just need to understand and plot it.
# 
# ```python
# DHrx = -6900. # J/mol
# Ea = 65700. # J/mol
# R = 8.314 # J/mol/K
# T0 = 330 # K
# FA0 = 0.9 * 163e3 # mol/hr
# CA0 = 9.3e3 # mol/m3
# 
# def Kc(T):
#     "Equilibrium constant, as a function of temperature"
#     return 3.03 * np.exp(( -1* DHrx / R) * (1./T - 1./(273.15+60)) )
# assert Kc(273.15+60) == 3.03
# assert Kc(273.15+65) < 3.03 # le Chatellier's principle
# 
# def k(T):
#     "Specific reaction rate in 1 hour, as a function of temperature"
#     return 31.1 * np.exp((-1*Ea/R) * (1./T - 1./360.))
# assert k(360) == 31.1
# assert k(365) > 31.1
# 
# def Teb(X):
#     "Temperature from energy balance, as a function of conversion."
#     Cps = np.array([141, 141, 161])
#     Thetas = np.array([1., 0., 0.1/0.9])
#     return T0 - DHrx * X / sum(Cps*Thetas)
# assert Teb(0) == T0
# assert Teb(0.5) > T0 # Exothermic
# 
# def dXdV(X,V):
#     """
#     dX/dV in a plug flow reactor is -rA/ FA0
#     """
#     T = Teb(X)
#     CA = CA0*(1.-X)
#     CB = CA0*X
#     rA = -k(T)*(CA - CB/Kc(T))
#     return -rA/FA0
# ```
# 
# ```
# Xeq = Kc(Teb(X))/ (1+Kc(Teb(X)))
# rate = k(Teb(X))*(CA0*(1.-X) - CA0*X/Kc(Teb(X)))
# 
# #For Levenspiel plot you may use this:
# 
# plt.plot(X, FA0/rate)
# plt.ylim(0,10)
# plt.title("Levenspiel Plot")
# plt.xlabel("X")
# plt.ylabel("$\\frac{F_{A0}}{-r_A}$")
# plt.plot((Xeq[-1],Xeq[-1]),(0,FA0/rate[-1]), ':')
# plt.show()
# ```

# In[ ]:


import matplotlib.pyplot as plt # plotting modules
get_ipython().run_line_magic('matplotlib', 'inline')
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
import scipy
from scipy.integrate import odeint

DHrx = -6900. # J/mol
Ea = 65700. # J/mol
R = 8.314 # J/mol/K
T0 = 330 # K
FA0 = 0.9 * 163e3 # mol/hr
CA0 = 9.3e3 # mol/m3

def Kc(T):
    "Equilibrium constant, as a function of temperature"
    return 3.03 * np.exp(( -1* DHrx / R) * (1./T - 1./(273.15+60)) )
assert Kc(273.15+60) == 3.03
assert Kc(273.15+65) < 3.03 # le Chatellier's principle

def k(T):
    "Specific reaction rate in 1/hour, as a function of temperature"
    return 31.1 * np.exp((-1*Ea/R) * (1./T - 1./360.))
assert k(360) == 31.1
assert k(365) > 31.1

def Teb(X):
    "Temperature from energy balance, as a function of conversion."
    Cps = np.array([141, 141, 161])
    Thetas = np.array([1., 0., 0.1/0.9])
    return T0 - DHrx * X / sum(Cps*Thetas)
assert Teb(0) == T0
assert Teb(0.5) > T0 # Exothermic

def dXdV(X,V):
    """
    dX/dV in a plug flow reactor is -rA/ FA0
    """
    T = Teb(X)
    CA = CA0*(1.-X)
    CB = CA0*X
    rA = -k(T)*(CA - CB/Kc(T))
    return -rA/FA0
    
volumes = np.linspace(0,5,501)
X = scipy.integrate.odeint(dXdV, 0., volumes)
plt.plot(volumes,X,label="Conversion")

Xeq = Kc(Teb(X))/ (1+Kc(Teb(X)))
plt.plot(volumes, Xeq,label="Equilibrium conversion")
plt.legend(loc="best")
plt.xlabel("Volume")
plt.ylabel("Conversion X")
plt.show()

rate = k(Teb(X))*(CA0*(1.-X) - CA0*X/Kc(Teb(X)))
plt.plot(volumes, rate, label="Rate")
plt.xlabel("Volume")
plt.ylabel("Rate")
plt.show()

for volume, conversion in zip(volumes, X):
    if conversion>0.7:
        break
print(volume)
print(conversion)


plt.plot(X, FA0/rate)
plt.ylim(0,10)
plt.title("Levenspiel Plot")
plt.xlabel("X")
plt.ylabel("$\\frac{F_{A0}}{-r_A}$")
plt.plot((Xeq[-1],Xeq[-1]),(0,FA0/rate[-1]), ':')
plt.show()