{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "bc68411114a596781fbbdee2115f17dd",
     "grade": false,
     "grade_id": "cell-3bbf4e304107dc5a",
     "locked": true,
     "schema_version": 1,
     "solution": false
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Tutorial 3 - Part 2\n",
    "_applying our knowledge to solve Chem Reactor problems_\n",
    "\n",
    "We now know ODEs\n",
    "\n",
    "It is time everything together with other cool python modules to solve engineering problems.\n",
    "\n",
    "If you don't know how to solve these problems it's ok, start small and make your way up.\n",
    "\n",
    "The easiest example we could do would be to analyze the chemical reaction kinetics in the simplest case:\n",
    "\n",
    "# Problem 1:\n",
    "\n",
    "*A -> B *\n",
    "\n",
    "We assume it is a first order reaction, and we use the design equation from the Batch reactor:\n",
    "\n",
    "So that means that $-r_A = k C_A$.\n",
    "\n",
    "Let's assume $k = 10^{-3} 1/s$ and $C_A(t=0) = 10$ moles\n",
    "\n",
    "We can expect the half-time of A to be around 10000 seconds, which is ~3 hours.\n",
    "\n",
    "Ok, let's see that through Python.\n",
    "\n",
    "We want to track the concentration, which we will call $y_A$.\n",
    "\n",
    "$dy_A/dt = - k C_A = -k y_A$\n",
    "\n",
    "An important tool we need to use here is the ODEINT module:\n",
    "\n",
    "`from scipy.integrate import odeint`\n",
    "\n",
    "In this module if you have an equation $\\frac{dy}{dt} = f(y,t)$ you need:\n",
    "\n",
    "- 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).\n",
    "\n",
    "\n",
    "Some additional commands that you may need in your script:\n",
    "\n",
    "```python\n",
    "#lets import the modules # plotting everything inline\n",
    "# %matplotlib inline \n",
    "import matplotlib.pyplot as plt # plotting modules\n",
    "plt.style.use('presentation') # just have in your script for prettier plotting\n",
    "import numpy as np # our matlab-like module\n",
    "from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself\n",
    "```\n",
    "\n",
    "To plot your data use these commands:\n",
    "\n",
    "```python\n",
    "plt.plot(times, y_A_exact,'k--', label='analytical')\n",
    "plt.xlabel('$t [s]$')\n",
    "plt.ylabel('$C_A [moles]$')\n",
    "plt.legend()\n",
    "plt.savefig('odeint.pdf')\n",
    "plt.show()\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#let's import the modules # plotting everything inline\n",
    "%matplotlib inline \n",
    "import matplotlib.pyplot as plt # plotting modules\n",
    "#plt.style.use('presentation') # just have in your script for prettier plotting\n",
    "import numpy as np # our matlab-like module\n",
    "from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself\n",
    "\n",
    "#reaction order\n",
    "def f(y, t):\n",
    "    return -k*y\n",
    "\n",
    "k = 0.001 # L/(mol*sec)\n",
    "time_start = 0 #s\n",
    "time_finish = 10000 #s\n",
    "N_points = 50\n",
    "\n",
    "times = np.linspace(time_start,time_finish,N_points)\n",
    "y_A0 = 10. # moles\n",
    "\n",
    "y_A_calc = odeint(f, y_A0, times)\n",
    "plt.plot(times, y_A_calc, 'ro-', label='integrated')\n",
    "\n",
    "\n",
    "y_A_exact = y_A0*np.exp(-k*times)\n",
    "plt.plot(times, y_A_exact,'k--', label='analytical')\n",
    "\n",
    "\n",
    "plt.xlabel('$t [s]$')\n",
    "plt.ylabel('$C_A [moles]$')\n",
    "plt.legend()\n",
    "plt.savefig('odeint.pdf')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "182158bfcdc905e6e2da8cb50ca92095",
     "grade": false,
     "grade_id": "cell-78c964bc077f121e",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "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:\n",
    "\n",
    "$X_A = \\dfrac{N_{A0} - N_A}{N_{A0}} = \\dfrac{C_{A0} - C_A}{C_{A0}}$ if the volume $V=const$\n",
    "\n",
    "So the design equation looks like:\n",
    "\n",
    "$\\dfrac{d X_A}{dt} = (-r_A/N_{A0}) V $\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reset\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt # plotting modules\n",
    "plt.style.use('ggplot') # just have in your script for prettier plotting\n",
    "import numpy as np # our matlab-like module\n",
    "from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself\n",
    "\n",
    "plt.clf()\n",
    "k=0.001\n",
    "y_A0 = 10.\n",
    "time_start = 0 #s\n",
    "time_finish = 10000 #s\n",
    "N_points = 50\n",
    "\n",
    "times = np.linspace(time_start,time_finish,N_points)\n",
    "\n",
    "y_A_exact = y_A0*np.exp(-k*times)\n",
    "\n",
    "x_A = ((y_A0 - y_A_exact)/y_A0)\n",
    "# print(y_A_exact)\n",
    "# print(x_A)\n",
    "plt.plot(times, x_A, 'ro-')\n",
    "\n",
    "plt.xlabel('$t [s]$')\n",
    "plt.ylabel('$X_A$')\n",
    "# plt.savefig('conversion factor.pdf')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "5ee0765159009802eff9ee396f3d1771",
     "grade": false,
     "grade_id": "cell-63d325154a978c19",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Problem 2:\n",
    "\n",
    "*2A ->B*\n",
    "\n",
    "We assume that it is a first order reaction and we use the design equation from the Batch reactor:\n",
    "\n",
    "So that means that $-r_A = k [C_A]^2$.\n",
    "\n",
    "Let's assume $k = 10^{-3} 1/s$ and $C_A(t=0) = 10$ moles\n",
    "\n",
    "So we can expect the half-time of A to be around 10000 seconds which is ~3 hours.\n",
    "\n",
    "Ok let's see that through Python.\n",
    "\n",
    "We want to track the concentration which we will call $y_A$.\n",
    "\n",
    "$dy_A/dt = - k C_A = -k y^2_A$\n",
    "\n",
    "The exact solution here would be\n",
    "\n",
    "$ 1/y_A = kt + C$\n",
    "$y_A(0) = 10$ => $C = 1/10$\n",
    "\n",
    "\n",
    "$1/y_A = kt + 1/y_A_0$\n",
    "\n",
    "So now it is easier to track $1/y_A$ rather than the function itself\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %matplotlib inline \n",
    "%reset\n",
    "import matplotlib.pyplot as plt # plotting modules\n",
    "#plt.style.use('presentation') # just have in your script for prettier plotting\n",
    "import numpy as np # our matlab-like module\n",
    "from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself\n",
    "\n",
    "#reaction order\n",
    "def f(y, t):\n",
    "    return -k*y**2\n",
    "\n",
    "k = 0.001 # 1/sec\n",
    "time_start = 0 #s\n",
    "time_finish = 10 #s\n",
    "N_points = 50\n",
    "\n",
    "times = np.linspace(time_start,time_finish,N_points)\n",
    "y_A0 = 10. # moles\n",
    "\n",
    "y_A_calc = odeint(f, y_A0, times)\n",
    "plt.plot(times, 1./y_A_calc, 'ro-', label='integrated')\n",
    "\n",
    "\n",
    "inverse_y_A_exact = (1/y_A0 + k*times)\n",
    "plt.plot(times, inverse_y_A_exact,'k--', label='analytical')\n",
    "\n",
    "\n",
    "plt.xlabel('$t [s]$')\n",
    "plt.ylabel('$1/C_A [moles]$')\n",
    "plt.legend()\n",
    "# plt.savefig('odeint.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "2d61628bcd471a9bbe8a5ae272d8934c",
     "grade": false,
     "grade_id": "cell-a62b63367a825b3d",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Problem 3: CSTR \n",
    "\n",
    "Having had solved ODEs, now we are going to solve \"Simple\" algebraic equations.\n",
    "\n",
    "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?\n",
    "\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Copyright (C) 2013 by John Kitchin\n",
    "\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "Fa0 = 5.0\n",
    "v0 = 10.\n",
    "\n",
    "V = 66000.0  # reactor volume L^3\n",
    "k = 3.0      # rate constant L/mol/h\n",
    "\n",
    "def func(Ca):\n",
    "    \"Mole balance for a CSTR. Solve this equation for func(Ca)=0\"\n",
    "    Fa = v0 * Ca     # exit molar flow of A\n",
    "    ra = -k * Ca**2  # rate of reaction of A L/mol/h\n",
    "    return Fa0 - Fa + V * ra\n",
    "\n",
    "# CA guess that 90 % is reacted away\n",
    "CA_guess = 0.1 * Fa0 / v0\n",
    "CA_sol, = fsolve(func, CA_guess)\n",
    "\n",
    "print('The exit concentration is {0} mol/L'.format(CA_sol))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's solve the same problem using Python's \n",
    "# built-in functionality through a Python package odeint\n",
    "\n",
    "# Step 1:\n",
    "# Import modules into Python\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.style.use('fivethirtyeight')\n",
    "# integration of ODEs so you don't have \n",
    "# to write your finite difference yourself\n",
    "from scipy.integrate import odeint \n",
    "\n",
    "# Step 2:\n",
    "# define the required parameters and functions\n",
    "def f(y, t, k):\n",
    "    dydt = -k*t\n",
    "    return dydt\n",
    "\n",
    "# define initial values, arrays, intervals\n",
    "t_0 = 0.\n",
    "t_f = 10.\n",
    "npoints = 10\n",
    "t_array = np.linspace(t_0, t_f, npoints + 1)\n",
    "y_0 = 1.\n",
    "k = 2.\n",
    "h = t_array[1] - t_array[0]\n",
    "y_num = np.zeros_like(t_array)\n",
    "\n",
    "# let's also define a numerical solution\n",
    "t_array_analyt = np.linspace(t_0, t_f, 1000 + 1)\n",
    "y_analyt = -k*t_array_analyt**2/2. + y_0\n",
    "\n",
    "# here the syntax is\n",
    "# module odeint(right hand side function,\n",
    "# initial value of y,\n",
    "# array for the time values\n",
    "# additional parameters like k in a tuple (k,))\n",
    "# if we needed two parameters like k1, k2\n",
    "# then args=(k1, k2)\n",
    "\n",
    "y_num = odeint(f, y_0, t_array, args=(k,))\n",
    "\n",
    "\n",
    "plt.xlabel('$t, s$')\n",
    "plt.ylabel('$y$')\n",
    "plt.plot(t_array_analyt, y_analyt, 'k--', lw=3,label='analytical')\n",
    "plt.plot(t_array, y_num, 'go',alpha=0.6, label='numerical')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convergence\n",
    "\n",
    "Let's compare a set number of intervals vs. 1001 intervals from the analytical, at a specific time. \n",
    "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.\n",
    "\n",
    "\n",
    "    When we have 5 intervals the difference in $y$ is ~4.13922407e-08\n",
    "\n",
    "    When we have 10 intervals the difference in $y$ is ~2.98024077e-08\n",
    "\n",
    "    When we have 20 intervals the difference in $y$ is ~7.4505806e-09\n",
    "\n",
    "    When we have 30 intervals the difference in $y$ is ~3.31137073e-09\n",
    "\n",
    "    When we have 100 intervals the difference in $y$ is ~2.98030045e-10\n",
    "\n",
    "\n",
    "In these cases, we have an overestimate, which gets more accurate as we increase the number of intervals.\n",
    "\n",
    "Although it is harder to see visually, most approximations are fairly accurate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#https://en.wikipedia.org/wiki/Euler_method\n",
    "#https://rosettacode.org/wiki/Euler_method\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.style.use('fivethirtyeight')\n",
    "from scipy.integrate import odeint \n",
    "\n",
    "def f(y, t, k):\n",
    "    dydt = -k*t\n",
    "    return dydt\n",
    "\n",
    "t_0 = 0.\n",
    "t_f = 10\n",
    "\n",
    "npoints0 = 5\n",
    "t_array0 = np.linspace(t_0, t_f, npoints0 + 1)\n",
    "\n",
    "npoints1 = 10\n",
    "t_array1 = np.linspace(t_0, t_f, npoints1 + 1)\n",
    "\n",
    "npoints2 = 20\n",
    "t_array2 = np.linspace(t_0, t_f, npoints2 + 1)\n",
    "\n",
    "npoints3 = 30\n",
    "t_array3 = np.linspace(t_0, t_f, npoints3 + 1)\n",
    "\n",
    "npoints4 = 100\n",
    "t_array4 = np.linspace(t_0, t_f, npoints4 + 1)\n",
    "\n",
    "y_0 = 1.\n",
    "k = 2.\n",
    "h = t_array1[1] - t_array1[0]\n",
    "\n",
    "t_array_analyt = np.linspace(t_0, t_f, 1000 + 1)\n",
    "y_analyt = -k*t_array_analyt**2/2. + y_0\n",
    "\n",
    "y_num0 = np.zeros_like(t_array0)\n",
    "y_num0 = odeint(f, y_0, t_array0, args=(k,))\n",
    "\n",
    "y_num1 = np.zeros_like(t_array1)\n",
    "y_num1 = odeint(f, y_0, t_array1, args=(k,))\n",
    "\n",
    "y_num2 = np.zeros_like(t_array2)\n",
    "y_num2 = odeint(f, y_0, t_array2, args=(k,))\n",
    "\n",
    "y_num3 = np.zeros_like(t_array3)\n",
    "y_num3 = odeint(f, y_0, t_array3, args=(k,))\n",
    "\n",
    "y_num4 = np.zeros_like(t_array4)\n",
    "y_num4 = odeint(f, y_0, t_array4, args=(k,))\n",
    "\n",
    "\n",
    "plt.xlabel('$t, s$')\n",
    "plt.ylabel('$y$')\n",
    "plt.plot(t_array_analyt, y_analyt, 'k--', lw=3,label='analytical')\n",
    "\n",
    "plt.plot(t_array0, y_num0, 'go',alpha=0.4, label='numerical0')\n",
    "plt.plot(t_array1, y_num1, 'bo',alpha=0.4, label='numerical1')\n",
    "plt.plot(t_array2, y_num2, 'co',alpha=0.4, label='numerical2')\n",
    "plt.plot(t_array3, y_num3, 'mo',alpha=0.4, label='numerical3')\n",
    "plt.plot(t_array4, y_num4, 'yo',alpha=0.4, label='numerical3')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "print(y_num0[5] - y_analyt[1000])\n",
    "print(y_num1[10] - y_analyt[1000])\n",
    "print(y_num2[20] - y_analyt[1000])\n",
    "print(y_num3[30] - y_analyt[1000])\n",
    "print(y_num4[100] - y_analyt[1000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 4\n",
    "\n",
    "This problem should be familiar, but we will throw in a twist.\n",
    "\n",
    "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:\n",
    "$\\dfrac{d C_A}{dt} = -k C_A(t)$\n",
    "$C_{A0} = 10$ mole/L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt # plotting modules\n",
    "%matplotlib inline \n",
    "plt.style.use('ggplot') # just have in your script for prettier plotting\n",
    "# instead of 'ggplot' you can use 'presentation'\n",
    "import numpy as np # our matlab-like module\n",
    "from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself\n",
    "\n",
    "#step 1 definition of the function:\n",
    "\n",
    "def f(C_A, t):\n",
    "    dcadt = - k*C_A\n",
    "    return dcadt\n",
    "\n",
    "k = 1\n",
    "C_A0 = 10.    \n",
    "time_start = 0.\n",
    "time_finish = 10000. \n",
    "N_points = 100\n",
    "time_array = np.linspace(0, 10., N_points)\n",
    "\n",
    "C_num = odeint(f,C_A0, time_array)\n",
    "\n",
    "plt.plot(time_array, C_num, 'r-')\n",
    "plt.plot(time_array, C_num, 'go', alpha=0.8, fillstyle='none') # alpha - 80% transparency\n",
    "# fillstyle - none - means no filling of the circles\n",
    "\n",
    "plt.xlabel('t, s')\n",
    "plt.ylabel('$C_A(t)$, mole/L')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now the same thing for B\n",
    "$ -(r_B) = -(-r_A) = - k C_A$\n",
    "so the design equation will have an additional equation:\n",
    "\n",
    "$\\dfrac{d C_B}{dt} = k C_A(t)$\n",
    "$C_{B0} = 0$ mol/L\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt # plotting modules\n",
    "%matplotlib inline \n",
    "plt.style.use('ggplot') # just have in your script for prettier plotting\n",
    "# instead of 'ggplot' you can use 'presentation'\n",
    "import numpy as np # our matlab-like module\n",
    "from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself\n",
    "\n",
    "#step 1 definition of the function:\n",
    "\n",
    "def f(C_list, t):\n",
    "    C_A = C_list[0]\n",
    "    C_B = C_list[1]\n",
    "    dcadt = - k*C_A\n",
    "    dcbdt = k*C_A\n",
    "    return [dcadt, dcbdt]\n",
    "\n",
    "k = 2\n",
    "C_A0 = 10.    \n",
    "C_B0 = 0.\n",
    "time_start = 0.\n",
    "time_finish = 10. \n",
    "N_points = 100\n",
    "time_array = np.linspace(0, 10., N_points)\n",
    "\n",
    "C0_list = [1, 0.]\n",
    "C_num_list = odeint(f, C0_list, time_array)\n",
    "\n",
    "print(C_num_list.shape)\n",
    "\n",
    "C_A_num = C_num_list[:,0]\n",
    "C_B_num = C_num_list[:,1]\n",
    "\n",
    "plt.plot(time_array, C_A_num, 'r--', fillstyle='none', label='$A$')\n",
    "plt.plot(time_array, C_B_num, 'g--', fillstyle='none', label='$B$') \n",
    "\n",
    "plt.xlabel('t, s')\n",
    "plt.ylabel('$C(t)$, mole/L')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see precisely what the values were from the integration, we can use \n",
    "\n",
    "```python\n",
    "print(C_num_list[:,0])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_num_list[:,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "e46b61f0c86e2680150507d49706a432",
     "grade": false,
     "grade_id": "cell-90a20d720bb72257",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "# Problem 5\n",
    "This one is more complicated:\n",
    "\n",
    "\n",
    "A system of ODES\n",
    "\n",
    "$ A + B \\to C$, $k_1 = 1$ L/(mol s)\n",
    "\n",
    "$ B + C \\to D$, $k_1 = 1.5$ L/(mol s)\n",
    "\n",
    "with initial concentrations \n",
    "$C_{A0} = 1$ mol/L\n",
    "$C_{B0} = 1$ mol/L\n",
    "$C_{C0} = 0$ mol/L\n",
    "$C_{D0} = 0$ mol/L\n",
    "\n",
    "1. Calculate concentration of each component\n",
    "2. Calculate the selectivity parameters $S = \\dfrac{C}{(C+D)}$, S(t=0) = 1.0\n",
    "\n",
    "r1 = k1*C_A*C_B\n",
    "r2 = k2*C_B*C_C\n",
    "\n",
    "dcAdt = -r1\n",
    "dcBdt = -r1 - r2\n",
    "dcCdt = r1 - r2\n",
    "dcDdt = r2\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt # plotting modules\n",
    "%matplotlib inline \n",
    "plt.style.use('ggplot') # just have in your script for prettier plotting\n",
    "# instead of 'ggplot' you can use 'presentation'\n",
    "import numpy as np # our matlab-like module\n",
    "from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself\n",
    "\n",
    "#step 1 definition of the function:\n",
    "\n",
    "def f(C_list, t):\n",
    "\n",
    "    k1 = 1.\n",
    "    k2 = 1.5\n",
    "    C_A = C_list[0]\n",
    "    C_B = C_list[1]\n",
    "    C_C = C_list[2]\n",
    "    \n",
    "    r1 = k1*C_A*C_B\n",
    "    r2 = k2*C_B*C_C\n",
    "\n",
    "    dcAdt = -r1\n",
    "    dcBdt = -r1 - r2\n",
    "    dcCdt = r1 - r2\n",
    "    dcDdt = r2\n",
    "    \n",
    "    return [dcAdt,dcBdt,dcCdt,dcDdt]\n",
    "\n",
    "k1 = 1.\n",
    "k2 = 1.5\n",
    "\n",
    "C0_list = [1.0, 1.0, 0.0, 0.0]\n",
    "\n",
    "time_start = 0.\n",
    "time_finish = 3. \n",
    "N_points = 100\n",
    "time_array = np.linspace(0, 3., N_points)\n",
    "\n",
    "\n",
    "C_num_list = odeint(f, C0_list, time_array)\n",
    "\n",
    "print(C_num_list.shape)\n",
    "\n",
    "C_A_num = C_num_list[:,0]\n",
    "C_B_num = C_num_list[:,1]\n",
    "C_C_num = C_num_list[:,2]\n",
    "C_D_num = C_num_list[:,3]\n",
    "\n",
    "plt.plot(time_array, C_A_num, 'r-', label='$A$')\n",
    "plt.plot(time_array, C_B_num, 'g-', label='$B$') \n",
    "plt.plot(time_array, C_C_num, 'm--', label='$C$') \n",
    "plt.plot(time_array, C_D_num, 'k--', label='$D$') \n",
    "\n",
    "plt.xlabel('t, s')\n",
    "plt.ylabel('$C(t)$, mole/L')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Additional problem : \n",
    "\n",
    "$$\\require{mhchem}$$\n",
    "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. \n",
    "\n",
    "Additional information: \n",
    "* $\\Delta H^\\circ_{Rx} = -6900$ J/mol $n$-butane\n",
    "* Activation energy $E_a = 65.7$ kJ/mol \n",
    "* $K_C = 3.03$ at $60^\\circ$ C  \n",
    "* $C_{P,n\\text{-Butane}} = 141$ J/mol/K \n",
    "* $C_{P,i\\text{-Butane}} = 141$ J/mol/K \n",
    "* $C_{P,n\\text{-Pentane}} = 161$ J/mol/K\n",
    "* $C_{A0} = 9.3$ mol/dm3 = 9.3 kmol/m3 feed concentration of n-butane\n",
    "\n",
    "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.\n",
    "\n",
    "b) Plot and analyze $X, X_e, T$, and $-r_A$ down the length of the reactor.\n",
    "\n",
    "\n",
    "You may use this code in your program.\n",
    "You just need to understand and plot it.\n",
    "\n",
    "```python\n",
    "DHrx = -6900. # J/mol\n",
    "Ea = 65700. # J/mol\n",
    "R = 8.314 # J/mol/K\n",
    "T0 = 330 # K\n",
    "FA0 = 0.9 * 163e3 # mol/hr\n",
    "CA0 = 9.3e3 # mol/m3\n",
    "\n",
    "def Kc(T):\n",
    "    \"Equilibrium constant, as a function of temperature\"\n",
    "    return 3.03 * np.exp(( -1* DHrx / R) * (1./T - 1./(273.15+60)) )\n",
    "assert Kc(273.15+60) == 3.03\n",
    "assert Kc(273.15+65) < 3.03 # le Chatellier's principle\n",
    "\n",
    "def k(T):\n",
    "    \"Specific reaction rate in 1 hour, as a function of temperature\"\n",
    "    return 31.1 * np.exp((-1*Ea/R) * (1./T - 1./360.))\n",
    "assert k(360) == 31.1\n",
    "assert k(365) > 31.1\n",
    "\n",
    "def Teb(X):\n",
    "    \"Temperature from energy balance, as a function of conversion.\"\n",
    "    Cps = np.array([141, 141, 161])\n",
    "    Thetas = np.array([1., 0., 0.1/0.9])\n",
    "    return T0 - DHrx * X / sum(Cps*Thetas)\n",
    "assert Teb(0) == T0\n",
    "assert Teb(0.5) > T0 # Exothermic\n",
    "\n",
    "def dXdV(X,V):\n",
    "    \"\"\"\n",
    "    dX/dV in a plug flow reactor is -rA/ FA0\n",
    "    \"\"\"\n",
    "    T = Teb(X)\n",
    "    CA = CA0*(1.-X)\n",
    "    CB = CA0*X\n",
    "    rA = -k(T)*(CA - CB/Kc(T))\n",
    "    return -rA/FA0\n",
    "```\n",
    "\n",
    "```\n",
    "Xeq = Kc(Teb(X))/ (1+Kc(Teb(X)))\n",
    "rate = k(Teb(X))*(CA0*(1.-X) - CA0*X/Kc(Teb(X)))\n",
    "\n",
    "#For Levenspiel plot you may use this:\n",
    "\n",
    "plt.plot(X, FA0/rate)\n",
    "plt.ylim(0,10)\n",
    "plt.title(\"Levenspiel Plot\")\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"$\\\\frac{F_{A0}}{-r_A}$\")\n",
    "plt.plot((Xeq[-1],Xeq[-1]),(0,FA0/rate[-1]), ':')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt # plotting modules\n",
    "%matplotlib inline \n",
    "plt.style.use('ggplot') # just have in your script for prettier plotting\n",
    "# instead of 'ggplot' you can use 'presentation'\n",
    "import numpy as np # our matlab-like module\n",
    "import scipy\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "DHrx = -6900. # J/mol\n",
    "Ea = 65700. # J/mol\n",
    "R = 8.314 # J/mol/K\n",
    "T0 = 330 # K\n",
    "FA0 = 0.9 * 163e3 # mol/hr\n",
    "CA0 = 9.3e3 # mol/m3\n",
    "\n",
    "def Kc(T):\n",
    "    \"Equilibrium constant, as a function of temperature\"\n",
    "    return 3.03 * np.exp(( -1* DHrx / R) * (1./T - 1./(273.15+60)) )\n",
    "assert Kc(273.15+60) == 3.03\n",
    "assert Kc(273.15+65) < 3.03 # le Chatellier's principle\n",
    "\n",
    "def k(T):\n",
    "    \"Specific reaction rate in 1/hour, as a function of temperature\"\n",
    "    return 31.1 * np.exp((-1*Ea/R) * (1./T - 1./360.))\n",
    "assert k(360) == 31.1\n",
    "assert k(365) > 31.1\n",
    "\n",
    "def Teb(X):\n",
    "    \"Temperature from energy balance, as a function of conversion.\"\n",
    "    Cps = np.array([141, 141, 161])\n",
    "    Thetas = np.array([1., 0., 0.1/0.9])\n",
    "    return T0 - DHrx * X / sum(Cps*Thetas)\n",
    "assert Teb(0) == T0\n",
    "assert Teb(0.5) > T0 # Exothermic\n",
    "\n",
    "def dXdV(X,V):\n",
    "    \"\"\"\n",
    "    dX/dV in a plug flow reactor is -rA/ FA0\n",
    "    \"\"\"\n",
    "    T = Teb(X)\n",
    "    CA = CA0*(1.-X)\n",
    "    CB = CA0*X\n",
    "    rA = -k(T)*(CA - CB/Kc(T))\n",
    "    return -rA/FA0\n",
    "    \n",
    "volumes = np.linspace(0,5,501)\n",
    "X = scipy.integrate.odeint(dXdV, 0., volumes)\n",
    "plt.plot(volumes,X,label=\"Conversion\")\n",
    "\n",
    "Xeq = Kc(Teb(X))/ (1+Kc(Teb(X)))\n",
    "plt.plot(volumes, Xeq,label=\"Equilibrium conversion\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xlabel(\"Volume\")\n",
    "plt.ylabel(\"Conversion X\")\n",
    "plt.show()\n",
    "\n",
    "rate = k(Teb(X))*(CA0*(1.-X) - CA0*X/Kc(Teb(X)))\n",
    "plt.plot(volumes, rate, label=\"Rate\")\n",
    "plt.xlabel(\"Volume\")\n",
    "plt.ylabel(\"Rate\")\n",
    "plt.show()\n",
    "\n",
    "for volume, conversion in zip(volumes, X):\n",
    "    if conversion>0.7:\n",
    "        break\n",
    "print(volume)\n",
    "print(conversion)\n",
    "\n",
    "\n",
    "plt.plot(X, FA0/rate)\n",
    "plt.ylim(0,10)\n",
    "plt.title(\"Levenspiel Plot\")\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"$\\\\frac{F_{A0}}{-r_A}$\")\n",
    "plt.plot((Xeq[-1],Xeq[-1]),(0,FA0/rate[-1]), ':')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}