membrane fluctuations gamma MCMC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['seaborn-dark', 'seaborn-darkgrid', 'seaborn-ticks', 'fivethirtyeight', 'seaborn-whitegrid', 'classic', '_classic_test', 'fast', 'seaborn-talk', 'seaborn-dark-palette', 'seaborn-bright', 'seaborn-pastel', 'grayscale', 'seaborn-notebook', 'ggplot', 'seaborn-colorblind', 'seaborn-muted', 'seaborn', 'Solarize_Light2', 'seaborn-paper', 'bmh', 'tableau-colorblind10', 'seaborn-white', 'dark_background', 'seaborn-poster', 'seaborn-deep']\n"
     ]
    }
   ],
   "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 lmfit as lm\n",
    "\n",
    "SMALL_SIZE = 16\n",
    "MEDIUM_SIZE = 18\n",
    "BIGGER_SIZE = 20\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "SEED = 350107321 # from random.org\n",
    "np.random.seed(SEED)\n",
    "\n",
    "print(plt.style.available)\n",
    "plt.style.use('seaborn-white')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we describe the system in a continuous way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "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": 5,
   "metadata": {},
   "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": 6,
   "metadata": {},
   "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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "time=np.logspace(-5,0,50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "acf=[quad(integrant,qmin,np.inf,args=(t))[0] for t in time]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc75211c828>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAECCAYAAADkaECYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU5d3/8fdkspGQhCWBEAKBBLhldwEVELRCiy2CaOsKuNDqo1V8XGofrUvV1uKCWq2i1Yq11W7QYvWRopJHEdkUVBSFmy1hB9kSwhKyze+PmeE3jJNkwCRnls/rurgmued7Tr7xeM0nZ7nPcXk8HkRERIIlON2AiIhEJgWEiIiEpIAQEZGQFBAiIhJSotMNNBVjTAowGNgO1DrcjohItHADnYCPrbVHAt+ImYDAGw4LnG5CRCRKDQc+DByIpYDYDvDaa6+Rm5vrdC8iIlFhx44dTJgwAXyfoYFiKSBqAXJzc8nPz3e6FxGRaPONQ/M6SS0iIiEpIEREJCQFhIiIhKSAEBGRkBQQIiISkgJCRERCiqXLXE/Ye6u/ZtH63SS5E0hyJ5CcmECS23XM9ymJCSS7E0hJSiDZ7SYlKYHURDepSQmkJrl9/xJoleQm0a3cFZHop4AA5q/ZxT+Wbaa6to7q2m//fIxkdwJpKW7SkxNplewmPdlNWnIi6SmJZKYmkpGaSOvURDJSk8hITSQzNYm2acm0SUuiTZr367RkNy6Xqwl+OxGRExN2QBhjrgV+DuQDnwG3WWsXN1DfD3gKOAPYCzwLPGqt9QTUDAemAf2BrcBUa+2MoPWc4qsZApQB/wTutNYeDLf3xtw/ri/3j+sLgMfjoabO4w2LGg9HamupqqnjSE0dVb5/R2rqOFJTS2V1HZXVtRyuruVItff7w9W1HKqq5VBVzTGvB4/UsLXsMKsrq6morKGispq6BrIo2Z1Am7QkslunkJORQnbrFLIzksnxfd8pqxWdslLJzUolSXssItIMwgoIY8yVwPPAg8DHwBTgbWPMQGttSYj6DsA8YCVwCXAq8BDemXrTfDW9gbnAm8Avge8BLxlj9ltrZ/lqegDz8d5j6QKgEHgYaA1cc2K/csNcLtfRw0skAyQ1x4/B4/FwqKqWisoa9ldWs+9gFWWHqyk/VM2+Q96v9x6oYveBI+w+cIS1OyvYfaCKqtq6oH6hgy8wOrdpRZd2aXTPTqOgfTrds9PpkJGiPREROSGNBoQxxoU3GF6w1j7gG3sXsMCtwM0hFrvRt+5x1tpDwBzf3VbvMsY8Za2tBu4ESoHLfXsVc40xOcB9wCzfen4JlAAXWGtrfD87EbjZGJPkW09UcrlcpKd4DzvlZqWGtYzH42H/4Rq+rqhke3kl28oOs628ku1lh9leXsmq7ft556sdxxwma5XkpqB9Gj07ZnBSbga9O2VwUm4mnbJSFRwi0qBw9iB6AAXAG/4Ba221MeYt4Lx6lhkFFPvCwe914B68d11d5Kt5NfCQk69mojEmD9iBd6/hHn84+H72s3gPV8Udl8tFVloSWWlJ9OyYEbKmpraO7eWVlOw+yMY9BynZfYjSPQf5dNM+3lyx7WhdRmoivXMz6ds5k1O6tuXUrm3o3KaVQkNEjgonIHr5XtcFjW8Aiowxbmtt8E2eegHvh6gH6GWMWQHk1bNO//KpQAawwxjzZ+BCvIeo/gzcHnzfcvFKdCfQpV0aXdqlATnHvLe/spo1OypYvaOC1Tv2s2p7BX/7aDMvLywFvIeqTu3allML2nBaQTsG5mfpiiyROBZOQGT6XiuCxivwzqNIB/aHWCZUvf+9htbpr/F/uj0F/Afv3kR/vOcy3MANYfQuATJTkxjUrR2DurU7OlZTW8fqHRV8smkfn2zcxyebypj75Q7Au5cxrCibEb1yGNErm/y2aU61LiIOCCcg/Mccgq+58Y/X8U2uEPV+dWGu0392eJW1drLv62LfOYhHjDEPWGt3NNa8NCzRnUC/zln065zFlUO6AbCr4ggflezlgzW7+GDtrqOBUZiTztm9cvh+v04MKmhLQoIOR4nEsnACotz3mgHsDBhvjfeDPNTlpuW++kAZAe/tDxoLXKe/5oDv67lBNe8CjwF98Z6nkCaWk5HCmAGdGDOgEx6Ph/W7DjB/zW4+WLOLvyzdxMsLS8nNTGXMgE6MHZjHwPwsnbsQiUHhBMRa32shx54zKARs0EnmwGUKg8b831tr7QFjzPYGatYAh/DuYSQH1fj3LL79jDZplMvlokeHDHp0yODHZ3XnwJEailft5M0V2/nz4o289GEJXdq1Ykz/PC4ZlE9hTuvGVyoiUSGcM5Brgc3AeP+AMSYJGAMU17NMMTDKGJMeMDYe2IN3kp2/Zqwxxh1Us9Jau9NaWwEsAX5ojAnscwxwBFgeRu/SxFqnJHLByZ35w1WD+PieUTz2owEUZrfmxQUbOPfx+Ux6aSnzvtpJbUOzAEUkKjS6B2Gt9RhjHgaeMcbsAxYCNwHZwJMAxpgiIMdau8S32HS8k+nmGGMeAwYCd+GdAV3lq5mGd9LdTGPMi3gve52Id2Kd393AO8A/jDG/B04BfgE8Ya0tRxyV1SqJiwd14eJBXdhVcYS/fbSJ15Zu4id/WkZ+21ZMPLOASwd1oW168E6giESDsK5htNZOB+4AJuGdxNYGGG2t9V+Wei+wOKB+O94P/ERf/XXA3dbaaQE1K4CxeA8rzfZ9fY21dmZAzXvAD4AueGdcT8E7ae/uE/hdpRnlZKQwZWRPFvzPd3huwqnkt23Fw/9ZzZlTi7l79hdsLz/sdIsicpxcHk9sHAowxnQDSoqLi8nPz3e6HQHsjgr+uKiUWcs343K5mHBGV244p4gOGeHNHBeR5rdlyxZGjhwJ0N1aWxr4nmZBSbMxuRlMvag/7/3sHC48uTN/WryREY++x9T/rGLvwarGVyAijlJASLPLb5vGIz8aQPFtZ/P9fp144YMNDH/k/3hq3loqq4Mn4YtIpFBASIvplp3Ok5eezDu3jGBErxyenLeG7z+1gEXrdjvdmoiEoICQFtezYwbPTTyNP//4dOo8Hq74w1Ju+8dn7Dmg22uJRBIFhDhmeM8c3r5lBDd9pwdvrtjGyCfm849lm4mVCydEop0CQhyVmuTmZ6MNb908nJ4dWvPzWZ9z2QtL2Famy2JFnKaAkIjQq2MGf79uCA9f1J+VW8sZ8/QC3rNfO92WSFxTQEjESEhwcdnpXXlzyll0zEzlmpc/5pG5q6mpDXXDYBFpbgoIiTiFOa15/cZhXH56V557fz1XvLiUHeWVTrclEncUEBKRUpPcTL2oP7+99GRWbivnB08vYP6aXU63JRJXFBAS0caf0pk3bjqLnNYpXP3yR/xhwYbGFxKRJqGAkIjXo4P3kNPoPrn8+q1VTJ2zijrdTlyk2SkgJCq0Snbz7IRTmXRmAb//YAM/m7mCap28FmlW4TxRTiQiuBNcPHhBXzpmpjDtnTXsPljFcxNOJT1F/xuLNAftQUhUcblc3HRuTx75YX8+XLuLy19cwm7dokOkWSggJCpdOrgrL0waxJqdFfzouUVs3nvI6ZZEYo4CQqLWqD4dee0nZ7LvUDUT/rCUnfs1V0KkKSkgJKqdVtCWVyafzp4DR5j4h6Xs04OIRJqMAkKi3sld2vDiVYPYuPcQV738ERWV1U63JBITFBASE4YWZTP9ilP5att+fvLKMj2pTqQJKCAkZozq05HHLxnIR6V7ufG1TzRPQuRbUkBITLng5M786oJ+FK/+mtv+sYJazbgWOWGaYSQxZ+KZBVRU1vDI3NW0T0/m/nF9nW5JJCopICQm3XBOEXsOHOEPH5bQu1MGlw7u6nRLIlFHh5gkZt35/ZMY3jObe1//kk827XO6HZGoo4CQmJXoTuB3l59Cx6wUbnh1OV9XaCKdyPFQQEhMa5OWzAuTBrH/cA03vPoJVTW6skkkXAoIiXm9O2Xy2MUDWL5xH/e/+aXT7YhEDZ2klrhw/oA8Vm7dz/Pz19MvL4srztBJa5HGaA9C4sYdow0jeuXwyzdWsnzjXqfbEYl4CgiJG+4EF7+77BTy2rTip699ohv7iTRCASFxJSstiWevOJW9B6u45/WVeDyaaS1SHwWExJ1+nbO4ZVQv3vpiO//+bJvT7YhELAWExKXrzy7itIK23PvvlWwrO+x0OyIRSQEhccmd4OKJSwZSW+fhZzNXUKeb+ol8gwJC4lZB+3TuO78Pi9bv4eVFpU63IxJxFBAS1y4d3IVRvTvwyNzVrNlZ4XQ7IhFFASFxzeVyMfWiAWSkJHLr3z/TrThEAiggJO7lZKTwm4v68+W2/TxVvMbpdkQihgJCBBjdN5dLBuXz3Pvr+XxLmdPtiEQEBYSIzz3n96Fdegr3vr5SjyoVQQEhclRmahL3nt+bFVvK+etHm5xuR8RxCgiRAOMG5jGksD2Pzl3N7gNHnG5HxFEKCJEALpeLX43vy+HqWqbOWe10OyKOUkCIBOnRIYNrhxfyz0+2sHTDHqfbEXGMAkIkhCnn9qRzm1bc+++VVNdqboTEJwWESAitkt3cP64va3YeYMaHJU63I+IIBYRIPb7bpyOjenfkt/PW6o6vEpcUECIN+OXYPnjw8OCbXzndikiLU0CINKBLuzSmnNuTuV/uYP6aXU63I9KiFBAijbh2eCFd26Uxdc4qzbCWuJIYbqEx5lrg50A+8Blwm7V2cQP1/YCngDOAvcCzwKPWWk9AzXBgGtAf2ApMtdbOqGd9LqAYSLDWnhNu3yLfVnJiAj8bbbj5r5/y+qdb+eFp+U63JNIiwtqDMMZcCTwPvAr8ECgD3jbGdK+nvgMwD/AAlwAvAA8BtwfU9AbmAiXARcCbwEvGmB/V08a1wHfC6VekqZ3fvxP9O2fx+DuWyupap9sRaRGNBoTvL/cHgRestQ9Ya+cA44DdwK31LHYj3r2TcdbaOdbaXwNTgbuMMUm+mjuBUuBya+1ca+1teAPovhA95AOPAnrCvDgiIcHFXT84iW3llbyip89JnAhnD6IHUAC84R+w1lYDbwHn1bPMKKDYWnsoYOx1oB0wOKDmfwMPOflq+htj8oLW9xwwG/g4jH5FmsXQomzOMTk8+946yg5VOd2OSLMLJyB6+V7XBY1vAIqMMe56lglVD9DLGJMO5DVU4x8wxkwATifg8JSIU/7nvJOoOFLD9PfXO92KSLMLJyAyfa/BD+yt8C2fXs8yoer97zW0zqM/0xiTg/dE983W2r1h9CrSrHp3yuSiU/L546JStuw71PgCIlEsnIBw+V6Dr+/zj4e6UY0rRL1f3XGs83fAYmvt38PoU6RF3P497w7uE+/o8aQS28IJiHLfa0bQeGu8H+QH61kmuD4j4L39DawToNwYcwEwBphijEk0xiTiDRCX73sXIg7Ia9OKa4Z1Y/ZnW/lq2/7GFxCJUuEExFrfa2HQeCFgg04yBy4Tqh7fMgeA7Q3UrAEuxBsYJUC17984YITv67PD6F2kWfz07B5kpibx8Fw9M0JiV7gBsRkY7x/wXao6Bu/EtVCKgVG+k9F+44E9eCfZ+WvGBp3kHg+stNbuBO7He8VT4L/5wCe+r5eH0btIs8hKS2LKuT34YM0uPly72+l2RJpFozOprbUeY8zDwDPGmH3AQuAmIBt4EsAYUwTkWGuX+BabDkwB5hhjHgMGAncBd1pr/dcHTsN72epMY8yLeC97nYh3Yh3W2lK88ySOMsaU+d5bdoK/r0iTmTSkgJcXlvL4u5ZhPdrjcumop8SWsGZSW2unA3cAk4BZQBtgtLXWf1nqvcDigPrteD/wE3311wF3W2unBdSsAMbiPaw02/f1Ndbamd/ydxJpESmJbm44p4hPN5Xx4TrtRUjscXk8sXHzMWNMN6CkuLiY/HzdK0daxpGaWs557H06t2nFzOuHaC9Cos6WLVsYOXIkQHffkZujdDdXkW/BvxexbOM+Fq3X86sltiggRL6lSwZ1oWNmCk8Vr228WCSKKCBEvqXUJDc3nF3ERyV7Way9CIkhCgiRJnDZ6V3pkJHCU8WaXS2xQwEh0gRSk9z819lFLNmwl6UbtBchsUEBIdJEJpzRlezWKTz9fzoXIbFBASHSRFKT3Fx/diEL1+1hWaluPizRTwEh0oSuOKMr7dOTdUWTxAQFhEgTSktO5LoRhSxYu5vlG/c53Y7It6KAEGlik4YU0C49md/pXIREOQWESBNLS05k8rBuvG93sWq7nhch0UsBIdIMJp5ZQFqymxc/2NB4sUiEUkCINIM2aclcOrgLb6zYxrayw063I3JCFBAizeTHZ3XHA8z4sMTpVkROiAJCpJnkt01j7IBO/PWjTZQfqna6HZHjpoAQaUbXjSjiYFUtry7d6HQrIsdNASHSjPrkZTKiVw4vLyylsrrW6XZEjosCQqSZXT+ikN0HjjD7061OtyJyXBQQIs1sSFF7+nXO5MUPNlBXFxuP+JX4oIAQaWYul4v/GlHEht0HeXfVTqfbEQmbAkKkBXy/Xy5d2rXi+fnr8Xi0FyHRQQEh0gIS3QlcO7yQTzeVsUw38ZMooYAQaSEXn9aFtmlJ/H7+eqdbEQmLAkKkhbRKdnPlkG7MW/U1G3YdcLodkUYpIERa0MQzC0h2J/DHRaVOtyLSKAWESAvKyUhh3Ml5zFy2RbffkIingBBpYZOHdedwdS1/+3iT062INEgBIdLC+uRlMqSwPa8sKqWmts7pdkTqpYAQccDks7qzrbySt7/UxDmJXAoIEQece1IHCtqn8dKHeuKcRC4FhIgD3Akurh7ajU82lfHpJk2ck8ikgBBxyMWDupCRksjLC0udbkUkJAWEiENapyRy6eAuzPliO9vL9dxqiTwKCBEHXTW0G3UeD39arCfOSeRRQIg4qEu7NEb3zeUvSzdxuEpPnJPIooAQcdjks7pTfriaf326xelWRI6hgBBx2KCCtvTvnMWMD0v0xDmJKAoIEYe5XC4mn9WN9bsO8sHaXU63I3KUAkIkAozpn0dORgqv6C6vEkEUECIRIDkxgQlndOU9u0vPipCIoYAQiRATziggye3SJa8SMRQQIhEiJyOFsQPymLlsMxWVelaEOE8BIRJBrh7WjYNVtcxcpktexXkKCJEIMiC/DacVtOWVxaW65FUcp4AQiTBXD+3Gxj2HeH/N1063InFOASESYc7rl0tuZqru8iqOU0CIRJgkdwKThhSwYO1u1u6scLodiWMKCJEIdNngLiQnJvDK4lKnW5E4poAQiUDtW6dwwcA8/rl8K+WHdcmrOEMBIRKhrh7WjcPVtfzj481OtyJxSgEhEqH65mVxevd2vLK4lFpd8ioOSAy30BhzLfBzIB/4DLjNWru4gfp+wFPAGcBe4FngUWutJ6BmODAN6A9sBaZaa2cErecq4DagB7ANeA34jbW2KtzeRaLVNUO7ccNrnzBv1U5G9811uh2JM2HtQRhjrgSeB14FfgiUAW8bY7rXU98BmAd4gEuAF4CHgNsDanoDc4ES4CLgTeAlY8yPAmquBl4G3gEuAJ4Dfgb89jh+R5Go9d0+HencphUzPixxuhWJQ40GhDHGBTwIvGCtfcBaOwcYB+wGbq1nsRvx7p2Ms9bOsdb+GpgK3GWMSfLV3AmUApdba+daa2/DG0D3BaznDuA1a+0d1tp51tongAeA640xrY/3lxWJNonuBK4e2o2lJXtZubXc6XYkzoSzB9EDKADe8A9Ya6uBt4Dz6llmFFBsrT0UMPY60A4YHFDzv4GHnHw1/Y0xecaYBOBt4E9B67aAy9eTSMy7ZHAX0pLdmjgnLS6cgOjle10XNL4BKDLGuOtZJlQ9QC9jTDqQ11CNtbbOWnubtfbdoJqxQCXevQ+RmJfVKomLT8vnzRXb+Lqi0ul2JI6EExCZvtfgKZ0VvuXT61kmVL3/vYbWGfgzj2GMGQ1cA0y31h5suG2R2HH1sO5U19Xx6pJNTrcicSScgHD5XoOvs/OP19WzTH3X5dWdyDqNMecC/wSWAHc30K9IzOmenc7Ikzrw2pKNVFbXOt2OxIlwAsJ/ZiwjaLw13g/yUH/Jl4eozwh4b38D6wz8mQAYYy4F5gCfA+dba7WfLXFn8rDu7DlYxRufbXO6FYkT4QTEWt9rYdB4IWCDTjIHLhOqHt8yB4DtDdSs8Q8YY64H/gp8AHzXWlsWRs8iMWdIUXtOys1gxsISPB5NnJPmF25AbAbG+wd8l6qOAYrrWaYYGOU7Ge03HtiDd5Kdv2Zs0Enu8cBKa+1O388ZD0wH/oV3z0HnHSRuuVwuJp/VndU7Kli8fo/T7UgcaHQmtbXWY4x5GHjGGLMPWAjcBGQDTwIYY4qAHGvtEt9i04EpwBxjzGPAQOAu4M6AGdDTgI+BmcaYF/Fe9joR78Q6jDGpeCfnbQeeBk41xgS29oUCQ+LNuIF5PPKf1cxYWMLQHtlOtyMxLqyZ1Nba6XgnrU0CZgFtgNHWWv9lqfcCiwPqt+P9wE/01V8H3G2tnRZQswLvJauFwGzf19dYa2f6Ss4EOuK9HHa+b/2B/3of/68rEt1Sk9xMOLOA4tVfU7Jbfx9J83LFyrFMY0w3oKS4uJj8/Hyn2xFpNl9XVHLWw+9x+eldeOCCfk63I1Fuy5YtjBw5EqC7tbY08D3dzVUkynTISGXswDxmLt+iZ0VIs1JAiESha4Z141BVLX/7SBPnpPkoIESiUL/OWQwpbM+MhSUcqdHEOWkeCgiRKHX9OUXs3H+Ef3+qiXPSPBQQIlFqRM9s+nTK5PkP1lOnJ85JM1BAiEQpl8vF9ecUsWHXQd5dtdPpdiQGKSBEotgP+uXStV0az72/XrffkCangBCJYonuBK4dUchnm8tYWrLX6XYkxiggRKLcxaflk906mefeX+90KxJjFBAiUS41yc01w7ozf80uvtq2v/EFRMKkgBCJARPPKCA92c3vP9BehDQdBYRIDMhKS2LCmQW8uWIbm/cecrodiREKCJEYMXlYd9wJLl5csKHxYpEwKCBEYkRuVioXnZLP3z/ezO4DR5xuR2KAAkIkhlx3diFVtXX8cWGp061IDFBAiMSQopzWnNc3l1cWlVJ2qKrxBUQaoIAQiTG3jOrFgaoanp+vcxHy7SggRGKMyc1g/Mmd+eOiEr7eX+l0OxLFFBAiMeiWUT2pqfXwzHvrnG5FopgCQiQGFbRP59LBXfjrR5s0L0JOmAJCJEZNObcnCS4Xv5231ulWJEopIERiVG5WKlcN7cbsT7ewdmeF0+1IFFJAiMSw688uIi05kSfeXeN0KxKFFBAiMaxdejI/Gd6d/6zcwedbypxuR6KMAkIkxv34rO60TUti2jvai5Djo4AQiXEZqUn89JwefLBmF0s27HG6HYkiCgiRODBpSAEdM1N47G2rZ1dL2BQQInEgNcnNraN6sXzjPv71yVan25EooYAQiROXDOrCaQVt+fVbX7H3oG7kJ41TQIjEiYQEF7+5sD8VlTX8Zs4qp9uRKKCAEIkjJjeD60YUMmv5Fhat3+10OxLhFBAicebmkT3p2i6Ne2avpLK61ul2JIIpIETiTGqSm4cu7MeG3QeZ/v56p9uRCKaAEIlDw3vmMP7kPJ57fx3rvj7gdDsSoRQQInHqnvP7kJacyC9mf0FdneZGyDcpIETiVHbrFH7xg5P4qGQvs5ZvcbodiUAKCJE4dvFpXTi9WzsemrOKnXo8qQRRQIjEsYQEF7+5qD/VtXVc/+pyXdUkx1BAiMS5Hh1a88QlA/l0Uxn3vL5S92qSoxQQIsJ5/Trx3yN7Mmv5FmYsLHW6HYkQCggRAeC/R/ZkdN+OPPTWV3y4VrOsRQEhIj4JCS6euORkenbI4Ma/fELp7oNOtyQOU0CIyFHpKYm8eOUgXC649k/LqKisdrolcZACQkSO0bV9GtOvOJUNuw9y699XaBJdHFNAiMg3DO2RzX3n92Heqp3cMetzqmrqnG5JHJDodAMiEpmuHFLAvkNV/HbeWrbsO8TvJ51Gm7Rkp9uSFqQ9CBEJyeVyccuoXvz20pP5dFMZF05fRIlOXMcVBYSINGj8KZ157dozKD9czYXTF7Jkwx6nW5IWooAQkUYN7taO2T8dSvv0ZCa9tJSZyzY73ZK0AAWEiISloH06/7phGIO7teOOWZ9z17++0A3+YpwCQkTClpWWxCuTT2fysO7MXLaZEY++x2/mrGLvwSqnW5NmEPZVTMaYa4GfA/nAZ8Bt1trFDdT3A54CzgD2As8Cj1prPQE1w4FpQH9gKzDVWjsjaD3jgV8BPYA1wN3W2v8Nt28RaVpJ7gTuG9uHq4YW8NS8tby4YAN/WbqJyWd15yfDu5OZmuR0i9JEwtqDMMZcCTwPvAr8ECgD3jbGdK+nvgMwD/AAlwAvAA8BtwfU9AbmAiXARcCbwEvGmB8F1JwLzALeBy4EPgdmG2POPJ5fUkSaXkH7dJ649GTeuWUEI3pl83TxWoY/8h6Pv2NZvnEvNbWaOxHtGt2DMMa4gAeBF6y1D/jG3gUscCtwc4jFbvSte5y19hAwxxiTAtxljHnKWlsN3AmUApf79irmGmNygPvwhgLAL4F3rbVTfN/PNcYUAL8Axp3ILywiTatnxwymTziNlVvLeeLdNTzz3jp+93/ryEhJZEhRe4b3zOasnjl0a5+Gy+Vyul05DuEcYuoBFABv+AestdXGmLeA8+pZZhRQ7AsHv9eBe4DBwCJfzauBh5x8NRONMXnAPmAo3wygfwO/Msa4rbV6uolIhOjXOYsZVw9m38EqFq3fw4frdrFg7W7e+WonAB0yUujSLo28Nq3Iy0qlU1Yqndq0olNWKukpibRKctMqyU1qkpuUxAQSEhQmTgsnIHr5XtcFjW8Aiur5oO6F97BQcD1AL2PMCiCvnnX6l9/l6y9UTSugC949EBGJIG3TkxkzoBNjBnTC4/Gwcc8hFqzbzWebythWdpgvtpTx9peVjd6+IyUxgSR3Ai4XuBNcJLhcJLi8E/iCs8PF/x+obyclluOmd6dMXrp6cJOvN5yAyPS9VgSNV+A9h5EO7A+xTKh6/3sNrdNfcySMGhGJYC6Xi27Z6XTLTmfSmQVHxz0eD3sOVpZfgIUAAAV5SURBVLG9rJId+ys5VFVDZXUth6tqOVxdx+HqWiqra6mt81Bb58Hj8VDngTqPhzqPh8CH3h3zNaFvLBjrD8nr2bF1s6w3nIDwB2/wf2L/eKg/A1wh6v3qwlznifxcEYkCLpeL7NYpZLdOoT9ZTrcj9QjnKqZy32tG0HhrvB/SoW7OUh6iPiPgvf1BY4Hr9Nc09HMD+xIRkWYQTkCs9b0WBo0XAjboJHPgMqHq8S1zANjeQM0avOca6uqpOQBsC6N3ERE5QeEGxGZgvH/AGJMEjAGK61mmGBhljEkPGBsP7ME7yc5fM9YY4w6qWWmt3WmtPYz3aqfxHOsC4H1dwSQi0rwaPQdhrfUYYx4GnjHG7AMWAjcB2cCTAMaYIiDHWrvEt9h0YAre+Q+PAQOBu4A7rbX+OfnTgI+BmcaYF/Fe9joR78Q6v6nAW8aYF4DZwBXAEGDEif/KIiISjrBmUltrpwN3AJPwTmJrA4y21vovS70XWBxQvx3vB36ir/46vLfImBZQswIYi/eQ0Wzf19dYa2cG1Mzx/cxzfDUDgPEN3eJDRESahssTI9d/GWO6ASXFxcXk5+c73Y6ISFTYsmULI0eOBOhurS0NfE93cxURkZBi6ZnUboAdO3Y43YeISNQI+Mx0B78XSwHRCWDChAlO9yEiEo06AesDB2IpID4GhuOdX6FLYEVEwuPGGw4fB78RMyepRUSkaekktYiIhKSAEBGRkGLpHIQ0M2PM40BHa+1Ep3uRpmWMuQ34Md67Jy8D/stae6ThpSQaGGMuBu4HUoA/+58MGg7tQUhYjDGjgSud7kOanjHmDOAa4HSgP96TlqEeJSxRxhiTCzwOfAfoDZxjjPleuMsrIKRRxpgOwK+Ah5zuRZrFPuAma+1B392ZPwe6OtyTNI3v4r256dfW2mrgz8Bl4S6sQ0zSIGOMC5gB3A4UOdyONANr7Rq8t9jHGNMJ+G/gaid7kibTmWMfjbANCPteRAoIAcAYcxnwTNDwTuAl4DNr7QLfXXslStW3ja21fX3vFwD/AZ6z1s5r6f6kWbg59qmcLo5jnpjmQUiDjDEfAa3w/k/VDu8T/f5qrb3R0cakSRljTgbeBB611v7O6X6kaRhjrgTOtdZe7ft+MjDUWvuTcJZXQEjYjDFXA6N0FVNsMcbkAF/gPQ8xy+l+pOkYY/KAJXgvQNgDvAM8ba2dHc7yOsQUY4wx44DXrLUZQePXAj/He/zxM+A2PVcjOjXDNr4VyATuMcbc4xt7y1p7dxO2Ld/CiW5za+02Y8ztwDwgGfh3uOEA2oOIKcaYocBcIMFa2zpg/ErgZeBBvPdbmQIMAwZaa0uc6FVOjLZx/HFym2sPIgYYY1LwXnnyK+Ag3r8U/O+58P4P9IJ/gowx5l3A4v3LUde7RwFt4/gTCdtc8yBiw/fxPvP7DiD4BGMPoAB4wz/gux76LeC8lmpQvjVt4/jj+DZXQMSGj/E+LvBpjr2kDaCX73Vd0PgGoMgY842HhEhE0jaOP45vcx1iigHW2q0NvJ3pe60IGq/A+wdCOrC/OfqSpqNtHH8iYZtrDyL2uXyvwX+B+MfrWrAXaR7axvGnRba5AiL2lfteM4LGW+P9n+hgy7YjzUDbOP60yDZXQMS+tb7XwqDxQsD6bs4m0U3bOP60yDZXQMS+tcBmYLx/wBiTBIwBip1qSpqUtnH8aZFtrpPUMc5a6zHGPAw8Y4zZBywEbgKygScdbU6ahLZx/Gmpba49iDhgrZ2O91rqScAsoA0w2lq7wdHGpMloG8efltjmutWGiIiEpD0IEREJSQEhIiIhKSBERCQkBYSIiISkgBARkZAUECIiEpICQkREQlJAiIhISAoIEREJ6f8B02TwAvuvERQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_list=[np.sqrt(l**2+m**2) for l in range(0,3) for m in range(0,3)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n",
      "[1.         1.         1.41421356 2.         2.         2.23606798\n",
      " 2.23606798 2.82842712]\n"
     ]
    }
   ],
   "source": [
    "q_list=np.array(sorted(q_list)[1:])\n",
    "print(q_list.size)\n",
    "print(q_list)\n",
    "q_list=q_list*qmin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "acf_d = np.array([corrfct(t,q_list) for t in time])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc728425d30>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAECCAYAAADkaECYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUZfrG8e9kMumNQGhJKAF5QXpHFBusZRWk2BF7r6urrq6ru7q7VnZdGyr+7Lo2WCwrohIElSJNmsKLdAIBpIYkhIRkfn+ciRvCJJlgkkkm9+e65ppwzvOePNlxc+e097i8Xi8iIiLlhQW7ARERqZ8UECIi4pcCQkRE/FJAiIiIX+HBbqCmGGMigf5ANlAc5HZERBoKN9AKWGCtPVh2RcgEBE44fBPsJkREGqghwLdlF4RSQGQDvP3227Rs2TLYvYiINAjbtm1j7Nix4PsdWlYoBUQxQMuWLUlLSwt2LyIiDc0Rh+Z1klpERPxSQIiIiF8KCBER8UsBISIifikgRETELwWEiIj4FUqXuR61pTPeJ9/OoCQsAm+Yh5Iwzy/vuD2UuCPBHYHXHYHXHQnhkeCOAo/zcnlicHmiwRONOyKGcE8EHreLcHcYEe4wPOEuPL6vI8PDiAgPIzLcTUR4GO4wV7B/fBERvxQQQMHKz+m54xPCOUSE69fP0nHQG04+UeQRxQFvJHlEstsbRR7R5BDNfm8Mub73PFcMB8LiyA9PpCA8gQJPIoURSeCJJSYynGiPm9jIcGIinHfn327iozzER4X/8p7g+zox2kOUx10D/6uISGMXcEAYY64B7gbSgCXAHdbauZXUdwOeAgYCu4HngMettd4yNUOA8UB3YAvwiLX2lXLb6e2rOQ7YC0wG7rHW5gXae1UG3vTy//7h9ULJISguxHuokENFBRQXHcRbdIDiooOUFB6kpKjA98qnpPAA3sJ8vEUH8BYV4C3Kx1WYD4V5uIvySCjKJ7EoH/ehXNxFuYQXZeMpysVzKJcwSv73fUuAQt8rD4rwsD8snt2uJHZ6E/nZm8i24ni2FSew3tuEbG8y2TRlm7cJh8p9jNEeN01iPCTGRNAkxkOTmAiaxkXQLC6SlPhIUuIiaRb/v68jwnWkUUSOFFBAGGMuBV4AHgIWALcAnxtjelpr1/upbw5MB1YA5wN9gL/j3Kk33lfTBZgGfAL8GTgNeNkYk2OtneSr6QjMwplj6RwgA3gUiAOuOLofuQouF7idQ0uuiFg8gKc2vo/XC4V5cDAHCvZB/m44sMf32o3nwB6S83aSnPczHXN3QO4ayNsBYYWHbwYXhdHNyYtqwf7Iluz0tCLb3ZpNtGJtcXM2Hkxg1bYcduUVsje/yO+PmxIXSeukaFKTommVGOV83SSajGaxtGkaQ2S49khEGqMqA8IY48IJhonW2gd9y74ELHA7cKufYTf5tj3CWpsPTPXNtnqvMeYpa20RcA+wAbjIt1cxzRiTAjwATPJt58/AeuAca+0h3/cOB241xnh822mYXC6IjHNeCa0DG+P1QsFe2L8dcrbAvixcOVuI3LeFyJwskveto+2Or6CkzP8snhhIzoC2hkMpx7IvoRM7ojuytSSZHbmFbM8pYOveA2zdW8DKbTlkrtpOQdH/9mzCXNA6KZr2zWLJaBZL+2axdG6VQJdWCSRG10p0ikg9EcgeREegLfBx6QJrbZEx5lPgjArGDAMyfeFQ6kPgTzizrs7x1bxV9pCTr+YSY0xrYBvOXsOfSsPB972fwzlc1fi4XBDdxHk17+y/pvgQ5GTBrrWwe53z2rUWshYQvmIyTYGmQJfIRGjRFVr1BNMP0gdAYjpeYE9+EZt257NhZx7ry7wmL95C7sFfPgpSk6Lp0iqeY32B0btNE1omRtXF/xIiUgcCCYhOvvc15ZavAzoYY9zW2vJndjsBM/3UA3QyxiwFWlewzdLxUUA8sM0Y8yYwCucQ1ZvA78vPWy4+7nBo0s55MfTwdQX7YMdK2P6D77UCFr8O3z3vrI9riSu9P8lpA0huM4hePfo42/Pxer3s2H+Qldk5rMzez4/ZOazMzmHGqh2U+GI+NSmafu2a0K9dMv3aNqFTi3hdqSXSQAUSEAm+9/3llu/HuY8iFsjxM8Zffem6yrZZWpPi+/op4DOcvYnuOOcy3MANAfQuZUUlQptBzqtU8SEnKLIWwOb5kDUfVn7irItMhIwTocNQ6DgUV1IbWiRE0SIhipNN8182UVBUzKpt+1m8cQ8LN+5m7tpdfLRkKwDxkeEMzGjKqZ2bc0rnFFolRtflTywiv0IgAVH655+3guUlHMnlp75USYDbLD3AvdJae6Xv60zfOYjHjDEPWmu3VdW8VMEdDq17Oa8B1zjLcnfAxtmwJhPWzvhfYDQ9BjoOg2NHQPogCHOuforyuOmVnkSv9CSuPKE9Xq+XzbsPsHDjbhZs2M3Xq3cyfeV2ADq3jOeUzs05tXNzeqcnEe7WFVQi9VUgAbHP9x4PbC+zPA7nF7m/y033+erLii+zLqfcsrLbLK3J9X09rVzNl8ATQFec8xRS0+KaQ9dRzsvrhZ2rfWGRCYtedQ5Jxbd21ncbA6l9nPMjPi6XizZNY2jTNIbRfdLwer38tCOXr1bt4Cu7g5e+XsfzM9fSNDaC4T1bM7pPKt1TE3G5dChKpD4JJCB+8r1ncPg5gwzAljvJXHZMRrllpf+21tpcY0x2JTWrgXycPYyIcjWlexYV7aFITXK5IMU4r+NuhIP7wU6DH/4DC16Cec9BUlsnLHqPg2Yd/WzCRacW8XRqEc91J3Ugp6CIb1bv5NPlW/n3d5t4bc4GOqTEMrpPGiN7p5KapMNQIvVBIPv3PwGbgZGlC4wxHuAsILOCMZnAMGNMbJllI4FdODfZldYMN8a4y9WssNZut9buB+YBY4wxZfs8CzgILAqgd6lpkfHQ4zy46B248yc4ZwI0OwbmPAPP9oU3R4H9DEoqviM9IcrDWT1aMWFsXxbcN4yHR3UnOTaCJz63nPDYDC5+aR6ZK7dTUqK/AUSCqco9CGut1xjzKPCsMWYPMBu4GWgGPAlgjOkApFhr5/mGTcC5mW6qMeYJoCdwL84d0KV3eo3HuenuA2PMSziXvV6Cc2NdqfuAL4D3jTEvAr2BPwL/tNbuQ4IrOgl6j3VeuTtg0euw8BV450JIagP9roI+l0JMcoWbSIzxcPHANlw8sA2bduXz4ZItvDt/E1e9vpAOKbFcMySDkb1TNX2ISBC4vN7A/kozxvweuA0nGJbgXGo617fuNeAya62rTH0/nCuQ+uKcu5hgrX2s3DZPBx4DOgObgIetta+Vq/kN8DeckPkZ547uR6y1JeXq2gHrMzMz9UzqYCouAjsV5r8EG76B8CjodTEMuRMSUwPaRFFxCVOXZ/PirHX8mJ1Ds7gILjuuHZcMakuT2PJHHEXk18jKymLo0KEA7a21G8quCzgg6jsFRD20/UeY/yJ8/za4wqDflXDC7RDfIqDhXq+XuWt38eLX65i1+meiPW6uGdKe607qQGyk5pkUqQmVBYSuMZTa0+JYGP4U3LoYepwP8yfCUz3hywcgb1eVw10uF4M7NuP1Kwfw+e9O5NQuzXl6xhpOGT+TDxZu1jkKkVqmgJDal9QGznkWbl7g3EMx+2l4qgfMfAyKCgLahGkZz3MX92HyDcfRKimauyYtY/iz3zJvXdVBIyJHRwEhdadpBxg9EW6cBx1OhZkPw/ODYd2sgDfRt20yU24YzFMX9mJPXiEXTpzHdW8uZPPu/KoHi0i1KCCk7jXvDBe8CeOmgLcE3hgBU66HvJ0BDQ8Lc3FOr1Rm3Hkyd57WiW9+2smZT33DpEVZhMo5NZH6QAEhwdPhVLhxrnOF0/JJ8Gw/+P4t5+7tAER53Nx86jF8/rsTObZVAnd+sJQb317M7rzCqgeLSJUUEBJcnmgYej9c/w2kdIaPboLXzoZ9WQFvIj05hneuHcQ9Z3Zm+srtnP6vr5lpd9Ri0yKNgwJC6ofmXeDyqTD8acheAi8MgZ++DHi4O8zF9Sd14MObjqdJjIfLX13A/R+u4EDhr3/GuEhjpYCQ+iMsDPpeBtfOcp6y9/a5MP0vzpTkAeraOpGPbz6Bq05oz5vzNjJqwmyy9x2ovZ5FQpgCQuqfZh3h6unQ93L49kl4fTjkbA14eJTHzf1nH8trV/Qna88BRj43mxVbNDOLSHUpIKR+8kQ7N9mNfgmyl8ILJ8Ca6dXaxMmmOZNvGEx4WBjnvTCX6T9ur3qQiPxCASH1W4/z4dqZENcC3joX5lbvceSmZTxTbhrMMS3iuObNhbw6e32ttCkSihQQUv+ldIKrM6HL2fD5H+GL+6HE34MM/WseH8V71x7Hace24MFPfuTPH63gUHHg40UaKwWENAwRMXDe69D/apjzNHx4gzNzbICiI9w8P7Yv156YwetzN3L9W4s4eEhXOIlURgEhDUeYG347Hk79Eyx7F/59ARzMrXpc6fAwF3/8bRceOqcr01fu4LZ3lmhPQqQSCghpWFwuOPEuGPEMrPsKXj8bcn+u1iYuPa4dD5x9LNN+2MadHyzVrLAiFVBASMPU51K48N+wYxW8chrs2VCt4Vee0J67Tjd8uGQr9324QnM4ifihgJCGy5wJl30M+bvhjXMgJ7taw286pSM3ndKBd+Zv4q//XamQEClHASENW/oAuOQ/zkywb450wqIa7jzNcPngdrwyez3//HJ1LTUp0jApIKThS+sLF70Du9fDW6OhICfgoS6XiwfOPpYL+qXzzIw1TJi5phYbFWlYFBASGtqfCOe/DtuWwzsXQVHg8y+Fhbl4eHR3RvRszePTLJ8sDXxaD5FQpoCQ0GHOhFEvwsbZ8MHl1bpPwh3mYvx5Penbtgl/mLwMu21/7fUp0kAoICS0dD8XzvoHrJ4GU66DksBvhosID2PC2D7ERoZz3ZsL2Xcg8IARCUUKCAk9/a+CYX+BFZNh2r3VGtoiIYoJY/uQtecAd7y3RPdISKOmgJDQdMLtcNzNMP9FWPxGtYb2b5fMA8OPJXPVDp6ZoZPW0ngpICR0DXsQMk6BT38PmxdUa+i4QW0Z3SeVf2WuZsYqTRMujZMCQkKXOxzOfQXiW8H742B/4L/oXS4XD4/qTpeWCdz27hI27MyrxUZF6icFhIS2mGRnSo6CfU5IHCoMeGiUx82L4/riDnNx3ZuLyC8M/NGnIqFAASGhr2U3OOc52PwdfHZ3tYamJ8fwzEW9Wb1jP49+tqqWGhSpnxQQ0jh0Gw3H/w4WvQoLX63W0CHHpHD54Ha8MXcj89btqqUGReofBYQ0HkMfgA5DYepdsOm7ag2963RDm+QY/jB5GQcK9aAhaRwUENJ4hLnh3JchMQ0+uKxaE/vFRITz2JgebNyVzxOf21psUqT+UEBI4xLdxJmzKW8n/Pd2qMYU38d1aMq4QW15dc56Fm6o3qyxIg2RAkIan1Y94ZR74ccPYfkH1Rp6z5mdSU2K5u5Jyygo0qEmCW0KCGmcjv8dpA+ET++EfVkBD4uNdA41rduZp+dHSMhTQEjjFOaGUS9AySH48AYoKQl46PEdm3HRgDb83zfr+H7TnlpsUiS4FBDSeCVnwBmPwPqv4bsXqjX0j7/tTMuEKO7SoSYJYQoIadz6XAqdzoTpf4EdKwMeFh/l4ZExPVizI5cJM9fWXn8iQaSAkMbN5YIRT0NkPPzn2mpNxXFSpxSG92zNxK/XsnVv4E+wE2koFBAicc1h+FOwbRnMerRaQ+8+3VDihfG6N0JCkAJCBKDL2dD7Evj2SdiyOOBh6ckxXHVCe/7z/RaWZe2txQZF6p4CQqTU6Q9DTDPn+RHVeFTpjSd3oGlsBH/7dCXeatx4J1LfKSBESkUlOiGxdTEsei3gYfFRHu44rRPz1+/m8x/0cCEJHQoIkbK6nwvthkDmg5D7c8DDLuiXzjHN43jks5UUHgr8ngqR+kwBIVKWywVn/QMK8+HLBwIeFu4O476zurBxVz5vzN1Qa+2J1CUFhEh5KQYG3wJL/w0bZgc87GTTnBM7pfB05k/syQv8clmR+koBIeLPiXdBYhvnhHVxUcDD7vttF3IPHuKpzJ9qsTmRuqGAEPEnIgbOfAx+XgnzJgQ8zLSM54L+bXhr3kbW/Zxbiw2K1D4FhEhFOv8WzG9h5qPVmvH1jt90IjI8jMem6RnW0rApIEQqc8ajzkOFpt0T8JCU+EiuHpLB5z9sZ2V2Ti02J1K7FBAilWnSFk66C1Z+AmumBzzsiuPbERvh5rmv1tRicyK1SwEhUpXjboEm7eDLPwd8h3VSTASXDm7Hp8uzWbND5yKkYQoPtNAYcw1wN5AGLAHusNbOraS+G/AUMBDYDTwHPG6t9ZapGQKMB7oDW4BHrLWvVLA9F5AJhFlrTw60b5FfLTwCTr0fJl8Fy96HXhcFNOzqE9rz6uz1TJi5hn+e36uWmxSpeQHtQRhjLgVeAN4CxgB7gc+NMe0rqG8OTAe8wPnARODvwO/L1HQBpgHrgdHAJ8DLxphzK2jjGuCUQPoVqXFdR0OrXjDjb1BUENCQpnGRjB3Ylo+WbGXTrvxablCk5lUZEL6/3B8CJlprH7TWTgVGADuB2ysYdhPO3skIa+1Ua+3fgEeAe40xHl/NPcAG4CJr7TRr7R04AXTE7avGmDTgcWBrdX44kRoTFgan/RVysmD+iwEPu/bEDNxhLp6fpXMR0vAEsgfREWgLfFy6wFpbBHwKnFHBmGFAprW27J9NHwLJQP8yNf8te8jJV9PdGNO63PaeB6YACwLoV6R2tD8ROv4GvvkH5O8OaEiLhCgu6JfOpEVZbNFDhaSBCSQgOvney/8JtA7oYIxxVzDGXz1AJ2NMLNC6sprSBcaYscAAyhyeEgmaYX+Bghz49p8BD7n+5A54vTBxlh5NKg1LIAGR4HvfX275ft/42ArG+KsvXVfZNn/5nsaYFJwT3bdaawP7k02kNrXsBj0vgu8mwt5NAQ1JTYpmTJ803lmwmR05gZ2/EKkPAgkIl++9/JNQSpf7m9vY5ae+VEk1tvkMMNda+14AfYrUjVPvc95n/D3gITee0oFDxSW89M26qotF6olAAmKf7z2+3PI4nF/keRWMKV8fX2ZdTrllZbcJsM8Ycw5wFnCLMSbcGBOOEyAu379diARDYhoMuh6WvQfblgc0pG3TWM7plcpb8zaxWzO9SgMRSECUTkuZUW55BmDLnWQuO8ZfPb4xuUB2JTWrgVE4gbEeKPK9RgAn+r4+KYDeRWrHCbc7T6D78s8BD7nplA4UHCrm5W+1FyENQ6ABsRkYWbrAd6nqWTg3rvmTCQzznYwuNRLYhXOTXWnN8HInuUcCK6y124G/4FzxVPY1C1js+3pRAL2L1I7oJs6U4GszYe1XAQ3p2Dye33ZrxRtzNpJ78FAtNyjy61V5J7W11muMeRR41hizB5gN3Aw0A54EMMZ0AFKstfN8wyYAtwBTjTFPAD2Be4F7rLWl+9fjcS5b/cAY8xLOZa+X4NxYh7V2A859Er8wxuz1rVt4lD+vSM0ZcA189wJ89XfIONl5Gl0Vrjkxg0+XZzN5URaXDW5X2x2K/CoB3UltrZ0A3AWMAyYBScDp1trSfeX7gbll6rNxfuGH++qvBe6z1o4vU7MUGI5zWGmK7+srrLUf/MqfSaRuhEc6h5qyFsC6wPYieqUn0Ss9idfmbKCkpKLrOETqB5fXGxr/kRpj2gHrMzMzSUtLC3Y70lgcOghP94bEdLhyWkB7ER8t2cJt7y7h1Sv6c4ppXgdNilQsKyuLoUOHArT3Hbn5hWZzFfk1SvciNs+D9bMCGnJmt1Y0j4/ktdkbarc3kV9JASHya/UeB/GtYNbjAZVHhIdxyaC2zFr9s6YCl3pNASHya3minL2IjbNh/TcBDbl4YBsi3GG8MXdDrbYm8msoIERqQp/LIK4lzHosoPJmcZEM79maSYuyyCkoquXmRI6OAkKkJnii4PjbYMM3sGF2QEOuOL4d+YXFvL9gcy03J3J0FBAiNaXfFRDbPOC9iG6pifRv14Q35m6kWJe8Sj2kgBCpKZ5oZy9i/SzYNK/qeuDywe3ZtDufGat21HJzItWngBCpSf2ugJhmMPPRgMpP69qCVolRvDZnfS03JlJ9CgiRmhQRC8ff6txZvXl+leUedxjjjmvL7DW7WL29/ONRRIJLASFS0/pfDTFNA74v4sL+bYgMD+NV3Tgn9YwCQqSmRcTCoBtgzZewbUWV5cmxEYzslcqU77PYm69nRUj9oYAQqQ39rgJPLMx5JqDySwe3paCohCnfb6nlxkQCp4AQqQ0xydDnUlgxCfZlVVnetXUiPdISeXf+ZkJlAk1p+BQQIrXluBvB64V5zwdUfmH/Ntjt+1myeW8tNyYSGAWESG1JagPdxsCi1+DAnirLR/RqTUyEm3fn685qqR8UECK16fhboTAXFr5SZWlcZDhn92jFJ8u26pGkUi8oIERqU8vu0GEozHsBigqqLL9wQBvyC4v5eMnWOmhOpHIKCJHadvxtkLcDlr1bZWnv9CRMi3jeW7CpDhoTqZwCQqS2tT8RWvV0LnktKam01OVyceGAdJZm7ePHrTl11KCIfwoIkdrmcjl7EbvWgJ1aZfmo3qlEhIfxrvYiJMgUECJ1ocs5kNQWZv/LufS1EkkxEZzZrSVTvt/CgcLiOmpQ5EgKCJG64A6HwbdA1oKApgK/sH8b9hccYury7DpoTsQ/BYRIXek1FqKTYfZTVZYOykimfbNY3tPT5iSIFBAidSUiBgZcC6s/g51rKi11uVxc0D+d+Rt2s2ZHbh01KHI4BYRIXep/Fbgj4LsXqiwd0yeN8DCXLnmVoFFAiNSluObQ/TxY8naV02+kxEfym2NbMHnxFg4e0slqqXsKCJG6NugGKMqHxW9UWXpB/3R25xUy/Uc9s1rqngJCpK617A7thsB3E6G48jmXhhyTQqvEKCYt0slqqXsKCJFgGHQj5GTBqk8qLXOHuRjVO5Wvf9rJjpyq53ISqUkKCJFg6HQ6NGkPcydUWTqmbxrFJV4+XKKnzUndUkCIBEOYGwZeD1nzIWthpaUdUuLo3SaJyYu26GlzUqcUECLB0nssRCYE9MS5c/umYbfvZ8UWTeAndUcBIRIskfHOc6t//BD2VX746OwerYkID2Py4qqfby1SUxQQIsE04FrwlsCClyotS4z2cNqxLfhoyRYKD1U+ZbhITVFAiARTk7bQ+WxY+CoU5ldaOqZvGnvyi5ixSvdESN1QQIgE26AboWBvlU+cG9KxGc3jI5m0SIeZpG4oIESCrc0gaNXLOVldyRPnwt1hjOqdyky7g525B+uwQWmsFBAiweZyOXsRO1fD2hmVlo7pm8ahEi8fLdlaR81JY6aAEKkPuo6CuBYw/8VKyzq1iKdHWiKTdZhJ6oACQqQ+CI+AflfCT19U+ayIc/um8WN2Dj9s3VdHzUljpYAQqS/6XQlhHpg/sdKy4T1a43G7mLxIU29I7VJAiNQXcc2h2xjnWREFFd8x3SQ2gmFdnHsiiop1T4TUHgWESH0y8DoozHVCohJj+qSxK6+QmfbnOmpMGiMFhEh9ktoH0gfCdy9WesnrSSaFZnEROlkttUoBIVLfDLwO9qyHNV9WWOJxhzGiZyqZq7azJ6+wDpuTxkQBIVLfdBkB8a2rnOV1TN9Uioq9fLJM90RI7VBAiNQ3bg/0vwrWfQU7VlVY1rV1Ip1bxjN5sa5mktqhgBCpj/peDu7IKi95PbdvGks372XNjty66UsaFQWESH0U2wy6nwdL34EDeyssO6dXKu4wl54TIbVCASFSXw28Dory4fs3KyxJiY/kpE4pTFm8heISPY5UapYCQqS+atUD2h7vHGYqKa6wbHSfVLblFDBn7c46bE4ag/BAC40x1wB3A2nAEuAOa+3cSuq7AU8BA4HdwHPA49Zab5maIcB4oDuwBXjEWvtKue1cBtwBdAS2Am8DD1trdW2fhL6B18H7l4L9DLqc7bdkWJcWJESFM3lRFkOOSanjBiWUBbQHYYy5FHgBeAsYA+wFPjfGtK+gvjkwHfAC5wMTgb8Dvy9T0wWYBqwHRgOfAC8bY84tU3M58CrwBXAO8DxwJ/CvavyMIg2XOQsS21R6yWuUx83ZPVsz7Ydt5B48VIfNSairMiCMMS7gIWCitfZBa+1UYASwE7i9gmE34eydjLDWTrXW/g14BLjXGOPx1dwDbAAustZOs9begRNAD5TZzl3A29bau6y10621/wQeBK43xsRV94cVaXDc4TDwWtj4LWQvrbBsTJ80CopKmLo8uw6bk1AXyB5ER6At8HHpAmttEfApcEYFY4YBmdbasg/Z/RBIBvqXqflv2UNOvpruxpjWxpgw4HPgjXLbtoDL15NI6Os9DjyxMO+FCkv6tEmifbNYTb0hNSqQgOjkey8/Sf06oIMxxl3BGH/1AJ2MMbFA68pqrLUl1to7rLXl5xsYDhTg7H2IhL7oJOg9FlZMgv3b/Za4XC5G907lu/W72bw732+NSHUFEhAJvvf95Zbv942PrWCMv/rSdZVts+z3PIwx5nTgCmCCtTav8rZFQsjA66G4CBa+XGHJqD6pAPxHd1ZLDQkkIFy+9/IXWZcu9zflpMtPfamSo9mmMeZUYDIwD7ivkn5FQk/TDtDpDFjwMhQV+C1JaxLDcRlN+c/3WXi9uidCfr1AAqL0uYbx5ZbH4fwi9/eX/D4/9fFl1uWUW1Z2m2W/JwDGmAuAqcAy4Gxrrf//h4iEskE3QP5OWP5BhSVj+qaxcVc+izbuqcPGJFQFEhA/+d4zyi3PAGy5k8xlx/irxzcmF8iupGZ16QJjzPXAO8DXwG+stRXPOyASytqfCC26OZe8VrCHcGa3lsREuDX1htSIQANiMzCydIHvUtWzgMwKxmQCw3wno0uNBHbh3GRXWjO83EnukcAKa+123/cZCUwA/oOz56DzDtJ4uVzOXsSOH2D9135LYiPDObNbK/67NJv8Qt0TIb9OlXdSW2u9xphHgWeNMXuA2cDNQDPgSQBjTAcgxVo7zzdsAnALMNUY8wTQE7gXuKfMHdDjgQ+3/08AAA/mSURBVAXAB8aYl3Aue70E58Y6jDFRODfnZQNPA32MMWVbW67AkEan27nw5Z+dvYiMk/yWXNA/ncmLs5i6fBvn9k2r4wYllAR0J7W1dgLOTWvjgElAEnC6tbb0stT7gbll6rNxfuGH++qvBe6z1o4vU7MU55LVDGCK7+srrLWlB1gHAS1wLoed5dt+2VeX6v+4Ig2cJ8p5VsTqabBrrd+S/u2akNEslvcWbKrj5iTUuELlagdjTDtgfWZmJmlp+qtJQtj+7fCvbs4zI377hN+SF2at5dHPVjH9jpPo2FyTDkjFsrKyGDp0KEB7a+2Gsus0m6tIQxPfwjnU9P3bFT4rYnSfVMLDXHywcHMdNyehRAEh0hANuh6K8mDx635XN4+PYmiX5kxenEXhIX+3KolUTQEh0hC16gnthjgnqw8d9FtyQf90duYWMmOV/+k5RKqigBBpqE74HezPhmXv+1194jEptEyI4r0FOswkR0cBIdJQdRgKLbvD7Keg5MjDSOHuMM7rl8as1T+zde+BIDQoDZ0CQqShcrng+N/Brp/ATvVbcl7fdEq8MEnTgMtRUECINGTHjoQm7eDbJ/1Ov9GmaQzHd2zK+ws3U1ISGpe0S91RQIg0ZO5wGHwLbFkIG2f7Lbmgfxuy9hxgztpdddycNHQKCJGGrtdYiE1x9iL8OO3YFiTFeHhXd1ZLNSkgRBo6T7TzQKE102Hb8iNWR3ncjOyVyhc/bGdPXqGfDYj4p4AQCQX9r4KIOOeKJj8u6J9OYXEJU77X0+YkcAoIkVAQ3QT6XQErJsOeDUes7tIqgZ7pSby7YJOeNicBU0CIhIpBN4LLDXOe9bv64gHprN6ey7x1u+u4MWmoFBAioSKhNfS8EL5/E3J/PmL1Ob1SSY6N4JXZ64PQnDRECgiRUHL8bc7cTN+9cMSqKI+bsQPbMH3ldjbs1LO2pGoKCJFQ0uwY6DIc5k+E/CMPJY0b1JbwMBevzdlQ971Jg6OAEAk1J98LB/f7vaKpeUIUw3u05v2Fm9l3oCgIzUlDooAQCTUtjoUe58N3L8L+bUesvvKE9uQXFvO+ZnmVKiggRELRyfdASRF8Pf6IVd1SExnQPpnX5mzgULEeJiQVU0CIhKLkDOhzKSx6ze99EVed0J4tew/wxY96mJBUTAEhEqpOvAvC3DDzsSNWDevSgjbJMbz8rS55lYopIERCVUJrGHANLHsXdqw6bJU7zMXlg9uxaOMelmzeG6QGpb5TQIiEsuNvB08sfPX3I1ad3z+d+MhwXtFehFRAASESymKbwuCbYeXHsGXxYaviIsO5oH86U5dnk71PjySVIykgRELdoBshOhlm/O2IVZcNbkeJ18sbczcGoTGp7xQQIqEuKgGG3AFrM2HDt4etSk+O4fSuLfn3d5vILzwUpAalvlJAiDQG/a+G+FaQ+dARz66+ekgG+w4U8ersDcHpTeotBYRIY+CJdqbg2PwdLH33sFV92zZhWJfmvDBzLbv1xDkpQwEh0lj0HgfpA+HzP0LersNW/eGMzuQVHuLZGWuC1JzURwoIkcYiLAzO/hcczIEv7z9s1TEt4jmvbzpvztvA5t35QWpQ6hsFhEhj0uJYGHwrLHkb1n992Krbf9MJd5iL8V/YIDUn9Y0CQqSxOeluaNIO/ns7FBX8srhlYhRXHt+ej5ZsZcWWfcHrT+oNBYRIY+OJhrOfhF1r4Nt/Hrbq+pM70CTGwyOfrcRb7monaXwUECKNUYdTofv58M0/4efVvyxOiPJwy6nHMHvNLr7+aWcQG5T6QAEh0lid/jBExMJ/fwcl/3suxNhBbUhPjubRz1ZRUqK9iMZMASHSWMWlwGl/hY2znZPWPpHhbu48zbAyO4cPl2wJYoMSbAoIkcas1yXQZjB88SfIyf5l8fAerememsg/vlhNQVFxEBuUYFJAiDRmYWEw/CkoLoL3x/1yVVNYmIt7zuzMlr0HmPCVbp5rrBQQIo1dSicY9QJkLYBP7/hlrqbjOzZjdJ9UnvlqDXPW6IR1Y6SAEBE4dgScdI9zLmLe878s/us53choFsut7y5hx/6CSjYgoUgBISKOk/4Anc+GL+6DtV8BEBsZzoSxfck9WMTt7y2hWFc1NSoKCBFxhIXBqBchpTN8cDnsWguAaRnPQyO6MXvNLk3m18goIETkfyLj4MJ/g8sF714MBTkAnNcvjVG9U/lX5mqdj2hEFBAicrjk9nDe67DzJ5hyHZSU4HK5+NvIbrRvFstt7y3h5/0Hg92l1AEFhIgcKeMkOONRsFPho5vgUKHvfEQfcg7ofERjoYAQEf8GXOM8hW7pv+HNUZC/m84tE3hwRFe+XbOTh6eu1FQcIU4BISL+uVxw8j0w+iXImg//Nwx2reWC/umMG9SWl79dzy3vfq87rUOYAkJEKtfjfLj0YyjYC/83FNfG2Tx0TlfuPbMzU5dnc+HEeTonEaIUECJStbbHwdXTITYF3hiJa8m/ue6kDjw/ti+rtuUw8rnZrNqWE+wupYYpIEQkMMkZcNUXTlh8dCN8chtntPHywXWDKSou4dzn5zLT7gh2l1KDFBAiErjoJnDJf2DQjfD9W/B0L7r/+AQfX9WZNskxXPnaAibMXMOBQp2XCAXhgRYaY64B7gbSgCXAHdbauZXUdwOeAgYCu4HngMettd4yNUOA8UB3YAvwiLX2lXLbGQn8FegIrAbus9b+N9C+RaSGuT1wxiPOVU4zH4M5z9Jy4Wt82P967ko8gcenWV76eh3jjmvHZce1pWlcZLA7lqMU0B6EMeZS4AXgLWAMsBf43BjTvoL65sB0wAucD0wE/g78vkxNF2AasB4YDXwCvGyMObdMzanAJGAmMApYBkwxxgyqzg8pIrUgOQNGvwg3zoOOpxIx+wn+te0y5gyYw/ktsnkucxWDH53BfVOWs35nXrC7laNQ5R6EMcYFPARMtNY+6Fv2JWCB24Fb/Qy7ybftEdbafGCqMSYSuNcY85S1tgi4B9gAXOTbq5hmjEkBHsAJBYA/A19aa2/x/XuaMaYt8EdgxNH8wCJSw5p3hvPfgK1LcH31MK2XPce9eLk7IZ5VUb34YPExXLmgK+06dmdARjO6pybSPTWRxBhPsDuXKgRyiKkj0Bb4uHSBtbbIGPMpcEYFY4YBmb5wKPUh8CegPzDHV/NW2UNOvppLjDGtgT3AYI4MoI+Avxpj3NZaHegUqS9a94Kx70P+blg/C/far+i67iu6ur8BN+zc3IQNG1LY6m3KO96mFES3JLpZW5q0akd8QhIxsfHExyUQFxdPYkIciTGRRIaH4XK5gv2TNVqBBEQn33v5aRzXAR0q+EXdCeewUPl6gE7GmKVA6wq2WTr+Z19//mqigXScPRARqU9ikqHrKOfl9cLudbDuK5plLSJxzya67tmMJ3cR4UWFkI3z8qPA6yGXcLy4KCYMLy68rjBKcOGt5Oi4F/+BUtHyULAtuiO9/vB5jW83kIBI8L3vL7d8P845jFig/AXQCRXUl66rbJulNQcDqBGR+szlgqYdnFf/q/EAHnCCI28n5GSRv2sz+ftzKDiQx8EDuRQeyKWoII/ig/mUlBRDSTFebzGUlOD1esF7CLzOCU7w4vJ6+eUwhLdxTv1R0qxT1UVHIZCAKI3d8v/Lly4vqWBMRZ9USYDbPJrvKyINgcsFcSkQl0JM697EBLsf8SuQq5j2+d7jyy2Pw/kl7e/yhH1+6uPLrMspt6zsNktrKvu+ZfsSEZFaEEhA/OR7zyi3PAOw5U4ylx3jrx7fmFycI48V1azGOddQUkFNLrA1gN5FROQoBRoQm4GRpQuMMR7gLCCzgjGZwDBjTGyZZSOBXTg32ZXWDDfGuMvVrLDWbrfWHsC52mkkhzsHmKkrmEREaleV5yCstV5jzKPAs8aYPcBs4GagGfAkgDGmA5BirZ3nGzYBuAXn/ocngJ7AvcA91tpCX814YAHwgTHmJZzLXi/BubGu1CPAp8aYicAU4GLgOODEo/+RRUQkEAHdSW2tnQDcBYzDuYktCTjdWlt6Wer9wNwy9dk4v/DDffXX4kyRMb5MzVJgOM4hoym+r6+w1n5Qpmaq73ue7KvpAYysbIoPERGpGS5viFwWZoxpB6zPzMwkLS0t2O2IiDQIWVlZDB06FKC9tXZD2XWazVVERPwKeDbXBsANsG3btmD3ISLSYJT5nekuvy6UAqIVwNixY4Pdh4hIQ9QKWFt2QSgFxAJgCM79FboEVkQkMG6ccFhQfkXInKQWEZGapZPUIiLilwJCRET8CqVzEFLLjDH/AFpYay8Jdi9Ss4wxdwBX4cyevBC4zlp7sPJR0hAYY84D/gJEAm+WPhk0ENqDkIAYY04HLg12H1LzjDEDgSuAAUB3nJOW/h4lLA2MMaYl8A/gFKALcLIx5rRAxysgpErGmObAX4G/B7sXqRV7gJuttXm+2ZmXAW2C3JPUjN/gTG66w1pbBLwJXBjoYB1ikkoZY1zAK8DvgQ5BbkdqgbV2Nc4U+xhjWgG3AZcHsyepMakc/miErUDAcxEpIAQAY8yFwLPlFm8HXgaWWGu/8c3aKw1URZ+xtbarb31b4DPgeWvt9LruT2qFm8OfyumiGveJ6T4IqZQxZj4QjfMfVTLOE/3esdbeFNTGpEYZY3oBnwCPW2ufCXY/UjOMMZcCp1prL/f9+0pgsLX26kDGKyAkYMaYy4FhuooptBhjUoDlOOchJgW7H6k5xpjWwDycCxB2AV8AT1trpwQyXoeYQowxZgTwtrU2vtzya4C7cY4/LgHu0HM1GqZa+IxvBxKAPxlj/uRb9qm19r4abFt+haP9zK21W40xvwemAxHAR4GGA2gPIqQYYwYD04Awa21cmeWXAq8CD+HMt3ILcDzQ01q7Phi9ytHRZ9z4BPMz1x5ECDDGROJcefJXIA/nL4XSdS6c/4Amlt4gY4z5ErA4fznqevcGQJ9x41MfPnPdBxEazsR55vddQPkTjB2BtsDHpQt810N/CpxRVw3Kr6bPuPEJ+meugAgNC3AeF/g0h1/SBtDJ976m3PJ1QAdjzBEPCZF6SZ9x4xP0z1yHmEKAtXZLJasTfO/7yy3fj/MHQiyQUxt9Sc3RZ9z41IfPXHsQoc/ley//F0jp8pI67EVqhz7jxqdOPnMFROjb53uPL7c8Duc/ory6bUdqgT7jxqdOPnMFROj7yfeeUW55BmB9k7NJw6bPuPGpk89cARH6fgI2AyNLFxhjPMBZQGawmpIapc+48amTz1wnqUOctdZrjHkUeNYYsweYDdwMNAOeDGpzUiP0GTc+dfWZaw+iEbDWTsC5lnocMAlIAk631q4LamNSY/QZNz518Zlrqg0REfFLexAiIuKXAkJERPxSQIiIiF8KCBER8UsBISIifikgRETELwWEiIj4pYAQERG/FBAiIuLX/wPrz1x1rpVGcAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 15,
   "metadata": {},
   "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(100000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEQCAYAAACqduMIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dedxV0/rAv2+hDJVCaaBQlkiZMnRlSrmVmduNiHIpP6TIkCmzaybDLXQNSVzzVJQplalEVCxpkOYSDSLVe35/7LNP+91n7Xk457zv+n4+7+c9Z589PHtaz1rPeoayTCaDRqPRaDR2qhVaAI1Go9EUJ1pBaDQajUaJVhAajUajUaIVhEaj0WiUaAWh0Wg0GiVbFFqAuBBC1ADaAouBTQUWR6PRaEqF6kBDYLKUcr31h0qjIDCUw4RCC6HRaDQlSntgonVBZVIQiwFGjhzJzjvvXGhZNBqNpiRYsmQJPXr0gGwbaqUyKYhNADvvvDNNmjQptCwajUZTauSZ5vUktUaj0WiUaAWh0Wg0GiVaQWg0Go1GiVYQGo1Go1GiFYRGo9FolGgFodFoNBolWkFoNAHo2rUrbdu2LbQYGk0qVKY4CI0mcUaPHl1oETSa1NAjCI1Go9Eo0QpCo9FoNEq0gtBoNBqNEq0gNBqNRqNEKwiNRqPRKNEKQqPRlAyzZ8/mjz/+KLQYVQatIDQaTUmwYcMGmjdvTvfu3QstSpVBKwiNRlMSrF9vVMN87733CixJ1UErCI0mBJlMptAiVDk2bNgAwJZbbllgSaoOWkFoNCFYs2ZNoUWocvz2228ArFu3rsCSVB20gtBoNCWBOfdgjiQ0yaMVhEajKQm++OKLQotQ5dAKQqPxyaZNm2u66zkITVVAKwhNIpSXl/Pll1/mLb/vvvu48sorCyBRdCZOnJj7/PLLLxdQEk3SzJkzR5uy0ApCkxD33nsvBx10EBMmTKiw/PLLL+fuu+8ukFThmTt3bs7NEuCrr74qoDTB6NatG2eeeWahxSgZli1bxh577MGll15aaFEKjlYQmkSYNm0aAPPmzSusIDHw73//m913353zzjsvt2z8+PEFlCgYL774IqNGjSq0GCXDihUrAPjoo48KK0gRoBWEJhHKysqAymGrHzRoEAALFizILZs7d26hxKmSTJkyJbVjmaalLbbQ9dS0gtAkwp9//gnAkiVLCixJMlQGxVdKzJkzJ7Vjbdy4EdABeaAVhCYhXnrpJWBz77uyUV5eXmgRqhSLFy9O/BiXXXYZvXr10iMIC/oKaBKlsva0dUbRdEnDxHT//fcD5Oaa9AhCjyA0CWPORWgKQ2VR0M8++2xqxzJHELoToBWEJmG0gigsn3/+eaFFKDlMc9bUqVMLLEnhCWRiEkLsDhwNNEStXDJSylviEExTObBGH1v566+/2GqrrVKWpuphOgto/LNq1apCi1A0+FYQQohuwDOA21udAbSC0CixNlZ6kjcdrHb0srIyvvvuO/baa68CSlT8/PXXX4UWoWgIMoK4GfgJ+D9gLqDuGmo0DjzyyCOFFqHKseuuu1b4Pm7cOK0gPNAKYjNBFMQuwOVSyveTEkZTubGOICrL5Gmx88knn1T47mTyKyWSNk3qHEybCTJJ/TWwW1KCaCon2267be5ztWqbH7cPPvigEOJUOeyN3dKlSwskSXwknVeqMqSHiYsgCmIA0FsIcb4QYqekBNJULn7//ffcZ6uCOP744wshTpXD7sv/1ltvFUiS+EjaM+7xxx9PdP9xMn78eFavXp3Y/oOYmOYA3wNDgaFCCNU6GSmlDr7TKLEqCE061KxZs8L36dOns3LlSurVq1cgiTRxsXz5co466ii6dOnC22+/ncgxgjTmjwN/AyYDPwAbE5FIU+no27cvhx56aMkqiAEDBuSibK1MmTKFnXfemSZNmhRAKn+0aNEib9luu+2mXTkrAebIYfTo0YkdI4iCOBZ4WErZLylhNJWTYcOGMWzYsJKsAwHOk6Jt27YFinvCXRUNnKRJIg2K+XqniZQy93nRokU0atQo9mME6dL9DsyIXQJNlaFUPWhKOWbj9NNP973uzJkzS6LxLQUZ0+CZZ57JfW7cuHEixwiiIEZgTFLr8FdNKEq1oTXTP5cifj1yPvnkE/bZZx+GDBmSrEAhWLt2bYXvuviRQRoZboOYmL4CugOzhBDvAEvJn4dINNWGEOJEYKSUslZSx9AkR6kqiKrQYzXrLUyePLnAkuRjV9BxB7KtXr06bzK/FEgjoC+IgrCmUzzfYZ3EUm0IIdplZdDZ30qUUm1oS1XuIBRzBcCkOxZ16tShQ4cOiR4jCdIY2QZREAUJkhNC1AAuxVA8v+OeC0pTxJRqAFLdunULLULiFLOCSGPu6v33Sy9BxAEHHJB4nQzfCkJK+VOSgrjQGRgEXAHsAFxeIDk0ERk+fHihRQhFVUhZXtUVRCmSRo2MINlce/pZT0r5jPdagZgM7Cal/E0IcWPM+9akSM2aNUsy/XQxNppJUYznWqpzV0mzbt26xI8RxMT0FMYcg6o7ZX2qYlUQUsqFce5PUzhKUTlAcTaaSVGMoyU9gigcQdxcjwaOyf43P/8dOAt4BcOj6Zy4BdTEwx9//EFZWRnPP/88AGPGjOGHH34osFSlQVVQELfddluhRXBEjyAKR5A5iPEuPz8nhBiFMVeQXvFYjW8+++wzAM444wy6d+9Oly5dgKrR+EWlKlyjGTOMGNgoI4gpU6bQrFkzdtxxx7jEAip663Tu3JmFC+MzKpRSYr5CEGdynPeB3WPcn0ZTFFQFBWEyatQojjnmmFDbtm3blkMPPTRmiSr6+2+55ZaxmsEuuOCC2PZVGYlTQXTBqBmh8cFvv/2WalH07bffHoA2bdpUaPA2bNjAoEGD+O2331KTpdSoSgoC4MMPPwy97ezZs2OUxMCqIKpVq1bl7oeKnXZKp+JCVC+mrTBcT7tiZHq9275eAl5NlYJjjz2WL7/8MrWHfc2aNQBMmzatwjFbtmzJ7Nmz+eWXX3jsscdiO17Xrl0TS0GcNqVoAy8vL+fjjz8utBixYC16VFZWVpL3I25WrFiRynHi8mIyudL2PUPMXk2VhS+//DLV440fr55CMnt869evj/V4lSm4rBR7rA888ACXX55cyFAmk+HWW2+lV69eiac71yOIiqxcuTK1YwVREEcnJoVPpJQ3AjcWWIxYyWQykWyqy5YtY+TIkfTv3991P9bKYg888EDo4/mlGN0lw1KKPdaZM2cmuv/p06dzww03MGbMmLy613GzfPny3Gc9goAff/wxtWPF5cWkCYmUkr322iv09meccQYffPABHTp0oHXr1o7rVa9ePfc5yZ6lSWVSEHfddVehRQhM3CNCO6bZJ43YlpNPPhmARo0a6REEFRVm0pRmia9KRNSH3RxueiXu2mIL975A3A16GqmINc4knYbBDF6zdjySZuDAgXoEwWaX9TTQCqLAHH/88ZG2NxWMVwNvL16fNOPGjUv1eJp0MRvpNMvILliwQI8gSLc+iVYQEVi5ciVvvPFGpH2YefjDsnTpUsD7RfXq6VX1l66qYfUMCoM5gkiztvV9992nRxCka77VCiICp5xyCieddBLLli0rmAxLliwBvHsV++yzj+vvaXtVpUV5eTnPPfdclW9U7PTv3z/S9ub1lFJyzz33pKYo5s6dy+zZs6t0hybNUZtWEBEwcxlFTSYWRzKyrbZyL5NhFjh3erjMVAuVjaFDh9KjRw+GDRsWeh9eyrUUlY+T27NfrM/sFVdcwSWXXBJVJF98+umnALmcYlWRNOd9gri5AiCEqIkRHKeUUko5P6pQpYL5kkTV6FtssQXLli2LFB3pNew0UwoUU2P266+/Jh4vYZrgoozy2rZt66pAN27c6Kmgi40oduwffviBo446qsKytCPxH3nkEc4444xUj1ksFKWJSQhRTwjxPLAKmA/MdfirMpgKwstDyA8LFiwIvM17770X+biF5Jdffim0CJ4MGzaMp556ynWdUkxHHUXmjz76KG9Z2h2PSZMmuf6+dOlSysrKlLKWOmma14K0bPcB3YB3MHIuJetoXQKYL9nkyZORUnLppZeG3leYm96xY8fc53333bdK22WdiHpN+vbt67lOVVMQqmsat4KoXbs2zZo145tvvgm1vWmKuvTSS5k2bVqcohWcNL2YgiiIE4EnpJQ6/WEW8yXr3LkzQCQFMX78eA444IBY5CoV0lBo1voXgwcPTuQYpa4gTEeHKMR9L9esWcO3334bevvVq1cDhFYwxUxUD7QgBDGeb4FR/lOTJWzDoCo0XkxzA5UJU3kfeeSRiR2jFBXEzz//nPscdcIaYP784pp6jMPsW6xYc1OZ1KlTJ5FjBVEQHwNHJSJFiWJvGMxeixdPPvlk3rLKlJrCbw8njRGEeU/iSpdxwgkn5C0rReVuvfZ+5G/cuDFXX3113rYmSeV+CqN8ly1bRo8ePRKQpjhQpVFJ6l0Komb7Ax8KIe4CXgSWA3lPVlX0YjLxm/9G5fXUsmXLWGTyombNmonnz7npppsS3X8Q4u7ZdurUiTfffLPCslIcQViZMGGC5zqLFi3izjvv5Nprr+WJJ55IQSqDMGmtv/jiiwQkKR7SVBBBRhDfAg2AgcBnwGy0F1OF736HtSo/5rjLNDqRRnK1efPmOf7Wvn37xI9v5YMPPoh1f+3atctbluakoR+CupxOnz7d97r9+vVLNagyzMi62O5H3KhMTMUwgrgTo76DxgG/8RAqBTFjxgzatm0bt0gF5/rrr+eWW27Jfbf2VkvR60rVYBXbCKJPnz6B1j/88MNdRxHW+2TGlbjx7bffsu+++waSwYk4FMTo0aNzNdgrA0WpILK1GDRZbr311rxlfm/Sfffdl7fss88+49xzz40qVlFgtWmbpU6LgY0bN+aN8qZOncp+++0XKdgxTa8SPwQtKFO7dm3X39etWxdof61bt+b7779HCBFoOxVhGj67gujatWtJdkacUCmI33//PZFjOb4VQohdhRBb2757/iUiZRFy/fXX5y2L8hAGnegs5gfempenmCKM7ekgJk2axIEHHsi9997rex+qHm2xKYigvW6vEZB1f2PGjPG1z7jyk1mLW9mjt52IMqLbf//9efnll0NvnwYqBZEUbt2mucAplu/zcI6erpJzEHaiNNpB86sUs+eMtUFp3Lix43ppK7mhQ4dW+G7OlXz11Ve+96HKxV/qCsJLfj/7s6eJiStf0B133JH7/L///c/XNlHuxzHHHMOpp54aevs0SFNBuJmYbga+sX0v3m5rERClwQs6sVZsdm8r1sahU6dOjuulMWHuhpnqI4gJZc2aNXnLklQQmUyGYcOG0a1bN+rVqxdqH61bt3YNGNtmm20qfN9uu+0CH8Ne5SyJ59NvrrKFCxeGPkYpxE8UhYKQUt5k+35j4tKUOFEURNAXqpg9NWbPnp37vO222zqut8MOO6QhjiPXXHMNAK+//rrvbVTKIEkFMXnyZC688ELGjRvn2/RhTy3hFU189NEVy83bn8UwE8Uffvhh6h5rJnPn+jdkbLHFFhXepbQLa4XBriD69++fWJ15ne67SKhMIwi/qcMLERxolS2MQt9jjz3yliWpIMz5nCD1Fvx4GlmxXwf7sxWmxxrFxBTV9BjkftiPddJJJ0U6dhrY4yBq1aqV2LG0goiROExMv/zyS4UeuNf6pYhpVw5zva644opIE99jx47NfQ7qnQNqr6wkFYRZx+P9999P7BjW+3D++efnKYjvvvsu8D5Vprgw8rjh1EkKotDsxyomrzsn7OeXZKp1rSBi5Keffgq9rdnItGzZkubNm/teP07SmjQ2bclhjnfPPfewYcOG0JP05eXllJeXc+WVV4baXiVzkgoizZTugwYNYuedd45ldHrnnXfSs2dPevXqFXhb69yUm6vsLrvsolwe5H7Yn6Owo9pHH320Qn6rJPnjjz8qfH/ooYeA+DzHrGgFESNRGljzobZP9jmRxAhi1qxZse9ThdtLOGPGDF89wLDXury8nG+//Za777471PaqiPckFUSrVq1i2c9pp53GOeeco/zN9KE/9NBDc6ahOLzkRowY4VlLQ8XHH3+c+7z33ns7rrd48WLlctX9UJ1PXD3vBx54gIsuuii1YDynORa/qX6CEKuCEEIUj9N7Adh5551Db2t/qL2G6MXmWhkGeyM/f/58WrVqRY0aNQJv65fy8vJIilwlW5L3Iq6X/qWXXnJsrDt06ADA8OHDcwrCOopIe67I2knaf//9A2+vuh833HBD3jJVQKHbuU6YMCGv9w4wYMAAIFjKklIhSEW5OUKIE11+PwNYFItUJUqUobl9RKB6EN3Wj4Ok7NwNGjSo8N18Ce0N9Y8//uh7n2F7uJs2bYrU4IUJlDv66KM5/PDDQx3PK8o5TtavX69UEGnWQIaKz0WLFi0q/Na9e3fP7VXvxujRo12PY2K/v+Y6s2fP5ogjjuDCCy/0PH6auOU9iwNHN1chxI6AdXzXDGgrhFCNy6phBNXVjFW6EiPKsNz+UHspG9VLkMlkIjV+w4YNS+UFME0D06ZNo1mzZrnlQXr2Xuv+8MMPPPjgg3nLy8vLmThxou/j2KlZM/8R91IQUcpeurkJx01ZWZlSQUStuR4U63tk98jq1KlTrgiUE1HSzdvfn40bN7Llllvm0sa7VafbdVfvRBKZTIbZs2f7mmd0olOnTowdO5axY8fStGnT0Pvxg9ud/xN4Dvgw+5cBrrF8t/69D5wOvKncUxUhioLYsGFDhUb/ww8/dF1fZT+NajdOqjSj/UUcMWIEsHlyzSRIQ+R1rqeffjqPPvpo3vJNmzYxdepU38exMmLECOWLnaSJafjw4Ynt286GDRuUCiLt4DHr89K/f3/H35xQzWGpOk5+FATAU089xWWXXQa4d9z8dM6GDh1KixYtInUaTE+8NNLYOL6RUsq1GGVGewPnAWXA49nv9r9zgM7A2QnLmxpr1671nXfGJIzbpMmGDRsqNPpXX321q5np0EMPzVtWrPmZ7N4V33//PZBv0gpS1MfrXJ1s96+//npoU+BZZ52lXJ6kgkiz9/7+++8XnYnJjp9OkN9sp34VRK9evXINepQyqAD/93//B1Ahw3EQrPMzdsWdxFyRa9dASjkVmAoghGgKvCylrHwzMQp69+7Niy++yKxZs2jevDl//vknHTp04JlnnlEGSwFcd911vPHGG6GON3PmzAov5c8//8yQIUO46qqrfO8j6Ahi6tSpFepgh0mxECfvvPOO73X9mJhUrF69OvYG3e/+Tj/9dJYuXeqrQI9JEAWxatWqSKPIzp07KxVE2nm/3I537LHHem4fxcT0zTffVDAVBel0+Wmg9957b2bOnBm6/rw1Ot1UEEOGDKFfv36JdCZ871FKeVNVUQ4AL774IgDvvvsuYPjff/LJJ662wzDVr0xWr16d1/MJGsEadARh9xBJ29YchbCNViaT8XQACIq9QVq3bp0yfuHll18ONP+RyWQ802RY2X777T3zNbnNaVSvXj13Llaf/mJQEGZ6/WbNmtG3b1/X7VUKwgw4tKJ6X8z8XG7rOOFHQRxxxBEA7L777r73a8V6HqaCMOfFkrhPgYyLQoi+QA9gZ0A17sxIKdXd6xJl0qRJXHTRRb7SR0Qdir/00kuRto/6gLRp0ybS9k7UrVuXX3/9NdZ9hjWnbdq0iQULFsQqy5IlSyp879OnD88++yxSSvbcc8/Q+w1q4vSDWyNWu3ZtLr/8cgAOOOCA3DVOO62L6t5avZe83jNVx0pl/lW9L3ZXda/nzOraao1PWLJkCTvuuKPj/E3Q5/e3335jyJAhFZaZ52l27JIwMQdxcx0MPIrh2bQc+EnxV+nqUZsX3U/j6zfIzYmoRT+iPiBBTB9umNlB99tvPyCZkUmUSGqVt0mUeANrSmrYnJrC9HwJS9jI2L/97W+Ov9lHGNYes31C2CTtEYQqnsiq2LzeM78u4KrG2z7x6/VOqUaKq1atomHDhlx66aV5v4VtzAcMGMDgwYMrLDMzN5jXJon7FOTN7Q18BDSSUraTUh6t+otdwgJjRhdbb6jKOwbUw9gg2B+aoJNOqgfk5JNPVu7nX//6VzDhAmD6rp94ohE2c+ONN1b43YyOjVJxLKwyzGQyyqyopikxDswGKqrCDhIX4he794w1rYWTV0zaIwhVeg7rM+xVF8LvHIR5vv/5z3/ylpl43UPVaMUMcjXN1FbCNuaqIEdzX0UxggDqA89JKeOP5y5izALtnTt3zi276KKLEjnWqFGjIm2vekCcUlknGR170EEHAZur7l188cW537bddttcD75OnTqhjxFlBKEiznTJprtwmDQTVsKmgrDb0a3stttujr85PRPmO5AWKnNkkOfVr4Iwn4Wtt84VzgysIFTOEKaXnmqk4xQkGoZjjjkG2KwgCj2CmA6EN6iWMC+//HLODzpJwmTNtFIsVeYymQxNmjRRDuFnzpyZm1SLUjAoaN1lq2wqvOJOrHTt2tXXek4jzaQxG6iglJWVKc2Bzz77bFSRIhNEQfh17jCfBes5mwrCdCP3ashfeeWVvGVu85UPP/ywr/36wZS1WExM1wEXCCHSyUhVRJx++umJptR1Ys6cOYHW//rrrwMfo379+oG38aK8vNxx3mHXXXfN2Wb91hhW4Vapzo04XqKRI0dG2t7v3IK1EYky2e2GVdlVq1Ytd0+sqUGcTEynn366r2PE0RgG2YdfBWE+C1blYza6p512mq/jquau/Di0+Jnb8Iq5MDtgSZqYgngx9QfWAm8KIdYBvwD2t63SeTEVkpkzZwZaf9WqVSxcuJBtttmGunXr+tpm6dKlsZubysvLXfdpepQMGTJEmQ7DD2Fz0MShIFTpNvr06cOee+6Z8wJyY9asWYko5jBYo+etIwjrdXK6ZvZSpU507dpVmQspCEHum1+Hg7feeguoWIzJVBBmhzCMgvAjq9d+O3bs6LmeqSCKZQRRE5gFfAxMAeZSBbyYCknQycGysjKaNGniq3ZxUnMQK1as4PXXX3f1XIqjpxPWpdjtmvqVS3XtHnvsMQYOHBh6exXWwEVVEakpU6bE6rIbVEH4PY843HWDNH5+3xvT+8za4zdLjt52222A2lxntSa0bt0673c/z1Ec74Apa1GMICqjh1IhOOywwxx/a9WqVQW/6k2bNrFy5Up22GEH3nzzTY4//njXfYd9QO6++26uuOKKUNva6dKlC6tWrXKt7eslp5/edVjPmpUrV1JWVqaUYdy4cRVMV1OmTFHuI6py9bv9uHHjcp/t5/vhhx/mJinjahjmzp2r7I3269dPWUsizTTgSXhSmdfNeh72jocqqLJu3bq5bVV11f0oszjOx1QMxTJJDYAQYgshxGFCiH8KIRoIIWoJIfzZMzR89tlnjr/tu+++Fb4vX76cr776CoD77rvPc9+qh9UJq+tpnDVtTbdg1cNqRvF6NWhJVMayYs0ga8Ueh9K2bVvleoWopW3HVA5RsMcL/PLLL8reqFkB0MrXX3+d6nVwim2Io8yvlaAjU1UchLXx//TTT5XbmZHhdv7880/22WefQHIVi4kJIcQ/MMxIEzEyve4DtAcWCCHi6YIWATfffHNBjmvGDZi0a9cu537px8vmv//9r+9jNWzYMPc5ymSxHfOFtXswffXVVznlUcgGtkaNGo6TmK+99pqvfaQ1gohrZOBUQ+GEE06o8L1JkyY5BWH1ElPJ0aZNm1TvYxIJEU1T0aJFm8vY2BviMI2udZt27dopOzxOQZSzZs1Szj2q5r3sxyuoghBCdAJGYcxDDMTI7grGXMS3wL+FEOp0lyWGPWIxLezVyv73v//lJtL84Ob/7kaQkYcXKtdBMKKqTaUUZ8Myf/581q5d63v94447zrGx8evOmVbDeNxxx/le102ZnHLKKcrl9uSIe++9d+6+Wec8iiFLsNM9C5qvTIX1vbEXaApjCrJfL9NMFWW+yM0l3OxEPvPMM6H370SQEcQNGJPTRwM5SaSU3wGHA59geDppQhL2RTR9tsM2XKo6y2ExzyFqXiq/6RKaNm0ayNzyzTffOPbe/PbAoioIv8cJkhvL7dnp1q2b7/0ESYuSRAqVa6+9Vrnc6Xmw58EKg/V+mp000/3XOtL2i9P9Pe+88zy3dbqPbu+TWVTJXlwpDoLc4f2BUVLKvLOXUm7EMDmFz51QSbAWXA9K2CGiaQ6ImvsnDsxzSEtBAEyePNn3uvPmzYsUoAdGg/Ldd9/RsmVL5f5N7CZDE7/3Ocg1iGPSs6yszHdhHYg3PYnJ7bffrlzuNIIw6yuYhAmgVLnrnnrqqQCuzhZOON3fKGWCGzdu7Pib+a4Veg7iL8Dtau0AJFc5pUQ48sgjQ2/rJ/J2+vTpjnVzP/jgA9dtVZONcWPK4taj9dMDT7uKWVD22msvBg0alLfcWivk5JNPrpDh08RvY25fz6kgVXl5eWwKQpVgzuleWlOCJ41TDRZ7fEWYuYrWrVvn7d98RsM06k7Xy9qAO5n9nHCrHmcqxSTqUwdREB8B5wkh8mZLhBANgf8D4kkHWgVQFQzxMxGtyq3kd6jvx9QQtRdivhxujYfbhJtJKdSmUOWSsl6/F154gZNPPjlvHb+NuT2exSmyfsKECbG5gdo96YoFv2a9MCPXatWqMWXKlJwThfV4Ya6r/R0y5yCs+3r11VeV26qUy8aNG/PuvRnpDZsVxCeffBJYVi+CvIXXAA2Bb4DrMWpUnyyEuA+YAdQBCjO7W2RYH4S1a9cydOjQvBvvFg8RFFUW2U2bNuUN1x955BHPfakynQbBzzyKfTJeRbHklXLDy/zw7rvvKidRw56b04Tspk2bYrteKsVsZictJG4N9fLly3PR0G6KxOnZXLduHdtvv32FYmBeCqJp06a5rMR27PfCTOjnls7fNA/bR0B//PEHJ5xwQt5ya+fBVOpuxaDCEqSi3HcYLq2LgH4YXkwXY0xM/wh0kFIGTwZUCbHezMsuu4wLL7ywQtATJN9Dfumllxwn/FR+1iZBPIJURLXvm0StjZEGfnqrqvvsN+7E3ji5Kd8k5yD69I9ex90AACAASURBVOkTy77t70AQrI2uvWGuX79+zkbvdo2sKTWsqDpOXrEF8+fPd0yu6XTfpk6dqlx/woQJ1KlTh7fffjuvJsfKlSs9S/H26NEDSCbLdKBWSkr5rZTyKGBH4BCgHUZ9iIOllM4RYFUMa8/STPlrb/DCesKMHz8+b5mqwVE11HvttReQfhF6O2YAmjUy3N47LrSMfvCj5FVzKU55sq699toKMTj2PFVO9vVMJhNoBPH3v//d8TfVOUXtNJiEzbsFVCjypIpX+f333/nrr79cFYRTnNCBBx6Yt8x8P8MotSCKHTY7thx//PFMmjTJdV92+WBzev04rRImobqxUsqVUsrJUsrPpJTR/cwqGdaXzCn/e5CgNiuqB9Ye8KQ6HpBLJNeli3NC3jR83rfaaiu23nrrnMKC4rN9+3EZ9KMgVIpOVUMADA8eawyOvTPgNGG6bt26QAqiQ4cOjr8lObJ9++23lcs3btzoeb2tDaKTjDVq1HCt+e30bP/jH/9wPF6Ymh72Rt1txA5www03OP7mdM+t16NDhw4sXrxYOd8VFd+uItl0GvcAnTDmIlRd4IyUMnb3EyHE+cCVQBPga+AyKaU6hr0IWL16dW4C00lBxNUrq1atmvLBVxW7b9SoEQsXLqRBgwaO+5s8eTK9e/eORTY37PmQ7I3m2rVrKySrSxtV2Us7bqU9TcxUKVZOOOEE2rVrl9db9MKpN3niiScGSk+iGr3279+fHXbYoSCmvQEDBvDwww/7TqnfpEkTx9/cEiYGSTpoLluxYoUvmazY75OXtcBNud90002+junneQ1DkO7Cf4BewDyMmIdnFH8jYpYPIURPYCjwLHAa8BvwrhDCuTRWBOLoQVu1vnnzo+7XKS+Q04Ph5BLZqFEjV/NNXHMIXjglzDO57rrrPPeRhKxB7pOfyXYnwnicuLlwBpFbFTF+//33O7pPJ41pMnKaI7Djdt3daqh8+umnlJWV5QXXuSmIMHVggpqY3HCKjk4rmj+IgugIPCylbC+l7Cml7KX6i1M4IUQZcDPwmJTyJinlaOBEYAUwIM5jmQSxk/bs2VNpGrE2wObDHzXi01oW0crTTz+tDPRRHc/PQ3XwwQcHFy4EXgrCqSe7++675z77zfsfBFX2zmLBzSd/2LBhvvfjVljK+owECUD0y8qVK2nfvj3z52+uDBBnGU43TNdSe3S0m4IIEwdh39+oUaMqnK/JNddcUxRpTNwIGigXrSZmcJoDTYE3zAVSyg3A24DzTFsEggz7t9xyS8eG287ixYvDigQ4T9rutNNOnH322XnLFy5cGOtx4sZLQTz//PPK5dbheBKybrvttrk6AYXGWssb3EcQqmy//fv3D5yfx+pgoQoEjMpzzz3HxIkTueuuu3LLkuoNSymZMME7NKtVq1Z5y0aMCG8MsSuVwYMHV6jQZ3LHHXcoXdT9UIwjiKeAnkKINENczTqLP9qWzwH2EELE3kIEufB169ZVNiZOwS5RcIosdmpkVZ5Nfs6tWBSEE9Ztkup9qYIRd9lll0SO5Ya9WI3bM6Qyhdx///3KzoMb1ojdJBoht9FCkPt5/fXXe66z5557KhtmO+3bt89bFqXIkUqROwWOJlHnIk6CNPbXY/TcfxBCjAGWYQTLWclIKW+JSzjATK1oj9RZg6HctgViTUAUJCho+vTp3HJL/umqHnSnGgR+cWq4mzVrppwI9ROVreL666/3lVQsKk6T616koSBUL3jQhjbqMVWNcxIpr91I4vqa9RGs+w5jYor6PlmJWxEGmV8q9oDQICOIM4EOQDPgQoyo6RsVf3Fi3jn7k2Muj/3qegWl+FnX+qCb3kz2pGJBcRpB2CNA3fDzIkQ1hfmlrKwscq59PzlvwmAPaGrfvr1j0GFc+JlPiToKDYpTbi8zMEvlAnrLLbfw2WefKc02sPmdsd6jMAoizoa1kDE3Yc+jGE1MgzFMPV2BFsBuir/dHbcOh+nWYLeXbIehHAoebluzZs28mAbrg271zFD19P3i9hD7NX8kkQ44LGFNTNYhuVPqibi9m9566y2lI0CcWM/L6bqkrSCcMAMcVaVhr7vuOg455BBHE9CeexpW46gjiDgVhKrz5eYKHifFPoIIYmJqiBF/EL0CuX/M7Fm7U3EeYndASikL6gJgNtq9evWqEDvg9KBHyWkTRwDT5MmT6dmzZ+T9xEEQBTFixAhq167NSSedVOGFGjJkiOO+4ySp4LFNmzblniFr4++UMt6MyvdDkGJDQTHrj1g9yuw4ZR81TUwjR44M5HllJ0rWZD/EUYjID0EKglkpxhHE1xjmpTSZBfwM5EIEhRBbYoxi3k9ZljzMF8WOU68gik03jh7s9ttvH3kfcRFEQfTs2TMXJWq9tnFV0JoxY4br70kF7FnrEltHEE4lYAcOHMh1113nS54kRxumwjRHA2DEVrRr185zWzOep1OnToDhAqpKie6FEMmWnjHNaEnjFkXtRjEqiIHA+UKIC4UQDYUQiedjzo4Q/g30FULcJoToAryOkQvq/qSP74VqiA3OisBvg6gqdB5Hj8YtB08cXHDBBb7XDWtisioIlW+5fR0/FCq1+AsvvJD77Neb5bbbbvMV7Rwl75EXP/5oDOatjVSPHj0quIg7XVOzpKepwM4888zcb5lMJhFlHKakbtOmTUMdyymgtVQJ8mYMw5gsfhhYAGwQQmyy/cXebZFSPgpcAZwNvARsDxwnpXQOmUwJpx6TU8Pnt2CP6uF88cUXgYopfYOG1yfd63j88cd9rxuHgnBCtd/p06cHPlbSWLOBxt3jd8v/YyZ3C4sfDyInG/6UKVMAePPNN/N+W79+vTJxXlTuvffewNuELcNrnl/SpDWCCDIH8SX53kSpIKW8Fwh+lxPG6SY5NXyqAjNB2XffffnsMyNxbtCoZ6csooUgSQVhrnPCCSfkGiK3BtMtPUNapOkPf+utt1YYTQYtqxmll++mCO2eT2PGjKFz586hjmPtSO23336Bty8mc2wh8a0gpJTnJihHpcJs+MI2PG6T2VaPiyDF6MePH6+soRyWTCbDnDlzHEtBehGHF5MT5n6POuoo3nrrLc455xzX9YuhvGmaCsLuERe0dKgfk5xpSrKz1VZb+TKR7bLLLpES0FmPEcaE6FbisyoR+MoJIXYUQpwhhLhKCNFMCLGTECK+lqeE8BpB2BtPs+fvhdsLFHZoecQRR4TazomHH36Y5s2bhx5Sp2FiWrFiBeXl5Tz55JOu64cpTB83aSoIs0SlSVCXTmuD27NnT2U80G67qXNp/vrrrwCuI4OjjjqK+fPnxzY35OYi7uRo4qcsblDcMs2a+E3dkxaB7oAQ4nJgPjASuB3D3bQdMF0I8Ug2uV6VwamxnjZtmnK53VvDKX7BTQlYfytkoi/TFTOMBwr4UxBnnXVWhe+bNm1yzFJrxYwi95tTqRgURJpR0o8++mjuc5gkktZn8Omnn1a61Ho1sE4OHtb9R7GzW2u+Oyma6tWrc8wxx3huHxd+TMxOitWOm4txnPhWEEKIM4G7gVeAf7A5mnkq8CrQF7gkbgGLGacH2Ckgztq47b///nlZJb32CxXtv34VxPnnn+9rvSCYtuSw5hk/CmLkyJEVvvs9lrX4vNfyLl26FIWCiHOS2iwM5cQpp5yS+xwmICyOnr3bM27mlYqiIKwmLqcRhDUOxY59+aJFiyK71vqJY5k5c6bnOu+88w79+vWLJItfgrq5jpNSngV8ZC6UUv4spTwdGA3E3xIVMU4PsNMLZNX6NWrUcPSUUO3X7E1b920t2emGnzQfH330Ue6zU8UzK2Za7LAKYuHChbzxxhue64XxE3dS0PaUJKNGjeK5554rimjWOEcQRx99tOvvXs4NrVu3dv09DgXhZpo0718UBWGtGaGS18xw4KQg7Ns0bNgwdOZVE6fAzqAcd9xxqblmBzlKS4wYBCfeJP5UG0WNqj40GLWF7XZeqDgnUV5ezksvvaTcXjU5Zw5PrQ+GX68kP7lmrKMbPxHf7777LuCcr8cPfqp1qZIheuGUKtxO9+7dqVOnju9CNUkSZ+Su0wSxidfz4BWU6bfhvuKKKxx/8+N2HEVBWAPdVOdrmiGdal6k5UZa7ARREGswYhCcaArEU0ezRHB7yK0+7ibWicjy8nJHO60qsMc0x4R5cK21n5344osv8o7lhzAlGYuNOL27whLnCMKr0p2XgvCaKPXbe7355psj5TSK0khbZVSd73PPPQe418nWBFMQ7wD/J4TIm10SQrQGLgLGxSVYKeBmu1aZLcaN23x5MplMIBOV2WiHCeDxY2O3ZhMN8nJEzYRperWkhaqecaHdXNevX++YXiMMXi6aXplj7QqiSZMmFXKN+VUQNWvWpH///srfdtppJ8+kiuZxokZXq55Rr3dCjyAMgiiIQRgT0zOAxzCC5voIIV4BJgMbgHCJRYoElVnIit3DwC1fi0pBWM0lbr101cNprp9UXiBro6BK42zFWpYzqoII6oMfFdVoqtAKwiunVNC8QF4KYtmyZa6/W4PMwHiWrUohSOPp1ENfvnw5e++9t699ODlz+E3YF+YZtV+DqopvBSGlXAgchFE06BgMZfEPoBPG3MQhxZD+IgpeBcrtLp1urnxeZpprrrnGv2BsHn0k5dpqtTt7mY2sL7bTyxc2VUHSqBq3QtYDAO8cVo0aNQq0Py8Tkxf33HNPhe+LFi3iiSeeyH0PoiAmTpzo+JtfF2mn47nJYX03VR0Ar1FyofJzOXHnnXcW5LiBroKUcnE2oroe0ABoBNSVUnYrdeUAwesluCkIr8Cn0047LdCxZs+eDSSnIKw9Ji+vnnnz5uU+u7kQ+sF+PkOHDvW1XVhUL36hRxBeBK3w5zWC8GrgnXKGmfc6rcbT61l3k8Mam6F6RovBc80N65wghEs4GAdBA+VOEEJ8BjSWUi6XUi4BHhFCfCmEcPetKwGsNtFTTz2V2bNnu2ZA/de//uX427HHHuvrmF9//XXOI8gP33zzje91g2BVEPbGPZPJOE6iOnli+Z1bsDcCSadxDqogunfvHvgYcadpCKrAvI7vpbxVDWq3bt1yLsdpK4gwIwirKVZ1PmlMQnft2jX0tvassIWaEwkSKHcS8BqwA2DtOk8EagBjhRDJVvFIGLstePfdd68wtLNnb42jIWjTpk2g+rpOrrVRsZol7Mrg9ttvZ6uttlK6g65YsYJPP/20wqgiCHYFEbbqnt8RS1ATk998V1YzkVP6hrTwei69KhDar0edOnVo2LAhN9xwA5lMJnXzi1Pj+P77/krCqBTsa6+9FkkmPzjNnfjlkEMOyX0uKytj9OjRqWWLNQlyp68FJgCtpJS56m5SymeA/YDPgOBO60XEyy+/nPv8yiuvABWDhqyeHHFiLbxSDHTs2LHC9+uuuw5wrlfdrl073ykC7NgVRFivJr8KQhXj4dZD92vSs+5DlT101apVqbkEe82pONWLNrErAPskdalRqDmmqOk6rIqxWrVqdO7cOZF06G4EDZR7TkqZ5yMnpdwIPAe0iUuwQuCUJM+cL6hduzaDBw8GjCjcOLE3yiqSzL1kfRijzisEIa5z8mtTXr16dd6yKA2IGR2/dq17CFDt2rVTsyN7jSC8Gnt7j71QCiKuZyOuOSavAEQ7UV2Xrfeh6E1MGIFybt3ERoC7g3WRY21krO555lCxdu3a3HjjjWQymVC26aikpSCc7LNBJvb8vkz29AVhXwS/sqnOze2YXtd88uTJfP3115HTMMSJV6K8oAqx0COIqI2jm+xBvAmDyhGn3KWgIMYAlwgh8gysQoj9MRL15ef9LVG6dOmS+3zXXXfx3//+N1dHNwn8PABppYR2UhBBju/3gbZP0Icti6oy31x77bV5y1TnEOXlq1evHm3atKnQa3dTKmZt7SAE6Rg4OQ1YCdLYz5o1iz/++CPPqyYN0shWvOuuu/pe1+k5sbsFe60f5niloCCuA34FJgkhvhBCjBRCjBBCfIoRKLcWI5iu0rH11lvTq1evRG+S08tgjTlw87ww6x7Y02H7SbNhZ8GCBY7H8Ivfa2VvrA477DDfx7CiKtOqShnhFUVspXPnzr4Vlt/ztZsdorhbqkZpjRs39twuiIIw3UXDOkd4TYi7ESW9jF/eeustz3XMuUmnOCmn7LlR5bbGkBRqBBckUG4x0Bq4B9gGOAk4HSMm4iHggGwwXclx/fXXFzy03qmhmDFjRu6zm4IoKyvj3HPPzWsUVTmhvHjooYeUy5OoWWC97nHb6FUvVZD6FaNHj/YdUWu1c7v1fO2eclGSHU6YMCFvWdzPcdREhldddVXobdNQEE71IKzE7UQSpiZ4odqnQLM3Uspfgauyf5WGW2+9tdAi+PLLTtLEZPW8ccrmmYTvuGka2m677ejZs2es+w5ia3fLjeWHSy65RNnYz58/v0JqErvy8KPAnRROq1atuOyyy7jvvvtyy4KcQ9++fT3XiXrPFy4M32dMQ0H4GWHHffygkfFJyOCXWMYtQohqQog6QohTvNfWqBpCe+/8wQcfzAuW+fzzzxOTqXnz5rm6EU5eMEn4YJsjp6gNtIk1bXYUG7YfrzIrTq6ju+yyS4UeqF0mP4VfnEZuZWVl3HvvvXnL/JDJZPjPf/7juV7UNORRTEymCU3lNqzCqw6GCrfrZc8Vtf32bsms83FKRhimip/K+y4NggTK1RFCDBdC/CaE2CiE2GT+YSTqWwl4z5CVAEEmrsJQr169vGX28pr9+vXLmxgMO4HrF9NMkslkWL58OWeddVYF198g5iq/DZV5zN9//5333nvP1zZudTCs9+6MM87wtT8VQVOAq+5pXDhFl6uucdw9Ta+Mq144pe3wQ7NmzZg4cSKPPfaY57rbbbed7+R/VoLY9oNeW6cOVZhcWVZTc5oEGUHcCfQCZgLvYSTrew54H0NB/AkESzBUpJx99tmJ7l+VWqPQ0bdWysrKuOGGGxg5ciRPP/10bnmQ3uRJJ53ka73XX389Z2bym0bEKZhu48aNFfJpqQov+SVoY+C3oQkzqgliKiv0XJqdoDnH7Pztb3/zrE8BRgCkWwrv3r17M2LEiLzqikEUxMCBAwF8KSwwZFcRRkEU/SQ1cDzwqpSyHWDmHx4ipewEHApsAoK7zBQhYQJcgngLqQJ39t9//8DHTBJTRqvXz5NPPul7e78vEeQnJvOiT58+yuX2if40G0tr49SrVy/H9YJ4LdWtW9czNbedYot4TvMemKNHVcM8fPhwzjrrrApVHSGY8h00aBBTp071XePdKR7FK05FRaGiwYM8TfWBsQBSyl+AhcAh2e9fA8OBM+MWsBCEeajHjx/P/fff72vduJO5JcHDDz8M5OescXKBteOnSJFJENdTcB4ZhM0HpeLBBx8MtP62227LK6+8wuLFiyOnWDD58MMPA5toim0EkSb9+vXjtddeU3p3mYTpRFgny/125NzmDMKMIEpBQay2rT8b2Nfy/VuMsqMlT5iXrH79+r4zuBZremmV+WPOnIpZ3JPwZApq53a6P3ZbfdqN5SmnnOJp1gpiYgozGihGBWFOhictW/Xq1TnppJNcj5NWmu9atWo5yhGmg1gKJqYvgH8KIUxVNgNoL4Qwr4KgxFNtmIS1XXsVHDJJWkEMGDAg0vbWB9uenyquB/Xxxx/PfY7iCllqBGmgwjhLFKOC6Nu3L5lMpmA1DazYFbSf51l1TWvVqhVqOwg2ujaxukqnSZC3/R7gcOAHIURd4CmMOYexQohHgUuBD2OXsAD4iUZV4ZTsz44f//MomHn7g2K+PNZJYPs5xdUAWRtKr0R3dvzKEFTWRx55JND6bjhNzgZREHXq1Al83GJUECZ+R09Jlvs88cQTK3x3y0JsjkhVcUEXXnih57FUpqT69et7jiD23XffvGVBTZ5xESSS+gOgCyCB36SUk4ELgcOAvsBXwGVJCJk2YV8yv6mqw+TjSQNr7QOzZ2X3XIprBGGNGwjq2RMm0MgPdg+XKLz00kvK8/JSEH56pm4kpSD8eBLFRZJFo+zWAbc4jaeffpoxY8Yolchll1Vs6qzR0Z9++ilgKAM7Xbp0oUOHDq4yNmjQwPX3NAlacnSslLKLlDKT/T4Mo4BQPSnlYaWaaiMu/KaicHqJM5lMLAnKwgY3WSO1ncxgQb1qnLA2OPZ5Di969+7NY4895nvEVkx4KYgwZiXrRHZSCiIO273fZzvJCdkgHZxatWo5xh7Z3w9rzIOby3p5eblnJmjTBGX1xrKPfNIicndQSrleSunP+F7JUb0A3bp1ywuCS9oMELaXb7WNWpWM1eQ2ffr08IJZsF6DZ5991nN960tXrVo1zj//fMeUICZpZAMNildD+9RTT1G/fn3GjRvne59nnrnZeTCpycw4Gm0/9+Mf//gHL774YuRjORHXuxe0NoTJjTfeCMCkSZP4/PPPlaMF0wRldYctVJyU49MkhCi3Rkv7/Eu+0GsCDBoUTxJalf/1Cy+8kPdyxfGQuj2gZlTvkCFDAu3TafLMOolcrVq1WIbAbg2ZygYdpoEqRm8xr0ayRo0aLF26NM8jzq1AlTXdRlwNoFkIySQOBXHllVd6rnPnnXcqM/PGRVzXJ8xEM2ye82jXrh0HH3wwn3zyieO+rSVVX3311VDHi4pbd+OZEH8jkhQ2KYqlxq5funTp4lj+E4zzyWQyXHLJJYH226JFi9xnpzrZZWVlsZgb3K65tUfsZ/0whEnLEAdeuYmcPLrcItOtjXdcDaA9oOyuu+6KvM+rr77ac52klXohAwlVSsVUxBdffHFu2Y8//pi33nbbbZecYC443g0p5bkpylFQ4tLOfrOthu19mPTt29fTvBIG67yAkwdN9erVYzHdBH1R42r4vvzySw488EDfBe/jpnnz5nz88ceOv7sl5vNDXNfJvp+wNceDknRAWCEVhDXrrhXzfTKDU810MXXr1s05vhSqbn1gdS2EqAMcCzTDSK/xA/CBlDJaVq8C8v3331f4HrYBdFIQ1v0NHz68YFGRXlgbBSclts8++3iOIPxEmxbKHfOAAw4o6NyE17Gj+s7HdV3tDWla9yvpdyOp/Xfq1ImxY8e6ruP3GprtyNFHH80rr7wCwPHHHx9NwJAEUqdCiKswUmz8D7gbuA94C1gkhDgnfvHSIa7oSj8Konfv3oH3O2zYsArfk/JosDYKbqmTM5mMa0ZVP/W63Xpyqka02HIMhcXLS8lJgVSvXp0DDzzQc/9x1Qyxl3BNy8SRtIKIM83NdtttR+vWrQHnGiphMKtCtmnTJrcsTH64OAiS7rsPcAcwBaOSXGtgf+CfwHfAf4UQhfHFipmwgTpOL2dUDwR7Pp6kUktbG2GnFymTyfDrr7+6Zrz18zAHbfAfeOAB5XI/jWYxce211/qqG63CGmDpdI3jqh1tD16My8ThlfAxjIIIU6EtDtasWcO0adMAfxmg/Y4gzPk/a031QlkdgpiY+gMfAB3NOIgs04QQr2BEUV8PvBGjfKnTqlWr0BNlTnZ7e53ooNhfVtUkbhxYG+3hw4cr1zG9Lpyytd57770cfPDBgY7lByfPKSEEX375ZaB9FZItt9wydArs8847jwYNGrDXXntVcCiw0r59+yji5WjQoEEF02uYDKQqzj33XNdst2HevTSD+Jw49dRTY9tX165dgYpKoRRyMTUFXrEpBwCklJuAF4DCuIbESJQb4eT5EzWVt71GdJhskH6w9nCcql6ZhUucEuz5Hd389NNPjr855dFXkUTywELi9vyVlZVxwgknOCoHiM+EYjchFrOJKWg24DSxjrz8jiCOPPLIvGWFctkO0hrOBNze3H0wMryWNFGHcqpUz+bL5rd0op3JkydX+F7IHpPbaGjo0KG+60q77eecc/Kns5xeLr/R60FxKiGaNFEng5PyYkqLyqYgpJRccMEFgP9rao9BgdIwMV0MvCOEuBf4t5RyOYAQYpvsb2cDneIXMV2iDuVU25s5WeKaaEpqBOEHJ9MTGCaQOIbC9hepdevWjvtNagThNBpMGlXvMQilriDCPD9R62bHiVucS9quynEQREGMAMqBAUB/IcRK4C+gAUb50TJgoi3RVkZKWXzhrC5E1dQtWrTIq0XbpEkTfvjhh9h8yf2W5oyCH7c9O0Fe7iDXeauttnJ084x7BFGzZk3+/PNPZS8uDaJmMtUKonD8+eefSvmDulWXqoKYBBRfcpsY2G233Zg7dy4QfQTRt29fZVoEN7uxF0cccUSF4Ko0KtIlrSCCzMvstttubL/99srf4lYQF198Mffcc0/olO9RCHq9VVRFBVEsJiavkX2lHkFU5sjqnj17ctNNNwHRRxBu8QFx7TON3m2Q2JDu3bvz/PPPJyaLW9nNuCOizcnAuOIJgtCxY8fI+ygFBbF48WJWrFihrHsQ5rjFoiCcCDqCCJsIMAnCRFJvI6Vcl/28A9Ad2Ai8KKVcGbN8qWBtDKOOIJIoaWg3r6RhHw/SQF577bWuyeSikuYEnXmsUvWOKgUFsfPOO+c9Xw899JBrHWcVtWrVYs2aNYHrdqeNtaa1H1R1JApFkEC57YUQ75CtGieEqA1MBYYA/wG+FUIUxnAbEauGX7VqVaR9JeGvbHdxO/fcc2M/hp0gZUDD9La9elUfffRR7nOaPuDnnHMONWvW5Iwzzkj0OH6izcNQCgoCKpoG69Wrx8UXX8w111wTaB9m/MuYMWNilS1uzHONmoOtEAR5824FjgHeyX7vDewCXAkcjTGBfWus0qXEwIEDc5+jTgCnMYJIo8E0E4f5Icw5q7axupZaK39Zs4DaU4DEbW5r0aIFf/zxB82bN491v3aSmkeKq2FPegRlfYYfffTRUPto0aIFmUwmVJGlNKkqCuJE4CEp5eDs91OAZVLKe6WU44FHMJL4lRxhav86YW34VDbWMBRjXQMrc/Z+bAAAEnRJREFUcTUm1hTc1gk/M4p6wYIFvP322xW2CVqNrlhIqrGIq7i9U2qTuLA26nFFaQfFLN6TNOb7UdkVRH1gOuQyuh4GWN0uVgDJVRsvEaw9z7hqJ1sVRNSo7CQI40mkMjFZRwfWczZ7xY0bNy6KtApx4CcdSRjicvls0qRJLPspZgYPHuy9UgyY70eYjt68efNYtGhR3CL5JoiCWAiY4/mTgeoYmVxN2gHzY5KrZLEWnb/tttti2ac111AxpgoPM4JQZb+0BhLuuOOOUUQKnRAvLc4//3zP4kFhiDqHVgiKya0zCQ477DAgvwiTH5o2bUrDhg3jFsk3QRTEmxgBckMwUn2vBN4UQjTKLusJJObrKISoJYT4SQhxelLHiJu4evtWN75iTHsdptdar169vPkea0MR1exQyJfKD2VlZUybNi01M4emcAwcOJAffvihQvruUiFIa3MlhgI4D/gV+KeU8g+gCXARMBL4d+wSYigH4HWguGejbMTVM7IOTYN6eqRB2GA1+xxNGOVnLdVoJa0KaFGoW7du7GaOQqUIiUJlH0FUq1YtUqBsIQkSKPcXcH72z8rXQGMppTr9Z0SEEEcCQzFSepQUcT347du3z+Wdd6tNXCjiimYOoyD22WefvGXr169PJdq8GClGE6QXlV1B+OWdd95JpJRwFCLbK6SUfyWlHLK8BnwL/D3BYxQ1/fv3L7QIrsTlxRSmoejTp08FN+UDDjig5JTDe++9l5fSPSzFaILU+OO4446LrZ5HXDiOIIQQc4D+Uso3LN+9yEgpg8/EuNNeSjldCNEs5v2WDIXM3urEb7/9lsuPFFcPMEw8RVlZWYXJ3ksvvTQWWdKkQ4cOdOjQIZZ96RGEJk7cTEw/Ab9bvs8nxmR9QogtATdlslRK+auUcnpcxyxVivEFql69OqNHj6ZLly65ClhR+emnn0IFqK1Zsyb3uVRTZMRFKSqIxYsXp3q8yZMnF316jmLBUUFIKY+2fT8q5mM3xqhl7cQAINlonRKhWBVE586dAycis/Pwww/nJprDurZaA5CKzYabNnEqiIYNG6bSeKcd7FioGtalSMFCdKWU8zBqSGg8KEYFEZfZq2vXrjkFETZlwu+/bx7oWlN0VEXinINo2bJlKgpi1qxZiR9DEw63OYgPwuxQSnlMeHE0KtJUEHvttVeFYvVOxNUQWd0yw45GrPmYStHEEgdNmjRhwYIFsT4ro0aNomXLltxxxx2x7VNFsaeSqcq43ZndyZ9zaADUxIiDmIXhBdUM2BH4BXeTkSYkaSqIjh07+lIQSRBWQVhHDXHm1SolJk2axBdffBHrPuvXr88vv/wS6z5VFCoXk8YbtzmIZtbvQoijMFJrnAs8K6Ust/x2BvAERsI+TcykqSCslevSol69eqxcuTK0grBGqO61115xiVVS7LrrrkWf1dSJuKsCauIjiJ1gCPCElPIZq3IAkFKOwlAON8cpnCZ9zIC8NHFSgAcccICv7aPWcdYUlmKpKa3JJ4jxrzkwzOX3BRieSYmQ9KT25MmTadu2LX//e/HF4xXjJHUS2EcQkyZNYt26dZ7bqTK/akoHPYIoXoKMICTQXQiRNwsohKiJUUAoWrWdAnLQQQexceNGRo8eHXlfcRebqeyNntP51axZk3r16nluX1UnpisLVT12pZgJMoL4NzAKmCiEeBKYgzFhvSdwIdAUiCdiqkDE1dB88cUXLFkSX/aRyq4gTMLOQVjjIKrKtapMhClZq0mHIMn6XhBCbI2hKIay2cOpDJgHnCKlHBe7hCVI3bp1qVu3bmz7q+yNXrNmzVixYkVoBV3Zr09lJ+6ysZr4COTMLqV8CmgEHAz8E+gGHCil3F1K+bbbtprwVPYG8K233uL5559nhx12iLyvyn6tKhOnnXYaAEceeWSBJdE4EThCJevBNCX7p0mByp6hs0GDBvzzn/8stBialDGdC/QcRPGiQxhLgGLrFccdkBUnxXatNM4MHjyYGTNmFKXnoMZAK4gSoNgavWKOfI2aPFCTHi1btuTbb78ttBgaFyq37aKSkKaJqWPHjp7rFLPJS7u8ajTxUbxvuiZHmiOIV155xXOdYlYQGo0mPvSbXgKk2SBvt912nnUZ1q9fn5I0Go2mkGgFUQKkPQfhVfpz0aJFKUmi0WgKiVYQJUDaJh2vyFY9EazRVA20gigB0h5B6NQHGo0GtIIoCdJWEHvuuafr78U4gihm11uNplTRCqIESNvEdMQRR+QtW7BgQe6ztUBPsTB//nxd21ijiRmtIEqAtEcQvXr1ylvWuPHmUh+1atVKUxxf7LTTTrGnWddoqjpaQWjyaN26de5zjx498n4vRhOTRqOJH60gSoTbb7+dr776KtVj7rrrrjz11FOsXr26wnKtIDSaqoHOxVQiDBo0KNXjTZ06lV122YUtttiiKE1KGo0mebSC0CjZf//9HX/TIwiNpmqgTUwajUajUaIVhCYwOuZAo6kaaAWhCcy2225baBE0Gk0KaAWh0Wg0GiV6klrjmz///LPQImg0mhTRCkLjmxo1ahRaBI1GkyLaxKTRaDQaJVpBaDQajUaJVhAajUajUaIVhEaj0WiUaAWh0Wg0GiVaQWg0Go1GSWVyc60OsGTJkkLLodFoNCWDpc2sbv+tMimIhqAucKPRaDQaTxoCs60LKpOCmAy0BxYDmwosi0aj0ZQK1TGUw2T7D2U6t79Go9FoVOhJao1Go9EoqUwmJleEEOcDVwJNgK+By6SUn7qs3wp4EDgEWAk8AtwlpSyZIVeIc24H3AbsD6wD3gOukFIuTUHcyAQ9X9u2NwKDpZRlyUkYPyHu8U7AvcDxGB3Ej4H+Uso5KYgbCyGf67uANsBy4GngdinlhhTEjRUhxInASCmlax3guNqvKjGCEEL0BIYCzwKnAb8B7wohdnNYvz5G45gBugGPYTScl6cicAyEOOeWwPvAGuAMYCDwt+w2W6YidASCnq9t21ZAukW/YyDEPd4SGAccDJwPnAvsAYwRQmyVhsxRCXHOewBjgbXZ9e8HrgLuSEXgGMkqumcB105MnO1XpVcQQogy4GbgMSnlTVLK0cCJwApggMNmF2GMrk6UUo6WUt6K8UANKpHGMsw5X4wxwX+alHKMlPI5oDtGr6tjCmKHJuT5mttWB4Zj9CxLhpDn3BPYE+gkpXxZSvka0AOoBeybgtiRCHnOp2NMwp4mpRwrpXwIeAC4ILu/okcIUUMIcSXwIbDRxyaxtV+VXkEAzYGmwBvmguzQ8m3g7w7bHAu8L6VcZ1n2GlAPaJuQnHES5pxnAPfaht0y+9+zF15gwpyvyQCgNvBQYtIlQ5hzPgV4R0o537LN11LKRlLKL5MUNibCnHMNYAPwh2XZL8B22d9Kgc4YI9wr8PecxtZ+VQUFsWf2/4+25XOAPbI9SNU2qvWt+ytmAp+zlPJRKeUjtsUnZP9/H7N8cRPmHiOEaA7ciGFuWZ+YdMkQ5pxbA98LIQYLIZYIIdYLId4WQuyaqKTxEeacR2K4vd8hhKgnhGgL9AdelVKWSgWsycBuUsohGGYjL2Jrv6qCgqid/b/GtnwNxvmrCizXdljfur9iJsw5V0AIsQtwDzAF+CBW6eIn8PlmzQtPACOklBOTFS8RwtzjnYBeGL3t84Czgb2Bt4UQpeCwEvicpZSzMebTBmKMHL4AlmFch5JASrlQSvlbgE1ia7+qgoIw7Yx2zWsuL3fYxklTq9YvNsKcc46scngf4/noXgKeW2HOtw+GyeKqpIRKmDDnvCWwFdBZSvm2lPJ/wD+AVsCpiUgZL4HPWQjxL4yOwBNABwylWBdDKZaKiSkosbVfVUFBrMr+t7uFbYdxsX532Ma+fi3Lb8VOmHMGch49n2D0NDpme2DFTqDzzSrAu4BLgXXZ3nO17G9bCCFK4b0Ic4/XAp9be6NSyikYnkBFP0lNuHO+GhgtpewjpfxASvks0AU4HGOCvjISW/tVCi9CVGZl/+9uW747IB16x7Mc1ofNE7fFTJhzRghxCIZf/CagvZTym+REjJWg59sB44V5CWMCcwNGbADZzzckJGechLnHP2KMIOxsgT/bdqEJc867AJ9ZF0gpv8cwN+0du4TFQWztV1VRED8DJ5sLsq5eXTHMKCreB44VQlhtmidjPFRfJyRnnAQ+ZyFEM2AMsBRoJ6WcpVqvSAl6vm9ieHNY/+7L/tYWw2+82AnzXI8F/iaEaGTZ5kiMHvgnyYkaG2HO+QeMeJ4cWeeEHYC5yYhZcGJrv0phYioSUsqMEOLfwMNCiF+BSRg+/ztiBM2YwTQ7SSnNnsajwCXAaCHE3RixAIOAq6WUf6V9DkEJec4PYpiVLgJ2tXm2/CSlXJzaCQQk6PlKKX/BeFlyCCEOz+5rSqrChyTkPb4f6I0RGDcY2Aa4G0M5jE35FAIT8pxvBv4nhHgCGAXsjOG5Ng94JtUTSIgk26+qMIJASvkohg/x2Rhmhe2B4yzpBa4HPrWsvxjDl3iL7PoXANdKKe9JU+4oBDnnbC+sC0ZA0XPZ5da/orfVBr3HlYEQz/VyjN70XGAE8DBGZHVXKWUpOF+EOecXMSKoDwBGYwSMfQwcIqW0e/qUKom1Xzqbq0aj0WiUVIkRhEaj0WiCoxWERqPRaJRoBaHRaDQaJVpBaDQajUaJVhAajUajUaIVhEaj0WiUaAWhKXmEEPWtUaNCiKeEEAX33y4WOTSasGgFoSlphBCdMfLL7GRZPAwjkKrQFIscGk0oKn2qDU2l5xCMaNoc2QL2BY+aLhY5NJqw6BGERqPRaJToVBuakkUI8RRwjmXReCnlUeZyKWWZZb2DMIoE3QPsBywGbsLIPXUzRoWxrTByE/1fNqGfeZy9gduAo7PrfAXcLKV81498NjkOxTA73YOROXYN8AJwlZTyD/WechXwrsfIi9UUI6//WOAaKeXPlvXqZs/nVIwkdnOAocAQazpsIUTt7Pmfll1vNvCglPIJt3PSVC30CEJTygwDXs1+HoDRiDvREHgLmABcDmwE/otR8P4Y4BYMZdENo/EGQAixL4aZaG/gduBajMpso4UQ/wwhc32Mhv17jIJFkzAyb97ksd01wGDgHYyMu49jpHAea9Zizk7Uf4yhgJ7GqL08HXgAIzGfeU5bZde7BOOaDMBQJI8LIfqFOCdNJUXPQWhKFinlp0KIb4BTgNeklPNcVq8HXCKlfBhACDEPQznsCQgp5frs8v2ATpbtHgKWAwdIKX/PrvMQRp3uB4UQrwZMoVwX6CelfCj7/XEhxEyMkcGVLtv1AMZIKS81FwghfgYuBJphjACuyJ7PQVLKb7Or/UcIcTswSAjxmJRyGkY96jZADynlc9l9PQaMz673iJRyU4Bz0lRStILQVCVetXz+Ift/jKkcsswF2gEIIXYAjsRQElsLIba27es+DDPRpIBy/M/2fRpGbWg3FgBHCyEuBZ6XUi6VUg7DGEWZnIYxYlgshNjRsvw1jHoAx2ePdTyG0htlrpCttXA2xuioJFJ/a5JHKwhNVWKp5fPG7P9ltnU2YRR9B9gj+/+S7J+KXQmuIJbbvq/HqMXhxkCMSngPAPcLIb4E3gAel1Iusci7tWL/VlkhO+Kwl+iUUv7kS3pNlUErCE2VQUq5UbHYzUvDbLQfweiFq5gRQo7APXQp5TdCiBbA34ETsv9vBi4TQhyWrbNcHZiI83zGouz/6uhRgsYHWkFoNM7My/7fKKV8z/pD1rNpN2Bd0kJkJ6HbAKullG9gjBwQQnTD8IA6H2PifR5QSyFrXaADRk1ngPlAa8VxOgPdgSullEvtv2uqHtqLSVPqmJOpsT/L2dKNU4BzhRCNzOXZEq3/xSjnmEYnqzrwIYZ5ycrn2f/mNXgDaCOE6Gpb7zrgRaBV9vtooIEQ4hTbegOArsCKOITWlD56BKEpdUx7+xVCiDHZHnac9MPwWPpSCPEo8AtwBkYE9yBrvERSSCn/EkIMAa4TQryK4eq6DUat4XUYygqMesunAa8IIYZimL8Ox3B7HZP9A2NiuzfwvBDiEYxUJV2BjkBv7cGkMdEjCE2p8zzwHkag251x7zybLuNvGCOJy4G7gW2Bc6WU/477eC4MBi4DmgP3Zr/PAY7Izj8gpVwJHAY8heEVNQQjMO8W4HRz7iMbkHcUMBxD2d0PNAa6SSmfTO2MNEWPjqTWaDQajRI9gtBoNBqNEq0gNBqNRqNEKwiNRqPRKNEKQqPRaDRKtILQaDQajRKtIDQajUajRCsIjUaj0SjRCkKj0Wg0SrSC0Gg0Go0SrSA0Go1Go+T/AZnAvnDqNC5+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(ts_all.size)*delta_t,ts_all,\"k\")\n",
    "plt.xlabel(\"time in sec\")\n",
    "plt.ylabel(\"displacement in µm\")\n",
    "plt.savefig(\"membrane_ts300.png\",format=\"png\",dpi=300,bbox_inches='tight',facecolor=\"white\",backgroundcolor=\"white\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc7305f13c8>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAECCAYAAADq7fyyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXRklEQVR4nO3dfYxU13nH8d96wTY1EIr8DrV5sx9XMQlITYyD3LolFgQaaxNHcRpsK1VFkyhuiklpQThysaGgOo5p6rqJo7ZSBIloaYrdQkzJ5rUUEkcyhaThKa/B4LemeG2C17As0z9mhgyzM3PPzM7MfZnvR0KDz9y7+3C0/t2z5557piuXywkAkB0XxV0AAKC5CHYAyBiCHQAyhmAHgIwZEXcBZnaJpHdJeknSYMzlAEBadEu6RtJz7n669I3Yg135UP9+3EUAQErdJuk/ShuSEOwvSdKGDRt09dVXx10LAKTCyy+/rIULF0qFDC2VhGAflKSrr75aEydOjLsWAEibIVPY3DwFgIwh2AEgYwh2AMgYgh0AMiYJN08bsvn543p0m+vFvn5dO26Uls419cycEHdZABC7VAb75uePa/nX96p/IH8z+Hhfv5Z/fa8kDSvcG7lYhJxT65hK70kKauNCBqCSrri37TWzSZIO9/b2Bi93nL32Wzre1z+kfcK4Udqx7HcaqqP8YiFJo0Z2a80Hp1cN0JBzah0jach7Iy/qkrqkgcFcZNvoS0eo780BXTtulH77piv07X3/qxf7+vW2USPV1SX1vTlQ9e9cHIB0O3bsmObMmSNJk939SOl7qRyxv1gh1Gu1h3h0m18QsJLUPzCoR7d51fALOafWMcW/lxo4N/RCW63ttTcHJOV/Y1m/6+j59/r6ByL/fryvX4s37taSf9ytczmpu6tLg7lc1dcJXAiA1EhlsF87blTFEfu140Y1/DUbuViEnNOKi1AzFa8Zg4Xf3Kq9Fi8EizfuvuD8i7ryX4PgB5IjlcG+dK5VnN4ozkU3opGLRcg5UcdUei9NiheGasEvSb/6KyP10PvfTugDbZLK5Y49MydozQena8K4UepSfrRYay48xNK5plEjuy9oi7pYhJxT65hK7428qEsju7si29LktTcHtHjjbk1atuX8n5kP/7s2P3887tKATErliF3Kh3szR4Clc+KhK09Czgk5pt5VMW8bNVKnzpy94GZq2hTDvjjCZ1QPNE8qV8Vg6DLJ0FUxZ84O6s2Bc3GXH+myi7u1+gPD+y0MyLLMrYrB8H5jKV4Ujvf1R66G6ZIUx6X/1JnB8yP6i7qkj95ynVb1TI+hEiB9CPYOVO9FYfPzx/Xnz/zkguWS7XQuJ63fdfT8kk5G80BtBDsi1boQxBH6paN55uaBoQh2DEvcoV96E/aeWUzXABLBjhYqD/1WB31xuqZL0kJCHh2MYEfbVBrdP7h5rzbsOtrUG7Q5/TLkGcWjExHsiNWqnukXBG+zg74Y8LOnjteGRbc26asCyZbKJ0+RXat6puvw2gU6snaB1t09Q6NGNudHdMfBE5q0bIsWfnlnU74ekGSM2JFYpVM3+e2P96h/mA9XFQOeETyyjBE7UqFn5gT99JH36cjaBbpn1nXD/nrFgH9w894mVAckC8GO1FnVM11HCtM1ww359buOavKyLWxIhkwh2JFqxZBfd/cMNTodn5O0eONu3bJ6e1NrA+JCsCMTemZO0P6/GN4N11dOnuEGKzKBYEemlM7Fz546vqGvUZx/Z3oGaUWwI7M2LLp1WAHP9AzSimBH5hUD/oYrL6v73OL0DKN3pAnBjo6xfcntWnf3jIbOXbxxN3PvSA2CHR2lZ+aEhpdJ7jh4Qjet2NqCqoDmItjRkYrLJOudnnlrMMeDTUg8gh0drdHpmfW7jnJjFYlFsKPjFadn6h29v3LyDFMzSCSCHShoZPRenJoBkoRgB0oUR+9jL+mu6zyWRCJJCHaggj0r59W9coYlkUgKgh2oorhyZkRX+Dk7Dp7QHZ//TstqAkIQ7ECEA2vqm5rZ/+opVswgVgQ7EGDPynl1rZphxQziRLADgepdNfPWYE7TlrNiBu1HsAN1KK6aCXU2J8IdbUewAw04snaBQu+pEu5oN4IdaNBhwh0JFRzsZrbIzPabWb+Z7TSzWyOOf4+ZfdvM+szsRTP7ipldNfySgeQ4XMdySMId7RIU7GZ2n6QvSlov6S5JfZK2mdnkKsf/uqReSScl/Z6kP5E0u3DOyCbUDSRGPcshz+bEahm0XGSwm1mXpIclPeXuK919q6Q7Jf1c0gNVTrtf0kuS7nL3b7j7VyV9RNI7Jd3RlMqBBNmzcp4u7Q4bur81mGOdO1oqZMQ+TdL1kp4pNrj7gKQtkuZVOecnkh4rHHf+tMJrxVE+kHb7Vs8PDvdXTp5hT3e0TEiw31h4PVDWfkjSVDMb8juouz/p7n9T1vz+wuu++koE0qOecF+/62iLq0GnCgn2sYXXk2XtJwvnRz6OZ2a/Julzkn4k6Vv1FAikzb7V84NvqLLlL1ohJNiLP6K5Ku3nap1cCPXewvf6iLuXfx0gcw6sCV8tQ7ij2UKC/fXC65iy9tHKh/qpaiea2c2S/lP5Uf8d7n6wkSKBNDqwJnyd+zseeraltaCzhAT7/sLrlLL2KZK82gjczG6R9D1Jg5Juc/c9DVcJpNThwO0H3jg9yF7uaJrQYH9BUk+xobAWfYHyUyxDmNkkSd+Q9Iqk97j7/krHAZ0gdG+ZHQdP8ClMaIoRUQe4e87M1kp6wsxek7RD+XXql0t6XJLMbKqkK9x9V+G0v1J++uVTkq4zs9KPovmZu7/UxH8DkHjr7p6hxRt3Rx63eONu9cyc0IaKkGVBT566+5OSlkq6V9ImSeMkzXX3Q4VDPitpp3R+ND9fUrekrxbaS/8sbGL9QCr0zJwQvJ87N1MxXF25XLyLVArTNod7e3s1ceLEWGsBWm3ysi1DlpdVU8/2wOg8x44d05w5cyRpsrsfKX2P3R2BNgq9mSqxUgaNI9iBNgsdib9xepBtB9AQgh2IQehH7LHtABpBsAMx4GYqWolgB2KyfcntwdsO3PH577S0FmQLwQ7E6MCasPn2/a9W3bkDGIJgB2IWejOVKRmEItiBBAi9mcqUDEIQ7EAC9MycEPS5qUzJIATBDiTEnpXVPmnyQkzJIArBDiRI6Hw7Dy6hFoIdSJiQ9e08uIRaCHYgYbYvuT3oOPaSQTUEO5BAIVMyb5webEMlSCOCHUiokFUy3EhFJQQ7kFChq2S4kYpyBDuQYLOnjo88hhupKEewAwm2YdGtQcct/PLOFleCNCHYgYQLuZG64+CJNlSCtCDYgRQIuZHKPjIoItiBFAi5kco+Migi2IGUuGrMxZHH3LJ6exsqQdIR7EBK/GDFHZHHvHLyTBsqQdIR7ECKhOwjw6gdBDuQIiH7yDBqB8EOpEzIqJ0VMp2NYAdSJmTUzgqZzkawAykUMmrnadTORbADKRQyaudp1M5FsAMpFbKunZ0fOxPBDqRUyLp2dn7sTAQ7kGIho3Z0HoIdSLGQUftkPmWp4xDsQMpd2t1V8/1cm+pAchDsQMrtWz0/8pibVmxtQyVICoIdyIDaY3bprUHG7Z2EYAcy4HDApyyxOVjnINiBDsHmYJ2DYAcy4p5Z10UewzYDnYFgBzJiVc/0yGPYZqAzEOxAhrDNACSCHcgUthmARLADmRP1wBKyj2AHMibkgaVJbDOQaQQ7kEGM2TsbwQ5kUMgDS4zas4tgB4CMIdiBjDoSMGpnS99sItiBDsbWYNlEsAMZtu7uGZHHTFvOqD1rCHYgw3pmTog85izD9swh2IGMmz11fOQxbOmbLcHBbmaLzGy/mfWb2U4zuzXwvDFm9jMz+1DjZQJo1IZF0f+rsqVvtgQFu5ndJ+mLktZLuktSn6RtZjY54rwxkp6WFL2fKICWueHKyyKP2fz88TZUgnaIDHYz65L0sKSn3H2lu2+VdKekn0t6oMZ5vyXph5Ki794AaKntS26PPGbxxt2tLwRtETJinybpeknPFBvcfUDSFknzapy3WdLeiGMAtAmbg3WOkGC/sfB6oKz9kKSpZtZd5bzb3P3Dkl5ttDgAzROyORg3UbMhJNjHFl5PlrWfLJxfcfLO3X88jLoAxICbqNkQEuzF39/KV7sW2881rxwArRSyzQDSLyTYXy+8jilrH618qJ9qakUAYsWTqOkXEuz7C69TytqnSHJ357k1IEWiPheVJ1HTLzTYX5DUU2wws5GSFkjqbVFdAFok5HNRkW6RwV4Yka+V9AkzW21m85V/6OhySY9LkplNNbNZLa0UQNvwIRzpFvTkqbs/KWmppHslbZI0TtJcdz9UOOSzkna2pEIATReyfwzSqyuXi3dCzcwmSTrc29uriRMnxloL0EmiRuWzp44P2mcG8Th27JjmzJkjSZPd/Ujpe+zuCKCiHQdPxF0CGkSwAx2KNe3ZRbADqOodDz0bdwloAMEOdLAREfuCvXF6sD2FoKkIdqCDHVjDdEwWEewAamLHx/Qh2IEON/aSajtv57HjY/oQ7ECH27OSz8LJGoIdQKQHN++NuwTUgWAHELnj4/pdR9tUCZqBYAfAjo8ZQ7ADQMYQ7AAkRT+sxFOo6UGwA5AU/bAST6GmB8EOABlDsAMItvn543GXgAAEO4Dzoj5ZafHG3W2qBMNBsAM4j09MygaCHQAyhmAHUJfJEZ+VivgR7AAuEPWRebk21YHGEewAkDEEO4C6TVvOdEySEewAhoha9niW+ZhEI9gBDMGyx3Qj2AE05KYVW+MuAVWMiLsAAMl0w5WXaf+rp6q+/9Yg8zG1TKqwLDRqxVGzMGIHUNH2JbfHXUJqVQr1Wu3NRrADaFi7gipLJi3bcv5PqxDsAKq6tDvi0zcwLK0Kd4IdQFX7Vs+PPIabqMlDsAMYFm6iJg/BDqCme2ZdF3cJqBPBDqCmVT3TI4/hJmqyEOwAkDEEO4BIIQ/WMGpPDoIdAJooCRc4gh1AkJCbqEkItbRpRZ8R7ACChNxERTIQ7ACC3XDlZZHHMGqPH8EOIBgbg6UD2/YCqMuRtQsiR+XveOhZ7Vk5r00VNV+cW+42AyN2AE33xunBuEtoWK0td9MyzUSwA6hb1GeiZlkawp1gB1C3kM9ETUMAlguteXLMH6QRhWAHgDolfT9Lgh1AQ9hm4EJJ+rcS7ABaauGXd8ZdQkuUTsckKdQlgh3AMKy7e0bkMTsOnmhDJcNXbzjnlA/3pIW6RLADGIaemROCjkti+DVDUufaCXYAw3LVmIuDjntw894WV5Jezb7wBQe7mS0ys/1m1m9mO82s5nonM7vZzHrN7BdmdtTM/szM+MhzIGN+sOKOoOPW7zra4kpQFBTsZnafpC9KWi/pLkl9kraZ2eQqx18p6ZvK/6byYUlPSVot6TNNqBlAwqTpcfs4tLt/IoO9MMp+WNJT7r7S3bdKulPSzyU9UOW0Tym/D82d7r7V3VdJWiNpuZmNbE7pANImq3PttcRx0QsZsU+TdL2kZ4oN7j4gaYukarv8vFdSr7u/WdK2WdJ4Se9qrFQAScaovbaoLY+bedELCfYbC68HytoPSZpqZt1Vzql0fOnXA5AxPLR0odL+CNnyuFl9ExLsYwuvJ8vaTxbOr3QZGlvl+NKvBwCJkLWLTUiwF1eylC/ZLLafq3JOtSWelY4HkBEhDy1lLUiTJiTYXy+8jilrH618SJ+qck758WNK3gOQUZ3+0FIShAT7/sLrlLL2KZLc3SuNzPdXOV6SPLw8AGl0aXfYIyuEe2uEBvsLknqKDYUliwsk9VY5p1fSe82sdP69R9L/SdrdWKkA0mLf6vnBxxLuv9SslUWRwV4Yka+V9AkzW21m8yU9LelySY9LkplNNbNZJac9KeliSVvN7HfNbIWk5ZLWuvuZplQOAG3QaNjGufwz6MlTd39S0lJJ90raJGmcpLnuXlzC+FlJO0uOf0n5tewjCsf/oaQV7v655pUOIMnqCbZOGbXX6pNmXghGhB7o7o9JeqzKex+T9LGyth9Jmj2M2gCk3JG1C4JDe9KyLbGMctt9UWnHv5HdHQG0VBZG7qE3g5OCYAfQcmkP93puBicBwQ4gcZIY7mkSPMcOAMNRz3x7q8T9/duFETuAtqlnSmba8uaF8KQGP5s0rTtWMmIH0FahI/ezueoj7LTP2bcaI3YAbTfckXAcYZ2m0TvBDiAW7Qj3ThytSwQ7gBRLanDHPbon2AHEphkBmNRwjxPBDiD1CPcLsSoGQKyatb6dcP8lRuwAYhf3nHQladyut4gRO4BEqBWISRmNJ+Hp2RCM2AEkXitGwUfWLqj6dZMw6h4ORuwAUqFVo+W0h3gljNgBpEazQrhVYZ6UiwTBDqCjDDd80zB9w1QMgFQZzpRM0kf8zUKwA0idYrCGBnzSg7jZmIoBkGmdFuoSwQ4gxaJCuxNDXSLYAaRcGm5mthtz7ABSr5NDvBJG7ACQMQQ7AGQMwQ4AGUOwA0DGJOHmabckvfzyy3HXAQCpUZKZ3eXvJSHYr5GkhQsXxl0HAKTRNZIOljYkIdifk3SbpJckDcZcCwCkRbfyof5c+RtduVyu/eUAAFqGm6cAkDEEOwBkTBLm2NvCzDZJmi6pv9C00t3/JcaSEsfMHpN0lbvfE3ctSWBmXZI+J+l9haYn3P3JGEtKDDNbIukPJOUk/UjSx939dLxVJYuZvU3SDkk97n6gnd+7Y4Jd0m9Ieqe7vx53IUlkZnMl3SdpW9y1JMidyg8GbpZ0qaT/NrPvuvtP4i0rXmZ2i6Tfl/RuSW9K+oqkT0t6NM66ksTMZkv6kqQb4/j+HTEVY2bXSBor6Z/MbI+Z/bmZdcS/PYSZXSnpEUmr464lSdz9aUkL3P2cpCuUHwj9It6qEuE1Sfe7+yl3z0naI+m6mGtKmo8X/rwYxzfvlHC7StI3JS2UdKuk35S0KNaKEqIw3fD3kj4jqS/mchLH3QfMbK2kn0rqlXQ05pJi5+7/4+7flc4Pmv5Y0tPxVpUs7n6fu++I6/tnairGzD4i6Ymy5lfc/e2SPlxy3BeUnx/8UhvLi1W1vpH0d5J2u/v3zWxq+yuLX8TPjdx9mZk9IulfJX1SUkfMs0f1i5ldL+kbkv7W3b/Z7vriFtU/ceqIdexmNkvS5e7+b4X//pCkj7r7B+OtLH5m9kNJo5R/OGy8pNGSvubun4q1sAQws+mScu7+48J/3y/pZnf/RLyVxc/MZih/oftLd//ruOtJKjM7Ium93DxtjYslrTOz70l6S/lR1z/EW1IyuPu7i383s48p/0PY8aFecLOkT5rZHOWnLT+g/G84Hc3MrpD0rPLz7JvirgdDJTLYzexOSRvcfUxZ+yJJfyppoqTdkpa4+86or+fu3zOzpyT9QPnHcP9Z0oamF94Gze6bLGnBz83XCiPT/1L+N5pN7v7V5lfeWi34mXlA+cUID5rZg4W2Le6+oollt00W/59K3FSMmb1H+dHARe4+uqT9PuVH2Q8rvzfCH0marfwSxsNx1Npu9E119E1l9EttWe2fxIzYzewS5e+uPyLplPLTJ8X3upTv4KfcfWWhbbskV3708Om2F9xG9E119E1l9EttWe+fJC13fJ+k5ZKWSiq/GTNN0vWSnik2uPuApC2S5rWrwBjRN9XRN5XRL7Vlun+SFOzPSZrs7l9Q/jHlUsWnt8rvLB+SNNXMhmw0nzH0TXX0TWX0S22Z7p/ETMW4+/Eab48tvJ4saz+p/MXpMklvtKKuJKBvqqNvKqNfast6/yRpxF5LV+G1/MpabD/XxlqShr6pjr6pjH6pLfX9k5ZgL27cNaasfbTynXyqveUkCn1THX1TGf1SW+r7Jy3Bvr/wOqWsfYokL2xE1Knom+rom8rol9pS3z9pCvYXJPUUG8xspKQFym/M1Mnom+rom8rol9pS3z+JuXlai7vnCjvsPWFmrym/ef39ki6X9HisxcWMvqmOvqmMfqktC/2TlhG7Cp9cs1TSvZI2SRonaa67H4q1sASgb6qjbyqjX2pLe/8kbksBAMDwpGbEDgAIQ7ADQMYQ7ACQMQQ7AGQMwQ4AGUOwA0DGEOwAkDEEOwBkDMEOABnz/8HA9fNQrleeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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[:20000]\n",
    "t = np.arange(autocorr.size)*delta_t\n",
    "plt.semilogx(t, autocorr,\"o\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "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": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Model]]\n",
      "    Model(fitcorr)\n",
      "[[Fit Statistics]]\n",
      "    # fitting method   = leastsq\n",
      "    # function evals   = 27\n",
      "    # data points      = 20000\n",
      "    # variables        = 1\n",
      "    chi-square         = 13.0893610\n",
      "    reduced chi-square = 6.5450e-04\n",
      "    Akaike info crit   = -146631.756\n",
      "    Bayesian info crit = -146623.852\n",
      "[[Variables]]\n",
      "    kappa:  10 (fixed)\n",
      "    sigma:  125 (fixed)\n",
      "    ga:     25 (fixed)\n",
      "    eta:    0.27417142 +/- 0.00229336 (0.84%) (init = 0.25)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAECCAYAAADq7fyyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f3/8ddkkrATBIki+SL70W9ZW2QRFfyBimAxpVpsoValWL5VW0RRKCqbfEXcrfrTutQFqqi/igLKImojshUFg6hHVtkpWwiyZZvfH5OJWSYzN8lk1vfz8chjwl1mPlzxfc+ce+65Lo/Hg4iIxI+kSBcgIiKhpWAXEYkzCnYRkTijYBcRiTPJkS7AGFMHuADYCxRGuBwRkVjhBloA/7bWni69IuLBjjfUP410ESIiMepiYHnpBdEQ7HsB5syZw9lnnx3pWkREYsK+ffsYMWIEFGdoadEQ7IUAZ599NhkZGZGuRUQk1lTowtbFUxGROKNgFxGJMwp2EZE4o2AXEYkz0XDxtFrmrdvNQ4ste3JOck6Teoy/wpDZvWWkyxIRibiYDPZ563Yz8Z8bOJnvvRi8O+ckE/+5AaBG4V6dk4WTfQJt428d4GiZTmQi4o8r0tP2GmNaA9uWLVvmeLhj35kfsTvnZIXlLZvU47MJ/6dadZQ/WQDUS3HzwLDOlQaok30CbQNUWJeS5AIX5Bd6gi5rWDeZnBP5nNOkHpee15yPvz3AnpyTpNVLweWCnBP5lf6uk4NIbNu1axcDBgwAaGOt3V56XUy22Pf4CfVAy514aLEtE7AAJ/MLeWixrTT8nOwTaBvf76XlF1U80Va27MiJfMD7jWX2qh0l63JO5gf9fXfOScbOXc+4N9dT5AG3y0Whx1Ppa0udCERiRkwG+zlN6vltsZ/TpF6137M6Jwsn+9TGSSiUfOeMwuJvbpW9+k4EY+euL7N/ksv7Hgp+kegRk8E+/grjt3vD1xddHdU5WTjZJ9g2/tbFEt+JobLgBzijfgqTf/4Thb5ImMTkcMfM7i15YFhnWjaphwtvazFQX7gT468w1Etxl1kW7GThZJ9A2/hbl5LkIsXtCroslhw5kc/YuetpPWFhyU/3aUuYt253pEsTiUsx2WIHb7iHsgVYuk/c6cgTJ/s42aaqo2LS6qVwPK+gzMXUWOMLe18LX616kdCJyVExUnGYpNNRMXkFhZzIL4p0+UE1SHUz4xc1+xYmEs/iblSM1Owbi++ksDvnZNDRMC4gEqf+43mFJS36JBf8plcr7s/sHIFKRGKPgj0BVfWkMG/dbqa8t7HMcMlwKvLA7FU7SoZ0qjUvEpiCXYIKdCKIROiXbs2rb16kIgW71EikQ7/0RdiRvdVdIwIKdqlF5UO/toPe113jAkYo5CWBKdglbPy17u+Zt4E5q3aE9AKthx9DXq14SUQKdomo+zM7lwneUAe9L+D7tmvKnNF9QvSuItEtJu88lfh1f2Znts0cwvaZQ3h8eDfqpYTmn+hnWw7TesJCRjy/MiTvJxLN1GKXqFW668Y7/XE2J2t4c5Uv4NWCl3imFrvEhMzuLflm+pVsnzmEkb1b1fj9fAF/z7wNIahOJLoo2CXm3J/Zme3F3TU1DfnZq3bQZsJCTUgmcUXBLjHNF/KPD+9GdbvjPcDYuevpNWNpSGsTiRQFu8SFzO4t2fS/Nbvguv9Yni6wSlxQsEtcKd0X37dd02q9h6//Xd0zEqsU7BK35ozuU6OAV/eMxCoFu8Q9X8B3SG9Q5X193TNqvUssUbBLwlg6rj+PD+9WrX3Hzl2vvneJGQp2SSiZ3VtWe5jkZ1sOc96k92uhKpHQUrBLQvINk6xq98ypQo9ubJKop2CXhFbd7pnZq3bowqpELQW7JDxf90xVW+/7j+Wpa0aikoJdpFh1Wu++rhmRaKJgFynF13pvXMddpf00JFKiiYJdxI/sqYOqPHJGQyIlWijYRSrhGzmT7HK+z2dbDnPZo5/UWk0iTijYRYLY/EDVumY2/ee4RsxIRCnYRRzInjqoSqNmNGJGIsnl8YTy+fBVZ4xpDWxLT08nOTmZSZMmMXDgQNavX8/tt99eYfsZM2Zw4YUXsmLFCiZNmgSA2+0mKSmJpKQkHnzwQbp27cqnn37KY489VrLO9zp9+nTatGlDVlYWc+bMwe12l9lmwoQJpKens3z5chYtWlRmX7fbzW233UajRo1YuXIlq1evLlmenJxMamoqI0aMIDU1la+//podO3aQmppKamoqKSkppKam0q1bN1wuF0eOHOH06dMl633buFxV+N4vYTdv3W7Gzl3vePtkl7fFLxJqu3btYsCAAQBtrLXbS6+LmmeeFhUVUVT04/MsPR5PmT+XXl56vcfjIT8/n6KiIgoLCyksLAQgNzeXTZs2lSz3vR4/fhyA7du38+6775ZZV1RUxB//+EfS09NZs2YNM2fOLHk/n5tuuolGjRqxaNEipk2bVqG+X/3qV6SmpvLCCy/w2GOPlVnncrlK3u/uu+/m+eefL7O+YcOGHDt2DIDRo0ezYMEC6tevT/369WnQoAEZGRm8/fbbADz55JNYa8usP/vss7n++usByM7OpqCggLS0tJKflJQUh/81pDK+57A6HeJY4IH2Excq3CWsoqbFvmzZMjIyMiJaS2U8Hk9J8Pta1adOneLUqVMlJ5OCggLy8vJo1aoVSUlJ7Nixgz179pCXl1fyk5+fz9VXXw1AVlYWX3/9dZn1SUlJTJgwAYAXX3yR1atXc+LEiZKftLQ05s6dC8CIESNYsmRJyTqATp06sWGD91b3Pn36sGrVqjJ/j4svvpisrCwArr/+eg4cOEDz5s1JT0+nefPmdOrUiSFDvAG0c+dO0tLSaNSokb5FVKLNhIU4/b9HLXcJtUAtdgV7HCgqKuLUqVPk5eXRpEkTANasWcOePXs4evRoyU96ejp/+MMfALjhhhv46quvOHDgAAcOHODkyZNcffXVzJs3D4Czzz6b/fv3k5qaSosWLcjIyGDYsGGMGzcOgPnz53PWWWeRkZHBWWedhdtdtXHf8ULhLpESkq4YY8xo4C4gA1gPjLPWVjpo1xhzITAD6A6cAD4Exltr91f1LyCBJSUllXTJ+PTs2TPgPi+//HKZPx8/fpzTp0+X/Pnhhx9m3759HDhwgL1797Jr1y5OnjwJwOnTpxk6dGjJtsnJyZx77rmMHTuWW2+9lYKCAhYvXkz79u1p06YNqampIfhbRqdtM4fQfuJCChyku7plJFwcBbsx5nrgWWAa8G/gNmCxMaartXabn+3PB5YBS4FfA2cA04v3ucBamx+i+iVEGjRoQIMGP476GDlyZKXbJicns27dOnbu3MmuXbvYuXMnW7dupVmzZgB8//33XHXVVYD3pNOhQwc6d+7MrbfeSr9+/UquM8RLK3/zA0PoMnkRuacLg25b4IHzJr3PtzMGh6EySVRBg90Y48Ib6H+z1k4tXrYUsMDtwJ/87HYrsBf4pS/EjTGbgDXAZYDGgcUwt9tNt27d6NbN/7wqLVu2ZMWKFWzevJnvvvuOjRs3sn79eg4fPgzA8uXLufLKK+nUqRM9e/akV69e9OzZkw4dOpCUFJsjcLOnDuK8Se9zqjB40/1UoYdeM5ayetJlYahMEpGTFnt74FzgPd8Ca22+MWYhMKiSfTYCX5drmdvi1zbVKVRiR926denTpw99+vTxuz49PZ0xY8awbt06XnnlFZ5++mkAVqxYQZ8+ffj222/Zs2cPffr0oV69euEsvUa+nTHYcbjvP5bHPfM2cH9m5zBUJonGSfOoY/Hr5nLLtwLtjDEVvk9ba5+x1j5dbvHPi1+/rVqJEm/OP/98Hn30UT7++GNycnLYsGEDL774It27dwfghRdeYMCAAZxxxhn079+fqVOnkpWVVWHoaTT6dsZg6rqdjSKavWpHLVcjicpJsDcufj1Wbvmx4v2D3o5njPkv4GFgLfBRVQqU+OZ2u+nUqRM33XQTdevWBeDee+9l/vz53HLLLeTm5jJ16lSGDRtWss/atWs5ePBgpEoO6tsZgx3PL6Mpf6U2OAl23z/R8t8vfcsr3kVUSnGoLyv+rOustZEdXylRLy0tjauuuopHHnmEL774gkOHDrFo0SLcbjcej4fhw4eTnp5O3759eeihh9ixI/pavpsfcD55mMJdQs1JsB8tfm1UbnlDvKF+vLIdjTGdgBV4W/2XWWu3VKdISWxnnHEGPXr0KPnz3Llzue+++zh16hR33XUX5557LlOmTIlcgZXY/MAQnN7a1WXyolqtRRKLk2DfVPzattzytoCtrAVujOkFZAGFwMXW2uxqVylSzOVy0aNHD6ZMmcLnn3/O5s2buf/+++nXrx8A1lquuuoq5s+fT0FBQYSr9Y5zdyL3dKHmcpeQcRrsO4FM3wJjTAowBG8XSwXFd5N+AOwHLrTWbvK3nUhNtWvXjkmTJnHppZcC3qkQPv/8c4YOHUrr1q2ZPHlyxLtqtjsM98+2HNZTmCQkggZ7cYt8JjDGGDPDGDMYeBc4E3gMwBjTzhjTu9RuT+DtfpkGtDLG9C710yLkfwuRYgMHDmTHjh288847dO7cmenTp3P++eeTm5sb0bqcPku1KjNHilTG0d0g1tpngPHAb4G3gSbAFdbarcWb3AushJLW/GDADfyjeHnpnxEhrF+kgpSUFDIzM/nggw/YunUrr7zyCo0bewd3jR07lvfee8/vzKG1KbN7S8fzuetiqtSUJgGThHHo0CF++tOfsmPHDrp06cK0adMYOnRoWGevrMqkYU67cCQxBZoELDbv3xaphmbNmrFlyxZee+01Tpw4QWZmJj179mTTpvBdAnJ6MRU0UkaqT8EuCSU5OZmRI0fyzTff8NJLL+FyuWjRwnvZxzd7ZW1z2hLPPV3IPfM21HI1Eo8U7JKQkpOTufHGG1mzZg0NGzYkPz+fHj16MGrUKPbt21frn+/0YqqmHZDqULCLAPn5+Vx55ZW89tprdOzYkUceeaRWx8HrYqrUJgW7CFC/fn0efvhhNm7cyCWXXMKdd95Jr1692Lt3b6195tJx/R1PO3DZo5/UWh0SfxTsIqV06NCB+fPn89Zbb9GiRQuaN29eq5/n9GlKm/5T6cwdIhUo2EXKcblcXHPNNSxYsIDk5GQOHTpEv379+OKLL2rl85xeTFWXjDilYBcJYteuXWzZsoXevXsza9asWrm5yenFVHXJiBMKdpEgunbtSnZ2NkOHDuXuu+9m4MCB7Nq1K6Sfkdm9JY3rBH8GrLpkxAkFu4gDTZs25a233uKll15izZo13HXXXSH/jOyplT1psix1yUgwCnYRh1wuFzfeeCNr167lr3/9KwC5ubmEcloOp/3tunFJAlGwi1TReeedR7NmzcjPz2fQoEFce+21/PDDDyF7fyfj23XjkgSiYBeppuTkZIYNG8Y777zDxRdfzO7doZlLfem4/o6201wyUhkFu0g1uVwu7rzzThYsWMDmzZvp1asXX375ZUje20mXTO7pwpB8lsQfBbtIDV155ZUsX74cgBtuuCFkwyGdjJLRhVTxR8EuEgJdu3Zl9erVvPnmmyQlJYXkgqrTUTK6kCrlKdhFQqRly5Z06NCBoqIibr75Zp566qkav2ffdk2DbqMLqVKegl0kxAoKCjh48CC33XYbM2bMqFHrfc7oPo62G/H8ymp/hsQfBbtIiKWmpvLWW28xcuRI7rnnHiZMmFCjcHdyIfWzLYer/f4SfxTsIrUgOTmZV155hTFjxjBr1iwmT55co/dzciFV88iIT3KkCxCJV0lJSTzzzDPUqVOHiy66qEbvlT11UNARMJpHRnzUYhepRS6Xi8cff5zLL78cgA0bqj+C5axGqUG36TVjabXfX+KHgl0kTD766CO6du3Ko48+Wq39V0+6LOg2+4/lVeu9Jb4o2EXC5JJLLuGXv/wld9xxB7Nnz67WeziZR0atdlGwi4RJcnIyr732Gv379+fGG29kyZIlVX4PJ/PIqNUuCnaRMKpbty7z5s3jJz/5CcOGDavWAzuctNo1QiaxKdhFwiwtLY0PPviARx99lIyMjCrv76TVrhEyiU3BLhIBLVq04OabbwZg48aN5OTkVGl/J6123Y2auBTsIhGUm5vLJZdcwnXXXUdBQYHj/Zy02nU3auJSsItEUOPGjZk1axaLFy/mjjvuqNK+Tsa1a+bHxKRgF4mwUaNGMW7cOJ588kmee+45x/s5GdeumR8Tk4JdJArMmjWLwYMHc+utt7J69WrH+zlptUviUbCLRAG3283rr7/OvffeS7du3Rzv56TV3kZPWUo4CnaRKNG4cWPuu+8+6tSpw5EjR8jPz3e0X123K+D6mj/LSWKNgl0kyhw9epSf/exnjB8/3tH2384YHHSb8ya9X9OyJIYo2EWiTFpaGkOHDuWJJ57gjTfecLRP4DY7nCpUuz2RKNhFotBDDz1E3759+f3vf88333wTdPttDp6ypMnBEoeCXSQKpaSk8Oabb1KvXj1+/etfc+rUqRq/pyYHSxwKdpEodc455/DKK6/QrFkzfvjhh6Dbj+zdKug2mmYgMSjYRaLY4MGD+fDDDznzzDODbnt/Zueg22iagcSgYBeJci6Xi3379jF8+HB2794dcFtNMyCgYBeJCUePHmXhwoWMHDmSwsLCSrfTNAMCCnaRmGCM4amnnuKTTz5h1qxZAbcNdsOSxD8Fu0iM+N3vfse1117L5MmT+fLLLyvdzskNS601zUBcU7CLxAiXy8UzzzxDs2bN+Mtf/hJ42zDVJNFJwS4SQ84880zmz5/PnDlzAm7n5IYltdrjl4JdJMb06NGDJk2akJeXx7Zt2yJdjkQhBbtIjBo+fDiXX345x4/7f3D1dgetdk3pG58U7CIx6s9//jObN2/m7rvvrvZ7aGqw+KRgF4lR/fv3Z+zYsTz99NN88sknfrd5fHjwh3a0n6hWe7xRsIvEsBkzZtC2bVtGjx7NiRMnKqzP7N4y6HsUqNkedxTsIjGsfv36vPDCC9SvX5/9+/f73aZvu6ZB30dT+saXZKcbGmNGA3cBGcB6YJy1NuhUccaYRsBXwB3W2rerW6iI+HfppZeybt06kpL8t9PmjO4TdGijpvSNL45a7MaY64FngdnAL4EcYLExpk2Q/RoB7wLB5xMVkWpLSkri6NGjTJ06lby8iiHdIb1B0PeYty7wBGMSO4IGuzHGBUwD/matnWqtfR8YChwEbg+wXz9gDeD8kesiUm3Lly9nypQpPPjggxXWLR3XP+j+Y+eur4WqJBKctNjbA+cC7/kWWGvzgYXAoAD7zQM2BNlGREJkyJAhXHfddUyfPp2vv/66wnpNDpY4nAR7x+LXzeWWbwXaGWPclex3sbX2V8B/qluciFTNE088QePGjRk1alSF6X2dTA6mi6jxwUmwNy5+PVZu+bHi/f123llrv6pBXSJSDenp6TzxxBOsWrWKZ599tsr76yJqfHAyKsb3/a38aFff8qLQlSMiNfWb3/yGTZs2MWRIxSkFts8cosm/EoCTFvvR4tdG5ZY3xBvq/ieqEJGIcLlcTJkyhdatW+PxVP3uI92JGvucBPum4te25Za3Bay1VvetiUShgwcPMnjwYBYuLBvUwZ6LqjtRY5/TYN8JZPoWGGNSgCHAslqqS0RqqHHjxuzYsYNbbrmlzAyQTp6LKrEtaLAXt8hnAmOMMTOMMYPx3nR0JvAYgDGmnTGmd61WKiJVkpqayrPPPsv333/PtGnTqrSv+uFjm6M7T621zwDjgd8CbwNNgCustVuLN7kXCDq9gIiE18UXX8yoUaN45JFHyM7OLlnuZP4YiV2u6lxcCSVjTGtg27Jly8jIyIhoLSLx6NChQ5x33nn07duXefPmlSwP1irv264pc0b3qe3ypJp27drFgAEDANpYa7eXXud4EjARiU3NmjVj/vz5nH/++VXa77Mth2upIqltmrZXJAH07t2btLQ08vPzOXrUO4LZyaPzJDYp2EUSRGFhIRdeeCFjxoxxvE+XyYtqsSKpLQp2kQThdrv5+c9/zhtvvMGiRd7ATg4yL1ju6cLAG0hUUrCLJJC7774bYwx//OMfOXHiBJsfUHdMPFKwiySQOnXq8Nxzz7Ft2zamT5/uaB/N+Bh7FOwiCaZfv37ceOONLF26lPz8fBrXqWzmbS/N+Bh7FOwiCeiJJ55g5cqVpKSkkD1Vz8KJNwp2kQTUqFEjUlJSyMnJISsrK+j298zbEIaqJFQU7CIJbMyYMQwdOpQz+CHgdrNX7QhTRRIKCnaRBDZt2jROnjzJed/PC76xxAwFu0gC69ixI5MmTeKNN97g5NbPI12OhIiCXSTB+ca2H17yDEX5pyrdTnehxg4Fu0iC841tv+KiHnjyKg923YUaOxTsIkK/fv1YsGAB7gZNIl2KhICCXURK5Ofs48jHL+HxFPldP2/d7jBXJNWhYBeREq3yvid3zT/54cvFftePnbs+zBVJdSjYRaTEihenUqdVF4588jKFPxyJdDlSTQp2ESnhcrlodsUteApOc/ij5yNdjlSTgl1Eykhp2pK0PsM58U0WJ7esrbC+TZBnpUrk6ZmnIlLG9plDOHd8Png81Mn47wrrPRGoSapGwS4iFbjcKTS56DcAeIoKcSUFntpXoou6YkSkUgVH/8Pev/+JE1v+XWZ5+4nqjolmCnYRqaBvu6YAuBucAcDhRX+l6NSPM0AWqD8mqinYRaSCOaP7AOBKTqHZ4LEUHs/h8EcvRLgqcUrBLiIB1WnRgca9r+H4hg/LjJI5b9L7EaxKAtHFUxHxq0N6Azb95zgATS78NSc3reLY+g+o164HAKcK1R8TSGs/w0K3zxwSls9Wi11E/Fo6rn/J767kFNKvnUrzzAmRKyiG+Av1QMtDTcEuIo4kN26Oy51C4cljnNy2DghfUMWT1hMWlvzUFgW7iFSqrttVYdmRj57nwD/vJ//QzghUFF9qK9wV7CJSqW9nDK6wrEm/G3Cl1OHg/IfxFObrImoUUrCLSJUkN2xKsyv/TN7+LRz55GVdRI1CCnYRCWhk71YVltXv0ItGP/s5x9a+ywm7IgJVSSAKdhEJ6P7Mzn6Xn3HpTTTqcTV1Ms7XRdQoo3HsIlItLncKTQeMBrwThZ04cYL69etHuCoBtdhFxIFAN9Z4PEX85837aP7TKygq8v+sVAkvBbuI1IjLlUTdNj/lhF3OlClTIl1OxEVDt5SCXUQc8XcR1adxz2E07HI506dP59VXXw1jVbGvNk4ECnYRcaSyi6jgfVZq08v/h7rnduH3v/89WVlZYaxMylOwi4hjHdIbVLrO5U7hzMy/4G7ZiSZNmoSxKilPwS4ijpWeGMwfd92GnDV8Ol26dMHj8XDw4MHwFCZlaLijiFTJ9plDgvYLd5m8iF+wkr///e8sX76cVq0q75+PRpGccjcU1GIXkZDLPV1IZmYmubm59O/fn61bt0a6JMcCTbkbDSNenFCwi0iV+Z6JGkj37t1ZunQpR48e5aKLLmLjxo1hqKz2xUK4K9hFpMp8z0QNpPWEhVxwwQX861//AmDgwIEcP368tkurEaeh3SbCD9IIRn3sIlKrOnXqxKeffspXX31FgwaVj6qJJdE+n6Va7CJSLU4uJvpasO3atePqq68G4B//+AcTJkygoKCgVusLt2hprYOCXURq2YjnV5b585o1a3jwwQfp378/33//fYSqqrnS3THRFOqgYBeRGnh8eLeg23y25XDZfR5/nNmzZ5OdnU3Xrl159dVX8Xgi37lR1XD24A33aAt1ULCLSA1kdm/paLvy4TdixAi+/PJLOnXqxO9+9ztWr15dG+XVusifjvxTsItIjZzVKNXRdvfM21Dmz23atCErK4slS5bQu3dvAF5//XV27ky8h2SHutXvONiNMaONMZuMMSeNMSuNMQHHOxljOhljlhljfjDG7DDG3G2MqfjIcxGJaasnXeZou9mrdlRYlpSUxGWXefc/evQoN998M+3bt2fMmDFs2LChwvbijKNgN8ZcDzwLzAZ+CeQAi40xbSrZPh34EO83lV8BfwNmAHeEoGYRiTKhuN0+LS2NjRs3csMNN/Dyyy/TpUsXLrnkErKzs0NQYWSFezqCoMFe3MqeBvzNWjvVWvs+MBQ4CNxeyW634B0jP9Ra+7619n7gAWCiMSYlNKWLSKwJ1uXQqlUrnnvuOXbv3s1DDz3E3r17adiwIQAff/wxM2fOZNWqVeTn54ej3JCIxBwzTlrs7YFzgfd8C6y1+cBCYFAl+wwElllrT5RaNg9oClxQvVJFJJqFMsCaNWvGnXfeyXfffUfbtm0BeO+995g4cSJ9+vQhLS2Nnj17cvPNN5eMqNm1axeHDx+OihE2/gSa8hhC28/u5M7TjsWvm8st3wq0M8a4rbWFfvb5xM/2vnUrqlKkiMQGJzM/tp6w0PFJwOX68bLcY489xsSJE8nKyuKzzz4jOzubjRs3lmwzatQolixZQr169WjatClNmjShW7duzJ49G4CpU6eya9cuUlNTqVOnDqmpqXTs2JGbbroJgNzP51N06gffB+PCRfIZ59Dg/IsBOPbFAoryTwGu4h9IaZZB/fY9vfuvfRdPYdkonDtxeMnvVxStZe3qb8qsr3NOR+r+V6dqHZtAnAR74+LXY+WWH8Pb4m8A5PrZx9/2pd9PRKRK0tPTueaaa7jmmmsqrLvzzjsZNGgQu3fvJicnh5ycHJo2/XGysjVr1rB+/XpOnz5NXl4ep0+f5sILL+Smm26i9YSFHPt8PgVH9pR5z3rtLigJ9qMr5lJ4/EiZ9fXP71cS7DmfzsaTd7LM+rc7uLn88ssBGD9+fIWaG/e6pkywh4qTYPedMst/v/Et9/dYcpef7X30GHOROPb48G6Mnbs+4DahapmWdtlll5WMsPFn4cLA3yTOGf0s+LpxPB68EfbjN4Zz/vB8car9GG2uJHfJ7xm3VHzW65P/e1XJ78eOHeO/71tUZr0rqXam63LSx360+LVRueUN8Ya0v+najvrZvlGpdSISp6p701KkuVxJuJLc3h93Mi53Ci73j8GblFKXpNS6JKXWK/lxJf84hs82TkAAAAclSURBVL/0ct9P3bp1S9Y3bNiwwnpXcu2MJXES7JuKX9uWW94WsNZafy3zTZVsD2Cdlycisaiu29ktK9EW7vHCabDvBDJ9C4qHLA4BllWyzzJgoDGm9GXgTOAQEPg7mojEvG9nDHa8rcL9R6Hqngoa7MUt8pnAGGPMDGPMYOBd4EzgMQBjTDtjTO9Suz0DpALvG2OuMsZMAiYCM621eSGpXEQkDKobtpF8RqqjO0+ttc8A44HfAm8DTYArrLW+IYz3AitLbb8X71j25OLtbwYmWWsfDl3pIhLNqhJsidJqD3RMQnkicHxJ1lr7CPBIJetuAG4ot2wt0LcGtYlIjHMyrt2nNkbKOP3ccArH31GzO4pIrYqHlrvTi8HRQsEuIrUu1sO9KheDo4GCXUSiTjSGeyypndueRETKqUp/e22J9OeHi1rsIhI2VemSaT8xdCHcuprPJo3kkMWaUItdRMLKacu9wFN5CzvW++xrm1rsIhJ2NW0JRyKsY6n1rmAXkYgIR7gnYmsdFOwiEsOiNbgj3bpXsItIxIQiAKM13CNJwS4iMU/hXpZGxYhIRIVqfLvC/UdqsYtIxEW6T9qfWJyu10ctdhGJCoECMVpa49Fw96wTarGLSNSrjVbw9plDKn3faGh114Ra7CISE2qrtRzrIe6PWuwiEjNCFcK1FebRcpJQsItIQqlp+MZC9426YkQkptSkSybaW/yhomAXkZjjC1anAR/tQRxq6ooRkbiWaKEOCnYRiWHBQjsRQx0U7CIS42LhYma4qY9dRGJeIoe4P2qxi4jEGQW7iEicUbCLiMQZBbuISJyJhounboB9+/ZFug4RkZhRKjPd5ddFQ7C3ABgxYkSk6xARiUUtgC2lF0RDsP8buBjYCxRGuBYRkVjhxhvq/y6/wuXxeMJfjoiI1BpdPBURiTMKdhGROBMNfexhYYx5G+gMnCxeNNVa+04ES4o6xphHgLOstSMjXUs0MMa4gIeBK4sXPWWtfSaCJUUNY8w4YBTgAdYCf7DWno5sVdHFGJMGfAZkWms3h/OzEybYgR5AV2vt0UgXEo2MMVcA1wOLI11LFBmKtzHQCagLfG2M+Ze1dmNky4osY0wv4EagJ3ACeBX4E/BQJOuKJsaYvsBzQMdIfH5CdMUYY1oAjYG3jDHZxpgpxpiE+Ls7YYxJB6YDMyJdSzSx1r4LDLHWFgHN8TaEfohsVVHhCHCrtfa4tdYDZAOtIlxTtPlD8c+eSHx4ooTbWcCHwAigD3AJMDqiFUWJ4u6Gl4A7gJwIlxN1rLX5xpiZwDfAMmBHhEuKOGvtd9baf0FJo+nPwLuRrSq6WGuvt9Z+FqnPj6uuGGPMdcBT5Rbvt9b+BPhVqe2exNs/+FwYy4uoyo4N8CKw3lr7qTGmXfgri7wg/26w1k4wxkwH5gP/AyREP3uw42KMORf4APi/1toPw11fpAU7PpGUEOPYjTG9gTOttQuK/3wN8Btr7bDIVhZ5xpg1QD28N4c1BRoCr1trb4loYVHAGNMZ8Fhrvyr+861AJ2vtmMhWFnnGmG54T3SzrLV/jXQ90coYsx0YqIuntSMVeNwYkwWcwtvq+ntkS4oO1tqevt+NMTfg/UeY8KFerBPwP8aYAXi7LX+B9xtOQjPGNAcW4e1nfzvS9UhFURnsxpihwBxrbaNyy0cDdwEZwHpgnLV2ZbD3s9ZmGWP+BqzGexvu/wPmhLzwMAj1sYkntfDv5vXilumXeL/RvG2t/UfoK69dtfBv5na8gxHuMcbcU7xsobV2UgjLDpt4/H8q6rpijDEX4m0NJFlrG5Zafj3eVvY0vHMj3Ab0xTuEcVskag03HZvK6dj4p+MSWLwen6hpsRtj6uC9uj4dOI63+8S3zoX3AP/NWju1eNlSwOJtPfwp7AWHkY5N5XRs/NNxCSzej080DXe8EpgIjAfKX4xpD5wLvOdbYK3NBxYCg8JVYATp2FROx8Y/HZfA4vr4RFOw/xtoY619Eu9tyqX57t4qf2V5K9DOGFNhovk4o2NTOR0b/3RcAovr4xM1XTHW2t0BVjcufj1WbvkxvCenBkBubdQVDXRsKqdj45+OS2DxfnyiqcUeiKv4tfyZ1be8KIy1RBsdm8rp2Pin4xJYzB+fWAl238Rdjcotb4j3IB8PbzlRRcemcjo2/um4BBbzxydWgn1T8WvbcsvbArZ4IqJEpWNTOR0b/3RcAov54xNLwb4TyPQtMMakAEPwTsyUyHRsKqdj45+OS2Axf3yi5uJpINZaT/EMe08ZY47gnbz+VuBM4LGIFhdhOjaV07HxT8clsHg4PrHSYqf4yTXjgd8CbwNNgCustVsjWlgU0LGpnI6NfzougcX68Ym6KQVERKRmYqbFLiIizijYRUTijIJdRCTOKNhFROKMgl1EJM4o2EVE4oyCXUQkzijYRUTijIJdRCTO/H9e5+sZp5hAMwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done with.. 0\n",
      "done with.. 1\n",
      "done with.. 2\n",
      "done with.. 3\n",
      "done with.. 4\n",
      "done with.. 5\n",
      "done with.. 6\n",
      "done with.. 7\n",
      "done with.. 8\n",
      "done with.. 9\n",
      "done with.. 10\n",
      "done with.. 11\n",
      "done with.. 12\n",
      "done with.. 13\n",
      "done with.. 14\n",
      "done with.. 15\n",
      "done with.. 16\n",
      "done with.. 17\n",
      "done with.. 18\n",
      "done with.. 19\n",
      "done with.. 20\n",
      "done with.. 21\n",
      "done with.. 22\n",
      "done with.. 23\n",
      "done with.. 24\n",
      "done with.. 25\n",
      "done with.. 26\n",
      "done with.. 27\n",
      "done with.. 28\n",
      "done with.. 29\n",
      "done with.. 30\n",
      "done with.. 31\n",
      "done with.. 32\n",
      "done with.. 33\n",
      "done with.. 34\n",
      "done with.. 35\n",
      "done with.. 36\n",
      "done with.. 37\n",
      "done with.. 38\n",
      "done with.. 39\n",
      "done with.. 40\n",
      "done with.. 41\n",
      "done with.. 42\n",
      "done with.. 43\n",
      "done with.. 44\n",
      "done with.. 45\n",
      "done with.. 46\n",
      "done with.. 47\n",
      "done with.. 48\n",
      "done with.. 49\n",
      "done with.. 50\n",
      "done with.. 51\n",
      "done with.. 52\n",
      "done with.. 53\n",
      "done with.. 54\n",
      "done with.. 55\n",
      "done with.. 56\n",
      "done with.. 57\n",
      "done with.. 58\n",
      "done with.. 59\n",
      "done with.. 60\n",
      "done with.. 61\n",
      "done with.. 62\n",
      "done with.. 63\n",
      "done with.. 64\n",
      "done with.. 65\n",
      "done with.. 66\n",
      "done with.. 67\n",
      "done with.. 68\n",
      "done with.. 69\n",
      "done with.. 70\n",
      "done with.. 71\n",
      "done with.. 72\n",
      "done with.. 73\n",
      "done with.. 74\n",
      "done with.. 75\n",
      "done with.. 76\n",
      "done with.. 77\n",
      "done with.. 78\n",
      "done with.. 79\n",
      "done with.. 80\n",
      "done with.. 81\n",
      "done with.. 82\n",
      "done with.. 83\n",
      "done with.. 84\n",
      "done with.. 85\n",
      "done with.. 86\n",
      "done with.. 87\n",
      "done with.. 88\n",
      "done with.. 89\n",
      "done with.. 90\n",
      "done with.. 91\n",
      "done with.. 92\n",
      "done with.. 93\n",
      "done with.. 94\n",
      "done with.. 95\n",
      "done with.. 96\n",
      "done with.. 97\n",
      "done with.. 98\n",
      "done with.. 99\n"
     ]
    }
   ],
   "source": [
    "gamma_list = []\n",
    "gamma_std_list = []\n",
    "N = 100000\n",
    "for i in range(100):\n",
    "    ts=simulate_membrane(N)\n",
    "    # calculate autocorrelation function\n",
    "    f = np.fft.rfft(ts)\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[:20000]\n",
    "    t = np.arange(autocorr.size)*delta_t\n",
    "\n",
    "    # fit the autocorrelation function\n",
    "    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=True)\n",
    "    paramsFC['eta'].set(value=eta,vary=False)\n",
    "\n",
    "    resultFC = fit_corr_model.fit(t=t,data=autocorr,params=paramsFC)\n",
    "\n",
    "    gamma_list.append(np.abs(resultFC.values['ga']))\n",
    "    gamma_std_list.append(resultFC.covar[0][0])\n",
    "    print(\"done with..\",i)\n",
    "    \n",
    "gamma_list = np.array(gamma_list)\n",
    "gamma_std_list = np.array(gamma_std_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/hstrey/anaconda3/envs/pymc3/lib/python3.6/site-packages/ipykernel_launcher.py:1: MatplotlibDeprecationWarning: \n",
      "The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.03, 0.18)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbl0lEQVR4nO3deZhkdX3v8XfD6AQdRJRRlhFZlK+gDEswGBVRwAf0Cg5CEr1j1IgLeRRHBRQQlEVQEZQdRGOIihHiBUTlIjgR2XMxBgHFLyggjjKgyC6y9v3jdwrOFNU93dWnuqp63q/n6aeYs9X3x+mqT//O7ywjo6OjSJLUhJX6XYAkaeYwVCRJjTFUJEmNMVQkSY0xVCRJjZnV7wJ6ISJmAy8HbgMe63M5kjQsVgbWAq7KzIe62cCMDBVKoFzS7yIkaUhtA1zazYozNVRuAzj99NNZc801+12LJA2FpUuXsnDhQqi+Q7sxU0PlMYA111yTefPm9bsWSRo2XQ8bOFAvSWqMoSJJaoyhIklqjKEiSWqMoSJJasxAnf0VESsDi4D3AusCvwFOAk7MzNGI2Aq4qsOqR2fmPtNXqSSpk4EKFeAgYD/gMOBKygU4xwDPAI4E5gMPADu0rff7aaxRkjSGgQmViFgJ+Cjw+cw8vJq8OCLmAvvwZKhcl5lX9qlMzSAjIyNT3oYPuZOWNUhjKqsBXwPOapuewNyIeCYlVK6Z7sIkSRMzMD2VzLwL+GCHWTsDSzLzgYjYFHgoIq4GNgFuBQ7LzH+bxlIlSWMYpJ7KU0TEeyjjJ0dGxNrAGsCLgcOBNwI/Bk6LiHf0r0pJUsvA9FTaRcRC4BTg28AJwCrATsA1mdm62dkPq7D5FOXQmSSpjwYyVCLiI8DRwLnAwswcBf4M/KDD4ucDO0XEnMy8fxrLlCS1GbhQiYgjgP0pPY89MvPRavpGwHbAv7Y9PGYV4EHKqcaSpD4aqDGViFhECZRjgXe1AqWyDnAyZSyltfwI8Bbgkqo3I0nqo4HpqUTEWsDngGuBbwFbR0R9kcspTyI7JSJWpzxE5v2U04xfPb3VSpI6GZhQAXYEZgObAld0mD8XeDNwBHAo8Fzgp8DrM/Mn01WkJGlsAxMqmXkacNoEFt2zt5VIkro1UGMqkqThZqhIkhpjqEiSGmOoSJIaY6hIkhpjqEiSGmOoSJIaY6hIkhpjqEiSGmOoSJIaY6hIkhpjqEiSGmOoSJIaY6hIkhpjqEiSGmOoSJIaY6hIkhpjqEiSGmOoSJIaY6hIkhpjqEiSGmOoSJIaY6hIkhpjqEiSGmOoSJIaY6hIkhpjqEiSGmOoSJIaM6vfBdRFxMrAIuC9wLrAb4CTgBMzczQiRoADgPcDawCXAXtl5i/7VLIkqWbQeioHAUcA3wB2Ac4EjgH2reZ/EjgQOAp4K7AasDgiVpv+UiVJ7QampxIRKwEfBT6fmYdXkxdHxFxgn4g4GdgHODgzj6vWuYTSm9kD+EIfypYk1QxST2U14GvAWW3TE5gLbAfMAc59YkbmXcCPgZ2mqUZJ0jgGpqdSBcQHO8zaGVgCzKv+/eu2+TcBb+5haZKkCRqYUOkkIt4D7AB8CHgW8FBmPty22H3VPElSnw3S4a9lRMRC4BTg28AJwAgw2mHREeDxaSxNkjSGgQyViPgI8HXge8DCzBwF7gFmR8TT2hafU82TJPXZwIVKRBxBOZPr68DutcNdN1J6Jeu3rbIBZTBfktRnAxUqEbEI2B84FnhXZj5am3058BdgQW351YFtgcXTWackqbOBGaiPiLWAzwHXAt8Cto6I+iI/AY4HPh0RjwM3AJ8A7gW+Mr3VSpI6GZhQAXYEZgObAld0mD+XcouWxykXQc6h9F7emZmOqUjSABiYUMnM04DTJrDoftWPJGnADNSYiiRpuBkqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMbM6ncBY4mIXYDTM3PV2rStgKs6LH50Zu4zbcVJkjoayFCJiFcC3wBG2mbNBx4Admib/vvpqEuSNL6BCpWImA0sAg6jhMfT2xaZD1yXmVdOd22SpOUbtDGVNwD7A/sCx3eYPx+4ZlorkiRN2KCFylXA+pl5HDDaYf6mwAsi4uqIeDgifhUR75zeEiVJY+n68FdEvBBYH1gDeAy4Hbg1M5d0u83M/N0477d29V4vpvRm7gLeBpwWEaOZ+bVu31eS1IxJhUpEvAz4IOUw1bxqcmswfbRa5lfAecC/ZOZ1DdUJcDewE3BNZt5WTfthFTafAgwVSeqzCYVKFSbHANsBfwLOp4xt3ATcSzmM9lxK0LyC0oP4UERcCByQmT+daqGZ+WfgBx1mnQ/sFBFzMvP+qb6PJKl7yw2ViDgK2BP4d2B74OLMfGw564xQTvt9B3BpRJycmXtPpdCI2IgSav+amQ/VZq0CPEg5W0waOiMj7WfOT87oaKfhx+mtoak6NPwm0lNZDdg4M38LEBEB5HgrZOYocCFwYUR8EjhoqoUC6wAnU8Zuzq5qGQHeAlxSvackqY+WGyqZ+d62SddGxPHAIZl57wTWvxl4d5f11V0MXAqcEhGrA7cB76ecZvzqBrYvSZqibk4p/hKwF3BjRLyn4XrGVB1yezOll3IocBYwF3h9Zv5kuuqQJI1tpJvjoBHxUsrA/fbAT4FFmXlZw7V1LSLWA25evHgx8+bNW97iWkENyjiCYyoaFEuWLGH77beHcr3gLd1so6uLHzPz55n5emA3YHXg4oj4ZkT4DS5JK7ApXVGfmWcDG1MG4ncGfhkRB1b38JIkrWCmfJuWzHw4M48AAvgv4BDg+ojYLyJ2jIjnTfU9JEnDoavbtETEi4DNgJfUfgKYUy2yHnAET15lvzQz15lqsZKkwTbpUImIYym3ammN7P0BuAE4k3L9SuvnEWBzYIvqVZI0w3XTU9kDuIgyjvKLzLx7nGVvopz6K0laAXQTKp8HLsjMK5ouRpqMJk6DnQn8/6BBMulQycxDelGIJGn4Lffsr+q5KVMSEetPdRuSpME3kVOKfxIRx1bPLZmUiNggIk6lPNFRkjTDTeTw1+aUW7LcHBGXUO69dV51o8hlVHcNfhmwDfBW4FXAd/HsL0laIUzkLsW/A/4uIl4DfIQSMMdFxAPALSz7kK61gWdQrk85F9gmMy/vTemSpEEz4YH6zLyYco+veZTHCb8K2IASJo8Dv6Xcmv4/gR9k5h+bL1eSNMi6OftrCfDl6keSpCd0dZuWlohYA3gh8Bhwc2be00hVkqSh1O29v7YBPgtszZO3a3ksIhYD+2bmdQ3VJ0kaIt3c++u1wA+AB4ATgRuBlYGNgIXAZRHxKoNFklY83fRUPk056+tV7YPxEXEocCXwGcrzVaSOZsqtRWZKO5owCE+wVP918zyVzYGTO53dlZm3AycBr5lqYZKk4dNNqNwOPH+c+X9FuXZFkrSC6SZUDgcWRcRTDm9FxNbAh4FDp1qYJGn4dDOm8rfAHcA5EfFL4BfAw8CGwMuBh4C3RcTbauuMZub2Uy1WkjTYugmVHSi3YbmVckuWrWrzbq1evSuxJK2Aurmi3sCQJHXUzZiKJEkdGSqSpMYYKpKkxhgqkqTGTOkuxb0UEbsAp2fmqrVpI8ABwPuBNYDLgL0y85f9qVKSVDeQPZWIeCXwDZ68A3LLJ4EDgaMojyteDVgcEatNb4WSpE4GqqcSEbOBRcBhlLsgP702b1VgH+DgzDyumnYJ8BtgD+AL016wJGkZg9ZTeQOwP7AvcHzbvFcAc4BzWxMy8y7gx8BO01WgJGlsgxYqVwHrVz2R9vtgb1S9/rpt+k21eZKkPhqow1+Z+btxZj8LeCgzH26bfl81T5LUZ4PWUxnPCE/tvbSmPz7NtUiSOhimULkHmB0RT2ubPqeaJ0nqs2EKlRspvZL2G1puAOT0lyNJajdMoXI58BdgQWtCRKwObAss7ldRkqQnDdRA/Xgy8/6IOB74dEQ8DtwAfILy6OKv9LU4SRIwRKFSOYAyKL8PZSzlcuCdmemYiiQNgIENlcw8GDi4bdqjwH7VjyRpwAzTmIokacAZKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxszqdwGafiMjI1Naf3R0tKFKJM009lQkSY0xVCRJjTFUJEmNMVQkSY0xVCRJjTFUJEmNMVQkSY0ZuutUIuK5wB87zPo/mbn7dNcjSXrS0IUKsFn1uiNwb236nX2oRZJUM4yhMh+4PTMv6HchkqRlDeOYynzgmn4XIUl6qmHtqfwlIi4HtqSMrxwHfD4zvSmVJPXRUIVKRKwEbAI8AOwD3Aq8EfgM8FfAof2rTpI0VKECjABvAm7NzF9V034UEXOAj0fEkZn5l/6VJ0krtqEKlcx8DPjPDrPOB/YEXgRcN61FSZKeMFShEhFrU3oqZ2fmH2qzVqleO12/IkmaJsN29tds4EvA29um7wbckJlLp78kSVLLUPVUMvPmiPh34LCIeBy4Hvg7Sqgs6GtxkqThCpXKHsBBwIeBtSjBsltmntvXqiRJwxcqmfkgcED1I0kaIMM2piJJGmBD11NR/42MjPS7BM1ATfxejY56U41+s6ciSWqMoSJJaoyhIklqjKEiSWqMoSJJaoyhIklqjKEiSWqMoSJJaoyhIklqjKEiSWqMoSJJaoyhIklqjKEiSWqMoSJJaoyhIklqjKEiSWqMoSJJaoyhIklqjKEiSWqMoSJJasysfhcgSU0ZGRmZ8jZGR0cbqGTFZU9FktQYQ0WS1BhDRZLUGENFktQYQ0WS1JihPPsrIt4LfAyYB1wNfDQzr+hvVZKkoeupRMQ7gFOAbwC7AXcDP4iI9ftamCRpuEIlIkaAQ4FTM/OQzDwP2AX4I/CRvhYnSRquUAFeBLwQOLc1ITMfAb4P7NSvoiRJxbCNqWxUvf6qbfpNwIYRsXJmPgasDLB06dLprG1ozJo1bLtdmj5Llizpdwl9U/vOXLnbbQzbt8uzqtf72qbfR+l1PRO4F1gLYOHChdNX2RDZYIMN+l2CNLC23377fpcwCNYCft3NisMWKq0b+7TfnKc1/fHq9SpgG+A24LFpqEuSZoKVKYFyVbcbGLZQuad6XRW4vTZ9DiVQHgDIzIeAS6e3NEmaEbrqobQM20D9jdVr+/GbDYDMTG8vKkl9NGw9lRuB3wILgAsAIuJpwP+inAHWN5O9IDMiXgYcC2wN/Ak4ETiyHowRcR3w0rZV78zMNWrLLAAOo5wZdwPwicz8XiONWrbeRtsXEeP9AfCuzPy36hTyeyg907r/zsytum7MU2vt6mLaiFgVuA7YOzO/3TZvG+AoYFPgd8BnMvOrbcsM5L6rrTde+94EHARsAtxJOSPzE5l5XzV/WvZd9V69aN/QfvZq63VsX68/e0PVU6m+kD4L7BkRh0fEG4HvAGsAX+xXXZO9IDMingf8kDI29PfAqcDhwN61ZZ5OOdttP+Bvaz871pbZDvg2cBGwK3ANcHZEvGLQ29fWptbPpcDvgf9bLbM+5Zf6nW3Lvatfbauttyrld2/dDvM2Bs4HbgbeAnwX+JeI2L22zEDuu9p647XvdZQQ+Xm1zcOBtwJn1Bbr+b6raulF+4b2s1dbb8z20ePP3rD1VMjMkyJiFWAR5YLHq4EdM/OmftTTfkFmNe1CIKv6PtRhtQ9Q/t/vkpl/Bs6LiNnA/hFxbHXtzSbA04DvZOYvx3j7TwEXZuZe1b/Pj4gXAgdQLgod2PZl5pVt77MAeDWwXWbeUU2eTxkr+3a1nUZ12TYiYlvKB/35Y2x6P+AW4G3VH0LnR8Rc4JOULyIY3H03kfbtA1yWme+urXM3cGZEbJKZv6DH+656z161b5g/e8ttX68/e0PVU2nJzKMzc93MfEZmvrLP9/3q5oLMHYDFbTvrHOA5wMurf88H/sKT40jLqIL1lfX3rXwH2CEiuj7PvE2v2veEKnC+CHwrM39UmzUf+HWvvpTo/mLac4Brx1lmB+B7bWN85wCbRsTaA77vWrWO174rKYcz67J6bf0F3et9B71r3zB/9mD57XtCLz57QxkqA2a5F2SOsU6n5evbm085Vn1GRNwbEfdExFeqbi2UkxNmjbGdVYAXTK4ZY+pV++r2BNYBPt42fT7wUERcEBF/jog/RMSR1ThaE7ppG8A2mfn3wB3tMyLimcDaY2yz9Z6DvO9gnPYBZOZhmfmttsk7V6+tv+x7ve+gR+1juD97sPz21TX+2TNUpm4iF2R2WqfT8vXtzQfWBH5GORHhQMox1XMm8L71+VPVq/YBEBErUbrxZ2TmrW3rzAc2pPylthNwDLAX8KVJ1D+ebtpGZl7X5TZb8wd53y2vfU8REZsB+wNnZWbrdNRe7zvoXfuG+bM34f3Xq8/e0I2pDKCJXpDZPm+sMzBay38cmF07/nlJRNwBfKs6s6h1Uedk3rcbvWpfyw6Uv/z+ocOy/wTcl5nXVP++OCIeBT4TEYdk5m/GrXz5umlbE9vsxft2W8uURMR8ypmYvwPeV5vV630HvWvfMH/2JqMnnz1DZeomdEFmh3XaT9VbtTaPzPyfDuudX71uBrSOf7ZvZ05bXVPVk/bVLKAcu/1J+0Yy87IO2z6fcgbgpsBUv5i6advy3FvbZl19v9wzgWWa0Iv2PSEiXkv56/12YIfMvLM1bxr2HfSofUP+2ZuMnnz2DJWpq1+QWT/2Od4FmTfS+QJOgIyIWcDbgZ+1/YKvUr3+kXJc9fExtnM/5fTAJjTevrbpO/HkGVFPiIjVgN2Bi2qHVGDZ/wdT1U3bxpWZ90fEbYzd/hsowTOo+25CImIX4EzgesrZl3fU5k3HvoMetG8GfPYmoyefPcdUpq5+QSawzAWZi8dYZzHlLJH6MdEFlMHBqzPzUeAQ4OC29XYDHgGuyMwHgcvr71t5M+WXoal7njXevtp21qCcLXQlT/Uw5Qyj9tMmdwPuopzdMlXdtG0iFgM7tw2kLgCuy8zbB3zfLVdE/A0lUK4Ctq0HSmU69h30oH0z4LM3Ib387NlTmaLMHI2IzwInRMRdwGXAB6ldkBkRGwJza8doT6IMep0XEZ+ndKn3B/bLzIerZQ4HvhQRx1Iunns55TqH42rHMz8DfD8iTgXOBv435QKl1wxB+wBe1nqbDu/7YER8AfhYRNxJ+RC/HvgosCgzp9r177ZtE3EU5Qv3PyLiy5Rj12+nXAjaMqj7biK+TPmCPQLYJCLq827IzD/1et9BT9s3zJ+9ierZZ8+eSgMy8yRgX+AfKd3JZ7PsBZkHAVfUlr+N8kUzq1r+fZRbPBxVW+ZUymDZ6yi/2O+j3BLiY7Vlzqve87WUX+z5wIKmr9vpRfsqz6te7x7jrQ+iXEi4EPge5S/Bf87ME6baplqtk2rbBLf5M8opthtQ9svOwD9l5n/UlhnIfbc8EbFeVesc4Lxq3frPdrXt9nTfQc/239B+9iahZ5+9kdFR78EoSWqGPRVJUmMMFUlSYwwVSVJjDBVJUmMMFUlSYwwVSVJjDBVJUmO8ol7SE6oHUF0LrJGZz+53PRo+9lQk1R0KLOl3ERpehooERMRFEXFLv+vop4jYEngj5b5WUlc8/CUVhzPGk/S6ERHPotxXaQT4f5m5dYdlVqPciXYu5Xb4z27dzjwi1gRuA76bmbu0rXci5c618xq4/Xlrm7MoN4r8QBPb04rLnooEZOaFmXnO8pecsC0pgfIg8NKIGOmwzME8+cjY/2kLiC1a0+srVNt5M/CdpgKlsjdwTWZe1OA2tQIyVKTe2LJ6PZvSA1rmgU4R8RJKr6AVZP/dtv7m1Wv7UwhfDqxTW2/KqtunfwDYp6ltasXl4S/NGBFxMOU24V/NzL06zD8e2BOYn5nXt827CFgvM9er/n0L5fGpl1KeBbMh5YFJx2TmiRMo56+r169SnrWxKVB/it4XKc91/xHlGeE/bVt/rFDZlfKY2dYjbYmIz1Fuyx6U523sRrlF+uXAuzPztxHxj8CHgY2rdny81jPbhnII7ufVs1GeDjwrIpYCb8nMyyfQXgmwp6KZ5QLgv4APRsRG9RkR8VJKoJzYHijjeANwHOU5Fh+hPBP8hIh44wTW3RK4FbiY8iS91kORiIg3UR7lujfw0mpyp57Kn2oPhWrZFfh+Zj5Sm7YF5TDbWcBzgU8D36Q80+bEiDiB0gs5g/KwqecCp0fE3Gr9MyhPAdy8+nkPcF/13095frk0HnsqmjEy8/KIOIjSu9iM8jz4lmMoA+cHT2KTLwA2z8xrACLibMrzxxdSHlDVUUTMATYCzs3MRyLiekpPpfU42KOBH2XmWRHxUcpzzW+orf9M4EXUeiPV9I0pvZED295yC8rzw4/IzG/Wlt+K8oCw7wNbtYIoIh6p/n/MBxZXj8d9sLben4DRzFw6kf9JUp2hopnmF9XrS1oTImJXyl/te2bmWE+66yRbgVL9Y2lE3A6suZz1NqccBWgd0roa+JvqvxdRDqXtHhErUcLv6sx8vLb+ZtX67Ye+FgAPUQ7Ltdo2j/J42XPrgVK5C3iUcgis3rO5t3p9tFPx1WC9Fz6qKx7+0oySmXcBf6AKlYiYTXlm/NWUU2Yn4w8dpj0ErLyc9VrjKfVQeXFErEt5TOuXMvNaSm9mDmOPp7RP3xW4MDPvr01rnRBwRoc6XkbpEd3RNn3j6vUpzyeXpsqeimaiG3iyp7I35cyrbdt6AxMx2eVbWl/09VCZRfnif5QyrlFfbrlnflU9kq2A97Yt2zr1+Mr6xIh4AWXw/Uqeakvg9x7eUi/YU9FMlEBExDrAAcCZmXnxNL7/XwNLM/O26t9XV6+vAD6VmXfWloPOPZUHWHZMaAEwCpzbtuwWwN2ZeVPb9PZga1+n03RpygwVzURJuTbk65Tf8X2n642rGzK+hNqXdjWO83HgE8DJtcW3pAyQX19bf2XKYauftfWsdgUuzcz2Q3JjBUQrsJbpBUXEesBzxlhHmjIPf2kmao0VvA44ODNvncb33pwy5rLMl3ZmHjnGsj/LzMdq015COZOrfuhrdeA1lGtRqE1/DrAucGaHbW8J/DEzf9thOu31SU2xp6KZqBUqtwKdvsx7aUJf2tVV7M/usFxrPKV+fcjOlD8A26+ib42njNVTGevQ13Lrk7o1Mjra5O2DpP6rvrB/BRyYmYf3u57JiIhTKYPx62fmLdW0s6t/bz7eutIg8PCXZqKxTskdaNUV7v8A/LgVKJUrmPzp0FJfGCqaica6b9ZAiojXUU4X/mfgaZRbwjxhjPEYaSAZKpqJtgBuH6LrME4Bng9cBuyemUMRhlInjqlIkhrj2V+SpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTGGCqSpMYYKpKkxhgqkqTG/H9p/Q5PmyfQ4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(gamma_list*4/1000, bins=12,color=\"k\",normed=True)\n",
    "plt.ylabel(r\"p($\\gamma$)\")\n",
    "plt.xlabel(r\"$\\gamma$ in $MJ/m^{4}$\")\n",
    "plt.xlim((0.03,0.18))\n",
    "#plt.savefig(\"membrane_gamma300.png\",format=\"png\",dpi=300,bbox_inches='tight',facecolor=\"white\",backgroundcolor=\"white\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.01072736905533967 4.98264007400951 26.027704926927523\n"
     ]
    }
   ],
   "source": [
    "print(gamma_std_list.mean(),gamma_list.std(),gamma_list.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.04290947622135868 19.93056029603804 104.11081970771009\n"
     ]
    }
   ],
   "source": [
    "print(4*gamma_std_list.mean(),4*gamma_list.std(),4*gamma_list.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "464.479225829314\n"
     ]
    }
   ],
   "source": [
    "print(gamma_list.std()/gamma_std_list.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "400.0\n"
     ]
    }
   ],
   "source": [
    "print(0.024/6e-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "simulated  0\n",
      "simulated  1\n",
      "simulated  2\n",
      "simulated  3\n",
      "simulated  4\n",
      "simulated  5\n",
      "simulated  6\n",
      "simulated  7\n",
      "simulated  8\n",
      "simulated  9\n",
      "simulated  10\n",
      "simulated  11\n",
      "simulated  12\n",
      "simulated  13\n",
      "simulated  14\n",
      "simulated  15\n",
      "simulated  16\n",
      "simulated  17\n",
      "simulated  18\n",
      "simulated  19\n",
      "simulated  20\n",
      "simulated  21\n",
      "simulated  22\n",
      "simulated  23\n",
      "simulated  24\n",
      "simulated  25\n",
      "simulated  26\n",
      "simulated  27\n",
      "simulated  28\n",
      "simulated  29\n",
      "simulated  30\n",
      "simulated  31\n",
      "simulated  32\n",
      "simulated  33\n",
      "simulated  34\n",
      "simulated  35\n",
      "simulated  36\n",
      "simulated  37\n",
      "simulated  38\n",
      "simulated  39\n",
      "simulated  40\n",
      "simulated  41\n",
      "simulated  42\n",
      "simulated  43\n",
      "simulated  44\n",
      "simulated  45\n",
      "simulated  46\n",
      "simulated  47\n",
      "simulated  48\n",
      "simulated  49\n",
      "simulated  50\n",
      "simulated  51\n",
      "simulated  52\n",
      "simulated  53\n",
      "simulated  54\n",
      "simulated  55\n",
      "simulated  56\n",
      "simulated  57\n",
      "simulated  58\n",
      "simulated  59\n",
      "simulated  60\n",
      "simulated  61\n",
      "simulated  62\n",
      "simulated  63\n",
      "simulated  64\n",
      "simulated  65\n",
      "simulated  66\n",
      "simulated  67\n",
      "simulated  68\n",
      "simulated  69\n",
      "simulated  70\n",
      "simulated  71\n",
      "simulated  72\n",
      "simulated  73\n",
      "simulated  74\n",
      "simulated  75\n",
      "simulated  76\n",
      "simulated  77\n",
      "simulated  78\n",
      "simulated  79\n",
      "simulated  80\n",
      "simulated  81\n",
      "simulated  82\n",
      "simulated  83\n",
      "simulated  84\n",
      "simulated  85\n",
      "simulated  86\n",
      "simulated  87\n",
      "simulated  88\n",
      "simulated  89\n",
      "simulated  90\n",
      "simulated  91\n",
      "simulated  92\n",
      "simulated  93\n",
      "simulated  94\n",
      "simulated  95\n",
      "simulated  96\n",
      "simulated  97\n",
      "simulated  98\n",
      "simulated  99\n"
     ]
    }
   ],
   "source": [
    "acf_all=None\n",
    "N = 100000\n",
    "for i in range(100):\n",
    "    print(\"simulated \",i)\n",
    "    ts=simulate_membrane(N)\n",
    "    # calculate autocorrelation function\n",
    "    f = np.fft.rfft(ts)\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[:20000]\n",
    "    t = np.arange(autocorr.size)*delta_t\n",
    "    \n",
    "    if acf_all is None:\n",
    "        acf_all = autocorr\n",
    "    else:\n",
    "        acf_all = np.vstack((acf_all, autocorr))\n",
    "        \n",
    "acf_std = np.std(acf_all, axis=0)\n",
    "acf_mean = np.mean(acf_all, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ErrorbarContainer object of 3 artists>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAZn0lEQVR4nO3db6xddb3n8feZUtEpJVxTFEkF2pJ+g9oSHhACpjEOMCU2mMYalCAEHoAPQBLggm0UQnEIjZaQGuxA/TP3ARoJTeSPxTZ4TOYyXFC8I1zqhK8FWgJOgcjfCmEscObBWls227PPXmefvc/e+6z3KzlZ7W/9fmv91jrrrM9ef/fYxMQEkqT6+U+D7oAkaTAMAEmqKQNAkmrKAJCkmjIAJKmmDABJqqlDqlaMiIuBa4DFwGPAlZn58BT1TwNuBE4C3gJ+DVydmS821dkNfLql6cuZuajyEkiSujJW5TmAiLgA+B/ADcCjwDeAzwInZubeSeqfAPxv4AHgvwP/BHwHOACcnJkHI+JDwF+Ba4H/2dT8YGb+e8v0DgVOBvYD705zGSWpruYBnwAezcz/1zqy4xFARIxR7Pi3ZebGsuwBIIErgMsnaXYZxc56XWYeLNvsAX4HnAncD3wKmA/ck5lPdujGycCDnfoqSZrUKuB/tRZWOQV0PHAscG+joPwEvwM4q02bPwL/p7HzbzQrh0vK4UrgbWBPhT7sB/jpT3/KUUcdVaG6JOmFF17gvPPOg3If2qpKACwvh0+1lD8DLIuIeZn5gdMymbl1kumcXQ4bn/ZXAi8Dd0bEfwUmgLuAKzLzQEvbdwGOOuooFi9eXKHLkqQmk546r3IX0OHlsHWnfKBsv6DTBCLik8Bm4PfAb8rilcBRwOPAGuDbwDrg7gp9kiTNUJUjgLFy2Hq1uFH+3lSNy53/OEVYfDUzG9P5JnBoZj5S/v/BiHgJ+HlErMpMz/lLUh9VCYDXy+FC4MWm8sModv5vtmsYEZ8BfkVxsffMzHy6MS4z/zBJk53l8ES86CtJfVUlABoXaZfywesAS4Fs+kT/ARFxCsXO/w3gv2TmnqZxhwBfAx5vCYKPlMO/VOu+JKlbVa4B7AGeA9Y2CiJiPsV5+/HJGkTEcRQ7/xeB05p3/gCZ+Q6wEbi+pek64CDQ9gEzSVJvdDwCyMyJiNgE3BoRrwIPUdznvwi4BSAilgFHNp3P30Jx8fhS4JiIOKZpks9m5n6Kp4Rvj4gtwH0U9/pfB3w/M5/tydJJktqq9C6g8rbOq4Hzge3AEcDqzHymrHIt5af28ujgCxRPoP2sLG/+Oa+c5jbgIuDzFAFwCcXTwtf0YLkkSR1UehXEoJWnlPaOj49P+zmA49bvAGDfpjW975gkDbHnn3+e008/HWBJZu5rHe/bQCWppgwASaopA0CSasoAkKSaMgAkqaZqEwArrt816C5I0lCpTQAcePudQXdBkoZKbQJAkvRBBoAk1ZQBIEk1ZQBIUk0ZAJJUUwaAJNWUASBJNWUASFJNGQCSVFMGgCTVlAEgSTVVqwBYtmHHoLsgSUOjVgHw7vB//bEkzZpaBYAk6X0GgCTVlAEgSTVlAEhSTR1StWJEXAxcAywGHgOuzMyHp6h/GnAjcBLwFvBr4OrMfLGpzipgM7AC+DNwU2b+pIvlkCRNU6UjgIi4ALgNuANYB7wG7IqIJW3qnwCMAweAc4F/Bj5btpnfVGcnsBf4EnAf8OOI+PJMFkiSVE3HI4CIGANuALZl5say7AEggSuAyydpdhmwH1iXmQfLNnuA3wFnAvcD64F9wLmZOQHsjIgjgeuA7TNbLElSJ1WOAI4HjgXubRSUO/UdwFlt2vwRuLmx8280K4eNo4YzgF+WO/+Gu4EVEXF0hX5JkmagyjWA5eXwqZbyZ4BlETEvM99tHpGZWyeZztnl8MmIWAAc3WaajXn+3wp9kyR1qcoRwOHl8EBL+YGy/YJOE4iIT1Jc7P098JsO02yepySpT6oEwFg5bH2RQqP8vakalzv/8XJeXy1P+cxomjPh+4AkqVAlAF4vhwtbyg+j2FG/2a5hRHwG+DeKT/RnZubT5ag3pphm8zx7zvcBSVKhSgDsKYdLW8qXAtlyEffvIuIU4F+Bd4FVmfkfjXGZ+VeKu4QmmybAnyr0S5I0A1UD4DlgbaOgvJd/DcWpnX8QEccBvwJeBE7LzD2TVBsHzo6IeU1la4HdzQ+LSZL6o+NdQJk5ERGbgFsj4lXgIYr7/BcBtwBExDLgyMx8pGy2heK0z6XAMRFxTNMkn83M/RQXhR8F7oqIH1LcFvo14JyeLJkkaUqVngQub+u8Gjif4iGtI4DVmdm4bfNa4GH4+9HBF4B5wM/K8uaf88ppPk5xa+hS4Bflvy/KzLt6sWCSpKlVfhdQZt4M3Nxm3IXAheW/DwLzK05zF7Crah8kSb3j20AlqaYMAEmqKQNAkmrKAJCkmqplABy33tdBSFItA0CSZABIUm0ZAJJUU3M+AOaNda4jSXU05wPg6ZvWDLoLkjSU5nwASJImZwBIUk0ZAJJUUwaAJNVUbQPAp4El1V1tA0CS6s4AkKSaMgAkqaYMAEmqKQNAkmrKAJCkmqpFACz88CGD7oIkDZ1aBMAT168edBckaejUIgDa+crtDw+6C5I0MLUOgN/ufWXQXZCkgal8cjwiLgauARYDjwFXZmbHj9ARsRDYDVyVmdtbxu0GPt3S5OXMXFS1X5Kk7lQ6AoiIC4DbgDuAdcBrwK6IWNKh3ULgHuCYScZ9CFgOrAdObfrxhL0kzYKORwARMQbcAGzLzI1l2QNAAlcAl7dp9zmK0Ph4m0l/CpgP3JOZT06/65KkmahyBHA8cCxwb6MgMw8CO4Czpmh3N/DEFHVWAm8Deyr1VJLUU1UCYHk5fKql/BlgWUTMa9NuVWaeA7zUZvxK4GXgzoh4IyJej4gflaeNJEl9ViUADi+HB1rKD5TtF0zWKDN3d5juSuAo4HFgDfBtiusLd1fokyRphqrcBTRWDifalL/X5by/CRyamY+U/38wIl4Cfh4RqzLzwS6nK0mqoMoRwOvlsPXUzGEUO/83u5lxZv6haeffsLMcntjNNLvhN4NJqqsqAdC4SLu0pXwpkJnZemTQUUQcEhEXRsRJLaM+Ug7/Mt1pdrJv05peT1KSRlrVAHgOWNsoiIj5FOftx7uZaWa+A2wErm8ZtQ44CPiOBknqs47XADJzIiI2AbdGxKvAQ8BlwCLgFoCIWAYcOckpnancCNweEVuA+4CTgeuA72fms9NbDEnSdFV6EjgztwJXA+cD24EjgNWZ+UxZ5Vqm+ak9M7cBFwGfpwiAS4DvULxuQpLUZ5XfBZSZNwM3txl3IXBhm3H7eP+OodZx/wL8S9U+SJJ6p9ZvA5WkOjMAJKmmDAB8FkBSPRkAklRTBoAk1ZQBIEk1VasA8HUQkvS+WgWAJOl9BoAk1ZQBUFpx/a5Bd0GSZpUBUDrw9juD7oIkzSoDQJJqygCQpJoyACSppgwASaqp2gXAwg9X/goESZrTahcAT1y/uu043woqqU5qFwCSpIIBIEk1ZQBIUk0ZAJJUUwaAJNWUAdDCO4Ek1YUBIEk1VfmpqIi4GLgGWAw8BlyZmQ9XaLcQ2A1clZnbW8atAjYDK4A/Azdl5k+qd787pyz5KL/d+0q/ZyNJQ63SEUBEXADcBtwBrANeA3ZFxJIO7RYC9wDHTDLuBGAnsBf4EnAf8OOI+PJ0FqAbd3791H7PQpKGXscjgIgYA24AtmXmxrLsASCBK4DL27T7HEVofLzNpNcD+4BzM3MC2BkRRwLXAdvbtJEk9UiVI4DjgWOBexsFmXkQ2AGcNUW7u4EnpqhzBvDLcuff3GZFRBxdoV9944VgSXVQJQCWl8OnWsqfAZZFxLw27VZl5jnAS60jImIBcHSbaTbPU5LUJ1UC4PByeKCl/EDZfsFkjTJzd5fTbB4vSeqTKgEwVg4n2pS/18V8+zFNSdI0VAmA18vhwpbywyh21G92Md83pphm8zwlSX1SJQD2lMOlLeVLgWy5iFtJZv4V2N9mmgB/mu40p6vTF8N4IVjSXFc1AJ4D1jYKImI+sAYYn8G8x4GzWy4irwV2Z+aLM5huJVN9MYwk1UHH5wAycyIiNgG3RsSrwEPAZcAi4BaAiFgGHJmZj0xj3puBR4G7IuKHFLeFfg04Z3qLIEnqRqUngTNzK3A1cD7FQ1pHAKszs3Hb5rVAx9dCtEzzceBsitM+vyj/fVFm3jWd6UiSulP5XUCZeTNwc5txFwIXthm3j/fv7mkdtwvYVbUPkqTe8W2gU/BCsKS5zACQpJqqdQDMm/TElCTVQ60D4Omb1gy6C5I0MLUOgCq8DiBprjIAJKmmDABJqikDQJJqygCowOsAkuai2gfAKUs+OuguSNJA1D4A7vz6qYPugiQNRO0DQJLqygCoyOsAkuYaA0CSasoAkKSaMgDo/P3ADcs2eBpI0txhAFD9+4HfnehzRyRpFhkAklRTBsA0feX2aX31sSQNLQNgmn6795VBd0GSesIAkKSaMgBK+zZV/3YwHwqTNBcYAJJUUwaAJNVUtSeggIi4GLgGWAw8BlyZmW1viYmIzwBbgFOAV4AfAN/NzImmOruBT7c0fTkzF1VeggE5bv2OaZ02kqRhU+kIICIuAG4D7gDWAa8BuyJiSZv6HwN+DUwA5wDbgBuBq5rqfAhYDqwHTm36qfZUliRpRjoeAUTEGHADsC0zN5ZlDwAJXAFcPkmzS8tpfzEz3wLuj4hDgQ0RsSUzDwKfAuYD92Tmkz1Zmhnat2nNtC7wehQgaZRVOQI4HjgWuLdRUO7AdwBntWlzBjBe7vwb7gY+Cpxc/n8l8DawZ5p9liT1QJUAWF4On2opfwZYFhHz2rSZrH7z9FYCLwN3RsQbEfF6RPwoIhZW6NPQ8JZQSaOqSgAcXg4PtJQfKNsvaNNmsvrN01sJHAU8DqwBvk1xfeHuCn2SJM1QlbuAxsph67swG+XvtWnT7t2ZjfrfBA7NzEfK/z8YES8BP4+IVZn5YIW+DQWvBUgaRVWOAF4vh62nZg6j2Jm/2aZNa/2FTePIzD807fwbdpbDEyv0qy/ckUuqiyoB0LhIu7SlfCmQzff1t7SZrD5ARsQhEXFhRJzUUucj5fAvFfo1VLwWIGnUVA2A54C1jYKImE9x3n68TZtx4IyIaL4+sJbiou9jmfkOsBG4vqXdOuAgMJLvXDYEJI2SjtcAMnMiIjYBt0bEq8BDwGXAIuAWgIhYBhzZdEpnK/ANivv/v0dxSmcDsD4z/1bWuRG4PSK2APdR3B56HfD9zHy2VwsoSZpcpSeBM3MrcDVwPrAdOAJYnZmNWzuvpelTe2bup3gW4JCy/iXAtzJzc1OdbcBFwOcpAuAS4DsUr5sYqJlcB/AoQNKoGJuYGP4vuo2I44C94+PjLF68eFbmOdMduReTJQ3a888/z+mnnw6wJDP3tY73baB94pGApGFnAPSR3x8saZgZAG304hSO3x8saZgZAH3mqSBJw8oAmAWGgKRhZABMoZd38hgCkoaNATCLDAFJw8QAmGWGgKRhYQB0MG+sc53pMgQkDQMDoIOnb+rPE73Hrd/Bsg0GgaTBMQAG6N0JjwYkDY4BUEG/3+tz3PodBoGkWWcADBGDQNJsqvKdwJplU4WAbxmV1CsGQEX7Nq0Zik/nk/XBUBhtzb9Tf5ewbMMO3q34lnrX18wYAHNA1WDq5o/Fo5GZW3H9Lg68/U6lunVe39NZTw2t62uurKPZWi4DoEZ6fQQzyC/N6Wbes71z6Pf6nqs7u35Ma97Y1Ld0z8aR9TCcQWjlN4JN0zD+Euea2Tjd1s+d51zcRjrtQKdjlNdP1e2m18vY7fba6RvBPALQ0JmNHURjHr0KglHeqVXhMyuFubYODIBpGpaLweqNKhdg/X1rrjIApJI7etWND4J1Ya5cfJNUbwZAlxZ+2IMnSaPNAOjSE9evHnQXJGlGKn+MjYiLgWuAxcBjwJWZ+fAU9T8DbAFOAV4BfgB8NzMnmuqsAjYDK4A/Azdl5k+6WI6B8IKwpFFW6QggIi4AbgPuANYBrwG7ImJJm/ofA34NTADnANuAG4GrmuqcAOwE9gJfAu4DfhwRX+52YSRJ1XUMgIgYA24AtmXmxsy8H/gi8BfgijbNLqU4uvhiZt6fmf8NuAnYEBHzyzrrgX3AuZm5MzOvpAiY62ayQLPNC8KSRlWVI4DjgWOBexsFmXkQ2AGc1abNGcB4Zr7VVHY38FHg5KY6v2w+JVTWWRERR1fr/nAwBCSNoioBsLwcPtVS/gywLCLmtWkzWX2A5RGxADh6qjoV+jVUDAFJo6ZKABxeDg+0lB8o2y9o02ay+o1xU02zeZ6SpD6pEgBj5bD1rXGN8vfatGn3lrn3upzm0PMoQFI/vPW36b0mu6oqAfB6OVzYUn4YxY76zTZtWusvbBr3xhTTbJ7nyDEEJPXaf/5Qfx48rRIAe8rh0pbypUC2XMRtbjNZfco2fwX2T1HnTxX6NbQMAUmjoGoAPAesbRSUt3KuAcbbtBkHzigv9jasBV6meIisUefslovIa4Hdmflite4PL0NA0rDreFyRmRMRsQm4NSJeBR4CLgMWAbcARMQy4MjMfKRsthX4BnB/RHwPOBHYAKzPzL+VdTYDjwJ3RcQPKW4L/RrFg2Nzgk8KSxpmlZ4EzsytwNXA+cB24AhgdWY2btu8Fni4qf5+ih36IWX9S4BvZebmpjqPA2dTnPb5RfnvizLzrhku01DZt2mNRwOShpJfCTmLPBqQ1I1+fSWkbwOdRR4NSBomBsAAGAKShoEBMCCGgKRB82utBqg5BLw+IGm2GQBDokoYVDlqmE6QVD0KMZykuckAGEIzOT3Uj1NLnaZpQAyXqX5fzb+rRj1/f/VlAGjGqu5w1Bu9/oDQqw8N3fyu+xFCvfwQNNe3XwNAfTVqnzJ7sfOosqxz8SaAYTty7YW5/uHGANCs6NUf+LINO3i3w7OLg96ZDHr+mh1z4SYOA0Aj5emb3Llq+OzbtIav3P4wv937Sk+mNVsMAEnqgTu/fuqguzBtPggmSTVlAEhSTRkAklRTBoAk1ZQBIEk1ZQBIUk0ZAJJUUwaAJNXUqDwINg/ghRdeGHQ/JGlkNO0z5002flQC4BMA55133qD7IUmj6BPA062FoxIAjwKrgP3AuwPuiySNinkUO/9HJxs5NjHR4dWKkqQ5yYvAklRTo3IK6O8i4mLgGmAx8BhwZWY+PEX9zwBbgFOAV4AfAN/NzImmOquAzcAK4M/ATZn5k74txCzq0/raDXy6penLmbmox92fddNdX03tFgK7gasyc3vLOLevf2w31fpy+3q//mnAjcBJwFvAr4GrM/PFpjpdb18jdQQQERcAtwF3AOuA14BdEbGkTf2PUaywCeAcYBvFyryqqc4JwE5gL/Al4D7gxxHx5f4tyezo0/r6ELAcWA+c2vSzum8LMkumu76a2i0E7gGOmWSc29c/tptqfbl9vV//BGAcOACcC/wz8NmyzfymOl1vXyNzBBARY8ANwLbM3FiWPQAkcAVw+STNLqVYxi9m5lvA/RFxKLAhIrZk5kGKDW0fcG75KXdnRBwJXAdsn2SaI6GP6+tTwHzgnsx8chYWZVZ0ub6IiM9R/FF/vM2k3b4+2K7T+nL7et9lFDe+rCv/9oiIPcDvgDOB+5nh9jVKRwDHA8cC9zYKypWyAzirTZszgPFyZ9ZwN/BR4OSmOr9sPsVR1lkREUf3qO+D0K/1tRJ4G9jT6w4PWDfrC4r188QUddy+PqjT+nL7et8fgZsbO/9Gs3LYOGqY0fY1SgGwvBw+1VL+DLAsIiZ70GF5m/oAyyNiAXD0VHW67Osw6Pn6KocrgZeBOyPijYh4PSJ+VB7Wj7Ju1hfAqsw8B3ipdYTb16Tarq+S21cpM7dm5g9ais8uh0/2YvsapQA4vBweaCk/QLEcC9q0max+Y9xU02ye5yjqx/qC4g/0KOBxYA3wbYrzmXfPsL+D1s36IjN3dznN5vGjqB/rC9y+2oqIT1Jc7P098JsO02yeZ1sjcw0AGCuHrQ8uNMrfa9Om3YMO73U5zVHRj/UF8E3g0Mx8pPz/gxHxEvDziFiVmQ922+EB68e24PY1fW5fkyh3/uMUYfHVzJworyt0PU0YrQB4vRwuBF5sKj+MYkHfbNOm9dBxYdO4N1rKmqfZPM9R1I/1RWb+YZJ2O8vhicCo/oF2s746cfuaJrevf1Temv0riovjZ2Zm45UOM96+RukUUOOi0NKW8qVAtlwEaW4zWX3KNn+luMrers6fuuzrMOj5+oqIQyLiwog4qaXOR8rhX7ru7eB1s76m5PY1PW5f/ygiTgH+leIVOKsy8z8a43qxfY1aADwHrG0UlPfCrqE4NJrMOHBGebGkYS3FRabHmuqc3XIRZi2wu/lhixHU8/WVme8AG4HrW9qtAw4CHR8AGmLdrK8q3L4qcvv6oIg4juKT/4vAaZk52Z1RM9q+RuYUUHnOaxNwa0S8CjxEcZ/sIuAWgIhYBhzZdP5wK/ANivvZv0dxCLkBWJ+ZfyvrbKZ4UdJdEfFDituqvkbxINTI6uP6uhG4PSK2UDx0cjLFPcffz8xnZ2fpeq/L9VWF29f01pfb1/vrawvFhdxLgWMiovnBuWczcz8z3L5G6QiAzNwKXA2cT/GQwxHA6sxs3PZ0LU2fEsoVdAZF0G0HLgG+lZmbm+o8TnFr1VLgF+W/L8rMu/q+QH3Wp/W1DbgI+DzFH+glwHcoHm8fadNdXxWn6fY1vWm6ffH3o4MvULzN82dlefPPeeU0Z7R9+TZQSaqpkToCkCT1jgEgSTVlAEhSTRkAklRTBoAk1ZQBIEk1ZQBIUk0ZAJJUUwaAJNXU/wec7zSTxpujLQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(t,acf_mean,yerr=acf_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [10:34:30<00:00, 38.07s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Variables]]\n",
      "    kappa:  10 (fixed)\n",
      "    sigma:  125 (fixed)\n",
      "    ga:     28.9101353 +/- 0.13163437 (0.46%) (init = 25)\n",
      "    eta:    0.25 (fixed)\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'dg_list2' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-27-24fd3db35144>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     36\u001b[0m     \u001b[0mg_list\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvaluesdict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'ga'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m     \u001b[0mdg_list\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcovar\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 38\u001b[0;31m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdg_list2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     39\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m     \u001b[0mgS_list\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mres\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatchain\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'ga'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'dg_list2' is not defined"
     ]
    }
   ],
   "source": [
    "g_list=[]\n",
    "dg_list=[]\n",
    "\n",
    "gS_list=[]\n",
    "dgS_list=[]\n",
    "\n",
    "def residual(pars, x, data=None, eps=None):\n",
    "    # unpack parameters: extract .value attribute for each parameter\n",
    "    parvals = pars.valuesdict()\n",
    "    ga = parvals['ga']\n",
    "    kappa = parvals['kappa']\n",
    "    sigma = parvals['sigma']\n",
    "    eta = parvals['eta']\n",
    "\n",
    "    model = fitcorr(x,kappa,sigma,ga,eta)\n",
    "\n",
    "    if data is None:\n",
    "        return model\n",
    "    if eps is None:\n",
    "        return model - data\n",
    "    return (model-data) / eps\n",
    "\n",
    "for y in acf_all:\n",
    "    pars = lm.Parameters()\n",
    "    pars.add('kappa', value=kappa, vary=False)\n",
    "    pars.add('sigma', value=sigma, vary=False)\n",
    "    pars.add('ga', value=ga, vary=True)\n",
    "    pars.add('eta', value=eta, vary=False)\n",
    " \n",
    "    result = lm.minimize(residual, pars,kws = {'x':t,'data':y,'eps':acf_std})\n",
    "    #result.params.add('__lnsigma', value=np.log(0.1), min=np.log(0.001), max=np.log(2))\n",
    "    res = lm.minimize(residual, kws = {'x':t,'data':y,'eps':acf_std},method='emcee', nan_policy='omit', burn=300, steps=1000, thin=20,\n",
    "                     params=result.params)\n",
    "\n",
    "    print(lm.fit_report(result.params))\n",
    "    g_list.append(result.params.valuesdict()['ga'])\n",
    "    dg_list.append(np.sqrt(result.covar[0,0]))\n",
    "    print(g_list[-1],dg_list[-1])\n",
    "\n",
    "    gS_list.append(res.flatchain['ga'].mean())\n",
    "    dgS_list.append(res.flatchain['ga'].mean())\n",
    "    print(gS_list[-1],dgS_list[-1])\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "gS_list.append(res.flatchain['ga'].mean())\n",
    "dgS_list.append(res.flatchain['ga'].mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28.91013532789121 0.13163437031611613\n"
     ]
    }
   ],
   "source": [
    "print(g_list[-1],dg_list[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28.916058902756763 28.916058902756763\n"
     ]
    }
   ],
   "source": [
    "print(gS_list[-1],dgS_list[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}