data_exploration/figure/Supplemental/motionAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c3f8996b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#dataframes\n",
    "import pandas as pd\n",
    "import h5py\n",
    "\n",
    "# user defined functions\n",
    "\n",
    "#speedup\n",
    "import multiprocessing as mp\n",
    "\n",
    "#suppress warnings\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "pd.TimeSeries = pd.Series \n",
    "\n",
    "#math\n",
    "import numpy as np\n",
    "import math\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from scipy import signal\n",
    "import scipy \n",
    "import scipy.stats as st\n",
    "from scipy.stats import kurtosis\n",
    "from scipy.cluster import hierarchy\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "from scipy import fft\n",
    "\n",
    "#plots\n",
    "import pylab as plt\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import MaxNLocator\n",
    "import matplotlib\n",
    "from figurefirst import FigureLayout,mpl_functions\n",
    "import figurefirst\n",
    "import matplotlib.image as mpimg\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable # for colorbar\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "sns.set_style(\"whitegrid\")\n",
    "import figurefirst as fifi\n",
    "from IPython.display import display,SVG\n",
    "\n",
    "import networkx as nx\n",
    "import community.community_louvain as community_louvain # pip install community\n",
    "\n",
    "#misc\n",
    "import cairosvg\n",
    "import time\n",
    "np.set_printoptions(suppress=True)\n",
    "import fly_plot_lib.plot as fpl\n",
    "# matplotlib.rcParams['font.family'] = [ 'sans-serif']\n",
    "FS=7\n",
    "dir = '~/DataAnalysis/data/Sprints/HighRes/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "297ceb16",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = pd.read_hdf(dir+'Windy/WindyMASigned.h5')\n",
    "df2 = pd.read_hdf(dir+'NotWindy/NotWindyMASigned.h5')\n",
    "df3 = pd.read_hdf(dir+'Forest/ForestMASigned.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3a5cd6a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fft_of_data(data,dt):\n",
    "    data = data-np.mean(data)\n",
    "    N = data.shape[0] # Length of data.\n",
    "    try:\n",
    "        yf = 2.0/N * np.abs(fft.fft(data)[0:N//2])\n",
    "    except ZeroDivisionError:\n",
    "        yf= 0\n",
    "\n",
    "    xf = fft.fftfreq(N,dt)[:N//2]# Only take the first half of the data.\n",
    "\n",
    "    return(xf, yf) # xf = frequency array. yf = amplitude or power.\n",
    "\n",
    "def calculate_rotation_matrix(angle):\n",
    "    c, s = np.cos(angle), np.sin(angle)\n",
    "    return np.array([[c, -s], [s, c]])\n",
    "\n",
    "def rotate_vector(vector, rotation_matrix):\n",
    "    return np.dot(rotation_matrix, vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "028d6b2e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ4AAALPCAYAAACT7S4pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABKr0lEQVR4nO3df0xV9f8H8CeXoAgssiHemzrC5Y3J11a4NSu2BJGLcYWMuk0wjfiytPVd1Mz7hyJXnOuy1SYM22oF/ViLzIlTyIjZhlhpMRcZgkWQJldNGHPDH1wv7+8ffLjjxo973tfDuff0eT42NjjnXH0et/P03st9n1eYEEKAiEiCIdgBiEh/WBxEJI3FQUTSWBxEJI3FQUTSWBxEJM1vcTidTqSlpcFsNuPMmTOTHuPxeOBwOLBixQpkZGRg7969ivYRkT7d5u+A9PR0vPDCC8jPz5/ymIMHD+Ls2bNoamrC4OAgcnNzsWzZMsybN2/afUSkT36LY+nSpX7/kMbGRjz77LMwGAyYPXs2VqxYgcOHD6OoqGjafUqMjIxgaGgIERERCAsLU/QYIpqeEAJutxvR0dEwGOTfsfBbHEq4XC6YTCbvz0ajERcuXPC7T4mhoaEpXyIR0a1ZtGgRZs2aJf04VYpjJkVERAAYPcHIyMggp/Hv1KlTSE5ODnYMRfSSVS85Af1kHR4expkzZ7zXlyxVisNoNKKvrw9LliwB4PssY7p9Soy9PImMjMTtt9+uRtwZp5ecgH6y6iUnoK+sgb78V+XXsRaLBXv37sXIyAgGBgbQ3NyMzMxMv/uISJ/8PuPYuXMnmpqacPnyZbz44ouIjY1FQ0MD/vd//xf/93//h//5n/9BTk4Ofv75Z6xcuRIA8Morr2D+/PkAMO0+ItInv8WxdetWbN26dcL2999/3/t9eHg4HA7HpI+fbh8R6RM/OUpE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCRN0R3Aenp6YLfbMTg4iNjYWDidTiQkJPgc8+abb6Krq8v7c1dXF6qrq5Geno6qqip89tlnmDNnDgDgkUcewfbt29U7CyLSlKLi2L59O9auXYucnBwcOHAApaWl+Pjjj32Oqaio8H7f2dmJ9evXIzU11bstNzcXW7ZsUSk2EQWT35cq/f396OjoQHZ2NgAgOzsbHR0dGBgYmPIxX375JaxWqy5uLkxE8vw+43C5XIiPj0d4eDiA0Tt6zZkzBy6XC7Nnz55w/PDwMA4ePIja2lqf7Q0NDWhtbUVcXBxeffVVPPzww1JBT506JXV8MLW1tQU7gmJ6yaqXnIC+sgZK9fEIzc3NMJlMSEpK8m57/vnn8fLLLyMiIgLHjh3Dpk2b0NjYiHvuuUfxn5ucnKyLu0e3tbUhJSUl2DEU0UtWveQE9JP1xo0bt/Sfsd+XKkajERcvXoTH4wEwOgv20qVLMBqNkx6/b98+PPPMMz7b4uLivPMbHn/8cRiNRvz2228Bhyai4PJbHPfeey+SkpJw6NAhAMChQ4eQlJQ06cuUCxcuoK2tDVar1Wf7xYsXvd+fPn0a58+fx/3333+r2YkoSBS9VCkrK4PdbseePXtw1113wel0AoDPiAQA2L9/P5YvX467777b5/HvvPMOfv31VxgMBkRERKCiogJxcXEqnwoRaUVRcSxcuBB79+6dsH38iAQA2Lhx46SPHysaIvp34CdHiUgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpCkqjp6eHthsNmRmZsJms6G3t3fCMVVVVVi2bBlycnKQk5MDh8Ph3Xft2jW89tpryMjIgMViwbfffqvaCRCR9lQbjwBMPQLhgw8+QExMDL755hv09vYiPz8fTU1NiI6OvvUzICLNzch4hH/66quvYLPZAAAJCQlITk5GS0tLgJGJKNhUHY8w1QiEvr4+3Hfffd7jjEYjLly4IBWU4xFmhl6y6iUnoK+sgVJtPIIaIxCmw/EI6tNLVr3kBPSTNWTGI0w3AsFkMuH8+fPeY10uF+bOnRtwaCIKLtXGI0w3AsFisaCurg4A0Nvbi19++cVnriwR6Ytq4xGmG4Hw0ksvwW63IyMjAwaDATt27EBMTMzMnRURzSjVxiNMNwLhzjvvRGVlZQDxiCgU8ZOjRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0hQtq+/p6YHdbsfg4CBiY2PhdDqRkJDgc0x1dTUaGxu99+MoKSnx3qzHbrfju+++895G0GKxYOPGjeqeCRFpRrXxCEuWLEFhYSGioqLQ2dmJgoICtLa24o477gAAFBcXo6CgQP0zICLNqTYeITU1FVFRUQAAs9kMIQQGBwfVT0xEQafqeIQx9fX1WLBggc8NiWtqalBXV4f58+fjjTfewMKFC6WCcjzCzNBLVr3kBPSVNVCqjUcYc+LECezevRsffvihd1tJSQni4uJgMBhQX1+PoqIiNDc3e8tICY5HUJ9esuolJ6CfrCEzHgEATp48ic2bN6O6uhqJiYne7fHx8TAYRv+q3NxcXL16VXogExGFDtXGI7S3t6OkpASVlZVYvHixz77xoxOOHj0Kg8GA+Ph4NfITURCoNh7B4XDg+vXrKC0t9T6uoqICZrMZW7ZsQX9/P8LCwhATE4N3330Xt92m+qskItKIauMR9u3bN+Xja2tr5ZMRUcjiJ0eJSBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4ikKSqOnp4e2Gw2ZGZmwmazobe3d8IxHo8HDocDK1asQEZGhs8y/On2EZH+KCqOsfEIX3/9NdauXetzs54xBw8exNmzZ9HU1IS6ujpUVVXhr7/+8ruPiPTH7418xsYj1NTUABgdj1BeXo6BgQGf2wc2Njbi2WefhcFgwOzZs7FixQocPnwYRUVF0+7zRwgBABgeHg70HDV348aNYEdQTC9Z9ZIT0EfWsetp7PqSpdp4BJfLBZPJ5P3ZaDR6b0g83T5/3G43AODMmTOKjg8FehrloJeseskJ6Cur2+32Dk2TEfI3/oyOjsaiRYsQERGBsLCwYMch+lcQQsDtdiM6Ojqgx/stjvHjEcLDw6ccj2A0GtHX14clS5YA8H2WMd0+fwwGA2bNmiV1UkTkXyDPNMaoNh7BYrFg7969GBkZwcDAAJqbm5GZmel3HxHpT5hQ8O5Id3c37HY7rly54h2PkJiY6DMewePxYMeOHTh27BiA0dEJNpsNAKbdR0T6o6g4iIjG4ydHiUgai4OIpLE4iEgai4OIpLE4iEhayBTHra7A1ZKSrNXV1XjqqadgtVqxZs0aHD16VPugUJZ1zB9//IGHHnoITqdTu4D/oTRnY2MjrFYrsrOzYbVacfnyZW2DQlnW/v5+FBcXw2q1IisrC2VlZbh586amOZ1OJ9LS0mA2m6dcshHwNSVCxLp160R9fb0QQoj6+nqxbt26Ccfs379fFBYWCo/HI/r7+0Vqaqo4d+6c1lEVZW1paRFXr14VQghx+vRpkZKSIq5du6ZpTiGUZRVCiJs3b4qCggLx+uuvi7feekvLiEIIZTnb29tFVlaWuHTpkhBCiCtXrojr169rmlMIZVl37tzp/XccHh4WeXl5oqGhQdOcP/74o+jr6xPLly8XXV1dkx4T6DUVEs84xlbgZmdnAxhdgdvR0YGBgQGf46ZaZRuKWVNTUxEVFQUAMJvNEEJgcHAwJLMCwHvvvYcnn3wSCQkJmmYElOesra1FYWEh4uLiAACzZs3C7bffHpJZw8LCMDQ0hJGREQwPD8PtdiM+Pl7TrEuXLp2wNOSfAr2mQqI4pluB+8/jAl1lq3XW8err67FgwQLMnTtXq5gAlGft7OxEa2srNmzYoGm+MUpzdnd349y5c8jPz8fTTz+NPXv2BLwsfKazbtq0CT09PXjiiSe8XykpKZpmVSLQayokiuPf7MSJE9i9ezfefvvtYEeZlNvtxrZt2+BwOLwXQ6jyeDzo6upCTU0NPvnkE7S0tODAgQPBjjWpw4cPw2w2o7W1FS0tLfjpp580f3Y8k0KiOMavwAXgdwXuGJfLpfn/4kqzAsDJkyexefNmVFdXIzExUdOcgLKsf//9N86ePYvi4mKkpaXho48+whdffIFt27aFVE4AMJlMsFgsiIyMRExMDNLT09He3q5ZTpmsn376KVavXu1d3Z2Wlobjx49rmlWJQK+pkCgONVbghlrW9vZ2lJSUoLKyEosXL9Y04xglWU0mE44fP44jR47gyJEjWL9+PZ577jmUl5eHVE5g9P2E1tZW770kfvjhBzz44IOa5ZTJOm/ePLS0tAAYvdvW999/jwceeEDTrEoEfE2p+jbuLfj9999FXl6eWLlypcjLyxPd3d1CCCGKiopEe3u7EGL0nf/S0lKRnp4u0tPTxeeffx6yWdesWSMeffRRsXr1au9XZ2dnSGYdr7KyMii/VVGS0+PxiF27dgmLxSJWrVoldu3aJTweT0hm/fPPP8WGDRtEdna2yMrKEmVlZcLtdmuas7y8XKSmpoqkpCTx2GOPiVWrVk3IGeg1xdWxRCQtJF6qEJG+sDiISBqLg4ik+S2OW/28e6isLyEi9fi9y3l6ejpeeOEF5OfnT3nM+Eltg4ODyM3NxbJlyzBv3rxp9xGRPvktjqVLl/r9Q2ZqihsAjIyMYGhoiHNViFQkxs1VMRjk37FQZSDTTE1xA4ChoSFdTXEj0pNFixYFNLco5Ce5RUREABg9wcjIyCCn8e/UqVNITk4OdgxF9JJVLzkB/WQdHh7GmTNnvNeXLFWKY6amuAHwvjyJjIzUfAl1oPSSE9BPVr3kBPSVNdCX/6r8OpZT3Ij+u/h9xrFz5040NTXh8uXLePHFFxEbG4uGhgafKW45OTn4+eefsXLlSgDAK6+8gvnz5wPAtPuISJ/8FsfWrVuxdevWCdvff/997/fh4eFwOByTPn66fUSkT/zkKBFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTRFdwDr6emB3W7H4OAgYmNj4XQ6kZCQ4HPMm2++ia6uLu/PXV1dqK6uRnp6OqqqqvDZZ59hzpw5AIBHHnkE27dvV+8siEhTiopj+/btWLt2LXJycnDgwAGUlpbi448/9jmmoqLC+31nZyfWr1+P1NRU77bc3Fxs2bJFpdhEFEx+X6r09/ejo6MD2dnZAIDs7Gx0dHRgYGBgysd8+eWXsFqturi5MBHJ8/uMw+VyIT4+HuHh4QBG7+g1Z84cuFwuzJ49e8Lxw8PDOHjwIGpra322NzQ0oLW1FXFxcXj11Vfx8MMPSwU9deqU1PHB1NbWFuwIiuklq15yAvrKGijVxyM0NzfDZDIhKSnJu+3555/Hyy+/jIiICBw7dgybNm1CY2Mj7rnnHsV/bnJysi7uHt3W1oaUlJRgx1BEL1n1khPQT9YbN27c0n/Gfl+qGI1GXLx4ER6PB8DoLNhLly7BaDROevy+ffvwzDPP+GyLi4vzzm94/PHHYTQa8dtvvwUcmoiCy29x3HvvvUhKSsKhQ4cAAIcOHUJSUtKkL1MuXLiAtrY2WK1Wn+0XL170fn/69GmcP38e999//61mJ6IgUfRSpaysDHa7HXv27MFdd90Fp9MJAD4jEgBg//79WL58Oe6++26fx7/zzjv49ddfYTAYEBERgYqKCsTFxal8KkSkFUXFsXDhQuzdu3fC9vEjEgBg48aNkz5+rGiI6N+BnxwlImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSpqg4enp6YLPZkJmZCZvNht7e3gnHVFVVYdmyZcjJyUFOTg4cDod337Vr1/Daa68hIyMDFosF3377rWonQETaU208AjD1CIQPPvgAMTEx+Oabb9Db24v8/Hw0NTUhOjr61s+AiDQ3I+MR/umrr76CzWYDACQkJCA5ORktLS0BRiaiYFN1PMJUIxD6+vpw3333eY8zGo24cOGCVFCOR5gZesmql5yAvrIGSrXxCGqMQJgOxyOoTy9Z9ZIT0E/WkBmPMN0IBJPJhPPnz3uPdblcmDt3bsChiSi4VBuPMN0IBIvFgrq6OgBAb28vfvnlF5+5skSkL6qNR5huBMJLL70Eu92OjIwMGAwG7NixAzExMTN3VkQ0o1QbjzDdCIQ777wTlZWVAcQjolDET44SkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJU7SsvqenB3a7HYODg4iNjYXT6URCQoLPMdXV1WhsbPTej6OkpMR7sx673Y7vvvvOextBi8WCjRs3qnsmRKQZ1cYjLFmyBIWFhYiKikJnZycKCgrQ2tqKO+64AwBQXFyMgoIC9c+AiDSn2niE1NRUREVFAQDMZjOEEBgcHFQ/MREFnarjEcbU19djwYIFPjckrqmpQV1dHebPn4833ngDCxculArK8QgzQy9Z9ZIT0FfWQKk2HmHMiRMnsHv3bnz44YfebSUlJYiLi4PBYEB9fT2KiorQ3NzsLSMlOB5BfXrJqpecgH6yhsx4BAA4efIkNm/ejOrqaiQmJnq3x8fHw2AY/atyc3Nx9epV6YFMRBQ6VBuP0N7ejpKSElRWVmLx4sU++8aPTjh69CgMBgPi4+PVyE9EQaDaeASHw4Hr16+jtLTU+7iKigqYzWZs2bIF/f39CAsLQ0xMDN59913cdpvqr5KISCOqjUfYt2/flI+vra2VT0ZEIYufHCUiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpKmqDh6enpgs9mQmZkJm82G3t7eCcd4PB44HA6sWLECGRkZPsvwp9tHRPqjqDjGxiN8/fXXWLt2rc/NesYcPHgQZ8+eRVNTE+rq6lBVVYW//vrL7z4i0h+/N/IZG49QU1MDYHQ8Qnl5OQYGBnxuH9jY2Ihnn30WBoMBs2fPxooVK3D48GEUFRVNu88fIQQAYHh4ONBz1NyNGzeCHUExvWTVS05AH1nHrqex60uWauMRXC4XTCaT92ej0ei9IfF0+/xxu90AgDNnzig6PhToaZSDXrLqJSegr6xut9s7NE1GyN/4Mzo6GosWLUJERATCwsKCHYfoX0EIAbfbjejo6IAe77c4xo9HCA8Pn3I8gtFoRF9fH5YsWQLA91nGdPv8MRgMmDVrltRJEZF/gTzTGKPaeASLxYK9e/diZGQEAwMDaG5uRmZmpt99RKQ/YULBuyPd3d2w2+24cuWKdzxCYmKiz3gEj8eDHTt24NixYwBGRyfYbDYAmHYfEemPouIgIhqPnxwlImksDiKSxuIgImksDiKSxuIgImkhUxy3ugJXS0qyVldX46mnnoLVasWaNWtw9OhR7YNCWdYxf/zxBx566CE4nU7tAv6H0pyNjY2wWq3Izs6G1WrF5cuXtQ0KZVn7+/tRXFwMq9WKrKwslJWV4ebNm5rmdDqdSEtLg9lsnnLJRsDXlAgR69atE/X19UIIIerr68W6desmHLN//35RWFgoPB6P6O/vF6mpqeLcuXNaR1WUtaWlRVy9elUIIcTp06dFSkqKuHbtmqY5hVCWVQghbt68KQoKCsTrr78u3nrrLS0jCiGU5WxvbxdZWVni0qVLQgghrly5Iq5fv65pTiGUZd25c6f333F4eFjk5eWJhoYGTXP++OOPoq+vTyxfvlx0dXVNekyg11RIPOMYW4GbnZ0NYHQFbkdHBwYGBnyOm2qVbShmTU1NRVRUFADAbDZDCIHBwcGQzAoA7733Hp588kkkJCRomhFQnrO2thaFhYWIi4sDAMyaNQu33357SGYNCwvD0NAQRkZGMDw8DLfbjfj4eE2zLl26dMLSkH8K9JoKieKYbgXuP48LdJWt1lnHq6+vx4IFCzB37lytYgJQnrWzsxOtra3YsGGDpvnGKM3Z3d2Nc+fOIT8/H08//TT27NkT8LLwmc66adMm9PT04IknnvB+paSkaJpViUCvqZAojn+zEydOYPfu3Xj77beDHWVSbrcb27Ztg8Ph8F4Mocrj8aCrqws1NTX45JNP0NLSggMHDgQ71qQOHz4Ms9mM1tZWtLS04KefftL82fFMConiGL8CF4DfFbhjXC6X5v+LK80KACdPnsTmzZtRXV2NxMRETXMCyrL+/fffOHv2LIqLi5GWloaPPvoIX3zxBbZt2xZSOQHAZDLBYrEgMjISMTExSE9PR3t7u2Y5ZbJ++umnWL16tXd1d1paGo4fP65pViUCvaZCojjUWIEbalnb29tRUlKCyspKLF68WNOMY5RkNZlMOH78OI4cOYIjR45g/fr1eO6551BeXh5SOYHR9xNaW1u995L44Ycf8OCDD2qWUybrvHnz0NLSAmD0blvff/89HnjgAU2zKhHwNaXq27i34Pfffxd5eXli5cqVIi8vT3R3dwshhCgqKhLt7e1CiNF3/ktLS0V6erpIT08Xn3/+echmXbNmjXj00UfF6tWrvV+dnZ0hmXW8ysrKoPxWRUlOj8cjdu3aJSwWi1i1apXYtWuX8Hg8IZn1zz//FBs2bBDZ2dkiKytLlJWVCbfbrWnO8vJykZqaKpKSksRjjz0mVq1aNSFnoNcUV8cSkbSQeKlCRPrC4iAiaSwOIpLmtzhu9fPuobK+hIjU4/cu5+np6XjhhReQn58/5THjJ7UNDg4iNzcXy5Ytw7x586bdR0T65Lc4li5d6vcPmakpbgAwMjKCoaEhzlUhUpEYN1fFYJB/x0KVgUwzNcUNAIaGhnQ1xY1ITxYtWhTQ3KKQn+QWEREBYPQEIyMjg5zGv1OnTiE5OTnYMRTRS1a95AT0k3V4eBhnzpzxXl+yVCmOmZriBsD78iQyMlLzJdSB0ktOQD9Z9ZIT0FfWQF/+q/LrWE5xI/rv4vcZx86dO9HU1ITLly/jxRdfRGxsLBoaGnymuOXk5ODnn3/GypUrAQCvvPIK5s+fDwDT7iMiffJbHFu3bsXWrVsnbH///fe934eHh8PhcEz6+On2EZE+8ZOjRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0hTdAaynpwd2ux2Dg4OIjY2F0+lEQkKCzzFvvvkmurq6vD93dXWhuroa6enpqKqqwmeffYY5c+YAAB555BFs375dvbMgIk0pKo7t27dj7dq1yMnJwYEDB1BaWoqPP/7Y55iKigrv952dnVi/fj1SU1O923Jzc7FlyxaVYhNRMPl9qdLf34+Ojg5kZ2cDALKzs9HR0YGBgYEpH/Pll1/CarXq4ubCRCTP7zMOl8uF+Ph4hIeHAxi9o9ecOXPgcrkwe/bsCccPDw/j4MGDqK2t9dne0NCA1tZWxMXF4dVXX8XDDz8sFfTUqVNSxwdTW1tbsCMoppeseskJ6CtroFQfj9Dc3AyTyYSkpCTvtueffx4vv/wyIiIicOzYMWzatAmNjY245557FP+5ycnJurh7dFtbG1JSUoIdQxG9ZNVLTkA/WW/cuHFL/xn7faliNBpx8eJFeDweAKOzYC9dugSj0Tjp8fv27cMzzzzjsy0uLs47v+Hxxx+H0WjEb7/9FnBoIgouv8Vx7733IikpCYcOHQIAHDp0CElJSZO+TLlw4QLa2tpgtVp9tl+8eNH7/enTp3H+/Hncf//9t5qdiIJE0UuVsrIy2O127NmzB3fddRecTicA+IxIAID9+/dj+fLluPvuu30e/8477+DXX3+FwWBAREQEKioqEBcXp/KpEJFWFBXHwoULsXfv3gnbx49IAICNGzdO+vixoiGifwd+cpSIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEiaouLo6emBzWZDZmYmbDYbent7JxxTVVWFZcuWIScnBzk5OXA4HN59165dw2uvvYaMjAxYLBZ8++23qp0AEWlPtfEIwNQjED744APExMTgm2++QW9vL/Lz89HU1ITo6OhbPwMi0tyMjEf4p6+++go2mw0AkJCQgOTkZLS0tAQYmYiCTdXxCFONQOjr68N9993nPc5oNOLChQtSQTkeYWboJatecgL6yhoo1cYjqDECYTocj6A+vWTVS05AP1lDZjzCdCMQTCYTzp8/7z3W5XJh7ty5AYcmouBSbTzCdCMQLBYL6urqAAC9vb345ZdffObKEpG+qDYeYboRCC+99BLsdjsyMjJgMBiwY8cOxMTEzNxZEdGMUm08wnQjEO68805UVlYGEI+IQhE/OUpE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCRN0bL6np4e2O12DA4OIjY2Fk6nEwkJCT7HVFdXo7Gx0Xs/jpKSEu/Neux2O7777jvvbQQtFgs2btyo7pkQkWZUG4+wZMkSFBYWIioqCp2dnSgoKEBrayvuuOMOAEBxcTEKCgrUPwMi0pxq4xFSU1MRFRUFADCbzRBCYHBwUP3ERBR0qo5HGFNfX48FCxb43JC4pqYGdXV1mD9/Pt544w0sXLhQKijHI8wMvWTVS05AX1kDpdp4hDEnTpzA7t278eGHH3q3lZSUIC4uDgaDAfX19SgqKkJzc7O3jJTgeAT16SWrXnIC+skaMuMRAODkyZPYvHkzqqurkZiY6N0eHx8Pg2H0r8rNzcXVq1elBzIRUehQbTxCe3s7SkpKUFlZicWLF/vsGz864ejRozAYDIiPj1cjPxEFgWrjERwOB65fv47S0lLv4yoqKmA2m7Flyxb09/cjLCwMMTExePfdd3Hbbaq/SiIijag2HmHfvn1TPr62tlY+GRGFLH5ylIiksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISJqi4ujp6YHNZkNmZiZsNht6e3snHOPxeOBwOLBixQpkZGT4LMOfbh8R6Y+i4hgbj/D1119j7dq1PjfrGXPw4EGcPXsWTU1NqKurQ1VVFf766y+/+4hIf/zeyGdsPEJNTQ2A0fEI5eXlGBgY8Ll9YGNjI5599lkYDAbMnj0bK1aswOHDh1FUVDTtPn+EEACA4eHhQM9Rczdu3Ah2BMX0klUvOQF9ZB27nsauL1mqjUdwuVwwmUzen41Go/eGxNPt88ftdgMAzpw5o+j4UKCnUQ56yaqXnIC+srrdbu/QNBkhf+PP6OhoLFq0CBEREQgLCwt2HKJ/BSEE3G43oqOjA3q83+IYPx4hPDx8yvEIRqMRfX19WLJkCQDfZxnT7fPHYDBg1qxZUidFRP4F8kxjjGrjESwWC/bu3YuRkREMDAygubkZmZmZfvcRkf6ECQXvjnR3d8Nut+PKlSve8QiJiYk+4xE8Hg927NiBY8eOARgdnWCz2QBg2n1EpD+KioOIaDx+cpSIpLE4iEgai4OIpLE4iEgai4OIpIVMcdzqClwtKclaXV2Np556ClarFWvWrMHRo0e1DwplWcf88ccfeOihh+B0OrUL+B9KczY2NsJqtSI7OxtWqxWXL1/WNiiUZe3v70dxcTGsViuysrJQVlaGmzdvaprT6XQiLS0NZrN5yiUbAV9TIkSsW7dO1NfXCyGEqK+vF+vWrZtwzP79+0VhYaHweDyiv79fpKaminPnzmkdVVHWlpYWcfXqVSGEEKdPnxYpKSni2rVrmuYUQllWIYS4efOmKCgoEK+//rp46623tIwohFCWs729XWRlZYlLly4JIYS4cuWKuH79uqY5hVCWdefOnd5/x+HhYZGXlycaGho0zfnjjz+Kvr4+sXz5ctHV1TXpMYFeUyHxjGNsBW52djaA0RW4HR0dGBgY8DluqlW2oZg1NTUVUVFRAACz2QwhBAYHB0MyKwC89957ePLJJ5GQkKBpRkB5ztraWhQWFiIuLg4AMGvWLNx+++0hmTUsLAxDQ0MYGRnB8PAw3G434uPjNc26dOnSCUtD/inQayokimO6Fbj/PC7QVbZaZx2vvr4eCxYswNy5c7WKCUB51s7OTrS2tmLDhg2a5hujNGd3dzfOnTuH/Px8PP3009izZ0/Ay8JnOuumTZvQ09ODJ554wvuVkpKiaVYlAr2mQqI4/s1OnDiB3bt34+233w52lEm53W5s27YNDofDezGEKo/Hg66uLtTU1OCTTz5BS0sLDhw4EOxYkzp8+DDMZjNaW1vR0tKCn376SfNnxzMpJIpj/ApcAH5X4I5xuVya/y+uNCsAnDx5Eps3b0Z1dTUSExM1zQkoy/r333/j7NmzKC4uRlpaGj766CN88cUX2LZtW0jlBACTyQSLxYLIyEjExMQgPT0d7e3tmuWUyfrpp59i9erV3tXdaWlpOH78uKZZlQj0mgqJ4lBjBW6oZW1vb0dJSQkqKyuxePFiTTOOUZLVZDLh+PHjOHLkCI4cOYL169fjueeeQ3l5eUjlBEbfT2htbfXeS+KHH37Agw8+qFlOmazz5s1DS0sLgNG7bX3//fd44IEHNM2qRMDXlKpv496C33//XeTl5YmVK1eKvLw80d3dLYQQoqioSLS3twshRt/5Ly0tFenp6SI9PV18/vnnIZt1zZo14tFHHxWrV6/2fnV2doZk1vEqKyuD8lsVJTk9Ho/YtWuXsFgsYtWqVWLXrl3C4/GEZNY///xTbNiwQWRnZ4usrCxRVlYm3G63pjnLy8tFamqqSEpKEo899phYtWrVhJyBXlNcHUtE0kLipQoR6QuLg4iksTiISJrf4rjVz7uHyvoSIlKP37ucp6en44UXXkB+fv6Ux4yf1DY4OIjc3FwsW7YM8+bNm3YfEemT3+JYunSp3z9kpqa4AcDIyAiGhoY4V4VIRWLcXBWDQf4dC1UGMs3UFDcAGBoa0tUUNyI9WbRoUUBzi0J+kltERASA0ROMjIwMchr/Tp06heTk5GDHUEQvWfWSE9BP1uHhYZw5c8Z7fclSpThmaoobAO/Lk8jISM2XUAdKLzkB/WTVS05AX1kDffmvyq9jOcWN6L+L32ccO3fuRFNTEy5fvowXX3wRsbGxaGho8JnilpOTg59//hkrV64EALzyyiuYP38+AEy7j4j0yW9xbN26FVu3bp2w/f333/d+Hx4eDofDMenjp9tHRPrET44SkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJU3QHsJ6eHtjtdgwODiI2NhZOpxMJCQk+x7z55pvo6ury/tzV1YXq6mqkp6ejqqoKn332GebMmQMAeOSRR7B9+3b1zoKINKWoOLZv3461a9ciJycHBw4cQGlpKT7++GOfYyoqKrzfd3Z2Yv369UhNTfVuy83NxZYtW1SKTUTB5PelSn9/Pzo6OpCdnQ0AyM7ORkdHBwYGBqZ8zJdffgmr1aqLmwsTkTy/zzhcLhfi4+MRHh4OYPSOXnPmzIHL5cLs2bMnHD88PIyDBw+itrbWZ3tDQwNaW1sRFxeHV199FQ8//LBU0FOnTkkdH0xtbW3BjqCYXrLqJSegr6yBUn08QnNzM0wmE5KSkrzbnn/+ebz88suIiIjAsWPHsGnTJjQ2NuKee+5R/OcmJyfr4u7RbW1tSElJCXYMRfSSVS85Af1kvXHjxi39Z+z3pYrRaMTFixfh8XgAjM6CvXTpEoxG46TH79u3D88884zPtri4OO/8hscffxxGoxG//fZbwKGJKLj8Fse9996LpKQkHDp0CABw6NAhJCUlTfoy5cKFC2hra4PVavXZfvHiRe/3p0+fxvnz53H//fffanYiChJFL1XKyspgt9uxZ88e3HXXXXA6nQDgMyIBAPbv34/ly5fj7rvv9nn8O++8g19//RUGgwERERGoqKhAXFycyqdCRFpRVBwLFy7E3r17J2wfPyIBADZu3Djp48eKhoj+HfjJUSKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImmKiqOnpwc2mw2ZmZmw2Wzo7e2dcExVVRWWLVuGnJwc5OTkwOFwePddu3YNr732GjIyMmCxWPDtt9+qdgJEpD3VxiMAU49A+OCDDxATE4NvvvkGvb29yM/PR1NTE6Kjo2/9DIhIczMyHuGfvvrqK9hsNgBAQkICkpOT0dLSEmBkIgo2VccjTDUCoa+vD/fdd5/3OKPRiAsXLkgF5XiEmaGXrHrJCegra6BUG4+gxgiE6XA8gvr0klUvOQH9ZA2Z8QjTjUAwmUw4f/6891iXy4W5c+cGHJqIgku18QjTjUCwWCyoq6sDAPT29uKXX37xmStLRPqi2niE6UYgvPTSS7Db7cjIyIDBYMCOHTsQExMzc2dFRDNKtfEI041AuPPOO1FZWRlAPCIKRfzkKBFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTRFy+p7enpgt9sxODiI2NhYOJ1OJCQk+BxTXV2NxsZG7/04SkpKvDfrsdvt+O6777y3EbRYLNi4caO6Z0JEmlFtPMKSJUtQWFiIqKgodHZ2oqCgAK2trbjjjjsAAMXFxSgoKFD/DIhIc6qNR0hNTUVUVBQAwGw2QwiBwcFB9RMTUdCpOh5hTH19PRYsWOBzQ+KamhrU1dVh/vz5eOONN7Bw4UKpoByPMDP0klUvOQF9ZQ2UauMRxpw4cQK7d+/Ghx9+6N1WUlKCuLg4GAwG1NfXo6ioCM3Nzd4yUoLjEdSnl6x6yQnoJ2vIjEcAgJMnT2Lz5s2orq5GYmKid3t8fDwMhtG/Kjc3F1evXpUeyEREoUO18Qjt7e0oKSlBZWUlFi9e7LNv/OiEo0ePwmAwID4+Xo38RBQEqo1HcDgcuH79OkpLS72Pq6iogNlsxpYtW9Df34+wsDDExMTg3XffxW23qf4qiYg0otp4hH379k35+NraWvlkRBSy+MlRIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaYqKo6enBzabDZmZmbDZbOjt7Z1wjMfjgcPhwIoVK5CRkeGzDH+6fUSkP4qKY2w8wtdff421a9f63KxnzMGDB3H27Fk0NTWhrq4OVVVV+Ouvv/zuIyL98Xsjn7HxCDU1NQBGxyOUl5djYGDA5/aBjY2NePbZZ2EwGDB79mysWLEChw8fRlFR0bT7/BFCAACGh4cDPUfN3bhxI9gRFNNLVr3kBPSRdex6Gru+ZKk2HsHlcsFkMnl/NhqN3hsST7fPH7fbDQA4c+aMouNDgZ5GOeglq15yAvrK6na7vUPTZIT8jT+jo6OxaNEiREREICwsLNhxiP4VhBBwu92Ijo4O6PF+i2P8eITw8PApxyMYjUb09fVhyZIlAHyfZUy3zx+DwYBZs2ZJnRQR+RfIM40xqo1HsFgs2Lt3L0ZGRjAwMIDm5mZkZmb63UdE+hMmFLw70t3dDbvdjitXrnjHIyQmJvqMR/B4PNixYweOHTsGYHR0gs1mA4Bp9xGR/igqDiKi8fjJUSKSxuIgImksDiKSxuIgImksDiKSFjLFcasrcLWkJGt1dTWeeuopWK1WrFmzBkePHtU+KJRlHfPHH3/goYcegtPp1C7gfyjN2djYCKvViuzsbFitVly+fFnboFCWtb+/H8XFxbBarcjKykJZWRlu3rypaU6n04m0tDSYzeYpl2wEfE2JELFu3TpRX18vhBCivr5erFu3bsIx+/fvF4WFhcLj8Yj+/n6Rmpoqzp07p3VURVlbWlrE1atXhRBCnD59WqSkpIhr165pmlMIZVmFEOLmzZuioKBAvP766+Ktt97SMqIQQlnO9vZ2kZWVJS5duiSEEOLKlSvi+vXrmuYUQlnWnTt3ev8dh4eHRV5enmhoaNA0548//ij6+vrE8uXLRVdX16THBHpNhcQzjrEVuNnZ2QBGV+B2dHRgYGDA57ipVtmGYtbU1FRERUUBAMxmM4QQGBwcDMmsAPDee+/hySefREJCgqYZAeU5a2trUVhYiLi4OADArFmzcPvtt4dk1rCwMAwNDWFkZATDw8Nwu92Ij4/XNOvSpUsnLA35p0CvqZAojulW4P7zuEBX2Wqddbz6+nosWLAAc+fO1SomAOVZOzs70draig0bNmiab4zSnN3d3Th37hzy8/Px9NNPY8+ePQEvC5/prJs2bUJPTw+eeOIJ71dKSoqmWZUI9JoKieL4Nztx4gR2796Nt99+O9hRJuV2u7Ft2zY4HA7vxRCqPB4Purq6UFNTg08++QQtLS04cOBAsGNN6vDhwzCbzWhtbUVLSwt++uknzZ8dz6SQKI7xK3AB+F2BO8blcmn+v7jSrABw8uRJbN68GdXV1UhMTNQ0J6As699//42zZ8+iuLgYaWlp+Oijj/DFF19g27ZtIZUTAEwmEywWCyIjIxETE4P09HS0t7drllMm66efforVq1d7V3enpaXh+PHjmmZVItBrKiSKQ40VuKGWtb29HSUlJaisrMTixYs1zThGSVaTyYTjx4/jyJEjOHLkCNavX4/nnnsO5eXlIZUTGH0/obW11XsviR9++AEPPvigZjllss6bNw8tLS0ARu+29f333+OBBx7QNKsSAV9Tqr6Newt+//13kZeXJ1auXCny8vJEd3e3EEKIoqIi0d7eLoQYfee/tLRUpKeni/T0dPH555+HbNY1a9aIRx99VKxevdr71dnZGZJZx6usrAzKb1WU5PR4PGLXrl3CYrGIVatWiV27dgmPxxOSWf/880+xYcMGkZ2dLbKyskRZWZlwu92a5iwvLxepqakiKSlJPPbYY2LVqlUTcgZ6TXF1LBFJC4mXKkSkLywOIpLG4iAiaX6L41Y/7x4q60uISD1+73Kenp6OF154Afn5+VMeM35S2+DgIHJzc7Fs2TLMmzdv2n1EpE9+i2Pp0qV+/5CZmuIGACMjIxgaGuJcFSIViXFzVQwG+XcsVBnINFNT3ABgaGhIV1PciPRk0aJFAc0tCvlJbhEREQBGTzAyMjLIafw7deoUkpOTgx1DEb1k1UtOQD9Zh4eHcebMGe/1JUuV4pipKW4AvC9PIiMjNV9CHSi95AT0k1UvOQF9ZQ305b8qv47lFDei/y5+n3Hs3LkTTU1NuHz5Ml588UXExsaioaHBZ4pbTk4Ofv75Z6xcuRIA8Morr2D+/PkAMO0+ItInv8WxdetWbN26dcL2999/3/t9eHg4HA7HpI+fbh8R6RM/OUpE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCRN0R3Aenp6YLfbMTg4iNjYWDidTiQkJPgc8+abb6Krq8v7c1dXF6qrq5Geno6qqip89tlnmDNnDgDgkUcewfbt29U7CyLSlKLi2L59O9auXYucnBwcOHAApaWl+Pjjj32Oqaio8H7f2dmJ9evXIzU11bstNzcXW7ZsUSk2EQWT35cq/f396OjoQHZ2NgAgOzsbHR0dGBgYmPIxX375JaxWqy5uLkxE8vw+43C5XIiPj0d4eDiA0Tt6zZkzBy6XC7Nnz55w/PDwMA4ePIja2lqf7Q0NDWhtbUVcXBxeffVVPPzww1JBT506JXV8MLW1tQU7gmJ6yaqXnIC+sgZK9fEIzc3NMJlMSEpK8m57/vnn8fLLLyMiIgLHjh3Dpk2b0NjYiHvuuUfxn5ucnKyLu0e3tbUhJSUl2DEU0UtWveQE9JP1xo0bt/Sfsd+XKkajERcvXoTH4wEwOgv20qVLMBqNkx6/b98+PPPMMz7b4uLivPMbHn/8cRiNRvz2228Bhyai4PJbHPfeey+SkpJw6NAhAMChQ4eQlJQ06cuUCxcuoK2tDVar1Wf7xYsXvd+fPn0a58+fx/3333+r2YkoSBS9VCkrK4PdbseePXtw1113wel0AoDPiAQA2L9/P5YvX467777b5/HvvPMOfv31VxgMBkRERKCiogJxcXEqnwoRaUVRcSxcuBB79+6dsH38iAQA2Lhx46SPHysaIvp34CdHiUgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpLE4iEgai4OIpCkqjp6eHthsNmRmZsJms6G3t3fCMVVVVVi2bBlycnKQk5MDh8Ph3Xft2jW89tpryMjIgMViwbfffqvaCRCR9lQbjwBMPQLhgw8+QExMDL755hv09vYiPz8fTU1NiI6OvvUzICLNzch4hH/66quvYLPZAAAJCQlITk5GS0tLgJGJKNhUHY8w1QiEvr4+3Hfffd7jjEYjLly4IBWU4xFmhl6y6iUnoK+sgVJtPIIaIxCmw/EI6tNLVr3kBPSTNWTGI0w3AsFkMuH8+fPeY10uF+bOnRtwaCIKLtXGI0w3AsFisaCurg4A0Nvbi19++cVnriwR6Ytq4xGmG4Hw0ksvwW63IyMjAwaDATt27EBMTMzMnRURzSjVxiNMNwLhzjvvRGVlZQDxiCgU8ZOjRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0lgcRCSNxUFE0hQtq+/p6YHdbsfg4CBiY2PhdDqRkJDgc0x1dTUaGxu99+MoKSnx3qzHbrfju+++895G0GKxYOPGjeqeCRFpRrXxCEuWLEFhYSGioqLQ2dmJgoICtLa24o477gAAFBcXo6CgQP0zICLNqTYeITU1FVFRUQAAs9kMIQQGBwfVT0xEQafqeIQx9fX1WLBggc8NiWtqalBXV4f58+fjjTfewMKFC6WCcjzCzNBLVr3kBPSVNVCqjUcYc+LECezevRsffvihd1tJSQni4uJgMBhQX1+PoqIiNDc3e8tICY5HUJ9esuolJ6CfrCEzHgEATp48ic2bN6O6uhqJiYne7fHx8TAYRv+q3NxcXL16VXogExGFDtXGI7S3t6OkpASVlZVYvHixz77xoxOOHj0Kg8GA+Ph4NfITURCoNh7B4XDg+vXrKC0t9T6uoqICZrMZW7ZsQX9/P8LCwhATE4N3330Xt92m+qskItKIauMR9u3bN+Xja2tr5ZMRUcjiJ0eJSBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4iksTiISBqLg4ikKSqOnp4e2Gw2ZGZmwmazobe3d8IxHo8HDocDK1asQEZGhs8y/On2EZH+KCqOsfEIX3/9NdauXetzs54xBw8exNmzZ9HU1IS6ujpUVVXhr7/+8ruPiPTH7418xsYj1NTUABgdj1BeXo6BgQGf2wc2Njbi2WefhcFgwOzZs7FixQocPnwYRUVF0+7zRwgBABgeHg70HDV348aNYEdQTC9Z9ZIT0EfWsetp7PqSpdp4BJfLBZPJ5P3ZaDR6b0g83T5/3G43AODMmTOKjg8FehrloJeseskJ6Cur2+32Dk2TEfI3/oyOjsaiRYsQERGBsLCwYMch+lcQQsDtdiM6Ojqgx/stjvHjEcLDw6ccj2A0GtHX14clS5YA8H2WMd0+fwwGA2bNmiV1UkTkXyDPNMaoNh7BYrFg7969GBkZwcDAAJqbm5GZmel3HxHpT5hQ8O5Id3c37HY7rly54h2PkJiY6DMewePxYMeOHTh27BiA0dEJNpsNAKbdR0T6o6g4iIjG4ydHiUgai4OIpLE4iEgai4OIpLE4iEhayBTHra7A1ZKSrNXV1XjqqadgtVqxZs0aHD16VPugUJZ1zB9//IGHHnoITqdTu4D/oTRnY2MjrFYrsrOzYbVacfnyZW2DQlnW/v5+FBcXw2q1IisrC2VlZbh586amOZ1OJ9LS0mA2m6dcshHwNSVCxLp160R9fb0QQoj6+nqxbt26Ccfs379fFBYWCo/HI/r7+0Vqaqo4d+6c1lEVZW1paRFXr14VQghx+vRpkZKSIq5du6ZpTiGUZRVCiJs3b4qCggLx+uuvi7feekvLiEIIZTnb29tFVlaWuHTpkhBCiCtXrojr169rmlMIZVl37tzp/XccHh4WeXl5oqGhQdOcP/74o+jr6xPLly8XXV1dkx4T6DUVEs84xlbgZmdnAxhdgdvR0YGBgQGf46ZaZRuKWVNTUxEVFQUAMJvNEEJgcHAwJLMCwHvvvYcnn3wSCQkJmmYElOesra1FYWEh4uLiAACzZs3C7bffHpJZw8LCMDQ0hJGREQwPD8PtdiM+Pl7TrEuXLp2wNOSfAr2mQqI4pluB+8/jAl1lq3XW8err67FgwQLMnTtXq5gAlGft7OxEa2srNmzYoGm+MUpzdnd349y5c8jPz8fTTz+NPXv2BLwsfKazbtq0CT09PXjiiSe8XykpKZpmVSLQayokiuPf7MSJE9i9ezfefvvtYEeZlNvtxrZt2+BwOLwXQ6jyeDzo6upCTU0NPvnkE7S0tODAgQPBjjWpw4cPw2w2o7W1FS0tLfjpp580f3Y8k0KiOMavwAXgdwXuGJfLpfn/4kqzAsDJkyexefNmVFdXIzExUdOcgLKsf//9N86ePYvi4mKkpaXho48+whdffIFt27aFVE4AMJlMsFgsiIyMRExMDNLT09He3q5ZTpmsn376KVavXu1d3Z2Wlobjx49rmlWJQK+pkCgONVbghlrW9vZ2lJSUoLKyEosXL9Y04xglWU0mE44fP44jR47gyJEjWL9+PZ577jmUl5eHVE5g9P2E1tZW770kfvjhBzz44IOa5ZTJOm/ePLS0tAAYvdvW999/jwceeEDTrEoEfE2p+jbuLfj9999FXl6eWLlypcjLyxPd3d1CCCGKiopEe3u7EGL0nf/S0lKRnp4u0tPTxeeffx6yWdesWSMeffRRsXr1au9XZ2dnSGYdr7KyMii/VVGS0+PxiF27dgmLxSJWrVoldu3aJTweT0hm/fPPP8WGDRtEdna2yMrKEmVlZcLtdmuas7y8XKSmpoqkpCTx2GOPiVWrVk3IGeg1xdWxRCQtJF6qEJG+sDiISBqLg4ik+S2OW/28e6isLyEi9fi9y3l6ejpeeOEF5OfnT3nM+Eltg4ODyM3NxbJlyzBv3rxp9xGRPvktjqVLl/r9Q2ZqihsAjIyMYGhoiHNViFQkxs1VMRjk37FQZSDTTE1xA4ChoSFdTXEj0pNFixYFNLco5Ce5RUREABg9wcjIyCCn8e/UqVNITk4OdgxF9JJVLzkB/WQdHh7GmTNnvNeXLFWKY6amuAHwvjyJjIzUfAl1oPSSE9BPVr3kBPSVNdCX/6r8OpZT3Ij+u/h9xrFz5040NTXh8uXLePHFFxEbG4uGhgafKW45OTn4+eefsXLlSgDAK6+8gvnz5wPAtPuISJ/8FsfWrVuxdevWCdvff/997/fh4eFwOByTPn66fUSkT/zkKBFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTRFdwDr6emB3W7H4OAgYmNj4XQ6kZCQ4HPMm2++ia6uLu/PXV1dqK6uRnp6OqqqqvDZZ59hzpw5AIBHHnkE27dvV+8siEhTiopj+/btWLt2LXJycnDgwAGUlpbi448/9jmmoqLC+31nZyfWr1+P1NRU77bc3Fxs2bJFpdhEFEx+X6r09/ejo6MD2dnZAIDs7Gx0dHRgYGBgysd8+eWXsFqturi5MBHJ8/uMw+VyIT4+HuHh4QBG7+g1Z84cuFwuzJ49e8Lxw8PDOHjwIGpra322NzQ0oLW1FXFxcXj11Vfx8MMPSwU9deqU1PHB1NbWFuwIiuklq15yAvrKGijVxyM0NzfDZDIhKSnJu+3555/Hyy+/jIiICBw7dgybNm1CY2Mj7rnnHsV/bnJysi7uHt3W1oaUlJRgx1BEL1n1khPQT9YbN27c0n/Gfl+qGI1GXLx4ER6PB8DoLNhLly7BaDROevy+ffvwzDPP+GyLi4vzzm94/PHHYTQa8dtvvwUcmoiCy29x3HvvvUhKSsKhQ4cAAIcOHUJSUtKkL1MuXLiAtrY2WK1Wn+0XL170fn/69GmcP38e999//61mJ6IgUfRSpaysDHa7HXv27MFdd90Fp9MJAD4jEgBg//79WL58Oe6++26fx7/zzjv49ddfYTAYEBERgYqKCsTFxal8KkSkFUXFsXDhQuzdu3fC9vEjEgBg48aNkz5+rGiI6N+BnxwlImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSxuIgImksDiKSpqg4enp6YLPZkJmZCZvNht7e3gnHVFVVYdmyZcjJyUFOTg4cDod337Vr1/Daa68hIyMDFosF3377rWonQETaU208AjD1CIQPPvgAMTEx+Oabb9Db24v8/Hw0NTUhOjr61s+AiDQ3I+MR/umrr76CzWYDACQkJCA5ORktLS0BRiaiYFN1PMJUIxD6+vpw3333eY8zGo24cOGCVFCOR5gZesmql5yAvrIGSrXxCGqMQJgOxyOoTy9Z9ZIT0E/WkBmPMN0IBJPJhPPnz3uPdblcmDt3bsChiSi4VBuPMN0IBIvFgrq6OgBAb28vfvnlF5+5skSkL6qNR5huBMJLL70Eu92OjIwMGAwG7NixAzExMTN3VkQ0o1QbjzDdCIQ777wTlZWVAcQjolDET44SkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJY3EQkTQWBxFJU7SsvqenB3a7HYODg4iNjYXT6URCQoLPMdXV1WhsbPTej6OkpMR7sx673Y7vvvvOextBi8WCjRs3qnsmRKQZ1cYjLFmyBIWFhYiKikJnZycKCgrQ2tqKO+64AwBQXFyMgoIC9c+AiDSn2niE1NRUREVFAQDMZjOEEBgcHFQ/MREFnarjEcbU19djwYIFPjckrqmpQV1dHebPn4833ngDCxculArK8QgzQy9Z9ZIT0FfWQKk2HmHMiRMnsHv3bnz44YfebSUlJYiLi4PBYEB9fT2KiorQ3NzsLSMlOB5BfXrJqpecgH6yhsx4BAA4efIkNm/ejOrqaiQmJnq3x8fHw2AY/atyc3Nx9epV6YFMRBQ6VBuP0N7ejpKSElRWVmLx4sU++8aPTjh69CgMBgPi4+PVyE9EQaDaeASHw4Hr16+jtLTU+7iKigqYzWZs2bIF/f39CAsLQ0xMDN59913cdpvqr5KISCOqjUfYt2/flI+vra2VT0ZEIYufHCUiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpKmqDh6enpgs9mQmZkJm82G3t7eCcd4PB44HA6sWLECGRkZPsvwp9tHRPqjqDjGxiN8/fXXWLt2rc/NesYcPHgQZ8+eRVNTE+rq6lBVVYW//vrL7z4i0h+/N/IZG49QU1MDYHQ8Qnl5OQYGBnxuH9jY2Ihnn30WBoMBs2fPxooVK3D48GEUFRVNu88fIQQAYHh4ONBz1NyNGzeCHUExvWTVS05AH1nHrqex60uWauMRXC4XTCaT92ej0ei9IfF0+/xxu90AgDNnzig6PhToaZSDXrLqJSegr6xut9s7NE1GyN/4Mzo6GosWLUJERATCwsKCHYfoX0EIAbfbjejo6IAe77c4xo9HCA8Pn3I8gtFoRF9fH5YsWQLA91nGdPv8MRgMmDVrltRJEZF/gTzTGKPaeASLxYK9e/diZGQEAwMDaG5uRmZmpt99RKQ/YULBuyPd3d2w2+24cuWKdzxCYmKiz3gEj8eDHTt24NixYwBGRyfYbDYAmHYfEemPouIgIhqPnxwlImksDiKSxuIgImksDiKSxuIgImkhUxy3ugJXS0qyVldX46mnnoLVasWaNWtw9OhR7YNCWdYxf/zxBx566CE4nU7tAv6H0pyNjY2wWq3Izs6G1WrF5cuXtQ0KZVn7+/tRXFwMq9WKrKwslJWV4ebNm5rmdDqdSEtLg9lsnnLJRsDXlAgR69atE/X19UIIIerr68W6desmHLN//35RWFgoPB6P6O/vF6mpqeLcuXNaR1WUtaWlRVy9elUIIcTp06dFSkqKuHbtmqY5hVCWVQghbt68KQoKCsTrr78u3nrrLS0jCiGU5WxvbxdZWVni0qVLQgghrly5Iq5fv65pTiGUZd25c6f333F4eFjk5eWJhoYGTXP++OOPoq+vTyxfvlx0dXVNekyg11RIPOMYW4GbnZ0NYHQFbkdHBwYGBnyOm2qVbShmTU1NRVRUFADAbDZDCIHBwcGQzAoA7733Hp588kkkJCRomhFQnrO2thaFhYWIi4sDAMyaNQu33357SGYNCwvD0NAQRkZGMDw8DLfbjfj4eE2zLl26dMLSkH8K9JoKieKYbgXuP48LdJWt1lnHq6+vx4IFCzB37lytYgJQnrWzsxOtra3YsGGDpvnGKM3Z3d2Nc+fOIT8/H08//TT27NkT8LLwmc66adMm9PT04IknnvB+paSkaJpViUCvqZAojn+zEydOYPfu3Xj77beDHWVSbrcb27Ztg8Ph8F4Mocrj8aCrqws1NTX45JNP0NLSggMHDgQ71qQOHz4Ms9mM1tZWtLS04KefftL82fFMConiGL8CF4DfFbhjXC6X5v+LK80KACdPnsTmzZtRXV2NxMRETXMCyrL+/fffOHv2LIqLi5GWloaPPvoIX3zxBbZt2xZSOQHAZDLBYrEgMjISMTExSE9PR3t7u2Y5ZbJ++umnWL16tXd1d1paGo4fP65pViUCvaZCojjUWIEbalnb29tRUlKCyspKLF68WNOMY5RkNZlMOH78OI4cOYIjR45g/fr1eO6551BeXh5SOYHR9xNaW1u995L44Ycf8OCDD2qWUybrvHnz0NLSAmD0blvff/89HnjgAU2zKhHwNaXq27i34Pfffxd5eXli5cqVIi8vT3R3dwshhCgqKhLt7e1CiNF3/ktLS0V6erpIT08Xn3/+echmXbNmjXj00UfF6tWrvV+dnZ0hmXW8ysrKoPxWRUlOj8cjdu3aJSwWi1i1apXYtWuX8Hg8IZn1zz//FBs2bBDZ2dkiKytLlJWVCbfbrWnO8vJykZqaKpKSksRjjz0mVq1aNSFnoNcUV8cSkbSQeKlCRPrC4iAiaSwOIpLG4iAiaSwOIpLG4iAiaSwOIpL2/w51aDv48ldZAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 595.276x841.89 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAAKfCAYAAAArTIf4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAru0lEQVR4nO3dT0jTAR/H8c+2Z1aIUYbaFMI6KCNPJQXGApuopKswxRIkEokSgjykRqKHPLRjih3rEEVR0UCzJA+h1KEIwUzTEoXwX6BIpKJrfp9DNOqpvr/V9tPF83lB0Ny/D/J7u6k/0CIiAiL6JetaDyCKZgyESMFAiBQMhEjBQIgUDIRIYRiI1+vFgQMHkJ6ejuHh4dXYRBQ1DANxu924efMmUlJSVmMPUVT5j9ENMjMzV2MHUVQyDCQcKysrmJ+fh91uh8ViMfOpiAAAIgK/34/Y2FhYreF/i21qIPPz8/y+hdZEWloa4uLiwn4cUwOx2+0Avo6NiYkx86n+Sn9/PzIyMtZ6xm9F875o3ba8vIzh4eHgsRcuUwP59rYqJiYG69atM/Op/lq07vommvdF87ZIvaU3fJPW1NSE/fv3Y2pqCidPnkRBQUFEnpjoX2D4ClJfX4/6+vrV2EIUdfibdCIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUoQUyOjoKEpLS5GXl4fS0lKMjY2ZPIsoOoQUSGNjI8rKytDZ2YmysjI0NDSYvYsoKhgGMjMzg4GBARQWFgIACgsLMTAwgNnZWdPHEa01w79yOzk5iaSkJNhsNgCAzWZDYmIiJicnER8fr95XRAB8/ePu0WppaWmtJ6iieV80bvt2rH079sJlGEg4/H4/AGB4eNjMpwlLf3//Wk9QRfO+aN7m9/uxfv36sB/HMBCHw4Hp6WkEAgHYbDYEAgF8/PgRDofD8MFjY2ORlpYGu90Oi8US9lgiIyICv9+P2NjYiDyeYSBbtmyB0+lEe3s7Dh8+jPb2djidTsO3VwBgtVoRFxcXkaFEoYrEK8c3FgnhzdrIyAjq6urw6dMnbNy4EV6vFzt27IjYCKJoFVIgRP+v+Jt0IgUDIVIwECIFAyFSMBAixV/9Jn10dBR1dXWYm5vDpk2b4PV6kZqa+sNtAoEAmpqa0NPTA4vFglOnTqGkpMTwukgIZV9rays6OjpgtVpht9tRXV0Nl8sFAKirq8Pz58+xefNmAEB+fj7OnDmzattaWlpw69YtJCYmAgB27dqFxsZGAMDi4iIuXLiAN2/ewGazoba2FtnZ2RHZFuq+mpoaDA0NBS8PDQ2htbUVbrdb3R4Or9eLzs5OjI+Po62tDWlpaT/dxpRjTv5CeXm5+Hw+ERHx+XxSXl7+020ePHggFRUVEggEZGZmRlwul3z48MHwukgIZV93d7csLCyIiMjg4KDs3r1bFhcXRUSktrZWbty4EbE9f7qtublZLl++/Mv7t7S0yMWLF0VEZHR0VLKysuTz58+ruu97g4ODsmfPHllaWjLcHo6XL1/KxMSEZGdny9DQ0C9vY8Yx98dvsUI9u7ejowMlJSWwWq2Ij49HTk4OHj9+bHhduELd53K5sGHDBgBAeno6RARzc3MR2RDuNs2jR49QWloKAEhNTUVGRga6u7vXbN+9e/fg8XgQExMTkQ2/k5mZaXh6kxnH3B8Hop3d+7+3S05ODl52OByYmpoyvC5coe77ns/nw7Zt27B169bgx65fvw6Px4OqqiqMjIys+raHDx/C4/GgoqICvb29wY9PTEwgJSUleHktP3fLy8toa2vD0aNHQ9puNjOOOVPP5v0XvHjxAleuXMG1a9eCH6uurkZCQgKsVit8Ph8qKyvR1dUVPHDMduzYMZw+fRp2ux3Pnj1DVVUVOjo6gt8TRYuuri4kJyfD6XQGP/avbA/VH7+CfH92L4Dfnt3rcDgwMTERvDw5ORn8Cq1dF65Q9wFAb28vzp8/j9bW1h/OLUtKSoLV+vVTc+TIESwsLETkq3So2xISEmC32wEA+/btg8PhwLt37wAAycnJGB8fD952rT53AHD//v2fXj207WYz45j740C+P7sXwG/P7s3Pz8fdu3exsrKC2dlZdHV1IS8vz/C6cIW6r6+vD9XV1WhubsbOnTt/uG56ejr4/56eHlitViQlJa3atu+ff3BwEOPj49i+fTuAr5+7O3fuAADGxsbw+vXr4E/fVmsfAExNTeHVq1fweDwhbzebKcfc3/xE4f3791JcXCy5ublSXFwsIyMjIiJSWVkpfX19IiLy5csXaWhoELfbLW63W27fvh28v3ZdJISyr6ioSPbu3SuHDh0K/nv79q2IiJw4cUIKCwvF4/HI8ePHpbe3d1W31dTUSEFBgXg8HikqKpKnT58G7z8/Py9nz56VnJwcyc3NlSdPnkRsW6j7RESuXr0q586d++n+2vZwXLp0SVwulzidTsnKypKDBw/+tMuMY45n8xIp+Jt0IgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSGAbi9Xpx4MABpKenY3h4eDU2EUUNw0Dcbjdu3ryJlJSU1dhDFFX+Y3SDzMzM1dhBFJUMAwnHysoK5ufnYbfbYbFYzHwqIgCAiMDv9yM2NhZWa/jfYpsayPz8PL9voTWRlpaGuLi4sB/H1EDsdjuAr2NjYmLMfKq/0t/fj4yMjLWe8VvRvC9aty0vL2N4eDh47IXL1EC+va2KiYnBunXrzHyqvxatu76J5n3RvC1Sb+kN36Q1NTVh//79mJqawsmTJ1FQUBCRJyb6Fxi+gtTX16O+vn41thBFHf4mnUjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFCEFMjo6itLSUuTl5aG0tBRjY2MmzyKKDiEF0tjYiLKyMnR2dqKsrAwNDQ1m7yKKCoaBzMzMYGBgAIWFhQCAwsJCDAwMYHZ21vRxRGvN8K/cTk5OIikpCTabDQBgs9mQmJiIyclJxMfHq/cVEQBf/7h7tFpaWlrrCapo3heN274da9+OvXAZBhIOv98PABgeHjbzacLS39+/1hNU0bwvmrf5/X6sX78+7McxDMThcGB6ehqBQAA2mw2BQAAfP36Ew+EwfPDY2FikpaXBbrfDYrGEPZbIiIjA7/cjNjY2Io9nGMiWLVvgdDrR3t6Ow4cPo729HU6n0/DtFQBYrVbExcVFZChRqCLxyvGNRUJ4szYyMoK6ujp8+vQJGzduhNfrxY4dOyI2gihahRQI0f8r/iadSMFAiBQMhEjBQIgUDIRI8Ve/SR8dHUVdXR3m5uawadMmeL1epKam/nCbQCCApqYm9PT0wGKx4NSpUygpKTG8LhJC2dfa2oqOjg5YrVbY7XZUV1fD5XIBAOrq6vD8+XNs3rwZAJCfn48zZ86s2raWlhbcunULiYmJAIBdu3ahsbERALC4uIgLFy7gzZs3sNlsqK2tRXZ2dkS2hbqvpqYGQ0NDwctDQ0NobW2F2+1Wt4fD6/Wis7MT4+PjaGtrQ1pa2k+3MeWYk79QXl4uPp9PRER8Pp+Ul5f/dJsHDx5IRUWFBAIBmZmZEZfLJR8+fDC8LhJC2dfd3S0LCwsiIjI4OCi7d++WxcVFERGpra2VGzduRGzPn25rbm6Wy5cv//L+LS0tcvHiRRERGR0dlaysLPn8+fOq7vve4OCg7NmzR5aWlgy3h+Ply5cyMTEh2dnZMjQ09MvbmHHM/fFbrFDP7u3o6EBJSQmsVivi4+ORk5ODx48fG14XrlD3uVwubNiwAQCQnp4OEcHc3FxENoS7TfPo0SOUlpYCAFJTU5GRkYHu7u4123fv3j14PB7ExMREZMPvZGZmGp7eZMYx98eBaGf3/u/tkpOTg5cdDgempqYMrwtXqPu+5/P5sG3bNmzdujX4sevXr8Pj8aCqqgojIyOrvu3hw4fweDyoqKhAb29v8OMTExNISUkJXl7Lz93y8jLa2tpw9OjRkLabzYxjztSzef8FL168wJUrV3Dt2rXgx6qrq5GQkACr1Qqfz4fKykp0dXUFDxyzHTt2DKdPn4bdbsezZ89QVVWFjo6O4PdE0aKrqwvJyclwOp3Bj/0r20P1x68g35/dC+C3Z/c6HA5MTEwEL09OTga/QmvXhSvUfQDQ29uL8+fPo7W19Ydzy5KSkmC1fv3UHDlyBAsLCxH5Kh3qtoSEBNjtdgDAvn374HA48O7dOwBAcnIyxsfHg7ddq88dANy/f/+nVw9tu9nMOOb+OJDvz+4F8Nuze/Pz83H37l2srKxgdnYWXV1dyMvLM7wuXKHu6+vrQ3V1NZqbm7Fz584frpueng7+v6enB1arFUlJSau27fvnHxwcxPj4OLZv3w7g6+fuzp07AICxsTG8fv06+NO31doHAFNTU3j16hU8Hk/I281myjH3Nz9ReP/+vRQXF0tubq4UFxfLyMiIiIhUVlZKX1+fiIh8+fJFGhoaxO12i9vtltu3bwfvr10XCaHsKyoqkr1798qhQ4eC/96+fSsiIidOnJDCwkLxeDxy/Phx6e3tXdVtNTU1UlBQIB6PR4qKiuTp06fB+8/Pz8vZs2clJydHcnNz5cmTJxHbFuo+EZGrV6/KuXPnfrq/tj0cly5dEpfLJU6nU7KysuTgwYM/7TLjmOPZvEQK/iadSMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUhoF4vV4cOHAA6enpGB4eXo1NRFHDMBC3242bN28iJSVlNfYQRZX/GN0gMzNzNXYQRSXDQMKxsrKC+fl52O12WCwWM5+KCAAgIvD7/YiNjYXVGv632KYGMj8/z+9baE2kpaUhLi4u7McxNRC73Q7g69iYmBgzn+qv9Pf3IyMjY61n/FY074vWbcvLyxgeHg4ee+EyNZBvb6tiYmKwbt06M5/qr0Xrrm+ieV80b4vUW3rDN2lNTU3Yv38/pqamcPLkSRQUFETkiYn+BYavIPX19aivr1+NLURRh79JJ1IwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiRUiBjI6OorS0FHl5eSgtLcXY2JjJs4iiQ0iBNDY2oqysDJ2dnSgrK0NDQ4PZu4iigmEgMzMzGBgYQGFhIQCgsLAQAwMDmJ2dNX0c0Voz/Cu3k5OTSEpKgs1mAwDYbDYkJiZicnIS8fHx6n1FBMDXP+4erZaWltZ6giqa90Xjtm/H2rdjL1yGgYTD7/cDAIaHh818mrD09/ev9QRVNO+L5m1+vx/r168P+3EMA3E4HJienkYgEIDNZkMgEMDHjx/hcDgMHzw2NhZpaWmw2+2wWCxhjyUyIiLw+/2IjY2NyOMZBrJlyxY4nU60t7fj8OHDaG9vh9PpNHx7BQBWqxVxcXERGUoUqki8cnxjkRDerI2MjKCurg6fPn3Cxo0b4fV6sWPHjoiNIIpWIQVC9P+Kv0knUjAQIgUDIVIwECIFAyFS/NVv0kdHR1FXV4e5uTls2rQJXq8XqampP9wmEAigqakJPT09sFgsOHXqFEpKSgyvi4RQ9rW2tqKjowNWqxV2ux3V1dVwuVwAgLq6Ojx//hybN28GAOTn5+PMmTOrtq2lpQW3bt1CYmIiAGDXrl1obGwEACwuLuLChQt48+YNbDYbamtrkZ2dHZFtoe6rqanB0NBQ8PLQ0BBaW1vhdrvV7eHwer3o7OzE+Pg42trakJaW9tNtTDnm5C+Ul5eLz+cTERGfzyfl5eU/3ebBgwdSUVEhgUBAZmZmxOVyyYcPHwyvi4RQ9nV3d8vCwoKIiAwODsru3btlcXFRRERqa2vlxo0bEdvzp9uam5vl8uXLv7x/S0uLXLx4UURERkdHJSsrSz5//ryq+743ODgoe/bskaWlJcPt4Xj58qVMTExIdna2DA0N/fI2Zhxzf/wWK9Szezs6OlBSUgKr1Yr4+Hjk5OTg8ePHhteFK9R9LpcLGzZsAACkp6dDRDA3NxeRDeFu0zx69AilpaUAgNTUVGRkZKC7u3vN9t27dw8ejwcxMTER2fA7mZmZhqc3mXHM/XEg2tm9/3u75OTk4GWHw4GpqSnD68IV6r7v+Xw+bNu2DVu3bg1+7Pr16/B4PKiqqsLIyMiqb3v48CE8Hg8qKirQ29sb/PjExARSUlKCl9fyc7e8vIy2tjYcPXo0pO1mM+OYM/Vs3n/BixcvcOXKFVy7di34serqaiQkJMBqtcLn86GyshJdXV3BA8dsx44dw+nTp2G32/Hs2TNUVVWho6Mj+D1RtOjq6kJycjKcTmfwY//K9lD98SvI92f3Avjt2b0OhwMTExPBy5OTk8Gv0Np14Qp1HwD09vbi/PnzaG1t/eHcsqSkJFitXz81R44cwcLCQkS+Soe6LSEhAXa7HQCwb98+OBwOvHv3DgCQnJyM8fHx4G3X6nMHAPfv3//p1UPbbjYzjrk/DuT7s3sB/Pbs3vz8fNy9excrKyuYnZ1FV1cX8vLyDK8LV6j7+vr6UF1djebmZuzcufOH66anp4P/7+npgdVqRVJS0qpt+/75BwcHMT4+ju3btwP4+rm7c+cOAGBsbAyvX78O/vRttfYBwNTUFF69egWPxxPydrOZcsz9zU8U3r9/L8XFxZKbmyvFxcUyMjIiIiKVlZXS19cnIiJfvnyRhoYGcbvd4na75fbt28H7a9dFQij7ioqKZO/evXLo0KHgv7dv34qIyIkTJ6SwsFA8Ho8cP35cent7V3VbTU2NFBQUiMfjkaKiInn69Gnw/vPz83L27FnJycmR3NxcefLkScS2hbpPROTq1aty7ty5n+6vbQ/HpUuXxOVyidPplKysLDl48OBPu8w45ng2L5GCv0knUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECKFYSBerxcHDhxAeno6hoeHV2MTUdQwDMTtduPmzZtISUlZjT1EUeU/RjfIzMxcjR1EUckwkHCsrKxgfn4edrsdFovFzKciAgCICPx+P2JjY2G1hv8ttqmBzM/P8/sWWhNpaWmIi4sL+3FMDcRutwP4OjYmJsbMp/or/f39yMjIWOsZvxXN+6J12/LyMoaHh4PHXrhMDeTb26qYmBisW7fOzKf6a9G665to3hfN2yL1lt7wTVpTUxP279+PqakpnDx5EgUFBRF5YqJ/geErSH19Perr61djC1HU4W/SiRQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIEVIgo6OjKC0tRV5eHkpLSzE2NmbyLKLoEFIgjY2NKCsrQ2dnJ8rKytDQ0GD2LqKoYBjIzMwMBgYGUFhYCAAoLCzEwMAAZmdnTR9HtNYM/8rt5OQkkpKSYLPZAAA2mw2JiYmYnJxEfHy8el8RAfD1j7tHq6WlpbWeoIrmfdG47dux9u3YC5dhIOHw+/0AgOHhYTOfJiz9/f1rPUEVzfuieZvf78f69evDfhzDQBwOB6anpxEIBGCz2RAIBPDx40c4HA7DB4+NjUVaWhrsdjssFkvYY4mMiAj8fj9iY2Mj8niGgWzZsgVOpxPt7e04fPgw2tvb4XQ6Dd9eAYDVakVcXFxEhhKFKhKvHN9YJIQ3ayMjI6irq8OnT5+wceNGeL1e7NixI2IjiKJVSIEQ/b/ib9KJFAyESMFAiBQMhEjBQIgUf/Wb9NHRUdTV1WFubg6bNm2C1+tFamrqD7cJBAJoampCT08PLBYLTp06hZKSEsPrIiGUfa2trejo6IDVaoXdbkd1dTVcLhcAoK6uDs+fP8fmzZsBAPn5+Thz5syqbWtpacGtW7eQmJgIANi1axcaGxsBAIuLi7hw4QLevHkDm82G2tpaZGdnR2RbqPtqamowNDQUvDw0NITW1la43W51ezi8Xi86OzsxPj6OtrY2pKWl/XQbU445+Qvl5eXi8/lERMTn80l5eflPt3nw4IFUVFRIIBCQmZkZcblc8uHDB8PrIiGUfd3d3bKwsCAiIoODg7J7925ZXFwUEZHa2lq5ceNGxPb86bbm5ma5fPnyL+/f0tIiFy9eFBGR0dFRycrKks+fP6/qvu8NDg7Knj17ZGlpyXB7OF6+fCkTExOSnZ0tQ0NDv7yNGcfcH7/FCvXs3o6ODpSUlMBqtSI+Ph45OTl4/Pix4XXhCnWfy+XChg0bAADp6ekQEczNzUVkQ7jbNI8ePUJpaSkAIDU1FRkZGeju7l6zfffu3YPH40FMTExENvxOZmam4elNZhxzfxyIdnbv/94uOTk5eNnhcGBqasrwunCFuu97Pp8P27Ztw9atW4Mfu379OjweD6qqqjAyMrLq2x4+fAiPx4OKigr09vYGPz4xMYGUlJTg5bX83C0vL6OtrQ1Hjx4NabvZzDjmTD2b91/w4sULXLlyBdeuXQt+rLq6GgkJCbBarfD5fKisrERXV1fwwDHbsWPHcPr0adjtdjx79gxVVVXo6OgIfk8ULbq6upCcnAyn0xn82L+yPVR//Ary/dm9AH57dq/D4cDExETw8uTkZPArtHZduELdBwC9vb04f/48Wltbfzi3LCkpCVbr10/NkSNHsLCwEJGv0qFuS0hIgN1uBwDs27cPDocD7969AwAkJydjfHw8eNu1+twBwP3793969dC2m82MY+6PA/n+7F4Avz27Nz8/H3fv3sXKygpmZ2fR1dWFvLw8w+vCFeq+vr4+VFdXo7m5GTt37vzhuunp6eD/e3p6YLVakZSUtGrbvn/+wcFBjI+PY/v27QC+fu7u3LkDABgbG8Pr16+DP31brX0AMDU1hVevXsHj8YS83WymHHN/8xOF9+/fS3FxseTm5kpxcbGMjIyIiEhlZaX09fWJiMiXL1+koaFB3G63uN1uuX37dvD+2nWREMq+oqIi2bt3rxw6dCj47+3btyIicuLECSksLBSPxyPHjx+X3t7eVd1WU1MjBQUF4vF4pKioSJ4+fRq8//z8vJw9e1ZycnIkNzdXnjx5ErFtoe4TEbl69aqcO3fup/tr28Nx6dIlcblc4nQ6JSsrSw4ePPjTLjOOOZ7NS6Tgb9KJFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEhhGIjX68WBAweQnp6O4eHh1dhEFDUMA3G73bh58yZSUlJWYw9RVPmP0Q0yMzNXYwdRVDIMJBwrKyuYn5+H3W6HxWIx86mIAAAiAr/fj9jYWFit4X+LbWog8/Pz/L6F1kRaWhri4uLCfhxTA7Hb7QC+jo2JiTHzqf5Kf38/MjIy1nrGb0Xzvmjdtry8jOHh4eCxFy5TA/n2tiomJgbr1q0z86n+WrTu+iaa90Xztki9pTd8k9bU1IT9+/djamoKJ0+eREFBQUSemOhfYPgKUl9fj/r6+tXYQhR1+Jt0IgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFShBTI6OgoSktLkZeXh9LSUoyNjZk8iyg6hBRIY2MjysrK0NnZibKyMjQ0NJi9iygqGAYyMzODgYEBFBYWAgAKCwsxMDCA2dlZ08cRrTXDv3I7OTmJpKQk2Gw2AIDNZkNiYiImJycRHx+v3ldEAHz94+7Ramlpaa0nqKJ5XzRu+3asfTv2wmUYSDj8fj8AYHh42MynCUt/f/9aT1BF875o3ub3+7F+/fqwH8cwEIfDgenpaQQCAdhsNgQCAXz8+BEOh8PwwWNjY5GWlga73Q6LxRL2WCIjIgK/34/Y2NiIPJ5hIFu2bIHT6UR7ezsOHz6M9vZ2OJ1Ow7dXAGC1WhEXFxeRoUShisQrxzcWCeHN2sjICOrq6vDp0yds3LgRXq8XO3bsiNgIomgVUiBE/6/4m3QiBQMhUjAQIgUDIVIwECLFX/0mfXR0FHV1dZibm8OmTZvg9XqRmpr6w20CgQCamprQ09MDi8WCU6dOoaSkxPC6SAhlX2trKzo6OmC1WmG321FdXQ2XywUAqKurw/Pnz7F582YAQH5+Ps6cObNq21paWnDr1i0kJiYCAHbt2oXGxkYAwOLiIi5cuIA3b97AZrOhtrYW2dnZEdkW6r6amhoMDQ0FLw8NDaG1tRVut1vdHg6v14vOzk6Mj4+jra0NaWlpP93GlGNO/kJ5ebn4fD4REfH5fFJeXv7TbR48eCAVFRUSCARkZmZGXC6XfPjwwfC6SAhlX3d3tywsLIiIyODgoOzevVsWFxdFRKS2tlZu3LgRsT1/uq25uVkuX778y/u3tLTIxYsXRURkdHRUsrKy5PPnz6u673uDg4OyZ88eWVpaMtwejpcvX8rExIRkZ2fL0NDQL29jxjH3x2+xQj27t6OjAyUlJbBarYiPj0dOTg4eP35seF24Qt3ncrmwYcMGAEB6ejpEBHNzcxHZEO42zaNHj1BaWgoASE1NRUZGBrq7u9ds37179+DxeBATExORDb+TmZlpeHqTGcfcHweind37v7dLTk4OXnY4HJiamjK8Llyh7vuez+fDtm3bsHXr1uDHrl+/Do/Hg6qqKoyMjKz6tocPH8Lj8aCiogK9vb3Bj09MTCAlJSV4eS0/d8vLy2hra8PRo0dD2m42M445U8/m/Re8ePECV65cwbVr14Ifq66uRkJCAqxWK3w+HyorK9HV1RU8cMx27NgxnD59Gna7Hc+ePUNVVRU6OjqC3xNFi66uLiQnJ8PpdAY/9q9sD9Ufv4J8f3YvgN+e3etwODAxMRG8PDk5GfwKrV0XrlD3AUBvby/Onz+P1tbWH84tS0pKgtX69VNz5MgRLCwsROSrdKjbEhISYLfbAQD79u2Dw+HAu3fvAADJyckYHx8P3natPncAcP/+/Z9ePbTtZjPjmPvjQL4/uxfAb8/uzc/Px927d7GysoLZ2Vl0dXUhLy/P8Lpwhbqvr68P1dXVaG5uxs6dO3+4bnp6Ovj/np4eWK1WJCUlrdq2759/cHAQ4+Pj2L59O4Cvn7s7d+4AAMbGxvD69evgT99Wax8ATE1N4dWrV/B4PCFvN5spx9zf/ETh/fv3UlxcLLm5uVJcXCwjIyMiIlJZWSl9fX0iIvLlyxdpaGgQt9stbrdbbt++Hby/dl0khLKvqKhI9u7dK4cOHQr+e/v2rYiInDhxQgoLC8Xj8cjx48elt7d3VbfV1NRIQUGBeDweKSoqkqdPnwbvPz8/L2fPnpWcnBzJzc2VJ0+eRGxbqPtERK5evSrnzp376f7a9nBcunRJXC6XOJ1OycrKkoMHD/60y4xjjmfzEin4m3QiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIwECIFAyFSMBAiBQMhUjAQIgUDIVIYBuL1enHgwAGkp6djeHh4NTYRRQ3DQNxuN27evImUlJTV2EMUVf5jdIPMzMzV2EEUlQwDCcfKygrm5+dht9thsVjMfCoiAICIwO/3IzY2FlZr+N9imxrI/Pw8v2+hNZGWloa4uLiwH8fUQOx2O4CvY2NiYsx8qr/S39+PjIyMtZ7xW9G8L1q3LS8vY3h4OHjshcvUQL69rYqJicG6devMfKq/Fq27vonmfdG8LVJv6Q3fpDU1NWH//v2YmprCyZMnUVBQEJEnJvoXGL6C1NfXo76+fjW2EEUd/iadSMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFAyESMFAiBQMhEjBQIgUIQUyOjqK0tJS5OXlobS0FGNjYybPIooOIQXS2NiIsrIydHZ2oqysDA0NDWbvIooKhoHMzMxgYGAAhYWFAIDCwkIMDAxgdnbW9HFEa83wr9xOTk4iKSkJNpsNAGCz2ZCYmIjJyUnEx8er9xURAF//uHu0WlpaWusJqmjeF43bvh1r3469cBkGEg6/3w8AGB4eNvNpwtLf37/WE1TRvC+at/n9fqxfvz7sxzEMxOFwYHp6GoFAADabDYFAAB8/foTD4TB88NjYWKSlpcFut8NisYQ9lsiIiMDv9yM2NjYij2cYyJYtW+B0OtHe3o7Dhw+jvb0dTqfT8O0VAFitVsTFxUVkKFGoIvHK8Y1FQnizNjIygrq6Onz69AkbN26E1+vFjh07IjaCKFqFFAjR/yv+Jp1IwUCIFAyESMFAiBQMhEjxV79JHx0dRV1dHebm5rBp0yZ4vV6kpqb+cJtAIICmpib09PTAYrHg1KlTKCkpMbwuEkLZ19raio6ODlitVtjtdlRXV8PlcgEA6urq8Pz5c2zevBkAkJ+fjzNnzqzatpaWFty6dQuJiYkAgF27dqGxsREAsLi4iAsXLuDNmzew2Wyora1FdnZ2RLaFuq+mpgZDQ0PBy0NDQ2htbYXb7Va3h8Pr9aKzsxPj4+Noa2tDWlraT7cx5ZiTv1BeXi4+n09ERHw+n5SXl/90mwcPHkhFRYUEAgGZmZkRl8slHz58MLwuEkLZ193dLQsLCyIiMjg4KLt375bFxUUREamtrZUbN25EbM+fbmtubpbLly//8v4tLS1y8eJFEREZHR2VrKws+fz586ru+97g4KDs2bNHlpaWDLeH4+XLlzIxMSHZ2dkyNDT0y9uYccz98VusUM/u7ejoQElJCaxWK+Lj45GTk4PHjx8bXheuUPe5XC5s2LABAJCeng4RwdzcXEQ2hLtN8+jRI5SWlgIAUlNTkZGRge7u7jXbd+/ePXg8HsTExERkw+9kZmYant5kxjH3x4FoZ/f+7+2Sk5ODlx0OB6ampgyvC1eo+77n8/mwbds2bN26Nfix69evw+PxoKqqCiMjI6u+7eHDh/B4PKioqEBvb2/w4xMTE0hJSQleXsvP3fLyMtra2nD06NGQtpvNjGPO1LN5/wUvXrzAlStXcO3ateDHqqurkZCQAKvVCp/Ph8rKSnR1dQUPHLMdO3YMp0+fht1ux7Nnz1BVVYWOjo7g90TRoqurC8nJyXA6ncGP/SvbQ/XHryDfn90L4Ldn9zocDkxMTAQvT05OBr9Ca9eFK9R9ANDb24vz58+jtbX1h3PLkpKSYLV+/dQcOXIECwsLEfkqHeq2hIQE2O12AMC+ffvgcDjw7t07AEBycjLGx8eDt12rzx0A3L9//6dXD2272cw45v44kO/P7gXw27N78/PzcffuXaysrGB2dhZdXV3Iy8szvC5coe7r6+tDdXU1mpubsXPnzh+um56eDv6/p6cHVqsVSUlJq7bt++cfHBzE+Pg4tm/fDuDr5+7OnTsAgLGxMbx+/Tr407fV2gcAU1NTePXqFTweT8jbzWbKMfc3P1F4//69FBcXS25urhQXF8vIyIiIiFRWVkpfX5+IiHz58kUaGhrE7XaL2+2W27dvB++vXRcJoewrKiqSvXv3yqFDh4L/3r59KyIiJ06ckMLCQvF4PHL8+HHp7e1d1W01NTVSUFAgHo9HioqK5OnTp8H7z8/Py9mzZyUnJ0dyc3PlyZMnEdsW6j4RkatXr8q5c+d+ur+2PRyXLl0Sl8slTqdTsrKy5ODBgz/tMuOY49m8RAr+Jp1IwUCIFAyESMFAiBQMhEjBQIgUDIRIwUCIFP8FS66PYT8rtRoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 595.276x841.89 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "layout = fifi.svg_to_axes.FigureLayout('../../../../Figure/motionAnalysis.svg',\n",
    "                                           autogenlayers=True, make_mplfigures=True, hide_layers=[])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "695aae86",
   "metadata": {},
   "source": [
    "## HWS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0c09fe02",
   "metadata": {},
   "outputs": [],
   "source": [
    "dat = pd.DataFrame() # Replace this with your actual data\n",
    "dat=df1\n",
    "heading='HWS'\n",
    "# Your specific points and velocities\n",
    "points = dat[['xsrc', 'ysrc']].values\n",
    "velocities = dat[['gps_linear_x', 'gps_linear_y']].values\n",
    "\n",
    "# Calculate the rotation angle for each point\n",
    "angles = np.arctan2(points[:, 1], points[:, 0])\n",
    "\n",
    "rotated_velocities = []\n",
    "\n",
    "# Rotate each velocity to the body-fixed frame\n",
    "for angle, velocity in zip(angles, velocities):\n",
    "    rotation_matrix = calculate_rotation_matrix(angle)\n",
    "    rotated_velocity = rotate_vector(velocity, rotation_matrix)\n",
    "    rotated_velocities.append(rotated_velocity)\n",
    "\n",
    "# Convert the list of rotated velocities to a numpy array\n",
    "rotated_velocities = np.array(rotated_velocities)\n",
    "\n",
    "# Extract the components of the velocities in the body-fixed frame\n",
    "parallel_components = rotated_velocities[:, 0]\n",
    "perpendicular_components = rotated_velocities[:, 1]\n",
    "\n",
    "# Add the results to the DataFrame\n",
    "dat['parallel'] = parallel_components\n",
    "dat['perpendicular'] = perpendicular_components\n",
    "\n",
    "dt = dat.time[1] - dat.time[0]\n",
    "\n",
    "f_parallel, parallel_fft = fft_of_data(np.array(dat.parallel), dt)\n",
    "f_perpendicular, perpendicular_fft = fft_of_data(np.array(dat.perpendicular), dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "47baa492",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig2','t1')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(dat.time,dat.parallel, linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,2000)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,2000],\n",
    "                            yticks = [-5.0,0,5.0],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "# ax.set_yticklabels([])\n",
    "\n",
    "ax.set_ylabel('Parallel, m/s')\n",
    "ax.yaxis.set_label_coords(-0.12, 0.5)\n",
    "\n",
    "ax.set_xlabel('2000 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig2'], 'fig2', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d3511972",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig1','a1')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(f_parallel, np.abs(parallel_fft), linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,1],\n",
    "                            yticks = [0,0.3],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "# ax.set_xlabel('1000 secs')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "ax.set_ylabel('Parallel \\nAmplitude, m/s')\n",
    "ax.yaxis.set_label_coords(-0.04, 0.5)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig1'], 'fig1', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ccbc6d3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig2','t2')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(dat.time,dat.perpendicular, linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,2000)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,2000],\n",
    "                            yticks = [-5.0,0,5.0],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "# ax.set_yticklabels([])\n",
    "\n",
    "ax.set_ylabel('Perpendicular, m/s')\n",
    "ax.yaxis.set_label_coords(-0.12, 0.5)\n",
    "ax.set_xlabel('2000 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig2'], 'fig2', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "43a8376d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig1','a2')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(f_perpendicular, np.abs(perpendicular_fft), linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,1],\n",
    "                            yticks = [0,0.3],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "# ax.set_xlabel('1000 secs')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "ax.set_ylabel('Perpendicular \\nAmplitude, m/s')\n",
    "ax.yaxis.set_label_coords(-0.04, 0.5)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig1'], 'fig1', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "600f64f2",
   "metadata": {},
   "source": [
    "## LWS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c1a08c49",
   "metadata": {},
   "outputs": [],
   "source": [
    "dat = pd.DataFrame() # Replace this with your actual data\n",
    "dat=df2\n",
    "heading='LWS'\n",
    "# Your specific points and velocities\n",
    "points = dat[['xsrc', 'ysrc']].values\n",
    "velocities = dat[['gps_linear_x', 'gps_linear_y']].values\n",
    "\n",
    "# Calculate the rotation angle for each point\n",
    "angles = np.arctan2(points[:, 1], points[:, 0])\n",
    "\n",
    "rotated_velocities = []\n",
    "\n",
    "# Rotate each velocity to the body-fixed frame\n",
    "for angle, velocity in zip(angles, velocities):\n",
    "    rotation_matrix = calculate_rotation_matrix(angle)\n",
    "    rotated_velocity = rotate_vector(velocity, rotation_matrix)\n",
    "    rotated_velocities.append(rotated_velocity)\n",
    "\n",
    "# Convert the list of rotated velocities to a numpy array\n",
    "rotated_velocities = np.array(rotated_velocities)\n",
    "\n",
    "# Extract the components of the velocities in the body-fixed frame\n",
    "parallel_components = rotated_velocities[:, 0]\n",
    "perpendicular_components = rotated_velocities[:, 1]\n",
    "\n",
    "# Add the results to the DataFrame\n",
    "dat['parallel'] = parallel_components\n",
    "dat['perpendicular'] = perpendicular_components\n",
    "\n",
    "dt = dat.time[1] - dat.time[0]\n",
    "\n",
    "f_parallel, parallel_fft = fft_of_data(np.array(dat.parallel), dt)\n",
    "f_perpendicular, perpendicular_fft = fft_of_data(np.array(dat.perpendicular), dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "02151979",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig2','t3')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(dat.time,dat.parallel, linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,2000)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,2000],\n",
    "                            yticks = [-5.0,0,5.0],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "ax.set_xlabel('2000 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig2'], 'fig2', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "770272e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig1','a3')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(f_parallel, np.abs(parallel_fft), linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,1],\n",
    "                            yticks = [0,0.1],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "# ax.set_xlabel('1000 secs')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "ax.set_ylabel('Parallel \\nAmplitude, m/s')\n",
    "ax.yaxis.set_label_coords(-0.04, 0.5)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig1'], 'fig1', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "df099562",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig2','t4')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(dat.time,dat.perpendicular, linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,2000)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,2000],\n",
    "                            yticks = [-5.0,0,5.0],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "ax.set_xlabel('2000 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig2'], 'fig2', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "91d0f9ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig1','a4')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(f_perpendicular, np.abs(perpendicular_fft), linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,1],\n",
    "                            yticks = [0,0.1],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "# ax.set_xlabel('1000 secs')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "ax.set_ylabel('Perpendicular \\nAmplitude, m/s')\n",
    "ax.yaxis.set_label_coords(-0.04, 0.5)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig1'], 'fig1', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a7b2daa",
   "metadata": {},
   "source": [
    "## Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "35ad9328",
   "metadata": {},
   "outputs": [],
   "source": [
    "dat = pd.DataFrame() # Replace this with your actual data\n",
    "dat=df3\n",
    "heading='Forest'\n",
    "# Your specific points and velocities\n",
    "points = dat[['xsrc', 'ysrc']].values\n",
    "velocities = dat[['gps_linear_x', 'gps_linear_y']].values\n",
    "\n",
    "# Calculate the rotation angle for each point\n",
    "angles = np.arctan2(points[:, 1], points[:, 0])\n",
    "\n",
    "rotated_velocities = []\n",
    "\n",
    "# Rotate each velocity to the body-fixed frame\n",
    "for angle, velocity in zip(angles, velocities):\n",
    "    rotation_matrix = calculate_rotation_matrix(angle)\n",
    "    rotated_velocity = rotate_vector(velocity, rotation_matrix)\n",
    "    rotated_velocities.append(rotated_velocity)\n",
    "\n",
    "# Convert the list of rotated velocities to a numpy array\n",
    "rotated_velocities = np.array(rotated_velocities)\n",
    "\n",
    "# Extract the components of the velocities in the body-fixed frame\n",
    "parallel_components = rotated_velocities[:, 0]\n",
    "perpendicular_components = rotated_velocities[:, 1]\n",
    "\n",
    "# Add the results to the DataFrame\n",
    "dat['parallel'] = parallel_components\n",
    "dat['perpendicular'] = perpendicular_components\n",
    "\n",
    "dt = dat.time[1] - dat.time[0]\n",
    "\n",
    "f_parallel, parallel_fft = fft_of_data(np.array(dat.parallel), dt)\n",
    "f_perpendicular, perpendicular_fft = fft_of_data(np.array(dat.perpendicular), dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "98ae175d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig2','t5')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(dat.time,dat.parallel, linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,2000)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,2000],\n",
    "                            yticks = [-5.0,0,5.0],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "# ax.set_ylabel('Parallel, m/s')\n",
    "# ax.yaxis.set_label_coords(-0.05, 0.5)\n",
    "\n",
    "ax.set_xlabel('2000 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig2'], 'fig2', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2d2da71b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig1','a5')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(f_parallel, np.abs(parallel_fft), linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,1],\n",
    "                            yticks = [0,0.1],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "# ax.set_xlabel('1000 secs')\n",
    "# ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "ax.set_ylabel('Parallel \\nAmplitude, m/s')\n",
    "ax.yaxis.set_label_coords(-0.04, 0.5)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig1'], 'fig1', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "ebf51ad9",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig2','t6')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(dat.time,dat.perpendicular, linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,2000)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,2000],\n",
    "                            yticks = [-5.0,0,5.0],\n",
    "                            linewidth=0.5)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "# ax.set_ylabel('Perpendicular, m/s')\n",
    "# ax.yaxis.set_label_coords(-0.05, 0.5)\n",
    "\n",
    "ax.set_xlabel('2000 secs')\n",
    "ax.xaxis.set_label_coords(0.5, -0.2)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig2'], 'fig2', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "f1e814f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax=layout.axes[('fig1','a6')]\n",
    "\n",
    "lw=0.5\n",
    "msize=3\n",
    "\n",
    "ax.plot(f_perpendicular, np.abs(perpendicular_fft), linewidth=lw,rasterized=True)\n",
    "\n",
    "ax.grid(False)\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_rasterization_zorder(-20)\n",
    "\n",
    "mpl_functions.adjust_spines(ax,['left','bottom'],\n",
    "                            spine_locations={'left':4,'bottom':5}, \n",
    "                            smart_bounds=True,\n",
    "                            xticks=[0,1],\n",
    "                            yticks = [0,0.1],\n",
    "                            linewidth=0.5)\n",
    "# ax.set_xticklabels([0,1])\n",
    "\n",
    "ax.set_xlabel('Frequency,Hz')\n",
    "ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "ax.set_ylabel('Perpendicular \\nAmplitude, m/s')\n",
    "ax.yaxis.set_label_coords(-0.04, 0.5)\n",
    "\n",
    "figurefirst.mpl_functions.set_fontsize(ax, FS)\n",
    "layout.append_figure_to_layer(layout.figures['fig1'], 'fig1', cleartarget=True)\n",
    "layout.write_svg('../../../../Figure/motionAnalysis.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a612e39c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "vscode38kernel",
   "language": "python",
   "name": "pywork38"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}