diff --git a/.github/workflows/publish_docker.yml b/.github/workflows/publish_docker.yml index 753bf033..bd09b935 100644 --- a/.github/workflows/publish_docker.yml +++ b/.github/workflows/publish_docker.yml @@ -1,4 +1,5 @@ -name: Publish tutorial docker image +# Recipe based on: https://docs.docker.com/build/ci/github-actions/multi-platform/#distribute-build-across-multiple-runners +name: Build and publish platform dependent docker image on: push: branches: @@ -15,8 +16,12 @@ env: IMAGE_NAME: ${{ github.repository }} jobs: - build-and-push-image: - runs-on: ubuntu-latest + build: + strategy: + matrix: + os: ["ubuntu-24.04", "ubuntu-24.04-arm"] + runs-on: ${{ matrix.os }} + permissions: contents: read packages: write @@ -25,12 +30,6 @@ jobs: - name: Checkout repository uses: actions/checkout@v4 - - name: Set up QEMU - uses: docker/setup-qemu-action@v3 - - - name: Set up Docker Buildx - uses: docker/setup-buildx-action@v3 - - name: Log in to the Container registry uses: docker/login-action@v3 with: @@ -38,30 +37,81 @@ jobs: username: ${{ github.actor }} password: ${{ secrets.GITHUB_TOKEN }} + - name: Set up Docker Buildx + uses: docker/setup-buildx-action@v3 + - name: Extract metadata (tags, labels) for Docker id: meta uses: docker/metadata-action@v5 with: images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} - - name: Build Docker image + - name: Set architecture tag (amd64) + if: ${{ matrix.os == 'ubuntu-24.04' }} + run: echo "ARCH_TAG=amd64" >> $GITHUB_ENV + + - name: Set architecture tag (arm) + if: ${{ contains(matrix.os, 'arm') }} + run: echo "ARCH_TAG=arm64" >> $GITHUB_ENV + + - name: Build and push by digest + id: build uses: docker/build-push-action@v6 with: - context: . - load: true - push: false file: docker/Dockerfile - platforms: linux/amd64 - tags: ${{ steps.meta.outputs.tags }} + platforms: ${{ env.ARCH_TAG }} labels: ${{ steps.meta.outputs.labels }} + outputs: type=image,"name=${{ env.REGISTRY }}/${{ env.IMAGE_NAME}}",push-by-digest=true,name-canonical=true,push=true - - name: Build (arm) and push (amd/arm) Docker image - uses: docker/build-push-action@v6 + - name: Export digest + run: | + mkdir -p ${{ runner.temp }}/digests + digest="${{ steps.build.outputs.digest }}" + touch "${{ runner.temp }}/digests/${digest#sha256:}" + + - name: Upload digest if: github.event_name == 'push' + uses: actions/upload-artifact@v4 with: - context: . - push: true - file: docker/Dockerfile - platforms: linux/amd64,linux/arm64 - tags: ${{ steps.meta.outputs.tags }} - labels: ${{ steps.meta.outputs.labels }} + name: digests-${{ env.ARCH_TAG }} + path: ${{ runner.temp }}/digests/* + if-no-files-found: error + retention-days: 1 + + merge-and-publish: + if: github.event_name == 'push' + runs-on: ubuntu-latest + needs: + - build + steps: + - name: Download digests + uses: actions/download-artifact@v4 + with: + path: ${{ runner.temp }}/digests + pattern: digests-* + merge-multiple: true + + - name: Log in to the Container registry + uses: docker/login-action@v3 + with: + registry: ${{ env.REGISTRY }} + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + + - name: Extract metadata (tags, labels) for Docker + id: meta + uses: docker/metadata-action@v5 + with: + images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} + + + - name: Create manifest list and push + working-directory: ${{ runner.temp }}/digests + run: | + docker buildx imagetools create $(jq -cr '.tags | map("-t " + .) | join(" ")' <<< "$DOCKER_METADATA_OUTPUT_JSON") \ + $(printf '${{ env.REGISTRY }}/${{ env.IMAGE_NAME}}@sha256:%s ' *) + + - name: Inspect image + run: | + docker buildx imagetools inspect ${{ env.REGISTRY }}/${{ env.IMAGE_NAME}}:${{ steps.meta.outputs.version }} + 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/elasticity_scaling.md b/chapter2/elasticity_scaling.md index da15c47d..e6a0b2f2 100644 --- a/chapter2/elasticity_scaling.md +++ b/chapter2/elasticity_scaling.md @@ -18,6 +18,6 @@ where $\beta = 1+\frac{\lambda}{\mu}$ is a dimensionless elasticity parameter an ``` is a dimensionless variable reflecting the ratio of the load $\rho g$ and the shear stress term $\mu \nabla^2u \sim \mu \frac{U}{L^2}$ in the PDE. -One option for the scaling is to chose $U$ such that $\gamma$ is of unit size ($U=\frac{\rho g L^2}{\mu}$). However, in elasticity, this leads to displacements of the size of the geometry. This can be achieved by choosing $U$ equal to the maximum deflection of a clamped beam, for which there actually exists a formula: $U=\frac{3}{2} \rho g L^2\frac{\delta^2}{E}$ where $\delta=\frac{L}{W}$ is a parameter reflecting how slender the beam is, and $E$ is the modulus of elasticity. Thus the dimensionless parameter $\delta$ is very important in the problem (as expected $\delta\gg 1$ is what gives beam theory!). Taking $E$ to be of the same order as $\mu$, which in this case and for many materials, we realize that $\gamma \sim \delta^{-2}$ is an appropriate choice. Experimenting with the code to find a displacement that "looks right" in the plots of the deformed geometry, points to $\gamma=0.4\delta^{-2}$ as our final choice of $\gamma$. +One option for the scaling is to choose $U$ such that $\gamma$ is of unit size ($U=\frac{\rho g L^2}{\mu}$). However, in elasticity, this leads to displacements of the size of the geometry. This can be achieved by choosing $U$ equal to the maximum deflection of a clamped beam, for which there actually exists a formula: $U=\frac{3}{2} \rho g L^2\frac{\delta^2}{E}$ where $\delta=\frac{L}{W}$ is a parameter reflecting how slender the beam is, and $E$ is the modulus of elasticity. Thus the dimensionless parameter $\delta$ is very important in the problem (as expected $\delta\gg 1$ is what gives beam theory!). Taking $E$ to be of the same order as $\mu$, which in this case and for many materials, we realize that $\gamma \sim \delta^{-2}$ is an appropriate choice. Experimenting with the code to find a displacement that "looks right" in the plots of the deformed geometry, points to $\gamma=0.4\delta^{-2}$ as our final choice of $\gamma$. -The simulation code implements the problem with dimensions and physical parameters $\lambda, \mu, \rho, g, L$ and $W$. However, we can easily reuse this code for a scaled problem: Just set $\mu=\rho=L=1$, $W$ as $W/L(\delta^{-1})$, $g=\gamma$ and $\lambda=\beta$. \ No newline at end of file +The simulation code implements the problem with dimensions and physical parameters $\lambda, \mu, \rho, g, L$ and $W$. However, we can easily reuse this code for a scaled problem: Just set $\mu=\rho=L=1$, $W$ as $W/L(\delta^{-1})$, $g=\gamma$ and $\lambda=\beta$. 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..c33bcace --- /dev/null +++ b/chapter2/helmholtz_code.ipynb @@ -0,0 +1,529 @@ +{ + "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 (\n", + " fem,\n", + " io,\n", + " default_scalar_type,\n", + " geometry,\n", + " __version__ as dolfinx_version,\n", + ")\n", + "from dolfinx.fem.petsc import LinearProblem\n", + "import ufl\n", + "import numpy as np\n", + "import numpy.typing as npt\n", + "from packaging.version import Version\n", + "\n", + "mesh_data = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)\n", + "if Version(dolfinx_version) > Version(\"0.9.0\"):\n", + " domain = mesh_data.mesh\n", + " assert mesh_data.facet_tags is not None\n", + " facet_tags = mesh_data.facet_tags\n", + "else:\n", + " domain, _, facet_tags = mesh_data\n" + ] + }, + { + "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..2eee4715 --- /dev/null +++ b/chapter2/helmholtz_code.py @@ -0,0 +1,301 @@ +# --- +# 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, + __version__ as dolfinx_version, +) +from dolfinx.fem.petsc import LinearProblem +import ufl +import numpy as np +import numpy.typing as npt +from packaging.version import Version + +mesh_data = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3) +if Version(dolfinx_version) > Version("0.9.0"): + domain = mesh_data.mesh + assert mesh_data.facet_tags is not None + facet_tags = mesh_data.facet_tags +else: + domain, _, facet_tags = mesh_data + +# - + +# 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() diff --git a/chapter2/navierstokes.md b/chapter2/navierstokes.md index dd23acc3..336562cc 100644 --- a/chapter2/navierstokes.md +++ b/chapter2/navierstokes.md @@ -59,7 +59,7 @@ $\langle -\nabla \cdot \sigma, v\rangle$. Just as for the [linear elasticity pro where $T=\sigma \cdot n$ is the boundary traction. If we solve a problem with a free boundary, we can take $T=0$ on the boundary. However, if we compute the flow through a channel or a pipe and want to model flow that continues into an "imaginary channel" at the outflow, we need to treat this term with some care. The assumption we then can make is that the derivative of the velocity in the direction of the channel is zero at the outflow, corresponding to that the flow is "fully developed" or doesn't change significantly downstream at the outflow. Doing so, the remaining boundary term at the outflow becomes -$pn - \mu \nabla u \cdot n$, which is the term appearing in the variational problem [](ipcs-one). Note that this argument and the implementation depends exact on the definition of $\nabla u$, as either the matrix with components $\frac{\partial u_i}{\partial x_j}$ or $\frac{\partial u_j}{\partial x_i}$. +$pn - \mu \nabla u \cdot n$, which is the term appearing in the variational problem [](ipcs-one). Note that this argument and the implementation depend exactly on the definition of $\nabla u$, as either the matrix with components $\frac{\partial u_i}{\partial x_j}$ or $\frac{\partial u_j}{\partial x_i}$. We here choose the latter, $\frac{\partial u_j}{\partial x_i}$, which means that we must use the UFL-operator `nabla_grad`. If we use the operator `grad` and the definition $\frac{\partial u_i}{\partial x_j}$, we must instead keep the terms $pn-\mu(\nabla u)^T \cdot n$.