{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 2: Classic Zeno effect\n",
    "\n",
    "The classic Zeno effect analyzes the intensity transmitted by a certain number of perfect polarizers with equally spaced azimuths between 0º and 90º. Using only two polarizers, as their transmission axes are perpendicular to each other, the output intensity of the system will be zero. As the number of polarizers increases, the output intensity will also increase. In the limit of infinite polarizers, if the incident light is linearly polarized with the same azimuth than the first polarizer, the output wave will have the same intensity than the incident light. \n",
    "\n",
    "This example will show the classic Zeno effect both in Jones and Mueller formalisms. We will also analyze the influence of using imperfect polarizers.\n",
    "\n",
    "## Jones formalism\n",
    "\n",
    "Let us start by Jones formalism. Jones formalism is easier to visualize and understand.\n",
    "\n",
    "### Perfect polarizers\n",
    "\n",
    "This is the easiest case. Start by importing all the required modules and classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from py_pol import degrees\n",
    "from py_pol.jones_vector import Jones_vector\n",
    "from py_pol.jones_matrix import Jones_matrix\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will work using linear polarization with the same azimuth as the first polarizer (0º)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Incident wave = \n",
      "[+1.000]   \n",
      "[+0.000]   \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Incident wave\n",
    "E0 = Jones_vector('Incident wave')\n",
    "E0.linear_light(intensity=1, azimuth=0*degrees)\n",
    "print(E0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will construct two arrays, one with the number of polarizers and other with the output intensity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nmax = 30\n",
    "Narray = np.arange(2,Nmax+1)\n",
    "J = Jones_matrix('System')\n",
    "I = np.zeros(Nmax-1)\n",
    "\n",
    "# For loop to calculate the intensities\n",
    "for ind, N in enumerate(Narray):\n",
    "    # Create the system of polarizers\n",
    "    angles = np.linspace(90*degrees, 0*degrees, N) # Start by 90º and end in 0º because the rightest  element is the first to be crossed by the light wave\n",
    "    J.diattenuator_perfect(azimuth=angles)\n",
    "    # Multiply them\n",
    "    J.prod()\n",
    "    # Multiply the system matrix by the incident light wave to obtain the output\n",
    "    E_final = J * E0\n",
    "    # Store the intensity\n",
    "    I[ind] = E_final.parameters.intensity(out_number=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, plot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEpCAYAAABr364UAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApqUlEQVR4nO3deZhcVZ3/8fcnAYRmVQmIgXSjooIKCJHFBQLjAqgDjjIK/QNhnGkRQXGZEY2OQY2KIzPuMq0TUekZ3ADRQRRHEGcQJSD7ogmkk8gWZadlSfL9/XFOkZtKVXff7ttdVd2f1/Pcp6rOXc65dbvvt+4595yriMDMzKyMGa0ugJmZdR4HDzMzK83Bw8zMSnPwMDOz0hw8zMysNAcPMzMrzcHD2oKSkyRdI2lIUuSpp9Vl6zSt/O4kXZrzPm6y87bJ5eAxjUg6q3BiKU4P5pP2v0jasUXF+xDwRWAPQMDdeVoz2QWRtI2kBZIWTHbeZp3CwWN6eoJ1J+d7gC1IJ+33A9dLenkLyvTu/PpeoCsinpGnFS0oyzbAR/Nk5SwHbgUeaHVBbGJt1OoCWEtcHhHzah8kdQFvBL5AOnF+T9KzIuIvk1EYSdsBs/LHr4WHPehYEXFsq8tgk8NXHkZEDEXEt4F35aRnAEdMYhE2K5Tl4UnM18zGyMHDir4LrM3v9y7OkDRD0jGSLpa0StLjku6Q9B1J+zbaWG43iNzWMiM3iP9W0v05/RRJASwrrFNsi1kw3jIU1t1c0vslXS7pXkmPSrpN0gWSeiVtnJe7FLi9SXk2KNMw+fXU1smfXybpx7ncQ7mN6SRJw/4PSvobSRfl9R6TtFLSgKS9RlOOum3NlHSQpM9LukrS3YXv8DxJBw+z7pMN4blN6HRJt+R9ub/Rco2+i1FMPU2+yy9KujXn91Au/wckbd6kvE9uT9Kukr4paYWkJySdX1huu9zWd4OkR/LfxYr8d/IxSd1lv+dpIyI8TZMJOAsI4NJhlrk7L9NfSNsSuDinBynAPFD4vAY4qcG2FuT53wTOz+9XA/fl9y8F7gJWFbZ1V2F6/3jLkNfdjRQQass+UShDberJy547THnWK9MI33VPYRtvzHlGzveJwrzzgI0arD8jf2+15YrfW21/39Ek7/X2qZD+wrp9fhR4uC7tQ022eWme/4/A0sL6DwL3N1juuELaTg2+x+J0/zBl/hvgL4X5Q8Bjhc/XAdsP8x0cAzyS3z+Yt3V+XqYbuKPuO743/23V0k5o9f9tu04tL4CnSTzYIwQPUvVR7R/nM4X083LatcBhwGY5fRvgg/mfeQ3wsrrtLcjrPZRPNO8gNYYDbAdsld/31P5Zhyn7WMvwNFIjbgC3AYcDm+R5WwIvBxYBOxbWGbE8o/iun9xGPjn+BNg5z9ucdBJeQ5MTNnAq64Lkh4Etc/ps0hViLYAc0GDdZifi5+Z1XwdsD6hwLD6cT55rgX0bbPPSwrFcDhwCzMjzntNgueNG+T1tDFyW17kG2LQw7yXA47lcnwbmkO7EmwnsC1yR1/vpMN/BQ7lML8zpAp6d3y/Ky/wBeEVhf55CCrQfB45o9f9tu04tL4CnSTzYIwePkwr/dG/Maa/Mn28HntZkvX/Ky/y4Ln1BYXt9w5Rr2JP1OMvwmZy+Cpg9yu9p2PKU3QZwA/CUBsvUvp8HgM0L6Zuz7qrqUw3Wmwn8Ks+/rMH8hsFjFGX+SF7vGw3m1YLC47UTcZNt1JY7bpR5npmXvwforpv3v3nee5qs+1Tgj3mZuU2+g6XkHxoN1r8pL/PmsR7n6Ty5zWOaU9Ij6f2kEy3AIPCj/P6t+fWsiLi3yWb+M78eJGlmg/l/Jv3KG6vxlOGY/PrZiPjjOMowHmdExGMN0v+VdEW2FfCqQvqrc9rjrDsmT4qINaRfxQCvkPSMispZO+YvG2aZn0TEDVVkJukdwNtJ1XhviojBwrxn53L8hRRgNhAR95Gu6GD976/oS9H8rsEH8+sOJYtu+Fbd6erAWkNuA3eSLtUfz59fml/fk//Zh9MFPJ30K7JocUSsHltRx16G3PhaO7FeOI78x+vSRokR8aCk3wH7A3uR2oXI7wGuzSfIRi4jVedslJcf1f5J2gw4gVR9txvp13v9eeCZw2zi16PJZxTlOAD4fP54ckRcVrdI7ZhvAtwuqdmmtsivOzWZP1x5LyRVf50uaRfg+8AVwwQbK3DwmJ6eIDUMQrpsf4TUHnAx8PW6E1btV9nWeRpJV4O0VWMs53jLsH0hbfk4yzAew13x1ObNKqTNqpu3gYh4VNKfSfs4q9lyRZJ2IAWy5xaSHyE1xK8lVYdtS6o2a2a8xxJJc0gn6o2Br0bEvzdYrHbMZ7L+cWym0d8dDF/e00l3Ff41cGKeVku6ktTG9rWIuH8UeU9LDh7T03qdBEdQq9o8PCIuGGN+4x1iZKxlaPpztY0MV8anVJzX50iB4zZSg/0lxR8KuapoyQjbGNexVOqQ+kNSwPsl60YWqFc75r+LiNK3JRc0LW+uSjxc0n7AG4ADSMFk/zz9o6RXRcS148h/ynKbh43k7vy6WweW4a7C++6KyjIWw1UD1X5hF38h1943LbOkTUnVc/XrNlt+E1JVFUBvRJzboEpsNL/wx+sbwJ6kvj1viognmixXO+a7SJrQH7kRcUVEfCAi9idV4x1FulKdBXx9IvPuZA4eNpJanfEbO60MEbGMdQHksBKr1jpKomEq20s4sFGipC1Z175xdWFW7f0ukmY32eYBrKs5uLrJMkXbsu5K5ndNlnnlKLYzZpLmA39Lqio7PCL+NMzitWO+BekGgkkREY9ExDlAX07au1lHxOnOwcNGclZ+nStp2HGLJD21Dcvw7fz6vmFOxPUeLLzfZpTrDOd9+Zd/vVOATXN+Pyuk/yynbUyqXlpPvpvsI/njryLirvplGniQ1L4F8KIG29wBOHkU2xkTSa8HPpbLcGxEXDfc8hFxC6kfB6QG7aYncEmbSSpdxdfkmNTUGs1FarS3Og4eNqyIuIjU6xpgkaTT8okGSCdrSYdL+iHp1tN2K8PppIbnbYFfSfrr2klD0haS5kk6R4Wh6HMj6R354/EV7MIc4Lza0BuSuiS9l3Wj9p4eEUOF/B8BPpk/vkvSfElb5HVnA/9F6txY60A4okhjhtVOxosk7Zm3N0PSX5HaHyakjUjSc4GzSeeb0yLi3BFWqTmZ1PnzhaRj98paFVYu9wskfZjUl2Mst9veIOmTkl5S+JuQpH1IjwcAuHKYO96mt1Z3NPE0eROjGJ6kyXqbs66Hd7HX9AN1ad+oW29BTj9rhO331LZRdRnyui8CVhSWeZwmw5MU1jmtMO9hUh39MuCUUX5nPYX1hxue5HwaD08ykw2HJykOnbEGOLFJ3s32aV/S8B7F/ap9/jOpTaThcWCUnf8aLQccV8jzboYfqmSnuu0dyvrDlzwG/Ckfw+Lx6x7Nd1C3THG7q/N3UNzuKmD3Vv/ftuvku61sRJF+Cb9B0muBvyOdhGaRTmRLgN+SrgwmrC/FeMoQEddLegHpl+wRpDuONiXddXQd8B1gZd1qHyPVzfcCz2Fd4/U2Yyj7DyQdRBpyZD/SiepGUmPsVyJibYN11gBvlXQBqf59b1LHwTtJVwlnRMRVJcvxG0n7k4L6gaSAfCdwEbCQFLAm2nYjzF+vDBHxk3zlcjKp3eo5pGNwP+m5IRcB34tCB8MSDgdeQ2o/mkO6YeBx4GbS39G/RUR9nyXLamPbmFmFchXV7QAR0Qm3DJuV4jYPMzMrzcHDzMxKc/AwM7PSHDzMzKy0adNgvu2220ZPT0+ri2Fm1lGuuuqqP0XEBoNvTptbdXt6eli8eHGri2Fm1lEkNbwN2tVWZmZWmoOHmZmV5uBhZmalOXiYmVlpDh5mZlaag4eZ2RQ0MDBAT08PM2bMoKenh4GBgUq37+BhZtZiZU70o1l2YGCAvr4+BgcHiQgGBwfp6+urNoC0ekz4yZr23nvvMDObLGeffXZ0d3eHpOju7o6zzz676XJdXV3FZ5NEV1dXw+VHu2x3d/d6y9Sm7u7u0vsBLI4G59SWn9Qna3LwMLPxmoiAUOZEP9plJTVcTlLpfXbwcPAwswZaHRDKnOhHu6yvPBw8zGwMOikgTMSVR5n9GomDh4OHWUebqgFhIto8ynxfI3HwcPAwazsOCOW+h7LLVsHBw8HDbNKM5gTngNAZHDwcPMzGpeqrBAeEzuDg4eBhtoFWVhs5IHQGBw8HD7P1tLrayAGhMzh4OHjYNDHaE2erq40cEDqDg4eDh3WwiaheaodqIweE9ufg4eBhHarVQ12ULUNteQeFqcHBw8HD2kyrq5ccEGw0HDwcPKyNtEP1Uq0cDgg2nGbBQ2ne1Dd37txYvHhxq4thBkBPTw+Dg4MbpHd3d7Ns2bIxL1t7jsPQ0NCTaV1dXfT399Pb21tJ2W16kXRVRMytT/fDoMwqNNqH+ixfvnzU6QsXLqSrq2u9tK6uLhYuXLjBsr29vfT399Pd3Y0kuru7HThsYjS6HJmKk6utbKJNVMN2bduuXrJWoEm1la88zEYw2quJ+fPnr1ddBDA0NMT8+fM3WLbM1QSkK4ply5axdu1ali1b5isJazkHD7NhlHkWdJmqKFcvWadzg7nZMCaqYdusU7jB3KzOaKqjJqph26zTOXjYtDTa6qg5c+Y0XL9RuquibDpxtZVNS6OtYnK/CZvuOqbaStIhkm6VtETSqQ3mby3pR5KulXSjpONbUU5rT1X3s/DVhFljbXXlIWkm8HvgVcBK4ErgqIi4qbDMh4CtI+IDkmYBtwLPiIjHh9u2rzymvjJXCW7cNhudTrny2AdYEhG35WBwDnB43TIBbClJwBbAvcDqyS2mtaOJ7GdhZutrt+AxG1hR+LwypxV9CdgVuAO4Hnh3RKxttDFJfZIWS1q8atWqiSivtRH3szCbPO0WPNQgrb5e7TXANcAzgT2BL0naqtHGIqI/IuZGxNxZs2ZVWU6bRKNtxyhzZxS417bZeLRb8FgJ7FT4vCPpCqPoeODcPOzKEuB24PmTVD6bZGV6eLsqymzytFvwuBLYRdLOkjYB3gJcULfMcuCvACRtDzwPuG1SS2mTpkw7hquizCZPW91tBSDpMOBzwExgUUQslHQCQEScKemZwFnADqRqrk9HxNkjbdd3W3WmGTNm0OhvVBJr1zZs6jKzCjW722qjVhRmOBFxIXBhXdqZhfd3AK+e7HJZa8yZM6fhLbXN2jHMbHK0W7WVTROjbQR3O4ZZe3LwsElXphHc7Rhm7ant2jwmits82od7d5t1jk7pYW7TQJnOfGbWnhw8bNKV7cxnZu3HwcMq40Zws+nDwcMq4UZws+nFDeZWCTeCm01NbjC3CeVGcLPpxcHDKuFGcLPpxcHDKuFGcLPpxcHDKuFGcLPpxcHDhjXa22/BD1cym07ablRdax+1229rz9Oo3X4LODCYTXO+8rCmyjyIycymFwcPa8q335pZMw4e1pRvvzWzZhw8rCnffmtmzTh4WFO+/dbMmvHYVmZm1pTHtjIzs8o4eExDZTr+mZk14k6C04w7/plZFXzlMc2445+ZVWFcVx6S/gf4XZ6uAW6OiLUVlMsmiDv+mVkVxlttdTfwWuAUQMBjkq4nBZJaULkuIv4yznysInPmzGn4xD93/DOzMsZVbRURR0fErsDWwCuBlcBewOuAzwGXAw9KukmSW2XbgDv+mVkVKmnziIhHgCOApcB2ETEb2Ap4C3AT0A3MqiIvGx93/DOzKlTWSVDSA8AbI+LndekbAz8FfhERn6gkszFwJ0Ezs/Imo5PgEPCM+sSIeAI4Azi+wrzMzKyFqgweA8BpkmY3mLcpsEOFeZmZWQtV2Unww8AewK2SvkaqqloOPA/4DHB9hXmZmVkLVRY8IuJRSa8B3gO8F3g3EKRbeG8H+qrKy8zMWqvSHuYRsTYizsh3W70AeDXp1t3nRcS1VeZl6/N4VWY2mSoNHpIOkHSWpEuATwF/iIhrImJ1iW0cIulWSUskndpkmXmSrpF0o6RfVlX+TlUbr2pwcJCIeHK8KgcQM5solQUPSb3ApcBOwA3A64Ft87wFkg4cxTZmAl8GDgV2A46StFvdMtsAXwH+OiJeABxZ1T50Ko9XZWaTrcorjw8Dn42IvyK1e6gw7wnSECYj2QdYEhG3RcTjwDnA4XXLHA2cGxHLASLinvEWvNN5vCozm2xVBo8e4KL8vr7n4Q2kto+RzAZWFD6vzGlFzwWeKulSSVdJOrbZxiT1SVosafGqVatGkX1najYulcerMrOJUmXwWAq8qMm8+xjd8CRqkFYfiDYC9iYNyPga4COSnttoYxHRHxFzI2LurFlTd3QUj1dlZpOtyuDxBdKJvNEVxi7AaH76ryS1mdTsCNzRYJmLIuKRiPgTcBmpf8m05fGqzGyyVdnPoz/387gc+D7piuEZkg4A5gMXj2IzVwK7SNoZ+CNpYMWj65b5IfAlSRsBmwD7Av9WzV50rt7eXgcLM5s0VT+G9k3AO4F/IlVB/Si/XkkKIMOKiNWSTiL1Tp8JLIqIGyWdkOefGRE3S7oIuA5YC3w9Im6oeD/MzGwYVY6qeyIwEBEP5M+7As8E7iQ9YbCajMbIo+qamZU3GaPqfpE0jhUAEXFzRPxPRNwEbCxpswrzMjOzFhpX8JD0NEmb1D4Os+iLgQfHk5eZmbWP8V55nAw8LOlGUgP5sZJeJ2nHuuU2A9aMMy8zM2sT420wX0QaMXcPYFfSnVEnAiHpPuAa4PfAfqTH0ZqZ2RQwruARESuAbwFIOpJ0t9VSUjDZA9iT1KFvFaO428rMzDrDuIKHpG0i4n6AiCiOhfGLPJmZ2RQ03jaPVZL2AJD0QUmHNnkMrY2Rn9NhZu1ovG0ebwZqo9q+G9iOde0d19ZNN0bEE+PMb1qpPaejNtx67TkdgHuTm1lLVdZJEEDSLNa1d+zOuob0jYEnIuIplWVWUid2Euzp6WFwcHCD9O7ubpYtWzb5BTKzaadZJ8FKhyeJiFXAz/NUy3gj0oOddq8yr+nAz+kws3ZV6WNoG4mI1RFxXUScPdF5TTV+ToeZtasqH0O7h6RfSLpT0s2Svifpw7nT4E4jb8Hq+TkdZtauqrzy+BawM+n54j8CtiSNsHsBsEzSnyvMa1rwczrMrF1VOaruw8CbI+K/69K3I3UW3D0iPltJZmPQiQ3mZmatNhkN5lcCW9QnRsQ9wM/yZGZmU0CV1VYfBE6VtHWF2zQzszZUZfD4LfBr4AZJH5D0EkmbVrh9MzNrE1UGj37gBGAr4DTgN8BDkm6UNCDp/RXmZWZmLVRl8HgzsDAitgY2J/UuPx74CbA96bnmZmY2BVTZYH4f8EuAiFgDXJ8ndw40M5tiqrzy+HfgiAq3Z2ZmbarK4LEx8DpJp0na4JZdMzObOqoMHm8H5gAfAe6R9EtJn5d0vKQXS9q4wrzMzKyFKmvziIgdGgzJfiDpDqyNgSeAlg3JbmZm1fGQ7GZmVlplwUNSH7A6IhYV0yNiNXBdnszMbAqoss3jXcCs2gdJ8yTdK+khSZ+qMB8zM2uxKoNHD3B14fOngKXAAuAdkv6uwrzMzKyFqmzz+AuwFkBSN7Av8NKIuELSWtKzPRYNs76ZmXWIKq88rgVeld/3AndHxBX589XALhXm1dEGBgbo6elhxowZ9PT0MDAw0OoimZmVUuWVx0LgQklzgZcBXyzM2x54vMK8OtbAwAB9fX0MDQ0BMDg4SF9fH4CfEGhmHaOyK4+IuAQ4GLgJ+Dzw0cLsecDvq8qrk82fP//JwFEzNDTE/PnzW1QiM7Pyqu7n8WvSMz3qzQS+W2VenWr58uWl0s3M2lGlwaOZiHj7ZOTTCebMmcPg4GDDdDOzTlFZtZWkp0r6pKTvS/qqpHdKermkLUtu5xBJt0paIunUYZZ7iaQ1kt40/tJPnoULF9LV1bVeWldXFwsXLmxRiczMyqvybquzgfcAmwF7kvp5XAbcL2mppO+PtAFJM4EvA4eShjQ5StJuTZY7HfhpZaWfJL29vfT399Pd3Y0kuru76e/vd2O5mXWUKqutDgTeWRyeRNJzSIHkxYxubKt9gCURcVte/xzgcFIjfNHJwA+Al4y/2JOvt7fXwcLMOlqVweNO4I/FhIhYAiwBRrzqyGYDKwqfV5I6Gz5J0mzgDaQ7u4YNHnm8rT5wm4KZWZWqrLb6KvDGcW5DDdKi7vPngA/kR90OKyL6I2JuRMydNWvWSIubmdkoVXnlsSlwsKRPAKdHxENj2MZKYKfC5x2BO+qWmQucIwlgW+AwSasj4vwx5GdmZmNQZfA4hXQy/xDwPklXAtcUphsiYqRe5lcCu0jamVQF9hbg6OICEbFz7b2ks4AfO3CYmU2uKp8kuJ2kHUgN47VpHunxtKN6kmBErJZ0EukuqpnAooi4UdIJef6ZVZXXzMzGruoe5neSGs6fvIU2P7t8N+BFo9zGhcCFdWkNg0ZEHDfWspqZ2dhNeA/ziHiCNOLutROdl5mZTY5Kg4ek55NurX2Y1Dfjloiov1vKzMw6XJXPMH8zqZd5sW3jUUk3kBrMf+c2CzOzqaHKfh4fJvX12IrUX+MfgW+S2jreBJxRYV5mZtZCVQaPZwPnkx9FC/wyIk4k9cu4F9ivwrzMzKyFqgweDwKKiLX5/dMBIuImUq/wT1eYl5mZtVCVweN60tUHwG3AXoV5NwOvqDAvMzNroarHttohv/8GcIqkg/NAhicA91eYl5mZtVAld1spDTT144g4Nyf9O/Ba4OekgQ3XAn9fRV5mZtZ6Vd2qOwN4SNLhEXFRHsPqEEn7kwY6vDoPz25mZlNAJcEjItZIuqdB+q+BX1eRh5mZtY8q2zwWAcdWuD0zM2tTVQYPgAMl9Ut6esXbNTOzNlLl2FZ9wPakhvHjJS0GFpOHJmF0z/MwM7MOUOXzPHaQNAvYI0+7k/p29DHK53mYmVlnqHJgxBOBgYj4OekW3Vr6RqTneexeVV5mZtZaVbZ5fBF4Xn1iRKwGbgF+UGFeZmbWQuMKHpKeJmmT2sdhFn0xabwrMzObAsZ75XEy8LCkG0k9yY+V9DpJO9YttxmwZpx5mZlZmxhvm8ci4HZSA/muwNHAiUBIuo90p9XvScOx3zTOvMzMrE2MK3hExArgWwCSjiQ99Gkp6+642hPYG1gFzB9PXmZm1j6qvFV3TuHjL/JkZmZT0IjBI3f2uzpPVwHXRcRjE10wMzNrX6NpMN8LeBvwZeAK0ui510r6hqSTJL1UUteElrIDDAwM0NPTw4wZM+jp6WFgYKDVRTIzmzCjqbZ6GymA7EXq6Lc58KI81QZCXCvpD6Qrk9pVyu8iYlrcnjswMEBfXx9DQ0MADA4O0tfXB0Bvb28ri2ZmNiEUEaNfOD306fmsCyZ7kRrFty4sVtzg0oh47viLOX5z586NxYsXT8i2e3p6GBwc3CC9u7ubZcuWTUieZmaTQdJVETG3Pr1Ug3mkSHNznp6sl5H0bNYPKC8GtmXdM82ntOXLl5dKNzPrdFU9DGop6Rbd79XSJO1ECiJT3pw5cxpeecyZM6fB0mZmna/q53k8KSJWRMQFE7X9drJw4UK6uta/Z6Crq4uFCxe2qERmZhNrwoLHdNLb20t/fz/d3d1Ioru7m/7+fjeWm9mUVarBvJNNZIO5mdlU1azB3FceZmZWmoOHmZmV5uBhZmaltV3wkHSIpFslLZF0aoP5vZKuy9PlkvZoRTnNzKaztgoekmaSxtA6lPTc86Mk7Va32O3AgRGxO/BxoH9yS2lmZm0VPIB9gCURcVtEPA6cAxxeXCAiLo+I+/LHK4D6pxaamdkEa7fgMRtYUfi8Mqc18zbgJ81mSuqTtFjS4lWrVlVURDMza7fgoQZpDTuiSDqIFDw+0GxjEdEfEXMjYu6sWbMqKqKZmVX2JMGKrAR2KnzeEbijfiFJuwNfBw6NiD9PUtnMzCxrtyuPK4FdJO0saRPgLcB642NJmgOcCxwTEb9vQRnNzKa9trryiIjVkk4CfgrMBBZFxI2STsjzzwT+GXg68JX0eBFWN+o6b2ZmE8djW5mZWVMe28rMzCrj4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBhZmalOXiYmVlpDh5mZlaag4eZmZXm4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBhZmalOXiYmVlpDh5mZlaag4eZmZXm4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBhZmaltV3wkHSIpFslLZF0aoP5kvSFPP86SXtNVFkGBgbo6elhxowZ9PT0MDAwMFFZmZl1lI1aXYAiSTOBLwOvAlYCV0q6ICJuKix2KLBLnvYFvppfKzUwMEBfXx9DQ0MADA4O0tfXB0Bvb2/V2ZmZdZR2u/LYB1gSEbdFxOPAOcDhdcscDnwrkiuAbSTtUHVB5s+f/2TgqBkaGmL+/PlVZ2Vm1nHaLXjMBlYUPq/MaWWXAUBSn6TFkhavWrWqVEGWL19eKt3MbDppt+ChBmkxhmVSYkR/RMyNiLmzZs0qVZA5c+aUSjczm07aLXisBHYqfN4RuGMMy4zbwoUL6erqWi+tq6uLhQsXVp2VmVnHabfgcSWwi6SdJW0CvAW4oG6ZC4Bj811X+wEPRMSdVRekt7eX/v5+uru7kUR3dzf9/f1uLDczo83utoqI1ZJOAn4KzAQWRcSNkk7I888ELgQOA5YAQ8DxE1We3t5eBwszswbaKngARMSFpABRTDuz8D6Ad052uczMbJ12q7YyM7MO4OBhZmalOXiYmVlpDh5mZlaaUvvz1CdpFTDY6nKMw7bAn1pdiAng/eo8U3XfvF+NdUfEBr2sp03w6HSSFkfE3FaXo2rer84zVffN+1WOq63MzKw0Bw8zMyvNwaNz9Le6ABPE+9V5puq+eb9KcJuHmZmV5isPMzMrzcHDzMxKc/Boc5KWSbpe0jWSFre6POMhaZGkeyTdUEh7mqSLJf0hvz61lWUciyb7tUDSH/Nxu0bSYa0s41hI2knSJZJulnSjpHfn9I4+ZsPs11Q4ZptK+q2ka/O+nZbTKz9mbvNoc5KWAXMjouM7L0k6AHiY9Az6F+a0zwD3RsSnJZ0KPDUiPtDKcpbVZL8WAA9HxGdbWbbxkLQDsENEXC1pS+Aq4AjgODr4mA2zX39L5x8zAZtHxMOSNgb+F3g38DdUfMx85WGTJiIuA+6tSz4c+GZ+/03SP3FHabJfHS8i7oyIq/P7h4Cbgdl0+DEbZr86XiQP548b5ymYgGPm4NH+AviZpKsk9bW6MBNg+9qTIPPrdi0uT5VOknRdrtbqqKqdepJ6gBcDv2EKHbO6/YIpcMwkzZR0DXAPcHFETMgxc/Bofy+LiL2AQ4F35ioSa39fBZ4N7AncCZzR0tKMg6QtgB8Ap0TEg60uT1Ua7NeUOGYRsSYi9gR2BPaR9MKJyMfBo81FxB359R7gPGCf1paocnfnOuhaXfQ9LS5PJSLi7vxPvBb4Gh163HK9+Q+AgYg4Nyd3/DFrtF9T5ZjVRMT9wKXAIUzAMXPwaGOSNs8NekjaHHg1cMPwa3WcC4C35vdvBX7YwrJUpvaPmr2BDjxuufH1P4CbI+JfC7M6+pg1268pcsxmSdomv98MeCVwCxNwzHy3VRuT9CzS1Qak583/Z0QsbGGRxkXSfwHzSENE3w18FDgf+C4wB1gOHBkRHdX43GS/5pGqPwJYBry9VufcKSS9HPgVcD2wNid/iNQ+0LHHbJj9OorOP2a7kxrEZ5IuDr4bER+T9HQqPmYOHmZmVpqrrczMrDQHDzMzK83Bw8zMSnPwMDOz0hw8zMysNAcPazuSdpYUkhYV0vbPaR9tZdlGI4+EHHnoi8nKMyT51kmbNA4e1o72y6+/LqS9tEGambXIRq0ugFkDteBxeSFtf1Lnrd9suLgBu7a6ADa9OHhYO9oPeAC4qZC2P3BTRDzQmiK1t4i4pdVlsOnF1VbWViRtQhoi4orIwx9I6gaeCVwxxm0uyG0CC3J7ytmS7pb0aH7a2vskNfwhpeQYSZdKui+vs1TSlyXtVLIc3ZI+mJ9it0LSY5LuzZ+PbrLOvFz2SyV1SfqEpFsk/SUPu11bboM2D0ln1dKHmS5tkOdOkj4v6dacz4OS/k/ScXlcqPrlL83bmifpAEn/LelPktZKOiIvs6mkUyVdLenhvO93Svp13qdNy3yX1nq+8rCWa9LQ+5oG6W+T9Lb8fjAiekpmtTOwGHiUNNroVsBBwGeBl0t6Yx5RtVYuAWcDRwNP5HXuJY22eiLwFkmHRMSVo8z/GODjwFLSYHX/Rxo2+xXAPEn7RsS7m6y7ac5/V+Ay4FpgkxHy+99h5r2cNPz4mmKipINI46ltDSwBLgK2IF0NfgM4GDi2yTaPBE4gXTFeTBrr6wlJM4D/zus+APwyv24PPA+YD3wJuGuE/bF2EhGePLV0Ar5emG4itW18v5C2PKd9o5D26RLbX5DXr21308K8XYCVed6JdeudmNPvAl5QSJ8JfIF1A+g9pW69ZXleT136S4rbqStDbR/3rZs3r1D235Ee6tNoHyP9O4/q+5gHPAYMFfMDdiAFx9WkkVdVmLdTzj+A4+q2d2mhjH0N8jsgz7uK9IjU4jwBLwO6Wv136Knc1PICePJUnEi/xh8GNimkLQduHcc2a8FjCNiuwfzj8/w/1KUvzen/0GCdTQon/N66eQ2Dxwhl/Ie8zr/UpReDx0uHWX9UwYN05XIf6YrjDXXzTs/bOb3JunNrQaAuvRY8ftZkvSPz/M+1+u/LU3WTq62sbUjailQldHFEPJ7TdiH96j2zgix+FumhWvUGSFczz5E0OyL+KGlH4FmkIbu/Xb9CRDwuaQA4lXSCHxhNAXLd/mtIVyGzgKfkWbVnSTy3yap3R8TlTeaNiqTtgQuBbYD3RMR5dYscll+/12QTV5EC+56SNo2IR+vmn9tgHYCrScHqbZJ+D/wgIu4uW35rLw4e1k7mkf4mf15IOzi//qKC7d/eKDEHgjuB2aQ2iD/m9wB3NjhJ1izNr7ObzF+PpP1Jz1TYcZjFtmqSPjiaPIbJuwv4EdADfCEiPtdgsWfl1ysbtIvXezrpeypqWMaIWCrpPaS2pS8DX5Z0G+lW7B8C50XEmkbrWvty8LCWkTQPuKTBrDMk1T8/+ruFE9o3I+K4CSpWrZFedZ8bGfEM++SC6eR9HqmR+D9Iz8teAjwUEWslvRr46TDb/Mto82qQ9wzgP0lXOz8E3tNk0Zn59TukmwqG81iDtKZljIgvSvoecASpsf7lwP/L0zWSDowp9Hz06cDBw1rpLtJTz2reRGqsPT9/FukOpbuAnxWWG+4uouH0NErMtwfXqo3uyK8r8+szJT0lIhqdLHfOr/W/wBs5gBQ4roqIv28w/zmj2MZY/RtwOPBb4Ogo3FFWZ0Uux8cj4saqCxERd5GqH88EkLQHqUpwT1L134eqztMmjoOHtUykjm3HASg9d/kYUn14Le1FpNtCz4mI91aQ5aslzYqIVXXpR5H6PC2NiJW5bCtz1cqzSL+O/6O4gqSNSbfwQmowHsnT8uuKJvMb9vMYL0mnAO8iVdm9PiKGhln8J8DJpAbuyoNHvYi4VtLnSe1Ne0x0flYtdxK0djGP9PdYrMY6KL9W0d4B0AV8SVKtkRpJzyb1vQD4fN3y/5pfPy7p+YV1ZgKfAbpJ9fzfH0XetR7gB9dta4akfybdrlqp3EHvDNLdVYc1uVmg6F+AB4EPSXpno46TkvaTdGTJchws6bD67eXvsdZIP642HZt8vvKwdlFrGK8PHmtIneKq8G3gtcBSSf9H6vx2MKkD3o9IjblFXyGd1I8CrpV0CelEvA/piuQ+4MgmVVrriYirJf0IeD2pjv8SUke5lwBzSMHon8a9h+s7gxSQB4FTmzSC3xIRn85lXJEDzvdJnfbmS7oR+DOph/+z8+t3aH5HViO7k6rOHpB0NXAnKZDvS6ouvIt0m7B1EAcPaxcHAXflqqxa7+4DgMUVNqTeRjpZf5IUNLbOaYtIfRDWawuIiJDUS6rO+QdSL+vNSO0iXwU+FRHNqqEaeROpsfoY0pXWw6RRgo/O2606eNQawPfMUyO/BD5d+xARl0h6Aan66rWkfd6YdIL/A/BFygUOSIF5G9LxfA5phOSHSf1kzgS+2qAq0dqcIvwIAJvaJC0APgqcFhELWlsas6nBbR5mZlaag4eZmZXm4GFmZqW5zcPMzErzlYeZmZXm4GFmZqU5eJiZWWkOHmZmVpqDh5mZlfb/ASLqkALiG3TNAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the result\n",
    "plt.figure()\n",
    "plt.plot(Narray, I, 'ko')\n",
    "plt.xlabel('# polarizers', fontsize=22)\n",
    "plt.ylabel('$I_{transmited}$', fontsize=22)\n",
    "plt.title('Perfect polarizers', fontsize=24);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to show how the transmitted intensity gets closer to 1 using logarithmic scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Npoins = 100\n",
    "Narray = np.logspace(1,5, Npoins)\n",
    "J = Jones_matrix('System')\n",
    "I = np.zeros(Npoins)\n",
    "\n",
    "# For loop to calculate the intensities\n",
    "for ind, N in enumerate(Narray):\n",
    "    # Create the system of polarizers\n",
    "    angles = np.linspace(90*degrees, 0*degrees, int(N)) # Start by 90º and end in 0º because the rightest  element is the first to be crossed by the light wave\n",
    "    J.diattenuator_perfect(azimuth=angles)\n",
    "    # Multiply them\n",
    "    J.prod()\n",
    "    # Multiply the system matrix by the incident light wave to obtain the output\n",
    "    E_final = J * E0\n",
    "    # Store the intensity\n",
    "    I[ind] = E_final.parameters.intensity(out_number=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEtCAYAAAAr9UYgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAt/UlEQVR4nO3deZhcVZ3/8fcnIRGaVZKACnQ3ICq4gBhZXDDghqiD6yhGFFwygOAybmg7P1En7oyiMmKPgyi04oriiAg6Io4rYYnsGiCbbCEBAgYNSb6/P86p5KZS1V3VfbuW7s/ree5TVecu59Tt7vvte7ariMDMzKwsU9pdADMzm1gcWMzMrFQOLGZmVioHFjMzK5UDi5mZlcqBxczMSuXAYh1NycmSrpG0RlLkpb/dZes27Tx3ki7LeR/X6ryt9RxYDEnnFC46xWV1vqB/RtLubSreB4EvAvsDAu7Ky/pWF0TSTpJOk3Raq/M26yYOLFb0MJsu3HcD25Eu6O8BrpX0rDaU6R359V+Bnoh4VF6WtaEsOwEfzos1ZylwM3B/uwti42+rdhfAOspvI2JO5YOkHuCVwBdIF9XvStorIh5qRWEk7QLMyh//KzxNRNeKiDe0uwzWOr5jsboiYk1EnAu8PSc9CnhZC4uwTaEsD7YwXzMbAwcWa8R3gA35/dOKKyRNkXSspEslrZC0VtLtkr4t6eBaB8vtFJHbdqbkxvk/Srovp79TUgCLC/sU235OG2sZCvtuK+k9kn4raZWkv0u6VdKFkuZKmpa3uwy4rU55tijTMPn1V/bJn58p6X9yudfkNq2TJQ37tynpFZIuzvv9Q9JySUOSDmykHFXHmirpcElnSLpS0l2Fc3iBpCOG2Xdjo3xug/qUpJvyd7mv1na1zkUDS3+dc/lFSTfn/B7I5X+/pG3rlHfj8STtK+nrkpZJeljSDwvb7ZLbFq+T9Lf8e7Es/558VFJfs+d5UokIL5N8Ac4BArhsmG3uytsMFtK2By7N6UEKPvcXPq8HTq5xrNPy+q8DP8zv1wH35vfPAO4EVhSOdWdhec9Yy5D33Y8ULCrbPlwoQ2Xpz9v+YJjybFamEc51f+EYr8x5Rs734cK6C4Ctauw/JZ+3ynbF81b5vifWyXuz71RIf1LVd/478GBV2gfrHPOyvP69wC2F/VcD99XY7rhC2h41zmNxuW+YMr8CeKiwfg3wj8LnPwG7DnMOjgX+lt+vzsf6Yd6mD7i96hyvyr9blbQT2v1328lL2wvgpf0LIwQWUpVU5Y/q04X0C3LaQuAoYJucvhPwgfyHvh54ZtXxTsv7PZAvQieSGuYBdgF2yO/7K3/Iw5R9tGXYmdSgHMCtwNHA9Lxue+BZwNnA7oV9RixPA+d64zHyhfOnwJ553bakC/R66lzMgVPZFEA/BGyf03cj3VlWgsthNfatd5F+XN73JcCugAo/iw/lC+sG4OAax7ys8LNcChwJTMnrHltju+MaPE/TgMvzPtcAWxfWPR1Ym8v1SaCX1GNwKnAw8Pu838+GOQcP5DI9KacL2Du/Pztv8xfg2YXv8whSEP4Y8LJ2/9128tL2Anhp/8LIgeXkwh/kK3Pa8/Ln24Cd6+z3vrzN/1Sln1Y43rxhyjXshXyMZfh0Tl8B7NbgeRq2PM0eA7gOeESNbSrn535g20L6tmy6G/tEjf2mAr/O6y+vsb5mYGmgzP+W9/tajXWVgLG2cpGuc4zKdsc1mOdZefu7gb6qdf+X172rzr6PBP6at5ld5xzcQv4npMb+N+RtXjPan/NkX9zGYjUp6Zf0HtJFGGAJ8OP8/o359ZyIWFXnMN/Mr4dLmlpj/UrSf4ejNZYyHJtfPxsRfx1DGcbi9Ij4R430/yDdye0APL+Q/oKctpZNP5ONImI96b9pgGdLelRJ5az8zJ85zDY/jYjryshM0onAv5CqBl8VEUsK6/bO5XiIFHy2EBH3ku4EYfPzV/SlqN+7cXV+fXSTRbfM3Y2t6DmVRuUa7iDd/q/Nn5+RX9+VLwTD6QFmkP77LFoQEetGV9TRlyE3BFcuuheNIf+xuqxWYkSslnQ1cChwIKkdivweYGG+eNZyOamKaKu8fUPfT9I2wAmkKsH9SP/1V18fHjPMIX7XSD4NlOMw4Iz88ZSIuLxqk8rPfDpwm6R6h9ouv+5RZ/1w5b2IVKX2KUn7AN8Dfj9MILIqDixW9DCpkRJSVcDfSO0PlwJfrbqYVf6b2zEvI+mpkbZilOUcaxl2LaQtHWMZxmK4O6XKulmFtFlV67YQEX+XtJL0HWfV265I0qNJQe5xheS/kToFbCBVsc0kVcXVM9afJZJ6SRfxacCXI+IrNTar/MynsvnPsZ5av3cwfHk/Rer9+E/ASXlZJ+kKUpvef0XEfQ3kPWk5sFjRZgMkR1CpRj06Ii4cZX5jnZZltGWo+29uBxmujI8oOa/Pk4LKraTOA78s/hORq58WjXCMMf0slQbj/ogUDH/FphkXqlV+5ldHRNNdqwvqljdXTx4t6RDg5cBhpEBzaF7eK+n5EbFwDPlPaG5jsdG6K7/u14VluLPwvq+ksozGcFVLlf/Mi/9ZV97XLbOkrUlVftX71tt+Oqn6C2BuRPygRjVbI3cGY/U14ADS2KVXRcTDdbar/Mz3kTSu/xhHxO8j4v0RcSipavAY0h3uLOCr45l3t3NgsdGq1FG/stvKEBGL2RRcjmpi18ogUTRM5X4TnlMrUdL2bGpPuaqwqvJ+H0m71TnmYWyqibiqzjZFM9l0B3R1nW2e18BxRk3SAPDPpOq3oyPinmE2r/zMtyN1ZmiJiPhbRJwPzMtJT6s3CNMcWGz0zsmvsyUNOw+UpEd2YBnOza/vHuYiXW114f1ODe4znHfnO4Zq7wS2zvldUki/JKdNI1VZbSb3evu3/PHXEXFn9TY1rCa1pwE8ucYxHw2c0sBxRkXSS4GP5jK8ISL+NNz2EXETaZwKpMb1uhd3SdtIarrasM7PpKLSgC9SBwKrwYHFRiUiLiaNRgc4W9JH8kUISBdySUdL+hGp+2ynleFTpEbwmcCvJf1T5YIiaTtJcySdr8LjAnKD7e354/ElfIVe4ILKdCWSeiT9K5tmT/5URKwp5P834OP549slDUjaLu+7G/At0sDOyuDJEUWag61yoT5b0gH5eFMkPZfU3jEubVKSHgecR7oOfSQifjDCLhWnkAa+Pon0s3tepVosl/uJkj5EGqsymi7D10n6uKSnF34nJOkg0iMcAK4YpmeetXsgjZf2LzQwpUud/bZl08j34mjy+6vSvla132k5/ZwRjt9fOUbZZcj7PhlYVthmLXWmdCns85HCugdJbQKLgXc2eM76C/sPN6XLD6k9pctUtpzSpTjdyHrgpDp51/tOB5OmRCl+r8rnlaQ2mJo/Bxoc+FhrO+C4Qp53Mfz0LntUHe9FbD7lyz+Ae/LPsPjz62vkHFRtUzzuunwOisddATyl3X+3nby4V5iNWqT/oF8u6cXAm0gXqFmki9wi4I+kO4pxGysyljJExLWSnkj6D/hlpJ5RW5N6R/0J+DawvGq3j5LaAuYCj2VTQ/pOoyj79yUdTpqm5RDSRex6UsPwf0bEhhr7rAfeKOlCUn3/00iDJu8g3V2cHhFXNlmOP0g6lBTwn0MK1ncAFwPzScFsvO0ywvrNyhARP813PKeQ2skeS/oZ3Ed67svFwHejMLiyCUcDLyS1V/WSOi+sBW4k/R59LiKqx2RZQWVOIDNrgVztdRtARHRDt2ezprmNxczMSuXAYmZmpXJgMTOzUjmwmJlZqdx4b2ZmpZr03Y1nzpwZ/f397S6GmVlXufLKK++JiJozaE/6wNLf38+CBQvaXQwzs64iqe4YIbexmJlZqRxYzMysVA4sZmZWKgcWMzMrVccFFklnS7pb0nV11kvSFyQtkvQnSQcW1h0p6ea87tTWldrMJruhoSH6+/uZMmUKM2fOZObMmS1/39/fz0knndRUOfr7+xkaGir3ZLR7euUaU1YfRnp63nV11h8F/JT0jIhDgD/k9Kmk5y/sRXoAz0Jgv5Hye9rTnhZmk8l5550XfX19ISlmzJgRM2bMaOv7vr6+OPHEEzuqTM2WFQhJxen6u2rp6emJ8847r6nfI2BB1LuO11vRzoX0zIp6geUrwDGFzzeTHuZzKPCzQvoHgA+MlJcDi3WqsgJA8WLY7RdAL+O39PX1NfX7yTCBpRvHsexGejhTxfKcViv94FoHkDSP/Ozq3t7e8SmlWZWhoSEGBgZYunQpO++8MwCrVq2q+X7lypVIqvyTxMqVKzcep9n3S5Ys4ctf/vLGz5VjmhUtXbq0tGN1Y2Cp9QyLGCZ9y8SIQWAQYPbs2f4rs1EZ70DhAGCtVOY/2d0YWJYDexQ+7056Dvn0OulmozZc8HjggQdYu3Yt4EBh3a2np4f58+eXdryO6xXWgAuBN+TeYYcA90fEHcAVwD6S9pQ0HXht3tZsRLV69Eji2GOPZcmSJUQEK1euZOXKlRvfV4KKWYWUKk5mzJjBjBkzkNTS9319fZx44on09fU1tc/g4CBz584t7Tx03B2LpG8Bc4CZkpYDHwamAUTEWaRnTh9Fep75GuD4vG6dpJOBn5F6iJ0dEde3/AtY16jcjSxZsqRuNdVEvcuofN8ZM2YA9avwWvG+t7eXo446iosuuqihasVOLWtvby/z588v9QLdrTousETEMSOsD+BtddZdRAo8ZhvVqs6qbvPo1ABSRgCovhj6AmjjreMCi1kZOvVupJlA4QBg3cqBxSaMesFkvALItGnT2GGHHZq6c3CgsMnAgcW6WiuCSa27DAcJs/ocWKzrjFcwcQAxK4cDi3WF8Q4mfX19DiBmJXFgsY5VZjDx3YhZ6ziwWMeo7hZcHNk+lmDiuxGz1nJgsY4wNDTEvHnzWLNmDbB5t+BmOJiYtV83TuliE0hlKpXXv/71G4NKsyrTaPT19XHuuecSESxevNhBxaxNfMdiLVev7aQZvjMx61wOLNZS1VVezQQVBxOz7uCqMGuJ0VR5TZs2bbMZWF3NZdYdfMdi4676LqURvisx614OLDZuim0pjerp6Sn92RBm1lquCrNxUblLaSSoFHt1OaiYdT8HFitVs20pbjsxm3hcFWalaaYtxVVeZhOX71isNAMDAw3fpTiomE1cDiw2ZpXqr5HaU3p6ejjvvPNc5WU2wTmw2Jg02kjvuxSzycNtLDYqjXYldluK2eTjOxZrmu9SzGw4vmOxpjXSSN/X18fixYtbUyAz6yi+Y7GGNdNIP3/+/BaVysw6jQOLNcTVX2bWKFeFWUNGqv5yI72ZVfiOxYbVSPWX71LMrMh3LFZXI1O0uJHezKr5jsXqaqT6y430ZlbNgcW24OovMxsLV4XZZlz9ZWZj5TsW24yrv8xsrDousEg6UtLNkhZJOrXG+kdKukDSnyT9UdKTCusWS7pW0jWSFrS25BPD0qVL665z9ZeZNaKjqsIkTQXOBJ4PLAeukHRhRNxQ2OyDwDUR8XJJT8jbP7ew/vCIuKdlhZ4gKpNKRkTN9a7+MrNGddody0HAooi4NSLWAucDR1dtsx/wC4CIuAnol7Rra4s5sYw0qt7VX2bWjE4LLLsBywqfl+e0ooXAKwAkHQT0AbvndQFcIulKSfPqZSJpnqQFkhasWLGitMJ3q+HaVVz9ZWbN6qiqMEA10qrrZj4JnCHpGuBa4GpgXV73zIi4XdIuwKWSboqIy7c4YMQgMAgwe/bs2nU/k0i9dhVJrv4ys6Z12h3LcmCPwufdgduLG0TE6og4PiIOAN4AzAJuy+tuz693AxeQqtasjsp4lXrtKr29vS0ukZlNBJ0WWK4A9pG0p6TpwGuBC4sbSNoprwN4C3B5RKyWtK2k7fM22wIvAK5rYdm7ittVzGy8dFRVWESsk3Qy8DNgKnB2RFwv6YS8/ixgX+AbktYDNwBvzrvvClwgCdL3+mZEXNzq79AtRmpXmT9/vttVzGxUVK8aZLKYPXt2LFgw+Ya8TJkypWYVmCQ2bNjQhhKZWTeRdGVEzK61rtOqwqxF6rWfuF3FzMbKgWWSKU4wmasNN3K7ipmVwYFlEqlusI+IjcHF41XMrCwd1Xhv46tWg31EeLoWMyuV71gmkXoDIYebeNLMrFkOLJOIG+zNrBUcWCa4SmP9lClTePDBB5k+ffpm691gb2Zlc2CZwIqN9RHBypUriQhmzJiBJDfYm9m4cOP9BFarsf7hhx9mu+224557/MgaMxsfYwoskn5Bml34auAa4MaI8LDtDuHGejNrh7HesdwFvBh4J2nK+39IupYUZCoB508R8dAY87FR6O3trTnJpBvrzWw8jamNJSJeFxH7AjsCzyNNe38g8BLg88BvgdWSbpA0NMayWoM8ut7M2qmUxvuI+BvwMuAWYJeI2A3YgTTt/Q2kpzzOKiMvG55H15tZu5XZeH8c8MqIuBcgP7P+u5J+SJoGf4snOVr5PLrezNqtzO7Ga4BHVSdGxMPA6cDxJeZldbjB3szarczAMgR8RNJuNdZtDTy6xLysDo+uN7N2KzOwfAi4FbhZ0uckHSlpP0kvBz4NXFtiXlbH/Pnz6enp2SzNDfZm1kqlBZaI+DvwQuDDwD8DF5GCyffzJvPKysu2VOkJduyxx7LNNtt4dL2ZtU2pI+/z4MjTgdMl7Qs8BlgJXBcR68rMyzap9ASrNNqvXLmSnp4ezj33XAcUM2u5Up95L+kw4E2k7sX3A2+PiI5uNZ4Iz7yvjFmp5p5gZjZeWvLMe0lzgcuAPYDrgJcCM/O60yQ9p6y8bHPuCWZmnaTsxvvPRsRzgXeRpnipeJg07YuNA/cEM7NOUmZg6Qcuzu+r69euI031YuPAPcHMrJOUGVhuAZ5cZ929eEqXcTN37lwGBwfp6+tzTzAza7sye4V9Afi4pF8DC6vW7QOsKDEvqzJ37lwHEjPrCGWOYxkEfkWa0fjrpOqwR+WeYgPApWXlZUnxscP9/f0MDXkCaTNrv7KfIPkq4G3A+0iN9z/Or1eQgouVpHrsypIlS5g3L41B9Z2LmbVTaeNYJJ0EDEXE/flzZYDkHaQnS5Y3YKZE3TqOxWNXzKydWjKOBfgi8PjKh4i4MSJ+ERE3ANMkbVNiXpOex66YWacaU2CRtLOk6ZWPw2z6VGD1WPKyzXnsipl1qrHesZwCPCjpelJj/RskvUTS7lXbbQOsH2NeVuCxK2bWqcbaeH82cBuwP7Av8DrgJCAk3QtcA/wZOIT0iGIbo6GhIQYGBli6dCk777wz22yzDatWraK3t5f58+e74d7M2m5MdywRsSwivhER7waWA0eSBkI+H5gP/BV4GmkMS0PT5ufnuNwsaZGkU2usf6SkCyT9SdIfJT2p0X27XfF59hHBypUreeihhzj33HNZvHixg4qZdYQx9QqTtFNE3FdaYaSppDuc55MC1RXAMbkDQGWbzwAPRsRHJD0BODMintvIvrV0U68w9wQzs04xnr3CVkjaP2fyAUkvqvNo4kYdBCyKiFsjYi1wPnB01Tb7Ab8AiIibgH5Juza4b1dzTzAz6wZjDSyvAe7O798B/ARYKukeSb+Q9B+S3ijpAEnTGjjebsCywuflOa1oIfAKAEkHkZ79snuD+3Y19wQzs24w1jaWH0TEHfn9o4BdSY8n/gTpwn4EMAhcBTzYwCFrdVmurqv7JPBISdeQeqVdDaxrcN+UiTRP0gJJC1as6J4pzNwTzMy6QdmPJl4B/DwvAEjailR99ZQGDrGc9KCwit2B26vyWA0cn48tUq+024CekfYtHGOQFPCYPXt2R84IUEulcb7SK8w9wcysE5U9V9gW8rPu/5SXkVwB7CNpT1KPsteSujBvJGknYE1uR3kLcHlErJY04r7dqtjF2MHEzDpdaYElN+J/jjSe5T7Sw70WksayLIyIZXV3ziJinaSTgZ8BU4GzI+J6SSfk9Wfl439D0nrS2Jg3D7dvWd+vXTzZpJl1mzInoVwI7EAaNLk9qeprf1K7SwD3RcSMUjIrUad3N3YXYzPrRMN1Ny6zKmxv4DUR8ZOqzHcBDqCxNhar4i7GZtZtypzd+Apgu+rEiLg7Ii6JiM+WmNek4S7GZtZtygwsHwBOlbRjicec9NzF2My6TZmB5Y/A74DrJL1f0tMlbV3i8SeluXPnMjg4SF9fH5Lo6+tjcHDQDfdm1rHKbLz/KvAm4AHgEcB00lT5fyb1DLu6E6vDOr3x3sysE7XqCZKvAeZHxI7AtqQeYccDPyX1DHtfiXlNeENDQ/T39zNlyhT6+/sZGhpqd5HMzBpSZq+we4FfAUTEeuDavJxXYh6TgseumFk3K/OO5SvAy0o83qQ1MDCwMahUrFmzhoGBgTaVyMyscWUGlmnASyR9RNIW3Y6tcR67YmbdrMzA8i9AL/BvwN2SfiXpDEnHS3pqg9PmGx67YmbdrbTAEhGPJjXSv4AUXBYDzwHOAq6ksWnzDY9dMbPu1mnT5hueHt/MuluZ41jmAesi4uxSDtgiHsdiZta8Vo1jeTswq5DpHEmrJD0g6RMl5mNmZh2szMDST3oEccUngFuA04ATJb2pxLzMzKxDlRlYHgI2AEjqAw4GTomI04GPAG8rMa8JyaPtzWwiKDOwLASen9/PBe6KiN/nz1cB+5SY14RTGW2/ZMkSImLjaHsHFzPrNmUGlvnAOyT9nNTd+NzCul2BtSXmNeF4tL2ZTRSldTeOiF9KOgI4BlhAqv6qmEOa5djq8Gh7M5soyh7H8jvSM1mqTQW+U2ZeE01vb2/NZ9t7tL2ZdZsyq8Lqioh/iYjPtyKvbuXR9mY2UZR2xyLpkcB7gccBK4DrSA36CyPigbLymag82t7MJooyR97/BDgC+F9gZ+CJwHZAkOYNuzoiXlVKZiXyyHszs+YNN/K+zDaW5wBvK07pIumxwAHAU/FcYWZmk0KZgeUO4K/FhIhYBCwCvldiPmZm1sHKbLz/MvDKEo9nZmZdqMzAsjVwhKR/l7R9icc1M7MuUmZgeSewF/BB0hMkL5f0BUlvknSgpOkl5jVheH4wM5toyhx5v4ukR5Ma6SvLHNIji6cBDwOPKCu/iaAyP1hlKpfK/GCAuxmbWdcqrbtx3QzSs+73A54cEeeNa2aj0M7uxv39/TVH2/f19bF48eLWF8jMrEGt6m5cU0Q8TB4oOd55dRvPD2ZmE1GpU7pIeoKkN0p6paR9JanM40809eYB8/xgZtbNSgsskl4DXEvqdvwd0pQuD0r6g6SvSDqhweMcKelmSYsknVpj/Y6SfixpoaTrJR1fWLdY0rWSrpHU8cPpPT+YmU1EZd6xfIgUVHYARJo37OvAk4FXAaePdABJU4EzgReR2mWOkbRf1WZvA26IiP1JnQNOr+pxdnhEHFCv7q+TzJ07l8HBQfr6+pBEX18fg4ODbrg3s65WZmDZG/gh+fHEwK8i4iRgNrAKOKSBYxwELIqIWyNiLXA+cHTVNgFsn6vZtsvHXjf24rfH3LlzWbx4MRs2bGDx4sUOKmbW9coMLKtJvcw25PczACLiBuDzwCcbOMZuwLLC5+U5rehLwL7A7aSqt3fkPCEFnUskXSlp3ii/h5mZjUGZgeVa0l0LwK3AgYV1NwLPbuAYtRr7q/tDvxC4BngMaYLLL0naIa97ZkQcSKpKe5ukw2pmIs2TtEDSghUrVjRQLDMza1TZc4U9Or//GvBOSUdI2g04AbivgWMsB/YofN6ddGdSdDzwg0gWAbcBTwCIiNvz693ABaSqtS1ExGBEzI6I2bNmzWrku5mZWYNKCSy5veN/IqLynPuvkO4qfg4sBV4O/FsDh7oC2EfSnrlB/rXAhVXbLAWem/PdFXg8cKukbStzlEnaFngBqWeamZm1UFkDJKcAD0g6OiIuzg3vR0o6lHQHclW+uxhWRKyTdDLwM2AqcHZEXF/pqhwRZwEfA86RdC2p6uz9EXGPpL2AC/LQma2Ab0bExSV9PzMza1ApgSUi1ku6u0b674DfNXmsi4CLqtLOKry/nXQ3Ur3frcD+zeRlZmblK7ON5WzgDSUez8zMulCpU7oAz5E0KGlGycedUDxVvplNZGVOQjkP2BV4C3B8nlJlAakR/2rgutz2Mql5qnwzm+hKnTZf0ixSO8f+pOex7E8azDgNeDgiOu55LK2eNt9T5ZvZRNCSafMlnQQMRcTPSd2MK+lbkeb9ekpZeXUzT5VvZhNdmW0sXySNKdlMRKwDbgK+X2JeXctT5ZvZRDemwCJp58LMwsM9e+WppPnDJj1PlW9mE91Y71hOIT1z5XrSnF5vkPQSSbtXbbcNsH6MeU0InirfzCa6MTXeS9oDOJzUSP8u0nxgO5GCzL2kHmF/Jk+ZnyeI7CjtfOa9mVm3GrfG+4hYBnwjZ/Jq0gO9bmFTz7ADgKcBK4CBseRlZmbdobReYRFRbH3+37yYmdkkM2IbS35uyaCkEyQ9XVLHjUXpBh5tb2aTRSN3LAeSenW9OX9eL+lG4Crgyvx6TUSsGZ8idj+PtjezyWTExntJx5OCy4GkQY7bFlZXdt4A/IVNgeYq4OqI6Pguxq1ovPdoezObaIZrvG+qV1h+oNcT2BRoDiQ10O9Y2Kx4wFsi4nHNFriVWhFYpkyZQq3zLIkNGzaMa95mZuOhtF5hka6ON+ZlYyOBpL3ZPNg8FZgJ7D3KMk8ovb29Ne9YPNrezCaiUqZ0iYhbIuK7EfGBiHhhROwC9JEeSTzpebS9mU0mZT+PZaOIWBYR1c+rn5Q82t7MJpNSp83vRh55b2bWvOHaWMbtjsXMzCYnBxYzMyuVA8s48mh7M5uMynzmvRV4tL2ZTVa+YxknAwMDG4NKxZo1axgY8CTPZjaxObCMEz/b3swmKweWceJn25vZZOXAMk482t7MJisHlnHi0fZmNll55L1H3puZNc0j783MrGUcWMzMrFQOLGZmVqqOCyySjpR0s6RFkk6tsX5HST+WtFDS9fnRyQ3ta2Zm46+jAoukqcCZwIuA/YBjJO1XtdnbgBsiYn9gDnC6pOkN7mtmZuOsowILcBCwKCJujYi1wPnA0VXbBLC9JAHbAauAdQ3uO+488aSZTXadNgnlbsCywuflwMFV23wJuBC4HdgeeE1EbJDUyL7jyhNPmpl13h2LaqRVD7R5IXAN8BjgAOBLknZocN+UiTRP0gJJC1asWDH60lbxxJNmZp0XWJYDexQ+7066Myk6HvhBJIuA24AnNLgvABExGBGzI2L2rFmzSiu8J540M+u8wHIFsI+kPSVNB15LqvYqWgo8F0DSrsDjgVsb3HdceeJJM7MOCywRsQ44GfgZcCPwnYi4XtIJkk7Im30MeIaka4FfAO+PiHvq7dvK8nviSTMzzxVW+lxhQ0NDDAwMsHTpUnp7e5k/f74b7s1swhlurjAHFk9CaWbWNE9CaWZmLePAYmZmpXJgMTOzUjmwmJlZqRxYzMysVA4sZmZWKgcWMzMrlQOLmZmVyoHFzMxK5cBiZmalcmAxM7NSObCYmVmpHFjMzKxUDixmZlYqBxYzMyuVA4uZmZXKgcXMzErlwGJmZqVyYDEzs1I5sJiZWakcWMzMrFQOLGZmVioHFjMzK5UDi5mZlcqBZRSGhobo7+9nypQpzJw5k5kzZzJlyhT6+/sZGhpqd/HMzNpqq3YXoNsMDQ0xb9481qxZA8DKlSs3rluyZAnz5s0DYO7cuW0pn5lZu/mOpUkDAwMbg0ota9asYWBgoIUlMjPrLA4sTVq6dGkp25iZTVQOLE3q7e0tZRszs4nKgaVJ8+fPp6enp+76np4e5s+f38ISmZl1FgeWJs2dO5fBwUH6+vqQxIwZM5gxYwaS6OvrY3Bw0A33ZjapKSLaXYbNSDoSOAOYCnw1Ij5Ztf69QOXKvRWwLzArIlZJWgw8AKwH1kXE7JHymz17dixYsKDEb2BmNvFJurLeNbajuhtLmgqcCTwfWA5cIenCiLihsk1EfAb4TN7+pcC7ImJV4TCHR8Q9LSy2mZkVdFpV2EHAooi4NSLWAucDRw+z/THAt1pSMjMza0inBZbdgGWFz8tz2hYk9QBHAt8vJAdwiaQrJc0bt1KamVldHVUVBqhGWr1GoJcCv6mqBntmRNwuaRfgUkk3RcTlW2SSgs48cNdgM7Oyddody3Jgj8Ln3YHb62z7WqqqwSLi9vx6N3ABqWptCxExGBGzI2L2rFmzxlxoMzPbpKN6hUnaCvgz8Fzgr8AVwOsi4vqq7XYEbgP2iIi/5bRtgSkR8UB+fynw0Yi4eIQ8VwBLCkk7Avc3+HkmMF4dBarzLXO/kbapt75Wus/X+J0vGL9z5vPVvNGcs4l8vvoiovZ/5hHRUQtwFCm43AIM5LQTgBMK2xwHnF+1317AwrxcX9l3FPkPNvoZWDCO52FwvPYbaZt662ul+3yN3/kaz3Pm89WaczZZz1entbEQERcBF1WlnVX1+RzgnKq0W4H9SyjCj5v8PF5Gm08j+420Tb31tdJ9vny+mt2mG8/XaPOalOero6rCuo2kBdHAIExLfL6a53PWHJ+v5ozX+eq0xvtuM9juAnQZn6/m+Zw1x+erOeNyvnzHYmZmpfIdi5mZlcqBxczMSuXAYmZmpXJgKZGkvST9t6Tvtbss3UDSyyT9l6QfSXpBu8vT6STtK+ksSd+TdGK7y9MNJG2b5w58SbvL0ukkzZH06/w7Nmcsx3JgGYGksyXdLem6qvQjJd0saZGkUyGNpYmIN7enpJ2hyfP1w4h4K2nA62vaUNy2a/J83RgRJwD/DEzKLrXNnK/s/cB3WlvKztHk+QrgQWBr0vRaozdeo1QnygIcBhwIXFdIm0qaGWAvYDpptP9+hfXfa3e5u+x8nQ4c2O6yd8P5Av4J+C1pqqO2l7+TzxfwPNKcgscBL2l32bvgfE3J63cFhsaSr+9YRhBpduRVVcnNPjdm0mjmfCn5FPDTiLiq1WXtBM3+fkXEhRHxDDY9RXVSafJ8HQ4cArwOeKukSXe9a+Z8RcSGvP5e4BFjybfjpnTpErWeG3OwpBnAfOCpkj4QEZ9oS+k6T83zBZxC+q9yR0mPjaqpeyaxer9fc4BXkP7oL9pyt0mr5vmKiJMBJB0H3FO4cE529X6/XgG8ENgJ+NJYMnBgGZ2az42JiJWkCTNtc/XO1xeAL7S6MF2g3vm6DListUXpCsM+xynS3IK2Sb3frx8APygjg0l3a1iSZp4bYz5fzfL5ao7PV3PG/Xw5sIzOFcA+kvaUNJ3UQHhhm8vUyXy+muPz1Ryfr+aM+/lyYBmBpG8BvwMeL2m5pDdHxDrgZOBnwI3Ad6LqYWSTlc9Xc3y+muPz1Zx2nS9PQmlmZqXyHYuZmZXKgcXMzErlwGJmZqVyYDEzs1I5sJiZWakcWMzMrFQOLNZV8qCukHR2Ie3QnPbhdpatEZIW57L2tzDPkORxBdYyDizWbQ7Jr78rpD2jRpqZtYknobRuUwksvy2kHUqadPAPrS9OV9i33QWwycWBxbrNIcD9wA2FtEOBGyLi/vYUqbNFxE3tLoNNLq4Ks66RJ8w7APh95LmIJPUBjwF+P8pjnpbbIE7L7TfnSbpL0t8lXS/p3ZJq/gOWH1R2rKTLJN2b97lF0pmS9qi1zzDl6JP0AUm/lLRM0j8krcqfX1dnnzm57JdJ6pH075JukvSQpGsK223RxiLpnEr6MMtlNfLcQ9IZ+bG2D0laLek3ko6TtMV07Llskct6mKSfSLpH0gZJL8vbbC3pVElXSXowf/c7JP0uf6etmzmX1n6+Y7GOVqfR+YU10t8s6c35/ZKI6G8yqz2BBcDfSc882YH0BMLPAs+S9Mrig6LyRfQ80tMJH877rCI9ne8k4LWSjoyIKxrM/1jgY6RHxt4E/IY0nfmzgTmSDo6Id9TZd+uc/77A5aRHzU4fIb//G2bds4C9gfXFREmHAxcAOwKLgIuB7Uh3kV8DjgDeUOeYryY9q+gG4FJgJvCw0lMdf5L3vR/4VX7dFXg8MEB66NSdI3wf6yTtfiazFy/DLcBXC8sNpLaU7xXSlua0rxXSPtnE8U/L+1eOu3Vh3T6kZ1cEcFLVfifl9DuBJxbSp5IeXhbAYuARVfstzuv6q9KfXjxOVRkq3/HgqnVzCmW/Gti1zneM9Kfe0PmYA/wDWFPMD3g0KXCuA95InsA2r9sj5x/AcVXHu6xQxnk18jssr7sS2LZqnYBnAj3t/j300tzS9gJ48dLoQvov/kFgeiFtKXDzGI5ZCSxrgF1qrD8+r/9LVfotOf2tNfaZXggGc6vW1QwsI5TxrXmfz1SlFwPLM4bZv6HAQrrjuZd0p/LyqnWfysf5VJ19Z1cCRFV6JbBcUme/V+f1n2/375eX8hZXhVlXkLQDqZrp0ohYm9P2If23fFYJWVwSEXfXSB8i3QU9VtJuEfFXSbsDewEbgHOrd4iItZKGgFNJF/+hRgqQ2xJeSLp7mUV6tj2kuwWAx9XZ9a6I+G2ddQ2RtCtwEel55++KiAuqNjkqv363ziGuJAX9AyRtHRF/r1pf75G3V5EC2Zsl/Rn4fkTc1Wz5rbM4sFi3mEP6ff15Ie2I/Pq/JRz/tlqJOUjcAexGavP4a34PcEeNC2jFLfl1tzrrNyPpUOA7OY96dqiTvqSRPIbJuwf4MdAPfCEiPl9js73y6xU12uirzSCdp6KaZYyIWyS9i9SWdSZwpqRbSd3JfwRcEBHra+1rncuBxTqSpDnAL2usOl3S6VVp3ylc7L4eEceNU7EqHQZU9bmWEa++GzdMF/YLSA3W/w18mdQ4/kBEbJD0AtLT/uod86FG86qR9xTgm6S7pB8B76qz6dT8+m1SB4fh/KNGWt0yRsQXJX0XeBmp48CzgNfn5RpJz4mI1SPkaR3EgcU61Z3A1wufX0VqOP5h/ixST6o7gUsK2w3X22k4/bUScxfnSlXU7fl1eX59jKRHREStC+me+bX6P/daDiMFlSsj4i011j+2gWOM1ueAo4E/Aq+LQs+3KstyOT4W4/DY34i4k1SleRaApP1J1YwHkKoUP1h2njZ+HFisI0Ua1HccgKSdSEHk+5W7EUlPJnVtPT8i/rWELF8gaVZErKhKP4Y03uuWiFiey7Y8V9fsRfqv+r+LO0iaRuqGDKnxeiQ759dlddbXHMcyVpLeCbydVA340ohYM8zmPwVOITW2j/vz5CNioaQzSO1b+493flYuD5C0bjCH9LtarBo7PL+W0b4C0AN8SVKlwRxJe5PGlgCcUbX9f+TXj0l6QmGfqcCngT5Su8L3Gsi7MjL+iKpjTZH0/0hdbkuVByeeTuoFdlSdjgtFnwFWAx+U9LZag0YlHSLp1U2W4whJR1UfL5/HSoeBMbUhWev5jsW6QaWRvjqwrCcNCCzDucCLgVsk/YY08O8I0uDDH5Malov+k3TBPwZYKOmXpIv0QaQ7mXuBV9epJttMRFwl6cfAS0ltCr8kDRJ8OtBLClTvG/M33NzppGC9BDi1ToP8TRHxyVzGZTkYfY80YHFA0vXAStLMB3vn129Tv+dYLU8hVcfdL+kq4A5SkD+YVAV5J6mrs3URBxbrBocDd+bqscqo98OABSU26t5KupB/nBRQdsxpZ5PGWGzW9hARIWkuqYroraTR59uQ2mG+DHwiIupVbdXyKlLD+bGkO7QHSbM1vy4ft+zAUmmMPyAvtfwK+GTlQ0T8UtITSVViLyZ952mki/9fgC/SXFCBFLR3Iv08H0uaqfpB0jigs4Av16ietA6nCD+mwSYvSacBHwY+EhGntbc0ZhOD21jMzKxUDixmZlYqBxYzMyuV21jMzKxUvmMxM7NSObCYmVmpHFjMzKxUDixmZlYqBxYzMyuVA4uZmZXq/wPC09b0qycWoQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the result\n",
    "plt.figure()\n",
    "plt.semilogx(Narray, I, 'ko')\n",
    "plt.xlabel('# polarizers', fontsize=22)\n",
    "plt.ylabel('$I_{transmited}$', fontsize=22)\n",
    "plt.title('Perfect polarizers', fontsize=24);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imperfect polarizers\n",
    "\n",
    "Imperfect polarizers will have a leak ($T_{min}$ > 0) and, much more important for the Zeno effect, will not perfectly transmit the transmission state ($T_{max}$ < 1). We will use some values of commercial polarizers with anti-reflecting coating: \n",
    "\n",
    "* $T_{max}$ = 0.99.\n",
    "* $T_{min}$ = 1e-4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nmax = 30\n",
    "Narray = np.arange(2,Nmax+1)\n",
    "J = Jones_matrix('System')\n",
    "I = np.zeros(Nmax-1)\n",
    "\n",
    "# For loop to calculate the intensities\n",
    "for ind, N in enumerate(Narray):\n",
    "    # Create the system of polarizers\n",
    "    angles = np.linspace(90*degrees, 0*degrees, N) # Start by 90º and end in 0º because the rightest  element is the first to be crossed by the light wave\n",
    "    J.diattenuator_linear(Tmax=0.99, Tmin=1e-4, azimuth=angles)\n",
    "    # Multiply them\n",
    "    J.prod()\n",
    "    # Multiply the system matrix by the incident light wave to obtain the output\n",
    "    E_final = J * E0\n",
    "    # Store the intensity\n",
    "    I[ind] = E_final.parameters.intensity(out_number=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum obtained intensity is 0.7201695688246019 for 17 polarizers.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEpCAYAAABr364UAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsTUlEQVR4nO3de5gcVZ3/8fcnAYQBEZWgGMgMIq6ACuII4gUDqy6yuvECKzI/EHV3FhEV3ItodEXdrJeVXRVRnNWou8wu6wUQFUVdBdwF1gwIyFUTTELkFrlDuCX5/v44p0ml093pStdMd08+r+epp7tPnapzqmumvl116pxSRGBmZlbGjG5XwMzM+o+Dh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBh04KSEyRdKWmVpMjTULfr1m+6+d1JujCXfexUl23lOHj0EUlfz/9YF3a7Lj3og8BpwD6AgNvztGaqKyJpB0mnSDplqss2myoOHjZdvDe/vg8YiIin5+nmLtRlB+AjebJylgM3Avd2uyLW2hbdroBZpyTtBMzKH/81PGxC34qIY7pdB2uPzzxsOtim9iYiHuhmRcw2Fw4e00SxoVHS9pI+LWmJpIck3STpY5K2LuT/Y0kXSPqDpAclXSzp5U3WfUpe99clzZB0kqSr8nJ3SjpP0v4bqd8MSUdL+omklZIelXSLpP+SdECb5Z4g6ZeS7snpJ0oKYGlhmShMp3Rah8Ky20r6G0mXSLpL0sP5ez1P0oikLWv7Afhdk/psUKcW5Q3VlsmfXyrp+7neq/KNASdIavk/LOmNkn6Ul3tE0gpJ45L2a6cedeuaKelgSZ+TdLmk2wvf4TmSDmmxbPHvcwdJn5J0Q96Wexrla/RdtDENNfkuT5N0Yy7v/lz/90vatkl9H1+fpD0lfUPSzZIek3RuId9Okv5J0jX5/+HhnO8Spf+5wbLfc9+ICE99MgFfBwK4sMG8C/O8k4Dr8/sHgEfz+wDOy3mPB9aSGpPvLcx/BHhpg3Wfkud/A/hOfv8YcE9h2dXAm5vU+4nATwp519aVuwY4YSPlnlso5+78/iXAbcDKwrpuK0x/02kd8rJ7kQJCLe9jhTrUpqGc9+wW9VmvThvZ10OFdbwplxm53McK884Btmiw/Iz8vRX3T7HOa4B3Nil7vW0qpD+3bpsfzn9jxbQPNlnnhXn+3wJLCsvfB9zTIN+xhbRdG3yPxemeFnV+I/BQYf4q0t957fPVwNNafAdHAw/m9/fldZ2b8wwCt9R9x3flv61a2nHdPm5M2vGo2xXwVGJntRc87gFuAF6W07cC/qJwwPkwKaD8I7BDzjMIXJLn/7LBuk8prHs1KUBtk+ftDvy48I+5e4Plz8nzrwIOKyy7A/CB/M+8hrrAVSj3/nygeSepMRxgJ2D7/H6o9s/a4rvb1Do8hdSIG8BNwDxgqzzvicDLgIXALoVlNlqfNvb14+vI3/sPgd3yvG1JB+E1NDlgAyezLkh+CHhiTp8NfJN1AeSgBss2OxA/Oy/7WuBpgAr74kP5b2MtcECLv8/78/d5KDAjz3tWg3zHtvk9bQlcnJe5Eti6MO9FpL/11cAngTmkO/FmAgcAl+XlLmjxHdyf6/TcnC7y33je7wH8Fnh5YXueQAq0Hwde3+3jxmRNXa+ApxI7q73g8Vjxn7Ew/6uFf4iFDeYPsu4X05y6eacUlp3fYNmtSQErgK/UzXtlTv8d8JQm2/V3Oc/3W5Q72uJ7GaLFwbrDOnw6p68EZre5n1rWp+w6gGuAJzTIU/t+7gW2LaRvy7qzqk80WG4m8Is8/+IG8xsGjzbq/OG83Nda/H0+WjsQN1lHLd+xbZZ5Rs5/BzBYN+9/8ryTmiz7ZOD3Oc9wk+9gCfmHRoPlr8t5Gp5xT/fJbR7Tz7ciYnGD9J8W3n+ifmZELANqyz23ybpXAZ9tsOzDwKn545skqTD7rfn16xFxV5P1/kd+PVjSzAbz7yT9yttUndTh6Pz6mYj4fQd16MSpEfFIg/R/Jp2RbQ+8qpD+6pz2KCn4rSci1pB+FQO8XNLTK6rn9/LrS1vk+WFEXFNFYZLeCfwV6QfT4flvuDZv91yPh0gBZgMRcTfpjA7W//6KvhARDzWZd19+3blk1acF36o7/fy6Sfod+fVh1gWJercDe5B+kTUyEREPNpl3UX7dAdiNdIkHUrsEwEn5n72VAeCphboWy129kWVb2aQ65MbX2oH1/A7K79SFjRIj4j5JvwIOBPYjtQuR3wNclQ+QjVxMupyzRc7f1vZJ2gY4jnT5bi/S30r9ceQZLVZxaTvltFGPg4DP5Y/vjoiL67LU9vlWwO/W/z2znu3y665N5req7/mky1+fkrQH8G3gshbBZlpx8Jh+bm2SXutpfXvkc+4WebZsMr/VL+/ivFmsCx61X2VPytPGDDRIW9nGcq1sah2eVkhb3mEdOtHO9z6rkDarbt4GIuJhSXeStnFWs3xFknYmBbJnF5IfJDXEryVdDtuRdNmsmU73JZLmkA7UWwJfiogvN8hW2+czWX8/NtPo7w5a1/dTwAuBPyPdhHI8sFrSIlIb279GxD1tlN2XfNnKqtLsp13tb2xeRKiNaWmDdXQ6xMim1qHpz9Ue0qqOT6i4rM+SAsdNpDvAnhIR20XEThHxdODFbayjo30paQD4LingXcS6kQXq1fb5r9rc58eWrW9EPBIR80hnfp9mXQN87fNvJO1TeiP7hIOHldHqckTxum/x19rt+XWv6qvTtk2tw22F94MV1WVTtPO9F7/z2vumdVbq8/PUBss2y78V6VIVwEhEnN3gklg7v/A79TVgX1LfnsMj4rEm+Wr7fA9Jk3qFJSIui4j3R8SBpMt4byGdqc4CvjKZZXeTg4eV8aL8y6+RV+TXeyh0kmPdNeM3TVal2rBJdchnILUAcliJRdfW3qjFxfYSXtEoUdITWde+cUVhVu39HpJmN1nnQay7bH1FkzxFO7LuTOZXTfK8so31bDJJ84E/J10qmxcRf2iRvbbPtyPdQDAlIuLBiDgLGM1JL2zWEbHfOXhYGQM0uEwg6QmkAQkBvl3XpvL1/DosqeW4RZKaNdR3qpM6/Ht+/esWB+J69xXe79DmMq38df7lX+9E0m3S95H62tT8OKdtSeoPsp58N9mH88dfRMRt9XkauI90SQbgeQ3WuTPw7jbWs0kkvQ74WK7DMRFxdav8EXED6TISpAbtpgdwSdvkv+GydWq0T2pqjeYiNdpPOw4eVsa9wMclvTffdYOkZ5KuQe9JupPrk8UFIuJHpF7XAAslfTQfaMjLP1nSPEnfJd16WrkO6/ApUsPzjsAvJP1Z7aAhaTtJcyWdJWmXQnn3kHoeA7ytgk2YA5xTG3pD0oCk97Fu1N5PRcSqQvkPkjqBArxH0nxJ2+VlZwP/SercWOtAuFGRxgyrHYwXSto3r2+GpD8mtT9MShuRpGcDZ5KOVx+NiLM3skjNu0mdP59L2nevrF3CyvXeW9KHSH05NuV222sk/aOkFxX+JqQ0VM9pOc+iFne89bdudzTx1P5Ee50Ej22y7Nw8f2mL9TdcB+sPE3I26zp73c36QzMc2WS927Kuh3dtuof1hwfZoHNZodyvb+R7Gaqto0WeTapDXvZ5wM2FPPXbvkGHOuCjhXkPkK7RLwVObHNfDxWWbzU8ybk0Hp5kJhsOT1IcOmMNcHyTsptt0wGkvj7F7ap9vpPUJtJwP2zs77NVPuDYQpm303qokl3r1vca1h++5BHgD6w/bE+wYQfDht9BXZ7ielfn76C43pXA87t5zJjMyWceVkYAR5AuUV1POh2/G/g+8JJI13o3XChdB34DaViLs0m/5LfJyy8mddA7nHSr4+RUvIM6RMSvgb1Jv9InSJcktibddXQuqYF0Rd1iHwPeTxo7SaTG60E24TJWRHwHOBj4Aemgv5o0zMq7gTdGgz4wEbEmIt6at+nHpAPddqRbuf8T2D8ivliyHv9HupPoXNJ+35LUJ+fLpEbsq8pu2ybYidQw32xar5NpRPyQdIfYP5Dadh4m7YP7SEPy/D2wZxQ6GJYwj9Th9n9JZ5rbkYLH1aQz8L1jI5fX+lltbBqzppRGgv0I8I1ofkujVShfovodQET0wy3DtpnxmYeZmZXm4GFmZqU5eJiZWWkOHmZmVtpm02C+4447xtDQULerYWbWVy6//PI/RMQGg2duNqPqDg0NMTEx0e1qmJn1FUkNb2P2ZSszMyvNwcPMzEpz8DAzs9IcPMzMrDQHDzMzK83Bw6xLxsfHGRoaYsaMGQwNDTE+Pj7lec02WbeH9Z2q6YUvfGGYbYozzzwzBgcHQ1IMDg7GmWee2XHeM888MwYGBopDgsfAwEDD/JOVd7K2zaYXYCIaDUnfKHE6Tg4etikm68A9ODi4Xr7aNDg4OGV5eyUoWW9z8HDwsDrtHOAm68AtqWFeSVOWt1eCkoNMb3PwcPDYLFR92ajfDty9EMDarYPPZvqDg4eDx7Q3GZeNeuXXeT9dOms30EzW91XL70BTjb4JHsChwI2kR4Oe3GD+3wJX5uka0mM5n7Kx9Tp4TH+TcYDrlXaBfmq0b3c/TNaZmgNNtfoieJCeP7wEeCbp2dJXAXu1yP864GftrNvBo3+1+889mQejyTjI94JuBqXJusTmQFOtfgkeBwIXFD5/APhAi/z/AfxlO+t28OhPk3UZpuxBw8pp5yA7WfvWgaZa/RI8Dge+Uvh8NPCFJnkHgLtocckKGAUmgIk5c+ZU+oXa1PA/9/Q2GZfYHGiq1S/B44gGweO0JnnfDHyv3XX7zKO3TMalqDLrtf7jQFPue6hKvwSPti9bAecAR7W7bgeP3jFZ/9xmNQ401emX4LEFcBOwG+sazPdukO9J+ZLVtu2u28Gjd/T6P4ttXqZzoKniDKUvgkeqJ4cBvyHddTU/px0HHFfIcyxwVpn1Onj0Dl+Ksn7VL4Gmyh9dfRM8Jmty8OgdvhRlm4NuBpoq/8ccPBw8Jt1k/LOYbQ66OazOxjh4OHhMql6/Y8Rsumjnf2cqzjyU5k1/w8PDMTEx0e1qTFtDQ0MsW7Zsg/TBwUGWLl069RUy24yNj48zOjrKqlWrHk8bGBhgbGyMkZGRUuuSdHlEDNen+0mCVonly5eXSjezyTMyMsLY2BiDg4NIYnBwcJMCRys+87BK+MzDbHrymYdNqgULFjAwMLBe2sDAAAsWLOhSjcxsMjl4WCWm4jTZzHqHg4e1ND4+ztDQEDNmzGBoaIjx8fGmeUdGRli6dClr165l6dKlDhxm09gW3a6A9a76OzaWLVvG6OgogAOD2WbOZx7W1Pz589e71Q9g1apVzJ8/v0s1MrNe4eBhTfn2WzNrxsHDmpozZ06pdDPbfDh4WFO+/dbMmnHwsKZ8+62ZNeMe5mZm1pR7mNvjyvTdMDNrxP08NjPuu2FmVfCZx2bGfTfMrAo9FzwkHSrpRkmLJZ3cJM9cSVdKulbSRVNdx37mvhtmVoWeCh6SZgKnA68B9gLeImmvujw7AF8E/iwi9gaOmOp69jP33TCzKvRU8AD2BxZHxE0R8ShwFjCvLs9RwNkRsRwgIu6Y4jr2NffdMLMq9FrwmA3cXPi8IqcVPRt4sqQLJV0u6ZhmK5M0KmlC0sTKlSsnobr9x303zKwKvXa3lRqk1XdE2QJ4IfDHwDbApZIui4jfbLBgxBgwBqmfR8V17VsjIyMOFmbWkV4LHiuAXQufdwFuaZDnDxHxIPCgpIuBfYANgoeZmU2OXrtstQjYQ9JukrYCjgTOq8vzXeDlkraQNAAcAFw/xfU0M9us9dSZR0SslnQCcAEwE1gYEddKOi7PPyMirpf0I+BqYC3wlYi4pnu1NjPb/HhsKzMza8pjW5mZWWUcPKYJD3ZoZlOpp9o8bNN4sEMzm2o+85gGPNihmU01B49pwIMdmtlUc/CYBjzYoZlNNQePacCDHZrZVHPwmAY82KGZTTV3EjQzs6bcSdDMzCrj4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWWkdDskv6b+BXeboSuD4i1lZQLzMz62GdPs/jduBPgRMBAY9I+jUpkNSCytUR8VCH5ZiZWQ/p6LJVRBwVEXsCTwJeCawA9gNeC3wWuAS4T9J1ktp6tJ2kQyXdKGmxpJMbzJ8r6V5JV+bp7zvZhl7nJwSaWS+q5EmCEfGgpNcDS4ADIuJuSVsB84APAc8iBZaWJM0ETgdelfMvknReRFxXl/UXEfHaKurey/yEQDPrVVU2mB8LnBoRdwNExKMR8S1gGPg/4OI21rE/sDgiboqIR4GzSAFos+QnBJpZr6oyeKwCnl6fGBGPAacCb2tjHbOBmwufV+S0egdKukrSDyXt3WxlkkYlTUiaWLlyZRvF9xY/IdDMelWVwWMc+KikRgf7rYGd21iHGqTVjxl/BTAYEfsApwHnNltZRIxFxHBEDM+aNauN4nuLnxBoZr2qyuDxIeAm4EZJ/5IbvveS9Abg08Cv21jHCmDXwuddgFuKGSLivoh4IL8/H9hS0o6VbEGP8RMCzaxXVRY8IuJh4E+AjwB/DpxPChjfyVlG21jNImAPSbvlBvcjgfOKGSQ9XZLy+/1J23BnJRvRY/yEQDPrVZP2JEFJewLPIB3Yr4mI1W0udxjpNt+ZwMKIWCDpOICIOEPSCcA7gdXAQ8D7IuKSja3XTxI0Myuv2ZMEKw0ekg4C3g4MAvcC74mInmjddfAwMytv0h9DK2kEuJDUZnEN8DpgxzzvFEmvqKosMzPrrqobzD8TEX8MnMT6d049RhrCxMzMpoEqg8cQ8KP8vv5a2DWkYUvMzGwaqDJ4LAGe12Te3UD/dbQwM7OGqgwenwc+LKnRGcYeQP918TYzs4aq7OcxBlxEGkn3G6RLV0/Pd2DNB35SVVlmZtZdVT9J8HDgb4CDSA3m3yPdgbWSFEDMzGwaqPJW3eOB7SPiCxExB9gbeDXwXODFEXF7VWWZmVl3VXnmcRrwR7UPEXF9RPx3fhbHlpK2qbAsMzProo6Ch6Sn5DGooPGIuDUvAO7rpCwzM+sdnZ55vBt4QNK1pAbyYyS9VtIudfm2AdZ0WJaZmfWITh9DuxD4HbAPsCdwFHA8EJLuBq4EfgO8GKh/lKyZmfWpjoJHRNwM/BuApCNId1stIQWTfYB9gRfiu63MzKaVjoKHpB0i4h6AfIdVzc/yZGZm01CnbR4rJe0DIOkDkl7T5DG0ZmY2jXTa5vFm4I78/r3ATqxr77iqbro2Ih7rsDwzM+sBnbZ5nF14/3RJs1jX3vF84BDgXcCWpGHZn9BJeWZm1hs6PfNYT0SsBH6aJwAkbQHsRQomZmY2DVQ9ttUGImJ1RFwdEWe2k1/SoZJulLRY0skt8r1I0hpJh1dX26kxPj7O0NAQM2bMYGhoiPHx8W5XycyslMrOPHLD+b+Q+nvcQ3oA1FWkvh5X5dt6N7aOmcDpwKuAFcAiSeflIU7q830KuKCq+k+V8fFxRkdHWbVqFQDLli1jdHQUgJGRkW5WzcysbVWeefwbsBvwRdJouk8ktXecByyVdGcb69gfWBwRN0XEo8BZwLwG+d4NfId1jfV9Y/78+Y8HjppVq1Yxf767wZhZ/6iyzWN34M0R8YNioqSdSJ0F22nzmA0Uz1BWAAfUrW828AZSY/yLWq1M0igwCjBnzpxWWafM8uXLS6WbmfWiKs88FgHb1SdGxB0R8eOI+Ewb62g0uGL989A/C7w/IjY6VlZEjEXEcEQMz5rVG0/BbRbEeiW4mZm1o8rg8QHgZElP6mAdK4BdC593AW6pyzMMnCVpKWk4lC9Ken0HZU6pBQsWMDAwsF7awMAACxYs6FKNzMzKqzJ4/BK4FLhG0vvz3VBbl1zHImAPSbvlod6PJLWZPC4idouIoYgYAr4NHB8R53Ze/akxMjLC2NgYg4ODSGJwcJCxsTE3lptZX1FE/VWhTVyR9BXg7cD9pM6AW5GGYf8N6Y6rX7Vz6UrSYaRLUzOBhRGxQNJxABFxRl3erwPfj4hvb2y9w8PDMTExUWKLzMxM0uURMbxBeoXB437gsxHx4Xwr7V6sG1l3X+D5EbFTJYVtAgcPM7PymgWPKu+2uhu4CCA3Zv86T211DjQzs/5RZZvHl4HXV7g+MzPrUVUGjy2B10r6qKQNbtk1M7Ppo8rg8VfAHODDwB2SLpL0OUlvk/QCSVtWWJaZmXVRZW0eEbFzgyHZXwEch4dkNzObVjwku5mZlVblqLqjwOqIWFhMj4jVwNV5MjOzaaDKNo/3AI8PICVprqS7JN0v6RMVlmNmZl1WZfAYAq4ofP4EsAQ4BXinpLdXWJaZmXVRlW0eDwFrASQNkoZSf0lEXCZpLenZHgtbLG9mZn2iyjOPq0hPAAQYAW6PiMvy5yuAPSosy8zMuqjKM48FwPmShoGXAqcV5j0NeLTCsszMrIuq7Ofxc0mHAG8BJoCPFmbPJY2ua2Zm00DV/TwuJT3To95M4JtVlmVmZt1TafBoJiL+airKMTOzqVFlJ8EnA38LPBtYCVxDakS/KiLur6ocMzPrvirPPM4EDgF+Rnr40wiwHRD5eeO/iojDKyzPzMy6pMrg8QrgXcXhSSQ9ixRIXoDHtjIzmzaq7OdxK/D7YkJELI6Ib0fE/Ih4XTsrkXSopBslLZZ0coP58yRdLelKSROSXlZR/c3MrE1VBo8vAW/qZAX52eenA68hjcT7Fkl71WX7b2CfiNgXeDvwlU7KNDOz8qoMHlsDh0j6B0lP3MR17A8sjoibIuJR4CxgXjFDRDwQEZE/bgsEZmY2paoMHicCzwQ+SHqS4MWSPi/p7ZL2k7RVG+uYDdxc+Lwip61H0hsk3QD8gHT2YWZmU6jKHuY7SdqZ1DBem+aSHk/b7pME1WjVDco6BzhH0kHAx4FXNlxZesbIKMCcOXPa2g4zM9u4qnuY30pqOL+glpafXb4X8Lw2VrEC2LXweRfglhblXSxpd0k7RsQfGswfA8YAhoeHfXnLzKwiVV62aigiHouIqyLizDayLwL2kLRbvsx1JHBeMYOkZ0lSfr8fsBVwZ9X1NjOz5io985D0HNJzPB4ArgNuKDRub1RErJZ0AunMZSawMCKulXRcnn8G6Y6uYyQ9RnqGyJvLlGFmZp1TVcddSW8m9TIvtm08TBqm5EpSD/MzKilsEwwPD8fExES3ijcz60uSLo+I4fr0Ki9bfYjU12N7UsP33wLfILV1HA6cWmFZZmbWRVUGj92Bc8mPogUuiojjgWHgLuDFFZZlZmZdVGXwuI90GWxtfv9UgIi4Dvgs8MkKyzIzsy6qMnj8mnT2AXATsF9h3vXAyyssy8zMuqjqsa12zu+/Bpwo6RBJs4HjgHsqLMvMzLqoklt1c7+L70fE2Tnpy8CfAj8l9RBfC/xFFWWZmVn3VdXPYwZwv6R5EfGjPKjhoZIOJPUYvyIiFldUlpmZdVklwSMi1ki6o0H6pcClVZRhZma9o8o2j4XAMRWuz8zMelTVY1u9QtKYpKdWvN6eNz4+ztDQEDNmzGBoaIjx8fFuV8nMbNJUObbVKPA0UsP42yRNABPkoUmAa3JbyLQzPj7O6Ogoq1atAmDZsmWMjo4CMDIy0s2qmZlNisrGtgKQNAvYJ0/Pz697kp/nEREbe57HpJnMsa2GhoZYtmzZBumDg4MsXbp0Uso0M5sKzca2quzMQ9LxwHhE/JR0i24tfQvS8zyeX1VZvWb58uWl0s3M+l2VbR6nAX9UnxgRq4EbgO9UWFZPafaUQj+90Mymq46Ch6SnFJ5N3ugRsjUvII13NS0tWLCAgYGB9dIGBgZYsGBBl2pkZja5Oj3zeDfwgKRrST3Jj5H0Wkm71OXbBljTYVk9a2RkhLGxMQYHB5HE4OAgY2Njbiw3s2mrowZzSbsCB5Maxk8ijV+1AymQ3E260+o35OHYI2K/BquZEn4YlJlZeZPSYB4RNwP/lgs4gvTQpyWsu+NqX+CFwEpgfidlmZlZ76jsbquIKLYO/yxPZmY2DW20zUPSRO41fpykF0ma1L4akg6VdKOkxZJObjB/RNLVebpE0j6TWR8zM9tQO2ce+5HulnpH/rxG0vXAFcDl+fXKiFjVaWUkzQROB14FrAAWSTovP42w5nfAKyLibkmvAcaAAzot28zM2tdO8HgHKYDsR+roty3wvDzVBkJcK+m3rAsmVwC/ioiyt+fuDyyOiJsAJJ0FzAMeDx4RcUkh/2VA/Z1dZmY2yTYaPCLia6QnA9Ye+vQc1gWT/UiN4k/K6c8BjqotK2lJRDy7RH1mAzcXPq+g9VnFO4AfNpspaZQ05pY77JmZVahUg3mk+3qvz9Pjw8ZK2p31A8oLgB1Z90zzdjXqaNjwXmJJB5OCx8ta1HeMdFmL4eHh6gbxMjPbzFX1MKglpFt0v1VLy31AXlByVStITx6s2QW4pT6TpOcDXwFeExF3lq6wmZl1pMoh2deT+4DcvNGM61sE7CFpN+D3wJEULoMBSJoDnA0cHRG/qaKuZmZWzqQFj00REaslnQBcAMwEFkbEtZKOy/PPAP4eeCrwxdQEw+pGvR/NzGzyVPo8j17m4UnMzMprNjxJ1Y+hNTOzzYCDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBhZmalOXiYmVlpDh5mZlaag4eZmZXm4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV1nPBQ9Khkm6UtFjSyQ3mP0fSpZIekfQ33aijmdnmbotuV6BI0kzgdOBVwApgkaTzIuK6Qra7gPcAr5/6GpqZGfTemcf+wOKIuCkiHgXOAuYVM0TEHRGxCHisGxU0M7PeCx6zgZsLn1fktE0iaVTShKSJlStXdlw5MzNLei14qEFabOrKImIsIoYjYnjWrFkdVMvMzIp6LXisAHYtfN4FuKVLdTEzsyZ6LXgsAvaQtJukrYAjgfO6XCczM6vTU3dbRcRqSScAFwAzgYURca2k4/L8MyQ9HZgAtgfWSjoR2Csi7utWvc3MNjc9FTwAIuJ84Py6tDMK728jXc4yM7Mu6bXLVmZm1gccPMzMrDQHDzMzK83Bw8zMSnPwMDOz0hw8zMysNAcPMzMrzcHDzMxKc/AwM7PSHDzMzKw0Bw8zMyvNwcPMzEpz8DAzs9IcPMzMrDQHDzMzK83Bw8zMSnPwMDOz0hw8zMysNAcPMzMrreeCh6RDJd0oabGkkxvMl6TP5/lXS9pvsuoyPj7O0NAQM2bMYGhoiPHx8ckqysysr2zR7QoUSZoJnA68ClgBLJJ0XkRcV8j2GmCPPB0AfCm/Vmp8fJzR0VFWrVoFwLJlyxgdHQVgZGSk6uLMzPpKr5157A8sjoibIuJR4CxgXl2eecC/RXIZsIOknauuyPz58x8PHDWrVq1i/vz5VRdlZtZ3ei14zAZuLnxekdPK5gFA0qikCUkTK1euLFWR5cuXl0o3M9uc9FrwUIO02IQ8KTFiLCKGI2J41qxZpSoyZ86cUulmZpuTXgseK4BdC593AW7ZhDwdW7BgAQMDA+ulDQwMsGDBgqqLMjPrO70WPBYBe0jaTdJWwJHAeXV5zgOOyXddvRi4NyJurboiIyMjjI2NMTg4iCQGBwcZGxtzY7mZGT12t1VErJZ0AnABMBNYGBHXSjouzz8DOB84DFgMrALeNln1GRkZcbAwM2ugp4IHQEScTwoQxbQzCu8DeNdU18vMzNbptctWZmbWBxw8zMysNAcPMzMrzcHDzMxKU2p/nv4krQSWdbseHdgR+EO3KzEJvF39Z7pum7erscGI2KCX9WYTPPqdpImIGO52Parm7eo/03XbvF3l+LKVmZmV5uBhZmalOXj0j7FuV2CSeLv6z3TdNm9XCW7zMDOz0nzmYWZmpTl4mJlZaQ4ePU7SUkm/lnSlpIlu16cTkhZKukPSNYW0p0j6iaTf5tcnd7OOm6LJdp0i6fd5v10p6bBu1nFTSNpV0s8lXS/pWknvzel9vc9abNd02GdbS/qlpKvytn00p1e+z9zm0eMkLQWGI6LvOy9JOgh4gPQM+ufmtE8Dd0XEJyWdDDw5It7fzXqW1WS7TgEeiIjPdLNunZC0M7BzRFwh6YnA5cDrgWPp433WYrv+nP7fZwK2jYgHJG0J/A/wXuCNVLzPfOZhUyYiLgbuqkueB3wjv/8G6Z+4rzTZrr4XEbdGxBX5/f3A9cBs+nyftdiuvhfJA/njlnkKJmGfOXj0vgB+LOlySaPdrswkeFrtSZD5dacu16dKJ0i6Ol/W6qtLO/UkDQEvAP6PabTP6rYLpsE+kzRT0pXAHcBPImJS9pmDR+97aUTsB7wGeFe+RGK970vA7sC+wK3AqV2tTQckbQd8BzgxIu7rdn2q0mC7psU+i4g1EbEvsAuwv6TnTkY5Dh49LiJuya93AOcA+3e3RpW7PV+Drl2LvqPL9alERNye/4nXAv9Kn+63fN38O8B4RJydk/t+nzXarumyz2oi4h7gQuBQJmGfOXj0MEnb5gY9JG0LvBq4pvVSfec84K35/VuB73axLpWp/aNmb6AP91tufP0qcH1E/HNhVl/vs2bbNU322SxJO+T32wCvBG5gEvaZ77bqYZKeSTrbgPS8+f+IiAVdrFJHJP0nMJc0RPTtwEeAc4FvAnOA5cAREdFXjc9Ntmsu6fJHAEuBv6pdc+4Xkl4G/AL4NbA2J3+Q1D7Qt/usxXa9hf7fZ88nNYjPJJ0cfDMiPibpqVS8zxw8zMysNF+2MjOz0hw8zMysNAcPMzMrzcHDzMxKc/AwM7PSHDys50jaTVJIWlhIOzCnfaSbdWtHHgk58tAXU1VmSPKtkzZlHDysF704v15aSHtJgzQz65Itul0BswZqweOSQtqBpM5b/7dhdgP27HYFbPPi4GG96MXAvcB1hbQDgesi4t7uVKm3RcQN3a6DbV582cp6iqStSENEXBZ5+ANJg8AzgMs2cZ2n5DaBU3J7ypmSbpf0cH7a2l9LavhDSsnRki6UdHdeZomk0yXtWrIeg5I+kJ9id7OkRyTdlT8f1WSZubnuF0oakPQPkm6Q9FAedruWb4M2D0lfr6W3mC5sUOaukj4n6cZczn2S/lfSsXlcqPr8F+Z1zZV0kKQfSPqDpLWSXp/zbC3pZElXSHogb/utki7N27R1me/Sus9nHtZ1TRp6/6RB+jskvSO/XxYRQyWL2g2YAB4mjTa6PXAw8BngZZLelEdUrdVLwJnAUcBjeZm7SKOtHg8cKenQiFjUZvlHAx8HlpAGq/tf0rDZLwfmSjogIt7bZNmtc/l7AhcDVwFbbaS8/2kx72Wk4cfXFBMlHUwaT+1JwGLgR8B2pLPBrwGHAMc0WecRwHGkM8afkMb6ekzSDOAHedl7gYvy69OAPwLmA18AbtvI9lgviQhPnro6AV8pTNeR2ja+XUhbntO+Vkj7ZIn1n5KXr61368K8PYAVed7xdcsdn9NvA/YupM8EPs+6AfSeULfc0jxvqC79RcX11NWhto0H1M2bW6j7r0gP9Wm0jZH+ndv6PuYCjwCriuUBO5OC42rSyKsqzNs1lx/AsXXru7BQx9EG5R2U511OekRqcZ6AlwID3f479FRu6noFPHkqTqRf4w8AWxXSlgM3drDOWvBYBezUYP7b8vzf1qUvyel/2WCZrQoH/JG6eQ2Dx0bq+Jd5mX+qSy8Gj5e0WL6t4EE6c7mbdMbxhrp5n8rr+VSTZYdrQaAuvRY8ftxkuSPy/M92++/LU3WTL1tZz5C0PemS0E8i4tGctgfpV+8ZFRTx40gP1ao3TjqbeZak2RHxe0m7AM8kDdn97/ULRMSjksaBk0kH+PF2KpCv7f8J6SxkFvCEPKv2LIlnN1n09oi4pMm8tkh6GnA+sANwUkScU5flsPz6rSaruJwU2PeVtHVEPFw3/+wGywBcQQpW75D0G+A7EXF72fpbb3HwsF4yl/Q3+dNC2iH59WcVrP93jRJzILgVmE1qg/h9fg9wa4ODZM2S/Dq7yfz1SDqQ9EyFXVpk275J+rJ2ymhR9gDwPWAI+HxEfLZBtmfm10UN2sXrPZX0PRU1rGNELJF0Eqlt6XTgdEk3kW7F/i5wTkSsabSs9S4HD+saSXOBnzeYdaqk+udHf7NwQPtGRBw7SdWqNdKr7nMjGz3CPp4xHbzPITUSf5X0vOzFwP0RsVbSq4ELWqzzoXbLalD2DOA/SGc73wVOapJ1Zn79L9JNBa080iCtaR0j4jRJ3wJeT2qsfxnw//J0paRXxDR6PvrmwMHDuuk20lPPag4nNdaemz+LdIfSbcCPC/la3UXUylCjxHx7cO2y0S35dUV+fYakJ0REo4Plbvm1/hd4IweRAsflEfEXDeY/q411bKp/AeYBvwSOisIdZXVuzvX4eERcW3UlIuI20uXHMwAk7UO6JLgv6fLfB6su0yaPg4d1TaSObccCKD13+WjS9fBa2vNIt4WeFRHvq6DIV0uaFREr69LfQurztCQiVuS6rciXVp5J+nX81eICkrYk3cILqcF4Y56SX29uMr9hP49OSToReA/pkt3rImJVi+w/BN5NauCuPHjUi4irJH2O1N60z2SXZ9VyJ0HrFXNJf4/Fy1gH59cq2jsABoAvSKo1UiNpd1LfC4DP1eX/5/z6cUnPKSwzE/g0MEi6zv/tNsqu9QA/pG5dMyT9Pel21UrlDnqnku6uOqzJzQJF/wTcB3xQ0rsadZyU9GJJR5SsxyGSDqtfX/4ea430HbXp2NTzmYf1ilrDeH3wWEPqFFeFfwf+FFgi6X9Jnd8OIXXA+x6pMbfoi6SD+luAqyT9nHQg3p90RnI3cESTS1rriYgrJH0PeB3pGv/PSR3lXgTMIQWjv+t4C9d3KikgLwNObtIIfkNEfDLX8eYccL5N6rQ3X9K1wJ2kHv6759f/ovkdWY08n3Tp7F5JVwC3kgL5AaTLhbeRbhO2PuLgYb3iYOC2fCmr1rv7IGCiwobUm0gH638kBY0n5bSFpD4I67UFRERIGiFdzvlLUi/rbUjtIl8CPhERzS5DNXI4qbH6aNKZ1gOkUYKPyuutOnjUGsD3zVMjFwGfrH2IiJ9L2pt0+epPSdu8JekA/1vgNMoFDkiBeQfS/nwWaYTkB0j9ZM4AvtTgUqL1OEX4EQA2vUk6BfgI8NGIOKW7tTGbHtzmYWZmpTl4mJlZaQ4eZmZWmts8zMysNJ95mJlZaQ4eZmZWmoOHmZmV5uBhZmalOXiYmVlp/x9C86qBZKpfNAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the result\n",
    "plt.figure()\n",
    "plt.plot(Narray, I, 'ko')\n",
    "plt.xlabel('# polarizers', fontsize=22)\n",
    "plt.ylabel('$I_{transmited}$', fontsize=22)\n",
    "plt.title('Imperfect polarizers', fontsize=24)\n",
    "\n",
    "# Find the maximum\n",
    "Imax = np.max(I)\n",
    "Nmax = np.argmax(I) + 2\n",
    "print('The maximum obtained intensity is {} for {} polarizers.'.format(Imax, Nmax))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the increase in transmission expected by the increase of the number of polarizers does not compensate the extra loss produced by the imperfect polarizer transmission. This produces a maximum of the transmitted intensity (using 17 polarizers for this data).\n",
    "\n",
    "## Mueller-Stokes formalism\n",
    "\n",
    "Mueller-Stokes formalism is more complicated than Jones formalism, but allows working with partial polarization states. This example uses completely polarized light, but it is a good example to start using *Mueller* and *Stokes* classes.\n",
    "\n",
    "### Perfect polarizers\n",
    "\n",
    "This is the easiest case. Start by importing all the required modules and classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from py_pol import degrees\n",
    "from py_pol.stokes import Stokes\n",
    "from py_pol.mueller import Mueller"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will work using linear polarization with the same azimuth as the first polarizer (0º)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Incident wave = \n",
      "[+1.000]   \n",
      "[+1.000]   \n",
      "[+0.000]   \n",
      "[+0.000]   \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Incident wave\n",
    "S0 = Stokes('Incident wave')\n",
    "S0.linear_light(intensity=1, azimuth=0*degrees)\n",
    "print(S0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will construct two arrays, one with the number of polarizers and other with the output intensity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nmax = 30\n",
    "Narray = np.arange(2,Nmax+1)\n",
    "J = Mueller('System')\n",
    "I = np.zeros(Nmax-1)\n",
    "\n",
    "# For loop to calculate the intensities\n",
    "for ind, N in enumerate(Narray):\n",
    "    # Create the system of polarizers\n",
    "    angles = np.linspace(90*degrees, 0*degrees, N) # Start by 90º and end in 0º because the rightest  element is the first to be crossed by the light wave\n",
    "    J.diattenuator_perfect(azimuth=angles)\n",
    "    # Multiply them\n",
    "    J.prod()\n",
    "    # Multiply the system matrix by the incident light wave to obtain the output\n",
    "    S_final = J * S0\n",
    "    # Store the intensity\n",
    "    I[ind] = S_final.parameters.intensity(out_number=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, plot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEpCAYAAABr364UAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApqUlEQVR4nO3deZhcVZ3/8fcnAYRmVQmIgXSjooIKCJHFBQLjAqgDjjIK/QNhnGkRQXGZEY2OQY2KIzPuMq0TUekZ3ADRQRRHEGcQJSD7ogmkk8gWZadlSfL9/XFOkZtKVXff7ttdVd2f1/Pcp6rOXc65dbvvt+4595yriMDMzKyMGa0ugJmZdR4HDzMzK83Bw8zMSnPwMDOz0hw8zMysNAcPMzMrzcHD2oKSkyRdI2lIUuSpp9Vl6zSt/O4kXZrzPm6y87bJ5eAxjUg6q3BiKU4P5pP2v0jasUXF+xDwRWAPQMDdeVoz2QWRtI2kBZIWTHbeZp3CwWN6eoJ1J+d7gC1IJ+33A9dLenkLyvTu/PpeoCsinpGnFS0oyzbAR/Nk5SwHbgUeaHVBbGJt1OoCWEtcHhHzah8kdQFvBL5AOnF+T9KzIuIvk1EYSdsBs/LHr4WHPehYEXFsq8tgk8NXHkZEDEXEt4F35aRnAEdMYhE2K5Tl4UnM18zGyMHDir4LrM3v9y7OkDRD0jGSLpa0StLjku6Q9B1J+zbaWG43iNzWMiM3iP9W0v05/RRJASwrrFNsi1kw3jIU1t1c0vslXS7pXkmPSrpN0gWSeiVtnJe7FLi9SXk2KNMw+fXU1smfXybpx7ncQ7mN6SRJw/4PSvobSRfl9R6TtFLSgKS9RlOOum3NlHSQpM9LukrS3YXv8DxJBw+z7pMN4blN6HRJt+R9ub/Rco2+i1FMPU2+yy9KujXn91Au/wckbd6kvE9uT9Kukr4paYWkJySdX1huu9zWd4OkR/LfxYr8d/IxSd1lv+dpIyI8TZMJOAsI4NJhlrk7L9NfSNsSuDinBynAPFD4vAY4qcG2FuT53wTOz+9XA/fl9y8F7gJWFbZ1V2F6/3jLkNfdjRQQass+UShDberJy547THnWK9MI33VPYRtvzHlGzveJwrzzgI0arD8jf2+15YrfW21/39Ek7/X2qZD+wrp9fhR4uC7tQ022eWme/4/A0sL6DwL3N1juuELaTg2+x+J0/zBl/hvgL4X5Q8Bjhc/XAdsP8x0cAzyS3z+Yt3V+XqYbuKPuO743/23V0k5o9f9tu04tL4CnSTzYIwQPUvVR7R/nM4X083LatcBhwGY5fRvgg/mfeQ3wsrrtLcjrPZRPNO8gNYYDbAdsld/31P5Zhyn7WMvwNFIjbgC3AYcDm+R5WwIvBxYBOxbWGbE8o/iun9xGPjn+BNg5z9ucdBJeQ5MTNnAq64Lkh4Etc/ps0hViLYAc0GDdZifi5+Z1XwdsD6hwLD6cT55rgX0bbPPSwrFcDhwCzMjzntNgueNG+T1tDFyW17kG2LQw7yXA47lcnwbmkO7EmwnsC1yR1/vpMN/BQ7lML8zpAp6d3y/Ky/wBeEVhf55CCrQfB45o9f9tu04tL4CnSTzYIwePkwr/dG/Maa/Mn28HntZkvX/Ky/y4Ln1BYXt9w5Rr2JP1OMvwmZy+Cpg9yu9p2PKU3QZwA/CUBsvUvp8HgM0L6Zuz7qrqUw3Wmwn8Ks+/rMH8hsFjFGX+SF7vGw3m1YLC47UTcZNt1JY7bpR5npmXvwforpv3v3nee5qs+1Tgj3mZuU2+g6XkHxoN1r8pL/PmsR7n6Ty5zWOaU9Ij6f2kEy3AIPCj/P6t+fWsiLi3yWb+M78eJGlmg/l/Jv3KG6vxlOGY/PrZiPjjOMowHmdExGMN0v+VdEW2FfCqQvqrc9rjrDsmT4qINaRfxQCvkPSMispZO+YvG2aZn0TEDVVkJukdwNtJ1XhviojBwrxn53L8hRRgNhAR95Gu6GD976/oS9H8rsEH8+sOJYtu+Fbd6erAWkNuA3eSLtUfz59fml/fk//Zh9MFPJ30K7JocUSsHltRx16G3PhaO7FeOI78x+vSRokR8aCk3wH7A3uR2oXI7wGuzSfIRi4jVedslJcf1f5J2gw4gVR9txvp13v9eeCZw2zi16PJZxTlOAD4fP54ckRcVrdI7ZhvAtwuqdmmtsivOzWZP1x5LyRVf50uaRfg+8AVwwQbK3DwmJ6eIDUMQrpsf4TUHnAx8PW6E1btV9nWeRpJV4O0VWMs53jLsH0hbfk4yzAew13x1ObNKqTNqpu3gYh4VNKfSfs4q9lyRZJ2IAWy5xaSHyE1xK8lVYdtS6o2a2a8xxJJc0gn6o2Br0bEvzdYrHbMZ7L+cWym0d8dDF/e00l3Ff41cGKeVku6ktTG9rWIuH8UeU9LDh7T03qdBEdQq9o8PCIuGGN+4x1iZKxlaPpztY0MV8anVJzX50iB4zZSg/0lxR8KuapoyQjbGNexVOqQ+kNSwPsl60YWqFc75r+LiNK3JRc0LW+uSjxc0n7AG4ADSMFk/zz9o6RXRcS148h/ynKbh43k7vy6WweW4a7C++6KyjIWw1UD1X5hF38h1943LbOkTUnVc/XrNlt+E1JVFUBvRJzboEpsNL/wx+sbwJ6kvj1viognmixXO+a7SJrQH7kRcUVEfCAi9idV4x1FulKdBXx9IvPuZA4eNpJanfEbO60MEbGMdQHksBKr1jpKomEq20s4sFGipC1Z175xdWFW7f0ukmY32eYBrKs5uLrJMkXbsu5K5ndNlnnlKLYzZpLmA39Lqio7PCL+NMzitWO+BekGgkkREY9ExDlAX07au1lHxOnOwcNGclZ+nStp2HGLJD21Dcvw7fz6vmFOxPUeLLzfZpTrDOd9+Zd/vVOATXN+Pyuk/yynbUyqXlpPvpvsI/njryLirvplGniQ1L4F8KIG29wBOHkU2xkTSa8HPpbLcGxEXDfc8hFxC6kfB6QG7aYncEmbSSpdxdfkmNTUGs1FarS3Og4eNqyIuIjU6xpgkaTT8okGSCdrSYdL+iHp1tN2K8PppIbnbYFfSfrr2klD0haS5kk6R4Wh6HMj6R354/EV7MIc4Lza0BuSuiS9l3Wj9p4eEUOF/B8BPpk/vkvSfElb5HVnA/9F6txY60A4okhjhtVOxosk7Zm3N0PSX5HaHyakjUjSc4GzSeeb0yLi3BFWqTmZ1PnzhaRj98paFVYu9wskfZjUl2Mst9veIOmTkl5S+JuQpH1IjwcAuHKYO96mt1Z3NPE0eROjGJ6kyXqbs66Hd7HX9AN1ad+oW29BTj9rhO331LZRdRnyui8CVhSWeZwmw5MU1jmtMO9hUh39MuCUUX5nPYX1hxue5HwaD08ykw2HJykOnbEGOLFJ3s32aV/S8B7F/ap9/jOpTaThcWCUnf8aLQccV8jzboYfqmSnuu0dyvrDlzwG/Ckfw+Lx6x7Nd1C3THG7q/N3UNzuKmD3Vv/ftuvku61sRJF+Cb9B0muBvyOdhGaRTmRLgN+SrgwmrC/FeMoQEddLegHpl+wRpDuONiXddXQd8B1gZd1qHyPVzfcCz2Fd4/U2Yyj7DyQdRBpyZD/SiepGUmPsVyJibYN11gBvlXQBqf59b1LHwTtJVwlnRMRVJcvxG0n7k4L6gaSAfCdwEbCQFLAm2nYjzF+vDBHxk3zlcjKp3eo5pGNwP+m5IRcB34tCB8MSDgdeQ2o/mkO6YeBx4GbS39G/RUR9nyXLamPbmFmFchXV7QAR0Qm3DJuV4jYPMzMrzcHDzMxKc/AwM7PSHDzMzKy0adNgvu2220ZPT0+ri2Fm1lGuuuqqP0XEBoNvTptbdXt6eli8eHGri2Fm1lEkNbwN2tVWZmZWmoOHmZmV5uBhZmalOXiYmVlpDh5mZlaag4eZ2RQ0MDBAT08PM2bMoKenh4GBgUq37+BhZtZiZU70o1l2YGCAvr4+BgcHiQgGBwfp6+urNoC0ekz4yZr23nvvMDObLGeffXZ0d3eHpOju7o6zzz676XJdXV3FZ5NEV1dXw+VHu2x3d/d6y9Sm7u7u0vsBLI4G59SWn9Qna3LwMLPxmoiAUOZEP9plJTVcTlLpfXbwcPAwswZaHRDKnOhHu6yvPBw8zGwMOikgTMSVR5n9GomDh4OHWUebqgFhIto8ynxfI3HwcPAwazsOCOW+h7LLVsHBw8HDbNKM5gTngNAZHDwcPMzGpeqrBAeEzuDg4eBhtoFWVhs5IHQGBw8HD7P1tLrayAGhMzh4OHjYNDHaE2erq40cEDqDg4eDh3WwiaheaodqIweE9ufg4eBhHarVQ12ULUNteQeFqcHBw8HD2kyrq5ccEGw0HDwcPKyNtEP1Uq0cDgg2nGbBQ2ne1Dd37txYvHhxq4thBkBPTw+Dg4MbpHd3d7Ns2bIxL1t7jsPQ0NCTaV1dXfT399Pb21tJ2W16kXRVRMytT/fDoMwqNNqH+ixfvnzU6QsXLqSrq2u9tK6uLhYuXLjBsr29vfT399Pd3Y0kuru7HThsYjS6HJmKk6utbKJNVMN2bduuXrJWoEm1la88zEYw2quJ+fPnr1ddBDA0NMT8+fM3WLbM1QSkK4ply5axdu1ali1b5isJazkHD7NhlHkWdJmqKFcvWadzg7nZMCaqYdusU7jB3KzOaKqjJqph26zTOXjYtDTa6qg5c+Y0XL9RuquibDpxtZVNS6OtYnK/CZvuOqbaStIhkm6VtETSqQ3mby3pR5KulXSjpONbUU5rT1X3s/DVhFljbXXlIWkm8HvgVcBK4ErgqIi4qbDMh4CtI+IDkmYBtwLPiIjHh9u2rzymvjJXCW7cNhudTrny2AdYEhG35WBwDnB43TIBbClJwBbAvcDqyS2mtaOJ7GdhZutrt+AxG1hR+LwypxV9CdgVuAO4Hnh3RKxttDFJfZIWS1q8atWqiSivtRH3szCbPO0WPNQgrb5e7TXANcAzgT2BL0naqtHGIqI/IuZGxNxZs2ZVWU6bRKNtxyhzZxS417bZeLRb8FgJ7FT4vCPpCqPoeODcPOzKEuB24PmTVD6bZGV6eLsqymzytFvwuBLYRdLOkjYB3gJcULfMcuCvACRtDzwPuG1SS2mTpkw7hquizCZPW91tBSDpMOBzwExgUUQslHQCQEScKemZwFnADqRqrk9HxNkjbdd3W3WmGTNm0OhvVBJr1zZs6jKzCjW722qjVhRmOBFxIXBhXdqZhfd3AK+e7HJZa8yZM6fhLbXN2jHMbHK0W7WVTROjbQR3O4ZZe3LwsElXphHc7Rhm7ant2jwmits82od7d5t1jk7pYW7TQJnOfGbWnhw8bNKV7cxnZu3HwcMq40Zws+nDwcMq4UZws+nFDeZWCTeCm01NbjC3CeVGcLPpxcHDKuFGcLPpxcHDKuFGcLPpxcHDKuFGcLPpxcHDhjXa22/BD1cym07ablRdax+1229rz9Oo3X4LODCYTXO+8rCmyjyIycymFwcPa8q335pZMw4e1pRvvzWzZhw8rCnffmtmzTh4WFO+/dbMmvHYVmZm1pTHtjIzs8o4eExDZTr+mZk14k6C04w7/plZFXzlMc2445+ZVWFcVx6S/gf4XZ6uAW6OiLUVlMsmiDv+mVkVxlttdTfwWuAUQMBjkq4nBZJaULkuIv4yznysInPmzGn4xD93/DOzMsZVbRURR0fErsDWwCuBlcBewOuAzwGXAw9KukmSW2XbgDv+mVkVKmnziIhHgCOApcB2ETEb2Ap4C3AT0A3MqiIvGx93/DOzKlTWSVDSA8AbI+LndekbAz8FfhERn6gkszFwJ0Ezs/Imo5PgEPCM+sSIeAI4Azi+wrzMzKyFqgweA8BpkmY3mLcpsEOFeZmZWQtV2Unww8AewK2SvkaqqloOPA/4DHB9hXmZmVkLVRY8IuJRSa8B3gO8F3g3EKRbeG8H+qrKy8zMWqvSHuYRsTYizsh3W70AeDXp1t3nRcS1VeZl6/N4VWY2mSoNHpIOkHSWpEuATwF/iIhrImJ1iW0cIulWSUskndpkmXmSrpF0o6RfVlX+TlUbr2pwcJCIeHK8KgcQM5solQUPSb3ApcBOwA3A64Ft87wFkg4cxTZmAl8GDgV2A46StFvdMtsAXwH+OiJeABxZ1T50Ko9XZWaTrcorjw8Dn42IvyK1e6gw7wnSECYj2QdYEhG3RcTjwDnA4XXLHA2cGxHLASLinvEWvNN5vCozm2xVBo8e4KL8vr7n4Q2kto+RzAZWFD6vzGlFzwWeKulSSVdJOrbZxiT1SVosafGqVatGkX1najYulcerMrOJUmXwWAq8qMm8+xjd8CRqkFYfiDYC9iYNyPga4COSnttoYxHRHxFzI2LurFlTd3QUj1dlZpOtyuDxBdKJvNEVxi7AaH76ryS1mdTsCNzRYJmLIuKRiPgTcBmpf8m05fGqzGyyVdnPoz/387gc+D7piuEZkg4A5gMXj2IzVwK7SNoZ+CNpYMWj65b5IfAlSRsBmwD7Av9WzV50rt7eXgcLM5s0VT+G9k3AO4F/IlVB/Si/XkkKIMOKiNWSTiL1Tp8JLIqIGyWdkOefGRE3S7oIuA5YC3w9Im6oeD/MzGwYVY6qeyIwEBEP5M+7As8E7iQ9YbCajMbIo+qamZU3GaPqfpE0jhUAEXFzRPxPRNwEbCxpswrzMjOzFhpX8JD0NEmb1D4Os+iLgQfHk5eZmbWP8V55nAw8LOlGUgP5sZJeJ2nHuuU2A9aMMy8zM2sT420wX0QaMXcPYFfSnVEnAiHpPuAa4PfAfqTH0ZqZ2RQwruARESuAbwFIOpJ0t9VSUjDZA9iT1KFvFaO428rMzDrDuIKHpG0i4n6AiCiOhfGLPJmZ2RQ03jaPVZL2AJD0QUmHNnkMrY2Rn9NhZu1ovG0ebwZqo9q+G9iOde0d19ZNN0bEE+PMb1qpPaejNtx67TkdgHuTm1lLVdZJEEDSLNa1d+zOuob0jYEnIuIplWVWUid2Euzp6WFwcHCD9O7ubpYtWzb5BTKzaadZJ8FKhyeJiFXAz/NUy3gj0oOddq8yr+nAz+kws3ZV6WNoG4mI1RFxXUScPdF5TTV+ToeZtasqH0O7h6RfSLpT0s2Svifpw7nT4E4jb8Hq+TkdZtauqrzy+BawM+n54j8CtiSNsHsBsEzSnyvMa1rwczrMrF1VOaruw8CbI+K/69K3I3UW3D0iPltJZmPQiQ3mZmatNhkN5lcCW9QnRsQ9wM/yZGZmU0CV1VYfBE6VtHWF2zQzszZUZfD4LfBr4AZJH5D0EkmbVrh9MzNrE1UGj37gBGAr4DTgN8BDkm6UNCDp/RXmZWZmLVRl8HgzsDAitgY2J/UuPx74CbA96bnmZmY2BVTZYH4f8EuAiFgDXJ8ndw40M5tiqrzy+HfgiAq3Z2ZmbarK4LEx8DpJp0na4JZdMzObOqoMHm8H5gAfAe6R9EtJn5d0vKQXS9q4wrzMzKyFKmvziIgdGgzJfiDpDqyNgSeAlg3JbmZm1fGQ7GZmVlplwUNSH7A6IhYV0yNiNXBdnszMbAqoss3jXcCs2gdJ8yTdK+khSZ+qMB8zM2uxKoNHD3B14fOngKXAAuAdkv6uwrzMzKyFqmzz+AuwFkBSN7Av8NKIuELSWtKzPRYNs76ZmXWIKq88rgVeld/3AndHxBX589XALhXm1dEGBgbo6elhxowZ9PT0MDAw0OoimZmVUuWVx0LgQklzgZcBXyzM2x54vMK8OtbAwAB9fX0MDQ0BMDg4SF9fH4CfEGhmHaOyK4+IuAQ4GLgJ+Dzw0cLsecDvq8qrk82fP//JwFEzNDTE/PnzW1QiM7Pyqu7n8WvSMz3qzQS+W2VenWr58uWl0s3M2lGlwaOZiHj7ZOTTCebMmcPg4GDDdDOzTlFZtZWkp0r6pKTvS/qqpHdKermkLUtu5xBJt0paIunUYZZ7iaQ1kt40/tJPnoULF9LV1bVeWldXFwsXLmxRiczMyqvybquzgfcAmwF7kvp5XAbcL2mppO+PtAFJM4EvA4eShjQ5StJuTZY7HfhpZaWfJL29vfT399Pd3Y0kuru76e/vd2O5mXWUKqutDgTeWRyeRNJzSIHkxYxubKt9gCURcVte/xzgcFIjfNHJwA+Al4y/2JOvt7fXwcLMOlqVweNO4I/FhIhYAiwBRrzqyGYDKwqfV5I6Gz5J0mzgDaQ7u4YNHnm8rT5wm4KZWZWqrLb6KvDGcW5DDdKi7vPngA/kR90OKyL6I2JuRMydNWvWSIubmdkoVXnlsSlwsKRPAKdHxENj2MZKYKfC5x2BO+qWmQucIwlgW+AwSasj4vwx5GdmZmNQZfA4hXQy/xDwPklXAtcUphsiYqRe5lcCu0jamVQF9hbg6OICEbFz7b2ks4AfO3CYmU2uKp8kuJ2kHUgN47VpHunxtKN6kmBErJZ0EukuqpnAooi4UdIJef6ZVZXXzMzGruoe5neSGs6fvIU2P7t8N+BFo9zGhcCFdWkNg0ZEHDfWspqZ2dhNeA/ziHiCNOLutROdl5mZTY5Kg4ek55NurX2Y1Dfjloiov1vKzMw6XJXPMH8zqZd5sW3jUUk3kBrMf+c2CzOzqaHKfh4fJvX12IrUX+MfgW+S2jreBJxRYV5mZtZCVQaPZwPnkx9FC/wyIk4k9cu4F9ivwrzMzKyFqgweDwKKiLX5/dMBIuImUq/wT1eYl5mZtVCVweN60tUHwG3AXoV5NwOvqDAvMzNroarHttohv/8GcIqkg/NAhicA91eYl5mZtVAld1spDTT144g4Nyf9O/Ba4OekgQ3XAn9fRV5mZtZ6Vd2qOwN4SNLhEXFRHsPqEEn7kwY6vDoPz25mZlNAJcEjItZIuqdB+q+BX1eRh5mZtY8q2zwWAcdWuD0zM2tTVQYPgAMl9Ut6esXbNTOzNlLl2FZ9wPakhvHjJS0GFpOHJmF0z/MwM7MOUOXzPHaQNAvYI0+7k/p29DHK53mYmVlnqHJgxBOBgYj4OekW3Vr6RqTneexeVV5mZtZaVbZ5fBF4Xn1iRKwGbgF+UGFeZmbWQuMKHpKeJmmT2sdhFn0xabwrMzObAsZ75XEy8LCkG0k9yY+V9DpJO9YttxmwZpx5mZlZmxhvm8ci4HZSA/muwNHAiUBIuo90p9XvScOx3zTOvMzMrE2MK3hExArgWwCSjiQ99Gkp6+642hPYG1gFzB9PXmZm1j6qvFV3TuHjL/JkZmZT0IjBI3f2uzpPVwHXRcRjE10wMzNrX6NpMN8LeBvwZeAK0ui510r6hqSTJL1UUteElrIDDAwM0NPTw4wZM+jp6WFgYKDVRTIzmzCjqbZ6GymA7EXq6Lc58KI81QZCXCvpD6Qrk9pVyu8iYlrcnjswMEBfXx9DQ0MADA4O0tfXB0Bvb28ri2ZmNiEUEaNfOD306fmsCyZ7kRrFty4sVtzg0oh47viLOX5z586NxYsXT8i2e3p6GBwc3CC9u7ubZcuWTUieZmaTQdJVETG3Pr1Ug3mkSHNznp6sl5H0bNYPKC8GtmXdM82ntOXLl5dKNzPrdFU9DGop6Rbd79XSJO1ECiJT3pw5cxpeecyZM6fB0mZmna/q53k8KSJWRMQFE7X9drJw4UK6uta/Z6Crq4uFCxe2qERmZhNrwoLHdNLb20t/fz/d3d1Ioru7m/7+fjeWm9mUVarBvJNNZIO5mdlU1azB3FceZmZWmoOHmZmV5uBhZmaltV3wkHSIpFslLZF0aoP5vZKuy9PlkvZoRTnNzKaztgoekmaSxtA6lPTc86Mk7Va32O3AgRGxO/BxoH9yS2lmZm0VPIB9gCURcVtEPA6cAxxeXCAiLo+I+/LHK4D6pxaamdkEa7fgMRtYUfi8Mqc18zbgJ81mSuqTtFjS4lWrVlVURDMza7fgoQZpDTuiSDqIFDw+0GxjEdEfEXMjYu6sWbMqKqKZmVX2JMGKrAR2KnzeEbijfiFJuwNfBw6NiD9PUtnMzCxrtyuPK4FdJO0saRPgLcB642NJmgOcCxwTEb9vQRnNzKa9trryiIjVkk4CfgrMBBZFxI2STsjzzwT+GXg68JX0eBFWN+o6b2ZmE8djW5mZWVMe28rMzCrj4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBhZmalOXiYmVlpDh5mZlaag4eZmZXm4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBhZmalOXiYmVlpDh5mZlaag4eZmZXm4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBhZmaltV3wkHSIpFslLZF0aoP5kvSFPP86SXtNVFkGBgbo6elhxowZ9PT0MDAwMFFZmZl1lI1aXYAiSTOBLwOvAlYCV0q6ICJuKix2KLBLnvYFvppfKzUwMEBfXx9DQ0MADA4O0tfXB0Bvb2/V2ZmZdZR2u/LYB1gSEbdFxOPAOcDhdcscDnwrkiuAbSTtUHVB5s+f/2TgqBkaGmL+/PlVZ2Vm1nHaLXjMBlYUPq/MaWWXAUBSn6TFkhavWrWqVEGWL19eKt3MbDppt+ChBmkxhmVSYkR/RMyNiLmzZs0qVZA5c+aUSjczm07aLXisBHYqfN4RuGMMy4zbwoUL6erqWi+tq6uLhQsXVp2VmVnHabfgcSWwi6SdJW0CvAW4oG6ZC4Bj811X+wEPRMSdVRekt7eX/v5+uru7kUR3dzf9/f1uLDczo83utoqI1ZJOAn4KzAQWRcSNkk7I888ELgQOA5YAQ8DxE1We3t5eBwszswbaKngARMSFpABRTDuz8D6Ad052uczMbJ12q7YyM7MO4OBhZmalOXiYmVlpDh5mZlaaUvvz1CdpFTDY6nKMw7bAn1pdiAng/eo8U3XfvF+NdUfEBr2sp03w6HSSFkfE3FaXo2rer84zVffN+1WOq63MzKw0Bw8zMyvNwaNz9Le6ABPE+9V5puq+eb9KcJuHmZmV5isPMzMrzcHDzMxKc/Boc5KWSbpe0jWSFre6POMhaZGkeyTdUEh7mqSLJf0hvz61lWUciyb7tUDSH/Nxu0bSYa0s41hI2knSJZJulnSjpHfn9I4+ZsPs11Q4ZptK+q2ka/O+nZbTKz9mbvNoc5KWAXMjouM7L0k6AHiY9Az6F+a0zwD3RsSnJZ0KPDUiPtDKcpbVZL8WAA9HxGdbWbbxkLQDsENEXC1pS+Aq4AjgODr4mA2zX39L5x8zAZtHxMOSNgb+F3g38DdUfMx85WGTJiIuA+6tSz4c+GZ+/03SP3FHabJfHS8i7oyIq/P7h4Cbgdl0+DEbZr86XiQP548b5ymYgGPm4NH+AviZpKsk9bW6MBNg+9qTIPPrdi0uT5VOknRdrtbqqKqdepJ6gBcDv2EKHbO6/YIpcMwkzZR0DXAPcHFETMgxc/Bofy+LiL2AQ4F35ioSa39fBZ4N7AncCZzR0tKMg6QtgB8Ap0TEg60uT1Ua7NeUOGYRsSYi9gR2BPaR9MKJyMfBo81FxB359R7gPGCf1paocnfnOuhaXfQ9LS5PJSLi7vxPvBb4Gh163HK9+Q+AgYg4Nyd3/DFrtF9T5ZjVRMT9wKXAIUzAMXPwaGOSNs8NekjaHHg1cMPwa3WcC4C35vdvBX7YwrJUpvaPmr2BDjxuufH1P4CbI+JfC7M6+pg1268pcsxmSdomv98MeCVwCxNwzHy3VRuT9CzS1Qak583/Z0QsbGGRxkXSfwHzSENE3w18FDgf+C4wB1gOHBkRHdX43GS/5pGqPwJYBry9VufcKSS9HPgVcD2wNid/iNQ+0LHHbJj9OorOP2a7kxrEZ5IuDr4bER+T9HQqPmYOHmZmVpqrrczMrDQHDzMzK83Bw8zMSnPwMDOz0hw8zMysNAcPazuSdpYUkhYV0vbPaR9tZdlGI4+EHHnoi8nKMyT51kmbNA4e1o72y6+/LqS9tEGambXIRq0ugFkDteBxeSFtf1Lnrd9suLgBu7a6ADa9OHhYO9oPeAC4qZC2P3BTRDzQmiK1t4i4pdVlsOnF1VbWViRtQhoi4orIwx9I6gaeCVwxxm0uyG0CC3J7ytmS7pb0aH7a2vskNfwhpeQYSZdKui+vs1TSlyXtVLIc3ZI+mJ9it0LSY5LuzZ+PbrLOvFz2SyV1SfqEpFsk/SUPu11bboM2D0ln1dKHmS5tkOdOkj4v6dacz4OS/k/ScXlcqPrlL83bmifpAEn/LelPktZKOiIvs6mkUyVdLenhvO93Svp13qdNy3yX1nq+8rCWa9LQ+5oG6W+T9Lb8fjAiekpmtTOwGHiUNNroVsBBwGeBl0t6Yx5RtVYuAWcDRwNP5HXuJY22eiLwFkmHRMSVo8z/GODjwFLSYHX/Rxo2+xXAPEn7RsS7m6y7ac5/V+Ay4FpgkxHy+99h5r2cNPz4mmKipINI46ltDSwBLgK2IF0NfgM4GDi2yTaPBE4gXTFeTBrr6wlJM4D/zus+APwyv24PPA+YD3wJuGuE/bF2EhGePLV0Ar5emG4itW18v5C2PKd9o5D26RLbX5DXr21308K8XYCVed6JdeudmNPvAl5QSJ8JfIF1A+g9pW69ZXleT136S4rbqStDbR/3rZs3r1D235Ee6tNoHyP9O4/q+5gHPAYMFfMDdiAFx9WkkVdVmLdTzj+A4+q2d2mhjH0N8jsgz7uK9IjU4jwBLwO6Wv136Knc1PICePJUnEi/xh8GNimkLQduHcc2a8FjCNiuwfzj8/w/1KUvzen/0GCdTQon/N66eQ2Dxwhl/Ie8zr/UpReDx0uHWX9UwYN05XIf6YrjDXXzTs/bOb3JunNrQaAuvRY8ftZkvSPz/M+1+u/LU3WTq62sbUjailQldHFEPJ7TdiH96j2zgix+FumhWvUGSFczz5E0OyL+KGlH4FmkIbu/Xb9CRDwuaQA4lXSCHxhNAXLd/mtIVyGzgKfkWbVnSTy3yap3R8TlTeaNiqTtgQuBbYD3RMR5dYscll+/12QTV5EC+56SNo2IR+vmn9tgHYCrScHqbZJ+D/wgIu4uW35rLw4e1k7mkf4mf15IOzi//qKC7d/eKDEHgjuB2aQ2iD/m9wB3NjhJ1izNr7ObzF+PpP1Jz1TYcZjFtmqSPjiaPIbJuwv4EdADfCEiPtdgsWfl1ysbtIvXezrpeypqWMaIWCrpPaS2pS8DX5Z0G+lW7B8C50XEmkbrWvty8LCWkTQPuKTBrDMk1T8/+ruFE9o3I+K4CSpWrZFedZ8bGfEM++SC6eR9HqmR+D9Iz8teAjwUEWslvRr46TDb/Mto82qQ9wzgP0lXOz8E3tNk0Zn59TukmwqG81iDtKZljIgvSvoecASpsf7lwP/L0zWSDowp9Hz06cDBw1rpLtJTz2reRGqsPT9/FukOpbuAnxWWG+4uouH0NErMtwfXqo3uyK8r8+szJT0lIhqdLHfOr/W/wBs5gBQ4roqIv28w/zmj2MZY/RtwOPBb4Ogo3FFWZ0Uux8cj4saqCxERd5GqH88EkLQHqUpwT1L134eqztMmjoOHtUykjm3HASg9d/kYUn14Le1FpNtCz4mI91aQ5aslzYqIVXXpR5H6PC2NiJW5bCtz1cqzSL+O/6O4gqSNSbfwQmowHsnT8uuKJvMb9vMYL0mnAO8iVdm9PiKGhln8J8DJpAbuyoNHvYi4VtLnSe1Ne0x0flYtdxK0djGP9PdYrMY6KL9W0d4B0AV8SVKtkRpJzyb1vQD4fN3y/5pfPy7p+YV1ZgKfAbpJ9fzfH0XetR7gB9dta4akfybdrlqp3EHvDNLdVYc1uVmg6F+AB4EPSXpno46TkvaTdGTJchws6bD67eXvsdZIP642HZt8vvKwdlFrGK8PHmtIneKq8G3gtcBSSf9H6vx2MKkD3o9IjblFXyGd1I8CrpV0CelEvA/piuQ+4MgmVVrriYirJf0IeD2pjv8SUke5lwBzSMHon8a9h+s7gxSQB4FTmzSC3xIRn85lXJEDzvdJnfbmS7oR+DOph/+z8+t3aH5HViO7k6rOHpB0NXAnKZDvS6ouvIt0m7B1EAcPaxcHAXflqqxa7+4DgMUVNqTeRjpZf5IUNLbOaYtIfRDWawuIiJDUS6rO+QdSL+vNSO0iXwU+FRHNqqEaeROpsfoY0pXWw6RRgo/O2606eNQawPfMUyO/BD5d+xARl0h6Aan66rWkfd6YdIL/A/BFygUOSIF5G9LxfA5phOSHSf1kzgS+2qAq0dqcIvwIAJvaJC0APgqcFhELWlsas6nBbR5mZlaag4eZmZXm4GFmZqW5zcPMzErzlYeZmZXm4GFmZqU5eJiZWWkOHmZmVpqDh5mZlfb/ASLqkALiG3TNAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the result\n",
    "plt.figure()\n",
    "plt.plot(Narray, I, 'ko')\n",
    "plt.xlabel('# polarizers', fontsize=22)\n",
    "plt.ylabel('$I_{transmited}$', fontsize=22)\n",
    "plt.title('Perfect polarizers', fontsize=24);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to show how the transmitted intensity gets closer to 1 using logarithmic scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "Npoins = 100\n",
    "Narray = np.logspace(1,5, Npoins)\n",
    "J = Mueller('System')\n",
    "I = np.zeros(Npoins)\n",
    "\n",
    "# For loop to calculate the intensities\n",
    "for ind, N in enumerate(Narray):\n",
    "    # Create the system of polarizers\n",
    "    angles = np.linspace(90*degrees, 0*degrees, int(N)) # Start by 90º and end in 0º because the rightest  element is the first to be crossed by the light wave\n",
    "    J.diattenuator_perfect(azimuth=angles)\n",
    "    # Multiply them\n",
    "    J.prod()\n",
    "    # Multiply the system matrix by the incident light wave to obtain the output\n",
    "    S_final = J * S0\n",
    "    # Store the intensity\n",
    "    I[ind] = S_final.parameters.intensity(out_number=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEtCAYAAAAr9UYgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAt/UlEQVR4nO3deZhcVZ3/8fcnIRGaVZKACnQ3ICq4gBhZXDDghqiD6yhGFFwygOAybmg7P1En7oyiMmKPgyi04oriiAg6Io4rYYnsGiCbbCEBAgYNSb6/P86p5KZS1V3VfbuW7s/ree5TVecu59Tt7vvte7ariMDMzKwsU9pdADMzm1gcWMzMrFQOLGZmVioHFjMzK5UDi5mZlcqBxczMSuXAYh1NycmSrpG0RlLkpb/dZes27Tx3ki7LeR/X6ryt9RxYDEnnFC46xWV1vqB/RtLubSreB4EvAvsDAu7Ky/pWF0TSTpJOk3Raq/M26yYOLFb0MJsu3HcD25Eu6O8BrpX0rDaU6R359V+Bnoh4VF6WtaEsOwEfzos1ZylwM3B/uwti42+rdhfAOspvI2JO5YOkHuCVwBdIF9XvStorIh5qRWEk7QLMyh//KzxNRNeKiDe0uwzWOr5jsboiYk1EnAu8PSc9CnhZC4uwTaEsD7YwXzMbAwcWa8R3gA35/dOKKyRNkXSspEslrZC0VtLtkr4t6eBaB8vtFJHbdqbkxvk/Srovp79TUgCLC/sU235OG2sZCvtuK+k9kn4raZWkv0u6VdKFkuZKmpa3uwy4rU55tijTMPn1V/bJn58p6X9yudfkNq2TJQ37tynpFZIuzvv9Q9JySUOSDmykHFXHmirpcElnSLpS0l2Fc3iBpCOG2Xdjo3xug/qUpJvyd7mv1na1zkUDS3+dc/lFSTfn/B7I5X+/pG3rlHfj8STtK+nrkpZJeljSDwvb7ZLbFq+T9Lf8e7Es/558VFJfs+d5UokIL5N8Ac4BArhsmG3uytsMFtK2By7N6UEKPvcXPq8HTq5xrNPy+q8DP8zv1wH35vfPAO4EVhSOdWdhec9Yy5D33Y8ULCrbPlwoQ2Xpz9v+YJjybFamEc51f+EYr8x5Rs734cK6C4Ctauw/JZ+3ynbF81b5vifWyXuz71RIf1LVd/478GBV2gfrHPOyvP69wC2F/VcD99XY7rhC2h41zmNxuW+YMr8CeKiwfg3wj8LnPwG7DnMOjgX+lt+vzsf6Yd6mD7i96hyvyr9blbQT2v1328lL2wvgpf0LIwQWUpVU5Y/q04X0C3LaQuAoYJucvhPwgfyHvh54ZtXxTsv7PZAvQieSGuYBdgF2yO/7K3/Iw5R9tGXYmdSgHMCtwNHA9Lxue+BZwNnA7oV9RixPA+d64zHyhfOnwJ553bakC/R66lzMgVPZFEA/BGyf03cj3VlWgsthNfatd5F+XN73JcCugAo/iw/lC+sG4OAax7ys8LNcChwJTMnrHltju+MaPE/TgMvzPtcAWxfWPR1Ym8v1SaCX1GNwKnAw8Pu838+GOQcP5DI9KacL2Du/Pztv8xfg2YXv8whSEP4Y8LJ2/9128tL2Anhp/8LIgeXkwh/kK3Pa8/Ln24Cd6+z3vrzN/1Sln1Y43rxhyjXshXyMZfh0Tl8B7NbgeRq2PM0eA7gOeESNbSrn535g20L6tmy6G/tEjf2mAr/O6y+vsb5mYGmgzP+W9/tajXWVgLG2cpGuc4zKdsc1mOdZefu7gb6qdf+X172rzr6PBP6at5ld5xzcQv4npMb+N+RtXjPan/NkX9zGYjUp6Zf0HtJFGGAJ8OP8/o359ZyIWFXnMN/Mr4dLmlpj/UrSf4ejNZYyHJtfPxsRfx1DGcbi9Ij4R430/yDdye0APL+Q/oKctpZNP5ONImI96b9pgGdLelRJ5az8zJ85zDY/jYjryshM0onAv5CqBl8VEUsK6/bO5XiIFHy2EBH3ku4EYfPzV/SlqN+7cXV+fXSTRbfM3Y2t6DmVRuUa7iDd/q/Nn5+RX9+VLwTD6QFmkP77LFoQEetGV9TRlyE3BFcuuheNIf+xuqxWYkSslnQ1cChwIKkdivweYGG+eNZyOamKaKu8fUPfT9I2wAmkKsH9SP/1V18fHjPMIX7XSD4NlOMw4Iz88ZSIuLxqk8rPfDpwm6R6h9ouv+5RZ/1w5b2IVKX2KUn7AN8Dfj9MILIqDixW9DCpkRJSVcDfSO0PlwJfrbqYVf6b2zEvI+mpkbZilOUcaxl2LaQtHWMZxmK4O6XKulmFtFlV67YQEX+XtJL0HWfV265I0qNJQe5xheS/kToFbCBVsc0kVcXVM9afJZJ6SRfxacCXI+IrNTar/MynsvnPsZ5av3cwfHk/Rer9+E/ASXlZJ+kKUpvef0XEfQ3kPWk5sFjRZgMkR1CpRj06Ii4cZX5jnZZltGWo+29uBxmujI8oOa/Pk4LKraTOA78s/hORq58WjXCMMf0slQbj/ogUDH/FphkXqlV+5ldHRNNdqwvqljdXTx4t6RDg5cBhpEBzaF7eK+n5EbFwDPlPaG5jsdG6K7/u14VluLPwvq+ksozGcFVLlf/Mi/9ZV97XLbOkrUlVftX71tt+Oqn6C2BuRPygRjVbI3cGY/U14ADS2KVXRcTDdbar/Mz3kTSu/xhHxO8j4v0RcSipavAY0h3uLOCr45l3t3NgsdGq1FG/stvKEBGL2RRcjmpi18ogUTRM5X4TnlMrUdL2bGpPuaqwqvJ+H0m71TnmYWyqibiqzjZFM9l0B3R1nW2e18BxRk3SAPDPpOq3oyPinmE2r/zMtyN1ZmiJiPhbRJwPzMtJT6s3CNMcWGz0zsmvsyUNOw+UpEd2YBnOza/vHuYiXW114f1ODe4znHfnO4Zq7wS2zvldUki/JKdNI1VZbSb3evu3/PHXEXFn9TY1rCa1pwE8ucYxHw2c0sBxRkXSS4GP5jK8ISL+NNz2EXETaZwKpMb1uhd3SdtIarrasM7PpKLSgC9SBwKrwYHFRiUiLiaNRgc4W9JH8kUISBdySUdL+hGp+2ynleFTpEbwmcCvJf1T5YIiaTtJcySdr8LjAnKD7e354/ElfIVe4ILKdCWSeiT9K5tmT/5URKwp5P834OP549slDUjaLu+7G/At0sDOyuDJEUWag61yoT5b0gH5eFMkPZfU3jEubVKSHgecR7oOfSQifjDCLhWnkAa+Pon0s3tepVosl/uJkj5EGqsymi7D10n6uKSnF34nJOkg0iMcAK4YpmeetXsgjZf2LzQwpUud/bZl08j34mjy+6vSvla132k5/ZwRjt9fOUbZZcj7PhlYVthmLXWmdCns85HCugdJbQKLgXc2eM76C/sPN6XLD6k9pctUtpzSpTjdyHrgpDp51/tOB5OmRCl+r8rnlaQ2mJo/Bxoc+FhrO+C4Qp53Mfz0LntUHe9FbD7lyz+Ae/LPsPjz62vkHFRtUzzuunwOisddATyl3X+3nby4V5iNWqT/oF8u6cXAm0gXqFmki9wi4I+kO4pxGysyljJExLWSnkj6D/hlpJ5RW5N6R/0J+DawvGq3j5LaAuYCj2VTQ/pOoyj79yUdTpqm5RDSRex6UsPwf0bEhhr7rAfeKOlCUn3/00iDJu8g3V2cHhFXNlmOP0g6lBTwn0MK1ncAFwPzScFsvO0ywvrNyhARP813PKeQ2skeS/oZ3Ed67svFwHejMLiyCUcDLyS1V/WSOi+sBW4k/R59LiKqx2RZQWVOIDNrgVztdRtARHRDt2ezprmNxczMSuXAYmZmpXJgMTOzUjmwmJlZqdx4b2ZmpZr03Y1nzpwZ/f397S6GmVlXufLKK++JiJozaE/6wNLf38+CBQvaXQwzs64iqe4YIbexmJlZqRxYzMysVA4sZmZWKgcWMzMrVccFFklnS7pb0nV11kvSFyQtkvQnSQcW1h0p6ea87tTWldrMJruhoSH6+/uZMmUKM2fOZObMmS1/39/fz0knndRUOfr7+xkaGir3ZLR7euUaU1YfRnp63nV11h8F/JT0jIhDgD/k9Kmk5y/sRXoAz0Jgv5Hye9rTnhZmk8l5550XfX19ISlmzJgRM2bMaOv7vr6+OPHEEzuqTM2WFQhJxen6u2rp6emJ8847r6nfI2BB1LuO11vRzoX0zIp6geUrwDGFzzeTHuZzKPCzQvoHgA+MlJcDi3WqsgJA8WLY7RdAL+O39PX1NfX7yTCBpRvHsexGejhTxfKcViv94FoHkDSP/Ozq3t7e8SmlWZWhoSEGBgZYunQpO++8MwCrVq2q+X7lypVIqvyTxMqVKzcep9n3S5Ys4ctf/vLGz5VjmhUtXbq0tGN1Y2Cp9QyLGCZ9y8SIQWAQYPbs2f4rs1EZ70DhAGCtVOY/2d0YWJYDexQ+7056Dvn0OulmozZc8HjggQdYu3Yt4EBh3a2np4f58+eXdryO6xXWgAuBN+TeYYcA90fEHcAVwD6S9pQ0HXht3tZsRLV69Eji2GOPZcmSJUQEK1euZOXKlRvfV4KKWYWUKk5mzJjBjBkzkNTS9319fZx44on09fU1tc/g4CBz584t7Tx03B2LpG8Bc4CZkpYDHwamAUTEWaRnTh9Fep75GuD4vG6dpJOBn5F6iJ0dEde3/AtY16jcjSxZsqRuNdVEvcuofN8ZM2YA9avwWvG+t7eXo446iosuuqihasVOLWtvby/z588v9QLdrTousETEMSOsD+BtddZdRAo8ZhvVqs6qbvPo1ABSRgCovhj6AmjjreMCi1kZOvVupJlA4QBg3cqBxSaMesFkvALItGnT2GGHHZq6c3CgsMnAgcW6WiuCSa27DAcJs/ocWKzrjFcwcQAxK4cDi3WF8Q4mfX19DiBmJXFgsY5VZjDx3YhZ6ziwWMeo7hZcHNk+lmDiuxGz1nJgsY4wNDTEvHnzWLNmDbB5t+BmOJiYtV83TuliE0hlKpXXv/71G4NKsyrTaPT19XHuuecSESxevNhBxaxNfMdiLVev7aQZvjMx61wOLNZS1VVezQQVBxOz7uCqMGuJ0VR5TZs2bbMZWF3NZdYdfMdi4676LqURvisx614OLDZuim0pjerp6Sn92RBm1lquCrNxUblLaSSoFHt1OaiYdT8HFitVs20pbjsxm3hcFWalaaYtxVVeZhOX71isNAMDAw3fpTiomE1cDiw2ZpXqr5HaU3p6ejjvvPNc5WU2wTmw2Jg02kjvuxSzycNtLDYqjXYldluK2eTjOxZrmu9SzGw4vmOxpjXSSN/X18fixYtbUyAz6yi+Y7GGNdNIP3/+/BaVysw6jQOLNcTVX2bWKFeFWUNGqv5yI72ZVfiOxYbVSPWX71LMrMh3LFZXI1O0uJHezKr5jsXqaqT6y430ZlbNgcW24OovMxsLV4XZZlz9ZWZj5TsW24yrv8xsrDousEg6UtLNkhZJOrXG+kdKukDSnyT9UdKTCusWS7pW0jWSFrS25BPD0qVL665z9ZeZNaKjqsIkTQXOBJ4PLAeukHRhRNxQ2OyDwDUR8XJJT8jbP7ew/vCIuKdlhZ4gKpNKRkTN9a7+MrNGddody0HAooi4NSLWAucDR1dtsx/wC4CIuAnol7Rra4s5sYw0qt7VX2bWjE4LLLsBywqfl+e0ooXAKwAkHQT0AbvndQFcIulKSfPqZSJpnqQFkhasWLGitMJ3q+HaVVz9ZWbN6qiqMEA10qrrZj4JnCHpGuBa4GpgXV73zIi4XdIuwKWSboqIy7c4YMQgMAgwe/bs2nU/k0i9dhVJrv4ys6Z12h3LcmCPwufdgduLG0TE6og4PiIOAN4AzAJuy+tuz693AxeQqtasjsp4lXrtKr29vS0ukZlNBJ0WWK4A9pG0p6TpwGuBC4sbSNoprwN4C3B5RKyWtK2k7fM22wIvAK5rYdm7ittVzGy8dFRVWESsk3Qy8DNgKnB2RFwv6YS8/ixgX+AbktYDNwBvzrvvClwgCdL3+mZEXNzq79AtRmpXmT9/vttVzGxUVK8aZLKYPXt2LFgw+Ya8TJkypWYVmCQ2bNjQhhKZWTeRdGVEzK61rtOqwqxF6rWfuF3FzMbKgWWSKU4wmasNN3K7ipmVwYFlEqlusI+IjcHF41XMrCwd1Xhv46tWg31EeLoWMyuV71gmkXoDIYebeNLMrFkOLJOIG+zNrBUcWCa4SmP9lClTePDBB5k+ffpm691gb2Zlc2CZwIqN9RHBypUriQhmzJiBJDfYm9m4cOP9BFarsf7hhx9mu+224557/MgaMxsfYwoskn5Bml34auAa4MaI8LDtDuHGejNrh7HesdwFvBh4J2nK+39IupYUZCoB508R8dAY87FR6O3trTnJpBvrzWw8jamNJSJeFxH7AjsCzyNNe38g8BLg88BvgdWSbpA0NMayWoM8ut7M2qmUxvuI+BvwMuAWYJeI2A3YgTTt/Q2kpzzOKiMvG55H15tZu5XZeH8c8MqIuBcgP7P+u5J+SJoGf4snOVr5PLrezNqtzO7Ga4BHVSdGxMPA6cDxJeZldbjB3szarczAMgR8RNJuNdZtDTy6xLysDo+uN7N2KzOwfAi4FbhZ0uckHSlpP0kvBz4NXFtiXlbH/Pnz6enp2SzNDfZm1kqlBZaI+DvwQuDDwD8DF5GCyffzJvPKysu2VOkJduyxx7LNNtt4dL2ZtU2pI+/z4MjTgdMl7Qs8BlgJXBcR68rMyzap9ASrNNqvXLmSnp4ezj33XAcUM2u5Up95L+kw4E2k7sX3A2+PiI5uNZ4Iz7yvjFmp5p5gZjZeWvLMe0lzgcuAPYDrgJcCM/O60yQ9p6y8bHPuCWZmnaTsxvvPRsRzgXeRpnipeJg07YuNA/cEM7NOUmZg6Qcuzu+r69euI031YuPAPcHMrJOUGVhuAZ5cZ929eEqXcTN37lwGBwfp6+tzTzAza7sye4V9Afi4pF8DC6vW7QOsKDEvqzJ37lwHEjPrCGWOYxkEfkWa0fjrpOqwR+WeYgPApWXlZUnxscP9/f0MDXkCaTNrv7KfIPkq4G3A+0iN9z/Or1eQgouVpHrsypIlS5g3L41B9Z2LmbVTaeNYJJ0EDEXE/flzZYDkHaQnS5Y3YKZE3TqOxWNXzKydWjKOBfgi8PjKh4i4MSJ+ERE3ANMkbVNiXpOex66YWacaU2CRtLOk6ZWPw2z6VGD1WPKyzXnsipl1qrHesZwCPCjpelJj/RskvUTS7lXbbQOsH2NeVuCxK2bWqcbaeH82cBuwP7Av8DrgJCAk3QtcA/wZOIT0iGIbo6GhIQYGBli6dCk777wz22yzDatWraK3t5f58+e74d7M2m5MdywRsSwivhER7waWA0eSBkI+H5gP/BV4GmkMS0PT5ufnuNwsaZGkU2usf6SkCyT9SdIfJT2p0X27XfF59hHBypUreeihhzj33HNZvHixg4qZdYQx9QqTtFNE3FdaYaSppDuc55MC1RXAMbkDQGWbzwAPRsRHJD0BODMintvIvrV0U68w9wQzs04xnr3CVkjaP2fyAUkvqvNo4kYdBCyKiFsjYi1wPnB01Tb7Ab8AiIibgH5Juza4b1dzTzAz6wZjDSyvAe7O798B/ARYKukeSb+Q9B+S3ijpAEnTGjjebsCywuflOa1oIfAKAEkHkZ79snuD+3Y19wQzs24w1jaWH0TEHfn9o4BdSY8n/gTpwn4EMAhcBTzYwCFrdVmurqv7JPBISdeQeqVdDaxrcN+UiTRP0gJJC1as6J4pzNwTzMy6QdmPJl4B/DwvAEjailR99ZQGDrGc9KCwit2B26vyWA0cn48tUq+024CekfYtHGOQFPCYPXt2R84IUEulcb7SK8w9wcysE5U9V9gW8rPu/5SXkVwB7CNpT1KPsteSujBvJGknYE1uR3kLcHlErJY04r7dqtjF2MHEzDpdaYElN+J/jjSe5T7Sw70WksayLIyIZXV3ziJinaSTgZ8BU4GzI+J6SSfk9Wfl439D0nrS2Jg3D7dvWd+vXTzZpJl1mzInoVwI7EAaNLk9qeprf1K7SwD3RcSMUjIrUad3N3YXYzPrRMN1Ny6zKmxv4DUR8ZOqzHcBDqCxNhar4i7GZtZtypzd+Apgu+rEiLg7Ii6JiM+WmNek4S7GZtZtygwsHwBOlbRjicec9NzF2My6TZmB5Y/A74DrJL1f0tMlbV3i8SeluXPnMjg4SF9fH5Lo6+tjcHDQDfdm1rHKbLz/KvAm4AHgEcB00lT5fyb1DLu6E6vDOr3x3sysE7XqCZKvAeZHxI7AtqQeYccDPyX1DHtfiXlNeENDQ/T39zNlyhT6+/sZGhpqd5HMzBpSZq+we4FfAUTEeuDavJxXYh6TgseumFk3K/OO5SvAy0o83qQ1MDCwMahUrFmzhoGBgTaVyMyscWUGlmnASyR9RNIW3Y6tcR67YmbdrMzA8i9AL/BvwN2SfiXpDEnHS3pqg9PmGx67YmbdrbTAEhGPJjXSv4AUXBYDzwHOAq6ksWnzDY9dMbPu1mnT5hueHt/MuluZ41jmAesi4uxSDtgiHsdiZta8Vo1jeTswq5DpHEmrJD0g6RMl5mNmZh2szMDST3oEccUngFuA04ATJb2pxLzMzKxDlRlYHgI2AEjqAw4GTomI04GPAG8rMa8JyaPtzWwiKDOwLASen9/PBe6KiN/nz1cB+5SY14RTGW2/ZMkSImLjaHsHFzPrNmUGlvnAOyT9nNTd+NzCul2BtSXmNeF4tL2ZTRSldTeOiF9KOgI4BlhAqv6qmEOa5djq8Gh7M5soyh7H8jvSM1mqTQW+U2ZeE01vb2/NZ9t7tL2ZdZsyq8Lqioh/iYjPtyKvbuXR9mY2UZR2xyLpkcB7gccBK4DrSA36CyPigbLymag82t7MJooyR97/BDgC+F9gZ+CJwHZAkOYNuzoiXlVKZiXyyHszs+YNN/K+zDaW5wBvK07pIumxwAHAU/FcYWZmk0KZgeUO4K/FhIhYBCwCvldiPmZm1sHKbLz/MvDKEo9nZmZdqMzAsjVwhKR/l7R9icc1M7MuUmZgeSewF/BB0hMkL5f0BUlvknSgpOkl5jVheH4wM5toyhx5v4ukR5Ma6SvLHNIji6cBDwOPKCu/iaAyP1hlKpfK/GCAuxmbWdcqrbtx3QzSs+73A54cEeeNa2aj0M7uxv39/TVH2/f19bF48eLWF8jMrEGt6m5cU0Q8TB4oOd55dRvPD2ZmE1GpU7pIeoKkN0p6paR9JanM40809eYB8/xgZtbNSgsskl4DXEvqdvwd0pQuD0r6g6SvSDqhweMcKelmSYsknVpj/Y6SfixpoaTrJR1fWLdY0rWSrpHU8cPpPT+YmU1EZd6xfIgUVHYARJo37OvAk4FXAaePdABJU4EzgReR2mWOkbRf1WZvA26IiP1JnQNOr+pxdnhEHFCv7q+TzJ07l8HBQfr6+pBEX18fg4ODbrg3s65WZmDZG/gh+fHEwK8i4iRgNrAKOKSBYxwELIqIWyNiLXA+cHTVNgFsn6vZtsvHXjf24rfH3LlzWbx4MRs2bGDx4sUOKmbW9coMLKtJvcw25PczACLiBuDzwCcbOMZuwLLC5+U5rehLwL7A7aSqt3fkPCEFnUskXSlp3ii/h5mZjUGZgeVa0l0LwK3AgYV1NwLPbuAYtRr7q/tDvxC4BngMaYLLL0naIa97ZkQcSKpKe5ukw2pmIs2TtEDSghUrVjRQLDMza1TZc4U9Or//GvBOSUdI2g04AbivgWMsB/YofN6ddGdSdDzwg0gWAbcBTwCIiNvz693ABaSqtS1ExGBEzI6I2bNmzWrku5mZWYNKCSy5veN/IqLynPuvkO4qfg4sBV4O/FsDh7oC2EfSnrlB/rXAhVXbLAWem/PdFXg8cKukbStzlEnaFngBqWeamZm1UFkDJKcAD0g6OiIuzg3vR0o6lHQHclW+uxhWRKyTdDLwM2AqcHZEXF/pqhwRZwEfA86RdC2p6uz9EXGPpL2AC/LQma2Ab0bExSV9PzMza1ApgSUi1ku6u0b674DfNXmsi4CLqtLOKry/nXQ3Ur3frcD+zeRlZmblK7ON5WzgDSUez8zMulCpU7oAz5E0KGlGycedUDxVvplNZGVOQjkP2BV4C3B8nlJlAakR/2rgutz2Mql5qnwzm+hKnTZf0ixSO8f+pOex7E8azDgNeDgiOu55LK2eNt9T5ZvZRNCSafMlnQQMRcTPSd2MK+lbkeb9ekpZeXUzT5VvZhNdmW0sXySNKdlMRKwDbgK+X2JeXctT5ZvZRDemwCJp58LMwsM9e+WppPnDJj1PlW9mE91Y71hOIT1z5XrSnF5vkPQSSbtXbbcNsH6MeU0InirfzCa6MTXeS9oDOJzUSP8u0nxgO5GCzL2kHmF/Jk+ZnyeI7CjtfOa9mVm3GrfG+4hYBnwjZ/Jq0gO9bmFTz7ADgKcBK4CBseRlZmbdobReYRFRbH3+37yYmdkkM2IbS35uyaCkEyQ9XVLHjUXpBh5tb2aTRSN3LAeSenW9OX9eL+lG4Crgyvx6TUSsGZ8idj+PtjezyWTExntJx5OCy4GkQY7bFlZXdt4A/IVNgeYq4OqI6Pguxq1ovPdoezObaIZrvG+qV1h+oNcT2BRoDiQ10O9Y2Kx4wFsi4nHNFriVWhFYpkyZQq3zLIkNGzaMa95mZuOhtF5hka6ON+ZlYyOBpL3ZPNg8FZgJ7D3KMk8ovb29Ne9YPNrezCaiUqZ0iYhbIuK7EfGBiHhhROwC9JEeSTzpebS9mU0mZT+PZaOIWBYR1c+rn5Q82t7MJpNSp83vRh55b2bWvOHaWMbtjsXMzCYnBxYzMyuVA8s48mh7M5uMynzmvRV4tL2ZTVa+YxknAwMDG4NKxZo1axgY8CTPZjaxObCMEz/b3swmKweWceJn25vZZOXAMk482t7MJisHlnHi0fZmNll55L1H3puZNc0j783MrGUcWMzMrFQOLGZmVqqOCyySjpR0s6RFkk6tsX5HST+WtFDS9fnRyQ3ta2Zm46+jAoukqcCZwIuA/YBjJO1XtdnbgBsiYn9gDnC6pOkN7mtmZuOsowILcBCwKCJujYi1wPnA0VXbBLC9JAHbAauAdQ3uO+488aSZTXadNgnlbsCywuflwMFV23wJuBC4HdgeeE1EbJDUyL7jyhNPmpl13h2LaqRVD7R5IXAN8BjgAOBLknZocN+UiTRP0gJJC1asWDH60lbxxJNmZp0XWJYDexQ+7066Myk6HvhBJIuA24AnNLgvABExGBGzI2L2rFmzSiu8J540M+u8wHIFsI+kPSVNB15LqvYqWgo8F0DSrsDjgVsb3HdceeJJM7MOCywRsQ44GfgZcCPwnYi4XtIJkk7Im30MeIaka4FfAO+PiHvq7dvK8nviSTMzzxVW+lxhQ0NDDAwMsHTpUnp7e5k/f74b7s1swhlurjAHFk9CaWbWNE9CaWZmLePAYmZmpXJgMTOzUjmwmJlZqRxYzMysVA4sZmZWKgcWMzMrlQOLmZmVyoHFzMxK5cBiZmalcmAxM7NSObCYmVmpHFjMzKxUDixmZlYqBxYzMyuVA4uZmZXKgcXMzErlwGJmZqVyYDEzs1I5sJiZWakcWMzMrFQOLGZmVioHFjMzK5UDi5mZlcqBZRSGhobo7+9nypQpzJw5k5kzZzJlyhT6+/sZGhpqd/HMzNpqq3YXoNsMDQ0xb9481qxZA8DKlSs3rluyZAnz5s0DYO7cuW0pn5lZu/mOpUkDAwMbg0ota9asYWBgoIUlMjPrLA4sTVq6dGkp25iZTVQOLE3q7e0tZRszs4nKgaVJ8+fPp6enp+76np4e5s+f38ISmZl1FgeWJs2dO5fBwUH6+vqQxIwZM5gxYwaS6OvrY3Bw0A33ZjapKSLaXYbNSDoSOAOYCnw1Ij5Ztf69QOXKvRWwLzArIlZJWgw8AKwH1kXE7JHymz17dixYsKDEb2BmNvFJurLeNbajuhtLmgqcCTwfWA5cIenCiLihsk1EfAb4TN7+pcC7ImJV4TCHR8Q9LSy2mZkVdFpV2EHAooi4NSLWAucDRw+z/THAt1pSMjMza0inBZbdgGWFz8tz2hYk9QBHAt8vJAdwiaQrJc0bt1KamVldHVUVBqhGWr1GoJcCv6mqBntmRNwuaRfgUkk3RcTlW2SSgs48cNdgM7Oyddody3Jgj8Ln3YHb62z7WqqqwSLi9vx6N3ABqWptCxExGBGzI2L2rFmzxlxoMzPbpKN6hUnaCvgz8Fzgr8AVwOsi4vqq7XYEbgP2iIi/5bRtgSkR8UB+fynw0Yi4eIQ8VwBLCkk7Avc3+HkmMF4dBarzLXO/kbapt75Wus/X+J0vGL9z5vPVvNGcs4l8vvoiovZ/5hHRUQtwFCm43AIM5LQTgBMK2xwHnF+1317AwrxcX9l3FPkPNvoZWDCO52FwvPYbaZt662ul+3yN3/kaz3Pm89WaczZZz1entbEQERcBF1WlnVX1+RzgnKq0W4H9SyjCj5v8PF5Gm08j+420Tb31tdJ9vny+mt2mG8/XaPOalOero6rCuo2kBdHAIExLfL6a53PWHJ+v5ozX+eq0xvtuM9juAnQZn6/m+Zw1x+erOeNyvnzHYmZmpfIdi5mZlcqBxczMSuXAYmZmpXJgKZGkvST9t6Tvtbss3UDSyyT9l6QfSXpBu8vT6STtK+ksSd+TdGK7y9MNJG2b5w58SbvL0ukkzZH06/w7Nmcsx3JgGYGksyXdLem6qvQjJd0saZGkUyGNpYmIN7enpJ2hyfP1w4h4K2nA62vaUNy2a/J83RgRJwD/DEzKLrXNnK/s/cB3WlvKztHk+QrgQWBr0vRaozdeo1QnygIcBhwIXFdIm0qaGWAvYDpptP9+hfXfa3e5u+x8nQ4c2O6yd8P5Av4J+C1pqqO2l7+TzxfwPNKcgscBL2l32bvgfE3J63cFhsaSr+9YRhBpduRVVcnNPjdm0mjmfCn5FPDTiLiq1WXtBM3+fkXEhRHxDDY9RXVSafJ8HQ4cArwOeKukSXe9a+Z8RcSGvP5e4BFjybfjpnTpErWeG3OwpBnAfOCpkj4QEZ9oS+k6T83zBZxC+q9yR0mPjaqpeyaxer9fc4BXkP7oL9pyt0mr5vmKiJMBJB0H3FO4cE529X6/XgG8ENgJ+NJYMnBgGZ2az42JiJWkCTNtc/XO1xeAL7S6MF2g3vm6DListUXpCsM+xynS3IK2Sb3frx8APygjg0l3a1iSZp4bYz5fzfL5ao7PV3PG/Xw5sIzOFcA+kvaUNJ3UQHhhm8vUyXy+muPz1Ryfr+aM+/lyYBmBpG8BvwMeL2m5pDdHxDrgZOBnwI3Ad6LqYWSTlc9Xc3y+muPz1Zx2nS9PQmlmZqXyHYuZmZXKgcXMzErlwGJmZqVyYDEzs1I5sJiZWakcWMzMrFQOLNZV8qCukHR2Ie3QnPbhdpatEZIW57L2tzDPkORxBdYyDizWbQ7Jr78rpD2jRpqZtYknobRuUwksvy2kHUqadPAPrS9OV9i33QWwycWBxbrNIcD9wA2FtEOBGyLi/vYUqbNFxE3tLoNNLq4Ks66RJ8w7APh95LmIJPUBjwF+P8pjnpbbIE7L7TfnSbpL0t8lXS/p3ZJq/gOWH1R2rKTLJN2b97lF0pmS9qi1zzDl6JP0AUm/lLRM0j8krcqfX1dnnzm57JdJ6pH075JukvSQpGsK223RxiLpnEr6MMtlNfLcQ9IZ+bG2D0laLek3ko6TtMV07Llskct6mKSfSLpH0gZJL8vbbC3pVElXSXowf/c7JP0uf6etmzmX1n6+Y7GOVqfR+YU10t8s6c35/ZKI6G8yqz2BBcDfSc882YH0BMLPAs+S9Mrig6LyRfQ80tMJH877rCI9ne8k4LWSjoyIKxrM/1jgY6RHxt4E/IY0nfmzgTmSDo6Id9TZd+uc/77A5aRHzU4fIb//G2bds4C9gfXFREmHAxcAOwKLgIuB7Uh3kV8DjgDeUOeYryY9q+gG4FJgJvCw0lMdf5L3vR/4VX7dFXg8MEB66NSdI3wf6yTtfiazFy/DLcBXC8sNpLaU7xXSlua0rxXSPtnE8U/L+1eOu3Vh3T6kZ1cEcFLVfifl9DuBJxbSp5IeXhbAYuARVfstzuv6q9KfXjxOVRkq3/HgqnVzCmW/Gti1zneM9Kfe0PmYA/wDWFPMD3g0KXCuA95InsA2r9sj5x/AcVXHu6xQxnk18jssr7sS2LZqnYBnAj3t/j300tzS9gJ48dLoQvov/kFgeiFtKXDzGI5ZCSxrgF1qrD8+r/9LVfotOf2tNfaZXggGc6vW1QwsI5TxrXmfz1SlFwPLM4bZv6HAQrrjuZd0p/LyqnWfysf5VJ19Z1cCRFV6JbBcUme/V+f1n2/375eX8hZXhVlXkLQDqZrp0ohYm9P2If23fFYJWVwSEXfXSB8i3QU9VtJuEfFXSbsDewEbgHOrd4iItZKGgFNJF/+hRgqQ2xJeSLp7mUV6tj2kuwWAx9XZ9a6I+G2ddQ2RtCtwEel55++KiAuqNjkqv363ziGuJAX9AyRtHRF/r1pf75G3V5EC2Zsl/Rn4fkTc1Wz5rbM4sFi3mEP6ff15Ie2I/Pq/JRz/tlqJOUjcAexGavP4a34PcEeNC2jFLfl1tzrrNyPpUOA7OY96dqiTvqSRPIbJuwf4MdAPfCEiPl9js73y6xU12uirzSCdp6KaZYyIWyS9i9SWdSZwpqRbSd3JfwRcEBHra+1rncuBxTqSpDnAL2usOl3S6VVp3ylc7L4eEceNU7EqHQZU9bmWEa++GzdMF/YLSA3W/w18mdQ4/kBEbJD0AtLT/uod86FG86qR9xTgm6S7pB8B76qz6dT8+m1SB4fh/KNGWt0yRsQXJX0XeBmp48CzgNfn5RpJz4mI1SPkaR3EgcU61Z3A1wufX0VqOP5h/ixST6o7gUsK2w3X22k4/bUScxfnSlXU7fl1eX59jKRHREStC+me+bX6P/daDiMFlSsj4i011j+2gWOM1ueAo4E/Aq+LQs+3KstyOT4W4/DY34i4k1SleRaApP1J1YwHkKoUP1h2njZ+HFisI0Ua1HccgKSdSEHk+5W7EUlPJnVtPT8i/rWELF8gaVZErKhKP4Y03uuWiFiey7Y8V9fsRfqv+r+LO0iaRuqGDKnxeiQ759dlddbXHMcyVpLeCbydVA340ohYM8zmPwVOITW2j/vz5CNioaQzSO1b+493flYuD5C0bjCH9LtarBo7PL+W0b4C0AN8SVKlwRxJe5PGlgCcUbX9f+TXj0l6QmGfqcCngT5Su8L3Gsi7MjL+iKpjTZH0/0hdbkuVByeeTuoFdlSdjgtFnwFWAx+U9LZag0YlHSLp1U2W4whJR1UfL5/HSoeBMbUhWev5jsW6QaWRvjqwrCcNCCzDucCLgVsk/YY08O8I0uDDH5Malov+k3TBPwZYKOmXpIv0QaQ7mXuBV9epJttMRFwl6cfAS0ltCr8kDRJ8OtBLClTvG/M33NzppGC9BDi1ToP8TRHxyVzGZTkYfY80YHFA0vXAStLMB3vn129Tv+dYLU8hVcfdL+kq4A5SkD+YVAV5J6mrs3URBxbrBocDd+bqscqo98OABSU26t5KupB/nBRQdsxpZ5PGWGzW9hARIWkuqYroraTR59uQ2mG+DHwiIupVbdXyKlLD+bGkO7QHSbM1vy4ft+zAUmmMPyAvtfwK+GTlQ0T8UtITSVViLyZ952mki/9fgC/SXFCBFLR3Iv08H0uaqfpB0jigs4Av16ietA6nCD+mwSYvSacBHwY+EhGntbc0ZhOD21jMzKxUDixmZlYqBxYzMyuV21jMzKxUvmMxM7NSObCYmVmpHFjMzKxUDixmZlYqBxYzMyuVA4uZmZXq/wPC09b0qycWoQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the result\n",
    "plt.figure()\n",
    "plt.semilogx(Narray, I, 'ko')\n",
    "plt.xlabel('# polarizers', fontsize=22)\n",
    "plt.ylabel('$I_{transmited}$', fontsize=22)\n",
    "plt.title('Perfect polarizers', fontsize=24);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imperfect polarizers\n",
    "\n",
    "Imperfect polarizers will have a leak ($T_{min}$ > 0) and, much more important for the Zeno effect, will not perfectly transmit the transmission state ($T_{max}$ < 1). We will use some values of commercial polarizers with anti-reflecting coating: \n",
    "\n",
    "* $T_{max}$ = 0.99.\n",
    "* $T_{min}$ = 1e-4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nmax = 30\n",
    "Narray = np.arange(2,Nmax+1)\n",
    "J = Mueller('System')\n",
    "I = np.zeros(Nmax-1)\n",
    "\n",
    "# For loop to calculate the intensities\n",
    "for ind, N in enumerate(Narray):\n",
    "    # Create the system of polarizers\n",
    "    angles = np.linspace(90*degrees, 0*degrees, N) # Start by 90º and end in 0º because the rightest  element is the first to be crossed by the light wave\n",
    "    J.diattenuator_linear(Tmax=0.99, Tmin=1e-4, azimuth=angles)\n",
    "    # Multiply them\n",
    "    J.prod()\n",
    "    # Multiply the system matrix by the incident light wave to obtain the output\n",
    "    S_final = J * S0\n",
    "    # Store the intensity\n",
    "    I[ind] = S_final.parameters.intensity(out_number=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum obtained intensity is 0.720169568824601 for 17 polarizers.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEpCAYAAABr364UAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsTUlEQVR4nO3de5gcVZ3/8fcnAYQBEZWgGMgMIq6ACuII4gUDqy6yuvECKzI/EHV3FhEV3ItodEXdrJeVXRVRnNWou8wu6wUQFUVdBdwF1gwIyFUTTELkFrlDuCX5/v44p0ml093pStdMd08+r+epp7tPnapzqmumvl116pxSRGBmZlbGjG5XwMzM+o+Dh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBh04KSEyRdKWmVpMjTULfr1m+6+d1JujCXfexUl23lOHj0EUlfz/9YF3a7Lj3og8BpwD6AgNvztGaqKyJpB0mnSDplqss2myoOHjZdvDe/vg8YiIin5+nmLtRlB+AjebJylgM3Avd2uyLW2hbdroBZpyTtBMzKH/81PGxC34qIY7pdB2uPzzxsOtim9iYiHuhmRcw2Fw4e00SxoVHS9pI+LWmJpIck3STpY5K2LuT/Y0kXSPqDpAclXSzp5U3WfUpe99clzZB0kqSr8nJ3SjpP0v4bqd8MSUdL+omklZIelXSLpP+SdECb5Z4g6ZeS7snpJ0oKYGlhmShMp3Rah8Ky20r6G0mXSLpL0sP5ez1P0oikLWv7Afhdk/psUKcW5Q3VlsmfXyrp+7neq/KNASdIavk/LOmNkn6Ul3tE0gpJ45L2a6cedeuaKelgSZ+TdLmk2wvf4TmSDmmxbPHvcwdJn5J0Q96Wexrla/RdtDENNfkuT5N0Yy7v/lz/90vatkl9H1+fpD0lfUPSzZIek3RuId9Okv5J0jX5/+HhnO8Spf+5wbLfc9+ICE99MgFfBwK4sMG8C/O8k4Dr8/sHgEfz+wDOy3mPB9aSGpPvLcx/BHhpg3Wfkud/A/hOfv8YcE9h2dXAm5vU+4nATwp519aVuwY4YSPlnlso5+78/iXAbcDKwrpuK0x/02kd8rJ7kQJCLe9jhTrUpqGc9+wW9VmvThvZ10OFdbwplxm53McK884Btmiw/Iz8vRX3T7HOa4B3Nil7vW0qpD+3bpsfzn9jxbQPNlnnhXn+3wJLCsvfB9zTIN+xhbRdG3yPxemeFnV+I/BQYf4q0t957fPVwNNafAdHAw/m9/fldZ2b8wwCt9R9x3flv61a2nHdPm5M2vGo2xXwVGJntRc87gFuAF6W07cC/qJwwPkwKaD8I7BDzjMIXJLn/7LBuk8prHs1KUBtk+ftDvy48I+5e4Plz8nzrwIOKyy7A/CB/M+8hrrAVSj3/nygeSepMRxgJ2D7/H6o9s/a4rvb1Do8hdSIG8BNwDxgqzzvicDLgIXALoVlNlqfNvb14+vI3/sPgd3yvG1JB+E1NDlgAyezLkh+CHhiTp8NfJN1AeSgBss2OxA/Oy/7WuBpgAr74kP5b2MtcECLv8/78/d5KDAjz3tWg3zHtvk9bQlcnJe5Eti6MO9FpL/11cAngTmkO/FmAgcAl+XlLmjxHdyf6/TcnC7y33je7wH8Fnh5YXueQAq0Hwde3+3jxmRNXa+ApxI7q73g8Vjxn7Ew/6uFf4iFDeYPsu4X05y6eacUlp3fYNmtSQErgK/UzXtlTv8d8JQm2/V3Oc/3W5Q72uJ7GaLFwbrDOnw6p68EZre5n1rWp+w6gGuAJzTIU/t+7gW2LaRvy7qzqk80WG4m8Is8/+IG8xsGjzbq/OG83Nda/H0+WjsQN1lHLd+xbZZ5Rs5/BzBYN+9/8ryTmiz7ZOD3Oc9wk+9gCfmHRoPlr8t5Gp5xT/fJbR7Tz7ciYnGD9J8W3n+ifmZELANqyz23ybpXAZ9tsOzDwKn545skqTD7rfn16xFxV5P1/kd+PVjSzAbz7yT9yttUndTh6Pz6mYj4fQd16MSpEfFIg/R/Jp2RbQ+8qpD+6pz2KCn4rSci1pB+FQO8XNLTK6rn9/LrS1vk+WFEXFNFYZLeCfwV6QfT4flvuDZv91yPh0gBZgMRcTfpjA7W//6KvhARDzWZd19+3blk1acF36o7/fy6Sfod+fVh1gWJercDe5B+kTUyEREPNpl3UX7dAdiNdIkHUrsEwEn5n72VAeCphboWy129kWVb2aQ65MbX2oH1/A7K79SFjRIj4j5JvwIOBPYjtQuR3wNclQ+QjVxMupyzRc7f1vZJ2gY4jnT5bi/S30r9ceQZLVZxaTvltFGPg4DP5Y/vjoiL67LU9vlWwO/W/z2znu3y665N5req7/mky1+fkrQH8G3gshbBZlpx8Jh+bm2SXutpfXvkc+4WebZsMr/VL+/ivFmsCx61X2VPytPGDDRIW9nGcq1sah2eVkhb3mEdOtHO9z6rkDarbt4GIuJhSXeStnFWs3xFknYmBbJnF5IfJDXEryVdDtuRdNmsmU73JZLmkA7UWwJfiogvN8hW2+czWX8/NtPo7w5a1/dTwAuBPyPdhHI8sFrSIlIb279GxD1tlN2XfNnKqtLsp13tb2xeRKiNaWmDdXQ6xMim1qHpz9Ue0qqOT6i4rM+SAsdNpDvAnhIR20XEThHxdODFbayjo30paQD4LingXcS6kQXq1fb5r9rc58eWrW9EPBIR80hnfp9mXQN87fNvJO1TeiP7hIOHldHqckTxum/x19rt+XWv6qvTtk2tw22F94MV1WVTtPO9F7/z2vumdVbq8/PUBss2y78V6VIVwEhEnN3gklg7v/A79TVgX1LfnsMj4rEm+Wr7fA9Jk3qFJSIui4j3R8SBpMt4byGdqc4CvjKZZXeTg4eV8aL8y6+RV+TXeyh0kmPdNeM3TVal2rBJdchnILUAcliJRdfW3qjFxfYSXtEoUdITWde+cUVhVu39HpJmN1nnQay7bH1FkzxFO7LuTOZXTfK8so31bDJJ84E/J10qmxcRf2iRvbbPtyPdQDAlIuLBiDgLGM1JL2zWEbHfOXhYGQM0uEwg6QmkAQkBvl3XpvL1/DosqeW4RZKaNdR3qpM6/Ht+/esWB+J69xXe79DmMq38df7lX+9E0m3S95H62tT8OKdtSeoPsp58N9mH88dfRMRt9XkauI90SQbgeQ3WuTPw7jbWs0kkvQ74WK7DMRFxdav8EXED6TISpAbtpgdwSdvkv+GydWq0T2pqjeYiNdpPOw4eVsa9wMclvTffdYOkZ5KuQe9JupPrk8UFIuJHpF7XAAslfTQfaMjLP1nSPEnfJd16WrkO6/ApUsPzjsAvJP1Z7aAhaTtJcyWdJWmXQnn3kHoeA7ytgk2YA5xTG3pD0oCk97Fu1N5PRcSqQvkPkjqBArxH0nxJ2+VlZwP/SercWOtAuFGRxgyrHYwXSto3r2+GpD8mtT9MShuRpGcDZ5KOVx+NiLM3skjNu0mdP59L2nevrF3CyvXeW9KHSH05NuV222sk/aOkFxX+JqQ0VM9pOc+iFne89bdudzTx1P5Ee50Ej22y7Nw8f2mL9TdcB+sPE3I26zp73c36QzMc2WS927Kuh3dtuof1hwfZoHNZodyvb+R7Gaqto0WeTapDXvZ5wM2FPPXbvkGHOuCjhXkPkK7RLwVObHNfDxWWbzU8ybk0Hp5kJhsOT1IcOmMNcHyTsptt0wGkvj7F7ap9vpPUJtJwP2zs77NVPuDYQpm303qokl3r1vca1h++5BHgD6w/bE+wYQfDht9BXZ7ielfn76C43pXA87t5zJjMyWceVkYAR5AuUV1POh2/G/g+8JJI13o3XChdB34DaViLs0m/5LfJyy8mddA7nHSr4+RUvIM6RMSvgb1Jv9InSJcktibddXQuqYF0Rd1iHwPeTxo7SaTG60E24TJWRHwHOBj4Aemgv5o0zMq7gTdGgz4wEbEmIt6at+nHpAPddqRbuf8T2D8ivliyHv9HupPoXNJ+35LUJ+fLpEbsq8pu2ybYidQw32xar5NpRPyQdIfYP5Dadh4m7YP7SEPy/D2wZxQ6GJYwj9Th9n9JZ5rbkYLH1aQz8L1jI5fX+lltbBqzppRGgv0I8I1ofkujVShfovodQET0wy3DtpnxmYeZmZXm4GFmZqU5eJiZWWkOHmZmVtpm02C+4447xtDQULerYWbWVy6//PI/RMQGg2duNqPqDg0NMTEx0e1qmJn1FUkNb2P2ZSszMyvNwcPMzEpz8DAzs9IcPMzMrDQHDzMzK83Bw6xLxsfHGRoaYsaMGQwNDTE+Pj7lec02WbeH9Z2q6YUvfGGYbYozzzwzBgcHQ1IMDg7GmWee2XHeM888MwYGBopDgsfAwEDD/JOVd7K2zaYXYCIaDUnfKHE6Tg4etikm68A9ODi4Xr7aNDg4OGV5eyUoWW9z8HDwsDrtHOAm68AtqWFeSVOWt1eCkoNMb3PwcPDYLFR92ajfDty9EMDarYPPZvqDg4eDx7Q3GZeNeuXXeT9dOms30EzW91XL70BTjb4JHsChwI2kR4Oe3GD+3wJX5uka0mM5n7Kx9Tp4TH+TcYDrlXaBfmq0b3c/TNaZmgNNtfoieJCeP7wEeCbp2dJXAXu1yP864GftrNvBo3+1+889mQejyTjI94JuBqXJusTmQFOtfgkeBwIXFD5/APhAi/z/AfxlO+t28OhPk3UZpuxBw8pp5yA7WfvWgaZa/RI8Dge+Uvh8NPCFJnkHgLtocckKGAUmgIk5c+ZU+oXa1PA/9/Q2GZfYHGiq1S/B44gGweO0JnnfDHyv3XX7zKO3TMalqDLrtf7jQFPue6hKvwSPti9bAecAR7W7bgeP3jFZ/9xmNQ401emX4LEFcBOwG+sazPdukO9J+ZLVtu2u28Gjd/T6P4ttXqZzoKniDKUvgkeqJ4cBvyHddTU/px0HHFfIcyxwVpn1Onj0Dl+Ksn7VL4Gmyh9dfRM8Jmty8OgdvhRlm4NuBpoq/8ccPBw8Jt1k/LOYbQ66OazOxjh4OHhMql6/Y8Rsumjnf2cqzjyU5k1/w8PDMTEx0e1qTFtDQ0MsW7Zsg/TBwUGWLl069RUy24yNj48zOjrKqlWrHk8bGBhgbGyMkZGRUuuSdHlEDNen+0mCVonly5eXSjezyTMyMsLY2BiDg4NIYnBwcJMCRys+87BK+MzDbHrymYdNqgULFjAwMLBe2sDAAAsWLOhSjcxsMjl4WCWm4jTZzHqHg4e1ND4+ztDQEDNmzGBoaIjx8fGmeUdGRli6dClr165l6dKlDhxm09gW3a6A9a76OzaWLVvG6OgogAOD2WbOZx7W1Pz589e71Q9g1apVzJ8/v0s1MrNe4eBhTfn2WzNrxsHDmpozZ06pdDPbfDh4WFO+/dbMmnHwsKZ8+62ZNeMe5mZm1pR7mNvjyvTdMDNrxP08NjPuu2FmVfCZx2bGfTfMrAo9FzwkHSrpRkmLJZ3cJM9cSVdKulbSRVNdx37mvhtmVoWeCh6SZgKnA68B9gLeImmvujw7AF8E/iwi9gaOmOp69jP33TCzKvRU8AD2BxZHxE0R8ShwFjCvLs9RwNkRsRwgIu6Y4jr2NffdMLMq9FrwmA3cXPi8IqcVPRt4sqQLJV0u6ZhmK5M0KmlC0sTKlSsnobr9x303zKwKvXa3lRqk1XdE2QJ4IfDHwDbApZIui4jfbLBgxBgwBqmfR8V17VsjIyMOFmbWkV4LHiuAXQufdwFuaZDnDxHxIPCgpIuBfYANgoeZmU2OXrtstQjYQ9JukrYCjgTOq8vzXeDlkraQNAAcAFw/xfU0M9us9dSZR0SslnQCcAEwE1gYEddKOi7PPyMirpf0I+BqYC3wlYi4pnu1NjPb/HhsKzMza8pjW5mZWWUcPKYJD3ZoZlOpp9o8bNN4sEMzm2o+85gGPNihmU01B49pwIMdmtlUc/CYBjzYoZlNNQePacCDHZrZVHPwmAY82KGZTTV3EjQzs6bcSdDMzCrj4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWWkdDskv6b+BXeboSuD4i1lZQLzMz62GdPs/jduBPgRMBAY9I+jUpkNSCytUR8VCH5ZiZWQ/p6LJVRBwVEXsCTwJeCawA9gNeC3wWuAS4T9J1ktp6tJ2kQyXdKGmxpJMbzJ8r6V5JV+bp7zvZhl7nJwSaWS+q5EmCEfGgpNcDS4ADIuJuSVsB84APAc8iBZaWJM0ETgdelfMvknReRFxXl/UXEfHaKurey/yEQDPrVVU2mB8LnBoRdwNExKMR8S1gGPg/4OI21rE/sDgiboqIR4GzSAFos+QnBJpZr6oyeKwCnl6fGBGPAacCb2tjHbOBmwufV+S0egdKukrSDyXt3WxlkkYlTUiaWLlyZRvF9xY/IdDMelWVwWMc+KikRgf7rYGd21iHGqTVjxl/BTAYEfsApwHnNltZRIxFxHBEDM+aNauN4nuLnxBoZr2qyuDxIeAm4EZJ/5IbvveS9Abg08Cv21jHCmDXwuddgFuKGSLivoh4IL8/H9hS0o6VbEGP8RMCzaxXVRY8IuJh4E+AjwB/DpxPChjfyVlG21jNImAPSbvlBvcjgfOKGSQ9XZLy+/1J23BnJRvRY/yEQDPrVZP2JEFJewLPIB3Yr4mI1W0udxjpNt+ZwMKIWCDpOICIOEPSCcA7gdXAQ8D7IuKSja3XTxI0Myuv2ZMEKw0ekg4C3g4MAvcC74mInmjddfAwMytv0h9DK2kEuJDUZnEN8DpgxzzvFEmvqKosMzPrrqobzD8TEX8MnMT6d049RhrCxMzMpoEqg8cQ8KP8vv5a2DWkYUvMzGwaqDJ4LAGe12Te3UD/dbQwM7OGqgwenwc+LKnRGcYeQP918TYzs4aq7OcxBlxEGkn3G6RLV0/Pd2DNB35SVVlmZtZdVT9J8HDgb4CDSA3m3yPdgbWSFEDMzGwaqPJW3eOB7SPiCxExB9gbeDXwXODFEXF7VWWZmVl3VXnmcRrwR7UPEXF9RPx3fhbHlpK2qbAsMzProo6Ch6Sn5DGooPGIuDUvAO7rpCwzM+sdnZ55vBt4QNK1pAbyYyS9VtIudfm2AdZ0WJaZmfWITh9DuxD4HbAPsCdwFHA8EJLuBq4EfgO8GKh/lKyZmfWpjoJHRNwM/BuApCNId1stIQWTfYB9gRfiu63MzKaVjoKHpB0i4h6AfIdVzc/yZGZm01CnbR4rJe0DIOkDkl7T5DG0ZmY2jXTa5vFm4I78/r3ATqxr77iqbro2Ih7rsDwzM+sBnbZ5nF14/3RJs1jX3vF84BDgXcCWpGHZn9BJeWZm1hs6PfNYT0SsBH6aJwAkbQHsRQomZmY2DVQ9ttUGImJ1RFwdEWe2k1/SoZJulLRY0skt8r1I0hpJh1dX26kxPj7O0NAQM2bMYGhoiPHx8W5XycyslMrOPHLD+b+Q+nvcQ3oA1FWkvh5X5dt6N7aOmcDpwKuAFcAiSeflIU7q830KuKCq+k+V8fFxRkdHWbVqFQDLli1jdHQUgJGRkW5WzcysbVWeefwbsBvwRdJouk8ktXecByyVdGcb69gfWBwRN0XEo8BZwLwG+d4NfId1jfV9Y/78+Y8HjppVq1Yxf767wZhZ/6iyzWN34M0R8YNioqSdSJ0F22nzmA0Uz1BWAAfUrW828AZSY/yLWq1M0igwCjBnzpxWWafM8uXLS6WbmfWiKs88FgHb1SdGxB0R8eOI+Ewb62g0uGL989A/C7w/IjY6VlZEjEXEcEQMz5rVG0/BbRbEeiW4mZm1o8rg8QHgZElP6mAdK4BdC593AW6pyzMMnCVpKWk4lC9Ken0HZU6pBQsWMDAwsF7awMAACxYs6FKNzMzKqzJ4/BK4FLhG0vvz3VBbl1zHImAPSbvlod6PJLWZPC4idouIoYgYAr4NHB8R53Ze/akxMjLC2NgYg4ODSGJwcJCxsTE3lptZX1FE/VWhTVyR9BXg7cD9pM6AW5GGYf8N6Y6rX7Vz6UrSYaRLUzOBhRGxQNJxABFxRl3erwPfj4hvb2y9w8PDMTExUWKLzMxM0uURMbxBeoXB437gsxHx4Xwr7V6sG1l3X+D5EbFTJYVtAgcPM7PymgWPKu+2uhu4CCA3Zv86T211DjQzs/5RZZvHl4HXV7g+MzPrUVUGjy2B10r6qKQNbtk1M7Ppo8rg8VfAHODDwB2SLpL0OUlvk/QCSVtWWJaZmXVRZW0eEbFzgyHZXwEch4dkNzObVjwku5mZlVblqLqjwOqIWFhMj4jVwNV5MjOzaaDKNo/3AI8PICVprqS7JN0v6RMVlmNmZl1WZfAYAq4ofP4EsAQ4BXinpLdXWJaZmXVRlW0eDwFrASQNkoZSf0lEXCZpLenZHgtbLG9mZn2iyjOPq0hPAAQYAW6PiMvy5yuAPSosy8zMuqjKM48FwPmShoGXAqcV5j0NeLTCsszMrIuq7Ofxc0mHAG8BJoCPFmbPJY2ua2Zm00DV/TwuJT3To95M4JtVlmVmZt1TafBoJiL+airKMTOzqVFlJ8EnA38LPBtYCVxDakS/KiLur6ocMzPrvirPPM4EDgF+Rnr40wiwHRD5eeO/iojDKyzPzMy6pMrg8QrgXcXhSSQ9ixRIXoDHtjIzmzaq7OdxK/D7YkJELI6Ib0fE/Ih4XTsrkXSopBslLZZ0coP58yRdLelKSROSXlZR/c3MrE1VBo8vAW/qZAX52eenA68hjcT7Fkl71WX7b2CfiNgXeDvwlU7KNDOz8qoMHlsDh0j6B0lP3MR17A8sjoibIuJR4CxgXjFDRDwQEZE/bgsEZmY2paoMHicCzwQ+SHqS4MWSPi/p7ZL2k7RVG+uYDdxc+Lwip61H0hsk3QD8gHT2YWZmU6jKHuY7SdqZ1DBem+aSHk/b7pME1WjVDco6BzhH0kHAx4FXNlxZesbIKMCcOXPa2g4zM9u4qnuY30pqOL+glpafXb4X8Lw2VrEC2LXweRfglhblXSxpd0k7RsQfGswfA8YAhoeHfXnLzKwiVV62aigiHouIqyLizDayLwL2kLRbvsx1JHBeMYOkZ0lSfr8fsBVwZ9X1NjOz5io985D0HNJzPB4ArgNuKDRub1RErJZ0AunMZSawMCKulXRcnn8G6Y6uYyQ9RnqGyJvLlGFmZp1TVcddSW8m9TIvtm08TBqm5EpSD/MzKilsEwwPD8fExES3ijcz60uSLo+I4fr0Ki9bfYjU12N7UsP33wLfILV1HA6cWmFZZmbWRVUGj92Bc8mPogUuiojjgWHgLuDFFZZlZmZdVGXwuI90GWxtfv9UgIi4Dvgs8MkKyzIzsy6qMnj8mnT2AXATsF9h3vXAyyssy8zMuqjqsa12zu+/Bpwo6RBJs4HjgHsqLMvMzLqoklt1c7+L70fE2Tnpy8CfAj8l9RBfC/xFFWWZmVn3VdXPYwZwv6R5EfGjPKjhoZIOJPUYvyIiFldUlpmZdVklwSMi1ki6o0H6pcClVZRhZma9o8o2j4XAMRWuz8zMelTVY1u9QtKYpKdWvN6eNz4+ztDQEDNmzGBoaIjx8fFuV8nMbNJUObbVKPA0UsP42yRNABPkoUmAa3JbyLQzPj7O6Ogoq1atAmDZsmWMjo4CMDIy0s2qmZlNisrGtgKQNAvYJ0/Pz697kp/nEREbe57HpJnMsa2GhoZYtmzZBumDg4MsXbp0Uso0M5sKzca2quzMQ9LxwHhE/JR0i24tfQvS8zyeX1VZvWb58uWl0s3M+l2VbR6nAX9UnxgRq4EbgO9UWFZPafaUQj+90Mymq46Ch6SnFJ5N3ugRsjUvII13NS0tWLCAgYGB9dIGBgZYsGBBl2pkZja5Oj3zeDfwgKRrST3Jj5H0Wkm71OXbBljTYVk9a2RkhLGxMQYHB5HE4OAgY2Njbiw3s2mrowZzSbsCB5Maxk8ijV+1AymQ3E260+o35OHYI2K/BquZEn4YlJlZeZPSYB4RNwP/lgs4gvTQpyWsu+NqX+CFwEpgfidlmZlZ76jsbquIKLYO/yxPZmY2DW20zUPSRO41fpykF0ma1L4akg6VdKOkxZJObjB/RNLVebpE0j6TWR8zM9tQO2ce+5HulnpH/rxG0vXAFcDl+fXKiFjVaWUkzQROB14FrAAWSTovP42w5nfAKyLibkmvAcaAAzot28zM2tdO8HgHKYDsR+roty3wvDzVBkJcK+m3rAsmVwC/ioiyt+fuDyyOiJsAJJ0FzAMeDx4RcUkh/2VA/Z1dZmY2yTYaPCLia6QnA9Ye+vQc1gWT/UiN4k/K6c8BjqotK2lJRDy7RH1mAzcXPq+g9VnFO4AfNpspaZQ05pY77JmZVahUg3mk+3qvz9Pjw8ZK2p31A8oLgB1Z90zzdjXqaNjwXmJJB5OCx8ta1HeMdFmL4eHh6gbxMjPbzFX1MKglpFt0v1VLy31AXlByVStITx6s2QW4pT6TpOcDXwFeExF3lq6wmZl1pMoh2deT+4DcvNGM61sE7CFpN+D3wJEULoMBSJoDnA0cHRG/qaKuZmZWzqQFj00REaslnQBcAMwEFkbEtZKOy/PPAP4eeCrwxdQEw+pGvR/NzGzyVPo8j17m4UnMzMprNjxJ1Y+hNTOzzYCDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV5uBhZmalOXiYmVlpDh5mZlaag4eZmZXm4GFmZqU5eJiZWWkOHmZmVpqDh5mZlebgYWZmpTl4mJlZaQ4eZmZWmoOHmZmV1nPBQ9Khkm6UtFjSyQ3mP0fSpZIekfQ33aijmdnmbotuV6BI0kzgdOBVwApgkaTzIuK6Qra7gPcAr5/6GpqZGfTemcf+wOKIuCkiHgXOAuYVM0TEHRGxCHisGxU0M7PeCx6zgZsLn1fktE0iaVTShKSJlStXdlw5MzNLei14qEFabOrKImIsIoYjYnjWrFkdVMvMzIp6LXisAHYtfN4FuKVLdTEzsyZ6LXgsAvaQtJukrYAjgfO6XCczM6vTU3dbRcRqSScAFwAzgYURca2k4/L8MyQ9HZgAtgfWSjoR2Csi7utWvc3MNjc9FTwAIuJ84Py6tDMK728jXc4yM7Mu6bXLVmZm1gccPMzMrDQHDzMzK83Bw8zMSnPwMDOz0hw8zMysNAcPMzMrzcHDzMxKc/AwM7PSHDzMzKw0Bw8zMyvNwcPMzEpz8DAzs9IcPMzMrDQHDzMzK83Bw8zMSnPwMDOz0hw8zMysNAcPMzMrreeCh6RDJd0oabGkkxvMl6TP5/lXS9pvsuoyPj7O0NAQM2bMYGhoiPHx8ckqysysr2zR7QoUSZoJnA68ClgBLJJ0XkRcV8j2GmCPPB0AfCm/Vmp8fJzR0VFWrVoFwLJlyxgdHQVgZGSk6uLMzPpKr5157A8sjoibIuJR4CxgXl2eecC/RXIZsIOknauuyPz58x8PHDWrVq1i/vz5VRdlZtZ3ei14zAZuLnxekdPK5gFA0qikCUkTK1euLFWR5cuXl0o3M9uc9FrwUIO02IQ8KTFiLCKGI2J41qxZpSoyZ86cUulmZpuTXgseK4BdC593AW7ZhDwdW7BgAQMDA+ulDQwMsGDBgqqLMjPrO70WPBYBe0jaTdJWwJHAeXV5zgOOyXddvRi4NyJurboiIyMjjI2NMTg4iCQGBwcZGxtzY7mZGT12t1VErJZ0AnABMBNYGBHXSjouzz8DOB84DFgMrALeNln1GRkZcbAwM2ugp4IHQEScTwoQxbQzCu8DeNdU18vMzNbptctWZmbWBxw8zMysNAcPMzMrzcHDzMxKU2p/nv4krQSWdbseHdgR+EO3KzEJvF39Z7pum7erscGI2KCX9WYTPPqdpImIGO52Parm7eo/03XbvF3l+LKVmZmV5uBhZmalOXj0j7FuV2CSeLv6z3TdNm9XCW7zMDOz0nzmYWZmpTl4mJlZaQ4ePU7SUkm/lnSlpIlu16cTkhZKukPSNYW0p0j6iaTf5tcnd7OOm6LJdp0i6fd5v10p6bBu1nFTSNpV0s8lXS/pWknvzel9vc9abNd02GdbS/qlpKvytn00p1e+z9zm0eMkLQWGI6LvOy9JOgh4gPQM+ufmtE8Dd0XEJyWdDDw5It7fzXqW1WS7TgEeiIjPdLNunZC0M7BzRFwh6YnA5cDrgWPp433WYrv+nP7fZwK2jYgHJG0J/A/wXuCNVLzPfOZhUyYiLgbuqkueB3wjv/8G6Z+4rzTZrr4XEbdGxBX5/f3A9cBs+nyftdiuvhfJA/njlnkKJmGfOXj0vgB+LOlySaPdrswkeFrtSZD5dacu16dKJ0i6Ol/W6qtLO/UkDQEvAP6PabTP6rYLpsE+kzRT0pXAHcBPImJS9pmDR+97aUTsB7wGeFe+RGK970vA7sC+wK3AqV2tTQckbQd8BzgxIu7rdn2q0mC7psU+i4g1EbEvsAuwv6TnTkY5Dh49LiJuya93AOcA+3e3RpW7PV+Drl2LvqPL9alERNye/4nXAv9Kn+63fN38O8B4RJydk/t+nzXarumyz2oi4h7gQuBQJmGfOXj0MEnb5gY9JG0LvBq4pvVSfec84K35/VuB73axLpWp/aNmb6AP91tufP0qcH1E/HNhVl/vs2bbNU322SxJO+T32wCvBG5gEvaZ77bqYZKeSTrbgPS8+f+IiAVdrFJHJP0nMJc0RPTtwEeAc4FvAnOA5cAREdFXjc9Ntmsu6fJHAEuBv6pdc+4Xkl4G/AL4NbA2J3+Q1D7Qt/usxXa9hf7fZ88nNYjPJJ0cfDMiPibpqVS8zxw8zMysNF+2MjOz0hw8zMysNAcPMzMrzcHDzMxKc/AwM7PSHDys50jaTVJIWlhIOzCnfaSbdWtHHgk58tAXU1VmSPKtkzZlHDysF704v15aSHtJgzQz65Itul0BswZqweOSQtqBpM5b/7dhdgP27HYFbPPi4GG96MXAvcB1hbQDgesi4t7uVKm3RcQN3a6DbV582cp6iqStSENEXBZ5+ANJg8AzgMs2cZ2n5DaBU3J7ypmSbpf0cH7a2l9LavhDSsnRki6UdHdeZomk0yXtWrIeg5I+kJ9id7OkRyTdlT8f1WSZubnuF0oakPQPkm6Q9FAedruWb4M2D0lfr6W3mC5sUOaukj4n6cZczn2S/lfSsXlcqPr8F+Z1zZV0kKQfSPqDpLWSXp/zbC3pZElXSHogb/utki7N27R1me/Sus9nHtZ1TRp6/6RB+jskvSO/XxYRQyWL2g2YAB4mjTa6PXAw8BngZZLelEdUrdVLwJnAUcBjeZm7SKOtHg8cKenQiFjUZvlHAx8HlpAGq/tf0rDZLwfmSjogIt7bZNmtc/l7AhcDVwFbbaS8/2kx72Wk4cfXFBMlHUwaT+1JwGLgR8B2pLPBrwGHAMc0WecRwHGkM8afkMb6ekzSDOAHedl7gYvy69OAPwLmA18AbtvI9lgviQhPnro6AV8pTNeR2ja+XUhbntO+Vkj7ZIn1n5KXr61368K8PYAVed7xdcsdn9NvA/YupM8EPs+6AfSeULfc0jxvqC79RcX11NWhto0H1M2bW6j7r0gP9Wm0jZH+ndv6PuYCjwCriuUBO5OC42rSyKsqzNs1lx/AsXXru7BQx9EG5R2U511OekRqcZ6AlwID3f479FRu6noFPHkqTqRf4w8AWxXSlgM3drDOWvBYBezUYP7b8vzf1qUvyel/2WCZrQoH/JG6eQ2Dx0bq+Jd5mX+qSy8Gj5e0WL6t4EE6c7mbdMbxhrp5n8rr+VSTZYdrQaAuvRY8ftxkuSPy/M92++/LU3WTL1tZz5C0PemS0E8i4tGctgfpV+8ZFRTx40gP1ao3TjqbeZak2RHxe0m7AM8kDdn97/ULRMSjksaBk0kH+PF2KpCv7f8J6SxkFvCEPKv2LIlnN1n09oi4pMm8tkh6GnA+sANwUkScU5flsPz6rSaruJwU2PeVtHVEPFw3/+wGywBcQQpW75D0G+A7EXF72fpbb3HwsF4yl/Q3+dNC2iH59WcVrP93jRJzILgVmE1qg/h9fg9wa4ODZM2S/Dq7yfz1SDqQ9EyFXVpk275J+rJ2ymhR9gDwPWAI+HxEfLZBtmfm10UN2sXrPZX0PRU1rGNELJF0Eqlt6XTgdEk3kW7F/i5wTkSsabSs9S4HD+saSXOBnzeYdaqk+udHf7NwQPtGRBw7SdWqNdKr7nMjGz3CPp4xHbzPITUSf5X0vOzFwP0RsVbSq4ELWqzzoXbLalD2DOA/SGc73wVOapJ1Zn79L9JNBa080iCtaR0j4jRJ3wJeT2qsfxnw//J0paRXxDR6PvrmwMHDuuk20lPPag4nNdaemz+LdIfSbcCPC/la3UXUylCjxHx7cO2y0S35dUV+fYakJ0REo4Plbvm1/hd4IweRAsflEfEXDeY/q411bKp/AeYBvwSOisIdZXVuzvX4eERcW3UlIuI20uXHMwAk7UO6JLgv6fLfB6su0yaPg4d1TaSObccCKD13+WjS9fBa2vNIt4WeFRHvq6DIV0uaFREr69LfQurztCQiVuS6rciXVp5J+nX81eICkrYk3cILqcF4Y56SX29uMr9hP49OSToReA/pkt3rImJVi+w/BN5NauCuPHjUi4irJH2O1N60z2SXZ9VyJ0HrFXNJf4/Fy1gH59cq2jsABoAvSKo1UiNpd1LfC4DP1eX/5/z6cUnPKSwzE/g0MEi6zv/tNsqu9QA/pG5dMyT9Pel21UrlDnqnku6uOqzJzQJF/wTcB3xQ0rsadZyU9GJJR5SsxyGSDqtfX/4ea430HbXp2NTzmYf1ilrDeH3wWEPqFFeFfwf+FFgi6X9Jnd8OIXXA+x6pMbfoi6SD+luAqyT9nHQg3p90RnI3cESTS1rriYgrJH0PeB3pGv/PSR3lXgTMIQWjv+t4C9d3KikgLwNObtIIfkNEfDLX8eYccL5N6rQ3X9K1wJ2kHv6759f/ovkdWY08n3Tp7F5JVwC3kgL5AaTLhbeRbhO2PuLgYb3iYOC2fCmr1rv7IGCiwobUm0gH638kBY0n5bSFpD4I67UFRERIGiFdzvlLUi/rbUjtIl8CPhERzS5DNXI4qbH6aNKZ1gOkUYKPyuutOnjUGsD3zVMjFwGfrH2IiJ9L2pt0+epPSdu8JekA/1vgNMoFDkiBeQfS/nwWaYTkB0j9ZM4AvtTgUqL1OEX4EQA2vUk6BfgI8NGIOKW7tTGbHtzmYWZmpTl4mJlZaQ4eZmZWmts8zMysNJ95mJlZaQ4eZmZWmoOHmZmV5uBhZmalOXiYmVlp/x9C86qBZKpfNAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the result\n",
    "plt.figure()\n",
    "plt.plot(Narray, I, 'ko')\n",
    "plt.xlabel('# polarizers', fontsize=22)\n",
    "plt.ylabel('$I_{transmited}$', fontsize=22)\n",
    "plt.title('Imperfect polarizers', fontsize=24)\n",
    "\n",
    "# Find the maximum\n",
    "Imax = np.max(I)\n",
    "Nmax = np.argmax(I) + 2\n",
    "print('The maximum obtained intensity is {} for {} polarizers.'.format(Imax, Nmax))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.7 64-bit",
   "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.7"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": true,
   "user_envs_cfg": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": true
  },
  "vscode": {
   "interpreter": {
    "hash": "26de051ba29f2982a8de78e945f0abaf191376122a1563185a90213a26c5da77"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}