Simple-diffusion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Why bother studying diffusion?\n",
    "# Simplest example of stochastic behavior.  ('Stochastic' = of random character.)\n",
    "# Random behavior and random numbers very important in science and computing - e.g., molecular dynamics, cryptography\n",
    "# Diffusion is key reference point for understanding more complex scientific phenomena."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the required modules\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# This makes the plots appear inside the notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7623467632724313 0.6813560666651842\n"
     ]
    }
   ],
   "source": [
    "## (Pseudo) random numbers\n",
    "# So where are 'real' random numbers?\n",
    "#seed = 12345  # a nmumber that tells the random number generator (\"rng\") where/how to start its sequence of \"random\" numbers\n",
    "#rng = np.random.default_rng(seed) # shortcut for function\n",
    "# The sequence is not truly random.  If we use the same seed, we will get the same sequence of numbers.\n",
    "rng = np.random.default_rng() # shortcut for function\n",
    "### Toggle with and without seed\n",
    "rfloat1 = rng.random() # a single random number \n",
    "rfloat2 = rng.random() # another \n",
    "print(rfloat1,rfloat2) # display the two numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the range function\n",
      "[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]\n",
      "\n",
      "trial loop\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "\n",
      "array test\n",
      "0.9090712819683793\n",
      "0.8887009172592641\n",
      "0.08216118357772395\n",
      "0.9021811151551513\n",
      "0.9732757372879874\n",
      "0.7442106217979246\n",
      "0.9849376131386149\n",
      "0.5147187078367013\n",
      "0.3697614324896277\n",
      "0.14845250493789042\n"
     ]
    }
   ],
   "source": [
    "# Looping to create a list of random numbers\n",
    "print(\"The range function generates this array:\")\n",
    "print(list(range(10)))\n",
    "\n",
    "print (\"\\nTrial loop\")   # \\n gives new line\n",
    "for i in range(10): # remind ourselves how loops work\n",
    "    print(i)\n",
    " \n",
    "print (\"\") # skip a line (another way, more clumsy than \\n)\n",
    "print (\"array test\")\n",
    "rng_array = []\n",
    "for i in range(10): # remind ourselves how loops work\n",
    "    rng_array.append(rng.random())\n",
    "\n",
    "for x in rng_array:\n",
    "    print(x)    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.08899193769704539\n",
      "-0.10682945579379166\n",
      "-0.6302516771588726\n",
      "-0.0453875403234536\n",
      "-0.14528095325246482\n",
      "-0.9584434295367956\n",
      "0.8689655725938519\n",
      "0.2731462401173592\n",
      "-0.9691924813803792\n",
      "-0.12443534436474613\n",
      "0.889299802141631\n",
      "-0.5003066518038239\n",
      "0.06761099670768567\n",
      "0.37556578432602983\n",
      "-0.3931379837199791\n",
      "-0.6028918390485678\n",
      "-0.02127591310436716\n",
      "-0.14815896474564916\n",
      "-0.565601863653888\n",
      "-0.15212494564861623\n"
     ]
    }
   ],
   "source": [
    "# Transform a (0,1) random number to a (-1,1) random number\n",
    "rng_array = [] # create an empty array (list) which we'll fill with (-1,1) random numbers\n",
    "\n",
    "for i in range(20):\n",
    "    rng_array.append(2.0*rng.random()-1.0) # need to shift and expand range (or expand, then shift)\n",
    "\n",
    "for dx in rng_array:\n",
    "    print(dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Number of occurrences in array')"
      ]
     },
     "execution_count": 184,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYCUlEQVR4nO3de5gldX3n8fcHBAHlGkZFYRw16K5XMI1BMRrRNQoCivdNFC+7oxsl+BgvY3AFg7tRvGRNdHVHxYXEBQWBFfCGKBBXRWYQEASj4JgQUAbDTQwgzHf/ONVyZuhL9TmnTvcc36/nOU9X1amq36ere75T/TtVv0pVIUmaPFssdgBJUjcs8JI0oSzwkjShLPCSNKEs8JI0oe6z2AH67brrrrVixYrFjiFJm421a9feUFXLZnpvSRX4FStWsGbNmsWOIUmbjSQ/ne09u2gkaUJZ4CVpQlngJWlCWeAlaUJZ4CVpQlngJWlCdVrgk+yU5JQkVya5IsmTu2xPknSPrq+D/zDw5ap6UZKtge06bk+S1OiswCfZAXga8CqAqroTuLOr9iRJG+vyDP7hwHrg00meAKwFjqiq2/pXSrISWAmwfPnyDuNoEqxYddaitLvuvQcuSrvSMLrsg78P8ETgY1W1N3AbsGrTlapqdVVNVdXUsmUzDqcgSRpAlwX+GuCaqrqgmT+FXsGXJI1BZwW+qn4G/HOSRzWLngn8oKv2JEkb6/oqmsOBzzRX0FwNvLrj9iRJjU4LfFVdDEx12YYkaWbeySpJE8oCL0kTygIvSRPKAi9JE8oCL0kTygIvSRPKAi9JE8oCL0kTygIvSRPKAi9JE8oCL0kTygIvSRPKAi9JE8oCL0kTygIvSRNq3gKfZE2SNyTZeRyBJEmj0eYM/mXAg4ELk5yU5I+SpONckqQhzVvgq+rHVXUk8Ejg/wDHAf+U5N1Jduk6oCRpMK364JM8Hvgg8H7g88CLgFuAr3cXTZI0jHmfyZpkLXAT8ClgVVXd0bx1QZL9ugwnSRpcm4duv7iqrp7pjao6dMR5JEkjMm+Br6qrkxwIPAbYpm/5X3YZTJI0nDaXSX4ceClwOBDgxcBDO84lSRpSmw9Zn1JVrwRurKp3A08G9ug2liRpWG0K/O3N118leTDwa+Bh3UWSJI1Cmw9Zz0iyE71LJC8CCvhEm50nWQfcCtwN3FVVUwPmlCQt0JwFPskWwDlVdRPw+SRnAttU1c0LaOMZVXXDMCElSQs3ZxdNVW2gd4PT9PwdCyzukqRF0qYP/qtJXjjg+DPVbL82ycqZVkiyshnQbM369esHaEKSNJM2ffBvBu4H3JXkdnqXSlZV7dBi2/2q6tokDwDOTnJlVZ3fv0JVrQZWA0xNTdXC4kuSZjPnGXzTB/+cqtqiqrauqh2qavuWxZ2qurb5ej1wGvCkoRNLklpp0wf/gUF2nOR+SbafngaeDVw2yL4kSQvXZR/8A4FvJrkE+C5wVlV9ecEJJUkD6awPvhmg7AnDR5QkDaLNYGPbjyOIJGm02pzB0zyPdU82Hk3y/Nm3kCQttjYP/PhPwBHA7sDFwL7At4H9u40mSRpGmw9ZjwD2AX5aVc8A9ga8I0mSlrhWo0lW1e0ASe5bVVcCj+o2liRpWG364K9pRpM8nd7dqDcC13YbS5I0rDZX0bygmTw6yTeAHQGvZ5ekJa7VVTTTquq8roJIkkarTR+8JGkzZIGXpAllgZekCTVvgU9yaJIfJbk5yS1Jbk1yyzjCSZIG1+ZD1mOBg6rqiq7DSJJGp00Xzc8t7pK0+WlzBr8myWfp3eh0x/TCqjq1s1SSpKG1KfA7AL+i90SmaQVY4CVpCWtzJ+urxxFEkjRasxb4JG+rqmOT/C29M/aNVNWfdZpMkjSUuc7gpz9YXTOOIJKk0Zq1wFfVGc3X48cXR5I0Kt7JKkkTygIvSRPKAi9JE6rNWDTHJtkhyVZJzklyQ5I/GUc4SdLg2pzBP7uqbgGeB1wDPBJ4a6epJElDa1Pgt2q+HgCcWFX/2mEeSdKItCnwZyS5EpgCzkmyDLi9bQNJtkzyvSRnDhpSkrRw8xb4qloFPBmYqqpf0xuX5pAFtHEE99w0JUkakzYfsm4HvAH4WLPowfTO5ueVZHfgQOCTgwaUJA2mTRfNp4E7gac089cA72m5//8BvA3YMNsKSVYmWZNkzfr161vuVpI0nzYF/hFVdSzwa4Cq+jcg822U5HnA9VW1dq71qmp1VU1V1dSyZcvaZJYktdCmwN+ZZFuaESWTPIK+B3/MYT/g4CTrgJOA/ZP8/aBBJUkL06bAHwV8GdgjyWeAc+h1u8ypqt5RVbtX1QrgZcDXq8obpCRpTNo88OPsJBcB+9Lrmjmiqm7oPJkkaShtrqJ5AXBXVZ1VVWcCdyV5/kIaqapzq+p5g4aUJC1cqy6aqrp5eqaqbqLXbSNJWsLaFPiZ1mnzsG5J0iJqU+DXJPlQkkckeXiSvwbmvPRRkrT42hT4w+nd6PRZ4GR649C8octQkqThtbmK5jZg1RiySJJGaN4Cn+SRwFuAFf3rV9X+3cWSJA2rzYelJwMfpzdg2N3dxpEkjUqbAn9XVX1s/tUkSUtJ2wd+/GmS3ZLsMv3qPJkkaShtzuAPa772P4e1gIePPo4kaVTaXEXzsHEEkSSNVqsnOiV5Z5LVzfyezVjvkqQlrOsnOkmSFklnT3SSJC2uLp/oJElaRG2uotn0iU77Aa/qMpQkaXhzFvgkWwA7A4fiE50kabMyZ4Gvqg1J3lhVnwPOGlMmSdIItOmDPzvJW5Ls4Z2skrT5aNMH/5rma/8Y8N7JKklLXJs++FVV9dkx5ZEkjcicXTRVtQGf3iRJmyX74CVpQtkHL0kTytEkJWlCtXkm6ytnWl5VJ8yz3TbA+cB9m3ZOqaqjBgkpSVq4Nl00+/RNbwM8E7gImLPA0xuvZv+q+mWSrYBvJvlSVX1nsKiSpIVo00VzeP98kh2Bv2uxXQG/bGa3al41QEZJ0gDanMFv6lfAnm1WTLIlsBb4XeCjVXXBDOusBFYCLF++fIA4UvdWrFqckTrWvffARWlXk6FNH/wZ3HPmvQXwaOBzbXZeVXcDeyXZCTgtyWOr6rJN1lkNrAaYmpryDF+SRqTNGfwH+qbvAn5aVdcspJGquinJucBzgMvmWV2SNAJtCvw/AddV1e0ASbZNsqKq1s21UZJlwK+b4r4t8CzgfcMGliS10+ZO1pOBDX3zdzfL5rMb8I0klwIXAmdX1ZkLjyhJGkSbM/j7VNWd0zNVdWeSrefbqKouBfYeJpwkaXBtzuDXJzl4eibJIYBPdJKkJa7NGfzrgc8k+Ugzfw0w492tkqSlo82NTlcB+ya5P5CqurX7WJKkYc3bRZPkvyfZqap+WVW3Jtk5yXvGEU6SNLg2ffDPraqbpmeq6kbggO4iSZJGoU2B3zLJfadnmmva7zvH+pKkJaDNh6x/D5yT5NP0hix4DXB8p6kkSUNr8yHrsc3NSs9qFh1TVV/pNpYkaVhtR5P8HvcM9/u97uJIkkalzVU0LwG+C7wIeAlwQZIXdR1MkjScNmfwRwL7VNX18JtBxL4GnNJlMEnScNpcRbPFdHFv/KLldpKkRdTmDP7LSb4CnNjMvxT4YneRJEmj0OYqmrcmORR4KhBgdVWd1nkySdJQWl1FU1WnAqd2nEWSNEL2pUvShLLAS9KEmrXAJzmn+epzVCVpMzRXH/xuSZ4OHJzkJHofsP5GVV3UaTJJ0lDmKvDvAlYBuwMf2uS9AvbvKpQkaXizFviqOgU4Jcl/rapjxphJkjQCba6DP6Z56PbTmkXnVtWZ3caSJA2rzWBjfwUcAfygeR3RLJMkLWFtbnQ6ENirqjYAJDme3pDB7+gymCRpOG2vg9+pb3rHLoJIkkarzRn8XwHfS/INepdKPg3P3iVpyWvzIeuJSc4F9qFX4N9eVT+bb7skewAnAA8CNtAbpOzDw8WVJLXVdrCx64AvLHDfdwF/XlUXJdkeWJvk7Kr6wUJDSpIWrrOxaKrquum7XavqVuAK4CFdtSdJ2ljbh24PJckKYG/gghneWwmsBFi+fPnAbaxYddbA22ph1r33wMWOIHVisepIV/+m5jyDT7JFksuGaSDJ/YHPA2+qqls2fb+qVlfVVFVNLVu2bJimJEl95izwzbXvlyQZ6NQ6yVb0ivtnmoeGSJLGpE0XzW7A5Um+C9w2vbCqDp5royQBPgVcUVWbDlYmSepYmwL/7gH3vR/wCuD7SS5ulv1FVfnAbkkagzbXwZ+X5KHAnlX1tSTbAVu22O6bbDKGvCRpfNoMNvafgVOA/9UseghwepehJEnDa3Md/BvodbfcAlBVPwIe0GUoSdLw2hT4O6rqzumZJPeh90QnSdIS1qbAn5fkL4Btk/wH4GTgjG5jSZKG1abArwLWA98HXgd8EXhnl6EkScNrcxXNhuYhHxfQ65r5YVXZRSNJS9y8BT7JgcDHgavoXfb4sCSvq6ovdR1OkjS4Njc6fRB4RlX9GCDJI4CzAAu8JC1hbfrgr58u7o2rges7yiNJGpFZz+CTHNpMXp7ki8Dn6PXBvxi4cAzZJElDmKuL5qC+6Z8DT2+m1wM7d5ZIkjQSsxb4qnr1OINIkkarzVU0DwMOB1b0rz/fcMGSpMXV5iqa0+mN634GsKHbOJKkUWlT4G+vqr/pPIkkaaTaFPgPJzkK+Cpwx/TCqrqos1SSpKG1KfCPo/dkpv25p4ummnlJ0hLVpsC/AHh4/5DBkqSlr82drJcAO3UdRJI0Wm3O4B8IXJnkQjbug/cySUlawtoU+KM6TyFJGrk248GfN44gkqTRanMn663c8wzWrYGtgNuqaocug0mShtPmDH77/vkkzwee1FkiSdJItLmKZiNVdTpeAy9JS16bLppD+2a3AKa4p8tmru2OA55H74Ehjx04oSRpIG2uoukfF/4uYB1wSIvt/jfwEeCEBaeSJA2tTR/8QOPCV9X5SVYMsq0kaXhzPbLvXXNsV1V1zCgCJFkJrARYvnz5KHapjq1YddZiR/it8dt4rNe998DFjjAx5vqQ9bYZXgCvBd4+qgBVtbqqpqpqatmyZaParST91pvrkX0fnJ5Osj1wBPBq4CTgg7NtJ0laGubsg0+yC/Bm4I+B44EnVtWN4wgmSRrOrF00Sd4PXAjcCjyuqo5eSHFPciLwbeBRSa5J8tqh00qSWpvrDP7P6Y0e+U7gyCTTy0PvQ9Y5hyqoqpePJKEkaSBz9cEv+C5XSdLSYRGXpAllgZekCWWBl6QJZYGXpAllgZekCWWBl6QJZYGXpAllgZekCWWBl6QJZYGXpAllgZekCWWBl6QJZYGXpAllgZekCWWBl6QJZYGXpAllgZekCWWBl6QJZYGXpAllgZekCWWBl6QJZYGXpAllgZekCWWBl6QJZYGXpAnVaYFP8pwkP0zy4ySrumxLkrSxzgp8ki2BjwLPBR4NvDzJo7tqT5K0sS7P4J8E/Liqrq6qO4GTgEM6bE+S1Oc+He77IcA/981fA/z+pislWQmsbGZ/meSHLfa9K3DD0Am7ZcbRMOPwlno+6MuY9y1yktl1dhyH/J4fOtsbXRb4zLCs7rWgajWwekE7TtZU1dSgwcbBjKNhxuEt9Xxgxq502UVzDbBH3/zuwLUdtidJ6tNlgb8Q2DPJw5JsDbwM+EKH7UmS+nTWRVNVdyV5I/AVYEvguKq6fES7X1CXziIx42iYcXhLPR+YsROpule3uCRpAngnqyRNKAu8JE2oJVngk7w4yeVJNiSZ9bKk2YZCSLJLkrOT/Kj5unMHGedtI8mjklzc97olyZua945O8i997x2wGBmb9dYl+X6TY81Ct+86Y5I9knwjyRXN78URfe91dhznG2ojPX/TvH9pkie23XaMGf+4yXZpkm8leULfezP+3Bch4x8mubnvZ/iuttuOMeNb+/JdluTuJLs0743lOA6kqpbcC/j3wKOAc4GpWdbZErgKeDiwNXAJ8OjmvWOBVc30KuB9HWRcUBtN3p8BD23mjwbe0vFxbJURWAfsOuz32FVGYDfgic309sA/9v2sOzmOc/1+9a1zAPAlevd87Atc0HbbMWZ8CrBzM/3c6Yxz/dwXIeMfAmcOsu24Mm6y/kHA18d5HAd9Lckz+Kq6oqrmu6N1rqEQDgGOb6aPB57fQcyFtvFM4Kqq+mkHWWYz7HFYEsexqq6rqoua6VuBK+jdKd2lNkNtHAKcUD3fAXZKslvLbceSsaq+VVU3NrPfoXc/yjgNcyyWzHHcxMuBEzvIMXJLssC3NNNQCNP/6B9YVddBrzgAD+ig/YW28TLu/UvxxuZP5+O66P5YQMYCvppkbXpDRyx0+3FkBCDJCmBv4IK+xV0cx7l+v+Zbp82248rY77X0/uKYNtvPfZTaZnxykkuSfCnJYxa47bgykmQ74DnA5/sWj+M4DqTLoQrmlORrwINmeOvIqvq/bXYxw7KRXvM5V8YF7mdr4GDgHX2LPwYcQy/zMcAHgdcsUsb9quraJA8Azk5yZVWdv9Assxnhcbw/vX9Yb6qqW5rFIzmOMzU3w7JNf79mW6fz38152r/3iskz6BX4p/Yt7vTnvoCMF9Hruvxl8xnK6cCeLbcdhYW0cxDw/6rqX/uWjeM4DmTRCnxVPWvIXcw1FMLPk+xWVdc1fzJfP+qMSRbSxnOBi6rq5337/s10kk8AZy5Wxqq6tvl6fZLT6P3Jej5L6Dgm2Ypecf9MVZ3at++RHMcZtBlqY7Z1tm6x7bgykuTxwCeB51bVL6aXz/FzH2vGvv+sqaovJvmfSXZts+24Mva511/iYzqOA9mcu2jmGgrhC8BhzfRhQJu/CBZqIW3cq8+uKWbTXgBcNtJ0PfNmTHK/JNtPTwPP7suyJI5jkgCfAq6oqg9t8l5Xx7HNUBtfAF7ZXE2zL3Bz0800rmE65m0nyXLgVOAVVfWPfcvn+rmPO+ODmp8xSZ5Ery79os2248rYZNsReDp9v6NjPI6DWexPeWd60fuHeg1wB/Bz4CvN8gcDX+xb7wB6V1RcRa9rZ3r57wDnAD9qvu7SQcYZ25gh43b0fll33GT7vwO+D1xK75dpt8XISO/KgUua1+VL8TjS61ao5lhd3LwO6Po4zvT7BbweeH0zHXoPtbmqyTA117Yd/VuZL+MngRv7jtua+X7ui5DxjU2GS+h9EPyUpXYcm/lXASdtst3YjuMgL4cqkKQJtTl30UiS5mCBl6QJZYGXpAllgZekCWWBl6QJZYHXokmyrvm6Ism/NaPx/SDJCc2NTaNo4+gkbxnRvtaNYj8LaO/czDGaqjQfC7yWiquqai/gcfTuJHzJIufZrCVZtLvUtXRY4LWY1m+6oKruBr5LM9hTkoOSXJDke0m+luSBzfKjm8HFzk1ydZI/m95HkiPTG9v7a/SGnZ5evleS7zQDk502PTBZs4+/TnJ+emPO75Pk1PTGqH/PXHmbvz6uSPKJ9Maq/2qSbfv2O9VM79r3F8urkpye5IwkP0nyxiRvbr7H76QZZ7zxJ+mN435Zc5fn9N2TxyW5sNnmkL79npzkDOCrA/w8NGEs8Fo0VbXPpsuSbAP8PvDlZtE3gX2ram96w7i+rW/1fwf8Eb2xP45KslWS36N3q/newKFAfxsnAG+vqsfTu/P0qL737qyqpwEfp3cr+huAxwKvSvI7s+Vt7Al8tKoeA9wEvLDFt/9Y4D822f8b8Kvme/w28Mq+9e5XVU8B/hQ4rll2JL3xyPcBngG8v7lNHuDJwGFVtX+LDJpw/hmnpeIRSS6mVyxPqapLm+W7A59txpzZGvhJ3zZnVdUdwB1JrgceCPwBcFpV/QogyRearzsCO1XVec22xwMn9+1reuyR7wOXVzOEcZKr6Q1E9Qtm95OquriZXgusaPH9fqN6Y9vfmuRm4Iy+9h/ft96JAFV1fpIdkuxEb7yTg/s+W9gGWN5Mn10bj3So32KewWupmO6D/11g3yQHN8v/FvhIVT0OeB29Yjbtjr7pu7nnhGWQ8Tem97Vhk/1uYP4Todly3MU9/8a2YWObttHffn97m34v08MRv7Cq9mpey6vqiub92+bJqt8iFngtKc2Z8yruGTt/R+BfmunDZtxoY+cDL0iybTPK30HNfm8GbkzyB816rwDOm2Ufo7IO+L1m+kUD7uOlAEmeSm+0ypuBrwCH943AuPeQOTWhLPBaik4HtmuK8dHAyUn+Abhhvg2r92i/z9IbOfHzwD/0vX0Yvf7qS4G9gL8cce5NfQD4L0m+Bew64D5ubLb/OL0HdkDvwSZbAZcmuayZl+7F0SQlaUJ5Bi9JE8oCL0kTygIvSRPKAi9JE8oCL0kTygIvSRPKAi9JE+r/AznxDlqHIowRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# A histogram is a great way to visualize data\n",
    "# matplotlib is a library of plotting tools\n",
    "plt.hist( rng_array )\n",
    "plt.xlabel(\"'Random' number\")\n",
    "plt.ylabel(\"Number of occurrences in array\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Interactively: What's a different way to convert to (-1,1) random numbers?\n",
    "# Make a new array and plot in histogram to see if it worked\n",
    "# Reminder: this is what we used before: rng_array.append(  2.0*rng.random() - 1.0  ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Position')"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29eXSjZ333/b203Npla/Nujz0znsw+mWRIQjZSmkAS0rAU+qYspzS04XlbWmhpKZSXp30OPPRAVwqlFMpWSqFQCDuBEEIgCUtmMpl983jsGa/a91263j/uxbIt2bK13LL0+5zjM/ItWbrmlnR/r9/OOOcgCIIgOg+N2gsgCIIg1IEEgCAIokMhASAIguhQSAAIgiA6FBIAgiCIDkWn9gI2gtvt5qOjo2ovgyAIYktx7NgxP+fcs/L4lhKA0dFRHD16VO1lEARBbCkYY9PljpMLiCAIokMhASAIguhQSAAIgiA6FBIAgiCIDoUEgCAIokMhASAIguhQSAAIgiA6FBIAgiA6hqcv+THhjau9jJaBBIAgiI6Ac44/+OIxfPTHl9ReSstAAkAQREdwLZhCNJ2HP55Reyktg+oCwBjTMsaOM8a+o/ZaCIJoX07PRQAAgXi2qa/ri2WQzhWa+prVoroAAHg7gHNqL4IgiPbm9KwoAMFE8wSAc477PvIzfOqnk017zY2gqgAwxoYAvALAv6u5DoIg2p/Tc1EAogA0axZ6PCO6nKYCyaa83kZR2wL4JwDvAlCs9ADG2COMsaOMsaM+n695KyMIom3gnOPMbAQaBuSLHNFUvimvK7ubAonWjDuoJgCMsQcAeDnnx9Z6HOf8k5zzI5zzIx7PqnbWBEEQ67IQTSOQyOLQcDeA5l2Q5YBzs+MO1aKmBXAbgAcZY1MAvgzgpYyx/1RxPQRBtCmnZ0X3z0t2iZvIZsUBlgSALIBlcM7fwzkf4pyPAngIwI85529Uaz0EQbQv04EEAOCmUScAwN+kHbn8Ov4mxh02gtoxAIIgiIbji2cgaDUY81gANN8CyOaLiGeaE3fYCC0hAJzzn3DOH1B7HQRBtCe+WAYemwFOiwAACDY5BgC0ZhygJQSAIAii3pyaieAdXz6OQpHDH8/CbRVg0GlhNegQaJIFUHrRb8VMIBIAgiDakifOL+IbL8xhNpRSLAAAcFqEprqAbEaddJssAIIgiKYgu19mwkn44xm4raIAuKxC09wx/ngW1/XaAJALiCAIomn4YpIABFMIxJcsAJdFaJoLyB/PYFwRAHIBEQRBrMmEN6b07akFWQBOzUZQ5FAsANEF1PiLcTpXQCydx2C3ETZj8+IOG4EEgCCIluJd/3MSf/bVE8rviUwer/vEszgr9fKpFtnnfmImDAAlMQBDU/oByRd8t9UAt9WwLCOoWGyNmgASAIIgWoZ4Jo8TMxHMhlLKsfMLMTw3FcKzl/1VPw/nXLEAzs2LwqHEACwCcgWOaLqxefmyy8dtNcBVEnh+8oIXN3zgcUz5Ew19/WogASC2HH//wwv46UVqDNiOPHcliEKRI5bJI5rOAQBmw6IYzIXTVT9PIltASurBnyuIu20lBmCVawEa65KRd/wuq7As8Hx6JoJwMof/+z31u+CTABBbCs45/u2nk/j2iTm1l0I0gJ9PBpTbshUwExJbKc+FU2X/phzy7r9HuugDgFu68DerGMwfW3IBuawGpQ5gISoK2eNnF/H0peqtmkZAAkBsKZLZArL5IkLJ1guoEbXz7GU/7FLevHzBl4VgPlK9AMi77+ul7p9GvQZWg/i8LosoCo1My+Sc42cTfmiYaHnILqBCkWMxmsYOjwUui4CvPT/TsDVUAwkAsaWQzfZQMqfySoh6E05mcWYuigevHwBQIgDKv9W7gGQLQG7/7LYawBgDADglS6CRWTl/8/3z+PaJObzt13bCqNeix2ZAkYvCtBBNY8RpxojLrKxTLUgAiC1FWLrwh1owpY6ojfMLMXAO3LO3D4JWo1zwZQvAH88gk69utq58YT0sCYCnxBXksjQ2BpDOFfCpn03iNYcH8Sf37AIAjLjEJnTTgSQWIhn0dRnhshhUH1BPAkBsKYJJ2QIgAWg3ZFHvsRnQ323EXDgFzjlmQik4zHoAwEKkOivAH89Aq2HYN9gFAPBYlwTAqNfCImgb5gLyxTLgHLhlu0uxOkZdZgDAJW8MgUQGvXYj3NbmFaRVggSA2FLIF4lwKodCi+RSE/VBFnenRcBAlwlz4RRCyRxSuQKOSH38Z6sMBPtiGTgtArpMevTYDBhymJfd77RWVww2H0nhygbTNX3Srr7U6hjsNkGnYXjuShCcA312I1xWMS6gZk0ACQCxpZDNds6BSGp1HCAQz+B93ziNdK46VwHROsji3m3WY6DbhNlwSnH/vGjUAaD6VFBfLKPs+r/6v16Mt989vux+p8VQ1e77fd84jQc/+jQmvLGq/x/+2GoB0Gk1GHKY8MsrQQBAr+QCKhR52c9xsyABILYU4RLXTzk30FMXffjCL6br0kqAaC7BRA5Wgw4GnRaD3UYsRtOYkiZ5yRbAfJUWgL+k9882lwVdJv2y+12W6hrCTfoSiGXyePhzR6u+UJezAOR1zEsuLNkCANRtE00CQGwpgqUCUGYHJwf/WrH1LrE2oWQWDot4oR7oNqHIgWPTIQDAdrcFbqsBc1WmgvpiS90/y1FNS+hikWMmnMKh4W5cDSbxVJXFh75YBowt1RvIjLktyu0+u1FZn5qfVRIAYksRSiztwsp9gb2KALRe50VibYKJLJxm8aI50G0CIPb0twhadJn0GOg2VpUKyjmHr8QCKIfsf1+rH5AvnkE2X8Tdu3sALLl21sMXy8BpFqDXLr+8bpMCwYJOg26zfskCIAEgiOoIJbPo7zICWEoJLcVHArBlES0A8aK4u88Gm1GHa8EUbtjmAGNMCQyvRySVQ67A1xYAi4BsoYjYGnN6rwXFCuT9g13QaVjVn6nS4TOljEqpoH12IxhjSwVpKrqAdKq9MkFsgmAii+0e0ZcaLBMDIAHYugQTWez0WAEAPXYjTv7Vy5AtFKHXiPvUEZcZT17woljk0GhYxefxK03YhIqPcUoX32A8C7tRX/YxM1IAethp3tAQmUrWh2wB9NrF+xxmPRgTXUB/8t8voL/LiHfdu7uq16gXZAEQW4pQMovBbhMEnUYJAp+4FsYDH/0ZoumcEoBrxelLxNqEEksWAAAwxmDQaZWL/Zjbgky+uG4cwFsmC2clcjHYWplAsgUw5DDBZTFUvVMvzUAqZchhhoYBvXbRgtVpNXCYBcyHU/juyXn8+9NXmjaqUoYEgNgycM4RSuTgsAhwmgUlCPyNF2ZxejaKs3NRsgC2KJl8AYlsYVXgtBQ5iDrpWzsvv1wjuJU4y1QDxzN5/M33ziGVFVOIr4WS8NgMMOq1cFkF+KrYVMhtqMuJj6DT4KGbRvDyfX3KMZdFwDMTfmQLRWTzRXzpV1fXfY16QgJAbBmS2QKyhSKcZgHdZj2CUkD4mQmxo+KFhZiSqkdZQFsLOZ7jMFcWgO0eUQDWK8yS33uP1VjxMXIA9uJiDB9+7Dyy+SKemfDj3346iV9eETuSXgumMOwwSc9lqGqkYyyTRyZfrGh9fPDVB/AbhwaWrWNOSg3d1WvFf/5iGrlCcd3XqRckAMSWQd6tOSwCnBYB4WQW3mgaFxfjAIDnpsQiG7OgrTpjg2gNlPfWXN4fD4gXYatBt64A+GIZCFoN7KbKIU45APvRH1/Cx39yGSdmwiUZZOJaZsJJpYLYZRXgj2fWnSLmq8L9tGwdkqvIbRXwtpeOYz6S3vDks1ogASC2DLLP32kW4DALCCazeEaaEmXSa/Erqcryuj4bYpk8VQNvIUIl4l4Jxhi2eyy47Iuv+VxiDYCg9OEph0nQwqTXIp0Td9sLkTR8Up9+XyyDfKGIuXAaw07RAnBbDUjnikhm1/5MKQKwRg1CKW7p/3toqFuxNpoZByABILYMpRaAw6JHKJHF05cCcFoE3D7uVnZwe/rtABrb7peoL6V9gNZizG1Z3wJYpwZAxmkRYBa0AIDFaBqLUfHz442lMR9Jo1DkGFYsgOpmCGzWAjg41K3835v5uSUBILYMS35iPZxmAeFUDt8/PY/bd7oV/zAA7OmzAai+cIdQH8UCWCMGAADb3VbMhlNrWnf+daqAZf7fu3bg7153CCa9FguRNLyxJQvgmjSFbNi55AIClto8VGLjAiBZAMNdigA0s9U5CQCxZZB3Rk6LAIdFAOdicc1f3r8HY64lAditWAAkAFsFOaDfvUYMAADGPBZwLvbVr0S1FsAbb9mG+w/0o6/LiIVoWrEgfbGMUgMwVBIEBkRX0fu/cxZXK7y+P56BTsNW9R6qxB07PXjgYD9uGnPCatBBr2VNtQCoEIzYMnhjaQg6DbpMety7vw/RVB6/d8cYLAadkiLotAjok/Ks5ZmsROuSzRfx/NUQ/PEM7EbdqvYJK9kuvc/Hr4ZwnWTplVIocgSqFACZXrsBi6UCEM9gJpiEhi21pJB36t8+MYfHzizg3HwUX/y9m1fFGeRq5rXiD6WMuMz42OtvUH53WoSmWgAkAMSWYTGSRq9dHO3X32Va1uJXFgCP1aCY/+uZ64S6HJsO4s+/ehKT/gT0WoZB6WK7Fjt7rNjuseA9j55CMJnFH9y1c9n9oWQWRV69CwYQWzM8NxVS0jxFF1AK/V0mRZDkrKEfn/cCAJ69HMB3T83jgYMDy56rtJ/RZnCYmzskhlxAxJZhIZpGr618brfHZoBF0MJjM8AkiBOfqBistXnvo6eRyhXwltvHkC/ydQPAgDjN61tvux0v2eXBPz5+EfkVOfOyD76aGIBMb5cRs+EUihzo7zIils5jwhvHoGNJkASdBnajDtlCES/e7sLefjv+6UeXVj2XWKhYnfunHK4qB9XUCxIAYsvgjWbQ21VeABhjeODgAO4YdwMA3DYDtYNoYbL5Iia8cbzq8CDe98Be/Odbbsa779tT1d9aDTrct78PuQJfNSBmo0FYAMs2FfsGxPjRufmokgEkI4vKXdd5cOcuD64GkqvqAkLJbFVCVgmHWUAomUMgnsHb/uv5ZfMvGoFqAsAYG2aMPckYO8cYO8MYe7taayFaH845FqJpxb9fjg+99iDe+pIdAMQvdbXzY4nmMxVIIF/k2NUrNn+7bacbN405q/57ubPmlcDylFDZ6qs2Dx8A+ko2FXsHxBnC+SJXagBkZAG4c5cHHpsB2UJx1ZCYUDKL7hpcQOKgmgyenvDjOyfn8dxUaNPPVQ1qWgB5AO/knO8BcAuAP2SM7VVxPUQTePK8d1P9TmKZPJLZgtJJcT2GnCYllY9oPS5J1dvjPasDudUgx3ymVtQEyIFc94aCwEsCsF+yAACsmiPc321Er92A3X02xcLwlaQaF4scoWSuthiARUA0nVfOz2yDP8OqCQDnfJ5z/rx0OwbgHIBBtdZDNIfP/3wKH31ite90PbxSlWbvGhZAKSNOMxaiaaoGblEuLsagYWJQdzPIMZ+VRWFT/gScFgFWQ/X5LcstgCUBGHYstwDec98eJfNHbjTnLRGAWDqPQpGvWc28HnKX0uPXxJ3/XIOt2JaIATDGRgEcBvDLMvc9whg7yhg76vNVN5KNaF0C8Sx88QyKxbV7qqxkISJ+0dZyAZUy4jSDc2C2yhmyRHO55I1hxGmGUa/d1N8zxrDNZVFmBstMeOPKTIFq6bEZwJhYYCgOaxGPy0VgMn1dRuyULJZyFoDcqmStfkbrIYvH8athAI3//KouAIwxK4CvAXgH53xVFyTO+Sc550c450c8Hk/zF0jUFX88g1yBlx3ovhYLm7AAAOBqkNxArUA4mcW//uSyIvwXF+MY792c+0dmzG1Z5gLinGPCF8fO3o0JgF6rgctiQI/NCJ1WA5dFgF7L1vysLVkASzt0uZ1FLRaAHECWew7NhtpYABhjeogX/y9yzr+u5lqIxsM5VzJzvBts07AoCUBfhSyglcgCcI0EoCX47ql5fOix8zi3EEU2X8SUP6EEgDfLqNuMa6GU0j45kMginMxt2AIAgFGXWdnxu60GDHSboF1j6pjVoINRr1luASSWmhVultIMIg1DVSMwa0G1QjAmlsp9GsA5zvk/qLUOonlEU3lkpS/rYjStNG2rhsVoGl0mfdUuA4/NAINOQwLQIsitFRajaeg0GikDqDYLYNRlQaHIMRNKYcxtUQKnm4krfPT1h5UL/h3jbuTXcVEyxuCxGZYJQLDKfkZrUSoAh4a7cfxqGJl8AQbd5lxl66FmJfBtAN4E4BRj7AXp2F9yzr+n4pqIBuIvKXDZqAWwIFUBVwtjDCNOM7mAWgRZABYiGaUF845N7NRLKc0EGnNbMOHbvAD0dy0FfN/7iuqSEXtsxmWfY6VZYQ2FYKXicftON45fDWMhksa2kl5X9UQ1AeCcPw2guoYZRFtQ2p3Tt1EXUCxTtf9fZthpxtUgBYFbATmdcSGaRjwjXihXBlk3ijIi0p/ArwG47I3DImjRX6WbsFY8VsOy2QTBZBZ6LdtQBtJK9Fqx4jiVK+BFo2JdxGw41TABUD0ITHQOpWMaZZ8+IPo518sKWoysXQRWjhGnGdeCSfzVN09vKvWUqB+KCyiSxkwoBbtRV3XHzEo4LQJsRp0SCJ7wxrGzx1p1I7Za6bEbllkAoYRYBFbr67usBow4zUoca2W1cz0hASDW5JeTATz8ueeQzdc+p1Su0nRaBHil4RuhRBZ3/e1P8K0TcxX/rlDk8MUz6NmACwgQd5jxTB6f//k0vnd6YfMLJ2oinSsoF8rFmCgAK4usNgNjTMwEklJBL3lj2LHJuoLN4LEaEEnlkMmLGTu1NoKTOTTUhVt3uNHfLW54GpkJRAJArMmjx2fx4/NenJoN1/xc/ngGjAHX9dqU9Lm5SArZQhFn5yvPQQ0lsygUOXoqNIKrxJhbvMgYdBr4YtQWQi3mS4qZFiJpzISSyxqt1cKoS5wQFknmsBjNbLqyeDPItQCyZRtO1tYITuafHjqM979qPww6sblhIzOBSACINTk2LVYk/vxyoObn8sfFHVJ/l1EZvydnTkwHKo/5k+MFPRso7weAO8Y9+PgbbsBbbh9DIJFd1TmSaA4zkv9/vMeKxahsAdRJANwWzIVTeF6qnN0/WH1mWa3IFqlcpR6ssRFcOQa7TevOQK4FEgCiIpFUDpe84ofvF5PBDf1tMptHIpNfdswfF0f1eexi+hznvEQAKmfreDfR4REQA2r3H+jHQLcJnC+PQRDNQ3ZhHBl1IJTMIZkt1MUFBIhWXpEDj50SXXz7pGZuzcBjFS1SeYMixwDqyT17e3F0OoTTs5G6Pq8MCQBRkeNXxV3V7j4bjk4HNxQH+OMvvYA/+tJxAMCEN4bTsxFRAGwCem1GpZOiLABXg6tb68pspsVvKeWqNonmMRNKQathODjUrRyrmwUgZcc8dmYB/V3Guu/A10JOS56PpFEscoRTtTWCK8ebXrwNdqMOH/vxRF2fV4YEgKjI89MhaBjwyJ3bkc4VcXKm+jjAhcUonrsSRLHI8e6vncKbP/scFiNpuK0GxXRejGYUAUhmCxUneG1myEcpPVL2kBx4JhrDG//9l3j0+Myq4zOhJPrsxmUTv+olAHIqaCSVa+ruHxA3JA6zHmfnopiPplEo8g3VqlSD3ajHm28bw2NnFnBhIVbX5wZIAIg1eP5qGHv67Xjp7h4wBvxisro4QLHIsRBJI5bJ47IvjpPS7n8uklZ6rgDijrx0/F2lQdveWBoWQQvLJvOr5S/lRovPiOpJZvN4esJfNlY0G05h0GFa1sajXi6gbrOgDJLfN9A8/z8gZiHtH+zCqdmIYi0fGu5e5682zsO3jeLOXR6l5UU9IQEgKnLiWhiHR7rRbRbQbzfiir+6qlp/Qmz4BgBfPTazzHXkti0NbZ+PpBGMZ2HQiR/DSnEAXyyj7OI3g9sqdnskF1B9SGULq9x1snW1uMLKyhWKuOSNY5vTrBTy2epQA1CK7AZqtgAAwMGhLlxcjOEXkwEYdJoNtTeplm6zgP94+CbsH6y/hUMCQJQlmc0jlskrOzW3zVD1jN35ksKVrxy9BgB4zWFx1IPbakB/txEaJvqGg8ks9g7YoWGVM4F8scyGJjytRK/VwGkWyAKoA8evhnDD+x/HV48td/XIhX2lBX4A8NOLPoSTObxsXx/sRh1Mem3ddv8y2yU30L4GXCDX48BgF/JFjm8en8PBoS5liPxWYWutlmgactdOOajmthqqbt8wHxGzPkx6LcLJHHpsBrzz5ddhV68V1w93Q6/VoL/LhJlgEsFEFn12Iwa6TZiu0LfHF89sOgAs47EZysYAFqNp3PD+x6t2b3UyM6Ekfv8/jiKVK+Ds3PK6jYUKAvD147NwWgS8ZJcHjDEMO01KfUa9uHtvL166uwcDTWoBUYq8K49l8rhhxNH0168VEgCiLLJv3qUIgFC1BSCXrt+5SxzQfnikG4PdJvzwT16idIAccogjG4OJLBwWAdtc5souoGjtAtBjN5YtBnvinBfBRFapdyAq84mnLiOeyaPXbljVZE+p7E4uVcZG0zk8fnYRv3GwH4Lk5vv4G27E+x6o7+TX+w/04zNvflHTWkCUMthtUjZJh0fq7/9vNCQARFmCiaW2DYBoAQQS2aomec1HUjDoNLhjXBzgc7jMzmjYacZUIIlQMguXRcCI01K2dXMqW0Ask69dAGzL+7Zck9JOn7zgBYBVowWJ1VxYiOHAYBcODXWveq9Kd/6yGHzz+Cyy+SJefcOQct/OHuuyzptbHTkQDJT/nLc6JABEWWQXkJx66bEZUJBynddjLpJGf5cRt+5wwWbQ4c7x1ZPchh1mqRhMFJk+uxGBRHZVpoNsddRDAHwxcRTlZV8cd3z4Sfz7z67g2Qk/gLUrkQlxmM8lbxw7e2wYdppxLbS8bmMxVtrqO41svohPPDWJwyPdODTUfN98M3nw0ABetrd3w91qWwE15wEQLYycn19qAQDiBXm9Ypv5cAr9XSZs91hx6v+8vOxjhp1Lu0CnRVAGXgTi2WXpgputAl5Jj82AfFEcRTkhVTf/zffPocjFNNGpNSqRiZJpWz1W6DQM6VxRbNAnpfQuRtOwG3WIpvNYjGYw4Z3BbDiFD7x6vyqumWby2huH8Nobh9Z/YAtCFgBRlkBCTM80C+KFWRGAKgLB85G00smwEqW94F0WQ9kh2+LvomuhliwgoKQYLJZRWhNoGIOg1eB1Nw7DF8usal1BLCFP2xrvsSrifa1k1oI3mlYqfRciafzbU5M4ONSFu3bRHO9WhgSAKEsgLvrm5d2bxybu+itV68rkC0UsRtMYWMfPO1ySCui0CEsCEF8eqN1sI7iVLJXtpzAXTsGo1+CDrz6AP/71ndgr5Y9PkRuoIhNesQp1vNe6at4y5xyL0Qyu67NB0GpwdDqISX8Cr71xqO13/1sdcgERZQkmMnBal1w9sgWwXiqoN5ZBkWNdC6DHZoCg0yCbL8JlFZQZrCuffyGahlbDau7xMuYW+8RP+hJiZWq3Cb/1omEAwJk5sdHWdCDZ9HYCW4UJbxxWgw59diMyUmGfLADRdB6pXAF9diN67Ab86KwYWL91h1u19RLVQRYAUZZgIguXZWnX3WXSQ69l63bUlGsA1rMANBqGIak3TLdZD7ckNitz9acDSQw5TNDVWGDjtAhwWgRclgWgxAKRx+2RBVCZS944dkjTtox6sU+9nAoqt0Pu7TKi1y42+uuxGbDD05gxhkT9IAEgyuKXXEAyjDG4LOtXAx+/KjaM29W3/mCOIacZNoMOBp0WBp0WXSb9KhfTdCCpuBxqZYfHgsu+OGZDKQyWWChWgw4emwHTVba66EQueeMYL5m2NSJlAgFL7R96bQbF1XbrDhe5f7YAJABEWYKJ1cMt3Lb1i8GenvBju9uyrPNjJV56nQcv3dOj/O6xra42ng4ksM1VLwGw4tx8FIFEdtX6Rl1mXCELoCyRVA6+WAY7SwRg2GFSgsByDUCv3ahkBd26k9w/WwESAGIVqWwBqVxhWQwAEDNx/PEMfnLBq6RSlpLJF/DLySDuGK/uy//m28bwkYcOL3v+UgEIJ7OIpvNKs69a2eGxIpYWM31WjiTc1WvDiWvhhg3e2MrIWVOlltiw04z5SAr5QlFpA9FjNyhtnm/d4Wr+QokN03EC8NRFHz74vXNqL6OlCUhVwK6VFoDVgElfAm/5/FF85IlLq/7u2HQIqVwBt5cp/KoGj82wzAUk5+bXzQXUsyQkK2MUb797HC6LgEf+42jVLS86hYWoKACl9Rn9XSYUuVgANh9Jodush1nQ4aGbRvDlR26pe8M3ojF0lABEkjn86X+/gE/9bBKFIsfp2Qj+/ocX1F5WyxFU+gAtT7102wxIZgsoFDmulnGXPH3JD62G4Zbtzk29bs8KF5BcnTvqrp8FILPSAuixGfHxN96IuUga3z81X5fXaxfk3k6lojkgxVDmwinMhZfSfq0GHW7ZTrv/rUJHCcCHf3AegUQWnIt+zW+fnMNHfzyBOBUALUPpBLrCBSTn4ntsBlwLpVb93TOXAzg83A2bcXO93j2SwMgFWVfrbAEMOcwQdBpoNUyZSVDKgcEuMAb4aHbwMhYiYipuaTX2gBRDEQUgpfxObC06RgAiyRz+61dXleBfKJlFUPqir2xh2+ms7AQq8xuHBvCBV+3Hm28dRTCRXSacqWwBZ2YjuGlsc7t/YKndg9z+YSogjhI06rWbfs5StBqG7W4L+uzGsmmlWg1Dt0mvNMIjROYiKfTaDNBqlrJ6+ruWhvrMhlJ1G/FINJeOEYAJXwycA/fs7QUAhBJZhJIkAOVY2QlUxm014I23bFOycko7Qp6cCSNf5Lhx2+Y7IsoC8NyVIL5xfBbTgQRG6pQBJPPAwX7ct7+v4v1Oi6C4wAiRhUh6mf8fAGxGPWxGHS4uxBDL5BWXELG16JhK4Ms+0Z984zYHPvfsFELJnLLTrXbQSacQSGQhaDWwVpjBK7dxuBZMKiPwjkkzUWtpiSsLwLu+dlI59ro6N9l620vH17yfBGA185G00oIvucgAACAASURBVC6jlIEuE56bDoq3yQW0JekYC2DSl4Bey3BA6t0dSmYRSixZAKdnI3jP109V1e++3QnGs3BZhYqFPHIjt9KhIMemQtjusdTUskHOIR/oMuJd914Ho16DG2qwKDYDCcByOOeYj6TQXyZm0t9tVGoBqqn7IFqPDrIA4tjmssAt7TJDiazyRfdGM/j2yTl86VdX8dY7t9ct62SrEihTBFaKw6yH1aDDTCiFP/3KCzDptXj+agh37+mt6XWdFgEfeeh6vGjUiYFuE37/ju3QaZpbTeq0GGg6WAnhZA7pXBH9ZS7wpbt+EoCtSccIwKQvjh0eKyyCFnotgy+WQVQqClqMZZDJiWPsLnnjJADrCABjDEMOE352yae41gDU5P+XeeX1g8ptNQZsuywCQskcikUOTZPFpxWZj4jxsf4y83blGbx6LVOaBRJbi6q+YYyx1zDGLjHGIoyxKGMsxhiLrv+XrUG+UMTVYBLbPWIzK4dZWDYC0BtNK7+Xq3DtNIKJzKoMoJUMO8247EtAw4BH7tyOLpMet7VB+b/DIqBQ5Iim15981gnIzf3KCYA82rG/y0RiuUWpdov1YQAPcs67OOd2zrmNc746KrRBGGP3MsYuMMYmGGPvrvX5KnEtlEKuwLFd6k7oMAu47BMv9IJWg4VoGtOSP/uS1Pe8kwnGs3Ba1t7Rybn5d4x78Jf378Hx992zbMjLVkUWvgDFAQCUWgCVXUDk/tm6VCsAi5zzuvZPYIxpAfwLgPsA7AXw24yxvfV8DZlJ6WIvt6ftNuuVAOaOHiumA0lkpR7nnW4BpHMFJLIFuKzrWABS3vdvSlk67bIDlF1fIRIAAKIFsLIITEZO/aQMoK1LtTGAo4yx/wbwDQBKziTn/Os1vPZNACY455MAwBj7MoBXAjhbw3OWZVLyU2+XhoI4LQLkZJ89fTacmxe9WXv67ZjwxsE579hWtpWKwFZy34F+zEfTePm+2gK/rYazxAK44k/AZRVg32RlczswH0mjZ0URmExflxFGvUaxrImtR7UCYAeQBPCykmMcQC0CMAjgWsnvMwBuXvkgxtgjAB4BgJGRkU290KQ/DqdFgEP6cnebly5uu/ttwHHx9j17e/HPT1zCXCTdsWatXB29Xjpnr92I99y3pxlLairy/3vCG8fb/ut5aDUMb7h5G973QEOM05bHF8so85RXYtBp8e233b6qrxKxdajKBcQ5/90yPw/X+NrlttirkvA555/knB/hnB/xeDbXZfJP77kOn/vdFym/O8xLO7pdveLgEougVVrYXlrs3DiA0gl0HRdQuyILwA/PLiJX4BhymPHpp68gXyiqvDJ18MUya85jHu+1wSx0TDJh21FtFtAQY+xRxpiXMbbIGPsaY6zWEs0ZAMMlvw8BmKvxOcvisRlwcKhb+V3+kneZ9EoPk1G3RRGDTo4DyLUR6wWB2xWjXguLoMWJa2EIWg0ekuYGh5KdmRXkjWXK+v+J9qDaIPBnAXwLwABE1823pWO18ByAccbYGGNMAPCQ9BoNR3YBOS0CPFL16ZjbosyN7WQBkDuBdqoFAEBxFR4Y6lKyXwId2CAuVygimMiuaQEQW5tqBcDDOf8s5zwv/XwOwOb8MRKc8zyAtwH4AYBzAL7COT9Ty3NWi+wCcpj1sBt1GHGa8aJRsYvlzh5rZwtAIgu9lsFWoQ9QJyAHwI9scyhCGOzAFtHyYBy5RQfRflT7Lfczxt4I4EvS778NIFDri3POvwfge7U+z0aRd3hOiwGMMTz153cp9433WPGdk/MdmwkUTGTgtFTuA9QJyC7CG7Y5FDHwd1Ba6GefuYLLvjh+64jo/iIXUPtSrQXwMIDfArAAYB7Aa6VjWxKH4gISLQHGmHLB29ljFYdgd+hYQHEYfGd/4eX//43bHHBJLQ6CHfR5+MkFH77+/CwWo7IF0Nmfh3amKguAc34VwIMNXkvTkF1A5S504z1LgeBONH0Diey6NQDtzj17e2HQa+C2GsSeQKyzKoPDySyS2QJOXAsDEIe9E+3JmgLAGHsX5/zDjLGPonyK5h83bGUNxG7U44GD/bhz1+reNTt7xGKxCW8ct+7Y+r1tNkognsXIyNZv6VAL9+7vw73S0BiNhsFpETpKAOSMp2cv+wGsng1NtA/rWQBy+4ejjV5IM9FoGD72+hvK3tdrN8Bm0HVsIDi4TifQTsRlMSDQQS4geVLeiZkInBYBgq5jxoZ0HGsKAOf829LNJOf8q6X3McZe17BVqQhjDDt7rbi02HkCkMoWEM/kqbXvCjppSEy+UERMapNeKHLy/7c51Ur7e6o81hbs9FhxqQMtgDmp9S/Nd12Oyyoo9RHtTji1vOCNMoDam/ViAPcBuB/AIGPsn0vusgPIN3JhajLea8VXj80glMgqKaOdwFxYEoAyrX87GVcHxQDCkvuny6RHJJUjAWhz1rMA5iD6/9MAjpX8fAvAyxu7NPWQB52fnd8yM2/qgiIAHdoIrxIuqwGRVE5pGd7OyAHgI9J0t07MhOsk1osBnABwgjH2RalytyPYPyAOjj89G2mLKVfVMhtOQ8PENr/EEsqMgGQWvRU6Y7YL8hyEI6NOPHHeSxZAm7OmBcAY+4p08zhj7GTJzynG2MkmrE8VHBYBg90mnJqNLDvOOcestEtuR2ZDKfTajarM4m1l3FI7iE6IA4QlC+DOXW7s8Fhww0j3On9BbGXWSwN9u/TvA41eSKuxf9COM3PLXUBPXfThzZ99Dl9+5Bbcst2l0soax1w4Re6fMsgFg53QEE5OAd3msuCJd96l7mKIhrPmVo9zPi/d9AO4xjmfBmAAcAgNat3cKhwY7MIVfwKxkuHgcm3Ax348odayGspchASgHEpDuA4IBIeSOQhaDSyCVu2lEE2gWlv/pwCMjLFBAE8A+F0An2vUolqBfYNiHKDUCpAHZD894cfxqyFV1tUoikWO+XCaUkDLoDSE6wgXUBbdZn1HNwPsJKoVAMY5TwJ4DYCPcs5fDXGQe9tSGgiWmY+kMNBlRJdJj889O6XSyhqDP5FBtlDs2FGYa9Fl0sOk1ypZUu1MKJlVmiUS7U+17aAZY+zFAN4A4C0b/NsticdmQH+XEcelhlgAMBdOY8xjgVajwRV/QsXV1Z/ZkHhxIwFYDWMM2z2WjmgPEkrm0F0yMpVob6q1AN4BsfL3Uc75GcbYdgBPNm5ZrcGtO9x4ZsKPQlHsgzcfSaG/y4Q+u0FxB7ULc2Hx/0MxgPJ0yqCgMFkAHUW1Q+Gf4pw/CODjjDEr53xyq3YC3Qh37nIjnMzh1GwEuUIR3lgGA11G9HWZ4I9nkGujQeFUBLY24z1WzIZTSGTauxwmlMzBYSELoFOodij8AcbYcQCnAZxljB1jjO1r7NLU5/adbjAG/PSiD95YBpwD/d0m9NmN4BzwxdonLfCZy354bAbYjW3t2ds0cpvwSV97uf5K4ZxLQWCyADqFal1A/wbgTznn2zjnIwDeCeBTjVtWa+CyGrB/oAs/vejDvLRD7u8yoq9LzAtvFzfQhYUYfnLBhzfdso2yPyogC8Alb0zllTSORLaAXIErA5OI9qdaAbBwzhWfP+f8JwAsDVlRi3HnLjeOXwvj/IL4xR/oNqHPLrpJFqPtIQCf+tkkTHot3nTLNrWX0rJsc1mg07C2jgPIbSDIAugcqhWAScbY+xhjo9LP/wfgSiMX1ircs7cPhSLHZ54R/7uiBSDmyi9UsAD88Qw4XzVArSVJZQv45guzeO2NQx3V+XSj6LUajLrbOxNILnRzkgB0DBsZCu8B8HXpxw2xGKztOTTUhfEeKyZ9CdgMOtiMejjMegg6TVkL4Fsn5nDkAz/C0xN+FVa7ccRgNsfBoS61l9Ly7PRYMeFrXwGQY1rUAK5zWK8ZnJEx9g4A7wdwBsDNnPMbOOfv4Jy3VylsBRhjeN2RIQBAv1QlyxhDn924KgZwfiGKv/gfsUfeC1fD2Aoouz7a/a/Lzh4rpgNJ5Nso+6sUX5wEoNNYzwL4PIAjAE4BuA/A3zZ8RS3Iqw4PQqth6C8ZlNJnN2JhhQXw9z+8CLOghcdmwPnFrREsDErNv8j9sz69XUYUirxtewLJFoDc+4hof9bL+dvLOT8AAIyxTwP4VeOX1Hr02Ix47/17MOI0K8f6uow4MbN8l38tmMThEXGQxoWFrSEAIfL7Vo1HujD64hn0tOFcAF8sg26zHgYdNYLrFNazAJRWmJ00EKYcD98+hrv39iq/93WJLqDSYK83lkGv3YDdfTZc8SeQyRfUWOqGkCdAUfXn+ritomukXZvC+WIZeKzk/ukk1hOAQ4yxqPQTA3BQvs0Y66x5iSvotRuRzReVARrZfBHBhDgx6ro+GwpFjsve1i8aCiWy0GoYbFQAti6KALRRAWApvniG/P8dxnrzALScc7v0Y+Oc60pu25u1yFakX0oFlQPBcgCtx2bAdX02AMCFxdbXyGAyC4dZD42GCsDWw22TLYA2FYAYCUCnQbP/Nsk2lxgPuCylBcopob12I8bcFui1DBcWWj9lMJSg5l/VYhG0MOo1bSkAnHNyAXUgJACbZGePFToNw7l5cZfvjUoWgN0AvVaDHR4rLiw0xgKIpnNI5+oTXwgmspQBVCWMMbithraMASSyBaRyBbIAOgwSgE1i0Gmxs8eqtIjwxkQLoMcmuoYOjzjw7OUApgP1jQNwzvGajz+Lv/rmmbo8XyiZpQygDSAKQPtZAFQE1pmQANTA7j6bYgEsRtPQapgyPvDtvz4OnYbhf3/zzLJMoS/8fApv/uzms2mfvxrChDeO03OR9R9cBdT+d2O4rYa26gIrQwLQmagiAIyxv2WMnWeMnWSMPcoY61ZjHbWyp9+O+Uga4WQW3mgGPTaDEkzt6zLinS+7Dk9d9OFnl5baQnz31Dx+etGnDJnZKN84PgcAmPInau43xDmnGMAG8diEtnQBkQB0JmpZAI8D2M85PwjgIsRpY1uO3f1iItS5+RgWY6IAlPL6m0eg0zD8YjIAQBy8fmY2iiIHAomN7yJzhSK+c3IOglaDRLagZB5tllgmj3yRUxuIDeC2GhBMZDYt4K2KT3JhUhC4s1BFADjnPywpLPsFgCE11lEre/rFdM9z81F4o+lV1aFGvRa7em04JQ2Wnw4mEZMmSslB443wzIQfoWQOv33TMABgyp+sZflKFTBZANXjthpQ5GLspJ3wxTPQahh9FjqMVogBPAzg+2ovYjP02IxwWwVRAKQq4JUcHOrCyZkIOOeKEABLQeONIA+if/UNol5O1TiYnhrBbZylauD2igP4Yhm4rQLVg3QYDRMAxtiPGGOny/y8suQx7wWQB/DFNZ7nEcbYUcbYUZ/P16jlbpobtznwgzMLCCaySgZQKQeGuhBJ5XAtmMKpkt5Bm7EA/NIubd+AHToNw5UaM4xC1Ahuw7ilfkD+WLatZkJ7qQisI2mYAHDO7+ac7y/z800AYIz9DoAHALyBrxHN5Jx/knN+hHN+xOPxNGq5m+Zd9+5GOideCMpZAIeGxPj2ydkwTs1GsEeKGyxuQgAC8SycFgF6rQYjTnMdLAC5DxBlAVWLXA38we+dw5EP/KhthsQvRNLKpDuic1ArC+heAH8B4EHOeW2ObJXZ4bHiD39tJwCgr2v1F2hXrw2CVoPjV8M4MxvFjdu64TDrN+UC8sczigtizG1RXEKbJUwWwIaRd8ln56OIpHI436Biv2YzF05hoLv9OpwSa6NWDOBjAGwAHmeMvcAY+4RK66gLf/BrO/Avr78Bt+1wrbpP0Gmwp9+GTz99BbFMHrftcKPHZoS3ilzyb5+YW9ZW2hfPKi6IUbcF04FkTamgwUQWOg2DzUCN4KrFZtCh26zHgUFxgtqZua0vAIlMHtF0ftm8C6IzUOWbzznfqcbrNgq9VoNXHOyveP/9B/qRyRfxJ/fswsv29uK/fnUV3ioGyv/F107i4FAXvvzIiwEAgXgG290WAKIApHIFLEYzyozijeKLZeCyCmCMAn/VwhjDd//4DritAm754BM4M7v1BWA+kgIAsgA6kFbIAmp73vqSHXjsHXfi5fv6wBirygJI5wpIZgv4xWQQ0wGx6Et0AYkWwA5JCC55Nz94Zj6SxkA37fo2ymC3CQadFvsGunB2vrIARJI5/MPjF5HNt3aweC4sbkbIAug8SABUoMcuthMorlFMJM8ZAID/OTaDRLaAdK6oxACUltM1TB4T/b70pd8sewfsuLAQq5gN9OjxGfzzE5fwqyvBJq9sY8gWQP8mLUli60ICoAK9NgPyRb5mMZF8n1GvwVePziguI1kAXFYDPDYDzs1vTgA455gNpzBIArBp9g3YkS0UMeEt3/b7uakQALR8oHgunAZj2LQrkdi6kACogFwxvFYqqCwAv3FwAAvRNH4utZMoHdi9u8+26aEzwUQWmXwRA/Sl3zR7pZTes2UCwZxz/GpK3Pmfb/H50PORFDxWsY050VnQO64Ccs+gtVJBZRfQS64Tax9+ckEsgnOX9GrZ3WfDxcU48psoSJoNy4E/sgA2y3aPFQadRukIW8p0IKk0WKvFTdcM5iNp9NPnoCMhAVCBXskCWCsQLFsAh0ccMOm1eHZC7ChaWq25u8+ObL6IqU1UBM+RANSMVsPQazeWbcon7/5fssuDi4uxTYl0s5gLp9BvJ0uwEyEBUAH5Ir4QWd8CcFkE7B2wI5EVJ4CV9u2RA8HPXw3jmy/MolAU2zu/5XPPrVtoNitlflAMoDacFkHpqVTK0akgus16PHBQTAGeCrRmvSPnXLIASAA6EaoAUgGjXgunRcDCGrUAoUQWJr0WRr0W+wfsODYdQrdZv8xPu7PHCq2G4b2PnkKuwDHYbUI8k8cT5714fjqMe/f3VXz+uXAKZkGLbmoDURMui4D5MkL+wrUwbhxxKK0/LizEsLPH2uzlrUs0lUcyW8AApYB2JGQBqESf3bimBRBK5pQePfsGxKpT94pe7Ua9FuM9VqU3vS+WUXaj5XalpcgpoFQEVhsOi1A2m2s+ksaw06yIdKtmAs3JKaBkAXQkJAAqMdBtLLtzlAkns+iWerPvGxR3ka4yPXs+8cYb8ZW3ipXCwWS2RADWLjSjGoD64LIICCSyy1pypLIFxNJ5eGwGGPVajLktLZsJdDUouqaGHWaVV0KoAQmASvR1GbEg7b7KEUpmlVm94z1iQzl3mXa9o24LDgyJFkIwnlV2o4F1LIDZcBqDtOurGadFQDZfVGI0wFJ2l5ztdV2fraIFEEvn8PDnnsO1oDoxArmh4JjHosrrE+pCAqAS/V0mhJI5pHOFsveHkznFAhB0Gvz5y6/D/3NkuOxjDTotbAYdAolsVS6gdK4AfzxDft86IHdSDZWcbzm7S6732NNnw7VgCvEyraPPzkXx4/NePDPhX3VfM5j0xeG2GmA3UiyoEyEBUIk+6eJQKQ4QSmaX9en//Tu3485dlechOK2iL7oaAZCH0VDlZ+3IbrlSi0s+v0sWwFIgeCWyWMhpuc3mij+hNBgkOg8SAJWQ+67MlXEDFYoc4VRuQ/NZ5XTEkDTkZS0BWJRcFL2U+10zclpuacxlpQtot5SuW84NJBeLyWm5zeaKP4ExEoCOhdJAVUKuvCxnAURTOXAOxQVUDU6zmI6YlQqOqrEAespMMCM2xpIALDXv88Yy0GuXBqwPOUywGnRrWgCz4ebHACKpHPzxLLaT/79jIQtAJWQXULlMIGVW7wZy9J2W5S6glZkppSxK9Qe9ZWYYExujrAUQzcBjNSgD1hljYiC4TOM+n+ICar4FII8UJQugcyEBUAmTVIRVzgIIJeVZvRuwAKwCAvEswsksjHrNqsyUUryxDASthorA6oDVoIOg1SyPAcTS8Kxwr+2WMoFWirLcRmI+klqzPXgjkDOAyALoXEgAVKTPvlQLkCsUccWfQL5QVGb1buQC7TQLyBaKKHJgu1usOA3Gy7uBvNE0PDYDFYHVAcYYHBb9snPtjWYU/7/M7j4bouk85iNpZPIF/OPjF5HI5JU237mCOPCnmUz64tAwYNhJNQCdCsUAVKS/y4iFqBgE/vyzU/jAd8/BqNcoF/CNBoFldvZYcXY+imAyixHX6i+3N5Yh/38dcVoMy6qBvbE0jow6lj1mh0d8T6f8CUz5E/jIE5ew3WOBP54Rq8KjacyEU0rqaDOY9Ccw5DDDoNM27TWJ1oIsABXp6zJhXvL9XvbFYTPq8OrDQ5jwijszp7V6ASidEyD3nKlUDeyNpVftUInNI1cDA0A2X0QomUPPiviK/F6GkjkEJbE4Ox9FIJHFoWGxkK+ZqaDpXAHPTPiVIkKiMyELQEUGu40IJLJI5wqYj6SxzWXG37zmAN720p24FkxuqDin1FoYlwQgUMEFtBjN4OYxV22LJxQcFgHXQmIWj+zT711hYcnvTyiZhezpf3YiAM6B64cd+MGZxaYKwDeOzyKUzOFNt2xr2msSrQcJgIrIvXjmwiksRNIYkvqxDHabNtym2WVZuuAsWQCiAPzskg/eaAa/eeMQ0rkCIqncqgsUsXlcJS2hZZ/+ShebHM8JJZYE4PRcBIAYhLUZdJgNNUcAOOf4zDNXsKffjpvHnE15TaI1IReQiiwJQBqL0TT6ujZ/US51Fw05zBC0GuWi9NlnpvAPj18EsJR2uNJFQWwep0VALJ1HNl9cagOx4vwadFpYBC2CyaV+TXJCkMdmwEC3qeHFYHKB2vFrYVxcjON3bx2lRIAOhwRAReRd/qQ/jlAyp9QGbAaLoIWg08Ck18IkiPMGAiVtIfzxDDjnykXAQxZA3ZAD8KFkVsnkWdm6GxBdReFkblnfIECsGB50mDDhjTUsFfTZy37c8sEnMOmL46JUkPbiHeQG7HRIAFSkr8sIDQOenw4BqK01A2MMTrOgXIxKJ1WFkuIA+Hgmr1QBUxFY/XBL1pcvloE/Jp5zV5kAvtKuI5lblrXlthrw4KEBTAWS+OIvpxuyxkuLcRS5GHi+GkxCp2FKOxKicyEBUBG9VoNeuxHHrooC0F9jd06nZUkAXNYlC0DecfpiGaUKmNJA64e82/fHMwgkMqsmt8l0mwWEJRfQvgE7zIIWXSY9jHotXnn9AG7b6cKHH7ugvEf1RJ4+N+VPYDqYxJDDBF2ZNRKdBX0CVGag24RrQTH4V0sMAAAeONSP+w/0AwA8VgP8sQxyhSKiabENsT+ehTeWgU4jWgtEfZAFIBDPIhDPlh3cAwBOs14Z2uO2GjDea1PmQzPG8IFXHUA8m8dXnrtW9zXKojLpT+BqIEnFXwQAygJSncFuE47VwQUEAH9w107ltsdmgC+eWVag5I9nsCBVAct9aojakQf1+OMZ+OIZuMr4/wEpBpDIgUPMCvqTu8eRLGnXMea24IYRB75/egF/9OvjdV3jYokFcDWYxKHh/ro+P7E1IQtAZeRMIKtBB1sdh3J4bAZk80VMB5a6TPrjGcwEUzT+r85YBC2Meo3oAopnlJjAShxmAbFMHvFMHk6zgLuu61EsNpn79veJfvpAfbuDyj2nzs3HEEnlsM1J/X8IEgDVkccy1jsvX3YtlLYg9scymA4myPyvM4wxuK0G+ONZBBLZZTUZpThKXEPdFdxEL9/XBwB47Mx8Xde4GBVbVKekCXT0GSAAEgDVGXSIFkC9p3PJAnBxcUkAZkIpLEYz2FamPxBRGy6rAQuRNMLJXNkUUADL4i6VYjDDTjP2Ddjx2OmFuq0tLlkd1w93K8foM0AAJACqI7uA+uz1nc/bs8ICcJj1OH4tDAAYod1f3fFYBUVsy6WAAsvnO6w16+FFo05cWozXbW2y++eW7Ut5//QZIACVBYAx9meMMc4Yc6u5DjWRi8HqnZPtsYrPJ1+UxnttSv/3ch1CidpwWw1K2m3FGECJ28dRwQUEiMkAsUweiTJD5DeDHAA+MuqEhonrsxgo/4NQMQuIMTYM4B4AV9VaQytgM+rxkYeux0117sliN+kg6DQIJXOwCNplvYVo91d/St0+FbOASl1AawqA+PfeWAZjdbhQyxbAsMOEIYe5okARnYea24B/BPAuAN9UcQ0twSuvH6z7czLG4LEaMBtOwWERlC+9WdBWzFMnNk/pRbVSDKB0wM9aw37kdODFaLou4xrlIrC+LiPe+bJdMAu0+ydEVPkkMMYeBDDLOT9Bzagah8cmCoDTIihB4RGnmRqANQC3rdQCKC+wRr3YEA7AmkNYZAugUkVwMpuHUaetupZjMZqGzaiDWdA1ZLNBbF0aJgCMsR8B6Ctz13sB/CWAl1X5PI8AeAQARkZG6ra+TkC+6DvMgrIrJfdPY5DPr6DVwLaG26a7igpseSqY3LeplGy+iDs+9CTecfc43vTi0arWthhN19RokGhfGiYAnPO7yx1njB0AMAZA3v0PAXieMXYT53xV7hvn/JMAPgkAR44cae7U7C2OnAnktCwJAKX/NQb5/LqswpoW1lq+fxmbQQeTXlvWApgKJBBIZHF6Nlr12haimbqnGRPtQdNdQJzzUwB65N8ZY1MAjnDO/c1eS7sjWwDdZv0yFxBRfzwlArAWrz48iPV2MYwx9NoNWIyttgDk9NDZDUwPuxZM4mV7e6t+PNE5UDSojZEv+k6zgOt6bXjnPbvwioMDKq+qPbGbdNBrWcUAsMzDt49V9Xw9dmNZC0BO650JVdcqIiw1n9vuodYPxGpUFwDO+ajaa2hX5F2pwyJAo2F1bzBGLMEYw7DDjCFHfQr6eu1GnJoJrzo+4RUtgLlwGsUiXzcQLNd+bHdb67Iuor1QXQCIxiHPF/DY6ttniCjPF37vZljrVGDVazPgR1FxiltpTOGSV7QAsoUifPHMuh1kJ32iAIyRBUCUgVpBtDH7B+341zfcgJfu7ln/wUTNDHab0GWqT0fXXrsRqVwBsZJq4FyhiCv+BHb32QCIvZ3W44o/lZ/hSwAACXNJREFUAa2GUeyHKAsJQBvDGMN9B/rLTqciWht5Ypu3JA4wHUggV+D4NUnQqwkET/rjGHGa6TNAlIU+FQTRgixVAy9lAskZQHft8gBYCgSHEln84Ez57qGTvkRdqomJ9oQEgCBaEFkA5D4+AHBiJgLGgIND3eg26zEruYC++MtpvPULx/Dkea/yWG8sjUQmj6lAAttJAIgKUBCYIFoQOXAfSIgWgC+WwRd+PoV79vTCJGgx5DApMYALkmXw/u+exe3jbhSKHK/456chaDVI54oUACYqQgJAEC2IRdBC0GkQiIstpv/xRxeRyRfx7vt2AxADzpelDJ9LizG4rQImfQl8/tkpuKwCfLEM5OQhSgElKkECQBAtCGMMbouAQCKLRCaPL//qKt5w8zZs94gX8yGHGU9d9CFXKGLSl8CbbxvFZW8cf/fDCxhymDHmtuAv7t2Nf3lyAvsG7Sr/b4hWhWIABNGiuKwGBOIZzEfSKHLgxm0O5b7dfTakc0U8fnYR2UIR4z1W/M1rDsCo12LCG8cbbh7Bvfv78O0/uh12Y31SU4n2gwSAIFoUp2QBeGNiIFhODQWAO8bFTKDPPH0FgDjxrcduxN+99hBuGnPidTcON3/BxJaDXEAE0aK4rAImvHH4pKZwPbalqt++LiN29VpxdDoEANjZI7qG7t7bi7up8RtRJWQBEESLIs4ZzihN4UotAGDJChjoMtatBQXRWZAAEESL4rQISOeKuOJPwqTXrho0c6dUEDbea1NjeUQbQAJAEC2KPLv53HwUPXbDqkEzN406YTXosG+AsnyIzUF2I0G0KPJsgQsLMewvk8ppErT4/tvvWHcIDUFUggSAIFoUeXxkKldQ5gSvZJi6fBI1QC4ggmhRSnf2PTTTgWgAJAAE0aK4LEsX/dIUUIKoFyQABNGimAQtzIIWANBrJwuAqD8kAATRwshuILIAiEZAAkAQLYzsBlpZBEYQ9YAEgCBaGLkWgILARCMgASCIFsZlFSDoNHUbNk8QpVAdAEG0ML990wj2D3atqgImiHpAAkAQLczhEQcOjzjWfyBBbAJyAREEQXQoJAAEQRAdCgkAQRBEh0ICQBAE0aGQABAEQXQoJAAEQRAdCgkAQRBEh0ICQBAE0aEwzrnaa6gaxpgPwPQm/9wNwF/H5dQLWtfGoHVtDFrXxmjVdQG1rW0b59yz8uCWEoBaYIwd5ZwfUXsdK6F1bQxa18agdW2MVl0X0Ji1kQuIIAiiQyEBIAiC6FA6SQA+qfYCKkDr2hi0ro1B69oYrbouoAFr65gYAEEQBLGcTrIACIIgiBJIAAiCIDqUjhAAxti9jLELjLEJxti7VVzHMGPsScbYOcbYGcbY26Xjf80Ym2WMvSD93K/C2qYYY6ek1z8qHXMyxh5njF2S/m3qZBLG2HUl5+QFxliUMfYONc4XY+wzjDEvY+x0ybGK54cx9h7p83aBMfbyJq/rbxlj5xljJxljjzLGuqXjo4yxVMl5+0ST11XxfVP5fP13yZqmGGMvSMebeb4qXRsa+xnjnLf1DwAtgMsAtgMQAJwAsFeltfQDuEG6bQNwEcBeAH8N4M9UPk9TANwrjn0YwLul2+8G8CGV38cFANvUOF8A7gRwA4DT650f6T09AcAAYEz6/GmbuK6XAdBJtz9Usq7R0sepcL7Kvm9qn68V9/89gP+twvmqdG1o6GesEyyAmwBMcM4nOedZAF8G8Eo1FsI5n+ecPy/djgE4B2BQjbVUySsBfF66/XkAr1JxLb8O4DLnfLOV4DXBOf8pgOCKw5XOzysBfJlznuGcXwEwAfFz2JR1cc5/yDnPS7/+AsBQI157o+taA1XPlwwTBy//FoAvNeK112KNa0NDP2OdIACDAK6V/D6DFrjoMsZGARwG8Evp0Nskk/0zzXa1SHAAP2SMHWOMPSId6+WczwPiBxRAjwrrknkIy7+Yap8voPL5aaXP3MMAvl/y+xhj7Dhj7CnG2B0qrKfc+9Yq5+sOAIuc80slx5p+vlZcGxr6GesEAWBljqma+8oYswL4GoB3cM6jAP4VwA4A1wOYh2iGNpvbOOc3ALgPwB8yxu5UYQ1lYYwJAB4E8FXpUCucr7Voic8cY+y9APIAvigdmgcwwjk/DOBPAfwXY8zexCVVet9a4nwB+G0s32Q0/XyVuTZUfGiZYxs+Z50gADMAhkt+HwIwp9JawBjTQ3yDv8g5/zoAcM4XOecFznkRwKfQIPN3LTjnc9K/XgCPSmtYZIz1S+vuB+Bt9rok7gPwPOd8UVqj6udLotL5Uf0zxxj7HQAPAHgDl5zGkrsgIN0+BtFvvKtZa1rjfWuF86UD8BoA/y0fa/b5KndtQIM/Y50gAM8BGGeMjUk7yYcAfEuNhUg+xk8DOMc5/4eS4/0lD3s1gNMr/7bB67IwxmzybYhBxNMQz9PvSA/7HQDfbOa6Sli2M1P7fJVQ6fx8C8BDjDEDY2wMwDiAXzVrUYyxewH8BYAHOefJkuMexphWur1dWtdkE9dV6X1T9XxJ3A3gPOd8Rj7QzPNV6dqARn/GmhHhVvsHwP0Qo+qXAbxXxXXcDtFMOwngBennfgBfAHBKOv4tAP1NXtd2iBkFJwCckc8RABeAJwBckv51qnDOzAACALpKjjX9fEEUoHkAOYi7r7esdX4AvFf6vF0AcF+T1zUB0T8sf8Y+IT32N6X39wSA5wH8RpPXVfF9U/N8Scc/B+B/rXhsM89XpWtDQz9j1AqCIAiiQ+kEFxBBEARRBhIAgiCIDoUEgCAIokMhASAIguhQSAAIgiA6FBIAgigDY8xV0gVyoaSLZZwx9nG110cQ9YDSQAliHRhjfw0gzjn/O7XXQhD1hCwAgtgAjLG7GGPfkW7/NWPs84yxH0p95F/DGPswE+cqPCaV9oMxdqPUTOwYY+wHKypiCUI1SAAIojZ2AHgFxPa8/wngSc75AQApAK+QROCjAF7LOb8RwGcA/F+1FksQpejUXgBBbHG+zznPMcZOQRxa85h0/BTEgSLXAdgP4HGx3Qu0EFsREITqkAAQRG1kAIBzXmSM5fhSUK0I8fvFAJzhnL9YrQUSRCXIBUQQjeUCAA9j7MWA2PKXMbZP5TURBAASAIJoKFwcQ/paAB9ijJ2A2OXxVnVXRRAilAZKEATRoZAFQBAE0aGQABAEQXQoJAAEQRAdCgkAQRBEh0ICQBAE0aGQABAEQXQoJAAEQRAdyv8PtAvaCMbM5GMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Simplest diffusion in one dimension (\"1D\") using a loop, and plot\n",
    "#seed = 987654\n",
    "#rng = np.random.default_rng(seed) # shortcut for function\n",
    "\n",
    "positions = []\n",
    "x = 0 # initial position taken to be zero\n",
    "positions.append(x)\n",
    "for t in range(200):\n",
    "    dx = 2.0*rng.random()-1.0 # random step in range (-1,1)\n",
    "    x = x + dx # new x position taken to be old position plus dx\n",
    "    positions.append(x) # this appends the *new* position\n",
    "    \n",
    "# print(positions)\n",
    "plt.plot(positions)\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Position\")\n",
    "\n",
    "# If you re-run this cell, it will produce a different diffusion trajectory every time.  Try it.  How could you prevent this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'y')"
      ]
     },
     "execution_count": 180,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZdrH8e+T3gvpCQkJNSShhy4ISLMidl2xi921vLru6rquZde1u6vr2rCLYqcJgiK9l0ASSAgppJLe68w87x8JMTEJTCBhJuH+XFcuMjNnztwnQ+aXc56mtNYIIYQQ5rCxdAFCCCF6DgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0hhBBms2hoKKUWKaXylVLxLe57SimVrZTa1/R1gSVrFEII8RtLn2l8CMxt5/5XtdYjm75WnuGahBBCdMCioaG13gAUW7IGIYQQ5rOzdAEduFcpdQOwC3hYa13y+w2UUguBhQCurq5jIiMjz3CJQgjRs+3evbtQa+3XmecoS08jopQKB5ZrrWOabgcAhYAGngGCtNa3nGgfsbGxeteuXd1cqRBC9C5Kqd1a69jOPMfSbRptaK2Paa2NWmsT8C4wztI1CSGEaGR1oaGUCmpxcz4Q39G2QgghziyLtmkopRYD0wBfpVQW8DdgmlJqJI2Xp9KBOyxWoBBCiFYsGhpa62vbufv9M16IEEIIs1jd5SkhhBDWS0JDCCGE2SQ0hDgJg9FEXlmtpcsQwipIaAhxEkt2ZTH5X7+w8kCupUsRwuIkNIQ4iUkDfDBpzX2L97J8f46lyxHCoiQ0hDiJcF9XZkcFYDRp/vjFPpbGnX3BkVlczSs/JWE0WXYGCWF51jr3lBAWobXmUF4FG5ILiMsqZXSYNxcND+b2Kf1ZnXAMo0nzwBd70Vozb2TIGa+toLKO1IKqpq9KUgurqKw1cPHIYC4bFYKrY/f8Sj/xfTzrkwuYHunPqDDvbnkN0TNIaIizXnFVPRsPF7AhuZCNhwvIr6gDINDDiZUH8nhu5UHGhfcBwMHWhpFhXjz45T5MWjN/VN9uq0trzTd7stmSUsiRwsaQqKg1ND/uaGdDhK8rWsNfv4/nXz8e4vLRISyY2I+B/u5dVsfmlELWJxcAsC21WELjLCehIc46BqOJvZmlbEguYENyAfuzy9AavFzsOWegL1MH+zF1kB+Bnk6kFlSyNC6n+ZJUvdHEjrRiHGxteGhJHEYTXDGm64Oj3mDiie8PsGRXFv7ujgz0d+PSkSH093Olv58b/X1dCfFyxsZGobVmz9FSPtmazuIdmXy0NYNJA3y4ZlwYYX1ccHO0xc3RHjcnO1zsbbGxUWbXYTJp/vnjQUK8nLG3VWxPK+KuaQO6/HhFzyGhIc4a+RW1PLfiIL8czKeizoCtjWJUqBcPzhzM1MF+DAvxxPZ3H6jBXs5EBrpzKNed1IKq5vvrjSYAHvk6Duja4CirbuDOT3ezNbWI+2cM5IGZgzv8oNdas+VIEW/9eoQLhgXxxEVRfLkzk8+3H+X+xXvbfY6box2ujrb0cXXkv38YTYSva4e1LNufQ3x2OS9fOYI9R0v4alcWCz/eRXSwJ9HBHsSEeBLg4YhS5geR6NkkNMRZYUdaMfd+vofy2gbmjwph6iA/Jg30xdPZvs22RpNmW2oRP+zL5sf4PCpqDfi6OXDTpHDmjQxmRF8vdh8tYem+HFYeyOXDLWldFhrphVXc8tFOsopreOWqEVw2+rf9aq3Jr6hj0+FC3lp/hJT8ylbP3ZRSyHXjw7hn+kDuPHcA+zJLKaupp6LWQFWdkcq6BirrjBRW1vH59qPU1BtxcbDtsJY6g5GXfkpiaJAHl44KYUSoF1tTi/gp8Rg/JR5r3s7H1YGopgCJDvZgbHgfAjycuuTnIayPxdfT6AqynoboiNaadzem8q9VSYT1ceGt60cTGejR7nbx2eV8vy+bZXE55FfU4eZox5zoQOaNDGbSAB/sbNt2NjQYTRi1xtGu4w9fc+1IK+aOT3ZR22DikTlD8HFzaGzwLqwirbCStIIqquqNbZ43OMCNu6cN5KLhQe3WCLAhuYDE3HIWTOjHX3+I59s92Xxw01imR/p3WM/7m9J4ZnkiH90yjnMHN67Tc/hYBbNe3dC8zf3nDSK3tIb4nHIOH6vAYNI429vy9V0TiQ72PM2fiOhup7KehoSG6LXKaxv4vyVx/JR4jPNjAnnhiuG4O7U9swD458qDvL0hFQBne1semDmIGyeF42R/+mFgjo2HC7j1w13Nl72OUwr6ejsT4dvYjrE/q5Q9R0txsLPhxz9OYYCf2wn3W28w8eLqQ7y7Ma3V/fefN4hbJ0eQXlRFZZ2BitoGKmoNVNQaqKwzUFRZx0dbMwCYPNCH0uoGSqsbKKmup7pFcL31h9GcP6xxNYM6g5HEnHLu+nQP9naKZfeeg5eLQ1f8eEQ3kdAQokliTjl3fbab7JIaHjs/klvPiTjhdffdGcUs2pzOtiNFFFXVA7Bwan/+csHQM1Lvqvg8Fm1KI8zHhQhfVwb4uRLh60Y/H5c2wbX1SBF3f7Ybk4b//mE0kwf6trvPo0XV3Ld4D3FZZTja2WBva0NlnaHdbU9kdJgX3i4OeLk44OVij5ezPS+vSQbgqzsnMrapZ9lxe4+WcMX/tuJib8sfZw5iTnQgoX1cOtz/gawy0ouquHhEcKdrE6dHQkMI4KtdmTzxfTxeLva8ed1oYn/3oXYiJpPmhkU72JRSyOvXjDzjYzHMdbSomts/3kVKQSV/vXAoN04KbxWKaYVVXPKfTVScJCS8Xex59tJhRPi64u5kh7uTHc4Otsx9bSNujnYsvXdyu2FbUdvAxf/ZRE2DkZX3T8HHzZH8ilo2JDd2z13WYgCkm6MdT10SzeWjQ1rtq6Sqnsv/t6W5g8HBp+fifII2FtH1JDTEWa22wchTSxP4Ymcmkwb48O9rR+Hr5tipfTz/4yH+t/4I/zd7MPfOGNRNlXaNyjoDD3yxj7UHj3HN2FCenheDg11jm8ax8lr+9eMhvt2b3bz9iFAv5kQHMCc6EH93R974JYVFm9NwtLPlvhkDuXlyRONlrwO53PXZHt68bjQXDg/q6OVJyClj/n+3MMjfrel2OQB+7o6cO9iP1Ql5rcaVAEQFeeBob8Peo6Wt7v/g5rFMH9Jx+4roHhIa4qxVUFHHwk92sfdoKfdMH8BDs4a06T57Mh9uTuOpZYksmNCPp+dF94hupCaT5uU1Sby57ghjw7156/ox+Lo5YjRp7l+8lxUHcnlo1mCuig0l0LNtj6bUgkqeW3GQnw/lE+HryhMXDuXfPx+mrKaBnx+edtKf4Zc7j/LE9/GMCvXm3CF+nDvYj6ggD2xsFHUGI098F0/ysQrisso63EfSs3O7pCOB6DwJDXFWSsqr4JYPd1JUVcerV41sbpjtjBX7c7l38R5mDQ3grevHNH9Yaq3JKqkhIaeccRF96ONqnQ27P+zL5tGv9+Pr5si7N8TyybYMFu84yuMXDOX2qf1P+vxfk/J5enli86Wi5+bHcO5gP+76dA9PXhzVpt2iJZNJmzVg8EBWGRe/san59gXDAnnzutE9Ipx7KwkNcVYpra5n8Y5M3lyXgouDLe/dGMvwvl6d3k9aYRVzXt3A8L6evH/jWJLzK9idUcKejBL2HC2lsLJxWpHHzo/kznOtdzT0/qxSbvtoV/M0KHdPG8CjcyPNfn6D0cRHW9LZm1nKy1eOYPn+XP7vqzjcHe34/PYJDOt76l1od2eUcPlbW5pvL7v3nNPan+gapxIaMrhP9Dgp+RV8sDmdb/ZkUdtgYsogX164YjhBns6ntL/Pt2dg0pr//mE0F/x7I9mlNQCE+7gwdZAvzg62fLb9KP1PMHLa0irrDKw4kNscGAB/nNm5Nhl7Wxtum9J4VpJeWMUbvxwGwNPFngWLtvPlwokMCez8nFZ/+e4An28/CsCwEE9+uGdyp6YyEdZFpkYXPYLJpPk1KZ8bFu1g5isb+Gp3FvNGhLDqgSl8cuv4Uw6MeoOJb/ZkM3NoAP4eTtw+JQIHWxu8XOx57PxIXrl6ZPM0GyPDOn8W09201ny/N5vzXv6Vt9enMjsqgOFNf8Gf/9pG1iXld2p/9QYTr689zOxXN5BeVE1MiAef3TYeB1sbrn9/O2mFVSffye8cX7zqo1vGsey+cyQweji5PCWsWnW9gW/2ZPPh5jSOFFTh7+7IDRP7ce24MHw62TOqPSv253LP53v48OaxTGvqvZN8rIKHluwjPrucS0cGU1RVT2pBFZsfm3Har9eejKIq/rnyEKF9nHn8wiizn5eQU8ZTSxPYmV7C8L6ePHVJNKObZqBddyifZ5YnklpYxYxIf/56UdQJ55iCxrEqj31zgMP5lVw8Ipg/jA8jxMuZ0D4uHD5WwdXvbMPZ3pYld04kxOvkIZ2QU0ZfLxfcnewkKKyUtGmIXuWDzWm8uiaZ8loDw/t6cus5EZwfE9TcrbQrLHh/O6kFVWx4dHqrnkINRhNvrkvhjV9SMJg07o52HPj7nC57XWjsIvzWr0d4a/0R6g0mfN0c2fXEzJM+r7S6npd+SuLz7UfxcnHg0TlDuCo2tM0Hc73BxIdb0vj3zynUGYzcck4E980YhNvv1twor23gxVVJfLo9g2BPZ569NKbd6UXis8u49t1t+Lk5suqBqSd8H77bm8VDS+IYEuDOkjsn4tHBSHxhWdKmIXqNitoGnltxkNFh3vzp/CGMDvPu8l42mcXVbDxcyAMzB7XpWmpva8MDMwczoq8XN3+4k4o6A499s5/HLxza4VQknfHLoWM8tTSRo8XVXDwiGF83Bz7YnE55bUOHH7DFVfWsPJDLSz8lUV7TwA0Tw3lw5mA8Xdrf3sHOhoVTB3DpqBBeWJXE2+tT+XZPNo/NjWT+qBBsbBSrE/J48od4CirquHlSBA/PHtzhQk5RQR6E9XGhuKqeE70V3+/N5uElccQEe3Iwt5w7P9nNhzeP69KwF5YjoSGs0pYjRRhMmodmD2ZMv467e2qtMWk6PSYjIaeMJ39IQCm4Kja0w+0amuaCGhfRhyW7Mtl4uJAXrxzOpAHtT91xMofyyvnXj4dYl1TAAD9XPrttPJMH+rIqPo8PNqeTVlCFr7sjKfmVzV9H8itJKaikuGl6k3ERffj7JdEMDWo78WJ7/N2deOnKEfxhfBhPLUvk4a/ieKVpGpDs0hoiA915Z0EsI0JP3Gbz9e4sEnLKef2akdh3MDHiD/uyeWjJPsZH+PDejbGsTsjjoSVxPPp1HK9ePVK61/YCEhrCKm1ILsDVwbb5Gv3vNRhN/LAvh7d+TSGrpIbRYd5M6O/DhP59GBnm1eFgsfyKWl75KZkvd2Xi5WzPi1eMIPgE1+dTmxp+F900lqS8Ch5eso/r3t3OTZPCuWJMX77fm83gQHdG9PVioL9bh+GVXVrDKz8l8+3eLNwc7fjz+ZHNI7ChcVJCgHlvbm71PC8Xewb6uTE7KoCB/m5EBXswsb/PKX34jgrz5ru7JvF/X8fx7Z7GkeKDA9xYdt85HYbAcRW1DbywOokx/by55HdzRGUWV7PyQC4rDuSyP6uMCf37MD3Sj9HPrKHO0Bi63+/L4ft9J19bfUakP28vGHPSeoTlSGgIq6O1Zn1yAZMG+ra5pFHbYGTJrkzeXp/a/FfytePC2JVRzGs/J6PXNi6DejxEJg7wYUSoJ452tizZlcnflyZQZzBx6+TG6/sdXdoBeGVNMv/+ubHbaczfVrd67MMt6Xy4Jb3Vfc72tsSEeDAsxIsRoZ4MC/HEy8WBt35NaZ4x9vYp/bl72oDm2V8NRhPf7m0cmHfcs5fGMNDfjYH+bvi4OnTpX+c2Noq50YHNoTHQ382sD+g31x2hsLKO92+MRSlFZnE1P8bnsmJ/bvNo72EhnsyNDiSvvJZ/rDzElEG+BHs68+WuTLPr++VQPg1Gk4SGFbNoQ7hSahFwEZCvtY5puq8P8CUQDqQDV2mtS060H2kI711SCyqZ8fJ6nrk0hgUT+gGNf+l+uu0o729KpbCynjH9vLl3+kCmDfFr/lAtq25gR3ox21KL2JZaRGJuOVo3hsiYft5sOVIEwC8Pn0v/k0wpDm0HpB3n5mjHyFAvNqUUcsPEftwwMZwD2aXEZZZxILuMhJwyahtaT3F+xZi+PDhrcHOvI5NJ882eLB5pERYA790Qy8yogM7/0E6ByaRRipOGUkZRFbNe2UBsuDfThvix4kAecZmNc0dFBXlgNGmSjlW0eZ6box39fFya56QCsFHw9oJYZrU4RqNJY2NGHaLr9bjeU0qpqUAl8HGL0HgBKNZaP6+Uegzw1lr/6UT7kdDoXT7YnMbflyWy4ZHpuDraNv9VX1FrYOpgP+6ZNoBxEX1O+iHTMkQ2HS5s/mBzcbDlwZmDmRtz4im7oXHK8js/3c1TF0eRX1HH+uSCVh+CK+4/p9ViQ0WVdYx5dm2rfbTszqu15pNtGTz5Q0KrbZ6bH8P8USG4OFjfyf8dn+xidcJvK/VFB3vg5+7Ir0kF7W4/LqIPs4YGkFlSzdHiajKKqluN7wjwcGT7X07eS0x0vx7Xe0prvUEpFf67u+cB05q+/wj4FThhaIjeJSGnHHcnOz7YksbiHUepM5iYGx3I3dMGdmrqCU8Xe2ZFBTArKgCjSfPLoXxu/3gX1fVGnlt5kOdWHmSQvxsjQ70YEerFyFAvhgS6t7o0crx7alSwJzdF9OHRuZEUVNSx8XAB1fXG5sAwmTR3fba71YfroptimRH521/Ur65J5vWmy13HfXbbeCYNOLU2ijMlo6ia6GAPAjycGlc4zClvFZwttVyUqaWU/EpmvrIeXzdH3rthbHeXLLqR9f1ZAwFa61wArXWuUkrmSz7LBHs6UVFr4OOtGVw6MoS7pvVnoH/np69o6e/LEvi4qV0B4PnLhlFW08C21CJ+PpTPV7uzgMZuqtHBHozo2xgixz/Lq1qsS+Hn7thq7W6Auz/b0xwY9583iAdnDmoVBEl5Fc2B4ePqwLd3T6Kfj/VOS9LSqgemYjRpxv9jLRW1BmZE+nPh8CA+3JzOrozfrhyPCPXqcLLIgf5uPDJnCC+uTiK3rEbmnerBLD64r+lMY3mLy1OlWmuvFo+XaK3bdKFRSi0EFgKEhYWNycjI+P0moofKKqnmuz3ZXDoq5KSXj8xxtKiaGS//yuSBvmxPK6K2wcSQAHdWPzgV+G0m27isUuIyf2ubqGn4bVnTN64bxUXDO15ZbkdaMcvicvjrRVHtjkcwmjRrEvOYOMAXT+eeOdBtf1YpA/zcqKg1MOGfP7d53MHWhpsnh3PPjIHtjjVpMJq4+D+bKKmuZ81D58qAPyvQ49o0oN3QSAKmNZ1lBAG/aq2HnGgf0qYhTuThJXEs35/Dxken42Bnw8in1wCQ/vyFHT7HYDSRUlBJXGYpaYXV3DYlotMLOvU2Wmvu/mwPP8bnNd/X19uZj28Zh4uDHS/9lMQ3e7Lo4+LAg7MGc83YUOxaXOorra7n6WWJfLs3mzvPHcBj55s/A6/oHj2uTaMDS4Ebgeeb/v3BsuWIniwlv5Lv9mZx6zkR+Hu0XoRIa91hW4KdrQ2RgR5EBpo3gK6301oT8eeVzbf/fkk0N04Kb7XNS1eO4KZJ4TyzPJEnvo/n463pPH5hFOcO9gNgya5Mvt2bTaCHE5GnMFuusA4W7QytlFoMbAWGKKWylFK30hgWs5RSh4FZTbeFOCWvrU3Gyd62eR2MmnojH9w8llsmR1h147M1umREMIP83Uh+9vw2gXFcTIgnXyycwNsLxlBvMHHjoh3cuGhH44SHsWEEezrh7GDL7Ogz061YdD2LX57qCnJ5SrRUVt3Asv05fLMnq3n510fmRFJaXc9lb20hxMuZT24db+kye716Q+Okj6//fBhHOxt2/GUmibnlXPfeNq4dF8Y/5g+zdIlnvd5yeUqIU1JSVc9ff4jnp8Rj1BsaG7sfv2AoN04Kp95g4s5Pd5NaUMXRomoq6wxtZnsVXWtrahFLmkaDTx/iz870YtKLqgj2dObz7Ue5ZEQwE/r7WLhK0VnyWyN6jdyyWtYdyqfeYEKpxnEa9UYTB3PL+WRbBttSi7l2XCiLd2Sy9UhRq1HJoms9tTSh1TQrqxLyWJXQ2IDu5WLP6DAv6T3VQ0loiF4jKtiDPU/OYnd6CRsOF7LxcAEvrk7ixdVJLbbxxNYmi42HCyQ0uonWmkN55QwL8STc15UIHxci/FwJ93Elwte1ed4t0TNJm4bo1T7dlsET38cDjYPyCprW0G45TkOIs5W0aYizxu6MYnJKa7l4RMcD7vZnlfLsikRGhXmx+PYJONrZkHyskg3JBfi6y1+7QpwKCQ3Ro+w5WsKra5LZeLiQCF9XLhoe1G7X2ZzSGm79aBe+bo68syAWJ/vG9TWGBLozRMYICHHKJDREj7Avs5RX1ySzPrlxZtUBfq58ePO4DsdaPPJ1HLX1Rj67bTx+7mf3SG4hupKEhrB6f/nuAJ9vP9p8e1x4H965YcwJG1QTcsqZNyqY/r6uJxz5LYToHAkNYfUONK0MB3DhsCBevmpE8+WmlgxGE8nHKtmZXkxpdeOiTYt3ZKK1xtXRDjdHu+Z/3Z3suCo29IRtIkKItiQ0hNUoq27gqWUJrIrP4/yYQG6dEsEAPzcOZDeGxk2TwnnyoihsfrcO97qkfN78JYX4dlbMu3FiOG6OtlTUGaisNVBVb6Ci1kBaYRX3Ld5LoKcTY8P7nLFjFKKnk9AQVuHXpHz+9M1+CivrmTnUn1UJeXy7N7vVNg/PHtwmMKCx0bvlug4tLdqchruTHf7ujvi7O+Htak9pdQNZJTUA3PXpHlbef06byQxbMpl0u68rxNlIxmkIi6qsM/CPlQf5fPtRBvm78cpVIxnW15OMoirOffHXVtu6O9rx93nRTOjvg6uDHRtTClh3qID1yfkUVtYDjVN1xwR7MqyvZ/O4jPzyWvIr6sivqKOkqp4Qb2fsbBTrmpYrje3nzee3T2heB6OqzsC+zFJ2phezM72YvUdLGRLozhvXjW5e41uI3qBHrqfRFSQ0eqbtqUX839dxZJXUcPuU/jw0a3BzW8XWI0Vc++423BztsLNVlFY3tLsPLxd7pg7yY3qkH1MH+eFj5poXWSXVnPOvdYzp583ujBICPZy4cHgQu9KLic8px2jSKAVDAz0Y3teT5ftzcbCz4c3rRjNxgMyXJHoHCQ3RI9Q2GHlxdRKLNqcR6u3Cy1eNaNWuUF1vYO5rG1EKVv1xKs4Otmit2XKkiIeXxJFXXgtAhK8rb1w3qnmd7s4Kf2xFm/vGRfRhXHgfYsO9Gd3Pu3l+pJT8Su74ZBfpRdX85YKh3DI5XHpkiR7vVELDoutpiLNPakElF/57I+9vSuP68f348Y9T2jREv7AqiaPF1bxw+XCcHRrPPJRSTB7oy7a/nMf6R6Zx06RwjpXXMu+NzcRnl7X3Uu0qrGy8XGUwmtp9fHZUAPefN4hpQ/xbTag30N+N7++ZzHmR/jyzPJEHv9xHTb2x3X1orVkVn8u6pHyz6xKip5AzDXFGfbI1nb/+kMAgfzc+uXU8gZ6tG6B3pBVz9TtbuWFCP/4+L+aE+yqqrGPu6xvxd3fkh3smt1patL1tX1t7mM93HMVoav1//soxfVkal0OdoTFIHO1suP+8QdwzfWCb/ZhMmjfXpfDK2mSGBnrw9oIxrdYxz6+o5fHv4lmTeIw50QG8vaBTf8QJcUbJmYawetdP6Mezl8aQVVLD7FfX88O+33pI1RmMPPp1HH29nXl0bsfrR6cWVPLyT0mkF1Xz9CXRJOSU8+7GtHa3rW0w8t9fUzj3xV/bDQyAr3ZnNQcGwJh+3ry4Oon/rT/SZlsbG8V95w3i39eM4lBeOZe8sYm4zFK01iyNy2H2qxtYk3gMgLkxgWb/XIToKeRMQ3SpHWnFvLcxFXs7G5zsbHG0t8HbxZ57pg/ExeG3Ht5phVU8tGQfe4+WctHwIJ69NIZ1Sfk8+GUci26KZUZk22nLCyrq+PfPrc8WLh/dl/SiKuKzy1j1wFQifF2BxjOCZftzeGFVEtmlNcwc6s9Af/d2g6AlpeD3vxIzhwZQXW+gqt5IdZ2B6nojuWU1mDTY2SgenTuEvUdL+TE+j5GhXjja2bA/q4xdT8zEVRZ6ElZMZrkVFldW08DmlEKqfne9P8TLhevGhzXfjvB15as7JnL3Z3tYvj+X5ftzmx+ztbFhc0ohSoGtUtQaTLyyJpm4zFKUguvH9+O2KREs3pHJ+5tSaTA2fso/8lUcX905kZ3pJTy3IpG4rDKigz148crhTBrgy8R//gxAsKcTc2OCmB7px7iIPigUD365jxUHctsEBsDag8eI7eeNp7M9wZ5OuDjYEezlxPgIH/LKa/nHyoNU1hp47PxIbpjYjwn/+Jm5MYESGKJXkjMN0eWq6gwsjcvh8+1Hm0dzA3x79yRGhXo19zoyGE3EPLW6zSjuE7lxYuu2jtSCSp5ensivTWMujgv0cOKROUOYPyqkeWBeakElNkoR3nQ2ApCQU8YjX+0nMbe8w9ecEx3AM5fGoDUsi8vB2cEWXzdHPt2WwcbDhQzv68nLV45gUIA7q+JzufPTPXx0yzjOHexn9nEJYQnS5VZYnQNZZVz8xqbm25GB7vxhfBjb04pZvj+Xf8wfxooDOWxOKTrhftyd7LhmbCgPzx7SZt4prTVrD+Zz+8e//R8YH9GHUWHeDA1yJzLQg/5+rti3aCivMxh545cU3vr1CF4uDsyK8mfxjkwWTm0cL7LuUD6fbs84aV3HX2tokAeDA9xZvj+Hw/mVbH1sxgkb5oWwBhIawiqlFlQy4+X19PdzxdneloSc9v+qf+3qkfT1dia1sIrd6SVklVazYEI4c6IDzBoTUV7bwOfbj5KUV8HB3HKOFFQ2X7pysLVhoL8bkUHu1BlMrGi6HHbZ6BDumzGIq97eio+rA/++dhSPf3eAnentT0sCMDLUi+p6AzEhnigUKfkVHM6vpLrpktztUyJ4/MKozv6YhDjjJDSE1broPxuxUYof7pnM/qwyrnlnGzUNbWByhe0AABaoSURBVMc5HPnHBdh20TxP9QYTqYWVHMqt4FBeBXGZpWxNbX3msOqBKby+9jBrDx7j+gn9+GBzepv9XDGmL0v35dC3jzPf3TUZTxf7NtuYTJrs0hoyiqoZFeYl7RmiR5CGcGG15o0I4bmVB0kvqmZEqBcHn5lLYk45l765mfoWA+0O5ZWf8gjv33OwsyEy0IPIQA8AFm1KaxMac1/b2Px9y8C4YWI/rhsfRqCHE5f9dwuujrZ8cNPYdgMDGrvihvZxaTVmQ4jeSEJDnLaKpstCqQVVVNYZfvuqbfy3oraByjoDABsPFzR3i40K9iDp2bks2pzOM8sTAfg1qYDBAe6t2h+6ytVjQ/Fzd+RY0wSG72xIbbPNs5fGMH9UCK6OdtQbTCx4fztZJTV8dvt4+vm4trNXIc4uEhrilNUZjHy27ShvrEuhuKoef3dH3J3scHOyx93RDj83R1ybFjxyc7TDw9mO82OCWu1DKcUtk8ObQ+PF1Um8uDqJ1Q9M7fK1vF0d7ZoXXXq7g/EaPx88xrmD/XBxsOXP3x5ge1oxr109UtbcEKKJhIboNKNJ88O+bF5Zk0xWSQ2TBvjwp7mRjAj1OqX9KaW4Y2p/fjmUz+H8SgA+3prOc/OHdWHVrX2zJwuAq2L78tQl0TjY2vDhlnReXZPMzFfWM2WQL2sP5vPAzEFcOiqk2+oQoqeRhnBhNq01vxzK58XVSRzKqyA62IM/zY1kyiDfLpvxVWtNea0Bd0e7bl34KLWgkoKKOsb3bz3NeW5ZDU8vS+TH+DzmjQzmtatHymy2otfqVb2nlFLpQAVgBAwnOjAJje63O6OY5388xM70EsJ9XHh49hAuHBbUa1e0S8mvJNzHRcZaiF6tN/aemq61LrR0EWez5GMVvLAqibUHj+Hn7sgzl8ZwzdjQbmmotiYD/d0sXYIQVsnaQ0NY0NYjRfzhvW2YNPi5O/LXi6KY2N+n1weGEKJj1hwaGvhJKaWBt7XW77R8UCm1EFgIEBYW1s7Txeka4OfKtePC2JdZyqG8Cu5fvBeAEC9nRoZ6MSLUk5Gh3sSEeLSawVYI0XtZc5tGsNY6RynlD6wB7tNab2hvW2nT6H419UYScsrYl1na/JVVUgOAjYLBAe7Mjg7k+vFh+Hs4nWRvQghr0KsawltSSj0FVGqtX2rvcQkNyyisrCMus5S4zFJ2ppewLa0IW6W4YFgQN04KZ3SYl/Q8EsKK9ZqGcKWUK2Cjta5o+n428LSFyxK/4+vmyHlDAzhvaOOCSemFVXyyLYMlOzNZGpfDsBBPbpwUzryRwdIOIkQvYZVnGkqp/sB3TTftgM+11s91tL2caViXqjoD3+7N5qMt6aTkVxLh68qf5g5hTnSgnHkIYUV67eWpk5HQsE7HBwM+/+MhDudXEtvPm2vHhfHUsgQqag1EBroT5OlEg1FTbzAxwN+Vf1423NJlC3HWkNAQVslgNPHV7ixe/imZwsq6DreztVHsfmImXi4OZ7A6Ic5eEhrCqlXVGXhnQypvrT9CvcGEi4Nt88JFLYX2cSYm2JOYkKavYA983BwtULEQvZuEhugRjpXX8uqaZL7Ymdl8388Pn8uxslrissqIzykjPruMjKLq5seDPZ2IDvEkxMuZiQN8mBMdaInShehVJDREjxL+2Irm710cbLkqNpRbJkcQ5tO4kFFZTQMJOWVsTy3m9Z8Pt3ru9r+cR4CMBxHitPSaLrfi7GBro7jr3AFcMCyI9zal8tn2DD7ems6c6EBum9KfMf28mTTAFwdbG17/+TAxIR6cM9CPqYN9JTCEsBAJDWExNgqMWhMV7MErV43k0TmRfLQ1nc+2ZfBjfB6jw7xYOLU/E/v7ohTMjgrk/vMGWbpsIc5qEhrCYmyU4ueDx6htMOLn7oifmyPnRfpz7/SBfLUrk/c3p3Hnp3sI6+OCg60NcZmlli5ZiLOehIawmMvH9GXbkSKW7MykqkUvqg2PTOemyREsmBjO6oQ83t2YytHiag7lVViwWiEESGiIM8hgNJFXXkt2SQ35FXXcP2MQo8O8eeOXw1QVVRMV5MGDswY3N4Tb2jTOY3XBsCB2Z5RQUdtg4SMQQkhoiC6ltWZXRgmpBZVkl9SQVVJDVmkN2SU15JXXYjS17a0XFeTBOwvGMCsqoMNpRsb08+7u0oUQZpDQEF1qf1YZV/5vKwBKQZCHEyHezoyL6EOIlzMh3s709XbG1kZx3bvbAUjMLWdokIfMSyVEDyDjNESX0lqzeEcmr6xJorCynvmjQnh07hCCPJ3bbFvbYCTqyVUcP/lYOLU/f7lg6BmuWIizlwzuE1ajoraBt349wnub0rBRsHBKf+44dwCujm1Pbj/bnsHj38U330775wXtnnU0GE1sSilkxf5cNiQXMCPSn79dHI2zg223HosQvZWEhrA6mcXVvLA6iWVxOfi7OzJ/VAhRwR5EB3sS4euKrU1jOJRVNzDi6Z8AiPvbbOoNJnzdHDCaNNtSi1m+P4dVCXmUVjfg7mRHfz834jJLGRzgxhcLJ9LHVSY5FKKzJDSE1dqdUcJLq5PYnVFCvdEEgLO9LZFB7kQHexAV5El0sAdDAt15b2MqL/2U3GYfAR6OBHo6YzCaSMwtR2sI9HBiyR0Tm3tcCSHMJ6EhrF69wcSRgkoScspJyCkjMaecxNxyKmoNZu9jRF9PZkQGMCPSn+hgD2xspAFdiFMhc08Jq+dgZ8PQIA+GBnlwxZi+QGPjeWZxDVNfXNdm+zH9vPnstvHkltWSXlRFdLAH/u4y75QQliILNwuLU0oR5uPCZ7eNx0bBxSOCuWvaAKDxstaDX+4j3MeF6UP8JTCEsDAJDWE1Jg/05eHZQ1gWl0OghxNPXRwFwI/xefx9WWKHz6ttaLuQk6XV1Bv5encWaxOPWboUIbqUXJ4SVuWucwew92gJz65I5IuFE3njulHc+/lePtySjp+7I/dMH9i8bUZRFU/+kMD65ILmRZoaV/xr7J0V4OF4xgcMphZU8tn2o3y1K5PyWgPnxwQyMyrgjNYgRHeShnBhdcqqG7jojY00GDQr7j+H+Jxybly0A4BnLo3hqti+vL0+lTfXpWBno7hufBj5FXXEZ5eRWljF8f/Svm4ORDeFyIT+PkwZ5Nct9TYYTaxNPMan2zPYnFLUfP/5MYG8ctVIGUcirJb0nhK9Rnx2GZe9tYWx4d58fMt4dqQVc+272wDo7+tKamEVFw4L4q8XRRHo+Vs7R1WdgYO55STklBOfXUZ8TjmHj1VgMGkWTu3PY3Mju6y3ldaaT7Zl8Oa6FI6V1+HqYNs8W+890wfw8Kwh0rNLWDXpPSV6jZgQT56dF8Oj3+zn1TXJ3DCpHzEhHsRnl2PUmg9vHsu0If5tnufqaEdseB9iw/s031fbYOSfKw/yzoZUjhZV8+rVp//Xf22DkT9/e4Dv9mYzoX8fbjunP59uz6ChtJbnLx/GZaP7ntb+hbBWcqYhrNqfvt7Pl7sycbSzQWu4c9oA7p42ACf7zn/of7A5jaeXJzI8xJN3b4w9pZ5YeWW1/OG9bRwpqALg4VmDGRPuzV2f7sHWRvH2gjGMbRFYQlgzuTwlep06g5GPtqSTVljF7VP609/P7bT2tybxGPcv3ksfVwcenj2YuTGBuDiYd8Idn13GRf/Z1Oq+cB8X0ouqGeTvxvs3jpWR6aJHkdAQwgwHssq4/4u9pBVW4eZox4XDgrgiti+x/bw77G21NvEYt3382/+xF68YziNf72++/a/Lh3FVbKhM7y56FAkNIcyktWZHWjFf785ixYFcquuNhPu4cMWYvlw2ui/BXs7N2727MZV/rDwENM6XNTs6gJ8SjlHTYGSQf+OZz+H8SsaF9+HJi6OICfG02HEJ0Rm9KjSUUnOB1wFb4D2t9fMdbSuhIU5HVZ2BlQdy+Xp3FtvTilEKzhnoy/xRIaw8kMvag/mttnewtWHeyGBuOSeCoUEeGE2aL3dm8vJPSRRX13PF6L48MmcI/h4yel1Yt14TGkopWyAZmAVkATuBa7XW7Q4LltAQXSWjqIpv9mTz8dZ0Sqtbr0nu7WLP9RP6sWBiv3Yb0ctrG3hzXQofbErHzlZx97QB3Dal/yk12gtxJpxKaJx0GhGl1L1KqTO9QPM4IEVrnaq1rge+AOad4RrEWaifjytXjw0l4Heh8Nz8GLY8dh4Pzx7SYa8rDyd7/nz+UNY8NJWpg/x46adkznt5PT8flKlERO9hTreRQGCnUmoPsAhYrbv/9CQEyGxxOwsY33IDpdRCYCFAWFhYN5cjziaPf3eApGMVTBrgw21TIpg22L9Tg/T6+bjyvwVj2HqkiGeWJzYP+BOiNzDr8pRq7BIyG7gZiAWWAO9rrY90S1FKXQnM0Vrf1nR7ATBOa31fe9vL5SnRlVLyG0eQRwZ6nPa+jCaNjUJ6VQmr1G0jwrXWWimVB+QBBsAb+FoptUZr/WjnSz2pLCC0xe2+QE43vI4QbQz0d++yfdnKNCKilzGnTeN+pdRu4AVgMzBMa30XMAa4vJvq2gkMUkpFKKUcgGuApd30WkIIIcxkzpmGL3CZ1jqj5Z1aa5NS6qLuKEprbVBK3QusprHL7SKtdUJ3vJYQQgjznTQ0tNZPnuCxg11bTqt9rwRWdtf+hRBCdJ6s3CeEEMJsEhpCCCHMJqEhhBDCbBIaQgghzCahIYQQwmwSGkIIIcwmoSGEEMJsEhpCCCHMJqEhhBDCbBIaQgghzCahIYQQwmwSGkIIIcwmoSGEEMJsEhpCCCHMJqEhhBDCbBIaQgghzCahIYQQwmwSGkIIIcwmoSGEEMJsEhpCCCHMJqEhhBDCbBIaQgghzCahIYQQwmwSGkIIIcwmoSGEEMJsEhpCCCHMJqEhhBDCbFYXGkqpp5RS2UqpfU1fF1i6JiGEEI3sLF1AB17VWr9k6SKEEEK0ZnVnGkIIIayXtYbGvUqp/UqpRUop7/Y2UEotVErtUkrtKigoONP1CSHEWUlprc/8iyq1Fghs56HHgW1AIaCBZ4AgrfUtJ9pfbGys3rVrV5fXKYQQvZlSarfWOrYzz7FIm4bWeqY52yml3gWWd3M5QgghzGR1l6eUUkEtbs4H4i1VixBCiNassffUC0qpkTRenkoH7rBsOUIIIY6zutDQWi+wdA1CCCHaZ3WXp4QQQlgvCQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0hhBBmk9AQQghhNouEhlLqSqVUglLKpJSK/d1jf1ZKpSilkpRScyxRnxBCiPbZWeh144HLgLdb3qmUigKuAaKBYGCtUmqw1tp45ksUQgjxexY509BaH9RaJ7Xz0DzgC611ndY6DUgBxp3Z6oQQQnTE2to0QoDMFrezmu4TQghhBbrt8pRSai0Q2M5Dj2utf+joae3cpzvY/0JgIUBYWNgp1SiEEKJzui00tNYzT+FpWUBoi9t9gZwO9v8O8A5AbGxsu8EihBCia1nb5amlwDVKKUelVAQwCNhh4ZqEEEI0sVSX2/lKqSxgIrBCKbUaQGudACwBEoFVwD3Sc0oIIayHRbrcaq2/A77r4LHngOfObEVCCCHMYW2Xp4QQQlgxCQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0hhBBmk9AQQghhNgkNIYQQZrNIaCilrlRKJSilTEqp2Bb3hyulapRS+5q+/meJ+oQQQrTPzkKvGw9cBrzdzmNHtNYjz3A9QgghzGCR0NBaHwRQSlni5YUQQpwiS51pnEiEUmovUA48obXe2N5GSqmFwMKmm3VKqfgzVaAF+AKFli6iG8nx9Wy9+fh687EBDOnsE7otNJRSa4HAdh56XGv9QwdPywXCtNZFSqkxwPdKqWitdfnvN9RavwO80/Rau7TWsb/fpreQ4+vZ5Ph6rt58bNB4fJ19TreFhtZ65ik8pw6oa/p+t1LqCDAY6PSBCSGE6HpW1eVWKeWnlLJt+r4/MAhItWxVQgghjrNUl9v5SqksYCKwQim1uumhqcB+pVQc8DVwp9a62IxdvtNNpVoLOb6eTY6v5+rNxwancHxKa90dhQghhOiFrOrylBBCCOsmoSGEEMJsPTo0evt0JB0dX9Njf1ZKpSilkpRScyxVY1dRSj2llMpu8Z5dYOmaTpdSam7T+5OilHrM0vV0NaVUulLqQNP71eN7OCqlFiml8luO+VJK9VFKrVFKHW7619uSNZ6ODo6v0793PTo0+G06kg3tPHZEaz2y6evOM1xXV2n3+JRSUcA1QDQwF/jv8V5nPdyrLd6zlZYu5nQ0vR9vAucDUcC1Te9bbzO96f3qDWMZPqTx96mlx4CftdaDgJ+bbvdUH9L2+KCTv3c9OjS01ge11kmWrqO7nOD45gFfaK3rtNZpQAow7sxWJ05iHJCitU7VWtcDX9D4vgkrpbXeAPy+t+Y84KOm7z8CLj2jRXWhDo6v03p0aJxEhFJqr1JqvVJqiqWL6WIhQGaL21lN9/V09yql9jedRvfYywBNeut71JIGflJK7W6a1qc3CtBa5wI0/etv4Xq6Q6d+76w+NJRSa5VS8e18neivtuPTkYwCHgI+V0p5nJmKO+cUj6+9mR6tvu/0SY71LWAAMJLG9+9lixZ7+nrke9RJk7XWo2m8BHePUmqqpQsSndbp3ztrnLCwld4+HcmpHB+Nf7WGtrjdF8jpmoq6j7nHqpR6F1jezeV0tx75HnWG1jqn6d98pdR3NF6Sa699sSc7ppQK0lrnKqWCgHxLF9SVtNbHjn9v7u+d1Z9pnIqzYDqSpcA1SilHpVQEjce3w8I1nZamX8jj5tPYCaAn2wkMUkpFKKUcaOy4sNTCNXUZpZSrUsr9+PfAbHr+e9aepcCNTd/fCHQ02WqPdCq/d1Z/pnEiSqn5wH8APxqnI9mntZ5D43QkTyulDIAR86cjsSodHZ/WOkEptQRIBAzAPVproyVr7QIvKKVG0ngJJx24w7LlnB6ttUEpdS+wGrAFFmmtEyxcVlcKAL5TjWvi2AGfa61XWbak06OUWgxMA3ybpjn6G/A8sEQpdStwFLjSchWeng6Ob1pnf+9kGhEhhBBm65WXp4QQQnQPCQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGaT0BBCCGE2CQ0huoFSamzTJHBOTaOnE5RSMZauS4jTJYP7hOgmSqlnASfAGcjSWv/TwiUJcdokNIToJk1zTu0EaoFJvWCqFyHk8pQQ3agP4Aa403jGIUSPJ2caQnQTpdRSGlfsiwCCtNb3WrgkIU5bj57lVghrpZS6ATBorT9vmqZ/i1Jqhtb6F0vXJsTpkDMNIYQQZpM2DSGEEGaT0BBCCGE2CQ0hhBBmk9AQQghhNgkNIYQQZpPQEEIIYTYJDSGEEGb7fwvtaCO0+rc4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Diffusion in two-dimensions (2D)\n",
    "xposns = [] # can write out as xpositions\n",
    "yposns = []\n",
    "x = 0 # initial position taken to be zero\n",
    "xposns.append(x)\n",
    "y = 0 # initial position taken to be zero\n",
    "yposns.append(y)\n",
    "for t in range(200):\n",
    "    dx = 2.0*rng.random()-1.0 # random step in range (-1,1)\n",
    "    x = x + dx # new x position taken to be old position plus dx\n",
    "    xposns.append(x) # this appends the *new* position\n",
    "    dy = 2.0*rng.random()-1.0 # random step in range (-1,1)\n",
    "    y = y + dy # new y position taken to be old position plus dx\n",
    "    yposns.append(y) # this appends the *new* position\n",
    "    \n",
    "# print(xposns)\n",
    "limit = 15\n",
    "plt.plot(xposns,yposns)\n",
    "plt.xlim([-limit, limit])\n",
    "plt.ylim([-limit, limit])\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interactively\n",
    "# write a program to determine the number of steps for a 1D particle to exceed a certain distance from the start\n",
    "# Same, but for 2D\n",
    "# What is the average distance at fixed time?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare to a YouTube video of diffusion of food color or bacteria.  What else is going on?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Play around more, learn more ...\n",
    "# Demonstrate that using the same seed always leads to the same sequence of random numbers\n",
    "# How would you transform a (0,1) random number to a (-5,-5) random number?\n",
    "# How would you change the time step in 1D diffusion?  How could you change the matplotlib command to plot with the new time step?\n",
    "# Use functions to do 1D diffusion: One function for (-1,1) random number, another for simulation based on number of steps chosen by user, with plot.\n",
    "# A loop of loops - Multiple repeats in 1D and some interesting properties of diffusion.  Plot all runs on the same graph for visual comparison. Plot the average of x and the average of x^2 vs. time\n",
    "# Measuring the diffusion constant: Use D = <x^2> / 2 t, where <x^2> is the average of x^2, averaging over all trajectories.  What causes diffusive behavior in real systems? \n",
    "# Use random numbers to implement diffusion on a grid."
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}