{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dual Rotating Retarder Mueller Polarimetry\n",
    "\n",
    "_written by Jaren N. Ashcraft_\n",
    "\n",
    "In order to measure how a system transforms the polarization state of light, we can measure the system's Mueller Matrix. There are a host of different methods by which we can do so. In this tutorial, we review the popular Dual Rotating Retarder Polarimeter (DRRP) for full Mueller matrix polarimetry. This consists of:\n",
    "\n",
    "- a source\n",
    "- a linear polarizer at a fixed angle\n",
    "- a quarter-wave plate that can rotate\n",
    "- the system you are trying to measure the Mueller matrix of\n",
    "- a quarter-wave plate that can rotate\n",
    "- a linear polarizer at a fixed angle\n",
    "- a detector to measure power\n",
    "\n",
    "The theory of this polarimeter is covered in Chipman [1], and also Azzam [2], but is reproduced in brief here. Our first polarizer/waveplate pair forms our *polarization state generator* (PSG). The latter pair forms our *polarization state analyzer* (PSA). The first waveplate must rotate $N_{meas} > 16$ times at some angular step size $\\delta\\theta$ to completely determine the Mueller matrix. The second waveplate must also rotate $N_{meas}$ times, but with a much greater step size. In this tutorial we use $5\\delta\\theta$.\n",
    "\n",
    "For each $\\delta\\theta$ we record the power on the detector. We next must understand that detectors can only measure power! This means, that it effectively measures the first row of the PSA matrix and the first column of the PSG matrix. This can be shown by working through the Mueller calculus, which we encourage the readers to try for themselves.\n",
    "\n",
    "$$\\mathbf{S}_{out} = [\\mathbf{PSA}] \\mathbf{M}_{system} [\\mathbf{PSG}] \\mathbf{S_{0}}$$\n",
    "\n",
    "For each measurement iteration we compute our PSA and PSG, and then store those results in our polarimetric data reduction matrix $\\mathbf{W}$. This is an $N_{meas} \\times 16$ matrix that contains what our PSA and PSG are doing. $\\mathbf{W}$ is given by:\n",
    "\n",
    "$$\\mathbf{W} = [\\mathbf{PSA}_{0,j}] \\otimes [\\mathbf{PSG}_{i,0}]$$\n",
    "\n",
    "Where $i,j$ are the row and column indices, respectively, and $\\otimes$ is the Kronecker product. For each measurement, we record the power on the detector and store it in a vector that is $N_{meas}$ long, called $\\mathbf{P}$, and we compute the unraveled Mueller matrix using the following relation:\n",
    "\n",
    "$$\\mathbf{M}_{meas} = \\mathbf{W}^{-1}\\mathbf{P}$$\n",
    "\n",
    "**References**\n",
    "- [1] Chipman, R., Lam, W.S.T., & Young, G. (2018). Polarized Light and Optical Systems (1st ed.). CRC Press. https://doi.org/10.1201/9781351129121\n",
    "- [2] Ellipsometry and Polarized Light: by R.M.A. Azzam and N.M. Bashara; published by North-Holland Publishing Company, Amsterdam, 1977; pp. xvii  529; price 196 Dfl. (n.d.). Surface Technology, 7(3), 269–269. https://doi.org/info:doi/\n",
    "\n",
    "For the case where $N_{meas} > 16$, The Moore-Penrose pseudo inverse is used to invert $\\mathbf{W}$. Now that the theory is out of the way, let's test it out! We begin by constructing a random Mueller matrix composed of a polarizer and retarder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.5       , 0.29399553, 0.40443371, 0.        ],\n",
       "       [0.28578716, 0.16804029, 0.23116392, 0.        ],\n",
       "       [0.40673081, 0.23915408, 0.3289913 , 0.        ],\n",
       "       [0.05381217, 0.03164107, 0.04352691, 0.        ]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from katsu.mueller import linear_polarizer,linear_retarder\n",
    "\n",
    "M_to_measure = linear_retarder(np.random.random(),np.random.random()) @ linear_polarizer(np.random.random())\n",
    "display(M_to_measure)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing a simulated measurement\n",
    "We begin by creating a simulated series of power measurements by setting up our polarization state generators and analyzers. We then propagate a Stokes vector through this system and grab the $I$ parameter in order to observe the power as we step through the experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "from katsu.mueller import stokes_from_parameters\n",
    "\n",
    "# How much more the PSA rotates than the PSG\n",
    "PSA_ANGULAR_INCREMENT = 5\n",
    "\n",
    "# katsu has a simulator built-in, so let's try use that\n",
    "thetas = np.linspace(0,2*np.pi,16) # 16 measurements between 0 and pi\n",
    "\n",
    "# set up polarization state generator\n",
    "psg_polarizer = linear_polarizer(0, shape=thetas.shape)\n",
    "psg_retarder  = linear_retarder(thetas, np.pi / 2, shape=thetas.shape)\n",
    "PSG = psg_retarder @ psg_polarizer\n",
    "\n",
    "# set up polarization state generator\n",
    "psa_polarizer = linear_polarizer(0, shape=thetas.shape)\n",
    "psa_retarder  = linear_retarder(PSA_ANGULAR_INCREMENT * thetas, np.pi / 2, shape=thetas.shape)\n",
    "PSA = psa_polarizer @ psa_retarder\n",
    "\n",
    "# set up system Mueller matrix\n",
    "Msys = PSA @ M_to_measure @ PSG\n",
    "\n",
    "# propagate Stokes vector to get power\n",
    "Sin = stokes_from_parameters(1, 0, 0, 0)\n",
    "Sout = Msys @ Sin\n",
    "power_measured = Sout[..., 0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we plot the power observed as the PSG and PSA step in angular increments of $\\delta\\theta$ and $5\\delta\\theta$, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAGwCAYAAABGogSnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB1AElEQVR4nO3de3yT9d3/8XeakPSUCqVH2tLzgRYBRyfjRsS5Ip5Q8ISoAzrFiaIicx4RlDFRh4o6J4o/8IjiFFGnA4EJ0wk6YKJSKG1pKS09UmiTHtImuX5/sGaG9vrSpmmuK1ffz8eD+x5pkn778jL9mOvKdekkSZJARERENEAEKL0AIiIiIl/i8ENEREQDCocfIiIiGlA4/BAREdGAwuGHiIiIBhQOP0RERDSgcPghIiKiAYXDDxEREQ0oBqUXoDZOpxPt7e3Q6/XQ6XRKL4eIiIh6QJIkOBwOGI1GBASI39vh8HOa9vZ2/Otf/1J6GUREROSBCRMmIDAwUHgfDj+n0ev1AICRI0e6/rc3HTx4EFlZWV5/Xi1gGzH2EWMfMfYRYx95/tLG4XDgxx9/7NHvbg4/p+nc1aXX6/tl+ImJiemX59UCthFjHzH2EWMfMfaR529tenLICg94JiIiogGFw4+P1dTUKL0E1WIbMfYRYx8x9hFjH3labMPhh4iIiAYUDj8+lpaWpvQSVIttxNhHjH3E2EeMfeRpsQ2HHx+rrKxUegmqxTZi7CPGPmLsI8Y+8rTYhsOPj7W2tiq9BNViGzH2EWMfMfYRYx95WmyjyuHnxRdfRFJSEgIDAzFu3Dh8++23svfdsGEDcnNzMXjwYISEhGDMmDF488033e4zZ84c6HQ6tz8XX3xxf/8Y3TrTiZcGMrYRYx8x9hFjH3mSJMEeYES1xYbGNjskSVJ6SaqixW1Hdef5Wb9+PRYuXIhVq1Zh3LhxWLlyJaZMmYLCwkJERUV1uX94eDgefvhhZGVlwWg04m9/+xvy8/MRFRWFKVOmuO538cUXY+3ata6/m0wmn/w8pxs+fLgi39cfsI0Y+4ixjxj7dGW12bGlqAEb99ehytIOoAAAEGs2YlpOJCanhyPUpLpfkz6nxW1Hde/8PPPMM5g7dy7y8/ORnZ2NVatWITg4GGvWrOn2/hdccAGmT5+OESNGIDU1FXfffTdGjRqFr776yu1+JpMJMTExrj9DhgzxxY/TxaFDhxT5vv6AbcTYR4x9xNjH3e6KJtzwzn6s2lWJaku729eqLe1YtasSN7yzH7srmhRaoXpocdtR1fDT3t6OPXv2IC8vz3VbQEAA8vLysHPnzjM+XpIkbNu2DYWFhTj//PPdvrZ9+3ZERUUhMzMT8+bNw/Hjx72+fiIiUr/dFU1YtLkENrsTEoDTd3J13mazO7FocwkHIA1S1ft59fX1cDgciI6Odrs9OjoaBw8elH1cY2Mj4uLiYLPZoNfr8Ze//AWTJ092ff3iiy/GVVddheTkZJSUlOChhx7CJZdcgp07d8qesttisbhdFdZkMnllV1lkZGSfn0Or2EaMfcTYR4x9TrHa7Fi6tRSS1HXoOZ303/+zdGsp1s3MGbC7wLS47Wjin6TZbMZ3330Hq9WKbdu2YeHChUhJScEFF1wAALj++utd9z377LMxatQopKamYvv27fjVr37V7XOOHDkSLS0trr/n5+fjzjvvRGxsLEpKSgCcGsokSUJtbS0AID09HRUVFWhtbUVgYCASEhJQVFQEAIiKikJAQACOHTuGuro6pKamorq6Gs3NzTCZTEhKSkJhYSEAICIiAkajEceOHQMAJCcno66uDlarFYMGDUJaWhoOHDgA4NQxT0FBQa6PIiYlJaGhoQFNTU3Q6/XIzMzEgQMHIEkSBg8eDLPZjKNHjwI4tR+3qakJJ0+ehE6nw4gRI1BYWAiHw4GwsDAMGTIER44cAQDEx8ejpaUFDQ0NAIDs7GwcOnQIdrsdZrMZERERKC0tBQAMGzYMNpvN9e5aVlYWDh8+jPb2doSEhCA6OhqHDx8GAMTGxsJut6Ourg4dHR0YMmQIysvL0dbWhqCgIMTFxaG4uNjVG/jf2UbT0tJQWVnp6j18+HDX27ORkZEwGAyoqqoCAKSkpKCmpgbNzc0wGo1ISUlxDdRDhw6FyWRy611fXw+LxQKDwYCMjAwUFBS4egcHB6OiogIAkJiYiBMnTsj2DgsLQ3l5OQAgISEBFotFtnd4eDjKysoAAHFxcWhtbXX1HjFiBGpqalBXV4fQ0FBERka69W5vb0d9fT0AIDMzE2VlZbDZbAgJCUFMTIxrm42JiYHT6XTbZo8ePerqHR8f77bN6nQ6V+/U1FRUVVWhpaUFJpMJiYmJwt61tbWwWq3d9g4MDOx2mz2995AhQxAaGuq2zTY2NqKxsREBAQHIysrCwYMH4XQ6YTQaERQU5NbbarXixIkTXbbZ7nq3tbV1u82GhoYiKiqq220WADIyMnDkyBHYbDYEBwf3+TWiurra1dubrxGdrz3+/BrR2bsvrxFvflmANrsTPdX5DtDH31diTEizq7caXyOKi4vR0dHh9deI+vp61NXVqf41ojcHquskFR3W3t7ejuDgYLz//vuYNm2a6/bZs2fj5MmT+Oijj3r0PLfccguOHj2KzZs3y94nMjISy5Ytw29/+1u32+12O3bs2IGUlJR+eeenoKAA2dnZfX4eLfLnNpIkocnmQGuHA0GD9Agz6Xt0cb3e8Oc+vsA+Yuxz6t/TOe8V/Pfg5p7TAYgxG/Haddle//faH/jLtuNwOLBv3z5MmjQJBoP4vR1VvfNjNBoxduxYbNu2zTX8OJ1ObNu2DfPnz+/x8zidTthsNtmvV1RU4Pjx44iNjZW9j9ls9qur2JIyun5a5BR+WoRIfZpsjl4PPsCpd3+qLO2w2BwIC+S/z1qgun+KCxcuxOzZs5Gbm4tzzz0XK1euRHNzM/Lz8wEAs2bNQlxcHJYvXw4AWL58OXJzc5GamgqbzYbPPvsMb775Jl566SUAgNVqxWOPPYarr77a9dbefffdh7S0NLePwvtKSkqKz7+nv/C3NrsrmrB0ayls3byF3vlpkbW7q7A4Lxm58WF9/n7+1sfX2EeMfYDWDkefHt/SMTCHHy1uO6r6tBcAzJgxAytWrMDixYsxZswYfPfdd9i0aZNrX255eblrnyEANDc34/bbb0dOTg4mTJiADz74AG+99RZuueUWAIBer8f333+PK664AhkZGbj55psxduxYfPnll4qc60eLV8f1Fn9qo8SnRfypjxLYR4x9gKBBfXs3P7iPj/dXWtx2VDnCzp8/X3Y31/bt293+vmzZMixbtkz2uYKCgoTH/vhac3Oz0ktQLX9po9SnRfylj1LYR4x9gDCTHrFmI6ot7Wf8d/enOo/5MZsG5vCjxW1Hde/8aJ3RaFR6CarlL222FDW43vHpic53gLYUNfTp+/pLH6Wwjxj7ADqdDtNyPPvY9rScyAF5sDOgzW2Hw4+PaXHfqbf4QxtJkrBxf12v/qux08b9dX26ZpA/9FES+4ixzymT08NhMgSgp2OMDoDJEIDJ6eH9uSxV0+K2w+HHx0Qnaxzo/KGNNz4t4il/6KMk9hFjn1NCTQYszkuGToczDkA6ADodsDgveUB/alOL2w6HH6Je8ManRYhIWbnxYVg2JdX1DtDpQ1DnbSZDAJZNSfXKpzVJXQbuKKuQoUOHKr0E1fKHNkp+WsQf+iiJfcTYx11ufBjWzczp9jxdMf89T9dFGUMRYhyYBzn/lBa3HQ4/PqbEx+v9hT+0UfLTIv7QR0nsI8Y+XYWaDJg+MgrTciJRUduAQcGhCB6kh7kfztDuz7S47XC3l491XouHuvKHNkp+WsQf+iiJfcTYR55Op4PleA1izCaEBRo4+JxGi9sOhx+iXuKnRYiI/BuHHx9LTk5Wegmq5S9tlPq0iL/0UQr7iLGPGPvI02IbDj8+Vl9fr/QSVMuf2ijxaRF/6qME9hFjHzH2kafFNjzg2ccsFotPv58kSWiyOdDa4UDQID3CVHwgn6/b9JWvPy3ib318jX3E2EeMfeRpsQ2HHx8zGHyT3Gqzd/tLOfa/v5Qnp4er7qRdvmrjTT/9tIjF5kBLh6PfPi3ij318iX3E2EeMfeRpsY1O6sv59jXIbrdjx44dGD16NPR6/zy/w+6KJizdWgqb3QnA/eKbnb+OTYYALM5L5sm7iIhIExwOB/bt24dJkyadcWDjMT8+VlBQ0K/Pv7uiCYs2l7guvHn6ZNt5m83uxKLNJdhd0dSv6+mN/m7j79hHjH3E2EeMfeRpsQ2HHw2x2uxYurUUktR16DmdBECSgKVbS2G12X2xPCIiIlXg8ONj4eH9d66XLUUNrnd8eqLzHaAtRQ39tqbe6M82WsA+Yuwjxj5i7CNPi204/PhYcHBwvzyvJEnYuL+uV5dc6LRxfx3UcOhXf7XRCvYRYx8x9hFjH3labMPhx8cqKir65XmbbA63T3X1lASgytIOi035q433VxutYB8x9hFjHzH2kafFNhx+NKK1o2/DS0sfH09EROQvOPz4WGJiYr88b9Cgvn0sP7iPj/eG/mqjFewjxj5i7CPGPvK02IbDj4+dOHGiX543zKRHrNnY44ttdtLh1IkPzSblh5/+aqMV7CPGPmLsI8Y+8rTYhsOPjzU19c95dXQ6HablRHr02Gk5kaq45EV/tdEK9hFjHzH2EWMfeVpsw+HHx/rzrNGT08NdF9rsic4Lb05OV8fHGP31jNq+wj5i7CPGPmLsI0+LbTj8+FhmZma/PXeoyYDFecnQ6bpeZfx0OgA6HbA4L1k11/jqzzZawD5i7CPGPmLsI0+LbTj8+NiBAwf69flz48OwbEqq6x2g04egzttMhgAsm5Kqqmt79Xcbf8c+Yuwjxj5i7CNPi23U8Z/8A4gvTiaYGx+GdTNzur2qe8x/r+p+UcZQhBjV9VamGk60qGbsI8Y+Yuwjxj7ytNiGw4+PDR482CffJ9RkwPSRUZiWEwmLzYGWDgeCB+lhNulVcXBzd3zVxl+xjxj7iLGPGPvI02IbDj8+Fhbm291MOp0OYYEGhAWq/x+1r9v4G/YRYx8x9hFjH3labMNjfnysvLxc6SWoFtuIsY8Y+4ixjxj7yNNiGw4/RERENKBw+PGxhIQEpZegWmwjxj5i7CPGPmLsI0+LbTj8+JjFYlF6CarFNmLsI8Y+Yuwjxj7ytNiGw4+PnTx5UuklqBbbiLGPGPuIsY8Y+8jTYhsOPz4iSRIa2+xosJ36/1o8b0JfqfUj+GrBPmLsI8Y+Yuwjz5ttOn8XVltsiv4u1En8LezGbrdjx44dGD16tFeuZ2K12bs92WDsf082ODk9XDWXlyAiIuoPvvhd6HA4sG/fPkyaNAkGg/i5+M5PP9pd0YQb3tmPVbsqUf2Tf9gAUG1px6pdlbjhnf3YXaG9K+Z6orCwUOklqBr7iLGPGPuIsY+8vrZR4+9CDj/9ZHdFExZtLoHN7oQE4PS31zpvs9mdWLS5hAMQTk3tJI99xNhHjH3E2EdeX9qo9Xchh59+YLXZsXRrKSSp6z/o00kAJAlYurUUVpvdF8tTLS2eRdSb2EeMfcTYR4x95HnaRs2/Czn89IMtRQ2uKbcnOqfeLUUN/bks1QsPD1d6CarGPmLsI8Y+Yuwjz9M2av5dyOHHyyRJwsb9dT3+h/1TG/fXDehPgZWVlSm9BFVjHzH2EWMfMfaR50kbtf8u5PDjZU02h9uR7D0lAaiytMNi435nIiLyb2r/Xcjhx8taO/r2D6ylj4/3Z3FxcUovQdXYR4x9xNhHjH3kedJG7b8LOfx4WdCgvp0bKLiPj/dnra2tSi9B1dhHjH3E2EeMfeR50kbtvws5/HhZmEmPWLMRvT0fpg6nTvZkNg3c4aehYWAf8H0m7CPGPmLsI8Y+8jxpo/bfhRx+vEyn02FaTqRHj52WE8lTrBMRkd9T++9CDj/9YHJ6OEyGgB5PvDoAJkMAJqcP7I9ajhgxQuklqBr7iLGPGPuIsY88T9uo+Xchh59+EGoyYHFeMnQ6nPEfug6ATgcszkse8Nf4Ki4uVnoJqsY+Yuwjxj5i7CPP0zZq/l3I4aef5MaHYdmUVNfUe/o/+M7bTIYALJuSitx4nl20o6ND6SWoGvuIsY8Y+4ixj7y+tFHr78KB/VZDP8uND8O6mTndXsk25r9Xsr0oYyhCjAP3IOefCg0NVXoJqsY+Yuwjxj5i7COvr23U+LuQw08/CzUZMH1kFKblRMJic6DB0oxwcwjMJj0Pbj5NZKRnB8cNFOwjxj5i7CPGPvK80eb034UtHQ4ED9Ir9ruQu718RKfTISzQgJa6SoQFGjj4dKO0tFTpJaga+4ixjxj7iLGPPG+26fxdGGM2Kfq7UJXDz4svvoikpCQEBgZi3Lhx+Pbbb2Xvu2HDBuTm5mLw4MEICQnBmDFj8Oabb7rdR5IkLF68GLGxsQgKCkJeXh6Kior6+8cgIiIiFVLd8LN+/XosXLgQS5Yswd69ezF69GhMmTIFtbW13d4/PDwcDz/8MHbu3Invv/8e+fn5yM/Px+bNm133eeqpp/D8889j1apV+OabbxASEoIpU6agra3NVz+Wy7Bhw3z+Pf0F24ixjxj7iLGPGPvI02Ib1Q0/zzzzDObOnYv8/HxkZ2dj1apVCA4Oxpo1a7q9/wUXXIDp06djxIgRSE1Nxd13341Ro0bhq6++AnDqXZ+VK1di0aJFuPLKKzFq1Ci88cYbOHbsGDZu3OjDn+yU9vbeX+htoGAbMfYRYx8x9hFjH3labKOq4ae9vR179uxBXl6e67aAgADk5eVh586dZ3y8JEnYtm0bCgsLcf755wM4ta+yurra7TnPOussjBs3TvicFosFTU1Nrj82m60PP9n/1NfXe+V5tIhtxNhHjH3E2EeMfeRpsY2qPu1VX18Ph8OB6Ohot9ujo6Nx8OBB2cc1NjYiLi4ONpsNer0ef/nLXzB58mQAQHV1tes5Tn/Ozq91Z+TIkWhpaXH9PT8/H3feeSdiY2NRUlLieg5Jkly75NLT01FRUYHW1lYEBgYiISHBdWxRVFQUAgICYLFYUFBQgNTUVFRXV6O5uRkmkwlJSUkoLCwEAERERMBoNOLYsWMAgOTkZNTV1cFqtWLQoEFIS0vDgQMHAJza7RcUFITKykoAQFJSEhoaGtDU1AS9Xo/MzEwcOHAAkiRh8ODBMJvNOHr0KABg+PDhaGpqwsmTJ6HT6TBixAgUFhbC4XAgLCwMQ4YMwZEjRwAA8fHxaGlpcV3jJTs7G4cOHYLdbofZbEZERITroLhhw4bBZrPh+PHjAICsrCwcPnwY7e3tCAkJQXR0NA4fPgwAiI2Nhd1uR11dHSwWC+x2O8rLy9HW1oagoCDExcW5TrDV+c+wpqYGAJCWlobKykpX7+HDh+PQoUMATn06wWAwoKqqCgCQkpKCmpoaNDc3w2g0IiUlxbVNDR06FCaTya13fX09LBYLDAYDMjIyUFBQ4OodHByMiooKAEBiYiJOnDgh2zssLAzl5eUAgISEBFgsFtne4eHhKCsrA3DqKsqtra2u3iNGjEBzczMKCgoQGhqKyMhIt97t7e2uF6jMzEyUlZXBZrMhJCQEMTExrm02JiYGTqfTbZs9evSoq3d8fLzbNqvT6Vy9U1NTUVVVhZaWFphMJiQmJgp719bWwmq1dts7MDCw22329N5DhgxBaGio2zbb2NiIxsZGBAQEICsrCwcPHoTT6URbWxusVqtbb6vVihMnTnTZZrvr3dbW1u02GxoaiqioqG63WQDIyMjAkSNHYLPZEBwc3OfXiM7XJW+/RnS+9vjza0Rn7/54jXA4HDhy5Ihfv0YUFxejo6PD668Rna89an+NkCQJPaWTenPvfnbs2DHExcXh66+/xvjx412333fffdixYwe++eabbh/ndDpx+PBhWK1WbNu2DX/4wx+wceNGXHDBBfj6668xYcIEHDt2DLGxsa7HXHfdddDpdFi/fr3bc9ntduzYsQMpKSkICPjfG2Mmkwkmk6nPP6PD4YBez/P6dIdtxNhHjH3E2EeMfeT5SxuHw4F9+/Zh0qRJMBjE7+2oardXREQE9Hq9a4rsVFNTg5iYGNnHBQQEIC0tDWPGjMHvfvc7XHPNNVi+fDkAuB7X2+c0m80ICwtz/fHG4APANbVTV2wjxj5i7CPGPmLsI0+LbVQ1/BiNRowdOxbbtm1z3eZ0OrFt2za3d4LOxOl0uo7RSU5ORkxMjNtzNjU14ZtvvunVc3qLt44d0iK2EWMfMfYRYx8x9pGnxTaqOuYHABYuXIjZs2cjNzcX5557LlauXInm5mbk5+cDAGbNmoW4uDjXOzvLly9Hbm4uUlNTYbPZ8Nlnn+HNN9/ESy+9BODUCZUWLFiAZcuWIT09HcnJyXjkkUcwbNgwTJs2zec/X0hIiM+/p79gGzH2EWMfMfYRYx95WmyjuuFnxowZqKurw+LFi1FdXY0xY8Zg06ZNrgPZysvL3Y7FaW5uxu23346KigoEBQUhKysLb731FmbMmOG6z3333Yfm5mbceuutOHnyJM477zxs2rQJgYGBPv/5RLvaBjq2EWMfMfYRYx8x9pGnxTaqOuBZDToPeB49enS/HOBVUFCA7Oxsrz+vFrCNGPuIsY8Y+4ixjzx/aeO3BzwTERER9TcOPz6mxbcPvYVtxNhHjH3E2EeMfeRpsQ2HHx9zOp1KL0G12EaMfcTYR4x9xNhHnhbbcPjxMbkLtBLbnAn7iLGPGPuIsY88Lbbh8ENEREQDCocfH0tPT1d6CarFNmLsI8Y+Yuwjxj7ytNiGw4+PdV4wkLpiGzH2EWMfMfYRYx95WmzD4cfH2tralF6CarGNGPuIsY8Y+4ixjzwttuHw42NBQUFKL0G12EaMfcTYR4x9xNhHnhbbcPjxsfj4eKWXoFpsI8Y+Yuwjxj5i7CNPi204/PhYUVGR0ktQLbYRYx8x9hFjHzH2kafFNhx+iIiIaEDh8ONjUVFRSi9BtdhGjH3E2EeMfcTYR54W23D48TGdTqf0ElSLbcTYR4x9xNhHjH3kabENhx8fq6mpUXoJqsU2Yuwjxj5i7CPGPvK02IbDDxEREQ0oHH58LDU1VeklqBbbiLGPGPuIsY8Y+8jTYhsOPz5WVVWl9BJUi23E2EeMfcTYR4x95GmxDYcfH2tpaVF6CarFNmLsI8Y+Yuwjxj7ytNiGw4+PmUwmpZegWmwjxj5i7CPGPmLsI0+LbTj8+FhiYqLSS1AtthFjHzH2EWMfMfaRp8U2HH587NChQ0ovQbXYRox9xNhHjH3E2EeeFttw+CEiIqIBhcOPj0VGRiq9BNViGzH2EWMfMfYRYx95WmzD4cfHDAaD0ktQLbYRYx8x9hFjHzH2kafFNhx+fEyL50vwFrYRYx8x9hFjHzH2kafFNhx+iIiIaEDh8ONjKSkpSi9BtdhGjH3E2EeMfcTYR54W23D48bHa2lqll6BabCPGPmLsI8Y+YuwjT4ttOPz4mNVqVXoJqsU2Yuwjxj5i7CPGPvK02IbDj48ZjUall6BabCPGPmLsI8Y+YuwjT4ttOPz4mBb3nXoL24ixjxj7iLGPGPvI02IbDj8+dvDgQaWXoFpsI8Y+Yuwjxj5i7CNPi204/BAREdGAwuHHx4YOHar0ElSLbcTYR4x9xNhHjH3kabENhx8fCwwMVHoJqsU2Yuwjxj5i7CPGPvK02IbDj49VVlYqvQTVYhsx9hFjHzH2EWMfeVpsw+GHiIiIBpQeXar1n//8p0dPfv7553v0OC1LSkpSegmqxTZi7CPGPmLsI8Y+8rTYpkfDzwUXXACdTuf6uyRJbn+X43A4PF+ZRjU0NCA4OFjpZagS24ixjxj7iLGPGPvI02KbHg0/X3zxhdvfbTYb7rvvPrS0tODWW29FZmYmgFPnAli9ejVCQkLw1FNPeX+1GtDU1KT0ElSLbcTYR4x9xNhHjH3kabFNj4afSZMmuf194cKFMBqN2LVrl9tR4FOnTsUdd9yBSZMmYdOmTZg8ebJ3V6sBBkOPkg9IbCPGPmLsI8Y+YuwjT4ttPDrg+e2338avf/3rbj/+FhwcjF//+td46623+rw4LcrIyFB6CarFNmLsI6bFPpIkobHNjmqLDY1tdkiS5PFzabGPN7GPPC228Wj4aW5uRlVVlezXq6qq0NLS4vGitKygoEDpJagW24ixj5iW+lhtdnz4Yy3mvFeAa9/6AbPWn/r/c94rwIc/1sJqs/f6ObXUpz+wjzwttvFo+MnLy8Nzzz2HDRs2dPnaBx98gOeeew55eXl9XhwR0UCzu6IJN7yzH6t2VaLa0u72tWpLO1btqsQN7+zH7grtHYdB5Cse7ch78cUXceGFF+Laa69FbGws0tLSAAAlJSU4duwYUlNT8cILL3h1oVoxZMgQpZegWmwjxj5iWuizu6IJizaXQJKA7nZwdd5mszuxaHMJlk1JRW58WI+eWwt9+hP7yNNiG4/e+YmLi8O+ffvwzDPPYOTIkaipqUFNTQ1ycnLw7LPPYt++fYiPj/f2WjUhNDRU6SWoFtuIsY+Yv/ex2uxYurVUdvD5KQmAJAFLt5b2eBeYv/fpb+wjT4ttPD7Dc2BgIO6++25s2rQJBw4cwIEDB7Bp0ybcddddCAoK8uYaNeXo0aNKL0G12EaMfcT8vc+WogbY7M4zDj6dJJx6B2hLUUOP7u/vffob+8jTYps+Xd7CZrNh586d+Oijj1BfX++tNRERDSiSJGHj/roeDz4/tXF/XZ8+BUY0EHk8/Dz//POIjY3FhAkTcNVVV+H7778HANTX1yMiIgJr1qzx2iK1ZPjw4UovQbXYRox9xPy5T5PNgarTDm7uCQlAlaUdFtuZz6bvz318gX3kabGNR8PP2rVrsWDBAlx88cVYs2aN2391RERE4MILL8S7777rtUVqSWNjo9JLUC22EWMfMX/u09rRt0sBtfTg8f7cxxfYR54W23g0/Dz99NO48sorsW7dOkydOrXL18eOHYv9+/d7vKgXX3wRSUlJCAwMxLhx4/Dtt9/K3nf16tWYOHEihgwZgiFDhiAvL6/L/efMmQOdTuf25+KLL/Z4fX2hxY3IW9hGjH3E/LlP0CB9nx4f3IPH+3MfX2AfeVps49HwU1xcjEsuuUT26+Hh4Th+/LhHC1q/fj0WLlyIJUuWYO/evRg9ejSmTJmC2trabu+/fft2zJw5E1988QV27tyJhIQEXHTRRaisrHS738UXX4yqqirXn3feecej9fVVQECfDrPSNLYRYx8xf+4TZtIj1mzEmS8X7U4HINZshNl05uHHn/v4AvvI02Ibj36iwYMHCw9wLigoQExMjEcLeuaZZzB37lzk5+cjOzsbq1atQnBwsOwxRG+//TZuv/12jBkzBllZWXj11VfhdDqxbds2t/uZTCbExMS4/ih13oKsrCxFvq8/YBsx9hHz5z46nQ7TciI9euy0nEjodGcem/y5jy+wjzwttvFo+Ln00kvxyiuv4OTJk12+tn//fqxevRpXXHFFr5+3vb0de/bscTs7dEBAAPLy8rBz584ePUdLSws6OjoQHh7udvv27dsRFRWFzMxMzJs374zvTFksFjQ1Nbn+2Gy2Xv883Tl48KBXnkeL2EaMfcT8vc/k9HCYDAE9fvdHB8BkCMDk9PAz3hfw/z79jX3kabGNR2d4XrZsGcaNG4eRI0di6tSp0Ol0eP3117FmzRp88MEHiI2NxeLFi3v9vPX19XA4HIiOjna7PTo6usfx77//fgwbNsxtgLr44otx1VVXITk5GSUlJXjooYdwySWXYOfOndDru3+7eOTIkW7XJ8vPz8edd96J2NhYlJSUuNYlSZJrl1x6ejoqKirQ2tqKwMBAJCQkoKioCAAQFRWFgIAANDY2oqCgAKmpqaiurkZzczNMJhOSkpJQWFgI4NRB40ajEceOHQMAJCcno66uDlarFYMGDUJaWhoOHDgA4NQuxqCgINduvqSkJDQ0NKCpqQl6vR6ZmZk4cOAAJEnC4MGDYTabXedsGD58OJqamnDy5EnodDqMGDEChYWFcDgcCAsLw5AhQ3DkyBEAQHx8PFpaWtDQcOqcItnZ2Th06BDsdjvMZjMiIiJQWloKABg2bBhsNptrwMzKysLhw4fR3t6OkJAQREdH4/DhwwCA2NhY2O121NXVwWKxwG63o7y8HG1tbQgKCkJcXByKi4tdvQGgpqYGAJCWlobKykpX7+HDh+PQoUMAgMjISBgMBtc16FJSUlBTU4Pm5mYYjUakpKS4tqmhQ4fCZDK59a6vr4fFYoHBYEBGRobr2jbh4eEIDg5GRUUFACAxMREnTpyQ7R0WFoby8nIAQEJCAiwWi2zv8PBwlJWVATh1ItHW1lZX7xEjRsBisaCgoAChoaGIjIx0693e3u56NzYzMxNlZWWw2WwICQlBTEyMa5uNiYmB0+l022aPHj3q6h0fH++2zep0Olfv1NRU13X7TCYTEhMThb1ra2thtVq77R0YGNjtNnt67yFDhiA0NNRtm21sbERjYyMCAgKQlZWFgwcPwul0oq2tDVar1a231WrFiRMnumyz3fVua2vrdpsNDQ1FVFRUt9sscOqij0eOHIHNZkNwcHCfXiPuGReBJ/5VCx3EJzrUAdDpgDkZepSXHOrRa0Tna48/v0Z09u6P1wiHw4EjR4749WtEcXExOjo6vP4a0fnao/bXiN6c8kEneXiCiNraWjz00EPYsGGD6x0gs9mMq6++Gk888QSioqJ6/ZzHjh1DXFwcvv76a4wfP951+3333YcdO3bgm2++ET7+iSeewFNPPYXt27dj1KhRsvc7fPgwUlNTsXXrVvzqV79y+5rdbseOHTuQkpLitp/TZDLBZDL1+mc6XWVlJeLi4vr8PFrENmLsI6aVPrsrmrB0aylsdicA9yGo810hkyEAi/OSe3xpC0A7ffoL+8jzlzYOhwP79u3DpEmTYDCI39vx6J0f4NTE9+qrr+LVV19FXV0dnE4nIiMj+3RgVEREBPR6vWuK7FRTU3PGY4hWrFiBJ554Alu3bhUOPsCpiTMiIgLFxcVdhp9OZrNZ9l2hvjjrrLO8/pxawTZi7COmlT658WFYNzMHW4oasHF/ndv5f2LMRkzLicRFGUMRYuzd65NW+vQX9pGnxTYeTSqfffYZHI7/nVciMjIS0dHRfT4i3Gg0YuzYsW4HK3cevPzTd4JO99RTT+EPf/gDNm3ahNzc3DN+n4qKChw/fhyxsbF9Wq8nOt/epK7YRox9xLTUJ9RkwPSRUXjtumy8f9PZeGPGqf//2nXZmD4yqteDD6CtPv2BfeRpsY1H08rll1+O6Oho3Hrrrdi2bRucTqfXFrRw4UKsXr0ar7/+Og4cOIB58+ahubkZ+fn5AIBZs2bhwQcfdN3/ySefxCOPPII1a9YgKSkJ1dXVqK6uhtVqBQBYrVb8/ve/x65du1BWVoZt27bhyiuvRFpaGqZMmeK1dRMReZtOp0NYoAExZhPCAg09+lQXEZ2ZR8PP3//+d1xxxRV4//33cdFFFyE2NhZ33HEHvvzyyz4vaMaMGVixYgUWL16MMWPG4LvvvsOmTZtcB7KVl5e7DpgCgJdeegnt7e245pprEBsb6/qzYsUKAIBer8f333+PK664AhkZGbj55psxduxYfPnll145hqe3EhISfP49/QXbiLGPGPuIsY8Y+8jTYhuPD3gGgI6ODmzevBnr16/HJ598AovFgtjYWFxzzTWYMWOGcFeVWnUe8Dx69Oh+OeanqqpKkd1t/oBtxNhHjH3E2EeMfeT5S5veHPDcp4N0Bg0ahMsvvxxvvvkmamtrsWHDBpx//vl49dVXMXHixL48tWZ1fuyWumIbMfYRYx8x9hFjH3labOO1c1ZbrVbU1taipqYGbW1tvfq8PREREZGv9Gm3V2NjIzZs2ID169fjiy++QEdHB84++2zMmDEDM2bMQGpqqjfX6hP9vduLiDwnSRKabA60djgQNEiPMJOeBwETEQAfnOfnzTffxHvvvYctW7agvb0dWVlZeOihhzBjxgxNXgPEmw4dOoSMjAyll6FKbCM2kPtYbfZuz3sT+9/z3kxOD8exI4cHbJ+eGMjbT0+wjzwttvFo+Jk9ezZSUlLwu9/9DjNmzDjjSQXpf+x2u9JLUC22ERuofU4/4/FPVVvasWpXJdbursLsdD009vrsVQN1++kp9pGnxTa9Hn4kScK///1vZGdnIygoqD/WpGlhYT0/Hf1AwzZiA7HP7oomLNpcAknq/lpXnbfZ7E6sPuBEYmJTry75MJAMxO2nN9hHnhbb9PqA5/b2dowbNw5//vOf+2M9mnf61ebpf9hGbKD1sdrsWLq1VHbw+Snpv3+Wbi2F1aa9/0r1hoG2/fQW+8jTYpteDz8mkwkxMTEwGo39sR7N67wiL3XFNmIDrc+WogbY7M4zDj6dJJx6B2hLUUN/LstvDbTtp7fYR54W23j0Ufc5c+bgjTfeQHt7+5nvTETUS5IkYeP+uh4PPj+1cX8dT7VBREIeHfB89tlnY+PGjcjJycGcOXOQlJTU7fE/V111VZ8XqDVxcXFKL0G12EZsIPVpsjncPtXVUxKAKks7LDYHwgI9ennTrIG0/XiCfeRpsY1Hrw4zZ850/e9HHnmk2/vodDq3K7/TKW1tbTjrrLOUXoYqsY3YQOrT2tG3146WDg4/pxtI248n2EeeFtt49OrwxRdfeHsdA8bx48ddF2kld2wjNpD6BA3q2wlGg/v4eC0aSNuPJ9hHnhbbeDT8TJo0ydvrICJyCTPpEWs2otrS3qvjfnQAYsxGmE0cfohIXp+u7WWz2bBz50589NFHqK+v99aaNI1nwJbHNmIDqY9Op8O0nEiPHjstJ5KXvOjGQNp+PME+8rTYxuPh5/nnn0dsbCzOO+88XHXVVfj+++8BAPX19YiIiMCaNWu8tkgtOXz4sNJLUC22ERtofSanh8NkCEBPxxgdAJMhAJPTtXdOEm8YaNtPb7GPPC228Wj4Wbt2LRYsWICLL74Y/+///T+3j5VGRETgwgsvxLvvvuu1RWoJTw8gj23EBlqfUJMBi/OSodPhjANQ59cX5yUj1MQDnbsz0Laf3mIfeVps49Hw8/TTT+PKK6/EunXrMHXq1C5fHzt2LPbv39/nxWlRaGio0ktQLbYRG4h9cuPDsGxKqusdoNOHoM7bTIYALPiZmZe2EBiI209vsI88Lbbx6D+RiouLcdddd8l+PTw8HMePH/d4UVoWFRWl9BJUi23EBmqf3PgwrJuZ0+1V3WP+e1X3izKGQu/sUHCV6jdQt5+eYh95Wmzj0fAzePBg4QHOBQUFiImJ8XhRWnb48GFkZ2crvQxVYhuxgdwn1GTA9JFRmJYTCYvNgZYOB4IH6WE26V0HNxcUFA7YPj0xkLefnmAfeVps49Fur0svvRSvvPIKTp482eVr+/fvx+rVq3HFFVf0dW1ERG50Oh3CAg2IMZsQFmjgp7qIyCMeDT/Lli2Dw+HAyJEjsWjRIuh0Orz++uu46aabkJubi6ioKCxevNjba9WE2NhYpZegWmwjxj5i7CPGPmLsI0+LbTwafoYNG4Y9e/bg4osvxvr16yFJEt5880188sknmDlzJnbt2oWIiAhvr1UT7Ha70ktQLbYRYx8x9hFjHzH2kafFNh6f5ycqKgqvvvoqGhoaUFNTg6qqKpw4cQJr1qzR5MFR3lJXV6f0ElSLbcTYR4x9xNhHjH3kabGNV06IERl56kys7e3t6OjoQEhIiDeeloiIiMjrPHrn591338U999zjdttjjz2G0NBQDB48GNOnT4fVavXKArUmIyND6SWoFtuIsY8Y+4ixjxj7yNNiG49Pctjc3Oz6+9dff43HHnsMU6ZMwT333INNmzbhj3/8o9cWqSVHjhxRegmqxTZi7CPGPmLsI8Y+8rTYxqPdXiUlJZg9e7br7+vWrUNMTAw+/PBDGAwGOJ1OfPDBB1i+fLnXFqoVNptN6SWoFtuIsY8Y+4ixj5hSfSRJQpPNgdYOB4IG6RH2k3NXqYUWtx2Phh+bzYbAwEDX3z///HNccsklMBhOPV12djb+8pe/eGeFGhMcHKz0ElSLbcTYR4x9xNhHzNd9rDZ7t2ctj/3vWcsnp4er5jp1Wtx2PNrtlZycjK1btwIAdu/ejeLiYlx88cWur9fU1GjyWiDeoMXzJXgL24ixjxj7iLGPmC/77K5owg3v7MeqXZWotrhfNLTa0o5Vuypxwzv7sbuiyWdrEtHituPR8PPb3/4W7733HkaNGoWLLroI8fHxuPzyy11f/9e//oWcnByvLVJLSkpKlF6CarGNGPuIsY8Y+4j5qs/uiiYs2lwCm90JCYB02tc7b7PZnVi0uUQVA5AWtx2Php8777wTL7/8MlJTU3HllVfi888/R1BQEACgoaEB1dXVuPHGG726UCIiIn9mtdmxdGspJKnr0HM6CYAkAUu3lsJq095JBpXm8Q7FuXPnYu7cuV1uDw8Px+7du/u0KC2Ljo5WegmqxTZi7CPGPmLsI+aLPluKGlzv+PRE5ztAW4oaMH2kcicP1uK206ejqRwOB/bs2YOysjIAQFJSEsaOHQu9Xu+NtWmSJPV0sx942EaMfcTYR4x9xPq7jyRJ2Li/rseDz09t3F+HaTmRin0KTIvbjseXt3jttdcQHx+P8ePH4/rrr8f111+P8ePHIy4uDmvWrPHmGjWltrZW6SWoFtuIsY8Y+4ixj1h/92myOdw+1dVTEoAqSzssNof3F9VDWtx2PHrn5+WXX8a8efMwZswYPProo66zPxYWFuLll1/G3Llz0d7ejttuu82riyUiIvJHrR19G15aOhwIC1THR9+1QCd58H5WSkoKEhISsHXrVgwaNMjtax0dHbjwwgtRWVmJw4cPe22hvmK327Fjxw6MHj26X3bfdXR0dGlGp7CNGPuIsY8Y+4j1d5/GNjuufesHjx///k1nKzb8+Mu243A4sG/fPkyaNMl13kE5Hu32qq6uxnXXXddtjEGDBuH6669HTU2NJ0+teRUVFUovQbXYRox9xNhHjH3E+rtPmEmPWLMRvT1qR4dTJz40m5Q7llaL245Hw88555yDQ4cOyX790KFDGDNmjKdr0rTW1lall6BabCPGPmLsI8Y+Yv3dR6fTYVpOpEePVfJgZ0Cb245Hw88LL7yA9957D88995xblNbWVjz77LN477338Oc//9lri9SSn14WhNyxjRj7iLGPGPuI+aLP5PRwmAwBPX73RwfAZAjA5PTw/lzWGWlx2+nRMT+jRo3qcltDQwOqqqpgMBgwbNgwAMCxY8dgt9sRGxuLoUOHYt++fd5fcT/jMT/KYRsx9hFjHzH2EfNVn84zPJ/pRIc6ADodsGxKKnLjw/p9XSL+su305pifHh09FR4e3uUtt6FDhyI9Pd3ttqSkpN6tdAAqKipCdna20stQJbYRYx8x9hFjHzFf9cmND8OyKalYurUUNrsTgPsQ1Pmb1mQIwOK8ZMUHH0Cb206Php/t27f38zKIiIgGhtz4MKybmdPtVd1j/ntV94syhiLEyBMG9xeeNMDHoqKUO0W52rGNGPuIsY8Y+4j5uk+oyYDpI6MwLScSFpsDLR0OBA/Sw2zSK3pwc3e0uO14PPw4HA689dZb+PTTT3HkyBEAQGJiIi6//HLceOONvMSFjIAAj0+qrXlsI8Y+Yuwjxj5iSvXR6XQICzSo+gSGWtx2PPqJGhsbMWHCBPzmN7/B559/jo6ODnR0dGDLli3Iz8/Heeedh6amJm+vVROqq6uVXoJqsY0Y+4ixjxj7iLGPPC228Wj4efjhh7Fnzx688MILqKurw969e7F3717U1tbiz3/+M3bv3o2HH37Y22slIiIi6jOPLm8RFxeHa665Bs8991y3X7/rrrvw/vvv49ixY31eoK/190fdbTYbTCaT159XC9hGjH3E2EeMfcTYR56/tOn3y1scP34cmZmZsl/PyspCQ0ODJ0+teVp8+9Bb2EaMfcTYR4x9xNhHnhbbeDT8pKWl4eOPP5b9+scff4zU1FSPF6Vlzc3NSi9BtdhGjH3E2EeMfcTYR54W23g0/Nx+++34/PPPcemll+Lzzz9HWVkZysrKsHnzZlx22WXYsmUL5s+f7+21aoI/vHWoFLYRYx8x9hFjHzH2kafFNh4PP4sXL8Y//vEPXHLJJUhNTUVqaiouvfRSbNu2DYsXL8a8efM8XtSLL76IpKQkBAYGYty4cfj2229l77t69WpMnDgRQ4YMwZAhQ5CXl9fl/pIkYfHixYiNjUVQUBDy8vJQVFTk8fr6gmfBlsc2Yuwjxj5i7CPGPvK02MbjD+8/+uijqKiowNtvv43HH38cjz/+ON5++21UVFRgyZIlHi9o/fr1WLhwIZYsWYK9e/di9OjRmDJlCmpra7u9//bt2zFz5kx88cUX2LlzJxISEnDRRRehsrLSdZ+nnnoKzz//PFatWoVvvvkGISEhmDJlCtra2jxep6cKCwt9/j39hTfaSJKExjY7qi02NLbZ4cHx/KrFbUeMfcTYR4x95GmxTZ/OqhQREYHrr7/eW2sBADzzzDOYO3cu8vPzAQCrVq3Cp59+ijVr1uCBBx7ocv+3337b7e+vvvoqPvjgA2zbtg2zZs2CJElYuXIlFi1ahCuvvBIA8MYbbyA6OhobN270+vpJGVabvdtTxcf+91Txk9PDEWpS70nEiIjId1R12sb29nbs2bMHeXl5rtsCAgKQl5eHnTt39ug5Wlpa0NHRgfDwcABAaWkpqqur3Z7zrLPOwrhx43r8nN4UERHh8+/pLzxts7uiCTe8sx+rdlWi+ieDDwBUW9qxalclbnhnP3ZX+PeJN7ntiLGPGPuIsY88LbZR1fBTX18Ph8OB6Ohot9ujo6N7/FG7+++/H8OGDXMNO52P6+1zWiwWNDU1uf7YbLbe/CiyjEajV55Hizxps7uiCYs2l8Bmd0KC+9WR8d+/SwBsdicWbS7x6wGI244Y+4ixjxj7yNNiG03tB3jiiSfw7rvvYvv27QgMDOzTc40cORItLS2uv+fn5+POO+9EbGwsSkpKAJwaoCRJch2PlJ6ejoqKCrS2tiIwMBAJCQmuA6ujoqIQEBCAoqIimM1mpKamorq6Gs3NzTCZTEhKSnLtV42IiIDRaHSdJDI5ORl1dXWwWq0YNGgQ0tLScODAAQBAeHg4goKCXMc4JSUloaGhAU1NTdDr9cjMzMSBAwcgSRIGDx4Ms9mMo0ePAgCGDx+OpqYmnDx5EjqdDiNGjEBhYSEcDgfCwsIwZMgQ13Xb4uPj0dLS4jp/U3Z2Ng4dOgS73Q6z2YyIiAiUlpYCAIYNGwabzYbjx48DOHXep8OHD6O9vR0hISGIjo7G4cOHAQCxsbGw2+2oq6uDxWLB2LFjUV5ejra2NgQFBSEuLg7FxcWu3gBQU1MDAIgZnozHtpRAkroOPaeT/vt/Hv28BI+PD8XIzFQcPHgQADB06FCYTCa33vX19bBYLDAYDMjIyEBBQYGrd3BwMCoqKgCcup7diRMnZHuHhYWhvLwcAJCQkACLxSLbOzw8HGVlZQBOnUi0tbXV1XvEiBEoKipCSEgIQkNDERkZ6da7vb0d9fX1AIDMzEyUlZXBZrMhJCQEMTExrm02JiYGTqfTbZs9evSoq3d8fLzbNqvT6Vy9U1NTUVVVhZaWFphMJiQmJuLQoUMAgMjISBgMBlRVVQEAUlJSUFtbC6vVCqPRiJSUFLfegYGB3W6zp/ceMmQIQkND3bbZxsZGNDY2IiAgAFlZWTh48CCcTifa2tqQkZHh1ttqteLEiRNdttnuere1tXW7zYaGhiIqKqrbbRYAMjIycOTIEdhsNgQHB/f5NaLzP8q8/RpRWFgIs9ns168Rnb17+hqRlpaGyspKV+/hw4fLbrMOhwONjY1obm7udpv1h9eI4uJidHR0eP01ovO1R+2vEb05xtOjMzz3l/b2dgQHB+P999/HtGnTXLfPnj0bJ0+exEcffST72BUrVmDZsmXYunUrcnNzXbcfPnwYqamp+M9//oMxY8a4bp80aRLGjBnT5SzVnWd4TklJcbuYm8lk8srH/QoKCpCdnd3n59Gi3rb58MdarNpVecbB56d0AG77RRymj/S/qxRz2xFjHzH2EWMfef7Spl/P8Gyz2fDxxx/j+++/93iBcoxGI8aOHYtt27a5bnM6ndi2bRvGjx8v+7innnoKf/jDH7Bp0ya3wQc4NaHHxMS4PWdTUxO++eYb4XOazWaEhYW5/njrPAfJycleeR4t6k0bSZKwcX9drwafThv31/nlp8C47Yixjxj7iLGPPC226fXwYzQace211+Lrr7/uj/Vg4cKFWL16NV5//XUcOHAA8+bNQ3Nzs+vTX7NmzcKDDz7ouv+TTz6JRx55BGvWrEFSUhKqq6tRXV0Nq9UKANDpdFiwYAGWLVuGjz/+GD/88ANmzZqFYcOGub275Cudb91SV71p02RzuH2qq6ckAFWWdlhsjl4/VmncdsTYR4x9xNhHnhbb9PqYH51Oh/T0dNe+Q2+bMWMG6urqsHjxYlRXV2PMmDHYtGmTa19ueXm52+6ol156Ce3t7bjmmmvcnmfJkiV49NFHAQD33Xcfmpubceutt+LkyZM477zzsGnTpj4fF+SJzqGMuupNm9aOvg0vLR0OhAX61yFv3HbE2EeMfcTYR54W23j06v/QQw9h4cKFuPbaa4UXOPXU/PnzZS+PsX37dre/dx78JaLT6bB06VIsXbrUC6vrm0GDBim9BNXqTZugQfo+fa/gPj5eCdx2xNhHjH3E2EeeFtt4NPzs2rULQ4cOxciRI3HBBRcgKSkJQUFBbvfR6XRdDiamU58+oO71pk2YSY9YsxHVlvZeH/AcYzbCbPK/4Yfbjhj7iLGPGPvI02Ibj4afP//5z67//dMDiX+Kw0/3Dhw44BdHzSuhN210Oh2m5URi1a7KM9/5NNNyIqHT6Xr9OKVx2xFjHzH2EWMfeVps49FJDp1O5xn/OBz+d0Ap+ZfJ6eEwGQLQ0zFGB8BkCMDk9PD+XBYREamcqs7wPBB0XnaDuuptm1CTAYvzkqHT4YwDkA6ATgcszkv222t8cdsRYx8x9hFjH3labNOn4WfXrl1Yvnw57rnnHtcZH1taWrB3715NHh3uDacfG0X/40mb3PgwLJuS6noH6PQhqPM2kyEAy6akIjc+zAsrVQa3HTH2EWMfMfaRp8U2Hg0/7e3tuOqqqzBhwgQ8/PDDeP75512nQw8ICMBFF13E431kdJ6um7rytE1ufBjWzczBbb+IQ4zZ/Ro0MWYjbvtFHN65YaRfDz4At50zYR8x9hFjH3labOPR+/+PPPII/va3v+Gll17CL3/5S7ePuwcGBuLaa6/FRx99hIcffthrCyUSCTUZMH1kFKblRMJic6Clw4HgQXqYTXq/PLiZiIj6j0fv/LzzzjuYN28ebr311m73BY4YMcJ1UTpyl5SUpPQSVMsbbXQ6HcICDYgxmxAWaNDU4MNtR4x9xNhHjH3kabGNR8NPbW0tzj77bNmv6/V6tyui0/90XoGXumIbMfYRYx8x9hFjH3labOPR8JOQkOC69Hx3/vWvf2nypEje0NTUpPQSVIttxNhHjH3E2EeMfeRpsY1Hw88NN9yAl19+GTt37nTd1rl7YfXq1Xjvvfcwa9Ys76xQY/R6/zuzsK+wjRj7iLGPGPuIsY88LbbRSZLUm6sDADj1aa+pU6fiH//4B0aMGIH9+/fj7LPPRkNDAyoqKnDppZfio48+8stgdrsdO3bswOjRo/1y/URERAORw+HAvn37MGnSJBgM4s9zefTOj9FoxKZNm7B27VqkpKQgKysLNpsNo0aNwmuvvYZPPvmEg4OMAwcOKL0E1WIbMfYRYx8x9hFjH3labOPxqW51Oh1uuukm3HTTTd5cj+Z58EbbgME2Yuwjxj5i7CPGPvK02Majd35efvllTU6CvjB48GCll6BabCPGPmLsI8Y+YuwjT4ttPHrnZ968edDpdAgPD8d5552HiRMnYuLEiRg7diwCAni5MBGz2az0ElSLbcTYR4x9xNhHjH3kabGNR5NKdXU11q9fjxtvvBHl5eW477778Itf/AKDBw/GRRddhD/84Q/Yvn27l5eqDZ2XAaGu2EaMfcTYR4x9xNhHnhbbePTOT1RUFK655hpcc801AACLxYKvv/4aX375Jd5//308+uij0Ol0sNvtXl0sERERUV95fMBzp5KSEnz55Zf48ssv8c9//hMlJSUICQnB+PHjvbE+zRk+fLjSS1AtthFjHzH2EWMfMfaRp8U2Hg0/f/7zn10DT3V1tevYn3nz5mHixIn42c9+xo+6y2hqakJoaKjSy1AlthFjHzH2EWMfMfaRp8U2Hg0/d911F/R6Pa6++mr8/ve/x9ixY729Ls06efIkhg0bpvQyVIltxNhHjH3E2EeMfeRpsY1HBzzfcccdGDlyJN5//31MmDAB5513Hh588EF89tlnaGxs9PYaNUVLVxn3NrYRYx8x9hFjHzH2kafFNh5d3qJTU1MTvvrqK3z11Vf48ssv8e9//xt2ux0jR47ExIkT8cILL3hzrT7By1sQERH5n36/vEWnsLAwXHrppXj88cfxxhtv4Pnnn0d6ejq+//57/OUvf+nLU2tWYWGh0ktQLbYRYx8x9hFjHzH2kafFNh5/2qugoMD1Ca8vv/wSlZWVAIBhw4bh+uuvx8SJE722SC1xOBxKL0G12EaMfcTYR4x9xNhHnhbbeDT8RERE4MSJE5AkCVlZWbjkkktcZ3pOSkry8hK1JSwsTOklqBbbiLGPGPuIsY8Y+8jTYhuPhp/Zs2dj4sSJOO+88xAREeHtNWnakCFDlF6CarGNGPuIsY8Y+4ixjzwttvHomJ+nn34a06ZN4+DjgSNHjii9BNViGzH2EWMfMfYRYx95WmzTpzM879ixA59++qkrTGJiIi677DJMmjTJK4sjIiIi8jaPhp/29nbMnDkTGzduhCRJrsvdnzx5Ek8//TSmT5+Od955B4MGDfLmWjUhPj5e6SWoFtuIsY8Y+4ixjxj7yNNiG492ez322GP48MMP8bvf/Q5VVVVoaGhAQ0MDqqurce+992LDhg1YunSpt9eqCS0tLUovQbXYRox9xNhHjH3E2EeeFtt4NPysW7cOs2fPxlNPPYXo6GjX7VFRUXjyyScxa9YsvPnmm15bpJY0NDQovQTVYhsx9hFjHzH2EWMfeVps49HwU1VVhXHjxsl+fdy4caiurvZ4UURERET9xaPhJz4+Htu3b5f9+o4dOzS5j9AbsrOzlV6CarGNGPuIsY8Y+4ixjzwttvFo+Jk9ezbee+893HbbbSgsLITD4YDT6URhYSHmzZuHv/71r5gzZ46Xl6oNhw4dUnoJqsU2YmrqI0kSGtvsqLbY0NhmRx8uEeg1auqjRuwjxj7ytNjGo097PfTQQygpKcErr7yC1atXIyDg1AzldDohSRJmz56Nhx56yKsL1Qq73a70ElSLbcTU0Mdqs2NLUQM27q9DlaXddXus2YhpOZGYnB6OUFOfzqDhMTX0UTP2EWMfeVps49GrlF6vx2uvvYaFCxfi008/RXl5OYBT5/m59NJLMWrUKK8uUkvMZrPSS1AtthFTus/uiiYs3VoKm93Z5WvVlnas2lWJtbursDgvGbnxvj8dvtJ91I59xNhHnhbb9Gr4aWtrw0cffYTS0lJERETgsssuw4MPPthfa9MknhVbHtuIKdlnd0UTFm0ugSQB3e3g6rzNZndi0eYSLJuS6vMBiNuPWG/6SJKEJpsDrR0OBA3SI8ykh06n68fVKY/bjzwttunx8FNbW4v/+7//Q2lpqWv/fnBwMDZu3Ii8vLx+W6DWlJaWavLgMW9gGzGl+lhtdizdWio7+PyU9N//s3RrKdbNzPHpLjBuP2I96aPm3Zr9jduPPC226fEBz3/4wx9QVlaGe+65B3/729+wcuVKBAUF4be//W1/ro+IFLalqAE2u/OMg08nCafeAdpSpL1zg2jZ7oom3PDOfqzaVYnqnww+wP92a97wzn7srmhSaIVE3tPjEf7zzz/HrFmzsGLFCtdt0dHRuOGGG1BYWIjMzMx+WaDWDBs2TOklqBbbiCnRR5IkbNxf1+PB56c27q/DtJxIn+0u6c8+WtgNJOrjD7s1+xtff+RpsU2Ph5/y8nLcf//9bredd955kCQJNTU1HH56yGazKb0E1WIbMSX6NNkcbrs/ekoCUGVph8XmQFigb3aT9EcfLe0GkuvjL7s1+xtff+RpsU2Pd3vZbDYEBga63db5dy1+DK6/HD9+XOklqBbbiCnRp7XD0afHt/Tx8b3h7T5a2w0k14e7NU/h6488Lbbp1dheVlaGvXv3uv7e2NgIACgqKnJd2f2nfvazn/VtdUSkqKBB+j49PriPj1fKQNkN5E+7NYm8SSf18NSsAQEB3W7kkiR1ub3zNofDd//V5y12ux07duzA6NGjodd7/4Xb6XS6TgpJ7thGTIk+kiRhznsFqLa09+oXpA5AjNmI167L9tkvR2/1sdrsuOGd/T1+N0QHwGQIUP1uoO76NLbZce1bP3j8nO/fdLbPdmv2N77+yPOXNg6HA/v27cOkSZNgMIi3yx5vtWvXru3zwgg4fPgw0tLSlF6GKrGNmBJ9dDodpuVEYtWuyl4/1tfvCnirT192A00fGdXn799fuuvjjd2aWhl++PojT4tterzVzp49uz/XMWC0t/f+4NGBgm3ElOozOT0ca3dX9fqdkMnp4f29NDfe6KPl3UDd9RmouzW7w9cfeVpso/73sTQmJCRE6SWoFtuIKdUn1GTA4rxk6HSnBhsRHQCdDlicl+zzXUDe6OONT7epVXd9wkx6xJqNZ/znejodTn3izWzSzvDD1x95WmzD4cfHoqOjlV6CarGNmJJ9cuPDsGxKKkyGgFMDzmlf77zNZAhQ7OBfb/Txp0+39VZ3fTp3a3pCze9yeYKvP/K02IbDj48dPnxY6SWoFtuIKd0nNz4M62bm4LZfxCHGbHT7WozZiNt+EYd3bhip2KeevNFHy7uB5PpMTg93DbU9odRuzf6m9L9faqbFNto4Uo2IfCLUZMD0kVGYlhMJi82Blg4HggfpYfbDMx53p3M3kKefbvPH3UCduzUXbS4BznCiQyV3axJ5k+re+XnxxReRlJSEwMBAjBs3Dt9++63sfffv34+rr74aSUlJ0Ol0WLlyZZf7PProo9DpdG5/srKy+vEnEIuNjVXse6sd24ipqY9Op0NYoAExZhPCAg2qGHy80UfLu4FEffxht2Z/U9O/X2qjxTaqGn7Wr1+PhQsXYsmSJdi7dy9Gjx6NKVOmoLa2ttv7t7S0ICUlBU888QRiYmJknzcnJwdVVVWuP1999VV//QhnxLNhy2MbMfYR81Yfre4GOlMfte/W7G/890ueFtuo6n3LZ555BnPnzkV+fj4AYNWqVfj000+xZs0aPPDAA13u//Of/xw///nPAaDbr3cyGAzC4ciX6urqEBnp2X9Zah3biLGPmLf6aHU3UE/6aH23pgj//ZKnxTaqeeenvb0de/bsQV5enuu2gIAA5OXlYefOnX167qKiIgwbNgwpKSm48cYbUV5efsbHWCwWNDU1uf5o8cJuRNS9gb4bSI27NYm8STX/qVJfXw+Hw9HlI3XR0dE4ePCgx887btw4vPbaa8jMzERVVRUee+wxTJw4ET/++CPMZrPs40aOHImWlhbX3/Pz83HnnXciNjYWJSUlrrVJkuTaLZeeno6Kigq0trYiMDAQCQkJKCoqAgBERUUhICAAkiShoKAAqampqK6uRnNzM0wmE5KSklBYWAgAiIiIgNFoxLFjxwAAycnJqKurg9VqxaBBg5CWloYDBw4AAMLDwxEUFITKylNn4E1KSkJDQwOampqg1+uRmZmJAwcOQJIkDB48GGazGUePHgUADB8+HE1NTTh58iR0Oh1GjBiBwsJCOBwOhIWFYciQIThy5AgAID4+Hi0tLWhoOHUxw+zsbBw6dAh2ux1msxkREREoLS0FAAwbNgw2m811MbysrCwcPnwY7e3tCAkJQXR0tOvTA7GxsbDb7airq4MkSbDb7SgvL0dbWxuCgoIQFxeH4uJiV28AqKmpAQCkpaWhsrLS1Xv48OE4dOgQACAyMhIGgwFVVVUAgJSUFNTU1KC5uRlGoxEpKSmu7Wro0KEwmUxuvevr62GxWGAwGJCRkYGCggJX7+DgYFRUVAAAEhMTceLECdneYWFhrmE7ISEBFotFtnd4eDjKysoAAHFxcWhtbXX1HjFiBPR6PQoKChAaGorIyEi33u3t7aivrwcAZGZmoqysDDabDSEhIYiJiXFtszExMXA6nW7b7NGjR1294+Pj3bZZnU7n6p2amoqqqiq0tLTAZDIhMTFR2Lu2thZWq7Xb3oGBgd1us6f3HjJkCEJDQ9222cbGRjQ2NiIgIABZWVk4ePAgnE4nzGYzrFarW2+r1YoTJ0502Wa7693W1ua2zQ5uq8XinxnwfaMB2ys7UG3tQKfIYD0mRAPnRuoxOiYYJSUlsNlsCA4O7vNrRHV1tau3N18jOl97/Pk1AgAyMjL65TUiMTERR44c8evXiOLiYnR0dHj9NcJoNKKgoED1rxE9vFoXgF5c26u/HTt2DHFxcfj6668xfvx41+333XcfduzYgW+++Ub4+KSkJCxYsAALFiwQ3u/kyZNITEzEM888g5tvvrnL1zuv7ZWSkuJ2LROTyQSTydS7H6obhw8fRkpKSp+fR4vYRox9xPqzjyRJfr8biNuPGPvI85c2/XJtr/4WEREBvV7vmiA71dTUePV4ncGDByMjI8P1XwpyzGZzv1zYtK2tzevPqRVsI8Y+Yv3Zp3M3kD9fx4rbjxj7yNNiG9Uc82M0GjF27Fhs27bNdZvT6cS2bdvc3gnqK6vVipKSEsU+uhcUFKTI9/UHbCPGPmLsI8Y+YuwjT4ttVPWfMQsXLsTs2bORm5uLc889FytXrkRzc7Pr01+zZs1CXFwcli9fDuDUQdKd+/3a29tRWVmJ7777DqGhoa4r0N57772YOnUqEhMTcezYMSxZsgR6vR4zZ85U5GeMi4tT5Pv6A7YRYx8x9hFjHzH2kafFNqp55wcAZsyYgRUrVmDx4sUYM2YMvvvuO2zatMl1EFt5ebnrYCng1HFC55xzDs455xxUVVVhxYoVOOecc3DLLbe47lNRUYGZM2ciMzMT1113HYYOHYpdu3Yp9rG9M+1uG8jYRox9xNhHjH3E2EeeFtuo6p0fAJg/fz7mz5/f7de2b9/u9vekpKQzHt397rvvemtpREREpAGqeudnINDi1XG9hW3E2EeMfcTYR4x95GmxDYcfIiIiGlA4/PjY6R/lp/9hGzH2EWMfMfYRYx95WmzD4YeIiIgGFA4/Ptb5EXzqim3E2EeMfcTYR4x95GmxDYcfH+u8Vgl1xTZi7CPGPmLsI8Y+8rTYhsOPj7W2tiq9BNViGzH2EWMfMfYRYx95WmzD4cfHAgMDlV6CarGNGPuIsY8Y+4ixjzwttuHw42PDhw9XegmqxTZi7CPGPmLsI8Y+8rTYhsOPjx06dEjpJagW24h110eSJDS22VFtsaGxzX7GM55rGbcfMfYRYx95WmyjustbENGZWW12bClqwMb9daiytLtujzUbMS0nEpPTwxFq4r/eRETd4aujjyl1QVV/wDZinX12VzRh6dZS2OzOLveptrRj1a5KrN1dhcV5yciND/P1MhXD7UeMfcTYR54W23C3l48ZDJw35bCNmMFgwO6KJizaXAKb3QkJwOk7uTpvs9mdWLS5BLsrmny/UIVw+xFjHzH2kafFNhx+fKyqqkrpJagW24iVHD2GpVtLIUldh57TSQAkCVi6tRRWm90Xy1Mctx8x9hFjH3labMPhh8hP/LvW4XrHpyc63wHaUtTQn8siIvI7HH58LCUlReklqBbbyJMkCTvrA3o8+PzUxv11A+JTYNx+xNhHjH3kabENhx8f0+LVcb2FbeQ12Ryotnb0+nESgCpLOyw2h/cXpTLcfsTYR4x95GmxDYcfH2tublZ6CarFNvJaO/o2vLT08fH+gNuPGPuIsY88Lbbh8ONjRqNR6SWoFtvICxqk79Pjg/v4eH/A7UeMfcTYR54W23D48TEt7jv1FraRF2bSI9ZshK6Xj9Ph1IkPzSbtDz/cfsTYR4x95GmxDYcfHzt48KDSS1AttpGn0+nwi6Ge7bqalhMJna63Y5P/4fYj5g99lLxciz/0UYoW22jvzEVEGvXzKD3+XiH1+OPuOgAmQwAmp4f399KI+oSXayFf4zs/PjZ06FCll6BabCOWEB2BxXnJ0Olwxt1fOgA6HbA4L3nA/NLg9iOm1j67K5pwwzv7sWpXJap/MvgA/7tcyw3v7O/3s5WrtY8aaLENhx8fM5lMSi9BtdhGzGQyITc+DMumpMJkCDg14Jx2n87bTIYALJuSOqCu7cXtR0yNfdR0uRY19lELLbbh8ONjx44dU3oJqsU2Yp19cuPDsG5mDm77RRxizO6fwogxG3HbL+Lwzg0jB9TgA3D7ORO19bHa7Kq6XIva+qiJFtsMjPfDiTQm1GTA9JFRmJYTCYvNgZYOB4IH6WE26QfEwc3k/7YUNXh8uZbpI6P6c2k0APCdHx9LTk5WegmqxTZi3fXR6XQICzQgxmxCWKBhQA8+3H7E1NRHkqRTl13x4LH9dbkWNfVRGy224fDjY/X19UovQbXYRox9xNhHTE19mmwOt0919VR/Xq5FTX3URottOPz4mMViUXoJqsU2Yuwjxj5iauqjxsu1qKmP2mixDYcfHzMYeJiVHLYRYx8x9hFTUx81Xq5FTX3URottOPz4WEZGhtJLUC22EWMfMfYRU1MfNV6uRU191EaLbTj8+FhBQYHSS1AtthFjHzH2EVNTH51Oh2k5kR49tr8u16KmPmqjxTYcfoiIyOcmp4e7TtbZE7xcC3kThx8fCw/nv7hy2EaMfcTYR0xtfUJNBlVdrkVtfdREi204/PhYcHCw0ktQLbYRYx8x9hFTYx81Xa5FjX3UQottOPz4WEVFhdJLUC22EWMfMfYRU2sftVyuRa191ECLbbT3+TUiIvIrvFwL+RqHHx9LTExUegmqxTZi7CPGPmL+0Kfzci1hgb7/1eQPfZSixTbc7eVjJ06cUHoJqsU2Yuwjxj5i7CPGPvK02IbDj481NTUpvQTVYhsx9hFjHzH2EWMfeVpsw+HHx/R675+ZVCvYRox9xNhHjH3E2EeeFtvoJEmSlF6EmtjtduzYsQOjR4/W5D9wIiIiLXI4HNi3bx8mTZp0xuuR8Z0fHztw4IDSS1AtthFjHzH2EWMfMfaRp8U2HH58jG+0yWMbMfYRYx8x9hFjH3labMPhx8cGDx6s9BJUi23E2EeMfcTYR4x95GmxDYcfHwsL69+zlPozthFjHzH2EWMfMfaRp8U2HH58rLy8XOklqBbbiLGPGPuIsY8Y+8jTYhsOP0RERDSgcPjxsYSEBKWXoFpsI8Y+Yuwjxj5i7CNPi204/PiYxWJRegmqxTZi7CPGPmLsI8Y+8rTYhsOPj508eVLpJagW24ixjxj7iLGPGPvI02Ib1Q0/L774IpKSkhAYGIhx48bh22+/lb3v/v37cfXVVyMpKQk6nQ4rV67s83P2N51Op9j3Vju2EWMfMfYRYx8x9pGnxTaqGn7Wr1+PhQsXYsmSJdi7dy9Gjx6NKVOmoLa2ttv7t7S0ICUlBU888QRiYmK88pz9bcSIEYp8X3/ANmLsI8Y+Yuwjxj7ytNhGVcPPM888g7lz5yI/Px/Z2dlYtWoVgoODsWbNmm7v//Of/xx/+tOfcP3118NkMnnlOftbYWGhIt/XH7CNGPuIsY8Y+4ixjzwttlHN8NPe3o49e/YgLy/PdVtAQADy8vKwc+dOnz+nxWJBU1OT64/NZvNoDadzOBxeeR4tYhsx9hFjHzH2EWMfeVpsI77sqQ/V19fD4XAgOjra7fbo6GgcPHjQ5885cuRItLS0uP6en5+PO++8E7GxsSgpKXE9jyRJrl1o6enpqKioQGtrKwIDA5GQkICioiIAQFRUFAICAtDa2oqCggKkpqaiuroazc3NMJlMSEpKck3XERERMBqNOHbsGAAgOTkZdXV1sFqtGDRoENLS0lwXmgsPD0dQUBAqKysBAElJSWhoaEBTUxP0ej0yMzNx4MABSJKEwYMHw2w24+jRowCA4cOHo6mpCSdPnoROp8OIESNQWFgIh8OBsLAwDBkyBEeOHAEAxMfHo6WlBQ0NDQCA7OxsHDp0CHa7HWazGRERESgtLQUADBs2DDabDcePHwcAZGVl4fDhw2hvb0dISAiio6Nx+PBhAEBsbCzsdjvq6urQ2toKu92O8vJytLW1ISgoCHFxcSguLnb1BoCamhoAQFpaGiorK129hw8fjkOHDgEAIiMjYTAYUFVVBQBISUlBTU0NmpubYTQakZKS4toGhg4dCpPJ5Na7vr4eFosFBoMBGRkZKCgocPUODg5GRUUFACAxMREnTpyQ7R0WFuY6QVhCQgIsFots7/DwcJSVlQEA4uLi0Nra6uo9YsQItLe3o6CgAKGhoYiMjHTr3d7ejvr6egBAZmYmysrKYLPZEBISgpiYGNc2GxMTA6fT6bbNHj161NU7Pj7ebZvV6XSu3qmpqaiqqkJLSwtMJhMSExOFvWtra2G1WrvtHRgY2O02e3rvIUOGIDQ01G2bbWxsRGNjIwICApCVlYWDBw/C6XRCkiRYrVa33larFSdOnOiyzXbXu62trdttNjQ0FFFRUd1uswCQkZGBI0eOwGazITg4uM+vEdXV1a7e3nyN6Hzt8efXiM7e/fEaERwcjCNHjvj1a0RxcTE6Ojq8/hpht9tRUFCg+teI3lyDTCep5Iplx44dQ1xcHL7++muMHz/edft9992HHTt24JtvvhE+PikpCQsWLMCCBQv69Jx2ux07duxASkoKAgL+98aYyWSS3bXWGy0tLQgODu7z82gR24ixjxj7iLGPGPvI85c2DocD+/btw6RJk2AwiN/bUc1ur4iICOj1etcE2ammpkb2YOb+fE6z2YywsDDXH28MPgBcUzt1xTZi7CPGPmLsI8Y+8rTYRjXDj9FoxNixY7Ft2zbXbU6nE9u2bXN710bp5yQiIiL/pppjfgBg4cKFmD17NnJzc3Huuedi5cqVaG5uRn5+PgBg1qxZiIuLw/LlywHAdQxE5/+urKzEd999h9DQUKSlpfXoOX0tLi5Oke/rD9hGjH3E2EeMfcTYR54W26hq+JkxYwbq6uqwePFiVFdXY8yYMdi0aZPrILby8nK343COHTuGc845x/X3FStWYMWKFZg0aRK2b9/eo+f0tdbWVpx11lmKfG+1Yxsx9hFjHzH2EWMfeVpso6rhBwDmz5+P+fPnd/u1zoGmU1JSUo+O7hY9p681NDR4fAyT1rGNGPuIsY8Y+4ixjzwttlHNMT9EREREvqCaj7qrRedH3UePHg29Xu/155ckSZPXSfEGthFjHzH2EWMfMfaR5y9t/PKj7gNF58m4qCu2EWMfMfYRYx8x9pGnxTYcfnyso6ND6SWoFtuIsY8Y+4ixjxj7yNNiGw4/PhYaGqr0ElSLbcTYR4x9xNhHjH3kabENhx8fi4yMVHoJqsU2Yuwjxj5i7CPGPvK02IbDj491XmiOumIbMfYRYx8x9hFjH3labMPhh4iIiAYUDj8+NmzYMKWXoFpsI8Y+Yuwjxj5i7CNPi21Ud4ZnrWtvb1d6CT4hSRKabA60djgQNEiPMJP+jOeJGChtPMU+Yuwjxj5i7CNPi204/PhYfX09oqKilF5Gv7Ha7NhS1ICN++tQZfnfvzCxZiOm5URicno4Qk3db3Zab9NX7CPGPmLsI8Y+8rTYhsMPec3uiiYs3VoKm93Z5WvVlnas2lWJtbursDgvGbnxYQqskIiIiJe36KK/L2/hcDj65XmVtruiCYs2l0CSANEGpQOg0wHLpqR2GYC02sZb2EeMfcTYR4x95PlLG17eQsXKysqUXoLXWW12LN1aesbBBzj1dUkClm4thdVmd/uaFtt4E/uIsY8Y+4ixjzwttuHw42M2m03pJXjdlqIG2OzOMw4+nSQANrsTW4oa3G7XYhtvYh8x9hFjHzH2kafFNhx+fCwkJETpJXiVJEnYuL+ux4PPT23cX4ef7nXVWhtvYx8x9hFjHzH2kafFNhx+fCwmJkbpJXhVk83h9qmunpIAVFnaYbE5XLdprY23sY8Y+4ixj9hA6yNJEhrb7Ki22NDYZofo8F8ttuHw42MlJSVKL8GrWjscZ76TQMtPHq+1Nt7GPmLsI8Y+YgOlj9Vmx4c/1mLOewW49q0fMGv9qf8/570CfPhjbZdjMQFttuFH3alPggb17RMAwX18PBER9QxPR/I/fOfHx7T29mGYSY9YsxHiczd3pcOpEx+aTf8bfrTWxtvYR4x9xNhHTOt9Ok9H0vnhlNN3cnXeZrM7sWhzCXZXNLm+psU2HH58zOnsOnH7M51Oh2k5kR49dlpOpNslL7TWxtvYR4x9xNhHTMt9+no6Ei224fDjY7W1tUovwesmp4fDZAjo8bs/OgAmQwAmp4e73a7FNt7EPmLsI8Y+Ylru09fTkWixDYcf6rNQkwGL85Kh0+GMA1DnGZ4X5yXLXuOLiIi8w5unI9ESDj8+lp6ervQS+kVufBiWTUl1vQN0+hDUeZvJENDtpS0A7bbxFvYRYx8x9hHTah9vnI5Ei204/PjY0aNHlV5Cv8mND8O6mTm47RdxiDEb3b4WYzbitl/E4Z0bRsp+ikDLbbyBfcTYR4x9xLTaxxunI9FiG+538LG2tjall9CvQk0GTB8ZhWk5kbDYHGjpcCB4kB5mk97t4ObuaL1NX7GPGPuIsY+YVvt443QkDRpsw+HHx4KCgpRegk/odDqEBRoQFtjzTWygtPEU+4ixjxj7iGm1T+fpSKot7b067keHU+/Ym016Tbbhbi8fi4+PV3oJqsU2Yuwjxj5i7COm1T7eOB2JFttw+PGxoqIipZegWmwjxj5i7CPGPmJa7tPX05FosQ2HHyIiIg3j6Ui64vDjY1FRUUovQbXYRox9xNhHjH3EtN6nL6cj0WIb7Y51KnWmTzwNZGwjxj5i7CPGPmIDoU/n6Ui2FDVg4/46t/P/xJiNmJYTiYsyhiLE6P4JMS224fDjYzU1NRg6dKjSy1AlthFjHzH2EWMfsYHSx5PTkWixDYcfIiKiAcaT05FoCY/58bHU1FSll6BabCPGPmLsI8Y+YuwjT4ttOPz4WFVVldJLUC22EWMfMfYRYx8x9pGnxTYcfnzIZrPhhRdegM1mU3opqsM2Yuwjxj5i7CPGPvK02obDjw/ZbDasXbtWcxuRN7CNGPuIsY8Y+4ixjzyttuHwQ0RERAMKhx8iIiIaUAbmZ9wEJOnUdW8dDofXn9vpdCI4OBhOp7Nfnt+fsY0Y+4ixjxj7iLGPPH9q07m+zt/jIjqpJ/caQNra2vCvf/1L6WUQERGRByZMmIDAwEDhfTj8nMbpdKK9vR16vfzZLomIiEhdJEmCw+GA0WhEQID4qB4OP0RERDSg8IBnIiIiGlA4/BAREdGAwuHHR1588UUkJSUhMDAQ48aNw7fffqv0klTjn//8J6ZOnYphw4ZBp9Nh48aNSi9JNZYvX46f//znMJvNiIqKwrRp01BYWKj0slTjpZdewqhRoxAWFoawsDCMHz8ef//735Velio98cQT0Ol0WLBggdJLUYVHH30UOp3O7U9WVpbSy1KVyspK3HTTTRg6dCiCgoJw9tlnY/fu3Uovyys4/PjA+vXrsXDhQixZsgR79+7F6NGjMWXKFNTW1iq9NFVobm7G6NGj8eKLLyq9FNXZsWMH7rjjDuzatQtbtmxBR0cHLrroIjQ3Nyu9NFWIj4/HE088gT179mD37t248MILceWVV2L//v1KL01V/v3vf+Pll1/GqFGjlF6KquTk5KCqqsr156uvvlJ6Sapx4sQJTJgwAYMGDcLf//53FBQU4Omnn8aQIUOUXpp3SNTvzj33XOmOO+5w/d3hcEjDhg2Tli9fruCq1AmA9OGHHyq9DNWqra2VAEg7duxQeimqNWTIEOnVV19VehmqYbFYpPT0dGnLli3SpEmTpLvvvlvpJanCkiVLpNGjRyu9DNW6//77pfPOO0/pZfQbvvPTz9rb27Fnzx7k5eW5bgsICEBeXh527typ4MrIHzU2NgIAwsPDFV6J+jgcDrz77rtobm7G+PHjlV6Oatxxxx247LLL3F6D6JSioiIMGzYMKSkpuPHGG1FeXq70klTj448/Rm5uLq699lpERUXhnHPOwerVq5Veltdw+Oln9fX1cDgciI6Odrs9Ojoa1dXVCq2K/JHT6cSCBQswYcIEjBw5UunlqMYPP/yA0NBQmEwm3Hbbbfjwww+RnZ2t9LJU4d1338XevXuxfPlypZeiOuPGjcNrr72GTZs24aWXXkJpaSkmTpwIi8Wi9NJU4fDhw3jppZeQnp6OzZs3Y968ebjrrrvw+uuvK700r+DlLYj8xB133IEff/yRxyWcJjMzE9999x0aGxvx/vvvY/bs2dixY8eAH4COHj2Ku+++G1u2bDnj2W4HoksuucT1v0eNGoVx48YhMTER7733Hm6++WYFV6YOTqcTubm5ePzxxwEA55xzDn788UesWrUKs2fPVnh1fcd3fvpZREQE9Ho9ampq3G6vqalBTEyMQqsifzN//nz87W9/wxdffIH4+Hill6MqRqMRaWlpGDt2LJYvX47Ro0fjueeeU3pZituzZw9qa2vxs5/9DAaDAQaDATt27MDzzz8Pg8Gg+us0+drgwYORkZGB4uJipZeiCrGxsV3+A2LEiBGa2TXI4aefGY1GjB07Ftu2bXPd5nQ6sW3bNh6XQGckSRLmz5+PDz/8EP/4xz+QnJys9JJUz+l0wmazKb0Mxf3qV7/CDz/8gO+++871Jzc3FzfeeCO+++476PV6pZeoKlarFSUlJYiNjVV6KaowYcKELqfVOHToEBITExVakXdxt5cPLFy4ELNnz0Zubi7OPfdcrFy5Es3NzcjPz1d6aapgtVrd/murtLQU3333HcLDwzF8+HAFV6a8O+64A+vWrcNHH30Es9nsOk7srLPOQlBQkMKrU96DDz6ISy65BMOHD4fFYsG6deuwfft2bN68WemlKc5sNnc5NiwkJARDhw7lMWMA7r33XkydOhWJiYk4duwYlixZAr1ej5kzZyq9NFW455578H//9394/PHHcd111+Hbb7/FK6+8gldeeUXppXmH0h83GyheeOEFafjw4ZLRaJTOPfdcadeuXUovSTW++OILCUCXP7Nnz1Z6aYrrrgsAae3atUovTRV+85vfSImJiZLRaJQiIyOlX/3qV9Lnn3+u9LJUix91/58ZM2ZIsbGxktFolOLi4qQZM2ZIxcXFSi9LVT755BNp5MiRkslkkrKysqRXXnlF6SV5DS9sSkRERAMKj/khIiKiAYXDDxEREQ0oHH6IiIhoQOHwQ0RERAMKhx8iIiIaUDj8EBER0YDC4YeIiIgGFA4/RERENKBw+CEizXn00Ueh0+mUXobHysrKoNPp8Nprr7lu8/efiUhNOPwQ+bHXXnsNOp3O9ScwMBAZGRmYP38+ampq3O5bVlaG/Px8pKamIjAwEDExMTj//POxZMmSbp/7k08+wdSpUxEdHQ2j0Yjw8HCcf/75ePrpp9HU1OSLH4+IqF/wwqZEGrB06VIkJyejra0NX331FV566SV89tln+PHHHxEcHIzi4mL8/Oc/R1BQEH7zm98gKSkJVVVV2Lt3L5588kk89thjrudyOp24+eab8dprr+Hss8/G7bffjoSEBFgsFuzcuROLFi3CZ599hm3btin4Ew88ixYtwgMPPKD0Mog0gcMPkQZccsklyM3NBQDccsstGDp0KJ555hl89NFHmDlzJp599llYrVZ89913SExMdHtsbW2t29+feuopvPbaa7jnnnvw9NNPu+1qufvuu1FVVYU33nij/38oP2S32+F0OmE0Gr3+3AaDAQYDX7KJvIG7vYg06MILLwQAlJaWAgBKSkoQHx/fZfABgKioKNf/bmlpwZNPPomcnBz86U9/6vYYk9jYWNx///1nXMOXX36Ja6+9FsOHD4fJZEJCQgLuuecetLa2ut1vzpw5CA0NRWVlJaZNm4bQ0FBERkbi3nvvhcPhcLvv8ePH8etf/xphYWEYPHgwZs+ejX379nU5PkbOW2+9hbFjxyIoKAjh4eG4/vrrcfTo0TM+rjudx+WsWLECK1euRGpqKkwmEwoKCtDe3o7Fixdj7NixOOussxASEoKJEyfiiy++6PI8J0+exJw5c3DWWWe5fqaTJ092uV93x/ysXbsWF154IaKiomAymZCdnY2XXnqpy2OTkpJw+eWX46uvvsK5556LwMBApKSkdBliOzo68NhjjyE9PR2BgYEYOnQozjvvPGzZssWjRkRqxf+MINKgkpISAMDQoUMBAImJidi6dSv+8Y9/uAaj7nz11Vc4efIk7r33Xuj1+j6t4a9//StaWlowb948DB06FN9++y1eeOEFVFRU4K9//avbfR0OB6ZMmYJx48ZhxYoV2Lp1K55++mmkpqZi3rx5AE7tjps6dSq+/fZbzJs3D1lZWfjoo48we/bsHq3nj3/8Ix555BFcd911uOWWW1BXV4cXXngB559/Pv7zn/9g8ODBHv2ca9euRVtbG2699VaYTCaEh4ejqakJr776KmbOnIm5c+fCYrHg//2//4cpU6bg22+/xZgxYwAAkiThyiuvxFdffYXbbrsNI0aMwIcfftjjn+mll15CTk4OrrjiChgMBnzyySe4/fbb4XQ6cccdd7jdt7i4GNdccw1uvvlmzJ49G2vWrMGcOXMwduxY5OTkADg1YC1fvhy33HILzj33XDQ1NWH37t3Yu3cvJk+e7FEfIlWSiMhvrV27VgIgbd26Vaqrq5OOHj0qvfvuu9LQoUOloKAgqaKiQpIkSfrxxx+loKAgCYA0ZswY6e6775Y2btwoNTc3uz3fc889JwGQNm7c6Ha73W6X6urq3P44nU7h2lpaWrrctnz5ckmn00lHjhxx3TZ79mwJgLR06VK3+55zzjnS2LFjXX//4IMPJADSypUrXbc5HA7pwgsvlABIa9eudd2+ZMkS6acvb2VlZZJer5f++Mc/un2PH374QTIYDF1u74nS0lIJgBQWFibV1ta6fc1ut0s2m83tthMnTkjR0dHSb37zG9dtGzdulABITz31lNtjJ06ceMafSZK6bzxlyhQpJSXF7bbExEQJgPTPf/7TdVttba1kMpmk3/3ud67bRo8eLV122WU9+OmJ/Bt3exFpQF5eHiIjI5GQkIDrr78eoaGh+PDDDxEXFwcAyMnJwXfffYebbroJZWVleO655zBt2jRER0dj9erVrufp/BRXaGio2/P/8MMPiIyMdPtz/Phx4ZqCgoJc/7u5uRn19fX4v//7P0iShP/85z9d7n/bbbe5/X3ixIk4fPiw6++bNm3CoEGDMHfuXNdtAQEBXd7h6M6GDRvgdDpx3XXXob6+3vUnJiYG6enp3e6O6qmrr74akZGRbrfp9XrXcT9OpxMNDQ2w2+3Izc3F3r17Xff77LPPYDAYXO9udT72zjvv7NH3/mnjxsZG1NfXY9KkSTh8+DAaGxvd7pudnY2JEye6/h4ZGYnMzEy3xoMHD8b+/ftRVFTUo+9P5K+424tIA1588UVkZGTAYDAgOjoamZmZCAhw/2+bjIwMvPnmm3A4HCgoKMDf/vY3PPXUU7j11luRnJyMvLw8mM1mAIDVanV7bFpamuu4jzfeeANvvvnmGddUXl6OxYsX4+OPP8aJEyfcvnb6L+bAwMAuA8SQIUPcHnfkyBHExsYiODi4y9rOpKioCJIkIT09vduvDxo06IzPISc5Obnb219//XU8/fTTOHjwIDo6Orq9f+fPdPqwmZmZ2aPv/a9//QtLlizBzp070dLS4va1xsZGnHXWWa6/Dx8+vMvjT2+8dOlSXHnllcjIyMDIkSNx8cUX49e//jVGjRrVo/UQ+QsOP0QacO6557o+7XUmer0eZ599Ns4++2yMHz8ev/zlL/H2228jLy8PWVlZAIAff/wRV155pesxoaGhyMvLA3DquKAzcTgcmDx5MhoaGnD//fcjKysLISEhqKysxJw5c+B0OrusqT85nU7odDr8/e9/7/Z7nT589MZP333p9NZbb2HOnDmYNm0afv/73yMqKgp6vR7Lly93HY/VVyUlJfjVr36FrKwsPPPMM0hISIDRaMRnn32GZ599tseNJUly/e/zzz8fJSUl+Oijj/D555/j1VdfxbPPPotVq1bhlltu8cq6idSAww/RANY5MFVVVQE4tavprLPOwrvvvosHH3ywy7tHPfXDDz/g0KFDeP311zFr1izX7X351FBiYiK++OILtLS0uL37U1xcfMbHpqamQpIkJCcnIyMjw+M19NT777+PlJQUbNiwwe0TWqefUDIxMRHbtm2D1Wp1G8AKCwvP+D0++eQT2Gw2fPzxx27v6vRlFx4AhIeHIz8/H/n5+bBarTj//PPx6KOPcvghTeExP0QDwJdffum266XTZ599BuB/u1mCg4Nx33334ccff8QDDzzg9q5Ap+5uO13nuww/va8kSXjuuec8Wj8ATJkyBR0dHW7HKDmdTrz44otnfOxVV10FvV6Pxx57rMv6JUk64/FLvdXdz//NN99g586dbve79NJLYbfb3T6e7nA48MILL3j0PRobG7F27VqP1316h9DQUKSlpcFms3n8nERqxHd+iAaAJ598Env27MFVV13lOn5j7969eOONNxAeHo4FCxa47vvAAw/gwIED+NOf/oTPP/8cV199NeLj43HixAns3bsXf/3rXxEVFYXAwEDZ75eVlYXU1FTce++9qKysRFhYGD744IMux/70xrRp03Duuefid7/7HYqLi5GVlYWPP/4YDQ0NACC87lVqaiqWLVuGBx98EGVlZZg2bRrMZjNKS0vx4Ycf4tZbb8W9994LANi+fTt++ctfYsmSJXj00Uc9Wuvll1+ODRs2YPr06bjssstQWlqKVatWITs72+14qqlTp2LChAl44IEHUFZWhuzsbGzYsKHLMVHdueiii2A0GjF16lT89re/hdVqxerVqxEVFeV6J6+3srOzccEFF2Ds2LEIDw/H7t278f7772P+/PkePR+RWnH4IRoAHnroIaxbtw47duzA22+/jZaWFsTGxuL666/HI4884nYQbkBAAN58801cffXVWL16NV544QWcOHECoaGhGDlyJP74xz9i7ty5wuNkBg0ahE8++QR33XUXli9fjsDAQEyfPh3z58/H6NGjPfoZ9Ho9Pv30U9x99914/fXXERAQgOnTp2PJkiWYMGGCcBgDTg11GRkZePbZZ12X80hISMBFF12EK664wnW/zuEkNjbWo3UCp07cWF1djZdffhmbN29GdnY23nrrLfz1r3/F9u3bXfcLCAjAxx9/jAULFuCtt96CTqfDFVdcgaeffhrnnHOO8HtkZmbi/fffx6JFi3DvvfciJiYG8+bNQ2RkJH7zm994tO677roLH3/8MT7//HPYbDYkJiZi2bJl+P3vf+/R8xGplU7qyXvYREQqtXHjRkyfPh1fffUVJkyY0Ofnu++++/DOO++guLgYJpPJCyskIrXhMT9E5DdOvzRG5/ExYWFh+NnPfuaV7/HFF1/gkUce4eBDpGHc7UVEfuPOO+9Ea2srxo8fD5vNhg0bNuDrr7/G448/3u1Hzj3x73//2yvPQ0Tqxd1eROQ31q1bh6effhrFxcVoa2tDWloa5s2bxwNyiahXOPwQERHRgMJjfoiIiGhA4fBDREREAwqHHyIiIhpQOPwQERHRgMLhh4iIiAYUDj9EREQ0oHD4ISIiogGFww8RERENKBx+iIiIaED5/+NDs6TbfKOOAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('bmh')\n",
    "\n",
    "plt.figure()\n",
    "plt.ylabel('Power observed')\n",
    "plt.xlabel('PSG angle, radians')\n",
    "plt.plot(thetas, power_measured, marker='o', linestyle='None', markersize=10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In katsu.polarimetry we have the data reduction available to convert this power measurement into a Mueller matrix measured. The full_mueller_polarimetry function takes dictionaries as inputs which store the state of the instrument. Below we show the defaults, but this freedom enables the flexibility necessary to use this technique accurately in the laboratory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "from katsu.polarimetry import full_mueller_polarimetry\n",
    "\n",
    "# The angles all optics begin at\n",
    "starting_angles={'psg_polarizer': 0,\n",
    "                 'psg_waveplate': 0,\n",
    "                 'psa_waveplate': 0,\n",
    "                 'psa_polarizer': 0}\n",
    "\n",
    "starting_polarization={'psg_Tmin': 0,\n",
    "                       'psg_ret': np.pi / 2,\n",
    "                       'psa_Tmin': 0,\n",
    "                       'psa_ret': np.pi / 2}\n",
    "\n",
    "M_measured = full_mueller_polarimetry(thetas, power_measured, PSA_ANGULAR_INCREMENT,\n",
    "                                      starting_angles=starting_angles,\n",
    "                                      starting_polarization=starting_polarization)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mueller Matrix Original\n",
      "[[0.5        0.29399553 0.40443371 0.        ]\n",
      " [0.28578716 0.16804029 0.23116392 0.        ]\n",
      " [0.40673081 0.23915408 0.3289913  0.        ]\n",
      " [0.05381217 0.03164107 0.04352691 0.        ]]\n",
      "------------------------------\n",
      "Mueller Matrix Measured\n",
      "[[ 5.00000000e-01  2.93995530e-01  4.04433713e-01 -5.96744876e-16]\n",
      " [ 2.85787156e-01  1.68040293e-01  2.31163922e-01  6.10622664e-16]\n",
      " [ 1.02871029e-01  6.04872455e-02  8.32090246e-02 -5.34294831e-16]\n",
      " [ 2.05742058e-01  1.20974491e-01  1.66418049e-01  2.49800181e-16]]\n",
      "------------------------------\n",
      "Percent Difference\n",
      "[[-2.22044605e-13 -4.72040776e-13  3.84318167e-13             inf]\n",
      " [ 6.21566365e-13  1.04058451e-12 -2.40137608e-13            -inf]\n",
      " [ 7.47078346e+01  7.47078346e+01  7.47078346e+01             inf]\n",
      " [-2.82333711e+02 -2.82333711e+02 -2.82333711e+02            -inf]]\n",
      "------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Work\\AppData\\Local\\Temp\\ipykernel_17996\\1858271723.py:8: RuntimeWarning: divide by zero encountered in divide\n",
      "  print(100*(M_to_measure-M_measured)/M_to_measure)\n"
     ]
    }
   ],
   "source": [
    "print('Mueller Matrix Original')\n",
    "print(M_to_measure)\n",
    "print('-'*30)\n",
    "print('Mueller Matrix Measured')\n",
    "print(M_measured)\n",
    "print('-'*30)\n",
    "print('Percent Difference')\n",
    "print(100*(M_to_measure-M_measured)/M_to_measure)\n",
    "print('-'*30)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pardon the rightmost column blowing up, but those were all zeros anyway. We successfully recover some of the information for the top two rows but not the bottom, can this be fixed with more measurements?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAGwCAYAAABGogSnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB4LUlEQVR4nO3de3yT9d3/8XeaNEkPCfR8pucDbQUcTIacnCtUp47ipoBToM4THiYydR4QPOBQJ54YP3F4A4qKKGLVzYHAhHkAHaCIbSn0RGnpkUKT9JA0yfX7o0skba9vmzTJdeXK5/l4cO9umqTfvLxy9dsr10HGcRwHQgghhBA/ESD0AAghhBBCvIkmP4QQQgjxKzT5IYQQQohfockPIYQQQvwKTX4IIYQQ4ldo8kMIIYQQv0KTH0IIIYT4FZr8EEIIIcSvKIQegNhYrVaYTCbI5XLIZDKhh0MIIYSQYeA4DhaLBUqlEgEB7G07NPnpx2Qy4auvvhJ6GIQQQghxwdSpU6FWq5n3oclPP3K5HACQn59v///d6fjx48jJyXH780oBtWGjPmzUh436sFEffr7SxmKx4McffxzW726a/PRj+6hLLpd7ZPITGxvrkeeVAmrDRn3YqA8b9WGjPvx8rc1wdlmhHZ4JIYQQ4ldo8uNlzc3NQg9BtKgNG/Vhoz5s1IeN+vCTYhua/BBCCCHEr9Dkx8syMjKEHoJoURs26sNGfdioDxv14SfFNjT58bKGhgahhyBa1IaN+rBRHzbqw0Z9+EmxDU1+vKy7u1voIYgWtWGjPmzUh436sFEfflJsQ4e6e9lQJ17yZ9SGjfqwuasPx3HQGS3o7rUgKFAOrUoaZ3un5YeN+vCTYhua/HjZmDFjhB6CaFEbNurDNtI+BqMZu0+2o6S0FY16k/32OI0SRXlRmJUZjlCV764yaflhoz78pNiGPvbyshMnTgg9BNGiNmzUh20kfQ7V63DD1lKsP9iApgsmPgDQpDdh/cEG3LC1FIfqdSMdpmBo+WGjPvyk2IYmP4QQv3aoXoflu6pgNFvBAeD6fd92m9FsxfJdVT49ASKE9KHJj5dFRUUJPQTRojZs1IfNlT4GoxlP7qkBxw2c9PTHAeA44Mk9NTAYzS6NUUi0/LBRH35SbEOTHy9TKHx3nwFPozZs1IfNlT67T7bbt/gMh20L0O6T7U7/LKHR8sNGffhJsQ1NfryssbFR6CGIFrVhoz5szvbhOA4lpa3DnvhcqKS0FRznyiOFQ8sPG/XhJ8U2NPkhhPglndHicFTXcHEAGvUm6I0W9w+KEOIVNPnxsrS0NKGHIFrUho36sDnbp7t3ZJOXrhE+3hUcx6Gjx4wmvREdPWantj7R8sNGffhJsY30PsgTuebmZiQnJws9DFGiNmzUh83ZPkGB8hH9vOARPt4Z7jgHES0/bNSHnxTb0JYfL+vs7BR6CKJFbdioD5uzfbQqOeI0Sjh77mYZ+iYdGpV3Jj/uOgcRLT9s1IefFNvQ5MfLlEql0EMQLWrDRn3YnO0jk8lQlOfaIbxFeVFeueSFO89BRMsPG/XhJ8U2NPnxMil+duou1IaN+rC50mdWZjhUioBhb/2RAVApAjArM9zpn+Usd5+DiJYfNurDT4ptaPLjZcePHxd6CKJFbdioD5srfUJVCqwoSIVMhiEnQDIAMhmwoiDVK9f4cvc5iGj5YaM+/KTYhiY/hBC/NilRi1WF6fYtQP0nQbbbVIoArCpMx6RErcfH5G/nICLE2+hoLy+LiIgQegiiRW3YqA/bSPpMStTinQV5gx5RFfu/I6pmZ0UgROmdnZzdcQ4irdpx9U7LDxv14SfFNjT58TKVSiX0EESL2rBRH7aR9glVKTA3PxpFeVHQGy3o6rUgOFAOjUrulZ2bL+SOcxD1n/zQ8sNGffhJsQ197OVlZ86cEXoIokVt2KgPm7v6yGQyaNUKxGpU0KoVXp/4AJ45BxEtP2zUh58U29DkhxBCRMZXzkFEiK+iyY+XpaamCj0E0aI2bNSHTUp9PHEOIin18QTqw0+KbWjy42VtbW1CD0G0qA0b9WGTWh93n4NIan3cjfrwk2Ibmvx4mV6vF3oIokVt2KgPm9T6uPscRFLr427Uh58U29Dkx8sUCjrAjg+1YaM+bFLs485zEEmxjztRH35SbCPj6GxYDsxmM/bv34/x48dDLqedBgkhwhvqqu7ePAcRIWJlsVhw9OhRzJw5c8gJm/SmcyJXVlaG3NxcoYchStSGjfqwSbmPO85BJOU+7kB9+EmxDU1+CCHER9jOQdT/BIaEEOfQPj9eFh7u+atB+ypqw0Z92KgPG/Vhoz78pNiGJj9eFhwcLPQQRIvasFEfNurDRn3YqA8/KbahyY+X1dfXCz0E0aI2bNSHjfqwUR826sNPim1o8kMIIYQQv0KTHy9LTk4WegiiRW3YqA8b9WGjPmzUh58U29Dkx8vOnTsn9BBEi9qwUR826sNGfdioDz8ptqHJj5fpdDqhhyBa1IaN+rBRHzbqw0Z9+EmxDU1+vIzOGs2P2rBRHzbqw0Z92KgPPym2octb9EOXtyCEEEJ8jzOXt6AtP15WXl4u9BBEi9qwUR826sNGfdioDz8ptqHJj5fRhjZ+1IaN+rBRHzbqw0Z9+EmxDV0gxstGjx4t9BA8juM46IwWdPdaEBQoh3aYF1/0hzYjQX3YqA8b9WGjPvyk2EaUW37WrVuHlJQUqNVqTJ48Gd9++y3vfXfs2IFJkyZh9OjRCAkJwYQJE7BlyxaH+yxevBgymczh3xVXXOHplzEorVYryM/1BoPRjA9/bMHi98pw3VvHsHBb3/8ufq8MH/7YAoPRzHy8lNu4A/Vhoz5s1IeN+vCTYhvRTX62bduGZcuWYeXKlThy5AjGjx+PwsJCtLS0DHr/8PBwPProozhw4AB++OEHFBcXo7i4GLt27XK43xVXXIHGxkb7v61bt3rj5QxQV1cnyM/1tEP1OtywtRTrDzagSW9y+F6T3oT1Bxtww9ZSHKrnP2RSqm3chfqwUR826sNGffhJsY3oJj8vvPACbr31VhQXFyM3Nxfr169HcHAwNm7cOOj9L7vsMsydOxdjx45Feno67r33XowbNw5ffvmlw/1UKhViY2Pt/8LCwrzxcvzCoXodlu+qgtFsBQeg/6fDttuMZiuW76piToAIIYQQTxPV5MdkMuHw4cMoKCiw3xYQEICCggIcOHBgyMdzHIe9e/eioqICM2bMcPjevn37EB0djezsbCxZsgRnz551+/iHIykpSZCf6ykGoxlP7qkBxw2c9PTHAeA44Mk9NYN+BCa1Nu5GfdioDxv1YaM+/KTYRlSTn7a2NlgsFsTExDjcHhMTg6amJt7HdXR0IDQ0FEqlEldddRXWrl2LWbNm2b9/xRVX4M0338TevXvx7LPPYv/+/bjyyithsVh4n1Ov10On09n/GY3Gkb/A/z2vlOw+2W7f4jMcti1Au0+2D/ie1Nq4G/Vhoz5s1IeN+vCTYhtJHO2l0Wjw/fffw2AwYO/evVi2bBnS0tJw2WWXAQDmz59vv+9FF12EcePGIT09Hfv27cOvfvWrQZ8zPz8fXV1d9q+Li4txzz33IC4uDlVVVQD6JmUcx9n3R8rMzER9fT26u7uhVquRlJSEkydPAgCio6MREBCA06dP4/z580hPT0dTUxM6OzuhUqmQkpKCiooKAEBkZCSUSiXOnDkDAEhNTUVraysMBgMCAwORkZFhP+9CeHg4goKC0NDQAABISUlBe3s7dDod5HI5srOzUV5eDo7jMHr0aGg0Gpw+fRoAMGbMGOh0Opw/fx4ymQxjx45FRUUFLBYLtFotwsLCcOrUKQBAYmIiurq60N7eN2nJzc1FRUUF3vuua9gTnwu9910DsmStiI+Ph9lsRmtrK/R6PaKjo1FXV4eenh4EBQUhISEBlZWV9t4A0NzcDADIyMhAQ0ODvfeYMWNw4sQJAEBUVBQUCgUaGxsBAGlpaWhubkZnZyeUSiXS0tJw/PhxAEBERARUKpVD77a2Nuj1eigUCmRlZaGsrMzeOzg4GPX19QD6Lvh37tw53t5ardb+eXlSUhL0ej1v7/DwcNTW1gIAEhIS0N3dbe89duxYNDQ04Pz58wgNDUVUVBRqamoAAPHx8TCZTGhrawMAZGdno7a2FkajESEhIYiNjbUvs7GxsbBarQ7L7OnTp+29ExMTHZZZmUxm752eno7GxkZ0dXVBpVIhOTmZ2bulpQUGg2HQ3mq1etBltn/vsLAwhIaGOiyzHR0d6OjoQEBAAHJycnD8+HFYrVb09PQM6G0wGOzXJcrNzcWJEydgNpsH7d3T02PfIpyTk4Pq6mqYTCaEhoYiOjoa1dXVAIC4uDj7MgsAWVlZOHXqFIxGI4KDg0e8jrD9oefudYRt3ePNdYStt0ajQWRkpMMyazQaB+0dEhKCmJgYZm9PrCMsFgt6e3t9eh1RWVmJ3t5et68jbOsesa8jnDkkX1RneDaZTAgODsb27dtRVFRkv33RokU4f/48Pvroo2E9zy233ILTp08P2On5QlFRUVi1ahVuv/12h9ttZ3hOS0tDQMBPG8ZUKhVUKpVzL2gQ5eXlGDt27IifRww6esy47q1jLj9++40XQav+af4tpTaeQH3YqA8b9WGjPvx8pY3PnuFZqVRi4sSJ2Lt3r/02q9WKvXv3YsqUKcN+HqvVyvyYqr6+HmfPnkVcXBzvfTQaDbRarf2fOyY+AHxiARqu7l7+jw2Ho6vf46XUxhOoDxv1YaM+bNSHnxTbiGryAwDLli3Dhg0b8MYbb6C8vBxLlixBZ2cniouLAQALFy7Eww8/bL//6tWrsXv3blRXV6O8vBxr1qzBli1bcOONNwIADAYDHnjgARw8eBC1tbXYu3cv5syZg4yMDBQWFnr99dk2W0tBUODIrn0W3O/xUmrjCdSHjfqwUR826sNPim1Et8/PvHnz0NraihUrVqCpqQkTJkzAzp077Z/l1tXVOXwc1dnZiTvvvBP19fUICgpCTk4O3nrrLcybNw9A39Vof/jhB7zxxhs4f/484uPjMXv2bDz11FNu25rjDNZO1r5Gq5IjTqNEk97k1H4/MgCxGiU0KsfJj5TaeAL1YaM+bNSHjfrwk2IbUe3zIwaevqp7fX09EhMT3f68QvnwxxasP9jg9OTnjl8kYG5+tMPtUmvjbtSHjfqwUR826sPPV9r47D4//iA8PFzoIbjVrMxwqBQBGPrKXX1kAFSKAMzKHNhBam3cjfqwUR826sNGffhJsQ1NfrzMdpiiVISqFFhRkAqZDENOgGQAZDJgRUEqQlUDZ+VSa+Nu1IeN+rBRHzbqw0+KbWjyQ0ZsUqIWqwrT7VuA+k+CbLepFAFYVZiOSYnSu0geIYQQ3yG6HZ6lLiEhQegheMSkRC3eWZCH3SfbUVLaisYLLm4aq1GiKC8Ks7MiEKLk349Kqm3chfqwUR826sNGffhJsQ1Nfrysu7sbo0aNEnoYHhGqUmBufjSK8qKgN1rQ1WtBcKAcGpUcMtnQewVJuY07UB826sNGfdioDz8ptqGPvbzMdipyKZPJZNCqFYjVqKBVK4Y18QH8o81IUB826sNGfdioDz8ptqHJDyGEEEL8Cp3npx9Pn+eH47hhbwnxN9SGjfqwUR826sNGffj5Shs6z4+I2a4+TAaiNmzUh436sFEfNurDT4ptaPLjZb29vUIPQbSoDRv1YaM+bNSHjfrwk2Ibmvx4WWhoqNBDEC1qw0Z92KgPG/Vhoz78pNiGJj9eFhUVJfQQRIvasFEfNurDRn3YqA8/KbahyY+X1dTUCD0E0aI2bNSHjfqwUR826sNPim1o8kMIIYQQv0KTHy+Lj48XegiiRW3YqA8b9WGjPmzUh58U29Dkx8tMJtPQd/JT1IaN+rBRHzbqw0Z9+EmxDU1+vKytrU3oIYgWtWGjPmzUh436sFEfflJsQ5MfQgghhPgVmvx4WXZ2ttBDEC1qw0Z92KgPG/Vhoz78pNiGJj9eVltbK/QQRIvasFEfNurDRn3YqA8/KbahyY+XGY1GoYcgWtSGjfqwUR826sNGffhJsQ1NfrwsJCRE6CGIFrVhoz5s1IeN+rBRH35SbEOTHy+LjY0VegiiRW3YqA8b9WGjPmzUh58U29Dkx8uqqqqEHoJoURs26sNGfdioDxv14SfFNjT5IYQQQohfocmPl0lx86G7UBs26sNGfdioDxv14SfFNjT58TKr1Sr0EESL2rBRHzbqw0Z92KgPPym2ocmPl7W0tAg9BNGiNmzUh436sFEfNurDT4ptaPLjJRzHoaPHjLM9ff/LcZzQQyKEEEK8yva7sElvFPR3oUKQn+pHDEYzdp9sR0lpKxr1/7sy7pFjiNMoUZQXhVmZ4QhV0X8GAMjMzBR6CKJGfdioDxv1YaM+/NzRZtDfhYBgvwtpy48HHarX4YatpVh/sAFNF/zHBoAmvQnrDzbghq2lOFSvE2iE4nL69GmhhyBq1IeN+rBRHzbqw2+kbcT4u5AmPx5yqF6H5buqYDRbwQHov2HPdpvRbMXyXVU0AQLQ09Mj9BBEjfqwUR826sNGffiNpI1YfxfS5McDDEYzntxTA44b+B+6Pw4AxwFP7qmBwWj2xvBEKygoSOghiBr1YaM+bNSHjfrwc7WNmH8X0uTHA3afbLfPcofDNuvdfbLdk8MSvcTERKGHIGrUh436sFEfNurDz9U2Yv5dSJMfN+M4DiWlrcP+j32hktJWvz4K7OTJk0IPQdSoDxv1YaM+bNSHnyttxP67kCY/bqYzWhz2ZB8uDkCj3gS90eL+QRFCCCFeJPbfhTT5cbPu3pH9B+sa4eN9WXR0tNBDEDXqw0Z92KgPG/Xh50obsf8upMmPmwUFykf0+OARPt6XyWQyoYcgatSHjfqwUR826sPPlTZi/11Ikx8306rkiNMo4eyiIkPfyZ40Kv+d/DQ3Nws9BFGjPmzUh436sFEffq60EfvvQpr8uJlMJkNRXpRLjy3Ki6K/PgghhPg8sf8upMmPB8zKDIdKETDsGa8MgEoRgFmZ4Z4cluilp6cLPQRRoz5s1IeN+rBRH36uthHz70Ka/HhAqEqBFQWpkMkw5H90GQCZDFhRkOr31/hqbGwUegiiRn3YqA8b9WGjPvxcbSPm34U0+fGQSYlarCpMt896+/+Ht92mUgRgVWE6JiVqvT9Ikenq6hJ6CKJGfdioDxv1YaM+/EbSRqy/C/17U4OHTUrU4p0FeYNeyTb2f1eynZ0VgRCl/+7kfCGVSiX0EESN+rBRHzbqw0Z9+I20jRh/F8o4fz6l8CDMZjP279+P8ePHQy53338IjuOgN1qg7zZBE9S3Jzvt3OzIbDZDoaD5OB/qw0Z92KgPG/Xh5842tt+FXb0WBAfK3fq70GKx4OjRo5g5c+aQ4xXlx17r1q1DSkoK1Go1Jk+ejG+//Zb3vjt27MCkSZMwevRohISEYMKECdiyZYvDfTiOw4oVKxAXF4egoCAUFBR4/VTmMpkMWrUCHY210KoVNPEZxIkTJ4QegqhRHzbqw0Z92KgPP3e2sf0ujNWoBP1dKLrJz7Zt27Bs2TKsXLkSR44cwfjx41FYWIiWlpZB7x8eHo5HH30UBw4cwA8//IDi4mIUFxdj165d9vs899xzeOWVV7B+/Xp88803CAkJQWFhIXp6erz1sgghhBAiEqKb/Lzwwgu49dZbUVxcjNzcXKxfvx7BwcHYuHHjoPe/7LLLMHfuXIwdOxbp6em49957MW7cOHz55ZcA+rb6vPTSS1i+fDnmzJmDcePG4c0338SZM2dQUlLixVfWJyrKtfMe+ANqw0Z92KgPG/Vhoz78pNhGVJMfk8mEw4cPo6CgwH5bQEAACgoKcODAgSEfz3Ec9u7di4qKCsyYMQMAUFNTg6amJofnHDVqFCZPnsx8Tr1eD51OZ/9nNBpH8Mp+Qp8p86M2bNSHjfqwUR826sNPim1E9Yra2tpgsVgQExPjcHtMTAyOHz/O+7iOjg4kJCTAaDRCLpfj//2//4dZs2YBAJqamuzP0f85bd8bTH5+vsPhfcXFxbjnnnsQFxeHqqoq+3NwHGf/SC4zMxP19fXo7u6GWq1GUlKSfd+i6OhoBAQE4OTJk9BoNEhPT0dTUxM6OzuhUqmQkpKCiooKAEBkZCSUSiXOnDkDAEhNTUVraysMBgMCAwORkZGB8vJyAH0f+wUFBaGhoQEAkJKSgvb2duh0OsjlcmRnZ6O8vBwcx2H06NHQaDQ4ffo0AGDMmDHQ6XQ4f/48ZDIZxo4di4qKClgsFmi1WoSFheHUqVMAgMTERHR1daG9vR0AkJubixMnTsBsNkOj0SAyMhI1NTUAgPj4eBiNRpw9exYAkJOTg+rqaphMJoSEhCAmJgbV1dUAgLi4OJjNZrS2tkKv12PixImoq6tDT08PgoKCkJCQgMrKSof/hrZTrWdkZKChocHee8yYMfbPpqOioqBQKOznp0hLS0NzczM6OzuhVCqRlpZmX6YiIiKgUqkcere1tUGv10OhUCArKwtlZWX23sHBwaivrwcAJCcn49y5c7y9tVot6urqAABJSUnQ6/W8vcPDw1FbWwsASEhIQHd3t7332LFjUVlZiZCQEISGhiIqKsqht8lkQltbGwAgOzsbtbW1MBqNCAkJQWxsrH2ZjY2NhdVqdVhmT58+be+dmJjosMzKZDJ77/T0dDQ2NqKrqwsqlQrJycnM3i0tLTAYDIP2VqvVgy6z/XuHhYUhNDTUYZnt6OhAR0cHAgICkJOTg+PHj8NqtaKnpweBgYEOvQ0GA86dOzdgmR2sd09Pz6DLbGhoKKKjowddZgEgKysLp06dgtFoRHBw8IjXEbb1krvXESdOnIBGo/HpdYSttyfWERaLBTqdzufXEb29vW5fR9jWPWJfRzhz/JaojvY6c+YMEhIS8PXXX2PKlCn22x988EHs378f33zzzaCPs1qtqK6uhsFgwN69e/HUU0+hpKQEl112Gb7++mtMnToVZ86cQVxcnP0x119/PWQyGbZt2+bwXLajvdLS0hAQ8NOGMZVK5ZZDIcvKypCbmzvi55EiasNGfdioDxv1YaM+/HyljTNHe4lqy09kZCTkcvmAi6g1NzcjNjaW93EBAQHIyMgAAEyYMAHl5eVYvXo1LrvsMvvjmpubHSY/zc3NmDBhAu9zajQatx7qbpOWlub255QKasNGfdioDxv1YaM+/KTYRlT7/CiVSkycOBF79+6132a1WrF3716HLUFDsVqt9n10UlNTERsb6/CcOp0O33zzjVPP6S58R60RajMU6sNGfdioDxv14SfFNqLa8gMAy5Ytw6JFizBp0iRccskleOmll9DZ2Yni4mIAwMKFC5GQkIDVq1cDAFavXo1JkyYhPT0dRqMRn376KbZs2YJXX30VQN85BZYuXYpVq1YhMzMTqampeOyxxxAfH4+ioiKvvz6DweD1n+krqA0b9WGjPmzUh4368JNiG9FNfubNm4fW1lasWLECTU1NmDBhAnbu3Gnfka2urs5hX5zOzk7ceeedqK+vR1BQEHJycvDWW29h3rx59vs8+OCD6OzsxG233Ybz589j2rRp2LlzJ9Rqtddfn1Kp9PrP9BXUho36sFEfNurDRn34SbGNqHZ4FgNPXd7Cxmq1OkzeyE+oDRv1YaM+bNSHjfrw85U2Pn95CyljHbLv76gNG/Vhoz5s1IeN+vCTYhua/BBCCCHEr9Dkx8siIiKEHoJoURs26sNGfdioDxv14SfFNjT58TIhdrL2FdSGjfqwUR826sNGffhJsQ1NfrzMdrpuMhC1YaM+bNSHjfqwUR9+UmxDkx9CCCGE+JVhnefnP//5j0tPbruyOvlJSkqK0EMQLWrDRn3YqA8b9WGjPvyk2GZYk5/LLrsMMpnM/jXHcQ5f87FYLK6PTKLa29sRHBws9DBEidqwUR826sNGfdioDz8pthnW5Ofzzz93+NpoNOLBBx9EV1cXbrvtNmRnZwPoOxfAhg0bEBISgueee879o5UAnU4n9BBEi9qwUR826sNGfdioDz8pthnW5GfmzJkOXy9btgxKpRIHDx502Av8mmuuwV133YWZM2di586dmDVrlntHKwFDnXXSn1EbNurDRn3YqA8b9eEnxTYu7fD89ttv46abbhr08Lfg4GDcdNNNeOutt0Y8OCnKysoSegiiRW3YqA8b9WGjPmzUh58U27g0+ens7ERjYyPv9xsbG9HV1eXyoKSsrKxM6CGIFrVhoz5s1IeN+rBRH35SbOPS5KegoAAvv/wyduzYMeB7H3zwAV5++WUUFBSMeHCEEEIIIe7m0gd569atw+WXX47rrrsOcXFxyMjIAABUVVXhzJkzSE9Px9q1a906UKkICwsTegiiRW3YqA8b9WGjPmzUh58U27i05SchIQFHjx7FCy+8gPz8fDQ3N6O5uRl5eXl48cUXcfToUSQmJrp7rJIQGhoq9BBEi9qwUR826sNGfdioDz8ptnH5DM9qtRr33nsvdu7cifLycpSXl2Pnzp344x//iKCgIHeOUVJOnz4t9BBEi9qwUR826sNGfdioDz8pthnR8WtGoxFHjhxBS0sLpk6disjISHeNixBCCCHEI1ze8vPKK68gLi4OU6dOxbXXXosffvgBANDW1obIyEhs3LjRbYOUkjFjxgg9BNGiNmzUh436sFEfNurDT4ptXJr8bNq0CUuXLsUVV1yBjRs3guM4+/ciIyNx+eWX491333XbIKWko6ND6CGIFrVhoz5s1IeN+rBRH35SbOPS5GfNmjWYM2cO3nnnHVxzzTUDvj9x4kSUlpaOeHBSJMWFyF2oDRv1YaM+bNSHjfrwk2IblyY/lZWVuPLKK3m/Hx4ejrNnz7o8KCkLCHD5k0bJozZs1IeN+rBRHzbqw0+KbVx6RaNHj0ZbWxvv98vKyhAbG+vyoKQsJydH6CGIFrVhoz5s1IeN+rBRH35SbOPS5OfXv/41/v73v+P8+fMDvldaWooNGzbgN7/5zUjHJknHjx8XegiiRW3YqA8b9WGjPmzUh58U27g0+Vm1ahUsFgvy8/OxfPlyyGQyvPHGG7jxxhsxadIkREdHY8WKFe4eqyRYrVahhyBa1IaN+rBRHzbqw0Z9+EmxjUuTn/j4eBw+fBhXXHEFtm3bBo7jsGXLFnzyySdYsGABDh48SOf84TFq1CihhyBa1IaN+rBRHzbqw0Z9+EmxjcsnOYyOjsbrr7+O119/Ha2trbBarYiKipLkjlHuJMWFyF2oDRv1YaM+bNSHjfrwk2Ibl2Yqn376KSwWi/3rqKgoxMTE0MRnGOrq6oQegmhRGzbqw0Z92KgPG/XhJ8U2Ls1Wrr76asTExOC2227D3r17Jfl5ICGEEEKkyaXJz7/+9S/85je/wfbt2zF79mzExcXhrrvuwhdffOHu8UlOUlKS0EMQLWrDRn3YqA8b9WGjPvyk2MalyU9hYSE2btyI5uZmfPTRR5g9ezbefvttXHbZZUhMTMTSpUtx4MABd49VEgwGg9BDEC1qw0Z92KgPG/Vhoz78pNhmRDvpBAYG4uqrr8aWLVvQ0tKCHTt2YMaMGXj99dcxffp0d41RUs6dOyf0EESL2rBRHzbqw0Z92KgPPym2cdseygaDAS0tLWhubkZPT4/DxU4JIYQQQsRCxo1gltLR0YEdO3Zg27Zt+Pzzz9Hb24uLLroI8+bNw7x585Cenu7OsXqF2WzG/v37MX78eMjlcqGHQwghhJBhsFgsOHr0KGbOnAmFgn0mH5e2/GzZsgXXXHMNYmJi8Ic//AF1dXV45JFHUFZWhqNHj+KRRx7xyYmPN5w4cULoIYgWtWGjPmzUh436sFEfflJs49JJDhctWoS0tDT86U9/wrx58zBu3Dh3j0uyzGaz0EMQLWrDRn3YqA8b9WGjPvyk2MbpyQ/Hcfjvf/+L3NxcBAUFeWJMkqbVaoUewgAcx0FntKC714KgQDm0KjlkMpnXx+GNNmJ5ra4Q47IjJtSHjfqwUR9+Umzj9OTHZDJh8uTJWL16NR544AFPjEnSwsPDhR6CncFoxu6T7SgpbUWj3mS/PU6jRFFeFGZlhiNU5fIVUJzmyTZie62uENOyI0bUh436sFEfflJs4/Q+PyqVCrGxsVAqlZ4Yj+TV1tYKPQQAwKF6HW7YWor1BxvQdMFkAACa9CasP9iAG7aW4lC9zmtj8lQbMb5WV4hl2REr6sNGfdioDz8ptnFph+fFixfjzTffhMlkGvrORHQO1euwfFcVjGYrOAD9D/ez3WY0W7F8V5XoJwUs/vRaCSGEDI9L2/kvuugilJSUIC8vD4sXL0ZKSsqg+/9ce+21Ix6g1CQkJAj68w1GM57cUwOOGzgR6I/73/95ck8N3lmQ5/GPhdzdRsyv1RVCLztiR33YqA8b9eEnxTYureEXLFhg//8fe+yxQe8jk8kcrvxO+vT09GDUqFGC/fzdJ9vtW0GGw7ZVZPfJdszNj/bk0NzeRsyv1RVCLztiR33YqA8b9eEnxTYuTX4+//xzd4/Db5w9exYxMTGC/GyO41BS2jrsycCFSkpbUZQX5dEjo9zZRuyv1RVCLju+gPqwUR826sNPim1cmvzMnDnT3eMgXqAzWhyOdBouDkCj3gS90QKtWnwfBw3Gn14rIYQQ54zo2l5GoxEHDhzARx99hLa2NneNSdJycnIE+9ndvSP7GLJrhI8fijvbiP21ukLIZccXUB826sNGffhJsY3Lk59XXnkFcXFxmDZtGq699lr88MMPAIC2tjZERkZi48aNLg9q3bp1SElJgVqtxuTJk/Htt9/y3nfDhg2YPn06wsLCEBYWhoKCggH3X7x4MWQymcO/K664wuXxjUR1dbUgPxcAggJHdq2y4BE+fijubCP21+oKIZcdX0B92KgPG/XhJ8U2Lk1+Nm3ahKVLl+KKK67A//3f/zlcwT0yMhKXX3453n33XZcGtG3bNixbtgwrV67EkSNHMH78eBQWFqKlpWXQ++/btw8LFizA559/jgMHDiApKQmzZ89GQ0ODw/2uuOIKNDY22v9t3brVpfGNlJCnB9Cq5IjTKOHsniwy9J0MUKPy7ITAnW3E/lpdQaeWYKM+bNSHjfrwk2IblyY/a9aswZw5c/DOO+/gmmuuGfD9iRMnorS01KUBvfDCC7j11ltRXFyM3NxcrF+/HsHBwbxbkt5++23ceeedmDBhAnJycvD666/DarVi7969DveznZzR9i8sLMyl8Y1UaGioID8X6DsCrygvyqXHemMHYHe2EftrdYWQy44voD5s1IeN+vCTYhuXJj+VlZW48soreb8fHh6Os2fPOv28JpMJhw8fRkFBwU8DDAhAQUEBDhw4MKzn6OrqQm9v74DTce/btw/R0dHIzs7GkiVLhhyfXq+HTqez/zMajU6/nsFERwt7CPWszHCoFAHD3iIiA6BSBGBWpudPb+7uNmJ+ra4QetkRO+rDRn3YqA8/KbZx6XCW0aNHM3dwLisrQ2xsrNPP29bWBovFMuCQupiYGBw/fnxYz/HnP/8Z8fHxDhOoK664Atdeey1SU1NRVVWFRx55BFdeeSUOHDgAuXzwjzfy8/PR1dVl/7q4uBj33HMP4uLiUFVVZR8Xx3H2j+QyMzNRX1+P7u5uqNVqJCUl4eTJkwD6Fp6AgACcPHkSGo0G6enpaGpqQmdnJ1QqFVJSUlBRUQGg76NDpVKJM2fOAABSU1PR2toKg8GAwMBAZGRkoLy8HEDfRDMoKMj+MV9KSgra29uh0+kgl8uRnZ2N8vJycByH0aNHQ6PRYFGmHH8vt0IG9sn/bJOGxVlyWI1daDrXhfb2dgBAbm4uTpw4AbPZDI1Gg8jISNTU1AAA4uPjYTQa7RPMnJwcVFdXw2QyISQkBDExMfbPkOPi4mA2m9Ha2gq9Xo+JEyeirq4OPT09CAoKQkJCAiorK+29AaC5uRkAkJGRgYaGBnvvMWPG4MSJEwCAqKgoKBQKl16r3NqL06cbodfroVAokJWVhbKyMnvv4OBg1NfXAwCSk5Nx7tw53t5arRZ1dXUAgKSkJOj1epw/fx4ymQxjx45FRUUFLBYLtFotwsPD7aeRT0hIQHd3t7332LFjcezYMYSEhCA0NBRRUVEOvU0mk/09mZ2djdraWhiNRoSEhCA2Nta+zMbGxsJqtToss6dPn7b3TkxMdFhmZTKZvXd6ejoaGxvR1dUFlUqF5OTkAb0bGxsBAGlpaWhpaYHBYIBSqURaWpr9PRwREQG1Wj3oMtu/d1hYGEJDQ3H69GkAwJgxY9DR0YGOjg4EBAQgJycHx48fh9VqRU9PD7Kyshx6GwwGnDt3bsAyO1jvnp6eQZfZ0NBQREdHD7rMAkBWVhZOnToFo9GI4ODgEa8jmpqa7L3duY44fvw4NBrNsNYRF/bW6XS8y2xYWBhOnToFAEhMTERXl2fXEbbe7lxH2JZZ22vq7OwcdJlVqVQOvdva2kS3jqisrOz7iEoVgpBRYWg5U48QRd99R7KOsK17xL6OuHAXnKHIOGfu/T8333wz/v3vf+P777+HxWJBVFQU9uzZg8svvxylpaWYPHkybr75ZrzyyitOPe+ZM2eQkJCAr7/+GlOmTLHf/uCDD2L//v345ptvmI9/5pln8Nxzz2Hfvn0YN24c7/2qq6uRnp6OPXv24Fe/+pXD98xmM/bv34+0tDQEBPy0YUylUkGlUjn1egZTVlaG3NzcET/PSB2q1+HJPTUwmq0AHCcGtomAShGAFQWpmJTonSv6eqqNGF+rK8Sy7IgV9WGjPmy+3seTF2/2lTYWiwVHjx7FzJkzoVCwX6tLJVatWoXJkycjPz8f11xzDWQyGd544w1s3LgRH3zwAeLi4rBixQqnnzcyMhJyudw+i7Rpbm4eckvS888/j2eeeQZ79uxhTnyAvhlnZGQkKisrB0x+bDQaDe9WoZGIi4tz+3O6YlKiFu8syBv0zRL7vzfL7KwIhCi9t+Ovp9qI8bW6QizLjlhRHzbqw+bLffr/gXch28WbNx1qdPkPPF9uw8elyU98fDwOHz6MRx55BNu2bQPHcdiyZQs0Gg0WLFiAZ555BpGRkU4/r1KpxMSJE7F3714UFRUBgH3n5bvvvpv3cc899xyefvpp7Nq1C5MmTRry59TX1+Ps2bOC/Ac1m81e/5l8QlUKzM2PRlFeFPRGC7p6LQgOlEOjkguyw68n24jttbpCTMuOGFEfNurD5qt9bBdv5ruGoe0228WbVxWmOz0B8tU2LC6f5yc6Ohqvv/462tvb0dzcjMbGRpw7dw4bN24c0c5Ry5Ytw4YNG/DGG2+gvLwcS5YsQWdnJ4qLiwEACxcuxMMPP2y//7PPPovHHnsMGzduREpKCpqamtDU1ASDwQAAMBgMeOCBB3Dw4EHU1tZi7969mDNnDjIyMlBYWOjyOF1l+9xaTGQyGbRqBWI1KmjVCsEmA95oI5bX6goxLjtiQn3YqA+bL/Zx9uLN3P8u3mwwOjeZ8cU2Q3HL+fujovoOKTaZTOjt7UVISIjLzzVv3jy0trZixYoVaGpqwoQJE7Bz5077jmx1dXUO++K8+uqrMJlM+N3vfufwPCtXrsTjjz8OuVyOH374AW+88QbOnz+P+Ph4zJ49G0899ZRb9uEhhBBChCC1izd7k0s7PL/77rv45ptv8OKLL9pve+KJJ/D000+D4zhcffXV2LJli0+eG8C2w/P48eM9ss+P2Wweckcsf0Vt2KgPG/Vhoz5svtaH4zgsfq/M6WsYytC3r+Pm63OHveXbV9o4s8Ozyyc57OzstH/99ddf44knnkBhYSHuu+8+7Ny5E08//bQrTy15tsNCR4LjOHT0mNGkN6Kjx+zU4X1i5o42UkZ92KgPm6f7+Pp6ydeWH3dcvHm4fK3NcLg0lauqqsKiRYvsX7/zzjuIjY3Fhx9+CIVCAavVig8++ACrV69220ClYiQnS/TkoYxi4K4TSUoV9WGjPmye6iOV9ZKvLT/uuHizVj28/y6+1mY4XNryYzQaoVar7V9/9tlnuPLKK+2bmXJzc+0neCKOgoODXXrcoXodbthaivUHG9DUb7ZvO5Txhq2lOFSvc8cwBeFqG39BfdioD5sn+khpveRry483L97sa22Gw6XJT2pqKvbs2QMAOHToECorKx2ukt7c3OyT+/t4gyuH19sOZbTt2NZ/Y7LtNtuhjL6wohmMFM8l4U7Uh436sLm7j9TWS762/Hjz4s2+1mY4XJr83H777Xjvvfcwbtw4zJ49G4mJibj66qvt3//qq6+Ql5fntkFKie0U4sPlrUMZxcDZNv6G+rBRHzZ39pHiesnXlh9vXrzZ19oMh0uTn3vuuQevvfYa0tPTMWfOHHz22WcICgoCALS3t6OpqQm///3v3TpQfzWSQxkJIcQTaL0kDlK7eLM3ubwX2q233opbb711wO3h4eE4dOjQiAYlZf0v2srCcRxKSluHvYK5UElpq9Oze6E508YfUR826sPmrj5SXS/54vITqlJgRUEqlu+qAobYCicDIJMBKwpSnd4B3RfbDGVEu+BbLBYcPnzYfpXZlJQUTJw40SPnx5EKZw7/dMehjMPdm18MfO3QWG+jPmzUh81dfaS6XvLV5WdSoharCtM9evFmX23D4vLlLTZv3ozExERMmTIF8+fPx/z58zFlyhQkJCRg48aN7hyjpLS0tAz7vu44lNGXONPGH1EfNurD5q4+Ul0v+fLyY7t48x2/SECsRunwvViNEnf8IgFbb8h3aeID+HYbPi5Nv1977TUsWbIEEyZMwOOPP46srCwAQEVFBV577TXceuutMJlMuOOOO9w6WH/jzUMZCSFkOGi9JE5SuHizN7l0eYu0tDQkJSVhz549CAwMdPheb28vLr/8cjQ0NKC6utptA/UWT1/eore3d0AzPrbTlzfpTU59vu7K6cvFwJk2/oj6sFEfNnf1kep6iZYffr7SxuOXt2hqasL1118/aIzAwEDMnz8fzc3Nrjy15Dlz8kdvHsooBnRiTDbqw0Z92NzVR6rrJVp++EmxjUuTn4svvhgnTpzg/f6JEycwYcIEV8ckad3d3U7d358OZXS2jb+hPmzUh82dfaS4XqLlh58U27g0+Vm7di3ee+89vPzyyw5Ruru78eKLL+K9997D3/72N7cNUkouvCzIcNgOZZTJMOSKZiSHMoqBs238DfVhoz5s7uwjxfUSLT/8pNhmWPv8jBs3bsBt7e3taGxshEKhQHx8PADgzJkzMJvNiIuLQ0REBI4ePer+EXuYmPb5udChep1HD2UUA1/5XFko1IeN+rB5oo+U1ku0/PDzlTbO7PMzrGl4eHj4gM9pIyIikJmZ6XBbSkqKcyP1QydPnkRubq7Tj7MdyjjY1ZNj/3f15NlZEQhR+u6RFK628RfUh436sHmij5TWS7T88JNim2FNfvbt2+fhYZDhoEMZCSFiQ+sl4ovE+wGsREVHR4/4OWQyGbRqhSjPkjoS7mgjZdSHjfqwebqPr6+XaPnhJ8U2Li+lFosFb731Fv75z3/i1KlTAIDk5GRcffXV+P3vf0+XuOAREODySbUlj9qwUR826sNGfdioDz8ptnHpFXV0dGDq1Km4+eab8dlnn6G3txe9vb3YvXs3iouLMW3aNOh0OnePVRKampqEHoJoURs26sNGfdioDxv14SfFNi5Nfh599FEcPnwYa9euRWtrK44cOYIjR46gpaUFf/vb33Do0CE8+uij7h4rIYQQQsiIuXR5i4SEBPzud7/Dyy+/POj3//jHP2L79u04c+bMiAfobZ4+1N1oNEKlUrn9eaWA2rBRHzbqw0Z92KgPP19p4/HLW5w9exbZ2dm838/JyUF7e7srTy15Utx86C7Uho36sFEfNurDRn34SbGNS5OfjIwMfPzxx7zf//jjj5Genu7yoKSss7NT6CGIFrVhoz5s1IeN+rBRH35SbOPS5OfOO+/EZ599hl//+tf47LPPUFtbi9raWuzatQtXXXUVdu/ejbvvvtvdY5UEX9h0KBRqw0Z92KgPG/Vhoz78pNjGpUPd77zzTrS0tOCZZ57Brl27HL4XGBiIFStWYMmSJW4ZoNTQWbD5URs26sNGfdioDxv14SfFNi4fvP/444+jvr4eb7/9Nv7yl7/gL3/5C95++23U19dj5cqV7hyjpFRUVAg9BNGiNmzUh436sFEfNurDT4ptRnQqzsjISMyfP99dYyGEEEII8TjpnbZR5CIjI4UegmhRGzbqw0Z92KgPG/XhJ8U2NPnxMqVSKfQQRIvasFEfNurDRn3YqA8/KbahyY+X+eKJH72F2rBRHzbqw0Z92KgPPym2ockPIYQQQvyK05Mfo9GIjz/+GD/88IMnxiN5qampQg9BtKgNG/Vhoz5s1IeN+vCTYhunJz9KpRLXXXcdvv76a0+MR/JaW1uFHoJoURs26sNGfdioDxv14SfFNk5PfmQyGTIzM9HW1uaJ8UiewWAQegiiRW3YqA8b9WGjPmzUh58U27i0z88jjzyCv/3tb5I88ZGnBQYGCj0E0aI2bNSHjfqwUR826sNPim1cOsnhwYMHERERgfz8fFx22WVISUlBUFCQw31kMhlefvlltwxSSjIyMoQegmhRGzbqw0Z92KgPG/XhJ8U2Mo7jOGcfFBAw9AYjmUwGi8Xi0qCEZDabsX//fowfPx5yudztz19WVobc3Fy3P68UUBs2KfXhOA46owXdvRYEBcqhVckhk8lG9JxS6uMJ1IeN+vDzlTYWiwVHjx7FzJkzoVCwt+24tOXHarW6NDBCiH8zGM3YfbIdJaWtaNSb7LfHaZQoyovCrMxwhKpGdNUdQggZEq1lvCw8PFzoIYgWtWHz9T6H6nV4ck8NjOaBfzw16U1Yf7ABmw41YkVBKiYlap1+fl/v42nUh4368JNimxGd5PDgwYNYvXo17rvvPpw8eRIA0NXVhSNHjkhy73B36L9vFPkJtWHz5T6H6nVYvqsKRrMVHID+n7XbbjOarVi+qwqH6nVO/wxf7uMN1IeN+vCTYhuXJj8mkwnXXnstpk6dikcffRSvvPIKTp8+3feEAQGYPXs27ezMo6GhQeghiBa1YfPVPgajGU/uqQHHDZz09McB4DjgyT01MBjNTv0cX+3jLdSHjfrwk2IblyY/jz32GP7xj3/g1VdfRUVFBS7cZ1qtVuO6667DRx995LZBEkJ81+6T7fYtPsNh2wK0+2S7J4dFCPFjLk1+tm7diiVLluC2224b9LPAsWPHorq6esSDk6KUlBShhyBa1IbNF/twHIeS0tZhT3wuVFLaCmcORhVrH47j0NFjRpPeiI4es1OvyZ3E2kcsqA8/KbZxafLT0tKCiy66iPf7crkcXV1dLg9q3bp1SElJgVqtxuTJk/Htt9/y3nfDhg2YPn06wsLCEBYWhoKCggH35zgOK1asQFxcHIKCglBQUGDfR8nb2tvpr1k+1IbNF/vojBaHo7qGiwPQqDdBbxz+6TLE1sdgNOPDH1uw+L0yXPfWMSzc1ve/i98rw4c/tjj9sd5Iia2P2FAfflJs49LkJykpCcePH+f9/ldffeXySZG2bduGZcuWYeXKlThy5AjGjx+PwsJCtLS0DHr/ffv2YcGCBfj8889x4MABJCUlYfbs2Q6fUT733HN45ZVXsH79enzzzTcICQlBYWEhenp6XBrjSOh0zu/I6S+oDZsv9unuHdm5vrqceLyY+hyq1+GGraVYf7ABTf0mf7Yj227YWurSjt2uElMfMaI+/KTYxqXJzw033IDXXnsNBw4csN9mO0HZhg0b8N5772HhwoUuDeiFF17ArbfeiuLiYuTm5mL9+vUIDg7Gxo0bB73/22+/jTvvvBMTJkxATk4OXn/9dVitVuzduxdA31afl156CcuXL8ecOXMwbtw4vPnmmzhz5gxKSkpcGuNIeOLEiVJBbdh8sU9Q4MjGHOzE48XSxxtHtrlCLH3Eivrwk2IblyY/jz76KC699FLMmDEDv/zlLyGTyXDfffdhzJgxuP3223HFFVfgvvvuc/p5TSYTDh8+jIKCgp8GGBCAgoICh4kWS1dXF3p7e+37ItXU1KCpqcnhOUeNGoXJkycP+zndKTs72+s/01dQGzZf7KNVyRGnUcLZczfL0HfiQ41q+CtdMfTx1pFtrhBDHzGjPvyk2MalyY9SqcTOnTuxadMmpKWlIScnB0ajEePGjcPmzZvxySefuDRTbGtrg8ViQUxMjMPtMTExaGpqGtZz/PnPf0Z8fLx9smN7nLPPqdfrodPp7P+MRqMzL4VXeXm5W55HiqgNmy/2kclkKMqLcumxRXlRTl3yQgx9xHxkmxj6iBn14SfFNi6f4Vkmk+HGG2/EjTfe6M7xjMgzzzyDd999F/v27YNarR7Rc+Xn5zvstF1cXIx77rkHcXFxqKqqAtA3geI4zr4/UmZmJurr69Hd3Q21Wo2kpCT7jtXR0dEICAiATqdDWVkZ0tPT0dTUhM7OTqhUKqSkpKCiogIAEBkZCaVSiTNnzgAAUlNT0draCoPBgMDAQGRkZNgXxvDwcAQFBdn3cUpJSUF7ezt0Oh3kcjmys7NRXl4OjuMwevRoaDQa+zmZxowZA51Oh/Pnz0Mmk2Hs2LGoqKiAxWKBVqtFWFgYTp06BQBITExEV1eXfce33NxcnDhxAmazGRqNBpGRkaipqQEAxMfHw2g04uzZswCAnJwcVFdXw2QyISQkBDExMfajAePi4mA2m9Ha2gq9Xg+z2Yy6ujr09PQgKCgICQkJqKystPcGgObmZgB9F9traGiw9x4zZgxOnDgBAIiKioJCoUBjYyMAIC0tDc3Nzejs7IRSqURaWpp9v7WIiAioVCqH3m1tbdDr9VAoFMjKykJZWZm9d3BwMOrr6wEAycnJOHfuHG9vrVaLuro6AH37yun1et7e4eHhqK2tBQAkJCSgu7vb3nvs2LEwGAwoKytDaGgooqKiHHqbTCa0tbUB6Psrrba2FkajESEhIYiNjbUvs7GxsbBarQ7L7OnTp+29ExMTHZZZmUxm752eno7GxkZ0dXVBpVIhOTmZ2bulpQUGgwHp8kCoFAHDnhTIAAQGAEnWVjQ2WhAaGuqwzHZ0dKCjowMBAQHIycnB8ePHYbVa0dPTA4PB4NDbYDDg3LlzA5bZwXr39PQMusyGhoYiOjp60GUWALKysnDq1Cn09PRg+1Gz00e2cQDe+64B02JkkMvl9j/K3L2OsK17fHkdYevtiXWExWLBqVOnfHodUVlZid7eXrevI2zrHk+tIwbrrVarB/291r93WFiYfR3hzJGULl3Y9LXXXsOMGTMwduxYZx/KZDKZEBwcjO3bt6OoqMh++6JFi3D+/HnmuYOef/55rFq1Cnv27MGkSZPst1dXVyM9PR3fffcdJkyYYL995syZmDBhwoCTMdoubJqWluZwAVeVSgWVSjXi13jmzBnEx8eP+HmkiNqw+XIf234wQ30cJAMgkwGrCtOdvsSF0H06esy47q1jLj9++40XQav23BWHhO4jdtSHn6+0cebCpi597LVkyRLk5+cjKioKc+fOxQsvvID//ve/I77gqVKpxMSJE+07KwOw77w8ZcoU3sc999xzeOqpp7Bz506HiQ/QN0OPjY11eE6dTodvvvmG+ZwajQZardb+zx0TH9vzksFRGzZf7jMpUYtVhelQKQL6Jjj9vm+7TaUIcGniAwjfx5tHtrlC6D5iR334SbGNS5OfpqYmbNu2Db///e9RV1eHBx98EL/4xS8wevRozJ49G0899RT27dvn0oCWLVuGDRs24I033kB5eTmWLFmCzs5OFBcXAwAWLlyIhx9+2H7/Z599Fo899hg2btyIlJQUNDU1oampyX5tMZlMhqVLl2LVqlX4+OOPcezYMSxcuBDx8fEOW5e8xbY5mQxEbdh8vc+kRC3eWZCHO36RgFiN0uF7sRol7vhFArbekO/SxAcQvo83j2xzhdB9xI768JNiG5e2sUZHR+N3v/sdfve73wHo2zn466+/xhdffIHt27fj8ccfh0wmg9ns/BEM8+bNQ2trK1asWIGmpiZMmDABO3futH+WW1dX5/Bx1KuvvgqTyWQfi83KlSvx+OOPAwAefPBBdHZ24rbbbsP58+cxbdo07Ny5c8T7BRFCnBOqUmBufjSK8qKgN1rQ1WtBcKAcGpXcqZ2bxch2ZFuT3uTUfj8y9E3+nDmyjRAyMi7t83OhqqoqfPHFF/jiiy/wn//8B1VVVQgJCcGUKVPw2WefuWucXmPb52f8+PEeObeBwWBAaGio259XCqgNG/VhE0OfD39swfqDDU5Pfu74RQLm5kd7algAxNHHGRzHQWe0oLvXgqBAObQeniD7Wh9v8pU2zuzz49KWn7/97W/2CU9TUxPCw8Mxbdo0LFmyBNOnT8fPfvYzSZ4UyR10Op1PLERCoDZs1IdNDH1mZYZj06FGp45sUykCMCtz4DUS3U0MfYbDYDRj98l2lJS2OlwaJU6jRFFeFGZlhiNU5f4dw32ljxCk2MalfX7++Mc/YseOHZgxYwb++9//oq2tDSUlJVi2bBl+/vOf08SH4fz580IPQbSoDRv1YRNDn1CVAisKUiGTDdypuz/bkW0rClI98su8PzH0GYqQlwXxhT5CkWIblyY/d911F/Lz87F9+3ZMnToV06ZNw8MPP4xPP/0UHR0d7h6jpPj6fg2eRG3YqA+bWPp448g2V4ilDx+hLwsi9j5CkmKbEe3zo9Pp8OWXX+LLL7/EF198gf/+978wm83Iz8/H9OnTsXbtWneO1Ss8vc8PIcQ/DPXxzeysCIQoaR0D9LW6YWup0x8XvrMgzytbzYhvcGafnxHv8GxTU1OD3bt348UXX0RFRQVkMhksFs+et8ITPD35qaiokOR1UtyB2rBRHzax9uE4ThRHtom1DyCOHcXF3EdovtLG4zs8A0BZWZn9CK8vvvjCfhrq+Ph4zJ8/H9OnT3f1qSXNFyeE3kJt2KgPm1j7yGQyaNUKj569eTjE2ofjOJSUtjp9WRAAKCltdfoacHzE2kcMpNjGpXdjZGQkzp07B47jkJOTgyuvvBLTpk3D9OnTkZKS4uYhSotW653P930RtWGjPmzUh02sfXRGi8PHgsPFAWjUm6A3WtwysRRrHzGQYhuXlphFixZh+vTpmDZtGiIjI909JkkLCwsTegiiRW3YqA8b9WETax93XBbEHZMfsfYRAym2celorzVr1qCoqIgmPi6wXQGZDERt2KgPG/VhE2sfsVwWRKx9xECKbUY0Xd6/fz/++c9/2sMkJyfjqquuwsyZM90yOEIIIdJGlwUhQnBp8mMymbBgwQKUlJSA4ziMHj0aQN+JkNasWYO5c+di69atCAwMdOdYJSExMVHoIYgWtWGjPmzUh02sfWQyGYryorD+YIPTj3XXzs6AePuIgRTbuPSx1xNPPIEPP/wQf/rTn9DY2Ij29na0t7ejqakJ999/P3bs2IEnn3zS3WOVhK6uLqGHIFrUho36sFEfNjH3mZUZbj8p5HB44rIgYu4jNCm2cWny884772DRokV47rnn7FdbB/qu9v7ss89i4cKF2LJli9sGKSXt7e1CD0G0qA0b9WGjPmxi7iOGy4KIuY/QpNjGpclPY2MjJk+ezPv9yZMno6mpyeVBEUII8S9ivSwIkSaXzvCckZGBSZMm4d133x30+/Pnz8ehQ4dQWVk54gF6G13eghBChEOXBSGucuYMzy5t+Vm0aBHee+893HHHHaioqIDFYoHVakVFRQWWLFmC999/H4sXL3blqSXvxIkTQg9BtKgNG/Vhoz5svtInVKXA3PxobL4+F9tvvAhvzuv7383X52JufrTHJj6+0kcIUmzj0gemjzzyCKqqqvD3v/8dGzZsQEBA3xzKarWC4zgsWrQIjzzyiFsHKhVms1noIYgWtWGjPmzUh83X+nj7siC+1sebpNjGpaVKLpdj8+bNWLZsGf75z3+irq4OQN95fn79619j3Lhxbh2klGg0GqGHIFrUho36sFEfNurDRn34SbGNU5Ofnp4efPTRR6ipqUFkZCSuuuoqPPzww54amyTRWbH5DdaG4zjojBZ091oQFCiHVqArYnvDUK+Vlh026sNGfdioDz8pthn25KelpQWXXnopampqYNtHOjg4GCUlJSgoKPDYAKWmpqYGubm5Qg9DlC5sM9ROj7Myw916mKuQhvtaadlhoz5s1IeN+vCTYpth7/D81FNPoba2Fvfddx/+8Y9/4KWXXkJQUBBuv/12T46P+KFD9TrcsLUU6w82oKnf1Z6b9CasP9iAG7aW4lC9TqARuo8/vVZCCBGLYf/p/Nlnn2HhwoV4/vnn7bfFxMTghhtuQEVFBbKzsz0yQKmJj48XegiiFR8fj0P1OizfVQWOw6DX+bHdZjRbsXxXlU+f78PZ1/rQ1BhI628v96L3Fhv1YaM+/KTYZthbfurq6jBt2jSH26ZNmwaO49Dc3Oz2gUmV0WgUegii1W7oxpN7angnAxfiAHAc8OSeGhiMvnckgsFodvq1rjnQ4pOv1VvovcVGfdioDz8pthn25MdoNEKtVjvcZvtaiofBecrZs2eFHoJo7TreCqPZOuwrO3Po2yqy+6TvnXp998l251+rhfPJ1+ot9N5ioz5s1IefFNs4tcdobW0tjhw5Yv+6o6MDAHDy5En7ld0v9LOf/WxkoyN+g+M4/KfRMuzJwIVKSlvdenVnT+M4DiWlrX7xWom4+dPRlIRcaNiXtwgICBj0TcFx3IDbbbdZLBb3jNKLPH15C6vVaj8pJPlJR48Z1711zOXHb7/xIq+dDG2k/Om1ehO9t9gu7ONPR1MOFy0//HyljTOXtxj20r1p06YRD4wA1dXVyMjIEHoYotPdO7KJclevxWcmBP70Wr2J3ltstj6H6nV4ck8NjGbrgPvYjjDcdKgRKwpSffZgAlfQ8sNPim2GvQZdtGiRJ8fhN0wm09B38kNBgSPbyhY8wsd7kz+9Vm+i9xabyWTyq6MpnUXLDz8pthH/diyJCQkJEXoIoqRVyREVFABn9zaQoW9TvUblOxMCrUqOOI3SL16rN9F7i02mDPKboyldQcsPPym2ocmPl8XExAg9BFGSyWSYkxfl0mN9bQdgmUyGIh97rRzHoaPHjCa9ER09ZgxzV0GvsI0NIWGiG5uYHDOo/OZoSlfQunlwHMchaFSEKN/7I0E7DnhZdXW15E4T7i5pAeegUgQMewUtA6BSBGBWZrinh+Z2szLDselQo1OvNTAAXn+tYt4xVsxjExuO4/DhsWY6wpCB1s2OpP7+oi0/RDSCFTKsKEiFTIYhPxKSAZDJgBUFqT75BgxVKZx+rTfnBHr1tYr50htiHpsY6YwWnHXhPHUcgEa9CXqj7x25S1znD+8vmvx4WVxcnNBDEK24uDhMStRiVWE6VIq+/X/6Twxst6kUAT6/M6azr/WXYxO8NjbbjrG2LVP9txjYbrPtGOvNlaCYxyZW7jjCUOqcWTeL+WPgkfKX95fv/cns4+hs2PxsbSYlavHOgrxBN7nG/m+T6+ysCIQofX/HX2dea2trq1fG5OylN/C/HWPfWZDn8S1TYh6bmNERhkMbzrpZ6h8F+dP7y7dGKwGtra2IinJtZ1epu7BNqEqBufnRKMqLgt5oQVevBcGBcmgkeAba4b5Wby07Ll164387xs7Nj/bk0EQ9NjHTquSIUAHtxqF/qV1Ihr5JuD8cYTjU+8sfzo/kT+8v+tiLiJpMJoNWrUCsRgWtWiG5ic+FxPBaR3rpDU9u/hfz2MROJpNhRpxrExh/2Nl5KP7wUZC/vb9o8uNlWVlZQg9BtKgNmzf66IwWh835w+WNHWPFPDZfsGBKtn3/suHw5aMpXcH3/nL2oyBfPT+Sv72/aPLjZXV1dUIPQbSoDZs3+oh5x1gxj80XtDc3+M3RlK7ge3+N5KMgX+Jv7y+a/HhZT0+P0EMQLWrD5o0+Yt4xVsxj8wU9PT1+dTSlswZ7f/nTR0H+9v6iyY+XBQUFCT0E0aI2bN7oI+ZLb4h5bL7AtvzYjjC84xcJiNUoHe4Tq1Hijl8kYOsN+X418QEGf3/500dB/vb+8o/tmSKSkOC9c7X4GmrD5o0+tktvrD/Y4PRjPb1jrJjH5gsuXH786WjK4Rrs/eWOj4K0at/4Netv7y/a8uNllZWVQg9BtKgNm7f6zMoMF+2OsWIem9gNtvyI4QhDsRisj799FORP7y+a/BBCHLhy6Q1v7Rgr5rER6fG3j4L86f1Fkx8voysH86M2bN7sI+YdY8U8NjGj9xfbYH1sHwW5whc/CgL85/3le9M1QohXiPkyI2IeG5GWWZnh2HSocdiHu/vyR0E2/vD+Et3kZ926dfjrX/+KpqYmjB8/HmvXrsUll1wy6H1LS0uxYsUKHD58GKdOncKLL76IpUuXOtzn8ccfxxNPPOFwW3Z2No4fP+6pl8DU3NyMiIgIQX622FEbNiH6iHnH2P5jO1pWgfG52aIYmxjR+4uNr4/to6Dlu6qAIU506OsfBV3owvfXoR/KkJSWIZr3vjuI6mOvbdu2YdmyZVi5ciWOHDmC8ePHo7CwEC0tLYPev6urC2lpaXjmmWcQGxvL+7x5eXlobGy0//vyyy899RIIkSQx7xhrG1uEWia6sRFp8JePggYjk8kQEigT5Xt/JGSciM7CNHnyZPz85z/H3/72NwCA1WpFUlIS7rnnHjz00EPMx6akpGDp0qWDbvkpKSnB999/P6wxmM1m7N+/H+PHj4dc7v5NeiaTCUqlcug7+iFqw0Z92KgPG/VhG06foa7q7usfBfHxlWXHYrHg6NGjmDlzJhQK9pY30WyXM5lMOHz4MB5++GH7bQEBASgoKMCBAwdG9NwnT55EfHw81Go1pkyZgtWrV2PMmDHMx+j1egQE/LRhTKVSQaVSjWgcANDQ0IDU1NQRP48UURs26sNGfdioD9tw+oj5Y2BPkuKyI5rJT1tbGywWy4A97mNiYka0f87kyZOxefNmZGdno7GxEU888QSmT5+OH3/8ERqNhvdx+fn56Orqsn9dXFyMe+65B3FxcaiqqrKPjeM4+8dymZmZqK+vR3d3N9RqNZKSknDy5EkAQHR0NAICAtDS0oLu7m6kp6ejqakJnZ2dUKlUSElJQUVFBQAgMjISSqUSZ86cAQCkpqaitbUVBoMBgYGByMjIQHl5OQAgPDwcQUFBaGjoOzFVSkoK2tvbodPpIJfLkZ2djfLycnAch9GjR0Oj0eD06dMAgDFjxkCn0+H8+fOQyWQYO3YsKioqYLFYoNVqERYWhlOnTgEAEhMT0dXVhfb2vuvV5Obm4sSJEzCbzdBoNIiMjERNTQ0AID4+HkajEWfPngUA5OTkoLq6GiaTCSEhIYiJiUF1dTUAIC4uDmazGa2trdDr9UhKSkJdXR16enoQFBSEhIQE+/k3bMtGc3MzACAjIwMNDQ323mPGjMGJEycAAFFRUVAoFGhsbAQApKWlobm5GZ2dnVAqlUhLS7MvVxEREVCpVA6929raoNfroVAokJWVhbKyMnvv4OBg1NfXAwCSk5Nx7tw53t5ardZ+zaCkpCTo9Xre3uHh4aitrQXQd8K17u5ue++xY8eira0N3d3dCA0NRVRUlENvk8mEtrY2AH37tNXW1sJoNCIkJASxsbH2ZTY2NhZWq9VhmT19+rS9d2JiosMyK5PJ7L3T09PR2NiIrq4uqFQqJCcnM3u3tLTAYDAM2lutVg+6zPbvHRYWhtDQUIdltqOjAx0dHQgICEBOTg6OHz8Oq9WKnp4eREVFOfQ2GAw4d+7cgGV2sN49PT2DLrOhoaGIjo4edJkF+i6IeerUKRiNRgQHB494HdHU1GTv7c51hG3d48vrCFtvT6wjLBYLTp06Nex1RMf/1hE6hQJakawjKisr0dvb6/Z1hG3dI/Z1hDMfZInmY68zZ84gISEBX3/9NaZMmWK//cEHH8T+/fvxzTffMB/P97FXf+fPn0dycjJeeOEF/OEPfxjwfdvHXmlpaR7Z8lNdXY20tLQRP48UURs26sNGfdioDxv14ecrbXzyY6/IyEjI5XL7DNKmubmZuTOzs0aPHo2srKwhz5ar0Wg8ss/PUB+3+TNqw0Z92KgPG/Vhoz78pNhGNEd7KZVKTJw4EXv37rXfZrVasXfvXoctQSNlMBhQVVWFuLg4tz2nM2ybAMlA1IaN+rBRHzbqw0Z9+EmxjWi2/ADAsmXLsGjRIkyaNAmXXHIJXnrpJXR2dqK4uBgAsHDhQiQkJGD16tUA+naStn3uZzKZ0NDQgO+//x6hoaHIyMgAANx///245pprkJycjDNnzmDlypWQy+VYsGCBMC+SEEIIIYIS1eRn3rx5aG1txYoVK9DU1IQJEyZg586d9p3Y6urqHPbDOXPmDC6++GL7188//zyef/55zJw5E/v27QMA1NfXY8GCBTh79iyioqIwbdo0HDx4EFFRrp2yfKSE+rm+gNqwUR826sNGfdioDz8pthHV5AcA7r77btx9992Dfs82obFJSUkZcu/ud999111Dc4uhdsLyZ9SGjfqwUR826sNGffhJsY1o9vnxF7bD/MhA1IaN+rBRHzbqw0Z9+EmxDU1+CCGEEOJXaPLjZb5wrgShUBs26sNGfdioDxv14SfFNjT58bL+5zEiP6E2bNSHjfqwUR826sNPim1o8uNlnZ2dQg9BtKgNG/Vhoz5s1IeN+vCTYhua/HiZL1wZVyjUho36sFEfNurDRn34SbENTX68TIqfnboLtWGjPmzUh436sFEfflJsQ5MfLxvJFeqljtqwUR826sNGfdioDz8ptqHJDyGEEEL8Ck1+vCwiIkLoIYgWtWGjPmzUh436sFEfflJsQ5MfL1OpVEIPQbSoDRv1YaM+bNSHjfrwk2Ibmvx42ZkzZ4QegmhRGzbqw0Z92KgPG/XhJ8U2NPkhhBBCiF+hyY+XpaamCj0E0aI2bNSHjfqwUR826sNPim1o8uNlbW1tQg9BtKgNG/Vhoz5s1IeN+vCTYhua/HiZXq8XegiiRW3YqA8b9WGjPmzUh58U29Dkx8sUCoXQQxAtasNGfdioDxv1YaM+/KTYRsZxHCf0IMTEbDZj//79GD9+PORyudDDIYQQQsgwWCwWHD16FDNnzhxywkZbfrysrKxM6CGIFrVhoz5s1IeN+rBRH35SbCO9bVmEEEJIPxzHQWe0oLvXgqBAObQqOWQymdDDIgKhyY+XhYeHCz0E0aI2bM728beVPS0/bP7ax2A0Y/fJdpSUtqJRb7LfHqdRoigvCrMywxGqUvhtn+GQYhua/HhZcHCw0EMQLWrDNtw+w13ZSw0tP2z+2OdQvQ5P7qmB0Wwd8L0mvQnrDzZg06FGrChIRZbW//oMlxSXHdrnx8vq6+uFHoJoURu24fQ5VK/DDVtLsf5gA5oumPgAP63sb9haikP1Ok8NUzC0/LD5W59D9Tos31UFo9kKDkD/I3tstxnNVizfVYXdx+q8P0gfIcVlhyY/hEiEsyt7KU6ACAH6tn4+uacGHDfwfdAfB4DjgE0VvTAYzd4YHhEBmvx4WXJystBDEC1qw8bq48rK/sk9NZJa2dPyw+ZPfXafbLf/ETAcHIBea9/jyEBSXHZo8uNl586dE3oIokVt2Fh9XFnZG81WSa3saflh85c+HMehpLR12O+FC5WUtoJOfTeQFJcdmvx4mU5HHzXwoTZsfH1oZd+Hlh82f+mjM1ocdvQfLg5Ao94EvdHi/kH5OCkuOzT58TI6azQ/asPG14dW9n1o+WHzlz7dvSNbnrtG+HgpkuKyQ5MfL8vOzhZ6CKJFbdj4+tDKvg8tP2z+0icocGS/qINH+HgpkuKyQ5MfLysvLxd6CKJFbdj4+tDKvg8tP2z+0kerkiNOo4Qrp/OM0yihUUnj/eBOUlx2aPLjZVLZv8ITqA0bXx9XV/YySGtlT8sPm7/0kclkKMqLcumxRXlRkj4LuqukuOzQ5MfLRo8eLfQQRIvasPH1oZV9H1p+2Pypz6zMcKgUAcP+g0AGQCWXYVam9C7j4A5SXHZo8uNlWq1W6CGIFrVhY/VxaWWvCJDUyp6WHzZ/6hOqUmBFQSpkMgz5npABkMmAP0+Lk+RlX9xBissOTX68rK6OTqHOh9qwsfq4srJfUZAqqZU9LT9s/tZnUqIWqwrT7X8U9H9f2G5TKQKwqjAd4b1nvT9IHyHFZYcmP4RIhLMr+0mJ0vtrjpALTUrU4p0FebjjFwmI1SgdvherUeKOXyRg6w359F7wQ9L5s89HJCUlCT0E0aI2bMPpY1vZD3ZV99j/XdV9dlYEQpTS2Mn5QrT8sPlrn1CVAnPzo1GUFwW90YKuXguCA+XQqOQO+7v5a5/hkGIbmvx4mV6vh0ajEXoYokRt2IbbZ7gre6mh5YfN3/vIZDJo1Qpo1YP/2vP3PixSbEMfe3nZ+fPnhR6CaFEbNmf72Fb2sRoVtGqFpCc+AC0/Q6E+bNSHnxTb0OTHy6T+C2gkqA0b9WGjPmzUh4368JNiGxknxbMXjYDZbMb+/fsxfvx4SV7PhBBCCJEii8WCo0ePYubMmVAo2Hv10JYfL6uoqBB6CKJFbdioDxv1YaM+bNSHnxTb0OTHyywWaVxE0hOoDRv1YaM+bNSHjfrwk2Ibmvx4mRTPlOku1IaN+rBRHzbqw0Z9+EmxDU1+vCw8XDqXE3A3asNGfdioDxv1YaM+/KTYRnSTn3Xr1iElJQVqtRqTJ0/Gt99+y3vf0tJS/Pa3v0VKSgpkMhleeumlET+np9XW1gr2s8WO2rBRHzbqw0Z92KgPPym2EdXkZ9u2bVi2bBlWrlyJI0eOYPz48SgsLERLS8ug9+/q6kJaWhqeeeYZxMbGuuU5CSGEECJtopr8vPDCC7j11ltRXFyM3NxcrF+/HsHBwdi4ceOg9//5z3+Ov/71r5g/fz5UKpVbntPTEhISBPm5voDasFEfNurDRn3YqA8/KbYRzeTHZDLh8OHDKCgosN8WEBCAgoICHDhwQDTPOVLd3d2C/FxfQG3YqA8b9WGjPmzUh58U24hm8tPW1gaLxYKYmBiH22NiYtDU1OT159Tr9dDpdPZ/RqPRpTH0197e7pbnkSJqw0Z92KgPG/Vhoz78pNiGLmzKIz8/H11dXfavi4uLcc899yAuLg5VVVUA+iZRHMfZ9x/KzMxEfX09uru7oVarkZSUhJMnTwIAoqOjERAQAL1ej7KyMqSnp6OpqQmdnZ1QqVRISUmxn0gqMjISSqUSZ86cAQCkpqaitbUVBoMBgYGByMjIQHl5OYC+vfCDgoLQ0NAAAEhJSUF7ezt0Oh3kcjmys7NRXl4OjuMwevRoaDQanD59GgAwZswY6HQ6nD9/HjKZDGPHjkVFRQUsFgu0Wi3CwsJw6tQpAEBiYiK6urrsb4Lc3FycOHECZrMZGo0GkZGRqKmpAQDEx8fDaDTi7NmzAICcnBxUV1fDZDIhJCQEMTExqK6uBgDExcXBbDajtbUVer0eZrMZdXV16OnpQVBQEBISElBZWWnvDQDNzc0AgIyMDDQ0NNh7jxkzBidOnAAAREVFQaFQoLGxEQCQlpaG5uZmdHZ2QqlUIi0tDcePHwcAREREQKVSOfRua2uDXq+HQqFAVlYWysrK7L2Dg4NRX18PAEhOTsa5c+d4e2u1WtTV1QHouzKyXq/n7R0eHm7fsTAhIQHd3d323mPHjkVnZyfKysoQGhqKqKgoh94mkwltbW0AgOzsbNTW1sJoNCIkJASxsbH2ZTY2NhZWq9VhmT19+rS9d2JiosMyK5PJ7L3T09PR2NiIrq4uqFQqJCcnM3u3tLTAYDAM2lutVg+6zPbvHRYWhtDQUIdltqOjAx0dHQgICEBOTg6OHz8Oq9WKnp4eGAwGh94GgwHnzp0bsMwO1runp2fQZTY0NBTR0dGDLrMAkJWVhVOnTsFoNCI4OHjE6wjbH2buXkfY1j2+vI6w9fbEOsJiseDUqVM+vY6orKxEb2+v29cRtnWP2NcRzlywQjSXtzCZTAgODsb27dtRVFRkv33RokU4f/48PvroI+bjU1JSsHTpUixdunREz2m7vEVaWhoCAn7aMKZSqXj3K3IGx3GSvE6KO1AbNurDRn3YqA8b9eHnK2188vIWSqUSEydOxN69e+23Wa1W7N27F1OmTPH6c2o0Gmi1Wvs/d0x8ANj/QiEDURs26sNGfdioDxv14SfFNqL62GvZsmVYtGgRJk2ahEsuuQQvvfQSOjs7UVxcDABYuHAhEhISsHr1agB9W3Zsm75MJhMaGhrw/fffIzQ0FBkZGcN6Tm/r7e0V5Of6AmrDRn3YqA8b9WGjPvyk2EZUk5958+ahtbUVK1asQFNTEyZMmICdO3faP8etq6tz+CjqzJkzuPjii+1fP//883j++ecxc+ZM7Nu3b1jP6W2hoaGC/FxfQG3YqA8b9WGjPmzUh58U24hmnx+xsO3zM378eMjlcrc/f3d3N4KCgtz+vFJAbdioDxv1YaM+bNSHn6+08cl9fvyFbe97MhC1YaM+bNSHjfqwUR9+UmxDkx9CCCGE+BWa/HhZfHy80EMQLWrDRn3YqA8b9WGjPvyk2IYmP15mMpmEHoJoURs26sNGfdioDxv14SfFNjT58TLbGTbJQL7WhuM4dPSY0aQ3oqPH7NTZRV3ha328jfqwUR826sNPim1Edag7Ib7AYDRj98l2lJS2olH/019EcRolivKiMCszHKEqemsRQohY0aHu/Xj6UHeLxeKR55UCX2hzqF6HJ/fUwGi2AgAufPPYTv6uUgRgRUEqJiVq3fqzfaGPkKgPG/Vhoz78fKUNHeouYraL0pGBxN7mUL0Oy3dVwWi2goPjxAf/+5oDYDRbsXxXFQ7V69z688XeR2jUh436sFEfflJsQ5MfLzMajUIPQbTE3MZgNOPJPTXguIGTnv44ABwHPLmnBgaj2W1jEHMfMaA+bNSHjfrwk2Ibmvx4WUhIiNBDEC0xt9l9st2+xWc4bFuAdp9sd9sYxNxHDKgPG/Vhoz78pNiGJj9eFhsbK/QQREusbTiOQ0lp67AnPhcqKW1121FgYu0jFtSHjfqwCdnH20eOOkuKyw5NfrysqqpK6CGIlljb6IwWh6O6hosD0Kg3QW+0uGUcYu0jFtSHjfqwCdHHYDTjwx9bsPi9Mlz31jEs3Nb3v4vfK8OHP7a49WPzkZDiskPH4xIyhO7ekU1eunot0KrprUYI+Un/I0cv1KQ3Yf3BBmw61OiRI0cJbfnxOiluPnQXsbYJChzZIZ7BI3y8jVj7iAX1YaM+bN7sI/SRo86S4rJDkx8vs1oHzvJJH7G20arkiNMo7efxGS4Z+k58qFG5Z/Ij1j5iQX3YqA+bt/qI4chRZ0lx2aHJj5e1tLQIPQTREmsbmUyGorwolx5blBcFmczZadPgxNpHLKgPG/Vh81YfMRw56iwpLjs0+SFkGGZlhkOlCBj21h8Z+s70PCsz3JPDIoT4ELEcOUpo8uN1mZmZQg9BtMTcJlSlwIqCVMhkGHICJAMgkwErClLdeo0vMfcRA+rDRn3YvNFHLEeOOkuKyw5Nfrzs9OnTQg9BtMTeZlKiFqsK0+1bgPpPgmy3qRQBWFWY7vYjNMTeR2jUh436sHmjjzuOHBWCFJcdOv7Wy3p6eoQegmj5QptJiVq8syBv0Ku6x/7vqu6zsyIQonT/RQB9oY+QqA8b9WHzRh+xHDnqLCkuOzT58bKgoCChhyBavtImVKXA3PxoFOVFQW+0oKvXguBAOTQqudt2bh6Mr/QRCvVhoz5s3uhjO3K0SW9yar8fGfr+uHLXkaPOkuKyQx97eVliYqLQQxAtX2sjk8mgVSsQq1FBq1Z4dOID+F4fb6M+bNSHzRt9xHLkqLOkuOzQ5MfLTp48KfQQRIvasFEfNurDRn3YvNXHF48cleKyQ5MfQgghxEvEcOQoocmP10VHRws9BNGiNmzUh436sFEfNm/2EfrIUWdJcdmhqaSXCfWZrS+gNmzUh436sFEfNm/3EfLIUWdJcdmhyY+XNTc3IyIiQuhhiBK1YaM+bNSHjfqwCdFHqCNHnSXFZYcmP4QQQoiAbEeOatX0K9lbaJ8fL0tPTxd6CKJFbdioDxv1YaM+bNSHnxTb0OTHyxobG4UegmhRGzbqw0Z92KgPG/XhJ8U2NPnxIqPRiLVr18JoNAo9FNGhNmzUh436sFEfNurDT6ptaPLjRUajEZs2bZLcQuQO1IaN+rBRHzbqw0Z9+Em1DU1+CCGEEOJXaPJDCCGEEL9Cx9X1w3F919q1WCxuf26r1Yrg4GBYrVaPPL8vozZs1IeN+rBRHzbqw8+X2tjGZ/s9ziLjhnMvP9LT04OvvvpK6GEQQgghxAVTp06FWq1m3ocmP/1YrVaYTCbI5eI6wyYhhBBC+HEcB4vFAqVSiYAA9l49NPkhhBBCiF+hHZ4JIYQQ4ldo8kMIIYQQv0KTHzdbt24dUlJSoFarMXnyZHz77bfM+7///vvIycmBWq3GRRddhE8//dRLI/U+Z9ps3rwZMpnM4d9QO7D5sv/85z+45pprEB8fD5lMhpKSkiEfs2/fPvzsZz+DSqVCRkYGNm/e7PFxCsXZPvv27Ruw/MhkMjQ1NXlnwF60evVq/PznP4dGo0F0dDSKiopQUVEx5OP8Zd3jSh9/Wf+8+uqrGDduHLRaLbRaLaZMmYJ//etfzMdIZbmhyY8bbdu2DcuWLcPKlStx5MgRjB8/HoWFhWhpaRn0/l9//TUWLFiAP/zhD/juu+9QVFSEoqIi/Pjjj14euec52wYAtFotGhsb7f9OnTrlxRF7V2dnJ8aPH49169YN6/41NTW46qqr8Mtf/hLff/89li5diltuuQW7du3y8EiF4Wwfm4qKCodlKDo62kMjFM7+/ftx11134eDBg9i9ezd6e3sxe/ZsdHZ28j7Gn9Y9rvQB/GP9k5iYiGeeeQaHDx/GoUOHcPnll2POnDkoLS0d9P6SWm444jaXXHIJd9ddd9m/tlgsXHx8PLd69epB73/99ddzV111lcNtkydP5m6//XaPjlMIzrbZtGkTN2rUKC+NTlwAcB9++CHzPg8++CCXl5fncNu8efO4wsJCD45MHIbT5/PPP+cAcOfOnfPKmMSkpaWFA8Dt37+f9z7+tO7pbzh9/Hn9ExYWxr3++uuDfk9Kyw1t+XETk8mEw4cPo6CgwH5bQEAACgoKcODAgUEfc+DAAYf7A0BhYSHv/X2VK20AwGAwIDk5GUlJScy/RvyRvyw7IzVhwgTExcVh1qxZfnP+ro6ODgBAeHg47338efkZTh/A/9Y/FosF7777Ljo7OzFlypRB7yOl5YYmP27S1tYGi8WCmJgYh9tjYmJ49zNoampy6v6+ypU22dnZ2LhxIz766CO89dZbsFqtuPTSS1FfX++NIYse37Kj0+nQ3d0t0KjEIy4uDuvXr8cHH3yADz74AElJSbjssstw5MgRoYfmUVarFUuXLsXUqVORn5/Pez9/Wff0N9w+/rT+OXbsGEJDQ6FSqXDHHXfgww8/RG5u7qD3ldJyQ5e3IKI0ZcoUh78+Lr30UowdOxavvfYannrqKQFHRnxBdnY2srOz7V9feumlqKqqwosvvogtW7YIODLPuuuuu/Djjz/iyy+/FHooojTcPv60/snOzsb333+Pjo4ObN++HYsWLcL+/ft5J0BSQVt+3CQyMhJyuRzNzc0Otzc3NyM2NnbQx8TGxjp1f1/lSpv+AgMDcfHFF6OystITQ/Q5fMuOVqtFUFCQQKMSt0suuUTSy8/dd9+Nf/zjH/j888+RmJjIvK+/rHsu5Eyf/qS8/lEqlcjIyMDEiROxevVqjB8/Hi+//PKg95XSckOTHzdRKpWYOHEi9u7da7/NarVi7969vJ+fTpkyxeH+ALB7927e+/sqV9r0Z7FYcOzYMcTFxXlqmD7FX5Ydd/r+++8lufxwHIe7774bH374If79738jNTV1yMf40/LjSp/+/Gn9Y7VaYTQaB/2epJYbofe4lpJ3332XU6lU3ObNm7mysjLutttu40aPHs01NTVxHMdxN910E/fQQw/Z7//VV19xCoWCe/7557ny8nJu5cqVXGBgIHfs2DGhXoLHONvmiSee4Hbt2sVVVVVxhw8f5ubPn8+p1WqutLRUqJfgUXq9nvvuu++47777jgPAvfDCC9x3333HnTp1iuM4jnvooYe4m266yX7/6upqLjg4mHvggQe48vJybt26dZxcLud27twp1EvwKGf7vPjii1xJSQl38uRJ7tixY9y9997LBQQEcHv27BHqJXjMkiVLuFGjRnH79u3jGhsb7f+6urrs9/HndY8rffxl/fPQQw9x+/fv52pqargffviBe+ihhziZTMZ99tlnHMdJe7mhyY+brV27lhszZgynVCq5Sy65hDt48KD9ezNnzuQWLVrkcP/33nuPy8rK4pRKJZeXl8f985//9PKIvceZNkuXLrXfNyYmhvv1r3/NHTlyRIBRe4ft0Oz+/2xNFi1axM2cOXPAYyZMmMAplUouLS2N27Rpk9fH7S3O9nn22We59PR0Tq1Wc+Hh4dxll13G/fvf/xZm8B42WBcADsuDP697XOnjL+ufm2++mUtOTuaUSiUXFRXF/epXv7JPfDhO2ssNXdiUEEIIIX6F9vkhhBBCiF+hyQ8hhBBC/ApNfgghhBDiV2jyQwghhBC/QpMfQgghhPgVmvwQQgghxK/Q5IcQQgghfoUmP4QQQgjxKzT5IYRIzuOPPw6ZTCb0MFxWW1sLmUyGzZs322/z9ddEiJjQ5IcQH7Z582bIZDL7P7VajaysLNx9990Drr5cW1uL4uJipKenQ61WIzY2FjNmzMDKlSsHfe5PPvkE11xzDWJiYqBUKhEeHo4ZM2ZgzZo10Ol03nh5hBDiEQqhB0AIGbknn3wSqamp6OnpwZdffolXX30Vn376KX788UcEBwejsrISP//5zxEUFISbb74ZKSkpaGxsxJEjR/Dss8/iiSeesD+X1WrFH/7wB2zevBkXXXQR7rzzTiQlJUGv1+PAgQNYvnw5Pv300wFXdyaetXz5cjz00ENCD4MQSaDJDyEScOWVV2LSpEkAgFtuuQURERF44YUX8NFHH2HBggV48cUXYTAY8P333yM5OdnhsS0tLQ5fP/fcc9i8eTPuu+8+rFmzxuGjlnvvvReNjY148803Pf+ifJDZbIbVaoVSqXT7cysUCigUtMomxB3oYy9CJOjyyy8HANTU1AAAqqqqkJiYOGDiAwDR0dH2/7+rqwvPPvss8vLy8Ne//nXQfUzi4uLw5z//ecgxfPHFF7juuuswZswYqFQqJCUl4b777kN3d7fD/RYvXozQ0FA0NDSgqKgIoaGhiIqKwv333w+LxeJw37Nnz+Kmm26CVqvF6NGjsWjRIhw9enTA/jF83nrrLUycOBFBQUEIDw/H/Pnzcfr06SEfNxjbfjnPP/88XnrpJaSnp0OlUqGsrAwmkwkrVqzAxIkTMWrUKISEhGD69On4/PPPBzzP+fPnsXjxYowaNcr+ms6fPz/gfoPt87Np0yZcfvnliI6OhkqlQm5uLl599dUBj01JScHVV1+NL7/8EpdccgnUajXS0tIGTGJ7e3vxxBNPIDMzE2q1GhEREZg2bRp2797tUiNCxIr+jCBEgqqqqgAAERERAIDk5GTs2bMH//73v+0To8F8+eWXOH/+PO6//37I5fIRjeH9999HV1cXlixZgoiICHz77bdYu3Yt6uvr8f777zvc12KxoLCwEJMnT8bzzz+PPXv2YM2aNUhPT8eSJUsA9H0cd8011+Dbb7/FkiVLkJOTg48++giLFi0a1niefvppPPbYY7j++utxyy23oLW1FWvXrsWMGTPw3XffYfTo0S69zk2bNqGnpwe33XYbVCoVwsPDodPp8Prrr2PBggW49dZbodfr8X//938oLCzEt99+iwkTJgAAOI7DnDlz8OWXX+KOO+7A2LFj8eGHHw77Nb366qvIy8vDb37zGygUCnzyySe48847YbVacddddznct7KyEr/73e/whz/8AYsWLcLGjRuxePFiTJw4EXl5eQD6JlirV6/GLbfcgksuuQQ6nQ6HDh3CkSNHMGvWLJf6ECJKHCHEZ23atIkDwO3Zs4drbW3lTp8+zb377rtcREQEFxQUxNXX13Mcx3E//vgjFxQUxAHgJkyYwN17771cSUkJ19nZ6fB8L7/8MgeAKykpcbjdbDZzra2tDv+sVitzbF1dXQNuW716NSeTybhTp07Zb1u0aBEHgHvyyScd7nvxxRdzEydOtH/9wQcfcAC4l156yX6bxWLhLr/8cg4At2nTJvvtK1eu5C5cvdXW1nJyuZx7+umnHX7GsWPHOIVCMeD24aipqeEAcFqtlmtpaXH4ntls5oxGo8Nt586d42JiYribb77ZfltJSQkHgHvuueccHjt9+vQhXxPHDd64sLCQS0tLc7gtOTmZA8D95z//sd/W0tLCqVQq7k9/+pP9tvHjx3NXXXXVMF49Ib6NPvYiRAIKCgoQFRWFpKQkzJ8/H6Ghofjwww+RkJAAAMjLy8P333+PG2+8EbW1tXj55ZdRVFSEmJgYbNiwwf48tqO4QkNDHZ7/2LFjiIqKcvh39uxZ5piCgoLs/39nZyfa2tpw6aWXguM4fPfddwPuf8cddzh8PX36dFRXV9u/3rlzJwIDA3HrrbfabwsICBiwhWMwO3bsgNVqxfXXX4+2tjb7v9jYWGRmZg76cdRw/fa3v0VUVJTDbXK53L7fj9VqRXt7O8xmMyZNmoQjR47Y7/fpp59CoVDYt27ZHnvPPfcM62df2LijowNtbW2YOXMmqqur0dHR4XDf3NxcTJ8+3f51VFQUsrOzHRqPHj0apaWlOHny5LB+PiG+ij72IkQC1q1bh6ysLCgUCsTExCA7OxsBAY5/22RlZWHLli2wWCwoKyvDP/7xDzz33HO47bbbkJqaioKCAmg0GgCAwWBweGxGRoZ9v48333wTW7ZsGXJMdXV1WLFiBT7++GOcO3fO4Xv9fzGr1eoBE4iwsDCHx506dQpxcXEIDg4eMLahnDx5EhzHITMzc9DvBwYGDvkcfFJTUwe9/Y033sCaNWtw/Phx9Pb2Dnp/22vqP9nMzs4e1s/+6quvsHLlShw4cABdXV0O3+vo6MCoUaPsX48ZM2bA4/s3fvLJJzFnzhxkZWUhPz8fV1xxBW666SaMGzduWOMhxFfQ5IcQCbjkkkvsR3sNRS6X46KLLsJFF12EKVOm4Je//CXefvttFBQUICcnBwDw448/Ys6cOfbHhIaGoqCgAEDffkFDsVgsmDVrFtrb2/HnP/8ZOTk5CAkJQUNDAxYvXgyr1TpgTJ5ktVohk8nwr3/9a9Cf1X/y4YwLt77YvPXWW1i8eDGKiorwwAMPIDo6GnK5HKtXr7bvjzVSVVVV+NWvfoWcnBy88MILSEpKglKpxKeffooXX3xx2I05jrP//zNmzEBVVRU++ugjfPbZZ3j99dfx4osvYv369bjlllvcMm5CxIAmP4T4MduEqbGxEUDfR02jRo3Cu+++i4cffnjA1qPhOnbsGE6cOIE33ngDCxcutN8+kqOGkpOT8fnnn6Orq8th609lZeWQj01PTwfHcUhNTUVWVpbLYxiu7du3Iy0tDTt27HA4Qqv/CSWTk5Oxd+9eGAwGhwlYRUXFkD/jk08+gdFoxMcff+ywVWckH+EBQHh4OIqLi1FcXAyDwYAZM2bg8ccfp8kPkRTa54cQP/DFF184fPRi8+mnnwL46WOW4OBgPPjgg/jxxx/x0EMPOWwVsBnstv5sWxkuvC/HcXj55ZddGj8AFBYWore312EfJavVinXr1g352GuvvRZyuRxPPPHEgPFzHDfk/kvOGuz1f/PNNzhw4IDD/X7961/DbDY7HJ5usViwdu1al35GR0cHNm3a5PK4+3cIDQ1FRkYGjEajy89JiBjRlh9C/MCzzz6Lw4cP49prr7Xvv3HkyBG8+eabCA8Px9KlS+33feihh1BeXo6//vWv+Oyzz/Db3/4WiYmJOHfuHI4cOYL3338f0dHRUKvVvD8vJycH6enpuP/++9HQ0ACtVosPPvhgwL4/zigqKsIll1yCP/3pT6isrEROTg4+/vhjtLe3AwDzulfp6elYtWoVHn74YdTW1qKoqAgajQY1NTX48MMPcdttt+H+++8HAOzbtw+//OUvsXLlSjz++OMujfXqq6/Gjh07MHfuXFx11VWoqanB+vXrkZub67A/1TXXXIOpU6fioYceQm1tLXJzc7Fjx44B+0QNZvbs2VAqlbjmmmtw++23w2AwYMOGDYiOjrZvyXNWbm4uLrvsMkycOBHh4eE4dOgQtm/fjrvvvtul5yNErGjyQ4gfeOSRR/DOO+9g//79ePvtt9HV1YW4uDjMnz8fjz32mMNOuAEBAdiyZQt++9vfYsOGDVi7di3OnTuH0NBQ5Ofn4+mnn8att97K3E8mMDAQn3zyCf74xz9i9erVUKvVmDt3Lu6++26MHz/epdcgl8vxz3/+E/feey/eeOMNBAQEYO7cuVi5ciWmTp3KnIwBfZO6rKwsvPjii/bLeSQlJWH27Nn4zW9+Y7+fbXISFxfn0jiBvhM3NjU14bXXXsOuXbuQm5uLt956C++//z727dtnv19AQAA+/vhjLF26FG+99RZkMhl+85vfYM2aNbj44ouZPyM7Oxvbt2/H8uXLcf/99yM2NhZLlixBVFQUbr75ZpfG/cc//hEff/wxPvvsMxiNRiQnJ2PVqlV44IEHXHo+QsRKxg1nGzYhhIhUSUkJ5s6diy+//BJTp04d8fM9+OCD2Lp1KyorK6FSqdwwQkKI2NA+P4QQn9H/0hi2/WO0Wi1+9rOfueVnfP7553jsscdo4kOIhNHHXoQQn3HPPfegu7sbU6ZMgdFoxI4dO/D111/jL3/5y6CHnLviv//9r1uehxAiXvSxFyHEZ7zzzjtYs2YNKisr0dPTg4yMDCxZsoR2yCWEOIUmP4QQQgjxK7TPDyGEEEL8Ck1+CCGEEOJXaPJDCCGEEL9Ckx9CCCGE+BWa/BBCCCHEr9DkhxBCCCF+hSY/hBBCCPErNPkhhBBCiF+hyQ8hhBBC/Mr/B2UolDUOWx93AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mueller Matrix Original\n",
      "[[0.5        0.29399553 0.40443371 0.        ]\n",
      " [0.28578716 0.16804029 0.23116392 0.        ]\n",
      " [0.40673081 0.23915408 0.3289913  0.        ]\n",
      " [0.05381217 0.03164107 0.04352691 0.        ]]\n",
      "------------------------------\n",
      "Mueller Matrix Measured\n",
      "[[ 5.00000000e-01  2.93995530e-01  4.04433713e-01 -1.38777878e-16]\n",
      " [ 2.85787156e-01  1.68040293e-01  2.31163922e-01 -5.55111512e-17]\n",
      " [ 4.06730810e-01  2.39154080e-01  3.28991304e-01 -2.01227923e-16]\n",
      " [ 5.38121679e-02  3.16410736e-02  4.35269098e-02 -1.11022302e-16]]\n",
      "------------------------------\n",
      "Percent Difference\n",
      "[[-8.88178420e-14  3.77632621e-14  2.47061679e-13             inf]\n",
      " [-1.16543693e-13 -8.25860725e-14 -2.88165130e-13             inf]\n",
      " [-4.09443911e-14 -8.12401065e-14 -1.68731363e-13             inf]\n",
      " [ 4.12628990e-13 -1.18422110e-12 -7.65197688e-13             inf]]\n",
      "------------------------------\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Work\\AppData\\Local\\Temp\\ipykernel_17996\\3633839214.py:40: RuntimeWarning: divide by zero encountered in divide\n",
      "  print(100*(M_to_measure-M_measured)/M_to_measure)\n"
     ]
    }
   ],
   "source": [
    "# katsu has a simulator built-in, so let's try use that\n",
    "thetas = np.linspace(0, np.pi, 24) # 24 measurements between 0 and pi\n",
    "\n",
    "# set up polarization state generator\n",
    "psg_polarizer = linear_polarizer(0, shape=thetas.shape)\n",
    "psg_retarder  = linear_retarder(thetas, np.pi / 2, shape=thetas.shape)\n",
    "PSG = psg_retarder @ psg_polarizer\n",
    "\n",
    "# set up polarization state generator\n",
    "psa_polarizer = linear_polarizer(0, shape=thetas.shape)\n",
    "psa_retarder  = linear_retarder(PSA_ANGULAR_INCREMENT * thetas, np.pi / 2, shape=thetas.shape)\n",
    "PSA = psa_polarizer @ psa_retarder\n",
    "\n",
    "# set up system Mueller matrix\n",
    "Msys = PSA @ M_to_measure @ PSG\n",
    "\n",
    "# propagate Stokes vector to get power\n",
    "Sin = stokes_from_parameters(1, 0, 0, 0)\n",
    "Sout = Msys @ Sin\n",
    "power_measured = Sout[..., 0, 0]\n",
    "\n",
    "plt.figure()\n",
    "plt.ylabel('Power observed')\n",
    "plt.xlabel('PSG angle, radians')\n",
    "plt.plot(thetas, power_measured, marker='o', linestyle='None', markersize=10)\n",
    "plt.show()\n",
    "\n",
    "# perform Mueller matrix measurement\n",
    "M_measured = full_mueller_polarimetry(thetas, power_measured, PSA_ANGULAR_INCREMENT,\n",
    "                                      starting_angles=starting_angles,\n",
    "                                      starting_polarization=starting_polarization)\n",
    "\n",
    "print('Mueller Matrix Original')\n",
    "print(M_to_measure)\n",
    "print('-'*30)\n",
    "print('Mueller Matrix Measured')\n",
    "print(M_measured)\n",
    "print('-'*30)\n",
    "print('Percent Difference')\n",
    "print(100*(M_to_measure-M_measured)/M_to_measure)\n",
    "print('-'*30)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed it can! Experimentally I've found that 24 is a good minimum number of measurements to get accurate measurements of the Mueller matrix in simulation. In the laboratory, the presence of noise further complicates things, so more measurements are needed for more accurate polarimetry. But this is the subject of future investigations!"
   ]
  }
 ],
 "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.10.14"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}