Melanin.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MONTE CARLO SIMULATION OF LIGHT-TISSUE INTERACTION: EFFECT OF MELANIN ON LIGHT ABSORPTION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Epidermis is the outermost layer of the skin and it acts as the physical and chemical barrier between the interior body and the exterior environment. It primarily consists of keratinocytes, melanocytes, Langerhans cells, and Markel cells. Melanocytes produce melanin- a natural pigment responsible for the color of our skin, hair, and eyes. Due to its light absorbing characteristics, melanin acts as a photoprotective substance of the human skin. Through a process called melanogenesis, melanin protects our skin by absorbing and reradiating harmful UV radiation and reduces the chance of melanoma.\n",
    "\n",
    "<img src=\"images/Epidermis.jpg\" width=\"400\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this project, I wish to simulate the light absorbing characteristics of melanin. However, modeling the light-tissue interaction in a turbid media like epidermal tissue using conventional deterministic method is difficult because of the inherent stochasticity of the process [1]. When light is incident on the tissue surface, it can randomly be reflected from the surface or refracted into the tissue. Light that finally enters the tissue are scattered randomly within the tissue. Some light may transmit through the tissue, some are reflected into incident medium while some are absorbed by the tissue.\n",
    "\n",
    "<img src=\"images/LightInTissue.png\" width=\"400\" />\n",
    "\n",
    "Due to this stochastic nature of light-tissue interaction, I will use Monte Carlo method for simulating the light absorbing behavior of melanin. Monte Carlo method is a computational technique, invented by Stanislaw Ulam in 1946, that uses repeated random sampling to estimate some functions of probability distribution. I will be using the random walk variant of the Monte Carlo method which relies on taking repeated random steps through a space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Walk of Photon in Tissue\n",
    "\n",
    "The random walk of light in epidermis follows the flowchart below. \n",
    "\n",
    "<img src=\"images/BlockDiagram.png\" width=\"700\" />\n",
    "\n",
    "A narrow-collimated photon is launched into the tissue and initially it is directed downwards. A random step size is then calculated using, and a check is done to see whether the photon hits the boundary or not. If the photon hits any of the top or bottom boundary, necessary adjustments to the step are made, and the photon is either reflected into the tissue or transmitted to the top or bottom medium. If the photon transmits to other medium, it is terminated. Otherwise, the photon is moved within the tissue. At the new position, a fraction of the photon energy is absorbed by the tissue and the photon is scattered in random direction. A random step size is calculated again and the photon wanders around the tissue in this manner until its energy has been fully absorbed. In that case, the photon is assumed to be dead and is terminated. The next photon is then launched. [2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I have implemented each step of the random walk in separate functions and then aggregated them in a single *RandomWalk* function. You will find these functions in their respective scripts in the [functions](functions/) directory of this repository. Each function is accompanied with a docstring, describing the function definition, input and output arguments."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook will walk you through the simulation. For avoiding verbosity, I have omitted the comments and docstrings in this notebook. The simulation can also be run just by executing [Melanin.py](Melanin.py) script."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import Necessary Libararies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import seed\n",
    "from random import random\n",
    "from math import floor\n",
    "from math import copysign\n",
    "from math import sqrt\n",
    "from math import log\n",
    "from math import pi\n",
    "from math import sin\n",
    "from math import cos\n",
    "from math import tan\n",
    "from math import asin\n",
    "from math import acos\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from tabulate import tabulate\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Launching Photon Into the Tissue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def LaunchPhoton():\n",
    "    \n",
    "    w = 1.0\n",
    "    \n",
    "    x = 0.0\n",
    "    y = 0.0\n",
    "    z = 0.0\n",
    "    \n",
    "    dpx = 0.0\n",
    "    dpy = 0.0\n",
    "    dpz = 1.0\n",
    "    \n",
    "    return w, x, y, z, dpx, dpy, dpz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate a Random Step Size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RandomStep(mu_att):\n",
    "    \n",
    "    s = - log(random()) / mu_att\n",
    "    \n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checking Whether the Photon Hits the Boundary or Not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BoundaryCheck(w, x, y, z, dpx, dpy, dpz, s, dDR, dDT, d, n_t, n_m, alpha_critical, epsilon = 1e-9):\n",
    "    \n",
    "    while True:\n",
    "        \n",
    "        if dpz < 0:\n",
    "            sbound = (0 - z) / dpz\n",
    "        elif dpz > 0:\n",
    "            sbound = (d - z) / dpz\n",
    "        else:\n",
    "            break\n",
    "        \n",
    "        if sbound > s:\n",
    "            break\n",
    "        else:\n",
    "            x += sbound * dpx\n",
    "            y += sbound * dpy\n",
    "            z += sbound * dpz\n",
    "            \n",
    "            s -= sbound\n",
    "            \n",
    "            alpha_i = acos(abs(dpz))\n",
    "            \n",
    "            if n_t > n_m:\n",
    "                if alpha_i > alpha_critical:\n",
    "                    dpz = - dpz\n",
    "                    break\n",
    "            else:\n",
    "                alpha_t = asin(n_t * sin(alpha_i) / n_m)\n",
    "                \n",
    "                R = 0.5 * ((sin(alpha_i - alpha_t))**2 / ((sin(alpha_i + alpha_t))**2 + epsilon) + (tan(alpha_i - alpha_t))**2 / ((tan(alpha_i + alpha_t))**2 + epsilon))\n",
    "                \n",
    "                if R < random():\n",
    "                    if dpz < 0:\n",
    "                        dDR += w\n",
    "                    else:\n",
    "                        dDT += w\n",
    "                    \n",
    "                    w = 0\n",
    "                \n",
    "                dpz = - dpz\n",
    "    \n",
    "    return w, x, y, z, dpx, dpy, dpz, s, dDR, dDT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Taking a Step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TakeStep(x, y, z, dpx, dpy, dpz, s):\n",
    "    \n",
    "    x += s * dpx\n",
    "    y += s * dpy\n",
    "    z += s * dpz\n",
    "    \n",
    "    return x, y, z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Absorbing Photon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AbsorbPhoton(w, z, pDW, absFact, dz):\n",
    "    \n",
    "    dw = w * absFact\n",
    "    \n",
    "    if (w - dw) < 0:\n",
    "        dw = w\n",
    "    w -= dw\n",
    "    \n",
    "    nz = floor(z / dz)\n",
    "    pDW[nz] += dw\n",
    "    \n",
    "    return w, pDW"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scattering Photon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ScatterPhoton(dpx, dpy, dpz, g): \n",
    "    \n",
    "    if g == 0:\n",
    "        theta = acos(2 * random() - 1)\n",
    "    else:\n",
    "        theta = acos((1 + g**2 - ((1 - g**2) / (1 - g + 2 * g * random()))**2) / (2 * g))\n",
    "    \n",
    "    psi = 2 * pi * random()\n",
    "    \n",
    "    \n",
    "    if abs(dpz) > 0.9999:\n",
    "        dpx = sin(theta) * cos(psi)\n",
    "        dpy = sin(theta) * sin(psi)\n",
    "        dpz = copysign(cos(theta), dpz)\n",
    "    else:\n",
    "        dpx = sin(theta) * (dpx * dpz * cos(psi) - dpy * sin(psi)) / sqrt(1 - dpz**2) + dpx * cos(theta)\n",
    "        dpy = sin(theta) * (dpy * dpz * cos(psi) - dpx * sin(psi)) / sqrt(1 - dpz**2) + dpy * cos(theta)\n",
    "        dpz = -sin(theta) * cos(psi) * sqrt(1 - dpz**2) + dpz * cos(theta)\n",
    "    \n",
    "    return dpx, dpy, dpz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregate the Functions Together to Take a Random Walk!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RandomWalk(mu_att, absFact, g, n_t, n_m, alpha_critical, d, dz, dDR, dDT, pDW):   \n",
    "    \n",
    "    w, x, y, z, dpx, dpy, dpz = LaunchPhoton()\n",
    "    \n",
    "    while True:\n",
    "        s = RandomStep(mu_att)\n",
    "        w, x, y, z, dpx, dpy, dpz, s, dDR, dDT = BoundaryCheck(w, x, y, z, dpx, dpy, dpz, s, dDR, dDT, d, n_t, n_m, alpha_critical)\n",
    "        \n",
    "        if w <= 0:\n",
    "            break\n",
    "        \n",
    "        x, y, z = TakeStep(x, y, z, dpx, dpy, dpz, s)\n",
    "        w, pDW = AbsorbPhoton(w, z, pDW, absFact, dz)\n",
    "        \n",
    "        if w <= 0:\n",
    "            break\n",
    "        \n",
    "        dpx, dpy, dpz = ScatterPhoton(dpx, dpy, dpz, g)\n",
    "    \n",
    "    return dDR, dDT, pDW"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effect of Melanin Concentration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first simulate the effect of melanin concentration on light absorption. I considere an in vitro setup where the epidermis is covered in a pair of coverslips.\n",
    "\n",
    "<img src=\"images/Setup.png\" width=\"400\" />\n",
    "\n",
    "The thickness of the tissue $d$ is $100 \\mu m$. For $550 nm$ light, the scattering coefficient $\\mu_s$, refractive index $n_t$, and anisotropy of scatter $g$ for epidermis is $74.7 cm^{-1}$, $1.447$, and $0.759$, respectively. The refractive index of the coverslip is $1.5$. The melanin concentration has been divided into three levels- light, medium, and high, corresponding to $1.3-6.3\\%$, $11-16\\%$, and $18-43\\%$ volume fraction of epidermis occupied by melanin, respectively. Their corresponding absorption coefficient $\\mu_a$ is given in table below. [3]\n",
    "\n",
    "|Melanin Concentration|Light | Medium | High  |\n",
    "|---------------------|------|--------|-------|\n",
    "|$\\mu_a (cm^{-1})$    |$1.63$|$6.85$  |$16.13$|\n",
    "\n",
    "\n",
    "For each level of melanin concentration, 5 million photons will be launched into the epidermis and the following metrics will be scored-\n",
    "- how much photon energy have escaped from the boundary of the tissue where the light was incident (diffuse reflectance)\n",
    "- how much photon energy have escaped from the opposite boundary (diffuse transmittance)\n",
    "- how much photon energy have been absorbed at various depths (absorption)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Melanin Level      Diffuse Transmittance    Diffuse Reflectance    Absorption\n",
      "---------------  -----------------------  ---------------------  ------------\n",
      "Low                             0.388032              0.332539       0.279429\n",
      "Medium                          0.159724              0.157917       0.682359\n",
      "High                            0.042999              0.0681065      0.888894\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fcfdf540970>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZzN1f/A8deZxWBGGHuEsa9jJzshEVIpO/FLK0VJUUqUUIoWhBZkTalos2/ZZedry5KdYWQ2s53fH+fOxhgzcz93m3k/H495uJ/PvffzeX/M8r7nc855H6W1RgghhLgTL1cHIIQQwr1JohBCCJEmSRRCCCHSJIlCCCFEmiRRCCGESJOPqwNwhIIFC+rSpUu7OgwhhPAoO3fuvKK1LnTr/iyZKEqXLs2OHTtcHYYQQngUpdSp1PbLrSchhBBpkkQhhBAiTZIohBBCpClL9lEIIbKfmJgYzpw5Q1RUlKtDcXs5c+akRIkS+Pr6puv1kiiEEFnCmTNnyJMnD6VLl0Yp5epw3JbWmpCQEM6cOUNQUFC63iO3noQQWUJUVBQFChSQJHEXSikKFCiQoZaXJAohRJYhSSJ9Mvr/JLeeMik21nzNmgU5c0K+fHDwoNlXuDA0agQBAVCiBFy/Dvnzg7e3q6MWQoiMy1KJQinVEehYrlw5y46ptfnjf/EijBgBc+ZA2bJw5gzcvJmxY3XrBjVrQv36EBgI1apJ8hAiKwkICCAsLMzVYVguSyUKrfVSYGndunUH2HusmBjYuRMaNrz9uePHkx57eUF8PNSpA6VLQ5UqULs2/PorzJwJ9eqZFsWRI7BggflKEBho3vvii6YF0qIF5M5tb+RCCGEt6aNIRUQEtGt3e5L44gv44w/45x84fdq0NuLizL87dsDixTB6NHTuDDNmmP3btsHhwxAWBpcvw59/JiUDf38IDYX334eHHzbbzZqZW1ibNmW8xSKEcD+nTp2iVatWBAcH06pVK06fPk1cXBxlypRBa01oaCheXl6sX78egKZNm3Ls2DEXR51SlmpRWGHCBHj99aTtfv3g66/tP66/v/l68EEID0/af+4cPPOMSTQXL8KGDVC1qnkuf37z9dln5n0+8t0SIn0GD4bdu609Zs2aMGlSht82cOBA+vTpQ9++ffn666956aWX+Omnn6hQoQIHDx7kxIkT1KlThw0bNtCgQQPOnDmDlbfPrSAtimTWrk1KEq+/DtHR1iSJtNx7LyxbBhcumNtda9eaFoe/v0kS//xjWhtly0LjxrBnD4SEODYmIYR1Nm/eTI8ePQDo3bs3GzduBEzLYf369axfv57hw4ezceNGtm/fTr169VwZbqrkM6qN1jBqFFSqBFu3wj33OD8GHx9o3jypxREdDXv3wi+/wLx55nZUzZrmuTZtTB9IyZLOj1MIt5eJT/7OkjA0tWnTpkybNo1z584xevRoPvzwQ9auXUuzZs1cHOHtpEVhoxQsXWo6oV2RJFKTIwfUrWv6Pfbvh++/T3puxQooVQpq1YI1a0yHuRDC/TRq1IgFtlEsc+fOpUmTJgA0aNCATZs24eXlRc6cOalZsyZffvklTZs2dWW4qZJEkUyePFCmjKujSF3OnNCli2n5XL4MEyea/bt3wwMPmHkcuXKZhCKEcI2IiAhKlCiR+PXxxx/z6aef8s033xAcHMycOXOYPHkyAH5+ftx3333cf//9gGlh3Lhxg+rVq7vyElKltNaujsFydevW1dll4aLvv4cnnzStoP/+M/sKFzZJZeJEk2CEyA4OHTpE5cqVXR2Gx0jt/0sptVNrXffW10qLwsM98YRpZVy/bm6dlSwJly7BlCmmhdGwoRmeK4QQmSWJIgvp0AFOnYLNm+Gpp0zS2LLFdNBXqADffGNmmQshREZIosiC7r/fJIUDB0zJEYCjR6F/fzMcd8kSWLXKtTEKITyHJIosLCAAevQwfRXHjplSIZcvw2OPQevWULw4/Pyzq6MUQrg7SRRZnJcXvPKKmbD3+edJJURy5jSzwjt3NkOD9+xxdaRCCHcliSKbSSghEhkJJ04k7U+oart/v/RjCCFSkkSRjZUubWZ/f/aZqWS7fTtUrw5NmyYNtRVCpJ9Sit69eydux8bGUqhQITp06JCh47Ro0YKEIf7t27cnNDTU0jgzShJFNufrCwMHmr6LP/4w5dK3bIFy5cytqvPnXR2hEJ7D39+f/fv3ExkZCcCKFSsoXry4Xcf87bffyJcvnxXhZZokCgGYvoy2bU0V26VLoVgxGDTIjJJq1MiUPhdC3F27du349ddfAZg/fz7du3dPfC48PJz+/ftTr149atWqxc+20SSRkZF069aN4OBgunbtmphoAEqXLs2VK1c4efIk1apVS9z/0UcfMWrUKMC0QIYMGUKzZs2oXLky27dv57HHHqN8+fK89dZbdl+TFAUUt+nQAR56yAyjnTrV1JKqWtXMy/jlF6hRw9URCpG2wX8MZvcFa8uM1yxak0kP3b3YYLdu3Rg9ejQdOnRg79699O/fnw0bNgDw/vvv88ADD/D1118TGhpK/fr1ad26NV9++SW5c+dm79697N27l9q1a2c4vhw5crB+/XomT57MI488ws6dOwkMDKRs2bIMGTKEAgUKZPiYCaRFIVLl42Nmfa9ebYbWFi5sFmuqWdOUB7l2zdURCuGegoODOXnyJPPnz6d9+/Ypnlu+fDnjxo2jZs2atGjRgqioKE6fPs369evp1atX4vuDg4MzfN5OnToBUL16dapWrUqxYsXw8/OjTJky/Pvvv3Zdk7QoxF2VLWv6Kr76CqZNgx9+MF8rVkCDBqaYohDuJD2f/B2pU6dODB06lLVr1xKSbAEZrTU//PADFStWvO09CeXH78THx4f4+PjE7aioqBTP+/n5AeDl5ZX4OGE71s6hjNKiEOni5QUDBpg+jEGDoEgRsyZG0aLw8cdm+VghhNG/f3/efvvt2yrBtm3bls8++4yEYqy7du0CoFmzZsydOxeA/fv3s3fv3tuOWaRIES5dukRISAg3b95k2bJlDr6KJJIoRIYoBZ9+am5HjRplEsWrr5oV+YYNg2QfeITItkqUKMHLL7982/6RI0cSExNDcHAw1apVY+TIkQA8//zzhIWFERwczIQJE6hfv/5t7/X19eXtt9+mQYMGdOjQgUqVKjn8OhJImXFhl7g4UyZk0aKkfT/9BI884rqYRPYkZcYzRsqMC6fx9oaFC80CSh07mn2dO5tRUkeOuDY2IYQ1JFEIS9SoYYbOXrwItWubeRcVK5pbVd9+a9bMEEJ4JkkUwlKFC8POnWaEVIJ+/VJuCyE8iyQK4RD9+5uO7cceM9sDBpjWxYABpl9DCOE5JFEIh1HKzLc4dixp38yZZjKfbSSgEMIDSKIQDle2rKlSe/QoJAyy6NULJk2S4bRCeAJJFMIpfH1NRdoDB2D5cmjXDoYMMUUHbfXThPB4AQEBKba//fZbBg4cCMC0adOYPXt2mu9P/np3IolCOJVSZkb30qWmH+PSJVOEsEIFWLdORkeJrOu5556jT58+rg4jUyRRCJfw9jYjoSIizPyLo0ehRQuzTOuVK66OTgjrjRo1io8++giA7du3ExwcTMOGDXnttddSlA8/d+4cDz30EOXLl2fYsGGuCjcFKQooXCpnTjP/YsUKs0xrVBQUKgTvvgtDh5rEIURGDR5sJoFaqWZN06+WlsjISGrWrJm4ffXq1cSqrsn169eP6dOn06hRI954440Uz+3evZtdu3bh5+dHxYoVGTRoEPfdd58l15BZ0qIQbqFNG7NW9zPPmO133oGWLc3Ke0J4ily5crF79+7Er9GjR9/2mtDQUG7cuEGjRo0A6NGjR4rnW7VqRd68ecmZMydVqlTh1KlTTok9LdKiEG7D2xu+/BImToSnnzalQQoXhilT4LnnTP+GEOlxt0/+rnS3+nrJS4R7e3vbXSLcCtKiEG4nIAAWLIBly6B0aXjhBVPm/PnnTatDCE+WP39+8uTJw5YtWwBYsGCBiyO6O0kUwm09/DAcPgwJt3CnTYOgINiyBWJiXBubEPb46quveOaZZ2jYsCFaa/LmzevqkNIkZcaFR4iIgLFj4f33k/bdvAk5crguJuFePKnMeFhYWOKci3HjxnH+/HkmT57s1BikzLjIcnLnhvfegx9/TNrn52dKggjhaX799Vdq1qxJtWrV2LBhA2+99ZarQ0qTdGa7ieNXj3OP3z2cun6Kned20rhkY4rnKc7p66fZc3EPQfmCyOWbC4B789xLREwEQfmC8PbydnHkzvXoo2Yk1P/9nxlWO2AAnDsHXbpAlSqujk6I9OnatStdu3Z1dRjpJonCwWLjY7kScQVfL19Co0JZdGARh0MOc+7GOaLjoll3ap1dxy+TvwwnQ09SMm9Jbty8wbjW4yiTvwzz9s2jTZk2NCjRgKIBRfH18s0ySaVgQfj5Zzh/HurXN0Np33nHzLv44ANTdFBkT1prlAyPu6uMdjlIH4WFtNacvXGWd9e+y8Xwi6z8ZyWRsZGpvtbXy5eY+KQeWX9ff8JjwgkuEkz3at0JyBFAWHQYR0OOsvDAQvLlzMfZG2epUKACR0Iyt3Rcz+o9ebDsg8TExdCnRh/O3jhL6XylM3Usd3H8OPz2G4webWZ058ljhtW2a+fqyISznThxgjx58lCgQAFJFmnQWhMSEsKNGzcICgpK8dyd+igkUVggXsez6MAiuv/Q/bbnvJU3Par3IE+OPEzZMYWPH/yYxyo/Rql8pRJfEx0XTQ7v9PfKxut4vJQXWmsiYiK4GH6R/27+R0RMBHsu7OFk6Ekal2zMowsfJV6nrzzrPX730KNaD+rcW4dewb3w8/bzuF+2GjVg717z+M03TZ+GyD5iYmI4c+YMUVFRrg7F7eXMmZMSJUrg6+ubYr/HJgqlVBngTSCv1rpLet7jjEShtWbO3jm8tuI1LoVfSvFcTp+cfPzgxzxX9zm3+GMbERPBhlMbqFyoMuM3jmfazmn0qdGHb3d/e8f3VChQgVcbvsqDZR8kIEcABXMXdF7AmXT5smlZfP652X7qKXjySTPrW25HCXF3bpUolFJfAx2AS1rrasn2PwRMBryBmVrrccmeW+wOieJU6ClazGrBydCTKfbP7DiTpqWaUuKeEuTyyeUWCSI9tNacCD3B2pNrOXzlMBM2TUj1dRUKVCAoXxANSzSkZ3BPigYUJSBHQKqvdbW9e6FWraS1Lho3hvXrzaQ9IcSduVuiaAaEAbMTEoVSyhs4ArQBzgDbge5a64O2512aKI6GHKXr4q7surArcV/L0i0Z3XI0VQtVJX+u/Jaez1W01mw9u5UKBSqw7ew2/jr9FyGRIUzdMTXV11cpVIXaxWozpuUYt+rv+O8/mD4dPvoILl6E6tXNqnrVq7s6MiHcl1slCgClVGlgWbJE0RAYpbVua9seDqC1/sC2nWaiUEo9AzwDULJkyTpWFtJadGARXRcnDWUbVH8QfWv0pc69dSw7hye4FnmNtSfX8tPhn5i95/YFWGoUqcF9ee9jdIvRVCpYKXE4ryvFx8PAgTB1KpQsCY8/Dvffb25JCSFS8oRE0QV4SGv9tG27N9AAeAd4H9PSmJmQONJiVYvi/I3z3PvxvYnbo1uMpmdwT8rkL2P3sbOCg5cPUipvKX49+iuvr3z9tttxANULV+fjth/Tukxr5weYzKZN0LYthIUlbTds6NKQhHA7njAzO7Wb+lprHaK1fk5rXTY9ScIqp0JPMfjPwQBULVSVnc/sZGTzkZIkkqlSqAr+Ofx5suqTnHj5BJFvRrK813ICcwUmvmbfpX20mdOGwPGBVJ9anQ//+pDw6HA2nNqQ7hFZVmjUCI4dM+sUALRqBRs3Ss0oIdLDnVoUad56ygh7WxR7L+6l+bfNCY0KpVPFTvzU9SeP6Zx2B/E6ngthFyjiX4S1J9ey9uRavtj+BdeirqV4nY+XD5v/bzPFAopR/J7iTotv61bTurh+3Ww//bTpz5BvscjuPOHWkw+mM7sVcBbTmd1Da30go8e2J1FcDr9M7em1iYuPY2WflVQqWAkv5U4NL88UERPB+Rvn+WjTR2z8dyP7L+1P8XyeHHkY2mgoCsUTVZ+gYoGKDk3Oly+bNQvGjk3ad+KEKWsuRHblVolCKTUfaAEUBC4C72itv1JKtQcmYYbHfq21fv/OR7mzzCaKm7E3uf+r+9l9YTer+qzigaAHMnN6kU77L+1n6eGljFg94rbnHqn4CEMbDaXxfY0dmjB27YLatc3jUqVMKfOHHnLY6YRwa26VKBwts4mi2TfN2HB6A5PaTuLl+192QGTiTmLiYjh+7Tgz/57JrD2zuBJxJcXzszvP5qFyD1Ewd0GHJI5XXoFPPjGP33kHRoyQEuYi+8kWiUIp1RHoWK5cuQFHjx7N8Pt/PPQjCkXnSp2lT8LFwqPDGbdxHAsPLOTo1aTvZRH/IjQr1YyB9QfSrFQzy86nNVy9Ct26wcqVkD8//P47NGhg2SmEcHuWJgqllD8QpbWOsyI4q8nCRVlHXHwcfxz7g7n75jJ///wUz7Uu05oaRWrQoUIHmpdqbklyj42Fl14y8y4SnDsHxYrZfWgh3J5diUIp5QV0A3oC9YCbgB9wGfgNmK61zvhHeAeRRJF1/XjoR0asGsGViCuERIakeO75us8zvvV48vjlsfs8b76Z1NGdI4dZIKl3b7sPK4RbszdRrANWAj8D+7U2A+CVUoFAS6AHsERr/Z2lUWeSJIrsISQihImbJ/LBxpQjqMsFluPRSo/yeOXHqVe8XqZHrYWEwDffwGuvme2OHWH2bMiXz97IhXBP9iYKX611mlOT0vMaZ5FEkb1ordl/aT+fbPmEqNio225Rvd3sbercW4cOFTpkKmls3gwjR8KqVWb75EkzQkqIrCZbdGYnkESRvYVHhzNsxTCm7JiSYn+h3IUYUHsALUq3oE3ZNhk6ptbQrx/MmmW2e/aEL78Ef3+rohbC9SxLFEqp17XW4y2LzAEkUQgwCePLnV9yKvQUn277NMVzPav35OnaT9OidIt0Hy95gUGAunVh8WJpXYisI9OJQim1KPkmUFNrXd7i+CwliULc6kLYBebtm8ery1+97bmnaz3N1A5T8Vbe6Ro5FRlpbkVNnAgFCsAvv5ikIfMuhKezJ1HMTKjoatueqrV+3gEx2s3eeRQie4iIieB/V/7Hw/Me5kLYhcT9BXIV4Jfuv3B/ifvv2pehtakZ1a0bJFS0378fqlZ1ZORCOJY9iSJIa30i2Xag1vqqA2K0jLQoRHrFxccxfed0XvjthcR9zUs158GyD/Jyg5fxz5F2J8Rff0GTJuZx7twwdKjpy5CaUcIT2d1HoZQqqLW+cvdXup4kCpFR8Tqe34/+zg+HfuCb3d8k7q9QoAKT2k6iZVBLcvrkTPW9sbGwfDl06WJuS/n7w5kzEBAga3ULz2JFovhFa93J8sgcQBKFsEd4dDjrTq1j6eGlzN03lxvRN/BW3rQu05p699ZjaKOh5M2Z97b3HT8O5colbTdpYuZdBAU5MXgh7GBFoliqte5oeWQOIIlCWOW/m//x1d9fseR/S9hwekPifj9vP2Z2mkn3at3x9vJO3H/tGrzwAixYYLbz5TM1pKR0mPAEVqxwl/UmXAhxF/f43cOQhkNY3289u57dxdgHTF2Pm3E36b2kNz5jfBi1dlTi6/Pnh/nz4exZsx0aCi++CLt3uyB4ISwiLQohMigyJpLNZzbT9ru2xMbHAlAwd0Geqf0M9YrXo3OlzgBERUH79rBmjXlfy5bw55/g6+uqyIVImxW3nqpprfff/ZWuJ4lCOEN0XDRTtk9hw+kN/Hjox8T9QfmCKBdYjjmPzqFQ7iJ88okZDZXgyhUz/0IIdyMlPIRwoA2nNrDsyDImbJpAQI4AwqLDEp9b1GURD5V8gscfhxUrzD5vbzMyqmhRFwUsRCosSRRKqbrAm0ApwAczU1trrYOtCtQeMuFOuIs1J9bwwOyUS+m+WO9Fal2cxNP9k8bMHjwIlSs7OzohUmdVojgMvAbsA+IT9mutT1kRpFWkRSHcgdaalf+s5Invn+D6zeuJ+xsUaU7B5UtZviyA/PkVP/8M99/vwkCFsLEqUWzUWjexNDIHkEQh3ElsfCyHrxxmzck1DPp9UOL++rGvsGvCBGKizfDao0dTzsMQwtmsShStgO7AKswqdwBorX+845tcQBKFcFcRMRGsPrGaiZsnsvbkWggtCZ8dhjgz63vzZmldCNexYh4FQD+gJvAQ0NH21cH+8ITIHnL75qZDhQ6s6buGiBERDGzTCd7MDe1NramGDaFGw8ucOh1/lyMJ4TwZTRQ1tNZ1tdZ9tdb9bF/9HRKZEFlcLt9cfNb+M+JHxfH6kHvg2ZqQ7wR7txSidCkvHhu5kCw4KFF4oIwmii1KqSoOiUSIbEopxbjW44iZsoMfNu3Cv8lXACx5ryul2/7EGyuGs/rEahdHKbKzjPZRHALKAicwfRRuNTw2gfRRCE83/8fr9HjcVnjQNxyqzyNHpyGMbjWSYY2HpWuBJSEyyqrO7FQXfZThsUJYb8cOeHXEddavsCWM/Mfh2dqULJKP8oHl+bDNh9QqVsu1QYosxarO7NHAda31KVty+A94x4oAhRAp1a0L65bn5c03bTuulYVx1zkdcoFVJ1ZRe3ptfMf4MmfPHJfGKbK+jCaKYK11aMKG1voaIB9phHCg996DpUuTth8+cJV3qy0kIEcAsfGx9PmpD+U/K8+2s9vIiiV5hOtlNFF4KaXyJ2wopQIxpTzcglKqo1Jq+vXr1+/+YiE8SIcOcPmyWW711yX+vNPlSUb7h/DDkz8AcOzqMRrMbIDXaC8+3fopcfFxLo5YZCUZ7aPoAwwHFmPWp3gSeF9r7VZtX+mjEFnVhQvQqRNs35607/2Pr3Cu4ihOXT/FsiPLEvf3Cu7FgNoDaFaqmQsiFZ7IsuqxtuGxD2BGPK3SWh+0JkTrSKIQWZnW8NFHMGxY0r4nn4S5cyE0+go9fujBin9WpHjPit4raFqyKX4+fk6OVngSuxKFUkrpu7wwPa9xFkkUIjuYNAnGjjW3pAAGDzYJxNsbjl89zsTNE5m6Y2qK9xx44QBl8pchp09OF0Qs3J29o57WKKUGKaVK3nLQHEqpB5RSs4C+VgQqhEifwYPh0iUICTH1oSZNgm7dICwMygaWZcrDU7jw6gVGNhuZ+J6qU6pyzwf38Nfpv1wYufA06W1R5AT6Az2BICAUyIVJNMuBL7TWbrMqsLQoRHZz4wa0aQNbt0L58vDtt9CoUdLzWmuOXj3KK3++wq9HfwWgVtFaTO84nXKB5QiNCqV0vtIuiV24Dyv7KHyBgkBk8qGy7kQShciO4uLgl1/g+efh4kWzLyQEAgNTvu6jTR/x2orXbnv/ldeuUCC3rNGanVk14Q6tdYzW+ry7Jgkhsitvb3j0Udi1yzwGszZ3eHjK1w1tNBT9jmbdU+sonqd44v6CHxZk8B+DiY6LdmLUwhNkOFEIIdxbsWIQEwMBAWa7USNYnUpNwWalmnHmlTPodzSLn1hMLp9cTN46Gb/3/Gj8dWOuR8l8JGFIohAiC1LK9FuMHQt790KrVvDpp9yxbPnjVR4nbEQYYx8YC8CmfzeRb3w+7p95P5fDLzsxcuGOMtxH4Qmkj0KIJKtXm0QBUKUKrF0LhQql/Z4F+xfQ/YfuidtVC1WlYYmGDKgzgPrF6zsuWOFSVlWP9QMeB0qTrHSH1nq0BTFaRhKFEClFR0PfvrBggdn+4AMYOhR80ijAE6/jGb5yOBM2TUixv3mp5sx9bC7F7yl+h3cKT2VVZ/bPwCNALBCe7MstSK0nIVKXIwd89x288YbZHj7c1I+KjLzze7yUF+PbjCfu7TjGtx6fuH/dqXWU+KQEg34bRFh0mIMjF+4goy2K/Vrrag6MxxLSohDizuLioHjxpCG0+/ZBtXT+Vu86v4sl/1vCmPVjACjsX5iXG7zM/9X6Pwr7F5YFlTycVS2KTUqp6hbFJIRwAW9vOH06abtWLVi+PH3vrVWsFqNbjiZ2ZCwzOs5Aa82bq9+k6MSiDFg6gJ3ndjomaOFSGU0UTYCdSqnDSqm9Sql9Sqm9jghMCOE4OXKYIbTTpkFsLLRtC6++eudRUbfy9vLm6dpPc/7V8zxS8REAvtr1FXVn1GX2ntlExqRxT0t4HFkKVYhsbvlykygSfPcd9OyZsWNsO7uNp356ikNXDiXua1C8AXMfm0vZwLIWRSoczZJbT7aEkA/oaPvK525JQgiRMQ8+CIcPJ03Q69ULnnkGzp9P/zHqF6/PwRcPEjIshOfqPAfA1rNbKfdZOfze82P72e13OYJwZxlKFEqpl4G5QGHb13dKqUGOCEwI4TwVKpgJejNmmO0ZM+Dee+HIkYwdJzBXIFM7TOXCqxdY2GUhJe4pQXRcNPVn1qfzgs6cCpXPlZ4oo7ee9gINtdbhtm1/YLPWOthB8WWK3HoSIvOOHIHgYLh502x//jm8+GLmjhUWHcayI8tSTN4b2nAoQxsNpUhAEQuiFVayatSTApIvxhtn2yeEyCIqVDDzK7p2NdsDB5qSIGGZmDIRkCOAbtW6cWbIGXoF9wLgo80fUXRiUdS7inUn11kYuXCUjCaKb4CtSqlRSqlRwBbgK8ujEkK4lFJmFndERNK+KlXg+PHMHa/4PcWZ8+gc1j+1nqdqPpW4v8WsFoxYNYKsWEooK8nMehR1gMaYlsR6rfUuRwRmD7n1JIR1Nm2Cxo3N45IlYcsWU6HWHr2X9Oa7vd+l2Df2gbEMbzrcvgMLu1i2cJEnkEQhhLXi4+Hxx+GnnyBvXpg1Cx55xP7jnr5+mlKTUo6671atG5+0/YSiAUXtP4HIELv6KJRSG23/3lBK/Zfs64ZS6j+rgxVCuBcvL/j+e7PEav780LkzjB5tFkWKi7vr2++oZN6SRL0ZxcIuC2lRugVgKtfWnFaTJYeWWBK7sJ+0KIQQGbJvnxkVlaBHD5g715pjbzu7jb9O/8XQFUOJ1/EA7H52NzWK1rDmBCJNlox6UkqNT88+IUTWVa0aTJ0KhQub7Xnz4IUXrM0gDHoAACAASURBVDl2/eL1GdJwCOdeOZe4r+aXNan1ZS0mbZmUmDyEc2V01FObVPa1syIQIYRnUAqeew7OnoVt28y+qVNhwgT4z6Ib0UUCinBp6CWGNRoGwO4Luxny5xC8R3sze89sGSXlZOm69aSUeh54ASgDJAyQU0AA8JfWupfDIswApVRHoGO5cuUGHD161NXhCJEtvPmmWXIVoH59WL8e/PysPceJayd46uenWH9qfeK+v/r/RaP7Gll7omzOrlFPSqm8QH7gA+CNZE/d0FpftSxKi0gfhRDOdeVKyuVVv/sOunSxNmFEx0Uzau0oPtj4QeK+coHl2PPcHnL75rbuRNmYXX0UWuvrWuuTQH/MUqifAB8DfZVSOa0MVAjheQoWNLO5H3rIbPfqBeXKWXcrCiCHdw7GthrLxaEXqVywMgDHrh7Df6w/8/bNk9tRDpTRPopZQFXgM+BzoDIwx+qghINcv24WH0jNnfYLkU45c8Lvv8PmzWYI7Zkz0KaNtckCzKp6B188yPJeSast9fyxJ16jvfjt6G/WnkwAGS8KuEdrXeNu+1wt2956unIFTp409aKXL4erVyF3bihSxNxEPnIESpc2v72rV5veyJIlzbqYf/1leinr1jUD42Njzf59+8zjKlXA3x9CQsyN6OhoOHcO7rkH+vUzK97Uq2cG1ufKBZcvw333ufp/RLiI1mbuBYCPDwwebDrAy1q8NMWxq8f46u+vmLh5IjHxMYBJJLuf3U3B3AXx9fa19oRZnCUzs5VS3wLTtNZbbNsNgL5aa4sGx1kjSyaK8+fNH2mtYedO80d6/XqzVNmePeZjW2ioNefy8bGmhVGzpvloefmyiX3gQHMje88ec4+ibVvwlV/krOqXX+CLL0xDdutW01+xcydUreqY843bOI7hq1KWAAkbHoZ/Dn/HnDALsipRHAIqAgkr7pYEDgHxgHaXcuMenyi0hu3bYcUK+OEH2JVGOa0KFSAqyvxRLlwYgoJMq6FgQahYEQ4ehGvXzMe7pk1NoZ7Chc1r8uc3rYSYGHPOf/8120WLmtbFqVPm2BERppWxd6957p9/TAujTBlYtsy8PyLCPD5+3CQGMEksOjrta82Xz7RKGjUyj2NioF07c4wqVUzMYOIRHklrU+5j6VKzffOm+dFwlDUn1vDA7AcAqF64OjM7zaR+8fporVHyc5QmuxOFMv/DTYE7rjziLqvdeVyiOH/e9P6tXp368/nymdZCu3YmEZw+bRLE22+bwjvu7PJl85dhwwbzuFo1OHYM3n/fXIefn0lc//6b9nFy5zb1I/LkgQIFTIunYkUIDHTKZQj7XL9ufowT/P47tGhh+jUcIS4+jmeXPct3e7/jZtzNxP0HXzhI5UKVHXPSLMCqFsVOrXUdSyNzALdPFHFxMHEifP212T58OOXzjRubT9C9esHTT4O3t/lYlpU/DV28CPv3m2u8csUkz8GD09cqqVXLJNCaNc1iz0WKmATk4+Oc2EW6hIaaRmyCevVMx7e3t+POefjKYSp9USnFvtV9VtMyqKXjTurBrEoUXwDfaq3degFct0sUCbdvzp0zf8j+97/b//g99RR07w7Nm1s/W8nTXb9ubnf99RccOgRz5piPp8WLw44d5mNpaqvq5MxpbrG9/LL5q1SnTtZPuG7uxg1o3970VURGmsl6773n2HNGxkRy9sZZ2s1tx7Grx/BSXoxuMZrXGr9GDm8H3gPzQFYlioOYPoqTQDhmdrbb9E0kcJtEER8Pr7wCkyen3F+vnhkR1K2bGYEUHm4WKJY/YJkTH2/6Yf74AxYtMr2oYIbYJF9pp2RJc7srwfDh8NJLpgUi//dOpTW0bAnr1pnG4MKF5tfCkWLjYzl+9Ti9lvRixznz96FfzX6Maz2Owv6FHXtyD2FVoiiV2n536ZtI4NJEERkJixeb4R5bt6Z87tVXTfW0MmVcE1t2ERJi+jHADBFeudLc3tu2zXS4pyYwEHr3hoYNTd9HmTKmk104zNat8NprZrxFSIhp/G3aZP+iSOnRYGYDtp01hapy++bm6KCj3JvnXsef2M1ZtnCRUqoGplMbYIPWeo8F8VnKJYkiNhbeegvGJyummy8ffPKJ6XMoX9658YjUxcfDzz/Dl1+ajvHo6KQWyK1y5DA30J97ztwyLFkyZZ0KYYmVK03DGkw307Ztjh81fePmDa5GXqXvT31Zd8qs2922bFtmdZ5FkYAijj25G7OqRfEyMAD40bbrUWC61vozS6K0iFMTRWys+Ul/6inTIQvmJuzkyebWh9zS8Axam36QefPMaLI7adXKfK8rVza3C319zVBkYZcrV6BSJdOyCAw0rYz4eMe3LrTWvPjbi0zdMTVxX5sybZj80ORsOTrKqkSxF2iotQ63bfsDm7NlH0V0tOl/WLPG/FSDude9a1fSPAThmSIjTb9RYCD89pvp+4iKgvnzzaCE1LRrB82amduLMokwU5LP5k6+zxmOhBxhyvYpTN6a1J+4pu+axFX3sgurEsU+oJ7WOsq2nRPYrrWublmkFnBoooiLMzOLu3SBEyfMvtq1TbnMytnvE0i2ExJiJgHMm2f+Tcu+feaWY1ycmQci7urMmZSVXyZNMoPWnGX+vvn0WtIrcYGk8a3H07J0S+oVd3BPu5uwKlG8AvQFlmBGPD2CGS47yapAreCwRHHhAnTsaIZkJtixwwy7FNlPVJQZgnv2rOmFnTHDzKZPrmBBkygGDoQnnzQTCYoVM++V5JGqpUvN0qoLF5rtJUvMzG5nNtIXHVjE878+z9VIs4rC5+0+54V6L2T5md1WdmbXBprYNtdrrXdbEJ+lHJIofv0VOnRI2h46FIYNk85NkVJsrClD8vPPMG2a2bd1q0kMyeXPD7Nnm1tW8fFyuyoVGzaYu3kAH39s5l868+90VGwU/X/uz/z98wEIyhfEjI4zaBnUEi+V0cLbnsGqFsUTwB9a6xtKqZFALWCM1jqNYkTOZ3mi+O4704EZF2e2z51zzhg+kTWcOWOG5d64YeZ5xMfD33/f/rqaNU2/V8eOKetdZGOrVkHr1ubxvfeauZbNmzt2NndyMXEx/Hr0Vz7f9jmrTqwCoFVQK/7o9Qc+Xllv5v+dEgVa63R/AXtt/zYB1mNuPW3NyDEc+QV0BKaXK1dOW2bMGK1Nn5rWM2dad1yRvUVEaP3GG1o3b6516dJJP2PJv954Q+s//9Q6MtLV0brU7t1aN2mS9N+ilNbr1zs3hujYaP351s81o9CMQpedXFavP+nkIJwA2KFT+dua0RbFLq11LaXUB8A+rfW8hH32ZjIrWdKiCAkxq8VPmGC2v//edGAL4Sjnz0OTJmaY7q1atjSDJsaMMet9ZLNSJOHhplG/eHHSPkdXob2Tz7Z+xkt/vARAtcLV2Pf8PucH4SB2LYWazFml1JfAk8BvSim/TBzDMwwcaJJE8eLmF1iShHC0YsVMZd0TJ8yM8iFDoEcP89yaNaaQZO7cJkF4eZne3lv7PrIof3/zWe3sWbNuFphakGvXOj+WQQ0GMavzLAD2X9pPvRn12HPB7eYdWyqjLYrcwEOY1sRRpVQxoLrWevld3upUdrcoJkyA1183tZhmz5aORuF6UVGmxMipU2atkuQee8wUTmzcGN591zXxOVFUlJkkn7DsyfjxZlyJs10Kv0Tr2a3Zd8m0KKa0n0Kbsm0oF1jO+cFYxLJRT57ArkSxZIn5xata1SzyExBgbXBCWCEsDPr3Nx+zlUo5M83Pz4yqWrXKLP6URc2eDX37msfVqsHu3c7r5E4Qr+OZt28evZf0TtxXtVBVfuv5GyXzlnRuMBaw5NaTUiqnUuoVpdSPSqkflFJDbJPusobNm02SKF/ePJYkIdxVQIAZQRUXl7T6YKdOZiGrmzfNnJ+qVc3oqYAAM2vt9Om7r+3hQXr3NiW7wCxl4uNjlm8JCXFeDF7Ki17BvTg9+DSDGwwG4MDlA7T9ri1XIq44LxAHy+itp0XADeA7267uQH6t9RMOiC3TMt2iGDAAZs40v3TV3WqyuRDpd/06fPWVuU01Z44pwZ7ciBHwxBPmZ9zZH8EdYNs2M8Up4VZUo0Zm6RJX6Lq4K4sOLErcPvDCAaoU8pxWnVXzKPZorWvcbZ+rZTpRjBljfqk+/tj6oIRwlZgY6NMHFixIuT842BSu7NMHGjQwrREPni3epElSghg2zHwlVJt3poiYCJ78/kl+PforAF93+pp+tfo5P5BMsCpRfAtM01pvsW03APpqrV+wKlAruM3CRUK4k+vXzUS/oUPNMKING25/zY8/mhWESpRwfnx2+ucfk/eSi4gwo4mdLSImggFLBzBv3zwA6hevz+89fycwl3uv8W5XorAVA9SAL2aFu4RlwkoCB7XW1SyM1W6SKIRIh7g40+Hdvn1S1YEE5cqZDoBSpcwiTh4yZyMmxqw7lVCvs3x5M7K4eHHXxBMZE8lTPz+VeDuqYYmGrOyzkty+7tlyszdRpLaynQJKACO01u3tD9E6kiiEyIRTp+DZZ+HPP1Pur1/f3KZq0sRM+nPz/rsbN8waYp9+mrRv2jRzaa6gtabUpFL8+9+/iftW91lNy6CWrgkoDVYWBawJ9MBMujsB/KC1/tySKC0iiUIIO5w9a4obrltnquHOnZty+G316jB6tPNLumZQdDQ0bWo6u8Gscrt8uesGM/78v5/pvLAzAN7Km3+H/EvRgKJuVZHW3hZFBaAbZpRTCLAQGKq1TnUNbVeTRCGEha5ehUuXzCJOr76a8rmXXjJJw9/fjE91Q99/byq8g0kWmza5LpaYuBhyvZ+LOG1u9Y1oMoLRLUfj7eUeo8/snUfxP6AV0FFr3USbpU/j7vIeIURWEBho1il95RVTzuaLL8wMNzD3d/LlM9ULqlUzJUjczBNPJFWg3bwZpky5fcSws/h6+7LoiaThs2M3jqX7D91x94nP6W1RPIppUTQC/gAWADO11kGODS9zpEUhhIPFx8PXX5tJfDNnmgSSoEULUyq9Vy+TZNygpREfDx98YHLc+fNmJFR4uGvvnF0Ov8zTS5/ml8O/ALCk6xI6V+rsuoCwbnisP9AZcwvqAWAWsCTL1XoSQmTMiRNmabpb1y2tXBl++skMPwKX92mcP2/WtQAzmOu778ztKFeJ1/H0WdKHufvmAnBk4BHKFyjvsngsKeGhtQ7XWs/VWnfAjHjaDbxhUYxCCE8VFGT6K+LjzZryTz9t9h86ZMarenmZv9CHD8PBgy4Ls1gxs45U8+Zm3kWLFrevXutMXsqLOY/OoUHxBgBU+LwC3x/4ngthF1wXVCoyXSJca31Va/2l1voBKwMSQngwpcxQ2hkzzEipw4dNFWYw9acqVTI1qJSCZ54xHeVOVry4KU9erJgZGfXgg6b7xVWUUmx5egtty7YF4MnFT/LyHy/f5V3OlTXXkhBCuIcKFWD+fJM0jh6F//u/pOdmzDA1NsaOhZMnnR7apElJjz/5xOSu3393ehiJfu72Mx8/+DG+Xr4sOrCID//60HXB3ELKjAshnG/OHFNjKoGvLwwaBIUKmRlzb7wBefI4JZTIyJQlrg4dMg0fV9l4eiNNv2kKwBuN36BTxU40vM85HSmyHoUQwv1cvWrWfXntNTO0NqEMeoEC5uN9vXpOCSM+3nRsJ6xvMWnS7f3yzvT3+b+pM70OAD5ePkS9GeWUuRZWLYUqhBDWCQw0taYOHDDraAwcaPaHhJjSIUqZv94O/kDr5ZWywO7IkbBxo0NPmabaxWqzsvdKulfrTmx8LD5jfPjj2B8ui0daFEII93Lpkrn/M39+0spE/v7w8MPmq2lTU93WQUsU79+fVM5q6FD40IVdBVGxUfiP9SdexwOwod8GmpRs4rDzya0nIYTnefZZWLkS7rvP1J5KEBxs7g+1dExhvd9/Nw0dgDffNCN8e/Vy3TSQMevG8PbatwG4Ouwq+XPld8h5JFEIITzbhQtmcbEpU1Lu/+knqFvX8lriFy+afHTpktles8bMu3CFkIgQCn5YEID7S9zP+NbjaVaqmeXnkT4KIYRnK1rU1OC4eTPlKpSdO5tbUSVLwrffWna6IkXM5LwEn3ySlDScrUDuAuh3NA+Xf5gtZ7bQ/NvmfLDhA6edXxKFEMKz5MgBQ4aYDu7du5P2//sv9Otnhth+840lHeC+vqbqeps28MsvJnk0agShoXYfOlO+aP9F4uMRq0dw7sY5p5xXEoUQwnPVqGESQlgYLFpkFle6cgX69wc/P0uGLt17b8q1nDZvNod3hVL5ShE+IpyqhaoCUPzj4lwKd3wzJ0slCqVUR6XU9OvXr7s6FCGEM/n7m3ri69aZWlNg1kVt2tRM3OvYEWbNyvThlTJTPN55B2rVgiVLYMQIi2LPoNy+uVnVZxWN72sMwIhVjg9EOrOFEFnTqVNmtvfSpUn7evc2K/M1aGD6NTLh9GmzlDjAzp1mdVhX6bq4K4sOLOLPXn/yYNkH7T6edGYLIbKXUqVMx8L16/C//5kkMWcOdOlihtu+9hpk4gNlyZKwerV5XKeOaW1cvGhx7On0ebvPqVigIl0Xd2XHOcd9OJZEIYTI2u65x0yEmD0bzp0zEyIAPvrIlAhp0AAuXzY1ptKpRYuk5VXBDMhas8basNOjkH8hfuz6I6FRofRZ0oeYuBiHnEcShRAi+yhWzLQqrl2Dt80ENrZtg8KFTUI5ciRdh1EKFi406zIlePHFpFJVzlSlUBU+b/c5h64cov8vjulll0QhhMh+8uWDd981ieGRR5L2V6xopmRv3pyuw+zfb+pDPfecqTri5+eSJTZ4sf6LPF/3eb7b+x0nrp2w/PjSmS2EEGFhZsm7WrVMKdkEe/aY6dl3ERdnahj+/beZxjFzJrRr57ByVKk6ff00M3bO4MX6L1I0oGimjiGd2UIIcScBASYhREWl3F+jhikPciLtT+ne3rBpk1lG4/Jl00h5/HEHxpuKknlLMuaBMZlOEmmRRCGEEAl8fSE83MzBSFgfdedOKFPGlJRNY0q2nx988IEpPQVmVG716nD8uBPidjBJFEIIkVzu3ODjAxMnmkKEbdqY/fv3Q/78Zlr2jz+apfFS8cgj5qlixcxbypUzk8U9mSQKIYS4kyJFYPly2LsXunc3+775xtxXyp0b9u1L9W05c5qRuAmlyidMcFK8DiKJQggh7qZ6dZg3z8z2TliFD0y/xgsvmPIgBw/e9rZly8w0jQ8/hKlT4fBhJ8ZsIRn1JIQQGaU1jB0Lb711+/5b7NqVssxHo0awYYNZftXdyKgnIYSwilJm6bu1a1Pur1TJdH4nU6tWUskPMKOjxoxxfIhWkkQhhBCZ1by5aUVcuGCmaR8+bIbTjhxp7jvFxgJmxdbYWDO3AmDUKNetaZEZkiiEEMJeRYrA1q1mIgXAe++Z0ua+voktDG9v+O23pEZIo0ZJFdHdnSQKIYSwQp48ZiJFTIxJFAnatDFDamfOBKBZM2jVypT8ePhhU7bc3UmiEEIIK/n4mP6L8+fNv9eumSG1AwbAH3+grl1l5UqzxtLZs65bACkjJFEIIYQjFC1qWhY7dsDzz5t97dpBgQLw3380a2aWxJg714yA+vdf14abFkkUQgjhSHXqwJQpKYs/5c0LISEMG2amYmhtZnS762wFSRRCCOEMixebMrMJChak4Lih7Ppb89RTZr6Fl5epDuKqFfPuRBKFEEI4i5eXqVA7frzZnjgRLx8vvnznXOJLHn/cDJhyJ5IohBDCmfz8YNgw08ldvToAOYKKc6lSM2pUMUvkbd8OZ864MsiUJFEIIYQr5MsHf/1l+iv8/Sn0vw3sOJibZV+eBcw8C1eslpcaSRRCCOEqefKYKdphYfDII/gQx8PPluDX8fv4918YPdrVARqSKIQQwh389BN8/z0A7V8PpnrBc0yeDEOHpuwDdwVJFEII4S66dDFTtvv25ecrjQGzflKVKhAS4rqwJFEIIYQ7qVQJJk0i6OGqLPbrCcCRIzDxI9dNspBEIYQQ7iZfPli2jMdnPMR/5AHgg3GKTye7JllIohBCCHfVqxd5PnyHAXnmAzByeAyrVjk/DEkUQgjhrpSCoUOZfr4TO+7tRO7IENo9GMu/J53buy2JQggh3J2/P3Vmv8wvdCIm3oeSQd78c9x5t6EkUQghhCdo1Yo6YesJxAx/ahx8g1OnnHNqSRRCCOEhvPxzERIVQAmvs1yIuIfSpZ0zx0IShRBCeBI/P76anSNxc0HXJQ4/pSQKIYTwMA/2LMS1v09QipN88sN9xK7d6NDzSaIQQggPlK9WEC/3DGEndZnVfqFDF7GQRCGEEB5q8Jw6FAqM5enIz3jpvh/hn38cch5JFEII4aGUgiFDfQD4LOZ5Dtfs6pD1VCVRCCGEBxs+HN4aboY+jbrxiln1yGI+lh9RCCGEU40Z681/12OZNqMrl4O8KGTx8SVRCCFEFvDMCz6cOgv//QeFLM4UkiiEECILqFrVrH3kCG6fKJRS/sAUIBpYq7We6+KQhBAiW3FJZ7ZS6mul1CWl1P5b9j+klDqslDqmlHrDtvsxYLHWegDQyenBCiFENueqUU/fAg8l36GU8ga+ANoBVYDuSqkqQAngX9vLXLxyrBBCZD8uSRRa6/XA1Vt21weOaa3/0VpHAwuAR4AzmGQBMpxXCCGczp3+8BYnqeUAJkEUB34EHldKTQWW3unNSqlnlFI7lFI7Ll++7NhIhRAiG3GnzmyVyj6ttQ4H+t3tzVrr6cB0gLp167puFXIhhMhi3KlFcQa4L9l2CeCci2IRQghh406JYjtQXikVpJTKAXQDfnFxTEIIke0p7YACUnc9qVLzgRZAQeAi8I7W+iulVHtgEuANfK21fj+Tx78MZHaRwILAlUy+11PJNWcPcs3Zgz3XXEprfdu8bpckCnemlNqhta7r6jicSa45e5Brzh4ccc3udOtJCCGEG5JEIYQQIk2SKG433dUBuIBcc/Yg15w9WH7N0kchhBAiTdKiEEIIkSZJFEIIIdKUrRLFHcqYJ3/eTym10Pb8VqVU6WTPDbftP6yUauvMuO2R2WtWSrVRSu1USu2z/fuAs2PPLHu+z7bnSyqlwpRSQ50Vsz3s/LkOVkptVkodsH2vczoz9syy4+faVyk1y3ath5RSw50de2al45qbKaX+VkrFKqW63PJcX6XUUdtX3wyfXGudLb4wk/iOA2WAHMAeoMotr3kBmGZ73A1YaHtcxfZ6PyDIdhxvV1+Tg6+5FnCv7XE14Kyrr8fR15zs+R+A74Ghrr4eB3+PfYC9QA3bdoFs8HPdA1hge5wbOAmUdvU1WXTNpYFgYDbQJdn+QOAf27/5bY/zZ+T82alFcacy5sk9AsyyPV4MtFJKKdv+BVrrm1rrE8Ax2/HcXaavWWu9S2udUGvrAJBTKeXnlKjtY8/3GaVUZ8wv0gEnxWsve673QWCv1noPgNY6RGvtCWu+2HPNGvBXSvkAuTArZ/7nnLDtctdr1lqf1FrvBeJveW9bYIXW+qrW+hqwglvWA7qb7JQo7lTGPNXXaK1jgeuYT1npea87sueak3sc2KW1vumgOK2U6Wu2Lbv7OvCuE+K0ij3f4wqAVkr9abtlMcwJ8VrBnmteDIQD54HTwEda61vXxnFH9vwNsvvvlzuVGXe0VMuYp/M16XmvO7Lnms2TSlUFxmM+fXoCe675XeATrXWYrYHhCey5Xh+gCVAPiABWKaV2aq1XWRui5ey55vqYlTLvxdyG2aCUWqm1/sfaEC1nz98gu/9+ZacWRXrKmCe+xtY0zYtZic9TS6Dbc80opUoAS4A+WuvjDo/WGvZccwNgglLqJDAYGKGUGujogO1k78/1Oq31Fa11BPAbUNvhEdvPnmvuAfyhtY7RWl8C/gI8oRaUPX+D7P/75epOGid2Bvlg7j0HkdQZVPWW17xIyg6wRbbHVUnZmf0PntHpZ88157O9/nFXX4ezrvmW14zCMzqz7fke5wf+xnTq+gArgYddfU0OvubXgW8wn7L9gYNAsKuvyYprTvbab7m9M/uE7fud3/Y4MEPnd/V/gJP/s9sDRzCjB9607RsNdLI9zokZ7XIM2AaUSfbeN23vOwy0c/W1OPqagbcw93J3J/sq7OrrcfT3OdkxPCJR2Hu9QC9Mx/1+YIKrr8XR1wwE2PYfsCWJ11x9LRZecz1M6yEcCAEOJHtvf9v/xTGgX0bPLSU8hBBCpCk79VEIIYTIBEkUQggh0iSJQgghRJokUQghhEiTJAohhBBpkkQhsj2lVJxSaretguoepdQrSqlM/24opUYke1xaKbU/ne8brJTqk9nzJjvOQKVUP3uPI0QCGR4rsj2lVJjWOsD2uDAwD/hLa/2OBccrDSzTWle7y3t8MJPfamtTmyjTlFK5MfHXsuc4QiSQFoUQyWhT1uEZYKAyvJVSHyqltiul9iqlngVQSrVQSq1XSi1RSh1USk1TSnkppcYBuWwtlLm2w3orpWbYWizLlVK5Ujn1A8DfCUlCKVVOKbXS1sL5WylV1nbOdUqpRUqpI0qpcUqpnkqpbbb1FcrariECOKmU8oQKx8IDSKIQ4hbaFIjzAgoD/wdc11rXw8x8HaCUCrK9tD7wKlAdKAs8prV+A4jUWtfUWve0va488IXWuioQiqnGe6vGwM5k23Nt76kBNMJUOwWoAbxsO2dvoILWuj4wExiU7P07gKaZ/C8QIgVJFEKkLqHi5oNAH6XUbmArplR1edtz27RZHyAOmI+pxJqaE1rr3bbHOzELzNyqGHAZQCmVByiutV4CoLWOsrUSALZrrc9rU/L9OLDctn/fLce9hKmQKoTdslOZcSHSRSlVBlOK+hImYQzSWv95y2tacHup5jt1+CVfxyMOs2DOrSIx9Ykg9bLQqR0rPtl2PCl/n3PajimE3aRFIUQySqlCwDTgc21GevwJPK+UmehUjQAAAN5JREFU8rU9X8G2wBFAfaVUkG2EVFdgo21/TMLrM+AQUA5Aa/0fcMa22l7C+s+5M3i8CphCf0LYTRKFEEmdzwcwpbaXk7TK3UxMldG/bcNcvyTpk/tmYBzmD/IJzNodANOBvck6s9Pjd6BZsu3ewEtKqb3AJqBoBq+pse1ahLCbDI8VIhNst56Gaq07WHjMJcAwrfVRO49TC3hFa93bmshEdictCiHcxxuYTm17FQRGWnAcIQBpUQghhLgLaVEIIYRIkyQKIYQQaZJEIYQQIk2SKIQQQqRJEoUQQog0/T/JjE7uA2Hg5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sMelaninLevel = ['Low', 'Medium', 'High']\n",
    "mu_as = np.array([1.63, 6.85, 16.13])\n",
    "mu_s = 74.7\n",
    "g = 0.759\n",
    "n_t = 1.447\n",
    "n_m = 1.5\n",
    "d = 0.1\n",
    "\n",
    "mu_atts = mu_as + mu_s\n",
    "absFacts = mu_as / mu_atts\n",
    "\n",
    "if n_t > n_m:\n",
    "        alpha_critical = asin(n_m / n_t)\n",
    "else:\n",
    "    alpha_critical = float('nan')\n",
    "\n",
    "N_bins = int(1e3)\n",
    "dz = d / N_bins\n",
    "\n",
    "pDW = np.zeros((3, N_bins))\n",
    "dDR = np.zeros(3)\n",
    "dDT = np.zeros(3)\n",
    "\n",
    "N_photons = int(5e6)\n",
    "\n",
    "for n, (mu_att, absFact) in enumerate(zip(mu_atts, absFacts)):\n",
    "    for nph in range(N_photons):\n",
    "        dDR[n], dDT[n], pDW[n] = RandomWalk(mu_att, absFact, g, n_t, n_m, alpha_critical, d, dz, dDR[n], dDT[n], pDW[n])\n",
    "\n",
    "\n",
    "dDR /= N_photons\n",
    "dDT /= N_photons\n",
    "dDW = np.sum(pDW, axis = 1) / N_photons\n",
    "\n",
    "print('\\n')\n",
    "print(tabulate({'Melanin Level': sMelaninLevel, 'Diffuse Transmittance': dDT, 'Diffuse Reflectance': dDR, 'Absorption': dDW}, headers = 'keys'))\n",
    "print('\\n')\n",
    "\n",
    "pz = np.arange(0, d, dz)\n",
    "pDW = pDW / (dz * N_photons)\n",
    "plt.figure()\n",
    "plt.semilogy(pz, pDW[0], 'r', label = sMelaninLevel[0])\n",
    "plt.semilogy(pz, pDW[1], 'g', label = sMelaninLevel[1])\n",
    "plt.semilogy(pz, pDW[2], 'b', label = sMelaninLevel[2])\n",
    "plt.xlabel('Depth (cm)')\n",
    "plt.ylabel('Absorption (cm$^{-1}$)')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is evident from the table above that the higher the melanin concentration, the higher the light absorption. The plot demonstrate the depth profile of absorption in the epidermis. Absorption rate is higher near the surface, and it drops off as we go down the epidermis. The rate drops off more quickly for highly pigmented epidermis than medium and lightly pigmented one. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effect of Wavelength"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, I will demonstrated the wavelength dependency of light absorption by melanin. I have only considered lightly pigmented epidermis for this part. The wavelengths of light $\\lambda$ being used, and the corresponding optical properties of epidermis are given in the table below. All the other parameters have the same value as stated above. [3]\n",
    "\n",
    "|$\\lambda (nm)$|$\\mu_a (cm^{-1})$|$\\mu_s (cm^{-1})$|$g$  |$n_t$|\n",
    "|--------------|-----------------|-----------------|-----|-----|\n",
    "|350           |30.16            |210.4            |0.702|1.488|\n",
    "|450           |13.5             |121.6            |0.728|1.448|\n",
    "|550           |6.85             |74.7             |0.759|1.447|\n",
    "|650           |3.68             |55.48            |0.787|1.433|\n",
    "\n",
    "For each wavelength of light, 5 million photons have been launched into the epidermis and their absorption depth profile has been recorded. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "  Wavelengths    Diffuse Transmittance    Diffuse Reflectance    Absorption\n",
      "-------------  -----------------------  ---------------------  ------------\n",
      "          350              0.000318329               0.129227      0.870455\n",
      "          450              0.024086                  0.150934      0.82498\n",
      "          550              0.159541                  0.15793       0.682529\n",
      "          650              0.362692                  0.166306      0.471001\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fcfb4b4c670>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd1zVZfvA8c/NXoKyFAEVFJy4cK80R27TslyZ+ahlPaaNX8/T1pa2s51Wmo+llZUrd64y98Q9EBVRQRFB2XD//rgVNRfjwGFc79eLF+f7Ped8z3XIznXudd1Ka40QQgiRVzbWDkAIIUTJJAlECCFEvkgCEUIIkS+SQIQQQuSLJBAhhBD5YmftAIqKt7e3rlatmrXDEEKIEmXr1q1ntdY+N7uvzCSQatWqsWXLFmuHIYQQJYpS6tit7pMuLCGEEPlS6hOIUqqXUmrKhQsXrB2KEEKUKqU+gWitF2itR3l4eFg7FCGEKFXKzBiIEKLsysjIIDo6mtTUVGuHUmw5OTkREBCAvb19rp8jCUQIUepFR0dTrlw5qlWrhlLK2uEUO1przp07R3R0NEFBQbl+XqnvwpIxECFEamoqXl5ekjxuQSmFl5dXnltopT6ByBiIEAKQ5HEH+fn7SBdWfmVkwKZNcOAAXLpkftLTISwMqlaFKlXA29vaUQohRKGRBJIXWVnw7bcwdy6sXQsXL97+8VWrQrNm0LkzdOgANWoUTZxCiGIlNTWVdu3akZaWRmZmJvfffz8TJkwAYNiwYaxZs4YrvSTTp0+nYcOGaK0ZO3YsixYtwsXFhenTp9O4cWNrvo0bSALJraNH4eGH4c8/ISQEHnoIOnaERo3A3R1cXSExEQ4ehDNnICoKtmyBNWvg55/NNUJCoFs3GDrUPM+m1PcgCiEAR0dHVq5ciZubGxkZGbRp04Zu3brRokULAN59913uv//+656zePFiDh06xKFDh9i4cSOjR49m48aN1gj/liSB5Ma+faYFkZIC06ebBHCz/kJnZ6hY8fpzWpuks3UrrFgBn34KH38MNWua6/TvbxKLEKLUUkrh5uYGmCnFGRkZdxxzmDdvHkOHDkUpRYsWLUhISODUqVP4+fnlPCYqKopu3brRpk0b/v77b/z9/Zk3bx7Ozs60b9+eRo0asXXrVuLi4pgxYwYTJ04kIiKCBx98kDfeeKPA76vUJxClVC+gV438dh9lZZnWRlYW/P031K2b1wCgXTvz89RTEBMDP/wAv/0GL75ofpo2hX//GwYOhDzMwRZC5MO4cbBjh2Wv2bAhfPTRbR+SlZVFeHg4hw8f5oknnqB58+Y597344ou89tprdOzYkUmTJuHo6MjJkycJDAzMeUxAQAAnT568LoEAHDp0iFmzZjF16lQeeOABfvnlF4YMGQKAg4MDa9euZfLkyfTp04etW7fi6elJ9erVeeqpp/Dy8irQ2y71fSgFnoV1+jSkpcFnn+U9edxM5crw7LOwbh0cPw7vvgsnT5ruMT8/eOklOHKk4K8jhChWbG1t2bFjB9HR0WzatIndu3cDMHHiRPbv38/mzZuJj4/n7bffBszajH+6WaslKCiIhg0bAhAeHk5UVFTOfb179wYgLCyMunXr4ufnh6OjI8HBwZw4caLA76nUt0AKzN8ftm0Du0L4UwUGmmTy1FOwZInp3nrrLXjzTXjgARg/HmrXtvzrClGW3aGlUNjKly9P+/btWbJkCfXq1ctpUTg6OvLII4/w3nvvAabFce2HfHR0NJUrV77heo6Ojjm3bW1tSUlJueE+Gxub6x5nY2NDZmZmgd9LqW+BWIS9/c3HPCzF1hZ69IDFi83g+/PPw6JFUK8eDB4MO3cW3msLIQpdXFwcCQkJAKSkpLBixQpq1aoFwKlTpwDT4pg7dy716tUDTOthxowZaK3ZsGEDHh4eN3RfWZskkOKmShXTComMhGeegfnzTf/qXXfBsmVmLEYIUaKcOnWKDh06UL9+fZo2bUrnzp3p2bMnAIMHDyYsLIywsDDOnj3LSy+9BED37t0JDg6mRo0ajBw5ks8//9yab+Gm1M362UqjJk2a6BK5oVR8PHzxBXzwgblduzZ8/TW0amXtyIQoMfbt20dt6Q6+o5v9nZRSW7XWTW72eGmBFHeenmamVkwMzJoFCQnQujXcdx/s32/t6IQQZVipTyClppiioyMMGGAWKk6YYLqz6tY1a0kOH7Z2dEKIMqjUJ5BSV0zRzQ1eecWMkTz5JPz0k1nV/sorpnUihBBFpNQnkFLLxwc+/NAUc+zaFV5/HYKCYOJEsMD0PCGEuBNJICVd1aqm1tb27dCmDbzwAgQHw7x51o5MCFHKSQIpLRo2NFN+f/nFLHq891645x7YvNnakQkhSilJIKWJUtCvn5mdNWmSKeLYurWp/XPunLWjE6LMy8rKolGjRjlrQMCUc79SjqRhw4bsuFynS2vNk08+SY0aNahfvz7btm2zVti3JAmkNHJwgP/8x9TUGjIEPvnErB/59FMZHxHCiiZPnnzT9SjvvvsuO3bsYMeOHTl1ra4t5z5lyhRGjx5d1OHekSSQ0szPz2yAtWaN2cxqzBjTIpFuLSGKXHR0NL///jsjRozI1eNvVc79WlFRUdSuXZuRI0dSt25dunTpklMLq3379jz11FO0a9eO2rVrs3nzZvr160dISEjOaveCkmKKFpKamcqR+COcSDzBvrh9JKQmsPfsXlIzU3G1dyUuOY6s7CwyszPJyM6gvFN53B3dcbV3JUtnUdG1Io62joR4heDl7IVSijo+dQh0D8TWxhYbVYBc36aNKUX/449m6m+zZjBihJm5VamS5f4IQpQA45aMY8dpy5Zzb1ipIR91vX2RxnHjxvHOO++QlJR0w30ltZy7JJB8ik+JZ1bELBYcXEBUQhRRCVGkZaXd8Dg/Nz9sbWzxcPQgS2dRya0SjnaOXEi9wIkLJ0hITSA5I5m0rDQysjLI0jfWunK1d6VhpYZU8aiCm4Mb3i7etKvajioeVXLO5cqDD5odEcePN1OAf/wR3nsPRo4s3GKRQpRxCxcuxNfXl/DwcFavXn3dfRMnTqRSpUqkp6czatQo3n77bV555ZVCK+cO5JRzlwRyBwXeUOoaGVkZzNo9i4UHF7L48GIupl+kpldNwiqG0a1GN4IqBFGtfDVaBLSgglMF7G3ztjlUZnYmkecjOZt8loTUBKITo4lJiiHuUhzbTm9jzbE1xF2KIyM7g4l/TTTvD0V45XBaBrSknEM5fF196RTciZreNbGzucl/Xnd3U1dr5EizidWjj5pdFj/8EK7Z4EaI0upOLYXCsG7dOubPn8+iRYtITU0lMTGRIUOGMHPmzBJdzr3UJxCt9QJgQZMmTUYW5DorIlcwasEojiYcxc/NjwfqPEC/2v3oHtL9jltT5padjR2hXqGEeoXe8jFXksue2D0cu3CMA2cPsCt2F99s/4bkjOScxznZOVHbuzZ1fetS06sm7aq2o65PXbxcLn/jqF0bli+Hr74yXVmtWpnurTffBBcXi7wfIYQxceJEJk40X/pWr17Ne++9x8yZMwFytqm9WTn3Tz/9lAEDBrBx48ZiWc691CeQgkpKS2LskrFM2zGN4ArB/NDvBwbUG2CxpJFX3i7eeLt4U8Pz+hZVWmYaO07vwEbZsOnkJiLPR7Inbg+ro1Yzc9fMnMeF+4XToGID2lRpQw3PGoSPeBiXwYPNrK3Jk+HXX82srctNXyFE4Ro8eDBxcXForWnYsCFffvklYMq5L1q0iBo1auDi4sK0adOsHOmNpJz7HaRnpeP5tid9avXhm97f4GTnVAjRFa4LqRdYFbWK3bG7WXpkKXvj9hKfEg+Ag60DrQNb0yqwFc1jbAh991tCd51E9expurYK2EcqRHEg5dxzJ6/l3CWB5ILW2motjsKQlZ1F5PlIDp47yKqoVSyPXM6uM7ty7q+iPWgfkUjtNHc69H6SZiPGo2xkxrcouSSB5E5eE4h0YeVCaUoeALY2toR4hRDiFUKP0B4AJGckszF6IzvP7GTdiXXMtV/MjMwLcOp1vF+ZRFP/ZjSt1ZFOwZ1oHtAcB1sHK78LIYS1SQtE3FS2zmb/qQj+nPEGG7b8xhbfLPb6KrKVxs7GjhYBLWhWuRmN/BrRt1ZfXB1crR2yELckLZDckRaIsAgbZUOdyg2o89+fefTMGXjiCRK//YWlrSqyelBLZpxewV/H/8p5fGO/xnQO7kyX6l1oW6VtnqcwCyFKHunYFndWsSLMmYP7nAX0j3Lhs0fnk3R8COljz/FN728Y0WgE2Tqbt9e9TccZHSk3sRwNvmzAC3+8wIboDWRkZVj7HQghCoF0YYm8OX8eXn4ZPv/cJJbPPzel45Ui7lIcq6JWMe/APLaf2s6+s/sAcLZzpql/U+6udjc9QnsQ7hde6saVRPEmXVi5k9cuLGmBiLypUMFU9d2wATw9Tfn4++6DuDh8XH14oO4DfN/ve/Y+sZf45+L5od8PDG80nGMJxxi/ZjxNpzbF+11vRi0YxcKDC0lMS7T2OxKiSFSrVo2wsDAaNmxIkyZXP4/Hjx+Pv79/Tjn3RYsW5dw3ceJEatSoQc2aNVm6dKk1wr4taYGI/EtLg48/hpdegvLlYerUWy5AzMrO4sC5Ayw/spyVUStZdXQVSemmqFwdnzr0qdmHu4PupmVASxmQFxZXHFog1apVY8uWLXh7e193fvz48bi5ufHss89ed37v3r0MHDiQTZs2ERMTQ6dOnTh48CC2traFFqO0QETRcXSE//s/2LoVKleGPn1MAjl27IaH2trYUsenDmNbjGXegHmcfe4sS4cs5aW2L+Fs58zb696m8/864zHJg74/9uXD9R8SlRBV9O9JiGJi3rx5DBgwAEdHR4KCgqhRowabNm264XFubm68+OKLNGjQgBYtWnDmzBnAbFQ1evRoOnToQHBwMGvWrGH48OHUrl2bYcOGWSRGmYUlCq5ePdi4ESZMgLfegrVrTctk8GC4xbclB1sHulTvQpfqXXj97teJToxmReQKFh1axF/H/2Lu/rk8vexpPJ096Rnak/51+tMyoOXVWl5C5NO4cbDDstXcadgQPrpDjUalFF26dEEpxaOPPsqoUaNy7vv000+ZMWMGTZo04f3336dChQqcPHmSFi1a5DzmSjn3f7p06RItWrTgzTff5LnnnmPq1Kk5+32cP3+elStXMn/+fHr16sW6dev4+uuvadq06XWbV+VXiU0gSqlg4EXAQ2t9v7XjyY+UFHByAq0hK8tsFhgdDfHxZnw6Lg6SkuDSJcjIMF/4nZzM7fLlzeMSE83trCzw9TXXTU0Fe3tzvaAgc62kJPO4Qhu7dnAwhRgffNDsgvjww2aAfeFC+EeT/WYC3AMY1nAYwxoOI1tns/7EepYdWcas3bP4ac9PzNg5A4C2VdoyoN4AetfsTYB7QCG9GSEsb926dVSuXJnY2Fg6d+5MrVq1aNeuHaNHj+bll19GKcXLL7/MM888w7fffpvrcu4ODg45W+SGh4ezfPnynPt69eqFUoqwsDAqVqxIWFgYAHXr1iUqKqp0JRCl1LdATyBWa13vmvNdgcmALfC11nqS1joS+JdSao51or01rc225OnpcPq0mbh08iQcPw5Hj5rzhw5BZKT5gp514xYghcLJyXzOX7pkEoxSEB4OsbGmOG/lyuDjYx7j42MSFUDdumZDw8xMCAgAu9v9q6lfH9avN4UZX33VZLB33jFl43NZDsVG2dC6SmtaV2nNhA4TSMlIYeXRlaw7sY4ZO2fwxKIneGLRE3g5e/FQ/YdoX609Xap3wdneueB/JFHq3amlUFiulGL39fWlb9++bNq0iXbt2lGxYsWcx4wcOTInGeS2nLu9vX1OYrG1tb2uTHtZK+c+HfgUmHHlhFLKFvgM6AxEA5uVUvO11nutEuE/JCXB4sWwZYtpEaxaZT6gL3dD3sDb29QnDA6Gvn3NB/nFi1CunLldrpx5vp+faVG4uZkP7JQUc31XV1NtPTXVXMvJybRC4uLA2dl8yMfEmO3QfXxMq+PUKTPefe6cSVgnT5qktmePea24OJPc7O3N+7kTLy9ISDCx+PqabUTs7EyiqVULbGxcSQx6gZDp91F/yr9xePxx+N//4IsvoEGDPP+Nne2d6RHagx6hPXitw2vsOrOL2btnszxyOZ9t/oyPNn5EOYdy3FXtLjoGdWRw2GC8XbxlqrAoNi5dukR2djblypXj0qVLLFu2jFdeeQW4Ws4d4LfffruunPugQYN4+umniYmJ4dChQzRr1sxq7+FmilUC0VqvVUpV+8fpZsDhyy0OlFKzgT6AVRJIUhIsWAC//276UQ8duvpN3d0d/P3Nt/V27aBmTZMIPD3NB62Li/mQL260vtq1lZZmksrZsyZJrF5t3sehQ+b+mBjTYjl0yCSxXbsgO9ucT03955VrAqY5XWXDcUIbHuBS4AnCewfQuo3C0dG0dnx9ISTEJLs7sbOxo7FfYxr7NQZMDa91x9fx896f+Xnvzyw8uJCnlj5F5XKVGVhvIF2qd6FDtQ6yMl5Y1ZkzZ+jbty8AmZmZDBo0iK5duwLw3HPPsWPHDpRSVKtWja+++gow3UwPPPAAderUwc7Ojs8++6xQZ2DlR7Gbxns5gSy80oWllLof6Kq1HnH5+CGgOfAq8CamZfK11nriTa41ChgFUKVKlfBjN5kdlBtam6QxfTosWmQ+ZF1cICzM9Nh07Gi+WNesWXZ3hs3IMEkkLs500cXHm+67vXtNC+nInlTO7znJ4aSKXOLmW/Da2UFoqGmdZWZCtWrmb1yjhmlt+fqaLrTbWXtsLcuOLOOPo3+wJWYLmdmZeDl7cW+te+lfpz9NKjeRgfgyqDhM4y0JSnw595skkP7APf9IIM201mPyct38rgOJiTELrTdvhkqVoH9/M07csmWuu/TFFVqT+cVULrzyPrGJTqQ9/hSRrR/iz79tOX3adKEdPw4XLsDhw6Zr75/8/U1Cr1fPJJby5U0Lr2ZN0wKsWdN03dnYmH1Q5u6fy++Hfmf+gfmkZaXhbOfMvbXu5a6qd9G/bn88nT2L/u8gipwkkNwpjcUUo4HAa44DgJjcPrmge6JXrGg+oL75BoYOvcMAsrg9pbB7fBRe/fvi9eijMPkRGq77jH7TppmMcI20NJME9u414zVnzpgWzfr1Zife6GiT1M+fv/lLVakCPj4e1K//MPd2epjGzsl41Y5g5fnpLD48m1m7ZzFm8Rh6hvakdWBrmgc0p7l/c+nqEiIPSkILxA44CHQETgKbgUFa6z15ua6sRC9mtIaff4YnnjCDLS+8YFa02+ftAzwjwySTnTtN99mFCzBrljnn7m4mEFy6dOPzGjW/yCXHg0TbrSGZs1AuhgpVT9LrrgAGhPWndZXWuDu6W+jNCmuTFkjulOguLKXULKA94A2cAV7VWn+jlOoOfISZxvut1vrNvF5bEkgxFRcHo0bB3LmmFbJwIVStarHLJyebhfGbNsHu3aYlExVlEkxWlukyu459MjifQ3kfxsfLjuAqjgzrXZ3m9b0ICTEzz0TJIwkkd0p0AikM13RhjTx0ZSqRKF6ys+HLL80SYRsbsyDx6aeLZEbClS6w/fshIgJ2RmSy60gcB4+kEXuw2nWPtXNMJzBQQ5YD1aop2rY1g/5+flC9uuk2K6uTKIo7SSC5IwnkFqQFUgJs22ZKxS9aBJ06wWefmU9oKzlzRrP+5F/8uGoPCzfu5eKxEIitC1F33/Txnp5mMWbNmmZtjJeXWaAZHn514aawDkkguSMJ5BYkgZQQWsOHH8Ibb5j+p9dfN62RYjD//VjCMaZsncIv+37lwNn9kGWPV0orqicPJMixKTZx9Ti434GICDOV+VpXpiT7+Ji1QIGB0LgxNGlSyCVmBFA8EkhCQgIjRoxg9+7dKKX49ttvadmyJePHj2fq1Kn4+PgA8NZbb9G9e3fAlHP/5ptvsLW15eOPP+aee+4p1BglgfyDdGGVUCdOwJNPmrGRpk1h2jSz1L2Y2H5qOz/u+ZEN0RtYc2wNAE52TvQM7Un3qvdTy+4eHHR5tmyBAwfMwstt28yamH9OT65UySSUoCCTVBo0MLP/wsJk1p+lFIcE8vDDD9O2bVtGjBhBeno6ycnJlC9fvkSXcy/1/zy11guABU2aNBlp7VhEHgQGwq+/wuzZMHas6Qd6+mkzPlIMvq438mtEI79GAJy5eIYpW6ew4eQG5uydw5y9c1Ao2lRpQ696vXiifz+CKwTnlFa5cMFMU/77b/juu6vlYxYvhp9+uv516taFOnWgdWvTJVajxtUCmaLkSExMZO3atUyfPh0wBRAdHBxu+5xblXNv2bLldY9zc3Nj7NixLFy4EGdnZ+bNm0fFihUZNmwYzs7O7N+/n2PHjjFt2jS+++471q9fT/PmzXNiKYhSn0BECaYUDBxolvoPHQoTJ8Jff5lFOSEh1o4uR0W3irx818sAnE0+y7Ijy1hyeAnLjizjuRXP8dyK5wj3C2ds87H0rtkbDw8PwCxQvffe6691+rSZLXbwoGm5REbCsmVmxvO1AgKgc2czK8zf3/w5unY1PX1OTkXxrkuucePGscPC9dwbNmzIR7ep0hgZGYmPjw+PPPIIO3fuJDw8nMmTJ+N6eVqflHMXhSY7OxutNVprLl26xNmzZ7Gzs8PGxoadO3eSkJBAZGQkFSpUwM/PjwsXLmBra8vp06fJzMwkMzMTFxcXkpOTsbW1JTk5GTc3N6KiooiLi6NOnTr4+vqSmJjIrl27cHR05M8//8THx4cBAwaQnJzMxo0bqV27Nv7+/qSnp/Pbb7/RpUsX7C+v23Bzc8up+FmpUiXs7OxyiscdPnyYqlWrUrt2bY4dO0aDBg1IS0sjJSUFR0dH3NzccHR0xNbW9uYFEH19TfGxSZNg/HjTxzNxIowZU+zKAXi7eDMobBCDwgYBsOnkJn7e8zO/7v+VoXOHAlCtfDVGNR5Ft5BuNKjY4Lr3XKmS+enU6eo1MzNNj97ff5vWyvHjZlry99+bhJGScn0Mfn5mYtugQeZP1aaNmSUmrCczM5Nt27bxySef0Lx5c8aOHcukSZN4/fXXpZx7cVbQleiWlpmZSVxcHFlZWURFRbFy5UpSUlI4cOAAGzZsoFKlSri6uhIfH4+DgwMxMTHExsbmNHfT/zk6awHz5s3Lue3l5cW5c+cAOHHiBNu2bcu5b/Hixdc972a7oxWUr68vFStWxMbGhvT0dOrXr49SirNnz5KZmYlb+/b4HTiA67hxOLz7LsFDh3Le3R1/f3+UUnTo0IH09HQ8PT1xcnIiISEBJyennG/9Ra2ZfzOa+Tfjjbvf4Oe9P7MxeiMro1bywsoXeGHlC9TxqUO3Gt24r/Z9NKnc5KYr4e3sTLdVUND157OyTAJJSoKVK+HPP81eXufOmZbLhx9efWzFiqY7rEEDs/iyQgXo2dPUcnN0LBa9gkXmdi2FwhIQEEBAQADNmzcH4P7772fSpEkAUs69OCvKMZD09HQOHz5MREQEp0+fZsmSJZQrV47z589z6NAhMjIyiIm5fRWWlJQUgoKCCAoKYt++ffj4+NC6dWvS0tKoV68e3t6mTPnFixdxdHQkMDCQ+vXrExMTg42NDZ6enpw7dw5/f388PT1zEs+VATsw/6gSExOJi4sjPT0dHx8f7OzsyMrKws/Pj/j4eNLT01FK4eTkxMWLF8nKysppxXh7exMfH4+Hhwdaa3bs2IGNjQ2+vr6cO3eOhIQE3NzcsLOz48SJE+zZs4cWLVpw/PjxnMcGBgZiZ2fHhQsXuHjxIpcuXeLAgQPExMQQGBiIvb09sbGxrF69mosXL2JnZ4eHhwdZWVksOXOGTDB16SfeUEPzpho2bEhcXByenp7Y2NjQvHlzEhMTqVGjBra2tsTGxuLv74+7uzt+fn6EhoZSvnx53NzcqFChAnDzb3+55WjnyJD6QxhSfwgAEWcimLJ1CquPreaD9R/w/vr3cbJzom+tvowKH0XbKm2xtbn9YOmVsdRy5cxuwn36XL0vK8skkq1bTfmX48dNWZgvvrhaNfn1168+fsgQk0wCAszvYjRfoVSoVKkSgYGBHDhwgJo1a/LHH39Qp04dQMq5lzlaayIiIpgxYwaxsbGcPHmSv//+m9Qb65nndOdkZGTQqlUrQkJC8PX1pVy5cjg4OGBjY0ObNm3IysqievXq2OSzS6Z+/fp5eryPj0/OtMF/8vK6vlpt+ZvUWXe9Zkn2XXfdlafXtoTz58+TdvIkF0ePZu9ff1Ghc2eO9OtHXFIS6enppKenk5WVxeHDh1FKcejQIWJjY0lOTiY1NZWIiAh8fX2ZPXt2rl+zVq1a+Pr64uLigpeXF3Z2dgQHm8HxjIwM/Pz8SEtLw93dnTp16hAUFHTdt8trhVUM45PunwCw8/RO09W192dm7Z7FrN2z8HHxYUj9IXQM6kj3kO55Tl62tqbnr1s383PFlS+d58/DihVmtvTevTBz5vXPt7c3M8Lq1YNGjcy+NG3bmvGVm3wJFrnwySefMHjwYNLT0wkODmbatGmAlHMvEQq6DuT06dPMnz+fDRs2MHfuXM5fU8XP3t4eV1dXmjdvTvfu3Vm6dCnPP/88tWvXxtPTUzY2Kkxaw2uvmR9bW7MT4mOP5bpP5siRI9ja2qK1Jj4+Hq01hw8fJjk5mfT0dGJiYli/fj0RERFUrlyZtLQ0YmNjSUpKIi0t7Y7Xd3V1xdvbGx8fHxITE6lUqRK1atXKGQOqWrUq4eHhpKen4+vri19NPxbvWcwvMb+w68wu0rLSCHQPpHNwZwaFDaKxX2MqOFco6F/tBsnJpsTLsWOm1RIRYQbyjx27sWBlaKgZp2nSxFSlbtTIrHMpZp9t1ykO03hLAlkHcgv5TSCnT59m+PDh1/X/29nZ0axZMx588EFGjBiBi4uLJUMV+fHnn/DAA2YaU48eZgCgEGdqZWdn5/xOTk5m+fLlVK1alaSkJFasWIG3tzfZ2dkcPHiQ2NhYTp06xeHDh/H09CS365EqVqqIq68rqc6pxGTHQDLYONlQr3E9ghyC6Hd3P2oF1CI0NPzq30oAACAASURBVBR3d/d8t15vR2tTN2z/fpNYEhNN4cqVK6/fQMzT00wxrl3bbKZWvbpJLO7FpB6lJJDcKZIEopRyBVK11kW0m3f+FXQhYVpaGm3btqVVq1Z0796dTp06oZSSVkVxlJVlWiD//a/pc3nzTRg92tpR3eDKrLpz586xdetWYmNj0Vrj6enJ2rVrc2bPRURE4OLiQmRkJPHx8be9ZsWKFalUqVJOSycgIABvb2+CgoKoUqUKYWFh+Pn54erqiq2t7R3XINxJZqZJIBs2wLp1sH272db5n7NMmzUzg//u7qa10ru3WXlf1K0VSSC5UygJRCllAwwABgNNgTTAEYgDFgFTtNbFepl3QbqwtNaSMEqS3bth8GCz3+5DD8E775g+lxIsIyODlJQUVq9ejbuXO0s2L2HDsQ1s3LPRjL1FgyuueDl7cerwKTLSM3ImPdyMr68vlSpVwtfXl8DAQMqXL4+NjQ3Vq1enUaNGVKxYETc3N9zc3HDOwz7Mp06ZFsuhQ2aM5eBBUwH58OGrj3F1NY3DevXg7ruhRQuoVatwZ4Lt27ePWrVqyf/Ht6G1Zv/+/YWSQNYAK4B5wG6tdfbl855AB2AQ8JvWeuatr2JdUgurjMnOhqeeMtOOXF1NSRQrDPYXtoysDFZFreLXfb+yPHI5kecjIQsaVG5Azxo96e3bm4ykDKKjozl+/DgJCQk5s+/27NnD6dOnOXny5A3TP69wcnLC1dWVmjVr0rZtW2xtbUlLSyM0NJSkpCS6du2Kt7c3FSpUuG2rJjXVbAA2b55ppezcaZJNQoK539nZdIMpZWaT+fjA8OGmIIElHD16lHLlyuHl5SVJ5CautIiTkpII+sd8cUskEHutdUZBH2NNkkDKqIUL4b77THXDRx6BKVNKbYEprTVrj61lxs4ZbI7ZTERsBE52TjSq1IiB9QbSuXpnannXuuF5Vwb0d+3aRUxMDEopFi9ejJubG/Hx8axcuZKsrCzOnj17y9cuV64cTZo0wd3dndDQULy8vMjMzKRFixZUqlQJPz8/srOz8fb2znlORoapD7Zhg2k07t1rFktef10z66tuXdNaGTDAlHTJq4wMk0RvNlNSGE5OTgQEBOQsDr5CBtGRBFKmJSaaVeszZpjjr76CESOK3Sp2S4s4E8Fnmz9j4cGFnEwygxPdQ7rTrko7OgZ3pLFfY2xU7v4GWmsuXryIi4sL+/fvZ9euXWzfvh0/Pz/mzJkDwNmzZzl48OBtr1OuXDn8/f3x8PCgS5cuhIaGcvLkSYYMGUJycjLVq9dg+3bF+vWmhRIbawbv9+y5vsJx9epmTCU4GO6/3xSeFIXDoglEKfUfrfXbFomsCEg1XgFcne77+efmU2nAALOJlZVWqBclrTW7Y3czZesUlkcu58C5AwCE+YbRIqAFj4Y/SnjlcIu9XkpKCikpKcTHxxMREcGxY8fIyMhg+/bt/P7772RnZ3PxnyWJr+Hu7p6zmDYkJITq1atTvXoNzp71Y9GiS2ze3ITjxxXXDu/Y25vV9QMGmGnG9eubIpTWGLAvbQqUQJRS19YHVUBDrXXxqWSXS9ICEYCZqTVmjBkb8faG6dOhe/cyVcsjOjGa73Z8x5dbvyQ6MRqAUK9Qwv3CebHti9TxqVPo4wSZmZmcPn2a06dPc+bMGXbt2kV8fDxnz55l27ZtORUNUlJScqZMX8vPrzLBwbVYt24lbds+S0BAH5Yvd+LsWTegJuajyhgyxLRWHBxMb2YFyy+jKdUKmkC+1lqPuOb4C6118ZsbeQeSQMR1fv0VHn0Uzp41ZW3nzy9zZWy11qyKWsXPe35my6ktbIkx/39UcKrA4LDB9KrZizZV2uBib711TtnZ2Zw8eZL169fz66+/5qyzqVat2i27y5ycXKhYsRYQRExMIFlZGWRntwTCgEBatChPq1aKli3NYkiloGrVonxXJUtBE0iQ1vroNceeWuvbT0ovhiSBiBtcuGAG1Z97zkz7mTr1+oJSZczh+MPM2TuHDzd8yNnks2SbyZaMaTaGcS3GEVQ+qNjNYLpSS+3PP//kf//7Hy4uLpw9e5bU1FR27959064ypVzQuhMQCiRha+tDixYdOXOmIsOH1+KBB5RUL76GRcZAlFLeWutbT8Mo5iSBiFuaP98MqsfFmXLxL79c6gfY7yQpLYlpO6bx9NKnybq8Xrha+Wr0DOnJfXXuo13VdrkegLeW7Oxs4uLiSExMZNu2baxdu5b09HQOHjzIxo0bb1OKxh0Hh2CCg9sSEOBFp06VefTR/pQvX75MrgmzVAKZr7XubdHIipAkEHFbhw/D//2fWS9Sp45Z0X7tphxl2O7Y3Sw7sowVkStYemRpTstkaIOh9KvVj641uuJo53iHqxQv2dnZpKamopRi7dq1XLx4kY0bT7Bo0fekpnoQF3eexMRdYOo+o5QzdnblyMiIxd3dm8aNG1CunAtPPvkkLVq0wM3NzbpvqBBZKoEs0Fr3smhkRUgSiLgjrc0U31dfhUuXzIythx4qUwPsdxKfEs+H6z/k400fk5iWmHN+QvsJPNH0Cbxc8rFIo5g6deosmzYd57ffolm/fgUxMalcvHgJWAmcvuXzmjZtyrPPPkurVq0ICAgosngLi7RAkAQi8iAmBvr3N6vaOnQwSaUYbaFbXBw6d4gP1n/Al1u/BEChqOJRhcFhgxkZPpJq5atZN0ALy8q6Wv/rt9+SWLZsA/HxaZw7tw74CYi84Tn+/v5UqFCB3bt306hRI3r16kW/fv2oW7cuNjY2hVIA09LKdAtE1oGIfMnMNFvovv+++eR45RUYO9YsOBDX0Vrz94m/WRG5gvXR61keuZxsnU2IZwjdQ7oztvlYAj0CsbMpfRUAMjJMIehvvjEr5OfNS+P48T8AhZPTZjw8juDmdp4jRxbc8Nw6derQuXNnfHx8aNeuHY0bN75un53iwlIJpJ7WerdFIytC0gIR+XL8uNlcfN066NoVpk0r8YUZC9vxC8f5dvu3/LrvVyJiIwCo4lGF3qG9ebzp49T2Kb1VcbU2+6rMnAmrVl0tgQ8ZODhcwsdnLqGh0URGzubSpTOkpqbmzBSzt7enVq1aRERE0KhRIxo2bMjjjz9O48aNrdpSkVImSAIRBfTee2aQvXx5swhxwABrR1QizNs/j6+3f010YjQ7Tu8AoEHFBrSt0panWz5NUIWgO1yh5IuPN91en3xiekUTLw8dXanxVa7cCpo2Pc25cxEsW7aEXbt23XCNDh060KJFCwC6d+9OmzZtiix+iyUQpVQT4EWgKmY7XAVorXXe9lO1AkkgosD+/hseftjM2Ore3awbkf1dc+3o+aNM2zGNDdEbWB21miydRVD5IJr5N2NQ2CC6h3Qv9lODC0pr0zKJjDSViVevhmuXqri5wd13n+f557NYuvQztNZMmDDhhuvY2dkxcOBAWrduTY8ePQp1sN6SCeQA8H9ABJBTX0BrfaygQRY2SSDCItLSYMIE+OgjMx7y1VfSGsmHk4kn+WjDR8w9MJfYS7EkpiUS5htGq8BWPNn8Ser41LF2iEXmt9/MP6e1a68/37Chqes1bBi0b69JTU3grbfeYufOnXh5eTF79uycx7q5uVGpUiWaN2/O22+/jb+/v8Xis2QC+UtrXXRtJwuSBCIsatcus/f6+vWmTPyXX5piSyLPjp4/ytglYzlw7gAHz5nyJAHuAbQObM27nd8l0MNCm4KUADt2mP3o160zxaNTUq6//7PPYOhQ01JJTk5my5YtzJkzh/nz53Ps2PXf4z09Pfn444/p2bMnHgUoGmrJBNIRGAj8gdmVEACt9a/5jq6ISAIRFpeYaGZnTZ4Mfn7wyy+map/It8Pxh/lk4ycsOLiAowlHsVW2dAruRIhnCGOajyHUK9TaIRapFSvMljbffHN9V9cjj0Dfvqacfd265lxKSgrvvfce+/fv54cffrjuOs8++yxvv/12vgbjLZlAZgK1gD1c7cLSWuvheY6qiEkCEYXml1/M/9EZGWaq7/PPl4ky8YUpKzuL7ae3M3XrVL7b+R1pWVfLjnSr0Y3/tP4P7aq2K1NlRWJi4McfzQ7Np69ZxzhunNkPZdiw6yvwpKSk8NtvvzFmzBiaNm3KkiVL8vW6lkwgEVrrErl1iyQQUaiOHIEXX4SffzZ7s37/PXTpYu2oSoXUzFR+3P0jb/75Jofir67lCvUKZVC9QTzc8OFSt2jxTi5ehCVLTAN43z5zztkZHn/czO/o0OH6AgpZWVnY5nNjFEsmkKnAh1rrvfmKxIokgYgisWwZ3HOPuf3qq/DMM2ZfVmERaZlp/H3ibwb9OojTF83XcHsbe+6tdS+9a/ZmQL0BpXLB4u1cuGB2J5g/35RyA7Nr84wZMHBgwa9vyQSyD6gOHMWMgRT7abyyEl0UuRMnzF4jixeDi4uZ/tuggbWjKnVSM1P589iffLzpY3ad2cXxC8cB6BzcmeGNhtOnZh+c7Z2tHGXRiomB0aNNMrnio4/g3//O/86MlkwgN912RabxCvEPWpv/i++91xx//LHZCVEUimydzdz9c5mxcwbLjiwjJTOFuj51Gdl4JM38m9HMvxm2NmVnb9vt26Fx46vHjRrBtm35u9btEkheh+RfAy5orY9dThqJwKv5C0uIUkwpsznVmjVmqsyTT5r1IqdOWTuyUslG2dCvdj/mDpjL4ScP83CDh0nJTGHc0nG0+rYVdq/b8ebaN3O6vUq7Ro0gO9uUVRkyBO6/v3BeJ68tkO1a60Z3OlccSQtEWE1GBrz1FkycCO7uZipNhw7WjqrU01pzOP4wDb9qSHJGMmDGS9pXa89jTR6jfbX2eDp7WjnK4s+SLRAbpVTOlvRKKU9MSRMhxK3Y25sB9e3bTcnWjh2haVOzakwUGqUUIV4hXHrhEvpVza7HdjG2+Vg2x2zmvp/uw+99P9QExbPLnqWs1AS0tLy2QIYCzwNzAA08ALyptf5f4YRnOdICEcVCUpJZKzJtmjl+5x1TpFEUmZSMFDad3MTUbVP5PuJ7AIIrBNMrtBc9Q3vSvlr7MjeT63YsWo1XKVUHuBszA+uPkjKlVxKIKFaOHIGePWH/frPr4RdfQDHcC6K0O33xNHP3z2XBwQX8EflHzoLF4Q2H81K7l8pEteA7KXACUUopfYcH5uYx1iQJRBQ7KSkwfDjMng2+vqZeRc+e1o6qzDqbfJYH5zzIyqMrc8652Lvw4T0fMqzhMBxsy2atM0uMgaxSSo1RSlX5x4UdlFJ3K6W+Ax4uaKBClCnOzjBrlil2lJoKvXpBnToQF2ftyMokbxdv/hj6B1Fjo/h90O80829GckYyjy58FMc3HLnvp/v4dV+xL/tXpHLbAnEChgODgSAgAXDGJKBlwGda62I9IigtEFGs7dkDbdvC+fOmhveyZVD1psuuRBHJ1tmkZaYxddtUFh1axJpja0jNTKWWdy0eCzezuOr61i314yWWHgOxB7yBFK11ggXiKxKSQESxp7VZNvzCC2bqb9++ZrDdzc3akQkgMzuTqVun8tXWr9h5ZmfO+YkdJ/JMy2ewt7W3YnSFR7a0RRKIKEGOHzer1ufPh3r14IMPoFOn66vjCavRWvPLvl/o/3P/nHO+rr70Du3Nf9r8hxqeNawYneVJAkESiChhtIYffoD//heio82G2v/+t7WjEv+QmZ3JrIhZzN4zm0WHFgHQKrAV91S/h9FNRuPj6mPlCAtOEgiSQEQJFR8P4eGmJkXt2iaJPP64taMSN7ElZguTN05m5q6ZOee61ehGy4CWPNvq2RJb2LFUJhCllCvwOZAOrNZaf3+7x0sCESVWcrIpxvjOO2aQ/eOPTSKRLq1iKTkjma0xW2k3vV3OuRqeNehWoxuTOk3Cxd7FitHlnSWr8ToC9wHVuKaEidb6tQLGeOX63wI9gVitdb1rzncFJgO2wNda60lKqYeABK31AqXUj1rrB293bUkgosSLjITOnc1vFxdTCiUkxNpRiVvI1tkciT/C5I2T+WzzZwA42znzUruXeKLpE3g4lYxdKy1ZC2se0AfIBC5d82Mp04Gu155QStkCnwHdgDrAwMur4QOAE5cflmXBGIQonoKDzfZzb75pWiUdO8Jff0GW/PMvjmyUDSFeIXza/VOOjzvO9/2+p02VNry48kXKv10e17dc+ev4X9YOs0Dy2gLZfW3LoDAopaoBC6+8jlKqJTBea33P5ePnLz80GjivtV6olJqttR5wk2uNAkYBVKlSJfzYsWK/bYkQufP55zBhAsTGmnUjP/4IDRtaOypxB1pr/jj6B5P+msQfR/8AwMfFh5n9ZtKlevHcAtmSLZC/lVJFvSe6P1dbGmAShz/wK3CfUuoLYMHNnqi1nqK1bqK1buLjU/JnQwiR4/HHTT2tgQPh4EFTAmXpUjN7SxRbSik6BXdixdAVnHn2DD1CemCjbLhn5j2oCYpRC0axN65ElBcE8p5A2gBblVIHlFK7lFIRSqldhRHYNW42Uqi11pe01o9orUffaQBdiFLJzc1M9V2+3CSOrl2hVSs4cMDakYlc8HX1ZeGghUSOjWREoxEEuAcwddtU6n5el/bT27Py6EpSM1OtHeZt5XUNfrdCieL2ooHAa44DgJjcPvmaPdEtHZcQxUOnTnD0KHz1ldn5sFYtWLfOJBNR7LnYuzC191QADpw9wOjfR7MqahUdZ3SknEM5Pu/xOUPqD7FylDeXpxbI5W1sywO9Lv+UL4L90DcDIUqpIKWUAzAAmH+H5+TQWi/QWo/y8CgZMx6EyBcHB7N6/ZdfzHGbNtC/v1nVLkqMmt41+WPoH5x46gQfd/2YpPQkHvrtIYInB/Pj7h+tHd4N8pRAlFJjge8B38s/M5VSYywVjFJqFrAeqKmUilZK/UtrnQn8G1gK7AN+0lrvsdRrClGq9OsHhw7Bww/DnDlQvz6sXWvtqEQeKKUIcA9gTPMxJP43kSebPcnRhKMM+GUAdT+vy5aY4rMcIa+zsHYBLbXWly4fuwLrtdb1Cym+ArumC2vkoUOHrB2OEEVnzRro3RsSE83v334Dm7wOe4riYO2xtTz+++Psibv63fnLHl8yKnwUqpAXlFpyFpbi+jUXWdx8kLvYkC4sUWbddZeZqRUYaAoz9u1rZmyJEqdd1XbsGr2L4+OOMzhsMACP/f4YL/zxAqujVlttT/e8JpBpwEal1Hil1HhgA/CNxaMSQliGtzccOwZvvw0rV0LNmqab6+JFa0cm8shG2RDoEcjMfjM58dQJAt0DmbRuEh2+60CvWb24lG7JNd25k5/9QMKB1piWx1qt9fbCCMxSpAtLiMvOnDHrRbZsMSVQli+XTatKMK01G6I3MHTuUA7HH6aqR1W+6PEF3UIsO1m2VBZTzCuphSUEkJQEP/9sZmwpZbbU7dXL2lGJAkjLTGPJ4SW8sPIF9sbtZUj9IUzrM81iOyUWeAxEKfXX5d9JSqnEa36SlFKJFolSCFH4ypWD4cMhIsKUh+/d2xzHx1s7MpFPjnaO9KnVh80jN9OgYgNm7pqJ0xtOLD+yvNDHRnKVQLTWbS7/Lqe1dr/mp5zW2r1QIxRCWF5wMKxYAXffbbbN9fKCTz+1dlSiAFzsXdjx2A6+u/c7snQWXWZ24cnFT5KtswvtNfO6DuTt3JwrTpRSvZRSUy5cuGDtUIQoXjw84I8/YMoUczxmDDz6KKSlWTcuUSBDGwxlzbA1AHy6+VP+Nf9fhfZaeZ2F1fkm56xR3iTXZBqvEHcwcqTZY6R/f5NMqlc3035FidWuajuyX8nmyWZPMn3HdMYtGVco3Vm5HQMZrZSKwKwQ33VNIcWjQGEXUxRCFLagIPjpJ/j9d3Pcpw8MG2YG3UWJpJRifPvx9K/Tn8jzkYWy4DBXs7CUUh5ABWAi8N9r7krSWpeI0TeZhSVELp04YVojGzea48WLTaVfUWJprfOdQAo8C0trfUFrHQUMx2xp+yHwAfCwUsopX1EVERkDESKPAgNhwwaYONEcd+sGzz4LJ09aNy6Rb4VV7iSvtbB+ApKAmZdPDQQqaK37F0JsFiUtECHy4fhx0/rYt88cjx8Pr75q1ZBE0bpdCySvK01qaq0bXHO8Sim1M/+hCSGKtSpVYM8es/jwwQdNArGzgxdftHZkohjI6yys7UqpFlcOlFLNgXWWDUkIUawoBQ88AHsvb7X60kumBMoumT9T1uU1gTTH7IsepZSKwuzdcVcRbW0rhLCm2rXNrKzBgyE6Gpo3h9dfh+zCW6gmirdcj4EoMwrTFrjlDoRFsDthvskYiBAWdOAAjB4Nq1aZAo0zZkCFCtaOShQCi+wHok2m+VBrfexWPxaL2IJkFpYQhaBmTVi2DO6/HxYuBE9PKYVSBuW1C2uDUqppoURSSGQluhCFxM7ODK5fmZU1ZgxUqyZjI2VIXhNIB0wSOXLNanT51yJEWTZ+PKSkmH3Yjx0zOyEuWgRlZKuIsiyvCaQbEAzcDfQCel7+LYQoy5ycYPp0OHoUKlaEHj2gZUuQTdxKtTwlkMvjHOUxSaMXUL64jn0IIaygWjXYuRM+/9yUQgkNhQ8+kNZIKZXXcu5jge8B38s/M5VSYwojMCFECeXoaGZoXSkT/8wzZjX7paLfs1sUrrx2Yf0LaK61fkVr/QrQAhhp+bCEECXeyJGQnAzPPWdmbLm5werV1o5KWFBeE4gCsq45zrp8rtiSabxCWJGzM7z9tunSAujQwVT6TZSdsEuDvCaQacBGpdR4pdQEYAPwjeXDshyZxitEMTB6NBw+bG7PmWM2rfrjD+vGJAosr4PoHwCPAPHAOWCY1vqjwghMCFHKVK9udj4cPNhsm9upk9lCNzXV2pGJfMrrIHp/4JDW+mPAA3hFKdWoUCITQpQ+QUEwcybExJiB9SlToFkzWcVeQuW1C+tlrXWSUqoNZn/074AvLR+WEKJUc3Mz2+d++CFERJhV7HfdBRcvWjsykQd5TSBXBtB7AF9qrecBDpYNSQhRJtjYwLhxZt2Iry+sXWt2Q5QJLyVGXhPISaXUV8ADwCKllGM+riGEEFfVr292Phw1ChISoHx5WL7c2lGJXMjrh/8DwFKgq9Y6AfAE/s/iUQkhyhZHR/jqq6tjIV26mC6tefOsG5e4rTztiV6SyX4gQpQQZ89C27awf7857tMHZs829bZEkbPIfiCXL+SklHpaKfWrUuoXpdRTSin5ryqEsBxvbzMusn69GWyfN88sSJQvgMVOXruwZgB1gU+AT4HawP8sHZQlyUp0IUogBwdo0cIMqPv5mXNNm5pFiCkp1o1N5MhrAqmptf6X1nrV5Z9RQGhhBGYpshJdiBLMxgY2b4YRI8xx//7mdkaGdeMSQN4TyHalVIsrB0qp5sA6y4YkhBDX8Pc3Cw43bYI6deCHH0wLZcMGa0dW5uUqgVyz82Bz4G+lVJRSKgpYD7QrxPiEEAKUMl1Yu3dDu8sfOS1bQo0a0hqxoty2QK7sPNgVCALuAtpjksfxQolMCCH+SSmzRmT2bHN85IhpjQwfLptWWUGuEojW+tiVH6AC8ASwGngdWFR44QkhxD84OMCDD5qS8IMGmXPTpsHYsdaNqwzKbRdWqFLqFaXUPszsqxOYNSQdtNZSBU0IUfTKlYPvv4fYWLj3XvjkE6hd25REEUUit11Y+4GOQC+tdRut9Sdcv7GUEEJYh48PfPst3H23WXx4111mB0RR6HKbQO4DTgOrlFJTlVIdKeY7EQohypAKFcwGVYsXm+N77oEJE2RcpJDldgzkN631g0AtzNjHU0BFpdQXSqkuhRifEELkXteuJonUqAHjx0O9enDokLWjKrXyuiPhJa3191rrnkAAsAP4b6FEJoQQ+dG1qymF0r497N0LdevCq6+aSr/CovJdil1rHa+1/kprfbclAxJCiAJzcYFVq8yGVdnZ8NprppvrzTetHVmpInt5CCFKr3HjzJ7rzz5rjl96yawZyZI5QJYgCUQIUbrZ2cG775q1ImB++/ub/UekyGqBlNgEopQKVkp9o5SaY+1YhBAlwLBhpjtr4kQ4cwYee8xU+pWZWvlmlQSilPpWKRWrlNr9j/NdlVIHlFKHlVK3HZzXWkdqrf9VuJEKIUoVpeC//4W5c81xSgo0ayZJJJ+s1QKZjqmrlUMpZQt8BnQD6gADlVJ1lFJhSqmF//jxLfqQhRClRp8+pjUybJjZqMrDA86ds3ZUJY5VEojWei0Q/4/TzYDDl1sW6cBsoI/WOkJr3fMfP7G5eR2l1Cil1Bal1Ja4uDgLvwshRImmlFnB/sorkJQEDRrAyy9DZqa1IysxitMYiD+mxtYV0ZfP3ZRSyksp9SXQSCn1/M0eo7WeorVuorVu4uPjY9lohRAln1Jmxfrvv8PJk/DGG9CoEcyYYWZvidsqTgnkZqVRbtkxqbU+p7V+TGtdXWs9sRDjEkKUdt27m4WGnTqZPUceftjsw37woLUjK9aKUwKJBgKvOQ4AYgp6UdkTXQiRKx4eZq+Rr76Cxo3NuZo1TUIRN1WcEshmIEQpFaSUcgAGAPMLelHZE10IkSejRsHWrWYHRICwMNPVlZJi3biKIWtN452F2Q63plIqWin1L611JvBvYCmwD/hJa73HGvEJIQTz5sFnn109dnGBo0etF08xpHQpn/+slOoF9KpRo8bIQ1KVUwiRV4cOQXi4mallawsjRpjBdm9va0dWJJRSW7XWTW52X3HqwioU0oUlhCiQkBBT8uSvv0wNra++giZNIDLS2pFZXalPIEIIUWBKQevWsHGjOT52DKpXhwULrBuXlZX6BCKzsIQQFtOs2dVdDwF69zYbV5VRpT6BSBeWEMKiunY1tbPefdccT5hgNqwqg9UuSn0CEUKIQvHss7Bpk7n92mvg6wuzZ1s3piImCUQIIfKraVMzuH7FwIHw5ZeQlma9mIpQqU8gBxJvygAACsJJREFUMgYihChUrVubLq3du83+IqNHQ+fOkJxc6svEl/oEImMgQogiUbeumZV1113w55/g6nq1bHwpVeoTiBBCFJnwcFi9Gj74wBwvWACtWpnFiKWQJBAhhLC0p56C+Hh48EGzdiQ0FGbNsnZUFicJRAghCkOFCmZW1tNPm+NBg+C776wbk4WV+gQig+hCCKt64w344gtze9gwuP9+M8BeCpT6BCKD6EIIq3J2hsceu9oS+eUXM8D+44/WjcsCSn0CEUKIYuH99+HcOXBzM8cDBpgaWyV4qq8kECGEKCqennD+PAwffvWcjQ0sXWq9mApAEogQQhQlOzv45hsYOfLqua5d4euvrRdTPpX6BCKD6EKIYmnKFLO/yBUjR8LkySVq4WGpTyAyiC6EKLZsbMyMLD8/czxuHAQElJjpvqU+gQghRLHm7AwxMZCQYI5PnTLTfSdPtmpYuSEJRAghigMPDzhy5OrxuHFmum8x3mdEEogQQhQXwcGQng7/+Y85Tk6GwEDIyLBuXLcgCUQIIYoTe3uYNAlOnjTHaWng4AATJ1o3rpuQBCKEEMVR5crw5JNXj194odhtVFXqE4hM4xVClFiTJ8OcOVcXHgYFwfz5kJlp3bguU7oEL6PPiyZNmugtW7ZYOwwhhMi77Gy45x5YseLquawsMw24kCmltmqtm9zsvlLfAhFCiBLPxgaWL4e+fa+eGzbMauFcIQlECCFKih9/hD17oEMH+N//4N13rRqOnVVfXQghRO7Z20OdOmajKj8/eO45Mx4yYIAZHyli0gIRQoiSxtcXjh83v194wawfWbu2yMOQBCKEECWRvz/8/bfZOhfgrruKfG8RSSBCCFFSVa8OUVFXjxcvLtKXlwQihBAlmbs7REaa2z16XC2DUgQkgQghREkXFAT9+5vb77wDEyYUyctKAhFCiNLgp59g1Spze/x46N690IswlvoEIqVMhBBlRvv28Ouv5vbixbBrV6G+XKlPILIjoRCiTOnbF3buNLe3bi3Ulyr1CUQIIcqcsDDz+9FHoVu3QuvKkgQihBCljVJXa2UtWQL/93+F8jKSQIQQojT69lto0MDcXrKkUF5CamEJIURppBTs2AHvvw+nThXKS0gCEUKI0uyZZwrt0tKFJYQQIl8kgQghhPj/9u41Rq6yjuP499cu9IaRKjQKqNuW9gVga4mtRgRJjQUNEaNNJGBL0KjRFCVYsQWN4gtT0URfYFKRN5igRTFNGm8LxQS1Elpa2oUFC9uLWCEpUIWUIpft3xfPs3EYd8PMeeay6/w+yWTPnMtz/v+dmfPf58zZ51TiAmJmZpW4gJiZWSUuIGZmVokLiJmZVeICYmZmlbiAmJlZJYoO30O3WyQ9DfytoIlTgGdaFM5k4Zx7Q6/l3Gv5QlnO74iIU8da0DMFpJSkByLi3d2Oo5Occ2/otZx7LV9oX84+hWVmZpW4gJiZWSUuII27pdsBdIFz7g29lnOv5QttytnfgZiZWSXugZiZWSUuIGZmVknPFxBJF0vaK2lY0roxlk+TdEdefr+k/ppl6/P8vZIu6mTcJarmLOlDknZKeij/XN7p2KsqeZ3z8rdLOippbadiLlX43l4k6T5JQ/n1nt7J2KsqeG+fIOm2nOujktZ3OvaqGsj5Akm7JL0qaWXdsislPZ4fVza984jo2QcwFdgHzANOBPYAZ9Wt80VgY56+DLgjT5+V158GzM3tTO12Tm3OeQlwWp4+B/hHt/Npd841y38F/BJY2+18OvA69wGDwOL8/M098N6+HNiUp2cCB4H+bufUopz7gUXAT4GVNfPfBOzPP2fn6dnN7L/XeyDLgOGI2B8RLwObgEvr1rkUuC1P3wl8UJLy/E0R8VJEHACGc3sTXeWcI+LBiHgyzx8Cpkua1pGoy5S8zkj6GOnDNdSheFuhJOcVwGBE7AGIiGcjYqRDcZcoyTmAWZL6gBnAy8DznQm7yOvmHBEHI2IQOF637UXA3RFxJCL+CdwNXNzMznu9gJwO/L3m+aE8b8x1IuJV4DnSX2SNbDsRleRc6xPAgxHxUpvibKXKOUuaBXwNuLEDcbZSyeu8EAhJA/nUx3UdiLcVSnK+E3gBeAp4Avh+RBxpd8AtUHIcKj6G9TWz8v8hjTGv/rrm8dZpZNuJqCTntFA6G/gu6S/VyaAk5xuBH0TE0dwhmSxKcu4D3g8sBY4B90jaGRH3tDbElivJeRkwApxGOp3zJ0lbI2J/a0NsuZLjUPExrNd7IIeAt9U8PwN4crx1cvf2jcCRBrediEpyRtIZwGZgdUTsa3u0rVGS83uAmyQdBK4Brpe0pt0Bt0Dpe/veiHgmIo4BvwXObXvE5Upyvhz4fUS8EhGHgW3AZBgvq+Q4VHwM6/UCsgNYIGmupBNJX6ptqVtnCzB6dcJK4A+RvoHaAlyWr+qYCywAtnco7hKVc5Z0MvAbYH1EbOtYxOUq5xwR50dEf0T0Az8EvhMRN3cq8AIl7+0BYJGkmfkg+wHgkQ7FXaIk5yeA5UpmAe8F/tqhuEs0kvN4BoAVkmZLmk06ozDQ1N67fRVBtx/AR4DHSFcy3JDnfRv4aJ6eTrr6ZphUIObVbHtD3m4v8OFu59LunIGvk84T7655zOl2Pu1+nWva+BaT5Cqs0pyBT5EuGngYuKnbubQ7Z+CkPH+IVCy/2u1cWpjzUlJv4wXgWWCoZttP59/FMHBVs/v2UCZmZlZJr5/CMjOzilxAzMysEhcQMzOrxAXEzMwqcQExM7NKXEDMxiFpRNLuPCLtHknXSqr8mZF0fc10v6SHG9zuGkmrq+63pp01kq4qbcdslC/jNRuHpKMRcVKengP8DNgWEd9sQXv9wK8j4pzX2aYP2AWcG2nspsokzSTFv6SkHbNR7oGYNSDS8BafA9bk/1aeKul7knZIGpT0eQBJF0r6o6TNkh6RtFHSFEkbgBm5R3N7bnaqpJ/kHs5dkmaMsevlwK7R4iHpTElbc49ol6T5eZ/3SvqFpMckbZB0haTt+f4W83MOx4CDkibDqNE2CbiAmDUo0sB6U4A5wGeA5yJiKek/fT+bh7SBNDDfV4B3AvOBj0fEOuDFiHhXRFyR11sA/Cgizgb+RRrhuN55wM6a57fnbRYD7yONHguwGPhy3ucqYGFELANuBa6u2f4B4PyKvwKz13ABMWvO6AimK4DVknYD95OGBF+Ql22PdH+GEeDnpJFtx3IgInbn6Z2kG//UeyvwNICkNwCnR8RmgIj4d+5VAOyIiKciDa+/D7grz3+ort3DpBFnzYr1+nDuZg2TNI805PdhUiG5OiIG6ta5kP8dEnu8Lxpr76UyQrqRUb0XSeM3wdjDb4/V1vGa58d57ed8em7TrJh7IGYNkHQqsBG4OdKVJwPAFySdkJcvzKO4AizLo6NOAT4J/DnPf2V0/SY8CpwJEBHPA4fyHRJH7+89s8n2FpIGSDQr5gJiNr7RL72HgK2k00Kjdya8lTRq6658Oe6P+e9f+vcBG0gH6gOk+6cA3AIM1nyJ3ojfARfUPF8FfEnSIPAX4C1N5nRezsWsmC/jNWuhfAprbURc0sI2NwPXRcTjhe0sAa6NiFWticx6nXsgZhPfOtKX6aVOAb7RgnbMAPdAzMysIvdAzMysEhcQMzOrxAXEzMwqcQExM7NKXEDMzKyS/wDnj2gyZOMd+QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xUZfbH8c9JKEEBAY2KFImKC1iIEnBta0dABHQtICroKjbsivhTV2XVtewu6oq62FBEEQuCZe2ydiWsoICidKKuRqyI9PP747mBIU7CBDK5meT7fr3mNXPbzLkZyMlzn+eex9wdERGR0rLiDkBERKonJQgREUlKCUJERJJSghARkaSUIEREJKk6cQdQWbbaaitv06ZN3GGIiGSUKVOmfOvuucm21ZgE0aZNGwoLC+MOQ0Qko5jZgrK26RKTiIgkpQQhIiJJKUGIiEhSNaYPQkRqp5UrV1JUVMSyZcviDqVay8nJoWXLltStWzflY5QgRCSjFRUV0ahRI9q0aYOZxR1OteTuLF68mKKiIvLy8lI+TpeYRCSjLVu2jC233FLJoRxmxpZbblnhVpYSxJgx0KYNZGWF5zFj4o5IRCpIyWHDNuZnVLsvMY0Zw8rTz+ajX3diGS1YtiCH5X96kmWTW7Gsyx9YvhyWLVv3WL4csrOhfn3IyQmPZK/r14cttoBttoGmTUPuERHJNLU7QVxxBd//Wp8Cpqxbtxy4rfI+ok4dyM2FrbcOCaP0c4sWsNNO0KpVSD4iknmKioo455xzmDlzJmvWrKFnz57ccsst1KtXb739vvzyS8477zyeeOKJct+vR48ePPLIIzRp0qTCsVxzzTU0bNiQSy65pMLHlla7E8TChTShDs/Qk/osJ4dl0fNycmZN+03roF49WLMmtCRKWheJrYzE1z/8AN98A19/vf7zZ5+F17/+un4o9evDDjtA27a/fbRooVaISKUZMwauuAIWLoTWreH666F//41+O3fn6KOP5qyzzmLChAmsXr2aQYMGccUVV3DLLbes3W/VqlVst912G0wOAM8///xGx1OZaneCaN2aegsW0JPn1l+//fawc/JDsrJCq2DzzTfto5csCQlj4UL4/PP1Hy++GJJNiZwc2Hln6NRp3aNjR2jQYNNiEKl1xoyBQYNg6dKwvGBBWIaNThKvvfYaOTk5nHLKKQBkZ2czfPhw8vLyyMvL4/XXX2fZsmX88ssv3H///fTs2ZPp06ezdOlSBg4cyKeffkr79u2ZP38+I0aMoKCgYG3poCVLltC9e3f2228/3nnnHVq0aMGECRNo0KAB99xzDyNHjmTFihXstNNOjB49ms0226wyfkpr1e4Ecf316/9jAdhss7A+zRo2DI8ddoADD1x/25o1UFS0ftKYMQOeeQYeeCDsk50Nu+yipCFSIVdcsf7/dwjLV1yx0QlixowZdOrUab11jRs3pnXr1qxatYp3332Xjz76iGbNmjF//vy1+9x55500bdqUjz76iOnTp5Ofn5/0/T///HMeffRR7rnnHo477jiefPJJTjzxRI4++mhOP/10AK688kruu+8+zj333I06h7LU7gRR8g+iEpublSErK4TSujUccsi69e6waBFMmbLu8eyz6yeNXXeFAw6Agw4Kz02bxnMOItXSwoUVW58Cd086Qqhk/WGHHUazZs1+s/2tt97i/PPPB2DXXXdl9913T/r+eXl5a5NHp06d1iaZ6dOnc+WVV/LDDz+wZMkSDj/88I0+h7LU7gQBIRnEnBBSZbYucRx1VFjnHlobU6ZAYSG8/z7ccw/cfnvYPz8/JIuDDoL99w+jq0Rqrdatw2WlZOs30i677MKTTz653rqffvqJRYsWkZ2dzeZlXI9295Tev379+mtfZ2dn82vUgTlw4ECefvppOnbsyKhRo5g0adLGnUA51PWZ4czCCKg+feC66+Dll+H77+GNN+Caa0JCGDECjjwSmjWDLl1gyBB44YXfdpSL1HjXXx8uIyfaxMvKhxxyCEuXLuWhhx4CYPXq1Vx88cUMHDiw3D6B/fbbj3HjxgEwc+ZMPv744wp97s8//0zz5s1ZuXIlY9J0/1ZaE4SZdTOzWWY228yGJtne2sxeN7MPzewjM+sRrW9jZr+a2dTocXc646xp6tcPrYU//xlefz2MqHrttXAlrX59uPVW6N4dttoqtERGjYLi4rijFqkC/fvDyJFhIIpZeB45cpOuIpgZ48eP5/HHH6dt27bsvPPO5OTkcMMNN5R73Nlnn01xcTG77747N910E7vvvjtbVKCJ/5e//IW99tqLww47jHbt2m10/OWxVJs5FX5js2zgM+AwoAiYDPRz95kJ+4wEPnT3u8ysA/C8u7cxszbAs+6+a6qfV1BQ4JowKDVLl8Kbb8LEieFRVBT6PfbZB3r1gt69w6gpkUzwySef0L59+7jDqLDVq1ezcuVKcnJymDNnDocccgifffbZb+6dqEzJflZmNsXdC5Ltn84WRBdgtrvPdfcVwFigd6l9HGgcvd4C+DKN8Uhks83g8MPDpaeFC0P/xZVXhqG3Q4bA734H7drBZZfB22+HUVUiUrmWLl3KfvvtR8eOHTnqqKO466670pocNkY6O6lbAIsSlouAvUrtcw3wkpmdC2wOHJqwLc/MPgR+Aq509zfTGGutZQZ77hke114b+u+eeQYmTIB//ANuvhlatoQTToATT4Tddos7YpGaoVGjRtV+muR0tiCSVYYqfT2rHzDK3VsCPYDRZpYFfAW0dvc9gIuAR8yscaljMbNBZlZoZoXFuoheKbbfHgYPDp3d334b7ivq2BH+/nfYfffw+uabw2UpEanZ0pkgioBWCcst+e0lpD8B4wDc/V0gB9jK3Ze7++Jo/RRgDknubXb3ke5e4O4Fubm5aTiF2m2LLULL4dln4auv4I47wuWpyy4LowIPPhjuuy90gotIzZPOBDEZaGtmeWZWD+gLTCy1z0LgEAAza09IEMVmlht1cmNmOwBtgblpjFU2IDcXzjkH3n033Nl9zTWhFXHaabDttnDMMfDcc7B6ddyRikhlSVuCcPdVwGDgReATYJy7zzCzYWbWK9rtYuB0M5sGPAoM9DCs6g/AR9H6J4Az3f27dMUqFbPTTmEI7axZ4ca8M84I91307Ak77gg33AD/+1/cUYrIpkrbMNeqpmGu8VqxInRs3313uOeiTh04+mg488xQa0rzuUi6VIdhrtnZ2eyWMILj6aef5ttvv+Whhx7i9ttvZ9KkSdSrV4999tln7fadd96ZDh06VOhzGjZsyJIlSzY6zooOc1WpDakU9erBsceGx6xZ8K9/hRvwxo0Lw2bPPBMGDFBtKKmZGjRowNSpU9db16ZNGwoKwu/dSZMm0bBhw/USRM+ePSucIKqaSm1Ipfvd78IQ2S++CEmiaVO48ELYbjs45ZRQM0qkpps0aRI9e/Zk/vz53H333QwfPpz8/Hz+85//MHHiRC699FLy8/OZM2cOc+bMoVu3bnTq1In999+fTz/9FIB58+ax995707lzZ6666qoqPwe1ICRtGjQIrYYBA2Dq1HD56eGHQ9I48EC49FLo1k2TIUnlueCC8G+tMuXnh/I05fn111/XVlzNy8tj/Pjxa7e1adOGM888c71Z3nr16kXPnj055phjgFDP6e6776Zt27a8//77nH322bz22mucf/75nHXWWZx88smMGDGick8sBUoQUiXy80OCuPnmUG321lvhiCOgQwe45JIwnDahaKVIRkl2iSlVS5Ys4Z133uHYY49du255NGPY22+/vbZS7EknncRll1226cFWgBKEVKnGjeHii+G882DsWPjb3+DUU0MhwfPPDyOimjxXuVNCSu2xob/0q6M1a9bQpEmTMhNMsrkmqooa9xKLunXhpJPC5YAXXwyz4w0dCq2br+TigYtZtGB1mOyiZErINJUzFqkKjRo14ueff0663LhxY/Ly8nj88ceBME/EtGnTANh3330ZO3YsQNpKepdHCUJiZQZdu4bSHv/9LxyZ9Ry3rTqbHZjLAEbxOTutmxJSJEMdeeSRjB8/nvz8fN5880369u3LLbfcwh577MGcOXMYM2YM9913Hx07dmSXXXZhwoQJANx2222MGDGCzp078+OPP1Z53LoPQqqXrCwWekuGcyH/4gyWU58TeZiruI6d/PO4o5NqqDrcB5EpqlO5b5GKa92a1ixiOBcxjzwu4FYe51ja8QmnnAJz5sQdoEjtoQQh1UvClJDb8A1/5xLmNtiVc7t9ztix4R6LP/0J5qoyl0jaKUFI9ZJkSsht7/kLw//dnrlzQ8HAMWNCojjtNJg/P+6ApTqoKZfK02ljfkZKEFL99O8ffvOvWROeoyGuzZvDbbeFy0xnngmjR0PbtmFo7Jeai7DWysnJYfHixUoS5XB3Fi9eTE5OToWOUye1ZKyiIrjxxtDgqFs33F9x6aXQqFHckUlVWrlyJUVFRSxbtizuUKq1nJwcWrZsSd26dddbX14ntRKEZLw5c8Io2Mceg222CXNVnHZaqCgrIuXTKCap0XbcMdyV/d57sPPOcNZZsOuuofx4Dfn7RyQWShBSY+y1F/znPyExmEGfPnDAAfDBB3FHJpKZ0pogzKybmc0ys9lmNjTJ9tZm9rqZfWhmH5lZj4Rtl0fHzTKzw9MZp9QcZtCrF3z8Mdx1V5ibYq+9oG9fDY0Vqai0JYhoTukRQHegA9DPzErPjnElYSrSPQhzVt8ZHdshWt4F6AbcWTJHtUgq6tQJI51mzw7Toz7zDLRrB5dfDpswIZdIrZLOFkQXYLa7z3X3FcBYoHepfRxoHL3eAigZrNgbGOvuy919HjA7ej+RCmnUCK69NiSKE04Io57atQsd2uqfEClfOhNEC2BRwnJRtC7RNcCJZlYEPA+cW4FjMbNBZlZoZoXFxcWVFbfUQM2bh4mK3nknjHTq2xcOPhimT487MpHqK50JIlkR89J/s/UDRrl7S6AHMNrMslI8Fncf6e4F7l6Qm5u7yQFLzbf33qHT+q67YNq0MJHRRRdBDIUyRaq9dCaIIqBVwnJL1l1CKvEnYByAu78L5ABbpXisyEbJzg79E599Fu6XuPXWULrjoYfCzdsiEqQzQUwG2ppZnpnVI3Q6Tyy1z0LgEAAza09IEMXRfn3NrL6Z5QFtAQ1WlEq11VZhGtTJk6FNmzB39v77w4cfxh2ZSPWQtgTh7quAwcCLwCeE0UozzGyYmfWKdrsYON3MpgGPAgM9mEFoWcwEXgDOcffV6YpVardOnULfxP33w+efh+ULLtBoJxGV2hBJ8MMPoWzHXXdBq1bhuUePDR8nkqlUakMkRU2awIgR8NZb0LAhHHEE9OsHX38dd2QiVU8JQiSJffYJc2Rfey089RS0bw8PPKB7J6R2UYIQKUP9+uEu7KlTYZdd4NRT4bDDwk13IrWBEoTIBrRvH4oAlox42m03uOkmWLky7shE0ksJQiQFWVlh5rpPPgmd1kOHQufO4TKUSE2lBCFSAdttB08+CePHQ3FxqBR77bVqTUjNpAQhshH69Al1nPr2DTPY/f73MGNG3FGJVC4lCJGN1LQpjB4dRjktWgR77gk33wyrdUun1BBKECKb6KijQuuhZ0+47LJQruPzz+OOSmTTKUGIVILcXHjiCRgzJnRkd+wI//yniv9JZlOCEKkkZmFSohkz4MAD4bzz4NBDYcGCuCMT2ThKECKVbLvt4Lnn4J571t03obuwJRMpQYikgVmYa+Ljj0Pn9amnhhFPP/wQd2QiqVOCEEmjNm3g1Vfhr38No506dgyFAEUygRKESJplZ4c7r99+G+rWhQMOgKuvhlWr4o5MpHxKECJVpEuXMFvdiSfCsGEhUcyfH3dUImVLa4Iws25mNsvMZpvZ0CTbh5vZ1OjxmZn9kLBtdcK20lOVimSkRo3gwQfDcNjp08Mlp7Fj445KJLm0JQgzywZGAN2BDkA/M+uQuI+7X+ju+e6eD/wTeCph868l29y9FyI1yAknrCsj3q8fnHIK/Pxz3FGJrC+dLYguwGx3n+vuK4CxQO9y9u9HmJdapFbIy4M33ghzTjz0EOyxRxgWK1JdpDNBtAAWJSwXRet+w8y2B/KA1xJW55hZoZm9Z2Z9yjhuULRPYXFxcWXFLVJl6tQJ1WAnTYIVK2DffcMd2LpnQqqDdCYIS7KurH/2fYEn3D2xzFnraCLtE4BbzWzH37yZ+0h3L3D3gtzc3E2PWCQm++8fLjkdfni4A/v44+Gnn+KOSmq7dCaIIqBVwnJL4Msy9u1LqctL7v5l9DwXmATsUfkhilQfzZrBhAlhtrqnnoJOnWDatLijktosnQliMtDWzPLMrB4hCfxmNJKZ/Q5oCrybsK6pmdWPXm8F7AvMTGOsItVCVhYMGQKvvw5Ll4Z5Ju69V5ecJB5pSxDuvgoYDLwIfAKMc/cZZjbMzBJHJfUDxrqv91+gPVBoZtOA14Eb3V0JQmqN/fcP90zstx+cfjoMGAC//BJ3VFLbmNeQP00KCgq8sLAw7jBEKtXq1XDddaEju337UFK8ffu4o5KaxMymRP29v6E7qUWqsezsUJbjpZfCHNgFBeEmO5GqoAQhkgEOPTRcctpzz1Cq44wzYPnyuKOSmk4JQiRDtGgROq+HDIGRI0Mtpy++iDsqqcmUIEQySJ06YRjs44+HWk6dOql8uKSPEoRIBjrmGHj//VD876CDYMQIDYWVyqcEIZKhdtkl1G7q2hUGDw6z1i1bFndUUpMoQYhksCZN4JlnQsG/UaPC/ROLFm3wMJGUKEGIZLisrHCfxNNPw6xZoV9i0qS4o5KaQAlCpIbo3Rs++CDUdDr0ULjtNvVLyKZRghCpQdq1C0miZ0+44AI4+WT49de4o5JMtcEEYWZHm9nnZvajmf1kZj+bmQoRi1RTjRuHarDDhsHDD8OBB8JXX8UdlWSiVFoQNwO93H0Ld2/s7o3cvXG6AxORjZeVBVddFRLF9OnQpUu4E1ukIlJJEF+7+ydpj0REKt1RR8Hbb4NZqAz75JNxRySZJJUEUWhmj5lZv+hy09FmdnTaIxORSpGfH/oldtst3GB33XXqvJbU1Elhn8bAUqBrwjoHnkpLRCJS6bbdNgx9Pe20cOlp5ky47z5o0CDuyKQ622CCcPdTqiIQEUmvnBwYPTrcgf1//wezZ4cpTps3jzsyqa5SGcXU0szGm9k3Zva1mT1pZi1TeXMz62Zms8xstpkNTbJ9uJlNjR6fmdkPCdsGRKOnPjezARU7LRFJxgwuvxzGjw+tiM6d4b//jTsqqa5S6YN4gDCX9HZAC+CZaF25zCwbGAF0BzoA/cysQ+I+7n6hu+e7ez7wT6LLVmbWDLga2AvoAlxtZk1TPSkRKV+fPqHzOisrdF4/8UTcEUl1lEqCyHX3B9x9VfQYBeSmcFwXYLa7z3X3FcBYoHc5+/cDHo1eHw687O7fufv3wMtAtxQ+U0RS1LFjKPaXnw/HHgs33qjOa1lfKgniWzM70cyyo8eJwOIUjmsBJJYNK4rW/YaZbQ/kAa9V5FgzG2RmhWZWWFxcnEJIIpJom23gtdegb99w6emMM2DlyrijkuoilQRxKnAc8D/gK+CYaN2GWJJ1Zf190hd4wt1XV+RYdx/p7gXuXpCbm0qjRkRKy8kJ81z/3//BPffAkUfCT6qVIKSQINx9obv3cvdcd9/a3fu4+4IU3rsIaJWw3BL4sox9+7Lu8lJFjxWRTZSVBddfHxLEK6+EsuFFRXFHJXErc5irmQ1x95vN7J8k/+v9vA2892SgrZnlAV8QksAJST7nd0BT4N2E1S8CNyR0THcFLt/A54nIJjrtNGjdOtxQt9de8NxzoY9CaqfyWhAl5TUKgSlJHuVy91XAYMIv+0+Ace4+w8yGmVmvhF37AWPd13WPuft3wF8ISWYyMCxaJyJp1rVrmOc6Kyu0JF54Ie6IJC7mGxi2YGbHuvvjG1oXt4KCAi8sLIw7DJEa44svQtnwjz+GO++EQZuPgSuugIULQzPj+uuhf/+4w5RNZGZT3L0g2bZUOqmTXdrR5R6RGq5FC3jjjdCiOOMMuHzgV6xZsDCMhV2wAAYNCr3bUmOV1wfRHegBtDCz2xM2NQZWpTswEYlfo0YwcSKc2/RhblxyCfNoySgGksNyWLo0tCjUiqixyqvF9CWh/6EX6/c5/AxcmM6gRKT6qFMH7lxyMjswjSHcwv/Ylgn0Zgt+CpebpMYqM0G4+zRgmpk9QrgvoR1hNNOs6M5oEaklbPvWXLrgb7SkiAE8yB94g3/Tne1a14s7NEmjVPogDgPmALcDdwCzo8tPIlJbXH89bLYZ/RjLcxzBXHZgH3uXWefcvuFjJWOlkiD+ARzk7ge6+wHAQcDw9IYlItVK//4wciRsvz2H2atM2rYfvzbamn1v6sX778cdnKRLKgniG3efnbA8F/gmTfGISHXVvz/Mnw9r1tDpq2d5578NaNIEDj4Ynn8+7uAkHVJJEDPM7HkzGxjNy/AMMFlTj4rUbjvuGEqGt2sHvXrBqFFxRySVLZUEkQN8DRwAHAgUA82AI4GeaYtMRKq9bbYJU5kedBCccopKhtc05U45Gk3685G7q89BRJJq1CjUbBo4MJQM//JLuPXWUKpDMlu5X2FUfrtXefuIiNSrBw8/DBdeCP/8J/TrB8uXxx2VbKpyWxCRd8zsDuAx4JeSle6umWxFZK2sLPjHP2C77eDSS+G778Lc1w0bxh2ZbKxUEsQ+0fOwhHUOHFz54YhIprvkEsjNhVNPhcMOC5efmjWLOyrZGBtMEO5+UFUEIiI1x4AB0LhxmMr0gAPgpZegefO4o5KK2mA3kpltYWb/KJn72cz+bmZbVEVwIpK5jjoq3B8xb16YV2LevLgjkopKZZzB/YQCfcdFj5+AB1J5czPrZmazzGy2mQ0tY5/jzGymmc2I6j6VrF9tZlOjx8RUPk9EqpdDDoFXXw39EfvuCzNmxB2RVEQqEwZNdff8Da1Lclw28BmhllMRYWa4fu4+M2GftsA44GB3/97Mtnb3b6JtS9w95e4tTRgkUn1Nnx7mlVi+PMxQ17lz3BFJiU2dMOhXM9sv4c32BX5N4bguwGx3nxtVfx0L9C61z+nACHf/HqAkOYhIzbLrrvDmm7DFFqE0x+uvxx2RpCKVBHEWMMLM5pvZAkJF1zNTOK4FsChhuShal2hnYGcze9vM3jOzbgnbcqI+j/fMrE+yDzCzQSV9I8XFxSmEJCJx2XHHMNf19ttD9+4wYULcEcmGbDBBuPtUd+8I7A7s5u57RHNFbIgle7tSy3WAtoQSHv2Ae82sSbStddTsOQG41cx2TBLbSHcvcPeC3NzcFEISkThttx385z/QsSP88Y8wenTcEUl5UhnFdL6ZNSZ0VP/DzP5rZl1TeO8ioFXCckvCLHWl95ng7ivdfR4wi5AwcPcvo+e5wCRgjxQ+U0SquS23hFdeCcNfTz4Z7rgj7oikLKlcYjrV3X8CugJbA6cAN6Zw3GSgrZnlmVk9oC9QejTS04T5JTCzrQiXnOaaWVMzq5+wfl9gJiJSI5TUb+rdG849F265Je6IJJlUEkTJpaIewAPR5aVkl4/W4+6rgMHAi8AnwDh3n2Fmw8yspL7Ti8BiM5sJvA5c6u6LgfZAoZlNi9bfmDj6SUQyX04OPP44HH88DBkCw4apEmx1k8ow1wcInct5QEcgG5jk7p3SH17qNMxVJDOtXh3Kcjz0EAwdCjfcALbBP0GlspQ3zDWVWkx/AvKBue6+1My2JFxmEhHZZNnZ8MAD0KBBmE/i119h+HAlieoglVpMa8ysDXCimTnwlruPT3dgIlJ7ZGXBXXeFy0633QbLlsGdd2pOibhtMEGY2Z3ATsCj0aozzOxQdz8nrZGJSK1iFloOiS2J++6DOqlc55C0SOVHfwCwq0edFWb2IPBxWqMSkVrJLPRBbLYZ/PnPoSXx8MNQt27ckdVOqSSIWUBrYEG03Ar4KG0RiUitZgZXXRUuNw0ZEuo3PfYY1K8fd2S1T5kJwsyeIdz5vAXwiZl9EG3qDLxbBbGJSC126aXhctO550KfPvDUU2FZqk55LYi/JVlnwH6EshgiImk1eHBoSQwaBD17wsSJsPnmcUdVe5SZINz9PyWvzSyfUBPpOGAecHf6QxMRgdNOC5eXBg4MSeLZZ5Ukqkp5l5h2JpTH6AcsBh4j3FinKUhFpEqddFIY8nryyXDEEaFMh5JE+pU3yvhT4BDgSHffz93/CayumrBERNbXv3+o/vrmm9CjByxZEndENV95CeKPwP+A183sHjM7hBRqMImIpMsJJ4Rhr2+9pSRRFcpMEO4+3t2PB9oRym1fCGxjZnelWO5bRKTS9esHjzwC77wTJh76+ee4I6q5Upkw6Bd3H+PuPQlzOkwFhqY9MhGRMhx/fEgS776rJJFOFap04u7fufu/3P3gdAUkIpKK446DRx+F996Dbt3gp5/ijqjmUSksEclYxx4b7rL+4AMliXRQghCRjPbHP4YkMXkyHH44/Phj3BHVHGlNEGbWzcxmmdlsM0vab2Fmx5nZTDObYWaPJKwfYGafR48B6YxTRDLb0UfDuHFQWBiShFoSlSNtCcLMsoERQHegA9DPzDqU2qctcDmwr7vvAlwQrW8GXA3sBXQBrjazpumKVUQy31FHhSlMp0zRENjKks4WRBdgtrvPdfcVwFigd6l9TgdGuPv3AO7+TbT+cODlqFP8e+BloFsaYxWRGqBPn3Ud1z17wi+/xB1RZktngmgBLEpYLorWJdoZ2NnM3jaz98ysWwWOxcwGmVmhmRUWFxdXYugikqmOOWbdHde9e4eJh2TjpDNBJLvr2kst1wHaAgcSaj7da2ZNUjwWdx/p7gXuXpCbm7uJ4YpITdGvH4waBa+9Fi49LVsWd0SZKZ0JoogwuVCJlsCXSfaZ4O4r3X0eYXKitikeKyJSppNOgnvvhRdfDK2KFSvijijzpDNBTAbamlmemdUjVIadWGqfp4GDAMxsK8Ilp7nAi0BXM2sadU53jdaJiKTs1FPh7rtD9dfjj4eVK+OOKLOkbTpwd19lZoMJv9izgfvdfYm94Q8AABC4SURBVIaZDQMK3X0i6xLBTEKl2EvdfTGAmf2FkGQAhrn7d+mKVURqrjPOCInh3HNDRdhHHoE6afvNV7OY+28u7WekgoICLywsjDsMEammhg+Hiy4K/ROjR0N2dtwRVQ9mNsXdC5JtUx4VkVrhwgtDS+Kyy6BuXbj/fiWJDVGCEJFaY8iQ0Fl91VUhSYwcGWaqk+SUIESkVrnyytCSGDYMNtsMbrsNTFOhJaUEISK1zjXXwNKl8Le/QcOGcMMNcUdUPSlBiEitYwY33xzqNf31r7D55nDFFXFHVf0oQYhIrWQGI0aEek1XXhmSxAUXxB1V9aIEISK1VlZWGM20dGkY5bT55nD66XFHVX2o/15EarU6dcLNc927h5vqHnlkw8fUFkoQIlLr1asHTz4JBxwAJ58MTz8dd0TVgxKEiAjQoAFMnAidO4e6TS+q+psShIhIiUaN4PnnoUOHUCb8jTfijiheShAiIgmaNoWXXoLttw+z0n3wQdwRxUcJQkSklNxceOUV2Gor6NYNPv447ojioQQhIpJEixbw6quhHEfXrjBnTtwRVT0lCBGRMuTlhctNK1bAYYfBl7VsXkslCBGRcnToAC+8AMXFoSXxXS2auiytCcLMupnZLDObbWZDk2wfaGbFZjY1epyWsG11wvrSU5WKiFSZzp3DENjZs6FHj1DDqTZIW4Iws2xgBNAd6AD0M7MOSXZ9zN3zo8e9Cet/TVjfK11xioik4qCD4LHHoLAQ+vSB5cvjjij90tmC6ALMdve57r4CGAv0TuPniYikVe/ecN99ofP6hBNg1aq4I0qvdCaIFsCihOWiaF1pfzSzj8zsCTNrlbA+x8wKzew9M+uT7APMbFC0T2FxcXElhi4iktyAAXDrrfDUUzBoEKxZE3dE6ZPOBJFsjiYvtfwM0MbddwdeAR5M2NY6mkj7BOBWM9vxN2/mPtLdC9y9IDc3t7LiFhEp1/nnw9VXwwMPwCWXgJf+zVZDpLPcdxGQ2CJoCaw3SMzdFycs3gPclLDty+h5rplNAvYAauFIZBGpjq6+OoxoGj4cttyyZk44lM4WxGSgrZnlmVk9oC+w3mgkM2uesNgL+CRa39TM6kevtwL2BWamMVYRkQoxC5eaTjopTDg0YkTcEVW+tLUg3H2VmQ0GXgSygfvdfYaZDQMK3X0icJ6Z9QJWAd8BA6PD2wP/MrM1hCR2o7srQYhItZKVFTqtf/wRBg+GZs2gX7+4o6o85jXk4llBQYEXFhbGHYaI1ELLloWaTe+8A88+G26oyxRmNiXq7/0N3UktIrKJcnJgwoRw1/XRR9ecCrBKECIilWCLLeDf/4attw53W8+aFXdEm04JQkSkkjRvHor7ZWeHy0xffBF3RJtGCUJEpBLttFNoSXz/feiX+P77uCPaeEoQIiKVbM894emn4bPP4MgjYenSuCPaOEoQIiJpcPDBMGZMGNl0/PGZWbdJCUJEJE2OOSbcQPfss6FuU6bdVZDOUhsiIrXeWWfB11/DtdeGEU433hh3RKlTghARSbOrrw5J4qabYJtt4MIL444oNUoQIiJpZgZ33AHffgsXXRRaEv37xx3VhilBiIhUgexsePjhkCROOSW0JA49NO6oyqdOahGRKlK/PowfD+3ahZIcU6fGHVH5lCBERKpQkybhRromTaB7d5g/P+6IyqYEISJSxVq0gBdeWFcFdvHiDR8TByUIEZEYdOgAEyeGFkR1vds6rQnCzLqZ2Swzm21mQ5NsH2hmxWY2NXqclrBtgJl9Hj0GpDNOEZE47L8/PPIIvPdemGiout1tnbYEYWbZwAigO9AB6GdmHZLs+pi750ePe6NjmwFXA3sBXYCrzaxpumIVEYnL0UfD7beH1sTgwdXrbut0tiC6ALPdfa67rwDGAr1TPPZw4GV3/87dvwdeBrqlKU4RkVgNHgxDh8K//gXXXx93NOukM0G0ABYlLBdF60r7o5l9ZGZPmFmrCh4rIlIj3HADnHQSXHUVPPBA3NEE6UwQlmRd6cbTM0Abd98deAV4sALHYmaDzKzQzAqLi4s3KVgRkTiZwX33hYmGTj89DIWNWzoTRBHQKmG5JfBl4g7uvtjdl0eL9wCdUj02On6kuxe4e0Fubm6lBS4iEoe6deGJJ2D33UMl2MmT440nnQliMtDWzPLMrB7QF5iYuIOZNU9Y7AV8Er1+EehqZk2jzumu0ToRkRqtUSN4/vlQiqNnT5g7N75Y0pYg3H0VMJjwi/0TYJy7zzCzYWbWK9rtPDObYWbTgPOAgdGx3wF/ISSZycCwaJ2ISI237bbhEtOqVdCjR3w30plXpzFVm6CgoMALCwvjDkNEpNK8+WYo6NelC7z8MuTkVP5nmNkUdy9Itk13UouIVFP77w8PPQRvvQUDBsCaNVX7+Sr3LSJSjR1/PCxcCEOGwPbbw803V91nK0GIiFRzl1wSajbdcgu0aQNnn101n6sEISJSzZnBbbfBokVw7rnQqlUo8Jdu6oMQEckAderAo4/CnntC375Vc4+EEoSISIbYfHN49tkwp3XPnjBvXno/TwlCRCSDbLNNuEdi5Urovu+PfNeqI2Rlhc6JMWMq9bOUIEREMky7djDhnJeY91UOvYvuYJnXgwULYNCgSk0SShAiIhlo/9GDeIiTeYv9GcCDrMHCtHRXXFFpn6FRTCIimWjhQo5nAQtpzRIaYiUFrxcurLSPUIIQEclErVvDggVcyt9+u76S6BKTiEgmuv562Gyz9ddttlmlTkmnBCEikon694eRI0P9DbPwPHJkWF9JdIlJRCRT9e9fqQmhNLUgREQkKSUIERFJSglCRESSUoIQEZGklCBERCSpGjMntZkVAws24S22Ar6tpHDiVFPOA3Qu1VVNOZeach6waeeyvbvnJttQYxLEpjKzwrIm7s4kNeU8QOdSXdWUc6kp5wHpOxddYhIRkaSUIEREJCkliHVGxh1AJakp5wE6l+qqppxLTTkPSNO5qA9CRESSUgtCRESSUoIQEZGkakWCMLMcM/vAzKaZ2QwzuzZaP8rM5pnZ1OiRH603M7vdzGab2Udmtme8Z7A+M8s2sw/N7NloOc/M3jezz83sMTOrF62vHy3Pjra3iTPuZJKcS6Z+J/PN7OMo5sJoXTMzezn6Xl42s6bR+kw8l2vM7IuE76VHwv6XR+cyy8wOjy/y3zKzJmb2hJl9amafmNnemfi9lHEe6f9O3L3GPwADGkav6wLvA78HRgHHJNm/B/Dv6LjfA+/HfQ6l4rsIeAR4NloeB/SNXt8NnBW9Phu4O3rdF3gs7thTOJdM/U7mA1uVWnczMDR6PRS4KYPP5RrgkiT7dgCmAfWBPGAOkB33OSTE9yBwWvS6HtAkE7+XMs4j7d9JrWhBeLAkWqwbPcrrne8NPBQd9x7QxMyapzvOVJhZS+AI4N5o2YCDgSeiXR4E+kSve0fLRNsPifavFkqfywZU2++kHIk//9LfS6adS1l6A2Pdfbm7zwNmA11ijgkAM2sM/AG4D8DdV7j7D2TY91LOeZSl0r6TWpEgYO2ljKnAN8DL7v5+tOn6qDk53MzqR+taAIsSDi+K1lUHtwJDgDXR8pbAD+6+KlpOjHXteUTbf4z2ry5Kn0uJTPtOIPzB8ZKZTTGzQdG6bdz9K4DoeetofSaeC8Dg6Hu5v+SyDNX7XHYAioEHosuY95rZ5mTe91LWeUCav5NakyDcfbW75wMtgS5mtitwOdAO6Aw0Ay6Ldk/2V3bs44HNrCfwjbtPSVydZFdPYVusyjgXyLDvJMG+7r4n0B04x8z+UM6+mXgudwE7AvnAV8Dfo32r87nUAfYE7nL3PYBfCJeUylJdz6Ws80j7d1JrEkSJqGk2Cejm7l9FzcnlwAOsa4YVAa0SDmsJfFmlgSa3L9DLzOYDYwmXlm4lNIVLpo9NjHXteUTbtwC+q8qAy/GbczGzhzPwOwHA3b+Mnr8BxhPi/rrkEkX0/E20e8adi7t/Hf2RtQa4h8z4XoqAooSrBU8QftFm2veS9Dyq4jupFQnCzHLNrEn0ugFwKPBpwj8SI1yHnB4dMhE4ORrV8Hvgx5ImaZzc/XJ3b+nubQidzq+5e3/gdeCYaLcBwITo9cRomWj7ax71YsWtjHM5MdO+EwAz29zMGpW8BroS4k78+Zf+XjLqXEpdiz+K9b+XvhZGzOUBbYEPqjLmsrj7/4BFZva7aNUhwEwy7Hsp6zyq4jups+FdaoTmwINmlk1IiuPc/Vkze83McglNsqnAmdH+zxNGNMwGlgKnxBBzRVwGjDWz64APiTqzoufRZjab0HLoG1N8FTEmA7+TbYDxUf9/HeARd3/BzCYD48zsT8BC4Nho/0w8l9EWhhw7YZTTGQDuPsPMxhF+8a4CznH31bFEnty5hH9T9YC5hJ91Fpn3vSQ7j9vT/Z2o1IaIiCRVKy4xiYhIxSlBiIhIUkoQIiKSlBKEiIgkpQQhIiJJKUFIRorKcFyQsPyimd2bsPx3M7uokj9zyYb3qvB75peqwnmNmV2SwnEWDdNuXAkx1DOzNxJuthQBlCAkc70D7ANgZlnAVsAuCdv3Ad6OIa6KyieMva+oHsA0d/9pUwNw9xXAq8Dxm/peUrMoQUimepsoQRASw3TgZzNrGhX4aw98aGYNzexVM/uvhTkOegOY2U1mdnbJm0V/uV8cvb7UzCZHRdCuTfbhyfYxszYWavXfY2HekZeiO/cxs87Rvu+a2S1mNj266WkYcLyFev4lv6A7mNkkM5trZueVcf79ie4A3sDnTopaW29E+3Q2s6cszIVwXcL7PR29p8haShCSkaJ6QavMrDUhUbxLmOdjb6AA+Cj6y3gZcFRUfO4g4O9RGY+xrP8X83HA42bWlVCaoAvhr/tOVqrw3gb2aQuMcPddgB+AP0brHwDOdPe9gdXROawA/kyYpyPf3R+L9m0HHB69/9VmVjfJj2BfILHQYVmfC7DC3f9AmCtkAnAOsCsw0MxKqvtOJxRIFFlLCUIyWUkroiRBvJuw/E60jwE3mNlHwCuEssfbuPuHwNZmtp2ZdQS+d/eFhNpDXQklS/5L+GXdttTnlrfPPHefGr2eArSxUAeskbuXxPTIBs7ruaiW/7eEQnLbJNmnmbv/nLD8m89N2DYxev4YmBEVRFxOKNnQCkK1Y2BFSR0mEag9tZikZirph9iN8BfwIuBi4Cfg/mif/kAu0MndV1qoHpsTbXuCUMRwW0KLAkJC+au7/6ucz026j4UpXZcnrFoNNCB5+eXylH6PZP9PV5lZVlTJM9kxDZK835pS+60p9d71CS0uEUAtCMlsbwM9ge+issffEaZi3JvQmoBQ4vybKDkcBGyfcPxYQgHDY1g3I9+LwKlm1hDAzFqY2dasL5V91nL37wn9I7+PViUWTfwZ2Ji/2mcRJpKpFNGlpmJ3X1lZ7ymZTwlCMtnHhNFL75Va92N0eQZgDFBgZoWE1sSnJTu6+wzCL+cvEmYYe4lwCehdM/uYkDjW+wWeyj5J/AkYaWbvEloUP0brXyd0Sid2UqfiOeDACuy/IQcRqpmKrKVqriJVwMwaejQvupkNBZq7+/mb8H7NCfMnH1ZJ8T0FXO7usyrj/aRmUB+ESNU4wswuJ/yfWwAM3JQ3c/evomGtjTf1XohouO3TSg5SmloQIiKSlPogREQkKSUIERFJSglCRESSUoIQEZGklCBERCSp/wcUk3OAq0gCHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nWavelengths = np.array([350, 450, 550, 650])\n",
    "mu_as = np.array([30.16, 13.5, 6.85, 3.68])\n",
    "mu_ss = np.array([210.4, 121.6, 74.7, 55.48])\n",
    "gs = np.array([0.702, 0.728, 0.759, 0.787])\n",
    "n_ts = np.array([1.488, 1.448, 1.447, 1.433])\n",
    "n_m = 1.5\n",
    "d = 0.1\n",
    "\n",
    "mu_atts = mu_as + mu_ss\n",
    "absFacts = mu_as / mu_atts\n",
    "\n",
    "alpha_criticals = np.zeros(4)\n",
    "for n, n_t in enumerate(n_ts):\n",
    "    if n_t > n_m:\n",
    "            alpha_criticals[n] = asin(n_m / n_t)\n",
    "    else:\n",
    "        alpha_criticals[n] = float('nan')\n",
    "\n",
    "N_bins = int(1e3)\n",
    "dz = d / N_bins\n",
    "\n",
    "pDW = np.zeros((4, N_bins))\n",
    "dDR = np.zeros(4)\n",
    "dDT = np.zeros(4)\n",
    "\n",
    "N_photons = int(5e6)\n",
    "\n",
    "for n, (mu_att, absFact, g, n_t, alpha_critical) in enumerate(zip(mu_atts, absFacts, gs, n_ts, alpha_criticals)):\n",
    "    for nph in range(N_photons):\n",
    "        dDR[n], dDT[n], pDW[n] = RandomWalk(mu_att, absFact, g, n_t, n_m, alpha_critical, d, dz, dDR[n], dDT[n], pDW[n])\n",
    "\n",
    "        \n",
    "dDR /= N_photons\n",
    "dDT /= N_photons\n",
    "dDW = np.sum(pDW, axis = 1) / N_photons\n",
    "\n",
    "print('\\n')\n",
    "print(tabulate({'Wavelengths': nWavelengths, 'Diffuse Transmittance': dDT, 'Diffuse Reflectance': dDR, 'Absorption': dDW}, headers = 'keys'))\n",
    "print('\\n')\n",
    "\n",
    "pz = np.arange(0, d, dz)\n",
    "pDW = pDW / (dz * N_photons)\n",
    "plt.figure()\n",
    "plt.semilogy(pz, pDW[0], 'r', label = str(nWavelengths[0]) + ' nm')\n",
    "plt.semilogy(pz, pDW[1], 'g', label = str(nWavelengths[1]) + ' nm')\n",
    "plt.semilogy(pz, pDW[2], 'b', label = str(nWavelengths[2]) + ' nm')\n",
    "plt.semilogy(pz, pDW[3], 'k', label = str(nWavelengths[3]) + ' nm')\n",
    "plt.xlabel('Depth (cm)')\n",
    "plt.ylabel('Absorption (cm$^{-1}$)')\n",
    "plt.legend()\n",
    "\n",
    "p = np.polyfit(nWavelengths, dDW, 3)\n",
    "nWaves = np.arange(350, 660, 10)\n",
    "fitted = np.polyval(p, nWaves)\n",
    "plt.figure()\n",
    "plt.plot(nWavelengths, dDW, 'ro', label = 'Original')\n",
    "plt.plot(nWaves, fitted, 'b', label = 'Fitted')\n",
    "plt.xlabel('Wavelength (nm)')\n",
    "plt.ylabel('Absorption')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the above table and plots, it is evident that the lower the wavelength, the higher the light absorption. The depth profile demonstrate that the absorption rate drops off more quickly for the lower wavelengths of light than the higher wavelengths."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "1. Zhu, C. and Liu, Q., 2013. Review of Monte Carlo modeling of light transport in tissues. Journal of biomedical optics, 18(5), p.050902.\n",
    "2. Wang, L., Jacques, S.L. and Zheng, L., 1995. MCML—Monte Carlo modeling of light transport in multi-layered tissues. Computer methods and programs in biomedicine, 47(2), pp.131-146.\n",
    "3. Bashkatov, A.N., Genina, E.A. and Tuchin, V.V., 2011. Optical properties of skin, subcutaneous, and muscle tissues: a review. Journal of Innovative Optical Health Sciences, 4(01), pp.9-38."
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}