{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exam of Numerical Methods, 26 June 2024\n",
    "\n",
    "### Exercise 1\n",
    "\n",
    "The Lane-Emden differential equation is written as:\n",
    "\n",
    "\\begin{equation}\\large\n",
    "\\frac{1}{\\xi^2}\\frac{d}{d \\xi} \\left(\\xi^2 \\frac{d \\Theta}{d \\xi} \\right) + \\Theta^n = 0 ,\n",
    "\\end{equation}\n",
    "\n",
    "with conditions $\\Theta(0) = 1$ and $\\Theta'(0) = 0$. This equation is used to study the density profile of stars, assuming a polytropic equation of state, with index $n$.\n",
    "\n",
    "1. Solve the equation for $n=1$ using a Runge-Kutta integrator at 4th order, in the range $0 < \\xi < 10$. Plot $\\Theta(\\xi)$ and show that it is consistent with the analytical solution: $\\Theta(\\xi) = j_0(\\xi) = {\\sin(\\xi)/\\xi}$\n",
    "\n",
    "2. Now repeat the same procedure for $n=5$, comparing with the analytical solution $\\Theta(\\xi) = \\frac{1}{\\sqrt{1+\\xi^2/3}}$.\n",
    "\n",
    "You are free to use python libraries for ODE integration, but you get bonus points if you write your own Runge-Kutta solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEmCAYAAABh8itbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3hElEQVR4nO3deXxU5dn/8c812VdCFgRCFiggooJA2BQBsX0Ul1K1UtTHrUXEqtWnfarW2upjW7fS2sUFcAPFn2itVbSIqCyCECAgyL4nEMISSEIgezLX74+ZpEkIEMgkM5m53q8Xr2TOOTnnPiH55p77XOc+oqoYY4zxfw5vN8AYY0zbsMA3xpgAYYFvjDEBwgLfGGMChAW+McYECAt8Y4wJEBb4xueJSLaILGqlfd8hIioio1tj/8b4Egt84zEi0kNEpovIFhEpFZFCEdkkIjNF5DIvtmu0iDwhInHeaoOvE5HxIvKGiKwTkSr3H8F0b7fLeFawtxtg/IOIZACLgSrgTWAjEAH0Bq4FjgELvdS80cDjwAygqNG6t4DZQGWbtsj3/BQYCqwDdgLnerc5pjVY4BtPeRyIBAao6tr6K0TkPqCzNxp1OqpaA9R4ux0+4DYgT1WrReQFLPD9kg3pGE/pBRxpHPYAqupU1bzGy0VkooisEZEyETkqIvNFZERzDuYecpjRxPIGY/LubR53r97tXqci8kRT29fbT6KIvCgie0Wk0v3xRRFJOMnxxojI/4rIThGpEJFtInJ7c87lFOf4hHvf54rIUyKS6973OhG5qiX7bkxV96hqtSf3aXyP9fCNp+wEzhWR61X1g9NtLCLPAg8BK4FHgRhgErBQRMap6lwPtWsaEAtcB/wPcNi9/NtTtK0DsAzoCbwOrAEGAPcAY0RkiKoea/RlT+EawpoGVLi3nSEiO1T16xaew0xcQ2VTgFDgQeBDEemtqtn12t0RCGrmPo+pakUL22XaGQt84ym/B74H/FNEtgNLgVXAIlXdXH9DETkX+CXwNTBGVSvdy18FNgEvich33MMtLaKqy0XkW1yB/2H9gDyFh3C9Y7lXVV+q1+61wAvu9b9p9DVhwOB65/I+sAu4D9d5tsRh4Fp1z3QoIgtx/aG8G/hVve2+AdKauc87cV3TMAHEAt94hDtYBwG/AMbiCpQ7AURkKXC7qu5ybz4OEOC52oB07yPPPQTzAK4edVbbnUED1wH5wPRGy6cBT7jXNw78lxqdyz4R2YbrD0dL/bU27N37XiUix5rY9y243mU0x0YPtMu0Mxb4xmNUdT1wB4CIpAGjgInApcBHIjLIHYrd3V/SVOhscH/sgfcCvzuQ1XhM231BcyswsImv2dXEsiM0v8d9Kk3tuwBocD3BA0NHxs9Z4JtWoao5wJsi8hawBLgEGIJrqEda8dDe+pk+2fCTJ861WfsWkSSaP4Z/VFXLWtQq0+5Y4JtWpaoqIitwBX6ye/FO98fz631eq6/7Y1O92voKgPgmlvdoqhnNaGp9u3BdgA6u38sXkWBc9xWcrm3esgobwzenYIFvPEJEvgcsbDwMIiIRwH+5X25yf5wDPAv8UkQ+VdUq97ZdcAVRDq4LkKeyDRguIpGqWur++o7ur2/suPtjPJDdjNP5EFfl0ERgar3ldwFJuMbyfZGN4ZtTssA3nvI8kCAic4D1QCmQAtyMq1f8pnuMH1XdKiJ/xFXt8pWIvMt/yjKjgVuaUaHzAjALWOAeNorDFcg5nHiTV6b747Mi8jZQDmxQ1Q007TngRuBFERmI64/PAOAnwFb3+rPirv9/HLhTVWec7X6a0pIxfBEZCYx0v8xwf7xPRIrc+/59y1pnfIEFvvGUn+OqvhkB3IArgI/iqnd/lkbDB6r6sIjswHVL/zO4pjZYAdysqktOdzBVfVtEuuIqe/wzrmGWJwEnrikC6m/7tYg8DEwGXsH1c/9//OcCceN9HxWRS9zbfB/Xu4aDuHr7jzdRg38mYtwf97VgH61hDP+5Qa3WL+p9boHvB8QeYm5M2xGRNbhuehrl7baYwGM9fGPaiIh0AvrT6B2IMW3FevjGGBMgbPI0Y4wJEBb4xhgTIHx6DD8xMVHT09O93QxjjGk3Vq9efVhVk5pa59OBn56eTlaWt6ZTMcaY9kdEck62zoZ0jDEmQFjgG2NMgLDAN8aYAGGBb4wxAcIC3xhjAoRHAl9EXheRQyLS5GRU4vI3EdkhIt+6ZyA0xhjThjxVljkD13S1b55k/Vhcz9/shWsekZdpxflEnE7lSEklidGhiLTmw5WMMUePHuXw4cNUVlaefmNz1kJDQ0lMTKRDhw5nvQ+PBL6qfiUi6afYZByu+dAVyBSROBHpoqr7PXH8+pxOZeRzC8k7WsZFKXG8N2k4hWVVJEaHoor9ITDGg8rLyzl48CDdunUjIiLCfq9aiapSVlZGbm4uYWFhhIeHn9V+2urGq2Rgb73Xue5lJwS+iEzC9SAMUlNTz/hA+4rKyC1yPapzzZ4iLvi/z6iocpIWH0lEaBDbDh4jIz2ed+4ahsNhP5zGtER+fj5JSUlERkZ6uyl+TUSIjIwkMTGR/Px8UlJSzmo/bXXRtqlkbXKaTlWdrqoZqpqRlNTk3cGn1K1jBEPSOxIk0LVDOOVVThTILihl84Fj1Cis3F3AquwC8o9VYLOFGnP2ysvLiY6O9nYzAkZMTAzl5eVn/fVt1cPPxfW4u1rdgLzWOJCIMHvScI6UVJIQFcJNr6wgK7uAC5JjKa10suOQ6/GmP5qeiQAXJnfgw3svsd6+MWehurqa4GCfnqHFrwQHB1NdXX36DU/29R5sy6nMwfV8zNm4LtYebY3x+1oOh5AUEwbAO3cNqxu3rx3DP3ysnKv/vhSnwrf7jvLbORuYOKIHaQmRNgZpzBmy35m209LvtUcCX0TeAUYDiSKSi+vZmCEAqjoVmAtcBezA9XDrOz1x3OaoH/4ikBQTRmJ0KIPT48nKLqBjVCizMvcwK3MP53eN5eP7Rlhv3xjjlzxVpXPTadYrcK8njuUJIlLX81dVhj/9JTUKG/OKeXnxTm4c1I2kmDDruRhj/ErA3mlb2/NPigkjIz2eIIEOESH88bOtDH3qS8ZPW47TaRd0jTH+I+CvttTv7dc4nVz8zAKcClnZheQWlpGaYOVmxhj/ELA9/Ppqe/vnxIYzOD0eh7hqRh989xt2HjpupZvGBKD58+czduxYEhISCA8Pp3fv3jz88MMUFhY2uf3999/Ptdde22DZggULEJET/iUnJwPw/PPP069fP5xOZ6ufD1jgN1Db21/x6Hd58aYBfLOniMv/vJgfvrzMhneMCSBPPfUUV1xxBeHh4bz66qt89tln3HPPPcycOZPBgwezd+/eBtvv3LmTadOm8fjjjzdYft555zF+/Hji4uKYOHEir732GvPnz2fevHkATJ48mUOHDjFz5sw2Oa+AH9JprLa3P6RHAiKgCqv3FLGvqIyUeBveMcbfLVy4kMcee4wHH3yQ559/vm75qFGjuO6667jooou47bbbWLhwYd26v/zlL/Tv35+MjIwG+5oyZQrr1q1j8+bNdO7c+YRjRUREcNtttzFlyhTuvLP1ixeth38StaWbtRWav/1oA/uLymx4xxg/99xzzxEfH8/TTz99wrr09HQeeOABFi1axIoVKwCoqKhg1qxZ3HzzzSdsP2PGDO6///4mw77WhAkT2LRpE8uWLfPcSZyE9fBPov7F3Hkb9vObjzZy8TMLGJzekdmThlutvjEn8X8fb2RTXrFX29C3ayyPX3v+GX9ddXU1ixcvZty4cYSGhjZ5V+t1113Hk08+yYIFCxg6dCiZmZkUFRVx6aWXnrBtXFwc06ZNIyUlhYEDB5KYmHjCxGcXXXQRsbGxzJs3j4svvviM23wmrId/CrXDO1de0KXuQm5WTiFHSmwaWGP80ZEjRygrKyM9PZ0nn3ySkJCQE/6lp6cD1I3jZ2ZmIiL069fvhP3NnTuXkJAQxo0bR0pKChERERw4cKDBNg6Hg379+pGZmdnq52c9/GZIjA4lI60jq7ILcSps2FfEZX3O8XazjPFJZ9Oz9hX1h2wnTZrENddcc9qvycvLIzY2ltDQ0AbL169fz4QJE6isrOSxxx5jwIABdOrUqcnhnaSkJLZt29byEzgNC/xmqJ2QbV9RGXe/tZqfzV7LjDuHMDA1zu7GNcaPJCYmEhERQXZ2Nl27dqVr164nbLN27VqAuimKy8vLCQsLa7BNaWkpV199NQMGDODdd9897fz1ERERlJWVeeYkTsGGdJrJ4RBS4iOZ+t8DKaus4YaXl3HjVLsb1xh/EhwczMiRI/n8889POg3xv/71LwDGjBkDQEJCwgm1+XPnzmXv3r0899xzzXpYSUFBAYmJiS1s/elZ4J+hiNBgnO63fattPN8Yv/PQQw9RUFDAo48+esK6nJwc/vrXvzJq1CiGDnU9pbVPnz5UVVWRm5tbt13t+H5zp47evXs35557rgdaf2oW+GeotlxTcF3Ezdx1xNtNMsZ40JgxY/j973/P888/z/XXX8+HH37I4sWLef755xk6dCgJCQm89dZbdduPHDkSgJUrV9YtGz16NMHBwVx77bVMnTqVL774gu3btzd5vKKiIrZt21a3n1alqj77b9CgQeqLamqcmldYqj94Yale8Pg8XbunUJ1Op7ebZUyb27Rpk7eb0GrmzZunV1xxhcbFxWloaKj27NlTH3roIS0oKDhh2yFDhugdd9zRYNmCBQv0qquu0qSkJA0KClJAJ02adMLXzpo1S8PCwvTw4cPNatfpvudAlp4kU0V9+EaijIwMzcrK8nYzTir7cAlj/rQIp8IQq883AWjz5s2cd9553m6G182YMYMHHniA/fv3N/l8X1Xl5z//OTNnzqSgoKDBurFjx5KYmNjgXcOpnO57LiKrVTWjqXU2pNMCUWH/GZ+z+nxjAtett95KcnIyL730UpPrd+7cycqVK/ne977XYPnatWtZuHDhCXPwtBYL/Baorc+vHc8/WmaBb0wgCgoK4vXXX2+ydw8wfvx4unfvztSpUxssP3DgAG+88QY9e/Zsi2bakE5LOZ3K1gPHmPBKJj07RfPe3cMJsmEdEyBsSKft2ZCOFzkcwnldY/ntNX1ZnVPIzGXZ5B+rsEnWjDE+xwLfQ64fmMyo3kn8/t+bGPbUF0yYnmk3ZRljfIoFvoeICP97xbk4FWrUbsoyxvgeC3wPuqBrLKnxEQD0SIwiMTr0NF9hjDFtxwLfg0SEz/9nFD2SoiiprKG8qm2eU2mMMc1hge9hYSFBPH3dhewrKuPvC7bbBVxjjM/wSOCLyJUislVEdojII02s7yAiH4vIOhHZKCKt//BGLxraI4HrBiTz8qKddgHXGOMzWhz4IhIEvAiMBfoCN4lI30ab3QtsUtX+wGjgTyLi1wPcd4/sgWIXcI0xvsMTPfwhwA5V3aWqlcBsYFyjbRSIEdfTQqKBAuDEh0X6kXM7x5DS0XUBt2enaLuAa4zxOk8EfjKwt97rXPey+l4AzgPygPXAA6ra5BVNEZkkIlkikpWfn++B5nmHiDD/wZF0jYtAARvRMcZ4mycCv6l5BBrH2xXAWqArcBHwgojENrUzVZ2uqhmqmpGUlOSB5nlPRFgwv77qPLYeOMa7q/ae/guMMaYVeSLwc4GUeq+74erJ13cn8IF7uuYdwG6gjweO7fOuurAzg9M78qf5W9l9uMQqdoxpJ7Zv346I8Omnnzb7a+6//36uvfbaE5YvWLAAETnhX3JyMs8//zz9+vXD6Wz9Mm5PBP4qoJeIdHdfiJ0AzGm0zR7gcgAROQc4F9jlgWP7PBHh11edx5GSSsZMWWQVO8a0EwkJCSxfvpzLL7+8Wdvv3LmTadOmNTnV8Xnnncf48eOJi4tj4sSJvPbaa8yfP5958+YxefJkDh06xMyZMz19Cido3gMXT0FVq0XkPuAzIAh4XVU3ishk9/qpwO+AGSKyHtcQ0MOqerilx24vkjtG1k2hnJVdwJGSSpJiwk73ZcYYL4qPj2fYsGHN3v4vf/kL/fv3JyPjxIkqp0yZwrp169i8eTOdO3c+Yf1tt93GlClTuPPO1q1Y90gdvqrOVdXeqvodVf2De9lUd9ijqnmq+l+qeqGqXqCqszxx3PYiMTqUft06AJAQHWYVO8a0AxdffDHjx49v1rYVFRXMmjWLm2++ucn1M2bM4P77728y7AEmTJjApk2bWLZs2Vm3tznsTts2ICL866eXcMPAbhSUVJJzpNTbTTLGnILT6eTbb79l4MCBzdo+MzOToqIiLr300ibXx8XFMW3aNObMmUNubi7l5eUN1l900UXExsYyb968Frf9VCzw24jDITw89lxCghxMmb/Vplww5iScTvX678eWLVsoKSk5o8AXEfr169fk+rlz5xISEsK4ceNISUkhIiKCAwcO1K13OBz069ePzMxMj7T/ZCzw21CnmHB+fEk6n3y736ZcMKYJTqdy0yuZDH/6S6/+fqxZswagQeDv3LmTESNG0Lt3bwYMGED9p/Hl5eURGxtLaOiJw7Xr16/n+uuvp7i4mMcee4x//vOfLFmy5IThnaSkJPLyGhc4elaLL9qaM/PDQSm8uGhngykX7AKuMS5HSipZnVNItVO9+vuxZs0aUlNTSUxMrFs2efJk7rjjDiZOnMjnn3/OLbfcwpYtWxARysvLCQs7sZ2lpaVcffXVDBgwgHfffZfw8PCTHjMiIoKysrJWOZ9a1sNvY+mJkXVTLvQ+J8Yu4BpTT2J0KIPSOhLsEAaldfTa78fq1asb9O7z8/PJzMzk9ttvB+B73/te3XbgKuEsLCw8YT9z585l7969PPfcc6cMe4CCgoIGf2BagwV+GxMR5v7sUuIiQ0iIDsU1vZAxBly/H+/cNYzlv7qc2ZOGeeX3Q1VZu3Ztg8Dfs2cPXbt2JSQkpG5ZWloae/bsAaBPnz5UVVWRm5vbYF9797rusA8OPv1gyu7duzn33HM9cQonZYHvBTERIfx09HdYsv0wq7ILvN0cY3yKwyEkxYR5rTO0Y8cOiouLT3vBtv5F5ZEjRwKwcuXKBtuMHj2a4OBgrr32WqZOncoXX3zB9u3bT9hXUVER27Ztq9tPa7HA95Jbh6WTGB3Gn+dv83pFgjHmP5q6YJuamkpeXh5VVVV1y3JyckhNTQUgPT2dIUOG8PHHHzfY14ABA5g/fz7du3fnt7/9LVdeeSW9e/fm7rvvbrDdv//9b0JDQ7nuuuta67QAC3yviQgNYvKoHizfdYShVrFjjM9YsmQJaWlpdOnSpW5ZUlISQ4YMYcaMGQB8/vnnqCqDBg2q2+aee+7hgw8+oLS04X02l112Gf/+9785dOgQVVVVPPjgg/zjH/9osM2sWbO48cYbSUhIaL0TwwLfq64431WW5bSHpBjjdbt27eKdd97hzTff5MYbbzxh/dSpU3njjTfo3bs3v/zlL3n77bcbDDvdeuutJCcn89JLL530GDt37mTlypV1F30B1q5dy8KFC5ucg8fTrCzTi7p1jCA1PpI9BaX0PscekmKMN913332sWrWKH/3oRzzxxBMnrO/Vq9cppz4ICgri9ddfrxsSasr48ePp27cvf//73+uWHThwgDfeeIOePXu2qP3NYYHvRSLCpz+7lBHPLSAx2nsXqYwxrhLKlho2bNgpJ1xr6o/BlVde2eLjNpcN6XhZVHgwk0Z+h6+2H2bt3iJvN8cY48cs8H3ArcPT6BARwgsLTizXMsYYT7HA9wHRYcH8+JLufLH5EEu3H7YSTWNMq7DA9xG3DU8jSIRbX1thJZqmXbEOSttp6ffaAt9HVDsVp2qDp2IZ4+tCQkJafcIv8x9lZWUNpnc4Uxb4PiIxOpSLUuIAiI+yp2KZ9qFTp07s27eP0tJS6+m3IlWltLSUffv20alTp7Pej5Vl+ggR4Z/3XMwv3/+WOev2cbC4gs4dTj27njHeFhsbC3DCtAPG80JCQjjnnHPqvudnwwLfhzgcwoPf7cWHa/fx2tJd/Prqvt5ukjGnFRsb26IQMm3HhnR8TEp8JNf068L/W7GHHYeO29tkY4zHWOD7oLtH9qCksob/en6xVewYYzzGAt8HJcWEI7gmVbOKHWOMp3gk8EXkShHZKiI7ROSRk2wzWkTWishGEVnsieP6q8ToUPp0jgFcQzxWsWOM8YQWB76IBAEvAmOBvsBNItK30TZxwEvA91X1fODEuUdNHRHhk/tHcH5X14UwG9ExxniCJ3r4Q4AdqrpLVSuB2cC4RtvcDHygqnsAVPWQB47r14KCHNx7WU+yj5Qyf+MBbzfHGOMHPBH4ycDeeq9z3cvq6w10FJFFIrJaRG472c5EZJKIZIlIVn5+vgea135dcX5n0hIimfbVLqvWMca0mCcCv6lJ3BunUzAwCLgauAL4jYj0bmpnqjpdVTNUNSMpKckDzWu/ghzCxBHdWbu3iM83HbTQN8a0iCcCPxdIqfe6G5DXxDbzVLVEVQ8DXwH9PXBsv3f9gG4EO4S731ptJZrGmBbxROCvAnqJSHcRCQUmAHMabfMRcKmIBItIJDAU2OyBY/u90qoaapw2qZoxpuVaHPiqWg3cB3yGK8TfU9WNIjJZRCa7t9kMzAO+BVYCr6rqhpYeOxAkRocyIDUOgIRom1TNGHP2xJfHhTMyMjQrK8vbzfA6p1P5n3fXMm/jAZb/6nLioyz0jTFNE5HVqprR1Dq707YdcDiE+8b0pKLayduZOd5ujjGmnbLAbyd6nRPD6HOTmLk8h/KqGm83xxjTDlngtyMTR/Tg8PEK3l6xx0o0jTFnzAK/HRneI57I0CB+98kmfjRtuZVoGmPOiAV+O1JQWlU3nJOVU2glmsaYM2KB344kRocyKK0jADHhIVaiaYw5Ixb47YiI8O6k4Uwa2YOjZVXszD/u7SYZY9oRC/x2xuEQJo3sQWiwg9e/zvZ2c4wx7YgFfjuUGB3GdRcl88GaXAptHN8Y00wW+O3Uj0d0p7zKyStLbOpkY0zzWOC3U706RRMbHsxLi3Yy3ko0jTHNYIHfTh0pqeR4RTUAq61E0xjTDBb47VRidCgZ7hLNiJAgEqJCvNwiY4yvs8Bvp0SE2ZOG89AV51JSWcM3e496u0nGGB9ngd+OORzC7RenExMezOtf7/Z2c4wxPs4Cv52LCgvmpiGpzNtwgLyiMm83xxjjwyzw/cBtw9NQVaYu3mklmsaYk7LA9wNdO0QQFxHCm8tzuHGqlWgaY5pmge8HjpRUcrSsCoA1e6xE0xjTNAt8P1B/Fs2wYCvRNMY0zQLfD9SWaD5xbV/Kqmr4eucRbzfJGOODLPD9hMMh3DQ0lcToMF5faiWaxpgTWeD7kbDgIP57WCoLt+azy+bKN8Y0YoHvZ24ZmkZokMNKNI0xJ/BI4IvIlSKyVUR2iMgjp9husIjUiMgPPXFcc6KEqFBiI4J5LyuXH05dZiWaxpg6LQ58EQkCXgTGAn2Bm0Sk70m2exb4rKXHNCd3pKSy7qEo3+wpshJNY0wdT/TwhwA7VHWXqlYCs4FxTWx3P/BP4JAHjmlOIjE6lIz0eABCghx0jLQSTWOMiycCPxnYW+91rntZHRFJBq4Dpp5uZyIySUSyRCQrPz/fA80LLCLCO3cN49kbLqSi2smXW+zvqzHGxROBL00sazxw/BfgYVWtOd3OVHW6qmaoakZSUpIHmhd4HA7hhoHdSI6L4A2bRdMY4+aJwM8FUuq97gbkNdomA5gtItnAD4GXROQHHji2OYngIAe3DU8jc1cBG/NsrnxjjGcCfxXQS0S6i0goMAGYU38DVe2uqumqmg68D/xUVT/0wLHNKUwYnEpESBBTF1mJpjHGA4GvqtXAfbiqbzYD76nqRhGZLCKTW7p/c/ZiwoOJCQ/m42/3c8PLVqJpTKAL9sROVHUuMLfRsiYv0KrqHZ44pjm9IyWVHDleAcDava4SzaSYMC+3yhjjLXanrR+rLdEUIMghxIQHebtJxhgvssD3Y7Ulmn+7aQBVNcrc9Qe83SRjjBdZ4Ps5h0O4pl8XenaK5rWlu+3irTEBzAI/AIgIP76kOxvzilmVXejt5hhjvMQCP0BcNyCZuMgQXlu6i/xjFdbTNyYAWeAHiIjQIG4anMpnGw8y7KkvmDA908o0jQkwFvgB5Jp+XQCoUVidYw87NybQWOAHkL5dY0mICgWgf0ocidGhXm6RMaYtWeAHEBHhtdszABh7QWdEmpr3zhjjryzwA8xFqR0Zkh7PjGXZVNc4vd0cY0wbssAPQD8e0Z3cwjLmbzro7aYYY9qQBX4A+l7fc0iNj+S1JbutRNOYAGKBH4CCHMIdF6ezek+hlWgaE0As8APUmD6dACvRNCaQWOAHqLSESLp0CAfggq6xVqJpTACwwA9QIsL7k4cT5BAGpHa0Ek1jAoAFfgBL7hjJ9/t35b2svRwtq/J2c4wxrcwCP8BNvLQ7JZU1vLNyj7ebYoxpZRb4Ae78rh24pGcCb3y9m7yiMivRNMaPWeAbfnJJdw4WVzDi2QVWommMH7PAN1yQ3AEAp0JWdoGVaBrjpyzwDUkxYfRIjAKgZ6doK9E0xk9Z4BtEhLk/u5ROMWF0jAq1Ek1j/JRHAl9ErhSRrSKyQ0QeaWL9LSLyrfvfMhHp74njGs8JDw3irkt7kLmrgHV7i7zdHGNMK2hx4ItIEPAiMBboC9wkIn0bbbYbGKWq/YDfAdNbelzjeROGpBATHsy0xTttUjVj/JAnevhDgB2quktVK4HZwLj6G6jqMlUtdL/MBLp54LjGw2LCQ7hlaCpzNxywSdWM8UOeCPxkYG+917nuZSfzE+DTk60UkUkikiUiWfn5+R5onjkT4/q7/utsUjVj/I8nAr+pK3xNdgtF5DJcgf/wyXamqtNVNUNVM5KSkjzQPHMm+nSJoVNMGAAXJnewih1j/IgnAj8XSKn3uhuQ13gjEekHvAqMU9UjHjiuaQUiwruThuEQGNw93ip2jPEjngj8VUAvEekuIqHABGBO/Q1EJBX4ALhVVbd54JimFXVPiubqfl15OzOHo6U2qZox/qLFga+q1cB9wGfAZuA9Vd0oIpNFZLJ7s98CCcBLIrJWRLJaelzTun46+juUVNYwc3m2t5tijPGQYE/sRFXnAnMbLZta7/OJwERPHMu0jfO6xDKmTyfe+Ho3P7gomZT4CBveMaadszttzUlNHtWDwtIqRk9ZaCWaxvgBC3xzUt0TowGbVM0Yf2GBb04qMTqUPp1jAEiJj7QSTWPaOQt8c1Iiwr/vH8GFyR2orHZSVWNDOsa0Zxb45pSCghz84r96k3e0nA/W5Hq7OaYdczq1bo6m+p+fal3j7UzLeKRKx/i3Ub2T6NetAy8u3MGo3kl07hBuFTumWZxOdV37UeXW11ey7eAxkqLDqKh2UlRWRXiwg5jwEApKKqlRJdghiEBVjRITHkyIw0FhaSXpiVH8YdwFdI4Lp3tiFKpwpKSSxGibzvtMiC//5czIyNCsLCvZ9wXzNx5g0lurXXfgpsfzzl3DcDjsF82cyOlU1u87yqa8o0yZv+2UF/sFuOL8zny28QDKf+ZpOVUqhQc7iAoLprC0kguTO/D+3cMpKq+28HcTkdWqmtHUOuvhm2a5KCUOaFixk+Sec8cYp1NZlV3Ash2Hmb5kN2VVNQ3WC/CdTlHszi9hYGoc4nCwJqeQQWkdeemWAdz0ygpW5xQyMDUORFhT7/PV2QX06RLD5v3HcCqUVzspr3b9EVmXe5Tzn/iMqhrl3M4x/OueizleWWPhfxIW+KZZkmLC6N0pmm2HjpOaEGUVOwanU8kpKGXptnyenreF0sqGIe8A+qV0YMO+YgaldeT/TRxKQWkVidGhJwzJvHPXsLrX9dfVfp4QFdLgj4ICa3IK6RoXwd7CMgC2HDjGBU/Mx6lKn84xfPTTSzhaYT3/+mxIxzRbTY2Ta174muKyKhb+72hCg+2af6CpHZM/VlbJhFdWcOhYRYP1DqB/ahzrc4+eEPItDd3aYzf1hyAru4C0hCh2Hy6pGw4KCRKqa5T+KR14/+6LKSzzTDt83amGdCzwzRlZuPUQd76xij9cdwG3DE3zdnNMG3I6le+/sJSNecUNxthP1ZNvi3Ct/UNQP/xT4iPJPlJat02HiGCOlVeTkdaR2ZOG+/X1JxvDNx4zuncSA1PjeGHBDkb37kTXOKvY8XdOp7J2bxFTF+9kQ14x4BqTPz85li37jzUZ8m15fcfh+M/xaoeG6od/YkwYB4td70RWZhcyc3k2V5zfmS4BWG1mPXxzxr7als9tr6+0ip0AcOR4BWP/uoRDxypwCHSODedgcTkZ6fFt3pM/U/V7/j+ansnqnELCgh2UVTkBSIuP5JP7RlBe4/TZczgb1sM3HlU73YJV7Pivqmonry7dzcuLd1BcVg24evX/uvcSHCJe6cmfqfo9/3cnDedISSVOp5PhzyzAqZBTUMqA33+OU5WBqR15727/HuoBu9PWnIWkmDD6dnGFfpe4CKvY8TPrc4vo/+R8np23BVW4oGsswQ4hIz2eTjFhJMWEtbvecG34d4oNZ3B6PMEOoXenaKqd6uq45BTytwXbyS0s9eu7em1Ix5wVp1O5c8Yq1uQU8tVDl9ExykK/vSutqOaZeVt4O3MPNe5cCBJY9qvLG/Tq27umLvJGhQVTXO56J9M9MZJP77+UY+20nt+qdEyr2H7wGFf85St+fEl3Hrumr7ebY1rgm5xCbno1k/IqJ0nRoaR0jOTbfa7SytmThrW70Guu2vBXVYY9/SW1j3wIC3ZQVeNsl1U9NoZvWkWvc2K4YWA33lyew7iLkrkgOdZvg8FfVVbVMGX+Nl5buruuV19YWsUnPxvkV736k6kd6lFVBqfHu+v5I9l12FXSuTK7kIVbD9GvW5xffC+sh29aZG9BKSP/uBAUhnS3ip32JOdwCVf/fQnHK2pIiAolLT4wevUnU3+oZ8L0TLJyCgkSocqpCNC/Wwc++OklPv/zbT1802rCQ4JAXZNdWcVO++B0Ku+vzuXJTzZxvMI1HcLRsipevjUwevUnU7+qZ7a7qud4eRVj/rQYBdbmHuX/PtnEz7/bi8oabZffJwt80yKJ0aEMTI1j9Z4iosKCSYgK8XaTzCmUV9YwaspCDhZXEBUaRL/kDmza77pDtlM7rL5pLbXhnxgdypDurqGejlGhzFyWzduZOThV2+X4vgW+aRER4R+TL+blRTv54/ytLNiSz3f7nuPtZpkmHDhazl1vZtXddVpeVcMrt2cEdK/+dBpP7LZ4Wz53vLEKcI3vf7bxABnp8e3m+2d1+KbFHA5h0qge9EiK4g9zN7O/qMyva5nbo+U7DnPV35awM/84PTtFt/u6+rZU29sXEUb1TmJIekccAqFBDu55ew1D/vAFP3jxa6qrnT7/dC6PXLQVkSuBvwJBwKuq+kyj9eJefxVQCtyhqmtOt1+7aNu+fL7pAHe9aQ9J8SWqysxl2Tzx8SYALkzuwL/uCZyZI1tD7cXdsspqRk9ZVFfKeU5sGIePVZDh5Z/9U120bXEPX0SCgBeBsUBf4CYRaVyUPRbo5f43CXi5pcc1vqd/tziEhlMuGO9wOpW8ojIe/dd6nvh4U92TpDbvL6awrMp69S1Q2+NPiY9kcHo8QQLxUSEcLK6gRmFVdgEHist9srfviTH8IcAOVd0FICKzgXHApnrbjAPeVNfZZ4pInIh0UdX9Hji+8RFJMWFckBzL+n3FxEeF2ZQLXuJ0Kj+cuow1e4oAmDyyB2v2FLJmTxGD0jra/4uH1B/fT4gK4doXvmZjXjFOhe/+eTHlVTU+907XE4GfDOyt9zoXGNqMbZIBC3w/IiJ8dO8IHvtoA++s3MO63KN1j0Y0bWd1TmFd2DsEfnJpDx6KCrWHfreC+qWcH983gsPHK/h800F+/eEGAFbuLmBVdgE9kqJ94nvviYu2TZ1B4/cxzdnGtaHIJBHJEpGs/Pz8FjfOtC2HQ/jV2D4kRYfx2482UOP0rbe0/szpVD75No8fz1xFsEMIcl9LSYwObXDh0bQOh0PoFBvOzUNTGZzmurDrEPjR9EyG/OELrn9pGU4v/z54ooefC6TUe90NyDuLbQBQ1enAdHBdtPVA+0wbiwkP4ddXn8cDs9fy6pJdTBrZw4KmlTmdyuV/XszuwyVEhATxxf+MJCo8xCd6lYFGRHj3bteNWwXHK7jyr0tQ4Ju9Rfz+35u4dXg66QmRXvl/8UQPfxXQS0S6i0goMAGY02ibOcBt4jIMOGrj9/7tmgu7EBMezNOfbuGGl73fs/FXTqdyqLic332yid2HSwCorK4hKjzEevReVPuOqnfnGIZ0d1/YjQzh9a+zuWzKIkb9cSHHy6ra/MJui3v4qlotIvcBn+Eqy3xdVTeKyGT3+qnAXFwlmTtwlWXe2dLjGt9WUFpFaYVrutlv9hTZlAutwOlUfjR9OVnZhSiQFB1GQUlF3Y1AxvvqX9hVVYY//SU1CnsKyhjw+8+pcSqD0jrybhvdseuRO21VdS6uUK+/bGq9zxW41xPHMu1DYnQoGenxrMouwKnwbW4Rl59nd+B6UvaRElZlFwKuseKP77+EIIfDhnF8TP0ZOTPS41mdU0j3xCi2HzoOwKrsQl5atIMfDEgmOS6iVf/vbLZM02qcTuVAcTl3vLGSY+XVvHPXMNK8NHbpT5xOZWNeMb94by3bDx1H3BdnA3GGy/amqRk5I0ODOe5+N5zSMYI5915CtXLWf7jtASjGq7KyC/jh1OUINoVySzmdyvdfWMqGvGKCRJhxx2D6dI21Xn07VBv+9Z+zC653a6qQkdaRd8/iObuteqetMaeTlhCF0HAKZXN2Pt1wgA15xe5XSp+usXZxtp1q6jm7fTpHo+7pxtfsKfT474rNlmlaXWJ0KIPSOroeKOFwEBps4XSmnE5l1oocfvfJJiJCgqisrrGLs36i8R27tc/ZbY3/XxvSMW3C6VQWbctn4sxVXDegG38a39/bTWo3amqcjJ6yiL2FZcSEB7P4F6Op4ezHeI1vqx3qaY0xfBvSMW3C4RDG9OnEvZf15J9rcnlv1V6fm1jKF9U4lUc+WM/ewjIASiuqqQEbxvFjrXlXtAW+aVP3XdaTyNAgHvrntz5xq7kvK62o5iczVvGP1bl06RBOkGDDOKZFbAzftKni8moqqlzPUf1mbxEHj5XTpUOEl1vle/KLy7nsT4s4XlFDWnwkX/58lM1hb1rMevimTdXekFVbafbW8hzvNsjHOJ1KVnYB1728rO4B4/uKymwOe+MR1sM3bap+RcKf5m/lpUU76X1ODOMu6hrwYeZ0Ktf8fQmb9h8j2CH07RzLtkPHbA574zEW+KbN1V6U+u01fZmzLo8H313LjK9388FPLwnoG7LeXrGHTfuPAa5HE77x48H2gHHjURb4xmtKKmvqxvPX5h5lT0Ep6YlRXm5V26upcfLsvK1MX7KLmPBgSiuq6x4wbkFvPMkC33hN4wnWnvxkE6/clkFQAPXySyuqufS5hXV111/972WUVNVYr960Crtoa7ymdjx/xaPf5cnvn8+CLYd4fM7GgKjPdzqVjfuOcsPUZXW3zxeWVFJSVWMXZ02rscA3XlU7nv/fw9I4JyaMWZk5jJ6yyK/r851O5doXlnL135ey9cAxenWKJtghVmNvWp0N6RifcKSkksPHKwDIOVLKW5k53H5xuncb1UpeXbqbje4J0ASYNXGoXZw1bcIC3/iE2vH8rOwCosKCefKTTXTpEM6A1I5+E4RlFdU88sF6PlqXR4eIEI6XV9nFWdOmLPCNT6hfnx8e4uCWV1dw96zVCK6He7T3OfR35R/nmr8vpbSyhuS4cBb8fDTFFdV+88fMtA82hm98Ru14fkx4CH+6sT+q4FRY1Y7n0Hc6lbeW59SFPcDB4gqKK6rt4qxpcxb4xif17BTNoNQ4wBX6C7YcJP9YRbuq4DlaWsmQp77gNx9twCFC/24dCHaI3TlrvMbmwzc+y+lUcgvLePRf37J0xxEc4nrs2+xJZ/7Yt7bkdCrzNh7gyY83caC4HIAggWW/utwuzppWZ/Phm3bJ4RBSEyJ59ob+CLXDO4XsKSj12d5+cVkVw5/5kp++vYbC0krO6xxTV3LZKSbMhnGMV9lFW+PzusaFMzjd9YhEp8LYvy2hssr1iD9fuZhbU+Pk/TW5/Gn+Ng4dc5WXVtc4mfHjIdarNz7DAt/4PBFh9qThHCmpZNnOwzwwey3guph78Fg5wQ6HVwN124Fj3DhtOUfLqogMCWowy6WVXBpf0qLAF5F44F0gHcgGxqtqYaNtUoA3gc6AE5iuqn9tyXFN4Kmt4Pl+/6688XU26/YW4VS44vmvKHFPNtbWvf28wjKe/WwLH6/Lo/bG4IrqGpvl0vislvbwHwG+VNVnROQR9+uHG21TDfxCVdeISAywWkQ+V9VNLTy2CUAiwgf3XMzh4xV8+M0+nvp0CwArdxewI/84HSNDWy1oax8uXVpZzatLdjNrRQ6qcE5MGF3jIli/76j16o1Pa1GVjohsBUar6n4R6QIsUtVzT/M1HwEvqOrnp9u/VemYU1FVrn9pGWv3FqFAkENwOpWLUuL4x93DPfJIwNqQjwsP5toXv2brgWMoEOwQapxa9/nXj4yxXr3xCaeq0mlp4Bepaly914Wq2vEU26cDXwEXqGrxSbaZBEwCSE1NHZSTY4/AMydXG8jbDx7jlldXUPvTHBcZQnFZFYNSO/LOXcPOKPzrh/wPXl7GprxiQoKEyhrX3gX4+P5L+N0nm1mdU8igtI7MnjTMgt74hBYFvoh8gWv8vbFfAzObG/giEg0sBv6gqh80p+HWwzfNpapMmJ5JVnYBCdFhdZUyALHhwRyvqKZXpxheuz2D0GAHSTFhqFI3D73TqWw7dJxDxeX85qON7C0oxeHuxddKS4gkt6CUjPR4Zk8a1uDrLeyNr2jNHn6zhnREJAT4BPhMVf/c3P1b4JszUdszT4gK4UfTMlm9p5CEqFDyj584LUOQQwgSobLGSVCjYK8lQGq9kP9/E4dSUNryYSJjWtOpAr+lF23nALcDz7g/ftTEwQV4Ddh8JmFvzJmqreQBePfu4XXhf9MrK8jKLqBHUhQ78ktcc/Q4Fad7AMjpVARQXHci9k2OZct+V1ll45Cv3b8x7VFLe/gJwHtAKrAHuFFVC0SkK/Cqql4lIiOAJcB6XGWZAI+q6tzT7d96+MYT6vf8b3plBatzChmYGgcirGn0eVMhb0x70mpDOq3NAt94Wm34J0aHNhiDt/F44y9ac0jHmHal/rCPCE1+boy/ssnTjDEmQFjgG2NMgLDAN8aYAGGBb4wxAcIC3xhjAoQFvjHGBAifrsMXkXzgbGdPSwQOe7A57UEgnjME5nkH4jlDYJ73mZ5zmqomNbXCpwO/JUQk62Q3H/irQDxnCMzzDsRzhsA8b0+esw3pGGNMgLDAN8aYAOHPgT/d2w3wgkA8ZwjM8w7Ec4bAPG+PnbPfjuEbY4xpyJ97+MYYY+qxwDfGmADhd4EvIleKyFYR2SEij3i7PW1BRFJEZKGIbBaRjSLygLfb1FZEJEhEvhGRT7zdlrYiInEi8r6IbHH/nw/3dptam4j8j/tne4OIvCMi4d5uU2sQkddF5JCIbKi3LF5EPheR7e6PTT43vDn8KvBFJAh4ERgL9AVuEpG+3m1Vm6gGfqGq5wHDgHsD5LwBHgA2e7sRbeyvwDxV7QP0x8/PX0SSgZ8BGap6ARAETPBuq1rNDODKRsseAb5U1V7Al+7XZ8WvAh8YAuxQ1V2qWgnMBsZ5uU2tTlX3q+oa9+fHcAVAsndb1fpEpBtwNfCqt9vSVkQkFhiJ6znRqGqlqhZ5tVFtIxiIEJFgIBLI83J7WoWqfgUUNFo8Dpjp/nwm8IOz3b+/BX4ysLfe61wCIPjqE5F0YACwwstNaQt/AR7iP89KDgQ9gHzgDfdQ1qsiEuXtRrUmVd0HTMH13Oz9wFFVne/dVrWpc1R1P7g6d0Cns92RvwV+Uw8jDZi6UxGJBv4JPKiqxd5uT2sSkWuAQ6q62tttaWPBwEDgZVUdAJTQgrf47YF7zHoc0B3oCkSJyH97t1Xtk78Ffi6QUu91N/z0rV9jIhKCK+zfVtUPvN2eNnAJ8H0RycY1dDdGRGZ5t0ltIhfIVdXad3Dv4/oD4M++C+xW1XxVrQI+AC72cpva0kER6QLg/njobHfkb4G/CuglIt1FJBTXhZ05Xm5TqxMRwTWmu1lV/+zt9rQFVf2VqnZT1XRc/88LVNXve32qegDYKyLnuhddDmzyYpPawh5gmIhEun/WL8fPL1Q3Mge43f357cBHZ7ujYI80x0eoarWI3Ad8hutK/uuqutHLzWoLlwC3AutFZK172aOqOtd7TTKt6H7gbXenZhdwp5fb06pUdYWIvA+swVWR9g1+OsWCiLwDjAYSRSQXeBx4BnhPRH6C64/fjWe9f5tawRhjAoO/DekYY4w5CQt8Y4wJEBb4xhgTICzwjTEmQFjgG2NMgLDAN8aYAGGBb4wxAcIC35gzICIHRUSb+FclIpHebp8xp+JXd9oa05rcd7b+BvgF0BtYBywCdgJ5qlrqvdYZc3oW+MacmTuBWOByVV3g7cYYcyYs8I1pvp/geqLYpaq61NuNMeZM2Ri+Mc13Oa4nqlnYm3bJAt+Y5osEwt1j+ca0Oxb4xjTfp7geqrNIRG4RkfNFpLN7jnZjfJ5Nj2zMGRCRX+OaozzEvagMiFXVau+1ypjmsYu2xjSDiHQCXgauBj4BlgI5wEELe9NeWOAbcxruIZtPgE5AP1Xd5uUmGXNWLPCNOb3RwGDgdgt7057ZRVtjTq+7+6NV55h2zS7aGnMaItID+BbXhdrZwCqgENigquu82TZjzoQFvjHNICLnA78ERgBdgQj3qsmqOs1rDTPmDFjgG3MWRCQWV6//oKoO9XZ7jGkOG8M35iyoajFQAuR7uy3GNJcFvjFnSER6iMgnQF/gRW+3x5jmssA35syNwFW5c5OqfurtxhjTXDaGb4wxAcJ6+MYYEyAs8I0xJkBY4BtjTICwwDfGmABhgW+MMQHCAt8YYwKEBb4xxgSI/w8E2cDnHe7JGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEmCAYAAAByJWuvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAyv0lEQVR4nO3deXhU5f3//+d7sofsJGwBEqygIKJAgABWEbVCFS1UrdrFYhGoS11aq7ULav26tWrtT1EsKrh8xKUu1SpqRVxBNkEFVFQIBJAlC0sSss39+2MmMWQhE7JMMnk9rosrzDlnztxHwyt33vd97mPOOUREpOPzBLsBIiLSMhToIiIhQoEuIhIiFOgiIiFCgS4iEiIU6CIiIUKBLkFlZpvMbHErnfuXZubMbFxrnF+kvVGgS0DM7Agze8jMPjezYjMrMLN1ZjbfzE4OYrvGmdmNZpYUrDa0d2a22P+Drb4/WcFun7Sc8GA3QNo//z/6d4By4DFgLRADDAAmAfuAt4PUvHHALGAeUFhr3+PAAqCsTVvUPu0Grq5n+zdt3RBpPQp0CcQsIBYY6pxbXXOHmV0O9AhGoxrjnKsEKoPdjnaiyDn3RLAbIa1LJRcJRH8gr3aYAzjnvM65bbW3m9k0M1tlZiVmtsfM3jCzEwL5MH8pYF492w+qifuPmeXfvbFGGeHG+o6vcZ5UM7vfzLaYWZn/6/1m1rWBzxtvZr8zs6/NrNTMvjSziwK5lkNc443+cx9lZreaWa7/3GvM7IfNOfchPtNjZglmZq1xfgk+9dAlEF8DR5nZFOfc840dbGZ3AL8HlgE3APHAdOBtMzvbOfdqC7VrDpAATMZXTtjt3/7JIdqWCHwIHAk8AqwChgK/Bsab2Ujn3L5ab7sVX4lpDlDqP3aemX3lnPugmdcwH18p6+9AJHAV8KKZDXDObarR7mQgLMBz7nPOldbalg7sx3cdxWb2OnCDc+7z5jVf2hMFugTiFuA04N9mtgF4H1gOLHbOra95oJkdBVwLfACMd86V+bfPBdYBs83se/5ySLM455aY2Sf4Av3FmgF4CL/H9xvHZc652TXavRq4z7//z7XeEwWMqHEtz+GrPV+O7zqbYzcwyflXyTOzt/H9IJwB/KHGcR8DGQGecyq+MYUqG/3t/ARfCWoUvrafYmYnOOc+bc4FSPuhQJdG+YNzOPBbYCK+wJgKYGbvAxc556oG184GDLizKgD959jmL5Fcia9HvKLtruAgk4FdwEO1ts8BbvTvrx3os2tdy1Yz+xLfD4bmurcqzP3nXm5m++o590/x9a4DsbbmC+fc1Fr7nzOz/wCLgbvx/bCWEKBAl4D4e3G/BDCzDOAkYBrwfeAlMxvuD71+/resrec0n/m/HkHwAr0fsMI5V1Fzo3Ouwsy+AIbV8576ZoLkEXiP+VDqO3c+cFA9vwVKOwdxzr1nZu8CJ5tZjHOupCXPL8GhQJcmc87lAI+Z2ePAe8BYYCS+UkxrDrgF6/u1ofJQS1xrQOc2szQCr6HvCTCgN+Gb9pkMKNBDgAJdDptzzpnZR/gCPd2/+Wv/12Nq/L3KIP/XxuY+5wMp9Ww/or5mBNDUmr7BN8AbXrOXbmbh+ObVt9d52cs5/Bp6Q/oDFfj+e0sIUKBLo8zsNODt2mUKM4sBfuB/uc7/9T/AHcC1Zvaac67cf2xPfEGTg2+A71C+BEabWaxzrtj//mT/+2vb7/+agq/H2ZgX8c28mQY8WGP7JUAavlp6e3RYNXT/rJ79tQehzewMfD+IX3POHWixVkpQKdAlEPcAXf0DaZ8CxUAf4EJ8vdrHqmZKOOe+MLO/4Zst8q6ZPc130xbjgJ8GMMPlPuAJYJG/rJOEL3BzqHsT01L/1zvM7EngAPCZc+4z6ncncC5wv5kNw/fDZSjwK+AL//7D4p//PguY6pybd7jnqU8zaugnA3eb2cv4fvuowFce+xm+GTZXtUgDpV1QoEsgrsE3e+UE4Mf4AnYPvmlwd1Dr13vn3HVm9hVwKXA7vlvvPwIudM6919iHOeeeNLNe+KbW3Y0viG4GvPim3NU89gMzuw6YCfwL3/f0TXw3AFv73HvMbKz/mLPw9fp34Outz6pnDnpTxPu/bm3GOVraF8BK4EygOxAB5OK73ludc+2prdJMpodEi7QMM1uF76aek4LdFumc1EMXaQFm1g04jlq/QYi0JfXQRURChBbnEhEJEQp0EZEQEbQaempqqsvMzAzWx4uIdEgrV67c7ZxLq29f0AI9MzOTFSuCtZyHiEjHZGY5De1TyUVEJEQo0EVEQoQCXUQkRCjQRURChAJdRCRENBroZvaIme00s3oXOzKff5rZV2b2iX8FOxERaWOBTFuch28508ca2D8R30L5/fGtY/EArbiehdfryCsqIzUuErPWfDiOSOsoLS0lPz+fffv2UVnZ7GdlS4iIjIwkNTWVxMTEwz5Ho4HunHvXzDIPccjZ+NbDdsBSM0sys57Oue2H3aoGeL2OKQ98wJotezi2dyLPzxxDQUm5wl06jNLSUjZv3kxycjKZmZlEREToe1dwzlFSUkJubi5RUVFER0cf1nla4saidGBLjde5/m11At3MpuN70AF9+/Zt8gflFZXxSe4eHPBJ7h6O++sblJRVcnyfJJ6bOab6GAW8tFf5+fkkJyeTmpoa7KZIO2JmxMbGkpqayq5du+jTp89hnaclAr2+5Kx3CUfn3EPAQwBZWVlNXuYxNS6SEZkprNiUT/fEaLYV+p6ctWpzIaf/4132l1awc+8BsjJTeOqSbDwehbq0L/v27UNLXkhD4uPjycvLO+z3t0Sg5+J7HFmV3sC2FjhvHWbGU5dkk1dURtcuEZz/0FJW5hTQMymGhJgINuz0PV5y2cZ8NucX0yUqXL11aVcqKyuJiIgIdjOknQoPD6eioqLxAxt6fwu04T/A5Wa2AN9g6J7WqJ9X8XiMtPgoABZMH11dYgGYfP8HrPGXZH5wz7tUeL2MUG9d2hl1MKQhzf3eaDTQzewpYByQama5+B6CGwHgnHsQeBX4IfAVvocH1/dk9lZRM9wBnr90LHlFZXy8uYDpj68EYPmmfHbuO0CYx6PeuoiEtEBmuVzQyH4HXNZiLWqGqoA/bVB3hmcksyqnAK+DSfd9QP7+UtXWRSSkheQzRc2MZ2eMZtf+Uh5Y/BXzPvStNrliUz55RWUH9epFREJFyN767/EY3ROimTXpGDK6xgKQGBtJUkxI/gwTEQndQK9iZrz923FceUp/8ovK+P2/P2XH3gPo4dgiLe+NN95g4sSJdO3alejoaAYMGMB1111HQUFBvcdfccUVTJo06aBtixYtwszq/ElPTwfgnnvuYciQIXi93la/no4m5AMdfL31q08bwG9PG8ALH28l+9a3OP+hpXi9CnWRlnLrrbdy+umnEx0dzdy5c3n99df59a9/zfz58xkxYgRbtmw56Pivv/6aOXPmMGvWrIO2Dxw4kPPOO4+kpCSmTZvGww8/zBtvvMHChQsBmDlzJjt37mT+/Pltdm0dhnMuKH+GDx/u2trOvQdc5nWvuIzrXnH9rn/F7dx7oM3bIJ3bunXrgt2EVrFo0SJnZu6qq66qs2/jxo0uMTHRjRs37qDtl19+ucvKyqpz/DXXXOOOOuoot3379gY/79prr3WDBg1qfsPboca+R4AVroFc7RQ99Cq+O02TMcA5+GbnPnbtK1X5RaSZ7rzzTlJSUrjtttvq7MvMzOTKK69k8eLFfPTRR4BvTZsnnniCCy+8sM7x8+bN44orrqBHjx4Nft7555/PunXr+PDDD1vuIkJApxohNDMWTB/Nxt1FTHtsBT99eBnOOU1nlKC76eW1rNu2N6htGNQrgVmTjmny+yoqKnjnnXc4++yziYyMrPdOx8mTJ3PzzTezaNEiRo0axdKlSyksLOT73/9+nWOTkpKYM2cOffr0YdiwYaSmptZZrOr4448nISGBhQsXMmbMmCa3OVR1qh46+Orp3+sWx53nDKHC66h0301nFJGmy8vLo6SkhMzMTG6++WYiIiLq/Klav6aqjr506VLMjCFDhtQ536uvvkpERARnn302ffr0ISYmhm+//fagYzweD0OGDGHp0qWtfn0dSafqodeUlZFM/25xbNi5n24J0dXLB4gEw+H0jNuLmiXL6dOnc+aZZzb6nm3btpGQkEBk5MH/7j799FPOP/98ysrK+NOf/sTQoUPp1q1bveWXtLQ0vvzyy+ZfQAjptIFuZrx+1Yn87tk1PP/xVt7dsJuTBqQFu1kiHU5qaioxMTFs2rSJXr160atXrzrHrF69GqB6WdgDBw4QFXXwDX7FxcWcccYZDB06lKeffrrRNcFjYmIoKSlpmYsIEZ2u5FKTx2PcOuVYBnSP45qnV7N++14NkIo0UXh4OCeeeCJvvvkmBw4cqPeYF154AYDx48cD0LVr1zpz01999VW2bNnCnXfeGdADHvLz87WufC2dOtABoiPC+Of5QykoLuOH977HT+Ys0fx0kSb6/e9/T35+PjfccEOdfTk5Odx7772cdNJJjBrlezrl0UcfTXl5Obm5udXHVdXXw8MDKxxs3LiRo446qgVaHzo6faADdI3z/erngBU5BRogFWmi8ePHc8stt3DPPfcwZcoUXnzxRd555x3uueceRo0aRdeuXXn88cerjz/xxBMBWLZsWfW2cePGER4ezqRJk3jwwQf53//+x4YNG+r9vMLCQr788svq84iPAh3f/PSsjOTqRy+VVujBvSJNdcMNN7Bw4UKKi4uZOnUqP/jBD5g9ezYXXXQRK1asOOixapmZmYwcOZKXX365etvQoUN544036NevH3/5y1+YMGECAwYMYMaMGXU+67///S+RkZFMnjy5Ta6to7Bg1YyzsrLcihUrgvLZ9fF6HZ9s3cOF/1rKsL7J3H3ecaTFR2n9dGlR69evZ+DAgcFuRrswb948rrzySrZv305sbGyd/c45rrnmGubPn09+fv5B+yZOnEhqaupBvf5Q0dj3iJmtdM5l1bdPPXQ/j8c4vk8S1084mve/2k32bVrvRaQ1/fznPyc9PZ3Zs2fXu//rr79m2bJlnHbaaQdtX716NW+//XadNWBEgV7H6YN7YIBXNxyJtKqwsDAeeeSRenvnAOeddx79+vXjwQcfPGj7t99+y6OPPsqRRx7ZFs3sUDrtPPSGdIuPYkjvRNbk7iEhJkI3HIm0ouzsbLKzs+vdt2rVqnq3T5gwoTWb1KGph16LmfHCpWO57OQjKSgu5411O4LdJBGRgCjQ6+HxGFed2p+je8Qz66W1bNxdpBuORKTdU6A3ICLMw62Tj+XbvQcY//fFGiAVkXZPgX4IfVJifWunowFSEWn/FOiHkBoXybC+SQBER4aREhsR3AaJiByCAv0QzIxnZ47hprOOoai0kmdW5jb+JhGRIFGgN8LjMX4xOoOR/VK447XP2bBjnwZIRaRdUqAHwMy4adIxFJaU84N73tUAqYi0Swr0AKXGR2mAVETaNQV6gFLjIhmWkQxogFSkIdOmTcPMuOaaa1r9s+bNm4eZsWnTpia9b/Hixdx44414vd6Dtm/atAkzY968eS3XyBrMjBtvvLFVzl1FgR4gM+PZGaOrB0j//fHWYDdJpF0pKSnh2WefBeDJJ5+koqIiyC2q3+LFi7npppvqBHrPnj1ZsmQJZ5xxRpBa1nwK9CaoGiAdnpHMHa99zsZduoNUpMoLL7zA3r17+eEPf8jOnTtZuHBhsJvUJFFRUWRnZ5OW1nGfLaxAbyIzY9aZg8grKmP8XbqDVKTK/PnzSU5OZt68ecTExPDYY48dtP/GG2/EzNiwYQNnnHEGcXFxZGRkcPPNNx/UWz5w4ABXX301gwcPJi4ujh49ejBp0iQ+//zzQ37+mWeeybBhw+ps37hxIx6Phzlz5nDjjTdy0003ARAREYGZVT/zoKGSyzvvvMNpp51GYmIiXbp04bjjjuPhhx+u3r9gwQLGjx9PWloacXFxDB06lPnz5zfpv11L0WqLh6FnUkydAdK0+KjG3iYSsrZt28b//vc/pk+fTlpaGj/60Y94/vnnKSgoIDk5+aBjJ0+ezNSpU7n66qt5+eWXmTVrFn369GHq1KkAlJaWsm/fPv70pz/Rs2dP8vPzmT17NtnZ2Xz++ef06NGj3jZceumlnHHGGSxbtoyRI0dWb3/ooYfo0qULF154IXv27CE3N5eHH36Y999/n7CwsENe10svvcSPf/xjxo4dy5w5c0hNTWXt2rXk5ORUH/PNN99wzjnncP311+PxeHj33XeZNm0aJSUlzJw583D/kx4e51xQ/gwfPtx1VF6v102+/32Xcd0r7thZC11lZWWwmyQdxLp161rlvJWVXrdz7wHn9Xpb5fyNuf322x3gPvzwQ+eccwsXLnSAe+CBB6qPmTVrlgPcI488ctB7Bw8e7E477bQGz11RUeGKiopcXFycu/vuu6u3P/roow5wGzdudM45V1lZ6Y444gh38cUXVx9TVlbmunfv7mbMmFGnHeXl5Qd9zsaNGx3gHn30Ueec7995RkaGGz58eMD/xisrK115ebmbNm2aGzJkyEH7ADdr1qxGz9HY9wiwwjWQqyq5HAYz47mZY/jdDwaw90AFb67fGewmSSfm9Tou+NdSRgfxKVuPPfYY/fv3Z/To0QCceuqp9OrVq07ZBagz6Dh48GA2b9580LZnnnmGUaNGkZSURHh4OF26dGH//v188cUXDbbB4/EwY8YMFixYwJ49ewB48cUX2bFjR73PJW3MF198QU5ODtOmTcPjaTgqN2zYwAUXXEB6ejoRERFEREQwd+7cQ7a1tSjQD5PHY8w86XsM6B7HLf9dx4FyPVhagiOvqIyVOQVUeB0rcwra/B6J5cuXs27dOqZMmUJhYSGFhYXs27ePKVOmsGTJEr788suDjk9JSTnodVRUFAcOHKh+/fLLL/OTn/yEgQMH8n//93989NFHLF++nLS0tIOOq8+vfvUrvF5v9bNGH3zwQUaOHMnQoUObfF15eXkA9O7du8Fj9u/fz2mnncaaNWu4/fbbee+991i+fDkXX3wxpaWlTf7M5lINvRnCwzzMmnQMP537Ef98awPXnn6UHiotbS41LpLhGcmszClgeEZymz9lq2oA8I477uCOO+6os/+xxx7jlltuCfh8CxYs4MgjjzxocLK8vLzOg6Lr07VrV84991zmzJnD6aefzttvv83cuXMD/uyaUlNTAdi6teEpykuWLCEnJ4f33nuPE044oXp7sKZsqofeTKOP6EpybASzF3/NlNkfasaLtDkz46lLslnyh1NYMD27TTsVZWVlLFiwgFGjRvH222/X+XP88cfz+OOPN2l6b3FxMeHhB/c1H3/8cSorA/st+NJLL+Wzzz5j2rRpJCQkcP755x+0PyrKN4GhpKTkkOcZMGAAmZmZzJ07t8H2FxcXA74ZM1UKCgp46aWXAmprS1MPvZnyisrYW1IOwOothZrxIkHh8VhQvu9eeeUV8vLyuOuuuxg3blyd/TNmzODXv/41ixcvDvicEyZM4MUXX+Tqq6/mzDPPZOXKlfzzn/8kKSkpoPdnZ2czbNgw3n33Xa644oo6D6EeNGgQAHfddRcTJ04kLCyMrKysOucxM/7xj38wZcoUxo8fz8yZM0lLS2P9+vXs3LmTm266iTFjxpCQkMBll13GTTfdRFFREbfccgupqanVdfy2pB56M6XGRZKVmVI9jTEnryjYTRJpM/Pnzyc+Pp5zzz233v0XXHABMTExTZqXfckll/DHP/6Rp59+mkmTJvHf//6Xl19+mcTExIDPcc455wDUOxh65plncumllzJ79mxGjx7NiBEjGjzP2WefzZtvvgn46vNnnXUWDz30EJmZmQCkpaXxwgsvUFlZyTnnnMMf/vAHpk2bxs9+9rOA29qSLJBfhcxsAnAvEAbMdc7dXmt/IvAE0Bdfr//vzrlHD3XOrKwst2LFisNtd7vi9Tq2FBRz3pwldE+I5sVLx+LxqJYuda1fv56BAwcGuxkhb+zYsXg8Ht57771gN6XJGvseMbOVzrm6v1IQQA/dzMKA+4GJwCDgAjMbVOuwy4B1zrnjgHHAXWbWtiMzQeTxGBldu/CHiQP5JHcPj36wUUsCiLSx0tJSlixZwl//+lc+/PBDrr322mA3qc0FUnIZCXzlnPvGOVcGLADOrnWMA+LNNxoTB+QD7XNlnlY0aUhP4qLC+et/13POgxogFWlL27dvZ8yYMdx9993ccMMNnHXWWcFuUpsLZFA0HdhS43UuMKrWMfcB/wG2AfHAT5xz3lrHYGbTgekAffv2PZz2tmv5xeWUlPl+jq3K0QCpSFvKzMzs9L8ZB9JDr68YXPu/2unAaqAXcDxwn5kl1HmTcw8557Kcc1kdeUWzhtQcIAXYU6KHYIhI2wkk0HOBPjVe98bXE69pKvC8f6mBr4CNwNEt08SOo2o+8MKrTvSVXl5Z3+l7DFKXviekIc393ggk0JcD/c2sn3+g83x85ZWaNgOnAJhZd+Ao4JtmtayD8niMo3rEc+Wp/Xnny1288PFW/QOWapGRkY3e0CKdV0lJyUE3KTVVo4HunKsALgdeB9YDzzjn1prZTDOrWhvyr8AYM/sUeAu4zjm3+7BbFQJ+np1BdISHa55Zw3lzlmiAVADf7eS5ubnk5+dTXl6uH/YC+HrmxcXFbN26lW7duh32eQK6U9Q59yrwaq1tD9b4+zbgB4fdihC090AFZRW+ceGqBZM0QCqJiYlERUWxa9cu8vLy2u1j2qTtRURE0L17dxIS6gw/Bky3/reS1LhIRmSmsGyjb0Gh8gDXoZDQFx0dTZ8+fRo/UKSJdOt/K6kaIH3hsrFEhHm47bW2XxtZRDoXBXor8niM4/skMfOk7/Hymm0s/SYv2E0SkRCmQG8Dvx73PdKTYvjLS5+xfU+JBsJEpFUo0NtAdEQYfzpjIF/u2M/Y2xcF7TFhIhLaFOhtZHhGMgZ4HazYlN/mjwkTkdCnQG8jafFRDOntW885KTayzR8TJiKhT4HeRsyMFy4dy8VjM8krKquezigi0lIU6G3I4zGuPf1o0pNi+PNLn7G9UAOkItJyFOhtLCYyjL+cOcg3QHqHBkhFpOUo0INgaN8kDZCKSItToAdBWnwUx/kHSBNiIjRAKiItQoEeBGbG85eO5fKTj6SguJy31u8MdpNEJAQo0IPE4zGuPLU/A7rHMes/a9mcV6wBUhFpFgV6EEWEefjr2YPZWljCSX97WwOkItIsCvQgOyItDsP3kFYNkIpIcyjQgyw1LpJhfZMA35ovybGH//gpEencFOhBZmY8O3MMt04+lqKySh79YBO79pWqni4iTaZAbwc8HuOCkX04+ag0bnttPdm3/k/1dBFpMgV6O2FmXHPaUXgdVOqGIxE5DAr0dmRwegIZXWMByOjaRTcciUiTKNDbETPjratP4vg+SezeX8qufaXBbpKIdCAK9HYmPNzD3ecdR2mFlz+9+Bk79x7QAKmIBESB3g4dkRbHVaf25411O8i+7S0NkIpIQBTo7dSPhqYDWpFRRAKnQG+neiREc2z6dysydu2iG45E5NAU6O2UmfHSZWP5zSn9KSgu56U124LdJBFp5xTo7ZjHY1x5Sn+yMpL5y0trWbt1jwZIRaRBCvR2Lsxj3PHjIRSVVnDm//c+P5mzRAOkIlIvBXoHkBDjq587YEVOgQZIRaReCvQOIDUukqyM5OpldguLFegiUpcCvQMwMxZMH82rV36f5NhIrn56NdsKS1RPF5GDKNA7CI/HGNgzgVsnD+azbXsZe/si3XAkIgdRoHcwwzNSqksvy3XDkYjUoEDvYFLjIhmekQz4ZsB4LMgNEpF2Q4HewZgZz8wYzZPTRmEY1z73iRbwEhFAgd4heTzG2CNT+cPEo1n0+U4t4CUigAK9QztjSE8MLeAlIj4K9A4sLT6KoX2TAAgP8xARpoK6SGcWUKCb2QQz+8LMvjKz6xs4ZpyZrTaztWb2Tss2U+pjZjw3cwyP/HIElV7Htc+qni7SmTUa6GYWBtwPTAQGAReY2aBaxyQBs4GznHPHAOe2fFOlPh6PMf7oblw/8WjeXK8HYoh0ZoH00EcCXznnvnHOlQELgLNrHXMh8LxzbjOAc25nyzZTGnPWcb1UTxfp5AIJ9HRgS43Xuf5tNQ0Aks1ssZmtNLNf1HciM5tuZivMbMWuXbsOr8VSr7T4KIZlJAG+Xnul1xvcBolImwsk0Osbaav9+3w4MBw4Azgd+LOZDajzJucecs5lOeey0tLSmtxYaZiZ8eyMMSyYnk1EmIfLnvxY672IdDKBBHou0KfG695A7cfn5AILnXNFzrndwLvAcS3TRAmUx2NkH9GVO6YMYeXmAq33ItLJBBLoy4H+ZtbPzCKB84H/1DrmJeD7ZhZuZrHAKGB9yzZVApX9va5a70WkE2o00J1zFcDlwOv4QvoZ59xaM5tpZjP9x6wHFgKfAMuAuc65z1qv2XIoqXGRjMj0r5/uYNPu/cFukoi0AQtWjTUrK8utWLEiKJ/dGXi9jo27i/jV/OXsO1DBvKkjGZyegJluPhLpyMxspXMuq759ulM0RHk8xve6xfGvn2dRWFzOpPve59wH9TxSkVCmQA9xSV0icf5JSStzCti5rzTILRKR1qJAD3G+enoKHvMNkv7rvW+C3SQRaSXhwW6AtC4z46lLsskrKuO+RRt4+P2NpMVFMeOkI1RPFwkx6qF3Ah6PkRYfxR9/OJCkmAhuX/g5E/7xrurpIiFGgd6J7DlQwb4D5QB8sWM/727Q8gsioUSB3omkxkWSlZlCmEFUuIdrnlnDqpwCLQ8gEiI0D72T8XodeUVlFBaVMvGf71PhdQzrm8RzM8fg0ROnRdo9zUOXalX19KQuUdU981WbC8nJKw5yy0SkuRTonVRV+cVjvuU0f//vNRQdqGDXvlKVYEQ6KJVcOrGq8svSb3bzmwWrSYyOYN+BcrIyU3jqkmyVYETaIZVcpF5V5ZdJx6Vz/YSjKSwpp1JPPBLpsBToAsD0E4+gb0oMAMldIkmJjQhyi0SkqRToAvjuKF38u5OZedIR7N5fxh9f/Iydew+oni7SgejWf6nm8RjXTxyIx4zZi7/m6eVbGJGZzILpo1VPF+kA1EOXOn45JrPGE48K2LVfKzSKdAQKdKkjLT7quyceAfe8+SU7VH4RafdUcpE6zIwF00eze38pj36wiQfeUflFpCNQD13q5fEY3RKimTo2s3ot9eWbCvh274FgN01EGqBAl0PylV++e0DGH1/4lC35xSq/iLRDulNUGlV1R+nCz7bz55fWAjA8I4lnZ2hBL5G2pjtFpVmq7iidMLgnVfm9MqeQz7/dq7VfRNoRDYpKwKqeT7p8Uz4Ak2d/SEWlV2u/iLQTCnQJWM3nk36aW8jF830ls6q1X9Lio4LcQpHOTSUXaZKq8svJR3djSHoi4BssXZWTr/KLSJCphy6Hxcx48bKxbNi5n2ufXcOMJ1bhMRih8otI0KiHLofN4zGO6hHPfRcOwwCvg2Ub89mhueoiQaFAl2brkxJz0FIBv3t2Dd/s2q/yi0gb0zx0aRFVc9UXrd/B9c9/igOGpCfy4mVjVX4RaUGahy6trmqwdPzA7tVz1T/ZuocXPt6qwVKRNqJBUWlRVQ+fXrEpn6iIMH777Bo8BlkZWthLpLUp0KVF1ZyrXlJWwUl/W4zX+Rb22ri7iO91iwt2E0VClkou0uKqyi99UmIZkZlcXYL55aPLeH/DLpVfRFqJBkWlVVUNlubs3s+Fc5dRVukls2ssb11zEmFh6k+INJUGRSVoqnrrGalxVHq9AGzKK2bGEyvZuKtIvXWRFqRAlzZRNVgaZtA7OYb/rd/JyXct5sx/vo/Xq1AXaQkaFJU2UXOw1DnH6NveotLB2u17ufP1z/nlmH50T4jCTLNgRA6XeujSZqrKL2nxUdW99ZQuETz4zjdk3/YWk+//QL11kWYIKNDNbIKZfWFmX5nZ9Yc4boSZVZrZOS3XRAk1Vb31pTecymu/+X71LJjVuXt4+ION7Nx7QLV1kcPQaKCbWRhwPzARGARcYGaDGjjuDuD1lm6khJ6q3nq3hGhG+Hvr8dHh/L//rmfUrW8xZfaH6q2LNFEgPfSRwFfOuW+cc2XAAuDseo67Avg3sLMF2ychrmZv/c2rT6x+GPXHWwp5bOkm9dZFmiCQQE8HttR4nevfVs3M0oHJwIMt1zTpLKp6691r9NbjosK48T/rGHXrW0xWb10kIIHMcqlv2kHtf13/AK5zzlUeapaCmU0HpgP07ds3wCZKZ1FzJkxFpZexdyzC62D1lkIeeOdrzhnem27xmgkj0pBAeui5QJ8ar3sD22odkwUsMLNNwDnAbDP7Ue0TOececs5lOeey0tLSDq/FEtKqeus9Er/rrSdEh/O3179g1K1vcdZ9mrcu0pBGb/03s3DgS+AUYCuwHLjQObe2gePnAa8455471Hl16780pmrZAK/Xy+jbfb11gF+MzuD3px9FSbmX1LhI9dilUznUrf+NllyccxVmdjm+2SthwCPOubVmNtO/X3VzaRVVvXXnHCP8S/KmdInisSU5LFi2hQqvV8vyitSgxbmkQ6jqrafGRfLW+p1Me8z3vWPAgunZHJEWp966dArN6qGLtAdVvXWAUwZ2Y2RmMityCjAzzn9oKQDD+ibx7Mwx6q1Lp6VAlw7HzFgwfTR5RWXk7y9lwr3v4YCVm32zYaYMS6dHQrR669LpaC0X6ZCqeuwDesQzst93d5r+7fUvGH3bIk7/x7uUl1fqeabSqaiGLh1eQ7Nh4qLCKCmrJCszhacuyVYpRkKCHnAhIa2+dWH6JMewv7SSSgfLNuazbFO+eusS8tRDl5BS1Vvv2iWC8+YsZdXmAt9255sRM6R3Iv+eOYaCknLNipEO6VA9dAW6hKyqcK85cArQNS6SwqIylWKkQ1LJRTql+gZO0+IiydtfVl2K+WRroUoxEjLUQ5dOoWYpZsoDH7Jmy56D9qsUIx2FSi4iNVSF+7d7Sjjrvg+qSzHJsRHsLSlXKUbaNZVcRGqoKsUMTk+sLsV0i4+ioLi8uhSzcO23eriGdDjqoUunVrMUc+6DS/h4SyEeMyr8k9mPTOvCK5efwL6ySpVipF1QyUUkAFXhXlJWwbi/L66+QSkizKj0Oob20VoxEnwquYgEoKoU0ycltvoGpSNSu1Be6fA631ox1z63huUb81WKkXZJPXSRetQsxVzwr49YsSmf5C6+KY8O31OU7pgyhOMzkrQQmLQplVxEmqEq3J1zjL7tLSpr/ZPpnRTDMzOyiQgPU51dWp3WQxdphppPTsrKTGFlTgHHpiewJncPXge5hSWMueNtDBjQPY5Xrvg+YR6rfiCHAl7ainroIk1QuxSzMqeAAd3jWL99X/V89h4J0YR5jG/3lJCVmcKTvxqlG5akxajkItIK6quz90uLo2uXSD7amF99XEbXWHLzi3XDkrQIlVxEWkHNx+I9dUl2dYkFYPL9H/DJ1j2Eezzk5BUDvhuW3lj3LUP7JtMtPkq9dWlx6qGLtIKq3ntyTDhn3vcBX3y7DzOq57b3Sozm4YuySI2PVilGmkQlF5Egqgr34tIKTr7ruxuWqvRJjuHxi0fSJTpC4S6NUslFJIiqZ8nERTLCP0tmYM94Ptu2F+dgS0EJ4+56B4C+KbE8cfFIYqLCFe7SZOqhi7Sh+mbJHNUjjnXb91H7n2LflBgenzqK2GiFu3xHJReRdqjBKZDf1g333skxPPyLLFLiohTunZwCXaSda0q490yM5t6fHE/frl3onqDZMp2NAl2kA6kv3I/uEc/a7XvrhHtaXCQ3TjqGwb0T6ZsSq3DvBBToIh1UfeE+qFc8n26tG+4J0eH8Znx/hmYkM6xvEs6h5QdCkAJdJATUF+4115SpKSrcQ2JMBLv3lzKsbzILLsnW8gMhQoEuEmLqC/chvRNZvaWwTrgDxEWFU1xWwaCeCTw/cwx7SisU7h2UAl0khNUX7sP6JuEFVuUUkBoXxc59pdXHh3kMr9eR2TWWR345gr4pseQXq/feUSjQRTqJqnBPjYusrqHXXDysd0osm/OKqfmvPjLcQ3mFl36pXZg/dQTRkeEHvV9B374o0EU6ufpWhjymVyKnH9Odv7/xJbVTIC0ukvjoCHLyirQEcDujQBeRajV78QDnP7SUFZvy6d89ji927K8zewYgOTaCPSXlDO6VyL9njqbwgGrwwaJAF5EG1VeDH9onkeJyL+u37yUhOoLCkvLq46tq8L2TY/jbj4fQu2ss6UkxKtG0EQW6iASkoRr8+Q8tZWVOAenJMWzJL6lToukSGUZcdDi79pUypHciz05XL761KNBFpFnq68UP7BnP2m17650mGRFmVFT6evG3Tj6WPimxZHSNVS++BSjQRaTFHGqa5Mc5BfRMimFrQd1efEyEhy5R4eQXlTGwZwLPXJJNcYVX4d5ECnQRaRWHmia5MqeAQf513+vrxRvggLS4KK4YfyR9u8aS3S+FyPAw9eIPodmBbmYTgHuBMGCuc+72Wvt/Clznf7kf+LVzbs2hzqlAFwldDfXiHb6bnfp2jWXT7uI6vXiAmIgwDpRX0jclltunHMvAXgmUVzrNjfdrVqCbWRjwJXAakAssBy5wzq2rccwYYL1zrsDMJgI3OudGHeq8CnSRzqGxXnzN9Wiqeu31SY6NICYijG/3HuDYdN/Aa9USBp0p6Jsb6KPxBfTp/td/AHDO3dbA8cnAZ8659EOdV4Eu0rk11IvHjFU5BRybnsgvx2Ry9TOr6y3ZePwP3U6LiyQmMpzcgmKO75PEM9NHV98EFYpB39xniqYDW2q8zgUO1fv+FfBaAw2ZDkwH6Nu3bwAfLSKhqupZqwBPXZJdpxdfdePT/y3bXB32lQ5Wb/5u+iTArv1lQBkAqzYXcvSshVRWOrolRBETEcaW/GKO75vE05eEdtBDYD30c4HTnXPT/K9/Dox0zl1Rz7EnA7OBE5xzeYc6r3roIhKIxko2Q/skUu6FT3IL6ZkYzbbCA/WWbaqmUnaLjyIm0hf0x3XAHn1ze+i5QJ8ar3sD2+r5kCHAXGBiY2EuIhKomj15Mw7Zq68d9KWVjs+27qFHQjTb9/iCfkeNlSdXbS7k6L8spMLrSIuLJDoijK2FJRybnsjTl2Szr6yywwQ9BNZDD8c3KHoKsBXfoOiFzrm1NY7pCywCfuGc+zCQD1YPXURaQyA9+rJKx6db9xyyR181QJsUE0F0RBg79x3gqO7xPHxRFuFhHtLio4IS9C0xbfGHwD/wTVt8xDn3/8xsJoBz7kEzmwv8GMjxv6WioQ+sokAXkbbUWNAP65tEudfxyZZCeifHsLmeJQ5qio0MIzo8jILiMvqlduHWyceS0TWWMP9vFDXDviWDXzcWiYg0IJCgr5o/3y+1C1/vLqp3RcqaYiI8REeEUVhcTt+UGKIiwvh6536GZSTz1LTmPQ5QgS4i0kSBBD1mrNyUz+D0RC4+oR9XP13/FMuaIsKMSq9jRGYKT12SjcfTtFBv7qCoiEin05TB2Ooplh9trhP2x6YnUlrp5Ytv99EtPoode0txwMqcAvKKyqrP2xIU6CIiTdBQ0EPDYV9fD394RnL1D4KWokAXEWkhDYV9Qz38lp4Zo0AXEWlDNUO/xc/dKmcVEZE2p0AXEQkRCnQRkRChQBcRCREKdBGREKFAFxEJEUG79d/MdvHdYl5NlQrsbsHmdBSd8bo74zVD57zuznjN0PTrznDOpdW3I2iB3hxmtqKx1RxDUWe87s54zdA5r7szXjO07HWr5CIiEiIU6CIiIaKjBvpDwW5AkHTG6+6M1wyd87o74zVDC153h6yhi4hIXR21hy4iIrUo0EVEQkSHC3Qzm2BmX5jZV2Z2fbDb09rMrI+ZvW1m681srZldGew2tSUzCzOzj83slWC3pS2YWZKZPWdmn/v/n48Odpvagpld7f/+/szMnjKz6GC3qTWY2SNmttPMPquxLcXM3jSzDf6vyYd7/g4V6GYWBtwPTAQGAReY2aDgtqrVVQC/dc4NBLKByzrBNdd0JbA+2I1oQ/cCC51zRwPH0Qmu3czSgd8AWc65wUAYcH5wW9Vq5gETam27HnjLOdcfeMv/+rB0qEAHRgJfOee+cc6VAQuAs4PcplblnNvunFvl//s+fP/A04PbqrZhZr2BM4C5wW5LWzCzBOBE4GEA51yZc64wqI1qO+FAjJmFA7HAtiC3p1U4594F8mttPhuY7//7fOBHh3v+jhbo6cCWGq9z6SThBmBmmcBQ4KMgN6Wt/AP4PeANcjvayhHALuBRf5lprpl1CXajWptzbivwd2AzsB3Y45x7I7italPdnXPbwdeBA7od7ok6WqDX9wC+TjHv0szigH8DVznn9ga7Pa3NzM4EdjrnVga7LW0oHBgGPOCcGwoU0YxfvzsKf834bKAf0AvoYmY/C26rOqaOFui5QJ8ar3sTor+a1WRmEfjC/Enn3PPBbk8bGQucZWab8JXWxpvZE8FtUqvLBXKdc1W/gT2HL+BD3anARufcLudcOfA8MCbIbWpLO8ysJ4D/687DPVFHC/TlQH8z62dmkfgGTv4T5Da1KvM9FvxhYL1z7u5gt6etOOf+4Jzr7ZzLxPf/eZFzLqR7bc65b4EtZnaUf9MpwLogNqmtbAayzSzW//1+Cp1gMLiG/wAX+f9+EfDS4Z4ovEWa00accxVmdjnwOr6R8Eecc2uD3KzWNhb4OfCpma32b7vBOfdq8JokregK4El/h+UbYGqQ29PqnHMfmdlzwCp8s7o+JkSXATCzp4BxQKqZ5QKzgNuBZ8zsV/h+uJ172OfXrf8iIqGho5VcRESkAQp0EZEQoUAXEQkRCnQRkRChQBcRCREKdBGREKFAFxEJEQp0ET8z22Fmrp4/5WYWG+z2iTSmQ90pKtJa/Hdm/hn4LTAAWAMsBr4GtjnnioPXOpHAKNBFvjMVSABOcc4tCnZjRJpKgS7i8yt8T4T6vnPu/WA3RuRwqIYu4nMKvqdhKcylw1Kgi/jEAtH+WrpIh6RAF/F5Dd8DUxab2U/N7Bgz6+Ffn1ukQ9DyuSJ+ZvZHfOtTR/g3lQAJzrmK4LVKJHAaFJVOz8y6AQ8AZwCvAO8DOcAOhbl0JAp06dT8JZVX8D1pfYhz7ssgN0nksCnQpbMbB4wALlKYS0enQVHp7Pr5v2p2i3R4GhSVTs3MjgA+wTcQugBYDhQAnznn1gSzbSJNpUCXTs/MjgGuBU4AegEx/l0znXNzgtYwkSZSoIvUYmYJ+HrtO5xzo4LdHpFAqYYuUotzbi9QBOwKdltEmkKBLlKDmR1hZq8Ag4D7g90ekaZQoIsc7AR8M18ucM69FuzGiDSFaugiIiFCPXQRkRChQBcRCREKdBGREKFAFxEJEQp0EZEQoUAXEQkRCnQRkRDx/wMoRmgWVQKqBgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate\n",
    "\n",
    "# RK4 integrator\n",
    "def rk4(f,specs):\n",
    "    # Solve ODE with 4th order Runge-Kutta\n",
    "    h = (specs['t1'] - specs['t0'])/specs['N']\n",
    "    x = np.array(specs['x0'],dtype='float')\n",
    "    Nequations = x.size\n",
    "    xt = np.zeros((specs['N']+1,Nequations),dtype='float')\n",
    "    xt[0,:] = x\n",
    "    ts = np.arange(specs['t0'],specs['t1'],h)\n",
    "    for i in range(len(ts)):\n",
    "        t = ts[i]\n",
    "        k1 = h*np.array(f(x,t))\n",
    "        k2 = h*np.array(f(x+0.5*k1,t+0.5*h))\n",
    "        k3 = h*np.array(f(x+0.5*k2,t+0.5*h))\n",
    "        k4 = h*np.array(f(x + k3, t + h))\n",
    "        x += (1./6.)*(k1 + 2.*k2 + 2.*k3 + k4)\n",
    "        xt[i+1,:] = x        \n",
    "    return ts, xt[:-1]\n",
    "\n",
    "# Function to calculate Theta and Theta'\n",
    "def g(theta,xi,n):\n",
    "    theta0 = theta[0]\n",
    "    theta1 = theta[1]\n",
    "    f0 = theta1  # theta' = y\n",
    "    if xi <= 1e-4:\n",
    "        f1 = -theta0**n\n",
    "    else:\n",
    "        f1 = -theta0**n -(2.0/xi)*theta1\n",
    "    return np.array([f0,f1])\n",
    "\n",
    "specs = {'x0':[1.,0.], 't0':0., 't1':10., 'N': 100}\n",
    "\n",
    "n = 1.\n",
    "# Function f(x,t) to pass to the ODE integrator, fixing n=1\n",
    "def f(theta,xi): return g(theta,xi,n)\n",
    "t, x1 = rk4(f,specs)\n",
    "\n",
    "# n=5 case\n",
    "n=5\n",
    "def f(theta,xi): return g(theta,xi,n)\n",
    "t, x5 = rk4(f,specs)\n",
    "\n",
    "# Plotting solutions for n=1 and n=5, and comparing with analytic expectations\n",
    "plt.plot(t,x1[:,0],label= r'$\\Theta(\\xi)$')\n",
    "analytical = [ ]\n",
    "def j0(t): \n",
    "    if t < 1e-20:\n",
    "        y = 1\n",
    "    else:\n",
    "        y = np.sin(t)/t\n",
    "    return y\n",
    "for t1 in t:\n",
    "    analytical.append(j0(t1))\n",
    "plt.scatter(t,analytical,s=5,label=r'$j_0(\\xi)$')\n",
    "plt.xlabel(r'$\\xi$',size=18)\n",
    "plt.legend(fontsize=16)\n",
    "plt.title('Solution, n=1',fontsize=18)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "def sol(t): return 1.0/np.sqrt(1.+t**2/3.0)\n",
    "analytical = [ ]\n",
    "for t1 in t:\n",
    "    analytical.append(sol(t1))\n",
    "plt.plot(t,x5[:,0],label= r'$\\Theta(\\xi)$')\n",
    "plt.scatter(t,analytical,s=5,label=r'Analytical')\n",
    "plt.title('Solution, n=5',fontsize=18)\n",
    "plt.xlabel(r'$\\xi$',size=18)\n",
    "plt.legend(fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "The Fermi-Dirac distribution is written as:\n",
    "\n",
    "\\begin{equation}\\large\n",
    "f(E) = \\frac{1}{\\exp{[(E-\\mu)/kT]} + 1}\n",
    "\\end{equation}\n",
    "\n",
    "The Fermi energy (chemical potential) $\\mu$, for fixed $T$, can be found by adjusting its value in such a way as to normalize the distribution to $1$ in the allowed energy range for the particle. That is, you impose $G(\\mu) = 1$, where:\n",
    "\n",
    "\\begin{equation}\\large\n",
    "G(\\mu) \\equiv \\int_{E_{min}}^{E_{max}} f dE\n",
    "\\end{equation}\n",
    "\n",
    "Assume your system is at room temperature $kT \\simeq {1/40}$ eV and that the energy of the particles is constrained between $0$ and $2$ eV. Find $\\mu$ for this case. Graphically verify that the result you found makes sense. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "root = 1.0000000000008642\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEOCAYAAACuOOGFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAj9klEQVR4nO3de5hV9X3v8feH4SYXRWG4K3ghChpRmKrRxrsEE9DmOe2ptifNyWlEW03Snp60ts95Amn/OO1pT05PGxvFxJr0NHraRBsgRIeYi6bGxIEoOgMoohEYZjOAchW5fc8fe6GbYQZmz6y919ozn9fz7Gdm/9Zv/fZ3B8OH9VuXnyICMzOz3hqQdQFmZtY3OFDMzCwVDhQzM0uFA8XMzFLhQDEzs1QMzLqALI0ZMyamTp2adRnla22FiROzrsLM+qGVK1dui4j6zrb160CZOnUqTU1NWZdRvpUrYfbsrKsws35I0i+72uYpLzMzS4UDpRY1NGRdgZnZcRwoZmaWCgeKmZmlwoFSixYuzLoCM7Pj5CZQJJ0p6YeS1khqlvS5TvpI0t9JWi9ptaRZJdvmSlqXbLu3utVX2aJFWVdgZnac3AQKcAj4o4iYDlwB3C1pRoc+NwPTktcC4CsAkuqA+5LtM4DbO9m37/A9KGaWQ7m5DyUitgBbkt93S1oDTAJaSrrdCnwjis/cf07SKEkTgKnA+ojYACDp0aRv6b59x5YtWVfwnvbd7/LMq+28sW1v1qWYWTcNGzKQu645N/VxcxMopSRNBS4FftZh0yRgY8n7TUlbZ+2XdzH2AopHN5x11lnpFNyPHDh0hKZf7uDpV7bx9CvttGzZ9d42KcPCzKzbxowY0j8CRdII4NvAH0TEro6bO9klTtB+fGPEYmAxQENDQ22uLjZr1sn7VMCXVrzC157ZwN4Dhxk4QMyacjqf/8j5XPOBemZMOJUBA5woZv1ZrgJF0iCKYfLPEfFYJ102AWeWvJ8MtAKDu2jvm1aurPpHbt/zLvf9cD1XnjuaT1wxhQ+dO5qRQwdVvQ4zy6/cnJSXJOBrwJqI+FIX3ZYAv5Nc7XUFsDM59/I8ME3S2ZIGA7clffumBQuq/pHfe7mNw0eCP/vodOZcON5hYmbHyU2gAFcBnwCul/RC8vqopLsk3ZX0WQ5sANYDDwK/DxARh4B7gCeBNcC/RERz1b9BtTz4YNU/cumLrZxbP5wLxo+s+mebWW3IzZRXRPyEzs+FlPYJ4O4uti2nGDiWssKu/fz8jR189vppyGfezawLeTpCsZxa/tIWImD+zAlZl2JmOeZAqUWbN1f145a+2MoF40dy3lhPd5lZ1xwotaiKV3ltfvsdVr35NvNn+u58MzsxB0otuuWWqn3Ud1cXr76ed7Gnu8zsxBwodkJLX9zCxZNPY8ro4VmXYmY550CxLr2xbS8vbd7poxMz6xYHSi164IGqfMx3Xyo+hPJjF/v8iZmdnAOlFlXpTvmlL7Yye8rpTBp1SlU+z8xqmwOlFlXh5sL1W3eztm23p7vMrNscKNappS9uQYKPfdCBYmbd40Cx40QEy1a3cvnZZzD21KFZl2NmNcKBUovmzavo8Gu27Oa19r3M88l4MyuDA6UWLV1a0eGXrW6lboC4+aLxFf0cM+tbHCi1aP78ig1dnO7awpXnjmb0iCEV+xwz63scKLVo2bKKDd2yZRdv7tjnq7vMrGwOFDtGY3OBAYIbp4/LuhQzqzG5WWBL0kPAPGBrRFzUyfbPA7+dvB0ITAfqI2KHpDeA3cBh4FBENFSn6r6nsaXA7Cmne7rLzMqWpyOUh4G5XW2MiL+OiEsi4hLgT4EfR8SOki7XJdv7fphEVGTYjTv2sWbLLubM8Ml4MytfbgIlIp4Gdpy0Y9HtwCMVLCffFi+uyLArWgoA3DTD011mVr7cBEp3SRpG8Ujm2yXNATRKWinphA+6krRAUpOkpvb29kqWWjl33lmRYRtb2jh/3EimjvGj6s2sfDUXKMB84N87THddFRGzgJuBuyVd3dXOEbE4IhoioqG+vr7StdaMt/Ye4Oev7/DRiZn1WC0Gym10mO6KiNbk51bgceCyDOqqaU+t3cqRgDkXOlDMrGdqKlAknQZcA3ynpG24pJFHfwfmAC9nU2GVLFmS+pArWtoYf+pQPjjptNTHNrP+IU+XDT8CXAuMkbQJWAgMAoiI+5NuHwcaI2Jvya7jgMdVfKT7QOCbEfFEterOxOzZqQ73zoHD/PiVdn5j9pmoCo/GN7O+KTeBEhG3d6PPwxQvLy5t2wDMrExVOTVpUqqXDv9k/Tb2Hzzi6S4z65WamvKyymhsbmPk0IFcfvborEsxsxrmQOnnDh8Jnlq7lesvGMvggf7Pwcx6zn+D1KI77khtqJW/fIsdew/4cmEz6zUHSi1K8U75xuY2BtcN4JoP+J4cM+sdB0otSukqr4igsaXAleeNZuTQQamMaWb9lwOlFq1alcowrxT28OaOfX4YpJmlwoHSjzU2twFw4/SxGVdiZn2BA6UWTUhnNcXGlgKXnjWKsacOTWU8M+vfHCi1qLW190O8/Q4vbd7p6S4zS40DpRYtWtTrIb6/xmufmFm6HCi16Itf7PUQjc0FzqkfznljR6RQkJmZA6Vf2rnvIM9t2O7pLjNLlQOlH/rhuq0cOhJ+GKSZpcqBUouamnq1+4qWAvUjh3DJ5FHp1GNmhgOl39l/8DA/WreVm2aMY8AAr31iZulxoNSihoYe7/rT17az98BhX91lZqnLTaBIekjSVkmdLt8r6VpJOyW9kLy+ULJtrqR1ktZLurd6VdeexpY2hg+u48pzvfaJmaUrN4FCcSXGuSfp80xEXJK8/hxAUh1wH3AzMAO4XdKMilZao44cCVa0bOXaC8YyZGBd1uWYWR+Tm0CJiKeBHT3Y9TJgfURsiIgDwKPArakWlzcLF/Zot19sfJtte95ljqe7zKwCchMo3fQhSS9K+p6kC5O2ScDGkj6bkrZOSVogqUlSU3t7eyVrrZwe3inf2NLGoDpx3QV+GKSZpa+WAmUVMCUiZgJ/D/xb0t7ZpUrR1SARsTgiGiKiob6+RheVmjixR7utaC5wxTmjOdVrn5hZBdRMoETErojYk/y+HBgkaQzFI5IzS7pOBnr/9MQ827Kl7F3Wb93Dhm17Pd1lZhVTM4EiabwkJb9fRrH27cDzwDRJZ0saDNwGLMmu0nxqbEnWPnGgmFmFDMy6gKMkPQJcC4yRtAlYCAwCiIj7gV8Hfk/SIeAd4LaICOCQpHuAJ4E64KGIaM7gK1TPrFll79LYXGDm5NOYcNopFSjIzCxHgRIRt59k+5eBL3exbTmwvBJ15dLKlWV1L+zazwsb3+bzHzm/QgWZmdXQlJeVWLCgrO5e+8TMqsGBUosefLCs7o3NBaaOHsY0r31iZhXkQOnjdu8/yLOvbWPOheNJrmkwM6sIB0of96N17Rw8HL5c2MwqzoFSizZv7nbXFS0FRg8fzKVnnV7BgszMHCi1qZtXeR04dIQfrt3KjdPHUee1T8yswhwoteiWW7rV7bkN29n97iEv9WtmVeFA6cMaW9oYNriOq84bk3UpZtYPOFD6qCNHgu+3bOXqafUMHeS1T8ys8hwoteiBB07a5aXNO2nbtd/TXWZWNQ6UWtSNO+UbW9qoGyCu99onZlYlDpRa1I0bFBubC1x+9hmMGja4CgWZmTlQ+qTXt+3l1a17/OwuM6sqB0oftCJZ+8SBYmbV5ECpRfPmnXBzY3OBCyeeyuTTh1WpIDMzB0ptWrq0y03tu99l5ZtvMWfG+CoWZGaWo0CR9JCkrZJe7mL7b0tanbyelTSzZNsbkl6S9IKkpupVnZH587vc9NSaAhGe7jKz6stNoAAPA3NPsP114JqIuBj4C2Bxh+3XRcQlEdFQofryY9myLjetaCkw+fRTmD5hZBULMjPLUaBExNPAjhNsfzYi3krePgdMrkphNWTvu4d4Zv025szw2idmVn25CZQy/S7wvZL3ATRKWinphHf9SVogqUlSU3t7e0WLrLanX2nnwKEjvjvezDIxMOsCyiXpOoqB8qslzVdFRKukscAKSWuTI57jRMRikumyhoaGqHjBlRCdl93YUmDUsEE0TPHaJ2ZWfTV1hCLpYuCrwK0Rsf1oe0S0Jj+3Ao8Dl2VTYZUs7nj6CA4ePsJTawrccME4BtbV1B+rmfURNfM3j6SzgMeAT0TEKyXtwyWNPPo7MAfo9EqxPuPOO49rev71Heza77VPzCw7uZnykvQIcC0wRtImYCEwCCAi7ge+AIwG/iE54XwouaJrHPB40jYQ+GZEPFH1L5CxxpYCQwYO4MPTvPaJmWUjN4ESEbefZPungU930r4BmHn8Hv1HRNDY3MaHp41h2ODc/JGaWT9TM1NeVmLJkmPeNrfuonXnfuZc6LvjzSw7PfrnrKQPABcCYylestsOvBwRr6ZYm3Vl9uxj3ja2FBgguMFrn5hZhrodKJKmA3cBv0HxvAXA0bvnIulTAP4FeCAi1qRYp5WaNOmYS4cbm9tomHIGo0cMybAoM+vvThooks4F/gr4OPAO8AzwU+A1YDvFUDkDOA+4guJ5js9Iegz4k+Qch1XIm9v3sbZtN//9Y9OzLsXM+rnuHKG0AC8B/xl4LCL2nqhzcunurwOfTfYd2ssa7QQavfaJmeVEdwLlP0bEd7o7YBI4Xwe+LunWHldmXbvjjvd+XdFS4ILxI5kyeniGBZmZdeMqr3LCJM197QSSO+V37D3A82/s8NGJmeWCLxuuRclVXk+tKXAk8GJaZpYLZQeKpL9JFruark6ekS5pkqRT0inPOrVqFVC8XHjCaUO5aNKpGRdkZtaz+1D+K8llwsA7klYDq4BfJK/bgOuAX0mlQuvUOwcO88yr7fxmw5le+8TMcqEngXIGMKvD6y6KRztHg2ZfKtVZ5yZM4JlX29l/8Ag3ebrLzHKi7ECJiLeBHyQvAJKn/X4UWASMAT6TTnnWqdZWGv/1RUYOHcjl55yRdTVmZkBKJ+UjYndE/D+K01x7KYaKVciRLyxM1j4ZyyCvfWJmOZHq30YRsQf4BvCHaY5rxxrwF3/OW/sOerrLzHKlEv+8bQcmVGBcKzG4bgDXnF+fdRlmZu8p+xyKpHaKV3OtOvqKiPXJNgEfA5rSLNLeFxEIuOq80YwY4rVPzCw/enKE8izwAeCPgUeBdZLelvQcsB64EvjfkgaVM6ikhyRtldTp8r0q+jtJ6yWtljSrZNtcSeuSbff24DvVjHWF3cz75N967RMzy52yAyUibo2IqRQvH74e+G/Ad4BTgDOBEcC3gD2SXpL0TUl/0o2hHwbmnmD7zcC05LUA+AqApDrgvmT7DOB2STPK/V61orG5gAQ3TPfaJ2aWLz2eM0kuH/5R8gJA0mDgIuCS5HUpxSmw36T4CPwTjfe0pKkn6HIr8I2ICOA5SaMkTQCmAuuPPiZf0qNJ35ayv1QNaGxpY9nDfwD/+LmsSzEzO0baV3kdiIhVEfFQRHw2Ij4cEadRnCLrrUnAxpL3m5K2rto7JWmBpCZJTfvWrwfp/dfKlcVXaduiRcUdJ058v+3oiokLFhzbt7UVli49ti15kOMxbfPnF9vmzz+2HYr9S9uWLi2Om7xf9tmr3/8ys2e/32/ixGLbokU1952Qip/r7+Tv5O+U/+90AoqSlf867SDdEBFPnXSkzve9MSK+X0b/qcCyiLiok23fBf5HRPwkef8UxfM45wAfiYhPJ+2fAC6LiJPeXNnQ0BBNTbVz/cDXn32DhUuaeeOv5h2zYqOZWbVIWhkRDZ1t684RyhOSfiBpXnK+4mQfNkjSxyU9DSwvt9gT2ETxHM1Rk4HWE7T3OY0tbZxbPxwWLsy6FDOz43TnHMqlwJeAJcA2SSuAn1NcAngHvLcE8DSKSwBfD5wONFI8j5KWJcA9yTmSy4GdEbEluYx5mqSzgc0UH075Wyl+bi7s3HeQ5zbsYMHV58DcRVmXY2Z2nJMGSkS8DMyR9CHg9yme8L6d9x8EeZSAXcBjwFci4vlyCpH0CHAtMEbSJmAhMCip4X6KRzsfpXhp8j7gU8m2Q5LuAZ4E6oCHIqK5nM+uBT9YV+DwkWDOjHHFec3WPnkQZmY1rNtXeUXET4GfJtNesyleoltPMVjagZeBX0TEkZ4UEhG3n2R7AHd3sW056U6v5c6KlgJjRw5h5uRRsGVL1uWYmR2nJ08bPkxxyuvn6Zdjndl/8DA/WtfOr106iQEDlHU5Zmad6tZlw5LOlLRD0qKStlMk3SXpVyQNqViFxrOvbWPfgcPF6S6AWbNOvIOZWQa6e4TyX4ChwP0lbcOAf6A45XVIUgvFZ3s9ERH/mmqV/Vxjc4ERQwbyoXNHFxtWrsy2IDOzTnT3xsaPAN+JiLZOtj0G/BCYSPFE+T8ld7BbCg4fCb6/psA159czZGBy1fbRm5bMzHKku4FyAfBcF9u+EhFzI2IcxUuHg+KjViwFL2x8i217Drw/3QXw4IPZFWRm1oXuBsoIYOfJOkXEa8BSYE5virL3NTYXGFQnrrvAD4M0s3zrbqDsBMZ1aNtB8b6Rjvd8/AL4YO/KMiiufdLYUuCKc0Zz6tCyVgMwM6u67gbKajocdUTR0xFR6NC3gNeUT8Vr7Xt4fdveY6e7ADZvzqYgM7MT6G6gfBu4VtIN3eg7Hjjc85LsqCebi1l9Y8dA8VVeZpZD3Q2Ur1F8dte/Jo9gOZE5FB+PYr3U2FJg5uTTmHDaKcduuOWWbAoyMzuBbgVKRLxL8RleB4AfS/pbSeeW9pE0QNIXgA8Dj6deaT9T2LWfFze+7aV+zaxmlPMsrzWSrgX+Cfgs8BlJb1A8GhlI8dleY4F1wP9Ku9D+ZkVLcbrrpo7TXWZmOVXWio0RsRa4jOKd808DU4CbgOuAUcA3gWsiYk+6ZfY/jS0Fpo4exrSxI47f+MAD1S/IzOwkevJwyAAeBh6WNIjicruDgI0RsT/d8vqnXfsP8tPXtvGpq85G6uRhkL5T3sxyqOxAKRURB4E30inFjvrRunYOHo6up7skLwFsZrlT1pSXVceKlgKjhw9m1lmnZ12KmVm35SpQJM2VtE7Sekn3drL985JeSF4vSzos6Yxk2xuSXkq2NVW/+nS8e+gwP1y7lRunj6POa5+YWQ3p1ZRXmpKVIO+jeJJ/E/C8pCUR0XK0T0T8NfDXSf/5wB9GxI6SYa6LiG1VLDt1z23YwZ53DzHnwhNc3TVvXvUKMjPrpjwdoVwGrI+IDRFxAHiU4r0vXbkdeKQqlVVRY3MbwwbXcdV5J3h6zdKl1SvIzKyb8hQok4CNJe83JW3HkTQMmEvxkTBHBdAoaaWkLi+DkrRAUpOkpvb29hTKTs+RZO2Tq6fVM3RQXdcd58+vXlFmZt2Up0Dp7IRBV5cyzQf+vcN011URMQu4Gbhb0tWd7RgRiyOiISIa6uvre1dxylZv3klh17snnu4CWLasOgWZmZUhT4GyCTiz5P1koLWLvrfRYborIlqTn1spPvrlsgrUWFGNzW3UDRDXe+0TM6tBeQqU54Fpks6WNJhiaCzp2EnSacA1wHdK2oZLGnn0d4oPqHy5KlWnaEVLgcvPPoNRwwZnXYqZWdlyc5VXRBySdA/wJFAHPBQRzZLuSrbfn3T9ONAYEXtLdh8HPJ7cVT4Q+GZEPFG96ntvQ/seXt26h9+6/KyTd/ZNjWaWQ7kJFICIWA4s79B2f4f3D1N89Etp2wZgZoXLq6iyHga5eLEfv2JmuZOnKa9+rbGlwIUTT2Xy6cNO3vnOOytfkJlZmRwoOdC++11WvfkWc2Z47RMzq10OlBx4ak2BCK99Yma1zYGSA40tBSaffgrTJ4zs3g5Ljrv4zcwscw6UjO199xA/Wb+NOTPGd772SWdmz65sUWZmPeBAydjTr7Rz4NCRk98dX2pSp0+kMTPLlAMlY40tBUYNG0TDFK99Yma1zYGSoYOHj/DUmgI3XDCOgXX+ozCz2ua/xTL089d3sGv/SdY+6cwdd1SmIDOzXnCgZGhFS4EhAwfw4WknWPukM4sXV6YgM7NecKBkJCJobG7jw9PqGTa4zCfg+CovM8shB0pGmlt30bpzf/nTXQCrVqVfkJlZLzlQMtLY3MYAwQ1e+8TM+ggHSkYaWwo0TDmD0SOGlL/zhAnpF2Rm1ksOlAy8uX0fa9t292y6C6C1q4Uszcyy40DJQGNLG9CLh0EuWpReMWZmKclVoEiaK2mdpPWS7u1k+7WSdkp6IXl9obv75kljS4ELxo9kyujhPRvgi19MtyAzsxTkZsVGSXXAfcBNwCbgeUlLIqKlQ9dnImJeD/fN3I69B2h6Ywd3X3de1qWYmaUqT0colwHrI2JDRBwAHgVurcK+VfXUmgJHAi+mZWZ9Tp4CZRKwseT9pqStow9JelHS9yRdWOa+SFogqUlSU3t7exp1l6WxpcCE04Zy0aRTez5IU1N6BZmZpSRPgdLZYiDR4f0qYEpEzAT+Hvi3MvYtNkYsjoiGiGior6/vaa098s6BwzzzajtzZozr/tonZmY1Ik+Bsgk4s+T9ZOCY62MjYldE7El+Xw4MkjSmO/vmwTOvtrP/4BFu6u10V0NDOgWZmaUoT4HyPDBN0tmSBgO3AcesdStpvJJ/2ku6jGL927uzbx40thQYOXQgl59zRtalmJmlLjdXeUXEIUn3AE8CdcBDEdEs6a5k+/3ArwO/J+kQ8A5wW0QE0Om+mXyRLhx6b+2TsQzy2idm1gflJlDgvWms5R3a7i/5/cvAl7u7b540/fIt3tp3sPfTXQALF/Z+DDOzlPmfylWyoqXA4LoBXHN+ChcC+E55M8shB0oVRASNLW1cdd5oRgxJ4aBw4sTej2FmljIHShWsbdvNxh3vMOfClG5m3LIlnXHMzFLkQKmCxuYCEtww3WufmFnf5UCpghVr2rj0zFGMHTk0nQFnzUpnHDOzFDlQKmzz2+/w8uZd6U13Aaxcmd5YZmYpcaBU2Irm4tonc3q69klnFixIbywzs5Q4UCpsxZoC59YP55z6EekN+uCD6Y1lZpYSB0oF7dx3kOc27Eh3usvMLKccKBX0g3UFDh+JdKe7zMxyyoFSQY3NBcaOHMLMyaPSHXjz5nTHMzNLgQOlQvYfPMyPX2nnxhnjGDAg5bVPfJWXmeWQA6VCnn1tG/sOHK7MdNctt6Q/pplZLzlQKqSxucCIIQP50Lmjsy7FzKwqHCgVcPhI8P01Ba49v54hA+uyLsfMrCocKBXwwsa32LbnADdV6uquBx6ozLhmZr2Qq0CRNFfSOknrJd3byfbflrQ6eT0raWbJtjckvSTpBUlN1a38WI3NBQbViesuqNDDIH2nvJnlUG5WbJRUB9wH3ARsAp6XtCQiWkq6vQ5cExFvSboZWAxcXrL9uojYVrWiO1Fc+6TAFeeM5tShgyrzIRJEVGZsM7MeytMRymXA+ojYEBEHgEeBW0s7RMSzEfFW8vY5YHKVazyp19r38Pq2vb473sz6nTwFyiRgY8n7TUlbV34X+F7J+wAaJa2U1OWckKQFkpokNbW3t/eq4M482VwA4KbpvjvezPqX3Ex5AZ3d/dfpvI6k6ygGyq+WNF8VEa2SxgIrJK2NiKePGzBiMcWpMhoaGlKfN2psKTBz8mmMPy2ltU86M29e5cY2M+uhPB2hbALOLHk/GWjt2EnSxcBXgVsjYvvR9ohoTX5uBR6nOIVWVYVd+3lx49uVn+5aurSy45uZ9UCeAuV5YJqksyUNBm4DlpR2kHQW8BjwiYh4paR9uKSRR38H5gAvV63yxIqW4nRXxR8GOX9+Zcc3M+uB3Ex5RcQhSfcATwJ1wEMR0SzprmT7/cAXgNHAP0gCOBQRDcA44PGkbSDwzYh4otrf4cnmNqaOHsZ5Y1Nc+6Qzy5ZVdnwzsx7ITaAARMRyYHmHtvtLfv808OlO9tsAzOzYXk3b97zLs69tZ8HV55AEm5lZv5KnKa+a9kRzG4ePBPMvnph1KWZmmXCgpGTpi62cUz+c6RNGVv7DfFOjmeWQAyUFW3ft52ev72DexROrM921eHHlP8PMrEwOlBQsf2kLETD/4gnV+cA776zO55iZlcGBkoKlq7dwwfiRTBtXhekuM7OccqD0Uuvb77Dyl28xr1pHJ2ZmOeVA6aXvrt4CwLxqXt21ZMnJ+5iZVZkDpZeWrm7lg5NOY+qY4dX70Nmzq/dZZmbd5EDphV9u38vqTTurP9016UQPYTYzy4YDpReWJdNdH/P5EzMzB0pvLH2xlVlnjWLy6cOyLsXMLHMOlB5av3UPa9t2V/dk/FF33FH9zzQzOwkHSg8tW92KlNF0l++UN7MccqD0QESw9MVWLpt6BuNOreDKjF3xVV5mlkMOlB5Y27ab19r3Mm9mRk8WXrUqm881MzuBXAWKpLmS1klaL+neTrZL0t8l21dLmtXdfdO0bHUrAwQ3X1ThpX7NzGpIbgJFUh1wH3AzMAO4XdKMDt1uBqYlrwXAV8rYNxXF6a4tXHXeGMaMGFKJjzi5Cb5M2czyJ08rNl4GrE9WX0TSo8CtQEtJn1uBb0REAM9JGiVpAjC1G/um4p2Dh7ny3NFced6YtIfuvtbW7D7bzKwLuTlCASYBG0veb0rautOnO/umYtjggfzlf7iYW7I6fwKwaFF2n21m1oU8BUpnK1N1XJqwqz7d2bc4gLRAUpOkpvb29jJLzIkvfjHrCszMjpOnQNkEnFnyfjLQcW6nqz7d2ReAiFgcEQ0R0VBfX9/ros3MrChPgfI8ME3S2ZIGA7cBHZ/TvgT4neRqryuAnRGxpZv7mplZBeXmpHxEHJJ0D/AkUAc8FBHNku5Ktt8PLAc+CqwH9gGfOtG+GXyN6mhqyroCM7Pj5CZQACJiOcXQKG27v+T3AO7u7r5mZlY9eZrysu5qaMi6AjOz4zhQzMwsFQ4UMzNLhYqnJfonSe3AL3u4+xhgW4rl1Ap/7/7F37t/6c73nhIRnd5z0a8DpTckNUVEvzuZ4e/dv/h79y+9/d6e8jIzs1Q4UMzMLBUOlJ7rr+vw+nv3L/7e/UuvvrfPoZiZWSp8hGJmZqlwoJiZWSocKGWq5tr1eSHpTEk/lLRGUrOkz2VdUzVJqpP0C0nLsq6lmpIVUb8laW3yZ/+hrGuqBkl/mPx3/rKkRyQNzbqmSpD0kKStkl4uaTtD0gpJryY/Ty9nTAdKGaq5dn3OHAL+KCKmA1cAd/eT733U54A1WReRgf8DPBERFwAz6Qf/G0iaBHwWaIiIiyg+vfy2bKuqmIeBuR3a7gWeiohpwFPJ+25zoJTnvXXvI+IAcHTt+j4tIrZExKrk990U/2KpyBLLeSNpMvAx4KtZ11JNkk4Frga+BhARByLi7UyLqp6BwCmSBgLD6GKxvloXEU8DOzo03wp8Pfn968CvlTOmA6U8VVu7Pq8kTQUuBX6WcSnV8rfAHwNHMq6j2s4B2oF/TKb7vippeNZFVVpEbAb+BngT2EJxEb/GbKuqqnHJooUkP8eWs7MDpTzdXru+L5I0Avg28AcRsSvreipN0jxga0SszLqWDAwEZgFfiYhLgb2UOf1Ri5JzBrcCZwMTgeGS/lO2VdUOB0p5ur12fV8jaRDFMPnniHgs63qq5CrgFklvUJzevF7S/822pKrZBGyKiKNHot+iGDB93Y3A6xHRHhEHgceAKzOuqZoKkiYAJD+3lrOzA6U8/XLtekmiOJe+JiK+lHU91RIRfxoRkyNiKsU/6x9ERL/412pEtAEbJZ2fNN0AtGRYUrW8CVwhaVjy3/0N9IOLEUosAT6Z/P5J4Dvl7JyrJYDzrt+tXf++q4BPAC9JeiFp+7Nk2WXruz4D/HPyj6cNwKcyrqfiIuJnkr4FrKJ4deMv6KOPYZH0CHAtMEbSJmAh8JfAv0j6XYrh+htljelHr5iZWRo85WVmZqlwoJiZWSocKGZmlgoHipmZpcKBYmZmqXCgmJlZKhwoZmaWCgeKWU5I+o6k7Z20ny4pJPWbpxRYbXKgmOXHJRTvzO7o6DO0OttmlhsOFLMcSJ5yexbFR350dDRQOttmlhsOFLN8uCT52VWgvAOsrVo1Zj3gQDHLh0uTn10FyuqIOFzFeszK5kAxy4dLgN3Aq6WNyVK85+HpLqsBDhSzfLgEWBvHP/77Zor/P/UJecs9B4pZxiQNAaYD9ZLqStpHAIuStw4Uyz0vsGWWvYso/n+xHlgu6bvAGIor5o1I+twmqRARGzOq0eykfIRilr1Lkp+fBIYD/5PiCplforhq4j7gemBHFsWZdZdXbDTLmKS/B+4CRkTEu1nXY9ZTPkIxy96lwCsOE6t1DhSzDEkScDHwUta1mPWWA8UsW+cCI4HVWRdi1ls+h2JmZqnwEYqZmaXCgWJmZqlwoJiZWSocKGZmlgoHipmZpcKBYmZmqXCgmJlZKv4/J3JWVAHZhZIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate\n",
    "import scipy.optimize as opt\n",
    "\n",
    "# Fermi Dirac distribution\n",
    "def FD(E,mu,kT): return 1.0/(np.exp((E-mu)/kT)+1.0)\n",
    "\n",
    "# Defining integral function G(mu), for given Emin, Emax boundaries and kT = 1/40\n",
    "def G(mu):\n",
    "    kT = 1.0/40.0\n",
    "    Emin = 0.0\n",
    "    Emax = 2.0\n",
    "    G = integrate.quad(FD,Emin,Emax,args=(mu,kT))\n",
    "    return G[0]\n",
    "\n",
    "\n",
    "# Find mu as the root of G=1, via bisection method\n",
    "a = 0.1\n",
    "b = 3.0\n",
    "def func(mu): return G(mu) - 1.0\n",
    "root = opt.root_scalar(func, method='bisect', bracket=[a,b])\n",
    "print('root =', root.root)\n",
    "\n",
    "# Display the solution\n",
    "muvec = [ ]\n",
    "Gvec = [ ]\n",
    "for mu in np.arange(0.,10.,0.2):\n",
    "    muvec.append(mu)\n",
    "    y = G(mu)\n",
    "    Gvec.append(y)\n",
    "plt.plot(muvec,Gvec)\n",
    "plt.axhline(y=1.0,linestyle='dashed',color = 'red', lw = 1)\n",
    "plt.axvline(x=root.root,ls='dashed',color='r',lw=1)\n",
    "plt.xlabel(r'$\\mu$',fontsize=18)\n",
    "plt.ylabel(r'$G(\\mu)$',fontsize=18)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
