P05-DemonAlgorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Demon Algorithm\n",
    "===========\n",
    "\n",
    "There are a number of approaches to complex problems involving large numbers of interactions where the objective is to find the \"average\" behavior of the system over a long period of time. We've seen that we can integrage Newton's 2nd Law to see the precise behavior of a multipartical system over time. When we have a handful of objects in a system this works well. However, if we have thousands or millions of particles, it's not practical. Looking at \"average\" behavior however glosses over the precision of following each interaction and attempts only to see what happens on a less fine-grained scale. This means we sacrifice the hope of getting a detailed pictured of a microscopic physical process, but achieve the reward of a more general understanding of the large scale consequences of that process. The demon algorithm is such an approach. It's a simple way to simulate the random exchange of energy between components of a system over time. Here's the basic idea:\n",
    "\n",
    "* Suppose we have a demon..\n",
    "\n",
    "    1 Make a small change to the system.\n",
    "\n",
    "    2 Compute $\\Delta E$. If $\\Delta E<0$ give it to the “demon” and accept the change.\n",
    "\n",
    "    3 If $\\Delta E>0$ and the demon has that much energy available, accept the change and take the energy from the demon.\n",
    "    \n",
    "    4 If the demon doesn’t have that much energy, then reject the change.\n",
    "    \n",
    "Example Problem\n",
    "---------------\n",
    "\n",
    "Compute the height distribution of nitrogen molecules near the Earth's surface. Assume T=const. and that the weight of a molecule is constant.\n",
    "\n",
    "$$ PE(y) = m g y $$\n",
    "\n",
    "so $\\Delta E$ is just $m g \\Delta y$.\n",
    "\n",
    "Below is a sample program that uses the demon algorithm to approach this problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as pl\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9420869702253908\n",
      "(array([ 15.,  51., 148., 256., 266., 175.,  68.,  18.,   2.,   1.]), array([-2.99631416, -2.26929282, -1.54227147, -0.81525012, -0.08822877,\n",
      "        0.63879258,  1.36581393,  2.09283528,  2.81985663,  3.54687798,\n",
      "        4.27389932]), <a list of 10 Patch objects>)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADjZJREFUeJzt3X+o3fV9x/Hna9atozqqeBdsEndLyQZpaWO5OIf7w9WttVoaHUyUzbpOSP9QUBBGbGHtGEJGVzvKNrd0ipbZdgEVpbpp6gQpzB/RZTYaXUMbMSGadG7VUnBE3/vjftOe6s09595zT7/nfvp8wOWe8znfc7/vaPLMN9/7PeemqpAktesX+h5AkjRZhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxb+t7AIDTTjutZmdn+x5DklaVJ5544vtVNTNsu6kI/ezsLLt27ep7DElaVZI8P8p2nrqRpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMZNxStjpWFmt97b2773b7uwt31LK8EjeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMZ5eaWWpM/LHCUtj0f0ktS4oaFPsj7JQ0meSfJ0kmu69c8lOZhkd/dxwcBzrk+yL8lzST4yyV+AJGlxo5y6OQpcV1VPJjkZeCLJzu6xL1bVXw1unGQjcCnwXuBdwDeT/HpVvb6Sg0uSRjP0iL6qDlXVk93tV4G9wNpFnrIZ+HpVvVZV3wP2AWetxLCSpKVb0jn6JLPAmcCj3dLVSZ5KckuSU7q1tcALA087wAJ/MSTZkmRXkl1HjhxZ8uCSpNGMHPokJwF3ANdW1SvATcB7gE3AIeALS9lxVW2vqrmqmpuZmVnKUyVJSzBS6JOcyHzkb6+qOwGq6qWqer2q3gC+zE9OzxwE1g88fV23JknqwShX3QS4GdhbVTcOrJ8+sNnFwJ7u9j3ApUl+Kcm7gQ3AYys3siRpKUa56uYc4HLg20l2d2ufBi5LsgkoYD/wKYCqejrJDuAZ5q/YucorbiSpP0NDX1XfArLAQ/ct8pwbgBvGmEuStEJ8ZawkNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNW5o6JOsT/JQkmeSPJ3kmm791CQ7k3yn+3xKt54kX0qyL8lTST446V+EJOn4RjmiPwpcV1UbgbOBq5JsBLYCD1bVBuDB7j7AR4EN3ccW4KYVn1qSNLKhoa+qQ1X1ZHf7VWAvsBbYDNzWbXYbcFF3ezPwlZr3CPDOJKev+OSSpJEs6Rx9klngTOBRYE1VHeoeehFY091eC7ww8LQD3ZokqQcjhz7JScAdwLVV9crgY1VVQC1lx0m2JNmVZNeRI0eW8lRJ0hKMFPokJzIf+dur6s5u+aVjp2S6z4e79YPA+oGnr+vWfkpVba+quaqam5mZWe78kqQhRrnqJsDNwN6qunHgoXuAK7rbVwB3D6x/orv65mzgBwOneCRJP2NvG2Gbc4DLgW8n2d2tfRrYBuxIciXwPHBJ99h9wAXAPuBHwCdXdGJJ0pIMDX1VfQvIcR4+b4HtC7hqzLkkSSvEV8ZKUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1ztBLUuMMvSQ1ztBLUuNG+cEj0s+12a339rLf/dsu7GW/ao9H9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0bGvoktyQ5nGTPwNrnkhxMsrv7uGDgseuT7EvyXJKPTGpwSdJoRjmivxU4f4H1L1bVpu7jPoAkG4FLgfd2z/m7JCes1LCSpKUbGvqqehh4ecSvtxn4elW9VlXfA/YBZ40xnyRpTOOco786yVPdqZ1TurW1wAsD2xzo1iRJPVlu6G8C3gNsAg4BX1jqF0iyJcmuJLuOHDmyzDEkScMsK/RV9VJVvV5VbwBf5ienZw4C6wc2XdetLfQ1tlfVXFXNzczMLGcMSdIIlvXDwZOcXlWHursXA8euyLkH+GqSG4F3ARuAx8aeUj+lrx9WLWl1Ghr6JF8DzgVOS3IA+CxwbpJNQAH7gU8BVNXTSXYAzwBHgauq6vXJjC5JGsXQ0FfVZQss37zI9jcAN4wzlCRp5fjKWElqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYZeklqnKGXpMYNDX2SW5IcTrJnYO3UJDuTfKf7fEq3niRfSrIvyVNJPjjJ4SVJw41yRH8rcP6b1rYCD1bVBuDB7j7AR4EN3ccW4KaVGVOStFxDQ19VDwMvv2l5M3Bbd/s24KKB9a/UvEeAdyY5faWGlSQt3XLP0a+pqkPd7ReBNd3ttcALA9sd6NYkST0Z+5uxVVVALfV5SbYk2ZVk15EjR8YdQ5J0HMsN/UvHTsl0nw936weB9QPbrevW3qKqtlfVXFXNzczMLHMMSdIwyw39PcAV3e0rgLsH1j/RXX1zNvCDgVM8kqQevG3YBkm+BpwLnJbkAPBZYBuwI8mVwPPAJd3m9wEXAPuAHwGfnMDMkqQlGBr6qrrsOA+dt8C2BVw17lCSpJXjK2MlqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaN/SHg0vqx+zWe3vZ7/5tF/ayX02OR/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1LixXhmbZD/wKvA6cLSq5pKcCvwzMAvsBy6pqv8Zb0xJ0nKtxBH971TVpqqa6+5vBR6sqg3Ag919SVJPJnHqZjNwW3f7NuCiCexDkjSicd/UrIAHkhTwD1W1HVhTVYe6x18E1oy5j6nV15tOSdJSjBv6366qg0l+FdiZ5NnBB6uqur8E3iLJFmALwBlnnDHmGJKk4xnr1E1VHew+HwbuAs4CXkpyOkD3+fBxnru9quaqam5mZmacMSRJi1h26JO8I8nJx24DHwb2APcAV3SbXQHcPe6QkqTlG+fUzRrgriTHvs5Xq+pfkzwO7EhyJfA8cMn4Y0qSlmvZoa+q7wIfWGD9v4HzxhlKkrRyfGWsJDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS45b9w8EltWl267297Xv/tgt723fLPKKXpMat+iP6Po8+JGk18Ihekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekho3sdAnOT/Jc0n2Jdk6qf1IkhY3kdAnOQH4W+CjwEbgsiQbJ7EvSdLiJvXK2LOAfVX1XYAkXwc2A89MaH+SGtDXK91bf4+dSYV+LfDCwP0DwG9OaF+SNJbW38itt/e6SbIF2NLd/WGS55b5pU4Dvr8yU02Uc64s51xZq2HO1TAjLHHO/OVY+/q1UTaaVOgPAusH7q/r1n6sqrYD28fdUZJdVTU37teZNOdcWc65slbDnKthRpjOOSd11c3jwIYk707yi8ClwD0T2pckaRETOaKvqqNJrgbuB04AbqmqpyexL0nS4iZ2jr6q7gPum9TXHzD26Z+fEedcWc65slbDnKthRpjCOVNVfc8gSZog3wJBkhrXROiT/EWSp5LsTvJAknf1PdNCknw+ybPdrHcleWffMy0kyR8keTrJG0mm6uqB1fLWGkluSXI4yZ6+ZzmeJOuTPJTkme7/9zV9z7SQJG9P8liS/+zm/PO+Z1pMkhOS/EeSb/Q9yzFNhB74fFW9v6o2Ad8A/qzvgY5jJ/C+qno/8F/A9T3Pczx7gN8HHu57kEGr7K01bgXO73uIIY4C11XVRuBs4Kop/e/5GvChqvoAsAk4P8nZPc+0mGuAvX0PMaiJ0FfVKwN33wFM5TcequqBqjra3X2E+dcXTJ2q2ltVy30B2yT9+K01qur/gGNvrTF1quph4OW+51hMVR2qqie7268yH6e1/U71VjXvh93dE7uPqfwznmQdcCHwj33PMqiJ0AMkuSHJC8AfMr1H9IP+BPiXvodYZRZ6a42pC9NqlGQWOBN4tN9JFtadDtkNHAZ2VtVUzgn8NfCnwBt9DzJo1YQ+yTeT7FngYzNAVX2mqtYDtwNXT+uc3TafYf6fzbdP85z6+ZDkJOAO4No3/et4alTV692p2XXAWUne1/dMb5bkY8Dhqnqi71nerLf3ulmqqvrdETe9nfnr9z87wXGOa9icSf4Y+BhwXvV4besS/ntOk6FvraGlSXIi85G/varu7HueYarqf5M8xPz3P6btG93nAB9PcgHwduBXkvxTVf1Rz3OtniP6xSTZMHB3M/BsX7MsJsn5zP+z7uNV9aO+51mFfGuNFZQkwM3A3qq6se95jifJzLEr1JL8MvB7TOGf8aq6vqrWVdUs8783/20aIg+NhB7Y1p12eAr4MPPf9Z5GfwOcDOzsLgX9+74HWkiSi5McAH4LuDfJ/X3PBPNvrcH8abn7mf/G4Y5pfWuNJF8D/h34jSQHklzZ90wLOAe4HPhQ9/txd3c0Om1OBx7q/nw/zvw5+qm5dHE18JWxktS4Vo7oJUnHYeglqXGGXpIaZ+glqXGGXpIaZ+glqXGGXpIaZ+glqXH/D1XXli1UlmCdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#\n",
    "# rand() returns a single random number:\n",
    "#\n",
    "\n",
    "print(np.random.rand())\n",
    "\n",
    "#\n",
    "# hist plots a histogram of an array of numbers\n",
    "#\n",
    "\n",
    "print(pl.hist(np.random.normal(size=1000)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 97., 102.,  99., 100.,  84.,  97., 111.,  94., 119.,  97.]),\n",
       " array([0.00389839, 0.20333071, 0.40276302, 0.60219534, 0.80162766,\n",
       "        1.00105997, 1.20049229, 1.39992461, 1.59935692, 1.79878924,\n",
       "        1.99822156]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEG9JREFUeJzt3X2sZHV9x/H3p6xAQSvg3lLkwYW4qQGjgd5Qn6IoJiI+LE0Ngdhm0W22tmi1NlUoSWmamGLa+NC0tdkAdU0oD0Ut1GIrBYxpza69IPIosC4guwH2Kg9KTVDst3/MWR3X3b33zpmZu/z6fiWTOed3fmfOd39z9nPPPWfm3FQVkqR2/cJyFyBJmiyDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktS4FctdAMDKlStr1apVy12GJD2r3Hzzzd+pqpmF+u0TQb9q1Srm5uaWuwxJelZJ8uBi+nnqRpIaZ9BLUuMMeklq3IJBn+TSJDuS3DHU9pdJvpnktiSfT3LI0LLzk2xJck+SN02qcEnS4izmiP7TwGm7tF0PvLSqXgbcC5wPkOR44CzghG6dv0uy39iqlSQt2YJBX1VfAR7bpe1LVfVMN7sJOKqbXgNcUVVPV9X9wBbg5DHWK0laonGco3838MVu+kjgoaFl27o2SdIy6RX0SS4AngEuG2Hd9UnmkszNz8/3KUOStBcjB32Sc4C3Au+sn/7h2e3A0UPdjurafk5Vbaiq2aqanZlZ8ItdkqQRjfTN2CSnAR8CXldVPxhadC3wj0k+BrwQWA18rXeVkv5fWHXevy7Ldh+46C3Lst1pWTDok1wOnAKsTLINuJDBp2wOAK5PArCpqt5TVXcmuQq4i8EpnXOr6seTKl6StLAFg76qzt5N8yV76f8R4CN9ipIkjY/fjJWkxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekho30t0rJU2ed3LUuHhEL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIat+BtipNcCrwV2FFVL+3aDgOuBFYBDwBnVtXjSQJ8Ejgd+AFwTlXdMpnSl5+3kZX0bLCYI/pPA6ft0nYecENVrQZu6OYB3gys7h7rgU+Np0xJ0qgWDPqq+grw2C7Na4CN3fRG4Iyh9s/UwCbgkCRHjKtYSdLSjXqO/vCqeribfgQ4vJs+EnhoqN+2rk2StEx6/ynBqqoktdT1kqxncHqHY445pm8ZmhKvS0jPPqMe0T+685RM97yja98OHD3U76iu7edU1Yaqmq2q2ZmZmRHLkCQtZNQj+muBtcBF3fM1Q+3vTXIF8OvAk0OneDQmy3VULenZaTEfr7wcOAVYmWQbcCGDgL8qyTrgQeDMrvt1DD5auYXBxyvfNYGaJUlLsGDQV9XZe1h06m76FnBu36IkSePT+2LscvM0hibNfax9y/keT+ODBt4CQZIaZ9BLUuMMeklqnEEvSY171l+MlTReXnxuj0f0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1Dj/wpSeFfyrR9LoPKKXpMYZ9JLUuF5Bn+QPk9yZ5I4klyc5MMmxSTYn2ZLkyiT7j6tYSdLSjRz0SY4E/gCYraqXAvsBZwEfBT5eVS8GHgfWjaNQSdJo+p66WQH8YpIVwEHAw8AbgKu75RuBM3puQ5LUw8hBX1Xbgb8Cvs0g4J8EbgaeqKpnum7bgCP7FilJGl2fUzeHAmuAY4EXAgcDpy1h/fVJ5pLMzc/Pj1qGJGkBfU7dvBG4v6rmq+pHwOeAVwOHdKdyAI4Ctu9u5araUFWzVTU7MzPTowxJ0t70CfpvA69IclCSAKcCdwE3Ae/o+qwFrulXoiSpjz7n6DczuOh6C3B791obgA8DH0yyBXgBcMkY6pQkjajXLRCq6kLgwl2atwIn93ldSdL4+M1YSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjesV9EkOSXJ1km8muTvJK5McluT6JPd1z4eOq1hJ0tL1PaL/JPBvVfUS4OXA3cB5wA1VtRq4oZuXJC2TkYM+yfOB1wKXAFTVD6vqCWANsLHrthE4o2+RkqTR9TmiPxaYB/4hydeTXJzkYODwqnq46/MIcPjuVk6yPslckrn5+fkeZUiS9qZP0K8ATgI+VVUnAv/DLqdpqqqA2t3KVbWhqmaranZmZqZHGZKkvekT9NuAbVW1uZu/mkHwP5rkCIDueUe/EiVJfYwc9FX1CPBQkl/tmk4F7gKuBdZ2bWuBa3pVKEnqZUXP9d8HXJZkf2Ar8C4GPzyuSrIOeBA4s+c2JEk99Ar6qroVmN3NolP7vK4kaXz8ZqwkNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TG9Q76JPsl+XqSL3TzxybZnGRLkiuT7N+/TEnSqMZxRP9+4O6h+Y8CH6+qFwOPA+vGsA1J0oh6BX2So4C3ABd38wHeAFzdddkInNFnG5Kkfvoe0X8C+BDwv938C4AnquqZbn4bcOTuVkyyPslckrn5+fmeZUiS9mTkoE/yVmBHVd08yvpVtaGqZqtqdmZmZtQyJEkLWNFj3VcDb09yOnAg8EvAJ4FDkqzojuqPArb3L1OSNKqRj+ir6vyqOqqqVgFnATdW1TuBm4B3dN3WAtf0rlKSNLJJfI7+w8AHk2xhcM7+kglsQ5K0SH1O3fxEVX0Z+HI3vRU4eRyvK0nqz2/GSlLjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjRg76JEcnuSnJXUnuTPL+rv2wJNcnua97PnR85UqSlqrPEf0zwB9V1fHAK4BzkxwPnAfcUFWrgRu6eUnSMhk56Kvq4aq6pZv+PnA3cCSwBtjYddsInNG3SEnS6MZyjj7JKuBEYDNweFU93C16BDh8HNuQJI2md9AneS7wWeADVfW94WVVVUDtYb31SeaSzM3Pz/ctQ5K0B72CPslzGIT8ZVX1ua750SRHdMuPAHbsbt2q2lBVs1U1OzMz06cMSdJe9PnUTYBLgLur6mNDi64F1nbTa4FrRi9PktTXih7rvhr4beD2JLd2bX8CXARclWQd8CBwZr8SJUl9jBz0VfWfQPaw+NRRX1eSNF5+M1aSGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjJhb0SU5Lck+SLUnOm9R2JEl7N5GgT7If8LfAm4HjgbOTHD+JbUmS9m5SR/QnA1uqamtV/RC4AlgzoW1JkvZiUkF/JPDQ0Py2rk2SNGUrlmvDSdYD67vZp5LcM+JLrQS+M56qxsq6lsa6lm5frc26liAf7VXXixbTaVJBvx04emj+qK7tJ6pqA7Ch74aSzFXVbN/XGTfrWhrrWrp9tTbrWppp1DWpUzf/DaxOcmyS/YGzgGsntC1J0l5M5Ii+qp5J8l7g34H9gEur6s5JbEuStHcTO0dfVdcB103q9Yf0Pv0zIda1NNa1dPtqbda1NBOvK1U16W1IkpaRt0CQpMbts0G/0C0UkhyQ5Mpu+eYkq4aWnd+135PkTVOu64NJ7kpyW5IbkrxoaNmPk9zaPcZ+cXoRtZ2TZH6oht8ZWrY2yX3dY+2U6/r4UE33JnliaNnExizJpUl2JLljD8uT5K+7um9LctLQskmO10J1vbOr5/YkX03y8qFlD3TttyaZm3JdpyR5cuj9+tOhZRO7Jcoi6vrjoZru6Papw7plkxyvo5Pc1OXBnUnev5s+09nHqmqfezC4gPst4Dhgf+AbwPG79Pl94O+76bOAK7vp47v+BwDHdq+z3xTrej1wUDf9ezvr6uafWuYxOwf4m92sexiwtXs+tJs+dFp17dL/fQwu3k9jzF4LnATcsYflpwNfBAK8Atg86fFaZF2v2rk9BrcZ2Ty07AFg5TKN1ynAF/ruA+Oua5e+bwNunNJ4HQGc1E0/D7h3N/8np7KP7atH9Iu5hcIaYGM3fTVwapJ07VdU1dNVdT+wpXu9qdRVVTdV1Q+62U0MvkMwDX1uO/Em4PqqeqyqHgeuB05bprrOBi4f07b3qqq+Ajy2ly5rgM/UwCbgkCRHMNnxWrCuqvpqt12Y4j62iPHak4neEmWJdU1z/3q4qm7ppr8P3M3P3yFgKvvYvhr0i7mFwk/6VNUzwJPACxa57iTrGraOwU/rnQ5MMpdkU5IzxlTTUmv7ze5XxKuT7PxS2z4xZt1prmOBG4eaJzlmC9lT7fvSLT523ccK+FKSmzP49vm0vTLJN5J8MckJXds+MV5JDmIQlp8dap7KeGVwavlEYPMui6ayjy3bLRBal+S3gFngdUPNL6qq7UmOA25McntVfWuKZf0LcHlVPZ3kdxn8RvSGKW5/IWcBV1fVj4falnvM9llJXs8g6F8z1Pyabrx+Gbg+yTe7I95puIXB+/VUktOBfwZWT2nbi/E24L+qavjof+LjleS5DH64fKCqvjfO116sffWIfsFbKAz3SbICeD7w3UWuO8m6SPJG4ALg7VX19M72qtrePW8FvszgJ/y4LOa2E98dqudi4NcWu+4k6xpyFrv8Wj3hMVvInmqf5HgtSpKXMXgP11TVd3e2D43XDuDzjO+05YKq6ntV9VQ3fR3wnCQr2QfGq7O3/Wsi45XkOQxC/rKq+txuukxnH5vERYi+Dwa/aWxl8Gv8zos3J+zS51x+9mLsVd30CfzsxditjO9i7GLqOpHBhafVu7QfChzQTa8E7mO8F6QWU9sRQ9O/AWyqn174ub+r8dBu+rBp1dX1ewmDC2OZ1ph1r7uKPV9cfAs/e6Hsa5Mer0XWdQyDa0+v2qX9YOB5Q9NfBU6bYl2/svP9YxCY3+7GblH7wKTq6pY/n8F5/IOnNV7dv/0zwCf20mcq+9jYBnrcDwZXo+9lEJoXdG1/zuAoGeBA4J+6Hf5rwHFD617QrXcP8OYp1/UfwKPArd3j2q79VcDt3U5+O7BuGcbsL4A7uxpuAl4ytO67u7HcArxrmnV1838GXLTLehMdMwZHdw8DP2JwDnQd8B7gPd3yMPgDOt/qtj87pfFaqK6LgceH9rG5rv24bqy+0b3PF0y5rvcO7V+bGPpBtLt9YFp1dX3OYfAhjeH1Jj1er2FwDeC2offq9OXYx/xmrCQ1bl89Ry9JGhODXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxv0fg4jvZO90QOQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m=28*1.67e-27  # mass of a molecule (e.g., Nitrogen)\n",
    "g=9.8          # grav field strength\n",
    "kb=1.67e-23    # boltzman constant\n",
    "demonE = 0.0   # initial demon energy\n",
    "N=10000        # number of molecules\n",
    "M=400000       # number of iterations\n",
    "h=20000.0      # height scale\n",
    "\n",
    "def setup(N=100,L=1.0):\n",
    "    y=L*np.random.rand(N)     # put N particles at random heights (y) between 0 and L\n",
    "    return y\n",
    "\n",
    "yarray = setup(N=1000,L=2.0)\n",
    "pl.hist(yarray)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "76.894 percent accepted\n",
      "Averge height=9966.415m\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYHVW57/HvjwBhnkzEkBA6SMADihFaHAAFURmFo94j4aAMopEDKCheCeDRyDmeGxkcuCoYBCEIYVYigzIcFDwK2IFABomEECQxJAEEwnAjSd77R60dqju7d1d17917d/L7PE89XXtV1ao3vXf63WutqlWKCMzMzMpYr9kBmJnZwOPkYWZmpTl5mJlZaU4eZmZWmpOHmZmV5uRhZmalOXlYS5F0saR/r1NdIyW9LGlQev1bSZ+rR92pvtslHVuv+kqc9z8lPSvpmSrb9pO0oJf17itpTsF9e30eWzus3+wAbN0haT6wLbACWAnMBiYDkyJiFUBEnFiirs9FxF3d7RMRfwU261vUq883AdgpIj6dq//getRdMo6RwOnADhGxpJ51R8R9wC71qEvS5cCCiPh6Peqz1uOWh/W3j0XE5sAOwETgDODSep9E0tr6xWgk8Fy9E4dZWU4e1hQR8WJETAWOBI6V9HbIvrFK+s+0PkTSLZJekPS8pPskrSfpSrI/or9K3VJfk9QmKSSdIOmvwH/nyvKJ5K2SHpT0kqSbJW2TzrVGN4yk+ZI+LOkg4CzgyHS+R9L21d1gKa6vS3pK0hJJkyVtmbZV4jhW0l9Tl9PZ3f1uJG2Zjl+a6vt6qv/DwJ3AdimOy2vUcXqKY5Gk43PlgyWdn+JYnLoJN672O5C0h6SHJS2TdL2kayvvTa3zSBoHHA18LcX5q1R+hqSFqb45kg7oLn5rfU4e1lQR8SCwANi3yubT07ahZN1dZ2WHxGeAv5K1YjaLiHNzx3wQ+CfgwG5OeQzwWWAYWffZhQVi/DXwX8C16XzvrLLbcWnZH9iRrLvsh1322YesW+gA4BuS/qmbU/5fYMtUzwdTzMenLrqDgb+lOI7r5vi3pOOHAycAP5K0ddo2EdgZGAPslPb5RtcKJG0I/AK4HNgGmAJ8vMh5ImIScBVwborzY5J2AU4B3p1angcC87uJ3wYAJw9rBX8j+wPV1etkf+R3iIjXI+K+6HkytgkR8UpEvNbN9isjYmZEvAL8O/CpyoB6Hx0NfDci5kXEy8CZwNgurZ5vRcRrEfEI8AiwRhJKsYwFzoyIZRExH7gA+EyJWF4Hzkm/s9uAl4FdJAkYB3w5Ip6PiGVkSXFslTreSzYmemGq5ybgwSLn6SamlcBgYFdJG0TE/Ih4osS/yVqMk4e1guHA81XKzwPmAndImidpfIG6ni6x/SlgA2BIoShr2y7Vl697fbIWU0X+6qhXqT6YPyTF1LWu4SVieS4iVlQ511BgE2Ba6gp8Afh1Ku9qO2Bhl2Td9Xfb3XnWEBFzgdOACcASSddI2q7Ev8lajJOHNZWkd5P9Yfx9123pm/fpEbEjcDjwlVw/eXctkJ5aJtvn1keSfXt+FniF7A9rJa5BdP6j2lO9fyO7CCBf9wpgcQ/HdfVsiqlrXQtL1tNd3a8Bu0XEVmnZMiKq/cFfBAxPrZWK7avs1501fl8RcXVE7EP2bwvgOyXqsxbj5GFNIWkLSYcB1wA/j4gZVfY5TNJO6Q/Yi2RdH6vS5sVkYwJlfVrSrpI2Ac4BboiIlcBfgI0kHSppA+DrZN0sFYuBNknd/Z+ZAnxZ0ihJm/HGGMmKbvavKsVyHfBtSZtL2gH4CvDzMvV0U/cq4BLge5LeDCBpuKRq40N/JPt9nyJpfUlHAHuVOF2n90fSLpI+JGkw8P/Iktiq7g621ufkYf3tV5KWkXWBnA18Fzi+m31HA3eR9aX/EfhxRNyTtv0f4Oup++WrJc5/Jdkg8DPARsCXILv6CzgJ+CnZt/xXyAbrK65PP5+T9FCVei9Ldd8LPEn2B/KLJeLK+2I6/zyyFtnVqf56OIOsK/B+SS+R/X7XGKeIiH8AnyAbCH8B+DRwC7C84HkuJRvfeEHSL8kS8USy1s8zwJvJxoVsgJIfBmVmRUh6ALg4In7W7Fis+dzyMLOqJH1Q0ltSt9WxwO5kA+xmnp7EzLq1C9n4y6ZkXWj/KyIWNTckaxXutjIzs9LcbWVmZqWttd1WQ4YMiba2tmaHYWY2YEybNu3ZiKh20+ga1trk0dbWRkdHR7PDMDMbMCQ91fNeGXdbmZlZaU4eZmZWmpOHmZmV5uRhZmalOXmYmVlpTh5mZlaak4eZmZXm5GFmZqU5eZiZWWlr7R3mfdE2/tZut82feGg/RmJm1poa1vKQdJmkJZJm5squlTQ9LfMlTU/lbZJey227OHfMnpJmSJor6cIuz1Q2M7MmaGTL43Lgh8DkSkFEHFlZl3QB2XOpK56IiDFV6rkI+DzwAHAbcBBwewPiNTOzghrW8oiIe4Hnq21LrYdPAVNq1SFpGLBFRNwf2YNHJgP/XO9YzcysnGYNmO8LLI6Ix3NloyQ9LOl3kvZNZcOBBbl9FqSyqiSNk9QhqWPp0qX1j9rMzIDmJY+j6NzqWASMjIh3AV8Brpa0RdlKI2JSRLRHRPvQoYWmpDczs17o96utJK0PfALYs1IWEcuB5Wl9mqQngJ2BhcCI3OEjUpmZmTVRM1oeHwYei4jV3VGShkoalNZ3BEYD8yJiEfCSpPemcZJjgJubELOZmeU08lLdKcAfgV0kLZB0Qto0ljUHyj8APJou3b0BODEiKoPtJwE/BeYCT+ArrczMmq5h3VYRcVQ35cdVKbsRuLGb/TuAt9c1ODMz6xNPT2JmZqU5eZiZWWlOHmZmVpqTh5mZlebkYWZmpTl5mJlZaU4eZmZWmpOHmZmV5uRhZmalOXmYmVlpTh5mZlaak4eZmZXm5GFmZqU5eZiZWWlOHmZmVpqTh5mZlebkYWZmpTl5mJlZaU4eZmZWmpOHmZmV1rDkIekySUskzcyVTZC0UNL0tByS23ampLmS5kg6MFd+UCqbK2l8o+I1M7PiGtnyuBw4qEr59yJiTFpuA5C0KzAW2C0d82NJgyQNAn4EHAzsChyV9jUzsyZav1EVR8S9ktoK7n4EcE1ELAeelDQX2CttmxsR8wAkXZP2nV3ncM3MrIRmjHmcIunR1K21dSobDjyd22dBKuuuvCpJ4yR1SOpYunRpveM2M7Okv5PHRcBbgTHAIuCCelYeEZMioj0i2ocOHVrPqs3MLKdh3VbVRMTiyrqkS4Bb0suFwPa5XUekMmqUm5lZk/Rry0PSsNzLjwOVK7GmAmMlDZY0ChgNPAj8CRgtaZSkDckG1af2Z8xmZramhrU8JE0B9gOGSFoAfBPYT9IYIID5wBcAImKWpOvIBsJXACdHxMpUzynAb4BBwGURMatRMZuZWTGKiGbH0BDt7e3R0dHRq2Pbxt/a6/POn3hor481M2smSdMior3Ivr7D3MzMSiucPCRt0shAzMxs4OgxeUh6v6TZwGPp9Tsl/bjhkZmZWcsq0vL4HnAg8BxARDwCfKCRQZmZWWsr1G0VEU93KVrZgFjMzGyAKHKp7tOS3g+EpA2AU4E/NzYsMzNrZUVaHicCJ5PNKbWQbGqRkxsZlJmZtbYeWx4R8SxwdD/EYmZmA0SPyUPShVWKXwQ6IuLm+odkZmatrki31UZkXVWPp2V3sgkKT5D0/QbGZmZmLarIgPnuwN65uaYuAu4D9gFmNDA2MzNrUUVaHlsDm+Vebwpsk5LJ8oZEZWZmLa1Iy+NcYLqk3wIiu0HwvyRtCtzVwNjMzKxFFbna6lJJt/HGM8XPioi/pfX/3bDIzMysZRWdGHE9YCnwd2AnSZ6exMxsHVbkUt3vAEcCs4BVqTiAexsYl5mZtbAiYx7/DOwSER4cNzMzoFi31Txgg0YHYmZmA0eRlserZFdb3U3u0tyI+FLDojIzs5ZWJHlMTYuZmRlQ7FLdK3pTsaTLgMOAJRHx9lR2HvAx4B/AE8DxEfGCpDayad7npMPvj4gT0zF7ApcDGwO3AadGRPQmJjMzq48ij6EdLekGSbMlzassBeq+HDioS9mdwNsjYnfgL8CZuW1PRMSYtJyYK78I+DwwOi1d6zQzs35WZMD8Z2R/wFcA+wOTgZ/3dFBE3As836XsjohYkV7eTzbBYrckDQO2iIj7U2tjMtnVX2Zm1kRFksfGEXE3oIh4KiImAIfW4dyfBW7PvR4l6WFJv5O0byobDizI7bMglVUlaZykDkkdS5curUOIZmZWTZEB8+WS1gMel3QK2dMEN+vhmJoknU3WkrkqFS0CRkbEc2mM45eSditbb0RMAiYBtLe3e1zEzKxBirQ8TgU2Ab4E7Al8Gji2tyeUdBzZQPrRlYHviFgeEc+l9Wlkg+k7kyWqfNfWiFRmZmZNVDN5SBoEHBkRL0fEgog4PiI+GRH39+Zkkg4CvgYcHhGv5sqHpnMhaUeygfF5EbEIeEnSeyUJOAbw0wvNzJqsZrdVRKyUtE9vKpY0BdgPGCJpAfBNsqurBgN3Zrlg9SW5HwDOkfQ62fxZJ0ZEZbD9JN64VPd2Oo+TmJlZExQZ83hY0lTgeuCVSmFE3FTroIg4qkrxpd3seyNwYzfbOoC3F4jTzMz6SZHksRHwHPChXFkANZOHmZmtvYrcYX58fwSytmgbf2vN7fMn1uMqZzOz5ir6MCgzM7PVnDzMzKy0InNbjSpSZmZm644iLY9qV0HdUO9AzMxs4Oh2wFzS24DdgC0lfSK3aQuyK7DMzGwdVetqq13IphHZiuwZHBXLyKZINzOzdVS3ySMibgZulvS+iPhjP8ZkZmYtrshNgnMlnQW05fePiM82KigzM2ttRZLHzcB9wF3AysaGY2ZmA0GR5LFJRJzR8EjMzGzAKHKp7i2SDml4JGZmNmDUulR3GdkEiALOkrQceD29jojYon9CNDOzVlPraqvN+zMQMzMbOHoc85C0R5XiF4GnImJF/UMyM7NWV2TA/MfAHsCM9PodwEyyO8//LSLuaFRwZmbWmooMmP8NeFdE7BkRewJjgHnAR4BzGxmcmZm1piLJY+eImFV5ERGzgbdFxLzGhWVmZq2sSLfVLEkXAdek10cCsyUNJrv6yszM1jFFWh7HAXOB09IyL5W9Duxf60BJl0laImlmrmwbSXdKejz93DqVS9KFkuZKejQ/UC/p2LT/45KOLfuPNDOz+uoxeUTEaxFxQUR8PC3nR8SrEbEqIl7u4fDLgYO6lI0H7o6I0cDd6TXAwcDotIwDLoIs2QDfBN4D7AV8s5JwzMysOWrdJHhdRHxK0gyymwU7iYjde6o8Iu6V1Nal+Ahgv7R+BfBb4IxUPjkiArhf0laShqV974yI51Ncd5IlpCk9nb8VtY2/tdtt8yce2o+RmJn1Xq0xj1PTz8PqfM5tI2JRWn8G2DatDweezu23IJV1V74GSePIWi2MHDmyjiGbmVlet91WlT/wEfFUKhqd1pcAz9fj5KmVsUarpg/1TYqI9ohoHzp0aL2qNTOzLnoc85D0ebJnlv8kFY0AftmHcy5O3VGkn0tS+UJg+9x+I1JZd+VmZtYkRa62OhnYG3gJICIeB97ch3NOBSpXTB1L9ryQSvkx6aqr9wIvptbPb4CPSto6DZR/NJWZmVmTFLnPY3lE/EMSAJLWp2BXk6QpZAPeQyQtILtqaiJwnaQTgKeAT6XdbwMOIbss+FXgeICIeF7SfwB/SvudUxk8NzOz5iiSPH6XHkO7saSPACcBvypSeUQc1c2mA6rsG2StnGr1XAZcVuScZmbWeEW6rcYDS8kmRvwCWQvh640MyszMWluPLY+IWAVckhYzM7NCz/PYG5gA7JD2rzxJcMfGhmZmZq2qyJjHpcCXgWnAysaGY2ZmA0GR5PFiRNze8EjMzGzAqDW3VWVW23sknQfcBCyvbI+Ihxocm5mZtahaLY8Lurxuz60H8KH6h2NmZgNBt8kjImo+q8PMzNZdRe7zMDMz68TJw8zMSnPyMDOz0opcqouk9wNt+f0jYnKDYjIzsxZX5A7zK4G3AtN54ybBAJw8zMzWUUVaHu3ArmnWWzMzs0JjHjOBtzQ6EDMzGziKtDyGALMlPUjnO8wPb1hUZmbW0ookjwmNDsIybeNvrbl9/sRD+ykSM7PaijzP43f9EYiZmQ0ctSZG/H1E7CNpGZ2fWV55nscWDY/OzMxaUq25rfZJPzfvv3DMzGwg6Pc7zCXtIml6bnlJ0mmSJkhamCs/JHfMmZLmSpoj6cD+jtnMzDordId5PUXEHGAMgKRBwELgF8DxwPci4vz8/pJ2BcYCuwHbAXdJ2jki/FRDM7MmafbcVgcAT0TEUzX2OQK4JiKWR8STwFxgr36JzszMqqqZPCQNknRPA88/FpiSe32KpEclXSZp61Q2HHg6t8+CVLYGSeMkdUjqWLp0aWMiNjOz2skjdQ2tkrRlvU8saUPgcOD6VHQR2RxaY4BFrPkkwx5FxKSIaI+I9qFDh9YtVjMz66zImMfLwAxJdwKvVAoj4kt9PPfBwEMRsTjVt7iyQdIlwC3p5UJg+9xxI1KZmZk1SZHkcVNa6u0ocl1WkoZFxKL08uNkc2oBTAWulvRdsgHz0cCDDYjHzMwKKnKH+RWSNgZGpiul+kzSpsBHgC/kis+VNIbshsT5lW0RMUvSdcBsYAVwsq+0MjNrriLP8/gYcD6wITAq/YE/py8TI0bEK8CbupR9psb+3wa+3dvzmZlZfRW5VHcC2aWxLwBExHRgxwbGZGZmLa5I8ng9Il7sUraqEcGYmdnAUGTAfJakfwUGSRoNfAn4Q2PDMjOzVlYkeXwROJvsQVBTgN8A/9HIoKy6Ws/78LM+zKw/Fbna6lXgbEnfyV7GssaHZWZmrazHMQ9J75Y0A3iU7GbBRyTt2fjQzMysVRXptroUOCki7gOQtA/wM2D3RgZmZmatq8jVVisriQMgIn5PdrOemZmto2o9hnaPtPo7ST8hGywP4Ejgt40PzczMWlWtbquus9p+M7cemJnZOqvWM8z3789AzMxs4Cgyt9VWwDFAW37/OkzJbmZmA1SRq61uA+4HZuBpSVpWrRsIwTcRmll9FUkeG0XEVxoeiZmZDRhFLtW9UtLnJQ2TtE1laXhkZmbWsoq0PP4BnEc2v1XlKqvA07Kbma2ziiSP04GdIuLZRgdjZmYDQ5Fuq7nAq40OxMzMBo4iLY9XgOmS7iGblh3wpbpmZuuyIsnjl2kxMzMDij3P44pGnFjSfGAZsBJYERHt6Squa8luSJwPfCoi/i5JwA+AQ8i60I6LiIcaEZeZmfWsyPM8npQ0r+tSp/PvHxFjIqI9vR4P3B0Ro4G702uAg4HRaRkHXFSn85uZWS8U6bZqz61vBPwL0Kj7PI4A9kvrV5DN3ntGKp8cEQHcL2krScMiYlGD4jAzsxp6bHlExHO5ZWFEfB+ox1wXAdwhaZqkcals21xCeAbYNq0PB57OHbsglXUiaZykDkkdS5curUOIZmZWTZGJEffIvVyPrCVSpMXSk30iYqGkNwN3SnosvzEiQlKpqd8jYhIwCaC9vd3TxpuZNUiRJJB/rscK0kB2X08cEQvTzyWSfgHsBSyudEdJGgYsSbsvBLbPHT4ilZmZWRMUudqq7s/1kLQpsF5ELEvrHwXOAaYCxwIT08+b0yFTgVMkXQO8B3jR4x1mZs1TpNtqMPBJ1nyexzl9OO+2wC+yK3BZH7g6In4t6U/AdZJOAJ7ijRbObWSX6Vbudj++D+c2M7M+KtJtdTPwIjCN3B3mfRER84B3Vil/DjigSnkAJ9fj3GZm1ndFkseIiDio4ZFYQ/lhUWZWT0UmRvyDpHc0PBIzMxswirQ89gGOk/QkWbeVyHqSdm9oZGZm1rKKJI+DGx6FmZkNKEUu1X2qPwIxM7OBo8iYh5mZWSdOHmZmVpqTh5mZlebkYWZmpTl5mJlZaU4eZmZWmpOHmZmVVo+HOtlaoNbcV573ysy6csvDzMxKc/IwM7PSnDzMzKw0Jw8zMyvNA+bWIz9Iysy6csvDzMxKc8vD+syX+Zqte/q95SFpe0n3SJotaZakU1P5BEkLJU1PyyG5Y86UNFfSHEkH9nfMZmbWWTNaHiuA0yPiIUmbA9Mk3Zm2fS8izs/vLGlXYCywG7AdcJeknSNiZb9GbWZmq/V7yyMiFkXEQ2l9GfBnYHiNQ44AromI5RHxJDAX2KvxkZqZWXeaOuYhqQ14F/AAsDdwiqRjgA6y1snfyRLL/bnDFtBNspE0DhgHMHLkyIbFbcX5Si2ztVPTrraStBlwI3BaRLwEXAS8FRgDLAIuKFtnREyKiPaIaB86dGhd4zUzszc0JXlI2oAscVwVETcBRMTiiFgZEauAS3ija2ohsH3u8BGpzMzMmqQZV1sJuBT4c0R8N1c+LLfbx4GZaX0qMFbSYEmjgNHAg/0Vr5mZrakZYx57A58BZkiansrOAo6SNAYIYD7wBYCImCXpOmA22ZVaJ/tKKzOz5ur35BERvwdUZdNtNY75NvDthgVlZmaleHoSMzMrzcnDzMxK89xW1lS+D8RsYHLLw8zMSnPyMDOz0pw8zMysNCcPMzMrzcnDzMxK89VW1tJ6uhqrFl+pZdY4bnmYmVlpTh5mZlaau61sreUbEM0axy0PMzMrzS0PW2fVapm4VWJWm1seZmZWmpOHmZmV5m4rsyo82G5Wm1seZmZWmlseZr3gwXZb1zl5mPUzd4nZ2mDAJA9JBwE/AAYBP42IiU0OyayqvszH1de6nXisvwyI5CFpEPAj4CPAAuBPkqZGxOzmRmZWf41MPn05rxOT5Q2I5AHsBcyNiHkAkq4BjgCcPMxymtXq6SmxeIxo7TNQksdw4Onc6wXAe7ruJGkcMC69fFnSnF6ebwjwbC+PbSTHVU4rxtWKMUEf49J3en/iHo5dK39fDdTXuHYouuNASR6FRMQkYFJf65HUERHtdQiprhxXOa0YVyvGBI6rLMc1cO7zWAhsn3s9IpWZmVkTDJTk8SdgtKRRkjYExgJTmxyTmdk6a0B0W0XECkmnAL8hu1T3soiY1cBT9rnrq0EcVzmtGFcrxgSOq6x1Pi5FRH+dy8zM1hIDpdvKzMxaiJOHmZmV5uSRI+kgSXMkzZU0vkHnuEzSEkkzc2XbSLpT0uPp59apXJIuTPE8KmmP3DHHpv0fl3RsrnxPSTPSMRdKUsG4tpd0j6TZkmZJOrUVYpO0kaQHJT2S4vpWKh8l6YFU17XpQgokDU6v56btbbm6zkzlcyQdmCvv1fsuaZCkhyXd0kIxzU+/4+mSOlJZK3y+tpJ0g6THJP1Z0vuaHZekXdLvqbK8JOm0ZseVjvty+rzPlDRF2f+Dpn++OokIL9m4zyDgCWBHYEPgEWDXBpznA8AewMxc2bnA+LQ+HvhOWj8EuB0Q8F7ggVS+DTAv/dw6rW+dtj2Y9lU69uCCcQ0D9kjrmwN/AXZtdmxp383S+gbAA6mO64Cxqfxi4N/S+knAxWl9LHBtWt81vaeDgVHpvR7Ul/cd+ApwNXBLet0KMc0HhnQpa4XP1xXA59L6hsBWrRBXl///z5DdJNfsz/xw4Elg49zn6rhW+Hx1irPsAWvrArwP+E3u9ZnAmQ06Vxudk8ccYFhaHwbMSes/AY7quh9wFPCTXPlPUtkw4LFceaf9SsZ4M9lcYi0TG7AJ8BDZ7ALPAut3fe/Irsh7X1pfP+2nru9nZb/evu9k9xrdDXwIuCWdo6kxpX3ns2byaOp7CGxJ9sdQrRRXl1g+CvxPK8TFGzNqbJM+L7cAB7bC5yu/uNvqDdWmQBneT+feNiIWpfVngG17iKlW+YIq5aWkZu+7yL7lNz02Zd1D04ElwJ1k35peiIgVVepaff60/UXgTb2ItyffB74GrEqv39QCMQEEcIekacqm64Hmv4ejgKXAz5R18/1U0qYtEFfeWGBKWm9qXBGxEDgf+CuwiOzzMo3W+Hyt5uTRYiL7KtC066clbQbcCJwWES/ltzUrtohYGRFjyL7t7wW8rb9jyJN0GLAkIqY1M45u7BMRewAHAydL+kB+Y5Pew/XJumovioh3Aa+QdQc1Oy4A0tjB4cD1Xbc1I640xnIEWdLdDtgUOKg/YyjCyeMNzZwCZbGkYQDp55IeYqpVPqJKeSGSNiBLHFdFxE2tFBtARLwA3EPW7N5KUuUm13xdq8+ftm8JPNeLeGvZGzhc0nzgGrKuqx80OSZg9bdWImIJ8AuyZNvs93ABsCAiHkivbyBLJs2Oq+Jg4KGIWJxeNzuuDwNPRsTSiHgduInsM9f0z1cnZfu51taF7NvRPLJsXxlE2q1B52qj85jHeXQeoDs3rR9K5wG6B1P5NmR9yFun5Ulgm7St6wDdIQVjEjAZ+H6X8qbGBgwFtkrrGwP3AYeRfUvMDx6elNZPpvPg4XVpfTc6Dx7OIxs47NP7DuzHGwPmTY2J7Bvq5rn1P5B9Y22Fz9d9wC5pfUKKqelxpWOvAY5voc/8e4BZZGN8IrvY4IvN/nytEWfZA9bmhexqir+Q9amf3aBzTCHrx3yd7BvZCWT9k3cDjwN35T54InsI1hPADKA9V89ngblpyX/w24GZ6Zgf0mWQskZc+5A1zx8FpqflkGbHBuwOPJzimgl8I5XvmP5jzk3/qQan8o3S67lp+465us5O555D7qqXvrzvdE4eTY0pnf+RtMyqHNfs9zAdNwboSO/jL8n+yLZCXJuSfUvfMlfWCnF9C3gsHXslWQJoic98ZfH0JGZmVprHPMzMrDQnDzMzK83Jw8zMSnPyMDOz0pw8zMysNCcPs0RSm3KzHRc85kRJx/Swz3GSftjNtrNqHCdJ/y1pixLxHCbpnKL7m/WWk4dZH0TExRExuQ9VdJs8yK7FfyS6TBPTg1uBj0napA8xmfXIycOss0GSLknPUrhD0sYAkt4q6ddpwsH7JL0tlU+Q9NW0/u70nIfpks7r0orZLh3/uKRz0/4TgY3T/ldVieVostmNK62ixyRdLukvkq6S9GFJ/5Pq3AtWz8X0W7K78M0axsnDrLPRwI9fVZ+GAAABqUlEQVQiYjfgBeCTqXwS8MWI2BP4KvDjKsf+DPhCZJM4ruyybQxwJPAO4EhJ20fEeOC1iBgTEUdXqW9vstlUK3YCLiCbGPJtwL+SzQzwVTq3YDqAfQv+e816Zf2edzFbpzwZEdPT+jSgLc00/H7g+tyD4AbnD5K0Fdm8Un9MRVfT+dv/3RHxYtp3NtlDh/LTYlezTUQs6xLbjFTHrFRnSJpBNl9axRKy2VjNGsbJw6yz5bn1lWSTMa5H9iyFMXWst8j/vRWS1ouIyjND8nWsyr1e1aW+jYDXehuoWRHutjLrQRqwflLSv8Dqq6De2WWfF4Blkt6TisYWrP71NBV+NXPIJsMra2eyCfXMGsbJw6yYo4ETJFVmrD2iyj4nAJekpx5uSvZEt55MAh7tZsD8VrJZe8vaPx1r1jCeVdesTiRtFhEvp/XxZM/BPrUP9Q0DJkfER0ocsy1wdUQc0NvzmhXhMQ+z+jlU0plk/6+eAo7rS2URsShdNrxFiXs9RgKn9+W8ZkW45WFmZqV5zMPMzEpz8jAzs9KcPMzMrDQnDzMzK83Jw8zMSvv/vDzDEZiW3gcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def shake(y, demonE, delta=0.1):\n",
    "    \"\"\"\n",
    "    Pass in the current demon energy as an argument.\n",
    "    delta is the size of change in y to generate, more or less.\n",
    "    randomly choose a particle, change it's position slightly (around delta)\n",
    "    return the new demon energy and a boolean (was the change accepted?)\n",
    "    \"\"\"\n",
    "    ix = int(np.random.rand()*len(y))\n",
    "    deltaY = delta*np.random.normal()\n",
    "    deltaE = deltaY*m*g\n",
    "    accept=False\n",
    "    if deltaE < demonE and (y[ix]+deltaY>0):\n",
    "        demonE -= deltaE  # take the energy from the demon, or give it if deltaE<0.\n",
    "        y[ix] += deltaY\n",
    "        accept=True\n",
    "        \n",
    "    return demonE, accept\n",
    "\n",
    "y = setup(N,L=h)\n",
    "\n",
    "acceptCount = 0\n",
    "\n",
    "demonList = []\n",
    "for i in range(M):\n",
    "    demonE,accept = shake(y, demonE, delta=0.2*h)\n",
    "    demonList.append(demonE)\n",
    "    if accept:\n",
    "        acceptCount += 1\n",
    "\n",
    "pl.title(\"Distribution of heights\")\n",
    "pl.xlabel(\"height (m)\")\n",
    "pl.ylabel(\"number in height range\")\n",
    "pl.hist(y,bins=40)\n",
    "print(100.0*acceptCount/M, \"percent accepted\")\n",
    "print(\"Averge height=%4.3fm\" % (y.sum()/len(y),))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xm4XFWZ7/HvLwnzYEACQkIISBRRkCENoaEFQTEICo2MjRDptFwfUbHFVlA0INod2osoKrmmmQKCiCgSJBBiGBxaIAmjATQxgglTAoFAGA2894+1CnZO6tTZJ2dXFZXz+zxPPWfvtXeteivnpN5aa+29liICMzOzKgxodwBmZrb6cFIxM7PKOKmYmVllnFTMzKwyTipmZlYZJxUzM6uMk4qtEkn/T9LXKqpruKRlkgbm/Vsk/VsVdef6rpc0tqr6evG635T0pKTHW/3a/Un+29mm3XFY4qRiK5H0kKQXJT0n6RlJ/yvpU5Je/3uJiE9FxJkl6/pAo3Mi4m8RsX5EvFpB7KdL+nGX+g+IiMl9rbuXcQwHTga2j4i31Tm+j6TX8gfiMkkLJV0p6R9aGeeqkDRCUhRirz2ObEc8+W9nfjte21bmpGLd+UhEbABsBUwAvgxcUPWLSBpUdZ1vEsOBpyJiUYNzHo2I9YENgNHAg8BvJe3XigArMDh/oNceP636BWqtV+scTirWUEQsjYgpwJHAWEnvAZB0saRv5u1NJP0qt2qWSPqtpAGSLiV9uF6bv8l+qfAtd5ykvwE3FcqKCebtku6Q9KykayRtnF9rH0kLizHWWkOSxgBfAY7Mr3dPPv56d1qO6zRJD0taJOkSSW/Jx2pxjJX0t9x19dXu/m0kvSU/f3Gu77Rc/weA6cAWOY6Le/g3johYGBFfB84Hziq8xnaSpud/1z9JOqJw7GJJ5+XuvWWSfi/pbZK+K+lpSQ9K2rlw/rvyv8UzkuZI+miXun4o6brcQr1d0tsbxd3g36VhXSXe00RJUyU9D7xf0lslXZv/FmYqdSv+rvCckLRt3l5L0v/Nv78nlLpp18nH6v6drsp7tO75H9RKiYg7gIXAP9U5fHI+NgTYjPTBHhFxLPA3Uqtn/Yj478Jz9gbeBXyom5c8DvhXYHNgOXBuiRhvAP4T+Gl+vffWOe0T+fF+YBtgfeAHXc7ZC3gnsB/wdUnv6uYlvw+8Jdezd475+Ij4NXAAuSUSEZ/oKfaCXwC7SFpP0nqk5HQ5sClwFHCepO0L5x8BnAZsArwM/AG4M+9fBXwHQNIawLXAjbmuzwKXSXpnoa6jgDOAjYB5wLd6EXdXdesq+Z7+JZ+/AfA74IfA88DbgLH50Z0JwDuAnYBtgaHA1/Oxun+nfXiPVoeTivXGo8DGdcr/Tvrw3yoi/h4Rv42eJ5U7PSKej4gXuzl+aUT8MSKeB74GHFFRV8gxwHciYn5ELANOBY7q0ko6IyJejIh7gHuAlZJTjuUo4NSIeC4iHgLOBo7tY3yPAgIGAwcBD0XERRGxPCLuAn4OHF44/+qImB0RLwFXAy9FxCV5fOqnQK2lMpqUQCdExCsRcRPwK+DoLnXdERHLgctIH8yNPJm/9dcexeTbXV1l3tM1EfH7iHiN9Lf1MWB8RLwQEfcDdcfHJAk4Afj3iFgSEc+RvmQclU9Zlb9T66XVtT/bmmMosKRO+beB04Eb0/9rJkXEhB7qWtCL4w8Da5C+fffVFrm+Yt2DSN9ca4pXa71A+jDuapMcU9e6hvYxvqGkb8/PkMazdpf0TOH4IODSwv4The0X6+zXYt8CWJA/qLuLt8z7LtokJ416uqurzHsq/u6H5OMLujlOl3PXBWbnv0NICbr2ZWRV/k6tl5xUrBSlq5KGkrojVpC/EZ4MnKw05nKTpJkRMYPuuxd6+oa4ZWF7OOlb5pOkbpB1C3ENJH2YlK33UdIHW7Hu5aQP42E9PLfoyRzTVsD9hboe6UUd9fwzcGdEPC9pAXBrRHywj3VCet9bShpQSCzDgT9XUHdvlHlPxd/hYtLvZxhvxLrlSs9IniQl0ndHxEq/hx7+Tq0i7v6yhiRtKOkg4ArgxxFxX51zDpK0be5+WAq8CtQ+uJ4gjTn01sclbS9pXeAbwFW5S+fPwNqSDszjBKcBaxWe9wQwosEA7E+Af5e0taT1eWMMprtv3HXlWK4EviVpA0lbAV8Aftz4mStTMlTSeODfSH39kLqn3iHpWElr5Mc/NBjjaeR2UovhS7mefYCPkH6vrdSr95T/nX8BnC5pXUnbkcau6p37GvA/wDmSNgXI/64fytuN/k6tIk4q1p1rJT1H+mb5VdKA7/HdnDsS+DWwjDRQfF5E3JyP/RdwWu5z/2IvXv9S4GJSN8rawOcgXY0GfJp0ldQjpJZL8Wqwn+WfT0m6s069F+a6fwP8FXiJNGi9Kj6bX38+qQV3ea6/rC0kLSP9u80EdgD2iYgb4fVv1vuTxgQeJf1bnMWKSbSUiHiFlEQOIH2jPw84LiIe7G1dBc9oxftUvlAijlV5T58hXRDxOOl39xPSRQn1fJl0YcBtkp4l/V3WLkZo9HdqFZHHqcysk0g6C3hbRLR8lgTrmVsqZvamlu9r2TF3E+4GjCNd6WZvQh6oN7M3uw1IXV5bkMbMzgauaWtE1i13f5mZWWXc/WVmZpXpd91fm2yySYwYMaLdYZiZdYzZs2c/GRFDej6zHyaVESNGMGvWrHaHYWbWMSQ93PNZibu/zMysMk4qZmZWGScVMzOrjJOKmZlVpqlJRdJgSVcprUD3gKQ9JG2cV32bm39ulM+VpHMlzZN0r6RdCvWMzefPlTS2UL6rpPvyc85VYb5rMzNrvWa3VL4H3BAR25EWOnoAOAWYEREjgRl5H9JEdyPz4wRgIoDSMrLjgd2B3YDxtUSUz/lk4Xljmvx+zMysgaYlFaV1v98HXABpltSIeAY4mDdWbpsMHJK3DwYuyet13wYMlrQ5abnZ6Xklt6dJS5GOycc2jIjb8uptlxTqMjOzNmhmS2Vr0gI7F0m6S9L5eX3qzSLisXzO47yx4t5QVlzRbWEua1S+sE75SiSdIGmWpFmLFy/u49syM7PuNDOpDAJ2ASZGxM6kdSdOKZ6QWxhNn3wsIiZFxKiIGDVkSKmbQs3MbBU08476hcDCiLg9719FSipPSNo8Ih7LXViL8vFHWHGZ0GG57BFgny7lt+TyYXXOb5oRp1y3UtlDEw5s5kuamXWUprVUIuJxYIGk2qpr+5HW8p4C1K7gGssbU1hPAY7LV4GNBpbmbrJpwP6SNsoD9PsD0/KxZyWNzld9HYenwzYza6tmz/31WeAySWuSllw9npTIrpQ0DngYOCKfOxX4MGkp0BfyuUTEEklnkpZbBfhGRCzJ258mLTm7DnB9fpiZWZs0NalExN3AqDqH9qtzbgAndlPPhdRZ+zsiZgHv6WOYZmZWEd9Rb2ZmlXFSMTOzyjipmJlZZZxUzMysMk4qZmZWGScVMzOrjJOKmZlVxknFzMwq46RiZmaVcVIxM7PKOKmYmVllnFTMzKwyTipmZlYZJxUzM6uMk4qZmVXGScXMzCrjpGJmZpVxUjEzs8o4qZiZWWWcVMzMrDJOKmZmVhknFTMzq4yTipmZVcZJxczMKuOkYmZmlXFSMTOzyjQ1qUh6SNJ9ku6WNCuXbSxpuqS5+edGuVySzpU0T9K9knYp1DM2nz9X0thC+a65/nn5uWrm+zEzs8Za0VJ5f0TsFBGj8v4pwIyIGAnMyPsABwAj8+MEYCKkJASMB3YHdgPG1xJRPueTheeNaf7bMTOz7rSj++tgYHLengwcUii/JJLbgMGSNgc+BEyPiCUR8TQwHRiTj20YEbdFRACXFOoyM7M2GNTk+gO4UVIAP4qIScBmEfFYPv44sFneHgosKDx3YS5rVL6wTvlKJJ1Aav0wfPjwvryflYw45bqVyh6acGClr2Fm1il6TCqS1gNejIjXJL0D2A64PiL+XqL+vSLiEUmbAtMlPVg8GBGRE05T5WQ2CWDUqFFNfz0zs/6qTPfXb4C1JQ0FbgSOBS4uU3lEPJJ/LgKuJo2JPJG7rsg/F+XTHwG2LDx9WC5rVD6sTrmZmbVJmaSiiHgBOBQ4LyIOB97d45Ok9SRtUNsG9gf+CEwBaldwjQWuydtTgOPyVWCjgaW5m2wasL+kjfIA/f7AtHzsWUmj81VfxxXqMjOzNigzpiJJewDHAONy2cASz9sMuDpf5TsIuDwibpA0E7hS0jjgYeCIfP5U4MPAPOAF4HiAiFgi6UxgZj7vGxGxJG9/mtRqWge4Pj/MzKxNyiSVzwOnAldHxBxJ2wA39/SkiJgPvLdO+VPAfnXKAzixm7ouBC6sUz4LeE9PsZiZWWv0mFQi4lbgVknr5v35wOeaHZiZmXWeHsdUJO0h6X7gwbz/XknnNT0yMzPrOGUG6r9LugHxKYCIuAd4XzODMjOzzlTqjvqIWNCl6NUmxGJmZh2uzED9Akn/CISkNYCTgAeaG5aZmXWiMi2VT5GuyhpKurlwJ7q5SsvMzPq3Mld/PUm6R8XMzKyhMnN/nVuneCkwKyJ8B7uZmb2uTPfX2qQur7n5sSNpnq1xkr7bxNjMzKzDlBmo3xHYMyJeBZA0EfgtsBdwXxNjMzOzDlOmpbIRsH5hfz1g45xkXm5KVGZm1pHKtFT+G7hb0i2ASDc+/meeefjXTYzNzMw6TJmrvy6QNJW0FgrAVyLi0bz9H02LzMzMOk7ZNeoHAIuBp4FtJXmaFjMzW0mZS4rPAo4E5gCv5eIgrQhpZmb2ujJjKocA74wID8qbmVlDZbq/5gNrNDsQMzPrfGVaKi+Qrv6aQeES4ojwQl1mZraCMkllSn6YmZk1VOaS4smtCMTMzDpfmau/RgL/BWxPmgcMgIjYpolxdbQRp1y3UtlDEw5sQyRmZq1VZqD+ImAisBx4P3AJ8ONmBmVmZp2pTFJZJyJmAIqIhyPidMBfu83MbCVlBupfljQAmCvpM6TVH9fv4TlmZtYPlWmpnASsC3wO2BU4FhjbzKDMzKwzlbn6a2beXAYcDyBpeDODMjOzztSwpSJpD0mHSdo07+8o6XLg92VfQNJASXdJ+lXe31rS7ZLmSfqppDVz+Vp5f14+PqJQx6m5/E+SPlQoH5PL5kk6pVfv3MzMKtdtUpH0beBC4GPAdZK+CdwI3A6M7MVrnAQ8UNg/CzgnIrYlzXo8LpePA57O5efk85C0PXAU8G5gDHBeTlQDgR8CB5Audz46n2tmZm3SqKVyILBzRBwN7A98HhgdEd+LiJfKVC5pWK7n/LwvYF/gqnzKZNKElQAH533y8f3y+QcDV0TEyxHxV2AeaW2X3YB5ETE/Il4BrsjnmplZmzRKKi/VkkdEPA3MjYiHeln/d4Ev8caU+W8FnomI5Xl/ITA0bw8FFuTXWw4szee/Xt7lOd2Vr0TSCZJmSZq1ePHiXr4FMzMrq9FA/TaSinN+bV3cj4iPNqpY0kHAooiYLWmfvoXZNxExCZgEMGrUqGhnLGZmq7NGSaVrV9LZvax7T+Cjkj5Mmt5lQ+B7wGBJg3JrZBjpvhfyzy2BhZIGAW8BniqU1xSf0125mZm1QbdJJSJu7UvFEXEqcCpAbql8MSKOkfQz4DDSGMhY4Jr8lCl5/w/5+E0REbl1dLmk7wBbkC4SuAMQMFLS1qRkchTwL32J2czM+qbMHfVV+zJwRb6a7C7gglx+AXCppHnAElKSICLmSLoSuJ80/9iJEfEqQL7DfxowELgwIua09J2YmdkKWpJUIuIW4Ja8PZ905VbXc14CDu/m+d8CvlWnfCowtcJQzcysD3qcpkXSDq0IxMzMOl+Zub/Ok3SHpE9LekvTIzIzs47VY1KJiH8CjiFdaTVb0uWSPtj0yMzMrOOUaakQEXOB00iD7HsD50p6UNKhzQzOzMw6S5kxlR0lnUOav2tf4CMR8a68fU6T4zMzsw5S5uqv75Mu9/1KRLxYK4yIRyWd1rTIzMys45RZT2XvBscurTYcMzPrZD0mFUn3AV3ny1oKzAK+GRFPNSOw1c2IU65bqeyhCQe2IRIzs+Yp0/11PfAqcHneP4q0vPDjwMXAR5oSmZmZdZwySeUDEbFLYf8+SXdGxC6SPt6swMzMrPOUuaR4oKTXp1WR9A+kubYgzcVlZmYGlGupjAMukrR+3n8OGCdpPeC/mhaZmZl1nIZJRdIAYJuI2KE2RUtELC2ccmUzgzMzs87SsPsrIl4jLQdMRCztklDMzMxWUGZM5deSvihpS0kb1x5Nj8zMzDpOmTGVI/PPEwtlAWxTfThmZtbJytxRv3UrAjEzs85XZkLJdSWdJmlS3h8p6aDmh2ZmZp2mzJjKRcArwD/m/UeAbzYtIjMz61hlksrbI+K/gb8DRMQLgJoalZmZdaQySeUVSeuQJ5WU9Hbg5aZGZWZmHanM1V/jgRuALSVdBuwJfKKZQZmZWWcqc/XXdEl3AqNJ3V4nRcSTTY/MzMw6TpmWCsDawNP5/O0lERG/aV5YZmbWicos0nUW6QbIOcBruTgAJxUzM1tBmZbKIcA7I6JXg/OS1iYlnrXy61wVEeMlbQ1cAbwVmA0cGxGvSFoLuATYFXgKODIiHsp1nUqaLflV4HMRMS2XjwG+R5qK//yImNCbGNvNq0Ga2eqmzNVf84E1VqHul4F9I+K9wE7AGEmjgbOAcyJiW1KX2rh8/jjg6Vx+Tj4PSduTVpt8NzAGOE/SQEkDgR8CBwDbA0fnc83MrE3KtFReAO6WNIPCpcQR8blGT4qIAJbl3TXyI4B9gX/J5ZOB04GJwMF5G+Aq4AeSlMuvyC2lv0qaB9QWDZsXEfMBJF2Rz72/xHsyM7MmKJNUpuRHr+XWxGxgW1Kr4i/AMxFRWzFyITA0bw8FFgBExHJJS0ldZEOB2wrVFp+zoEv57qsSp5mZVaPbpCJpw4h4NiIm1zk2vEzlEfEqsJOkwcDVwHarHGkfSDoBOAFg+PBSoZuZ2SpoNKZyS20jd30V/bI3LxIRzwA3A3sAgyXVktkw0lxi5J9b5tcbBLyFNGD/enmX53RXXu/1J0XEqIgYNWTIkN6EbmZmvdAoqRTn9+q6KFePc39JGpJbKORpXj4IPEBKLofl08YC1+TtKXmffPymPC4zBThK0lr5yrGRwB3ATGCkpK0lrUkazF+lbjozM6tGozGV6Ga73n49mwOT87jKAODKiPiVpPuBKyR9E7gLuCCffwFwaR6IX0JKEkTEHElXkgbglwMn5m41JH0GmEa6pPjCiJhTIi4zM2uSRkllU0lfILVKatvk/R77kCLiXmDnOuXzeePqrWL5S8Dh3dT1LeBbdcqnAlN7isXMzFqjUVL5H2CDOtsA5zctIjMz61jdJpWIOKOVgZiZWecrc0e9mZlZKU4qZmZWGScVMzOrTJmp79cCPgaMKJ4fEd9oXlhmZtaJysz9dQ2wlDSHl9emNzOzbpVJKsMiYkzTIzEzs45XJqn8r6QdIuK+pkdjXrjLzDpamaSyF/AJSX8ldX+JtFzKjk2NzMzMOk6ZpHJA06MwM7PVQo/rqQDPtTAeMzPrYI1aKpcDB5Gu+gpWnO4+gG2aGJeZmXWgRnN/HZR/bt26cMzMrJP5jnozM6uMk4qZmVXGScXMzCrTMKlIGijpwVYFY2Zmna1hUslrwf9J0vAWxWNmZh2szM2PGwFzJN0BPF8rjIiPNi0qMzPrSGWSyteaHoWZma0WekwqEXGrpK2AkRHxa0nrAgObH5rVeJJJM+sUPV79JemTwFXAj3LRUOCXzQzKzMw6U5lLik8E9gSeBYiIucCmzQzKzMw6U5mk8nJEvFLbkTSINPeXmZnZCsoklVslfQVYR9IHgZ8B1zY3LDMz60RlksopwGLgPuD/AFOB03p6kqQtJd0s6X5JcySdlMs3ljRd0tz8c6NcLknnSpon6V5JuxTqGpvPnytpbKF8V0n35eecK0krR2JmZq3SY1KJiNeAycCZwBnA5Igo0/21HDg5IrYHRgMnStqelKRmRMRIYEbeh7QY2Mj8OAGYCCkJAeOB3YHdgPG1RJTP+WTheWNKxGVmZk1S5uqvA4G/AOcCPwDmSepxNciIeCwi7szbzwEPkK4cO5iUpMg/D8nbBwOXRHIbMFjS5sCHgOkRsSQingamA2PysQ0j4rac5C4p1GVmZm1Q5ubHs4H3R8Q8AElvB64Dri/7IpJGADsDtwObRcRj+dDjwGZ5eyiwoPC0hbmsUfnCOuX1Xv8EUuuH4cM944yZWbOUGVN5rpZQsvn0YolhSesDPwc+n5cnfl1uYTT9SrKImBQRoyJi1JAhQ5r9cmZm/VajNeoPzZuzJE0FriQlgMOBmWUql7QGKaFcFhG/yMVPSNo8Ih7LXViLcvkjwJaFpw/LZY8A+3QpvyWXD6tzvpmZtUmj7q+PFLafAPbO24uBdXqqOF+JdQHwQER8p3BoCjAWmJB/XlMo/4ykK0iD8ktz4pkG/GdhcH5/4NSIWCLpWUmjSd1qxwHf7ymu1YWnbjGzN6NGa9Qf38e69wSOBe6TdHcu+wopmVwpaRzwMHBEPjYV+DAwD3gBOD7HsUTSmbzROvpGRCzJ258GLiYluevpxTiPmZlVr8eBeklbA58FRhTP72nq+4j4HdDdfSP71Tk/SFPC1KvrQuDCOuWzgPc0isPMzFqnzNVfvyR1Y10LvNbccMzMrJOVSSovRcS5TY/EzMw6Xpmk8j1J44EbgZdrhbUbG83MzGrKJJUdSAPu+/JG91fkfTMzs9eVSSqHA9sUp783MzOrp8wd9X8EBjc7EDMz63xlWiqDgQclzWTFMZWGlxRb6/mGSDNrtzJJZXzTozAzs9VCj0klIm5tRSBmZtb5ytxR/xxvzCS8JrAG8HxEbNjMwMzMrPOUaalsUNvOk0QeTFrJ0czMbAVlrv56XV6V8Zek1RjNzMxWUKb769DC7gBgFPBS0yIyM7OOVebqr+K6KsuBh0hdYGZmZisoM6bS13VVzMysn2i0nPDXGzwvIuLMJsRjFat3QyT4pkgza45GLZXn65StB4wD3go4qZiZ2QoaLSd8dm1b0gbASaQlfq8Azu7ueWZm1n81HFORtDHwBeAYYDKwS0Q83YrAzMys8zQaU/k2cCgwCdghIpa1LCozM+tIjW5+PBnYAjgNeFTSs/nxnKRnWxOemZl1kkZjKr262946i6fJN7NmcOIwM7PKOKmYmVllnFTMzKwyTipmZlaZpiUVSRdKWiTpj4WyjSVNlzQ3/9wol0vSuZLmSbpX0i6F54zN58+VNLZQvquk+/Jzzs1rvZiZWRuVmaV4VV0M/AC4pFB2CjAjIiZIOiXvfxk4ABiZH7sDE4Hd882X40nT7QcwW9KUfAPmROCTwO3AVGAMcH0T389qz1eEmVlfNa2lEhG/AZZ0KT6YdGc++echhfJL8iJgtwGDJW1OWgxsekQsyYlkOjAmH9swIm6LiCAlrkMwM7O2avWYymYR8VjefhzYLG8PBRYUzluYyxqVL6xTXpekEyTNkjRr8eLFfXsHZmbWrbYN1OcWRrTotSZFxKiIGDVkyJBWvKSZWb/U6qTyRO66Iv9clMsfAbYsnDcslzUqH1an3MzM2qjVSWUKULuCayxwTaH8uHwV2Ghgae4mmwbsL2mjfKXY/sC0fOxZSaPzVV/HFeoyM7M2adrVX5J+AuwDbCJpIekqrgnAlZLGAQ8DR+TTpwIfBuYBL5DWbSEilkg6E5iZz/tGRNQG/z9NusJsHdJVX77yy8yszZSGNvqPUaNGxaxZs1bpud0tzdvf+DJjs/5F0uyIGFXmXN9Rb2ZmlXFSMTOzyjTzjnpbTfnOezPrjlsqZmZWGScVMzOrjJOKmZlVxmMqVgmPs5gZuKViZmYVclIxM7PKOKmYmVllPKZiTeNxFrP+xy0VMzOrjFsq1lJuvZit3txSMTOzyjipmJlZZdz9ZW3nLjGz1YdbKmZmVhm3VOxNya0Xs87kloqZmVXGLRXrGG69mL35uaViZmaVcUvFOlq91ks9btGYtYZbKmZmVhm3VKxf8HiMWWs4qVi/VbbrDJyAzMpyUjErwWM3ZuV0fFKRNAb4HjAQOD8iJrQ5JOvHnHysv+vopCJpIPBD4IPAQmCmpCkRcX97IzNrrDddb2U4SdmbRUcnFWA3YF5EzAeQdAVwMOCkYv1K1UmqGZz4+odOTypDgQWF/YXA7l1PknQCcELeXSbpT6v4epsAT67iczuZ33f/0pT3rbOqrrFy/n13b6uylXV6UiklIiYBk/paj6RZETGqgpA6it93/+L33b9U/b47/ebHR4AtC/vDcpmZmbVBpyeVmcBISVtLWhM4CpjS5pjMzPqtju7+iojlkj4DTCNdUnxhRMxp4kv2uQutQ/l99y9+3/1Lpe9bEVFlfWZm1o91eveXmZm9iTipmJlZZZxUSpA0RtKfJM2TdEq742kVSRdKWiTpj+2OpVUkbSnpZkn3S5oj6aR2x9QqktaWdIeke/J7P6PdMbWKpIGS7pL0q3bH0kqSHpJ0n6S7Jc2qpE6PqTSWp4L5M4WpYICj+8NUMJLeBywDLomI97Q7nlaQtDmweUTcKWkDYDZwSD/5fQtYLyKWSVoD+B1wUkTc1ubQmk7SF4BRwIYRcVC742kVSQ8BoyKisps+3VLp2etTwUTEK0BtKpjVXkT8BljS7jhaKSIei4g78/ZzwAOkmRtWe5Esy7tr5Mdq/61T0jDgQOD8dseyOnBS6Vm9qWD6xYdMfydpBLAzcHt7I2md3A10N7AImB4R/eG9fxf4EvBauwNpgwBulDQ7T2fVZ04qZnVIWh/4OfD5iHi23fG0SkS8GhE7kWan2E3Sat3tKekgYFFEzG53LG2yV0TsAhwAnJi7vPvESaVnngqmn8njCT8HLouIX7Q7nnaIiGeAm4Ex7Y6lyfYEPprHFq4A9pX04/aG1DoR8Uj+uQi4mtTd3ydOKj3zVDD9SB6svgB4ICK+0+54WknSEEmD8/Y6pItTHmxvVM0VEadGxLCIGEH6v31TRHy8zWG1hKT18sUoSFoP2B/o85WeTio9iIjlQG0qmAeAK5s8FcybhqT48wYdAAAFZUlEQVSfAH8A3ilpoaRx7Y6pBfYEjiV9Y707Pz7c7qBaZHPgZkn3kr5MTY+IfnWJbT+zGfA7SfcAdwDXRcQNfa3UlxSbmVll3FIxM7PKOKmYmVllnFTMzKwyTipmZlYZJxUzsze5Kid3lbSTpD/kSUPvlXRk4djWkm7Pk+f+NN9G0StOKrbakvRq4bLgu98sM0xLuiXPen2PpJmSdmp3TF1J+ryk4/L2xZIOy9tXSBrZ3uj6pYup7kbUF4DjIuLduc7v1u5PAs4CzomIbYGngV7fRuCkYquzFyNip8JjQl8rlFTVEtzHRMR7gfOAb1dUZyXye/xX4PI6hyeS5smyFqo3uaukt0u6Ic/b9VtJ25Ws688RMTdvP0qa521IvvF3X+CqfOpk4JDexuqkYv1OXkPiDEl35rUktsvl6+Vuhjvy2hoH5/JPSJoi6SZghqQBks6T9KCk6ZKmSjpM0r6Sfll4nQ9KurqHcP5AYYJSSRMlzeq6nkmDmIfkGOZIOl/Sw5I2ycc+nt/L3ZJ+lCeLHJhbHn/M9fx7nZj2Be7MN/529VvgAxUmV1t1k4DPRsSuwBdJX1B6RdJuwJrAX4C3As8Ufu+rNHmuk4qtztbp0v11ZOHYk3kivYmk/5AAXyVN07Eb8H7g23n6CoBdgMMiYm/gUGAEsD3p7vs98jk3A9tJGpL3jwcu7CHGMcAvC/tfjYhRwI7A3pJ27CHm8Tnmd5O+YQ4HkPQu4EhgzzxB5KvAMcBOwNCIeE9E7ABcVCemPUnryKwkIl4D5gHv7eF9WRPlCU//EfiZ0qzSPyLNiICkQ/OXhq6PaV3q2By4FDg+/14r4W8btjp7MX+g1lObKHI2KUlAmvvoo5JqH9hrkz+kSVOW1Lof9gJ+lv8jPi7pZkjrkUi6FPi4pItIyea4bl7/sjwIuj7pg77mCKUpyAeRPiS2B+5tEPNewD/n179B0tO5fD9gV2Bm6tVgHVI3x7XANpK+D1wH3Fgnts1JUxJ1ZxGwBd0kHmuJAaRWxUp/33kS1IYToUrakPT7/2phEbangMGSBuXWyipNnuuWivVXL+efr/LGlysBHyuMwQyPiNqH6/Ml670I+DhwNCnx1OtCgtRq2IbUb/19SFfekFog+0XEjqT/9Gv3EHN3BEwuvJd3RsTpEfE0qZVxC/Ap6i9M9WKX1+1q7XyOtUlejuGvkg6HNBGqpFKtx/xl5mrSiq618RMizdl1M3BYLhoLXNPb2JxUzN4wDfhsHrBE0s7dnPd74GN5bGUzYJ/agTzw+ShwGvW7liicG8DXgNF5jGRDUvJamus9oETMvweOyPHuD2yUy2cAh0naNB/bWNJWebxlQET8PMe4S506HwC2bfCa76CC2WytPNWf3PUYYJzShJBzKL8i7RHA+4BPFLqGay2eLwNfkDSPNMZyQW9jdfeXrc7Wyf3NNTdERKPLis8krQJ4r6QBwF+BeuuV/5zUvXQ/aVXQO4GlheOXAUMKrZxuRcSLks4G/iMixkm6izTd/AJSwujJGcBPJB1L+tB5HHguIp6UdBppVb8BwN+BE0ktjItyGcCpdeq8ntTXXjOI3ErKye7FiHi8RGxWkYg4uptDvb7MOCJ+DNRdMyYi5tPHNVU8S7HZKpC0fkQsk/RW0rThe9Y+aCX9ALgrInr9LW8V4lgLeDUilkvaA5jYYBypN/VeTbp0+C+kafCPjYj789Viz7bivVlnckvFbNX8SumGsTWBMwsJZTapC+vkFsUxHLgytzxeAT5ZUb2nkAbsrwFujoj7c/kzrNiKMVuBWypmZlYZD9SbmVllnFTMzKwyTipmZlYZJxUzM6uMk4qZmVXm/wPfFWFqK2RZQwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#\n",
    "# Build a histogram of Demon Energies\n",
    "#\n",
    "\n",
    "pl.title(\"Distribution of Demon Energies\")\n",
    "pl.xlabel(\"Energy Ranges (J)\")\n",
    "pl.ylabel(\"Number in Energy Ranges\")\n",
    "ns, bins, patches = pl.hist(demonList, bins=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonic Thermometer\n",
    "====================\n",
    "\n",
    "You can easily see that the demon acts like an small thermometer. According to the Maxwell-Boltzmann distribution the energy distribution of the demon's energy should go like:\n",
    "\n",
    "$$P(E) = P_0 e^{-E/k_B T}$$\n",
    "\n",
    "Where $P_0$ is the basically the probability of having an energy of zero. (Actually, maybe a better way to think of it is as a normalization constant that's determined by the requirement that the total probability to have *any* energy is 1.0). The histogram of demon energies tells us the number of times the demon have various values of energy during the calculation. This is proportional to the probability that the demon had various energies. We can fit that probability to an exponential curve (or the log of the probability to a straight line) and from the slope of the line deduce the temperature!\n",
    "\n",
    "See below how the code does exactly this.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "slope= -1.9387548999149866e+20 +/- 4.470166272915833e+17\n",
      "intercept= 11.058280298449063 +/- 0.003256417693895689\n",
      "Temperature= 308.85925561598225 +/- 0.7121334561645698 K\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x113268ad0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VGX2wPHvyaTRlRiVTqRJiSBFCQgCFkQBXcGuP3V1g66KriLq2ntjlRUbrC7YFhWwV4pEFGIBREJROkpRMNJNz/n9cW/iEEJyk0xJMufzPPOQuXPn3jOgOfO284qqYowxJnJFhTsAY4wx4WWJwBhjIpwlAmOMiXCWCIwxJsJZIjDGmAhnicAYYyKcJQJjqhkRuUhEZgbwestFZID78z0i8moAr/1PEXkhUNcz4WGJIIKJyAYRyRKRPSKyU0QWiMhVIlLt/7sQkTQRyRaRvX6P98MdV3lEZIqI5Lp/53tEZJmIPCwijYrOUdXXVPVUj9d6oLzzVLWzqqZVMXREZICIbCpx7YdU9cqqXtuEV7X/H94E3TBVbQC0Ah4BbgFeDG9Inl2rqvX9HsMCfQMRiQ70NYHH3L/zROByoDcwX0TqBfImQYrd1EKWCAwAqrpLVd8DzgMuFZEuACISJyLjROQnEflVRJ4XkTruawNEZJOIjBWRbSKyVUTOEpHTRWSViPwuIv8suod7rfEissV9jBeRuBLXusnvWpdX5rOUdy2Pn+kWEfkFmOweH+teZ4uIXCkiKiJtRaSXew2f3/XPFpHvPfydZ6vqt8BwIAEnKSAil4nIl+7PIiJPup9jt4hkiEgXEUkFLgLG+reG3FbeLSKyFNgnItHusZP9bh0vIm+4LZLFItLVL3YVkbZ+z6eIyANukvoYaOrXAmtasqtJRIa7XVE73VZbR7/XNojIGBFZKiK73BjivfybmuCyRGD2o6rfAJuAfu6hR4D2QDegLdAMuMvvLUcC8X7H/wNcDPRwr3GniCS5596O8+23G9AVOA64o8S1GrnXugJ4RkQOreRHKetaXj5TY5xWUqqInAbcCJzsnj+g6ET3F3km4N+VcwnwstdAVXUPMIs//879nQr0d+NtBJwLZKrqJOA1nNZFydbQBcAZwCGqml/KNc8Eprmf8X/AOyISU06M+4AhwBa/FtgW/3NEpD0wFbgBp7XzEfC+iMT6nXYucBqQBBwDXFbWfU1oWCIwpdkCNBYRAVKBf6jq7+4vrIeA8/3OzQMeVNU84HXgMODfqrpHVZcDK3B+6YPzDfY+Vd2mqtuBe3F+afpf6z5VzVPVj4C9QIcy4nzK/eZZ9Li/vGt5/EyFwN2qmqOqWTi/vCar6nJV/QO4p0QcL+EkP0SkMTAY5xdsRWzB+cVcUh7QADgaEFVdqapby7nWU6r6sxt7aRap6nT33+wJnETeu4LxluY84ENVneVeexxQB+hTIrYtqvo78D5OMjZhZn2IpjTNgN9xvtXVBRY5vz8BEMDnd26mqha4Pxf94vnV7/UsoL77c1Ngo99rG91j/tfy/wb7h997SzNaVQ82Y+Vg1/Lymbararbf86bAQr/nP5e416vASrf75FzgCw+/rEsq+jvfj6p+JiJPA88ArUTkLWCMqu4u41ol4zvo66pa6A4ANy3jfK/2+/d1r/0zzmcr8ovfz38E6L6miqxFYPYjIr1w/sf9EvgN5xd5Z1U9xH00UtWyfjmXZQtOd0uRlu6xUPLymUqW5N0KNPd73sL/RVXdDKQDZ+O0cF6pSEAiUh+n2+mL0l5X1adUtQfQCaeL6OaDxHmw+Esqjl+cGWLN+fPf4Q+cRFnkyApcd79/X7f11QLYXM77TJhZIjAAiEhDERmK073zqqpmqGohTp//kyJyuHteMxEZXMnbTAXuEJFEETkMp18+YHPavajkZ3oTuFxEOopIXeDOUs55GRgLJANveYnFHbTuAbwD7MAdmC5xTi8ROd7tw98HZON0XYHT8jrKy71K6OEOaEfj9OfnAF+5ry0BLhQRnzs2cqLf+34FEsRvqmsJbwJniMhJbrw3uddeUIkYTQhZIjDvi8genO6C23H6jP1n69wCrAG+EpHdwGzK7rcvywM4XSxLgQxgsXussp6W/dcRLPL4vgp9JlX9GHgKmFv0PvelHL/T3sb5Nvy2O45QlrHu33kmTgJZBPRxB2RLaoiTuHbgdLtkAo+7r70IdHLHR94p557+3sXpz9+B04I52+3TB7geGAbsxBnTKb6uqv6Ak8zXuffcr1tHVX/EGSuZgNPyGoYzPTm3ArGZMBDbmMaYinGnRC4D4vzHIURkLTBKVWeHLThjKsFaBMZ4ICJ/cbtyDgUeBd4vkQRG4PShfxauGI2pLEsExngzCtgGrAUKgKuLXhCRNOA54Bp3DMKYGsW6howxJsJZi8AYYyJcjVhQdthhh2nr1q3DHYYxxtQoixYt+k1VE8s7r0YkgtatW7Nw4cLyTzTGGFNMRDaWf5Z1DRljTMSzRGCMMRHOEoExxkS4GjFGYIwJjLy8PDZt2kR2dnb5J5saIz4+nubNmxMTU+a2EgdlicCYCLJp0yYaNGhA69at8SvDbWowVSUzM5NNmzaRlJRU/htKYV1DxkSQ7OxsEhISLAnUIiJCQkJClVp5tToRZExKJ23ww2RMSg93KMZUG5YEap+q/puW2TXkbiw9FGcv1aY4G3osw9mObnmV7hxkGZPSaTtqEB3JI3dmLBnMITk1JdxhGWNMtXPQFoGI3AvMB1KAr4GJOBtP5AOPiMgsETkmJFFWQuaMNOLJJpoC4shm99QPwx2SMQbw+Xx069aNzp0707VrV/71r39RWFg9avUVxVb0eOSRR4J+zyuvvJIVK1YE/T5lKatF8I2q3n2Q155wd3dqGYSYAiJhxACyZ8YTRw5RKL3Tn4DblRWH9WfbJ4tJGDHAWgjGhEGdOnVYsmQJANu2bePCCy9k9+7d3HvvvWGObP/YAiU/P5/o6IP/qn3hhYNtux1CqlrmA0gu75xgP3r06KGVsXTiAp176kP64z2vqZ57ripoIWgBon9QR5dOXFCp6xpTU61YsaLC71mwQPWhh5w/A6FevXr7PV+7dq02btxYCwsLNT8/X8eMGaM9e/bU5ORkff7551VVde7cudq/f38dPny4JiUl6S233KKvvvqq9urVS7t06aJr1qxRVdX169frwIEDNTk5WQcNGqQbN25UVdVLL71Ur7vuOk1JSdGkpCSdNm2ap9iKtGrVSu+66y499thjtUuXLrpy5UpVVd27d69efvnl2qtXL+3WrZu+8847qqo6efJkHTZsmA4cOFD79++vBQUFevXVV2uHDh305JNP1iFDhhTHcOKJJ+q3336rqqqffvqp9u7dW4899lgdOXKk7tmzR1VVb7nlFu3YsaMmJyfrTTfdVGqMpf3bAgvVw+9YL4ngC+Ab4O9AIy8XDfSjsomgpG9SRmsBFCeEdW1OUt25MyDXNqYmqGgiWLBAtU4dVZ/P+TMQyaC0X7aNGjXSX375RSdOnKj333+/qqpmZ2drjx49dN26dTp37lxt1KiRbtmyRbOzs7Vp06Z61113qarq+PHj9frrr1dV1aFDh+qUKVNUVfXFF1/UM888U1WdRDBy5EgtKCjQ5cuXa5s2bUqNLSoqSrt27Vr8eP3111XVSQRPPfWUqqo+88wzesUVV6iq6m233aavvPKKqqru2LFD27Vrp3v37tXJkydrs2bNNDMzU1VVp02bpkOGDNGCggLdunWrHnLIIQckgu3bt2u/fv107969qqr6yCOP6L333qu//fabtm/fXgsLC4vvU5qqJIJy1xGoaj8RaQf8FVgkIt8Ak1V1VrBaKcESf9n5ZKf/h1hyEZSktXMgKYlf+p/D6j1NOOS8wdZdZIyftDTIzYWCAufPtDRICeL/IjNnzmTp0qVMnz4dgF27drF69WpiY2Pp1asXTZo0AaBNmzaceuqpACQnJzN37lwA0tPTeeuttwC45JJLGDt2bPG1zzrrLKKioujUqRO//vprqfcvq2vo7LPPBqBHjx7F95g5cybvvfce48aNA5zpuT/99BMAp5xyCo0bNwbgyy+/5JxzziEqKoojjzySgQMHHnD9r776ihUrVtC3b18AcnNzSUlJoVGjRsTHx3PFFVcwdOhQhg4d6unvsiI8LShT1dUicgfOxuNPAceKM1/pn6r6VsCjCpLk1BQymEPmjDRnjKBnHLsvG82R707iCCDvs4dYnvMxna87KdyhGlMtDBgAsbFOEoiNdZ4H2rp16/D5fBx++OGoKhMmTGDw4MH7nZOWlkZcXFzx86ioqOLnUVFR5OfnUx7/92slNuQqer/P5yu+n6oyY8YMOnTosN+5X3/9NfXq1avQ9VWVU045halTpx7w2jfffMOcOXOYPn06Tz/9NJ99FtgdUctdRyAix4jIk8BKYBAwTFU7uj8/GdBoQiA5NYUBn97mfPPv3p3FTc4gnygEiCWPdjcNh8cfh337wh2qMWGXkgJz5sD99zt/Bro1sH37dq666iquvfZaRITBgwfz3HPPkZeXB8CqVavYV4H/F/v06cPrr78OwGuvvUa/fv0CG3AJgwcPZsKECcWJ5bvvviv1vL59+zJjxgwKCwv59ddfSUtLO+Cc3r17M3/+fNasWQPAvn37WLVqFXv37mXXrl2cfvrpPPnkk3z//fcB/xxeWgQTgBdxvv1nFR1U1S1uK6FGSxgxgNyZcSi5FBBNXrtkYseOhXHj2Nr/PFbvSODQc0+1LiMTsVJSApsAsrKy6NatG3l5eURHR3PJJZdw4403As5Uyg0bNtC9e3dUlcTERN555x3P154wYQKXX345jz/+OImJiUyePLlSsRU57bTTypxCeuedd3LDDTdwzDHHUFhYSFJSEh988MEB540YMYI5c+bQqVMnWrRoQffu3WnUqNF+5yQmJjJlyhQuuOACcnJyAHjggQdo0KABZ555JtnZ2agqTzzxRIU+kxc1Ys/inj17ajA3psmYlP5nd1FqCixYwJ4rr6fByoUokEcMq5+eSedrBgQtBmNCYeXKlXTs2DHcYUSkvXv3Ur9+fTIzMznuuOOYP38+Rx55ZMCuX9q/rYgsUtWe5b233BaBiGQAJbPFLpzxggdUNbMCsVZLyakp4P+Nv08fFrU4m34rF+OjkFjyaHvjcNCHWCZd+O29dFuHYIypkKFDh7Jz505yc3O58847A5oEqspL19DHQAHwP/f5+UBd4BdgCjAsKJGFWcKIAeTMjCPG7TIqaHkUcdddRycEBXJmxlvZCmOMZ6WNC1QXXhLByara3e95hogsVtXuInJxsAILtwNmGP2tN0t6XUnXRf8lCqhDFrH/fgwue8OZTmGMMTWUl0TgE5HjVPUbABHpBfjc18qfs1WDlewy8qVeSdaoqcSSgwAdVrwDHTqwqd8FrN1al8bnnGQtBGNMjeMlEVwBTBaR+u7zPcAVIlIPeDhokVVD+7cSTiS5xS7+uPpGmr/yMM2AvNn3k1Ewm+SrgztlzRhjAqm8MtRRwFGqmiwijQBUdZffKW+W8d7/4pSw3qaqXdxjjYE3gNbABuBcVd1RlQ8QaiVbCd90+I5+G+90B5VzaTPmLGjwbzL2tibz7S9sUNkYU+2VuaBMVQuBse7Pu0okgfJMAU4rcexWYI6qtgPmuM9rtIQRA8khjjx85BCLNG4Ml1xC56v702/mHbQZdZJtjGOMn5Klnjds2MDChQsZPXo04AyqLliwIMxRRhYvXUOzRWQMzjf54iV+qvp7WW9S1Xki0rrE4TOBAe7PLwFpwC3eQq2eDhhUvvJ4lnX/Pzp//xpRKHXIIuqZp+DK4yGqVm8IZ4wnpdXzad26NT17OtPd09LSqF+/Pn369AlHeBGp3AVlIrK+lMOqqkeVe3EnEXzg1zW0U1UPcX8WYEfR81LemwqkArRs2bLHxo0by7tdtZExKZ02o04ijhwEJQqFLl3Y2PcC1q+DhJEDrbvIhEV1WFBWv3599u7du9+xtLQ0xo0bx9NPP03v3r3x+XwkJiYyYcKEoJeJqC2CuqBMVZOqEFtZ11UROWgWUtVJwCRwVhYHI4Zg2a+V8Jd+JDf6mZx/3EqribfTEsidFUsGc0lOtW88JoxuuAECvAkL3brB+PFlnuJfxiEpKYm33367+LXWrVtz1VVXUb9+fcaMGRPY2MxBeVlZXBe4EWipqqluSeoOqnpgQY3y/SoiTVR1q4g0AbZV4ho1QslB5fTJ6+g36y58FBJHLq3GnANNnidjawKZMz63QWUTMYKxC5ipGi9jBJOBRUDR19fNwDSgMongPeBS4BH3z3crcY0aKWHkIHJmPUgMuSg+4uKA4cPpjKAIOTPjbKWyCa1yvrmbyOFl9LKNqj4G5AGo6h+AlPcmEZkKpAMdRGSTiFyBkwBOEZHVwMnu84iQnJrC2olzmH/q/fw4MY24LRv4ofPZCIqPQuqQReGk/zj7pxkTwRo0aMCePXvCHUZE8dIiyBWROriF50SkDZBT3ptU9YKDvBSxu76U7C7KGz2GrFEfFw8qd100GU74kXXHn89Py/eQMMIGlU3kGTZsGCNHjuTdd9+1weIQ8ZII7gY+AVqIyGtAX+CyYAYVKfwHlQ87sy9dfCvJve1OjlowmiRwu4vmWjIwtUrJGUMAAwYMYIC7/Vn79u1ZunRpiKOKbF5mDc0SkcVAb5wuoetV9begRxYh9m8l9Cd9+jZOmH03PpR4cmh+68XQcQr063fgvgnGGBMAnvYsBuKBHe75nUQEVZ0XvLAiV+NzTiZn9sPuoHIU9Qt2Qv/+7Dm6J+1+yKAj+eTOjLWBZWNMwHiZPvoocB6wHCh0DytgiSAIDlipfHFXeP554m69i9jioZkcMmek7b+ZjjEeqSrOek5TW1R1p0kvLYKzcNYNlDtAbALjgB3TbryR1THdaDf6NGLII5pCum95HxaeQsbiPOsuMp7Fx8eTmZlJQkKCJYNaQlXJzMwkPj6+0tfwkgjWATF4mClkgqfzdYPIiPucHW98SruGv9Jk3pvQqxednAIW5No6BONB8+bN2bRpE9u3bw93KCaA4uPjad68eaXf76XW0AygK0610OJkoKqjK33XCgr25vU10u7drOsxkqQ1sxCcPrtFKdfRa8FT4Y7MGFNNBKzWEM5q4PeqHpIJqIYN2XfzvWSP+pJYsolC6ZU+AUZuYVXns9jy1c/WXWSM8cTL9NGX3AVlLVX1xxDEZDzyH1g+fEgPOu2YT/5j42g3YwZtgdyZ8WTwmSUDY0yZvHQNDQPGAbGqmiQi3YD7VHV4KAIE6xqqiC8H3UmfuQ8ShaLAtiOP4Yi5b5Axb4cNKhsTYQLZNXQPcBzOJjKo6hIRKXcvAhMejc4/ney5/yKGXEA47PdVaKdOdNQoQG1Q2RhzAC9F5/JK2aKysNQzTdj5F7f7YeI8fJt+4udWJ+CjgGgKiSebfa/MCHeYxphqxEuLYLmIXAj43L0IRgO2oWg1VnIdwq7bHiVx1KDi4nbHL3gSrtjBj21OZ+vnq6y7yJgI5yURXAfcjjN19H/Ap8ADwQzKBJYzqPwZmTPSOOLkLnTcPIfCZ56lff5/aVe8F4INKhsTqcodLK4ObLA48BaceBu95z1aPKi8pflxNFswHVq0CHdoxpgA8TpY7GWMwNRCDS4aTjbx5OOjAB9Nti6Ctm3JHDiSBSfeRsak9HCHaIwJEa/VR00tc0Bxu8FN+f3i0TROm0EKkD9vHCt3vUXHm4eFO1RjTJBZIohgJQeVl9btzQl8QDSFRJNP+1vPhl9Gw9ixZLy7ztYhGFNLeSlD3RyYAJyAU376C5zNaTYFOTYTYgkjBpA7Mw4ll3xiyDp+II3Hj6dgwtMcnVeIoLYXgjG1kJcxgsk4tYaaAE2B991jppbxX4OwZuJnNF7wEaxcyfbDOhFNPtEUEEc2u1//KNyhGmMCyEuJiSWq2q28Y8Fks4bCK2NSOm391iEUxtXBN+ZGViT0Y9sni627yJhqKpAlJjJF5GJgqvv8AiCzKsGZmsV/HULTlJa0X/kePPggHYGjEXJmxlt3kTE1mJcWQSucMYIUnDGCBcB1qvpz8MNzWIug+vm2z/X0SH8Kp4IRbGhzEkmLZkCjRuEOzRjjCuQ6guaqOlxVE1X1cFU9C2hZxeD+ISLLRWSZiEwVkcrvsWbCIv6y88mmDvn4KCSKpLVzICmJX84axRcn3WPrEIypQby0CBaravfyjnm+oUgz4Eugk6pmicibwEeqOuVg77EWQfWUMSn9zymlPePYffloGi6djwJ5xLD6qY/pfN1J4Q7TmIhV5TECEUkB+gCJInKj30sNAV8V44sG6ohIHlAX2FLF65kwKLkOYfGRZ3DC0nSiKSSWPNrdNByy74G//x3q1QtfoMaYMpXVNRQL1Mf5pd3A77EbGFnZG6rqZpyNbn4CtgK7VHVmZa9nqo+EEQPIJY48fGQTR067ZBg7lrwmLVjT4QyWP5MW7hCNMaXwNFisqhsDdkORQ4EZwHnATmAaMF1VXy1xXiqQCtCyZcseGzcGLAQTRPt1F6WmsPaWSbR+7Gp8FKLA1nNvoOlLD0O8DQsZE2wBGywOZBJwnQysV9XtqpoHvIXTBVXyvpNUtaeq9kxMTAxwCCZYklNTGPDpbcVTSX9ekokigDO7qOmb46FNG3jmGcjJCWOkxpgi4ag++hPQW0TqiogAJwErwxCHCQGnuyjW7S6qw7obJziJ4NpryW50OJsvGgu5ueEO05iIVm4iEJG+Xo55papfA9OBxUCGG8Okyl7PVG/+ZSvWTpzDUf+6loyLHyGbWOJydtPsf4+Te2RL+M9/WPbcF6QNftimnhoTYiGfPloZNn20dkkb/DAnzLyTaArIJ4p9jZrSaNcmnLJ2kEM8ayfaSmVjqqrKYwQikiIiN+FOH/V73EPVp4+aCObfXZRLHD899gZLu1+KoPhQ6pBF9IQnID8/3KEaExHKqjVUcvpokSpNHzXmgE1xUlPIQMga9Sax5BCF0nHZdOj0Pdx1Fxl7W5P59hdW3M6YIAn59NHKsK6hyPDn1NMTST58G9xzD3z/vdtlJOQQZ11GxlSA166hgyYCERmvqjeIyPs4M//2o6rDqx6mN5YIIlRhIcu6X0Ln7/+H4PxHuOKY8+n83WsQZdttG1OeQJShfsX9c1xgQjKmgqKi0L9fS9aot4v3Qui89HXouoyNfS9g/TpIGDnQWgjGVFG5XUP7neysCm6hqkuDF9KBrEUQ2Yq7jP7Sj+RGP5Pzj1uJ+/UnFMglllUT55KcesCaRGMiXsA2phGRNGC4e+4iYJuIzFfVG8t8ozEBUrK4XfrkdfSbdRc+Cokjl1ZjzoEmz5OxNYHMGZ/boLIxFeRlh7JGqrpbRK4EXlbVu0UkpC0CY/wljBxEzqwHiSEXxUdcHDB8OJ2LBpVnxtmOacZUgJcRt2gRaQKcC3wQ5HiMKZf/auUfJ6YRt2UDP3Q+212HUEgdsiicOAlUyZiUbquVjSmHlxbBfcCnwJeq+q2IHAWsDm5YxpStZHdR3ugxZI36uHhQueviKexru4h261bRkXxyZ8ZaK8GYgyg3EajqNJxS0UXP1wEjghmUMRXlv0jtsDP70EVW4PvHrcRTVOE0h8wZafslD2OMo6wdysaq6mMiMoHS1xGMDmpkxlTQ/q2EE1lT2JH2155KDHlEU0jPDdPhixPIWBm936pmYyJdWS2CotLQNm/T1EhdrhlARszn7HhzFm0TfqfpvDegf386EeVMPbVBZWMAbyUmuqvq4hDFUypbR2AC4o8/WHPsObRZ9RECFALfHX81Pb569oCd1YypDQK2jgD4l4gcibOHwBuquqzK0RkTDnXrknXTHWSNmksc2USh9Pj6OXYnL6HdssU2qGwilpetKgcCA4HtwEQRyRCRO4IemTFBUDT19ItTH2TFU7Ph4Yep88MS4skhmgJiyXUGlY2JIBUtMZEMjAXOU9XYoEVVgnUNmWBa8dRs2l0/hGjyEWB3t/40nPIUdO0a7tCMqZJAlpjoCJyHM2U0E3gDuKnKERpTTXQafTIZ8fPY+cYntK+7mSO+mA7durHr2AGsij+G+MvOt64iU6t5GSxOB14HpqnqlpBEVYK1CExI7djBtktuIvHDyQiQj49197xM+7svDHdkxlRIlbeqLKKqKar673AlAWNC7tBDWZHXjgJ3R1YfBbS75yK46CL48UcrW2FqHS+zhoyJOAkjBpA7MxYll3xi2HvaOSS++xY6dSodNQpQW4dgag3b5smYUvgXtlsz8TMSP34Z1q/n51Yn4KOAaAqJJ5t9r8wAKLWVYC0HU1NUaNZQuNgYgakuMial03bUoOLidkRFsSPldOrMn00MeeQSy9qJcwBoM+okYsktPmYtBxNqgZw11B64GWjlf76qDqpShMbUQE5xu8/InJHGESd3oePmOTSa8DRRFCAAfusQOpJLNAVo0TFLBKaa8jJGMA14HvgPUBCIm4rIIcALQBecgnZ/VVVrP5saYf/idsNY3XQQbW4ZQTT5+CggOftbtp9+afEYQx6xJIwYEM6QjSmTl0SQr6rPBfi+/wY+UdWRIhIL1A3w9Y0JmaPHDifjkHnsfe0d2rOKhPQPSPjqQ34bOIxVhe1ocOFw6xYy1ZqXRPC+iPwdeBuKi7ujqr9X5oYi0gjoD1zmXicXyK3MtYypLvZrJWzcCA8+yGGTJ3OYzwfdsmBra2jSxIrbmWrJy4Ky9aUcVlU9qlI3FOkGTAJWAF2BRcD1qrrvYO+xwWJTI61fDw88AC+9BDExbD/hLOrPfme/QWVLBiaYArmgLKmUR6WSgCsa6A48p6rHAvuAW0ueJCKpIrJQRBZu3769CrczJkySkuDFF+GHH+Dcc0mY/QbxZBNNATFW3M5UI+UmAhGJEZHRIjLdfVwrIjFVuOcmYJOqfu0+n46TGPajqpNUtaeq9kxMTKzC7YwJs7Zt4aWXWH3fVArwoUA0BXSOWgmZmeGOzhhPYwTPATHAs+7zS9xjV1bmhqr6i4j8LCIdVPVH4CScbiJjarUOd55HxhEtyXrpTdrnLiPx01eh9dts6z+SVVktaHT+EOtAsQ7LAAAaTklEQVQqMmHhZYzge1XtWt6xCt3UGSd4AYgF1gGXq+qOg51vYwSmVlq+nJ2XjuaQRZ+hQD7RrH7yQzrdcGq4IzO1RCB3KCsQkTaquta98FFUcT2Bqi4Byg3OmFqtc2eWJJzMCaQRTSEx5NPu5jNhzz9Z0SiFbR9+a7OLTEh4SQQ3A3NFZB0gOCuMLw9qVMZECKe4XRxKLgVEk9u5FzF33UVHoANC7sx4K2xngs7LrKE5QDtgNHAd0EFV5wY7MGMigX9xu9UT59JwyTwW9r4GBXwo8WRRd9y9sHcvYIXsTHB4GSPwAWcArdm/1tATQY3Mj40RmEiSMSm9uGCdoPgohMMOY+uJ59FoxovE2joE41EgxwjeB7KBDKCwqoEZY8rmFLab8+cK5GTgnntoMuMZFKd/1grZmUDykgiaq+oxQY/EGFNs/8J2wKefsvaWibR+7O/4KCSaAtofsg2ys8l4+TsrW2GqxEsi+FhETlXVmUGPxhhzUG0eHUVGm2PIf2Ey7XYtpOmb48n79FU67NpFFIXkzoy1gWVTKV52KPsKeFtEskRkt4jsEZHdwQ7MGHOg5NQUjv1mEvV/XAyffcY+qU8seURTQBw5/D5tTrhDNDWQl0TwBJAC1FXVhqraQFUbBjkuY0x5Bg7k50dfI5s4CgEfhfRe/Ay88ALk5YU7OlODeEkEPwPLtCbsaWlMhElO7cPqiXOZd+qDrB/9BHHtWsHf/kZO01b80GUEy577ItwhmhrAy/TRKcBRwMfsvx+BTR81prpRZcPoJ2j59M1EoRQibLnsdpr/524y/vutDSpHmEBOH13vPmLdhzGmuhJhw6pcmhPl7qOsNJ/yADkfTaH9tl/x2aCyKUW5iUBV7wUQkbqq+kfwQzLGVIVTtuLP/ZK3XXUXCa/+mwY44wZCDpnTP7M1CKZYuYlARFKAF4H6QEsR6QqMUtW/Bzs4Y0zFHbAgLTWFjG79aH/VIGLJdQaVl70Ab7YjY0czMt+aZ91FEc7LGMHXwEjgPXdHMURkmap2CUF8gI0RGBMIGZPSyZw+l6NaFdBy/lRYuZJCBEXIIc5KVtRCAduqEkBVfy5xqEplqI0xoZecmsKAmf+k5X/uhIwMViSfV1zLqA5Z8NyzoEpG+7NJa/wXK2wXQTxNHxWRPoC621aOAVYGOS5jTDD5fBRcez1Z1KGAKBQhecmr/NH6aNqv/pATdrxPm1EnWTKIEF4SwVXANUAzYDPQzX1ujKnBikpgf3HqAyx/9nOYMgXZvp04commgFhyyJxhFecjgZdZQ78BF4UgFmNMiO1f3K4fa/84ivZ/P4kY8oimkF6rX4OZPeGUU0AkrLGa4DloIhCRCcBBR5JVdXRQIjLGhE2Xq/uR4fuc36fNoc0Re2g+byoMHgx9+sB998GgQZYQaqGyuoYWAouAeKA7sNp9dMMWlhlTayWnpnDirDto/uqjsHo1PPssbNwIJ5/MzthE1t30TLhDNAF20ESgqi+p6kvAMcAAVZ2gqhOAk3CSgTGmtouLg6uvZtltr5JLDI3yMznqiWvZ26E7fPGFzTCqJbyUmDgUaAj87j6v7x4zxkSI395L52gKEaCAKGI2rIb+/elEFArkjvrUylbUYF5mDT0CfCciU0TkJWAx8FBwwzLGVCcJIwaQSyx5+MghjjX/eo817YcQRSHRFBJPFrn/fSXcYZpK8jJraLKIfAwc7x66RVV/CW5YxpjqpNSyFbHxZI2aSxw5RKH0+Po5OGMjazL+YNPeQ0h4ZKy1EGqIcktMBO3GIj6cAenNqjq0rHOtxIQx1VPGpHQyZ6SReEYvOu/7lvwHHib6jz0UArnEsXriXEsGYRTIMtTBcj3OCmXb7cyYGmr/dQgnkz47mz6f3Y8PJZ4cjrzzSug9lYyR95D5m1oroZo66BiBiCQF66Yi0hw4A3ghWPcwxoTeIeedRg7x5OEjj2gO3bMRunal4+r3OGHHe/uXrRgwwHn4sVlI4VHWYPF0ABEJxm7Y44GxQGEQrm2MCZOishXzT72fHybOI3rzz2w4aiA+CtxB5WyyXnoTgIwtjUlbemjxL/2MSem0Wf2J1TkKg7K6hqJE5J9AexG5seSLld2qUkSGAttUdZGIDCjjvFQgFaBly5aVuZUxJgz27y6CPbc8SPaoQcS6g8q9Foxnx/HLabt6Hh3JL556mjkjjY5unSMll8wZabZ5ToiUlQjOB85yz2kQwHv2BYaLyOk4q5Ybisirqnqx/0mqOgmYBM5gcQDvb4wJIWfG0Wdkzkjj8MHd6LTtcxo8/gQ+8nCKVeQWz0by31ktYcSA8AYeQbxsTDNEVT8Oys2dFsEYmzVkTGRZMe5D2t18FtHkA7AjZQiNX5lAxpCbbVA5gAI5a2iBiDwB9Heffw7cp6q7qhKgMSZydRpzBhkN57Hnf+/TLnodifPfhQ4dSL70UrjjDkgqMVelaFA5LS3UoUYELyuL/wvsAc51H7uByYG4uaqmldcaMMbUTsmpKfRJe4jE2a/DunVwzTXw2mvQvj2kpjqF7lwlB5ZNYHnpGlqiqt3KOxZM1jVkTITYvBkeegheeAFU4cor+aHVYFreegGx5JJLrO2tXAGB3LM4S0RO8LtwXyCrKsEZY0ypmjWDZ56BNWvgiivghRdod9sI4skmmgJiimYTmYDyulXlMyKyQUQ2AE8Do4IalTEmsrVoAc89B6tXs6vvGQjq7pIlHHFKcpiDq3081xoSkYYAqro7qBGVwrqGjIlsP7YeTMOtP3Jk/s+Iu0cCY8fCEUeEO7RqLZBdQ4CTAMKRBIwxpsOGT2mSswH58Uc491wYPx6OOspJBtu3hzu8Gs9zIjDGmLBr2xamTIGVK+Evf4Fx45ypprfdBpmZpb+nlJpGZn9lJgIRiRKRPqEKxhhjPGnfHjZtgh49YNgwePRRaN0aWrWCvn33O9WmnpavzESgqoWA7VRtjKme6tWDqVMhIwOGDIGffoKvv4Z77oGdO62QnUdeuobmiMgIEZGgR2OMMR7t902/c2d4801WtxzEtqjD4d57ISmJ+o/dRaxbyK7cqacR3IXkpcTEKOBGoEBEsgABVFVtQxljTFgUfdPvSG5x9VKANj+lE0su2cSR2yqZpO9no0ABQr4VsjuoclsEqtpAVaNUNUZVG7rPLQkYY8Imc0baAd/0/Y/5yGfxEUNg4UL21D0SH0pMvRiSd8yDvXtLvWYkjyWUmwjEcbGI3Ok+byEixwU/NGOMKV3CiAHkEuvuhOZ80y/tGD160HDfVvj6a6L794Vbb3WmnbZpA/36FV8v0scSvIwRPAukABe6z/diA8jGmDBKTk1hbbvTmH/osOLaQ6UdK3bccfDRR7BgAXTr5hS5+/prePJJyMoqtYURSbwUnVusqt1F5DtVPdY99r2qdg1JhNjKYmNMYK1tfiKNf13Bofm/wZFHsqX/+Rzy5kRi3E1xakthu0CuLM4TER84pT5EJBHba9gYU0NlTEqnyeZvaZC/g2zi2NfgSJq+OZ4Yn7K+bmfWPfNxrUgCFeElETwFvA0cISIPAl8CDwU1KmOMCZKSg8rfJp0Dc+YQUy+W9n8spcsjl8DEiZCbG+5QQ8bLrKHXgLE4v/y3AGep6rRgB2aMMcFw4KDyQBg0CHbuhJkzoXlzuOoqZ/XyCy9AXl64Qw46r7WG6gI+9/w6wQvHGGOC66CDyiJwyikwfz588olT2fRvf4MOHWDyZMjPD2/gQeRl+uhdwEtAY+AwYLKI3BHswIwxJliSm/7OgGN2lD4WIAKDB8NXX8GHH0LjxvDXv8LRR8PLL9fKhOClRXAR0EtV71HVu4HewCXBDcsYY4IoLc15lEUEHnvMqWf03nvQsCFceil06uTsrVxQEIpIQ8JLItgCxPs9jwM2ByccY4ypZkScCqeLFsHbb0OdOnDxxdClC7z+eq1ICAdNBCIyQUSeAnYBy0VkiohMBpYBO0MVoDHGVAsicNZZ8N13MG0a+HxwwQVwzDHO88KaO6u+rBbBQmARztTRfwJzgTTgduDdoEdmjDFhVmr9oagoGDkSli51WgSqzq5p3brBW2/VyITgec/icLKVxcaYUMuYlE6bUScRSy65Za02LiiAN95wSl+vWgVduzo/Dx/utCLCKGAri0VkqIh8JyK/i8huEdkjIpXeu9gtWjdXRFaIyHIRub6y1zLGmGAps/6Q/94FPh9ceCEsX+7MKtq3z+lC6tkTPvjAaTFUc14Gi8cDlwIJASpDnQ/cpKqdcGYgXSMinapwPWOMCbhSq5mWJToaLrnE2U958mRngdqwYXD88fDxx9U6IXhJBD8DyzRAfUiqulVVF7s/7wFWAs0CcW1jjAmUMquZliU6Gi67DH74wVmZvH07nH469OnjrFyuhgnBSyIYC3wkIreJyI1Fj0DcXERaA8cCX5fyWqqILBSRhdu3bw/E7YwxpkLKXHhWnpgYuOIK+PFHp3bR5s3OQrV+/WDOnGqVELwkggeBP3DWEjTwe1SJiNQHZgA3qOoBYw6qOklVe6pqz8TExKrezhhjwiM2FlJTYfVqePZZ2LABTj7ZGWP4/PNwRwd4SwRNVfVsVb1bVe8telTlpiISg5MEXlPVt6pyLWOMCbVKbWsZFwdXXw1r1sCECU5iGDDAKXj35ZdBi9ULL4ngIxE5NVA3FBEBXgRWquoTgbquMcYEXCmlKKq8rWV8PFx7Laxd6+yQtmKF0110yimQHp4tMr0kgquBT0QkKxDTR4G+OLWKBonIEvdxehWuZ4wxIROwbS3r1IEbbnC2zRw3Dr7/3hlQHjIEvvkmoDGXx8t+BA1UNUpV6wRi+qiqfqmqoqrHqGo39/FRZa9njDGhVOFppeWpWxduugnWr4dHH4Vvv3WmnA4d6tQ3CgEvexb3L+24qs4LSkSlsJXFxpjqJKP92WT+piQ8MvbPGUVFC8zKq2pa3nl79sDTT8Nddzklr7/6ykkMleB1ZXG0h2vd7PdzPHAcTg2iQZWKzBhjarjkpr9DUyAYexs3aAC33ebshbBtGxx3XODvUUK5iUBVh/k/F5EWOKuNjTHGBEt0NDRtGpJ6RV63qvS3CegY6ECMMcaER7ktAhGZABQNJEQB3YDFwQzKGGNqBa/jBmHmZYzAf5Q2H5iqqvODFI8xxlR/1fwXe0V5GSN4KRSBGGOMCQ8v+xH0FZFZIrJKRNaJyHoRWReK4IwxpqbwWnYi0OcFgpeuoReBf+BMGa35uzQbY0yAFZWd6EguuaM+JYPSy1YH+rxA8TJraJeqfqyq21Q1s+gRtIiMMaaG8Vp2ItDnBYqXRDBXRB4XkRQR6V70CGpUxhhTg3gtOxHo8wLFS4mJuaUcVlUN2cpiKzFhjKnuvJadqNJ5FRSwEhOqOrBSERhjTATxWnYi0OcFwkG7hkTkYhEp6/U2InJCcMIyxhgTKmW1CBKA70RkEc6Moe04RefaAicCvwG3Bj1CY4ypoTK2NHa6dyalB3XWT1UdNBGo6r9F5GmcKqN9gWOALGAlcImq/hSaEI0xpuYJ9RTQqihzjEBVC4BZ7sMYY4xHmTPS6OhOAdWiKaDVNBFUpvqoMcaYcoR6CmhVWCIwxpggSE5NYW2705h/6DDWTqy+3ULgrcSEMcaY8pRSkbTUKaBeK5eGsMKpl/0IegL9cD5OFrAMmKWqO4IcmzHGmBAoa53A5SKyGLgNqAP8CGwDTgBmi8hLItIyNGEaY4wJlrJaBHWBvqqaVdqLItINaAfYNFJjjKnBylpH8ExZb1TVJYEPxxhjTKh5GSNIBP4GtPY/X1X/WtmbishpwL8BH/CCqj5S2WsZY4ypGi+zht4FvgBmE4CNaUTEBzwDnAJsAr4VkfdUdUVVr22MMTVRuEtReEkEdVX1lgDe8zhgjaquAxCR14EzAUsExpiIUx1KUXhZUPaBiJwewHs2A372e77JPbYfEUkVkYUisnD79u0BvL0xxlQfod6NrDReEsH1OMkgS0R2i8geEdkd7MBUdZKq9lTVnomJicG+nTHGBF5aWrkLw6pDKQovG9M0CPA9NwMt/J43d48ZY0zESU5NIWPcaVXejawqDpoIRKS1qm4o43UBmqnqpgre81ugnYgk4SSA84ELK3gNY4ypNUK5G1lpymoRPO7uUPYuB25MMxA4Cbgbp4/fM1XNF5FrgU9xpo/+V1WXVyJ2Y4wxAVDWgrJzRKQTcBHwV6AJf25M8yHwoKpmV+amqvoR8FFl3muMMSawytuYZgVwe4hiMcYYEwZeVhafXcrhXUCGqm4LfEjGGGNCycuCsiuAFGCu+3wAzphBkojcp6qvBCk2Y4wxIeAlEUQDHVX1VwAROQJ4GTgemAdYIjDGmBrMy4KyFkVJwLXNPfY7kBecsIwxxoSKlxZBmoh8AExzn490j9UDdgYtMmOMMSHhJRFcA5yNszMZwEvADFVVnPUExhhjajAvJSZURL4EcgEFvnGTgDHGmFqg3DECETkX+AanS+hc4GsRGRnswIwxxoSGl66h24FeRWsG3B3LZgPTgxmYMcaY0PAyayiqxMKxTI/vM8YYUwN4aRF8IiKfAlPd5+dhdYKMMSZwytmzINi8DBbfLCIjgL7uoUmq+nZwwzLGGBMqXloEqOoMYEaQYzHGGBMGZW1MswdnuugBL+HMKm0YtKiMMcaETFn7EQR6i0pjjDHVkM3+McaYCGeJwBhjIpwlAmOMiXCWCIwxJsJZIjDGmAgnNaGQqIhsBzZW4C2HAb8FKZyaIJI/v332yGSfvXStVDWxvAvUiERQUSKyUFV7hjuOcInkz2+f3T57pAnEZ7euIWOMiXCWCIwxJsLV1kQwKdwBhFkkf3777JHJPnsV1MoxAmOMMd7V1haBMcYYjywRGGNMhKt1iUBEThORH0VkjYjcGu54QklE/isi20RkWbhjCSURaSEic0VkhYgsF5Hrwx1TqIhIvIh8IyLfu5/93nDHFGoi4hOR70Tkg3DHEmoiskFEMkRkiYgsrPR1atMYgYj4gFXAKcAm4FvgAlVdEdbAQkRE+gN7gZdVtUu44wkVEWkCNFHVxSLSAFgEnBUJ/+4iIkA9Vd0rIjHAl8D1qvpVmEMLGRG5EegJNFTVoeGOJ5REZAPQU1WrtJiutrUIjgPWqOo6Vc0FXgfODHNMIaOq84Dfwx1HqKnqVlVd7P68B1gJNAtvVKGhjr3u0xj3UXu+3ZVDRJoDZwAvhDuWmqy2JYJmwM9+zzcRIb8QjENEWgPHAl+HN5LQcbtGlgDbgFmqGjGfHRgPjAUKwx1ImCgwU0QWiUhqZS9S2xKBiWAiUh9nb+0bVHV3uOMJFVUtUNVuQHPgOBGJiG5BERkKbFPVReGOJYxOUNXuwBDgGrd7uMJqWyLYDLTwe97cPWZqObd/fAbwmqq+Fe54wkFVdwJzgdPCHUuI9AWGu/3krwODROTV8IYUWqq62f1zG/A2Tvd4hdW2RPAt0E5EkkQkFjgfeC/MMZkgcwdMXwRWquoT4Y4nlEQkUUQOcX+ugzNR4ofwRhUaqnqbqjZX1dY4/69/pqoXhzmskBGReu7kCESkHnAqUKkZg7UqEahqPnAt8CnOgOGbqro8vFGFjohMBdKBDiKySUSuCHdMIdIXuATnG+ES93F6uIMKkSbAXBFZivNFaJaqRtw0ygh1BPCliHwPfAN8qKqfVOZCtWr6qDHGmIqrVS0CY4wxFWeJwBhjIpwlAmOMiXCWCIwxJsJZIjDGmCAIZBFIEekmIuluYcGlInKe32tJIvK1W2jzDXfqfIVYIjC1nogU+E0rXVKdqtKKyHQROcr9eYOIHCYisSIyT0Siwx2fqZIpBG5x3x/A/6lqZ/ea44vWjwCPAk+qaltgB1DhaeOWCEwkyFLVbn6PR6p6wUD8khaRzoBPVdf5H3cLJs4Bziv1jaZGKK0IpIi0EZFP3NpAX4jI0R6vtUpVV7s/b8GpK5XoLqYcBEx3T30JOKuisVoiMBHL/QZ+r4gsdmu6H+0er+c2679x69yf6R6/TETeE5HPgDkiEiUiz4rIDyIyS0Q+EpGRIjJIRN7xu88pIvJ2KSFcBLx7kPDecV83tcsk4DpV7QGMAZ6t6AVE5DggFlgLJAA73cW0UMlCm9b0NJGgjluds8jDqvqG+/NvqtpdRP6O8z/mlcDtOOUK/uo2v78Rkdnu+d2BY1T1dxEZCbQGOgGH46xm/y9OvZ9nRSRRVbcDl7vHS+oLTD1IzMuAXpX8vKYacosi9gGmOV/kAYhzXzsbuK+Ut21W1cF+12gCvAJcqqqFftepEksEJhJkudU5S1NUoG4RcLb786k4xczGuM/jgZbuz7NUtai5fwIwTVULgV9EZC44ewSIyCvAxSIyGUgB/q+UezcBtpcWlKoWiEiuiDRw91gwNV8Uzrf3A/5bdAslllksUUQaAh8Ct/ttPJQJHCIi0W6roFKFNq1ryES6HPfPAv78YiTACL8xhZaqutJ9bZ/H604GLgYuwEkW+aWck4WTZA4mDsj2eD9Tzbml0deLyDngFEsUka5e3uvOBHobZ/fBovEA1KkRNBcY6R66lIN3Nx6UJQJjDvQpcJ07EIeIHHuQ8+YDI9yxgiOAAUUvuAN6W4A7cJJCaVYCbUt7QUQScLqt8ir1CUzYHaQI5EXAFW6huOV430HxXKA/cJnf7LeilsUtwI0isgZnzODFisZqXUMmEpQcI/hEVcuaQno/zs5XS0UkClgPlLYX7gzgJGAFzs54i4Fdfq+/BiT6tSZK+hAneRSNP0TzZwtloPu6qaFU9YKDvFThKaWq+ipQ6l4L7qyzSu1DUMSqjxpTBSJS3904PgGnFHBfVf3Ffe1p4DtVLfUbmrt/wFycQePGwBJVbea+9hZwq6quCsXnMJHNWgTGVM0H7syiWOB+vySwCGc84aaDvVFVs0Tkbpzugn8Ct7nvjQXesSRgQsVaBMYYE+FssNgYYyKcJQJjjIlwlgiMMSbCWSIwxpgIZ4nAGGMi3P8DeRQ+bNPiuDMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#\n",
    "# Use a \"curve fit\" to find the temperature of the demon\n",
    "#\n",
    "\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def fLinear(x, m, b):\n",
    "    return m*x + b\n",
    "\n",
    "energies = (bins[:-1]+bins[1:])/2.0\n",
    "xvals = np.array(energies)  # fit log(n) vs. energy\n",
    "yvals = np.log(np.array(ns))\n",
    "sig = 1.0/np.sqrt(np.array(ns))\n",
    "\n",
    "#\n",
    "# make initial estimates of slope and intercept.\n",
    "#\n",
    "\n",
    "m0 = (yvals[-1]-yvals[0])/(xvals[-1]-xvals[0])\n",
    "b0 = yvals[0]-m0*xvals[0]\n",
    "\n",
    "popt, pcov = curve_fit(fLinear, xvals, yvals, p0=(m0, b0), sigma=sig)\n",
    "\n",
    "m=popt[0]          # slope\n",
    "dm=np.sqrt(pcov[0,0]) # sqrt(variance(slope))\n",
    "b=popt[1]          # int\n",
    "db=np.sqrt(pcov[1,1]) # sqrt(variance(int))\n",
    "Temp=-1.0/(m*kb)   # temperature\n",
    "dT = abs(dm*Temp/m)# approx uncertainty in temp\n",
    "\n",
    "print(\"slope=\", m, \"+/-\", dm )\n",
    "print(\"intercept=\", b, \"+/-\", db)\n",
    "print(\"Temperature=\", Temp, \"+/-\", dT, \"K\")\n",
    "pl.title(\"Demon Energy Distribution\")\n",
    "pl.xlabel(\"Energy (J)\")\n",
    "pl.ylabel(\"log(n) (number of demon visit to energy)\")\n",
    "pl.errorbar(xvals, yvals, sig, fmt='r.')\n",
    "pl.plot(xvals,yvals,'b.',label=\"Demon Energies\")\n",
    "pl.plot(xvals,fLinear(xvals, m, b),'r-', label=\"Fit\")\n",
    "pl.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Proj 5: (option 1) \n",
    "==========\n",
    "\n",
    "Velocity distribution in a 1-D gas\n",
    "--------------------------------\n",
    "\n",
    "For project 5 use the above example as a starting point to investigate the distribution of velocities of nitrogen molecules in a 1D gas. Start the velocities out in some random distribution, then each 'step' should bump one molecule up or down in velocity randomly. Use the demon algorithm to exchange energy between the system and the demon. Use the same strategy to show a final distribution of velocities and energies to determine the  temperature of the system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Proj 5: (option 2) \n",
    "==========\n",
    "\n",
    "Monte-Carlo Simulation\n",
    "---------------------\n",
    "\n",
    "The main point of this project is to use the generation of random data to answer non-trivial questions. If you don't love the idea of the velocity distribution of particles in a gas, you can attack any problem you care to dream up that could be modeled using a large number of random numbers in some way. If you choose to do this, please let me know your plans!\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}