diff --git a/_toc.yml b/_toc.yml index 3f46be09..42eb9934 100644 --- a/_toc.yml +++ b/_toc.yml @@ -35,6 +35,10 @@ parts: - file: chapter2/ns_code1 - file: chapter2/ns_code2 - file: chapter2/hyperelasticity + - file: chapter2/helmholtz + sections: + - file: chapter2/helmholtz_code + - caption: Subdomains and boundary conditions chapters: - file: chapter3/neumann_dirichlet_code @@ -48,4 +52,4 @@ parts: - file: chapter4/solvers - file: chapter4/compiler_parameters - file: chapter4/convergence - - file: chapter4/newton-solver \ No newline at end of file + - file: chapter4/newton-solver diff --git a/chapter2/helmholtz.md b/chapter2/helmholtz.md new file mode 100644 index 00000000..00d6d1ac --- /dev/null +++ b/chapter2/helmholtz.md @@ -0,0 +1,125 @@ +# The Helmholtz equation +Author: Antonio Baiano Svizzero + +The study of computational acoustics is fundamental in fields such as noise, vibration, and harshness (NVH), noise control, and acoustic design. In this chapter, we focus on the theoretical foundations of the Helmholtz equation - valid for noise problems with harmonic time dependency - and its implementation in FEniCSx to compute the sound pressure for any acoustic system. + +## The PDE problem +The acoustic Helmholtz equation in its general form reads + +$$ +\begin{align} +\nabla^2 p + k^2 p = -j \omega \rho_0 q \qquad\text{in } \Omega, +\end{align} +$$ + +where $k$ is the acoustic wavenumber, $\omega$ is the angular frequency, $j$ the imaginary unit and $q$ is the volume velocity ($m^3/s$) of a generic source field. +In case of a monopole source, we can write $q=Q \delta(x_s,y_s,z_s)$, where $\delta(x_s,y_s,z_s)$ is the 3D Dirac Delta centered at the monopole location. + +This equation is coupled with the following boundary conditions: + +- Dirichlet BC: + + $$ + \begin{align} + p = \bar{p} \qquad \text{on } \partial\Omega_p, + \end{align} + $$ + +- Neumann BC: + + $$ + \begin{align} + \frac{\partial p}{\partial n} = - j \omega \rho_0 \bar{v}_n\qquad \text{on } \partial\Omega_v, + \end{align} + $$ + +- Robin BC: + + $$ + \begin{align} + \frac{\partial p}{\partial n} = - \frac{j \omega \rho_0 }{\bar{Z}} p \qquad \text{on } \partial\Omega_Z, + \end{align} + $$ + +where we prescribe, respectively, an acoustic pressure $\bar{p}$ on the boundary $\partial\Omega_p$, +a sound particle velocity $\bar{v}_n$ on the boundary $\partial\Omega_v$ and +an acoustic impedance $\bar{Z}$ on the boundary $\partial\Omega_Z$ where $n$ is the outward normal. +In general, any BC can also be frequency dependant, as it happens in real-world applications. + +## The variational formulation +Now we have to turn the equation in its weak formulation. +The first step is to multiplicate the equation by a *test function* $v\in \hat V$, +where $\hat V$ is the *test function space*, after which we integrate over the whole domain, $\Omega$: + +$$ +\begin{align} +\int_{\Omega}\left(\nabla^2 p + k^2 p \right) \bar v ~\mathrm{d}x = -\int_{\Omega} j \omega \rho_0 q \bar v ~\mathrm{d}x. +\end{align} +$$ + +Here, the unknown function $p$ is referred to as *trial function* and the $\bar{\cdot}$ is the complex conjugate operator. + +In order to keep the order of derivatives as low as possible, we use integration by parts on the Laplacian term: + +$$ +\begin{align} +\int_{\Omega}(\nabla^2 p) \bar v ~\mathrm{d}x = +-\int_{\Omega} \nabla p \cdot \nabla \bar v ~\mathrm{d}x ++ \int_{\partial \Omega} \frac{\partial p}{\partial n} \bar v ~\mathrm{d}s. +\end{align} +$$ + +Substituting in the original version and rearranging we get: + +$$ +\begin{align} +\int_{\Omega} \nabla p \cdot \nabla \bar v ~\mathrm{d}x +- k^2 \int_{\Omega} p \bar v ~\mathrm{d} x = \int_{\Omega} j \omega \rho_0 q \bar v ~\mathrm{d}x ++ \int_{\partial \Omega} \frac{\partial p}{\partial n} \bar v ~\mathrm{d}s. +\end{align} +$$ + +Since we are dealing with complex values, the inner product in the first equation is *sesquilinear*, +meaning it is linear in one argument and conjugate-linear in the other, +as explained in [The Poisson problem with complex numbers](../chapter1/complex_mode). + +The last term can be written using the Neumann and Robin BCs, that is: + +$$ +\begin{align} +\int_{\partial \Omega} \frac{\partial p}{\partial n} \bar v ~\mathrm{d}s = +-\int_{\partial \Omega_v} j \omega \rho_0 \bar v ~\mathrm{d}s +- \int_{\partial \Omega_Z} \frac{j \omega \rho_0 \bar{v}_n}{\bar{Z}} p \bar v ~\mathrm{d}s. +\end{align} +$$ + +Substituting, rearranging and taking out of integrals the terms with $j$ and $\omega$ we get the variational formulation of the Helmholtz. +Find $u \in V$ such that: + +$$ +\begin{align} +\int_{\Omega} \nabla p \cdot \nabla \bar v ~\mathrm{d}x ++ \frac{j \omega }{\bar{Z}} \int_{\partial \Omega_Z} \rho_0 p \bar v ~\mathrm{d}s +- k^2 \int_{\Omega} p \bar v ~\mathrm{d}x += j \omega \int_{\Omega} \rho_0 q \bar v ~\mathrm{d}x +-j \omega\int_{\partial \Omega_v} \rho_0 \bar{v}_n \bar v ~\mathrm{d}s \qquad \forall v \in \hat{V}. +\end{align} +$$ + +We define the sesquilinear form $a(p,v)$ is + +$$ +\begin{align} +a(p,v) = \int_{\Omega} \nabla p \cdot \nabla \bar v ~\mathrm{d}x ++ \frac{j \omega }{\bar{Z}} \int_{\partial \Omega_Z} \rho_0 p \bar v ~\mathrm{d}s, +- k^2 \int_{\Omega} p \bar v ~\mathrm{d}x +\end{align} +$$ + +and the linear form $L(v)$ reads + +$$ +\begin{align} +L(v) = j \omega \int_{\Omega}\rho_0 q \bar v ~\mathrm{d}x - j \omega \int_{\partial \Omega_v} \rho_0 \bar{v}_n \bar v ~\mathrm{d}s. +\end{align} +$$ diff --git a/chapter2/helmholtz_code.ipynb b/chapter2/helmholtz_code.ipynb new file mode 100644 index 00000000..f85949bb --- /dev/null +++ b/chapter2/helmholtz_code.ipynb @@ -0,0 +1,519 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Implementation\n", + "Author: Antonio Baiano Svizzero and Jørgen S. Dokken\n", + "\n", + "In this tutorial, you will learn how to:\n", + "- Define acoustic velocity and impedance boundary conditions\n", + "- Compute acoustic sound pressure for multiple frequencies\n", + "- Compute the Sound Pressure Level (SPL) at a given microphone position\n", + "\n", + "## Test problem\n", + "As an example, we will model a plane wave propagating in a tube.\n", + "While it is a basic test case, the code can be adapted to way more complex problems where velocity and impedance boundary conditions are needed.\n", + "We will apply a velocity boundary condition $v_n = 0.001$ to one end of the tube and an impedance $Z$ computed with the Delaney-Bazley model,\n", + "supposing that a layer of thickness $d = 0.02$ and flow resistivity $\\sigma = 1e4$ is placed at the second end of the tube.\n", + "The choice of such impedance (the one of a plane wave propagating in free field) will give, as a result, a solution with no reflections.\n", + "\n", + "First, we create the mesh with gmsh, also setting the physical group for velocity and impedance boundary conditions and the respective tags." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Info : Meshing 1D...\n", + "Info : [ 0%] Meshing curve 1 (Line)\n", + "Info : [ 10%] Meshing curve 2 (Line)\n", + "Info : [ 20%] Meshing curve 3 (Line)\n", + "Info : [ 30%] Meshing curve 4 (Line)\n", + "Info : [ 40%] Meshing curve 5 (Line)\n", + "Info : [ 50%] Meshing curve 6 (Line)\n", + "Info : [ 60%] Meshing curve 7 (Line)\n", + "Info : [ 60%] Meshing curve 8 (Line)\n", + "Info : [ 70%] Meshing curve 9 (Line)\n", + "Info : [ 80%] Meshing curve 10 (Line)\n", + "Info : [ 90%] Meshing curve 11 (Line)\n", + "Info : [100%] Meshing curve 12 (Line)\n", + "Info : Done meshing 1D (Wall 0.000926869s, CPU 0.000589s)\n", + "Info : Meshing 2D...\n", + "Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay)\n", + "Info : [ 20%] Meshing surface 2 (Plane, Frontal-Delaunay)\n", + "Info : [ 40%] Meshing surface 3 (Plane, Frontal-Delaunay)\n", + "Info : [ 60%] Meshing surface 4 (Plane, Frontal-Delaunay)\n", + "Info : [ 70%] Meshing surface 5 (Plane, Frontal-Delaunay)\n", + "Info : [ 90%] Meshing surface 6 (Plane, Frontal-Delaunay)\n", + "Info : Done meshing 2D (Wall 0.0311526s, CPU 0.030342s)\n", + "Info : Meshing 3D...\n", + "Info : 3D Meshing 1 volume with 1 connected component\n", + "Info : Tetrahedrizing 1285 nodes...\n", + "Info : Done tetrahedrizing 1293 nodes (Wall 0.0127086s, CPU 0.011753s)\n", + "Info : Reconstructing mesh...\n", + "Info : - Creating surface mesh\n", + "Info : - Identifying boundary edges\n", + "Info : - Recovering boundary\n", + "Info : Done reconstructing mesh (Wall 0.0264354s, CPU 0.02449s)\n", + "Info : Found volume 1\n", + "Info : It. 0 - 0 nodes created - worst tet radius 2.85527 (nodes removed 0 0)\n", + "Info : 3D refinement terminated (1746 nodes total):\n", + "Info : - 0 Delaunay cavities modified for star shapeness\n", + "Info : - 0 nodes could not be inserted\n", + "Info : - 6588 tetrahedra created in 0.0199287 sec. (330577 tets/s)\n", + "Info : 0 node relocations\n", + "Info : Done meshing 3D (Wall 0.0772706s, CPU 0.07553s)\n", + "Info : Optimizing mesh...\n", + "Info : Optimizing volume 1\n", + "Info : Optimization starts (volume = 0.01) with worst = 0.0138257 / average = 0.769909:\n", + "Info : 0.00 < quality < 0.10 : 18 elements\n", + "Info : 0.10 < quality < 0.20 : 69 elements\n", + "Info : 0.20 < quality < 0.30 : 66 elements\n", + "Info : 0.30 < quality < 0.40 : 67 elements\n", + "Info : 0.40 < quality < 0.50 : 150 elements\n", + "Info : 0.50 < quality < 0.60 : 273 elements\n", + "Info : 0.60 < quality < 0.70 : 913 elements\n", + "Info : 0.70 < quality < 0.80 : 1842 elements\n", + "Info : 0.80 < quality < 0.90 : 2090 elements\n", + "Info : 0.90 < quality < 1.00 : 1095 elements\n", + "Info : 152 edge swaps, 0 node relocations (volume = 0.01): worst = 0.300511 / average = 0.784453 (Wall 0.00193693s, CPU 0.00201s)\n", + "Info : No ill-shaped tets in the mesh :-)\n", + "Info : 0.00 < quality < 0.10 : 0 elements\n", + "Info : 0.10 < quality < 0.20 : 0 elements\n", + "Info : 0.20 < quality < 0.30 : 0 elements\n", + "Info : 0.30 < quality < 0.40 : 67 elements\n", + "Info : 0.40 < quality < 0.50 : 135 elements\n", + "Info : 0.50 < quality < 0.60 : 271 elements\n", + "Info : 0.60 < quality < 0.70 : 907 elements\n", + "Info : 0.70 < quality < 0.80 : 1872 elements\n", + "Info : 0.80 < quality < 0.90 : 2113 elements\n", + "Info : 0.90 < quality < 1.00 : 1081 elements\n", + "Info : Done optimizing mesh (Wall 0.00743052s, CPU 0.007641s)\n", + "Info : 1746 nodes 9265 elements\n" + ] + } + ], + "source": [ + "import gmsh\n", + "\n", + "gmsh.initialize()\n", + "\n", + "# meshsize settings\n", + "meshsize = 0.02\n", + "gmsh.option.setNumber(\"Mesh.MeshSizeMax\", meshsize)\n", + "gmsh.option.setNumber(\"Mesh.MeshSizeMax\", meshsize)\n", + "\n", + "\n", + "# create geometry\n", + "L = 1\n", + "W = 0.1\n", + "\n", + "gmsh.model.occ.addBox(0, 0, 0, L, W, W)\n", + "gmsh.model.occ.synchronize()\n", + "\n", + "# setup physical groups\n", + "v_bc_tag = 2\n", + "Z_bc_tag = 3\n", + "gmsh.model.addPhysicalGroup(3, [1], 1, \"air_volume\")\n", + "gmsh.model.addPhysicalGroup(2, [1], v_bc_tag, \"velocity_BC\")\n", + "gmsh.model.addPhysicalGroup(2, [2], Z_bc_tag, \"impedance\")\n", + "\n", + "# mesh generation\n", + "gmsh.model.mesh.generate(3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we import the gmsh mesh with the ```dolfinx.io.gmshio``` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from mpi4py import MPI\n", + "from dolfinx import fem, io, default_scalar_type, geometry\n", + "from dolfinx.fem.petsc import LinearProblem\n", + "import ufl\n", + "import numpy as np\n", + "import numpy.typing as npt\n", + "\n", + "mesh_data = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)\n", + "domain = mesh_data.mesh\n", + "assert mesh_data.facet_tags is not None\n", + "facet_tags = mesh_data.facet_tags" + ] + }, + { + "cell_type": "markdown", + "id": "92d8cd12", + "metadata": {}, + "source": [ + "We define the function space for our unknown $p$ and define the range of frequencies we want to solve the Helmholtz equation for." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06cebc38", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "V = fem.functionspace(domain, (\"Lagrange\", 1))\n", + "\n", + "# Discrete frequency range\n", + "freq = np.arange(10, 1000, 5) # Hz\n", + "\n", + "# Air parameters\n", + "rho0 = 1.225 # kg/m^3\n", + "c = 340 # m/s\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "## Boundary conditions\n", + "\n", + "The Delaney-Bazley model is used to compute the characteristic impedance and wavenumber of the porous layer,\n", + "treated as an equivalent fluid with complex valued properties\n", + "\n", + "$$\n", + "\\begin{align}\n", + "Z_c(\\omega) &= \\rho_0 c_0 \\left[1 + 0.0571 X^{-0.754} - j 0.087 X^{-0.732}\\right],\\\\\n", + "k_c(\\omega) &= \\frac{\\omega}{c_0} \\left[1 + 0.0978 X^{-0.700} - j 0.189 X^{-0.595}\\right],\\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "where $X = \\frac{\\rho_0 f}{\\sigma}$.\n", + "\n", + "With these, we can compute the surface impedance, that in the case of a rigid passive absorber placed on a rigid wall is given by the formula\n", + "$$\n", + "\\begin{align}\n", + "Z_s = -j Z_c cot(k_c d).\n", + "\\end{align}\n", + "$$\n", + "\n", + "Let's create a function to compute it.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Impedance calculation\n", + "def delany_bazley_layer(f, rho0, c, sigma):\n", + " X = rho0 * f / sigma\n", + " Zc = rho0 * c * (1 + 0.0571 * X**-0.754 - 1j * 0.087 * X**-0.732)\n", + " kc = 2 * np.pi * f / c * (1 + 0.0978 * (X**-0.700) - 1j * 0.189 * (X**-0.595))\n", + " Z_s = -1j * Zc * (1 / np.tan(kc * d))\n", + " return Z_s\n", + "\n", + "\n", + "sigma = 1.5e4\n", + "d = 0.01\n", + "Z_s = delany_bazley_layer(freq, rho0, c, sigma)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since we are going to compute a sound pressure spectrum, all the variables that depend on frequency (that are $\\omega$, $k$ and $Z$) need to be updated in the frequency loop.\n", + "To make this possible, we will initialize them as dolfinx constants.\n", + "Then, we define the value for the normal velocity on the first end of the tube" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "omega = fem.Constant(domain, default_scalar_type(0))\n", + "k = fem.Constant(domain, default_scalar_type(0))\n", + "Z = fem.Constant(domain, default_scalar_type(0))\n", + "v_n = 1e-5" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We also need to specify the integration measure $ds$, by using ```ufl```, and its built in integration measures" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "ds = ufl.Measure(\"ds\", domain=domain, subdomain_data=facet_tags)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Variational Formulation\n", + "We can now write the variational formulation." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "p = ufl.TrialFunction(V)\n", + "v = ufl.TestFunction(V)\n", + "\n", + "a = (\n", + " ufl.inner(ufl.grad(p), ufl.grad(v)) * ufl.dx\n", + " + 1j * rho0 * omega / Z * ufl.inner(p, v) * ds(Z_bc_tag)\n", + " - k**2 * ufl.inner(p, v) * ufl.dx\n", + ")\n", + "L = -1j * omega * rho0 * ufl.inner(v_n, v) * ds(v_bc_tag)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The class ```LinearProblem``` is used to setup the PETSc backend and assemble the system vector and matrices.\n", + "The solution will be stored in a `dolfinx.fem.Function`, ```p_a```." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "lines_to_end_of_cell_marker": 2 + }, + "outputs": [], + "source": [ + "p_a = fem.Function(V)\n", + "p_a.name = \"pressure\"\n", + "\n", + "problem = LinearProblem(\n", + " a,\n", + " L,\n", + " u=p_a,\n", + " petsc_options={\n", + " \"ksp_type\": \"preonly\",\n", + " \"pc_type\": \"lu\",\n", + " \"pc_factor_mat_solver_type\": \"mumps\",\n", + " },\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "## Computing the pressure at a given location\n", + "Before starting our frequency loop, we can build a function that, given a microphone position,\n", + "computes the sound pressure at its location.\n", + "We will use the a similar method as in [Deflection of a membrane](../chapter1/membrane_code).\n", + "However, as the domain doesn't deform in time, we cache the collision detection" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "class MicrophonePressure:\n", + " def __init__(self, domain, microphone_position):\n", + " \"\"\"Initialize microphone(s).\n", + "\n", + " Args:\n", + " domain: The domain to insert microphones on\n", + " microphone_position: Position of the microphone(s).\n", + " Assumed to be ordered as ``(mic0_x, mic1_x, ..., mic0_y, mic1_y, ..., mic0_z, mic1_z, ...)``\n", + "\n", + " \"\"\"\n", + " self._domain = domain\n", + " self._position = np.asarray(\n", + " microphone_position, dtype=domain.geometry.x.dtype\n", + " ).reshape(3, -1)\n", + " self._local_cells, self._local_position = self.compute_local_microphones()\n", + "\n", + " def compute_local_microphones(\n", + " self,\n", + " ) -> tuple[npt.NDArray[np.int32], npt.NDArray[np.floating]]:\n", + " \"\"\"\n", + " Compute the local microphone positions for a distributed mesh\n", + "\n", + " Returns:\n", + " Two lists (local_cells, local_points) containing the local cell indices and the local points\n", + " \"\"\"\n", + " points = self._position.T\n", + " bb_tree = geometry.bb_tree(self._domain, self._domain.topology.dim)\n", + "\n", + " cells = []\n", + " points_on_proc = []\n", + "\n", + " cell_candidates = geometry.compute_collisions_points(bb_tree, points)\n", + " colliding_cells = geometry.compute_colliding_cells(\n", + " domain, cell_candidates, points\n", + " )\n", + "\n", + " for i, point in enumerate(points):\n", + " if len(colliding_cells.links(i)) > 0:\n", + " points_on_proc.append(point)\n", + " cells.append(colliding_cells.links(i)[0])\n", + "\n", + " return np.asarray(cells, dtype=np.int32), np.asarray(\n", + " points_on_proc, dtype=domain.geometry.x.dtype\n", + " )\n", + "\n", + " def listen(\n", + " self, recompute_collisions: bool = False\n", + " ) -> npt.NDArray[np.complexfloating]:\n", + " if recompute_collisions:\n", + " self._local_cells, self._local_position = self.compute_local_microphones()\n", + " if len(self._local_cells) > 0:\n", + " return p_a.eval(self._local_position, self._local_cells)\n", + " else:\n", + " return np.zeros(0, dtype=default_scalar_type)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The pressure spectrum is initialized as a numpy array and the microphone location is assigned" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "p_mic = np.zeros((len(freq), 1), dtype=complex)\n", + "\n", + "mic = np.array([0.5, 0.05, 0.05])\n", + "microphone = MicrophonePressure(domain, mic)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Frequency loop\n", + "\n", + "Finally, we can write the frequency loop, where we update the values of the frequency-dependent variables and solve the system for each frequency" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "for nf in range(0, len(freq)):\n", + " k.value = 2 * np.pi * freq[nf] / c\n", + " omega.value = 2 * np.pi * freq[nf]\n", + " Z.value = Z_s[nf]\n", + "\n", + " problem.solve()\n", + " p_a.x.scatter_forward()\n", + "\n", + " p_f = microphone.listen()\n", + " p_f = domain.comm.gather(p_f, root=0)\n", + "\n", + " if domain.comm.rank == 0:\n", + " assert p_f is not None\n", + " p_mic[nf] = np.hstack(p_f)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## SPL spectrum\n", + "After the computation, the pressure spectrum at the prescribed location is available.\n", + "Such a spectrum is usually shown using the decibel (dB) scale to obtain the SPL, with the RMS pressure as input,\n", + "defined as $p_{rms} = \\frac{p}{\\sqrt{2}}$." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAB9QAAAKsCAYAAAC9C4YeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAD6DElEQVR4nOzdd3hUddrG8Xtm0jsJpEFC772DgBVR7L279vKqq4iuK7u6rrrq6hZ7XRXsvWEvKChKk947pBfSeyYz8/4xyWQCARJSzpTv57q8rjP9Acmck3Of3/OYHA6HQwAAAAAAAAAAAAAAoAmz0QUAAAAAAAAAAAAAAOCJCNQBAAAAAAAAAAAAAGgGgToAAAAAAAAAAAAAAM0gUAcAAAAAAAAAAAAAoBkE6gAAAAAAAAAAAAAANINAHQAAAAAAAAAAAACAZhCoAwAAAAAAAAAAAADQDAJ1AAAAAAAAAAAAAACaQaAOAAAAAAAAAAAAAEAzCNQBAAAAAAAAAAAAAGiGoYF6WVmZZs2apZ49eyo0NFRHHXWUVqxY4Xrc4XDob3/7m5KSkhQaGqrp06dr+/btBlYMAAAAAAAAAAAAAPAXhgbq1157rb7//nu98cYbWr9+vWbMmKHp06crMzNTkvTYY4/pqaee0gsvvKBly5YpPDxcJ510kqqrq40sGwAAAAAAAAAAAADgB0wOh8NhxAdXVVUpMjJSn332mU499VTX/WPHjtXMmTP14IMPKjk5WXfccYfuvPNOSVJJSYkSEhI0b948XXTRRUaUDQAAAAAAAAAAAADwEwFGfXBdXZ1sNptCQkKa3B8aGqrFixdr9+7dysnJ0fTp012PRUdHa+LEiVqyZMlBA/WamhrV1NS4btvtdhUWFiouLk4mk6lj/jAAAAAAAAAAAAAAAK/gcDhUVlam5ORkmc2HbupuWKAeGRmpyZMn68EHH9TgwYOVkJCgd955R0uWLFG/fv2Uk5MjSUpISGjyuoSEBNdjzXnkkUd0//33d2jtAAAAAAAAAAAAAADvlp6erh49ehzyOYYF6pL0xhtv6Oqrr1b37t1lsVg0ZswYXXzxxVq5cuURv+ecOXM0e/Zs1+2SkhKlpqZq9+7dioyMbI+yAQAGslqt+umnn3TccccpMDDQ6HIAAB6EfQQA4GDYRwAADoX9BAD4n7KyMvXu3btF+bGhgXrfvn21aNEiVVRUqLS0VElJSbrwwgvVp08fJSYmSpJyc3OVlJTkek1ubq5GjRp10PcMDg5WcHDwAffHxsYqKiqq3f8MAIDOZbVaFRYWpri4OH7BAQA0wT4CAHAw7CMAAIfCfgIA/E/D931LRoYfuiF8JwkPD1dSUpKKior07bff6swzz1Tv3r2VmJioBQsWuJ5XWlqqZcuWafLkyQZWCwAAAAAAAAAAAADwB4auUP/222/lcDg0cOBA7dixQ3/60580aNAgXXXVVTKZTJo1a5b+8Y9/qH///urdu7fuvfdeJScn66yzzjKybAAAAAAAAAAAAACAHzA0UC8pKdGcOXOUkZGh2NhYnXvuuXrooYdcS+zvuusuVVRU6Prrr1dxcbGmTp2qb775RiEhIUaWDQAAAAAAAAAAAADwA4YG6hdccIEuuOCCgz5uMpn0wAMP6IEHHujEqgAAAAAAAAAAAAAAnsrhcKiurk42m+2gzwkMDJTFYmnzZxkaqAMAAAAAAAAAAAAA0FK1tbXKzs5WZWXlIZ9nMpnUo0cPRUREtOnzCNQBAAAAAAAAAAAAAB7Pbrdr9+7dslgsSk5OVlBQkEwm0wHPczgcys/PV0ZGhvr379+mleoE6gAAAAAAAAAAAAAAj1dbWyu73a6UlBSFhYUd8rndunXTnj17ZLVa2xSom4/4lQAAAAAAAAAAAAAAdDKz+fAxd3Mr14/os9rlXQAAAAAAAAAAAAAA8DEE6gAAAAAAAAAAAAAANINAHQAAAAAAAAAAAACAZhCoAwAAAAAAAAAAAADQDAJ1AAAAAAAAAAAAAIDXcDgc7fKcliBQBwAAAAAAAAAAAAB4vMDAQElSZWXlYZ9bW1srSbJYLG36zIA2vRoAAAAAAAAAAAAAgE5gsVgUExOjvLw8SVJYWJhMJtMBz7Pb7crPz1dYWJgCAtoWiROoAwAAAAAAAAAAAAC8QmJioiS5QvWDMZvNSk1NbTZwbw0CdQAAAAAAAAAAAACAVzCZTEpKSlJ8fLysVutBnxcUFCSzue0T0AnUAQAAAAAAAAAAAABexWKxtHk+eku0PZIHAAAAAAAAAAAAAMAHEagDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJpBoA4AAAAAAAAAAAAAQDMI1AEAAAAAAAAAAAAAaAaBOgAAAAAAAAAAAAAAzSBQBwAAAAAAAAAAAACgGQTqAAAAAAAAAAAAAAA0g0AdAAAAAAAAAAAAAIBmEKgDAAAAAAAAAAAAANAMAnUAAAAAAAAAAAAAAJphaKBus9l07733qnfv3goNDVXfvn314IMPyuFwuJ7jcDj0t7/9TUlJSQoNDdX06dO1fft2A6sGAAAAAAAAAAAAAPgDQwP1Rx99VM8//7yeeeYZbd68WY8++qgee+wxPf30067nPPbYY3rqqaf0wgsvaNmyZQoPD9dJJ52k6upqAysHAAAAAAAAAAAAAPi6ACM//LffftOZZ56pU089VZLUq1cvvfPOO1q+fLkk5+r0J554Qvfcc4/OPPNMSdLrr7+uhIQEffrpp7rooosMqx0AAAAAAAAAAAAA4NsMDdSPOuoovfTSS9q2bZsGDBigtWvXavHixfrvf/8rSdq9e7dycnI0ffp012uio6M1ceJELVmypNlAvaamRjU1Na7bpaWlkiSr1Sqr1drBfyIAQEdr+C7nOx0AsD/2EQCAg2EfAQA4FPYTAOB/WvOdb2igfvfdd6u0tFSDBg2SxWKRzWbTQw89pEsvvVSSlJOTI0lKSEho8rqEhATXY/t75JFHdP/99x9w/3fffaewsLB2/hMAAIzy/fffG10CAMBDsY8AABwM+wgAwKGwnwAA/1FZWdni5xoaqL///vt666239Pbbb2vo0KFas2aNZs2apeTkZF1xxRVH9J5z5szR7NmzXbdLS0uVkpKiGTNmKCoqqr1KBwAYxGq16vvvv9eJJ56owMBAo8sBAHgQ9hEAgINhHwEAOBT2EwDgfxq6nLeEoYH6n/70J919992u1u3Dhw/X3r179cgjj+iKK65QYmKiJCk3N1dJSUmu1+Xm5mrUqFHNvmdwcLCCg4MPuD8wMJAdIQD4EL7XAQAHwz4CAHAw7CMAAIfCfgIA/Edrvu/NHVjHYVVWVspsblqCxWKR3W6XJPXu3VuJiYlasGCB6/HS0lItW7ZMkydP7tRaAQAAAAAAAAAAAAD+xdAV6qeffroeeughpaamaujQoVq9erX++9//6uqrr5YkmUwmzZo1S//4xz/Uv39/9e7dW/fee6+Sk5N11llnGVk6AAAAAAAAAAAAAMDHGRqoP/3007r33nt10003KS8vT8nJybrhhhv0t7/9zfWcu+66SxUVFbr++utVXFysqVOn6ptvvlFISIiBlQMAAAAAAAAAAAAAfJ2hgXpkZKSeeOIJPfHEEwd9jslk0gMPPKAHHnig8woDAAAAAAAAAAAAAPg9Q2eoAwAAAAAAAAAAAADgqQjUAQAAAAAAAAAAAABoBoE6AAAAAAAAAAAAAADNIFAHAAAAAAAAAAAAAKAZBOoAAAAAAAAAAAAAADSDQB0AAAAAAAAAAAAAgGYQqAMAAAAAAAAAAAAA0AwCdQAAAAAAAAAAAAAAmkGgDgAAAAAAAAAAAABAMwjUAQAAAAAAAAAAAABoBoE6AAAAAAAAAAAAAADNIFAHAAAAAAAAAAAAAKAZBOoAAAAAAAAAAAAAADSDQB0AAAA+weFw6NZ3VuvkJ37Wjrwyo8sBAAAAAAAA4AMI1AEAAOAT1mWUaP7aLG3JKdMbS/YaXQ4AAAAAAAAAH0CgDgAAAJ+QVVzl2t5dUGlgJQAAAAAAAAB8BYE6AAAAfEJeWY1rO72QQB0AAAAAAABA2xGoAwAAwCfklVW7tjOKKmWzOwysBgAAAAAAAIAvIFAHAACAT8grbVyhbrU5lFNafYhnAwAAAAAAAMDhEagDAADAJ+S6tXyXpDTmqAMAAAAAAABoIwJ1AAAA+IS8/VakM0cdAAAAAAAAQFsRqAMAAMAn5O+/Qp1AHQAAAAAAAEAbEagDAADA61ltdhVU1Da5j0AdAAAAAAAAQFsRqAMAAMDr7SuvOeC+vQTqAAAAAAAAANqIQB0AAABeL6/0wECdGeoAAAAAAAAA2opAHQAAAF4vr+zAQL2wolZl1VYDqgEAAAAAAADgKwjUAQAA4PXyyqpd24EWk2s7vbDKiHIAAAAAAAAA+AgCdQAAAHg995bvQ5OjXdtptH0HAAAAAAAA0AYE6gAAAPB67ivUx/Xs4tpmjjoAAAAAAACAtiBQBwAAgNdzX6E+rlesa5sV6gAAAAAAAADagkAdAAAAXi+vzBmom03SmNQY1/0E6gAAAAAAAADagkAdAAAAXq+h5XvXiGB1iwxWWJBFEi3fAQAAAAAAALQNgToAAAC8ms3u0L7yWklSfFSwTCaTUmPDJEkZRVWy2R1GlgcAAAAfUlNn098+26BHvtrMcSYAAICfIFAHAACAVyusqHWdzIyPDJEkpdQH6rU2u3JLqw2rDQAAAL7lo5WZen3JXr348y4t3JpndDkAAADoBATqAAAA8GoN7d4lKT4yWJJcK9Ql5qgDAACg/axJL3Jtb88rN7ASAAAAdBYCdQAAAHi1vLIa1zaBOgAAADrS5uwy13ZWcZWBlQAAAKCzEKgDAADAq+W5tXTvFuVs+e4eqKcTqAMAAKAd1Nns2ppLoA4AAOBvCNQBAADg1fJKD1yhnuIWqO8tIFAHAABA2+3aV6HaOrvrdkYRgToAAIA/IFAHAACAV3Nv+Z5Qv0K9R5dQ1320fAcAAEB72Jxd2uQ2K9QBAAD8A4E6AAAAvFpeWWPL94YV6iGBFiXWh+u0fAcAAEB72LRfoF5aXaeyaqtB1QAAAKCzEKgDAADAq7mvUO8aEezabpijXlBRq/Kauk6vCwAAAL5lc3bZAfdlFVc380wAAAD4EgJ1AAAAeLWGGeqx4UEKCmg8vHWfo84qdQAAALTV/i3fJdq+AwAA+AMCdQAAAHgth8Oh/PoV6g3t3hukugXqzFEHAABAW+wrr3Edd7rLJFAHAADweQTqAAAA8FolVVbV2uySpG77B+pxoa5tVqgDAACgLdxXp/ePj3BtE6gDAAD4PgJ1AAAAeC33+enxkSFNHkuNDXdts0IdAAAAbeEeqE8fkuDapuU7AACA7yNQBwAAgNfKLa12bcdH0fIdAAAAHWNzdplr+4RB8a5tAnUAAADfR6AOAAAAr5VX6r5CvWmg3jUiSKGBFkkE6gAAAGibhhXqAWaThveIVmx4kCQpq7j6UC8DAACADyBQBwAAgNdyb/meENW05bvJZHKtUs8orJLd7ujU2gAAAOAbaups2pFXLknq2y1CwQEWdY8JlSTllFarzmY3sjwAAAB0MAJ1AAAAeK28MreW7/utUJeklPpAvdZmV24Zq4cAAADQettzy1VXf3Hm4KRISVJyjPNiTpvdoVy3izwBAADgewjUAQAA4LXcV6jHR4Yc8HiTOeoFtH0HAABA6zW0e5ekwUlRkqTk+hXqEnPUAQAAfB2BOgAAALxWvvsM9agDV6inxjae6GSOOgAAAI7E5uwy13ZDoN7dLVDPLCJQBwAA8GUE6gAAAPBaDS3fI0MCFBJoOeDx1Di3FeoE6gAAADgCza1QbxKos0IdAADApxGoAwAAwCs5HA5Xy/fm5qdL+7V8J1AHAABAKzkcDm3OcQbqXSOC1a3+uJOW7wAAAP6DQB0AAABeqbymTpW1NknNz0+XpB5dCNQBAABw5HJKq1VcaZUkDUmOct3fvQsr1AEAAPwFgToAAAC8UsPqdKn5+emSFBJoUUL9Y+kE6gAAAGilpu3eI13bceFBCgpwnlplhToAAIBvI1AHAACAV8ordQvUD9LyXWps+76vvFYVNXUdXhcAAAB8x+bsMtf2kKTGFeomk8k1Rz2zqEoOh6PTawMAAEDnIFAHAACAV8orq3ZtJ0Q13/JdklLc5qinF7FKHQAAAC23qckK9agmjyXHOI9BK2ptKq3iwk0AAABfRaAOAAAAr5Tv1vK9WwtWqEtSWgGBOgAAAFquoeV7UIBZfbqGN3msYYW6xBx1AAAAX2ZooN6rVy+ZTKYD/rv55pslSdXV1br55psVFxeniIgInXvuucrNzTWyZAAAAHiIJjPUIw++Qr1JoM4cdQAAALRQVa1Ne/ZVSJIGJEQowNL0VGqyW6DOHHUAAADfZWigvmLFCmVnZ7v++/777yVJ559/viTp9ttv1+eff64PPvhAixYtUlZWls455xwjSwYAAICHyCttbPkeH3XwFeo949xavhOoAwAAoIW25pbJXj8afXBi1AGPJ7NCHQAAwC8EGPnh3bp1a3L7n//8p/r27atjjjlGJSUleuWVV/T222/r+OOPlyTNnTtXgwcP1tKlSzVp0iQjSgYAAICHaLpC/eCBegor1AEAAHAENmUdfH66JPVghToAAIBfMDRQd1dbW6s333xTs2fPlslk0sqVK2W1WjV9+nTXcwYNGqTU1FQtWbLkoIF6TU2NamoaT66WljoPfK1Wq6xWa8f+IQAAHa7hu5zvdAC59SvUQwPNCjY7Dvq9EBNsVkigWdVWu/YWVPL94cPYRwAADoZ9BI7Exsxi1/aA+LAD/v3ERwS6tjMKOc4EvBn7CQDwP635zveYQP3TTz9VcXGxrrzySklSTk6OgoKCFBMT0+R5CQkJysnJOej7PPLII7r//vsPuP+7775TWFhYM68AAHijhjEhAPxXVqFFkknhZpu+/vrrQz43JsCiHKtJaQXl+uLLr2Q2dU6NMAb7CADAwbCPQGv8ttl5vClJ6euXqmBz08etdqnh9OrGPdn66quMTq0PQPtjPwEA/qOysuWdLD0mUH/llVc0c+ZMJScnt+l95syZo9mzZ7tul5aWKiUlRTNmzFBU1IGtmQAA3sVqter777/XiSeeqMDAwMO/AIBPqrbaVLVkgSSpV2IXnXLKhEM+/7PC1crZmq86h0njph2vxKiQzigTnYx9BADgYNhHoLXsdof+supHSTYlR4fovDOObvZ5j25cqPzyWlWZQnTKKcd0bpEA2g37CQDwPw1dzlvCIwL1vXv36ocfftDHH3/sui8xMVG1tbUqLi5usko9NzdXiYmJB32v4OBgBQcfOEMzMDCQHSEA+BC+1wH/ll3a2JIpISr0sN8HPbuGS1vzXa9NiYvs0PpgLPYRAICDYR+BlkorqFRFjU2SNCQ56qD/bpK7hCm/vFZ55TVymCwKCjB3ZpkA2hn7CQDwH635vveII7y5c+cqPj5ep556quu+sWPHKjAwUAsWLHDdt3XrVqWlpWny5MlGlAkAAAAPkVdW7dqOjzrwYsr9pcY2jv7ZW1DRITUBAADAd2zKblyxNDjp4F0vu8c4Ox85HFJOSfVBnwcAAADvZfgKdbvdrrlz5+qKK65QQEBjOdHR0brmmms0e/ZsxcbGKioqSn/84x81efJkTZo0ycCKAQAAYLS8shrXdnzk4du3uwfq6YUtn48EAAAA/7S5xYF6qGs7s7hKqXFhB30uAAAAvJPhgfoPP/ygtLQ0XX311Qc89vjjj8tsNuvcc89VTU2NTjrpJD333HMGVAkAAABPklfqtkI9snUr1NMI1AEAAHAYLQ3Uk90C9aziqg6tCQAAAMYwPFCfMWOGHA5Hs4+FhITo2Wef1bPPPtvJVQEAAMCTNVmh3oKW7z26EKgDAACg5TbnOAP1sCCLesYefNU5gToAAIDv84gZ6gAAAEBrtLble2iQxbWSPa2QE50AAAA4uLJqq9LrjxkHJkbKbDYd9Ln7t3wHAACA7yFQBwAAgNfJbWXLd6mx7fu+8hpV1tZ1SF0AAADwfltyylzbh2r3LhGoAwAA+AMCdQAAAHid/PoV6kEWs2LCAlv0Gvc56umsUgcAAMBBbMpq2fx0SYoJC1RooEUSLd8BAAB8FYE6AAAAvE5Dy/dukcEymQ7egtNdSixz1AEAAHB4m7MbA/UhhwnUTSaTkmOcI4gyi6vkcDg6tDYAAAB0PgJ1AAAAeJXaOrsKK2olOQP1luoZR6AOAACAw2sI1E0maVBi5GGf372L8ziz2mpXUaW1Q2sDAABA5yNQBwAAgFfZV17j2k6Ianmg3rTlO4E6AAAADmSzO7Q11zlDvWdsmMKDAw77mu71K9Ql2r4DAAD4IgJ1AAAAeJWGdu+SFB8ZcohnNpVKy3cAAAAcxu59Faq22iUdfn56g+ToUNd2RhGBOgAAgK8hUAcAAIBXySutdm3Ht6Lle7fIYAUHOA9/CdQBAADQHPf56S0N1Lt3aQzUWaEOAADgewjUAQAA4FWarFBvRct3k8nkWqWeXlgpu93R7rUBAADAux1JoJ4cQ6AOAADgywjUAQAA4FWOtOW71Nj2vabOrny3WewAAACAtH+gHtmi13R3C9QzCdQBAAB8DoE6AAAAvIp7y/durWj5LkkpbnPU9xbQ9h0AAABNbc4ukyRFhQQ0CcoPJTE6RCaTc5sV6gAAAL6HQB0AAABe5UhbvkuNK9Ql5qgDAACgqaKKWuXUX7w5KClKpoaU/DACLWYl1HdOyiyuPsyzAQAA4G0I1AEAAOBV8sqcJynNJikunEAdAAAA7cO93fuQFs5Pb5Ac4wzU95XXqNpqa9e6AAAAYCwCdQAAAHiVvFLnCvWuEcGymFu2aqhBalxjoJ5OoA4AAAA3m9oQqHfv0nicmV3CKnUAAABfQqAOAAAAr2GzO7Sv3BmoJ0SFtPr1KV1YoQ4AAIDmuQfqg49whbrEHHUAAABfQ6AOAAAAr1FQUSO7w7kdH9m6du+SFBpkUbf61xGoAwAAwN3m7DJJksVsUv+EiFa9tntMqGs7s4hAHQAAwJcQqAMAAMBrNLR7l6T4qNYH6lLjHPX8shpV1TLfEgAAAFJtnV078pyBep+u4QoJtLTq9cnRboE6K9QBAAB8CoE6AAAAvEZ+WWOg3i2y9S3fpcZAXZLSi1ilDgAAAGlnfrmsNmcrpNa2e5ek7l0aA3VavgMAAPgWAnUAAAB4jdzSatf2kbR8l6QUt0A9rYBAHQAAANKGzBLX9pDk1gfqyW4t37NKCNQBAAB8CYE6AAAAvEae2wr1Iw3Ue8U1Buo78svbXBMAAAC8n3ugPqJ7dKtfHxUSoIjgAEnMUAcAAPA1BOoAAADwGnllbivUo46s5bv7iqP1bidOAQAA4L/WuR0XDj2CQN1kMql7/Sr1rJJq2e2OdqsNAAAAxiJQBwAAgNfIK237CvV+3SIUEug8DF6fQaAOAADg76w2uzZllUpydjOKDg08ovdJjnFe8FlbZ1dBRW271QcAAABjEagDAADAa7i3fO8acWSBeoDFrCFJzlXqaYWVKqm0tkttAAAA8E7bc8tVU2eXJA3vEXPE7+M+Rz2zmLbvAAAAvoJAHQAAAF4jvz5QjwsPUlDAkR/KDndr47khi1XqAAAA/qyt89MbdO/SGKhnEagDAAD4DAJ1AAAAeAWHw+EK1LsdYbv3Bu4rj9bR9h0AAMCvrcssdm0P79GGQD2GQB0AAMAXEagDAADAKxRXWlVrc7bijI8KadN7NVmhnkmgDgAA4M/Wu11gOTQ56ojfx73le0YRgToAAICvIFAHAACAV3Cfnx7fxhXqfbuFKyTQeSi8nkAdAADAb9XW2bU5u0yS1KdbuCJDAo/4vVihDgAA4JsI1AEAAOAVckurXdttDdQDLGYNTXauUk8rrFRxZW2b3g8AAADeaVtumasLUlvmp0vOY1SL2SRJyiohUAcAAPAVBOoAAADwCu25Ql3av+17aZvfDwAAAN7HvVvR8B4xbXqvAItZifWjiTJp+Q4AAOAzCNQBAADgFfLK3Faot3GGuiQNcwvU12UWt/n9AAAA4H3Wuc1PH97GFeqSlBzjPE4tqrSqsrauze8HAAAA4xGoAwAAwCvklbbvCvURPdxXqDNHHQAAwB+tr7+w0mSShiZHtfn9ms5Rrz7EMwEAAOAtCNQBAADgFfKbtHxv+wr1vt0iFBpokdS01ScAAAD8Q02dTVtzyiRJ/bpFKDw4oM3vmewWqGcW0/YdAADAFxCoAwAAwCs0bfne9hXqFrNJQ+pXIaUXVqmoorbN7wkAAADvsTWnTFabQ5I0vEfb271LTQP1LAJ1AAAAn0CgDgAAAK+QV79CPSokQCH1K8vbyn1O5oYsVqkDAAD4E/f56SPaYX66JHXvQqAOAADgawjUAQAA4PEcDodrhnp8VNvbvTdwD9Rp+w4AAOBf1rsF6u21Qt19hnpmEYE6AACALyBQBwAAgMcrq6lTldUmSYqPbHu79wYj3E6cup9QBQAAgO9bV39BpdkkDUlq/5bvzFAHAADwDQTqAAAA8HgNq9Ol9g3U+3SLUFiQs308K9QBAAD8R7XVpu25ZZKkAQmRCg1qn5FCEcEBig4NlCRllRCoAwAA+AICdQAAAHi8vLJq13Z7tny3mE0akhQlScooqlJRRW27vTcAAAA81+bsUtXZHZKajgFqDw2r1LOLq2Wr/wwAAAB4LwJ1AAAAeLz8so5ZoS41nZfJKnUAAAD/4H7c117z0xt0j3FeAFpndzQ5jgUAAIB3IlAHAACAx3Nv+d6tvQP17gTqAAAA/mZdhlug3s4r1LszRx0AAMCnEKgDAADA4zVp+R7Zfi3fpf0C9QwCdQAAAH+wof5CygCzSYPrRwC1l2S3QD2LQB0AAMDrEagDAADA4+W5tcpMiGrfFep9ukUoLMgiiRXqAAAA/qCq1qZtuWWSpAEJkQoJtLTr+yezQh0AAMCnEKgDAADA47m3fI+Pat8V6hazSUOTnauSMourVFhR267vDwAAAM+yKbtEdodzu73bvUtS9y6sUAcAAPAlBOoAAADweA0t38OCLIoIDmj39x/GHHUAAAC/0WR+eo8OCNRp+Q4AAOBTCNQBAADg8RpWqMdHtm+79wYj3E6kbiBQBwAA8GnuF1CO6IBAvVtEsAItJklSRhGBOgAAgLcjUAcAAIBHK6m0qqymTpIUH9m+7d4buLf6XJ9BoA4AAODLGo73Ai0mDUyMbPf3N5tNSox2HreyQh0AAMD7EagDAADAo/26c59re1gHzLiUpN5dIxQeZJFEy3cAAABfVlFTpx355ZKkQYlRCg6wdMjnNLR9L62uU1m1tUM+AwAAAJ2DQB0AAAAe7edt+a7towd07ZDPsJhNGprsDOszi6tUWFHbIZ8DAAAAY23MKpXD4dzuqIs1JSnZbY76nn2VHfY5AAAA6HgE6gAAAPBYDofDFagHBZg1sXdch32W+wlVVqkDAAD4po6en+56b7djy8U79h3imQAAAPB0BOoAAADwWDvyypVVUi1Jmtg7VqFBHdOSU2p6QnV9RnGHfQ4AAACM436cN7wDV6gfOzDetb1wa16HfQ4AAAA6HoE6AAAAPNYi93bv/bt16GexQh0AAMD3ras/zgsKMGtAQmSHfU6vruHqFRcmSVq5t4g56gAAAF6MQB0AAAAe6+ftje0xjxnYsYF6n67hCq9fAb8+g0AdAADA15RVW7Urv0KSNDgxUkEBHXtq9JgBzuPXOrtDv+4o6NDPAgAAQMchUAcAAIBHqrbatGyX88RjYlSI+sdHdOjnmc0mDa1fpZ5VUq2C8poO/TwAAAB0ro1Zpa7t4R04P72Be9v3Rdto+w4AAOCtCNQBAADgkZbvLlRNnV2SdPSArjKZTB3+mcNp+w4AAOCz3LsQjege0+GfN6lPnGsV/MKt+XI4HB3+mQAAAGh/BOoAAADwSD+7z08f0LHt3hs0CdRp+w4AAOBT1rldMNkZK9RDgyya1CdOkpRdUq1tueUd/pkAAABofwTqAAAA8EiL6gN1s0ma2q9rp3ym+4lVVqgDAAD4lvUZxZKk4ABzh48TanCM24WhtH0HAADwTgTqAAAA8DhZxVXanudcwTMyJUYxYUGd8rm948IVERwgSdpAoA4AAOAzSqqs2lNQKUkakhylAEvnnBY9dmBjoL5wa/4hngkAAABPZXignpmZqcsuu0xxcXEKDQ3V8OHD9fvvv7sedzgc+tvf/qakpCSFhoZq+vTp2r59u4EVAwAAoKP9st2t3Xv/zmn3Lklms0lDk6MkSVkl1dpXXtNpnw0AAICOszHTfX56x7d7b9Cna7hSYkMlSSv2FKq8pq7TPhsAAADtw9BAvaioSFOmTFFgYKC+/vprbdq0Sf/5z3/UpUsX13Mee+wxPfXUU3rhhRe0bNkyhYeH66STTlJ1dbWBlQMAAKAj/bxtn2u7s+anN2gyR51V6gAAAD6h6fz0mE77XJPJpGMHxEuSrDaHluws6LTPBgAAQPswNFB/9NFHlZKSorlz52rChAnq3bu3ZsyYob59+0pyrk5/4okndM899+jMM8/UiBEj9PrrrysrK0uffvqpkaUDAACgg9jsDi3e4QzUo0ICNLJH560gkprOUd+QQaAOAADgC9a7HdeN6OTjS/c56gu3MkcdAADA2wQY+eHz58/XSSedpPPPP1+LFi1S9+7dddNNN+m6666TJO3evVs5OTmaPn266zXR0dGaOHGilixZoosuuuiA96ypqVFNTWNrztLSUkmS1WqV1Wrt4D8RAKCjNXyX850O+K7V6cUqqXL+jB/VN04Ou01Wu63TPn9wQrhre216Ed83XoR9BADgYNhHYG1GsSQpNNCs1JjgTv23ML5nlAItJlltDi3cmqfa2lqZTKZO+3wAh8d+AgD8T2u+8w0N1Hft2qXnn39es2fP1l/+8hetWLFCt956q4KCgnTFFVcoJydHkpSQkNDkdQkJCa7H9vfII4/o/vvvP+D+7777TmFhYe3/hwAAGOL77783ugQAHeTrdJMkiyQppipLX32V2amfb3dIwRaLamwm/b4rT1999VWnfj7ajn0EAOBg2Ef4pwqrlFHkPA2aGGLTt9983ek19I4wa1uJWZnF1Zr70ddK5DQl4JHYTwCA/6isrGzxcw0N1O12u8aNG6eHH35YkjR69Ght2LBBL7zwgq644oojes85c+Zo9uzZrtulpaVKSUnRjBkzFBUV1S51AwCMY7Va9f333+vEE09UYGCg0eUA6ABzX1omydmS86ZzjlNSdEin1/BOzgot31Ok4lqTJhx9grpGBHd6DWg99hEAgINhH+HfFu8okH5fKUk6elhPnXLKoE6vITt6j/75zTZJkil5qE45qmen1wDg4NhPAID/aehy3hKGBupJSUkaMmRIk/sGDx6sjz76SJKUmJgoScrNzVVSUpLrObm5uRo1alSz7xkcHKzg4ANPeAYGBrIjBAAfwvc64JtKKq1aVz/fsn98hFK7RhpSx8iUGC3fUyRJ2pJbqeO6RBhSB44M+wgAwMGwj/BPm3LKXdujUrsY8m/ghMGJrkD9lx0Fuv6Yfp1eA4DDYz8BAP6jNd/35g6s47CmTJmirVu3Nrlv27Zt6tnTeYVm7969lZiYqAULFrgeLy0t1bJlyzR58uROrRUAAAAdb/GOfbI7nNtHD+hmWB3Duke7ttdnlhhWBwAAANpufUbj8dxwt+O8ztQvPkLdY0IlSct2Faqyts6QOgAAANB6hgbqt99+u5YuXaqHH35YO3bs0Ntvv62XXnpJN998syTJZDJp1qxZ+sc//qH58+dr/fr1+sMf/qDk5GSdddZZRpYOAACADvDztnzXtpGBuvuJ1nUZBOoAAADerOECyfAgi3p3NabzkMlk0jEDnce3tTa7lu4qMKQOAAAAtJ6hgfr48eP1ySef6J133tGwYcP04IMP6oknntCll17qes5dd92lP/7xj7r++us1fvx4lZeX65tvvlFISOfP0gQAAEDHcTgc+nm7M1APDjBrYu9Yw2rpFReuyGDndKQNrFAHAADwWpnFVcosrpIkDe0eLYvZZFgtx7hdMLpwa/4hngkAAABPYugMdUk67bTTdNpppx30cZPJpAceeEAPPPBAJ1YFAACAzrYjr1zZJdWSpAm9YxUSaDGsFrPZpKHdo7R0V6FySquVV1at+Egu6AQAAPA2n67OdG1P7dfVwEqkKf26KtBiktXm0MKt+XI4HDKZjAv4AQAA0DKGrlAHAAAAGixya/d+jIHt3huM7BHj2mYFEQAAgPdxOBz6cGWG6/bZo7sbWI0UERygcT2dXZjSCiu1p6DS0HoAAADQMgTqAAAA8AieFqifNCzRtf2R24lYAAAAeIdVacXava9CkjS5T5xSYsMMrkg6dqB72/c8AysBAABASxGoAwAAwHDVVpuW7y6UJCVFh6hffITBFUmjU2LUp2u4JGnZ7kKlF7KCCAAAwJu4r04/d2wPAytpdMxA5qgDAAB4GwJ1AAAAGG7Z7kLV1NklSUf37+YRsyRNJlOTE68fr8o8xLMBAADgSaqtNn2xLkuSFBZk0Uy37kNGGpgQqcSoEEnS0l0FqrbaDK4IAAAAh0OgDgAAAMP97Nbu/WgPaPfe4OzR3dWQ7X+8OkMOh8PYggAAANAi323KVVl1nSRp5rAkhQcHGFyRk8lkcrV9r6mza+muAoMrAgAAwOEQqAMAAMBwDYG62SRN7dfV4GoaJceE6qi+cZKkvQWV+n1vkcEVAQAAoCU+cmv3fp6HtHtvcCxt3wEAALyK3wTqW3NKjS4BAAAAzcgqrtL2vHJJ0qiUGEWHBRpcUVPnjmk8Aet+YhYAAACeKbe0Wr9sdwbV3WNCNbF3rMEVNXVUv64KMDvbIC3aRqAOAADg6fwmUH9vRbrRJQAAAKAZntruvcHJwxIVHmSRJH25Lps5lwAAAB7uk9WZstdP6jl3bA+Z68NrTxEVEqgxPbtIknbvq9DeggqDKwIAAMCh+E2g/sW6LJXX1BldBgAAAPbz83bPDtTDggI0c3iSJKmspk7fbswxuCIAAAAcjMPh0IduXYXOHdPdwGoOzr3tO6vUAQAAPJvfBOqVtXZ9ujrT6DIAAADgps5m1+Lt+yRJ0aGBGtkjxtiCDsK97fvHqzimBAAA8FTrMkq0o36c0IReseoZF25wRc07ZgBz1AEAALyF3wTqkvTm0r1yOBxGlwEAAIB6azNKVFrt7CI0tV9XWTysHWeDib1j1T0mVJL0y/Z85ZZWG1wRALSPqlqbvtuYo6KKWqNLAYB20WR1+ljPXJ0uSUOSohQfGSxJWrKzgLFCAAAAHsyvAvUtOWValVZkdBkAAACo13R+elcDKzk0s9nkahdqd4jORwB8xl8/Wa/r31ipa15bYXQpANBmNXU2zV+bJUkKCTTrlPqxPZ7IZDK5VqlXWW1asafQ4IoAAABwMH4VqEvSm0vTjC4BAAAA9X7amufa9sT56e7OcWv7/tGqDDofAfAJq9OLJUmr0opVU8fqSADebcHmPJVUWSVJJw9NVGRIoMEVHdqxA+Nd27R9BwAA8Fx+E6hHhwZIkr5cl61CWtkBAAAYbnVakdZllEiSBidFKSk61OCKDq1X13CN69lFkrQtt1wbMksNrggA2q68ps61nVdaY2AlANB2HzVp997jEM/0DFP7dVXDxKOftuZxwSYAAICH8ptA/azRzhadtTa7Pvg93eBqAAAA8PIvu13bVx7V08BKWs79xOxHqzIO8UwA8A4VboF6dkm1gZUAQNvkl9VoYf04oaToEB3V13PHCTWIDgvU2PoLNnflV2jlXkZVAgAAeCK/CdTPH5fi2n57eZrsdq74BAAAMEp6YaW+3pAtSeoaEaQzR3U3uKKWOXVEkoIDnIfQn63JVG2d3eCKAODI2e0OVdY2tnnPLqkysBoAaJvP1mTKVn++7+zR3WVpWPrt4S4an+ra/t8vuwysBAAAAAfjN4F6r7hwTe3nvDJ1b0GlFu/YZ3BFAAAA/uvVX3er4frGP0zupZBAi7EFtVBUSKBmDE2UJBVVWpvMgAcAb1NRW9fkdg4r1AF4KYfDoQ+9rN17g9NHJishKliS9N2mXO3ZV2FwRQAAANif3wTqknTZpMYrPt9cutfASgAAAPxXSZVV769wjuAJCTTrskne0e69wbljGlfTu8/pBABvU1Fja3Kblu8AvNXGrFJtySmTJI1OjVHfbhEGV9RyQQFmXXFUL0mSwyG9snj3oV8AAACATudXgfr0wQmuKz5/2JxLOzsAAAADvLM8TRX1LYbPHdNDseFBBlfUOtP6d1N8pPOY8qeteSqsqDW4IgA4MuU1rFAH4BvcV6ef50Wr0xtcOqGnwoKcHZs+WJmuIo4vAQAAPIpfBeoBFrNrLpHdIb2zPN3gigAAAPxLbZ1d837dI0kymaRrpvY2tqAjYDGbdPZo5yp1q82h+WsyDa4IAI5Mxf6BeimBOgDvU1tn1/y1WZKcq71PG5FscEWtFx0WqAvGpUiSqq12vbWMzpoAAACexK8CdUm6aEKKLGaTJOnd5Wmy2uwGVwQAAOA/vlyf5Qpspg9OUB8vasfpzn0u50erCNQBeKcDAnVWqAPwQu4dg2YMSVB0aKDBFR2Za6b2Vv0pS837ba9q6myHfgEAAAA6jd8F6knRoTphULwkKa+sRgs25xpcEQAAgH9wOBz638+NMyGvm9bHwGraZkBCpIZ3j5Ykrc8s0bbcMoMrAoDWK9svUM8rq1YdF50D8DIfubV7P9cL2703SIkN08nDEiVJ+8pr9NnqLIMrAgAAQAO/C9Ql6bJJPV3bby5NM7ASAAAA/7FkZ4E2ZZdKkkb2iNb4Xl0Mrqhtzh3T3bXtfiIXALzF/ivU7Q4pv7zGoGoAoPUKymv045Y8SVJ8ZLCm9etqcEVtc63bBacvL94lh8NhYDUAAABo4JeB+tR+XdUzLkyStHjHPu3eV2FwRQAAAL7vf7/scm1fO62PTCaTgdW03RmjuivQ4vwzfLI6k1WdALzO/oG6JGXT9h2AF5m/Nkt1dmfofPbo7gqwePepzjGpXTSup/Oi02255Vq0Ld/gigAAACD5aaBuNpt06cRU1+23lu41sBoAAADftyOvTD9tdZ4Q7B4Tqpn17Sy9WWx4kI4b2DhKaPGOfQZXBACtU15z4Hxe5qgD8BZ2u0PvrUh33fbmdu/umqxS/2X3IZ4JAACAzuKXgboknTc2RUEBzj/+ByszVG098EQCAAAA2of7ycCrpvTy+tVDDdxP3H60KtPASgCg9VihDsCbfbE+W1tyyiRJI1NiNCAh0uCK2seJQxKadNbcmFVicEUAAADwjTOZRyA2PEinDk+SJJVUWfXFumyDKwIAAPBN+WU1+ni1M2yODA7QheNTDK6o/Rw3MF5dwgIlSd9uzFFeKUEUAO9R3kygnlNSZUAlANA6NXU2/evbLa7bs08cYGA17ctiNumaqb1dt19hlToAAIDh/DZQl6TLJrm1fV9G23cAAICO8MbSvaqtc84Xv3hiqiJDAg2uqP0EBZh1Xv0q9do6u55csN3gigCg5VihDsBbvbU0TemFzguApvSL09H9uxpcUfs6b2wPRYc6j5nnr81iHAcAAIDB/DpQH5PaRYMSne2gVqcV00IJAACgnVXV2vTmUueFiwFmk648qpexBXWAG4/pq4jgAEnSuyvStTO/3OCKAKBlKmqbW6FOaAPAs5VWW/X0j40XMc6ZOVgmk8nAitpfWFCALp/UU5JUZ3do3m97jC0IAADAz/l1oG4ymXRZ/cGpJL25NM3AagAAAHzPR6syVFhRK0k6dUSSkmNCDa6o/cVFBOuGo/tIkmx2h/797VaDKwKAlimvsbm2gyzO0wM5jK4A4OFeWLhTRZVWSdKZo5I1rHu0wRV1jD8c1dP13fz2sr3NjukAAABA5/DrQF2SzhrdXeFBFknSZ2syVVBeY3BFAAAAvsFud+jVxY0zH6+b1sfAajrWNdN6q2tEsCTp6w05WpVWZHBFAHB47i3fe3UNkyTlllbLbncYVRIAHFJ2SZVeqT++DLKYdeeMgQZX1HHiI0N05qhkSVJpdZ3eX5FucEUAAAD+y+8D9YjgAJ1bP/eystam/3y/zeCKAAAAfMOCLXnata9CkjSpT6zPrh6SnG05Z03v77r9z6+3yOEgkALg2cqrnYF6SKBZPbo4A3WrzaGC+s4iAOBpHv9+m2rq7JKkP0zuqZTYMIMr6ljXul2Q+uqvu1VnsxtYDQAAgP/y+0Bdkm45rp9rlfo7y9OYpQ4AANAO/vfLLte2L69Ob3Dh+BT17houSVq+u1ALt+YbXBEAHFpD++CI4AAlRoe47meOOgBPtDWnTB+uzJAkRYYE6Obj+hlcUccbmBipowd0kyRlFFXp2425BlcEAADgnwjUJcVHheiPJzhXFDkc0gOfb2JFEQAAQBusyyjW8t2FkqQ+3cJ13MB4gyvqeIEWs/50UmPb0Ue/2SIbbZMBeLCKWmegHh4coKSoxkA9u6TKqJIA4KAe/WaLGg6tbj6un7qEBxlbUCe53u3C1Jd+2cU5SwAAAAMQqNe7akov9YxztolatrtQX63PMbgiAAAA7/XUgu2u7Wun9pHZbDKwms4zc1iiRqbESJK25JTp09WZxhYEAIfQMEM9PGi/FeqlrFAH4FmW7CzQj1vyJElJ0SG68qhexhbUiab0i9OgxEhJ0tr0Yv2+t8jgigAAAPwPgXq94ACL7jl1iOv2w19tVrXVZmBFAAAA3um7jTn6YbPzhGe3yGCdM6a7wRV1HpPJpLtPHuS6/d/vt3FMCcAj1dTZZLU5VzlGBAcoKTrU9Vg2Ld8BeBC73aFHvt7sun3HjIEKCbQYWFHnMplMTcYnPb9wp4HVAAAA+CcCdTfTB8drWv+ukqTM4iq99POuw7wCAAAA7ipq6vT3+Rtdt+89bYhfnfCUpMl943TsQOesy8ziKr25dK/BFQHAgSpqGi/2CQ+2MEMdgMf6cn221mWUSJIGJUbq7NH+c7Fmg9NHJiuxfjTHj1vytGAzs9QBAAA6E4G6G5PJpL+dNkSW+pakzy3coaxiZscBAAC01OPfb1NWfRAzrX9XnT4iyeCKjHHXSYNkqu9y/8xPO1RSZTW2IADYT0O7d8k5Q909UGeGOgBPUVtn17++3eq6fffMQa7zdv4kKMCsOac0dkG6b/5GVdXSBQkAAKCzEKjvp39CpC6f1FOSVG21659fbzG4IgAAAO+wIbNEr/66W5IUHGDWP84aJpPJ/054StKQ5CidPcq5eqq40qoXF9GaE4BnKXcL1COCAxQRHKDI4ABJUm5pjVFlAUATby3bq7TCSknOWeLHDOhmcEXGOWNksqb0i5MkZRRV6ekftxtcEQC03Psr0nXl3OXakFlidCkAcEQI1Jtx+/QB6hIWKEmavzZLK/YUGlwRAACAZ7PZHfrrJ+tld47j1R+P76eeceHGFmWw208coCCL83D71V9300IZgEfZf4W6JNcq9eySKjkcDkPqAoAGpdVWPf3jDtftu08e7LcXa0rOzpoPnDnMdXz50s+7tC23zOCqAODwKmvrdM9nG7Rwa74e/YYFjAC8E4F6M6LDAnXHjIGu2/d/vlE2OycTAAAADuatZXu1tn62Zb/4CF1/dF+DKzJeSmyYLp/c2PnoyQXbDK4IABrtv0JdagzUq612RlUAMNyLi3aqsKJWknTmqGQN7xFtcEXG69stQjce6zzOrrM7dM+nG7gACoDHSy+sUm2dXZK0Nr2Y7y0AXolA/SAunpCqQYmRkqQNmaX6cGW6wRUBAAB4ptzSav3rm8bZlg+fPVxBARxmStLNx/VztVB+b0W6duSVG1wRADg1F6gnNZmjTlcNAMbJKanWK4udo4SCLGbd6bbwxd/ddGxf9YwLkyQt312oj1ZlGlwRABxaw+gOSSqtrlNGUZWB1QDAkeFM50FYzCbdd/pQ1+1/fbtVpdVcoQ8AALC/B77YpLL6YOaCcT00oXeswRV5jtjwINcqIrtDeoz2dgA8RPMt30Nd9zGmAoCR/vn1ZlVbnasZL5/cUymxYQZX5DlCAi164MxhrtsPf7VZRfUr+QHAE+0tqGhye2NWqUGVAEBTVbW2Fj+XQP0QJveN0ynDEyVJ+8pr9Yzb3CYAAABIP23N05frsiU5w+M5MwcbXJHnuWpKL8VHBkuSvtuUq5V7Cw2uCACk8prGEwfhwRZJrFAH4Bk+W5OpT9dkSZIiQwJ0y3H9DK7I8xwzoJtOHZEkSSqsqNVj33LRJgDPle62Ql2SNmUTqAPwDP/4clOLn0ugfhhzZg52tSyd++tu7cqnTScAAIDkvIrz3k83uG7/5ZTB6hIeZGBFniksKECzpg9w3X7wi82qs9kNrAgAmq5Q33+GuiTllNCKE0DnSy+s1D2fNB5fPnjmMI4vD+Jvpw1xfX+/szydizYBeKy0/QN1VqgD8AAf/J6uz+ov4mwJAvXDSIkN0w1H95EkWW0O/ePLzQZXBAAA4Bme+nG7a/bZpD6xOndMd4Mr8lwXjOuhPt3CJUlr0ov13MKdBlcEwN811/KdFeoAjFRns+u2d1e7RgmdM7q7zhrN8eXBJESF6M4ZjRdt/vWTDbJy0SYAD3RgoF5iUCUA4LQ1p0z3frbh8E90Q6DeAv93bF8lRjlPLPy4JU8/bc0zuCIAAABjbc0p0/9+3iVJCrKY9dDZw2UymQyuynMFWMz613kjZK7/K3pywXatTisytigAfq28mRXqSVFuM9RLCdQBdK6nFmzXqrRiSVJqbJjuP3OosQV5gcsn99Kw7lGSpC05ZZr36x5jCwKA/djtDqUXNe18lFVSraKKWoMqAuDvKmrqdNNbK1Vtbd2FiATqLRAWFKA5pwxy3f77/I1NTj4AAAD4E7vdob9+sl51dock6cZj+6pvtwiDq/J8Y3vG6pbj+0uSbHaHbn9vTZMVogDQmZpboR4VGqCQQOdpAlaoA+hMy3YV6JmfdkiSAswmPXXxaEWGBBpcleezmE166Kzhariu9fEftimzmJEdADxHXlmNausODK2Yow7ACA6HQ/d8ukE78yskSQMSWn4+k0C9hc4YmaxxPbtIkvYWVGrOx+vlcDgMrgoAAKDzvf97un7f61xd3btruG46tq/BFXmPW4/vp1EpMZKkPQWVeuDzTcYWBMBvldfYXNvhwRZJkslkUlK0c5V6LoE6gE5SXFmrWe+tUf21mrr9xAGu4yUc3siUGF02sackqbLWpvvnbzS4IgBotLegwrUdE9Z4oRRz1AEY4b0V6fpkdaYkKTzIov9cMLLFryVQbyGTyaT/XDBSkfVX7n++NktvL08zuCoAAIDOta+8Ro98vcV1+8Ezhykk0GJgRd4lwGLWExeOUliQ8+/svd/T9c2GbIOrAuCPKppp+S7JNe6srKZOZdXWTq8LgH9xOBy6+6P1rq4Yk/vE6cZjuFizte48aaC6RgRLkr7blKsfNuUaXBEAOLnPTz9xcIJreyNz1AF0ss3ZpbrP7cLDf547Qr27skK9Q/SMC9ej541w3b7/80188QMAAL9hszs0+/21KqlyBixnj+6uqf27GlyV9+nVNVx/P6NxJujdH69XDitBAXSyilpnoG42SaFuF0YlRYe4tnOZow6gg727Il3fbMyR5Fy5+N8LR8piNhlclfeJDg3UvacNdt2+b/5GVdYyWgiA8dLdAvUTBscr0OL8jqflO4DOVFZt1U1vrVJN/QiKyyf11Okjk1v1HgTqrXTK8CRdMdnZRqm2zq5b3l7NVfsAAMAv/Oe7rfp5W74kKS48SH89dfBhXoGDOX9sD80clihJKq606s4P1spuZ5wQgM5TXu0MWsKDA2QyNYZXiW6BOnPUAXSkHXnluv/zxlVCj547wjV2Aq13xshkTekXJ0nKLK7SY99sNbgiAGi6Qr1vtwj1j4+UJO3Mr1C11XawlwFAu3E4HJrz8Xrt3uccQTGse5TuOa315zQJ1I/AX04drOHdoyVJu/dVME8dAAD4vK/WZ+u5hTslSRazSc9cMsbVVhKtZzKZ9PDZw5UQ5fw7XLxjn179dbfBVQHwJ+X1Ld/d271LTVeoE6gD6Cg1dTbd+s5qVVudq4QunZiqk4YmGlyVdzOZTHrwzGEKsjhP9877bY8+X5tlcFUA/J17oJ4SG6YhyVGSnB3wtuaUGVUWAD/y5rI0fbHOOW4xMjhAz14yRsEBrR9fSaB+BIIDLHr2kjGueepfrMvWW8uYpw4AAHzTttwy3fnBWtftv54yWJP7xhlYkW/oEh6k/5w/ynX7sW+2ajNt7wB0koYZ6uH7BeqJbqtDGUcBoKM89s1WV7vf/vERuufUIQZX5Bv6dIvQvac3/l3++aN12pZLYAXAOGmFVZKkhKhghQRaNLQ+UJdo+w6g423ILNGDn29y3X7svBHqGRd+RO9FoH6EUuPC9JjbPPUHvtikDZnMUwcAAL6lpMqq61//XZW1zlZsZ4/urqum9DK2KB8ytX9XXTettySp1mbXrHfX0PYOQIez2x2qqP9e3z9QZ4U6gI62cGueXlns7MwTFGDWUxePVmhQ61cJoXmXTUzVOWO6S5Iqa2268Y2VjKsEYIiKmjrtK6+RJKXGhkmShiQ1Buobs8hTAHSc0vq56bU2Z0ekK4/qpZnDk474/QjU22Dm8CRdeVQvSQ3z1FdxgAoAAHyG3e7QrHdXa0+Bs0Xb0OQoPXz28CazdtF2d540UIPrTypszS3TP7/eYnBFAHxdpduFOxHBTUMs9xnqOSVVnVYTAP+QX1bTpPPRnJmDXMdBaB8No4UaQqtd+yp05wdrGVcJoNOlFzVt9y5Jg91XqGexQh1Ax3A4HPrzh+tcYydG9ojWX05p/dx0dwTqbTTnlEEa0cM5T31PQaXuZp46AADwEU/8sE0/bc2XJHUJC9QLl41l9VAHCA6w6MmLRik4oHHe5aJt+QZXBcCXNbR7l6TwoKYr1GPDghRocV44xQp1AO2psrZO177+u/aV10qSjhvYzbVQBe0rJNCiFy4bq6gQ53f8txtz9cKiXQZXBcDfpBU0BuoNK9SjQgJd25uzy2Szk6UAaH+vLN6trzfkSJKiQgL0zCVjFBTQtkicQL2NggMseubiMYqsP0D9cl223mSeOgAA8HLfbszRUz/ukCSZTdIzl4xxXVGO9jcgIbLJlbJ3frBWBfWt8QCgvZW7BeoR+7V8N5tNSohyrlLPKSVQB9A+rDa7bnprldamF0tyztL91/kj6XzUgVLjwvTkRaPV8Ff8r2+36Ncd+4wtCoBfaVgZKjUG6lJj2/cqq017Cio6vS4Avu3T1Zl66KvNrtv/Pn9ku5zTNDRQ//vf/y6TydTkv0GDBrker66u1s0336y4uDhFRETo3HPPVW5uroEVNy81Lkz/cpun/uDnzFMHAADea0deme54370V52BN6dfVwIr8wx8m99SxA7tJcrZDve3dNbLWz3kCgPbUZIX6foG61DhHvbjSqmq39vAAcCQcDofu/mi9FtZ3PooMCdC8qyaoa0SwwZX5vuMGxeu2E/pLkuwO6Y/vrFZWMeM8AHSOdLdAvWdcY5g1NNl9jjpt3wG0nx825eqOD9aqoZH4Lcf104yhie3y3oavUB86dKiys7Nd/y1evNj12O23367PP/9cH3zwgRYtWqSsrCydc845BlZ7cCcPc5unbrPrZuapAwAAL1RabdX1b6x0rV48fWSyrp3W2+Cq/IPJZNJj541QXHiQJGnxjn26+yPGCQFof01WqIccGKgnRoe6tnNo+w6gjR77dqs+WpUhSQqymPW/P4xjbnonuvX4/jqu/qLNwopa/d9bq1RTx8VSADre3sIDZ6hL0hDmqAPoAL/t3Keb3l7lGiVx6cRU3TFjQLu9/4G/OXeygIAAJSYeeHVASUmJXnnlFb399ts6/vjjJUlz587V4MGDtXTpUk2aNKnZ96upqVFNTWN7zNJS5xey1WqV1dqxAfefTuynVXsLtS6zVHsLKnXrO6v03MWjFGAx/LoFAPAZDd/lHf2dDvgju92h299do135zpZrgxIi9I8zBqmuru4wr0R76RJi0bMXj9Qf5q1UbZ1dH63KUFJUkG47oZ/RpXkF9hFAy5RWNP7OHBpgOuBnJj4i0LWdUViu7tFBnVYb0FHYRxjjtSV79fzCnZIkk0n693nDNDYliv8Pneyxc4bp7BeWKqOoSmvTi3XfZxv04BlDjC4L8CjsJ9pfWn0795BAs2KCza6/2wHxjeH6hsxi/s4BtNnajBJd99rvqq1zdno8bXii7j1l4GHPabbm+8fkMHDJy9///nf961//UnR0tEJCQjR58mQ98sgjSk1N1Y8//qgTTjhBRUVFiomJcb2mZ8+emjVrlm6//faDvuf9999/wP1vv/22wsI6fu5nQbX0r3UWVdmcA4rGd7Prkr52mRkJBQAAPNw36SZ9nWGRJIVZHLpjhE1dQwwuyk+tKTBp3jazHHIeRF7Ux6bJCaxUB9A+VuSb9OYO5/f9ub1sOjqp6ffLwmyTPtnjfPyyfjaN78b3D4DWW73PpNe2Nx7PnNfbpmmJfJ8YJaNCemK9RVaH8//HxX1tmhTP/w8AHcPukO5cZpHNYVJiqENzRjV2xnA4pL/+blFFnUkRAQ79Y5xNJvITAEcou1J6aqNFlXXOL5IhMXZdO9Culqx1rqys1CWXXKKSkhJFRR26g5KhK9QnTpyoefPmaeDAgcrOztb999+vadOmacOGDcrJyVFQUFCTMF2SEhISlJOTc9D3nDNnjmbPnu26XVpaqpSUFM2YMeOwfxntpe/IAl37xipZbQ6tyDdraL9e+svMgTKxVwCANrNarfr+++914oknKjAw8PAvANAin63N1tdL1kuSzCbpmcvGahpz0w1ziqSk3/bq4a+3SpI+2BOgE6eO1tH9+X9yKOwjgJYpWp4u7dgsSRo/eoROGdO9yePmjbn6ZM9aSVJi70E65WhGf8D7sY/oXEt2Feit11fJIWdg+3/H9Nbs6f0Nrgrd+mXpro83SJI+2huoC2dMaDLLGPBn7CfaV3ZJtWxLf5YkDUmN1ymnjG7y+Af5K/XrzgKV15k0/ugTFB8ZbESZALxcWmGlHnp5hSrrnF3YJvTqolf+MEYhgZYWvb6hy3lLGBqoz5w507U9YsQITZw4UT179tT777+v0NDQQ7zy4IKDgxUcfOCXb2BgYKftCI8ZlKinLx6tm95aJbtDmrckTXERIfrjCfziAADtpTO/1wFf98W6LN310XrX7TtPGqjjBycZWBEk6fpj+imntFav/rpbNrtDt767Vu/dMFnDukcbXZrHYx8BHFp1XeOKxOiw4AN+XnrEhru288pr+XmCT2Ef0fE2ZpXoprfXympzftdcMK6H7jp5MAtNPMAFE3pqXVap3lyapto6u255d60+v2WquoQz2gNowH6ifWS7hVQ9u4Yf8Hc6rHu0ft1ZIEnallep7rERnVofAO+XW1qtK19bqbwyZ5g+oke0XrlyvCJDWv4d3prve48a7h0TE6MBAwZox44dSkxMVG1trYqLi5s8Jzc3t9mZ657m5GFJeuSc4a7b//l+m95YutfAigAAAA70zYYc3fbuGtnrs5XLJqXq/47pa2xRcLnn1MGaOcx57FtRa9NV81Yoo6jS4KoAeLuKmsY5cuHBB15nnxTdeIF7dkl1p9QEwDekF1bqyrkrVF7/PXPCoHg9fPZwwnQPcu9pQzQqJUaSlFFUpSvnNf7/AoD2klbY+Htrz9gDR/EOceuOsSm75StEAUCSiipqddnLy5ReWCVJ6h8foXlXTWhVmN5aHhWol5eXa+fOnUpKStLYsWMVGBioBQsWuB7funWr0tLSNHnyZAOrbLkLx6dqzsxBrtt/+2yD5q/NMrAiAACARj9sytUf31klW32aftH4FD1wxjBOeHoQs9mkxy8cpbE9u0iS8stqdOXcFSqptBpcGQBvVn6YQL1rRJDM9buCHAJ1AC1UUF6jP7y6XPn1q4RGp8bomUvGKKAlAyzRaYIDLHr+sjHqVt9eeW16sa6et0JVtbbDvBIAWs49UE+NOzBQdx83sTGrpFNqAuAbyqqtumLucm3PK5ckpcSG6o1rJiq2gzvuGHpEe+edd2rRokXas2ePfvvtN5199tmyWCy6+OKLFR0drWuuuUazZ8/WTz/9pJUrV+qqq67S5MmTNWnSJCPLbpUbjumrG+tXeTkc0uz31mjh1jyDqwIAAP5u4dY83fTWKlcrznPH9NDDZw+X2UyY7mlCAi363x/GqXdXZwvmHXnluv6N31VTx0lPAEfGfYV6RDOBeoDFrPjIEEmsUAfQMuU1dbp63grt3lchSerbLVyvXjFeoUEtm1+JzpUUHao3r5momDDnKq7luwt145srOb4E0G6aBOrNrFDv3TVCIYHOeGpTFivUAbRMtdWm617/XesynBfixEcG661rJikxOqTDP9vQQD0jI0MXX3yxBg4cqAsuuEBxcXFaunSpunXrJkl6/PHHddppp+ncc8/V0UcfrcTERH388cdGlnxE/nzyQF08IUWSVGd36MY3V2rl3kKDqwIAAP7q1x37dP0bK1Vrs0uSzhiZrMfOG0GY7sFiw4M076rxiqu/2nbZ7kLd+cE62e2Ow7wSAA5UUdMYmIQHNx92NZyQKKioUW2dvVPqAuCdCsprdMn/lmpt/YnNhKhgvXb1BOZye7iBiZF6/eoJrgurFm3L123vrFGdje98AG3nHqj36HJgoG4xmzQo0blKfU9Bpcqq6cIG4NBKqqy64tXlWrrLma/GhAXqjWsmNtsFoyMYGqi/++67ysrKUk1NjTIyMvTuu++qb9/GmZ0hISF69tlnVVhYqIqKCn388cdeMT99fyaTSf84a7hOGe6svdpq11VzV2gzs0EAAEAnW7qrQNe8tsIVjswclqj/XjBSFsJ0j9czLlyvXDnedRX/52uz9Oi3WwyuCoA3Kj/MCnVJSqoP1B0OKa+MVeoAmpdZXKXzX1ziWiUUHRqo166e0Gx4As8zokeMXnU7vvxmY47u+pCLNgG0XXp9oJ4QFayQwOYv4HSfo74lp6xT6gLgnbKKq3T+C79p2W5nmB4eZNG8qyZoYGJkp9XQ/G/O+5k9e3ar3/iee+5RbGxsq1/nqyz18y/Lqn/XL9v3qbS6Tn94dbk+uvGoTrt6AgAA+Lff9xTq6nkrVG11huknDknQUxePZq6lFxmVEqNnLh6j69/4XXaH9OKiXYoJDdL/Hdv38C8GgHoVh5mhLqlJy7yckmrCMQAH2JFXpstfWe4aDZEQFaw3rpmoAQmdd2ITbTehd6xeunycrn3td9Xa7Pp4dabCgi168MxhMpm46BZA61XU1Glfea2k5tu9N2gyRz2zRON7kScBONCWnFJd+eoK5ZQ6jzljw4P06pXjNSolplPraFGg/sQTT2jy5MkKCmpZq6bFixfrlltuIVDfT3CARS9cNlaXvrxMa9KLlV9Wo8teWaYPbpyshKiO7+8PAAD815r0Yl05d4Uqa51tfo8b2E3PXDJagYTpXmf6kATdf+Yw3fvpBknSo99sUWVtnWafOICTngBapGGFenCA+aD7gSS3QJ056gD2tya9WFfNXa6iSmeL3t5dw/X61ROUcojgBJ7r6AHd9PQlo3XTW6tkszv05tI0hQcF6O6Zgzi+BNBq6UXu89PDD/q8IUmNgfomuvkCaMZvO/fphtdXqqz+d9iecWF67aoJ6tX14N8tHaVFgbokffLJJ4qPj2/RcyMjuRL1YMKDAzT3yvG64MUl2p5XrrTCSp37/G9645qJ6m3APwAAAOD7NmSW6PJXlrkClGn9u+r5y8YqOKD5tmvwfJdP6qmyaqse+2arJOnpH3eoosame08bzElPAIfVsD84WLt3SUqMDnVt5xCoA3CzePs+Xf/G764LNYcmR+m1qyeoa0SwwZWhLU4amqj/nD9St7+/Rg6H9OLPuxQeHKBbT+hvdGkAvMzeAvdA/eAXWg1KjJLZJNkdBOoADjR/bZbueH+NrDbnKJqRPaL1ypXjDTvmbNGSpLlz5yo6OrrFb/riiy8qISHhiIvydV3Cg/TGNROVEus8QZFRVKXznv9N6zKKjS0MAAD4nOW7C3Xpy8tUVu0MTyb3idNLl4876AwzeI+bju2nv58+xHX71V93a87H62Vj5iWAw2ho+X6wdu8SK9QBNO+r9dm6at5yV5g+sXes3r1+EmG6jzhrdHc9dNZw1+3/fr9NryzebWBFALxRw/x0SUqNCz3o80KDLOrTLUKStC2nXFabvcNrA+D5HA6HXvp5p259Z7UrTD9+ULzeMfiYs0WB+hVXXKHg4JYXeckllyg8nNXWh5IYHaIPbzxKgxKdq/kLKmp10UtL9fO2fIMrAwAAvuLT1Zm67OVlKqlytuKc0CtWr1w5TqFBhOm+4sopvfXYuSNkrl+U/u6KdN3+3hpORAA4pIoaZxB2qEA90W0sWU5pVYfXBMDzvbVsr25+e5XrxOaMIQl67eoJigwJNLgytKdLJqbqnlMHu24/+MUmvbs8zcCKAHibtMKWrVCXGueo19rs2pFX3qF1AfB8NrtD93++SQ9/tcV138UTUvTS5WMVFtTipusdos1DM3ft2qWNGzfKbuekXWslRIXovRsma0Jv56z5ylqbrp63Qp+tyTS4MgAA4M0cDoee/GG7Zr23RrX1werRA7rp1avGG37wifZ3wfgUPXnRaAXUp+rz12bpprdWqdpqM7gyAJ6ots7u2jdEBB/8AquEKFaoA3ByOBx69qcd+usnG+Sob4Rz/tgeeu7SMXQ98lHXTuujWdMbW73P+WS9Xvttj3EFAfAq7oF6ymECdfc56huzaPsO+LNqq023vL1K89yOOW6fPkAPnz1cAZY2x9lt1uIKrFar7rvvPp1++ul66KGHZLPZdPHFF6t///4aMWKEhg0bpj179nRgqb4pOjRQr189QScNdbbIr7M7dNu7a2inBAAAjkhtnV13frBOj/+wzXXfxRNS9coV4w45Kxfe7fSRyXrhsrEKCnAe3n+/KVfXvf67KmvrDK4MgKdpaPcuHXqFelCAWV0jgiQxQx3wZza7Qw9+sVn/+nar674bju6jx84b4REnNtFxbjuhv66b1luS5HBI983fqEe/2SKHg/FCAA6tIVAPDbSo22HaMw9JbgzUNxGoA36roLxGf3hlub7ekCNJsphNeuzcEbpten+ZTCaDq3Nq8ZHv3Xffreeff16JiYl69dVXdc4552j16tV6++239e677yogIEB//etfO7JWnxUSaNFzl47VxRNSXfc9+MUm/fNrDlIBAEDLlVRa9YdXl+mjVRmu++bMHKSHzx6mQE54+rzpQxI098rxCqtv6f/L9n264tXlKq22GlwZAE9S3sJAXXKOKpOkvLIa2ez8bgr4m6KKWl05d7le/bVx0cfdMwdpzimDPebEJjqOyWTSX04ZrJuP6+u67/mFO3XH+2sZLwTgoOx2hzIKneOCUmPDDru/aLpCvaRDawPgmVbsKdSpTy3W8j2FkqSwIItevmKcLhifYnBlTbV4mdKHH36oefPm6ZRTTtG2bds0aNAgffnll5o5c6YkKT4+XpdeemmHFerrLGaTHj57mOIjg/Xkgu2SpBcW7VR+WY3+ee5wToIDAIBDSiuo1FXzlmtnfoUkKTjArMcvHKVThicZXBk605R+XfXGNRN05dwVKquu04o9Rbr0f8v02tUTFBseZHR5ADxAhVvniojDjAFJjArVhsxS2ewO7SuvadIGHoBvW59RohvfXKnMYmco0nDe6sLxqYd5JXyJyWTSn04apMSoEP1t/kY5HNLHqzOVX16j5y8bSwcsAAfIKa12jRc6XLt3SYqLCFZiVIhySqu1KbtUDoeDi7YAP2G3O/TSL7v0r2+3ui7g7hoRrFevHKcRPWKMLa4ZLU5ps7KyNHLkSEnSgAEDFBwcrH79+rkeHzBggHJyctq/Qj9iMpl0+4kD9I+zhqlhn/HRqgzd8MZKVdUyAxMAADRvVVqRzn7uV1eYHhcepHeun0SY7qfG9ozVO9dNUpewQEnS+swSXfTSEtcJcQD+raUt3yUpKZo56oA/em9Fms594TfXsUPXiCC9de1EwnQ/dvnkXnr+0sbxQr9s36eLXlqivDL2DQCacp+fntqCQF2Shta3fS+rrlNGEb+3Av6gqKJW177+u/759RZXmD6xd6y+unWqR4bpUisCdZvNpsDAQNftgIAAWSyWxjcym2lP3k4um9RTz10yRkH1q9J/3JKnS15eqvyyGoMrAwAAnuar9dm6+KWlKqiolST17RauT26aojGpXQyuDEYa1j1a798wWfGRznl123LLdcbTi7V8d6HBlQEwWnlN48XaESEta/kuSTklnNwEfF211aa7P1qnP3+0XrV1ztWFY1Jj9MUfp2lSnziDq4PRTh6WqLeunajoUOf54Q2ZpTr3+d+0K7/c4MoAeJKmgXpoi17jPkd9I3PUAZ+3Kq1Ipz29WD9uyZMkmUzSLcf101vXTlS8B3dFa1Uf8W+//Vbz58/X/PnzZbfbtWDBAtftb7/9tqNq9EszhyfptasnKLJ+xcDqtGKd9vQvWrm3yODKAACAJ7DbHXpu4Q7d9NYq1dSf8JzUJ1Yf/98Upca17Cpw+Lb+CZH64MbJ6ln/76GgolaXvrxUby9LM7gyAEZyX6EeEWw5xDNZoQ74k4yiSl3w4hK9uyLddd8Vk3vq3esnN7m4Bv5tfK9YfXjjZCXX/5tIL6zSeS8s0eo0zlcCcEp3D9RbeG5iqFugvok56oDPcjgcevmXXbrghcYuirHhQZp31QTdedJABXj46OtWDbq54oormty+4YYbmtxmtkX7mtw3Tu/dMFlXzl2uvLIa5ZbW6MIXl+ieUwfriqN68fcNAICfKqyo1R3vr9FPW/Nd950zprv+ec4IVxtGQJJ6xoXrs5un6Ja3V2vxjn2y2hz6yyfrtSm7RPedPlSBHv7LCoD2V17d8pbvTVeoE6gDvuqX7fm69Z3VKqq0SpJCAs165JzhOnt0D4MrgyfqnxCpj2+aoivnLteWnDIVVtTqkv8t07OXjtbxgxKMLg+AwY6k5fuQpGjX9qZsVqgDvqikyqq7Plyrbzfmuu4b36uLnrp4tJKiW9bNwmgtPoNmt9sP+5/Nxpzv9jYkOUpf3DpVE3vHSpLq7A79/fNNuu3dNaqsrTvMqwEAgK9ZuqtAM5/8uUmYfvv0AfrP+SMJ09GsmLAgzbtqvK6Z2tt135tL03Tpy8tUUM5IIcDflDdZoX64GeqNJzZYoQ74HrvdoWd/2qErXl3uCtNTY8P08f9NIUzHISVGh+j9GydrUh/n+coqq03Xvb5S7yynExLg79wD9R5dWhao9+gS6urUu4mW74DPWZfh7MDtHqbfcEwfvX3dJK8J06VWtnyHMeIjQ/TWtRN1/dF9XPfNX5uls579lTlFAAD4CZvdoacWbNcl/1uq3FJnCBoXHqTXrp6g26b3p3MNDinAYta9pw3Rv88fqaD6VenLdxfqjGd+1UZa6gF+xb3le3jQYVaoR7FCHfBV+WU1uu713/Wvb7fK7nDed8KgeH1+y9Qms2yBg4kKCdRrV0/QaSOSJDl/X5nz8Xr99ZP1qqlj0RXgr9IKnIF6YlSIQgIPPV6ogdls0uD6fU9WSbWKKmo7rD4Anae2zq4nftimc5//TemFzhbv0aGBeuWKcZozc7DXdU1sUcv3+fPnt/gNzzjjjCMuBgcXYDHrL6cM1uiUGP3pw3Uqr6nTttxynfHMr/r3+SN08rAko0sEAAAdJK+0WrPeW6Pfdha47juqb5yeuHCU4qOYaYmWO29sD/XtFq4b3lipvLIaZRZX6bznl+jf54/UqSM4ngT8QXlty1u+hwZZFB0aqJIqq7JLqzq6NACd5It1Wbr30w2uVekmkzR7+gDdfFw/mc1cpImWCw6w6KmLRishKkSvLN4tSXprWZo2ZJbo2UvHtHh1KgDfUF5Tp4L6MLyl7d4bDEmK0vLdhZKcbd+n9Ova7vUB6Dxr04v154/WaUtOmeu+0akxeuaSMeoe4z2r0t21KFA/66yzmtw2mUxyOBxNbjeg7XvHmjk8SQMSI3XjGyu1Pa9c5TV1uvHNVbrh6D7600kDFeBlV3QAAIBDW7QtX7PfW+P6pdRskmbVn/C0cMITR2B0ahd9/sepuv6NlVqbXqwqq003v71Km7P7afaJAziRDvi4ila0fJekpOgQlVRZlVtSI7vdwXcE4MUKK2p172cb9OW6bNd9XSOC9O/zR+rYgfEGVgZvZjabdO9pQzQwMVL3frpBNXV2rc0o0WlPL9aTF43WMQO6GV0igE6S7tbuPaWVgfpQt+4oG7NKCNQBL1Vttenx77fpf7/scnVBsphNuvGYPrrthAFePa6yRZW7z0n/7rvvNGrUKH399dcqLi5WcXGxvvrqK40ZM0bffPNNR9cLSX27RejTm6fo9JHJrvte/HmXLntlmfLLmIMJAIAvsNrsevSbLbri1eWuMD0hKlhvXzdJt57QnzAdbZIQFaL3rp+kc8Z0d933zE87dN3rv6uQ9nqAT6uoabwIPjz48G04E6OdnVBqbXYVVvL9AHir7zbmaMbji5qE6acOT9K3s44mTEe7uGBcij6+6SjXqtTiSquunLtcT/6wXXa74zCvBuAL3Oent3qFulugzhx1wDst312omU/+ohd/bgzThyRF6bObp+hPJw3y6jBdauEKdXezZs3SCy+8oKlTp7ruO+mkkxQWFqbrr79emzdvbtcC0bzw4AA9ddEojUmN0UNfblad3aGluwp16lO/6NFzR+i4QfwyBACAt8ooqtRt767Ryr1FrvuOG9hN/z5/pOIigg2sDL4kJNCi/5w/UkOTo/XQl5tkd0gLtuTp5Cd+1n8uGKlp/VlNBPii8iNYod4gp6RaXdkPAV6lpNKq+z/fqI9XZ7ruiwkL1INnDmuyUANoD0OTo/X5LVN1xwdr9MPmPDkc0uM/bNPq9CI9fsEodQkPMrpEAB3IfYV6alzrWjr3j49UoMUkq82hjQTqgFcpr6nTo19v0RtL97ruC7KYddv0/rr+6D5eNyv9YFr9p9i5c6diYmIOuD86Olp79uxph5LQUiaTSVdN6a13r5+khCjnSY28shpdNW+F7v5oncqqrQZXCAAAWsNud+j1JXt00uM/u8L0ALNJfz1lsF65YjxhOtqdyWTSNVN767WrJ6hLWKAk5/Hk5a8s1z++2KSaOsY5Ab6mScv3kMMH6olRjSdDc0qqO6QmAB1j4dY8zXhiUZMwffrgBH13+9GE6egw0WGBeunycfrTSQPV0FRr4dZ8nfb0Yq3LKDa0NgAdqy0r1IMCzOofHylJ2plfrmorv4sC3mDRtnyd9PjPTcL0Makx+uq2qbr5uH4+E6ZLRxCojx8/XrNnz1Zubq7rvtzcXP3pT3/ShAkT2rU4tMy4XrH64o/TdLTbTKJ3V6Tr5Cd+0W879hlYGQAAaKkdeeW64MUl+ttnG1VR6/zFsUeXUH1w42Rdd3QfZtaiQ03r303fzDpa0/o3zql7efFunfnMr9qWW2ZgZQDaW0OgbjZJoYGHb/nuvkI9u5RAHfAGZdVW3f3ROl05d4VyS52jASNDAvTv80fqf38Yq/jIkMO8A9A2ZrNJNx/XT29cM1Fx9avSM4urdN7zS/T2sjQ5HLSAB3zR3gL3QD281a9vmKNud0hbcvg9FPBkuaXVmv3eGl3x6nJlFldJcv5++bfThuiDG49Sv/oLZHxJqwP1V199VdnZ2UpNTVW/fv3Ur18/paamKjMzU6+88kpH1IgW6BYZrNeuGq9Hzhmu8CDnSZHM4ipd8vIy/X3+RlXW1h3mHQAAgBGsNrue+XG7TnnyF/3u1uL94gmp+uq2aRqd2sXA6uBPEqJC9NpVE3TvaUMUVH8F8ZacMp3+9GK99tseTnwCPqKh5Xt4UIBMpsNfrJXYpOV7VYfVBaDtHA6HPl2dqRP+s0jvrkh33X/0gG767vajdd7YHi36uQfay5R+XfXFrVM1OjVGklRrs+svn6zX7e+tUUkVnTUBX9PQ8j000KKuEa0f8cAcdcDzVVtteubH7Tru3wubdEE6qm+cvp11tK6e2lsWH10U1OoZ6v369dO6dev0/fffa8uWLZKkwYMHa/r06RyUG8xkMuniCama2q+r/vThWi3dVShJmvfbHi3cmqf/XDBSY3vGGlwlAABosC6jWHd9uK7Jlde94sL0yDkjNLlvnIGVwV+Zzc4W8Ef1jdOsd9doa26Zaursum/+Ri3cmqfHzhupbpGMHgC8mStQb8H8dGm/Feq0fAc81qasUv19/kYt31Poui88yKJ7Thuii8ancM4OhkmKDtV710/Ww19t1rzf9kiSPl2TpWW7C/XYeSM0rX+3Q78BAK9gszuUUeS8+DI1NuyI9jtDk6Nd21+sy9LFE9h/AZ7C4XDoi3XZ+ufXW1wr0iUpKiRAfzllsC70g+PNVgfqkjO4nTFjhmbMmNHe9aAdpMSG6e1rJ+m1JXv06DdbVG21a09Bpc57YYmun9ZHt584QCEtaO0HAAA6RlWtTY//sE0v/7JL9vpFv2aTdN3RfXT7dPbTMN7gpCh9dssU/fPrLa4Tnz9tzdfJT/ysf50/QscPSjC2QABHrKLGOVYkPLhl+5qmK9QJ1AFPU1xZq/9+v01vLt3rOq6UnLPS7zt9iFJaOcMW6AhBAWb9/YyhGtOzi/768XqV1dQpu6Ral7+yXJdNStVfThmssKAjOk0NwEPkllar1maXpCPe94xOjVFqbJjSCiv1284CLdqWr2MHxrdnmQCOwNr0Yj3wxSatdOusaTGbdOnEVN0+fYC6hLe+I4U3alHL96eeekrV1S3/xfmFF15QWRkzLoxkNpt01ZTe+urWaa62Sg6H9OLPu3T604u1LqPY0PoAAPBXv+3Yp5Of/Fkv/dwYpg9OitJnN0/VnJmDCdPhMUICLfr7GUM176rx6hrhXJVeUFGrq+f9rjkfr6dNJ+CFHA6HKurHgUW0cIV6ZEiga6wYgTrgOWx2h95Znqbj/r1Qry9pDNN7dw3X3KvG6+UrxhGmw+OcMTJZ39x+tKb0a+zG9ebSNM188hf97tZdAYD3SSt0n59+ZPufQItZfzppoOv2P7/eIpud0WOAUXJKqjX7/TU689lfm4Tp0/p31de3TdMDZw7zmzBdamGgfvvtt7cqIL/rrruUn59/xEWh/fTpFqEPbzxKd88c5JqFuT2vXGc++6vu+XS9iitrDa4QAAD/kF1SpdvfW6NLXl6mvQXOXzSDApy/LM6/ZYqG94g+zDsAxjh2YLy+nTVN0wc3rgx4Z3mapv93kb5cl81sdcCLVNba1PAj29KW71LjKvXskmp+5gEPsCqtSGc9+6vmfLxeRZXOC9zCgiz688mD9M2saTqO1XzwYN1jQvXG1RP1wJlDFRLoPFe5t6BS57+4RI98tVnVVpvBFQI4EmkF7oF66BG/z6nDkzSy/vzIlpwyfbwqo821AWidaqtNTy2on5O+qnFOep9u4Zp75Xi9fvUEDUiINLBCY7ToN2iHw6ETTjhBAQEt+4W7qqrq8E9Cp7GYTbrxmL46bmC87vhgjTZklsrhcF4B+tX6HP355IE6f2yKzGbfnm8AAIARKmvr9OKiXXrx552qttpd94/v1UX/PHeE+naLMLA6oGXiIoL1vz+M01vL0vTwV5tVWWtTflmNbn57lY4fFK8HzhyqHl1YBQd4uor6+elS6wL1pOhQ7cyvUJXVptKqOkWHBXZEeQAOI6+0Wo99u1UfrmwaLpw+Mll/OWWQkqKPPMAAOpPZbNIfJvfStP7ddMf7a7QqrdjVWfOnrXn67wWjNKw7FxwD3sR9hXrPuPAjfh+z2aS7Zw7Wxf9bKkn6z3fbdPrIZLr5AZ2g2mrTO8vT9PzCncorq3HdHx0aqFnT++uyST0VaGnROm2f1KLfoO+7775WvemZZ56p2NjYIyoIHWdgYqQ+uWmKXl28W08u2K7KWpsKK2r154/W653l6XrwzGGsjgMAoJ3Y7Q59tjZTj369VTmljS1yY8ICdceMgbp0QioXs8GrmEwmXTapp44bFK/7PtugHzbnSZJ+3JKnJTsLNPvEAbpqSi8F+PEvV4CnK3cL1Fva8l3ab456aTWBOtDJiipq9cLPO/Xab3uaXKA5MCFS9585VJP6xB3i1YDn6t01XB/ceJT+98su/fe7baq12bUtt1xnPfurbjm+n24+rp9fn7gHvIl7oN7WkSOT+8bphEHxWrAlTzml1Xr119266dh+bS0RwEFUW216b0W6nlu4Q7mljUG6xWzSZRNTNcuP5qQfSocE6vBcgRazbjimr84c1V0PfbVZn6/NkiStSS/WGc8u1sUTUvWnGQP54QAAoA1W7i3SA19s0tr0Ytd9AfWrMG47oT9BBLxa95hQ/e8P4/TtxhzdN3+jcktrVGW16aGvNuuT1Zl65JzhGpkSY3SZAJpRUdPYRjc8uOWrfJLcAvXskioNTPS/9n6AEcpr6vTq4t3638+7VOZ2QUxkSIDuOHGALpvUkwvZ4PUaOmseO7CbZr+3VpuyS1Vnd+iJH7brmw05+sdZwzSuFwu3AE/nHqj36NL2jil/njlIP23Nk90hPf/TTl00PlWxZBZAu6qps+n9Fel69qedTRYDSdLJQxN1x4wB6u+Hrd0PpuWXpMOnJEaH6OmLR+viCSm677ON2p5XLodDentZmr5en627Th6kC8fRBh4AgNbILK7So19v0fz6C9YanDAoXn85dTDt3eEzTCaTTh6WpKP6ddV/vt2q15fulcMhbcou1VnP/aorJvfSHTMGKDKEi0cAT9J0hXrLfz6brFAvqT7EMwG0h2qrTW8u3avnFu5UYUWt6/6gALMum9hTNx3XV10jgg2sEGh/gxKj9OnNU/TMj9v17MKdstkd2pJTpvNeWKLzx/bQ3TMHKY5/94DHSq8P1BOjQtqlPfuAhEhdMC5F765IV1lNnZ5asF1/P2Nom98XgDNI/+D3DD330w5l7ff73YwhCbpten8NTaab9f4I1P3cUX276qvbpmner3v0xA/bVFFrU1GlVXM+Xq93l6fpvjOGakxqF6PLBADAo5VVW/XSz7v00s+7VFPX2IZzQEKE7j1tiKb172ZgdUDHiQoJ1P1nDtNZo7trzsfrtSWnTA6HNO+3PfpmQ47uOW2wTh2eJJOJizQBT1DRJFA/0hXqBOpAR7Ha7Prg9ww9tWB7k1VCFrNJF4zroT8e31/JMcxJh+8KCjBr9oyBmj4kQXM+Xq+NWaWSpA9WZui7Tbm66+SBumh8qiwsAAI8SnlNnQrqLwBLbWO7d3e3nzhAn63JUlX9hWZXHtVLvboe+Xx2wN/V1tn1wcp0PfvjgUH69MEJmjW9v4Z1J0g/GAJ1KNBi1nVH99EZo5L10JebXavq1maU6JznftNJQxP0p5MGql88rR0AAHBXXlOn137bo5d+3qWSKqvr/tjwIM0+cYAuGp9CG074hdGpXfT5H6fq1cW79fgP21RttSuntFq3vL1ac3vu0T2nDtZoLtIEDFdR2xioh7dmhnpUY4DHCnWg/dnsDn2+NkuP/7BNewsqmzx2xshk3X7iAPUmQIAfGdEjRvNvmaq3lu3Vv77dqrLqOpVUWfXXTzbo/RXp+sdZwzW8Byf8AU+RVtB+89PdJUSF6LppvfXUjztUZ3foX99t1bOXjGm39wf8RXFlrd5alqbXftujvLKaJo9NHxyv204YwH61BQjU4ZIQFaKnLh6tiyek6r75G7Qtt1yS9O3GXH2/KVfnje2hWdMHcDU0AMDvVdbW6fUle/Xiop0qqmwM0gMtJl15VC/dcnx/RYfS6hr+JdBi1g3H9NUpw5N0z6cbtGhbviRp5d4inf3cbzp9ZLLuOmlgu55gAdA6ZdVHFqi7r1DPLK5q15oAf1ZttenDlRn63y+7DgjSpw9O0B0zBmhwUpRB1QHGsphN+sPkXjp5WKIe+WqLPlmdKcm5AOiMZxfr8kk9dceMgfzeBXgA9/npPePa9/e964/pq7eWpamgolZfrsvWtVOLuFgbaKHd+yr06uLd+nBlhqqstiaPHT8oXrOm99eIHjHGFOeF2i1Q37Vrl2688UZ999137fWWMMjkvnH68tZpev/3dD35w3blldXI7pDe/z1Dn67J0hWTe+qmY/upS3iQ0aUCANCpqmqdbcZeWLTT1c5Mkswm6ezRPXTrCf3UM47VQ/BvKbFhmnfVeP20NU8PfblZO/MrJEmfr83StxtzdM3U3rrp2L7MVwcM0LTle8tPB8SEBSomLFDFlVYt2VWg7bll6p9ABzPgSJVUWfXm0r2a++tu7SuvbfLYUX3jdOdJAxm/B9SLjwzR4xeO0gXjUvS3zzZoe165HA7p9SV79dX6bM2ZOVjnjOnOiCHAQOlugXp7tnyXnMess6b3172fbZQkPfL1Fr13/SR+5oGDcDgcWr67UC8v3q0fNufK4Wh8zGySThqaqBuO6atRKTGG1eit2i1QLysr04IFC9rr7WCwQItZl07sqXNG99Dc33br+YU7VVZdp9o6u/73y269uzxdNxzTR1dP7a2wIBodAAB8W7XVpreXpem5hTu1r7yxNZLJJJ01qrv+eHw/9ekWYWCFgGcxmUw6flCCpvXvpneXp+nxH7arsKJWtXV2Pb9wp95fka7bGYsAdDr3QL01K9RNJpOundpb//5um2x2hx7+arPmXjWhI0oEfFpOSbVe/XW33l6WpnK3n0dJmtIvTjcf209H9etqUHWAZ5vcN05f3TZNry7erSd+2K4qq037ymt1xwdr9frSvfrrKYM1oXes0WUCfsl9hXpHdCS7aEKqXv11j3bvq9Dy3YVasDlP04cktPvnAN7MarPrq/XZevmX3VqfWdLksbAgiy4Yl6Krp/RWajt3kfAnJKE4pNAgi246tp8umZCq5xft1Lxf96imzq6ymjr9+7ttmvfbXt12Qj9dOD5VQQGcDAUA+JaqWpveW5Gm5xftVG5p0yD9tBHJuu2EfuoXzwo94GACLWZdPrmXzhjVXc/9tENzf92jWptdBRW1uufTDXrttz36yymDdezAbqwwADpBeU1jm7+IYEurXnvttD56e1maskqq9dPWfP28LV9HD+jW3iUCPmlHXrle+nmnPlmdKautcZmQ2STNHJakG4/py9xKoAUaRgydPjJZD36xSV9vyJEkrU0v1gUvLtH0wQm6e+Yg9YvnYmegM6V14Ap1yfmz/+eTB+rGN1dJkh75erOOHdiNi7MBOS/YfP/3dL2zPE3ZJdVNHkuMCtGVU3rp4gmpjEhpBwTqaJGYsCDNmTlYVx7VS0/+sF3v/54uu0PaV16jez/bqOcX7tT1R/fRRRNSFRLYuhMzANpm0bZ8PfzlZp0xKlk3H9fP6HIAn5BfVqM3luzRG0v3NpmRLkmnDk/SbdP7awCtboEWiw4N1JxTBuuyST316Ddb9MW6bEnS9rxyXTVvhSb0jtXsEwdoUp84gysFfNuRrlCXpJBAi/48c5Bue3eNJOmhLzfrqL5xnMgEDsLhcOi3nQWa99ueA9ptBgWYdd7YHrp+Wh/16sq4IKC1kmNC9fxlY13nQ7bmlkmSfticq5+25unC8SmaNb2/4iNDDK4U8A8NLd9DAy3qGtExY2JPGpqosT27aOXeIu3Mr9D7v2fokompHfJZgKez2R1auDVP7yxP049b8mR3NH18WPcoXTetj04ZnqRAfl9rNwTqaJWk6FD989wRunZaH/3nu62uK0GzSqr198836Zmfdujqqb11+aSezMUEOsmzP+3Q1twy/evbrTpjZHKHtFYC/MWOvHK9sniXPlqVqdo6e5PHThqaoFnTB2hwUpRB1QHeLyU2TM9cMkZXTSnSP77cpNVpxZKk5bsLddFLSzW5T5xuP3EA7TqBDlJe6xaoH8HortNHJOvVX/dobXqxtuaWcSITaEZJlVUfr8rQG0v3ald+RZPHIkMCdPmknrpySi+CPqAdHDOgm6b266qPVmboP99vVW5pjWx2h95elqZPV2fquml9dP3RfVp9ERmAlrPZHUovcgbqqbFhHdZ5zGQy6S+nDNK5zy+RJD3+wzadOSqZn2/4lcziKr2/Il3v/55+wGp0s0k6flC8rp3WRxN7x9IFsAO0+Ntm9OjRh/wfUFlZedDH4Hv6xUfo+cvGak16sZ5esF0LtuRJkvaV1+qxb7bq+YU7deVRvXTVlN6KDe+Yq9IAOGWXVLm2P1qVoVnTBxhYDeB9HA6HVuwp0ks/79QPm/OaPBZgNumMkcm6dlofDUkmSAfay9ieXfTx/x2lL9Zl6/EftrkChyW7CrTkxSWa0i9Ot08foHG9CNaB9uS+Qj3iCE4+ms0m3XvqYJ33gvNE5n+/36rTRyZxMTUgaVNWqd5Yukefrs5SldXW5LGEqGBdM7W3Lp6Qys8L0M4sZpMuGJ+i00cm69Vfd+v5hTtVXlOnylqbnlywXW8tS9PtJ/bXheNS6KoCdICc0mrXOJOOns08tmesTh6aqG825ii/rEYv/7Jbt03v36GfCRitzmbXj1ucq9EXbcs/YDV6UnSILhiXogvGp6h7TKgxRfqJFv8GfdZZZ3VgGfBWo1Ji9MqV47Upq1TPLdyhL9dny+GQyqrr9PSPO/TyL7t16cRUXXd0HyVEcfU10N4cDofyyxrnOn+4MkO3Ht9fZjNXoAGHU2ez69uNuXrpl11am17c5LHI4ABdMjFVV07ppaRoDkaBjmAymXT6yGSdMjxJ89dm6skftmtPgfMi3V93FOjXHUs0rX9X3X7iAI1J7WJwtYBvaEvL9wbjesXq1BFJ+nJdtvaV1+q5hTv155MHtVeJgFepqbPpmw05en3JXq3cW3TA45P6xOrySb00Y2gC7TaBDhYaZNHNx/XTReNT9FR9kF5nd2hfeY3++skGvbJ4t2ZNH6BThyfJwjkToN2kFXTs/PT93XXyQH2/OVc2u0Mv/rxT54zpTrdO+ByHw6ENmaX6bE2m5q/NUp7b+X+pYTV6gi6ZmKJjBsSzX+skLf4N+r777uvIOuDlhiRH6ZlLxmh2frmeX7hTn6zOVJ3doSqrTS8v3q3Xl+zVeeN66JqpvdW3W4TR5QI+o7ymTtXWxrbUGUVVWrq7QEf17WpgVYBnyy+r0fu/p+vtZWnKLK5q8lhSdIiuntJbF05IURSrh4BOYTGbdPboHjp9RLI+XZOlp3/crr31J2V+2b5Pv2zfp2MGdNPtJw7QqJQYY4sFvFx5jXPVbFCAWUEBRx7u3X3yIH2/MVe1NrteWbxbl0xI5UQm/MqOvHJ9tCpD769IV0FFbZPHIoIDdM6Y7rpsUk8NSIg0qELAf8VFBOv+M4fpyim99a9vt+ir9c5xlbvyK3TrO6v15A/bdOsJ/XXaiGQCCKAdNMxPlzonUO/TLUKXTEjVG0v3qrLWpvNfWKLXrp6ggYnsc+H99hZU6LM1Wfp0TeYBo4MkqXtMqC4an6Lzx6UoMZoFrJ2tVZek/z979x3eZn3v//8lyfLeeyfOdvZOTEJCBgkJm0BLS1vK4bQ95xdoIW1PS09bSlsOpf2eltIyTiktXSkUWvYIIYGEQPZetpM4ie3YlrdlyUPz94dsxSbKIrHl8Xxcly9Jt27JH0Gse7zu9/uzZcsWvf7663I4HFq0aJGuueaanhoX+qlhKdH6xW2TdN/Vo/S7Dcf0/PYytbs8crg9Wr21VKu3luqq0Sn68hVDNW9kClW0wCWq+cTVaZL00o5yAnXgE7xer7aU1OuvW09qzYEquT7RH2lsRqy+Om+Yrp2YQfUQECQhJqNunZatGydn6uXdp/Sb9UdUVu+76GVDcY02FNfoiuFJ+tr84Zo3Mpn5wIBPobNC/dO0e+8qJzFSd80dqv/bUCKHy6NH3ynUbz8/9XIMEeizmlqcem1fhf65s1x7PtHdSJJGpUXriwVDdfOUrEv+GwNw6fKSo/TkHdO0q7RBP3urUNtO1EuSjtXY9Y3n9+jxdUd078KRun4SwTpwKUp7OVCXpPuvHqUPj9ToRF2Lqqxtuu3pj/Xsl2doBlOGoR+qtbXrjb0VenVvhXaXNp7xvNlk0MIxqfrczFxdOTKFbVYQXfAe/ksvvaTPfvazioiIkNls1i9/+Us9+uij+ta3vtWT40M/lRUfoYduHK97Fo7Us5uO669bTsrWcfLmg6IafVBUo+EpUfryFUN1y9TsT91uEBjsAgXqbx2o1EM3jmNuPkC+E5//3FWuv209qWOfuLLTYJCuGpWif79ymK4YnkQ4B/QRZpNRn5meo5unZOmfO8v1m/VH/d0kPj5Wp4+P1WlMeoy+Om+Yrp+UyUUwwEXoPCaLCjNd8nutXDBCL+0oV53doTf2VequOfWaNoSTmBhYXG6PNh6p0T93ntLaQ76uDF2FGA26Zny6vjh7iGbmJbI/CfRBU3MT9MLXZuvjY3X69XtHugXr973QEawvGqHrJ2YyxzrwKXQN1HurY1FiVKhe+s8r9G/Pbde+8iZZ21z6wu+36vHPTdHScem9MgbgUjS1OPXeYYte21uhTUdr5f7kxOiSZuUl6qYpWVo+PkNxkZzn7wsMXq/3zP9TAUybNk0zZszQE088IZPJpEceeUS/+MUvVF9f39NjvCRWq1VxcXFqampSbGxssIczaDW1OvWP7WX60+YTKm/o3l43JjxEt8/I0ZcKhtImELhIb+yr0D2rd0uSws1Gf/v3R1dM0Gdn5AZzaD3G6XTqrbfe0vLly2U2szOBM3m9Xu0tb9LftpzU6/squk2LIEnJ0aH6zPQcfY72tEC/4HB59NLOcv1u4zH/HOudMuLCdffcPN0+M1fRYSFsI4DzGPX9t+VweTQmPUbv3Dfvkt/vr1tO6vuvHJAkTcqJ18v/eQVdyNBnXcw2orDKqn/uLNfLuytUazvzIub8jFitmJqlGydnKSUmrKeGDOAy83q92lxSp8feO6Jtx7uf085LjtI9C0boxskE64MVxxKfzo1PfKS9ZY0yGKTDP75G4eZLv3DzQtnbXfqPv+7Uh0dqJfnmlf7pTRP0+VkD85wo+jeLtU3vHqzSmoMWbSmpO6N7piSNSY/RTVOydMOkTGXGRwRhlIPPxWTIF1wWXFRUpBdeeEEmk+8L8Zvf/KZ++MMfqrq6WqmpqZc2Ygx4cRFmfWXeMP3b3DytPWTRHz86rq0dO67NbS498+FxPbvpuBbnp+nLc4aqYBiVgsCF6FqhfvuMXD338QlJ0os7ygdsoA6cTXVzm17bU6GXdparsKr5jOdn5SXqjtlDdM249EuaNxZA7woNMerzs3L12Rk5WnuoSk9vKPG32q1satNP3zysX687ojtmDdEXZmYFd7BAH+Z0e+Rw+S4yu1ztqG+fkaM/bz6hYotNe8sa9fq+Ct04mb9D9E8lNTa9tb9Sb+6v0uFK6xnPJ0WF6sbJWVoxLUvjMuOCMEIAl8pgMOiK4cm6YniyNh+r02PvFfvPTx6vteubL+7V4+uP6N+vHKZbp2YrIrT3gkGgvzrVUTyXGhPWq2G6JEWFhejZO2fov17aq1f2VMjjlb738n7VNLfr64tGkC8g6I7X2rXmYJXeOVAVcMogydft+YbJmbppcpZGp8f07gBxUS74KLqlpaVbOh8aGqrw8HDZbDYCdVwwU0c7tGvGp+tgRZP+9PEJvbKnQg6XRx6v9O4hi949ZNGw5Ch9dkaOVkzLVnI0V3sDZ9O1WmLhmFR9fKxWxRabdpxsUEmNTcNSooM4OqDntTndeveQRf/aVa6NxTX65MWdMeEhWjE1W3fMytXINHZKgf7Mtx+ZoaXj0rX9RIN+t/GY3jtcLcl3gebTG47p2U0lmpJo1JAKqyYPSQryiIG+pXP+dEmXbcqtEJNR/33tWN35h22SpEffLtTScem9fjIV+LSO19r11v5KvbGvMmCIbjYZtGhMmlZMy9ZVo1OYZgQYQAqGJ6lgeIG2lPhawW8uqZMknaxr0Q9eOaBfrS3WF2cP0ZcKhiiJc5NAQC63R3V237nJtNjwoIwhNMSoX35mslJiwvTMh8clSb96r1jVzW368Y3jmW8avcrj8Wr/qSa9d9iiNQerVGyxBVwvOyFCS8f5crJpuQl0+eonLuoo+ve//72io0+HMy6XS88995ySk5P9y77+9a9fvtFhQBuXGaef3zpJ37lmjP6+rVR/2XJSFqtvA1xSa9cjbxfqF2uKdPXYNN0+M1dzRySzAQQ+oWuFempsmG6blqOH3zosSXppZ7n+65oxwRoa0GO8Xq+2n2jQv3aV6819lWruEhB0mpwTr8/PzNV1kzIUGXp5QgMAfYPBYNDMvETNzEvUEUuznvmwRK/srpDD7ZHT7dW2GqNuemqLpubG60sFQ7VsQrrCQgj3AFuX7eXlqlCXpPmjUjR/VIo2FNeooqlNz246rpULRly29wcut5N1LVpzuEZv7qvUoQAhuuSbwuCWjnabCVGhvTxCAL1p9rAkzf5qkraW1Ok3649q01Ff6+h6u0O/XndET284plunZevfrxymvOSoII8W6FtqbQ51TiicGhOcQF2SjEaD/vvasUqNCfefF/3b1lLV2Rx67PbJXOyJHlVra9eHR2q0oahGG4/Uqt7uCLjemPQYLRmXrqXj0jQ2I5YOCv3QBc+hPnTo0PP+DzYYDCopKbksA7tcmEO9/3C6PXr7QJX+vrXUf1VoV1nxEfrM9BzdNj2b+SOADnf9cZveL6qRJO38/mJ5vNLsR9bJ7fEqPTZcH3134YC7EIU5rQavYzU2vbanQv/aXa6y+tYzns+Kj9DNU7J089QsDac7AzCoVFvb9MePT+ivW06qua37RTZJUaG6fWaOPj9riLLYh8QgVlTVrKWPbZQkfWZ6tn5+66TL9t5HLM265tcfyu3xKirUpPe/fVVQT6oCXXm9Xh2ssGrtwUr9c+tRldsDHx9NyonXtRPStWx8hnISI3t5lAD6igOnmvTMhyV6Y1+l3F1aoBkM0pKxafrqvGGaNiQxiCNET+F808XbV96oG377kSTpjlm5evjmCUEekfTy7nJ9+8V9/vmpZ+Yl6pkvTVdcBP9PcXm43B7tLmvUhqIabSiu0f5TTQHXMxikqbkJWjouTUvGpmsoF2X1ST0yh/qJEycudVzAOZlNRt0wKVM3TMrUiVq7XthRphd3lPtbWp9qbNWv3ivWr9cVa/6oFH12Ro4WjEml4giDWk3H34fJaFBCZKiMRoMWjE7Ve4ctqrK26cMjNbpqNNNyoP86VmPTW/sq9eb+yoDzokeFmrRsQoZumZql2XlJtEgCBqnU2HB955ox+trcIXr4b2u11x6n4mpfa7U6u0NPvH9MT31wTIvz0/SlgqGaMyKJq8Ex6HSvUL+8JxRHpsXoczNz9NctpbI73Prlu8X62YqJl/V3ABejzenW5mN1eu+wResLq1XZ1NbxTPfv/knZcbp2YgYhOgC/8Vlx+vXtU/Rf14zRHzYd1/PbfNs2r1dac9CiNQctmjYkQf8+N09Xj01TCFNBYBDr7DYrBbdCvaubp2QrMSpM//nXnWpxuLXteL1ufvIj/ffyfC0ck8pxIC6a1+tVSa1dW0rq9NHRWn14pPaMC/k7RYeFaM6IJM0flarFY1P7zN8FLg96oKJPGpocpe9cM0arrh6l9YXVen5bqTZ0zI3r8UrvF9Xo/aIaxYaH6NqJGbpxcpZmDk0kSMGg09nyPTk61P/v/7bp2XrvsEWSr+07gTr6m/OF6AaDNHdEsm6ZmqWl49Jp6Q7ALyosRHPTvXp4WYF2lzfrz1tOas2BKrk8Xnm80ruHLHr3kEXDUqL0hVlDdPOULFr5YtCwdwvUL/9FyfcvHqVXd1eoud2lF3aU6faZuZqcE3/Zfw9wNtXWNq0vrNZ7h6v10dFatTrdAdebmBWraydmavkEQnQAZ5cVH6EfXDdWX180Uqu3luqPHx1Xdcc5mJ0nG7TzZIMy4sJ1x6xcfXZGrlJimGcdg4/F2ua/nxbbd/4G5o9K0d+/Mlt3Pbdd9XaHSmrsuvtPOzQ1N17fXjpGBcOTgj1E9GFdA/QtJfXaUlLXbdrVTxqbEav5o1N01agUTR2SIDMXWg1YF3wGevPmzaqrq9N1113nX/bnP/9ZDz74oOx2u2666Sb95je/UVhY3/niRP9nNhm1dFy6lo5LV0Vjq17cUa5/7CjTqUZfq19rm0t/31amv28rU2ZcuG6YnKWbp2RpdHpMkEcO9DyPx6tam29Olq4HbgvHpCopKlR1dofePWRRU4tTcZG0NULfdr4QXfLNi37thAxdNylDGXG0bQZwdgaDQbOGJWnWsCRZrG36+7ZSrd5a6j8JWlJj14/fOKSfvV2oq8em6bbp2bpyZMqAmyYF6KproB51GedQ75QUHaZ7Fo7QI28XyuuVPve7LfrpTeO1Ylr2Zf9dgCQ5XB7tKm3QR0drtbG4RnvLA7fbDA0x6orhSbpqZJK8FQf0hZtn08oXwAWLizDrP68arrvn5unVPaf0zIclKrb4OiFVNrXp/71brF+vO6LlEzL0pYIhmpqbQAUsBo3qLiFjah8K1CXfVC7//M8rdN/zu/37CLtKG/W5Z7boypHJ+taS0ZrExZ/QxQfo8ZFmXTkyRfNHpWjeyGSlxlKFPlhc8FH0j3/8Y1111VX+QH3//v26++679eUvf1n5+fn6xS9+oczMTP3oRz/qqbFikMuMj9A3Fo/UPQtH6KOjtXpl9ym9c7BKLQ7fVecVTW16esMxPb3hmPIzYnXT5EzdMDmT0AUDVkOLwz+fV0r06Z1Ws8moGydn6Q8fHZfD5dFre0/piwVDgzRKIDCX26NdpY1ad9ii9w5bdKzGHnC9zhB92YR0ZSdQQQTg4qXFhuu+xaO0csEIvXvQoj9vPqGtx+slSQ63R2/u913MkxEXrhVTs3Xb9GwNSWJuMww8th4O1CXpy3OG6s39ldpX3qRWp1vffHGvth6v00M3jFdEKFN14dJ4vV4VW2zadLRWm47UaOvxev/5gE9Kjg7TwjEpWpSfprkjkhUVFtIxN+6BXh41gIEiNMSo26bn6NZp2dp4pFZ/2XxC6wqr5fVKTrdXr+6p0Kt7KjQ2I1Z3XjFEN0zKYtuHAa+6S4V6X2xtnZccpVdWztGagxb977tFOtIxLdiHR3xtu68Zl65vLhmlkWkU5w0mzW1O7Stv0q6TDdpd1qjdpQ1qaHGedf2oUJNm5CVq9rAkzR6WpAlZcVyMP0hd8FH0nj179JOf/MT/+Pnnn9esWbP0zDPPSJJycnL04IMPEqijx5mMBs0blaJ5o1L0U4dLaw9Z9MruU9p4pNYfLh6utOpwpVU/e6dQM4cmavmEDF0zPl1pXC2EAaRz/nRJZ7QWu216tv7w0XFJ0os7ywnU0SdY25zaWFyjdYer9X5RtRrPsrNKiA6gJ5hNRl07MUPXTszQEUuzXtxZrn/tKvd3e6lsatNv3z+q375/VLOHJeq2aTlaNoFpJTBwdJ9DvWf+XYeFmPTCVwv04GsH9I8d5ZKkf+wo156yRj3x+amcrMRFq7a2dQTotdp0tLZbJdwn5WfEanF+qhblp2liVhxTwgHoEQaDQfNH+SoTy+pb9LetpXphe6k/jDlUadV3/rlf//NWoW6blq07Zg9RXjIXa2Jg6rpd7qvn3Q0Gg64Zn66rx6bp1T2n9Kv3ilVW7+t++87BKq05VKWbp2Tp/sWjmApmAPJ4vCqptWnXyUbtLmvQ7tJGFVma5fWe/TWfDNDHZ8YqhDbu0EUE6g0NDUpLS/M/3rBhg5YtW+Z/PGPGDJWVlV3e0QHnERkaohsnZ+nGyVmqtbXrzX2Venn3Ke0pa5Qkeb3S1uP12nq8Xg++dlDThiRo2fh0LZuQoax4KtfRv3VtPZMc3T1Qz8+I1fisWB04ZdW+8iYVVTUzFQKCoqy+Re8dtmjd4WptPV4np/vMPVajQZo2JEFLxqYTogPoFSPTYvS95fn69tLRWl9YrRd3lOn9ohr/xZm+Nm++/cdrJ2ToxsmZmjUsiavQ0a/1dMv3ThGhJv381kmalZek779yQK1Ot4otNt3w249oAY9z8nq9Kqtv1dbjddp+ol7bTzToeG3gLkaS76LiK0cka86IZM0dmdxnT+QDGLhyEiP13WVjdN/ikXpjX6X+svmEv7V0U6tTv990XL/fdFwz8xJ1+4wcLRufQdU6BpTOOdRNRoOSokKDPJpzMxkNumVqtq6bmKkXdpTpN+uOqLq5XV6v9K9dp/T63gpdOyFDyydkaN6oFIWb+Vvtb1xuj47V2HWwokkHK6z+2+Y21zlflxBp1uSceM3MS9LsYYkanxXHPOgI6IKPotPS0nT8+HHl5OTI4XBo165deuihh/zPNzc3MwcVgio5Okx3XjFUd14xVMdr7Xp1zym9tqdCJV0OwHeebNDOkw366ZuHNSknXsvHp2vZ+AzlJhHeoP+pPUeFuiTdNi1HB04dlCS9uKNM379ubK+NDYNXc5tTW0rq9eGRGn14pPasJ0Gjw0I0b1SyFuen6arRqUrs4wdeAAYms8mopePStXRcuqqtbfrX7lP6x44ylXRMQ2Frd+mFHWV6YUeZ0mLDdP3ETN00JUvjMmOZGxP9jq39dGvsqLCeP0G4Ylq2JuXE6f/72y4VW2y0gMcZPB6viqubtb3jIvjtJ+plsZ69Aj3CbNLsYYmaMyJZV45M0ai0aL6LAfQJ4WaTbp2WrVunZWtvWaP+vPmkXt9XIYfLI0nadrxe247X68FXD+qGyZm6fUauxmexP4n+r7NCPSU6rN90hgkNMeqLs4fo1qnZ+vPmE3pqwzE1tjjldHv1yp4KvbKnQpGhJi0ck6pl4zO0YEwKXcv6oDanW4VVzV3Cc6sKK61q7/jePRuT0aAx6TGakhuvKTkJmjokQUOTIvk+xgW54G+C5cuX67vf/a4effRRvfLKK4qMjNSVV17pf37fvn0aPnx4jwwSuFh5yVG6b/EofWPRSBVbbHprf6XePlCpYovNv87eskbtLWvUI28XalxmrJaMTdei/FROkKLf6FqhHihQv3Fyph5+87Acbo9e2XNK31k2hqvrcNm5PV7tK2/smH+qRrtLG+XyBO6blJ0QocX5aVqUn6pZeUkKDeHfI4C+IzU2XP8xf7i+Nm+YdpU26B/by/XGvgrZO+bntVjb/VVGw1KidNPkLN04OZP51tFv2Huh5fsnjUiN0asr5+pHrx3UCzt8He1oAT94NbU4tbe8UXvKfD87TzaoqfXs81WaTQZNyo7X7GFJmjsyWVNzE9h/BNDnTcqJ1//mxOu/r83XSzvL9ML2Mh3ruFizud2lv20t1d+2lio/I1a3z8jRTZOzFBdJkRr6H5fb4y/2SY0987xkXxcRatLX5g/X52bl6vcfHtcfPzrur2Rucbj1xr5KvbGvUuFmo+aPStHyCRlaOCZVMeH8vfamVodbx2psOlLdrGKLTUcsNh2tblZpfYvOcvqxm7TYME3MjteU3HhNzU3QxOw4LpDAp3bB/3J+8pOf6JZbbtH8+fMVHR2tP/3pTwoNPV1N9oc//EFLlizpkUECn5bBYNDo9BiNTo/R/VeP0tFqm945UKm39lfpUKXVv17nVUy/eq9Y6bHhWpifqkVjUjVnRDLtXdBndQvUo8/ccY2PDNXVY9P05v5K1doc+qCoRlePTTtjPeBieL1eldTatbWkXpuO1uijo3VnPREaYjRoam6C5o9O0aL8VI1Oi+GCJQB9nsFg0LQhiZo2JFE/umGc3jts0at7KrShuNo/bUVJjV2/XFusX64t1uSceN04OVPXTshQKu2G0Yf1Vsv3T4oINenRWydq1rBE/ffL3VvAP3TjON06NbvfVDThwrW73Dpc2aw9pQ3aW96kvWWN3brHBRIZatK0IQmaOTRRM/ISNTknnuNxAP1WYlSovjpvuL5ype9izee3lemNfZVqdfou1jxcadWDrx3Uw28d1jXj0nXrtGzNGZHMFEPoN2ptDv881Kkx/fc4KDbcrFVXj9I9C0boo2O1ent/pd49ZFFji+9cV5vTozUHLVpz0KJQk1FzRyarYFiSJmbHaXxWXK/uVw9UHo9X1c3tOlFnV2ldS0eA7gvRyxtazznfeVdDkyI1LjNO47JifbeZsWdMkwpcigv+a09OTtbGjRvV1NSk6OhomUzdD2pefPFFRUdHf+qB/OxnP9MDDzygb3zjG3rsscckSW1tbfrmN7+p559/Xu3t7Vq6dKmefPLJbnO5AxdjRGq07lk4UvcsHKmTdXa9faBKb++v9M9vJElV1jat3lqq1VtLFW42au6IZC0c46uoZE429CXnq1CXpFunZ+vN/ZWSfG3fCdRxsbxer45U27S1pE5bjtdra0l9t+kGPmlYcpSuHOlrwzl7eFKvVcABQE+ICDXp+kmZun5SphrsDr19oEqv7jmlrcfr/et0Vlr++I1Dmj4kQcvGZ+ia8enKjI8I4siBM9mCUKHe1S1TszUxO04r/7ZbRZZmtTrd+q+X9unxdUd0+4wc3TY9h+OtfqrV4VaRpVmHK606VGHVvlNNOlxhlcN97pabCZFmzRiaqJl5iZoxNFHjMmMVQkctAANM14s1f3j9WL2xr1IvbC/TnrJGSZLD5dFreyv02t4KpcaE6YZJmbp5apbGZtBBE31bdXOb/35aP6xQ/6TQEKMWjE7VgtGpetjt0daSer11oFLvHqxSrc0hSXK4PVpfWK31hdWSJINBGpESrYnZ8ZqUE6eJ2fHKz4hRWAgXBH5Su8utysY2X2he36KTdS06WWfXyboWlda3nLdVe1cRZpNGpEZrdHqMxmX6wvP8jBi6B6DHXfRRdFxcXMDliYmJn3oQ27dv1//93/9p4sSJ3Zbff//9evPNN/Xiiy8qLi5O99xzj2655RZ99NFHn/p3AZ2GJEXpP+YP13/MH67KplatO1ytdYct+uhYnX+OozanR+8drtZ7h6ull6UJWXGaPypFV45M1hTazSHIas4zh7okXTkiWakxYapubtf6wmrV2tq5Mg/n5PF4VWRp1taSOm3tmM+y3u446/pxEWbNHZGsK0cma+7IZGUnRPbiaAGg9yREherzs3L1+Vm5qmhs1et7ffPrHe7oeuT1SttPNGj7iQb9+I1DmpwTr+UT0rVsfIZyEvluRPDZHcEN1CVfC/hXVs7p1gK+vKFV/+/dYv3qvSNaOCZVn5+Zq3mjUqjQ66Nqmtt1qCM4P1xp1aFKq0pqbOdtuRlqMmpsZqwm58Rrck68JuXEM18lgEEnJtysz83M1edm5qqoqlkvbC/Ty7vL1dBRCVvdfHqKoVFp0bppSpZumpzFhZrokyzW0+cl+3OFeiDmjkr0uSOT9ZMbx2v7iXq9c6BKbx+o7Pa5vV51VFLb9M9d5R2vNWhMeqzGZ8VpaFKkshMilZ0QoeyECCVGhQ7IfZ92l1uWpnZVNLWqqqlNlU1tqmpqVUVTm//xuYpzziYq1KQRaTEamRrt+0mL1sjUGGXFR9DhCkER9LIxm82mO+64Q88884x++tOf+pc3NTXp2Wef1erVq7Vw4UJJ0h//+Efl5+dry5Ytmj17drCGjAEoIy5CX5g9RF+YPUQtDpc+OlqndYctWldY3a0KeP+pJu0/1aTfvn9UUaEmFQxP1rxRvkpMTgagt3X+2ww3G896UjTEZNQtU7P19IZjcnm8emX3Kf37lcN6c5jo45pandpT1qhdJxu0u6xRe0obZG1znXX96LAQTR+aoNnDkjR7WJImZMVxwhvAoJMZH6GvzR+ur80frmJLs97YW6G3D1TpSLXNv05n5fr/vFWo8VmxWjY+Q8snZCgvmTnXERy2dl+LWYPB11o7WDpbwF8zPl3PfXxCG4/UyOuV3B6v1h6yaO0hizLjwvWZGTn6zPQcQoQg8Hq9qrM7dLTjBPGxapuOVttUZGnudnx8LsNSojQ5O16Tc+M1KTte+RmxXJAOAF2MTo/RD68fq+8sG633C6v1r12n9H7R6SmGii02/fydIv1iTZFm5SXqlinZumZCumKpwEQfMdAq1M/GZDT4z4H98LqxKqxq1r7yRu0tb9K+8kYVVTXL1eXKQqfb688QPinCbFJ2QoRyEk+H7OlxEYqPMCuuy09shDko59q8Xq9anW7Z2l2yt7vV3OZUvd3h/6mzO1Rv67i1t/uXNZ/jPOL5hJqMykmM0JCkKOUmRmpoUqSGJkdpZFqMMuPCyVvQpwQ9UF+5cqWuvfZaLV68uFugvnPnTjmdTi1evNi/bMyYMcrNzdXmzZvPGqi3t7ervf30AZ7V6qsYcTqdcjoDz/EKdGU2SFeNTNRVIxP10HVjdLDSqveLarS+qEYHK5r969kdbr132KL3DlskSdnx4b4r14YnqWBYomIj2MFFz+o8mZUcHSaX6+w7LjdPStfTG45J8rV9/9Ks7H69M9L5Xc53+sXzeLw6VmP3BeflTdpd2qhjtfZzzkUUGx6i6UMSNDMvQTOHJig/PaZbK06P2yWPuxcGDwAXIBjbiLzEcN27YJjuXTBMR6pt/vn1Ci2nw/UDp6w6cMqqX6wp0oiUKC0ck6JFY1I1KZuLktB7bG2+v4vIUNM59x17y9zhCZo7PEHlDa16adcpvbTzlCwd+7cVTW167L0jenzdEc0bmawbJmVo+pAEZcQNrOqnYHO5PTrV1KaTdS06VmPXsRqbjtXYdbTarsbWC/seNZsMGpkarfyMGI1Jj1F+x88Zx8Net5zO4O40chwBoC8ySlo0OlmLRierscWptw5U6bW9ldpZ2ijJVwG7paReW0rq9YNXD2j+qGQtH5+uBaOTFRka9FP7AwrbiYtT2dDiv58YGTJo/ruNTInQyJQIrZiSIUlqd7p1uKpZ+09Ztf9Uk/adsqrkLOfaWp1uf0X7+USHhSguIkSx4WbFRYQoMjREISaDzEaj79bUcWs0KMRkVIjRoBCTQfJKTo9Xbo9XLren+313x32PV3aHLzT3hecu2R1u2dtd5+06dLGMBik1JkzpceFKjw3XkMRI5SZGKDcxUkOSIpUWE3bWavO+cMyCge9ivrsMXu+5TqP3rOeff14PP/ywtm/frvDwcF111VWaPHmyHnvsMa1evVp33XVXt3BckmbOnKkFCxbo0UcfDfieP/rRj/TQQw+dsXz16tWKjKTVIi5Nk0MqbjKosNGgwiaDbM7AX/YGeZUdJY2M82pErFfDY70KZ+oUXEYuj/TNrb4Dp6HRXt0/4dwnp36136QTNt+/129NcCknuseHiCDzeqW6dqnMZlCZ3aAyu+9+q/vcwU2M2au8GN/31ohYrzIjfTu/AICLU90q7a03aG+dUWX2wF+k0SFejUvwanyiV6PjvApjfxE96MGdJjU6DIoze/Xj6X3vaji3VzrcYNDH1QYdajDIqzP/buJDffspnT9ZkRJTbp9bu1uqbZNq2wyqa/fddj5uaJc8Af47n01kiFdZkV5lRUnZUV5lRnqVFiFReA4Al19tm7SjxqAdtUbVtJ35XR1q9O1HTknyKj/eqyA2n8Eg9fwxozZX+3YCvj3RpWwacfm1uaWqFqm+3bf/Vd9uUH1bx2275PIOjBNtESavos1StNm3nx4fKsWHddx2PI4JlUwD4+NigGppadHnP/95NTU1KTY29pzrBu0ytrKyMn3jG9/Q2rVrFR5++a4yf+CBB7Rq1Sr/Y6vVqpycHC1ZsuS8/zGAi+HxeFVoadamo3XadLROO042+NsyedURXtkNWl/haw0zLjNGs/MSNTsvUVNz4xUVpHkLMTBUNrVJWzdKkkbmpGn58snnXL85tVzff/WQJGm3M0OfWzC+33ZRcDqdWrt2ra6++mqZzf3zM1xuXq9XpfWtOlBh1YEKqw52/JyrdbskhRgNys+I8c1lmR2nKblxyo6P6NcdDAAMbn11G1HW0KI1B6v13uFq7Spr9Fcr2FwGba0xaGuNFBpiVMGwRC0ak6IFo1OUHkslLi6v/961XpJLSXFRWr58brCHE9D1HbeVTW16adcpvbjzlG+/t0Ojw6DddQbtrvM9DjcbNTErTlNz4zUlN16j06KVHhs+aDo/tLs8qm7unKeyXVVW3zyVVdZ2VXbMV1lnd1z0+6bFhml4SpSGp0RreEqURnT89Pd5P/vqNgIAzuZL8h3v7y1v0qt7K/X2AYv/e93hOb1NjAo1adGYVC2fkKa5I5IVxpVOnwrbiYvz8l92SdW1kqRbli1ScvTAbft+OXk8vql1yhtbdaqhVdXN7Wpqdamp1ammVqesbU41tbpkbXWqqc0pa6urW0v5yy0q1KSosBD/bXRY5+MQRYWZlBAZqsQos5KiQpXY5Sch0iwzV7ZiAOjscn4hgpbo7dy5U9XV1Zo6dap/mdvt1saNG/Xb3/5Wa9askcPhUGNjo+Lj4/3rWCwWpaenn/V9w8LCFBZ25pe32WxmQ4jLblJukiblJmnlQqnF4dLWknptKK7R5mN1KrKcbg/v9ni1r9yqfeVW/e7DEwoxGjQxO06zhiVpxtAETctNVFwk/z5x4ZraTrdVSo0NP+/3241TsvXTtwrV5vRoXWGNFj+2SfctHqXPz8rttzs/g/V7vbnNqWJLswqrmlVY2ayiqmYdrrJe0HxFKTFhmpITr6lDEjQ1N0ETsuIUwWXsAAagvraNGJYap/9MjdN/LhipOlu71hdWa93ham08UqMWh69S2OHyaENxrTYU10o6rPyMWM0flaKrRqdo2pCEfru9Rt/g9XrV4vDtK0SH962/j0Byk81atWSMvrF4tDYfq9PW43XaebJBe8oa/X8zktTm9GjbiQZtO9HgXxZiNCgjPlw5CZG+n8QIZXfc5iREKjn67G0l+4JWh1u1tnb/3JS1NofqbL77dR1zVtba2mWxtqnWdvFheaeYsBDlJvlabeYmRvmC89RoDU+NHvDz8/a1bQQAnM+MYSmaMSxFP7rBo63H6/XGvgq9c6BKDS2+NrV2h1uv7avUa/sqFRMeoiVj07VsfLrmjkxWuJlj/ovFduLC1HTsh5iMBqXFRfXp/au+JjMsVJmJF9Y+1Lcf71aLwy2XxyOX2yun2yOXp+PW7ZXL42vn7nJ7ZTDI3/49pKM9fOet2WiUqaNNfESoSVGhIfx/w6B3Md/3QQvUFy1apP3793dbdtddd2nMmDH6zne+o5ycHJnNZq1bt04rVqyQJBUVFam0tFQFBQXBGDJwTpGhIVowJlULxqRKkups7dp6vF5bSuq0+Vhdt7lRXB6vdpU2aldpo57qWDYqLVrThiRq+pAEzRiaqJxEqkRxdjW205U6KTHnvwI0Jtys/7l5gr77r/1yuDxqaHHqwdcO6k+bT+iBZflanJ/Kv7c+xuHy6GSdXYVVvtC8sMqqwqpmlTe0XtDrU2PCNCErTuOz4jQhK04TsuOURrUjAARdUnSYbpueo9um56jN6dbmkjq9d8ii9w5bZLGenu7qcKVVhyutenrDMUWHhWjOiCRdNTpV80elKDM+IoifAP1Rq9Ptnw8xqh/Nt2oyGjR3ZLLmjkyW5Jvzu7CqWbtKG7TzpO/nk/tGLo9XZfWtKqtvlVR3xnuGhhiVEGlWbLhZsRFmxYaHKDbCrJjwkC7LzIoOD5HZaJCp44Sk0eA7GWnqWGYyGhTScQLS5fHI4fLK4fbI6fLI6fb47nec8HS6PWp1uNXc5pKt3SVrm1O2jvvN/lunmttcand5Lst/O998leHKSYxQbmKUhvjD80gNSYpSQqSZ/X8A6GdCTEbNGZGsOSOS9eMbx+vjY3V6Y2+F1hys8neoa25z6Z+7yvXPXeWKCjXpqjGpWjouXQtGpyhmgF8whd5V3ew7dknp4xcr9ncGg8FXMU63WyDogvZXGBMTo/Hjx3dbFhUVpaSkJP/yu+++W6tWrVJiYqJiY2N17733qqCgQLNnzw7GkIGLkhQdpuUTMrR8QoYkqaa5XVuP1/kD9mM19m7rF1tsKrbY9PdtpZJ8Ien0IQmaPjRR04YkKD8jRmEhXFUKn5rm0yfcLyRQl6RbpmZr1rAk/eKdQr2yp0KSVFJj11f+vEMFw5L039fma3xWXI+MF2fX1OLU0RqbjnX+VNtVUmPTyfoWuS+wpVN6bLjGZ8WeDs+z4pRKeA4AfV642aQFo1O1YHSqfnrTeB2ssOrdQxZ9UFStfeVN/vVs7S6tOWjRmoMWSb4LMTvD9WlDEqg8wnnZ2k93sunPJ+NCTEaN77hg8EsFQyVJFmubdp1s0O6yRp2ss/vC9IaWs3bvcbg8sljbu13A0t+YjAalxoQpPS5cGXHhyoiLUEZceMdj3/3UmDCF0NkCAAYss8mo+aNSNH9Uih6+eYI2Ha3RG3sr9e4hi3+7b3e49ea+Sr25r1KhJqPmjkzWNePStXhsmhKjQoP8CdCfudwe1dp8+1KpsbR6BzA49Okj6V/96lcyGo1asWKF2tvbtXTpUj355JPBHhbwqaTEhOm6iZm6bmKmJKm6uU07TjT4fk7W62CFtVt4VtPcrrcPVOntA1WSpFCTUfmZsZqSE6/JOfGalBOvoUmRVBUMUt0C9YuYoygrPkKP3T5FX56Tp4ffPKTtHe0xN5fU6frfbtItU7L17aWjlR5HGHs5WducKq1r0cm6FpXWt6i03q5jNb7g/GJadUaFmjQ6PUaj02OVnxGj0WkxGpMey5QRADAAGAwGf1C46upRqrW168MjNfqgqEYbi2v8LT2l0xdi/m5jicJCjJqZl6g5I5I1d0SyxmbEUiGCM9jbT7dJjw4bWBdgpMWGa9mEDC3ruJC5U1OrU2X1LSpvaFF5Q6vK6ltU1tCqisZW3/yUrU7Zu7SPDxajQYoOC1FMeEeVfIRZydGhSooKU1J0qJKiQpUUHdbtNi7CzN85AMAvNMSohWPStHBMmtqcbn18rFbvHKjS2kMW/z6kw+3R+sJqrS+slvFf0sy8RC0dl67F+WnKSYwM8idAf1Nrc8jbcRo7NYZziAAGhz4VqH/wwQfdHoeHh+uJJ57QE088EZwBAT0oNSa8WwV7i8OlPaWN2nGyQdtP1Gt3aWO3ShKH26O9ZY3aW9boXxYfadakbF/APjknXhOz45R0EeEq+q9PU6He1eSceP3jawV650CVfvZOoU7Wtcjrlf65q1xv7q/Qv88dpmUT0jU6LYbKlgvQ5nSrqqlNFU2tKq9v1cl6u0rrW1VaZ9fJ+hY1dglBLkS42ajhKdEalhKt0WnRGp0eqzHpMcqKj+DkKQAMEsnRYbp5SrZunpItt8er/aea9EFRtT4oqtHe8kb/Cax2l0cfHqnVh0dqJUkJkWZd0RGuzx2RzAlSSJLsA6RC/WLERZgV13GRytk43R7Z2nwt2K2tnbdOWTvar7s9Xrm9XrndXrk8Xnm8vlt3lx+P16tQk1HmEKPMJqPCQowymwwym3yPfc8ZFB5iUkxHK/nosBDFhocoOjxEEWYTF0kDAC6bcLPJH6673B5tO1GvNQeqtOagRVVW3/SBHq+0paReW0rq9dDrhzQmPUZXj03Tovw0TcyK47wDzqu6+fRUlGlUqAMYJAbHkTTQD0SGhuiKEcm6YoRvfkC3x6vCKqt2nmzQntJG7SlrVElt9zbxjS1ObSiu0YbiGv+yzLhwjeto+Tw+K1bjM2n9PBDV2C4tUJd8lXDLJmRoYX6q/rL5pB5fd0TWNpfanB799v2j+u37RxUZatKk7HhNyY3X1NwETcmNH3QXbbQ53aq1tctibVNFY5sqm1pV0dimisZWVTb5Hl9MlXlXKTFhGp4SpeEp0b6f1GiNSI1WRmw4B7AAAD+T0eC/gPK+xaPUYHdo45EabTpSq01Ha1XZdPqEVkOL09/aU5KGJEXqiuFJmj0sSbPykuhCM0h1vVA3OpzTAJ3MJqMSokKVQNtbAMAAFGIy6orhybpieLIevH6c9pY36p2DVVpzoEon6lr86xVWNauwqlm/WX9UqTFhWpSfpqvHpuqK4clMLYSAuk6dQ4U6gMGCI2mgjzIZDRqXGadxmXH6UoFvWVOLU3vLfeH63jLfbZ29e5BX0dSmiqY2rT1k8S9LjQnztxAdnxmrsZmxyoqPoBKiH+taoZ58iQF3WIhJ/37lMK2Ymq3H1x/RXzaflKtj+oEWh1ubS+q0uaTOv/6QpEh/uD4uM07pceFKjg5VWEj/OMjyeLxqbnOpsdWhxhan6lscqmluV62tXTXNXX5s7aptbpf1LPNvXgiDQcqMi1BOYoSGJEYpNylSuYmRGpIUqSFJUYqLoFU7AODiJUSF6sbJWbpxcpa8Xq9Kau366GitNh2p1eZjdWruEp6e7Jhy5O/byiT5tuOz8hI1Ky9Js4YlKjuBCvbBoGuFenQopwEAABhsjEaDpuQmaEpugr57zRgVW2x677BFaw9ZtKdLN8zq5nb9fVup/r6tVBFmk+aOTNbi/FRdNTpVaRTsoAMV6gAGI46kgX4kLtKseaNSNG9UiiTJ6/WqvKHVH7DvP9WkgxXWbhUokm9nuHOepE4x4SHK75iHOT8jVmMyYjU6LUYRof0jFB3sOgP1mPCQy3a1cEJUqB68fpzuuiJP7x6q0u6yRu0+2aCKLlVv0ukT8y/vPtVteVyEWakxYUqJCfPf+u6HKy7CrNAQXwtM362py/3TyyTJ4/W1z3R5vPJ0tNnsvHV7vGpzOFVhl3aXNcrhNqjF4VKr0y17u9t33+GW3eG739TqVGOLU42tTjW1OHy3rU5/m9xLZTT45u3MjI9QRpzvNis+QrlJkRqSGKmshIh+c6EBAKB/MhgM/k4nXyoYKpfbo32nmvRRR/X6rtIGOd2nN3yd2/F/7CiXJGXFR2jWsETNzvNVseckctHlQGQbhC3fAQBAYAaDQaPTYzQ6PUYrF4xQtbVN6wur9d5hiz48Uqt2l0eS1Op0a+0hi79oZ2xGrK4anaIFY1I1JSeeKQIHsW4V6gTqAAYJjqSBfsxgMCgnMVI5iZG6flKmJF/17cn6Fl+4fqpJ+0816cCppjOqbJvbXNp2ol7bTtT7lxkN0tDkKH/QPjItRiNTo5WbGMlOch/TGah/2nbv55KbFKl/v3KY/3FVU5t2lzZoV2mDdpX6LtxwdBxcddXUEVYfqbZd9jGdKUTat61Hf0NUqMl/UUDnhQEZceHKiI9QVny4MuIilBoTxt8GAKBPCTEZNTU3QVNzE3TvopFqcbi062Sjth6v09aSeu0pa5TDfXo7fqqxVf/adUr/2uW7UC4jLtxXwT4sSbPyEpWXHEXAPgDY293++9EE6gAAoIvU2HDdPjNXt8/MVavDrU1Ha/XeIYvWFVq6TXF3qNKqQ5VWPfnBMcWGh2jeqBRdNTpV80el9Mj5KfRdNV0q1Gn5DmCw4EgaGGCMRoPykqOUlxylGzpCdq/Xq7L6Vh2o8IXrhyutOlzZrCpr98pjj1cqqbGrpMauN/dX+peHmowalhLlD9hHpkZrZFq0hiRFyUyY2OtaHC7ZHb6Toim9MJ95ely4lk3I0LIJGZIkh8ujQ5VW7TrZoBN1dtU0t6u6o016dXOb2pxnhu3BZjBIseFmxUeaFR9hVlxkqOIjTj/uGpynRIcrOSZUkbRDBQAMAJGhIZo7MllzRyZLktqcbu0qbdDWknptPV6n3aWN/iokSapsatMreyr0yp4KSb6pg2Z2BOyz8xI1IjWagL0fsrU7/fepUAcAAGcTEWrS1WPTdPXYNHk8Xu0pb9QHRTX6oKha+8qb/OtZ21x6Y1+l3tjnO384MTtOV41O1VWjUzQpO14mI/uLA1nXCnWmAgAwWHAkDQwCBoPBN3dzUqSWd4SiktRgd+hwlS9cL6y06nCVVcUW2xnVxw63R4VVzSqsau62PKQjvB+ZFq0RqTH+oD0vOYo21z2otvn01cHBuAI4NMSoyTnxmpwTf8ZzXq9XtnbXJ0L2dtnbXWp3udXu9Mjh9py+dbnlcHnU7vItkyST0SCT0SCj0SCTQTIZjTIZfcuNBoMM8qq2qkKjhw9VdHioIkJNigw1KSo0xH8/ouOxLzAPVUx4iIwczAEAoHCzSVcMT9YVw30Be7vLrb1lTdpaUqetx+u182SDWp2nq5mrm9u7nSxNigrVjKGJmj40QTOGJmpsZiwXWPYDti4V6lFh7KcDAIDzMxoN/s5Hq64epZrmdm0ortH7RdX6sLimWzfMfeVN2lfepMfXHVFCpFnzR/law88bmaKEqNAgfgr0BEtHkZbJaFAS/38BDBIE6sAglhAV2u2EqiS53B4dr7XrcFWzjlbbdLS6WcUWm07U2uXydJ942uXx6ki1raPFd5V/uclo0JCkyI5q9hiNTPPN6zksJYqq38ugxna6s0Bfa6llMBgUE25WTLhZw1Kie+R3OJ1OvfVWuZYvHyOz2dwjvwMAgMEiLMSkmXmJmpmXqHvl60Sz/1STtnQG7Cfq/Z1xJKnO7tA7B6v0zkHfvl+E2aQpufGaPiRB04cmakpuvGLC2T73NfYuc6jT8h0AAHwaKTFhunVatm6dli2X26PdZY16v7Ba7xfV6HCl1b9eQ4vT3/HIYJAm58RrwehULRidqnGZsRQ8DADVnVNRRofx/xPAoMGRNIBuQkxGX2v3tJhuyx0uj07W2X0BusWmI9W+wL2kxt5tHk5Jcnu8/tbxaw5auj2XEReuvOQoDUuJ0rBkX8g+PCVamfERtIO6QJ3zp0t9L1AHAAD9W2iIUdOGJGjakAStXOC72PJAhdVfwb79eL2au4SzrU63Pj5Wp4+P1UmSjAYpPyPWX8U+fUii0uNoAxlsXQN1Wr4DAIBLFWIyasbQRM0Ymqj/umaMqpratKG4Wu8X1mjT0VrZOvY9vF5pd2mjdpc26pdri5UcHaZ5o5I1f1SKrhyZokSqm/sdl9ujWpvv3GRqLOclAQweHEkDuCChIV2C9gmnl7vcHpXWt3QE7c3+wP1Yja3bfJydKpvaVNnU5j/p2vX9hyZF+kP2YSm+1vHDU6IUH8nOdVfdAvVemEMdAAAMXiGm01O9fG3+cLk9XhVbmrXjRL22n2jQzpMNOtXY6l/f45UOVlh1sMKq5z4+IUnKTojo1iZ+REo0lSy9zEaFOgAA6EHpceH67IxcfXZGrhwuj3aebNAHRdV6v6haxRabf71aW7v+teuU/rXrlAwGaWJ2vK4alaL5zL3eb9TZHfJ2NDFNjeHCWQCDB0fSAC5JiMmoYSnRGpYSraXj0v3L3R6vyhtadMRiU3F1s45V23W81qaSWrsaW5xnvI/D5VGxxdZtJ7tTYlSohnVUtef5q9qjlJsYpdCQwTdnJxXqAAAgWExGg/IzYpWfEasvFgyVJJ1qbNWOE/XacaJB20/Uq8jS7D/JJknlDa0qbzill3efkiTFRZj9LeJnDE3Q+Kw4hZuZ17snUaEOAAB6S2iIUQXDk1QwPEkPLM/XqcZWX7heWKOPjtaq1embTsjrlfaWNWpvWaN+ve6I4iLMunKkr3p9/qgUpcYS1vZFnfOnS1IaFeoABhGOpAH0CN886lEakhSlxWPTuj1Xb3eopMYXrvtaw/vun6yzy+n2nvFe9XaH6u0O7TjZcMbvyEmI6GghH+1vIz88JUopMWEyGAbmVa01NgJ1AADQd2TFRyhrcpZunJwlSWpqdWpXaYN2dgTse8oau3Uuamp1al1htdYVVkuSQk1GTcyO8wfs04Yk0KHoMrO3u/33o8K4eAEAAPSerPgI3TFriO6YNUTtLrd2nmjQhuIabSiuUWFVs3+9plan3thXqTf2VUryTSPUGa5PG5IwKItq+iKL9fR5SSrUAQwmBOoAel1iVKgSoxI1fWhit+Uut0enGltVUmPXsY6Q/XiNXSW1tm47a53cHq9O1LXoRF2L3i+q6fZcdFiI8pKjNDItWqPSYjQ6LUaj0mOUGRfe74N2Wr4DAIC+LC7CrAWjU7VgdKokXyeiAxVN/jbxO07Uq6FLxyKH26MdJxu042SDnt7gWzYqLVrThyZq+hBfm/jshIh+vw8XTJ0t30NNRoWFEKgDAIDgCAsx6YoRybpiRLIeWJ6vyqZWbewI1z88UqvmttNddQ5XWnW40qqnNxxTdFiIrhiepPmjUzRvZIpyEiOD+CkGt+pmKtQBDE4E6gD6jBCT0V/VvmBMarfnbO0uf7heUmPvqG636XitXS0O9xnvZWt3af+pJu0/1dRteXRYiEalRWt0eoxGpsZodHqMRqXF9KtK785A3WDwXZwAAADQl4WGGDU1N0FTcxP01XmS1+vVsRq7dp48HbCfqGvp9prOqYBWby2V5DtZN31oomZ0tIrPz4hljs2LYHf4Tk5TnQ4AAPqSjLgI/9zrLrdHe8oa/dXr+8pPn9Oztbv07iGL3j1kkSQNT4nS/FGpmj86RbPyEpk+qBd1q1AnUAcwiBCoA+gXosNCNCE7ThOy47ot93q9sljbVVJj07GOkL2kxq7jtXaVN7TI84kO8rZ2l3aVNmpXaWO35Zlx4ZqcG68pOQmanBuvCX14Ls9am0OSlBQVqhAT7a4AAED/YjAYNCI1WiNSo/XZGbmSfJUuvhbxDdpxsl4HK6xyd9mRs1jb9ea+Sr3Z0QI0OixEU3LjNX2Ir0385Nx4RYZyeHs2nXOoM386AADoq0JMRl+HoqGJ+uaS0aq1tWvTkVptKK7RxuIa1dkd/nWP1dh1rOa4/vDRcUWYTZozIkkLxvg6JGXGRwTxUwx8NV0q1Gn5DmAw4WgaQL9mMBiUHheu9LhwXTEiudtzbU63jtXYVGxpVlGV77bY0qzyhtYz3qeiqU0V+6v01v4qSVKI0aAxGTG+gD0nXlNy45WXHBX0VqNer9dfoZ5Mu3cAADBApMaEa9mEDC2bkCHJFwDvLWv0B+y7TjbI3qUrka3dpQ+P1OrDI7WSJJPRoPGZsZqZl6g5I5I1My+RgL2Lzvap0QTqAACgn0iODtNNU7J005QseTxeHaywakNxtTYU12hXaaP/4stWp1vvHa7We4erJUlj0mO0YEyqFo5J1ZSceIpRLrOuFeppsQTqAAYPjqYBDFjhZpPGZcZpXGb3qnZbu0tHOsL1oiqbDlU2aV95U7fW8S6PVwdOWXXglFV/2XJSkpQQadbVY9N03cRMFQxPkjkIO+TWVpccbo8k9as29QAAABcjKizEP7+mJLncHhVWNfvmYT/ZoO3H61XdfPpkntvj1d7yJu0tb9IzHx6X2WTQlNwEzR2RrDkjkjUpO27Qnkx1uT1qd/n2H6lQBwAA/ZHRaPB3rrxn4Ug1tTr10dFavV9YrfeLalRrO71fWFjVrMKqZj31wTHFRZg1f1SKFuWn6qrRqYqLMAfxUwwMFquvQt1kNCiJqSgBDCIcTQMYdHwtQhM0JTfBv8zt8epIdbN2lzZqT2mjdpc16Ei1Td4uLeMbWpz6x45y/WNHuRIizbpmfIaum5ihWXmJvXaCtsZ2uq0SgToAABgsQkxGjc+K0/isOH15Tp68Xq/KG1q1/cTpediPVNv86zvdXm07Xq9tx+v1y7XFigkL0axhSZozIklzRyRrRGp00DsP9RZ7++mLRgnUAQDAQBAXYdbyCRlaPiHDX72+vrBa64uqta+80X8+r6nVqdf2Vui1vRUKMRpUMDxJS8am6eqx6UqPo7r60+i8qDUlOkxG4+DYnwYAiUAdACT5rqockx6rMemx+txM31yezW1O7Stv0p6yRu0ubdDmY3X+VqMNLU79fVup/r6tVMnRobpmfLqum5ipGUMTZerBncmulVgE6gAAYLAyGAzKSYxUTmKkbpmaLUlqsDu0uaROm47W6uOjtTpR1+Jfv7ndpfcOW/TeYYskKSs+QlePTdPVY9M0My8xKJ2HeovN4fLfjw4zBXEkAAAAl1/X6vVvLB6pWlu7Piiq0fuF1dpYXKPmdt++kMvj9U8Z9INXD2pSTryWjE3T0nFpGp4SHeRP0T+43B5/N4DUWM5LAhhcCNQB4Cxiws2a09EmVPLNyf5BUbVe31ep9Yer1er0heu1Nof+uqVUf91SqtSYMN0wKVP3LhrZI22karoG6syhDgAA4JcQFeqvVJKksvoWfXysVpuO1unjo7Wqszv8655qbNVzH5/Qcx+fUGx4iBaOSdWScemaNyplwM0zbm8/HahHMa88AAAY4JKjw3TrtGzdOi1bTrdHO040aO0hi9YcrNKpxlb/envLGrW3rFG/WFOkYclRWjQmRTHNkrdru0p0U2d3+Kv/U2Oo8AcwuHA0DQAXKNxs0jXjM3TN+Ay1OFxaX1itN/dVan1htX9eyurmdv1+03GtOVSl335uqiblxF/WMdRQoQ4AAHBBchIj9dnEXH12Rq48Hq8Kq5r10dFabTxSoy0ldXK6fWcDrW0uvbKnQq/sqVCoyag5I5J09dh0LR6bOiBOFNq6BOrR4ZwCAAAAg4fZZFTB8CQVDE/SD67L16FKq949aNG7hyw6XGn1r1dSa1fJJrukEL1UsUk3Tc7SDZMzNSI1JniD74M650+XpDQq1AEMMhxNA8CnEBkaousmZuq6iZmytbu07rBFb+yr1IaiGjncHpXVt+rWpz/Wfy/P151XDL1sc3TW2KhQBwAAuFhGo0FjM2M1NjNWX5k3TNY2pz4oqtHaQxZ9UFjtbwXqcHv0flGN3i+q0X+/Il0xPEmfm5mrJWPTFRrSP9vCd61QH2jV9wAAABfKYDBoXGacxmXG6f6rR6msvkXvdlSu7zhRL09H5XVpfaseX39Uj68/qrEZsbpxcqaun5SpzPiI4H6APqDaevq85EC48BQALgZH0wBwiaLDQnTj5CzdODlLpxpbde/qXdpV2iin26sfvX5IW4/X69FbJyo2/NJbwFOhDgAAcOliw826YVKmbpiUKYfLoy0ldVp7yKK1hyyq6qi88Xqlj47W6aOjdUqKCtWt07J1+8xc5SVHBXn0F6dby3cCdQAAAEm+bkZ3z83T3XPzVGdr15oDFXru/QM6ajX6w/VDlVYdqrTqkbcLNTMvUTdOztTy8RlKiAoN7uCDxNJMhTqAwYujaQC4jLLiI/TC1wr0izVF+t3GEknS2weqdLDCqic+P1UTsuMu6f0J1AEAAC6v0BCj5o1K0bxRKfrxjeO0/1ST1h6y6LW9FTpZ1yLJN1/k/20s0f9tLFHBsCR9blaulo5LU1iIKcijPz9bu9t/n0AdAADgTEnRYbptWraiLPs0/coFWnOoRq/urdDeskb/OtuO12vb8Xo9+OpBLcpP1RdnD9WcEUmXrStlf2DpWqFOoA5gkOFoGgAuM7PJqO8tz9fMoYn65ot71dTqVGl9i1Y89bG+f12+vjh7yKfe2a61OTp+h0FxEZde8Q4AAIDTDAaDJmbHa2J2vO5fPEpbSuq0elup1hys8s+5vrmkTptL6pQQafZXrQ9PiQ7yyM+ue8v3vn8BAAAAQDClxoTp3+bm6d/m5ulErV2v7a3Qq3tO6ViNXZLk8ni15qBFaw5aNCwlSl+cPUQrpmVfls6UfV1Nlwp1Wr4DGGz65yRwANAPLB6bpje/PleTc+Il+ebk/OGrB3XP6t2ytjk/1Xt2VqinRIcNqitgAQAAepvRaNAVI5L1289P1ZYHFul7y8d0a/fe0OLUMx8e16L/3aCVq3epvKEliKM9O1vXlu+hXFMPAABwoYYmR+nri0bqvVXz9ebX5+pr84YptUvHyJIaux56/ZBmPbxOD/xrvw5XWoM42p7XtUI9LZZAHcDgQqAOAD0oOyFS//hagf59bp5/2Zv7K3X9bzbpwKmmi3ovt8erentHoE67dwAAgF6TFB2mr84brvXfnK+/f2W2bpiUqVDT6cPpN/dVatH/btCv1har1eE+xzv1Plu3CnUCdQAAgItlMBg0LjNODyzP10ffXagnPj9Vs/IS/c+3Ot36+7ZSLfv1h7rt6Y/12t4KOVyeII64Z1R3VKibjAYlDdJ55AEMXhxNA0APCw0x6vvXjdXMvER968W9sra5dLKuRZ/5v81au2q+suIjLuh96uzt8vg6jRKoAwAABIHBYFDB8CQVDE9Svd2hf+4s11Mbjqne7lC7y6Nfrzuil3aW64HlY3TthIw+0VGoa8t35lAHAAC4NGaTUddOzNC1EzNUbGnWXzaf1L92lcvecVHl9hMN2n6iQcnRYbp7bp7umjNU4eaBMe1OZ4V6SnSYjMbg7+cCQG+iQh0AesmScel68+tXalJ2nCSpxeHWOweqLvj1ne3eJQJ1AACAYEuMCtVX5g3T+9+6SnfPzVNIx0nFU42tumf1bn32d1t0sOLiOhL1BBuBOgAAQI8YlRajn9w0Xlu+t0g/vnGcRqRG+5+rtbXr0XcKNf8X7+vv20rlcvfvinWX26Nam+/cZGos5yUBDD4E6gDQi3ISI/XwzRP8j/eVN17wa7sG6snR7LgCAAD0BXERZv3gurF6574rNW9Uin/5tuP1uv43m/S9l/erztZ+jnfoWXZavgMAAPSomHCzvlQwVGvvn6e/f2W2lk9IV2cBt8Xargf+tV9LHtuodw5Uyuv1Bnewn1Kd3aHOoafGMH86gMGHQB0Aetno9BiFhvi+fveVX3jVEhXqAAAAfdeI1Bj96a4ZevbO6RqaFClJ8nil1VtLteD/faA/bDout6f3T6Da20/P6R4dTqAOAADQUzqnB3ryjml65755Wpyf5n+upMau//jrLt385MfafKwuiKP8dCzWNv/9NCrUAQxCBOoA0MvMJqPGZsRKko7X2tXU6ryg19V0qWxKoUIdAACgzzEYDFqUn6Y198/TA8vG+CvCrW0u/fiNQ/rG87t7vd1n15bvkQNk/k4AAIC+blRajH5/53S99B8FmjE0wb98T1mjPvfMFt35h219YnqgC1VtPX1ekgp1AIMRgToABEHnPOqStP8Cq9SpUAcAAOgfwkJM+tr84Vr/rfm6bVq2f/kb+yr1Xy/t69VK9c6W71GhJhk7e48CAACgV0wfmqh/fK1Az945XaPTYvzLNxTX6NrHN+kbz+/uVv3dV1maqVAHMLgRqANAEEzMjvff33eq8YJeU2tz+O8TqAMAAPR9qTHh+sVtk/TsndNlNvnC7H/tPqUH/rVPnl4K1f2BOvOnAwAABEVnF6O3vnGl/ve2ScqKj/A/9+qeCi3/9Yf68EhNEEd4fpauFeoE6gAGIQJ1AAiCSTmnK9T3lV1ohfrpK0GTafkOAADQbyzKT9MTn5+qkI4K8X/sKNf3Xz0gr7fnQ/XOlu/RBOoAAABBZTIatGJattZ/a75+cN1YJUSaJUl1doe+9Idt+tXa4l7tZHQxup6XpOU7gMGIQB0AgiAvOVpRob45LPeVN17QazpbvkeFmqgwAgAA6GeWjEvX45+bIlNHqL56a6keev1Qj4bqXq/XH6iz/wgAANA3hIWYdPfcPK3/5lVaMDpFkuT1Sr9ed0R3/mGbam3t53mH3te1Qj0tlkAdwOBDoA4AQWAyGjQ+y1elXtHU1m1+9LPpXId27wAAAP3T8gkZ+uVnJqlzKvPnPj6hh9883GOhepvTo84ip6gwU4/8DgAAAHw6CVGhevbOGfr20tH+/cNNR2t17eMfatvx+uAO7hOqOyrUTUaDkqJCgzwaAOh9BOoAECSTcuL99/efZx71Nqdb1jZfdRGBOgAAQP914+Qs/eLWSTJ0nDT9/abjevSdoh4J1Tur0yVavgMAAPRFRqNBKxeM0OqvzPaf87NY2/W5Z7bo6Q3HemWKoAvRWaGeEh0mY2f6DwCDCIE6AATJxOzT86jvPc886l1bPTF/OgAAQP+2Ylq2fnbLBP/jpzcc06/eO3LZf4+9S6BOy3cAAIC+a/awJL319St1xfAkSZLb49XP3i7UV/68Q00tzqCOzeX2qK7j3GRqLOclAQxOBOoAECQTs+L99883j3rXlvBUqAMAAPR/n52Rq5/eNN7/+PF1R/SbdZc3VLcRqAMAAPQbKTFh+svds3TvwhH+Ze8drta1v/lQe8sagzauOrvDP41QagzzpwMYnAjUASBIchIjlBBpliTtK286ZwunboE6FeoAAAADwhdmD9GD14/1P/7ftcV6esOxy/b+XSvUYwjUAQAA+jyT0aBvLhmt5+6a4T9vWN7Qqluf/ljrDluCMiaLtc1/P40KdQCDFIE6AASJwWDQhOx4Sb4rPSua2s66bo2NCnUAAICB6K45efrv5fn+xz97u1DvF1Vflve2O6hQBwAA6I+uGp2qN79+pabmxkuSnG6vvvH8Hh2tbu71sVRbT5+XpEIdwGBFoA4AQTSpyzzq+87Ruqm22eG/T6AOAAAwsHxl3jB9a8ko/+PnPjpxWd7X1u723ydQBwAA6F8y4yP0wtcKdO2EDEm+6Xy+8uedamrt3TnVLc1UqAMAgToABNHEjgp1Sdpb3nTW9Wpsp3dcCdQBAAAGnv/vqhHKToiQJG08UqOKxtZLfs+uLd+jw0yX/H4AAADoXWaTUb+4baLyM2IlScdr7brv+d1ye84+deTl1q1CnUAdwCBFoA4AQTSxa4V6eeNZ1+s2hzqBOgAAwIBjNBr0mek5kiSvV3pxR/klv2fXQJ0KdQAAgP4pMjREv/viNP+c6u8X1eiXa4t67fdXd6lQp+U7gMGKQB0AgigtNtzfKmn/qSZ5znJ1addAPSmKQB0AAGAgunVatgwG3/0Xd5addd/wQjW3EagDAAAMBDmJkfrt56fKZPTtLD7x/jG9ua+yV363pUuFelosgTqAwYlAHQCCrLPte3ObSyfq7AHXqbH5dlzjI80KDeGrGwAAYCDKjI/QvJEpkqTyhlZ9fKzukt6ve8t3AnUAAID+bM6IZH1veb7/8bde3KvDldYe/72dFeomo0FJUaE9/vsAoC8ilQGAIJvUre37mfOoe71ef4V6SjTV6QAAAAPZ7TNy/Pef3156Se9ld3SpUA8lUAcAAOjv/m3OUN0yNUuS1Op066t/2aEGu6NHf2dnhXpKdJiMHRXyADDYEKgDQJBN6KhQl6S9AeZRt7W71Ob0SGL+dAAAgIFuUX6av/Ln3YOWSzpBamt3++9ToQ4AAND/GQwG/c/NEzSxo0CnrL5V9/59t1xuT4/8Ppfbo7qOzpmpsZyXBDB4EagDQJBNzDp3hXrX+dMJ1AEAAAa20BCjbp7iqzpyuD16Zc+pT/1e3Vq+hxOoAwAADAThZpOe/sI0JUf7LsLcdLRWj75T2CO/q87ukMfru58aw/zpAAYvAnUACLKEqFDlJkZKkg5WNJ1xRWm3QJ2W7wAAAAPeZ7u0fX9he5m8Xu+neh9bl0A9Ksx0yeMCAABA35AZH6En75imkI4W7M98eFyv7P70F2KejcXa5r+fRoU6gEGMQB0A+oDONk1tTo+OVNu6PVdrO93mkwp1AACAgW9kWoym5sZLkgqrmgN2MboQnRXqZpNBYSEE6gAAAAPJzLxEPXjDOP/j7/xznw6c+nT7jWdTbT1d6EOFOoDBjEAdAPqAzkBdkvZ9Yh71mubTV4ISqAMAAAwO3arUd5R9qvfoDNSjmD8dAABgQPrCrFzd3rHf2O7y6Kt/3qEGu+M8r7pwlmYq1AFAIlAHgD5hYna8//7eT1Qg1diYQx0AAGCwuW5ipqJCfVXlr+2pUIvDdZ5XnMnW7pYkRYUSqAMAAAxEBoNBD904zt/dqKKpTX/8+MRle/9uFeoE6gAGMQJ1AOgDxmfFyeCb8kj7PxmoNxOoAwAADDZRYSG6bmKmJN9c6G/tr7ro9+isUI+mQh0AAGDACgsx6fHPTZGpYz7157eVyun2XJb3ru5SoU7LdwCDGYE6APQB0WEhGpESLUkqrLKq3eX2P9c1UE+OJlAHAAAYLD7Tpe37P7ZfXNt3l9ujVmdHhXoY86cDAAAMZNkJkVqcnypJqm5u19pDlsvyvl0r1NNiCdQBDF4E6gDQR0zomEfd6fbqcGWzf3lny3eT0aCEyNCgjA0AAAC9b2puvEak+i663HaiXsdqbBf8Wrvj9AWazKEOAAAw8H1x9lD//b9sPnlZ3rNzDnWT0aCkKM5LAhi8CNQBoI+Y1GUe9X3ljf77nRXqSVGh/tZNAAAAGPgMBoNu71qlvuPCq9Q7271LtHwHAAAYDK4YnqRhyVGSpM0ldTpa3XyeV5yfpaNCPSU6TEbOSwIYxAjUAaCPmNhRoS5J+zrmUfd4vKq1OSQxfzoAAMBgdPOULJlNvpOX/9x56oLnw+waqFOhDgAAMPAZjQbdMXuI//Fft5Re0vu53B7VdXTOTI3lvCSAwS2ogfpTTz2liRMnKjY2VrGxsSooKNDbb7/tf76trU0rV65UUlKSoqOjtWLFClksl2fuDwDoa/IzYhXScaVnZ4V6Q4tDbo9XEoE6AADAYJQUHabF+WmSpFpbu9YXVl/Q62xUqAMAAAw6t07NVrjZF/v8c2d5t4ssL1ad3aGO05JKjWH+dACDW1AD9ezsbP3sZz/Tzp07tWPHDi1cuFA33nijDh48KEm6//779frrr+vFF1/Uhg0bVFFRoVtuuSWYQwaAHhNuNmlMRowk6Wi1TfZ2l786XfK1VgIAAMDg89mubd+3X1jbd3v76TnUCdQBAAAGh7hIs26YlClJam536bW9FZ/6vao72r1LUhoV6gAGuaAeVV9//fXdHj/88MN66qmntGXLFmVnZ+vZZ5/V6tWrtXDhQknSH//4R+Xn52vLli2aPXt2wPdsb29Xe/vpL3qr1SpJcjqdcjqdPfRJAODyGJcRqwOnrPJ4pT2ldXK6vf7nkqLMfI9J/v8G/LcAAHwS2wgMVLOHxisjLlyVTW16v6haZXXNSo89d5VQU0ub/354iIG/Cwx6bCMAAOcykLYTn5uerX/sKJck/fnjE1oxOV0Gw8XPf36qwea/z3lJAAPRxXyv9ZnL1N1ut1588UXZ7XYVFBRo586dcjqdWrx4sX+dMWPGKDc3V5s3bz5roP7II4/ooYceOmP5u+++q8jIyB4bPwBcDoZ6gySTJOnF97Yq2iz/4+rSo3rrrSNBG1tfs3bt2mAPAQDQR7GNwEA0MdqoyiajPF7pZ8+/ryXZ3nOuv63m9H7liaOFest2uBdGCfR9bCMAAOcyULYTQ6JNOmkz6HBVs578x9vKi7n49/jIcnp/svpEsd56q+jyDhIAgqylpeWC1w16oL5//34VFBSora1N0dHRevnllzV27Fjt2bNHoaGhio+P77Z+Wlqaqqqqzvp+DzzwgFatWuV/bLValZOToyVLlig2NranPgYAXBZ5lc16/snNkiRnbJayMmOlo8WSpPkzp2j5hPRgDq9PcDqdWrt2ra6++mqZzeZgDwcA0IewjcBANrGhVe/+6kN5vdJ+e7R+ec1cGY1nrzSq31oqHS2UJM2aOknLJ2f21lCBPoltBADgXAbadqIt45S+8y/f1LrHTTlauXzCRb/H0fVHpZISSdLCK6ZrweiUyzpGAAi2zi7nFyLogfro0aO1Z88eNTU16aWXXtKdd96pDRs2fOr3CwsLU1jYmfN5mM3mAbEhBDCwjc2KV1iIUe0ujw5UWJUZH+F/Li0+ku+xLvheBwCcDdsIDER5qWbNGZ6sTUdrVVrfql3lzSoYnnTW9VtdpyvYYyPD+JsAOrCNAACcy0DZTtw4JUePvFOsxhan3j5g0Q+vH6ek6IubB73WfroVcmZC1ID47wIAXV3M95qxB8dxQUJDQzVixAhNmzZNjzzyiCZNmqRf//rXSk9Pl8PhUGNjY7f1LRaL0tOp0AQwMIWYjBqX6eumcbKuRUeqT89VlBJzcTu9AAAAGFg+MyPHf/+F7aXnXNfW5vLfjw4L+rX0AAAA6EXhZpM+M9237+hwe/xzql+Mamu7/35abPhlGxsA9Ed97qja4/Govb1d06ZNk9ls1rp167RixQpJUlFRkUpLS1VQUBDkUQJAz5mYHa9dpY2SpM3H6vzLCdQBAAAGtyVj0xQXYVZTq1Ov76vU/lNNSogMVXxkqBIizUqIClV8pFkJkaE6VHm6dV0UgToAAMCg8/mZufrdRl/L9tXbTuqr84bJdI4pgz7J0twmSTIZDUqKCu2RMQJAfxHUo+oHHnhAy5YtU25urpqbm7V69Wp98MEHWrNmjeLi4nT33Xdr1apVSkxMVGxsrO69914VFBRo9uzZwRw2APSoSTlx/vvtLo8kKSzEqBhOhAIAAAxq4WaTbp6Spec+PiG3x6tjNXZJ9vO+LjrM1PODAwAAQJ8yNDlK80alaGNxjcrqW7WxuEYLxqRe8Os7K9RTosNkvIggHgAGoqCmM9XV1frSl76kyspKxcXFaeLEiVqzZo2uvvpqSdKvfvUrGY1GrVixQu3t7Vq6dKmefPLJYA4ZAHrcxOz4M5alxITJYGDHFQAAYLC7d+EIlda3qNjSrMYWp2ztrnOuHxMeooy4iF4aHQAAAPqSL84eoo3FNZKkv2w5ecGBusvtUa3NF6inxtI1EwCCGqg/++yz53w+PDxcTzzxhJ544oleGhEABF9eUpRiwkLU3OXkKO3eAQAAIElJ0WH6w5dn+B87XB41tjrU2OJUg92hhhanGlt8t/Z2lxbmp9LyHQAAYJBaOCZVWfEROtXYqveLqlVW36KcxMjzvq7O7pDH67ufGsP86QDAUTUA9DFGo0Hjs+K0uaTL/OnRBOoAAAA4U2iIUakx4ZzoBAAAwBlMRoM+PytXv1hTJK9X+tvWUn132Zjzvq6z3bskpVGhDgAyBnsAAIAzTewyj7pEhToAAAAAAACAi/eZ6Tkym3xTSf5jR5naXe7zvsZibfPf58JNACBQB4A+adIn5lFPpkIdAAAAAAAAwEVKiQnTNeMzJEn1dofe3l913tdUN1OhDgBdEagDQB80MZsKdQAAAAAAAACX7ouzh/jv/2XLyfOu361CnUAdAAjUAaAvyoqPUGJUqP8xgToAAAAAAACAT2PG0ASNTouRJO082aCDFU1nrNPc5tSu0ga9sL1UHxRV+5fT8h0ApJBgDwAAcCaDwaCJ2XH6oKhGEoE6AAAAAAAAgE/HYDDoCwVD9INXDkiSnnz/mK4anaIj1TYVW5pVXNWsiqa2gK9NjyNQBwACdQDoo26anKUPimqUFhumsRmxwR4OAAAAAAAAgH7q5ilZ+tlbh2V3uPXm/kq9ub/yvK+5ZUqWkqMp9AEAAnUA6KNumpKlCdlxSo0JU7jZFOzhAAAAAAAAAOinosNCtGJatv68+cw51GPCQzQqLUaj0qI1MjXGfz81lup0AJAI1AGgTxueEh3sIQAAAAAAAAAYAL69dLQMklqdbo1Ki9HItBiNTotRWmyYDAZDsIcHAH0WgToAAAAAAAAAAMAAFxNu1kM3jg/2MACg3zEGewAAAAAAAAAAAAAAAPRFBOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAEEN1B955BHNmDFDMTExSk1N1U033aSioqJu67S1tWnlypVKSkpSdHS0VqxYIYvFEqQRAwAAAAAAAAAAAAAGi6AG6hs2bNDKlSu1ZcsWrV27Vk6nU0uWLJHdbvevc//99+v111/Xiy++qA0bNqiiokK33HJLEEcNAAAAAAAAAAAAABgMQoL5y995551uj5977jmlpqZq586dmjdvnpqamvTss89q9erVWrhwoSTpj3/8o/Lz87VlyxbNnj07GMMGAAAAAAAAAAAAAAwCQQ3UP6mpqUmSlJiYKEnauXOnnE6nFi9e7F9nzJgxys3N1ebNmwMG6u3t7Wpvb/c/tlqtkiSn0ymn09mTwwcA9ILO73K+0wEAn8Q2AgBwNmwjAADnwnYCAAafi/nO7zOBusfj0X333ac5c+Zo/PjxkqSqqiqFhoYqPj6+27ppaWmqqqoK+D6PPPKIHnrooTOWv/vuu4qMjLzs4wYABMfatWuDPQQAQB/FNgIAcDZsIwAA58J2AgAGj5aWlgtet88E6itXrtSBAwe0adOmS3qfBx54QKtWrfI/tlqtysnJ0ZIlSxQbG3upwwQABJnT6dTatWt19dVXy2w2B3s4AIA+hG0EAOBs2EYAAM6F7QQADD6dXc4vRJ8I1O+55x698cYb2rhxo7Kzs/3L09PT5XA41NjY2K1K3WKxKD09PeB7hYWFKSws7IzlZrOZDSEADCB8rwMAzoZtBADgbNhGAADOhe0EAAweF/N9b+zBcZyX1+vVPffco5dfflnr169XXl5et+enTZsms9msdevW+ZcVFRWptLRUBQUFvT1cAAAAAAAAAAAAAMAgEtQK9ZUrV2r16tV69dVXFRMT458XPS4uThEREYqLi9Pdd9+tVatWKTExUbGxsbr33ntVUFCg2bNnB3PoAAAAAAAAAAAAAIABLqiB+lNPPSVJuuqqq7ot/+Mf/6gvf/nLkqRf/epXMhqNWrFihdrb27V06VI9+eSTvTxSAAAAAAAAAAAAAMBgE9RA3ev1nned8PBwPfHEE3riiSd6YUQAAAAAAAAAAAAAAPgEdQ51AAAAAAAAAAAAAAD6KgJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIgUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AAAAAAAAAAAAAAACIFAHAAAAAAAAAAAAACAAAnUAAAAAAAAAAAAAAAIIaqC+ceNGXX/99crMzJTBYNArr7zS7Xmv16sf/vCHysjIUEREhBYvXqwjR44EZ7AAAAAAAAAAAAAAgEElqIG63W7XpEmT9MQTTwR8/uc//7kef/xxPf3009q6dauioqK0dOlStbW19fJIAQAAAAAAAAAAAACDTUgwf/myZcu0bNmygM95vV499thj+v73v68bb7xRkvTnP/9ZaWlpeuWVV3T77bcHfF17e7va29v9j61WqyTJ6XTK6XRe5k8AAOhtnd/lfKcDAD6JbQQA4GzYRgAAzoXtBAAMPhfznR/UQP1cjh8/rqqqKi1evNi/LC4uTrNmzdLmzZvPGqg/8sgjeuihh85Y/u677yoyMrLHxgsA6F1r164N9hAAAH0U2wgAwNmwjQAAnAvbCQAYPFpaWi543T4bqFdVVUmS0tLSui1PS0vzPxfIAw88oFWrVvkfW61W5eTkaMmSJYqNje2ZwQIAeo3T6dTatWt19dVXy2w2B3s4AIA+hG0EAOBs2EYAAM6F7QQADD6dXc4vRJ8N1D+tsLAwhYWFnbHcbDazIQSAAYTvdQDA2bCNAACcDdsIAMC5sJ0AgMHjYr7vjT04jkuSnp4uSbJYLN2WWywW/3MAAAAAAAAAAAAAAPSUPhuo5+XlKT09XevWrfMvs1qt2rp1qwoKCoI4MgAAAAAAAAAAAADAYBDUlu82m01Hjx71Pz5+/Lj27NmjxMRE5ebm6r777tNPf/pTjRw5Unl5efrBD36gzMxM3XTTTcEbNAAAAAAAAAAAAABgUAhqoL5jxw4tWLDA/3jVqlWSpDvvvFPPPfec/uu//kt2u11f/epX1djYqLlz5+qdd95ReHh4sIYMAAAAAAAAAAAAABgkghqoX3XVVfJ6vWd93mAw6Mc//rF+/OMf9+KoAAAAAAAAAAAAAADow3OoAwAAAAAAAAAAAAAQTATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAABOoAAAAAAAAAAAAAAARAoA4AAAAAAAAAAAAAQAAE6gAAAAAAAAAAAAAABECgDgAAAAAAAAAAAABAAATqAAAAAAAAAAAAAAAEQKAOAAAAAAAAAAAAAEAA/SJQf+KJJzR06FCFh4dr1qxZ2rZtW7CHBAAAAAAAAAAAAAAY4Pp8oP7CCy9o1apVevDBB7Vr1y5NmjRJS5cuVXV1dbCHBgAAAAAAAAAAAAAYwPp8oP7LX/5SX/nKV3TXXXdp7NixevrppxUZGak//OEPwR4aAAAAAAAAAAAAAGAACwn2AM7F4XBo586deuCBB/zLjEajFi9erM2bNwd8TXt7u9rb2/2Pm5qaJEn19fVyOp09O2AAQI9zOp1qaWlRXV2dzGZzsIcDAOhD2EYAAM6GbQQA4FzYTgDA4NPc3CxJ8nq95123TwfqtbW1crvdSktL67Y8LS1NhYWFAV/zyCOP6KGHHjpjeV5eXo+MEQAAAAAAAAAAAADQ/zQ3NysuLu6c6/TpQP3TeOCBB7Rq1Sr/Y4/Ho/r6eiUlJclgMARxZACAy8FqtSonJ0dlZWWKjY0N9nAAAH0I2wgAwNmwjQAAnAvbCQAYfLxer5qbm5WZmXnedft0oJ6cnCyTySSLxdJtucViUXp6esDXhIWFKSwsrNuy+Pj4nhoiACBIYmNjOcABAATENgIAcDZsIwAA58J2AgAGl/NVpncy9vA4LkloaKimTZumdevW+Zd5PB6tW7dOBQUFQRwZAAAAAAAAAAAAAGCg69MV6pK0atUq3XnnnZo+fbpmzpypxx57THa7XXfddVewhwYAAAAAAAAAAAAAGMD6fKD+2c9+VjU1NfrhD3+oqqoqTZ48We+8847S0tKCPTQAQBCEhYXpwQcfPGN6DwAA2EYAAM6GbQQA4FzYTgAAzsXg9Xq9wR4EAAAAAAAAAAAAAAB9TZ+eQx0AAAAAAAAAAAAAgGAhUAcAAAAAAAAAAAAAIAACdQAAAAAAAAAAAAAAAiBQBwAAAAAAAAAAAAAgAAJ1AEDQPfLII5oxY4ZiYmKUmpqqm266SUVFRd3WaWtr08qVK5WUlKTo6GitWLFCFoul2zqlpaW69tprFRkZqdTUVH3729+Wy+XqzY8CAOhhP/vZz2QwGHTffff5l7GNAIDB69SpU/rCF76gpKQkRUREaMKECdqxY4f/ea/Xqx/+8IfKyMhQRESEFi9erCNHjnR7j/r6et1xxx2KjY1VfHy87r77btlstt7+KACAy8ztdusHP/iB8vLyFBERoeHDh+snP/mJvF6vfx22EwCAC0GgDgAIug0bNmjlypXasmWL1q5dK6fTqSVLlshut/vXuf/++/X666/rxRdf1IYNG1RRUaFbbrnF/7zb7da1114rh8Ohjz/+WH/605/03HPP6Yc//GEwPhIAoAds375d//d//6eJEyd2W842AgAGp4aGBs2ZM0dms1lvv/22Dh06pP/93/9VQkKCf52f//znevzxx/X0009r69atioqK0tKlS9XW1uZf54477tDBgwe1du1avfHGG9q4caO++tWvBuMjAQAuo0cffVRPPfWUfvvb3+rw4cN69NFH9fOf/1y/+c1v/OuwnQAAXAiDt+vlWAAA9AE1NTVKTU3Vhg0bNG/ePDU1NSklJUWrV6/WrbfeKkkqLCxUfn6+Nm/erNmzZ+vtt9/Wddddp4qKCqWlpUmSnn76aX3nO99RTU2NQkNDg/mRAACXyGazaerUqXryySf105/+VJMnT9Zjjz3GNgIABrHvfve7+uijj/Thhx8GfN7r9SozM1Pf/OY39a1vfUuS1NTUpLS0ND333HO6/fbbdfjwYY0dO1bbt2/X9OnTJUnvvPOOli9frvLycmVmZvba5wEAXF7XXXed0tLS9Oyzz/qXrVixQhEREfrrX//KdgIAcMGoUAcA9DlNTU2SpMTEREnSzp075XQ6tXjxYv86Y8aMUW5urjZv3ixJ2rx5syZMmOAPSiRp6dKlslqtOnjwYC+OHgDQE1auXKlrr72227ZAYhsBAIPZa6+9punTp+u2225TamqqpkyZomeeecb//PHjx1VVVdVtGxEXF6dZs2Z120bEx8f7QxJJWrx4sYxGo7Zu3dp7HwYAcNldccUVWrdunYqLiyVJe/fu1aZNm7Rs2TJJbCcAABcuJNgDAACgK4/Ho/vuu09z5szR+PHjJUlVVVUKDQ1VfHx8t3XT0tJUVVXlX6drUNL5fOdzAID+6/nnn9euXbu0ffv2M55jGwEAg1dJSYmeeuoprVq1St/73ve0fft2ff3rX1doaKjuvPNO/3d8oG1A121Eampqt+dDQkKUmJjINgIA+rnvfve7slqtGjNmjEwmk9xutx5++GHdcccdksR2AgBwwQjUAQB9ysqVK3XgwAFt2rQp2EMBAPQBZWVl+sY3vqG1a9cqPDw82MMBAPQhHo9H06dP1//8z/9IkqZMmaIDBw7o6aef1p133hnk0QEAgu0f//iH/va3v2n16tUaN26c9uzZo/vuu0+ZmZlsJwAAF4WW7wCAPuOee+7RG2+8offff1/Z2dn+5enp6XI4HGpsbOy2vsViUXp6un8di8VyxvOdzwEA+qedO3equrpaU6dOVUhIiEJCQrRhwwY9/vjjCgkJUVpaGtsIABikMjIyNHbs2G7L8vPzVVpaKun0d3ygbUDXbUR1dXW3510ul+rr69lGAEA/9+1vf1vf/e53dfvtt2vChAn64he/qPvvv1+PPPKIJLYTAIALR6AOAAg6r9ere+65Ry+//LLWr1+vvLy8bs9PmzZNZrNZ69at8y8rKipSaWmpCgoKJEkFBQXav39/t4OctWvXKjY29oyTbACA/mPRokXav3+/9uzZ4/+ZPn267rjjDv99thEAMDjNmTNHRUVF3ZYVFxdryJAhkqS8vDylp6d320ZYrVZt3bq12zaisbFRO3fu9K+zfv16eTwezZo1qxc+BQCgp7S0tMho7B6BmEwmeTweSWwnAAAXjpbvAICgW7lypVavXq1XX31VMTEx/jmo4uLiFBERobi4ON19991atWqVEhMTFRsbq3vvvVcFBQWaPXu2JGnJkiUaO3asvvjFL+rnP/+5qqqq9P3vf18rV65UWFhYMD8eAOASxMTEaPz48d2WRUVFKSkpyb+cbQQADE7333+/rrjiCv3P//yPPvOZz2jbtm363e9+p9/97neSJIPBoPvuu08//elPNXLkSOXl5ekHP/iBMjMzddNNN0nyVbRfc801+spXvqKnn35aTqdT99xzj26//XZlZmYG8dMBAC7V9ddfr4cffli5ubkaN26cdu/+/9u705Covj+O459pyjCvYYvlSEaLbYRgC0XRMi2kJlFUFhbZpAlthLQgQRsh0V5EZEVlRSUVkaHQRjklPqgki4TAkooCK8KENkfR+T8IL01zS6sp4/97v+DC3HPPnPOd+2QYPnPOLdWuXbuUmpoqie8JAEDz2bxer7eliwAA/LfZbDbL9pycHLlcLklSTU2NVq5cqdzcXHk8HsXFxWn//v0+22s9f/5cixcvltvtVkhIiObPn68tW7aodWv+PwYA/0+cTqdiY2O1Z88eSXxHAMB/WUFBgdasWaPHjx+rZ8+eWrFihdLT083rXq9XGzZs0KFDh1RdXa1Ro0Zp//796tu3r9mnqqpKy5YtU35+vlq1aqUZM2Zo7969MgyjJT4SACBA3r9/r3Xr1unChQt68+aNIiMjlZycrPXr1ysoKEgS3xMAgOYhUAcAAAAAAAAAAAAAwALPUAcAAAAAAAAAAAAAwAKBOgAAAAAAAAAAAAAAFgjUAQAAAAAAAAAAAACwQKAOAAAAAAAAAAAAAIAFAnUAAAAAAAAAAAAAACwQqAMAAAAAAAAAAAAAYIFAHQAAAAAAAAAAAAAACwTqAAAAAAAAAAAAAABYIFAHAAAAAAAB43K5ZLPZZLPZlJeXF9Cx3W63Ofa0adMCOjYAAAAAAFYI1AEAAAAA+IGvA+KvjydPnrR0af+s+Ph4VVZWKiEhwWz7XsDucrmaHY6PHDlSlZWVmjVrVoAqBQAAAADgx1q3dAEAAAAAAPzr4uPjlZOT49MWHh7u16+2tlZBQUF/q6x/Vtu2bRURERHwcYOCghQREaHg4GB5PJ6Ajw8AAAAAwLdYoQ4AAAAAQBMaA+KvD7vdLqfTqWXLlikjI0OdO3dWXFycJKmsrEwJCQkyDENdu3bVvHnz9PbtW3O8jx8/KiUlRYZhyOFwaOfOnXI6ncrIyDD7WK3oDgsL07Fjx8zzFy9eaNasWQoLC1PHjh01depUPXv2zLzeuPp7x44dcjgc6tSpk5YuXaq6ujqzj8fjUWZmpqKiotS2bVtFR0fryJEj8nq9io6O1o4dO3xquH///h9bof/s2TPL3QCcTmfA5wIAAAAAoDkI1AEAAAAA+A3Hjx9XUFCQiouLdeDAAVVXV2v8+PEaNGiQSkpKdPnyZb1+/dpnm/LVq1fr5s2bunjxoq5evSq326179+791Lx1dXWKi4tTaGioioqKVFxcLMMwFB8fr9raWrNfYWGhKioqVFhYqOPHj+vYsWM+oXxKSopyc3O1d+9ePXr0SAcPHpRhGLLZbEpNTfVbmZ+Tk6MxY8YoOjr6127YD0RFRamystI8SktL1alTJ40ZMybgcwEAAAAA0Bxs+Q4AAAAAQBMKCgpkGIZ5npCQoHPnzkmS+vTpo23btpnXsrKyNGjQIG3evNlsO3r0qKKiolReXq7IyEgdOXJEJ0+e1IQJEyR9CeW7dev2UzWdOXNGDQ0NOnz4sGw2m6QvYXdYWJjcbrcmTZokSerQoYP27dsnu92u/v37KzExUdevX1d6errKy8t19uxZXbt2TRMnTpQk9erVy5zD5XJp/fr1unPnjoYNG6a6ujqdPn3ab9V6cyUnJ8tut/u0eTweJSYmSpLsdru5VXxNTY2mTZumESNGaOPGjb80HwAAAAAAv4tAHQAAAACAJowbN07Z2dnmeUhIiPl6yJAhPn0fPHigwsJCnwC+UUVFhT5//qza2loNHz7cbO/YsaP69ev3UzU9ePBAT548UWhoqE97TU2NKioqzPOBAwf6hNgOh0MPHz6U9GX7drvdrrFjx1rOERkZqcTERB09elTDhg1Tfn6+PB6PkpKSfqrWRrt37zaD+0aZmZmqr6/365uamqr379/r2rVratWKDfYAAAAAAC2DQB0AAAAAgCaEhIR8d4vzr8N1Sfrw4YOmTJmirVu3+vV1OBzNfva4zWaT1+v1afv62ecfPnzQkCFDdOrUKb/3hoeHm6/btGnjN25DQ4MkKTg4uMk6Fi5cqHnz5mn37t3KycnR7Nmz1a5du2Z9hm9FRET43cfQ0FBVV1f7tGVlZenKlSu6c+eO3x8GAAAAAAD4mwjUAQAAAAAIoMGDB+v8+fPq0aOHWrf2/9ndu3dvtWnTRrdv31b37t0lSe/evVN5ebnPSvHw8HBVVlaa548fP9anT5985jlz5oy6dOmi9u3b/1KtMTExamho0M2bN/1WjjeaPHmyQkJClJ2drcuXL+vWrVu/NFdznT9/Xps2bdKlS5fUu3fvPzoXAAAAAABNYc80AAAAAAACaOnSpaqqqlJycrLu3r2riooKXblyRQsWLFB9fb0Mw1BaWppWr16tGzduqKysTC6Xy29b8/Hjx2vfvn0qLS1VSUmJFi1a5LPafO7cuercubOmTp2qoqIiPX36VG63W8uXL9fLly+bVWuPHj00f/58paamKi8vzxzj7NmzZh+73S6Xy6U1a9aoT58+GjFiRGBulIWysjKlpKQoMzNTAwcO1KtXr/Tq1StVVVX9sTkBAAAAAPgRAnUAAAAAAAIoMjJSxcXFqq+v16RJkxQTE6OMjAyFhYWZofn27ds1evRoTZkyRRMnTtSoUaP8nsW+c+dORUVFafTo0ZozZ45WrVrls9V6u3btdOvWLXXv3l3Tp0/XgAEDlJaWppqamp9asZ6dna2ZM2dqyZIl6t+/v9LT0/Xx40efPmlpaaqtrdWCBQt+4840raSkRJ8+fVJWVpYcDod5TJ8+/Y/OCwAAAADA99i83z6QDQAAAAAA/HVOp1OxsbHas2dPS5fip6ioSBMmTNCLFy/UtWvXH/Z1uVyqrq5WXl7eH6vnb8wBAAAAAIDECnUAAAAAAPAdHo9HL1++1MaNG5WUlNRkmN6ooKBAhmGooKAgoPUUFRXJMAydOnUqoOMCAAAAAPA9rVu6AAAAAAAA8G/Kzc1VWlqaYmNjdeLEiWa9Z9u2bVq7dq0kyeFwBLSeoUOH6v79+5IkwzACOjYAAAAAAFbY8h0AAAAAAAAAAAAAAAts+Q4AAAAAAAAAAAAAgAUCdQAAAAAAAAAAAAAALBCoAwAAAAAAAAAAAABggUAdAAAAAAAAAAAAAAALBOoAAAAAAAAAAAAAAFggUAcAAAAAAAAAAAAAwAKBOgAAAAAAAAAAAAAAFgjUAQAAAAAAAAAAAACw8D83rQHH5VkF4AAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "if domain.comm.rank == 0:\n", + " import matplotlib.pyplot as plt\n", + "\n", + " fig = plt.figure(figsize=(25, 8))\n", + " plt.plot(freq, 20 * np.log10(np.abs(p_mic) / np.sqrt(2) / 2e-5), linewidth=2)\n", + " plt.grid(True)\n", + " plt.xlabel(\"Frequency [Hz]\")\n", + " plt.ylabel(\"SPL [dB]\")\n", + " plt.xlim([freq[0], freq[-1]])\n", + " plt.ylim([0, 90])\n", + " plt.legend()\n", + " plt.show()" + ] + } + ], + "metadata": { + "jupytext": { + "formats": "ipynb,py:light" + }, + "kernelspec": { + "display_name": "Python 3 (DOLFINx complex)", + "language": "python", + "name": "python3-complex" + }, + "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.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/chapter2/helmholtz_code.py b/chapter2/helmholtz_code.py new file mode 100644 index 00000000..09694ebc --- /dev/null +++ b/chapter2/helmholtz_code.py @@ -0,0 +1,290 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,py:light +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.16.6 +# kernelspec: +# display_name: Python 3 (DOLFINx complex) +# language: python +# name: python3-complex +# --- + +# # Implementation +# Author: Antonio Baiano Svizzero and Jørgen S. Dokken +# +# In this tutorial, you will learn how to: +# - Define acoustic velocity and impedance boundary conditions +# - Compute acoustic sound pressure for multiple frequencies +# - Compute the Sound Pressure Level (SPL) at a given microphone position +# +# ## Test problem +# As an example, we will model a plane wave propagating in a tube. +# While it is a basic test case, the code can be adapted to way more complex problems where velocity and impedance boundary conditions are needed. +# We will apply a velocity boundary condition $v_n = 0.001$ to one end of the tube and an impedance $Z$ computed with the Delaney-Bazley model, +# supposing that a layer of thickness $d = 0.02$ and flow resistivity $\sigma = 1e4$ is placed at the second end of the tube. +# The choice of such impedance (the one of a plane wave propagating in free field) will give, as a result, a solution with no reflections. +# +# First, we create the mesh with gmsh, also setting the physical group for velocity and impedance boundary conditions and the respective tags. + +# + +import gmsh + +gmsh.initialize() + +# meshsize settings +meshsize = 0.02 +gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) +gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize) + + +# create geometry +L = 1 +W = 0.1 + +gmsh.model.occ.addBox(0, 0, 0, L, W, W) +gmsh.model.occ.synchronize() + +# setup physical groups +v_bc_tag = 2 +Z_bc_tag = 3 +gmsh.model.addPhysicalGroup(3, [1], 1, "air_volume") +gmsh.model.addPhysicalGroup(2, [1], v_bc_tag, "velocity_BC") +gmsh.model.addPhysicalGroup(2, [2], Z_bc_tag, "impedance") + +# mesh generation +gmsh.model.mesh.generate(3) +# - + +# Then we import the gmsh mesh with the ```dolfinx.io.gmshio``` function. + +# + +from mpi4py import MPI +from dolfinx import fem, io, default_scalar_type, geometry +from dolfinx.fem.petsc import LinearProblem +import ufl +import numpy as np +import numpy.typing as npt + +mesh_data = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3) +domain = mesh_data.mesh +assert mesh_data.facet_tags is not None +facet_tags = mesh_data.facet_tags +# - + +# We define the function space for our unknown $p$ and define the range of frequencies we want to solve the Helmholtz equation for. + +# + + +V = fem.functionspace(domain, ("Lagrange", 1)) + +# Discrete frequency range +freq = np.arange(10, 1000, 5) # Hz + +# Air parameters +rho0 = 1.225 # kg/m^3 +c = 340 # m/s + +# - + +# ## Boundary conditions +# +# The Delaney-Bazley model is used to compute the characteristic impedance and wavenumber of the porous layer, +# treated as an equivalent fluid with complex valued properties +# +# $$ +# \begin{align} +# Z_c(\omega) &= \rho_0 c_0 \left[1 + 0.0571 X^{-0.754} - j 0.087 X^{-0.732}\right],\\ +# k_c(\omega) &= \frac{\omega}{c_0} \left[1 + 0.0978 X^{-0.700} - j 0.189 X^{-0.595}\right],\\ +# \end{align} +# $$ +# +# where $X = \frac{\rho_0 f}{\sigma}$. +# +# With these, we can compute the surface impedance, that in the case of a rigid passive absorber placed on a rigid wall is given by the formula +# $$ +# \begin{align} +# Z_s = -j Z_c cot(k_c d). +# \end{align} +# $$ +# +# Let's create a function to compute it. +# +# + + +# + +# Impedance calculation +def delany_bazley_layer(f, rho0, c, sigma): + X = rho0 * f / sigma + Zc = rho0 * c * (1 + 0.0571 * X**-0.754 - 1j * 0.087 * X**-0.732) + kc = 2 * np.pi * f / c * (1 + 0.0978 * (X**-0.700) - 1j * 0.189 * (X**-0.595)) + Z_s = -1j * Zc * (1 / np.tan(kc * d)) + return Z_s + + +sigma = 1.5e4 +d = 0.01 +Z_s = delany_bazley_layer(freq, rho0, c, sigma) +# - + +# Since we are going to compute a sound pressure spectrum, all the variables that depend on frequency (that are $\omega$, $k$ and $Z$) need to be updated in the frequency loop. +# To make this possible, we will initialize them as dolfinx constants. +# Then, we define the value for the normal velocity on the first end of the tube + +omega = fem.Constant(domain, default_scalar_type(0)) +k = fem.Constant(domain, default_scalar_type(0)) +Z = fem.Constant(domain, default_scalar_type(0)) +v_n = 1e-5 + +# We also need to specify the integration measure $ds$, by using ```ufl```, and its built in integration measures + +ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags) + +# ## Variational Formulation +# We can now write the variational formulation. + +# + +p = ufl.TrialFunction(V) +v = ufl.TestFunction(V) + +a = ( + ufl.inner(ufl.grad(p), ufl.grad(v)) * ufl.dx + + 1j * rho0 * omega / Z * ufl.inner(p, v) * ds(Z_bc_tag) + - k**2 * ufl.inner(p, v) * ufl.dx +) +L = -1j * omega * rho0 * ufl.inner(v_n, v) * ds(v_bc_tag) +# - + +# The class ```LinearProblem``` is used to setup the PETSc backend and assemble the system vector and matrices. +# The solution will be stored in a `dolfinx.fem.Function`, ```p_a```. + +# + +p_a = fem.Function(V) +p_a.name = "pressure" + +problem = LinearProblem( + a, + L, + u=p_a, + petsc_options={ + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps", + }, +) + + +# - + +# ## Computing the pressure at a given location +# Before starting our frequency loop, we can build a function that, given a microphone position, +# computes the sound pressure at its location. +# We will use the a similar method as in [Deflection of a membrane](../chapter1/membrane_code). +# However, as the domain doesn't deform in time, we cache the collision detection + + +class MicrophonePressure: + def __init__(self, domain, microphone_position): + """Initialize microphone(s). + + Args: + domain: The domain to insert microphones on + microphone_position: Position of the microphone(s). + Assumed to be ordered as ``(mic0_x, mic1_x, ..., mic0_y, mic1_y, ..., mic0_z, mic1_z, ...)`` + + """ + self._domain = domain + self._position = np.asarray( + microphone_position, dtype=domain.geometry.x.dtype + ).reshape(3, -1) + self._local_cells, self._local_position = self.compute_local_microphones() + + def compute_local_microphones( + self, + ) -> tuple[npt.NDArray[np.int32], npt.NDArray[np.floating]]: + """ + Compute the local microphone positions for a distributed mesh + + Returns: + Two lists (local_cells, local_points) containing the local cell indices and the local points + """ + points = self._position.T + bb_tree = geometry.bb_tree(self._domain, self._domain.topology.dim) + + cells = [] + points_on_proc = [] + + cell_candidates = geometry.compute_collisions_points(bb_tree, points) + colliding_cells = geometry.compute_colliding_cells( + domain, cell_candidates, points + ) + + for i, point in enumerate(points): + if len(colliding_cells.links(i)) > 0: + points_on_proc.append(point) + cells.append(colliding_cells.links(i)[0]) + + return np.asarray(cells, dtype=np.int32), np.asarray( + points_on_proc, dtype=domain.geometry.x.dtype + ) + + def listen( + self, recompute_collisions: bool = False + ) -> npt.NDArray[np.complexfloating]: + if recompute_collisions: + self._local_cells, self._local_position = self.compute_local_microphones() + if len(self._local_cells) > 0: + return p_a.eval(self._local_position, self._local_cells) + else: + return np.zeros(0, dtype=default_scalar_type) + + +# The pressure spectrum is initialized as a numpy array and the microphone location is assigned + +# + +p_mic = np.zeros((len(freq), 1), dtype=complex) + +mic = np.array([0.5, 0.05, 0.05]) +microphone = MicrophonePressure(domain, mic) +# - + +# ## Frequency loop +# +# Finally, we can write the frequency loop, where we update the values of the frequency-dependent variables and solve the system for each frequency + +for nf in range(0, len(freq)): + k.value = 2 * np.pi * freq[nf] / c + omega.value = 2 * np.pi * freq[nf] + Z.value = Z_s[nf] + + problem.solve() + p_a.x.scatter_forward() + + p_f = microphone.listen() + p_f = domain.comm.gather(p_f, root=0) + + if domain.comm.rank == 0: + assert p_f is not None + p_mic[nf] = np.hstack(p_f) + +# ## SPL spectrum +# After the computation, the pressure spectrum at the prescribed location is available. +# Such a spectrum is usually shown using the decibel (dB) scale to obtain the SPL, with the RMS pressure as input, +# defined as $p_{rms} = \frac{p}{\sqrt{2}}$. + +if domain.comm.rank == 0: + import matplotlib.pyplot as plt + + fig = plt.figure(figsize=(25, 8)) + plt.plot(freq, 20 * np.log10(np.abs(p_mic) / np.sqrt(2) / 2e-5), linewidth=2) + plt.grid(True) + plt.xlabel("Frequency [Hz]") + plt.ylabel("SPL [dB]") + plt.xlim([freq[0], freq[-1]]) + plt.ylim([0, 90]) + plt.legend() + plt.show()