{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1fed9816",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Before: 1488\n",
      "After: 4\n",
      "0.03769629626957275 0.010397710724117578\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Dysthe/opt/anaconda3/lib/python3.9/site-packages/trackpy/motion.py:235: FutureWarning: Using the level keyword in DataFrame and Series aggregations is deprecated and will be removed in a future version. Use groupby instead. df.median(level=1) should use df.groupby(level=1).median().\n",
      "  results = msds.mul(msds['N'], axis=0).mean(level=1)  # weighted average\n",
      "/Users/Dysthe/opt/anaconda3/lib/python3.9/site-packages/trackpy/motion.py:236: FutureWarning: Using the level keyword in DataFrame and Series aggregations is deprecated and will be removed in a future version. Use groupby instead. df.median(level=1) should use df.groupby(level=1).median().\n",
      "  results = results.div(msds['N'].mean(level=1), axis=0)  # weights normalized\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAG0CAYAAADdM0axAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPsklEQVR4nO3de1iUZf4/8PcwylEY5YyIgGcQRYVV8VB2ELMyyd+WpZKm1eqltsjXtlwz02wxt1r9tsFqpa5Z6q6m5jejuNb1FJWJoiJmHlA8DKJoHDRAZ57fH9NMDjMDM8zheWae9+u6uGoe7pm5kcO85z58boUgCAKIiIiIZMJL7A4QERERuRLDDxEREckKww8RERHJCsMPERERyQrDDxEREckKww8RERHJCsMPERERyUobsTsgNVqtFpcvX0ZgYCAUCoXY3SEiIiIrCIKA2tpadOzYEV5ezY/tMPw0cfnyZcTExIjdDSIiImqFCxcuoFOnTs22YfhpIjAwEIDuHy8oKEjk3hAREZE1ampqEBMTY3gdbw7DTxP6qa6goCCGHyIiIjdjzZIVLngmIiIiWWH4ISIiIllh+CEiIiJZ4ZofIiKSBK1Wi8bGRrG7QRLVtm1bKJVKhzwWww8REYmusbERZWVl0Gq1YneFJKx9+/aIjIy0uw4fww8REYlKEASo1WoolUrExMS0WKCO5EcQBNy6dQuVlZUAgKioKLsej+GHiIhEdefOHdy6dQsdO3aEv7+/2N0hifLz8wMAVFZWIjw83K4pMMZrIiISlUajAQB4e3uL3BOSOn04vn37tl2Pw/BDRESSwPMUqSWO+hnhtBcROZRGK+BA2XVU1tYjPNAXA+ODofTiixoRSQfDDxE5TH6JGot2lEJdXW+4FqXyxcIxiXgoyb4FikREjsJpLyJyiPwSNWasP2QUfACgoroeM9YfQn6JWqSekVxotAK+PVOF7cWX8O2ZKmi0gthdcqhz585BoVCguLhY7K64PYYfIrKbRitg0Y5SmHup0V9btKPU416MSDryS9QY9tYuPP3Bd/jjxmI8/cF3GPbWLqeFboVC0ezHlClTnPK8LYmLizP0wd/fH0lJSVi5cqVdjzlixAgoFAosXbrU5HMPP/wwFAoFXn/9dcO1s2fP4umnn0bHjh3h6+uLTp06YezYsfjpp58Mbe7+twoICED37t0xZcoUFBUV2dVXazH8EJHdDpRdNxnxuZsAQF1djwNl113XKZINMUYd1Wq14WP58uUICgoyurZixQqj9vbuTrLF4sWLoVarcfToUWRkZGD69OnYtGmTXY8ZExODNWvWGF27fPkydu3aZVRzp7GxESNHjkRNTQ0+++wznDx5Eps2bUJSUhKqq6uN7r9mzRqo1WocP34c77//Purq6jBo0CCsW7fOrr5ag+GHiOxWWWs5+LSmHZG1xBp1jIyMNHyoVCooFArD7fr6erRv3x7/+te/MGLECPj6+mL9+vWoqqrC008/jU6dOsHf3x99+vTBhg0bjB5Xq9XirbfeQrdu3eDj44POnTvjzTffNNsHrVaL559/Hj169MD58+cN1wMDAxEZGYlu3bphyZIl6N69O7Zt24Zz587By8sLBw8eNHqc9957D7GxsRAEy/9Gjz76KKqqqvDNN98Yrq1duxbp6ekIDw83XCstLcXZs2eRm5uLwYMHIzY2FkOHDsWbb76J3/3ud0aPqa/WHBcXh/T0dGzevBkTJ07ErFmzcOPGjZa/CXZg+CEiu4UH+jq0HZG1pDzq+PLLL+PFF1/EiRMnMGrUKNTX1yMlJQX/93//h5KSErzwwgvIzMzE999/b7jPvHnz8NZbb2HBggUoLS3Fp59+ioiICJPHbmxsxJNPPomDBw9i//79iI2NtdgPX19f3L59G3FxcXjwwQdNRnDWrFmDKVOmNLuN3NvbGxMnTjS679q1azF16lSjdmFhYfDy8sLmzZsN9ZtsMWfOHNTW1qKgoMDm+9qC4YeI7DYwPhhRKl9Y+tOpgG7X18D4YFd2i2RAyqOOWVlZGDduHOLj49GxY0dER0dj7ty56NevH7p06YLZs2dj1KhR+Pe//w0AqK2txYoVK7Bs2TJMnjwZXbt2xbBhw/Dcc88ZPW5dXR0eeeQRVFRUYPfu3UYjL3e7c+cO1q5di2PHjuGBBx4AADz33HPYsGEDGhoaAABHjhxBcXExnn322Ra/nmnTpuFf//oXbt68ib1796K6uhqPPPKIUZvo6Gj87//+L1577TV06NAB999/P9544w2cPXvWqn+zXr16AdAt7nYmhh8ispvSS4GFYxIBwCQA6W8vHJPIej/kcFIedUxNTTW6rdFo8Oabb6Jv374ICQlBu3bt8PXXX6O8vBwAcOLECTQ0NBiCiiVPP/006urq8PXXX0OlUpl8/uWXX0a7du3g5+eHmTNn4qWXXsIf/vAHAEBGRgbatGmDrVu3AgBWr16N++67D3FxcS1+PX379kX37t2xefNmrF69GpmZmWjbtq1Ju5kzZ6KiogLr169HWloa/v3vf6N3795Wjebop96cXfCS4YeIHOKhpCjkTRqASJXxi0ykyhd5kwawzg85hZRHHQMCAoxuv/POO/jb3/6GP/3pT9i1axeKi4sxatQoNDY2Avjt7KqWPPzwwzh69Ci+++47s59/6aWXUFxcjPPnz6Ourg7Lli0zHBbr7e2NzMxMrFmzBo2Njfj0009Npq6aM3XqVLz//vvYvHlzs/cLDAzEY489hjfffBNHjhzB8OHDsWTJkhYf/8SJEwCA+Ph4q/vUGixySEQO81BSFEYmRrLCM7mMftRxxvpDUABGC5+lNuq4b98+jB07FpMmTQKgW7B86tQpJCQkAAC6d+8OPz8//Oc//zGZ6rrbjBkzkJSUhMceewxffPEF7r33XqPPh4aGolu3bhbv/9xzzyEpKQm5ubm4ffs2xo0bZ/XXMGHCBMydOxfJyclITEy06j4KhQK9evVCYWFhi231O+cefPBBq/vUGgw/RORQSi8F0rqGiN0NkhH9qGPT6uKREqsu3q1bN2zZsgWFhYXo0KED3n33XVRUVBjCj6+vL15++WX86U9/gre3N4YOHYqrV6/i+PHjmDZtmtFjzZ49GxqNBo8++ii+/PJLDBs2zOp+JCQkYPDgwXj55ZcxdepUq0ecAKBDhw5Qq9Vmp7sAoLi4GAsXLkRmZiYSExPh7e2NPXv2YPXq1Xj55ZeN2v7888+oqKhAQ0MDfvrpJ6xcuRLbtm3DunXr0L59e6v71BoMP0RE5PbcYdRxwYIFKCsrw6hRo+Dv748XXngBGRkZRvVvFixYgDZt2uC1117D5cuXERUVhenTp5t9vKysLGi1Wjz88MPIz8/HkCFDrO7LtGnTUFhYaNOUl15zwaRTp06Ii4vDokWLDBWp9bfnzJlj1Fa/yNrX1xfR0dEYNmwYDhw4gAEDBtjcJ1sphOY29stQTU0NVCoVqqurERQUJHZ3iIg8Xn19PcrKyhAfHw9fX5ZDcIU333wTGzduxLFjx8Tuik2a+1mx5fWbC56JiIhkoq6uDj/88APee+89vPjii2J3RzQMP0RERDIxa9YsDBs2DPfee2+rprw8Bdf8EBERycTatWuxdu1asbshOo78EBERkaww/BAREZGsMPwQERGRrDD8EBERkaww/BAREZGsMPwQERG1giAIeOGFFxAcHAyFQoHi4mKMGDECWVlZLuvD2rVrm624rK+yXFxcDADYvXs3FAoFfv75Z5f0T6oYfoiIiFohPz8fa9euxf/93/9BrVYjKSkJn332Gd544w1Dm7i4OCxfvtzofi0FFmcaMmQI1Go1VCqVKM8vFazzQ0RE1ApnzpxBVFSU0ZlawcHBIvaoZd7e3oiMjBS7G6KT/MhPbm6u4QyPlJQU7Nu3z2LbKVOmQKFQmHz07t3bhT0mIiJPN2XKFMyePRvl5eWGwzsBGE17jRgxAufPn8ecOXMMr0e7d+/Gs88+i+rqasO1119/HQDQ2NiIP/3pT4iOjkZAQAAGDRqE3bt3Gz3v2rVr0blzZ/j7++Pxxx9HVVWVTf1uOu2lH4X66quvkJCQgHbt2uGhhx6CWq02ut+aNWuQkJAAX19f9OrVC7m5ubb+k0mKpEd+Nm3ahKysLOTm5mLo0KFYuXIlRo8ejdLSUnTu3Nmk/YoVK7B06VLD7Tt37iA5ORlPPPGEK7tNRET2EATg1i1xntvfH1C0fBL8ihUr0LVrV6xatQo//PADlEqlSZvPPvsMycnJeOGFF/D8888D0I0MLV++HK+99hpOnjwJAGjXrh0A3Snn586dw8aNG9GxY0ds3boVDz30EI4dO4bu3bvj+++/x9SpU/GXv/wF48aNQ35+PhYuXGj3l3zr1i28/fbb+Pjjj+Hl5YVJkyZh7ty5+OSTTwAAH3zwARYuXIi///3v6N+/Pw4fPoznn38eAQEBmDx5st3PLwZJh593330X06ZNw3PPPQcAWL58Ob766ivk5eUhJyfHpL1KpTKax9y2bRtu3LiBZ5991uJzNDQ0oKGhwXC7pqbGgV8BERHZ7NYt4NdA4HJ1dUBAQIvNVCoVAgMDoVQqLU4jBQcHQ6lUIjAw0KiNSqWCQqEwunbmzBls2LABFy9eRMeOHQEAc+fORX5+PtasWYO//OUvWLFiBUaNGoVXXnkFANCjRw8UFhYiPz/fnq8Yt2/fxj/+8Q907doVgO78r8WLFxs+/8Ybb+Cdd97BuHHjAADx8fEoLS3FypUr3Tb8SHbaq7GxEUVFRUhPTze6np6ejsLCQqse46OPPsKDDz6I2NhYi21ycnIMoUmlUiEmJsaufhMREdnq0KFDEAQBPXr0QLt27Qwfe/bswZkzZwAAJ06cQFpamtH9mt5uDX9/f0PwAYCoqChUVlYCAK5evYoLFy5g2rRpRv1asmSJoV/uSLIjP9euXYNGo0FERITR9YiICFRUVLR4f7VajS+//BKffvpps+3mzZuH7Oxsw+2amhoGICIiMfn760ZgxHpuEWi1WiiVShQVFZlMoemnxQRBcMpzt23b1ui2QqEwPJdWqwWgm/oaNGiQUTtzU33uQrLhR0/RZO5VEASTa+boF3FlZGQ0287Hxwc+Pj72dJGIiBxJobBq6skdeHt7Q6PRtHitf//+0Gg0qKysxPDhw80+VmJiIr777juja01vO1pERASio6Nx9uxZTJw40anP5UqSDT+hoaFQKpUmozyVlZUmo0FNCYKA1atXIzMzE97e3s7sJhERkUVxcXHYu3cvnnrqKfj4+CA0NBRxcXGoq6vDf/7zHyQnJ8Pf3x89evTAxIkT8cwzz+Cdd95B//79ce3aNezatQt9+vTBww8/jBdffBFDhgzBsmXLkJGRga+//tru9T7WeP311/Hiiy8iKCgIo0ePRkNDAw4ePIgbN24YzZy4E8mu+fH29kZKSgoKCgqMrhcUFBjVVDBnz549OH36NKZNm+bMLhIRETVr8eLFOHfuHLp27YqwsDAAukKD06dPx/jx4xEWFoZly5YB0G0nf+aZZ/A///M/6NmzJx577DF8//33hqUYgwcPxocffoj33nsP/fr1w9dff41XX33V6V/Dc889hw8//BBr165Fnz59cO+992Lt2rWIj493+nM7i0Jw1iSiA2zatAmZmZn4xz/+gbS0NKxatQoffPABjh8/jtjYWMybNw+XLl3CunXrjO6XmZmJU6dOtWo4sKamBiqVCtXV1QgKCnLUl0JERBbU19ejrKzMUNONyJLmflZsef2W7LQXAIwfPx5VVVVYvHixoXT4zp07Dbu31Go1ysvLje5TXV2NLVu2YMWKFWJ0mYiIiCRO0iM/YuDIDxGRa3Hkh6zlqJEfya75ISIiInIGhh8iIiKSFYYfIiKSBK7CoJY46meE4YeIiESlrxTc2Ngock9I6m79euBt06rUtpL0bi8iIvJ8bdq0gb+/P65evYq2bdvCy4vvy8mYIAi4desWKisr0b59e7uP1mD4ISIiUSkUCkRFRaGsrAznz58XuzskYe3bt0dkZKTdj8PwQ0REovP29kb37t059UUWtW3b1mGHqTL8EBGRJHh5ebHOD7kEJ1aJiIhIVhh+iIiISFYYfoiIiEhWGH6IiIhIVhh+iIiISFYYfoiIiEhWGH6IiIhIVhh+iIiISFYYfoiIiEhWGH6IiIhIVhh+iIiISFYYfoiIiEhWGH6IiIhIVhh+iIiISFYYfoiIiEhWGH6IiIhIVhh+iIiISFYYfoiIiEhWGH6IiIhIVhh+iIiISFYYfoiIiEhWGH6IiIhIVhh+iIiISFYYfoiIiEhW2ojdASIiZ9FoBRwou47K2nqEB/piYHwwlF4KsbtFRCJj+CEij5RfosaiHaVQV9cbrkWpfLFwTCIeSooSsWdEJDZOexGRx8kvUWPG+kNGwQcAKqrrMWP9IeSXqEXqGRFJAcMPEXkUjVbAoh2lEMx8Tn9t0Y5SaLTmWhCRHDD8EJFHOVB23WTE524CAHV1PQ6UXXddp4hIUiQffnJzcxEfHw9fX1+kpKRg3759zbZvaGjA/PnzERsbCx8fH3Tt2hWrV692UW+JpEejFfDtmSpsL76Eb89UefyIR2Wt5eDTmnZE5HkkveB506ZNyMrKQm5uLoYOHYqVK1di9OjRKC0tRefOnc3e58knn8SVK1fw0UcfoVu3bqisrMSdO3dc3HMiaZDjot/wQF+HtiMiz6MQBEGybwMHDRqEAQMGIC8vz3AtISEBGRkZyMnJMWmfn5+Pp556CmfPnkVwcHCrnrOmpgYqlQrV1dUICgpqdd+JxKZf9Nv0F1y/0Ttv0gCPDEAarYBhb+1CRXW92XU/CgCRKl/sf/l+bnsn8iC2vH5LdtqrsbERRUVFSE9PN7qenp6OwsJCs/f5/PPPkZqaimXLliE6Oho9evTA3Llz8csvv1h8noaGBtTU1Bh9ELk7OS/6VXopsHBMIoDfgp6e/vbCMYkMPkQyJtnwc+3aNWg0GkRERBhdj4iIQEVFhdn7nD17Fvv370dJSQm2bt2K5cuXY/PmzZg5c6bF58nJyYFKpTJ8xMTEOPTrIBKD3Bf9PpQUhbxJAxCpMp7ailT5euyIFxFZT9JrfgBAoTB+dyYIgsk1Pa1WC4VCgU8++QQqlQoA8O677+L3v/893n//ffj5+ZncZ968ecjOzjbcrqmpYQAit8dFv7oANDIxkhWeiciEZMNPaGgolEqlyShPZWWlyWiQXlRUFKKjow3BB9CtERIEARcvXkT37t1N7uPj4wMfHx/Hdp5IZFz0q6P0UiCta4jY3SAiiZHstJe3tzdSUlJQUFBgdL2goABDhgwxe5+hQ4fi8uXLqKurM1z76aef4OXlhU6dOjm1v0RSMjA+GFEqX5M1L3oK6HZ9DYxv3cYAIiJ3JtnwAwDZ2dn48MMPsXr1apw4cQJz5sxBeXk5pk+fDkA3ZfXMM88Y2k+YMAEhISF49tlnUVpair179+Kll17C1KlTzU55EXkqLvolIrJMstNeADB+/HhUVVVh8eLFUKvVSEpKws6dOxEbGwsAUKvVKC8vN7Rv164dCgoKMHv2bKSmpiIkJARPPvkklixZItaXQCQa/aLfpnV+Iu+q88NTz4lIjiRd50cMrPNDnsZSwJFjAUQi8ly2vH4z/DTB8ENyINcCiETkuTyiyCEROYecCyASEQEMP0SyI/cCiEREDD9EMsMCiEQkdww/RDLDAohEJHcMP0QywwKIRCR3DD9EMsMCiEQkdww/RDLEU8+JSM4kXeGZiJyHp54TkVwx/BDJGE89JyI54rQXERERyQrDDxEREckKww8RERHJCsMPERERyQrDDxEREckKww8RERHJCre6E5HNNFqB9YGIyG0x/BCRTfJL1Fi0oxTq6t9OfY9S+WLhmERWhiYit8BpLyKyWn6JGjPWHzIKPgBQUV2PGesPIb9ELVLPiIisx/BDRFbRaAUs2lEKwczn9NcW7SiFRmuuBRGRdDD8kNvRaAV8e6YK24sv4dszVXyxdZEDZddNRnzuJgBQV9fjQNl113WKiKgVuOaH3ArXm4instZy8GlNOyIisXDkh9wG15uIKzzQ16HtiIjEwvBDboHrTcQ3MD4YUSpfWNrQroBuFG5gfLAru0VEZDOGH3ILXG8iPqWXAgvHJAKASQDS3144JpH1fohI8hh+yC1wvYk0PJQUhbxJAxCpMp7ailT5Im/SAK67IiK3wAXP5Ba43kQ6HkqKwsjESFZ4JiK3xfBDbkG/3qSiut7suh8FdKMPXG/iGkovBdK6hojdDSKiVuG0F7kFrjchIiJHYfght8H1JkRE5Aic9iK3wvUmRERkL4Yfcjtcb0JERPZg+CGSEI1W4KgWkQfg77K0MfwQSQTPLSPyDPxdlj4ueCaSAJ5b1noarYBvz1Rhe/ElfHumikeckKj4u+weOPJDJLKWzi1TQHdu2cjESA6bN8F32CQl/F12Hxz5IRIZzy1rHb7DJqnh77L7kHz4yc3NRXx8PHx9fZGSkoJ9+/ZZbLt7924oFAqTjx9//NGFPSayDc8ts11L77AB3TtsToGRMzWdcq2o4e+yu5D0tNemTZuQlZWF3NxcDB06FCtXrsTo0aNRWlqKzp07W7zfyZMnERQUZLgdFhbmiu4StQrPLbOdLe+wWRaBWtKanVnmplyDA9pa9Xz8XRafpMPPu+++i2nTpuG5554DACxfvhxfffUV8vLykJOTY/F+4eHhaN++vYt6SWQfnltmO46WkaO0Zt2Yfsq16e/r9Zu3m30u/i5Lh2SnvRobG1FUVIT09HSj6+np6SgsLGz2vv3790dUVBQeeOAB/Pe//222bUNDA2pqaow+iFyJ55bZjqNlZCtzuwJbs26suSnXuzX3uwyAOxRFJtmRn2vXrkGj0SAiIsLoekREBCoqKszeJyoqCqtWrUJKSgoaGhrw8ccf44EHHsDu3btxzz33mL1PTk4OFi1a5PD+E9lCf25Z03egkdy5ZBZHy8gW5kZ3IoN8UH9Ha/POrJamXPU6BHjj+s3G357v199lABj21i7uUBSZZMOPnkJhnJ8FQTC5ptezZ0/07NnTcDstLQ0XLlzA22+/bTH8zJs3D9nZ2YbbNTU1iImJcUDPqSWsgGqM55ZZTz9aNmP9ISgAoxcwjpbR3SxNUVXUNDR7P0vrxqydSl3wSAIiVX5Gv8sFpRXm+/LrSBMPaHYdyYaf0NBQKJVKk1GeyspKk9Gg5gwePBjr16+3+HkfHx/4+Pi0up/UOqzPYh7PLbMeR8uoJdZOUTWnadixdio1UuVn9LvMGkDSItnw4+3tjZSUFBQUFODxxx83XC8oKMDYsWOtfpzDhw8jKop/BKXE4jsxvvshG3G0jJpj7RRVc0IDfPDtmSrDz1dKbIdWTblyh6K0SDb8AEB2djYyMzORmpqKtLQ0rFq1CuXl5Zg+fToA3ZTVpUuXsG7dOgC63WBxcXHo3bs3GhsbsX79emzZsgVbtmwR88ugu/DdDzkaR8vIEnt2+ykAqPzb4n/+fcSofk+UyhePJUdh1d4ym6ZcuUNRWiQdfsaPH4+qqiosXrwYarUaSUlJ2LlzJ2JjYwEAarUa5eXlhvaNjY2YO3cuLl26BD8/P/Tu3RtffPEFHn74YbG+BGqC736IyFmariMMbde6JQ36UPPzrdsAjLevV1TXY9XeMrxwTzw+P6K2esqVOxSlRSEIAvfY3aWmpgYqlQrV1dVGhRLJMbYXX8IfNxa32G7FU/0wtl+08ztEksEF8GSP5nZ0Vd+6bXGKSuXfFr5tlEajO/r76cKP+ftFqnyx56X7UHT+hlU/sxqtgGFv7Wpxumz/y/fz576VbHn9lvTID3kevvshc7gAnuxhaR3hlZoGwzVLU1RLx/UxWTemFQRM/PB7i8+nH6EuOn/D6hFq7lCUFskWOSTPpK/PYunXWwHdix7rs9jGXAE3d8EDSske1qwjbO/fFhFBxm+oIlW+hs0V+nVjY/tFI61rCK7VNb8NXs/W9Tn6HYqRKst98UiCAKjVwNdfA++8A0yZAnz1lahd4sgPuRTf/bSepWkhdx414QJ4spc16wh/vnUbn0wbAC8vhVVTVM4cofb4HYq1tUBJCXDs2G8fJSVAVZVxu6goYNQocfoIhh8SAeuz2M5SwNHvOnHXsgFcAE/2snb05drNBqvXETq7grhH7FC8fRv46SfjkHPsGHDunPn2Xl5A9+5Anz66j5EjXdrdploVfj7//HOb7zNy5Ej4+fm15unIA3n8ux8HsrSeQV1dj5V7y8zex11GTbj9l+zljFEajlDfRRCACxdMQ86PP+oCkDlRUb+FHP1HQgIgoQzQqvCTkZFhU3uFQoFTp06hS5curXk68lAe8e7HyeypUOsOoyZcAE/2ctYojSxHqG/cMA05JSWApQO/AwOBpCTjkJOUBIRI8+/N3Vo97VVRUYHw8HCr2gYGBrb2aYhkzREVaqU8asIDSslezhyl8dgR6oYG4MQJ06Bz6ZL59m3aAL16mY7mdO4MWDhrU+paFX4mT55s0xTWpEmTWDOHqBUcEVykPGrC6QVyBGeO0rj1CLVWC5SVmYacU6cAjcb8fTp3Ng05PXsC3t6u7buTschhEyxySFLy7ZkqPP3Bd626rzsVTXPnHWskHbIulFlZaTpddfw4cPOm+fYdOhhPVen/q1K5tt8OxCKHRB6ipWkhPXcZNbH04uSx0wvkUm49SmOtW7d0oabpaE5lpfn2Pj5AYqJp0OnY0W2nrBzB5vDzyy+/4Pr164iONt4yePz4cfTu3dthHSMi66aFbD1jSCwtje7I4oWLyFp37gCnT5uO5pw5o9uB1ZRCAXTpYjpl1a2bbs0OGbFp2mvz5s2YM2cOgoODIQgCPvjgAwwaNAgAMGDAABw6dMhpHXUVTnuRFLUUHJwx3O/Ix7S0XV//aFKvR0TkNPrqx01HckpLdQuTzQkPNx3J6d0bCAhwbd8lxmnTXkuWLMGhQ4cQFhaGgwcPYvLkyZg/fz4mTJgALh0icp6WpoUcPWriyDU4rOJM9KuaGvPVj69fN9/e3994K7n+/63caU2W2RR+bt++jbCwMABAamoq9u7di3HjxuH06dNQyHjukMgVXDUtZGmUprVVo1nFmWSnsRE4edI44Bw7Bpw/b769lxfQo4fplFV8vO5z5HA2hZ/w8HAcPXoUffv2BQCEhISgoKAAkydPxtGjR53SQSJyHWeM0rCKM3ksQdAFmqYh58cfdWt2zImONh3JSUgAfO0rSdHcNLVUPiel3Xg2hZ+PP/4YbZosnPL29saGDRswa9Ysh3aMiFzPGaM0rOJMHuH6dfPVj2trzbcPCjIOOPr/D3Z8wc7mpqkBSOJzjyVHmWzMELOchUPq/NTX1+Po0aOorKyEVqs1+txjjz1m78O7FBc8k5xtL76EP24sbrHdiqf6WX1IpEYrYNhbu1qs4uwO9YhIBn755bfqx3evz7l82Xz7tm3NVz+OiXHJVvLmNhNYenF39ecscfSGB5fW+cnPz0dmZiaqmh5XD92ZXhpLVSSJSHJ4SCTJhlYLnD1rvvpxkzfxBnFxplNWPXqIVv24pWlqS1z9uebuI9aGB7vDz6xZs/Dkk0/itddeQ0REhCP6REQi4SGR5JGuXDEdyTl+XFcw0JyQEPNbySU2G+CIs//EJtaGB7vDT2VlJbKzsxl8ZExKi9jIPjwkkqSu2b83dXW6UNN0O/nVq+YfzNfXuPqx/iMy0i2qH3vSJgFXfy12h5/f//732L17N7p27eqI/pCb4ZlMnoeHRJJU6f/eVN64ifjrl9Dz6nn8WHsRD+MaIs6f0k1jmaNQ6Cod3734WF/9WKl07RfhQJ60ScDVX4vdC55v3bqFJ554AmFhYejTpw/atm1r9PkXX3zRrg66Ghc8W49Vez0bR/RIdIIAXLoEHDuGkwXfoPTrQvS8eg5dqy7AR2NhK3lEhOlITmKirmCgh2lpM4E7cOSGB1tev+0OPx9++CGmT58OPz8/hISEGBU7VCgUOGspiUsUw4919L90luabuYOHiGxSXW26LqekBLhxw2zzm219cSq0M34Mi8PJsFhcjeuBFW8+A2WEvKof69+EAqbT1IKZ/xfjc5a49W6vV199FYsXL8Yrr7wCL1ailA1W7SWiVmls1BUBbLrL6sIF8+2VStyK64pdbSPwY1gsTobF4cewOFxUhUNQGL/mTKxTIk1my09bmqYGTOvuuPpzlur8iLnhwe7w09jYiPHjxzP4yAyr9hJRs7RaXfXjpouPT560XP04JsZ0XU6vXig4cc2q+lNy/XvT0mYCqXzuTw8lSGYq3e7wM3nyZGzatAl//vOfHdEfchOs2ktEBteumT+ws67OfHuVynRdTlIS0L692eb8e9Oy5jYTSOVzUtrwYHf40Wg0WLZsGb766iv07dvXZMHzu+++a+9TkAQ5qx6MO+GCYJKdX34BSktNp6wqKsy39/bWnVvVdDSnUyebtpLz7w05mt3h59ixY+jfvz8AoKSkxOhzPOndc8m9ai+3+JNH02iAM2dMQ86ZM5arH8fHm47mdO+uO/7BTnL/e0OO55CzvTwJd3vZRo4hgFv8yWMIwm/Vj+/+KC3VjfKYExpqemBn795AYKDTuyvHvzdkPZdudb9w4QJiYmLseQhJYfixnZymf7jFn9xWXZ3pupxjxwAz5zIC0FU/7t3bdF2OyNWP5fT3hmzj0q3usbGx6NChA5KTk5GcnIx+/fohOTkZDQ0NeP/997Fu3Tp7n4IkTkqL2JyNW/xJ8m7f1h3O2TTklJWZb+/lpat03PQsq65dJVn9WE5/b8h57A4/Z8+eRXFxMYqLi3H48GFs3rwZly9fBgCOnJDHvUvjFn+SDEEALl40DTk//qirpWNOZKT56sd+fq7tO5HI7A4/cXFxiIuLQ0ZGhuHat99+i8mTJ+Ott96y9+HJjXni/Ly1W2lDA3zw7Zkqjwl9JLIbN8xvJa+uNt++XTvTHVZJSbr1OiRpnvaGUaqctuD5yy+/xKuvvoqioiJnPLzTcM2PY3jqouCWztJRAFD5t4VvGyUqajwn9JGLNDQAJ06YBp2LF823b9MG6NnTNOTExuqms8iteOIbRldy6YLn27dvm9T2AYDTp08jKSkJ9fXuNfzP8GM/T18UbM1ZOk25e+gjB9NqgXPnTKesfvpJt83cnJgY0ymrXr10tXTI7XnqG0ZXcumC54CAACQmJqJ///7o168f+vfvj44dO+K9995Denq6vQ9PbsjTFwVbOksnIsgH9Xe0+PnWbZP7CND9EVu0oxQjEyPdMvRRK129ajpdVVIC3Lxpvn379qYjOc1UPyb3p9EKWLSj1OybJ/7tcA67w8+uXbtw5MgRHDlyBJ988gn+/Oc/45df60Okp6dj/vz56Nu3L/r27YuEhAS7O0zSJ4dFwebO0tEKAiZ++L3F+7h76KMW3LplvvrxlSvm2+urHzcdzYmOtmorOdeGeA5Pf8MoRXaHn2HDhmHYsGGG21qtFidPnjTsACsqKsLq1atRWVkJjaXhXPIocjmHp+mW2+3Fl6y6nzuHPoJuWur0adPRnNOndTuwzOnSxXz14zat+xPMtSGeRQ5vGKXG7vDTlJeXFxISEpCQkICnn37acP2KpXc/LcjNzcVf//pXqNVq9O7dG8uXL8fw4cNbvN8333yDe++9F0lJSSguLm7Vc1PryPUcHrmEPtkQBECtNg05paWApbWMYWGmU1a9e+t2XzmIpbUhFdX1mLH+ENeGuCH+7XC9VoWfo0ePIikpCV5W7iY4fvw4evbsafPzbNq0CVlZWcjNzcXQoUOxcuVKjB49GqWlpejcubPF+1VXV+OZZ57BAw880OrQRa0n13N45Br6PEJtrfnqx9evm2/v729c/Vi/rTwiwqnd5NoQz2Tt346U2A4soeEgrdrtpVQqUVFRgbCwMKvaBwUFobi4GF26dLHpeQYNGoQBAwYgLy/PcC0hIQEZGRnIycmxeL+nnnoK3bt3h1KpxLZt25od+WloaEBDQ4Phdk1NDWJiYrjbywHkODTf3E4wgDs2RHf7NnDy5G+jOPqQc+6c+fZeXrrpqaZTVl26iLKV/NszVXj6g+9abLfh+cFuvzZEbmuaWvrb8cI98fj8iFpWf09t5fTdXoIgYMGCBfD397eqfaOlaqMt3KeoqAivvPKK0fX09HQUFhZavN+aNWtw5swZrF+/HkuWLGnxeXJycrBo0SKb+0ctM7co2NP/gFnaCRbJP1KuJQhAeblpyPnxR10AMqdjR9MDOxMSJFX9WC5rQ+T4xqm5vx2PJUdh1d4yTnU6UKvCzz333IOTJ09a3T4tLQ1+Nv4BuXbtGjQaDSKaDCNHRESgoqLC7H1OnTqFV155Bfv27UMbKxcSzps3D9nZ2Ybb+pEfcgw5nsMjx9Anqhs3TKerSkqAmhrz7QMDzVc/DpH+z6kc1obIeU2Tub8dKbEdcO9f/8upTgdrVfjZvXu3g7thmaLJlk9BEEyuAYBGo8GECROwaNEi9OjRw+rH9/HxgY+Pj939JLqbHEOf09XXm69+fMnCLrs2bXRFAJsGndhYUU8lt4enryvjmibTvx3fnqniNngncPhuL0cJDQ01rC26W2VlpcloEADU1tbi4MGDOHz4MGbNmgVAt+1eEAS0adMGX3/9Ne6//36X9J2I7KDVAmfPmk5ZnTplufpx586m63J69nT76sfm1r148mYC1rsxJZepTleTbPjx9vZGSkoKCgoK8PjjjxuuFxQUYOzYsSbtg4KCcOzYMaNrubm52LVrFzZv3oz4+Hin95mIbFRZaTpldfy4rmCgOcHBputykpIAD9yc0Ny6F09dV8YXelNymOoUg2TDDwBkZ2cjMzMTqampSEtLw6pVq1BeXo7p06cD0K3XuXTpEtatWwcvLy8kJSUZ3T88PBy+vr4m14nIxW7e1IWaplNWlZXm2/v4AImJpqM5UVFuO2VlC2vWvex/+X6L68rcdacUX+hNefpUp1gkHX7Gjx+PqqoqLF68GGq1GklJSdi5cydiY2MBAGq1GuXl5SL3kogM7twxrX587JhuGstcVQ2FAuja1XRdTrdura5+7O5sWfdiburHnXdK8YXelFzrpjmb3ae63+3GjRsICgqCUql01EO6HE91J7KCIACXL5uGnBMngLvqZhkJDzcdyUlMBAICXNt3ibOnlo8nnAzOWlnmuXOodRWXnup++vRpbN++Hdu3b8e3336LwMBAPPzwwxg7dixGjx6Ndg4s605EIqiu/u0k8ruDzo0b5tv7+/92EvndQSc83LX9dlOtXffiKTulWCvLPJbQcKxWhZ8LFy4gLy8P27dvR3l5OR588EFMnToVW7duxcWLF/H5559j6dKleOaZZzBixAg89thjmDFjhqP7TkSO1Nj4W/Xjuz8sTS0rlearH8fHi1L92FO0dt2LO+6UsrQ2iS/05rGEhuO0KvwcOnQIFRUVyMnJQXp6Onx9f/slDAkJQXJyMhYsWICLFy9i69at2Lp1K8MPkVQIAnD+vGnIOXlSt2bHnOho05DTqxfgK5+Fp67S2nUv7rZTqqVpHL7QkzM5dM2PJ+CaH/IoVVWm01UlJbqDPM0JCjINOUlJQIcOru23zLVm3Ys7nfvlCWuTSHpcuuaHiCTgl190i42bjuao1ebbt22rG7lpGnRiYmSxlVxKzE39tGbdi7vslPKUtUnk3mwOP7/88guuX7+O6Ohoo+vHjx9H7969HdYxIjJDo/mt+vHdIzmnTukqI5sTF2daGLBHD7evfuwJWpr6sWXdi1S3RDcNd1pBcLu1SeR5bJr22rx5M+bMmYPg4GAIgoAPPvgAgwYNAgAMGDAAhw4dclpHXYXTXiQZV66YjuSUllqufhwSYjpd1bu3R1Y/9gTOmvqR0pZoc31p79cWP/9yu8X7rniqH8b2i26xHZGe06a9lixZgkOHDiEsLAwHDx7E5MmTMX/+fEyYMAFcOkTUSnV1uurHTUdzrl41397X13z148hITlm5CWdO/Uhlp5SlcGdN8AHkVcWZXM+m8HP79m2EhYUBAFJTU7F3716MGzcOp0+fNnvSOhHd5c4d4KefTA/sPHvWfHuFQlfpuOloTrduum3m5LacvS1d7J1SzYW7lkhlbRJ5NpvCT3h4OI4ePYq+ffsC0G1rLygowOTJk3H06FGndJDI7QgCcOmS+erHjY3m7xMRYb76sb+/a/tOLuFu29Jt1VK4s4THNZCr2BR+Pv74Y7Rpct6Ot7c3NmzYgFmzZjm0Y0Ru4eefzW8l//ln8+0DAkxPJO/TB/h1RJXkwdMP8LQ2tDVd/yP3Ks7kOjaFn06dOhnd1mg02LZtG9LT0zF06FCHdoxIUhobgR9/NB3NuXDBfHulEujZ0/SIh7g4Vj8mUbelu+LEd2tD2/sTBsDLS8EqzuRydtX5USqVmDRpEo4fP47AwEBH9YlIPFqtafXjkpLmqx936mS++rGPj2v7TpJlLnCIsS3dVTvBrA13g7uGMOyQKOwucjhw4ECUlZWhS5cujugPketcu2YccPT/rasz316lMp2uYvVjakFzgcOVB3ha2n1VUV2PGesP2VVVWSrhjshadh9vsWXLFixbtgybN29GTEyMo/olGtb58UC3bunq4zTdZVVRYb69tzeQkGBaGLBTJ24lJ5tYU8vHFdvSNVoBw97aZXERsn4kZv/L99v83M2FOwCSqTlEns+W12+7w4/Xr+sX2rVrh8ceewwjRoxA//790adPH3i7YQVZhh83ptEAZ86Yrss5c8Zy9eP4eNMpq+7ddcc/ENnBmYHDVs4690sq4Y4IcPHZXmVlZSguLsaRI0dQXFyMnJwcnDt3DkqlEr169eIWeHI8QdCN2jSdsiot1Z1xZU5oqOlITu/eANeqkZM4u5aPLZyxtd6WQo08poKkxu7wExsbi9jYWIwdO9Zwrba2FkeOHMHu3bvtfXiSu7o6063kx47pTis3x89PF2qaBp2ICMlMWblitw25XtPva0WNdGr5OGNrvZTCHZGtHH6qe3V1NT755BN89NFHKC4uxquvvuropyBPdPu2rvpx06BTVma+vUKhm55qOmXVpYukqx9L6dwlchxz39fgAOumTl1Ry8cZW+s9vVAjeTaHhZ9du3Zh9erV2LJlCwIDAzF8+HAUFxc76uHJUwgCcPGi6UjOjz9arn4cFWU6kpOYqBvlcSPO3G1D4rH0fb1+s/kzrFx5jIMzTnz39EKN5NnsCj8XL17E2rVrsWbNGly5cgVjx47Fli1bMGrUKJw4cQLbtm1zUDfJLf38s2nIKSkBqqvNtw8M1AWcphWQQ0Nd2m1ncOZBluQ6Tae2UmI7WHWGlRS2ez+UFNXi1vrmpmTNfe1iFWokslerw8/DDz+M//73v7j//vuxePFiZGRkICAgwPB5HnQqIw0NunOrmm4lv3jRfPs2bXTVj5tOWcXGSmZdjqNxfYT7sBQALE1ttTTCAwAdArxx/eZvI5tiHePQ3Invrdmy/lhyFFbtLZNEuCOyRavDT35+PiZMmICsrCykpqY6sk8kVVqtbg1O03U5P/2k22ZuTufOpkc89Owpu+rHXB/hHiwFAP2LvK1TW3oLHklApMpPEovczZ343tyU7PT1h8w+TkV1PVbtLcML98Tj8yNqlxRqJHKUVoefb775BqtXr8b999+PqKgoTJw4ERMmTEC3bt0c2T8Sy9WrplNWx48DN2+ab9++velITlKSriqyRIi5y4rrI6TF3M9CQWmF2QCgrq7Hyr0WFt5bKVLlJ9kRvZamZC3RT9d+fkSNPS/dh6LzNyQR7ois0erwk5aWhrS0NKxYsQIbN27E6tWrsWjRIvzud7/DxIkT0bt3b0f2k5zl7urHd39cuWK+vbe3brFx09Gc6GhJT1mJvctKzIMsyZi5n4XIIB/U39G2uHbHVu7wfW1pSrY5+unaovM3JBvuiMyxe7eXv78/pk6diqlTp+LkyZP46KOP8Je//AVXrlzhuh8p0WiA06fNVz+2VOS7Sxfz1Y/bOLxCglNJYZeVM3bbkO0s/izUNDj8udzl++qIqVZO15K7ceirWM+ePbFs2TLk5ORgx44dWL16tSMfnqzRtPqx/qO0FKi38AcqLMw05CQmAu3aubbvTiClXVbW7LYh52nuZ8ERgiWyqNlWjphq5XQtuZtWhZ/Kykr4+fkh0MLRAEqlEhkZGcjIyAAA/PTTT+jRo0erO0kW1NaaLj4uKbFc/djf/7fqx3d/hIe7tt8uJLVdVs3ttiHnsmd6pzn6qS13XffS0pRsc9xhWo/InFaFnz179iAzMxMjRozA2LFjMXbsWHTs2NGozffff4/t27dj69atuHDhAurq6hzSYVlqbNTtqGo6mnP+vPn2Xl6Wqx//ehCtXEhxl5W53TbkfI74Hjc3Zendxsstv68tTckKZv5ffxuQ/rQekTmtCj9PPPEEBg8ejO3bt2PLli3IyspCcnIyxowZg0uXLuHzzz+HRqPBI488gqVLlyI9Pd3R/fZMggCUl5uGnJMndcc/mNOxo2nISUgAfB07DO2u51FxlxXptfZ7rP8p9+Qt3S1NyQKmdX485WsneVIIgqXVrtarrq7Gjh078MUXX6Bjx47IyMjA0KFD4eWGoww1NTVQqVSorq5GUFCQ857oxg3z1Y9rasy3Dww0PeKhTx8g2PnDzWLvlLKHRitg2Fu7Wtxltf/l+90izFHrWfOzoPJvC982SqNDSaOsrIDsCWyp8OxpXzu5P1tevx0SfjyJ08LP7t3AF1/8FnQuXzbfrk0boFcv03o5IlU/trQ7Rt8TdziPSv81AOaH7d3hayDHsOZngWuyiNwTw48dnBZ+FiwAliwxvhYbazpl1aOHrpaOBOjfKVtaJOpOoyb2jF7xHa9nceeRTCKyzJbXb/cq2OLO7r9fd9CnPuT07g04c1rNAaS2U8oerd1lxRdKz8Mdd0TE8OMq992n+3AjUtwpZY/mdlnZctyBK4sjknNwxx2RvDH8kEHTABDazrrDR919p5Stxx24ujgiERE5FsMPAbAcANr7t0X1rdseex5Va487cKcpPyIiMib5vei5ubmIj4+Hr68vUlJSsG/fPott9+/fj6FDhyIkJAR+fn7o1asX/va3v7mwt+5JHwCaru+5UtOAn38NPk3HNjyhwJkjjjtwlyk/IiL6jaRHfjZt2oSsrCzk5uZi6NChWLlyJUaPHo3S0lJ07tzZpH1AQABmzZqFvn37IiAgAPv378cf/vAHBAQE4IUXXhDhK5A+a86+Mlf7xBMKnDniuAN3n/IjIpIjSW91HzRoEAYMGIC8vDzDtYSEBGRkZCAnJ8eqxxg3bhwCAgLw8ccfW9XeZUUOJeLbM1V4+oPvWmz3ybRB8PJSeNTumO3Fl/DHjcWtuq87bfMnIpIDW16/JTvt1djYiKKiIpOjMdLT01FYWGjVYxw+fBiFhYW49957LbZpaGhATU2N0YecWDttc+1mA9K6hmBsv2ikdQ3xiBd8e487cOcpPyIiOZNs+Ll27Ro0Gg0iIiKMrkdERKCioqLZ+3bq1Ak+Pj5ITU3FzJkz8dxzz1lsm5OTA5VKZfiIiYlxSP/dhZzPvtKfZm0pvigAtPdvi8gg4689UuXLbe5ERG5M0mt+AEDR5EgHQRBMrjW1b98+1NXV4bvvvsMrr7yCbt264emnnzbbdt68ecjOzjbcrqmpkVUA0geAls6+cucdXZa0dJo1ACwd14cF8YiIPIxkw09oaCiUSqXJKE9lZaXJaFBT8fHxAIA+ffrgypUreP311y2GHx8fH/j4WFfPxhNZEwA8eXqnpdOs9aM73M5OROQ5JBt+vL29kZKSgoKCAjz++OOG6wUFBRg7dqzVjyMIAhoamq/ZInfWBgB30JpzuHjcARGRvEg2/ABAdnY2MjMzkZqairS0NKxatQrl5eWYPn06AN2U1aVLl7Bu3ToAwPvvv4/OnTujV69eAHR1f95++23Mnj1btK/BXXhCALDnHC4ed0BEJB+SDj/jx49HVVUVFi9eDLVajaSkJOzcuROxsbEAALVajfLyckN7rVaLefPmoaysDG3atEHXrl2xdOlS/OEPfxDrS3Ar7hwALFZq5jlcRETUhKTr/IjB0+v8tGZaSOo0WgHD3tplsWAha/IQEXk+W16/JT3yQ45lz7SQlLVUqZnncBER0d0kW+eHHMvS+V36aaH8ErVIPbOftYUaeQ4XEREBDD+y0NL5XQCwaEcpNFr3nAGVc6FGIiKyHcOPDNgyLeSOrKnUHOWhhRqJiMh2DD8y4OnTQvpCjQBMApAcCjUSEZFtGH5kQA7TQvpCjZEqnsNFRETN424vGZDL+V2eUKiRiIicj+FHBuR0fpc7F2okIiLX4LSXTHBaiIiISIcjPzLCaSEiIiKGH9nhtBAREckdww8RSZ4nnklHROJh+CEiSfPUM+mISDxc8ExEkuXJZ9IRkXgYfohIkjz9TDoiEg/DDxFJkqefSUdE4mH4ISJJ8vQz6YhIPAw/RCRJcjiTjojEwfBDRJKkP5PO0oZ2BXS7vtz9TDoicj2GHyKSJP2ZdABMApCnnUlHRK7F8ENEksUz6YjIGVjkkCSJFX1Jj2fSEZGjMfyQ5LCiLzXFM+mIyJE47UWSwoq+RETkbAw/JBms6EtERK7A8EOSwYq+RETkCgw/JBms6EtERK7A8EOSwYq+RETkCgw/JBms6EtERK7A8EOSwYq+RETkCgw/JCms6EtERM7GIockOazoS0REzsTwQ5LEir5EROQsnPYiIiIiWeHIjwfioaBERESWMfx4GB4KSkRE1DxOe3kQHgpKRETUMoYfD8FDQYmIiKwj+fCTm5uL+Ph4+Pr6IiUlBfv27bPY9rPPPsPIkSMRFhaGoKAgpKWl4auvvnJhb8XDQ0GJiIisI+nws2nTJmRlZWH+/Pk4fPgwhg8fjtGjR6O8vNxs+71792LkyJHYuXMnioqKcN9992HMmDE4fPiwi3vufBqtgG/PVGF78SV8e6YKFTU8FJSIiMgaCkEQJDsPMmjQIAwYMAB5eXmGawkJCcjIyEBOTo5Vj9G7d2+MHz8er732mtnPNzQ0oKGhwXC7pqYGMTExqK6uRlBQkH1fgJOYW9QcHNAW12/ebvG+G54fzPo5RETkcWpqaqBSqax6/ZbsyE9jYyOKioqQnp5udD09PR2FhYVWPYZWq0VtbS2Cgy0fhJmTkwOVSmX4iImJsavfzmZpUXNLwYeHghIREelINvxcu3YNGo0GERERRtcjIiJQUVFh1WO88847uHnzJp588kmLbebNm4fq6mrDx4ULF+zqtzM1t6j5bjwUlIiIyDLJ1/lRKIxfrAVBMLlmzoYNG/D6669j+/btCA8Pt9jOx8cHPj4+dvfTFVpa1KzXIcAb1282Gm5Hss4PERGRgWTDT2hoKJRKpckoT2VlpcloUFObNm3CtGnT8O9//xsPPvigM7vpUtYuVl7wSAIiVX6s8ExERGSGZKe9vL29kZKSgoKCAqPrBQUFGDJkiMX7bdiwAVOmTMGnn36KRx55xNnddKnwQF+r2kWq/JDWNQRj+0UjrWsIgw8REdFdJDvyAwDZ2dnIzMxEamoq0tLSsGrVKpSXl2P69OkAdOt1Ll26hHXr1gHQBZ9nnnkGK1aswODBgw2jRn5+flCpVKJ9HY4yMD4YUSpfVFTXm133o4BuiouLmomIiCyT7MgPAIwfPx7Lly/H4sWL0a9fP+zduxc7d+5EbGwsAECtVhvV/Fm5ciXu3LmDmTNnIioqyvDxxz/+UawvwaGUXgosHJMIgIuaiYiIWkvSdX7EYEudALHw8FIiIiJjtrx+S3rai8x7KCkKIxMjcaDsOhc1ExER2YjhR8I0WsFiwFF6KVipmYiIqBUYfiSKU1tERETOIekFz3Jl6QiLiup6zFh/CPklapF6RkRE5P4YfiSmuSMs9NcW7SiFRst16kRERK3B8CMxLR1hIQBQV9fjQNl113WKiIjIgzD8SIy1R1hY246IiIiMMfxIjLVHWFjbjoiIiIwx/EiM/ggLSxV7FNDt+uIRFkRERK3D8CMxPMKCiIjIuRh+JOihpCjkTRqASJXx1Fakyhd5kwawzg8REZEdWORQoniEBRERkXMw/EgYj7AgIiJyPE57ERERkaww/BAREZGsMPwQERGRrHDNj8g0WoGLmomIiFyI4UdE+SVqLNpRanSWV5TKFwvHJHI7OxERkZNw2ksk+SVqzFh/yOQQ04rqesxYfwj5JWqRekZEROTZGH5EoNEKWLSjFIKZz+mvLdpRCo3WXAsiIiKyB8OPCA6UXTcZ8bmbAEBdXY8DZddd1ykiIiKZYPgRQWWt5eDTmnZERERkPYYfEYQH+rbcyIZ2REREZD2GHxEMjA9GlMrX5NR2PQV0u74Gxge7sltERESywPAjAqWXAgvHJAKASQDS3144JpH1foiIiJyA4UckDyVFIW/SAESqjKe2IlW+yJs0gHV+iIiInIRFDkX0UFIURiZGssIzERGRCzH8iEzppUBa1xCxu0FERCQbnPYiIiIiWWH4ISIiIllh+CEiIiJZYfghIiIiWWH4ISIiIllh+CEiIiJZYfghIiIiWWH4ISIiIllh+CEiIiJZkXz4yc3NRXx8PHx9fZGSkoJ9+/ZZbKtWqzFhwgT07NkTXl5eyMrKcl1HiYiIyC1IOvxs2rQJWVlZmD9/Pg4fPozhw4dj9OjRKC8vN9u+oaEBYWFhmD9/PpKTk13cWyIiInIHCkEQBLE7YcmgQYMwYMAA5OXlGa4lJCQgIyMDOTk5zd53xIgR6NevH5YvX27Tc9bU1EClUqG6uhpBQUGt6TYRERG5mC2v35Id+WlsbERRURHS09ONrqenp6OwsNBhz9PQ0ICamhqjDyIiIvJckg0/165dg0ajQUREhNH1iIgIVFRUOOx5cnJyoFKpDB8xMTEOe2wiIiKSHsmGHz2FQmF0WxAEk2v2mDdvHqqrqw0fFy5ccNhjExERkfS0EbsDloSGhkKpVJqM8lRWVpqMBtnDx8cHPj4+Dns8IiIikjbJjvx4e3sjJSUFBQUFRtcLCgowZMgQkXpFRERE7k6yIz8AkJ2djczMTKSmpiItLQ2rVq1CeXk5pk+fDkA3ZXXp0iWsW7fOcJ/i4mIAQF1dHa5evYri4mJ4e3sjMTFRjC+BiIiIJEbS4Wf8+PGoqqrC4sWLoVarkZSUhJ07dyI2NhaArqhh05o//fv3N/x/UVERPv30U8TGxuLcuXOu7DoRERFJlKTr/IiBdX6IiIjcj0fU+SEiIiJyBoYfIiIikhWGHyIiIpIVhh8iIiKSFYYfIiIikhWGHyIiIpIVhh8iIiKSFYYfIiIikhWGHyIiIpIVhh8iIiKSFUmf7eVJNFoBB8quo7K2HuGBvhgYHwyll0LsbhEREckOw48L5JeosWhHKdTV9YZrUSpfLByTiIeSokTsGRERkfxw2svJ8kvUmLH+kFHwAYCK6nrMWH8I+SVqkXpGREQkTww/TqTRCli0oxSCmc/pry3aUQqN1lwLIiIicgaGHyc6UHbdZMTnbgIAdXU9DpRdd12niIiIZI7hx4kqay0Hn9a0IyIiIvsx/DhReKCvQ9sRERGR/Rh+nGhgfDCiVL6wtKFdAd2ur4Hxwa7sFhERkawx/DiR0kuBhWMSAcAkAOlvLxyTyHo/RERELsTw42QPJUUhb9IARKqMp7YiVb7ImzSAdX6IiIhcjEUOXeChpCiMTIxkhWciIiIJYPhxEaWXAmldQ8TuBhERkexx2ouIiIhkheGHiIiIZIXhh4iIiGSF4YeIiIhkheGHiIiIZIXhh4iIiGSF4YeIiIhkheGHiIiIZIXhh4iIiGSFFZ6bEAQBAFBTUyNyT4iIiMha+tdt/et4cxh+mqitrQUAxMTEiNwTIiIislVtbS1UKlWzbRSCNRFJRrRaLS5fvozAwEAoFLYfPFpTU4OYmBhcuHABQUFBTugh2YPfH2nj90fa+P2RNrl/fwRBQG1tLTp27Agvr+ZX9XDkpwkvLy906tTJ7scJCgqS5Q+fu+D3R9r4/ZE2fn+kTc7fn5ZGfPS44JmIiIhkheGHiIiIZIXhx8F8fHywcOFC+Pj4iN0VMoPfH2nj90fa+P2RNn5/rMcFz0RERCQrHPkhIiIiWWH4ISIiIllh+CEiIiJZYfghIiIiWWH4caDc3FzEx8fD19cXKSkp2Ldvn9hdol/t3bsXY8aMQceOHaFQKLBt2zaxu0S/ysnJwe9+9zsEBgYiPDwcGRkZOHnypNjdol/l5eWhb9++hsJ5aWlp+PLLL8XuFlmQk5MDhUKBrKwssbsiaQw/DrJp0yZkZWVh/vz5OHz4MIYPH47Ro0ejvLxc7K4RgJs3byI5ORl///vfxe4KNbFnzx7MnDkT3333HQoKCnDnzh2kp6fj5s2bYneNAHTq1AlLly7FwYMHcfDgQdx///0YO3Ysjh8/LnbXqIkffvgBq1atQt++fcXuiuRxq7uDDBo0CAMGDEBeXp7hWkJCAjIyMpCTkyNiz6gphUKBrVu3IiMjQ+yukBlXr15FeHg49uzZg3vuuUfs7pAZwcHB+Otf/4pp06aJ3RX6VV1dHQYMGIDc3FwsWbIE/fr1w/Lly8XulmRx5McBGhsbUVRUhPT0dKPr6enpKCwsFKlXRO6puroagO4FlqRFo9Fg48aNuHnzJtLS0sTuDt1l5syZeOSRR/Dggw+K3RW3wINNHeDatWvQaDSIiIgwuh4REYGKigqRekXkfgRBQHZ2NoYNG4akpCSxu0O/OnbsGNLS0lBfX4927dph69atSExMFLtb9KuNGzfi0KFD+OGHH8Tuittg+HEghUJhdFsQBJNrRGTZrFmzcPToUezfv1/srtBdevbsieLiYvz888/YsmULJk+ejD179jAAScCFCxfwxz/+EV9//TV8fX3F7o7bYPhxgNDQUCiVSpNRnsrKSpPRICIyb/bs2fj888+xd+9edOrUSezu0F28vb3RrVs3AEBqaip++OEHrFixAitXrhS5Z1RUVITKykqkpKQYrmk0Guzduxd///vf0dDQAKVSKWIPpYlrfhzA29sbKSkpKCgoMLpeUFCAIUOGiNQrIvcgCAJmzZqFzz77DLt27UJ8fLzYXaIWCIKAhoYGsbtBAB544AEcO3YMxcXFho/U1FRMnDgRxcXFDD4WcOTHQbKzs5GZmYnU1FSkpaVh1apVKC8vx/Tp08XuGkG3E+L06dOG22VlZSguLkZwcDA6d+4sYs9o5syZ+PTTT7F9+3YEBgYaRlBVKhX8/PxE7h39+c9/xujRoxETE4Pa2lps3LgRu3fvRn5+vthdIwCBgYEm6+MCAgIQEhLCdXPNYPhxkPHjx6OqqgqLFy+GWq1GUlISdu7cidjYWLG7RgAOHjyI++67z3A7OzsbADB58mSsXbtWpF4RAEN5iBEjRhhdX7NmDaZMmeL6DpGRK1euIDMzE2q1GiqVCn379kV+fj5GjhwpdteIWo11foiIiEhWuOaHiIiIZIXhh4iIiGSF4YeIiIhkheGHiIiIZIXhh4iIiGSF4YeIiIhkheGHiIiIZIXhh4iIiGSF4YeIiIhkheGHiIiIZIXhh4hcZsSIEcjKypLdcxORtPBgUyLyKCNGjEC/fv2wfPlyo+ufffYZ2rZtK0qfsrKycO7cOWzbtk2U5yciYxz5ISJZCA4ORmBgoCjP/cMPP2DgwIGiPDcRmWL4ISJR5OfnY9iwYWjfvj1CQkLw6KOP4syZM0ZtamtrMXHiRAQEBCAqKgp/+9vfmp2+mjJlCvbs2YMVK1ZAoVBAoVDg3LlzAIynvUaMGIHZs2cjKysLHTp0QEREBFatWoWbN2/i2WefRWBgILp27Yovv/zS6PEFQcCyZcvQpUsX+Pn5ITk5GZs3b7b4Nd6+fRve3t4oLCzE/PnzoVAoMGjQoFb/mxGRYzD8EJEobt68iezsbPzwww/4z3/+Ay8vLzz++OPQarWGNtnZ2fjmm2/w+eefo6CgAPv27cOhQ4csPuaKFSuQlpaG559/Hmq1Gmq1GjExMWbb/vOf/0RoaCgOHDiA2bNnY8aMGXjiiScwZMgQHDp0CKNGjUJmZiZu3bpluM+rr76KNWvWIC8vD8ePH8ecOXMwadIk7Nmzx+xzKJVK7N+/HwBQXFwMtVqNr776qjX/XETkQFzzQ0Si+H//7/8Z3f7oo48QHh6O0tJSJCUloba2Fv/85z/x6aef4oEHHgAArFmzBh07drT4mCqVCt7e3vD390dkZGSzz5+cnIxXX30VADBv3jwsXboUoaGheP755wEAr732GvLy8nD06FEMHjwYN2/exLvvvotdu3YhLS0NANClSxfs378fK1euxL333mvyHF5eXrh8+TJCQkKQnJxs/T8OETkVww8RieLMmTNYsGABvvvuO1y7ds0w4lNeXo6kpCScPXsWt2/fNloro1Kp0LNnT4c8f9++fQ3/r1QqERISgj59+hiuRUREAAAqKysBAKWlpaivr8fIkSONHqexsRH9+/e3+DyHDx9m8CGSGIYfIhLFmDFjEBMTgw8++AAdO3aEVqtFUlISGhsbAejW1wCAQqEwup/+ur2a7vxSKBRG1/TPqw9l+v9+8cUXiI6ONrqvj4+PxecpLi5m+CGSGIYfInK5qqoqnDhxAitXrsTw4cMBwLA2Rq9r165o27YtDhw4YFi3U1NTg1OnTpmdYtLz9vaGRqNxeJ8TExPh4+OD8vLyZp+/qWPHjuHxxx93eH+IqPUYfojI5Tp06ICQkBCsWrUKUVFRKC8vxyuvvGLUJjAwEJMnT8ZLL72E4OBghIeHY+HChfDy8jIZDbpbXFwcvv/+e5w7dw7t2rVDcHAwvLzs39sRGBiIuXPnYs6cOdBqtRg2bBhqampQWFiIdu3aYfLkyWbvp9VqcfToUVy+fBkBAQFQqVR294WI7MPdXkTkcl5eXti4cSOKioqQlJSEOXPm4K9//atJu3fffRdpaWl49NFH8eCDD2Lo0KFISEiAr6+vxceeO3culEolEhMTERYWhvLycof1+4033sBrr72GnJwcJCQkYNSoUdixYwfi4+Mt3mfJkiXYtGkToqOjsXjxYof1hYhaTyE4agKdiMjJbt68iejoaLzzzjuYNm2a2N0hIjfFaS8ikqzDhw/jxx9/xMCBA1FdXW0YORk7dqzIPSMid8bwQ0SS9vbbb+PkyZPw9vZGSkoK9u3bh9DQULG7RURujNNeREREJCtc8ExERESywvBDREREssLwQ0RERLLC8ENERESywvBDREREssLwQ0RERLLC8ENERESywvBDREREssLwQ0RERLLC8ENERESy8v8Bftki3OGM438AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt \n",
    "from scipy import stats\n",
    "import pandas as pd \n",
    "from pandas import DataFrame, Series\n",
    "# from PIL import Image\n",
    "from multiprocessing import freeze_support\n",
    "import pims \n",
    "import trackpy as tp \n",
    "\n",
    "R = 0.5*(10**(-6)) #radius of your particle; half a micro meter\n",
    "# mu = 0.0010016 #dynamic viscosity \n",
    "eta = 0.9532*(10**(-3)) #Pa*s\n",
    "T = 295.15 #K aka 22 C\n",
    "k = 1.380649*(10**(-23)) #Boltzmann constant\n",
    "\n",
    "D = (k*T)/(6*np.pi*eta*R)\n",
    "# print(D)\n",
    "\n",
    "\n",
    "@pims.pipeline \n",
    "def gray(image):\n",
    "    return image[:, :, 1]\n",
    "\n",
    "\n",
    "def main():\n",
    "    # frames = pims.open('../uke2/no_contrast/*.tif')\n",
    "    frames = gray(pims.open('Stack2/brownian*.png')) #use pims to open the video\n",
    "    tp.quiet() \n",
    "    f = tp.batch(frames, 11, minmass = 100, threshold  = 10, processes = 'auto') #100 might have also worked \n",
    "    # tp.annotate(f, frames[0]);\n",
    "\n",
    "    t = tp.link(f, 10, memory=3)\n",
    "    t1 = tp.filter_stubs(t, 7)\n",
    "    # t2 = t1[((t1['mass'] > 260) & (t1['size'] < 2.6) & (t1['ecc'] < 0.3))]\n",
    "# # Compare the number of particles in the unfiltered and filtered data.\n",
    "    print('Before:', t['particle'].nunique())\n",
    "    print('After:', t1['particle'].nunique())\n",
    "    # plt.figure()\n",
    "    # tp.mass_size(t1.groupby('particle').mean()); # convenience function -- just plots size vs. mass\n",
    "    # tp.annotate(t1[t1['frame'] == 0], frames[0]);\n",
    "    # plt.figure()\n",
    "    # tp.plot_traj(t1);\n",
    "    d = tp.compute_drift(t1)\n",
    "    # d.plot()\n",
    "    tm = tp.subtract_drift(t1.copy(), d)\n",
    "    # plt.show()\n",
    "    # ax = tp.plot_traj(tm)\n",
    "    # plt.show()\n",
    "    em = tp.emsd(tm, 260/1902., 15)\n",
    "    time = em.index\n",
    "    res = stats.linregress(em.index,em)\n",
    "    print(res.slope, res.stderr)\n",
    "   # microns per pixel = 100/285., frames per second = 24\n",
    "    plt.plot(em.index, em, 'o', label = 'TrackPy MSD') #plot the MSD\n",
    "    plt.plot(em.index, res.intercept + res.slope*em.index, 'r', label='fitted line')\n",
    "    \n",
    "    # # plt.xscale('log')\n",
    "    # # plt.yscale('log')\n",
    "    plt.ylabel(r'$\\langle \\Delta r^2 \\rangle$ [$\\mu$m$^2$]')\n",
    "    plt.xlabel('lag time $t$')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    freeze_support()\n",
    "    main()\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# f = tp.locate(frames[0], 9, invert = True)\n",
    "# tp.annotate(f, frames[0])\n",
    "\n",
    "# fig, ax = plt.subplots()\n",
    "# ax.hist(f['mass'], bins=20)\n",
    "\n",
    "# Optionally, label the axes.\n",
    "# ax.set(xlabel='mass', ylabel='count')\n",
    "# plt.show()\n",
    "\n",
    "\n",
    "# f = tp.batch(frames, 9, minmass = 200, threshold  = 5, processes = 1)#100 might have also worked \n",
    "\n",
    "\n",
    "# tp.annotate(f, frames)\n",
    "\n",
    "# fig, ax = plt.subplots()\n",
    "# ax.hist(f['mass'], bins=20)\n",
    "# ax.set(xlabel='mass', ylabel='count')\n",
    "# plt.show()\n",
    "\n",
    "# a = tp.subpx_bias(f)\n",
    "# plt.show()\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d3cbfe9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}