Auralizacja - symulacje.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Auralizacaja - Symulacje odpowiedzi impulsowych\n",
    "## 1. Wstęp\n",
    "\n",
    "Celem dzisiejszych zajęć będzie zapoznanie się z różnymi metodami generowania odpowiedzi impulsowych pomeiszczeń na potrzeby auralizacji. Podstawą zajęć jest kod przygotowany w języku programowania Python. \n",
    "\n",
    "## 2. Ray Tracing\n",
    "Podstawą dzisiejszych zajęć będzie algorytm wykorzystujący metodę promieniową do modelowania akustyki pomieszczenia. Cały algorytm podzielony jest na bloki, które wykonują różne akcje. Poszczególne bloki są opisane, aby ułatwić zrozumienie działania algorytmu.\n",
    "\n",
    "Pierwszy blok kodu odpoiwada za importy wykorzystywanywch przez program bibliotek. Najważniejszą z punktu widzenia dzisiejszych zajęć jest biblioteka numpy. Zapewnia ona obsługę macierzy w sposób zbliżony do pakietu Matlab. Dokładne informacje o odpowiednikach komend i instukcji z języka Matlab dostępne są na stronie: https://docs.scipy.org/doc/numpy-1.15.0/user/numpy-for-matlab-users.html."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from ipywidgets import interact, widgets, IntProgress\n",
    "#from IPython import display\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import cython"
   ]
  },
  {
   "attachments": {
    "tree_struct.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAASoAAAEhCAYAAADI/e3sAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAGdYAABnWARjRyu0AAFl/SURBVHhe7b0JWFvnmff9B4HQAhJIIBaJfQcb7BizGOxgh6x23rbB7Vv6JePOm2s8nbmmddtpv5m+7fW+7aT9MtNpJ3HaTho3bdo4TpPWZLGdOIsTHDs2xg7YeGWzWYRYLbELIbbveY6EkcSOgcjm/uU6wTrn6Jwj6Xnuc99HR8/Po76+/uvj4+O7QBAE4Z5c8BgYGPixp6fn/7XPIIgVp7+/X/jr6+sr/CUIR1giddyD/e/H7N8UqIjPje7ubuGvv7+/8JcgHOGBytP+b4IgCLeFAhVBEG4PBSqCINweClQEQbg9FKgIgnB7KFARBOH2iH784x/ns798IojPBYvFIvyVSCTCX2KFMVWheP9hvH/NAl18CEYMVahpNGFYEgClxBOj3ddRdvgV/PpXh3HZOAYvj36YOvrh6auCXGzfxnxw2k84lN72+XPTsICMygL9mQN45vvfwS/+dAoNtnv0BCz6Uhx49iXsK66C0T6PIO5mlrLNDzRfwAeH3sFHFQbWy6aBd/B9P0PR3/4bnj90BZ2j9vkCA2g49TJ+8Z0i/GxfMapM9tkLYagXzTcaUNtohHnUjI4rn6Lk2AlcaR9gC7tRU3IEH55shH9mPgo2qtF3/mMc//AsGrqGbc+fL077sc+bJwsIVKMwmxpRe/F9/PntIzha3sDeIvsSsxGNNQ240dwLq30ecXdy4sQJ/OxnP8PJkyftc1YnS9nmBzrKUXb5E1R2DLBeNg1CB7+Ksk+Lcfz0p6htdQgQ3U24ePkqSsuq0XqjGb1D9vmLRobwjQ/jS489hA06KXtsQZ+hDUPDSiRuzEBKQhLuuf+L+B+P3YuEQJHtKSvAAkq/YRhry1Bm7MXIuBV+w1bIw1MQoRRj2FiDsrIWDAUmIDs7BKMN5/HRX/+A53/3El4/UoqmQV8EasOgGjeg9PU/4P1jFaiuv4qDf3gJhypM8PEVwXTpXbzw6xfx1ulWQBWGsBBfsC2ju+4MDu3/DZ7700F8WtMPmToMoWoZVu4tIiY4fvw4fv7zn+Ott95CU1MTdDodIiMj7UsXz51Y+jm3eS387PN5x+68VoZ3Xvktnv/Di3ij5BJMnmqEhATDd7gTVafewh+f/w1+95f3cdHkBY+us/j4rRdx8MMaVJ8vQ7VxFEptAsKVXvbtMfoNbF/NMLGIKApWQaKOQlqUP7xYWOu89iFOVlaCxRKEhcYjKTsbWnYwA21VOPXmH/Hb53+H199j2Y9ZjgBNMNQy3nMG0HbpOA4+/2u8+FoJ6tpu4npTL8b8I9lr0aD348M4ebYVniFyNB99DS+9/C7ONNbjcrkRlrFxdDecwOXqPvjHpyJQOlsfnW0/CQicf9nYsIhAJUHCug3QtrSgvlOO8OQISHsdP7RAWDuMGPRQI2H9WoSMVuHcuSp0yWOQFD6M6vffxeGyDoiTM3CPdgQ1n32Kd0ubgMAYrI1XovfSWVy/KUJYYhIkLSfw58PHUeWZhPvy10LRcQ4X6jrgFZKISP+FFMfE7fLxxx8LQeqTTz4RfvLS3NyMtra2JQlWd1eg6sfNtl4Me2uRmhEPX9MVnLjcguHgCCgM7+Ot45dgib4POzbFwl8RAKUqEiHerajvF0GX+SUUPpCN+FA1fMUOxY4QqOph8QtBtJcIY/0SRKyNhUp0E1ffqUBLqwhSuQdEPsFIYIFK038BRw68gZJ6H6zduhUZYSOoO12Gq23eCI2Phkh/Aq//5SguWmOw+f6NCB2tx6Xy6xBp09hrCUbvJ+/hXI0FARk5SI30gVXfBIunDhu/+BgeXCdB69mPUd0uRmRGNuRtn8zYR5WdpbPsZ2GBauHf+nlK4J+wBZtYQOpoPo53Luhhti+yIUZgeALW5eQidxOraTdnQjMmQlt9z636OygsHOkbcpF/3xbkaL2hYBVoWGoOtj2Yj5xoX4wZGtHW04Gm8jK0sjPFmpz7cf/m+5GbrkJ/zzVUsBp32lqeWBYOHz6MV199FQqFAtu2bcOOHTuwlXUAmUyGQ4cOrfoy0Bk5gmOSkbElF5tyH0J+VgSULBC1GY3oMhrQ2tkPb00C1rI2fV9OGtYksZNuVCRCwwLZSX8DNqTEQCOfvl4QB8YhPYUFqsEKXGnqxnBrHcqtnhhKTkNCsI99LQva607hnKkbypz7sP2Bbdi2/QHkJYyhq5EFq5YG1q/OoZklHOsefBSPbNuCzTkZiA4Ptj/fAZEvNDFxSIwMRWhQOJIy1iA+Usl6+ATds/TR66id737mwcIDFcND7I+EvHxkqy0wHP8IFc1mh9pahPH+Tlw9+Ve88NxP8O//XYxTLGUcHx7HuH0Nha8caqUf5HJfKCQ+iNQEIDI0CL5svkIihs/oMIb7b6JN34qrn3yA/f/ft/C/dv0v/L//9Q6OX7FgwGyZvpYnloUNGzZgz549eOqpp/DLX/7y1vSTn/wETz75JJKTk+1rEhxrdwMqPngFz//i3/CffyxBad0YxsckiFz3ADYnBaDqr/+Jf/nxizh6sR0j3l7w9PSEp4cHPDxFEIn4v+0bcsVLiZi1aVAOi1B9sQoXqs7BOjKI+JQoTFaKfeisaoVXvx+S4sKhkrH+xDKxMJ0Sw+JetBjb0NnMMh4fNeKidPBjWayfUg25r8L+fGc82THx4/Pgx8ePzfHgBo0z91FTG1oXsJ+5WFSgYqEKPqok3JufiUBzPT7+sBJN9m8BRzuv4OiBt/FuxRiiNn0RX9qej7VxKttCO/ylOn4WYpEI3mzi3FomFkMuEiMmbiP+x9/8k9BRvv+jn2Pvj/bg8Uwd6IvslSMsLAxr165FQkLClCkpKQmBgYH2NYmBhjP465tH8HG7EukFO/GlB9chVceXeEEWlo77dv49vvsPX0Gu6gaOHzuCklrbyBHzwwPyqBSs1Y7C+/JhvPvRCCsxk5EU6evQn3wgDfTBoGwIN81DGBPmWTFssUI25gN/lgyIPb0wPjiEIav9yvs4W4tPC2W2PpoXDcVS7YexyEDF8PBBQHIeS21VsDZ+gnMXDcJsa28r6m7Uo83qh9BwHQKkYxizOheH80IUCO3aOIiDx9FuESEoNhVJUf5QsbRY5EkX0wl3YBTDwyx7GBiwTZZh9HbeQEPDdfR7qaGNULPOCowJzb8bN2oa0Njujag1a5AaLsaouQM3+60QeftD1MPWYBmIiW1jtmrBwycMsclB6G07g09b+6FKSUGkzLEbyxGWkIbQ4X5c/eg0LuqN6Ky5gNJLA/BUprL9RiA4MQSDw/X49Mx5XGel4IVzpbhytd7+/AUwWx/1CkXYUu2HsfhAxfCQBCApZzPyUuIQZB0U5ok1McjYEATRtQP41X/8Dh9U92NMPnmpcf5IEJb5AB7cHIeR0hfwf/b8I7771J9w7HIHWFlOEG7ANZx862l872//Fn/Lph/95jVUIwEZsSq0ffR7/Py/DqG81QJvpW1ta9t5vP/b/4NvPfm/8eKJcSTdU4C8uAAEaNOwXitDw6Ff4j//+zV81jzbFVhvhLCAEBOUjrWhsUhJCoOPY3nCTuG+sXl4dMcWxPZ/gud++C189z+L0azJwUOP5iMhwA8RWQ/hgWwdut57Fj99ah9ONIzBJzDI/vyFMFsflS3hflismf/AeaMwGw1o7huBTKWDTmG/pDZmQXdrM5pbhiELZfO1Elh5fWroQD984e8vwcjgOERyNXShIvQ1N6NvVAa1jtWtIjOMjo/FVmG5cdALflot1LJxDPYY0d7aCmPfMIuCfggMDUWwWgkJpVR3DXfiwHn8PiqDoRktRtsJmiNVh0EXooaXuRPNrV2sLGPZhd8oLKNjkChCofLsR0dzG0wDo/DyDUBwWAgCA+TwHh1EV3srmtuMGFcEQ6fTQiV1aODWPjQ3G2EW+bFlasjGB3Cz+SYGWPYUqAuE3MO1H41j2NwDY3sb2kx9GBHJoArRIiRQCZkXi2rjw+if6KNjbJnaDyM8kxMr2Pb9MMr682QfHHXukz5DLvsambmPes62H/Y65tmHaYRPwi2gET6J2aARPgmCuCOgQEUQhNtDgYogCLeHAhVBEG4PBSqCINwej3Pnzv3Yw8ODvvUjPjdGRkaEv15eDiMGEISdsbGx4x61tbV0ewJBEG4J/40h47hHT08PZVQEQbglPFAJN3wODw9TRkUQhNsiBKrR0VHKqAiCcFvoJzQEQbg99BMagiDuCBaVUfX09KCsrAxVVVU82tnnAiqVCllZWcKAagRBEEvBojMqPrB/RUUFLl265DSVlpbi2LFjuHHjhn1NgiCI22dRGVV5eTnefvttIWC5wk0ifGhauVxun2MjIiICeXl5CA8Pt89xxYiq4jPQD6qQVJiDcK4Umzf8ucWoMqqRtLMQQe3FOKO3QJVciJxw10GL+bpcGsnWLUxi/yeIOxUuBeY6OgkishfaZ+4cFp1RCYPRs2k6uPqIB7Dq6mqnqaamBkbjbE5ZK3qbW9C2CIsq4AdtVgE2FmRB68uHQ25GS3sjjNNuiO/nJm6SLJVwFwb1KH3lGXyvqAhFRd/DM6+UshO2fdmscClwJ4ydi+kzdxaLClS+vr5ChrR582Z89atfnXPi64WGhi6jt00MP10MtDF8xEH7LIK4I2BZUcV7aOyzYOMTe7DniY2w9DXivQo9KeEcWICAdJK+vj6YTCbExMRg48aNQtCabRodHRUyLS6r5Bfcp6cfhjKWdfV0o6n+U7zyX8/hpTI9hhQ6JOiU8OJnnU8vo67dA6H8sZD2nsTl2jZ4hMjQeqgE5yuMEMeHw7O9DDXGYfhps5EQKIbxWjFeeeEp/OL5YlSb+jHcoYREyrIwJ3EkQXwedOH6e83oG4nHxi9vwZoYXwx91oERkxi6DaFwvoDCwpq+FK+/tBdP730JZ5q7YTWPQyqVIzA+AYGjrI+8/gfsfZr1nWLe1j2giY7BYOXrOHuxFsPqFATy8pD3pdfP4mLtMOR+tThs315xlQkeQVFICZLZdjaBsL7rdh23ZV/G+qupV4+27g4YESrYnh2Pt5hf2rH4Qcf6qNLbtul5sggBqR1e+vn4+AgSyrkmsVhs84LZfrczCy04e60Obf5p+PKePXgiw4Jrl4tx6BorGUfNQorbaZpwCDqmvTOXjfyNOnT1JnpSnhB0PvdqB3FdfxxX7MsJ4vPFBz7qDgyJ9ejlyrn+XujFQ+hQs/m2FSYxVeFQWSMa1fn4BmvLX05ToK32Gs7aLxVb2jvhERyL/N08M3sQ4zda8c771zEk9cKlmzdRqrdderG016JFZMSQpBMlDtvbU1iALO6Dd2H67fLrvEZUHT6Dxka1bdnONMjryvDJ0Qu2yy78eN+4imuDG/EE3/7ufKgbG3HmMH/uwlh0oFoIjrcwzI4KcTHpyNucg4zsbBTk5yP5phRdZwwLfmE2TLh+vBaj9UHIyy1ANttm5uYHkBWdipku6RPEyiJHdP7DiNd04K2fFKHoJ2+hQxOPh/OjXbIpVkFcLsFAjwkbsrchj7XljOx85KXHIc7+jZA3e96a7K3YtiUb2QVcEKyBpmMAEl0itGI1vDu62FYsaKvvgshbhphkf0jqh+E/GIYUtr3s5JhJaYsD02+3F736yyjp74ZpQ7ZtWUYO8vPSESd4PIdsxysawKbtj6CAb3/LNmRvYFnbeCkMJtu258uKBqq5MyoJVP4qBPnbvH3iQB2ixEFQ91oWeeF7CAMmGVSewdDZRfdihQaBiiAsztdKEEvMYAsqj5bjcmMUcr/Gso6v5SKq8TLKj1aixemCOqsg+mVQyjTQaWz9QyTzR5DaHyr7t30imRXNH/wZzzzOAt6uf8azh0uh5/P9I5DIop5vbxPaTG24YRJhxCcCEbpYbC2IwJjhdXyj6GfYJ3wbPpWZtjtq7oMsQA5NeIjNKCOSwT9IBX/hgOzHqw5DZKSfTQPPloeEawCpN9q6FnYFbkUC1QRzByoXrCzQaLwAnXxKrT5/BjBk5ZP94eAAugfNWIifliCWC16GnRvsQXfieuRk2LKS9Ynd6Bk8h9p2185sYu2YT/aHMGOgeYBN/N/8tptPcKVWjo28DNvzJHZm2SsHkRIJSRp4KcZx9eJVdHmNQBwRAaVYBnXqPfjirq/j+0+koMd4BYfOuF7En2W7DJNlSJgmMPe3YKDfJiOeerwsdbCY4DXeB7nPwi5SrVhGNe9rVA2X2MTjOv82pATXe+shjguFXCpDgLELw7UGdLIzjUV/ASdKz9+qz6dHDU1KH+p8ylF63XauMF4vRfn102gQHhHE54u3nwrBPVJIWy0YFLKSQVhaTRhvG4GPl2NnliBIG4bW3n6UXbMFE4v+KspqrqGSZ16DXWiq7UO3OArJvAxL0GFcBiHzYRuFLDICAWMWNH/cDJGXBknxSpttXOwHXUwysnMTEN1jxuDFTvQJz7Ezy3YlrOIJM/Siv+yq7XaKQT2ultXg2kX+wH68V9tRcsyeqQnLvdHTFg+tRtj7vHG/a1QeZlS+sRdFRbvwzYONGIrfglz+poqCkbDOD+0N+/HtXUX4wf4T6JFqbtXn0yNGWNZD2BIxgspf72bbLMLe92vh45uOGPsaBPF5IvJPQMF9URhuOIDdwn1Uu3GgwQ8JD38FaU6dWQRlfAG2JnnC8OY3sYut+4P9pzGgCkWSji0WhyE1T4mulpfxbV6iPXsAZ/VDk9di2XKtbzMMcgPaWSkYxWs1UxWK9/1M6BdFf/cbVFr9kLE1Hnaxs41Ztms79hR4tryNb7I+WfTD/Tg9GIDQJC1b6ikc744H42Auf8b22nb9FqctMdiwPQthC7yNaFF3pre3t+O9995DXV2dcE/VXPDbGaKjo/Hwww8jLCzMPtcVbkk2wjw4hB4MwmTsAZSsHtdOWpkd7bRStYbV5uxMJBZDpuZ2V/6NH7+fSg0vlgsbB8fYyYKbXlk6yg3PhhaYzHyTGijYfz5iHyh09tqZID5PBBNyM5o7WJtnKDWszQsWYuGhE/xm5mbWB/iqzn1ADZ8hB3uzUskCjgJBfqy129t5a9leHGuwIjRrDwqi2Byn/Uqh4qbzaQzGTlZo1+06bkOmgqj9XTSafRC6YTdydSIXo7QSGva6uI15If2OxajFDfMyPDyMxsZG6PV64d9zwW9l4D+diYqKEm5pIAhihWFl1ycflqJTFI57789B0HKcoVmG9upLx9Dgm4Sv/W0BeCxcChYdqDh8QH4epNjz7XNmhl+b4gP3e3sv7AIaQRC3C78xuhgH/3IKJt0j+OJXtmGdzvat4W3Db/YsPojid85CuHxuVUCb/gAeLXoAOfFLV63cVqAiCOJOYNR+6aOLVV7Tl3aLht+Ezcu+VhN4YTdb+Xg7UKAiCMLt4YFqRe+jIgiCWAwUqAiCcHsoUBEE4fZQoCIIwu2hQEUQhNtDgYogCLeHAhVBEG4PBSqCINyeRd3wSQJSgiBWikXf8Ml/LU0CUoIgVopFZVTLIyCdBj5eDhc7qJJQmEyqUOJuxlmimzSTrGkVwjOqRQWq8+fPC4GqqanJPmduIiMj8YUvfAHr1q2zz5kHrWex92wrEJKJPVmh9pnTwAMaH0KVB7TscMxqD1zIusTqxWVkgDVbd6JQCCC2gFJ8sASXbWva0GViRyFbZ9Ftio/H1oxeqxiKGcaimh+20RKqTBKEZ98dAW/Rgaq2thYnT54Uhm/RavlofrNjMBiEIWG2bNmCpKQk+9x5MN9Axdd7txoDoWy9RxJnH199IesSqxbj5QrUssAx6K+BlA+he1KPcV0OHijaAH9jM6smOmAb5m4I+rKDqBoMQ2bhHhTEL9EQKotmANXv7sXZVjkSH9mDzFm6zZ3CogNVQ0ODUP7xIJWWlmafOzOVlZXCQHsZGRmIi4uzz50Ox7PVGmzNDUdzuBJqe6DiItFidpYr4VI++xlsu7YZ7/z5j3jhaC1GxQo8WPgEvsLOfEHt81n3y8jhw7iyc+Cts49wJq2CSSqDP7phNg6gi51Tz5awM6vrWdPprMuOl+23sDAJaqfXoUXm9kLsLMxBuN0WQrg/1t4+FoIAH4UfxKMmXCjehwqjCut27sY9QbZ1BExXsP/Ti+jzTcETeenOmZC9LelvNrP2cQVnyyC0hR2ZwNl3i3HE/tjWNnibOQP9oApJ7LGsoRhn9EYM9LITfVkJzjZPti+5vhSlegsk4TnICWctke/nTBMsUjkCDKX4675XUGoVQZb2JRQV2drrgEPf0WbtYNspFJ7LvZfFrK0eKWN5Y+pWdixsHy6XWZzWcWrnDMc+wPrH1iwttGFBUNmPzfm5i+sLPFAtypTc3d0tDEfMr0Xxa098QLzZpo6ODnR1dQmBbWZTsi1lvdDQBc3mr6BgK9vu9c9QdtGE4LXrka2zou4G+9QkKdh6/yZoR66hua0DXtoMpAeLYBhWIChtG776cDaiAofR1DifdTdBhxYcbzLhpjQa60NksBjKUWG4idFwHeTXPsYbH9/A6Po8bH9oM+LldTjfXItuaTxSgwZR9eZJVN/wRsIXt2NbihQ17R3QD3vDt/lDXL71OjYgPTUaIQEyeNPNIHcMIh8fYTRaITvy7EfLKRMsIzrE5MbA/1bKxNrsuc/Q2itBYkY2YpUuudRQGyrffwMnG6yI2fooCqKHca3sbXx4YwyxGx9APn9saEWHKASpEeNo+rAcjZ0iaDIT4MXa0OGDH+HGSDrytj+MVOUNXLhZj35FIhJGanChuQcjijjBBC7sp1KPnlEVUjdEQzpigLcmCjkPfRX5KSGQGD/DGxfq0aXZjK+w/qAZvITO9k54eYtworwa9d4JeOyRbchdvxYp2iAofBxfhwXNde2wemmQdd82JPtPHsda1z6QrkJ/xXGUXTVDk3wPEjxv4I3XK3BlMBo7vvogtqZrMHi5E+03vRCSGggXH/NsLN6UvBBYRLT/axZM13G8bhSNmjzkFWQjOzsTmx/MQnTqxMV3P0QnZ2BzAReJ5uG+rBQEB4jRy/KeiMho6CKiEBqfhoQYXt/Pd90I6GJnkDMmqCCeRYhqcJUvFqzHeqsHJBe70TfahW7JGOSJ69j+05EQubQDiRErCT+BHsF5cwfEG9KdpQS8zTb1odVHh0S7N3Iq4QgPz2DtIAN5W7OQ4peKCK9ErJ94POyD4Ssu5hc7KlZ9pOflCxqtgvxsqKW+0Dd3CZne9HBnQCQi40MRFaFDSlICItSDaHCS8OZhU1osC75WNPUNwHNOAak3NPFrkJ1vk54KxxGsQWevlWVLq1FAOjQAk4S9wCAdVMJ7xd70oED4aSZUoWJ491zGoT/8b5bOPo5/+L9/wLFLM73a+a87o5xROG3OLEQdNPehq/Y4Xvv1D/DkhJjxL4dwxdQFSdz/wFrlKE4+swvfe/aAkKYTdyC8rHnlTXxUJkPy9i8hh2UMjt3YqD8LsaQHSSwguCZTkyjgKwuAHztTiXykkPsEI8jx8SirOgZHMGJf2xGJyh+qIH/hJCcOCEW0bxBChkcxal8+P6wYMNWj8s0X8B/f3HWrP7x4sgENI4HzEJCyEnKkGe+/9gweLyrCrn9+FodLbRKu1SsgNQ9hmE/2hxZzNywDNlUor3X/UlGL9pjt2LNnD/7x8YdxjxC5p7KQdWeUM07X8FyEqNqkTDz29W8J+9iz5/v40c+/g78p2sBKz2isyXscT+z+R+SrG9FYdghVCzyDEJ8zPEidqkaLOBK5XyzAOpZ9q53SYiMMFTJWViUgMXIFLqCLvCBh0+Iuc8YgfVMRdn+Xt9M9+M6//gQ/3f01PJIUOreAlH9LfuIK6vw24hu8nT+5E1m3qhy2eNUJSNUapPbWQVJeiutCpzbi+ulyXC+1qUL7Oi6hv7cDEbEbWPq6DpFqT3gNTt/7F7KucMaYSc44ixA1kGVXAV1e8DF5QZvOS1U+TZR5bJvqCCSk5SE9RoIBcz2aFngGIT5P+Gd9CTdNFkRtWoeY6Up3kwHlMgnaYuIW7Ki7XSTyABi7h3GjpZMdKTvWyhMoq7QLFqbAJbxmdIQb0a2Mt7dTNk2UeXMISC1dTajr7YZ3RLJQ+sVrxyEdt2VUq1NAKg5D1oNboLNW4rm/Y6VU0V68V+cDeZpNFaqMz0Iye8NO//Dv2LIf4HeHbsAktmdJ6mhkqjvRWfwv+OW+YrSq57+ukOlMJ2cUmFmI6u2fgNx7EzDUVGyTL3LB4k/3ofiaAfozB/DM9/m8x/Htt7sxGPEANoSSgefOoQ+dVyrw8b5f4anvPimUPJOfr604MurLYR1tgDrMuRxcCTw1CVjn1476V76NXax97z/ZA4kmjrVWjgTh8RvQZ6jCv/37Mzhwph2qrPuQqzWjYq9NwltU9D0880op9GyduQSk3qGpyFV2oWX/t4X3Ye8BFhCHbBnVKhWQMpysFkpo+PUphQ/EYgWL/mw7fJlw/4oUSqUUUlbD+/mpnSSjo77cgsHq4I75rmu7uW6KnBGtOLv3CKqNvgjbuQ5ycw+Lls5CVFdxpG15KFTDnbeEp1OeQ9wB2G6+nLxXyo7DZ8lFoJOS2xmyA6E998EqkkHNpZuzPvbBEBfwTivRZdvn6xrNGOXrqtm6t6S6Uqg0KkilUviI+TJWhk60y34Rwtjxap3W5wdmt8VoxDC1TrTfmQwyExYbu8BXycKYIoj1Jd4nWZsmAekKwlLSqXJGe6DqD0Xmnh1InDseE8TqhgSkywX/6nkmOSMFKoKYFf6Fw60bnhkkIF0uZpMz8vR/MhWfKbsniFWL06UaDglICYJYpfBAtaL3UREEQSwGClQEQbg9FKgIgnB7KFARBOH2UKAiCMLtoUBFEITbQ4GKIAi3hwIVQRBuz6Ju+CQBKUEQK8Wib/gkASlBECvJojKq5RGQOls4FmKp4HaaM3oLVMk2swZB3HnYzEUkIJ3KojMqPmwLn6bDYrEIAay6utppqqmpgdE4dUTmSazobW5BWyP/IbB91jzhYwO1tDfCuNAnEsRi4KMGvPIMvicMQleEn00MwHhb+EGbVYCNBVnQ3tZIHXxEkAP48N2lOCb3YVEZ1fIISBc/rAof9O5IzQBCM/dgBzc1EMQyYmlgJ152MjZ6ySE1VeHYiX5I0wvw9SK76+5zhQSkt1geAak9UHWOYFRrweUTNqGho/DTUUDqKFF0DFQ5Y++xMtBR3OgqPZxBDgo+ts6kLNIsy0Hh9kx0lUwnGCVWM6NmM8v/uVRFBrG1Acf2nkWrOAoFezLhFBeE8ZrufAGp67hTk3p7l2WrQ0Dazz6Q4zha3grPe7jwMxXKhovoMlgQkJIOxU32gqeTKIakQmIqQ41xGH7abAT3HHcSN26Ol6PuvAG1PTLEp/rC6Cg5zXYQIsaNoMZRFnmPBn0nK9HQaJcr5q3H2hQdghR2KSWxavG0t2uhHfQ04RNjHywRsdgU6+/cNu4KASnLIA3X0T7kBc3G+7AtxR83LphQ369A4loPNK06AamAo/CzAPmbgyH1b4Ghy4DrM0kUp7G7OIobswu2IT/VCGl3KQy1LpLTLZuQFs/qX2sTuoSRvyZlkckx0ZBLGzGssiBsDTd3JCNmgWM9E3crPPvZh5/xa1TPfojesGRszwqboW3c6QJSsKAXjzXZW+2i3XxkqzXQdPSid1UKSAUchZ9iqMKi4RscAuvI4MwSxe6pH7GjuBFiFXQRXghS9cLSPwDTjUoUP/8f+CduzXj8H/Dj37+IU40N6O3nz5yURQJyROffh4jRZrz+jaW6YErcHdgufD/GPXcPJLM6rgR/PVRuU0ZN4U4XkPIS14rmD/6MZyZEu4dLwYVZq1dA6oLISwIvr4lidgaJYryj3Gc6hjFk1sCL1cdyvqnodGwq2o3vCNLQ7+Bff/JT7C56BFM3w918qbjni1/H17//BFJ6buLKW+wDmrYxEqsLnrnEIJlnCZmbkZ/UhwBPVk4JJ7tl4vMSkArZ4ye4UivHxt38+U9iZ1YqyxNtrD4B6azMIVF0oeVsAy6xSdCG6stRcqMH9T5xCI3UINXcgUhjN/wT7NtgU3KMTZk1lYkGmYeEmD6Yh5a5MRJ3ABbo9QY22XOP0W50N3rC3KGCz/S135KzkgJSDHahqbYP3eIoJPPyLkGHcRmEjGp1CkhnRYywmSSK02Q4qjgW2SuLsZett+ubxWi0JGDLpgQopWHI2paL0IEKPCNITtn0fS5pdD2LcByuQxT9HX5zwQq/DdumybyI1QXLBDpKUPy8vR3uehrHusKx7qGtrDSzr7LMrKSAlMt5U/OU6Gp5Gd/mpd+zB3BWPyRkVKtUQDrV+sLlhcbBUYikzuJQJ4kiW9dR1mi99GscqTHBV7sT62QD6HGVHrqaM2Qqm6TRH+h2lEMKx8PWs0sopazW1mlnkU0SqwZ+g3GzoRk27+ws5pVZhaOuj91VQGrrh7ckokol61EKBPkpoOCvgQSki4Nu/iSIzwkSkM4fClQEsUK43AhKAtIFwNNyp5SZIIjlwfUyymxl8G1wVwYqgiDuLnigWtH7qAiCIBYDBSqCINweClQEQbg9FKgIgnB7KFARBOH2UKAiCMLtoUBFEITbQ4GKIAi3Z1E3fJKAlCCIlWLRN3ySgJQgiJVkURnVSghIuYWj+JoR6uRC5OM4Pg/BKIlNiZWDt38SkE7HojOqlRCQ+mmzUJBVgCyt3woJRqeKG0lsuorhIwNMKxl1HFDRYZpxAMb5QgLS2ViULqu/v1+4ThUZGYlNmzZhzZo1s04BAQGQyWSIjY0Vsq3p4bosFsys7APLTkCwnwIBigBB3dNvmNRhCXqgZWEIbRffR3WrBX4x2WDxcYX2S7gjxurrMI14CYKPvHQtLFeb0NQCBCVFQaMKhDZlLdYJw/quh9ajAR7eUsRv2Ip4tfciv6ESwYe39wAFXGxVC2RqO74LaFhU6besAlK7KTlQP1l2RbW8gCM1nRgZZw3m0gkXqai9ZNR34DJacLbssk2QuDUVV0qKcbDE/ng6YSJ7aFvG4nTJX7DfRdy4yeMNHJtxv8JBE3cp1t4+QUvlo+AjcZpwgWVRFUYV1u3cjXuCbOsImK5g/6cX0eebgify0p3H3hfa2p0vIHWWiLqIeB370+oRkE5mVFzAWN7UCZEmE2GWChw/+hlaPdZPIxX1RtuHh/HXig5Ic+/HV9ercOPCSfzhvBFBa3JRmM4fTwgTA+E5rUxRhbRNMfBzETd6mSpwwmG/qcoGXDQaYAlIQbpmAfpE4o5D5OMjjEYrJDee/Wg5ZYJlRIeY3Bj438p4WJl1jrWPXgkSM7IRq3RJhe4KAakFzXXtsHppkHXfNiT7Tx7H2qBBVK1OAenMzCgVFepwFeJ0qcjPyrQLEmOQ5RmCzIyJx2L4Ntoko9PLFC3wmSJu5G5BVwnkFgT7qtDS2nUb1yKIOwt+zecIzps7IN6Q7iwlMF3H8aY+tProkDjjpYE7XUDqDU38GmTnb0MeW0c4jmANOnutLFtatQLSmZlRKip8ahKo5Eo2sTdZ7AOF2B9RTo89IBu1YnSUnSlnkCnOxNJIIIk7EuGC+pv4qEyG5O1fQg7LGBy7sVF/FmJJD5JYQHBNpia50wWkrIQcacb7rz2Dx4uKsOufn8XhUluPIQHpnDhIRRdUgfHrADPLFOfktiSQxB0FD1KnqtEijkTuFwuwLibCZWhrIwwVMlZWJSAx0pZ9Lyufl4DUVIXiE1dQ57cR3+DS3id3Iit1sseQgNSFGaWiC/kqdxaZIkFMwsq9iku4abIgatM6xEROMwa4yYBymQRtMXELdtTdLispILV0NaGutxveEclC6RevHYd03NZjSEA6DTNKRRfyemeRKfI31lnceDv3xBB3Nn3ovFKBj/f9Ck9990mh5BHulfrpPuEmZI6RnSytow1QhzmXgyvBSgpIvUNTkavsQsv+bwvvw94DLCAO2XoMCUhdBIx+I63s32YM9bKgbDS5SEVdxaWzP/YZmodM0S5uDPI2oXtGCeQKpPvE54SzfPYWSg10rF3wTGRetqNZhaOuj91VQDrqJP5Vsj4DRRD8WJ8RMjISkBIEcUdBAlKCINwOx5s9+WMSkBIE4XYIpSsr+0hAShDEaocHqhW9j4ogCGIxUKAiCMLtoUBFEITbQ4GKIAi3hwIVQRBuDwUqgiDcHgpUBEG4PRSoCIJwe5YsUHX3DuFi1U0cO9WE81c6MTA4928ACYIg5sOS3Jl+sfomXjtSg4vXOmGxjkIiFiE/S4cvPxKPSK3CvhZBEMTCWbI7089f6cD7JxpQdaMLDc29wt+yyjY0tUw3EvQSIAwR+xf8hY+pY/uR0Yxwiei+NybHESIId4e32Xc+OCCYZggbi7LQTNDVY4G+tY9FPCAlVoUtG7Xw8vKEvq0f/n4+CA2SC+sMDI4gQCGBSLREQxF7ekMWoEZIZAjUATJ4zxJuvXwUCAyKRGRgEBRjLSh9/Swu1g5DnbIgCwZB3DY8AL3ywlP4xfPFKC7TY0ihQ4JOCS/78gl6bkwamMgnKdCw6EB1/monfvvqJbx6uBpnzrfhap0J166b0GjoQ795GEEqKcTeIpwsbxGuXYVqZAgLvi0F7CQ8UPkrofSfPUhxRCxQTYhMbfqiarRa/BCbrcXd4WYk7gyMaOkYhVixDrn3ZkA7cg3NbR0QBaciQukcqkh8O4XFCUg5B4/W4rk/XWBv/oB9jjNZ6SGIiVCi7EIbLNYRfGvXOux8KN6+9DYRxsGpgkkajuzCIHRycaNxAF0w4GzJWRgcRIr8LCaITFmDaDn9V+xj5eKQSIa0x4pQZJcwLhRBqnjNZBvE/mqJIHWckJwGtU+KU23bdhZLkrh0tWKF2cz9MWLIWCpvqjyId6otCEjfiR2JctsqdlrL9s4i3LW3vxmknry9c5W7mrU/QbjL5Qz8socqCdt9a9i/9ehoYc+7YoY8uxDbM7tw/G4VkHIuXOvE6YpWmC0j0LJM6ZH8KHyhIAb3Zupwb5YOG1mgau0cQPmlDnh6eCBjbTDWJs6kc18gTpmRLzq5gPTjGxhdn4ftD6VCef0CTA39UCSthXerPY0O3YyMKCmGm70RFJWDR4rykRLsD9lcKdk0DLVV4v23/oxLNwOQfd9XsHX9pPI7ONwTn127CcOYBunR/vAy1eFwXTvaVdHYlKCeMwMk7lZEdiEvbwBWtF+7ji6LHBGpKQiROTcKnlE5CnedxLdeTbNKPa2sbKzt6IU4xK5z776BD2s70OsdgvWi66zdnkSDJQZbHy3AhuA+nLzYgKa7WUDKx/Xkk5+vGFsytcK3e9X1XbjKyr9L1Tfx5gd1+Oi0ngUyfpvCxNrLhQpxMenI25yDjOwC5GerIZbr0eToDuO2jchIxIVGIUKXgqR4V/3RAtGlIjrzAWzO5GLF7cjL9MTw6Dk0iXRIGlAhqNYk2Dy4BEAmMSIhKWpJBxMj7jx4dnHg2e+hqOgHePWiNyI35iJhUrvsxPTi2zbU3q7UUxeOcLbNnIxkxETKIbnbBaTeXp7CFOgvwdj4OM5caMWHnzbhw1NNQoA6w0q+NpZR8Qvt3t78bLKcvVQClT8LDv5cuCCGKiwavsEhsI4soyaUD3AfqIEf/1xZKRkSroMiUAnzeDDiEjshU1ai1WR3v7XFIm6lnUqE2+GtTsDm7U9gz56vIiu0AUePFeNozfTfRk8vHx1C3+1KPRW+kAX42U6avtF3v4A0MSYAKXFq+Ph4YWhoFDdNg+gbGIaEPdaG+CImXImkWJUw8Xuq1sQ7173LichLAi+vlb0YxMWKMk8zAuQKhKWlQRQQhLbKi5+b+41wP0QyNSLi0wW9e86me6DzlGKwuWtuLZuTfHTppJ48cNz1AtKE6AD8y99n4OnvbcLGtGChsPNkW4sI9UX2uhDkrA/B//pyCn7+L3n4+6K1woX1Oxr+gb37oeD7Ez6Wq7XQl9fY7uNiyz48bsVVu1hRHJgI3ZgvTOXlQIAIaWlhSzbQPXGHYtJDb2CTcN+fCIPd9TBY9eiT+bDiar5IoJlD6in31aJhQIIrHXypEVXHT6KKf9szE3e7gFTKMqdYFnz4BXJ+K4LI0wNjY+Ooa+wWyr/T5a0YHR1nWZcKOpZh+YiXs/SbJ9JwJGzoQ8u1f8PTzy7whrqhXjS3tqHRZGZJLUOngLnzGPb/kMsbn0HFQChyt9nFiiIlC9iDaBkfQ9eYDon0FTMhMqP23f14mos6WQm1++VaeMXej0fXBi/ADymaU+opj05DUlcbzv6Qy0b34p0qD0gcMiAn+Mn3bhaQusKvRR06dgMH368TSkCNWiZcYP/yw/GIj/K3r7WEzCJu5LUylx4aB0chkqrt8tJJkSMXRzYbmtHnFQaddhZ5pCtctGg0Y4ydgQJufoxfnW2GR+A6FLLPRJBAutg3BqqP4GBlFySJ27EzXbWAxkjclbgaWxxkpq5MkZs6yUdlrL6aTerpKE9lyzQKQOEDsViBUK8+GPqs8GT9Yoo89G4UkLoyyjIpfgd6W6cZ1uFRIXtSB0iFC+38TvX5UFVcjOKDB3HZ/piTmMIjM1B9daLmZegysaNwJwqzw1lyeRs43v9hn8Xv80hJ4H8NuFojzLDhcn8JD0J7z7ZCnrADe7JChXnOGHFl/ynUmv2R8kQuEuYbDAniToZlaG4nIF1qbGeBZieFtp/Cdvmwr1c4B9lgUdumqb5NpfoUJxlHCoVwu/ogeh0LdZez32yBit9wV1z8Kq4PP4BHC7+GzLSlEzEShFvherInAal7IaSzjim0AxOl5ZAPT3F1ttsXCOJuhASkBEEQNnigWvS3fgRBECsFBSqCINweClQEQbg9FKgIgnB7KFARBOH2UKAiCMLtoUBFEITbQ4GKIAi3Z8kCFQlICYJYLpbkznQSkBIEsVws2Z3pKy4gXQR8vOqS0yVzjkG1ePmjBfozJThVMrcUlSCWhtXT5m4rUPGhXWrquwQLzbf+Zh1+9I+Z2JYTLgztMsQyqxtNPSgp1eNyjRHDw2P2Z30+8B8Sdxo7YRSURTPDf1Tc0t5oX4+rrvbhZ3zQLzZ9b9bB9kZhNrHtd/KxseyzCGIeTEofeDubaezy6Vg9bW7RgYoLSP/9t5/he/9+Env/eAEH36tD8ft1QlDiI3ty6hp78PJbVYKotLKqU5h3R2HqwqguERt378Ge3Y8h3tSKC0fLKWMilhAjrnd4IDBmO/bs2YNvFIaiXX9q6tjlq5xFe/1OfdaCYhacbuh7cbNr8NY0cRGdyx0C/CW4VmeCoaP/lgxiZnj2UoLzH/AP6V08t/d5VJs8oFF44pPiF/DUL+yPo1MQaPc2OCmyi6th8tAgakLVzsfKef0P2Pv0c3ipTI/uYTPGpTLI1QmCfZafxV5/aS+e3vsSiqtM8AiKQkqQzNlSG6yAr0aLiJhoREeGwsvYAjOLweqkBPDRhR23caa5G1a2UCqVIzDettyG/XV9egHl5R/gpV/YjkfQeat7cc5FMS9ss+wiaq1q4XiIux0vSH1VCI+KRmxsJMKDPdB+tg29vQrEbQiFs5p0jjY36tDmi4vt/SUGg5Wv4+zFWgyr7X1H6Bu2dif3q8XhafqBE4596dZ2Hbc12c9MvXq0dXewVh+KcKWXcz/jIl6LH3Tx4VAuzO2weKX70gtIe3CDi0QrOiDNvR9fXa/CjQsn8YfzRgStyUVhOn9sQn2/AolrA+HJ3oDiTy6jVZGNwh0PIDdBipr2DuiH/ZCsG0GdixjRfKkMlxs8EJm2wSZG/LQa9XbxYrKkBm2degz7JcO/r9whUMlt0kiedw614+r1Lgz5RmBNSgik3VVO20hXmXG57DLqxyORvsExUNlf1y1B6mbEy+twvrkW3fJYhJhOor7rJuTR6xEstaC5ogKGm6MIT1mLcF8aGfTuRwQfHx822T5rz4EWnDZaMKKNQW6sv/PgkKbZ25xv53W0D3lBs/E+bEvxt/cXJeLCbuJkvQk3pdFYHyKDxVCOCsNNjMpFqG3QzyEgZcHRMN12eT/0QNMKCUgXHaj4N338gjm/HnV/XgQiwhSobexGp2mQHYiZlYYd+IwFqe5eC2QSL2StC5kjUPXDUFYPsygSuYUPICdVDevFVowPhWLzYxOPGzDQ5YHQDf7ofPcMmtrU2Pil7chNiUZEJAteLXXo7xmC1OsmSuoMGEzMw5e2pSM6LARqTxOMHlL46WIhrT+C6pZBJOZ+CfnpLFuSDuJcqxVGbx0Sxy45eP/HoD/zOv6w92k8d6wDPnEZKNiUiCD5KAzn3nTaRliwGp4mIyD1Y0HGMVBNvq7srzyILHas0UFe6DhtQm9PEJLWDqNjCBhVJiBW2oaKix2wsga1MV0L6ZLdPELcEQgjZp7CDZMSmQXrEOV0orJAP0ebC/L1hUobjZi4WERH+7P+0o1xixzxecHobmSF5Lg/EmJlaDt/kbU5CaLTQtF7ugWjoljkf2EjEoMCpgQpjod4uu1KEaxtwbGqOfpZpwXJ938FBRtiERmhg3SwFB2WNojD+InZvoO5WbwpeXkEpBKo5Eo2sV4u9oFC7I8op8cekI1aMcqmAZMMKnEYIkPtQxKLVdAFdUHuWYebHXOJEbtQffw1/Opfn7wlVfzLkSvoMvW4XBfwhjp+M7Y/sQd7/mcWgm8cxQfFR1FrsssVZRroNLb9i2T+CFL7Q9jFFBwFqexQA3WIEgdB3TsMiS4RsjFf9DW1obv9Brq8RiCOiIByPm8XcdfAS6QDhz7COd9kPPqlHKRPnunszN3mRDIrmj/4M555nBtf/hnPHi4FF1uJ/COQyGpI394mtJnaWCAUYcSHJRcsmMwtIJ15uyQgnRcDGLLyyf4Qwxgy+7H3RwUf9q7MLkbUIinjMez6JxaAuFTx+z/Cz7/zNyi6J8zFsSaCTB2BhDSbrjpvgw6eskE0dfHtusoVzRhoHmCT/eFsWIcwoPECdHIo/BOQxP7tN3oVF652QTQsZg1IaQu+xKqAB6lTNS3wisjFF+9bh5jImYbyna3N8Wuhn+BKrdz25c+eJ7EzKxWC2EqkREKSBl6KcVy9eHXyZCieh4B0tu0ySEA6K2poUvpQN3QSJRX2N9Z0HaXtHqj3iUN01FxixAD0ePugW6zFOu7E51NaAiJchBFGvR4GNgnbFw2iu96AYT1/k/1scsXefpRdsy236K+irOYaKu37K33lQ7xz6wzVgrMNl9jEH7EUvqIE13vrIY4LhZydZaIiAjBqacaHehG8NElI8KcwtWpgbaX80k0YLVHIS4+Z0gYnsQs9Z2xzXWiq7WNtOgrJW1h7TtBhXAYh8xFOuJERCBizoPnjZoi8NEiKt58M5xCQzrZdEpDOiRhhWV/EQ5kyNB78Jnbx+09+UowW30Q8vCkBAeq5xIi52BI/NPlcNv1sXzGqTLatT+BtrsXR/U/b19mNl+vYa77/UaRpvAW54tYkTxjetG3jB/tPY0AViiQdeyIf9L6xDS3NvbCdTFSI8zCj8o29bDu78M2DjRiK34Jce2MRh2mhMBigbGlHQETUkg6MT7g5/Z24cv5j7HvuKXz3ycft91JN1x5t8tEZ25w4DKl5SnS1vIxv8xLt2QM4qx+6lfnw5VrfZhjkBrSzUjCKN7J5CEhn2+4qFZBycaKjSHSux5PGF8Gf6KrRchEjalgwlUp9IJaqnUSkwnMZSo3NGuPDMpsJ+aM/uh3kiZPrTJhlHLchVWugkrJ9iMWQ+Usx2DohSO3Hhb1HUG30RdjOdZCb2cpT5JOtOLv3GBqsocjcs3Q+NOIOwEkCOolrW5tgxjanVsNnyEH2qVSygKNAkJ8CCrvws7VsL441WBGatQcFvJHNS0DKz7uzbHe1CUjvXngQYoGqnwehHUj0tc92wKL/BB+UdsIz/F48kBO0ZD40grgFK7s++bAUnaJw3Ht/DoKWo5Eto4B0SaIIL/sCWWBak6DGPakapMarERIooyA1F8K1rGfwg/84huvyaKxPV1GQIpYY/nvAA3jmh0/jWJ0M0enpUC1VI7O33+/Zy9Wiv3sGFQOhyN228NJuLsjrt+xMLVlvMSFwZFlzwFLYnwliCqMwGw2sXOxildf0pd2iIQEpQRCEjSUr/QiCIJYTClQEQbg9FKgIgnB7KFARBOH2UKAiCMLtoUBFEITbQ4GKIAi3hwIVQRBuz5IFKhKQEgSxXJCAlCAIt2bJ7ky/EwSkBHH3sXoEpIuWO3D40C761j4W8YCUWBW2bNQKIybo2/rh7+eD0CC5sM7A4AgCFBKIRB72ZxIEcQs+CsEtHdVClFJDaLt4AfruESjjHIUidx0Niy79uID0tcPVuHbDhHEHCbKx2yKYaJJjA5CeFITm9n7Ipd74emGyoMwiCMIRdrI/tR+nLt6ER+RWRKIRJY3jCEzLxRO54ZDY15qeAVS/+w6q+2VIfGT6sc7uBnjpt+hAdfBoLZ770wW0dAzY5ziTlR4ijJNedqENFusIvrVrHXY+FG9fOh18EHl2NtF34DJaYEAUCgvzEXD2HRQfLGHzAG3WDjavEDnyBhTzQehVSSjMtn2YXEZ6Rm+BKpktD7d9vLZ5Rgz0AoayEpxt1iJzeyF2FuZA1uC6vn3/gyokseXh81f5EMRt0ImKfXx01yDcs7uAtfoGHNtXgU5xFAp234Mg+1oTcBFEcXExjpQZWH/YgjVSCVQBQbZAJeK6rYMofucs6z/Amq07UbhzO+Q1fFhjCcKzC5GkYgt4BldcBZM0HOsyu3Dcvj2kbmV9oxCFyS4iFmF91+06bsu+TJeJrVlaaMOCoArPEfqV4/FyqcpE/1tI/+KBym0FpDs3pSHa+ybaBj0hSdmKgmwthq81o7VDhJC1KvR/cgE39WPQpEfD39uIuiN1aG9XIXpTPNSCMZRt8caHOHzwI9wYSUfe9oexOV6OuvMG1PbIkBxlxWl2FjOMaZAe7Q8vUx0O17WjXRWNTQlqm3SUIJadUZgaTqHFOsTKt/UItjbhk+pWmNRaZKRo4GNfS4AEpO4nIM2IULHn+UIREomY+FjERjso1VPXQTNYhm4W3HzCNyB0qA7v6vvRr0vFVhbmJ95mrmevN3sjMvureDArBdHRwfC6eRqmXhNUKRsgrexm0cwHug2hsNS8C31/P8LXbEWiyvmDIojlwwu+wWqMGctx+I+/wW8/McJvTTYevTcBwTJvh/uHWIlIAtKFs+wCUo7YG92XD+HFHxah6PF/wI9/fwyXBO+OGGHxSehk61a2GGHUl7Og1obY+LApQ/lKVP5QBfnbRhzkktIILwSpetnHHoa4xE7IlJVoZR+2oUIGSVssO/vcvVckCTdksAWVR8txuTEKuV/bgz1fy0VU42WUH60Ei0kOkIB0USy/gJR/9foXVLByLPYRLj78R/w/j9wDtksBcVg60rzUCOq4iIvzDjJcUqph5zAt5DIW7NLSIGL1fVvlRZTLJGiLiVvysZ4JYjYs7bU4N9iD7sT1yMmwiW7XJ7JaYfAcattdOzMJSBfM8gtI+9B5sR+9rRGIvYd9gOmRULMQMzgR8ll2lKgbhi9LmcvH1RClpQtBhl9Af+eDAyjV297ulrMNuMQm/jQLy7xKbvQIklJ2mBAHJkI35gtTeTkQIEJa2tSMjCCWE28/FYJ7pJC2WjAoZCWDsLSaMN42Ah8vx85MAtJFsfwCUiXis5PR13MaP/w7lnL+8Hc40miCz63ETASlLhSDzWMYMw1Dl2AzuHDvWUt7I4zmUWEtVRyL8pXF2FtUhF3fLEajJQFbNiVAKTQKJQusg2gZH0PXmA6Jd/GNKIR7YpN4RmG44YBN0lm0Gwca/JDw8FeQ5tSZSUB62z+hWR4BKZ/lKDdUQimVQqX0g5/avry/GkeKL6BLmojthevAr4HzQDUhELVe+jWO1Jjgq92JdbIB9EwjQByoPoKDlV2QJG7HzvTJC/EEsWK4SEhnko9ySEB6G3xeAlLj5f04VWuG/5onkBs/VTXFP5gjNQMIzdyDHfxq4hSMuLL/FGrN/kh5IhcJrp8QQdwtsLKLBKQrLSC1p6y7X65Fa2A21kcu3IfHr2Xt++luvFzXisCc9bZUmCDuOviXUiQg/XyYSDeHfKDTslRSMf274lgGql0C0UQKPeQzc5pNEHc+JCAlCIJYdpas9CMIglhOKFARBOH2UKAiCMLtoUBFEITbQ4GKIAi3hwIVQRBuDwUqgiDcHgpUBEG4PUsWqEhAShDEcrEkd6aTgJQgiOViye5MvzMEpHykwn14ax83cthnrTDcyFFyuuTWoH4EcXvwHxyTgHRO7iwBqRd8FIEIjIxEUKAC04xhvyj4KAyvvPAUfvF8MYrL9BhS6JCgU7K9TWWorRIXmnswoohDAg3SR9jhJ7DXX9qLp/e+hOLiapg8NIhKmY+lZfUISBcdqLiA9LevXsKrh6tx5nwbrtaZcO26CY2GPvSbh4VRP8XeIpwsbxGuXYVqZAgL/jwNiSIWqAKgCFi6IMWztJaOUYgV65B7bwa0I9fQ3NYBUXAqIpRTQ9WwsQa1xmF4qxMoUBF2LKzN9MFTpsOmzXlI8L+BCzfr0a9IxNqguULVMIy1tTBavREYT6bkaVkZAWkh1uP8pMDQRZDoJDfUZWJH4U5BSDrAspzi4oMoueIgLQ0fmBSMbg3A2cqpAlNeEqqTCxE14LBdh33ahKZN6LgKtLBFUYWP4uE0LXu+GDLWpkyVB/FOtQUB6TttA/W5yBm3rJVAEqBGUMJ2p4H8bNslUerqZBRmM7cfiFgbEsPacAx7y1ohjirAnqxQ2yoOkIB0gayMgLQN75XV3BIuJktq0Napx7BfMoJ6z+EvJSztlaXblgWrEKaWQjlci0MX6tGl2Yyv3L8JmsFL6GzvhFeIDr2flqOxUwRNThiGTkwVmLa0BkAXZ8XHpZOSR8d9+t0sweGD59EhzcX9OzchLUqLIKUMPoKt1Ir2a9fRZZEjIjUFIbIuVLnIGc2XynC5wQORaRucMioSpa5mPOHt7X1LJ9fd+AmM/RZExm5CjL9L6k8CUncUkMpgvPC2k3AxWjqIc61WGL0DEdxahZY+KdK2PojcFLZMF4Jg9Qgaj5ZD3xyAjB0PIScpCoFow/XuYViVYZDU3GBpsh+0udMLTHvCYhDVU4K6afepQ+LYJdSbvRCZuxMPZERAJfOG9db1hQ/R4R2LjPtykaiRY9RwDm/OKGdMcQpUJEpd5fARa/e/gKd+8TxqLUFYn/sFbAgPcLlEQQLSRbH8AlIuMOxC9fHX8Kt/fRKPFxVh1z8/i78cuYIuUzNaakWQ94chMtQ2DLGY1V5SmRcGTfWofPMF/Mc3d6Go6HH8w//9A1482YCGbsdvIB0FpgboLzdALuK6rGCMzrjPHkFRJFHJoGTTRJjh15s2b38Ce/Z8FVmhDTh6rBhHa4xzyBmnQqLUVYyvFlkFj7E2tAcPJAFlJX/FoQpXvx4JSBfF8gtIOVokZTyGXf/ExYds+v6P8PPv/A2K7tHAc7AHAz0DTnJDGzFI31SE3d+1Pec7//oT/HT31/BIvLNXUBwWh9g2CWSnynGmcwyePhO6rJn2GYbplIkimRoR8enIzs5AzqZ7oPOUYrC5C1zJOLuccTZIlLqqmHDrZWcjc3M+kvsD4Onq1xMgAemCWX4BKRcYBqDH2wfdYi3WsQ+Rf5DZaQmIUIcjJKUPdUMnUWI/8xj1ehj0A1CkmNERbkS3Mt62Pp+SY6aOqy5mmUpMGyRDH7J029f2wYlm2+c0AgmTHnoDm4R7WEQY7K6HwapHn8wHfrPKGW0X0EmUSvB2IbQh+719Yz3daBSxNqz2gY9tlp3VLSBd9DUqXvaplBIEB8rQYTTj3MV24S//6UxzWz/aO81ITw4SsimFrxheorliIr9GVWO7hpTNv2r1hFgRCMlABUrffh6//9PbeKO4GNUmD2ii1yIiIgx+Y9dQdvh3wrK363ohiVyLnLQIiDrP45NXf4eXDryFYv4NmYV9GPH+6C933D778ETd+PRMExqsEci/NxVhLH+deZ8pkJjKUGMchp8223aNydqG8rdew75f/jdeeaMYbzWwQL3xYRTeE4sARQBC5RYWfI7g+d//EW9/1garVIIAVqerdSkI7jmO8qZOiDSZCLNUoN7cie4aA8pe+i1ePKSHKDwbD29fj0i5NzsBiCEeqMbJ6+2waNajYK0W0kWfYgi3wrMX1w8XY//P9+LFt4rxZkU3lGmb8YWtyQiWsc/evhrPKcSKUMgGL7E2b2ubZ9uskEgDWAmnRnhqGnQSA86dK8aLf3oVb9U2wdjsizUhWsRka+HHyq6x9hM40dYCz4R8bGWVjsjh2ljx25cwHJSCrTsyEMWC5K0w4imBxHP67cZtSUb8jG08FWuiIxAk6cKFE6/hdy+x575xHdaQDOR/aSOiAhz2MTeLvz3BkWUTkPK5DsJFjqOc0WmZktXugpFGZLdutMBk5s+YMGP4YbR1OsEpmydigWwe+/SxuFhtXC0ct47Bnu+4yBk1rHaXSn0glqrhN9JKolSCwds9ayPNHexzZ7B2EsbakHa6DJ7h2DZJQLpAPi8B6d0CiVKJZYeVXSQgXWkB6SqCX8siUSqxeEhASiwRPJ0nUSqxPJCAlCAIYtlZstKPIAhiOaFARRCE20OBiiAIt4cCFUEQbg8FKoIg3B4KVARBuD0UqAiCcHsoUBEE4fYsWaAiASlBEMvFktyZTgJSgiCWiyW7M305BaTcYnHg4F9wYMrIg7eH68B182W5jmd2Vo9oklh827ybcXsBqYe3DAHqEEQGq+HvNJDY7cHNLxMD183m2OONpuT8eRjFCQhXei3h8fAhXl/BC0/9As8XF+NM8xD8tLZ9TGXViCbvanhbmo+sdr5tcxWxeLkDF5D++28/w/f+/ST2/vECDr5Xh+L363CZiw1G+cDEQF1jD15+q0oQlVZWdQrzFoowJrkufPqhgFcAPnpBS3sjjOZR4fGSHY+pC6O6RNs41LsfQ7ypFRe4QYcyprsUI7rGwpGY+Q3s2fMNPBbfjdbK91BOWdO8WHRGdeqzFhSz4HRD34ubXYO3pomL6FzuEOAvwbU6Ewwd/bdkELPBy6pbausqEzyCohBjqcTrZRdRa1UjxdeI0tfPovTjE/jgzAE8t/cDYZjhQEkz3tvPn2d7rIsPh3JEj9JPL6Ou3QOhwlmLl08ncbm2DR4h4fBsdxhWWN7OtvsH7H36Obx0a+jhGAxWvo5X/rgfb3xUjfKTHRj0DkaIXy0OTxxPkMz5mG8Ne8z27z2RjX2KC+fL8cFLv8BzLznoukUS+Gq0iIiJRnRkKLyMLTCz+K5OsmVMjts909wNK1solcpdjLg8KyvB+U8voLz8A7z0C3b8t87UPmg98zrOXqzFsDoFgVwMwscPYu/fxdphqOelDCeWDi9IfIOgjWCfd3QEQr060TLAP/DkKVkT16fVGHvQ3VyPT1/5L9ZunNvVXG3ufM01mH3tnzkfbriiGtcG5Pa+VIqPTn6AMwc+xqU+ja09u/Q53q4dcd6fi3JeaFP2vsPanqlXj7buDtYyQ4XqYLZjXQCL9/otuYDURa44IUP06bqM92taYfGNRXZgPyrffwMnG6yI2fooCqKHca3sbXx4YwyxGx9APn9saEWHKASpGjOuVOrRM6pEXHwgxC7lk0/nZKCKsDRNK25MzYyHatzAOn4Q0u77Kh7OjoLCfA3HJo5HZphVsGh1EYumKh103SFKm3iS57RD7bh6vQtDvhFYkxICaffsosnJdm2Xtn58A6Pr87D9oc2Il9fhfHMtuqWJSPJsxAfTiSfDU7E2XE7DGa8oXBnHRaP8A3eV1ToXNjxQHT/6GVo91tvbTQMuGg2wBKQg3atpzjZX29ELcUg2tH5sY9038GFtB3q9Q7BedB3vv3USDZYYbH20ABuC+3DyYsMcAlILmuvaYfXSIOu+bUh2Us4PTpHsLpeAdNGlHy/u+OTnKxbGR+ff7lXXd7EOZ8Kl6pt484M6we9ntvAMa2LtWRCNwrN+GP6DYUiZyRwjEI7w8AxwPVXe1iyk+KUiwisR6yceD/tg+Mp0qqGZ8dbEY032Vmzjlo2CfGSrNdB0DMBLHYHIKB0iokIRtzYBMU5jPbMM7XIJBkQD2LT9ERTwY96yDdkb2FlpvBQGu1VEFReH9Lx85GRkoyA/G2qpL/TNXezZ9pEXv1+Eoh++isviSGRtSoBSZN9ujwkbsrchj203IzsfeelxmD4hVSEuhgXCzTlsPb6PfCTflKLrTKtQWmrFanh32PbXxj4fkbcMMQmqBY1ZTSwNPLs48Oz3UFT0A7x60RuRG3OR4GpDtuPcbrYg2FeFltY21M6jzc2KjvUfts2cjGTERMohmbPPeUMTvwbZ+ba2KLThYA06e63s9VxGSX83TBuybX0nIwf5eemIi+Ou96F59Y/5suhAteQCUt/oeckQAQV8ZQHwk4kg8pFC7hOMIMfHo+ysNTiCEfva82EmweLszE+w6CgWFQeEIpql/yHDo+zZ3lDHb8b2J/Zgz//MQvCNo/ig+ChqTXOLJp2RQOWvQpC/bV1xoA5R4iCoey0YmUk8OUPnIJaXmWS10zF9uxlC3+1KPRW+kAX42UbgnFefE0E20oz3X3vmlpD3cKmtd6xOASl7EdPJECfVhsuFaVbB4uwsULAo8oKETbZ4wxoAy9gS0mxnorwNOnjKBtHUJahLXbbrKJqcA+sQBjRe7Mwph3wm8STFqc+FmWS1c3ZZp3azwDY3GzP0Oafj4de4TlxBnd9GfIMLeZ/ciazUyd6xOgWk08gQe+2LFoxUhgBjF4ZrDegc5Gn3BZwoPY+zrp19VnHjbNyeYNEkCFPtjUI0iO56A4b1/AP0m100yebwkvHDd4tRJaTPLTjbcIlN/AjYsooSXO+thzguFPwq1IziSWJlmUVWO/8uK4FmjjYn99WiYUCCKx18qRFVx0+iquQKX3N65hCQWrqaUNfbDe+IZKH0i9eOQzpu6x0rKSBddKCSsswplgUffoE8iKV63EQzNjaOusZuofw7Xd4q3KaQwupVHcuwuEJrVnjk3vczu83iN6i0+iFjazwW6le+hSgYCev80N6wH9/eVYQf7D+BHqlm6nUecRhS85ToankZ3+al37MHcFY/dCujUkdlQt3eieJ/+aWQGk+W1iIo4wuw48E4mMufwW5+3Lt+i9OWGGzYPreFw8tci6P7n8YuweDBLTPs/bz/UaRpvIXtbk3yhOHNbwrLf7D/NAZUoUjS8WeylNrUiLbWZvQKJzIV4jzMqHxjL9vOLnzzYCOG4rcg95YJNwxa32YY5Aa0s1KQLDafEyIzat/dj6dZW+RtfPfLtfCKvR+Prg1ewIlj7jYnj05DUlcbzv5wN9vPXrxT5QGJQwbkxDz6nHdoKnKVXWjZ/22h9Nt74CwMQ7btifwTUHBfCjxb3sY3+ev64X6cHgxAaJKWLfW8rf7hivsISGeQIQpSxT4rPKVqVs2MwNjcBytPWfmFbcGCMfNjkYP80FXW6MXqqAnziz+6ZxYs3jJtjMJPp0OoyorOieNRiGcVLE6xy/BtGc0YZccXIB1E663nsWc6iFU5M4sm2fvZbUCf1ZP9W4Rrzx9BtdEXYTvXQW5mK7tKUBlTxJPEyuNqbJnmc5pgtnajZn0Ls0o9HYWmbJlGASh8IBYrEOrVN9mX+H7nJSCdsNjYhL5K1j+gCIIf6x9TtkECUmJ6WnF2LwtU/aHI3LMDib722Y6wdHvZxZMEwWEZGglIiQViv/1hOcSTBMFhJ0ESkBLzgKf5RphHxawsnSllXwbxJEFwXMtZEpASBLFaWbLSjyAIYjmhQEUQhNtDgYogCLeHAhVBEG4PBSqCINweClQEQbg9FKgIgnB7KFARBOH2LFmgIgEpQRDLBQlICYJwa5bszvTlFJDeafBxsUtOl5A8klgB+A/PV4eY1u0FpO7AfMWRnKG2Slxo7sGIIo7kkcS8mFVHNSurRky7eF0WF5Byseirh6tx5nwbrtaZcO26CY2GPvSbh4VRP8XeIpwsbxGuXYVqZAgLnm7AJHfHiJaOUYgV65B7bwa0I9fQ3NYBUXAqIqaxGg8ba1BrHBYG8qdARcyNhbWnPnjKdNi0OQ8JTjqquULVMIy1tTBavV18j3cdDYu+RnXwaC2e+9MF1okH7HOcyUoPEcZJL7vQBot1BN/atQ47H4q3L50efmYpLi7GkTIDkLoVOwsLUZisdp4PLTK3F7JlOQiX2rIdPna4OrkQSdzSw4dXvWYEVEnY7lvD/q1HRwt73hUz5Nlse4VRMJ8pxsHiI8L46dqsHWxeIXLCB1DF9lF8sASXnfZhhVmwJIshY+3GVHkQ71RbEJC+Ezu44oWPyVN8EMXvnIVBl4ktayWQBKgRlLDdttwOP84zeiMGegFDWQnb9+Q+ZA18mQUq9hpywiV8bXYsZ1g6r0KS/XUSdyujrH1x+4GItS8xrA3HsLesFeKoAuzJCrWt4oBjX9BmbcEaqQSqgCAkPrIDiSKHtsjWXbN1Jwp3boe8hvcRCcJZ+xf6iNBmq2CShmNdZheOT9PnnHBs4+yhbbuO25ps/1uztNCGBUEVniO05dn67nzh16jcXkCqGKhdOuFiciAGr76NTyr1kKc9hu1bk6FRhUEtVWK47hAuNHRBs/krKMh23EcwlDOKI7umCBjNl8pwucEDkWkbnDKqHhcZ6eZ4OerOG1DbI0NylBWnL96EYUyD9Gh/eJnqcLiuHe2qaGxKUNskpcRdiqddTGobwKm78RMY+y2IjN2EGFetmUsfcRXT+nZen1akGxd2EyenE9HKRaht0M8hIGXB0TDddhVIXOuBJhKQchni2O1LDB2Fi+oBNBm8IVJvwX1b85CdXYBt2euQIDfgeN0oGjV5yCvg+9iENJb8eVib0DVoO4tNJ46cXcA4FUepZHbBNuSnGiHtLkWnTxySBlQIqjUJBhCjvhwyiREJSVE02N1qwEGysPdYL7TJ25E1ZYhMy5xi2plEupKZRLTJ/vMQkM603V70koB0Qoa4BBJDR+Fivwl1nnL0h0aybIjP4OWclC0bhOlGJYqf/w/8E7dpPP4P+PHvX8Spxgb09rPjn0EcObuAcSqOUkmIVdBFeCFI1cuaThjiEjshU1ai1WSEoUIGSVssOxPevRcdCAd8tcgqeIy1rz14IAkoK/krDlW4+PUm+sIsYtqZRLqimUS0uth5SX9n2i4JSJ1kiEsoXBwZwWDfAHrNQ5hyO2p0OjYV7cZ3uGRxz3fwrz/5KXYXPYJ4JT+06cWRgip0VgHjbAxjyKwBe7cgl4kRlpYGUUAQ2iovolwmQVtM3JKPO024KRNuPZZ1ZG7OR3J/ADxd/Ho2XPuCo5iWX9ecQaQ7k4hWPFOfc2SW7TJIQCrIEPsgWUrhIktZU3vrYDlRgnLhPicj9HoD9GYFUs0diDR2wz+Bpag8TWVTcowOfv0ziyP9ZhUwsq1fK8Y7Hxy4dU9Vy9kGXGITP1ILK+9KbvSgnpV97C2DODARujFfmMrLgQAR0tLCFqQUIu5QWJsR2pe9FBrr6UajyIwOtQ98bLPs2IWeM4lpZxXpziKinUNAOtt2V1JAuuiL6bzsUyklCA6UocNoxrmL7cJf/tOZ5rZ+tHeakZ4cJGRTCl8xvERzxERep+9/AU/94nkUv30Jw0Ep2LojEwlRrEyTdOHCidfwu5dexVtvXIc1JAP5X9qIqAD2Yfp6obf0DD564UUcOHwDA+NSmONC4R8cg3skrSgz9EPMSrdsnR9PjRAU5oexa2U48rvf449vv43rfRJErs1BWrgIHRWf4MDvXsKrbzncL6Uxo/Kt17Dvl/+NV94oxlsNLBhvfBiF98QiQBGAULmFBZ8jeP73f8Tbn7XBKpUgIEIHtS4FwT3HUd7UCZEmE2GWCtSbO9FdY0DZS7/Fi4f0EIVn4+Ht6xEp92ZBXgzxQDVOXm+HRbMeBWu1kC76NELcMXj24vrhYuz/+V68yNrdmxXdUKZtxhe2JiNYxtqFfTWeU4gVoZANXkLZ4efx+z+9jbNtVkikAayEUyM8NQ06iQHnzhXjxT+xflLbBGOzL9aEaBGTrYUfa/tj7Sdwoq0Fngn52MoqHdG0fS4DUSxI3gojnhJIPKffbtyWZMTP2P5TsSY6AkGz9N0FhKrF357gyHIKSHntO7vEcAHCxSnrs61NiD/ZWczJqCFTIUyrg5b7Ph3nu4ojXQSMGlafS6Vs/2yffiOtt0SS1ku/xpEaE3y1O7FONmA7VhcZ40D1ERys7IIkcTt2prOGZJ9P3M04t8db7Y71oek+/5nFtHZZ70Q/cRXpsudOEdHOS0Dq0v9ct0sC0rsL3kiO1AwgNHOP0/1VkxhxZf8p1Jr9kfJELhJcWwtB3A6s7Fp2ES3L0EhAehfDr2Xt++luvFzXisCc9YiiIEUsGcsoomXBjwSkdxk8ZZ8oA9UugWginR/ysZehy3G2I1YpyyiiJQEpQRCEjSUr/QiCIJYTClQEQbg9FKgIgnB7KFARBOH2UKAiCMLtoUBFEITbQ4GKIAi3hwIVQRBuDwUqgiDcHgpUBEG4PRSoCIJweyhQEQTh9lCgIgjC7fEYGxv7e/aXTwRBEG4Iyv9/X2P83+dyoVsAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Następnie importowana jest biblioteka SimulationRT napisana na potrzeby dzisiejszych zajęć. Zawiera ona prosty algorytm śledzenia promieni. Pliki związane z tą biblioteką znajdują się w folderze simpleRT. Biblioteka ta wykorzystuje pakiet Cython, który pozwala na automatyczne tłumaczenie kodu napisanego w Pythonie na język C. Wykorzystanie biblioteki Cython pozwoliło na przyśpieszenie obliczeń około 1000-krotnie, kosztem przejrzystości i łatwości edycji kodu programu.\n",
    "\n",
    "![tree_struct.png](attachment:tree_struct.png)\n",
    "\n",
    "Po wczytaniu biblioteki tworzony jest podstawowy obiekt który będzie zawierał informacje o prowadzonych symulacjach - SimulationRT(). Obiekt ten przechowuje informacje o modelu, odbiornikach i nadajnikach oraz zawiera funkcie odpowiedzialne za wykonywanie obliczeń."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!python _CythonSetup.py build_ext --inplace    # Uncomment to build the library first\n",
    "from simpleRT import SimulationRT\n",
    "simulation = SimulationRT()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kolejnym etapem jest wczytanie modelu, dla którego będą przeprowadzane obliczenia. Modele są importowane z popularnego formatu obj. Format ten przypomina sposób opisu geometrii w programie CATT-Acoustic. Na początku pliku znajdują się informacje o wierzchołkach modelu, następnie definicje normalnych powierzchni i na końcu powierzchnie utworzone z wierzchołków.\n",
    "\n",
    "Po wczytaniu model jest przypisywany do symulacji."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simpleRT import Model3D\n",
    "import os\n",
    "\n",
    "folder = 'modele'\n",
    "basename = 'korytarz'\n",
    "filename = basename+'.obj'\n",
    "model3D = Model3D()\n",
    "\n",
    "os.chdir(folder)\n",
    "model3D.loadOBJ(filename)\n",
    "os.chdir('..')\n",
    "\n",
    "simulation.model3d = model3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Następne bloki kodu odpowiadają za wyświetlenie modelu sali i zadanie niektórych parametrów. Ich szczegółowa analiza nie jest konieczna, ponieważ odpowiadają one tylko za interakcję z użytkownikiem.\n",
    "\n",
    "Prawy przycisk myszy przesuwa model, lewy przycisk obraca model, kółko przybliża i oddala.\n",
    "\n",
    "UWAGA! W celu zatwierdzenia wprowadzanych wartości konieczne jest wciśnięcie klawisza enter!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "model_display_holder = widgets.Output(layout={'border': '1px solid black'})\n",
    "model_display_holder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from notebook_helpers.helpers import *\n",
    "display_model(model_display_holder, model3D, simulation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from notebook_helpers.helpers import *\n",
    "\n",
    "#### Create source, assign it to model, handle interactions with GUI\n",
    "from simpleRT import Source3D\n",
    "source = Source3D()\n",
    "simulation.source = source\n",
    "\n",
    "##### Create receiver, assign it to model, handle interactions with GUI\n",
    "from simpleRT import Receiver3D\n",
    "receiver = Receiver3D()\n",
    "simulation.receiver = receiver\n",
    "\n",
    "display_position_controls(simulation ,model3D, model_display_holder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "### Create table with material absorption and diffusion coefficients, handle GUI interactions\n",
    "from simpleRT.helpers import matloader as mat\n",
    "import ipysheet as ip\n",
    "\n",
    "mat_filename = os.path.join(folder,basename+'.mat')\n",
    "try:\n",
    "    model3D = mat.load_previous_materials(model3D, mat_filename)\n",
    "except:\n",
    "    print('Material file not found: %s' % s)\n",
    "    pass\n",
    "    \n",
    "df = model3D.export_materials_as_dict()\n",
    "\n",
    "sheet = ip.from_dataframe(df)\n",
    "sheet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ustal pozycje źródła i nadajnika zwryfikuj przyjęte wszpółczynniki pochłaniania i rozproszenia dźwięku. Po zadaniu parametrów modelu można przejść do wykonywania obliczeń. W poniższym bloku wykonywane są symulacje. Na początek uruchom symulacje z domyślnymi parametrami. \n",
    "\n",
    "Parametry:\n",
    " 1. no_rays - liczba promieni\n",
    " 2. fs - fs\n",
    " 3. sim_length - dlugość symulacji\n",
    " 4. max_reflection_order - określa liczbę odbić po której promień jest odrzucany,\n",
    " 5. brute_force - zmienia sposób przeszukiwania geometrii, jeżeli wartość jest ustawiona na false wykorzystywane jest drzewo BVH do szybszego wykrywania przecięć (więcej informacji na wykładzie), \n",
    " 6. reflection_model - zmienia sposób modelowania odbicia dźwięku, dostępne opcje `ReflectionModel.Phong` i `ReflectionModel.Base` (więcej informacji na wykładach)\n",
    " 7. intersection_model - zmienia sposób wykrywania przecięć promieni z modelem `RaySphereIntersection.Full` i `RaySphereIntersection.Simple`. W przypadku modelu Full uwzględniany jest czas propagacji promienia przez odbiornik, w przypadku modelu Simple uwzględniany jest tylko czas dotarcia do środka odbiornika. Korzystając z modelu simple nie jest konieczne wykrywanie maksimów na etapie auralizacji (więcej informacji na wykładach).\n",
    "\n",
    "W czasie trwania obliczeń zapoznaj się z opisem działania algorytmu znajdującym się poniżej."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simpleRT import ReflectionModel, RaySphereIntersection\n",
    "\n",
    "no_rays = 15000\n",
    "fs = 44100\n",
    "sim_length = 2\n",
    "max_reflection_order = 40\n",
    "brute_force = False\n",
    "reflection_model = ReflectionModel.Phong\n",
    "intersection_model = RaySphereIntersection.Full\n",
    "\n",
    "# Update and optionaly save materials\n",
    "updated_mat = ip.to_dataframe(sheet)\n",
    "model3D.import_materials_from_dict(updated_mat)\n",
    "#mat.save_materials(model3D, mat_filename)\n",
    "    \n",
    "# Define callback to display progress bar\n",
    "global f\n",
    "f = IntProgress(min=0, max=no_rays) # instantiate the bar\n",
    "display(f) # display the bar\n",
    "def callback(no_rays):\n",
    "    if no_rays%100==0:\n",
    "        global f\n",
    "        f.value += 100\n",
    "\n",
    "# Perform all required operations to start simulations\n",
    "def run_simulations(sim):\n",
    "    sim.initialize_simulation_parameters(no_rays, fs, sim_length, \\\n",
    "                                        max_reflection_order, brute_force, \n",
    "                                        reflection_model,intersection_model)\n",
    "    sim.initialize_rays() #create array of rays, set ray energy and position in the source, generate random ray direction\n",
    "    sim.run_simulation(callback)\n",
    "    return sim\n",
    "\n",
    "simulation = run_simulations(simulation)\n",
    "print(simulation.no_lost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Główny algorytm programu można przedstawić w postaci pseudokodu jako:\n",
    "\n",
    "---\n",
    "```\n",
    "Jeżeli liczba odbić jest mniejsza od maksymalnego rzędu odbicia:\n",
    "    Sprawdź przecięcie promienia z odbiornikiem\n",
    "    Wyszukaj z którą powierzchnią przecina się promień.\n",
    "    Jeżeli promień nie przecina się z powierzchnią\n",
    "        Uznaj promień za stracony\n",
    "    Sprawdż czy przecięcie z odbiornikiem jest przed przecięciem z powierzchnią\n",
    "        Jeżeli tak dodaj energię promienia do odbiornika\n",
    "    Przesuń promień do punktu następnego odbicia\n",
    "    Oblicz energię i kierunek promienia po odbiciu\n",
    "```\n",
    "---\n",
    "\n",
    "Szczegółowo przedstawia się on w sposób następujący:\n",
    "\n",
    "---\n",
    "```cython \n",
    "last_face = None\n",
    "for i in range(self.max_refl_ord):\n",
    "\n",
    "    #Find colision with receiver \n",
    "    t1 = np.Inf\n",
    "    rec_hit = libmath.ray_sphere_intersection(ray.position, ray.direction ,self.receiver.position, self.receiver.radius)\n",
    "    if rec_hit is not -1:\n",
    "        t1 = rec_hit[0] #entry point\n",
    "        t2 = rec_hit[1] #out_point\n",
    "\n",
    "    #Find colision with model geometry\n",
    "    if self.brute_force:\n",
    "        tmin, hit_face = self.brute_force_triangle_intersection(ray, self.model3d.faces, last_face)\n",
    "    else:         \n",
    "        tmin, hit_face = self.kd_tree_triang_intersection(ray, last_face)\n",
    "\n",
    "\n",
    "    # Test if model hit any face, if not it is lost\n",
    "    if hit_face is None:\n",
    "        ray.ray_lost()\n",
    "        self.no_lost+=1\n",
    "        return\n",
    "\n",
    "    # Test if colision with model is after colision with receiver\n",
    "    if tmin>t1 and t1>0:              \n",
    "        entry_sample = int(min((ray.total_dist + t1)/self.dv, self.fs*self.sim_len))\n",
    "        exit_sample = int(min((ray.total_dist + t2)/self.dv + 1, self.fs*self.sim_len))\n",
    "        ds = exit_sample-entry_sample\n",
    "\n",
    "        # If simlified intersection model is selected - using only middle intersection sample\n",
    "        if self.ray_sphere_int == RaySphereIntersection.Simple:\n",
    "            entry_sample = int((exit_sample-entry_sample)/2)\n",
    "            exit_sample = entry_sample+1\n",
    "            ds=1\n",
    "\n",
    "        # add energy to result\n",
    "        self.result[:,entry_sample:exit_sample] += np.transpose(np.tile(ray.energy,(ds,1)))\n",
    "        ray.receiver_hit = True\n",
    "\n",
    "    # Move ray to face hit position, update energy\n",
    "    dif_coeff = self.model3d.materials[hit_face.mat_name]['diffusion_coeff']\n",
    "    abs_coeff = self.model3d.materials[hit_face.mat_name]['absorption_coeff']\n",
    "    last_face = hit_face\n",
    "    ray.position = libmath.move_ray(ray.position, ray.direction, tmin)\n",
    "    ray.total_dist += tmin\n",
    "    ray.energy*=(1-abs_coeff)\n",
    "\n",
    "    # calculate reflection based on chosen reflection model\n",
    "    if self.reflection_model == ReflectionModel.Base:\n",
    "        ray.direction = libmath.base_reflection(ray.direction, self.model3d.normals[hit_face.normal_idx], dif_coeff)\n",
    "    elif self.reflection_model == ReflectionModel.Phong:\n",
    "        ray.direction = libmath.phong_reflection(ray.direction, self.model3d.normals[hit_face.normal_idx], dif_coeff)\n",
    "    ray.refl_order+=1\n",
    "```\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "W przypadku kiedy odbiornik ma określoną objętość (`intersection_model = RaySphereIntersection.Full`, konieczne jest podzielenie wyniku symulacji przez objętość odbiornika. W ten sposób wyniki są niezależne od objętości odbiornika."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "receiver_volume = (4/3*np.pi*simulation.receiver.radius**3)\n",
    "simulation.result = simulation.result/receiver_volume\n",
    "simulation.result_spatial = simulation.result_spatial/receiver_volume"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nastęnie możemy wyświetlić wyniki symulacji. Zmienna simulation.results zawiera echogramy dla poszczególnych częstotliwości. Wynik dla określonej częstotliwości można otrzymać przez `simulations.result[idx]`, \n",
    "gdzie `idx`to index częstotliwości. Pierwszy wykres przedstawia zarejestrowaną energię docierającą do odbiornika (kwadrat ciśnienia akustycznego), drugi przedstawia echogram w skali logarytmicznej."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e_tot = np.sum(simulation.result, axis=0)\n",
    "\n",
    "plt.figure(figsize=(10, 3),num='Sumaryczny echogram, skala liniowa')\n",
    "plt.plot(simulation.time,e_tot)\n",
    "\n",
    "plt.figure(figsize=(10, 3),num='Sumaryczny echogram, skala logarytmiczna')\n",
    "zeros = np.zeros(len(e_tot))\n",
    "Lp = 10*np.log10(e_tot/4e-10)\n",
    "Lp[Lp==-np.inf]=0\n",
    "plt.plot(simulation.time,Lp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Na koniec możemy zapisać wyniki. Wynik symulacji zapisywany jest do pliku, którego nazwa składa się z nazwy modelu, liczby promieni w symulacjach oraz rozszerzenia \".res\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pickle.dump( simulation, open( os.path.join(folder,basename+str(no_rays)+\".res\"), \"wb\" ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Auralizacjia\n",
    "\n",
    "Pierwszym krokiem jest wczytanie wyników symulacji. Jeżeli w danej sesji zostały przeprowadzone symulajce, wybierany jest wynik przechowywany w pamięci, w przeciwnym razie wyniki wczytywane są z dysku. Konieczne może być też ponowne zaimportowanie wykorzystywanych bibliotek."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import cython\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "from simpleRT import SimulationRT\n",
    "\n",
    "#del simulation\n",
    "\n",
    "fs=44100\n",
    "if 'simulation' not in locals():\n",
    "    no_rays = 1500000\n",
    "    folder='modele'\n",
    "    basename = 'halaWidowiskowa'\n",
    "    simulation = pickle.load( open( os.path.join(folder,basename+str(no_rays)+\".res\"), \"rb\" ) )\n",
    "    \n",
    "    #simulation = SimulationRT()\n",
    "    #spos = [0.8,4,15,1]\n",
    "    #rpos = [2.9,1.1,3]\n",
    "    #fs=44100\n",
    "    #simulation.import_ISIMPA_results('ISimpa/Sound level.csv', 'ISimpa/Punctual receiver intensity.csv',spos,rpos,fs)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Następnie na podstawie echogramów chcemy uzyskać odpowiedzi impulsowe. Poniżej przedstawione są dwa podejścia do przekształcenia echogramu w odpowiedź impulsową nadającą się do auralizacji. \n",
    "\n",
    "### Wariant 1 - wykrywanie maksimów i bezpośredni splot\n",
    "Ponieważ symulowane echogramy są próbkowane z częstotliwością 44100Hz, przejście promienia przez odbiornik trwa więcej niż jedną próbkę. Na potwierdzenie powiększ prążek odpowiadający za dźwięk bezpośredni na wykresach powyżej. Zauważ, że energia jest rozłożona na wiele próbek. Dla poprawnej auralizacji w tym wariancie konieczne jest założenie, że odbiornik jest punktem, a poszczególne promienie docierają do niego tylko w pojedynczych, określonych chwilach czasu. Wtedy wynik symulacji możemy potraktować jako ciąg delt diraca, z odpowiednią amplitudą i opóźnieniami (jedna delta - jedno docierające do odbiornika odbicie). Taki ciąg możemy splatać z naszym sygnałem.\n",
    "\n",
    "Przeprowadź ponowne obliczenia dla odbiornika o promieniu 5cm, porównaj wyniki z otrzymanymi w poprzedniej symulacji. Nie zapisuj wyników!\n",
    "\n",
    "Ok, jak widać znaczniej mniej promieni trafiło do odbiornika, a wyniki są wyraźnie mocniej zaszumione. W celu zwększenia dokładności możliwe jest zwiększenie liczby promieni, ale to wiąże się z dłuższym czasem obliczen. Zamiast tego można sprowadzić wyniki do punktu przez wykrywanie maksimów.\n",
    "\n",
    "Funkcja `peak_detection()` w pierwszej kolejności znajduje indeksy maksimów i zapisuje je do zmiennej `peaks`, nastęnie tworzony jest wektor zer o długości odpowiadającej długości echogramu, na koniec wykryte maksima są przepisywane w odpowiednich miejscach do wektora zer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import find_peaks\n",
    "import numpy as np\n",
    "\n",
    "def peak_detection(result):\n",
    "    peaks, _ = find_peaks(result, height=0)\n",
    "    pf_res = np.zeros(len(result))\n",
    "    pf_res[peaks] = result[peaks]\n",
    "    return pf_res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mając echogramy po wykrywaniu impulsów koniecznie jest uwzględnienie jakiego pasma częstotliwości one dotyczą. Wtym celu wykonywana jest filtracja echogramu z wkorzystaniem filtrów oktawowych. Poniższy blok służy do wygenerowania banku filtrów."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import find_peaks\n",
    "from pyfilterbank import FractionalOctaveFilterbank\n",
    "\n",
    "#Initialize octave filter bank\n",
    "ofb = FractionalOctaveFilterbank(fs, order=4,nth_oct=1,  start_band=-3, end_band=2)\n",
    "print(ofb.center_frequencies)\n",
    "\n",
    "def filter_octave(data, ofb, idx):\n",
    "    y, states = ofb.filter(data)\n",
    "    return y[:,idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mając gotowe filtry można przefiltrować echogramy dla każdej częstotliwości i je ze sobą zsumować, żeby otrzymać modelowaną odpowiedź impulsową."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ZMIENNA KONTROLUJĄCA CZY WYKRYWANIE MAKSIMÓW JEST AKTYWNE\n",
    "do_peak_detection = True\n",
    "\n",
    "# Variable creation\n",
    "pf_filter = np.empty(simulation.result.T.shape, dtype=float)\n",
    "peak_det = np.empty(simulation.result.T.shape, dtype=float)\n",
    "\n",
    "# Convert energy to pressure\n",
    "result_press = np.sqrt(simulation.result)\n",
    "\n",
    "# For each freq range detect peaks anf filter the echogram\n",
    "for idx in range(len(simulation.result)):\n",
    "    if do_peak_detection:\n",
    "        peak_det[:,idx] = peak_detection(result_press[idx])\n",
    "        pf_filter[:,idx] = filter_octave(peak_det[:,idx], ofb, idx)\n",
    "    else:\n",
    "        pf_filter[:,idx] = filter_octave(result_press[idx], ofb, idx)\n",
    "        \n",
    "# Sumowanie wszystkich częstotliwości\n",
    "ir_w1 = np.sum(pf_filter,axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poniższy blok kodu wyświetla etapy obliczania odpowiedzi impulsowej oraz ostateczny wynik"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# Wyświetlanie wyników\n",
    "plt.figure(figsize=(10, 5),num='Proces wyznaczania odpowiedzi impulsowej')\n",
    "ax = plt.subplot(311)\n",
    "ax.set_title(\"Echogram\")\n",
    "plt.plot(result_press[3,2000:8000])\n",
    "ax = plt.subplot(312)\n",
    "ax.set_title(\"Echogram po wykrywaniu impulsów\")\n",
    "plt.plot(peak_det[2000:8000,3])\n",
    "ax = plt.subplot(313)\n",
    "ax.set_title(\"Echogram po wykrywaniu impulsów i filtracji\")\n",
    "plt.plot(pf_filter[2000:8000,3])\n",
    "plt.tight_layout()\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10, 5),num='Obliczona odpowiedź impulsowa w skali liniowej, wariant 1')\n",
    "plt.plot(simulation.time,ir_w1)\n",
    "\n",
    "plt.figure(figsize=(10, 5),num='Obliczona odpowiedź impulsowa w skali logarytmicznej, wariant 1')\n",
    "Lp = 10*np.log10(ir_w1**2/(2e-5)**2)\n",
    "Lp[Lp==-np.inf]=0\n",
    "plt.plot(simulation.time,Lp)\n",
    "plt.ylim(np.max(Lp)-90,np.max(Lp)+5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wariant 2 - Echogram jako obwiednia przebiegu kwadratu ciśnienia akustycznego\n",
    "\n",
    "W symulacjach metodą geometryczną śledzimy rozchodzenie się energii w modelu. Energia jest proporcjonalna do kwadratu ciśnienia akustycznego. Ponieważ w symulacjach zaniedbujemy fazę otrzymane echogramy możemy rozpatrywać jako obwiednię przebiegu kwadratu ciśnienia akustycznego. W takim razie aby odtworzyć przebieg ciśnienia możemy wygenerować szum biały (płaskie widmo), na który nałożymy obwiednię otrzymaną z symulacji. W celu pozbycia się gwałtownych skoków w echogramie wykonywany jest splot z krótkim oknem Gaussa. \n",
    "\n",
    "Uwaga - szerokość okna Gaussa regulowana parametrem `std_ms` ma duży wpływ na ostatecznie uzyskiwane brzmienie!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal\n",
    "from scipy.signal import fftconvolve\n",
    "\n",
    "# Convert energy to ~pressure\n",
    "result_press = np.sqrt(simulation.result)\n",
    "\n",
    "# Variable creation\n",
    "ir_gaussian = np.empty(simulation.result.T.shape, dtype=float)\n",
    "noise_filter = np.empty(simulation.result.T.shape, dtype=float)\n",
    "ir_bands = np.empty(simulation.result.T.shape, dtype=float)\n",
    "\n",
    "# Generate white noise\n",
    "noise = np.random.rand(len(simulation.time))*2 - 1\n",
    "noise = np.multiply(noise, signal.tukey(len(simulation.time), alpha=0.005))\n",
    "\n",
    "# Generate gaussian window to smooth the echogram\n",
    "glen_ms = 10\n",
    "gstd_ms = 1\n",
    "glen_sample = int(fs*glen_ms/1000)\n",
    "gstd_sample = int(fs*gstd_ms/1000)\n",
    "gaussian = signal.gaussian(glen_sample, std=gstd_sample)\n",
    "\n",
    "# For each freq range filter noise and apply envelope\n",
    "for idx in range(len(simulation.result)):\n",
    "    ir_gaussian[:,idx] = fftconvolve(result_press[idx], gaussian, mode = 'same') # Convovle ir with gaussian to smooth the result\n",
    "    noise_filter[:,idx] = filter_octave(noise, ofb, idx)                         # Filter noise using octave filters\n",
    "    ir_bands = np.multiply(ir_gaussian,noise_filter)                             # Apply the envelope by signals multiplication (element by element)\n",
    "    \n",
    "ir_w2 = np.sum(ir_bands,axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poniższy blok kodu wyświetla etapy obliczania odpowiedzi impulsowej oraz ostateczny wynik"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# Wyświetlanie wyników\n",
    "plt.figure(figsize=(10, 5),num='Proces wyznaczania odpowiedzi impulsowej')\n",
    "ax = plt.subplot(411)\n",
    "ax.set_title(\"Echogram\")\n",
    "plt.plot(result_press[3,2000:8000])\n",
    "ax = plt.subplot(412)\n",
    "ax.set_title(\"Echogram spleciony z oknem Gaussa\")\n",
    "plt.plot(ir_gaussian[2000:8000,3])\n",
    "ax = plt.subplot(413)\n",
    "ax.set_title(\"Przefiltrowany szum\")\n",
    "plt.plot(noise_filter[2000:8000,3])\n",
    "ax = plt.subplot(414)\n",
    "ax.set_title(\"Szum z nałożoną obwiednią\")\n",
    "plt.plot(ir_bands[2000:8000,3])\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.figure(figsize=(10, 3),num='Obliczona odpowiedź impulsowa w skali liniowej, wariant 2')\n",
    "plt.plot(simulation.time,ir_w2)\n",
    "\n",
    "plt.figure(figsize=(10, 3),num='Obliczona odpowiedź impulsowa w skali logarytmicznej, wariant 2')\n",
    "Lp = 10*np.log10(ir_w2**2/(2e-5)**2)\n",
    "Lp[Lp==-np.inf]=0\n",
    "plt.plot(simulation.time,Lp)\n",
    "plt.ylim(np.max(Lp)-90,np.max(Lp)+5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Splot odpowiedzi impulsowej z sygnałem audio\n",
    "\n",
    "Mając wyznaczone odpowiedzi impulsowe można dokonać splotu z wybranym sygnałem.\n",
    "Funkcja fftconvolve wykonuje szybki splot z wykorzytaniem transformacji Fouriera."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import soundfile as sf\n",
    "from scipy.signal import fftconvolve\n",
    "\n",
    "#Load audio sample and convolve\n",
    "audio_folder = 'audio'\n",
    "audio_file = 'mowa.wav'\n",
    "data, samplerate = sf.read(os.path.join(audio_folder,audio_file))\n",
    "convolved_w1 = fftconvolve(data[:], ir_w1)\n",
    "convolved_w2 = fftconvolve(data[:], ir_w2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Następnie można przesłuchać wyniki auralizacj."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import sounddevice as sd\n",
    "import IPython.display as ipd\n",
    "\n",
    "#sd.play(convolved, fs)\n",
    "print('Wariant 1:')\n",
    "display(ipd.Audio(convolved_w1, rate=fs))\n",
    "print('Wariant 2:')\n",
    "display(ipd.Audio(convolved_w2, rate=fs))\n",
    "\n",
    "print('Wariant 2 IR:')\n",
    "display(ipd.Audio(ir_w2, rate=48000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Uprzestrzennienie\n",
    "\n",
    "Metoda promieniowa oprócz przebiegu energii może dostarczać informacji na temat kierunku z którego promienie docierają do odbiornika. Na tej podstawie można odtworzyć wrażenia przestrzenne związane z odsłuchem w symulowanym pomieszczeniu. Poniważ docelową formą odsłuchu na naszych zajęciach jest wykorzstyanie słuchawek, w celu nadanie przestrzenności symulowanej odpowiedzi impulsowej konieczne jest wykorzystanie funkcji HRTF.\n",
    "\n",
    "Poniższy fragment kodu otwiera pliki HRTF zapisane w otwartym formacie .sofa, oraz wyświetla zapisane w pliku pozycje odbiorników i przykładową odpowiedź impulsową dla prawego i lewego ucha."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits import mplot3d\n",
    "\n",
    "import numpy as np\n",
    "import sofa\n",
    "\n",
    "def plot_coordinates(coords, title):\n",
    "    x0 = coords\n",
    "    n0 = coords\n",
    "    fig = plt.figure(figsize=(5, 5))\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    q = ax.quiver(x0[:, 0], x0[:, 1], x0[:, 2], n0[:, 0],\n",
    "                  n0[:, 1], n0[:, 2], length=0.1)\n",
    "    plt.xlabel('x (m)')\n",
    "    plt.ylabel('y (m)')\n",
    "    plt.title(title)\n",
    "    return q\n",
    "\n",
    "# Open and print sofa Metadata\n",
    "HRTF_path = \"./hrtfs/mit_kemar_normal_pinna.sofa\"\n",
    "HRTF = sofa.Database.open(HRTF_path)\n",
    "\n",
    "# plot Source positions\n",
    "source_positions = HRTF.Source.Position.get_values(system=\"cartesian\")\n",
    "receiver_positions = HRTF.Receiver.Position.get_values(system=\"cartesian\")\n",
    "emitter_positions = HRTF.Emitter.Position.get_values(system=\"cartesian\")\n",
    "plot_coordinates(source_positions, 'Source positions');\n",
    "\n",
    "source_positions_sph = HRTF.Source.Position.get_values(system=\"spherical\")\n",
    "\n",
    "# plot Data.IR at M=5 for E=0\n",
    "measurement = 5\n",
    "emitter = 0\n",
    "legend = []\n",
    "\n",
    "t = np.arange(0,HRTF.Dimensions.N)*(1/HRTF.Data.SamplingRate.get_values(indices={\"M\":measurement}))\n",
    "\n",
    "plt.figure(figsize=(10, 2))\n",
    "for receiver in np.arange(HRTF.Dimensions.R):\n",
    "    plt.plot(t, HRTF.Data.IR.get_values(indices={\"M\":measurement, \"R\":receiver, \"E\":emitter}))\n",
    "    legend.append('Receiver {0}'.format(receiver))\n",
    "plt.title('HRIR at M={0} for emitter {1}'.format(measurement, emitter))\n",
    "plt.legend(legend)\n",
    "plt.xlabel('$t$ in s')\n",
    "plt.ylabel(r'$h(t)$')\n",
    "plt.grid()\n",
    "\n",
    "HRTF.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poniższy blok kodu przygotowuje dane do dalszych obliczeń. Wykonywana są konwersje układów współrzędnych i przeliczanie wyników z radianów na stopnie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from aura.helpers import cart2sph, rotation_relative2source\n",
    "\n",
    "# Convert coordinates from [x z -y] to [x y z]\n",
    "spos=simulation.source.position.asArray()\n",
    "spos=[spos[0], -spos[2], spos[1]]\n",
    "rpos=simulation.receiver.position.asArray()\n",
    "rpos=[rpos[0], -rpos[2], rpos[1]]\n",
    "\n",
    "# Calculate vector pointing from the receiver to the source\n",
    "source_vector = np.subtract(spos,rpos)\n",
    "source_vector = source_vector/np.linalg.norm(source_vector)\n",
    "sv_spherical = np.array(cart2sph(*source_vector))\n",
    "sv_spherical = np.rad2deg(sv_spherical)\n",
    "\n",
    "# Rotate simulation results so that source is at 0,0 or angles given in rotation_relative2source\n",
    "spatial_info = np.sum(simulation.result_spatial,axis=0).T\n",
    "spatial_info = np.vstack((spatial_info[:,0], -spatial_info[:,2], spatial_info[:,1])).T\n",
    "#spatial_spherical = np.array([ for vector in spatial_info.T])\n",
    "#ssa = spatial_spherical_aligned = rotation_relative2source(spatial_spherical,sv_spherical,0,0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "W kolejnym bloku obliczane jest faktyczne uprzestrzennienie odpowiedzi impulsowej. Parametry `azim` i `elev` odpowiadają za ustawienie głowy względem odbiornika. Wartości `0,0` oznaczają, że słuchacz skierowany jest bezpośrednio na odbiornik. Następnie wybierana jest odpowiedź impulsowa wyznaczona w poprzednim kroku analiz (domyślenie 2-gi wariant) i wczytywane są dane z pliku HRTF. W razie konieczniści odwpowiedź impulsowa jest przepróbkowana.\n",
    "\n",
    "Właściwe obliczenia zaczynają się od ustalenia długości okna analizy - domyślnie 20 ms. Okno to określa czas przez który sygnał będzie przypisany do konkrentego kierunku - (krótsze okno - częstsza zmiana kierunków). Na podstawie długości okna wyznaczane są niezbędne do obliczeń parametry i inicjalizowany jest obekt przechowujący wynik `binaural_ir`. \n",
    "\n",
    "Dla każdego okna (pętla `for`) wyznaczany jest średni kierunek z którego dochodził sygnał. Kierunek ten jest przeliczany na współrzędne sferyczne i stopnie oraz przekształcany w odniesieniu do pozycji źródła (tak aby kąty 0,0 oznaczały, że dźwięk dochodzi ze źródła). Znając kierunek z którego dochodzi sygnał wybierany jest najbliższy HRTF i przeprowadzany jest splot danego okna odpowiedzi impulsowej z wybranym filtrm HRTF. \n",
    "\n",
    "Schemat filtracji jest identyczny jak w algorytmie szybkiego splotu overlap-add, z tą różnicą, że dla każdego okna wykorzystywany jest inny filtr. Wykorzystanie pełnego wyniku splotu (overlap) i dodawanie zachodzącej na siebie części pozwala na uniknięcie problemów z nagłą zmianą filtra i skokami sygnału."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import resample_poly, convolve\n",
    "from scipy import spatial\n",
    "from fractions import Fraction\n",
    "\n",
    "# Set head rotation\n",
    "azim=0\n",
    "elev=0\n",
    "\n",
    "# Select impulse response to spacialize\n",
    "ir = ir_w2\n",
    "\n",
    "# Load required info from HRTF\n",
    "HRTF_path = \"./hrtfs/mit_kemar_normal_pinna.sofa\"\n",
    "HRTF = sofa.Database.open(HRTF_path)\n",
    "hrtf_len = HRTF.Dimensions.N\n",
    "hrtf_fs = HRTF.Data.SamplingRate.get_values(indices={\"M\":0})\n",
    "hrtf_source_positions = HRTF.Source.Position.get_values(system=\"spherical\")\n",
    "hrtf_data = HRTF.Data\n",
    "\n",
    "# Build kd tre for fast nearest point searching\n",
    "tree = spatial.KDTree(hrtf_source_positions[:,0:-2])\n",
    "\n",
    "# Check sampling rates, resample if needed\n",
    "if simulation.fs != hrtf_fs:\n",
    "    Print('Resample required, resampling simulation results')\n",
    "    ratio = fraction(hrtf_fs,simulation.fs)\n",
    "    ir = resample_poly(ir,ratio.nominator,ratio.denominator)\n",
    "    spatial_info = resample_poly(spatial_info,ratio.nominator,ratio.denominator)\n",
    "\n",
    "# Set windowing parameters and initialize result vector\n",
    "wlen_ms = 15\n",
    "wlen_sample = int(hrtf_fs*wlen_ms/1000)\n",
    "no_loops = int(len(ir)/wlen_sample)\n",
    "binaural_ir = np.zeros(((no_loops*wlen_sample)+hrtf_len-1,2))\n",
    "\n",
    "# Calculate spatial IR\n",
    "for i in range(no_loops):\n",
    "    window_data = ir[i*wlen_sample:(i+1)*wlen_sample]\n",
    "    avr_dir = np.sum(spatial_info[i*wlen_sample:(i+1)*wlen_sample],axis=0)\n",
    "    avr_dir_sph = np.rad2deg(cart2sph(*avr_dir))\n",
    "    avr_dir_rel = rotation_relative2source([avr_dir_sph],sv_spherical,azim,elev)\n",
    "    _, index = tree.query(avr_dir_rel[0,0:-2])\n",
    "    #print(avr_dir_rel)\n",
    "    \n",
    "    hrtf_l = HRTF.Data.IR.get_values(indices={\"M\":index, \"R\":1, \"E\":0})\n",
    "    hrtf_r = HRTF.Data.IR.get_values(indices={\"M\":index, \"R\":0, \"E\":0})\n",
    "    \n",
    "    # Overlap add\n",
    "    binaural_ir[i*wlen_sample:(i+1)*wlen_sample + hrtf_len-1,0]+=convolve(window_data, hrtf_l, mode='full')\n",
    "    binaural_ir[i*wlen_sample:(i+1)*wlen_sample + hrtf_len-1,1]+=convolve(window_data, hrtf_r, mode='full')\n",
    "    \n",
    "HRTF.close()\n",
    "    \n",
    "# Plot spatial IR\n",
    "t = np.arange(0,len(binaural_ir))*(1/hrtf_fs)\n",
    "plt.figure(figsize=(10, 5),num='Binauralna odpowiedź impulsowa')\n",
    "plt.plot(t,binaural_ir[:,0])\n",
    "plt.plot(t,binaural_ir[:,1])\n",
    "\n",
    "#import sounddevice as sd\n",
    "import IPython.display as ipd\n",
    "\n",
    "print('Binauralna odpowiedź impulsowa')\n",
    "display(ipd.Audio([binaural_ir[:,0], binaural_ir[:,1]] , rate=48000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nastęnie można wykonać splot binauralnej odpowiedzi impulsowej z wybranym sygnałem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import soundfile as sf\n",
    "from scipy.signal import fftconvolve\n",
    "\n",
    "#Load audio sample and convolve\n",
    "audio_folder = 'audio'\n",
    "audio_file = 'mowa.wav'\n",
    "audio, samplerate = sf.read(os.path.join(audio_folder,audio_file))\n",
    "\n",
    "# Check sampling rates, resample if needed\n",
    "if samplerate != hrtf_fs:\n",
    "    Print('Resample required, resampling simulation results')\n",
    "    ratio = fraction(samplerate,hrtf_fs.fs)\n",
    "    binaural_ir = resample_poly(ir,ratio.nominator,ratio.denominator)\n",
    "\n",
    "convolved_l = fftconvolve(audio[:], binaural_ir[:,0])\n",
    "convolved_r = fftconvolve(audio[:], binaural_ir[:,1])\n",
    "\n",
    "#import sounddevice as sd\n",
    "import IPython.display as ipd\n",
    "\n",
    "#sd.play(convolved, fs)\n",
    "print('Wynik splotu z binauralną odpowiedzią impulsową')\n",
    "display(ipd.Audio([convolved_l, convolved_r] , rate=samplerate))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Zadania\n",
    " ### Zadanie 1\n",
    "   Porównaj wyniki auralizacji z wykorzystaniem wariantu 1 i 2. Które podejście daje lepsze rezultaty, dlaczego? Wykorzystaj modele \"korytarz\" i \"halaWidowiskowa\". (W folderze z modelami są już wyniki dla dużych liczb promieni).\n",
    "\n",
    " ### Zadanie 2\n",
    "   Porównaj wyniki auralizacji dla Korytarza z wynikiem splotu ze zmierzoną odpowiedzią impulsową (2gie zajęcia) i rzeczywistym nagraniem z tego pomieszczenia.\n",
    "    \n",
    " ### Zadanie 3\n",
    "   Jaki wpłw na wynik auralizacji ma wykrywanie maksimów - wykonaj auralizację w wariance 1 zmieniając wartość zmiennej `do_peak_detection` na False (możesz pominąć etap uprzestrzennienia)\n",
    "    \n",
    " ### Zadanie 4\n",
    "   Sprawdź jak zmieniają się wyniki auralizacji w zaleźności od liczby promieni i liczby odbić w symulacjach. (W folderze z plikami są już wyniki dla dużych liczb promieni)\n",
    "    \n",
    " ### Zadanie 5\n",
    "   Sprawdź jak zmieniają się wyniki auralizacji w zaleźności od szerokości okna Gaussa w wariancie drugim (zmienna `gstd_sample`).\n",
    "    \n",
    " ### Zadanie 6*\n",
    "   Zastanów się jak korzystając z powyższego programu wykonać auralizację stereofonicznego systemu nagłośnienia, czy możliwe jest uwzględnienie funkcji HRTF?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "e54e0e078e955b1b0746e886d91677b8823d0ccc982f8f55da62263e84df05a1"
  },
  "kernelspec": {
   "display_name": "Python 3.8.8 64-bit ('base': conda)",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}