{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example 5: Optical elements polarimetry\n",
"\n",
"Optical elements polarimetry consists on measuring the influence of an optical element upon a light wave polarization, i.e., its Jones or Mueller matrix of the optical element.\n",
"\n",
"There are different ways of implementing a polarimeter for optical elements. In this example we will describe one of the common ones: the four-axis polarimeter, a polarimeter composed of a circularly-polarized light source, two rotating quarter-wave plates (linear retarders with 90º retardance), two rotating perfect polarizers and a photodetector:\n",
"\n",
"\n",
"\n",
"Most optical element polarimeters, like this one, are usually divided into two parts: the polarization states generator (PSG) and the polarization states analyzer (PSA). The polarization states generator is the part of the polarimeter capable of generating many polarization states. In our case, it is composed by the first rotating polarizer and waveplate. If both of them are perfect elements, this PSG can generate all the possible polarization states. If the incident light wave is circularly polarized, all the generated states will have the same intensity. The second part is the polarization states analyzer (PSA). This part is the same as the light polarimeter described in the previous example. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.optimize import least_squares\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from py_pol import degrees\n",
"from py_pol.jones_vector import Jones_vector\n",
"from py_pol.jones_matrix import Jones_matrix, create_Jones_matrices\n",
"from py_pol.stokes import Stokes, create_Stokes\n",
"from py_pol.mueller import Mueller, create_Mueller"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Jones formalism\n",
"\n",
"The most general Jones matrix can be described using 4 parameters. One selection of parameters is:\n",
"\n",
"$J=\\exp\\left(i\\Phi\\right)\\left[\\begin{array}{cc}\n",
"J_{0} & J_{1}\\exp\\left(i\\delta_{1}\\right)\\\\\n",
"J_{2}\\exp\\left(i\\delta_{2}\\right) & J_{3}\\exp\\left(i\\delta_{3}\\right)\n",
"\\end{array}\\right]$,\n",
"\n",
"where the modules $J_i$ are real and positive and the phases $\\delta_i$ and the global phase $\\Phi$ are real and range between 0 and 360º. In principle, it is not possible to measure the global phase using an intensity detector without using interferometric techniques, so first we will concentrate in measuring the other parameters:\n",
"\n",
"$J=\\left[\\begin{array}{cc}\n",
"J_{0} & J_{1}\\exp\\left(i\\delta_{1}\\right)\\\\\n",
"J_{2}\\exp\\left(i\\delta_{2}\\right) & J_{3}\\exp\\left(i\\delta_{3}\\right)\n",
"\\end{array}\\right]$.\n",
"\n",
"There are several techniques to determine the 7 remaining parameters of the Jones matrix. Here we are going to implement the method described in \"Calibration method to determine the complete Jones matrix of SLMs\", J. Hoyo, L.M. Sanchez-Brea and A. Soria-Garcia, Optics and Lasers in Engineering, Volume 151, April 2022, 106914. This method requires performing 10 intensity measurements, varying the rotation angles of the elements of the PSG and PSA, which we can characterize using the bra-ket notation [1]:\n",
"\n",
"1. $I_{1}=\\left|\\left\\langle 0\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=J_{0}^{2}$,\n",
"1. $I_{2}=\\left|\\left\\langle 0\\text{º},0\\text{º}\\right|J_{sample}\\left|90\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=J_{1}^{2}$,\n",
"1. $I_{3}=\\left|\\left\\langle 90\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=J_{2}^{2}$,\n",
"1. $I_{4}=\\left|\\left\\langle 90\\text{º},0\\text{º}\\right|J_{sample}\\left|90\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=J_{3}^{2}$,\n",
"1. $I_{5}=\\left|\\left\\langle 0\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},-45\\text{º}\\right\\rangle \\right|^{2}=\\left(J_{0}^{2}+J_{1}^{2}+2J_{0}J_{1}\\sin\\delta_{1}\\right)/2$,\n",
"1. $I_{6}=\\left|\\left\\langle 0\\text{º},0\\text{º}\\right|J_{sample}\\left|45\\text{º, 0º}\\right\\rangle \\right|^{2}=\\left(J_{0}^{2}+J_{1}^{2}+2J_{0}J_{1}\\cos\\delta_{1}\\right)/2$,\n",
"1. $I_{7}=\\left|\\left\\langle 0\\text{º},45\\text{º}\\right|J_{sample}\\left|0\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=\\left(J_{0}^{2}+J_{2}^{2}+2J_{0}J_{2}\\sin\\delta_{2}\\right)/2$,\n",
"1. $I_{8}=\\left|\\left\\langle 45\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=\\left(J_{0}^{2}+J_{2}^{2}+2J_{0}J_{2}\\cos\\delta_{2}\\right)/2$,\n",
"1. $I_{9}=\\left|\\left\\langle 90\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},-45\\text{º}\\right\\rangle \\right|^{2}=\\left[J_{2}^{2}+J_{3}^{2}+2J_{2}J_{3}\\sin\\left(\\delta_{3}-\\delta_{2}\\right)\\right]/2$,\n",
"1. $I_{10}=\\left|\\left\\langle 90\\text{º},0\\text{º}\\right|J_{sample}\\left|45\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=\\left[J_{2}^{2}+J_{3}^{2}+2J_{2}J_{3}\\cos\\left(\\delta_{3}-\\delta_{2}\\right)\\right]/2$,\n",
"\n",
"where $\\left|\\phi,\\chi\\right\\rangle$ is the state generated by the PSG, $\\left\\langle \\phi,\\chi\\right|$ the state analyzed by the PSA (the state with maximum transmission), and $\\phi$ and $\\chi$ are the azimuth and ellipticity angle. Then, 7 parameters can be calculated from the intensity measurements:\n",
"\n",
"1. $J_{0}=\\sqrt{I_{1}}$,\n",
"1. $J_{1}=\\sqrt{I_{2}}$,\n",
"1. $J_{2}=\\sqrt{I_{3}}$,\n",
"1. $J_{3}=\\sqrt{I_{4}}$,\n",
"1. $\\delta_{1}=\\arctan\\left(\\frac{2I_{5}-I_{1}-I_{2}}{2I_{6}-I_{1}-I_{2}}\\right)$,\n",
"1. $\\delta_{2}=\\arctan\\left(\\frac{2I_{7}-I_{1}-I_{3}}{2I_{8}-I_{1}-I_{3}}\\right)$,\n",
"1. $\\delta_{3}=\\arctan\\left(\\frac{2I_{9}-I_{3}-I_{4}}{2I_{10}-I_{3}-I_{4}}\\right)+\\delta_{2}$.\n",
"\n",
"Let's simulate the experiment. The first stepp is implementing a function that calculates the intensity depending on the elements angles and another to calculate the matrix parameters from the intensities."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def sample_polarimeter_measurement(Jsample, angleDG, angleRG, angleRA, angleDA):\n",
" \"\"\"This function simulates the light polarimeter measurements.\n",
" \n",
" Parameters:\n",
" E: Incident light wave Jones vector.\n",
" angleDG: Rotating angle of the diattenuator (polarizer) of the PSG.\n",
" angleRG: Rotation angle of the quarter-wave plate of the PSG.\n",
" angleRA: Rotation angle of the quarter-wave plate of the PSA.\n",
" angleDA: Rotating angle of the diattenuator (polarizer) of the PSA.\n",
" \n",
" Returns:\n",
" I: Measured intensity.\"\"\"\n",
" # Rotate the optical elements\n",
" E = Jones_vector().circular_light(intensity=2) # Intensity = 2 produces a normalized light state of intensity 1 after the PSG\n",
" Jdg = Jones_matrix().diattenuator_perfect(azimuth=angleDG)\n",
" Jrg = Jones_matrix().quarter_waveplate(azimuth=angleRG)\n",
" Jra = Jones_matrix().quarter_waveplate(azimuth=angleRA)\n",
" Jda = Jones_matrix().diattenuator_perfect(azimuth=angleDA)\n",
" # Multiply all objects\n",
" E_final = Jda * Jra * Jsample * Jrg * Jdg * E\n",
" # Measure the intensity\n",
" I = E_final.parameters.intensity()\n",
" # Return\n",
" return I\n",
"\n",
"def matrix_parameters(I):\n",
" \"\"\"This function calculates the parameters of the Jones matrix from the intensity measurements.\n",
" \n",
" Parameters:\n",
" I: Array of intensities.\n",
" \n",
" Returns:\n",
" J0, J1, J2, J3, d1, d2, d3: Matrix parameters.\"\"\"\n",
" # Calculate modules\n",
" J0, J1, J2, J3 = np.sqrt(I[:4])\n",
" # Calculate phases\n",
" d1 = np.arctan2(2*I[4] - I[1] - I[0], 2*I[5] - I[1] - I[0]) % (360*degrees) \n",
" d2 = np.arctan2(2*I[6] - I[2] - I[0], 2*I[7] - I[2] - I[0]) % (360*degrees) \n",
" d3 = (np.arctan2(2*I[8] - I[2] - I[3], 2*I[9] - I[2] - I[3]) + d2) % (360*degrees) \n",
" # Return\n",
" return J0, J1, J2, J3, d1, d2, d3 "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can check if it works."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Comparison between original and calculated parameters:\n",
"- J0: Original = 0.0666; Calculated = 0.0666\n",
"- J1: Original = 0.6734; Calculated = 0.6734\n",
"- J2: Original = 0.6419; Calculated = 0.6419\n",
"- J3: Original = 0.5885; Calculated = 0.5885\n",
"- d1 (deg): Original = 134.6173; Calculated = 134.6173\n",
"- d2 (deg): Original = 181.8143; Calculated = 181.8143\n",
"- d3 (deg): Original = 274.9632; Calculated = 274.9632\n"
]
}
],
"source": [
"# Generate a random Jones matrix for the sample\n",
"J0, J1, J2, J3 = np.random.rand(4)\n",
"d1, d2, d3 = np.random.rand(3) * 360*degrees\n",
"Jsample = Jones_matrix().from_components([J0, J1*np.exp(1j*d1), J2*np.exp(1j*d2), J3*np.exp(1j*d3)])\n",
"\n",
"# Calculate the intensities\n",
"angleDG = np.array([0, 90, 0, 90, 45, 45, 0, 0, 45, 45]) * degrees\n",
"angleRG = np.array([0, 90, 0, 90, 0, 45, 0, 0, 0, 45]) * degrees\n",
"angleRA = np.array([0, 0, 90, 90, 0, 0, 0, 45, 90, 90]) * degrees\n",
"angleDA = np.array([0, 0, 90, 90, 0, 0, 45, 45, 90, 90]) * degrees\n",
"I = sample_polarimeter_measurement(Jsample, angleDG, angleRG, angleRA, angleDA)\n",
"\n",
"# Calculate the matrix parameters\n",
"J0_calc, J1_calc, J2_calc, J3_calc, d1_calc, d2_calc, d3_calc = matrix_parameters(I)\n",
"\n",
"# Compare results\n",
"print(\"Comparison between original and calculated parameters:\")\n",
"original = [J0, J1, J2, J3, d1/degrees, d2/degrees, d3/degrees]\n",
"calculated = [J0_calc, J1_calc, J2_calc, J3_calc, d1_calc/degrees, d2_calc/degrees, d3_calc/degrees]\n",
"parameter = [\"J0\", \"J1\", \"J2\", \"J3\", \"d1 (deg)\", \"d2 (deg)\", \"d3 (deg)\"]\n",
"for ind in range(7):\n",
" print(\"- {}: Original = {:.4f}; Calculated = {:.4f}\".format(parameter[ind], original[ind], calculated[ind]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The match is perfect."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mueller-Stokes formalism\n",
"\n",
"The method for determining the Mueller matrix of an optical element is similar to the one for Stokes light polarimetry. The k-th measurement of the light polarimeter experiment is:\n",
"\n",
"$S_{f}^{k}=\\left[\\begin{array}{cccc}\n",
"A_{00}^{k} & A_{01}^{k} & A_{02}^{k} & A_{03}^{k}\\\\\n",
"\\bullet & \\bullet & \\bullet & \\bullet\\\\\n",
"\\bullet & \\bullet & \\bullet & \\bullet\\\\\n",
"\\bullet & \\bullet & \\bullet & \\bullet\n",
"\\end{array}\\right]\\left[\\begin{array}{cccc}\n",
"M_{00} & M_{01} & M_{02} & M_{03}\\\\\n",
"M_{10} & M_{11} & M_{12} & M_{13}\\\\\n",
"M_{20} & M_{21} & M_{22} & M_{23}\\\\\n",
"M_{30} & M_{31} & M_{32} & M_{33}\n",
"\\end{array}\\right]\\left[\\begin{array}{c}\n",
"S_{0}^{k}\\\\\n",
"S_{1}^{k}\\\\\n",
"S_{2}^{k}\\\\\n",
"S_{3}^{k}\n",
"\\end{array}\\right]=\\left[\\begin{array}{c}\n",
"I^{k}\\\\\n",
"\\bullet\\\\\n",
"\\bullet\\\\\n",
"\\bullet\n",
"\\end{array}\\right]$.\n",
"\n",
"Being $A^{k}=M_{D2}(\\theta_{D2}^{k})*M_{R2}(\\theta_{R2}^{k})$ the Mueller matrix of the PSA, $S^k=M_{R1}(\\theta_{R1}^{k})*M_{D1}(\\theta_{D1}^{k})*S_i$ the Stokes vector generated by the PSG and $M_{ij}$ the elements of the Mueller matrix of the sample. Due to the fact that the detector only measures the intensity, only some parts of the Mueller matrix of the PSA are important (the rest of the parameters are replaced by $\\bullet$). Then, if we consider the 16 required measurements:\n",
"\n",
"$\\left[\\begin{array}{ccc}\n",
"A_{00}^{0}S_{0}^{0} & \\ldots & A_{03}^{0}S_{3}^{0}\\\\\n",
"\\vdots & \\ddots & \\vdots\\\\\n",
"A_{00}^{16}S_{0}^{16} & \\ldots & A_{03}^{16}S_{3}^{16}\n",
"\\end{array}\\right]\\left[\\begin{array}{c}\n",
"M_{00}\\\\\n",
"\\vdots\\\\\n",
"M_{44}\n",
"\\end{array}\\right]=\\left[\\begin{array}{ccc}\n",
"b_{0,0} & \\ldots & b_{0,16}\\\\\n",
"\\vdots & \\ddots & \\vdots\\\\\n",
"b_{16,0} & \\ldots & b_{16,16}\n",
"\\end{array}\\right]\\left[\\begin{array}{c}\n",
"M_{00}\\\\\n",
"\\vdots\\\\\n",
"M_{44}\n",
"\\end{array}\\right]=\\left[\\begin{array}{c}\n",
"I^{0}\\\\\n",
"\\vdots\\\\\n",
"I^{16}\n",
"\\end{array}\\right]$.\n",
"\n",
"If $b$ is invertible ($\\det(b)\\neq0$) then $M_{ij}$ can be easily calculated as $M_{ij}=b^{-1}I'$. This condition is satisfied if there are no repeated set of four angles.\n",
"\n",
"Let us reproduce this numerically. Let us define a function that simulates the experiment and calculates the $b$ matrix. Also, we will include an error amplitude which we will use in the future."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def Mueller_polarimeter_measurement(Msample, Mpol, S_ini, angles, Ierror=0):\n",
" \"\"\"This function simulates the light polarimeter measurements.\n",
" \n",
" Parameters:\n",
" Msample (Mueller): Mueller matrix of the sample.\n",
" Mpol (list): List with the four Mueller objects of the polarimeter: Md1, Mr1, Mr2, Md2.\n",
" Sini (Stokes): Stokes object of the illumination.\n",
" angles (numpy.ndarray): 4xN array with the angles used to perform the measurement.\n",
" Ierror (float): Intensity error amplitude.\n",
" \n",
" Returns:\n",
" Ik: Measured intensities vector.\n",
" b: b matrix.\"\"\"\n",
" # Extract the information from the input parameters\n",
" N = angles.shape[1]\n",
" Md1, Mr1, Mr2, Md2 = Mpol\n",
" # Rotate the optical elements\n",
" Md1_rotated = Md1.rotate(angle=angles[0,:], keep=True)\n",
" Mr1_rotated = Mr1.rotate(angle=angles[1,:], keep=True)\n",
" Mr2_rotated = Mr2.rotate(angle=angles[2,:], keep=True)\n",
" Md2_rotated = Md2.rotate(angle=angles[3,:], keep=True)\n",
" # Multiply the objects as required\n",
" Sg = Mr1_rotated * Md1_rotated * S_ini\n",
" A = Md2_rotated * Mr2_rotated\n",
" S_final = A * Msample * Sg\n",
" # Measure the intensity\n",
" Ik = S_final.parameters.intensity() + np.random.randn(N) * Ierror\n",
" # Calculate the b matrix\n",
" Ak = A.parameters.matrix() # Extract the 4x4xN matrix of the Mueller object\n",
" Ak = Ak[0,:,:] # Extract all the first rows\n",
"# Ak = np.transpose(Ak) # Transpose so Ak rows corrrespond to different measurements\n",
" Sk = Sg.parameters.matrix() # Extract the 4xN matrix of the Stokes object\n",
"# Sk = np.transpose(Sk) # Transpose so Sk rows corrrespond to different measurements\n",
" b = np.zeros((N,16))\n",
" for ind in range(N):\n",
" b[ind,:] = np.outer(Ak[:,ind], Sk[:,ind]).flatten()\n",
" # Return\n",
" return Ik, b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we can perform the experiment and see if the calculated Mueller matrix is correct. We will use as sample the product of a random diattenuator, a random polarizer and a random diagonal diattenuator. We will use a circularly polarized light wave for simplicity."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Sample = \n",
"[+0.508 -0.018 +0.093 -0.179] \n",
"[-0.407 +0.022 -0.091 +0.230] \n",
"[-0.038 +0.026 -0.069 -0.070] \n",
"[+0.056 +0.016 +0.136 -0.004] \n",
" Calculated = \n",
"[+0.508 -0.018 +0.093 -0.179] \n",
"[-0.407 +0.022 -0.091 +0.230] \n",
"[-0.038 +0.026 -0.069 -0.070] \n",
"[+0.056 +0.016 +0.136 -0.004] \n",
"\n"
]
}
],
"source": [
"# Start by defining the polarimeter system\n",
"S_ini = Stokes('Light source')\n",
"S_ini.circular_light(intensity=1)\n",
"Mpol = create_Mueller(('Diattenuator 1', 'Retarder 1', 'Retarder 2', 'Diattenuator 2'))\n",
"Mpol[0].diattenuator_perfect(azimuth=0)\n",
"Mpol[1].quarter_waveplate(azimuth=0)\n",
"Mpol[2].quarter_waveplate(azimuth=0)\n",
"Mpol[3].diattenuator_perfect(azimuth=0)\n",
"\n",
"# Create the sample parameters\n",
"p1 = np.random.rand() * 0.5 + 0.5\n",
"p2 = np.random.rand() * 0.5\n",
"alphaD = np.random.rand() * 90*degrees\n",
"delayD = np.random.rand() * 360*degrees\n",
"\n",
"R = np.random.rand() * 180*degrees\n",
"alphaR = np.random.rand() * 90*degrees\n",
"delayR = np.random.rand() * 360*degrees\n",
"\n",
"d = [np.random.rand(), np.random.rand(), np.random.rand()]\n",
"\n",
"# Create the sample \n",
"Md, Mr, Mp = create_Mueller(N=3)\n",
"Md.diattenuator_charac_angles(p1=p1, p2=p2, alpha=alphaD, delay=delayD)\n",
"Mr.retarder_charac_angles(R=R, alpha=alphaR, delay=delayR)\n",
"Mp.depolarizer_diagonal(d=d)\n",
"\n",
"Msample = Md*Mr*Mp\n",
"Msample.name = 'Sample'\n",
"\n",
"# Now, define the rotation angles for the optical elements\n",
"angles = np.random.rand(4,16) * 180*degrees\n",
"\n",
"# Make the measurements\n",
"Ik, b = Mueller_polarimeter_measurement(Msample=Msample, Mpol=Mpol, S_ini=S_ini, angles=angles)\n",
"\n",
"# Calculate the original Mueller matrix\n",
"b_inv = np.linalg.inv(b)\n",
"Ik = np.reshape(Ik, (16,1))\n",
"M_calc = Mueller('Calculated')\n",
"M_calc.from_matrix(M=b_inv@Ik)\n",
"\n",
"# Compare the result\n",
"print(Msample, M_calc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As explained in the previous example, experiments are subject to errors. This means that the calculated Mueller matrices will be wrong. However, the error in the calculated matrix can be reduced by oversampling, i.e., taking more than the minimum 16 measurements. In this case, $b$ will be a Nx16 matrix and $I'$ a Nx1 vector. Then, $M_{ij}$ can be calculated using the pseudoinverse of $b$: \n",
"\n",
"$b^+=(b^T*b)^{-1}*b^T$.\n",
"\n",
"Let us define a function that performs one polarimetry experiment."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def Mueller_polariemter_experiment(N, Ierror):\n",
" \"\"\"This funcion simulates a polarimetry experiment with errors in the detection.\n",
" \n",
" Parameters:\n",
" N (int): Number of measurements.\n",
" Ierror (float): Intensity error amplitude.\n",
" \n",
" Returns:\n",
" dif (numpy.ndarray): Difference between the calculated and the original array.\"\"\"\n",
" # Create the sample parameters\n",
" p1 = np.random.rand() * 0.5 + 0.5\n",
" p2 = np.random.rand() * 0.5\n",
" alphaD = np.random.rand() * 90*degrees\n",
" delayD = np.random.rand() * 360*degrees\n",
"\n",
" R = np.random.rand() * 180*degrees\n",
" alphaR = np.random.rand() * 90*degrees\n",
" delayR = np.random.rand() * 360*degrees\n",
"\n",
" d = [np.random.rand(), np.random.rand(), np.random.rand()]\n",
"\n",
" # Create the sample \n",
" Md, Mr, Mp = create_Mueller(N=3)\n",
" Md.diattenuator_charac_angles(p1=p1, p2=p2, alpha=alphaD, delay=delayD)\n",
" Mr.retarder_charac_angles(R=R, alpha=alphaR, delay=delayR)\n",
" Mp.depolarizer_diagonal(d=d)\n",
"\n",
" Msample = Md*Mr*Mp\n",
" Msample.name = 'Sample'\n",
"\n",
" # Now, define the rotation angles for the optical elements\n",
" angles = np.random.rand(4,N) * 180*degrees\n",
"\n",
" # Make the measurements\n",
" Ik, b = Mueller_polarimeter_measurement(Msample=Msample, Mpol=Mpol, S_ini=S_ini, angles=angles, Ierror=Ierror)\n",
"\n",
" # Calculate the Stokes vector\n",
" b_T = np.transpose(b)\n",
" b_inv = np.linalg.inv(b_T @ b)\n",
" b_inv = b_inv @ b_T\n",
" Ik = np.reshape(Ik, (N,1))\n",
" M_calc = b_inv @ Ik\n",
" \n",
" # Calculate the difference\n",
" dif = np.linalg.norm(Msample.M.flatten() - M_calc.flatten())\n",
"\n",
" # Return\n",
" return dif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, lets do the experiment many times, changing the number of measurements and the error amplitude. "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# Define some variables\n",
"Nstat = 100 # Number of experiments done under the same condition.\n",
"I_errors = np.array([0, 0.005, 0.01, 0.02]) # Array of error factors\n",
"Ns = np.arange(20, 100, 10) # Array of the number of measurements\n",
"tol = 1e-12\n",
"I_wave = 1 # We will maintain this parameter constant so we can compare the errors easily\n",
"\n",
"# Create the result variables\n",
"error = np.zeros((Ns.size, I_errors.size))\n",
"legend = []\n",
"\n",
"# Make the loops changing the variables\n",
"for indI, Ierror in enumerate(I_errors): # Loop in the error amplitude\n",
" for indN, N in enumerate(Ns): # Loop in the number of measurements\n",
" # Create the temporal variables\n",
" aux = np.zeros((Nstat))\n",
" # Repeat the same conditions several times\n",
" for indS in range(Nstat):\n",
" # Do the experiment\n",
" dif = Mueller_polariemter_experiment(N=N, Ierror=Ierror);\n",
" # Calculate the errors\n",
" aux[indS] = dif\n",
" # Calculate the mean square error\n",
" error[indN, indI] = np.linalg.norm(aux) / Nstat\n",
" # Create the legend for the plots\n",
" legend += ['Ierror = {} %'.format(Ierror*100)]"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'MSE')"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAHgCAYAAABJrX+JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABXTUlEQVR4nO3dd5xU9b3/8ddnZmfLDAssvTepFkREYBeDsUWxRqOyaKIxJsZo6v2lx3tzTa/3GnOTGGPXUGyxRYEYCyKggAoWQJC6S6/b+/f3x5ldtrIL7OyZmX0/H4957Mw537PzOaC8z/mec75fc84hIiIiySXgdwEiIiLS/hTwIiIiSUgBLyIikoQU8CIiIklIAS8iIpKEFPAiIiJJKMXvAtpTr1693LBhw/wuQ0REpEOsXLlyr3Oud3Prkirghw0bxooVK/wuQ0REpEOY2ZaW1qmLXkREJAkp4EVERJKQAl5ERCQJJdU1eBEROXqVlZXk5eVRVlbmdynSgvT0dAYNGkQoFGrzNgp4EZFOLi8vj8zMTIYNG4aZ+V2ONOKcY9++feTl5TF8+PA2b6cuehGRTq6srIyePXsq3OOUmdGzZ8+j7mGJacCb2YVmts7MNpjZ95tZb2Z2V3T9ajObWG9ddzN7wszWmtkaM8uOZa0iIp2Zwj2+HcvfT8wC3syCwJ+AGcCJwCwzO7FRsxnAqOjrZuAv9db9AZjvnBsLnAqsiVWtIiLiry5duvhdQqucc3z9619n5MiRjB8/nrfffrvZdps2bWLKlCmMGjWKmTNnUlFR0aTNG2+8wfjx4znjjDPYsGEDAAcPHuSCCy7AOdcu9cbyDH4ysME5t9E5VwHMBS5v1OZy4GHnWQZ0N7P+ZtYVmA7cB+Ccq3DOHYxhrSIikiCcc9TU1LT4uSXV1dXH9b0vvvgi69evZ/369dxzzz185Stfabbd9773Pb71rW+xfv16srKyuO+++5q0+f3vf8+TTz7JL37xC/7yF+/c9qc//Sk//OEP2603JZYBPxDYVu9zXnRZW9qMAPYAD5jZO2Z2r5lFmvsSM7vZzFaY2Yo9e/a0X/UiIuKL3/72t5xxxhmMHz+eH//4xwBs3ryZcePGceuttzJx4kRef/31Bp+3bdvGd77zHU4++WROOeUU5s2bB8Crr77K2WefzbXXXsspp5xyXHU988wzXH/99ZgZU6dO5eDBg+zYsaNBG+ccL7/8MldddRUAN9xwA08//XST3xUKhSgtLaWkpIRQKMTHH39Mfn4+Z5111nHVWF8s76Jv7hCkcb9DS21SgInA15xzb5rZH4DvA//ZpLFz9wD3AEyaNKl9+jVERDqpO577gA+3F7Tr7zxxQFd+fOlJbWq7cOFC1q9fz1tvvYVzjssuu4xFixYxZMgQ1q1bxwMPPMCf//xnNm/e3ODzk08+ybvvvsuqVavYu3cvZ5xxBtOnTwfgrbfe4v3332/2DvSZM2eybt26Jsv/4z/+g+uvv77Bsvz8fAYPHlz3edCgQeTn59O/f/+6Zfv27aN79+6kpKQ0aNPYD37wA26++WYyMjJ45JFH+Pa3v81Pf/rTNv0ZtVUsAz4PGFzv8yBgexvbOCDPOfdmdPkTeAEvIiJJbOHChSxcuJDTTjsNgKKiItavX8+QIUMYOnQoU6dOrWtb//PixYuZNWsWwWCQvn37ctZZZ7F8+XK6du3K5MmTW3y8rPZMvy2auzbeuDu9LW0AJkyYwLJlywBYtGgRAwYMwDnHzJkzCYVC/P73v6dv375trq05sQz45cAoMxsO5AO5wLWN2jwLfNXM5gJTgEPOuR0AZrbNzMY459YB5wIfxrBWERGBNp9px4pzjh/84Ad8+ctfbrB88+bNRCINr9TW/3ykG9Mab1ff0ZzBDxo0iG3bDl9VzsvLY8CAAQ3a9OrVi4MHD1JVVUVKSkqzbepzzvGzn/2MefPm8dWvfpU77riDzZs3c9ddd/Hzn/+8xe3aImbX4J1zVcBXgQV4d8A/5pz7wMxuMbNbos1eADYCG4C/AbfW+xVfA/5uZquBCcAvYlWriIjEhwsuuID777+foqIiwOsW3717d6vbTZ8+nXnz5lFdXc2ePXtYtGgRkydPbnW7efPm8e677zZ5NQ53gMsuu4yHH34Y5xzLli2jW7duDbrnwTtbP/vss3niiScAeOihh7j88sb3lx/20EMPcfHFF5OVlUVJSQmBQIBAIEBJSUmrtbcmpiPZOedewAvx+svurvfeAbe1sO27wKRY1iciIvHlU5/6FGvWrCE72xv6pEuXLjz66KMEg8EjbnfFFVewdOlSTj31VMyM3/zmN/Tr14+1a9e2W20XXXQRL7zwAiNHjiQcDvPAAw80WHfvvfcyYMAAfv3rX5Obm8vtt9/Oaaedxk033dTs7yspKeGhhx5i4cKFgNdr8JnPfIbU1FTmzJlz3PVaez1vFw8mTZrkNB+8iMjRWbNmDePGjfO7DGlFc39PZrbSOdfsybCGqm1BTXk5VQcO+F2GiIjIMVHAN8PV1LDhnHPZ84c/+F2KiIjIMVHAN8MCATJOOYXiJUv9LkVEROSYKOBbEMnJoXLrViry8vwuRURE5Kgp4FsQmZYDQPGSJT5XIiIicvQU8C1IHTGClD591E0vIiIJSQHfAjMjkpNDydKluDbMUiQiIscumaaL/fznP8/w4cOZMGECEyZM4N13323SJtGni014kWk5VB86RNmHmopeRCRexPt0seDNiFc7Kt6ECROarE/06WITXiQ6iYGuw4uIdJxEni62rRJ9utiEl9K7N2mjR1O8dAm9bv6S3+WIiMTei9+Hne+17+/sdwrM+FWbmib6dLG1fvSjH/GTn/yEc889l1/96lekpaU1WJ/o08UmhUhODgdmz6amrIxAerrf5YiIJLVEny4W4Je//CX9+vWjoqKCm2++mV//+tf813/9V4M2iT5dbFKI5GSz/8EHKVmxki5nTvO7HBGR2GrjmXasJPp0sUDdGX1aWho33ngjv/vd71r8/oScLjZZhCdNwkIhipfqOryISKwl+nSxQN11eeccTz/9NCeffHKL35+w08Umg0A4TMZpp+l5eBGRDpAM08Ved9117NmzB+ccEyZM4O67727292m62KMQq+li9979V/bceSej3lhMSs+e7f77RUT8pOliE4Omi42BSI53JFkcvSFCREQk3ing2yD9pJMIdOum5+FFRCRhKODbwIJBIlOmULxkabsNISgiIhJLCvg2iuRkU7VjBxWbN/tdioiISKsU8G0UydH0sSIikjgU8G2UOmQIoUGD9LiciIgkBAX8UYhkZ1Py5pu4qiq/SxERSSqJMF3s2rVryc7OJi0t7Yij023atIkpU6YwatQoZs6cSUVFRZM2mi42zkSm5VBTVETpe+08EYOIiLSZX9PF9ujRg7vuuotvf/vbR2z3ve99j29961usX7+erKws7rvvviZtNF1snAlPmQJmug4vIhJD8TpdbJ8+fTjjjDMIhUIttnHO8fLLL3PVVVcBcMMNN/D00083aafpYuNMSlYW6SeeSPHSpfS+7Ta/yxERaXe/fuvXrN3ffsO7AoztMZbvTf5em9rG83SxbbFv3z66d+9OSooXr7VTyjam6WLjUCQnh30PPEB1UTHBLi3PUCQiIkcvnqeLbYu2Timr6WLjUGRaDvv+9jdKlr9F5tln+12OiEi7auuZdqzE83SxbdGrVy8OHjxIVVUVKSkpLU4pW0vTxcaRjNNOw9LSKF6qx+VERNpbPE8X2xZmxtlnn80TTzwBeNPBXn755S2213SxcSSQlkZ40iTdaCciEgPxPF3szp07mTRpEgUFBQQCAe68804+/PBDunbt2mC62F//+tfk5uZy++23c9ppp3HTTTc1+/s0XexRiNV0sY3tu+9+dv/2t4x87VVCx3mNRETEb5ouNjFoutgOUDd9rLrpRUQkTingj0HamDEEe/RQN72IiMQtBfwxsECASHY2xUs1fayIiMQnBfwxiuRkU71nL+Xr1/tdioiISBMK+GOk6WNFRCSeKeCPUah/f1KHD1fAi4hIXFLAH4dIdjYly1fgmpkKUERE2i4Rpov9+9//zvjx4xk/fjw5OTmsWrWq2XaaLjYJRKbl4EpLKXn3Xb9LERHpNPyaLnb48OG89tprrF69mv/8z//k5ptvbradpotNAuHJkyEYVDe9iEg7itfpYnNycsjKygJg6tSp5OXlNWmj6WKTRDAzk4xTTqF4yVL45jf9LkdE5Ljt/MUvKF/TvtPFpo0bS78f/rBNbRNlutj77ruPGTNmNFmu6WKTSCQnh7133031oUMEu3XzuxwRkYSWCNPFvvLKK9x3330sXry4yTpNF5tEItNy2PvnP1P85pt0/dSn/C5HROS4tPVMO1bifbrY1atX88UvfpEXX3yRnj17Nlmv6WKTSMb48QTCYV2HFxFpB/E8XezWrVu58soreeSRRxg9enSzv0/TxSYRC4UIT56siWdERNpBPE8X+5Of/IR9+/Zx6623ApCSkkLtDKaaLjbGOmq62Mb2P/wIu37xC0546SVSBw3s8O8XETkemi42MWi6WB/UTR+75A2fKxEREfEo4NtB6gknkNKnj7rpRUQkbijg24GZEcnJoWTpMlwbRlMSERGJNQV8O4nkZFN98CBlH67xuxQRkaOWTPdjJaNj+ftRwLeTSPSOz+KlelxORBJLeno6+/btU8jHKecc+/btIz09/ai202Ny7SSld2/SRo+meMkSen3pS36XIyLSZoMGDSIvL489e/b4XYq0ID09nUGDBh3VNgr4dhTJzubAnDnUlJUROMojLRERv4RCoRaHcpXEpS76dhSZloOrqKBk5Uq/SxERkU5OAd+OwpMmQSikYWtFRMR3Cvh2FAiHCU+Y4E0fKyIi4iMFfDuLTMuhfM0aqvbv97sUERHpxBTw7SySkwOgUe1ERMRXMQ14M7vQzNaZ2QYz+34z683M7oquX21mE+ut22xm75nZu2bW8TPIHKP0k04i0LWrrsOLiIivYvaYnJkFgT8B5wN5wHIze9Y592G9ZjOAUdHXFOAv0Z+1znbO7Y1VjbFgwSCRKVMoXrIU5xxm5ndJIiLSCcXyDH4ysME5t9E5VwHMBRrPen858LDzLAO6m1n/GNbUISLTcqjasYOKzZv9LkVERDqpWAb8QGBbvc950WVtbeOAhWa20sxujlmVMVA3bK266UVExCexDPjm+qYbD3R8pDbTnHMT8brxbzOz6c1+idnNZrbCzFbEyzCLoSFDCA0cqBvtRETEN7EM+DxgcL3Pg4DtbW3jnKv9uRv4B16XfxPOuXucc5Occ5N69+7dTqUfn7rpY5e9iauq8rscERHphGIZ8MuBUWY23MxSgVzg2UZtngWuj95NPxU45JzbYWYRM8sEMLMI8Cng/RjW2u4iOdnUFBVR+t57fpciIiKdUMzuonfOVZnZV4EFQBC43zn3gZndEl1/N/ACcBGwASgBboxu3hf4R/QO9BRgtnNufqxqjYXw1KlgRvHSpYRPO83vckREpJOxZJr/d9KkSW7Fivh5ZH7TZ67CMtIZ9uijfpciIiJJyMxWOucmNbdOI9nFUCQnm9J3V1FdVOx3KSIi0sko4GMokpMDVVWUrFjudykiItLJKOBjKGPiRCwtTc/Di4hIh1PAx1AgLY3w6acr4EVEpMMp4GMsMi2Hig0fU7lrt9+liIhIJ6KAj7HD08fqLF5ERDqOAj7G0saMIdijh7rpRUSkQyngY8wCASJTp1K81Js+VkREpCMo4DtAZFoO1Xv2Ur5+vd+liIhIJ6GA7wCaPlZERDqaAr4DhAYMIHXYME0fKyIiHUYB30EiOTmULF+Bq6jwuxQREekEFPAdJDItB1dSQsm77/pdioiIdAIK+A4SnjwZgkF104uISIdQwHeQYGYmGaecohvtRESkQyjgO1AkJ4ey996n+tAhv0sREZEkp4DvQJGcbKipofitt/wuRUREkpwCvgNlnHoqgXBY3fQiIhJzCvgOZKEQ4cmTFfAiIhJzCvgOFsnJpnLLViry8v0uRUREkpgCvoNp+lgREekICvgOlnrCCaT06aNuehERiSkFfAczMyLZ2ZQsXYarqfG7HBERSVIKeB9EpuVQffAgZWvW+F2KiIgkKQW8DzR9rIiIxJoC3gcpvXuTNmoUJRqXXkREYkQB75NITg4lK1ZSU1bmdykiIpKEFPA+iUzLwVVUULJypd+liIhIElLA+yQ8aRKEQuqmFxGRmFDA+yQQDhOeMIEi3WgnIiIxoID3UWRaDuUfrqFq/36/SxERkSSjgPdR7eNyJcuW+VyJiIgkGwW8j9JPPplA167qphcRkXangPeRBYNEpkyheMkSnHN+lyMiIklEAe+zSE42Vdt3ULlli9+liIhIElHA+6x2+lh104uISHtSwPssNGQIoYEDNS69iIi0KwW8z8yMSE42JW++hauq8rscERFJEgr4OBDJyaGmsJCy99/3uxQREUkSCvg4EJ46Fcx0HV5ERNqNAj4OpGRlkT5uHCVLNC69iIi0DwV8nIhMy6Fk1Spqiov9LkVERJKAAj5ORHJyoLKS4uXL/S5FRESSgAI+TmRMnIilpelxORERaRcK+DgRSEsjfPrpmh9eRETahQI+jkSm5VC+fgOVu3b7XYqIiCQ4BXwcqZ0+tnipuulFROT4KODjSNrYsQR79FA3vYiIHDcFfByxQIDI1KkUL1mq6WNFROS4KODjTCQnm6o9eyhfv97vUkREJIEp4ONM7fSx6qYXEZHjoYCPM6EBA0gdNkzj0ouIyHFRwMehSE42JctX4Coq/C5FREQSlAI+DkVycnAlJZSuWuV3KSIikqAU8HEoPGUKBALqphcRkWOmgI9DwcxMMk45RePSi4jIMVPAx6nItBzK3nuf6oICv0sREZEEpICPU5GcHKipofjNN/0uRUREElBMA97MLjSzdWa2wcy+38x6M7O7outXm9nERuuDZvaOmT0fyzrjUcb48Vg4rG56ERE5JjELeDMLAn8CZgAnArPM7MRGzWYAo6Kvm4G/NFr/DWBNrGqMZ5aaSuSMMyhZogFvRETk6MXyDH4ysME5t9E5VwHMBS5v1OZy4GHnWQZ0N7P+AGY2CLgYuDeGNca1yLQcKrZsoTI/3+9SREQkwcQy4AcC2+p9zosua2ubO4HvAjUxqi/u1Q5bq8flRETkaMUy4K2ZZY2nSGu2jZldAux2zq1s9UvMbjazFWa2Ys+ePcdSZ9xKPeEEUvr00bj0IiJy1GIZ8HnA4HqfBwHb29hmGnCZmW3G69o/x8webe5LnHP3OOcmOecm9e7du71qjwtmRiQ7m+Kly3A1nbYjQ0REjkEsA345MMrMhptZKpALPNuozbPA9dG76acCh5xzO5xzP3DODXLODYtu97Jz7rMxrDVuRablUH3gAGVrOuW9hiIicoxiFvDOuSrgq8ACvDvhH3POfWBmt5jZLdFmLwAbgQ3A34BbY1VPogpPnQpo+lgRETk65lzjy+KJa9KkSW7FihV+l9HuNl56GSm9ezHk/vv9LkVEROKIma10zk1qbp1GsksAkZwcSlaspKaszO9SREQkQSjgE0AkJxtXUUHp22/7XYqIiCQIBXwCCJ9xBoRCGrZWRETaTAGfAALhMOEJEzTgjYiItJkCPkFEcrIp/3ANVQcO+F2KiIgkAAV8gqgdtlaPy4mISFso4BNE+sknE8jMVDe9iIi0iQI+QVgwSGTqFIqXLCGZxi4QEZHYUMAnkEhODlXbd1C5ZYvfpYiISJxTwCcQTR8rIiJtpYBPIKEhQwgNGKAb7UREpFUK+ARiZkSm5VC87E1cVZXf5YiISBxTwCeYSE4ONYWFlL3/vt+liIhIHFPAJ5jw1KlgRrG66UVE5AgU8AkmJSuL9HHjKH5DN9qJiEjLFPAJKDIth5JVq6gpLva7FBERiVMK+AQUyc6GykpKVqzwuxQREYlTCvgElHH66VhamqaPFRGRFingE1AgLY3w6acr4EVEpEUK+AQVycmmfP0GKnfv9rsUERGJQwr4BKXpY0VE5EgU8AkqbexYgllZ6qYXEZFmKeATlAUCRLKnUrxkqaaPFRGRJhTwCSySk0PVnj1UbNjgdykiIhJnFPAJrPY6vLrpRUSkMQV8AgsNGEDq0KEUL9GNdiIi0pACPsFFpuVQvHw5rqLC71JERCSOKOATXCQnB1dSQumqVX6XIiIicUQBn+DCkydDIKDpY0VEpAEFfIILdu1KximnaPpYERFpQAGfBCLTcih97z2qCwr8LkVEROKEAj4JRLKzoaaGkrfe8rsUERGJEwr4JJBx6qlYOKzn4UVEpI4CPglYaiqRM87QdXgREamjgE8SkZxsKrZsoTI/3+9SREQkDijgk0TdsLV6XE5ERFDAJ43UkSNJ6d1b1+FFRARQwB9RIk3DamZEcrIpXroMV1PjdzkiIuIzBXwzKqsr+cpLX+Ge1ff4XcpRieTkUH3gAOVr1/pdioiI+EwB34xQMIRzjsfWPUZlTaXf5bRZODsb0PSxIiKigG9R7thcdpfu5uWtL/tdSpuF+vQhbdRITR8rIiIK+JZ8YuAnGBAZwNy1c/0u5ahEcnIoWbmSmvJyv0sREREfKeBbEAwEuWbMNazYtYL1B9b7XU6bRXJycOXllK5c6XcpIiLiIwX8EVw56kpSA6nMWzfP71LaLDxpEoRCeh5eRKSTU8AfQVZ6FhcOv5BnP36WwopCv8tpk0AkQvjUUzVsrYhIJ6eAb8WssbMorSrl2Y+f9buUNotMy6FszRqqDhzwuxQREfGJAr4VJ/c6mZN7nsy8dfMSZuCbSHY2OEfJsmV+lyIiIj5RwLdB7thcNh3axJs73/S7lDZJP/lkApmZeh5eRKQTU8C3wYXDL6R7WveEeWTOUlKITJ1C8RtLEqbXQURE2pcCvg3SgmlcOepKXtn2CjuKdvhdTpuEs7Op3L6dyi1b/C5FRER8oIBvo2vGXINzjsc/etzvUtqki6aPFRHp1BTwbTSwy0DOGnQWT65/korqCr/LaVVo6FBCAwboOryISCelgD8KuWNz2V+2n4VbFvpdSqvMjMi0HIqXvYmrqvK7HBER6WAK+KOQPSCboV2HJszNdpHsbGoKCyn74AO/SxERkQ6mgD8KAQswc8xMVu1ZxYf7PvS7nFaFs7PBTN30IiKdkAL+KF0+8nIyUjIS4iw+JSuL9HHjNGytiEgnpIA/Sl1Tu3LR8It4YdMLHCo/5Hc5rYrkZFOyahU1xcV+lyIiIh1IAX8MZo2dRXl1OU9veNrvUloVycmBykpKVqzwuxQREelARwx4M/tsvffTGq37amu/3MwuNLN1ZrbBzL7fzHozs7ui61eb2cTo8nQze8vMVpnZB2Z2R9t3KfbG9BjDxD4Tmbt2LjWuxu9yjijj9NOxtDRdhxcR6WRaO4P/j3rv/9ho3ReOtKGZBYE/ATOAE4FZZnZio2YzgFHR183AX6LLy4FznHOnAhOAC81saiu1dqjcsbnkFeWxOH+x36UcUSAtjfDpEyleogFvREQ6k9YC3lp439znxiYDG5xzG51zFcBc4PJGbS4HHnaeZUB3M+sf/VwUbROKvuJqUPXzhpxHz/SeCXGzXSQnh/L166ncvdvvUkREpIO0FvCuhffNfW5sILCt3ue86LI2tTGzoJm9C+wG/uWci6up3ELBEFeNvorF+YvZVrit9Q18FIkOW1uiYWtFRDqN1gJ+bPTa+Hv13td+HtPKts2d4Tc+KGixjXOu2jk3ARgETDazk5v9ErObzWyFma3Ys2dPKyW1r6tHX03AAjy27rEO/d6jlTZ2LMGsLHXTi4h0IimtrB93HL87Dxhc7/MgYPvRtnHOHTSzV4ELgfcbf4lz7h7gHoBJkyZ1aDd+30hfzhlyDk+tf4pbJ9xKRkpGR359m1kgQCR7KsVLvOljzVq7uiIiIonuiGfwzrkt9V9AETAR6BX9fCTLgVFmNtzMUoFc4NlGbZ4Fro/eTT8VOOSc22Fmvc2sO4CZZQDnAWuPeu86wKyxsyioKGD+pvl+l3JEkZwcqvbsoWLDBr9LERGRDtDaY3LP13aNm1l/vDPoLwCPmNk3j7Stc64K+CqwAFgDPOac+8DMbjGzW6LNXgA2AhuAvwG3Rpf3B14xs9V4Bwr/cs49fwz7F3OT+k5iZPeRzFk7B+fi6j7ABiLZ2YCmjxUR6Sxa66If7pyr7Ra/ES9orzezTOAN4M4jbeycewEvxOsvu7veewfc1sx2q4HTWq0+DpgZM8fM5Odv/pzVe1dzau9T/S6pWaGBA0kdOpTiN5bQ4/rr/S5HRERirLWb7CrrvT+XaFg75wqB+B7hpQNdesKlREKRuH9kLjIth+Lly3EV8T+fvYiIHJ/WAn6bmX3NzK7Au/Y+H+qui4diXVyiiIQiXHbCZSzYvIB9pfv8LqdF4exsXEkJpatX+12KiIjEWGsBfxNwEvB5YKZz7mB0+VTggdiVlXhyx+RSWVPJU+uf8ruUFkWmTIFAQMPWioh0Aq3dRb/bOXeLc+5y59zCestfcc79LvblJY4R3Ucwpd8UHvvoMapqqvwup1nBrl3JOOUUTR8rItIJtHYX/bNHenVUkYli1thZ7CzeyWt5r/ldSovCOdmUvvce1YWFfpciIiIx1FoXfTbe4DOvA78Dft/oJfWcNfgs+kX6xfXNdl1ycqCmhpI342rkXxERaWetBXw/4IfAycAfgPOBvc6515xz8Xua6pOUQApXj76aZTuWsfHQRr/LaVbGqadi4bCuw4uIJLnWrsFXO+fmO+duwLuxbgPwqpl9rUOqS0BXjrqSlEAK89bO87uUZllqKuEzJmlcehGRJNfaGTxmlmZmVwKP4g1KcxcQv7eK+6xXRi8+NfRTPPvxs5RUlvhdTrO65ORQsXkzldsbTw0gIiLJorWb7B4CluA9A3+Hc+4M59xPnXP5HVJdgpo1dhZFlUU8vzEuR9etmz5W3fQiIsmrtTP4zwGjgW8AS8ysIPoqNLOC2JeXmE7tfSrjeoyL2/HpU0eOJKV3b3XTi4gksdauwQecc5nRV9d6r0znXNeOKjLRmBm5Y3PZcHADK3at8LucJsyMSE42xUuX4mo04rCISDJq9Rq8HJsZw2fQNbVr3D4yF8nJofrAAcrXxuUsvCIicpwU8DGSkZLBp0d+mpe3vszukt1+l9NEeKqmjxURSWYK+BiaOWYm1a6aJz56wu9Smgj17UPaqJEatlZEJEkp4GNoSNchTBs4jcc/epzK6srWN+hgkZwcSlaupKa83O9SRESknSngY2zW2FnsLd3Lv7f+2+9SmghnZ+PKyyl9+22/SxERkXamgI+xaQOmMbDLQOasneN3KU1EzjgDQiEOzJ2Hq672uxwREWlHCvgYCwaC5I7J5e3db/PRgY/8LqeBQCRCr1u+TOGCBeR/81vqqhcRSSIK+A5wxagrSAumxeUjc71vu42+P/g+hf/6F9tu/jLVRUV+lyQiIu1AAd8BuqV1Y8bwGTy/8XkKKuJvAMAeN9zAgN/+hpKVK9ly/fVU7d3rd0kiInKcFPAdJHdsLqVVpTy74Vm/S2lWt0svZfCf/0TFxk1svu46KvLy/C5JRESOgwK+g5zU8yTG9xrPvHXzqHHxOTxsl+nTGfLA/VQfPMTmWbMoW7fO75JEROQYKeA7UO7YXDYXbGbZjmV+l9Ki8GmnMezRR7BAkC2f/RwlK+JvLH0REWmdAr4DXTDsAnqk94jLm+3qSxs1imFzZpPSqxdbb/oihS+/7HdJIiJylBTwHSg1mMqVo67ktbzX2F603e9yjig0YABDZ/+dtNGjyfva1zn45FN+lyQiIkdBAd/Brhl9DQCPrXvM50pal5KVxdAHHyAyZQo7fvQj9t17r98liYhIGyngO1j/Lv05a9BZPLX+Kcqr439gmUAkwuC7/0LXi2aw+3e/Z9dvfqs55EVEEoAC3gezxs7iQPkBFm5e6HcpbWKpqQz43e/Iuu469t9/Pzt++CNcZfxNniMiIocp4H0wtf9UhnUdFvc329VngQB9b/8Rvb7+NQ49/TR5X/0aNaWlfpclIiItUMD7wMzIHZvL6r2r+WDvB36X02ZmRu9bb6Xff/+YokWL2PqFm6g+dMjvskREpBkKeJ9cdsJlZKRkxOUsc63Jys1l4P/+L2Xvv8+Wz36Oyl27/C5JREQaUcD7JDM1k0tGXML8zfM5WHbQ73KOWtcLL2Dw3+6hMj+fLbOupXzTJr9LEhGRehTwPsodm0t5dTn/2PAPv0s5JpGpUxny8MPUlJWx5drrKH3vfb9LEhGRKAW8j0Znjeb0vqczb908qmuq/S7nmGScfBLDZv+dQDjM1htuoHjpUr9LEhERFPC+yx2bS35RPovzF/tdyjFLHTaMobNnExo4kG03f5mC+fP9LklEpNNTwPvs3CHn0jujN3PWJd7NdvWF+vZh6KOPkD5+PPnf+g8OzEns/RERSXQKeJ+FAiGuHn01b+S/wdaCrX6Xc1yC3box5L576fLJT7Lzjp+w5//+hHPO77JERDolBXwcuGr0VaRYCvPWzfO7lOMWSE9n0B/votsVV7D3//6PXT/9Ka46Me8vEBFJZAr4ONA73Jtzh57LPzb8g9KqxB8dzlJS6P+Ln9Pjpi9wYPYc8r/9bWoqKvwuS0SkU1HAx4ncMbkUVhTywsYX/C6lXZgZfb/zHfp859sUvjifvFtuobqo2O+yREQ6DQV8nDi97+mM7D6SuevmJtV165433UT/X/yC4jffYuuNN1K1f7/fJYmIdAoK+DhhZswaO4u1+9eyas8qv8tpV92vvIJBf/wj5R99xJZrr6MyP9/vkkREkp4CPo5cMuISuoS6JOT49K3JPOdshtx/H1X797N51rWUr1/vd0kiIklNAR9HwqEwl4+8nIVbFrK3dK/f5bS78OmnM/SRR8A5Nn/2c5S8/Y7fJYmIJC0FfJyZOWYmVTVVPPnRk36XEhPpY0YzdM4cUrp3Z+sXvkDRa6/5XZKISFJSwMeZ4d2GM7X/VB7/6HGqaqr8LicmUgcNZOjsv5M2YgTbbr2NQ88843dJIiJJRwEfh2aNncWukl28uu1Vv0uJmZSePRny8EOEzziD7d/7PvsefNDvkkREkooCPg6dNegs+kf6M3ftXL9Lialgly4MvuevZF5wAbt/9Wt2//5/kuoRQRERPyng41AwEOSaMdfw5s432Xhwo9/lxFQgNZWB//N7uufOZN/f/saO22/HVSXnpQkRkY6kgI9TV466klAglJSPzDVmwSD9fvxjet12G4eefIq8r3+DmrIyv8sSEUloCvg41SO9BxcOu5DnNj5HcWXyD/FqZvT+2lfpe/vtFL3yClu/+EWqCwr8LktEJGEp4ONY7thciiuLee7j5/wupcP0+Ox1DPjdbyldtZotn7ueyt27/S5JRCQhKeDj2Cm9TuHEnicyd21yjU/fmm4XX8zgv/yFim3b2HLdZ6nYutXvkkREEo4CPo6ZGbljcvn40Mes2LXC73I6VJczpzH0wQeoKSxk87XXUfbhh36XJCKSUBTwcW7G8Bl0S+vWKW62ayxj/HiGzv47lhpiy/U3UPzmW36XJCKSMBTwcS49JZ0rR17Jy1tfZmfxTr/L6XBpI0YwbPZsUvr1ZduXvkTBwoV+lyQikhAU8AngmjHXUONqeOKjJ/wuxRehfv0Y9uijpJ94Ivnf/BYHHnvM75JEROKeAj4BDMocxCcGfYInPnqCyupKv8vxRbB7d4bcfx+RM6ex879+zN67/9qpbjwUETlaMQ14M7vQzNaZ2QYz+34z683M7oquX21mE6PLB5vZK2a2xsw+MLNvxLLORJA7Jpd9Zft4aetLfpfim0A4zOA//Ymul13KnjvvZNcvf4mrqfG7LBGRuBSzgDezIPAnYAZwIjDLzE5s1GwGMCr6uhn4S3R5FfD/nHPjgKnAbc1s26lMGziNwZmDO+XNdvVZKMSAX/2KHjfcwIGHH2H7d7+Hq6jwuywRkbgTyzP4ycAG59xG51wFMBe4vFGby4GHnWcZ0N3M+jvndjjn3gZwzhUCa4CBMaw17gUswMwxM3ln9zus27/O73J8ZYEAfb7/PXr/v/+g4Pnn2XbbV6kpKfG7LBGRuBLLgB8IbKv3OY+mId1qGzMbBpwGvNncl5jZzWa2wsxW7Nmz53hrjmufHvlp0oPpnf4sHrwxAnp96Uv0/9lPKX7jDbbceCNVBw74XZaISNyIZcBbM8sa3xV1xDZm1gV4Evimc67Zgcmdc/c45yY55yb17t37mItNBN3SunHRiIt4YdMLHCo/5Hc5caH7VVcx6K4/UL5mLVs++zkqd+zwuyQRkbgQy4DPAwbX+zwI2N7WNmYWwgv3vzvnnophnQkld0wupVWlPLPhGb9LiRuZ553H4Hv/RtWuXWyedS3lH3/sd0kiIr6LZcAvB0aZ2XAzSwVygWcbtXkWuD56N/1U4JBzboeZGXAfsMY59z8xrDHhjOs5jgm9JzBv3TxqnO4grxWZPJmhjzyMq6piy7XXUbpqld8liYj4KmYB75yrAr4KLMC7Se4x59wHZnaLmd0SbfYCsBHYAPwNuDW6fBrwOeAcM3s3+rooVrUmmtyxuWwt3MrS7Uv9LiWupI8bx7DZfyfQrRtbPn8jRa8v9rskERHfWDINFjJp0iS3YkXyT8pSUV3B+U+cz/he4/njuX/0u5y4U7VnD1tv/jLlGzYw4Je/pNslF/tdkohITJjZSufcpObWaSS7BJQaTOUzoz7Da3mvkV+U73c5cSeld2+GPvwQ4QkT2P6d77D/kUf9LklEpMMp4BPUNWOuwcyYt26e36XEpWBmJoPv/Rtdzj2HXT//Obv/8AcNbSsinYoCPkH1i/TjnMHn8I/1/6C8utzvcuJSIC2NQXfeSferr2LfX+5m02c+w567/kjpqlW46mq/yxMRiSkFfALLHZvLwfKDzN803+9S4palpNDvJz+h73/9J4G0dPbefTebZ+ay/sxPkP/d73Lo+X9qgBwRSUq6yS6BOef49DOfJpwSZs4lGt2uLaoOHKD4jSUULXqN4tcXU33gAAQCZJx6Kl3Omk6X6dNJGzcO70lNEZH4dqSb7BTwCW72mtn88q1fMvui2ZzS+xS/y0korrqasg8+oOi1RRQtWkTZe+8B3k16kU98gi7TpxOZlkMwM9PnSkVEmqeAT2JFFUWc+/i5nDf0PH5+5s/9LiehVe3dS9HixRQvWkTR4jeoKSiAlBTCp51Gl7OmE5k+nbRRo3R2LyJxQwGf5H627Gf8Y/0/eOnql8hKz/K7nKTgqqooXbWq7uy+fO1aAFL696fL9Ole4E+ZQiAS8blSEenMFPBJbsOBDVzx7BV8c+I3uemUm/wuJylV7tpF8euvU/TaIoqXLKGmuBgLhQifcUbd2X3qsGE6uxeRDqWA7wS+sOAL5Bfm88KVLxAMBP0uJ6m5igpK3n6HokWLKFr0GhUbvMltQkOG0OUTn6DLWdMJT55MID3d50pFJNkp4DuBhZsX8v9e+3/88Zw/8snBn/S7nE6lIi+f4tcXeWf3y5bhysqwtDTCU6dEu/PPInXQIL/LFJEkpIDvBCprKrnwiQsZmTWSv57/V7/L6bRqysspeWs5Ra8vovi1RVRs2QJA6ogRddfuM04/nUBqqs+VikgyUMB3Enevups/vfsnnvv0cwzrNszvcgSo2LyZokWvU7RoESVvvYWrqCAQDhPOyabLJ6bTZfonCPXv73eZIpKgFPCdxN7SvZz/xPnkjsnle5O/53c50khNSQnFb75J0SLv7L5y+3YA0kaPrhtkJ2PCBCwU8rlSEUkUCvhO5LuvfZfF+Yt56eqXCIfCfpcjLXDOUfHxx95jeK+/TsmKFVBVRSAzk8i0aV53/ifOJKV3b79LFZE4dqSAT+noYiS2csfm8uLmF3lh0wtcNfoqv8uRFpgZaSNHkjZyJD1v+gLVRUUUL13qDbLz2iIK53vzC6SfdFLd2X36KadgQT0hISJtozP4JOOc4+rnrsbheOLSJ/RcdgJyzlG+bl3dIDul77wDNTUEu3ePDqH7CSJnnklKlgY1EunsdAbfiZgZuWNzuWPpHbyz+x0m9p3od0lylMyM9LFjSR87ll5fvpnqQ4cofuMNL/AXL6bguefAjIzx44mcNZ0u088i/cRxWECTQ4rIYTqDT0IllSWc9/h5nDnwTH5z1m/8LkfakaupoeyDDyla9Jo3Qc7q98A5gr161Q2yEznzTIJduvhdqoh0AJ3BdzLhUJjLR17O3HVz2Vu6l14ZvfwuSdqJBQJknHIyGaecTO/bbqNq/36KFy/2rtu//DKH/vEPLDWVyJln0vXCC+hyzjkKe5FOSmfwLSk9AKmZEEzMY6AtBVu45B+XcNuE27jl1Fv8Lkc6gKuqonT1agoXLKRgwQKqdu7EQqGGYa+pb0WSih6TO1oVJXDPJ2HgRLj8z5Cg1zZv+dctrD+wnvlXzScU0LPVnYmrqaF01SoK5y+gYOFCqnbsUNiLJKEjBXxiJlespYbhlKtg1RyY/z1I0IOg3LG57C7dzStbX/G7FOlgFggQPu00+v7g+4z890sMmzuHrOuuo2ztWrZ/7/usz5nGtlu+wsGnn6a6oMDvckUkBnQG3xLnYOHtsPT/YPp34Jzb2+f3dqDqmmoueuoiBmYO5P4L7ve7HIkDzjnKVq+m4MX5FCxcQNX2HRAK0SUnh8wLLyTz3HMIdu3qd5ki0ka6ye5YmMGnfgblBbDot5DWFaZ93e+qjkowEOSaMddw59t3suHABkZmjfS7JPGZmZFx6qlknHoqfb73XS/s5y+gYMF8il57jR0Ke5GkoTP41tRUw5M3wQf/gEvuhEk3tu/vj7EDZQc47/HzuGLUFdw+NfF6IaRjOOcoe+89Cl6cT+GCBd44+aEQkZxsul4QDftu3fwuU0Qa0U12x6uqAuZeCxtegs/c612fTyA/WvwjXtryEv+++t90SdUjU3JkdWE/fwGF8+cfDvvsqV7Yn3euwl4kTijg20NFCfz9Ktj2JuTOhtEXxOZ7YuD9ve8z65+z+MHkH3DtuGv9LkcSiHOOsvffp2D+fArnL6AyPx9SUohkZ9P1wgvIPPdcgt27+12mSKelgG8vZQXw0KWwZy1c9wQM/0TsvqudzXp+FiVVJTx9+dMan16OiRf2H1Aw/0WFvUicUMC3p+J98OBFcCgPrn8WBp0e2+9rJ89seIbb37idez91L1P6T/G7HElwtWFfuGA+BfMXUJmX54X91Knec/bnnqvJcEQ6gAK+vRXsgPsv8O6w//wL0PfE2H/ncSqvLue8x8/j9L6nc+fZd/pdjiQR5xxlH3xI4fwXG4b9lClkXngBmeedp7AXiREFfCzs3wQPzABXA1+YDz1GdMz3Hof/Xfm/PPjBgyz4zAL6Rfr5XY4kobqwrz2z37YNgkEiU6cq7EViQAEfK7vXeiGf1gVunA/dBnbcdx+D/KJ8Zjw5gy+e8kW+PjGxnumXxOOco+zDD73hchcsoHLrVi/sa8/szz9fYS9ynBTwsZT/Njx0GXTtDze+CJH4nrnta//+Gqv3ruZfV/2L1GCq3+VIJ+Gco3zNGm8EvQZhP5nMCy4k8/zzSOnRw+8yRRKOAj7WNr8Bj14JvcfADc9Bevw+I/xG/hvc8tIt/OoTv+LiERf7XY50QnVhHx1Br3KLwl7kWCngO8JHC2HuLBh0Bnz2KW/CmjhU42q47OnL6J7WnUcvetTvcqSTc85RvnZt9Mz+cNiHJ5/hDarzqfMV9iJHoIDvKO8/CU/cBCPPhdw5kBKfXeCPfPgIv1n+Gx675DHG9RzndzkiQL2wj46gV7FlCwQChCdP9p6zP/98Unr29LtMkbiigO9IKx+C574OJ34arrofAkF/62lGQUUB5z1+HjOGz+COnDv8LkekCecc5evW1Y2gV7F5s8JepBkK+I625P9g4Y/gtM/CpX+EQMDvipr47yX/zT83/pOXrn6Jbmnxe8+AiHOO8o8+ouDFFxuEfeqwYaSOGE7a8BGkjhhB2ojhpI4YQTAz0++SRTqMpovtaDlfhbJDsOg33jSzF/zCm342jswaO4sn1z/J0xue5oaTbvC7HJEWmRnpY8aQPmYMvb/xDco/+ojCf71E+bq1lG/cRNGrr0FVVV37lN69SR0xokn4p/Trh8XhwbZIrCjgY+XsH3oj3S37s3dX/Se/73dFDYzpMYaJfSYyb908Pnfi5wiY/uGT+Fc/7Gu5ykoq8vKo2LSJio0bKf94IxUbN1Lw/D+pKSw8vG1GBqnDh5E24gQv/EeMIHX4CFKHDSWQlubH7ojElAI+Vszggl9CeSG8+kvvTD77Vr+raiB3bC7fXfRdXt76MucNPc/vckSOiYVCpA0fTtrw4XDOOXXLnXNU79tH+caNVGzcRMWmjZRv3ETpO+9Q8Pzz9X6BERo0qN4Z/3DSTjiB1BEjNBCPJDQFfCwFAnDpXV7IL/gBpGXCxM/5XVWd84acx8AuA/nWq99iSv8pXDv2Ws4adBbBOLwxUORomRkpvXqR0qsXkcmTG6yrKS2lYssWyj/+uEH4lyx7E1deXtcu2L17o+5+78w/NGgQFtT/JxLfdJNdR6gqhzm5sPFV7876k67wu6I6B8oO8OT6J5m3bh47i3cyIDKAa8Zcw2dGfYbu6d39Lk+kQ7maGiq376Bik9fNX75xExUff0z5pk1U79tX185CoehNfiMadPenDR9GIBLxcQ+ks9Fd9PGgohgeuRLyV8KsuTAqvrrEq2qqeHXbq8xZO4e3dr5FWjCNGcNnMGvsLE7sGf+z5YnEWvXBg5Rv2tTgjL9i40Yqtm2D6uq6din9+kXv6G94rT+lT28szm62lcSngI8XZYfgwUtg73r43FMwNMfvipq1/sB65q6dy3Mbn6O0qpQJvScwa+wszh96PqFgyO/yROKKq6igYuvWJtf6KzZupKa4uK5dIBI5/Dhf/Wv9gwdjqfE5KJbEPwV8PCna481AV7gTPv8cDDjN74paVFBRwDMbnmHu2rlsLdxKr4xeXD36aq4afRV9wn38Lk8krjnnqNq9Jxr40fDf6L2v2rnzcMNgkNTBg5uG/4gRBLtpjAo5MgV8vDmUD/dfCJXF3gx0vce0vo2PalwNS7YvYfaa2SzOX0zQgpw/9HxmjZvFhN4T1O0ocpSqi4qp2Ly5SfhXbN6Mq6ysaxfs2ZOUvn1IyepBMCuLYFYWKT2you97EMzqTkqP6Lru3XXjXyekgI9H+z72Qj4QhC/Mh6xhflfUJlsLtjJ33VyeXv80hZWFjO0xlmvHXsuM4TNIT0n3uzyRhOaqq6nMz/dC/+ONVGzeRNWevVQd2E/1gYNU799PTVFR8xubEezalWBt4PfIIiUri2D3rOiyegcDWT1IyeqOhcM6QE9wCvh4tesDeOAiyOgON8735pRPECWVJfxz0z+ZvWY2Gw5uoFtaN64ceSUzx85kYJeBfpcnkrRcRQVVBw5SffAA1fv3U33gAFUHDlC9/0D0/eGDgdp19Uf6q8/S0g4fDNQdCGQdPhjoXu9AoUcPgt26YSl6ujqeKODjWd4KeOgy6D4EbnwBwok1NaZzjhW7VjBn7Rxe3voyNa6GswafxbVjr2Vq/6k6OxDxmXOOmqKiw4EfPRCoPrC/xQODVnsJshodDGQ16jWIrkvJylIvQYwp4OPdpkXw6FXQ9yS44VlvQJwEtLN4J4+te4wn1z/J/rL9DO82nNwxuVw+8nIiIT0bLJIoXEUFVQcPRg8EvJ6CY+4lSE2tF/jdo/cO1DsYyOpBsEcWqYMHk9K3r+YLOEoK+ESw7kWYex0MyYbPPgGhDL8rOmYV1RUs2LyA2Wtm8/6+94mEIlx2wmXkjs1lRLcRfpcnIu2s1V6CRgcD1QcONJgnoFbdfAH1niRIHTGC1KFDCaTrHp/mKOATxerH4akvwahPQe7fIQmeOX9vz3vMWTuH+ZvnU1lTSXb/bGaNncX0QdM1JK5IJ9agl2DfPiq2bqVi06a6MQQqt2+H2nwyIzRwIKnDhzd9lLBnz059CUABn0hW3A/PfwtO/gxc+TfvLvsksK90X92QuLtLdjOwy0BmjpnJFSOv0JC4ItJE7XwBdUMGb9pE+aaNVGzajCstrWsX6NqVtOHDSR0+/PBYAiNGeAMIhRL/JKk1CvhEs/hOeOnHMPEGuPQPcTeX/PGoqqnilW2vMHvNbFbsWkFaMI2LR1zMrLGzGNtjrN/liUicczU1VO3cWXemX76pdgTBTVTt3n24YUrK4QGEhg9L2gGEFPCJ6KU7YPH/QM7X4PyfJlXI1/rowEfMWTuHf278J6VVpZzW5zSuHXst5w49l1Ag+Y+8RaR9VRcVUbFp0+Gz/o3eWALlm7dAowGE0qJn/IfnCxhOaMCAhBssSAGfiJyDF74Dy/8G59wO07/jd0Uxc6j8kDck7rq5bCvcRu+M3lw9+mquHnM1vTJ6+V2eiCQ4V1V1eAChjZsanPVXHzhQ185SU+tmCWxwrX9Y/M4S6FvAm9mFwB+AIHCvc+5XjdZbdP1FQAnweefc29F19wOXALudcye35fuSKuABamrg6Vtg9TyY8RuY8mW/K4qpGlfD4vzFzFk7h8X5i0kJpHD+0PO5duy1nNr71E59I42IxEbVgQNNzvrLN22kclue929wVN0sgY3u8E/p08fXf5t8CXgzCwIfAecDecByYJZz7sN6bS4CvoYX8FOAPzjnpkTXTQeKgIc7bcADVFfBY9fDun/Cp/8CE671u6IOsaVgC3PXzuXpDU9TVFnEuB7jmDV2lobEFZEOUVNRQWVbZgkMhxve4Bc9AEgdOpRAWlrM6/Qr4LOB/3bOXRD9/AMA59wv67X5K/Cqc25O9PM64JPOuR3Rz8OA5zt1wANUlsHsa2Dz63D1Q3DiZX5X1GFKKkt4fuPzzFk7hw0HN9A9rTtXjrqSmWNmMqDLAL/LE5FOxjlH1Z49TUK/fNNGqrbvONwwECA0aFCT5/rTRo1q15v8jhTwsRxUeCCwrd7nPLyz9NbaDAR2IIeF0iF3NjzyaXjyJkibByec43dVHSIcCnPNmGu4evTVLN+5nDlr5/DgBw/y4AcP8slBn2TWuFlM6TdF3fci0iHMjFCfPoT69CEytWGk1ZSUULFlS5Oz/pJlb+LKywHo893v0vMLN3ZIrbEM+Ob+xW3cXdCWNkf+ErObgZsBhgwZcjSbJpa0LnDd4/DgJd6Id5/7BwyZ6ndVHcbMmNx/MpP7T2ZH0Q4e++gxnvzoSV7e9jIjuo1g1thZXHrCpRoSV0R8EwiHSR83jvRx4xosdzU1VG7fQcWmTaQOG9ph9aiLPtEU7fammS3eC59/HvqP97si35RXlzN/03xmr53Nh/s+pEuoC5ePvJyZY2YyvNtwv8sTEYm5I3XRx3JU/+XAKDMbbmapQC7wbKM2zwLXm2cqcKg23KUFXfrA9c94E9I8cgXsXe93Rb5JC6Zx+cjLmXvxXB696FHOGnwW89bN47KnL+PL//oyr217jeqaar/LFBHxRawfk7sIuBPvMbn7nXM/N7NbAJxzd0cfk/s/4EK8x+RudM6tiG47B/gk0AvYBfzYOXffkb6vU5zB19q7Ae6/AFLS4Avzvelmhb2le3nyoyd5bN1j7C71hsTNHZPLFaOuoFta8oxeJSICGugmee1Y7V2Tj/SELyzwzu4FgMqaSl7e+jJz1s5h5a6VpAfTGd97PIMzBzM4czBDug5hSOYQBmcOJhwK+12uiMgxUcAns61venfX9xjhXZPPyPK7orizbv86Hv/ocdbuX8u2wm3sL9vfYH3P9J4M6eqFfW3o137WWb+IxDMFfLL7+GWYPRP6nwqfe9q7415aVFRRxLbCbWwt3Mq2wm3e+4KtbC3cyu6S3Q3adkvrxuAugxnc1Qv/2uAfnDmYnumde5pKEfGfAr4zWPMcPHYDDJsG1z7uPTsvR620qpT8wvy68N9asLXuYGBH8Q5q3OGhK8Mp4QaBX/8AoE+4DwGL5T2sIiIK+M7j3Tne2PVjLoZrHoKgZmRrT5XVleQX5Tc4+689AMgryqOqpqqubWog1Qv+rvXCP3MIg7sOpn+kPymBWA5BISKdhV8j2UlHmzALygvhxe/AM7fBp++GgM4i20soGGJYt2EM6zasybrqmmp2luysC/z63f7Lti+jrLqsrm2KpTCgy4C6bv+6a/9dBzOoyyBSg6kduFcikqwU8Mlmys1Qfghe/hmkdoGLf5+Uc8nHm2AgyMAuAxnYZSDZZDdY55xjT+mehuFfuJWtBVtZtXsVRZVFdW0No1+kX13g1575D8ocpDv+ReSoKOCT0Se+DWUFsOQuSO8G5/3Y74o6NTOjT7gPfcJ9mNSvYU+ac46D5QfrAr/+AcC/t/ybA+UHGrTvndG76aN+0csAXVO7duRuiUicU8AnIzM4/yded/3i/4H0rnDmt/yuSpphZmSlZ5GVnsWpvU9tsr6wovDwNf+Cw9f+l25fyjMfP9Ogbfe07gzJHEL/Lv3pE+5D33Bf+kb60i/cj77hvvQK9yIU0H0ZIp2FAj5ZmXnd8+WF8NJ/e0PbnvFFv6uSo5SZmsmJPU/kxJ4nNllXUllCXlGed9YfDf+thVtZt38di/IWUVpV2qC9YfTM6OkFfzT8aw8E+kW8g4A+4T6kp+gJDJFkoIBPZoEgXHE3VBTBP78NqZlw6ky/q5J2Eg6FGZ01mtFZo5usc85RUFHA7pLd7CrZxa7iXd7P6Gtr4VaW71pOYUVhk227pXVrchBQ2wvQN+It75KqsRZE4p0CPtkFQ3D1g/D3q+Hpr3hn8mMv8rsqiTEzo1taN7qldWNU1qgW25VUlhw+CGh8IFC8iw/2fdBk5D/wxgCoDfv6wV//ffe07hoISMRHeg6+sygvhIcug10fwHWPwYhP+l2RJIiK6gp2l+xusTdgV/Eu9pTuaTAIEHhjAfQJ9zkc/vUPAqKfe6b3JBgI+rRnIolPA92Ip2Q/PHgxHNjiTTk7+Ay/K5IkUVVTxb7SfU16A3aW7PSWRQ8KKmsqG2wXtCC9Mno1CP9+kX4NbhLsk9GHkAZtEmmWAl4OK9wJ918Ipfvh8y9Av5P9rkg6CeccB8oPNAj8ncU72VWyq+7AYGfxziY3BwL0SO/RoBegX6QfAyID6h4Z7JbWTZcDpFNSwEtDB7Z4IV9T5c0l3/MEvysSAbyDgKLKInYV1wv9kp0NPu8q2cWh8kMNtssMZdYNBlT/NaTrEM0LIElNAS9N7VkHD8yAUBiufsibiS6oey4lMdROClQ7MNC2wm1sK/IeF9xetJ0q13BegIGZA5uE/6DMQRoaWBKeAl6at/1deOhSKC+AlAwv5AdNgoGne6/uQzTMrSScqpoqdhbvrAv+vMK8BgcCJVUldW0No2+kb5Pwr31lpmb6uCcirVPAS8sKd8Hm1yF/JeStgB2roLrcWxfpHQ37STDodBgwETK6+1quyPFwzrG/bH/DM/96r8aPBHZP6153tt84/Htn9NZ1f/GdAl7arroSdr0fDfyVkL8C9n50eH3PUV7o157p9z0ZUtTFKcmhuLK4yRn/1sKt5BXmsaN4R4NHAdOD6c0G/+DMwfTv0l/DAkuHUMDL8Sk7BPlve6Ffe6ZfvNtbF0yD/uMbnulnDVfXviSdyupKthdvb3LWn1eYR15hXoMpgYMWpF+kX4td/5oVUNqLAl7al3NwKM87u68909/xLlRGr21m9Dh8Hb/2TD/cw9eSRWKpxtWwt3RvgxkB63oCirY1ueu/Z3rPJjf81b7vkd5DXf/SZgp4ib3qKtizxju7rz3T370GiP73lTU8GvbRwO93CoQ0qYl0DgUVBU2Cv/ZgYHfJbhyH/x0Op4QbhH/39O5kpGQQTgkTDoUJp4S9z9H34ZD3OSMlQ48DdkIKePFHeaF3p379M/3C7d66QMgbZGfgpMNn+T1OgID+gZLOpby6vOkjf9FXflF+k9H/jqQ26BscALRwQND4fUao6bJwKEx6MF09CnFMAS/xo2D74ev4+Sth+zvebHcA6d28O/Xrn+l36e1vvSI+cs5RVl1GSWUJpVWllFSVUFJZQklV9HNlo5/NrC+pin6uLK1b3txogS0x7OgOEJpp09xBRigQ0oFDOzhSwGtkE+lYXQd4r3GXep9rqr1Bd/JXHj7Tf/1/wFV767sPqXcD3yToNx5SdYOSdA5mVndG3p5qXA1lVWUtHxC0csBQWllKYUUhu0p2NVhfXvuIbRsELdjgwCA9JZ3UYCppwTTvZyCt7n2D5fV+HtOyQGqnObDQGbzEn4pi73n8ujP9t+HQVm+dBaHvSfVu4JsEvUara18kDlTVVLWpR6G5NuXV5VRUVzT42dyyipqK464zNdD2A4bjOpjogAMMddFL4ivcBdvfjgb+Csh/B2rvTE7NhIGnNTzTz+znb70iEhM1robKmspWDwSOZ1lL6ypqKo6ql6I535n0Ha4/6fp2+tNQF70kg8y+MGaG9wKoqYF9Gw537eetgCV/9CbQAeg6sOFjev0nQFoX38oXkfYRsEDdWbEfnHN1BxjHchAxse/EDqtVAS+JKRCA3qO914RZ3rLKMti5ul7X/gpY86y3zgJe6Gf2h679672P3hOQ2d976dE9ETkCM6u7LyCT+J6rQAEvySOUDoMne69axfsOP5d/cIt3F/+uD2HDvw/fvV9fuCdkDogeBAxo/n16d43UJyJxTwEvyS3SE0Z/yns1VlbgBX7hdijY0eh9vvcIX/GeptuFwk3P/uveRw8CuvSFQDD2+yci0gIFvHRe6V29V5+xLbepqoDCHd6rYHv0IKDe+61LvQOCxoORWNAL+SY9AfUuDWT21yN/IhIzCniRI0lJhayh3qslNTVQsi969t/MQcDe9bDxNSgvaLptevcj9wR0HQgZWbokICJHTQEvcrwCAW/EvS69of+pLbcrL4oGf7531l93QBB9v/M9KNoNNHp0NSXde+zvSDcJdukLQU1PKiKHKeBFOkpaF0gbBb1GtdymuhKKdjXqCcg/fI9A/gpYswOaPItr0KWPF/rhnpCW6b3Sux1+n5YJaV0P/0zvenh5KKLBgkSSjAJeJJ4EQ9BtkPdqiXNQsr/hDYH1LwmUHYRD27zJfsoKoLK4DV9s9cI/s2H41x0YdG1mXdeG26VGdDlBJE4o4EUSjZn3dECkpzftbmtqqr2wLy/07gOo/76soJl10Z8l++HAlsPLKkvaUFugaU9Bk4OGltbV71EI60BB5Dgp4EWSXSAIGd291/GoroKKwsM9Aw0ODOp9Lmt8oLAXDmw6vL5NBwrB5i8lND4wSI14BwOpXbz3qZFG76OvYKoOGKTTUcCLSNsEU7w7+jOyju/3VFe20KNQCGWHWu5tKNoN+z4+/LmqrO3fGUhpGv6hRgcBdeuOcMDQeBuNdSBxTAEvIh0rGIJwD+91PKqrvPsLKoqhosQbmbCi9nP0fWULy2u3KdpV73N0fe1UxW2RktG0t6C5XoTQEdY1Xp6Srt4GaRcKeBFJTMEUCHbznhRoL85BdUUzBwNFRz6QqCiud7BRDMV7G27TphsdoyzghX0ofPgSRCgj+gp7QzLXX5ZSf11Go1fYO2CoW1dve/U+JD0FvIhILTNISfNex9vDUF9NTbQ34QgHBc0eSBR5kyhVlkBlKZQe8C5NVJYeXtaWexqaE0w9hgOERuua3b7RAYXGZ/CNAl5EJNYCgeg4CF2Avu37u52DqvJ6gV8KVaWNDgIafa4qa3iAUP8gorLUm4OhyfYl4GqOvj4LtnKAED0gSEn1DjqCqd5BQTCt3vvUI79PaWV93fs0r+eik1wCUcCLiCQys+hZc4ynOnbOu0Gy7iCh8YFDo4OGqrJmDjAaHTSUHqjXvsK7PFJdGf1ZcXT3Q7SZtXIQEPJ6cI60/lgORGoPQnqcAN0GxmC/mlLAi4hI68y8kEpJPf5HLtuqprpe4NcL/gavypbfV5U3Wn6EttXNtK0s9Z7sqN/2eA9ELvwVTP1K7P7M6lHAi4hIfAoEvVeseyeO19EciGQN67CyFPAiIiLHI04PRDS7hIiISBJSwIuIiCQhBbyIiEgSUsCLiIgkIQW8iIhIElLAi4iIJCEFvIiISBJSwIuIiCQhBbyIiEgSUsCLiIgkIQW8iIhIEoppwJvZhWa2zsw2mNn3m1lvZnZXdP1qM5vY1m1FRESkZTELeDMLAn8CZgAnArPM7MRGzWYAo6Kvm4G/HMW2IiIi0oJYnsFPBjY45zY65yqAucDljdpcDjzsPMuA7mbWv43bioiISAtiOV3sQGBbvc95wJQ2tBnYxm1j6o7nPuDD7QUd+ZUiIpLkThzQlR9felKHfFcsz+CtmWWujW3asq33C8xuNrMVZrZiz549R1miiIhIcorlGXweMLje50HA9ja2SW3DtgA45+4B7gGYNGlSswcBx6KjjrBERERiIZZn8MuBUWY23MxSgVzg2UZtngWuj95NPxU45Jzb0cZtRUREpAUxO4N3zlWZ2VeBBUAQuN8594GZ3RJdfzfwAnARsAEoAW480raxqlVERCTZmHPt1qvtu0mTJrkVK1b4XYaIiEiHMLOVzrlJza3TSHYiIiJJSAEvIiKShBTwIiIiSUgBLyIikoQU8CIiIklIAS8iIpKEFPAiIiJJSAEvIiKShBTwIiIiSUgBLyIikoQU8CIiIklIAS8iIpKEFPAiIiJJSAEvIiKShBTwIiIiSSip5oM3sz3Alnb8lb2Ave34++KV9jO5aD+Ti/YzubT3fg51zvVubkVSBXx7M7MVzrlJftcRa9rP5KL9TC7az+TSkfupLnoREZEkpIAXERFJQgr4I7vH7wI6iPYzuWg/k4v2M7l02H7qGryIiEgS0hm8iIhIElLAA2Y22MxeMbM1ZvaBmX0juryHmf3LzNZHf2b5XevxMLN0M3vLzFZF9/OO6PKk2s9aZhY0s3fM7Pno56TbTzPbbGbvmdm7ZrYiuiwZ97O7mT1hZmuj/59mJ9t+mtmY6N9j7avAzL6ZbPsJYGbfiv4b9L6ZzYn+25SM+/mN6D5+YGbfjC7rsP1UwHuqgP/nnBsHTAVuM7MTge8D/3bOjQL+Hf2cyMqBc5xzpwITgAvNbCrJt5+1vgGsqfc5WffzbOfchHqP3iTjfv4BmO+cGwucivf3mlT76ZxbF/17nACcDpQA/yDJ9tPMBgJfByY5504GgkAuybefJwNfAibj/Td7iZmNoiP30zmnV6MX8AxwPrAO6B9d1h9Y53dt7biPYeBtYEoy7icwKPo/zznA89Flybifm4FejZYl1X4CXYFNRO8ZStb9bLRvnwLeSMb9BAYC24AeQArwfHR/k20/rwburff5P4HvduR+6gy+ETMbBpwGvAn0dc7tAIj+7ONjae0i2m39LrAb+JdzLin3E7gT73+mmnrLknE/HbDQzFaa2c3RZcm2nyOAPcAD0Usu95pZhOTbz/pygTnR90m1n865fOB3wFZgB3DIObeQJNtP4H1gupn1NLMwcBEwmA7cTwV8PWbWBXgS+KZzrsDvemLBOVftvC7AQcDkaDdSUjGzS4DdzrmVftfSAaY55yYCM/AuLU33u6AYSAEmAn9xzp0GFJPg3bdHYmapwGXA437XEgvRa86XA8OBAUDEzD7rb1Xtzzm3Bvg18C9gPrAK73Jwh1HAR5lZCC/c/+6ceyq6eJeZ9Y+u74931psUnHMHgVeBC0m+/ZwGXGZmm4G5wDlm9ijJt58457ZHf+7Gu147meTbzzwgL9rbBPAEXuAn237WmgG87ZzbFf2cbPt5HrDJObfHOVcJPAXkkHz7iXPuPufcROfcdGA/sJ4O3E8FPGBmBtwHrHHO/U+9Vc8CN0Tf34B3bT5hmVlvM+sefZ+B9z/aWpJsP51zP3DODXLODcPr6nzZOfdZkmw/zSxiZpm17/GuY75Pku2nc24nsM3MxkQXnQt8SJLtZz2zONw9D8m3n1uBqWYWjv7bey7eTZPJtp+YWZ/ozyHAlXh/rx22nxroBjCzM4HXgfc4fM32h3jX4R8DhuD9R3m1c26/L0W2AzMbDzyEd9dqAHjMOfcTM+tJEu1nfWb2SeDbzrlLkm0/zWwE3lk7eN3Ys51zP0+2/QQwswnAvUAqsBG4keh/wyTXfobxbkAb4Zw7FF2WjH+fdwAz8bqs3wG+CHQh+fbzdaAnUAn8h3Pu3x3596mAFxERSULqohcREUlCCngREZEkpIAXERFJQgp4ERGRJKSAFxERSUIKeJHjZGbOzH5f7/O3zey/2+l3P2hmV7XH72rle66OztL2Sqy/KxGZ2aejE1CJJAwFvMjxKweuNLNefhdSn5kFj6L5TcCtzrmzY1VPezGzFB++9tOAAl4SigJe5PhVAfcA32q8ovEZuJkVRX9+0sxeM7PHzOwjM/uVmV1nZm+ZN7/7CfV+zXlm9nq03SXR7YNm9lszW25mq83sy/V+7ytmNhtv4KbG9cyK/v73zezX0WX/BZwJ3G1mv23Uvk11RkdJfDJaz3IzmxZdPtnMlkQniVlSOxqdmZ0U/R3vRusfZWbDzOz9et9d1xNiZq+a2S/M7DXgG2Z2erSulWa2oN7Qn6+a2f+a2aJoj8QZZvaUeXNv/6ze7/5sve//a+3BkJkVmdnPzWyVmS0zs75mloM3Nvxvo+1PMLOvm9mH0drntuU/EpEO5/eUenrplegvoAhvStPNQDfg28B/R9c9CFxVv2305yeBg3jTRaYB+cAd0XXfAO6st/18vIPxUXjjsqcDNwO3R9ukASvwJu/4JN5kLMObqXMA3shZvfFGvnsZ+HR03at483M33qatdc4Gzoy+H4I37DPRP5eU6PvzgCej7/8IXBd9nwpkAMOA9+t9d/0/x1eBP0ffh4AlQO/o55nA/fXa/bpefdvr1Z6HN6rYOOA5IBRt92fg+uh7B1waff+ben/Gjf8etwNp0ffd/f5vUC+9mnv50dUlknSccwVm9jDwdaC0jZstd9FpI83sY2BhdPl7QP2u8secczXAejPbCIzFG3d+fL3egW54BwAVwFvOuU3NfN8ZwKvOuT3R7/w7MB14uh3qPA840cxqt+lq3jj53YCHzGwUXniGouuXAj8ys0HAU8659fW2bcm86M8xwMnAv6LbBPGmHa31bL36PqhX+0a86TrPBE4Hlke3z+DwhB8VePOTA6wEzm+hltXA383saVr/8xPxhQJepP3cCbwNPFBvWRXRS2HmpUlqvXXl9d7X1PtcQ8P/NxuPJ+0AA77mnFtQf4V5Y+8Xt1BfqwnagrbUGQCynXMNDm7M7I/AK865K8xsGN4ZNs652Wb2JnAxsMDMvgh8RMPLhumN6qjdL8ML7uxW6q1fa/16DXjIOfeDZratdM7V/nlX0/K/kRfjHRxdBvynmZ3knOvQqUBFWqNr8CLtxHkTRjyGd8Narc14Z4vgzYEd4uhdbWaB6PXuEcA6YAHwFfOmOcbMRps3o9yRvAmcZWa9otecZwGvHUM9zVkIfLX2g3mTw4B3Bp8fff/5eutHABudc3fhnXGPB3YBfcysp5mlAZe08F3rgN5mlh39XSEzO+koav03cJUdnumrh5kNbWWbQqB25r4AMNg59wrwXaA73kQpInFFAS/Svn4P1L+b/m94ofoWMIWWz66PZB1eEL8I3OKcK8ObWe1D4O3ojWl/pZUeuWhX9Q+AV4BVeHOOt9dUlV8HJkVvOvsQuCW6/DfAL83sDbyu9FozgffN7F28Sw4PO29u8J/gHYg8jzeVcXP7UQFcBfzazFYB7+LNJ94mzrkPgduBhWa2GvgX3nX6I5kLfMfM3sG7FPKomb2HNxPa/zrnDrb1+0U6imaTExERSUI6gxcREUlCCngREZEkpIAXERFJQgp4ERGRJKSAFxERSUIKeBERkSSkgBcREUlCCngREZEk9P8Bn+tO82uoBLkAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Plot the results\n",
"plt.figure(figsize=(8,8))\n",
"plt.plot(Ns, error[:,:])\n",
"plt.legend(legend)\n",
"plt.xlabel('Number of measurements')\n",
"plt.ylabel('MSE')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These plots show that the error in the determiantion of the Stokes vectors is proportional to the initial error amplitude, and can be decreased increasing the number of measurements."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.12"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
![](\"Optical_elements_polarimeter.png\")
\n",
"\n",
"Most optical element polarimeters, like this one, are usually divided into two parts: the polarization states generator (PSG) and the polarization states analyzer (PSA). The polarization states generator is the part of the polarimeter capable of generating many polarization states. In our case, it is composed by the first rotating polarizer and waveplate. If both of them are perfect elements, this PSG can generate all the possible polarization states. If the incident light wave is circularly polarized, all the generated states will have the same intensity. The second part is the polarization states analyzer (PSA). This part is the same as the light polarimeter described in the previous example. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.optimize import least_squares\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from py_pol import degrees\n",
"from py_pol.jones_vector import Jones_vector\n",
"from py_pol.jones_matrix import Jones_matrix, create_Jones_matrices\n",
"from py_pol.stokes import Stokes, create_Stokes\n",
"from py_pol.mueller import Mueller, create_Mueller"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Jones formalism\n",
"\n",
"The most general Jones matrix can be described using 4 parameters. One selection of parameters is:\n",
"\n",
"$J=\\exp\\left(i\\Phi\\right)\\left[\\begin{array}{cc}\n",
"J_{0} & J_{1}\\exp\\left(i\\delta_{1}\\right)\\\\\n",
"J_{2}\\exp\\left(i\\delta_{2}\\right) & J_{3}\\exp\\left(i\\delta_{3}\\right)\n",
"\\end{array}\\right]$,\n",
"\n",
"where the modules $J_i$ are real and positive and the phases $\\delta_i$ and the global phase $\\Phi$ are real and range between 0 and 360º. In principle, it is not possible to measure the global phase using an intensity detector without using interferometric techniques, so first we will concentrate in measuring the other parameters:\n",
"\n",
"$J=\\left[\\begin{array}{cc}\n",
"J_{0} & J_{1}\\exp\\left(i\\delta_{1}\\right)\\\\\n",
"J_{2}\\exp\\left(i\\delta_{2}\\right) & J_{3}\\exp\\left(i\\delta_{3}\\right)\n",
"\\end{array}\\right]$.\n",
"\n",
"There are several techniques to determine the 7 remaining parameters of the Jones matrix. Here we are going to implement the method described in \"Calibration method to determine the complete Jones matrix of SLMs\", J. Hoyo, L.M. Sanchez-Brea and A. Soria-Garcia, Optics and Lasers in Engineering, Volume 151, April 2022, 106914. This method requires performing 10 intensity measurements, varying the rotation angles of the elements of the PSG and PSA, which we can characterize using the bra-ket notation [1]:\n",
"\n",
"1. $I_{1}=\\left|\\left\\langle 0\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=J_{0}^{2}$,\n",
"1. $I_{2}=\\left|\\left\\langle 0\\text{º},0\\text{º}\\right|J_{sample}\\left|90\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=J_{1}^{2}$,\n",
"1. $I_{3}=\\left|\\left\\langle 90\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=J_{2}^{2}$,\n",
"1. $I_{4}=\\left|\\left\\langle 90\\text{º},0\\text{º}\\right|J_{sample}\\left|90\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=J_{3}^{2}$,\n",
"1. $I_{5}=\\left|\\left\\langle 0\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},-45\\text{º}\\right\\rangle \\right|^{2}=\\left(J_{0}^{2}+J_{1}^{2}+2J_{0}J_{1}\\sin\\delta_{1}\\right)/2$,\n",
"1. $I_{6}=\\left|\\left\\langle 0\\text{º},0\\text{º}\\right|J_{sample}\\left|45\\text{º, 0º}\\right\\rangle \\right|^{2}=\\left(J_{0}^{2}+J_{1}^{2}+2J_{0}J_{1}\\cos\\delta_{1}\\right)/2$,\n",
"1. $I_{7}=\\left|\\left\\langle 0\\text{º},45\\text{º}\\right|J_{sample}\\left|0\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=\\left(J_{0}^{2}+J_{2}^{2}+2J_{0}J_{2}\\sin\\delta_{2}\\right)/2$,\n",
"1. $I_{8}=\\left|\\left\\langle 45\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=\\left(J_{0}^{2}+J_{2}^{2}+2J_{0}J_{2}\\cos\\delta_{2}\\right)/2$,\n",
"1. $I_{9}=\\left|\\left\\langle 90\\text{º},0\\text{º}\\right|J_{sample}\\left|0\\text{º},-45\\text{º}\\right\\rangle \\right|^{2}=\\left[J_{2}^{2}+J_{3}^{2}+2J_{2}J_{3}\\sin\\left(\\delta_{3}-\\delta_{2}\\right)\\right]/2$,\n",
"1. $I_{10}=\\left|\\left\\langle 90\\text{º},0\\text{º}\\right|J_{sample}\\left|45\\text{º},0\\text{º}\\right\\rangle \\right|^{2}=\\left[J_{2}^{2}+J_{3}^{2}+2J_{2}J_{3}\\cos\\left(\\delta_{3}-\\delta_{2}\\right)\\right]/2$,\n",
"\n",
"where $\\left|\\phi,\\chi\\right\\rangle$ is the state generated by the PSG, $\\left\\langle \\phi,\\chi\\right|$ the state analyzed by the PSA (the state with maximum transmission), and $\\phi$ and $\\chi$ are the azimuth and ellipticity angle. Then, 7 parameters can be calculated from the intensity measurements:\n",
"\n",
"1. $J_{0}=\\sqrt{I_{1}}$,\n",
"1. $J_{1}=\\sqrt{I_{2}}$,\n",
"1. $J_{2}=\\sqrt{I_{3}}$,\n",
"1. $J_{3}=\\sqrt{I_{4}}$,\n",
"1. $\\delta_{1}=\\arctan\\left(\\frac{2I_{5}-I_{1}-I_{2}}{2I_{6}-I_{1}-I_{2}}\\right)$,\n",
"1. $\\delta_{2}=\\arctan\\left(\\frac{2I_{7}-I_{1}-I_{3}}{2I_{8}-I_{1}-I_{3}}\\right)$,\n",
"1. $\\delta_{3}=\\arctan\\left(\\frac{2I_{9}-I_{3}-I_{4}}{2I_{10}-I_{3}-I_{4}}\\right)+\\delta_{2}$.\n",
"\n",
"Let's simulate the experiment. The first stepp is implementing a function that calculates the intensity depending on the elements angles and another to calculate the matrix parameters from the intensities."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def sample_polarimeter_measurement(Jsample, angleDG, angleRG, angleRA, angleDA):\n",
" \"\"\"This function simulates the light polarimeter measurements.\n",
" \n",
" Parameters:\n",
" E: Incident light wave Jones vector.\n",
" angleDG: Rotating angle of the diattenuator (polarizer) of the PSG.\n",
" angleRG: Rotation angle of the quarter-wave plate of the PSG.\n",
" angleRA: Rotation angle of the quarter-wave plate of the PSA.\n",
" angleDA: Rotating angle of the diattenuator (polarizer) of the PSA.\n",
" \n",
" Returns:\n",
" I: Measured intensity.\"\"\"\n",
" # Rotate the optical elements\n",
" E = Jones_vector().circular_light(intensity=2) # Intensity = 2 produces a normalized light state of intensity 1 after the PSG\n",
" Jdg = Jones_matrix().diattenuator_perfect(azimuth=angleDG)\n",
" Jrg = Jones_matrix().quarter_waveplate(azimuth=angleRG)\n",
" Jra = Jones_matrix().quarter_waveplate(azimuth=angleRA)\n",
" Jda = Jones_matrix().diattenuator_perfect(azimuth=angleDA)\n",
" # Multiply all objects\n",
" E_final = Jda * Jra * Jsample * Jrg * Jdg * E\n",
" # Measure the intensity\n",
" I = E_final.parameters.intensity()\n",
" # Return\n",
" return I\n",
"\n",
"def matrix_parameters(I):\n",
" \"\"\"This function calculates the parameters of the Jones matrix from the intensity measurements.\n",
" \n",
" Parameters:\n",
" I: Array of intensities.\n",
" \n",
" Returns:\n",
" J0, J1, J2, J3, d1, d2, d3: Matrix parameters.\"\"\"\n",
" # Calculate modules\n",
" J0, J1, J2, J3 = np.sqrt(I[:4])\n",
" # Calculate phases\n",
" d1 = np.arctan2(2*I[4] - I[1] - I[0], 2*I[5] - I[1] - I[0]) % (360*degrees) \n",
" d2 = np.arctan2(2*I[6] - I[2] - I[0], 2*I[7] - I[2] - I[0]) % (360*degrees) \n",
" d3 = (np.arctan2(2*I[8] - I[2] - I[3], 2*I[9] - I[2] - I[3]) + d2) % (360*degrees) \n",
" # Return\n",
" return J0, J1, J2, J3, d1, d2, d3 "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can check if it works."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Comparison between original and calculated parameters:\n",
"- J0: Original = 0.0666; Calculated = 0.0666\n",
"- J1: Original = 0.6734; Calculated = 0.6734\n",
"- J2: Original = 0.6419; Calculated = 0.6419\n",
"- J3: Original = 0.5885; Calculated = 0.5885\n",
"- d1 (deg): Original = 134.6173; Calculated = 134.6173\n",
"- d2 (deg): Original = 181.8143; Calculated = 181.8143\n",
"- d3 (deg): Original = 274.9632; Calculated = 274.9632\n"
]
}
],
"source": [
"# Generate a random Jones matrix for the sample\n",
"J0, J1, J2, J3 = np.random.rand(4)\n",
"d1, d2, d3 = np.random.rand(3) * 360*degrees\n",
"Jsample = Jones_matrix().from_components([J0, J1*np.exp(1j*d1), J2*np.exp(1j*d2), J3*np.exp(1j*d3)])\n",
"\n",
"# Calculate the intensities\n",
"angleDG = np.array([0, 90, 0, 90, 45, 45, 0, 0, 45, 45]) * degrees\n",
"angleRG = np.array([0, 90, 0, 90, 0, 45, 0, 0, 0, 45]) * degrees\n",
"angleRA = np.array([0, 0, 90, 90, 0, 0, 0, 45, 90, 90]) * degrees\n",
"angleDA = np.array([0, 0, 90, 90, 0, 0, 45, 45, 90, 90]) * degrees\n",
"I = sample_polarimeter_measurement(Jsample, angleDG, angleRG, angleRA, angleDA)\n",
"\n",
"# Calculate the matrix parameters\n",
"J0_calc, J1_calc, J2_calc, J3_calc, d1_calc, d2_calc, d3_calc = matrix_parameters(I)\n",
"\n",
"# Compare results\n",
"print(\"Comparison between original and calculated parameters:\")\n",
"original = [J0, J1, J2, J3, d1/degrees, d2/degrees, d3/degrees]\n",
"calculated = [J0_calc, J1_calc, J2_calc, J3_calc, d1_calc/degrees, d2_calc/degrees, d3_calc/degrees]\n",
"parameter = [\"J0\", \"J1\", \"J2\", \"J3\", \"d1 (deg)\", \"d2 (deg)\", \"d3 (deg)\"]\n",
"for ind in range(7):\n",
" print(\"- {}: Original = {:.4f}; Calculated = {:.4f}\".format(parameter[ind], original[ind], calculated[ind]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The match is perfect."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mueller-Stokes formalism\n",
"\n",
"The method for determining the Mueller matrix of an optical element is similar to the one for Stokes light polarimetry. The k-th measurement of the light polarimeter experiment is:\n",
"\n",
"$S_{f}^{k}=\\left[\\begin{array}{cccc}\n",
"A_{00}^{k} & A_{01}^{k} & A_{02}^{k} & A_{03}^{k}\\\\\n",
"\\bullet & \\bullet & \\bullet & \\bullet\\\\\n",
"\\bullet & \\bullet & \\bullet & \\bullet\\\\\n",
"\\bullet & \\bullet & \\bullet & \\bullet\n",
"\\end{array}\\right]\\left[\\begin{array}{cccc}\n",
"M_{00} & M_{01} & M_{02} & M_{03}\\\\\n",
"M_{10} & M_{11} & M_{12} & M_{13}\\\\\n",
"M_{20} & M_{21} & M_{22} & M_{23}\\\\\n",
"M_{30} & M_{31} & M_{32} & M_{33}\n",
"\\end{array}\\right]\\left[\\begin{array}{c}\n",
"S_{0}^{k}\\\\\n",
"S_{1}^{k}\\\\\n",
"S_{2}^{k}\\\\\n",
"S_{3}^{k}\n",
"\\end{array}\\right]=\\left[\\begin{array}{c}\n",
"I^{k}\\\\\n",
"\\bullet\\\\\n",
"\\bullet\\\\\n",
"\\bullet\n",
"\\end{array}\\right]$.\n",
"\n",
"Being $A^{k}=M_{D2}(\\theta_{D2}^{k})*M_{R2}(\\theta_{R2}^{k})$ the Mueller matrix of the PSA, $S^k=M_{R1}(\\theta_{R1}^{k})*M_{D1}(\\theta_{D1}^{k})*S_i$ the Stokes vector generated by the PSG and $M_{ij}$ the elements of the Mueller matrix of the sample. Due to the fact that the detector only measures the intensity, only some parts of the Mueller matrix of the PSA are important (the rest of the parameters are replaced by $\\bullet$). Then, if we consider the 16 required measurements:\n",
"\n",
"$\\left[\\begin{array}{ccc}\n",
"A_{00}^{0}S_{0}^{0} & \\ldots & A_{03}^{0}S_{3}^{0}\\\\\n",
"\\vdots & \\ddots & \\vdots\\\\\n",
"A_{00}^{16}S_{0}^{16} & \\ldots & A_{03}^{16}S_{3}^{16}\n",
"\\end{array}\\right]\\left[\\begin{array}{c}\n",
"M_{00}\\\\\n",
"\\vdots\\\\\n",
"M_{44}\n",
"\\end{array}\\right]=\\left[\\begin{array}{ccc}\n",
"b_{0,0} & \\ldots & b_{0,16}\\\\\n",
"\\vdots & \\ddots & \\vdots\\\\\n",
"b_{16,0} & \\ldots & b_{16,16}\n",
"\\end{array}\\right]\\left[\\begin{array}{c}\n",
"M_{00}\\\\\n",
"\\vdots\\\\\n",
"M_{44}\n",
"\\end{array}\\right]=\\left[\\begin{array}{c}\n",
"I^{0}\\\\\n",
"\\vdots\\\\\n",
"I^{16}\n",
"\\end{array}\\right]$.\n",
"\n",
"If $b$ is invertible ($\\det(b)\\neq0$) then $M_{ij}$ can be easily calculated as $M_{ij}=b^{-1}I'$. This condition is satisfied if there are no repeated set of four angles.\n",
"\n",
"Let us reproduce this numerically. Let us define a function that simulates the experiment and calculates the $b$ matrix. Also, we will include an error amplitude which we will use in the future."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def Mueller_polarimeter_measurement(Msample, Mpol, S_ini, angles, Ierror=0):\n",
" \"\"\"This function simulates the light polarimeter measurements.\n",
" \n",
" Parameters:\n",
" Msample (Mueller): Mueller matrix of the sample.\n",
" Mpol (list): List with the four Mueller objects of the polarimeter: Md1, Mr1, Mr2, Md2.\n",
" Sini (Stokes): Stokes object of the illumination.\n",
" angles (numpy.ndarray): 4xN array with the angles used to perform the measurement.\n",
" Ierror (float): Intensity error amplitude.\n",
" \n",
" Returns:\n",
" Ik: Measured intensities vector.\n",
" b: b matrix.\"\"\"\n",
" # Extract the information from the input parameters\n",
" N = angles.shape[1]\n",
" Md1, Mr1, Mr2, Md2 = Mpol\n",
" # Rotate the optical elements\n",
" Md1_rotated = Md1.rotate(angle=angles[0,:], keep=True)\n",
" Mr1_rotated = Mr1.rotate(angle=angles[1,:], keep=True)\n",
" Mr2_rotated = Mr2.rotate(angle=angles[2,:], keep=True)\n",
" Md2_rotated = Md2.rotate(angle=angles[3,:], keep=True)\n",
" # Multiply the objects as required\n",
" Sg = Mr1_rotated * Md1_rotated * S_ini\n",
" A = Md2_rotated * Mr2_rotated\n",
" S_final = A * Msample * Sg\n",
" # Measure the intensity\n",
" Ik = S_final.parameters.intensity() + np.random.randn(N) * Ierror\n",
" # Calculate the b matrix\n",
" Ak = A.parameters.matrix() # Extract the 4x4xN matrix of the Mueller object\n",
" Ak = Ak[0,:,:] # Extract all the first rows\n",
"# Ak = np.transpose(Ak) # Transpose so Ak rows corrrespond to different measurements\n",
" Sk = Sg.parameters.matrix() # Extract the 4xN matrix of the Stokes object\n",
"# Sk = np.transpose(Sk) # Transpose so Sk rows corrrespond to different measurements\n",
" b = np.zeros((N,16))\n",
" for ind in range(N):\n",
" b[ind,:] = np.outer(Ak[:,ind], Sk[:,ind]).flatten()\n",
" # Return\n",
" return Ik, b"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we can perform the experiment and see if the calculated Mueller matrix is correct. We will use as sample the product of a random diattenuator, a random polarizer and a random diagonal diattenuator. We will use a circularly polarized light wave for simplicity."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Sample = \n",
"[+0.508 -0.018 +0.093 -0.179] \n",
"[-0.407 +0.022 -0.091 +0.230] \n",
"[-0.038 +0.026 -0.069 -0.070] \n",
"[+0.056 +0.016 +0.136 -0.004] \n",
" Calculated = \n",
"[+0.508 -0.018 +0.093 -0.179] \n",
"[-0.407 +0.022 -0.091 +0.230] \n",
"[-0.038 +0.026 -0.069 -0.070] \n",
"[+0.056 +0.016 +0.136 -0.004] \n",
"\n"
]
}
],
"source": [
"# Start by defining the polarimeter system\n",
"S_ini = Stokes('Light source')\n",
"S_ini.circular_light(intensity=1)\n",
"Mpol = create_Mueller(('Diattenuator 1', 'Retarder 1', 'Retarder 2', 'Diattenuator 2'))\n",
"Mpol[0].diattenuator_perfect(azimuth=0)\n",
"Mpol[1].quarter_waveplate(azimuth=0)\n",
"Mpol[2].quarter_waveplate(azimuth=0)\n",
"Mpol[3].diattenuator_perfect(azimuth=0)\n",
"\n",
"# Create the sample parameters\n",
"p1 = np.random.rand() * 0.5 + 0.5\n",
"p2 = np.random.rand() * 0.5\n",
"alphaD = np.random.rand() * 90*degrees\n",
"delayD = np.random.rand() * 360*degrees\n",
"\n",
"R = np.random.rand() * 180*degrees\n",
"alphaR = np.random.rand() * 90*degrees\n",
"delayR = np.random.rand() * 360*degrees\n",
"\n",
"d = [np.random.rand(), np.random.rand(), np.random.rand()]\n",
"\n",
"# Create the sample \n",
"Md, Mr, Mp = create_Mueller(N=3)\n",
"Md.diattenuator_charac_angles(p1=p1, p2=p2, alpha=alphaD, delay=delayD)\n",
"Mr.retarder_charac_angles(R=R, alpha=alphaR, delay=delayR)\n",
"Mp.depolarizer_diagonal(d=d)\n",
"\n",
"Msample = Md*Mr*Mp\n",
"Msample.name = 'Sample'\n",
"\n",
"# Now, define the rotation angles for the optical elements\n",
"angles = np.random.rand(4,16) * 180*degrees\n",
"\n",
"# Make the measurements\n",
"Ik, b = Mueller_polarimeter_measurement(Msample=Msample, Mpol=Mpol, S_ini=S_ini, angles=angles)\n",
"\n",
"# Calculate the original Mueller matrix\n",
"b_inv = np.linalg.inv(b)\n",
"Ik = np.reshape(Ik, (16,1))\n",
"M_calc = Mueller('Calculated')\n",
"M_calc.from_matrix(M=b_inv@Ik)\n",
"\n",
"# Compare the result\n",
"print(Msample, M_calc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As explained in the previous example, experiments are subject to errors. This means that the calculated Mueller matrices will be wrong. However, the error in the calculated matrix can be reduced by oversampling, i.e., taking more than the minimum 16 measurements. In this case, $b$ will be a Nx16 matrix and $I'$ a Nx1 vector. Then, $M_{ij}$ can be calculated using the pseudoinverse of $b$: \n",
"\n",
"$b^+=(b^T*b)^{-1}*b^T$.\n",
"\n",
"Let us define a function that performs one polarimetry experiment."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def Mueller_polariemter_experiment(N, Ierror):\n",
" \"\"\"This funcion simulates a polarimetry experiment with errors in the detection.\n",
" \n",
" Parameters:\n",
" N (int): Number of measurements.\n",
" Ierror (float): Intensity error amplitude.\n",
" \n",
" Returns:\n",
" dif (numpy.ndarray): Difference between the calculated and the original array.\"\"\"\n",
" # Create the sample parameters\n",
" p1 = np.random.rand() * 0.5 + 0.5\n",
" p2 = np.random.rand() * 0.5\n",
" alphaD = np.random.rand() * 90*degrees\n",
" delayD = np.random.rand() * 360*degrees\n",
"\n",
" R = np.random.rand() * 180*degrees\n",
" alphaR = np.random.rand() * 90*degrees\n",
" delayR = np.random.rand() * 360*degrees\n",
"\n",
" d = [np.random.rand(), np.random.rand(), np.random.rand()]\n",
"\n",
" # Create the sample \n",
" Md, Mr, Mp = create_Mueller(N=3)\n",
" Md.diattenuator_charac_angles(p1=p1, p2=p2, alpha=alphaD, delay=delayD)\n",
" Mr.retarder_charac_angles(R=R, alpha=alphaR, delay=delayR)\n",
" Mp.depolarizer_diagonal(d=d)\n",
"\n",
" Msample = Md*Mr*Mp\n",
" Msample.name = 'Sample'\n",
"\n",
" # Now, define the rotation angles for the optical elements\n",
" angles = np.random.rand(4,N) * 180*degrees\n",
"\n",
" # Make the measurements\n",
" Ik, b = Mueller_polarimeter_measurement(Msample=Msample, Mpol=Mpol, S_ini=S_ini, angles=angles, Ierror=Ierror)\n",
"\n",
" # Calculate the Stokes vector\n",
" b_T = np.transpose(b)\n",
" b_inv = np.linalg.inv(b_T @ b)\n",
" b_inv = b_inv @ b_T\n",
" Ik = np.reshape(Ik, (N,1))\n",
" M_calc = b_inv @ Ik\n",
" \n",
" # Calculate the difference\n",
" dif = np.linalg.norm(Msample.M.flatten() - M_calc.flatten())\n",
"\n",
" # Return\n",
" return dif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, lets do the experiment many times, changing the number of measurements and the error amplitude. "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# Define some variables\n",
"Nstat = 100 # Number of experiments done under the same condition.\n",
"I_errors = np.array([0, 0.005, 0.01, 0.02]) # Array of error factors\n",
"Ns = np.arange(20, 100, 10) # Array of the number of measurements\n",
"tol = 1e-12\n",
"I_wave = 1 # We will maintain this parameter constant so we can compare the errors easily\n",
"\n",
"# Create the result variables\n",
"error = np.zeros((Ns.size, I_errors.size))\n",
"legend = []\n",
"\n",
"# Make the loops changing the variables\n",
"for indI, Ierror in enumerate(I_errors): # Loop in the error amplitude\n",
" for indN, N in enumerate(Ns): # Loop in the number of measurements\n",
" # Create the temporal variables\n",
" aux = np.zeros((Nstat))\n",
" # Repeat the same conditions several times\n",
" for indS in range(Nstat):\n",
" # Do the experiment\n",
" dif = Mueller_polariemter_experiment(N=N, Ierror=Ierror);\n",
" # Calculate the errors\n",
" aux[indS] = dif\n",
" # Calculate the mean square error\n",
" error[indN, indI] = np.linalg.norm(aux) / Nstat\n",
" # Create the legend for the plots\n",
" legend += ['Ierror = {} %'.format(Ierror*100)]"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'MSE')"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAHgCAYAAABJrX+JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABXTUlEQVR4nO3dd5xU9b3/8ddnZmfLDAssvTepFkREYBeDsUWxRqOyaKIxJsZo6v2lx3tzTa/3GnOTGGPXUGyxRYEYCyKggAoWQJC6S6/b+/f3x5ldtrIL7OyZmX0/H4957Mw537PzOaC8z/mec75fc84hIiIiySXgdwEiIiLS/hTwIiIiSUgBLyIikoQU8CIiIklIAS8iIpKEFPAiIiJJKMXvAtpTr1693LBhw/wuQ0REpEOsXLlyr3Oud3Prkirghw0bxooVK/wuQ0REpEOY2ZaW1qmLXkREJAkp4EVERJKQAl5ERCQJJdU1eBEROXqVlZXk5eVRVlbmdynSgvT0dAYNGkQoFGrzNgp4EZFOLi8vj8zMTIYNG4aZ+V2ONOKcY9++feTl5TF8+PA2b6cuehGRTq6srIyePXsq3OOUmdGzZ8+j7mGJacCb2YVmts7MNpjZ95tZb2Z2V3T9ajObWG9ddzN7wszWmtkaM8uOZa0iIp2Zwj2+HcvfT8wC3syCwJ+AGcCJwCwzO7FRsxnAqOjrZuAv9db9AZjvnBsLnAqsiVWtIiLiry5duvhdQqucc3z9619n5MiRjB8/nrfffrvZdps2bWLKlCmMGjWKmTNnUlFR0aTNG2+8wfjx4znjjDPYsGEDAAcPHuSCCy7AOdcu9cbyDH4ysME5t9E5VwHMBS5v1OZy4GHnWQZ0N7P+ZtYVmA7cB+Ccq3DOHYxhrSIikiCcc9TU1LT4uSXV1dXH9b0vvvgi69evZ/369dxzzz185Stfabbd9773Pb71rW+xfv16srKyuO+++5q0+f3vf8+TTz7JL37xC/7yF+/c9qc//Sk//OEP2603JZYBPxDYVu9zXnRZW9qMAPYAD5jZO2Z2r5lFmvsSM7vZzFaY2Yo9e/a0X/UiIuKL3/72t5xxxhmMHz+eH//4xwBs3ryZcePGceuttzJx4kRef/31Bp+3bdvGd77zHU4++WROOeUU5s2bB8Crr77K2WefzbXXXsspp5xyXHU988wzXH/99ZgZU6dO5eDBg+zYsaNBG+ccL7/8MldddRUAN9xwA08//XST3xUKhSgtLaWkpIRQKMTHH39Mfn4+Z5111nHVWF8s76Jv7hCkcb9DS21SgInA15xzb5rZH4DvA//ZpLFz9wD3AEyaNKl9+jVERDqpO577gA+3F7Tr7zxxQFd+fOlJbWq7cOFC1q9fz1tvvYVzjssuu4xFixYxZMgQ1q1bxwMPPMCf//xnNm/e3ODzk08+ybvvvsuqVavYu3cvZ5xxBtOnTwfgrbfe4v3332/2DvSZM2eybt26Jsv/4z/+g+uvv77Bsvz8fAYPHlz3edCgQeTn59O/f/+6Zfv27aN79+6kpKQ0aNPYD37wA26++WYyMjJ45JFH+Pa3v81Pf/rTNv0ZtVUsAz4PGFzv8yBgexvbOCDPOfdmdPkTeAEvIiJJbOHChSxcuJDTTjsNgKKiItavX8+QIUMYOnQoU6dOrWtb//PixYuZNWsWwWCQvn37ctZZZ7F8+XK6du3K5MmTW3y8rPZMvy2auzbeuDu9LW0AJkyYwLJlywBYtGgRAwYMwDnHzJkzCYVC/P73v6dv375trq05sQz45cAoMxsO5AO5wLWN2jwLfNXM5gJTgEPOuR0AZrbNzMY459YB5wIfxrBWERGBNp9px4pzjh/84Ad8+ctfbrB88+bNRCINr9TW/3ykG9Mab1ff0ZzBDxo0iG3bDl9VzsvLY8CAAQ3a9OrVi4MHD1JVVUVKSkqzbepzzvGzn/2MefPm8dWvfpU77riDzZs3c9ddd/Hzn/+8xe3aImbX4J1zVcBXgQV4d8A/5pz7wMxuMbNbos1eADYCG4C/AbfW+xVfA/5uZquBCcAvYlWriIjEhwsuuID777+foqIiwOsW3717d6vbTZ8+nXnz5lFdXc2ePXtYtGgRkydPbnW7efPm8e677zZ5NQ53gMsuu4yHH34Y5xzLli2jW7duDbrnwTtbP/vss3niiScAeOihh7j88sb3lx/20EMPcfHFF5OVlUVJSQmBQIBAIEBJSUmrtbcmpiPZOedewAvx+svurvfeAbe1sO27wKRY1iciIvHlU5/6FGvWrCE72xv6pEuXLjz66KMEg8EjbnfFFVewdOlSTj31VMyM3/zmN/Tr14+1a9e2W20XXXQRL7zwAiNHjiQcDvPAAw80WHfvvfcyYMAAfv3rX5Obm8vtt9/Oaaedxk033dTs7yspKeGhhx5i4cKFgNdr8JnPfIbU1FTmzJlz3PVaez1vFw8mTZrkNB+8iMjRWbNmDePGjfO7DGlFc39PZrbSOdfsybCGqm1BTXk5VQcO+F2GiIjIMVHAN8PV1LDhnHPZ84c/+F2KiIjIMVHAN8MCATJOOYXiJUv9LkVEROSYKOBbEMnJoXLrViry8vwuRURE5Kgp4FsQmZYDQPGSJT5XIiIicvQU8C1IHTGClD591E0vIiIJSQHfAjMjkpNDydKluDbMUiQiIscumaaL/fznP8/w4cOZMGECEyZM4N13323SJtGni014kWk5VB86RNmHmopeRCRexPt0seDNiFc7Kt6ECROarE/06WITXiQ6iYGuw4uIdJxEni62rRJ9utiEl9K7N2mjR1O8dAm9bv6S3+WIiMTei9+Hne+17+/sdwrM+FWbmib6dLG1fvSjH/GTn/yEc889l1/96lekpaU1WJ/o08UmhUhODgdmz6amrIxAerrf5YiIJLVEny4W4Je//CX9+vWjoqKCm2++mV//+tf813/9V4M2iT5dbFKI5GSz/8EHKVmxki5nTvO7HBGR2GrjmXasJPp0sUDdGX1aWho33ngjv/vd71r8/oScLjZZhCdNwkIhipfqOryISKwl+nSxQN11eeccTz/9NCeffHKL35+w08Umg0A4TMZpp+l5eBGRDpAM08Ved9117NmzB+ccEyZM4O67727292m62KMQq+li9979V/bceSej3lhMSs+e7f77RUT8pOliE4Omi42BSI53JFkcvSFCREQk3ing2yD9pJMIdOum5+FFRCRhKODbwIJBIlOmULxkabsNISgiIhJLCvg2iuRkU7VjBxWbN/tdioiISKsU8G0UydH0sSIikjgU8G2UOmQIoUGD9LiciIgkBAX8UYhkZ1Py5pu4qiq/SxERSSqJMF3s2rVryc7OJi0t7Yij023atIkpU6YwatQoZs6cSUVFRZM2mi42zkSm5VBTVETpe+08EYOIiLSZX9PF9ujRg7vuuotvf/vbR2z3ve99j29961usX7+erKws7rvvviZtNF1snAlPmQJmug4vIhJD8TpdbJ8+fTjjjDMIhUIttnHO8fLLL3PVVVcBcMMNN/D00083aafpYuNMSlYW6SeeSPHSpfS+7Ta/yxERaXe/fuvXrN3ffsO7AoztMZbvTf5em9rG83SxbbFv3z66d+9OSooXr7VTyjam6WLjUCQnh30PPEB1UTHBLi3PUCQiIkcvnqeLbYu2Timr6WLjUGRaDvv+9jdKlr9F5tln+12OiEi7auuZdqzE83SxbdGrVy8OHjxIVVUVKSkpLU4pW0vTxcaRjNNOw9LSKF6qx+VERNpbPE8X2xZmxtlnn80TTzwBeNPBXn755S2213SxcSSQlkZ40iTdaCciEgPxPF3szp07mTRpEgUFBQQCAe68804+/PBDunbt2mC62F//+tfk5uZy++23c9ppp3HTTTc1+/s0XexRiNV0sY3tu+9+dv/2t4x87VVCx3mNRETEb5ouNjFoutgOUDd9rLrpRUQkTingj0HamDEEe/RQN72IiMQtBfwxsECASHY2xUs1fayIiMQnBfwxiuRkU71nL+Xr1/tdioiISBMK+GOk6WNFRCSeKeCPUah/f1KHD1fAi4hIXFLAH4dIdjYly1fgmpkKUERE2i4Rpov9+9//zvjx4xk/fjw5OTmsWrWq2XaaLjYJRKbl4EpLKXn3Xb9LERHpNPyaLnb48OG89tprrF69mv/8z//k5ptvbradpotNAuHJkyEYVDe9iEg7itfpYnNycsjKygJg6tSp5OXlNWmj6WKTRDAzk4xTTqF4yVL45jf9LkdE5Ljt/MUvKF/TvtPFpo0bS78f/rBNbRNlutj77ruPGTNmNFmu6WKTSCQnh7133031oUMEu3XzuxwRkYSWCNPFvvLKK9x3330sXry4yTpNF5tEItNy2PvnP1P85pt0/dSn/C5HROS4tPVMO1bifbrY1atX88UvfpEXX3yRnj17Nlmv6WKTSMb48QTCYV2HFxFpB/E8XezWrVu58soreeSRRxg9enSzv0/TxSYRC4UIT56siWdERNpBPE8X+5Of/IR9+/Zx6623ApCSkkLtDKaaLjbGOmq62Mb2P/wIu37xC0546SVSBw3s8O8XETkemi42MWi6WB/UTR+75A2fKxEREfEo4NtB6gknkNKnj7rpRUQkbijg24GZEcnJoWTpMlwbRlMSERGJNQV8O4nkZFN98CBlH67xuxQRkaOWTPdjJaNj+ftRwLeTSPSOz+KlelxORBJLeno6+/btU8jHKecc+/btIz09/ai202Ny7SSld2/SRo+meMkSen3pS36XIyLSZoMGDSIvL489e/b4XYq0ID09nUGDBh3VNgr4dhTJzubAnDnUlJUROMojLRERv4RCoRaHcpXEpS76dhSZloOrqKBk5Uq/SxERkU5OAd+OwpMmQSikYWtFRMR3Cvh2FAiHCU+Y4E0fKyIi4iMFfDuLTMuhfM0aqvbv97sUERHpxBTw7SySkwOgUe1ERMRXMQ14M7vQzNaZ2QYz+34z683M7oquX21mE+ut22xm75nZu2bW8TPIHKP0k04i0LWrrsOLiIivYvaYnJkFgT8B5wN5wHIze9Y592G9ZjOAUdHXFOAv0Z+1znbO7Y1VjbFgwSCRKVMoXrIU5xxm5ndJIiLSCcXyDH4ysME5t9E5VwHMBRrPen858LDzLAO6m1n/GNbUISLTcqjasYOKzZv9LkVERDqpWAb8QGBbvc950WVtbeOAhWa20sxujlmVMVA3bK266UVExCexDPjm+qYbD3R8pDbTnHMT8brxbzOz6c1+idnNZrbCzFbEyzCLoSFDCA0cqBvtRETEN7EM+DxgcL3Pg4DtbW3jnKv9uRv4B16XfxPOuXucc5Occ5N69+7dTqUfn7rpY5e9iauq8rscERHphGIZ8MuBUWY23MxSgVzg2UZtngWuj95NPxU45JzbYWYRM8sEMLMI8Cng/RjW2u4iOdnUFBVR+t57fpciIiKdUMzuonfOVZnZV4EFQBC43zn3gZndEl1/N/ACcBGwASgBboxu3hf4R/QO9BRgtnNufqxqjYXw1KlgRvHSpYRPO83vckREpJOxZJr/d9KkSW7Fivh5ZH7TZ67CMtIZ9uijfpciIiJJyMxWOucmNbdOI9nFUCQnm9J3V1FdVOx3KSIi0sko4GMokpMDVVWUrFjudykiItLJKOBjKGPiRCwtTc/Di4hIh1PAx1AgLY3w6acr4EVEpMMp4GMsMi2Hig0fU7lrt9+liIhIJ6KAj7HD08fqLF5ERDqOAj7G0saMIdijh7rpRUSkQyngY8wCASJTp1K81Js+VkREpCMo4DtAZFoO1Xv2Ur5+vd+liIhIJ6GA7wCaPlZERDqaAr4DhAYMIHXYME0fKyIiHUYB30EiOTmULF+Bq6jwuxQREekEFPAdJDItB1dSQsm77/pdioiIdAIK+A4SnjwZgkF104uISIdQwHeQYGYmGaecohvtRESkQyjgO1AkJ4ey996n+tAhv0sREZEkp4DvQJGcbKipofitt/wuRUREkpwCvgNlnHoqgXBY3fQiIhJzCvgOZKEQ4cmTFfAiIhJzCvgOFsnJpnLLViry8v0uRUREkpgCvoNp+lgREekICvgOlnrCCaT06aNuehERiSkFfAczMyLZ2ZQsXYarqfG7HBERSVIKeB9EpuVQffAgZWvW+F2KiIgkKQW8DzR9rIiIxJoC3gcpvXuTNmoUJRqXXkREYkQB75NITg4lK1ZSU1bmdykiIpKEFPA+iUzLwVVUULJypd+liIhIElLA+yQ8aRKEQuqmFxGRmFDA+yQQDhOeMIEi3WgnIiIxoID3UWRaDuUfrqFq/36/SxERkSSjgPdR7eNyJcuW+VyJiIgkGwW8j9JPPplA167qphcRkXangPeRBYNEpkyheMkSnHN+lyMiIklEAe+zSE42Vdt3ULlli9+liIhIElHA+6x2+lh104uISHtSwPssNGQIoYEDNS69iIi0KwW8z8yMSE42JW++hauq8rscERFJEgr4OBDJyaGmsJCy99/3uxQREUkSCvg4EJ46Fcx0HV5ERNqNAj4OpGRlkT5uHCVLNC69iIi0DwV8nIhMy6Fk1Spqiov9LkVERJKAAj5ORHJyoLKS4uXL/S5FRESSgAI+TmRMnIilpelxORERaRcK+DgRSEsjfPrpmh9eRETahQI+jkSm5VC+fgOVu3b7XYqIiCQ4BXwcqZ0+tnipuulFROT4KODjSNrYsQR79FA3vYiIHDcFfByxQIDI1KkUL1mq6WNFROS4KODjTCQnm6o9eyhfv97vUkREJIEp4ONM7fSx6qYXEZHjoYCPM6EBA0gdNkzj0ouIyHFRwMehSE42JctX4Coq/C5FREQSlAI+DkVycnAlJZSuWuV3KSIikqAU8HEoPGUKBALqphcRkWOmgI9DwcxMMk45RePSi4jIMVPAx6nItBzK3nuf6oICv0sREZEEpICPU5GcHKipofjNN/0uRUREElBMA97MLjSzdWa2wcy+38x6M7O7outXm9nERuuDZvaOmT0fyzrjUcb48Vg4rG56ERE5JjELeDMLAn8CZgAnArPM7MRGzWYAo6Kvm4G/NFr/DWBNrGqMZ5aaSuSMMyhZogFvRETk6MXyDH4ysME5t9E5VwHMBS5v1OZy4GHnWQZ0N7P+AGY2CLgYuDeGNca1yLQcKrZsoTI/3+9SREQkwcQy4AcC2+p9zosua2ubO4HvAjUxqi/u1Q5bq8flRETkaMUy4K2ZZY2nSGu2jZldAux2zq1s9UvMbjazFWa2Ys+ePcdSZ9xKPeEEUvr00bj0IiJy1GIZ8HnA4HqfBwHb29hmGnCZmW3G69o/x8webe5LnHP3OOcmOecm9e7du71qjwtmRiQ7m+Kly3A1nbYjQ0REjkEsA345MMrMhptZKpALPNuozbPA9dG76acCh5xzO5xzP3DODXLODYtu97Jz7rMxrDVuRablUH3gAGVrOuW9hiIicoxiFvDOuSrgq8ACvDvhH3POfWBmt5jZLdFmLwAbgQ3A34BbY1VPogpPnQpo+lgRETk65lzjy+KJa9KkSW7FihV+l9HuNl56GSm9ezHk/vv9LkVEROKIma10zk1qbp1GsksAkZwcSlaspKaszO9SREQkQSjgE0AkJxtXUUHp22/7XYqIiCQIBXwCCJ9xBoRCGrZWRETaTAGfAALhMOEJEzTgjYiItJkCPkFEcrIp/3ANVQcO+F2KiIgkAAV8gqgdtlaPy4mISFso4BNE+sknE8jMVDe9iIi0iQI+QVgwSGTqFIqXLCGZxi4QEZHYUMAnkEhODlXbd1C5ZYvfpYiISJxTwCcQTR8rIiJtpYBPIKEhQwgNGKAb7UREpFUK+ARiZkSm5VC87E1cVZXf5YiISBxTwCeYSE4ONYWFlL3/vt+liIhIHFPAJ5jw1KlgRrG66UVE5AgU8AkmJSuL9HHjKH5DN9qJiEjLFPAJKDIth5JVq6gpLva7FBERiVMK+AQUyc6GykpKVqzwuxQREYlTCvgElHH66VhamqaPFRGRFingE1AgLY3w6acr4EVEpEUK+AQVycmmfP0GKnfv9rsUERGJQwr4BKXpY0VE5EgU8AkqbexYgllZ6qYXEZFmKeATlAUCRLKnUrxkqaaPFRGRJhTwCSySk0PVnj1UbNjgdykiIhJnFPAJrPY6vLrpRUSkMQV8AgsNGEDq0KEUL9GNdiIi0pACPsFFpuVQvHw5rqLC71JERCSOKOATXCQnB1dSQumqVX6XIiIicUQBn+DCkydDIKDpY0VEpAEFfIILdu1KximnaPpYERFpQAGfBCLTcih97z2qCwr8LkVEROKEAj4JRLKzoaaGkrfe8rsUERGJEwr4JJBx6qlYOKzn4UVEpI4CPglYaiqRM87QdXgREamjgE8SkZxsKrZsoTI/3+9SREQkDijgk0TdsLV6XE5ERFDAJ43UkSNJ6d1b1+FFRARQwB9RIk3DamZEcrIpXroMV1PjdzkiIuIzBXwzKqsr+cpLX+Ge1ff4XcpRieTkUH3gAOVr1/pdioiI+EwB34xQMIRzjsfWPUZlTaXf5bRZODsb0PSxIiKigG9R7thcdpfu5uWtL/tdSpuF+vQhbdRITR8rIiIK+JZ8YuAnGBAZwNy1c/0u5ahEcnIoWbmSmvJyv0sREREfKeBbEAwEuWbMNazYtYL1B9b7XU6bRXJycOXllK5c6XcpIiLiIwX8EVw56kpSA6nMWzfP71LaLDxpEoRCeh5eRKSTU8AfQVZ6FhcOv5BnP36WwopCv8tpk0AkQvjUUzVsrYhIJ6eAb8WssbMorSrl2Y+f9buUNotMy6FszRqqDhzwuxQREfGJAr4VJ/c6mZN7nsy8dfMSZuCbSHY2OEfJsmV+lyIiIj5RwLdB7thcNh3axJs73/S7lDZJP/lkApmZeh5eRKQTU8C3wYXDL6R7WveEeWTOUlKITJ1C8RtLEqbXQURE2pcCvg3SgmlcOepKXtn2CjuKdvhdTpuEs7Op3L6dyi1b/C5FRER8oIBvo2vGXINzjsc/etzvUtqki6aPFRHp1BTwbTSwy0DOGnQWT65/korqCr/LaVVo6FBCAwboOryISCelgD8KuWNz2V+2n4VbFvpdSqvMjMi0HIqXvYmrqvK7HBER6WAK+KOQPSCboV2HJszNdpHsbGoKCyn74AO/SxERkQ6mgD8KAQswc8xMVu1ZxYf7PvS7nFaFs7PBTN30IiKdkAL+KF0+8nIyUjIS4iw+JSuL9HHjNGytiEgnpIA/Sl1Tu3LR8It4YdMLHCo/5Hc5rYrkZFOyahU1xcV+lyIiIh1IAX8MZo2dRXl1OU9veNrvUloVycmBykpKVqzwuxQREelARwx4M/tsvffTGq37amu/3MwuNLN1ZrbBzL7fzHozs7ui61eb2cTo8nQze8vMVpnZB2Z2R9t3KfbG9BjDxD4Tmbt2LjWuxu9yjijj9NOxtDRdhxcR6WRaO4P/j3rv/9ho3ReOtKGZBYE/ATOAE4FZZnZio2YzgFHR183AX6LLy4FznHOnAhOAC81saiu1dqjcsbnkFeWxOH+x36UcUSAtjfDpEyleogFvREQ6k9YC3lp439znxiYDG5xzG51zFcBc4PJGbS4HHnaeZUB3M+sf/VwUbROKvuJqUPXzhpxHz/SeCXGzXSQnh/L166ncvdvvUkREpIO0FvCuhffNfW5sILCt3ue86LI2tTGzoJm9C+wG/uWci6up3ELBEFeNvorF+YvZVrit9Q18FIkOW1uiYWtFRDqN1gJ+bPTa+Hv13td+HtPKts2d4Tc+KGixjXOu2jk3ARgETDazk5v9ErObzWyFma3Ys2dPKyW1r6tHX03AAjy27rEO/d6jlTZ2LMGsLHXTi4h0IimtrB93HL87Dxhc7/MgYPvRtnHOHTSzV4ELgfcbf4lz7h7gHoBJkyZ1aDd+30hfzhlyDk+tf4pbJ9xKRkpGR359m1kgQCR7KsVLvOljzVq7uiIiIonuiGfwzrkt9V9AETAR6BX9fCTLgVFmNtzMUoFc4NlGbZ4Fro/eTT8VOOSc22Fmvc2sO4CZZQDnAWuPeu86wKyxsyioKGD+pvl+l3JEkZwcqvbsoWLDBr9LERGRDtDaY3LP13aNm1l/vDPoLwCPmNk3j7Stc64K+CqwAFgDPOac+8DMbjGzW6LNXgA2AhuAvwG3Rpf3B14xs9V4Bwr/cs49fwz7F3OT+k5iZPeRzFk7B+fi6j7ABiLZ2YCmjxUR6Sxa66If7pyr7Ra/ES9orzezTOAN4M4jbeycewEvxOsvu7veewfc1sx2q4HTWq0+DpgZM8fM5Odv/pzVe1dzau9T/S6pWaGBA0kdOpTiN5bQ4/rr/S5HRERirLWb7CrrvT+XaFg75wqB+B7hpQNdesKlREKRuH9kLjIth+Lly3EV8T+fvYiIHJ/WAn6bmX3NzK7Au/Y+H+qui4diXVyiiIQiXHbCZSzYvIB9pfv8LqdF4exsXEkJpatX+12KiIjEWGsBfxNwEvB5YKZz7mB0+VTggdiVlXhyx+RSWVPJU+uf8ruUFkWmTIFAQMPWioh0Aq3dRb/bOXeLc+5y59zCestfcc79LvblJY4R3Ucwpd8UHvvoMapqqvwup1nBrl3JOOUUTR8rItIJtHYX/bNHenVUkYli1thZ7CzeyWt5r/ldSovCOdmUvvce1YWFfpciIiIx1FoXfTbe4DOvA78Dft/oJfWcNfgs+kX6xfXNdl1ycqCmhpI342rkXxERaWetBXw/4IfAycAfgPOBvc6515xz8Xua6pOUQApXj76aZTuWsfHQRr/LaVbGqadi4bCuw4uIJLnWrsFXO+fmO+duwLuxbgPwqpl9rUOqS0BXjrqSlEAK89bO87uUZllqKuEzJmlcehGRJNfaGTxmlmZmVwKP4g1KcxcQv7eK+6xXRi8+NfRTPPvxs5RUlvhdTrO65ORQsXkzldsbTw0gIiLJorWb7B4CluA9A3+Hc+4M59xPnXP5HVJdgpo1dhZFlUU8vzEuR9etmz5W3fQiIsmrtTP4zwGjgW8AS8ysIPoqNLOC2JeXmE7tfSrjeoyL2/HpU0eOJKV3b3XTi4gksdauwQecc5nRV9d6r0znXNeOKjLRmBm5Y3PZcHADK3at8LucJsyMSE42xUuX4mo04rCISDJq9Rq8HJsZw2fQNbVr3D4yF8nJofrAAcrXxuUsvCIicpwU8DGSkZLBp0d+mpe3vszukt1+l9NEeKqmjxURSWYK+BiaOWYm1a6aJz56wu9Smgj17UPaqJEatlZEJEkp4GNoSNchTBs4jcc/epzK6srWN+hgkZwcSlaupKa83O9SRESknSngY2zW2FnsLd3Lv7f+2+9SmghnZ+PKyyl9+22/SxERkXamgI+xaQOmMbDLQOasneN3KU1EzjgDQiEOzJ2Hq672uxwREWlHCvgYCwaC5I7J5e3db/PRgY/8LqeBQCRCr1u+TOGCBeR/81vqqhcRSSIK+A5wxagrSAumxeUjc71vu42+P/g+hf/6F9tu/jLVRUV+lyQiIu1AAd8BuqV1Y8bwGTy/8XkKKuJvAMAeN9zAgN/+hpKVK9ly/fVU7d3rd0kiInKcFPAdJHdsLqVVpTy74Vm/S2lWt0svZfCf/0TFxk1svu46KvLy/C5JRESOgwK+g5zU8yTG9xrPvHXzqHHxOTxsl+nTGfLA/VQfPMTmWbMoW7fO75JEROQYKeA7UO7YXDYXbGbZjmV+l9Ki8GmnMezRR7BAkC2f/RwlK+JvLH0REWmdAr4DXTDsAnqk94jLm+3qSxs1imFzZpPSqxdbb/oihS+/7HdJIiJylBTwHSg1mMqVo67ktbzX2F603e9yjig0YABDZ/+dtNGjyfva1zn45FN+lyQiIkdBAd/Brhl9DQCPrXvM50pal5KVxdAHHyAyZQo7fvQj9t17r98liYhIGyngO1j/Lv05a9BZPLX+Kcqr439gmUAkwuC7/0LXi2aw+3e/Z9dvfqs55EVEEoAC3gezxs7iQPkBFm5e6HcpbWKpqQz43e/Iuu469t9/Pzt++CNcZfxNniMiIocp4H0wtf9UhnUdFvc329VngQB9b/8Rvb7+NQ49/TR5X/0aNaWlfpclIiItUMD7wMzIHZvL6r2r+WDvB36X02ZmRu9bb6Xff/+YokWL2PqFm6g+dMjvskREpBkKeJ9cdsJlZKRkxOUsc63Jys1l4P/+L2Xvv8+Wz36Oyl27/C5JREQaUcD7JDM1k0tGXML8zfM5WHbQ73KOWtcLL2Dw3+6hMj+fLbOupXzTJr9LEhGRehTwPsodm0t5dTn/2PAPv0s5JpGpUxny8MPUlJWx5drrKH3vfb9LEhGRKAW8j0Znjeb0vqczb908qmuq/S7nmGScfBLDZv+dQDjM1htuoHjpUr9LEhERFPC+yx2bS35RPovzF/tdyjFLHTaMobNnExo4kG03f5mC+fP9LklEpNNTwPvs3CHn0jujN3PWJd7NdvWF+vZh6KOPkD5+PPnf+g8OzEns/RERSXQKeJ+FAiGuHn01b+S/wdaCrX6Xc1yC3box5L576fLJT7Lzjp+w5//+hHPO77JERDolBXwcuGr0VaRYCvPWzfO7lOMWSE9n0B/votsVV7D3//6PXT/9Ka46Me8vEBFJZAr4ONA73Jtzh57LPzb8g9KqxB8dzlJS6P+Ln9Pjpi9wYPYc8r/9bWoqKvwuS0SkU1HAx4ncMbkUVhTywsYX/C6lXZgZfb/zHfp859sUvjifvFtuobqo2O+yREQ6DQV8nDi97+mM7D6SuevmJtV165433UT/X/yC4jffYuuNN1K1f7/fJYmIdAoK+DhhZswaO4u1+9eyas8qv8tpV92vvIJBf/wj5R99xJZrr6MyP9/vkkREkp4CPo5cMuISuoS6JOT49K3JPOdshtx/H1X797N51rWUr1/vd0kiIklNAR9HwqEwl4+8nIVbFrK3dK/f5bS78OmnM/SRR8A5Nn/2c5S8/Y7fJYmIJC0FfJyZOWYmVTVVPPnRk36XEhPpY0YzdM4cUrp3Z+sXvkDRa6/5XZKISFJSwMeZ4d2GM7X/VB7/6HGqaqr8LicmUgcNZOjsv5M2YgTbbr2NQ88843dJIiJJRwEfh2aNncWukl28uu1Vv0uJmZSePRny8EOEzziD7d/7PvsefNDvkkREkooCPg6dNegs+kf6M3ftXL9Lialgly4MvuevZF5wAbt/9Wt2//5/kuoRQRERPyng41AwEOSaMdfw5s432Xhwo9/lxFQgNZWB//N7uufOZN/f/saO22/HVSXnpQkRkY6kgI9TV466klAglJSPzDVmwSD9fvxjet12G4eefIq8r3+DmrIyv8sSEUloCvg41SO9BxcOu5DnNj5HcWXyD/FqZvT+2lfpe/vtFL3yClu/+EWqCwr8LktEJGEp4ONY7thciiuLee7j5/wupcP0+Ox1DPjdbyldtZotn7ueyt27/S5JRCQhKeDj2Cm9TuHEnicyd21yjU/fmm4XX8zgv/yFim3b2HLdZ6nYutXvkkREEo4CPo6ZGbljcvn40Mes2LXC73I6VJczpzH0wQeoKSxk87XXUfbhh36XJCKSUBTwcW7G8Bl0S+vWKW62ayxj/HiGzv47lhpiy/U3UPzmW36XJCKSMBTwcS49JZ0rR17Jy1tfZmfxTr/L6XBpI0YwbPZsUvr1ZduXvkTBwoV+lyQikhAU8AngmjHXUONqeOKjJ/wuxRehfv0Y9uijpJ94Ivnf/BYHHnvM75JEROKeAj4BDMocxCcGfYInPnqCyupKv8vxRbB7d4bcfx+RM6ex879+zN67/9qpbjwUETlaMQ14M7vQzNaZ2QYz+34z683M7oquX21mE6PLB5vZK2a2xsw+MLNvxLLORJA7Jpd9Zft4aetLfpfim0A4zOA//Ymul13KnjvvZNcvf4mrqfG7LBGRuBSzgDezIPAnYAZwIjDLzE5s1GwGMCr6uhn4S3R5FfD/nHPjgKnAbc1s26lMGziNwZmDO+XNdvVZKMSAX/2KHjfcwIGHH2H7d7+Hq6jwuywRkbgTyzP4ycAG59xG51wFMBe4vFGby4GHnWcZ0N3M+jvndjjn3gZwzhUCa4CBMaw17gUswMwxM3ln9zus27/O73J8ZYEAfb7/PXr/v/+g4Pnn2XbbV6kpKfG7LBGRuBLLgB8IbKv3OY+mId1qGzMbBpwGvNncl5jZzWa2wsxW7Nmz53hrjmufHvlp0oPpnf4sHrwxAnp96Uv0/9lPKX7jDbbceCNVBw74XZaISNyIZcBbM8sa3xV1xDZm1gV4Evimc67Zgcmdc/c45yY55yb17t37mItNBN3SunHRiIt4YdMLHCo/5Hc5caH7VVcx6K4/UL5mLVs++zkqd+zwuyQRkbgQy4DPAwbX+zwI2N7WNmYWwgv3vzvnnophnQkld0wupVWlPLPhGb9LiRuZ553H4Hv/RtWuXWyedS3lH3/sd0kiIr6LZcAvB0aZ2XAzSwVygWcbtXkWuD56N/1U4JBzboeZGXAfsMY59z8xrDHhjOs5jgm9JzBv3TxqnO4grxWZPJmhjzyMq6piy7XXUbpqld8liYj4KmYB75yrAr4KLMC7Se4x59wHZnaLmd0SbfYCsBHYAPwNuDW6fBrwOeAcM3s3+rooVrUmmtyxuWwt3MrS7Uv9LiWupI8bx7DZfyfQrRtbPn8jRa8v9rskERHfWDINFjJp0iS3YkXyT8pSUV3B+U+cz/he4/njuX/0u5y4U7VnD1tv/jLlGzYw4Je/pNslF/tdkohITJjZSufcpObWaSS7BJQaTOUzoz7Da3mvkV+U73c5cSeld2+GPvwQ4QkT2P6d77D/kUf9LklEpMMp4BPUNWOuwcyYt26e36XEpWBmJoPv/Rtdzj2HXT//Obv/8AcNbSsinYoCPkH1i/TjnMHn8I/1/6C8utzvcuJSIC2NQXfeSferr2LfX+5m02c+w567/kjpqlW46mq/yxMRiSkFfALLHZvLwfKDzN803+9S4palpNDvJz+h73/9J4G0dPbefTebZ+ay/sxPkP/d73Lo+X9qgBwRSUq6yS6BOef49DOfJpwSZs4lGt2uLaoOHKD4jSUULXqN4tcXU33gAAQCZJx6Kl3Omk6X6dNJGzcO70lNEZH4dqSb7BTwCW72mtn88q1fMvui2ZzS+xS/y0korrqasg8+oOi1RRQtWkTZe+8B3k16kU98gi7TpxOZlkMwM9PnSkVEmqeAT2JFFUWc+/i5nDf0PH5+5s/9LiehVe3dS9HixRQvWkTR4jeoKSiAlBTCp51Gl7OmE5k+nbRRo3R2LyJxQwGf5H627Gf8Y/0/eOnql8hKz/K7nKTgqqooXbWq7uy+fO1aAFL696fL9Ole4E+ZQiAS8blSEenMFPBJbsOBDVzx7BV8c+I3uemUm/wuJylV7tpF8euvU/TaIoqXLKGmuBgLhQifcUbd2X3qsGE6uxeRDqWA7wS+sOAL5Bfm88KVLxAMBP0uJ6m5igpK3n6HokWLKFr0GhUbvMltQkOG0OUTn6DLWdMJT55MID3d50pFJNkp4DuBhZsX8v9e+3/88Zw/8snBn/S7nE6lIi+f4tcXeWf3y5bhysqwtDTCU6dEu/PPInXQIL/LFJEkpIDvBCprKrnwiQsZmTWSv57/V7/L6bRqysspeWs5Ra8vovi1RVRs2QJA6ogRddfuM04/nUBqqs+VikgyUMB3Enevups/vfsnnvv0cwzrNszvcgSo2LyZokWvU7RoESVvvYWrqCAQDhPOyabLJ6bTZfonCPXv73eZIpKgFPCdxN7SvZz/xPnkjsnle5O/53c50khNSQnFb75J0SLv7L5y+3YA0kaPrhtkJ2PCBCwU8rlSEUkUCvhO5LuvfZfF+Yt56eqXCIfCfpcjLXDOUfHxx95jeK+/TsmKFVBVRSAzk8i0aV53/ifOJKV3b79LFZE4dqSAT+noYiS2csfm8uLmF3lh0wtcNfoqv8uRFpgZaSNHkjZyJD1v+gLVRUUUL13qDbLz2iIK53vzC6SfdFLd2X36KadgQT0hISJtozP4JOOc4+rnrsbheOLSJ/RcdgJyzlG+bl3dIDul77wDNTUEu3ePDqH7CSJnnklKlgY1EunsdAbfiZgZuWNzuWPpHbyz+x0m9p3od0lylMyM9LFjSR87ll5fvpnqQ4cofuMNL/AXL6bguefAjIzx44mcNZ0u088i/cRxWECTQ4rIYTqDT0IllSWc9/h5nDnwTH5z1m/8LkfakaupoeyDDyla9Jo3Qc7q98A5gr161Q2yEznzTIJduvhdqoh0AJ3BdzLhUJjLR17O3HVz2Vu6l14ZvfwuSdqJBQJknHIyGaecTO/bbqNq/36KFy/2rtu//DKH/vEPLDWVyJln0vXCC+hyzjkKe5FOSmfwLSk9AKmZEEzMY6AtBVu45B+XcNuE27jl1Fv8Lkc6gKuqonT1agoXLKRgwQKqdu7EQqGGYa+pb0WSih6TO1oVJXDPJ2HgRLj8z5Cg1zZv+dctrD+wnvlXzScU0LPVnYmrqaF01SoK5y+gYOFCqnbsUNiLJKEjBXxiJlespYbhlKtg1RyY/z1I0IOg3LG57C7dzStbX/G7FOlgFggQPu00+v7g+4z890sMmzuHrOuuo2ztWrZ/7/usz5nGtlu+wsGnn6a6oMDvckUkBnQG3xLnYOHtsPT/YPp34Jzb2+f3dqDqmmoueuoiBmYO5P4L7ve7HIkDzjnKVq+m4MX5FCxcQNX2HRAK0SUnh8wLLyTz3HMIdu3qd5ki0ka6ye5YmMGnfgblBbDot5DWFaZ93e+qjkowEOSaMddw59t3suHABkZmjfS7JPGZmZFx6qlknHoqfb73XS/s5y+gYMF8il57jR0Ke5GkoTP41tRUw5M3wQf/gEvuhEk3tu/vj7EDZQc47/HzuGLUFdw+NfF6IaRjOOcoe+89Cl6cT+GCBd44+aEQkZxsul4QDftu3fwuU0Qa0U12x6uqAuZeCxtegs/c612fTyA/WvwjXtryEv+++t90SdUjU3JkdWE/fwGF8+cfDvvsqV7Yn3euwl4kTijg20NFCfz9Ktj2JuTOhtEXxOZ7YuD9ve8z65+z+MHkH3DtuGv9LkcSiHOOsvffp2D+fArnL6AyPx9SUohkZ9P1wgvIPPdcgt27+12mSKelgG8vZQXw0KWwZy1c9wQM/0TsvqudzXp+FiVVJTx9+dMan16OiRf2H1Aw/0WFvUicUMC3p+J98OBFcCgPrn8WBp0e2+9rJ89seIbb37idez91L1P6T/G7HElwtWFfuGA+BfMXUJmX54X91Knec/bnnqvJcEQ6gAK+vRXsgPsv8O6w//wL0PfE2H/ncSqvLue8x8/j9L6nc+fZd/pdjiQR5xxlH3xI4fwXG4b9lClkXngBmeedp7AXiREFfCzs3wQPzABXA1+YDz1GdMz3Hof/Xfm/PPjBgyz4zAL6Rfr5XY4kobqwrz2z37YNgkEiU6cq7EViQAEfK7vXeiGf1gVunA/dBnbcdx+D/KJ8Zjw5gy+e8kW+PjGxnumXxOOco+zDD73hchcsoHLrVi/sa8/szz9fYS9ynBTwsZT/Njx0GXTtDze+CJH4nrnta//+Gqv3ruZfV/2L1GCq3+VIJ+Gco3zNGm8EvQZhP5nMCy4k8/zzSOnRw+8yRRKOAj7WNr8Bj14JvcfADc9Bevw+I/xG/hvc8tIt/OoTv+LiERf7XY50QnVhHx1Br3KLwl7kWCngO8JHC2HuLBh0Bnz2KW/CmjhU42q47OnL6J7WnUcvetTvcqSTc85RvnZt9Mz+cNiHJ5/hDarzqfMV9iJHoIDvKO8/CU/cBCPPhdw5kBKfXeCPfPgIv1n+Gx675DHG9RzndzkiQL2wj46gV7FlCwQChCdP9p6zP/98Unr29LtMkbiigO9IKx+C574OJ34arrofAkF/62lGQUUB5z1+HjOGz+COnDv8LkekCecc5evW1Y2gV7F5s8JepBkK+I625P9g4Y/gtM/CpX+EQMDvipr47yX/zT83/pOXrn6Jbmnxe8+AiHOO8o8+ouDFFxuEfeqwYaSOGE7a8BGkjhhB2ojhpI4YQTAz0++SRTqMpovtaDlfhbJDsOg33jSzF/zCm342jswaO4sn1z/J0xue5oaTbvC7HJEWmRnpY8aQPmYMvb/xDco/+ojCf71E+bq1lG/cRNGrr0FVVV37lN69SR0xokn4p/Trh8XhwbZIrCjgY+XsH3oj3S37s3dX/Se/73dFDYzpMYaJfSYyb908Pnfi5wiY/uGT+Fc/7Gu5ykoq8vKo2LSJio0bKf94IxUbN1Lw/D+pKSw8vG1GBqnDh5E24gQv/EeMIHX4CFKHDSWQlubH7ojElAI+Vszggl9CeSG8+kvvTD77Vr+raiB3bC7fXfRdXt76MucNPc/vckSOiYVCpA0fTtrw4XDOOXXLnXNU79tH+caNVGzcRMWmjZRv3ETpO+9Q8Pzz9X6BERo0qN4Z/3DSTjiB1BEjNBCPJDQFfCwFAnDpXV7IL/gBpGXCxM/5XVWd84acx8AuA/nWq99iSv8pXDv2Ws4adBbBOLwxUORomRkpvXqR0qsXkcmTG6yrKS2lYssWyj/+uEH4lyx7E1deXtcu2L17o+5+78w/NGgQFtT/JxLfdJNdR6gqhzm5sPFV7876k67wu6I6B8oO8OT6J5m3bh47i3cyIDKAa8Zcw2dGfYbu6d39Lk+kQ7maGiq376Bik9fNX75xExUff0z5pk1U79tX185CoehNfiMadPenDR9GIBLxcQ+ks9Fd9PGgohgeuRLyV8KsuTAqvrrEq2qqeHXbq8xZO4e3dr5FWjCNGcNnMGvsLE7sGf+z5YnEWvXBg5Rv2tTgjL9i40Yqtm2D6uq6din9+kXv6G94rT+lT28szm62lcSngI8XZYfgwUtg73r43FMwNMfvipq1/sB65q6dy3Mbn6O0qpQJvScwa+wszh96PqFgyO/yROKKq6igYuvWJtf6KzZupKa4uK5dIBI5/Dhf/Wv9gwdjqfE5KJbEPwV8PCna481AV7gTPv8cDDjN74paVFBRwDMbnmHu2rlsLdxKr4xeXD36aq4afRV9wn38Lk8krjnnqNq9Jxr40fDf6L2v2rnzcMNgkNTBg5uG/4gRBLtpjAo5MgV8vDmUD/dfCJXF3gx0vce0vo2PalwNS7YvYfaa2SzOX0zQgpw/9HxmjZvFhN4T1O0ocpSqi4qp2Ly5SfhXbN6Mq6ysaxfs2ZOUvn1IyepBMCuLYFYWKT2you97EMzqTkqP6Lru3XXjXyekgI9H+z72Qj4QhC/Mh6xhflfUJlsLtjJ33VyeXv80hZWFjO0xlmvHXsuM4TNIT0n3uzyRhOaqq6nMz/dC/+ONVGzeRNWevVQd2E/1gYNU799PTVFR8xubEezalWBt4PfIIiUri2D3rOiyegcDWT1IyeqOhcM6QE9wCvh4tesDeOAiyOgON8735pRPECWVJfxz0z+ZvWY2Gw5uoFtaN64ceSUzx85kYJeBfpcnkrRcRQVVBw5SffAA1fv3U33gAFUHDlC9/0D0/eGDgdp19Uf6q8/S0g4fDNQdCGQdPhjoXu9AoUcPgt26YSl6ujqeKODjWd4KeOgy6D4EbnwBwok1NaZzjhW7VjBn7Rxe3voyNa6GswafxbVjr2Vq/6k6OxDxmXOOmqKiw4EfPRCoPrC/xQODVnsJshodDGQ16jWIrkvJylIvQYwp4OPdpkXw6FXQ9yS44VlvQJwEtLN4J4+te4wn1z/J/rL9DO82nNwxuVw+8nIiIT0bLJIoXEUFVQcPRg8EvJ6CY+4lSE2tF/jdo/cO1DsYyOpBsEcWqYMHk9K3r+YLOEoK+ESw7kWYex0MyYbPPgGhDL8rOmYV1RUs2LyA2Wtm8/6+94mEIlx2wmXkjs1lRLcRfpcnIu2s1V6CRgcD1QcONJgnoFbdfAH1niRIHTGC1KFDCaTrHp/mKOATxerH4akvwahPQe7fIQmeOX9vz3vMWTuH+ZvnU1lTSXb/bGaNncX0QdM1JK5IJ9agl2DfPiq2bqVi06a6MQQqt2+H2nwyIzRwIKnDhzd9lLBnz059CUABn0hW3A/PfwtO/gxc+TfvLvsksK90X92QuLtLdjOwy0BmjpnJFSOv0JC4ItJE7XwBdUMGb9pE+aaNVGzajCstrWsX6NqVtOHDSR0+/PBYAiNGeAMIhRL/JKk1CvhEs/hOeOnHMPEGuPQPcTeX/PGoqqnilW2vMHvNbFbsWkFaMI2LR1zMrLGzGNtjrN/liUicczU1VO3cWXemX76pdgTBTVTt3n24YUrK4QGEhg9L2gGEFPCJ6KU7YPH/QM7X4PyfJlXI1/rowEfMWTuHf278J6VVpZzW5zSuHXst5w49l1Ag+Y+8RaR9VRcVUbFp0+Gz/o3eWALlm7dAowGE0qJn/IfnCxhOaMCAhBssSAGfiJyDF74Dy/8G59wO07/jd0Uxc6j8kDck7rq5bCvcRu+M3lw9+mquHnM1vTJ6+V2eiCQ4V1V1eAChjZsanPVXHzhQ185SU+tmCWxwrX9Y/M4S6FvAm9mFwB+AIHCvc+5XjdZbdP1FQAnweefc29F19wOXALudcye35fuSKuABamrg6Vtg9TyY8RuY8mW/K4qpGlfD4vzFzFk7h8X5i0kJpHD+0PO5duy1nNr71E59I42IxEbVgQNNzvrLN22kclue929wVN0sgY3u8E/p08fXf5t8CXgzCwIfAecDecByYJZz7sN6bS4CvoYX8FOAPzjnpkTXTQeKgIc7bcADVFfBY9fDun/Cp/8CE671u6IOsaVgC3PXzuXpDU9TVFnEuB7jmDV2lobEFZEOUVNRQWVbZgkMhxve4Bc9AEgdOpRAWlrM6/Qr4LOB/3bOXRD9/AMA59wv67X5K/Cqc25O9PM64JPOuR3Rz8OA5zt1wANUlsHsa2Dz63D1Q3DiZX5X1GFKKkt4fuPzzFk7hw0HN9A9rTtXjrqSmWNmMqDLAL/LE5FOxjlH1Z49TUK/fNNGqrbvONwwECA0aFCT5/rTRo1q15v8jhTwsRxUeCCwrd7nPLyz9NbaDAR2IIeF0iF3NjzyaXjyJkibByec43dVHSIcCnPNmGu4evTVLN+5nDlr5/DgBw/y4AcP8slBn2TWuFlM6TdF3fci0iHMjFCfPoT69CEytWGk1ZSUULFlS5Oz/pJlb+LKywHo893v0vMLN3ZIrbEM+Ob+xW3cXdCWNkf+ErObgZsBhgwZcjSbJpa0LnDd4/DgJd6Id5/7BwyZ6ndVHcbMmNx/MpP7T2ZH0Q4e++gxnvzoSV7e9jIjuo1g1thZXHrCpRoSV0R8EwiHSR83jvRx4xosdzU1VG7fQcWmTaQOG9ph9aiLPtEU7fammS3eC59/HvqP97si35RXlzN/03xmr53Nh/s+pEuoC5ePvJyZY2YyvNtwv8sTEYm5I3XRx3JU/+XAKDMbbmapQC7wbKM2zwLXm2cqcKg23KUFXfrA9c94E9I8cgXsXe93Rb5JC6Zx+cjLmXvxXB696FHOGnwW89bN47KnL+PL//oyr217jeqaar/LFBHxRawfk7sIuBPvMbn7nXM/N7NbAJxzd0cfk/s/4EK8x+RudM6tiG47B/gk0AvYBfzYOXffkb6vU5zB19q7Ae6/AFLS4Avzvelmhb2le3nyoyd5bN1j7C71hsTNHZPLFaOuoFta8oxeJSICGugmee1Y7V2Tj/SELyzwzu4FgMqaSl7e+jJz1s5h5a6VpAfTGd97PIMzBzM4czBDug5hSOYQBmcOJhwK+12uiMgxUcAns61venfX9xjhXZPPyPK7orizbv86Hv/ocdbuX8u2wm3sL9vfYH3P9J4M6eqFfW3o137WWb+IxDMFfLL7+GWYPRP6nwqfe9q7415aVFRRxLbCbWwt3Mq2wm3e+4KtbC3cyu6S3Q3adkvrxuAugxnc1Qv/2uAfnDmYnumde5pKEfGfAr4zWPMcPHYDDJsG1z7uPTsvR620qpT8wvy68N9asLXuYGBH8Q5q3OGhK8Mp4QaBX/8AoE+4DwGL5T2sIiIK+M7j3Tne2PVjLoZrHoKgZmRrT5XVleQX5Tc4+689AMgryqOqpqqubWog1Qv+rvXCP3MIg7sOpn+kPymBWA5BISKdhV8j2UlHmzALygvhxe/AM7fBp++GgM4i20soGGJYt2EM6zasybrqmmp2luysC/z63f7Lti+jrLqsrm2KpTCgy4C6bv+6a/9dBzOoyyBSg6kduFcikqwU8Mlmys1Qfghe/hmkdoGLf5+Uc8nHm2AgyMAuAxnYZSDZZDdY55xjT+mehuFfuJWtBVtZtXsVRZVFdW0No1+kX13g1575D8ocpDv+ReSoKOCT0Se+DWUFsOQuSO8G5/3Y74o6NTOjT7gPfcJ9mNSvYU+ac46D5QfrAr/+AcC/t/ybA+UHGrTvndG76aN+0csAXVO7duRuiUicU8AnIzM4/yded/3i/4H0rnDmt/yuSpphZmSlZ5GVnsWpvU9tsr6wovDwNf+Cw9f+l25fyjMfP9Ogbfe07gzJHEL/Lv3pE+5D33Bf+kb60i/cj77hvvQK9yIU0H0ZIp2FAj5ZmXnd8+WF8NJ/e0PbnvFFv6uSo5SZmsmJPU/kxJ4nNllXUllCXlGed9YfDf+thVtZt38di/IWUVpV2qC9YfTM6OkFfzT8aw8E+kW8g4A+4T6kp+gJDJFkoIBPZoEgXHE3VBTBP78NqZlw6ky/q5J2Eg6FGZ01mtFZo5usc85RUFHA7pLd7CrZxa7iXd7P6Gtr4VaW71pOYUVhk227pXVrchBQ2wvQN+It75KqsRZE4p0CPtkFQ3D1g/D3q+Hpr3hn8mMv8rsqiTEzo1taN7qldWNU1qgW25VUlhw+CGh8IFC8iw/2fdBk5D/wxgCoDfv6wV//ffe07hoISMRHeg6+sygvhIcug10fwHWPwYhP+l2RJIiK6gp2l+xusTdgV/Eu9pTuaTAIEHhjAfQJ9zkc/vUPAqKfe6b3JBgI+rRnIolPA92Ip2Q/PHgxHNjiTTk7+Ay/K5IkUVVTxb7SfU16A3aW7PSWRQ8KKmsqG2wXtCC9Mno1CP9+kX4NbhLsk9GHkAZtEmmWAl4OK9wJ918Ipfvh8y9Av5P9rkg6CeccB8oPNAj8ncU72VWyq+7AYGfxziY3BwL0SO/RoBegX6QfAyID6h4Z7JbWTZcDpFNSwEtDB7Z4IV9T5c0l3/MEvysSAbyDgKLKInYV1wv9kp0NPu8q2cWh8kMNtssMZdYNBlT/NaTrEM0LIElNAS9N7VkHD8yAUBiufsibiS6oey4lMdROClQ7MNC2wm1sK/IeF9xetJ0q13BegIGZA5uE/6DMQRoaWBKeAl6at/1deOhSKC+AlAwv5AdNgoGne6/uQzTMrSScqpoqdhbvrAv+vMK8BgcCJVUldW0No2+kb5Pwr31lpmb6uCcirVPAS8sKd8Hm1yF/JeStgB2roLrcWxfpHQ37STDodBgwETK6+1quyPFwzrG/bH/DM/96r8aPBHZP6153tt84/Htn9NZ1f/GdAl7arroSdr0fDfyVkL8C9n50eH3PUV7o157p9z0ZUtTFKcmhuLK4yRn/1sKt5BXmsaN4R4NHAdOD6c0G/+DMwfTv0l/DAkuHUMDL8Sk7BPlve6Ffe6ZfvNtbF0yD/uMbnulnDVfXviSdyupKthdvb3LWn1eYR15hXoMpgYMWpF+kX4td/5oVUNqLAl7al3NwKM87u68909/xLlRGr21m9Dh8Hb/2TD/cw9eSRWKpxtWwt3RvgxkB63oCirY1ueu/Z3rPJjf81b7vkd5DXf/SZgp4ib3qKtizxju7rz3T370GiP73lTU8GvbRwO93CoQ0qYl0DgUVBU2Cv/ZgYHfJbhyH/x0Op4QbhH/39O5kpGQQTgkTDoUJp4S9z9H34ZD3OSMlQ48DdkIKePFHeaF3p379M/3C7d66QMgbZGfgpMNn+T1OgID+gZLOpby6vOkjf9FXflF+k9H/jqQ26BscALRwQND4fUao6bJwKEx6MF09CnFMAS/xo2D74ev4+Sth+zvebHcA6d28O/Xrn+l36e1vvSI+cs5RVl1GSWUJpVWllFSVUFJZQklV9HNlo5/NrC+pin6uLK1b3txogS0x7OgOEJpp09xBRigQ0oFDOzhSwGtkE+lYXQd4r3GXep9rqr1Bd/JXHj7Tf/1/wFV767sPqXcD3yToNx5SdYOSdA5mVndG3p5qXA1lVWUtHxC0csBQWllKYUUhu0p2NVhfXvuIbRsELdjgwCA9JZ3UYCppwTTvZyCt7n2D5fV+HtOyQGqnObDQGbzEn4pi73n8ujP9t+HQVm+dBaHvSfVu4JsEvUara18kDlTVVLWpR6G5NuXV5VRUVzT42dyyipqK464zNdD2A4bjOpjogAMMddFL4ivcBdvfjgb+Csh/B2rvTE7NhIGnNTzTz+znb70iEhM1robKmspWDwSOZ1lL6ypqKo6ql6I535n0Ha4/6fp2+tNQF70kg8y+MGaG9wKoqYF9Gw537eetgCV/9CbQAeg6sOFjev0nQFoX38oXkfYRsEDdWbEfnHN1BxjHchAxse/EDqtVAS+JKRCA3qO914RZ3rLKMti5ul7X/gpY86y3zgJe6Gf2h679672P3hOQ2d976dE9ETkCM6u7LyCT+J6rQAEvySOUDoMne69axfsOP5d/cIt3F/+uD2HDvw/fvV9fuCdkDogeBAxo/n16d43UJyJxTwEvyS3SE0Z/yns1VlbgBX7hdijY0eh9vvcIX/GeptuFwk3P/uveRw8CuvSFQDD2+yci0gIFvHRe6V29V5+xLbepqoDCHd6rYHv0IKDe+61LvQOCxoORWNAL+SY9AfUuDWT21yN/IhIzCniRI0lJhayh3qslNTVQsi969t/MQcDe9bDxNSgvaLptevcj9wR0HQgZWbokICJHTQEvcrwCAW/EvS69of+pLbcrL4oGf7531l93QBB9v/M9KNoNNHp0NSXde+zvSDcJdukLQU1PKiKHKeBFOkpaF0gbBb1GtdymuhKKdjXqCcg/fI9A/gpYswOaPItr0KWPF/rhnpCW6b3Sux1+n5YJaV0P/0zvenh5KKLBgkSSjAJeJJ4EQ9BtkPdqiXNQsr/hDYH1LwmUHYRD27zJfsoKoLK4DV9s9cI/s2H41x0YdG1mXdeG26VGdDlBJE4o4EUSjZn3dECkpzftbmtqqr2wLy/07gOo/76soJl10Z8l++HAlsPLKkvaUFugaU9Bk4OGltbV71EI60BB5Dgp4EWSXSAIGd291/GoroKKwsM9Aw0ODOp9Lmt8oLAXDmw6vL5NBwrB5i8lND4wSI14BwOpXbz3qZFG76OvYKoOGKTTUcCLSNsEU7w7+jOyju/3VFe20KNQCGWHWu5tKNoN+z4+/LmqrO3fGUhpGv6hRgcBdeuOcMDQeBuNdSBxTAEvIh0rGIJwD+91PKqrvPsLKoqhosQbmbCi9nP0fWULy2u3KdpV73N0fe1UxW2RktG0t6C5XoTQEdY1Xp6Srt4GaRcKeBFJTMEUCHbznhRoL85BdUUzBwNFRz6QqCiud7BRDMV7G27TphsdoyzghX0ofPgSRCgj+gp7QzLXX5ZSf11Go1fYO2CoW1dve/U+JD0FvIhILTNISfNex9vDUF9NTbQ34QgHBc0eSBR5kyhVlkBlKZQe8C5NVJYeXtaWexqaE0w9hgOERuua3b7RAYXGZ/CNAl5EJNYCgeg4CF2Avu37u52DqvJ6gV8KVaWNDgIafa4qa3iAUP8gorLUm4OhyfYl4GqOvj4LtnKAED0gSEn1DjqCqd5BQTCt3vvUI79PaWV93fs0r+eik1wCUcCLiCQys+hZc4ynOnbOu0Gy7iCh8YFDo4OGqrJmDjAaHTSUHqjXvsK7PFJdGf1ZcXT3Q7SZtXIQEPJ6cI60/lgORGoPQnqcAN0GxmC/mlLAi4hI68y8kEpJPf5HLtuqprpe4NcL/gavypbfV5U3Wn6EttXNtK0s9Z7sqN/2eA9ELvwVTP1K7P7M6lHAi4hIfAoEvVeseyeO19EciGQN67CyFPAiIiLHI04PRDS7hIiISBJSwIuIiCQhBbyIiEgSUsCLiIgkIQW8iIhIElLAi4iIJCEFvIiISBJSwIuIiCQhBbyIiEgSUsCLiIgkIQW8iIhIEoppwJvZhWa2zsw2mNn3m1lvZnZXdP1qM5vY1m1FRESkZTELeDMLAn8CZgAnArPM7MRGzWYAo6Kvm4G/HMW2IiIi0oJYnsFPBjY45zY65yqAucDljdpcDjzsPMuA7mbWv43bioiISAtiOV3sQGBbvc95wJQ2tBnYxm1j6o7nPuDD7QUd+ZUiIpLkThzQlR9felKHfFcsz+CtmWWujW3asq33C8xuNrMVZrZiz549R1miiIhIcorlGXweMLje50HA9ja2SW3DtgA45+4B7gGYNGlSswcBx6KjjrBERERiIZZn8MuBUWY23MxSgVzg2UZtngWuj95NPxU45Jzb0cZtRUREpAUxO4N3zlWZ2VeBBUAQuN8594GZ3RJdfzfwAnARsAEoAW480raxqlVERCTZmHPt1qvtu0mTJrkVK1b4XYaIiEiHMLOVzrlJza3TSHYiIiJJSAEvIiKShBTwIiIiSUgBLyIikoQU8CIiIklIAS8iIpKEFPAiIiJJSAEvIiKShBTwIiIiSUgBLyIikoQU8CIiIklIAS8iIpKEFPAiIiJJSAEvIiKShBTwIiIiSSip5oM3sz3Alnb8lb2Ave34++KV9jO5aD+Ti/YzubT3fg51zvVubkVSBXx7M7MVzrlJftcRa9rP5KL9TC7az+TSkfupLnoREZEkpIAXERFJQgr4I7vH7wI6iPYzuWg/k4v2M7l02H7qGryIiEgS0hm8iIhIElLAA2Y22MxeMbM1ZvaBmX0juryHmf3LzNZHf2b5XevxMLN0M3vLzFZF9/OO6PKk2s9aZhY0s3fM7Pno56TbTzPbbGbvmdm7ZrYiuiwZ97O7mT1hZmuj/59mJ9t+mtmY6N9j7avAzL6ZbPsJYGbfiv4b9L6ZzYn+25SM+/mN6D5+YGbfjC7rsP1UwHuqgP/nnBsHTAVuM7MTge8D/3bOjQL+Hf2cyMqBc5xzpwITgAvNbCrJt5+1vgGsqfc5WffzbOfchHqP3iTjfv4BmO+cGwucivf3mlT76ZxbF/17nACcDpQA/yDJ9tPMBgJfByY5504GgkAuybefJwNfAibj/Td7iZmNoiP30zmnV6MX8AxwPrAO6B9d1h9Y53dt7biPYeBtYEoy7icwKPo/zznA89Flybifm4FejZYl1X4CXYFNRO8ZStb9bLRvnwLeSMb9BAYC24AeQArwfHR/k20/rwburff5P4HvduR+6gy+ETMbBpwGvAn0dc7tAIj+7ONjae0i2m39LrAb+JdzLin3E7gT73+mmnrLknE/HbDQzFaa2c3RZcm2nyOAPcAD0Usu95pZhOTbz/pygTnR90m1n865fOB3wFZgB3DIObeQJNtP4H1gupn1NLMwcBEwmA7cTwV8PWbWBXgS+KZzrsDvemLBOVftvC7AQcDkaDdSUjGzS4DdzrmVftfSAaY55yYCM/AuLU33u6AYSAEmAn9xzp0GFJPg3bdHYmapwGXA437XEgvRa86XA8OBAUDEzD7rb1Xtzzm3Bvg18C9gPrAK73Jwh1HAR5lZCC/c/+6ceyq6eJeZ9Y+u74931psUnHMHgVeBC0m+/ZwGXGZmm4G5wDlm9ijJt58457ZHf+7Gu147meTbzzwgL9rbBPAEXuAn237WmgG87ZzbFf2cbPt5HrDJObfHOVcJPAXkkHz7iXPuPufcROfcdGA/sJ4O3E8FPGBmBtwHrHHO/U+9Vc8CN0Tf34B3bT5hmVlvM+sefZ+B9z/aWpJsP51zP3DODXLODcPr6nzZOfdZkmw/zSxiZpm17/GuY75Pku2nc24nsM3MxkQXnQt8SJLtZz2zONw9D8m3n1uBqWYWjv7bey7eTZPJtp+YWZ/ozyHAlXh/rx22nxroBjCzM4HXgfc4fM32h3jX4R8DhuD9R3m1c26/L0W2AzMbDzyEd9dqAHjMOfcTM+tJEu1nfWb2SeDbzrlLkm0/zWwE3lk7eN3Ys51zP0+2/QQwswnAvUAqsBG4keh/wyTXfobxbkAb4Zw7FF2WjH+fdwAz8bqs3wG+CHQh+fbzdaAnUAn8h3Pu3x3596mAFxERSULqohcREUlCCngREZEkpIAXERFJQgp4ERGRJKSAFxERSUIKeJHjZGbOzH5f7/O3zey/2+l3P2hmV7XH72rle66OztL2Sqy/KxGZ2aejE1CJJAwFvMjxKweuNLNefhdSn5kFj6L5TcCtzrmzY1VPezGzFB++9tOAAl4SigJe5PhVAfcA32q8ovEZuJkVRX9+0sxeM7PHzOwjM/uVmV1nZm+ZN7/7CfV+zXlm9nq03SXR7YNm9lszW25mq83sy/V+7ytmNhtv4KbG9cyK/v73zezX0WX/BZwJ3G1mv23Uvk11RkdJfDJaz3IzmxZdPtnMlkQniVlSOxqdmZ0U/R3vRusfZWbDzOz9et9d1xNiZq+a2S/M7DXgG2Z2erSulWa2oN7Qn6+a2f+a2aJoj8QZZvaUeXNv/6ze7/5sve//a+3BkJkVmdnPzWyVmS0zs75mloM3Nvxvo+1PMLOvm9mH0drntuU/EpEO5/eUenrplegvoAhvStPNQDfg28B/R9c9CFxVv2305yeBg3jTRaYB+cAd0XXfAO6st/18vIPxUXjjsqcDNwO3R9ukASvwJu/4JN5kLMObqXMA3shZvfFGvnsZ+HR03at483M33qatdc4Gzoy+H4I37DPRP5eU6PvzgCej7/8IXBd9nwpkAMOA9+t9d/0/x1eBP0ffh4AlQO/o55nA/fXa/bpefdvr1Z6HN6rYOOA5IBRt92fg+uh7B1waff+ben/Gjf8etwNp0ffd/f5vUC+9mnv50dUlknSccwVm9jDwdaC0jZstd9FpI83sY2BhdPl7QP2u8secczXAejPbCIzFG3d+fL3egW54BwAVwFvOuU3NfN8ZwKvOuT3R7/w7MB14uh3qPA840cxqt+lq3jj53YCHzGwUXniGouuXAj8ys0HAU8659fW2bcm86M8xwMnAv6LbBPGmHa31bL36PqhX+0a86TrPBE4Hlke3z+DwhB8VePOTA6wEzm+hltXA383saVr/8xPxhQJepP3cCbwNPFBvWRXRS2HmpUlqvXXl9d7X1PtcQ8P/NxuPJ+0AA77mnFtQf4V5Y+8Xt1BfqwnagrbUGQCynXMNDm7M7I/AK865K8xsGN4ZNs652Wb2JnAxsMDMvgh8RMPLhumN6qjdL8ML7uxW6q1fa/16DXjIOfeDZratdM7V/nlX0/K/kRfjHRxdBvynmZ3knOvQqUBFWqNr8CLtxHkTRjyGd8Narc14Z4vgzYEd4uhdbWaB6PXuEcA6YAHwFfOmOcbMRps3o9yRvAmcZWa9otecZwGvHUM9zVkIfLX2g3mTw4B3Bp8fff/5eutHABudc3fhnXGPB3YBfcysp5mlAZe08F3rgN5mlh39XSEzO+koav03cJUdnumrh5kNbWWbQqB25r4AMNg59wrwXaA73kQpInFFAS/Svn4P1L+b/m94ofoWMIWWz66PZB1eEL8I3OKcK8ObWe1D4O3ojWl/pZUeuWhX9Q+AV4BVeHOOt9dUlV8HJkVvOvsQuCW6/DfAL83sDbyu9FozgffN7F28Sw4PO29u8J/gHYg8jzeVcXP7UQFcBfzazFYB7+LNJ94mzrkPgduBhWa2GvgX3nX6I5kLfMfM3sG7FPKomb2HNxPa/zrnDrb1+0U6imaTExERSUI6gxcREUlCCngREZEkpIAXERFJQgp4ERGRJKSAFxERSUIKeBERkSSkgBcREUlCCngREZEk9P8Bn+tO82uoBLkAAAAASUVORK5CYII=\n",
"text/plain": [
"