Generate iPOT structures.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 168,
   "id": "1e572222",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "804a8563",
   "metadata": {},
   "source": [
    "# 2D projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "id": "cb4f0290",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t0) lose ends? True\n",
      "\t1) lose ends? True\n",
      "\t2) lose ends? True\n",
      "\t3) lose ends? True\n",
      "\t4) lose ends? True\n",
      "\t5) lose ends? True\n",
      "\t6) lose ends? True\n",
      "\t7) lose ends? True\n",
      "\t8) lose ends? True\n",
      "\t9) lose ends? True\n",
      "1) n=11 stable? True\n"
     ]
    }
   ],
   "source": [
    "# defining the number of steps (!= number of points you will have in the end)\n",
    "n0 = 20\n",
    "n = n0\n",
    "# maximum number of attempt to find a stable structures\n",
    "max_attempt = 100\n",
    "str_attempt = 0\n",
    "stable = False\n",
    "\n",
    "while not stable and str_attempt < max_attempt:\n",
    "    n = n0 # reset number of steps\n",
    "    str_attempt += 1\n",
    "    # random walk trajectory\n",
    "    traj_xy = np.zeros((n,2))\n",
    "\n",
    "    # compute RW steps\n",
    "    for i in range(1, n):\n",
    "        test_xyz = traj_xy[i - 1]\n",
    "        attempt = 0\n",
    "        # self-avoiding random walk\n",
    "        while np.min(np.linalg.norm(traj_xy-test_xyz, axis=1)) == 0 and attempt<max_attempt:\n",
    "            attempt += 1\n",
    "            val = random.randint(1, 4)\n",
    "            test_xyz = traj_xy[i - 1].copy()\n",
    "\n",
    "            if val == 1:\n",
    "                test_xyz[0] = traj_xy[i - 1, 0] + 1\n",
    "                test_xyz[1] = traj_xy[i - 1, 1] \n",
    "            elif val == 2:\n",
    "                test_xyz[0] = traj_xy[i - 1, 0] \n",
    "                test_xyz[1] = traj_xy[i - 1, 1] + 1\n",
    "            elif val == 3:\n",
    "                test_xyz[0] = traj_xy[i - 1, 0] - 1\n",
    "                test_xyz[1] = traj_xy[i - 1, 1] \n",
    "            else:\n",
    "                test_xyz[0] = traj_xy[i - 1, 0] \n",
    "                test_xyz[1] = traj_xy[i - 1, 1] - 1\n",
    "        traj_xy[i] = test_xyz\n",
    "    # shift to 0 height\n",
    "    traj_xy[:,1] -= np.min(traj_xy[:,1])\n",
    "    # delete overlapping points\n",
    "    traj_xy = np.unique(traj_xy, axis=0)\n",
    "    n = traj_xy.shape[0]\n",
    "    \n",
    "    # prune loose ends\n",
    "    connect = []\n",
    "    connect_ids = []\n",
    "    # loop over points\n",
    "    for i in range(len(traj_xy)):\n",
    "        for j in range(len(traj_xy)):\n",
    "            if i == j: continue # no self link\n",
    "            # check if two points are closer\n",
    "            dist = np.linalg.norm(traj_xy[i]-traj_xy[j])\n",
    "            if np.isclose(dist-1,0):\n",
    "                connect.append([traj_xy[i], traj_xy[j]])\n",
    "                connect_ids.append([i, j])\n",
    "    connect = np.array(connect)\n",
    "    connect_ids = np.array(connect_ids)\n",
    "    # compute nearest neighbours\n",
    "    ids, id_counts = np.unique(connect_ids[:,0], return_counts=True)\n",
    "\n",
    "    # iteratively delete nodes with only 1 neighbour\n",
    "    ii = 0\n",
    "    while np.any(id_counts<=1):\n",
    "        print(f'\\t{ii}) lose ends?', np.any(id_counts<=1))\n",
    "        ii+=1\n",
    "        ids, id_counts = np.unique(connect_ids[:,0], return_counts=True)\n",
    "        # update trajectory\n",
    "        traj_xy = traj_xy[ids[id_counts>1]]\n",
    "        traj_xy[:,1] -= np.min(traj_xy[:,1])\n",
    "        traj_xy = np.unique(traj_xy, axis=0)\n",
    "        n = traj_xy.shape[0]\n",
    "        # recompute connectivity\n",
    "        connect = []\n",
    "        connect_ids = []\n",
    "        for i in range(len(traj_xy)):\n",
    "            for j in range(len(traj_xy)):\n",
    "                if i == j: continue\n",
    "                dist = np.linalg.norm(traj_xy[i]-traj_xy[j])\n",
    "\n",
    "                if np.isclose(dist-1,0):\n",
    "                    connect.append([traj_xy[i], traj_xy[j]])\n",
    "                    connect_ids.append([i, j])\n",
    "        connect = np.array(connect)\n",
    "        connect_ids = np.array(connect_ids)\n",
    "\n",
    "    # center of mass (only of points, while weight is probabl more in lines)\n",
    "    xcm = np.mean(traj_xy, axis=0)\n",
    "    # identify base\n",
    "    base_mask = np.where(traj_xy[:,1] == np.min(traj_xy[:,1]))\n",
    "    base_mask, traj_xy[base_mask]\n",
    "    # compute stability\n",
    "    stable = xcm[0] > np.min(traj_xy[base_mask, 0]) and xcm[0]<np.max(traj_xy[base_mask, 0])\n",
    "    \n",
    "    print(f'{str_attempt}) n={n} stable?', stable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "id": "0253de36",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASoAAAGdCAYAAACl9T7iAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAmCklEQVR4nO3dfXQTZd438G+SJmnRpsJik5amFSjy/lJAauu64FroYvWm66MgtzegIrDe5Ra2uyKsKAe53ayKiruyAtuFHmVZARH0AIKlUBCouLx0RWRZEIFWmyLPAwlvTZrkev6oRCNN20k7zTXw/Zwzh2RyTeY3XMm3M1dmEp0QQoCISGL6aBdARNQUBhURSY9BRUTSY1ARkfQYVEQkPQYVEUmPQUVE0mNQEZH0YqJdQHMEAgF88803iI+Ph06ni3Y5RNQKhBA4f/48kpOTodc3vs+kiaD65ptvYLfbo10GEamgsrISKSkpjbbRRFDFx8cDqN8gi8US5WqIqDW43W7Y7fbg+7sxmgiqK4d7FouFQUV0jWnOcA4H04lIegwqIpIeg4qIpMegIiLpMaiISHoMKiKSHoOKiKTHoCIi6TGoiEh6ioLqzTffRL9+/YJniGdlZeHDDz9sdJnVq1ejR48eiI2NRd++fbFx48YWFUxE1x9FQZWSkoI//OEP2LdvH/bu3Yuf//znGDVqFA4dOtRg+927d2Ps2LGYOHEiDhw4gPz8fOTn5+Pzzz9vleIbUlVVhW3btqGqqkq1dZB62H/UINFC7du3F0VFRQ0+Nnr0aJGXlxcyLzMzU0yZMkXROlwulwAgXC5Xo+2KioqEXq8XAIRerw9bF8mJ/Xd9ae77WgghIh6j8vv9eOedd3Dx4kVkZWU12Ka8vBw5OTkh83Jzc1FeXt7oc3s8Hrjd7pCpKVVVVZg8eTICgQCA+u+wmjJlCv8yawT7jxqj+NsTDh48iKysLNTW1uLGG2/E2rVr0atXrwbbOp1OWK3WkHlWqxVOp7PRdTgcDsydO1dRXUePHg2+yK/w+/3o168fTCaToueituf1ehvsv2PHjjX5XUV07VMcVN27d0dFRQVcLhfeffddTJgwAdu3bw8bVpGYNWsWCgsLg/evfG9NY7p16wa9Xh/yYjcYDPjss8/4QteAqqoqpKWlXdV/6enpUayKZKH40M9kMiE9PR2DBg2Cw+FA//798frrrzfY1mazoaamJmReTU0NbDZbo+swm83BTxab+x1UKSkpWLJkCQwGA4D6F/nixYsZUhrB/qPGtPg8qkAgAI/H0+BjWVlZKC0tDZlXUlISdkyrpSZOnIgTJ06gffv2OHHiBCZOnKjKekgd7D8KR9Gh36xZszBy5Eikpqbi/PnzWLFiBcrKyrB582YAwPjx49GpUyc4HA4AwLRp0zB06FC88soryMvLwzvvvIO9e/diyZIlrb8l30lJSYHJZOJfYo1i/1FDFAXV6dOnMX78eFRXVyMhIQH9+vXD5s2bMXz4cADAqVOnQn5NIjs7GytWrMDs2bPxu9/9Dt26dcO6devQp0+f1t0KIrqm6YQQItpFNMXtdiMhIQEul6tZ41U2m63JTxZJXuy/64OS9zWv9SMi6TGoiEh6DCoikh6Dioikx6AiIukxqIhIegwqIpIeg4qIpMegIiLpMaiISHoMKiKSHoOKiKTHoCIi6TGoiEh6DCoikh6Dioikx6AiIukxqIhIegwqIpIeg4qIpMegIiLpMaiISHoMKiKSHoOKiKTHoCIi6TGoiEh6DCoikh6Dioikx6AiIukxqIhIegwqIpIeg4qIpMegIiLpMaiISHoMKiKSHoOKiKTHoCIi6TGoiEh6DCoikh6Dioikx6AiIukxqIhIegwqIpIeg4qIpMegIiLpKQoqh8OB2267DfHx8UhMTER+fj6OHDnS6DLFxcXQ6XQhU2xsbIuKJqLri6Kg2r59OwoKCvDJJ5+gpKQEdXV1GDFiBC5evNjochaLBdXV1cHp5MmTLSqaiK4vMUoab9q0KeR+cXExEhMTsW/fPvzsZz8Lu5xOp4PNZousQiK67rVojMrlcgEAOnTo0Gi7CxcuIC0tDXa7HaNGjcKhQ4cabe/xeOB2u0MmIrp+RRxUgUAA06dPxx133IE+ffqEbde9e3csXboU77//PpYvX45AIIDs7GxUVVWFXcbhcCAhISE42e32SMskomuATgghIlnwiSeewIcffoidO3ciJSWl2cvV1dWhZ8+eGDt2LObNm9dgG4/HA4/HE7zvdrtht9vhcrlgsViaXIfNZoPT6Wx2TSQX9t/1we12IyEhoVnva0VjVFdMnToV69evx44dOxSFFAAYjUZkZGTg2LFjYduYzWaYzeZISiOia5CiQz8hBKZOnYq1a9di69at6Ny5s+IV+v1+HDx4EElJSYqXJaLrk6I9qoKCAqxYsQLvv/8+4uPjg7vnCQkJiIuLAwCMHz8enTp1gsPhAAA8//zzuP3225Geno5z587h5ZdfxsmTJ/H444+38qYQ0bVKUVC9+eabAIBhw4aFzF+2bBkeeeQRAMCpU6eg13+/o3b27FlMmjQJTqcT7du3x6BBg7B792706tWrZZUT0XUj4sH0tqRk0A3gYKzWsf+uD0re17zWj4ikx6AiIukxqIhIegwqIpIeg4qIpMegIiLpMaiISHoMKiKSHoOKiKTHoCIi6TGoiEh6DCoikh6Dioikx6AiIukxqIhIegwqIpIeg4qIpMegIiLpMaiISHoMKiKSHoOKiKTHoCIi6TGoiEh6DCoikh6Dioikx6AiIukxqIhIegwqIpIeg4qIpMegIiLpMaiISHoMKiKSHoOKiKTHoCIi6TGoiEh6DCoikh6Dioikx6AiIukxqIhIegwqIpIeg4qIpMegIiLpMaiISHoMKiKSHoOKiNQjBODztPhpFAWVw+HAbbfdhvj4eCQmJiI/Px9HjhxpcrnVq1ejR48eiI2NRd++fbFx48aICyYiDRACOLYF+MtdwGt9AFdVi55OUVBt374dBQUF+OSTT1BSUoK6ujqMGDECFy9eDLvM7t27MXbsWEycOBEHDhxAfn4+8vPz8fnnn7eo8HCqqqrg9XpRVdWy/xiKDvaftlVVVuKfa16Fd+FPgeX/B/jmn8DF08DFMy17YtECp0+fFgDE9u3bw7YZPXq0yMvLC5mXmZkppkyZ0uz1uFwuAUC4XK5G2xUVFQm9Xi8ACL1eL4qKipq9Doo+9p+GBQLiwz/9Vnz6+A1CzLGIumfjhZhj+X76+sBVizT3fS2EEC0ao3K5XACADh06hG1TXl6OnJyckHm5ubkoLy8Pu4zH44Hb7Q6ZmlJVVYXJkycjEAgAAAKBAKZMmcK/zBrB/tOwL7fBu/Cn+MWZJchIqo+UGL2uVVcRE+mCgUAA06dPxx133IE+ffqEbed0OmG1WkPmWa1WOJ3OsMs4HA7MnTtXUT1Hjx4Nvsiv8Pv96NevH0wmk6Lnorbn9XrZfxq1Y3Qtbu0gALR+QF0RcVAVFBTg888/x86dO1uzHgDArFmzUFhYGLzvdrtht9sbXaZbt27Q6/UhL3aDwYDPPvsMKSkprV4jta6qqiqkpaWx/7Toy23wfjgbpjOfwxcQqoRVRId+U6dOxfr167Ft27YmX0Q2mw01NTUh82pqamCz2cIuYzabYbFYQqampKSkYMmSJTAYDADqX+SLFy/mi1wj2H8a1vUumAp2YlPHyTjgrP9D4wuIVl2FoqASQmDq1KlYu3Yttm7dis6dOze5TFZWFkpLS0PmlZSUICsrS1mlzTBx4kScOHEC7du3x4kTJzBx4sRWXweph/2nYTodfjH1ZSQ9exgPrNMjkNj3uwda6VRNJQP7TzzxhEhISBBlZWWiuro6OF26dCnYZty4cWLmzJnB+7t27RIxMTFi/vz54vDhw2LOnDnCaDSKgwcPNnu9Sj4dEEIIq9Xa/I0i6bD/tM1qtQoRCAhxtESIxUOFeDldiHNVV7VT8r5WFFQAGpyWLVsWbDN06FAxYcKEkOVWrVolbr31VmEymUTv3r3Fhg0blKyWQXWdYf9pW0j/BQJC1NU22E7J+1rRYLoQTR93lpWVXTXvwQcfxIMPPqhkVUR0LdDpgBhzi5+G1/oRkfQYVEQkPQYVEUmPQUVE0mNQEZH0GFREJD0GFRFJj0FFRNJjUBGR9BhURCQ9BhURSY9BRUTSY1ARkfQYVEQkPQYVEUmPQUVE0mNQEZH0GFREJD0GFRFJj0FFRNJjUBGR9BhURCQ9BhURSY9BRUTSY1ARkfQYVEQkPQYVEUmPQUVE0mNQEZH0GFREJD0GFRFJj0FFRNJjUBGR9BhURCQ9BhURSY9BRUTSY1ARkfQYVEQkPQYVEUmPQUVE0mNQEZH0GFREJD0GFRFJj0FFRNJjUBGR9BQH1Y4dO3DfffchOTkZOp0O69ata7R9WVkZdDrdVZPT6Yy0ZiK6zigOqosXL6J///5YuHChouWOHDmC6urq4JSYmKh01UR0nYpRusDIkSMxcuRIxStKTEzETTfdpHg5IqI2G6MaMGAAkpKSMHz4cOzatavRth6PB263O2QiouuX6kGVlJSERYsWYc2aNVizZg3sdjuGDRuG/fv3h13G4XAgISEhONntdrXLJCKJ6YQQIuKFdTqsXbsW+fn5ipYbOnQoUlNT8fbbbzf4uMfjgcfjCd53u92w2+1wuVywWCxNPr/NZuNgvYax/7Stuf3ndruRkJDQrPe14jGq1jBkyBDs3Lkz7ONmsxlms7kNKyIimUXlPKqKigokJSVFY9VEpEGK96guXLiAY8eOBe9/9dVXqKioQIcOHZCamopZs2bh66+/xltvvQUAWLBgATp37ozevXujtrYWRUVF2Lp1Kz766KPW2woiuqYpDqq9e/firrvuCt4vLCwEAEyYMAHFxcWorq7GqVOngo97vV785je/wddff4127dqhX79+2LJlS8hzEBE1pkWD6W1FyaAbwMFYrWP/aZsag+m81o+IpMegIiLpMaiISHoMKiKSXlRO+KTrxKlTwJkz4R/v2BFITW27ekizGFSkjlOngO7dgdra8G1iY4EjRxhW1CQe+pE6zpxpPKSA+scb2+Mi+g6Dioikx6AiIukxqIhIegwqIpIeg4qIpMegIiLpMahIHR071p8n1ZjY2Pp2RE3gCZ+kjtTU+pM5eWY6tQIGFaknNRVITYXX78Xc8rkAgDlZc2AymKJcGGkND/1Idb6ADx98+QE++PID+AK+aJdDGsQ9KlKdUW9E4aDC4G0ipRhUpDqjwYhH+zwa7TJIw3joR0TS4x4VqS4gAvj20rcAgJvb3Qy9jn8fSRkGFamu1leLnHdzAAB7/nMP2hnbRbki0hoGFbWJGB1fahQ5vnpIde2M7XBg/IFol0EaxsECIpIeg4qIpMdDP1Kd1+/FS/94CQAw47YZvISGFOMeFanOF/Bh5ZGVWHlkJS+hoYhwj4pUZ9Qb8UT/J4K3iZRiUJHqjAYj/nvAf0e7DNIwHvoRkfS4R0WqE0LgfN15AEC8MR46nS7KFZHWMKhIdZd9l3HH3+8AwEtoKDI89CMi6XGPilQXFxOH/eP2A+A1fxQZvmpIdTqdDkYdT0ugyPHQj4ikxz0qUl2dvw5/PPBHAMCTGU/CaODeFSnDPSpSXV2gDsWHilF8qBh1gbpol0MaxD0qUp1Rb8QjvR8J3iZSikFFqjMajPjN4N9EuwzSMB76EZH0uEdFqhNCwCfqv94lRhfDS2hIMe5Rkeou+y5j4NsDMfDtgbjsuxztckiDGFREJD3FQbVjxw7cd999SE5Ohk6nw7p165pcpqysDAMHDoTZbEZ6ejqKi4sjKJW0Ki4mDrvG7sKusbsQFxMX7XJIgxQH1cWLF9G/f38sXLiwWe2/+uor5OXl4a677kJFRQWmT5+Oxx9/HJs3b1ZcbHNUVVXB6/WiqqpKlecn5XQ6HSwmCywmS5PjU+w/bVOt/0QLABBr165ttM2MGTNE7969Q+aNGTNG5ObmNns9LpdLABAul6vRdkVFRUKv1wsAQq/Xi6Kiomavg6KP/adtSvuvue9rIYRQPajuvPNOMW3atJB5S5cuFRaLJewytbW1wuVyBafKysomN6iysjL4n3RlMhgMorKyUskmkQq8Pq9YeGChWHhgofD6vA22Yf9pWyT9pySoVD89wel0wmq1hsyzWq1wu924fPky4uKuHrNwOByYO3euovUcPXoUgUAgZJ7f70e/fv1gMvHnmaLKCNz8vzcDAJ6/73mggatovF4v+0/DwvXfsWPHkJKS0uLnl/I8qlmzZqGwsDB43+12w263N7pMt27doNfrQ/6zDAYDPvvss1b5j6LIhfyu34mGf9evqqoKaWlp7D+NCtd/6enprfL8qp+eYLPZUFNTEzKvpqYGFoulwb0pADCbzbBYLCFTU1JSUrBkyRIYDAYA9f9Jixcv5otcAiaDCbNvn43Zt88O++Oj7D9tU73/WnJcimYOpvfp0ydk3tixY1UZTBei/li5ffv2HNvQKPaftinpPyXva8V7VBcuXEBFRQUqKioA1J9+UFFRgVOnTgGoP2wbP358sP2vfvUrHD9+HDNmzMC//vUv/PnPf8aqVavw61//uhVi9mopKSkwmUz8S6xR7D9tU6v/FAfV3r17kZGRgYyMDABAYWEhMjIy8NxzzwEAqqurg6EFAJ07d8aGDRtQUlKC/v3745VXXkFRURFyc3NbaRNIdpfqLiHjrQxkvJWBS3WXol0OaZDiwfRhw4ZBCBH28YbOOh82bBgOHDigdFV0DblyUTJRJKT81I+uLbExsdjywJbgbSKlGFSkOr1OD+sN1qYbEoXBb08gIulxj4pUV+evw/LDywEA/9Xzv/grNKQYg4pUVxeow6v7XgUAjOk+hkFFijGoSHUx+hj8R9f/CN4mUoqvGlKdyWDCCz99IdplkIZxMJ2IpMegIiLpMahIdZfqLiF7RTayV2TzEhqKCMeoqE2crzsf7RJIwxhUpLrYmFis/+X64G0ipRhUpDq9To80S1q0yyAN4xgVEUmPe1SkurpAHd7997sAgAdufQBGPc9MJ2UYVKS6On8dfr/n9wCAUV1HMahIMQYVqc6gN2B42vDgbSKlGFSkOrPBjFeHvRrtMkjDOJhORNJjUBGR9BhUpLrLvsu4e9XduHvV3bjsuxztckiDOEZFqhNC4PTl08HbREoxqEh1ZoMZq+9bHbxNpBSDilRn0BvQo0OPaJdBGsYxKiKSHveoSHV1gTpsOL4BAJDXJY9nppNiDCpSXZ2/Ds/uehYAMCJtBIOKFGNQkeoMegPu7HRn8DaRUgwqUp3ZYMafc/4c7TJIwziYTkTSY1ARkfQYVKS6y77LyHsvD3nv5fESGooIx6hIdUIInDp/KnibSCkGFanObDDjrZFvBW8TKcWgItUZ9AZkJGZEuwzSMI5REZH0uEdFqvMFfCg9VQoAuDv1bsTo+bIjZfiKIdV5/V78dvtvAQB7/nMPg4oU4yuGVKfX6THYOjh4m0gpBhWpLjYmFst+sSzaZZCG8c8bEUmPQUVE0mNQkepqfbV44IMH8MAHD6DWVxvtckiDOEZFqguIAI6cPRK8TaQUg4pUZzaYsXj44uBtIqUiOvRbuHAhbrnlFsTGxiIzMxOffvpp2LbFxcXQ6XQhU2xsbMQFk/YY9AZkJ2cjOzmb3/BJEVEcVCtXrkRhYSHmzJmD/fv3o3///sjNzcXp06fDLmOxWFBdXR2cTp482aKiiej6ojioXn31VUyaNAmPPvooevXqhUWLFqFdu3ZYunRp2GV0Oh1sNltwslqtLSqatMUX8GFH1Q7sqNoBX8AX7XJIgxQFldfrxb59+5CTk/P9E+j1yMnJQXl5edjlLly4gLS0NNjtdowaNQqHDh1qdD0ejwdutztkIu3y+r0oKC1AQWkBvH5vtMshDVIUVGfOnIHf779qj8hqtcLpdDa4TPfu3bF06VK8//77WL58OQKBALKzs1FVVRV2PQ6HAwkJCcHJbrcrKZMko9fp0fsnvdH7J715CQ1FRPVP/bKyspCVlRW8n52djZ49e2Lx4sWYN29eg8vMmjULhYWFwftut5thpWGxMbF45953ol0GaZiioOrYsSMMBgNqampC5tfU1MBmszXrOYxGIzIyMnDs2LGwbcxmM8xmfoxNRPUU7YebTCYMGjQIpaWlwXmBQAClpaUhe02N8fv9OHjwIJKSkpRVSkTXLcWHfoWFhZgwYQIGDx6MIUOGYMGCBbh48SIeffRRAMD48ePRqVMnOBwOAMDzzz+P22+/Henp6Th37hxefvllnDx5Eo8//njrbglJq9ZXi0kfTQIA/GXEXxAbw/PoSBnFQTVmzBh8++23eO655+B0OjFgwABs2rQpOMB+6tQp6PXf76idPXsWkyZNgtPpRPv27TFo0CDs3r0bvXr1ar2tIKkFRAAV31YEbxMppRMa+P0it9uNhIQEuFwuWCyWJtvbbLawn0JS2/MFfNhetR0AMDRlaJPf8Mn+07bm9p+S9zWv9SPVxehjcHfq3dEugzSMJ7UQkfS4R0Wq8wf82H96PwBgYOJAXphMijGoSHUevwePbX4MQP2v0LTTt4tyRaQ1DCpSnU6nQ9eErsHbREoxqEh1cTFxWJe/LtplkIZxMJ2IpMegIiLpMahIdVcuoZn00ST+Cg1FhGNUpLqACOCT6k+Ct4mUYlCR6kwGExx3OoK3iZRiUJHqYvQxuLfLvdEugzSMY1REJD3uUZHq/AE/Dv+/wwCAnh168hIaUoxBRarz+D0Yu2EsAF5CQ5FhUJHqdDodkm9IDt4mUopBRaqLi4nD5gc2R7sM0jAOphOR9BhURCQ9BhWpzuP34MmtT+LJrU/C4/dEuxzSII5Rker8AT+2VW4L3gbPTiCFGFSkOqPBiDlZc4K3iZRiUJHqjHojHrj1gWiXQRrGMSoikh73qEh1ARHA8XPHAQBdbuoCvY5/H0kZBhWprtZXi19+8EsA311CY+QlNKQMg4raRHtz+2iXQBrGoCLVtTO2w46HdkS7DNIwDhYQkfQYVEQkPQYVqc7j9+DpHU/j6R1P8xIaigiDilTnD/ix8auN2PjVxvpLaIgU4mA6qc5oMGLGbTOCt4mUYlCR6ox6I8b1GhftMkjDeOhHRNLjHhWpLiACqL5YDQBIuiGJl9CQYgwqUl2trxa/WPMLALyEhiLDoKI2ERcTF+0SSMMYVKS6dsZ2+PThT6NdBmkYBwuISHoMKiKSHg/9SHVevxe/3/N7AMDvMn8Hk8EU5YpIa7hHRarzBXxYc3QN1hxdA1/AF+1ySIO4R0WqM+qN+J+M/wneJlLqmtij8vj8EELA7/ejrKwMly9fRllZGfx+XgArA6PBiMn9JmNyv8mNXuvH/tM2VftPROCNN94QaWlpwmw2iyFDhog9e/Y02n7VqlWie/fuwmw2iz59+ogNGzYoWp/L5RIAhMvluuqxr89eEoPmfSQGT/xfcbMtWQAITikpKWLNmjWK1kUqKSkRomfP+n8bsGbNGpGSksL+06hI+q+x9/WP6YQQQkmwrVy5EuPHj8eiRYuQmZmJBQsWYPXq1Thy5AgSExOvar9792787Gc/g8PhwL333osVK1bgxRdfxP79+9GnT59mrdPtdiMhIQEulwsWiyXksc+/duHnU1/Ct+t+H3b5mJgY6PXXxM6jNgmBXQEfBvsF9up0uCMmBtDpgg8HAgH4fOHHrth/cmuq/9asWYP777//qvmNva+vojQ5hwwZIgoKCoL3/X6/SE5OFg6Ho8H2o0ePFnl5eSHzMjMzxZQpU5q9zsaSt+Lk/xWG+I4hSc5Jrik3BkLg+2mEBDVxarvJbrcLn8+n6H39Y4r+THm9Xuzbtw85OTnBeXq9Hjk5OSgvL29wmfLy8pD2AJCbmxu2PQB4PB643e6QKZz9e3bDf/6Mks2gNva8H/B9twPlAzAvqtVQW6usrMTHH3/coudQFFRnzpyB3++H1WoNmW+1WuF0Ohtcxul0KmoPAA6HAwkJCcHJbreHbfvt6RoFW0BtbQSAIQKIEfX3YwAM+W4+XT+qq6tbtLyUB/6zZs2Cy+UKTpWVlWHb3pxoDfsYRd881O9F/RD3qq4/SUlJLVpe0XlUHTt2hMFgQE1N6F5MTU0NbDZbg8vYbDZF7QHAbDbDbDY3q6aBmdkwxHfk4Z+ERqB+7+nHfrhX9VGbVkTRYLfbceedd7boORQFlclkwqBBg1BaWor8/HwA9SP+paWlmDp1aoPLZGVlobS0FNOnTw/OKykpQVZWVsRF/5DBYECHuyfzUz/ZCIEXfD74hYChgYf9AF7Q6VAWE4OAEPzUT8Oa+tRvwYIFMBgaehUo0OyP3r7zzjvvCLPZLIqLi8UXX3whJk+eLG666SbhdDqFEEKMGzdOzJw5M9h+165dIiYmRsyfP18cPnxYzJkzRxiNRnHw4MFmr7OxTwe+OXdJDJ5X0uB5VHa7nefhREttrRBWa8infVdNNlt9O9HweTjsP+2IpP9UPY8KAN544w28/PLLcDqdGDBgAP74xz8iMzMTADBs2DDccsstKC4uDrZfvXo1Zs+ejRMnTqBbt2546aWXcM899zR7fU2db+Hx+WEy6BEIBPDxxx+juroaSUlJuPPOO1ue5BS5ykrg22/DP56YCKSkBO/6/X72n4Yp7T8l51FFFFRtTdGJYUSkCUre1zzwJyLpMaiISHoMKiKSHoOKiKTHoCIi6TGoiEh6DCoikh6Dioikx6AiIulp4ldorpw839gX6BGRtlx5Pzfn4hhNBNX58+cBoNEv0CMibTp//jwSEhIabaOJa/0CgQC++eYbxMfHQ/eDHwVoiNvtht1uR2Vl5TV5XSC3T9u4fd8TQuD8+fNITk5u8mt8NLFHpdfrkfKDq+ybw2KxXJMvhCu4fdrG7avX1J7UFRxMJyLpMaiISHrXXFCZzWbMmTOn2d+5rjXcPm3j9kVGE4PpRHR9u+b2qIjo2sOgIiLpMaiISHoMKiKSnqaC6pZbboFOp7tqKigoCLvM6tWr0aNHD8TGxqJv377YuHFjG1asjNLtKy4uvqptbGxsG1fdfH6/H88++yw6d+6MuLg4dO3aFfPmzWvyWq+ysjIMHDgQZrMZ6enpIT/FJpNItq+srKzBPnc6nW1YefOdP38e06dPR1paGuLi4pCdnY1//OMfjS7TKv3Xwt8dbFOnT58W1dXVwamkpEQAENu2bWuw/a5du4TBYBAvvfSS+OKLL8Ts2bMV//hpW1K6fcuWLRMWiyVkmSs/BCujF154QfzkJz8R69evF1999ZVYvXq1uPHGG8Xrr78edpnjx4+Ldu3aicLCQvHFF1+IP/3pT8JgMIhNmza1YeXNE8n2bdu2TQAQR44cCelHv9/fhpU33+jRo0WvXr3E9u3bxdGjR8WcOXOExWIRVVVVDbZvrf7TVFD92LRp00TXrl1FIBBo8PHRo0eLvLy8kHmZmZliypQpbVFeizW1fcuWLRMJCQltW1QL5OXlicceeyxk3v333y8efvjhsMvMmDFD9O7dO2TemDFjRG5urio1tkQk23clqM6ePatydS136dIlYTAYxPr160PmDxw4UDzzzDMNLtNa/aepQ78f8nq9WL58OR577LGwFyqXl5cjJycnZF5ubi7Ky8vbosQWac72AcCFCxeQlpYGu92OUaNG4dChQ21YpTLZ2dkoLS3Fv//9bwDAP//5T+zcuRMjR44Mu4yW+jCS7btiwIABSEpKwvDhw7Fr1y61S42Iz+eD3++/anghLi4OO3fubHCZ1uo/TVyU3JB169bh3LlzeOSRR8K2cTqdsFqtIfOsVqu0x/8/1Jzt6969O5YuXYp+/frB5XJh/vz5yM7OxqFDhxRfxN0WZs6cCbfbjR49esBgMMDv9+OFF17Aww8/HHaZcH3odrtx+fJlxMXFqV12s0WyfUlJSVi0aBEGDx4Mj8eDoqIiDBs2DHv27MHAgQPbsPqmxcfHIysrC/PmzUPPnj1htVrx97//HeXl5UhPT29wmdbqP80G1V//+leMHDkSycnJ0S5FFc3ZvqysLGRlZQXvZ2dno2fPnli8eDHmzZvXFmUqsmrVKvztb3/DihUr0Lt3b1RUVGD69OlITk7GhAkTol1ei0Wyfd27d0f37t2D97Ozs/Hll1/itddew9tvv91WpTfb22+/jcceewydOnWCwWDAwIEDMXbsWOzbt0/V9WoyqE6ePIktW7bgvffea7SdzWZDTU1NyLyamhrYbDY1y2ux5m7fjxmNRmRkZODYsWMqVdYyTz31FGbOnImHHnoIANC3b1+cPHkSDocj7Bs5XB9aLBap9qaAyLavIUOGDAl7KBVtXbt2xfbt23Hx4kW43W4kJSVhzJgx6NKlS4PtW6v/NDlGtWzZMiQmJiIvL6/RdllZWSgtLQ2ZV1JSErIXIqPmbt+P+f1+HDx4EElJSSpV1jKXLl266gvSDAYDAoFA2GW01IeRbF9DKioqpO3DK2644QYkJSXh7Nmz2Lx5M0aNGtVgu1brP2Xj/tHn9/tFamqqePrpp696bNy4cWLmzJnB+7t27RIxMTFi/vz54vDhw2LOnDlSn54ghLLtmzt3rti8ebP48ssvxb59+8RDDz0kYmNjxaFDh9qy5GabMGGC6NSpU/Dj+/fee0907NhRzJgxI9hm5syZYty4ccH7Vz7efuqpp8Thw4fFwoULpT09IZLte+2118S6devE0aNHxcGDB8W0adOEXq8XW7ZsicYmNGnTpk3iww8/FMePHxcfffSR6N+/v8jMzBRer1cIoV7/aS6oNm/eHDzv5MeGDh0qJkyYEDJv1apV4tZbbxUmk0n07t1bbNiwoY0qjYyS7Zs+fbpITU0VJpNJWK1Wcc8994j9+/e3YbXKuN1uMW3aNJGamipiY2NFly5dxDPPPCM8Hk+wzYQJE8TQoUNDltu2bZsYMGCAMJlMokuXLmLZsmVtW3gzRbJ9L774oujatauIjY0VHTp0EMOGDRNbt26NQvXNs3LlStGlSxdhMpmEzWYTBQUF4ty5c8HH1eo/fs0LEUlPk2NURHR9YVARkfQYVEQkPQYVEUmPQUVE0mNQEZH0GFREJD0GFRFJj0FFRNJjUBGR9BhURCQ9BhURSe//AxD2J6PnDWsEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# interactive (works bad)\n",
    "# %matplotlib notebook \n",
    "# static\n",
    "%matplotlib inline\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "for line in connect:\n",
    "    ax.plot(line[:,0], line[:,1], '-k', lw=0.5)\n",
    "    \n",
    "ax.plot(traj_xy[:,0], traj_xy[:,1], '.k')\n",
    "ax.plot(traj_xy[np.argmin(traj_xy[:,1]),0], traj_xy[np.argmin(traj_xy[:,1]),1], '>')\n",
    "ax.plot(traj_xy[-1,0], traj_xy[-1,1], '<')\n",
    "\n",
    "\n",
    "ax.plot(traj_xy[base_mask, 0], traj_xy[base_mask, 1], 'ok', lw=5)\n",
    "ax.plot([np.min(traj_xy[base_mask, 0]), np.max(traj_xy[base_mask, 0])], \n",
    "        2*[np.min(traj_xy[base_mask, 1])], '-k', lw=5, zorder=-1)\n",
    "\n",
    "ax.plot([xcm[0], xcm[0]], [0, xcm[1]], ':')\n",
    "ax.plot(xcm[0], xcm[1], 'sr')\n",
    "ax.plot(xcm[0], 0, '^r')\n",
    "\n",
    "ax.set_aspect('equal')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae5c79e6",
   "metadata": {},
   "source": [
    "### make 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "id": "594fbee1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(82, 2, 3)\n"
     ]
    }
   ],
   "source": [
    "# make 3D by addin a layer behind\n",
    "\n",
    "N = traj_xy.shape[0]*2\n",
    "traj_xyz = np.zeros((N, 3))\n",
    "\n",
    "traj_xyz[:N//2,0] = traj_xy[:,0]\n",
    "traj_xyz[N//2:,0] = traj_xy[:,0]\n",
    "\n",
    "traj_xyz[:N//2,2] = traj_xy[:,1]\n",
    "traj_xyz[N//2:,2] = traj_xy[:,1]\n",
    "\n",
    "traj_xyz[:N//2,1] = 0\n",
    "traj_xyz[N//2:,1] = 1\n",
    "\n",
    "connect = []\n",
    "connect_ids = []\n",
    "for i in range(len(traj_xyz)):\n",
    "    for j in range(len(traj_xyz)):\n",
    "        if i == j: continue\n",
    "        dist = np.linalg.norm(traj_xyz[i]-traj_xyz[j])\n",
    "        if np.isclose(dist-1,0):\n",
    "            connect.append([traj_xyz[i], traj_xyz[j]])\n",
    "            connect_ids.append([i, j])\n",
    "connect = np.array(connect)\n",
    "connect_ids = np.array(connect_ids)\n",
    "\n",
    "xcm = np.mean(traj_xyz, axis=0)\n",
    "base_mask = np.where(traj_xyz[:,2] == np.min(traj_xyz[:,2]))[0]\n",
    "body_mask = [i for i in range(traj_xyz.shape[0]) if i not in base_mask]\n",
    "\n",
    "print(connect.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "c77176d8",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAGbCAYAAAAr/4yjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABc8UlEQVR4nO3dd3gU1f4G8Hc3ZTe9kBASUkggBAi9XkpI6NKlC6iAWFFsXBtgQVERG15QBEVEBSGErqCAgKEJgqJSEkISAqGE9GR7O78/8svIkrbpAd7P8/jcy+zM7Nlssu/OnHO+RyaEECAiIgIgr+8GEBFRw8FQICIiCUOBiIgkDAUiIpIwFIiISMJQICIiCUOBiIgkDAUiIpIwFIiISHJHhcL06dPRrFmz+m5GnYuJiUFMTIz07wMHDkAmkyEuLq7+GlVPZDIZ3njjjfpuRp1r1qwZpk+fXt/NoFu88cYbkMlk9d2MSqlSKHz99deQyWTSf/b29mjatCmmT5+OK1eu1HQbb0uxsbGQyWTYsmVLicc6dOgAmUyG/fv3l3gsODgYvXr1qosmlury5ctYsGABunfvDi8vL/j4+CAmJgZ79+4tse8vv/yChx56CC1btoSzszPCwsLw8MMP49q1a/XQ8up75513sHXr1vpuRq3S6XR499130aZNGzg7O6Np06aYMGECzpw5Y7Xfnfbe1iaNRoM33ngDBw4cqO+m1Aj76hz85ptvIjQ0FDqdDr/99hu+/vprHDp0CKdPn4ZSqaypNt6W+vTpAwA4dOgQxowZI20vKCjA6dOnYW9vj8OHD6Nfv37SY5cvX8bly5dx33331Xl7i23btg3vvfce7r33XkybNg0mkwnffPMNBg0ahK+++gozZsyQ9n3ppZeQk5ODCRMmIDw8HCkpKVi2bBl++OEHnDp1Ck2aNKnz9mu1WtjbV+3X+p133sH48eNx77331myjGpCpU6di+/bteOSRR9C5c2dcvXoVn376KXr27Il//vkHISEhABrme9tQaTQaLFiwAACsrtgBYP78+Xj55ZfroVXVIKpg9erVAoD4/fffrba/9NJLAoDYsGFDVU5bbdOmTRMhISH18tylCQ0NFd27d7fa9tNPPwmZTCYmT54shgwZYvXYunXrBACxbdu2Sj1PdHS0iI6Olv69f/9+AUBs3Lix0m0+ffq0yMzMtNqm0+lEq1atRGBgoNX2X3/9VZjN5hLbAIh58+ZV+rnrm4uLi5g2bVp9N6NKQkJCKmx7enq6ACD++9//Wm3ft2+fACA++ugjadud9t7WpszMTAFAvP766/XdlBpRo30KUVFRAIDk5GRpm8FgwGuvvYYuXbrAw8MDLi4uiIqKKnHr5OLFi5DJZPjggw+wcuVKNG/eHAqFAt26dcPvv/9e4rm2bt2Ktm3bQqlUom3btqXepgEAtVqNOXPmICgoCAqFAhEREfjggw8gbikOK5PJ8NRTT2Hjxo1o06YNnJycpG9PALBixQq0aNECSqUSMTExuHjxYoU/jz59+uDPP/+EVquVth0+fBiRkZEYOnQofvvtN1gsFqvHZDIZevfuDQBYvXo1+vfvj8aNG0OhUKBNmzZYvnx5hc9bGr1ejxEjRsDDwwNHjhwpc7/IyEj4+PhYbVMoFBg2bBjS09NRWFgobe/bty/kcutfob59+8Lb2xvnzp2rsE0xMTFo27YtTp48iV69esHJyQmhoaH4/PPPS+x748YNzJw5E35+flAqlejQoQPWrFlTYr9b+xSK7+leuHAB06dPh6enJzw8PDBjxgxoNBqr49RqNdasWSPdFi2+R19YWIhnn30WzZo1g0KhQOPGjTFo0CD88ccf5b6+tLQ0zJo1CxEREXByckKjRo0wYcKEEr87xbdjDx8+jOeffx6+vr5wcXHBmDFjkJmZabWvEAILFy5EYGAgnJ2d0a9fvxK3fspS/N75+flZbff39wcAODk5Sduq+94CwLFjxzBs2DB4eXnBxcUF7du3xyeffGK1z759+xAVFQUXFxd4enpi9OjRJc5v63sI/Pt3XPz5oFAoEBkZiZ9++qlE+65cuYKHHnoIfn5+0n5fffVVif10Oh3eeOMNtGzZEkqlEv7+/hg7diySk5Nx8eJF+Pr6AgAWLFgg/e4U/w6W1qdgMpnw1ltvSZ9xzZo1w9y5c6HX6632a9asGUaMGIFDhw6he/fuUCqVCAsLwzfffGO1n9FoxIIFCxAeHg6lUolGjRqhT58+2LNnTznvTjmqkiRlXSksW7ZMABDLly+XtmVmZgp/f3/x/PPPi+XLl4vFixeLiIgI4eDgIP78809pv9TUVAFAdOrUSbRo0UK89957YvHixcLHx0cEBgYKg8Eg7fvzzz8LuVwu2rZtKz766CMxb9484eHhISIjI62uFCwWi+jfv7+QyWTi4YcfFsuWLRMjR44UAMSzzz5r1XYAon379iIoKEgsWrRILFq0SHh4eIjg4GCxbNky0aZNG/Hhhx+K+fPnC0dHR9GvX78Kf04rVqwQAMT+/fulbf379xePPvqouHDhggAg/vrrL+mxjh07itatW0v/7tatm5g+fbr4+OOPxdKlS8XgwYMFALFs2TKr56noSkGj0YhBgwYJLy8vcfz48QrbXZopU6YIZ2dnYTKZyt2vsLBQODo6ikcffbTCc0ZHR4uAgADRuHFj8dRTT4n//e9/ok+fPgKAWLVqlbSfRqMRrVu3Fg4ODuK5554T//vf/0RUVJQAIJYsWWJ1Ttzyje3111+Xfq/Gjh0rPvvsM/Hwww8LAOLFF1+U9vv222+FQqEQUVFR4ttvvxXffvutOHLkiPTaHR0dxfPPPy++/PJL8d5774mRI0eK7777rtzXt3HjRtGhQwfx2muviZUrV4q5c+cKLy8vERISItRqtbRf8d9Tp06dRP/+/cXSpUvFnDlzhJ2dnZg4caLVOefPny8AiGHDholly5aJhx56SAQEBAgfH58KrxQMBoMIDAwUTZo0Edu3bxeXL18Wx44dE9HR0SI0NFTk5uaWe3xl3tvdu3cLR0dHERISIl5//XWxfPly8fTTT4uBAwdK++zZs0fY29uLli1bisWLF4sFCxYIHx8f4eXlJVJTU6X9bH0PhSh6/zt06CD8/f3FW2+9JZYsWSLCwsKEs7OzyMrKkva7fv26CAwMFEFBQeLNN98Uy5cvF6NGjRIAxMcffyztZzKZxIABAwQAcd9994lly5aJd999V/Tv319s3bpVqFQqsXz5cgFAjBkzRvrdKf67Lm77zaZNmyYAiPHjx4tPP/1UPPjggwKAuPfee632CwkJEREREcLPz0/MnTtXLFu2THTu3FnIZDJx+vRpab+5c+cKmUwmHnnkEfHFF1+IDz/8UEyePFksWrSowvepNNUKhb1794rMzExx+fJlERcXJ3x9fYVCoRCXL1+W9jWZTEKv11sdn5ubK/z8/MRDDz0kbSsOhUaNGomcnBxp+7Zt2wQAsWPHDmlbx44dhb+/v8jLy5O27d69WwCwCoWtW7cKAGLhwoVWzz9+/Hghk8nEhQsX/v1BAEKhUFj9MhZ/qDdp0kQUFBRI21955RUBwGrf0pw5c0YAEG+99ZYQQgij0ShcXFzEmjVrhBBC+Pn5iU8//VQIIURBQYGws7MTjzzyiHS8RqMpcc4hQ4aIsLAwq23lhUJhYaGIjo4WPj4+ViFcGUlJSUKpVIoHHnigwn3feustAUD88ssvFe4bHR0tAIgPP/xQ2qbX60XHjh1F48aNpS8CS5YsEQCsPoQNBoPo2bOncHV1tXpvygqFm3/XhBBizJgxolGjRlbbyrp95OHhIZ588skKX8+tSnv/jh49KgCIb775RtpW/Pc0cOBAYbFYpO3PPfecsLOzk37Pb9y4IRwdHcXw4cOt9ps7d64AYNOtr2PHjonmzZsLANJ/Xbp0EdeuXavwWFvfW5PJJEJDQ0VISEiJoLm53cXvc3Z2trTtr7/+EnK5XDz44IPStsq8hwCEo6Oj1d/2X3/9JQCIpUuXSttmzpwp/P39rYJCCCHuu+8+4eHhIb13X331VYlba7e+lvJuH90aCqdOnRIAxMMPP2y133//+18BQOzbt0/aFhISIgCI+Ph4aduNGzeEQqEQc+bMkbZ16NBBDB8+vMRzV1W1bh8NHDgQvr6+CAoKwvjx4+Hi4oLt27cjMDBQ2sfOzg6Ojo4AAIvFgpycHJhMJnTt2rXUy+9JkybBy8tL+nfxLamUlBQAwLVr13Dq1ClMmzYNHh4e0n6DBg1CmzZtrM61c+dO2NnZ4emnn7baPmfOHAghsGvXLqvtAwYMsBrS2qNHDwDAuHHj4ObmVmJ7cZvK0rp1azRq1AiHDh0CAPz1119Qq9XS6KJevXrh8OHDAICjR4/CbDZLHdSA9eV8fn4+srKyEB0djZSUFOTn55f73MXHDB48GAkJCThw4AA6duxY4TG30mg0mDBhApycnLBo0aJy942Pj8eCBQswceJE9O/f36bz29vb47HHHpP+7ejoiMceeww3btzAyZMnARS9j02aNMHkyZOl/RwcHPD0009DpVLh119/rfB5Hn/8cat/R0VFITs7GwUFBRUe6+npiWPHjuHq1as2vaZiN79/RqMR2dnZaNGiBTw9PUv93X/00UetbjVERUXBbDYjLS0NALB3714YDAbMnj3bar9nn33W5jZ5eXmhY8eOePnll7F161Z88MEHuHjxIiZMmACdTlfmcZV5b//880+kpqbi2Wefhaenp9Vjxe0u/juePn06vL29pcfbt2+PQYMGYefOnSXOa+t7OHDgQDRv3tzqnO7u7tLfqxACmzZtwsiRIyGEQFZWlvTfkCFDkJ+fL70/mzZtgo+PD2bPnl2iPVUZalr8up5//nmr7XPmzAEA/Pjjj1bb27RpI30GAoCvry8iIiKsPns8PT1x5swZJCUlVbo9palWKHz66afYs2cP4uLiMGzYMGRlZUGhUJTYb82aNWjfvr10v8vX1xc//vhjqR9swcHBVv8uDojc3FwAkP5AwsPDSxwbERFh9e+0tDQEBARYfaADRR/WN5+rrOcuDp2goKBStxe3qSwymQy9evWS+g4OHz6Mxo0bo0WLFgCsQ6H4f28OhcOHD2PgwIHS/VZfX1/MnTsXAGwKhWeffRa///479u7di8jIyAr3v5XZbMZ9992Hs2fPIi4uDgEBAWXum5CQgDFjxqBt27b48ssvbX6OgIAAuLi4WG1r2bIlAEj33tPS0hAeHl7iHndZ72NpKvq9Ks/ixYtx+vRpBAUFoXv37njjjTcq/EIAFI2Eeu2116T+LB8fH/j6+iIvL69Gf/d9fX2tvkiVJT8/H1FRUejZsyfeffddjB49GnPmzMGmTZtw6NAhrF69utTjKvveFvcptm3btsx9il/LrX+zQNH7mpWVBbVabbXd1vfw1v2K9y3eLzMzE3l5eVi5ciV8fX2t/iseXXfjxg3ptURERFR5RNut0tLSIJfLpc+AYk2aNIGnp2eFn0m3vhagaBRoXl4eWrZsiXbt2uGFF17A33//XeU2VisUunfvjoEDB2LcuHHYvn072rZtiylTpkClUkn7fPfdd5g+fTqaN2+OVatW4aeffsKePXvQv39/q07WYnZ2dqU+l6iDVUPLeu7qtKlPnz7Iz8/HP//8g8OHD1vNQejVqxfS0tJw5coVHDp0CAEBAQgLCwNQ9Ms4YMAAZGVl4aOPPsKPP/6IPXv24LnnngOAUn92txo9ejSEEFi0aJFN+9/qkUcewQ8//ICvv/663G+Hly9fxuDBg+Hh4YGdO3eWCOGGoDrv4cSJE5GSkoKlS5ciICAA77//PiIjI0tcad5q9uzZePvttzFx4kTExsZi9+7d2LNnDxo1alQvv/ubNm1CRkYGRo0aZbU9Ojoa7u7u0heTmzWk99bWn09F+xX/7O+//37s2bOn1P+KB3vUFluvMmx5zX379kVycjK++uorKbg7d+5cqS9nN6uZ+ENR4999913069cPy5Ytk8bmxsXFISwsDJs3b7b6Qbz++utVep7icdSlXSolJiaW2Hfv3r0oLCy0+mVOSEiwOldtunm+wuHDh60u9bt06QKFQoEDBw5IIzWK7dixA3q9Htu3b7f6tlDahLey3HvvvRg8eDCmT58ONze3So1ceuGFF7B69WosWbLE6rbNrbKzszF48GDo9Xr88ssv0kgWW129ehVqtdrqauH8+fMAIN3KCwkJwd9//w2LxWJ1tVDT72N5f6j+/v6YNWsWZs2ahRs3bqBz5854++23MXTo0DKPiYuLw7Rp0/Dhhx9K23Q6HfLy8qrUvpt/94u/PABF33xtueLJyMgAUHQFeDMhBMxmM0wmk9X2qr63xbduTp8+jYEDB5b7Wm79mwWK3lcfH58SV5A1xdfXF25ubjCbzWW2r1jz5s1x7NgxGI1GODg4lLpPZW4jhYSEwGKxICkpSbrSBYrem7y8vCr/Lnt7e2PGjBmYMWMGVCoV+vbtizfeeAMPP/xwpc9Vo0NSY2Ji0L17dyxZskS6P1mcdDcn27Fjx3D06NEqPYe/vz86duyINWvWWF2C79mzB2fPnrXad9iwYTCbzVi2bJnV9o8//hgymazcP+ia0rVrVyiVSqxduxZXrlyxulJQKBTo3LkzPv30U6jVaqtbR6X93PLz88u8xC/Lgw8+iP/973/4/PPP8dJLL9l0zPvvv48PPvgAc+fOxTPPPFPmfmq1GsOGDcOVK1ewc+fOUm/pVcRkMmHFihXSvw0GA1asWAFfX1906dIFQNH7eP36dWzYsMHquKVLl8LV1RXR0dGVft7SuLi4lPjANpvNJW71NG7cGAEBASWGEN7Kzs6uxLfYpUuXlvhQttXAgQPh4OCApUuXWp13yZIlNh1ffFtu/fr1Vtu3b98OtVqNTp06Sduq89527twZoaGhWLJkSYmfZ3G7b/47vnmf06dPY/fu3VZfkGqanZ0dxo0bh02bNuH06dMlHr95GPC4ceOQlZVV4jME+Pe1ODs7A4BNYV/8um59zz766CMAwPDhw216DTfLzs62+rerqytatGhR4e9nWWrsSqHYCy+8gAkTJuDrr7/G448/jhEjRmDz5s0YM2YMhg8fjtTUVHz++edo06aN1W2mynj33XcxfPhw9OnTBw899BBycnKwdOlSREZGWp1z5MiR6NevH+bNm4eLFy+iQ4cO2L17N7Zt24Znn33WqjOqtjg6OqJbt244ePAgFAqF9EFXrFevXtI3yZtDYfDgwXB0dMTIkSPx2GOPQaVS4YsvvkDjxo0rXWrgqaeeQkFBAebNmwcPDw+pX6I0W7ZswYsvvojw8HC0bt0a3333ndXjgwYNksa5T506FcePH8dDDz2Ec+fOWY0vd3V1tWlmcEBAAN577z1cvHgRLVu2xIYNG3Dq1CmsXLlS+mb26KOPYsWKFZg+fTpOnjyJZs2aIS4uDocPH8aSJUtq7JZGly5dsHfvXnz00UcICAhAaGgoIiIiEBgYiPHjx6NDhw5wdXXF3r178fvvv1tdAZRmxIgR+Pbbb+Hh4YE2bdrg6NGj2Lt3Lxo1alSl9vn6+uK///0v3n33XYwYMQLDhg3Dn3/+iV27dpWYW1KakSNHIjIyEm+++SbS0tLwn//8BxcuXMCyZcvg7++PmTNnSvtW572Vy+VYvnw5Ro4ciY4dO2LGjBnw9/dHQkICzpw5g59//hlA0ZePoUOHomfPnpg5cya0Wi2WLl0KDw+PWq9ftWjRIuzfvx89evTAI488gjZt2iAnJwd//PEH9u7di5ycHABFX6q++eYbPP/88zh+/DiioqKgVquxd+9ezJo1C6NHj4aTkxPatGmDDRs2oGXLlvD29kbbtm1L7VPp0KEDpk2bhpUrVyIvLw/R0dE4fvw41qxZg3vvvdeqwoGt2rRpg5iYGHTp0gXe3t44ceIE4uLi8NRTT1Xth1OVIUtlzVMQQgiz2SyaN28umjdvLkwmk7BYLOKdd94RISEhQqFQiE6dOokffvihxOzj4iGp77//folzopThXps2bRKtW7cWCoVCtGnTRmzevLnUGc2FhYXiueeeEwEBAcLBwUGEh4eL999/32poXPFz3DrssKw2VXbGcPEQ1l69epV4bPPmzQKAcHNzKzEHYPv27aJ9+/ZCqVSKZs2aiffee08aInfzcFhbZzS/+OKLpc5zuFnxELqy/rt5zkXxkLnS/rNlZnl0dLSIjIwUJ06cED179hRKpVKEhISU2r6MjAwxY8YM4ePjIxwdHUW7du3E6tWrS+x36+9K8eu5dZZ28e/wzT/HhIQE0bdvX+Hk5CQN8dTr9eKFF14QHTp0EG5ubsLFxUV06NBBfPbZZxW+vtzcXKnNrq6uYsiQISIhIaHE7OOy/p6K38ebf+Zms1ksWLBA+Pv7CycnJxETEyNOnz5t04xmIYTIyckRzz33nGjZsqVQKBTCx8dH3HfffSIlJcVqv+q+t0IIcejQITFo0CDp59a+fXurYaFCCLF3717Ru3dv4eTkJNzd3cXIkSPF2bNnrfapzHtY2t9x8eu59eeTkZEhnnzySREUFCQcHBxEkyZNxIABA8TKlSut9tNoNGLevHkiNDRU2m/8+PEiOTlZ2ufIkSOiS5cuwtHR0ep3sLR5CkajUSxYsEA6X1BQkHjllVeETqcr0ebShpre+ve+cOFC0b17d+Hp6SmcnJxEq1atxNtvv201t6syZELUQQ8uUSliYmKQlZVV6iU8EdWPO6p0NhERVQ9DgYiIJAwFIiKSsE+BiIgkvFIgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiIhIwlAgIiIJQ4GIiCQMBSIikjAUiKhOpaenY//+/UhPT6/vplApGApEVGdWrVqFkJAQ9O/fHyEhIVi1alV9N4luIRNCiPpuBBHd+dLT0xESEgKLxSJts7Ozw8WLFxEYGFiPLaOb8UqBiOpEUlKSVSAAgNlsxsGDB+upRVQahgIR1Ql/f/9St8+ZMweNGzdGdHQ0li9fDp1OV8cto5sxFIio1l29ehX9+vXD8OHDIZPJABTdOvryyy9x9epVXLp0CQ888ABiY2MRHByMoKAgTJkyBfHx8fXc8rsP+xSIqFadPn0a0dHRWLhwIZ544gk88cQTSEtLw8qVK8vsS0hLS8Onn36KH374ARkZGQgICMCoUaPw5JNPIiAgoI5fwd2FVwpEVGv27duHqKgorF69Gk888QQAwMHBAaGhoeV2LoeEhGDx4sU4e/YsMjMz8dFHHyEpKQndu3eHn58fBg0ahDVr1sBgMNTVS7lrMBSIqFasXbsW48ePx+7duzFq1Chpu06ng1KptPk8crkcgwYNQmxsLNLT05GcnIxRo0bhiy++QNOmTdGsWTPMmDEDv//+e228jLsOQ4GIatzixYvx3HPP4eTJk+jWrZvVY5UNhVu5urpi9uzZOHToEDIzM/Hzzz/D3d0dU6dOhY+PDzp37oyFCxciKyur1OM5ea58DAUiqlFPPvkkli1bhoSEBISGhpZ4XKfTwdnZucaeLyIiAp988gnOnz+PGzdu4PXXX8eJEyfQrl07NGnSBMOHD0dsbCwsFgsnz9mAHc1EVGNGjx6N5ORk/PHHH3B0dCx1n3vvvRd9+/bF888/X+vtycvLw8qVKxEXF4cLFy4gNzfX6nFOniuJVwpEVG0mkwndunVDYWEh/v777zIDASi6UnBycqqTdnl6euLFF1/E8ePHsWnTphKPm81mXLhwoU7acrtgKBBRtWg0GrRu3RoRERHYt28f5PLyP1YMBgNcXFzqqHX/unHjRoltdnZ2aNGiRZ23pSFjKBBRlV2/fh0tWrTA+PHj8d1339l0jMFgqNE+BVvExcXhiSeewNy5c6VtdnZ2WLFiBW8d3cK+vhtARLenM2fOoG/fvnjrrbcwa9Ysm4+r61BYuXIl5s6di6NHjyIiIgIHDhxAVFQUnnrqKQZCKRgKRGSz9PR0JCUlITMzE4899hi+/vprjB49ulLnMBgMcHV1raUWWnv33XexZMkSnDp1SgoAIQR69OjBQCgDQ4GIbLJq1So8+uijUqXTuXPnVjoQAMBoNNbJlcKcOXMQGxuLxMREeHp6Stt1Oh08PDxq/flvV+xTIKIKpaenWwUCUPQtvCoTwIxGI9zc3GqyeSVMnz4dO3bsQFJSklUgAIBer4e7u3utPv/tjKFARBUqbS0EIQQiIyPRu3dvLF26FBqNxqZzGY3GWh19NGLECJw6dQpnz54tdea0Xq+Hl5dXrT3/7Y6hQEQVCg8PLzHU1M7ODn/99ReeeOIJbN++HaGhoQgKCsLUqVNx6NChMs9lMplqpU/BYrGgV69eUKvV+OOPP2BvX/rdcYPBwFAoB0OBiCoUGBiIlStXWq2FsGLFCjRr1gz3338/9uzZg4yMDBw+fBhNmzbFI488Ah8fH3Tq1AlvvvmmVR2i2ggFk8mEDh06oEmTJti/f3+5cyWMRiNvH5WDZS6IyGbz58/Hr7/+iu+//77C0TsWiwU7d+7EqlWrcOzYMQgh0L17d8THxyM7O7vCSW620mg0aNeuHQYMGICVK1dWuH+TJk1w/fr1GnnuOxGvFIjIZkqlEv7+/jYN55TL5RgxYgS2bNmCq1ev4syZM+jevTtUKhUaN26MVq1a4YUXXkBaWlqV25OVlYXw8HBMmjTJpkCgijEUiMhmer2+3LpG5fH29sa8efPg4+ODrKwsrFq1CteuXUNUVBQaN26MAQMGYM2aNTCZTDad79KlS2jTpg3mzJmDd955p0ptopI4T4GIbFadUChWfMe6d+/e6N27N4CiW0Bff/01vvjiC7zwwgtwcXHBoEGD8NRTT6F9+/bSscWT54QQmDBhAj766CNMmzatWu0hawwFIrKZXq+HQqGo8fM6Oztj1qxZUrmMM2fO4NNPP8WYMWOQn5+P8PBwBAcHIzY2VjrmqaeeqlQgFBoKse3CNmhNWsSnxyPYLRhNXZvCwc6hxl/P7YyhQEQ2MxgM1b5SKB7BVJ7IyEh89tlnAIpGFn3++eeYPXu21T7Lly/HSy+9ZFP/Rp4uD7HnY3Em6wxMFhM2J22Gwk4BL4UXWnq1RKhnKILcgtBI2cim9t3JGApEZDO9Xl/t2ciVHfBob2+PyMjIEtuL10KoKBSytFnYkLABSXlJCPUIhZ3MDuGe4dCb9cjX5+PQ1UM4ePUgXB1c0cS5CVp5t0KQexCC3ILgZF836z40JAwFIrKZwWColdtHFSmePHfzrGpb1kK4rr6O9QnrcbHgIkLdQ6VbRTKZDEp7JZT2SvjBDxZhgcqgQroqHRfyLsBObgdPhSfCPMLQ3LM5gt2C0di5MezkdrX6OhsChgIR2awmQqEqt2eKJ8898sgjEELYtBbC5cLLiE2MRXphOsI8wmAvty9RqqOYXCaHu8Id7oqiSW0GswEFhgL8kfEHjl8/Dmd7Z/g4+aCVdyuEuIcgyD0I7o535gQ4hgIR2cxoNFYrFNLT06HX65Genl7p0tUzZ86El5cXXnjhBfz666/lHp+Sn4LYhFjc0NxAmGcY7GRF3/A1BRrYOVT8bd/RzhE+Tj7wcfKBEAJqkxo5uhzsSdsDuUwON0c3hLiHoIVnCwS7ByPANQAO8jujw5qhQEQ20+v1VV5f+ebS2yEhIVi5ciVmzpxZqXMolUp4e3uXGwiJOYmIOx+HHF0OwjzDIJf9Ox1Lp9bB3qFyH3symQyuDq5wdSgqzWGymFBgKMC57HP4K/MvKO2U8FR6IsIrAqEeoQh0C7ytO6wZCkRkM6PRWGrl0YrcWnrbYrHgsccew5AhQyp1xaBSqcp9/jNZZxCXFAeVQYUwj7ASH8y2XimUx15uD2+lN7yV3hBCQGfWIV+fj4PpBxGfHg83Rzc0cfn/Dmu3IAS6Bd5WHdYMBSKyWVVDobTS22azGXFxcXj22WdtPk9hYWGZt6/+yPgDWy9shd6sRzP3ZqV+U9cWauHgWHO3eWQyGZzsnYo+9F0Ai7Cg0FCISwWXcD73PBzkDkUd1p5haO7RHMHuRR3WN1+9NDQMBSKyWVVDITw8vMQ2mUyGzz77DAsXLkRgYCDGjh2LWbNmwcfHp8zzaDSaEs8vhMCx68ew48IOQAYEuwWXeetGp9HVaCjcSi6Tw0PhAQ9F0cpuBrMB+fp8nMg4gePXjsPZwRm+Tr6I8I4o6rB2C4KbY+0uOFRZDAUisllVQsFisWDKlCkIDw9HUlISgH9Lb8+cORMWiwU//fQTvvjiC7Rv3x5CCPznP//BI488gnvuuceqmuqtt4+EEDiYfhA7U3fC0c4RTVyalNsWbaEW9oq6+9hztHOEr7MvfOFb1GFtVCNTm4m0i2lFI54c3Ys6rL1aINgtGP6u/vXeYc1QICKbmUymSnU0m0wmdO7cGc2bN0d8fDz27duHadOm4ejRo1Jfglwux7BhwzBs2DAARZVPP/vsM8ydOxcPPvgg/Pz8MHr0aDz11FPQaDTS+s4WYcEvl37B7ou74ergCl9n3wrbo1Pr4Kio3ozsqpLJZHB1dIWro3WH9ZnsMziVeQpKOyW8ld5o6d0SzdybIdg9GF4KrzrvsGYoEJHNTCaTzVcKKpUKHTp0sFrnwM3NDR4eHuV2Lvv4+OC1117Da6+9BgA4cOAAPv/8c3Tv3h05OTnw8fHBhtgNcO7kjANXDsBT4YlGTo1sapNerYeDsmEMHb21w1pr0qLAUID49Hj8Kn6Fm6Mb/F390crr3w5rpX3lb91Vul21/gxEdMcwGo02XSncuHEDHTt2xIwZM/D2229L2wsLC+HgULkP5ZiYGMTExAAomquQmJiI+QvnIz0tHS4eLugQ3QHR90XDp2nZfRHF9Fo9HJX1c6VQHplMBmcHZzg7FF0FmYUZKoMKaflpSMxJlDqsW3i1QJhHGILcgmqtw5qhQEQ2M5vNFYZCamoqevTogXnz5uGZZ56xekyj0VQ6FG6mN+rhG+mLrvd1RRPnJsg+n4342Hh8PPNjmAwmNA1viu4juqPr0K6lrtFckFUAbaEWuRm58PJruOs028nsrDqs9WY9CvQFOH7tOH67+pvUYd3auzWC3YMR5BYk3ZaqLoYCEdnMZDLBxcWlzMdPnTqF/v3749NPP8XkyZNLPK7RaKo8I1pj1CDpRhLkAXL0dO0JVwdXuLdzR2i7UACAQWvAkW1HcHjTYWz+aDOULkq06dUG0fdFwz/MH0e2HMHRbUcBAKcPncaU+VPQa0yvKrWlrinsFEUd1s5FHdYqo6pEh3Uzj2YI9wxHkFsQ/F39YS+v2sc7Q4GIbGYymaSO3lvt27cP48aNw/r16zFkyJBS91Gr1VUqvV1oKMSm85two+AG2rVpJ80uvpmjkyNi7otBzH0xAID08+mI3xCPz2Z/Bp1KB02BRtpXWAS+X/g9Wvdq3aCvGEojk8ng5ugmDWU1Wowo0BfgdNZpnLrx/x3WTt5o6dUSzTyaIdgtGJ4KT5s7rBkKRGSzskIhLi4Ojz76KHbv3o1u3bqVebxWq610KNy8FoKjxRGuLrbdJglsGYgpr04BAJw9ehafzvrU6nGLxYLTh04jalxUpdrT0DjIHdDIqREaOTWy6rA+cPkAkA64O7gjwDUAEd4RCHYLRphnWLnnYygQkc0sFkuJ0UfLly/Hq6++iqNHjyIiIqLc47VabaXmOWRpsxCbGIvE3ESEeoTCbDBD4Vz5209yu9I7ZHd/tRs7lu6Ab5AvutzTBb3G9ILSufZH+NSWmzusm7g0gdliRqGxECn5KTiXcw6Ockcsjl5c7jkYCkRUKTdPJnvzzTexfPly/P333wgICKjwWLVabXOfwnX1dWxI3IDU/FSEuYfBwc4BZqMZCqfKhcKlc5fw5Qtfou/EvoiPjZdew+T5k9FrTC+YTCac+uUUftv6G3Z/tRtyOznCu4Sj76S+aN6xeaWeq6EpXhfCU+EJoKjDuiIMBSKqktmzZ2P79u1ITEyEu7ttawvodDqbQiG9MB0bEjdYrYUAAEaDEUpX27/JX/jzAj5/5nPMfG8mWvdsjeQ/k9Hlni7oPry71Jdgb2+PrkO6ouuQrgCAnGs5+DX2V6xdsBaqXBW8/L3QeWBnRE2IgrN76f0ptwuFXcU/e4YCEdnMbDZj//79+Pjjj5GUlISkpKRK9RHodLqKh7TmpyI2MRYZmgyEeYRZrXZmNtp+++jskbP46uWv8MT/npC+8Qsh0KJzi3I7l739vTHmmTEY88yYon6H+NM4vOUw9n+/HzKZDGEdwtBnXB9E9Iiwumq6UzAUiMgmq1atQlZWFvr37w8AWLlyZaU7jbVabblXCudzz2Nj4saitRA8wkpMzjIZTVC6VHylcGrfKXy34Ds8veJpBLcOtjre2c32b/tyuRztY9qjfUx7AEBBTgHi18dj80ebkZ+ZDw9fD3To1wF97+sLd+87YyU2hgIRVah4PYSbPfHEExg6dGil1kPQ6XRl3mo6k30Gm85vQqGhsNS1EADAbDJXGArHfzyOjYs3Ys7qOfAP87d6zGQwVer2063cvd0xYtYIjJg1AgBw7ug5HNp0CIvuW1S0eFBkCHqN6YV2fduVehWRm5GLzEuZ8A32bbBDYRkKRFShstZDePLJJ7F48eIKRx0V0+v1pQ5p/fPGn9iStAUGs6HMtRAAwGK2wMm17NtP8Rvj8cNnP+CldS+VWvbCbDLDxb3syXeV1bpna7Tu2RpA0QI+BzcexK4Vu7B2wVq4erqiXXQ7RE+Khre/N45sOYJ1C9dBWARkclmDnTzHUCCiCoWHh0Mul1sFg0wmg6urK4YMGQKNRoPIyEhMnz4dU6dOLbXEBFCyT0EIgePXj2N78nYIIRDkFlTuJCuL2VJm7aK93+zF3jV7MTd2Ljx9Pcs+3ql2ah85uztjyMwhGDKzaOLehT8v4ODGg/hwxocw6oxQ56ulfRvy5Lk7r5eEiGpcYGCgVOkUKFoP4YsvvsDatWtx8eJFXLp0CRMnTsTKlSvRpEkTNG/eHE8//TSSk5OtznPzGs9CCBy6cghbkrZALpMj0C2wymWif/z8R+z7bh/mx80vMxDqWotOLTDjnRl4+6e38eBbD5Z43GKxIPNyZj20rHwMBSKyycyZM+Hj44P9+/fj4sWLmDlzpvSYUqnEE088gcOHDyMrKwubNm2CSqVCv3790LhxYwwaNAhr166Vbh9ZhAX7Lu3DjuQdUNor4e/iX84zl2/Th5twdNtRvLr5Vbh61UxRuJrWKLBkaW+5XA7foIrXgKhrvH1ERDazs7OTyliXp2PHjvjqq68AFBXBW7VqFZYuXYqTJ0/i9OnT+OHID2g0qBFCQkJsXguhNOveWofE44l4bctrDbIkNgBoVVosfWwp2ka1xemDpwH8O3muod06AhgKRFTLnJ2dMXv2bMyePRs9evRAREwETpw6gSvPXIEwCQS2DETPMT3RaWCnSo37Xz13NS4nXMarm16FvWPD/CjTFGiwcPxC9B7bG8MfH473H3wfEd0jEDUhqkEGAsBQIKI6ojfrkanKRNPGTfHgmw/CXeEOnUaHw5sOY++avdjwzga4eLqgfUx7xEyOKfGhmZuRC7PRjNyMXMQuikXO9RzMj5vfYCeQFQdC1IQoDH1kKICiju7mnZo32EAAGApEVAlV7QjWGDXYcmELCjQFaOLZBO6KorkKSmclBjwwAAMeGAAASD6VjPjYeLz/wPuwmC0IbhOMXmN6QZWrwvp31kNYBObfMx8+gT54fdvrDTYQVHkqvD3+bcRMjpFGIwGAQWeAi0fNDYmtDQwFIqpVKoMKcefjcCrzFORCDk83zzL3bd6xuVSSQqvS4uDGg/jhsx9wLfma1X45V3OQn5nfIL9xq3JVWDhhIfpP7Y/BMwZbPWbUGxtsZ3ixhhmzRNTg3Dp5zRb5+nx8n/A9/rzxJ4LdgiHMwqYyFQDg5OqEwTMGY+JLE0tty6qXVuGf+H+q1K7aospVYeH4hRj44MASgQAUzahu6KHAKwUisonBYKjU7aNsbTY2JG5AYk4imnk0g8JOAbPJXOnJY77BvpDJZRAWYbW9ecfm2PHpDnz72rdw93FHx/4dy61BZNAZylxXoSYU5BTgnQnvYPBDg9F/av9S97GYLQ1+vQaGAhHZRKPRlDlT+VYZ6gxsSNyA5PxkhHkUrYUA/P+HolPlPhS9/LwwZf4UrH1zLQDrtRDGPDsGAJBwLAEHNx7Eu5PeBQQQ2iEUfcb1Qav/tJL6HTSFGtjZ25X5PNWRn5WPdya+g3seuQf9JverleeoKwwFIrKJSqWCnV3FH6pXVFewPmE9LhdeRnOP5lYLyJuNZqT8nQK/Zn6V6g/oNaYXmrVvhk8e+QQvf/9yiWNb9WiFVj1aFbUzV4VfN/yKLR9vkSqZdhrYCRHdI2DnUPOhkJeZh3cnvYuhjw6V1oe+nTEUiMgmWq22wiuFi/kXEZsYi2vqa2ju0dxqLYQjW45AlafC0seXVqkgnMJJAYWzosIwcfVyxfDHh2P448NhsVikSqb71+2HTq3Dly98iagJUYjoblsRv/LkZuRi0eRFGP7EcPSd0LfiA6o2eKtOMRSIyCYVhUJSbhI2nt+IbG02mns2t1oLITcjV7r9A1StIJxOo6v07R+5XI7I3pGI7B2JM4fPYNsn2+DXzA+x78WiMLsQXk280HlQ1VZVy7mWg0VTFmHUU6PQZ1yfSh3bkDEUiMgm5YXCmewz2Hx+MwoMBaWuhZCekF7iGIvFghXPrcCABwagy5AuFc450Kl0sHewh+u1HCjzVGXv5+kKlb93yfartFC6KDHyyZEY+eRI61XV1u2H3E6O5p2aI/q+aLTo1KLUcxevh+CgdMDyp5dj9OzR6D22d7nttiIq3qW+MRSIyCZarRYODg4ltp+6cQpbkrZAZ9aVuhaCKleF7xZ8V+o52/VthwPrDmDjexvLnc0MAHqNHsEyYMq9C2BvMJXZTpOjPdZtfb1EMOjVejgo/23/rauq5WXm4dfvf8W6t9ZBlaOCd4A3ugzpgj7j+sDJ1clqPQQA6DGqR+UC4TbBUCAim+h0OqsrBSEEfr/+O7Ylb4MQAsFuwSUCIedaDt6b+h6GzBwCpbOy1BFEwx8fDgBI+SsFv67/VZrNHNI2BH3G9UFkn0jI5XLo1Xr4QlZuIACAvcEEZZ6qRCjo1Do4KEqGWjFPX0+Mfno0Rj89GhaLBX/t/wtHtxzFL9/8AplMhoLsAqv9f//hd4ycNbJSt79qc0hsTWEoEJFNbr5SKF4LYWfKTtjL7eHvVrL0dUZaBj6Y9gHGPDsGve4t6lDe9r9tmPn+TPgGlVyOMqxDGMI6hAH4dxWzHZ/uwDevfQNPX0/4BPogpBodtTq1zuZKqnK5HJ0GdEKnAZ0AACd+OoHVr6y22qd4PQRbQ0GVq2qwhftu1vBbSEQNgk6ng4ODAyzCggOXDuCntJ/gbO+Mxs6NS+x7OeEyPnnsE0ydPxWdBnWStsvkMrTs2rLC57p1FbNzR89hy8dboEi9XuX2G3SGKq+61rxT8xIT6Cq7HoIqV1XulUpD0fCvZYioQSjuaP4p9SfsTN0JNwe3UgMh+VQyPnn0E8x4Z4ZVIACo8pDM1j1bo+vQrgjrFFa1E6CoT6Kqay4UT6ArVpX1ENT56ga75sPNeKVARDZRa9XINeZi76W9aKRsBC9lyQ/Es0fO4quXv8Ljnzxe+gieaoy+MWgNcHO07Zv2tZRrcGgdbLWtMKcQ9o72yM3IrVIhvV5jemH7p9tx/4L70bRF00qfg6FARHcMvVmPk+knkWfKQ2OnxlLp65v9uedPrFu4Ds+sfAZBrYJqvA0GncHme/Lb/7cdx97fBL9mfug2vBsA4NQvpwAAJ38+WemJcxIBtO3dtvLHoaifROGkqNKxdYmhQETl0pq02JK0BReyLsBF6VJqIBzddhSbP96M/675L/ya+ZV9smp0FBt0BjjYeKXw2JLHMDo8ACd2ncChuENI/TtVeqwqE+dqgrZQC0fnhn+lwD4FIiqTyqBCbGIsjl07BheZCxSKkt9096/dj62fbMUr379SbiBUt8S1UW+Ext0ZpgquFkyO9tB5usLe3h7/GfkfjHpqVKlt+fzZz3H8x+Mwmcof4lpTtIVam8uG1ydeKRBRqfL1+dh4fiP+yfoHIe4hOGs+C3sH64+MnSt24uDGg5gXN6/MktXFDDpDtVZKM+gMyA/0wbqtrxfNaBaAU04hLA720LsqAXnRZcitM5qP/XCs1PN1GtgJBzceRNwHcXB2d0a7vu3Qb0o/eJcyG7omaFVaOLk61cq5axJDgYhKyNZmIzYxFgk5CdJaCCaDyeqe/qYPN+GPPX/gta2v2fRhp9foq1W62mQwwVHpCJW/N1T+3rDTGzFx6nsAgC8OfQhjKd/Cv3/7eyT/mYxJcydhwzsbAFhPnLtn5j0AgLSzaTjw/QF8OONDmI1mBLcORq+xvdA+pn2NLfmpU+sa5Epxt2IoEJGVG5ob2JBQci0Ek8EkXSmsfXMtzv9+Hq9ved3msf86deUL2t3MqDNC4fzv7Ss7gwlmB3vYGU0wl3JL6bs3vkPyqWTMj5sPe0d77Pt2HwY8OABto9qW+HAOaROCaW9NK2qnRoeDGw/ipy9/wrq31sHN2w0d+3dE9OToKrcdKApFJ3deKRDRbeSq6irWJ6zHpcJLCPMIs1oLQVOgAQAsf3o5sq5k4dUtr9q86A5QNKRUbl/1b91Go9EqgAxuTlhx/BPAYgFuKa/x9byvcTnhMubFzZPaKIRAqx6tKvy2rnRWYtC0QRg0bRAA4PyJ84jfEI9F9y2CKleFlXNWou/EvtL6DbbSa/VwcXep1DH1gaFARACAtII0bEjYUOZaCH/s/kP69+T5kysVCMD/Dymt5DE3M+lNpS9lecvtnS9f/BIZqRl4ZcMrVs9nMpqqdE+/ZdeWaNm1JSwWC14e8DKahjdF3PtxyM/Kh3cTb3Qe0hl9J/St8NwGraHS5bnrA0OBiHAh9wJiz8ciS5tV6loI6xaus9r/+4XfI7JPZKXukes1+mqtfGYymaxuH5Vm5ZyVyL6ajVc2vFKiL8BsMkPpVvXRP3qNHg6ODlYL+Jw5dKZoAZ/vikpvh3cNR8zkGIS2Cy1xvE6tQ/bV7CpPnqsrDAWiu9zZ7LPYdH4T8g35CPMIswoEAMi8lGlV86fY2xPeRmi7UPQe19umDlmDzlDtjmal678f6q7XctB59W5ovd3w++PDsfzp5SjILsBLa18qvS0C1bpSUeerrTra5XI52vVth3Z92wEoCs8D3x/AN69+A3W+Gr6Bvug6rCt6j+mNE7tO4EbaDcQuisXGxRurPnmuDjAUiO5iN6+FEOoeWqL0NQD4BvuWWgxuzpo5OB1/Gru+2IW1b66Fp68nOg3shJjJMaXeJjFoDdW6UjCbzFa3j5yz8tF240EUBDTCtL9SodPo8MK3L9TYaKFbaQo05c6o9vLzwphnx2DMs2NgsVhw8ueTOLrtKHat2AV1vlrar74mz9mKoUB0FxJC4ETGCWy7sA0WYSl1LYRixcXgbl0LwT/UH/6h/lKH7NkjZ3Fw40G8OeZNyO3kaNmtJfpN7YeQNiEA/n9GcimL9NjKbDRbXSlofD1x/PFh2PfDcRgNRjy/+vlaCwSgKBRsrXIql8vRbWg3dBvaDed/P49PHv3E6vHKlt2uSwwForuMEAJHrh7BDyk/wE5mh0C3wAqP6TWmF45sPYLwLuHoO6lvqR9mbXq1QZtebQAULa6zf91+rH55NTQFGviF+sGzsWe1Rh+ZTWY4ufzbmVvQ2ANTjiXCzt8bz618psrntZW2QFul0ten9p0qsa2yZbfrEstcEN1FLMKC/Zf3Y9uFbXC0c0SAa4DNx8rkMgS2CrTp2623vzfGzRmHN7a/gXf2voNe9/bCxX8uIuFYAl4f+To2fbgJOddyKtd2swX2iqLvsRaLBR9N/wiOSkc8UweBAACaQk2lq5z++PmP+Gv/X5jw8gRpW1XKbtclXikQ3SXMFjP2pO3B3kt74e7oDh8nn0odbzFZbC5IdzN7e3v0HN0T2VeykXUlC/3v74/96/bjw+kfwmwyI6RtCPpO7IvWPVuXf/tHVvSBarFY8P4D78PH2xVPvDMDJoMRliq0q7K0Km2lFunZ9cUuHN58GK9ufhVOrk7Y9fmuMleda0gYCkR3AaPFiF0pu3Ag/UCZayFUxGwyV2vlMKPBCHtHewS3DpZmD2sKNIiPjce2T7bh63lflzvuX1gEEo4lYOPijfBr5ofF0e0xoO8LSOvVBj9++mSV22UrnUpnc+nrn1f/jIMbD2J+3Hzpddi66lx9YygQ3eEMZgO2J2/H4SuH0di5MTwUHlU6j9lkrvJylkBRmYpbrzSc3Z1xz8P34J6H75HG/R/ceBD7v9sPOwc7tOrRCv2m9sOlM5egylVh6eNLAQAD7h8AO4OxqF11tMSlVq2tcJ4EAOxZswf71+7H/Lj5t8VktVsxFIjuYMVrIRy/fhwBLgFwdXSt8rnM5upfKZQ36/fWcf9ZV7Kwf+1+fP7M58i9nmu17/cLv0frHxcg4fgnkJkrLsld3bLdQNHkNReP8stU7Fu7D7988wvmx82Hq2fVf9b1iR3NRHcotVGN2MRYHL92HIGugdUKBKCoo7c6oWAymCp1vE9TH0x4cQIefPPBkm2xWLB5yRbk5KlgtqHzV1OgqdYcCaAoFMpbD+HA+gPY/dXuokDwsv5ZG7QGyO1uj49bXikQ3YGktRAy/0GwezCc7KtfndNitlRpjWFhEhBCQJgEFEoFLMaib+0ymQwy+4qXYitt8hxQVAtp8dTFsFgsCOsQVm6ROm2htsRaEJVl0Bqs5kncLH5jPHat3IW5sXNLBAIAFOYW2ryUaH27PVpJRDbL0eVgQ8IGq7UQaoLFbIGDsnJXCsIkoL+mh8Vkgb+XP3xdfaFL1wEA5PZyKPwVFQZDWZPnJnu7ISCkMS60CsLXFzOkInWNAhqh29BuiBofJfWBaFXaal8pGLSGUm9/Hdp0CD8u/xGvbHgFHj6l99eo8lRVGrlVHxgKRHeQG5ob2JC4ARdyL1ithVATqnKlIISAxWSBzE4GnU4HmZ0Mcgc5hKVouxACMhsWbu41phe2/W+b1ZDOph9sQse1+yCmD8LwZ+6VitT9tf8vHI47jN2rd8NB4YDWPVsjtH1otT+UjXojnN2sO46PbjuK7cu245X1r8DT17PMY9V56koHan1hKBDdIa6qrmJD4gakFaShuWdzq7UQakJVbx8Vu3T+EjKvZMIszIjoGgF5Jbs0bx3SeaVbSwg7Ga52biFtk8vl6DSgEzoN6AQAyLiYgX1r92HrJ1uhLdTif4/9Dz3v7YkuQ7pUuiSGQW+As8e/ofDbjt+w5eMtePn7lyucd6DOV1frZ1eXGApEd4BLBZewIXEDrqlKroVQk6pUW0gAOz7bgYDmAYjoEYGEowk4tfcU3NzcYOdnh+gp0fBpWv5EutJGD6VFt0NadLtyj/Nr5ofJ8yYjtH0ojmw9gnYx7bB/7X5sfG8jXL1c0XFAR8RMialwfWmg6EqheJGc4zuPY9OHm/DyupdtWtNZU6CxeY5DfWMoEN3mkvOSEZtYtBZCmGfJ0tf1yWKx4OdVP8PZwxmDZhQVzgvvEg5hEshJz0H87nh8PPNjmI1mhESGIGpCFNr0blMifKq7vrNOrYPSSYl+k/uh3+R+AICkk0n4df2veHfSuwCA5p2aI2ZyDFp0alHqOQw6A66nXkfSH0nYsWwHXvzuRTQKaGTT82sLbZvj0BAwFIhuY+eyzyHufBzyDfkI9QhtUIFgMpnw0fSP0LlbZ/Qe3xvC9O/oIWER8PTzxMRXJuI+h/ugVWkRvzEe25duxzevfgMvfy90vacroiZEQemshLZQW6KYnp3eCCGXw2IvL7Ec5630an2Je/rhXcIR3iUcAJCflY8D3x/AujfXoTC3EI2DGqPb8G7odW8vOCodcWTLEeRl5OHLF74EAIx8cmSlCtppVdpyh7M2JAwFotvU35l/Y3PSZmhN2jLXQqgvJoMJb098Gy06tkCfCX2KOpXN1kNK5fZyqc1Ork4YMmMIhswYAovFgtPxp3Eo7hB++eYX2DvaIyQyBLhlnZ97/vsFQg6dwS8L7kfiqJ7ltkev1Zd7T9/DxwOjZ4/G6NmjYbFYcOKnEzi69Sh2rtgJB4UD8jLyrPb/cfmP6DGyh801jLQqLZzcqj8suC4wFIhuMzevhWAW5nLXQqgPBp0Bb094G63+0wqT502W5incqqx5CnK5HO1j2qN9THsAQOblTGxZsgV5mXl4qf9LaNqyKXqP7Y2R+v8vc2HDqCK9Rm9ziQ65XI7uw7qj+7DuAIDDWw9j3QLr5Ugrux6CXq0vc7hqQ8NQILqNCCFw9OpR7EjZUbQWgmtgwwoErQFvjXsL7WPaY8KLReWiZfYym4adlsU3yBdR46OQcy0Hz3/1PA5vOoy9a/Zia3oWGjfzQ8SFq+iVpyq1rITFYAEsgDAKuLq6wqL7/w5rOSB3tO1WW5uebUpdea4yt490Gl2J4awNFUOB6DZhERb8evlX/HTxJyjtlPBz8avvJlnRqrRYOH4hug7tijHPjKnRc+s0Otg72MNR6Yh+U/uh39SizuKEYwmIj43HvvELIZdbr/ZmMVigOa+BxWhBmG8YnNycoDqnAgDIHeRwbulsUzCUNXmuMuWv9Ro9nNx5+4iIaojZYsbeS3uxJ21PldZCqG2qPBXemfgOeo/tjeGPD6/x8+s1+lInn7Xq0UoqbZGbkYv9a/9/tbdCDZq1aobB9wyGb4gvdDodlB5KyB3lEGZRVGqjEjXySps8VxkGraHCYnoNBUOBqIEzWoz4KfUnHLh8AN5K7yqthVCbVLkqLBy/EP2m9MOQmUNq5TkMOkOJ2kEd1+yBvc6Ic2N6Qd3YE15+Xhj7/FiMfX4sTAYTTvxwAhkXM3Du+DmoVWrk5+TDv6U/FI6KEp3etqjOeggGHUOBiGqAwWzAjuQdOHT1EBo7VX0thNqSl5mHRfctwuCHBqP/1P619jx6Tckhpe3X7odrZj4uRreDurGn1WP2jvboPqw7VOdUSDmdgqRTSVC6KnEw9iDs7e3h19QPXkYvtOhe+pyEmmbUG+HiyVAgomrQmrTYemErjl07Bn8Xf7g5utV3k6zkXMvBoimLMGLWCPSd0LdWn8ugNZS4fZQ46j9Q5qmh9S7753I54TKuXLiCwTMGS1VS9So9riZcxZYlW5B5NRO+Qb7oPqI7eo7uWeaw1equx2AymODm1bDev7IwFIgaILVRjU3nN+FkxkkEugXCxaFhfcvMupKFxVMXY8xzY9BzdPlzBGqCQWcocaVw7KlR5R6z95u9sM+yR69RvWAns5Mmzzk4OCC4TTDmjJsDOAInfz6JI1uOYOeKnVA4K9A2qi3639/fqvSGTqWrVpVVs8lcZtnthoahQNTAFBgKsDFxI/7O/LvG1kKoKovOIg3FdHNzg1ljRlZ6Fj575jOMf3G8NJa/thl1RjgqrL/FB/6WgD6LN+LQixOQ/h/rdRR+/PxHnPjxBJ558xnIIS8amnoTuYO8aFiqXI5uQ7uh29BuAID08+k4sO5AidIbvsG+1a6yWqW6UfWAoUDUgOTqcrEhcQPOZZ+r0bUQqsKis6Dw70IIY1EoDBkwBFlHsnDp7CU89ORDCOofVGdtMegMcPG66WrJYkGPpdvgnXodPZZtQ3qPCKnUxc4VO3F482G8uvlVKBwVpY8yKmOeQmDLQNz/xv0AYFV6I/taNkwGE/as2SOV3qiUhjOVpEIMBaIGIlOTifWJ63Eh9wJCPULhaFe/pZaFRUAYBWR2MkAOaLQaJP6RiNDIULh5upVYCa02GfXWVwrNfv0bfmcvAQD8zlxC0NFzuNyrDXZ9sQuHNh3Cq5tfLXc9aFvcXHrj9KHT2Lh4I5J+T5JKb7Tu2Rr9p/aHf3P/Uo8368yAuej/e7h7wKz+/3/YAXbK2qliWxMYCkQNwDXVNaxPXI+0gjSEeYbBQd5wFmSR2cmg1+uh0+jgHeANJ4+6v51l0Bv+rTIqBLqt2CU9ZpHL0P2zHfjiTBoObjyI+XHzqx0It9KpdfD09cSsZbMAFJXe2PfdPnz29GcwaA1S6Y1OAztBLpfDrDOj8I9CqR9j2OBhKDhZAKBohrdbZ7cGGwwMBaJ6drngMtYnrsdV1dVaXQuhOhTOCrTs1hLZV7OR8FsClEolkn5OQs+JPdGsbbNaf36T3iTVLgo6eg6+ienSY3KLgN+ZS7BPz8L87Qvg7F7z5SS0BVqrkUm+Qb6Y9MokTMIkGHQGqfTGhnc2wN3HHf8Z8h90bN4Rdo52/6465yCDMBetVV18BdEQMRSI6lFKXgo2JG5ApjYTYZ5hsJM1vEAAiorXNQpohEYBjSBMAga1ATdO38DX876GpkAD/zB/9B7bG12Hdq2VDlWjwVh0H18IdP9sByxyGeQ33b4yAfi0iTe21FIlUq1KW+ZymqWV3jix/QSumK9Ar9dD6a6EwWCQiv/V5W23qmAoENWThJwEbEzciHxDPsI8GtbiODe7dfavMAvYO9hjyENDMOypYTAZTDi67SjiN8Qj7v04uHm7ofPgzoi+L7rUInVVYTIUXSkEHT0HvzOXSjxuD8A/MV3qW6hpOrXO5kVyWvVohfC24Sg4WQCTxYQLpy5Ue55DXWIoENWDfzL/waakTdAaG95aCMVkclnRLQ+jKBEMMgcZZPKiNts72iNqQhSiJkQBAM6fOI8D3x/AwnELYWdvh4geERjw4AA0bdG0ym0xGUxQOivQfdl2CJkMslJKcQtZUd/C5Z6tK1x0p7J0al3lRxwByLyUCZlcho79O9Zoe2oTQ4GoDgkhcDLjJLZd2AaTMCHYvWGthXAzuVIOt/aljzKSyWWQK0u/smnZtaVUIyj7ajZ++fYXrHh2BfQaPQIjAtFnfB906NehUreZTEYTXBQOcLuWW2ogAIBMCLhez4PcaIKlmnMKbqXX6OHeqOJ1nG+WfSUberUeEV0ipFCtSs2lusZQIKojQggcvXYUO5Ib5loIpSnrg99WjQIaYeJLEzHxpYkwaA04tPkQfl71M9a9tQ5ejb3QZWgXRE+KLvNbuFljhjALuDq7wl7hgnVfvojEbUfx5y8n8cCrk0qsZqb1dqvxQACKQkEZbPuVwq7Vu9AETRDWNgwwW4eBzF4GNMyuIwCATJS2JBIR1SiLsCA+PR67UndBYadAE5cm9d2kemWxWHD28FnEb4jHpXOXYO9oj8g+keh/f3/4hRStE2HWmJF/NB/CLJCfmQ83bzcYdAbo1Dp4NPaAVx8v2DnXzafrZ099hsioSERPiq5w350rduLwlsN4df2rcHAoJaDqeZ7CrI6zyn2cVwpEtcxsMeOXS79gT9oeuDq4wtfZ9hW77lRyuRxto9qibVRbAEDGxQz88u0vWPr4UpgMJoREhqDfuH7wNfsCMkBv0MNebQ91gRo+TX0gg6xOb8XodXqb1lje+81eHIw7iNe2vHbb1Dq6FUOBqBaZLCbsSt2FA5cPwEvpBW+ld303qUHya+aHKa9OAfD/5SU2xGPvN3sxNGYozMIMo96IAlMBmoQ1gRzyOr83b9QZK1xOM35jPPau2Vsrk+dqii03hhgKRLXEYDbgh+QfcPDqwQa5FkJD5eTqhCEzh2DAxAHIP5IPnU4HnUYHiKIrChc3F7i41W3VWKPeWO6kuN+2/4YfPvsB8zbOg6tXzQzDrWlCCFwqKDmc91YMBaJaoDPpsO3CNhy9drRBroVw25AVhURQq6LiewadAepcNQqyC7D6odXwDPZEvyn9EN4lvFabYTQYy1w57eTuk9j00Sa88v0r8PBpuMF/RXUFSvuKb2kxFIhqmMaowaakorUQmro2bXBrIdxubh4S6+DgUPTBK4DHlzyO+G3x+H7h91DlqtA4pDF6ju6JHqN6wN6+Zj/aTAZTqaHwT/w/WP/2evz3m//C27/h3hq8pr4GmUyGsS3HVrgvQ4GoBjWktRBudzI7GWR2/18v6NbJc3YyuHi5YMSsERgxawRMJhNO7DqBI1uOYNvSbXB2d0anAZ3Q7/5+cPeu3PyC0piMJijdrL9lJxxLwJpX1+C5Vc9JI6YaokxNJoxmI8aGj0UH3w4V7s8hqUQ1JFeXi9jEWJzNPotm7s2gsK+/tRDuFMXzFG4ls5OVOxw15a8U7F+3Hxf+uACZTIbwLuHo/0B/hLQJqVI7Xh7wMhb9skj6d/KpZCx/ejlmfz67yuesCzm6HOTr8zEybCSiAqNsmhfDKwWiGpClzcL6hPVIyk1qEGsh3CmqOg8hrEMYwjqEAQDyMvOw79t9WP3yamgKNQhoEYBeY3qh6z1VK96XdjYNy59ejsc/ebxBB0K+Ph+5ulwMaTbE5kAAeKVAVG3X1dexPmE9LhZcRKhHaINaC4GsmQwmHNl2BMe2H0Pm5Uy4NXJDl8FdEDM5psToIpPKJK2HsPTxpZj9+WxkpWfhuze+w6hnRiEiKqI+XoJNVAYVrqmvoX9wfwwLHVapcuwMBaJquFx4GesTitZCCPMIa5BrIVDZEo4l4Nf1vyL171TYOdih9X9aY8ADA+Dr54vcX3LLPdZrgBfsXRvezRaNUYN0VTr6NO2D0S1GV/pLCkOBqIpS8lOwIWEDMjWZCPUMbbBrIZBtsq5kYd+3+/BP/D/wdPPEg1MeLHd/z2hPOHg2rKtCnUmHtMI0dG/SHeNbjq/SGt8MBaIqSMxJxMbzG5Grz0Woe2iDXQuBqkaboYXqN1W5+zS0UDCYDUjNT0XHxh0xKWISnB2qtgJdw7v2IWrgTmedRtz5OGiMGoS5hzX4SqdUefaK2+uj0WQx4WL+RbRu1BrjWo6rciAADAWiSjmZcRJbL2yF0WxEiHsIA4HqndliRkp+Cpp7NcekiElwd6zevAyGApENhBD47dpv2JG8AzKZDEFuQQwEqncWYUFKfgqC3YIxKWISvJRe1T4nQ4GoAkIIxKfHY2fqTq6FcJeQ2Vcc+LbsU5uEEEjNT0UTlyaY1GoSGjs3rpHzsqOZqBxcC+HudfM8hVvJ7GX1OhxVCIG0gjR4KDwwpfUUhHqE1ti5eaVAVAaTxYSfU3/Gvsv7uBbCXaghzkEAigIhXZUOZwdnjG85vkYDAWAoEJXKaDbih5QfEJ8eD19nX3gqPOu7SUQAiiqe2snsMDZ8LCK8a35WNUOB6BbFayH8du03NHFpwrUQqMG4obkBs8WMsS3Hor1v+1p5DoYC0U2K10I4cf0EAt0CuRYCNRjZ2myojWqMaj4KXf261trzMBSI/l+hoRBx5+Nw6sYpBLkFVWsCEFFNytfnI0+fh3tC70Gfpn1qdTg0Q4EIQJ4uD7HnY3Em6wzXQqAGpdBQiBvaG+gf1B/9g/rX+vwYhgLd9bK0WdiQsAFJeVwLgRoWjVGDa+priGoahaGhQ+ukCi9Dge5qXAuBGiqdSYd0VTq6N+mOkc1Hwl5eNx/XDAW6a10uvIzYxFikq9IR5hFWZ390RBUxmA1IK0hDp8adcG+Le+v06pV/BXRXSslPQWxCLG5obiDMM4xrIVCDYbQYkZqfikifyGpXPK0KhgLddRJzEhF3Pg45uhyEeYZxLQRqMMwWM1LzUxHuFY4JLSfUyxwZhgLdVc5knUFcUhzUBjXCPLgWAjUcxRVPQ9xDaqziaVUwFOiu8UfGH9h6YSsMZgPXQqAGpbjiqb+LPyZFTKrXwosMBbrjCSFw7Pox7LiwA5CBayFQgyKEwMWCi/BWemNixEQEuAbUa3sYCnRHE0LgYPpB/Jj6I9dCoAZHCIHLhZfh6uiKCRET0MyjWX03iaFAdy6LsOCXS79g98XdXAuBGqRr6mtwkDtgbIuxaOnVsr6bAwDgsAu6Yx1MP4hdqbsAAI2cGtVza4isZWgyYBZmjG4xGu1829V3cyQMBbpj+bn4IdInEk72TkjOS0ZKXgpydDkwW8z13TS6y2Vrs6ExajAibAS6+HWp7+ZY4XKcdEcTQiBHl4PU/FQk5CTgQt4F5OvzIZPJ4OHoAU+lJ0tbUJ3K0+chW5uNYaHD0D+49gvcVRZDge4q+fp8pOan4nzueZzPPY9cXS6EEHBTuMFL4cVieFSrCg2FuK6+joHBAzE0bGiDnDjJUKC7lsaoQWpBKpJyk3A2+yxytDkwCzNcHVzhpfSC0l5Z302kO4jaqMYV1RX0bdoXo1qMarC1thgKRAD0Zj0u5l9Ecl4yzmSfQZY2C3qzHi4OLvBSeMHJ3qnBXebT7UNn0uFSwSX0COiBceHjGvQVKUOB6BZGixGXCy4jJT8Fp7NOI0OdAa1ZC6WdEl5KL7g6uDIgyGZ6sx4XCy6ic+POmBgxEU72TvXdpHIxFIjKYbaYcVV1VQqIq+qrUBvVcLRzhJfCC26Obg3yvjA1DEazESkFKWjbqC0mt5oMV0fX+m5ShRgKRDayCAsy1BlIzU/FmewzuFR4CYWGQjjIHeCp8IS7wp0luElitpiRnJ+McM9wTG09FZ5Kz/pukk0YCkRVIIRAti4bKXkpOJdzDqn5qcg35EMGGbwUXvBQeDTYjkSqfWZhRkpeUcXT+9vcDx8nn/puks0YCkQ1IE+Xh9SCVCTmJOJ87nnk6fIgIOCh8ICXwgsOdpwLcbcoLoHt7+KP+1vfD39X//puUqUwFIhqmMqgQmp+0VDXhNwEaRa1m2PRXAiFvaK+m0i1pLjiqZfSC1NbT0WIe0h9N6nSGApEtUhr0pYY6mqymODi4AJPhSeHut5BhBC4VHgJzg7OmNJqCsK9wuu7SVXCUCCqIwazAZcKLiE5Pxmns04jU5MJnUkHZwdneCo94WLvwoC4jV1RXYGdzA6TIiYh0ieyvptTZQwFonpgspiQXpiO1PxU/JP1D66rr0Nj1EBhr5CGujIgbh8Z6gwYLAaMbzm+wRW4qyyGAlE9swgLrqmvISUvBaezT+NK4RWojCo4yB3gpfSCu6M750I0YFnaLKgMKtzb4l70DOh524c5Q4GoARFC4IbmBlILUnE2+ywu5l9EgaEA9jJ7eCg84KnwhJ2ccyEailxdLnJ1uRgWNgz9gvrd9oEAMBSIGiwhBHL1uUjJS0FCTgKS85KRb8gHAJb9bgAKDAW4obmBgcEDcU/oPXfM1RxDgeg2UWAoQEpeCi7kXUBCTtFQVwgUDXVVsux3XVIb1biquoq+gX0xsvnIO2qiIkOB6DZUXPb7Qu4FnM0+i2xdNkwWE9wc3Fj2u5ZpTVpcLryMnv49MTZ87B03MZGhQHSbKy77XVy0r7jst7O9M7yUXnC2d74j7nU3BMU/6y5+XTAxYuIdGb4MBaI7iNFixOXCy0UjmbJOI0OTAa2JZb9rgtFsREp+Ctr6tMWU1lPg4uBS302qFQwFojtUcdnv1IJU/JP5D8t+V4PJYkJKfgpaerXE1NZT4aHwqO8m1RqGAtFdQAiBDE0GUvJScCb7DC4XXkahoRB2MjtpLgSHupbudq54WhUMBaK7THHZ79T8VJzLOYeUvBSp7LenwhOeCs87ajRNddxc8fSBNg+giUuT+m5SrWMoEN3l8vX5SMlP+bfstz4PAODu6H5Xl/0WQiC1IBWNlI0wtfVUBLsH13eT6gRDgYgkxWW/L+RdwLmcc3dt2W8hBC4XXoaLgwsmt5qMFl4t6rtJdYahQESl0pq0SCtIk+ZCZGmzYLAY4OrgeseX/ZYqnraahMhGt2/F06pgKBBRhYxmI9IK0qS5EDc0N6Az6+Bs7wxPhSdcHO6cst/X1ddhspgwIWICOjXuVN/NqXMMBSKqFJPFhCuqK9JciGvqa9CYNFDY3f5lv7O0WVAb1VLF07sRQ4GIqqy47HdqfirOZBUNdS0u++2p8ISHwuO2mQuRq8tFrj4Xw0OHIyYo5rYNtupiKBBRjRBCIFObiZT8FJzLPieV/ZbL5fB09GzQZb8L9AW4ob2BwSGDMbjZ4NsmyGoDQ4GIakWOLgcpeSlIzE3EhdwLyDPkQQZZUdlvhWeDGeqqMqpwVXUV0UHRGBU2qsEGV11hKBBRrSswFCA1PxVJuUlIyElAri4XFmGBu6M7PJWeUNjVz1BXrUmLSwWX0KtpL4xtcedVPK0KhgIR1anist/Juck4k30G2bpsmC3moqGuyqKhrnVBb9YjNT8V3Zp0w4SWE+7IiqdVwVAgonqjN+uRVpCG5LzkOi37XVzxtJ1vO0xuNfmOrXhaFQwFImoQist+p+al4p+sf2qt7LfJYkJyfjJaebXClNZT7uiKp1XBUCCiBsdsMeOq+ipS81NxOus0rqiuFJX9ljvCS1n1st/FFU+beTTD/a3vRyOnRrXQ+tsbQ4GIGjSp7Hd+ijQXorjst6fSEx6OHjaNGLIIC1LyUtDUrSmmtp56V1Q8rQqGAhHdNkot+63Ph1wmh4fCo8yy3zdXPL2/zf0Icguqh9bfHhgKRHTbytfnIzU/FYk5iUjMTbQq++2p8ISjnSOEELhUeAmuDq6Y2noqwjzD6rfRDRxDgYjuCGqjWpoLUVz222QxwU5mB1cHV9zX6j60btS6vpvZ4DEUiOiOozPpcLHgIpJzk3E+7zz6BfVDx8Yd67tZtwWGAhHd0YQQd21xu6q4e6s+EdFdgYFQOQwFIiKSMBSIiEjCUCAiIglDgYiIJAwFIiKSMBSIiEjCUCAiIglDgYiIJAwFIiKSMBSIiEjCUCAiIglDgYiIJAwFIiKSMBSIiEjCUCAiIglDgYiIJAwFIiKSMBSIiEjCUCAiIglDgYiIJAwFIiKSMBSIiEjCUCAiIglDgYiIJAwFIiKSMBSIiEjCUCAiIglDgYiIJAwFIiKSMBSIiEjCUCAiIglDgYiIJAwFIiKSMBSIiEjCUCAiIglDgYiIJAwFIiKSMBSIiEjCUCAiIsn/AU3TFhmZHOJPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# %matplotlib notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# plotting stuff:\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d', computed_zorder=False)\n",
    "\n",
    "plt.title(f\"Random Walk {N} points and {connect.shape[0]} connections\")\n",
    "\n",
    "# add a green carpet to understand better the 3D view\n",
    "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n",
    "x = [np.min(traj_xyz[:,0])-1, np.max(traj_xyz[:,0])+1,np.max(traj_xyz[:,0])+1,np.min(traj_xyz[:,0])-1]\n",
    "y = [-3,-3,5,5]\n",
    "z = [0,0,0,0]\n",
    "verts = [list(zip(x,y,z))]\n",
    "ax.add_collection3d(Poly3DCollection(verts, zorder=-20, alpha=0.5, color='tab:green'))\n",
    "\n",
    "# plot points\n",
    "ax.plot3D(traj_xyz[body_mask,0], traj_xyz[body_mask,1], traj_xyz[body_mask,2], '.k')\n",
    "# plot lines\n",
    "for line in connect:\n",
    "    ax.plot3D(line[:,0], line[:,1], line[:,2], '-k', lw=0.5)\n",
    "    \n",
    "# show CM\n",
    "ax.plot3D([xcm[0], xcm[0]], [xcm[1], xcm[1]], [0, xcm[2]], ':r')\n",
    "ax.plot3D(xcm[0], xcm[1], xcm[2], 'sr')\n",
    "ax.plot3D(xcm[0], xcm[1], 0, '^r')\n",
    "\n",
    "# highlight base\n",
    "ax.scatter3D(traj_xyz[base_mask, 0], traj_xyz[base_mask, 1], traj_xyz[base_mask, 2], marker='s', color='gray', zorder=10)\n",
    "\n",
    "# nice view\n",
    "ax.set_aspect('equal')\n",
    "ax.axis('off')\n",
    "\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}