{"cells":[{"cell_type":"code","execution_count":1,"metadata":{},"outputs":[],"source":["import numpy as np\n","import matplotlib.pyplot as plt"]},{"cell_type":"markdown","metadata":{},"source":["# Intro to Dual Channel Polarimetry"]},{"cell_type":"markdown","metadata":{},"source":["_Written by Manxuan (Rebecca) Zhang_\n","\n","Similar to the single rotating retarder Stokes polarimeter example, we can also use a dual channel polarimeter instead to measure the polarization state of an incoming light \n","source. The components are as follows:\n","\n","- a normalized Stokes vector with some polarization state and no circular polarization\n","- a half wave-plate (HWP) that can rotate\n","- a Wollaston prism to split two orthogonal polarization states\n","- a detector to measure the power of both the ordinary ($I_o$) and extraordindary ($I_e$) beams\n","\n","Using this set-up, we can either use single differencing or double differencing to retrieve the Q and U of $S_{in}$. Note that we are unable to retrieve V as the HWP is not able to rotate circular polarization states."]},{"cell_type":"markdown","metadata":{},"source":["# Single Differencing"]},{"cell_type":"markdown","metadata":{},"source":["As with the the single channel polarimeter, we rotate the half wave plate in steps of $\\delta\\theta$ to sample possible power responses on the detector. A two measurements is necessary completely determine input Stokes Q and U.\n","\n","The simplest way to extract input Stokes vectors is to calculate the single difference, $I_o - I_e$, for each HWP angle $\\theta$, and use that as the observable in retrieving the input Stokes vector.\n","\n","From multiplying the input generalized Stokes vector with the Mueller matrices of the components, we see that the power of the single differences follows a sinusoid of two frequencies:\n","\n","$$P \\propto Qcos(4\\theta) + Usin(4\\theta)$$"]},{"cell_type":"code","execution_count":2,"metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["Input normalized Stokes vector\n","[1.         0.32205895 0.61911011 0.        ]\n"]}],"source":["# We begin with a simulated dataset and a randomized normalized Stokes vector with no Stokes V\n","from katsu.polarimetry import dual_channel_polarimeter\n","from katsu.mueller import wollaston, linear_retarder\n","\n","thetas = np.linspace(0, np.pi, 10)\n","# For testing normalized Stokes vectors\n","I = 1.0  # Total intensity\n","Q = np.random.random()\n","U = np.random.random()\n","V = 0.0  # Circular polarization\n","\n","S_to_measure = np.array([I, Q, U, V])\n","print(\"Input normalized Stokes vector\")\n","print(S_to_measure)"]},{"cell_type":"markdown","metadata":{},"source":["Let's test how well katsu.dual_channel_polarimeter performs at reconstructing the input Stokes vector"]},{"cell_type":"code","execution_count":3,"metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["Stokes Vector Measured\n","[1.         0.32205895 0.6191101  0.        ]\n","Percent Difference (Q, U only)\n","[2.35802600e-06 1.15021692e-06]\n"]}],"source":["S_out = dual_channel_polarimeter(thetas, S_in = S_to_measure, \n","    sub_method = \"single_difference\", normalized = False)\n","print('Stokes Vector Measured')\n","print(S_out)\n","print('Percent Difference (Q, U only)')\n","print(100 * (S_to_measure[1 : -1] - S_out[1 : -1]) / S_to_measure[1 : -1])"]},{"cell_type":"markdown","metadata":{},"source":["Great, it seems like we can retrieve Stokes Q and U up to a certain precision. We now simulate a real dataset, give as input the measured o and e beam powers, and see how well we can reconstruct the Stokes parameters."]},{"cell_type":"code","execution_count":4,"metadata":{},"outputs":[{"data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAYkAAAEXCAYAAABYsbiOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9JElEQVR4nO2deZycVZnvv093V6Urne5Ob+kQCAkK2qJxxYCDItroIDMDTIKCc0WZCXNxDKjjuDDXi5cBr+K43OGqMzKiFxRnAAUhYhQwiAsGA0Q2IUpkaBtC7031kl6rn/vH+3bnTaWqu7q2c076fD+f+tS7nDrv7/eeqveps7zvEVXF4/F4PJ5MVJgW4PF4PB578UHC4/F4PFnxQcLj8Xg8WfFBwuPxeDxZ8UHC4/F4PFnxQcLj8Xg8WfFBwjOHiFwuIntLkO96EVEReWOx857nmId4EZF3isgfRCQlIteF204VkcdFZEpE7i2XPpsJz4mKyFGmtRSTw9VXqfFBYokgIgkRuVJEnhKRMREZEJEHROSDkWRfAE4ypXEhwgu/hq+UiAyKyC4RuUJEmtOSH+RFRCqBbwI3A0cDHwp3/RuwG3gRsKn0LswhIjER+biIPBp+B4ZE5Ociclj79hSGDxJLh38D3gt8DDgeeAvwVWDlbAJVHVHVPiPqcucZ4AjgKOBPCDxsBh4XkZfOJsrg5QhgBbBdVZ9T1WS4/TjgblXtVNWBfARJQCyfz5aLUN+PgH8A/oXgO3ASsAO4SUQuNyYuRETipjV4MqCq/rUEXsALwMULpLkc2Ju+DpwF7AFGgXuB49I+927gD8A48CvgzwEF3hjuXx9dD7e1AtcBvcAwcB9wymL0RbbXhjp/miktcEF4/Ojr1AzbLgjTHwvcEp6zQeAuYEMk7wuAaYJA+xtgEngHEAuP+1/hufgtcFGaVgU+AHw79P0s8I9paaqA/xWe0wngOeDLkf0rgKvD7ftDDZsWOHcfCY99YoZ9nwj3vS5cnz03fwHsCr08Drw18pkY8KVQ/wTwPHBjWr7nAQ+Hn38mTF8T2X8v8A3gyvDzXcD/Bn6XQeO/Ab+MrL8uLJeR8Dt0K7Au7TOXhPr2A3cS/ElS4CjTv0eXXsYF+FeZChqeBO4AGudJM3dhjayPAj8Of5SvAh4CfhFJ8zpgBvg08FLgbIILdtYgASSAJwguxCcQXJQ/GV5sXparvrR9/xDqaElPGx7v9aGGM4HVQDx8V2BruJwgCF5d4UVpQ+jpy0B/JO8LwmPtIggULwJaCILeo8DbgWOAcwkCzZaITgW6gb8FXhweW4H2SJrrgR7g/DDNScDfh/sE+CnBBfaN4bH/O0Ggap/n3P0G+EmWfdVhOX8pXD811PQUQcB/GcHFfBQ4IkzzEYIL8KkEzXevBz4cyfMCggB7fqjxlPDcfDuS5l6CQPk1gprNBuAlpAUzYBkwAPz3cP14guDwT0Bb+LnvAr8HqsM0ZxEE8o+EeW4Jz7sPEou9dpgW4F9lKmg4GegAUuGP9d8JLugSSXM5hwaJ6dmLY7jt3PACOftj/A6RoBFuez/zB4kLwgtMVdrn7gH+ZR4PB+lL23d6eIyNWbwcpCGyXYH3pB3j/rQ0QvCv/sMR/Qq8KZLmmPC8tKV99lPAw2nH+79paZ4EPhsuHxumOSeLz1MJ/pnXp23/JnDbPOduP3D1PPsfBX4YOYZycHCrCr8/V4brV4flJVnyewZ4f9q2U8J8G8L1ewku7BVp6e4HvhpZPwcYA1aG69dxaK1lWejx7HD9l8B30tJ8AR8kFv2qwrMkUNX7ROTFwEbgDQQ/2O8BPxKRMzX8FWVgn6r2RtcJLpqrgD8S/Kv7Sdpndi4g5/UE/9xfEJHo9mUEF4N8mM0om49ceT3wOhEZSdueIOi/iPJAZPmEUMODaZ6qCAJzlIfT1vcR1GAAXhu+3zWPvjjwXNpx4gT//IvJXDmq6rSI7AJeHm76f8DdwF4RuTtc/oGqTopIC7AO+JKIfCGS36zgYzlw7h5S1Zm0414PXCkiH1bVKYJmom2q+kK4//XAsRnKqJoDZXQ88J9p+39JUOP0LAIfJJYQqjpN0GfwK+CLIvIegrbxU4CfZfnYZHo24XtFhm25UkHw7/kvM+zbv8i8Znl5qOO/8vz8LBUEnbkXZ9iXjCynVHU87XMQdKane0g/P5nOaa6DSCpCHa/PsC893yi/B16RaYeIVBM0a6UH+6yo6sMicgzwNoImt6sJLuwnccDLhwiaxtJ5NrI8mmH/jQSd638mIvcR1BLPjuyvIPjeXpXhs/25evDkhg8SS5snw/dVBeTxBEHNJMpCw2gfJPh3OKSqPQUcGwARqQX+DrhXCx+d9SBhc1haEFiIh8L3o1X1jgKOvzt8fztBTS+dBwlGpFWr6uOLyPcG4PMicqKq/jpt34eA5QRNh1FOIihfRKSKoBb67dmdqjoCfB/4voh8hqDz+c2q+gMR6QReqqpfX4TG2XwHReQHBP0ZRxP0R9wZSfIg8ErgD/PUgJ/gwOi3WU5erBaPDxJLBhH5GUH1+0GC0SDHAp8h6FjN9G8vV74EPCAiVxBciNo4UKXP9gP+DvD3wA9F5JME/3JbgbcCT6rqbfMcr1JEVhM0XdQTXLg+AdQQBIpC+QpBJ+ftIvJpoJNguO07CNrsf5XpQ6q6V0S+CXxdRD5O0FRTQ9Cx36Kqn8vl4GE+3wH+NfyHvxNoBP5EVWf7AX4C3Boe51GggeCCOD7PRflq4M+AbSJyKUF/QDXwLuB/Aleo6kNpn7lURLoIamcfIeic/1cAEfkYQTPZwwQ1p3cTNKv9PvzsJ4FviMggcDswRdAB/g5VvSiHU/Etgs7olxH0LUSb7D5DMGjgBhG5muD7vJ6gtnG1qj4NfBH4bthEtp2gk//8HI7rScd0p4h/lecFXAr8gmDUzDhBf8INwPGRNJeTYQhsWj5vJLj4r49smx0CO0FwUXsXBw+pXM+hQ2CbCEYQPUfQTPIcwb/S18zj4XIODFdNEQS4B4ArgOYMaRfdcR1uW0cQyHpDTx3huTom3H8BMJ1BXyXwcYLhwpNAH0Ez3jsXON5PgOsi6zGCYaHPhPk8S6RDn6B/5CqCi/ckwWisHxMZoprl/MXD78Hj4XdgOPxObE5LdyoHRoI9FJ6DJ4C3RdJcFO4bIhhp9ABwVlo+Z4ffh/1huoeBT0X23wtcm0VrjOC7qsCrMuzfQBB8Bgn6sfYSDMZojKT5UPi9GgvP8fvwHdeLfkl4Mj2eoiEi7yXo2GzSA52NHo/HQXxzk6dgROSjBE1WAwQdqp8DvusDhMfjPj5IeIrBKwn6IRoJ2vBvILhj2OPxOI5vbvJ4PB5PVvwD/jwej8eTFR8kPB6Px5OVw65P4t5779Vly5aZluHxeDxOsX///r729vaW9O2HXZBYtmwZbW1teX22o6ODdevWFVlReXHdg9dvHtc9uK4fzHjYvXt3R6btvrkpQixm9bwxOeG6B6/fPK57cF0/2OXBB4kI9fX1piUUjOsevH7zuO7Bdf1glwcfJCL09dk+c+fCuO7B6zeP6x5c1w92eTjs+iQKwWT03jc0wS2P9bBj7wBjUzMkYhW0H9vI5g2rWFOXe0e8Tf9A8sHrN4/rHlzXD3Z58EEiwuTkfI/jLx27OpNcueMZplMzpMJ7G/dPzbB9Tx93PTXAZe3r2bg2ty+NKQ/Fwus3j+seXNcPdnnwzU0RxsbynRQtf/YNTXDljmeYmD4QIGZJKUxMz3DljmfYNzSRU34mPBQTr988rntwXT/Y5cHXJCKsXr267Me85bEeplPpszcezHRqhlse6+GSk9cumJ8JD8XE6z+UYjVF5oovA/PY5MHXJCJ0dXWV/Zg79g4cUoNIJ6VBulww4aGYeP0Hs6szyUW37mH7nj72T82gHGiKvOjWPezqTC6Yx2LxZWAemzz4IBEhHo+X/ZhjU/PXIhabzoSHYuL1H6DYTZG54svAPDZ58EEiQm1tbdmPmYjlVgS5pjPhoZh4/QdYTFNkMfFlYB6bPPggEaG/v7/sx2w/tpFKmT9NpQTpcsGEh2Li9R+g2E2RueLLwDw2efBBIkJDQ0PZj7l5wyqqKucvhqrKCjZvWJVTfiY8FBOv/wDFborMFV8G5rHJgw8SEUwMO1tTt4zL2tezrKrikBpFpcCyqgoua1+f8ygWm4bO5YPXf4BiN0Xmii8D89jkwQeJCOPj40aOu3FtPddsauOMtmaWxyoQYHmsgjPamrlmU1vON9KBOQ/Fwus/QLGbInPFl4F5bPJgzX0SInI6cDVQCVyrqldlSPMu4HJAgUdU9a+KqcHk2OQ1dcu45OS1Od0LMR82ja/OB6//AJs3rOKupwZITWdvTlpMU2Su+DIwj00erKhJiEgl8FXgHcDxwLtF5Pi0NMcB/wicrKovBz5cbB02jU3OF9c9eP0HKHZTZK74MjCPTR6sCBLARmCvqj6tqpPAjcBZaWn+Fviqqg4CqGpxx/0B1dXVxc6y7Ljuwes/mGI2ReaKLwPz2OTBluamI4HOyPqzwIlpaV4CICL3ETRJXa6qP07PqKenhy1btlBVVUUqlWLTpk1s3bqVrq4uampqqKysZGhoiJaWFgYGBlBVWlpa6O7uprKykv7+fkZGRmhtbaW3txcRobGxkd7eXurq6kilUoyOjrJ69Wq6urqIxWLU19fT19dHfX09k5OTjI2Nze2Px+PU1tbS399PQ0MDY2NjjI+Pz+2vrq4mkUgwODhIU1MTw8PDTE5Ozu1PJBLE43GSySTNzc0kk0mmpqbm9qd7GhkZYXJycs7TihUrAJzxlEwmqaysnLecbPe0YsWKBctpMZ6mBnt5z/F1vPuliYM8MdLPaGXxPa1YsYKOjo5Ff/dsKaeZmRk6OjqK8nsy5SmZTLJy5cqSXCOyecqGqC4wELsMiMg5wOmqemG4fj5woqpeHElzBzAFvAs4Cvg5sEFVX4jmtXPnTvXTl7rrwes3j+seXNcPxqYvfai9vf2E9O22NDc9B0R7bI8Kt0V5FtimqlOq+l/A74HjiimiqampmNkZwXUPXr95XPfgun6wy4MtQeIB4DgROUZE4sB5wLa0NLcBpwKISDNB89PTxRQxPDxczOyM4LoHr988rntwXT/Y5cGKIKGq08DFwJ3Ak8DNqvpbEblCRM4Mk90J9IvIE8BPgY+palHvXbdpoo98cd2D128e1z24rh/s8mBLxzWquh3YnrbtU5FlBT4SvkqCTWOT88V1D16/eVz34Lp+sMuDFTUJW7BpbHK+uO7B6zeP6x5c1w92efBBIkIikTAtoWBc9+D1m8d1D67rB7s8+CARwaaJPvLFdQ9ev3lc9+C6frDLgw8SEZLJ4k8FWW5c9+D1m8d1D67rB7s8+CARobm52bSEgnHdg9dvHtc9uK4f7PLgg0QEm6J3vrjuwes3j+seXNcPdnnwQSLC1NSUaQkF47oHr988rntwXT/Y5cEHiQg2jU3OF9c9eP3mcd2D6/rBLg8+SESwaWxyvrjuwes3j+seXNcPdnnwQSJCTU2NaQkF47oHr988rntwXT/Y5cEHiQiVlZWmJRSM6x68fvO47sF1/WCXBx8kIgwNDZmWUDCue/D6zeO6B9f1g10efJCI0NLSYlpCwbjuwes3j+seXNcPdnnwQSLCwMCAaQkF47oHr988rntwXT/Y5cEHiQg2TOVaKK578PrN47oH1/WDXR58kIhgUxUvX1z34PWbx3UPrusHuzz4IBGhu7vbtISCcd2D128e1z24rh/s8uCDRIQVK1aYllAwrnvw+s3jugfX9YNdHnyQ8Hg8Hk9WfJCIMDIyYlpCwbjuwes3j+seXNcPdnmwJkiIyOki8jsR2Ssil86TbrOIqIicUGwNra2txc6y7Ljuwes3j+seXNcPdnmwIkiISCXwVeAdwPHAu0Xk+AzpaoEPAb8uhY7e3t5SZFtWXPfg9ZvHdQ+u6we7PFSZFhCyEdirqk8DiMiNwFnAE2nprgQ+B3ysFCJEpBTZlhXXPXj95snmYd/QBLc81sOOvQOMTc2QiFXQfmwjmzesYk3dsjKrzM7hXAYmsCVIHAl0RtafBU6MJhCR1wJrVfWHIpI1SPT09LBlyxaqqqpIpVJs2rSJrVu30tXVRU1NDZWVlQwNDdHS0sLAwACqSktLC93d3cRiMfr7+xkZGaG1tZXe3l5EhMbGRnp7e6mrqyOVSjE6Osrq1avp6uoiFotRX19PX18f9fX1TE5OMjY2Nrc/Ho9TW1tLf38/DQ0NjI2NMT4+Pre/urqaRCLB4OAgTU1NDA8PMzk5Obc/kUgQj8dJJpM0NzeTTCaZmpqa25/uaWpqis7OzjlPs6MkXPE0Pj5OX1/fvOVks6fx8XFGR0cXLCebPdXX19PR0XFQOXVMJvjnXz5HagZS4X1e+6dm2L6njzt/18cnTjmSY6onrfAUj8fp6Ogoyu/JZDlNTEyU5BqRzVM2xIY7+0TkHOB0Vb0wXD8fOFFVLw7XK4B7gAtU9RkRuRf4qKo+mJ7Xzp07ta2tLS8dHR0drFu3Lk8XduC6B6/fPOke9g1NcNGte5iYnsn6mWVVFVyzqc2KGsXhWAblYPfu3Q+1t7cf0tdrRZ8E8BywNrJ+VLhtllrgFcC9IvIMcBKwrdid13V1dcXMzgiue/D6zZPu4ZbHephOZQ8QANOpGW55rKeUsnLmcCwDk9gSJB4AjhORY0QkDpwHbJvdqapJVW1W1fWquh64HzgzU02iEFKpVDGzM4LrHrx+86R72LF3YK6JKetnNEhnA4djGZjEiiChqtPAxcCdwJPAzar6WxG5QkTOLJeO0dHRch2qZLjuwes3T7qHsan5axGLTVdqDscyMIktHdeo6nZge9q2T2VJe2opNNg0+Xi+uO7B6zdPuodErIL9OQSARMyK/5yHZRmYxI5StQSbJh/PF9c9eP3mSffQfmwjlQuMyKyUIJ0NHI5lYBJrahKmiI793j81w/LYC1aO/c6VWCxmWkJBeP3mSfewecMq7npqgNQ8o5uqKivYvGFVqaXlxOFYBvNR6vtXlnRNYldnkotu3cP2PX1z1enZsd8X3bqHXZ1JwwoXT319vWkJBeH1myfdw5q6ZVzWvp5lVRWH1CgqJRj+eln7emv+VB2OZZCN9GuYUvxr2JINEvuGJrhyxzNMTM8cMnIjpTAxPcOVO55h39CEGYF50tfXZ1pCQXj95snkYePaeq7Z1MYZbc0sj1UgwPJYBWe0NXPNpjY2rrXnwny4lkE65bqGLdnmpsWM/b7k5LXzprMJ1/9Fef3myeZhTd0yLjl5rfW/h8O5DKKU6xq2ZGsSro39zpXJyUnTEgrC6zeP6x5c1w+5eSjXNWzJBgnXxn7nytjYmGkJBeH1m8d1D67rh9w8lOsatmSDRK5jum0Z+50rNo2vzgev3zyue3BdP+TmoVzXMLeugEXEtbHfuWLT+Op88PrN47oH1/VDbh7KdQ1bskFi84ZVVFXOb9+msd+5Eo/HTUsoCK/fPK57cF0/5OahXNewJRskXBv7nSu1tbWmJRSE128e1z24rh9y81Cua9iSDRLg1tjvXOnv7zctoSC8fvO47sF1/ZC7h3Jcw5bsfRKzRMd+Dw0NWfUc93xoaGgwLaEgvH7zuO7Bdf2wOA+lvn9lSdck0lkqQ+dsxus3j+seXNcPdnnwQSLC+Pi4aQkF47oHr988rntwXT/Y5cEHiQhLZXy1zXj95nHdg+v6wS4PPkhEWCrjq23G6zeP6x5c1w92efBBIkJ1dbVpCQXjugev3zyue3BdP9jlwQeJCIlEwrSEgnHdg9dvHtc9uK4f7PLgg0SEwcFB0xIKxnUPXr95XPfgun6wy8O890mIyM9zzGdcVd9eiBAROR24GqgErlXVq9L2fwS4EJgGeoG/UdWOQo6ZTlNTUzGzM4LrHrx+87juwXX9YJeHhW6mez3w/gXSCMHFPW9EpBL4KvA24FngARHZpqpPRJL9BjhBVfeLyN8B/wycW8hx0xkeHmbFihXFzLLsuO7B6zeP6x5c1w92eVgoSPxKVa9fKBMR+asCdWwE9qrq02F+NwJnAXNBQlV/Gkl/P/CeAo95CEtlshKb8frN47oH1/WDXR7mDRKq2p5LJoU2NQFHAp2R9WeBE+dJvwX4UaYdPT09bNmyhaqqKlKpFJs2bWLr1q10dXVRU1NDZWUlQ0NDtLS0MDAwgKrS0tJCd3c31dXV9Pf3MzIyQmtrK729vYgIjY2N9Pb2UldXRyqVYnR0lNWrV9PV1UUsFqO+vp6+vj7q6+uZnJxkbGxsbn88Hqe2tpb+/n4aGhoYGxtjfHx8bn91dTWJRILBwUGampoYHh5mcnJybn8ikSAej5NMJmlubiaZTDI1NTW3P90TQGdn55yn2X8jrnhKpVL09fXNW042e5o97kLlZLOnxsZGOjo6Fv3ds8VTTU0NHR0dRfk9mfKUSqWYmJgoyTUim6dsiOoC89+VARE5BzhdVS8M188HTlTVizOkfQ9wMfBmVT1khu+dO3dqW1tbXjo6OjpYt25dXp+1Bdc9eP3mcd2D6/rBjIfdu3c/1N7efkL69oJHN4nIDwvNA3gOiD6d6qhwW/qxTgM+CZyZKUAUik3DzvLFdQ9ev3lc9+C6frDLQzGGwP6yCHk8ABwnIseISBw4D9gWTSAirwGuIQgQPUU45iEslclKbMbrN4/rHlzXD3Z5KDhIqOpni5DHNEET0p3Ak8DNqvpbEblCRM4Mk30eWAF8V0QeFpFtWbLLm2QyWewsy47rHrx+87juwXX9YJeHBeeTEJG3Ztg8BXSo6h+LJURVtwPb07Z9KrJ8WrGOlY3m5uZSH6LkuO7B6zeP6x5c1w92echl0qFvZNgWA1aJyAPAu1T1kP4DF0kmk9TU1JiWURCue/D6zeO6B9f1g10eFgwSqnpMpu0ishy4Cvg/wLuKrMsIU1NTpiUUjOsevH7zuO7Bdf1gl4e8py8N73z+R2BvEfUYxaZnuOeL6x68fvO47sF1/WCXh0I7rqc5jObJtukZ7vniugev3zyue3BdP9jlodAg8WHgoSLosAJb2gALwXUPXr95XPfgun6wy0Muo5t+AaTflh0DjgbGgT8rgS4jVFZWmpZQMK578PrN47oH1/WDXR5yaSq6NsO2aeCPwK9V1Z4nURXI0NAQDQ0NpmUUhOsevH7zuO7Bdf1gl4dcRjct+BTYw4XZB+S5jOsevH7zuO7Bdf1gl4diPLvp3cUQYgMDAwOmJRSM6x68fvO47sF1/WCXh2I8u+mTRcjDCmx4Im6huO7B6zeP6x5c1w92eSjGs5teUQwhNmBTFS9fXPfg9ZvHdQ+u6we7POQdJESkUUS2isiuYgoySXd3t2kJBeO6B6/fPK57cF0/2OVhUTfCiUgVwZDX9wFnEMwgd00JdBnBljllC8F1D16/eVz34Lp+sMtDTkFCRF5HEBjeHX7mVoJ7JP6kVHM7eDwej8c8CzY3icjjwC+ANcBFwGpV3QKMlVhb2RkZGTEtoWBc9+D1m8d1D67rB7s85NInsRxIEQSF/cBhc/NcOq2traYlFIzrHrx+87juwXX9YJeHBYOEqr6IoB9iErgJ6BKRLwPVHPq4Dqfp7e01LaFgXPfg9ZvHdQ+u6we7POQ0uklVfx42Ma0G/gF4KVAL3CsiHyihvrIiIqYlFIzrHrx+87juwXX9YJeHRQ2BVdUxVb1BVd8OrANuIJib+rCgsbHRtISCcd2D128e1z24rh/s8pD3fRKq+pyqflZVjy+mIJPYVMXLF9c9eP3mcd2D6/rBLg/zBgkRyenhfiLy/woVIiKni8jvRGSviFyaYf8yEbkp3P9rEVlf6DHTqaurK3aWZcd1D16/eVz34Lp+sMvDQvdJnCMi1wELNZD9JfDX+YoQkUrgq8DbCG7Qe0BEtqnqE5FkW4BBVT1WRM4DPgecm+8xM5FKpYqZnRFc9+D1m8d1D67rB7s8LBQkeoFv5pBPoXPtbQT2qurTACJyI3AWEA0SZwGXh8vfA74iIqJFfBLW6Ogozc3NxcrOCK578PrN47oH1/WDXR7mDRKqur5MOo4EOiPrzwInZkujqtMikgSagL5oop6eHrZs2UJVVRWpVIpNmzaxdetWurq6qKmpobKykqGhIVpaWhgYGEBVaWlpobu7m+rqavr7+xkZGaG1tZXe3l5EhMbGRnp7e6mrqyOVSjE6Osrq1avp6uoiFotRX19PX18f9fX1TE5OMjY2Nrc/Ho9TW1tLf38/DQ0NjI2NMT4+Pre/urqaRCLB4OAgTU1NDA8PMzk5Obc/kUgQj8dJJpM0NzeTTCaZmpqa25/uCaCzs3PO0+zt/a54SqVS9PX1zVtONnuaPe5C5WSzp8bGRjo6Ohb93bPFU01NDR0dHUX5PZnylEqlmJiYKMk1IpunbIgNj6QVkXOA01X1wnD9fOBEVb04kubxMM2z4fofwjQHBYmdO3dqW1tbXjo6OjpYt25dni7swHUPXr95XPfgun4w42H37t0Ptbe3n5C+vRjzSRSD54C1kfWjwm0Z04QPGqwH+ospIhaLFTM7I7juwes3j+seXNcPdnmwJUg8ABwnIseISBw4D9iWlmYbwUMGAc4B7ilmfwRAfX19MbMzgusevH7zuO7Bdf1gl4ecgoSIVIjIW8MLeNFR1WmCm/LuBJ4EblbV34rIFSJyZpjsG0CTiOwFPgIcMky2UPr6+hZOZDmue/D6zeO6B9f1g10ecnpUuKrOiMjtqlpbKiGquh3YnrbtU5HlceCdpTo+2BW988V1D16/eVz34Lp+sMvDYpqbfi4iJ5VMiQVMTrr/gFvXPXj95nHdg+v6wS4Pi5mZrgP4kYjcTjAUda4/IPqP32XGxtyfIsN1D16/eVz34Lp+sMvDYoJEArgtXD6q+FLMs3r1atMSCqZQD/uGJrjlsR527B1gbGqGRKyC9mMb2bxhFWvqlhVJZXZcLwPX9YP7HlzXD3Z5yLm5SVX/OturlALLSVdXoTeOm6cQD7s6k1x06x627+lj/9QMCuyfmmH7nj4uunUPuzqTxROaBdfLwHX94L4H1/WDXR4WNQRWRNpE5DIR+Uq4/lIReWVppJWfeLwkg7fKSr4e9g1NcOWOZ5iYniGVNrA4pTAxPcOVO55h39BEEVRmx/UycF0/uO/Bdf1gl4ecg4SIvJNgrusjgfeGm2uBL5VAlxFqa0s2eKts5Ovhlsd6mE7NzJtmOjXDLY/15JV/rrheBq7rB/c9uK4f7PKwmJrEFcBpqvp+gjmvAR4BXlV0VYbo7y/qDdxGyNfDjr0Dh9Qg0klpkK6UuF4GrusH9z24rh/s8rCYILEKeDRc1si7+Yc/FYmGhgbTEgomXw9jU/PXIhabLl9cLwPX9YP7HlzXD3Z5WEyQeAg4P23becCu4skxi03DzvIlXw+JWG5fhVzT5YvrZeC6fnDfg+v6wS4PixkC+0HgLhHZAtSIyJ3AS4C3l0SZAcbHx01LKJh8PbQf28j2PX3zNjlVSpCulLheBq7rB/c9uK4f7PKQc5BQ1T0i0gb8OXAHwQ11d6jqSKnElRubxibnS74eNm9YxV1PDZCazt6cVFVZweYNq/KVlhOul4Hr+sF9D67rB7s8LGZ00ytVdb+q3qyqn1fVGw+nAAF2jU3Ol3w9rKlbxmXt61lWVUFl2mS1lQLLqiq4rH19yW+os6kM9g1N8OX7Ojn7+kf402t/w9nXP8KX7+ucdxiwTfrzxXUPrusHuzwsprnpDhGpIRgG+7Pw9ZtiP67bJNXV1aYlFEwhHjaureeaTW1G77i2pQx2dSa5csczTKcO3Dcye2PhXU8NcFn7ejauPfQhbLboLwTXPbiuH+zysJjmpqNF5EXAKcCbCR7t3SQiv1TVPy+VwHKSSCRMSyiYQj2sqVvGJSev5ZKT1y6cuATYUAbRGwvTSSmkwhsLr9nUdkjgtEF/objuwXX9YJeHRQ1VUdWngV8BO4H7Ce6XKG0jdRkZHBw0LaFgXPdgg/5Cbiy0QX+huO7Bdf1gl4fF9EncJCJ/BL4FvAj4DrBeVTeWSly5aWpqMi2hYFz3YIP+Qm4stEF/objuwXX9YJeHxdQkXgvMENxl/QjwsKoOl0SVIYaH3bfjugcb9BdyY6EN+gvFdQ+u6we7PCzmKbDHAW8A7gHeSDC3xO9F5NpSiSs3Nk30kS+ue7BBfyE3Ftqgv1Bc9+C6frDLw2L7JJ4HfgfsBZ4BVgPvKL4sM9g0NjlfXPdgg/72YxsPGQacTrYbC23QXyiue3BdP9jlYTF9EttEZAC4naDp6QfA61T1yFKJKzc2jU3OF9c92KB/84ZVVFXO/9PIdmOhDfoLxXUPrusHuzwspiZxK0FQWKeq56vqtar6VKECRKRRRO4WkafC90OebCUirxaRnSLyWxF5VETOLfS4mbBp2Fm+uO7BBv2F3Fhog/5Ccd2D6/rBLg+L6ZO4DugUkVNE5N3h+2JuxsvGpcCOsM9jR7iezn7gvar6cuB04F9EZGURjn0QNk30kS+ue7BF/+yNhWe0NbM8VoEAy2MVnNHWzDWb2jLeSAf26C8E1z24rh/s8pDzRT58btMPCOa67gTWAuMi8heq+mQBGs4CTg2XrwfuBT4RTaCqv48s7xORHqAFeKGA4x5CMplk5cqVxcyy7LjuwSb9+dxYaJP+fHHdg+v6wS4Pi6kJ/Cvw78AXZh/FISIfDbe/pQANrWGHOEAX0DpfYhHZCMSBP2Ta39PTw5YtW6iqqiKVSrFp0ya2bt1KV1cXNTU1VFZWMjQ0REtLCwMDA6gqLS0tdHd3E4vF6O/vZ2RkhNbWVnp7exERGhsb6e3tpa6ujlQqxejoKKtXr6arq4tYLEZ9fT19fX3U19czOTnJ2NjY3P54PE5tbS39/f00NDQwNjbG+Pj43P7q6moSiQSDg4M0NTUxPDzM5OTk3P5EIkE8HieZTNLc3EwymWRqampuf7qnVCpFZ2fnnKcVK1YAOONpamqKvr6+ecvJZk9TU1OMjo4uWE42e6qvr6ejo2PR3z1bPFVXV9PR0VGU35MpT1NTU0xMTJTkGpHNU9Zrbq6PXgo7rVtUNRXZVgX0quq8M2SIyE8IRkKl80ngelVdGUk7mC0/ETmCoKbxPlW9P1OanTt3altb2wJuMrNv3z7WrFmT12dtwXUPXr95XPfgun4w42H37t0Ptbe3n5C+fTE1iX0Ez2y6J7LtTeH2eVHV07LtE5FuETlCVZ8Pg0DGSZRFpA74IfDJbAGiUKampkqRbVlx3YPXbx7XPbiuH+zysJgg8T+AbSJyB9ABrAP+DHhPgRq2Ae8Drgrfb09PICJx4PvAt1T1ewUeLys2jU3OF9c9eP3mcd2D6/rBLg+LGd20DXgN8DhQG76/TlUPuagvkquAt4nIU8Bp4ToickLkbu53ETx99gIReTh8vbrA4x6CTWOT88V1D16/eVz34Lp+sMvDgjUJEVkO/E/gFcBu4LOqmn3WlUWiqv1Ae4btDwIXhss3ADcU65jZqKmpKfUhSo7rHrx+87juwXX9YJeHXGoSXwX+AtgDnAN8oaSKDFJZWWlaQsG47sHrN4/rHlzXD3Z5yCVInA68XVU/TvCcpsNigqFMDA0NmZZQMK578PrN47oH1/WDXR5yCRI1s/cxqGonkPlW08OAlpYW0xIKxnUPXr95XPfgun6wy0Muo5uqROQtgGRZR1XvyfhJxxgYGGD58uWmZRSE6x68fvO47sF1/WCXh1yCRA/wzch6f9q6EsxU5zy53lhoM6578PrN47oH1/WDXR4WDBKqur4MOqzApipevrjuwes3j+seXNcPdnlY1KRDhzvd3d2mJRSM6x68fvO47sF1/WCXBx8kIiz0oCsXcN2D128e1z24rh/s8uCDhMfj8Xiy4oNEhJGREdMSCsZ1D16/eVz34Lp+sMuDDxIRWlvnncrCCVz34PWbx3UPrusHuzz4IBGht7fXtISCcd2D128e1z24rh/s8uCDRAQRWTiR5bjuwes3j+seXNcPdnnwQSJCY2OjaQkF47oHr988rntwXT/Y5cEHiQg2VfHyxXUPXr95XPfgun6wy4MPEhHq6upMSygY1z14/eZx3YPr+sEuDz5IREilUqYlFIzrHrx+87juwXX9YJcHHyQijI6OmpZQMK578PrN47oH1/WDXR58kIhg0+Tj+eK6B6/fPK57cF0/2OXBB4kINk0+ni+ue/D6zeO6B9f1g10ejAcJEWkUkbtF5KnwvWGetHUi8qyIfKUUWmKxWCmyLSuue/D6zeO6B9f1g10ejAcJ4FJgh6oeB+wI17NxJfDzUgmpr3d/ZlbXPXj95nHdg+v6wS4PNgSJs4Drw+XrgbMzJRKR1wGtwF2lEtLX11eqrMuG6x68fvO47sF1/WCXh1ymLy01rar6fLjcRRAIDkJEKoAvAu8BTpsvs56eHrZs2UJVVRWpVIpNmzaxdetWurq6qKmpobKykqGhIVpaWhgYGEBVaWlpobu7GxGhv7+fkZERWltb6e3tRURobGykt7eXuro6UqkUo6OjrF69mq6uLmKxGPX19fT19VFfX8/k5CRjY2Nz++PxOLW1tfT399PQ0MDY2Bjj4+Nz+6urq0kkEgwODtLU1MTw8DCTk5Nz+xOJBPF4nGQySXNzM8lkkqmpqbn96Z6mpqbo7Oyc8zT7XHpXPI2Pj9PX1zdvOdnsaXx8nNHR0QXLyWZPy5cvp6OjY9HfPVs8xWIxOjo6ivJ7MulpYmKiJNeIbJ6yIeWYS1VEfgJk6q7/JHC9qq6MpB1U1YP6JUTkYmC5qv6ziFwAnKCqF2c61s6dO7WtrS0vnT09PaxatSqvz9qC6x68fvO47sF1/WDGw+7dux9qb28/IX17WWoSqpr137+IdIvIEar6vIgcAfRkSPYG4E0i8gFgBRAXkRFVna//YtGMjY0VMzsjuO7B6zeP6x5c1w92ebChuWkb8D7gqvD99vQEqvrfZpcjNYmiBgiwa2xyvrjuwes3j+seXNcPdnmwoeP6KuBtIvIUQX/DVQAicoKIXFtOITaNTc4X1z14/eZx3YPr+sEuD8ZrEqraD7Rn2P4gcGGG7dcB15VCSzweL0W2ZcV1D16/eVz34Lp+sMuDDTUJa6itrTUtoWBc9+D1m8d1D67rB7s8+CARob+/37SEgnHdg9dvHtc9uK4f7PLgg0SEhoasTwRxBtc9eP3mcd2D6/rBLg8+SESwadhZvrjuwes3j+seXNcPdnnwQSLC+Pi4aQkF47oHr988rntwXT/Y5cEHiQg2jU3OF9c9eP3mcd2D6/rBLg8+SESwaWxyvrjuwes3j+0e9g1N8OX7Ojn7+kf402t/w9nXP8KX7+tk39AEYL/+XLDJg/H7JGyiurratISCcd2D128emz3s6kxy5Y5nmE7NkAofO7d/aobte/q466kBLmtfzzqL9eeKTWXgaxIREomEaQkF47oHr988tnrYNzTBlTueYWL6QICYJaUwMT3DlTueIZly/7+vTWXgg0SEwcFB0xIKxnUPXr95bPVwy2M9TKdm5k0znZrhe492l0lR6bCpDHyQiNDU1GRaQsG47sHrN4+tHnbsHTikBpFOSuHX3VPlEVRCbCoDHyQiDA8Pm5ZQMK578PrNY6uHsan5axGzjOeYzmZsKgP3G++KyOTkpGkJBeO6B6/fPLZ6SMQq2J9DAFhWKWVQE/SR3PJYDzv2DjA2NUMiVkH7sY1s3rCKNXXLCsrbpjLwNYkINo1NzhfXPXj95rHVQ/uxjSx0/a8UaH/xypJr2dWZ5KJb97B9Tx/7p2ZQDoyyuujWPezqTBaUv01l4INEBJvGJueL6x68fvPY6mHzhlVUVc5/yaqqrOBPWkrb3JTrKKvZ+zbywaYy8EEigk3DzvLFdQ9ev3ls9bCmbhmXta9nWVXFITWKSoFlVRXBfRJNK0qqI9dRVrc8lmkm5tywqQx8kIhg00Qf+eK6B6/fPDZ72Li2nms2tXFGWzPLYxUIsDxWwRltzVyzqY2Na+tLrj/XUVY79g7kfQybysB3XEdIJpOsXLnStIyCcN2D128e2z2sqVvGJSev5ZKT12bcX2r9uY6yyjVdJmwqA1+TiNDc3GxaQsG47sHrN4/rHkqtPxHL7bKZa7pM2FQGPkhESCYLG5FgA6578PrN47qHUuvPeZTVsY15H8OmMjAeJESkUUTuFpGnwveMUzKJyNEicpeIPCkiT4jI+mJrmZpy/05N1z14/eZx3UOp9ec6ymrzhlV5H8OmMjAeJIBLgR2qehywI1zPxLeAz6vqy4CNQP5DB7Jg09jkfHHdg9dvHtc9lFp/rqOsCrmhzqYysCFInAVcHy5fD5ydnkBEjgeqVPVuAFUdUdX9xRZi09jkfHHdg9dvHtc9lEN/LqOsCsGmMrBhdFOrqj4fLncBrRnSvAR4QURuBY4BfgJcqqqp9IQ9PT1s2bKFqqoqUqkUmzZtYuvWrXR1dVFTU0NlZSVDQ0O0tLQwMDCAqtLS0kJ3dzczMzP09/czMjJCa2srvb29iAiNjY309vZSV1dHKpVidHSU1atX09XVRSwWo76+nr6+Purr65mcnGRsbGxufzwep7a2lv7+fhoaGhgbG2N8fHxuf3V1NYlEgsHBQZqamhgeHmZycnJufyKRIB6Pk0wmaW5uJplMMjU1Nbc/3dP4+DidnZ1znlasCMaMu+Jp//799PX1zVtONnvav38/o6OjC5aTzZ7i8TgdHR2L/u7Z4klE6OjoKMrvaSFPZx8Nf7muKeJpBanJYTo6ugrytH//fiYmJkpyjcjmKRuiusCA3yIgIj8BMtWfPglcr6orI2kHVfWgfgkROQf4BvAa4I/ATcB2Vf1GeoY7d+7Utra2vHQODg7S0JCxS8QZXPfg9ZvHdQ+u6wczHnbv3v1Qe3v7Cenby9LcpKqnqeorMrxuB7pF5AiA8D1TX8OzwMOq+rSqTgO3Aa8tts6hoaFiZ1l2XPfg9ZvHdQ+u6we7PNjQJ7ENeF+4/D7g9gxpHgBWikhLuP5W4IliC2lpaVk4keW47sHrN4/rHlzXD3Z5sCFIXAW8TUSeAk4L1xGRE0TkWoCw7+GjwA4ReQwQ4OvFFjIwkP9t9Lbgugev3zyue3BdP9jlwXjHtar2A+0Ztj8IXBhZvxt4ZYm1lDL7suC6B6/fPK57cF0/2OXBhpqENdhUxcsX1z14/eZx3YPr+sEuDz5IROjudn8Cddc9eP3mcd2D6/rBLg8+SERYaLywC7juwes3j+seXNcPdnkw3ifh8WSilPMHezye3PE1iQgjIyOmJRSM6x5GRkZKPn9wKXH9/IP7HlzXD3Z58EEiQmtrpieCuIXrHmaWN5R8/uBS4vr5B/c9uK4f7PLgg0SE3t5e0xIKxnUP/7n72ZLPH1xKXD//4L4H1/WDXR58kIggssBMIg7guof7902WfP7gUuL6+Qf3PbiuH+zy4INEhMbG/GeSsgXXPUwsFCFCCpk/uJS4fv7BfQ+u6we7PPggEcGmKl6+uO5hWWVu6QqZP7iUuH7+wX0PrusHuzzY+UszRF1dnWkJBeO6h1PW1ZZ8/uBS4vr5B/c9uK4f7PLgg0SEVOqQOYycw3UP73jxipLPH1xKXD//4L4H1/WDXR58kIgwOjpqWkLBuO5hBRMlnz+4lLh+/sF9D67rB7s8+DuuI9g0+Xi+uO5h9erVrFu2jGs2tTl5x7Xr5x/c9+C6frDLgw8SEbq6uli3bp1pGQXhuodZ/WvqlnHJyWu55OS1piUtCtfPP7jvwXX9YJcH39wU4bbbbjMtoWBc9+D1m8d1D67rB7s8+CAR4dZbbzUtoWBc9+D1m8d1D67rB7s8+CARYXp62rSEgnHdg9dvHtc9uK4f7PIgNk2TVwx27NjRC3Tk89mBgYHmxsbGviJLKiuue/D6zeO6B9f1gzEP69rb2w+ZEu+wCxIej8fjKR6+ucnj8Xg8WfFBwuPxeDxZWZJBQkROF5HficheEbk0w/5lInJTuP/XIrLegMys5KD/AhHpFZGHw9eFJnRmQ0S+KSI9IvJ4lv0iIv839PeoiLy23BoXIgcPp4pIMlIGnyq3xvkQkbUi8lMReUJEfisiH8qQxtpyyFG/7WVQLSK7ROSR0MM/ZUhj/lqkqkvqBVQCfwBeBMSBR4Dj09J8APhauHwecJNp3YvUfwHwFdNa5/FwCvBa4PEs+88AfgQIcBLwa9Oa8/BwKnCHaZ3z6D8CeG24XAv8PsP3yNpyyFG/7WUgwIpwOQb8GjgpLY3xa9FSrElsBPaq6tOqOgncCJyVluYs4Ppw+XtAu9gzC0gu+q1GVX8OzDdr0FnAtzTgfmCliBxRHnW5kYMHq1HV51V1d7g8DDwJHJmWzNpyyFG/1YTndXYy61j4Sh9JZPxatBSDxJFAZ2T9WQ79cs2lUdVpIAk0lUXdwuSiH2Bz2ETwPRFx69kWuXu0nTeETQk/EpGXmxaTjbAJ4zUE/2SjOFEO8+gHy8tARCpF5GGgB7hbVbOWgalr0VIMEkuBHwDrVfWVwN0c+CfiKR+7gXWq+irgy8BtZuVkRkRWALcAH1bVIdN6FssC+q0vA1VNqeqrgaOAjSLyCsOSDmEpBonngOg/66PCbRnTiEgVUA/0l0XdwiyoX1X7VXUiXL0WeF2ZtBWLXMrIalR1aLYpQVW3AzERaTYs6yBEJEZwgf2OqmZ6DoTV5bCQfhfKYBZVfQH4KXB62i7j16KlGCQeAI4TkWNEJE7QGbQtLc024H3h8jnAPRr2HFnAgvrT2o3PJGivdYltwHvD0TUnAUlVfd60qMUgIqtn245FZCPBb82WPxqE2r4BPKmqX8qSzNpyyEW/A2XQIiIrw+UE8DZgT1oy49eiJfeocFWdFpGLgTsJRgp9U1V/KyJXAA+q6jaCL9+3RWQvQefkeeYUH0yO+j8oImcC0wT6LzAmOAMi8p8EI0+aReRZ4H8RdNqhql8DthOMrNkL7Af+2ozS7OTg4Rzg70RkGhgDzrPojwbAycD5wGNhmzjA/wCOBifKIRf9tpfBEcD1IlJJEMBuVtU7bLsW+cdyeDwejycrS7G5yePxeDw54oOEx+PxeLLig4TH4/F4suKDhMfj8Xiy4oOEx+PxeLLig4RnySPBU3N/aVpHoYjIvZLHE3/Dp6XOiMiIiKTfzFWInutE5NPh8kvC/FP5aPSYY8ndJ+HxFIKIXA4cq6rvMa2lyOxT1aNKlbmq/h5YISL3luoYntLgaxIej2dBwkdCeJYgPkh4So6I/LWI/CCy/pSIfDey3ikirw6Xrw7Xh0TkIRF5U7h9jYiMiUhj5HOvEZG+8Bk+iMjfiMiTIjIoIneKyLpIWhWRD4rI0+FnPi8iGb//82g4neCu3nPDppNHwu31IvINEXleRJ4TkU+Hd9FmynujiOwUkRfC9F8JH68S1fn+8By9ICJfjTxaolJEvhjq/y8RuThMn/ECPt/5WIiwCe4+Efk/ItIPXC4iLxaRe0SkP9TwndnHSkTKY7eIDIvITUB1rsfz2IsPEp5y8DPgTSJSISJrCCZLegOAiLwIWAE8GqZ9AHg10Aj8B/BdEalW1X3ATmBzJN+/Ar6nqlMichbBBXwT0AL8AvjPNB1/CZxAMFnQWcDfZNGbTcOPgc8QTPyyIny6KMB1BI9AOZbgkdVvB7K1u6eAvweaw3PQTjCxTJQ/B14PvBJ4F/Cn4fa/Bd4RanstcHaWY5Dj+ViIE4GngVbgfxNMkvNZYA3wMoIHz10eHi9O8JTVbxOct+9ycFl5XKXcsxz519J8ETwT/7UEz575d2AX0EbwPKBt83xuEHhVuHwhwQPOILhgdQKnhOs/ArZEPldB8LyhdeG6AqdH9n8A2BEuXwD8MkcNlwM3RPa1AhNAIrLt3cBPczwvHwa+H1lX4I2R9ZuBS8Ple4CLIvtOC9NXhev3Ahfmcj7SNJwKPJu27QLgjwtoPxv4Tbh8CrCP8FE/4bZfAZ9O+8ycRv9y4+VrEp5y8TOCi9Ep4fK9wJvD189mE4nIR8MmkqSIvEDwaOTZxzvfQjCJzBFhPjME/5AB1gFXh000LxA8DE04eJKc6AQ6HQT/iA9hAQ3prCN4sN/zkWNfA6zKkvdLROQOEekSkSGCmkl63l2R5f0ENS1CvVEP0eVMuhY6HwtxUP4i0ioiN4ZNakPADRHta4DnNIwEIR2LOJbHUnyQ8JSL2SDxpnD5Z6QFibDt/+METSwNqrqSYCYuAVDVQeAu4FyCpqYbIxelToJ/2Ssjr4Sq/iqiITo3wtEE/3wPYiENHDq9ZCdBTaI5ctw6Vc02C9q/ETwO+jhVrSNoEsp1OsrnCeZ0yOQnnVzOx0Kke/1MuG1DqP09HND+PHDkbP9JyNGLOJbHUnyQ8JSLnwFvIWiWeZagBnA6wVSMvwnT1BK07fcCVSLyKaAuLZ//AN5L8Bjo/4hs/xrwjxJOURl2Jr8z7bMfE5EGCaZz/RBwUwadC2noBtbPdnprML/CXcAXRaQu7Hd5sYi8Oct5qAWGgBERaQP+Lku6TNwMfEhEjgw7jD8xT9pczsdiqQVGgKSIHAl8LLJvJ8F5+6CIxERkE8F87B7H8UHCUxY0GCc/Qtg8pMFUk08D96lqKkx2J/Bj4PcETRXjHNqksg04DuhS1Uci+X8f+BxwY9gU8jhBJ2+U24GHgIeBHxI8qz+dhTTMjsrqF5Hd4fJ7CTrjnyDov/gewVwBmfgoQS1oGPg6mQNVNr5OEJAeJQis2wkuzKn0hDmej8XyTwT9SkmC8zc3G5yqThJ0kl9A0LR1bnS/x138fBKeJYGIKEETz17TWoqFiLwD+Jqq5jy0NUs+pxAExwngXFW9sxj60o5xHMGosTjwAVW9rtjH8JQGHyQ8S4LDIUhIMMXlWwhqE60EHfn3q+qHTeryHN745iaPxx2EoMlnkKC56UngU0YVeQ57fE3C4/F4PFnxNQmPx+PxZMUHCY/H4/FkxQcJj8fj8WTFBwmPx+PxZMUHCY/H4/FkxQcJj8fj8WTl/wNlM7wbEgRPqwAAAABJRU5ErkJggg==","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["power_o_matrix = np.zeros_like(thetas)\n","power_e_matrix = np.zeros_like(thetas)\n","single_diff_matrix = np.zeros_like(thetas)\n","\n","for i,angle in enumerate(thetas):\n","    # The polarization state analyzer for both beams of the Wollaston prism\n","    PSA_o = wollaston(0) @ linear_retarder(angle, np.pi)\n","    PSA_e = wollaston(1) @ linear_retarder(angle, np.pi)\n","    power_o_matrix[i] = (PSA_o[0, :] @ S_to_measure)\n","    power_e_matrix[i] = (PSA_e[0, :] @ S_to_measure)\n","    single_diff_matrix[i] = power_o_matrix[i] - power_e_matrix[i]\n","\n","plt.style.use('bmh')\n","plt.figure()\n","plt.title('Single Difference Observed')\n","plt.plot(thetas, single_diff_matrix, marker = 'o', linestyle = 'None', \n","    markersize = 10)\n","plt.xlabel('waveplate angle [rad]')\n","plt.ylabel('Power [A.U.]')\n","plt.show()"]},{"cell_type":"markdown","metadata":{},"source":["We can now pass the full list of powers to `dual_channel_polarimeter` instead of a single Stokes vector to perform a real measurement "]},{"cell_type":"code","execution_count":5,"metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["Stokes Vector Measured\n","[1.         0.32205895 0.6191101  0.        ]\n","Percent Difference\n","[2.35802600e-06 1.15021692e-06]\n"]}],"source":["S_out = dual_channel_polarimeter(thetas, power_o = power_o_matrix, \n","    power_e = power_e_matrix, sub_method = \"single_difference\", normalized = False)\n","print('Stokes Vector Measured')\n","print(S_out)\n","print('Percent Difference')\n","print(100 * (S_to_measure[1 : -1] - S_out[1 : -1]) / S_to_measure[1 : -1])"]},{"cell_type":"markdown","metadata":{},"source":["# Double Differencing"]},{"cell_type":"markdown","metadata":{},"source":["Let's do the same with the double difference, where we use $(I_o(\\theta_2) - I_e(\\theta_2)) - (I_o(\\theta_1) - I_e(\\theta_1))$ as the observable instead, with $\\theta_1$ and $\\theta_2$ corresponding to two different HWP angles. \n","\n","The double difference power also follows a sinusoid with two frequencies, where:\n","\n","$$P \\propto Q(cos(4\\theta_2) - cos(4\\theta_1)) + U(sin(4\\theta_2) - sin(4\\theta_1))$$\n","\n","Let's test the double difference method with the same randomized normalized input Stokes vector"]},{"cell_type":"code","execution_count":6,"metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["Stokes Vector Measured\n","[1.         0.32205896 0.6191101  0.        ]\n","Percent Difference (Q, U only)\n","[-1.59943387e-06  2.25577236e-06]\n"]}],"source":["S_out = dual_channel_polarimeter(thetas, S_in = S_to_measure, \n","    sub_method = \"double_difference\", normalized = False)\n","print('Stokes Vector Measured')\n","print(S_out)\n","print('Percent Difference (Q, U only)')\n","print(100 * (S_to_measure[1 : -1] - S_out[1 : -1]) / S_to_measure[1 : -1])"]},{"cell_type":"markdown","metadata":{},"source":["Great, let's test it with a list of measured powers as well!"]},{"cell_type":"code","execution_count":7,"metadata":{},"outputs":[{"data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAZAAAAEXCAYAAACDChKsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+OklEQVR4nO29e5hcVZmo/37pO01307d0EwgJCNLKREExI6KINgh6/AmTeME5KjjhDGpEPR5v6MBwxAsef0fHUWdkBAa8ggJqZKKgweCFKJdwCWCACImBpPpOX5LqW/V3/ti7m6JS3V3dVbv2qi/rfZ56UnvvtXd9716d+mqvvfZaoqp4PB6Px7NQlsQdgMfj8XhKE59APB6Px7MofALxeDwez6LwCcTj8Xg8i8InEI/H4/EsCp9APB6Px7MofALxFAUR2Ski/zRPmc0icnURY1opIioir05bt1xENonIPhHRcF29iPxERAbD8iuLFaPL5FKnpYhVryjwCeQgRESuC78IVUQmRKRXRH4vIp8Qkdq448uHtKQw/donIk+IyHdF5FUZxXcDhwN/Slv3aWApcGK4DeD9wCnAq8N1uyOViBkROV1EfiEi/SIyJiKPi8gXRKQu7tg8buETyMHL7wi+DFcArwO+D3wQ2CoibXEGViDOIfA7gSABCPB7EfnodAFVTalqQlUn0vY7DrhbVZ9Q1UTaukdUdVtYPrWYgESkclEmRURE1gGbgB1AJ/BCgqT6duAPIlIfY3hIQEWcMXjSUFX/OshewHXAr7OsPwLoB/4zbV0FcCXwDDAOPAr8fcZ+CrwrY92vgevSlncCnweuBoaAXuALwJK0MpuBqzOOczGwHRgFngA+A5TP4bYyjOfVWbb9H2ACeEG2suH79Nd1Ydzp6zannZfLgafC2B4BLspyXj4E/AAYBG4M158J/AFIhuf1P4HmzPoB/hHYFZ6vDUBbxvHPIPghsD88/p3TbuH284AHwvh2Al8Bauc4d8vCsv+WZduKMN5/XWCdngPcH8b4LHA3cFLa9mOBm8NtA8DtwKq07RcAkwQ/cu4n+Bu8ODy3r8qI8W/D9ceFy4cCXwvP8f5w/zUZ+7wUuAsYI/j7envo9U9x/z8thVfsAfhXDJU+SwIJt309/DJaEi5/GegD3sZzv0angM60fXJNIEPAZ4HjgXcD+4APp5XZTFoCIfiC3gX8HXA08Cbgr8AVc7itZPYE0hLG/rFsZYH28Mvk++H7BqAVuBH4bbiuKe0cPgS8IYztHeGX4LqM89JHcGX3AoIrmdeHX2YXh8uvAH5D8OUvacceBH4I/A1B89lTwHfTjn0GkAL+JfwS7ADWAR3h9gsIvpDfDRwDnBbG+905zt2Hw5iPnGX7NQRJYjrOOes0PF/jwCfCc/Qi4O8JEwTQBiSAfwdWhcf4enjOWtM8pggSz+tCl9awnv49I75/A+4K30t4XjcTND0eQ5CQxwn/doEaguSyMTyHpwD3hPXjE0gOr9gD8K8YKn3uBPK+8EtkKXAIwS+zD2SU+QlwR9pyrgnkdxllvgDsTlveTJhAws/eD5ydsc97gGfncFvJLAkk3J4g/IWdrSzZr4Ked77CL8Mpwi/rtPWXAQ9knJdrMspsBq7MWHdUWPbEtM/rBqrSynwS2Ju2/Dvg1jnOw07gfRnrTgs/p3GWff4NGJzjmB8N929N+4xZ6xQ4KSy/cpbjXQ78MWOdAH8BPhIuXxAe4zVZ/k77gcpwuZIg8VwULp9OcDXVkLHftcBPw/cXAiPp54MgYSs+geT08vdAPJlI+K8SNC9UEvz6TudOgnsLC2VLxvIfgCNnaVc/geAX4s0iMjL9Aq4CGkSkdRGfD4GfLnLfaU4Oj3NvRmyfJriqSOfujOVXAB/J2O/RcFv6vttVdSxteQ/BL/ZpXk7Q3HMA4blZAXwl43N+ERY5NjfNnJirTh8CbgMeDnuxfVhElqeVfQXw8owYhwkSe+Z5vCdj+UaCHxlvDpffDNSG66ePXQk8k3H8d6Ud+8XAn1V1YPqgqvowwdWfJwfK4w7A4xwnEPwH6uO5XkjzoTyXeKbJ90bn9I+btwGPZ9nev9ADhl+srcCTecQFz8X2KoKrpHQyk9O+LPt+CfhuluMm0t6PZzlu5jmeL74PEzTjZPL0LPs9DtSLyHJVzdbT7ASCv4veXIJQ1ZSIvJHgy/wMYC1wpYi8TVVvDePcRNDEl0n6l3hKVUczjj0gIj8nuCK9Jfx3g6o+GxZZEh7jFVmOnXluPYvEJxDPDCJyBPDfgVtUdUpEdhA0YZ0GPJxW9LUZy90EN2Cnj1NF8OvuqYyPeGXG8quAZ1R1KEs4jxA0QRyjqhsXoZONjxPcN/hJnse5L/z3qPCLcCHcC5ygqjsKEMMbgH/N3KCqXSKyGzheVb+9gGP+mKDDxKcJeq7NICIrCO5ffFvDtp6QOes0LHt3+PqCiPwSeC9wK8G5uAB4OjNB5Mj1wC0icjzB/bE1advuBQ4DqsOrimw8CvyjiBw2nXhE5ASCe1+eHPAJ5OClUkTaCX6pNRPcaLyEIBlcAqCq+0XkX4ErRKQHeBB4K0HPmjPTjvVr4H0i8luCJojPEDQfZHKiiFxO0CvpZIJfyJdmC05VR0TkCwRfOhp+RjnBzdaTVPWT8/g1hX5VBE02FxAkx/+lqnldgajqDhG5Fvi2iHyCoBmnlqBZqVVVvzTH7pcBt4vIV4DvEJyv4wiutD6oqskcw7gC+IWI/AtBu/4YwU3gLar6GEEdXCMiA8DPCHqfvQh4o6peNIvXMyLyIeDfRWSSoHfVAMGv+CsJeillPmA3a52Gz910EjS17Q09X0JwMx7gGwQ3/n8mIp8jeL7mSOCNwH+p6l3znINfhvHdEP77y7RtdxD8zdwS1tFDQCNBghsNE+sPwvP4PRH5DEGT6dcIept5ciHumzD+VfwXwU1aDV+TBM0SvyfoLVObUTaXbrztwM8JeuTsJvj1Ols33v8My/WFx52vG++FPNcVdYDgob/3z+G2Ms1NCZqYdhA0GWV2+5wuu6Cb6OG6svB8bQ/PSy/BvaG3pZU5oHNBuP414fkZJmji+jNBb6ryOT7vXYQ/6NPWnUWQvJIEzTW/Ibhim95+brh9f3jOHwAuy+Hvo5Pg3sWzodsTBDfH6zLKzVmnBE1eGwma5sYIetR9mfDGd1hmBUGvt560Mt8Djg63XwBMzhHrV8Pz/NUs22rCeJ4KPRIESeb1aWVOCs/RGMHN+/Pw3Xhzfk13x/N4PB6PZ0H4Xlgej8fjWRQ+gXg8Ho9nUfgE4vF4PJ5F4ROIx+PxeBaFTyAej8fjWRQHzXMgmzdv1qqqqrjD8Hg8npJi//79vZ2dnVmHDjpoEkhVVRUdHR2xff6uXbtYsWJFbJ8fFVa9wK6bVS+w6xan19atW3fNts03YRWJigqbc+BY9QK7bla9wK6bq14+gRSJhgabw+tY9QK7bla9wK6bq14+gRSJ3t6cBjAtOax6gV03q15g181Vr4PmHkjcuPoLIl+seoFdN6teYNdtMV57hsa4eVs3m3b0k5yYoqZiCZ3HNrF21VKW1RemQ5FPIEVifNzmFARWvcCum1UvsOu2UK+7dw9yxaadTKamSIXDHe6fmGLj9l5uf6KfSztXsnp5/snWN2EViWTS5gjRVr3ArptVL7DrthCvPUNjXLFpJ2OTzyWPaVIKY5NTXLFpJ3uGxrIfYAH4BFIk2tvb4w4hEqx6gV03q15g120hXjdv62YyNTVnmcnUFDdv6843LJ9AikUikZi/UAli1Qvsuln1ArtuC/HatKP/gCuPTFIalMsXn0CKRGVltgn6Sh+rXmDXzaoX2HVbiFdyYu6rj4WWmwufQIpEXV1d3CFEglUvsOtm1Qvsui3Eq6Yit6/1XMvNhU8gRaKvry/uECLBqhfYdbPqBXbdFuLVeWwTZTJ3mTIJyuWLTyBForGxMe4QIsGqF9h1s+oFdt0W4rV21VLKy+b+ai8vW8LaVUvzDcs/B1Iskskk9fX1cYdRcKx6QfRuxXjQKxu+zkqPhXgtq6/i0s6VBzwHAsGVR3nZEi7tXFmQvzGfQIrE6Oho3CFEglUviNatWA96ZcPXWemxUK/Vyxu4ak1H5D9QRHWe/l5G2LJli8Y5nPvY2BgW5yOx6gXRue0ZGuOiW7YzNjl7L5iq8iVctaYjkisRX2elR5xeW7duva+zs/PkbNv8PZAi4funlx5RuRXzQa9s+DorPVz18gmkSFRXV8cdQiRY9YLo3Ir5oFc2fJ2VHq56+QRSJGpqauIOIRKsekF0bsV80Csbvs5KD1e9nEkgInK2iDwmIjtE5FNZtn9VRB4IX4+LyLNp21Jp2zYUNfAcGRgYiDuESLDqBdG5FfNBr2z4Ois9XPVyoheWiJQB3wTOBJ4G7hGRDar66HQZVf2faeUvBk5KO0RSVU8sUriLorm5Oe4QIsGqF0Tn1nlsExu3987ZjFWoB72y4eus9HDVy5UrkNXADlV9UlXHgRuAc+Yo/07gh0WJrEAMDw/HHUIkWPWC6NyK+aBXNnydlR6uermSQI4AdqctPx2uOwARWQEcDdyRtrpaRO4VkT+KyLmRRZkHfqKb0iMqt+kHvarKlxww5ESZBF14C/WgVzZ8nZUerno50YS1QM4DblLVVNq6Far6jIgcA9whIttU9S/pO3V3d7Nu3TrKy8tJpVKsWbOG9evXk0gkqK2tpaysjKGhIVpbW+nv70dVaW1tpauri0MPPRSAkZER2tra6OnpQURoamqip6eH+vp6UqkU+/bto729nUQiQUVFBQ0NDfT29tLQ0EBFRQW7du2a2V5ZWUldXR19fX00NjaSTCYZHR2d2V5dXU1NTQ0DAwM0NzczPDzM+Pj4zPaamhoqKysZHBykpaWFwcFBJiYmZrYXw2l8fJxUKsXY2Jgpp2QySXt7O6lUir1790bi9NL2di772zo275nk97v3Mzqp1FQs4ZWHV3LmihpOaK1i165dBXdKJBIsWbKEkZERM/WU/reXSqUYGhoy5dTX10dtbS1dXV2xOM2FEw8SisgpwOWqela4fAmAqn4xS9n7gfWqetcsx7oOuFVVb0pfH/eDhLt27WLFihWxfX5UWPUCu25WvcCuW5xepfAg4T3AcSJytIhUElxlHNCbSkQ6gEZgS9q6RhGpCt+3AKcCj2buGzeudsPLF6teYNfNqhfYdXPVy4kmLFWdFJEPArcBZcC1qvqIiHwWuFdVp5PJecAN+vzLphcBV4nIFEFCvDK995Yr+IluSg+rbla9wK6bq15OJBAAVd0IbMxYd1nG8uVZ9rsLWBVpcAVgcHCQww47LO4wCo5VL7DrZtUL7Lq56uVKE5Z5Wlpa4g4hEqx6gV03q15g181VL59AisTg4GDcIUSCVS+w62bVC+y6uerlE0iRmJiYiDuESLDqBXbdrHqBXTdXvXwCKRLt7e1xhxAJVr3ArptVL7Dr5qqXTyBFwtXx/PPFqhfYdbPqBXbdXPXyCaRI1NbWxh1CJFj1ArtuVr3ArpurXj6BFImysrK4Q4gEq15g182qF9h1c9XLJ5AiMTQ0FHcIkWDVC+y6WfUCu26uevkEUiRaW1vjDiESrHqBXTerXmDXzVUvn0CKRH9/NPNbx41VL7DrZtUL7Lq56uUTSJFwYdTjKLDqBXbdrHqBXTdXvXwCKRKuXoLmi1UvsOtm1Qvsurnq5RNIkejq6oo7hEiw6gV23ax6gV03V72cGY3XOrnM7lWKlILXnqExbt7WzaYd/SQnpqipWELnsU2sXbV0zmljS8FtMVj1Arturnr5BOIxzd27B7li004mU1Okwmbk/RNTbNzey+1P9HNp50pWL2+IN0iPp0TxTVhFYmRkJO4QIsFlrz1DY1yxaSdjk88lj2lSCmOTU1yxaSd7hsay7u+yWz5Y9QK7bq56+QRSJNra2uIOIRJc9rp5WzeTqak5y0ymprh5W3fWbS675YNVL7Dr5qqXTyBFoqenJ+4QIsFlr007+g+48sgkpUG5bLjslg9WvcCum6tePoEUCRGJO4RIcNkrOTH31cd85Vx2ywerXmDXzVUvZxKIiJwtIo+JyA4R+VSW7ReISI+IPBC+Lkzbdr6IPBG+zi9u5LnR1NQUdwiR4LJXTUVuf96zlXPZLR+seoFdN1e9nEggIlIGfBN4I/Bi4J0i8uIsRW9U1RPD19Xhvk3APwN/C6wG/llEGosUes64egmaLy57dR7bRNk8P9zKJCiXDZfd8sGqF9h1c9XLiQRC8MW/Q1WfVNVx4AbgnBz3PQv4lar2q+oA8Cvg7IjiXDT19fVxhxAJLnutXbWU8rK5/8TLy5awdtXSrNtcdssHq15g181VL1cSyBHA7rTlp8N1mawVkYdE5CYRWb7AfWMllUrFHUIkuOy1rL6KSztXUlW+5IArkTKBqvIlXNq5ctaHCV12ywerXmDXzVWvUnqQ8OfAD1V1TEQuAq4HXp/rzt3d3axbt47y8nJSqRRr1qxh/fr1JBIJamtrKSsrY2hoiNbWVvr7+1FVWltb6erqmnkKdGRkhLa2Nnp6ehARmpqa6Onpob6+nlQqxb59+2hvbyeRSFBRUUFDQwO9vb00NDTQ3d39vO2VlZXU1dXR19dHY2MjyWSS0dHRme3V1dXU1NQwMDBAc3Mzw8PDjI+Pz2yvqamhsrKSwcFBWlpaGBwcZGJiYmZ7MZzGx8dJJBLU1dU567RMhH954wpuvH8Pf0qMMzqpVJUJnS84jFe1TnFE2T727St/nlMymaS9vZ09e/YwMTHhnNNi6mnaKZFIMDIyQnV1tSmn6b+9PXv2UFlZacqpr6+PsbExUqlULE5zIS6M8igipwCXq+pZ4fIlAKr6xVnKlwH9qtogIu8ETlfVi8JtVwGbVfWH6fts2bJFOzo6otSYk7GxMaqqZh82o1Sx6gV23ax6gV23OL22bt16X2dn58nZtrnShHUPcJyIHC0ilcB5wIb0AiJyeNriW4A/h+9vA94gIo3hzfM3hOucIpFIxB1CJFj1ArtuVr3ArpurXk40YanqpIh8kOCLvwy4VlUfEZHPAveq6gbgQyLyFmAS6AcuCPftF5ErCJIQwGdV1bnZVyoqKuIOIRKseoFdN6teYNfNVS8nEgiAqm4ENmasuyzt/SXAJbPsey1wbaQB5klDg80B+6x6gV03q15g181VL1easMzT29sbdwiRYNUL7LpZ9QK7bq56+QRSJFz9BZEvVr3ArptVL7Dr5qqXTyBFYnx8PO4QIsGqF9h1s+oFdt1c9fIJpEgkk8m4Q4gEq15g182qF9h1c9XLJ5Ai0d7eHncIkWDVC+y6WfUCu26uevkEUiRc7cedL1a9wK6bVS+w6+aql08gRaKysjLuECLBqhfYdbPqBXbdXPVy5jkQF9kzNMbN27rZtKOf5MQUNRVL6Dy2ibWrls46AN9s1NXVRRRlvFj1ArtuVr3ArpurXv4KZBbu3j3IRbdsZ+P2XvZPTKHA/okpNm7v5aJbtnP37sEFHa+vry+aQGPGqhfYdbPqBXbdXPXyCSQLe4bGuGLTTsYmpw6YUzulMDY5xRWbdrJnaCznYzY2OjfHVUGw6gV23ax6gV03V718AsnCzdu6mUzNPZ/2ZGqKm7d153xMV7vh5YtVL7DrZtUL7Lq56uUTSBY27eg/4Mojk5QG5XJldHQ0z6jcxKoX2HWz6gV23Vz18gkkC8mJua8+FloO3O3HnS9WvcCum1UvsOvmqpdPIFmoqcjttORaDtztx50vVr3ArptVL7Dr5qqXTyBZ6Dy26YA5tDMpk6BcrlRXV+cZlZtY9QK7bla9wK6bq14+gWRh7aqllJfNfWrKy5awdtXSnI9ZU1OTb1hOYtUL7LpZ9QK7bq56+QSShWX1VVzauZKq8iUHXImUCVSVL+HSzpULephwYGCgwFG6gVUvsOtm1QvsurnqNeeT6CLy2xyPM6qqbyhAPM6wenkDV63pKNiT6M3NzRFFGi9WvcCum1UvsOvmqtd8Q5m8AnjfPGUE+FphwnGLZfVVXHzqci4+dXnexxoeHubQQw8tQFRuYdUL7LpZ9QK7bq56zZdA7lLV6+c7iIj8fb6BiMjZBImoDLhaVa/M2P5R4EJgEugB/kFVd4XbUsC2sOhfVfUt+cZTaFydECZfrHqBXTerXmDXzVWvOROIqnbmcpB8m69EpAz4JnAm8DRwj4hsUNVH04rdD5ysqvtF5P3A/wHeEW5LquqJ+cQQNa72484Xq15g182qF9h1c9XLlZvoq4Edqvqkqo4DNwDnpBdQ1d+o6v5w8Y/AkUWOMS9c7cedL1a9wK6bVS+w6+aqV94JRET+qwBxHAHsTlt+Olw3G+uAX6QtV4vIvSLyRxE5twDxFBxXu+Hli1UvsOtm1QvsurnqVYj5QH5fgGPkjIi8CzgZeG3a6hWq+oyIHAPcISLbVPUv6ft1d3ezbt06ysvLSaVSrFmzhvXr15NIJKitraWsrIyhoSFaW1vp7+9HVWltbaWrq2vm5tXIyAhtbW309PQgIjQ1NdHT00N9fT2pVIp9+/bR3t5OIpGgoqKChoYGent7aWhoYP/+/ezatWtme2VlJXV1dfT19dHY2EgymWR0dHRme3V1NTU1NQwMDNDc3Mzw8DDj4+Mz22tqaqisrGRwcJCWlhYGBweZmJiY2V4Mp/HxcQYGBmhoaDDllEwmaW9vZ2BggFQqZcopkUiQSqU45JBDTDlN/+0NDAxQXV1tyqmvr4+Kigq6urpicZrz+1h1nlEDi4CInAJcrqpnhcuXAKjqFzPKnQF8HXitqmYdCldErgNuVdWb0tdv2bJFOzo6Iog+N3bt2sWKFSti+/yosOoFdt2seoFdtzi9tm7del9nZ+fJ2bbNewUiIq/PsnoC2KWqf803uJB7gONE5GjgGeA84Hk9u0TkJOAq4Oz05CEijcB+VR0TkRbgVIIb7E7R0tISdwiRYNUL7LpZ9QK7bq565dKEdU2WdRXAUhG5B3i7qj6TTxCqOikiHwRuI+jGe62qPiIinwXuVdUNwJeBQ4Efiwg81133RcBVIjJFcE/nyozeW04wODhIbW1t3GEUHKteYNfNqhfYdXPVa94EoqpHZ1svIocAVwJfBd6ebyCquhHYmLHusrT3Z8yy313Aqnw/P2omJibiDiESrHqBXTerXmDXzVWvRd9ED5/HuATYUcB4zOJqP+58seoFdt2seoFdN1e98u3GO0lhenKZx9V+3Pli1Qvsuln1Arturnrlm0A+AtxXgDjM42L7ZSGw6gV23ax6gV03V71y6YX1OyCzr28FcBQwCvy3COIyR1lZWdwhRIJVL7DrZtUL7Lq56pVL89PVWdZNAn8F/hQOPeKZh6GhIRobG+MOo+BY9QK7bla9wK6bq1659MKadzRez/y0trbGHUIkWPUCu25WvcCum6tehRgL652FCMQ6/f39cYcQCVa9wK6bVS+w6+aqVyFG4/1MAY5hHheGjIkCq15g182qF9h1c9Ur7wSiqn9TiECs4+olaL5Y9QK7bla9wK6bq16LTiAi0iQi60Xk7kIGZJWurq64Q4gEq15g182qF9h1c9VrQQ8Bikg5Qbfd84E3EczbcVUEcZnDxfmMC4FVL7DrZtUL7Lq56pVTAhGRlxMkjXeG+9xC8AzIq2YbVt3j8Xg8tpm3CUtEHgZ+BywDLgLaVXUdkIw4NlOMjIzEHUIkWPUCu25WvcCum6teudwDOQRIESSM/YB/cHARtLW1xR1CJFj1ArtuVr3ArpurXvMmEFU9huC+xzhwI5AQka8D1Rw4xIlnFnp6euIOIRKseoFdN6teYNfNVa+cemGp6m/DZqt24H8BxwN1wGYR+UCE8ZkhnATLHFa9wK6bVS+w6+aq14K68apqUlW/p6pvAFYA3wM+GElkxmhqaoo7hEiw6gV23ax6gV03V70W/RyIqj6jql9U1RcXMiCruHoJmi9WvcCum1UvsOvmqtecCUREchpIUUT+szDh2KW+vj7uECLBqhfYdbPqBXbdXPWa7zmQt4rIdcB8DXB/B7w3n0BE5Gzga0AZcLWqXpmxvQr4DvByoA94h6ruDLddAqwj6C32IVW9LZ9YoiCVSsUdQiRY9QK7bla9wK6bq17zJZAe4NocjpPXfIsiUgZ8EziT4On2e0Rkg6o+mlZsHTCgqseKyHnAl4B3iMiLgfOAEwieVfm1iLxQVZ064/v27aOlpSXuMAqOVS+w62bVC+y6ueo1ZwJR1ZVFimM1sENVnwQQkRuAc4D0BHIOcHn4/ibgGxJ0TTgHuEFVx4CnRGRHeLwtRYo9J9rb2+MOIRKseoFdN6teYNfNVa9CDOdeCI4AdqctPx2uy1pGVSeBQaA5x31jJ5HI6yLNWax6gV03q14Qj9ueoTG+/ofdnHv9g5x19f2ce/2DfP0Pu9kzNFawz3C1zhY0mGIp093dzbp16ygvLyeVSrFmzRrWr19PIpGgtraWsrIyhoaGaG1tpb+/H1WltbWVrq6umYHMRkZGaGtro6enBxGhqamJnp4e6uvrSaVS7Nu3j/b2dhKJBBUVFTQ0NNDb20tDQwNjY2Ps2rVrZntlZSV1dXX09fXR2NhIMplkdHR0Znt1dTU1NTUMDAzQ3NzM8PAw4+PjM9tramqorKxkcHCQlpYWBgcHmZiYmNleDKfx8XGGh4cZGxsz5ZRMJmlvb2d4eJi9e/eackokEiSTSUZGRkw5Tf/tDQ8PMzQ0VDSnv47X8MU7d5NSSIWPVe+fmGLj9l5ue7yXT77mCI6qSObl1NfXRyqVoqurK5Z6mgtxYaISETkFuFxVzwqXLwFQ1S+mlbktLLMlHBU4AbQCn0ovm14u/TO2bNmiHR0dxdDJyr59+6itrY3t86PCqhfYdbPqBcV12zM0xkW3bGdscmrWMlXlS7hqTQfL6qvy+qw462zr1q33dXZ2npxtW05NWCKyREReLyKVhQ1thnuA40Tk6PAzzgM2ZJTZQDAiMMBbgTs0yH4bgPNEpEpEjgaOA5ybo6S3tzfuECLBqhfYdbPqBcV1u3lbN5Op2ZMHwGRqipu35T9guat1lutQJlPAz1Q1koEUw3saHwRuA/4M/EhVHxGRz4rIW8Ji1wDN4U3yj/LclccjwI8Ibrj/EljvWg8sgIaGhrhDiASrXmDXzaoXFNdt047+mWar2UhpUC5fXK2zhdwD+a2IvFJV/xhFIKq6EdiYse6ytPejwNtm2ffzwOejiKtQjI/bHMTYqhfYdbPqBcV1S07MffWx0HJz4WqdLSSB7AJ+ISI/I+j1NJN707/oPdlJJm1On2LVC+y6WfWC4rrVVCxhfw7JoaYi/86urtbZQsxqgJ8SJI4jgeVpL888uNqPO1+seoFdN6teUFy3zmObKJtnjI4yCcrli6t1lnMCUdX3zvaKMkAruNqPO1+seoFdN6teUFy3tauWUl4291doedkS1q5amvdnuVpnC7q2EpEOEblURL4RLh8vIi+JJjRbVFZG1YEtXqx6gV03q15QXLdl9VVc2rmSqvIlB1yJlEnQhffSzpV5d+EFd+ss5wQiIm8jmBv9COA94eo64CsRxGWOurq6uEOIBKteYNfNqhcU32318gauWtPBmzpaOKRiCQIcUrGEN3W0cNWaDlYvL0zvKVfrbCE30T8LnKGqD4rIO8J1DwIvLXxY9ujr68vpyc5Sw6oX2HWz6gXxuC2rr+LiU5dz8anR3Q52tc4W0oS1FHgofK9p/8b/KHsJ0NjYGHcIkWDVC+y6WfUCu26uei0kgdwHvDtj3Xk4+NS3i7jaDS9frHqBXTerXmDXzVWvhTRhfQi4XUTWAbXhmFMvBN4QSWTGGB0djTuESLDqBXbdrHqBXTdXvXJOIKq6XUQ6gDcDtxI8THirqo5EFZwlXO3HnS9WvcCum1UvsOvmqtdCemG9RFX3q+qPVPXLqnqDTx6542o/7nyx6gV23ax6gV03V70W0oR1q4jUEnTlvTN83a8ujAdfAlRXV8cdQiRY9QK7bla9wK6bq14LeRL9KOAVBMOZvAT4MTAgIrdGE5otampq4g4hEqx6gV03q15g181VrwU9iR7OWX4XwXzjfwRSBN17PfMwMDAQdwiRYNUL7LpZ9QK7bq56LeQeyI0i8lfgO8AxwPeBlaq6OqrgLNHc3Bx3CJFg1Qvsuln1Arturnot5ArkZcAUwdPnDwIPqOpwJFEZZHjY5qmy6gV23ax6gV03V70Wcg/kOOAU4A7g1QRzgzwuIldHFZwlXJ0QJl+seoFdN6teYNfNVa+F3gPZCzwG7AB2Au3AGwsflj1c7cedL1a9wK6bVS+w6+aq10LugWwQkX7gZwTNWT8HXq6qR0QVnCVc7cedL1a9wK6bVS+w6+aq10KuQG4hSBgrVPXdqnq1qj6RbwAi0iQivxKRJ8J/Dxg1TEROFJEtIvKIiDyUNhowInKdiDwlIg+ErxPzjSkKXO2Gly9WvcCum1UvsOvmqtdChjK5TkTKReQ0gjlBngHuUtXJPGP4FLBJVa8UkU+Fy5/MKLMfeI+qPiEiy4D7ROQ2VX023P5xVb0pzzgixdUJYfLFqhfYdZvPa8/QGDdv62bTjn6SE1PUVCyh89gm1q5aWpDJkaLkYK2zuFhIE1YH8GfgBwQDK/4A2C4iL8ozhnOA68P31wPnZhZQ1cenr3ZUdQ/QDbTm+blFZXBwMO4QIsGqF9h1m8vr7t2DXHTLdjZu72X/xBQK7J+YYuP2Xi66ZTt373b7nByMdRYnC2nC+jfgP4DlqnqKqh4JfCtcnw9t4c15gATQNldhEVkNVAJ/SVv9+bBp66si4uRPpJaWlrhDiASrXmDXbTavPUNjXLFpJ2OTU6QyBihKKYxNTnHFpp3sGRorQpSL42Crs7hZyFhYJwJnZox99S/AZ+bbUUR+TdBjK5Pn7auqKiKzjq0lIocD3wXOV9WpcPUlBImnkiDBfZJg9sTn0d3dzbp16ygvLyeVSrFmzRrWr19PIpGgtraWsrIyhoaGaG1tpb+/H1WltbWVrq6umZnARkZGaGtro6enBxGhqamJnp4e6uvrSaVS7Nu3j/b2dhKJBBUVFTQ0NNDb20tDQwPd3d1UVVXNbK+srKSuro6+vj4aGxtJJpOMjo7ObK+urqampoaBgQGam5sZHh5mfHx8ZntNTQ2VlZUMDg7S0tLC4OAgExMTM9uL4TQ+Pk53dzfHHXecKadkMkl7eztPPfUUzc3NppwSiQTJZJIjjzzyAKdbdimTqanM/zrPYzI1xffu3sU/nrzUKafpv72dO3dy1FFHmain9P9PqVSK2traWJzmQnIdC1FEHgY+pKp3pK17HfANVT0hp4NkP+5jwOmqujdMEJtV9fgs5eqBzcAXZrvfISKnAx9T1TdnbtuyZYt2dHQsNsy82bVrFytWrIjt86PCqhfYdZvN69zrH2T/xNwJBII5v396vpszWR9sdVYMtm7del9nZ+fJ2bYt5Ark08CGcPDEXcAK4L8B78ozvg3A+cCV4b8/yywgIpXAT4DvZCYPETk8TD5CcP/k4TzjiQRX+3Hni1UvsOs2m1cyh+SxkHJxcLDVWdws5En0DcBJBF/QdeG/L1fVA77wF8iVwJki8gRwRriMiJyc9pT724HTgAuydNf9vohsA7YBLcDn8ownElztx50vVr3ArttsXjUVuX0d5FouDg62Ooubea9AROQQ4J+AvwG2Al9U1YLdRVPVPqAzy/p7gQvD998DvjfL/q8vVCxRUltbG3cIkWDVC+y6zebVeWwTG7f3HnADPZ0yCcq5ysFWZ3GTy0+JbwL/H7AdeCvw/0cakVHKysriDiESrHqBXbfZvNauWkp52dxfCeVlS1i7yt0ZHA62OoubXBLI2cAbVPUTBONeHXCD2jM/Q0NDcYcQCVa9wK7bbF7L6qu4tHMlVeVLKJPnbysTqCpfwqWdK51+mPBgq7O4yeUmeu30cxqqultEGiKOySStrSX13GPOWPUCu25zea1e3sBVazpK9kn0g7HO4iSXBFIedteVWZZJ79rryU5/fz+HHHJI3GEUHKteYNdtPq9l9VVcfOpyLj51eRGjKgwHa53FRS4JpBu4Nm25L2NZCWYo9MxBrs/blBpWvcCum1UvsOvmqte8CURVVxYhDvO4egmaL1a9wK6bVS+w6+aql7sduo3R1dUVdwiRYNUL7LpZ9QK7bq56+QRSJHIZV6YUseoFdt2seoFdN1e9fALxeDwez6LwCaRIjIyMxB1CJFj1ArtuVr3ArpurXj6BFIm2tjmnOSlZrHqBXTerXmDXzVUvn0CKRE9PT9whRIJVL7DrZtUL7Lq56uUTSJEIRpu3h1UvsOtm1Qvsurnq5RNIkWhqcncE03yw6gV23ax6gV03V718AikSrl6C5otVL7DrZtUL7Lq56uUTSJGor6+PO4RIsOoFdt2seoFdN1e9fAIpEqlUKu4QIsGqF9h1s+oFdt1c9fIJpEjs27cv7hAiwaoX2HWz6gV23Vz18gmkSLS3t8cdQiRY9QK7bla9wK6bq16xJxARaRKRX4nIE+G/jbOUS4nIA+FrQ9r6o0XkTyKyQ0RuFJHK4kWfO4lEIu4QIsGqF9h1s+oFdt1c9Yo9gQCfAjap6nHApnA5G0lVPTF8vSVt/ZeAr6rqscAAsC7acBdHRUVF3CFEglUvsOtm1Qvsurnq5UICOQe4Pnx/PXBurjtK8HTN64GbFrN/MWlosDkTsFUvsOtm1Qvsurnq5UICaZuecx1IALMN+lItIveKyB9F5NxwXTPwrKpOhstPA0dEF+ri6e3tjTuESLDqBXbdrHqBXTdXvXKZ0jZvROTXQLa7QJ9JX1BVFZHZ5m5coarPiMgxwB0isg0YzDWG7u5u1q1bR3l5OalUijVr1rB+/XoSiQS1tbWUlZUxNDREa2sr/f39qCqtra10dXXNjMU/MjJCW1sbPT09iAhNTU309PRQX19PKpVi3759tLe3k0gkqKiooKGhgd7eXhoaGlBVdu3aNbO9srKSuro6+vr6aGxsJJlMMjo6OrO9urqampoaBgYGaG5uZnh4mPHx8ZntNTU1VFZWMjg4SEtLC4ODg0xMTMxsL4bT+Pg4o6OjjI2NmXJKJpO0t7czOjrK3r17TTklEglSqRQjIyOmnKb/9kZHRxkaGjLl1NfXR0VFBV1dXbE4zYXEPdeuiDwGnK6qe0XkcGCzqh4/zz7XAbcCNwM9QLuqTorIKcDlqnpW5j5btmzRjo6OwgvkSHd3N0uXLo3t86PCqhfYdbPqBXbd4vTaunXrfZ2dnSdn2+ZCE9YG4Pzw/fnAzzILiEijiFSF71uAU4FHNch+vwHeOtf+LpBMJuMOIRKseoFdN6teYNfNVS8XEsiVwJki8gRwRriMiJwsIleHZV4E3CsiDxIkjCtV9dFw2yeBj4rIDoJ7ItcUNfoccbUfd75Y9QK7bla9wK6bq16xJxBV7VPVTlU9TlXPUNX+cP29qnph+P4uVV2lqi8N/70mbf8nVXW1qh6rqm9T1bG4XObC1X7c+WLVC+y6WfUCu26uesWeQA4WKiudfL4xb6x6gV03q15g181VL59AikRdXV3cIUSCVS+w62bVC+y6uepVlG68Hujr68upW1ypke61Z2iMm7d1s2lHP8mJKWoqltB5bBNrVy1lWX1VzJEunIOhzqxh1c1VL59AikRjY9Yhvkqeaa+7dw9yxaadTKamSIU9w/dPTLFxey+3P9HPpZ0rWb3czadpZ8N6nVnEqpurXr4Jq0i42g0vX5LJJHuGxrhi007GJp9LHtOkFMYmp7hi0072DDnZv2FWLNeZVay6uerlE0iRGB0djTuESBgdHeXmbd1MpqbmLDeZmuLmbd1FiqowWK4zq1h1c9XLJ5Ai4Wo/7nxpb29n047+A648MkkpbNrRX5ygCoTlOrOKVTdXvXwCKRKu9uPOl0QiQXJi7quPaXIt5wqW68wqVt1c9fIJpEhUV1fHHUIkVFdXU1OR259RruVcwXKdWcWqm6tepfU/uoSpqamJO4RIqKmpofPYJspk7nJlAp3HNhUnqAJhuc6sYtXNVS+fQIrEwMBA3CFEwsDAAGtXLaW8bO4/pfKyJaxdVVqjpFquM6tYdXPVyyeQItHc3Bx3CJHQ3NzMsvoqLu1cSVX5kgOuRMoEqsqXcGnnypJ7mNBynVnFqpurXj6BFInh4eG4Q4iEaa/Vyxu4ak0Hb+po4ZCKJQhwSMUS3tTRwlVrOkruIUKwX2cWsermqpd/Er1IjI+Pxx1CJKR7Lauv4uJTl3PxqctjjKhwHAx1Zg2rbq56+SuQIuFqP+58seoFdt2seoFdN1e9fAIpEq72484Xq15g182qF9h1c9XLJ5Ai4Wo3vHyx6gV23ax6gV03V718AikSrk4Iky9WvcCum1UvsOvmqpdPIEVicHAw7hAiwaoX2HWz6gV23Vz1ir0Xlog0ATcCK4GdwNtVdSCjzOuAr6at6gDOU9Wfish1wGuB6TN8gao+EG3UC6elpaVgx3Jp4qZCermGVTerXmDXzVUvF65APgVsUtXjgE3h8vNQ1d+o6omqeiLwemA/cHtakY9Pb3cxeUDhfkHcvXuQi27ZzsbtveyfmEJ5buKmi27Zzt27i/tLxdVfRoXAqptVL7Dr5qqXCwnkHOD68P31wLnzlH8r8AtV3R9lUIVmYmIi72O4OHFTIbxcxaqbVS+w6+aqlwsJpE1V94bvE0DbPOXPA36Yse7zIvKQiHxVRJwcL6MQ/bhdnLjJ1f7phcCqm1UvsOvmqldR7oGIyK+BbGfgM+kLqqoiMuvURCJyOLAKuC1t9SUEiacS+A/gk8BnM/ft7u5m3bp1lJeXk0qlWLNmDevXryeRSFBbW0tZWRlDQ0O0trbS39+PqtLa2kpXV9fMZPYjIyO0tbXR09ODiNDU1ERPTw/19fWkUin27dtHe3s7iUSCiooKGhoa6O3tpaGhgT179lBbWzuzvbKykrq6Ovr6+mhsbCSZTDI6Ojqzvbq6mpqaGgYGBmhubmZ4eJhfPd6b08RNv3q8j3ceXxO50/j4OIlEguOPP37RTuPj4zPba2pqqKysZHBwkJaWFgYHB5mYmJjZXox6Gh8fJ5lM0t7ezmOPPUZra6spp0QiwcjICCtWrDDlNP239+STT7Jy5UpTTn19fYyNjdHQ0BCL01yI6jzfSBEjIo8Bp6vq3jBBbFbV42cp+2HgBFX9x1m2nw58TFXfnLlty5Yt2tHRUbjAF0hvb2/eN8LOuvp+cqktAW678KS8PitXCuHlKlbdrHqBXbc4vbZu3XpfZ2fnydm2udCEtQE4P3x/PvCzOcq+k4zmqzDpICJCcP/k4cKHmD9lZWV5H8PFiZsK4eUqVt2seoFdN1e9XEggVwJnisgTwBnhMiJysohcPV1IRFYCy4E7M/b/vohsA7YBLcDnihH0QhkaGsr7GC5O3FQIL1ex6mbVC+y6ueoV+3MgqtoHdGZZfy9wYdryTuCILOVeH2V8haK1tTXvY6xdtZTbn+gnNTn7jfRiT9xUCC9Xsepm1Qvsurnq5cIVyEFBf39/3sdwceKmQni5ilU3q15g181Vr9ivQA4WCtVZYXriJleeRI+7E0aUWHWz6gV23Vz18gmkSBTyEtSliZtcvbQuBFbdrHqBXTdXvXwTVpHo6uqKO4RIsOoFdt2seoFdN1e9fAIpErk8lFOKWPUCu25WvcCum6tevgnL4/GUBC6NQu0J8FcgRWJkZCTuECLBqhfYdStFr1xHoS5Ft1xw1csnkCLR1jbfGJGliVUvsOtWal4LGYW61NxyxVUvn0CKRE9PT9whRIJVL7DrVmpeCxmFutTccsVVL59AikQwVJc9rHqBXbdS89q0oz+nUag37egvObdccdXLJ5Ai0dRUvPGpiolVL7DrVmpeyYm5rz7Sy5WaW6646uUTSJFw9RI0X6x6gV23UvNayCjUpeaWK656+QRSJOrr6+MOIRKseoFdt1LzWsgo1KXmliuuevkEUiRSqVTcIUSCVS+w61ZqXmtXLaW8bO6vqulRqEvNLVdc9fIJpEjs27cv7hAiwaoX2HUrNa+FjEJdam654qqXfxK9SLS3Z5sSvvSx6gV23UrRK9dRqEvRLRdc9fIJpEgkEglWrFgRdxgFx6oX2HUrVa9cRqEuVbf5cNXLN2EViZ/+9KdxhxAJVr3ArptVL7Dr5qqXTyBF4pZbbok7hEiw6gV23ax6gV03V718AikSk5OTcYcQCVa9wK6bVS+w6+aql7g6VWKh2bRpUw+wK67P7+/vb2lqauqN6/OjwqoX2HWz6gV23WL2WtHZ2Zl1SsSDJoF4PB6Pp7D4JiyPx+PxLAqfQDwej8ezKHwCKTAicraIPCYiO0TkU1m2XyAiPSLyQPi6MI44F4qIXCsi3SLy8CzbRUT+NfR+SEReVuwYF0MOXqeLyGBafV1W7BgXg4gsF5HfiMijIvKIiHw4S5lSrbNc3Equ3kSkWkTuFpEHQ6//naVMlYjcGNbZn0RkZQyhPoeq+leBXkAZ8BfgGKASeBB4cUaZC4BvxB3rItxOA14GPDzL9jcBvwAEeCXwp7hjLpDX6cCtcce5CK/DgZeF7+uAx7P8LZZqneXiVnL1FtbDoeH7CuBPwCszynwA+Fb4/jzgxjhj9lcghWU1sENVn1TVceAG4JyYYyoIqvpboH+OIucA39GAPwKHicjhxYlu8eTgVZKo6l5V3Rq+Hwb+DByRUaxU6ywXt5IjrIfpyc8rwldmL6dzgOvD9zcBnRLjbFM+gRSWI4DdactPk/0Pe23YZHCTiMw+LkNpkat7KXJK2KzwCxE5Ie5gFkrYzHESwS/adEq+zuZwgxKsNxEpE5EHgG7gV6o6a52p6iQwCDQXNcg0fAIpPj8HVqrqS4Bf8dyvCY+bbAVWqOpLga8DP403nIUhIocCNwMfUdWhuOMpJPO4lWS9qWpKVU8EjgRWi8jfxBzSnPgEUlieAdKvKI4M182gqn2qOhYuXg28vEixRc287qWIqg5NNyuo6kagQkRaYg4rJ0SkguAL9vuqmm0sjJKts/ncSrneAFT1WeA3wNkZm2bqTETKgQagr6jBpeETSGG5BzhORI4WkUqCm1wb0gtktDG/haD91gIbgPeEPXteCQyq6t64g8oXEWmfbmMWkdUE/2di+w+bK2HM1wB/VtWvzFKsJOssF7dSrDcRaRWRw8L3NcCZwPaMYhuA88P3bwXu0PCOehz44dwLiKpOisgHgdsIemRdq6qPiMhngXtVdQPwIRF5CzBJcPP2gtgCXgAi8kOCni0tIvI08M8EN/lQ1W8BGwl69ewA9gPvjSfShZGD11uB94vIJJAEzovzP+wCOBV4N7AtbFMH+DRwFJR2nZGbWynW2+HA9SJSRpDwfqSqt2Z8f1wDfFdEdhB8f5wXX7h+KBOPx+PxLBLfhOXxeDyeReETiMfj8XgWhU8gHo/H41kUPoF4PB6PZ1H4BOLxeDyeReETiMczCxKMnPz7uOPIFxHZLIsY9Tkc0XZKREZEJPOBtnziuU5EPhe+f2F4/NRiYvTEi38OxOMpACJyOXCsqr4r7lgKzB5VPTKqg6vq48ChIrI5qs/wRIe/AvF4PIsmHE7Dc5DiE4gnNkTkvSLy87TlJ0Tkx2nLu0XkxPD918LlIRG5T0ReE65fJiJJEWlK2+8kEekNx0tCRP5BRP4sIgMicpuIrEgrqyLyIRF5MtznyyKS9f/FHDGcTfAk9DvC5pgHw/UNInKNiOwVkWdE5HPhU8bZjr1aRLaIyLNh+W+Ew+Gkx/m+8Bw9KyLfTBuqo0xE/m8Y/1Mi8sGwfNYv97nOx3yEzXp/EJGvikgfcLmIvEBE7hCRvjCG708PyZFWH1tFZFhEbgSqc/08j9v4BOKJkzuB14jIEhFZRjAJ1ykAInIMcCjwUFj2HuBEoAn4AfBjEalW1T3AFmBt2nH/HrhJVSdE5ByCL/c1QCvwO+CHGXH8HXAywcRS5wD/MEu8s8XwS+ALBJP7HBqOAAtwHcGQNccSDDn+BmC2dv4U8D+BlvAcdBJMHpTOm4FXAC8B3g6cFa7/H8Abw9heBpw7y2eQ4/mYj78FngTagM8TTIT0RWAZ8CKCwf4uDz+vkmAk3O8SnLcf8/y68pQycc5m5V/+RTC3wcsIxvT5D+BuoINgXKYNc+w3ALw0fH8hwaByEHyZ7QZOC5d/AaxL228JwbhPK8JlBc5O2/4BYFP4/gLg9znGcDnwvbRtbcAYUJO27p3Ab3I8Lx8BfpK2rMCr05Z/BHwqfH8HcFHatjPC8uXh8mbgwlzOR0YMpwNPZ6y7APjrPLGfC9wfvj8N2EM4bFK47i7gcxn7zMToX6Xz8lcgnri5k+CL6rTw/WbgteHrzulCIvKxsNllUESeJRjGenp47psJJg86PDzOFMEva4AVwNfCZp9nCQagE54/cVL6pEq7CH5JH8A8MWSygmBQxr1pn30VsHSWY79QRG4VkYSIDBFc0WQeO5H2fj/BFRphvOkO6e+zxTXf+ZiP5x1fRNpE5IawmW4I+F5a7MuAZzTMEiG7FvBZHofxCcQTN9MJ5DXh+zvJSCDhvYZPEDTbNKrqYQQzsQmAqg4AtwPvIGi+uiHtC2s3wa/zw9JeNap6V1oM6XNiHEXwi/l5zBcDB049upvgCqQl7XPrVXW2mfH+nWDo7uNUtZ6gmSnXqUr3Eszlkc0nk1zOx3xkun4hXLcqjP1dPBf7XuCI6fs1IUct4LM8DuMTiCdu7gReR9DU8zTBlcPZBNN03h+WqSO4l9ADlIvIZUB9xnF+ALyHYBjvH6St/xZwiYRTmoY3tt+Wse/HRaRRgumFPwzcmCXO+WLoAlZO34DXYF6N24H/KyL14X2eF4jIa2c5D3XAEDAiIh3A+2cpl40fAR8WkSPCm9efnKNsLudjodQBI8CgiBwBfDxt2xaC8/YhEakQkTXA6jw/z+MIPoF4YkWD5wBGCJucNJia9EngD6qaCovdBvwSeJyg+WOUA5tpNgDHAQlVfTDt+D8BvgTcEDavPExwwzmdnwH3AQ8A/0Uw50Im88Uw3XusT0S2hu/fQ9Ax4FGC+yU3Ecz5kI2PEVw9DQPfJnsSm41vEySrhwiS7kaCL+1UZsEcz8dC+d8E97EGCc7fzAyBqjpOcMP+AoLmsnekb/eUNn4+EM9BjYgoQbPRjrhjKRQi8kbgW6qac/fcWY5zGkHiHAPeoaq3FSK+jM84jqB3WyXwAVW9rtCf4YkOn0A8BzUWEogE05++juAqpI2gU8EfVfUjccblsY9vwvJ4Sh8haEYaIGjC+jNwWawReQ4K/BWIx+PxeBaFvwLxeDwez6LwCcTj8Xg8i8InEI/H4/EsCp9APB6Px7MofALxeDwez6LwCcTj8Xg8i+L/AeU+ZFcS3mjpAAAAAElFTkSuQmCC","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["power_o_matrix = np.zeros_like(thetas)\n","power_e_matrix = np.zeros_like(thetas)\n","double_diff_matrix = np.zeros_like(thetas[1 : ])\n","\n","for i , angle in enumerate(thetas):\n","    # The polarization state analyzer for both beams of the Wollaston prism\n","    PSA_o = wollaston(0) @ linear_retarder(angle, np.pi)\n","    PSA_e = wollaston(1) @ linear_retarder(angle, np.pi)\n","    power_o_matrix[i] = (PSA_o[0, :] @ S_to_measure)\n","    power_e_matrix[i] = (PSA_e[0, :] @ S_to_measure)\n","    if i > 0 and i < (len(thetas)):\n","        double_diff_matrix[i - 1] = (power_o_matrix[i] - power_e_matrix[i]) - \\\n","            (power_o_matrix[i - 1] - power_e_matrix[i - 1])\n","\n","plt.style.use('bmh')\n","plt.figure()\n","plt.title('Double Difference Observed')\n","plt.plot(thetas[1 :], double_diff_matrix, marker = 'o', linestyle = 'None', \n","    markersize = 10)\n","plt.xlabel('waveplate angle [rad]')\n","plt.ylabel('Power [A.U.]')\n","plt.show()"]},{"cell_type":"code","execution_count":8,"metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["Stokes Vector Measured\n","[1.         0.32205896 0.6191101  0.        ]\n","Percent Difference\n","[-1.59943387e-06  2.25577236e-06]\n"]}],"source":["S_out = dual_channel_polarimeter(thetas, power_o = power_o_matrix, \n","    power_e = power_e_matrix, sub_method = \"double_difference\", normalized = False)\n","print('Stokes Vector Measured')\n","print(S_out)\n","print('Percent Difference (Q, U only)')\n","print(100 * (S_to_measure[1 : -1] - S_out[1 : -1]) / S_to_measure[1 : -1])"]}],"metadata":{"kernelspec":{"display_name":"base","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.9.7"}},"nbformat":4,"nbformat_minor":2}