examples/002 Materials.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "logging.getLogger('trimesh').disabled = True\n",
    "logging.getLogger('shapely.geos').disabled = True\n",
    "logging.getLogger('matplotlib').setLevel(logging.CRITICAL)\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from pvtrace import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Materials\n",
    "\n",
    "In *pvtrace* material must have a refractive index,\n",
    "\n",
    "```python\n",
    "air = Material(refractive_index=1.0)\n",
    "plastic = Material(refractive_index=1.5)\n",
    "```\n",
    "\n",
    "Multiple **components** can be added to materials to change the optical properties of the *volume*.\n",
    "\n",
    "There are three components types available: absorber, scatterer, luminophore.\n",
    "\n",
    "#### Absorber\n",
    "\n",
    "The `Absorber` uses the Beer-Lambert law, $I = I_o e^{-\\alpha x}$ to determine the propability that rays are non-radiatively absorbed as they travel along a path. Here the `coefficient` parameter it the absorption coefficient per unit length $\\alpha$.\n",
    "\n",
    "```python\n",
    "opaque = Material(\n",
    "    refractive_index=1.5\n",
    "    components=[\n",
    "        Absorber(coefficient=10.0)  # cm-1\n",
    "    ]\n",
    ")\n",
    "```\n",
    "\n",
    "The `components` attribute is a `list` of optical properties within the materials volume, such as absorption of light using the `Absorber` component type.\n",
    "\n",
    "#### Scatterer\n",
    "\n",
    "Change to `Scatterer` for a material which volume scatters light,\n",
    "\n",
    "```python\n",
    "opaque = Material(\n",
    "    refractive_index=1.5\n",
    "    components=[\n",
    "        Scatterer(\n",
    "            coefficient=1.0,  # e.g. units: cm-1\n",
    "            quantum_yield=1.0\n",
    "       )\n",
    "    ]\n",
    ")\n",
    "```\n",
    "\n",
    "The scattering coefficient is per unit length. To make scattering events non-radiative using the `quantum_yield` keyword. Scattering is assumed to be `isotropic`, change the phase function for different behaviour,\n",
    "\n",
    "```python\n",
    "opaque = Material(\n",
    "    refractive_index=1.5\n",
    "    components=[\n",
    "        Scatterer(\n",
    "            coefficient=1.0, # cm-1\n",
    "            quantum_yield=1.0,\n",
    "            phase_function=lambda : henyey_greenstein(g=0.8),\n",
    "        )  \n",
    "    ]\n",
    ")\n",
    "```\n",
    "\n",
    "Henyey-Greenstein is a useful one-parameter phase function which moves between favoring back scattering (g < 0.0), isotropic (g = 0.0) and forward scattering (g > 0.0).\n",
    "\n",
    "\n",
    "#### Luminophore\n",
    "\n",
    "*pvtrace* trace was originally designed for the study of luminescent solar concentrators, in which organic dyes or nanocrystals (luminophores) are embedded in plastic sheets. They absorb rays and (possibly) remit them at longer wavelength. To use the `Lumiophore` component you will need to provide an absorption and emission spectrum.\n",
    "\n",
    "\n",
    "```python\n",
    "\n",
    "# Arrays of data\n",
    "absorption_spectrum = ...  # like np.column_stack((x, y)), spectrum cm-1 nm-1\n",
    "emission_spectrum = ...  # as above, spectrum nm-1\n",
    "\n",
    "opaque = Material(\n",
    "    refractive_index=1.5\n",
    "    components=[\n",
    "        Luminophore(\n",
    "            coefficient=absorption_spectrum,\n",
    "            emission=emission_spectrum,\n",
    "            quantum_yield=1.0\n",
    "            phase_function=isotropic  # can also have phase functions\n",
    "        )\n",
    "    ]\n",
    ")\n",
    "```\n",
    "\n",
    "---\n",
    "**Units**\n",
    "\n",
    "*pvtrace* does *not* assume units, they could be metres, kilometres or lightyears. But you **must** pick one. Define absorption and scattering coefficient as inverses of the same unit that defines the dimensions of objects in the scene. For example, if the scene is defined in units of $\\textit{cm}$ then attenuation coefficients should be $\\textit{cm}^{\\textit{-1}}$.\n",
    "\n",
    "**Using experimental data**\n",
    "\n",
    "Don't use experimental absorption coefficient data directly with *pvtrace*. For best results fit an equation to the long wavelength tail of the absorption and use a combination of real data (at short wavelength) and the fit (at long wavelengths).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example with luminescent materials\n",
    "\n",
    "*pvtrace* comes with a default spectrum of the dye Lumogen F Red,\n",
    "\n",
    "```python\n",
    "from pvtrace.data import lumogen_f_red_305\n",
    "```\n",
    "\n",
    "Use the `absorption` and `emission` functions,\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXl8XNV5979H+75vlmVbsmx5320MYSlLEkwoJCmkQJM0ZCEf0oa3LU1CyNJspKQJffOmbVraQEN2IGQj7CSGAGazLSTLkmVLsmTt+zIaLSPNzHn/uBpZy4zmzsy9d+69zPfz8Wc8M3fOfX56zn3mzLnPeY6QUhIjRowYMexFXLQNiBEjRowY2hML7jFixIhhQ2LBPUaMGDFsSCy4x4gRI4YNiQX3GDFixLAhseAeI0aMGDYkFtxjxIgRw4bEgnuMGDFi2JBYcI8RI0YMG5IQrRMXFBTI8vLyaJ0+RowYMSzJ8ePHB6WUhcGOi1pwLy8v59ixY2F9tqWlhcrKSo0tig4xLebDLjogpsWsRKJFCHFOzXGWnJbJy8uLtgmaEdNiPuyiA2JazIoRWiwZ3CcnJ6NtgmbEtJgPu+iAmBazYoQWSwb3uDhLmu2XmBbzYRcdENNiVozQErU590hITEyMtgmaEdNiPuyiA/TXMjs7S2dnJ9PT07qeB8Dj8TA6Oqr7eYxAjZaUlBTKysrC9qElg7vT6aSgoCDaZmhCTIv5sIsO0F9LZ2cnmZmZlJeXI4TQ7TwALpeL5ORkXc9hFMG0SCkZGhqis7OTioqKsM4R9LeBEOJ/hRD9QoiTAd4XQoh/E0I0CyFOCCH2hmVJCNjlwoOYFjNiFx2gv5bp6Wny8/N1D+wACQmWHIv6JZgWIQT5+fkR/SJSM/HzEHBohfevATbO/fsk8F9hW6OSzs5OvU9hGDEt5sMuOsAYLUYEdoCZmRlDzmMEarRE+ncNGtyllC8Bwysc8l7gx1LhdSBHCLEqIquCsGHDBj2bN5SYFvNhVh3j07P88lgHM26vug84B9joOQNT9pinTklJibYJmmGEFi1u2a4GOhY875x7bRlCiE8KIY4JIY719PQwODhIT08PXV1djIyM0NLSwtTUFA0NDXi9XqqrqwE4fvw4ANXV1Xi9Xl555RWmpqZoaWlhZGSErq4ufO21tbXhdDppbGzE7XZTW1u7qA3fY11dHS6Xi6amJhwOB+3t7fT399Pf3097ezsOh4OmpiZcLhd1dXV+26itrcXtdtPY2IjT6aStrS1kTS+++CJer5eGhgbLa3rxxRcX+cmqmk6ePBmw70VL06tvHuPjDx3js4+d4JfHO4Jqmnz+XrhvA/G//DATP/wLvO5ZXTT5RqCTk5N4vV6mp6fxeDy4XC5mZ2eZnZ3F5XLh8XiYnp7G6/XOpwFOTEwsepycnERKydTU1KI2ZmZmmJmZYXJycr6NqakppJTzn3344YcRQlBdXT3fxuHDh3nPe96D2+2eb2OhPb42Atnz7W9/m4ceeigkTV/4whe47777gmq68sor6evrw+12B9Tke1zqJ7UINRtkCyHKgSeklNv9vPcE8C0p5Stzz/8I3CWlXHH56f79+2W4K1RjxHi7cecjNfympouCjGTy0pJ45u8vDfyzfawT/n0flF8Kay+Ew9+Ayz4LV35Jc7tOnTrFli1bNG83VG666Sa6u7u58sor+drXvgYog4377ruPJ554IuT23G43e/fupbq6WvVcv9vt5p577iEjI4PPfOYzKx77ox/9iM7OTr74xS+ueJy/v68Q4riUcn8we7QYuXcBaxY8L5t7TTd8owY7ENNiPsymw+X28NTJHm65YC2fffcmTveN8/rZFWZKX/hnkF748//L8fQrYNdfwcv/CuO9xhmtA76R7FKcTievvPIKDz74IA8//PCi9xwOB9deey2bNm3i9ttvx+v14vF4uPXWW9m+fTs7duzgu9/97rI2Dx8+zN69e+cD+w9+8AMOHDjArl27uOGGG+ZH6rfeeiu33347Bw8e5HOf+xyg/FK86KKL2LhxIz/4wQ/82nzVVVfxi1/8Iuy/hRq0uP38OPBpIcTDwEFgTErZo0G7Adm3b5+ezRtKTIv5MJuOE51jTM96uWxjIZdvKuTep0/x0zfOcVFl/vKDRzug5udw0d9Czlr27VsLazOg9udQ/xu48FO62fm139fT0O3QtM2tpVl85bptAKSnp/s95ne/+x2HDh2iqqqK/Px8jh8/Pu/DN998k4aGBtatW8ehQ4f49a9/TUVFBV1dXZw8qSQA+ss3P3LkyKJ+8Bd/8RfcdtttAHzpS1/iwQcf5I477gCUm9avvvoq8fHxfPWrX+XEiRO8/vrrTExMsGfPHq699lpKS0sXtV9WVobL5WJoaIj8fD9+1AA1qZC/AF4DNgkhOoUQHxdC3C6EuH3ukKeAs0Az8APgb3SxdAFmG1lFQkyL+TCbjtdbhhACLlyfR0piPJdVFVLTHuAm6ZlnAAl7PwLMaSncBCU7oe6XxhmtA4FG7r/4xS+4+eabAbj55psXjYgvuOAC1q9fT3x8PLfccguvvPIK69ev5+zZs9xxxx0888wzZGVlLWuzp6eHwsLzhRdPnjzJpZdeyo4dO/jZz35GfX39/Hsf+MAHiI+Pn3/+3ve+l9TUVAoKCrjiiit48803/WopKiqiu7s79D+ESoKO3KWUtwR5XwJ/q5lFKjDbyCoSYlrMh9l0vHZ2iM0lWeSkJQFQVZzJ72q6cbrcZCQvuYRPPw15lVCwEVigZccH4Pkvw1AL5OtTWdE3wtYLfyP34eFhDh8+TF1dHUIIPB4PQgi+853vAMvTCYUQ5ObmUltby7PPPsv999/Po48+yv/+7/8uOi41NXVRjvmtt97Kb3/7W3bt2sVDDz00nzzgzy5/5/SnZXp6mtTUVHXiw8CSxRp8GQR2IKbFfJhJh8vt4fi5ES5af/6n+8aiDACa+saXHDwObS/DpmtgLqDMa9lxIyCg/tdGmK0L/optPfbYY3z4wx/m3LlztLW10dHRQUVFBS+//DKgTMu0trbi9Xp55JFHuOSSSxgcHMTr9XLDDTdwzz33zGcRLWTLli00NzfPPx8fH2fVqlXMzs7ys5/9bEU7f/e73zE9Pc3Q0BAvvvgiBw4cWHbMxMQEvb296LmnhSWDe1VVVbRN0IyYFvNhJh21HWO43F4uXH++ROymkkwAziwN7i0vgGcGqs6vOZzXklUKxduh7YjuNuuFv9zwX/ziF7z//e9f9NoNN9wwPzVz4MABPv3pT7NlyxYqKip4//vfT1dXF5dffjm7d+/mQx/6EPfee++ydq+55hpeeuml+eff+MY3OHjwIBdffDGbN29e0c6dO3dyxRVXcOGFF/LlL395fr599+7d88fU19dz4YUX6rvqVkoZlX/79u2T4XLmzJmwP2s2YlrMh5l0/O8rZ+W6u56Q/Y7p+dc8Hq/c9KWn5Nd/X7/44N99Wsp710jpnpl/aZGWJ/5Rym+WSume1cy+hoYGzdoKxtTUlGHnklLK973vfbr1hb/5m7+Rf/jDH4Ie5+/vCxyTKmKsJUfuxcXF0TZBM2JazIeZdJzpc5KblkhBRtL8a3Fxgo1FmctH7h1HYc1BiD9fRXCRlrUXwowT+swz7RQKRlfr/Na3vkVPjz6Jfzt37uSqq67SpW0flgzudin7CTEtZsRMOpr7x9lYlLnsptzG4ozFwX3aAQONsHrx2pZFWtZepDy2v66XubridrsNPd+mTZu47LLLdGn71ltv1aXdhVgyuMdqTJgTu2gxiw4pJWf6nGwozlj2XlVxJn0OF2NTs8oL3W8BElYvzvRZpCV7NeSshfbXdLRaP2KbdYR4Dt3PECNGjLAYcCrBu6poeXDfVKzcVJ3PmOmay81fHaTi9tqL4NxroKLsSAxrY8ngbsSuL0YR02I+zKKjuc8JwMa5QL6QtflpAHSOzBWS6joOeeshbfHGy8u0lB2AiX6l/ozF8HpVVsO0AEZosWRwz8nJibYJmhHTYj7MoqOpfy64+xm5l2Qp0y09Y9PKKLzz2LL5dvCjpWSn8thXv+xYs/N22qxDCywZ3Pv6+qJtgmbEtJgPs+g40zdOdmoihZnLt2NLT04gKyWB3rEpcHSDs3fZfDv40VI0V2Gwz+/GaqZmdnbW7+vx8fHs3r17/t+3vvUtzc994403cvbs2ZA+c/nllxOo8q1Py8DAAIcOrbQXUvhY8qtw7dq10TZBM2JazIdZdDT1O9lYlBGwtO+q7FRl5O4bha/aueyYZVpSsiBnnSWDe1JSkt/XU1NTqamp0e289fX1eDwe1q9fr/ozHo9nxfd9WgoLC1m1ahVHjhzh4osvjsjOpVhy5H7mzJlom6AZMS3mwyw6zg44qSxcPiXjoyQ7hV7HNAycUl4oXL5y0q+Wkh2WnJYJ9V5IeXk5d999N7t372b//v1UV1dz9dVXU1lZyf333w8oBcIuu+wydu/ezfbt2+fLFizkZz/7Ge9973vnn3/qU59i//79bNu2ja985SuLznfXXXexd+9efvlLpUjbT37yk/m2FxYQW6jlfe97X9CSBuFgyZH7jh07om2CZsS0mA8z6HC63Aw6Zygv8F/mFmBVdgoNPQ7ob4SMkmU3UyGAluJtcPopmJ2CRA0LVz39eejVeIFUyQ64RplmSUtL83vI1NTUoqX9d999NzfddBOg/HKpqanhH/7hH7j11ls5cuQI09PTbN++ndtvv52f//znXH311Xzxi1/E4/H4rV9z5MgRbrnlfP3Eb37zm+Tl5eHxeLjqqqs4ceIEO3cqv5ry8/Pna9Xcf//9TE5OUlNTw0svvcTHPvax+TLDC7Xs37+fL31J+41ULDlyN1tJ1kiIaTEfZtBxbkgpb1ue7z+ggTJyH3S68PY3QJH/eid+tRRvUzbz6D+lia1GEajkr29axvfPF9gBrr/+ekD5kjt48CCZmZkUFhaSnJzM6OgoBw4c4Ic//CFf/epXqaurIzNzeWbS0vK/jz76KHv37mXPnj3U19fT0NAw/97CcwPzXwqXXXYZDodjflHZQi16lf615MjdbCVZIyGmxXyYQce5IWUEuS5/5ZE70gsDp2H/R/0e41dL8dxumX31wfPiQ+Ea7W9kLiTQZh0rkZys3IyOi4ub/7/vudvt5rLLLuOll17iySef5NZbb+XOO+/kr//6rxe1sbD8b2trK/fddx9Hjx4lNzeXW2+9ddEUi9ryvwuP06v0b2zkHmViWsyHGXS0zY3c160wci/OSqFMDBDnnvI73w4BtORWQGK65ebdA43cI+HcuXMUFxdz22238YlPfCJo+V+Hw0F6ejrZ2dn09fXx9NNPr9j+I488AsArr7xCdnY22dnZwGItZ86cYfv2ZdtTR0xs5B5lYlrMhxl0nBucpCAjmfSlm3EsYFV2KpvE3GKkoq1+j/GrJS4OCjbAoDluHKsl0Mh96Zz7oUOHVKdDvvjii3znO98hMTGRjIwMfvzjHy875tprr+XFF1/kne98J7t27WLPnj1s3ryZNWvWBM1wSUlJYc+ePczOzs5vCHLs2DHuv/9+HnjgAQBeeOEFrr32WlX2hoSa0pF6/Iuk5G9NTU3YnzUbMS3mwww6/vL+V+UN/3lkxWNGJ2fkv3zhNim/kiXl1JjfYwJq+eVHpfzu9kjNNLTk78TEhGHnWsjk5KQ8ePCgdLvdmrW5UMull14qh4eH/R73tiv5u22bvtt5GUlMi/kwg45zQ5MrzrcDZKUksCWhi7HEYiV/3Q8BteRvVDbTnjVHqQU16LklXbDzfu1rX6Orq0vTNkFZxHTnnXeSm5urWds+LBncF25/ZXViWsxHtHVMz3rodUyvmCkDys25jfG9dCeWBTwmoJaCjYCE4dBWXUaTaNb8ufrqqzVd3ObTUlhYyPve9z7N2l2IJYN7WVngzmw1YlrMR7R1tA/PZcqskOMOgJSskT20eksCHhJQS/4G5XGoKRwTl5hhTIXJQCtUrYgaLZH+XS0Z3AcHB6NtgmbEtJiPaOtoGwye4w7A5BDpcoImT+CdowJq8QX3wciCe0pKCkNDQ4YEeKM369CTYFqklAwNDUW0t4Als2UyMgIvybYaMS3mI9o65nPc84KM3IdaADg1XRjwkIBakjMgc9V8G+FSVlZGZ2cnAwMDEbWjBo/HQ3x8vO7nMQI1WlJSUiL6FWnJ4B6oOpwViWkxH9HW0TY0QU5aItlpQfYMHVYC82l3EZMzbtKSll/OK2rJ3xDxtExiYiIVFRURtaGWnp4eVq1aZci59MYILZaclokV7TcndtESbR1qMmUAGGrBK+LpkIUMOWf8HrKiloKNyrSMRXZlirZftCS2WUcAAhUQsiIxLeYj2jrahiaCz7cDDLcwnb4aNwkMTfgP7itqyd8I06MwORSmpcYSbb9oiRFaLBnch4eHo22CZsS0mI9o6phxe+kenVI9cnfnKFMiQ06X30NW1JI3V598pC1EK6ODXfoXGKPFksG9tLQ02iZoRkyL+Yimjs6RSbwS1uUFGdlJJUddzGW9BJqWWVFLbrnyONwahqXGY5f+BcZosWRwb221RmdUQ0yL+YimDl+mTHlBkODu7IcZJ8lFc8E9wLTMilpy1ymPFhm526V/gTFaLBncN2/2XwHPisS0mI9o6jhfDTLItMzcytKkoipSE+MDTsusqCUxVUmHtEhwt0v/AmO0qAruQohDQojTQohmIcTn/by/VgjxghDiLSHECSHEe7Q39Tx67pdoNDEt5iOaOs4NTZKRnEB+epAVjKPnlMfccvIzkgKO3INqyS2HEWuMiO3Sv8AYLUGDuxAiHvg+cA2wFbhFCLG0vuiXgEellHuAm4H/1NrQhezdq+EGA1EmpsV8RFNH29AE6/LTAm6KPc9IGyAgZw35GckMBhi5B9WSW26Zkbtd+hcYo0XNyP0CoFlKeVZKOQM8DLx3yTES8JWlywa03zNqAWbYTEErYlrMRzR1tA9NUq4mU2bknDKlkpBMQXpSwBuqQbXkVoCj2xLVIe3Sv8AYLWqC+2qgY8HzzrnXFvJV4ENCiE7gKeAOTawLgBk2U9CKmBbzES0dsx4v7cOTK+6+NM/oufkbonnpSQwHmJYJqiW3HJAw2h6asVHALv0LjNGi1Q3VW4CHpJRlwHuAnwghlrUthPikEOKYEOJYT08Pg4OD9PT00NXVxcjICC0tLUxNTdHQ0IDX653f8sr3LVddXY3X6+Xw4cNMTU3R0tLCyMgIXV1d+Npra2vD6XTS2NiI2+2mtrZ2URu+x7q6OlwuF01NTTgcDtrb2+nv76e/v5/29nYcDgdNTU24XC7q6ur8tlFbW4vb7aaxsRGn00lbW1vImp5//nm8Xi8NDQ2W1/T8888v8pNVNR0/fjxg39NTU01TB26vpDQzPqimmf4mvNlraWxsJCtZMOh0MTAwsEzTm2++ufL1lKOUsXX1nTa9n1577TXVMcLsfe/IkSMhxz2fJtUE280DuAh4dsHzu4G7lxxTD6xZ8PwsULRSu5HsxOTxeML+rNmIaTEf0dLx7Mkeue6uJ+Rb7SMrHzg7LeVXsqU8/M9SSil/8FKLXHfXE3J0cmbZoUG1OHqVnZxe/+9wzTYMu/QvKSPTgoY7MR0FNgohKoQQSSg3TB9fckw7cBWAEGILkALoViausbFRr6YNJ6bFfERLR8uAkga5vjDInPtYJyDnp2XyM5TMGn/pkEG1ZBRBYpolbqrapX+BMVqCBncppRv4NPAscAolK6ZeCPF1IcT1c4f9I3CbEKIW+AVw69w3jC4YVYXOCGJazEe0dLQMOCnKTCYrJUg1SF8gzpkL7unJgP+FTEG1CKG0Y4Hgbpf+BcZoUTXnLqV8SkpZJaWslFJ+c+61f5JSPj73/wYp5cVSyl1Syt1Syuf0NLq7W9dkHEOJaTEf0dLRMuCkslBFLfkFOe6g3FAF/yN3VVpy1sBYR/Djooxd+hcYo8WSK1Tz8vKibYJmxLSYj2jokFJydmCCyiI1aZBtEJ+kpEJyPriPTC6v3a5KS7Y1grtd+hcYo8WSwX1ycjLaJmhGTIv5iIaOoYkZxqZm1Y3cR84pATlOuXxz03zBffm0jCot2WUwNQIuZ0g2G41d+hcYo8WSwT0uzpJm+yWmxXxEQ0dLvxJY1U3LtMNcCiNASmIcSQlxjPkZuavS4mvL5KN3u/QvMEaLJf9aiYlBbjhZiJgW8xENHaozZUDJlslZM/9UCEFuWqLfkbsqLdlzbY2aO7jbpX+BMVosGdydTnP/fAyFmBbzEQ0dDT1jZCYnUJqduvKBs9Mw0X8+IM+Rm5bEqJ+Ruyotvi+KMXOvUrVL/wJjtFgyuBcUFETbBM2IaTEf0dBR1+Vg2+os4uKCFAxzdCmP2WWLXs5OTfQb3FVpySiBuETTj9zt0r/AGC2WDO6dnZ3RNkEzYlrMh9E6Zj1eTvU42F6aHfzgsTnblgT33LQkRqeWT8uo0hIXB9mrz7dtUuzSv8AYLZYM7hs2bIi2CZoR02I+jNbR3O9kxu1lR1kowX3xtExOWqLfVEjVWiyQDmmX/gXGaLFkcK+vr4+2CZoR02I+jNZxsmsMgG2qR+4CshbvwZmTlsTo5AxLF4ar1pKz1vTTMnbpX2CMFksG9127dkXbBM2IaTEfRus42TVGelI86wvUZMp0QEYxJCQvejknLZFZj2RyxrPoddVasstgvAfc/ksHmwG79C8wRoslg3usaL85sYsWo3Wc7HawtVTFzVRQgvuS+XaA3DQltW5pOqRqLdlrAHn+hq0JsUv/AvNs1mE6YkX7zYldtBipY8btpaHbwfbVKqZkQJmW8RPcc+ZWqS7NmFGtZT4d0rxTM3bpX2CtzToMJfYNbk7sosVIHcfahpma9fCOShWpcVIGDu6pysh9aXAPbeSOqefd7dK/IDZyD0jsG9yc2EWLkTpeON1PUnwc76jMD37w5BC4p5dlygDkzhUPW5oOqVqL7wsjNnI3hNjIPQC+La3sQEyL+TBSx4unB7igIo/05ITgB/sC7woj96XpkKq1JCQri5lMHNzt0r/AGC2WDO5VVVXRNkEztNDyyNF2vvp4PV2jIeyvqAN28YtROjpHJmnqd3L5pkJ1HwiwgAkge+6G6uiSDTtC0pKzxtTTMnbpX2CMFksG9/Z2c9fACIVItEgp+eaTDdz1qzoeerWNK77zIi+e7tfQutCwi1+M0nG4UfHV5ZuK1H3AF9wXVIT0kZwQT1pSPKNTi0fuIWkx+UImu/QvMEaLJYN7cXFxtE3QjEi0HG7s5wcvt/LhC9fx8ueuYE1eKl9/ogG3x6uhheqxi1+M0OH2eHnwlVa2rsqiUk0lSFCCe2IapOb6fTs3LWlZKmRIWnLWKOfwRqf/BMMu/QuM0WLJ4D46OhptEzQjEi3/89JZSrNT+KfrtrImL43PHdrM2YEJHjsenRocdvGLETp+/VYX54Ym+Yd3VSGEivx2OJ/jHuD4nLTEZTXdQ9KSvQY8M0rVSRNil/4FxmixZHBPSUmJtgmaEa6WE52jvNE6zMcuqSAxXnHju7cWs2dtDv/2xya8Xt32Jw+IXfyit47OkUm+94cmdpZl884tKqdkQJkP9zPf7iPHT033kLSYPB3SLv0LjNGi4hZ9DDPy0JE2MpMTuOnA4k0bPnJROX//SA1vdYyyb53/n+8xjMMxPcszdb20DDhJSYzH6XLzq+pOPB7Jv92yW/2oHZQpk5LtAd/OSUuiZ8wRvrEL67qvORB+OzFMgSWD+/T0dLRN0IxwtLg9Xv7Y2M+7t5WQmbJ4R5crNheRGC94rr7X8OBuF79opaOmY5QPPfAGTpebpPg4ZjxeUhPj2VqaxX0f2EWFmloyPgJs0rGQHD813UPS4mvbpKV/7dK/wBgtlgzuOTk50TZBM8LRUtMxytjULFduXv6TPjs1kYsqC3imvpfPX7M5tJFhhNjFL1ro6Bmb4rYfHyM3PZGffPwCdq/JwSshXk39GH8E2KRjIblzlSG9XjlfpyYkLSlZkJJt2mkZu/QvMEaLJefc+/r6om2CZoSj5XBjP/Fxgkur/C9ZP7SthHNDk5zuG4/UvJCwi1+00PHF35xk0uXmwY8cYM/aXIQQ4Qd2WDHH3UdOWiJeCeMu9/xrIWvJXmvadEi79C8wRoslg/vatcvzfK1KOFpeOD3A/nW5ZKX432T3XVuVNKs/njI268EufolUx9kBJ4cb+/nkZZVUFWdqY1SATToWcr542PmbqiFryS4z7bSMXfoXGKPFksH9zJkz0TZBM0LV0js2zakeB1f4mZLxUZiZTFVxBkfbhiM1LyTs4pdIdfz4tXMkxgtuORg4EIdMgE06FuIr+7tw3j1kLdllph2526V/gTFaLBncd+zYEW0TNCNULW/OBexLNqxcRXB/eR7Hz43gMTAl0i5+iUSH0+XmseOd/PnOUooyNUx3C7BJx0Jy/NR0D1lLdhlMj8F0BFk3OmGX/gXGaLFkcH87l/6sPjdCamI8m0tW/rl/oDyX8Wk3Zwycd7eLXyLR8eLpfpwuN7dcoPHP7gCbdCzENy0ztqAEQchafOcw4aYddulfECv5G5C3c+nPt9pH2FmWTUL8yq7bvy4PUOqFG4Vd/BKJjpfODJCVksDetRpnQwSo476Q+cqQC4qHhazFxOmQdulfECv5G5C36zf49KyH+m4He9YGz18vy02lJCuFo20jkZgXEnbxS7g6pJS8dGaQSzcWBv3yDbFhVcE920/Z35C1mHhHJrv0L4iN3APydv0GP9k1htsrVY0KhRDsL8+NjdzDIFwdZ/qc9DqmuSxAimrYrLBJx0IS4uPITElYNC0TspaMYohLiI3cdcY0I3chxCEhxGkhRLMQ4vMBjvlLIUSDEKJeCPFzbc1cTG1trZ7NG0ooWt5qV4oN7Vb5k3/P2ly6x6YZGHeFZVuo2MUv4ep46cwAAJdVqazPrpYVNulYytLKkCFriYtXMnJMGNzt0r/AGC1Bg7sQIh74PnANsBW4RQixdckxG4G7gYullNuAv9fB1nm2bdumZ/OGEoqWtzpGKMtNVZ2Fsa00C4D67rGwbAsVu/glXB0vNw+ysSiDVdnmq52WAAAgAElEQVSp2ho0X8c9eGplbtriEgRhacleY8rgbpf+BcZoUTNyvwBollKelVLOAA8D711yzG3A96WUIwBSSl1XzzQ3N+vZvKGEoqWua4xdZepv1G2dD+7GpLXZxS/h6JBSUtM+wv7yPO0NUrGAyUf2XAkCH2H5xKS57nbpX2CMFjXBfTWw0NOdc68tpAqoEkIcEUK8LoQ45K8hIcQnhRDHhBDHenp6GBwcpKenh66uLkZGRmhpaWFqaoqGhga8Xi/V1dXA+ZsP1dXVeL1epqenmZqaoqWlhZGREbq6uvC119bWhtPppLGxEbfbPf/zx9eG77Gurg6Xy0VTUxMOh4P29nb6+/vp7++nvb0dh8NBU1MTLpdrfr/DpW3U1tbidrtpbGzE6XTS1tYWsiaHw4HX66WhoWFFTb1Do3QMT7GpOF21ptH+HlZnJ3P8rDGaxsbGFvkpmCaz+qm0tDRg3wuk6a3mLhzTbjYXpmivaawTT3wKpOYG1ZQi3AyOT81rKigoCHo9LdU0HpeNdHQz2N9nKj9lZmaqjhFm73vJyckhxz2fJtVIKVf8B9wIPLDg+YeB/1hyzBPAb4BEoALlyyBnpXb37dsnw6W1tTXsz5oNtVrebB2S6+56Qv6hoTek9j/102Pysm8fDsOy0LGLX8LR8ZvqTrnurifkqZ4x7Q16+ENS/vt+VYf+02/r5I6vPDP/PCyfHH1Qyq9kSTnWFfpndcQu/UvKyLQAx2SQuC2lVDVy7wIW/h4sm3ttIZ3A41LKWSllK3AG2Kj+KyY0MjIy9GracNRqOdWjTK34plrUsq00m3NDkzimZ4MfHCF28Us4Omo6RklLimdjkUa1ZBaiIg3SR05aEo5p9/xWi2H5xKS57nbpX2CMFjXB/SiwUQhRIYRIAm4GHl9yzG+BywGEEAUo0zRnNbRzEbOz+gcqo1Cr5VSPg5y0REqyQlvS7vsyaDBg3t0ufglHR23nKNtLsyOr/BiIkIK7kuvumFYqQ4blE9+5TDbvbpf+BcZoCRrcpZRu4NPAs8Ap4FEpZb0Q4utCiOvnDnsWGBJCNAAvAJ+VUg7pZbTXpBv4hoNaLQ0942wpyQq5Pvs2A2+q2sUvoeqY9Xip73awa0229sao2KRjIblzJQh86ZBh+SRr7paayUbudulfYIwWVZt1SCmfAp5a8to/Lfi/BO6c+6c7aWlpRpzGENRo8Xglp3sd/NUF60JuvygzhYKMJE736h/c7eKXUHWc7h1nxu1l1xodNmBQsUnHQnKWVIYMyye+TTtMFtzt0r/AGC2WXKE6PGxsKVs9UaOlbWiC6VkvW1aFN5+7sSiTM33OsD4bCnbxS6g6fFNe20p1GLmHkAYJy2u6h+0TE+a626V/gTFaLBncS0sD17S2Gmq0nO5VKjtuWRXazVQfm0oyaeob92U26YZd/BKqjsbecVIS41ibp8NoTMUOTAvJTVtcXyZsn5gw190u/QuM0WLJ4N7a2hptEzRDjZamPidCQGVheHfYNxZnMDHjoWs0hBzZMLCLX0LVcaZvnKriTP1upgbZpGMhOamLR+5h+8SEOzLZpX+BMVosGdw3b94cbRM0Q42Wpv5x1uSmkZoUH9Y5fFu9Nek8NWMXv4Sqo7F3nE1abae3FBWbdCwkMyWBOHG+pnvYPskug6kRcOk/nacWu/QvMEaLJYN7TU1NtE3QDDVamvudbCwKPy+2ai73Wu+NO+zil1B0DDldDDpdbAqyeUrYqNikYyFxcYLs1MT5bJmwfWLCXHe79C8wRoslg/vevXujbYJmBNPi9ng5OzDBhuLwg3t2WiLFWcmc1jm428Uvoejw/U31C+7qc9x95KYlzWfLhO2T+Vx38wR3u/QvMEaLJYP726lof/vwJDMeb8QrH6uKM3WflrGLX0LR4bvZrcu0jMpNOpaSvaAyZNg+MeFCJrv0L4ht1hGQt1PR/qZ+JSBHMi2jfD6Tpv5xvDpumG0Xv4Si43TvOLlpiRRmqpsTDwmVm3QsJTctidEpZVombJ9klICIN9XI3S79C0y0WYfZ8FVNswPBtDTPBffKSIN7cQbTs15dM2bs4pdQdJyey5QJdeWwKnyjZhV13BeSk3p+5B62T+ITTLdph136FxijxZLBfffu3dE2QTOCaWnqG2d1TioZyaoWEwdkfUE6AGcHJyJqZyXs4he1OqSUNPc757ORNCfEHHcfC6dlIvKJydIh7dK/wBgtlgzujY2N0TZBM4Jpaep3Rjxqh/Mj/5Z+/ebd7eIXtToGxl2MT7upLEzXx5AQV6f6yE1LwulyM+vxRuYTky1kskv/AmO0WDK4V1RURNsEzVhJi5SS1sEJTYJHfnoSWSkJnB3UL7jbxS9qdTQPKH/LDXqU+QUluCemQWpuSB9bWF8mIp9kl4GjG7ye8NvQELv0LzBGiyWDe3d3d7RN0IyVtPQ5XEzOeFgf5srUhQghWF+YwdkB/aZl7OIXtTpa5u+H6DRyH21XAmyI8/m++jJjUzOR+SS7DLyz4NR110zV2KV/gTFaLBnc8/J02KcySqykxTfK9s2XR0plYQYtA/qN3O3iF7U6WgYmyEhOCLnGvmrCSIME5YYqKPVlIvKJyRYy2aV/gTFaLBncJycno22CZqykxTfKrtAouK8vTKfP4cLpcmvS3lLs4he1Opr7nVQWpuuTKQMhr071kTtfGXI2Mp+YLNfdLv0LjNFiyeAeF2dJs/2ykpbWwQlSE+M1Gxn65u5bdZqasYtf1OpQgrtO26XNTMLEAOSEXsM/Z74y5ExkPjHZKlW79C8wRosl/1qJiYnRNkEzVtJydsBJeUE6cRpVG/QFIr2mZuziFzU6nC43vY5pTTKZ/DKf4x56cM+eC+5jk7OR+SQlG5KzTBPc7dK/wBgtlgzuTqd5KtVFykpaWgcnWK9hmt3a/DTihPKloQd28YsaHfM3U/UauY+2K485a0P+aGZyAvFxgtGpmch9YqJNO+zSv8AYLZYM7gUFBdE2QTMCaZlxe+kYmdLsZipAckI8a/LSaNFpIZNd/KJGh2/l8Aa9Ru6j55THMIK7EIKc1ERGJmcj94mJct3t0r/AGC2WDO6dneYYSWhBIC3tw5N4vFLTkTsomTd6LWSyi1/U6GgecJIQJ1iXr9NemKPtEJ+k1HIPg5y0RMYmZyP3iYlWqdqlf4ExWiwZ3Dds2BBtEzQjkBbf1ElFgbYjw8rCDFoHJ3QpIGYXv6jR0dyv3A9JjNfpEhptV6ZEwrzxlpOWxMjkTOQ+yS6DqWGY0W99hFrs0r/AGC2WDO719fXRNkEzAmlpHdQ2DdLH+sIMXG59CojZxS9qdLT0O9mg13w7KME9jCkZH7lz9WUi9sl8rntXZO1ogF36FxijxZLBfdeuXdE2QTMCaWkdnKAgI4nsVG3vqvumefQoIGYXvwTTMeP2cm54Ur/5dog4uGenJjE6ORO5T+bTIdsja0cD7NK/wBgtlgzub4ei/WcHJjQftcP57A49Mmbs4pdgOtqGJvB4pX7BfT7HPbKR+8jkbOQ+8dkwci6ydjTALv0LYpt1BOTtULT/7OAE6zWebwcoyEgiMyVBl1x3u/glmA7dM2UiyHH3kZuexNSsh207Iywtm7lKubE7Gv3gbpf+BbHNOgJi929wx/Qsg04XFTqUktWzgJhd/BJMhy+4a53JNE8EOe4+8tKVEgQvvR6hT+LiFDtG2iJrRwPs0r8gNnIPiN2/wX3lAbTMcV9IZWG6LsHdLn5RM3JfnZNKWlJkG6gEJIIcdx+++jKrKzdFbk/OOlNMy9ilf0Fs5B6Qurq6aJugGf60zFeD1GlkWFmYQa9jWvMCYnbxSzAdzRptoBKQCHPc4fzIvfrkmcjtyV1nimkZu/QvMEaLJYN7VVVVtE3QDH9aWgcmiBOwNk+f4O77RaB1ATG7+GUlHV6v5OygAWmQEeS4A+SlK1lW6XlFkduTWw5TIzA9FnlbEWCX/gXGaLFkcG9vj35allb409IyOMGavDSSEvRxj2/zD613ZbKLX1bS0TU6xfSs19RpkHB+Wqalqy9ye3w3dqM8NWOX/gXGaFEVPYQQh4QQp4UQzUKIz69w3A1CCCmE2K+dicspLg7/56rZ8KeldWBCt/l2gHVzBcRaNB6528UvK+nQPVMGNAnu2amJCAGe+NTI7cmdC+5RnpqxS/8CY7QEDe5CiHjg+8A1wFbgFiHEVj/HZQJ/B7yhtZFLGR0d1fsUhrFUi2/fVK3LDiwkJTGestw0zdMh7eKXlXToHtw1yHEHSIiPIzs1kb5RDXycW648Rjljxi79C4zRombkfgHQLKU8K6WcAR4G3uvnuG8A/wJMa2ifX1JSdNrWLAos1dLrmGZq1qNfmt0c63XImLGLX1bS0dzvJC89af6GpeZokOPuIy8tCedsxM0oG3QnZ0d9WsYu/QuM0aImuK8GFtb87Jx7bR4hxF5gjZTySQ1te1uidxqkj/UFGbQOOnUpIGZnmgcMuJkKEY/cQVnINDqlUUZUrjly3WOoJ+I7dkKIOOD/Av+o4thPCiGOCSGO9fT0MDg4SE9PD11dXYyMjNDS0sLU1BQNDQ14vV6qq6uB8wn/1dXVeL1ezp49y9TUFC0tLYyMjNDV1YWvvba2NpxOJ42Njbjdbmpraxe14Xusq6vD5XLR1NSEw+Ggvb2d/v5++vv7aW9vx+Fw0NTUhMvlmk9bWtpGbW0tbrebxsZGnE4nbW1tIWs6ffo0Xq+XhoYGpqameONUGwCpnnFdNaXMjDI96+W5I0c109TY2LjITz5NVvPT1NSU377n8Xg43TNGeV6KbpqGz9YAcHbEE7GmZGYZdrqCXk9q/DSekI8caYuqn4aGhlTHCLP3vb6+vpDjnk+TaqSUK/4DLgKeXfD8buDuBc+zgUGgbe7fNNAN7F+p3X379slwGRsbC/uzZmOplq89Xi83f+lp6fV6dT3vq82Dct1dT8g/ne7XrE27+CWQjn7HtFx31xPygZfP6nfy574s5dcLpPR4Im7qs7+skQfueU4Do6SUz3xBym8UaWJXuNilf0kZmRbgmAwSt6WUqkbuR4GNQogKIUQScDPw+IIvhzEpZYGUslxKWQ68DlwvpTym/ismNPr6NEjvMglLtZwddFJRkI4Q2uybGgjfZtlaFhCzi18C6TAsUybCHHcfuelJjEzO+gZhETZWDu5pcEbPx3bpX2CMlqA9SErpBj4NPAucAh6VUtYLIb4uhLhebwP9sXZt5PORZmGpFq33TQ1EYWYyGckJmpb+tYtfAulo7h8HzJ8G6SMvLYlZj2RyxhN5Y76MmSimQ9qlf4ExWlQND6SUT0kpq6SUlVLKb8699k9Sysf9HHu5nqN2gDNnNFhSbRIWanG5PXQMT+p+MxV8BcS0zZixi18C6ajvdpCdmkhpto6ZDqPtkLNGk6Zy5zJ6hidmIm9sfiFTW+RthYld+hcYo8WSK1R37NgRbRM0Y6GWjuFJvPL8ClK9WV+Qrum0jF38EkhHfbeDbaVZ+k2ZucaVHPfcCk2ay5tbpToyqUVwj35dd7v0LzBGiyWDu11Lf/pG0Xps0uGPysIMusemmZzRJl3OLn7xp2PW4+V07zjbSrP0O/Fwq/KYt16T5jQduSemKLXdozgtY5f+BcZo0almqb7YtfSnrxyAHnXc/TFfY2Zggu2rsyNuTyu/HD83zGPHu2gfnuCKTUV86MJ1pCTGa9K2GvzpaO53MuPxsq008r9TQEZ8wV2jkXu6hiN3mCv926ZNW2Fg1+teL2Ij9yizUEvLgJPirGSyUrTdNzUQWu+nqoVfHj3awU3//TpPnOim3+HinidPcc33XmbI6dLAQnX401Hf7QDQeeR+VnnUeFpmyKlRcM8tj+q0jF2ve72wZHC36zd4c79zfo9TI1BSLrVLh4zUL3881cfnfnWCiyrzOfL5K3n+zj/jRx+7gO7RKT75k+NMz2qQ9aECfzrqu8dISYzT937IcCukFUCKNl8gmSkJJMQJhrSYlgGlgJijC9watRcidr3u9cKSwd23oswO+LRIKWkZMDa4pyTGszonVbPqkJH4ZXRyhs//uo7NJZk8+JED879e/qyqkO/etJvj50b47vPGZEv401Hf7WDLqizi43RcfzB8VrP5doC4OEFWstDuV0/OOkCer39jMHa87vXEksF927Zt0TZBM3xaBsZdjE+79c2h9oOyn6o2I/dI/HLvU42MTMxw3wd2Latj/54dq7hxXxk/PNJGx/BkpGYGZakOt8dLQ7eD7XrOt4Myn63RfLuPkpx0BrWcloHz9wYMxo7XvZ5YMrg3NzdH2wTN8GlpnguwRo7cQUmHbB2c0GQVY7h+aR+a5LHqTj580bqAN3Y/8+5NxMcJvvVMYyQmqmKpjlM94zhdbvaX5+p30tlpGOvUdOQOkCrcDGo1cvd98QxHJ7jb8brXE0sG97KysmiboBk+LS1GLG33Q2VhOpMzHnrGIq/UHK5f/utPzcTHCW7/s8qAx5Rkp/CJSyt48kSP5nXol7JUxxutQwAcrMjX76Sj5wCp2c1UH6vzsxgc1yi4Z66CxLTzN34Nxo7XvZ5YMrgPDg5G2wTN8GlpGZggIzmB4qxkQ8/v2+jZVzclEsLxS8/YFI8d7+Sm/Wsozlp55edH3lFOUnwcP3q1LUwL1bFUx+tnh1mXn0aJnitTNc5x95EW52ZwYkab+jJCQF4lDEVnBG3H615PLBncMzKMHd3qiU+Lkimjf8GwpWwsygSgSYPgHo5fHjnawaxHctulwYNaQUYy1+0q5bHjnYxNabELhX8W6vB6JUfbhjlYkafb+YDzo2Ed5txn3F7GXRrVdc+PXnC343WvJ5YM7rOz+l3YRuPTYnSmjI+CjCRy0xLni2JFQqh+8Xgljx7t4NKNBazNT1P1mY9eXM7kjIfHjneGY6IqFuo43TfO2NSsvlMyAENNkJINadqeJytZGSxoNjWTX6nkunuMvwbteN3riSWDu9frjbYJmuH1enG63PSMTc9PkRiJEIKNRZmc6Yt85B6qX15qGqB7bJqbD6ivkLd9dTa7yrL5lY7BfaGO11qU+fYL9B65DzZB/kZl6kNDclKVlb2aZczkbwDpicpiJrtd93pjyeCelqZulGcF0tLS5lMRozFyB9hYnEFT33jE87Kh+uXRox3kpyfxrq2h7QT//j2raehx0NjrCOlzalmo46m6HjYVZ7ImT+c+N9QMBRs1b7Y0T5l20yzXPX+D8hiFqRm7Xfd6Y8ngPjw8HG0TNGN4eNiYTSBWYGNRBo5pNwMR/nQPxS+O6Vn+2NjPdbtKl+W1B+O6XaUkxAl+U90Vqomq8OnoHp3i2LkR/nznKl3OM49rHMZ7zgdODYmfURaoaZYOGcXgbrfrXm8sGdxLS0ujbYJmlJaW0jLgJCFOsE7lvLPWVBUro7tIp2ZC8cvz9X3MuL1cvzt0X+ZnJPNnVYX8tqYLjw4bfPt0PFXXA8Cf79K5v/kCpQ4j983ryxACBrSalknLg5QcGG7Rpr0QsNt1rzeWDO6trdFZRKEHra2tNPc7WZefRmJ8dNyxoVj5xdAU4U3VUPzyeG03Zbmp7FmTE9a5rt9dSp/DRXX7SFifXwmfjt/XdrN9dZb+JZgH54J7vvbBvePcOXLTkrQbuYMyeo/CyN1u173eWDK4b968OdomaMbmzZtpGZiI2nw7QGFGMtmpiRGP3NX6ZXhihleaB7luV2nYqZ9Xbi4iKT6OZ0/2hvX5ldi8eTOvtQxR2znG+/cYsHBmqAkQmue4g6KlICNJ26qa+RvOfyEZiN2ue72xZHCvqamJtgmacaz6LdoGJ6KSKeNDCMGm4kxOR3iDUq1fnm/oxeOVXLsj/LnszJRELt6QzzP1vdos0FnAW2+9xbeePsWq7BQ+eNCAfTsHm5SdjhK1XyRVU1NDQUaydtkyAIVVMN4N0/rc0A6Ena57I7RYMrjv3bs32iZoRv66Tbi9kg1RHLkDbF6VyenecbwRzGGr9ctz9X2U5aZGXBv90PYSOkem5muta0VP4ipqO8f4h3dVGbNJiE6ZMqD4RAnuGo7cC+dGnYPG7mlqp+veCC2WDO52Ktr/7GtK6c9ojtwBtqzKYmLGQ+fIVNhtqPHLhMvNy82DvGtrccSrcd+5pZg4Ac/Vazc1M+vxcs/vTrCpOJMb9howJeP1wlCLLvPtoPikICNZu0VMcD64D+hfxG0hdrruY5t1BMBORftlppLjXWnQ1nqB2LJKGUU39IQ/Clbjl5ebBphxe3n31pKwz+MjPyOZA+V5PKNhcP/Fm+30Tni465pN+tZu9zHWAbMTULhJl+b37dtHYWYyEzMeJrQqQZCzDuKTDQ/udrruY5t1BKC6ujraJmjGaw3nWJ2TSqZBW+sFoqo4AyGIaGGQGr88V99HTloiBzQqn3toewln+pya1KR3utx87w9NbCtM4opNRRpYp4L+U8pj0VZdmq+urmbVXMGzXkfklT8BiE9QppEGTmvTnkrsdN0bocWSwX337t3RNkEz+l0JbC7JjLYZpCUlUJGfzqkIRu7B/DLr8fLHxn6u2lxMgkZpn1dvU34BPFvfF3Fb//PSWYYmZrjnxn3GFXAbmAvuOo3cd+/ePV9ts0+Dss7zFG4yfORup+veCC2WDO6NjcZ2Kr2YcXtpGXCyyQTBHZSpmVM94ee6B/PL0bZhxqZmQy43sBKlOansLMuOeGqmf3yaB14+y7U7VpE8oX16ZeATn4Ks1ZAaXr5/MBobG+dLFWtRs3+ews0w2g4z2mzRqAa7XPdgjBZLBveKCm3LokaLlgEnHolpgvvmkkzahydxhjk3G8wvz9X3kZwQx2VVBWG1H4irt5VQ2zFK92j4N4P/84UWZtxePnP1JmP7V3/D+RuUOlBRUUFJlsbTMnD+l8Zgk3ZtBsEu1z0Yo8WSwb27uzvaJmiCb37bdzMz2vjsaAxzamYlv0gpeb6hj0s3FpKWlBBW+4E4tF2Zmgk3a2ZscpZHj3Vw/e5SKgrSjetfXg8MnIGiLbqdoru7m9SkeLJTE+nTMrgXzAV3A6dm7HLdgzFaLBnc8/J0Lr9qEI294yTECf2Xt6tkR5myf+mJzrGwPr+SX+q7HXSNTvHubdpNyfioLMxgY1FG2FMzP3+znckZD5+4RFkhalj/Gm4Fj0u3m6lwXktJVoq20zL5lRCfBH0ntWszCHa57sEYLZYM7pOTk9E2QRNO945TnpsctZoySynOSqE4K5kTnaNhfX4lvzzf0EecgKs265OFcs32Et5sHQ55mf2sx8tDr7ZyyYYCts4tqjKsf/luphbpNy3j01KcnaLtyD0+UfnF0WtccLfLdQ/GaDFHVAmRuDhLmr2M073jVBakRtuMRewsywl75L6SX55r6GN/eR75GfrsEXv19hK8UvkSCYUXGvvpc7i49R3l868Z1r98aZAF+mTKwHktq7JS6NVy5A5QvAN660Dj8g+BsMt1D8ZoseRfKzExujnhWjA8MUPP2DRVxebaF3Ln6mzODk7gmA59G7BAfukYnuRUj4N3a5gls5Stq7JYk5ca8tTMY8c7KchI5vJNhfOvGda/emqVDaeT9esDPi3F2SkMOF3MejTcAahkO0wOgjPyNFQ12OG692GEFlXBXQhxSAhxWgjRLIT4vJ/37xRCNAghTggh/iiEWKe9qedxOiNfsBJt6ruV0XF5lrm+X3fOleA9GcboPZBfnpsbTWuxKjUQQgiu2b6KI82Dqr+YhpwuDjf28/49pYvy7g3rX70nYNUuXU/h01KSlYKURLwhyyJKdiiPBk3N2OG692GElqCRRQgRD3wfuAbYCtwihFh6B+gtYL+UcifwGPBtrQ1dSEGBtql00eBkl5KRcuGm1VG2ZDE7Vys3VWvDCO6B/PJcfS+bSzJVb4IdLldvK2HWI3mhsV/V8Y/XduP2Sm7Yt7iGjCH9a3JYyRNftVPX0/i0aL5KFaB4m/LYV6ddmytgh+vehxFa1AwbLwCapZRnpZQzwMPAexceIKV8QUrpu0PwOqBrxaXOTv02RzaKk91jrMlLZXzImJ+0aslNT2JNXmpYN1X9+WV4YoajbcO8e5t+o3Yfe9bkUJSZzNN16qZmflXdyfbVWWwuWZyKakj/6p0LiDqP3H1adFmlmpoL2WsNG7nb4br3YYQWNcF9NdCx4Hnn3GuB+DjwtL83hBCfFEIcE0Ic6+npYXBwkJ6eHrq6uhgZGaGlpYWpqSkaGhrwer3z9Rd8FdSqq6vxer243W6mpqZoaWlhZGSErq4ufO21tbXhdDppbGzE7XZTW1u7qA3fY11dHS6Xi6amJhwOB+3t7fT399Pf3097ezsOh4OmpiZcLhd1dXV+26itrcXtdtPY2IjT6aStrU29ppZ+tpdmMzExgdfrpaGhwTSa9pTl8EbLAOPj4yFp8v3U9PmpoaGBZ0504pVw0ZpU3TXV15/k6m0lHG7sZXLGvaKfXj/VzskuB++szFymaf369QH7nlZ+cncq/x/PWK9r31u9ejUNDQ0UZSYByipVLTWNppRBb50h11NBQYHqGGGm68mfpszMzJDjnk+TaqSUK/4DbgQeWPD8w8B/BDj2Qygj9+Rg7e7bt0+GS01NTdifNQOjkzNy3V1PyP843GRKLT99vU2uu+sJ2TrgDOlz/rR84kdH5Tvu/aP0er1ambcir7UMynV3PSF/Xd2x4nH3PFEvN3zhSTnkdC17zxCfPPZxKf91q+6n8Wnxer2y6otPyXueqNf2BC/cK+VXc6ScHte2XT+Y8VoJl0i0AMdkkPgqpVQ1cu8C1ix4Xjb32iKEEO8EvghcL6XU8K7Ncnbt0venrN40zG0usX11tim1HKxQFli82RraDu1LtUzNeHi5aUCT2u1qOViRR3l+Gg+/2RHwGLfHy2/e6ubKzUXkpScte98Qn/TU6j7fDue1CCFYk5dG+7DG+dWr987T9ZEAAByKSURBVIH0Qo/+OwuZ8VoJFyO0qAnuR4GNQogKIUQScDPw+MIDhBB7gP9GCezq7mZFgNWL9p/sUm5Wbi/NMqWWysIM8tKTeL11KKTPLdXypzMDTM96dU2BXIoQgg/sX8MbrcMBywA/W9/HoNPFB/at8fu+7j5xOZWaLCX6B/eFWtbmpdE+HH79Hb+Uzu0o1KV/PzbjtRIuptisQ0rpBj4NPAucAh6VUtYLIb4uhLh+7rDvABnAL4UQNUKIxwM0pwlWL9pf0zHK6pxU8jOSTalFCMEF5Xkhj9yXanmyroe89CQOVBi7bPzGfWXExwkeOep/9P6/R1pZm5fGFQFWy+ruk67jgISyA/qeh8Va1ual0TE8qe2es+n5kFsOnce0azMAZrxWwsU0m3VIKZ+SUlZJKSullN+ce+2fpJSPz/3/nVLKYinl7rl/16/cYmRY/Ru8un2EveuUzSrMquWCijw6R6boCqHS4kItkzNu/tDQx3t2lBheXqE4K4VD20r46evnlu0dWtMxyvFzI3z04vKAOy3p7pPON5XHMv0v8KUjd6fLzfCEhptlgzI106X/5hNmvVbCwRQjdzNi5W/w7tEpesam2btWWSxkVi0Xrs8H4EjzoOrPLNTyh1P9TM16uG5nqea2qeHOd1cx7fbyH4eb51+TUvK9P5whMzmBD+z3PyUDBvik46hSciBVm92oVmLpyB3QZ97d0Qnj+tbBN+u1Eg6mGbmbDV/akRWpbh8BYO9a5cI2q5YtqzIpyUpRvSAIFmt5vKab4ixlj9NoUFmYwV/uX8PP3jjH0TZleumnr5/jhdMD/P27qshIDlx2WFefSAmdR2GN/lMysFiLbxGZ9sF9v/Ko87y7Wa+VcDBCiyWDe1VVVbRNCJvqc6MkJ8TN1043qxYhBFduKeKlM8qG1mrwaRlyuvjTmX6u21lKnBGbTAfgzndVUZabxgcfeIPbf3Kcbzx5iss3FfLRBUXC/KGrT4ZaYGoYyi7Q7xwLWKhlTe5ccB/SOLiv2glxidDxhrbtLsGs10o4GKHFksG9vb092iaETXX7CDvLsklKUP70ZtZy1eYiJmY8qm+s+rT85q0uZj2SvzwQeOrDCAozk/nVp97BBeV51HSMct3OUv71A7uCfuHo6hPffPsaY4L7Qi2pSfEUZSZrP3JPTIXVe6HtiLbtLsHM10qoGKFF2y1xDKK42LjUOi2ZnvVQ3z3Gxy4+v8WWmbW8o7KA5IQ4/tjYxyUbg9fCKC4uRkrJw0c72Ls2h6ri6G8fmJeexE8/cTCkz+jqk/bXITlb1zK/C1mqZa0eue4A5ZfAK/9PSfPUqcqlma+VUDFCiyVH7qOj4W0mEW2qz40w65EcXH9+HtrMWlKT4rlkQwHPnuzF4w2ePjc6Okp1+wjN/U5uPrDWAAv1QVeftP5JCYQG1SZfqsWXDqk56y4G6dF1asbM10qoGKHFksE9JSUl2iaExastQ8THiUU3Gc2u5S/2ltE9Ns0rKrJmUlJSePCVVjJTErh25yoDrNMH3Xwycg5G2mD9n+nTvh+Walmbn0aPY5rpWY+2J1pzEOISoO0VbdtdgNmvlVAwQoslg7tVebVlkJ1l2WSmWGfTgXduVZboP3I0+BzhueEpnj7Zy19ftI70FbJR3ra0/kl5rLgsaiZUFmYgJbQEWL0bNskZULoHzuk77x5DPZYM7tPTGm8XZgBOl5vazjEurlw8d212LckJ8bx/z2qeb+gLuj/pQ290kRQfx0cX3FOwIrr55OyfIKMYCvXbM3UpS7VsLlHug5zuHdf+ZOWXKOmQLh3axvzXSigYocWSwT0nJyfaJoTM0dZhPF7JOyrzF71uBS03H1jDrEfy4CutAY9p6hvn2dOj3HRgDQU67ZNqFLr4REpofUkZtRtURA2WaykvSCcpPk6f4L7hneB1Q8sL2reNNa4VtRihxZLBva/PXBtcqOGV5kGSEuLmyw74sIKWjcWZvH/Pah58pdVvOQIpJV/8zUlSEwT/56qNUbBQW3TxSe8JmOiH9Zdr3/YKLNWSGB9HZVEGjXoE9zUHlUygpme1bxtrXCtqMUKLJYP72rXWysSQUvKHU328ozKflMT4Re9ZRctnr1ZS9/75yVPLCk/96NU23mwb5q5Dmyw/agedfNL4FCBg49Xat70C/rRsLsmksdeh/cniE2HDldD0PHg13Ih7DqtcK2owQoslg/uZM2eibUJINPc7OTc0yTu3LM9ttYqW0pxU7rhyA0/W9XDPXICXUvLLYx187YkGrtxcxI4Me2xgrItPGp+EtRdCRqH2ba+APy2bSjLpc7gYndS4gBgoX17OPuit1bxpq1wrajBCiyVTGnbs2BFtE0Li+VPKTzB/wd1KWv72ig0MOmd48JVWnm/oIy0pnsbecS7ekM9/fnDvsl8lVkVzn4ycUzaRfvc92rarAn9aNs3dVG3sHZ8vEKcZG98FCDj9jJI9oyFWulaCYYQWS47crVb68/mGPnaWZVOSvTy31UpahBB85bqtfPvGnWwoyiAjOYF73redBz9ygJTEeEtpWQnNdZx+Snnc9B5t21WBPy26ZsykFyhZMyd/pdxE1hC79C8wRoslR+5WKv3Z55impmOUv7/Kf6EgK2kBJcD/5f41/KWfkrlW0xIIzXWc/DUUbYX8Sm3bVYE/LSVZKeSkJVLfPabPSXfcCL//O2UrwdLdmjVrl/4FsZK/AbHSN/jvarqQEq7b5X/FppW0BMMuWjTVMdikFAvbdYt2bYaAPy1CCPatzeVY24g+J91yvVIlsu6XmjZrl/4Fsc06AmKVb3ApJb863sWetTmsL/RfTMkqWtRgFy2a6qj5OYh42HmTdm2GQCAtF1TkcXZwgv5xHRbTpOUpc+8nfwVe7coc2KV/QWzkHpDaWu3vxOtBfbeD033j3LC3LOAxVtGiBrto0UyH1wO1DyuBLjM6FQ0DafHta3u0VafR+65bYLwHzjyjWZN26V9gjBZLBvdt27ZF2wRVPHa8k6T4OP58hSJaVtGiBrto0UzH6adgvBt2f1Cb9sIgkJbtpdmkJsbP71KlOZveA1ll8Mb9mjVpl/4FxmixZHBvbm4OflCUGZuc5dFjHVy7cxU5aUkBj7OCFrXYRYsmOqRU6pvnlkclS8ZHIC3Kaukc3lC5EUvIxCfAgY8rJRf6GjRp0i79C4zRYsngXlYWeJrDLPz8zXYmZzx84tKVi2hZQYta7KJFEx3tr0HXMbjo00qgixIrabmgPJ/GXgcjEzosZgLY+xFISIFX/12T5uzSv8AYLZYM7oODwWuLRxOX28MPj7RyyYYCtpVmr3is2bWEgl20RKxDSvjTtyEtP6pTMrCylis3FyElPNfQq8/J0/PhwCfgxMPQ3xhxc3bpX2CMFksG94wMfbbx0oofv3qO/nEXt/9Z8Lxms2sJBbtoiVjHmWfh7Atw6WcgKU0bo8JkJS3bV2dRnp/G72t79DPgkjshKQMOfyPipuzSv8AYLZYM7rOzs9E2ISBDThf/driJyzcVqtp31MxaQsUuWiLS4Z6B574I+Rvhgtu0MypMVtIihOC6XaW82jKoT0okKKP3d/wfaHwCWg5H1JRd+hcYo8WSwd2rQ8U5rfjOs6eZnPHwpWu3qDrezFpCxS5aItLxwjdhqBkO3atUSYwywbRct6sUr4SnTug4en/HHVBQBb+7A6bDr0Zpl/4FxmixZHBPS4vuT91APHOyl4ePdvCJSyvYUJSp6jNm1RIOdtESto6zL8KR7yk3Eje+S1ObwiWYlqriTLaVZvHj187h9ugUcBJT4H3/paSFPvmPYdecsUv/AmO0WDK4Dw/rlL4VAW2DE9z1qxPsLMvmH9+1SfXnzKglXOyiJSwdA6fhl7dCwUZl1G4S1Gi548qNnB2c4Lc13foZUrYfLv8C1D0KL98XVhN26V9gjBZLBvfS0tJom7CInrEpPvjAG8QJ+Leb95CUoP7PajYtkWAXLSHrGGyGn7xfqafyV49AUro+hoWBGi1Xbytm++osvvfHM7jc2pULWMZln1HKMBy+J6z0SLv0LzBGiyWDe2tr4L08jaa+e4wb/+s1HFOz/PhjBykvCO3CNpOWSLGLlpB0NP8RHrgS3NPw4V9D3nr9DAsDNVqEEHzu6s10DE/x9d9rs+AowIng+n+Hre+D574ET30WZtXfyLVL/wJjtIilW6b5PUiIQ8D3gHjgASnlt5a8nwz8GNgHDAE3SSnbVmpz//798tixY2EZ7fV6iYuL7vfSjNvLD4+08t0/nCEnNYn/+et97CwLfdNbM2jRCrtoUaVjvE9J73vrJ1C4Bf7qYWU1qskIxSfferqR+//Uwleu28pHL1558V1kRnnguS/D69+Hws1w9T9D5ZVBNw63S/+CyLQIIY5LKfcHOy5o60KIeOD7wDXAVuAWIcTWJYd9HBiRUm4Avgv8S+gmq6empkbP5lek3zHNAy+f5cp/fZF7n27kkg0F/P6OS8IK7BBdLVpjFy0BdXhmofVlePwO+H87lIqPF/8d3HbYlIEdQvPJZ6/exFWbi/ja7xu4+9d1jE3qlK4XFw+H/hk++CuYmYSf/gU8+C6o/glMBp6Ltkv/AmO0BB25CyEuAr4qpbx67vndAFLKexcc8+zcMa8JIRKAXqBQrtB4JCN3PfF6JVOzHiZcbkanZukcmaRjeIqWASfV7SOc7FJSuS4oz+NTV1RyxaaiKFscQzO8XpidgJkJmBoFRyc4upXUxo6j0P0WuKcgMQ12fEAJ7FHYgENP3B4v9z13hvv/1EJaUjzX7yrlHRsKWF+QzqrsFPLSkxBBRtihndAF1T+GN/8HBs8AAoq3Q/E2KN6qrBfIKIL0QqWUcEJqVMs5mAG1I3c1wf1G4JCU8hNzzz8MHJRSfnrBMSfnjumce94yd0zANbbhBvejv/4eBSf+GxEnYInpvi4nl76x8D2pHLH8M0r9da9c+JnzT+KEIDkh7v+3d+4xdhX3Hf9873Of3t27u7g2htiUV0xdG5t3UFVKQxsaEbVFhCRSkj6UPqgaSh+CRkpLVPWd9CFFzauK2ipJTSglEVEFARJFAkLiNTaLjXEMBmMb7N29+/Tu3uevf8zc9WW9u17Wy733XM1HujpzZubMme+9M787Z845v6ElFac9FScVX6KBL/qdnhmfy+VIp9NneTxsibQlf75FEld8rqWPyxfypJLznKS9Q+d6+8ed5ZjCjDPcCxFLwrqtcME1cMG17jHHBrppuhQDAwMr8h3+4hsTfOkHr/D4/hNM5opz8al4jNZUnGRcJOMxkvEYiZhOd6QlWDKLGZeVD3FtaYB3F1/iovKr9NnCo/gicXKkyCtFiThlPwFhiDKCylZ6S5wtp5I1Yv8FH+Z9v/HnKzp2uca9pn+Bkj4BfALc3eLh4WEKhQLlcpm2tjay2Szr16/n8OHDXH755ezZs4ft27fPNdDdu3eT6OzjZHojiUScUqlELB4HKwNCcoOvWDxGqVQimUhSKORJpVLk83lS6TS5XJ6WtNtPJlMUi0XiiTjlcplEzDVUYaRTScqlAt2d7ZRzM/T39jA6mqWnJ0N2NEumJ8Po2BhdXd1MTk3S1tZGLpcjkUhg5TIGxOMJ8vkcra1tTE1N0dW1huzoGL2ZDCMjI/T29jKVzZLq6WF8YoKOjg5mZmZIpdKUikWQiMVEsVgknW5henqazs5OxsbGyFSVMTySpa+3l9GxUdasWcOpU6doaWkhn88TT7ifuFQqkUqlmJ2Zpb29nfGJCTI9GYZHhunr7ZvbZrNZunt6mJiYoL29ndncLMlEknK5jJmRSCbJ5XK0tbUxOTlJd3c3IyMj9PX2Mj4yQn9fH8PDrl5jY2N0dnYyPT1NuqWFYrGIgFgsRqFYpKWlhalTp+ha08XoaJbe3l6Ghobp7+/z235Gslm6u7uYmpyitfIdx+OYRKlUIp1KMT0zQ2dHB2Nj4/T2Zjg5NMx5/f2cHBrivP5+hkdGyGR6GB+foL2jg9nZWZLJJOVSGcNIJBLMFqG9Zy3ZyVn6N2zitTdH2bT1RgaPZNlywy3s3jvItm3bOHDgAJtKMY6//DKZTIbp6WlisRjJZJKpqSn6+vo4evQoF198Mfv27WPr1q1z7beyHRwc5NJLL+XIkSOsXbuWsbExWlrc+rqzs7N0d3dz4sQJLrzwQg4ePMiWLVvOKGPv3r1cccUVHDp0iA0bNjA8PExHR8cZ/Wnz5s3s379/0f40p2nTJo4fPz6nqTsW41M3b+Cuq7sYLbfw/MvHoK2Hfa8cZU1PH8ffPEFXTy8nhoZZ093D+Pi4+61PTdPSkiZfKBCPx+faXjKZZHZ2lra2diYnJ+nq7mI0O0pPpmduOzb60/yw60oen5qirbWNxMwQ62yIzlKWztI4nTZFrDhDW6IE+VO0J6CYmyGdTlHIzZJOJSnkc6RSSUqFPIlEnHKxSDwmsFLFBmFlczaiWCKRSFAoFkgmkxTyeZKp1OltoUAikaBUKhGPxSlbGfk/CMOIKUapXCIej1MsFhctI5lInrYzpTKKiTX953Ps2LFl273q32nZmNmSH+B64NGq/fuA++bleRS43ocTwDD+qmCxz44dO2ylDAwMrPjYRiNoaTyaRYdZ0NKonIsWYJedxW6b2bKmZRLAQeBm4BjwY+DDZravKs9dwBYz+11JdwK/ZmZ3LFVu1J+WWS2ClsajWXRA0NKoNMTTMmZWBP4ANzp/EXjAzPZJ+oyk23y2fwd6JR0C7gHuXVGtl8mBA+fuPrRRCFoaj2bRAUFLo1ILLct6zv2d4FxG7jMzM7S2tq5yjepD0NJ4NIsOCFoalXPRsmoj90bk+PF30AdGjQlaGo9m0QFBS6NSCy2RNO6ZTKbeVVg1gpbGo1l0QNDSqNRCSySN+/T0dL2rsGoELY1Hs+iAoKVRqYWWSBr3ZrljDkFLI9IsOiBoaVRqoSWS31YyWf8VblaLoKXxaBYdELQ0KrXQUrenZSQNAa+t8PA+3ItSzUDQ0ng0iw4IWhqVc9HyLjPrP1umuhn3c0HSruU8ChQFgpbGo1l0QNDSqNRCSySnZQKBQCCwNMG4BwKBQBMSVeP+pXpXYBUJWhqPZtEBQUuj8o5rieSceyAQCASWJqoj90AgEAgsQUMad0kXSPqepP2S9kn6pI/PSPqupJ/4bY+Pl6R/lXRI0vOSttdXgUNSi6QfSdrrddzv4zdJetbXd6eklI9P+/1DPn1jPeu/EJLikp6T9Ijfj6QWSa9KGpS0R9IuHxep9lVBUrekByUdkPSipOujpkXSZf63qHwmJN0dNR0VJP2R7/MvSPqGtwW17SvLcfpe6w+wDtjuw504f/Kbgb8H7vXx9wJ/58O3Av+HW8nrOuDZemvw9RLQ4cNJ4FlfvweAO338F4Df8+HfB77gw3cCO+utYQFN9wBfBx7x+5HUArwK9M2Li1T7qqr3fwC/7cMpoDuqWnwd47h1mN8VRR3A+cBhoNXvPwB8vNZ9pe5fxDK/rG8B7wVeAtb5uHXASz78ReBDVfnn8jXKB2gDdgPX4l5eSPj4uZWuWMGKVjXWsAF4AvgF4BHfsaKqZSHjHrn2BXR5Q6J58ZHTUlWnW4CnoqrDG/fXgYxv+48Av1TrvtKQ0zLV+EuUK3Gj3rVm9oZPehNY68OVL7PCUR9Xd/w0xh7gJPBd4GVgzNwiKPDWus7p8OnjQG9ta7wk/wz8GVD2+71EV4sBj0kakFvbFyLYvoBNwBDwVT9d9hVJ7URTS4U7gW/4cOR0mNkx4B+BI8AbuLY/QI37SkMbd0kdwP8Ad5vZRHWaub+5hn/Ux8xKZrYNN+q9Bri8zlVaEZLeD5w0s4F612WVuNHMtgPvA+6S9HPViVFpX7iR3nbg38zsSuAU81ZCi5AW/Dz0bcA356dFRYe/L/AB3B/veqAd+OVa16NhjbukJM6wf83MHvLRJySt8+nrcKNhcGu7XlB1+AYf1zCY2RjwPdzlWLfc2rTw1rrO6fDpXcBIjau6GO8BbpP0KvDfuKmZfyGaWiqjK8zsJPC/uD/eKLavo8BRM3vW7z+IM/ZR1ALuz3a3mZ3w+1HU8YvAYTMbMrMC8BCu/9S0rzSkcZck3LqsL5rZ56qSvg18zIc/hpuLr8R/1N9Bvw4Yr7qUqxuS+iV1+3Ar7r7Bizgjf7vPNl9HRd/twJN+tFJ3zOw+M9tgZhtxl81PmtlHiKAWSe2SOith3BzvC0SsfQGY2ZvA65Iu81E3A/uJoBbPhzg9JQPR1HEEuE5Sm7dlld+ktn2l3jcfFrkhcSPu8ut5YI//3Iqbh3oC+AnwOJDx+QV8HjefPQhcVW8Nvl4/CzzndbwAfNrHXwT8CDiEu/xM+/gWv3/Ip19Ubw2L6Pp5Tj8tEzktvs57/Wcf8CkfH6n2VaVnG7DLt7OHgZ4oasFNX4wAXVVxkdPh63c/cMD3+/8C0rXuK+EN1UAgEGhCGnJaJhAIBALnRjDugUAg0IQE4x4IBAJNSDDugUAg0IQE4x4IBAJNSDDugZog6Z8k3V21/6ikr1Ttf1bSPat8zqnVLM+XuU3SrVX7fynpT5ZxnCQ9KWnNKtQhJekHVS/EBAJnEIx7oFY8BdwAICmGW/39iqr0G4Cn61Cvt8s23DsXb5dbgb02z43GSjCzPO7Z7w+ea1mB5iUY90CteBrnegGcUX8BmJTUIykNvBvYLalD0hOSdsv5W/8AgKS/lXRXpbDqEbOkP5X0Y+/X+/6FTr5QHkkb5fyff9n73n7Mv0mMpKt93j2S/sH75U4BnwE+6OMrxnWzpO9LekXSHy6i/yP4NxLPct7v+6ucXT7P1ZIekvNn/ldV5T3sywwEFiQY90BNMLPjQFHShbhR+jM4T5/XA1cBg35EOgv8qjmnXjcBn/WvcO8E7qgq8g5gp6RbgEtwvmG2ATvmOwE7S55LgM+b2RXAGPDrPv6rwO+Yc/pW8hrywKdx/ra3mdlOn/dynEvXa4C/8H6R5vMenGfACoudFyBvZlfhfH5/C7gL+Bng45Iq3gJfAK5e4DyBABCMe6C2PI0z7BXj/kzV/lM+j4C/lvQ87nXz83FuX58DzpO0XtJWYNTMXsf5hbkF5+ZhN87QXjLvvEvlOWxme3x4ANjo/QF1mtkzPv7rZ9H1HTPLmdkwzrHV2gXyZMxssmr/jPNWpX3bbweBfWb2hpnlgFfwDqbMrATkKz5yAoH5hBsygVpSmXffght5vg78MTCBGymDm2roB3aYWcF7oWzxad/EOVb6KdxIHtyfwd+Y2ReXOO+CeeTWCshVRZWA1hXoml/GQv2qKClmZuVFjqk+byWtPC9feV7ZadyVTiBwBmHkHqglTwPvB7Lm/NxncUvCXc/pm6ldOL/xBUk34ZZaq7AT55Hydk77+34U+E053/9IOl/SefPOu5w8c5hzzzwp6VofdWdV8iRu6ce3y0s4x1Grgp+eGTbnUjYQOINg3AO1ZBD3lMwP58WN+ykNgK8BV0kaBD6K86wHgJntwxnWY+bdu5rZY7hpk2f8MQ8yz/guJ88C/BbwZblVtNpxq+OAc9u6ed4N1eXwHZw3zdXiJl9mILAgwStkILAAkjrMbMqH78Wtz/nJcyhvHfCfZvbeVarfQ7iFow+uRnmB5iPMuQcCC/Mrku7D9ZHXcKvXrxgze8M/+rjmXJ91949kPhwMe2Apwsg9EAgEmpAw5x4IBAJNSDDugUAg0IQE4x4IBAJNSDDugUAg0IQE4x4IBAJNSDDugUAg0IT8PxT/zXXbm71+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(200, 800, 200)  # wavelength, units: nm\n",
    "absorption_spectrum = lumogen_f_red_305.absorption(x)  # units: nm-1\n",
    "emission_spectrum = lumogen_f_red_305.emission(x)      # units: nm-1\n",
    "plt.plot(x, absorption_spectrum, label=\"Abs (arb.)\")\n",
    "plt.plot(x, emission_spectrum, label=\"Ems (arb.)\")\n",
    "plt.xlabel(\"Wavelength (nm)\")\n",
    "plt.legend()\n",
    "plt.grid(linestyle=\"dotted\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The absorption spectrum is normalised to one so that it is easy to convert to unit you wish. Let's say units of length is now in centimetres as we construct the scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "You can open the visualizer by visiting the following URL:\n",
      "http://127.0.0.1:7016/static/\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<div style=\"height: 400px; width: 600px; overflow-x: auto; overflow-y: hidden; resize: both\">\n",
       "<iframe src=\"http://127.0.0.1:7016/static/\" style=\"width: 100%; height: 100%; border: none\"></iframe>\n",
       "</div>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "world = Node(\n",
    "    name=\"world (air)\",\n",
    "    geometry=Sphere(\n",
    "        radius=10.0,\n",
    "        material=Material(refractive_index=1.0)\n",
    "    )\n",
    ")\n",
    "\n",
    "sphere = Node(\n",
    "    name=\"sphere (glass)\",\n",
    "    geometry=Sphere(\n",
    "        radius=1.0,\n",
    "        material=Material(\n",
    "            refractive_index=1.5,\n",
    "            components=[\n",
    "                Luminophore(\n",
    "                    coefficient=np.column_stack((x, absorption_spectrum)),\n",
    "                    emission=np.column_stack((x, emission_spectrum)),\n",
    "                    quantum_yield=1.0,\n",
    "                    phase_function=isotropic\n",
    "                ),\n",
    "                Absorber(coefficient=0.1)\n",
    "            ]\n",
    "        ),\n",
    "    ),\n",
    "    parent=world\n",
    ")\n",
    "ray = Ray(\n",
    "    position=(-1.0, 0.0, 0.9),\n",
    "    direction=(1.0, 0.0, 0.0),\n",
    "    wavelength=555.0\n",
    ")\n",
    "scene = Scene(world)\n",
    "np.random.seed(0)\n",
    "vis = MeshcatRenderer(wireframe=True)\n",
    "vis.render(scene)\n",
    "for _ in range(10):\n",
    "    steps = photon_tracer.follow(scene, ray)\n",
    "    path, decisions = zip(*steps)\n",
    "    vis.add_ray_path(path)\n",
    "vis.vis.jupyter_cell()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiple components can be added simultaneously,\n",
    "\n",
    "```python\n",
    "components=[\n",
    "    Luminophore(\n",
    "        coefficient=np.column_stack((x, absorption_spectrum)),\n",
    "        emission=np.column_stack((x, emission_spectrum)),\n",
    "        quantum_yield=1.0,\n",
    "        phase_function=isotropic\n",
    "    ),\n",
    "    Absorber(coefficient=0.1)\n",
    "]\n",
    "```\n",
    "\n",
    "this material simulates and radiative and non-radiative absorption. For example, non-radiative absorption occurs in the host material of glass of PMMA. Non-radiative absorption by the luminophore itself can be include by reducing the `quantum_yield` keyword."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happened above:\n",
    "\n",
    "* The green ray enters the luminophore sphere. One ray refracts as it exists, unabsorbed (it still remains green).\n",
    "\n",
    "* Other rays are absorbed an remitted at longer (orange) wavelengths and eventually escape.\n",
    "\n",
    "* One ray is emitted into a trapped optical mode (red)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}