brandon thermostate.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "1CzRnSUbYV2E"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from tqdm import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "Q0fw2bTbYV2G"
   },
   "outputs": [],
   "source": [
    "rcut = 3.0  ##### Cut-off distance. If distance between atoms is greater than this, interaction between them is ignored.\n",
    "rcutsquare = rcut**2\n",
    "g2 = 1./rcutsquare\n",
    "g6 = g2*g2*g2\n",
    "potcut = 4.*g6*(g6-1) ##### Term to be subtracted from pair potential energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "i51VhivYYV2G"
   },
   "outputs": [],
   "source": [
    "##### Function to compute acceleration for a pair of atoms. r12 is relative displacement of atoms.\n",
    "\n",
    "def acceleration(r12): \n",
    "    r12square = np.dot(r12,r12)\n",
    "    f2 = 1./r12square #1/r^2\n",
    "    f6 = f2*f2*f2 #1/r^6\n",
    "    acc = 48.*f2*f6*(f6-0.5)*r12\n",
    "    return acc\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "n8B7Y76QYV2H"
   },
   "outputs": [],
   "source": [
    "#### Function to compute potential energy for the system of particles. `pos' is the array of position vectors.\n",
    "def potentialenergy(pos):\n",
    "    potential = 0.\n",
    "    for i in range(Natoms-1):\n",
    "        for j in range(i+1,Natoms):\n",
    "            rij = pos[i] - pos[j] ## Relative position vector of the pair\n",
    "            for l in range(2):  ### Calculating the correct separation between particles because of periodic boundary conditions.\n",
    "                if abs(rij[l])>0.5*L: rij[l] -= L*np.sign(rij[l])\n",
    "            rijsquare = np.dot(rij,rij)\n",
    "            if rijsquare < rcutsquare: # Imposing interaction cut-off\n",
    "                f2 = 1./rijsquare\n",
    "                f6 = f2*f2*f2\n",
    "                potential += 4.*f6*(f6-1.) - potcut\n",
    "    return potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "id": "tT49TSmdYV2H"
   },
   "outputs": [],
   "source": [
    "number = 5\n",
    "Natoms = number**2  # No. of atoms\n",
    "Ratom = 0.5  ### Radius of atom used to illustrate the atom as a sphere.\n",
    "rho = 0.5 ### Number Density\n",
    "L = math.pow(Natoms/rho,0.5) # Length of a side of the square containing the gas.\n",
    "T0 = 120. # Natural temperature scale, T0 = epsilon/k.\n",
    "# Temp = 0.01 #dimensionless temperature\n",
    "T = 0.1   # Temperature in Kelvin\n",
    "dt = 1e-2  # Time step for Verlet Algorithm\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "id": "12feYni7YV2I"
   },
   "outputs": [],
   "source": [
    "###### Atoms are arranged initially in a square lattice configuration with spacing between neighbors equal to equilibrium separation of the potential energy.\n",
    "\n",
    "equilib_sep = math.pow(2.,1./6.)  ## Equilibrium separation of atoms r_0 = 2^{1/6}\n",
    "wall_spacing = (L-(number-1)*equilib_sep)/2.  ## Distance betwee the walls and the initial array of atoms.\n",
    "poslist = []  # List for positions of atoms\n",
    "vlist = []   # List for velocities\n",
    "\n",
    "################  Initialize positions and velocities ##############################################\n",
    "x, y = wall_spacing, wall_spacing\n",
    "\n",
    "for i in range(number):\n",
    "    for j in range(number):\n",
    "        xi, yi = x + equilib_sep*j, y + equilib_sep*i\n",
    "        poslist.append((xi,yi))\n",
    "\n",
    "\n",
    "####### Initially, all particles have the same speed but random velocity directions. ################\n",
    "v0 = math.sqrt(2.*T/T0)\n",
    "\n",
    "for i in range(Natoms):\n",
    "    phi = 2*np.pi*np.random.random()\n",
    "    vx = v0*np.cos(phi)\n",
    "    vy = v0*np.sin(phi)\n",
    "    vlist.append((vx,vy))\n",
    "###################################################################################################\n",
    "\n",
    "pos = np.array(poslist)  ### Converts lists to numpy arrays\n",
    "v = np.array(vlist) ## Scaled velocity\n",
    "\n",
    "\n",
    "V_cm = np.sum(v, axis = 0)/Natoms ## Correcting for CM velocity.\n",
    "V = np.array([V_cm,]*Natoms)\n",
    "v -= V  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "dWLHWFPEYV2J",
    "outputId": "26f61300-d63b-4425-9b3d-beba394105e6"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████| 5000/5000 [00:37<00:00, 134.55it/s]\n"
     ]
    }
   ],
   "source": [
    "time = 0. # Initial time.\n",
    "t_final = 50. # Time upto which simulation is carried out.\n",
    "\n",
    "potential_energy = potentialenergy(pos)\n",
    "kinetic_energy = 0.5*sum(np.square(v).sum(axis=1))\n",
    "energy = kinetic_energy + potential_energy\n",
    "Time_List = [time]\n",
    "Energy_List = [energy]\n",
    "PotentialEnergy_List = [potential_energy]\n",
    "KineticEnergy_List = [kinetic_energy]\n",
    "speed_list = [] #speed list for plotting histogram\n",
    "\n",
    "# Temp = 0.1 #temperatue that we can vary\n",
    "K_mean = T\n",
    "\n",
    "iterations = int((t_final-time)/dt) # Number of iterations of the Verlet Algorithm.\n",
    "\n",
    "##################################### The Verlet Algorithm Loop ###################################\n",
    "for step in tqdm(range(iterations)):\n",
    "    pos += v*dt/2. ## Updating positions.\n",
    "    for i in range(Natoms): ## Imposing periodic boundary conditions.\n",
    "        for j in range(2):\n",
    "            if pos[i][j] > L: pos[i][j] -= L\n",
    "            elif pos[i][j] < 0: pos[i][j] += L\n",
    "    accel = np.zeros((Natoms,2)) ## Initialising acceleration array.\n",
    "    for i in range(Natoms-1):\n",
    "        for j in range(i+1,Natoms):\n",
    "            rij = pos[i] - pos[j]\n",
    "            for l in range(2):  ### Calculating the correct separation between particles.\n",
    "                if abs(rij[l])>0.5*L: rij[l] -= L*np.sign(rij[l])\n",
    "            if np.dot(rij,rij) < rcutsquare:\n",
    "                acc = acceleration(rij) # Computing acceleration for a pair of particles.\n",
    "                accel[i] += acc\n",
    "                accel[j] -= acc\n",
    "    v += accel*dt  ## Updating velocities.\n",
    "    pos += v*dt/2. ## Final updating of positions.\n",
    "    for i in range(Natoms): ## Imposing periodic boundary conditions.\n",
    "        for j in range(2):\n",
    "            if pos[i][j] > L: pos[i][j] -= L\n",
    "            elif pos[i][j] < 0: pos[i][j] += L\n",
    "                \n",
    "    potential_energy = potentialenergy(pos)\n",
    "    kinetic_energy = 0.5*sum(np.square(v).sum(axis=1))\n",
    "    \n",
    "    #scaling the velocity for keeping the temperature constant\n",
    "    Km = np.mean(kinetic_energy)/Natoms\n",
    "    lamda = np.sqrt(K_mean/Km)\n",
    "    v *= lamda\n",
    "    energy = kinetic_energy + potential_energy\n",
    "    Energy_List.append(energy)\n",
    "    PotentialEnergy_List.append(potential_energy)\n",
    "    KineticEnergy_List.append(kinetic_energy)\n",
    "    time += dt\n",
    "    Time_List.append(time)\n",
    "    if time>10:\n",
    "        for atom in range(Natoms):\n",
    "            speeds = np.sqrt(np.dot(v[atom],v[atom]))\n",
    "        speed_list.append(speeds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARkAAAD8CAYAAABD7tCSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAATYElEQVR4nO3de3BV5b3G8edNIHKpNy6CTbQgFwteENmWaulUtIyKglYRZcp4Wi9oBbHTi4WhU0dqWztOeyyUjgWUUdpK8QpKbSmi1qroSaqoR0eFigNoUAGxXEwI+Z0/3p3DTtYGdmLedyWb72dmjdlv1l7rSXb2w7okvs7MBAChlKQdAEBxo2QABEXJAAiKkgEQFCUDIChKBkBQBywZ59zxzrmXc5ZPnHPfbbKOc87Ncs6tcc694pw7NTs+sslzP3XOXZT93JTs+uac6xHiiwOQPtec35NxzpVK2ihpuJm9mzM+WtINkkZLGi7pN2Y2vMlzu0laI6nCzHY654ZK2irpKUkZM/voM34tANqgDs1c/2xJa3MLJutCSfeab6xVzrkjnHNHm9n7OeuMk/S4me2UJDN7SZKccy2MDqA9aG7JXC7pvjzj5ZLW5zzekB3LLZnLJf26mfuTc26SpEmS1LVr12Ff/OIXm7sJAAWqqqr6yMx6tuY2Cy4Z51yZpLGSpjd3J865oyWdJOlvzX2umc2VNFeSMpmMVVZWNncTAArknGt6lvKZNefu0nmS/mVmm/J8bqOkY3IeV2THGoyX9LCZ7W5+RADtWXNKZoLynypJ0lJJV2TvMn1Z0rYm12P291wARaygknHOdZU0StJDOWPXOeeuyz78i6R/y989mifp+pz1+sgf5TzdZJtTnXMb5I96XnHOzW/5lwGgrWrWLey0cU0GCMs5V2VmmdbcJr/xCyAoSgZAUJQMgKAoGQBBUTIAgqJkAARFyQAIipIBEBQlAyAoSgZAUJQMgKAoGQBBUTIAgqJkAARFyQAIipIBEBQlAyAoSgZAUJQMgKAoGQBBUTIAgqJkAARFyQAIipIBEBQlAyAoSgZAUJQMgKAoGQBBUTIAgqJkAARFyQAIipIBEBQlAyAoSgZAUJQMgKAoGQBBUTIAgqJkAARFyQAIipIBEBQlAyAoSgbhbNniFxzUKJli8+mn0qJF0uzZ0tq16WTYskW66CKpZ0+pRw9pzBhp8+Z0siB1lEwxWb9eGjxYmjBBmjpVGjBAuvPO+Dmuu05askSqr5fMpMcekyZNip+jwXPPSZMn++9JVVU6GerqpB//WOrdWzrySOk735G2b08nS2xm1m6WYcOGWZtSX292xx1mJ51kNmiQ2U9/arZ7d3p5vv1tM/+23rt07my2dWu8DLt2mZWWJnOUlppt3x4vR4P58xvnKCkxW7w4fo7p05Pfk/Hj4+c4AEmV1srv29SLozlLmyuZW25J/uBcf316eQYOTOaRzFasiJehpsbskEOSGcrKzHbujJfDzBd+r17JLP37x81hZta9ezJHSYnZRx/Fz7IfIUqG06XPYtas5Nhdd0k7dsTPIvnTo6ZKSqR+/eJlKCuTJk5Mjk+YIHXuHC+H5K8NbdqUHF+zRqqtjZtl587kWH29VFMTN0cKKJmWqq+XPv44OV5Tk17JzJghderUeOzKK6U+feLmmDXLX5fp0sUXy6RJ0pw5cTNI/qJz377J8SFDfBnGNH58cmzECOnzn4+bIw2tfWgUcmlzp0tjxiQPgYcPTzfTq6/6U7ZLLzW75x6zPXvSy1Jf75c0LV3qT9UaXp8uXcxWroyfY+tW//PinM9xxhlm774bP8cBKMDpkvPbbR8ymYxVVlamHWOvDRuksWOll17yjwcOlJYulY4/Pt1caGzdOun++6XSUumyy6Ty8vSyVFf7U7Vjj00vw34456rMLNOq26RkWsFLL/lblJmM5FzaaYAWC1EyHVpzYwetoUPTTgC0WVz4BRAUJQMgKEoGQFCUDICgKBkAQVEyAIKiZAAERckACIqSARAUJQMgKEoGQFCUDICgKBkAQVEyAIKiZAAERckACIqSARAUJQMgKEoGQFCUDICgKBkk1dUdFDMbIg5KpsGHH0r33SetWOFnh4xt3To/6+Lpp0uTJ0vr18fPsHu3dOON0uGH+9kfv/GN/NO8xrJ5s3TrrX6upNtvlz75JL0sEsXbUq09W1zIJdgMkn/8Y+NJ4k85xezDD8PsK58PPjDr3bvxTJTl5WZbtsTLYGY2bVpyRsyRI+NmaLB1q1m/fo2zDBlitmtX/CyLF+/NcvLJZitWxM8QiQLMIJl6cTRnCVIy27aZde2afHPdcEPr72tfbr89uX/J7De/iZfBLFl0DUsa06n+6lf5s9xzT9wcVVVmJSWNM3TubLZ+fdwckYQoGU6XXnhB2rEjOb5iRbwM+zo12rAhXgbJv4XySeP08c0384+/9VbcHPfem/z6d+2SFi2Km6Mdo2SOOSb/eMy5ikeNat54KBMnJsdGjJD69ImbQ5K+8pX842ecETfHnj35x9Mo3vaqtQ+NQi7BrslcfHHjw+EOHcxWrgyzr3zq682uvXbv/p2Le7rW4NNPzSZN8tennDM791yzjRvj5zAzq601O+ecxq/L+PH+exXT88/770VujrIys3feiZsjEgU4XXK2r0PkNiiTyVhlZWXrb7i2Vpo9W1q2TDrqKGnq1Pj/YkrSG29Ir74qnXKKNHBg/P032LXL32k67LD0Mkj+Lb18uf+eZDLSmWemk2PBAmnGDOn996XjjpPuuEMaMyadLIE556rMLNOq26RkgALU1Ulbt0rdu0slxXuVIUTJdGjNjQFFq0MHqWfPtFO0S8VbyQDaBEoGQFCUDICgKBkAQVEyAIKiZAAERckACIqSARAUJQMgKEoGQFCUDICgKBkAQVEyAIKiZAAERckACIqSARAUJQMgKEoGQFCUDICgKJnWsm2btHFj2inQnmzbJs2cKZ1zjjRlirR2bdqJgmi/JbNkiTR0qJ+24/zz/XQiadi9W7r2Wj+VSkWFNGyY9Prr6WSRpMWL/Q/tmWdKv//9vmeFRLrq6qSzz5ZuvtlP+zJnjvSlL0nr1qWdrPW19kROIZf/n9zt2WeT8xP37m22fftnmNaqhW65xRLzNQ8YEH8SMjOzOXOSWW66KX6OBtu2mc2bZ3bbbWavvZZejrbokUeSr5Vk9oMfpBpLzIWdNXducprQ6mrpkUfiZ7nvvuTY229LacwPddttybHf/lbavj1+lrVrpUGDpGuukaZNk048UZo1K36OBp984r8/F1/sJ2qrrk4viyS9807zxtux9jnv0s6d+cd37IibQ5LKypo3HtJ77yXHdu70k5J97nNxs9x8czLPtGl+vu1u3eJmqa2VvvY16eWX/eOHH5YWLpSqqtKbS2nkyOaNt2Pt80jm0kuTY4ccIo0dGz/L1VcnxzIZaciQ+FlGjUqODRokHXNM/CyrViXHdu2SVq+On+XBB/cWTIP166V58+JnaTBkiDR9euOxUaOkq65KJ09A7bdkZsyQOnf2j3v3lhYt8v+NbcoU6Wc/k3r1kjp29IfjS5bEzyH505G+ffc+7tFDuvvudLIMGpQcKy2VBgyIn2XNmvzjb78dN0dTP/+5v2Exf7709NP+AnCnTulmCqB9z4XdcNt4wAD/Boe/a7FypVRTI33963uLOLYXX/SH/rmntjfe6Cerj+3JJ6WzzkqO33WXdOWV8fO0YSHmwm7fJYO27e23/UX6jz6SxozxR3lpufJKacGCvY9HjZIeeyyda2dtGCVDyeCzeOEFf63ohBP876g4l3aiNidEybTPu0tASwwf7hdE1T4v/AJoNygZAEFRMgCComQABEXJAAiKkgEQFCUDIChKBkBQlAyAoCgZAEFRMgCComQABEXJAAiKkgEQFCUDIChKBkBQlAyAoCgZAEFRMkB7t2mT9Oyz0scfp50kL0oGjT3zjJ/a9rnn0k7Sdpn5N/WyZenMWppr+nQ/ed+IEVJ5uTR7drp58mntybVDLsOGDftMk4m3Wbt3m61ebVZdnV6GujqzSy5pPPn7ZZeZ7dmTXiYzs5oas5kzzU4+2eyMM8z+8Id082zaZHbqqXu/R0ccYbZ8eTpZli1r/Ho1LK+80uJNSqq0Vn7fpl4czVmKsmT+/nez8nL/UpSWml19tS+d2B54IP8P7JIl8bPkmjgxmWnevPTyXHVVMk+vXma1tfGzXHNN/tfs1ltbvMkQJcPpUpr+8x9p3Dg/C6Yk7dnjpyydMyd+lmeeyT/+j3/EzZGrulr605+S47/+dfwsDf761+TYpk3JubZj6NYt//iRR8bNcQCUTJqeeMJPtdvUgw/Gz9KvX/7x/v3j5si1ebNUX58c/+CD+FkaHH10csw5Pxd6bFdfLXXp0nisZ09pwoT4WfaDkknToYc2bzykK65IFs3AgdI3vxk/S4NBg6S+fZPj558fP0uDH/4wOXbZZdKxx8bP0r+/n/f83HP9/seN80eebexIhmlq01RfL514ovTGG43HH31UuuCC+Hk+/NDfnVi9Who6VLrhBql79/g5cj3/vJ9Du7raP85k/F2do45KL9Ojj/o7cB9/LF10kfT97xfNnNrMhV1sJSNJGzZIN90kLV8uVVRIP/pRmzvcTV1trb9lfOihvmQQDHNhF6OKivwXN7FXWZk0cmTaKdBCXJMBEBQlAyAoSgZAUJQMgKAoGQBBUTIAgqJkAARFyQAIipIBEBQlAyAoSgZAUJQMgKAoGQBBUTIAgqJkAARFyQAIipIBEBQlAyAoSgZAUJQMwvjgA+mBB6RVq9JOgpRRMmh98+f7SeAvvVQ6/XTprLPSn5i+wXvvSevXp53ioELJFKvaWmnr1vj73bRJuv56v/8GTz6Z7tSykrRli3TeeVJ5uZ8I7atf9dPRpOX116XHH0/nNYqMkilGM2f6yc+6dZNOOy3uPM1PPSXt3p0cX748XoZ8Jk9uPI/1P/8pfetb8XPU1kqXXCKdcII0erQvvQUL4ueIiJIpNgsWSDffvHeO7cpK/8NcUxNn/+Xl+ccrKuLsP5/6+vzziz/xhJ9vO6bf/U566KG9j3ftkiZN8qdxRYqSKTYLFybH3n9fWrEizv5HjPCnIrnKyqTvfS/O/vNxTurUKTneoYPUsWPcLH/5S3Ksri7e65MCSqbYlOzjJS0tjZdh2TLpJz/xF33Hj5eeecaftqXFOemaa5LjEyZIhx0WN0uvXs0bLwZm1m6WYcOGGQ5g4UIzqfFSUWFWW5t2snTV1ppNm2Z21FFmRx5pNmWK2Y4d8XO8+KJZx46NX5+TTzbbsyd+ljwkVVorv285kik2EydKv/yl1KOHfzxihL/gGfu0oK3p2FH6xS/83a8tW6TZs6UuXeLnOO00f7ftggukE0+Upk7114b2dQRaBJwvr/Yhk8lYZWVl2jHah/p6f7G3c+e0k6Adcc5VmVmmNbdZvPV5sCspoWDQJlAyAIKiZAAERckACIqSARAUJQMgKEoGQFCUDICgKBkAQVEyAIKiZAAERckACIqSARAUJQMgKEoGQFCUDICgKBkAQVEyAIKiZAAERckACIqSARAUJQMgKEoGQFCUDICgKBkAQVEyAIKiZAAERckACIqSARAUJQMgKEoGQFCUDICgKBkAQVEyAIKiZAAERckACIqSARAUJQMgKEoGQFCUDICgCioZ59y5zrk3nXNrnHPT8nz+EOfcn7Off8E51yfnc9Oz4286584pdJsAisMBS8Y5VyppjqTzJA2WNME5N7jJaldJ2mpm/SX9t6RfZp87WNLlkk6QdK6k3znnSgvcJoAiUMiRzJckrTGzf5tZraRFki5sss6Fku7JfvyApLOdcy47vsjMaszsHUlrstsrZJsAikCHAtYpl7Q+5/EGScP3tY6Z1Tnntknqnh1f1eS55dmPD7RNSZJzbpKkSdmHNc651wrIDKBljm/tDRZSMqkys7mS5kqSc67SzDIpRwKKlnOusrW3Wcjp0kZJx+Q8rsiO5V3HOddB0uGSNu/nuYVsE0ARKKRk/kfSAOdcX+dcmfyF3KVN1lkq6b+yH4+TtNLMLDt+efbuU19JAyS9WOA2ARSBA54uZa+xTJH0N0mlku42s/91zs2UVGlmSyXdJWmhc26NpC3ypaHseoslvS6pTtJkM9sjSfm2WUDeuc3+CgE0R6u/x5w/4ACAMPiNXwBBUTIAgkqtZPhTBSAs59zxzrmXc5ZPnHPfbbKOc87Nyr5nXnHOnZodH9nkuZ865y7Kfm5Kdn1zzvU4YBAzi77IX+xdK+k4SWWSVksa3GSd6yXdmf34ckl/zn48OLv+IZL6ZrdTWsg2WVgO1iX7/qiW9IUm46MlPS7JSfqypBfyPLeb/A2dLtnHQyX1kbROUo8D7TutIxn+VAGI62xJa83s3SbjF0q617xVko5wzh3dZJ1xkh43s52SZGYvmdm6QnecVsnk+1OF8n2tY2Z1knL/VCHfcwvZJnCwulzSfXnGC3nf7Ou5BeHCL1Dksr/wOlbS/S147tGSTpL/nbYWSatk+FMFIJ7zJP3LzDbl+dyB3jfjJT1sZrtbuvO0SoY/VQDimaB9n+4slXRF9i7TlyVtM7P3C3xuYVK82j1a0lvyd4RmZMdmShqb/biT/OHdGvkSOS7nuTOyz3tT0nn72yYLy8G8SOoqfwZweM7YdZKuy37s5P8HcmslvSopk7NeH/mjmpIm25wqf+2mTtJ7kubvLwN/VgAgKC78AgiKkgEQFCUDIChKBkBQlAyAoCgZAEFRMgCC+j9Gs2HBTIjlfwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def create_picture(positions):\n",
    "    plt.cla()\n",
    "    plt.gca().set_aspect('equal')\n",
    "    plt.axis([0, L, 0, L])\n",
    "    plt.setp(plt.gca(), xticks=[0, L], yticks=[0, L])\n",
    "    for x,y in positions:\n",
    "        atom = plt.Circle((x, y), 0.1, fc='r')\n",
    "        plt.gca().add_patch(atom)\n",
    "create_picture(pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 1\n",
    "K_mean = np.mean(np.array(KineticEnergy_List))/Natoms\n",
    "T = K_mean\n",
    "\n",
    "def prob_dist(v):\n",
    "    return (v/T)*(np.exp(-v**2/(2*T)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x28ccb566b30>]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnVklEQVR4nO3dd3hU153/8fd3ZlRRQUISCHV6x4DAVBfsZDGOIQlu2LjFcYnX2WST3Y1/m8S/PE5+qbubJ4kdl9iOa+y4xAmOsbHXJTYY2YgOEkUIgSoqqHdpzu8PDURWBBpgRnfmzvf1PDyemXuY+V4N/nA559xzxBiDUkqp4OewugCllFK+oYGulFI2oYGulFI2oYGulFI2oYGulFI24bLqg5OSkkx2drZVH6+UUkFp27ZttcaY5MGOWRbo2dnZ5OfnW/XxSikVlETk6OmOaZeLUkrZhAa6UkrZhAa6UkrZhAa6UkrZhAa6UkrZhAa6UkrZhAa6UkrZhGXz0NWZHW/q4K29VdS1dJ56bXxKDMunpBAbGWZhZUqpQKWBHkDcbsNrO8p5eVspnxw5gTEg0nfs5LL14S4HF01M5saFmVw6OcW6YpVSAUcDPUBUN3Xw7Zd38dGhWsYlj+Abl03kqtljGZ8cA/SF/fZj9WzYU8WGPZXc9vutXDkrlR9cNZ3k2AiLq1dKBQIZasciEXkS+AJQbYyZMcjxG4HvAAI0A18zxuwa6oNzc3ON3vrf5+19VXzn1d20d/fy/S9M44YFmcjJS/NBdPW4eezDw/z63SKiwp08sHo6qy9IG8aKlVJWEZFtxpjcwY55Myj6FLDiDMePABcbY2YCPwQeO+sKQ9ir28q489ltpCdE88a/LOPGC7POGObQ1+1y7/KJvPnNZUxMieEbL+7k0b8dHqaKlVKBashAN8Z8CJw4w/GPjTH1nqd5QLqParO9dwqO8x+v7mbJhFG88rVFp7pXvDU+OYY/3LGQL8xK5Sdv7ucnbxaie8QqFbp83Yd+O/Dm6Q6KyJ3AnQCZmZk+/ujgsuVwHf/8h+3MSIvn0ZtyiXA5z+l9wl0OfnX9HOKjwnj0b8W0dPTwoy/OGPIqXyllPz4LdBG5lL5AX3q6NsaYx/B0yeTm5obspeSR2lbueCafrMRonrp1PjER5/c1OB3Cj744g5hIF4/+rZi0hCjuuWSCj6pVSgULnwS6iMwCHgeuMMbU+eI97arXbfj2SztxCDz9lQUkjAj3yfuKCPetmEJFQwc/f+sA45JiWDFjjE/eWykVHM77TlERyQT+BNxkjDl4/iXZ26MfHmb7sQZ++MUZjB0Z5dP3FhF+cfUsLsgYyb/+cSd7yxt9+v5KqcA2ZKCLyAvAFmCyiJSJyO0icreI3O1pcj8wCvitiOwUEZ2LeBoFFU388p2DrJw5hlWzx/rlMyLDnDx28zwSosO445l8Gtq6/PI5SqnAM+Q8dH8JtXnoXT1uVj24idqWLt7+14tI9FFXy+nsLmvgy7/9mBUzxvCbtXN0kFQpmzjfeejKB57NO8r+qmZ+/KUZfg9zgFnpI/nm5RP56+5K/rKzwu+fp5Syngb6MGjq6ObB9w6xdEISn58+fAOVd188nnlZCXz/L3spb2gfts9VSllDA30YPPLBYerburnviinD+rkup4NfXnsBbs/MGrc7ZGeKKhUSNND9rKqxgyc3H2H1BWOZkRY/7J+fOSqa739hGnnFJ3hle9mwf75SavhooPvZL985SK/b8G+fn2xZDdfmZpCblcBPNhRS36qzXpSyKw10PyqqbuHlbaWsW5hFRmK0ZXU4HMIPvziDpo4efr7xgGV1KKX8SwPdj373YXHfyoiXWn8b/tTUOG5bnM2LW4+x41j90L9BKRV0NND9pKa5k9d2lLNmbjqjYgJjA4pvfm4SKbERfO/Pe+nVAVKlbEcD3U+e3VJCt9vN7UtzrC7llJgIF9+7chr7Kpp4VQdIlbIdDXQ/aO/q5dm8o1w2ZTTjznKNc3/7wqxULsgYyf+8fZD2rl6ry1FK+ZAGuh+8ur2M+rZu7lgWOFfnJ4kI/7lyKlVNfdMplVL2oYHuY2634clNR5iVHs+CnESryxnUgpxELp86moc/OExdS6fV5SilfEQD3cfeP1BNcW0rX102LqAXxLrvism0dfXwm/eKrC5FKeUjGug+9uLWUpJiIrgiwDeXmJASy3XzM3ku7yjH6tqsLkcp5QMa6D5U3dzBe/urWTM3jTBn4P9ov3n5RBwO4cH3D1ldilLKBwI/dYLIn3eU0+s2XJObYXUpXhkdF8kNCzJ5dXs5R+tarS5HKXWeNNB9xBjDS/llzMtKYEJKYE1VPJN7LhmPyyE8qH3pSgU9DXQf2VHaQFF1C9fmpltdyllJiYvkxguz+NOOckpq9SpdqWCmge4jL+eXEhXm5MpZ/tkr1J/uvngcLofojBelgpwGug+0dfXw+q5KrpyVSkyEy+pyzlpKXCTrFmbx2o4yvUpXKohpoPvAW3uraOns4br5wTEYOpi7Lh6Hy+ng0Q8PW12KUuocaaD7wOu7KkhPiCI3K8HqUs5ZSmwk18xL59Vt5Rxv6rC6HKXUOdBAP0+Nbd18dKiWK2emBvSdod6466Lx9LjdPLFJ13hRKhgNGegi8qSIVIvI3tMcFxH5tYgUichuEZnr+zID18aCKnrchitnpVpdynnLHBXNVbPH8nzeURrbuq0uRyl1lry5Qn8KWHGG41cAEz2/7gQePv+ygscbuyvJSIxipgUbQPvD3RePp7Wrl2e2lFhdilLqLA0Z6MaYD4ETZ2iyGnjG9MkDRopI8F+ueqG+tYvNRbVcOXNs0He3nDQ1NY7lU1L4/cclul66UkHGF33oaUBpv+dlntf+gYjcKSL5IpJfU1Pjg4+21tsnu1tm2uvvr69dMp4TrV28vK106MZKqYAxrIOixpjHjDG5xpjc5OTk4fxov/jr7koyE6OZkRZndSk+lZuVwAUZI3li0xHde1SpIOKLQC8H+k/ATve8ZmsnWrv4+HAdV84K/tktA4kIdywbx9G6Nv638LjV5SilvOSLQF8P3OyZ7bIQaDTGVPrgfQPaxn1V9Nqwu+Wkf5o+mvSEKB7/qNjqUpRSXvJm2uILwBZgsoiUicjtInK3iNztabIBKAaKgN8B9/it2gDyTsFx0hOimD7WXt0tJ7mcDr6yJIetJfXsLG2wuhyllBeGXHjEGLN2iOMG+GefVRQE2rt62VxUy9oFmbbrbunv2vkZ/PJ/D/L4R8U8eENI3V6gVFDSO0XPweaiWjp73Fw2NcXqUvwqJsLFDQsyeXNvFWX1uk2dUoFOA/0cvLv/ODERLi7MGWV1KX5365JsAJ7+uMTSOpRSQ9NAP0tut+HdwmoumpREuMv+P77U+CiumDGGP24tpa2rx+pylFJnYP9E8rF9FU1UN3eyfMpoq0sZNrctyaapo4dXt9t+NqpSQU0D/Sz9b+FxRODSycF/Y5S35mYmMCs9nqc2H6FvDFwpFYg00M/Se/urmZuZwKiYCKtLGTYiwq2Lszlc08pHh2qtLkcpdRoa6GfheFMHe8obWT7F3rNbBnPlrFSSYiJ4SgdHlQpYGuhn4b391QBcPjV0+s9PinA5ufHCTN7bX80R3XdUqYCkgX4W3t9fTdrIKCaNjrG6FEvcuDCTMKfoFEalApQGupd6et1sOVzHRZOSbH136JmkxEaycmYqr24ro7VTpzAqFWg00L20q6yB5s4elk0Mndktg7l5UTbNnT28tkOnMCoVaDTQvfTRoVpEYPF4+98deiZzM0cyIy2OZ7aU6BRGpQKMBrqXPjpUy6z0kYyMDre6FEuJCDcvyubg8Rbyis+0M6FSarhpoHuhqaObnaUNLJuQZHUpAWHV7LGMjA7TjaSVCjAa6F7IO1xHr9uwbKIGOkBkmJPrcjN4u+A4lY3tVpejlPLQQPfCR4dqiQ53MiczwepSAsa6hVm4jeH5vGNWl6KU8tBA98KmoloWjRsVEqsreisjMZrlk1N4cWspXT1uq8tRSqGBPqTSE20cqW1lqXa3/IN1i7KobenkrX1VVpeilEIDfUibivoWowr1+eeDuXhiMpmJ0Ty35ajVpSil0EAf0qZDtaTGRzI+eYTVpQQch0NYtzCTT0tOsL+qyepylAp5GuhnYIwhr7iOReNHhezt/kO5Zl4G4S4Hz+XpVbpSVtNAP4Oi6hbqWrtYOC607w49k4QR4Vw1ayyvbS+nuaPb6nKUCmka6GeQV1wHwCIN9DO6aVEWrV29ur6LUhbzKtBFZIWIHBCRIhG5b5DjmSLyvojsEJHdIrLS96UOv7ziE4yNjyQ9IcrqUgLa7PR4ZqbF81zeUV3fRSkLDRnoIuIEHgKuAKYBa0Vk2oBm3wNeMsbMAa4HfuvrQofbyf7zheO0/3woIsJNC7M4eLyFrSX1VpejVMjy5gp9AVBkjCk2xnQBLwKrB7QxQJzncTxQ4bsSraH952fnqtljiYt08awOjiplGW8CPQ0o7fe8zPNafz8A1olIGbAB+PpgbyQid4pIvojk19TUnEO5wyfvSN9Kghro3okKd3L1vAze2ltJTXOn1eUoFZJ8NSi6FnjKGJMOrASeFZF/eG9jzGPGmFxjTG5ycmDfqJNXXMfY+EgyErX/3Fs3Lsyku9fwUn7p0I2VUj7nTaCXAxn9nqd7XuvvduAlAGPMFiASCNp75Y0xfKL952dtfHIMSyaM4g+fHKPXrYOjSg03bwJ9KzBRRHJEJJy+Qc/1A9ocAy4DEJGp9AV6YPepnMHhmhZqW7T//FysuzCL8oZ23t9fbXUpSoWcIQPdGNMD3AtsBArpm82yT0QeEJFVnmbfBu4QkV3AC8CtJojnr20p1v7zc3X5tNGMjovQwVGlLODyppExZgN9g539X7u/3+MCYIlvS7OO9p+fuzCng+vnZ/Lr9w5xtK6VrFG6Bo5Sw0XvFB3AGMPWIyeYn5Oo/efnaO2CTBwi/OET3fxCqeGkgT5AWX071c2d5GYnWl1K0BoTH8nnp43mpfxSOrp7rS5HqZChgT7A1pK+/vPcLN1u7nysW5hFfVs3G/ZUWl2KUiFDA32A/KP1xEa6mDQ61upSgtri8aMYlzxCB0eVGkYa6APkl5xgbmYCTof2n58PEWHdhVnsONbA3vJGq8tRKiRooPfT2NbNweMtzM/W7hZfWDMvncgw3fxCqeGigd7PtmN9/efzsnRA1Bfio8JYPTuNv+ysoLFdN79Qyt800PvJL6nH5RAuyBhpdSm2cdOiLNq7e3llW5nVpShlexro/eSX1DM9LZ6ocKfVpdjGjLR45mSO5Lm8o7h1fRel/EoD3aOzp5ddZQ06XdEPbl6UxZHaVjYfrrW6FKVsTQPdY295E509bh0Q9YOVM1NJHBHOM1t0cFQpf9JA99h2VAdE/SXC5eS6+Rm8W3ic8oZ2q8tRyrY00D22ltSTPSqa5NgIq0uxpRsvzATgD5/oVbpS/qKBTt+CXDuO1TNX+8/9Jj0hmuVTRvPip6V09uj6Lkr5gwY6fQty1bZ0MTdTA92fblmcRV1rl67vopSfaKAD24/VAzAnc6S1hdjckvFJjEsewVMfa7eLUv6ggQ7sONZAVJiTybogl185HMIti7LZVdrAztIGq8tRynY00IEdpQ3MSo/H5dQfh799eW4aI8KdPLOlxOpSlLKdkE+wju5eCioamaP958MiNjKMNfPS+euuSmpbOq0uRylbCflA31fRRHev0f7zYXTzomy6et38cWup1aUoZSshH+g7Tg6I6oJcw2ZCSgxLJyTxXN5RunvdVpejlG1ooJc2kDYyipS4SKtLCSm3LM6msrGDjfuqrC5FKdsI+UDfeaxBu1sssHxKCpmJ0Ty1ucTqUpSyDa8CXURWiMgBESkSkftO0+ZaESkQkX0i8gfflukfx5s6KG9o1wFRCzgdwi2Ls8k/Ws+eMt2iTilfGDLQRcQJPARcAUwD1orItAFtJgL/B1hijJkOfNP3pfrejmMNgN5QZJVrctMZEe7k95uPWF2KUrbgzRX6AqDIGFNsjOkCXgRWD2hzB/CQMaYewBhT7dsy/WNnaQPhTgfTx8ZZXUpIiosM45rcDF7fXUF1c4fV5SgV9LwJ9DSg//yyMs9r/U0CJonIZhHJE5EVvirQn3Ycq2fa2DgiXLpDkVVuWZxNj9vwfN4xq0tRKuj5alDUBUwELgHWAr8TkZEDG4nInSKSLyL5NTU1Pvroc9PT62Z3WaPuH2qxnKQRXDo5hec/OaqrMCp1nrwJ9HIgo9/zdM9r/ZUB640x3caYI8BB+gL+M4wxjxljco0xucnJyedas08U1bTQ3t3L7Ix4S+tQ8JUlOdS2dLF+Z4XVpSgV1LwJ9K3ARBHJEZFw4Hpg/YA2f6bv6hwRSaKvC6bYd2X63u7SvpkVs9JHWluIYsmEUUwZE8sTm45gjG4krdS5GjLQjTE9wL3ARqAQeMkYs09EHhCRVZ5mG4E6ESkA3gf+3RhT56+ifWFXWQOxkS5yRo2wupSQJyLcvjSH/VXNbCrSjaSVOlde9aEbYzYYYyYZY8YbY/6f57X7jTHrPY+NMeZbxphpxpiZxpgX/Vm0L+wua2RWejwOh1hdigJWXTCW5NgIHv9IpzAqda5C8k7Rju5eCiubtLslgES4nNyyKIu/Hazh4PFmq8tRKiiFZKAXVjbR4zbMTtcB0UByw4VZRIY5eHKTXqUrdS5CMtB3e241n61TFgNK4ohw1sxN5087yqlp1rXSlTpbIRnou8oaSI6NYIyusBhwbl+aQ3evm6c/LrG6FKWCTmgGemkDs9PjEdEB0UAzLjmGf5o2hme2lNDa2WN1OUoFlZAL9OaOboprW3VANIDddfE4mjp6eFF3NFLqrIRcoO8pb8QY7T8PZHMyE1iQk8iTm47ojkZKnYWQC/STA6Kz0nSGSyC7++JxlDe088buSqtLUSpohFyg7yptIDMxmoQR4VaXos7gkkkpTBodwyN/O6zLASjlpZAL9JN3iKrA5nAId140nv1VzXxwwNqVOZUKFiEV6HUtnZQ3tGugB4nVF4wlbWQUD75fpFfpSnkhpAJ9T3lf//nMtJHWFqK8EuZ0cNfF49h2tJ5PjpywuhylAl5IBfpeT6BPT9Mt54LFtbkZJMVE8ND7RVaXolTAC6lA31PeSE7SCOIiw6wuRXkpMszJV5fl8NGhWnaVNlhdjlIBLaQCfW95EzN0umLQWbcwi/ioMB7Uq3SlzihkAv1EaxflDe3M1O6WoBMT4eLWxdm8U3Cc/VVNVpejVMAKmUA/OSCqV+jB6bYl2cREuPj1u4esLkWpgBUygb5XAz2ojYwO57Yl2WzYU0VhpV6lKzWYkAn0PWWNZI+K1gHRIHb70hxi9SpdqdMKnUAvb9Sr8yB38ir9zb1VFFToVbpSA4VEoNefGhDVQA92ty8dp1fpSp1GSAT63+8Q1UAPdvHRYdy2NIe39lWxr6LR6nKUCighFejTNdBt4falOcRFuvivjQesLkWpgBISgb63vJGsUdHER+mAqB3ER4Vx9yXjef9ADVtLdI0XpU7yKtBFZIWIHBCRIhG57wzt1oiIEZFc35V4/nRA1H5uW5xDcmwEv3jrgK7EqJTHkIEuIk7gIeAKYBqwVkSmDdIuFvgG8ImvizwfDW1dlNW3M2OsBrqdRIU7+ZflE/i05AQfHNT10pUC767QFwBFxphiY0wX8CKwepB2PwR+BnT4sL7zts8zvU0HRO3nuvmZZCRG8Yu3DuB261W6Ut4EehrQf/v1Ms9rp4jIXCDDGPPGmd5IRO4UkXwRya+pGZ6rqlNL5o7VNVzsJtzl4Fufm0RBZROv766wuhylLHfeg6Ii4gD+B/j2UG2NMY8ZY3KNMbnJycnn+9Fe2VfRxNj4SN1D1KZWzU5jWmocP3/rAB3dvVaXo5SlvAn0ciCj3/N0z2snxQIzgA9EpARYCKwPlIHRfRWNOl3RxpwO4btXTqW8oZ2nPy6xuhylLOVNoG8FJopIjoiEA9cD608eNMY0GmOSjDHZxphsIA9YZYzJ90vFZ6G1s4fi2lbtbrG5JROSuHRyMg++X8SJ1i6ry1HKMkMGujGmB7gX2AgUAi8ZY/aJyAMissrfBZ6P/VVNGIPOcAkB/7lyKq2dPbokgAppLm8aGWM2ABsGvHb/adpecv5l+cbe8r4ZLrqHqP1NHB3L9QsyeS7vKDcvymJccozVJSk17Gx9p+i+ikZGjQhnTFyk1aWoYfCvl08iMszJj94otLoUpSxh60DfW97EtLFxiIjVpahhkBwbwTcum8h7+6t5t/C41eUoNexsG+hdPW4OVTczXfvPQ8oti7MZnzyCB/5aoNMYVcixbaAfPN5Md69hhvafh5Rwl4MfrJrO0bo2nth0xOpylBpWtg30k2tl6xV66Fk2MZkV08fw4HtFVDS0W12OUsPGxoHeREyEi6zEaKtLURb43hem4jaGB14vsLoUpYaNbQN9b3kj01LjcDh0QDQUpSdE843LJ/LWvireKdABUhUabBnovW5DYWWzzj8PcXcsG8eUMbHc/5e9tHT2WF2OUn5ny0A/UttKe3cv01I10ENZmNPBj788k6qmDv77bd2uTtmfLQO9oNJzh6gOiIa8uZkJrLswi6c/LmFXaYPV5SjlV7YM9H0VjYQ7HUxI0du/Ffz7iskkx0bwH6/sprNH56Yr+7JloBdUNDFxdAzhLluenjpLcZFh/OTLMzlwvFkX71K2ZrvEM8ZQUNGk/efqM5ZPGc0189J5+IPD2vWibMt2gV7T3EldaxfTdA10NcD3r5rG6LhIvv3yLl0WQNmS7QL95KbQOiCqBoqLDONna2ZRVN2is16ULdku0E/OcJmSGmtxJSoQXTQpmZsWZvG7j47w4cHh2ahcqeFiv0CvaCIzMZq4yDCrS1EB6rtXTmXS6Bi+9dIuals6rS5HKZ+xXaDvq2jUAVF1RpFhTn6zdi7NHd3828u7cLuN1SUp5RO2CvSWzh5K6tp0U2g1pMljYvnelVP54EANT27WZXaVPdgq0Pd7+s91hovyxrqFWXx+2mh++uZ+8ktOWF2OUufNVoFeoIGuzoKI8ItrZpOeEMU9z2+nurnD6pKUOi+2CvR95U0kRIfpptDKa/FRYTxy0zyaO3q49/kddPe6rS5JqXNmq0AvqGxi+th43RRanZUpY+L46ZqZfFpygh9vKLS6HKXOmVeBLiIrROSAiBSJyH2DHP+WiBSIyG4ReVdEsnxf6pl197o5cLxZu1vUOVl9QRq3Lcnm95tLeOHTY1aXo9Q5GTLQRcQJPARcAUwD1orItAHNdgC5xphZwCvAz31d6FCKa1rp6nEzVW8oUufouyuncsnkZL7/571sLqq1uhylzpo3V+gLgCJjTLExpgt4EVjdv4Ex5n1jTJvnaR6Q7tsyh1ZQqZtCq/Pjcjr4zdo5jE+O4e7ntlFU3WJ1SUqdFW8CPQ0o7fe8zPPa6dwOvDnYARG5U0TyRSS/psa3t10XVDQR7nIwLmmET99XhZbYyDCeuDWXCJeD2576lOomnfmigodPB0VFZB2QC/xisOPGmMeMMbnGmNzk5GRffjSFlc1MHh2Ly2mrcV5lgfSEaJ64ZT51LV3c/OSnNLZ3W12SUl7xJv3KgYx+z9M9r32GiFwOfBdYZYwZ1gUyjDEUVDZp/7nymdkZI3nsplwO17Tw1ae30t6ly+2qwOdNoG8FJopIjoiEA9cD6/s3EJE5wKP0hXm178s8s+rmTk60dukaLsqnlk5M4lfXzyH/aD33PL9Nt69TAW/IQDfG9AD3AhuBQuAlY8w+EXlARFZ5mv0CiAFeFpGdIrL+NG/nFwWeNdCnaqArH1s5M5Uff2km7x+o4Z7ntmuoq4Dm8qaRMWYDsGHAa/f3e3y5j+s6Kydv+Z+qc9CVH6xdkInbGL772l7ufnYbD6+bR2SY0+qylPoHthhBLKhsIj0hStdAV35z44VZp67U73p2m/apq4Bki0AvrNRNoZX/3XBhJj9bM5MPD9Vw4+N5NLR1WV2SUp8R9IHe1tXDkdpW7T9Xw+K6+Zn89oa57C1v4upHtlDR0G51SUqdEvSBfqCqGWN0yVw1fK6YmcrTX1nA8cYO1jz8MXvLG60uSSnABoF+ag10vUJXw2jR+FH88a5FAFzzyBY27Km0uCKlbBDohZVNxEa4SE+IsroUFWKmjY3jL/cuYWpqLPc8v51fvnNQ9ydVlgr6QC+oaGJqapyuga4skRIbyQt3LmTN3HR+9e4hbn1qK7Utw3qjtFKnBHWgu92G/VW6BrqyVoTLyX9dM4sffXEGecV1rPzVR+QV11ldlgpBQR3ox0600dbVq2u4KMuJCOsWZvHaPYsZEeHiht/l8bO39uudpWpYBXWgn7pDVAdEVYCYPjae17++lKvnpfPwB4e56jeb2F3WYHVZKkQEdaAXVjbhEJg0Wq/QVeCIiXDx86tn8/vb5tPY3s2XfvsxP/prAS2dPVaXpmwu6AN9fHKMrquhAtKlk1N4+5sXc828dB7fdITL/vsDXt9VgTE6E0b5R5AHerN2t6iAFh8dxk/XzOJP9ywmKSaCr7+wg2sf3cK2o/VWl6ZsKGgDvbGtm/KGdg10FRTmZiaw/t6l/OiLMzhS28aahz/mrmfzOVDVbHVpyka8Wj43EP19QFT7z1VwcDr6ZsJ8aU4aT246wqMfFrNx34f80/TR3HvpRGam6wbn6vwEbaAX6i3/KkiNiHDx9csmsm5hFr//uISnNh9h477jLB4/ituW5LB8SgpOh94op85eUAf6qBHhJMdGWF2KUuckYUQ43/rcJO5YlsNzecd4ZksJdzyTT0ZiFDcsyGLN3DRS4iKtLlMFkeAN9Cq95V/ZQ2xkGF+7ZDx3LMvh7YLjPLW5hJ+9tZ//evsAF09K5stz01g+JYXo8KD931UNk6D8E9LT6+bg8RZuWZRldSlK+YzL6WDlzFRWzkyluKaFl7eV8aftZby3v5qoMCfLp6ZwxYwxXDQpWXfnUoMKykAvrm2lq8eta7go2xqXHMN3Vkzh3z4/mU+PnOCNPRW8uaeKN3ZX4nII87MTuXRKMovHJzEtNQ6H9rkrgjTQC/WWfxUinA5h0fhRLBo/ih9cNZ2dpQ28u7+adwuP8+MN+wEYGR3GhTmJ5GYlMjcrgRlpcUS49Ga7UBSUgV5Q2US408H45BirS1Fq2LicDnKzE8nNTuQ7K6ZQ1djBluJaNhfV8emRE2zcdxyAMKcweUwsM9PimT42niljYpk4Opb4KO2msbvgDPSKJiakxBDmDNr7opQ6b2PiI/nSnHS+NCcdgOrmDrYfrWfHsQb2VjTyxu5KXvi09FT71PhIxiWPICdpBDlJMWQkRJGRGE1GYjQxEUEZBWoAr75FEVkB/ApwAo8bY3464HgE8AwwD6gDrjPGlPi21L8rrGzm4knJ/np7pYJSSmwkK2aksmJGKgDGGMrq2zlU3cyBqhYOHW+muLaV9TsraOr47EJhcZEuUuOjSB0ZSUpsBMmxESTHRJAYE8GoEeEken7FR4Xp2kkBbMhAFxEn8BDwOaAM2Coi640xBf2a3Q7UG2MmiMj1wM+A6/xRcE1zJ7UtnXqHqFJDEJFTV+DLp4w+9boxhvq2bkpPtFFa30bpiXaqGtupaOygsrGdwsomalu66D3NdnoRLgdxUWHERrqIjQwjNsLFiAgnI8JdREc4iQ53ERXmJDrcSYTLQWSYk8iwvscRYQ7CnU7CnEK4y0GY8+QvIczpwOkQXE7B5fA8dghOh+CQk/9FpyqfgTdX6AuAImNMMYCIvAisBvoH+mrgB57HrwAPiogYPywrp3eIKnV+ROTUFffsjJGDtnG7DSfaujjR2kVdS99/G9q7aGzvpqGtm6b2bpo7emjq6Kals4ea5k5aOnto6+qhvbuXjm633+p3CDhEcHgC3iF9gS+nHvedo3jOVQTPYxD6ngOnjnPyWL+/Jwa2O/V6v0af+WvlNH/HnO6vnrULMvnqsnFen7O3vAn0NKC03/My4MLTtTHG9IhIIzAKqO3fSETuBO4EyMzMPKeCo8KdXD41RWe4KOVHDoeQFBNBUkwEjB66/UBut6Gjpy/YO7p76ejupavXTWe3m84eN929brp63HT1uunpNfS4+573ug3dbkNPb99jtzH0uA1ut6HXDb3GYIyh123oNQYMuE3fMYPBGOh1m1OP+/6R0fe4r3nfNabnt/79tX6Xnn2vm1OPT70+oM3fXx/8uvVMV7NJMf65w31YR0KMMY8BjwHk5uae09X7/OxE5mcn+rQupZRvORxCdLiL6HCrKwkt3kwTKQcy+j1P97w2aBsRcQHx9A2OKqWUGibeBPpWYKKI5IhIOHA9sH5Am/XALZ7HVwPv+aP/XCml1OkN2eXi6RO/F9hI37TFJ40x+0TkASDfGLMeeAJ4VkSKgBP0hb5SSqlh5FUfujFmA7BhwGv393vcAVzj29KUUkqdDb3VUimlbEIDXSmlbEIDXSmlbEIDXSmlbEKsml0oIjXA0XP87UkMuAs1BOg5hwY959BwPuecZYwZdHVCywL9fIhIvjEm1+o6hpOec2jQcw4N/jpn7XJRSimb0EBXSimbCNZAf8zqAiyg5xwa9JxDg1/OOSj70JVSSv2jYL1CV0opNYAGulJK2URAB7qIrBCRAyJSJCL3DXI8QkT+6Dn+iYhkW1CmT3lxzreKSI2I7PT8+qoVdfqKiDwpItUisvc0x0VEfu35eewWkbnDXaOveXHOl4hIY7/v+P7B2gUTEckQkfdFpEBE9onINwZpY6vv2stz9u13bTxbOgXaL/qW6j0MjAPCgV3AtAFt7gEe8Ty+Hvij1XUPwznfCjxoda0+POeLgLnA3tMcXwm8Sd/2jAuBT6yueRjO+RLgr1bX6eNzTgXmeh7HAgcH+bNtq+/ay3P26XcdyFfopzanNsZ0ASc3p+5vNfC05/ErwGUS3FuCe3POtmKM+ZC+NfRPZzXwjOmTB4wUkdThqc4/vDhn2zHGVBpjtnseNwOF9O1F3J+tvmsvz9mnAjnQB9uceuAP4zObUwMnN6cOVt6cM8Aazz9JXxGRjEGO24m3PxO7WSQiu0TkTRGZbnUxvuTpGp0DfDLgkG2/6zOcM/jwuw7kQFeDex3INsbMAt7h7/9CUfaxnb71OmYDvwH+bG05viMiMcCrwDeNMU1W1zMchjhnn37XgRzoobg59ZDnbIypM8Z0ep4+Dswbptqs4s2fA1sxxjQZY1o8jzcAYSKSZHFZ501EwugLtueNMX8apIntvuuhztnX33UgB3oobk495DkP6FNcRV+/nJ2tB272zIBYCDQaYyqtLsqfRGTMybEgEVlA3/+nwXyhgud8ngAKjTH/c5pmtvquvTlnX3/XXu0pagUTgptTe3nO/yIiq4Ae+s75VssK9gEReYG+kf4kESkD/i8QBmCMeYS+vWxXAkVAG3CbNZX6jhfnfDXwNRHpAdqB64P8QgVgCXATsEdEdnpe+08gE2z7XXtzzj79rvXWf6WUsolA7nJRSil1FjTQlVLKJjTQlVLKJjTQlVLKJjTQlVLKJjTQlVLKJjTQlVLKJv4/PwwDougOVOcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "vs = np.linspace(0,2.5,100)\n",
    "prob = prob_dist(vs)\n",
    "plt.plot(vs,prob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "id": "CcZlX78AYoh5"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxRklEQVR4nO3dd3xUVfrH8c+TRq8hSCBUQZokASIgoKACgrogiiCrlFVA/NnYtaCuEEa3iN1FEVBZRKUJgqygIK6AiJQAIVJEEEESWkB6TTm/P+4kG2JIBjKTO3Pneb9e85py79z53gw8OTnn3nPFGINSSqnAF2J3AKWUUt6hBV0ppRxCC7pSSjmEFnSllHIILehKKeUQWtCVUsohtKArRxCXLBWXDCnG+58Vl7znzUxKlbQwuwOowCYu6Qi8BDQHsoCtwAiTaNbaGiwPcckYoKFJNPfme90AjUyi2WESzT883NZS4COTaLT4K7+jLXR12cQlFYHPgXFAVaAW4ALO2ZnLTuISrzaSvL095Wz6j0UVx1UAJtFMdz8/AyzOWSguGQwMBTYAA4B9wEMm0XztXl4JeA24BcgG/g0kmkST5V5+H/AkUANYAwwziWa3e1lXrF8k0cCHgBRnR/K24sUlpYH3gB5AKLAduA14FLgOaCcueQOYYhLNw+6W/sPACKz/U/XFJUOBkVi/6FYAw02i2ev+rG7u7DWAj7H+uvnQJJr38vzM1gADgXfEJf8G3gXiAAMscv8cj7q3twt42/0zvhKYATwLTAE6AquBu0yiOVKcn5Hyf9pCV8XxE5AlLvlAXNJDXFKlgHXaAj8D1YBE4FNxSVX3silAJtAQaAl0A4YAiEt6YRWlO4Ao4FtguntZNeBT4Dn3dn8GOnhxvwYBlYDaQCQwHDhjEs1f3TkeNommvEk0D+d5z+3ufW0mLrkR+CfQF+sXzm6sIpuTfTbwjHvb24D2+T6/LbATuAL4O9Yvq38CNYGm7lxj8r3nTqAr1i/ZPwBfYP38orD+nz96mT8LFUC0ha4um0k0x9196COxWpA1xCULgaEm0Rxwr3YQeMMkGgPMFJc8DtwqLlmM1TKvbBLNGeCUuOR1YBgwEauI/tMkmq0A4pJ/AM+KS+oCnYDNJtHMdi97A3i8iLh9xSW3ebhrGVjFtqFJNCnAOg/e80+TaH5z57kHmGwSzXr382eAI+KSesD17uyfupf9C3gi37b2mkQzzv04E9jhvgGki0tew/rlmNe4nJ+5uORb4KBJNBvcz+cCN3m05yqgaUFXxeIuuIMBxCVNgI+AN4D+7lXS3MU8x26slmZdIBzYJ67c3pIQYI/7cV3gTXHJq3neK1j99DXzrIdJNEZcsofCzbrIoGhBPsRqBc8Ql1R279NfTaLJKGT7eT+/JrA+T76T4pLDhWRPLWRbiEuuAN7E6u6pgPVzyt99ciDP4zMFPC9fSHblEFrQldeYRPOjuGQK8ECel2uJSyRPUa8DzMcqWueAaibRZBawuT3A302i+Tj/AnFJI6yCm/Nc8j4vLnfhdgEud6t6IVbXyPtYfdgFvi3P471Yv5By8pXDavGnYY0jxOTLHsOF8n/GP9yvtTCJ5jdxye3AW5e0UyooaEFXl83dIr8VmGkSTaq4pDZWy3xVntWqA4+KS8Zj9TM3BRaaRHPY3e3yqrhkFHASqA/EmESzDJgAvCAuSTaJZrN7ALWbSTSfAAuAt8Qld2D9cngIa4DRW/t1A3AI2AIcx+qCyXYvPgA0KGIT04Hp4pJpWIdx/gNYbRLNLnHJSXf227GOEBruQfYKwDHgmLikFtZAsVK/o4OiqjhOYA3grRaXnMIq5Ju4sD97NdAIq0D+HehjEs1h97KBQARW4TyCNVgYDWASzVxgLFa3x3H3dnu4lx0C7gJeBA67t/+dF/erhjvLcayCvAyrGwasro8+4pIj7v7v3zGJZgkwCpiD1SK/Erg7X/aX3NmbAUkUfqinC2iFVdQXYA0IK/U7ohe4UL7iPgRviEk0He3O4q/EJSFAKnCPSTTf2J1HBTbtclGqhIlLbsb6y+UMVveJcGE3lVKXRbtclCp512IdO38I65jx292HbipVLNrlopRSDqEtdKWUcgjb+tCrVatm6tWrZ9fHK6VUQFq3bt0hY0xUQctsK+j16tUjKSnJro9XSqmAJCK7L7ZMu1yUUsohtKArpZRDaEFXSimH0BOLlHKAjIwMUlNTOXv2rN1RlJeULl2amJgYwsPDPX6PFnSlHCA1NZUKFSpQr149RIp18SblB4wxHD58mNTUVOrXr+/x+7TLRSkHOHv2LJGRkVrMHUJEiIyMvOS/uLSgK+UQWsyd5XK+Ty3oSinlEFrQlVLFdvToUcaPH5/7fOnSpdx2m6eXcL00EyZMYOrUqQBMmTKFvXv35i4bMmQIW7ZsueRt+jJvSdKC7mDRMVmIcMEtOibL7ljKgfIXdF/JzMxk+PDhDBw4EPh9QX/vvfdo1qyZz3N4IjOzoCsr+pYWdAfbnxZK3ZELLrjtTwu1O5ZyoKeffpqff/6Z+Ph4nnzSukLeyZMn6dOnD02aNOGee+4hZ2bXdevW0alTJ1q3bs3NN9/Mvn37AEhOTqZdu3bExsbSu3dvjhyxroPduXNnRowYQUJCAm+++SZjxozhlVdeYfbs2SQlJXHPPfcQHx/PmTNn6Ny5c+6UIl9++SWtWrUiLi6Om266CYA1a9Zw7bXX0rJlS9q3b8+2bdsK3a+srCyeeOIJrr76amJjYxk3bhxgTV1y6NAhAJKSkujcuTMAY8aMYcCAAXTo0IEBAwbQrl07Nm/enLu9nHynTp3ivvvuo02bNrRs2ZLPPvvMG1+DHrbo14yB7Oz/PQ/VYqw8MGIEJCd7d5vx8fDGGxdd/OKLL7Jp0yaS3Z+7dOlSNmzYwObNm6lZsyYdOnTgu+++o23btjzyyCN89tlnREVFMXPmTP76178yefJkBg4cyLhx4+jUqROjR4/G5XLxhvszz58/n1uox4wZA0CfPn146623eOWVV0hISLggT3p6OkOHDmX58uXUr1+f3377DYAmTZrw7bffEhYWxpIlS3j22WeZM2fORfdr0qRJ7Nq1i+TkZMLCwnK3U5gtW7awYsUKypQpw+uvv86sWbNwuVzs27ePffv2kZCQwLPPPsuNN97I5MmTOXr0KG3atKFLly6UK1euyO0XRgu6v8nOhpUr4eOP4ZNP4PDh3EU/0YhZ9OXr6n34Zn+c1YeilJ9q06YNMTExAMTHx7Nr1y4qV67Mpk2b6Nq1K2C1gKOjozl27BhHjx6lU6dOAAwaNIi77rord1v9+vW7pM9etWoV119/fe4x3FWrVgXg2LFjDBo0iO3btyMiZGRkFLqdJUuWMHz4cMLCwi7YTmF69uxJmTJlAOjbty/dunXD5XIxa9Ys+vTpA8DixYuZP38+r7zyCmAddvrrr7/StGnTS9rP/LSg+5M1a2DQIPjxRyhTBnr1YtSMZlTuuI1Qk03r1K088+s/ee7g36HTdfDuu9C4sd2plb8ppCVdkkqVKpX7ODQ0lMzMTIwxNG/enO+///6CdY8dO1botorbcs0xatQobrjhBubOncuuXbtyu0ouVVhYGNnuv57zHyueN2utWrWIjIwkJSWFmTNnMmHCBMA6cWjOnDk09vL/X+1D9wfZ2fDSS9ChA5w+DVOnwsGDMH06f2MU4zr0542O9zDg7r/R5qGpPMqbsGkTxMVZ77Nh8EWpvCpUqMCJEyeKXK9x48akp6fnFvSMjAw2b95MpUqVqFKlCt9++y0AH374YW5r/XI+t127dixfvpxffvkFILer5NixY9SqVQuwBlSL0rVrVyZOnJg7wJmznXr16rFu3TqAQrtswPrr4qWXXuLYsWPExsYCcPPNNzNu3LjccYUNGzYUmcUTRRZ0EZksIgdFZNNFlt8jIiki8oOIrBSROK8kCwLRMVmUk1N8GXoLjBzJ52E9rb7PAQOgfPkC33O4XGXG8Shs2QK3WO/jxhuhiBaOUr4UGRlJhw4duPrqq3MHRQsSERHB7NmzGTlyJHFxccTHx7Ny5UoAPvjgA5588kliY2NJTk5m9OjRRX7u4MGDGT58eO6gaI6oqCgmTZrEHXfcQVxcXG6XzVNPPcUzzzxDy5YtPToKZciQIdSpU4fY2Fji4uKYNm0aAImJiTz22GMkJCQQWsTYVp8+fZgxYwZ9+/bNfW3UqFFkZGQQGxtL8+bNGTVqVJFZPFHkNUVF5HrgJDDVGHN1AcvbA1uNMUdEpAcwxhjTtqgPTkhIMMF+gYtwyWBRg/Z0+mU9o7sO5x9L3oTs3/eC1R254ILnu8feijFYg6Yffwz33QctWsCiRVCtWu56IoW8VznK1q1bi93/qvxPQd+riKwzxiQUtH6RfejGmOUiUq+Q5SvzPF0FxHgWNcgZw7sM5cadSTxz88NMj+8Oi8MKLMAXJQL33gtVq8Kdd0KnTrBkCURH+zi8UsofebsP/X7gi4stFJFhIpIkIknp6ele/ugA8+yzDOYDXu/wR6uYF8ctt8DChbB7N9xwA7iP3y1Q6O9PNgorlaUnHynlAF47ykVEbsAq6B0vto4xZhIwCawuF299dsCZNw9efJEJPMCbHbxzunH0gOtpdGohS7Z14b9V+3MrC4AC+vayQgv8KyDva4X+VaD8ljFGJ+hykKK6wwvilRa6iMQC7wG9jDGHi1o/qB0+DMOHQ3w8jzDu8o4lL6CVvT8tlF9HnmD0zQ/QnUW8c81dRW9HOUbp0qU5fPjwZRUB5X9y5kMvXbr0Jb2v2C10EakDfAoMMMb8VNztOd4jj1hFfdEiMuM9vxLJBS7SygaYEd+dpum/MGztXFYwleVEFjexCgAxMTGkpqYS9F2ZDpJzxaJLUWRBF5HpQGegmoikAolAOIAxZgIwGogExrv/3Mu82Ahs0Js7F6ZPB5fLOobcR164cSiNDu1h4q8PcNvh1/k5srbPPkv5h/Dw8Eu6so1ypiK7XIwx/Y0x0caYcGNMjDHmfWPMBHcxxxgzxBhTxRgT775pMS/IkSPw4IPWnBjPPOPTj8oMDeOxPzzBacry6oLXCc3WQU6lgoGeKVpC3qrzD7IPHKRl8mQkItzn07Ckl6/Kg7xD/L6feGB14WeyKaWcQedyKQlpadx/8i3mNu/Mkdv2Uhdr/mZfH03yCX35T5O3GLFiGt9cmcDW6g18+nlKKXtpC70kvPACoWTxesd7SvyjR3V7kKNlKvDa568RlqVzvijlZFrQfW3HDnj/fSYxjNTKNUr844+Wqchfb36Ipum7GLT+8xL/fKVUydGC7mujR0NEBH/jOdsifNWwLd80aM2IFR9zBftty6GU8i0t6L70ww/WYYqPPcYBSr51nksE103DiMjKYCwj7cuhlPIpLei+9PrrULYsPPGE3UnYVbUW713Tm0FMpVXqVrvjKKV8QAu6rxw8aE1tO3iwNRuiH3jr2n7sIYbnl0wgRI9NV8pxtKD7yoQJcP48PPqo3UlynYkozZO8zNUHfub2LUvtjqOU8jIt6L5w7hyMH29Na+tn1/ycRV9SajTkL99+TERm4RfIVUoFFi3ovjBzJhw4ACNG2J3kdwwhvHT9IGKOH+SPyRedul4pFYC0oHubMdZV15s1gy5d7E5ToBX14vmubiwPfz+TcudO2x1HKeUlWtC9bcUK2LABHnvs8uY6LwkivHT9IKqdPsb9SZ/ZnUYp5SVa0L3t/fehQgXrWp9+bGPNxnxxVXuGrvmUSA7ZHUcp5QVa0L3p5EmYPRv69bOOP/dzr153L+XOn+XPvG53FKWUF2hB96Y5c+DUKevY8wCwo1odFjbuwCOMo+LZk3bHUUoVkxZ0b5oyBRo2hPbt7U7isbfa96MiJ7hP+9KVCnha0L1l1y5YuhQGDfLfwdAC/Fi9PvPoxZ+S5lNej3hRKqBpQfeWqVOt+wED7M1xGV5gFJXOnWKgTq+rVEDTgu4NxsAHH8CNN0LdunanuWTrac1/GyQwZO08yp4/Y3ccpdRl0oLuDStWwM6dVndLgBrX/m6qnjnO3RsX2x1FKXWZtKB7w/Tp1mGKd9xhd5LLtqFWE1bHNOf+tfMIQ+d4USoQaUEvrqws+PRTayKu8uXtTlMsE9veSa0T6fRllt1RlFKXQQt6ca1caU3Edeeddicptm+uTGB7ZG2e5GVrXEApFVCKLOgiMllEDorIpossFxH5l4jsEJEUEWnl/Zh+bPZsKFUKbr0196XomCxEuOAWCIyEMKlNb+LZCF99ZXccpdQl8qSFPgXoXsjyHkAj920Y8E7xYwWI7Gyru6V7d2v+Frf9aaHUHbngglug+KzZDewlGl5+2e4oSqlLVGRBN8YsB34rZJVewFRjWQVUFpFobwX0a2vWQGqqI7pbcpwPC+dNHoMlS6xZI5VSAcMbfei1gD15nqe6X/sdERkmIkkikpSenu6Fj7bZnDkQHg5/+IPdSbxqIg9AuXLWvO5KqYBRooOixphJxpgEY0xCVFRUSX609xlj9Z937QqVK9udxquOURn+9CeYMQP277c7jlLKQ94o6GlA7TzPY9yvOdv69db8LQ7qbrnAI49YF7meMMHuJEopD3mjoM8HBrqPdmkHHDPG7PPCdv3b3LkQGgq9etmdxDeuuso6cuedd6yLXiul/J4nhy1OB74HGotIqojcLyLDRWS4e5WFwE5gB/Au8H8+S+tPPv8cOnSAyEi7k/jOiBFw8KDV9aKU8nthRa1gjOlfxHIDPOS1RIEgLQ02boSxY+1O4ls33QTNm1uDowMHBs4B9UoFKT1T9HIsXGjd33KLvTl8JdR9YlSIMHTzY5CczB3Vv7E7lVKqCEW20FUBFiyAOnWgeXOiY7LYnxZqdyLvygrNPRlqeUYkR8ZXoP+hd4Ab7c2llCqUFvRLde6cddKNuwsi56zQvHaPvfUibw4858JLMSu2K/evmWt1NdUq8BQDpZQf0C6XS7V8uXUhaKd2txTgo5a3EEI2TJxodxSlVCG0oF+qBQugdGnr6kRBYk/lGizgVpg0yTo2XSnll7SgX6qFC+GGG6wLWgSRt3nImiZ4zhy7oyilLkIL+qXYvt26BVF3S47FdIOGDeGtt+yOopS6CC3ol2KBe/DzVucMenrKEAIPPWRd0ENnYVTKL2lBvxSLF0PjxlC/vt1J7DF4sNXVNH683UmUUgXQgu6p8+dh2TLo0sXuJPapXBn694dp0+DYMbvTKKXy0YLuqdWr4fTp4C7oAA8+aP0cpk61O4lSKh8t6J5asgRCQqBzZ7uT2CNnOoCE1qzhGrY8+g7RtTLtTqWUykMLuqeWLIGEBMddzMJjWf+7TuonPdrTjK1ctfc7u1MppfLQgu6J48etLpebbrI7iV/4T9PrOFaqHA8G0fXAlQoEWtA9sXw5ZGVp/7nb2fDSzG7RhTv41DrZSCnlF7Sge2LJEut0//bt7U7iNz6O70EEGfD++3ZHUUq5aUH3xNdfQ8eOVlFXAOyMjOG/dOaXv75LiGRbA6YC0TFZdkdTKmhpQS/K/v2waZN2txRgAg9Sn10MuOtvuQOmjpsbXqkAogW9KF9/bd3rgOjvzON20stW5p7kL+yOopRCC3rRvv4aqlSBli3tTuJ3Mohgdosu3LhjDTWOH7I7jlJBTwt6UZYtg06dIFS7EgoyLb47YSabfimL7Y6iVNDTgl6Y1FTYuROuv97uJH5rT+UaLK/Xkn4piwnN1gFRpeykBb0w335r3WtBL9TH8T2oeeIQnXcm2R1FqaCmBb0wy5ZBhQoQH293Er/2dcM2HChflXs26OCoUnbyqKCLSHcR2SYiO0Tk6QKW1xGRb0Rkg4ikiIgzLumzfLl1/Ln2nxcqMzSMWS260nnnOmrzq91xlApaRRZ0EQkF3gZ6AM2A/iLSLN9qzwGzjDEtgbuBwL8CwsGDsHXrBd0t0TFZuSfQ5NyUZUbczQAM4T2bkygVvMI8WKcNsMMYsxNARGYAvYAtedYxQEX340rAXm+GtMWKFdZ9noK+P82acTCv3WOD73J0BUmrVJ1lDVpx/873IXM0hHnyT0sp5U2edLnUAvbkeZ7qfi2vMcC9IpIKLAQeKWhDIjJMRJJEJCk9Pf0y4pag5cuhTBlrylzlkelx3anF3v9de1UpVaK8NSjaH5hijIkBbgE+FJHfbdsYM8kYk2CMSYiKivLSR/vIsmVw7bUQEWF3koDxdcM2pFETJk60O4pSQcmTgp4G1M7zPMb9Wl73A7MAjDHfA6WBat4IaIujR2HjRj1c8RJlhYTyPvfDl1/Crl12x1Eq6HhS0NcCjUSkvohEYA16zs+3zq/ATQAi0hSroPt5n0ohvvsOjNGCfhneY4j7gQ6OKlXSiizoxphM4GFgEbAV62iWzSLyvIj0dK/2ODBURDYC04HBxhjjq9A+t3w5hIdDu3Z2Jwk4e6gDPXrA5MmQkWF3HKWCikeHIhhjFmINduZ9bXSex1uADt6NZqNvv4VrrrEGRdWle+AB6NULPv8ceve2O41SQUPPFM3v7FlYtw46OOf3U4kKzSKs1y2kUosv75ioF71QqgRpQc9v/Xo4f14vN3e5skKJGbmI2R2upxuL6fDA+3rRC6VKiBb0/FautO6vvdbeHAFuZmw3jAj9Ny6yO4pSQUMLen4rV8KVV8IVV9idJKDtqxjFNw1a0/eHrwhDB0eVKgla0PMyxiro2t3iFdPiexB16ii9Qub+bg4c7VdXyvt0wo28fvkFDhzQgu4lSxu0Jq1CFA+ceJekkeUuWKZz4CjlfdpCz+v77617LehekR0Syoy4bnRlCXWPBP58bUr5Oy3oea1caV3QonlznSrXS2bGdiOTUPonf2l3FKUcT7tc8lq5Etq2hdBQ9qehU+V6wcEKkXxGL+76YQmvXTeA82HhdkdSyrG0hZ7jxAlISdHuFh+YwHAizxyn+08r7Y6ilKNpQc+xZg1kZ2tB94GvuYndlWtwT7Jec1QpX9KCnmPlSqujvG1bu5M4jiGEafHdabtnEw0P6TVHlfIVLeg5Vq6E5s2hcmW7kzjSJy26cj4kTFvpSvmQFnSwTihavVqny/Wh38pW4ovGHbhz038pc/6s3XGUciQt6AA7dsCRI9rd4mMftexBxXOn+MPW5XZHUcqRtKCD1ToHLeg+tjamOT9Wq8uADQuAwL3+iVL+Sgs6WEe4lCsHzZrZncTZRPio5S20OPAz17DW7jRKOY4WdLBa6AkJEKrzdvvavOY3cDKiDP/HeLujKOU4WtDPnYPkZGjTxu4kQeFkqbLMa9aZfsyEw4ftjqOUo2hB37jRukKR9p+XmI9a3kIZzsKUKXZHUcpRtKCvWWPdawu9xPxYvT4raM+OJ94hRLJ1jnSlvEQn51q9GqKjISbG7iRB5W0eZjp/ZMBdf2NZg9aATn6mVHFpC33NGqt1rvPjlqg53El6ucoMXP+53VGUcozgLuhHjsBPP2n/uQ0yiGBaXHdu+DmJ2kf32x1HKUfwqKCLSHcR2SYiO0Tk6Yus01dEtojIZhGZ5t2YPrLWOhb6rldb68UsbDAtvjvZIty7YaHdUZRyhCL70EUkFHgb6AqkAmtFZL4xZkuedRoBzwAdjDFHRKS6rwJ71erVIMLiw231YhY2OFChGl9e1Z5+KYt5veMf7Y6jVMDzpIXeBthhjNlpjDkPzAB65VtnKPC2MeYIgDHmoHdj+siaNdCkCcepZHeSoDW19W1UPnuSnlt0fhelisuTgl4L2JPnear7tbyuAq4Ske9EZJWIdC9oQyIyTESSRCQpPT398hJ7izH/GxBVtlkT05ytUfUYvP4/6PwuShWPtwZFw4BGQGegP/CuiFTOv5IxZpIxJsEYkxAVFeWlj75Mqalw8CBcc429OYKdCFNb3Uazg7/QkRUXLCroQt16rLpSF+fJcehpQO08z2Pcr+WVCqw2xmQAv4jIT1gF3n9nYEpKsu4TEuzNoZjbvDNPLfuAEWffAK7LfX1/WqiObSh1CTxpoa8FGolIfRGJAO4G5udbZx5W6xwRqYbVBbPTezF9YN06azKu2Fi7kwS9s+GlmRF3M7czD3btsjuOUgGryIJujMkEHgYWAVuBWcaYzSLyvIj0dK+2CDgsIluAb4AnjTH+PfNSUhJcfTWUKWN3EgVMbXUrBoG337Y7ilIBy6M+dGPMQmPMVcaYK40xf3e/NtoYM9/92Bhj/mKMaWaMaWGMmeHL0MVmjFXQtbvFb+yrGMWn0psjr7xHOTml5wIodRmC80zR3butqVu1oPuVN80IqnCUP3d76nd950qpogVnQc8ZEG3d2t4c6gIrac/GGo0YnDQfMdl2x1Eq4ARnQV+3DsLDdUDU7wiTE3rS8LdUOu1cb3cYpQJOcBb0pCRo0QJKlbI7icpnYZOO7CsfyZC1c+2OolTACb6CrgOifi0jNJwPWv+Bjrs3Ekey3XGUCijBV9B37oSjR7X/3I9Ni+/OqfDS/IXX7I6iVEAJvoK+bp11ry10v3W8dHlmxXalP9O54sQhu+MoFTCCr6AnJUFEhHVSkfJb7yf0IoRsBukVjZTyWHAW9Lg4q6grv5VauQafcgf3bPiCsufP2B1HqYAQXAU9O9vqctH+84DwKo9T6dwp+m9cZHcUpQJCcBX0nTvh+HEt6AFiNe1YVftq7l87j/CsDLvjKOX3gqugb9hg3bdsaW8O5bEJbftQ88Qhem1ZZncUpfxecBX09eshLEwHRAPI0gat2RpVjwdWz9HpAJQqQvAV9KuvJvrKsAuugqP8mAjvtOtDo8N76LJjjd1plPJrwVPQjbG6XFq1yr0STs5N+bcFTa5jT6UreHDVJxCSqZelU+oiPLkEnTOkpUF6uvafB6CskFAmtenNC19NoBMr2DXy1AXL9bJ0SlmCp4W+3j17X6tW9uZQl2VWi66kl6vMc/zN7ihK+a3gKegbNlh/n8fF2Z1EXYZz4aWY2OYOuvA1rdK22h1HKb8UPAV9/Xpo0gTKlbM7ibpMH8ffQjrVeGSlf1/hUCm7BFdB1/7zgHYmojSv8Rdu2LmOFvu22x1HKb8THAU9PR1SU7X/3AHe5iGOli7PI9/PtDuKUn4nOAp6zhmiWtAD3gkq8u/WPem2fRVND+60O45SfiU4CnrOES7x8bbGUN7x74SeHC9Vjj+vmGZ3FKX8SnAU9A0boH59qFLF7iTKC46XLs+kNr3ptn0Vsft+gtAsPdlIKTws6CLSXUS2icgOEXm6kPXuFBEjIv51OaD167W7xWH+3bonv5WpyF++/RiyLjzzt+7IBexPC7U7olIlrsiCLiKhwNtAD6AZ0F9EmhWwXgXgMWC1t0MWy/HjsGOHHuHiMKdKlWVi2zvo/Ms62vOd3XGU8guetNDbADuMMTuNMeeBGUCvAtZ7ARgLnPVivuLbuNG614LuOFNb3kZ6ucq8wCi7oyjlFzwp6LWAPXmep7pfyyUirYDaxphCZ7oSkWEikiQiSenp6Zcc9rIkJ1v3OiDqOGciSjO+3V3cyDe035VsdxylbFfsQVERCQFeAx4val1jzCRjTIIxJiEqKqq4H+2Z5GSIioLo6JL5PFWipsX34FdqM3LZBzpfugp6nhT0NKB2nucx7tdyVACuBpaKyC6gHTDfXwZGUz5MZnF6PBIiOv+5A50Li2AULxC3fzu3/rjC7jhK2cqTgr4WaCQi9UUkArgbmJ+z0BhzzBhTzRhTzxhTD1gF9DTGJPkk8aXIyKBxxiZ2tqmg85872Efcy9aoejy5fKpee1QFtSILujEmE3gYWARsBWYZYzaLyPMi0tPXAYvlxx8pxXm2XHGl3UmUD2UTythOg6l7dD9/TP7S7jhK2cajC1wYYxYCC/O9Nvoi63YufiwvcZ/yv7l6A5uDKF9b2qA139dpwaPfTWfO1TfZHUcpWzj7TNHkZE5Thl+q1rQ7ifI1Ef7Z+U9EnjnO8NWz9exRFZScfQm65GR+oAXZIXrWYDBIib6Kuc06M3TNXF7N+hUzcssFy/VSdcrpnNtCNwaSk0km3u4kqgSN7TSYrJAQXuIpu6MoVeKcW9D37IEjR7SgB5n9FasxoW0f7mI2bX/9we44SpUo5xZ09xmiWtCDz6Q2vfmV2oz++l1CsrXfXAUP5xZ090Whf6CF3UlUCTsbXponeZnmB3fSL+Uru+MoVWKcW9CTk+GqqzhFebuTKBvMoi+ra1/NU8s+oMrpY3bHUapEOLug64RcQUx4ruuDlD9/mmeW/tvuMEqVCGcW9KNHYdcuiIuzO4my0faourx3TW/6/rCEhNTNdsdRyuecWdBz5kDXFnrQ+1f7u0mtGMXfFo0nDJ3nRTmbswu6ttCD3pmI0ri6PECTQ7sZIa/p2aPK0Zx5pujGjVCtms6BrgD4qlE7Fjdqx/Pbx7B2aBV2Vf3f9Vn07FHlJM5tocfF6eTnKtdzXR/kLKUZ+8W/9EIYyrGcV9AzM2HzZu1uURc4WCGSP/M6bVM3c++GhUW/QakA5LyCvn07nD2rBV39zgcMYln9Vjy9dAoxxw7YHUcpr3NeQdcBUXVRwrM3P4wR4aWFb2rXi3IcZxb08HBo2tTuJMoPpVWqzvM3DqX9rykMWTPP7jhKeZUzC3rTphARYXcS5admxXZlUaN2PLl8KrFstDuOUl7jzIKu3S2qMCI83f0RjpStyDT+CGfO2J1IKa9wVEFvHn0A9u7l8Q/jck8cUaogR8pW4skej9GcLYwr+5SebKQcwVEnFtXYvwmA/f3OUrfeAkBPHFEXt7xBa17jz/yF19naqzwLm3QE9N+MClyOaqHHuftDt1avb3MSFSie5kXW12zM2C/epP5vaXbHUapYHFfQD5Svym9lK9kdRQWIDCJ4uNdIMkPCGD/vn5TKOGd3JKUum+MK+tYobZ2rS7O3YnX+fNvjNE3fxQtfvQMYuyMpdVk8Kugi0l1EtonIDhF5uoDlfxGRLSKSIiJfi0hd70ctwvnzNGMLP1avV+IfrQLf0isTeLP93fT9YQmP8i+74yh1WYos6CISCrwN9ACaAf1FpFm+1TYACcaYWGA28JK3gxbpxx+JIIMt2n+uLtMbHf/IokbteI2/wOLFdsdR6pJ50kJvA+wwxuw0xpwHZgC98q5gjPnGGHPa/XQVEOPdmB5ISQHQLhd12YyE8OfbHmczzaFfP/jpJ7sjKXVJPCnotYA9eZ6nul+7mPuBLwpaICLDRCRJRJLS09M9T+mJlBTOEcEvVQuLplThTkeUoSfzISwMbrsNDh2yO5JSHvPqoKiI3AskAC8XtNwYM8kYk2CMSYiKivLmR0NKCptpTmaoow6tVzbYTT2YNw/27LGK+qlTdkdSyiOeFPQ0oHae5zHu1y4gIl2AvwI9jTElf+xXSgopxJb4xyqH6tABpk+HtWut7pfMTLsTKVUkTwr6WqCRiNQXkQjgbmB+3hVEpCUwEauYH/R+zCKkp8O+fVrQlXeEZlnTAPS+neHZb8OCBTBsGGTrdLvKvxXZP2GMyRSRh4FFQCgw2RizWUSeB5KMMfOxuljKA5+INYHKr8aYnj7MfSH3gKhV0PXEEFVMWaHUHWlNHfEltXGNHU3iv5+HMmXgrbd0kiDltzzqcDbGLAQW5nttdJ7HXbyc69K4C/pG4ijHGlujKOcZwxgSnzgNr7xizbX/+uta1JVfcsaZoikpUKMGh/DyQKtSAKHZyCsv8QaPwZtv8nLIU0TX0j515X+ccUhISgrExsJ+u4MoR3J3wbxhulBlyQ6eXP8KZfeehuxxEOKMNpFyhsD/15iZCZs3WwVdKV8SIbHLcCa0uYOHGA8DBkBGht2plMoV+AX9p5/g3Dkt6KpkiPDiDfcxkhdh2jTo3VuPU1d+I2ALenSMdWjZ3c2tAdHYgXrZOVVyXmIkTJwIX3wB110Hqal2R1IqcAv6/jSrX/O6dp+SERLKmSd22h1JBZthw+A//4EdO6BtW1i3zu5EKsgFbEHP0fTgL+yIrE1GaLjdUVQwuuUW+O4763DG666Djz6yO5EKYgFf0Juk7+LHqHp2x1DBrEULWL0arrnGGigdPhzOnrU7lQpCAV3QK505Qc0Th/Qaoqrk5UwPkHOrcQUxOxbDU09ZfesdO+r0u6rEBfRx6E3SdwGwTVvoqqTlmR4gx+6xt8LYsdC+PQweDC1bwssvw4MP6pmlqkQEdAs9p6DrVYqUX+nVCzZtslrpDz0EPXrAr7/anUoFgcAu6Ad/4XCZiqSXq2J3FKUuVKsWfPklT1cax+lFyzlZtxmPy6uESSbRMVl2p1MOFdAFvWn6Lqu7Rf+cVf4gf796iDD22MN0e2Acaxo05VWeIKX6VVyZ9r3dSZVDBWxBF7K56tBuPcJF+Q93v3reG0Bq5Rrc1yeR4bc/Q5XTx1nBdXDnnbB9u82BldMEbEFvwE7KZpxja/V6dkdRqmgifNm4AzcMm8hzvACLFkGzZtYhjrt3251OOUTAFvRYrFP+f4zSAVEVOM6Gl+bvPGedXTp0KEyeDI0awQMPwE4921kVT0AX9CwJYXu12kWvrJS/qVEDxo+Hn3+2CvuUKdCwIdxxB3z7LRhjd0IVgAK6oO+qUpOz4aXtjqLUpck7eFqnNjL+bWqbHfzdPMPhucvg+utJDmnJs5X/BUeO2J1WBZCAPbEolhS26oCoCkQXOSnp3ZHt+TCjFb03L6X/xi/5x/7HOFv1KeZxO9PpT3LNruxOK2tTaBUIArOFfvIkDfmZbVF17U6ilFedDS/N9Pju9Bz0Bi1Zz8xWXehW5gs+43Y27q0JgwbBp5/CyZN2R1V+KDAL+qZNAPyoZ4gqB0umJYldH6TNQ1MZdJeLz6Qnv039D9x5J+cqRLJIbsZV6SVITobsbLvjKj8QmF0uP/wAoF0uKihkhoaxrEFrpprRNHhyPgmpW+iyYzWddq7n5sOLoeVIiIyEDh2s6Qbat7fmkSmr3TPBJjALekoKx6lAWqXqdidRqkRlhYSyuk4LVtdpwd9vhPNjW7L3gyWwdCmsWAHz51srhoRA8+bQurV1ecbYWGua36goPbPawQK2oP9AC4wEZo+RUt6yL/QKZNBAYCAA1TlAj6ormfLQBli7FhYutA6JzFG5MjRpAo0bw5VXQoMGUL8+1KkD0dEQGmrHbigv8aigi0h34E0gFHjPGPNivuWlgKlAa+Aw0M8Ys8u7Ud2MgZQUUrjbJ5tXKqAUcMTMB6/05IMXeuc+r84B4sOSaZq5lcZHt9F41TauWrWEGD644H2ZhHIgtCa1EmpaxT06Gq64wmrVR0VBtWpQtap1q1IFypXT1r6fKbKgi0go8DbQFUgF1orIfGPMljyr3Q8cMcY0FJG7gbFAP18EJjUVjh4lhVifbF6pgFdAkV889la2jcwEGgG3AbB/7E10GDKF2kf3U/PEIaKPH6LSqnCiV+8nmu1Es5xIfrv454SGQqVKULEiVKhg3cqXtwp9zq1MGasvv0wZKF2aZ/8WQfqxMpwngnOU4jwRlKsayofTSkNEBISFWZfzy7kPDbUeh4ZeeAsJ+f19SIj1Cybvfc7j/Ddw5C8jT1robYAdxpidACIyA+gF5C3ovYAx7sezgbdERIzxweluKdYp/z/QAjjm9c0rFSzOUZqfI2vzc+T/zrbe/f2tF/wyCMvK5OQrbYm7bw5Vzhyn0pmTZM+rx7svHYVjx+DoUTh+nC8+OU7E2ROU5yjlSKMcpyjLacpwhrKcJgxryuB/FBTkN6C7L/e0CHkL/MUe51/X08cX8/jj8Pzzl5e3EFJUzRWRPkB3Y8wQ9/MBQFtjzMN51tnkXifV/fxn9zqH8m1rGDDM/bQxsO0yc1cDDhW5lrPoPgcH3efgUJx9rmuMiSpoQYkOihpjJgGTirsdEUkyxiR4IVLA0H0ODrrPwcFX++zJYSJpQN4ZsGLcrxW4joiEAZWwBkeVUkqVEE8K+lqgkYjUF5EI4G5gfr515gOD3I/7AP/1Sf+5Ukqpiyqyy8UYkykiDwOLsA5bnGyM2SwizwNJxpj5wPvAhyKyA2uIw9fHFBa72yYA6T4HB93n4OCTfS5yUFQppVRg0FMtlVLKIbSgK6WUQ/h1QReR7iKyTUR2iMjTBSwvJSIz3ctXi0g9G2J6lQf7PFhE0kUk2X0bYkdObxGRySJy0H0uQ0HLRUT+5f55pIhIq5LO6G0e7HNnETmW5zseXdIZvU1EaovINyKyRUQ2i8hjBazjqO/aw3327ndtjPHLG9YA7M9AAyAC2Ag0y7fO/wET3I/vBmbanbsE9nkw8JbdWb24z9cDrYBNF1l+C/AFIEA7YLXdmUtgnzsDn9ud08v7HA20cj+uAPxUwL9tR33XHu6zV79rf26h5045YIw5D+RMOZBXL8idYWg2cJNIQE/Q4Mk+O4oxZjkUNmEIvYCpxrIKqCwi0SWTzjc82GfHMcbsM8asdz8+AWwFauVbzVHftYf77FX+XNBrAXvyPE/l9z+M3HWMMZlYk7tElkg63/BknwHudP9JOltEahew3Ek8/Zk4zbUislFEvhCR5naH8SZ312hLYHW+RY79rgvZZ/Did+3PBV0V7D9APWNMLPAV5JsDVTnBeqz5OuKAccA8e+N4j4iUB+YAI4wxx+3OUxKK2Gevftf+XNCDccqBIvfZGHPYGHPO/fQ9rDnoncyTfweOYow5bow56X68EAgXkWo2xyo2EQnHKmwfG2M+LWAVx33XRe2zt79rfy7owTjlQJH7nK9PsSdWv5yTzQcGuo+AaAccM8bsszuUL4lIjZyxIBFpg/X/NJAbKrj3531gqzHmtYus5qjv2pN99vZ37beXoDP+OeWAT3m4z4+KSE8gE2ufB9sW2AtEZDrWSH81EUkFEoFwAGPMBGAh1tEPO4DTwJ/sSeo9HuxzH+BBEckEzgB3B3hDBaADMAD4QUSS3a89C9QBx37XnuyzV79rPfVfKaUcwp+7XJRSSl0CLehKKeUQWtCVUsohtKArpZRDaEFXSimH0IKulFIOoQVdKaUc4v8BxA2pVEEoRaQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(speed_list,bins=50,density=1, edgecolor='b' )\n",
    "plt.title('Speed Histrogram', color='green')\n",
    "plt.plot(vs,prob, label='theoritical curve',color='r')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x28cc9ae5090>]"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvr0lEQVR4nO3dd5gUVdYG8Pd0TyANeQiShoyoIDIiSUAwgDmy5lXXRdd116yYP7NrWrOYMK+CCoqCoCAoGYfMkHOGIQ8wqafv90dX1VR1V3WY6Z6xi/f3PDx0qK6+1VN17q1zb90SpRSIiCj5eaq6AEREFB8M6ERELsGATkTkEgzoREQuwYBOROQSKVX1xQ0bNlRZWVlV9fVERElp/vz5e5RSmXbvVVlAz8rKQk5OTlV9PRFRUhKRTU7vMeVCROQSDOhERC7BgE5E5BIM6ERELsGATkTkEgzoREQuwYBOROQSDOhEREFmrduDNbvyq7oYMWNAJ6JKsTu/EFnDx2Pqqt1VXZSIrn5/Ls767++2763dfRiLthyo3AJFiQGdiAx5+UX4fM4mJOLGN8u3HwIAjJyxwfL6o98tw6y1e+L+fYly5iu/4eK3ZlZ1MWwxoBOR4aGxS/HId8uwYkf80w0eEQCAua5QSuGzOZtw9Qdz4/59xyIGdArrYEEJtuw7WtXFoCCfz9mEJVsPxH29eflFAICCktK4r1sP6H5TRPfzDphxxYBOYV345gyc/sLUqi5GpTpS5MOTPyzHviPFVV0UR498twwXvhn/034t5gKIf6T1aOs2t9BLGdHjigGdwtq098/fOl+fdxj/m7s5buubvGIXRs7cgLenro3bOpOFXVpEp5TCws37I+bXJyzdgQWb94e+oQV0cws9UQFdKYXpa/JQ7PMnZP36dzi9/u8vF2Lehn3Ga3d8tRDdnvw5YWXRHfMBfcu+oxids6Wqi0EVcOPHf+ChsUtxtNgXl/UdKigBABT6QtMOs9buQdbw8X/q1ntF6A10uzg7KXcnLnl7FiYu22n72Vlr92Dm2j247YsFuPTtWaELaOu0tNAdguLc9XuRNXw8dh4sjKH0ZRZtOYDrPpyHNxNYKZeUWste6ld49LtlmLdhH8Yt3o6/jpxnvPf9ou3Yf7QkYWXRJX1A91ewhr/uw7m4/5slcQsGbhXrqIfd+YX4ZNbGhIyWCLZtfwEA4GhxnPK+Rq439K0Rv68HgHLlr+dt2GdptQXz+1WF9+eKEiMtEloOfZTKKofx2Vd/MBfXhOnc9GnbphC5hf7pnMCU3498twyDX/095v1IrwhW7TwU0+diURRU4e86VIjP5mzCw98tA2DfIEi0pA7oz/+0Ev1fqlh+d8/hQEvrQDlrz99X55W7FZFMYj01fmjMMjw+Lhe52xN3QOn0U/jCCnbkFfv82J1faLRS7WKI/t4dXy2Kef1D352Noe/Odnz/0e+XYchr0x3fv+jNGRi7cKvltYEvTYtrJSDaFtqtMS0lEC5KSsuXxvD5A5/zx5BDn7xiF1buzI9p/9t9qNAIponM0RcFpXP071q7+zAA+/0n0ZI6oI/4bR227CtAfmH5T2X0jpojRZFb6G/+uibk6rHrR87D+W/MKPf3l9f2AwX4deUu+P2qwmcXhSWl2J0fqJQ27z2K+79ZHJJ79JkOjILiUjzxQ67tbzZ73V78vjoPe48UGetONL1owQdYrIZ/uwQ9npliyvGGHpH6/nKwoCTuweKLuZttW79KKSilsHjrQdw1arElgK/fc8R2RMqCzfsx+NXf8caUNbEVwibPrUv16gG9fNvt0z6nlMKuQ4XYsu9o1Mfui5NWRZ0P7/HsFNw1anHgO2P8Gx0u8qFflIMAgve34igqOl85K8NoJWVAn79pP7KGjzeen/R/PxtBbfLyXcgvLEHXJ37GmAVbnVZh0P/cRyKcrh8t9uGln1fjyvfm4PJ3ZuE/E1ca7+05XBT7Rmhmr9uLd39bF/L6jDV7sHW/c4fktR/OxU0f5+CZCSvQ+bFJKKhAuuHerxcHAplf4d5vFmN0ztaQlIK5Vfb5nE34aOZGjLAp91Xvz8H1I+eVda6Vu1ShIlUO0VYepX5le1Y1KTeQGz5SFFiPOaY9M345Rs7YACkbBmI55S4sKcUtn+Vg2baDUZUhHHPgGvHbOmQ/PdkSRIPzzsGppsKSUlz69iys3JmPl39ZHVPFo2+d/pkiX6nx2KvVZr5yBnR9PQs2H8Bpz07B6S9MxWXvWM9YlAqknYIbC+/+vh5fzgt0fPv9KuoUjP6dvlK/bb/H/iPFlt9nxY5D2BzlMN2ioP0tmt+loo2OSJIyoAefdgLAY9/nYuHm/bj50xzcNWoRDhaU4N6vF0f9hz8aoYV+WHs/v8iHnE378c60dTGd6paU+m2Xv+r9OXjup5Uhr1/74dywwwXX5x0BEAiu5vIBgVO+UX9EP+rjxyU7AATGHusdguOX7sB5r5ed/tsFhXA56+AhauvzDkdMTe0+VIjez03B2t2hrdQt+46i06MTMfoP5w7swhK/9n+pMZ7azkczN6Dnc1OwYU/gN3x18mrMXLvHqIT0gGredd6fvgFP/rjcsj+V+Moe524/iEm5u/DcTyuwYPP+kMplfd5h43HPZ6c4lg2wjgF//qeV2Huk2EhXAKF/C70yz8svwsGCErw4aZXlfbsUib4dL05aib+Y0kD6b6C3bDs+MhF3fLUQQFlAL/VHF5TMDZLdhwpxoCC0NR7cGHp6/Aq0eWgCpq3KC1lW39/6/udXPDR2aVRl0H+rx8fl4pSnfgmphLs99QvemVbWceoxVdiRBAfnaFJRq3bl4/mfVlqO13hKyoCe4gkt9jfzt2JS7i4AwOQVgbki/CpwmW40Nuw9Yvt6qV/htGcn43WbU9fglpJSCrd8loPJy3cZry3acgAjZ2xA+4d/wmPjljl+v13FoxQidizqO5H5gD//jel44NuluO/rxfhmfmjlt3Dzftuzl8KSUuQXBna0j2ZutOS/zS3EVG9gpz9aXIrNDsMa9Zasfoo58OXf0Pt5+0B2/zeLkTV8PCbm7sT2g4X4aObGkGXWaEF+Ym7oCAu98li05QCUUnhwzFL0fM45aE5eEfj77DhQgLz8Irw6eY2lM0//TZXN+YW5EVZUWhYc9Bi7Ie8ILn17Fu7/ZonlcwNfLtsPdx4KX7EFt/wAa+UZnA5ZufMQ/H6FU5+ZjN7PTQlpiQYHnmXbDqL1gxMwa90evDV1HeaaOmr1Q6u0tKwV/OOSHdh/pBgvaRVFtGmMvv8pa5D0eHYKHhwTPggX+/z4MGhaADMt44PtBwvx5bwt2Hu4CGe8NM3IWdvx+RVe+XkVvtCGtR40VSr6vv7Gr2UB3S6eHzhajFF/bA5pkJUnoF/69iyM+G0dXvl5dcRlyyNJA7p9LWqXAliXZx+oDdrf6OGxgWC7eMsByzDGgpJS7DpUhM/nbA757uCW0tHiUkzK3YWbP80xXrv4rZl48sflAGCsw86pz0y2DdyPj8vF53M2heTexMh1Bv4fu3Cb8Z7eUv16/lbc+/XikNPXS96ehbtHL8Y1H8yxbENBSallhzfTK4ydBwuN4VdfztuMfi9Otb2SVP+Vinx+o8XrFAdG5wQqF33zzQfVxGU7MGf93rL3TJ8b8tp0SwB46sfl+HzOJoxduA2lfoVinx8Hj5bgcJEPT/yQa7SK9HI8+v0yazDQVq4fmHbltbTQTdFd/43179Bb6Vv2HcU9oxeHrCfc2Z1dTvyQ6e8SvN8N+2y+EZSOFJcaFa4uOPe8UBsjPm7R9pDv0TtFfX5lCVhPj19hpCVX78pH7vboUktPa/t+NF76eVXY939cssNy5vP2tHXYsOcI3tdGHgGhDaN5G/bhdVPANqcm9bRikc8PpRTuHrUIM9eEzinz/vT1eODbpcjZZB1bH5JyieGMvSL9fuGkJGStCeb1Rn9aBARantVSvbbv5ZuC3ZEiHy7SJt0Zmt0CB44W44Pp1haD+Y9mrpE37T2C/i9Oi6lcZnsOF6PI50e1VG/Iwf7o97nYe6QYd57ZwXhNYM1PvzBxFW4b0M523YcKS1AzPfRPPXPtXsspb2GJH0ccOljzC33wegptW7678wvRon4Ny2v6qWuRz4/Vpo6+wpJSpHk92H6wAM3rWT/zlZZOSfOW/a1u/XwBAOD967MBwJLDXrHjEJ76cXkgFaAdyMtNc5AcOFqMHqb0Rq30FNxzdkfjh1uXd8RIuwBllUWxEdAD+VyPqRI3xwtzoNTz7now9ojgryPnWVq/ZkdLSvHI2KVoUqc6Fm7ej//9vafpN/Lj7lGLcPYJjY3XzKfo//0l9GxxxtrQFIVRzlI/fs7diWGfzcf8R85EnRppAKyt1cKSUhSX+jFDmyTL5/dbGgLmQPrHxv047/UZ2Pj8ecZryx1GM30wYwMeOb+zY9nM3jMFZjtLth7EsxNWGM/1yrxEa2yMztkScmYUTD/TWbHjEP72SVnD6+ZPcjBlpf0skMu2BbbtwNFirDOlzswV3l/enY2NDmf5dhKVS0/KgB5LngsAth0oQNvMWgCA/0xciV2HCvHK0JNDWqOHTLWmUgovTFoVcgWiuQVgvnjC6WILO/uPFGP+pv04s3Njy+uDXv4NN5/eGlf1aBnymVcnr8HmvUfRJrMmRMSx9Sg2v82hAh+a1rEvi7kDKFyn4j+/WIBNDp1FdmXRW5HFpX68ZkpXZT89GU9dfALuGrUYY2/rjW4t6xnvrdgROHD0FubsdXuN9/Sc+L4jRSgsKTVGXJi/K6Ds8bYDBZYy2Y2D3muq0MRUCQHAmAXbMGbBNkvgMjNX6Hrw01vtHoFjMAcCY5a/M7WQ9VFBQKADfszCbRhjOus6VFAWXEfODE1LmP/u+hmPrtjnN1IOS7aWtax/Mu2znR6diKwGZRXslBW7LRWuXfopa/h49GrTALPX7w15zyyeo4GWhul0fjlCCx8oC+gb91iDr1MwL/b5jZTesxNWYKMpxWgOyuH+1nZ8UfZDxCopUy4OGRdH+oUnAPDOtHUYsyBwoOwOymWaa/cin9925Eh6SlnrcY3pdN1pl03zhv7ET/yQi5s/zcGOg9aAs+1AAZ74YbljTnDMwm146efVIZ1e5jLbGZ2zBU85nPpeMaKsQ6zI53esLNfsPuw4bGz7gQJkDR9v6bDUO5+KfX7UMp0dHC7y4XstkG3YY9+i8SuF3YcKcdX7c4zX9E6wBZsPoNOjE3HDR/NsP/vlvLIyBI9W+Hz2Jm39Za/tsQT0wP/Blfi1pvy6eVih/nv4Sv0hnVx2FatZ8L5lvg7CrsPscFH4U/Rwx0Sxz490bQz5jR//4TgAwBysxi7chsveKWuwOHXjRArmAOJ60d7CzQdCXttzuBjnvzEdew9Hvnr3izmb8Pa0tVENMQSADo/8ZPSbbAzqL3pn2tqwF4oBgWHAz4wPPfbKO1IokqRsoXtNB8sNvbOw42CB0SFq56dlO3FqVn1UTysLxt8t3IZGtdMBAI0y0rE7vwjTTfmzoqBApHM6boJ3+Jlr96BPu4ZIT/VYdh7zcMsdDqM+ytsDfv2H8zD61l4hr+unpo+e3xm/rXY+NS/ylTpuXzj66fb931orRCAQTMy/OwBjBEORz287zK/Y57ecWtuZbpPrDBbcQs8v8uFwkc+SNsuLYsjpDNNc3ebRM/rf9bIRs7E46IYHkX7H4Mpml6lxYdePoXfgOZEw37j3cJHl+7buL3Bc1szcso4mcDuJ9vvK6/cw+3Qw/aynY+OMqD/jlO9esPkArv1wLi7v3tzxs/1etB+pFuv4+GhFDOgi0gLApwAaI9AQfU8p9VrQMgMAfA9APxcco5R6Mq4ltX6f8Tg91YMaaeE348t5m7HvSBHeuvoU47U7Ry3CRScfBwAY1q8Nnh5vDSBFJaW2eed8h2BrHpcOANd8MBcnNasT9kB0Glr3q8PpXyTzNoZvLXw2eyMe/T7X8f3CktJy7WjbbSqmlTsDLdliX6lji/DT2ZuMNItZQUlp2Ao6WuttOsTz8osswTeaVp2TDXlH8HXOlpBgDtiPljC77YsFlufmTm27gB6pkre7EEj32Lhcy1lfuCGdTsp7JTWAsFe/VhWn6QvshLvosNjnL9fEcL3aNIj5M9GIpoXuA3CPUmqBiGQAmC8ivyilgs8jpiulzo9/EUOZD5ZqKV7bDs/3ruuOk1vUNTrFJuXuQruHf7IsowfOzIz0kM8X+fyomWbfkRqtcPk+wDnnF6lzKBxz3jlYuGAOlAXhWIW792Jxqd9xvLrdeHMAKCjxO45kioXdkM0zXppmeW4+S4o1aL3z2zrH9Fis/Tx6GhAAVtv8HezGZZuFC1DBZVycgHnU3Wz1LudhkeV1VlD/WbxEzKErpXYopRZoj/MBrADQLCGliZJ5KGJ6qgdpNqNeaqSloFHtaqhTPdVxPXrruVFGtZD3inylCZgR2ioRO4o57xyrFyZG7lSyszpcQPc5B3SnS8gLikstI0sSKdqrAu2EG/9ckVEMn2i5frNwqTIgtsqovBU3xU9KjCP1ohVTp6iIZAHoBsBuSrVeIrJYRH4SkRMcPj9MRHJEJCcvL/q8l9mMNXvww+Ky0QHpKV60a1QrZLnqaYFNe/XKkyOuMzMjLeS1wLhbzsAYzC7OhsvSFPmch0KaNahZ9jcoLCmNueP7z6Y8FUW9Gs6ND3KXVJvBEvEQ9VpFpBaAbwHcqZQKTnwuANBKKdUVwBsAvrNbh1LqPaVUtlIqOzMzs1wFDr6iMz3Fg2tOa4X/3Xya5XU9DXNGx0YR11m7WuiBNG1VHt79LTT10a1lXbxmU0kMzXbuGGlWt7oxjjrZ1a8Zmp4Kp9hhtFCw7q3Khi/OWLunwnNHX3lqiwp9vio0q1e9qotAlSQeKUU7UQV0EUlFIJh/oZQaE/y+UuqQUuqw9ngCgFQRaRjXkmqCf4g0rwcej6B3O+vXReooNUt3uOjIztjb+tjm3FsGXVhjdmKz2rafiVV7mzORyhZpP3z2kpMsz9/9fT32RnEziHgHsyZ1rGm04UM6xXX9idCktrXM4fYpqnw9surHbV0pVdVCl8CQkg8BrFBKveKwTBNtOYhID2295R/nFEZwQDf3/k+5p7/x2Jw7f+HyLmHXWTPNi7G39Y66DHY595YNalqe925b1ot94GgJaler+AjRWjGso5V2kYg+h3XzOAVMp86+fw9sh7kPDULPNqE7fV5+UcTvr18jNO01NLs5nrrINnsXUY2gDu1EpTNqpafgS9NVnhXRsJa10r+ga9OQZaqlViwQ5D5xToU+7+SCrsclZL2JcH2vVuX63Khb7P/O5Qn0wdMzxEs0e0cfANcBGCgii7R/54rIrSJyq7bM5QCWichiAK8DuFIl6FY1wbkn89V6bTNr4aMbT0XPNvVR1xTQh2a3cLza796zOyDF60Gj2qFB2ok+ft2sf4dM/PcvXY0LOM7rUnYw7jxUiNo2nbMNa1mDWL0aqTjhuNrG809v6mF5P8MmNeREHxbVSmvl3TGofVx2Iq9DE71to1poXLtayJhzXYugy/zrm3LmV/VoaXuR0eXdW6BVUEUZreys+mhYK82o2J2mfqiogpJS9GhdHw8O6YR/D7SfeiFa5t8EABrYpLcm/Pv0Cn2H3VBcJ2ceH91IjN5tG1TJ2eMZHTPxjc11F5FceopzetTJfed0dLxYLPjvFg27CQbjIZpRLjOUUqKU6qKUOln7N0EpNUIpNUJb5k2l1AlKqa5KqZ5KKZsbCsZHcO9wz6DxnGd0bISvhvWKepREB+0CAz0QRyMj6KD4399PQ53qqbikW3NjbLtHBENObAIAKCrxI8OmdW2+6hQAZg0fhLG39TGe10y3vl/d1Dqbck//sGkE/WAcdHxjbHz+PFyR3QK/3XcG/hVF0Ll/cEfH94L3Q/1nNgJnin3gDG6h33x6a+Pxc5eeZDsSJqthDXRvVc+SD4+2ddUoIx05j5yFvu0DqbhwV24O6hS5n8VJqV/B6xHc0r8tmtukSEZce4rNp+wFn1U0qBUaKNJi2E8r6qmLT8BVPSL3RSjlXNHbWfDoWRUplsHr8eD4poEGULhW8k19AvtaRnoKNj5/Hjo3re24rJP+HZz7/P5m2pcD5Yr8W1RlC/1PJbiF3rVF3QqtT78aVD+YOjSO3NIQETxy3vHGc3Ntq58weAR4YHAg4BaX+kOCNxCa562e5rUcsCkej+WKNr2Ved85HdE2sxZu7d/W8nk9rdCxcQbO7NwYk+/uj7vOam+8f1zd6mGvatOFOxMITrlkaS1o/TPmlvDIG8o6goMn7+rUxHqlXqppu2cOH4g3ruqGRhnVUDM9Bc9fVpYy69K8rvH4pSu6AoDtKCe9HPpVxX6/ckwLvG0Kum0zrWcEgzo1Ql+tf+bNq7uFfNa8HdVtzgJ6tXHuSprz4CDL8xSvB9ecVjaPj56CqZnmNR6np3htp5MAgFv7t8VPd0RuwX99ay98+4/wKcZOTTLQsFa6cQW0XYNEF+tosHAt2lgGD6R4BDXTU/Djv/rioxtPdVzuuLqB40yfOC3aSvEc0+RoTg2+5y49Cadm1ccrQwP7YpPa1YxGTmObM3ldpKkhyivpArrTzhyN4+qEplVqGgE9BWNu642Pb+wRsgwQ6JA053NvPr2NMZmR+azhCm20y8kt6qGelhdu3dA+bfDm1d1wQ++skNf12rtaqtdygDptu955JiKY99AgfPfPQCu/XaNaIRVJqwY1Me3eAbbrCadF/UALOzigvzS0K5686ASc0rIuAOuO3699WatGryhPa10fP/6rL05rHTiz0juL7zqzrOJpVre6Y/C9+OTjjIqya/M6+Pmufhh9S+hpt14OvbXkVwpvXNUNfdqFXqFn/o2m3DPA8t7Vp7VEHa2iLPUrzBo+0Hjv6YtPxMvagWz+LjOnfo/re7VCkzrV8H8XdEY/rfVXKz0FT198orGM/tu0zqyJT246Fbf0a4OGtdIcA9LwIZ1wfNPa+CJoxFfwmdypWfVxstYQMvf1vGLalol39kOq12ME9IfPPR5Olmw9aHtl7F+ynVv3Tn0akS64MVeg+qyrJzarEzaVpP+OdldBm0dXBe9Hp5v2X6ff/KRmgVnv9H2oZf0axjHy1EUn4qY+rfHMJSfafjYRkm4ul4oMyP/+9r6YvGKXZaL9uqYd65SW9Rwv8x11S6+QloWe1kk1tdB7tmmAdc+eaxzcL1/RFX2CRuD8Y0BbnN25MZrWqY6bT2+Nj2dttLyf6vWgpLQU1VO9ltRRupZyMU+vu/rpIdhzuAi9n/8VAKLqC4iUT840neqPu70PfH6FRhnp6PufqejXvqEl331cneo4xTRjorm8KV4PNjx3LoCy+VBSvIITtYPg3rM7YGCnwAHcJrMW2mTWxIAO4dMfKV4P/jGgLYb1axP21FZvLesHlz4vSXCF9MlN9hW4rmZ6Cm7snYXxS3agZ5sGaGz6fa/taU3/BF9+371VvZAytm5YE1PvHWDM2nlDn9a4rHtzjPhtHS7v3hwigrQUD4p9fjTKSMfb15yC7Kx6aJRRDSccF/jd0lI8gHb1/guXdbHMoQPAsr9tfP48KKXwfNBdsbwewfT7z0C1VC9OfWYygEBu+e6gudv1bTL/bkNObILpa/bgpSu64tbP5yM9xWOZXynVK3jm4pNwySnNMOSkJvhs9iZc07Mlbvq4bLra167shqXbDjpONOfkk5t64Jv5W/HipFW2Q//q1khFs7rVLTdn0WdaNbeYp947AHWrp6JezTSMztmCjo0z0LVFXbx5dTfc/j/rHZqAsmOmQc00y6gt/cxFb0CIlP1WjWpXw2MXdMbRYh/2HS7Gy78k5qYWZkkX0CsyID8zI90yFCzVK5YDFHA+tbKrSPSdOPg9845wmSnF8fNd/VAjzWuZltRue/QdNT1oRIOenzY3NNJSPMb3141yJIe5tXFKy7q4oU9rNKldDXVrpGLnwUKcrLW2AWuK45e7+qFNZi3LlYx2V+J2apKBS08JXEysn1rqv4m5NXz7wPaWz/0a1DoOJ1KeUh8Wpp9ZNNRaac9f1gVvTV1rzL8RLjcKBFrNJzar49ipbqZXGs3rVcfW/QWWU3ad3jI2n3JnVEvFfeeUtaJfuKwLXpi4EhnVUnHuSaEjXczbXrt6Km7p1ybsNBNOp/ct6tcwyny2Q8vY2NVMq3jn2u7G48fO74yebRpg1rqyCcxSPB4M1fo9BnRshAEdG4U0lPp1yES/DplGQP/un30cR0L9+K++xo3YM2ul4/perbB060Hce7a1r+ePh89ErfQUTFm5ywjKQOAYG3tbbzStU7Z+81nzUNOZxPldjsOYBdvw68rdliG6+rH3wOBOuP/bJbikWzMM69fG6LTXjymPiPE5PT7USEvBvwa1Z0C3Y+5MmP3gwDBL2ksznYqvfnpIyM6uB4LBJzSx3O6shk2rVj+woh3P08Fmhje7gP7kRSfi0e+XGSkbnR7gg+emzqyVjvsHd8QFXaIbOmautMaYOmH1Mprvu2jWPqj8aV6P7TC6iXf2C3mtR1Z93Ny3NYb1bxNVGYM9en5n1K/pXGG1alADm2xuh3fbgHbo3LQ2BmiBu1nd6nj2kpOinlApluGefds1RJPa1fDudd3hEbGd0S+awV8Xd2uGi7s5z65hTr2legUPhkmHROL1CGY8cEbIkEmd3kJ3qj5v6hvoEJy3oWyUsl39Yde/YJaZke5YBv2MDgicAWZUS8WI67qHLKenVoLPwlK8Ypl3PxL98yKCmmleHCkuNY49vfGklDI6ZIGy21GmeMU4Sw0eQDD7wYHo9dyvUZejPJIuoJv/WOYaN1r6wSBwbrn88fCZqFcj1ZjMa9ztfWwvBNDLEm6mu0jserudDmivtocEf52ION6tyI5eqQWP1jHej/IsaPUzQ6L+zhSvJ+o719j5W9/WYd//6IZTjft2Tr677HqEtBQPzj6hScjyqV6xjJC6sU+WMb/5OSc0xqTcXZg1fCDq2oyPd9KgVjrmPDQo5PV+HTKNKV7jMZbXvM/YTfEMBAJoZ9MQ2Bcu7xLS4asznzH2aF3fUokN6tQYYxZsQ9cWdTGgY6YlvWZmTvN4bY6rSKPO4jnqIzigx9rvpp8he0Tw/l+zMXn5biPloseB4Hy8/p1dmtcxbiISfBbZtE51fHB9tu3NsuMl6QJ6RQ8IvaYNNxte8FWd5rSDWdld0MtfqpppKaiVnoIHwgxBfP2qbpiyYpdxKlfR3yDV68H9gzs65qtFBAM7NTKGYCaDNpllI13sRr0EW/PMuZbnj19Q1uE94truUCpyEIrWpzf1wLjF2/HvLxfGZfyxflZ3Y58sZDsM1wu+gGhomA5Ks+COwfO6NMWg4wejWqrXccAAEDh7m3JPfwx6+bdy/W7moJuRnuI4TXU0gg/tWMvjNQI60LttQ/RuW1ZZtdAqu+ODhj72bdcQD597PP7So4Vx9mdXsQXfpSzeki6gV5Sx40T5Nw43XOuR847HA98uQfsohjo68XgEyyJcvXdh1+NwYdfj8Lp+K7c4XLMVqUU/8gbnYWBuJyIR5zOP1ZATm+CmPq1xewUvPgKAvu0bYs3uw7j5dOeO4XjOVhntRVn6GV8sY9J15jPgqfcNQEFxKU5/wf7mEBHXFfT9dunScDymFnqwbi3r4cd/9Q0J6F6P4O/92hiPzeupTMdeQNfSDalR/NizHxyIGqnOP9FpbRpg2n1nxK1skehD/7IchkFWliu6N8ecDQmZ2aFCRt/SC/uPlv+GFYmU6vXgsQvKn3Iye2BwJ/zl1BZoVvfPNZmXfpXwJWHy/07MKRenXHq0+nfIxG0D2uLanq2w/2gx6sV4JadeFKdUqjmnb8cYCJCgsebhHLMBvWEUk2WVJ0efSOec0ATf/qOXYx6zsrx4RdfIC1WBHq3jN3lSOA8M7hRVWidRqqV60alJ7Fc7JlpGtVRMvrt/TBXNC5d3waezN9peeFdeKV4P7teuVTiuHJWe3ldV3tvEGaNc2EJPvIY103FVj5bGBUDJRETQvVXlBC1y9o8BbSMvdIyKtaIbmt0i6vx+ZTEa1uXMbHqNUTLxKU8sjrmA7vEInrv0pMgLElFcnda6PvZEcVNuO1d0b441Ye4QFU+eCCmXSPRPJWZ6wvCOuYBORFVjlM0UDdGqzDRf2XDk8n3+6h4t8fIvqyvcF1AeSRvQP/yrO+4ARER/LlLB60tuH9gOt53Rjjn0WFTFj0VElWfx42dHdWVtvBnXe5Tzu0UECZodN6KkDehE5G528wRVhoqmXKpS0k2fS0SUSFLBFnpVSrqAnoS/MRElEf36k1gvSPozSNqUS6Lu+EFEx7a/n94azetVx/ldQqcu/rNL2oBORJQIKV6P4x2z/uySLuVCRET2GNCJiFyCAZ2IyCUY0ImIXIIBnYjIJRjQiYhcggGdiMglIgZ0EWkhIlNFZLmI5IrIHTbLiIi8LiJrRWSJiJySmOISEZGTaC4s8gG4Rym1QEQyAMwXkV+UUstNywwB0F77dxqAd7T/iYiokkRsoSuldiilFmiP8wGsABB8F9iLAHyqAuYAqCsiCbpulpO5EBHZiSmHLiJZALoBmBv0VjMAW0zPtyI06ENEholIjojk5OXlxVjUoHVV6NNERO4TdUAXkVoAvgVwp1LqUHm+TCn1nlIqWymVnZmZWZ5VEBGRg6gCuoikIhDMv1BKjbFZZBsA8627m2uvERFRJYlmlIsA+BDACqXUKw6LjQNwvTbapSeAg0qpHXEsJxERRRDNKJc+AK4DsFREFmmvPQSgJQAopUYAmADgXABrARwFcGPcS0pERGFFDOhKqRmI0AepAvdq+me8CkVERLHjlaJERC7BgE5E5BIM6ERELsGATkTkEgzoREQukXQBXXEqFyIiW0kX0HXCyVyIiCySNqATEZEVAzoRkUswoBMRuQQDOhGRSzCgExG5BAM6EZFLMKATEbkEAzoRkUswoBMRuQQDOhGRSyRdQOdULkRE9pIuoOsk/F3xiIiOOUkb0ImIyIoBnYjIJRjQiYhcggGdiMglGNCJiFyCAZ2IyCUY0ImIXCJiQBeRkSKyW0SWObw/QEQOisgi7d9j8S8mERFFkhLFMh8DeBPAp2GWma6UOj8uJSIionKJ2EJXSv0OYF8llIWIiCogXjn0XiKyWER+EpET4rROW4qTuRAR2Yom5RLJAgCtlFKHReRcAN8BaG+3oIgMAzAMAFq2bFmhLxVO5UJEZFHhFrpS6pBS6rD2eAKAVBFp6LDse0qpbKVUdmZmZkW/moiITCoc0EWkiUigvSwiPbR17q3oeomIKDYRUy4i8iWAAQAaishWAI8DSAUApdQIAJcD+IeI+AAUALhSKWa6iYgqW8SArpS6KsL7byIwrJGIiKoQrxQlInIJBnQiIpdgQCcicgkGdCIil2BAJyJyiaQL6BwRSURkL+kCuo5X/hMRWSVtQCciIisGdCIil2BAJyJyCQZ0IiKXYEAnInIJBnQiIpdgQCcicgkGdCIil2BAJyJyCQZ0IiKXSLqAzplciIjsJV1AN3AyFyIii+QN6EREZMGATkTkEgzoREQuwYBOROQSDOhERC7BgE5E5BIM6ERELsGATkTkEhEDuoiMFJHdIrLM4X0RkddFZK2ILBGRU+JfTCIiiiSaFvrHAAaHeX8IgPbav2EA3ql4sYiIKFYRA7pS6ncA+8IschGAT1XAHAB1RaRpvAoYWp5ErZmIKLnFI4feDMAW0/Ot2mshRGSYiOSISE5eXl6FvlQ4mQsRkUWldooqpd5TSmUrpbIzMzMr86uJiFwvHgF9G4AWpufNtdeIiKgSxSOgjwNwvTbapSeAg0qpHXFYLxERxSAl0gIi8iWAAQAaishWAI8DSAUApdQIABMAnAtgLYCjAG5MVGGJiMhZxICulLoqwvsKwD/jViIiIioXXilKROQSDOhERC7BgE5E5BIM6ERELsGATkTkEkkX0BU4mQsRkZ2kC+g64VQuREQWSRvQiYjIigGdiMglGNCJiFyCAZ2IyCUY0ImIXIIBnYjIJRjQiYhcggGdiMglGNCJiFwi+QI6r/wnIrKVfAFdwyv/iYiskjagExGRFQM6EZFLMKATEbkEAzoRkUswoBMRuQQDOhGRSzCgExG5BAM6EZFLMKATEblEVAFdRAaLyCoRWSsiw23ev0FE8kRkkfbv5vgXlYiIwkmJtICIeAG8BeAsAFsB/CEi45RSy4MWHaWUuj0BZbTgVC5ERPaiaaH3ALBWKbVeKVUM4CsAFyW2WJGJcDYXIiKzaAJ6MwBbTM+3aq8Fu0xElojINyLSwm5FIjJMRHJEJCcvL68cxSUiIifx6hT9AUCWUqoLgF8AfGK3kFLqPaVUtlIqOzMzM05fTUREQHQBfRsAc4u7ufaaQSm1VylVpD39AED3+BSPiIiiFU1A/wNAexFpLSJpAK4EMM68gIg0NT29EMCK+BWRiIiiEXGUi1LKJyK3A5gEwAtgpFIqV0SeBJCjlBoH4N8iciEAH4B9AG5IYJmJiMhGxIAOAEqpCQAmBL32mOnxgwAejG/RiIgoFrxSlIjIJRjQiYhcggGdiMglGNCJiFwi6QK64mQuRES2ki6g6ziVCxGRVdIGdCIismJAJyJyCQZ0IiKXYEAnInIJBnQiIpdgQCcicgkGdCIil2BAJyJyCQZ0IiKXYEAnInKJpAvoCpzMhYjITtIFdB2nciEiskragE5ERFYM6ERELsGATkTkEgzoREQuwYBOROQSDOhERC7BgE5E5BIM6ERELhFVQBeRwSKySkTWishwm/fTRWSU9v5cEcmKe0mJiCisiAFdRLwA3gIwBEBnAFeJSOegxf4GYL9Sqh2A/wL4T7wLqtt/tEQrV6K+gYgoOUXTQu8BYK1Sar1SqhjAVwAuClrmIgCfaI+/ATBIJDEhN6tBDZzeviGOb1o7EasnIkpa0QT0ZgC2mJ5v1V6zXUYp5QNwEECD4BWJyDARyRGRnLy8vHIVuEvzuvjsb6ehRlpKuT5PRORWldopqpR6TymVrZTKzszMrMyvJiJyvWgC+jYALUzPm2uv2S4jIikA6gDYG48CEhFRdKIJ6H8AaC8irUUkDcCVAMYFLTMOwF+1x5cD+FUpxYnLiYgqUcREtFLKJyK3A5gEwAtgpFIqV0SeBJCjlBoH4EMAn4nIWgD7EAj6RERUiaLqWVRKTQAwIei1x0yPCwFcEd+iERFRLHilKBGRSzCgExG5BAM6EZFLSFUNRhGRPACbyvnxhgD2xLE4yYDbfGzgNh8bKrLNrZRSthfyVFlArwgRyVFKZVd1OSoTt/nYwG0+NiRqm5lyISJyCQZ0IiKXSNaA/l5VF6AKcJuPDdzmY0NCtjkpc+hERBQqWVvoREQUhAGdiMglki6gR7q/qRuIyEgR2S0iy0yv1ReRX0RkjfZ/vaosY7yJSAsRmSoiy0UkV0Tu0F537XaLSDURmScii7VtfkJ7vbV2b9612r1606q6rPEkIl4RWSgiP2rP3b69G0VkqYgsEpEc7bWE7NdJFdCjvL+pG3wMYHDQa8MBTFFKtQcwRXvuJj4A9yilOgPoCeCf2t/WzdtdBGCgUqorgJMBDBaRngjck/e/2j169yNwz143uQPACtNzt28vAJyhlDrZNPY8Ift1UgV0RHd/06SnlPodgWmIzcz3bf0EwMWVWaZEU0rtUEot0B7nI3DAN4OLt1sFHNaepmr/FICBCNybF3DZNotIcwDnAfhAey5w8faGkZD9OtkCejT3N3WrxkqpHdrjnQAaV2VhEklEsgB0AzAXLt9uLf2wCMBuAL8AWAfggHZvXsB9+/irAO4H4NeeN4C7txcIVNI/i8h8ERmmvZaQ/Zp3Wk5CSiklIq4cbyoitQB8C+BOpdShQAMuwI3brZQqBXCyiNQFMBZAp6otUeKIyPkAdiul5ovIgCouTmXqq5TaJiKNAPwiIivNb8Zzv062Fno09zd1q10i0hQAtP93V3F54k5EUhEI5l8opcZoL7t+uwFAKXUAwFQAvQDU1e7NC7hrH+8D4EIR2YhAunQggNfg3u0FACiltmn/70ag0u6BBO3XyRbQo7m/qVuZ79v6VwDfV2FZ4k7LpX4IYIVS6hXTW67dbhHJ1FrmEJHqAM5CoO9gKgL35gVctM1KqQeVUs2VUlkIHLu/KqWugUu3FwBEpKaIZOiPAZwNYBkStF8n3ZWiInIuAnk4/f6mz1RtieJPRL4EMACBKTZ3AXgcwHcARgNoicC0w0OVUsEdp0lLRPoCmA5gKcryqw8hkEd35XaLSBcEOsS8CDSuRiulnhSRNgi0YOsDWAjgWqVUUdWVNP60lMu9Sqnz3by92raN1Z6mAPifUuoZEWmABOzXSRfQiYjIXrKlXIiIyAEDOhGRSzCgExG5BAM6EZFLMKATEbkEAzoRkUswoBMRucT/A8FMLDxo38fRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Time_List,KineticEnergy_List)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "Lennard_Jones_Fluid.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}