statmek.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "from sympy import *\n",
    "from scipy.integrate import quad\n",
    "plt.rcParams.update({'font.size': 16, 'figure.figsize': [8.0, 6.0]})\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N,V,T,Q,L,q,k,h,m,beta = symbols(\"N V T Q Lambda q k h m beta\", positive=True, real=True)\n",
    "E,S,F,A = symbols(\"E S G A\", real=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ideal monatomic gas\n",
    "\n",
    "Here we show how to use `sympy` to symbolically manage expressions, take the derivative etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Partition function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = h/sqrt(2*pi*m/beta) # de Broglie wavelength\n",
    "q = V/L**3              # translational partition function\n",
    "Q = q**N / factorial(N) # molecular part. func.; independent particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADgAAAArCAYAAAAtxEsrAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEZklEQVRoBe2a7VEbMRCGzxkKcJIOTAd8dAAdAB0EOgiTX/Y/BjqAVJCBDnAHfHSASyDuwHkfoRWSTgdOxnfBHu+M0Gq10u27u9LpZHqz2ayCRqPRQNWJa1RVXzXtc8nHXpZUkn+XYF9lz3dsS/YYK6lN37nKlspU5V7lUHL4SjXPeYAXTVScXDX6UCzDHvQvNe6CTj/+WuyOyrPam8hj6gHQKwLGACI7kCKDMegmHhTz6rtVGyA34g/jPuMlfxKPAwyA61IbJ+2qfIv7xP+WrGaw5Nhzq/pKdSC1cdK96mC/dX7yzLHqYykAysgi98MEeS19gF2qoHugNl4uEeATcF5pV/IQUWR+DiJlz/eqrvqlv0Q1JyKOo2tkAEktDAhGNBiUT0CKYghpCNU8qHlwwp3rjf54ICWj0IdKfTiQNA/k56ddcki1QY+U6PwMbySZRZMINVFfejhlrBrPkgmnmTJOOMtkNAcak6Sa10EfKhn86J/3ovHyF/1pQe56HcBYG17KeJGonIgvGYFaTuizAZDq8RhzQqIvnRIAdHj2RP0hm2xgwxj0w1zSITCs677KqaWoeAdsSwosfFKNtE3SQe1A0mPikEZqG6g8gjVDwyQZozlIQQwLBmcqSVP66Ib1p3acdUfqO8oBkgIXKuyGLOiHaJCaCdn6i4WAJPUAX/k6OCFWbODdOPXNO8b0WSL2THZ8gLNkxglACQJpAIp4/1o8A3IqpZ5tNhZFnNCYBfmEar+1/grqTh8bAcdrBZsr1QSK19LEARRDatrLFR0jM848ZXJqJk6ICSUgvfbED1SXnJCMyRo8p7j+Mj1rOn01tlXYL2oYLIK8KEnHUqRsslBLj4lrW79XiKNYc0KYJGPMIRL/zfrDiWcay57BIQAMCUgDiCGll/GO5FD+0NL6c4p6ALr2yph3LTHWsmTeMWwikNlm2QboSnawWVYGkDWTTCwFdiQiSuhDJMQj49QSZGrn5KIoHXt43l9q2yFh3jE4ufRexLnQV/64syiMjMGD8VkST4TDtvoB9lMFPXheI6SHW9jiE5KcM6NtGkmfNdSPI3gO81kEAYjz7tTvDtXia8T8EvLeDc8XT9Q4cHP2vVJ7GgBKsJJkKbqS4AC1BrjsoV1HcB3BD+6BdYp+8AC9a15vOBy+3Bu+q7qcCuuTzHLG7dXq9Sbz6ovl5FY+ghtdxEXfZQM9xy6iuCV4pi0535StUuu7qAfHh2n4+BXPhy4fp/vi5/2C/ydHdJGigLHrCGekQBFNvtq5KGqVugDIVcSTQHEtEROR41qR9G2NugAIkLfuOnPgCwXb+hpsslaR42KIa/5ek84i5F1EsGanQHE5S2razlrTWZTgvwCU8WwuXDQ3XgsuCmDnKSpQ/KBKaobXxqLAlObpNIICxS/AX7oCB+DOAAoUPwVsqg635+KJJGuxNeoEoECwqfAfFfmmAmiOba1R62vQR4jfEUpHMn5HrP3zziLRdnHYBhxpyPrLqfXD9h+9DaG8jnme/QAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\frac{3 N T k}{2}$"
      ],
      "text/plain": [
       "3⋅N⋅T⋅k\n",
       "───────\n",
       "   2   "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E = -diff(ln(Q), beta)\n",
    "E.subs(beta, 1/(k*T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Free energy (experimental section)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAA/CAYAAAA8PsvmAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAXe0lEQVR4Ae2d67XdttGGt7RUwLHSgd2BbFUgp4PEqsBKB/HyP//TSjpIVIESdeC4AtvqIPkq+JTTwcn7cHNokpsgSAIEL2ewFjcvAGYGLwYzuJH7ycPDw+WxhB9++OFvKuu/df7rYymzl3N7BFzv5teBMLtTrp90vNL1/XwKnsMR2A6Bufr75LE4YgHzZ1XLa52/3K56nPNjQ8D17nIRBl+r3r/TwflPuv+76UEd95f6/r3um06yrmmzv9f595bez47AURCYo79Pj1KoFDkFCAaAxv7HFDqe1xGYg4Dr3RUt4fAvXf2o44MOc7pVZB33Xjd/1HXjhIm0e507eaqM/uMI7ByBOfp7ekcsMJji+qcOeuL/2XnduXgnQcD17qYiv9CTb3XcCZs3vdjfjbRNOs9vFP+HXh6/dQSOgMAk/T29I1ZN0Zv+jxpyMx12hNpzGQ+PgOtdtwrv1QZZ66UdMk09KdR5cODvdE2n2oMjcBgEpurvqR2xQHihGqP3TUP24AgUQcD1rgtzjcfP9VOc8Od6Vo1we3HdjPWd0jCl/UmHT1EPIuQP94zAFP09tSNW5bzT8UFAfNxzRblsp0Mgm95Jd3/UQYeyeBBfppHhnzoSZY8G68QX0bJRsTnVJi5SQBw4U9SfR9J5tCOwRwRG9fe0jlgNlh43BmzyNNgea89lOhYCOfVOtNjb8E+dN+lIii9Ok/bzq65TnPEXNS2rTJwwo2KcMOvD8BkNSsOoGBx4BdGDI3AoBGL6e1pHrFqisf9LAPgGrUOp7OGFzaJ30luccPLeBtF5oYN35xc5UuXD+VXOOKFmOo5WNGmTNzuoJ9B/qzRfK/8mMwQT5PMkjsAYAkH9PaUjrhvq50LEe89jauFxWRHIpXeiw74GRoyLZnOUjyllRtLoP9PktIXFQXRwmizx0DmYFZQHp2nrw+28lC0U107XXNdy4NT/1Dz0C0fgIAiM6e8pHbHq5Xvqpi74QarJxTwBAsl6J53FaTKqXvzOu2iwQ5n3cnFY73PgKlqV49S5/+pRkLzS8kEOnDevDnZGsbq3UXG1dhwkchvxDz2aLMNtdn/iCGyKwKD+ns4Rq4EzBcf6sL+utKm+PS7mGfWOD1+8rR3V3kDEGf+lLmtUNqX7qw7Wh/k61s06t57RWehMW0eJXh37RfncGU8Ay5PsDoFqVqmvv6dzxIL9mxp6n5benQ6eWqBkvasb53OdO1+Y2gtqkospakayjNg3CZKBEbRPT2+CvjNNRSCkv2d0xNWUngp80wNPBdHzOwIjCOTQOxzc3mdyGBVv/RoR03tsQmMa34MjcDQEbvT3jI546nuJR6s8l3ffCCTpnZwK66ksq7Czcreh7tHTyd1sVCzeTN8TwNyDI3A0BG7091SOWEbCGqaPho+mmgeWN5PesbGKnclz10y3QI7Owh8kKx2HLYK1b/9Xpi3Qd56pCNzo76kcsdCxnabW40gFzPM7AlMQSNI7OTR2FDPNuuvRcAsI2+m8yYYp4cU6NR0W63i3RPNLR2DfCAzp79kc8Vd1Ffyy76pw6U6GQKre8doTrxxZT3nX8EhOnCDOmFH8VoE2zvvSW43Ktyq38z0HAh39PZsjZmTB14gwFB4cgVIIpOpd0vpyqUL2+PAaBh8dYSS/RbBZLx8Vb4G+80xFoKO/p3HELYPAtJUHR6AIAql6p/w4cUZ11jCLyJ2BCTs/CVs5Qmvnvk58rQf/PRYCHf09jSNWHWDQCFbA690j/205iixI5KaXRahtiaTqnTkyW3fNXZrf1QSf5yQsPWDWiWMrR2jtfKsReU44s9Eq0T5L8MgGyEJCsTLG4iew7ejvmRyxNch/TwDhUSSRsvBKjDmKXGVmOhK6Hq4IpOodjoz1YWuYWXAVPb41zfSxbaiq7vXM7nPwYZ0rt35NlcvwMvyn5jttOtXtGu29g1cJHh2G293E7FwsPiZ5R3+fxVLH4lUxTKv9WqeDOL1kgjXQ9jMaDen/pnzV14Pq/BgMNrx80v0XOi8Jls8KOIuG+KLEGEUboXypZ53NM7onjvcnKRvlxBB1PtOnNMS900FZ+fcn21Gr23JBfPnMJ38xV+EMZ11Txtc6rG4YhVEXH3Rugu6Jr9YAdQZPylFtzNGZaz6m8EZH8Y9PwFvybI6vZLCQpHcigt6jR1mDcCqhd7QPay9Z5Y8RU/novJBsU0e8F32UHJ32XsvFJkDay3OAUmjPumCH0T3C35WeD7WMhj6PfmLFZ7Ghfbpr3dcYYeew9ea3Kna6H7VzsfiYzPDTQbJKf588PDzE8ozGixjgv9TxLcQtsa7/q+sbx6rnFJw/G+8Ycd3jzH/RedFOTOVjjQ2jcONA9WxyaNHhnc5BY6bnjLpvKq/NRGkoDxvHBmm00+a+Fk8a2U86f9mnXcdRNx+H4i294sCSj/UPyq/nlO+Vzk2dW94S55r/Jvi2yyc5Fuud8lJP1AWGcJHet2UpfS2ZGV3zKVm+Jd028kVEEU/aIYbsM11voodWUPHfa3unHaOjg/ZMclN/fFZ1sJ23yhe0KZbGzqJlbWKQJ+mUJmpDjV7Os/hSjnc6PumgI0JHJag/Sj9q52Lxoh0MNQaV/j4Nppoe8VIE+6NCiFPgocb5Xs+HRq0AQgUuDdbrW9wgVQ6UFsVEbj5YQDmGAgoW40NFbxUYtVOOm1DLjeyhslmeoBOuE0A/6etK4K1j6VTplvgaRpxT9M7q4KjLKTaSp+1uEawNWh1sIYPx3FIfg+1dwg3ZWpP5ovZHB/CueRC+GOPR5KJN6yaXDW3oti/goWOR3VA+RqL4K8r9vk03cB2zc7H4ANnqcaO/T8dSxeJUGAzJkPOkMghDceSxBlwlEh1LP+S4qzQTfqBLSGkQTE0jgzmYm1FKLevPMNpx+EZydmYcerLSOIONT3mZ5Xjby9O5renDJ0ink2H4hrwp+Yepln2aoneWd9RYli3OLG4mt03Pz8qcIbG1dcMxA8lDkoi191ihrB7H0k3lUcKGFrMbMTsXix8DVHGN/j6NJIxFs2A9ZPCpDMKQY2VK9P4a3fySnp5K/3mTYMJFZdBTaZBfB3KjnEO9LlO0CSLlTyLZRo2O4lkrijWsKn6IVv2MteWPE6SHzjcT0h0myRAmEeFT9M5GcrH6ioiwTbSwsvZq5SgtiPGv6qA08xL8Yvqo+CntPSbq0ICpyTOTxx16oWO3NrQp2PSLmJ2LxYc4Nfr7LJRiyvMa7KGkjHAHP6wRyEP6xmkrDcr1UgcN7Dvdm8C6LRaqqRjxfqOj3dmoFG2pFKJFmaBt05GMJtjR2pQf2roHE0bkVDLOl8ZCXnBljT00CpnSUTDDD1271mUV+L/Z0fUiS6gzMsOvjU8ruvylZF8b35yFQlbC/fV02F8rx2ELsJbgBfRxSnsfLZ5k7GzWHEi8lMcqNnRAvrUfxexcLD4qX+qI+IaBKhXjTsPsOJabhPWDWlFZY6p6ZbrHCROYe2e0VWTEJb44vqZnqHtzLv3dhIuNpmhSThb/cXb8aToHzpbOBtPBVdC14cEGOOJwjCg1PU0azc2UeZXx+vOVTubkW487lxZPXTVBtG3zTfMscgGdDo1I+lWjJX8JfHOWwd7xtSmqnLRL0aI9PC/FrMfHcNuKf0+c7m0hfZzS3ruCzb+bxEPlXd2Gzhc9S46YnYvFh4Ro9De7IxZHKoPQOLXrbfDX0rNdvLrWGWeDM2e0FnXoSk/a1DDU68MZM/1ucnGeWq4hedgxzkav/igUZ49zxpEQvtfRn1EAh6pToHRjmICFVbAuB4Pxb0bVoolD/SJCu08MPrtxxJKlBL4NBsIqVe9S8zeybHiBDpyhHGtAWEIfp7T31LJN5VHChqaWZUn+mJ2LxUd5PquNyU9KOacxsevsY4A6lUEYcxbXFNdf0t/rwMn9Yo6gpn/z+s01S/AXOkvD0JQzo1BGiThAyoOso5uYFD8YVB6cLE7rZqMXZdVBvtc6QrgSPyUwOojhYJvlzPFDd86UtMmBQ4/qjcrG7EbVmbGM9RlZL4ofGuGDydQpcmiUwrdXhOo2hvdQHp5Z+ZfmD9Et/bwqR2mm4me4RXXQZJOekDanvTPSnXNBfQR7w6EjQ8abqTyy2lBhuLrdmIhRzM7F4kNsrN7ucMTczHV4IcI8x+D2R3PR9EqADLybywcmUp3RGL9QnIHSxEsOyoED/lpnnOiQojXpIxc2crzh08r3or5GAVkHbvMjjo5BchDdex3QqWTS9dwpaZNhUgMV/SFHe9FzliGYcWg+OmKEF5x3g+8C2WdlEV5pL/9HuIn+k0iSw0arbLntXQiLR6OPLQBubJvwXmxDlbeE3WiJH7yM2blYfJCwRTyzixxnAYfy3emwD8KPklV60pKHETZTtjjxX3VO+ijHKNOByJrvzUi1TorzQy5GxTeKVqeZcrLpYMocCpaGESvT4u8kG8/AiBGrrVvrNhiYJhnjYRkruqJJWqakp9C2vHYmL/z2EAy7sbJbmhR8c5a1wo460HE/lbDS7s1R7kUHpkJYIp3p2tr6OLW9p5Q5ykM6iY1c24amlCElb8zOxeKjvLM6YnGjMghT11FtI5ZNY9uUKY6Hqck/64iOlpTGRnhjSl8JFvgJTjmLNmvXNCpGjTbtHiATfiwa9toWND60UyrOcGNNicD9zU7qKib+g6zgFwuWDj5Ly0VP0AxOjN+q8QXxbcohnql6d18TS+5RN0KVv9hSB6y9G47lSx/gWFAfrR0HJMnyeAqP1W1olpIsIxLT8Vh8iGujv7kdsU0lmGMNCWDPqTxzUPaMsxl321Xajst6rQYDGHxFq9oIFSDOqJgp86nlCpC5vFIEn57sT79Dnx3URh/Dwgj4W53NyHzSvV3rcTB8VAyvfsWCYZwy3c1yAvz2Ekrgm7Os/18TQwePHEyXjlyGNWQvoY+x9m66Zecl5RzlIbsE7VI2dIn8oTzmX2Id4Zidi8WH+DfPkx2xKgFDzgiMynhRU2aUda/rn3UeG9GS522d50IeHThEPq/IlvAmztLkOos+vN/pYPTJ1CCvFb3VuTNa1bOLnvEt4Ekbh5QODNj1zJZ/8jLK5RUk1krodFBp7XejwY57c8IXrnXocfWqE+cq6BkGDyc+hul7pbGR9TXj8C/lpVwN3+Fko0/BzjpfowlzRUrerfHNVRTo0EaKBuGHvlFvn3RwTacNvV/aoaIdQWuLgAElbMX/sgN9HGzvLbmoYwL7XMwuzK3vEI9VbOhV3PV+WzjQDgj4K2zrzX8gVLHX9jJm55bawUZ/k//0oRZ085OA5OP5ONS9rZ/NxkZlYGofh1t1CnRG4ak0nBAdn1FnrPR0YsZ2tis6LYgHDXzswyJRBqKRc7NWlJ8lEN8kfI0OZ9FarHd1+TGOWf40QfQqg6BzsNOoOPjRKaxmgHTm/fXvdf5M51lBedBLys9sztiM0iy6UxOLJ2VBh6o9JlPz7S2dypGkj8pfor2vzmNqvai8xeyGeI3auVj8WJmUt9HfZ2MJDxZHrxhHzHF/MNkbcSV71UvTuRmZ1+WhTBhQHDKGdmxUjLOmBzfWi1N0UsDwwiclUKaidZUJ33aZU/TOpnRt1NKmu+Sa+sZIjYW3vUim55aOKE1um2LvkV799q7mYDiuzjA3g0z6WKK9l+AxFd6SdiNm52LxY2Vq9PfpWKqDxVljNONwMPEbcX/RFa/0MPodChhbXm8KBuW1D5GsgoXoQzf0nfGgXP0I0WEKfslu7T6pOffJ+PaYLdY7lf1jTYvlihzBOnFBZwzPFl94knZph8r0yzDIUYY5NIx/0c7cHAEnpE3Wx7oN0R4Njwls5yUpwWOqRJKliN2o8QzauVj8hPJYfd2faURsxoAR42GDKvdeB4aZ6cLXOttog5ELPSg+e2kGXLfBwKi5mvYMplgeQUdgzdH2cskiOTPia5xS9Y78XxmxpWeVCyeMUeeM3nzQEQxKz5Q06VjmWNoZMkMyRR+DsiREVG1d8lsdJJDaJqtkP0J7N3DWtCnGY0/nmJ2LxcfK0ujvmRwxaxgEMw7XuwP+0jglNlMei0PdwNn0NukVsKmMoKe00D208VMZkvBt4ZWqdzgxnGdqwEjSOUKe4IjYmKj+WNfFAbNR5cK9xc04vySt8m6lC3RMjzwarqAWfrtt75WA9Q9y6shuU9o89nKtco7auVj8xHI0+nsmR2wNsvl+8kQwTptMyoKBXGJgg5iIZlZ6QUbHiUjVOzbq8OpH6t6G59S3DkbC0GOX7OiOeMVjWOmQ8BEdpvvmjmzpQIzyUPwqQbJah/vTKgwOSFSYZG/vfRhK8Ojz3OJe5Ry1c7H4mMzK39HfMzli65VbAWNYeLwjkAOBVL1jOpmAUxudTq5SDfyoUbOf4Oc66r3OjIgZIXecpNLRA/8/Ha90bU7XOhJMj9szXY6H2pBAj47EFsHa+WSZtxDSeToCAQQ6+vs0kOiIj82gWQGPWAaX+XgIJOld7RBxhqzXLg3krZy46Jkz/6ZPTHHwYQRpnQeS0AHg+T+4mRHIRzB+17tyv9bOrQNSjrNzcgTSEejo72kccW1kMDCMDjw4AkUQyKR3OEFzbEvkZmdn27kyEmaqe6gtMFJmIyB7B1gH4+t2fNsdZzwnkI+p8DbfOflT09pO886oP5Wo53cECiHQ0d9nhZiWYkOjfCPj0DdMpfg7n8eJQKresbsdvY2u6/bhRdf1rD8qhB6OnZFyZ+pW6bnvPNP9ksD0N+vLWwWm0i91ebaSwfk6AksR6OjvaUbENRp8tpGQMrq4UvBfR2A6Akl6J2eCI2dkyWh1bsAh9qeHbZqZuOxB8vIHKISlrz1dc6f9MtrfajSeJrnndgSuM7eN/p7NEds0lQ37vcIdgRII5NA7PqrBqPhupsAvladp0OTVPdPMjHrX+sgDI2H+tnTudLaypQfxtSn3fgcknbhTcARWRmBIf0/liFVADBLGwUfEKyuTk/8NgRx6JxqMLtHd73+jPH6lPExLd5xwKwe7pwlZP7winrQt+G45LW3t28pIOT04AkdB4EZ/T+WI61rAoDESuDtKrbicp0Agh97h3NhENVV3GUXbl9f6INpoMff0tI2GQx2Avhxr3LNRjHegc6x1ryGf03QExhC40d8zOmI+O0awdazrnf86AusikKx3ciw4c5zL6KhY6V7oYF0aJ8v/Vnfe5dU9DpoNWwQ6pXywgx3SSUE06MlzfJtEKCGzZKCTggy2HJBAzbM6AmURCOnvs7JirM9NBeWVCozZax2jX0dZXxrn8FgQyKh3OLmfRI9vQA+OOmv9Du6DUPxa08Z0NvjLQabQtwr2fvTbrQRwvo5AAgKD+nvGETEYYTAYNbCW5cERKIVAst7VThYnA63dBMnFKJt/brIp761kY83bOttbyeB8HYGlCAzq71kdsb2+QaE9OAKlEMiid3J2zOSwBorz2zxIjmpKWuclr1dlk1/86VizY3oXuGQrmBN6FAiM6e8pHbEKzNQZ622+TvwoVHwfhcypd6KF0+MDH7k3W80CS/zN8b2alXGdxFXHWjLRtj04AkdDIKi/p3TEde2wTsZn/pI3qRyttl3eTRHIpnfSXdaB+ds5e2+2aMHEl41R73TwJxFbrgtfalloy2utfxfF1pk9LgRi+vvk4eHhtIio8DRcvqv72WkL6QXbHQKud/mrRJiyZs4Mgf/NaX54neLKCMT098wj4osKz1rbJ519VLyyojn53xBwvfsNixxXwvNz0WGZyfd85ADUaRRFYIr+ntoR12jTeHnXksbswREohYDrXT6k2ZzFJzX93eF8mDqlcghE9ff0jrhuvIyM7QMH5eB3To8WAde7PFUvHBkJb/oRkTwlcSqPEYGp+nvqNeJ2xQsQvj70o844ZQ+OQBEEXO+WwyzsmMXiC2JsFvPPWS6H0nNugMAc/T39iLiFP6+DbLYDtSWHXz4uBFzvltc3s1jfuRNeDqDn3BSByfr7aByxGjOvX/CxbQ+OQDEEXO+WQS3ceHUKJ+zvDC+D0HNtiMBc/f0feuC6wTs8OYkAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle - T k \\left(N \\log{\\left(V \\right)} + \\log{\\left(\\left(\\frac{1}{\\Lambda^{3}}\\right)^{N} \\right)} - \\log{\\left(\\Gamma\\left(N + 1\\right) \\right)}\\right)$"
      ],
      "text/plain": [
       "     ⎛              ⎛    N⎞                ⎞\n",
       "     ⎜              ⎜⎛1 ⎞ ⎟                ⎟\n",
       "-T⋅k⋅⎜N⋅log(V) + log⎜⎜──⎟ ⎟ - log(Γ(N + 1))⎟\n",
       "     ⎜              ⎜⎜ 3⎟ ⎟                ⎟\n",
       "     ⎝              ⎝⎝Λ ⎠ ⎠                ⎠"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = symbols('Lambda')\n",
    "Q = (V/L**3)**N / factorial(N)\n",
    "F = -k*T*ln(Q)\n",
    "F = simplify(F)\n",
    "F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMwAAABACAYAAABWUJhRAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAQKElEQVR4Ae2d65EVtxLHDxQBLOsMIAMDEYAz4BEBkAEuvvGNsjPwJYLFZGATAdgZmBvBXW8G3P9Pq1bNQ9JoRjPnnN2VqsSMpFar1VK/pDnLre/fv+9aahy4yhx49+7dE9H/szLP1yr/x+bj237x5TOVf7W2Jc9bTWCWsG1eHy3Sb+rxT+1izRv1ZkGLt28040fKT/R+tzt73/ZJz2/det5Vd6LHZ+XHer+gLpdu5xpbWz0HtAgs5EM9qzRbPSXXHsN9zfCl8ol4/Wow2x9UNxIWYFSPkJwp/055KjWBmeJQRbsWAxcBd+BZBZrWtYwDF37z447hnhUn9XPKTE9z3ZJ9m8AkWVPXIOZj6tFa+NRR7VY3wnH3Zv7Kr5T/2ppSjfGjxvjix0FY7qnuKeVBmweJPlBq0Ov6RSFU2QQmxZn6erTVNy1ACEDrUaYxaBx89z+Uvyv3XBLf9peeZFzEWQl84CjtJFg2MDSgNMizkvoT881J0PYnHdQXFwuem7UIbbSnku+HS/dB70mam8CkOFhRL4bbhmEB9pI0JhvmD+VPyrZZ3Ni+7UyFZ3qfFUsJHo37wONw+Kb+EezffpyllvUX9Wcupem+4BEUS8wfK4OwEL902wxm9BQcvDtX7vGvC9gEpsuN9d4/CBWnMn+vh7II06LAN4VZ9KNp3+r5OgWzRb3GQ9B+17PUGvYEwvcfKY5CWnHpsKj3YvB3YpWtbjkHxGg0MhZmMtD3G+KFh2dQrMRvqmexQ1IZfMRDLCKb6U/VxTaxC3zVZoFvrTvImHPdI3WpT8xBmaN4nj2B6GJXG7yx+KXbxMb/R/l9t3LqXfhM0THvn4bwzcIMOVJfxpyzoSfdEcHgHj32Q+LG/MSCDUlQHZYKAWExcT9GwqK67sZZGviGoYUP4eQ4vFboAs4FL2xarHU0iTYsEELNwQrzD0ll+A8vXWwTGspeEDJiwh5OujaBKWNgEZRnMButWCurD9qTTL9cYlPkrFYIbj3O2YHvYHCE7pDCAjmM/1TzOaEwTKr/VRkFgqIZub+qI2ZLWqchPiurD4JGv5FiagJjXFrn+RY0nuFzMKINo5vC40OTovVyaZXAtzMAp1xzAu9O13Ve/WZHEJ6vg3EWlo+C7p020rsJzCwepoG1uGx44pclWtm5b8IxsjK+jpOekQYdUNPTpII3lyR54jPoH4rqa67I11B5uBdoyFnWrSjD1duJFz2haQKzHrtNCxa7Y52hncCoPBIY1XHEinuUTH6DpwLfbmyTxDFowL3jDqknhAOYZFH9ONJFULG47l1lrOSSROD+cEnHmj6il9iH+ffcsnZKVsPVfl+nBcXoKUvQ73VZYlOQegIjXGi3rAD6jciissHJYXxfXhL4PhK+RcKifjvG1SMr5MAVJubDVwPkxTQVjjUEc26ZxkXonVJrAjNk0fJyCLoXoDALwz2KSyySXohLsi6e2jlpI0eT2pe4MydCdh5FuP9KowN+BGWwJzKI4VBarK1bh9t7GvhaD6NNCUNJSxfUYgWLHcA16YoBtFE6Fd59a/PUVIwOaNp3svUM9zHNwqyzBKbFF50qSeC4cIQS55LpfdIVW4fsJJZJCyMaV/3lofDdSlBjFibRvF21aLI4zhTirgnMOvy2oNQsxRKsuGX4ymzWSVdsyQAz+qDVsxo9s8Gjwwge62l8eqB3vmgwDR7t4yuNjkMJDmv6hHVRvmguWW6pytvYDKaNynv1IS2O4RuqtQLm/ggbl9hUmSE4NfsqGGIBLPGHDGy3yXAixIdI5jU4K9MEpnIJtAGcGyU0tuGXYrT+s+9Nlg6Y6QcttlEzYKOm/4ofKI9Y4tadz3/g1wvl7OlfB8GhLYyti4tjmsB0Vmbhq20QY+xCNDt+aMWHhku+fVo6ZqofdzrmPqVgevWiGw2MkPXuLQxI7cRpCAtwuFel7it96HsoC2PrCh3tph8mVCbHSOGwu5RF6LQhEJboZluEsK4TsQU++xwrgwbmeNsucHsUCJe7y2CeasCK8mO2EvzEO4dUIk1geitZX7C7E2NsPcYDY9BGtlvuWVZGZONmIWhPu1PwgtFVKCWCYiigweIIq9vbU7SbZXOKsfqUzDPDfrfNprEBuq6K1TEozOKExF22+f58twNjzlW2DahiPAkG3ASN4ONTejvWjXfYttYszLURGM+uj3piNSa1u/iPgLCm9qUBlvKTskuqx6Xia2uOy1l/viTgC2PbFw5u+I/agWWtoeWQibV1J5jVAiNEMAET/rLLAL3/qzoEIFz6qLxTGeEIjFKZd5iH0BX5tYJlPH42a4Kq4sHSqR85zOlglKw7MNbis3LJid0jrYUJCP0QHixN4Ines18sJEjHveM3QAFPAm7rahv/9PYKI8Gs3u8OVDZLEtNOZxozpo3RJHNNL8HjoZNZmGOgZTVeaA1RSljvnns1HEDtJ6rrultmDaKxzLD/RBlhLRHYCTTVzba296oExgtGbJNzEkKKtbHBepZEeAw+JmAO0RH/w4bZaQ6mhY6Y1NmkvVSPtxO98DBMSIwPWJuqAwzx033toGdMuU6QtHqzre1JlcCILPy6mKk1NywmAJzFGwE2M+Dxc4f11r7Xp+gwq7HXcY9tML8e7/XM3Q0Nf7jGNHDLflzKR98Pr+XXY+NJVQyjCcUEgjliMaI334k+wAdcgsENIDA8Uf5Z5VmCJHj6scjmKnCQwA16GEPlncqMiyZ0QZ2eWET6Qjt/42vyAEKw1zqJB8QQp8p8HhL4p/dw8KJ3joid4tQTZWPWBZ4jcBbflPIKN+yQBzlDOs0lO60SmCFWyp5hbLpgpmNwVid4YGE+Won+5jNTJqgnx6yYqsdJ/cHFwQIHCcGc6x0BQOs5rcW7YBCQu3q/0HOnJwcVCGhRoCk4aL/2SfMc8V91xDjckfSS6uF51WYXDhO4Hu5jKNS6ZLE5oLVJsfjlsqX/r8ETYLp3PdFIbEaYH7Rav1uyhLCw4YOweEi0Fp/MIygkfPOhFWQsF2QKbs64TuBA2tK15ICt78kdbQw2JseHPEsT/iUaJpZy8UsKHoIQlq+2UT3+kQaLIbA69UEYcAm+WJ09wadM8YVyinba95JEy6qfx++F6Gs0iPif+jlBdpYIDJt11sbMYszEL4l+CArWABq4Wyn97DuGDmEhMadUMguDy4eb1r0voC0X4KZwzq5fumCzB2odVuXAqi6ZNgEbFktV5M6wWQVLHwJD/FbcKQJI29QqzkrmhoE3lQyGo218c/74NK4aY/Mc+espRII1wcyNl+re6q8OB2x9L+6sTDPWglQav9jllgkYm5iEEOFCvVEuPloUrB1Z4xb2TmbUZrQhGCTKo5Mz19L+aRxIcGBVC6Mx7HTDBCAxbKhmY9smD5V6MSvwQ7ey8P2x4J5LQIZWCleLv5RotGEd3CGA6rhPct8KFY5xpcCY25UieAGxG8/x1JN0Xv1/XIpQNiILgtkyLc6mZEN+UXvSQqgNS0TMEqyB3t+ojvsP7lCif4haMAgDp1zd8TgOdoKmJ/Rw2gUNJMqMY8LiKlXmyHooWOAANkm36+z/ERxH0cRBi4LILq4t3kUX/OQ0MPB4i3EOjVPzYy+Ea4M16RFuvBKuO55VC8yahO0Tl5jQE1aVEXg0CQKEEigSGvVDsBHIcJ+j99lJeEwJ8IQOUlfAoe+hq71UJO7425ejD+F0F8B6RmFVjzBh5U3xcOjSO0FUmTb4AV0oINzm3reDKu8liRZoYPNCpynDMLbqmA9fjBTHoaFz5kX42Cvw4cGNFBgxgMljkewIvMcu1fMdE5si2t4FFkxgpt57m60LV/ouHNAGTu6SRheAquN0j5v3UVt3DLUjYJ/1nDwBFYzNITomeAWDYohuVNq3ShqXeXxQPldGYSA0SeUkeLyGx3qOBEr1i5KfO0rx/u1FGK5+J7QkMQvMjyViMTZmSbJ467QEuADG8EVBRTO0sYmmkrOSU0DCh4AyV6wZfymfjRFLCNOiTcgYyiih2YkxlVFezPusAAFzYe5rJre2ouHbnTWxXhVcfhHQvG/1zkXm/zztHDKwGfltT6m1QPOSUhvtsnXdf7NC5Yfi4MMOYXKjY0XfKyMMCA99ei6c8FD/RXlpgqfkzZNoJe79V3n2N4gZ4qAd/uxupMAwcTEUBvQ2BvULkmOk+u3zQ01cqGTS3IhdSoQKHHZxa/8JFJZgyBcTKuCvQmLuz5WrYxnx0hQhLmH7IxgwoTLZxjTGVqKb7q5FnDrxYoN3DwymkV5C4MogQEP3yYSqFM+h4Zj7ZPxZSKStq/M4bhd2amBpDhAPkYyxl6XD/ktwbK5ikhIJBq5WsFYqm0YeWhizoklcR9ZgJ5drkGXr6lzSJjCVLNUmYzNhZVIHCJUjLOp+ol7OhZjoHbNECA0HIgjTzj+DUE3gO5Zm5m4bvZYmO2V0FvvGxjC1XBz0h5mvtLnYaOaiDUD2WuRUp8QqxFwt3DKLY8y14VBgMmnunFA5QRsA2ylT7BCCLz2yR+QDXCVF1gClsUbCWu9Eo3PJmsCswdLLH7mBic1ibs06mLfFMhIqbQy+CkBQnuiJlo4JVZQqwccEYqd6DiFQJkVfT0SRz6ssVRglWPEcghJsLlkJy6ZhnLkWmJnv6R7bQuCSZDWsNi/CnToqtnsMYpmRUG1L+irYmXuJS5odTDxCWEjhkKUJzCVDqv4VY9FAbKyYO1KFe2Fn6ME65FIsfnHwmg8KABy4ZlctfmEOWJhgFahYmGw9z6x/ExjjRP3TguWsZi8Yxvrbs6DLCAR/+9Go1ldIIMDNrf5FCkb1zsp44cmA7bWJi2USApFLWHoXc+SACtpQKhfiQcB1p6BTAynjAAEvH/+hlWf76loUzP9bZbMMxBB8aEjiB3Zh0S6rsv+iEa1vABQOBOWDMpqT2ITvrqJ/1UVt3JivHYxruPlJdNhcTOPzOyYsCL+YjcWMwEXjqdLRhRdegSe4Y/S9kR9fMvEtkpjMBtzpefBYRjRwF5H72wtbsCCKU7TsLejXWCgcBKnqywv1R/GhBHtfcDeXLLrEiythML/JMCuxGNEKHXGpqrTsCjQYCly/nPtncGs8Oahw7mQlMnjHiWHPsjeBqeTqoPtHXz74RtVCW0x1cOEVLXynFnOdBuyrK2oM5srxddVYHg8u8kjwmsDUrVGvtxiNFmWxMOfHkIhBsHo3JTHXNZSVwxETvCYw628lXAICag4ADpq8AL8+Blq2ZoSfI3OtOk5Wf4J91m74PZ2bQgv6N1hJv3j81ubuBugbyg05oDXDSnFCGT00aBZmA+aL2Rwrn+t5cCuzwfSuLUqt1z1NDnc66dY1gdlu+WE6f8aJRWjpanCAIJ+fYtunTiOqm8CMWLJOhWc6lsYu3dZB3LBswgGtF5aFi8qXuQFaDJPjzgptWgi+xeIibfbt/wrDNxQFHNDa4AXYX5vp3bsMuzcLM+TI+mWOdjm94Vy/pePkAF4AfzQjKyyQ3gRm4wXUInA3w0d8LR0hB7Q+HCMjLEWXnf8HUfGWzbEuH2wAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle - T k \\log{\\left(\\frac{V^{N} \\left(\\frac{1}{\\Lambda^{3}}\\right)^{N}}{\\Gamma\\left(N + 1\\right)} \\right)}$"
      ],
      "text/plain": [
       "        ⎛       N⎞\n",
       "        ⎜ N ⎛1 ⎞ ⎟\n",
       "        ⎜V ⋅⎜──⎟ ⎟\n",
       "        ⎜   ⎜ 3⎟ ⎟\n",
       "        ⎜   ⎝Λ ⎠ ⎟\n",
       "-T⋅k⋅log⎜────────⎟\n",
       "        ⎝Γ(N + 1)⎠"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = F.subs( log(factorial(N)), N*log(N)-N ) # stirling's approximation\n",
    "logcombine(F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rotational Partition Function\n",
    "\n",
    "The following plots the rotational partition function and how each quantum number, $J$, contributions. We use `ipywidgets` to add an interactive slider to control the temperature.\n",
    "The manual for `ipywidgets` can be found [here](https://ipywidgets.readthedocs.io/en/latest/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "88f3bd63eca8448a97bb4edc79642f25",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=525, description='T', max=1000, min=50), Output()), _dom_classes=('widge…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.constants import h, c, k, N_A\n",
    "from scipy.misc import derivative \n",
    "\n",
    "B = 10.591 * 100 # 1/m # rotational wave-number (related to moment of inertia)\n",
    "J = np.arange(0,15,1)  # quantum numbers (called \"l\" in Dill's book)\n",
    "\n",
    "def q(J,beta):\n",
    "    ''' rotational partition function for a given J '''\n",
    "    return (2*J+1) * np.exp( -beta*h*c*B*J*(J+1) )\n",
    "\n",
    "def qrot(beta):\n",
    "    ''' partition function for \"inifinite\" number of quantum numbers '''\n",
    "    J = np.arange(0,500,1)\n",
    "    return np.sum( q(J,beta) )\n",
    "\n",
    "def qrotHT(beta):\n",
    "    ''' partition function in the high temperature limit, i.e. summation --> integral '''\n",
    "    return 1/(beta*h*c*B)\n",
    "\n",
    "def plot(T):\n",
    "    ''' plots partition function at given temperature '''\n",
    "    beta = 1/(k*T)\n",
    "    plt.plot(J, q(J,beta).cumsum(), 'b--o', alpha=0.7, lw=2, ms=10,\n",
    "             label=r'Explicit summation <$E$>={:.2f} kT'\n",
    "             .format(-derivative( qrot, beta, dx=0.001*beta ) / qrot(beta) * beta) )\n",
    "    plt.plot( [0,1000], [qrotHT(beta),qrotHT(beta)], 'r-', alpha=0.5, lw=4,\n",
    "             label=r'High temperature limit <$E$>={:.2f} kT'\n",
    "             .format(-derivative( qrotHT, beta, dx=0.001*beta ) / qrotHT(beta) * beta))\n",
    "    plt.xlim( 0, 14 )\n",
    "    plt.xlabel( 'Quantum number, $J$' )\n",
    "    plt.ylabel( 'Summed partition function, $q$' )\n",
    "    plt.legend(loc=0, frameon=False)\n",
    "    plt.title('Rotational Partition Function, $q_r$')\n",
    "\n",
    "# this part sets up a temperature slider to control the plot\n",
    "i = interact(plot, T=(50, 1000), continuous_update=False)\n",
    "i.widget.children[0].description=r'$T$ (Kelvin)';"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two-level system\n",
    "\n",
    "Here we consider a system where two energy level are possible; the ground state (0) and the first excited state, $\\epsilon$. The partition function is\n",
    "\n",
    "$$\n",
    "q = 1+e^{-\\epsilon / k_BT}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAGDCAYAAADK03I6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABlAUlEQVR4nO3deXhU5f3//+c7k0z2fSGsIeyyK3EBqbgiFgva1q1KbV272FZbl7b6rVattdaf2larRbtoUVv9VAWrBREBF0RZDCAICIQlbNn3dWbu3x9nEiaTSTKTZbLM+3Fdc03mnPucc080vM65z33uW4wxKKWUUio0hPV2BZRSSikVPBr8SimlVAjR4FdKKaVCiAa/UkopFUI0+JVSSqkQosGvlFJKhZAOg19Eviki/xGRAyJSKyK7ROS3IhLvx7ZRIvJ7ETnq3vZjETnLR7kwEfmFiOwXkToR2SIi3+jsl1JKKaWUb/5c8d8OOIFfAvOAp4HvAytFpKPt/wrcCPwKuBg4CqwQkele5R4A7gOeBC4C1gOvishX/foWSimllPKLdDSAj4ikG2MKvZZ9G3geOM8Y814b200DcoHrjDF/dy8LB7YDu4wxC9zLMoBDwMPGmHs9tl8FpBtjpnbyuymllFLKS4dX/N6h77bB/T60nU0XAI3Avz325QD+BVwoIpHuxRcCdmCJ1/ZLgCkikt1RHZVSSinln8527pvjfv+inTKTgDxjTI3X8u1YQT/Go1w9sMdHOYCJnayjUkoppbwEHPwiMhS4H3jXGLOxnaIpQKmP5SUe65vey0zrew7e5ZRSSinVReGBFBaROGAp4AC+21FxwFcHAulkOV/1uQm4CSA2NnbGhAkTOtpEKaWUGhA2bdpUZIxJD3Q7v4NfRKKAZcAoYI4xJr+DTUqAET6WJ3usb3pPFhHxuur3LteKMWYxsBggJyfHbNzYXgOEUkopNXCIyIHObOdXU7+IRAD/AU4DvmqM2ebHZtuBbBGJ8Vo+EWjgxD397UAkMNpHOYAd/tRRKaWUUh3zZwCfMOBF4DxgoTFmvZ/7XgZEAJd57CscuAJ4xxhT7168HOtE4Gqv7a8BPjfG5Pl5PKWUUkp1wJ+m/qewwvs3QLWInOGxLt8Yky8iWcBe4H5jzP0AxphcEfk38IS7xSAPa+CfbDxC3hhTICKPA78QkUpgM9bJwbnAwi5/Q6WUUko18yf4L3K/3+1+efo11oh7Atho3YLwXawThgeBJGALMM8Ys9mr3N1AFfATIBPYBVxujHnTny+hlFJKKf90OHJff6Gd+5RSSoUSEdlkjMkJdDudnU8ppZQKIRr8SimlVAjR4FdKKaVCiAa/UkopFUI0+JVSSqkQosGvlFJKhRANfqWUUiqEaPArpZRSIUSDXymllAohGvxKKaVUCNHgV0oppUKIBr9SSikVQjT4leqnHnvsMUaPHk1cXBxnnnkmu3bt6u0qKaX6AQ1+pfqhxYsX8/vf/54nn3yS3NxcYmNjefDBB3u7WkqpfkCDX6l+6IknnuDBBx/koosuYsyYMXzlK1+hsrKyt6ullOoHNPiV6meOHj3KF198wbx585qXrVy5kmnTpjV/fuKJJ5gwYQJRUVEMGjSIn/3sZ71RVRVicnJyOPXUUzu17euvv46IsGLFim6ulfKmwa96zZIlS7j55puZMWMGkZGRiAhvvPFGv6vP8uXLg1r33Nxc0tLSGDp0KGVlZfz85z8nNzeXG264AYB//OMfPPLIIzz22GN8+eWX/O9//+PrX/96UOqm/PPBBx9w4403Mm7cOBISEoiOjmbcuHF873vfY+vWrb1dvU7Jz89n06ZNLFy4sFPbX3jhhURHR7N06dJ2y3Xm783lchEXF4eI+PUaNWpUp75DfxHe2xVQoeuee+7hwIEDpKenM2jQIA4dOtQv67Ns2TKioqK44IILeriGlq1btzJt2jRefvllvvWtbwHw0EMPMXz4cAC+/PJLJkyYwJw5c4iNjW1ernpfbW0tN9xwAy+99BKRkZHMmzePBQsW4HA42LBhA4sXL+bZZ5/l448/5rTTTuvt6uJ0Olm5ciX79u0jKiqK2bNnM27cOJ9lly1bBtDp4I+JieH8889n2bJlPPXUU4hIm8fx/nvrqJ6VlZXcfvvtLfazZ88eXnzxRaZPn96qzmPGjOnUd+g3jDED4jVjxgyj+peVK1eaAwcOGGOMuffeew1gXn/99S7tc86cOSYrKyuo9Rk+fLiZP39+p45pjDF33XWXAdp8eX+fq666yvzsZz8zRUVF5qOPPjI///nPjYiYFStWGGOMOXz4sBk3bpwJCwszqampxul0drpuqvvU1dWZc845xwBmwYIFpqCgoFWZTZs2mWnTppljx471Qg1bO3r0aIv/F0XE3HXXXT7Lzp0712RnZ3fpeM8995wBzIYNG9os4+vvLZB6NvnjH/9oAPOb3/ymS3XuTcBG04m81Kb+EFRZWck999xDdnY2UVFRnHTSSSxZsoR//vOfiAgvvvhiUOpx/vnnM2LEiKAcyx+dqc/mzZs5dOgQCxYsaLWusbGRP//5z8ycOZP4+HhiYmKYM2cOH374YYtyv/zlLzl06FCbr/Xr17cov2XLFqZPn05qaiqzZs3it7/9LTNmzGDDhg00NDRw1VVXcckll7B582Y2bdpEWJj+mfcFv/jFL1i9ejWXXnopb7zxBunp6a3KnHLKKaxfv55Bgwb1Qg1bS0pK4rPPPqOuro4dO3Zw7rnn8rvf/Y7//Oc/LcpVVFSwZs2aNq/2161bx+WXX86IESOIjIwkNTWVU089lb/85S8tyn3ta18jLCyszeb+tv7e/K2np88++wyAk08+ucPfw4DTmbOFvvjSK37/FBcXmylTphgRMV/96lfNnXfeaS644AIDmPPPP98AZtu2bUGvV1+44u9Mfe69914jIubIkSMtlldVVZnZs2cbwJxxxhnmtttuM9/97ndNbGysiYiIMOvXr+9Uverq6kx4eLj59NNPm5fV1NSYpKQks3r1avPiiy+a6dOnd2rfqufs3LnThIWFmYyMDFNSUtLb1em04uJiExUVZc4777wWy//1r38ZwKxevbrVNi+++KIJCwszQ4cONTfddJP55S9/ab7//e+bs846yzz22GOtys+aNctMmTLF5/Hb+nvzt56epk+fbgBz9OjRdvfVl9HJK369xx9ivv3tb7Nz506WLVvGxRdf3Lz85ptvZvHixURGRjJhwgSf2953330BHSvQ8v3RsmXLyMnJYfDgwS2WX3PNNaxbt46///3vfOc732le/pOf/IRTTz2V+++/n7feeivg423fvh2Hw8G9997Lo48+isPh4Fe/+hUTJ07krLPO4vnnn2fPnj289tprzJgxg5KSEj7//HMWLVrU1a+quuDxxx/H5XJx6623kpyc3GPH6em/0ZSUFKZPn86WLVtaLF+6dCkpKSnMnj271TZ33303o0aNYvPmzcTHx3d4jIULF3LXXXeRl5dHdnZ2i3Vt/b35W88mDQ0N7Nixg8zMTDIzMzus00CjwQ9c/48NvV0Fv/z1O517TKbJihUreOutt/jhD3/YIvQBFi1axOLFi5k8eTLh4b7/t/j1r38d0PEGevDn5+fz2Wef8cADD7RY/t///pc33niDO++8s0XoA0ybNo2TTz6ZjRs3duqYW7ZsITs7m4qKCnJyckhJSWH+/Pk899xzhIWFcfXVV7Nx40ZuueUWiouLyczM5LrrruvsVwyKW1bd0ttV8MuT5z3Z6W2bTvKuuuqq7qqOT939N/p///d/3H333eTn53Pqqafy2GOPER0dTUVFRXOZxsZG3n77bRYsWODz347i4mLi4uLYsWMHJ598Mna7vd1jNgX/0qVLufXWW5uXt/X35m89PW3fvp2GhobQbOZHgz+k/PnPf8Zms3HnnXe2WpeSkgLQ4llwb1bLUt/RVq/fttZ1d/2bejF732/885//DEBNTY3Pf1gLCws7XZetW7cye/ZsXnjhBZ/r7XY7Tz31FE899VSn9q+6X2FhIfn5+QwaNIiRI0f26LG68//xDRs2cOWVV7Jw4UKWLFnCmjVruOCCC4iJiWnRB2HLli2Ul5czZ84cn/v5wx/+wG233cYZZ5zRvCwvL6/N38X48eMZNGgQa9eubRH8bf29+VtPTyF9fx8N/pBhjGHVqlVMmTLFZwe2I0eOADB9+vQg16zz7r333lbL/vGPf1BWVtbiH4yesmzZMrKyspg6dWrzMmMMa9asAeDJJ9u+QuzsPzhbtmxpMXCP6vuKiooA+kyHPX/97ne/IyYmhn/84x/Ex8dz6qmnUlRUxCOPPMKsWbOay02ePJnY2Fg++eQTrr/++hb7cDgclJSUMHLkSCZOnMioUaOw2+1kZWW1edwDBw5w/PhxTj/99BbLff29BVJPTxr8KiQUFxdTXV3d5jPdq1atAtoP/r52j9/X/tesWcP+/ft7/NhVVVWsWbOGm266qcXy4uJiamtrufrqq1myZEm3H3fr1q2tnkdWfVtMTAwAx44dC2i7hIQE/v73v/ONb3yjedlTTz3Fn/70J7Zv347NZmu1TXf+ja5cuZLzzjuvxX35efPm8cgjjzBz5szmZVFRUcydO5dly5bxzDPPtHiK5LrrrmPlypV8/PHHfrd2NA3M4/mEQFt/b4HU05MGv+ryvfP+oOm+WklJSat1RUVFLF68GBFpt6lf7/GfsHz5curr61s1O7pcLuDEVV5366n99qau3DvvD0aMGEF6ejoFBQXs3Lmzzc6zYP3/0xScEydO5IsvvmheV1dXx0MPPcTjjz/uM/Sh+/5Gy8rKqKioYPz48S2Wv/rqqwBccsklLZZfcsklvP7663zyySfNYXvw4EH++c9/cscddwR0i2Pp0qWMHTuWk046qXlZW39vgdYTrFa5LVu2kJCQMOBH6GuLPuAbIhISEhgyZAiffvope/bsaV5eUVHBZZdd1twcl5CQ0OY+An1kZCBbtmwZiYmJre5rZmRkMHLkSFavXs2nn37aaruqqio2b94crGqqPkBE+MEPfgDA9ddfT0FBQasyJSUl3Hvvvc2BBTBp0qQWwf/UU0+RkZHBZZdd1uaxuutvNCYmhoiICKqrq5uX7d+/n+eff55TTz21Ve/9+fPnY7PZWjx/3/Q9c3NzWx3L6XRSWlra6rilpaV88MEHrQK7rb+3QOsJ1siWVVVVTJ8+vd1+QgOZX1f8IjIMuAvIAaYB0UC2MWZ/B9vdB7S+EWupN8ZEeZTdD/i68XOpMeYNf+qp2vfTn/6U22+/nVmzZnHFFVdgjOG1117jlFNOAYJ/f/+5555rHswmNzcXgD/+8Y/NTX033HCDzz/c3q6P0+nk7bffZt68eURERLTaz+9+9zuuvPJKzjzzTBYsWMDYsWMpLy/n0KFDrFmzhp///OfNv3MVGu6++242btzIW2+9xbhx45g7dy6jRo2irq6O7du38+GHH1JXV8e2bduat5k0aVLzYFpVVVU8/PDDPP/880EJK7vdzoIFC3jllVe4+eabiYmJ4etf/zo2m43nnnuuVfnU1FRmz57N0qVLefjhh5vrn5WVxcqVK5k8eTLnnXceUVFR5Ofn89577/HCCy8wd+7cFvt56623cDgcLZr52/t7C7SeoM38gH9niMDZwHHgbWAF1rCII/3YbhhwhtfrPKAReMWr7H5guY/yyf7UUQfw6ZjT6TT33HOPyczMNFFRUWbatGnmmWeeMatWrTKA+fWvfx3U+lx77bXtDlX797//PeB9dmUAH3/rs3btWgOYF198sc19vfvuu+bCCy80ycnJJiIiwmRmZppZs2aZ++67z+Tn53eqfqp/czqd5vnnnzfnnXeeSUlJMeHh4Wbw4MHm1FNPNT/72c/M+++/b1wuV3P55cuXm5iYGONyucwDDzxgzjzzzKDWt7Cw0MyfP7/5//+pU6e2O5TuY489ZgCza9eu5mV79uwxV199tRkyZIgJDw83CQkJZsqUKeb73/++KSsra7WPb37zmyY9Pb3FMNMd/b0FWs+mIbI78+9LX0MnB/DxN/jDPH6+wd/gb2Nfi9zbz/davh9Y0pl9Gg3+Lnn88ccNYJYuXdrbVekXfvazn5nw8PB+PQKb6vsOHTpkAJObm2uSkpLM2rVre6UeZWVlfp2s7t271wDmkUce6dRx6urqTHx8vLnuuutaLPf3783feg4knQ1+v+7xG2Nc/pTz07VYrQc66XIf0TS6VXsd+9QJb775Jl/5yld6dAQ2pYYNG0ZiYiI33XQTp59+OmeddVav1CMxMZGhQ4d2WG7UqFFMnjy5w2l12/Lee+9RWVnZarx/f//e/K2nCnKvfndfgXOAJ4wxDh9FviYiNYAN+Ax42Oj9/R63ZcsWkpOT2322Vp2wa9eu3q6CChETJ07k448/7vRIj8H24YcfUl9f36lt58yZw9GjR1tNXqR/b90v2I/zLcJ6kuB5H+veBDYAecAg4BbgdRFZZIzx+UC0iNwE3AT0qVne+hOHw8GOHTvafN5VKdV71q1b19tVCEhiYmKnt42JiWke80D1LDEBPnYlIjcAz+JHr34f234B1BpjOuzSLCI2YD2QaYzxPeqMh5ycHNNfzoqVUkqprhKRTcaYnEC3C9pz/CJyGjAB31f7rRhjnMCrwDARaX8qJqWUUkr5JZgD+FwLOICXAtim6YHVgT0ajFJKKRUkQQl+EbEDVwJvG2MK/dwmHLgMOGiMCWyQa6WUUkr55HfnPhH5pvvHGe73i0SkECg0xqwVkSxgL3C/MeZ+r80vBlJoo5lfRK4CFmINEHQIq3PfD93H6tkJrJVSSqkQEkiv/le9Pv/Z/b4Wa2Q/wXoMz1crwrVACfDfNvadB2QAv8c6QajB6uE/zxijz/srpZRS3cTv4DfGtDtAtLuHv88yxpiFvpZ7rF8PnOtvXZRSSinVOTo7n1JKKRVCNPiVUkqpEKLBr5RSSoUQDX6llFIqhGjwK6WUUkES6DD5PSHYk/QopZRSA5Ixhsp6B8VVDZRUN1Babb2X1Fg/F1c3MGlIAt89M7tX66nBr5RSSvmhpsEK9dIad6B7vJqWOZztX9EXVXVu2uLupMGvlFIq5DmcLkprGimqqqfEfXVeUlVPSU0jJdX1lFY3Utfo7PJxSqobuqG2XaPBr1QnTZw4kby8PCZPnsyzzz7L9OnTe7tKSqk2OJwuSqobKKpqoLi6nuKqBoqq6ilyv5fVNBCM2++l1Y0YYxBpd0y8HqXBr5Tba6+9xjPPPMOmTZsoKSkhLy+PkSNHtls+MzOThx9+mEsuuYQ9e/YQHm79SW3atInvfe97bNu2jenTp3PzzTdz2223UVZWFpwvo1SIaWwOdivMi6tahnt5bXCCPcpuIyXGTkqs9UqOtZMa6/E5xt6roQ8a/Eo1q66uZvbs2XzjG9/ge9/7XoflJ0yYAMBvfvMbnnzySTZs2MDMmTPJzc3lrLPO4r777uOVV17hjTfe4Mc//rG2CCjVBcYYymsbKaisp7CynoLKOgorT1yxl9c09ngdImxhJMfaSYmNICU2kpTYCJJj7KTGRpIcG0FKrJ0Ye9+P1b5fQ6WCZNGiRQDs2bOnw7Lr16/nlltuoaGhgYceeojExESKiooA+NGPfsQtt9zCHXfcAcBtt93GE088wdSpU3uu8koNAI1OF8VVDRRU1lFQUU9RVT0F7pAvqmyg0enq0eMnxkSQFhfZ4grd88o9PjK816/Wu4MGv1IBcjgcXH755dxzzz1ceumlXHPNNRQWFpKRkcHevXv58MMPeemll1psY7fbNfiVAqrrHc1X7Z5X7gWVPXufXQQSo+2kxdlJjbNbAR8XaX2OjSQl1o49PDSGttHgD1E5OTmICBs2bAhou9dff52vf/3rLF++nAsvvLCHate3rVu3joqKCm688UZEhJ/+9KesWrWKCRMm8O6775KUlMTw4cOby1dXV3PgwAENfhUyahucHK+o43hFHcfc78crrHCvqXf0yDFFICnGCvX0uMjmcE+JtT4nx9qJsIVGsHdEgz8E5efns2nTJh544IGAt73wwguJjo5m6dKlIRv8x44dY+zYsc1Nfu+88w5z5swhMTEREcHpdOJyuQgLs/6R+dOf/oTD4WDy5Mm9WW2lupXD6aKwqp5j5e6AL6/jeGU9x8vrKK/tmfvt0XYb6fGRZMRHud8jSYuLJC3eTkqMnXANdr9o8IegZcuWAbBw4cKAt42JieH8889n2bJlPPXUUwPiflegJk2axOHDh6murmbnzp08/fTTvPPOO4DVklJbW8uDDz7IokWLWLFiBY8++igjR44kPj6+l2uuVGCMMZRUN3Cswrrnfswd8E3N8z3RLJ8ca28O9fT4SNLjIslIsII+LlIjqzvobzEELV26lOzsbKZMmdKp7RcuXMibb77Jpk2byMnJ6eba9X2TJk3i1ltvZfr06URGRvLyyy8ze/ZsAEaMGMGTTz7Jfffdx+OPP86VV17JFVdcweHDh3u51kq1rcHh4nhFHUfKajlSXsvR8jqOl1vN893doS7cJu5Aj2oZ8O6r91C5z96bNPhDTEVFBWvWrOEHP/iBz/Xr1q3jiSeeYP369Rw/fpy4uDhGjRrFDTfcwM033wzA1772NcLCwli6dOmACv6SkhIOHjzIoUOHANixYwdlZWWMGTOGuLi4FmXvvPNO7rzzTp/7ufnmm5t/VwDz5s3jtNNO67mKK+WnukYnx8rrOFJey5GyOo66g767r97Dwqxwz0yIIjMhioyESDITo8iIjyI5JiIkWwr7Eg3+EPO///2PhoYGn838L730EosWLWLw4MHMnz+ftLQ0SktL2b59OzU1Nc3lMjIyOOOMM1i6dGmn+gn0VcuWLeO73/1u8+f58+cDsHr1as4+++xO73fbtm3ccMMNXa2eUn6ra3RaV+9lJ67ij5TVUlzVvcPFJsfaGZRgBfwg9yszMYrUWL3f3pdJX5gisDvk5OSYjRs3dm7jl67o3sr0lG/9u+u7+Na3WLFiBcePH28eZa5JdnY24eHhbN68ucP70Y888gh33XUX+/btIzu7/Zmm7rvvvoDqGGj5vqy0tJSUlBR27drFuHHjers6aoCpaXB4Bbz1XtqN48FH220MTmwZ7IPirav4qAhbtx1HBU5ENhljAm521eCHkAn+xsZG0tPTWbBgAS+88EKr9QkJCcTFxfH6669z8sknY7fb29zXrl27mDBhAo8//ji33npru8cNtFlvoPw/qVR3cbkMBZX1HCqt4VBJDYdKaskvrenWCV9S4+wMSYpmSGI0g5OimsM+boAMWjMQdTb4tak/hGzZsoXy8nLmzJnjc/0f/vAHbrvtNs4444zmZW2NVz9+/HgGDRrE2rVrOwz+ngzyQP9BamxsbNXSoVRfUl3vIL+01gr40hryS2s5XFrbLZ3sRCA9PtId7tEMSYpiSGI0mYlRevUeQvRfwBAyefJkYmNj+eSTT7j++utbrHM4HJSUlDBy5EgmTpzIqFGjsNvtZGVl+dzXgQMHOH78OKeffnowqt4mz050/mh6tl6p3uZyGY5X1p0I+ZJaDpXWdEszvYgwKCHSuoJ3h/uQpGgGJURpr3mlwQ90y73z/iAqKoq5c+eybNkynnnmmRYheN1117Fy5Uo+/vjjdmeka/LGG28A/o0F0JP3+J955pmA9t0Wf1sOums6Tb2dEVpqG5wcLKnp9qv4sDAhMyHqRMAnuQM+PlI716k26T3+EPPCCy9w7bXXsm7dOmbOnAnAwYMHycrK4o477uCRRx7xaz/nnnsu+fn57N69u8Oyeo9fhZK6RicHims4UFzN/uJq9hfXcLy8rsv7jYsKZ3hyDMNTohmWHMPw5BgGJ0XpMLQhTO/xK7/Mnz8fm83G0qVLm4O/oKAAgNzc3FZXtE6nk4qKCpKTk5uXlZaW8sEHH3Dbbbf5dcz+FOSbNm3ie9/7Htu2bWP69OncfPPN3HbbbZSVlQW8r8cee4ynnnqK48ePM23aNP72t78xfvz47q+06jV1jU7yS2vIKzoR9MfK67r0THxYmDA4MYrhyTEMS45meIr1nhitz7+r7qHBH2JSU1OZPXs2S5cu5eGHHwaskeiysrJYuXIlkydP5rzzziMqKor8/Hzee+89XnjhBebOndu8j7feeguHw9GpIX/7stzcXM466yzuu+8+XnnlFd544w1+/OMfM3369ID3tXjxYn7/+9/zt7/9jbFjx/KDH/yABx98kH/+85/dX3EVFPUOJ4dKatlfVM2Bkhr2F1VztLy2SyEfHxXeHOxW0OtVvOp5GvwhaOHChfz0pz9l9+7djBs3jujoaFatWsW9997L6tWrefrpp4mJiSErK4uvf/3rrTrwLV26lPT09OYWg4HiRz/6Ebfccgt33HEHALfddhtPPPFEm7PqTZw4kby8PCZPnsyzzz7b4gThiSee4MEHH+Siiy4C4Ctf+QqbNm3q8e+gukeDw8Wh0hoOFteQV1TNgeJqDpfVdbr1SgQyE6PISolt0VSfEK2Pyqng8yv4RWQYcBeQA0wDooFsY8x+P7Zt6y/lZGNMrke5MPcxbgYygV3A/caY//hTR+W/puBfunRpc8iNHj2aJUuWdLhtfX09K1as4LLLLhtQPeT37t3Lhx9+yEsvvdRiud1ubzP4X3vtNTIzM3n44Ye55JJL2LNnD+Hh4Rw9epQvvviCefPmNZdduXIl55xzTo9+B9U5xhiKqxvYW1DFnsIq9hRUkV9ai8vV+ZAflBDFyNRYslJjyE6LZXhKjD4up/oMf6/4xwCXA5uAD4C57Rdv5R/AX7yWefcKewC4HbjbfZwrgVdF5GJjzNsBHk+1Y9SoUUyePLlF8Pvrvffeo7KyckA28yclJTF8+PDmZdXV1Rw4cKDN4J8wYQIAv/nNb3jyySfZsGEDM2fOJDc3l7S0NIYOHUpZWRkPP/wwubm5vPjii0H5Lqp9jU4XB4pr2FNQxd5C61Ve0/lpZDMSohiZGsPItFhGpsYyIiWGaLuGvOq7/A3+940xgwBE5AYCD/7Dxpj1ba0UkQys0H/YGPOoe/FqERkDPAxo8HezDz/8kPr6+oC3mzNnDkePHiU9Pb0HatV7RASn04nL5WpuyfjTn/6Ew+Fg8uTJrcqvX7+eW265hYaGBh566CESExMpKioCYOvWrUybNo2XX36Zb33rWwA89NBDLU4qVPCU1TRYAV9QzZ7CKvYXVePs5NV8RkIkWamxzVfzWakxxNj1jqnqX/z6P9YY073zMrZ2IWAHvNualwB/E5FsY0xeD9chpCQmJnZqu5iYGGJiYrq5Nr0vJyeH2tpaHnzwQRYtWsSKFSt49NFHGTlyZKt5CxwOB5dffjn33HMPl156Kddccw2FhYVkZGQA1giJ06dPZ+7cuXz00Ue8+eab3H333cyYMaNFJ0nV/ZwuQ35pDXvdTfZ7Cqo6PTFNWlwkI9OsgG8K+lidD14NAMH6v/j7InIH4ATWA/caYz7wWD8JqAf2eG233f0+EdDgVz1mxIgRPPnkk9x33308/vjjXHnllVxxxRUcPny4Vdl169ZRUVHBjTfeiIjw05/+lFWrVjU3/W/ZsoWvfvWrpKamMmvWLGbNmsW7777Lhg0bNPi7WVW9g70eTfb7CqtpcAR+nRJhCyM7PZYx6XGMzohjVHosCVERPVBjpXpfMIJ/CfBf4AiQBdwBvCciFxhj1rjLpABlpnWX2RKP9a2IyE3ATWD9w61UV9x8880thgCeN28ep512Wqtyx44dY+zYsc29sd955x3mzJlDYmIi9fX17N69u8Xz+rW1tezZs4czzzyz57/EAFdd72DX8Up2Hatk59EK8ktrO7WflFg7YzLiGJ0ex5iMOIYlR+tIdypk9HjwG2MWeXz8QESWAp8DDwKz3csF8HXTrd3nXIwxi4HFYI3c1/XaKnXCtm3buOGGG1otnzRpEocPH6a6upqdO3fy9NNP88477wCwfft2HA4H9957L48++igOh4Nf/epXTJw4kbPOOivYX6Hfq2lwsPt4FTuPVrDzWCX5pTUBPzdvCxOyUmOaQ350ehzJsW3PPKnUQBf0G1bGmEoReQvwnCWmBEgWEfG66k/2WK9U0JSWlnLkyBGfPfonTZrErbfeyvTp04mMjOTll19m9mzrHHbLli1kZ2dTUVFBTk4OKSkpzJ8/n+eee25APf7YU2obnHxZUMnOo5XsPFbJwZLqgIM+ITrCHfCxjMmIY0RKrE5Mo5SH3uqp4n2Fvx2IBEbT8j7/RPf7jiDVSykAkpOT2x2s5c477+TOO+9stXzr1q3Mnj2bF154oSerN2DUNTr58ngVu45bTff7i2sCGiRHBIYlxzA6I859fz6W9LhIHRRHqXYEPfhFJAGYD3zisXg50ABcDfzaY/k1wOfao1/1F1u2bGkxcI9qqd7hZE9BFTuPVrLreCV5RdUBD5QzLDmaCYMTGJ8Zz7hB8cRpT3ulAuL3X4yIfNP94wz3+0UiUggUGmPWikgWsBdrtL373dvcDowHVnOic9/tWCPzXd20b2NMgYg8DvxCRCqBzcAVwLnAwBopRg1oW7du5fbbb+/tavQZxhj2F9ew7XA524+Uk1cY+DP0Q5KimTA4ngnuoI/X3vZKdUkgp8qven3+s/t9LXA2VvO9DfC8mbYLuNT9SgQqgI+A640xn3rt726gCvgJJ4bsvdwY82YAdVSqVzUN4hPKKusa+fxwBZ8fLufzI+VU1TkC2j4zMYoJgxOYkBnP+Mx4faxOqW4m/WnK1Pbk5OSYjRs39nY1lAo5Lpchr7iazw+Xsy2/nP3FgXXIG5QYZYX8oHgmZCaQGKNBr5Q/RGSTMSYn0O305phSKmAVdY3WFf3hcj4/XEF1vf9X9RkJkYwfFM/4TOuqXh+tUyq4NPiVUh1quqrfll/OtsPlHAjgqj7abmPikASmDE1k4uAEUuMie7aySql2afArpXyqqGvkc3fQbz8S2FX98JQYpgxNZMqwREalxeqoeEr1IRr8SinA6oGfX1rLpgOlbM23rur9FW23MWlIIlOGJjJ5aAJJMdp8r1RfpcGvVAgzxrCvqJpNB0r57GApBRX+T9U8PCWGqcOssB+VHoctTAfNUao/0OBXKsQ4XYYvCyrZdKCUzQfKKKvxb9ramMhwJrnv1U8ekqi975XqpzT4lQoBjU4XO49WsvFACbmHyvx+tj4rNZYpw6ywz07Tq3qlBgINfqUGqHqHk88Pl7P5QBm5+WXUNTg73MYWJkwaksiMrGSmDEskMVqv6pUaaDT4lRpAahocbDlUzuaDpWzLL6fR6epwm8iIMKYMTeKUEUlMHZZEtN0WhJoqpXqLBr9S/VxFXSO5B8vYdKCUL45W+DUWfrTdxvThSczISmbSkESdtlapEKLBr1Q/VNPgYOP+UtbvK2b38Uq/BtOJjwrnlKxkThmRzITMeH22XqkQpcGvVD/hcLrYfqSCdXuLyT1UisPZcdonx9qZ4Q77sRlxhGnnPKVCnga/Un1Y07S2H+8t5tO8Yir96I2fkRDFjKxkZmQlMzI1BhENe6XUCRr8SvVBxVX1rN9Xwrq9RRwrr+uw/LDkaE5xh/3QpGgNe6VUmzT4leojahucbDpQysf7ith1rOP79mlxkcwcncoZo1LJTIwKTiWVUv2eBr9SvcjpMuw4UsG6vUV8drCsw8fvou02TstOYeaoVMZkxOmVvVIqYBr8SgWZMYZDJbWs21vEJ3klVNQ2tls+LEyYOjSRmaNTmTosSR+9U0p1iQa/UkFSWt3A+n3FfLyvmMOltR2Wz06LZeboVE7NTiEhSkfQU0p1Dw1+pXqQy2XYkl/Gml2FbD9S3uF9+9Q4O2eMSmXW6DS9b6+U6hEa/Er1gPLaRj74spA1uwoprW5/9rsou42crGRmjU5j3CC9b6+U6lka/Ep1E2MMu49XsXpXAZsPlLY7dK6IMMV93376cL1vr5QKHg1+pbqotsHJx/uKWL2zkCNl7d+7H5Eaw6zRaZw+Su/bK6V6hwa/Up10qKSG1bsKWL+vmPrGth/Di7CFcfqoFM4Zn8HItNgg1lAppVrT4FcqAI1OFxv3l7JmVwF7CqraLTsoMYpzxmcwa3QqsZH6p6aU6hv0XyOl/FBYWc/a3YV88GUhVe2Mly8inDwiiXMnZDAhM1476iml+hwNfqXa4HIZth0uZ/WuAj4/3P6jeIkxEcwZl85ZY9NJjrUHr5JKKRUgDX6lvFTUNfLhl0Ws2VVAcVX7j+KdNDiBcyakM21Yks5vr5TqFzT4lXI7Vl7Hiu3H+GhPUbuP4kXbbcwek8bZ4zN0kB2lVL/jV/CLyDDgLiAHmAZEA9nGmP0dbJcD3AScBYwAioAPgHuMMXleZfcDWT52c6kx5g1/6qlUZ+wtrGL558f47GBpu835WamxnDMhndOyU4gMtwWvgkop1Y38veIfA1wObMIK7rl+bnclMAn4I7AdGAr8P2CjiEw3xhzyKr8CuM9r2S4/j6WU34wxbM0vZ/n2Y+w+VtlmuQhbGKdmp3DuhAyy9VE8pdQA4G/wv2+MGQQgIjfgf/D/zhhT6LlARD4C8oAbgV95lS8yxqz3c99KBczhdPFJXgnLPz/W7mA7qXF2zjtpEGeOSSNOH8VTSg0gfv2LZoxpf5Lwtrcr9LHsgIgUYl39KxUUdY1O1u4uZOWO4+2OnT88JYavThnMjKxkbGH6KJ5SauAJ+qWMiJwEZABf+Fj9NRGpAWzAZ8DDen9fdUV5TSPvfnGc1bsKqG1wtllu4pAE5k3OZOLgBH32Xik1oAU1+EUkHHgGKAT+6rX6TWAD1m2AQcAtwOsissgYs6SN/d2E1XmQESNG9FS1VT/kTw99ETh1ZAoXTR7MiNSYINdQKaV6R7Cv+J8EZgHzjTGlniuMMT/y/CwirwPrgd8CPoPfGLMYWAyQk5PTwUznKhT400M/whbGV8alMXdiJunxkcGtoFJK9bKgBb+I/Bbr6vxaY8w7HZU3xjhF5FXgdyIy2BhztMcrqfolf3vox0aGc95JGZw7IYN4nRlPKRWighL8InI38HPgx8aYfwayqftdr+ZVKw6ni0/zSvifHz30503O5Mwxafr8vVIq5PV48IvIj4EHgbuNMX8KYLtw4DLgoDHmWE/VT/U/DqeLj/cV8+aWI+0OqTs8JYaLJmeSMzJFe+grpZSb38EvIt90/zjD/X6R+7G8QmPMWhHJAvYC9xtj7ndvcyXwBLAceE9EzvDYZYUxZoe73FXAQuBt4BBW574fuo91VSe/mxpgmgL/v1uOUlRV32Y57aGvlFJtC+SK/1Wvz392v68FzsZqlrcBnjOVzHMvn+d+eWraDqye/BnA74EUoAarh/88Y8yKAOqoBiCny/Dx3mL+u/UIhZW+A7+ph/68yZlkpeoIe0op1Ra/g98Y0+6lk3vcfvFa9h3gO37sez1wrr91UaHB6TKs32cFfkGF78C3hQlfGZfOvEnaQ18ppfyhY5GqPsfpMnyyr5g3Owr8sWnMnzqElFh7kGuolFL9lwa/6jNcLsP6vGLe3HKUgoo6n2VsYcLssWl8dcpg0uL0Cl8ppQKlwa96nctl+CSvhDe3HuF4ue/ADwsTZo9JY/5UDXyllOoKDX7Va1wuw4b9JSzbcoRj7QT+maNTuXjaEA18pZTqBhr8KuiaAv/NrUc4WuY78EWEM8ekcvHUIdppTymlupEGvwoaYwwbD5SyNPdwu4E/a3QqF08bTEZ8VJBrqJRSA58Gv+pxxhg2HShlae6RNofWFRFmjk7la1MHk5Ggga+UUj1Fg1/1GGMMnx0qY+lnh8kvbSvw4YxRqSyYNkQDXymlgkCDX3W7ptny3sg9zMHiGp9lmgL/4qlDyEzUwFdKqWDR4FfdxhjD9iMVLM09zL7Cap9lROD07FS+Nk0DXymleoMGv+oWXxyt4I3cw+w5XuVzfdNY+gumD2FwYnSQa6eUUqqJBr/qki+PV/JG7mF2Hq1ss8wpWcksnD6EYckxQayZUkopXzT4VafsLaxi6WeH2X6kos0y04cnsXD6UEakauArpVRfocGvArK/qJo3cg+zLb+8zTKThyZyyclDyU7T6XGVUqqv0eBXfjlUUsMbnx0m91BZm2UmDklg4fShjMmIC17FlFJKBUSDX7XrcFktS3MPs2l/aZtlxmXGc8n0oYzPjA9izZRSSnWGBr/y6Vh5HUtzD7NhfwnG+C4zJiOOS04eyoTMeEQkuBVUSinVKRr8qoVDJTWs2H6M9fuK2wz8kWmxXHryUCYNSdDAV0qpfkaDX+F0GT47WMqqnQXsPtb2Y3nDU2K49OShTB2WqIGvlFL9lAZ/CKuoa+T93YWs2VVIaXVDm+WGJkezcPpQThmRpIGvlFL9nAZ/CNpfVM2qnQV8mleMw9lGez6QmRjFJScPJScrWQNfKaUGCA3+EOFwuth0wGrO31vge1jdJuMy4zn/pAxOHp5MWJgGvlJKDSQa/ANceU0ja3YXsHZXIeW1jW2Wi7CFMWtMKueMz2B4io60p5RSA5UG/wBkjGFfUTWrvjjOxv2lOF1tN+enxUVyzoQMvjI2jdhI/d9BKaUGOv2XfgBpdLrYkFfCqp0F7C/yPS1uk4lDEjjvpEFMHZqozflKKRVCNPgHgJLqBtbsKuD93YVU1jnaLBcZEcas0Wmcd1KGTo2rlFIhSoO/H2p0ujhSVsvBkho+P1zBpgOlmLZG2wEyEqI4b0IGZ45JI9puC2JNlVJK9TV+Bb+IDAPuAnKAaUA0kG2M2e/HtlHAA8A1QBKQC9xljHnfq1yY+xg3A5nALuB+Y8x//PsqA1Ndo5ODJTUcLK7hQEkNh0pqOFJW2+59ewARa5a88yYMYvJQHWFPKaWUxd8r/jHA5cAm4ANgbgDH+CswH7gD2Af8EFghIjONMbke5R4Abgfudh/nSuBVEbnYGPN2AMfrtyrqGjlYXMPBkhoOuN8LKuoC2ke03cbsMWmcOyGDjISoHqqpUkqp/srf4H/fGDMIQERuwM/gF5FpwLeA64wxf3cvWwtsB+4HFriXZWCF/sPGmEfdm68WkTHAw8CACn5jDEVVDRwsqXZfzddyoKSa8pq2H7fryOCkKM47aRAzR6USFaHN+UoppXzzK/iNMa5O7n8B0Aj822NfDhH5F/BzEYk0xtQDFwJ2YInX9kuAv4lItjEmr5N18Jsxhh1HKzAGXMY0v7uMtc4ALpfXZ8/1XuVbvGOoqnNYQV9SQ22Ds0t1TY61k5USw4jUGCZkJjBuUJw25yullOpQT3fumwTkGWNqvJZvxwr6Me6fJwH1wB4f5QAmAj0e/E6X4bF3dvf0YQIiAoMSohiREkNWagzDU2LISo0lTp+5V0op1Qk9nR4pQKmP5SUe65vey0zrrune5XpUWC9fMdvChGHJMYxIiSYrNZbhKTEMS47WpnullFLdpqeDXwBf3c+9E9bfci1XitwE3AQwYsSIztTPa39d3oXfouw2RqTEtHgNTowi3BYWvEoopZQKOT0d/CWAr0RO9ljf9J4sIuJ11e9drgVjzGJgMUBOTk77z7f5QUSYOCQBcf8cJoIIhEnrz9bPguD+HOb5WZrLICc+h9vCGJJoNdunx0fqPXmllFJB19PBvx24VERivO7zTwQaOHFPfzsQCYym5X3+ie73HT1cz2Y/mzs+WIdSSik10DXWQvlhKD8E5fmQMgpGntmrVerp4F8G/Bq4DHgeQETCgSuAd9w9+gGWY50IXO0u3+Qa4PNg9OhXSimlOq2xDircAV/mDvnyfKgpalkua1b/CX4R+ab7xxnu94tEpBAoNMasFZEsYC/WaHv3AxhjckXk38ATIhKB1TP/+0A2VsjjLlcgIo8DvxCRSmAz1snBucDCLn1DpZRSqrs0B3y+x+sgVBd1vC1Y5XtZIFf8r3p9/rP7fS1wNlZHPBvg3Tvtu8BvgAexhuzdAswzxmz2Knc3UAX8hBND9l5ujHkzgDoqpZRSXedogAqPcC87ZF3NVxfhuy+6nyqOgMsJYb33tJa0N7lLf5KTk2M2btzY29VQSinVnxgDVQVQdgBKD0DZfivkqwroUsA3E4jPhMRhkDjceh9+Oti6fqddRDYZY3IC3U5HgVFKKRUaGuug7KD7td8d9AfBEdicKL4JxA9qGfCJwyB+CITbu2H/3UeDXyml1MBijNUkX3YASvefuJrvlqt4gbiMluGeOAwShva5gG+LBr9SSqn+y1FvNc03X8G7r+Iba7u4Y4G4dK+AH96vAr4tGvxKKaX6h4Ya6wq+ZN+JV+UxunwVHxEDyVmQlHXiPXEYhEd2R637HA1+pZRSfU+PhLy7o513yMekBnfM9l6mwa+UUqp39UTIR0RD0giPgB9pXcVHRHVTpfsvDX6llFLB0xzye90hnweVR7u2z7hBHlfxI6332LSQuooPhAa/UkqpnuF0WCFf/CUUfem+ku9KyAskDLbGu296JY+0ru6V3zT4lVJKdZ0xUFNsBXzRbivsS/LA5ejkDgUShrgDPltDvhtp8CullApcY511Bd90NV+8B2pLO7kz75Af7Q55vR/fEzT4lVJKtc8Yq7NdsftqvuhL61l54+rEzjxD3rO5XkM+WDT4lVJKtdRQDcV73U32e6ygb6jq3L5iUiFtLKSN0yv5PkKDXymlQl1NCRTuhIIvrPeyQ3TqUTpbhBXuaWMhdaz1HpPS7dVVXaPBr5RSocQYqC6Egh1QsNN6rzreuX3FDbKu5NPGQdoYSBzRLbPOqZ6l/4WUUmogMwYqDltX801hX1sS+H7Cozyu5MdY71EJ3V9f1eM0+JVSaiBxuawJawq+ONF0X18Z4E4EEoe6r+TdYZ8wFMLCeqLGKsg0+JVSqj9zNlqP1TVdzRfuDHx+ebFZj9FlTISMCZA+AeyxPVNf1es0+JVSqj9xuaygP7YVjn9u9bx3Nga2D1uEdTWfPgEyTrKu6LWnfcjQ4FdKqb6s6Rn6Y9vcYb8dGmsC20d41ImQzzjJenbeFtEz9VV9nga/Ukr1NXXlcOxzd9hvg5qiwLaPjLcCPt0d9ElZen9eNdPgV0qp3uaot+7NH3U335fuD2z7mFSPoJ9gdcTTmelUGzT4lVIq2FwuKM2zmu6PbYPCXYFNZmOPg8zJkDkFBk2BuAwNeuU3DX6llAoG7/v0DdX+bxsWbl3RZ06xXsnZGvSq0zT4lVKqJzgdVvP94U1wZLMV/H4T6/G6zCkwaLLVMS/c3mNVVaFFg18ppbpLXTkcybWC/ugWaKz1f9vYNMic5g77SToqnuoxGvxKKdVZxljT0x7ZDIc3W7PY+Tu5jT3WuprPnGqFffygHq2qUk00+JVSKhCOBijYbjXhH/7M/0ftwsIhfbz7Pv1U6z69PmKneoEGv1JKdaSmxH1V/5nVOc/Z4N92kQkw5GQYOsMKfHtMz9ZTKT/4FfwiMhx4HLgAEOBd4FZjzMEOtrsPuLeN1fXGmCiPsvuBLB/lLjXGvOFPPZVSqlsYYw2Le3iT1YRfmuf/tskjYcgpMPQUSB2jve9Vn9Nh8ItIDPAeUA9ci3UD60FgtYhMNca090zKc8Byr2Wx7mXLfJRfAdzntWxXR3VUSqkuczqsJvyDn1iBX1fm33a2COtZ+qEzrKv72NQeraZSXeXPFf+NwChgvDFmD4CIbAW+BG4GHmtrQ2NMPpDvuUxEFrmP+7yPTYqMMev9q7pSSnWR0wHHt1lhn78BGqr82y4m1d2Ef4rVQS88smfrqVQ38if4FwDrm0IfwBiTJyIfAQtpJ/jbcC1wHOvqXimlgqs57Ne7w96fgXQEUkdbV/VDT7HGvtcmfNVP+RP8k4ClPpZvBy4L5GAiMgw4B3jCGONrfMqviUgNYAM+Ax7W+/tKqS5zOuDYlhNX9v7MbhceBYOnWUE/5GSISuz5eioVBP4EfwpQ6mN5CZAc4PEWAWH4buZ/E9gA5AGDgFuA10VkkTFmia+dichNwE0AI0aMCLAqSqkBzdlo9cA/uB7yN/oX9vY4GHYqjDjdum9v0wef1MDj7//Vvkak6Ew717eBz4wxW1sdwJgftdi5yOvAeuC3gM/gN8YsBhYD5OTk+DlqhlJqwHI0nAj7wxv9GznPHgfDT4MRZ0DGJA17NeD58394KdZVv7dkfLcE+CQipwETgFv9KW+McYrIq8DvRGSwMeaov8dSSoWQ5rD/2Lqyd9R1vE1kPAxrCvuJGvYqpPjzf/t2rPv83iYCOwI41rWAA3gpgG2aWhX0al4pdYLLZd2z3/9hAGGfYF3ZDz/dGgs/zNbz9VSqD/In+JcBj4rIKGPMPgARGQmcCfzcn4OIiB24EnjbGFPo5zbhWJ0HDxpjApnWSik1UJUdgry1VuDX+tHg2BT2I2Za09pq2CvlV/A/i9XRbqmI3IN19f0AcAj4S1MhEckC9gL3G2Pu99rHxVi3C3x16kNErsJ6NPBt934HAT8EZgBXBfB9lFIDTX0l7P8I8t6Hkr0dl49KdF/Zu5vxdTx8pVroMPiNMdUici7WkL3/xGp+X4U1ZK/naBeC9Rier7+ya7GeAvhvG4fJAzKA32OdINRg9fCfZ4zR5/2VCjVOBxzNhX1r4Mhn4PL19K+HprAfMRPST9KwV6odfvVocY/J/40OyuynjZ7+xpiFHWy7HjjXn7oopQYoY6B0v7sp/yOor2i/fHikdVWffZZe2SsVAO3KqpTqXbVl1j37vLXW3PbtEqtj3qg5Vq/8iKgOyiulvGnwK6WCz9FgTYST977VpG9c7ZePz4TsOdbVfWxaUKqo1EClwa+UCg5joHiPdWV/YF3HY+RHxEDWLCvw08bq2PhKdRMNfqVUz6qvtDrp7X0PKo50UFis8fFHzYGhORBuD0YNlQopGvxKqZ5Rsg92r4ADH1nj5rcncZh1ZT9yNsT4GihUKdVdNPiVUt3H2WgNnbt7hdWs3x57nBX02WdByihtylcqSDT4lVJdV10Ee96FPavafwxPbNYUt6PmwJBTdIx8pXqB/tUppTrHGDj+OexeDvmbaHdKjZg0GHsBjD5H57VXqpdp8CulAtNQY/XM//KdjjvrZU6BcfOsq3sdYEepPkGDXynln7JD8OUK69l7R33b5SKirY564y6EhCHBq59Syi8a/EqptjkdcHij1Zxf8EX7ZROHW1f3I2friHpK9WEa/Eqp1mpLrY56e95tf/pbscHwU2Hshda0t9ozX6k+T4NfKXVCSR58sQwOfdr+jHjRyTD6PBhznj53r1Q/o8GvlLKa8be/YY2b356Mk6yr+2Gn6qN4SvVT+perVKgyxprrfscbULir7XI2uzXIzti5kJwVtOoppXqGBr9SocblskbX27EUyg60XS5+sNUzP/sssMcGr35KqR6lwa9UqHA2Ws/f71gGVcfbLpc2FiZeCkNP0c56Sg1AGvxKDXSNdVbv/J3/bb+H/uBpMPES7Z2v1ACnwa/UQFVfCbuWW8/gN1S1UUhg+GkwcSGkjg5q9ZRSvUODX6mBpqbEurrf827bI+yJDbK/YgW+jq6nVEjR4FdqoKg4aj2Dn/d+28/g2yJgzPkw4WsQmxrc+iml+gQNfqX6u5I8q4f+wfW0OUNeRIw1nO74eTo7nlIhToNfqf6qeC9se9V6Fr8tUUkwYb51lW+PCVrVlFJ9lwa/Uv1NxRHY8i849EnbZWLTrfv32XMg3B68uiml+jwNfqX6i5oS+Pw/sPc9MC7fZRKHw6RLYMRMCLMFtXpKqf5Bg1+pvq6hGr540+qp72z0XUYH3VFK+UmDX6m+ytEAX66wJs9p6zn8lFEw7SrInKKBr5Tyiwa/Un2NywX734etr0BNse8ycYOswB9xhga+UiogfgW/iAwHHgcuAAR4F7jVGHPQj23beL6Ik40xuR7lwoC7gJuBTGAXcL8x5j/+1FGpfs8YOLwZtrwE5fm+y0QlwuRvwuhzdVpcpVSndPgvh4jEAO8B9cC1WA8KPwisFpGpxphqP47zD+AvXst2e31+ALgduBvYBFwJvCoiFxtj3vbjGEr1X4W74LMlUOT9Z+EWHgUTF8D4+RARFdy6KaUGFH8uGW4ERgHjjTF7AERkK/Al1tX5Y37s47AxZn1bK0UkAyv0HzbGPOpevFpExgAPAxr8amAqz4fcl+HwRt/rw8Jh7FyYdClEJQS3bkqpAcmf4F8ArG8KfQBjTJ6IfAQsxL/g78iFgB1Y4rV8CfA3Eck2xuR1w3GU6huqi2HbK7BvLb5H2xNrLP0pl0NcerBrp5QawPwJ/knAUh/LtwOX+Xmc74vIHYATWA/ca4z5wOsY9cAer+22u98nAhr8qv+rr7R66e9e3vZ4+kNOgWlXQnJWUKumlAoN/gR/CuBrEu8SINmP7ZcA/wWOAFnAHcB7InKBMWaNxzHKjDHelz4lHutbEZGbgJsARowY4UdVlOoljnrY9T9rTP3GGt9lUsfA9Kth0MTg1k0pFVL87RbcRlukHxsas8jj4wcishT4HKuD4GyPfQV8DGPMYmAxQE5OTltPDyjVe4yBQ5/C5ufbfjQvYQhM+xYMy9FH85RSPc6f4C/F9xV3Mr5bAtpljKkUkbeA6z0WlwDJIiJeV/3JHuuV6l8qj8HGv8PRXN/ro1Ngyjdh1Nk6vK5SKmj8Cf7tWPfgvU0EdnTyuN5X+NuBSGA0Le/zN7V5dvY4SgWfowF2vGE16/u6jx8RY42nP+4inUBHKRV0/gT/MuBRERlljNkHICIjgTOBnwd6QBFJAOYDnlOLLQcagKuBX3ssvwb4XHv0q37j8GbY9HeoKmi9Tmww/iLr0bzIuODXTSml8C/4nwVuAZaKyD1YV+oPAIfwGJRHRLKAvVij7d3vXnY7MB5YzYnOfbdjjcx3ddO2xpgCEXkc+IWIVAKbgSuAc7EeGVSqb6sqhM3/gPw2nsfPOAlOvQEShwW1Wkop5a3D4DfGVIvIuVhD9v4Tq5l+FdaQvZ4zhwhgA8I8lu0CLnW/EoEK4CPgemPMp16HuhuoAn7CiSF7LzfGvNmJ76VUcDgdsPNNa7pcXzPnRSXCyYtg5GztuKeU6hOk9RN0/VNOTo7ZuLGNqy2lesKxbbDhr1B51MdKgXEXwtTLwR4b9KoppQY+EdlkjMkJdDud5UOpQNWUwGf/hAPrfK9PHWM166dkB7deSinlBw1+pfzlcloj7m19BRx1rdfb4+Dkq2HUOdqsr5TqszT4lfJH4S7Y8ByUtTET9ehzYdpVOpGOUqrP0+BXqj115fDZi5C31vf65JFWs37a2KBWSymlOkuDXylfXC7Yuwq2vAwN1a3XR0TD1CutKXPDwlqvV0qpPkqDXylvxXut3vole32vH/kV615+tD9zVCmlVN+iwa9UE0cDbP037HwLn3NGJQyFU6+HQb5GsFZKqf5Bg18psK7yP34KKg63XhceCZO/CeO/Cjb9k1FK9W/6r5gKbU4HbH8Ntr8OxtV6/fDT4JTvQGxq0KumlFI9QYNfha7SA7D+z1C6v/W6mFQ47UYYcnLQq6WUUj1Jg1+FHpcLvlgG2171PW3uqLPhlG/rULtKqQFJg1+Flooj1r384j2t10Ulwek3wdAZQa+WUkoFiwa/Cg3GWMPt5r7oexa9rFmQcx1Exge/bkopFUQa/GrgqyqET56G49tbr7PHWSPvZc0Mfr2UUqoXaPCrgcsY2LcaNj3ve1KdoTlw2g06EI9SKqRo8KuBqaYEPl0MRz5rvS4iGmZ8B7Ln6Cx6SqmQo8GvBhZj4MBHsPFvvsfYz5wCp39fn8tXSoUsDX41cNSVW1PnHvq09TqbHU5eBGMv0Kt8pVRI0+BXA8OhDVbTfn1F63Xp4+GMH0B8ZvDrpZRSfYwGv+rfGqph499h/wet14WFw7SrrDH2depcpZQCNPhVf3Z0K6x/GmpLWq9LGQ0zfwCJw4JfL6WU6sM0+FX/43TAtldgx9LW68QGk78OEy/RmfSUUsoH/ZdR9S/VRfDRH6Bod+t1SSOse/kp2cGvl1JK9RMa/Kr/yN9ozabX6jE9gYkLYco3wRbRK1VTSqn+QoNf9X3ORvhsiTXWvrfoZJj1Yxg0Mfj1UkqpfkiDX/VtlcfgwyegNK/1uiEnwxnfh6jEoFdLKaX6Kw1+1Xft/8h6Nt97nH2xwfSrYMLFOhiPUkoFSINf9T2Oetj0D9j7Xut1sWlw5q2QNjbYtVJKqQHBr1FNRGS4iPyfiJSLSIWIvCYiI/zYLkdEFovIThGpEZGDIvKiiLTqdi0i+0XE+Hhd0onvpfqrskOw4pe+Q3/46XDRIxr6SinVBR1e8YtIDPAeUA9cCxjgQWC1iEw1xviYCaXZlcAk4I/AdmAo8P+AjSIy3RhzyKv8CuA+r2W7/Pgeqr9rmkJ349+sznyewsLhlGt1nH2llOoG/jT13wiMAsYbY/YAiMhW4EvgZuCxdrb9nTGm0HOBiHwE5Ln3+yuv8kXGmPV+1l0NFA01sOFZOLCu9br4wTD7VkgeGexaKaXUgORP8C8A1jeFPoAxJs8d4AtpJ/i9Q9+97ICIFGJd/atQV7zXGpCn6njrddlnQc71EBEV/HoppdQA5c89/knA5z6WbwcCfnhaRE4CMoAvfKz+mrsvQL2IrNf7+wOYMbDrf7DyV61D32a3HtOb+UMNfaWU6mb+XPGnAKU+lpcAyYEcTETCgWeAQuCvXqvfBDZg3QYYBNwCvC4ii4wxS9rY303ATQAjRnTY11D1FfWVsP4ZOLyx9bqkEVav/URtEFJKqZ7g7+N8xseyzvSyehKYBcw3xrQ4mTDG/KjFzkVeB9YDvwV8Br8xZjGwGCAnJ8dXHVVfU7jLatqvKW69bsz5Vie+cHvw66WUUiHCn+Avxbrq95aM75YAn0Tkt1hX59caY97pqLwxxikirwK/E5HBxpij/h5L9UHGwI43YOsrYFwt10VEw2k3Q9bMXqmaUkqFEn+CfzvWfX5vE4Ed/hxERO4Gfg782BjzT/+r19yqoFfz/VldOax7Eo5tbb0uZTSc+ROIHxT8eimlVAjyJ/iXAY+KyChjzD4AERkJnIkV5u0SkR9jPfd/tzHmT/5WzN0f4DLgoDHmmL/bqT7m+A5Y90eo9dE4NOFimHYV2HQASaWUChZ//sV9Fquj3VIRuQfr6vsB4BDwl6ZCIpIF7AXuN8bc7152JfAEsBx4T0TO8NhvhTFmh7vcVViPBr7t3u8g4IfADOCqLnw/1VuMge2vwdZXadVgY4+DmT+AoTN6pWpKKRXKOgx+Y0y1iJwLPA78E6v5fRVwqzGmyqOoADZaPiI4z718nvvlaS1wtvvnPKxH/H6P1Z+gBquH/zxjzIrAvpLqdXXlsO5PcGxbi8VOY7BlnGRNoxub2kuVU0qp0CbGDIzb5zk5OWbjRh+Ph6mgaHQ2UlpfSumRTZR+9jxl9WWUGielOCjBSZlxUhefQURSFrERccRExBAbEUtsRCzR4dHNP8eGxxITEdO8PibcerfbtKe/Ukp5EpFNxpicQLfTm6uqQw6Xg/L6civM60qtgK+zXmX1ZZTUlVDdWAWVx6DiSOsdhIVDyiiISqTR5aCsvoyy+rKA6hARFmGdELhPBDxPDBIjE0mLTiMtOo3U6FQibZHd88WVUmoA0uBXLZTVlbH+6HoOVx1uDvmK+gpMew9WuBqhZD/UV7ReZ4+DlGxrNL4uaHQ1Ul5fTnl9eYdl4+3xLU4E0qLTSIuyPidGJiI60Y9SKoRp8CsAjlcf592D7/Lp0U9xGqf/G9ZXQmle6xn1AOIzIX5I0GfUq2yopLKhkrzyvFbrwsPCSY1KbXlS4PGzthYopQY6Df4Qd6DiACsPrGRLwZb2r+pbMS2a9gVIxEayhJMUHkPKmLkkZ04nKTKJ5KhkkqOSiY+Ip95ZT3VjNTWNNVQ73O+N1VQ7qqltrLV+dn+uaaxpXh/QyUg7HC4Hx2uOc7zGx6RAQJw9rrl1ID0mnfTodNJj0smIySA2IrZb6qCUUr1Jgz8EGWPYVbqLd/a/w+7S3X5tE2ePIznSCvBkWxQpeR+TVFVPsm0QydhIwIZNBNInWAPyxPga7BGiwqOICo8iNdr/Xv3GGBpcDdQ01lDVWGWdEDisE4KqxipK60opqi2iuLaY4rpiutJhtaqhiqqGKvZX7G+1LiY8hrSYNDKiM5rf9aRAKdXfaPCHEJdxsaVwC+/sf4dDlYfaLJcWncacYXMYGj+U5MhkkiKTiLBFWCs9B+QJ8wq7SZfClMsgzNat9RYRIm2RRNoiSY5qf14op8tJaX0pxbXFFNUWNb+K66zPNY01na5HjaOGgxUHOVhxsNW69k4KYsJjtF+BUqrP0OAPAY2uRjYc3cDKgysprClss9yw+GFckHUB09OnY/MOb2Ng++vWWPvetwQi460pdIec3P2VD5AtzNZ8334841utr2mssU4K6opanhjUFlNSV4LLex4BP7V3UhAdHt1826Cpbk2vpMgkPSlQSgWVBv8AVuuoZd3hdbx36L12e8OPTR7LBVkXcFLKSb5DqL2x9tPHw6yf9JsBeZrGCBieMLzVOqfLSVl9GUW1RRTWFlJYU9j8XlRbRKPLRwdGP9Q6ats8KbCJrbljoa9Oh1HhUZ06plJKtUWDfwCqbKhkzaE1vJ//PrWO2jbLTU2fygVZF5CdmN32ztoba3/iJTD18m5v2u8ttjArhFOjU1u1FhhjKKsvo6CmgKLaolbvnT0pcBonBTUFFNQU+Fzf1Nmw+eTA46QgKTKJMAnzuZ1SSrVFg38AKaot4r2D7/HxkY/bDKIwCePUzFO5IOsCMmMz295ZP2naDxYRaX46IVgnBdB+Z0Ob2EiJTiEl6sQrOSqZ5MhkUqJSSIpKIiIsotPHVkoNTBr8A8DhqsOs3L+STQWb2uzRbrfZOXPImZw74twOO8gNpKb9YPDnpKCwtpCimqJWfQu60tnQaZzW7Yh2+m0k2BNIjkpueWIQlUxKpPVzbESs9jFQKsRo8PdjByoO8Pa+t9levL3NMjERMZw97GzmDJ/j3yNnR3Lhk2dComk/GDxPCsYlj2u1vqaxhuK6YutRxKYnEdwnByW1JV0ev6CioYKKhgoOVBzwud5us594TNN9gtDUapAUmURiZKL2M1BqgNHg74dqGmtYuncp6w6va3PQnaTIJM4bcR6zhs7ybzS62lLY/AIcWNd6nT0OZt0SUk37wdLc2TC+dWdDl3E1dzZsevLAs7WgurG6y8dvcDa0O6ARWCcHTScBSZFJJNgTmn9OjEy0XvbEE498KqX6NA3+fsQYw6fHPuX1Pa9T1VDls0xmbCYXZF3AjEEzCA/z4z+vMbBnFeS+CL6anbVpv9eESVhzE72v1oJaRy3FtcXNcyoU1xY3T6BUUlfS8RwLfmpwNrTbAbFJTHgMCZFeJwX2xBMnB5GJxNvjtd+BUr1Mg7+fOFp1lH/v+jd7yvb4XD8yYSRzR85lStoU/+/Zlh2ET5+FIl+j9wlMXABTr9Cm/T4qOjyaYfHDGBY/zOf6plkVS+pKrJOB+pLmWRVL6kooqSuhwdnQbfWpcVgjKh6rPtZhvePsccTb40mwJxAfEd/ic1yE9XO8PZ7o8Gjtg6BUN9Pg7+PqnfUsz1vOqoOrfA4ukx6TzmXjLmv7GXxfHA3w+X/gizfB1z3k5JFw2k2QOrprlVe9KjwsvPnxRF+MMdQ4appPDJpOCErrSymrs6ZOLq8v77Z5EprUOmqpddS22ymxiU1sxNnjrBMEe3zzSUHzZ3sc8RHxxETEEGePwx5m1xMFpTqgwd+HbS3cyqu7X6W0rnVHu/CwcC4ceSHnjzg/sHurR7fChmehykezbXikdYU/bp5e5YcAESE2IpbYiFiffQzgxMlB00lAeX055Q3llNe53z2WdWWOhLY4jdPv6ZjhxIlC0/dqesVFxBETEdP8s+c6HVJZhRoN/j6ouLaYV3e/yudFn/tcPzF1IpeNu4z0mHT/d1pXDpuehwMf+V4/dAbkXAexaZ2osRqoPE8OhsYNbbOcy7ioaqw6cSLgcUJQVl9GRX0F5fXlVDZW9sgJQpNATxQABCE6Ipq4iLgWJwgx4VbHy+jw6BY/N32OjojWFgbVL2nw9yEOl4NVB1exPG+5z0FfEiMTuWzcZUxLn+b/PzbGwN73rM57DT56gUcnw4zvwvDTQP8BU50UJmEk2BNIsCe02XoAVgtC06yK5fXlVDVWUdlQ2frVaL13Zx+ENuuEaZ4CuoD2OzB6s4mN6Aj3iUB4dOuTA/eyphOHKJs1O2V0eDRR4VF64qB6hQZ/H7G7dDev7HrFZ8coEeGc4efw1eyvBvZMdXk+fLoYCnf5WCkwbi5MvRLsMZ2vuFIBEBHi7HHE2ePaHznSrd5Z3+KEoKqhiorGiuafKxsqqWqsorqxmurGahwuRxC+xQlO42weXbEzBGmeqjrKduKEoNVnrxMGPYFQXaHB38sqGip4/cvX2XBsg8/12YnZXDnhynabWVtxNFjD7X6xDHz9Q5iUZXXeSxvTyVorFRyRtkgioyNJi+74FpQxhgZXA1UNVdQ01lDtqLZ+dtQ0tzJUN1RT7ahuPlGobqymzlEXhG/SRp0xzZ0du0IQImwRRNmisNvsRIVHNU9lHWmLtE4ObPYT621RRIZb6+w2O9G2aOw2e4uyejIxcGnw9xKXcfHh4Q95c++bPv/oYyJiuGTMJcwcPDOwP76jW2HDc1DlY0AWm90aeW/8V7XznhpwRKT5RKGtJxl8cbgcVDdWU9NY09x6UNNYQ62j1npEsdF6TLHWUdtqeXc/8dBZBkODs6Hbb43Ybfbmk4CmnyNtka0+R4RFNJ9kRNgiTpTz2LapnN1mJyIsgoiwiNbTf6ug0ODvBQcrDvKvXf/yOU0rwMwhM1k4eiFx9jj/d1pXDpv/Cfs/8L1+yMmQcz3EBdAhUKkQEB4W3jzAUKAanY0tTgp8nSB4fq511FLnrKPWUUu9o75LEzgFQ0+cTHiyic06UQizE2GLaD4h8Dw5aDrBaPrcVN57WfPPYRGEh4Vjt9kJDws/8TnM+qytGBr8QVXTWMN/9/2XD/I/8Dmi2pC4IVw5/kpGJY3yf6fGwL7V8NmL4Os+Y3QyzPgODD9dO+8p1c0ibBEk2jp30gBWa0O9s/7ESYGjrvlznaOOOmcddY66Vp+bTh5qHbU9Hs49yWmcOB1O6gje7Zamkw3PkwTvE4jwsHDCJAyb2AiTsBavpmWe60SkxXLPz97bxkTEMCl1UtC+ry8a/EFQXl/OtqJtvL3vbSoaKlqtt9vszB81nznD5vg3zC6AsxEKvoDtr1nvrQiMPR+mXQV2PybnUUoFXXhYOOFh4f5NoNUOl3FR76ynwdlAnaPOendaJxH1znrqHfXUOT2WO6zlDa6G5pOHemd984lHg7Ohz7dGdFZvnGx4Ghw7WIN/IKpz1LGnbA9flHzBrpJd7Q5hOj1jOt8Y+42Op8oFqDwGR7dYM+gd/xzaOstPGuHuvDe2c19AKdWvhElY82OEnW198OYyLqs1wdXQfGLQ1LrQ/Nl14rN3uUZXY/OJh+d2ja5G6+Vs7Ja5JPqbvtCvQYO/GzhcDg5WHGRnyU52lu5kf/l+n8PrekqLTuOy8Ze1f+bnqIfj2+ForhX2vjrsebJFwJTLYPx8sOl/WqVU54VJmPXoID0zLbMxBodx0Oi0TgSaTh6aP3v+7D5h8C7TtMzhctDoamx+bzqx8FzW9N7bwiSst6vgX/CLyHDgceACQIB3gVuNMb57p7XcNgp4ALgGSAJygbuMMe97lQsD7gJuBjKBXcD9xpj/+PldgsYYw/Ga41bQl+zky9IvqXfW+7WtTWxcMPICLsy6sPVQu8ZAxWEr5I/mWk34/j6XPHg6nHo9xGUE8lWUUqpXiAgREhHU2RqNMS1ODlqcKLhPFppORlzGdeKFC5fLhdM4cRnr3RhjffZa1942LuPy69HUntZh8ItIDPAeUA9cCxjgQWC1iEw1xnQ0KfhfgfnAHcA+4IfAChGZaYzJ9Sj3AHA7cDewCbgSeFVELjbGvB3Qt+oB5fXl7C7d3Rz2AQ0JKsLIhJGMTxnP6Zmntxxqt6EGjm2zmvCP5kJNsf+Vik62An/E6da7dt5TSqk2iVjjHQQ0v8kA5M8V/43AKGC8MWYPgIhsBb7Eujp/rK0NRWQa8C3gOmPM393L1gLbgfuBBe5lGVih/7Ax5lH35qtFZAzwMBD04K931rOndA87S3ayq3QXR6qOBLR9RkwGE1ImMCFlAmOSxhAT4R4dzxgoybNC/ugWKNzte4Y8X8LCIX28FfKDp1n38jXslVJKBcCf4F8ArG8KfQBjTJ6IfAQspJ3gd2/bCPzbY1uHiPwL+LmIRBpj6oELATuwxGv7JcDfRCTbGJPn1zfqgqKaIjbu+x87y/eSV5mP0zisoLYqjtXY4fXZABjibJFMiB3K+OjBTIgZTHJ4FLhcULAPju8B47Lu0R/dYj1z76/YdBgyHQafDIMmQUTP3G9TSikVGvwJ/knAUh/LtwOX+bFtnjGmxse2dmCM++dJWLcS9vgoBzAR6PngrzrKfzf/2a+yEQhjJZLxEsUEiWIITqRsD62/QoBsEZAxyR320yB+sF7VK6WU6jb+BH8K0HpCeCgBOnoGrb1tm9Y3vZeZ1vN1epfrUaOTxxKO4PDxiIkAI8TOBIlivESRTSQR3RXICUOtkB8yHdJPgnB79+xXKaWU8uLvM1++Hrb0J/XEz239LddypchNwE3uj1Ui4msaOn+lAUVd2F61pr/TnqG/1+6nv9Pup7/TnuH5e83qzA78Cf5SfF9xJ+P7at5TCTCijW2b1je9J4uIeF31e5drwRizGFjcQR38IiIbjTE53bEvZdHfac/Q32v3099p99Pfac/ojt+rPyMJNN2D9zYR2OHHttnuRwK9t23gxA3x7UAkMNpHOfw4jlJKKaX84E/wLwPOEJHmmWNEZCRwpntdR9tG4NEJUETCgSuAd9w9+gGWY50IXO21/TXA58Ho0a+UUkqFAn+a+p8FbgGWisg9WPfiHwAOAX9pKiQiWcBerNH27gcwxuSKyL+BJ0QkAqtn/veBbDxC3hhTICKPA78QkUpgM9bJwblYjwwGQ7fcMlAt6O+0Z+jvtfvp77T76e+0Z3T59yqtO9L7KCQygpZD9q7CGrJ3v0eZkVjB/mtjzH0ey6OB32AN5JMEbMEasneN1zFswC+wBgzyHLL3/zr53ZRSSinlxa/gV0oppdTA0PvTBPUiERkuIv8nIuUiUiEir7lbN1QniMg3ReQ/InJARGpFZJeI/FZE4nu7bgOJiCwXESMiD/Z2Xfo7EfmqiLwvIlXufwM2isi5vV2v/kxEzhSRd0SkwP073Swi1/V2vfoLERkmIn8SkY9FpMb9tz7SR7lkEXlORIpEpFpE3hWRKf4cI2SD32PyoQlYkw8tAsZizREQ25t168duB5zAL4F5wNNYfTpWumdfVF0kIlcB03q7HgOBiNyMNSrpJuBSrE7IrwLeTyEpP4nIVKzZWyOwbtt+A9gA/FVEvt+bdetHxgCXYz0u/4GvAiIiWJ3n5wE/wvo9R2Dl17AOj2CMCckX8BOskBrjsSwbcAA/7e369ccXkO5j2bexOoSe29v16+8vrD4yx4Cr3L/TB3u7Tv31BYwEarH6KvV6fQbKC3gI6wmtOK/l64GPe7t+/eEFhHn8fIP7b32kV5mF7uXneCxLxBrz5o8dHSOUr8J8Tj4ENE0+pAJkjCn0sXiD+31oMOsyQD0CbDfGvNzbFRkArgNcwDO9XZEBxo41MVut1/IyQriFORDGGJcfxRYAR4wxqz22KwfexI/8CuX/EJOAz30s386JgYNU181xv3/Rq7Xo50RkNlbryQ96uy4DxGxgJ3CliOwVEYeI7BGRH/Z2xfq5f7jf/ygiQ0QkSURuBM7DejJMdY/28muEiMS1t7G/Y/UPRF2ZfEj5QUSGAvcD7xpjNvZ2ffor9xgYfwEeNcZ0ZT4KdcIQ9+v3WH1S9mLd439SRMKNMX/ozcr1V8aYz0XkbOB1TpykNgLfM8b8q7fqNQClAPt9LG8a3j4ZqGpr41AOfuj85EOqA+4zzqVYfSa+28vV6e/uAprGw1DdIwyIB75jjHnNvew9d+/pX4jIH437xqnyn4iMBf6DdeX5Pawm/4XAMyJSZ4x5sTfrN4B0amK7JqEc/F2ZfEi1Q0SisHqcjgLmGGPye7lK/Zb78dK7sTr5RIpIpMfqSBFJAiqNMc7eqF8/Voz1FM9Kr+XvYPWUHgwcCXalBoCHsK7wLzbGNLqXrRKRVOAPIvKyn/ewVftKaDu/oIMMC+V7/F2ZfEi1wd0s/R/gNOCrxphtvVyl/m4UEAUswfpjbnqB9fhkKeDXs7uqhe1tLG+6YtJw6pwpwBaP0G/yKZAKZAS/SgNSe/l10BjTZjM/hHbwd2XyIeWD+1n9F7E68iw0xqzv5SoNBLnAOT5eYJ0MnMOJWS6V/153v1/otfxCIN8YcyzI9RkojgHTRcTutfx0oI42plhXAVsGDBWRps7TiEgC8DX8yK9Qbur3a/IhFZCnsDpI/QaoFpEzPNbla5N/4IwxZcAa7+XW+B0cMF5zXii/vQ2sBv4iImnAPuCbwFy0T0pXPIk1CNKbIvJnrHv8C7DGnnjcGNPQm5XrL0Tkm+4fZ7jfLxKRQqDQGLMWK9w/BpaIyB1YLX+/wGqxeqTD/Ydy/xV/Jh9S/hOR/UBWG6tbTN6kukZEDPAbY8w9vV2X/sp9hfRbrMBPxnq872FjzEu9WrF+TkQuwuqQOgnrNtVerBnl/qJ9Ufzj/vv2Za0x5mx3mRTgUeASrN/zx1iDz23pcP+hHPxKKaVUqAnle/xKKaVUyNHgV0oppUKIBr9SSikVQjT4lVJKqRCiwa+UUkqFEA1+pZRSKoRo8CvVC0TE+PHa39v1DDYRGSki93mOqKmU6l6hPHKfUr1pptfn14EtwH0ey+qDVpu+YyRwL/Ah1mh6SqlupsGvVC/wnsdAROqBooE4v4GIRBpjevUkpi/UQam+Qpv6leqjRCRbRF4UkUIRqReRXBG51KvMfe7bAhNEZIWIVIvIQRH5rnv9IhHZKSJVIrJaREZ7bb9fRJaIyI0iskdE6kRks4icgxcRmSMiq0Sk0n2cFSIy2avMGhH5UES+JiKfuU9ofuBed4uIfCwiJSJSJiLrRWS+x7ZnY42fD7DS45bH2e71RkTu8zreSPfy73gs+4eI5IvITBFZJyK1uMcvF5E0EXlaRA67f6c7ReQm//+rKNX/afAr1QeJyHDgE2AacBvWRCebgf+IyAIfm7wKvIU1bvcm4G8i8hDwfeDnWBPPjAd8jUM/B/gpcDdwJdYthv+JyHiP+szHmsuiCrgG+BYQD3zgrqunccAfgT9hzXa3yr18JPAc1kROVwAbgf+6x3bH/f1+6P75x1i3Q2a6lwcqEfgX8DJwEfCSe2z+j4D5WLdU5gNvAk+LyI86cQyl+iVt6leqb7oPa+KoOcaYYveyFe6QvZ/WU2/+3hjzAoCIbMSanvNmINsYU+FePhj4g4hkGWMOeGw7CDjTGHPQXW4VcAC4B1jkLvMHrAlCFjZtJCKrse7D/wy41WN/acBcY0yuZwWNMbd7bBuGdUIwDvge8D9jTIWI7HAX+aKLtz3igGuMMUs9jvn/sCaRmmKM+dK9+F0RSQLuFZGnjTGOLhxTqX5Br/iV6pvmYU0dWy4i4U0vYAUwzX316ul/TT8YY0qBAmB9U+i77XS/e1+hr28Kfff2lVitBzMBRGQsMBp40asuNVgzgp3ltb/93qHv3s8MEfmviBwHHEAj1syY473LdgMH8F+vZfOwWlHyfPxOU4GJPVAPpfocveJXqm/KAL7tfvmSCniGeqnX+oY2loE1haen4z72fxwY6lEXgL+6X94Oen0+6l3A3VKxCtgB/Mi9jQN4ADjJxz67qsDHFLAZwBisEw5fUnugHkr1ORr8SvVNxcAHwO/aWH+kG481qI1lhz3qAvAL4F0fZRu8Pvua63se1n33y40x+U0LRSQmgHrWA3avZW2Fta86FGO1hPykjW12BVAXpfotDX6l+qblWE3t240xtT18rDNEZLgx5hCAiMRjdXx7y71+F7AfmGSMebiTx2gK+OarbREZB5wJ5HuUa3rkLtrHPg4Ak72WzfdRri3Lcbc2GGMKAthOqQFFg1+pvulXwKfA+yLyJFbwJmMF3yhjzHXdeKzjwDvuR+XqgbuAWKxmeIwxRkR+CCwVETvwClCE1SowCytIH+vgGO9iNe2/ICL/HzAY+DVWk79nX6Pd7nLXiUiJuz673P0O/gXcIyJ3A+uBrwBXBfA9H8d6muADEXkc64QmFpgAfMWz46JSA5l27lOqD3J3tsvBGs3vIWAl8DTWo3fvdfPh1gL/n/s4/8bqA3CRMWa3R33exurEF4v1SN4KrGfjM7E6+LXLGLMduBqrV/0y4E6sxwzf9ypXDNyC9RjjWmADMMO9+rfAk+71b2D1DViEn4wx5VgnKm9jndysAP4GLOTE+AFKDXhijK9bYUqpUOCeD+BDY8w1vV0XpVRw6BW/UkopFUI0+JVSSqkQok39SimlVAjRK36llFIqhGjwK6WUUiFEg18ppZQKIRr8SimlVAjR4FdKKaVCiAa/UkopFUL+f1q7xh88vpERAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.constants import h, c, k, N_A # import physical constants\n",
    "\n",
    "def q(T): return 1+np.exp(-epsilon/(k*T))\n",
    "def energy(T): return epsilon / (np.exp(epsilon/(k*T)+1) ) # atkins p625 (10th)\n",
    "def Cv(T): return np.gradient( energy(T)/k, T[1]-T[0] )\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "epsilon=k*4\n",
    "T=np.arange(0.1, 10, 0.5)\n",
    "plt.plot(T,q(T), label=r'$q = 1+1 / e^{-\\beta \\epsilon}$', lw=4, alpha=0.7)\n",
    "plt.plot(T, energy(T)/k, label=r'$\\langle \\epsilon \\rangle = -\\frac{1}{q}\\frac{\\partial q}{\\partial \\beta}$',\n",
    "         lw=4, alpha=0.7)\n",
    "plt.plot(T, Cv(T), label=r'$C_v=\\partial \\langle \\epsilon \\rangle / \\partial T$', lw=4, alpha=0.7)\n",
    "plt.legend(frameon=False, loc=0, fontsize='large', ncol=2)\n",
    "plt.ylim(0, 2)\n",
    "plt.xlabel('Temperature');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bragg-Williams Model for a Two-component System\n",
    "\n",
    "This example demonstrates how the Bragg-Williams (BW) model predicts the partial pressure of compoments A and B over a liquid mixture of the two. The interaction parameter $\\alpha$ includes the relative A-B neighbor interaction, compared to the A-B and B-B interactions, as well as the number of neighbors, $z$. The partial pressure is given by,\n",
    "$$\n",
    "p_A = p_A^0 \\chi_A e^{\\alpha \\chi_B^2}\n",
    "$$\n",
    "\n",
    "where $p_A^0$ is the vapor pressure of pure A. If either $\\alpha$ or the mole fraction, $\\chi_B$, of the other component approach zero, ideal behavior, i.e. Rault's law, is obtained:\n",
    "$$\n",
    "p_A = p_A^0 \\chi_A\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6deef06586d940f993f02b52b515aaa9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.0, description='alpha', max=3.0, min=-3.0, step=0.5), FloatSlider(va…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "chiA = np.arange(0,1.1,0.05) # mole fraction of A\n",
    "chiB = 1-chiA                # mole fraction of B\n",
    "\n",
    "def plotbw(alpha=1, p0A=0.7, p0B=0.2):\n",
    "    plt.plot( chiA, chiA * p0A * np.exp(alpha*chiB**2),'b-', label='Bragg-Williams'  )\n",
    "    plt.plot( chiA, chiA * p0A, 'b--', label=\"Rault's law\")    \n",
    "    plt.plot( chiA, chiB * p0B * np.exp(alpha*chiA**2),'r-')\n",
    "    plt.plot( chiA, chiB * p0B, 'r--')\n",
    "    plt.xlim(0,1)\n",
    "\n",
    "    plt.xlabel(r'$\\chi_A=1-\\chi_B$')\n",
    "    plt.ylabel(r'$p$', color='k')\n",
    "    \n",
    "    plt.annotate(r'$p_A^0$', (1,p0A), color='blue')\n",
    "    plt.annotate(r'$p_B^0$', (0,p0B), color='red')\n",
    "    plt.legend(loc=0, frameon=False)\n",
    "\n",
    "i = interact(plotbw, alpha=(-3,3, 0.5), p0A=(0,1,0.1), p0B=(0,1,0.1), continuous_update=False )\n",
    "\n",
    "i.widget.children[0].description=r'$\\alpha$'\n",
    "i.widget.children[1].description=r'$p_A^0$'\n",
    "i.widget.children[2].description=r'$p_B^0$'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Equation of State using Generalized van der Waals Theory\n",
    "\n",
    "This notebook illustrates the generalized van der Waals theory (gvdW) for the equation of state for interacting particles. Based on the lecture notes, _Properties of Molecular Fluids in Equilibrium_ by Sture Nordholm.\n",
    "\n",
    "### Pair potential\n",
    "\n",
    "The particles are here assumed to interact via a _Lennard-Jones_ potential,\n",
    "$$\n",
    "\\beta w(r) =\n",
    "4\\beta \\varepsilon_{LJ} \\left ( \\left ( \\frac{\\sigma}{r}\\right )^{12} - \\left ( \\frac{\\sigma}{r}\\right )^{6}\\right )\n",
    "$$\n",
    "\n",
    "Any potential may in principle be given and must return the energy in units of $k_BT=\\beta^{-1}$.\n",
    "\n",
    "### Interaction parameter\n",
    "\n",
    "Here we integrate the above pair potential to get the average interaction energy per particle, assuming that the pair correlation function, $g(r)$, can be described by a simple step function, zero when $r<\\sigma$, unity otherwise:\n",
    "$$ \n",
    "\\hat{a} = -\\frac{1}{2} \\int_{\\sigma}^{\\infty} 4\\pi w(r) r^2 dr\n",
    "$$\n",
    "\n",
    "In this Notebook we simply do the integration numerically so that we can use arbitrary pair potentials.\n",
    "From this we calculate the pressure, $p$, versus density, $n$, using,\n",
    "$$\n",
    "\\beta p_{gvdW} = \\frac{1}{v-v_0} - \\frac{\\hat{a}}{v^2}\n",
    "$$\n",
    "\n",
    "where $v=1/n$ and $v_0=2\\pi\\sigma^3/3$ is the particle volume.\n",
    "For reference we'll also plot EOS for an ideal system (van 't Hoff), $\\beta p_{ideal}=n$,\n",
    "where $\\beta = 1/k_BT$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2d0fbe6ae76543b7b15921be73aa2b84",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='eps', max=2.0), FloatSlider(value=4.0, description='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def pairpotential(r, eps, sigma): # Lennard-Jones potential\n",
    "    return 4 * eps * ( (sigma/r)**12 - (sigma/r)**6 )\n",
    "\n",
    "def Pgvdw(n, eps, sigma):\n",
    "    a  = -2*pi*quad(lambda r: pairpotential(r, eps, sigma)*r**2, sigma, np.infty, limit=20)[0]\n",
    "    v0 = 2*pi*sigma**3 / 3\n",
    "    v  = 1 / n\n",
    "    return 1/(v-v0) - a/v**2\n",
    "\n",
    "def Pideal(n): # ideal gas law\n",
    "    return n\n",
    "\n",
    "def plot( eps=0, sigma=4.0 ):\n",
    "    plt.plot(n, Pideal(n),  'k--', label='ideal', lw=2)\n",
    "    plt.plot(n, Pgvdw(n, eps, sigma),  'r-', label='gvdW', lw=2)\n",
    "    plt.legend(loc=0, frameon=False)\n",
    "    \n",
    "    plt.xlabel(r'Density, $n$')\n",
    "    plt.ylabel(r'Pressure, $\\beta p$')\n",
    "\n",
    "n = np.linspace(1e-7, 6e-3, 50)    \n",
    "i = interact(plot, eps=(0.0, 2.0, 0.1), sigma=(0, 5, 0.1), continuous_update=False )\n",
    " \n",
    "i.widget.children[0].description=r'$\\beta\\varepsilon_{LJ}$'\n",
    "i.widget.children[1].description=r'$\\sigma_{LJ}$'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variation of Dill problem 5.1 (page 91)\n",
    "\n",
    "Let's assume that each corner (n, s, e, w) are associated with an energy, $\\epsilon_i$. Specifically we set an attraction for north ($\\beta\\epsilon_n=-0.5$) and a repulsion for south ($\\beta\\epsilon_s = 0.5$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p_i = [0.38745562 0.14253696 0.23500371 0.23500371]\n",
      "sum(p_i) = 1.0\n",
      "S/k = 1.326\n"
     ]
    }
   ],
   "source": [
    "energy = np.array([-0.5, 0.5, 0.0, 0.0])\n",
    "p = np.exp(-energy) / np.exp(-energy).sum()\n",
    "entropy = -(p*np.log(p)).sum()\n",
    "print(\"p_i =\", p)\n",
    "print(\"sum(p_i) =\", p.sum())\n",
    "print(\"S/k = {:.3f}\".format(entropy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.5"
  },
  "widgets": {
   "state": {
    "87180b4c08944079bd38f6e9c5fe21a8": {
     "views": [
      {
       "cell_index": 15
      }
     ]
    },
    "fbba64368cd84ce3a14a3022b2eb8f39": {
     "views": [
      {
       "cell_index": 17
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}