membrane fluctuations MCMC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import quad\n",
    "import langevin\n",
    "import langevin_cached_model as lcm\n",
    "import pymc3 as pm\n",
    "import lmfit as lm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# here we define the values of the physical parameters\n",
    "# all energies are in kT = 4e-21J and lengths are in micrometers\n",
    "# kappa = 10kT, sigma = 0.5 microJ/m^2, gamma = 0.1MJ/m^4, eta = 1.0mPa s\n",
    "kappa = 10\n",
    "sigma = 125\n",
    "ga = 25\n",
    "eta= 0.25\n",
    "qmin=np.sqrt(3)/10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we describe the system in a continuous way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gamma(q, kappa=kappa, sigma=sigma, ga=ga, eta=eta ):\n",
    "    return (kappa*q**4+sigma*q**2+ga)/(4*eta*q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def avg(q,resolution):\n",
    "    return (np.exp(-gamma(q)*resolution)-1.0+gamma(q)*resolution)/(gamma(q)**2*resolution**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def integrant(q,t):\n",
    "    return q*np.exp(-gamma(q)*t)/(kappa*q**4+sigma*q**2+ga)*avg(q,1e-6)*np.exp(-0.25*0.5**2*q**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "time=np.logspace(-5,0,50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "acf=[quad(integrant,qmin,np.inf,args=(t))[0] for t in time]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11e0dcfd0>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEACAYAAACznAEdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8leWZ//HPdXKyQAIJkACBAGEJKIILpoCoaF0QrFOq\ntYpaabUVscVOOzqOnc6vM7+ZcXSmtk6pjhZbWulo1bHWotWqVXGjKAEVUbbIGowQ9iVkv+aPnHbS\nFJITcpLnLN/363VeyfM8931y3RzIl/tZzd0REREJBV2AiIjEBwWCiIgACgQREYlQIIiICKBAEBGR\nCAWCiIgACgQREYlQIIiICKBAEBGRCAWCiIgAEA66gI7Iz8/34uLioMsQEUkoK1as2OXuBe21S6hA\nKC4upqysLOgyREQSipltiaaddhmJiAigQBARkQgFgoiIAAoEERGJUCCIiAigQBARkYiEOu30eH34\n8QE+3neE9HCI9JA1f00LkZ5mZKSFyAiHyAynRb42L4dDhpkFXbqISLdJiUB45O0t/PeyrR3qEzLI\nDKeRlR6iR3oaWelpZKY3L/fMSKNnRpjsjDR6RL72zAyTk5lGTmY6vbLCkVc6vbPC9O6RTm6PdLLS\n07pohCIinZcSgfD1T4/iytKh1DU2Uf9nL6euoan51dhEbX0jdY3Ny7UNTdTUN1JT3/z1SKvvdx+q\n5kh9I4drG6mua6C6rrHdOrLSQ/TpmUFuj3T69MygT3Y6/bIzKeiVSX5OJvk5GX/6fkDvLDLC2qMn\nIt0nJQKhMLcHhbk9uvRnNDU51fWNHKyp51BNAwdqGjhYU8/BmgYO1NSzr7qefdV17KuuZ2/k+7Wf\nHGT3od3sP1L/F+9nBvk5mQzKzaIwtweD8nowKC+Loj49GZ6fzbB+PTXjEJGYiioQzGw68EMgDfiJ\nu9/VartFtl8MVANfdveV7fU1s5uBrwONwG/d/bZOjyggoZCRkxkmJzMMuR3rW9vQyO5Ddew6VMuu\nQ7VUHaylcn8Nlftq+Hj/ETbsPMhrG6r+YhZSmJvFsH7NATGyIIcTBvZmzMBeFPTKjOHIRCRVtBsI\nZpYG3AdcCFQAy81ssbt/2KLZDKAk8poE3A9MaquvmX0amAmc4u61ZtY/lgNLJJnhtMgM4NizGHdn\n/5F6tu6pZvPuajbvOszmXYfZtPswv1v9CXur/2+WkZ+TwQkDe3PCwF6cWNib04bmMTw/WwfJRaRN\n0cwQJgLl7r4RwMwepfkXectAmAkscncHlplZnpkVAsVt9L0JuMvdawHcfWdshpSczIy8nhnk9czg\n5KK8v9i++1At6z45yJpPDrK28gDrdhzkF8u2UNvQBECfnumcNrQPE4bmMWFoH04Zkkd2ZkrsMRSR\nKEXzG2EwsK3FcgXNs4D22gxup+9o4GwzuwOoAW519+XRly4t9cvJZMqoTKaMyv/TusYm56OqQ6zc\nspeVW/eycus+Xl7bnLtpIePUIXlMLSlg6uh8Ti7KIy2kGYRIKgvyv4hhoC8wGfgU8LiZjYjMMv7E\nzOYAcwCGDh3a7UUmsrSQMXpAL0YP6MWsic1/dvur63ln217KNu/l9Q1V/OdL67nn9+vJ7ZHOWaPy\nmTo6n/NPHEB+jo5DiKSaaAJhOzCkxXJRZF00bdLb6FsBPBkJgLfNrAnIB6pavrG7LwAWAJSWlv5Z\nWEjH5fZM59wx/Tl3TH9uvWgMew7X8Ub5Ll5bX8XrG6r47fuVhOx9pozM55KTC5k+biB5PTOCLltE\nuoG1+g/5XzYwCwPrgfNp/mW+HLja3T9o0eYzwDyazzKaBMx394lt9TWzucAgd/+umY0GXgKGtp4h\ntFRaWup6QE7XcXfWVB7k2fcreWbVx2zeXU04ZJxdks8lJw/ionEDm8+iEpGEYmYr3L203XbtBULk\nzS4G/pPmU0cXuvsdkV/ouPsDkdNO7wWm03za6XXuXnasvpH1GcBC4FSgjuZjCC+3VYcCofu4O6u3\nH+CZVR/zzKpKtu87Qk5mmMsmDGb2GcMY1b9X0CWKSJRiGgjxQoEQDHdn5da9PPzWVp55r5K6xiam\njOzH7DOKueDE/oTTdEW1SDxTIEiX2H2olsfKtvHwsq1s33eEQblZzJ5SzLWTh+k0VpE4pUCQLtXQ\n2MRLa3ey6A+bebN8N/2yM7jp3JF8cfIw3VJDJM4oEKTbrNy6lx+8sJ43yncxoHcm8z49iis+NYTM\nsIJBJB4oEKTb/eGj3fzgxXUs37yXwXk9+OvzS7j89CJCuuBNJFDRBoKOBkrMnDGyH4/feAYPXT+R\n/JwMbvvVKj7/wFLWVB4IujQRiYICQWLKzDhndAFPff1M7rnyFLbsruaSH73Bnc+tobquIejyRKQN\nCgTpEmbGpacV8dLfnMPlE4r48asbmXbPa7yyVvcwFIlXCgTpUn2yM/j3y0/msTmTyUpP47qfL+fr\nD69k96HaoEsTkVYUCNItJo3ox7PfOJtbp43mxTU7+Mz8N1i+eU/QZYlICwoE6TYZ4RDzzivhyZum\nkJUeYtaCZdy/5COamhLnTDeRZKZAkG43bnAuT998FtNPGsi//24tX3loOXsP1wVdlkjKUyBIIHpl\npXPv1afxzzNP4s3y3Vw8/3VWbNEuJJEgKRAkMGbG7DOK+dVNUwinGVf+eBm/WLYl6LJEUpYCQQI3\nviiXZ24+m6mjC/h/T63m7ufXkUhX0IskCwWCxIXcHuksuPZ0riwdwr2vlHP7r96nobEp6LJEUoru\nVyxxI5wW4q7Pj6d/70x+9HI5uw/X8qOrJtAjQzfJE+kOmiFIXDEzbpk2hn+ZeRIvrd3JNT9ZpjOQ\nRLqJAkHi0rVnFHP/NRNY/fEBLn9gKRV7q4MuSSTpKRAkbk0fV8gvrp/IzoO1zFqwjMr9R4IuSSSp\nKRAkrk0a0Y9HvjqZfdX1fPEnb+keSCJdSIEgcW98US4//VIpFXuP8KWfvc2BmvqgSxJJSgoESQiT\nRvTjgS+eztrKg3z152UcqWsMuiSRpKNAkITx6RP6c8+Vp7J8yx5uengFdQ26TkEklqIKBDObbmbr\nzKzczG4/ynYzs/mR7avMbEJ7fc3sn8xsu5m9G3ldHJshSTL7q1MGceel41myropvPfYujbpTqkjM\ntHthmpmlAfcBFwIVwHIzW+zuH7ZoNgMoibwmAfcDk6Loe4+73x2z0UhKmDVxKAdrGrjj2TX0ygpz\n52XjMbOgyxJJeNFcqTwRKHf3jQBm9igwE2gZCDOBRd58A5plZpZnZoVAcRR9RTrshqkj2Hekjvte\n+YjRA3px/VnDgy5JJOFFs8toMLCtxXJFZF00bdrre3NkF9NCM+sTddUiwC0XjmHa2AHc8ewaln60\nK+hyRBJekAeV7wdGAKcClcD3j9bIzOaYWZmZlVVVVXVnfRLnQiHj+1ecwvD8bOY98o6uZhbppGgC\nYTswpMVyUWRdNG2O2dfdd7h7o7s3AQ/SvGvqL7j7AncvdffSgoKCKMqVVNIrq/kuqfWNTdz4ixU6\nHVWkE6IJhOVAiZkNN7MMYBawuFWbxcDsyNlGk4H97l7ZVt/IMYY/uhRY3cmxSIoaUZDDD2edyoeV\nB/j2k6v0LAWR49TuQWV3bzCzecDzQBqw0N0/MLO5ke0PAM8CFwPlQDVwXVt9I2/9H2Z2KuDAZuDG\nWA5MUst5JwzglgtHc/cL6xk3OJevnj0i6JJEEo4l0v+mSktLvaysLOgyJE65Ozf990pe+PATFl0/\nibNK8oMuSSQumNkKdy9tr52uVJakYWbcfcUpjOqfw7xfrtTdUUU6SIEgSSUnM8yPry2lrqGJW//n\nPZp0JbNI1BQIknSG52fzD58Zy5vlu3noD5uDLkckYSgQJCldNXEI553Qn7ueW0v5zoNBlyOSEBQI\nkpTMjLs+P57szDDffOxd3RlVJAoKBEla/Xtl8W+Xjmf19gP86OUNQZcjEvcUCJLUpo8byOWnF3Hf\nK+Ws2LI36HJE4poCQZLeP/7VWApze3DL4+9yuLYh6HJE4pYCQZJer6x0fnDFKWzZU80dz64JuhyR\nuKVAkJQwaUQ/5pw9gkfe2sqSdTuDLkckLikQJGX8zbTRjCzI5ru/+YCaet0VVaQ1BYKkjMxwGv8y\ncxxb91Rz/5KPgi5HJO4oECSlTBmVz8xTB3H/ko/YtOtw0OWIxBUFgqSc71x8IpnhEN/9zWo9O0Gk\nBQWCpJz+vbO4ZdpoXt+wi2ff/yTockTihgJBUtIXJw/jpEG9+ednPuCQrk0QARQIkqLCaSH+9XPj\n2Hmwlv98cX3Q5YjEBQWCpKzThvZh1qeG8rOlm1lTeSDockQCp0CQlHbbRWPI7ZHOPzy1Wg/TkZSn\nQJCU1ic7g9tnnMCKLXt5YmVF0OWIBEqBICnv8glFnDY0j++/sI7qOh1gltSlQJCUFwoZf3/xiew4\nUMvCNzYFXY5IYBQIIsCnivty4dgBPPDqRnYfqg26HJFARBUIZjbdzNaZWbmZ3X6U7WZm8yPbV5nZ\nhA70vcXM3MzyOzcUkc75u+ljqK5r4Ecvlwddikgg2g0EM0sD7gNmAGOBq8xsbKtmM4CSyGsOcH80\nfc1sCDAN2NrpkYh00qj+vbjyU0N5+K0tbNmt+xxJ6olmhjARKHf3je5eBzwKzGzVZiawyJstA/LM\nrDCKvvcAtwE630/iwrcuKCEcCvG959cFXYpIt4smEAYD21osV0TWRdPmmH3NbCaw3d3fa+uHm9kc\nMyszs7KqqqooyhU5fv17Z/HVs4fzzKpK3tu2L+hyRLpVIAeVzawn8PfAd9tr6+4L3L3U3UsLCgq6\nvjhJeXOmjqBfdgZ3PrdGd0OVlBJNIGwHhrRYLoqsi6bNsdaPBIYD75nZ5sj6lWY2sCPFi3SFXlnp\nfOP8EpZt3MOSdZqVSuqIJhCWAyVmNtzMMoBZwOJWbRYDsyNnG00G9rt75bH6uvv77t7f3YvdvZjm\nXUkT3F33Ipa4cNXEoQzr15O7nltLo25pISmi3UBw9wZgHvA8sAZ43N0/MLO5ZjY30uxZYCNQDjwI\nfK2tvjEfhUiMZYRD3HbRCazbcZBf6ZYWkiIskfaRlpaWellZWdBlSIpwdz73X0upOlDDkr/9NBlh\nXccpicnMVrh7aXvt9Ddc5BjMjG9dUMLH+2t4YoVmCZL8FAgibThndAGnDMnjvlfKqW9sCrockS6l\nQBBpg5nxzfNL2L7vCE/qWIIkOQWCSDvOHVPAyUW53KtZgiQ5BYJIO8yMb5xXwrY9R3jqndaX4Igk\nDwWCSBTOP7E/Jw3qzb2vlNOgWYIkKQWCSBTMjG+cX8KW3dX85t2Pgy5HpEsoEESiNG3sAE4s1CxB\nkpcCQSRKZsZfnz+KTbsO88yqyqDLEYk5BYJIB0wbO5AxA3ox/+UNuseRJB0FgkgHhELNxxI2Vh3m\nmVU6liDJRYEg0kEzxg2kpH8OP3q5nCbNEiSJKBBEOigUMuadN4rynYf4/ZodQZcjEjMKBJHj8Jnx\nhQzO68GC1zYGXYpIzCgQRI5DOC3EDWcPp2zLXso27wm6HJGYUCCIHKcrPjWEvJ7p/FizBEkSCgSR\n49QzI8zsycN48cMdlO88FHQ5Ip2mQBDphNlTiskMh/jJ65olSOJTIIh0Qn5OJl8oLeLJldvZeaAm\n6HJEOkWBINJJXz1rBA1NTfxs6eagSxHpFAWCSCcV52czY1wh/71sCwdr6oMuR+S4KRBEYmDO1BEc\nrGng0be3BV2KyHFTIIjEwClD8pg8oi8L39xEXYNujS2JKapAMLPpZrbOzMrN7PajbDczmx/ZvsrM\nJrTX18z+JdL2XTN7wcwGxWZIIsG48ZyRVO6v4en3dNM7SUztBoKZpQH3ATOAscBVZja2VbMZQEnk\nNQe4P4q+33P3k939VOAZ4LudH45IcM4dXcCYAb348Wsf4a6b3kniiWaGMBEod/eN7l4HPArMbNVm\nJrDImy0D8syssK2+7n6gRf9sQP+CJKGZGTeeM4L1Ow6xZH1V0OWIdFg0gTAYaHmkrCKyLpo2bfY1\nszvMbBtwDceYIZjZHDMrM7Oyqir9I5P4dsnJg+jfK5Ofvbk56FJEOizQg8ru/h13HwI8DMw7RpsF\n7l7q7qUFBQXdW6BIB2WEQ8w+Yxivra9iw46DQZcj0iHRBMJ2YEiL5aLIumjaRNMXmgPh81HUIhL3\nrp40jMxwiIWaJUiCiSYQlgMlZjbczDKAWcDiVm0WA7MjZxtNBva7e2Vbfc2spEX/mcDaTo5FJC70\nzc7gsgmDeXJlBXsP1wVdjkjU2g0Ed2+geXfO88Aa4HF3/8DM5prZ3EizZ4GNQDnwIPC1tvpG+txl\nZqvNbBUwDfjr2A1LJFjXnTmc2oYmHnl7a9CliETNEun0uNLSUi8rKwu6DJGoXPvTt1i/4yCv33Ye\nGWFdAyrBMbMV7l7aXjv9LRXpItefOZwdB2p5bnVl0KWIREWBINJFzhldwIiCbH76xiZdqCYJQYEg\n0kVCIeO6M4ezqmI/K7fuDbockXYpEES60OcnDCa3Rzo/fWNT0KWItEuBINKFemaEuWriUH63+hMq\n9lYHXY5ImxQIIl1s9hnDMDMW/WFL0KWItEmBINLFBuX1YMa4gfzy7a0crm0IuhyRY1IgiHSD688a\nzsGaBp5YURF0KSLHpEAQ6QYThvbh1CF5PLR0M01NOgVV4pMCQaSbXHdmMRt3Hea1DbqNu8QnBYJI\nN5kxrpCCXpn8fOnmoEsROSoFgkg3yQiH+OKkYSxZV8VHVYeCLkfkLygQRLrR1ZOGkpEWYpFmCRKH\nFAgi3aigVyaXnFLIEysqOFBTH3Q5In9GgSDSza6bMpzDdY38T5lOQZX4okAQ6Wbji3I5fVgfHlq6\nmUadgipxRIEgEoDrzixm655qlqzbGXQpIn+iQBAJwEUnDWRg7yx+9ubmoEsR+RMFgkgA0tNCXHvG\nMN4o38WGHQeDLkcEUCCIBOaqiUPJCId0oZrEDQWCSED6ZmfwuVMH8eTK7eyv1imoEjwFgkiAvjSl\nmCP1jTxWtjXoUkSiCwQzm25m68ys3MxuP8p2M7P5ke2rzGxCe33N7HtmtjbS/tdmlhebIYkkjpMG\n5TJpeF8eWrqFhsamoMuRFNduIJhZGnAfMAMYC1xlZmNbNZsBlERec4D7o+j7IjDO3U8G1gPf7vRo\nRBLQ9WcNZ/u+I7zw4Y6gS5EUF80MYSJQ7u4b3b0OeBSY2arNTGCRN1sG5JlZYVt93f0Fd//j46OW\nAUUxGI9IwrngxAEM7duThW9sCroUSXHRBMJgYFuL5YrIumjaRNMX4HrguShqEUk6aSHjy1OKKduy\nl/e27Qu6HElhgR9UNrPvAA3Aw8fYPsfMysysrKpKDxaR5PSF0iJyMsMsfFOzBAlONIGwHRjSYrko\nsi6aNm32NbMvA5cA17j7UW/q4u4L3L3U3UsLCgqiKFck8fTKSueK0iH8dlUln+yvCbocSVHRBMJy\noMTMhptZBjALWNyqzWJgduRso8nAfnevbKuvmU0HbgM+6+7VMRqPSML68pRiGt35xbLNQZciKard\nQIgc+J0HPA+sAR539w/MbK6ZzY00exbYCJQDDwJfa6tvpM+9QC/gRTN718weiN2wRBLP0H49mTZ2\nAA+/tZUjdY1BlyMpyI6xpyYulZaWellZWdBliHSZtzbu5soFy/i3S8dz9aShQZcjScLMVrh7aXvt\nAj+oLCL/Z+Lwvpw0qDcL39xEIv1nTZKDAkEkjpgZXzlrOOU7D/Hahl1BlyMpRoEgEmc+c3IhBb0y\ndaGadDsFgkicyQynce3kYby6vorynXpWgnQfBYJIHLpmUvOzEhbqiWrSjRQIInGoX04ml546mF+t\nqGDXodqgy5EUoUAQiVM3TB1BXWMTP9csQbqJAkEkTo3qn8O0sQNY9IfNHKptaLe9SGcpEETi2Nxz\nRnKgpoFH39YT1aTrKRBE4thpQ/sweURffvL6Juoa9EQ16VoKBJE4d9O5o/jkQA1Pvdv6JsMisaVA\nEIlzU0vyGVvYmwde/YimJt3OQrqOAkEkzpkZc88dycaqw7y4Rs9dlq6jQBBJABePG8iQvj24f8lH\nuumddBkFgkgCCKeFmDN1JO9u28dbm/YEXY4kKQWCSIL4wulF5Odk8MCrHwVdiiQpBYJIgshKT+PL\nU4pZsq6KDz8+EHQ5koQUCCIJ5NrJxWRnpPHj1zRLkNhTIIgkkNye6VwzeRhPv/cxG6sOBV2OJBkF\ngkiCueHsEWSG07jn9xuCLkWSjAJBJMEU9Mrk+rOKefq9j3UsQWJKgSCSgOacPZLeWWG+/8K6oEuR\nJKJAEElAuT3TufGckby0dicrtui6BImNqALBzKab2TozKzez24+y3cxsfmT7KjOb0F5fM/uCmX1g\nZk1mVhqb4YikjuvOLCY/J5P/+N06Xb0sMdFuIJhZGnAfMAMYC1xlZmNbNZsBlERec4D7o+i7GrgM\neK3zwxBJPT0zwtx83ije2rSHN8p3BV2OJIFoZggTgXJ33+judcCjwMxWbWYCi7zZMiDPzArb6uvu\na9xdO0BFOmHWxCEMzuvB957XLEE6L5pAGAxsa7FcEVkXTZto+orIccoMp/HNC0pYVbGf5z/QnVCl\nc+L+oLKZzTGzMjMrq6qqCrockbhz6WmDGVmQzfdfWEejnpcgnRBNIGwHhrRYLoqsi6ZNNH3b5O4L\n3L3U3UsLCgo60lUkJYTTQtwybQwbdh7iN3qqmnRCNIGwHCgxs+FmlgHMAha3arMYmB0522gysN/d\nK6PsKyKdNP2kgYwb3Jt7fr9ez16W49ZuILh7AzAPeB5YAzzu7h+Y2Vwzmxtp9iywESgHHgS+1lZf\nADO71MwqgDOA35rZ8zEdmUgKCYWMW6eNYdueIyx8c1PQ5UiCskQ6M6G0tNTLysqCLkMkbs1ZVMZr\nG6p48VvnMKRvz6DLkThhZivcvd3rveL+oLKIRO+fPnsSaWb8w1OrdRqqdJgCQSSJDMrrwa0XjeHV\n9VU8s6oy6HIkwSgQRJLM7DOKObkol///9Ifsr64PuhxJIAoEkSSTFjL+7dLx7K2u467frQ26HEkg\nCgSRJDRucC7Xn1nML9/eStlm3Q1VoqNAEElS37xgNIPzevDtJ9/XtQkSFQWCSJLKzgzzzzNPYsPO\nQzz4+sagy5EEoEAQSWLnnziAi8cP5IcvbWDTrsNBlyNxToEgkuT+8a9OIisc4msPr6S6riHociSO\nKRBEktyA3ln88KrTWPvJAf72f1bpgjU5JgWCSAr49Jj+3D79BH77fiX3vVIedDkSp8JBFyAi3WPO\n1BGs/eQgd7+wnjEDe3Ph2AFBlyRxRjMEkRRhZtx52XhOLsrlm4++w/odB4MuSeKMAkEkhWSlp7Hg\n2lJ6Zoa5YVEZ+6rrgi5J4ogCQSTFDMzN4oEvnk7lvhrmPfIODY26aE2aKRBEUtDpw/rwr5eO443y\nXfz9r99XKAigg8oiKeuK0iFU7Klm/svlfHKglnuvPo3eWelBlyUB0gxBJIX9zbQx3HXZeJaW7+Ly\n+5eybU910CVJgBQIIilu1sShPHT9RCr313Dpf73Jyq17gy5JAqJAEBHOHJXPr792Jj0zwsxasIyn\n3/s46JIkAAoEEQFgVP8cnvr6mZw8OJebf/kOdz+/jiN1jUGXJd1IgSAif9I3O4OHb5jEZRMGc+8r\n5Uz93is8tHQztQ0KhlSgQBCRP5MZTuMHV5zK4zeewfD8bP5x8Qecd/erPL58m05PTXJRBYKZTTez\ndWZWbma3H2W7mdn8yPZVZjahvb5m1tfMXjSzDZGvfWIzJBGJhYnD+/LYnMksun4i/XIyuO1Xq5h2\nz2v8+p0KDtTUB12edAFr71a4ZpYGrAcuBCqA5cBV7v5hizYXAzcDFwOTgB+6+6S2+prZfwB73P2u\nSFD0cfe/a6uW0tJSLysrO86hisjxcnde+HAHP3hhPet2HCQtZJxSlMtZJQWcXZLPqUPySE/TDod4\nZWYr3L20vXbRXJg2ESh3942RN34UmAl82KLNTGCRN6fLMjPLM7NCoLiNvjOBcyP9HwKWAG0GgogE\nw8y46KSBXHjiAN7evIc3Nuzi9fJd3PvyBua/tIGczDCThvdlREE2hbk9GJTXg0F5WRTm9qBfdgah\nkAU9BIlCNIEwGNjWYrmC5llAe20Gt9N3gLtXRr7/BNC9eEXiXChkTB7Rj8kj+nHrRWPYX13P0o+a\nw+Gtjbt5o3wXtQ1/fpwhIy1E7x5hstLTyEpPo0fklZkeIiMtRChkhAxCZpHvDQMskiF/jBKLrEjV\naLlh6ghOLOzdpT8jLm5d4e5uZkfdd2Vmc4A5AEOHDu3WukSkbbk905kxvpAZ4wuB5l1Lew7XUbm/\nho/3HWn+uv8Ih2oaOFLfSE19IzX1TRypa+RATQMNjU00eXO/xianyR13aIrsyv7jL4U/7tl2Uvdp\nb1ccGdLlPyOaQNgOtKykKLIumjbpbfTdYWaF7l4Z2b2082g/3N0XAAug+RhCFPWKSEDMjH45mfTL\nyWTc4Nygy5EOiuYo0HKgxMyGm1kGMAtY3KrNYmB25GyjycD+yO6gtvouBr4U+f5LwG86ORYREemE\ndmcI7t5gZvOA54E0YKG7f2BmcyPbHwCepfkMo3KgGriurb6Rt74LeNzMvgJsAa6I6chERKRD2j3t\nNJ7otFMRkY6L9rRTnTgsIiKAAkFERCIUCCIiAigQREQkQoEgIiJAgp1lZGZVNJ+imgvsb7Gp5fKx\nvs8HdsWgjNY/uzNtj7X9aOtTbcwdWU7EMXf0M269HM9jjtXf69bLGvPxj3mYuxe028rdE+4FLDjW\nchvfl3XFz+5M22NtP9r6VBtzR5YTccwd/YwTacyx+nutMXfNmNt6Jeouo6fbWD7W9131szvT9ljb\nj7Y+1cbckeVEHHNHP+PWy/E85lj9vW69rDF3sYTaZdQZZlbmUVyYkUw05tSgMaeG7hhzos4QjseC\noAsIgMacGjTm1NDlY06ZGYKIiLQtlWYIIiLSBgWCiIgACgQREYlQIABmdq6ZvW5mD5jZuUHX013M\nLNvMysw1ujXXAAACyUlEQVTskqBr6Q5mdmLkM37CzG4Kup7uYGafM7MHzewxM5sWdD1dzcxGmNlP\nzeyJoGvpSpF/uw9FPttrYvW+CR8IZrbQzHaa2epW66eb2TozKzez29t5GwcOAVlARVfVGisxGjPA\n3wGPd02VsRWLMbv7GnefS/PDmM7synpjIUZjfsrdbwDmAld2Zb2dFaPxbnT3r3RtpV2jg+O/DHgi\n8tl+NmY1JPpZRmY2leZf5ovcfVxkXRqwHriQ5l/wy4GraH5q252t3uJ6YJe7N5nZAOAH7h6zxO0K\nMRrzKUA/mkNwl7s/0z3VH59YjNndd5rZZ4GbgF+4+yPdVf/xiNWYI/2+Dzzs7iu7qfwOi/F4n3D3\ny7ur9ljo4PhnAs+5+7tm9oi7Xx2LGtp9hGa8c/fXzKy41eqJQLm7bwQws0eBme5+J9DW7pG9QGZX\n1BlLsRhzZNdYNjAWOGJmz7p7U1fW3Rmx+pzdfTGw2Mx+C8R1IMToczaaH1f7XDyHAcT833LC6cj4\naQ6HIuBdYrinJ+ED4RgGA9taLFcAk47V2MwuAy4C8oB7u7a0LtOhMbv7dwDM7MtEZkhdWl3X6Ojn\nfC7NU+1Mmp8Dnog6NGbgZuACINfMRnnzM9ATSUc/437AHcBpZvbtSHAksmONfz5wr5l9hhje3iJZ\nA6FD3P1J4Mmg6wiCu/886Bq6i7svAZYEXEa3cvf5NP/ySAnuvpvm4yVJzd0PA9fF+n0T/qDyMWwH\nhrRYLoqsS2Yas8acjFJtvK116/iTNRCWAyVmNtzMMoBZwOKAa+pqGrPGnIxSbbytdev4Ez4QzOyX\nwB+AMWZWYWZfcfcGYB7wPLAGeNzdPwiyzljSmDVmknDMqTbe1uJh/Al/2qmIiMRGws8QREQkNhQI\nIiICKBBERCRCgSAiIoACQUREIhQIIiICKBBERCRCgSAiIoACQUREIv4XN+noCuTm6EkAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x117be6f60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogx(time,acf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "lets try the discrete one.  Here we assume that we have a $q_{min}$ and that $q^2=q_{x}^2+q_{y}^2$ where each $q_{x}=lq_{min}$ and $q_{y}=mq_{min}$ which results in $q_{lm}=q_{min}\\sqrt{l^2+m^2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0,\n",
       " 1.0,\n",
       " 1.4142135623730951,\n",
       " 2.0,\n",
       " 2.0,\n",
       " 2.2360679774997898,\n",
       " 2.2360679774997898,\n",
       " 2.8284271247461903,\n",
       " 3.0,\n",
       " 3.0,\n",
       " 3.1622776601683795,\n",
       " 3.1622776601683795,\n",
       " 3.6055512754639891,\n",
       " 3.6055512754639891,\n",
       " 4.2426406871192848]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q_list=[np.sqrt(l**2+m**2) for l in range(0,4) for m in range(0,4)]\n",
    "sorted(q_list)[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15\n",
      "1.0 160.0\n",
      "1.0 160.0\n",
      "1.41421356237 222.738636074\n",
      "2.0 342.5\n",
      "2.0 342.5\n",
      "2.2360679775 402.49223595\n",
      "2.2360679775 402.49223595\n",
      "2.82842712475 588.666395338\n",
      "3.0 653.333333333\n",
      "3.0 653.333333333\n",
      "3.16227766017 719.418167688\n",
      "3.16227766017 719.418167688\n",
      "3.60555127546 926.349327696\n",
      "3.60555127546 926.349327696\n",
      "4.24264068712 1299.89796608\n"
     ]
    }
   ],
   "source": [
    "q_list=np.array(sorted(q_list)[1:])\n",
    "print(q_list.size)\n",
    "for q in q_list:\n",
    "    print(q,gamma(q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def corrfct(t,q_list):\n",
    "    corr=np.exp(-gamma(q_list)*t)/(kappa*q_list**4+sigma*q_list**2+ga)*avg(q_list,1e-6)*np.exp(-0.25*0.5**2*q_list**2)\n",
    "    return corr.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "acf_d = np.array([corrfct(t,q_list) for t in time])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11e4efc50>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEACAYAAACznAEdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOXd//H3dyZ7AgmQAIFAwi4oe2QVFS0KuKBWLSju\nSnFfalv72J/t87RVazfl0WqxLsW11K2oqHXfWYLIJqABWcIatrCE7Pfvj5n2SVNIJjDJmZl8XteV\ni5wz95187iuaT86cM2fMOYeIiIjP6wAiIhIZVAgiIgKoEEREJEiFICIigApBRESCVAgiIgKoEERE\nJEiFICIigApBRESCVAgiIgJAnNcBGiMzM9Pl5eV5HUNEJKosWrRoh3Muq6FxUVUIeXl5FBQUeB1D\nRCSqmNn6UMbpKSMREQFUCCIiEqRCEBERQIUgIiJBKgQREQFUCCIiEhRVl50eqTXL5rF327f44hLw\nx8Xji0vEH5+APz4Rf1wicYlJxCUmEx+fRHxSMomJySTEJ+Dzqy9FpOVoEYWw44OHGb7zlUbNqXZG\nKQmUWwLlJFJhCVRYIlW+BCp8yVT6U6jyp1ATl0x1fAouPg0S0rCkVviT0/GnZBCfkk5iWhuSW7ch\nLaM96a1akRCnkhGRyNQiCiHv3Lv4ZsfV1FSVU11RQU1VReDzqkpcVTk1lWW4qnJcZRmuqgKqy6Gq\nDKv14a8uw1ddRlx1GYnVZbSqKCGh5iBJroxkykihrMEcB10CW2jFfl8aB/zplMWlU5HYlqqU9lha\nFv7WHUjM6ERq24607ZhL+4xWxOkoRUSaSYsohA45PeiQ06Npv0lNDVVl+yjdt4fSfbso27+Hiv17\nqCzdQ+WB3bjS3bjSXfjKdhNXsYfUihKyKr+lVdliWpfs/88v54wdpFPsy2RvfHtKk7OpatUJf9tu\npGb3pn1uHzpntSMp3t+06xKRFiOkQjCz8cADgB/4s3Pu3jqPW/DxiUApcLlz7ouG5prZjcD1QDXw\nunPuR0e9Iq/4fMSlpNM6JZ3WHXIbN7eqnLKSrezbsZn9OzdTtnsLlbs3YXuLSDywma7lG2lb8gXJ\nJWVQBCwNTNvs2rLF14k9yV0oy+hBQqf+tO02iB7dupORkhD2JYpIbGuwEMzMDzwEjCPw62ihmc1x\nzn1Va9gEoFfwYzjwMDC8vrlmNhaYBAx0zpWbWftwLiyqxCWS1C6XpHa5HPbuU87hDu6mZPM37C5a\nRdnWb3A719Bu/zp6HfyU1qVvwGagAIpda+b78tiZ2ouq9seS3nMkvfsNIjsjpRkXJSLRJpQjhGFA\noXNuLYCZPU/gF3ntQpgEzHLOOWCemWWYWTaQV8/ca4F7nXPlAM657eFZUowyw1LaktFzOBk9h//H\nw25/MbvXfcmutYup3rKc7N2rGbT/VRL3vwhrYddbaXzi68PONgPxdR1Gp2NPoH+3zjrJLSL/Ekoh\ndAY21touInAU0NCYzg3M7Q2MMbNfAWXA7c65haFHl9osLYu2x42j7XHj/m9nTTWV21axdcVHlH07\nj57Fizlh1+Ow63GqFvtYSk/WtxlJXO/v0HfISfTo0JrAs38i0hJ5eVI5DmgLjACOB2abWffgUca/\nmNk0YBpA165dmz1kVPP5ic8+li7ZxxI4IAMO7mZv4TyKv/qA9us/ZNCeWfgW/IU981N5zz+Q3R3H\nkJV/DiMG9CExTiesRVqSUAphE9Cl1nZOcF8oY+LrmVsEvBQsgAVmVgNkAsW1v7BzbiYwEyA/P//f\nykKOQHIbWvefQOv+EwLbB3ayY9lblCx7k6FbPyFj82dU//0+Fsw5jo2dJtB55IUM69eDeF3+KhLz\nrM4f5P85wCwO+Bo4lcAv84XARc65FbXGnAHcQOAqo+HADOfcsPrmmtl0oJNz7i4z6w28C3Ste4RQ\nW35+vtMb5DQh56jcvJRNnz1H6jdzyKrYRKXzM88GsKXLRHqeOJnBPbvoaSWRKGNmi5xz+Q2Oa6gQ\ngl9sInA/gUtHH3fO/Sr4Cx3n3CPBy04fBMYTuOz0CudcweHmBvcnAI8Dg4AKAucQ3qsvhwqhGTlH\n+cYv2PzpM6SveY22VdvY55L5IOkU/MOnMfaEMSQn6CklkWgQ1kKIFCoEjzhH2befs+Xdh+m8aS4J\nVDGf4yjqOZVh4y+mS2ZrrxOKSD1UCNIk3P5iNr77J9KWzaJt1TY2uXbMa3cex537A/p0yfY6nogc\nggpBmlZ1Fbu+nMO+j/5IbslCdrjWfNzhEoacdxu5HTO9TicitYRaCLp0RI6MP462Q88j99Z32Df1\nDfZnHMO52x8i8eF8/j7zv9myc4/XCUWkkVQIctRa9RxF3q1vs/vClyhP68qkzb+nZsZQXn3y1+w7\nWO51PBEJkQpBwqZNv1PJvf1Dis95luqUTM5adzfr7juB+Z9/6HU0EQmBCkHCy4ysQWfQ9UfzWHfi\nH+jCVoa+eQ5vz5jOjl27vE4nIvVQIUjTMCPvlCtJuXUxqzuexbhdz1E+Yxifzn2WaLqQQaQlUSFI\nk0poncmx185i46QXqPEnMXrBtSy472y2by3yOpqI1KFCkGbRZfA4Ot+xiMU9b2BQ6ee4R8aw4vM3\nvY4lIrWoEKTZ+OITGTz1V2y94FUqLZE+b05h4dN34WqqvY4mIqgQxAO5x42k9S2f8WXaGI4vfICv\nfjeR0j16fyQRr6kQxBOt09sy5LZX+LDnj+m5v4ADD4xk87IPvI4l0qKpEMQzPr+Pk6b+F19NeIFy\n5yPrhfP4+vX7vY4l0mKpEMRzg0eMxaZ/xBcJQ+i98GcUPv9j0KWpIs1OhSARoXPHbI655VXeSTqN\nnqseYd0TV0F1ldexRFoUFYJEjPTUZEbe8iwvt7qIvA0vUvTIeVBR6nUskRZDhSARJTUpngk3PcjT\nbW+k0/aP2Prg6VCqW16INAcVgkScpHg/37v+f/hzp5/TpmQlO/93LG73eq9jicQ8FYJEpHi/j6uu\nuZnHuv2O+NLtlDxyOpRs8jqWSExTIUjE8vuM6ZdeyqzeM/CX7WHPzDPgwA6vY4nELBWCRDSfz7h2\nyvk83OlukvYXUfLo2VBW4nUskZikQpCI5/cZN15xGfel/5SU3avY98T5uvpIpAmoECQqJCf4uW7a\ntfwq8VZSty2k9JmLoarC61giMSWkQjCz8Wa22swKzeyOQzxuZjYj+PhSMxvS0Fwz+7mZbTKzL4Mf\nE8OzJIlVmWmJXHrNrfzSvk/K+vcon30V6E6pImHTYCGYmR94CJgA9AOmmFm/OsMmAL2CH9OAh0Oc\n+wfn3KDgx9yjXYzEvu5ZaZxx+Y+5t3oqiV/PoWrOzbrNhUiYhHKEMAwodM6tdc5VAM8Dk+qMmQTM\ncgHzgAwzyw5xrkijDM1ty8ALf8pDVZOI+/Ip3LyHvY4kEhNCKYTOwMZa20XBfaGMaWjujcGnmB43\nszYhp5YWb0L/bBLG3cVb1fm4f/wUvv3I60giUc/Lk8oPA92BQcAW4HeHGmRm08yswMwKiouLmzOf\nRLirT+zBaz3uYm1NR6r+ehns2eB1JJGoFkohbAK61NrOCe4LZcxh5zrntjnnqp1zNcCjBJ5e+g/O\nuZnOuXznXH5WVlYIcaWlMDN+ceFIfpr4Ew6WlVH93EW6HFXkKIRSCAuBXmbWzcwSgMnAnDpj5gCX\nBq82GgGUOOe21Dc3eI7hn84Flh/lWqQFykhJ4IcXn8Etlddj25bjXtVJZpEj1WAhOOeqgBuAt4CV\nwGzn3Aozm25m04PD5gJrgUICf+1fV9/c4Jz7zGyZmS0FxgK3hm9Z0pIMzW1L/rgp/K7yAmzZbPj8\nIa8jiUQlc1H011R+fr4rKCjwOoZEoJoax+VPLODi9T/lNH8BNvUl6DHW61giEcHMFjnn8hsap1cq\nS0zw+Yzff28Q9yTcxDrLwf3tCt0dVaSRVAgSMzLTErl7ykiuLLuZivIyeOVaqKnxOpZI1FAhSEwZ\n1SOTs8eeyM8rLoJvP4QFM72OJBI1VAgSc248pSdLss7hExuKe+dnULza60giUUGFIDEnzu/j3vMH\ncFvZVZS6JHjpGt0ZVSQEKgSJSQNyMjhr9GBuO3gFbFkCH93ndSSRiKdCkJh127jeLG99Im/Gn4r7\n+HewcYHXkUQimgpBYlZqYhy/PPc4bt83hX0JHeDl70P5fq9jiUQsFYLEtLF92nPKwJ5MPzANt+tb\n+MdPvY4kErFUCBLz7jqrHyvij2NO6ndh0RPwzTteRxKJSCoEiXmZaYnceUZffrjzLEpSu8HcH0Dl\nQa9jiUQcFYK0CBcMzWFo947cduAS2L0OPrnf60giEUeFIC2CmXH3ef35uLIfX6R/Bz75Pexc43Us\nkYiiQpAWo1tmKpePzmP69nOp9iXC3Nv13gkitagQpEW5/uSeVCS355nUS2DNe/DVK15HEokYKgRp\nUdJT4rnplF7899aR7GvTD978CZTv8zqWSERQIUiLM3VELjntWnFnxZW4fVvhg3u9jiQSEVQI0uIk\nxPm4Y/wxzNnZibVdvwvzHoatektvERWCtEjjj+vI0Nw2TN98Ji4pA16/TW+mIy2eCkFaJDPjvyb2\n5Zv9Cbydcx1snA9LnvU6loinVAjSYg3NbcMZ/bO5ZVU/KrOHwnu/hIoDXscS8YwKQVq0H43vQ2UN\n/Dn5Sti3Beb90etIIp5RIUiLltsulUtH5vGblW3Yl3cafPIAHNjhdSwRT4RUCGY23sxWm1mhmd1x\niMfNzGYEH19qZkMaMfcHZubMLPPoliJyZG48pSdpiXHcVzkZKg/AR7/xOpKIJxosBDPzAw8BE4B+\nwBQz61dn2ASgV/BjGvBwKHPNrAtwGrDhqFcicoQyUhK4ekx3nlqTxK4+k2HhY7BrrdexRJpdKEcI\nw4BC59xa51wF8Dwwqc6YScAsFzAPyDCz7BDm/gH4EaAbyoinLhuVR6ukOO49eC744+HdX3gdSaTZ\nhVIInYGNtbaLgvtCGXPYuWY2CdjknFtS3zc3s2lmVmBmBcXFxSHEFWm89OR4rhjdjdmrK9nR/2pY\n8RJsWuR1LJFm5clJZTNLAf4LuKuhsc65mc65fOdcflZWVtOHkxbrqtHdSEuM4549p0FKJrz9M90N\nVVqUUAphE9Cl1nZOcF8oYw63vwfQDVhiZuuC+78ws46NCS8STukp8Vw+Ko+XVu5l+5CbYd3H8M3b\nXscSaTahFMJCoJeZdTOzBGAyMKfOmDnApcGrjUYAJc65LYeb65xb5pxr75zLc87lEXgqaYhzbmu4\nFiZyJK46oRsp8X7u3j4C2nSDd34GNdVexxJpFg0WgnOuCrgBeAtYCcx2zq0ws+lmNj04bC6wFigE\nHgWuq29u2FchEiZtUhO4dFQef19WzJbjfwzbv4Ilz3kdS6RZmIui50jz8/NdQUGB1zEkxu3cX84J\nv36f0/u15/79t8O+bXDTYohL8DqayBExs0XOufyGxumVyiJ1tEtL5NKRucxZuoUtg2+FvUXw5TNe\nxxJpcioEkUO4ekx3EuJ8/GZNDnQeCh//HqorvY4l0qRUCCKHkNUqkanDc/n7ki1sG3wLlGzQuQSJ\neSoEkcOYdlJ34nzG777tCp0Gw0e/1VGCxDQVgshhtG+VxOTju/Dyl5vZffytsGc9LJ3tdSyRJqNC\nEKnHVSd0p7rG8ejW3tBxQOBOqNVVXscSaRIqBJF6dG2XwunHduTp+RsoG3077P4Wlv3N61giTUKF\nINKAq8d0Z29ZFc+X9IcO/XWUIDFLhSDSgKG5bRia24bHPltH9Yk/hF1rAndDFYkxKgSREFwzphsb\ndx3krep8aN8PPrxP9ziSmKNCEAnBuH4dyW2XwqOfrIOTfgQ7v4EVL3sdSySsVAgiIfD7jKtO6Mbi\nDXtYlDoGso4JHiXUeB1NJGxUCCIhOn9oDunJ8cz8eB2MuR12rIav3/A6lkjYqBBEQpSSEMfUEV35\nx1fbWNfxNEjvCp/O8DqWSNioEEQa4bKRecT7fDz++UYYdQNsnAcb5nkdSyQsVAgijdC+dRKTBnVi\ndsFGdve+EJLb6ChBYoYKQaSRrh7TnbLKGp5ZvAOOvwZWvw7FX3sdS+SoqRBEGqlPx1ac2DuLJz9b\nT/nQqyAuCT7/X69jiRw1FYLIEbjqhG7s2F/Om99Ww6CLYcnzsG+r17FEjooKQeQIjOmZSbfMVP7y\n2ToYeT3UVMH8R7yOJXJUVAgiR8DnMy4ZkcsXG/aw7GAm9D0bFj4OZXu9jiZyxFQIIkfou0NzSEnw\nM+vzdTD6JigvgS/+4nUskSOmQhA5QunJ8Zw7uDN/X7KZ3Rn9IW8MzHsYqiq8jiZyREIqBDMbb2ar\nzazQzO44xONmZjOCjy81syENzTWzXwTHfmlm/zCzTuFZkkjzuXRkHhVVNfy1YCOMvhn2boLlL3od\nS+SINFgIZuYHHgImAP2AKWbWr86wCUCv4Mc04OEQ5v7GOTfAOTcIeA246+iXI9K8+nRsxYjubXnq\n8/VUdz81cGvsTx8A57yOJtJooRwhDAMKnXNrnXMVwPPApDpjJgGzXMA8IMPMsuub65yrffYtFdD/\nQRKVLhuZx6Y9B3l31fbAUULxSih8x+tYIo0WSiF0BjbW2i4K7gtlTL1zzexXZrYRuJjDHCGY2TQz\nKzCzguLi4hDiijSvcf06kJ2exKzP18Ox50Fax8C5BJEo4+lJZefcnc65LsAzwA2HGTPTOZfvnMvP\nyspq3oAiIYjz+7h4eFc+KdxB4a5yGHY1rHkXtq/yOppIo4RSCJuALrW2c4L7QhkTylwIFMJ3Q8gi\nEpEmD+tKgt/HU5+vh6FXBm5nMV9HCRJdQimEhUAvM+tmZgnAZGBOnTFzgEuDVxuNAEqcc1vqm2tm\nvWrNnwTozymJWplpiZwxIJsXFhWxz98aBnwvcDuL0l1eRxMJWYOF4JyrIvB0zlvASmC2c26FmU03\ns+nBYXOBtUAh8ChwXX1zg3PuNbPlZrYUOA24OXzLEml+l43K40BFNS8v3gQjroWqMlj0hNexREJm\nLoouj8vPz3cFBQVexxA5rEkPfsL+8ireue0k7KlzoXgV3LwU4hK8jiYtmJktcs7lNzROr1QWCaNL\nRuaxpvgAn63ZCSOug31b4Ku/ex1LJCQqBJEwOnNANhkp8Twzfz30/A606wXzHtIL1SQqqBBEwigp\n3s8FQ3P4x4ptbN9fASOmw+bFsHGB19FEGqRCEAmzi4bnUlXj+OvCjTBwCiRlBI4SRCKcCkEkzLpl\npnJCz0yeW7CB6rgUGHo5rHwV9mzwOppIvVQIIk1g6oiubC4p4/1V22HYNYDBgplexxKplwpBpAmc\n2rcD7Vsl8vT89ZCeA/0mwaJZUL7f62gih6VCEGkC8X4fk4d15cOvi9m4qzRwCWp5CXz5rNfRRA5L\nhSDSRKYM64IBzy7YAF2Oh875sOBPUFPjdTSRQ1IhiDSR7PRkTu3bgdkLN1JeVQ3Dp8POQljzntfR\nRA5JhSDShKaOyGXngQreXL41cB4hrQPMf8TrWCKHpEIQaUJjemaS2y6FZ+ZvCNzPKP8qKHwbdnzj\ndTSR/6BCEGlCPp9x0bCuLPh2F19v2wf5V4A/QZegSkRSIYg0sfOH5pDg9/Hs/A2Q1h6O+27gaqOy\nEq+jifwbFYJIE2uXlsjE/h15cVERpRVVMPz7ULEfFj/jdTSRf6NCEGkGU0fksq+8ijlfboZOg6HL\n8OAlqNVeRxP5FxWCSDMYmtuGYzq24un563HOBS5B3b0Ovnnb62gi/6JCEGkGZsbFI3JZvmkvS4pK\noO9Z0KoTzH/Y62gi/6JCEGkm5w7uTGqCn6fnrQd/PBx/Faz9ALav8jqaCKBCEGk2aYlxnDO4M68u\n2cye0goYegX4EwPnEkQigApBpBlNHZFLeVUNLywqgtR2MOACWPI8HNztdTQRFYJIc+qb3ZqhuW14\nZv4GamocDPs+VJbCF095HU0ktEIws/FmttrMCs3sjkM8bmY2I/j4UjMb0tBcM/uNma0Kjn/ZzDLC\nsySRyHbJiFy+3XGAz9bshOwBkHtC4JXL1VVeR5MWrsFCMDM/8BAwAegHTDGzfnWGTQB6BT+mAQ+H\nMPdt4Djn3ADga+AnR70akSgwoX9H2qYm8NS8dYEdI66Fko2w6jVPc4mEcoQwDCh0zq11zlUAzwOT\n6oyZBMxyAfOADDPLrm+uc+4fzrl//kk0D8gJw3pEIl5inJ8L8nN4Z+V2tpaUQZ8J0CYP5ukSVPFW\nKIXQGdhYa7souC+UMaHMBbgSeCOELCIx4eJhudQ4x3MLNoDPH3ih2sZ5sGmR19GkBfP8pLKZ3QlU\nAYe8sYuZTTOzAjMrKC4ubt5wIk2ka7sUTuyVxfMLN1BZXQODLoaEVjpKEE+FUgibgC61tnOC+0IZ\nU+9cM7scOBO42DnnDvXNnXMznXP5zrn8rKysEOKKRIdLRuSybW85767cBkmtYcglsOJl2LvZ62jS\nQoVSCAuBXmbWzcwSgMnAnDpj5gCXBq82GgGUOOe21DfXzMYDPwLOds6Vhmk9IlFj7DHt6ZyRzFPz\n1gd2DJsWuNndwj97G0xarAYLIXji9wbgLWAlMNs5t8LMppvZ9OCwucBaoBB4FLiuvrnBOQ8CrYC3\nzexLM9P7CkqL4vcZU4Z14dPCnawt3g9tu8ExZ0DB41Chv5Gk+dlhnqmJSPn5+a6goMDrGCJhs31f\nGaPueY9LR+Zx11n9YN2n8OREOPP+wLuriYSBmS1yzuU3NM7zk8oiLVn7VklM7J/N3wo2cqC8CnJH\nQccBgZPLUfTHmsQGFYKIxy4fnce+8ipe/KIIzGDk9bBjNax51+to0sKoEEQ8NqRrGwZ2yeDJz9YF\n7m907LmQ1kGXoEqzUyGIRIArRuWxtvgAH31TDHGJcPzVUPgOFK/2Opq0ICoEkQgwsX82Wa0SeeLT\ndYEd+VcG3itBRwnSjFQIIhEgIc7H1OG5fPh1MWuK90NqJgy4EJY8B/v1Cn1pHioEkQhx0fCuJPh9\n/OWzdYEdo26CqnKYr5foSPNQIYhEiKxWiZw5MJsXFhWxt6wSsnoHXqi28FEo3+d1PGkBVAgiEeSK\nUd0orahm9sLgTYJPuBXKSmDRX7wNJi2CCkEkgvTPSSc/tw2zPl9PdY2DnHzIGwOfPwRVFV7Hkxin\nQhCJMFeM7saGXaW8t2p7YMcJt8C+zbBstrfBJOapEEQizGnHdiA7PYknP/s2sKPHqdCxP3xyP9TU\neBtOYpoKQSTCxPt9XDIyl08Ld7J6677A7SxG3wI7v4HVc72OJzFMhSASgaYc35XEON//HSX0Owcy\ncuGTP+imd9JkVAgiEahNagLnDenMi19sYvu+MvDHweibYFMBrP/U63gSo1QIIhFq2ok9qKqu4bGP\ng0cJgy6G1KzAuQSRJqBCEIlQ3TJTOWNAJ56et549pRUQnwzDvw+Fb8PWZV7HkxikQhCJYNeP7cGB\nimqe/OftLI6/GhLS4NMHPM0lsUmFIBLBjunYmu/07cATn65jf3kVJLcJvLXm8hdhR6HX8STGqBBE\nItz1Y3tQcrCSZ+evD+wYdRPEJcEHd3sbTGKOCkEkwg3u2obRPdvx6MffUlZZDWntYcS1gaMEnUuQ\nMFIhiESB60/uSfG+cv62qCiwY9SNkJQO7/3S22ASU1QIIlFgZI92DO6awSMfrKGyuiZwLmH0zfD1\nm7BhvtfxJEaEVAhmNt7MVptZoZndcYjHzcxmBB9famZDGpprZheY2QozqzGz/PAsRyQ2mRk3jO3J\npj0HmfPl5sDO4dMhtT28+z969bKERYOFYGZ+4CFgAtAPmGJm/eoMmwD0Cn5MAx4OYe5y4Dzgo6Nf\nhkjsO+WY9hzTsRV//KCQmhoHCalw4g9h/Sew9n2v40kMCOUIYRhQ6Jxb65yrAJ4HJtUZMwmY5QLm\nARlmll3fXOfcSufc6rCtRCTGmRnXj+3JmuIDvLlia2Dn0MsgvauOEiQsQimEzsDGWttFwX2hjAll\nroiEaGL/bLplpvLQ+4U45yAuEU6+AzYvhlWveR1PolzEn1Q2s2lmVmBmBcXFxV7HEfGU32dcd3IP\nVmzey6tLtwR2DvgeZPYOXHFUU+1tQIlqoRTCJqBLre2c4L5QxoQyt17OuZnOuXznXH5WVlZjporE\npPOG5NAvuzW/fmNV4HUJ/jgYeycUr4Jlf/M6nkSxUAphIdDLzLqZWQIwGZhTZ8wc4NLg1UYjgBLn\n3JYQ54pII/h9xk/P7MumPQd57JPgnVD7ng3ZA+H9u/Xey3LEGiwE51wVcAPwFrASmO2cW2Fm081s\nenDYXGAtUAg8ClxX31wAMzvXzIqAkcDrZvZWWFcmEsNG9cjktH4d+OP7hWzfWwY+H5xyF+xZD/P+\n6HU8iVLmoujKhPz8fFdQUOB1DJGIsG7HAcb94UPOHdyZ+84fGNj5/MVQ+C5cPw/a5HmaTyKHmS1y\nzjX4eq+IP6ksIoeWl5nKZSPz+NuiIpZvKgnsnPBr8Pnh9R/oMlRpNBWCSBS78dReZCTH88vXvwpc\nhpqeA6f8Pyh8B1a85HU8iTIqBJEolp4cz23jejNv7S7+8dW2wM5h10CnwfDGHXBwt7cBJaqoEESi\n3JRhXenVPo27566kvKo68JTRWQ9A6U545+dex5MookIQiXJxfh93ntGX9TtLmfVZ8E10sgcG3jNh\n0ZOwYZ6n+SR6qBBEYsDJfdpzUu8sZrz3DTv2lwd3/gTSu8CrN+u1CRISFYJIjPh/Z/alvKqGH/5t\nSeBuqIlpMPG3gVcwfzbD63gSBVQIIjGiZ/tW3DmxL++vLubxT4OvYO4zHvpNgg/vg51rvA0oEU+F\nIBJDLh2Zy7h+Hfj1m6tYVhR8bcL4X0N8Esy+DCoOeBtQIpoKQSSGmBm/OX8AmWmJ3PjcF+wvr4LW\n2fDdx2DbcnjlOr1gTQ5LhSASYzJSEnhg8mA27CrlrleWB3b2Ggfj/hu+egU+/q23ASViqRBEYtCw\nbm25+dTevLR4Ey99URTYOeqmwHsnvPdLWDXX24ASkVQIIjHqhlN6MrxbW376ynLWFu8Hs8AL1joN\nhpeuge0r23ioAAAInklEQVQrvY4oEUaFIBKj/D7j/smDSIjzceNziwOvYo5PhsnPQkIqPDcFSnd5\nHVMiiApBJIZlpyfzm/MHsmLzXm6bvYTK6hpo3Qm+9zTs3QQvXAHVVV7HlAihQhCJceP6deDOiX15\nfekWrn16UeBtN7sMgzP/AGs/gNduVikIoEIQaRGuObE7vzjnON5ZuZ1rZhVQWlEFg6fCiT+CxU/D\nsxdCWYnXMcVjKgSRFuKSEbn89oKBfFq4g8seX8Desko45U44awZ8+yE8djrsXu91TPGQCkGkBTl/\naA7/O2UIizfsYeqf57P7QAUMvQymvgh7N8OfT4WNC72OKR5RIYi0MGcMyOZPlwxl1dZ9TJ45j+37\nyqD7yXD1O4Grj548A5a/6HVM8YAKQaQFOrVvB564/Hg27Cpl/P0f8+z8DVS36wVXvxd4ncILV8K7\nv4CKUq+jSjNSIYi0UKN7ZvLy9aPo2T6N/3p5GWc/+AkLiw0umwMDpwRucTFjEMyfCVXlXseVZmAu\nim50lZ+f7woKCryOIRJTnHO8tnQL98xdyeaSMs4e2ImfTDyG7D2LA7e5WP9p4I12TvpxoCj8cV5H\nlkYys0XOufyGxoV0hGBm481stZkVmtkdh3jczGxG8PGlZjakoblm1tbM3jazb4L/tgl1cSISPmbG\nWQM78e4PTuamU3vx1oqtnPLbD/n1yrYs+86zuKkvQWomzLkB/jgclvxVl6jGqAaPEMzMD3wNjAOK\ngIXAFOfcV7XGTARuBCYCw4EHnHPD65trZvcBu5xz9waLoo1z7sf1ZdERgkjT27irlHveWMkby7fi\nHHRoncgpfdozufVS+n/9IL7ilWB+6DwUeoyF7mMhJx/88V5Hl8MI9QghlEIYCfzcOXd6cPsnAM65\ne2qN+RPwgXPuueD2auBkIO9wc/85xjm3xcyyg/P71JdFhSDSfHbuL+f91cW8u3IbH31dzIGKapLj\nYWr2Zsb4lnFM6SKy9n6FUYNLSIO80Vi7XpCeE/ho3Tnwb0om+HS60kuhFkIoTwZ2BjbW2i4icBTQ\n0JjODczt4JzbEvx8K9AhhCwi0kzapSVy/tAczh+aQ3lVNQu+3cW7K7fz8dpWPL+7K/vKT6c1+xnl\n+4oxVcsYsXoZne19kqj4t69TSTwHfalU+JKotEQqfIlUBj+vtjic+ajBhzPD4Q/+a4AFv0Lg33/+\n6erMaInajbud7sfV/dUbXhFxdsg558zskIcqZjYNmAbQtWvXZs0lIgGJcX7G9MpiTK+sf+0rOVhJ\n0e5SinafRNHugzy35yAHyqvwl+0i+eAWUsu20apiGxmV20iqLiXBlZPoykioriCxqpxEt4cEqjFq\n8DlHsBYwarDgr////Lfl2r13Z5N/j1AKYRPQpdZ2TnBfKGPi65m7zcyyaz1ltP1Q39w5NxOYCYGn\njELIKyLNID05nvTkdI7tlO51lBahUzN8j1Ce2FsI9DKzbmaWAEwG5tQZMwe4NHi10QigJPh0UH1z\n5wCXBT+/DPj7Ua5FRESOQoNHCM65KjO7AXgL8AOPO+dWmNn04OOPAHMJXGFUCJQCV9Q3N/il7wVm\nm9lVwHrgwrCuTEREGkUvTBMRiXFhfWGaiIjEPhWCiIgAKgQREQlSIYiICKBCEBGRoKi6ysjMiglc\nopoO1L7dYu3tw32eCewIQ4y63/toxh7u8UPtb2lrbsx2NK65sT/jutuRvOZw/Xddd1trPvI15zrn\nshoc5ZyLug9g5uG26/m8oCm+99GMPdzjh9rf0tbcmO1oXHNjf8bRtOZw/XetNTfNmuv7iNanjF6t\nZ/twnzfV9z6asYd7/FD7W9qaG7MdjWtu7M+47nYkrzlc/13X3daam1hUPWV0NMyswIXwwoxYojW3\nDFpzy9Aca47WI4QjMdPrAB7QmlsGrbllaPI1t5gjBBERqV9LOkIQEZF6qBBERARQIYiISJAKATCz\nk83sYzN7xMxO9jpPczGzVDMrMLMzvc7SHMysb/Bn/IKZXet1nuZgZueY2aNm9lczO83rPE3NzLqb\n2WNm9oLXWZpS8P/dvwR/theH6+tGfSGY2eNmtt3MltfZP97MVptZoZnd0cCXccB+IAkoaqqs4RKm\nNQP8GJjdNCnDKxxrds6tdM5NJ/BmTKObMm84hGnNrzjnrgGmA99ryrxHK0zrXeucu6ppkzaNRq7/\nPOCF4M/27LBliParjMzsRAK/zGc5544L7vMDXwPjCPyCXwhMIfCubffU+RJXAjucczVm1gH4vXMu\nbI3bFMK05oFAOwIluMM591rzpD8y4Vizc267mZ0NXAs85Zx7trnyH4lwrTk473fAM865L5opfqOF\neb0vOOfOb67s4dDI9U8C3nDOfWlmzzrnLgpHhgbfQjPSOec+MrO8OruHAYXOubUAZvY8MMk5dw9Q\n39Mju4HEpsgZTuFYc/CpsVSgH3DQzOY652qaMvfRCNfP2Tk3B5hjZq8DEV0IYfo5G4G3q30jkssA\nwv7/ctRpzPoJlEMO8CVhfKYn6gvhMDoDG2ttFwHDDzfYzM4DTgcygAebNlqTadSanXN3ApjZ5QSP\nkJo0XdNo7M/5ZAKH2okE3gc8GjVqzcCNwHeAdDPr6QLvgR5NGvszbgf8ChhsZj8JFkc0O9z6ZwAP\nmtkZhPH2FrFaCI3inHsJeMnrHF5wzj3pdYbm4pz7APjA4xjNyjk3g8AvjxbBObeTwPmSmOacOwBc\nEe6vG/UnlQ9jE9Cl1nZOcF8s05q15ljU0tZbV7OuP1YLYSHQy8y6mVkCMBmY43GmpqY1a82xqKWt\nt65mXX/UF4KZPQd8DvQxsyIzu8o5VwXcALwFrARmO+dWeJkznLRmrZkYXHNLW29dkbD+qL/sVERE\nwiPqjxBERCQ8VAgiIgKoEEREJEiFICIigApBRESCVAgiIgKoEEREJEiFICIigApBRESC/j+vrOWe\nc/DBOgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11e4efcf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot correlation functions and normalize them to the same msd\n",
    "plt.semilogx(time,acf_d*acf[0]/acf_d[0])\n",
    "plt.semilogx(time,acf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "delta_t = 0.00001\n",
    "def simulate_membrane(N):\n",
    "    ts_all = None\n",
    "    for q in q_list:\n",
    "        ts = langevin.time_series2(1/(kappa*q**4+sigma*q**2+ga),gamma(q),delta_t,N)\n",
    "    #    print(1/(kappa*q**4+sigma*q**2+ga),gamma(q))\n",
    "        if ts_all is None:\n",
    "            ts_all = ts\n",
    "        else:\n",
    "            ts_all = ts_all + ts\n",
    "    return ts_all\n",
    "\n",
    "ts_all=simulate_membrane(1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11e484940>]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXeYFtX1x79nd2GBZQGBpZdFWXp3QQRBBVQQFGPFikZD\nNNbEXxLsxoqJMTGiIKLGJBo1aoIRgoUuvYp0Flh6WXpZypb7+2PfWWZnZ+adcqe875zP8+yzb7nv\n3PPOO3PPveecew4JIcAwDMNEj5SgBWAYhmGCgRUAwzBMRGEFwDAME1FYATAMw0QUVgAMwzARhRUA\nwzBMRGEFwDAME1FYATAMw0QUVgAMwzARJS1oAcyoX7++yM7ODloMhmGYhGHp0qX7hRBZVtqGWgFk\nZ2djyZIlQYvBMAyTMBDRVqtt2QTEMAwTUVgBMAzDRBRWAAzDMBGFFQDDMExEYQXAMAwTUVgBMAzD\nRBRWAAzDMBGFFQDDRIh5m/Zjc8HxoMVgQkKoN4IxDCOXW95ZCADIHzM0YEmYMMArAIZhmIjCCoBh\nGCaiSFEARDSYiNYTUR4RjTZocwkRrSCi1UQ0S0a/DMMwjHNc+wCIKBXAmwAuA7ADwGIi+lIIsUbV\npg6AtwAMFkJsI6IGbvtlGIZh3CFjBdALQJ4QYrMQ4gyAjwEM17S5BcAXQohtACCE2CehX4ZhkpSx\n0zcie/TkoMVIemQogKYAtque74i9pqYNgHOIaCYRLSWiO4wORkSjiGgJES0pKCiQIB7DMInGq99s\nCFqESOCXEzgNwPkAhgK4AsBTRNRGr6EQYoIQIlcIkZuVZammAcMwDOMAGQpgJ4DmqufNYq+p2QHg\nayHECSHEfgCzAXSV0DfDOGbn4ZPIHj0ZP2w/HLQoDBMIMhTAYgA5RNSKiKoCGAHgS02bSQAuIqI0\nIqoB4AIAayX0zTC2KCopReGZYgDA7A1lJsZ/LtoWpEiBUFRSihe+WoMjhUVBi8IEiGsFIIQoBvAA\ngK9RNqh/KoRYTUT3EtG9sTZrAUwFsBLAIgAThRCr3PbNMHa56/3F6PD010GLETiTV+7GxO+34KUp\nPA+LMlJSQQghpgCYonltvOb5HwD8QUZ/DOOU7/P2By1CKCguFQCAotLSgCVhgoR3AjORR4igJfAf\n4cOXnrdpP7YdKNR97/jpYhw5Gbz5ac+RU1i960jQYgQGJ4NjIgsFLUAIIA/PglniufOf/xani0sD\nT0rX++VpAKKbHI9XAEzS8ff5+Xj44+VBi8GYcLqYTU9hgBUAk3Q8NWk1Jq3YFbQYoSaCVi9GB1YA\nTOQRER4Oie1gkSZSCmDZtkMYO31j0GIwDGMRP5zVUSZSCuDat+ZxjpGQkz16MibO2exrn146Qhkm\nzERKATCJwQuTeXMSw/gBKwAA36zeg5NnSoIWgwkIPR/AkcIiLN16MABpfIItKwxYAWDVziMY9fel\neHoSZ6aIGmYO0GFj5+C6cfNRXGIervif5Ttx1RvfS5bMHruPnMRbM/Mc2cvZ+BVtIq8Ajp4q2424\n/ZD+jsVkZMv+E0GLoIufDr9bJy7AtLXGdYm2HzwJIP5E+ZFPVuDHne53kgoh8MrUdVjl4Fj3/n0p\nfj91PTbb+F0TJfJJxiWxdvdR7D16yv2BkpDIK4CoMXP9Plz66kxMWqHN2G2fySt3I3v05HIlmkjM\nzTuAb9bsDVqMcs6UlGLczE249q15tj97/HRZdlNHK4AkWAIcLjyDRz/9ASdi50HLkNfnoM+Y6T5L\nlRiwAogY6/ccAwCs3nXU9bHGzcoDAGzdn9irJ7NxM1mjEA+eKFPa//1hd8CSuOcv0/Lw+bIdpmm9\nS0qT9Id0CSuAOIybuQmdn+X0wWbIMick62DrNU5OW96+4wCAk0XBBz8cO1WEzs98jbmcqdV3Iq0A\nth8sxB9j+wKMBp9Xpq7DsVP6S8tERoa9PVni562YQZZuPYhxMzd5L4wrEvP3WLv7GI6dLsafv+M9\nOn4TaQVw7z+WYunWQ0GL4RlTV+3BNW/OrTDYe2Hz9XLmvmjLwUCX78rq5rpx8/HK1HWByWEHIQRu\nHD8f35n4OBLFCWyFZPBjBEUkFcDsDQVYtfMIziR5RsIHPlqGFdsPlxf/sMKRk0U4ZtGp6/WNN3/T\nAdz49ny8NSPP036SzfR0urgUi/IP4v6PlgUtimuS7KcJHZGsB3DHe4sAADkNaro6zsLNBzBn435c\n3a0J2jTMlCGabxgNel1/9w0Ae/nRnd6k6/ccQ06DmkhJId3j7DlaFoqZV3DcYQ/+smrnEZwsKkHP\n7LqOj+FoZu7gI7sPJ09YZGFsE2eyKXI/iOQKQA8ns9mbJizA2Bl5uPxPs+ULJBH1jRG03X7j3mP4\n2/x8LN92CFf8eTbesZD3J0iJ7Qwqw974HjeMn++oHxm/i51reO+x8CkAvXMdz1e1dOvB8ugfOyvd\nIDh2qgjzQuboZgUQMJNX7saQ1+f4tgkq/0DZZqGC46d96U/LlX+Zg6cnrcaOQ2Wz+5UmG5+SeUb3\n+L9/RPboyeXP/bbJe6VUJ63Yif02ry03psRlWw+XP/7b/HznB/KBBz5ajlsmLsSBgO49PSKtANQX\n3oLNweR9efjj5Vi7+2jc2curX6/HRa+438zy4cKy2dKUH93Hfzu5b4tKnGxWit9T3r7jjp3FWw3q\n1gLAwD/OcnTMeHy0UD9m3c1KQKswzc7GkZPyI9sKjp3Gwx+vwD0fLLH8menrzDfj2flFDxda35C4\nfNuh8g10frFhb9kenHjV0P734268Mc2ftPVSFAARDSai9USUR0SjTdr1JKJiIrpeRr+yCUORaiPG\nzsjDjkMn4xaw3rj3WHliOzOlEm92vf1gYXmseDxkrF6KSkoxe2OBo+NvKjiOQa/NwnmPT8HOwydt\n970o/yCKDHL+WD2eWsb8OCkZ1PmFThWVYPvBQoyYsKBSuzPFpfjr3C2eREFZnaW/P3eL5aLpyjnU\npl0Y8MeZhp/5/dT1lo5theJSa0EdhWeK8ZO35uG+fyyV1rdM7vtwGf74rT8hsa4VABGlAngTwBAA\nHQDcTEQdDNq9AuAbt31aZeWOw4Y3ti4JYHIY+hfjxGOni0tw2Z9m4xcfxr+w4604+v1+Bga9Fmf2\nKyMMKCbGn7/bgLveX1zxLaF0Y96PesD5k8GNs/1gId6aaRxNVCrR3nTJqzNN339pytlw0p/9bQn+\n9N0GLN92uFK7t2dtwrP/XYNPFm+31K+sqKzl2w4he/RkbNh7DL/77xrTa04P7ancXOA895Sdn8Xq\n6rKouKzdD9srn/OoIWMF0AtAnhBisxDiDICPAQzXafcggM8BGGfgkkjevuO4euxcvDTF/9zyB46f\nRuvHp2Bxvr9mJWWm6Lc5S3vbHSkswsod5jeXMmgrMpslqIs7rlm47+98f1Hc2WZRSSlWWBgUdhwq\ndGU+mKNa6czZWNEpqPYFKCtSoxw3XjF5ZZl5cNZ6/RWZEYoC2mMz8ZrZCseub+SlKWtDn/ht8srd\nvv+mRshQAE0BqKcoO2KvlUNETQH8BMA4Cf1Z4uCJMwBQnl2xzZP/q9Smkr1V0gxqcf4hFJcKvDPb\n38pWTiaxt05cgNvfXeioP6PTNeKdBbh67FzTzypFXxZsOWDYRuaCrNBCvYeXp6zDNW/OxcaYrdaI\ni16Zges0Sdu8dFg7dRAnyv4oJRW7jCieCbM344KXpsWdgASBknbjxSlr8cS/fwxYmjL8cgL/GcBv\nhRBx7TFENIqIlhDRkoICezMQNVq7cTJt+oo3+9SaAswGkLl5ByrNQp3wxrSN5Ttl1+62nmhOxiCl\n/nYzbc5a1ayK2boPxCYPZqyPoyTMkL2Bzivd41dk0oa9Zb4mPTOcU8X6ts+TLyuondS7joRjlSJD\nAewE0Fz1vFnsNTW5AD4monwA1wN4i4iu0TuYEGKCECJXCJGblZXlWjg7URXybkxvb5xOz+gnp7Pb\n6ymJicD++O0GV7ly7Ea/FJeUolRnxmjk3NwdkhsOqPxd1XUJikpEpTBBq+dG28rpVWh2H3yyeBu+\nlzBh0MPvsF833Z0qKsF/lu9M+KL1MhTAYgA5RNSKiKoCGAHgS3UDIUQrIUS2ECIbwGcAfiGE+I+E\nvg1RfpZF+QcxddUe3TZuZnFWsKJQ9C6fbQcKsWiLczu+Ubcz1lV0vzzw0XLHfahxew/sO3oKWw9W\n9gEoN9fy7Ycr2XVbP/E/PPqvH6T0D5TZofdIVhJHThbhsS9+jFtuVBt9dv4L31V4rszEzxSXVgpq\nOHmmRHqBH7Pz+dvPf8RtFkyGEy1s8EtkXpm6Do98ssLV6vnd77fg/blbJEplH9cKQAhRDOABAF8D\nWAvgUyHEaiK6l4judXt8GdwrOdxr6daDlWaf+46dQnFJKV6ashb7j8c3IWhR33T9/zADN75tf0ep\n0WxEefmuv1aMsvlurbuCKLJWTL1emoZVO43NRlv2n9CNx//3cvdFbRQOFxZh28Gy/QBmTkT15i01\nemf+zRl5+Oeibfhw4VZXsilO/bZP/Q99xkxHSakonyColcf6PccchcECZY7mBZuN/TFWUK9UFB9P\nPNbtOTsJc2JycnMNurl8lWtEyRQshDCMKtp1+KTudfP8V2vwu/+ucSGFe6T4AIQQU4QQbYQQ5wkh\nXoy9Nl4IMV6n7Z1CiM9k9GuGE5s/Afhq5S7sMrmJ5m86gOvGzce4WRXNHb1enIZZGwowYfZmPGWj\nvrAfjjpl161bhBAoKRU4XVyiiWJwPgW3sskLqOj3uOZNcwezE9Qpv7W2aCdlGgHgUMyXoPUp2B20\npq/bh8X5ByFE2War8bM24ca351dKK3Dfh8vQ12Hlq0c//QEjJizA3qOnLMsnhDC9V9TEq61s3Iez\n9wllyrrVY/oK2yl/m5+P7NGTcaa4Ysfvzc3H8Dfn6tY0cKtYvSRpdwLbcUSqeeCj5bhunHFZvt1H\nYgnKdDZJKeFsdkwSflgQrdj6TxeX4HSxebtxszbhvMenYMCrs9Dxma/lOHBNTpbRO9pQzXgzRysb\nqdRRGVqb+zCHRd//tXQHgLKlvlv2HzvrF1Cuvd1HTklz1P4Qi5qxEi2l8O73W9BnzHSs2xP/Xntd\n0s7Wk2dKLCkTZbet2b141EGdj9di+0yOn65otlsfOwcP/XN5QlUfS1oF4AanDkPtz27HsSkz4sLq\nrFpN52e+QY/nvjVt8/GismhfrZnBjQ3eiaxq9h07FTd8cPys+M7pJR7WhXD6DY1OjfK60bfW21QW\nD+WaX73rbJr0eBu4FNPUNpNUGgpW/BTbD1ZeTWjvi/ZPT8XdFlJNfOtRvWclkudoLJWGIp9yrx84\ncQbT11nY6mRyyV7dtQla1c9wJ6hFklYBBKGD9x2rGL2xYMsB9B0z3TBqB5BrAnLznc+UlOKEwexv\nbt5+3dXByh1lphGndmcgzvfXfKGjp4oqpRXo9eI0rIgz4P3ha3npBowwS0Ogzv2yce+xCnZvJyiD\nze7DJ7Fxb+WVqJsNiIcLi7A4v0wZfmxxB/I3FgZbK9em1VQsszY4D/UFoHvRqc18M9bvQ1FJKU6c\nLsYni7fprlKVpIp6aFcAYS5YE8l6AEaoL8Cjp4pQq1oVW59/6j8Vbf+HC4viJqjyQlFVCgd00cnq\nXUdw68SFuL13y0rvKTPv2RvO2j3/4mESqy7P+pNFxO5OVuDs7tl4yNiZrgwoRvlijH7vr1buct23\nHp8t3YFXb+jqWbW5LftPoFqVFDSuXb3S+4761Dk/w974HvljhmLepv246/3F+MUl52H/8dP4dMkO\ntKibgZb1aqBJnbP9m5tV3d3Vfk5ek3cF4OAs7j16dgZfYpBXRDmum1mv2XHVZI+ebBh1ombinM3Y\nfrDQkzhqRYFtMinKol6mv2YziZXZDWx1Z6iRfXnczE2Wzp8WJ+YDbYbHlyUM9F/+oD9gxxvz/hSr\nrasNhLAa9mt1UPVrZnvpqzNx4cvuM+Fa4UAsgm/rwUIUxFb0/1qyHX3GTK/gdFcuTcf3nMG5W7Xz\nCP77wy7fdnEnrQJwgjaiQW/p98mSsmWxmzh9mew7dgovTF6Lke8twvxNZRfoMc1OYc93dJoc/r8/\n7MIni/VTHwMwDJl9c0YeHne5Xd7PGr5jp1dMNGe0E9WOz0M9IfGTMJksgnSnbow525dtKzOJ3TKx\n8v6HSSsqhyJbUgoGbZwGHDglaRWAk0Fv3qaz4VpE+j9kWAZ+BUXGHYdOSknpoEVZAajPTSUZTD7/\n4D+X47ef2x/I/bDby8TqilDG2OpkgLaa/lnAox25Fo+pDLaeYvH8KeHT+SZO7rW7y/w5er9JSamw\nlI14Xt5+zNsUTKWwpFUAzjh7lQZdOtEqSqTGmZJSeBF9pk6qpWyU0pLo2+G1WP3l2z811TMZzM6p\n3Wtz//HTyNXsLjY+tjMF88WyHbjgpWmG71udkF371rwKkyy/ry2llsNRiw5pM8U/9C9zkPNE5SSU\nWm6ZuBC3vOMsIaNbklYBOBvAz35m7Z6jiVAeAPtUtV2NZnlu7iErpSP9rqzkNVZDQk86yKVkdXBd\nvcs4tt7uAH3QQnK7Csd3kHtIZtp1Za8NUHFfwouTvd81+8aMMlOejNW0Eu2ldz4XxYvU8mn+mbQK\nwAk1qqaWPx4xYYG9YjIeor4hzDBzXjrdzaq3s1HLd2t9KfGQ0Exftzd2PVm7s81swVbDMwGU7yA2\nolGtahWer9tzFNsPxY/rBypadbyaqC9UrQbemSMpb46JrMmUNdgKrABUqBUAUHEjjLIU1bZxi5Xd\nw24jIE6cKXbsXJLliLzp7flSs48mGj/96xLc/+EySysqmcyJEzPftXntCs//sWCb5dq6dqOllm49\nhAmz7WWMlWkC8srO/vXqPRWUcqUNoSG2JietAsisZn+LQ42qFT+jtlu2emwKAKBb8zqO5BFC4O1Z\nm3C48Oxy/FWVo9PKTNspQdkX1SzcctBxeo5k4Zs1ewMpQ2hmf5fl64pXQ+HbNXtx3bh5Fcph+k35\nfSB5QP753ysmm1ynuc7D7CJLWgVwXlZN259JTzM/HW/P2mTbnqrwfd5+vPy/dRitioh5R5Uy956/\nxd/erkfNdN7LxxizdNshPPLxCsP3/Sr6YrVeL+DtgPnMpFV47Rv9CLNXJUWe/XPx9vJEgID1CCw1\nfi0aknb0kJFaVnshvvw/57OXZyatBiB/GZppc7eyn3T9nT87d6PAcocrh7l55pkoEyhvmRQ+mG+c\nmnvsjDzD9+xQcOw0uj9/Nq+W1dTYQZC0CsAJXs48NivhZaoMhFqFc+RkEaqm2luUrdnlzLnrB1Zz\nuzDxmeBRicOw/0bbDUKPFcJsXkkEWAGo0CZzc5NUywldf/cNmtapnO/EjGcNCkrolUtkGC37jwWz\n29gM9er91W/0U4tkj56M7397qV8iJS1J6wNwgjZtcBDVemTlGLrgZeNNOQyjsHn/CUv5ko6dCt9K\n4bNYvYVkxG2adKskrwJIgAnwqSLvYo4LQjizYxKXzhIzsTpJ0KfHn7/zLvNsVGATUAAIIfCmJIdT\nIjFphTfpiBl/mGGl0IkEomrXz9vnrk6EE5J3BRDizRcFx08b2jaTmb/Oyw9aBMYFd/11cfnjsOR/\n0pYHTWReDCBaKHkVQDiuT112HXZWcpJhwoJXJRftMvlHa4V4Eo2EqgdARIOJaD0R5RHRaJ33byWi\nlUT0IxHNI6KuMvpNVO77x9L4jeJgpcYqw3iFlwkAlVKjjPe4VgBElArgTQBDAHQAcDMRddA02wLg\nYiFEZwDPA5jgtt9EptCg9q4d/r28ciEKhkkGLBVVZ6QgYwXQC0CeEGKzEOIMgI8BDFc3EELME0Io\nOXYXAGgmod+EJczJoRgmaIzqTjDykaEAmgJQ56fdEXvNiLsBxK+S4JIQuwAsF5tgGCaa+DVJ9DUM\nlIguRZkCuMikzSgAowCgRYsWPknmL7xJl2EYLX5t/lIjYwWwE0Bz1fNmsdcqQERdAEwEMFwIYZih\nSggxQQiRK4TIzcrKkiAewzBM+AkitFaGAlgMIIeIWhFRVQAjAHypbkBELQB8AeB2IUT0AuA9gN0I\nTJCs3+v/pqVkZ2sAvg/XCkAIUQzgAQBfA1gL4FMhxGoiupeI7o01expAPQBvEdEKInKW/J4p5/Vp\nvA2eCY4Dx53VxWCMUVcglFWoJx5SfABCiCkApmheG696fA+Ae2T0ZRWrZe0YhrFP1GrnJitJuxP4\n/o+WBS0CwyQtJSFJBcG4I2kVQBQIY4peJhpMXpmcKRjCgl8+FlYACYzMFL0Mw0QPVgAMwzARhRUA\nwzBMRGEFwDAME1FYATAMw0QUVgAMwzARhRUAwzBMRElKBfDe91uCFoFhGCb0JKUCeO6rNUGLwDAM\nE3qSUgEwDMMw8WEFwDAME1FYATAMw0QUVgAMwzARhRUAwzBMRGEFwDAME1FYATAMw0QUVgAMwzAR\nhRUAwzBMRGEFwDAME1FYATAMw0QUKQqAiAYT0XoiyiOi0TrvExH9Jfb+SiLqIaNfhmEYxjmuFQAR\npQJ4E8AQAB0A3ExEHTTNhgDIif2NAjDObb8MwzCMO2SsAHoByBNCbBZCnAHwMYDhmjbDAfxNlLEA\nQB0iaiyhb4ZhGMYhMhRAUwDbVc93xF6z2wYAQESjiGgJES0pKCiQIB7DMAyjR+icwEKICUKIXCFE\nblZWVtDiMAzDJC0yFMBOAM1Vz5vFXrPbhmEYhvERGQpgMYAcImpFRFUBjADwpabNlwDuiEUD9QZw\nRAixW0LfDMMwjEPS3B5ACFFMRA8A+BpAKoD3hBCrieje2PvjAUwBcCWAPACFAO5y2y/DMAzjDtcK\nAACEEFNQNsirXxuveiwA3C+jL4ZhGEYOoXMCMwzDMP7ACoBhGCaisAJgGIaJKKwAGIZhIgorAIZh\nmIjCCoBhGCaisAJgGIaJKKwAGIZhIgorAIZhmIjCCoBhGCaisAJgGIaJKKwAGIZhIgorAIZhmIjC\nCoBhGCaisAJgGIaJKKwAGIZhIgorAIZhmIjCCoBhGCaisAJgGIaJKKwAGIZhIgorAIZhGB94bnjH\noEWohCsFQER1iehbItoY+3+OTpvmRDSDiNYQ0WoiethNnwzDVKZjk1pBi8AkIG5XAKMBTBNC5ACY\nFnuupRjAo0KIDgB6A7ifiDq47JdhEpoOjeUO2LKPx0QDtwpgOIAPYo8/AHCNtoEQYrcQYlns8TEA\nawE0ddkvw1Ti+vObBS1CYOw/fjpoEZg4UNAC6OBWATQUQuyOPd4DoKFZYyLKBtAdwEKX/TJMJVLC\neIcZkJGeKvV4qSkp+Pcv+kg9JiOXzGpVghahEmnxGhDRdwAa6bz1hPqJEEIQkTA5Tk0AnwN4RAhx\n1KTdKACjAKBFixbxxGOYciiUcyx90tPkKoDTxSXo3qKSC86Uc+tnYPP+E1LlYIzJykwPWoRKxF0B\nCCEGCSE66fxNArCXiBoDQOz/Pr1jEFEVlA3+HwohvojT3wQhRK4QIjcrK8v+N2J8oVd23cD6bt2g\npu7rlDjjv3TmbNxv+zPX53prMsuuV8PT44eNJrWrmb4fxsvTrQnoSwAjY49HApikbUBEBOBdAGuF\nEK+57I8JAb++oi1SA7S3tKqfEVjfsgiDskrxWIiZv77U0+OHjRrp5gaV2jXCZwJyqwDGALiMiDYC\nGBR7DiJqQkRTYm36ArgdwAAiWhH7u9Jlv0yA3H9p66BF0CUMg6pVhKGxlHFCIvz2NeMoiCBwJZEQ\n4gCAgTqv7wJwZezx9wjn6ic0dGxSC6t3GbpFLJFZLQ3HThVLkihRSZzL7LbeLfF9nn2zjUzCoISy\nMtNRcMx9BFP1KqkoPFMiQSLvOCejatAiVIJ3AocAGeaU81vacwAmI5d1aBC0CJYZ3EkvrqIiaR6Z\n2c7NCo8JbUTP5lKP1+e8elKPZwcRBo1qE1YASUKphWtP5oCS4tGVU62K8wMPaGcahZxwND2nuifH\n/dVlbTw9vh1kq7iXr+0s+YjyCOP6lBWAj3i5Xd/K7MNrp58fZKan6d5IP+lecW9h/ZrhW26HhWFd\nmiB/zFBUryI3FNUJp4pLXR/j7otalT9OwEl4oLAC8JGqad6d7tIku/J/fUVb3dcXPD4QVr7pqzd0\nlStQACTZT6rLidPu/VZdm9cJxezaKDxZIYw/JyuAJKHEig0oDHeJRaoZzE69VKKyeW54RzStE7yZ\nxQl+hfnKcNyqJQ1ykP3jjd0C7N0ZiXM3hYTpj17s+LNezuhuuaCllOOkWxxgvdp1a+W4iaLH7rgw\nG3NHD3D8eTsWu57ZiRcEcGefbFzUur7r46jPU5CO2DCGecaDFYBNzs0yX+Zp6dqsdvljt5dm73ON\nd98O7dw47udlDpxeuBPaNcosf6y9kT+/rw9euKYTqqQm9iX79LAO+O8DF1lqa3aKJz9U8RiXdbDv\nADcbLGX9vBlVjf0Mz17dEUUl7n0AnZvWjt8IQN/W9iOEwuAn8ZLEvptCxt/v7oVvf9m/wmuTLNzs\nK3ccsXT8N27ugRsMMl5auWH/73J9u3qF43g8vTazLNzUszlevaErWjeoWW4CuuH8Zpg3egDOb3kO\nbuttvMpx42C/uZfcUEQz6mZURedm1gYsM2p5nFjMyJ7drXkdW8fpESc8ed2eY7aOp0eLujXw+ND2\nAIDqJgpnbt6BSq/1ahVcSpMwECkFkOtxrHy/nCzkNMw0buByeZqVmY7nr+lkuhIwIn/MUPys/7lx\n23mdUG3h44NM+gaGdmmM7351cXnEUloqoYnGjv7UsMrlJNz4Bn4ZC4sMG2SijbVpBepmyE00ZpRa\nu1Et83w3Wsbe3MPUOXrkZJGt4+lBRLj1gpbIHzPUNNJNb6VwVRfzlXMSBM6ZEikF0NNnbT/r15cA\nQLkjsF5N9zdptSqpuLZ75ZtT1oXq9QVfT8JuyOZ1zyYZszsj1cPPLKLLtx2SchztCsCKCVAGwqYh\ns3aNKhjc0XjTm182+waZ6fjTTZUjw9S9m8mZrERKAbjlmavsFTJrWa9sx+WsX1+C8bf18PQCM5st\nhgkzMZ1waiGXAAAW8UlEQVR8ByOTkp1j+XnqPpi/1XLbeGKpZ9Zmpo+gMTu/DeNk0JTFuyN7onUD\nk9U5gPG3ny+tv4cGtEa/HPcObq+JlAKQPdmwGuKXlpqCwZ0ao0MC1G2Nd45+f30XAOFROEqOda00\nejPLd+7ILX98V9/s8sfqz74+ohteuS4ku0njnOLHr2zn6vB+xcvUl7DytYr6Z9crkzmofcMK5sKj\nEkxQevzq8rZoa2YODgmRUgA7DhVKPd6FNvOOdLIYrRAkZn6CDo1r4cZcdw5TM8VRIZzP4vD0ynVd\ncFvvFrghjlyj+p9b/ntlVE3FI4PO2v3VMg3v1hQ39XReiOhcyamqv3rQWsSQbGQqeKuhxbLR+woT\nR+ZiwwtDyp9v2W8+Jsic5qSGZNKkJlIK4KuVu+M30uHOPtlyBQkx3ZvXwWND9GeWUx7uV/5YfSnL\nSrugv/rQv2mUmOta1arghWs6V9o4pjeA1aiSir6t6+HNW3tU2IEq87aUWfUpo2oaMqsZx5bXrm4e\nCfTLQebObSffW5tyIwjM/B0Na9k7/2bnF5C7SspIT8PbOmamIBM5RkoBOMXIURU+fS4BB1/q4YE5\nFZ4/f00nScIY85/7++DJoe2RYmPHakoK4cN7euOStg1wWpWDJoQTMwDAiz8xP4/ntzQParjvkvNM\n37c7uI2/7XwM7uTO2axVkE7i7M02vakVvxWTr+zazPG4QscPeLmDPRyyiIwCIDIuUffEle0NP1fD\nxLnmx8Dx0IBwFF/pqoldV3/3VE1q0NtN4vXNsHM+WzfIxD394oe1lh9b81y9AclqFNDYW7pjzXNX\nWO7TLfFm+PGwGxqrdlrqnRErKaz1UP+ui5+oGAacE8cxK4sWBvd+zfTwVenyk8goAEB/t+SCxwaa\n2r2rpKYE5vDMHzMUv7q8LYZ3a4JOTVUOLY/F6XNe5eiFi9sY12e2GxoYBiokIbN4PlvWzUCNqmmm\nOewFgFpxzApW8To8VTtDzj9gv0D8bb1bxHWaa/dxqIlngtFCZN0/oW6mVaZ39slGu0aZaBlTDJe0\nLbu+n72qQyh2/zrZteyESCkA5YK/5YKzTr5GOmFo2s0uxiYguTfotT3O2lfVdsHXR3THVw/20/uI\nayrN7AF0bla7gixAZXOBF0OTzGPuO3rKer9U8b+am3udvVaUHbyXtDUvPDP5oX4Yf1sPy/17iZ0a\nEKUOsjKkpaQYJu5T6JdjPHloYzNSZsvLQy3vHTBr9uzVHTH1kf7lbZSV/p19W2Ht84PL27m5JpWE\nevdc1Ap/vaunaVvtZrlxt8kLSTUjMgpAiLNFU+JtRmpVPwPLn7oMQMVBQXsxXNW1iUQJy2aYCndc\nKCe5W1yMZlOam6dHC38dVW5DdrcfNI/uSE87O2iZ3eQv/aQThnZujNdHnM30aOr0FmUb1azYyod3\nk3v96NHQ5s5dBauLXrcbufQmYDL5/fVdMPUR48mTkkY93mTuyaHtTfMaKbxxc3eMvaU7AOD+Aa1x\ne++WePTythUmDX/7aa9Kn8vQJJLzOtWHQuKlr3OBYqqoXb0Kvnygr2Eekn5t6iOzWhoaZKbjyWEd\nsDT/YIX301IIxaVCeh4Ry5YmDy0uyvJa20WbRpm67SR3buUlKVRJPXtg5bvodUVEePPWHprXvJFJ\nD21qgxn/d4mtcObXR3TD9ePnW2prNtkxw+8EnN1tTEbihS2Xix7nC4/o1QL39DsX2aMnm7ZTTwpr\nVauiGxDRX2NObVFX3z/hB5FZAdTNqFp+oaYQoUuzOroXx9/v7oX7Lj4PaakpWPTEIFytM8tXLhon\nA8GNufo5VgD3JpAeLaylRbi0rfGSvGmdarH/5pvc1LZbtwOA4l/o1sx9WgdDNCe3To2qRm/ZP5gK\nmf6Qg4Vn0ExTtrFV/QxTk4qSJ0lxAOdme5v+RPbY38VCoryuzetg7XOD47azgrKCCapa3uf3XYgh\nnRsHFlGYlApg/QuDK6TLHXNtZ/z7F33KY8fNQr/65WRZdzI5kE1KiKRBx1/8om/crJhbXr4S791p\nbI9UtssP62puwvhp31a6r+ePGWr4GW0EjRKV9cTQ9lj3/GApWTKtUlWVVtruvS9rrNCGz2o5erLI\n9kpLUdyXmDjtFXI14ZR1LeZpeuW6zhjQ7qxJQ6bSs/p9ZaW+UCYvQQ3ASihvUGEUrhQAEdUlom+J\naGPsv+HajIhSiWg5EX3lpk8rpKelomOTs4PJiF4t0LJeBh4Y0BqPDWmH63oYz8L16NWqzCPfUWcn\nb50a9mx1atuzGWaz6naNjB1n6nQHehCRFPONOi+/1aiJGlUrWhzVx4jnSLRNnO+ojly168xXt171\nO+dhoedm1cS7I41/rw17y0yUZhFYbqhfM72Cws5pkFm+QtVeI+oSnTf1bGEq0+sjuuFpnYytVmhk\ncyOXGVYGVUV5+VQAzUSQYFSA2xXAaADThBA5AKbFnhvxMIC1LvtzRbUqqfh5zLyjR45B2tqhXRpj\nyZOD0FNnOW0nyiIeVsflLs3qYMmT+mmVm9SpXsk38enPL3QrWiXUsiqORm1EkRF/ubl73F2qrolz\nQ6kHOPsrgLMf0FaBenpYR1vHGtjeeBPQpoLjAOxuVnI+kDQ7p3r5BEV7TgxjBUTlUz28W1P89CL9\nFaIZCx8fiGbn+GsPVyKfrE6KZN7vRjTxKUEe4F4BDAfwQezxBwCu0WtERM0ADAUw0WV/npH34hBM\nfaS/4fvqhFbqyAcrpXi1/Kzf2Zujz3n1ym8uM9uumTxa1HUP6mZUteysNio2YxVtJIMRV3dtgocH\nmZs/XM+HHKxyZJj+nJix4s3wnUwOnSzyiM7uVNVOdrSrJPXxZeW4chqxZISVlZNVf55yzy9/+jL8\n3+UeTF5UAsS7N2TiVgE0FEIoCXb2ADCazvwZwG8AxI00JqJRRLSEiJYUFBS4FM86aakplgthXxAz\nCRERzrFpAgKAJ4Z2QOOYln/1hq7lCcTUMz2zIhrxUF/MdmYsfhUCt4Ndib771cWYN3oADhee0RxH\n3ndra2KCc8IHOmGBAJDpUyigAoFwUU595I8ZivY6mTT1EBBo0zCzPMFa1RCU7FSygN7WO35Sv3hh\noNpJQWa1Kq53aMejuY+roLi/FhF9R0SrdP6Gq9uJMhVZaa5CRMMA7BNCLLUikBBighAiVwiRm5Xl\nje3TLRNH5mLqI/2QmkJ4wGGqhmeu6ogGmemoXzO9wklrXrfMiScrDtjO5FE907TqEFSjN4syS0fw\n2o3dMKh9A7SSmEGzdYOaaFKnOmauN588qFdxdiNAFEXptTlA2ZB4aTvzjWdq3JiSzWs1aJ5r3q+a\nloLnhnfElIfdZS+VcUbfvTMXTw5tb8mcNKBdA2RlpuNn/fVNVn4VrOnXOpjaAXEVgBBikBCik87f\nJAB7iagxAMT+79M5RF8AVxNRPoCPAQwgon9I/A6+k5GehnaNymYZTouUD+7UCIueGFQ2QFoNRvaR\nBpnVsPDxgZayK5pFgSx4bKDhe52b1cbEkT11z+Gwzo3RsUkt/Ly/eUIzI7bF2QimoLXhhwllRVi3\nhgNl7NO1pB4f77gwO27RFQWjTVUyhtvGtatbzhNVv2Y6Fj8xqPx+NsLIPNhAUvbXoMqSul2vfQlg\nZOzxSACTtA2EEI8JIZoJIbIBjAAwXQhxm8t+Q4PMCYKTfPgy0fbZsFY1wxzmei/rDTpOVhIAcE5G\nVUx+qJ9hEi+7WJnkl1h06Cjn5LFYEkGv8rY0rVP23RvX8ccpaEtluIgk+/HZy7HYIIgh0ZgZK/vq\nlqDMr26nP2MAfEpEdwPYCuBGACCiJgAmCiGudHl8JmCMhkQ/6+jKoFMTfUelEMK2sk1JoQrhk16f\ni44GsuvRLmb/HtJZP3Onmc/KzpiubDoc2N66eUrBb9+GDNSmIPXVog1ttsrXj/RH4Zni+A09xpUC\nEEIcAFBpjS+E2AWg0uAvhJgJYKabPpMRveEniAHWboUzhYBCmC0z+9eXVlpJBJXhdZDBgLn4iUHo\n+eJ3ro/fqn4GNr44xNA0afZT2TknHZvUxqaXrpQ+c1WONrRzY0z+0VkBJ5nonRMZPiuzQAI/b6fg\nXfYJjgxHoCiPRFC95uIyUCsPO4PzT7pbDwPVNQGFdFEgy4zklAmqKlBGv4deJTGn14BTv9Smfcdt\ntffSbNGteR3bmyz9QlEAdquPhZHwesASBOUmkLFbk4h8m/l/du+FOCejKn716Q+W2mulUm6CZ65y\ntuPTa4jCszK5XFUFykyk87IysKlAPyf/N7/sj702UlwbcUUH46IuX6/eY/ieX7pdfX6U32/BYwOx\nZvcRZNeTW2/ZLWkp3syf/ZxHsQJwiXLBqmddb93aw3BXcViwmyRMO3BVSU0pt4PP3uDffg2rNKld\nHTsPnwxajEo4zfzYpmGm7dz5WhY9MRDnOIgoCpr0tBQMaBdc2UQtSmhzE5+c817CJiAPuLJzY+TY\nuFllT1T9NsWEZKJdgZ9fbL1cpB+ZIG+KZZ6VvYnMDg0yq5mah7T27vwxQzHKpFqeF6glCKtJsUFm\nNYy9pTvevt0871YiwArAJTKuUacZCbX1Vb3Ermx39c3GNT4UPDFCmS0bJc5TzNfqUNffDG6r21YG\n3WNRM+dlhXdlGNLxNpQM69LEcYhzmGATUAhoVKsath0stF3EW89x6BVms/ze59bFsC6N8Zsr2pW/\n9sxV9pKieUUtg237Naqm4fUR3XBBq3qVwjq94KaezdH73HrIthFBYmU3+Le/7I8SWc6OEGgAPR8A\n4x2sAELAuNt6YPbGggrFs+Vd/NYO9Mmo3o57SE9LxdhbwlED1w7DuzWN30gSRGRr8AeAemalJ2PY\nMTXGQ2/89ysVghlhMAXFOwtSU3b7+H2TWgHM+c2l2H3EfeSEGTJuj3o108tDMIO62OOtPkJwD/qO\nH6l/w4TZPgC/rstEPONbXk7c/a5JrQCa162B5gHW23SDxBpLUo4S/DzQX+aOHmC50I1sJt6Ri5U7\nDsctyykbs0mA1wuBET0rl2cNxerDSpswLFEcktQKwA/C+NOrZfroZxeYtr04pz5+2H4YDSTnYk90\n/B581bSsVwODOvgf9lhDUplFu+S9OCSUqciBYCY+fmYB4CigkCLjEnhkUE7c2PFHBrXBgscGxh3w\nwnl7GqMUzOnW3MNC8x7h94RyeCxaa+SF2YZtvJQpLTWlfBatFNRp0ygz7sz66WEdfIvESbTr3yq8\nAnCJVzMEV8e1cbempBAa+ViCzi9aN6iJKQ/1Q5uG4Q27DAvKzF+vLrPfVpjh3Zqia7M6lhzmP72o\nlaPSk04I3hjlDbwCCBnJOtMIgg5NahnWf/aCWtUSez5llnvIT7OEMvgnig9AFlZLt8qEFYBLwjhg\nKwNRIqbdTWReH9E9aBGSkkRLPZ5IJPaUJQlpXrcG8g8UIt3iprDfX9cFXZpXzBc/sk820lIIt/Vu\n6YWIjAFpqck3UCmTiBrp7hzEd1/UCgs2H5AhUtLjpw+IFUDIGHtzD8zfvL/CpjAzbtQJn6uSmoI7\n+/pjG7XCF7/og1NFJUGLEXqcpnH2knsvORcZ6anluYyc8tQw+1ljs+tnYOWOI/Ao6aYlWjesieXb\nDhtWxpPJk0Pb47ef/4guzawXAHILK4CQUbtGFQzu1DhoMaTSo8U5QYuQEEy4PRcfLdoWWL4gPVNL\nelqq5fq6snn/zp5Ytu1woKbM9+/siR93HkGGD7WjuzSrg/893M/zftSwAmAYl1SJmX7c2qpb1KuB\n0UPaxW/oEUHUoTajXs10XBbAfgg1dWpURb8ciWkeQkb41pxMqOhzXn0AQIbD2qdR4JXru+Dn/c91\nXFKTYYKC72rGlJeu7YQHBrRG7ZCW5wsDDTKr4bEr2wctBsPYhlcAjCnpaalSimAzDBM+XCkAIqpL\nRN8S0cbYf11vHxHVIaLPiGgdEa0logvd9BsmLji3HprXrY5HBuUELQrDMIwt3K4ARgOYJoTIATAt\n9lyP1wFMFUK0A9AVwFqX/YaG2tWrYM5vBqBTU/9CtxiGYWTgVgEMB/BB7PEHAK7RNiCi2gD6A3gX\nAIQQZ4QQh132yzAMw7jErRO4oRBid+zxHgB6MVutABQAeJ+IugJYCuBhIcQJl30zDGOBaY9ejCMn\niwzfVzKncuqQ6BFXARDRdwAa6bz1hPqJEEIQkV4gcRqAHgAeFEIsJKLXUWYqesqgv1EARgFAixYt\n4onHMEwc4m0se3BADpqdUx1XdUmuDYhMfOIqACHEIKP3iGgvETUWQuwmosYA9uk02wFghxBiYez5\nZzD2FUAIMQHABADIzc0N184UhklCqqal4KaePNmKIm59AF8CGBl7PBLAJG0DIcQeANuJqG3spYEA\n1rjsl2EYhnGJWwUwBsBlRLQRwKDYcxBREyKaomr3IIAPiWglgG4AXnLZL8MwDOMSV05gIcQBlM3o\nta/vAnCl6vkKALlu+mIYhmHkwjuBGYZhIgorAIZhmIjCCoBhGCaisAJgGIaJKJwOmmEYRsOE288H\n+VmcNyBYATAMw2i4vKNe8oPkg01ADMMwEYUVAMMwTERhBcAwDBNRWAEwDMNEFFYADMMwEYUVAMMw\nTERhBcAwDBNRWAEwDMNEFBIivEW3iKgAwFaHH68PYL9EcRIB/s7JT9S+L8Df2S4thRBZVhqGWgG4\ngYiWCCEiVYOAv3PyE7XvC/B39hI2ATEMw0QUVgAMwzARJZkVwISgBQgA/s7JT9S+L8Df2TOS1gfA\nMAzDmJPMKwCGYRjGhKRTAEQ0mIjWE1EeEY0OWh6vIaLmRDSDiNYQ0WoiejhomfyCiFKJaDkRfRW0\nLH5ARHWI6DMiWkdEa4nowqBl8hoi+mXsul5FRP8kompByyQbInqPiPYR0SrVa3WJ6Fsi2hj7f44X\nfSeVAiCiVABvAhgCoAOAm4moQ7BSeU4xgEeFEB0A9AZwfwS+s8LDANYGLYSPvA5gqhCiHYCuSPLv\nTkRNATwEIFcI0QlAKoARwUrlCX8FMFjz2mgA04QQOQCmxZ5LJ6kUAIBeAPKEEJuFEGcAfAxgeMAy\neYoQYrcQYlns8TGUDQpNg5XKe4ioGYChACYGLYsfEFFtAP0BvAsAQogzQojDwUrlC2kAqhNRGoAa\nAHYFLI90hBCzARzUvDwcwAexxx8AuMaLvpNNATQFsF31fAciMBgqEFE2gO4AFgYriS/8GcBvAJQG\nLYhPtAJQAOD9mNlrIhFlBC2UlwghdgJ4FcA2ALsBHBFCfBOsVL7RUAixO/Z4D4CGXnSSbAogshBR\nTQCfA3hECHE0aHm8hIiGAdgnhFgatCw+kgagB4BxQojuAE7AI7NAWIjZvYejTPk1AZBBRLcFK5X/\niLJQTU/CNZNNAewE0Fz1vFnstaSGiKqgbPD/UAjxRdDy+EBfAFcTUT7KzHwDiOgfwYrkOTsA7BBC\nKKu7z1CmEJKZQQC2CCEKhBBFAL4A0CdgmfxiLxE1BoDY/31edJJsCmAxgBwiakVEVVHmMPoyYJk8\nhYgIZXbhtUKI14KWxw+EEI8JIZoJIbJR9htPF0Ik9cxQCLEHwHYiaht7aSCANQGK5AfbAPQmohqx\n63wgktzxreJLACNjj0cCmORFJ2leHDQohBDFRPQAgK9RFjHwnhBidcBieU1fALcD+JGIVsRee1wI\nMSVAmRhveBDAh7HJzWYAdwUsj6cIIRYS0WcAlqEs2m05knBXMBH9E8AlAOoT0Q4AzwAYA+BTIrob\nZRmRb/Skb94JzDAME02SzQTEMAzDWIQVAMMwTERhBcAwDBNRWAEwDMNEFFYADMMwEYUVAMMwTERh\nBcAwDBNRWAEwDMNElP8HdqsE7qK7sh4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11eade358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(ts_all.size)*delta_t,ts_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11eb9b400>]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEACAYAAACznAEdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF4BJREFUeJzt3X+MHOV9x/HPh8MG56ehXC1yQE0Si9TpJRidwAipgipu\nMLT1BakJCAptIyyHEjWhQjUFFdIexUqaH0JFRKCiQHFDSIMctzi4FkmUCmHEORAbSlwfDgEuju3I\nMVBw8a9v/7g5OJ/3bp+9ndndmX2/pJVvZ57Ze3ZY9nPzfWaecUQIAIBj2t0BAEBnIBAAAJIIBABA\nhkAAAEgiEAAAGQIBACCJQAAAZAgEAIAkAgEAkCEQAACSpGPb3YFGnHTSSTF//vx2dwMASmXTpk2/\niojeeu1KFQjz58/X8PBwu7sBAKVi++cp7SgZAQAkEQgAgAyBAACQRCAAADIEAgBAUsnOMpqJNU+N\n6kvrt+oXe/fpfXPn6PqPn6HBRX3t7hYAdJxKB8Kap0Z1w0NbtO/AIUnS6N59uuGhLZJUWCg0G0AE\nGIB2qXQgfGn91rfCYNy+A4f0pfVbC/mSbTaAZrJ9SoAQMgBSVHoM4Rd79zW0vFnTBVAR248HyOje\nfQq9HSBrnhpNbrPmqVGdt+r7On3lwzpv1feP2BZAd6n0EcL75s7RaI0v//fNnVPI72s2gBrdPuUI\nqF7ITHVEMr7txKOKWss40gCqo9KBcP3HzzjiC0+S5szqeevLLW/NBlCj26cEyHRtpgqLW9Y+qzcP\nHj4iKK7/9k8kSwcOxVvLbnhoi4Z/vkc/+Onut0Ligg/1HvGc0ADKo9Ilo8FFfbrtkn71zZ0jS+qb\nO0e3XdJf2BfU9R8/Q3Nm9RyxrJEAanT7qYJi4vLp2kwVFnv3HTgqKA4cjrfCYNy+A4e0euOLR5Sj\n7p/0/Ppv/0SL/u4/KUkBJeCIqN+qQwwMDESnT27XyrOMJg9CS2MBMjH0pmvzpfVbax6RFGnWMda7\njj9We984oPfOmSVb2vvGAY4mgALZ3hQRA3XbEQjl1sxZRlOFxfGzjtGv3zjQ6rciSwpJcwkKIFcE\nApLUCgtJRwXFrGN8xBiC9PYXeCu8Y9YxOm5WDyEBzACBgKZMFRQTl13woV59Z9PoUeMNrXLCO2bp\n5j/8MMEA1EEgoCUmB8fEs4zeO2eWXt9/8KjB6LwdY+lwjJ00wJEDcDQCAR1hYmAQEEB7EAjoSJMD\nwpZ+/caBwscjKC+hmxEIKJXJQbH/4CG9ceBw7r+Howd0IwIBpdeqkCAcUHUEAippzVOjumXts9q7\nL9/rJCzp8sWnaWiwP9fXBToBgYBKGz96GN27L/fxB44YUDUEArpKEQHxztk9uvUTxc19BbQKgYCu\nlndpiaMGlBmBAGQmHj00i6MGlBGBANSQVzgQDCgTAgGo46Y1W7R644tNjTcQDCiD1ECo9A1ygOkM\nDfbrq586U31N3FL19f2H9LlvPa0P/+0j3PwHpccRApAZuz/EZu1r4uI3jhjQiThCABo0uKhPz/39\nUn2tiaOG8SOGm9Zsybl3QPE4QgCm0cxRA0cL6BQMKgM5IhhQZrmWjGxfaHur7RHbK2ust+3bs/Wb\nbZ+VLT/V9g9s/7ftZ23/5YRtTrS9wfa27N8TGnmDQCtNLCfNmdVYpZUyEsqi7ifbdo+kOyQtlbRQ\n0mW2F05qtlTSguyxXNKd2fKDkv4qIhZKWizpLyZsu1LSoxGxQNKj2XOgozUTDPdvfJGzkdDRUj7R\nZ0saiYjtEbFf0gOSlk1qs0zSfTFmo6S5tk+OiB0R8WNJiojXJD0nqW/CNvdmP98rabDJ9wK0zHgw\nXLH4tIa24zRVdLKUQOiT9NKE5y/r7S/15Da250taJOmJbNG8iNiR/fxLSfOSegx0kKHBfspIqIyW\nnHZq+12SviPpcxHx6uT1MTayXXN02/Zy28O2h3fv3l1wT4HGNVtGuvzuxwvqGdCYlE/vqKRTJzw/\nJVuW1Mb2LI2FweqIeGhCm522T87anCxpV61fHhF3RcRARAz09vYmdBdoj5mWkR57fo/OuOl7lJDQ\ndimB8KSkBbZPtz1b0qWS1k5qs1bSldnZRoslvRIRO2xb0j9Lei4ivlJjm6uyn6+S9N0Zvwugg8yk\njPTmwcOUkNB2dT+xEXFQ0rWS1mtsUPjBiHjW9grbK7Jm6yRtlzQi6W5J12TLz5P0J5J+z/bT2eOi\nbN0qSUtsb5P0sew5UAkzPVqghIR24sI0oGAzuajt2GOsf/zjj3IxG3LBXEZAh5g46JxaRTp4OCgh\noeUIBKBFBhf1ads/XKzzPnBi8jaUkNBKBALQYquvPrehsYXHnt+jD/7NOs5CQuEIBKANxs9EooSE\nTkIgAG0yXkJa8JvvTN7m/o0vEgooDIEAtNmG685nXAEdgUAAOsDqq89tqIT02PN7tOQrPyy0T+g+\nBALQIRotIW3b9TqhgFwRCECHaaSERCggTwQC0IEaOTV1267XOS0VuSAQgA41NNifHArjp6USCmgG\ngQB0sPHrFVIRCmgGgQB0uMFFfYQCWoJAAEpgPBRS/4clFDATBAJQEoOL+rR9VfppqUx1gUYRCEDJ\nbLju/ORQuH/jixwpIBmBAJRQI6HwuW89XXBvUBUEAlBSjYTC+294mCMF1EUgACWWGgqHgzEF1Ecg\nACXHmALyQiAAFcCYAvJAIAAVseG68zXv3bOT2n7k5kcK7g3KiEAAKuSJG5ckhcKrbx7SObduaEGP\nUCYEAlAxT9y4JKl8tPO1/UydjSMQCEAFpY4pcD8FTEQgABWVOqawbdfr3KMZkggEoNKeuHGJnNDu\nsef3cI0CCASg6r6aOHU21yiAQAAqbnBRX/Kd17hGobsRCEAXaOR2nFyj0L0IBKBLDA3267wPnFi3\nHdcodC8CAegiq68+l2sUMCUCAegynI6KqRAIQBd64sYles9xPXXbPfb8Hs486iJJgWD7QttbbY/Y\nXlljvW3fnq3fbPusCevusb3L9jOTtrnF9qjtp7PHRc2/HQCpNn/hwqRrFDjzqHvUDQTbPZLukLRU\n0kJJl9leOKnZUkkLssdySXdOWPcNSRdO8fJfjYgzs8e6BvsOoEmp1ygwyNwdUo4QzpY0EhHbI2K/\npAckLZvUZpmk+2LMRklzbZ8sSRHxI0l78uw0gHwMLupLOvNo52v7GU/oAimB0CfppQnPX86WNdqm\nls9mJaZ7bJ+Q0B5AzlLPPGI8ofraOah8p6T3SzpT0g5JX67VyPZy28O2h3fv3t3K/gFdI/XMI8YT\nqi0lEEYlnTrh+SnZskbbHCEidkbEoYg4LOlujZWmarW7KyIGImKgt7c3obsAZiJ1IjzGE6orJRCe\nlLTA9um2Z0u6VNLaSW3WSroyO9tosaRXImLHdC86PsaQ+YSkZ6ZqC6A1UgaZGU+orrqBEBEHJV0r\nab2k5yQ9GBHP2l5he0XWbJ2k7ZJGNPbX/jXj29v+pqTHJZ1h+2Xbn85WfdH2FtubJV0g6fN5vSkA\nM5M6yMx4QjU5Itrdh2QDAwMxPDzc7m4AlfeRmx/Rq28eqtvuhVUXt6A3aJbtTRExUK8dVyoDOMrm\nL0x16dCRGE+oFgIBQE1fYzyh6xAIAGoaXNTH9QldhkAAMKUN153PfEddhEAAMK3U+Y64f0L5EQgA\nppV6T+Ztu16ndFRyBAKAuoYG+5PGEygdlRuBACBJ6ngCp6KWF4EAIFnq1BY3rdnSgt4gbwQCgGSp\np6Lev/HFFvQGeSMQADQktXT0kZsfKbwvyBeBAKBhKaWjV988ROmoZAgEAA1LnRWV0lG5EAgAZmT1\n1efqPcf11G3HBWvlQSAAmLGUWVG5YK08CAQATUm5ipkL1sqBQADQlKHBfh2TcNoRpaPORyAAaNpX\nPln/rCNKR52PQADQtNQL1igddTYCAUAuNlx3flI7Skedi0AAkJuU225SOupcBAKA3FA6KjcCAUCu\nUktHTGvReQgEALlLuTaBaS06D4EAIHdDg/2a9+7ZddtxM53OQiAAKMQTNy6p22bna/sZYO4gBAKA\nwjCtRbkQCAAKMzTYr+N76s9rwbUJnYFAAFCon956Ud02XJvQGQgEAIXj2oRyIBAAFI5pLcqBQADQ\nEikDzJSO2otAANASQ4P9SbfcpHTUPgQCgJZJueWmROmoXZICwfaFtrfaHrG9ssZ62749W7/Z9lkT\n1t1je5ftZyZtc6LtDba3Zf+e0PzbAdDpUktHaL26gWC7R9IdkpZKWijpMtsLJzVbKmlB9lgu6c4J\n674hqdafBSslPRoRCyQ9mj0HUHGp1ybMX/lwC3qDiVKOEM6WNBIR2yNiv6QHJC2b1GaZpPtizEZJ\nc22fLEkR8SNJe2q87jJJ92Y/3ytpcCZvAED5pFybIFE6arWUQOiT9NKE5y9nyxptM9m8iNiR/fxL\nSfMS+gKgIigddZ6OGFSOiJAUtdbZXm572Pbw7t27W9wzAEVJLR2dTumoZVICYVTSqROen5Ita7TN\nZDvHy0rZv7tqNYqIuyJiICIGent7E7oLoCxSSkch6fK7Hy++M0gKhCclLbB9uu3Zki6VtHZSm7WS\nrszONlos6ZUJ5aCprJV0VfbzVZK+20C/AVRESunosedrDUMib3UDISIOSrpW0npJz0l6MCKetb3C\n9oqs2TpJ2yWNSLpb0jXj29v+pqTHJZ1h+2Xbn85WrZK0xPY2SR/LngPoMqmlow/eQOmoaB4r35fD\nwMBADA8Pt7sbAAqQcprpFYtP09Bgfwt6Uy22N0XEQL12HTGoDADnfeDEum24D3OxCAQAHWH11eeq\nfuGI+zAXiUAA0DF+turium24D3NxCAQAHYWb6bQPgQCgo6TeTIfSUf4IBAAdJ+XaBEpH+SMQAHSc\nocF+HZswwkzpKF8EAoCONHJb/QFmiRlR80QgAOhYKdcmMCNqfggEAB0r9doEbqaTDwIBQEdLuTZB\nonSUBwIBQMejdNQaBAKAjrf66nOT2lE6ag6BAKAUXkgsHXEznZkjEACURkrpiJvpzByBAKA0UktH\n3Id5ZggEAKWSUjoKSTet2VJ8ZyqGQABQOikzonIzncYRCABKJ3VG1A/duK7YjlQMgQCglL72qTPr\ntvm/Q8GMqA0gEACU0uCiPr3nuJ667ZgRNR2BAKC0Nn/hwqR2H7n5kYJ7Ug0EAoBSS7k24dU3D1E6\nSkAgACi11GsTKB3VRyAAKL3UaS2YEXV6BAKASki5NoEZUadHIACohNRrE5gRdWoEAoDKSC0dnXPr\nhoJ7Uk4EAoBKmffu2XXb7Hxtfwt6Uj4EAoBKeeLGJUntKB0djUAAUDmcdTQzBAKASkopHXHW0ZEI\nBACVROmocQQCgMriPsyNSQoE2xfa3mp7xPbKGutt+/Zs/WbbZ9Xb1vYttkdtP509LsrnLQHA21Jm\nROU+zGPqBoLtHkl3SFoqaaGky2wvnNRsqaQF2WO5pDsTt/1qRJyZPbiTBYDcpc6I+sEbKB2lHCGc\nLWkkIrZHxH5JD0haNqnNMkn3xZiNkubaPjlxWwAo1BWLT6vb5mCo62dETQmEPkkvTXj+crYspU29\nbT+blZjusX1Ccq8BoAFDg/1yQrtunxG1nYPKd0p6v6QzJe2Q9OVajWwvtz1se3j37t2t7B+ACvlZ\n4gBzN99MJyUQRiWdOuH5KdmylDZTbhsROyPiUEQclnS3xspLR4mIuyJiICIGent7E7oLALWlzIj6\n6puHWtCTzpQSCE9KWmD7dNuzJV0qae2kNmslXZmdbbRY0isRsWO6bbMxhnGfkPRMk+8FAKbFjKjT\nqxsIEXFQ0rWS1kt6TtKDEfGs7RW2V2TN1knaLmlEY3/tXzPdttk2X7S9xfZmSRdI+nx+bwsAaku9\nNuH0LgwFR0S7+5BsYGAghoeH290NACV3zq0bkmY8vWLxaRoa7G9Bj4ple1NEDNRrx5XKALpO6rQW\n9298seCedBYCAUBXSi0dddN4AoEAoGulTGshdU8oEAgAulbqtBZSd4QCgQCgq6WWjqTqz3dEIADo\neilzHUnVn++IQADQ9YYG+3V8T8psR9We74hAAABJP701/ZYsVR1PIBAAINPIeEIVQ4FAAIAJujkU\nCAQAmKRbQ4FAAIAaUqbKHleVUCAQAKCG1Kmyx1UhFAgEAJhCI6UjqfyhQCAAwDS6KRQIBACoo1tC\ngUAAgATdEAoEAgAkmkkolGlCPAIBABrwwqqLdd4HTkxufzDKc7RAIABAg1ZffW7yDKnjyhAKBAIA\nzMDQYL/S5kd9W6eHAoEAADP0swbHFKTODgUCAQCa0OhAs9S5oeCIaHcfkg0MDMTw8HC7uwEAR7n8\n7sf12PN7ZrTtTEKlEbY3RcRA3XYEAgDkZ6Z//R9raeS2YoIhNRAoGQFAjmb6134nnJ5KIABAzpop\nAc1f+XDbgoFAAIACNDsuMH/lwzrn1g059SYNgQAABWk2FHa+tr+lRwsEAgAUKI8ziFpVRuIsIwBo\noWa/2GcSMJxlBAAd6IVVFzc96FwUAgEA2qCZYLhpzZacezOGQACANppJKHzziZcK6EliINi+0PZW\n2yO2V9ZYb9u3Z+s32z6r3ra2T7S9wfa27N8T8nlLAFAujR4tHCpo7LduINjukXSHpKWSFkq6zPbC\nSc2WSlqQPZZLujNh25WSHo2IBZIezZ4DQNdKDYUeNzrxdpqUI4SzJY1ExPaI2C/pAUnLJrVZJum+\nGLNR0lzbJ9fZdpmke7Of75U02OR7AYDSSzlauOycUwv53SmB0CdpYsHq5WxZSpvptp0XETuyn38p\naV5inwGg8qYKhSsWn6ahwf5CfuexhbxqgyIibNcsitlerrEylE47rbFb1gFAmRU9LfZkKUcIo5Im\nHp+cki1LaTPdtjuzspKyf3fV+uURcVdEDETEQG9vb0J3AQAzkRIIT0paYPt027MlXSpp7aQ2ayVd\nmZ1ttFjSK1k5aLpt10q6Kvv5KknfbfK9AACaULdkFBEHbV8rab2kHkn3RMSztldk678uaZ2kiySN\nSHpD0p9Nt2320qskPWj705J+LumTub4zAEBDmMsIACqOuYwAAA0p1RGC7d2S9kp6ZdKq9yYsO0nS\nr4rr3VFq9anI7VPaT9dmqnUp+7bWMvZ3fvu71vJa7Vq5z5vd342+RlH7e6p1VfuM/1ZE1D8rJyJK\n9ZB010yWSRpudz+L3D6l/XRtplrH/m7//p5i/9b6b9Cyfd7s/m70NYra39Psy8p9xlMeZSwZ/XsT\ny1qp2d/f6PYp7adrM9U69vfM2+e1v2stL/v+bvQ1itrfU62r4me8rlKVjJphezgSBlWQD/Z367HP\nW6uK+7uMRwgzdVe7O9Bl2N+txz5vrcrt7645QgAATK+bjhAAANMgEAAAkggEAECGQJBk+3zb/2X7\n67bPb3d/uoHtd9oetv0H7e5L1dn+7eyz/W+2P9Pu/nQD24O277b9Ldu/3+7+pCp9INi+x/Yu289M\nWj7tfaAnCUn/K+l4jd3EB1PIaX9L0l9LerCYXlZHHvs7Ip6LiBUam0DyvCL7WwU57fM1EXG1pBWS\nPlVkf/NU+rOMbP+uxr7M74uI38mW9Uj6H0lLNPYF/6SkyzQ24+ptk17izyX9KiIO254n6SsRcXmr\n+l82Oe3vj0r6DY0F8K8i4j9a0/vyyWN/R8Qu238k6TOS/iUi/rVV/S+jvPZ5tt2XJa2OiB+3qPtN\n6Yg7pjUjIn5ke/6kxW/dy1mSbD8gaVlE3CZpuhLFryUdV0Q/qyKP/Z2V5d4paaGkfbbXRcThIvtd\nVnl9viNiraS1th+WRCBMI6fPuDU2xf/3yhIGUgUCYQq17uV8zlSNbV8i6eOS5kr6p2K7VkkN7e+I\nuFGSbP+psqOzQntXPY1+vs+XdInG/thZV2jPqquhfS7ps5I+Jum9tj8YY/eN6XhVDYSGRMRDkh5q\ndz+6TUR8o9196AYR8UNJP2xzN7pKRNwu6fZ296NRpR9UnkLKfaCRH/Z3a7G/W68r9nlVAyHlPtDI\nD/u7tdjfrdcV+7z0gWD7m5Iel3SG7ZdtfzoiDkoav5fzc5IejLfv5YwmsL9bi/3det28z0t/2ikA\nIB+lP0IAAOSDQAAASCIQAAAZAgEAIIlAAABkCAQAgCQCAQCQIRAAAJIIBABA5v8BaHIW8BHqD9kA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11e9e40b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = ts_all.size\n",
    "# calculate autocorrelation function\n",
    "f = np.fft.rfft(ts_all)\n",
    "acf = np.fft.irfft(f * np.conjugate(f))\n",
    "acf = np.fft.fftshift(acf) / N\n",
    "autocorr = acf[int(N / 2):]\n",
    "\n",
    "autocorr = autocorr[:2000]\n",
    "t = np.arange(autocorr.size)*delta_t\n",
    "plt.semilogx(t, autocorr,\"o\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fitcorr(t,kappa,sigma,ga,eta):\n",
    "    corr=[np.sum(np.exp(-gamma(q_list,kappa=kappa, sigma=sigma, ga=ga, eta=eta)*tt)/(kappa*q_list**4+sigma*q_list**2+ga)) for tt in t]\n",
    "    return corr\n",
    "\n",
    "fit_corr_model=lm.Model(fitcorr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Model]]\n",
      "    Model(fitcorr)\n",
      "[[Fit Statistics]]\n",
      "    # function evals   = 9\n",
      "    # data points      = 2000\n",
      "    # variables        = 1\n",
      "    chi-square         = 0.000\n",
      "    reduced chi-square = 0.000\n",
      "    Akaike info crit   = -34792.507\n",
      "    Bayesian info crit = -34786.906\n",
      "[[Variables]]\n",
      "    kappa:   10 (fixed)\n",
      "    sigma:   125 (fixed)\n",
      "    ga:      25 (fixed)\n",
      "    eta:     0.24079124 +/- 0.000175 (0.07%) (init= 0.25)\n",
      "[[Correlations]] (unreported correlations are <  0.250)\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEACAYAAACznAEdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//HXJyshgYRA2IIURMqiUNCoILZolQJKS7R1\n4VutVivVVq1fa1uU+q38flRta939UqVasVLBVkVUlCK4FAtoKP5YxAVRZF8SQliyTs7vj0zGmckk\nmcBkMjN5Px+PeTBz7rkzZ64x79xzzj3XnHOIiIgktXUDREQkNigQREQEUCCIiIiXAkFERAAFgoiI\neCkQREQEUCCIiIiXAkFERIAwA8HMJpjZR2a2ycymhdhuZvagd/taMzvZW36cmb1hZh+Y2QYz+5nf\nPneY2XYze9/7OC9yX0tERFoqpbkKZpYMPAKMA7YB75nZQufcB37VJgIDvY/TgVnef2uAnzvn/mNm\nnYDVZrbEb9/7nHP3RO7riIjI0Wo2EIDTgE3Ouc0AZjYPmAz4B8Jk4ClXtw7GSjPLMbNezrmdwE4A\n59xBM9sI5AftG7Zu3bq5fv36Hc2uIiLt1urVq/c55/KaqxdOIOQDW/1eb6Pur//m6uTjDQMAM+sH\njARW+dW7wcx+ABRRdyaxv6mG9OvXj6KiojCaLCIi9cxsSzj1ojKobGZZwHPATc65Mm/xLOB4YAR1\nwfHHRvadamZFZla0d+/eaDRXRKRdCicQtgPH+b3u4y0Lq46ZpVIXBnOdc8/XV3DO7XbOeZxztcBs\n6rqmGnDOPeacK3DOFeTlNXvGIyIiRymcQHgPGGhm/c0sDbgUWBhUZyHwA+9so1HAAefcTjMz4HFg\no3PuXv8dzKyX38sLgPVH/S1EROSYNTuG4JyrMbPrgcVAMvCEc26DmV3r3f4nYBFwHrAJOAL80Lv7\nGOByYJ2Zve8tu805twj4vZmNABzwOfDjiH0rERFpMYunG+QUFBQ4DSqLiLSMma12zhU0V09XKouI\nCKBAEBERr3CuQ0gIN954Izt37iQ1NZW0tDRSU1MZPnw4N9xwAwD33HMPhw8f9m1LS0tj0KBBjB8/\nHoAFCxZQW1sbsD0/P5+vfvWrAHzwwQekpqb6tnXo0IGsrCzS0tLa7DuLiLREuwmErVu38vHHH1NV\nVUV1dTVVVVUcOnTIFwgPPfQQX3zxRcA+F198sS8QrrjiCsrKygK2X3XVVTz++OMADBs2jNra2oDt\nN954Iw888AAVFRXk5+fTsWNHMjIyfP/+6Ec/4uqrr6asrIwbb7zRt63+ce6553Lqqady5MgR/v3v\nf9O5c2c6depEp06d6Ny5M1lZWSQl6SRPRCKj3QTCCy+80OT2LVu2UFtb6wuL6upqkpOTfdtXrVpF\nZWWlb3tVVRU9evQAwDnHM888E7CtsrKS4cOH+7ZPmTKFI0eOUF5e7vs3JaXu8B8+fJg33niD8vJy\n3/ba2lo6dOjAqaeeypYtWxg3blyDNj/66KNMnTqV9evXc+mll5Kbm0vXrl19j8suu4xhw4ZRUlLC\n+vXr6dq1q6+OzlxEJFi7CYRwJCUlkZ6eTnp6eoNtgwcPbnQ/M+Piiy9udHtGRgYPP/xwo9t79erF\nli1fXlnunKOystL3um/fvrz99tuUlZVx8OBB37+jR48GIDk5mUGDBlFcXMymTZtYtWoVxcXFjBkz\nhmHDhrFixQomTZoU8JlZWVm8/PLLjB07luXLl/Poo4/Su3dvevfuTX5+Pr1792bEiBF07Nix0XaL\nSGJJ+EBYsGY7f1j8ETtKy+mdk8Evxg+icGR+TH+emdGhQwff68zMTL7+9a83Wn/IkCE899xzAWXO\nOeqnFI8aNYolS5ZQXFxMSUkJxcXFFBcX07dvXwD27dvH8uXL2bFjB1VVVb73WLduHSeddBJPPvkk\n9913ny8o8vPz6devHxdddBFZWVkt+m4iErsS+jqEBWu2c+vz6yiv9vjKMlKTuevCYa0SCpH4vJYG\nSiQDzzlHSUkJ27dvZ8eOHYwdO5aMjAxefPFFHn/8cXbs2MGOHTvYtWsXzjlKS0vJzs5mxowZPP30\n0xx//PEcf/zx9O/fnwEDBlBYWBjQ7SYibSPc6xASOhDG3L2M7aXlDcrzczJ4Z9o3I9m0iHxeSwMl\n3PpNhcbRBEpVVRVbt25lwIABADzzzDO88MILfPbZZ2zevJmSkhI6derEgQMHMDOmT5/Ohg0bGDx4\nMEOGDGHIkCGceOKJZGZmNntMROTYhRsICd1ltCPEL+emytv68/6w+KOAX+4A5dUe/rD4o5C/pMOp\nHxwa20vLufX5db76TW1rLCjS0tJ8YQAwZcoUpkyZ4nt94MABdu7cSd1SVnVnHp988gmLFi2iuroa\nqBuT2bhxIwDz5s0jOzubESNG0LNnT99+IhJdCR0IvXMyQv7F3jsnIyY/r6WBEk55U6FR/zx42x0L\nN1BZU9sgKIq2lPDGh3sDQqL+M4KDIzs72/eed955J3feeSc1NTVs3ryZDz74gJqaGt/2m2++mZ07\n626d0bt3b0aNGsX3vve9gJARkdaX0IHwi/GDQnap1P8ii7XPa2mghFP/aM5aSsurG5SVV3uYu/IL\n6jsYt5eW84u//z8wqPa4gLIZL22g9Eg1vXMyOHtwXoMQKSwsDHjvjRs3snbtWtasWcO7777LypUr\n6devH1OmTKGqqoqxY8dyyimnMHr0aMaOHUufPn0abbuIHL2EHkOA+Jpl1BpjCE2NawAht7Wm1CQj\nq0MKpUeqyc5IxQxfePgfK4/HQ3JyMjt27ODyyy/n3Xff5dChQwAMHDiQe++9t8FUWhEJTYPKcSrS\ns4yaCg0g5LYOqUnsP9LwLKG1GXVroeeH+B4ej4f169ezbNkyli1bxvTp0xk1ahSvvvoq06ZN49xz\nz2XSpEmceeaZpKamRr3tIrFMgSA+LZ1lBA2Dov6XdbTUf16yGR7nQoYEwNKlS7nzzjtZvnw5VVVV\nZGdnM2HCBGbNmkWXLl2i2GKR2KVAkGMSHBRnD87judXbA0IiNckCxhBaW1MhcejQIV5//XVeeukl\nioqKWLNmDUlJSTz66KMAfPe736Vbt25RaadIrFEgSMQ1djZRX5adkcrhqpqoBUS9JINaF7qr6Zxz\nzmHZsmUkJydz7rnncskll3DBBReQk5MT1TaKtCUFgrSJUGcW9bOMohUY/gFxy7e+Sv+kfcyfP5/5\n8+fz2WefceWVV/KXv/wFgJqaGt8igyKJSoEgMck/MOpnGe0/Ut3qYxRdOqbyP5OG0sezk8zMTIYO\nHcratWsZP348V155JVdddRUDBw5sxRaItB0FgsSV+qDYXlrequHgf/ZwUX8Pb/ztYV555RU8Hg9j\nx47luuuu48ILL9RMJUkoCgSJW/7hUD+A3BohYcD3R/Xlp6d3Y86cOcyePZs9e/awbds2srOzqa2t\n1Q2IJCEoECShtPYZRJeOqdx+/hAGdSjjxBNPxDnHmDFjOOmkk7jlllt8t0oViUcKBElorRkQ+TkZ\n3PiNvix98vf89a9/paqqiosvvpjbbruNYcOGRfCTRKJDgSDtSmsERGZaMr/4Rk8+fH0ejzzyCIcO\nHeK1117z3WdbJF4oEKRdW7BmO3cs3BByob6j0T2tmkGlq/jzH+4gPT2dxYsXk5+fz0knnRSR9xdp\nTQoEES//s4djlZmWzMzCk/j15ePZuHEjV1xxBTNmzPDdjlQkFikQREKI1JlDB88RBu9eymvPPgnA\n9ddfz/Tp07V+ksQkBYJIE+pWgV1LeXXtMb1PWnkxx32+iDdf/jsLFizg29/+doRaKBI54QaCJllL\nu1Q4Mp+N/3ci918ywndviKNRldGVT4dczvHXzaYmfyQATzzxBO+++26kmioSNTpDEPGKxFlDx+Ra\nSuZcz+7tX3D11Vfzu9/9jq5du0awlSItpzMEkRaKxFnDEU8SaRfdQ9fR3+XJOXMYMmQIzz77LPH0\nh5e0XzpDEGnCsZw1VO39nP2vPkDV7k/ZuHGjrnaWNqMzBJEI8D9ryEht2f8uaXn96H7ZPeRdMpMZ\nbxUDsGLFCp0tSMwK6yfczCaY2UdmtsnMpoXYbmb2oHf7WjM72Vt+nJm9YWYfmNkGM/uZ3z65ZrbE\nzD7x/qv5ehKzjjYYLCmZDn2H886nJfT94f2cccYZXHTRRZSUlLRia0WOTrM/2WaWDDwCTASGAlPM\nbGhQtYnAQO9jKjDLW14D/Nw5NxQYBfzUb99pwFLn3EBgqfe1SEw7ljMG6z6AnLN+yAsLXmT48OEs\nW7aslVopcnTC+Yk+DdjknNvsnKsC5gGTg+pMBp5ydVYCOWbWyzm30zn3HwDn3EFgI5Dvt88c7/M5\nQOExfheRqKkPhstGhX+FslkS2ad/l+6X3cOeCuOcc87h9ttvb8VWirRMOIGQD2z1e72NL3+ph13H\nzPoBI4FV3qIezrmd3ue7gB5htVgkhswsHMb9l4wgJyP8G+qk9zyBnlfcT9aIifx59X4WrNneii0U\nCV9UbiZrZlnAc8BNzrmy4O3OOWdmIUfazGwqdd1QWi9GYlLhyHwKR9b9/bNgzXZ+8ff3aW5SUlJq\nB7qO/ykAN81/n6fnPcvUc07kW9/6Vms3V6RR4ZwhbAeO83vdx1sWVh0zS6UuDOY65573q7PbzHp5\n6/QC9oT6cOfcY865AudcQV5eXhjNFWk7hSPz+eTO8xkzIDfsfZxzvPT0o4yfMIGZM2dSW3tsy2mI\nHK1wAuE9YKCZ9TezNOBSYGFQnYXAD7yzjUYBB5xzO83MgMeBjc65e0Psc4X3+RXAi0f9LURizNxr\nRnP/JSMIZ9zZzOhxyW/JHDKW22+/ndPGjqO0tLT1GykSJKwL08zsPOB+IBl4wjn3WzO7FsA59yfv\nL/6HgQnAEeCHzrkiMzsT+BewDqj/s+c259wiM+sKPAv0BbYAFzvnmpyLpwvTJB59f/YK3vm0+Wmm\nzjkOrVlEydLH6Ny9D1s/Wkvnzp2j0EJJdFrtVCSG/HrBOp5e+UVYdSu+WEfF5+9z6veuZcnNZ7Vu\nw6RdUCCIxJijWQbD7fuM7/Sp4qE7bm7Flkmi09IVIjHG/6K2cK9pK373RR6e8XNOPO8KLXkhrU6B\nIBJlLZmJ1HXCDWR9bTwfvPoU+aefR01NTRRaKO2VAkGkjcy9ZnSzVzpbUjK5468ne8wUdr73GtmD\nz2D+vzdFqYXS3igQRNpQ/ZXOTf2PaGbknPl9cr/1E2pra/nlC+v59YJ1UWujtB8aVBaJEeFMT3Wu\nFrMkPOUH+e6IHtz/w7Oj1DqJZxpUFokzYXUhWRLOOfYuuJNHbrmcGx9/PUqtk/ZAgSASQ8LuQvr6\nZXgOlTDrlsspvPv5JmqLhE+BIBJjCkfms/nu8xnYPbPROh36nEiPS2ZSW36Ql++6ljG3zY1iCyVR\nKRBEYtSSm89qcmpqeu9BdL/0t7jqSlb/9S5Om/nPKLZOEpECQSSGNTeukN7zBHpMuZNuk25mz6Fq\nTv/tkii2ThKNAkEkxs0sHNbkmUJaXj9Ssnvgaj18tOgJRv7qmSi2ThKJAkEkDsy9ZnSzVzbXlO6i\nrOhF1v/5l3zlxrm6E5u0mAJBJE7U32Ohsf9pU3Pz6f693+Ap28vOeb/mhif/pQvYpEUUCCJxpH4G\nUo9OaSG3dzjuJPK+ezvVJdvZ8/f/4am3P9KZgoRNgSASh1ZNH9doKGT0G0He5GlUF2+jas9n3DT/\n/Si3TuKVAkEkTjUVCh0Hnk7+tY/Toc9QAI6/9RWdKUizFAgicaypUEjOqLv95sH3X6P4rb9y0/z3\nNaYgTVIgiMS5VdPHNXlVc9XuTZStmE/Z6pd4euUXCgVplAJBJAEsufmsRkMhd9x1ZAwcxf7XH+Pw\nh8sVCtIoBYJIglhy81khu48sKZlu3/4F6fmDKX7lXip3fszTK7/QmII0oEAQSSCNjSkkpaaTd8F0\nkjO7ULltI4BmH0kDCgSRBNPYmEJyZg69rnqYzqdO9pUN/81r0WyaxDgFgkgCamxMISktA4CKrevZ\nt+gBDlRUKxTER4EgkqCaGmiu2vMZh9ctofStpyir9GiVVAEUCCIJrbGB5k4nTyJr5HmUrfoHh9Yv\nY/fBKsbd+2b0GygxRYEgkuBWTR9H5/TkgDIzI/ecqaT3HU7J4oep3LWJT/Yc5vuzV7RRKyUWKBBE\n2oG1MyY0DIXkFPIm/4qkjGwOr3sdgHc+LdE1Cu2YAkGknVg7YwIWVJbcMZuel/+BLuf+2FemaxTa\nLwWCSDty3yUjGpSldOqGmVFzYDdlq18CdI1Ce6VAEGlHCkfmN3qP5oP/eYX9rz/K4Q+XAzB4+qJo\nNk1igAJBpJ1p7B7NOd+4nPTegyledD/V+7ZS4XGajtrOKBBE2qG514xucI2CJafSrfBWLDWdvS/e\nRW11BbsPVmnmUTsSViCY2QQz+8jMNpnZtBDbzcwe9G5fa2Yn+217wsz2mNn6oH3uMLPtZva+93He\nsX8dEQlXqGsUUjp1pdukW6jet5UDK/4OaOZRe9JsIJhZMvAIMBEYCkwxs6FB1SYCA72PqcAsv21P\nAhMaefv7nHMjvA91WIpE2arp4+iQHDj3KKP/SPIuuJXsUd/zlWnmUfsQzhnCacAm59xm51wVMA+Y\nHFRnMvCUq7MSyDGzXgDOubeBkkg2WkQi58PfNjw57/jVM0hKy6C2qoKaA3sAzTxqD8IJhHxgq9/r\nbd6yltYJ5QZvF9MTZtYlVAUzm2pmRWZWtHfv3jDeUkRa6v4Q01Gdc+x5bgZ7/nEHtVUVABpkTnBt\nOag8CzgeGAHsBP4YqpJz7jHnXIFzriAvLy+a7RNpNwpH5jeYeWRmZI+6iOp9Wyl5/U8AGmROcOEE\nwnbgOL/XfbxlLa0TwDm32znncc7VArOp65oSkTYy95rRDQaZM/qfTPboizm87nUObXgDqBtk1nhC\nYgonEN4DBppZfzNLAy4FFgbVWQj8wDvbaBRwwDm3s6k3rR9j8LoAWN9YXRGJjlAL4WWf+V+k9xlK\nyT9nUV26C9B4QqJqNhCcczXA9cBiYCPwrHNug5lda2bXeqstAjYDm6j7a/8n9fub2TPACmCQmW0z\ns6u9m35vZuvMbC1wNvDfkfpSInL01s4InBRoScl0m/Rz0nsPCijXeELiMedcW7chbAUFBa6oqKit\nmyGS8H69YB1Pr/yi2XpjBuQy95rRUWiRHAszW+2cK2iunq5UFpEGZhYOC3m3NU/FIfYuuIuKbRsB\njSckGgWCiIS05OazQt5Yp3LXJopfvofayiOAxhMSiQJBRBoVPJ6QlJ5Jt0m3UFO2l5IlXy5IMPw3\nr0W7adIKFAgi0qTg5bI79BlC9hmXcHjDGxz+4C0Ayio9Wu8oASgQRKRJocYTss+4lLTegyh952+4\nWg9AWIPQEtsUCCLSrODxBEtKJu87v6Lnf/0OS/qyXDfViW8KBBEJS/B4Qkp2d5Izc3C1Hqr2bAbQ\nTXXinAJBRMIWahG8/W/+hV1zf+W7inn3wSqNJ8QpBYKIhK1wZH6D9Y46n/IdAIoX3U/d0mQaT4hX\nCgQRaZFV08cFvE7J7k7uOddQuXU9B1e/5CvXVNT4o0AQkRYL7jrKHDaOjAGnUvrWHKqLtwGaihqP\nFAgi0mLB908wM3In3EBqXn9qqyt85eo6ii8KBBE5KnOvGR0wFTUlK5eel99Des8TAuqp6yh+KBBE\n5Kg1WCrbjNqqCkqW/ZnqfXV31VXXUfxQIIjIMQle2sJVl3N43VKKX31AVzHHGQWCiByTmYXDArqO\nkjO70OXcqVTu+JCD/3nZV66uo9inQBCRYxbcdZQ59Cwyji+g9O2nfBesqeso9ikQRCQi/LuOzIzc\n8T8FS2L/kj/5ytV1FNsUCCISEcFdRymd8+j27Vvocs41AfXUdRS7FAgiEjHBXUcdTzid1Nx8nHPU\nVlcC6jqKZQoEEYmo4FlHAMWv3Mu+F+/GOQeo6yhWKRBEJKJmFg7DgsrSegyg/NP3OPzBm74yLZMd\nexQIIhJx9wWtddTplG+T1nsQ+5fOxnO4FKhbJnvBmu1t0TxphAJBRCKuwVpHScl0m3gTtVVHKFk6\n21d+0/z326J50ggFgoi0irnXjA7oOkrtdhzZoy6m4vM1eA7v95V/f/aK6DdOQlIgiEirCe46yh51\nEb1/NIvkzC6+snc+LYl2s6QRCgQRaTWFI/MZ2D3T99pSUknumI2r9VCxdb2vvP+0V9qieRJEgSAi\nrWrJzWc1KCtb9Ry7n7mNyl2bAHBo1lEsUCCISKsLvjah08jzSO6YTclrD/lWRNWso7anQBCRVhe8\nrEVShyy6nPtjqnZ/ysGiF33lmnXUthQIIhIVDZa1GDSGjBNOo3T5XN+KqKC1jtqSAkFEoqbBiqjj\nriO1W19c5WFfudY6ajthBYKZTTCzj8xsk5lNC7HdzOxB7/a1Znay37YnzGyPma0P2ifXzJaY2Sfe\nf7sEv6+IJJZQK6L2vPxe0noMCKintY7aRrOBYGbJwCPARGAoMMXMhgZVmwgM9D6mArP8tj0JTKCh\nacBS59xAYKn3tYgkuND3YS5n/1tz8JQf9JWr6yj6wjlDOA3Y5Jzb7JyrAuYBk4PqTAaecnVWAjlm\n1gvAOfc2EOrKk8nAHO/zOUDh0XwBEYk/wbOOakp3Ufbu8+x/4wlfmbqOoi+cQMgHtvq93uYta2md\nYD2cczu9z3cBPcJoi4gkgOCuo7Tu/el82gUcXreEii1rfeXqOoqumBhUdnWLpLtQ28xsqpkVmVnR\n3r17o9wyEWktwV1H2WdMISWnF8WLH/LdTAdg8PRF0W5auxVOIGwHjvN73cdb1tI6wXbXdyt5/90T\nqpJz7jHnXIFzriAvLy+M5opIvPDvOkpKTSd3/E+p2b+TshXP+sorPE5dR1ESTiC8Bww0s/5mlgZc\nCiwMqrMQ+IF3ttEo4IBfd1BjFgJXeJ9fAbzYRF0RSUDBXUcZ/UbQ5eyryDzpmwH11HUUHc0GgnOu\nBrgeWAxsBJ51zm0ws2vN7FpvtUXAZmATMBv4Sf3+ZvYMsAIYZGbbzOxq76a7gXFm9glwrve1iLQz\nwV1HnU+7kNTcuiHI+ltugmYdRYP5H/BYV1BQ4IqKitq6GSISYQvWbA9YtqK2qoJ9L99DRr8RdDp5\nkq/8/ktGUDiyufkqEszMVjvnCpqrFxODyiLSvhWOzKdHpzTfa0tNx1VXsv+tOdSU7fOVa62j1qVA\nEJGYsGr6ON9zMyN3/E+htpaS1/8UUG/cvW9GuWXthwJBRGLG/X53WEvN6Un2mVMo/2QlRz7+t6/8\nkz2HQ+0qEaBAEJGYUTgyP2DWUeeCQlK796f0nXkBA8z9dIe1VqFAEJGY4j/ryJJTyJs8jR6XzsTM\nAurpDmuRp0AQkZjjf8Faam4+yRmdcbWegAFm3WEt8hQIIhJzZhYOIyXwhIC9L97Nnmf/B+ep9pVp\n1lFkKRBEJCZtuuv8gNdZw8ZRXfwFB1Y9F1CuC9YiR4EgIjFrzIBc3/OOJ5xGx8Ff58C/51FdvM1X\nXlbpUddRhCgQRCRmzb1mNP49R7nnTCUpJZ3ixQ8HzDpS11FkKBBEJKZ9dveXXUfJWV3IOfsqPIdL\n8RzeH1DvhFs1FfVYKRBEJOb5zzrKGj6O3j98iJSs3IA6NQ4tk32MFAgiEvP8Zx2ZJWEpqdRWHuHQ\n+qUB9bRM9rFRIIhIXAiedXTw/dcofuU+yj8NXAFZs46OngJBROKGf9dR51O+TWrX4yj+5/9SW1Xu\nK9eso6OnQBCRuBHQdZSSSu6E6/GU7aF0+dyAepp1dHQUCCISV/y7jjr0OZGsERM4WLSQyl2bAuqp\n66jlFAgiEnf8u466jL2SjoPGkJTeMaCOuo5aToEgInFnZuEw3wVrSR2yyJv8K1K79G5QT11HLaNA\nEJG45H/BGkDNwWL2vvg7qvfvDCjXMtnhUyCISNzyX+sI5yjfvJriVx/AuVpf8e6DVW3QsvikQBCR\nuDX3mtG+5ymdu5F7zo+o3LqeQ2sWBdTTHdbCo0AQkbj2uV/XUeawcXTofzL733yS6tJdAfU066h5\nCgQRiXsDu2cCYGZ0nXA9mFH61pyAOmWVnrZoWlxRIIhI3Fty81m+5ymdu5N34a/J/dZ1Deqp66hp\nCgQRSQj+XUcZX/ma7z7MnvKDAfUGT18UvKt4KRBEJGHUdx0BOOfYPf929i38fcDNdCo8ThesNUKB\nICIJw7/ryMzIHHwmFZ+v4dDafwbU0wVroSkQRCSh+HcdZY2YQIevDGf/sj9TU7YnoJ7GExpSIIhI\nwvly1lESuRNuBKB40f0BF6yB7rAWTIEgIgnHv+soNacnXb55DZ5D+/EcLg2opzusBVIgiEhCCug6\nGj6OXlc+0OA+zAD91XXko0AQkYRVv9aRmWEpadRWlVP6zjO4mmpfHYe6juqFFQhmNsHMPjKzTWY2\nLcR2M7MHvdvXmtnJze1rZneY2XYze9/7OC8yX0lEpI7/WkcAlds+4MDyuZQufzqgXF1HdZoNBDNL\nBh4BJgJDgSlmNjSo2kRgoPcxFZgV5r73OedGeB+6WkREIi7ggrXjTyHraxMoW/U8FVvXB9TTBWvh\nnSGcBmxyzm12zlUB84DJQXUmA0+5OiuBHDPrFea+IiKtyv+CtS7fvJqULj3Z9/K91FYe9pXrgrXw\nAiEf2Or3epu3LJw6ze17g7eL6Qkz6xJ2q0VEWsB/1lFSWgbdJt2C5+A+9r/xREC99n7BWlsOKs8C\njgdGADuBP4aqZGZTzazIzIr27t0bzfaJSALx7zpK7z2IrhN/RudRFzWo156XyQ4nELYDx/m97uMt\nC6dOo/s653Y75zyu7kqR2dR1LzXgnHvMOVfgnCvIy8sLo7kiIqH16JTme5417BxSc3rinMNz5ICv\nvD0vkx1OILwHDDSz/maWBlwKLAyqsxD4gXe20SjggHNuZ1P7escY6l0ArEdEpBWtmj6uQVnJaw+x\n+5nbqK2u8JW112Utmg0E51wNcD2wGNgIPOuc22Bm15rZtd5qi4DNwCbq/tr/SVP7evf5vZmtM7O1\nwNnAf0dejBz2AAAJfUlEQVTua4mIhObfdQTQcdAYqvdtYf/S2QHl7TEUzH9Z2FhXUFDgioqK2roZ\nIhLnhv/mtYCuof1vPknZqn/Q7Tu/JHPIN3zlYwbkNriWIR6Z2WrnXEFz9XSlsoi0O2tnTAh4nfP1\ny0jvPZji1x6iuuTLIdJ3Pi2JdtPalAJBRNol/64jS06h2+RfktIpD8+RsoB67anrSIEgIu1W5/Rk\n3/OUzt3pdfXDdOgzpEG99hIKCgQRabeCu47MknCultLlczmw6vmAbe0hFBQIItKuBc86AqN631ZK\n33qS8s8Dr1xO9FBQIIhIuxcwnmBG1/N+RmrXPuxb+HtqDuwOqHvCrYkbCgoEERG+vHcC1K13lHfB\ndKj1sOcfMwIWwatxiXv/BAWCiAh1905IsS9fp+bm063wVmoO7KZy+4cBdRP1/gkKBBERr013BY4n\nZPQbQf6PHyfj+FMa1E3E8QQFgoiIn+BB5uTMHAAOf7icstUvBWxLtFBQIIiIBAkOBeccRz78F/tf\nf4zDH70TsK3ftFc4/bdLotm8VqNAEBEJocHMo/P/m/T8wex76Q8NpqPuPliVEGcLCgQRkUYEzDxK\n7UDe935Dam4f9j4/k8odHzWoH++hoEAQEWlE8EqnyR2y6H7x/yG5Yzbln/0n5D7xHAoKBBGRJgSP\nJ6Rk5dLrygfIGTMFqBtfCBavoaBAEBFpRnAoJHXIAqBq7xZ2Pnkj1cXbGuwTj6GgQBARCcPnd58f\ncE/mep5D+9n9zK1U79vaYFu8hYICQUQkTKumjws4W0jL+wo9ptyJw7Fr3q1U7f28wT7xFAoKBBGR\nFgoIhW596TnlLsyS2DX3V1Tt+axB/XgJBQWCiMhR8A+F1K7H0fPye8gcNIbU3PyQ9eMhFBQIIiJH\nyT8UUjp3p+vEG7GUNDwVhzi0bmmDGUixHgoKBBGRY9DwBjtwsOhFihfdR/GrD+JqqgO2xXIoWKg5\ntLGqoKDAFRUVtXUzREQa8P9F71wtB/41lwMr5pPWaxB5F9xGSqeuIfcLFSiRZmarnXMFzdXTGYKI\nSAR8fvf5vvspmCWR843L6VZ4K9X7trBzzs+o3PlxyP36TXslZs4aUtq6ASIiiaL+fgr1v+DrB5n3\nL51NSufuTe5bv080zhgaozMEEZEIu2xUX9/ztLx+9Lj0tyRn5uBqPRQvfpiq3Zsb3bctzxgUCCIi\nETazcFjA7Tjr1ezfwZGPV7Jzzk3sf+MJaqsrGn2PtggFDSqLiLSSwdMXUeEJ/B3rKT9I6Zt/4dDa\nf5KS3YMuZ19NxldHYxYiQbyOtRsp3EFlBYKISCtasGY7N81/v0F5xdb1lCz+X5yrpffVj2BJyc2+\n19EGgwJBRCTG/HrBOp5e+YXvtav14DlYTEp2d2orj1D86oN0OmUS6X1ObPSM4WhCQdNORURizMzC\nYXx+9/m+QWdLSiYlu272UfW+L6jYuo7df5vG7qd/weEPl+M81Q3eozXHFnSGICLSRoJ/uddWV3J4\n3escePd5PAd2k9Qxm94/mkVyRueAepeN6svMwmFhf466jERE4kRwMLhaDxWfraFi2wa6jL2iQf1k\nMz6967yw3z+iXUZmNsHMPjKzTWY2LcR2M7MHvdvXmtnJze1rZrlmtsTMPvH+2yXcLycikkg+v/t8\nxgzI9b22pGQyBhSEDAMATyv9Id9sIJhZMvAIMBEYCkwxs6FB1SYCA72PqcCsMPadBix1zg0Elnpf\ni4i0S3OvGR32gHFyE1NUj0U4ZwinAZucc5udc1XAPGByUJ3JwFOuzkogx8x6NbPvZGCO9/kcoPAY\nv4uISNz7/O7zmw2GKacf1yqfHU4g5AP+Nwvd5i0Lp05T+/Zwzu30Pt8F9AizzSIiCa+xUGjpgHJL\nxMTids45Z2YhO8XMbCp13VD07ds3VBURkYQU7YXuwjlD2A74n5/08ZaFU6epfXd7u5Xw/rsn1Ic7\n5x5zzhU45wry8vLCaK6IiByNcALhPWCgmfU3szTgUmBhUJ2FwA+8s41GAQe83UFN7bsQqB9CvwJ4\n8Ri/i4iIHINmu4ycczVmdj2wGEgGnnDObTCza73b/wQsAs4DNgFHgB82ta/3re8GnjWzq4EtwMUR\n/WYiItIiujBNRCTBaS0jERFpkbg6QzCzvUApcCBoU3YYZd2Afa3XugZCtak19w+nflN1GtsWzrEN\nVabjHbnjHao8VL1oHvNjPd4tfY/WOt6NbUu0n/GvOOean5XjnIurB/DY0ZQBRW3dztbcP5z6TdVp\nbJuOd9sf70aOb6j/BlE75sd6vFv6Hq11vJs4lgn3Mx7OIx67jF46hrJoOtbPb+n+4dRvqk5j23S8\nj75+pI53qPJ4P94tfY/WOt6NbUvEn/FmxVWX0bEwsyIXxqCKRIaOd/TpmEdXIh7veDxDOFqPtXUD\n2hkd7+jTMY+uhDve7eYMQUREmtaezhBERKQJCgQREQEUCCIi4qVAAMzsLDP7l5n9yczOauv2tAdm\nlmlmRWY2qa3bkujMbIj3Z/sfZnZdW7enPTCzQjObbWbzzexbbd2ecMV9IJjZE2a2x8zWB5U3eR/o\nIA44BHSg7iY+0ogIHW+AXwHPtk4rE0ckjrdzbqNz7lrqFpAc05rtTQQROuYLnHPXANcCl7RmeyMp\n7mcZmdk3qPtl/pRz7iRvWTLwMTCOul/w7wFTqFtx9a6gt7gK2OecqzWzHsC9zrnvR6v98SZCx/tr\nQFfqAnifc+7l6LQ+/kTieDvn9pjZd4DrgL865/4WrfbHo0gdc+9+fwTmOuf+E6XmH5OYuGPasXDO\nvW1m/YKKffdyBjCzecBk59xdQFNdFPuB9NZoZ6KIxPH2dstlAkOBcjNb5Jyrbc12x6tI/Xw75xYC\nC83sFUCB0IQI/YwbdUv8vxovYQAJEAiNCHUv59Mbq2xmFwLjgRzg4dZtWkJq0fF2zk0HMLMr8Z6d\ntWrrEk9Lf77PAi6k7o+dRa3assTVomMO3ACcC2Sb2Qmu7r4xMS9RA6FFnHPPA8+3dTvaG+fck23d\nhvbAOfcm8GYbN6Ndcc49CDzY1u1oqbgfVG5EOPeBlsjR8Y4uHe/oaxfHPFEDIZz7QEvk6HhHl453\n9LWLYx73gWBmzwArgEFmts3MrnbO1QD193LeCDzrvryXsxwDHe/o0vGOvvZ8zON+2qmIiERG3J8h\niIhIZCgQREQEUCCIiIiXAkFERAAFgoiIeCkQREQEUCCIiIiXAkFERAAFgoiIeP1//KvDiUdBS7kA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11e1f6da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "paramsFC = fit_corr_model.make_params()\n",
    "paramsFC['kappa'].set(value=kappa,vary=False)\n",
    "paramsFC['sigma'].set(value=sigma,vary=False)\n",
    "paramsFC['ga'].set(value=ga,vary=False)\n",
    "paramsFC['eta'].set(value=eta,vary=True)\n",
    "\n",
    "resultFC = fit_corr_model.fit(t=t,data=autocorr,params=paramsFC)\n",
    "\n",
    "plt.figure()\n",
    "plt.semilogx(t,autocorr,\"o\")\n",
    "plt.semilogx(t,resultFC.best_fit, 'k--')\n",
    "\n",
    "print(resultFC.fit_report(min_correl=0.25))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# initial prior, A inv gamma, D gamma\n",
    "# both D and A have mean 1 and std 10\n",
    "alpha_A=0.01\n",
    "beta_A=0.01\n",
    "alpha_D=2.01\n",
    "beta_D=1.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    ga = pm.Gamma('ga', alpha=aD, beta=bD)\n",
    "    A = pm.InverseGamma('A', alpha=aA, beta=bA)\n",
    "\n",
    "    B = pm.Deterministic('B', pm.math.exp(-delta_t * ga))\n",
    "\n",
    "    path = Ornstein_Uhlenbeck('path',A=A, B=B, observed=ts_all)\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}