physical/haber_process.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bcadaf0f-ea22-4840-b56d-84b2628d6f9f",
   "metadata": {},
   "source": [
    "$$\\require{mhchem}$$\n",
    "# Thermodynamics of the [Haber-Bosch process](https://en.wikipedia.org/wiki/Haber_process)\n",
    "\n",
    "$\\ce{\\frac{1}{2}N2(g) + \\frac{3}{2}H2(g) -> NH3(g)}$\n",
    "\n",
    "(This is a solution to KFKA05's Exercise K17.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "44a7d68e-2fd5-4ffa-b1e6-a884b4965c65",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.constants import R\n",
    "from math import exp, sqrt\n",
    "import scipy.integrate as integrate\n",
    "from scipy.optimize import fsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ab80fb85-ddbe-4926-ac35-8499f7299dfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Heat capacity coefficients a, b, c\n",
    "table = {\n",
    "    \"NH3\": np.array([29.75, 25.10e-3, -1.55e5]),\n",
    "    \"N2\": np.array([28.58, 3.77e-3, -0.5e5]),\n",
    "    \"H2\": np.array([27.28, 3.26e-3, 0.5e5]),\n",
    "}\n",
    "\n",
    "\n",
    "def Cp(T):\n",
    "    \"\"\"heat capacity of reaction as a function of temperature (J/K/mol)\"\"\"\n",
    "    a, b, c = table[\"NH3\"] - 1 / 2 * table[\"N2\"] - 3 / 2 * table[\"H2\"]\n",
    "    return a + b * T + c / T**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee632ab3-c21b-4233-afde-65d1d4287324",
   "metadata": {},
   "source": [
    "## Equilibrium constant, $K(T)$\n",
    "Here we calculate the equilibrium constant at different temperatures by numerically integrating\n",
    "the heat capacity, $C_p(T)$, to obtain contributions from the enthalpy and entropy:\n",
    "$$\n",
    "\\Delta H(T) = \\Delta H(T_0) + \\int_{T_0}^T C_p(T) dT\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\Delta S(T) = \\Delta S(T_0) + \\int_{T_0}^T \\frac{C_p(T)}{T} dT\n",
    "$$\n",
    "The reaction free energy is $\\Delta G(T) = \\Delta H(T) - T\\Delta S(T)$ and the equilibrium constant,\n",
    "$$\n",
    "K(T) = e^{-\\Delta G(T)/RT}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "429ba141-b164-40d6-be35-7de69cfccc6e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dG(298.15) = -16.48 kJ/mol   K = 771.1\n",
      "dG(800) = 38.78 kJ/mol   K = 0.002939\n",
      "dG(900) = 50.47 kJ/mol   K = 0.001177\n"
     ]
    }
   ],
   "source": [
    "def reaction_free_energy(T):\n",
    "    \"\"\"calculates the reaction free energy at a different temperature\"\"\"\n",
    "    T0 = 298.15  # reference temperature\n",
    "    #dG0 = -16.48e3  # reference dG, J/mol\n",
    "    dS0 = 192.45 - 1 / 2 * 191.61 - 3 / 2 * 130.684  # reference dS, J/mol\n",
    "    dH0 = -46.11e3  # reference dH, J/mol\n",
    "    dH = dH0 + integrate.quad(Cp, T0, T)[0]\n",
    "    dS = dS0 + integrate.quad(lambda T: Cp(T) / T, T0, T)[0]\n",
    "    dG = dH - T * dS\n",
    "    return dG\n",
    "\n",
    "\n",
    "for T in [298.15, 800, 900]:\n",
    "    dG = reaction_free_energy(T)\n",
    "    K = exp(-dG / R / T)\n",
    "    print(\"dG({}) = {:.4} kJ/mol   K = {:.4}\".format(T, dG / 1000, K))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ced86af6-48da-4c82-a8a6-92929d0b6fdd",
   "metadata": {},
   "source": [
    "## Composition at elevated pressure and temperature\n",
    "\n",
    "We now fix the temperature to 800 K and the pressure to $p$=80 bar.\n",
    "What is the equilibrium composition if the mole-fraction of NH3 is 0.03 and with an excess of hydrogen?\n",
    "$K$ is already known from the previous calculation.\n",
    "The total pressure, $p=\\sum_i p_i$ and $\\sum_i x_i = 1$.\n",
    "With $x_{H2}=1-x_{N2}-x_{NH3}$ we write\n",
    "\n",
    "$$ K(800~K) = 0.002939 = \\frac{\\overbrace{x_{NH3}p}^{p_{NH3}}/p^{\\circ}}{ (x_{N2}p/p^{\\circ})^{1/2} \\cdot (1-x_{N2}-x_{NH3})^{3/2} (p/p^{\\circ})^{3/2}  } $$\n",
    "\n",
    "and it's now only a matter of solving for $x_{N2}$.\n",
    "Although this is in priciple easily done, Python's numerical schemes have trouble finding all roots in one go.\n",
    "We will therefore need to start the solver with initial _guesses_ and as can be seen below, we obtain two roots.\n",
    "As hydrogen is in _excess_, only the first found root is valid, and thus $x_{N2}$ = 0.0189."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "261fbd28-14a6-47df-aa55-f89fc99b8a00",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0189242204053502, 0.6820767732127475)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K = 0.002939  # at 800 K (see above)\n",
    "p = 80  # (bar divided by 1 bar, i.e. the standard state)\n",
    "xNH3 = 0.03  # NH3 mole fraction\n",
    "\n",
    "def func(xN2):\n",
    "    xNH3 * p - K * (sqrt(xN2 * p) * ((1 - xN2 - xNH3) * p) ** (3 / 2))\n",
    "fsolve(func, 0.0)[0], fsolve(func, 0.95)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e0d89e67-b92f-4ad0-b85c-19fac0facc6d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.LineCollection at 0x108aefa00>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnxElEQVR4nO3dd3xUZdr/8c+VBiSkAAkl9A7SSUCwgW0tu2tbuyhYKLq6q27xcatr2fL4e3Rdd1dFQIqi4tr7yrqIFUgQIXSkJ0BCCyUJJJn790eCyyrCJGTmzMz5vl+vvMxkhpzvSeK55tznPtdtzjlERMR/4rwOICIi3lABEBHxKRUAERGfUgEQEfEpFQAREZ9SARAR8akErwPURWZmpuvUqZPXMUREokp+fv5251zW178eVQWgU6dO5OXleR1DRCSqmNmGI31dQ0AiIj6lAiAi4lMqACIiPqUCICLiUyoAIiI+pQIgIuJTvigAK7fu5d8rir2OISLyX7bvO8CriwoJBLxpyx9V9wHU19OfbeCNxUV8/pvveB1FRHwsEHAUFJXy/opi/r2yhMWbd+McdM1qSt+26WHP44sCkJGcSGl5JYGAIy7OvI4jIj5TUFjKtE/W8++VJWzfdwAzGNg+gzvP6sHpvVpyQps0T3L5pAAkEXCwt6KK9OREr+OIiE+sKd7HQ++t5K0lW0ltlMDIXi05o1cWp3XPokXTRl7H80kBaFJz0N9dflAFQERCbvOuMh6ZvZoXF26mSWI8PzqzOzed2pm0xpF1/PFHAag96O8qq6RjC4/DiEjM2rX/II/8azUz520Eg+tP7swtI7tGxLv9I/FJAUgCYHfZQY+TiEisenfpVn75cgG7yg5yeW47bjujO9kZTbyOdVQ+KQA1ZwCl5ZUeJxGRWLNr/0HueX0pry4qonebNKbdMIQ+2eGf0VMf/igAh64BlKkAiEjDOfSuf3fZQW4/qzu3jOxGUkL03F7liwKQ3uTQNQANAYnI8Sstq+Q3rxVE5bv+w/miACTEx5HaOEFnACJy3AoKS5nwdD5bSyui8l3/4XxRAOA/N4OJiNTX8ws28utXl9IiJYlZE4YzuEMzryMdF88KgJm1B6YDrYEAMNE590iotpfRJElDQCJSLxWV1fz21aU8n7eJU7pl8siVAyN2amddeHkGUAX8xDm30MxSgXwze885tywUG8tITtQQkIjU2aadZdz8TD4FhXu49fRu3HF2D+JjpKWMZwXAObcF2FL7+V4zWw60BUJUAJLYvKs8FN9aRGLU3FUl3DpzIQ6YdF0uZ53QyutIDSoirgGYWSdgEDDvCM+NA8YBdOjQod7byGiSqBvBRCRoz87fyK9eKaB7y6Y8cW0OHVukeB2pwXl+6drMmgIvArc75/Z8/Xnn3ETnXK5zLjcrK6ve2zm8I6iIyLcJBBx/fHsFd7+0hFO6ZfLChOExefAHj88AzCyRmoP/M865l0K5LXUEFZFjqais5iezvuDNJVu4+sQO3HtBHxLiPX+fHDJezgIyYDKw3Dn3UKi3p46gInI0O/YdYOz0PBZu3M0vzu/F2FO7UHOYil1elraTgWuBM8xsUe3H+aHa2OEdQUVEDre2ZB8X//0Tlhbt4bFrBjPutK4xf/AHb2cBfQSE7SesjqAiciQFhaWMnjIfBzw7bljU39xVFxExCygc1BFURL4ub/1Orp+6gNRGCcy46US6ZjX1OlJY+acAqCOoiBxmzspiJjydT3Z6E2bcdCJtI7x3fyj4pgCoI6iIHPLm4i3c/vzndG+ZyvQbh5IZA20d6sM3BUAdQUUEahq63f3SEnI6NmPS6CFfvTn0I98UAFBHUBG/m/LROu59YxkjemTx+KgcmiTFex3JU/4qAOoIKuJbkz5cy/1vLue8vq155MpBUdvDvyH5qwCoI6iILz05dy0PvLWc7/Zrw5+vHEhiDN/dWxe++ilkJCdpCEjEZybO/bLm4N9fB/+v89cZgDqCivjKEx98yR/eXsH3+rfhz1cMjOm+PvXhq59GM3UEFfGNx+bUHPy/PyBbB/9v4aufSPphHUFFJHZNnPslf3pnBRcMyObhywfo4P8tfPVTObwjqIjEpumfruf3b9UM+zykg/9R+eono46gIrFtVt4mfvPqUs7q3YqHNexzTL766agjqEjsev2LIv7nxcWc2j2Tv149SLN9guCrn5A6gorEptnLtnHH84vI7dicidfm0jjR33f4BstfBUAdQUVizoerS7jlmYX0aZvO5DG5vm/vUBe+KgDqCCoSWxas38nY6Xl0yUph2vVDSG3s38Zu9eGrAqCOoCKxY1nRHm6YuoDsjCY8fdOJX13jk+D5qgCAOoKKxIKNO8oY/dR8mjZKYMaNJ/q2n//x8l8BUEdQkahWvLeCUZPnUVkdYMaNQ325kldD8V8BUEdQkahVWl7J6CkL2L7vAE+NGUK3lqleR4pqPiwA6ggqEo0qKqsZOy2PNcV7eXxUDoM6NPM6UtTzVTdQqJkKqiEgkehSVR3g1pkLWbBhJ3+5chCn9cjyOlJM8N0ZgDqCikQX5xx3v7SE2cuLuffCvnx/QLbXkWKG7wpAenISTh1BRaLGw++t4oX8zfzozO5cO6yj13FiiqcFwMymmFmxmRWEa5vqCCoSPZ6Zt4G/vL+GK3Lbc8dZ3b2OE3O8PgOYCpwbzg02S1FHUJFoMHvZNn79SgGn98zigYv7YmZeR4o5nhYA59xcYGc4t5neRB1BRSLdwo27uPXZhfRrm87frhmsts4h4rufqjqCikS2tSX7uHHqAlqlNWbymCEkJ/lusmLYRHwBMLNxZpZnZnklJSXH/f0OXQPYtV9nACKRpmTvAUY/NZ84M6ZdP1QtHkIs4guAc26icy7XOZeblXX8c3/Tv7oIrDMAkUhSdrCKG6ctYPveg0weM4ROmSleR4p5EV8AGpo6gopEnuqA40fPLqKgsJRHrxrEwPYZXkfyBa+ngT4LfAr0NLPNZnZjOLarjqAikeW+N5Yxe/k27rmgD2ed0MrrOL7h6dUV59xVXmy3WbI6gopEiikfrWPqJ+u56ZTOXDe8k9dxfMV3Q0BQcx1AQ0Ai3nt36Vbue3MZ5/ZpzS/O7+11HN/xZQFQR1AR7y3atJsfP/c5A9pl8PAVA4mL041e4ebPAqCOoCKe2rSzjJumLSArtRGTRmshd6/48g6LwzuC6l2HSHiVlldy/dQFVFY7nhujuf5e8uUZgDqCinijsrav/4Yd+3l8VA7dWjb1OpKv+bIAqCOoSPg55/jta0v5cPV2Hri4H8O7tvA6ku/5sgCoI6hI+E3+aB0z523k5pFduTy3vddxBJ8WAHUEFQmv95Zt44G3lnNe39b87Ds9vY4jtXxZANQRVCR8CgpL+fFzn9O/bToPXa7pnpHEnwVAHUFFwmJraQU3Tcsjo0kiT16n6Z6RxpfTQNURVCT0yg9WM3Z6HnsrKvnHzSfRMq2x15Hka3xZANQRVCS0AgHHT15YREFRKZOuy6V3mzSvI8kR+HIICNQRVCSUHp69ireWbOWX5/fmzN7q7hmpfFsA1BFUJDReXVTIo++v4Yrc9tx4Smev48hR+LYAqCOoSMPL37CLn/1jMSd2bs59F/XFTDN+IplvC4A6goo0rM27yhg/I4826Y15fFQOSQm+PbxEDV9eBIaahnAaAhJpGPsOVHHTtDwOVAV4blwuzVKSvI4kQfBtic5o8p+OoCJSf4GA4/bnFrG6eB9/u3ow3Vqmeh1JguTbAqCOoCIN48F/rmT28m38+ru9Oa1HltdxpA58WwDUEVTk+L38+WYem/MlV5/YgdEndfI6jtSRbwuAOoKKHJ+FG3dx14tLGNalOb+7oI9m/EQh3xYAdQQVqb+i3eWMm55P67TGPHZNDonxvj2URDXfzgJSR1CR+ik7WMXY6XlUVFbz7NgTNeMnivm2ADRLrvmjVUdQkeAFAo6fzPqCZVv2MGX0ELq30oyfaObb87a0xjW1Tx1BRYL3yL9W83bBVn5xXm9O79XS6zhynHxbANQRVKRu3ly8hUf+tZpLc9px06nq8RMLPC0AZnauma00szVm9j/h3r46gooEp6CwlJ+8sIicjs144GL1+IkVnhUAM4sH/gacB5wAXGVmJ4QzgzqCihxbyd4DjJueR/PkJB4flUOjBK3qFSu8vAg8FFjjnFsLYGbPARcCy0KxsZEjR37ja8u37KE64Bg5LT0UmxSJegHnWL5lD/sPVNMnO43L3vXtvBHPzZkzp8G/p5dDQG2BTYc93lz7tf9iZuPMLM/M8kpKSho0QEJ8HFXqBSTyrdZt38/eiiq6tmxKSiMd/GONl7/RIw0ifuNo7JybCEwEyM3NrffR+kjV8zevFvDK54XMueec+n5bkZg16cO13P/mcn5/ZnfuPLuH13EkBI5ZAMysMfA94FQgGygHCoA3nXNLj2Pbm4H2hz1uBxQdx/ers86ZKeypqKJ4T4UWrBY5zAerSvj9W8s5r29rbj+zu9dxJESOOgRkZvcAHwPDgXnAE8AsoAr4o5m9Z2b967ntBUB3M+tsZknAlcBr9fxe9dInu2bsv6CoNJybFYloX5bs49aZC+nZOo3/u3wAcXGa8ROrjnUGsMA5d8+3PPeQmbUEOtRnw865KjO7FXgXiAemHOcZRZ2dkJ0GwNLCPZzRSwtXi5SWVTJ2Wh5J8XE8eV0OyUka949lR/3tOufehJphIOdcxeHPmVmmc64YKK7vxp1zbwFv1fffH6+mjRLonJnC0qI9XkUQiRhV1QFufXYhm3aVMXPsMNo1S/Y6koRYsLOAFpjZsEMPzOwHwCehiRReJ2SnaQhIBPjD2yv4cPV27r+oL0M6Nfc6joRBsOd3VwNTzGwONReCWwBnhCpUOPXNTufNxVsoLaskvbZDqIjfzFqwickfrWPMSZ24Yki9RnUlCgV1BuCcWwI8AEwATgdudc5tDmWwcOlz6DqAzgLEp/LW7+SXryzhlG6Z/Oq7vb2OI2EUVAEws8nA7UB/4HrgdTP7YQhzhc1/CoCuA4j/FO4uZ8LT+bTNaMJfrx5EghZ28ZVgf9sFwOnOuXXOuXeBYcDg0MUKnxZNG9EmvbHOAMR3yg5WMW56HgcqA0wanUtGshZ28ZugrgE45x7+2uNS4MaQJPJAn+w0CnQGID4SCDh++sJ/Fnbp1lILu/jRsW4Ee93Mvm9m37g6amZdzOxeM7shdPHCo092OmtL9lF2sMrrKCJh8ej7a3hryVbuPq+XFnbxsWMNAY2lpgXECjNbYGZvmdn7ZraWmruC851zU0KeMsT6ZKcRcLB8y16vo4iE3NtLtvDw7FVcMqgtY0/t4nUc8dCxbgTbCvwc+LmZdQLaUNMLaJVzriz08cKjb9ualhDLikrJ6djM4zQiobO0qJQ7Z33BoA4Z/P6SflrYxeeCugZgZp2BLc659bWPG5tZp0OPo12b9MY0S06koFDXASR2bd93gHHT80lvksgTo3JonKiFXfwu2FlALwCBwx4Har8WE8yMPtnpLN2imUASmw5UVTNhRj479h/gyety1f1WgOALQIJz7qu1E2s/j6k5Y33aprFq6z4OVgWO/WKRKOKc49evFJC3YRcPXjqAfu20Ap7UCLYAlJjZBYcemNmFwPbQRPJGn+x0DlYHWF2sC8ESW6Z8vJ5ZeZu57YxufH9AttdxJIIEWwAmAL8ws41mtgm4Cxgfuljh11d3BEsMmrOymAfeXMY5fVpxx1la1Uv+W7A3gn0JDDOzpoA552LubXKnFimkJMWztLAUctsf+x+IRLg1xfu4bebn9GydxkOXD9TCLvINwc4CagT8AOgEJByaOuacuzdkycIsLs7o3SZNZwASE3aXHeSmaQtolBjHpNG5WtBdjijYIaBXgQupWQpy/2EfMaVv23SWbdlDdaDea8+LeK6yOsAPZy6kaHcFT1ybQ9uMJl5HkggV7NuCds65c0OaJAKckJ1G2cFq1u/YT9espl7HEamX+95YxsdrdvD/LhtATkct7CLfLtgzgE/MrF9Ik0SAvrWLxGsYSKLVjM82MP3TDYw/rQuX5rTzOo5EuGALwClAvpmtNLPFZrbEzBaHMpgXurdqSlJ8XM2FYJEo89Hq7dzz2lLO6NWSn5/by+s4EgWCHQI6L6QpIkRifBw9W6fqDECiztqSfdzyTD5ds1J45MqBxGvGjwQh2DMA9y0fMadP7SLxzsXk7kkMKi2r5MZpeSTExzF59BBSG2ttawlOsGcAb1JzwDegMdAZWAn0CVEuz/TJTuO5BZsoKq3Q7AmJeJXVAW6Zmc/mXWXMHDuM9s2TvY4kUSTYG8H+6wKwmQ0mxu4EPqRPbWvopYWlKgAS0Zxz3PPa0q9m/AzppBk/Ujf1WgHaObcQGNLAWSJC79ZpxBlaIlIi3vRPN/DMvI1MGNFVM36kXoK9E/jOwx7GATlASX03amaXAfcAvYGhzrm8+n6vhtYkKZ6uWU1ZpkXiJYLNXVXC715fytkntOLn5/T0Oo5EqWOtCTyj9tPfAKm1H42AN6i5M7i+CoBLgLnH8T1Cpk92GksKdSFYItPqbXv54TML6dEqlT9foR4/Un/HOgPIMbOOwEbg0a89lwxU1GejzrnlQMQuRze4YzNeWVTE+h1ldM5M8TqOyFd27DvADdMW0CgxnsljhqjHjxyXY/31PA68Q82sn8OHaYyaWUEhX1HazMYB4wA6dOgQ6s0BMKJHFlBzmq0CIJGiorKacTPyKd5zgOfHD9ckBTluRx0Ccs79xTnXG3jKOdflsI/OzrmjHvzNbLaZFRzho05DR865ic65XOdcblZWVl3+ab11bJFCpxbJfLCq3pc5RBqUc467XlxM/oZdPHT5QAa2z/A6ksSAYKeB3lzXb+ycO6vucSLHiB5ZzMrbTEVltRbPFs/95V9reHVRET/9Tg++27+N13EkRtRrGqgfjOiZRXllNXnrd3kdRXzutS+KeHj2Ki4Z3JYfnt7N6zgSQzwpAGZ2sZltBoYDb5rZu17kOJphXVqQFB/HB6uKvY4iPpa/YRc/feELhnZqzh8u6RexEyckOnlSAJxzLzvn2jnnGjnnWjnnzvEix9EkJyUwpHMz5q7a7nUU8alNO8sYPyOPNumNefzaHBolaChSGpaGgI5iRI8sVm7by5bScq+jiM+UllUy5qn5VFY7Jo8eQvOUJK8jSQxSATiKET1aAjXTQUXC5WBVgAlP57NxZxlPXJtDt5ZanU5CQwXgKHq0akrrtMaaDiph45zj7peW8OnaHfzvpf0Z1qWF15EkhqkAHIWZMaJHFh+u3k5VdcDrOOIDf31/DS8u3MztZ3Xn4kFq8CahpQJwDCN6ZrG3oopFm3Z7HUVi3KuLCvm/92qme/74zO5exxEfUAE4hpO7ZhJnug4goTV/3U5+9sJihnVpzh8v6a/pnhIWKgDHkJ6cyKAOzXQdQEJmTfE+xs3Io13zJjwxKpekBP1vKeGhv7QgjOiRxeLCUnbsO+B1FIkxxXsrGPPUfBLijKljhpKerPV8JXxUAIIwokcWzsFHa3RTmDSc/QequGHqAnbsO8iUMUPo0ELr+Up4qQAEoV/bdJqnJPHBSg0DScOoqg7ww5kLWb5lL3+/ZjD922V4HUl8SAUgCHFxxqndM5m7uoRAQKuEyfFxzvGrVwqYs7KE+y/qy+m9WnodSXxKBSBIp3XPYvu+gyzbosXi5fg8+v4anluwidvO6MZVQ8OzyJHIkagABOnUHpkAmg0kx+Uf+Zt5qHau/51n9/A6jvicCkCQWqY2pk92mgqA1Nu/VxRz14uLOaVbpub6S0RQAaiDkT2zyN+wi+2aDip19PnGXdzyzEJ6t0nlsVGDNddfIoL+CuvgwoFtqQ44XltU5HUUiSJrivdxw9QFtExrxFNjhpLaWHP9JTKoANRBj1apDGiXzj/yN3sdRaLE1tIKRk+ZT3ycMf2GoWSlNvI6kshXVADq6NKcdizbsoelRaVeR5EIV1pWyegp8yktr2Tq9UPp2CLF60gi/0UFoI6+PyCbpPg4Xswv9DqKRLCKympumr6Atdv3MfHaHPq2Tfc6ksg3qADUUUZyEmef0IpXFhVysEprBMg3VVUHuO3Zz8nbsIuHrxjISd0yvY4kckQqAPVwaU47du4/yJyVxV5HkQgTCDh+/uJi3lu2jXu+34fv9c/2OpLIt1IBqIdTu2eSldpIF4PlvzjnuPeNZby0sJA7zurB6JM6eR1J5KhUAOohIT6OSwa15f0VxWoRLV/58+zVTP1kPTec3JkfndnN6zgix6QCUE8/yGlHVcDxqu4JEGDKR+t45F+ruTSnHb/6bm/d5StRQQWgnnRPgBzyQt4m7n1jGef2ac0fL+lHXJwO/hIdPCkAZvagma0ws8Vm9rKZZXiR43jpngB5p2Ard724mFO7Z/LIVQNJiNd7KokeXv21vgf0dc71B1YBd3uU47jongB/m7OymB89+zkD2mfw+KgcGiXEex1JpE48KQDOuX8656pqH34GtPMix/HSPQH+9fGa7YyfkU/3Vk2ZOmYoKY0SvI4kUmeRcL56A/C21yHqS/cE+M/8dTu5aVoenTNTmHHjiVrIXaJWyAqAmc02s4IjfFx42Gt+CVQBzxzl+4wzszwzyyspibxe/LonwF8WbtzF9U/NJzujMTNuPJHmKUleRxKpt5Cdtzrnzjra82Y2GvgecKZz7lsX2nXOTQQmAuTm5kbcgryH7gmY/NE6CneX0zajideRJESWbC5l9JT5ZKY2YubYYersKVHPq1lA5wJ3ARc458q8yNCQrjupE2bwxAdfeh1FQmRZ0R6unTKP9CaJzBw7jFZpjb2OJHLcvLoG8FcgFXjPzBaZ2eMe5WgQbTOa8IPB7XhuwSaK91R4HUca2PItexg1eR5NEuN5duwwneVJzPBqFlA351x759zA2o8JXuRoSDeP7Ep1wDFx7lqvo0gDWlpUytVPfkZSfBwzxw6jffNkryOJNJhImAUUEzq2SOHCAdk8M2+j+gPFiILCUq6ZVPPO//nxw+icqQVdJLaoADSgW07vRkVVNVM+Xud1FDlOizfv5uonPyMlKYHnxw/Xal4Sk1QAGlC3lk05v18bpn2ygdKySq/jSD19vnEX10yaR1qTRJ4bp2EfiV0qAA3s1tO7se9AFVM/We91FKmH/A27uG7yfJolJ/H8+OE6+EtMUwFoYL3bpHH2Ca2Y8vE69lboLCCafPrlDq6bPI8WTZN4frxm+0jsUwEIgdvO6EZpeSVPf7bR6ygSpNnLtjH6qflkZzTh+fHDaZOug7/EPhWAEOjfLoMRPbKY9OFayg5WHfsfiKde+byQ8U/n07t1KrPGD9dNXuIbKgAhctsZ3dix/yDPzt/kdRQ5ihmfrueOWYsY0qkZz4wdRjP19hEfUQEIkdxOzRnepQVPfPClzgIikHOOv/17Db9+dSln9mrJ1OuH0lQtncVnVABC6Kfn9KB47wEemb3a6yhyGOccf3x7BQ++u5KLBmbz2KgcGidqMRfxHxWAEMrp2JyrhrZn0kfrWFa0x+s4AhysCnDH84t4Yu5aRg3rwEOXDyRRyziKT+kvP8TuOrcXzZIT+cXLS6gORFw3a1/ZU1HJmKfm88qiIn76nR7cd2FfLeAuvqYCEGIZyUn8+nsnsGjTbmbO2+B1HN/aUlrO5Y9/yvx1O/m/ywZw6xndMdPBX/xNBSAMLhiQzSndMvnfd1ayTe2iw275lj1c/LdP2LyrnKeuH8IPcqJyCWqRBqcCEAZmxv0X9eVAdYB7X1/mdRxf+XjNdi5//FMcjlnjh3Nq9yyvI4lEDBWAMOmUmcKPzujGm0u28P6KbV7H8YVn5m1g9JSau3tfvuVkTshO8zqSSERRAQijcad1pVvLpvz6laW6NyCEKqsD/OqVJfzy5QJO7pbJrAnDyVZfH5FvUAEIo6SEOH5/cT8Kd5fr3oAQ2bn/INdOnsfTn21k/GldmDJmCOlNEr2OJRKRdOtjmA3t3Jwrh9TcGzCyZ0uGd23hdaSYsXzLHsZOz6N47wEevmIAFw/SxV6Ro9EZgAd++d3edM5M4YczF7J5V5nXcWLCOwVb+MFjn1BZHWDW+OE6+IsEQQXAA6mNE5l4bQ6V1QHGz8in/GC115GiVmV1gN+/tZwJTy+ke6tUXrv1FAa2z/A6lkhUUAHwSJespvzlykEs27KHu15cjHO6S7iuCneXc/kTnzKxtq3D8+OGqZWzSB3oGoCHTu/Vkp+d05P/fWclfbLTGD+iq9eRosa/lm/jzllfUB1wPHrVIL4/INvrSCJRRwXAYzeP6MrSwj386Z0V9GqTxogeulHpaCqrAzz47komzl3LCW3S+Ns1g+mcmeJ1LJGopCEgj5kZD17Wnx6tUrlt5kLWb9/vdaSItX77fq44bMjnpVtO0sFf5DioAESA5KQEnrwul7g4Y+z0PHbuP+h1pIgSCDimfbKe8x75kNXF+3j0qkHcf1E/9fAXOU6eFAAzu8/MFpvZIjP7p5n5fgC3ffNk/n7NYDbuLOOKJz5V07ham3aWMWryPH772lKGdG7OP+84TeP9Ig3EqzOAB51z/Z1zA4E3gN94lCOinNQ1k6nXD6VodzmXPv4JG3f49x4B5xzPzt/IuX+eyxebdvOHS/ox7fohtElXSweRhuJJAXDOHb48VgqgOZC1hndtwTNjh7GnvIrLnviE1dv2eh0p7DbuKGPMUwu4+6Ul9G+XwTu3n8ZVQzuof79IAzOv5p+b2QPAdUApcLpzruRbXjcOGAfQoUOHnA0b/LGoysqtexk1eR5V1QGm33Ai/dqlex0p5MoPVvPYnDU8PnctiXHGz8/txbXDOmrVLpHjZGb5zrncb3w9VAXAzGYDrY/w1C+dc68e9rq7gcbOud8e63vm5ua6vLy8BkwZ2dZv3881k+ZRWl7J5NG5nNglNvsGOed4u2ArD7y5nMLd5Vw4MJu7z+tN63Td1CXSEMJeAIJlZh2BN51zfY/1Wr8VAKhZyvCaSfMo3FXO3ef14rrhnWLqHfHqbXu55/WlfLxmB71ap/K7C/rEbKET8cq3FQBPbgQzs+7OuUP9kC8AVniRIxq0SW/CC+OH85MXvuCe15cxe3kxD17WP+ovhm7cUcZf/72aFxcWkpIUz+8u6MM1J3YgIV4zk0XCxZMzADN7EegJBIANwATnXOGx/p0fzwAOcc4xc/5G7n9jOYnxxn0X9eWCAdlRd2F0084yHn2/5sAfH2dcPbQDt57RjcymjbyOJhKzInYIqC78XAAOWb99P3fOWsTCjbv5bv82PHBRXzKSk7yOdUybdpbx1/fX8OLCzcTVHvhvHtlVzdtEwiCihoCk/jplpjBr/HCemLuWh99bxby1O7l5ZFeuGtqe5KTI+nUGAo65q0uYOW8j/1pRTHycMWpYRyaM6KoLvCIRQGcAUaygsJT73ljGvHU7aZacyJiTOjP6pI6enxEU763ghbzNPDt/I5t3ldM8JYnLctpx/cmddeAX8YCGgGJY/oZdPDZnDbOXF5OSFM/VJ3bgplO7hHV4Zef+g8xdVcI7BVuZvXwbVQHH8C4tuPrEDnynTysaJahvj4hXVAB8YMXWPTw250te/6IIMyOnQzNG9MxiRI8s+mSnNegFY+ccS4v2MGdlMe+vKGbRpt0EHGQ2TeKSwe24ckh7umQ1bbDtiUj9qQD4yMYdZTy3YCMfrCphaVFN142s1Eac1j2LU7tn0ikzhez0xmQ2bRTUPQXlB6tZU7yPldv2snrbXlZu28vSoj2U7D0AwIB26Yzs2ZIzerWkX9v0mLpPQSQWqAD4VPGeCuau3s6clcV8uHo7peWVXz2XGG+0Tm9Mm/QmtExtRFW1o6KqmorKaioqA1RUVrPvQBWFu8s59GeSlBBH16ym9GqdysndMhnRI4usVE3hFIlkKgBCdcCxatteCneVU1RaTtHuCraUllO0u5ySvQdIjI+jcWI8jRMP/TeelKR4Omc2pUerpvRonUrH5sm6WUskymgaqBAfZ/Ruk0bvNmleRxGRCKC3ciIiPqUCICLiUyoAIiI+pQIgIuJTKgAiIj6lAiAi4lMqACIiPqUCICLiU1F1J7CZlVCzglh9ZALbGzBONNA++4P22R+OZ587Oueyvv7FqCoAx8PM8o50K3Qs0z77g/bZH0KxzxoCEhHxKRUAERGf8lMBmOh1AA9on/1B++wPDb7PvrkGICIi/81PZwAiInIYFQAREZ+KuQJgZuea2UozW2Nm/3OE583M/lL7/GIzG+xFzoYUxD5fU7uvi83sEzMb4EXOhnSsfT7sdUPMrNrMLg1nvoYWzP6a2UgzW2RmS83sg3BnbGhB/F2nm9nrZvZF7T5f70XOhmRmU8ys2MwKvuX5hj1+Oedi5gOIB74EugBJwBfACV97zfnA24ABw4B5XucOwz6fBDSr/fw8P+zzYa97H3gLuNTr3CH+HWcAy4AOtY9bep07DPv8C+BPtZ9nATuBJK+zH+d+nwYMBgq+5fkGPX7F2hnAUGCNc26tc+4g8Bxw4ddecyEw3dX4DMgwszbhDtqAjrnPzrlPnHO7ah9+BrQLc8aGFszvGeA24EWgOJzhQiCY/b0aeMk5txHAOeeHfXZAqpkZ0JSaAlAV3pgNyzk3l5r9+DYNevyKtQLQFth02OPNtV+r62uiSV3350Zq3kFEs2Pus5m1BS4GHg9jrlAJ5nfcA2hmZnPMLN/MrgtbutAIZp//CvQGioAlwI+dc4HwxPNMgx6/Ym1ReDvC174+zzWY10SToPfHzE6npgCcEtJEoRfMPv8ZuMs5V13zBjGqBbO/CUAOcCbQBPjUzD5zzq0KdbgQCWafzwEWAWcAXYH3zOxD59yeEGfzUoMev2KtAGwG2h/2uB017w7q+ppoEtT+mFl/YBJwnnNuR5iyhUow+5wLPFd78M8EzjezKufcK2FJ2LCC/bve7pzbD+w3s7nAACBaC0Aw+3w98EdXMzi+xszWAb2A+eGJ6IkGPX7F2hDQAqC7mXU2syTgSuC1r73mNeC62qvpw4BS59yWcAdtQMfcZzPrALwEXBvF7wgPd8x9ds51ds51cs51Av4B3BKlB38I7u/6VeBUM0sws2TgRGB5mHM2pGD2eSM1ZzyYWSugJ7A2rCnDr0GPXzF1BuCcqzKzW4F3qZlFMMU5t9TMJtQ+/zg1M0LOB9YAZdS8i4haQe7zb4AWwN9r3xFXuSjupBjkPseMYPbXObfczN4BFgMBYJJz7ohTCaNBkL/j+4CpZraEmqGRu5xzUd0i2syeBUYCmWa2GfgtkAihOX6pFYSIiE/F2hCQiIgESQVARMSnVABERHxKBUBExKdUAEREfEoFQETEp1QAROrAzP5tZmfXfn6/mf3F60wi9RVTN4KJhMFvgXvNrCUwCLjA4zwi9aYzAJE6qG3Xa8CdwJXOuWoAM9t02JnBC2Z2sZk9aWavmtl3PIws8q1UAETqwMz6AW2AA865vbVfaw+8A3zv0Muccy8758YCY4ArvMgqciwqACJBql144xlqFuXYb2bn1D6VA8wBks2sG7DhsH/2K+Bv4cwpEiwVAJEg1HbYfAn4iXNuOTWNyO6pfToHyAfmUjM0tLC2W+OfgLedcws9iCxyTGoGJ3KczOx54CogjZp+7UOAs4HR1LQ1XhRrHUolNqgAiIj4lIaARER8SgVARMSnVABERHxKBUBExKdUAEREfEoFQETEp1QARER8SgVARMSnVABERHzq/wMNfdvK1CFP0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = np.vectorize(func)\n",
    "x = np.linspace(0, 0.97, 50)\n",
    "plt.plot(x, f(x))\n",
    "plt.xlabel(\"$x_{N2}$\")\n",
    "plt.ylabel(\"func(x)\")\n",
    "plt.hlines(0, 0, 1, color=\"black\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}