{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Harmonic Oscillator using HMC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import the Library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''import section'''\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math as  math\n",
    "import random as random\n",
    "import seaborn as sns\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to apply HMC to collection of 100 independent Harmonic Oscillator to get equilibrium configuration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hamiltonian of Harmonic Oscillator in 1D is:\n",
    "\n",
    "\\\\( H = \\frac{1}{2} p^{2} + \\frac{1}{2}q^{2}\\\\)  with \\\\( m = 1,k = 1\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function calculates the  total Hamiltonian of the configuration\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hamiltonian(x,p,np):\n",
    "    '''x,p: x and p are list of position and momentum'''\n",
    "    '''np : number of particles in the system '''\n",
    "    H = 0.0\n",
    "    for k in range(np):\n",
    "            H = H + ((x[k]*x[k])/2.0 + (p[k]*p[k])/2.0 )\n",
    "    return H      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating Random Momentum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to generate random momentum we use \"random.gauss\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drawp(np):\n",
    "    '''this function returns a list of random numbers'''\n",
    "    t = [0.0 for k in range(np)]\n",
    "    for k in range(np):\n",
    "        r = random.gauss(0.0,1.0)\n",
    "        t[k] = r\n",
    "    return(t)       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can check whether the generated numbers are normally distributed or not by doing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAFoCAYAAABgwz7vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dbWxTZ57+8ctxnpoENoXaCeIN2mUE2oEMS6MlE6G02oEEAm5aNp3loZt2C9kCU2UmUqkiWg0PbTUdNCXsztBqStvtbkkqukOcEi1N05nRdLVKNNuk2ha2w+7wYt6U4rgJMLFrg2N7X/DHf9wEnBDfPnb8/UhInN85h/MzB4dL57592xaNRqMCAACAETlWNwAAADCbEbYAAAAMImwBAAAYRNgCAAAwiLAFAABgEGELAADAIMIWAACAQblWN3A7ly75FYmkZhmw+fNLNDLiS8m1YB3uc/bgXmcH7nN2SPf7nJNj0913F99yf1qHrUgkmrKwdeN6mP24z9mDe50duM/ZIZPvM8OIAAAABhG2AAAADCJsAQAAGETYAgAAMIiwBQAAYBBhCwAAwCDCFgAAgEGELQAAAIMIWwAAAAZNKWz19PSovr5etbW16ujouOVxv/nNb/RXf/VXse0//vGP+vu//3utX79e27Ztk9frnXnHAAAAGSRh2PJ4PGpvb1dnZ6e6u7t14sQJnT9/fsJxX375pX784x/H1Y4cOaLKykq99957evjhh/XCCy8kr3MAAIAMkDBs9ff3q6qqSqWlpSoqKlJdXZ16e3snHPfss8/qySefjKv95je/kcvlkiRt3LhR//7v/65QKJSk1gEAANJfwrA1PDwsh8MR23Y6nfJ4PHHH/Mu//Iv+/M//XN/61rdueW5ubq5KSko0OjqajL4BAAAyQm6iAyKRiGw2W2w7Go3Gbf/v//6v+vr69Oabb+rixYu3/bOi0ahycqY+J3/+/JIpH5sMDseclF4P1uA+Zw97SUh/vPrHuNrcgrmad9c8izqCCbyns0Mm3+eEYau8vFyDg4Oxba/XK6fTGdvu7e2V1+vVX//1XysUCml4eFhbt25VZ2ennE6nvvzyS5WXl2t8fFx+v1+lpaVTbm5kxKdIJDrNl3RnHI458nrHUnItWIf7nD0cjjkavjKi14feiKtvv/dxhX15FnWFZOM9nR3S/T7n5Nhu+4Ao4WOm6upqDQwMaHR0VIFAQH19faqpqYntb2lp0fvvv693331Xr776qpxOpzo7OyVJ9913n7q7uyVJp0+fVmVlpfLy+CEHAACyR8KwVVZWptbWVjU1NenBBx/Uxo0bVVFRoebmZp05c+a2537/+9/Xf/3Xf2nDhg3q7OzUD3/4w6Q1DgAAkAls0Wg0NeN0d4BhRCQb9zl7OBxzdO7i7ycdRpxvL7eoKyQb7+nskO73ecbDiAAAALhzCSfIA0C6Cdp88o/74mrFuSUqjKb2E8wAMBWELQAZxz/umzA82PyXj8sf/v8BzH95RGFNbRFlwhsAkwhbAGaFwHhAxz9+O7ZdXFygh5ZsmtK5k4W37fc+rkI7YQvAzDFnCwAAwCDCFgAAgEGELQAAAIOYswUgrU02eX2qE98BIB0QtgCktckmrz+ycotF3QDA9BG2AKSNVDzFyrFLI+GLRq8BADcjbAFIG6l4ivX1JSJMXAMAbsYEeQAAAIN4sgXAEkx8B5AtCFsALMHEdwDZgmFEAAAAg3iyBQCTmOxTixJfUA1g+ghbADCJyT61KPEF1QCmj2FEAAAAgwhbAAAABhG2AAAADCJsAQAAGETYAgAAMIiwBQAAYBBLPwAwjq/mAZDNCFsAjJtNX80z2WKnBbn5ujp+La7G4qcAbiBsAcA0TLbY6SMrt0yosfgpgBuYswUAAGAQYQsAAMAgwhYAAIBBhC0AAACDphS2enp6VF9fr9raWnV0dEzY/8EHH8jlcmnDhg1qa2vTtWvXP5Xjdru1evVqNTQ0qKGhQe3t7cntHgAAIM0l/DSix+NRe3u7urq6lJ+fr82bN2vVqlVavHixJOmrr77SwYMH5Xa7dc8996i1tVVut1t/8zd/o7Nnz6qtrU0bN240/kIAAADSUcInW/39/aqqqlJpaamKiopUV1en3t7e2P6ioiL9+te/1j333KNAIKCRkRHNnTtXknTmzBm53W65XC499dRTunLlirlXAgAAkIYShq3h4WE5HI7YttPplMfjiTsmLy9PH374oe6//35dunRJq1evliQ5HA7t3r1bp06d0oIFC3Tw4MEktw8AAJDeEg4jRiIR2Wy22HY0Go3bvuG+++7Tb3/7Wx0+fFj79+/XSy+9pKNHj8b279ixQ2vXrp1Wc/Pnp3ZBQIdjTkqvB2twn1PPf3lExcUFcTW7PSfjatM5trAwT45S/q2lAu/p7JDJ9zlh2CovL9fg4GBs2+v1yul0xrYvX76ss2fPxp5muVwutba2amxsTCdPntRjjz0m6XpIs9vt02puZMSnSCQ6rXPulMMxR17vWEquBetwn60RDIfk91+Nq4XDEaO14uKCpF9jOscGgyH+raUA7+nskO73OSfHdtsHRAmHEaurqzUwMKDR0VEFAgH19fWppqYmtj8ajWrPnj26cOGCJKm3t1crV65UUVGRXnvtNX3yySeSpOPHj0/7yRYAZKob36F486+gzZf4RACzTsInW2VlZWptbVVTU5NCoZAaGxtVUVGh5uZmtbS0aPny5Xruuef0xBNPyGazafHixTpw4IDsdruOHDmi/fv3KxgMatGiRTp06FAqXhMAWG6y71Dk+xKB7DSlL6J2uVxyuVxxtWPHjsV+v2bNGq1Zs2bCeZWVlXK73TNsEQAAIHOxgjwAAIBBhC0AAACDCFsAAAAGEbYAAAAMmtIEeQCYqqDNJ/94/BIHYYUs6gYArEfYApBU/nGfXh96I672yMotFnUDANZjGBEAAMAgwhYAAIBBhC0AAACDCFsAAAAGEbYAAAAM4tOIAO4YyzwAQGKELQB3jGUeACAxwhYApEiOXRoJX4yrFeeWqDBaYlFHAFKBsAUAKRIYD+j4x2/H1bbf+7gK7YQtYDZjgjwAAIBBPNkCMCVMhgeAO0PYAjAlTIYHgDvDMCIAAIBBhC0AAACDCFsAAAAGEbYAAAAMImwBAAAYRNgCAAAwiLAFAABgEGELAADAIMIWAACAQYQtAAAAgwhbAAAABhG2AAAADJpS2Orp6VF9fb1qa2vV0dExYf8HH3wgl8ulDRs2qK2tTdeuXZMkXbhwQdu2bdO6deu0a9cu+f3+5HYPAACQ5hKGLY/Ho/b2dnV2dqq7u1snTpzQ+fPnY/u/+uorHTx4UP/0T/+kf/u3f9PVq1fldrslSQcOHNDWrVvV29urZcuW6eWXXzb3SgAAANJQwrDV39+vqqoqlZaWqqioSHV1dert7Y3tLyoq0q9//Wvdc889CgQCGhkZ0dy5cxUKhfTRRx+prq5OkrRp06a48wAAALJBwrA1PDwsh8MR23Y6nfJ4PHHH5OXl6cMPP9T999+vS5cuafXq1bp06ZJKSkqUm5srSXI4HBPOAwAAmO1yEx0QiURks9li29FoNG77hvvuu0+//e1vdfjwYe3fv19PP/30hOMmO+925s8vmdbxM+VwzEnp9WAN7vOd8V8eUXFxQVzNbs/JqpqJP7OwME+OUv5NzgTv6eyQyfc5YdgqLy/X4OBgbNvr9crpdMa2L1++rLNnz2r16tWSJJfLpdbWVs2bN09jY2MKh8Oy2+0TzpuKkRGfIpHotM65Uw7HHHm9Yym5FqzDfb5zwXBIfv/VuFo4HEnbWnFxQdKvYaLvYDDEv8kZ4D2dHdL9Pufk2G77gCjhMGJ1dbUGBgY0OjqqQCCgvr4+1dTUxPZHo1Ht2bNHFy5ckCT19vZq5cqVysvLU2VlpU6fPi1J6u7ujjsPACDl2KWR8MW4X0Gbz+q2ACRRwidbZWVlam1tVVNTk0KhkBobG1VRUaHm5ma1tLRo+fLleu655/TEE0/IZrNp8eLFOnDggCRp3759amtr0yuvvKIFCxbo8OHDxl8QAGSSwHhAxz9+O662/d7HVWhP7TQKAOYkDFvS9aFBl8sVVzt27Fjs92vWrNGaNWsmnLdw4UK99dZbM2wRAAAgc7GCPAAAgEGELQAAAIMIWwAAAAYRtgAAAAwibAEAABhE2AIAADBoSks/AJidgjaf/OMTF9Aszi1RYZR1ngAgGQhbQBbzj/v0+tAbE+osqgkAyUPYAjDBja+QuVlYIYu6AYDMRtgCMMFkXyHzyMotFnUDAJmNCfIAAAAGEbYAAAAMImwBAAAYRNgCAAAwiLAFAABgEJ9GBIA0M9nSGyw0C2QuwhYApJnJlt5goVkgczGMCAAAYBBPtgAgAzC0CGQuwhYAZACGFoHMxTAiAACAQYQtAAAAgwhbAAAABhG2AAAADGKCPJAlgjaf/OO+uFpYIYu6AYDsQdgCsoR/3KfXh96Iqz2ycotF3QBA9mAYEQAAwCDCFgAAgEEMIwKzEPOzsgOrygOZgbAFzELMz8oOrCoPZIYpDSP29PSovr5etbW16ujomLD/l7/8pRoaGvTAAw9o9+7dunLliiTJ7XZr9erVamhoUENDg9rb25PbPQAAQJpL+GTL4/Govb1dXV1dys/P1+bNm7Vq1SotXrxYkuTz+bR//36dPHlSZWVl+od/+Af99Kc/1bPPPquzZ8+qra1NGzduNP5CAAAA0lHCJ1v9/f2qqqpSaWmpioqKVFdXp97e3tj+UCikffv2qaysTJK0ZMkSffHFF5KkM2fOyO12y+Vy6amnnoo98QIAAMgWCcPW8PCwHA5HbNvpdMrj8cS27777bq1du1aSFAwG9eqrr2rNmjWSJIfDod27d+vUqVNasGCBDh48mOz+AQAA0lrCYcRIJCKbzRbbjkajcds3jI2N6Xvf+56WLl2qhx56SJJ09OjR2P4dO3bEQtlUzZ+f2kmeDseclF4P1siG++y/PKLi4oK4mt2eM6XadI7Ntlq69TNZrbAwT47S2f9v/GbZ8J5GZt/nhGGrvLxcg4ODsW2v1yun0xl3zPDwsLZv366qqirt3btX0vXwdfLkST322GOSroc0u90+reZGRnyKRKLTOudOORxz5PWOpeRasE623OdgOCS//2pcLRyOTKk2nWPTuVZcXJD0a6TT67tVLRgMZcW/8Ruy5T2d7dL9Pufk2G77gCjhMGJ1dbUGBgY0OjqqQCCgvr4+1dTUxPaHw2Ht3LlT69ev1zPPPBN76lVUVKTXXntNn3zyiSTp+PHj036yBQAAkOkSPtkqKytTa2urmpqaFAqF1NjYqIqKCjU3N6ulpUUXL17UZ599pnA4rPfff1+StGzZMr3wwgs6cuSI9u/fr2AwqEWLFunQoUPGXxAAAEA6mdKipi6XSy6XK6527NgxSdLy5ct17ty5Sc+rrKyU2+2eYYsAAACZi+9GBAAAMIiwBQAAYBBhCwAAwCDCFgAAgEGELQAAAIMIWwAAAAYRtgAAAAwibAEAABhE2AIAADCIsAUAAGAQYQsAAMAgwhYAAIBBhC0AAACDcq1uAMDMBG0++cd9cbWwQhZ1AwD4OsIWkOH84z69PvRGXO2RlVss6gYA8HUMIwIAABhE2AIAADCIYUQAmEVy7NJI+GJcrTi3RIXREos6AkDYAoBZJDAe0PGP346rbb/3cRXaCVuAVRhGBAAAMIiwBQAAYBBhCwAAwCDCFgAAgEGELQAAAIMIWwAAAAYRtgAAAAwibAEAABhE2AIAADCIsAUAAGDQlMJWT0+P6uvrVVtbq46Ojgn7f/nLX6qhoUEPPPCAdu/erStXrkiSLly4oG3btmndunXatWuX/H5/crsHskzQ5tNI+GLcr7BCVrcFALiNhGHL4/Govb1dnZ2d6u7u1okTJ3T+/PnYfp/Pp/379+vVV1/VqVOntGTJEv30pz+VJB04cEBbt25Vb2+vli1bppdfftncKwGygH/cp9eH3oj7NR4dt7otAMBtJAxb/f39qqqqUmlpqYqKilRXV6fe3t7Y/lAopH379qmsrEyStGTJEn3xxRcKhUL66KOPVFdXJ0natGlT3HkAgNTIsWvCE9GgzWd1W0DWyE10wPDwsBwOR2zb6XTq008/jW3ffffdWrt2rSQpGAzq1Vdf1d/+7d/q0qVLKikpUW7u9Us4HA55PJ5k9w8ASCAwHtDxj9+Oq22/93EV2kss6gjILgnDViQSkc1mi21Ho9G47RvGxsb0ve99T0uXLtVDDz0kj8cz4bjJzrud+fNT+4PA4ZiT0uvBGpl8n/2XR1RcXBBXs9tzkloz8WfOllq69TOTWmFhnhylmfteuFkmv6cxdZl8nxOGrfLycg0ODsa2vV6vnE5n3DHDw8Pavn27qqqqtHfvXknSvHnzNDY2pnA4LLvdPul5iYyM+BSJRKd1zp1yOObI6x1LybVgnUy/z8FwSH7/1bhaOBxJas3En2lFrbi4gL+b29SCwVBGvxduyPT3NKYm3e9zTo7ttg+IEs7Zqq6u1sDAgEZHRxUIBNTX16eamprY/nA4rJ07d2r9+vV65plnYk+v8vLyVFlZqdOnT0uSuru7484DAADIBgmfbJWVlam1tVVNTU0KhUJqbGxURUWFmpub1dLSoosXL+qzzz5TOBzW+++/L0latmyZXnjhBe3bt09tbW165ZVXtGDBAh0+fNj4CwIAAEgnCcOWJLlcLrlcrrjasWPHJEnLly/XuXPnJj1v4cKFeuutt2bYIgAAQOZiBXkAAACDCFsAAAAGEbYAAAAMmtKcLQDA7HJjVfmbFeeWqDDKQqdAshG2ACALsao8kDoMIwIAABhE2AIAADCIsAUAAGAQYQsAAMAgwhYAAIBBhC0AAACDCFsAAAAGsc4WkKaCNp/84764Wlghi7oBANwpwhaQpvzjPr0+9EZc7ZGVWyzqBgBwpxhGBAAAMIiwBQAAYBBhCwAAwCDCFgAAgEFMkAfSAJ88BIDZi7AFpAE+eQgAsxfDiAAAAAYRtgAAAAwibAEAABjEnC0AgCQpxy6NhC/G1YpzS1QYLbGoI2B2IGwBACRJgfGAjn/8dlxt+72Pq9BO2AJmgmFEAAAAgwhbAAAABhG2AAAADCJsAQAAGETYAgAAMGhKYaunp0f19fWqra1VR0fHLY97+umn1dXVFdt2u91avXq1Ghoa1NDQoPb29pl3DAAAkEESLv3g8XjU3t6urq4u5efna/PmzVq1apUWL14cd8y+ffs0MDCgqqqqWP3s2bNqa2vTxo0bzXQPAACQ5hI+2erv71dVVZVKS0tVVFSkuro69fb2xh3T09Oj73znO1q/fn1c/cyZM3K73XK5XHrqqad05cqV5HYPAACQ5hKGreHhYTkcjti20+mUx+OJO2bHjh16+OGHJ5zrcDi0e/dunTp1SgsWLNDBgweT0DIAAEDmSDiMGIlEZLPZYtvRaDRu+3aOHj0a+/2OHTu0du3aaTU3f35qVy12OOak9HqwRjreZ//lERUXF8TV7PYcS2pWXjvda+nWTypq+QV2+W0jcbW5BXM17655Shfp+J5G8mXyfU4YtsrLyzU4OBjb9nq9cjqdCf/gsbExnTx5Uo899pik6yHNbrdPq7mREZ8ikei0zrlTDscceb1jKbkWrJOu9zkYDsnvvxpXC4cjltSsvHYya8XFBfzdJKE2FvRN+hU+YV+e0kG6vqeRXOl+n3NybLd9QJRwGLG6uloDAwMaHR1VIBBQX1+fampqEl64qKhIr732mj755BNJ0vHjx6f9ZAsAACDTJXyyVVZWptbWVjU1NSkUCqmxsVEVFRVqbm5WS0uLli9fPul5drtdR44c0f79+xUMBrVo0SIdOnQo6S8AAAAgnSUMW5LkcrnkcrniaseOHZtw3Isvvhi3XVlZKbfbPYP2AAAAMhsryAMAABhE2AIAADCIsAUAAGDQlOZsAUieoM0n/7gvrhZWyKJuAACmEbaAFPOP+/T60BtxtUdWbrGoGwCAaYQtAMC05NilkfDFCfXi3BIVRlP7zR9AJiBsAQCmJTAemLCqvHR9ZflCO2EL+DomyAMAABhE2AIAADCIYUTAID55CAAgbAEG8clDAADDiAAAAAYRtgAAAAwibAEAABhE2AIAADCIsAUAAGAQYQsAAMAgwhYAAIBBhC0AAACDCFsAAAAGEbYAAAAMImwBAAAYxHcjAgCSIscujYQvxtWKc0tUGC2xqCMgPRC2AABJERgP6PjHb8fVtt/7uArthC1kN8IWAMAYnnYBhC0AgEE87QIIW0DSBG0++cd9cbWwQhZ1AwBIF4QtIEn84z69PvRGXO2RlVss6gYAkC5Y+gEAAMAgwhYAAIBBUwpbPT09qq+vV21trTo6Om553NNPP62urq7Y9oULF7Rt2zatW7dOu3btkt/vn3nHAAAAGSRh2PJ4PGpvb1dnZ6e6u7t14sQJnT9/fsIxO3fu1Pvvvx9XP3DggLZu3are3l4tW7ZML7/8cnK7BwAASHMJw1Z/f7+qqqpUWlqqoqIi1dXVqbe3N+6Ynp4efec739H69etjtVAopI8++kh1dXWSpE2bNk04DwAAYLZL+GnE4eFhORyO2LbT6dSnn34ad8yOHTskSUNDQ7HapUuXVFJSotzc65dwOBzyeDxJaRoAACBTJAxbkUhENpstth2NRuO2b2Wy46Zy3s3mz0/toncOx5yUXg/WMHWf/ZdHVFxcEFez23MyrpZu/aRTLd36SafadI4tLMyTozR570N+dmeHTL7PCcNWeXm5BgcHY9ter1dOpzPhHzxv3jyNjY0pHA7LbrdP+bybjYz4FIlEp3XOnXI45sjrHUvJtWAdk/c5GA7J778aVwuHIxlXS7d+7rRWXFzA300Ka9M5NhgMJe19yM/u7JDu9zknx3bbB0QJ52xVV1drYGBAo6OjCgQC6uvrU01NTcIL5+XlqbKyUqdPn5YkdXd3T+k8AACA2SRh2CorK1Nra6uampr04IMPauPGjaqoqFBzc7POnDlz23P37dund955R/X19RocHNQPfvCDpDUOWClo82kkfDHuF1/NAwCYzJS+rsflcsnlcsXVjh07NuG4F198MW574cKFeuutt2bQHpCe+GoeAMBUsYI8AACAQYQtAAAAgwhbAAAABhG2AAAADCJsAQAAGDSlTyMCAJAsOXZpJHwxrlacW6LCaGq/NQRIFcIWACClAuMBHf/47bja9nsfV6GdsIXZiWFEAAAAg3iyBSQQtPnkH/fF1VgtHgAwVYQtIAFWiwcAzARhCwBgOSbNYzYjbAEALMekecxmTJAHAAAwiLAFAABgEGELAADAIMIWAACAQYQtAAAAg/g0InATFjAFACQbYQu4CQuYAgCSjWFEAAAAg3iyBQBIS6wqj9mCsAUASEusKo/ZgrCFrMVkeABAKhC2kLWYDA8ASAUmyAMAABhE2AIAADCIsAUAAGAQYQsAAMAgwhYAAIBBfBoRAJAxvr7Qqf/yiOy2AhY6RVqbUtjq6enRK6+8ovHxcT366KPatm1b3P7f/e53euaZZ+T3+1VZWakDBw4oNzdXbrdbL730kubPny9Juv/++9Xa2pr8VwHcxs3rafkvjygYvr6WFmtqAZnn6wudFhcXaPPSbSx0irSWMGx5PB61t7erq6tL+fn52rx5s1atWqXFixfHjtmzZ4+ef/55rVixQnv37tU777yjrVu36uzZs2pra9PGjRuNvgjgdm5eT6u4uEB+/1VJrKkFAEiNhHO2+vv7VVVVpdLSUhUVFamurk69vb2x/Z9//rmCwaBWrFghSdq0aVNs/5kzZ+R2u+VyufTUU0/pypUrhl4GAABAekoYtoaHh+VwOGLbTqdTHo/nlvsdDkdsv8Ph0O7du3Xq1CktWLBABw8eTGbvAAAAaS/hMGIkEpHNZottR6PRuO3b7T969GisvmPHDq1du3Zazc2fn9oxeIdjTkqvh9TwXx5RcXFBbPvG7+32nLg6Neuvne61dOsnnWpWXruwME+OUn5+z3aZ/H90wrBVXl6uwcHB2LbX65XT6Yzb7/V6Y9tffvmlnE6nxsbGdPLkST322GOSrocwu90+reZGRnyKRKLTOudOORxz5PWOpeRaSK1gOBSbp3XznK1wOBL7/Q3U0q+fO60VFxfwd5PCmlXXLi4uUDAY4uf3LJfu/0fn5Nhu+4Ao4TBidXW1BgYGNDo6qkAgoL6+PtXU1MT2L1y4UAUFBRoaGpIkvfvuu6qpqVFRUZFee+01ffLJJ5Kk48ePT/vJFgAAidxYDuLmX0Gbz+q2gJiET7bKysrU2tqqpqYmhUIhNTY2qqKiQs3NzWppadHy5cv1k5/8RM8++6x8Pp+++c1vqqmpSXa7XUeOHNH+/fsVDAa1aNEiHTp0KBWvCQCQRb6+HIQkbb/3cZaDQNqY0jpbLpdLLpcrrnbs2LHY75cuXapf/OIXE86rrKyU2+2eYYsAAACZixXkMavcvIDpDSxeCgCwEmELs8rNC5jewOKlAAAr8UXUAAAABhG2AAAADGIYEQAw69xYDuJmxbklKozyCUWkHmELGYvJ8ABuheUgkE4IW8hYTIYHAGQC5mwBAAAYxJMtZASGDAHMFPO4YBXCFjICQ4YAZop5XLAKw4gAAAAG8WQLAJC1GFpEKhC2AABZi6FFpALDiAAAAAYRtgAAAAxiGBGWmmxJh4LcfF0dvxZXY5kHAECmImzBUrda0uHrcyhY5gFAqjBpHslG2AIA4CZMmkeyEbYAAEiAp12YCcIWUoav3AGQqXjahZkgbCFl+ModAEA2ImzBCJ5iAQBwHWELM3arYPXm0FtxNZ5iAZhNJpvHJTGXCxMRtjBjDA8CyEaTzeOSmMuFiQhbAAAk0WRPvCZbrJknYNmDsAUAQBJN9sRrssWaeQKWPQhbmBYmvgNAcrB2V/YgbOGWmPgOAOawdlf2IGzhlpj4DgDAzOVM5aCenh7V19ertrZWHR0dE/b/7ne/06ZNm1RXV6dnnnlG4+PjkqQLFy5o27ZtWrdunXbt2iW/35/c7pE0QZtPI+GLcb8YHgSA1LoxtHjzr6DNl/hEpLWET7Y8Hm1LR+gAAAYXSURBVI/a29vV1dWl/Px8bd68WatWrdLixYtjx+zZs0fPP/+8VqxYob179+qdd97R1q1bdeDAAW3dulUbNmzQ0aNH9fLLL2vPnj1GX1C2mmzI71Zj/wwPAkB6mmxosfkvH5c/PLWf70hPCcNWf3+/qqqqVFpaKkmqq6tTb2+vnnzySUnS559/rmAwqBUrVkiSNm3apH/8x3/Uww8/rI8++khHjx6N1R955BHCliGTDflN9gaVCFYAkEmmGsAmW16CJSfSQ8KwNTw8LIfDEdt2Op369NNPb7nf4XDI4/Ho0qVLKikpUW5ublx9OnJybNM6fqZSfb2vu2r7SoFw/FBrvj1P18KhhLUcu3T3XaVxtfHouLo/OzXhOg8ue2DCsXn2vKyoFRUWKD9yNW36ScdauvVzp7WiwgL+blJYs+raJu5zJvw9TPbz/cFlD0yptvlb39W1aDCuNtn/K3fZi1UQLZrw+qxi9f/Rt5OoN1s0Go3e7oBXXnlFV69e1Q9+8ANJ0jvvvKOzZ8/q4MGDkqShoSG99NJL6uzslCT94Q9/0M6dO/XP//zP+u53v6sPP/xQkjQ+Pq6/+Iu/0JkzZ2b8ogAAADJFwgny5eXl8nq9sW2v1yun03nL/V9++aWcTqfmzZunsbExhcPhSc8DAADIBgnDVnV1tQYGBjQ6OqpAIKC+vj7V1NTE9i9cuFAFBQUaGhqSJL377ruqqalRXl6eKisrdfr0aUlSd3d33HkAAADZIOEwonR96Yef//znCoVCamxsVHNzs5qbm9XS0qLly5fr3LlzevbZZ+Xz+fTNb35TP/rRj5Sfn6/PP/9cbW1tGhkZ0YIFC3T48GH9yZ/8SSpeFwAAQFqYUtgCAADAnZnSoqYAAAC4M4QtAAAAgwhbAAAABhG2AAAADCJsAQAAGETY+prPPvtMy5Yts7oNGDQ0NKTGxkY1NDTo0Ucf1eeff251S0iinp4e1dfXq7a2Vh0dHVa3A0N+9rOfacOGDdqwYYMOHTpkdTsw7Mc//rHa2tqsbuOOEbZuEggE9NxzzykUCiU+GBlrz549ev755/Xuu+/K5XLp+eeft7olJInH41F7e7s6OzvV3d2tEydO6Pz581a3hSTr7+/Xf/zHf8jtdqu7u1v//d//rQ8++MDqtmDIwMCA3G631W3MCGHrJi+++KIeffRRq9uAQdeuXdP3v/99LV26VJK0ZMkSffHFFxZ3hWTp7+9XVVWVSktLVVRUpLq6OvX29lrdFpLM4XCora1N+fn5ysvL05/92Z/pwoULVrcFAy5fvqz29nbt3LnT6lZmhLD1//zqV79SMBjUunXrrG4FBuXn56uhoUGSFIlE9LOf/Uxr1qyxuCsky/DwsBwOR2zb6XTK4/FY2BFM+MY3vqEVK1ZIkv7whz/ovffe03333WdxVzDhhz/8oVpbWzV37lyrW5mRXKsbSLX33ntPP/rRj+Jqf/qnfyqfz6c333zTmqZgxK3u9Ztvvqlr166pra1N4+PjeuKJJyzqEMkWiURks9li29FoNG4bs8vvf/97PfHEE3r66ae1aNEiq9tBkv3rv/6rFixYoG9/+9vq6uqyup0Z4et6dP2G/vznP1dxcbEk6dy5c1q6dKk6OjpUUlJicXdINr/fr127dqm0tFQ/+clPlJ+fb3VLSBK3263BwUG98MILkqSjR48qGo3qySeftLgzJNvQ0JBaWlq0d+9ebdiwwep2YMDf/d3fyev1ym6368qVK/rqq6/04IMPau/evVa3Nm2ErUksWbJE//M//2N1GzBk9+7dmj9/vg4cOKCcHEbSZxOPx6MtW7boF7/4he666y5t3rxZzz33nCoqKqxuDUn0xRdf6KGHHlJ7e7u+/e1vW90OUqCrq0v/+Z//qRdffNHqVu5I1g0jIrt99tln+tWvfqXFixfroYceknR9Xs+xY8cs7gzJUFZWptbWVjU1NSkUCqmxsZGgNQu9/vrrunr1atx/vJs3b9aWLVss7Aq4NZ5sAQAAGMQYCgAAgEGELQAAAIMIWwAAAAYRtgAAAAwibAEAABhE2AIAADCIsAUAAGAQYQsAAMCg/wPwrArQTrub7gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 100000 \n",
    "p = [0.0 for k in range(N)]\n",
    "p = drawp(N)\n",
    "num_bins = 100\n",
    "plt.figure(figsize = [10,6])\n",
    "plt.hist(p,num_bins, density= 1.0, facecolor='green', alpha = 0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Leap Frog "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use leap frog approximation to evolve the system according to time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def leap_frog(N,dt,ix,ip,np):\n",
    "    \n",
    "        ''' N : number of steps to evolve\n",
    "            dt:  fraction of time ie T  = dt*N\n",
    "            ix,ip : initial position and momentum\n",
    "            np : number of the particles in the system\n",
    "        '''\n",
    "        ''' Returns\n",
    "          x,p : final position and momentum'''   \n",
    "    \n",
    "    \n",
    "        x = ix\n",
    "        p = ip\n",
    "        k = 0\n",
    "        while k < N:\n",
    "            if k == 0:\n",
    "                    for i in range(np):\n",
    "                        p[i] = p[i] - ((dt/2.0)*x[i])\n",
    "            elif k > 0 :\n",
    "                    if k < N - 1:\n",
    "                        for i in range(np):   \n",
    "                            x[i] = x[i] + (dt*p[i])\n",
    "                            p[i] =   p[i] - (dt*x[i])\n",
    "                            #S1 = hamiltonian(x,p,np)\n",
    "                            #print \"k =\",k,\"S1=\",S1\n",
    "                            \n",
    "            elif k == N - 1:\n",
    "                          for i in range(np): \n",
    "                            p[i] = (p[i] - (dt/2.0)*x[i])\n",
    "                       \n",
    "            k = k+1\n",
    "        return x,p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HMC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we run the HMC - simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def HMC(np,N,dt,steps,x0):\n",
    "    \n",
    "    ''' np : number of particles in the system\n",
    "        N = number of steps in Leap - Frog\n",
    "        dt = fraction of time in Leap - Frog\n",
    "        steps: total steps in HMC '''\n",
    "    \n",
    "    \n",
    "       \n",
    "    xt = [0.0 for k in range(np)]\n",
    "    pt = [0.0 for k in range(np)]\n",
    "\n",
    "    \n",
    "    p0 = drawp(np)\n",
    "    H = [0.0 for k in range(steps)]\n",
    "       \n",
    "    S0 = hamiltonian(x0,p0,np)\n",
    "    #print (\"=======>\", 0,\"S0=\", S0)\n",
    "\n",
    "\n",
    "\n",
    "    chain = 1\n",
    "    total_frac =  0.0\n",
    "    while chain < steps:\n",
    "        s_stor = [0.0]\n",
    "        xt,pt = leap_frog(N,dt,x0,p0,np)\n",
    "        S1 = hamiltonian(xt,pt,np)\n",
    "        frac = math.exp(-(S1-S0))\n",
    "        #print frac\n",
    "        a = min(1,frac)\n",
    "        b = random.uniform(0.0,1.0)\n",
    "\n",
    "        if b < a:\n",
    "            #print(\"=======>\", chain, \"S1=\",S1,frac,a,b)\n",
    "            H[chain] = S1\n",
    "            x0 = xt\n",
    "            p0 = drawp(np)\n",
    "            S0 = hamiltonian(x0,p0,np)\n",
    "        else:\n",
    "            H[chain] = S0\n",
    "            p0 = drawp(np)\n",
    "            \n",
    "        chain = chain+1\n",
    "            \n",
    "    return H   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seting Constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "np = 1000\n",
    "N = 1000\n",
    "dt = 0.001\n",
    "steps = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Call HMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = [1.0 for k in range(np)]\n",
    "x0 = [random.uniform(0.0,1.0) for k in range(np)]\n",
    "H = HMC(np,N,dt,steps,x0)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3UAAAD7CAYAAADNXL0VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df2xUdb7/8deUli5uIdo6YwmyZEXy5futSpPFH1XSolFaqF3vzrpZRFM3xq/iNeiS3bosEAiJLuBtbLJfrZvcEK9L/GPrilSbMmiCNkJJVprdJSLcGAPIIpahoFBoS9s53z+8nYjQOdM5c358Zp6PxMj8aOdzzrzP+ZzX+XzOaciyLEsAAAAAACMV+N0AAAAAAEDmCHUAAAAAYDBCHQAAAAAYjFAHAAAAAAYj1AEAAACAwQh1AAAAAGAwQh0AAAAAGKzQ7wak68yZ80okgvUn9crKStTX1+93M5DjqDN4gTqD26gxeIE6gxf8qLOCgpCuueaH475uTKhLJKzAhTpJgWwTcg91Bi9QZ3AbNQYvUGfwQtDqjOmXAAAAAGAwQh0AAAAAGIxQBwAAAAAGI9QBAAAAgMGMuVEKAAAAkE/2HvhK27o+V9/ZIZVNK1a0ZraqKsr9bhYCiFAHAAAABMzeA1/p9R2HdHEkIUnqOzuk13cckiSCHS5DqIMnONMEAACQvm1dnycD3ZiLIwlt6/qcYyhchlAH13GmCcD3BflEj13bgtx2wGRsW5fqOzs0oedNwnedfYQ6uI4zTQC+K8gneuzaFuS2AybL120rVbgpm1Z8xQBXNq3Y62ZmVb5+124j1MF16ZxpcvOMDWeDkMtM3HaCfKLHrm1BbjtgsnzctuzCTbRm9iWvS9LkwgJFa2b70t5sycfv2guEOrjO7kyTm2ds0vndhD6YytRtJ8hTiuzaFuS2A9ngV5+Yj9uWXbgZW++5doySj9+1Fwh1cJ3dmSY3z9jY/W6mAGQfIXniMr2Gy9RtJ50pRX7VkV3bcnU6FCD5Oy0unRPAuda3pBNuvhvuTJKP00r9RqiD6+zONLl5xsbudzMFILucju7kYqdtx8k1XNnYdsbW+emzQyr9zjp3c9uxO9Hj5wi7XdtydToUIKW3Xfux7eXqCViTw41dX56P00r9RqiDJ1KdaXJzp2b3u50eFOdjCEnFyeiOpJzstO04uYbL6baT6vtwc9uxO9Hj5wi7XdtydToUINlv1063vVR9Zqptq6l1j9EnYMdb7myEGz+OQ+zqIMjTSnP5uI1QB9+5ecbG7nc7OSjO1TOHTjgZ3Rn795Vey+UdvZNruP5vw/9xtO2k+j7c3Hak1Cd6/B5ht5vuZOp0KMCO3XbtZNtLp88cb9uayA3Xvj/rwG/pLHemfY/bxyGZTv13Oq3Urf4414/b0gp1/f39Wrp0qf70pz/p+uuvV3d3tzZu3KihoSEtXrxYK1eulCQdPHhQa9as0fnz5zV//nxt2LBBhYWF+vLLL9XU1KS+vj79+Mc/VnNzs374wx+6umAwh5tnbOx+t5NAydTNy7kxuuPFXVKDfB1Jqtedbjupvg83tx07bo+wA0EQxBEDu+3aybbnpM/084ZrTqUzapVpG908DnEy9T+oJ8xz/bjNNtT985//1Nq1a3XkyBFJ0uDgoFavXq2tW7dq+vTpevLJJ9XV1aWamho1NTXp+eefV2VlpVavXq22tjYtW7ZMGzZs0LJly1RfX69XXnlFra2tampqcnvZYBA3z3yn+t12B66pOl2/DyxNPCDI5GDdi07bzx2902u4nGw7dt+Hk23HCbdHCZE7grgfHOPkmiO/2G3XTm5m4qTPdPuGa27WkZ/3Dcj0JlxS6nVqVwd+nzAfb7n8Pm5zm22oa2tr0/r16/Xcc89Jkvbv369Zs2Zp5syZkqSGhgbFYjHdeOONGhwcVGVlpSQpGo3qj3/8o37xi1/o448/1iuvvJJ8/pFHHiHUuSTIHVym3F6m8Q5c7TpdPw8sTT0gsNvR+3WXVLd39JleR5LO6044HW1z62SM0zpyciADcwT5T9Y4vebIT6m2ayc3M3HSZ9rtE5zsw93uT/26b4CTm3DZBSC7qf9O+i2n/bGTa8VNZxvqXnjhhUsenzx5UuFwOPk4Eomot7f3sufD4bB6e3t15swZlZSUqLCw8JLnJ6qsrGTCP+OFcHiq301I+rDnmP4c+28NDY9K+raQ/xz7b02b+gMt/MlMn1uXGT+XafvuvVfsdLfvPqyfLpyjX91foZff/GeybZJUXDRJv7q/Iut18f3fZ9c2P/104dRx2/DThVM1beoP9OcdB3XqzICuvWaKGhf/7+R3meq10+Ps0E+fHUp7fX/Yc+yKvz98zRTFzwxc9v7wNVMcf5fp1HCqdZbO65my+z6cGm99p9u2TOrIbn3n4n4y6JzUQSp2+8Eg9x/Z2J/5IdW299jz77naZ6baJ6SzDx+vDt3uT908Vkj1u/+842DK5bJb7lTr9KcL59j2HZn2W07741TLle3vImjb6oRvlJJIJBQKhZKPLctSKBQa9/mx/3/X9x+no6+vX4mENeGfc1M4PFXx+Dm/m5H0Xx0HLilUSRoaHtV/dRxQxY+u9qlVzvi5TFfaqYw9H4+fU8WPrlZj3f+67ExUxY+uzmpdXKnO7NrmlJtntyt+dLU2P1l1yXNjbU71Wuk4Z9hKpxWntczfP3sXPzOg/9f2D509N6h/W/DjK551/LcFP3a8Pv3eLu2+y7F1/t06y1YNjbe+s1FL49WK3fr2+/uw4+fIkhfXq2azDuz2g0HuP5zuz/w03rbnZ59ptw9PVYdu96duLnfK322zXHav263TVP21E07741TLlc3vwo8MUFAQSjnINeFQV15erng8nnwcj8cViUQue/7UqVOKRCIqLS3VuXPnNDo6qkmTJiXfj+zLxbnCfi5TOsP0dtPP3DpQcnMKQVCndjqdKphqutN//Ptdyfdk+7vys4b9/C79ml5mt76DvJ/06/sy9XrVIN9Qx+k1R26Gez/7JS+mbF/p7pdu/mmYdNvn9X0DnNyEa+z3St7/2QGnn+vkWnHTTTjUzZs3T4cPH9bRo0d1/fXXq6OjQz//+c81Y8YMFRcXq6enRz/5yU/U3t6u6upqFRUVaf78+ers7FRDQ4O2b9+u6upqN5Yl7+XiXOF0lsmvP0Jsx80DJTevKUrnb4T5MZKQzo7eyUX6bu3o/dwu/bxux68DaqcHMk452T7c/r4yvT25k+Vysw6CfEOddG5wJF15f+bmtYJ+9ktuG9uHX2kExcn1YZKZ1+E6vQmX5F8AcvK52Th2M+27HjPhUFdcXKxNmzZpxYoVGhoaUk1Njerq6iRJzc3NWrt2rfr7+1VRUaHGxkZJ0vr167Vq1Sq9+uqrmj59ul566aXsLkUeSVVsuVjI6YQXv/4IsR03D9CcHBA4uTja71G8VDt6v29s4+Yfl81UkEcq3JKNA5lMOd0+3Py+nNye3MlyuVkHTm+okw67O1g6ucHRePuzdE6sZfp9+NUv+S1VHdq12+9+L1N+3oTLT06Wy9TvekzaoW7Xrl3Jf1dVVemdd9657D1z587VX//618uenzFjhrZu3ZphEzHGrtj8LGS3AqHdMrl9ZtvJ2SK3D6gzPSCwez1V5xfkO7bZtc3vg3k/Os4gj1S4xc8DGafbh5Pbxdu97mT6mZORvHTqwEn/kWof7XR0P9V2Lcl2m8+0/7DrO5x8H371S35LZ+R0vHa7OZLtNrvvI6jflx279Z3pcgX5GCcdEx6pg3/SKTa7kQw3wpHbZzZSLVOQr4/x64Da6TVFqTq//3z30wl9ppfSmV4p+XMw71fH6ecooZ9h1q8DGaf7Iye3i3fz9uTZGMkb71ont/8sgZPR/VTb9di/r/Sa09pyeq1gPt/SfTxO9kdujmRj4txc30E+pkwHoc4gTorNSYdvx88zG0HuoLI1HfZKB0KpuHlx9Nhz4/2sn/y8SD+oHYHf02tMPQucKafXAKf6vppa97g2Am9XJ05H8sb+u9K1Tm5ONbRj99mZbNfZ2OadXiuYarn8vu7NT5nuj7Ixko3s8fPmS0FHqDOIk2Jz0uGP8Ws6RypB7qDSmafvxhx+Ny+ODvL69rNtQe4I8i1Y+Skb1wCP9325OQKf6nPT+VknfUA2phpmyu6zMxkxC8K1gqmWy+8TPSZys/4xcW6u7yAf46SDUGcQJ8XmtMMP6nSOoHdQ4x0oOZn2Y7dsbl5TFOT17WfbTO8IkB12Nehku3ZzBN7pcjnpA5xONXTC7rPttms3t3kn1wqmUwtB2Gebws36x8T5efOloCPUBYzTu2mNx2mHH+TpHCZ2UG5M+/kuN68pCvL69qttpncEyB63rgF2cwQ+HU5G8lJxOtXQiXRGL6XU23UQpzb73SfnItZ3cLi9voN8jGOHUBcgTqbm2HHa4TOdI7ucTvtB8JjcEcAbTrZrN0fgnXJzFNDNA7h02m03YhbEbZ4+2Vusb2+xvscXsizL8rsR6ejr61ciEaymXumibyeaWveM2+H/x7/f5fj3O7mDmNttyzd26/P7AV/69kDm0cVz2XH5JKi3rM6WbO/PcLl8364zrbFc3/aQXezL4AU/6qygIKSyspJxX2ekLkDcvtjWyVlFphdk10Sm/Uz07pfIPm5ZjWzgDHNmgjoiBgBBQqgLkCBPueNgJLsmMu2Hs47+45bVyBYCCgDADYS6AAn6aBgHI9nF+jQHt6wGAABBRqjzwXjXBzAaBgRTkEfRAQAACHUes7s2h9EbIHiCPooOAADyW4HfDcg3qa7NARBMVRXlenTx3Ev+0HO+3LEQAAAEHyN1HuPaHMBMjKIDAICgYqTOY+Ndg8O1OQAAAAAyQajzWLRmtiYXXrrauTYHAAAAQKaYfukx7nAJAAAAIJsIdT7g2hwAAAAA2cL0SwAAAAAwGKEOAAAAAAxGqAMAAAAAgxHqAAAAAMBghDoAAAAAMBihDgAAAAAMRqgDAAAAAIMR6gAAAADAYIQ6AAAAADAYoQ4AAAAADOYo1LW3t6u+vl719fXavHmzJOngwYOKRqOqra3VmjVrNDIyIkn68ssv9fDDD6uurk5PPfWUzp8/77z1AAAAAJDnMg51AwMDeuGFF7R161a1t7dr37596u7uVlNTk9atW6edO3fKsiy1tbVJkjZs2KBly5YpFovppptuUmtra9YWAgAAAADyVcahbnR0VIlEQgMDAxoZGdHIyIgKCws1ODioyspKSVI0GlUsFtPw8LA+/vhj1dbWXvI8AAAAAMCZwkx/sKSkRM8++6wWL16sKVOm6NZbb1VRUZHC4XDyPeFwWL29vTpz5oxKSkpUWFh4yfMAAAAAAGcyDnWHDh3SW2+9pQ8++EBTp07Vb3/7W+3Zs0ehUCj5HsuyFAqFkv//ru8/tlNWVpJpU10VDk/1uwnIA9QZvECdwW3UGLxAncELQauzjEPd7t27VVVVpbKyMknfTqncsmWL4vF48j2nTp1SJBJRaWmpzp07p9HRUU2aNEnxeFyRSGRCn9fX169Ewsq0ua4Ih6cqHj/ndzOQ46gzeIE6g9uoMXiBOoMX/KizgoJQykGujK+pmzt3rrq7u3XhwgVZlqVdu3bptttuU3FxsXp6eiR9e3fM6upqFRUVaf78+ers7JQkbd++XdXV1Zl+NAAAAADgf2Q8UrdgwQJ9+umnikajKioq0s0336wnnnhC9913n9auXav+/n5VVFSosbFRkrR+/XqtWrVKr776qqZPn66XXnopawsBAAAAAPkqZFlWsOY0joPpl8hX1Bm8QJ3BbdQYvECdwQs5Nf0SAAAAAOA/Qh0AAAAAGIxQBwAAAAAGI9QBAAAAgMEIdQAAAABgMEIdAAAAABiMUAcAAAAABiPUAQAAAIDBCHUAAAAAYDBCHQAAAAAYjFAHAAAAAAYj1AEAAACAwQr9bkAu2nvgK23r+lx9Z4dUNq1Y0ZrZqqoo97tZAAAAAHIQoS7L9h74Sq/vOKSLIwlJUt/ZIb2+45AkEewAAAAAZB3TL7NsW9fnyUA35uJIQtu6PvepRQAAAAByGaEuy/rODk3oeQAAAABwglCXZWXTiif0PAAAAAA4QajLsmjNbE0uvHS1Ti4sULRmtk8tAgAAAJDLuFFKlo3dDIW7XwIAAADwAqHOBVUV5YQ4AAAAAJ5g+iUAAAAAGIxQBwAAAAAGI9QBAAAAgMEIdQAAAABgMEIdAAAAABiMUAcAAAAABiPUAQAAAIDBCHUAAAAAYDBHoW7Xrl2KRqNavHixnn/+eUlSd3e3GhoatGjRIrW0tCTfe/DgQUWjUdXW1mrNmjUaGRlx1nIAAAAAQOah7tixY1q/fr1aW1v1zjvv6NNPP1VXV5dWr16t1tZWdXZ26pNPPlFXV5ckqampSevWrdPOnTtlWZba2tqythAAAAAAkK8yDnXvv/++lixZovLychUVFamlpUVTpkzRrFmzNHPmTBUWFqqhoUGxWEzHjx/X4OCgKisrJUnRaFSxWCxrCwEAAAAA+aow0x88evSoioqKtHz5cp04cUILFy7UnDlzFA6Hk++JRCLq7e3VyZMnL3k+HA6rt7fXWcsBAAAAAJmHutHRUe3bt09bt27VVVddpaeeeko/+MEPFAqFku+xLEuhUEiJROKKz09EWVlJpk11VTg81e8mIA9QZ/ACdQa3UWPwAnUGLwStzjIOdddee62qqqpUWloqSbr33nsVi8U0adKk5Hvi8bgikYjKy8sVj8eTz586dUqRSGRCn9fX169Ewsq0ua4Ih6cqHj/ndzOQ46gzeIE6g9uoMXiBOoMX/KizgoJQykGujK+pu/vuu7V7926dPXtWo6Oj+uijj1RXV6fDhw/r6NGjGh0dVUdHh6qrqzVjxgwVFxerp6dHktTe3q7q6upMPxoAAAAA8D8yHqmbN2+eHn/8cS1btkzDw8O666679NBDD+mGG27QihUrNDQ0pJqaGtXV1UmSmpubtXbtWvX396uiokKNjY1ZWwgAAAAAyFchy7KCNadxHEy/RL6izuAF6gxuo8bgBeoMXsip6ZcAAAAAAP8R6gAAAADAYIQ6AAAAADAYoQ4AAAAADEaoAwAAAACDEeoAAAAAwGCEOgAAAAAwGKEOAAAAAAxGqAMAAAAAgxHqAAAAAMBghDoAAAAAMBihDgAAAAAMRqgDAAAAAIMR6gAAAADAYIQ6AAAAADAYoQ4AAAAADEaoAwAAAACDEeoAAAAAwGCEOgAAAAAwGKEOAAAAAAxGqAMAAAAAgxHqAAAAAMBghDoAAAAAMBihDgAAAAAMRqgDAAAAAIMR6gAAAADAYIQ6AAAAADAYoQ4AAAAADOY41G3evFmrVq2SJB08eFDRaFS1tbVas2aNRkZGJElffvmlHn74YdXV1empp57S+fPnnX4sAAAAAEAOQ93evXv19ttvJx83NTVp3bp12rlzpyzLUltbmyRpw4YNWrZsmWKxmG666Sa1trY6azUAAAAAQJKDUPf111+rpaVFy5cvlyQdP35cg4ODqqyslCRFo1HFYjENDw/r448/Vm1t7SXPAwAAAACcyzjUrVu3TitXrtS0adMkSSdPnlQ4HE6+Hg6H1dvbqzNnzqikpESFhYWXPA8AAAAAcK4wkx968803NX36dFVVVWnbtm2SpEQioVAolHyPZVkKhULJ/3/X9x+no6ysJJOmui4cnup3E5AHqDN4gTqD26gxeIE6gxeCVmcZhbrOzk7F43E98MAD+uabb3ThwgWFQiHF4/Hke06dOqVIJKLS0lKdO3dOo6OjmjRpkuLxuCKRyIQ/s6+vX4mElUlzXRMOT1U8fs7vZiDHUWfwAnUGt1Fj8AJ1Bi/4UWcFBaGUg1wZTb987bXX1NHRofb2dj3zzDO65557tHHjRhUXF6unp0eS1N7erurqahUVFWn+/Pnq7OyUJG3fvl3V1dWZfCwAAAAA4Huy+nfqmpubtXHjRtXV1enChQtqbGyUJK1fv15tbW1asmSJ9u3bp1//+tfZ/FgAAAAAyFshy7KCNadxHEy/RL6izuAF6gxuo8bgBeoMXsiZ6ZcAAAAAgGAg1AEAAACAwQh1AAAAAGAwQh0AAAAAGIxQBwAAAAAGI9QBAAAAgMEIdQAAAABgMEIdAAAAABiMUAcAAAAABiPUAQAAAIDBCHUAAAAAYDBCHQAAAAAYjFAHAAAAAAYj1AEAAACAwQh1AAAAAGAwQh0AAAAAGIxQBwAAAAAGI9QBAAAAgMEIdQAAAABgMEIdAAAAABiMUAcAAAAABiPUAQAAAIDBCHUAAAAAYDBCHQAAAAAYjFAHAAAAAAYj1AEAAACAwQh1AAAAAGAwQh0AAAAAGIxQBwAAAAAGcxTqXn75ZdXX16u+vl4vvviiJKm7u1sNDQ1atGiRWlpaku89ePCgotGoamtrtWbNGo2MjDhrOQAAAAAg81DX3d2t3bt36+2339b27dt14MABdXR0aPXq1WptbVVnZ6c++eQTdXV1SZKampq0bt067dy5U5Zlqa2tLWsLAQAAAAD5KuNQFw6HtWrVKk2ePFlFRUWaPXu2jhw5olmzZmnmzJkqLCxUQ0ODYrGYjh8/rsHBQVVWVkqSotGoYrFY1hYCAAAAAPJVYaY/OGfOnOS/jxw5oh07duiRRx5ROBxOPh+JRNTb26uTJ09e8nw4HFZvb++EPq+srCTTproqHJ7qdxOQB6gzeIE6g9uoMXiBOoMXglZnGYe6MZ999pmefPJJPffcc5o0aZKOHDmSfM2yLIVCISUSCYVCocuen4i+vn4lEpbT5mZVODxV8fg5v5uBHEedwQvUGdxGjcEL1Bm84EedFRSEUg5yObpRSk9Pj371q1/pN7/5jX72s5+pvLxc8Xg8+Xo8HlckErns+VOnTikSiTj5aAAAAACAHIS6EydO6Omnn1Zzc7Pq6+slSfPmzdPhw4d19OhRjY6OqqOjQ9XV1ZoxY4aKi4vV09MjSWpvb1d1dXV2lgAAAAAA8ljG0y+3bNmioaEhbdq0Kfnc0qVLtWnTJq1YsUJDQ0OqqalRXV2dJKm5uVlr165Vf3+/Kioq1NjY6Lz1AAAAAJDnQpZlBetCtXFwTR3yFXUGL1BncBs1Bi9QZ/BCzl1TBwAAAADwF6EOAAAAAAxGqAMAAAAAgxHqAAAAAMBghDoAAAAAMBihDgAAAAAMRqgDAAAAAIMR6gAAAADAYIQ6AAAAADAYoQ4AAAAADEaoAwAAAACDEeoAAAAAwGCEOgAAAAAwGKEOAAAAAAxGqAMAAAAAgxHqAAAAAMBghDoAAAAAMBihDgAAAAAMRqgDAAAAAIMR6gAAAADAYIQ6AAAAADAYoQ4AAAAADEaoAwAAAACDEeoAAAAAwGCEOgAAAAAwGKEOAAAAAAxGqAMAAAAAgxHqAAAAAMBghV5+2LvvvqtXX31VIyMjevTRR/Xwww97+fFZs/fAV9rW9blOnx1S6bRiRWtmq6qi3O9mAQAAAMhDnoW63t5etbS0aNu2bZo8ebKWLl2q22+/XTfeeKNXTciKvQe+0us7DuniSEKS1Hd2SK/vOCRJBDsAAAAAnvNs+mV3d7fuuOMOXX311brqqqtUW1urWCzm1cdnzbauz5OBbszFkYS2dX3uU4sAAAAA5DPPRupOnjypcDicfByJRLR///60f76srMSNZk3Y6bND4z4fDk/1uDXIF9QWvECdwW3UGLxAncELQaszz0JdIpFQKBRKPrYs65LHdvr6+pVIWG40bUJKpxWr7wrBrnRaseLxcz60CLkuHJ5KbcF11BncRo3BC9QZvOBHnRUUhFIOcnk2/bK8vFzxeDz5OB6PKxKJePXxWROtma3JhZeutsmFBYrWzPapRQAAAADymWeh7s4779TevXt1+vRpDQwM6L333lN1dbVXH581VRXlenTxXJVNK1ZIUtm0Yj26eC43SQEAAADgC8+mX1533XVauXKlGhsbNTw8rAcffFC33HKLVx+fVVUV5aqqKGeIHwAAAIDvPP07dQ0NDWpoaPDyIwEAAAAgp3k2/RIAAAAAkH2EOgAAAAAwGKEOAAAAAAzm6TV1ThQUpP837bwU1HYht1Bn8AJ1BrdRY/ACdQYveF1ndp8XsizL/7/oDQAAAADICNMvAQAAAMBghDoAAAAAMBihDgAAAAAMRqgDAAAAAIMR6gAAAADAYIQ6AAAAADAYoQ4AAAAADEaoAwAAAACDEeoAAAAAwGCEugy9++67WrJkiRYtWqQ33njD7+YgR7z88suqr69XfX29XnzxRUlSd3e3GhoatGjRIrW0tPjcQuSSzZs3a9WqVZKkgwcPKhqNqra2VmvWrNHIyIjPrYPpdu3apWg0qsWLF+v555+XxP4M2dXe3p7sMzdv3iyJfRmyp7+/X/fff7/+9a9/SRp//xWYmrMwYV999ZV19913W2fOnLHOnz9vNTQ0WJ999pnfzYLh9uzZY/3yl7+0hoaGrIsXL1qNjY3Wu+++a9XU1FhffPGFNTw8bD322GPWhx9+6HdTkQO6u7ut22+/3frd735nWZZl1dfXW3//+98ty7Ks3//+99Ybb7zhZ/NguC+++MJasGCBdeLECevixYvWQw89ZH344Yfsz5A1Fy5csG699Varr6/PGh4eth588EFrz5497MuQFf/4xz+s+++/36qoqLCOHTtmDQwMjLv/CkrNMVKXge7ubt1xxx26+uqrddVVV6m2tlaxWMzvZsFw4XBYq1at0uTJk1VUVKTZs2fryJEjmjVrlmbOnKnCwkI1NDRQa3Ds66+/VktLi5YvXy5JOn78uAYHB1VZWSlJikaj1Bkcef/997VkyRKVl5erqKhILS0tmjJlCvszZM3o6KgSiYQGBgY0MjKikZERFRYWsi9DVrS1tWn9+vWKRCKSpP37919x/xWk/rPQl0813MmTJxUOh5OPI5GI9u/f72OLkAvmzJmT/PeRI0e0Y8cOPfLII5fVWm9vrx/NQw5Zt26dVq5cqRMnTki6fJ8WDoepMzhy9OhRFRUVafny5Tpx4oQWLlyoOXPmsD9D1pSUlOjZZ5/V4sWLNWXKFN16660qKipiX4aseOGFFy55fKVj/97e3rqe7o4AAAKYSURBVED1n4zUZSCRSCgUCiUfW5Z1yWPAic8++0yPPfaYnnvuOc2cOZNaQ1a9+eabmj59uqqqqpLPsU9Dto2Ojmrv3r36wx/+oL/85S/av3+/jh07Rp0haw4dOqS33npLH3zwgT766CMVFBRoz5491BhcMV4/GaT+k5G6DJSXl2vfvn3Jx/F4PDk8CzjR09OjZ555RqtXr1Z9fb3+9re/KR6PJ1+n1uBUZ2en4vG4HnjgAX3zzTe6cOGCQqHQJXV26tQp6gyOXHvttaqqqlJpaakk6d5771UsFtOkSZOS72F/Bid2796tqqoqlZWVSfp22tuWLVvYl8EV5eXlVzwe+/7zftYcI3UZuPPOO7V3716dPn1aAwMDeu+991RdXe13s2C4EydO6Omnn1Zzc7Pq6+slSfPmzdPhw4d19OhRjY6OqqOjg1qDI6+99po6OjrU3t6uZ555Rvfcc482btyo4uJi9fT0SPr2jnLUGZy4++67tXv3bp09e1ajo6P66KOPVFdXx/4MWTN37lx1d3frwoULsixLu3bt0m233ca+DK4Y73hsxowZgak5RuoycN1112nlypVqbGzU8PCwHnzwQd1yyy1+NwuG27Jli4aGhrRp06bkc0uXLtWmTZu0YsUKDQ0NqaamRnV1dT62ErmqublZa9euVX9/vyoqKtTY2Oh3k2CwefPm6fHHH9eyZcs0PDysu+66Sw899JBuuOEG9mfIigULFujTTz9VNBpVUVGRbr75Zj3xxBO677772Jch64qLi8c9HgtK/xmyLMvy5ZMBAAAAAI4x/RIAAAAADEaoAwAAAACDEeoAAAAAwGCEOgAAAAAwGKEOAAAAAAxGqAMAAAAAgxHqAAAAAMBghDoAAAAAMNj/BzdUzcZvNZtWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = [1.0*k for k in range (steps)] \n",
    "plt.figure(figsize = [15,4])\n",
    "plt.scatter(t,H)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}