membrane fluctuations gamma.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 0x7fc68867eba8>]"
      ]
     },
     "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,4) for m in range(0,4)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15\n",
      "[1.         1.         1.41421356 2.         2.         2.23606798\n",
      " 2.23606798 2.82842712 3.         3.         3.16227766 3.16227766\n",
      " 3.60555128 3.60555128 4.24264069]\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 0x7fc679a34940>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAECCAYAAADkaECYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUVf7H8ffNpDdCCCGE0BLkgNJBQFAsYFtFUJG1997L6k/Xbe6ua2PXxbWgrq513RUVK6IQFBVFOojl0JHQaxoEUub3x524IUySiWQyyeTzep55Bu79npsTRvPJvfeccx2v14uIiEh1EaHugIiINE0KCBER8UsBISIifikgRETEr8hQd6ChGGNigCOBTUB5iLsjItJceID2wDxr7b6qO8ImIHDD4fNQd0JEpJk6Bvii6oZwCohNAK+++ioZGRmh7ouISLOwefNmLrjgAvD9DK0qnAKiHCAjI4OsrKxQ90VEpLk56NK8blKLiIhfCggREfFLASEiIn4pIERExC8FhIiI+KWAEBERv8JpmOvPtmTm6+yxM6mIiMYbEUVFRNT/3j1ReD0xeD3R4ImmwhMDnhi8kTEQGQtRcTiRcRAdhxMVR0RUPJHR0URFOERFRhAZ4RDliSA6MuKn92jfe4zvzxERTqj/CUREDqKAAEq+/4i+W98jkjKinUNfpWOfN5I9xFJMLHu9MRQTQ4E3lmLiKCCOQm88Rb73QuLZ4yRQHJnM3shWlEQmsy8qBaISiI32EB8dSXy0h7hoD/HRHhKiI4mPjiQxNpKkGPc9MSaSpFj31SoumlZxUURH6uRQRA5NwAFhjLkKuAvIAhYDt1trv6qlvhcwERgC7ASeAB621nqr1BwDTAB6AxuAB6y1z1c7Tn9fzVHAbuBN4G5rbXGgfa/LkBue+99fvF6oKIPy/XjL9lNeWkJ56T4qSvdSUbrP9ypxX/v34t2/l4rSPXhL9+It3Qule2D/HthfTFTpHqJLi0kp3UNEaTGe0iIiSzcTVVpIVFkREVQc2JEy36sESomiMCKJ3U4K20lhu7cVWyqS2VyezKryVmzyprKJNmz2tqbMz8eYEO0hJd4Ni5T4KFIToklLjKFtUgxpidG0SYghLSmGdN8r0qNAEZEDBRQQxpiLgUnAH4F5wE3AR8aYvtbaNX7q04EZwDJgPDAAuB93pt4EX01PYBrwHvB74CTgOWNMgbX2DV9NN2AW7hpLY4Bs4EEgEbjs533LdXAc8ESBJwonOoFIgnSa5fXC/mLYVwAl+bBnJ+zd5XvtJGrvLlKLt5NavI3soq1QtAKKtwL73aW1Kg+DQ2l8OnvjMiiKyWBXTCabIzuwwWnPGm8G6/fHsGtvKcs25LOjaD+F+8oO6oonwiE9KYb2rWJp3yqO9q1iyUyJo0taPJ3bJNCxdbzOSERaoDp/9hljHNxgeMZae59v23TAArcBN/tpdoPv2GdYa/cAU32rrd5jjJlorS0F7gbWAuf5ziqmGWPaAr8D3vAd5/fAGmCMtbbM97UjgZuNMVG+4zRPjgMxie4rOTOwNl4vlOyGwi1QsAHy83AKNhCdv4Hogjxa5a+iw+aZ9Kqo8s8SFQ+p2dDVQLsj2N/mcHYkdmMLaWwv2s+WwhI255ewcXcJm/L38t2mAmZ8v4V9Zf87u4lwoEPrOLq0SaBLmwS6ZyTRIyMJk5FEcmxUA//DiEhTEcgvx92AzsC7lRustaXGmA+AU2poMwrI9YVDpbeB3+Cuuvqlr+aVqpecfDUXGmMygc24Zw2/qQwH39d+AvdyVcvjOBDX2n2l9/BfU14GBXmwYxXsXO2+dqyCvHmw7E2icdf1bR/TCtodAe37QtYgGDwYWvUBx8Hr9bKjeD/rduxh7fZi1u0oZu2OPazbUczbizYccBbSISWOHhlJ9GifRJ+sFAZ0ak3bpJhG+ecQkeAKJCC6+95XVtu+GsgxxnistdXv7HYHPvVTD9DdGLMEyKzhmJXtY4EkYLMx5mXgTNxLVC8Dd1Rft1x8PJHQuov7YuSB+0ryYev3sOVb32sZLHwRvn7K3Z+YAR2PxMkaTFqnoaRlDWBg59YHHMLr9bIxvwS7uYDvNxViNxfyw+YCZi3fRlmFm/WdUuMZ0CmFAZ1bM6BTa3pkJOkeh0gzFEhAJPveC6ttL8SdR5EAFPhp46++cl9tx6ysaev780TgQ9yzid649zI8wHUB9F2qim0FnYa6r0rlZW5Q5M2D9XMhby58/567L6YVZI+AnJHQbSSkdMJxHDqkxNEhJY4TerT76TAlpeV8uzGfhet2s2DdLr5ctYO3F28EIDk2kmMOa8ux3dsyontbMlrFNuZ3LSI/UyABUTlI31vD9mpDcX7aV72+UkWAx6y8uP29tfZy359zffcgHjLG3Get3VxX56UOnkjI7Oe+Bl/lbivaCutmw8pcWDXzf4HR5jDoNgoOPwM6DoWI/50VxEZ5GNg5lYGdU7kK90xjw+69LFi3i9krtzNr+TY++MZdbt60S+JY05YTeqQzuEuq5oGINFGBBES+7z0J2FJleyLuD3J/w03zffVVJVXZV1BtW9VjVtYU+f48rVrNdOAR4Ajc+xTS0BLT4Ygz3ZfXC9uX+8IiFxb8y70klZTp7u91NnQY4N4fqcJxHLJax5PVOp4x/Trg9XqxWwr5bPk2Zi3fxguz1/LMZ6vJSI7l9D7tGdOvA706JOM4CguRpiKQgFjhe8/mwHsG2YCtdpO5apvsatsq/26ttUXGmE211CwH9uCeYURXq6k8s6jpDEUakuNAW+O+jroe9hWCnQbfvgXznoU5T0BKZzcs+l8Ead1qOIxDj4xkemQkc/WIHIr3lTHzh628s3gjL361ln9+sYbstATO6JfJmH4d6JqW0Ljfp4gcJJA7hyuA9cDYyg3GmCjgNCC3hja5wChjTNX/y8cCO3An2VXWjDbGeKrVLLPWbrHWFgJzgLONMVX7eRqwD1gQQN+locUkQZ9z4LzX4FcrYMyTkHYYfPkPeHwgvHwm2A+hovYZ6QkxkYzum8k/LxnEvHtH8cBZvUlPjmFi7gqOn/Aplzw/l1nLt+H16vcAkVCp8wzCWus1xjwIPG6M2QXMBm4E0oBHAYwxOUBba+0cX7MncSfTTTXGPAL0Be7BnQG931czAXfS3WRjzLO4w14vxJ1YV+le4GPgdWPM00B/4NfA36y1+UhoxaVA/wvcV9FWWPAizH8eXjsXUjrBoCtgwMUQn1rrYVLiozlvcCfOG9yJzfklvD5/PS/PWcclz8+lW3oilw/vypn9OxAX7an1OCLSsJxAf0MzxtwB3IIbDItxh5p+5dv3AnCJtdapUj8IdwTSQNx7F09aax+qdsyTgYeAHsCPwF+stS9UqzkR+DNuyGzDndH9gLW2olpdF2BNbm6unkkdSuWlYKfC3Gdh7efugob9zodjfgWtOgR8mH1l5XywdBPPfbGGbzcWkBIfxfmDO3HlMdmkJlS/6igiP1deXh4jR44E6GqtXVt1X8AB0dQpIJqgLd/B3Kdh0avgRMCgy+Ho2yCpXd1tfbxeL/PW7uL5L9bw8XebSYiO5Nrjcrh8eFedUYg0gNoCQrOXJHjaHQ6jJ8LNC6HPeJj7DEzsC9N/B8U7AjqE4zgM7prKpIsG8tGtIxiS3YZHPrIcP+FTXp+3nvKK8PgFR6QpUkBI8KV0gjGPw43z3DkUsx+DiX3g04egtCTgwxzWLol/XjKI1685ioxWsdz15lJOnfgZM3/YopvZIkGggJDG0yYHznoGrp8DOSfAp3+Bp4bB6ln1OszgrqlMuX4YT14wgNJyL5e/MJ8rX5zP1oLAw0ZE6qaAkMaX3gN++TJcNAW8FfDSGTDlWijeHvAhHMfhF73b8/FtI/jNaT35YuV2Tnz0M95ZvEFnEyINRAEhoZNzAlz/lTvC6Zs34PFBsOgVd/Z2gKI8EVx5TDYf3nIMOW0TuOU/i7nulYVsL9JajiKHSgEhoRUVByN/C9d+Dm17wDs3wAunQ35evQ6T3TaRydcO455TezDzh62c9OhnTPWt/SQiP48CQpqG9J5w6VQY/RhsWgyTjoEV0+t1CE+EwzXH5vDBzUeT1TqO619dyB2vL6Gk9NCfMy7SEikgpOmIiICBl8DVs9yn7L06Dmb8wV2SvB4Oa5fEW9cN4+aRh/HmwjzOe3YOWwt1A1ukvhQQ0vSkdYMrZ8DAS+GLR+HF0VCwsV6HiPREcPuJ3Zl04QB+2FTImMdns2yDVmcRqQ8FhDRNUXHuJLuznoVNS2DS0bByRr0Pc0qv9rxx3VE4wDmTvuJD3ZcQCZgCQpq2PuPh6k8hsR28Mg6+qv/jyI/IbMXbNw6nR/skrnt1IY/lrtBQWJEAKCCk6WvbHa7MhZ6nw0e/ho9/CxX+HmRYs/SkWF67aihn9e/A36Yv59b/Lqa0vH7HEGlpAnlgkEjoRcfDOS/Ch3fBl4+5y4uPeRw8UXW39YmN8vDX8X3JSU/kkY8s5RVe/v7LfkR69HuSiD8KCGk+IjzwiwmQlAEz/wzF22D8SxCTWHdbH8dxuOH4bkR5HP4y9Qc8EQ5/G98Pj56LLXIQBYQ0L44DI+5070m8dwu8eDqcPxkS29brMFePyKGswsvD0ywex+GRc/oqJESqUUBI8zTgYkhoC5Mvg+dPctd1at2lXoe4/rhuVFR4mfDxciIiHB4+uw8RCgmRn+jiqzRf5lS45F3YsxNeGgMF9R/CeuMJh3HrqMN4Y0Ee97z1DRV6voTITxQQ0rx1HAwXvuWuBPvyWDcs6unWUd25+YRu/Hf+eu59e5mGwIr4KCCk+csaCOe9BjvXwCtnQUlBvQ9x24nduf64HF6b+yNPzVoVhE6KND8KCAkPXUfA+Bdh8zfw2nlQurdezR3H4c6TDaP7ZvLIR5aZP2wJUkdFmg8FhIQPcyqc+TSsmw2TL4Xy0no1dxz3RvURmcnc8tpiVm4tCk4/RZoJBYSEl97j4LS/wvJpMOUaqKjfUt9x0R6evmgQMVERXP3SfPL31i9kRMKJAkLCz5FXwKg/wLI3Ydo99W7eISWOpy4cyPpde7j5tUWUa2STtFAKCAlPR98GR90Ic5+GhS/Vu/mRXVL545hezFq+jYc/+iEIHRRp+hQQEr5G3QfZx8MHd8D6efVuft7gTlw0tDNPz1rN24s2BKGDIk2bAkLClycSxj0PSe3h9YugsP4jk343+nCGdE3l/95cyrcb9cAhaVkUEBLe4lPh3H9DSb4bEmX769U8yhPBkxcMoFVcFHe8voR9ZXq+tbQcCggJfxm9YMwTsP5rd7nwemqTGMNDZ/fhh82FTJyxIggdFGmaFBDSMvQ6C4bfCgv+BfP/Ve/mx/dI59wjOzJp1ioW/rgrCB0UaXoUENJyjPwd5IyEqXfCj1/Xu/m9p/Wkfas4fvX6Evbu16UmCX8KCGk5Ijww7jlolQWTL6n3wn5JsVE8ck4fVm8v1tBXaREUENKyxLV212wq3g7v3wb1XLl1WE4alw7rwr9mr+XLVduD1EmRpkEBIS1P+75w/D3w3dvwzeR6N/+/U3rQNS2BOycvpWhfWRA6KNI0KCCkZRp+K3QcAh/8CvLz6tU0LtrDhHP6sil/L/d/8F2QOigSegoIaZkiPHDmJKgog7evg4qKejUf2Lk11xybw2tz1/OJ3RqkToqElgJCWq7UbDjlAVjzGXw9qd7Nbx11GN3SE7nv3W81gU7CkgJCWrYBF0P3U2HGH2Dr9/VqGhPp4benH87aHXt4YfbaoHRPJJQUENKyOQ6c8RjEJMFbV9d7KY5ju7dlZI90/jFzJdsK9wWpkyKhoYAQSUyH0RNh81KY9WC9m997Wk/2lZUz4SMbhM6JhI4CQgSg5+nQ/0L44lHYsLBeTbPbJnLpsC68vmA9yzZoxVcJHwoIkUon/wXi09znR9TzUaU3jTyM1Pho7nvvW7z1nHwn0lQpIEQqxbZyQ2LjQljwQr2aJsdG8auTDfPW7uKDbzYFp38ijUwBIVJV73HQ5RjIvQ+KttWr6fhBHenZPpkHpv5ASamGvUrzp4AQqcpx4LS/wv49MP139WrqiXD4/ejD2bB7L898tjpIHRRpPAoIkeraGhh2Eyz5N6ydXa+mQ7Pb8IveGTz16So25e8NUgdFGocCQsSfEXdCq07uDevy0no1vefUnpR7vTwyTcNepXlTQIj4Ex0Ppz4E276HOU/Wq2nH1HguG9aFtxdvYNW2oiB1UCT4FBAiNenxCzC/gE8frPeKr1eNyCYm0sPjM1cGqXMiwaeAEKnNKQ+6DxWadne9mqUlxnDh0E68s3gDq3UWIc2UAkKkNq07w7F3wvfvwcoZ9Wp69YgcoiMjePwTnUVI86SAEKnLUTdB6y4w/ff1mmHdNimGC4d05p3FG1m7vTh4/RMJkshAC40xVwF3AVnAYuB2a+1XtdT3AiYCQ4CdwBPAw9Zab5WaY4AJQG9gA/CAtfb5Go7nALlAhLX2uED7LXLIIqPhhN/Cm1fA0teh33kBN7362GxenrOOxz9ZyYRz+gaxkyINL6AzCGPMxcAk4BXgbGA38JExpmsN9enADMALjAeeAe4H7qhS0xOYBqwBzgLeA54zxoyroRtXAccH0l+RBnfEWdC+H8z8M5SWBNwsPSmWC4Z0ZsqiDazbobMIaV7qDAjfb+5/BJ6x1t5nrZ0KnAFsB26rodkNuGcnZ1hrp1pr/ww8ANxjjIny1dwNrAXOs9ZOs9bejhtAB01fNcZkAQ8DG+vzzYk0mIgIOOlPUJAHc5+uV9Nrj80mMsLRiCZpdgI5g+gGdAberdxgrS0FPgBOqaHNKCDXWrunyra3gVTgyCo171e95OSr6W2Myax2vKeAKcC8APorEhxdR0C3E+Hzv8KenQE3S0+O5bzBnXhr0QZ+3LGn7gYiTUQgAdHd917915/VQI4xxlNDG3/1AN2NMQlAZm01lRuMMRcAg6lyeUokZEb9AUoK4Iu/1avZdcfl4IlweEIjmqQZCSQgkn3vhdW2F/raJ9TQxl995b7ajvnT1zTGtMW90X2ztTbwX9lEgiWjF/Q9D75+Bnb/GHCzdsmxnD+4E28uzGP9Tp1FSPMQSEA4vvfqT0Gp3F5RQ5uanppSUY9j/gP4ylr73wD6KdI4TrjXfZ95f72aXXtsDhGOw5Of6ixCmodAAqLyGYpJ1bYn4v4g9zc0I99PfVKVfQW1HBMg3xgzBjgNuMkYE2mMicQNEMf3dweRUGiVBUOvhaX/hc3fBNwso1Us5w7uyBsL8thSEPhIKJFQCSQgVvjes6ttzwZstZvMVdv4q8fXpgjYVEvNcuBM3MBYA5T6XmcAI3x/PjaAvosEx9G3uU+gm/77ejW74uiulFV4eWXOuiB1TKThBBoQ64GxlRt8Q1VPw5245k8uMMp3M7rSWGAH7iS7yprR1W5yjwWWWWu3AH/AHfFU9TULWOj784IA+i4SHHGt3SXBV+XCqk8Cbta5TQKjerbj1a9/1FPnpMmrcya1tdZrjHkQeNwYswuYDdwIpAGPAhhjcoC21to5vmZPAjcBU40xjwB9gXuAu621+301E3CHrU42xjyLO+z1QtyJdVhr1+LOk/iJMWa3b9/8n/n9ijScwVfB15Pgk/sh+zj3aXQBuHx4V6Z/t4V3Fm/gl0d2CmoXRQ5FQDOprbVPAncCFwFvACnAydbaymGpvwW+qlK/CfcHfqSv/mrgXmvthCo1S4DRuJeVpvj+fJm1dvIhfk8ijSMyxr3UlDcPVgd+FjE0O5We7ZN57os1eL01jeUQCT0nXP4DNcZ0Adbk5uaSlZUV6u5IS1G2Dx7rD606wuXTAj6LmDx/PXe+sZRXrhjC0YelBbmTIjXLy8tj5MiRAF19V25+otVcRQ5F5VnE+jmwZlbAzc7ol0laYjTPz14TxM6JHBoFhMih6n8RJLWHWQ8H3CQm0sOFQzsz84eteqCQNFkKCJFDFRXrnkWsmw1rPg+42QVDOhPtieCFL9cGr28ih0ABIdIQBlwCiRkw66GAm7RNiuGMfplMnp9H/p7SIHZO5OdRQIg0hKhYGH4LrP0c1s4OuNllw7uwt7Sc/8wLfF0nkcaigBBpKIMug4T0ep1FHJHZiqHZqbz45VrKyv0tayYSOgoIkYYSFeeeRayZBT/Oqbve54qjs9mYX8JH324JYudE6k8BIdKQBl0G8Wnw6YMBNzmhRzqd28RryKs0OQoIkYYUnQDDb3ZnVq+fG1ATT4TDpcO6sGDdLr7Jy6+7gUgjUUCINLQjr4T4NvWaF3H2wCxioyL491zdrJamQwEh0tCiE2DodbByOmxeFlCT5NgoTu+TybuLN1C0ryzIHRQJjAJCJBgGXQFRCfDlPwJuct7gThTvL+e9JRuD2DGRwCkgRIIhPhUGXAzL3oD8vICaDOiUgmmXxH90mUmaCAWESLAcdT14vTDnqYDKHcfhvMEdWZKXz7INulktoaeAEAmWlE7Q62xY8ALs3RVQkzP7ZxETGaGZ1dIkKCBEgmn4zbC/COY/H1B5q/goTuvTnrcXbWTPft2sltBSQIgEU0ZvyBkJcyZBaUlATc4f3ImifWW8v2RTkDsnUjsFhEiwDb8FirfC0v8EVD6wc2u6pSdqToSEnAJCJNi6joD2fd0hrxV1L8jn3qzuxOL1u/l+U0EjdFDEPwWESLA5jnsWsWMl2KkBNTmrfweiIyM05FVCSgEh0hh6joGUzjD77+7Q1zq0TojmF70yeGvRBvbuL2+EDoocTAEh0hg8kTDsJsibF/BS4OcN7kRhSRkffKOb1RIaCgiRxtLvAohLhdkTAyof3DWV7LYJvKbLTBIiCgiRxhIdD4OvhuUfwvaVdZY7jsP5gzuxYN0u7ObCRuigyIEUECKN6cgrwBMNX08KqPysAVlEeRwmz18f5I6JHEwBIdKYEtOh9zmw+NWAlt9ITYhmZI92vL14A6V6ZrU0MgWESGMbeh2U7oGFLwVUPm5gFtuL9jPLbgtyx0QOpIAQaWwZvaHLMfD1M1Be93pLx5q2pCVG8+bCwJYNF2koCgiRUBh6PRTkwQ/v1Vka5YlgTL8OzPh+C7uK9zdC50RcCgiRUOh+MrTuCl89GVD5uIFZlJZ7eVdPm5NGpIAQCYUIDwy5FvLmQt78Ost7tk/miMxk3ligy0zSeBQQIqHS/wKISQ74iXPjBmbxzYZ8ftisBfykcSggREIlJsl9bvV3b0P+hjrLx/TrQJTH4U2dRUgjUUCIhNLgq8FbAfOerbM0NSGa4006UxZtpExzIqQRKCBEQql1Z+hxOsz/F+zfU2e5OydiH5+t0JwICT4FhEioDb0eSnYH9MS543uk0yYhWjerpVEoIERCrdNQaN/PvVldxxPnfpoT8d1WzYmQoFNAiISa47hnEduXw6qZdZaPG5jF/vIK3luqORESXAoIkabgiDMhsR3MfbrO0sMzkzm8veZESPApIESagshoGHQ5rPg4oGdFnD0wi6V5+XpOhASVAkKkqRh0OUREwdxn6iwd0y+TyAhHC/hJUCkgRJqKxHTodbb7rIiS2mdLpyXGcJxJZ8qiDZoTIUGjgBBpSoZcA/uL3JCow7iBHdhWuI/Zq3Y0QsekJVJAiDQlHQZAxyHw9dN1Dnk9vkc6reKitPSGBI0CQqSpGXIN7FoDK6fXWhYT6WF03/Z89O1mCktKG6lz0pIoIESamp5nQFJmQKu8nj0gi31lFXz4zeZG6Ji0NAoIkabGEwVHXgGrP4GtP9Ra2q9jCtlpCbyh0UwSBAoIkaZo4KXgialzyKvjOJw9MIu5a3ayfmfdi/2J1IcCQqQpSkiD3ufAktdg7+5aS8f27wDAlEV1P1NCpD4UECJN1ZBroHQPLHq51rIOKXEcld2Gtxbm4fV6G6lz0hIoIESaqvZ9oPNw9zJTRXmtpWcPzGLtjj0s/HFXI3VOWoLIQAuNMVcBdwFZwGLgdmvtV7XU9wImAkOAncATwMPWWm+VmmOACUBvYAPwgLX2+WrHuQS4HegGbAReBf5irdVaxxL+hlwDr18M9kPoeXqNZaf0yuC3by/jzYUbGNg5tRE7KOEsoDMIY8zFwCTgFeBsYDfwkTGmaw316cAMwAuMB54B7gfuqFLTE5gGrAHOAt4DnjPGjKtScynwL+BjYAzwFPAr4O/1+B5Fmi9zGrTqVOeQ18SYSE7plcH7SzZSUlr72YZIoOoMCGOMA/wReMZae5+1dipwBrAduK2GZjfgnp2cYa2daq39M/AAcI8xJspXczewFjjPWjvNWns7bgD9rspx7gRetdbeaa2dYa39G3AfcK0xJrG+36xIs+OJhCFXw7ovYNOSWkvPHpBFQUkZud9vbaTOSbgL5AyiG9AZeLdyg7W2FPgAOKWGNqOAXGtt1XF3bwOpwJFVat6vesnJV9PbGJNpjIkAPgJeqnZsCzi+PomEv/4XQVQCzJlUa9lROW3ISI7lLc2JkAYSSEB0971XX6R+NZBjjPHU0MZfPUB3Y0wCkFlbjbW2wlp7u7W2+noDo4ES3LMPkfAXlwL9L4Blb0DhlhrLPBEOY/t34NPl29hWuK8ROyjhKpCASPa9V38ySaGvfUINbfzVV+6r7ZhVv+YBjDEnA5cBT1pri2vvtkgYGXItlJfC/OdqLTt7QAfKK7y8u0SPI5VDF0hAOL736gOsK7f7W3LS8VNfqeLnHNMYcwLwJjAHuLeW/oqEnzY50P0UmPcclJbUWHZYuyT6ZLXS40ilQQQSEPm+96Rq2xNxf5D7+00+3099UpV9BdW2VT1m1a8JgDHml8BUYClwurW25v9DRMLV0Otgz3b4ZnKtZecMzOL7TQUs25Bfa51IXQIJiBW+9+xq27MBW+0mc9U2/urxtSkCNtVSs7xygzHmWuA14DPgRGtt7esOiISrriOgXS93yGstM6bP6NuB6MgIJs9f34idk3AUaECsB8ZWbvANVT0NyK2hTS4wynczutJYYAfuJLvKmtHVbnKPBZZZa7f4vs5Y4EngLdwzB913kJbLcdyziK3fwprPaixrFR/FKUdk8PZizYmQQ1PnTGprrdcY8yDwuDFmFzAbuBFIAx4FMMbkAG2ttXN8zZ4EbgKmGhATFcQAAA/vSURBVGMeAfoC9wB3V5kBPQGYB0w2xjyLO+z1QtyJdRhjYnEn520CHgMGGGOqdu0bBYa0OL3GwfTfu2cR2cfWWDZ+UEfeXbKR6d9tYXTfzEbsoISTgGZSW2ufxJ20dhHwBpACnGytrRyW+lvgqyr1m3B/4Ef66q8G7rXWTqhSswR3yGo2MMX358ustZUXWIcC7XCHw87yHb/qq2f9v12RZi4q1n1WxPJpsGNVjWXDctrQISWO13WZSQ6BEy6rPxpjugBrcnNzycrKCnV3RIKncAv8vZf7zIhfPFJj2d9nLGdi7go+v+t4slrHN17/pFnJy8tj5MiRAF2ttWur7tNqriLNTVI791LToldrfVbEuIHuL0pvLtBzIuTnUUCINEdDr4XSYlj4Yo0lWa3jGZ6TxuQF66moCI8rBdK4FBAizVH7vtDlGPdmdVnNy2qcMyiLvF17mbN6RyN2TsKFAkKkuTr6VijcBEtfr7Hk5CMySI6N1M1q+VkUECLNVc5IyOgNsydChb8VbyA2ysOYfh34cNlm8veWNnIHpblTQIg0V44Dw2+FHSvATq2xbPygjuwrq9ACflJvCgiR5uzwsdC6C3zxaI3Lb/TqkEyPjCQtvSH1poAQac48kTDsJtgwH9bN9lviOA6/PLIjS/Py+X5Tgd8aEX8UECLNXb8LIKGtexZRg7H9OhDtiWDyfC0DLoFTQIg0d1Fx7gOFVs6Azd/4LWmdEM2Jh7djyqI8LeAnAVNAiISDI6+A6ER3RFMNzhvciV17Svlg6aZG7Jg0ZwoIkXAQ1xoGXQbL3oRda/2WDO/Whm7pibzw5VrCZQ02CS4FhEi4GHo9OB748nG/ux3H4dJhXfhmQz4Lf9zVyJ2T5kgBIRIukjOh77mw6GUo2ua35Mz+HUiKjeSFL9c1cuekOVJAiIST4be4azN9Pcnv7oSYSH45qCMffrOJLQV6tLvUTgEhEk7SDoOeo2HuM7Bnp9+Si4/qQrnXy6tzdBYhtVNAiISb4+6BfYU1jmjq1CaekT3S+ffcH9lXpiGvUjMFhEi4aXc49BkPXz8NhZv9llwyrAvbi/ZryKvUSgEhEo6OuxsqSuGzCX53H90tTUNepU4KCJFwlJoNAy6GBS/4nRfhOA6XDOvC0rx8Fq2v+bGl0rIpIETC1Yg7IcIDnz7kd/dZlUNeZ69t3H5Js6GAEAlXyZkw+CpY+h/Y+sNBuxNiIhk/qCNTNeRVaqCAEAlnw2+DqAT45H6/uy8+qrM75PXrHxu5Y9IcKCBEwllCGxh2I3z/LmxYeNDuzm0SOMGk8++v12nIqxxEASES7oZeD3GpMPPPfndfNrwr24v261kRchAFhEi4i02GY26HVbmw9ouDdg/v1oZBnVvzj5kr9KwIOYACQqQlOPJKSGoPuX886NnVjuPwq5MNWwr28fJXWn5D/kcBIdISRMW5S3Cs/xqW/Oeg3UOz23DMYWk8NWsVRfvKQtBBaYoUECItRf+LoOMQ+OjXULzjoN2/Osmws3g/z3+xJgSdk6ZIASHSUkREwOl/h30FMP23B+3u2zGFkw5vx7OfrWb3nv0h6KA0NQoIkZak3eEw7GZY/Cqs+eyg3XecZCjaX8akWatD0DlpahQQIi3NsXdB6y7w/m1QeuAMapORxJi+mbzw5Rq2Fmp2dUungBBpaaLi4PRHYcdK+OJvB+2+dVR3Ssu9PDFzZQg6J02JAkKkJco5AXqPh8//BtuWH7CrS1oC4wdl8e+5P5K3a0+IOihNgQJCpKU6+S8QnQDv3woVFQfsuumEw3Ach8dyV4Soc9IUKCBEWqrEtnDSn2DdbPemdRWZKXFcOKQzbyzIY+XWohB1UEJNASHSkvW7EDoNg49/AwUHPn70+uNzSIiO5O43l1JeoafOtUQKCJGWLCICRk+E8lJ4/aIDRjWlJcbwx7FHMH/dLp76VDesWyIFhEhL17Y7nDkJ8ubBB7cfsFbT2H4dGN03k7/PWMESPZq0xVFAiAgcfgYce7d7L2LOUz9tdhyHP4/tRXpSDLf+dzHFWqepRVFAiIjr2P+DHqfDx/fCqk9+2twqLoq/ju/H2h3F/PmD70PYQWlsCggRcUVEwJlPQ9seMPlS2LHqp11H5bTh6hHZvDb3Rz7+dnPo+iiNSgEhIv8Tkwjn/hscB/5zPpQU/LTrjhMNR2Qmc/db32gZjhZCASEiB0rtCue8CNtXwJRrfppEFx0ZwcRz+1G8r4w7Jy/F69XQ13CngBCRg2UfC6c8CHYqvHMDlLnLf3dLT+Le03oya/k2nv5MK76Gu8hQd0BEmqjBV8HenfDpA7D7R/jlyxCfykVDOzNn9Q4e/PAHSkrLuWWkuyyHhB+dQYiIf44Dx90NZz0LeXPhn6Ngxyp3jaZz+zNuYBZ/n7GC37/7LRWaaR2WFBAiUrs+4+Hid6FkN/xzJKz9gkhPBI+M68PVI7J56at13PyfRewvq6j7WNKsKCBEpG6dj4IrZ0BCW3hpLCx6Fcdx+PUvenLPqT14f+kmrnhxnibShRkFhIgEJjUbrvjYDYt3rof3boGCTVxzbA4Pj+vDl6t2cP6zc9hZrOdZhwsFhIgELq41XPgWDL0eFr0Cj/WDj3/D+J7xTLpwID9sLuS0xz7n1a/X6ZJTGHACHctsjLkKuAvIAhYDt1trv6qlvhcwERgC7ASeAB621nqr1BwDTAB6AxuAB6y1z1c7zljgT0A3YDlwr7X2fT9frwuwJjc3l6ysrIC+JxE5BDtXw6cPwdL/QnQiHHU9SzpewB8+zmPRj7vpkBLH9cfncM7AjkRH6nfRpiovL4+RI0cCdLXWrq26L6BPzRhzMTAJeAU4G9gNfGSM6VpDfTowA/AC44FngPuBO6rU9ASmAWuAs4D3gOeMMeOq1JwAvAF8CpwJLAWmGGOGBtJvEQmi1Gw462m4fg50OwFmPUTfN0bwlpnJlNGRtEuM5N4pyzh+wqf8++sfdUbRDNU5D8IY4wB/BJ6x1t7n2zYdsMBtwM1+mt3gO/YZ1to9wFRjTAxwjzFmorW2FLgbWAuc5zurmGaMaQv8DjcUAH4PTLfW3uT7+zRjTGfg18AZP+cbFpEGlt4Dxr8EGxfDJ3/B+fyv9MfLmzHJ7Og2hDd3H8Yzb6/m0emdGZLdhgGdWjOgc2sOb5+sM4smLpCJct2AzsC7lRustaXGmA+AU2poMwrI9YVDpbeB3wBHAl/6al6pesnJV3OhMSYT2AUM4+AAegf4kzHGY60tD6D/ItIYMvvBBa/Dnp2wZhbOqk9IW/0J1xRN55oY2O1tw48r0lj7fSpzvG14z0kjpk1H0jK7kpTcmsSkJBKTkmmVlEzr5CRSEmKIi/LgiXA0ES9EAgmI7r736o+UWg3k1PCDujvuZaHq9QDdjTFLgMwajlnZfpuvf/5q4oCOuGcgItKUxKfCEWe6L6/XvVex+hNS8haQkr+ew3fn4RQswFOx371YXcNziEq8UewlEi8O5UTgxcHrRFCBg7eWq+Ne/IdJTdvDwea4bvT7v48a/LiBBESy772w2vZC3HsYCUBBtX3JNdRX7qvtmJU1+wKoEZGmzHGgTY77OvJKwPdDx+uF4u1QkEdF/kZKigso3lPE3uIi9u0pYn9JEWX7iqkoL8dbUY63ogKvtxxvhRe8ZeB1b3CCF8fr5afLEC10AcGKtO51F/0MgQREZexW/5ev3O7vzpPjp75SRYDH/DlfV0SaA8eBxLaQ2JaIzP7EA/Gh7pMcJJA7RPm+96Rq2xNxf0gX19Cmen1SlX0F1bZVPWZlTW1ft2q/REQkCAIJiBW+9+xq27MBW+0mc9U2/urxtSkCNtVSsxz3XkNFDTVFwMYA+i4iIj9ToAGxHhhbucEYEwWcBuTW0CYXGGWMSaiybSywA3eSXWXNaGOMp1rNMmvtFmvtXtzRTmM50BjgU41gEhEJrjrvQVhrvcaYB4HHjTG7gNnAjUAa8CiAMSYHaGutneNr9iRwE+78h0eAvsA9wN3W2sqFWiYA84DJxphncYe9Xog7sa7SA8AHxphngCnA+cBRwIif/y2LiEggApqlYq19ErgTuAh3ElsKcLK1tnJY6m+Br6rUb8L9gR/pq78ad4mMCVVqlgCjcS8ZTfH9+TJr7eQqNVN9X/M4X00fYGxtS3yIiEjDCHgtpqZOazGJiNTfIa/FJCIiLU84PZPaA7B58+ZQ90NEpNmo8jPTU31fOAVEe4ALLrgg1P0QEWmO2gOrqm4Ip4CYBxyDO79CQ2BFRALjwQ2HedV3hM1NahERaVi6SS0iIn4pIERExK9wugchQWaM+SvQzlp7Yaj7Ig3LGHM7cAXu6snzgWustftqbyXNgTHmHOAPQAzwcuWTQQOhMwgJiDHmZODiUPdDGp4xZghwGTAY6I1709Lfo4SlmTHGZAB/BY4HegLHGWNOCrS9AkLqZIxJB/4E3B/qvkhQ7AJutNYW+1ZnXgp0CnGfpGGciLu46VZrbSnwMnBuoI11iUlqZYxxgOeBO4CcEHdHgsBauxx3iX2MMe2BW4BLQ9knaTAdOPDRCBuBgNciUkAIAMaYc4HHq23eAjwHLLbWfu5btVeaqZo+Y2vtEb79nYEPgaestTMau38SFB4OfCqnQz3miWkehNTKGDMXiMP9jyoV94l+r1lrbwhpx6RBGWP6Ae8BD1tr/xHq/kjDMMZcDJxgrb3U9/fLgWHW2isDaa+AkIAZYy4FRmkUU3gxxrQFvsG9D/FGqPsjDccYkwnMwR2AsAP4GHjMWjslkPa6xBRmjDFnAK9aa5Oqbb8KuAv3+uNi4HY9V6N5CsJnfBuQDPzGGPMb37YPrLX3NmC35RD83M/cWrvRGHMHMAOIBt4JNBxAZxBhxRgzDJgGRFhrE6tsvxj4F/BH3PVWbgKGA32ttWtC0Vf5efQZtzyh/Mx1BhEGjDExuCNP/gQU4/6mULnPwf0P6JnKCTLGmOmAxf3NUePdmwF9xi1PU/jMNQ8iPJyK+8zvO4HqNxi7AZ2Bdys3+MZDfwCc0lgdlEOmz7jlCflnroAID/NwHxf4GAcOaQPo7ntfWW37aiDHGHPQQ0KkSdJn3PKE/DPXJaYwYK3dUMvuZN97YbXthbi/ICQABcHolzQcfcYtT1P4zHUGEf4c33v130Aqt1c0Yl8kOPQZtzyN8pkrIMJfvu89qdr2RNz/iIobtzsSBPqMW55G+cwVEOFvhe89u9r2bMD6FmeT5k2fccvTKJ+5AiL8rQDWA2MrNxhjooDTgNxQdUoalD7jlqdRPnPdpA5z1lqvMeZB4HFjzC5gNnAjkAY8GtLOSYPQZ9zyNNZnrjOIFsBa+yTuWOqLgDeAFOBka+3qkHZMGow+45anMT5zLbUhIiJ+6QxCRET8UkCIiIhfCggREfFLASEiIn4pIERExC8FhIiI+KWAEBERvxQQIiLilwJCRET8+n+o3GyRZUKTNAAAAABJRU5ErkJggg==\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dd5gURdrAf7uKBAVFwQAKGM4Cz4TxRFHBnHNExQSeigkF05lQEdPhqXCo4CHmCN+dCRVOUcQA5kCJCdRDwIAICMgy3x89vfb0VPd09VRP9+zW73n22ZmeDm+neqveekNNLpfDYrFYLBY/tWkLYLFYLJZsYhWExWKxWJRYBWGxWCwWJVZBWCwWi0WJVRAWi8ViUbJy2gKYQgjRFNgemA3UpSyOxWKxVAsrAesBb0spl3p/yLyCEEKsBJwH9AE6ADOB4cAwKaXXR3d74NXKS2ixWCwNgu7Aa94FmVcQwBXAJcC1wBs4J3Eb0AK4ybPebIAHH3yQddddt9IyWiwWS1Xy/fff06tXL8i3oV4yrSCEELVAf+BmKeX1+cUThBBtgYsoVBB1AOuuuy7rr79+ZQW1WCyW6qfINJ/1SerVgTHAU77lEmgrhFi18iJZLBZL4yDTIwgp5c9AP8VPBwHfSikXVVgki8ViaTRkWkGoEEKcDuwJnJu2LBaLxdKQybqJqQAhRC9gBPAEcGfK4lgsFkuDpmoUhBDiAuB+4Gmgl8/F1WKxWCyGqQoFIYQYDPwdR0EcKaVclrJIFovF0uDJvIIQQpwHXAr8AzhZSrk8ZZEsFksKrFixgpqaGv7+97+nLUqjIdOT1EKI9YAbgQ+BR4AdhRDeVaZahWGxNA6WLXMMB5deein9+/dPWZrGQaYVBLAP0BTYApii+L0t8ENFJbJYLKlSU1OTtgiNhkwrCCnlaGB0ymJYLJYMYMsjV57Mz0FYLBYL/KEg7AiiclgFYbFYqoIVK1YAVkFUEqsgLBZLVWBNTJXHKgiLxVIV2BFE5bEKwmKxVAULFy4EYPHixSlL0niwCsJisVQFX331VdoiNDqsgrBYLFVBhw4d0hah0WEVhMXSCFiyZAlDhw6lrq6oaFjV0KxZMwBatmyZsiSNB6sgLJZGwKBBg+jfvz8PPPBA2qLExvViqmYlV21YBWGxNAJ+/PFHAH777beUJYnPzJkzATtJ7XLKKadQU1OT6D21CsJiaQS4LqIrrbRSypLE55tvvklbhEwxevRoAFq0aMEHH3yQyDGsgrBYGgGuWaa2tnpfeVfJWYp55513Etlv9T4tFoslMq6CqOYRhJ17CCapKHOrICyWRkBDUBCdOnVKW4TMYhWExWKJzYMPPgjA77//nrIk8XHdW9u0aZOyJNljzpw5iezXKgiLpRFRzdHI7ijI5mIqxioIi8VSNtU8SS2lBGDevHkpS5I9XDdm01Tv02KxWLSp5t73ww8/nOj+3377bV577bVEj5EUSQVAWgVhsTQixowZk7YIsVm+fHmi+99hhx3o3r17oseoNqyCsFgaEXPnzk1bhNhYN9fKs7LOykKIjYAewHqolUtOSnmtCcEsFot5dt1118Dfcrlcpk1Q7dq1S1uERkdkBSGEOBoYA6wSsloOsArCkinmzZtH27Zt0xYjE7Ro0UK5/JtvvqFDhw6MHj2a3r17V1iqaBxzzDGMHDkybTEaFTompkHATGAvYBNgQ8XfRqYFtFjK4fnnn2fttdfm+eefT1uU1GjdunX95yZNmijXmT59OpDcZKcJbE3qyqNjYtoAuFBKOSEpYSwW07z11lsAvPTSS+y7774pS5MO8+fPr/+8+uqrK9dx3V+zbOf/6KOP0hah0aEzgngPZ5RgsVQNX3zxBQC33nprypJkg91220253E3BkeWEeP3796/IcXK5HI899hhLliypyPHi8Oqrr1bkODoK4gLgVCFEHyGENehaqgJv7/ndd99NUZJssGjRIuXyahhBVIpJkyZxzDHHMGDAgLRFUTJnzpxQZwOT6JiYvgSmAyOAEUII1To5KaWWZ5TFkiTeHvE222zT6O3YQXMQ33//PUDVBoqZZNasWQB8/vnnALz55pv85S9/QUrJpptumqZoQHJR0yp0GvN7gJ2Bt4HPgGSjViwVY9GiRVx33XVcddVV9XV/GwqmFcLChQupq6sLtOVnnaDr8cYbb1RYkuxyxhlnADB58mQA7rvvPsCZx8qCgliwYEHFjqWjIPYE7pRSnpuUMJZ0uOmmmxgyZAht27atmJ23Upg2mbRp04alS5eSy+VYsGABTZs2pWnTpkaPkSRejyYvWc7RlMvl6Nu3b8WO55bwdM1xS5cuBWCVVcI8/CvHDz/8ULFj6TwVi4CPkxLEkh7ffvstUN31ioP49NNPje7PbSzA8Qjaeeedje7fNP70FEGBcFmuEzF//vxU4h9c8+S9994LwOWXX15xGVRUMmW7joK4H2eSOhtq1GIM9wX429/+lrIk5nEL3SfFtGnTEt1/ufiV/vz586mpqSlKfJflEURWJs6zkqakkgpCx8T0LnAsMEMI8Twwh+J5CJtqw5JZVl111bRFqDj+EcMrr7wCwM0338xxxx1Xv3zllbPrW6JyvZ09ezbrrbdeCtKkTyXdb3WeCm+IZZ+AdWyqDUtm2WWXXdIWoeL4J6Xd0aJfcWR5BOE167m0a9eu0XqkvfzyyxU7ls5ToUqtYVNtWKqGAw88MG0RKsYGG2zA8OHDAxtRv0LIsoJYtmxZ2iJkim233bZix4o8gpBSJmvMtVgsxvj22285++yz6dWrl/J3/wgiy5PUWZmDyApJ18XwopPN9aQo60kpE6tIIoQ4GHhQStkyqWNYkmP+/Pmcfvrp3HXXXay11loVP36WU1knRdAIIkxBLF++PFNzEroN4ooVK7j++us588wzadOmTdnH33DDDTNVy9stvVoJdJ6C0ThzDKq3zPsUJqIghBDdcOZBMvOWjxw5kj59+jB79mzWXXfdtMXJPMOGDePJJ59ECMH111+ftjjaeJPFPfXUU9rbz58/n7q6uooqxyAF4feE8ZqYsmbb11UQL7/8MldeeSXvvvturPvkJyvxDy5PP/10xY6lY3jsAfTM/3c/7wucADyF49FkPJG8EKKpEGIg8F8yFr09cOBA4I+Iy2okjcZg8ODBFT+mCdwUDECsHmXr1q2N9GhL4b2nn332mXIdf14qr4LIWsI+XQXx3XffATB27FgjHj9ZM79ttdVWFTuWzhzEKyE/PySEeBi4lEJvJxPsl9/vAGAt4ELD+4/Nzz//DMC4ceM44ogjUpYmHuPGjUtbhIoxaNAg+vbtG5iPqBTemAKdRqPSk6zPPPNM/ef//e9/kbbxnk/WbP66fv/efFLNmzcvuxP0ySeflLW9abI6gijFBJLxYnob2FBKeTuFpqzMUMnAFdN4s50mTdqmi7lz5/Kvf/0r9vZuzxT0IrTbtGlT0dxNd955Z/3nqLm1vAqi2kcQWevxm6Zjx44VO5ZJBbE/Ts0Io0gpv5NSVq4Vi0E1u+FVstGeNGlSxY4VhMqnPireEcTdd98debtff/21osFN48ePr/8ctbGsRhPT5ptvrlwex2XXdEqWJNlggw0qdqxyvZhWwTH7HICT6fVm/3pJejVlhUq6nfmZNm0a7733Hqeddlqs7SupINziPWnin4z997//zf777x9qdlq0aBFjx46tylxVU6ZMibTeCy+8UP+5WhTExx+rU8PFGUH88ssv2ttUkrlz53LcccfxyCOPFFyPLl26JKrcTHkxuQz0fc+RkFdTlvjPf/6T2rG32247gKpQECYankmTJjF58mQuvfTSSOv36NGD//73v/XfvQVznn76aQ499FD22GMPXnrppcB9nHXWWYwZM4addtopvuAZx2vXrvY5iDgKQnXObdvGr4s2btw4unfvbsxjbdiwYUycOJFhw4YVXI9hw4bRs2dPI8dQUY4XU5S/5CRvoCxYsKCiufmT6i3edNNN1NTUFCTLO+igg8re72677cZll10WeX2vcnDlchkyZAgAEyaEl1kfM8bp40TtjVc7WTO36I7Q45iYVAoibsTynDlzOOywwzj88MNjba/CjUupq6sruB49evQwdgzlcaOuWMKLyWIIdzJz4cKFFUkul9QI4uKLLwbgscceqy/dmIV0Dt6G4Jtvvom0TYcOHQpcXKuJNdZYQ3ubiRMnZipvVSUUhMk5Ineey2Rwnasgli9fXnQ9unfvnliN6vTf2Crlp59+SnT/uj17t2RkGE8//XSRrdUb/JUEXgX0j3/8I9Fj6RJVQbhmvGokqlns9ttvr/9cyYI0UYh6n1zimJhU3ny5XC7WCDuJTpd7TsuXLy8a7ey0006JFa2yCiImSUfD6qaFKPVSzJw5k4MOOogTTjihYLnXJTIJLr744opmnyyFG82sQ9YnMMOI2lh5R6vluAInwbnn6hWxjKMgXn/99aJluVyuYM4qKu41N5naxTWDfvLJJ0UjiG+++aYs77wwqkpBSCmvllKulrYclUC3F1JqBOHmb/n666/jihQblZ103rx5Ze2znF6abmGkUnMUWSZqD9i73uLFi5MSpyLEURAbbaQO4fKaq/bdd99I+0pCQbheZi+//HL9JPX+++8PUFT8ySRVpSAaE7oNYKmXYp999gGSNykF4ZY1dQlyUYxKOQrCnZxuDES9Tl4FkQU319deey1yFLifOHMQqnxLuVyu4PrpJjBMIjnkkiVL6kcQ/oJJSZi2rIJoIGQ9U6l/Ei3NEUTa/PjjjxU7VtTGvpKec1Ho3r077du3j7VtHAURxZW2ZctoSaSTeDbd93vFihX1CsJVWMcee2z9b6bRvpJCiGZCiPZCiA6qP+MSNmCefPJJFixYoPxtwIABWv7f99xzjymxEuHEE08s+H700UeXtb8s9HLjMnv27IodKyiP0MEHH1zwPWvzDuUQx8Skep78I4ggM5RqOzDbafMGcvoVhBtRnqqCEEKsKYR4BPgFmAV8FfDXIMjlcrz11luJ7X/69OkceeSRnHLKKcrf7777bu6///7I+0tjbkGHcoOvxowZw+eff17/vZpHEEnL3rVr1/rPI0eOVK4TxeutGlDNl8QZQTzxxBNFy/z3KaqzRRIKwpvOx68g3PNNewTxd+BonKR8NwKDAv4aBKNGjWLHHXdMLNvpwoULAQoCyfzoeCaMHTu2bJmyTO/evQvcTcOuW9ZJWkF4r02QOctU56euro4+ffpUtIiNlzPPPLNoWRwFEcWVdvLkyVrxMEmZfV2vOnek5J5vEhHwOrMuBwMjpZR9jUuRQdxoUm+vNWn8DccHH3yQ6PFU5/b777/HToedNF5300p62pxwwgk88IC5LPZJm8d22mmn+pTfSRe7efPNNxk5ciSjR49OJavx9OnTi5bFURADBgygX79+Bcv8JiZwzIMdOoRb0t1tgszHpnBHEEOHDgUcJSeEMHoMnSu5Mk7q7UaBq/2HDx/Ol19+aXz/ql6kfxSQdJlDb/pql0ceeSTRY8ZBda3KTWOg07uL2uB89913keaCklYQf/rTn+o/BxUM6t69u9Y+FyxYoExW6HqEmUxY6b3fpSb0VSMh3TmI+fPnKy0FqkJgUQJk3XXKdcQohduRmzNnDpDMqFpHQUwCdjcuQUZxG5CvvvqKjTfeOLHjeKNW/cP0pBsSVW6hLAW1uagURBJKO4ioymTvvfemb9++JSORk76v3usVdKzVVtMLJ1p99dVp0aJF0fIkggi9HlVxFI+ugujTp48yWePSpUuLnr1SowcILzC0YsUKYxUo/W63abu5ng/sKoS4SQixvRCiU0P2YtJ9MHVztLvmI6/Wf/bZZwvWSdpWrcqIOnXq1ESPGYe0J6SjKIh77703cuUx3cZZlygKKE6EMDh27pEjR9a/H0nY2cPMh65LZxh+BVHKTKyaoA4iivlVFZXtcuONN7LLLrvw73//O/Ixg8iagvgQWAe4CHgD+IIG7MV02223aa3vDwQrxemnn160zJ8wrHPnzqH7SCLrZtqNsYo0XVp79uzJ6NGjS67nTbdeKqAq6fPx15tWoVO8yftMjBw5kj59+tTn1UpCQYRdn/vuu6/k9q4DiIvX5FYuqvdj6dKlBccMm4txMxG7ySzLwX3O1l9/fQDWXHPNsvdZdAyNdW8koyU/Gwr+B7tUYM5mm21mXIYsKohKyLRixQrlXIM/XXgUSpk4kp449tZkNoH3+rtJ7ebOnQskU97T643jVxZRrt0ll1xiTBb/s/fAAw9w7bXXFizbcccdef/99+vXjeJ9GHcE58VVENdddx0nn3xyIvnhdNJ9X2386FVKUuUj/bbrNCvVZYmgHmUulzPWgx03bpyx/P2lGs00lPAhhxzC//3f/yl/K5V+xXv9vVlFIZkU7qXmUGpraxMbhU2dOpWWLVvWewP9+uuvBb8PHTq0SEG8//77Bd9dD7IwTCTXc/fhKookrkmggsjPJ8yTUv7m+V4SKWV1Js7XIGgSsq6urqwelVUQaoIa1OXLlxtzyTVZTrTUM5CGyezKK69k8eLFvPjii/XLXAWrchX14pXX7fkmWXXOe79Vx3nyySc57LDDjBzL9QByqa2tZdNNN63/7veiitIh8SsVFSYUhBs8l1ag3FeA9y58TXD0dIOagyjVwwt6OUz7ge+4445G9wdw3nnnhf6eRRNTkKI0eb3LGYn4lYv35VddzzSu8SeffMLWW2+tlOOVV8JrgXnldUfP7j3xKhxTlFIQJucU/F5J/hFRq1atCr77zcAqttpqq5LrmPD+atasGZCeghgEfOD7HuWv6ikVVRmkILzh8GEEVX+64oorCr77h64m8L4A6667rvH9J0HQ9Yp6vaNQTqPtNzl640vefPPNovXTGEF89913Rcd1z7lUTRCVvKp3YN999+XBBx8sQ8pivMfZcMMNgWieRFHxP0N+BaGbSSGXyxW8t0kWFquEggg0MUkpr/F9v9r40TNKKbtq0I2YP39+UY9DxZVXXqlc3rx584LvN9xwA4MHDy65Px28snfu3LkoJ09UV81K8tRTTymXSymNjbLKURD+bb0mpnfeeado/TQURG1tbeyqY95G+p///CcAI0aMqP/sMn78eMaPH0+vXr3iC0qhy7i3sT3nnHMAswrCj38keeGFF2pt78+fNmPGjLKf0SAToHs/s5KLqdFQynQR1Jh07NixpD0XgoPRKuWt43LUUUelIoNL1LTXQYnlTNa2MPlyzZ8/n+bNmwcGRKVhYqqpqSnqvESVw9sLrnR1Pe99cTMAJ+E55VLupLvfEcA/oe0l6rXs0qWLcnklFIReBYxGQqlCJWE3okuXLrEbgEr0LHXr+ybJ/PnzI7nmBc1BmMzHVM6kq/9+v/nmmyxZsoQhQ4YoU3unMYJo0aJF0eRp1Oe00qnAvXK57rRAfX2Itm3bJnbschSEqt0IG7WtscYa/PDDD8p34KeffuK3334LrYnhuvzaEUSFCbqpRx11FNddd11iL3glGg6vTVW3QpZpok4MBzUIJhVEnHgHF79y8V5X1SgpDQWx1157FS2LU22uEtx44431n1VzAM2bNy/wNDJJKQXhNwN7ad++fVFuKHeeIIggj8i11lqrPgAuiLQnqRstQQ3nE088wRVXXGE8EKmmpoZZs2ZpRUKaMFOcccYZZe+jHKIqiKACO2HzPUOHDuXuu++OLItu8jov/ufFTbXw8ccfK58l3Xs3a9assu93TU1N4CR1KZJ0aVXhVao66bVNUMrs89tvv/Hkk08G/u7PqPDQQw+VJU9YvWn3mU0y3bdRBSGESDZEtEKU6kX07Ws+4/mjjz6qNXSu9EubJhMmTKj/7L03nTp1YsSIEXz44YdF2/Tv319LAUad+LzpppuKlvmfFzc1+FdffaXsQer09CZNmkTHjh059dRTI2+joqampkhZRfHXh+Bnbdttty1LpiBat25d/zloPjDoGpYbOxTFc/DII48s6xheSnWSwka2a6yxBvDHnEzaFeW+FEIcHPL7cUC8KuMZI8lJsCDWXHNNLZ/ychWE6VGQl7BcQOuss0795zixB9tss03951wux5lnnsmWW26pvZ+4qPLth/XEVSMInRfZdUGNkg+qFP369SvIARTVBdNbl8R7/VUeWibwKoUgD6Cg5z9I0UfN/puk2TUokj2MsLYo1UlqIUQbwJvspxOwvRBivmL1WpygunCDW5WguilJe5506dJFqxaDKj1xFFxb9L777htr+yhcf/31gb+VssmWwttbj5J6OSpRc+PoPhuqBkvnWTL50jdv3pwhQ4bU2/ij9ra9aehVKb9N41UQbibWESNGFKxz7LHHcsMNNxQsC+twbLzxxpGuu6sgNtxwQ+P1WAYOHFi0rJwATXfbtOYglgAPAf/N/+WAyzzfvX8TgCOB/xiXMAVUjYButlZdcrkcJ510UuT1DzzwwFjHcT0m/EN3k0FzYaObuL74LlFqHcTB39gE8cUXX4TK5Ec1Cakjt9fcUg4qD5uoCsJ7zm3atDEiTxje58dNPz9x4sSCda677rr6zyNGjDAW0e2++z169ChYvtFGGxV8/+KLL7QViKp4k0pBeCPz/YpRRSoKQkq5EKfM6KnAaUANcE/+u/+vN7AfcKJxCVNApSC8D2Qp4ph/VOUN3dTAJnDtlblcjiVLlhRFkJpM2xzW8HizccaxFyelIFTV9VR4S4+69TT8jZcXVSOiI7e/oYqL97q5z0LU0YDXY2j8+PEl1y+31KZq3sEfLOkdSZ555pnsvffeZR3TxX33/Y4Re+65Z8H3TTbZhI022shITiU/uvW9U4uDkFK+A7wDIIToCDwppTQXnZRRonqenHTSSYwdO7Zosm/RokWRIqq9TJo0qcjnefXVV9faRxjuw7NixQpltknvOS9cuLCsojZff/114G9bbbVVfYBbnIhY70uQtPulap7G23lo164dgHbksI6JSTUBH4Y/klfFLbfcwumnnx457fgGG2xQ3/vt2bNnyWylq6++Ol988UVRrzsqKgVRqcSV7v31X5ugZ7XcGBHVM6ybQiYTbq5Symsag3IIQnXx99prL+UEaalI7GOOOaZo2eWXX150DJO5mNxRTS6XU8rn7ZEF5T6KSliE8+67717/Oc6EoLdxjZMWRMfzRuX66s0i6o66dJMG6rzIjz76qNa+g2piezss7uR01HkXb7qJKKmsIbgWdhRUysBbkClJ3GfSrxCC4h+imJ51R1S6o5JMKAgAIcRfhRCvCiFm5L2a/H/FBtoGgqqxrq2tVaaBKNVgBL30/p5luT7UXtzGIJfLKQu9e4sPlWsiCOPkk0+u751FcbP0+8F7X4ITT9S3aJZrSuvWrRsbb7wxvXv3LrmvIHdpnRc5rOdcV1dX1BAH7dvbwLmeTP7ny2SivXKcOlTnXK5zQ1SCFETQCCLK3JVbZElFUIU6HTKhIIQQVwHDcTyb5gEzFX8NohZEx44di5apajWrfMshfpbRcj2lgrb32+1VroN//vOf6z+XMyFf6hxWXnnl+pd9+PDhJffnT69c7jUqt+Z2bW0tn3/+OaNHjy6pIMIKHUUlLHjy6quvRghRYLOOsm/XjOKfKwtKUx3nmsf1sgO1Db5SUf/utfErhP79+yvXj9Iohyl51bXVrTYnhGDXXXcteIdNoTOCOBV4GWgnpewmpeyh+jMuYYapra1VPrj+egBBkcB+ym38gnLVexuCCRMmKI9jKkNmlHKPrpx33HFHyXX9vdy061V4lULc0UiQGUjFqquuqlz+22+/1Re+945ig66PdwThzp20aNGiYL4oKE9XnGv+8ccfK5cvX76cG2+8MbRA088//1y0rFKxSe777Pe2K8d7K8xpZddddy1apvKUC2OttdbilVdeCc3bFBcdBbE28JCU0vy0fZVSU1OjdNv02vAvueQS2rVrF2myMWwoGoWgXpa3l7NgwQLlC2/qBRw5cmTJdXSGwv6YiijXUbV/U14uXrNRqUI7QYSlT/ATdK1atGhRH8DmvXfuvfXnKvLGjLhxEPvvv39B5yVIEXjrlESdeB4/fjyHHnpo0TW6//77ueSSS7jmmmuU2wV5k5kqbfrJJ5/UX1PV+brvUNQJ/L/+9a8F33faaaeidcIUxLx58wq+33nnndopxpNE56p/BCSTIatKefzxx9lhhx2Klnsn1NzUDG5vL4wBAwYULZs5cyY1NTX1+X3CUD3wUkqeffbZgmX+7xMnTjRm402yQEpUVEP68ePHk8vllD02HbyjBp3Axrj85S9/KbmOdzLdNW+GTRK77q25XK7gWq1YsYLjjjuuaH1vZHzUiGRwIoe9TgnwR3GloJxHQQnqTHRgpk6dyp///Gc6d+4MqJWv+/xGfR+8mQH+/Oc/07Jly6J1dDyw/va3v0Vaz6+YkkJHQfwN6CuE2D8pYaqN1q1bR84q6a86FpUzzzwTgJ133rnkuqoHvnPnzkX1e/3Jxnr06FFQxjGNbKMmCXMSCDMLRbGbm+rJvvPOO1x22WUl64f07NkTILL5IMoIy70G/upndXV19aMbk9X6vLgNva7bahwF4c/F5eY1mjFjBqB+zt0YkahzHmPGjKn//PHHHyvldM81SuGgNNL8hKHztJ8PLAT+I4T4VQjxdUP0YtKxtx5wwAGRg43CbK5h6LhQRpE9yF3Pq0SynAhw4403LrlOXAWhSokdtr1uymmvOXLbbbflhhtuKBkI5zZiUec73HP3mjoOP/zwgnW8CsKt0gaFjXZScz1uA6irgOI0nP4Eh373a5WCcHNtRT1/fyCkqgPhzkn6R1PlYDKwNQwdBdEMmAFMAqYCX9EAvZhGjRqltf6hhx4aab0oxc5V6Li8Ren5t2vXTrmed5La1AjirLPO0lr/l19+oVmzZqFpE6LIphuXoIO3AXAjqaPQrFkzZSbYUg1l3IbazWEExZPPXgWhWu7/bBLXi8zb845CHAXhfw78+1A9S65pKe51V1039x1WmZ/8ZMFE6yWy71iaHkpCiD7AQGB94D2gv5RyShLHCkub4Kempibyi9StW7dY8ugErUV5qKN4SJiKWtXN7/TBBx+wdOlSrrnmmsDefJRzDJO/3IbP6yadvzIAACAASURBVEigMwdx3333KWUvNVqLq6y9NRXefvvtgt/ca3DzzTcXLPde86RGEFFyC6mI4+bqv9f+3n1YR8mkgnA7AWuvvXasfaaJtkFVCLGyEGInIcQxQoh1hBAthRBmMoqpj3cSMAJ4ADgCmA+MF0JsmMTxdMwrX375ZeSUGpWw65s6hqkeuDe1dBTcRkDXb9xPWK+8XAUxePDg+s86/uqdOnVSxhmUet7iKBUIt3e71+CJJ54oWF5OepUwPv3007L3UW7UPRR7j6muo6sgguIewujTp49yufs8htV2iFqfHRzTdtTJ7HLRjaQ+CseM9BpOptc/A92Bb4UQxS44ZSKEqAEGAXfnU308i5NA8AfgAtPHAz0F8dtvv9GsWbNIjVaSdn3XV143wCYIUxOUq6yyijK5m8pTBv4wAQQF6tXW1ka61mEOAeUqCG9DpbOv2tpapetkqefi/PPPLzqWKn4kl8sVxIwExU/49xVE1GJCUfBG6cfFxOStv0StqiPi3t/11ltPe/+qhJvwR4dryhS10eOZZ56hTZs27L9/NP+fp59+uj6WJWl0Iqn3Bh7GmYe4CCe7KzhzER8CQ4QQJxiWbxOgI1DvIyql/B14BkikoIFOQ65jirn11lvjiBOJTp06AWilCw/DpAeLqiC71zXQi5viI8gXPugF9JOkgvC6P+p4NAVF3Zd63lTKMshM6E0aGGbei3INknxe4xBHQZQ6zzAFEYdcLqccxbsKQuXGDn+4wD/33HOxj50UOiOIK3Emp3sA9TNMUspPgV2A13E8nUziuol87lv+JbCxEMK4T5iOgtBpSE0Ms4PYb7/9AHjjjTeM7M9kCuMuXboULQvqKalqFnjJ5XKB0b5ewhSE35VVV6nGtSPX1NQoCxyFPUP+eJUpU6bw6aef1leZC8O/rV+WUiT5vMYhjoIo1ZlQvevldCAWLFigvJ/usqCJ+agmXSllpOffJDoKoivwsJSySEVKKZfjmJyK6zGWh2vg9493f8WRPXgcHZOo2SrBbENaDq77pKm5A5MjCJUbcFB+IVPZa3ViTnRrPe+zzz664gDRU7dvv/329UkI/Q1Kt27dAs01/onoMPfZKI2g3xyTNnEURKnG1LTZ9/HHH1fG0rjv05tvvqncLmrK8E033TQwkDApdBTEMiAsYc9agGn/QvdJ9ncF3OXGZ351Jnor+RL17ds38Dd3WFxXV2ckE2tSQVIuQT07U95TOgpCN/DNG8ugk2ojSvwGOG6gblEi77FKNXalAu68RFEQSTtVqEyPYegqiH333bdkNHSp52233XbTOmYQSbpdJ43O2/EycJoQouiqCyHWA84CyiskUIwbj+93IF4NRzmYmZWNyeOPP17W9joJ8u6++27l8s6dOxfsx0SRIXdkVMrkE5egxseUgtAZ2YVlS1URNUdPHPxOBjrpT6L4+LuoFPQtt9wS+9hR8R5Xt8HXXX/AgAEl722p989U1HzSHa4k0bkClwHrAR8AV+D06g8VQvwd+BhYHbjKsHwz8v/9GcI2AqSUMtXUnrojCH/v/thjjy1bhqOPPtpYJlaXRYsWMWbMGNq3b29sXsNLpRVEWI9Z9dvs2bMDt0lSQfjdTHWqyfkTPYYpCFUwqBth7+YpMn2ey5YtK+jkeCeD582bx9ixY5k5c2bBNt563LqTxzU1NSXTXw8dOhSAXXbZRfm7KU+hRqEg8pPR3YH/AefimHn64UxMfw7sIaV8z7B8M4BvgPpwZSFEE+AAYILhY2mjOwfhfwGilId0Ofroo5XLmzZtWhA16yWs9KeLKnPkM888w8svvwzomS6iEtT4hlWi08Gf1sR73aOMGMIS3ZlWxmEEuUWquOuuuwq+hykIVaI815ziKgods0iU3n3Tpk0LEsx5R6cHHXQQhx9+OLfddlvBNl4lFVVBDB8+nG233ZYddtihpBnL9dRSlZYF+Oc//xnpmKVoLCYmpJQfSil3B9oAOwLdcOpD7CClNN7VzI8QhgB/FUJcn08U+H/54w81fbxKojuEDzIdrbLKKkyePLlg2bRp04A/kpIFcf7553Pttdcqf3Mf6iQKtWy//fZAcUW4KBlro3DZZZcVfN9ggw3qPx955JGh25byfPHet6C8Vmngj1sIG92q5G7bti3wh1nl6aefNihdOG4+I39P29uRiGpi2mmnnZg6dWpoHEhUoqTGiELcEYT3uU2LWEY2KeVPUsq3pZRvSCmLa24aREo5HBgAnAg8AawB7COljJ53OCF0ewYvvPAC4CQEO+SQQ7S2DWq4VlpppaLazNtttx1QeqLxwgsvLGgs3EYC/ujh6faYoxRWqa2tpVmzZrGCkaIQFEfhHtuL/7rmcrnQ6Hiv2eLAAw+MKWFpdFM9+EeLEyYED7DDOife61MpJwy38ffnIfIWQoqqIPz1FaZNmxZJkU+ePDmwxnkcc9sVV1xRr+DiKogsZFXWCZRrLYQYJYT4RgixXAhRp/gzY0T2IaW8VUrZQUrZIl/NLpE8TLq4NsyouA1XXV1dWRNg3sYjzL5eyo3PL8PWW29d/9nNSaU7glA1bHvssUfRspqamlReAP/18iv577//PjBD76uvvlrgWWRqElOFqfkYoMgWr6pw6JrevOcUtRceN3bATVviyhOW2yqqgnjssccKvm+zzTbcfvvtJbfr1q2bMmbnl19+0UqD4TJo0KD6zpX/GYvqwZV29UTQG0H8EzgF+Bon5mGM4i+6Ub0K+Pbbb0N7ibomBveBqauri5064IcffigwzYS9nKUaYP/oYM6cOUXr6PaeVA/1NttsU7SsVNqMX3/9tcgmHYVSozr/9fJ7ai1YsCDwuvknM5PM3V+uj753PsevbDbcsDiNmdsRSCIgLYjLL7888rrlXOty0my3atUqdo4q9/3yjyBKmTldXIeBNNHpHu4F3CmlPDcpYbJG+/btlSVFXXTtnOussw7PPfccX331VWhlsx9//LGol+G+hF4zEIRPupZqZPz72mefferLWLqYmJS95pprijKHlhpB9OvXTzsl9BNPPMFRRx1VsMzfe/Y3Zn6Fsnjx4sjzQ0mOIKSUZW3vPW9/sSmVgnBJ8pzKIaqCUBXWSsJl189WW21VFOi5yiqrsGjRoqJnLOqIKwv3QjdQLlvx9xUgrNBPUNK5IOrq6urTTITdfB3f/KC0Ff3799c2U6jcbnVMTD///LMyn7070vI2TKVGELrKAWDs2LFFy/yeUX7PFH8HYPHixZFde/35jlSpNOKim3gxLHJ6+PDhBd/DlL7/ufSWG/XinRSvhCkkqoJQdegq0dCqcl+5sui4XXupNgUxGjhJCGHerSXDhNXg1R32envMQa51QQS9hEH2zKFDh/LOO+9oHUP1cqka3SDCFNvMmTMLGhvTcxArVqyIlfrE37scM2ZM5DxERxxxRMF3VcF6KF2PvEePHkXH1J37cRM2qvDf17B9+xulK6+8UrneaqutRu/evYHKTKZGvR6q9+Hnn382LU4RqkY/aA4iqoLwrpdW8kSdp/AKnCyqnwkhngPmUpwCIyelVPtNVikmk2N5TT4zZsygtrY28ssVp2Sprh1bNTkb5Nmhi793HTV1d1RWrFgRK915p06daNeuXf1cxLfffssRRxxRFFmswvsCv/baa3Tt2pVHH320aL2DDjoodD8vv/xy0VyX7mhE9RwdeeSRykJV3bt3D9yPv9MT1jC7519TU5P4KCKqgvAG17mYSoN/yy23cNFFFyl/Uykhr4LwdgjjKAjd1CSm0BlBHA/sAXQCzsSJmr5a8degcIuYm8Bv8olazxripb0Iy+ipQjUZFzYHE0aplB/ljiDc3qvLihUrtBupBQsW0KlTpwK79Ztvvsn1118faXtvTeqdd95Z63768Y9Ur7jiirK2BycVzAUXFJdN2XbbbQu+e3MO7bnnngW/jRs3jpYtW7LvvsXZ9d0GLGjkFIWgBtfvMBF1tO66eHsxFcm80Ub+hA5/oErE55YZnjVrVoFSjmo68q5XqRrURTJorHsVTsT0AcCfgA0Vf8FXsEpRVQGLir9R9yuIMLOAnzj2SJ1UDVA8aQ3xUy6Uaqxra2tZtmxZ5BGa1yFgr732YvTo0UXHU0UIh6EKhPrpp5+K0lZEkcnFVIqKkSNHaq2vU5rW39h6Gx9VxH6zZs1Cn9VyPIyCTCf+tOpRj6FqSDfZZJOiZd4R+Y033hhp37rMnTsXgEmTJhUsjzJpPmLEiIJzSWs+Queo6wG3SSmfk1J+IaWcqfpLStC06NmzZ6ztpk2bVtRD9pt8dJL9hT0gUTOFqnqTpUgqtcSSJUu455576NChQyQb8VVX/ZHmSzXyyOVyxvJGRVWsqsZIp051Wvjl9n5XKfZcLlewjut+6S6rRONVjhJSTSB7LQMmoq51uPLKK0uW4z3jjDMKvlfDCOI9HPNSo8KbP0aHhQsXFt1Ufy+vY8eOkfenilFwiZKWeM899+T444+PfDyXpEYQ7vAb/ki1EIbXdKDadznmqrj2c1XDGGbfh/JGpEnhPY8oCsLtEXvnIJLGtBLymp2iym9Kkay22moMGTIk8PeDDz64aFk1KIiLgD5CiDOFEOsJIdL3wUoQt6C796EIyvqoolmzZiUVhI6nSlgiuyA7rpcXXnghVl6lJLOXuoQpP5dSvVxVaol+/fpFOv5f/vKXSOuFyRSEf24gSiqSSuM9j6DRmXcdd17KbbSz4I5ZDlHNc3vttVekan7l4o6WvJ5t1aAg7sLxWroT+Bb4vVKpNtLA78II6vKZQXz66aclb6qpSNwoctXU1MRKLlcJBRElXsNrnlNdN1Wvy00KWAqV50sUVPfXv8wflObNL5QVvFXQVJ5vP/30U8F5+UcOlVAQSbrSRnUuqKmp4eyzz468X7cU8J/+9CetBt5dN8zFvlLo3NlpwH9Qp9hokKk2/IS9CP6I3fnz5zNu3LjY+0sCIfQrwkZJGV4upa4TFCrBqCUat9hii0jrRbkPqnQNpRSEOwrVPVaaRGmI3bomlTQxRUmMecMNN8Tad1JzEFtuuSWgryBUz0hQudKk0akHcbKU8pRSf0kKmzZhN/ndd9/lvvvuq//etGlTzjrrrETliZrTxcuf/vQnrfVfeuklfv7559CJW9V8i45d/9577y25jrd4S9S6vF27do20XpRGO6zkaxCqUWhapgI/gwYNUi4Paoi9crteOG7nQXX9hg0bVqaEhXgTSQYR181Yd7vnnnsu0npuzrTevXsrFe8JJ5yg3E71jKSV2VW7OyOEaCOEOE4IcbEQopMQoq0QIrrtJeOoMo+6eCOF/bUMmjRpUnATdUtZ+vHnanLt5N4GPk4hksMOO0x7m27durHlllsGRn/HyXapQznJ1qIQ5h7bokULcrlc5LQqpby+sqIgdObTAOUcxPjx4wF1nIGOOTYKUXr5cdPH644g/DEhb731VsF3d+TomkKDGvcgr8IopstKoaUghBAXArOAB4HBOHEP3YCPhBDDhBDZePrLoHXr1my22WbK34466ihyuRyvvvqqska016RRqthIqR6W33fbfeG8Lq1xAoDiZKZ0q8qV8tDxousZFPYC+BPwuahMOFHxXsewyXtV4FUYLVu25NZbby1K3OaSFQURlLspaM7Jez/910s1gjA9dxXFwSJOqhWIP/Jw8c91tW/fHvjjurz44ovK7VRZjqFKFYQQ4njgZuAp4CickqMA7wBjgb8C55gWsNKsWLGi6IFftGgR48aNqx/m7rLLLspgF69Jw31IggiKUHajUv2Nv/uCel9UVebKUrj+33369NHeNi2CXg6VCScqXtfgsOymccyE/fv3r7c/+zGVuiQpghT7HXfcUd+5UAXa+VNBxJnv0t3Wn6hSx23cS7kKwo+rHN3nNooJ1UtVKggcN9cXpZQnAC+7C6WU30gpjwSeBaqn1QlApSBatGgRqQJcbW1tfXK7qAnf/LgxAcccc0yRXFD4Al966aWxjgHxU2ioUPUgVcFJcQl7OXRfbjfnkTcmRHeUUA4mC9ifeOKJkWzzOoTZuqdMmcLgwYOLzGi5XK6ozGs5KbajjnL95XJ1RrheTJeOdUc7cR0S3Ofd28msBgXRBacedBD/oQGk2lApCB3cGxnVZj106NCCB911g/SnGnAVg1u2FLLjEaOS47///a+x/Ye9HOeff77WvkaOHMmll15aYGIJmzfIiklIxZgxYwLTccclrOOw+eabKzslixcv5oILLiioVLfaaqvx0ksvxRqpRr3mpmJKomYiCMObO8t9nuI+Ow8//DBAgam7GhTErzj1oIPoCCwM+b0qmDdvHkuWLIm9vW6jff755ysbOb9dMwv1aaGwDkAYJguue6/pOuusU/Cb7ghinXXWYfDgwQX7DLNvxwkujIrJRJCm8CYgjMo777xDTU1N/ajRTcWxxx57cM8992jvL+o99Sr2ckq0xjHV+vGOmF25yu3Aeedx3nvvvbL2FRedM3geOEsIsbb/ByHElsDZgHo2poqYMmVK/aRsHOI8FN4XwnXn9PcYylEQql5mUI9k8803D92XyuMn6VTPrqwff/xxaER5XIKUwIABAyKZFnXw2smDai2kjd+8WQrv5PC0adOYPHlyWcd3e9Cl8CqIcoJOTYzEvbVXyh1BuHi999wa8ZVG58pcijMx/TFwN05U9RlCiKeAt4HfgWw+8RUkzkNRW1vL6NGj2XrrrQNNM+UoCB079R133BH6u8oEMWrUKG2Zgthrr72KlrnXdLPNNkskVUWQgrjpppuM15329tAXL15sdN+mKKfB3Gabbcpy8e7SpUtBzEsYSY7uSuFXgt731tQI4sILLyxrexPoBMp9B2yHUzSoJ46yOArYG2duYkcpZfqx4SkT9aHwT2L37t2bd999N3CIrypEkwSlFJxKQVx99dXGjq9yMQ6TycTopZJzOd5UFt99910ix3BTPJQiqChRmnNbN910U+R14ygIU7mU/B5T3rQYOiOIsLm6LMx/aT0JUsrZUsqTgTWBdYB2QGsp5dFWOThUeoJNF5W7rJc4CsI0/h5onAZLJ0aiki+i11bupqwoha6Xjb86XRD+mhouaSoInUbfXVcn0K1Udb+ohHlp6SiIsCJEXvyBuZVCN1DuICHEG0B7KeU8KeX3wDAhxDQhRHDV9EZEVjyLgjj88MPp3Lkz/fv3V/5eyuXTpJumin/+859FJSJ1a2tD/KjapPGOIB588MFI2+g2alGVuH/C3yXqM5xEZ0FHQTRr1oyrrrqK119/PfI2puqbhAUC6hwjqgnTXwWwUugEyh0CjAPWArzq8zWgKfCCEKJ0YYIGTtYVxPrrr8+nn34aGBxWyoPETc3922+/MWrUKOMT1MuWLSuKiA0rKBSUIypuVG3SxFGwumkxfvrpp0jrBXn+RH2Gr7nmGiB6WpnTTjuNs88+O7QWgq7Z6Oqrrw4MSvRywgknsMkmmxibtwhTjjrzVnHKj1YSnaNeDrwKbC6l/NxdKKUcA2wNvAFcG7BtoyFJc8WOO+6Y2L6j4r4YAwcO5PTTT4+cuKwcXnnllcDfPvjgA+Vy3YY4biZQXVTlL0uh2+uN+gwGKfeoDZybOSCqu+4555zDnXfeGVpNLamJ5/vvv58ZM2YYG0GE7cc9hyj3Ieq1rgYF0QV4SEpZ1DWTUi4HHgKyVy6rwkS9kXEUiT8pmJdKzA3AHw+0O8Hqre2bFGF1q4MSFurKdckll2itHxfXdq1Tc1q30YzSo4Zgz7ioz7D/WVDhdZmNMlfgPddyorGD8DfsbsyGLmHvr3v9ooyuo17rPffcM5pghtENlAtOWuNMWGdzXB8RE+YSkyMIHbvjtGnTjB03qJA8/HGN3IZZ1ZPaZ599jMlSiiAzSSUUVxzc6xf0nKhSQOsqiKhBikH7jdpouVH9Yea8zz77rP6z95yDju1drnJ5Lhf/hH+5WZdV6CiIVq1aRdqnycBTHXQUxHPAOUKIovqMQoiuOIn6njclWBqYiFaO+nJFeXiCssqqMJleOeyhda+RV0H4XzI3DXQYPXv2LEPCP5g5c6ZyuYno2CRw5xOC7q2ql61rFunUqVOk9YImqY8++uhI20cxLwalAglS7F4FYcoc5MX/fqYdzb7KKqtEaguqIdXG34CfgclCiLeEEA8KIe4XQkzBCZRbiBNMV7WYUBBRb6SbtTWMESNG8NRTT0XaX21tLZ9//rmRkpZ+LyIvdXV1LF68uL5nuMoqq8Tq3ZiMnVCRdHR3XE499VRmzZoVWAe7e/furL766gXLdIP1oj6D/nxfLj169Ih0/ZKYG/GalZJQEH781zouKqcPk89g5ucgpJSzgS2BW4AWwCHAkTgxEXcA2+SD6aoWVT1eXb799ttI63nTTQfRokULrQI/G2+8cWCvUIcwRblixQoOOeSQ+qyzTZo0iVW4KAn7spc4MlWCmpqaeoV6wAEHFP3erFkz2rZtW/YxTOBP4e3nr3/9q5HjePEGoFVCQZQz0txjjz047bTTAPWoW+VqHdQxKEVaIwgt46aU8mfg4vxfg8OEgpgxY0bRsg4dOjBr1qyy950FVqxYUVDkvkmTJkU1DqLUaUi6h590vIYJ/vWvfxX14ps0aVI0YtB1QDB1bUtVCgxzmlDhbeQmTZpUVDXRv04lFIS/lrwO3vdAFfiqmmeJYn5VUQ0mpkCEELVCiNWFEPr1LDOEiYlN1VDwgQceKPiumwwtS/hHF6qXeO+99y65n6QVRKVGEG+++WZBmgUd2rZtWxSzsPvuuxcpCF0TU6VyFC1cWDp5s7dErveeB81zVVpBmGp4o9bgjjop7SfzCiKvAEYJIeYLIZYLIercP5xEfT8B8WtAZgDXbFIOKgXhDz7r27dv2ceJyoABA7S3CXsY/aMsVcK5008/veQxdBSEzstxyimnANEnaoOI6iCwww47hFakK0Xr1q0Lvrdq1aroGVI9U0FzOF999ZUxBeE196gSMp588skl9+HeD4jWAYvi6ZRFyp3LuOiii0J/T+ta6IwgbgROAT4BXsJJ1vcQMAFHQSwB4teATJFly5Yxe/bssnLKu/gbs6lTp9K1a1fOPffc2PtU1b+OSljUahADBw4M/M2f48g/l9KkSZNIE2o6CuLMM8+MvO7xxx9PLpcr+4XVjV42yd///veC767bsDsn8Pvvv3PVVVcptzU5mek10+2www5Fv0fpDU+YMEHrmNU6gih3PzfffLMROUyj8zQdCIyVUnYDeuWX3S6l3Bv4C1AHxIs6SZlevXrRrl07I3MQ/hd02223pba2ln/84x/1y8ICv1SU86LEaTDCzDPeqmEqoo7CdJRxmInF3wMvp4GcP39+/edSE7RJ4Fa58/v/N2nShFwux9y5c1m0aFFob1LXHPXwww8HRqN777WqTkip2iGgl2oeCp/1SigIU2Qh82oS6LxNawMvAEgpfwS+A3bMf38PGAUcH7h1mQghWgohZgohjjS9b7dX7JYeLCfQy/ugBLkR6mK6JkE5hCnRO++8s6CObhg6CiKs0fd7eZWjILyjDhOjSV3CUlCAc26lcmXpnv+xxx7LFltsobWNS8uWLQECRzMQHNUdxTW6EgpCJ9YojCAF4aYjqVZ0nqYFvvW/ALxP1oc4ZUeNI4RoiVNzQp3A3jA6rqV+vOVKg2yuuuaPalEQUWzSUfbTq1evgu86jV6czK8qVB421YDJZyVKKupcLhca0+Id7XjNiv7gyl9++aWoRkpSCsIb4Gcq7X7QM3r//fdr70tlzksLHQXxFnCMEMJ9Aj8GugshXNUpSCDVRj5D7Fs4CQETxfXQKWe46M3xH7V+cymmTJliZD8mCCvwrtOQh/XQ/Y1cWM/aP5fxv//9L7IMYUStqWASldJ0M6ZGxeQchPsevPrqq7H34ZUnLACzVatWRXmRklIQ++67r/F9BuWZiuNGe/vtt5crjjF0nqZbgF2Az4QQrYHROHMOLwghhgPnAcHlkeIzDmd0Yv6u+nAbm3JeMm+QUzneLV507bhJEtazNqUg/IQF//mDyqrZFqya+wmrO6DCpIJw5yDKMZV674c/XqYU1TQHYTLwM+kgUh10IqknAvsDEpgvpXwbOBPYCfgr8C6grkJTHt2llEcDcxPYdwGueaicRsZb+cnUcL/SE6beHky3bt0KfgvzPtK5bmEKQmc//h52NQTIBaFKeqfb4Js0MbkuzEH5rqLglV/X9biaFEQUorq3e5//u+++m8cffzwpkUqiG0n9AvmJ6vz3u4QQo4HmUsr5gRsqEEI0AYLtFTBHSvmzlPIjnf2Ww9y5jg4qR0Gsu+669Z+DGivd/ateelMT4Cq8tme/rGENsE5jFjYHoeMC26xZM1588cV6z5+sZnGNgnf+yiVKttGzzjqL4cOHA8nk7Ikb3AV/1L9o1qxZvZdWVKopDiIKcZwB+vTpk4Ak0Sn7aZJSLtVVDnnaA5+G/PUuVzZd3JernJfM26D6e8lu7p1ye4UDBw40EtQXxH777cepp57K119/XSTr9ddfH7idjuIrVdpUB2/DU82NildBuBO/3kCzILxRvCZHEK4HUjm1Rrp27cqMGTOUAZWlqGZzoYqobvRZOu/At0kIsQLQzYeQk1JGekOllF/jBNtlBjcTqqkb5Pbq/JQ7gjj//PNLujuWQ21tbX3krF/WMHODjuILMh/EaYy816fc/P5dunThhx9+KGsfcfGamK666qpQ99EgTDYup512Gueddx7rr79+WfuJU0UPsuW9F4Xp06eHFkWKmv6lKhQEMAZ9BVHVuHnwTQ3Tg3LN6+Yh8j8wqiyRXp544gmOPNJMuIj/Wpiag/AXbnF57733yir/GSWNehi6E6kmMVFH22Tjcs4559CvX7/EU02fddZZyuVZr+/uRwgR+nuDUhBSypMrKEcmcG+gqRtkapJN180wSjbVqPjNZKaS7LVu3Zrp06dz7bXX8nJNbwAAFQ9JREFU8uCDDwJOeovOnTvTo0cPxowZE2u/DcXEFBddr6cwampqEm2s2rZty7x58wJt80ke+9lnnw0tlZoEDUpBBCGEWB3YE+iEk17jM2CilLL8pzsjmLpB/pKh/fv355lnntEOhEmz+M1rr72W2L6FEKy22mpFx+rduzf77bdfwYR/VKqt1+nFhIKopvM//PDDueuuuwJ/T/Jc9ttvv8T2HURU1+6qVRBCiIuBK4DmFM4fzBdCXCClvM+kcGlh6sH021B79uwZq7HP0gNjmokTJxYtq6mpiV34qJoaSD9ZinepBG6Z3A4d1AkSqvleqoiqILJ03jrpvs8AbgCm4lSS2xLoChyD43V0rxDi4CSEBGdSW0pZI6U0nlL8qKOOKvhuqkE2daOzpCBOOOEEo/vz5r0qlfI4Cll6uaKy5557An/kNmosnHPOObz66qvsv//+yt/TrhdtmqgmJu+oOm103qbzgYlADynlWCnlR1LK96WUjwO7ApNxRhdVh3/CdOrUqUb2m5VUwibZYostlF4tqvKZUfD2msPSeEQlaPI7y4waNYpTTz2Vnj17xt5HNTamtbW1oWnVDz/8cLp27Zqq44BJSo0gFi5cyCuvvBI54WUl0FEQHYGnpJRFNhIpZR3wKGAmNWKF8VdJu/XWW43s11TDbiplhwl+/PHHyHW3o+D18zfR+w/KHpplOnTowKhRo8pyavjoo4+YNGmSQanSp6amhnfeeafeFFXtlBpBrLrqqplLEqnzRn4ChFX4/jNOhteqw68gskY5PUvTBBXvmTNnTqz91dbW1hdyN+H3Xo0mJhO0b9+e7t27py2GJQQT9WYqjc4kdT/geSHErcAQKeU8ACFEi/xvJwKlixFnEBPeI1ljzJgxsQOU4lDOqOKNN94AiosNvfvuu0gptfbVWBWEJbtstNFGfPnllxx33HFpi6KNjoK4H1gBXACcL4T4CVgGrIPj0VQDvOYLFokcWZ0mWY/YjGOq8iYNrAQmJtbuu+8+rrvuuvrvW2+9tbZnj1UQlqyx9tpr8+WXX1Zl8kGdxnsyDTSyOutDvyxNUgdhIse+iUysWVf2lsaHG7yZRpXCcomsIBpyZLV/DqJdu3YpSaJmxowZ9Z932223FCUJxkTvyM2mG4fff/+dxYsXV4UytTQu3E5Lg1YQLkKIFlLKxfnPawHHAsuBx6WUPxmWryL4g9fSiLIMw5sJM6sT6nED20yx8sorl5WW2mJJCncEkXVLhYrICkIIsQbwCNAa2FEI0Qp4B1gfZ/7hSiFEdynll4lImiD+RjfMNzsNvHb1rCoIN9GhxWIppJpHEDozetcBPYHn899PBTYABgI9cCawr1Nvmm38ja7JhGcm6Nq1a/3nNPMyhVFuSmiLpaGSpchoXXRMTAcDd0gp3ST1hwFzpZS3AgghhpFMydHE8Wv2rCkIr33frZyWNUzMQRxzzDEGJLFYssVdd92FECKz724YOiOItYGPoD6j6054yo8CPwDB1TIyzPjx4wu+l1NBKwm8E69XXnllipIkS5JFkCyWtGjTpg2DBw+uSg87HQXxHeAWKz4UWAl42vN7N2CWIblSJWsjCK+CaMh+/g353CyWakTHxPQfnAC51XE8l34C/iOEaAdcApwEXGtexMpTrrlk0KBB3Hefuczn1nXTYrGkgU6XbSCOF9NpwM/AMVLK33C8mM4GHgSGGJcwBTp27FjW9ldccQWff/65IWkaT8969uzZaYtgsVg86ATKLQP65P+8vAe0l1J+b1KwNGnTpk3aIhTQWEYQ5daTtlgsZim7ayqlXFbtyqFz584F37PWIKcpz7HHHpv4MQ455BAAjj/++MSPZbFYohM4ghBCfAmcL6X8t+d7KXJSyvKrvlSYbt26MX369PrvWVMQaZqYunXrxiOPPJLoMcaMGcMrr7zCRhttVHpli8VSMcJanpnAIs/3WfllYX9V6cXkTxORNQWRpjx9+/YN/G2rrbYycoxWrVpx0EEHGdmXxWIxR+AIQkrZw/d998SlSQmrIIJp2rQpN954I927d6dbt24Fv73//vspSWWxWCpB43CPKUG/fv0KvlsFUcjAgQPZbrvtUpXBYrFUnrA5iIlxdiilzE59zIj4IxyzlnUxC26u/tiQ2trazCYOtFgsZghzc92I4gJB6wDNcOIgZuCMQDoBbYAfgU/Ni1h5mjdvnrYIBaQ9glCx8sorGynwY7FYskvYHEQn73chxO44qTVOBh6QUq7w/HYcMBIYloSQlcbN354VsqogTjnlFO666660RbFYLAmhY7u4HRgppRzjVQ4AUsqHcZTDIJPCpUEW01ZnwcTkZ+WVV+aOO+5IWwyLxZIgOi3PJjhmpSC+BdqXJ076fPvtt2mLUERWRxDVWITdYrFER0dBSOBYIURRzlohRDOcAkIfmBLMkm1OO+00ACZOnMhDDz2UsjQWiyUJdIztQ4CHgdeEEP8CvsSZsN4UOBPoCBxgXEJLPVdffXXaIgCw8847M2SIk5exR48eJda2WCzVik6yvkeFEM1xFMUI/vBwqgG+Bg6TUr5oXEILkK1So82aNcvkvIjFYjGLlruOlHK0EGIMsA2Oe2sO+FJK+W4CslkySu/evdMWwWKxVABtf868B9PU/J+lEbLBBhukLYLFYqkA1k5g0SZrkeYWiyUZrIKwaGNTbFgsjQOrICza/P7772mLYLFYKoBVEHmeeuqptEWoGpYuXZq2CBaLpQJkK+mQAiFEN+B6oCuwGHgJGCClnGPyOFJKk7tr0NgkfRZL4yDTIwghRBdgAvArcBxwEbAzMF4IYTTPwzbbbGNydw0aqyAslsZB1kcQ/YDZwBFSyt8BhBAzgLeAvYBnTR2oWbNmpnbV4LEmJoulcZB1BfEx8ImrHPK4tqANTR5oyZIlJnfXoFljjTXSFsFisVSATCsIKeVwxWK3uv10k8d6/fXXTe6uQbLFFlvw4YcfsuGGRnWzxWLJKKkpiPwcwsYhq8yRUv7s22YD4BacKO5YJVGDyGJK7azRoUMHPvzwQ3utLJZGQpojiPaElyi9ALjN/ZJXDhNwJtaPlVJmJ3tdI2HMmDE89NBDdO3aNW1RLBZLBUhNQUgpv8bJBFsSIcTmwHNAE2AvKeUXpuWx0cGlWXPNNenXr1/aYlgslgqRaTdXACHEjsAkoA7oLqVMpChRltJpWywWSxbI9CS1EKITzshhDrCHlPJ/SR2rffuqr5ZqsVgsRsm0ggD+AbQCzgY6CCE6eH6bKaWcbepAO+ywg6ldWSwWS4Mgswoi7+W0P7ASoCp6PADHo8kIdg7CYrFYCsmsgsgHxxlNpxFGkyYVO5TFYrFUBZmfpK4UzZs3B2D99ddPWRKLxWLJBlZB5GnTpg0Ap5xySsqSWCwWSzbIrImp0rRu3ZoFCxaw6qqrpi2KxWKxZAKrIDy0bNkybREsFoslM1gTk8VisViUWAVhsVgsFiVWQVgsFotFiVUQFovFYlFiFYTFYrFYlFgFYbFYLBYlDcnNdSWA77//Pm05LBaLpWrwtJkr+X9rSApiPYBevXqlLYfFYrFUI+sBBcXYGpKCeBvoDszGKS5ksVgsltKshKMc3vb/UGMrqVksFotFhZ2ktlgsFouShmRiCkUI0QcYCKwPvAf0l1JOCVl/c5yKdjsCPwHDgJuklFUz5Ipxzt2A64GuwGLgJWCAlHJOBcQtG93z9W17NXCVlLImOQnNE+MetwVuBQ7E6SBOAs6XUn5ZAXGNEPO5vgnYCpgH3AcMztecqSqEEAcDD0opQxPHmWq/GsUIQghxEjACeAA4ApgPjBdCbBiw/to4jWMOOBq4G6fhvLAiAhsgxjl3ASYAvwLHARcBO+e3yXw1Jd3z9W27OXBpshKaJ8Y9bgK8COwA9AFOBjYGnhNCrFIJmcslxjlvDLwALMyvPxS4GLihIgIbJK/oHgBCOzEm268GryCEEDXAIOBuKeU1UspngYOBH4ALAjY7G2d0dbCU8lkp5XU4D9SlVdJYxjnnfjgT/EdIKZ+TUj4EHIvT69qrAmLHJub5utuuBIzC6VlWDTHP+SRgU2BvKeWTUspxQC+gJbBFBcQui5jnfCTOJOwRUsoXpJR3ALcBffP7yzxCiKZCiIHAf4HlETYx1n41eAUBbAJ0BP7tLsgPLZ8B9g3YZk9ggpRysWfZOGBNYPuE5DRJnHP+GLjVN+yW+f8le+EpE+d8XS4AWgF3JCZdMsQ558OA56WUszzbvCelbCelnJaksIaIc85Ngd+B3zzLfgRWy/9WDeyHM8IdQLTn1Fj71RgUxKb5/5/7ln8JbJzvQaq2Ua3v3V+W0T5nKeVwKeUw3+KD8v+nG5bPNHHuMUKITYCrccwtSxOTLhninPOWwHQhxFVCiO+FEEuFEM8IITokKqk54pzzgzhu7zcIIdYUQmwPnA+MlVIuSU5Uo7wNbCilvB3HbFQKY+1XY1AQrfL/f/Ut/xXn/FUl5FoFrO/dX5aJc84FCCE2AG4BpgITjUpnHu3zzZsXRgL3SylfS1a8RIhzj9sCp+D0tk8DTgQ2A54RQlSDw4r2OUspv8CZT7sIZ+TwFjAX5zpUBVLK76SU8zU2MdZ+NQYF4doZ/ZrXXb4iYJsgTa1aP2vEOed68sphAs7zcWwVeG7FOd8zcEwWFyclVMLEOecmwCrAflLKZ6SUjwFHAZsDhycipVm0z1kIcTpOR2AksAeOUmyNoxSrxcSki7H2qzEoiF/y//1uYavhXKxFAdv412/p+S3rxDlnoN6j53WcnsZe+R5Y1tE637wCvAk4D1ic7z3X5n9bWQhRDe9FnHu8EHjT2xuVUk7F8QTK/CQ18c75EuBZKeUZUsqJUsoHgP2BXXAm6BsixtqvangRymVG/v9GvuUbATKgdzwjYH34Y+I2y8Q5Z4QQO+L4xdcB3aWUHyQnolF0z3cPnBfmCZwJzN9xYgPIf74yITlNEucef44zgvCzMtFs22kT55w3AN7wLpBSTscxN21mXMJsYKz9aiwK4hvgUHdB3tXrABwziooJwJ5CCK9N81Cch+q9hOQ0ifY5CyE6Ac8Bc4BuUsoZqvUyiu75/gfHm8P79/f8b9vj+I1nnTjP9QvAzkKIdp5tdsPpgb+enKjGiHPOn+HE89STd05YC/gqGTFTx1j7VQ0TU2UhpcwJIYYAdwohfgYm4/j8t8EJmnGDadpKKd2exnDgHOBZIcTNOLEAlwKXSCmXVfocdIl5zv/AMSudDXTwebbMlFLOrtgJaKJ7vlLKH3FelnqEELvk9zW1osLHJOY9HgqcihMYdxXQArgZRzm8UOFT0CbmOQ8CHhNCjAQeBtbF8Vz7GhhT0RNIiCTbr8YwgkBKORzHh/hEHLPCGsA+nvQCVwBTPOvPxvElXjm/fl/gcinlLZWUuxx0zjnfC9sfJ6Doofxy71/mbbW697ghEOO5nofTm/4KuB+4Eyey+gApZTU4X8Q558dxIqi3AZ7FCRibBOwopfR7+lQribVfNpurxWKxWJQ0ihGExWKxWPSxCsJisVgsSqyCsFgsFosSqyAsFovFosQqCIvFYrEosQrCYrFYLEqsgrBUPUKItb1Ro0KI0UKI1P23syKHxRIXqyAsVY0QYj+c/DJtPYvvwgmkSpusyGGxxKLBp9qwNHh2xImmrSdfwD71qOmsyGGxxMWOICwWi8WixKbasFQtQojRQG/PoleklLu7y6WUNZ71tsMpEnQLsDUwG7gGJ/fUIJwKY6vg5CY6K5/Qzz3OZsD1QI/8Ou8Cg6SU46PI55PjLzhmp1twMsf+CjwKXCyl/E29p/oKeFfg5MXqiJPX/wXgMinlN571WufP53CcJHZfAiOA273psIUQrfLnf0R+vS+Af0gpR4adk6VxYUcQlmrmLmBs/vMFOI14EOsBTwOvAhcCy4F7cQre9wSuxVEWR+M03gAIIbbAMRNtBgwGLsepzPasEOKYGDKvjdOwT8cpWDQZJ/PmNSW2uwy4CngeJ+PuPTgpnF9wazHnJ+on4Sig+3BqL38E3IaTmM89p1Xy652Dc00uwFEk9wghzo1xTpYGip2DsFQtUsopQogPgMOAcVLKr0NWXxM4R0p5J4AQ4msc5bApIKSUS/PLtwb29mx3BzAP2EZKuSi/zh04dbr/IYQYq5lCuTVwrpTyjvz3e4QQn+CMDAaGbNcLeE5KeZ67QAjxDXAm0AlnBDAgfz7bSSk/zK/2TyHEYOBSIcTdUsr3cepRbwX0klI+lN/X3cAr+fWGSSnrNM7J0kCxCsLSmBjr+fxZ/v9zrnLI8xXQDUAIsRawG46SaC6EaO7b199xzESTNeV4zPf9fZza0GF8C/QQQpwHPCKlnCOlvAtnFOVyBM6IYbYQoo1n+TicegAH5o91II7Se9hdIV9r4USc0VFVpP62JI9VEJbGxBzP5+X5/3N969ThFH0H2Dj//5z8n4oO6CuIeb7vS3FqcYRxEU4lvNuAoUKIacC/gXuklN975G2u2L9XVsiPOPwlOqWUMyNJb2k0WAVhaTRIKZcrFod5abiN9jCcXriKj2PIod1Dl1J+IIT4E7AvcFD+/yCgvxBip3yd5ZWA1wiez/hf/v9K2FGCJQJWQVgswXyd/79cSvmS94e8Z9OGwOKkhchPQm8FLJBS/htn5IAQ4mgcD6g+OBPvXwMtFbK2BvbAqekMMAvYUnGc/YBjgYFSyjn+3y2ND+vFZKl23MlU489yvnTjVOBkIUQ7d3m+ROu9OOUcK9HJWgn4L455ycub+f/uNfg3sJUQ4gDfen8DHgc2z39/FlhHCHGYb70LgAOAH0wIbal+7AjCUu249vYBQojn8j1sk5yL47E0TQgxHPgROA4ngvtSb7xEUkgplwkhbgf+JoQYi+Pq2gKn1vBiHGUFTr3lI4CnhBAjcMxfu+C4vT6X/wNnYvtU4BEhxDCcVCUHAHsBp1oPJouLHUFYqp1HgJdwAt1uNL3zfLqMnXFGEhcCNwOrAidLKYeYPl4IVwH9gU2AW/PfvwR2zc8/IKX8CdgJGI3jFXU7TmDetcCR7txHPiBvd2AUjrIbCrQHjpZS/qtiZ2TJPDaS2mKxWCxK7AjCYrFYLEqsgrBYLBaLEqsgLBaLxaLEKgiLxWKxKLEKwmKxWCxKrIKwWCwWixKrICwWi8WixCoIi8VisSixCsJisVgsSqyCsFgsFouS/wfH5OuCwx6K7wAAAABJRU5ErkJggg==\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 0x7fc638a00898>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAECCAYAAADq7fyyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAY/UlEQVR4nO3df7CU1X3H8ff1ikoEQhx/BorAJX7tRFQ6TYQwpibEwXgn9ia2MQnq2HZoktGkiMGBqmNQGO9EiSYxxui0nbaQDA1J0QQiRfKzFBLtSPzR+O1FIChRqRWUIur1sv1jn8Vl2R9n9+7u82M/rxlm4ex59p575vJ5zj3Pec7TlcvlEBGR7Dgq7gaIiEhzKdhFRDJGwS4ikjEKdhGRjDk67gaY2bHA+4DngaGYmyMikhbdwGnAI+7+RvEbsQc7+VD/ZdyNEBFJqfOBfy8uSEKwPw+wYsUKTj311LjbIiKSCi+88AJz5syBKEOLJSHYhwBOPfVUxo8fH3dbRETS5ogpbF08FRHJGAW7iEjGKNhFRDJGwS4ikjFJuHjakNWP7eL2dc7v9x7g3WNHsmC20TdtXNzNEhGJXSqDffVju1j0gyc4MJi/GLxr7wEW/eAJgGGFeyMni5BjqtUp9x6gk5aINCyVwX77Oj8U6gUHBoe4fZ03HICNnCxCjqlWBzjivQXf+w10weBQ7rD6j/7uZX769P8cCvsPnXlSxX+/c+QIurpg72uDh/1dJwmRztAV937sZjYR2L5hw4bgdeyTFq6hXKu7gO39vQ21Y2b/T9i198AR5ePGjmTjwg83fEy1OkDZ98rpgrLfc70Kn9Pd1cVQLlfxdZxOAiKJ9txzzzFr1iyASe6+o/i9VI7Y3z12ZNlAfHcUlo34fYWArVQeekwjn1tOs06/hc8Zik7olV537T3AvJVbmP8vWziYO/zEclQXHMyh8BdJqFQG+4LZdtgUBsDIEd2H5qcb0cjJIuSYWnVCR+xxORilea5MWSH8563cctgx73rHCG7+2HsV+CIxSeVyx75p47jtE1MZN3YkXeRHjrd9YuqwgmTBbGPkiO7DymqdLEKOqVan3HsjjupiRHfXYWWH/yv59rw2yLyVW5i4cM1hf6bd8m+sfmxX3M0TybxUjtghH+7NHBEWPque1Sghx4TUqbUq5kNnnsT3/3PXEReM06YQ+IUR/vHHdLP048M7IYvIkVJ58bQTlS6LrHdVzJ7XBpt2AbZVjuqCz5w3gSV9U+NuikjiZe7iaSdqxm8ohZPDrr0HKq6GiTP8D+Zg+eadLN+8E9BcvUijFOwdJPTkUO4EEEfgF0/daNpGJJyCXY5Q6wRQHPztCvz9bw4dCvnLp2u6RqQaBbvUrVLwtyvwC9M1GsWLlKdgl6apFvhffvAp9h4YbOrX0yhepDwFu7RcaeC3Iug1ihd5m4Jd2q5c0C/6weMcGDw47M8ujOK/9L3fcMefn6OAl46UyjtPJVv6po3jt7d+lB39vdx12bmMHTli2J/51sEc81ZuYcrfrtXdrtJxNGKXRCkdzd+4+olD69obUQh4jeClk2jELom2pG8qO/p7uXz6hGF9TiHgL/zqz5rTMJEEU7BLKhQCfkd/LzN7Tmj4cwZ272fiwjXcuPqJ2pVFUkrBLqmzYu6MYY/il2/eyeRFazT/LpmkYJfUKozi77rs3IZ+kA/m0PSMZJKCXVKvb9o4tg0j4Ad272fSQo3eJTsU7JIZxQFfrxz50fuc+zc1v2EibaZgl8zpmzau4Tn4jc+8zNk3P9SCVom0j4JdMqswB/+ek4+v67hX3xhioqZmJMUU7JJ56+df0ND0jC6sSloFB7uZzTWzATM7YGabzGxGjfqzzewRM9sfHfcFM0vbc5klIwrTM/WugR/YvV9TM5I6QcFuZlcC9wLLgUuBvcA6M5tUof4M4EfAk8CfAvcDXwXmNaHNIg1bMXdG3aP3V98Y4swb1raoRSLNVzPYo1H2LcB97r7Y3dcClwAvAddWOOxa4CngL939YXf/CvmTwtXNabZI4wqj93rm3l8fymneXVIjZMQ+BTgdeLBQ4O6DwBrgogrHXAd82t2LH6LzJnBsg+0UabpG5t7nrdyi7Qgk8UKC/YzodWtJ+Tagx8y6Sw9w92fd/bcAZjY2msopTOeIJEZh9H5cd/jln+WbdyrcJdFCgn1M9LqvpHxfdHzF32fN7HRgD/CP5Ofbv9VAG0Va7umlF3PK6GOC6y/fvFM3M0lihQR7YShT+mziQnm1x968CnwY+AzwLmCTmb2jrhaKtMmvbriwrlUzG595WcshJZFCgv2V6HV0Sfko8qG+v9KB7r7H3X/q7t8FPk5+WufSRhoq0g4r5s6o647Vgd37Fe6SOCHBPhC9Ti4pnwx4yQVSAMysz8zeV1L8JDAI6BE2kmhL+qbWdVFV4S5JExrszwJ9hQIzGwH0AhsqHLMQuKOk7EPACEBXnSTxChdVjw68pqpwlySpGezRiLwf+JyZLTWzi4EHgBOBOwHMrMfMphcdthT4oJl928xmmdnngZXAzwDd6SGpsfU2hbukT9Cdp+5+D7AAuAJYBYwFZrv7tqjKTcCmovo/JH/H6R8BP4ze/2egt9zUjUiS1RvuWi0jcevK5eLNWTObCGzfsGED48ePj7UtItWcecNaXh8K+/9y12Xn0jdNl5OkdZ577jlmzZoFMMnddxS/p90dRQI9vfTi4BuZ5q3c0uLWiFSmYBepw9NLLw6elpm4cE1rGyNSgYJdpE71zLkr3CUOCnaRBtQT7tryV9pNwS7SoK239QbVe30op5Uy0lYKdpFh2NEfFu4bn3lZe7lL2yjYRYYpNNy1UkbaRcEu0gShe8tMWaSLqdJ6CnaRJuibNi7oUXtv5dC2A9JyCnaRJlk//4KglTIDu/drvl1aSsEu0kShK2U03y6tpGAXabLQ+Xatb5dWUbCLNFnofPvrQzk9FFtaQsEu0gKh8+3LN+9sfWOk4yjYRVokdL5dUzLSbAp2kRYKmW/XlIw0m4JdpIX6po3jlNHH1KynKRlpJgW7SIv96oYLg+qdt3R9i1sinULBLtIGIVMyL+57UzcuSVMo2EXaIHRKRjcuSTMo2EXaJHRK5uybH2pxSyTrFOwibXT59Ak167z6xpCmZGRYFOwibbSkb2rQjUuakpHhULCLtFnojUuakpFGKdhFYqApGWklBbtIDDQlI62kYBeJSeiUzJz7N7W4JZI1CnaRGM3sOaFmnY3PvNyGlkiWKNhFYrRi7gwCZmS0A6TURcEuErPt/bWnZF4fyulCqgRTsIskQMgTl3QhVUIp2EUSYP38C4Lqad92CaFgF0mIkB0gtW+7hFCwiyRE37RxHNdd+1Kqlj9KLQp2kQR5eunFNeto+aPUcnRoRTObC1wPjAe2APPdveLQwcw+ACwFpgGvAQ8DC9z9xWG1WCTjThl9DC/ue7NqnRtXP8GSvqltapGkTdCI3cyuBO4FlgOXAnuBdWY2qUL9PwQ2APuATwNfAmZGx4xoQrtFMitk33bNtUs1NYPdzLqAW4D73H2xu68FLgFeAq6tcNg1wPPApe7+Y3f/DvAp4Bwg7GkDIh0sZPmj5tqlkpAR+xTgdODBQoG7DwJrgIsqHPMUsCyqd+iw6LXsKF9E3hay/FFz7VJJSLCfEb1uLSnfBvSYWXfpAe5+j7t/s6T4Y9Hr0/U1UaQzhYzata5dygkJ9jHR676S8n3R8TV/+szsD4A7gEeBn9TTQJFOFTJq11y7lBMS7IWFtbkK5QerHRyF+oboa33K3Us/R0QqOGX0MTXraNQupUKC/ZXodXRJ+Sjyob6/0oFmdhbwH+RH/Re6+zONNFKkU2mFjDQiJNgHotfJJeWTAa80Ajez84BfAEPA+e7+eMOtFOlgIaN2kWKhwf4s0FcoiNai95KfYjmCmU0Efgy8CHzA3QfK1ROR2kJG7dqvXYrVvPPU3XNm1g/cbWZ7gI3k16mfCNwJYGY9wEnuvjk67Gvkp1+uBiaYWfGTe3/n7s838XsQybxad6O+PqRLV/K2oDtP3f0eYAFwBbAKGAvMdvdtUZWbgE1waDR/MdANfCcqL/4zp4ntF+kIIaP285aub0NLJA2C94px92XAsgrvXQVcFf19ENC2ASJNdlx3V9WRea39ZaRzaHdHkZQI2flR2wwIKNhFMkXbDAgo2EVS5fLpE2rW0UOvRcEukiIhe7DrodeiYBdJmZDNwaSzKdhFUiZkczAtfexsCnaRFBpz7BG7ZR9GSx87m4JdJIUeX1zpGTdv066PnUvBLpJR2vWxcynYRVJqZs8JcTdBEkrBLpJSK+bOqFnn7JsfakNLJGkU7CIpVusi6qtvDLWpJZIkCnaRFNNFVClHwS6ScbqI2nkU7CIpF7J/jHQWBbtIyoXsH6OLqJ1FwS6SAbqIKsUU7CIZEHIRVQ/h6BwKdpEOoYdwdA4Fu0hG6CKqFCjYRTIi5CLqlEVr2tASiZuCXSRDal1EfSvXpoZIrBTsIhkSchH1wq/+rPUNkVgp2EU6zMDu/XE3QVpMwS6SMbqIKgp2kYwJuYg6caEuomaZgl0kg2pdRJVsU7CLZFDIRVQtfcwuBbtIh9LSx+xSsItk1I7+3pp1tPQxmxTsIh1MSx+zScEukmEze06oWUePzsseBbtIhq2YO6NmHT06L3sU7CIZp6WPnUfBLpJxIUsfdcNSthwdWtHM5gLXA+OBLcB8d6/5SBYzGw08CVzn7qsabaiIiIQJGrGb2ZXAvcBy4FJgL7DOzCbVOG408ACgzStEYqSlj52lZrCbWRdwC3Cfuy9297XAJcBLwLVVjvsT4NfAuU1qq4i0kJY+ZkfIiH0KcDrwYKHA3QeBNUC1ybvVwBM16ohIm2jXx84REuxnRK9bS8q3AT1mVumS+/nu/klgd6ONE5Hm0a6PnSMk2MdEr/tKyvdFxx9f7iB3f3IY7RKRFtDSx84QEuxd0WvplkGF8oPNa46ItJIendcZQoL9leh1dEn5KPKhrisuIhmii6jpFxLsA9Hr5JLyyYC7uzb/FEkRXUTNvtBgfxboKxSY2QigF9jQonaJSIvoImr21bzz1N1zZtYP3G1me4CNwDXAicCdAGbWA5zk7ptb2VgRaY6ju/SgjSwLuvPU3e8BFgBXAKuAscBsd98WVbkJqLm9gIgkw9bbat+Jet7S9W1oibRC8F4x7r4MWFbhvauAqyq8t4O3V9CISEq8uO/NuJsgDdLujiIdShdRs0vBLtKhdBE1uxTsIh3suG7NkmaRgl2kgz299OKadebcr3URaaNgF5GqNj7zctxNkDop2EU63F2X6ZEJWaNgF+lwfdPG1ayjjcHSRcEuIjVvNNHGYOmiYBcRtgc8E1XSQ8EuIkG0OiY9FOwiAtRe067VMemhYBcRIGxNu6SDgl1Egt24+om4myABFOwicsjRNZbHLN+8sz0NkWEJ3rZXRLJv62292virBSr16Y4WrUbSiF1E6qIHcNSn2omyVSdRBbuI1EUP4Eg+BbuIHKZV0wPSPgp2Eamb5uHDxNVPCnYRkRi1IvwV7CJyhJDpGG0x0DzNDnctdxSRhqRli4F2LzVMAo3YRaSsy6dPqFnn7JsfakNLGhfHUsMkULCLSFlL+qbWrPPqG0NtaEnrZDXcFewiMixJDcdJCW1XO2iOXUQq2tEftsVAoU6S5q1zTf68cv2QpO+3mEbsItI0SR29V1PPiaueY+vti2bunKlgF5Gq6h2VpjHcq2nX99PMnTMV7CLSdHGvlsnayaVemmMXkZpC59oL0rpaJisnBI3YRSRIlqdkVj+2K1XtrUXBLiLBkroKZLjmrdzS8LFJPCEo2EWkLjv6e4MDPo7Q65SvWY2CXUQakuRwj0Px9xn396xgF5GWizvo2mXiwjWJ+F4V7CLSsHrm3NsReEkI1SQIXu5oZnOB64HxwBZgvrtX3JDZzM4CvgacB7wMfBP4irs3+05fEYnRe04+noHd+4PqTly4JrMXYJMkaMRuZlcC9wLLgUuBvcA6M5tUof7JwMPkt2v4JHAfsBS4rgltFpEEWT//grrqN/PWeSmvZrCbWRdwC3Cfuy9297XAJcBLwLUVDrua/G8Dl7j7WndfAtwGLDKzEc1puogkRT2j8OWbd2rKpMVCRuxTgNOBBwsF7j4IrAEuqnDMR4AN7v5aUdlq4ATgfY01VUSSrJEbmBTwrRES7GdEr1tLyrcBPWbWXeGYcvWLP09EMqaR+fNmhXutz+mkuf2QYB8Tve4rKd8XHX98hWPK1S/+PBHJoDjDXfJCgr0rei1dzVIoP1jhmEqrX8rVF5EMSWq4d8qoPSTYX4leR5eUjyIf0uXWOb1Spv7oovdERKREs05uIcE+EL1OLimfDHiFdekDFeoDeHjzRCStkjRqf8/J5WaMsys02J8F+goF0ZLFXmBDhWM2AB8xs+Le7AP+l/zNTSLSAQobhnXVrnpIK8K9eK39cKdj0jCdU/POU3fPmVk/cLeZ7QE2AtcAJwJ3AphZD3CSu2+ODrsH+AKw1sxuB84BFgEL3f3N5n8bIpJk26MwTPtF0jSEOgTeeeru9wALgCuAVcBYYLa7F5Yw3gRsKqr/PPm17EdH9f8auMHd72he00UkbUKDccqiZJwACr9xlG5VnPSAD94rxt2XAcsqvHcVcFVJ2aPAzGG0TUQyKOQxe2/VuaNUmn4TqPb9N+uEoWeeiog0oN7nwJYe20ratldE2i7pUxmhyk3T1KrfDgp2EUmkNE2vJI2CXURi0a7R6+XTJ8T2teOiOXYRSazSUXu5QK41sl/SN7WpbUoDjdhFJDU0PRNGwS4isZnZc0Ldx7Qr3NM8XaNgF5HYrJg7o6HjCo/X0wi+PAW7iMSqnn1kCkIfr1ftN4I0j8hrUbCLSKy2tzBga/1GUCnc0x76WhUjIrEbzl2czfjazfqcpEwNacQuIonQ7FFyI1M8WaERu4gkRmm4D2cE3MopnqTTiF1EMiftc+TDpWAXkUzp9FAHBbuIJJhCujEKdhFJtDHHdgfX1YkgT8EuIon2+OKLguolPdTb2T4Fu4ikXlJCPSk3PGm5o4gkXjueE9osSWiPgl1EUiEJgZkWmooREckYBbuISMYo2EVEMkbBLiKSMUm4eNoN8MILL8TdDhGR1CjKzCPu4EpCsJ8GMGfOnLjbISKSRqcBzxQXJCHYHwHOB54HhmJui4hIWnSTD/VHSt/oyuVy7W+OiIi0jC6eiohkjIJdRCRjkjDH3hZmtgqYChyIiha7+7/G2KTEMbNlwCnufnncbUkCM+sC7gA+GhXd7e73xNikxDCz+cBfATngUeCz7v5GvK1KFjN7J7AR6HP3re382h0T7MAfA+e4+ytxNySJzGw2cCWwLu62JMgl5AcDZwHHAf9lZj9396fibVa8zOw84C+A9wOvAf8EfBG4Pc52JYmZzQS+DZwRx9fviKkYMzsNGAN8z8weN7Mvm1lHfO8hzOxk4FZgadxtSRJ3fwDodfeDwEnkB0L/F2+rEmEPcI2773f3HPA4MCHmNiXNZ6M/v4/ji3dKuJ0CPAzMAWYAHwTmxtqihIimG/4euA7YG3NzEsfdB82sH/gtsAHYGXOTYufu/+3uP4dDg6a/AR6It1XJ4u5XuvvGuL5+pqZizOxTwN0lxS+6+3uBTxbV+zr5+cFvt7F5sarUN8DfAVvc/Zdm1tP+lsWvxs8N7r7QzG4Ffgh8HuiIefZa/WJmpwM/Br7l7g+3u31xq9U/ceqIdexmNh040d1/FP37z4DPuPsn4m1Z/Mzs18BI8jeHnQCMAr7r7lfH2rAEMLOpQM7dn4z+fQ1wlrt/Lt6Wxc/MziV/ovuKu38j7vYklZntAD6ii6etcQxwl5n9Anid/KjrH+JtUjK4+/sLfzezq8j/EHZ8qEfOAj5vZrPIT1t+nPxvOB3NzE4CHiI/z74q7vbIkRIZ7GZ2CbDC3UeXlM8FrgfGA1uA+e6+qdbnufsvzOw+4Ffkb8P9PrCi6Q1vg2b3TZa04Ofmu9HI9Dfkf6NZ5e7faX7LW6sFPzPXkl+McKOZ3RiVrXH3G5rY7LbJ4v+pxE3FmNkHyI8GjnL3UUXlV5IfZd9Cfm+ELwAzyS9h3B5HW9tNfVOZ+qY89Ut1We2fxIzYzexY8lfXbwX2k58+KbzXRb6D73P3xVHZesDJjx6+2PYGt5H6pjL1TXnql+qy3j9JWu74UWARsAAovRgzBTgdeLBQ4O6DwBrgonY1MEbqm8rUN+WpX6rLdP8kKdgfASa5+9fJ36ZcrHD3VumV5W1Aj5kdsdF8xqhvKlPflKd+qS7T/ZOYqRh331Xl7THR676S8n3kT07HA6+2ol1JoL6pTH1Tnvqluqz3T5JG7NV0Ra+lZ9ZC+cE2tiVp1DeVqW/KU79Ul/r+SUuwFzbuGl1SPop8J+9vb3MSRX1TmfqmPPVLdanvn7QE+0D0OrmkfDLg0UZEnUp9U5n6pjz1S3Wp7580BfuzQF+hwMxGAL3kN2bqZOqbytQ35alfqkt9/yTm4mk17p6Ldti728z2kN+8/hrgRODOWBsXM/VNZeqb8tQv1WWhf9IyYid6cs0C4ApgFTAWmO3u22JtWAKobypT35Snfqku7f2TuC0FRERkeFIzYhcRkTAKdhGRjFGwi4hkjIJdRCRjFOwiIhmjYBcRyRgFu4hIxijYRUQyRsEuIpIx/w+a1MZTTo4TewAAAABJRU5ErkJggg==\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   = 17\n",
      "    # data points      = 20000\n",
      "    # variables        = 1\n",
      "    chi-square         = 13.5013288\n",
      "    reduced chi-square = 6.7510e-04\n",
      "    Akaike info crit   = -146011.989\n",
      "    Bayesian info crit = -146004.085\n",
      "[[Variables]]\n",
      "    kappa:  10 (fixed)\n",
      "    sigma:  125 (fixed)\n",
      "    ga:     25 (fixed)\n",
      "    eta:    0.22270477 +/- 0.00139195 (0.63%) (init = 0.25)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAECCAYAAADq7fyyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f3/8ddkB0JIKAiIhiXAsSgKBWSxfi1uIFgawaqICy7UryguCBaKPlQaBFuplZ8i6rf1W4taKlrALwgialVExQoiqIewI8giJCGELcv8/phJyDIzuUkmmSXv5+PBY+DecyefXON7Ts4991yX2+1GRESiR0yoCxARkeBSsIuIRBkFu4hIlFGwi4hEmbhQF2CMSQT6Aj8AxSEuR0QkUsQC7YA11toT5XeEPNjxhPpHoS5CRCRCXQh8XH5DOAT7DwCvvPIKbdu2DXUtIiIRYe/evYwePRq8GVpeOAR7MUDbtm0544wzQl2LiEikqTKErYunIiJRRsEuIhJlFOwiIlFGwS4iEmXC4eJprSxcu5s/LrfsyT3G6alNmDTYkNmrfajLEhEJuYgM9oVrdzPlza85Vui5GLw79xhT3vwaoE7hXpsPCyfHBGrjax+gDy0RqbWIDPY/LrdloV7qWGExf1xuax2AtfmwcHJMoDZAlX2TXv8KXFBY7K7Q/osdh3j/uwNlYT/orNZ+/92iSTwuF+QeLazwd31IiDQOrlCvx26M6QhsW7lypeN57J0mL8FX1S5g28xhtarjgpnvsTv3WJXt7VObsGryxbU+JlAbwOc+X1zg83uuqdL3iXW5KHa7/b6214eASFj7/vvvueSSSwA6WWu3l98XkT3201Ob+AzE071hWRt7/ASsv+1Oj6nN+/oSrI/f0vcp9n6g+3vdnXuM++avY8I/11HirvjBEuOCEjcKf5EwFZHBPmmwqTCEAdAkPrZsfLo2avNh4eSY6to47bGHSok3zd0+tpWG/33z11U4Jq1pPI/88mwFvkiIROR0x8xe7ZkxogftU5vgwtNznDGiR52CZNJgQ5P42ArbqvuwcHJMoDa+9sXHuIiPdVXYVvFf4S/naCH3zV9Hx8lLKvzpNe0dFq7dHeryRKJeRPbYwRPuwewRlr5XTWajODnGSZvqZsUMOqs1b/xnd5ULxpGmNPBLe/jNEmKZflXdPpBFpKqIvHjaGFWeFlnTWTE5RwuDdgG2vsS44Pp+6WRl9gh1KSJhL+ounjZGwfgNpfTDYXfuMb+zYUIZ/iVumPfpTuZ9uhPQWL1IbSnYGxGnHw6+PgBCEfjlh240bCPinIJdqqjuA6B88DdU4BecLC4L+Rv6a7hGJBAFu9SYv+BvqMAvHa5RL17ENwW7BE2gwH908UZyjxUG9eupFy/im4Jd6l3lwK+PoFcvXuQUBbs0OF9BP+XN9RwrLKnze5f24ie+/hVP/vo8Bbw0ShF556lEl8xe7fn291ewfeYw/nxtT1KbxNf5PYtK3Nw3fx1dfrdUd7tKo6Meu4SVyr35hxZ+XTavvTZKA149eGlM1GOXsJaV2YPtM4dxQ//0Or1PacBf9qcPglOYSBhTsEtEKA347TOHcUFGy1q/T/b+AjpOXsJDC7+uvrFIhFKwS8R5ZeyAOvfi5326k85Tlmj8XaKSgl0iVmkv/s/X9qzVD3KJGw3PSFRSsEvEy+zVnq11CPjs/QV0mqzeu0SPsFm299JLL6Vp06Zcf/319OrViy1btjB37twq7W+55Ra6d+/ON998w0svvQRAXFwcsbGxxMXFMWbMGDp27Mi3337L0qVLiY2NLdsXGxvLVVddRevWrdm8eTNr1qypcGxsbCwXXXQRycnJ7Nq1i23btlXYFxsbyznnnEN8fDwHDx4kNze3bH98fDyJiYm0aNEClyvSHo0RXRau3V3lqU5OXZDRklfGDghyRSLBFxHL9v7973+nuLiYPn360KtXL3bv3s2cOXOqtBs0aBDdu3dn+/btzJkzB7fbTXFxMUVFRZSUlDBo0CA6duzIl19+ycSJE6sc36dPH1q3bs27777LnXfeWWX/pk2b6Nq1K/Pnz2fSpElV9u/Zs4d27doxe/Zspk2bVmX/4cOHad68Ob/97W+ZM2cOiYmJZX+SkpLYuHEjLpeLWbNmsXz5cpKSkkhMTKRZs2akpaXx1FNPAbB8+XJ27NhBcnIyzZo1Izk5mdTUVHr37g1Afn4+cXFxJCUl6YOkktIpk7WZKrlqyyHOfWQZ6x8bUk/VidS/sOmxB+NBG6Xfi8vlorCwkGPHjlFcXFz2p6ioiNatW5OQkEBubi579+4t2176eu6555KUlMTOnTvJzs6ucGxxcTFDhgwhKSmJr776ivXr11NUVERRURGFhYWcOHGC8ePHExcXx+LFi/nggw84ceJE2Z+ioiJee+01AGbOnMmiRYs4ceIEx48f5+jRo8THx5OdnQ3AyJEjefPNNyt8f2eeeSY7d3qCasiQISxfvpy4uDjS0tJIS0ujZ8+ezJ8/H4Cnn36aAwcO0LJlS9LS0mjdujUdO3bknHPOqdM5jkSX/ekDsvcX1Pi4P1/bU/PeJWwF6rFHVbBHkyNHjpCXl0dBQQFHjhzhyJEjuN1uLrroIgDeeOMNsrOzycvLIycnh5ycHNq0acPs2bMBuOiii/j4448pKTl1m/5FF13EBx98AMCAAQPIycmhbdu2tGnThjZt2tCvXz9Gjx4NeH4zad26NfHxdb8LNBzUdnim62nNWDHhF8EvSKSOgjIUY4wZCzwInAGsAyZYa1cHaD8YyAK6A3uA2cAz1tpwfjpb2EhOTiY5Odnv/pEjRwY8/t///jclJSXk5+dz6NAhDhw4QGzsqQdnDxo0iM2bN7N3717WrVvH3r17ycnJYfTo0bjdbowxHD16lNNPP50OHTqQnp5OZmYm11xzDeAZckpJSQnON9sASodnRr+4mlVbDjk+Lnt/gYZmJOI4CnZjzE3AXGAasAYYDyw3xpxnrd3mo/0A4P+AecAU4GfAn7xf76nglC7ViYmJoUWLFrRo0YJOnTpV2Pf4449XaV9c7HlYdklJCU899RQ7duxg586d7Nixg9WrV9O1a1cA8vLySE1N5bTTTqNbt24YYzj77LO54oorOOuss+r/G6uDV8YOqHHv/fCJYs6aupTvpg+tx8pEgqfaYDfGuPAE+gvW2se821YAFrgfuMfHYfcDG4FbvT30d40xPwXuQsEetkp79LGxsdx+++1V9pcO27ndbp544gk2bdrEpk2beOutt/jLX/5C06ZNOeuss9i8eTMTJkzgvPPOo2fPnpx//vmceeaZDfq9BFLae6/J2PvxYjcdJy/RuLtEBCc99i5AB2Bx6QZrbaExZgng7/fTB4DkSsMuJ4HE2hYqoVc6+yY1NZUHH3ywwr59+/aRmOj5z3vgwAG2bt3K0qVLy34LOOOMM1iwYAH9+vWjoKCA+Ph4EhISGvYbqGTFhF/UuPd+3/x1fLHjkB7qIWHNyf0c3byvmytt3wpkGGNiK23HWrvLWvstgDEm1TuUUzqcI1GoTZs2pKamAp4Lsxs2bCA/P5/PP/+c2bNnc+GFF5Ke7lkC4IUXXiA1NZUhQ4Ywa9YsvvrqqwoXeRtSZq/2bJ85jKRY51NG5326U2vNSFhzEuylV8jyK23P9x7fzN+BxpgOQA7wN2AD8FwtapQI1aRJE/r27cv48eN59dVXadeuHQD9+/fn9ttvZ+fOnUycOJGePXty5plncuLECYCQhPx304fSprnz3yDmfbqT0S/6nTsgElJOgr20K1N5Nkvp9kD/Fx4GLgauB9KA1caYpjWqUKLOgAEDmD17Nt988w27du3if//3fxk7dmzZUM7QoUO5/PLLee6559izZ0+D1fXZ1MtqtHLkqi2HtM6MhCUnwZ7nfW1eaXsynlD3e/XJWptjrX3fWvsacBWeYZ3A8/SkUTnjjDO4+eabefTRRwHPhdk+ffqwfft2xo0bR/v27bngggv417/+1SD1vDJ2QI1WjczeX6Bwl7DjJNizva+dK23vDFhf89KNMZnGmL6VNm8ACgFNKRC/XC4XWVlZWGvZuHEjWVlZ5Obmsm2bZ1btkSNHWLRoEYWFwXsQdmVZmT3487U9HbdXuEu4cRrsu4DM0g3GmHhgGLDSzzGTgScrbRsExAO66iTVcrlcdO/enalTp7Jhwwbuucczq3bRokVkZmbSsWNHpk2bxg8//FAvX7/0omqcw2uqCncJJ9UGu7dHPhP4b2PMdGPMUGAR0ArvnHRjTIYxpn+5w6YD/2WMed4Yc4kx5k5gPvABsDTI34NEOZfLRVycZ2butddey6JFi+jRowePPPII6enpXHPNNeTnV762HxybZyjcJfI4Wr7aWjsHmATcCCwAUoHB1tqt3iYPA6vLtX8L+BWeO07f8u7/OzBMSwpIXcTFxTF8+HCWLVtGdnY29957L/n5+WXLL3z33XcEe/2jmoa7ZstIqGkRMIkaP/74Ix06dKBr16787ne/4+qrryYmJnjPkjlr6lKOFzv7/0V3qEp9C7QImJ6gJFEjJSWFZ599lhMnTnDttdfSp08fli9fHrQe/HfThzq+kam2D/oQCQYFu0SNhIQExowZw4YNG3j55ZfJyclhyJAhZWvcB8N304c6HpbpOHlJ0L6uSE0o2CXqxMbGcuONN2Kt5e2336ZbN8+qGH/729/Yt29fnd+/JmPuCncJBQW7RK2EhASGDPGsU/fDDz9wxx130K1bN55++uk6z4OvSbifNVUTwaRhKdilUWjXrh1fffUV/fv357777qNPnz58+eWXdXrPzTOGOWp3vNitmTLSoBTs0mgYY1i2bBlvvPEGBw4cYNCgQeTl5VV/YADbZzoL91VbDrFw7e46fS0RpxTs0qi4XC5GjBjBxo0bef3112nRogVut5utW7dWf7AfTsNdM2WkoSjYpVFKS0vj8ssvB+Af//gHZ511Fn/6059qvWSw07VlukzRxVSpfwp2afQGDx7MsGHDeOCBB7jyyivZv39/jd8js1d7up7m99EEZYrcaNkBqXcKdmn0WrZsyZtvvsmzzz7Le++9x3nnncdHH31U4/dZMeEXjmbKZO8v0Hi71CsFuwiesfdx48axZs0aWrRowYEDB2r1Pk5nymi8XeqTgl2knB49erB+/XpGjBgBwPvvv1/2yD6nnI63a3671BcFu0glCQmeZ5/u3LmTwYMHc/HFF/Pjjz86Pt7pePvxYrceii31QsEu4kd6ejqvvvoqX375JQMHDmTLli2Oj3U63j7v0511qFDENwW7SABXX301K1eu5ODBgwwYMIDPP//c8bFOx9s1JCPBpmAXqcbAgQP55JNPSE5OZtmyZTU61sl4u4ZkJNgU7CIOGGP48ssvefjhhwEoKChwdFxmr/a0aZ5QbTsNyUgwKdhFHEpNTcXlcpGdnU23bt34xz/+4ei4z6Ze5qhdv+kr6lKeSBkFu0gNtW3bli5dunD99dfz0ksvOTrGyZDMvvyTunFJgkLBLlJDzZs35+233+bSSy/l1ltv5eWXX672GKdDMrpxSYJBwS5SC02bNmXx4sVcfPHF3HLLLaxcubLaY5wOyZz7SM0u0IpUpmAXqaWkpCQWLVrExIkTGTBggKNjbuifXm2bwyeKNSQjdaJgF6mD5ORknnjiCZo2bcrhw4erneeeldnD0Y1LGpKRulCwiwTJXXfdxWWXXcbXXweek+70xiUNyUhtKdhFgmT69Ok0a9aMoUOHsnt34KEUDclIfVKwiwRJeno6S5cuJTc3l6FDh3L48GG/bTUkI/VJwS4SRD179mTBggVs3LiRu+66K2Bbp0Myo19cHYzSpBGJC3UBItFm8ODBzJ8/nz59+lTb9oKMlqzacihgm+r2i1SmHrtIPRg5ciQdOnSgpKSEL774wm+7V8YOwMGIjFaAlBpRsIvUoxkzZjBw4EBWrVrlt822mdUPyRwvdutCqjimYBepR+PGjaNDhw6MHDmS77//3m87J09c0oVUcUrBLlKP0tLSWLRoEUePHuWqq67i+PHjPtutmPALR++nddvFCQW7SD3r3r078+bN44svvmDixIl+2zlZAVLrtosTCnaRBjB8+HBmzJjBdddd57dNZq/2JMVWfylV0x+lOgp2kQYyefJkfv7znwP4HZL5bvrQat9H0x+lOo6D3Rgz1hiTbYw5ZoxZbYwJuJydMWagMeZ9Y0yuMWaPMeZlY0ybupcsEtkef/xxBg4c6DfcnazbrrF2CcRRsBtjbgLmAvOAkUAusNwY08lP+58CK4F8YBQwEbjAe0x8EOoWiVjnnXcea9eu5f777/e538m67Rprl0CqDXZjjAuYBrxgrX3MWrsUGA78CPj+yYS7gR+Akdbat621rwLXAecBzp42IBKlhg0bxsSJE5k7dy5vvfWWzzZOpj9qrF38cdJj7wJ0ABaXbrDWFgJLgCF+jtkIzPK2KzvM++qzly/SmGRlZdGzZ09uu+029u3bV2W/k+mPGmsXf5wEezfv6+ZK27cCGcaY2MoHWGvnWGufrbT5l97X72pWokj0SUxMZN68eRQWFrJmzRqfbZz02jXWLr44CfYU72t+pe353uOr/ekzxpwJPAl8AbxXkwJFotXZZ5/Njh07uPLKK33ud9Jr11i7+OIk2Esn1rr9bC8JdLA31Fd6v9Z11trK7yPSaKWkePpN8+fPZ9OmTVX2a4aM1IaTYM/zvjavtD0ZT6gX+DvQGHMO8AmeXv9l1tottSlSJJodOnSIO++8k5tvvpni4uIK+zRDRmrDSbBne187V9reGbD+euDGmH7Ah0AxcKG1dn2tqxSJYi1btuSZZ57h008/Zfbs2VX2O+m1i5TnNNh3AZmlG7xz0YfhGWKpwhjTEXgb2AcMtNZm+2onIh6jRo3il7/8JVOnTmXz5orzFJz02rVeu5RX7ROUrLVuY8xM4BljTA6wCs889VbAUwDGmAygtbX2U+9hT+MZfrkLSDfGlH9y7w5r7Q9B/B5EIp7L5eK5557j7LPP5rbbbuP9998nJuZUv6tN8wT25Z/0e/zxYl26klMc3XlqrZ0DTAJuBBYAqcBga+1Wb5OHgdVQ1psfCsQCr3q3l/8zOoj1i0SN9u3bM3fuXO6++25croqLgTnptfebvqK+SpMI4/iZp9baWcAsP/vGAGO8fy8EtGyASC2UX/3R7XZXCPikWFfAnnmgHr00LlrdUSQMPfvss4wYMQK3+1SQO1n5UcsMCCjYRcJSSUkJCxcu5J///GeNjtMyAwIKdpGwNG7cOHr37s19991HXl5e2fYb+qcHOMpDD70WBbtIGIqNjeX5559n//79PPTQQ2XbszJ7VHusHnotCnaRMNW7d2/uuusunnvuOXbuPHV3qZPFwaRxU7CLhLGsrCw+/vhj0tNPDcE4WRxMUx8bNwW7SBhLSUmhf//+ABw8ePDU9sQqq2VXoKmPjZuCXSQCzJ07l4yMDPbs2QPA+sf8PePmFK362Hgp2EUiwGWXXcbx48f57W9/6/gYrfrYeCnYRSJARkYGEydOZN68eaxatQqACzJahrgqCVcKdpEIMWXKFNq3b88999xDcXExr4wdUO0x5z6yrAEqk3CjYBeJEM2aNePJJ5/k22+/Zf16z+MNqruIevhEccD9Ep0U7CIR5Nprr2XLli306tUL0EVU8U3BLhJBXC4X7dq1w+12l/Xaq6OLqI2Pgl0kAs2ePZuf/exnbNiwwdH6MdK4KNhFItANN9xASkoK9957L7//1TnVttdF1MZFwS4SgX7yk5+QlZXFe++9x8KFC3URVSpQsItEqN/85jd0796dBx98kC+mXlxtez2Eo/FQsItEqLi4OJ588kkKCgrYvHlzte31EI7GQ8EuEsGGDBnCli1b6N69uy6iShkFu0gEc7lcNGnShJMnT3JZWvU98i5TljRAVRJqCnaRKPDoo49y6aWXknjkh4DtitwBd0uUULCLRIF77rmHpKQkuu+uvkd+2Z8+qP+CJKQU7CJRoG3btkyZMoWFCxdyfGfgJQSy9xc0UFUSKgp2kShx//33c+aZZ+L67GXc7pJQlyMhpGAXiRJNmjRhxowZnJ7ahJKCvIBtO07WRdRopmAXiSKjRo1izZo1pP2kVahLkRBSsItEkZiYGGJiYvjg3vM5tvU/Adtq6mP0UrCLRKEJEyZwYOHjFOUf9NtGUx+jl4JdJAo9/PDDxFFC7kfzArbT1MfopGAXiUIZGRmMHz+egq/f5eT+rX7baepjdFKwi0Sphx56iIRmzcl576+43f7HXfTovOijYBeJUmlpafzx8d/jSmyCu+iE33Z6dF70iQt1ASJSf8aPH8//HOymB200Muqxi0Qxl8vF+seGUHjw+4DTH3XDUnRx3GM3xowFHgTOANYBE6y11T6SxRjTHNgAPGCtXVDbQkWk9g6teI6TB7bT/jcvEJPYLNTlSD1z1GM3xtwEzAXmASOBXGC5MaZTNcc1BxYBegKASAh9/MZfKDmaR97q1/220dTH6FFtsBtjXMA04AVr7WPW2qXAcOBH4P4Ax10EfA70DFKtIlJLvXv3ptnZgzj8xSKK8vb5bKOpj9HDSY+9C9ABWFy6wVpbCCwBhgQ4biHwdTVtRKSB3HrfFFyuGHL+/XKoS5F65iTYu3lfKz8tdyuQYYyJ9XPchdbaa4D9tS1ORIJn9u2XkTLg18Qlt/S7rK8uokYHJxdPU7yv+ZW25+P5YGgGHK58kLV2Q91KE5FgSx80WlMfGwEnPXaX97XyrWul27Wiv0iEWP+YZ2T02PZ1HNu+zmcbXUSNfE6CvXTF/uaVtifjCXVdcRGJIG53CTnv/Q+Hlj+Du6iwyn5dRI18ToI92/vaudL2zoC11mrxT5EIcuOAjqQNupWi3L3kf/l/oS5H6oHTYN8FZJZuMMbEA8OAlfVUl4jUk6zMHjTp9DOSOvUm75N/UHysyiUyXUSNcNUGu7dHPhP4b2PMdGPMUDw3HbUCngIwxmQYY/rXa6UiEjRxLkgbdCslJ4+Rt+q1UJcjQebozlNr7RxgEnAjsABIBQZba0sXen4YqHZ5AREJD5tnDCOhdQdS+o0k/idn+mzTb/qKBq5KgsXxWjHW2lnALD/7xgBj/OzbzqkZNCISRtIuutnvvn35JxuwEgkmre4o0kjd0N+zhJO7pJj8r5ZzfJduPYkWCnaRRiorswcA7uIi8j6Zz6F3nsNdXFShjS6iRiYFu0gjlhTrIiY+kZaXjKXwxx2a/hglFOwijdh304cC0KRrf5I69Sb341coOnKoQpvRL2peRKRRsIsILpeLlpf+BndxIbnv/7XCvlVbDvk5SsKVgl2kkfvztZ5HJsS3bE/aoNto1v0XoS1I6kzBLtLIZfZqX/b3lN6/pElGnypttDBYZFGwi0iFG03cJcXkfPgy+WuXlm3TwmCRRcEuImybOezUP1wxnNy7hZwPXqLo8IHQFSW1pmAXkQpcLhctL78T3CUcemcObrdnAVfNjokcCnYRATxz2kvFp7Yl9ec3cGzLGo5+9xGg2TGRRMEuIsCpOe2lmvcZTkK7ruS89xfcxVUfyCHhS8EuIj65YmL5ydD7Oe2ax3DFxgPw0MKvQ1yVOKFgF5EycZXWYU1olU5C644AFB/NY96nOxu+KKkxx8v2ikj02zxjmM+Fv3I/foUjX71Du1ufCUFVkc/fYmrby89GCiL12EWkWk269KP4aC6H3pnD+VnvhLqciBJohcz6Wj1TwS4i1Ups24XUn4/m6HcfsfUzBXu4U7CLSAX+hgdS+o0ksX13Dq2Yy86dGmsPZwp2EXHEFRPLT66cQExCE/pOejnU5USEUD2oRMEuIo7Fp7al/R0v0qRTr1CXEjXqI/wV7CJSRaDZGq7YeNxuN31GPcDixYsbsKroFexw13RHEam5kiK+/mAxNy/7K//5z3/o3LlzqCvyq6GnGoYD9dhFxKcb+qf73eeKjafVryZz+HgRI0aM4MiRIw1YmXOhmGoYDhTsIuJTVmaPgPvjU9vS6peT2LBhA6NGjaK4uLiBKgueaA13BbuI1FqTzr1pcckdLFmyhA8//DDU5VTQKUpD2wmNsYuIX9tn+l5ioLzmvYaS2P6n3LL8KNsHNVBhDriD/H6+zkO4jtOrxy4idZZwWicA2lwzjb/97W8hrqZmnAzH+GsTzDH8YK6cqWAXkYBq0ivNX7uEMbfcyvz58+uxoobVUOPwwVw5U8EuIkHTaviDJLb/KdddP5rXXnstZHVE60VRpzTGLiLVcjLWDhATn8RpVz/C/jemMXr0aHJychg3blwDVBgc0fKBoB67iDjidEgmJrEpba6ZRlLG+axfv77sYdjhbOHa3VET6qBgF5EacBrurrgEWl/1O5599llcLhebNm3i8OHD9Vxd7d03f12tjw3HDwQFu4jUyPaZwxwFvCsmloypyygsLGTYsGH079+f9evX13t9oQjacAt3BbuI1IrT3nvXh9/hxRdf5NChQ/Tp04cnnngiIu9SrU75cA910CvYRaTejVlWwIYNGxg+fDiTJ0+mf//+HDhwINRlBV3HyUtCHuqgYBeROqjJHPc+T37G66+/zquvvkqHDh1o1aoVALm5uUGrJxxCNRw4nu5ojBkLPAicAawDJlhrVwdofw7wNNAPOAQ8C/zBWhv+l8hFxLGupzUje3+Bo7adpixl+8xRjBo1CoCDBw/SuXNnrrjiCiZMmEDfvn1xuVz1WW6j4KjHboy5CZgLzANGArnAcmNMJz/tTwPexbNcwzXAC8B04IEg1CwiYWTFhF/UqH35W+djYmIYN24cS5cupV+/fvTo0YNZs2Zx8ODBIFfZuFQb7MYYFzANeMFa+5i1dikwHPgRuN/PYXfh+W1guLV2qbU2C5gBTDHGxAendBEJFzUZkpn36c6yIZO0tDRmzJjB999/zwsvvEDz5s2ZOHEihw4dAuCTTz5hwYIF7Nq1KyLmw4cLJz32LkAHoOwZWNbaQmAJMMTPMZcCK621R8ttWwi0BPrWrlQRCWc1Xemw/IXGlJQUxo4dy+rVq9m6dStdu3YF4MUXX+TXv/416enptG3blgsvvJA77rijLOS//kytm0wAAAjpSURBVPpr1q5dy+7duzl8+DDukuibbVMbTsbYu3lfN1favhXIMMbEWmsrn81uwAc+2pfu+6QmRYpIZHC69EB5HScvqfCh0KnTqRHe559/nnHjxvHZZ5+xdu1asrOz+eabb8rG4R944AFWrFhR4f0S2mTQbszTABz41+OcPLAdXDGYdi3YtL+AhLZdaDX0Ps/+t/5Icd7+isefbmh58e0A7H8zi5KjeRX2J6WfS+p/3QjAvn8+gvvksQr7m2T0ocWAawDY++pkcJdU2N/UXEBKn1/hLirkyNcraN5rqPOT5ZCTYE/xvuZX2p6Pp8ffDKh8S1mKn/bl309EolAwwr1UQkICffv2pW9f37/o/+EPf2Dbtm3s37+fKfM/p6TwBDFJyWX741ulQ2wclJRgzGlsPbmH2KapZftdsQm44hIrvKcr9tRosSsuscp+YuPK7Y+HkqKK+2PK7Y+JAyoGO64Y338PIifBXnqJuvIAV+n2SlWX7fM3IOarvYhEkWCGeyA9e/akZ8+eAMzYdkaV/akX3lD29ze8712+rlZD7w34/q2HTwq4/7QRDwXc3+a6LL/7XHHxNO/pbzS7bpx8XJT+HtK80vZkPCHta55Tno/2zcvtExGRSoI1D99JsGd7XztX2t4ZsH7mpWf7aQ9gnZcnIpGqNo+Nq68bjLqe1qxe3jdcOQ32XUBm6QbvlMVhwEo/x6wELjXGlD+bmcBBPDc3iUgjULpgWE1uOaqPcC8/176uzykN1+ecllftGLu11m2MmQk8Y4zJAVYBdwOtgKcAjDEZQGtr7afew+YA44Glxpg/AucBU4DJ1tqTwf82RCScbfMxvh2JIiHUweGdp9baOcAk4EZgAZAKDLbWlk5hfBhYXa79D3jmssd52/8GmGqtfTJ4pYtIpHEajF2mhMcHQOlvHJWXKg73gHe8Voy1dhYwy8++McCYStu+AC6oQ20iEoWczJgpquFNppH0m0Cg7z9YHxh65qmISC3UZkpn+WPrk5btFZEGF+5DGU75Gqaprn1DULCLSFiKpOGVcKNgF5GQaKje6w3900P2tUNFY+wiErYq99p9BXJ1PfuszB5BrSkSqMcuIhFDwzPOKNhFJGQuyGhZ42MaKtwjebhGwS4iIfPK2AG1Oq708XrqwfumYBeRkKrNo6vLP14vkEC/EURyj7w6CnYRCalt9Riw1f1G4C/cIz30NStGREKuLndxBuNrB+t9wmVoSD12EQkLwe4l12aIJ1qoxy4iYaNyuNelB1yfQzzhTj12EYk6kT5GXlcKdhGJKo091EHBLiJhTCFdOwp2EQlrKYmxjtvqg8BDwS4iYW39Y0MctQv3UG/I+hTsIhLxwiXUw+WGJ013FJGw1xDPCQ2WcKhHwS4iESEcAjNSaChGRCTKKNhFRKKMgl1EJMoo2EVEokw4XDyNBdi7d2+o6xARiRjlMrPKHVzhEOztAEaPHh3qOkREIlE7YEv5DeEQ7GuAC4EfgOIQ1yIiEili8YT6mso7XG63u+HLERGReqOLpyIiUUbBLiISZcJhjL1BGGMWAD2AY95Nj1lr/xXCksKOMWYW0MZae0OoawkHxhgX8CRwhXfTM9baOSEsKWwYYyYAtwFu4AvgDmvtidBWFV6MMS2AVUCmtXZzQ37tRhPsQB/gPGttXqgLCUfGmMHATcDyUNcSRobj6QycAyQB3xhj/m2t3RjaskLLGNMPuAU4HzgKvAzcA/wxlHWFE2PMBcDzQLdQfP1GMRRjjGkHpACvG2PWG2MeNcY0iu/dCWPMacDvgemhriWcWGsXAcOstSVAazwdoSOhrSos5AB3W2sLrLVuYD2QHuKaws0d3j97QvHFG0u4tQHeBUYDA4D/AsaGtKIw4R1u+CvwAJAb4nLCjrW20BgzE/gWWAnsDHFJIWet3WSt/TeUdZruBRaFtqrwYq29yVq7KlRfP6qGYowx1wHPVNq8z1p7NnBNuXaz8YwPPt+A5YWUv3MD/AVYZ639yBiT0fCVhV41PzdYaycbY34PvAXcCTSKcfbqzosxpgPwNvCctfbdhq4v1Ko7P6HUKOaxG2P6A62stf/n/ffVwPXW2hGhrSz0jDGfA03w3BzWEkgGXrPW3hXSwsKAMaYH4LbWbvD++27gHGvtf4e2stAzxvTE80H3B2vt/wt1PeHKGLMduFQXT+tHAvBnY8yHwHE8va6XQltSeLDWnl/6d2PMGDw/hI0+1L3OAe40xlyCZ9jyKjy/4TRqxpjWwDI84+wLQl2PVBWWwW6MGQ68Yq1tXmn7WOBB4AxgHTDBWru6uvez1n5ojHkB+AzPbbhvAK8EvfAGEOxzE03q4efmNW/P9Cs8v9EssNa+GvzK61c9/Mzcj2cywkPGmIe825ZYa6cGsewGE43/T4XdUIwxZiCe3kCMtTa53Pab8PSyp+FZG2E8cAGeKYzbQlFrQ9O58U/nxjedl8Ci9fyETY/dGJOI5+r674ECPMMnpftceE7wC9bax7zbVgAWT+/hngYvuAHp3Pinc+Obzktg0X5+wmm64xXAFGASUPliTBegA7C4dIO1thBYAgxpqAJDSOfGP50b33ReAovq8xNOwb4G6GStnY3nNuXySu/eqnxleSuQYYypstB8lNG58U/nxjedl8Ci+vyEzVCMtXZ3gN0p3tf8Stvz8Xw4NQMO10dd4UDnxj+dG990XgKL9vMTTj32QFze18qfrKXbSxqwlnCjc+Ofzo1vOi+BRfz5iZRgL124q3ml7cl4TnJBw5YTVnRu/NO58U3nJbCIPz+REuzZ3tfOlbZ3Bqx3IaLGSufGP50b33ReAov48xNJwb4LyCzdYIyJB4bhWZipMdO58U/nxjedl8Ai/vyEzcXTQKy1bu8Ke88YY3LwLF5/N9AKeCqkxYWYzo1/Oje+6bwEFg3nJ1J67HifXDMJuBFYAKQCg621W0NaWBjQufFP58Y3nZfAIv38hN2SAiIiUjcR02MXERFnFOwiIlFGwS4iEmUU7CIiUUbBLiISZRTsIiJRRsEuIhJlFOwiIlFGwS4iEmX+PxsBoWF8t+ecAAAAAElFTkSuQmCC\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": 21,
   "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": 22,
   "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": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYhklEQVR4nO3deZhlVXnv8W9JKyE2OACKgMoUXsZmCEQjwQHwEbmCjZCBi8ZEoyGPElRABiEyiHgdoqCoYOKdxAjxAqJyEUIgyJSLSQBRfAGBMEgzzyJDU/ePtQ/sPpzqqjq1Tp061d/P89Rzuve4Vu/e51dr7732GhsfH0eSpBpeMOwCSJLmD0NFklSNoSJJqsZQkSRVY6hIkqpZMOwCDEJErAxsD9wJLB1ycSRpVKwEvAq4MjOf6GcD8zJUKIHy42EXQpJG1I7AJf2sOF9D5U6AU089lbXWWmvYZZGkkbBkyRL23XdfaL5D+zFfQ2UpwFprrcW666477LJI0qjp+7aBN+olSdUYKpKkagwVSVI1hookqRpDRZJUjaEiSarGUJEkVTNf+6lonhgbGxvIdh2cThoMWyqSpGoMFUlSNYaKJKkaQ0WSVI2hIkmqxlCRJFVjqEiSqjFUJEnVGCqSpGoMFUlSNYaKJKkaQ0WSVI2hIkmqxlCRJFVjqEiSqjFUJEnVGCqSpGoMFUlSNYaKJKkaQ0WSVI2hIkmqxlCRJFVjqEiSqlkw7AJIwzA2NlZ9m+Pj49W3KY0aWyqSpGoMFUlSNXP28ldE7AGcmpmrtqZtB1zZY/EvZOZBs1Y4SVJPczJUIuINwLeA7gvfi4DHgF26pv9qNsolSVq+ORUqEbEycABwLCU8XtS1yCLg2sy8YrbLJkma3Fy7p/J24DDgYODLPeYvAq6Z1RJJkqZsroXKlcD6mXki0Ov5zC2BV0fEVRHxZETcGBHvnd0iSpImMqcuf2XmHRPNi4i1gTWA36G0Zh4A9gH+R0SMZ+b/mp1SSpImMqdCZRIPArsC12Tmnc20f2rC5pOAoSJJQzYyoZKZvwZ+1GPWucCuEbEwMx+d5WJJklpGJlQiYmNgJ+C/Z+YTrVmrAI9TnhaTJA3RXLtRvzzrAF8DdutMiIgx4F3AjzPTFy9J0pCNTEsFuBi4BPh6RLwMuBP4S8pjxn8wzIJJ4EsqJRihlkpmLgXeCZwJHAOcAawJvDUzfzLMskmSijnbUsnMo4CjuqbdD+w3jPJIkiY3Mi0VSdLcZ6hIkqoxVCRJ1RgqkqRqDBVJUjWGiiSpGkNFklSNoSJJqsZQkSRVM2d71EsazPvEwHeKaXBsqUiSqjFUJEnVGCqSpGoMFUlSNYaKJKkaQ0WSVI2hIkmqxlCRJFVjqEiSqjFUJEnVGCqSpGoMFUlSNYaKJKkaQ0WSVI2hIkmqxlCRJFVjqEiSqjFUJEnVGCqSpGoMFUlSNYaKJKkaQ0WSVI2hIkmqZsGwCyBp9o2NjVXf5vj4ePVtavTYUpEkVWOoSJKqMVQkSdUYKpKkagwVSVI1hookqRpDRZJUjaEiSarGUJEkVWOoSJKq6fs1LRHxWmB9YA1gKXAXcGtm3l6pbJKkETOtUImILYAPA28H1m0md14iNN4scyNwDvD3mXltvwWLiD2AUzNz1da0MeBw4C8pYXYpsH9m/qLf/UiS6plSqDRh8iVgJ+B+4FzgGuAm4GHKZbTVKUHzemAf4K8j4nzg8Mz89+kUKiLeAHyL5wKr42+AQ4FDgFuAI4ALImKzzHxoOvuQJNU3aahExOeB/YB/AHYGLs7MpZOsMwbsAvwpcElEfC0zD5zCvlYGDgCOBR4DXtSatypwEHBUZp7YTPsx8J/A+4G/nWz7kqTBmsqN+pcAm2bmBzLzQmCjyVbIzPHMPD8z3wNsDrxsiuV5O3AYcDDw5a55rwcWAme39vMA8C/ArlPcviRpgCYNlSZMbmtN+mlEfCEiVpvKDjLz5sx83xTLcyWwftMS6R6cYePm85dd029qzZMkDVE/jxSfDOwP3BARf1GzMJl5R2Y+OMHs1YAnMvPJrumPNPMkSUM27VDJzP2BbSg36k+JiJ9ExA7VS/Z8Yzy/9dKZ/sws7F+SNIm+Oj9m5s8y863AXpT7JRdHxLcjYt1JVp2Jh4CVI+KFXdMXNvMkSUM2ox71mXkmsClwJLA78IuIOKJ5iqu2GyitkvW7pm8A5AD2J0maphm/piUzn8zMTwMB/CtwNHBdRBwaEW+LiFfMdB+Ny4DfAIs7EyLiZcCbgAsq7UOSNAN9vaYlIjYCtgI2af0E5VIUwHrAp3mul/2SzFxnJgXNzEcj4svApyLiGeB64BOUzpd/N5NtS5LqmHaoRMQJlFe1dHq730P5gj+dchmq8/MUsDXlpv7WNQpLeUXLM5ROkAsprZf32ptekuaGfloq7wcuotxH+flyHgGG0ofkjD72QWYeBRzVNe1pymtaDu1nm5KkweonVD4HnJeZl9cujCRptE07VDLz6EEURJI0+iZ9+qsZN2VGIqL7MWBJ0jw0lUeKfxIRJ0TE2tPdeERsEBGnUN7pJUma56Zy+WtrylgqNzevmj8TOCczb+5esHnl/RbAjsCfADsA36fe01+SpDls0lDJzDuAP4yINwIfpQTMiRHxGGWgrPYgXWsDv03pn3I2sGNmXjaYokuS5pop36jPzIsp7/halzLuyQ6UV6SsTuk7chtwCfDPwI8y8976xZUkzWX9PP11O/CN5keSpGf19ZqWjohYA3gtsBS42Z7tkrRi6/fdXzsCnwFex3Ova1kaERcAB2fmtZXKJ0kaIf28++vNwI+Ax4CTKK+kX4kypO++wKURsYPBIkkrnn5aKp+iPPW1Q/fN+Ig4BrgCOJ4yvookaQXSz3gqWwNf6/V0V2beBXwVeONMCyZJY2NjA/nR4PQTKncBr1zO/N+i9F2RJK1g+gmV44ADIuJ5l7ci4nXAR4BjZlowSdLo6eeeyu8DdwNnRcQvgJ8DTwIbAtsDTwD7RMQ+rXXGM3PnmRZWkjS39RMqu1Bew3Ir5ZUs27Xm3dp8+lZiSVoB9dOj3sCQJPXUzz0VSZJ6MlQkSdUYKpKkagwVSVI1hookqRpDRZJUzYzGU5HafKeSJFsqkqRqDBVJUjWGiiSpGkNFklSNoSJJqsZQkSRVY6hIkqoxVCRJ1RgqkqRqDBVJUjWGiiSpGkNFklSNoSJJqsZQkSRVY6hIkqoxVCRJ1RgqkqRqDBVJUjWGiiSpGkNFklTNgmEXYLoiYnXg3h6z/k9m7j3b5ZEkPWfkQgXYqvl8G/Bwa/p9QyiLJKllFENlEXBXZp437IJIkpY1qqFyzbALIWlZY2Njwy6C5oBRDZXfRMRlwLaU+ysnAp/LzPGhlkySVnAjFSoR8QJgM+Ax4CDgVmA34Hjgt4Bjhlc6SdJIhQowBrwDuDUzb2ymXRgRC4FDIuKzmfmb4RVPklZsIxUqmbkU+Oces84F9gM2Aq6d1UJJkp41UqESEWtTWipnZuY9rVmrNJ+9+q9IkmbJqPWoXxk4GXh31/S9gOszc8nsF0mS1DFSLZXMvDki/gE4NiKeAa4D/pASKouHWjhJ0miFSuP9wJHAR4BXUYJlr8w8e6ilkiSNXqhk5uPA4c2PJGkOGbV7KpKkOcxQkSRVY6hIkqoxVCRJ1RgqkqRqDBVJUjWGiiSpGkNFklSNoSJJqsZQkSRVY6hIkqoxVCRJ1RgqkqRqDBVJUjWGiiSpmpEbT2VFMzY2NuwiSNKU2VKRJFVjqEiSqjFUJEnVGCqSpGoMFUlSNYaKJKkaQ0WSVI2hIkmqxlCRJFVjqEiSqjFUJEnV+O6vinxPlzQaBnGujo+PV9/mKLKlIkmqxlCRJFVjqEiSqjFUJEnVGCqSpGoMFUlSNYaKJKkaQ0WSVI2hIkmqxlCRJFVjqEiSqjFUJEnVGCqSpGoMFUlSNYaKJKkaQ0WSVI2hIkmqxlCRJFUzksMJR8QHgI8D6wJXAR/LzMuHWypJ0si1VCLiT4GvA98C9gIeBH4UEesPtWCSpNEKlYgYA44BTsnMozPzHGAP4F7go0MtnCRptEIF2Ah4LXB2Z0JmPgX8ENh1WIWSJBWjdk9l4+bzxq7pNwEbRsRKmbkUWAlgyZIls1k2FiwYtX9OSbXcfvvtwy7CjLW+M1fqdxuj9i24WvP5SNf0RyitrhcDDwOvAth3331nr2TABhtsMKv7kzR37LzzzsMuQk2vAn7Zz4qjFipjzef4BNOfaT6vBHYE7gSWzkK5JGk+WIkSKFf2u4FRC5WHms9Vgbta0xdSAuUxgMx8ArhkdosmSfNCXy2UjlG7UX9D89l9nWkDIDOzuwUjSZpFo9ZSuQG4DVgMnAcQES8E/gvlCbChmW6HzIjYAjgBeB1wP3AS8Nl2MEbEtcDmXavel5lrtJZZDBxLeTLueuATmfmDKpVatrxV6xcRy/sF4M8y8382j5A/RGmZtv1bZm7Xd2WeX9a+OtNGxKrAtcCBmfndrnk7Ap8HtgTuAI7PzG92LTMnj11rveXV7x3AkcBmwH2UJzI/kZmPNPNn5dg1+xpE/Ub23Gut17N+gz73Rqql0nwhfQbYLyKOi4jdgO8BawBfHFa5ptshMyJeAfwT5d7QHwGnAMcBB7aWeRHlabdDgd9v/byttcxOwHeBi4A9gWuAMyPi9XO9fl116vxcAvwK+L/NMutT/lO/t2u5PxtW3VrrrUr5v/eaHvM2Bc4FbgbeBXwf+PuI2Lu1zJw8dq31lle/t1BC5GfNNo8D/gQ4rbXYwI9dU5ZB1G9kz73WehPWjwGfe6PWUiEzvxoRqwAHUDo8XgW8LTNvGkZ5ujtkNtPOB7Ip31/3WO1DlH/7PTLz18A5EbEycFhEnND0vdkMeCHwvcz8xQS7/yRwfmbu3/z93Ih4LXA4pVPonK1fZl7RtZ/FwB8AO2Xm3c3kRZR7Zd9ttlNVn3UjIt5EOdFfOcGmDwVuAfZpfhE6NyLWBP6G8kUEc/fYTaV+BwGXZub7Wus8CJweEZtl5s8Z8LFr9jmo+o3yuTdp/QZ97o1US6UjM7+Qma/JzN/OzDcM+b1f/XTI3AW4oOtgnQW8HNi++fsi4Dc8dx9pGU2wvqG938b3gF0iou/nzLsMqn7PagLni8B3MvPC1qxFwC8H9aVE/51pzwJ+upxldgF+0HWP7yxgy4hYe44fu05Zl1e/KyiXM9uy+ez8Bj3oYweDq98on3swef2eNYhzbyRDZY6ZtEPmBOv0Wr69vUWUa9WnRcTDEfFQRPxd06yF8nDCggm2swrw6ulVY0KDql/bfsA6wCFd0xcBT0TEeRHx64i4JyI+29xHq6GfugHsmJl/BNzdPSMiXgysPcE2O/ucy8cOllM/gMw8NjO/0zV59+az85v9oI8dDKh+jPa5B5PXr636uWeozNxUOmT2WqfX8u3tLQLWAq6mPIhwBOWa6llT2G97/kwNqn4ARMQLKM340zLz1q51FgEbUn5T2xX4ErA/cPI0yr88/dSNzLy2z2125s/lYzdZ/Z4nIrYCDgPOyMzO46iDPnYwuPqN8rk35eM3qHNv5O6pzEFT7ZDZPW+iJzA6yx8CrNy6/vnjiLgb+E7zZFGnU+d09tuPQdWvYxfKb35/3GPZPwceycxrmr9fHBFPA8dHxNGZ+Z/LLfnk+qlbjW0OYr/9lmVGImIR5UnMO4APtmYN+tjB4Oo3yufedAzk3DNUZm5KHTJ7rNP9qN6qrXlk5n/0WO/c5nMroHP9s3s7C7vKNVMDqV/LYsq12590byQzL+2x7XMpTwBuCcz0i6mfuk3m4dY229rH5aEpLFPDIOr3rIh4M+W397uAXTLzvs68WTh2MKD6jfi5Nx0DOfcMlZlrd8hsX/tcXofMG+jdgRMgI2IB8G7g6q7/4Ks0n/dSrqs+M8F2HqU8HlhD9fp1Td+V556IelZEvATYG7iodUkFlv03mKl+6rZcmfloRNzJxPW/nhI8c/XYTUlE7AGcDlxHefry7ta82Th2MID6zYNzbzoGcu55T2Xm2h0ygWU6ZF4wwToXUJ4SaV8TXUy5OXhVZj4NHA0c1bXeXsBTwOWZ+ThwWXu/jXdS/jPUeudZ9fq1trMG5WmhK3i+JylPGHU/NrkX8ADl6ZaZ6qduU3EBsHvXjdTFwLWZedccP3aTiojfowTKlcCb2oHSmI1jBwOo3zw496ZkkOeeLZUZyszxiPgM8JWIeAC4FPgwrQ6ZEbEhsGbrGu1XKTe9zomIz1Ga1IcBh2bmk80yxwEnR8QJlM5z21P6OZzYup55PPDDiDgFOBP4r5QOSm8cgfoBbNHZTY/9Ph4Rfwt8PCLuo5zEbwU+BhyQmTNt+vdbt6n4POUL9x8j4huUa9fvpnQE7Zirx24qvkH5gv00sFlEtOddn5n3D/rYwUDrN8rn3lQN7NyzpVJBZn4VOBh4D6U5+VKW7ZB5JHB5a/k7KV80C5rlP0h5xcPnW8ucQrlZ9hbKf+wPUl4J8fHWMuc0+3wz5T/2ImBx7X47g6hf4xXN54MT7PpISkfCfYEfUH4T/KvM/MpM69Qq67TqNsVtXk15xHYDynHZHfjzzPzH1jJz8thNJiLWa8q6EDinWbf9s1NruwM9djCw4zey5940DOzcGxsf9x2MkqQ6bKlIkqoxVCRJ1RgqkqRqDBVJUjWGiiSpGkNFklSNoSJJqsYe9ZKe1QxA9VNgjcx86bDLo9FjS0VS2zHA7cMuhEaXoSIBEXFRRNwy7HIMU0RsC+xGea+V1Bcvf0nFcUwwkl4/ImI1ynuVxoD/l5mv67HMSyhvol2T8jr8l3ZeZx4RawF3At/PzD261juJ8ubadSu8/ryzzQWUF0V+qMb2tOKypSIBmXl+Zp41+ZJTti0lUB4HNo+IsR7LHMVzQ8b+R1dAbNOZ3l6h2c47ge/VCpTGgcA1mXlRxW1qBWSoSIOxbfN5JqUFtMyAThGxCaVV0Amyf+taf+vms3sUwu2BdVrrzVjz+vQPAQfV2qZWXF7+0rwREUdRXhP+zczcv8f8LwP7AYsy87queRcB62Xmes3fb6EMn3oJZSyYDSkDJn0pM0+aQnF+t/n8JmWsjS2B9ih6X6SM634hZYzwf+9af6JQ2ZMyzGxnSFsi4r9RXsselPE29qK8Iv0y4H2ZeVtEvAf4CLBpU49DWi2zHSmX4H7WjI3yImC1iFgCvCszL5tCfSXAlorml/OAfwU+HBEbt2dExOaUQDmpO1CW4+3AiZRxLD5KGRP8KxGx2xTW3Ra4FbiYMpJeZ1AkIuIdlKFcDwQ2byb3aqnc3xoUqmNP4IeZ+VRr2jaUy2xnAKsDnwK+TRnT5qSI+AqlFXIaZbCp1YFTI2LNZv3TKKMAbt38/AXwSPPn541fLi2PLRXNG5l5WUQcSWldbEUZD77jS5Qb50dNY5OvBrbOzGsAIuJMyvjj+1IGqOopIhYCGwNnZ+ZTEXEdpaXSGQ72C8CFmXlGRHyMMq759a31XwxsRKs10kzflNIaOaJrl9tQxg//dGZ+u7X8dpQBwn4IbNcJooh4qvn3WARc0AyP+3hrvfuB8cxcMpV/JKnNUNF88/Pmc5POhIjYk/Jb+36ZOdFId71kJ1CavyyJiLuAtSZZb2vKVYDOJa2rgN9r/nwA5VLa3hHxAkr4XZWZz7TW36pZv/vS12LgCcpluU7d1qUML3t2O1AaDwBPUy6BtVs2DzefT/cqfHOz3o6P6ouXvzSvZOYDwD00oRIRK1PGjL+K8sjsdNzTY9oTwEqTrNe5n9IOld+JiNdQhmk9OTN/SmnNLGTi+ynd0/cEzs/MR1vTOg8EnNajHFtQWkR3d03ftPl83vjk0kzZUtF8dD3PtVQOpDx59aau1sBUTHf5js4XfTtUFlC++J+m3NdoLzfpk19Ni2Q74ANdy3YePb6iPTEiXk25+X4Fz7ct8Csvb2kQbKloPkogImId4HDg9My8eBb3/7vAksy8s/n7Vc3n64FPZuZ9reWgd0vlMZa9J7QYGAfO7lp2G+DBzLypa3p3sHWv02u6NGOGiuajpPQN+d+U/+MHz9aOmxcybkLrS7u5j3MI8Anga63Ft6XcIL+utf5KlMtWV3e1rPYELsnM7ktyEwVEJ7CWaQVFxHrAyydYR5oxL39pPurcK3gLcFRm3jqL+96acs9lmS/tzPzsBMtenZlLW9M2oTzJ1b709TLgjZS+KLSmvxx4DXB6j21vC9ybmbf1mE53+aRabKloPuqEyq1Ary/zQZrSl3bTi/2lPZbr3E9p9w/ZnfILYHcv+s79lIlaKhNd+pq0fFK/xsbHa74+SBq+5gv7RuCIzDxu2OWZjog4hXIzfv3MvKWZdmbz962Xt640F3j5S/PRRI/kzmlND/c/Bv6lEyiNy5n+49DSUBgqmo8mem/WnBQRb6E8LvxXwAspr4R51gT3Y6Q5yVDRfLQNcNcI9cP4OvBK4FJg78wciTCUevGeiiSpGp/+kiRVY6hIkqoxVCRJ1RgqkqRqDBVJUjWGiiSpGkNFklSNoSJJqsZQkSRV8/8BbK0koWZ0X/sAAAAASUVORK5CYII=\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": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.015286631875968612 5.809850837172658 25.841321918303993\n"
     ]
    }
   ],
   "source": [
    "print(gamma_std_list.mean(),gamma_list.std(),gamma_list.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.06114652750387445 23.239403348690633 103.36528767321597\n"
     ]
    }
   ],
   "source": [
    "print(4*gamma_std_list.mean(),4*gamma_list.std(),4*gamma_list.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "380.06088485103436\n"
     ]
    }
   ],
   "source": [
    "print(gamma_list.std()/gamma_std_list.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "400.0\n"
     ]
    }
   ],
   "source": [
    "print(0.024/6e-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 1/1000 [00:37<10:29:11, 37.79s/it]Traceback (most recent call last):\n",
      "  File \"/Users/hstrey/anaconda3/envs/pymc3/lib/python3.6/site-packages/emcee/ensemble.py\", line 546, in __call__\n",
      "    return self.f(x, *self.args, **self.kwargs)\n",
      "  File \"/Users/hstrey/anaconda3/envs/pymc3/lib/python3.6/site-packages/lmfit/minimizer.py\", line 2206, in _lnpost\n",
      "    out = userfcn(params, *userargs, **userkwargs)\n",
      "  File \"<ipython-input-31-24fd3db35144>\", line 15, in residual\n",
      "    model = fitcorr(x,kappa,sigma,ga,eta)\n",
      "  File \"<ipython-input-18-ea6971c635aa>\", line 2, in fitcorr\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",
      "  File \"<ipython-input-18-ea6971c635aa>\", line 2, in <listcomp>\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",
      "KeyboardInterrupt\n",
      "  0%|          | 1/1000 [01:05<18:09:27, 65.43s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "emcee: Exception while calling your likelihood function:\n",
      "  params: [27.62722113]\n",
      "  args: (<function residual at 0x7fc650aab9d8>, Parameters([('kappa', <Parameter 'kappa', value=10 (fixed), bounds=[-inf:inf]>), ('sigma', <Parameter 'sigma', value=125 (fixed), bounds=[-inf:inf]>), ('ga', <Parameter 'ga', 27.62722112702247, bounds=[-inf:inf]>), ('eta', <Parameter 'eta', value=0.25 (fixed), bounds=[-inf:inf]>)]), ['ga'], array([[-inf,  inf]]))\n",
      "  kwargs: {'is_weighted': True, 'float_behavior': 'posterior', 'userargs': [], 'userkws': {'x': array([0.0000e+00, 1.0000e-05, 2.0000e-05, ..., 1.9997e-01, 1.9998e-01,\n",
      "       1.9999e-01]), 'data': array([0.29121829, 0.29080743, 0.29039779, ..., 0.00037773, 0.00037372,\n",
      "       0.00036959]), 'eps': array([0.04125825, 0.04125783, 0.04125701, ..., 0.03005998, 0.03006058,\n",
      "       0.03006088])}, 'nan_policy': 'omit'}\n",
      "  exception:\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-31-24fd3db35144>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     31\u001b[0m     \u001b[0;31m#result.params.add('__lnsigma', value=np.log(0.1), min=np.log(0.001), max=np.log(2))\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     32\u001b[0m     res = lm.minimize(residual, kws = {'x':t,'data':y,'eps':acf_std},method='emcee', nan_policy='omit', burn=300, steps=1000, thin=20,\n\u001b[0;32m---> 33\u001b[0;31m                      params=result.params)\n\u001b[0m\u001b[1;32m     34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     35\u001b[0m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit_report\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[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pymc3/lib/python3.6/site-packages/lmfit/minimizer.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(fcn, params, method, args, kws, iter_cb, scale_covar, nan_policy, reduce_fcn, calc_covar, **fit_kws)\u001b[0m\n\u001b[1;32m   2446\u001b[0m                        \u001b[0mnan_policy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnan_policy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreduce_fcn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mreduce_fcn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2447\u001b[0m                        calc_covar=calc_covar, **fit_kws)\n\u001b[0;32m-> 2448\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mfitter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mminimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/envs/pymc3/lib/python3.6/site-packages/lmfit/minimizer.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(self, method, params, **kws)\u001b[0m\n\u001b[1;32m   2115\u001b[0m                         val.lower().startswith(user_method)):\n\u001b[1;32m   2116\u001b[0m                     \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'method'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mval\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2117\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mkwargs\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   2118\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2119\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pymc3/lib/python3.6/site-packages/lmfit/minimizer.py\u001b[0m in \u001b[0;36memcee\u001b[0;34m(self, params, steps, nwalkers, burn, thin, ntemps, pos, reuse_sampler, workers, float_behavior, is_weighted, seed, progress)\u001b[0m\n\u001b[1;32m   1291\u001b[0m         \u001b[0;31m# now do a production run, sampling all the time\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1292\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mEMCEE_VERSION\u001b[0m \u001b[0;34m>=\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1293\u001b[0;31m             \u001b[0moutput\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msampler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_mcmc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msteps\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprogress\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mprogress\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   1294\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_lastpos\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcoords\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1295\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pymc3/lib/python3.6/site-packages/emcee/ensemble.py\u001b[0m in \u001b[0;36mrun_mcmc\u001b[0;34m(self, initial_state, nsteps, **kwargs)\u001b[0m\n\u001b[1;32m    384\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    385\u001b[0m         \u001b[0mresults\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 386\u001b[0;31m         \u001b[0;32mfor\u001b[0m \u001b[0mresults\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minitial_state\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0miterations\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnsteps\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\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    387\u001b[0m             \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    388\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pymc3/lib/python3.6/site-packages/emcee/ensemble.py\u001b[0m in \u001b[0;36msample\u001b[0;34m(self, initial_state, log_prob0, rstate0, blobs0, iterations, tune, skip_initial_state_check, thin_by, thin, store, progress)\u001b[0m\n\u001b[1;32m    343\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    344\u001b[0m                     \u001b[0;31m# Propose\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 345\u001b[0;31m                     \u001b[0mstate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maccepted\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmove\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpropose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstate\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    346\u001b[0m                     \u001b[0mstate\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom_state\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom_state\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    347\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pymc3/lib/python3.6/site-packages/emcee/moves/red_blue.py\u001b[0m in \u001b[0;36mpropose\u001b[0;34m(self, model, state)\u001b[0m\n\u001b[1;32m     90\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     91\u001b[0m             \u001b[0;31m# Compute the lnprobs of the proposed position.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 92\u001b[0;31m             \u001b[0mnew_log_probs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_blobs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompute_log_prob_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq\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     93\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     94\u001b[0m             \u001b[0;31m# Loop over the walkers and update them accordingly.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pymc3/lib/python3.6/site-packages/emcee/ensemble.py\u001b[0m in \u001b[0;36mcompute_log_prob\u001b[0;34m(self, coords)\u001b[0m\n\u001b[1;32m    427\u001b[0m                 \u001b[0mmap_func\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmap\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    428\u001b[0m             results = list(\n\u001b[0;32m--> 429\u001b[0;31m                 \u001b[0mmap_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_prob_fn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\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[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    430\u001b[0m             )\n\u001b[1;32m    431\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pymc3/lib/python3.6/site-packages/emcee/ensemble.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    544\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\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    545\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 546\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkwargs\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    547\u001b[0m         \u001b[0;32mexcept\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# pragma: no cover\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    548\u001b[0m             \u001b[0;32mimport\u001b[0m \u001b[0mtraceback\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/pymc3/lib/python3.6/site-packages/lmfit/minimizer.py\u001b[0m in \u001b[0;36m_lnpost\u001b[0;34m(theta, userfcn, params, var_names, bounds, userargs, userkws, float_behavior, is_weighted, nan_policy)\u001b[0m\n\u001b[1;32m   2204\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2205\u001b[0m     \u001b[0;31m# now calculate the log-likelihood\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2206\u001b[0;31m     \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0muserfcn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0muserargs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0muserkwargs\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   2207\u001b[0m     \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_nan_policy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnan_policy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnan_policy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhandle_inf\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2208\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-31-24fd3db35144>\u001b[0m in \u001b[0;36mresidual\u001b[0;34m(pars, x, data, eps)\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0meta\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mparvals\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'eta'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m     \u001b[0mmodel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfitcorr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mkappa\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mga\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0meta\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     16\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mdata\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-18-ea6971c635aa>\u001b[0m in \u001b[0;36mfitcorr\u001b[0;34m(t, kappa, sigma, ga, eta)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfitcorr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mkappa\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mga\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0meta\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----> 2\u001b[0;31m     \u001b[0mcorr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mgamma\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq_list\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mkappa\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mkappa\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msigma\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mga\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mga\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0meta\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mtt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkappa\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mq_list\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mq_list\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mga\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mtt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mt\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      3\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mcorr\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mfit_corr_model\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mModel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfitcorr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-18-ea6971c635aa>\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfitcorr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mkappa\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mga\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0meta\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----> 2\u001b[0;31m     \u001b[0mcorr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mgamma\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq_list\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mkappa\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mkappa\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msigma\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mga\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mga\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0meta\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mtt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkappa\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mq_list\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0msigma\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mq_list\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mga\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mtt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mt\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      3\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mcorr\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mfit_corr_model\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mModel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfitcorr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "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_list2[-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": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g_list"
   ]
  },
  {
   "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
}