{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## QuTiP Lecture: Single-photon Interference"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Authors: Amélie Orban (Editor), Inga Schöyen, Alessandro Delmonte, Alexander Rihane. \n",
" Based on the project Simulating Single Photon Interference using QuTiP (Quantum Toolbox in Python) carried out in June 2021 at Maastricht University, Maastricht Science Programme.\n",
" Last updated: 09/11/21."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import qutip\n",
"from numpy import pi, exp, cos, sin, sqrt, random\n",
"from qutip import Qobj, basis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Table of Contents\n",
"* [Introduction](#section0)\n",
"* [1. Theoretical Background](#section1)\n",
"* [2. Coding the constitutive elements](#section2)\n",
"* [3. Simulation of Single-Photon Interference Experiment](#section3)\n",
"* [4. Variations of Single-Photon Interference Experiment](#section4)\n",
"* [Conclusion and Takeaways](#section5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Introduction "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wave-particle duality is one of the key concepts of quantum mechanics and it illustrates the fact that every quantum particle or object can be described using both particle and wave properties. Neither of these classical concepts can fully describe quantum objects on its own.\n",
"\n",
"In a similar manner to Thomas Young's double-slit experiment, originally carried out in 1801, the wave-particle duality of light can be observed by performing a single-photon interference experiment, or in this case, simulation. The individual *particle*-like photons interfere with themselves, which is an intrinsically *wave*-like property, thus exhibiting both particle and wave characteristics at the same time."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This lecture investigates the phenomenon of single-photon interference by creating a simulation of a quantum optical experiment in which single photons go through a polarisation interferometer.\n",
"* As mentionned, the experiment to be simulated is similar to Young‘s, in that it also produces interference by allowing a photon to take two separate paths. In a classical sense, a beam of light is split into two separate waves, one of which will travel a different path length, to then be recombined into a single wave by using beam splitters and a polarization analyzer.\n",
"* Varying the path lengths results in a phase shift, which creates an interference pattern.\n",
"* Applying this classical understanding to the case under study, a single photon, which is an indivisible packet of light, would have to be ‘split’ into two ‘waves’ for an interference pattern to occur. Because the photons are sent through the path individually and cannot be split any further, it seems that they cannot interfere with each other. \n",
"* But they do interfere! The canonical understanding of the single-photon interference is that the photon’s probability wave interferes with itself, in opposition to interference happening between two distinct states."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1. Theoretical Background "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This section presents the setup of the single-photon interference experiment from a theoretical perspective, and details the different optical elements that must be simulated. The mathematical expression of each element will also be presented, in the form of Jones matrices acting as operators on the given quantum state.\n",
"\n",
"The setup and corresponding theory of this single photon interference experiment closely follows the work of Mark Beck in \"Quantum Mechanics, Theory and Experiment\" (*Beck, M. (2012). Quantum mechanics, theory and experiment. Oxford University Press.*), more specifically Experiment 6 of Section 3."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Half-wave plate $\\frac{\\lambda}{2}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A wave plate is used to modify the polarization of a wave, using the fact that its effective index of refraction depends on the polarization of the incident wave.\n",
"* A wave plate has two orthogonal axes. The fast axis represents the direction with the lower index of refraction as light polarized along that direction propagates faster, and the slow axis, orthogonal to the fast axis, represents the direction with the higher index of refraction.\n",
"* This difference in indices of refraction leads to the orthogonal components of the wave (polarization of a wave can always be decomposed into orthogonal components) acquiring different phase shifts when propagating through the material, resulting in an overall relative phase shift between the two components.\n",
"* A wave plate makes use of this relative phase shift between the components of the wave to modify its overall polarization.\n",
"\n",
"In this case, a half-wave plate is used, for which the relative phase shift between the fast and slow axes corresponds to a half wavelength shift.\n",
"* The fast axis of the half-wave plate is positioned in order to create a $45^\\circ$ angle with the horizontal plane. \n",
"* Half of the incident wave is thus polarized along the fast axis and result in a $|+45\\rangle$ polarization state, while the other half is polarized along the slow axis into a $|-45\\rangle$ state, with a relative phase shift of $\\pi$.\n",
"* The purpose of a half-wave plate is to rotate the linear polarization of some wave by an arbitrary angle (wave does *not* get split into its orthogonal components), depending on the orientation of the wave plate. \n",
"\n",
"The mathematical expression of a half-wave plate at $45^\\circ$ is the following:\n",
"\n",
"\\begin{equation}\n",
"J_{\\lambda/2 \\hspace{1mm} 45^\\circ} = \\left[\\begin{array}{cc} \\cos(2\\cdot45) & \\sin(2\\cdot45) \\\\ \\sin(2\\cdot45) & -\\cos(2\\cdot45) \\end{array}\\right] = \\left[\\begin{array}{cc} 0 & 1 \\\\ 1 & 0 \\end{array}\\right]\n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Polarization analyzers PA$_{HV}$ and PA$_{45}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A polarization analyzer, also called beam displacing polarizer or polarizing beam splitter, has the property that it can both *split* a beam into two orthogonal components, and *displace* each of these components differently. \n",
"* A general wave composed of orthogonal polarization components gets split into two waves travelling in different directions. When the faces of the polarization analyzer are parallel, the outgoing waves will emerge parallel, with a displacement separating them.\n",
"\n",
"If the polarization analyzer is positioned perpendicular to the horizontal plane, in such a way that waves propagating horizontally are normally incident, the general waves get split into vertical and horizontal components and it is referred to as a PA$_{HV}$. \n",
"* In this case, the $|V\\rangle$ polarization state gets transmitted without bending, while the $|H\\rangle$ gets bent when entering and leaving the PA$_{HV}$ such as to emerge parallel to $|V\\rangle$, but displaced and phase-shifted. \n",
"* This polarization analyzer can also be used to recombine orthogonal components of a beam. The process simply gets reversed. \n",
"* The Jones matrix representing a PA$_{HV}$ is composed of three other Jones matrices: the horizontal polarizer, vertical polarizer and phase-shifting matrices. This combination is made such that it represents the splitting and the phase-shifting caused by this element, and is as follows: \n",
"\n",
"\\begin{equation}\n",
"J_{PA_{HV}} = J_{V} + J_{\\phi}J_{H} = \\left[\\begin{array}{cc} 0 & 0 \\\\ 0 & 1 \\end{array}\\right] + \\left[\\begin{array}{cc} e^{i \\phi} & 0 \\\\ 0 & 1 \\end{array}\\right]\\left[\\begin{array}{cc} 1 & 0 \\\\ 0 & 0 \\end{array}\\right] = \\left[\\begin{array}{cc} e^{i \\phi} & 0 \\\\ 0 & 1 \\end{array}\\right]\n",
"\\end{equation}\n",
"\n",
"If one were to rotate a PA$_{HV}$ by $45^\\circ$, a general wave incident on such a polarization analyzer would be analyzed into $|+45\\rangle$ and $|-45\\rangle$ polarization components, making it a PA$_{45}$. \n",
"* In this case, the $|+45\\rangle$ polarization state gets transmitted without bending, while the $|-45\\rangle$ gets bent when entering and leaving the PA$_{45}$ such as to emerge displaced and parallel to $|+45\\rangle$. \n",
"* The Jones matrix for PA$_{45}$ is also a combination of other matrices: the linear polarizer with transmission axis at $45^\\circ$, linear polarizer with transmission axis at $-45^\\circ$ and phase-shifting matrices. The mathematical expression is thus: \n",
"\n",
"\\begin{equation}\n",
"\\begin{split}\n",
"J_{PA_{45}} &= J_{+45} + J_{\\phi}J_{-45} \\\\\n",
"&= \\left[\\begin{array}{cc} \\cos^2(45) & \\cos(45)\\sin(45) \\\\ \\cos(45)\\sin(45) & \\sin^2(45) \\end{array}\\right]\n",
"+ \\left[\\begin{array}{cc} e^{i \\phi} & 0 \\\\ 0 & 1 \\end{array}\\right]\\left[\\begin{array}{cc} \\cos^2(-45) & \\cos(-45)\\sin(-45) \\\\ \\cos(-45)\\sin(-45) & \\sin^2(-45) \\end{array}\\right] \\\\\n",
"&= \\left[\\begin{array}{cc} \\frac{1}{2}(e^{i \\phi} +1) & \\frac{1}{2}(1-e^{i \\phi}) \\\\ 0 & 1 \\end{array}\\right]\n",
"\\end{split}\n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Complete Setup"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The setup of the experiment, represented in the figure below, is as follows.\n",
"* A photon is prepared in the $|+45\\rangle$ state and is sent through an interferometer, composed of a PA$_{HV}$, a half-wave plate, another PA$_{HV}$, and finally a PA$_{45}$.\n",
"* The first polarization analyzer PA$_{HV}$ is used to split the incident photon into two orthogonal polarization components $|H\\rangle$ and $|V\\rangle$, and displace them accordingly.\n",
"* Each of the components emerges parallel to the other and then goes through the half-wave plate. Inserting the half-wave plate in the polarization interferometer allows for the equalization of the path lengths of the two arms.\n",
"* As already mentioned, a half-wave plate with its fast axis rotated to $45^\\circ$ causes half of the incident wave to be polarized into a $|+45\\rangle$ state (fast axis) and the other half into a $|-45\\rangle$ state (slow axis), with a phase shift of $\\pi$. In this specific setup, this implies that the polarization of the two components gets \"switched\" from $|H\\rangle$ to $|V\\rangle$ and vice-versa as the linear polarizations are rotated by $90^\\circ$, making the overall behaviour of the two arms symmetric with equalized path lengths and no relative phase shift (for now).\n",
"* The second polarization analyzer PA$_{HV}$ is then used to recombine the components.\n",
"* Finally, the use of a PA$_{45}$ at the end of the interferometer allows for the splitting of the final photon into $|+45\\rangle$ and $|-45\\rangle$ polarization states, providing two different paths and two different ports for the output. In this way, it is easy to measure how the intensity of each of the outputs gets modulated.\n",
"\n",
"This apparatus acts as an interferometer, and the final super-positioned states should interfere producing an interference pattern and output intensities that depend on the relative phase between them (equivalently, on the path length difference).\n",
"* This relative phase (initially zero) can be modified using the PA$_{HV}$'s.\n",
"* One (or both) of the two polarization analyzers can be tilted in order to vary the relative phase between the two arms, as the relative phase shift is proportional to the tilt angle of the PA$_{HV}$ (because it modifies the path lengths).\n",
"* When the phase is $\\phi = 0$, there should be constructive interference, with all of the output coming from the $+45^\\circ$ port, while when the phase is $\\phi = \\pi$, there should be destructive interference, with all of the output coming from the $-45^\\circ$ port.\n",
"* Note that this interference is visible because of the last element, the PA$_{45}$, as it projects the two, initially orthogonal (no possibility of displaying interference), polarization states onto the $+45^\\circ$ and $-45^\\circ$ axes where they do interfere."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that a single photon is sent through these optical elements, making the process quantum mechanical in nature. This requires for an important distinction, when it comes to the splitting of the state into orthogonal components.\n",
"* A classical wave is split into its orthogonal components deterministically, with the wave actually getting divided into two components, each *in proportion* to the coefficient of the corresponding polarization state in the general polarization ($\\psi = c_H |H\\rangle + c_V |V\\rangle$).\n",
"* On the other hand, when it comes to single photons, as they cannot be broken down any further, they are 'split' randomly and it is not possible to determine with certainty which port any photon will emerge from.\n",
"* The proportion of the basis states $|H\\rangle$ and $|V\\rangle$ in a general state $\\psi$ describes the *probability* of the photon taking that corresponding path.\n",
"* Note also that the state of the photon never actually gets 'split' and it stays as one single state throughout the entire interferometer."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. Coding the constitutive elements "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Basis vectors\n",
"Code the polarization basis vectors, for the bases $H/V$ and $+/-45$. They can then be used to code any polarization state, since any state can be expressed as a linear combination of basis vectors."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# HV basis\n",
"\n",
"# horizontal polarization\n",
"H = basis(2, 0)\n",
"# vertical polarization\n",
"V = basis(2, 1)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# +45/-45 basis (in terms of HV basis)\n",
"\n",
"# +45 polarization\n",
"p45 = 1 / sqrt(2) * (H + V)\n",
"# -45 polarization\n",
"n45 = 1 / sqrt(2) * (H - V)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Polarization analyzer (HV) n°1\n",
"It induces a phase shift on the H polarization relative to V polarization, because of the difference in the path lengths of the two arms of the interferometer. \n",
" In the case of a single photon, the state doesn't actually get split and there is only a single state throughout the apparatus, so the PA$_{HV}$ only affects the phase of the components."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False\\begin{equation*}\\left(\\begin{array}{*{11}c}(0.707+0.707j) & 0.0\\\\0.0 & 1.0\\\\\\end{array}\\right)\\end{equation*}"
],
"text/plain": [
"Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False\n",
"Qobj data =\n",
"[[0.70710678+0.70710678j 0. +0.j ]\n",
" [0. +0.j 1. +0.j ]]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Polarization analyzer (HV) n°1\n",
"\n",
"phaseshift1 = pi / 4 # CONSTANT\n",
" # should depend on real size of the setup (here: arbitrarily chosen)\n",
"PA_HV1 = Qobj([[exp(1j * phaseshift1), 0], [0, 1]])\n",
"PA_HV1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Half-wave plate $\\frac{\\lambda}{2}$\n",
"In this case, the fast axis makes an angle of $45^\\circ$ with the horizontal.\n",
" For a single photon, the effect of going through such a half-wave plate is that the $|H\\rangle$ and $|V\\rangle$ polarization components are switched."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\\begin{equation*}\\left(\\begin{array}{*{11}c}0.0 & 1.0\\\\1.0 & 0.0\\\\\\end{array}\\right)\\end{equation*}"
],
"text/plain": [
"Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\n",
"Qobj data =\n",
"[[0. 1.]\n",
" [1. 0.]]"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Half-wave plate\n",
"\n",
"θ = pi / 4 # fast axis orientation\n",
" # (!) numpy calculates with rad\n",
"halfwave = Qobj([[cos(2 * θ), sin(2 * θ)], [sin(2 * θ), -cos(2 * θ)]])\n",
"\n",
"halfwave.tidyup() # removes very small elements (numerical artifacts from the finite precision of the computer)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Polarization analyzer (HV) n°2\n",
"As before, it induces a phase shift on the H polarization relative to V polarization.\n",
"* Initially, the value of the phase shift induced by PA$_{HV2}$ should be set to the same value as for PA$_{HV1}$. The initial state of our setup has both of the PA$_{HV}$s in the exact same position (both tilted in the same way), with the overall relative phase shift being zero.\n",
"* The phase shift induced by PA$_{HV2}$ can then be changed in order to observe different levels of interference.\n",
"* A difference in the phase shifts induced by the two polarization analyzers causes the *relative* phase shift between the components to change, and it is this *relative* phase shift that is responsible for the interference."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False\\begin{equation*}\\left(\\begin{array}{*{11}c}(0.707+0.707j) & 0.0\\\\0.0 & 1.0\\\\\\end{array}\\right)\\end{equation*}"
],
"text/plain": [
"Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = False\n",
"Qobj data =\n",
"[[0.70710678+0.70710678j 0. +0.j ]\n",
" [0. +0.j 1. +0.j ]]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Polarization analyzer (HV) n°2\n",
"\n",
"phaseshift2 = pi / 4 # CHANGE TO CHANGE INTERFERENCE\n",
" # should depend on real size of the setup\n",
"PA_HV2 = Qobj([[exp(1j * phaseshift2), 0], [0, 1]])\n",
"\n",
"PA_HV2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Polarization analyzer (45)\n",
"The PA$_{45}$ acts on the photon such that it has to come out as a $|+45\\rangle$ polarized photon or as a $|-45\\rangle$ polarized photon.\n",
"* Based on the coefficients of each polarization basis components constituting the general polarization state, the probability of the photon getting out in that corresponding state can be computed.\n",
"* It should also induce a phase shift on $|-45\\rangle$ component with respect to $|+45\\rangle$ as a result of different path lengths (as seen in the theory section), but this effect is ignored as this is the last optical element, and the goal is to determine whether the final state of the photon is $|-45\\rangle$ or $|+45\\rangle$, not which one \"gets out first\".\n",
"* The role of this polarization analyzer is to determine the outcome of the experiment.\n",
"* The single photon interference happens just before the photon goes through this last polarization analyzer, e.g. when the states are 'recombined' after the second PA$_{HV}$ (this is not to be interpreted too literally because, as previously mentioned, the state of the photon never actually gets split and it stays as one state throughout the apparatus), and thus phase shift does not have an importance anymore."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# Polarization analyzer (45)\n",
"\n",
"# linear Polarizer, transmission axis +45 wrt horizontal\n",
"θ = pi / 4\n",
"Pp45 = Qobj([[cos(θ)**2, cos(θ) * sin(θ)], [cos(θ) * sin(θ), sin(θ)**2]])\n",
"# linear Polarizer, transmission axis -45 wrt horizontal\n",
"θ = -pi / 4\n",
"Pn45 = Qobj([[cos(θ)**2, cos(θ)*sin(θ)], [cos(θ)*sin(θ), sin(θ)**2]])\n",
"\n",
"def PA_45(vector):\n",
" p45_comp = Pp45 * vector # retrieve only +45 component\n",
" n45_comp = Pn45 * vector # retrieve only -45 component\n",
" return p45_comp, n45_comp"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3. Simulation of Single-Photon Interference Experiment "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Define initial variables:**\n",
"* One has to choose the initial minimal and maximal values of the phase shift induced by PA$_{HV2}$, along with the number of steps to be effectuated between these values. These different values will be run through during the simulation and the experiment will be carried out for each of them in order to collect data to determine what the relationship is between the interference and the relative phase shift. The higher the number of steps is, the more precise the final result will be.\n",
"* The number of iterations, for one value of the phase shift, also has to be defined. Since the photon getting out of the interferometer in a certain polarization state is governed by probabilities, it is important to carry out the same process numerous times to get a good sample and a good representation of the results."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"psi_0 = p45 # define the initial state (+45 vector)\n",
"\n",
"phaseshift2_init = pi / 4 # initial value\n",
"phaseshift2_max = 8 * pi\n",
"n = 100 # resolution of φ (amount of steps)\n",
"step = (phaseshift2_max - phaseshift2_init)/n # interval divided by number of small steps we want\n",
"N_init = 1000 # number of iterations (range(N) -> 0 to N-1, both included)\n",
"\n",
"# # create x- and y- coords. arrays to store the values needed to plot output graph\n",
"x_coords = [] # relative phase shift\n",
"y1_coords = [] # amount of photons in +45\n",
"y2_coords = [] # amount of photons in -45"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Create a *for loop* to run through all of the possible values for the phase shift considered for PA$_{HV2}$:**\n",
" In this loop, the passage of the photon through the interferometer is simulated using an effective matrix. The PA$_{45}$ is used after that to compute the probability related to each of these polarization states, and eventually to get the results of the experiment.\n",
"* 'Output' values are first defined to keep track of the results of the numerous iterations that are performed to measure the final state of the photon for a given phase shift. They are used to count the amount of photons that come out as $|+45\\rangle$ and the amount that come out as $|-45\\rangle$.\n",
"* The phase shift is also defined and directly added to the array of x-coordinates.\n",
"* With this phase shift, PA$_{HV2}$ can be defined, which allows for the definition of the final effective matrix. The effective matrix regroups all of the operators acting on the photon (in the correct order) in a single matrix.\n",
"* This is then multiplied with the initial state of the photon, to get the final state that will be passed through the PA$_{45}$.\n",
"* The $|+45\\rangle$ and $|-45\\rangle$ components are then retrieved and used to compute the final probabilities. To do so, the norm of the vector of each component is computed and then squared.\n",
"\n",
"Inside the first *for loop* is another loop performing what can be considered as the final measurement. Whether the photon comes out in a $|+45\\rangle$ or $|-45\\rangle$ polarization state is a random process that depends on the corresponding probabilities.\n",
"* To simulate this, a random number is generated between 1 and 100 (100 because $100\\%$ of the photons need to come out in either one of the states).\n",
"* If the number is smaller or equal to the probability (in percentage form) related to the $|+45\\rangle$ component, the photon is measured to have come out in the $|+45\\rangle$ state. Otherwise, it is measured in the $|-45\\rangle$ state.\n",
"* These measurements are then accounted for and stored in the 'output' values that were previously defined.\n",
"* When all of the iterations are performed, the output values are added to the arrays corresponding to different y-coordinates (one for $|+45\\rangle$ photons and one for $|-45\\rangle$ photons).\n",
"* They are then reinitialized at the beginning of the first *for loop*, before the entire process starts anew for another value of the phase shift."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"for i in range(n + 1):\n",
" output_p45 = 0\n",
" output_n45 = 0\n",
"\n",
" phaseshift2 = phaseshift2_init + i * step\n",
" x_coords.append((phaseshift2 - phaseshift1) / pi) # add realtive phase shift to x coords\n",
" PA_HV2 = Qobj([[exp(1j * phaseshift2), 0], [0, 1]]) # create corresponding PA_HV2\n",
" EffM = PA_HV2 * halfwave * PA_HV1 # define the effective matrix\n",
"\n",
" psi_final = EffM * psi_0 # apply the effective matrix to the initial state to get the final state\n",
"\n",
" psi_p45 = PA_45(psi_final)[0] # retrieve +45 and -45 components\n",
" psi_n45 = PA_45(psi_final)[1]\n",
"\n",
" proba_p45 = round(psi_p45.norm()**2, 5) # probability is rounded up to five decimals to avoid machine precision artifacts\n",
" proba_n45 = round(psi_n45.norm()**2, 5)\n",
"\n",
" for j in range(N_init):\n",
" a = random.randint(1, 100) # generates random number between 1 and 100 (both included), 100 because 100% of the photons need to come out in either +45 or -45 state\n",
" if a <= proba_p45 * 100:\n",
" output_p45 = output_p45 + 1\n",
" else:\n",
" output_n45 = output_n45 + 1\n",
"\n",
" y1_coords.append(output_p45)\n",
" y2_coords.append(output_n45)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Create output plot:**\n",
" Using the arrays (which are collections of coordinates) that are created as the simulation goes on for different values and setups, a plot can be created to visualize the dependency of the interference on the relative phase shift, or similarly, on the different in path lengths of the two arms of the interferometer."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABS0ElEQVR4nO2deZhUxdW438MMwwwwDCDLyDaLImowKBARFeEDjZioxCWCW1xAUb4vilvUX2I0MYkriWsUN4gkUTQxUVyDioDiAiguI4oIqMAMIEIQ2Znz+6OqZ3qa3qanu2811Ps89fTtunVvnVu3bp2qU5uoKh6Px+PxxKNZ0AJ4PB6Px328svB4PB5PQryy8Hg8Hk9CvLLweDweT0K8svB4PB5PQryy8Hg8Hk9CvLLwpIwYJonIOhF5J4XrXxORMZmQzTVE5AUROSeF6waJyKcuyWSvvUFEbkizSB6H8coig9jCcJ2ItAhalmiIyLki8noTbnEkcAzQTVUPTZNYSSEiy0Tk6GzG2RRU9ThV/UuicCKiIrJv2HWzVbVXkDKlGxHpIiLLkwj3ik2P/DC/ZSKyWUQ2WvefDMo5JBk5I65p8P52J7yyyBAiUg4MAhQ4MVhpMkYZsExVvwtaEI87JNHq+BHwYoJ7nAk0j3H6BFVtbd0Pw65pLyKxrvE0Ea8sMsfPgLeAyUCDpr6ITBaRP1szwEYReUNESkXkDtsS+UREDgkLf4BtpawXkSoROTHsXANTTmRrwdZ0LhKRz+z191rz0QHA/cBAK8P6aA9ha4HPiMg3IrJYRC6w/qOBh8Ku/02Ua8+1z3aPiPzXPtewiGBlNsy3IvIfEekQdv2J9nnX2+c8wPpPAXoA02zcv4gX3p5bJiJXisgHVpapIlJoz3UQkWftdd+IyGwRifptiMj+IjLdhvtURE6z/vtYv75h6bZGRIZEvicR2VdEZlo5vhaRqdZ/lo3mfftcIyNrt/Gew57/hYhUi8hKERkTr6YbIdO5IvK6iNxu8+BSETku2nVp4EfA87FOikgJcD3wi0be9xhguYhMEJHeyV4kIj8SkY9tHlxh07cV8ALQRepbMV1E5FARedPmlWqbtwvsfXZ5f9b/eBFZYK+ZIyLfD4v7ahvntzY/RX4f7qCq3mXAAYuBcUA/YDvQOezcZOBre64QeBVYilEwecDvgBk2bHN7r/8HFABDgW+BXvb8a8CYsHufC7we9l+BZ4G2mAJ2DTA8WtgYzzEL+LOV82B7/dBkrrfndwCX2ecYCfwXaB8m++fAfkCR/X+zPbcf8B2mAGiOKTgWAwX2/DLg6LC4kgn/DtAFaA8sBC6y527CKM7m1g0CJMrztAK+As4D8oFD7Hs80J6/APgYaAm8BNwedm3dewIeA36JqawVAkdGvK99w/4PAZaH/Y/3HMOBGuB7Voa/Rt4v4nnCZToXk08vwOTBi4GV0dLBhr8BuCGFc81tmhXHyTf32jxTbuXPj3j+VZh8+B+gT8S1vYHbrOxzMd9guwR5vBoYZI/bAX2jpb316wccZt9/uU3/8XHe3yHAamCATddz7DO0AHrZ/NTFhi0H9gm67IrlfMsiA4jIkRgTzROqOh9TIJ4REexfqjpfVbcA/wK2qOqjqroTmIrJZGAyZmtMIbpNVV/FFP6nN0Kkm1V1vap+CczAFPrJPEd34AjgalXdoqoLMK2JnzUi7tXAHaq6XVWnAp8CPw47P0lVF6nqZuCJMNlGAs+p6nRV3Q7cjlEoh8eIJ5nwd6nqSlX9BpgWFtd2YG+gzMo5W+3XG8HxGLPbJFXdoarvAf8Efgqgqg9iFNTb9n6/jCHrdkz+6GLTtbH9RrGe4zRMelap6iZMod0YvlDVB20e/It9hs6NvEcijgLeV9Vvo50Ukf6YPHd3jOvPxBSqZZi8/JKItA2dVNWPVPUqoDvm+YcAS0XkcRFpE+Oe24EDRaSNqq5T1XdjCW+/2bfs+18GTAQGxwoPXAhMVNW3VXWnmj6irZjveidGaRwoIs1VdZmqfh7nXoHilUVmOAf4j6p+bf//nQhTFKZ2FGJzlP+t7XEX4CtVrQ07/wXQtRHy1IQdbwq7dyK6AN9EfNiNjXtFRMH7hb1vItm62LAA2Of/Kk7cyYSPFddtmEL+PyKyRESuiRFHGTDAmhPWizHdnQmUhoV5EFO7vVtVt8a4zy8AAd6xZrPzY4SLRbw0+yrsXPhxo+5rlQ0kmVfCzHjrgWuAa8LS6dmwoHUmKBE5M8zE84I1/f0ZuFRVd0SLR1XfUNXNqrpJVW8C1mNagpHhdgIfAu8D32DeSaz+jFOsXF9Y8+DAOM+5n33WGhHZAPwB6BArPCbPXBGRZ7pjKgqLgfEYpbbaKrQuMe8UMF5ZpBkRKcLU8AbbDFWDaVL3EZE+KdxyJdA9wobeA1hhj7/DmBxChBdciUi05PBKoL2IFMeIOxm6iohEXL8yietWYj40wAzTxXxkobgjZU8UPiaq+q2qXqGqlZjBCJfHsB1/BcxU1bZhrrWqXmzjbA3cATwM3CAi7WPEV6OqF6hqF2As8OdY/QqNpBroFva/exrumRSqenwoTYCbMa3ZUBodHxa0Tlmo6t+0vqP6OKAN0B+Yar+bufaa5SKyi0IIRY1RvIB5B7b/5VXgXUxlYaSq9lbVtTFkn6uqI4BOwL8xLdzQvSO5D/gE6KmqbTDmYYkSLsRXwO8j8kxLVX3Mxv13VQ1ZIhS4Jc69AsUri/TzE0zz8kCMeeBg4ABgNo0z34R4G1N7/IWINBfTYXoC8Lg9vwA4WURa2gJndCPuvQroFuqgi0RVvwLmADeJSKHtmBuNsYUnSyfgEiv7TzFpEbNzM4wngB+LyDAxI1yuwDTf54TJXtmI8DGxHZD7WgXzX8z7q40S9FlgPxE52z5PcxH5gdR3pN8JzFPVMcBzmH6QaPH9VERChfo6TCERii/yuRrDE8B5YgZEtASuS/E+GUFEKoAWqrowRpD/YlpHB1v3I+vfD3hbRHqIyBEiUmDz41WYWv0b9v7DMZWGkRjzUFdVHaeqc4mBvdeZIlJizZcbaPgu9hLT4R6i2IbZKCL7Y/p2wol8fw8CF4nIADG0EpEfi0ixiPQSkaFihtZvwVgUouU7J/DKIv2cg7Ebf2lrkDWqWgPcA5wpYWPGk0FVt2GUw3GYjsE/Az9T1U9skD8B2zCZ9C/A3xpx+1eBKqBGRL6OEeZ0jI14JaZv5XpVfbkRcbwN9LSy/x44NVYNLxxV/RQ4C2O7/hqTBifY9ADTKf0r27S/Monw8egJvAxsBN4E/qyqM6LI9C3wQ2AUJj1qMDXBFiIyAtPBHCo8Lgf6ihkCGskPMIXfRuAZjNlliT13A/AX+1ynJSF7uHwvAHdhbPmLMaPxwChNF/gxcSoKagj/ZtbYU6vseyzG1OzXYVqMw4HjwvLTp8D+auaPTI1jBozkbGCZNStdhDEtYr+xx4Al9n10Aa7E9D9+i1EEUyPudQNh709V52EGDdxj5V6MGUwApr/iZkx+rcFUrK5NUuasI9H78TyepiMi52JG2xwZtCx7IrbF8xGmNh+1D6AJ974BQFVvaMQ1zwP3qGoyLUuPY/iWhcezGyEiJ4lICxFph2n1TEu3omgCr2FaPZ4cxCsLj2f3YixmuPLnmL6XSJt6unjNuqRR1VvVDJH25CDeDOXxeDyehPiWhcfj8XgS0qiROblChw4dtLy8PGgxPB6PJ6eYP3/+16raMdq53VJZlJeXM2/evKDF8Hg8npxCRL6Idc6boTwej8eTEK8sPB6Px5MQryw8Ho/HkxCvLDwej8eTEK8sPB6Px5MQryw8Ho/Hk5Ddcuisx1BbW8vXX3/N+vXr2blzZ9DieDxOkpeXR9u2benQoQPNmvn6cyy8stiNWb58OSJCeXk5zZs3p+EeRB6PR1XZvn07q1atYvny5fTo0SNokZzFq9HdmO+++46uXbtSUFDgFYXHEwURoaCggK5du/Ldd98FLY7TZExZiMgjIrJaRD6Kcu4KEVER6WD/i4jcJSKLReQDEekbFvYcEfnMush9rD0J8M1qjycx/jtJTCZTaDJmJ6sGiEh3zG5jX4Z5H4fZrawncCFmNyzsHsbXAwOAQ4Hr7Tr9Ho/H48kiGVMWqjoL+CbKqT8Bv6DhZugjgEfttopvAW1FZG/gWGC6qn6jquuA6URRQB6Px+PJLFlte9l9ileo6vsRp7oCX4X9X279YvlHu/eFIjJPROatWbMmWhDPbsSyZcsQEV5//fWgRUkrkydPJj/fjzvxuEfWlIWItAT+H/DrTNxfVR9Q1f6q2r9jx6gr7HpyiHPPPRcRQUTIz8+nrKyMiy66iLVr16Y1ntdffx0RYdmyZWm9b6qMHDmSFStWZDye/Px8Jk+e3OjrxowZw5AhQ9IuT2OpqamhtLQUEWH58uV1/qFKRKT71a9+FaC0uwfZbFnsA1QA74vIMqAb8K6IlAIrgO5hYbtZv1j+nixTUwNvvml+s8WgQYOorq5m2bJl3HXXXfzzn//kZz/7WfYECICioiI6d+4ctBiBMHny5KQUUW1tLWeeeSaHHnpozDBPP/001dXVde6aa65Jo6R7JllTFqr6oap2UtVyVS3HmJT6qmoN8AzwMzsq6jDgv6paDbwE/FBE2tmO7R9aP08WmTIFKith+HDzO2VKduItKCigtLSUbt26MWLECMaPH8+LL77I5s312zivXLmS448/npYtW1JZWblLbbm6uppRo0bRtm1bioqKGDJkSN1eJ8uWLWPQoEEAVFRUICJ1hZWqcvvtt1NZWUlBQQH77LMPd9xxR4N7l5eX8+tf/5pLL72U9u3b07lzZy677DJ27NhRF+b111/niCOOoLi4mOLiYvr06cNLL8XOwpFmqND/N954g759+9KyZUv69evH3Llz46ZdVVUVxx57LG3btqVVq1YccMABTLEvrry8nJ07d3LeeefV1bwB1q1bx1lnnUWPHj0oKiqiV69eTJgwgdDWyzfccAMPP/wwM2fOrLsulN4bN27k0ksvpWvXrrRs2ZJDDjmEp556Kq6MqXLjjTdSUFDAZZddFjNM+/btKS0trXOtW7fOiCx7FKqaEQc8BlQD2zGKYXTE+WVAB3sswL2YTeY/BPqHhTsfWGzdecnE3a9fP/Wofvzxx02+R3W1alGRKtS7oiLjn0nOOeccHTZsWAO/CRMmKKAbNmzQpUuXKqAVFRU6depU/eyzz/Taa6/VvLw8/fTTT1VVtba2Vg899FDt06ePzp49Wz/44AM97bTTtG3btrpmzRrdsWOHPv300wroO++8o9XV1bp27VpVVb3nnnu0sLBQJ06cqIsWLdL77rtPW7RooQ899FCdPGVlZdq2bVu96aabdNGiRTp16lTNz8+vC7N9+3Zt166dXnbZZbpo0SJdtGiRPvXUUzpr1qyYzz1p0iTNy8tr8F9EdNCgQTpr1ixduHChDh8+XMvLy3X79u0x73PQQQfp6aefrlVVVfr555/r888/r9OmTVNV1dWrV2teXp7ecccdWl1drdX2ZVZXV+tNN92k8+fP1yVLluiUKVO0VatW+sgjj6iq6rfffqtnnHGGDhw4sO66TZs2aW1trQ4ZMkQHDx6ss2fP1s8//1wnTpyozZs315dffjnpdz5p0iQdPHhw3DCvvvqqdu3aVVetWqUzZsxQQL/66qu686F80b17d91rr720X79+OmHCBN22bVvC+NPxveQ6wDyNVabHOpHLzisLQzoy/5w5qm3aNFQWJSXGP5NEKouqqiqtrKzUAQMGqGp9oTBhwoS6MDt27NDWrVvr/fffr6qqL7/8sgJaVVVVF2bLli1aWlqqv/nNb1RVdfbs2Qro0qVLG8TfrVs3veqqqxr4jR8/XisqKur+l5WV6QknnNAgzPDhw3XUqFGqqvrNN98ooDNmzEj6uaMpC0Dnz59f5/fWW28poJ988knM+7Rp00YnTZoU83xeXl7c8yEuueQSPfroo+v+jx49epcCfcaMGdqiRQtdv359A//zzjtPR4wYkTCOEImURU1NjXbp0kWnT59eF2+kslizZo3eeuut+sYbb+h7772nd955p7Zp00bPOuushPF7ZRFfWfhhF564VFTA9u0N/bZtM/6Z5rXXXqN169bs3LmTrVu3MmzYMCZOnNggzMEHH1x3nJeXR6dOnVi1ahVgTDF77bUXBx54YF2YFi1aMGDAAKqqqmLGu2HDBpYvX85RRx3VwH/w4MHceeedbNq0iZYtW+4SP0CXLl1YunQpAO3atWPMmDEce+yxDB06lMGDB3PSSSfRq1evRqWDiNCnT58GcQCsWrUq5r2uvPJKxowZU9cPcOKJJ9K3b9+oYUPU1tZy66238vjjj7N8+XK2bNnC9u3bKSsri3vd3Llz2bZtG127NhyouG3bNnr27BnzutmzZ3PcccfV/d+xYwfbt29vYDIaNGgQL7zwAgBnnnkmP/vZzzj66KNj3rNDhw5cddVVdf8PPvhgiouLOf/887n55pt3kdGTPH7aoicupaUwcSIUFUFJifmdONH4Z5oBAwawYMECFi5cyJYtW5g+fTqVlZUNwhQUFDT4LyLU1tZmXrgk43/wwQeZP38+xxxzDDNnzqR37967KLxENGvWjLy8vAZxAHGf87rrrmPRokWcdtppfPTRRxx22GEJRwRNmDCBm266iUsuuYTp06ezYMECxowZw7Zt2+JeV1tbS0lJCQsWLGjgPv7447qCPhr9+/dvEP63v/3tLn4PPfRQXfhXXnmF2267jfz8fPLz8xk2bBhg+mDGjh0bM57DDz8cgC++iLm9tCcJfMvCk5Czz4ZjjoGlS02LIhuKAszIoH333Tfl67/3ve+xdu1aPv7447rWxdatW3n77bcZN24cUF/Yh6/K26ZNG7p168asWbM4/vjj6/xnzpxJRUVFXasiWXr37k3v3r25/PLLueiii3jggQfiFm7porKyknHjxjFu3DhuvvlmbrvtNn73u98B5rkjVyKeNWsWw4cP5/zzz6/z++yzzxqEiXZd//79Wb9+PVu2bKF3795Jyxf5fjt16hT3nX/44YcN/s+dO5fzzz+fl156iQMOOCBmPO+++y4A3bp1S1o2z654ZeFJitLS7CmJdDF06FAOPfRQzjjjDO69915KSkq48cYb2bJlCxdffDEAZWVlNGvWjOeff56RI0fSokULSkpKuPbaa7niiivo2bMnQ4YM4dVXX+W+++7j3nvvTTr+xYsX8+CDD3LCCSfQvXt3Vq5cyezZsxOag5rKxo0bufrqqznllFOoqKhg/fr1vPjiiw3McRUVFcyYMYPjjjuOgoICOnToQK9evZgyZQozZsyga9euPProo7z99tu0a9euwXVPPvkkVVVVdO7cmeLiYoYOHcrRRx/NySefzK233sr3v/991q1bx5w5cygsLOSCCy5Iy3NFKqKvv/4agF69etWZ5iZPnkxeXh59+/alsLCQ2bNnc9VVV3Hqqaf6FWWbiDdDeXZbRIR///vf7L///vz4xz/mBz/4ATU1NUyfPp0OHToA0LlzZ2666SZuvvlm9t57b0aMGAHAxRdfzG9/+1v+8Ic/cOCBB3LLLbdw8803M3r06KTjb9WqFZ999hmjRo1iv/3245RTTuHwww/nnnvuycjzhsjPz2fdunWMHj2aAw44gGOPPZbOnTvz97//vS7MhAkTmD9/PuXl5YQmsV533XUMHjyYESNGMHDgQNatW8cll1zS4N6jR4/mBz/4AYcffjgdO3bkscceQ0R45plnOPnkk7nsssvq0vu5555jn332yeizRtKsWTNuvfVWBgwYwPe//31uv/12rrrqKv76179mVY7dEVHVxKFyjP79+2toLP2ezMKFC+M2zz0eTz3+ewERma+q/aOd8y0Lj8fj8STEKwuPx+PxJMQrC4/H4/EkxCsLj8fj8STEKwuPx+PxJMQrC4/H4/EkxCsLj8fj8STEKwuPx+PxJMQrC4/H4/EkxCsLT04S2mv59ddfD1qUtBK5U57H4wpeWXic5Nxzz63bujM/P5+ysjIuuugi1q5dm9Z4Xn/9dUSEZcuWpfW+qTJy5EhWrMj8NvP5+fm7bEGbDGPGjElqn+x0UFtby7BhwxCRXdZ2Ki8vr8sfIXfkkUdmRa49FV+F8SRHTU3W1ygfNGgQTzzxBDt27GD+/PmMGTOGr776iueeey4r8QdBUVERRUVFQYvhBL/97W9p1apVzPNXX30148ePr/sfubeIJ71krGUhIo+IyGoR+SjM7zYR+UREPhCRf4lI27Bz14rIYhH5VESODfMfbv0Wi8g1mZLXE4cpU6CyEoYPN79TpmQl2oKCAkpLS+nWrRsjRoxg/PjxvPjii2zevLkuzMqVKzn++ONp2bIllZWVu9SWq6urGTVqFG3btqWoqIghQ4YQWmRy2bJlDBo0CDBLb4tIXa1ZVbn99tuprKykoKCAffbZhzvuuKPBvcvLy/n1r3/NpZdeSvv27encuTOXXXYZO3bsqAvz+uuvc8QRR1BcXExxcTF9+vThpZdeivnMkWao0P833niDvn370rJlS/r168fcuXPjpl1VVRXHHnssbdu2pVWrVhxwwAFMse+tvLycnTt3ct5559XVygHWrVvHWWedRY8ePSgqKqJXr15MmDCB0GKjN9xwAw8//DAzZ86suy6U3hs3buTSSy+la9eutGzZkkMOOYSnnnoqrozxePXVV3nkkUeYNGlSzDCtW7emtLS0zrVv3z7l+DxJEGu/1aY64CigL/BRmN8PgXx7fAtwiz0+EHgfaAFUAJ8DedZ9DlQCBTbMgYni9ntwG9Kyp3B1tWpRUcNNuIuKjH8GidyDW1V1woQJCuiGDRvq9uCuqKjQqVOn6meffabXXnut5uXl6aeffqqqqrW1tXrooYdqnz59dPbs2frBBx/oaaedpm3bttU1a9bojh079Omnn1ZA33nnHa2urta1a9eqquo999yjhYWFOnHiRF20aJHed9992qJFC33ooYfq5CkrK9O2bdvqTTfdpIsWLdKpU6dqfn5+XZjt27dru3bt9LLLLtNFixbpokWL9KmnntJZs2bFfO5oe3CLiA4aNEhnzZqlCxcu1OHDh2t5eblu37495n0OOuggPf3007Wqqko///xzff7553XatGmqqrp69WrNy8vTO+64Q6urq7Xavsvq6mq96aabdP78+bpkyRKdMmWKtmrVSh955BFVVf3222/1jDPO0IEDB9Zdt2nTJq2trdUhQ4bo4MGDdfbs2fr555/rxIkTtXnz5vryyy8n/c5D1NTUaNeuXev2Lgd0ypQpDcKUlZVp586dtX379nrggQfqz3/+c/36668bHVc4fg/u+HtwxyvspwHPxHKxrou4R3m4sog4dxLwN3t8LXBt2LmXgIHWvRTm3yBcLOeVhSEtmX/OHNU2bRoqi5IS459BIpVFVVWVVlZW6oABA1RV65TFhAkT6sLs2LFDW7durffff7+qqr788ssKaFVVVV2YLVu2aGlpqf7mN79RVdXZs2croEuXLm0Qf7du3fSqq65q4Dd+/HitqKio+19WVqYnnHBCgzDDhw/XUaNGqarqN998o0BdoZcM0ZQFoPPnz6/ze+uttxTQTz75JOZ92rRpo5MmTYp5Pi8vL+75EJdccokeffTRdf9Hjx6tgwcPbhBmxowZ2qJFC12/fn0D//POO09HjBiRMI5wdu7cqcOGDdPrrruuzi+asrj99tt1+vTp+sEHH+gTTzyh++23n+633366adOmRsUXjlcW8ZVFvD6L2+3vyUApEOphOh1YFbe5khznA1PtcVfgrbBzy60fwFcR/gOi3UxELgQuBPyOWOmkogK2b2/ot22b8c8wr732Gq1bt2bnzp1s3bqVYcOG7bJ/9cEHH1x3nJeXR6dOnVi1ymTPqqoq9tprrwY7xLVo0YIBAwZQVVUVM94NGzawfPlyjjrqqAb+gwcP5s4772TTpk11W6uGxw/QpUsXli5dCkC7du0YM2YMxx57LEOHDmXw4MGcdNJJ9OrVq1HpICL06dOnQRwAq1atinmvK6+8kjFjxjB58mSGDBnCiSeemHCHvtraWm699VYef/xxli9fzpYtW9i+fTtlZWVxr5s7dy7btm2ja9euDfy3bdtGz549Y17XunXruuNBgwbxwgsv8Ic//IGtW7dy/fXXx43ziiuuqDs+6KCD6NevH/vuuy//+te/OOOMM+Je60mNmH0WqjpTVWcCR6jqSFWdZt0ZwKCmRCoivwR2AH9ryn3CUdUHVLW/qvYP7fzlSQOlpTBxIhQVQUmJ+Z04MSud3AMGDGDBggUsXLiQLVu2MH36dCorKxuEiezUFBFqa2szLluy8T/44IPMnz+fY445hpkzZ9K7d+9dFF4imjVrRl5eXoM4gLjPed1117Fo0SJOO+00PvroIw477DB+9atfxY1nwoQJ3HTTTVxyySVMnz6dBQsWMGbMGLZt2xb3utraWkpKSliwYEED9/HHH/PCCy/EvC487EMPPQTAyy+/zJw5c2jRogX5+fl1/TfnnHMO+++/f8x7VVZW0rlzZ2dGte2OJDMaqpWIVKrqEgARqQBiD1FIgIicCxwPDLPNHoAVQPewYN2sH3H8Pdni7LPhmGOyPhqqqKiIfffdN+Xrv/e977F27Vo+/vjjutbF1q1befvttxk3bhxQX9jv3Lmz7ro2bdrQrVs3Zs2axfHHH1/nP3PmTCoqKupaFcnSu3dvevfuzeWXX85FF13EAw88wNixY1N+rmSprKxk3LhxjBs3jptvvpnbbruN3/3ud4B57vBnBpg1axbDhw/n/PPPr/P77LPPGoSJdl3//v1Zv349W7Zs2WWf7HhEe7eTJk3iu+++a+B30EEH8fvf/55TTjkl5r1WrFjB6tWr6d69e8wwnqaRjLK4DHhNRJYAApQBKeV0ERkO/AIYrKqbwk49A/xdRP4IdAF6Au/Y+HpaBbUCGAX4NmYQlJZmTUmki6FDh3LooYdyxhlncO+991JSUsKNN97Ili1buPjiiwEoKyujWbNmPP/884wcOZIWLVpQUlLCtddeyxVXXEHPnj0ZMmQIr776Kvfddx/33ntv0vEvXryYBx98kBNOOIHu3buzcuVKZs+endAc1FQ2btzI1VdfzSmnnEJFRQXr16/nxRdfbGCOq6ioYMaMGRx33HEUFBTQoUMHevXqxZQpU5gxYwZdu3bl0Ucf5e2336Zdu3YNrnvyySepqqqic+fOFBcXM3ToUI4++mhOPvlkbr31Vr7//e+zbt065syZQ2FhIRdccEHSslfEMG9269atzqT15ptv8sYbbzB06FD22msvPvnkE6655hp69OjBSSedlGKqeRISqzMj3GFGKfWxrkWS1zwGVAPbMX0No4HFmD6IBdbdHxb+l5iRT58Cx4X5/whYZM/9Mpm4fQe3IZc77KKNhgon1ME9e/bsBv777LOPXn/99XX/V65cqSNHjtSSkhItLCzUo446SufOndvgmltuuUW7dOmizZo1q+u8ra2t1VtvvVXLy8s1Pz9fKyoq9E9/+lOD68rKyvTGG29s4BfeAbxy5Uo96aSTtGvXrlpQUKB77723jhkzZpeO4HCidXCH/1dV/eqrr+J2nG/evFlPP/10LS8v1xYtWmjHjh31tNNO0y+//LIuzAsvvKD777+/Nm/eXE0xoLp+/Xr96U9/qsXFxdq+fXsdN26c/upXv9KysrK669auXavHHXectmnTRoG6TvJNmzbp1VdfreXl5dq8eXPt3LmzHnvssfrKK6/EfNZkIaKDe/78+Tpw4EBt166dFhQUaGVlpV500UV1o7pSJZe/l3RBnA5u0TpLUHREpCVwOVCmqheISE+gl6o+23RVlRn69++vobH0ezJ+A3qPJ3n89wIiMl9V+0c7l8ykvEnANswwVjDmoN+lSTaPx+Px5ADJKIt9VPVWjDkJNX0NklGpPB6Px+MUySiLbSJSBCiAiOwDbM2oVB6Px+NximRGQ90AvAh0F5G/AUcA52VSKI/H4/G4RUJloar/EZH5wGEY89Olqvp1xiXzpAVVrZvE5fF4opNooI8nCTOUiLyiqmtV9TlVfVZVvxaRV7IhnKdpNG/evMEKrR6PJzqbN2+mefPmQYvhNDGVhYgUikh7oIOItBOR9taVU79uk8dhOnXqxIoVK9i0aZOvOXk8UVBVNm3axIoVK+jUqVPQ4jhNPDPUWGA8Zkb1fOpHQG0A7smsWJ500KZNG8Ds+bA9cjFAj8cDmBZ4586d674XT3RiKgtVvRO4U0R+rqp3Z1EmTxpp06aN/wg8Hk+TSWbobG3EjnbtRGRc5kTyeDwej2skoywuUNX1oT+qug5IfmUwj8fj8eQ8ySiLPAkbeykieZgtTj0ej8ezh5DMpLwXgakiEtqxZaz183g8Hs8eQjLK4mqMgrjY/p8OPJQxiTwej8fjHMnM4K4VkcnAq6r6aeZF8ng8Ho9rJDOD+0TMRkUv2v8Hi8gzGZbL4/F4PA6RTAf39cChwHoAVV0ARN/70OPxeDy7Jckoi+2q+t8IP792hMfj8exBJKMsqkTkDMwQ2p4icjcwJ9FFIvKIiKwWkY/C/NqLyHQR+cz+trP+IiJ3ichiEflARPqGXXOODf+ZiJyTwjM2ipoaePNN8+scDgoXLpJz4jktnOPiRQrkmIBOpx3kgIApEGtz7pADWgK/B+Za9zugRRLXHQX0BT4K87sVuMYeXwPcYo9/BLyAWX/qMOBt698eWGJ/29njdoni7tevX0qblT/6qGpRkWqbNub30UdTuk1mcFC4cJGaNzfOGfGcFs5x8SLz2oUXOpX3nE471RwQMDbAPI1Vpsc6ofUF/E+T8YtxbXmEsvgU2Nse7w18ao8nAqdHhgNOByaG+TcIF8uloiyqq827hHpXWKg6bZo5FxjV1UaIwsKGwhUVBSpYtPRyRjynhXNcvETCBSyg02mXEwLGJ56ySMYMdW2SfsnQWVWr7XEN0NkedwW+Cgu33PrF8t8FEblQROaJyLw1a9Y0WrClSyFyOfstW2DUKKishClTGn3LpjNliol85EgjTDj5+fDMM4E1baOlVziBiue0cI6Ll0g4gIICEy7L1NSYdMmPM+A/INEMzgvYNOLtZ3Gc7Z/oavsTQm4ysKOpEVstlraOclV9QFX7q2r/jh07Nvr6igqItor3d9/B5s0wdmyWP96aGhPp5s2wadOu57/9Fq68MjBNFiu9QgQmXk0NrF0L27bFDuPTLjaJhAPYutWkcRY/iFC96YorTPrEYts28whZx3kBm068lsVKYB6wBbOfRcg9AxybYnyrRGRvAPu72vqvALqHhetm/WL5p53SUpg4EYqKoFWrXc9nvUIQq4ZXVFR//O23gWiymhoj3i23GHFKSoyoBQVQXBygeKEP9swzYedOI5AzwhmcTbtYwhUVwYUXNhR2506TxlnSZuH1po0b6/2Li+vTLiTqLbeYRwisYpdIwIkTTWGTi8SyT4Uc0BwoAnolChvl2nIa9lncRsMO7lvt8Y9p2MH9jtZ3cC/FdG63s8ftE8Wbage3qkNdBLE6UW67TbW4uKF/SYnqnDlZESuy7/Ouu0zU1dXGTZwYkHiJOp0CFc7gbNolEk410A9jzhwjVni0xcUmvUJpN2eOETmQfvhkBYw8dhCa2MF9AqbDean9fzDwTBLXPQZUA9sxfQ2jgb2AV4DPgJdDBb9VEvcCnwMfAv3D7nM+sNi68xLFq01UFiFC305JSYCDGKIJEa1QzJImSybqwMSL9sFGlrQ+7VIXTjW5NA5IvABfbfKROziiMZKmKov5QAnwXpjfh4muC9KlQ1moBlgJSFQTCUiTJVtWRIoXWUnNCI39YHMk7bIiXrLCBVgiJ0qXgPRY8gIGqs2Sp6nK4i37+16Y3weJrgvSpUtZhJM1xZFs7SOApm1j8nsgpoFkS9ocSbuslSONES7AJne8dAm0VZZMXgpcmyVHU5XFw8AZwAdAT+Bu4P5E1wXp0q0sstZ6TCXHZ7lp25iyIpAPuDElrcNpl3Ua+2IdtLtnPX0bk3/2kJZF+Azuefa4MNF1Qbp0KousvuPG1j4CyoDJlhVOV6YcT7tAcFq45MjaIzSlYudkbcEQT1kks5/FJuCX1u1xhEawbt5c7xcaRpv2EXDRxrjHG5edVeHqKS1N7vaNfZys4njahQiNZq2oyMKIy8YKBxkXsLG3T+URUiKV/HP22XDMMVl8oekl3qS8aSLyTCyXTSGDJKsFXvhkj2TGZTtdGjf+cbKK42kH9dNGhg8PcBWBeGRYQKefP9X8U1oKAwea31xbYDBWkwMYbN2dwFTMENoTgL8Df4p1nQsuU30WWWs9pmJ3z7BwTWneO2vdcDjtsmYlS/XlZFjAdNw+4/muKfnH0WG0NLHPYpeL493QBZeTo6EcLo0dzdf17IZpl5X+nqa82AwL2NTbZ3VQirM1gcbTVGWxEKgM+18BLEx0XZAuE8oiozhcGjucrw27adplPN2bGoHDLQvn86zDIz/iKYtkVp29DHhNRF4TkZnADODSdJnB9njC15XZsCGgVQtjE22JKmcWztyN0y7j/T1NfbEZFrApt3c6z0JO9JdFI5nRUC+KSE9gf+v1iapuzaxYbpPWASDpHpWT5tEpTufrgEY0JUtT0y6jg2fS8WIzPLon1ds7nWehXhOOHWvy67ZtDo38iEOsJkcuu0yaodJu9UhnmzlDJhlnh4dnwt6Q5j4MZ9NO1XHhmkZOPFoG+suaekua0meRiy5TyiJjttB05Ows2JDTka/T/n2ks1TIkLJ1djSYquPCNY2MPJrD6ZWO7OuVRZrIaL9UUzOhw51mITLWF52OD9j5XlFP4OymgynCiacsEnZwi8gRItLKHp8lIn8UkbIMW8ecJKO20PDJOqnguKE2o33RTU07yIFe0RyYw5UmAZ18zt14MEWyJDMa6j5gk4j0Aa7A7DnxaPpEyB2cnpGcAeHS+dE6XxY7rmzTOps5E6VxmgR0dta24xk4K9k3VpMj5IB37e+vgdHhfq66TM+zcNhsmTbh0t3izgkrj6O9omlNu0yYUtIkoNN5xGnhDOnIvjRxUt5M4FpgEVCKaY3sEZsf7alk6rtwtCxuiIM1gbR1R2XqxaZJwEx2u6XltWYyA6cp32VyNFQyZqiRwFbbqqgBumH20vbspmSqxX322bBkCbzwgvk9++ym3S8jpKP/g/RaetJmYsjUi02TgJkypaTNtJWpDJxG21uasm90YmmRTDrMrPAq4CPMXt2FmGVE3sbstT0VKLBhW9j/i+358kT3z5mWhYO1WNWcaHE7m3aqmbH0pKVSm8kXm6Zad7or787nZccEpIlmqJOBz4D/AhuAb4ENia6Lc7+uwFKgyP5/AjjX/o6yfvcDF9vjcdid+YBRwNREceSEsnB4GJ6q4yYjh9Muk9++N6U0HudHlDsmYFOVxWLggEThknVWWXwFtMcsN/IscCzwNZBvwwwEXrLHLwED7XG+DSfx4nBeWThWm4iFk5V3x9POsW8/Ok6+2MzgeHZxTsB4yiKZPotVqrowiXBJoaorgNuBL4FqTItlPrBeVXfYYMsxSgXqlQv2/H+BvSLvKyIXisg8EZm3Zs2adImbkJRs044PwwuRUfsnu2faOT4C15DpF+sQTg93hxwQsJ5klMU8EZkqIqeLyMkhl2qEItIOGIHpo+gCtAKGp3q/EKr6gKr2V9X+HTt2bOrtkiLlfqmcKFEyy+6adjn07e8xOD+wwnkBDckoizbAJuCH1O+Wd3wT4jwaWKqqa1R1O/AUcATQVkRCq+B2A1bY4xVAdwB7vgRY24T400KTJnRmq0RxcipsjqRdE8iRb3+PokmNqWx8RznQ2ktmifLz0hznl8BhItIS2AwMA+Zh9sk4FXgcOAd42oZ/xv5/055/1drWAqXJq2NnevP2KVNMCdy8uamJT5zoTKnlfNqlgdJSJ8XyNBaHv6NsI4nKXRHpBtyNqf0DzAYuVdXlKUcq8hvM/I0dwHvAGEzfxOOYju/3gLNUdauIFAJTgEOAbzAjppbEu3///v113rx5qYqXFDU1xnwSXuAVFZmaZOCFhNPCOS+eJx2keV+VQNgDM6qIzFfV/tHOJWOGmoSp3Xexbpr1SxlVvV5V91fV3qp6tqpuVdUlqnqoqu6rqj9Vu8GSqm6x//e15+MqimzhtDXE8U5gp9NudyWbJskUOqSctJg6/h2FyFraxRomFXLAgmT8XHLZHDrr5ChEx4bjxcLJtAvHeQGTJJvzUlLIe85Om8mB7yjdaUcT51m8ApwF5Fl3FvBKouuCdM7Ps8gGTs+qywGcLcEMSeuxbBd4jZxo4nx57PB3lIm0i6cskjFDnQ+cBtRg5kWciplx7XGZ3XFITrba247vXdAoK0+2TSmNHNochKWnUdnI4e8o22mXjLLopqonqmpHVe2kqj8BemRGHE9ayYHheEmTzY0OHLZVN1qPZXteSiM7pLItXkrZyNHvKOtTjmI1OUKOKHtXRPNzyTlthnLcDu6keNm2VThsG0lpOZEgTCmNyEjZEs/h19qQANOOVPosMOszXYFZauPyMHcD8H6s61xwzioLx+3gzooXxIJLKXyF2VC0KRd4TtYC6smGeDmxblcKH2E60y6esog5z0JEBgNDgIswq8CG+BaYpqqfpb2ZkyayMc+i0Tg+Zttp8YISrhFzBbI5dysUV0GBMTvswfPEGoXTeRycEDDePIuYM7hVdSYwU0Qmq+oXItLa+m/MkJy7N02etpxZnBYvZAePLCEzLViS07DD+xFC6Td2rJlkngkRc2ACu5MElY2SxumPMInlPoBiEXkPM7MaEfkaOEdVP8qoZDlI3Iqo4wvgOS6e0yVkEN+4X04kNRzORs5/hMmMhnoAuFxVy1S1DNOP8UBmxco9Eo6yCHracoLxgkGLlxR+VIonDTiajZz/CJNZG+p9Ve2TyM8lst1n0ShTYxBr5jTCoB7Ukj65vpSQ70fwpI0AP4Z4fRbJKIt/Ae9iFvMDM4O7n6qelFYp00i2lcWbb5oWxYYN9X4lJWYez8CBWRMjOg50miVid1nYM9cVnsfT1IUEzwc6YvadeMoen58+8XIfp80QDk8wA+cnSzcKZ80bHk8aSKgsVHWdql6iqn2tu1RV12VDuFzBaVOj05rMeV3m2VNxchncYEmoLERkPxF5QET+IyKvhlw2hMslnF1CxmlN5rwuq8cXHnsO2VxaJodIqoMbMylvPrAz5K+q8zMrWuo4OSkvaBw2qDvfOby7dKp4EpMDfXyZpKkd3PNVtV9GJMsQXlnkHs7qsgSFhwtyuyDDboPTo1XqydQ7T6mDW0Tai0h7YJqIjBORvUN+1t/jSRvOdg7H6VRxwVoRVwZXTGeuyJEMOWAXDSzfxVo0ClgKLLG/kW5JrOtccM4uJOjJPWKs3Lfq/erAVzCNu6igK6tCxpHD2bUN97ANj8Ihlc2PVLVCVSvtb6SrbIqCEpG2IvIPEflERBaKyEDbYpkuIp/Z33Y2rIjIXSKyWEQ+EJG+TYnb42kUMQYIfP5daeCjuGI1epbPc2Q8cpxx0S60ykIi7tLocWm0SoSAQY4eTGY0VKGIXC4iT4nIP0VkvIgUNjHeO4EXVXV/oA+wELgGs11rT8xWrtfYsMcBPa27ELiviXF7PI0jSuHhgrUipgw4Mh45Rsm2dt5SJ3RZXIXlgl00ioCB5rtYTY6QA54AHgb+x7oHgScTXRfnfiUYU5ZE+H8K7G2P9wY+tccTgdOjhYvlvBnKkw1csFZElcGVXX5iyDF3WnXg+0q4kkQxiSNgJvMdccxQyaw621tVDwz7P0NEPm6CfqoA1gCTRKQPZkjupUBnVa22YWqAzva4K2YDphDLrV91mB8iciGm5UGPHn7X190Gh4f6uLCCaXQZHFmLO8aa4N36lwbeKnN8NfC4Ap59dmkw+S6WFgk54K/AYWH/BwCPJrouzv36AzuAAfb/ncCNwPqIcOvs77PAkWH+rwD948XhWxa7Ca500uYqrvQgR5Ej6FZZLrcsMgmpdHCH0Q+YIyLLRGQZ8CbwAxH5UEQ+SEE/LQeWq+rb9v8/gL7AKhHZG8D+rrbnVwDdw67vZv1yh1waOugKu9OiUUHhgt09hhxB9yE7vrCBkwImY4Yans4IVbVGRL4SkV6q+ikwDPjYunOAm+3v0/aSZ4D/E5HHMa2a/2q9ucpZQtaTA+ZNoe3VfvZvo3HeTuBpKkFv4OSCGTEujgmYcAZ3RiIVORh4CCjAzOU4DzMy6wmgB/AFcJqqfiMiAtyDUVqbgPNUNe707KBncIdWh+iaV8P7GytpiWNLB0T0AzjZLbCHL7vg8QRBSntwZxJVXYDpu4hkWJSwCvxvpmVKF+HWkw4sZRvNGyqLoGvHEescvX72RH445Wz3Gj7Ob5js8exZxGxZiEgLVd2aZXnSQpAti/ClZTpTwxIcallEqa1voohKlrCK0sDFi4qLzR4XZfJ40kCqmx+9aS/26/M2gvBJM6soZSwT2UQRtcUOdFJFmSS1nQIzicviwl4SDcYDuNJJG8KVqcceT5aJpywKROQM4HAROTnSZUvAXCNyEMM/i87mubuW0OwlB5YOiDL9sznbWEr9APeg10xzuiz2I7Q8ezDx+iwuAs4E2gInRJxTzBarnijsOoihFHCgZhylH+DdsyeyYUopJQ50C4SXxSFL2dixJi2daFjkwggt101krsvnMIEnXawJGFo/CW50ojCuOT8pLwERk6Rcmbs1Z44GvgxEXFyfyeX6JEbH5XPlO4hGtpKOOJPyklEWBcAlmMlz/wB+DjRPdF2QziuL3MT1slhVg596HAvXE89x+R59VLW8sFqHtZyj5YXVzrxW1fqk60y1HsYc7Ux1xpIunrJIZgb3nzGzuP9sXV/8yq+eDODgpNVdCXrqcSyWmmHa4WzDgdEKIRyWr6YGZoyeQtWWSv6xaThVWyqZMXqKM11RS5fCWTqFJVTyAsNZQiVnMSXrSZfMPIsfqGqfsP+v2n25PZ60c/bZcGyfGla9tZTOh1XQ6fsuaQpL0FOPo7C6VQWtN2+nIMxvx+ZtrG9VQafApKrHZfmWz6vhnu1jacnmumHu92wfy8fzjqH0+ODf8z6tarhjS0P57tg8lo2tjiGbfaHJtCx2isg+oT8iUgnszJxInj2aKVPodFglB101nE6HuTYcyuDiUl+ff1fK+EIzTHs9JWyiiPFFZpMmF3BZvgqWsp34Q8qDpNN3S8kvaihfflEBnb7LrnwJl/sQkWHAJMyyHAKUYZbcmJF58VIj6OU+PCmSA0t8REyAd2bGeyjp2myuoYKlLKWCDUWlziSd0/LV1LCjRyX52+vz3Y6CIvK/cEE4svpdpDopDwBVfQWzS90lmM7tXi4rCk9iXKwZA8HuGZkELk+zCPX3bCgqZWHJQDYUlTrV3+O0fKWl5D88ES0sYkerErSwiPyHXBEOZzrzAllIMNM40bIIfFB0dKLWjI9xRFbHWxbhS7mEKCkxfd0DBwYnVziOZrs6nJbPaeHIinzxWhaBD3PNhAt86Kyj48mjjV48r/mjWlvokKyuDk1V50d/ejxNhjhDZ33LIt1EqR1rYRHzn1xCt/6lgVZYImvGnalhKZUUubLQYQiHa3ihlln4Qrgu9Fl4POmgSX0WIvJKMn4eSxS7+4YtBVwzcmngax1FLg1Vwa5j353oI3Bt8cAwXJ1m4fFkmpjKQkQKRaQ90EFE2olIe+vKga5ZkzDXiLFY30ebKgLvEI3sJ6sprKBV84ayBr6SYAQudsY7rMs8nowRr2UxFpgP7G9/Q+5pzM51nmiElcg7Wpnx5GOZWLdfRNAV9/Ca8ZtLzSiQoEdZxMLpFWg9nj2MZOZZ/FxV786SPGnBldFQa+ctpf9PK1i2pb7wdaFLYBcc7CNwfGCUx5MVsv1pNnWexd0icriInCEiPwu5NAiVJyLviciz9n+FiLwtIotFZKqIFFj/Fvb/Ynu+vKlxZ4XSUvY6fiC/faDU1Yp7PQ7aVRyfcuHZXXDRzmlxrWWdTAf3FOB24EjgB9ZFH4fbOC4FFob9vwX4k6ruC6wDRlv/0cA66/8nGy5n8B2iqRGl68e17hTncbgcNAQtoGulcRhOTgCNNaY25DAFuiQK1xgHdANeAYYCz2KWEfkayLfnBwIv2eOXgIH2ON+GiytP4PMsPGnB4SkX9biwCUIUGVya6hM1iYIW0PFJMwn3dslQvqOJ+1k8CeydKFxjHGZfjH7AEKssOgCLw853Bz6yxx8B3cLOfQ50iHLPC4F5wLwePXqkNQE9weFCWRyToAu8GDK4VA5GTSIXBHR8p624SZTBfNdUZTEDYxZ6CXgm5BJdF+d+xwN/tsdpUxbhzrcsPBnHhQIvhgxzp1U7UQ7GSqKvpzlQULvw/hIQtWWdYbnjKYtk9rO4IYkwjeEI4EQR+RFQCLQB7gTaiki+qu7AmKlW2PArMMpjuYjkAyXA2jTL5PHEZZdRKS7sxx1DhgqWsn17QxmC6O+JmURUsFfQHVJR9qMPH33iwgDBs882+883kOPNAPNdLC2SDYdtWdjjJ4FR9vh+YJw9/l/gfns8Cngi0X19y8KTTpw1pcSRwYX+nqRMKUF3SDne37MLAbYskinQvwU2WLcFs/HRhkTXJeMilEUl8A6w2CqOFta/0P5fbM9XJrqvVxaedOF8gRdHBhf6e+ImkQsCRuBCHSAhGcx38ZRFoxYSFBEBRgCHqeo1jW3FZAsnJuV5dgsSLkvugr3CBRni4Lh4DciFZeiBjCVqvEl5yfRZ1GE1z79F5HrAWWXhyUEcLVESzvdwYT9uF2SIg+PiNSBn5vcEkKjJTMo7OcydKiI3Y8xRnhwi6PlPcXF4cpQjm5R5soR/37FJZm2oSWF/dwDLgAdVdXUG5WoSgZmhHK0du7pvNJAzi0A5+mo9GWJPfd9NMkOp6nnpF2k3xNESOXzZgFB5PHasGZLnxEfgwhDUJMglU4qn6fj3vSvJmKG6ici/RGS1df8UkW7ZEC5ncHIhF4PzC/LljJHY49mzSagsgEmYWdtdrJtm/TwhHC6RnS+LdycjsdMdQx5P00hGWXRU1UmqusO6yUDHDMuVWzSiRM52edLosjiIAm93WJrX4U56T27hbJ0j1gSMkMOsDnsWkGfdWcAria4L0gUyKS+JiTJBzgxNav6T01NXHSYnZnJ56nBwMmCIoD9BmjiDuwxjhloDrAb+DfRIdF2QLrAZ3HEyofPlifMCOvyNO76CqSeMoEvjOLjwCcZTFsnslPeFqp6oqh1VtZOq/kRVv8xAIyf3ibPjnMPdGgbHBXTayhNEx5CztgqHcXggCjj/CSY1GqpCRP4oIk+JyDMhlw3hdiec72h2WMCgvvGky+Nsd9KnoDmd1y3ZENDx0tjhT9AQq8kRcsD7wCXA/wCDQy7RdUE6VxcSdGHdubg4KmAQVp6UrBXZsJOlYKsIwvLSqKTIloAu2HkSkNInmMZ8RxP7LN5OFMY156qyUHXY7h7CQQGz/Y07XaY0UnMG8SyNKvuzLaCjFaJwglS0TVUWZwDXY/bF7htyia4L0rmsLLKJg+V+ymTzG3e6v7qRhWu2n6XRZX8Qib27fBgZULTxlEUyq84eBJwNDAVqQ9Yr+9/jKI6uPpIyUXcNyxBO244T7PAWSbafpdGrtwSR2I1cy8PZdaKyvVROLC0ScphNhwoShXPJ7ektC6fNKDmC89aKRtSOs/ksKeU9hxPb4ZG2WW9ZJKMs/g10ShTOJZdVZeFgk9ZpM0o4DqZdOI6L1yiy+SxBd9Kmi5yodKVZ0cZTFsmYodoCn4jIXGBrWIvkxDQ3cnIPR209TptRQjiaduHsTiuPZvNZUjIZOpjYObEgchbts8nsZzE4mr+qzkwpQpHuwKNAZ0zfxwOqeqeItAemAuWYPTNOU9V1divXO4EfAZuAc1X13XhxZGU/C8f3YQiVxeFmbWfKYsfTzuOBPTObxtvPIpkZ3DPDHbATOK0J8uwArlDVA4HDgP8VkQMx27S+oqo9MetRhbZtPQ7oad2FwH1NiDt9pGGCTybnITm9Np/jaZcWnBfQk4jdaUHkdJDMqrOIyCEicpuILANuBBamGqGqVodaBqr6rb1XV2AE8Bcb7C/AT+zxCCBkiHsLaCsie6caf9pooq0nG8tXxFl9JFhyIO2ahPMCepLF6UpXtonVmQHsh5lf8QnwOvBz4ItY4VNxGJPTl0AbYH2Yv4T+A88CR4adewXoH+VeFwLzgHk9evRoUidP0qTYuZQTHWeZZndNO+cF3ENxsAM9HFfEI8WFBD/BzKU4XlWPVNW7MSaotCAirYF/AuNVdUP4OSt0/M6UCFT1AVXtr6r9O3bM0nYbKVY7HF+iph4H7WTOp53zAu6BON7Sc1y8OuIpi5OBamCGiDwoIsMwNf4mIyLNMYrib6r6lPVeFTIv2d/V1n8F0D3s8m7Wzw1SsPXkzGglB+1kzqed8wJmrg7gZDeN4yvNOi5eA2IqC1X9t6qOAvYHZgDjgU4icp+I/DDVCO3opoeBhar6x7BTzwDn2ONzgKfD/H8mhsOA/6pqdarxu4DzHWcO52Dn0y6TAqahNM5UHSCt902n1nG8pee4eA2JZZ+K5oB2mL6BlHfKA47EmJg+ABZY9yNgL0x/xGfAy0B7re+/uBf4HPiQKP0VkS5XZnC7YqfchRyY1eds2oVIt4BpmEqcqe6UtN433VOm0/zQ6X6trnVx0ZQZ3LnockVZOEsmc7DjpbyT4qXpfWSqDpC2+2Yq36VplnOmlv5wabUTryw8jScTOdjphXYcFi9NpbHzLYtMtmibWAvIdAsgbZWUJt7IK4t042T105BW0dJ5M9fa2xE4LV4ahctULTYt93X4JeSAZTYttR2vLNKJs9VPp0Vz/mtzXLy0lvKZquuk5b4u2WTCcFiPGdIkoFcW6cLhHOOwaIYMCpiOQsr59FN1ukWbVhx9Tkf1mCFNtZ14yiKp5T48lgyOc2vqaEHnh+BlaEhpuoZsOj8kFzKyfouTcyPS9Zxpfjinl/7IxvyeWFokl12utSzSYT7KeM04XbW9NNYaM/HMGanUOl5TdtJs2VR264eLQRqaPngzVBpJc1s0nQVexprJjn54zvczqDY97TKkaNKV7zKuB1OJICdsim5WTLyySDdpfMvpLvDSngEd/vAcFs3QVAEzqKTTke8yXodINYIcqEU4Wv/yysJlnC/wHP/wdttOxwxnjKbePitmz1QjcPyjclm8eMrCd3AnS4Z6AtPRsZrRTkrHF8Zraqejs2mX4RELTc13GR9Q0ZQIsjhaIZX84/xglFjE0iK57NLesshCmzFV81FWmrNNrb7vyR28qaZdlqqfqb4ap1sW4ffIYL5LNf/kassi8II9Ey6tyiKAN5tsHs+qaE5rs8aTU2nnpI0tC+I5/Pzp6o5y7dG8smgKWbbZN6Zsdbw7wekqlPNpF8LRVlkIJ0dDZYF05B8XHy2essgP2ArmPlm02YdvI7F5s/EbOxaOOSa6udXx7oR642zoYaDeOBvwbDfn0y5EaWngaRWPjIvn6POnmn9qakz2r6hw9tFi4ju4E5HFzrLGdnwFNus42V69AErkZEVzPu0CIFnRHH6ErAmXSv7Jle1TYxKryZHLLiNDZ7PQZkzVapPV5mxj+yCyaJxNVrTw9HI67bJIsqIF+giJXlYAwjnZv9gE8H0WuUOyZWsg9s7G5Pgsl8jJihZYYZdq2jkkWqAFXqIX50BpHO+15UofWTxl4c1Q8QigvR05b+CYY3YVIbDmbLJ2skgBp09P+wJ4qYgW6NbiqaZdFl5usqIFNj8gmRcX8OSFRK8tZ/rI4pAzykJEhovIpyKyWESuyVhEIQVx992BGRhDi25On95QhLvvhmefhQsvDKjAi5bjt26FtWuNADU1gQmYjGjPPAP5EUM6slae5HDagfldu9YUcOFkpcCLpgjy880LDaVdYMJF12UXXmheZ0i8pUvhllscX9U4EbGaHC45IA/4HKgECoD3gQNjhU/ZDBVq6rZu3bC9GFCTNrJVDaotW+7ql9XmbLidrHlz49q0qT8OUMBEogX+WnMw7YqKVC+8sN4C1Ly5akFBlucHxPoYiot3TcesCxfdxASqrVrtmpZ33eXecNlwyPU+C2Ag8FLY/2uBa2OFT0lZxMqQARkYY2XAaC7rHWXV1arTpqkWFjonYLKiFRcH1Me8G6RdYaEJl9U8F9JmxcXOCZeo6HC9UzuceMoiV8xQXYGvwv4vt351iMiFIjJPROatWbOm8TFEa+qGk2UDYzTTQCStWgXUnC0thb32MjaceAQgYDKiFRfD7bcHtIFNjqcdQIsWJlxW81yoM+/2280LjEUAwoUPo23VKn7YnFgDKga5oiwSoqoPqGp/Ve3fsWPHxt8gVulcXBxIiRyeAaN9G4WF8PjjAe7YlUibBShgItF27IATTwzQXpzDaQcBdsyWlpoXt2NH7DABCRfSZY8/bl5fLHKtU7sBsZocLjmyYYZS3XXcqgMGxtBwvLvucnAtmUhDdwD24hwUzeCwgJGfQajPwhHxnE47VefFiwtxzFBizruNiOQDi4BhwApgLnCGqlZFC9+/f3+dN29eapFFzsd3CCdFCxcKnBLQYdEMDgsYmdecy3sOpx04L15MRGS+qvaPei4XlAWAiPwIuAMzMuoRVf19rLBNUhYej8ezhxJPWeTMQoKq+jzwfNByeDwez57IbtPB7fF4PJ7M4ZWFx+PxeBLilYXH4/F4EuKVhcfj8XgS4pWFx+PxeBLilYXH4/F4EuKVhcfj8XgS4pWFx+PxeBKSMzO4G4OIrAG+aMItOgBfp0mcdONlSw2XZQO35fOypY7L8kWTrUxVo67Eulsqi6YiIvNiTXkPGi9bargsG7gtn5ctdVyWr7GyeTOUx+PxeBLilYXH4/F4EuKVRXQeCFqAOHjZUsNl2cBt+bxsqeOyfI2SzfdZeDwejychvmXh8Xg8noR4ZeHxeDyehHhlEYaIDBeRT0VksYhcE7Q84YjIIyKyWkQ+ClqWSESku4jMEJGPRaRKRC4NWqYQIlIoIu+IyPtWtt8ELVMkIpInIu+JyLNByxKJiCwTkQ9FZIGIOLX9pIi0FZF/iMgnIrJQRAYGLROAiPSy6RVyG0RkfNByhRCRy+y38JGIPCYihUld5/ssDCKSh9nn+xhgOWaf79NV9eNABbOIyFHARuBRVe0dtDzhiMjewN6q+q6IFAPzgZ+4kHYiIkArVd0oIs2B14FLVfWtgEWrQ0QuB/oDbVT1+KDlCUdElgH9VdW5iWUi8hdgtqo+JCIFQEtVXR+wWA2w5coKYICqNmWicLrk6Yr5Bg5U1c0i8gTwvKpOTnStb1nUcyiwWFWXqOo24HFgRMAy1aGqs4BvgpYjGqpararv2uNvgYVA12ClMqhho/3b3Dpnakgi0g34MfBQ0LLkEiJSAhwFPAygqttcUxSWYcDnLiiKMPKBIhHJB1oCK5O5yCuLeroCX4X9X44jBV4uISLlwCHA2wGLUoc18ywAVgPTVdUZ2YA7gF8AtQHLEQsF/iMi80XkwqCFCaMCWANMsia8h0SkVdBCRWEU8FjQQoRQ1RXA7cCXQDXwX1X9TzLXemXhSRsi0hr4JzBeVTcELU8IVd2pqgcD3YBDRcQJM56IHA+sVtX5QcsShyNVtS9wHPC/1hzqAvlAX+A+VT0E+A5wrZ+xADgReDJoWUKISDuMxaQC6AK0EpGzkrnWK4t6VgDdw/53s36eJLD9Af8E/qaqTwUtTzSsmWIGMDxgUUIcAZxo+wUeB4aKyF+DFakhtiaKqq4G/oUx17rAcmB5WCvxHxjl4RLHAe+q6qqgBQnjaGCpqq5R1e3AU8DhyVzolUU9c4GeIlJhawSjgGcCliknsJ3IDwMLVfWPQcsTjoh0FJG29rgIM4Dhk0CFsqjqtaraTVXLMfntVVVNqpaXDUSklR2wgDXx/BBwYjSeqtYAX4lIL+s1DAh8QEUEp+OQCcryJXCYiLS03+0wTB9jQvIzKlYOoao7ROT/gJeAPOARVa0KWKw6ROQxYAjQQUSWA9er6sPBSlXHEcDZwIe2bwDg/6nq88GJVMfewF/sqJRmwBOq6twQVUfpDPzLlCnkA39X1ReDFakBPwf+Zit3S4DzApanDqtcjwHGBi1LOKr6toj8A3gX2AG8R5LLfvihsx6Px+NJiDdDeTwejychXll4PB6PJyFeWXg8Ho8nIV5ZeDwejychXll4PB6PJyFeWXh2QUR22tUyPxKRaaF5CnHC3yAiVyYI8xMROTDs/29F5Og0iZySTJlERIbEWkXWLk1xoD3+qV0xdYa9JuYEKZuGv06TfHXpIyLnikiXaPLFuX6yiJyaDlli3P8xEflARC5LMvz/icj5mZLH45WFJzqbVfVgu7rtN8D/puGePwHqCiBV/bWqvpyG++YcqjombEXe0cAFqvo/mHk08WbT/gL4cwZEOhez9EM0+bKOiJQCP1DV76vqn5K87BHMvAtPhvDKwpOIN7ELKorIPiLyol1UbraI7B8ZWEQuEJG5YvaP+KedKXo4Zo2c22yLZZ9QzVTMHiJPhl1fVyMXkR+KyJsi8q6IPGnXnoqM7zURuTOsJRS+HMWB9vwSEbkk7Jp/22eoCi2OZxcbnGzv8WGoRpvkMw+W+r0L3gvNegZaS/1+C3+zM2ZDMve3rYQjgYdtGlwEXGbvMygijv2AraGlwq2s94nIW/b5hojZ82ShiEwOu25j2PGp4edCfpjl0f9m4y0KyRe6XkT+ZNPqFRHpGOX5+4nITJtGL4lZsh4RuUTMHicfiMjjUa4rFJFJNr3fE5H/saf+A3SNTAcROVBEakRkhT23RkS+FJFiVd0ELIt4/550oqreedfAARvtbx5mEbTh9v8rQE97PACzPAXADcCV9nivsPv8Dvi5PZ4MnBp2bjJwKmZm8JeYPScA7gPOAjoAs8L8rwZ+HUXW14AH7fFRwEdhMs0BWth7rQWa23Pt7W8RZvmKvYB+mBVpQ/dtG++ZI2SYBhxhj1vbZxoC/BezxlgzjNI9Mkzm/lGO69IxShznARMi0u9xQDALw20ADrJxzQcODn+X9vhUYHKUd1YnQxSZFDjTHv8auCfi/TW36dzR+o/ErH4AZunrFuHpGfFMV4SF3d/mg0KgPPQeo1wTLvc9wLlh534JXBH097O7Or/chycaRWKW7eiKWTdmuq3VHw48aSvIYAriSHqLyO+AtpiC86V4EalZZuVF4AQxyxD8GGNuGYwxW71h4yvAFLjReMzea5aItJH6PpbnVHUrsFVEVmOWr1gOXCIiJ9kw3YGewKdApYjcDTyHWZY72Wd+A/ijiPwNeEpVl9vw76jqcgCbnuWYjWdSYW/MktzhTFNVFZEPgVWq+qGNq8rGtSDFuMKpBaba479iFp4LpxfQG5NHwFQwqu25DzAtln8D/45y7yOBuwFU9RMR+QLYD6P4UmE1Rul4MoBXFp5obFbVg0WkJaaw/19MTXK9mqW+4zEZs0ve+yJyLqaGnYjHgf/D9I/MU9VvrclmuqqensT1kWvWhP5vDfPbCeSLyBDMypsDVXWTiLwGFKrqOhHpAxyLMQedBowniWdW1ZtF5DngRxjldmys+JN4llhsBkoi/EL3r42IqzYsrvC0SWr7zAREprUAVaoabUvTH2NaeycAvxSRg1R1RxpkiEUhJp08GcD3WXhiosYOfAnGXLAJWCoiPwWz0qwtXCMpBqrFLFl+Zpj/t/ZcNGZilpe+AKM4AN4CjhCRfW18razdPhojbZgjMZu5/DfOY5UA66yi2B84zF7bAWimqv8EfgX0VbMnR8JnFpF9VPVDVb0Fs3pxqrXbeGm0ENg3hXuuEpEDRKQZcFKMMPHibYYxNwGcwa4to0+BjmL3vxaR5iLyPRtfd1WdgTEhlmBamuHMxuYR+2572Pslyw5MizPEfjiyKu7uiFcWnrio6nsYc8LpmA97tIi8D1QRfdvZ6zC75L1Bw6XAHweush2Z+0TEsRN4FrP+/7PWbw1mlM5jIvIBxgQVqxDeIiLvAfdjRhfF40VMC2MhcDNGKYExub1mzUV/Ba61/sk883gxHeMfANuBFxLIEItpwEnROrgx/TeHSJg9LEmuwaTpHOrNQ5FMBu4PdXBHnPsOs2HUR8BQ4LfhJ9VsQXwqcItNowUY010e8FdrInsPuEt33fb0z0AzG2Yqpv9hK8nzFvA7Efme/X8EML0R13sagV911pPTWDPSlao6L2hZMo2I3Inpp8jakGMR2aiqu4xCcw0ROQS4XFXPDlqW3RXfsvB4coc/AC2DFsJROmBatZ4M4VsWHo/H40mIb1l4PB6PJyFeWXg8Ho8nIV5ZeDwejychXll4PB6PJyFeWXg8Ho8nIf8fMPNF84OqlcsAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(x_coords, y1_coords, 'b.', markersize=9, label='Photons in state +45')\n",
"plt.plot(x_coords, y2_coords, 'r.', markersize=9, label='Photons in state -45')\n",
"legend = plt.legend(loc='upper center', fontsize='x-large')\n",
"plt.ylim([-100, N_init + 500])\n",
"plt.xlabel('Relative phase shift (multiples of π)')\n",
"plt.ylabel('Amount of photons detected')\n",
"plt.title('Amount of photons existing in |+/-45> states');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The data shows a clear interference pattern. The measurements oscillate between only measuring photons in the state $|+45\\rangle$ or $|-45\\rangle$ as the relative phase shift increases."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4. Variations of Single-Photon Interference Experiment "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alternative cases can be investigated to get a better understanding of what is happening, why this interference pattern appears, and of the quantum nature of this experiment.\n",
"\n",
"In this case, it is shown that, when the superpositon of states collapses (here the superposition refers to the photon taking both paths \"at the same time\"), the interference pattern disappears. This is done by blocking either the vertical or the horizontal output port of the very first polarization analyzer, thus preventing the photon to be in the previously mentionned superposition of states as it can take only one path. This comes with some modifications, with respect to how the simulation was implemented previously:\n",
"* The main element that has to be modified is the first polarization analyzer. Depending on which port is blocked, the corresponding polarization component as to be put to zero for any further computation. \n",
"* The other big difference with the full interference simulation concerns the amount of photons that goes through the entire interferometer. When the initial $|+45\\rangle$ polarized photon goes through PA$_{HV1}$, it has $50\\%$ chance of coming out $|V\\rangle$ polarized and $50\\%$ chance of coming out $|H\\rangle$ polarized. Under the assumption that the experiment is carried out enough times to actually get half of the photons to come out from each of the ports, it implies that only half of the photons will get to the other end of the interferometer, as the other half will have been blocked. Where, for the original interference simulation, 1000 photons got out the interferometer, here, only 500 will do so. Because of how it was coded in this simulation, this requires that the number of iterations (for one phase shift value) is halved. This is directly linked to the final probabilities of the photon coming out in state $|+45\\rangle$ or $|-45\\rangle$. The sum of these probabilities is not 1 anymore (or $100\\%$), but 0.5 because half of the photons are not concerned by this. The missing 0.5 is thus the probability that the photon was blocked at the first polarization analyzer.\n",
"* The rest of the code is otherwise unchanged."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Output port 'V' of the first PA$_{HV}$ is blocked"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Polarization analyzer (HV) n°1, with V output port BLOCKED\n",
"\n",
"phaseshift1 = pi / 4 # CONSTANT\n",
" # should depend on real size of the setup (here: arbitrarily chosen)\n",
"PA_HV1vb = Qobj([[exp(1j * phaseshift1), 0], [0, 0]])"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEWCAYAAABWn/G6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABMYklEQVR4nO2deXhVxfn4Py9ZyEIIIAqyJQERpbiC4oZQwIordaniVkVwoy2Lu7/W1lorFOVbpVrFDYTWtaUqdUVFQHFFcQEUULCiAQVBBAQCeX9/zNybk5t7b05C7k2I7+d5znPvmTPLO3Nm5p39iKpiGIZhGOmiSX0LYBiGYfy4MMVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPPWAOCaLyDoReasW7l8RkeGpkK2hISLPisj5tXDXR0Q+aUgyebc3iMgNdSySkQQROVZEnqhvOcIiIheIyKtJnu90+ReRKSJy08744f2JyioiTUXkYxHZvTp3DVLx+IRdJyJN61uWeFSXMUJwFHAM0EFVD60jsUIhIitEZGA6w9wZVPU4VX2wOnsioiKyV8DdXFXtVp8y1TUi0k5EVoaw95JPj8yA2QoR+UFENvrrhRTK2S+MnDFuKr2/OubPwDgfzscicmGc8EeJyDspCj82rFTGtd5Q1a3AA8C11dltcIpHRIqBPoACJ9evNCmjCFihqpvqWxCj4RCiN3Q88Fw1fpwDZCV4fJKqNvPXzwJuWolIIje7NCJyCFCoqm94oweBX8axep5/lkpZMqu3tcvzEHB+dZ2GBqd4cJniDWAKUGk4w3cP/+6HOjaKyGsi0lZEbvM9pI9F5KCA/X1972m9iCwUkZMDzyp1V2N7Mb5VcqmILPXu7/RDZPsCdwOHexnWx4uEb50+JSLfisgyEbnImw8D7gu4/2Mctxf4uN0hIt/5eA2IsVbk7XwvIi+ISOuA+5N9fNf7eO7rzacBnYAZPuyrk9n3z1aIyJUi8oGX5VERyfHPWovIf727b0VkrojEzVMiso+IzPT2PhGRM7x5F292cCDdvhGRfrHvSUT2EpHZXo41IvKoN5/jg3nfx+vM2FZ3snj451eLSKmIfCUiw5O1SmNkukBEXhWRW30eXC4ix8VzVwccDzyT6KGIFAJ/AK6uob/HACtFZIKI9AjrSESOF5FFPg9+6dM3H3gWaCcVvat2InKoiLzu80qpz9vZ3p8q78+bnygiC7ybeSKyfyDsa3yY3/v8FFs+IhwHzA7cTwOOEpGigF/dgf2BhxPEc4WIXOfjuk7cMHkw71wkrox/K67Mtws8UxH5lYgsBZYmimv8YJOW/4ilJiLyOxH5XES+FpGpPh9Enh/l0269iHwhIhfE8aNARGaJyERxxC2r3u5uPo4bxE0TdAn6paorgXXAYQniFbXYoC5gGTAC6AmUAW0Cz6YAa/yzHOBlYDlOWWUANwGzvN0s79f/A7KB/sD3QDf//BVgeMDvC4BXA/cK/BdogausvwEGxbObIB5zgL97OQ/07vuHce+fbwfG+HicCXwHtArI/imwN5Dr78f5Z3sDm3CVSRauEloGZPvnK4CBgbDC2H8LaAe0AhYDl/pnY3FKOMtffQCJE5984AtgKJAJHOTfY3f//CJgEZAHPA/cGnAbfU+4iuG3uAZTDnBUzPvaK3DfD1gZuE8Wj0HAKuAnXoZ/xPoXE5+gTBfg8ulFuDx4GfBVvHTw9m8AbqjFsyyfZgVJ8s2dPs8Ue/kzY+K/GpcPXwAOiHHbA7jFy/42rgy2rCaPlwJ9/P+WwMHx0t6b9cRVRplevsXA6CTv7yDga6C3T9fzfRyaAt18fmrn7RYDXRLI+DhwVYzZTOB3gfuxwBNJ4rkC+Ajo6PPOa8BN/ll//14O9rL9DZgTE6+Z3l1uvLjWsvxH8t+FuPLaGWgGTAem+WdFuDrvLO/PbsCBgbr0Jm/2ViA+1ZXVR4DHvL0ewJfE1GXAU8DIZHmnQfV4ROQoXGI9pqrzcZXr2THW/qOq81V1C/AfYIuqTlXVHcCjuIQCl8mb4Srkbar6Mk6RnFUDkcap6npV/R8wC6dAwsSjI3AkcI2qblHVBbheTrwufiK+Bm5T1TJVfRT4BDgh8Hyyqi5R1R9wGSEi25nA06o6U1XLgFtxyumIBOGEsT9RVb9S1W+BGYGwyoA9gSIv51z1OS+GE3FDi5NVdbuqvgf8G/gFgKreiys8b3r/fptA1jJc/mjn07Wm82yJ4nEGLj0XqupmnAKoCZ+r6r0+Dz7o49Cmhn5Ux9HA+6r6fbyHItILl+f+lsD9ObgKugiXl58XkRaRh6r6kapehatcb8Apj+Ui8oiINE/gZxnQXUSaq+o6VX03kfC+zL7h3/8KYBLQN5F94GJgkqq+qao71M2pbcWV6x24Sr67iGSp6gpV/TSBPy1wlW+QB3FDa4jroZ9D9cNsd6jqFz7v/JmKeuQc4AFVfVfdHMd1uNGM4oDbsar6rS+rYamu/Ec4B/g/Vf1MVTf68IeIG9Y7G3hRVR/2/qz1dVGEdrje4OOq+jtvlrCsikgGcBrwe1XdpKofET/dvsele0IalOLBtWpeUNU1/v4hYobbcK22CD/EuW/m/7cDvlDV8sDzz4H2NZBnVeD/5oDf1dEO+Damkqhp2F/GVOKfe3+rk62dtwuAj/8XScIOYz9RWLfgFMYLIvKZiCSaVCwCevvu/npxw5PnAG0Ddu7FtaD+5gtwPK4GBHhL3NBglUniakiWZl8EngX/18hfr7ggZF6RiqHK9bhJ2WsD6fTfgNXoMJuInBMYxnrWV55/B0ap6vZ44ajqa6r6g6puVtWxwHpcDzXW3g7gQ+B94FvcO0k0/3Oal+tzcUOghyeJ594+rqtEZANwM9A6kX1cnrkiJs90xDU6lgGjcQrya68c2yXwZx1QEGM2HdhTRA7DKdg84OkkskDlPBEsi7HlZyOwlsrlp6b5Caov/xEqhe//Z+IaPh1xjfdEnIBrZN4dMEtWVnf3fsemRSwFuPyVkAajeEQkF9fy7Osz5ypcV/MAETmgFl5+BXSUynMOnXBdQ3DDS3mBZ8FKsDqqO9L7K6CViAQzfDDsMLQXEYlx/1UId1/hMg/gBopxGTASdqzs1dlPiKp+r6pXqGpn3EKQyxOMRX8BzFbVFoGrmape5sNsBtwG3A/cICKtEoS3SlUvUtV2wCXA36VuVgeVAh0C9x3rwM9QqOqJkTTBrbwaF0ijEwNWo4pHVf+pFYsEjgOaA72AR325edu7WSkiVZRLJGicEgfcOxA3X/Uy8C6u4jxTVXuo6toEsr+tqoOBPYAncD3viN+x3AV8DHRV1ea4IXCJYy/CF8CfY/JMnqo+7MN+SFUjIyQK/CWBPx/ghpODcm8G/oUbgTgPeERVtyWRBSrniWBZjC0/+bjhq2D5qc0nAMKW/0rhe3vbcQ3yL4iZg4nhXtxilWe83JC8rH7j/Y5Ni1j2xTVcEtJgFA/wc1wXujtuCORAXATmUrMhqghv4lq1V4tIlrjJ6pNwY5QAC4BTRSTPV17DauD3aqCD+MnRWFT1C2AeMFZEcsRNig7DzR2EZQ9gpJf9F7i0SDixHOAx4AQRGSBupdIVuCGKeQHZO9fAfkLETf7u5QvId7j3Vx7H6n+BvUXkPB+fLBE5RCoWMdwOvKOqw3Etz7vj+IGI/EJEIgpiHa5AR8KLjVdNeAwYKm4xSh5wfS39SQkiUgI0VdXFCax8h2v5Huiv4715T+BNEekkIkeKSLbPj1fhehuvef8H4SqwM3FDYO1VdYSqvk0CvF/niEihH6LdQOV3sZsEJrlxreANwEYR2Qc3FxYk9v3dC1wqIr39hHe+iJwgbiK8m4j0F7dyagtupCNevgNXZuIN6T3o43sa4Vaz/UpEOvhG0W9xw/rg5h2HisiBXp6bgTf9cGIiwuTVsOX/YWCMiJT4BtzNwKO+5/tPYKCInCEimeIWBhwY4/7XuGG8Gb7xn7Cs+t7wdFzjME/coozYBWDtcfNZb5CEhqR4zseNs//Pt2xXqeoq4A7gHKnhUkTfgjkJt6plDW4o4peq+rG38ldgGy4TPIh7SWF5GVgIrBKRNQnsnIUbU/8KNxf1B1V9sQZhvAl09bL/GTg9UcsziKp+ApyLG+tfg0uDkwIturHA73w3+soQ9pPRFXgR2Ai8DvxdVWfFkel74GfAEFx6rMK1UJuKyGDc5H6kIrocOFjcsuBYDsFVpBtxE5ijVPUz/+wG4EEfrzPiuE2Iqj4LTMTNfSyjotAkGvJLNyeQpNGhjmCZ+cY/Wu3fYwGux7EO1xIfBBwXyE+fAPuo25/0aJKhzljOA1b4obNLcUMy+DL2MPCZfx/tgCtxcw7f45TKozF+3UDg/anqO7gFG3d4uZfhJt3Bze+Mw+XXVbhK+roEafMu8J2I9I55NAensFcmU7ABHsItyvgMN3x1k/f/RVxD5d+4nnMXXD5PRqW4JrATtvw/gFupNwe30GoL8Bsv2/9wjZArcMOmC4BKo0d+OO9iYCXwJG7eLm5Z9U5+jRtGXoVboDA5Rp6zgQery0OiceeCjfpE3JLH4X4owUgzvif2Ea6XEXfOZCf8vgFAVW+ogZtncJPbYXq8Rgwi8jNghKr+vJbuV+DKY00ajj86fI/vfeBoVf06md2G1OMxjHpDRE4Rd+RHS1wLb0ZdK52d4BVcb8yoBar6Qm2VjhEeVd2qqvtUp3QghYpHRB4Qt6Hpo4BZK3Ebk5b635beXMRtXlomboPfwQE353v7S6WW52MZRgguwS1h/RQ3VxU7B1FXvOKv0KjqeK3ZUlzDaNCkbKhNRI7Gjf1PVdUe3mw8bpnxOHFLb1uq6jUicjxuXPJ43Iax21W1t5/Iewe3YkeB+UBPVV2XEqENwzCMlJOyHo+qzsFNaAUZTMUKkgdxK9ki5lP9JOkbQAsR2RM4FpipbvPVOtwO4EGpktkwDMNIPek+tK6Nqpb6/6uo2N3dnsqbklZ6s0TmVRCRi3GrM8jPz++5zz771KHYhmEYjZ/58+evUdVqP2uws9TbaamqqiJSZ+N8qnoPcA9Ar1699J130nLCuWEYRqNBROKdRFDnpHtV22o/hIb/jax++JLKu2E7eLNE5oZhGMYuSroVz1NU7HQ9H7dhKWL+S7+67TDgOz8k9zzwMxFp6VfA/cybGYZhGLsoKRtqE5GHcQfwtRb3XZQ/4HYbPybumzSf485mA7cr+3jc7uTNuCO5UdVvReRPVJw9daO602ENwzCMXZRGeXKBzfEYhmHUHBGZr6q9Uh2OnVxgGIZhpBVTPIZhGEZaqbfl1EbqKS8vZ82aNaxfv54dO3bUtziG0SDJycmhQ4cOZGUl+t6dUdeY4mnErFy5EhGhuLiYrKwsRJJ9d8swfnyoKmvXrmXlypWUlJTUtzg/GmyorRGzadMm2rdvT3Z2tikdw4iDiLDbbruxZcuW+hblR4UpnkZOkyb2ig0jGdYoSz9WKxmGYRhpxRSPYRiGkVZM8Ri7JCtWrEBEePXVV+tblDplypQpZGbamh+jcWOKx2iQXHDBBYgIIkJmZiZFRUVceumlrF27tk7DefXVVxERVqxYUaf+1pYzzzyTL79M/Tm4mZmZTJkypcbuhg8fTr9+/epcnpqyatUq2rZti4iwcuXKqHmkQRJ7/e53v6tHaY1YTPEYoVi1Cl5/3f2miz59+lBaWsqKFSuYOHEi//73v/nlL3+ZPgHqgdzcXNq0aVO9xUbIlClTQim18vJyzjnnHA499NCEdp588klKS0uj17XXXluHkho7iykeo1qmTYPOnWHQIPc7bVp6ws3OzqZt27Z06NCBwYMHM3r0aJ577jl++OGHqJ2vvvqKE088kby8PDp37lylFV9aWsqQIUNo0aIFubm59OvXj8g5fitWrKBPnz4AlJSUICLRik9VufXWW+ncuTPZ2dl06dKF2267rZLfxcXF/P73v2fUqFG0atWKNm3aMGbMGLZv3x618+qrr3LkkUdSUFBAQUEBBxxwAM8/n/iA9dihtsj9a6+9xsEHH0xeXh49e/bk7bffTugHwMKFCzn22GNp0aIF+fn57LvvvkzzL664uJgdO3YwdOjQaI8AYN26dZx77rl06tSJ3NxcunXrxoQJE4ic53jDDTdw//33M3v27Ki7SHpv3LiRUaNG0b59e/Ly8jjooIOYPn16Uhlry5/+9Ceys7MZM2ZMQjutWrWibdu20atZs2YpkcWoJara6K6ePXuqobpo0aKd9qO0VDU3VxUqrtxcZ55Kzj//fB0wYEAlswkTJiigGzZs0OXLlyugJSUl+uijj+rSpUv1uuuu04yMDP3kk09UVbW8vFwPPfRQPeCAA3Tu3Ln6wQcf6BlnnKEtWrTQb775Rrdv365PPvmkAvrWW29paWmprl27VlVV77jjDs3JydFJkybpkiVL9K677tKmTZvqfffdF5WnqKhIW7RooWPHjtUlS5boo48+qpmZmVE7ZWVl2rJlSx0zZowuWbJElyxZotOnT9c5c+YkjPfkyZM1IyOj0r2IaJ8+fXTOnDm6ePFiHTRokBYXF2tZWVlCf/bbbz8966yzdOHChfrpp5/qM888ozNmzFBV1a+//lozMjL0tttu09LSUi31L7O0tFTHjh2r8+fP188++0ynTZum+fn5+sADD6iq6vfff69nn322Hn744VF3mzdv1vLycu3Xr5/27dtX586dq59++qlOmjRJs7Ky9MUXXwz9zidPnqx9+/ZNaufll1/W9u3b6+rVq3XWrFkK6BdffBF9HskXHTt21N1220179uypEyZM0G3btiX1ty7KSmMAeEfTUEfXu5JIxWWKx1EXhWnePNXmzSsrnsJCZ55KYhXPwoULtXPnztq7d29VrahgJkyYELWzfft2bdasmd59992qqvriiy8qoAsXLoza2bJli7Zt21b/+Mc/qqrq3LlzFdDly5dXCr9Dhw561VVXVTIbPXq0lpSURO+Lior0pJNOqmRn0KBBOmTIEFVV/fbbbxXQWbNmhY53PMUD6Pz586Nmb7zxhgL68ccfJ/SnefPmOnny5ITPMzIykj6PMHLkSB04cGD0ftiwYVWUw6xZs7Rp06a6fv36SuZDhw7VwYMHVxtGhOoUz6pVq7Rdu3Y6c+bMaLixiuebb77R8ePH62uvvabvvfee3n777dq8eXM999xzk4ZtiseRLsVjy2eMpJSUQFlZZbNt25x5qnnllVdo1qwZO3bsYOvWrQwYMIBJkyZVsnPggQdG/2dkZLDHHnuwevVqwA037bbbbnTv3j1qp2nTpvTu3ZuFCxcmDHfDhg2sXLmSo48+upJ53759uf3229m8eTN5eXlVwgdo164dy5cvB6Bly5YMHz6cY489lv79+9O3b19OOeUUunXrVqN0EBEOOOCASmEArF69OqFfV155JcOHD4/Om5x88skcfPDBScMpLy9n/PjxPPLII6xcuZItW7ZQVlZGUVFRUndvv/0227Zto3379pXMt23bRteuXRO6mzt3Lscdd1z0fvv27ZSVlVUaFuvTpw/PPvssAOeccw6//OUvGThwYEI/W7duzVVXXRW9P/DAAykoKODCCy9k3LhxVWQ06geb4zGS0rYtTJoEublQWOh+J01y5qmmd+/eLFiwgMWLF7NlyxZmzpxJ586dK9nJzs6udC8ilJeXp164kOHfe++9zJ8/n2OOOYbZs2fTo0ePKsqzOpo0aUJGRkalMICk8bz++utZsmQJZ5xxBh999BGHHXZYtSu7JkyYwNixYxk5ciQzZ85kwYIFDB8+nG3btiV1V15eTmFhIQsWLKh0LVq0KKo04tGrV69K9m+88cYqZvfdd1/U/ksvvcQtt9xCZmYmmZmZDBgwAHBzVpdccknCcI444ggAPv/886TxMNKH9XiMajnvPDjmGFi+3PV00qF0wK3w2muvvWrt/ic/+Qlr165l0aJF0V7P1q1befPNNxkxYgRQoTiCp3c3b96cDh06MGfOHE488cSo+ezZsykpKYn2dsLSo0cPevToweWXX86ll17KPffck7SirCs6d+7MiBEjGDFiBOPGjeOWW27hpptuAly8Y08snzNnDoMGDeLCCy+Mmi1durSSnXjuevXqxfr169myZQs9evQILV/s+91jjz2SvvMPP/yw0v3bb7/NhRdeyPPPP8++++6bMJx3330XgA4dOoSWzUgtpniMULRtmz6FU1f079+fQw89lLPPPps777yTwsJC/vSnP7FlyxYuu+wyAIqKimjSpAnPPPMMZ555Jk2bNqWwsJDrrruOK664gq5du9KvXz9efvll7rrrLu68887Q4S9btox7772Xk046iY4dO/LVV18xd+7caoe8dpaNGzdyzTXXcNppp1FSUsL69et57rnnKg05lpSUMGvWLI477jiys7Np3bo13bp1Y9q0acyaNYv27dszdepU3nzzTVq2bFnJ3eOPP87ChQtp06YNBQUF9O/fn4EDB3Lqqacyfvx49t9/f9atW8e8efPIycnhoosuqpN4xSq1NWvWANCtW7fo8OOUKVPIyMjg4IMPJicnh7lz53LVVVdx+umn06lTpzqRw9h5bKjNaLSICE888QT77LMPJ5xwAocccgirVq1i5syZtG7dGoA2bdowduxYxo0bx5577sngwYMBuOyyy7jxxhu5+eab6d69O3/5y18YN24cw4YNCx1+fn4+S5cuZciQIey9996cdtppHHHEEdxxxx0piW+EzMxM1q1bx7Bhw9h333059thjadOmDQ899FDUzoQJE5g/fz7FxcXsvvvugBue69u3L4MHD+bwww9n3bp1jBw5spLfw4YN45BDDuGII45g99135+GHH0ZEeOqppzj11FMZM2ZMNL2ffvppunTpktK4xtKkSRPGjx9P79692X///bn11lu56qqr+Mc//pFWOYzkiFvI0Ljo1auXRvZq/JhZvHhx0iEIwzAcVlYcIjJfVXulOhzr8RiGYRhpxRSPYRiGkVZM8RiGYRhpxRSPYRiGkVZM8RiGYRhpxRSPYRiGkVZM8RiGYRhpxRSPYRiGkVZM8RiGYRhpxRSPsUuyYsUKRIRXX321vkWpU2K/QGoYjRFTPEaD5IILLoh+XjkzM5OioiIuvfRS1q5dW6fhvPrqq4gIK1asqFN/a8uZZ57Jl19+mfJwMjMzq3wmPAzDhw+Pfh481ZSXlzNgwABEpMpZa8XFxdH8EbmOOuqotMhl7DzWtDLCsWpV2r+L0KdPHx577DG2b9/O/PnzGT58OF988QVPP/10WsKvD3Jzc8nNza1vMRoEN954I/n5+QmfX3PNNYwePTp6H/ttJKPhYj0eo3qmTYPOnWHQIPc7bVpags3OzqZt27Z06NCBwYMHM3r0aJ577jl++OGHqJ2vvvqKE088kby8PDp37lylFV9aWsqQIUNo0aIFubm59OvXj8gBsitWrKBPnz6AO+5fRKKteVXl1ltvpXPnzmRnZ9OlSxduu+22Sn4XFxfz+9//nlGjRtGqVSvatGnDmDFj2L59e9TOq6++ypFHHklBQQEFBQUccMABPP/88wnjHDvUFrl/7bXXOPjgg8nLy6Nnz568/fbbSdNu4cKFHHvssbRo0YL8/Hz23Xdfpvn3VlxczI4dOxg6dGi0twCwbt06zj33XDp16kRubi7dunVjwoQJRA4SvuGGG7j//vuZPXt21F0kvTdu3MioUaNo3749eXl5HHTQQUyfPj2pjMl4+eWXeeCBB5g8eXJCO82aNaNt27bRq1WrVrUOz0gz6fi+drqvnj171vRT442SOvmOfGmpam6uKlRcubnOPIWcf/75OmDAgEpmEyZMUEA3bNigy5cvV0BLSkr00Ucf1aVLl+p1112nGRkZ+sknn6iqanl5uR566KF6wAEH6Ny5c/WDDz7QM844Q1u0aKHffPONbt++XZ988kkF9K233tLS0lJdu3atqqrecccdmpOTo5MmTdIlS5boXXfdpU2bNtX77rsvKk9RUZG2aNFCx44dq0uWLNFHH31UMzMzo3bKysq0ZcuWOmbMGF2yZIkuWbJEp0+frnPmzEkY78mTJ2tGRkalexHRPn366Jw5c3Tx4sU6aNAgLS4u1rKysoT+7LfffnrWWWfpwoUL9dNPP9VnnnlGZ8yYoaqqX3/9tWZkZOhtt92mpaWlWurfZWlpqY4dO1bnz5+vn332mU6bNk3z8/P1gQceUFXV77//Xs8++2w9/PDDo+42b96s5eXl2q9fP+3bt6/OnTtXP/30U500aZJmZWXpiy++GPqdR1i1apW2b99eZ82apaqqgE6bNq2SnaKiIm3Tpo22atVKu3fvrr/5zW90zZo1NQ4rQp2UlUYA8I6moY5O/ABmAE8lutIhXG0vUzyOOilM8+apNm9eWfEUFjrzFBKreBYuXKidO3fW3r17q6pGFc+ECROidrZv367NmjXTu+++W1VVX3zxRQV04cKFUTtbtmzRtm3b6h//+EdVVZ07d64Cunz58krhd+jQQa+66qpKZqNHj9aSkpLofVFRkZ500kmV7AwaNEiHDBmiqqrffvutAtEKNAzxFA+g8+fPj5q98cYbCujHH3+c0J/mzZvr5MmTEz7PyMhI+jzCyJEjdeDAgdH7YcOGad++fSvZmTVrljZt2lTXr19fyXzo0KE6ePDgasMIsmPHDh0wYIBef/31UbN4iufWW2/VmTNn6gcffKCPPfaY7r333rr33nvr5s2baxReBFM8jnQpnmRzPLf631OBtkBkdu8sYPXO9LJEZAwwHFDgQ2AosCfwCLAbMB84T1W3iUhTYCrQE1gLnKmqK3YmfKMGlJRAWVlls23bnHmKeeWVV2jWrBk7duxg69atDBgwgEmTJlWyc+CBB0b/Z2RksMcee7B6tcueCxcuZLfddqv05c2mTZvSu3dvFi5cmDDcDRs2sHLlSo4++uhK5n379uX2229n8+bN0c9fB8MHaNeuHcuXLwegZcuWDB8+nGOPPZb+/fvTt29fTjnlFLp161ajdBARDjjggEphAKxevTqhX1deeSXDhw9nypQp9OvXj5NPPrnaL5+Wl5czfvx4HnnkEVauXMmWLVsoKyujqKgoqbu3336bbdu20b59+0rm27Zto2vXrgndNWvWLPq/T58+PPvss9x8881s3bqVP/zhD0nDvOKKK6L/99tvP3r27Mlee+3Ff/7zH84+++ykbo36J+Ecj6rOVtXZwJGqeqaqzvDX2UCf2gYoIu2BkUAvVe0BZABDgL8Af1XVvYB1QORTj8OAdd78r96ekS7atoVJkyA3FwoL3e+kSWlZYNC7d28WLFjA4sWL2bJlCzNnzqRz586V7MROKIsI5eXlKZctbPj33nsv8+fP55hjjmH27Nn06NGjivKsjiZNmpCRkVEpDCBpPK+//nqWLFnCGWecwUcffcRhhx3G7373u6ThTJgwgbFjxzJy5EhmzpzJggULGD58ONu2bUvqrry8nMLCQhYsWFDpWrRoEc8++2xCd0G79913HwAvvvgi8+bNo2nTpmRmZkbnu84//3z22WefhH517tyZNm3aNJjViUZywqxqyxeRzqr6GYCIlACJl5qEDzdXRMqAPKAU6A9EmioPAjcAdwGD/X+AfwF3iIj4bqGRDs47D445Ju2r2nJzc9lrr71q7f4nP/kJa9euZdGiRdFez9atW3nzzTcZMWIEUKE4duzYEXXXvHlzOnTowJw5czjxxBOj5rNnz6akpCTa2wlLjx496NGjB5dffjmXXnop99xzD5dcckmt4xWWzp07M2LECEaMGMG4ceO45ZZbuOmmmwAX72CcAebMmcOgQYO48MILo2ZLly6tZCeeu169erF+/Xq2bNlCjx49QssX791OnjyZTZs2VTLbb7/9+POf/8xpp52W0K8vv/ySr7/+mo4dO4YO36g/wiieMcArIvIZIEARUOtSo6pfisitwP+AH4AXcENr61U1shxoJRDpt7cHvvBut4vId7jhuDVBf0XkYuBigE6dOtVWPCMRbdumTeHUFf379+fQQw/l7LPP5s4776SwsJA//elPbNmyhcsuuwyAoqIimjRpwjPPPMOZZ55J06ZNKSws5LrrruOKK66ga9eu9OvXj5dffpm77rqLO++8M3T4y5Yt49577+Wkk06iY8eOfPXVV8ydO7faIa+dZePGjVxzzTWcdtpplJSUsH79ep577rlKQ44lJSXMmjWL4447juzsbFq3bk23bt2YNm0as2bNon379kydOpU333yTli1bVnL3+OOPs3DhQtq0aUNBQQH9+/dn4MCBnHrqqYwfP57999+fdevWMW/ePHJycrjoootCy16SYAi3Q4cO0WG7119/nddee43+/fuz22678fHHH3PttdfSqVMnTjnllFqmmpFOql1OrarPAV2BUbghsm6qmng9aDWISEtcL6YEaIfrPQ2qrX8BOe9R1V6q2mv33XffWe+MRoCI8MQTT7DPPvtwwgkncMghh7Bq1SpmzpxJ69atAWjTpg1jx45l3Lhx7LnnngwePBiAyy67jBtvvJGbb76Z7t2785e//IVx48YxbNiwZEFWIj8/n6VLlzJkyBD23ntvTjvtNI444gjuuOOOlMQ3QmZmJuvWrWPYsGHsu+++HHvssbRp04aHHnooamfChAnMnz+f4uJiIuXl+uuvp2/fvgwePJjDDz+cdevWMXLkyEp+Dxs2jEMOOYQjjjiC3XffnYcffhgR4amnnuLUU09lzJgx0fR++umn6dKlS53Hr2nTpkyfPp2BAwey9957M2LECA477DBef/31SvNGRsNFqhuxEpE84HKgSFUvEpGuOOXz31oFKPILYJCqDvP3vwQOB34BtPW9msOBG1T1WBF53v9/XUQygVXA7smG2nr16qWRvRo/ZhYvXsy+++5b32IYRoPHyopDROaraq9UhxNmA+lkYBtOOQB8Cdy0E2H+DzhMRPLEzZIOABYBs4DTvZ3zgSf9/6f8Pf75yza/YxiGsesSRvF0UdXxQBmAqm7GzfXUClV9E7dI4F3cUuomwD3ANcDlIrIMN4dzv3dyP7CbN78cuLa2YRuGYRj1T5jFBdtEJBe35wYR6QJs3ZlAVfUPQOxC/c+AQ+PY3YIbhjMMwzAaAWEUzw3Ac0BHEfkncCRuw6dhGIZh1JhqFY+qviAi84HDcENso1R1TTXOjAaCqkY3HBqGURWbMk4/1c7xiMhLqrpWVZ9W1f+q6hoReSkdwhk7R1ZWVqWTnA3DqEpZWZl9fC/NJFQ8IpIjIq2A1iLSUkRa+auYis2dRgNmjz324Msvv2Tz5s3WqjOMOJSXl7N69WoKCwvrW5QfFcnU/CXAaNwmz/lUrGTbAKR2B5xRJzRv3hxw36wpiz3o0zAMwG30jWwoNtJDQsWjqrcDt4vIb1T1b2mUyahDmjdvHlVAhmEYDYEw+3jKRaRF5MYPu41InUiGYRhGYyaM4rlIVddHblR1HRD+1D/DMAzDCBBG8WRIYD2uiGQA2UnsG4ZhGEZCwqwhfA54VEQiX6+6xJsZhmEYRo0Jo3iuwSmby/z9TOC+lElkGIZhNGrCnFxQLiJTcKdCf5J6kQzDMIzGTJiTC04GFuCH10TkQBF5KsVyGYZhGI2UMIsL/oA7NXo9gKouwH091DAMwzBqTBjFU6aq38WY2fkrhmEYRq0Is7hgoYicjVtW3RUYCcxLrViGYRhGYyVMj+c3wE9wH397CPgOGJVKoQzDMIzGS5gezwmq+lvgtxEDEfkF8HjKpDIMwzAaLWF6PNeFNDMMwzCMaknY4xGR44DjgfYiMjHwqDmwPdWCGYZhGI2TZENtXwHvACfjvscT4XtgTCqFMgzDMBovyb7H8z7wvog85O11spMLDMMwjJ0lzBzPIOzkAsMwDKOOCKN4bsBOLjAMwzDqCDu5wDAMw0grdnKBYRiGkVZqenLBw8AGYHQKZTIMwzAaMWG+x7MZd2rBb6uzaxiGYRjVkWwD6QySzOWo6skpkcgwDMNo1CTr8dzqf08F2gL/8PdnAatTKZRhGIbReEm2gXQ2gIhMUNVegUczROSdlEtmGIZhNErCLC7IF5HOkRsRKQHyUyeSYRiG0ZgJs5x6DPCKiHwGCFAEXJxSqQzDMIxGS5hVbc/5/Tv7eKOPVXVrasUyDMMwGithhtpQ1a2q+r6/dlrpiEgLEfmXiHwsIotF5HARaSUiM0Vkqf9t6e2KiEwUkWUi8oGIHLyz4RuGYRj1RyjFkwJuB55T1X2AA4DFwLXAS6raFXjJ3wMcB3T118XAXekX1zAMw6gr0q54RKQQOBq4H0BVt6nqemAw8KC39iDwc/9/MDBVHW8ALURkz7QKbRiGYdQZ1SoeETlSRPL9/3NF5P9EpGgnwiwBvgEmi8h7InKf97+NqpZ6O6uANv5/e+CLgPuV3ixWzotF5B0Reeebb77ZCfEMwzCMVBKmx3MXsFlEDgCuAD4Fpu5EmJnAwcBdqnoQsImKYTUAVFWp4QnYqnqPqvZS1V677777TohnGIZhpJIwime7VwSDgTtU9U6gYCfCXAmsVNU3/f2/cIpodWQIzf9+7Z9/CXQMuO/gzQzDMIxdkDCK53sRuQ44F3haRJoAWbUNUFVXAV+ISDdvNABYBDwFnO/Nzgee9P+fAn7pV7cdBnwXGJIzDMMwdjHCbCA9EzgbGKaqq0SkE3DLTob7G+CfIpINfAYMxSnBx0RkGPA5cIa3+wxwPLAM2OztGoZhGLso4kbRGhe9evXSd96x4+QMwzBqgojMjzmbMyWEWdV2qt/U+Z2IbBCR70VkQ6oFMwzDMBonYYbaxgMnqeriVAtjGIZhNH7CLC5YbUrHMAzDqCvC9HjeEZFHgSeA6Dltqjo9VUIZhmEYjZcwiqc5bjXZzwJmCpjiMQzDMGpMmM8i2PJlwzAMo84Is6qtg4j8R0S+9te/RaRDOoQzDMMwGh9hFhdMxp0e0M5fM7yZYRiGYdSYMIpnd1WdrKrb/TUFsFM4DcMwjFoRRvGs9Z9DyPDXucDaVAtmGIZhNE7CKJ4LceemrQJKgdOBC1Iok2EYhtGICbOcuoOqnhw0EJEjqfxxNsMwDMMIRZgez99CmhmGYRhGtSTs8YjI4cARwO4icnngUXMgI9WCGYZhGI2TZENt2UAzbyf4xdENuHkewzAMw6gxCRWPqs4GZovIFFX9XESaefONaZPOMAzDaHSEWVxQICLvAa0ARGQNcL6qfpRSyQzDMIxGSZjFBfcAl6tqkaoWAVd4M8MwDMOoMWEUT76qzorcqOorQH7KJDIMwzAaNWGG2j4TkeuBaf7+XOCz1IlkGIZhNGbCnlywO+77O9P9/wtTKZRhGIbReAnzPZ51wMg0yGIYhmH8CKhW8YjI3sCVQHHQvqr2T51YhmEYRmMlzBzP48DdwH3AjtSKYxiGYTR2wiie7ap6V8olMQzDMH4UJDurrZX/O0NERgD/AbZGnqvqtymWzTAMw2iEJOvxzAcUEH9/VeCZAp1TJZRhGIbReEl2VltJOgUxDMMwfhyEWdWWA4wAjsL1dOYCd6vqlhTLZhiGYTRCwiwumAp8T8XH387GnWLwi1QJZRiGYTRewiieHqraPXA/S0QWpUogwzAMo3ET5sicd0XksMiNiPQG3kmdSIZhGEZjJkyPpycwT0T+5+87AZ+IyIeAqur+KZPOMAzDaHSEUTyDUi6FYRiG8aMhzCGhn6dDEMMwDOPHQZg5npQgIhki8p6I/Nffl4jImyKyTEQeFZFsb97U3y/zz4vrS2bDMAxj50moeESkaYrDHgUsDtz/Bfirqu4FrAOGefNhwDpv/ldvzzAMw9hFSdbjeR1ARKYlsVMrRKQDcALuxGtERID+wL+8lQeBn/v/g/09/vkAb98wDMPYBUk2x5MtImcDR4jIqbEPVXX6ToR7G3A1UODvdwPWq+p2f78SaO//twe+8GFuF5HvvP01QQ9F5GLgYoBOnTrthGiGYRhGKkmmeC4FzgFaACfFPFPcZ7BrjIicCHytqvNFpF9t/IiHqt4D3APQq1cvrSt/DcMwjLol2SGhrwKvisg7qnp/HYZ5JHCyiBwP5ADNgduBFiKS6Xs9HYAvvf0vgY7AShHJBAqBtXUoj2EYhpFGwqxqmyYiI0XkX/76jYhk1TZAVb1OVTuoajEwBHhZVc8BZgGne2vnA0/6/0/5e/zzl1XVejSGYRi7KGEUz99xpxf83V8HA6n4Iuk1wOUisgw3hxPpZd0P7ObNLweuTUHYhmEYRpoIc3LBIap6QOD+ZRF5vy4CV9VXgFf8/8+AQ+PY2YKdhG0YhtFoCNPj2SEiXSI3ItIZ2JE6kQzDMIzGTJgez1W4TyF8hvsMdhEwNKVSGYZhGI2WMGe1vSQiXYFu3ugTVd2aWrEMwzCMxkqYHg9e0XyQYlkMwzCMHwH1dkioYRiG8ePEFI9hGIaRVqpVPCLyUhgzwzAMwwhDwjkeEckB8oDWItISt6IN3BE37RO5MwzDMIxkJFtccAkwGmgHzKdC8WwA7kitWIZhGEZjJdkhobcDt4vIb1T1b2mUyTAMw2jEhNnH8zcROQIoDtpX1akplMswDMNopFSrePwXSLsAC6g4KkcBUzyGYRhGjQmzgbQX0N0+RWAYhmHUBWH28XwEtE21IIZhGMaPgzA9ntbAIhF5C4ie0aaqJ6dMKsMwDKPREkbx3JBqIQzDMIwfD2FWtc1OhyCGYRipYtUqWL4cSkqgrU0c1Dthjsz5XkQ2+GuLiOwQkQ3pEM4wjIbFqlXw+uvud1dh2jTo3BkGDXK/06bVt0SpY1d5P9UqHlUtUNXmqtocyAVOA/6ecskaGLvKCzWMVLErVuCrVsEll8APP8CGDe73kkvqsRyHqEhqW9fsSu+nRqdTq+MJ4NjUiFP/xHvpu9IL3VVJVNgausIPJXcdRaI+06LBVeAhWb4csrIqm2VnO/O0E6IiqW1ds8u9H1VNegGnBq7TgXHA69W5q8+rZ8+eWhumTlXNzVVt3tz9Tp2qWlrq/kPFlZvrzOuU0lLVefNCe5zIeg29qVNqG3a8dE9mnmp5whJG7qFZU7UsK2Bp4sRaCVXfaTFvngs7WA4KC515SgKsCUkKw5oZ87Q4p3Snym+dRCVeRZKTozpjRtTj6uqaoByxMtX4/SQAeEfTUEeHUTyTA9e9wG+BPdIhXG2v2iieRPnilltUCwoqm3ctKNUPJiXPiYkySaUME7mZOLFyrZKocvL2H5tYmpKKukaJFSNfbNjBKCQruIkK2/vv11LhJ0ijxyZWCJGsACeJYo3lbkOpbiLGErgMlZur6yZOrT7sBJVnTJ0V13Gi7BXMF3HzY5xI16gBFpMZEsUzEaEq+gSRi4S1bmKFDGVZuTo0a6oWFlbNC9WRrEyFVkilpaqTJlWtSEA1Pz/qcTzlUVDgnAajmZXlrmDergsFq9qAFM+ueNVG8cR76ZGXF7w/l6m6iVzdURCnUMWp9IKZJPg/2hJu1qxqoIHKKVar7ChorpvI1XOZWm1FHaZyCk2SWixepRSJQmwhiVWGiVpqkyYlLoSJKupIZRObRsF3FqmE4hXgiLJMVlEnk7trQan+5+p5ulczVwEcxjxdR5xM5a9N5OpezUorVYbBvBPJI2X5Vd95sM6KVqQBwSPxjJe9IvkiGM8qPbM4kY5UwpEKPNZKaanq2zNKtTwnN2E8q+vwhWo8RSzFidwmcvWwvPerKPzynFx9e0ZpJYUUG0Bsg2TGDJdW1aVdFWUW9CiJrLEer506o4ryiLznRM7i5e1E7ycMDUbxAB2A/wBf++vfQId0CFfbq656PLFXl/yqLdhIoYoU3HiKIfZK2BKOc5Xn5Op7U9+PW5jbUBq3wgtbOcXWBKvfdz251e9XdP0rtR7jFJ7ynFydOr40bmMuQfmqUIYJWvPFOaW6dGpV86A+jlUSezWL/256ULUSCqZdmEJeRYHHkTtSAWxvVvH+q3vP6yjUw5gX2m08uaMVT378SjhRPIN5PV5Y5U1zdNEtM3T1+4l7iPHq1wF583R9jLKNxDPRO0yWF4qalurcW3yeTKQNYsIazqSqCr+w0LlN0DqL11jMy4sfTNLGaLA1k5Pj/sd6kMjj/Pyo8ghTnuK+N69gazss2JAUz0xgKG7PTyZwATAzHcLV9trZOZ54lU9Bgep/rp7nMldMRj+eGTWq3KprCQev7yjQW7Ou1u8oqBJusNLamptY4UXsbMttpuUxD8tBy/IKtKxJlm4lS9fj/Lm/31QtzinV03Nn6GYSF/TvKNBfZ09KGNd4heWnOfN0dOZE3eZb89sCLbWhWVOrmMcrhBF/IuHGS9NElVBsRRhW5uKcUp178VQtz6ksX9eC+Eqia4FvkGTnxh1mSaYYj2dGFbnXU6g/zZlXSa4wiq26+MVLu3LQ78nXTeTqSJmoA/Jc/IMd8OKcUh2QN0+LmpZG69ewCjNid0BeRV7YkhO/p7qe5rqFLN2ekaVlOXlV8nCYNNXcXKd4mleNZ1nT/Gobi8nyRthGpOIrkvHjkyrPsI25uPVIbSZ2AjQkxbMgjFlDumqreFQTN6pyc9W1unJjC1WOXs4tVRTDegp0OK5CbkOpHkZFJZkws8bReOWgG8ivUtg2kauH58evtEpyK8I5nuSKI9G1hSzdTI5uIHlBL8cpn0jBDTMssJ6qCnBbluvZbcvKrWIeKYSRdBzBRN1Erq7zSvJXTIwbz2RplKgiDL6noMwurBzdSlYV+RbdMqNKg2RHQaF+MGle3Lm8HQWF0fRKpDDjNWbKc3L1haml0bwZZigvGJdIQzt0fozzjodmTdX333cNhOA7iPTSgu9nW15hlQo91k68vJBQeSS5djQriKu0dhQUJl8pFFOWj2dG3DIbSfNY85o0IqMVSbCLGK/AFBbqmhnzkg5fZ2dr3AbPzq58akiK5yXgXCDDX+cCL6VDuNpeO6N4IkRadAPzK7f0dKpr8W7LK9QtZOkWsuIWnkhh3ULlnsTIJhN1YP48HSmRgucK5xvnVowdlefk6pacgoR+biJXnz93qs69peqwxnoK9Znr5+nQTFfwqlMcySqbZM+2No2vDKeOL9WJEyvSrkNmqXbMKtVTmyZXgOsp1C+unxQ3PkunzotWdNWndbZuyyvU8pzKE82ul1Ko27JydVTmxEqyDciLr8xOya4qc2zY6ynUtVNnuPmRgHlZdoIKIDAPWJxTqr+I06PcRI7+Im9GoFfo4hPJhJG8Gc9tpAcbjGdxjnsnM2ZoNF8E4zkgryI/fh+nkRP7jh+7vqpSiDRUIv6Oypyoa2ZUxDOYxvHeYeU0LdBbMqr28mPjGemRRcpOJKxInOMuIvDltyxO/o34GSnX0bhkTNRnrp+nozOq5pHqRgS2Z2RpedOcKu8wmhdmzNCyrMruI3kntg6Kt2AnmLer+F8LGpLiKQKeAr7xczxPAJ3SIVxtr7pQPNEMmt88bqGPl+ES9U5i7WzLLYhm3EjrKThpWZxTqr/OnpS0F1WcU6pvXB+/0vrk+qlVMnNt5ExU0H/FRP1VHPm+o0DX3zIpqjzL8ptreWaWlmdlaVnT6odI4lVom3A9odiKPbE/OXp67oxKFU+w8gsO8UVlSzQEmV1V5njK9oWppYEeQGG0Z5BE70QrjG15zX1Fl63rKIw2VLblNY8qj9hhrorGT2W3QUVSKZ45fkIlTiUXmx+r6yEnaiDE68FGJvPLc3LjpnGyvLcluzoFWNE7iZSdikq4eaWGR4ze0eKc0riNiqRlIT++otqWW1Dl/QXfx/lNpmpR09Kq71ArRlcqGgMxeSe2DoqjeYJ5O9b/2tBgFM+ueO2U4kky1rb6/dJo9zdeF3s9BXozyVtqsZVWcCikugnfiP1zmaqbSVxp7YgzeRmrOCoPWVUUGNcDy9GyJlkxYVcU9ETyRQpomLhH7Ed6cKMyJ+rSqa7CDBbCUZkTnTILuXqhphP2Nb0qWvYV8s29xS3sCA7DBJfcxy4IjL8QIkfPYWrc4cJImhfnlOraqVXzZph3o5B4UjtO/tpErn5H1V53WbZbQhnbEIjXEzw9t+pwYbK8EK8xFOzNBiv02PmYktzEC38SjbRF4lldIyxcWU48TBe8Ykfa4g1/di0o1UXjEyyiCCwVTbTCdWf2GzUYxQOUAP8HTPc9n6eAp9IhXG2vWiue2NwQ88KXXT0punIskWKoydh0dZO/kYLxQ05h0spzEzk6PLdqpRVr5/ScGVUKQ2SCtzinVO+7KbCqLTCkWN2ihXjDgskqmO25+ZWGgkZnTowu5S2LZx5vNR1oWW784b6aTNjXRObynFyde/HU6PBHRL4wy7eDUajNQohIY6Msp2rejM1HNZ53iONPZDXZ+psmxh/KmRoY4mladfVWmPSOHRb8TU7VXvSOfFfu7rupomUfr0JPlKaR+CRYW6BtqH4YuC7KcuSqbm1BNO9U00iITetImuzk2gJtSIrnfWAk8FOgb+RKh3C1vWqleEKsp042gVmWXXllVlm23/AQmQlMsKopshAg0VWcUxodK8/NVR2YX3VeZ0dBoS67elKVCW6F6HrqyLBDpNUdWeufdF9FzL6keJ2O4pzSGvVIgmuTE+37SLSEPFpy46ynru2Efdz0SiJzJFkS7VcJs5ghUaMlkdvqGjObyNUu+cn9r3LFyeubcCvxquwBqcF66rL85I2kuO8wSV4INuFje47BV1XdarpqVlNXbBmopswmyiORtIs4LSys+WrqeL22MFdQ6UXqi9p2exqS4nkzHYLU5VUrxRN2B2lMAU20Iz5uAY2p9ddNnJpwW0Ls/lHVxBVeTXeQ1sX+0Uqb1KrbQRopiSneQRpRkKGWNwcrmFgtHDeS1eeXRA2AeK3h2FVXVXbcF1Y0ZuI1NqIVYIJGRbJl3FV2QRZWliFhvqhud2eS/U3R1WVhd5BWswsy3oh4xfBq1V56bAMurvchymyiPBJMu3j7m6prl+XkqM69pep2jUR1UMK8Xc0m4OpoSIrnbOAPwOHAwZErHcLV9qqzHk+CM3MqLZWtTTgJjpsJ1QuJ5yD2zJyd2bpcuygkj0QyTZfoLJZanplTZdNrvEo1Uau92kjWTu7YHk+kUZHw6JZA2AkbG2EaFUlbCyHjWV1847iLzQo1OaKmRjIlCSuoYGIr5Joe4VPrPBKwFuLEnMR1UGwjIabBtG7i1Gp7i2FpSIpnLLASmA3M8tfL6RCuttdOz/EEC2gNCtzOUONeyE4WhpRQ27BTpUhTnRYh5A4OwYZqVNQ0rLDsbFrU8BTKdGbDZEWhrirkupAx1IBEsvecrMFUw/eTiIakeJYB2XUWIHT0ymsRsBAY5c1b+VMSlvrflt5cgIlejg/C9LZ2elVbvDWYaehJ/KhpiIo0DCHkrrMo1HejIg0NsDqnjirkuiB0NVLbsfA6eD8NSfE8UZenUQN7RpQHUAAsAboD44Frvfm1wF/8/+OBZ70COizMnFOd7OOJpaFXgIaRanbFBlgDU5gprUbq4P2kS/GICysxIvIKsD/wNrA1Yq6qJyd1GBIReRK4w1/9VLVURPYEXlHVbiIyyf9/2Nv/JGIvkZ+9evXSd955py7EMwwjyKpV7itqJSXQtm19SxOOadPcV9Gys2HbNpg0Cc47r76lSg07+X5EZL6q9kqBZJXIDGHnD6kKXESKgYOAN4E2AWWyCmjj/7cHvgg4W+nNKikeEbkYuBigU6dOqRLZMH7ctG276yicCOedB8ccs+spzNqwi7yfahWPqs4O3ovIUcBZuMUGtUZEmuE+sTBaVTeISDBMFZHkXbGqct4D3AOux7MzshmG0cjYRSrkHwthejyIyEG4ZdW/AJbjFEatEZEs78c/VXW6N14tInsGhtq+9uZf4hYkROjgzQzDMIxdkCaJHojI3iLyBxH5GPgb8D9AVPWnqnpHbQMU17W5H1isqv8XePQUcL7/fz7wZMD8l+I4DPgu2fyOYRiG0bBJ1uP5GJgLnKiqywBEZEwdhHkkcB7woYgs8Gb/DxgHPCYiw4DPgTP8s2dwK9uWAZtxH6UzDMMwdlGSKZ5TgSHALBF5DngEt6R5p1DVV5P4MyCOfQV+tbPhGoZhGA2DhENtqvqEqg4B9sFt+BwN7CEid4nIz9Ikn2EYhtHISKh4IqjqJlV9SFVPwk3svwdck3LJDMMwjEZJtYoniKquU9V7VLXKkJhhGIZhhKFGiscwDMMwdhZTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZa2WUUj4gMEpFPRGSZiFxb3/IYhmEYtWOXUDwikgHcCRwHdAfOEpHu9SuVYRiGURt2CcUDHAosU9XPVHUb8AgwuJ5lMgzDMGpBZn0LEJL2wBeB+5VA76AFEbkYuNjfbhSRT3YivNbAmp1wn0pMttrRkGWDhi2fyVY7GrJsEF++onQEvKsonmpR1XuAe+rCLxF5R1V71YVfdY3JVjsasmzQsOUz2WpHQ5YN6le+XWWo7UugY+C+gzczDMMwdjF2FcXzNtBVREpEJBsYAjxVzzIZhmEYtWCXGGpT1e0i8mvgeSADeEBVF6YwyDoZsksRJlvtaMiyQcOWz2SrHQ1ZNqhH+URV6ytswzAM40fIrjLUZhiGYTQSTPEYhmEYacUUT4CGfCyPiDwgIl+LyEf1LUssItJRRGaJyCIRWSgio+pbpggikiMib4nI+162P9a3TLGISIaIvCci/61vWWIRkRUi8qGILBCRd+pbniAi0kJE/iUiH4vIYhE5vL5lAhCRbj69ItcGERld33JFEJExvix8JCIPi0hO2mWwOR6HP5ZnCXAMboPq28BZqrqoXgXziMjRwEZgqqr2qG95gojInsCeqvquiBQA84GfN4S0ExEB8lV1o4hkAa8Co1T1jXoWLYqIXA70Apqr6on1LU8QEVkB9FLVBrcRUkQeBOaq6n1+tWueqq6vZ7Eq4euVL4Heqvp5A5CnPa4MdFfVH0TkMeAZVZ2STjmsx1NBgz6WR1XnAN/WtxzxUNVSVX3X//8eWIw7baLeUcdGf5vlrwbT2hKRDsAJwH31LcuuhIgUAkcD9wOo6raGpnQ8A4BPG4LSCZAJ5IpIJpAHfJVuAUzxVBDvWJ4GUXnuSohIMXAQ8GY9ixLFD2UtAL4GZqpqg5ENuA24GiivZzkSocALIjLfH0vVUCgBvgEm+2HK+0Qkv76FisMQ4OH6FiKCqn4J3Ar8DygFvlPVF9Ithykeo84QkWbAv4HRqrqhvuWJoKo7VPVA3IkXh4pIgxiqFJETga9VdX59y5KEo1T1YNzJ8L/yQ74NgUzgYOAuVT0I2AQ0tHnZbOBk4PH6liWCiLTEjeSUAO2AfBE5N91ymOKpwI7l2Qn8/Mm/gX+q6vT6licefihmFjConkWJcCRwsp9HeQToLyL/qF+RKuNbyKjq18B/cEPSDYGVwMpA7/VfOEXUkDgOeFdVV9e3IAEGAstV9RtVLQOmA0ekWwhTPBXYsTy1xE/g3w8sVtX/q295gojI7iLSwv/PxS0e+bhehfKo6nWq2kFVi3H57WVVTXvrMxEiku8Xi+CHsX4GNIhVlaq6CvhCRLp5owFAvS9mieEsGtAwm+d/wGEikufL7QDcnGxa2SWOzEkH9XAsT40QkYeBfkBrEVkJ/EFV769fqaIcCZwHfOjnUgD+n6o+U38iRdkTeNCvLmoCPKaqDW7ZcgOlDfAfVz+RCTykqs/Vr0iV+A3wT99Q/AwYWs/yRPGK+hjgkvqWJYiqviki/wLeBbYD71EPR+fYcmrDMAwjrdhQm2EYhpFWTPEYhmEYacUUj2EYhpFWTPEYhmEYacUUj2EYhpFWTPEYVRCRHf5U3Y9EZEZkH0wS+zeIyJXV2Pm5iHQP3N8oIgPrSORayZRKRKRfotOm/fEu3f3/X/iTlWd5Nwk38/k0/H0dyRdNHxG5QETaxZMvifspInJ6XciSwP+HReQDERkT0v6vReTCVMlj1C2meIx4/KCqB/pTsL8FflUHfv4ciFZmqvp7VX2xDvzd5VDV4YGTu4cBF6nqT3H7tJLtIr8a+HsKRLoAd3xKPPnSjoi0BQ5R1f1V9a8hnT2A29dj7AKY4jGq43X8Yaki0kVEnvMHRs4VkX1iLYvIRSLytrjv3/zb75A+Andm1S2+J9Ul0mIW9w2kxwPuoz0FEfmZiLwuIu+KyOP+LLjY8F4RkdsDPbTgkS7d/fPPRGRkwM0TPg4LIwdf+oNEp3g/Poy0tEPGua9UfHvlvchuf6CZVHwv5p9+p3hE5l6+93IUcL9Pg0uBMd6fPjFh7A1sjXyewMt6l4i84ePXT9w3mxaLyJSAu42B/6cHn0XMcJ9k+KcPNzciX8S9iPzVp9VLIrJ7nPj3FJHZPo2eF/eZDERkpLhvNH0gIo/EcZcjIpN9er8nIj/1j14A2semg4h0F5FVIvKlf/aNiPxPRApUdTOwIub9Gw0VVbXLrkoXsNH/ZuAOOBzk718Cuvr/vXFHvADcAFzp/+8W8Ocm4Df+/xTg9MCzKcDpuB3x/8N9MwfgLuBcoDUwJ2B+DfD7OLK+Atzr/x8NfBSQaR7Q1Pu1Fsjyz1r531zcETC7AT1xJ1dH/G2RLM4xMswAjvT/m/k49QO+w5351wSnwI8KyNwrzv9oOsYJYygwISb9HgEEd+jjBmA/H9Z84MDgu/T/TwemxHlnURniyKTAOf7/74E7Yt5flk/n3b35mbhTP8Adt980mJ4xcboiYHcfnw9ygOLIe4zjJij3HcAFgWe/Ba6o7/JjV/WXHZljxCNX3NE37XHnOM30vY0jgMd9wx1cpR5LDxG5CWiBq4SfTxaQuqOKngNOEneUxwm4IaW+uKG513x42bjKOx4Pe7/miEhzqZiTelpVtwJbReRr3BEwK4GRInKKt9MR6Ap8AnQWkb8BT+M+BRA2zq8B/yci/wSmq+pKb/8tVV0J4NOzGPcRrtqwJ+4zAEFmqKqKyIfAalX90Ie10Ie1oJZhBSkHHvX//4E7VDJIN6AHLo+Aa6yU+mcf4HpSTwBPxPH7KOBvAKr6sYh8DuyNU6K14WucAjMaOKZ4jHj8oKoHikgeTnH8CtfCXa/u8wLJmIL7+uj7InIBruVfHY8Av8bNJ72jqt/7YamZqnpWCPex5z5F7rcGzHYAmSLSD3dC7+GqullEXgFyVHWdiBwAHIsb8joDGE2IOKvqOBF5GjgepyiPTRR+iLgk4gegMMYs4n95TFjlgbCCaVMXnziOTWsBFqpqvM9On4DrhZ4E/FZE9lPV7XUgQyJycOlkNHBsjsdIiLpx85G4IZHNwHIR+QW4E6l9RR1LAVAq7jMJ5wTMv/fP4jEbd6T9RTglBPAGcKSI7OXDy/fzHPE409s5Cvdhq++SRKsQWOeVzj7AYd5ta6CJqv4b+B1wsLpvClUbZxHpoqofqupfcKec17bVnSyNFgN71cLP1SKyr4g0AU5JYCdZuE1wQ2oAZ1O1x/YJsLuIHA7u8xgi8hMfXkdVnYUbJi3E9YCDzMXnEf9uO3n/wrId1xOOsDcN5PRsIzmmeIykqOp7uCGTs3CVxDAReR9YSPxPg1+P+/roa1T+/MAjwFV+ErlLTBg7gP/ivl/yX2/2DW611cMi8gFumC1Rhb5FRN4D7satEkvGc7iez2JgHE7BgRtWfMUPif0DuM6bh4nzaHGLEj4AyoBnq5EhETOAU+ItLsDNdx0kgTG/kFyLS9N5VAyBxTIFuDuyuCDm2Sbcx/M+AvoDNwYfqvtM/OnAX3waLcANT2YA//DDgO8BE7Xqp6n/DjTxdh7FzddsJTxvADeJyE/8/ZHAzBq4N+oJO53a2KXxQ2VXquo79S1LqhGR23HzOmlbhi4iG1W1ymrChoaIHARcrqrn1bcsRvVYj8cwdh1uBvLqW4gGSmtcb9vYBbAej2EYhpFWrMdjGIZhpBVTPIZhGEZaMcVjGIZhpBVTPIZhGEZaMcVjGIZhpJX/DxTTg0q+zpqpAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"psi_0 = p45 # Defining the initial state (+45 vector)\n",
"\n",
"phaseshift2_init = pi / 4 # initial value\n",
"phaseshift2_max = 8 * pi\n",
"n = 100 # resolution of φ (amount of steps)\n",
"step = (phaseshift2_max - phaseshift2_init)/n # interval divided by number of small steps we want\n",
"N_init = 1000 # number of iterations (range(N) -> 0 to N-1, both included)\n",
"\n",
"x_coords = [] # create x- and y- coords. arrays (x = phase shift of 2nd PA_HV, y1 = amount of photons in +45, y2 = amount of photons in -45)\n",
"y1_coords = []\n",
"y2_coords = []\n",
"\n",
"for i in range(n + 1):\n",
" output_p45 = 0\n",
" output_n45 = 0\n",
"\n",
" phaseshift2 = phaseshift2_init + i * step\n",
" x_coords.append((phaseshift2 - phaseshift1) / pi) # add realtive phase shift to x coords\n",
" PA_HV2 = Qobj([[exp(1j * phaseshift2), 0], [0, 1]]) # create corresponding PA_HV2\n",
" EffM = PA_HV2 * halfwave * PA_HV1vb # Defining the effective matrix\n",
"\n",
" psi_final = EffM * psi_0 # Applying the effective matrix to the initial state to get the final state\n",
" psi_p45 = PA_45(psi_final)[0] # Determining the probabilities\n",
" psi_n45 = PA_45(psi_final)[1]\n",
" proba_p45 = round(psi_p45.norm()**2,5)\n",
" proba_n45 = round(psi_n45.norm()**2,5)\n",
"\n",
" if (proba_p45 + proba_n45) == 1:\n",
" N = N_init # all of the photons get to the end\n",
" else:\n",
" N = int(N_init / 2) # half of the photons are blocked (should only get N/2 in the ouput) -> total proba should be 0.5\n",
"\n",
" for j in range(N):\n",
" a = random.randint(1, 50) # generates random number between 1 and 50 (both included), 50 because 50% of the photons need to come out in either +45 or -45 state (since other 50% was blocked)\n",
" if a <= proba_p45 * 100:\n",
" output_p45 = output_p45 + 1\n",
" else:\n",
" output_n45 = output_n45 + 1\n",
"\n",
" y1_coords.append(output_p45)\n",
" y2_coords.append(output_n45)\n",
"\n",
"plt.plot(x_coords, y1_coords, 'b.', markersize=9, label='Photons in state +45')\n",
"plt.plot(x_coords, y2_coords, 'r.', markersize=9, label='Photons in state -45')\n",
"legend = plt.legend(loc='upper center', fontsize='x-large')\n",
"plt.ylim([0, 1000])\n",
"plt.xlabel('Relative phase shift (multiples of π)')\n",
"plt.ylabel('Amount of photons detected')\n",
"plt.title('Amount of photons existing in |+/-45> states (V port blocked)');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Output port 'H' of the first PA$_{HV}$ is blocked"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# Polarization analyzer (HV) n°1, with H output port BLOCKED\n",
"\n",
"PA_HV1hb = Qobj([[0, 0], [0, 1]])"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAEWCAYAAAC5XZqEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABLXUlEQVR4nO2deZhUxfWw38PsDMMAgqAsM4MCSnAFxQ0hoBFX4hLFXQQ3kiC4RP0So1ETCMovSjQGNwgkKpoYlbiiIouoERQXRHEBI9qAICObwMCc74+qbrp7unvuDNM9i+d9nvt037p1q07VrapTdapuXVFVDMMwDCOTNKtvAQzDMIwfHqZ8DMMwjIxjyscwDMPIOKZ8DMMwjIxjyscwDMPIOKZ8DMMwjIxjyqceEMdkEVknIv+txf2visiIdMjW0BCR50Tkwlrc109EPm5IMvl7bxaRm+tYJCMFInKciDxZ33LUBBGZIiK3pbiuIrL3LsaxXESO2ZUwfDgRWUVkfxGZH+S+Bql8fOO6TkTy6luWRIjIRSIybxeCOAo4FuikqofWkViBqKsClylU9XhV/Vt1/uIro6rOVdUe9SlTXSMie4rIigD+Xvb5kR3ltlxEvheRjf54MY1yDggiZ9w9u9yYpuD3wLhUcflOwd/TFH8MtcmfxoKqvgeUi8jJ1fltcMpHREqBfoACp9SvNGmjBFiuqpvqWxCj4RBgVHQC8Hw1YZwL5CS5fLKqtvDHT6LuaSMiye5p1IjIIUCxqr5R37IARHcImjD/AC6rzlODUz7ABcAbwBQgxrThh3d/8WaPjSLymoh0EJE7/UjpIxE5KMr/vn4UVS4ii0XklKhrMaar+NGM7x1dLiKf+Pvv8eayfYG/Aod7GcoTJcL3Up8WkW9F5FMRucS7DwceiLr/dwnuvcin7W4R+c6na1CctxLvZ4OIvCgibaPuP8Wnt9ync1/vPg3oAszwcf8qlX9/bbmIXCMi73lZpotIvr/WVkT+4+/7VkTmikjCMiUi+4jITO/vYxE507vv5d0Ojsq3b0RkQPxzEpG9RWS2l2ONiEz37nN8NO/6dJ0V37tMlQ5//VciEhKRr0VkRKLecZTfaJkuEpF5InKHL4PLROT4RPfVAScAzya7KCLFwE3Ar2oY7rHAChGZICK9gt4kIieIyIe+DH7l87cQeA7YU3aOsvYUkUNF5HVfVkK+bOf6cKo8P+9+kogs8vfMF5H9o+K+zse5wZen+PoR5nhgdg3zIz6dpb48XOrLR0hErom6nieuDfraH3eKt9qEy6GXdyXwSKL8SRJ1W19nNvhyX5JEvmIRmerrzRci8pvoeigil4jIEh/Oh+G6FhfGvr7snu3PU+X9QSLytg9vOpAfF9yrwCCpznKlqg3qAD4FRgK9gQqgfdS1KcAafy0feAVYhlNYWcBtwCzvN8eH9f+AXGAgsAHo4a+/CoyICvsiYF7UuQL/AVrhGuxvgMGJ/CZJxxzgL17OA/39A4Pc769vB8b4dJwFfAe0iZL9M6A7UODPx/lr3YFNuAYlB9cQfQrk+uvLgWOi4gri/7/AnkAbYAlwub82FqeIc/zRD5AE6SkEvgSGAdnAQf459vTXLwE+BJoDLwB3RN0beU64ivtrXKcpHzgq7nntHXU+AFgRdZ4qHYOBlcCPvAx/jw8vLj3RMl2EK6eX4MrgFcDXifLB+78ZuLkW13J8nhWlKDf3+DJT6uXPjkv/Klw5fBE4IO7eXsDtXva3cHWwdTVlPAT08/9bAwcnynvv1hs4zD//Up//o1M8v4OA1UBfn68X+jTkAT18edrT+y0F9koi4+PAtXFuVZ6tz/u/JwkjnJ+P4Mryfj4fj/HXb8F1mHcH2gHzgVuj8mI78Ecve0Gi/EkQ5xRce3W0v+8uqrZPe/v/U4GngCIv61JguL/2M+Ar4BBAgL2Bkui2ADgY+B9wUoC8zwW+YGfbdAau/N8WJ/96YP9UaWxQIx8ROQpnknpMVRfiGthz4rz9W1UXquoW4N/AFlWdqqo7gOm4jANX0FvgGuVtqvoKTpmcXQORxqlquar+D5iFUyJB0tEZOBK4TlW3qOoi3GjnghrEvRq4U1UrVHU68DFwYtT1yaq6VFW/Bx6Lku0s4BlVnamqFcAduAJ/RJJ4gvifqKpfq+q3wIyouCqAPXCFuULdPEuizQJPwpkZJ6vqdlV9B/gXrmKgqvfjFN6bPrxfJ5G1Alc+9vT5WtN5t2TpOBOXn4tVdTOuIaoJX6jq/b4M/s2noX0Nw6iOo4F3VXVDoosi0gdX5v6c5P5zcQ1TCa4svyAircIXVfUDVb0W6IxL/wBgmYg8KiItk4RZAfQUkZaquk5V304mvK+zb/jnvxyYBPRP5h+4FJikqm+q6g51c2xbcfV6B64h7CkiOaq6XFU/SxJOK1wjHs/bvldfLs56cX0KWcL8TlU3qer7wGR2tiXnAreo6mpV/Qb4HXB+1H2VwE2qutXX16A8o6pzVHUrrk4c7tuWCCKSBQwFblDVDT5vJ0TFPwIYr6pvqeNTVf0iKoh+wNPABar6H++WKu8PwymdcNv0T1xnJZ4NuLxPSoNSPjgN+6KqrvHnDxNnesP13sJ8n+C8hf+/J/ClqlZGXf8C6FgDeVZG/d8cFXZ17Al8G9dQ1DTur+Ia8i98uNXJtqf3C4BP/5cp4g7iP1lct+OUxosi8rmIJKvAJUDfuMp+LtAhys/9uN73n31lS8SvcL23/4ozE16cxF8yUuXZl1HXov/XKFyvvCBgWZGdZstyXAN4fVQ+/SfKa8TkJiLnRplsnvMmlr8AV6rq9kTxqOprqvq9qm5W1bFAOa7hife3A3gfeBf4FvdMks0Hne7l+sKbhQ5Pkc7uPq0rRWQ98AegbTL/uDJzdVyZ6YzreHwKjMYpydVeQSYzXa3DjQjiOVhVW4UPohYkpCC6XETXx5g6RNW6+o3vLNeUSHyquhH3POLT2Rb3fOLjD9ffzrhOfDIuB+ar6qtRbknz3h+J2qZ4inBlLCkNRvmISAGuB9rfF9CVuKHdASJyQC2C/BroLLFzEF1wQ1BwpqbmUdeiG8LqqG4r8K+BNiISXeij4w5CRxGRuPu/DnDf17jCA7hl3biCE447Xvbq/CfF97SuVtWuuMUhVyWxvX8JzI6u7Oomva/wcbYA7gQeBG4WkTZJ4lupqpeo6p64Cc2/SN2skAoBnaLOOyfzWNeo6klxDeC4qDw6KcprRPmo6j9058KB44GWQB9guq834Z7oChGpomDCUeMUOeCegbj5q1eAt3GN11mq2ktV1yaR/S1VHYIzNz2JG4GHw47nXuAjoJuqtsSZwyWBvzBfAr+PKzPNVfURH/fDqhq2lCjOrJWI93Cm5bogulxE18eYOkTVuhqfH0E/JRCJz9eRNlRtA9aw0yIQHX+4/n4J7JUijsuBLiLypyi3VHkfInHbFEFEOuLMcylfdWgwygf4KW443RNnDjkQ2BeYS83MVWHexPVufyUiOeImsE8GHvXXFwGniUhz34ANr0HYq4BO4idM41HVL3F237Eiku8n64bj5hKCsjswysv+M1xeJJ1sjuIx4EQRGSRuBdPVuCFzeO39KqBrDfwnxU9K7u0L4ne451eZwOt/gO4icr5PT46IHCI7FzbcBSxQ1RHAM7h5pETx/UxEwkpiHa4Sh+OLT1dNeAwY5iddmwM31jKctCAiZUCeqi5J4uU7XI/0QH+c4N17A2+KSBcROVJEcn15vBbXY37Nhz8Y16idhTOHdVTVkaqayJwSlinXj8CK1Zlr1xP7LHYTtwAiTJH3s1FE9sHNjUUT//zuBy4Xkb7iKBSRE0WkSER6iMhAP6G9BWfxSFTuwNWZVOa9mnCjby9+hJu/nO7dHwF+IyLtxC38+S2p63qi/EnECSJylG9nbgXe8G1LBD9SfQz4vc+bEuCqqPgfAK4Rkd4+H/eW2IULG3BznkeLSHj0lzTvgddxc1jhtuk0IP51kf7AKyksGEDDUj4X4uzu//M93JWquhK4GzhXarhEUVW34ZTN8bjewV9wds2PvJc/AdtwBeFvuOWBQXkFWAysFJE1SfycjbOxf42bm7pJVV+qQRxvAt287L8HzkjWA41GVT8GzsPZ/tfg8uBknx/gFgn8xg+nrwngPxXdgJeAjbhC+RdVnZVApg3AT3C26a9xZqo/AnkiMgRX+MON0VXAweKWDMdzCK4x3YizU1+pqp/7azcDf/PpOjOA7NHyPQdMxM2FfIqbPAanhBsCJ5Ki4+Ft+dF15ht/aZV/jkW4kcc6XI94MHB8VHn6GNhH3ftL06trNKI4H1juzWiX40yp+Dr2CPC5fx57Atfg5m834Bq36XFh3UzU81PVBbhFHHd7uT/FLe4AN98zDldeV+I6ajckyZu3ge9EpG/ANKVitpfjZdyimPC7UrcBC3CjrPdxI8ekL4gmyZ9EPIxbvfgtriNxXhJ/v8RZcj4H5vn7HvJxPY5rPx7G5f2TuBFUtDzluAVHx4vIrany3pen0/z5t7gOyxNx8pxLkg5kNKIJ54eN+kRELsKtpjqqvmX5IeJHZB/gRhsJ51B2IeybAVT15hrc8yxwt6oGGfkacYjIT4CRqvrTWt5filtVm1PX5aGp4a08k1Q16fxfmIY08jGMekNEThX3vkZr3KhsRgNqaF7FjcqMWqCqL9ZW8Rg1Q1XfC6J4II3KR0QeEpHVIvJBlFsbcS9NfeJ/W3t3EZGJ4l7GfE+iXoISkQu9/0+klvtpGUYALsMtb/8MN3cVPydRV7zqj8Co6nit2RJdw2jwpM3sJiJH4+YCpqpqL+82HrcEeZy4ZbmtVfU6ETkBZ7c8Afdi012q2lfcqqcFuJU8CiwEeqvqurQIbRiGYWSEtI18VHUObkIqmiG4yX3870+j3Kf6idM3gFYisgdwHDBTVb/1CmcmbrLUMAzDaMRkepO79qoa8v9XsvMt8I7EvsC1wrslc6+CiFyKezOXwsLC3vvss08dim0YhtH0Wbhw4RpVbZeJuOpth1VVVRGpM5ufqt4H3AfQp08fXbBgQV0FbRiG8YNARBLtVpAWMr3abZU3p+F/V3v3r4h9e7iTd0vmbhiGYTRiMq18nmbnXm0X4nZiDbtf4Fe9HQZ8581zLwA/EZHWfmXcT7ybYRiG0YhJm9lNRB7B7YzbVtx3VW7CvZX8mLhv2nyB28sN3NvbJ+DepN2M27oCVf1WRG5l515Vt6jbkdgwDMNoxDTJHQ5szscwDKPmiMhCVe2TibhshwPDMAwj45jyMQzDMDJOvS21NtJPZWUla9asoby8nB07dtS3OIbRIMnPz6dTp07k5CT7Zp6RDkz5NGFWrFiBiFBaWkpOTg4iqb7dZRg/PFSVtWvXsmLFCsrKyupbnB8UZnZrwmzatImOHTuSm5triscwEiAi7LbbbmzZUpuvXBu7gimfJk6zZvaIDSMV1jGrH6xlMgzDMDKOKR/DMAwj45jyMRoly5cvR0SYN29efYtSp0yZMoXsbFsHZDR9TPkYDZKLLroIEUFEyM7OpqSkhMsvv5y1a9fWaTzz5s1DRFi+fHmdhltbzjrrLL76Kv1752ZnZzNlypQa3zdixAgGDBhQ5/LUlJUrV9KhQwdEhBUrVkTcw52S+OM3v/lNPUprJMKUjxGIlSvh9dfdb6bo168foVCI5cuXM3HiRP71r39xwQUXZE6AeqCgoID27dtX77EJMmXKlECKrbKyknPPPZdDDz00qZ+nnnqKUCgUOa6//vo6lNSoC0z5GNUybRp07QqDB7vfadMyE29ubi4dOnSgU6dODBkyhNGjR/P888/z/fffR/x8/fXXnHTSSTRv3pyuXbtW6c2HQiGGDh1Kq1atKCgoYMCAAYT3/Vu+fDn9+vUDoKysDBGJNH6qyh133EHXrl3Jzc1lr7324s4774wJu7S0lN/+9rdceeWVtGnThvbt2zNmzBi2b98e8TNv3jyOPPJIioqKKCoq4oADDuCFF5JvzB5vdgufv/baaxx88ME0b96c3r1789ZbbyUNA2Dx4sUcd9xxtGrVisLCQvbdd1+m+QdXWlrKjh07GDZsWGRkALBu3TrOO+88unTpQkFBAT169GDChAmE93+8+eabefDBB5k9e3bkvnB+b9y4kSuvvJKOHTvSvHlzDjroIJ544omUMtaWW2+9ldzcXMaMGZPUT5s2bejQoUPkaNGiRVpkMXYBVW1yR+/evdVQ/fDDD3c5jFBItaBAFXYeBQXOPZ1ceOGFOmjQoBi3CRMmKKDr16/XZcuWKaBlZWU6ffp0/eSTT/SGG27QrKws/fjjj1VVtbKyUg899FA94IADdO7cufree+/pmWeeqa1atdJvvvlGt2/frk899ZQC+t///ldDoZCuXbtWVVXvvvtuzc/P10mTJunSpUv13nvv1by8PH3ggQci8pSUlGirVq107NixunTpUp0+fbpmZ2dH/FRUVGjr1q11zJgxunTpUl26dKk+8cQTOmfOnKTpnjx5smZlZcWci4j269dP58yZo0uWLNHBgwdraWmpVlRUJA1nv/3207PPPlsXL16sn332mT777LM6Y8YMVVVdvXq1ZmVl6Z133qmhUEhD/mGGQiEdO3asLly4UD///HOdNm2aFhYW6kMPPaSqqhs2bNBzzjlHDz/88Mh9mzdv1srKSh0wYID2799f586dq5999plOmjRJc3Jy9KWXXgr8zCdPnqz9+/dP6eeVV17Rjh076qpVq3TWrFkK6Jdffhm5Hi4XnTt31t1220179+6tEyZM0G3btqUMty7qSlMAWKAZaqfrXVGk4zDl46iLCjV/vmrLlrHKp7jYuaeTeOWzePFi7dq1q/bt21dVdzYyEyZMiPjZvn27tmjRQv/617+qqupLL72kgC5evDjiZ8uWLdqhQwf93e9+p6qqc+fOVUCXLVsWE3+nTp302muvjXEbPXq0lpWVRc5LSkr05JNPjvEzePBgHTp0qKqqfvvttwrorFmzAqc7kfIBdOHChRG3N954QwH96KOPkobTsmVLnTx5ctLrWVlZKa+HGTVqlB5zzDGR8+HDh1dRELNmzdK8vDwtLy+PcR82bJgOGTKk2jjCVKd8Vq5cqXvuuafOnDkzEm+88vnmm290/Pjx+tprr+k777yjd911l7Zs2VLPO++8lHGb8nFkUvnYshojJWVlUFER67Ztm3NPN6+++iotWrRgx44dbN26lUGDBjFp0qQYPwceeGDkf1ZWFrvvvjurVq0CnOlpt912o2fPnhE/eXl59O3bl8WLFyeNd/369axYsYKjjz46xr1///7cddddbN68mebNm1eJH2DPPfdk2bJlALRu3ZoRI0Zw3HHHMXDgQPr378+pp55Kjx49apQPIsIBBxwQEwfAqlWrkoZ1zTXXMGLEiMg8yimnnMLBBx+cMp7KykrGjx/Po48+yooVK9iyZQsVFRWUlJSkvO+tt95i27ZtdOzYMcZ927ZtdOvWLel9c+fO5fjjj4+cb9++nYqKihgTWb9+/XjuuecAOPfcc7ngggs45phjkobZtm1brr322sj5gQceSFFRERdffDHjxo2rIqNRf9icj5GSDh1g0iQoKIDiYvc7aZJzTzd9+/Zl0aJFLFmyhC1btjBz5ky6du0a4yc3NzfmXESorKxMv3AB47///vtZuHAhxx57LLNnz6ZXr15VFGh1NGvWjKysrJg4gJTpvPHGG1m6dClnnnkmH3zwAYcddli1K74mTJjA2LFjGTVqFDNnzmTRokWMGDGCbdu2pbyvsrKS4uJiFi1aFHN8+OGHEcWRiD59+sT4v+WWW6q4PfDAAxH/L7/8MrfffjvZ2dlkZ2czaNAgwM1hXXbZZUnjOeKIIwD44osvUqbDyCw28jGq5fzz4dhjYdkyN+LJhOIBt/Jr7733rvX9P/rRj1i7di0ffvhhZPSzdetW3nzzTUaOHAnsVB7Ru363bNmSTp06MWfOHE466aSI++zZsykrK4uMeoLSq1cvevXqxVVXXcXll1/Offfdl7KxrCu6du3KyJEjGTlyJOPGjeP222/ntttuA1y643c6nzNnDoMHD+biiy+OuH3yyScxfhLd16dPH8rLy9myZQu9evUKLF/88919991TPvP3338/5vytt97i4osv5oUXXmDfffdNGs/bb78NQKdOnQLLZqQfUz5GIDp0yJzSqSsGDhzIoYceyjnnnMM999xDcXExt956K1u2bOGKK64AoKSkhGbNmvHss89y1llnkZeXR3FxMTfccANXX3013bp1Y8CAAbzyyivce++93HPPPYHj//TTT7n//vs5+eST6dy5M19//TVz586t1vy1q2zcuJHrrruO008/nbKyMsrLy3n++edjzI9lZWXMmjWL448/ntzcXNq2bUuPHj2YNm0as2bNomPHjkydOpU333yT1q1bx9z3+OOPs3jxYtq3b09RUREDBw7kmGOO4bTTTmP8+PHsv//+rFu3jvnz55Ofn88ll1xSJ+mKV2xr1qwBoEePHhFT5JQpU8jKyuLggw8mPz+fuXPncu2113LGGWfQpUuXOpHDqBvM7GY0WUSEJ598kn322YcTTzyRQw45hJUrVzJz5kzatm0LQPv27Rk7dizjxo1jjz32YMiQIQBcccUV3HLLLfzhD3+gZ8+e/PGPf2TcuHEMHz48cPyFhYV88sknDB06lO7du3P66adzxBFHcPfdd6clvWGys7NZt24dw4cPZ9999+W4446jffv2PPzwwxE/EyZMYOHChZSWltKuXTvAmer69+/PkCFDOPzww1m3bh2jRo2KCXv48OEccsghHHHEEbRr145HHnkEEeHpp5/mtNNOY8yYMZH8fuaZZ9hrr73SmtZ4mjVrxvjx4+nbty/7778/d9xxB9deey1///vfMyqHUT3iFjg0Lfr06aPhdzl+yCxZsiSlOcIwDIfVFYeILFTVPpmIy0Y+hmEYRsYx5WMYhmFkHFM+hmEYRsYx5WMYhmFkHFM+hmEYRsYx5WMYhmFkHFM+hmEYRsYx5WMYhmFkHFM+hmEYRsYx5WM0SpYvX46IMG/evPoWpU6J/5KpYTRVTPkYDZKLLroo8qnm7OxsSkpKuPzyy1m7dm2dxjNv3jxEhOXLl9dpuLXlrLPO4quvvkp7PNnZ2VU+OR6EESNGRD41nm4qKysZNGgQIlJlb7bS0tJI+QgfRx11VEbkMuoG62IZwVi5MuPfVOjXrx+PPfYY27dvZ+HChYwYMYIvv/ySZ555JiPx1wcFBQUUFBTUtxgNgltuuYXCwsKk16+77jpGjx4dOY//tpLRsLGRj1E906ZB164weLD7nTYtI9Hm5ubSoUMHOnXqxJAhQxg9ejTPP/8833//fcTP119/zUknnUTz5s3p2rVrld58KBRi6NChtGrVioKCAgYMGEB409nly5fTr18/wH0qQEQivXpV5Y477qBr167k5uay1157ceedd8aEXVpaym9/+1uuvPJK2rRpQ/v27RkzZgzbt2+P+Jk3bx5HHnkkRUVFFBUVccABB/DCCy8kTXO82S18/tprr3HwwQfTvHlzevfuzVtvvZUy7xYvXsxxxx1Hq1atKCwsZN9992Waf26lpaXs2LGDYcOGRUYNAOvWreO8886jS5cuFBQU0KNHDyZMmEB48+Gbb76ZBx98kNmzZ0fuC+f3xo0bufLKK+nYsSPNmzfnoIMO4oknnkgpYypeeeUVHnroISZPnpzUT4sWLejQoUPkaNOmTa3jM+qBTH2vO5NH7969a/jl8qZJnXyXPhRSLShQhZ1HQYFzTyMXXnihDho0KMZtwoQJCuj69et12bJlCmhZWZlOnz5dP/nkE73hhhs0KytLP/74Y1VVrays1EMPPVQPOOAAnTt3rr733nt65plnaqtWrfSbb77R7du361NPPaWA/ve//9VQKKRr165VVdW7775b8/PzddKkSbp06VK99957NS8vTx944IGIPCUlJdqqVSsdO3asLl26VKdPn67Z2dkRPxUVFdq6dWsdM2aMLl26VJcuXapPPPGEzpkzJ2m6J0+erFlZWTHnIqL9+vXTOXPm6JIlS3Tw4MFaWlqqFRUVScPZb7/99Oyzz9bFixfrZ599ps8++6zOmDFDVVVXr16tWVlZeuedd2ooFNKQf5ahUEjHjh2rCxcu1M8//1ynTZumhYWF+tBDD6mq6oYNG/Scc87Rww8/PHLf5s2btbKyUgcMGKD9+/fXuXPn6meffaaTJk3SnJwcfemllwI/8zArV67Ujh076qxZs1RVFdBp06bF+CkpKdH27dtrmzZttGfPnvrLX/5S16xZU+O4wtRJXWkCAAs0Q+108gswA3g62ZEpAWtzmPJx1EmFmj9ftWXLWOVTXOzc00i88lm8eLF27dpV+/btq6oaUT4TJkyI+Nm+fbu2aNFC//rXv6qq6ksvvaSALl68OOJny5Yt2qFDB/3d736nqqpz585VQJctWxYTf6dOnfTaa6+NcRs9erSWlZVFzktKSvTkk0+O8TN48GAdOnSoqqp+++23CkQa0SAkUj6ALly4MOL2xhtvKKAfffRR0nBatmypkydPTno9Kysr5fUwo0aN0mOOOSZyPnz4cO3fv3+Mn1mzZmleXp6Wl5fHuA8bNkyHDBlSbRzR7NixQwcNGqQ33nhjxC2R8rnjjjt05syZ+t577+ljjz2m3bt31+7du+vmzZtrFF8YUz6OTCqfVHM+d/jf04AOQHjG72xg1a6MtkRkDDACUOB9YBiwB/AosBuwEDhfVbeJSB4wFegNrAXOUtXluxK/UQPKyqCiItZt2zbnnmZeffVVWrRowY4dO9i6dSuDBg1i0qRJMX4OPPDAyP+srCx23313Vq1yxXPx4sXstttuMV/wzMvLo2/fvixevDhpvOvXr2fFihUcffTRMe79+/fnrrvuYvPmzZFPaUfHD7DnnnuybNkyAFq3bs2IESM47rjjGDhwIP379+fUU0+lR48eNcoHEeGAAw6IiQNg1apVScO65pprGDFiBFOmTGHAgAGccsop1X5BtbKykvHjx/Poo4+yYsUKtmzZQkVFBSUlJSnve+utt9i2bRsdO3aMcd+2bRvdunVLel+LFi0i//v168dzzz3HH/7wB7Zu3cpNN92UMs6rr7468n+//fajd+/e7L333vz73//mnHPOSXmv0TBIOuejqrNVdTZwpKqepaoz/HEO0K+2EYpIR2AU0EdVewFZwFDgj8CfVHVvYB0Q/mTkcGCdd/+T92dkig4dYNIkKCiA4mL3O2lSRhYd9O3bl0WLFrFkyRK2bNnCzJkz6dq1a4yf+ElmEaGysjLtsgWN//7772fhwoUce+yxzJ49m169elVRoNXRrFkzsrKyYuIAUqbzxhtvZOnSpZx55pl88MEHHHbYYfzmN79JGc+ECRMYO3Yso0aNYubMmSxatIgRI0awbdu2lPdVVlZSXFzMokWLYo4PP/yQ5557Lul90X4feOABAF566SXmz59PXl4e2dnZkfmvCy+8kH322SdpWF27dqV9+/YNZtWiUT1BVrsVikhXVf0cQETKgORLUILHWyAiFUBzIAQMBMJdlr8BNwP3AkP8f4B/AneLiPghopEJzj8fjj0246vdCgoK2HvvvWt9/49+9CPWrl3Lhx9+GBn9bN26lTfffJORI0cCO5XHjh07Ive1bNmSTp06MWfOHE466aSI++zZsykrK4uMeoLSq1cvevXqxVVXXcXll1/Offfdx2WXXVbrdAWla9eujBw5kpEjRzJu3Dhuv/12brvtNsClOzrNAHPmzGHw4MFcfPHFEbdPPvkkxk+i+/r06UN5eTlbtmyhV69egeVL9GwnT57Mpk2bYtz2228/fv/733P66acnDeurr75i9erVdO7cOXD8Rv0SRPmMAV4Vkc8BAUqAWtccVf1KRO4A/gd8D7yIM7OVq2p4mdAKIDyG7wh86e/dLiLf4Uxza6LDFZFLgUsBunTpUlvxjGR06JAxpVNXDBw4kEMPPZRzzjmHe+65h+LiYm699Va2bNnCFVdcAUBJSQnNmjXj2Wef5ayzziIvL4/i4mJuuOEGrr76arp168aAAQN45ZVXuPfee7nnnnsCx//pp59y//33c/LJJ9O5c2e+/vpr5s6dW635a1fZuHEj1113HaeffjplZWWUl5fz/PPPx5gfy8rKmDVrFscffzy5ubm0bduWHj16MG3aNGbNmkXHjh2ZOnUqb775Jq1bt4657/HHH2fx4sW0b9+eoqIiBg4cyDHHHMNpp53G+PHj2X///Vm3bh3z588nPz+fSy65JLDsZUnMuZ06dYqY8F5//XVee+01Bg4cyG677cZHH33E9ddfT5cuXTj11FNrmWtGpql2qbWqPg90A67Emct6qGrytaLVICKtcaOZMmBP3ChqcG3Di5LzPlXto6p92rVrt6vBGU0AEeHJJ59kn3324cQTT+SQQw5h5cqVzJw5k7Zt2wLQvn17xo4dy7hx49hjjz0YMmQIAFdccQW33HILf/jDH+jZsyd//OMfGTduHMOHD08VZQyFhYV88sknDB06lO7du3P66adzxBFHcPfdd6clvWGys7NZt24dw4cPZ9999+W4446jffv2PPzwwxE/EyZMYOHChZSWlhKuLzfeeCP9+/dnyJAhHH744axbt45Ro0bFhD18+HAOOeQQjjjiCNq1a8cjjzyCiPD0009z2mmnMWbMmEh+P/PMM+y11151nr68vDyeeOIJjjnmGLp3787IkSM57LDDeP3112PmkYyGjVRnvRKR5sBVQImqXiIi3XAK6D+1ilDkZ8BgVR3uzy8ADgd+BnTwo5vDgZtV9TgRecH/f11EsoGVQLtUZrc+ffpo+F2OHzJLlixh3333rW8xDKPBY3XFISILVbVPJuIK8pLpZGAbTkEAfAXctgtx/g84TESai5s5HQR8CMwCzvB+LgSe8v+f9uf466/YfI9hGEbjJojy2UtVxwMVAKq6GTf3UytU9U3cwoG3ccusmwH3AdcBV4nIp7g5nQf9LQ8Cu3n3q4Draxu3YRiG0TAIsuBgm4gU4N7JQUT2ArbuSqSqehMQv5D/c+DQBH634ExyhmEYRhMhiPK5GXge6Cwi/wCOxL0UahiGYRi1olrlo6ovishC4DCcue1KVV1TzW1GA0FVIy8lGoZRFZtCrh+qnfMRkZdVda2qPqOq/1HVNSLyciaEM3aNnJycmB2gDcOoSkVFhX3Arx5IqnxEJF9E2gBtRaS1iLTxRyk7XwA1GjC77747X331FZs3b7benWEkoLKyklWrVlFcXFzfovzgSKXuLwNG414EXcjOFW7rgfS+JWfUCS1btgTcN28q4jcHNQwDcC8Dh186NjJHUuWjqncBd4nIL1X1zxmUyahDWrZsGVFChmEYDYUg7/lUikir8Ik3wY1Mn0iGYRhGUyeI8rlEVcvDJ6q6Dgi+U6BhGIZhxBFE+WRJ1FpdEckCclP4NwzDMIyUBFlf+DwwXUTCX8C6zLsZhmEYRq0IonyuwymcK/z5TOCBtElkGIZhNHmC7HBQKSJTcLtJf5x+kQzDMIymTpAdDk4BFuFNbSJyoIg8nWa5DMMwjCZMkAUHN+F2my4HUNVFuK+QGoZhGEatCKJ8KlT1uzg326vFMAzDqDVBFhwsFpFzcEuuuwGjgPnpFcswDMNoygQZ+fwS+BHuA3IPA98BV6ZTKMMwDKNpE2Tkc6Kq/hr4ddhBRH4GPJ42qQzDMIwmTZCRzw0B3QzDMAwjEElHPiJyPHAC0FFEJkZdaglsT7dghmEYRtMlldnta2ABcAruez5hNgBj0imUYRiG0bRJ9T2fd4F3ReRh76+L7XBgGIZh1AVB5nwGYzscGIZhGHVIEOVzM7bDgWEYhlGH2A4HhmEYRsaxHQ4MwzCMjFPTHQ4eAdYDo9Mok2EYhtHECfI9n8243Q1+XZ1fwzAMwwhCqpdMZ5BibkdVT0mLRIZhGEaTJ9XI5w7/exrQAfi7Pz8bWJVOoQzDMIymTaqXTGcDiMgEVe0TdWmGiCxIu2SGYRhGkyXIgoNCEekaPhGRMqAwfSIZhmEYTZ0gS63HAK+KyOeAACXApWmVyjAMw2jSBFnt9rx/v2cf7/SRqm5Nr1iGYRhGUyaI2Q1V3aqq7/pjlxWPiLQSkX+KyEciskREDheRNiIyU0Q+8b+tvV8RkYki8qmIvCciB+9q/IZhGEb9Ekj5pIG7gOdVdR/gAGAJcD3wsqp2A1725wDHA938cSlwb+bFNQzDMOqSjCsfESkGjgYeBFDVbapaDgwB/ua9/Q34qf8/BJiqjjeAViKyR0aFNgzDMOqUapWPiBwpIoX+/3ki8n8iUrILcZYB3wCTReQdEXnAh99eVUPez0qgvf/fEfgy6v4V3i1ezktFZIGILPjmm292QTzDMAwj3QQZ+dwLbBaRA4Crgc+AqbsQZzZwMHCvqh4EbGKniQ0AVVVquHO2qt6nqn1UtU+7du12QTzDMAwj3QRRPtu9MhgC3K2q9wBFuxDnCmCFqr7pz/+JU0arwuY0/7vaX/8K6Bx1fyfvZhiGYTRSgiifDSJyA3Ae8IyINANyahuhqq4EvhSRHt5pEPAh8DRwoXe7EHjK/38auMCvejsM+C7KPGcYhmE0QoK8ZHoWcA4wXFVXikgX4PZdjPeXwD9EJBf4HBiGU4SPichw4AvgTO/3WeAE4FNgs/drGIZhNGLEWdSaFn369NEFC2z7OcMwjJogIgvj9vJMG0FWu53mX/z8TkTWi8gGEVmfCeEMwzCMpkkQs9t44GRVXZJuYQzDMIwfBkEWHKwyxWMYhmHUJUFGPgtEZDrwJBDZ101Vn0iXUIZhGEbTJojyaYlbZfaTKDcFTPkYhmEYtSLIJxVsabNhGIZRpwRZ7dZJRP4tIqv98S8R6ZQJ4QzDMIymSZAFB5Nxuwzs6Y8Z3s0wDMMwakUQ5dNOVSer6nZ/TAFs507DMAyj1gRRPmv9pxSy/HEesDbdghmGYRhNlyDK52LcPmsrgRBwBnBRGmUyDMMwmjhBllp3UtVToh1E5EhiP/BmGIZhGIEJMvL5c0A3wzAMwwhE0pGPiBwOHAG0E5Groi61BLLSLZhhGIbRdElldssFWng/0V8uXY+b9zEMwzCMWpFU+ajqbGC2iExR1S9EpIV335gx6QzDMIwmSZAFB0Ui8g7QBkBE1gAXquoHaZXMMAzDaLIEWXBwH3CVqpaoaglwtXczDMMwjFoRRPkUquqs8ImqvgoUpk0iwzAMo8kTxOz2uYjcCEzz5+cBn6dPJMMwDKOpE3SHg3a47/c84f9fnE6hDMMwjKZNkO/5rANGZUAWwzAM4wdCtcpHRLoD1wCl0f5VdWD6xDIMwzCaMkHmfB4H/go8AOxIrziGYRjGD4Egyme7qt6bdkkMwzCMHwyp9nZr4//OEJGRwL+BreHrqvptmmUzDMMwmiipRj4LAQXEn18bdU2BrukSyjAMw2japNrbrSyTghiGYRg/HIKsdssHRgJH4UY8c4G/quqWNMtmGIZhNFGCLDiYCmxg5wfkzsHtdvCzdAllGIZhNG2CKJ9eqtoz6nyWiHyYLoEMwzCMpk+Q7XXeFpHDwici0hdYkD6RDMMwjKZOkJFPb2C+iPzPn3cBPhaR9wFV1f3TJp1hGIbRJAmifAanXQrDMAzjB0WQjUW/yIQghmEYxg+HIHM+aUFEskTkHRH5jz8vE5E3ReRTEZkuIrnePc+ff+qvl9aXzIZhGEbdkFT5iEhemuO+ElgSdf5H4E+qujewDhju3YcD67z7n7w/wzAMoxGTauTzOoCITEvhp1aISCfgRNxO2YiIAAOBf3ovfwN+6v8P8ef464O8f8MwDKORkmrOJ1dEzgGOEJHT4i+q6hO7EO+dwK+AIn++G1Cuqtv9+Qqgo//fEfjSx7ldRL7z/tdEBygilwKXAnTp0mUXRDMMwzDSTSrlczlwLtAKODnumuI+qV1jROQkYLWqLhSRAbUJIxGqeh9wH0CfPn20rsI1DMMw6p5UG4vOA+aJyAJVfbAO4zwSOEVETgDygZbAXUArEcn2o59OwFfe/1dAZ2CFiGQDxcDaOpTHMAzDyDBBVrtNE5FRIvJPf/xSRHJqG6Gq3qCqnVS1FBgKvKKq5wKzgDO8twuBp/z/p/05/vorqmojG8MwjEZMEOXzF9wuB3/xx8FAOr5seh1wlYh8ipvTCY+2HgR28+5XAdenIW7DMAwjgwTZ4eAQVT0g6vwVEXm3LiJX1VeBV/3/z4FDE/jZgu2gbRiG0aQIMvLZISJ7hU9EpCuwI30iGYZhGE2dICOfa3GfUfgc90ntEmBYWqUyDMMwmjRB9nZ7WUS6AT2808equjW9YhmGYRhNmSAjH7yyeS/NshiGYRg/EOptY1HDMAzjh4spH8MwDCPjVKt8ROTlIG6GYRiGEZSkcz4ikg80B9qKSGvcSjdw2+F0THafYRiGYVRHqgUHlwGjgT2BhexUPuuBu9MrlmEYhtGUSbWx6F3AXSLyS1X9cwZlMgzDMJo4Qd7z+bOIHAGURvtX1alplMswDMNowlSrfPyXTPcCFrFzWx0FTPkYhmEYtSLIS6Z9gJ72GQPDMAyjrgjyns8HQId0C2IYhmH8cAgy8mkLfCgi/wUie7qp6ilpk8owDMNo0gRRPjenWwjDMAzjh0WQ1W6zMyGIYTQlVq6EZcugrAw6mNHaMKoQZHudDSKy3h9bRGSHiKzPhHCGsXIlvP66+20sTJsGXbvC4MHud9q0+pbIaGo0xnoRT7XKR1WLVLWlqrYECoDTgb+kXTLjB0u4Yv35z42vEV+5Ei67DL7/Htavd7+XXda4GwmjYdFUOjc12tVaHU8Cx6VHHKOhUlc9rerCCVesn/wERo1qfI34smWQkxPrlpvr3I1dJEAhbAojglQ0pc5NELPbaVHHGSIyDtiSAdnqhQZReBuEEDupq55WdeFEV6yNG6ven50NTz9dNVsaTHatXEn3ta/TZlusIFu3wtq1qeXLlHKvTThBw0zrc0hReBrzSLmmNKnOjaqmPIDJUcf9wK+B3au7rz6P3r17a22YOlW1oEC1ZUv3O3VqYn+hkOr8+e43iHsQwveumxhQiJpGHMBPIi+hkBMDdh4FBYmDSRVFkHDmz3fJjvYTfxQVxWZL/DObOLFmzyBa5qTyB8nfKEEqcgp0WM5ULS5WzclxR6rHGbTcVceu5kWicILIX5dpiCcUUn1rRkgr8+MKT36+6owZ+tjEkBYUqLZoUbWsJCundVVfdsF7yntSlcma1MfaACzQDLXT9a4o0nHURvkEfajJKlmqylddIZs40d2zd4uQbiJxJYu/OYiyqolCSyR/KKQ6aZJr8KNFKi524SZKQ8IoQiF9b9J83btFKGk4yZ4BqBYWVnXLz3dx5OdXr6Din0EimZM2sgkyJuZ5hkLu+cQJUplfoC9ODVWRL6ZMhUK6ZsZ8Lc0PJX3kQRu1ZHmXKC9qE06qOlFt3allQx5+PoOaz9dyqvZKdjQv1E0U6HlMTShrTPlKUlDXTZxapVwEqVPJylHCvE6Q/iDtSLIyGfZTXOx+H5sYpAcVjAalfIBOwL+B1f74F9ApUwLW5qiN8knU6y4qco1v+JkmaGO0ND+kn0yt2oCU5IV07u3z9YHbQgl7o+ECG91jO4z5ui5BJdPCwpiKkkpZVeYX6FszQoH8pEpXTo5zS9SjDDeOidIQnS9rZuwUdkdRyyoNRaJGNhxmuGJNnJhYAYbvTzVKipczukInkjle/rVTq2ZMRW6BluaHtGVL1WE5U7Uip0C1efOqARQX63uT5ictU+EGrqKwar6EH3l84xPdyKx61yn0Ve+GknYSEj6P+AYq6n+QcKLrRFiGF6eGqk1nouFYfBqS1Yv2JOiURR2bKND2hKpcCj//iAwJHvomCnTvFqFIXqfqAIZHWdWVo5i8jtMy6yZOTVjfUnWmEoWbUFEGHaamoKEpn5nAMNw7QdnARcDMTAlYm6OuRj7tCemg5vO1U3ZIO2WHdEDe/JhCfh5TdRMFujk3tgEJu5dTtWFpT0h/nD8/YWWpvpLl6+l5MyL3JlJW31Gkv8idVL2fvElJ05Wq4cnJUe2UHdKBBcnvCad/S14LrUxQ2bsVhSLhROfvoOZOicebi6rrjSc72hPS/rnB0pZI/q25zavIX06xHsb8ap+VFhToqndDCcvUGQUzdDOxrUyyBjRepm3NW2pFsxzdSk6kfI2SicGeR35LrczO0cqcHK0ojP2/LadAr8xOHU50nRglE2PK+IXNpgZKZyVoRfOiWqehIr/qMHidfyaJymkiGVLdm6i+VIJW5KUeZSXM67x8J0hcHT4td0aVdFZXVmPCzY8yTdR0mFoNDU35LAri1pCOXZ3zKSra+bDX0VK3kKNbyNF1UcokUeOziQLtxbsJ3PP1BGboSCZGwtxEgf6ciZGGLL6Q7WhetZJVgm6gMHLvCSSu3N9RFMhPonTFV4iwfEVFquPHqw7LnpoyDdU1yjuKinXu7fNjwomXY1jO1KTmzrAJLlq2cD1LVFkTyRl/b3yaq+tl9+JdHcEk/Y4EQwQ/Sg33OhOVqfUkVmo/zp9fJW3tCSV8hsme+XlMjTFTVqskU4QTbsRzc2PlLydxp6IkL1RtOoPEnaheQPLRaLhDEx5YzZixs5xWJ0O4bgbN63j/1T2n+Lij63B8ZzVVPYx/hhW5BS6hqSZJ4+3aAWhoyudl4Dwgyx/nAS9nSsDaHLVVPqp+pDy++hHIddm3V2l8yinS27N+VcW9EnQjBQkLYnSFC/f8H5vobGEVOdU3OFvI0a3k6pa8oqThp/ITn67TcmdoSV5IL2wWWyGG5UzVF6cmMOHFpeHnuUka5ahK887Uqgo6vjF5a0bVHtuqd50pc3RWrBIflj01Yl4rzQ/paXnBFW443ztlh7RzTkh/kVdV/krQivxC3ZFXoJOaXZq0Aa7IST4/V12Z2pbj5ogSKeWaNOKbKNCp40ORvEiUnpqEE7bMVSd/Oa5TUZ2/IEeyehG2Is29NKwEi3UTBTrr9IkR852qW6CQSoZK0G0FsWWh6v9c3UBhUuWR6N4NCZ5TdfVtWPbUBKNgVw9TlclyinXt1BnO7JuirjX2kU8J8DTwjZ/zeRLokikBa3PsivLRUEjLb09dYStBt+YkViZbcqsW2KAVrqJ5kVbm77QL7+y9pQ5zE/l6Q05VZVhTP2HzQmVevm7PijUXVOQW6NqpMxJO/MbkS4L0Rzcmw3Km6idTE08gR1esNTPiemxTd9rtkzb6EydqZX6BVuTVrMdd0bwoYoJKbCp0jUGfvHd1a1Zy5RtWgolWyyUqU9G94CuzJ/rGpJqJrGqO7yjS8tsnRfIiUXpqFI7XPtXViU0U6Dcvv5s0nYk6XzWtF+HFHgUFO0eHYWvCjqKd8xxrZlQtX/FWg+rMcZvI13OpqhhqemwhRzeTn7AOV/qym9A9r7CaMuk7KznhzkpxRGmu80o5kQWhOhqU8mmMR62VT4oGLlXlSFawNtagxxpfsMoKdpqxqjMFrKNYT2BGyt5eED8pD293C9o4JjOjdCsK6YfjUzeyVXpsQSZ9ajMpFFD+sBnkMKo2auUU6Qh2zrHFWd1Slqloc+xmCrQiP8HChTiZdva6c/U7Eo92Ey0PjB355cY0VuXJwgnb3fzMeqpOxUcDLg2UzviGMlEakqY/343GwgsiEpoU/cx9fPmKNpeFn2XCxT1R9eUw5keZEIN3KuNNa0HqcE3L5LCcqRGrW7z5L/y/FlY3bVDKBygD/g94wo+AngaezpSAtTlqpXwSNHCxveJcrcivWgC/o0gn5FY1te0oLNJlI8e7Scf4MAtSF+T4SVCImnAuqFpZw5PVOyclA/gpKE6ZroSHb4wqc3K1IoEc8fnyi7xJMXb7yHxWUcud4RRGyVFY7Hq48SsOgrwAlKSybs+vXY/7O4r057mx8rcnpJsTzPMlmjsqKHBmwkRlaltBUUzDFLRDsIl8/VnzGVqaH9IHbnMrxcpv86O9AM/jl/mTIqacYwrnx/wfnT1Rt+VUH070SKT8Nmfu+ubld6tNZ3Qehs1oNU1DOcU6sCD14gCFyFLBypxc3da8OPA8SrLnmmzxRHXPqSQvFLPeILJoIi9xZ3V7btVRUKJnGDZBpmG9gTY05fMuMAr4MdA/fGRKwNoctVI+tVxrXZnv5jCqvAQXfvLxi/Lj11onWNca36CFe9ORpbaR5ctVK1bVJc4p/ARZa51sjXP0GuYEfuKXfHcrquYdpmQvTYRnkBOtQU20xDmVnOG3PsMz6Kka2STv6QzLmeqedXGxVuS6Hmii95C6FYX0018lWLfsy1R42e4xhUlMkOG11rm5kbIT/U5KDAHWSMcvr0/4WkhN11oHqDvhdMYX/6Qv8tagXuxVWI3ijlsiHR/kzs5QcUxeh59rvMyRZc3Fsf5TPaeE7+TUcq119DMM51d0/saJUasXfRua8nkzU8LU1VFXI5+Ub5nGP+Fk7uGwA7xlmkhRJHnHNGHhS/RyZ7V+UqWrNm+ZJogkFFJ9b9J8N+JJFk6yZxBuyOJrViqlVNO3TFPU3ISPNe7dmHgxIo1aYerX7kOh1G/vN5W3TGv83mPc86korFovqrxHlKgHEFW+khXTZC9oJpU5meZOkciEl4K0I8m0SYJ3h+roHdMGp3zOAW4CDgcODh+ZErA2xy7P+VTXdahOmdTmyddEUQQUp8aiJfIUVClXF0mQcKozryXSxMlGlkGfQYAGJEjehcVIOMKrTgkELXfVsat5kSicoN3pukpDPKHEu0DEFJ1ko/ck5XRXG+g6I0g7Eu+nJvWxFjQ05TMWWAHMBmb545VMCVibY1dXu9V3yWwAIsSSrsYxPpzqet3JZlAbSIYlHeElMlUlurku0pCOcGoy+krTcwhUBNOlABsSiTpotVlZkIRMKh9x8SVHRD4FeqrqtpQeAyIinYGpQHtAgftU9S4RaQNMB0qB5cCZqrpORAS4CzgB2AxcpKpvp4qjT58+umDBgroQ1whTV5/mrC6cadPc1tbZ2bBhQ+y1ggL4/POG/WnQlSvdlsrff7/TrTHI3QgIVASb+idk01y+RGShqvbZ5YACEOR7Ph8Areowzu3A1araEzgM+LmI9ASux7282g33Yuv13v/xQDd/XArcW4eyGEHp0AEOP3zXC3h14Zx/vqtIL7wAEye6ilVc7H4nTWr4DUqHDk7OxiZ3IyBQEayrctpQaULlK8jI51Vgf+AtYGvYXVVPqRMBRJ4C7vbHAFUNicgewKuq2kNEJvn/j3j/H4f9JQvTRj5NiMbak22schuNgzSVr0yOfLID+LkpXZGLSClwEPAm0D5KoazEmeUAOgJfRt22wrvFKB8RuRQ3MqJLly7pEtnINB06NM7Gu7HKbTQOmkD5qtbspqqzow9gB3DmrkYsIi1wn2cYrarr4+JU3HxQYFT1PlXto6p92rVrt6viGYZhGGkkyMgHETkIt+T6Z8AynNKoNSKS48P4h6o+4Z1XicgeUWa31d79K6Bz1O2dvJthGIbRSEk68hGR7iJyk4h8BPwZ+B9ujujHqnp3bSP0q9ceBJao6v9FXXoauND/vxB4Ksr9AnEcBnyXar7HMAzDaPikGvl8BMwFTlLVTwFEZEwdxHkkcD7wvogs8m7/DxgHPCYiw4Ev2Gnaexa3zPpT3FLrYXUgg2EYhlGPpFI+pwFDgVki8jzwKCC7GqGqzksRzqAE/hX4+a7GaxiGYTQckprdVPVJVR0K7IPb1WA0sLuI3CsiP8mQfIZhGEYTJMhqt02q+rCqnoyb7H8HuC7tkhmGYRhNliA7HERQ1XV+SXMV85hhGIZhBKVGyscwDMMw6gJTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGaTTKR0QGi8jHIvKpiFxf3/IYhmEYtadRKB8RyQLuAY4HegJni0jP+pXKMAzDqC2NQvkAhwKfqurnqroNeBQYUs8yGYZhGLUku74FCEhH4Muo8xVA32gPInIpcKk/3SgiH+9CfG2BNbtwfzox2WpHQ5YNGrZ8JlvtaMiyQWL5SjIVeWNRPtWiqvcB99VFWCKyQFX71EVYdY3JVjsasmzQsOUz2WpHQ5YN6l++xmJ2+wroHHXeybsZhmEYjZDGonzeArqJSJmI5AJDgafrWSbDMAyjljQKs5uqbheRXwAvAFnAQ6q6OI1R1on5Lk2YbLWjIcsGDVs+k612NGTZoJ7lE1Wtz/gNwzCMHyCNxexmGIZhNCFM+RiGYRgZx5RPFA15Cx8ReUhEVovIB/UtSzwi0llEZonIhyKyWESurG+ZwohIvoj8V0Te9bL9rr5likdEskTkHRH5T33LEo+ILBeR90VkkYgsqG95ohGRViLyTxH5SESWiMjh9S0TgIj08PkVPtaLyOj6liuMiIzxdeEDEXlERPLrRQ6b83H4LXyWAsfiXmJ9CzhbVT+sV8E8InI0sBGYqqq96lueaERkD2APVX1bRIqAhcBPG0LeiYgAhaq6UURygHnAlar6Rj2LFkFErgL6AC1V9aT6licaEVkO9FHVBveypIj8DZirqg/4VbDNVbW8nsWKwbcrXwF9VfWLBiBPR1wd6Kmq34vIY8Czqjol07LYyGcnDXoLH1WdA3xb33IkQlVDqvq2/78BWILblaLeUcdGf5rjjwbT4xKRTsCJwAP1LUtjQkSKgaOBBwFUdVtDUzyeQcBnDUHxRJENFIhINtAc+Lo+hDDls5NEW/g0iAa0MSEipcBBwJv1LEoEb9ZaBKwGZqpqg5ENuBP4FVBZz3IkQ4EXRWSh38KqoVAGfANM9ibLB0SksL6FSsBQ4JH6FiKMqn4F3AH8DwgB36nqi/Uhiykfo84QkRbAv4DRqrq+vuUJo6o7VPVA3M4Yh4pIgzBbishJwGpVXVjfsqTgKFU9GLej/M+9+bchkA0cDNyrqgcBm4CGNk+bC5wCPF7fsoQRkdY4i04ZsCdQKCLn1Ycspnx2Ylv47AJ+PuVfwD9U9Yn6licR3iwzCxhcz6KEORI4xc+rPAoMFJG/169IsfieMqq6Gvg3zjzdEFgBrIgaxf4Tp4waEscDb6vqqvoWJIpjgGWq+o2qVgBPAEfUhyCmfHZiW/jUEj+p/yCwRFX/r77liUZE2olIK/+/ALeg5KN6FcqjqjeoaidVLcWVt1dUtV56oYkQkUK/gARv0voJ0CBWW6rqSuBLEenhnQYB9b7AJY6zaUAmN8//gMNEpLmvt4Nwc7QZp1Fsr5MJ6mELnxohIo8AA4C2IrICuElVH6xfqSIcCZwPvO/nVgD+n6o+W38iRdgD+JtfddQMeExVG9yS5gZKe+Dfro0iG3hYVZ+vX5Fi+CXwD99Z/BwYVs/yRPDK+ljgsvqWJRpVfVNE/gm8DWwH3qGettmxpdaGYRhGxjGzm2EYhpFxTPkYhmEYGceUj2EYhpFxTPkYhmEYGceUj2EYhpFxTPkYVRCRHX433g9EZEb4PZkU/m8WkWuq8fNTEekZdX6LiBxTRyLXSqZ0IiIDku1S7beC6en//8zvyDzL35P0hT+fh7+tI/ki+SMiF4nInonkS3H/FBE5oy5kSRL+IyLynoiMCej/FyJycbrkMeoeUz5GIr5X1QP97tnfAj+vgzB/CkQaNFX9raq+VAfhNjpUdUTUjt/DgUtU9ce497hSvW3+K+AvaRDpItxWK4nkyzgi0gE4RFX3V9U/BbztIdx7P0YjwZSPUR2v4zdYFZG9ROR5v8nkXBHZJ96ziFwiIm+J+37Ov/yb1Efg9ri63Y+o9gr3nMV9Q+nxqPsjIwYR+YmIvC4ib4vI437vuPj4XhWRu6JGatHbv/T01z8XkVFR9zzp07A4vFmm33x0ig/j/XCPO2Ca+8vOb7e8E94VAGghO7838w//RnlY5j5+FHMU8KDPg8uBMT6cfnFxdAe2hj9t4GW9V0Te8OkbIO6bT0tEZErUfRuj/p8RfS3shvucwz98vAVh+cL3i8iffF69LCLtEqS/t4jM9nn0grhPbCAio8R94+k9EXk0wX35IjLZ5/c7IvJjf+lFoGN8PohITxFZKSJf+WvfiMj/RKRIVTcDy+Oev9GQUVU77Ig5gI3+Nwu3KeJgf/4y0M3/74vbDgbgZuAa/3+3qHBuA37p/08Bzoi6NgU4A/fm/P9w39wBuBc4D2gLzIlyvw74bQJZXwXu9/+PBj6Ikmk+kOfDWgvk+Gtt/G8BbruY3YDeuB2vw+G2SpXmOBlmAEf6/y18mgYA3+H2CGyGU+JHRcncJ8H/SD4miGMYMCEu/x4FBLdR5HpgPx/XQuDA6Gfp/58BTEnwzCIyJJBJgXP9/98Cd8c9vxyfz+28+1m43UHAbdWfF52fcWm6OsrvPr4c5AOl4eeY4J5oue8GLoq69mvg6vquP3YEO2x7HSMRBeK2yemI2/dpph91HAE87jvw4Br2eHqJyG1AK1xD/EKqiNRta/Q8cLK4bT9OxJmX+uPMdK/5+HJxDXgiHvFhzRGRlrJzjuoZVd0KbBWR1bjtYlYAo0TkVO+nM9AN+BjoKiJ/Bp7BfUYgaJpfA/5PRP4BPKGqK7z//6rqCgCfn6W4D3nVhj1wnxCIZoaqqoi8D6xS1fd9XIt9XItqGVc0lcB0///vuI0oo+kB9MKVEXAdlpC/9h5uRPUk8GSCsI8C/gygqh+JyBdAd5wirQ2rcUrMaASY8jES8b2qHigizXHK4+e4nm65uk8TpGIK7ium74rIRbgRQHU8CvwCN7+0QFU3eBPVTFU9O8D98XtEhc+3RrntALJFZABuZ9/DVXWziLwK5KvqOhE5ADgOZ/46ExhNgDSr6jgReQY4Aacsj0sWf4C0JON7oDjOLRx+ZVxclVFxRedNXXwuOT6vBVisqok+YX0ibjR6MvBrEdlPVbfXgQzJyMflk9EIsDkfIynq7OijcOaRzcAyEfkZuJ2sfWMdTxEQEveJhXOj3Df4a4mYjdsO/xKcIgJ4AzhSRPb28RX6eY9EnOX9HIX7ONZ3KZJVDKzzimcf4DB/b1ugmar+C/gNcLC6bxJVm2YR2UtV31fVP+J2R69t7ztVHi0B9q5FmKtEZF8RaQacmsRPqnib4cxrAOdQdeT2MdBORA4H92kNEfmRj6+zqs7CmUyLcSPhaObiy4h/tl18eEHZjhsRh+lOA9l126geUz5GSlT1HZz55GxcQzFcRN4FFpP4M+M34r5i+hqxny54FLjWTyzvFRfHDuA/uO+f/Me7fYNbhfWIiLyHM7kla9S3iMg7wF9xq8dS8TxuBLQEGIdTcuBMjK9689jfgRu8e5A0jxa3UOE9oAJ4rhoZkjEDODXRggPc/NdBEmX/C8j1uDydz05zWDxTgL+GFxzEXduE+wDfB8BA4Jboi+o+OX8G8EefR4twpsos4O/eJPgOMFGrfub6L0Az72c6bv5mK8F5A7hNRH7kz48EZtbgfqMesV2tjUaNN5tdo6oL6luWdCMid+HmeTK2RF1ENqpqlVWGDQ0ROQi4SlXPr29ZjGDYyMcwGg9/AJrXtxANlLa4UbfRSLCRj2EYhpFxbORjGIZhZBxTPoZhGEbGMeVjGIZhZBxTPoZhGEbGMeVjGIZhZJz/D9qmf4960qAKAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"psi_0 = p45 # Defining the initial state (+45 vector)\n",
"\n",
"phaseshift2_init = pi / 4 # initial value\n",
"phaseshift2_max = 8 * pi\n",
"n = 100 # resolution of φ (amount of steps)\n",
"step = (phaseshift2_max - phaseshift2_init) / n # interval divided by number of small steps we want\n",
"N_init = 1000 # number of iterations (range(N) -> 0 to N-1, both included)\n",
"\n",
"x_coords = [] # create x- and y- coords. arrays (x = phase shift of 2nd PA_HV, y1 = amount of photons in +45, y2 = amount of photons in -45)\n",
"y1_coords = []\n",
"y2_coords = []\n",
"\n",
"for i in range(n + 1):\n",
" output_p45 = 0\n",
" output_n45 = 0\n",
"\n",
" phaseshift2 = phaseshift2_init + i * step\n",
" x_coords.append((phaseshift2 - phaseshift1) / pi) # add realtive phase shift to x coords\n",
" PA_HV2 = Qobj([[exp(1j * phaseshift2), 0], [0, 1]]) # create corresponding PA_HV2\n",
" EffM = PA_HV2 * halfwave * PA_HV1vb # Defining the effective matrix\n",
"\n",
" psi_final = EffM * psi_0 # Applying the effective matrix to the initial state to get the final state\n",
" psi_p45 = PA_45(psi_final)[0] # Determining the probabilities\n",
" psi_n45 = PA_45(psi_final)[1]\n",
" proba_p45 = round(psi_p45.norm()**2,5)\n",
" proba_n45 = round(psi_n45.norm()**2,5)\n",
"\n",
" if (proba_p45 + proba_n45) == 1:\n",
" N = N_init # all of the photons get to the end\n",
" else:\n",
" N = int(N_init / 2) # half of the photons are blocked (should only get N/2 in the ouput) -> total proba should be 0.5\n",
"\n",
" for j in range(N):\n",
" a = random.randint(1, 50) # generates random number between 1 and 50 (both included), 50 because 50% of the photons need to come out in either +45 or -45 state (since other 50% was blocked)\n",
" if a <= proba_p45 * 100:\n",
" output_p45 = output_p45 + 1\n",
" else:\n",
" output_n45 = output_n45 + 1\n",
"\n",
" y1_coords.append(output_p45)\n",
" y2_coords.append(output_n45)\n",
"\n",
"plt.plot(x_coords, y1_coords, 'b.', markersize=9, label='Photons in state +45')\n",
"plt.plot(x_coords, y2_coords, 'r.', markersize=9, label='Photons in state -45')\n",
"legend = plt.legend(loc='upper center', fontsize='x-large')\n",
"plt.ylim([0, 1000])\n",
"plt.xlabel('Relative phase shift (multiples of π)')\n",
"plt.ylabel('Amount of photons detected')\n",
"plt.title('Amount of photons existing in |+/-45> states (H port blocked)');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Conclusion and Takeaways "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The results show that, when both pathways are unobstructed, an interference pattern emerges. As the relative phase varies, the detected photons oscillate between the $|+45\\rangle$ and $|-45\\rangle$ states, clearly displaying interference. Note that, while inside the interferometer, the beams are in orthogonal polarizations, so no interference should occur there, though the components get affected and phase-shifted differently. Interference \"happens\", or rather becomes visible, at the PA$_{45}$ where both the horizontal and vertical polarisation are projected onto $+45^\\circ$ and $-45^\\circ$ axes. The degree of interference oscillates with increasing phase difference between the two paths.\n",
"\n",
"When one of the output port from the first polarization analyzer is blocked, the amount of photons detected is constant (apart from some random fluctuations) under a changing relative phase shift and there is no interference pattern. The results are identical for both setups, with the V beam and the H beam blocked. This illustrates that when the superposition collapses, no interference is displayed."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Important takeaways:**\n",
"* This experiment demonstrates that photons can behave as particles and waves at the same time. This illustrates that the wave-particle duality is not an either-or concept, and that the transition is fluid.\n",
"* *So, which path did the photon take before being 'reunited' with its other half?* It was mentioned that the photon never actually gets split nor reunited and remains in one (evolving) state throughout. However, the interference pattern is a function of phase length variance between the two arms of the interferometer, and thus a function of both paths. Hence there is no way to describe the measured interference other than by considering that the photon travels through both paths at the same time. This is done by considering the *wavefunction* of the photon. Its wavefunction is in a superposition of travelling down both paths simultaneously, and will remain in this state until it collapses due to a measurement or an other external factor that causes quantum decoherence. This superposition results in the observed interference pattern.\n",
"* By blocking one of the beam’s paths, the photon is prevented from being in superposition and therefore can't interfere with itself. It further results in only half the amount of photons exiting the interferometer. This is due to the photons originally having a 50/50 chance of passing through either path. We still observe photons in both polarization states with a 50/50 distribution, demonstrating that a quantum state can collapse from one superposition and still remain in another superposition at the same time. While the obstruction of one path permanently determines which path the photons are taking (thus collapsing the superposition that allowed for the interference), their polarisation is still in a superposition of two basis states afterwards.\n",
"* The simulation shows that single-particle quantum interference can be described using an *event-by-event* basis. The use of frameworks such as wave fields or time-dependent system evolutions is not necessary to demonstrate the basic workings of the physical mechanism. This is possible, as the mathematical essence of single photon interference can be described with simple vector and matrix functions. However, when intending to use the simulation to make actual predictions about physical mechanisms, it is important to take into account physical limitations to ensure the results of the simulation are meaningful. One example of such a limitation is coherence length, which puts a real limit on the relative phase differences that would still result in interference."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"QuTiP: Quantum Toolbox in Python\n",
"================================\n",
"Copyright (c) QuTiP team 2011 and later.\n",
"Current admin team: Alexander Pitchford, Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, Eric Giguère, Boxi Li, Jake Lishman and Simon Cross.\n",
"Board members: Daniel Burgarth, Robert Johansson, Anton F. Kockum, Franco Nori and Will Zeng.\n",
"Original developers: R. J. Johansson & P. D. Nation.\n",
"Previous lead developers: Chris Granade & A. Grimsmo.\n",
"Currently developed through wide collaboration. See https://github.com/qutip for details.\n",
"\n",
"QuTiP Version: 4.6.2\n",
"Numpy Version: 1.20.3\n",
"Scipy Version: 1.6.3\n",
"Cython Version: None\n",
"Matplotlib Version: 3.4.3\n",
"Python Version: 3.8.10\n",
"Number of CPUs: 4\n",
"BLAS Info: OPENBLAS\n",
"OPENMP Installed: False\n",
"INTEL MKL Ext: False\n",
"Platform Info: Linux (x86_64)\n",
"Installation path: /home/simon/venvs/py3/qutip-notebooks/lib/python3.8/site-packages/qutip\n",
"================================================================================\n",
"Please cite QuTiP in your publication.\n",
"================================================================================\n",
"For your convenience a bibtex reference can be easily generated using `qutip.cite()`\n"
]
}
],
"source": [
"qutip.about()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
![](\"./images/single-photon-interference-setup.jpg\")
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that a single photon is sent through these optical elements, making the process quantum mechanical in nature. This requires for an important distinction, when it comes to the splitting of the state into orthogonal components.\n",
"* A classical wave is split into its orthogonal components deterministically, with the wave actually getting divided into two components, each *in proportion* to the coefficient of the corresponding polarization state in the general polarization ($\\psi = c_H |H\\rangle + c_V |V\\rangle$).\n",
"* On the other hand, when it comes to single photons, as they cannot be broken down any further, they are 'split' randomly and it is not possible to determine with certainty which port any photon will emerge from.\n",
"* The proportion of the basis states $|H\\rangle$ and $|V\\rangle$ in a general state $\\psi$ describes the *probability* of the photon taking that corresponding path.\n",
"* Note also that the state of the photon never actually gets 'split' and it stays as one single state throughout the entire interferometer."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. Coding the constitutive elements "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Basis vectors\n",
"Code the polarization basis vectors, for the bases $H/V$ and $+/-45$. They can then be used to code any polarization state, since any state can be expressed as a linear combination of basis vectors."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# HV basis\n",
"\n",
"# horizontal polarization\n",
"H = basis(2, 0)\n",
"# vertical polarization\n",
"V = basis(2, 1)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# +45/-45 basis (in terms of HV basis)\n",
"\n",
"# +45 polarization\n",
"p45 = 1 / sqrt(2) * (H + V)\n",
"# -45 polarization\n",
"n45 = 1 / sqrt(2) * (H - V)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Polarization analyzer (HV) n°1\n",
"It induces a phase shift on the H polarization relative to V polarization, because of the difference in the path lengths of the two arms of the interferometer. \n",
"