{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Start-to-Finish Example: Numerical Solution of the Scalar Wave Equation, in Cartesian Coordinates\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module solves the scalar wave equation in Cartesian coordinates, using the [Method of Lines](Tutorial-Method_of_Lines-C_Code_Generation.ipynb), and validates the solution against the exact one, confirming its expected convergence behavior.\n",
    "\n",
    "**Notebook Status:** <font color='green'><b>Validated</b></font>\n",
    "\n",
    "**Validation Notes:** This module has been validated to converge at the expected order to the exact solution (see [plot](#convergence) at bottom).\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [ScalarWave/ScalarWave_RHSs.py](../edit/ScalarWave/ScalarWave_RHSs.py) [\\[**tutorial**\\]](Tutorial-ScalarWave.ipynb) generates the right-hand side for the Scalar Wave Equation in Cartesian coordinates.\n",
    "* [ScalarWave/InitialData.py](../edit/ScalarWave/InitialData.py) [\\[**tutorial**\\]](Tutorial-ScalarWave.ipynb) generates C code for plane wave or spherical Gaussian initial data for the scalar wave equation.\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "As outlined in the [previous NRPy+ tutorial notebook](Tutorial-ScalarWave.ipynb), we first use NRPy+ to generate initial data for the scalar wave equation, and then we use it to generate the RHS expressions for [Method of Lines](https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html) time integration based on the [explicit Runge-Kutta fourth-order scheme](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (RK4).\n",
    "\n",
    "The entire algorithm is outlined as follows, with links to the relevant NRPy+ tutorial notebooks listed at each step.\n",
    "\n",
    "1. Allocate memory for gridfunctions, including temporary storage for the Method of Lines time integration.\n",
    "    * [**NRPy+ tutorial notebook on Method of Lines algorithm**](Tutorial-Method_of_Lines-C_Code_Generation.ipynb)\n",
    "1. Set gridfunction values to initial data.\n",
    "    * [**NRPy+ tutorial notebook section on plane-wave solution to scalar wave equation**](Tutorial-ScalarWave.ipynb#planewave)\n",
    "1. Next, integrate the initial data forward in time using the Method of Lines coupled to a Runge-Kutta explicit timestepping algorithm.\n",
    "    1. At the start of each iteration in time, output the difference between the numerical and exact solution.\n",
    "        * [**NRPy+ tutorial notebook section on plane-wave solution to scalar wave equation**](Tutorial-ScalarWave.ipynb#planewave)\n",
    "    1. At each RK time substep, do the following.\n",
    "        1. Evaluate scalar wave RHS expressions.\n",
    "            * [**NRPy+ tutorial notebook section on right-hand sides of scalar wave equation, in 3 spatial dimensions**](Tutorial-ScalarWave.ipynb#rhss3d)\n",
    "        1. Apply boundary conditions [*a la* the SENR/NRPy+ paper](https://arxiv.org/abs/1712.07658).\n",
    "1. Repeat above steps at two numerical resolutions to confirm convergence to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows.\n",
    "\n",
    "1. [Step 1](#setup): Set up core functions and parameters for solving scalar wave equation\n",
    "    1. [Step 1.a](#applybcs) `apply_bcs()`: outer boundary condition driver function\n",
    "    1. [Step 1.b](#mol) Generate Method of Lines timestepping code\n",
    "    1. [Step 1.c](#freeparams) Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h`\n",
    "1. [Step 2](#mainc): `ScalarWave_Playground.c`: The Main C Code\n",
    "1. [Step 3](#convergence): Code validation: Verify that relative error in numerical solution converges to zero at the expected order\n",
    "1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='setup'></a>\n",
    "\n",
    "# Step 1: Set up core functions and parameters for solving scalar wave equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{setup}$$\n",
    "\n",
    "Let's pick up where we left off in the [previous module](Tutorial-ScalarWave.ipynb), interfacing with the [ScalarWave/InitialData](../edit/ScalarWave/InitialData.py) and [ScalarWave/ScalarWave_RHSs](../edit/ScalarWave/ScalarWave_RHSs.py) NRPy+ modules to generate\n",
    "* monochromatic (single-wavelength) plane wave scalar wave initial data, and\n",
    "* the scalar wave equation RHSs at **4th** finite difference order in **3 spatial dimensions**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.410919Z",
     "iopub.status.busy": "2021-10-04T19:30:07.410301Z",
     "iopub.status.idle": "2021-10-04T19:30:07.642502Z",
     "shell.execute_reply": "2021-10-04T19:30:07.642912Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step P1: Import needed NRPy+ core modules:\n",
    "from outputC import lhrh, add_to_Cfunction_dict  # NRPy+: Core C code output module\n",
    "import finite_difference as fin       # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par        # NRPy+: Parameter interface\n",
    "import grid as gri                    # NRPy+: Functions having to do with numerical grids\n",
    "import cmdline_helper as cmd          # NRPy+: Multi-platform Python command-line interface\n",
    "import shutil, os                     # Standard Python modules for multiplatform OS-level functions\n",
    "\n",
    "# Step P2: Create C code output directory:\n",
    "Ccodesrootdir = os.path.join(\"ScalarWave_Ccodes\")\n",
    "# First remove C code output directory if it exists\n",
    "shutil.rmtree(Ccodesrootdir, ignore_errors=True)\n",
    "# Then create a fresh directory\n",
    "cmd.mkdir(Ccodesrootdir)\n",
    "\n",
    "# Step P3: Create executable output directory:\n",
    "outdir = os.path.join(Ccodesrootdir, \"output\")\n",
    "cmd.mkdir(outdir)\n",
    "\n",
    "# Step P4: Set domain_size, the physical extent of numerical grid;\n",
    "#          in Cartesian coordinates xmin=ymin=zmin=-domain_size,\n",
    "#          and xmax=ymax=zmax=+domain_size\n",
    "domain_size = 10.0\n",
    "\n",
    "# Step P5: Set timestepping algorithm (we adopt the Method of Lines)\n",
    "RK_method = \"RK4\"\n",
    "\n",
    "# Step P6: Set the finite differencing order to 4.\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 4)\n",
    "\n",
    "# Step P7: Enable/disable SIMD. If enabled, code should run ~2x faster on most CPUs.\n",
    "enable_SIMD = True\n",
    "\n",
    "# Step 1: Import the ScalarWave.InitialData module.\n",
    "#         This command only declares ScalarWave initial data\n",
    "#         parameters and the InitialData() function.\n",
    "import ScalarWave.InitialData as swid\n",
    "\n",
    "# Step 2: Import ScalarWave_RHSs module.\n",
    "#         This command only declares ScalarWave RHS parameters\n",
    "#         and the ScalarWave_RHSs function (called later)\n",
    "import ScalarWave.ScalarWave_RHSs as swrhs\n",
    "\n",
    "# Step 3: The spatial dimension parameter defaults to 3;\n",
    "#         no need to set!\n",
    "# par.set_parval_from_str(\"grid::DIM\", 3)\n",
    "\n",
    "# Step 4: Call the InitialData() function to set up initial data.\n",
    "#         Options include:\n",
    "#    \"PlaneWave\": monochromatic (single frequency/wavelength) plane wave\n",
    "#    \"SphericalGaussian\": spherically symmetric Gaussian, with default stdev=3\n",
    "swid.InitialData(WaveType=\"PlaneWave\")\n",
    "\n",
    "# Step 5: Generate SymPy symbolic expressions for\n",
    "#         uu_rhs and vv_rhs; the ScalarWave RHSs.\n",
    "#         This function also declares the uu and vv\n",
    "#         gridfunctions, which need to be declared\n",
    "#         to output even the initial data to C file.\n",
    "swrhs.ScalarWave_RHSs()\n",
    "\n",
    "# Step 6: Enable \"FD functions\". In other words, all finite-difference stencils\n",
    "#         will be output as inlined static functions. This is essential for\n",
    "#         compiling highly complex FD kernels with using certain versions of GCC;\n",
    "#         GCC 10-ish will choke on BSSN FD kernels at high FD order, sometimes\n",
    "#         taking *hours* to compile. Unaffected GCC versions compile these kernels\n",
    "#         in seconds. FD functions do not slow the code performance, but do add\n",
    "#         another header file to the C source tree.\n",
    "enable_FD_functions = True\n",
    "par.set_parval_from_str(\"finite_difference::enable_FD_functions\", enable_FD_functions)\n",
    "\n",
    "# Step 7: If enable_SIMD, then copy SIMD/SIMD_intrinsics.h to $Ccodesrootdir/SIMD/SIMD_intrinsics.h\n",
    "cmd.mkdir(os.path.join(Ccodesrootdir,\"SIMD\"))\n",
    "if enable_SIMD:\n",
    "    shutil.copy(os.path.join(\"SIMD\", \"SIMD_intrinsics.h\"), os.path.join(Ccodesrootdir, \"SIMD\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.647237Z",
     "iopub.status.busy": "2021-10-04T19:30:07.646715Z",
     "iopub.status.idle": "2021-10-04T19:30:07.648540Z",
     "shell.execute_reply": "2021-10-04T19:30:07.648204Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_exact_solution_single_point():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc = \"Exact solution at a single point. params.time==0 corresponds to the initial data.\"\n",
    "    c_type = \"void\"\n",
    "    name = \"exact_solution_single_point\"\n",
    "    params = \"\"\"const paramstruct *restrict params,\n",
    "                const REAL xx0, const REAL xx1, const REAL xx2,\n",
    "                REAL *uu_exact, REAL *vv_exact\"\"\"\n",
    "    body = fin.FD_outputC(\"returnstring\",[lhrh(lhs=\"*uu_exact\",rhs=swid.uu_ID),\n",
    "                                          lhrh(lhs=\"*vv_exact\",rhs=swid.vv_ID)],\n",
    "                          params=\"includebraces=False,preindent=1,outCverbose=False\")\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.652260Z",
     "iopub.status.busy": "2021-10-04T19:30:07.651854Z",
     "iopub.status.idle": "2021-10-04T19:30:07.653570Z",
     "shell.execute_reply": "2021-10-04T19:30:07.653889Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_exact_solution_all_points():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc = \"Exact solution at all points. params.time==0 corresponds to the initial data.\"\n",
    "    c_type = \"void\"\n",
    "    name = \"exact_solution_all_points\"\n",
    "    params = \"const paramstruct *restrict params,REAL *restrict xx[3], REAL *restrict in_gfs\"\n",
    "    body = \"\"\"exact_solution_single_point(params, xx0, xx1, xx2,\n",
    "           &in_gfs[IDX4S(UUGF,i0,i1,i2)], &in_gfs[IDX4S(VVGF,i0,i1,i2)]);\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"), loopopts = \"AllPoints,Read_xxs\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.659176Z",
     "iopub.status.busy": "2021-10-04T19:30:07.658769Z",
     "iopub.status.idle": "2021-10-04T19:30:07.660491Z",
     "shell.execute_reply": "2021-10-04T19:30:07.660807Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_rhs_eval():\n",
    "    desc=\"Evaluate the scalar wave RHSs\"\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    if enable_FD_functions:\n",
    "        includes += [\"finite_difference_functions.h\"]\n",
    "    if enable_SIMD:\n",
    "        includes += [\"SIMD/SIMD_intrinsics.h\"]\n",
    "    c_type = \"void\"\n",
    "    name = \"rhs_eval\"\n",
    "    params = \"const paramstruct *restrict params, const REAL *restrict in_gfs, REAL *restrict rhs_gfs\"\n",
    "    body = fin.FD_outputC(\"returnstring\",[lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"uu\"),rhs=swrhs.uu_rhs),\n",
    "                                          lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"vv\"),rhs=swrhs.vv_rhs)],\n",
    "                          params=\"enable_SIMD=\"+str(enable_SIMD))\n",
    "    loopopts = \"InteriorPoints\"\n",
    "    if enable_SIMD:\n",
    "        loopopts += \",enable_SIMD\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"), loopopts = loopopts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.665143Z",
     "iopub.status.busy": "2021-10-04T19:30:07.664729Z",
     "iopub.status.idle": "2021-10-04T19:30:07.666695Z",
     "shell.execute_reply": "2021-10-04T19:30:07.666343Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_diagnostic_output_2D_xy_plane():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc=\"As a diagnostic, output to file 2D data closest to the xy plane (z approx 0 plane).\"\n",
    "    c_type = \"void\"\n",
    "    name = \"diagnostic_output_2D_xy_plane\"\n",
    "    params = \"\"\"const paramstruct *restrict params, REAL *xx[3],\n",
    "               const REAL *numerical_gridfunction_data,REAL *gridfunction_to_store_exact\"\"\"\n",
    "    body = r\"\"\"  char filename[100];\n",
    "  sprintf(filename, \"output/out2D__resolution_%dx%dx%d__iter_%d.txt\", Nxx0,Nxx1,Nxx2,params->n);\n",
    "  FILE *out2D = fopen(filename, \"w\");\n",
    "  // Output on z=midpoint plane (points closest to z=0):\n",
    "  const int i2  = (int)((Nxx2+ 2*NGHOSTS)*0.5);\n",
    "  const REAL xx2 = xx[2][i2];\n",
    "  for(int i0=0;i0<Nxx0+2*NGHOSTS;i0++) {\n",
    "    for(int i1=0;i1<Nxx1+2*NGHOSTS;i1++) {\n",
    "      // Zoom in; do not output data near outer boundary.\n",
    "      if(i0> (Nxx0+2*NGHOSTS)*.25 && i0< (Nxx0+2*NGHOSTS)*.75 &&\n",
    "         i1> (Nxx1+2*NGHOSTS)*.25 && i1< (Nxx1+2*NGHOSTS)*.75) {\n",
    "        const REAL xx0 = xx[0][i0];\n",
    "        const REAL xx1 = xx[1][i1];\n",
    "        REAL uu_exact,vv_exact; exact_solution_single_point(params,xx0,xx1,xx2, &uu_exact,&vv_exact);\n",
    "        fprintf(out2D,\"%e %e %e %e\\n\", xx0, xx1,\n",
    "                numerical_gridfunction_data[IDX4S(0,i0,i1, (int)((Nxx2+ 2*NGHOSTS)*0.5))], uu_exact);\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "  fclose(out2D);\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='applybcs'></a>\n",
    "\n",
    "## Step 1.a: `apply_bcs()`: outer boundary condition driver function \\[Back to [top](#toc)\\]\n",
    "$$\\label{applybcs}$$\n",
    "\n",
    "When solving the wave equation on a 3D Cartesian numerical grid cube (or, if you like, rectangular prism), at each step in time, we first evaluate the right-hand sides (RHSs) of the $\\partial_t u$ and $\\partial_t v$ equations. \n",
    "\n",
    "These RHSs generally contain spatial derivatives, which we evaluate using finite-difference differentiation ([**tutorial**](Tutorial-Finite_Difference_Derivatives.ipynb)). Each finite-difference derivative depends on neighboring points on the left and right, so the RHSs can only be evaluated in the grid interior. For example, a standard fourth-order centered finite difference derivative depends on two points to the left and right of the point at which the derivative is being evaluated. In order for the same interior to be filled at the next time step, we need to fill in the data at the boundaries; i.e., we need to apply boundary conditions.\n",
    "\n",
    "Here we quadratically extrapolate data to the outer boundary using the `FACE_UPDATE()` C macro defined below. The C code function `apply_bcs()` below updates all 6 faces of the cube. To ensure that all gridpoints on the outer boundary (also known as \"ghost cells\") are filled, the following algorithm is implemented, starting at the innermost ghost cells (i.e., the ghost cells closest to the grid interior):\n",
    "\n",
    "1. The lower $x$ face is updated on only the interior points of the face.\n",
    "1. The upper $x$ face is updated on only the interior points of the face.\n",
    "1. The lower $y$ face is updated on the interior points of that face, plus the lower and upper $x$ boundary points.\n",
    "1. The upper $y$ face is updated on the interior points of that face, plus the lower and upper $x$ boundary points.\n",
    "1. The lower $z$ face is updated on the interior points of that face, plus the lower and upper $x$ boundary points, plus the lower and upper $y$ boundary points.\n",
    "1. The upper $z$ face is updated on the interior points of that face, plus the lower and upper $x$ boundary points, plus the lower and upper $y$ boundary points.\n",
    "1. The above is repeated on the next outer ghost cell until all outer boundary points are filled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.671571Z",
     "iopub.status.busy": "2021-10-04T19:30:07.671158Z",
     "iopub.status.idle": "2021-10-04T19:30:07.673048Z",
     "shell.execute_reply": "2021-10-04T19:30:07.672713Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_apply_bcs():\n",
    "    desc=\"\"\"Apply (quadratic extrapolation) spatial boundary conditions to the scalar wave gridfunctions.\n",
    "BCs are applied to all six boundary faces of the cube, filling in the innermost\n",
    "ghost zone first, and moving outward.\"\"\"\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    c_type = \"void\"\n",
    "    name = \"apply_bcs\"\n",
    "    params = \"const paramstruct *restrict params,REAL *restrict gfs\"\n",
    "    prefunc = r\"\"\"\n",
    "// Declare boundary condition FACE_UPDATE macro,\n",
    "//          which updates a single face of the 3D grid cube\n",
    "//          using quadratic polynomial extrapolation.\n",
    "const int MAXFACE = -1;\n",
    "const int NUL     = +0;\n",
    "const int MINFACE = +1;\n",
    "#define  FACE_UPDATE(which_gf, i0min,i0max, i1min,i1max, i2min,i2max, FACEX0,FACEX1,FACEX2) \\\n",
    "  for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++) { \\\n",
    "        gfs[IDX4S(which_gf,i0,i1,i2)] =                                  \\\n",
    "          +3.0*gfs[IDX4S(which_gf,i0+1*FACEX0,i1+1*FACEX1,i2+1*FACEX2)]  \\\n",
    "          -3.0*gfs[IDX4S(which_gf,i0+2*FACEX0,i1+2*FACEX1,i2+2*FACEX2)]  \\\n",
    "          +1.0*gfs[IDX4S(which_gf,i0+3*FACEX0,i1+3*FACEX1,i2+3*FACEX2)]; \\\n",
    "      }\n",
    "\"\"\"\n",
    "    body = r\"\"\"\n",
    "#pragma omp parallel for\n",
    "    for(int which_gf=0;which_gf<NUM_EVOL_GFS;which_gf++) {\n",
    "#include \"set_Cparameters.h\"\n",
    "      int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };\n",
    "      int imax[3] = { Nxx_plus_2NGHOSTS0-NGHOSTS, Nxx_plus_2NGHOSTS1-NGHOSTS, Nxx_plus_2NGHOSTS2-NGHOSTS };\n",
    "      for(int which_gz = 0; which_gz < NGHOSTS; which_gz++) {\n",
    "        // After updating each face, adjust imin[] and imax[]\n",
    "        //   to reflect the newly-updated face extents.\n",
    "        FACE_UPDATE(which_gf, imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL); imin[0]--;\n",
    "        FACE_UPDATE(which_gf, imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL); imax[0]++;\n",
    "\n",
    "        FACE_UPDATE(which_gf, imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL); imin[1]--;\n",
    "        FACE_UPDATE(which_gf, imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL); imax[1]++;\n",
    "\n",
    "        FACE_UPDATE(which_gf, imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE); imin[2]--;\n",
    "        FACE_UPDATE(which_gf, imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE); imax[2]++;\n",
    "    }\n",
    "  }\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        prefunc=prefunc, body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mol'></a>\n",
    "\n",
    "## Step 1.b: Generate Method of Lines timestepping code \\[Back to [top](#toc)\\]\n",
    "$$\\label{mol}$$\n",
    "\n",
    "The Method of Lines algorithm is described in detail in the [**NRPy+ tutorial notebook on Method of Lines algorithm**](Tutorial-Method_of_Lines-C_Code_Generation.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.675995Z",
     "iopub.status.busy": "2021-10-04T19:30:07.675602Z",
     "iopub.status.idle": "2021-10-04T19:30:07.699606Z",
     "shell.execute_reply": "2021-10-04T19:30:07.699248Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1.b: Generate Runge-Kutta-based (RK-based) timestepping code.\n",
    "#       As described above the Table of Contents, this is a 2-step process:\n",
    "#       1.b.A: Evaluate RHSs (RHS_string)\n",
    "#       1.b.B: Apply boundary conditions (post_RHS_string, pt 1)\n",
    "import MoLtimestepping.MoL_simple as MoL\n",
    "# from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict\n",
    "# RK_order  = Butcher_dict[RK_method][1]\n",
    "# cmd.mkdir(os.path.join(Ccodesrootdir,\"MoLtimestepping/\"))\n",
    "MoL.register_C_functions_and_NRPy_basic_defines(RK_method,\n",
    "    RHS_string      =  \"rhs_eval(params,   RK_INPUT_GFS, RK_OUTPUT_GFS);\",\n",
    "    post_RHS_string = \"apply_bcs(params,  RK_OUTPUT_GFS);\", enable_SIMD=enable_SIMD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='freeparams'></a>\n",
    "\n",
    "## Step 1.c: Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h` \\[Back to [top](#toc)\\]\n",
    "$$\\label{freeparams}$$\n",
    "\n",
    "Here we output `free_parameters.h`, which sets initial data parameters, as well as grid domain & reference metric parameters, applying `domain_size` and `sinh_width`/`SymTP_bScale` (if applicable) as set above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.702625Z",
     "iopub.status.busy": "2021-10-04T19:30:07.702219Z",
     "iopub.status.idle": "2021-10-04T19:30:07.704039Z",
     "shell.execute_reply": "2021-10-04T19:30:07.703653Z"
    }
   },
   "outputs": [],
   "source": [
    "domain_size_str=str(domain_size)\n",
    "# Step 3.d: Set free_parameters.h\n",
    "with open(os.path.join(Ccodesrootdir,\"free_parameters.h\"),\"w\") as file:\n",
    "    file.write(r\"\"\"\n",
    "// Set free-parameter values.\n",
    "\n",
    "// Set free-parameter values for the initial data.\n",
    "params.time = 0.0; params.wavespeed = 1.0;\n",
    "//params.kk0 = 1.0; params.kk1 = 1.0; params.kk2 = 1.0;\n",
    "\n",
    "const REAL domain_size    = \"\"\"+str(domain_size)+r\"\"\";\n",
    "\n",
    "// Override parameter defaults with values based on command line arguments and NGHOSTS.\n",
    "const int Nx0x1x2 = atoi(argv[1]);\n",
    "params.Nxx0 = Nx0x1x2;\n",
    "params.Nxx1 = Nx0x1x2;\n",
    "params.Nxx2 = Nx0x1x2;\n",
    "params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*NGHOSTS;\n",
    "params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*NGHOSTS;\n",
    "params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*NGHOSTS;\n",
    "// Step 0d: Set up space and time coordinates\n",
    "// Step 0d.i: Declare \\Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:\n",
    "const REAL xxmin[3] = {-\"\"\"+domain_size_str+\"\"\",-\"\"\"+domain_size_str+\"\"\",-\"\"\"+domain_size_str+\"\"\" };\n",
    "const REAL xxmax[3] = {+\"\"\"+domain_size_str+\"\"\",+\"\"\"+domain_size_str+\"\"\",+\"\"\"+domain_size_str+r\"\"\" };\n",
    "\n",
    "params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0);\n",
    "params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1);\n",
    "params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2);\n",
    "params.invdx0 = 1.0 / params.dxx0;\n",
    "params.invdx1 = 1.0 / params.dxx1;\n",
    "params.invdx2 = 1.0 / params.dxx2;\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mainc'></a>\n",
    "\n",
    "# Step 2: `ScalarWave_Playground`: The Main C Code \\[Back to [top](#toc)\\]\n",
    "$$\\label{mainc}$$\n",
    "\n",
    "Next we will write the C code infrastructure necessary to make use of the above NRPy+-generated codes. Again, we'll be using RK4 time integration via the Method of Lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.710026Z",
     "iopub.status.busy": "2021-10-04T19:30:07.705672Z",
     "iopub.status.idle": "2021-10-04T19:30:07.711810Z",
     "shell.execute_reply": "2021-10-04T19:30:07.711480Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_main__ScalarWave_Playground():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    prefunc = \"\"\"\n",
    "// Frequency of output.\n",
    "const int NSKIP_0D_OUTPUT = 1;\n",
    "const int NSKIP_2D_OUTPUT = 10;\n",
    "\"\"\"\n",
    "    desc = \"\"\"\n",
    "// main() function:\n",
    "// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates\n",
    "// Step 1: Set up scalar wave initial data\n",
    "// Step 2: Evolve scalar wave initial data forward in time using Method of Lines with RK4 algorithm,\n",
    "//         applying quadratic extrapolation outer boundary conditions.\n",
    "// Step 3: Output relative error between numerical and exact solution.\n",
    "// Step 4: Free all allocated memory\n",
    "\"\"\"\n",
    "    c_type = \"int\"\n",
    "    name = \"main\"\n",
    "    params = \"int argc, const char *argv[]\"\n",
    "    body = r\"\"\"\n",
    "  paramstruct params;\n",
    "  set_Cparameters_to_default(&params);\n",
    "\n",
    "  // Step 0a: Read command-line input, error out if nonconformant\n",
    "  if(argc != 2 || atoi(argv[1]) < NGHOSTS) {\n",
    "    printf(\"Error: Expected one command-line argument: ./ScalarWave_Playground [Nx(=Ny=Nz)],\\n\");\n",
    "    printf(\"where Nx is the number of grid points in the x,y, and z directions.\\n\");\n",
    "    printf(\"Nx MUST BE larger than NGHOSTS (= %d)\\n\",NGHOSTS);\n",
    "    exit(1);\n",
    "  }\n",
    "  if(atoi(argv[1])%2 != 0) {\n",
    "    printf(\"Error: Algorithm for setting up cell-centered grids here requires Nx, Ny, and Nz to be a multiple of 2 .\\n\");\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  // Step 0b: Set free parameters, overwriting Cparameters defaults\n",
    "  //          by hand or with command-line input, as desired.\n",
    "#include \"free_parameters.h\"\n",
    "  //          ... and then set up the numerical grid structure in time:\n",
    "  const REAL CFL_FACTOR = 0.5; // Set the CFL Factor\n",
    "  #define MIN(A, B) ( ((A) < (B)) ? (A) : (B) )\n",
    "  params.dt = CFL_FACTOR * MIN(params.dxx0,MIN(params.dxx1,params.dxx2)); // CFL condition\n",
    "\n",
    "  // Step 0c: Now that params struct has been properly set up, create\n",
    "  //   list of const's containing each parameter. E.g.,\n",
    "  //  const REAL dxx0 = params.dxx0;\n",
    "#include \"set_Cparameters-nopointer.h\"\n",
    "\n",
    "  // Step 0d: Declare struct for gridfunctions and allocate memory for gridfunctions\n",
    "  MoL_gridfunctions_struct gridfuncs;\n",
    "  MoL_malloc_y_n_gfs(&params, &gridfuncs);\n",
    "  MoL_malloc_non_y_n_gfs(&params, &gridfuncs);\n",
    "\n",
    "  // Step 0e: Set t_final, and number of timesteps based on t_final\n",
    "  params.t_final = xxmax[0]*0.8; /* Final time is set so that at t=t_final,\n",
    "                                    data at the origin have not been corrupted\n",
    "                                    by the approximate outer boundary condition */\n",
    "  int Nt = (int)(params.t_final / params.dt + 0.5); // The number of points in time.\n",
    "                                                    //Add 0.5 to account for C rounding down integers.\n",
    "\n",
    "  // Step 0f: Set up cell-centered Cartesian coordinate grids\n",
    "  REAL *xx[3];\n",
    "  xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);\n",
    "  xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);\n",
    "  xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);\n",
    "  for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] = xxmin[0] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*params.dxx0;\n",
    "  for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*params.dxx1;\n",
    "  for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*params.dxx2;\n",
    "\n",
    "  // Step 1: Set up initial data to be exact solution at time=0:\n",
    "  params.time = 0.0; exact_solution_all_points(&params, xx, gridfuncs.y_n_gfs);\n",
    "\n",
    "  while(params.time < params.t_final) { // Main loop to progress forward in time.\n",
    "    // Step 1a: Set current time to correct value & compute exact solution\n",
    "    params.time = ((REAL)params.n)*params.dt;\n",
    "\n",
    "    /* Step 2: Validation: Output relative error between numerical and exact solution, */\n",
    "    // Step 2b: Output to 2D grid (region of x-y plane near origin)\n",
    "    //          every NSKIP_2D_OUTPUT iterations.\n",
    "    if((params.n)%NSKIP_2D_OUTPUT ==0) {\n",
    "      // Step 2a: Output data on gridpoints closest to xy-plane, ignoring points \"too close\" to outer boundary\n",
    "      diagnostic_output_2D_xy_plane(&params, xx, gridfuncs.y_n_gfs, gridfuncs.diagnostic_output_gfs);\n",
    "    }\n",
    "\n",
    "    if((n)%NSKIP_0D_OUTPUT ==0) {\n",
    "      // Step 2c: Output relative error between exact & numerical at center of grid.\n",
    "      const int i0mid=Nxx_plus_2NGHOSTS0/2;\n",
    "      const int i1mid=Nxx_plus_2NGHOSTS1/2;\n",
    "      const int i2mid=Nxx_plus_2NGHOSTS2/2;\n",
    "      REAL uu_exact,vv_exact; exact_solution_single_point(&params, xx[0][i0mid],xx[1][i1mid],xx[2][i2mid],\n",
    "                                                          &uu_exact,&vv_exact);\n",
    "\n",
    "      const REAL numerical = gridfuncs.y_n_gfs[IDX4S(UUGF,i0mid,i1mid,i2mid)];\n",
    "      const REAL relative_error = fabs((uu_exact-numerical)/uu_exact) + 1e-16; // Add 1e-16 in case relerror\n",
    "                                                                               // is exactly zero.\n",
    "      printf(\"%e %e || %e %e %e: %e %e\\n\",params.time, log10(relative_error),\n",
    "             xx[0][i0mid],xx[1][i1mid],xx[2][i2mid], numerical,uu_exact);\n",
    "    }\n",
    "\n",
    "    // Step 3: Evolve scalar wave initial data forward in time using Method of Lines with RK4 algorithm,\n",
    "    //         applying quadratic extrapolation outer boundary conditions.\n",
    "    // Step 3.b: Step forward one timestep (t -> t+dt) in time using\n",
    "    //           chosen RK-like MoL timestepping algorithm\n",
    "    MoL_step_forward_in_time(&params, xx, &gridfuncs);\n",
    "  } // End main loop to progress forward in time.\n",
    "\n",
    "  // Step 4: Free all allocated memory\n",
    "  MoL_free_memory_y_n_gfs(&params, &gridfuncs);\n",
    "  MoL_free_memory_non_y_n_gfs(&params, &gridfuncs);\n",
    "  for(int i=0;i<3;i++) free(xx[i]);\n",
    "  return 0;\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes, prefunc=prefunc,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"), enableCparameters=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.716463Z",
     "iopub.status.busy": "2021-10-04T19:30:07.716077Z",
     "iopub.status.idle": "2021-10-04T19:30:07.795372Z",
     "shell.execute_reply": "2021-10-04T19:30:07.794991Z"
    }
   },
   "outputs": [],
   "source": [
    "def register_C_code_functions_scalarwave():\n",
    "    add_to_Cfunction_dict_exact_solution_single_point()\n",
    "    add_to_Cfunction_dict_exact_solution_all_points()\n",
    "    add_to_Cfunction_dict_rhs_eval()\n",
    "    add_to_Cfunction_dict_diagnostic_output_2D_xy_plane()\n",
    "    add_to_Cfunction_dict_apply_bcs()\n",
    "    add_to_Cfunction_dict_main__ScalarWave_Playground()\n",
    "\n",
    "import outputC as outC\n",
    "outC.outputC_register_C_functions_and_NRPy_basic_defines()  # #define M_PI,  etc.\n",
    "# Declare paramstruct, register set_Cparameters_to_default(),\n",
    "#   and output declare_Cparameters_struct.h and set_Cparameters[].h:\n",
    "outC.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(os.path.join(Ccodesrootdir))\n",
    "\n",
    "gri.register_C_functions_and_NRPy_basic_defines(enable_griddata_struct=False)  # #define IDX3S(),  etc.\n",
    "fin.register_C_functions_and_NRPy_basic_defines(NGHOSTS_account_for_onezone_upwind=False)  # #define NGHOSTS, etc.\n",
    "\n",
    "# all functions needed for scalar wave:\n",
    "register_C_code_functions_scalarwave()\n",
    "\n",
    "# Output functions for computing all finite-difference stencils.\n",
    "#   Must be called after defining all functions depending on FD stencils.\n",
    "if enable_FD_functions:\n",
    "    fin.output_finite_difference_functions_h(path=Ccodesrootdir)\n",
    "\n",
    "# Call this last: Set up NRPy_basic_defines.h and NRPy_function_prototypes.h.\n",
    "outC.construct_NRPy_basic_defines_h(Ccodesrootdir, enable_SIMD=enable_SIMD)\n",
    "outC.construct_NRPy_function_prototypes_h(Ccodesrootdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:07.798471Z",
     "iopub.status.busy": "2021-10-04T19:30:07.798087Z",
     "iopub.status.idle": "2021-10-04T19:30:09.478368Z",
     "shell.execute_reply": "2021-10-04T19:30:09.477865Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(EXEC): Executing `make -j18`...\n",
      "(BENCH): Finished executing in 0.80 seconds.\n",
      "Finished compilation.\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./ScalarWave_Playground 48`...\n",
      "(BENCH): Finished executing in 0.20 seconds.\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./ScalarWave_Playground 64`...\n",
      "(BENCH): Finished executing in 0.60 seconds.\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./ScalarWave_Playground 96`...\n",
      "(BENCH): Finished executing in 1.61 seconds.\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd\n",
    "cmd.new_C_compile(Ccodesrootdir, \"ScalarWave_Playground\",\n",
    "                  uses_free_parameters_h=True, compiler_opt_option=\"fast\") # fastdebug or debug also supported\n",
    "os.chdir(Ccodesrootdir)\n",
    "for res in [\"48\", \"64\", \"96\"]:\n",
    "    cmd.Execute(\"ScalarWave_Playground\", res, os.path.join(\"output\", \"out\"+res+\".txt\"))\n",
    "os.chdir(\"..\")\n",
    "# import sys\n",
    "# sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convergence'></a>\n",
    "\n",
    "# Step 3: Code Validation: Verify that relative error in numerical solution converges to zero at the expected order \\[Back to [top](#toc)\\]\n",
    "$$\\label{convergence}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:09.490232Z",
     "iopub.status.busy": "2021-10-04T19:30:09.489648Z",
     "iopub.status.idle": "2021-10-04T19:30:09.817378Z",
     "shell.execute_reply": "2021-10-04T19:30:09.817159Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHHCAYAAABTMjf2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAACZDElEQVR4nOzddXhT59sH8O+Jp97SFmuRFnd3t+Ky4WNQKPryGwzbYBu2DRnOGAwdMsZwGwwdxYZDcQpFajgUammjz/tHdkJDkjZtk6Zyf64r18bJycl9ktPkzv0YxxhjIIQQQggpgASODoAQQgghxFEoESKEEEJIgUWJECGEEEIKLEqECCGEEFJgUSJECCGEkAKLEiFCCCGEFFiUCBFCCCGkwKJEiBBCCCEFFiVChBBCCCmwKBEidnXy5ElwHIeTJ086OhSSgziOw4wZMxwdRpbw1+zOnTsdHYpZpUqVQnBwsKPDICTfoESIZMmGDRvAcZzhJpPJUK5cOfzvf//Dy5cvbfIcf//9d6a+TFu0aGGIRyAQwM3NDeXLl8fnn3+OY8eO2SSmvEShUGDGjBl2S0Iz+/44wtmzZw3XxJs3b4zu27JlC5YsWeKYwPIQrVaL9evXo0WLFvDy8oJUKkWpUqUwePBgXLlyxdHhEZJtIkcHQPK277//HqVLl0ZqairOnj2LX3/9FX///Tdu374NJyenbB3777//xvLlyzP1Zevn54c5c+YAAJKTk/Hw4UPs3r0bmzdvRu/evbF582aIxeJsxZVXKBQKzJw5E4A+SbS19N6flJQUiESO/XjR6XT44osv4OzsjOTkZJP7t2zZgtu3b+PLL7/M+eDyiJSUFHzyySc4fPgwmjVrhm+++QZeXl6IjIzE9u3bsXHjRkRHR8PPz8/RoRKSZZQIkWzp0KED6tSpAwAYOnQoChUqhEWLFmHfvn3o169fjsfj7u6OAQMGGG2bO3cuxowZgxUrVqBUqVL46aefcjyuvCA5ORnOzs42OZZMJrPJcbJj9erViImJwdChQ7F06VJHh2M1W74PGdFoNNDpdJBIJGbvnzRpEg4fPozFixebJIzTp0/H4sWLcyBK+1IoFNn+0UbyNmoaIzbVqlUrAMCTJ0/S3W/Hjh2oXbs25HI5vL29MWDAADx9+tRwf3BwMJYvXw4ARk1wWSEUCvHzzz+jUqVK+OWXXxAfH290/+bNmw2xeHl5oW/fvoiJiTHap0WLFqhSpQpu3ryJ5s2bw8nJCWXKlDH0Izl16hTq168PuVyO8uXL4/jx4yZxhIWFoUOHDnBzc4OLiwtat26NCxcuGO3DNzn++++/GD9+PHx8fODs7IwePXrg9evXRvteuXIFQUFB8Pb2hlwuR+nSpTFkyBAAQGRkJHx8fAAAM2fONLx+fPUmODgYLi4uePToETp27AhXV1d89tlnAIAzZ86gV69eKFGiBKRSKfz9/TFu3DikpKQYnjuj9+fjPkIzZswAx3F4+PAhgoOD4eHhAXd3dwwePBgKhcLovFJSUjBmzBh4e3vD1dUVXbt2xdOnTzPV7yguLg7fffcdvv/+e3h4eJjc36JFCxw8eBBRUVGG2EuVKmW0j06nw6xZs+Dn5weZTIbWrVvj4cOHVj0/AJw4cQJNmzaFs7MzPDw80K1bN9y7d89oH/51uXv3Lvr37w9PT080adIEAMAYw48//gg/Pz84OTmhZcuWuHPnjtnnev/+Pb788kv4+/tDKpWiTJky+Omnn6DT6Qz7REZGguM4LFiwAEuWLEFgYCCkUinu3r1r9pixsbFYtWoV2rZta7ZqJhQKMXHiRKNqkC2v8c6dOyMgIMBsbA0bNjT8AONl5u/46tWraNasGZycnPDNN98AAN6+fYvPP/8cbm5u8PDwwKBBg3Djxg1wHIcNGzYYHSc8PBw9e/aEl5cXZDIZ6tSpg/3792fpPHmHDh1C8+bN4erqCjc3N9StWxdbtmwx2ufixYto37493N3d4eTkhObNm+Pff/81+xoR61FFiNjUo0ePAACFChWyuM+GDRswePBg1K1bF3PmzMHLly+xdOlS/PvvvwgLC4OHhwdGjBiBZ8+e4dixY/j999+zHZdQKES/fv0wdepUnD17Fp06dQIAzJo1C1OnTkXv3r0xdOhQvH79GsuWLUOzZs0MsfDevXuHzp07o2/fvujVqxd+/fVX9O3bF3/88Qe+/PJLjBw5Ev3798f8+fPRs2dPxMTEwNXVFQBw584dNG3aFG5ubvjqq68gFouxatUqtGjRwpBEpfXFF1/A09MT06dPR2RkJJYsWYL//e9/2LZtGwDg1atXaNeuHXx8fDB58mR4eHggMjISu3fvBgD4+Pjg119/xahRo9CjRw988sknAIBq1aoZnkOj0SAoKAhNmjTBggULDL+Kd+zYAYVCgVGjRqFQoUK4dOkSli1bhtjYWOzYsQMAsvz+9O7dG6VLl8acOXNw7do1rF27Fr6+vkZVuuDgYGzfvh2ff/45GjRogFOnThneL2tNnToVRYoUwYgRI/DDDz+Y3P/tt98iPj4esbGxhqqGi4uL0T5z586FQCDAxIkTER8fj3nz5uGzzz7DxYsXM3z+48ePo0OHDggICMCMGTOQkpKCZcuWoXHjxrh27ZpJ0tWrVy+ULVsWs2fPBmMMADBt2jT8+OOP6NixIzp27Ihr166hXbt2UKlURo9VKBRo3rw5nj59ihEjRqBEiRI4d+4cpkyZgufPn5v0g1q/fj1SU1MxfPhwSKVSeHl5mT2HQ4cOQaPR4PPPP8/wfAHbX+N9+vTBwIEDcfnyZdStW9fwuKioKFy4cAHz5883bMvM3/Hbt2/RoUMH9O3bFwMGDEDhwoWh0+nQpUsXXLp0CaNGjUKFChWwb98+DBo0yOx5Nm7cGMWLF8fkyZPh7OyM7du3o3v37ti1axd69OiRqfME9J+JQ4YMQeXKlTFlyhR4eHggLCwMhw8fRv/+/QHoE+sOHTqgdu3amD59OgQCAdavX49WrVrhzJkzqFevnlXvEzGDEZIF69evZwDY8ePH2evXr1lMTAzbunUrK1SoEJPL5Sw2NpYxxlhoaCgDwEJDQxljjKlUKubr68uqVKnCUlJSDMc7cOAAA8CmTZtm2DZ69GiWmUu0efPmrHLlyhbv37NnDwPAli5dyhhjLDIykgmFQjZr1iyj/W7dusVEIpHR9ubNmzMAbMuWLYZt4eHhDAATCATswoULhu1HjhxhANj69esN27p3784kEgl79OiRYduzZ8+Yq6sra9asmWEb/7q2adOG6XQ6w/Zx48YxoVDI3r9/b3Quly9ftni+r1+/ZgDY9OnTTe4bNGgQA8AmT55scp9CoTDZNmfOHMZxHIuKijJsS+/9+fh5p0+fzgCwIUOGGO3Xo0cPVqhQIcO/r169ygCwL7/80mi/4OBgi+fysRs3bjChUMiOHDli9NyvX7822q9Tp06sZMmSJo/nr9mKFSsypVJp2L506VIGgN26dSvDGGrUqMF8fX3Z27dvjeISCARs4MCBhm18bP369TN6/KtXr5hEImGdOnUyug6++eYbBoANGjTIsO2HH35gzs7O7MGDB0bHmDx5MhMKhSw6OpoxxtiTJ08YAObm5sZevXqV4TmMGzeOAWBhYWEZ7suY7a/x+Ph4JpVK2YQJE4yeZ968eUbXYlb+jleuXGm0765duxgAtmTJEsM2rVbLWrVqZfK33Lp1a1a1alWWmppq2KbT6VijRo1Y2bJlM32e79+/Z66urqx+/fpGn4n8cfn/li1blgUFBRkdS6FQsNKlS7O2bdsyknXUNEaypU2bNvDx8YG/vz/69u0LFxcX7NmzB8WLFze7/5UrV/Dq1Sv83//9n1E/kk6dOqFChQo4ePCg3WLlf/EnJiYCAHbv3g2dTofevXvjzZs3hluRIkVQtmxZhIaGmjy+b9++hn+XL18eHh4eqFixotGvXf7/Hz9+DEA/6ubo0aPo3r27Uam/aNGi6N+/P86ePYuEhASj5xo+fLhRU1PTpk2h1WoRFRUFAIZfuAcOHIBarc7yazJq1CiTbXK53PD/ycnJePPmDRo1agTGGMLCwrL8XAAwcuRIo383bdoUb9++NZz/4cOHAQD/93//Z7TfF198YfVzjBkzBh06dEC7du2yFevgwYON+s40bdoUwIf31ZLnz5/j+vXrCA4ONqq2VKtWDW3btsXff/9t8piPX5fjx49DpVLhiy++MLoOzDVR7dixA02bNoWnp6fRddymTRtotVqcPn3aaP9PP/3U0GyaHv494aua6bHHNe7m5oYOHTpg+/bthioZAGzbtg0NGjRAiRIlAGT+71gqlWLw4MFG2w4fPgyxWIxhw4YZtgkEAowePdpov7i4OJw4cQK9e/dGYmKi4bnevn2LoKAgREREGDXxW3Oex44dQ2JiIiZPnmzSt45/3PXr1xEREYH+/fvj7du3hudNTk5G69atcfr0aaNmUJI51DRGsmX58uUoV64cRCIRChcujPLly0MgsJxf83/85cuXN7mvQoUKOHv2rN1iTUpKAvDhgz0iIgKMMZQtW9bs/h+PLvPz8zPpp+Tu7g5/f3+TbYC+KQ0AXr9+DYVCYfacK1asCJ1Oh5iYGFSuXNmwnf+Q53l6ehods3nz5vj0008xc+ZMLF68GC1atED37t3Rv39/SKXSdF6FD0QikdnRPtHR0Zg2bRr2799veD7ex/2rMiu983Jzc0NUVBQEAgFKly5ttF+ZMmWsOv62bdtw7tw53L59O1txZhQroL+e+GsK0De/+vj4pHuNV6xYEUeOHDHpEP3x+fLH+Pja9PHxMcTBi4iIwM2bNy0mN69evTL698fPZYmbmxuADz8c0mOPaxzQN4/t3bsX58+fR6NGjfDo0SNcvXrVqLkvs3/HxYsXN+kcHhUVhaJFi5p0mv74unv48CEYY5g6dSqmTp1q9vlevXpl9EMwo/PkuxNUqVLF7PEA/TkCMNtUx4uPjze5Noh1KBEi2VKvXj2TTou5Ff/lyH+46XQ6cByHQ4cOQSgUmuz/cZ8Rc/uktz3tr9jMyuiY/IR/Fy5cwF9//YUjR45gyJAhWLhwIS5cuGASuzlSqdQkadVqtWjbti3i4uLw9ddfo0KFCnB2dsbTp08RHByc7V+d9nit0po0aRJ69eoFiUSCyMhIAPqOxAAQExMDlUqFYsWKWXWsjGJdsGCBYXoCAChZsqThOTMrbRUus3Q6Hdq2bYuvvvrK7P3lypXL0nNVqFABAHDr1i3UqFEjy/FZYs210KVLFzg5OWH79u1o1KgRtm/fDoFAgF69ehn2yezfcXZfawCYOHEigoKCzO7zcfJki2uef9758+dbfC+s+Zsn5lEiRHJUyZIlAQD37983jDDj3b9/33A/gCyPEjNHq9Viy5YtcHJyMozKCQwMBGMMpUuXNvmysCUfHx84OTnh/v37JveFh4dDIBCYVJWs1aBBAzRo0ACzZs3Cli1b8Nlnn2Hr1q0YOnRoll6/W7du4cGDB9i4cSMGDhxo2G5uQkpbvj+8kiVLQqfT4cmTJ0a/8K0drRUTE4MtW7aYjLYBgFq1aqF69eq4fv06gOzHP3DgQMO1BHz4gk17jX8sPDwc3t7eGQ6P548RERFh1NT0+vVrkypdYGAgkpKS0KZNm6ydiAUdOnSAUCjE5s2bM+wwba9r3NnZGZ07d8aOHTuwaNEibNu2DU2bNjVKZm3xd1yyZEmEhoaaDKX/+Lrj3wuxWGyz1zswMBCA/oeapconv4+bm5vN32dCw+dJDqtTpw58fX2xcuVKKJVKw/ZDhw7h3r17RqOD+C8L/hd9Vmm1WowZMwb37t3DmDFjDCX/Tz75BEKhEDNnzjT5dcYYw9u3b7P1vDyhUIh27dph3759RhWDly9fYsuWLWjSpIkhJmu9e/fOJGb+lyL/uvIf6Jl5/fhfr2mPzRgzOw+Prd6ftPhf2StWrDDavmzZMqsev2fPHpNbnz59AACbNm0ymvfG2dk5W019AQEBaNOmjeHWuHFjAPp+MTVq1MDGjRuNXpvbt2/j6NGj6NixY4bHbtOmDcRiMZYtW2b0XpibCbt37944f/48jhw5YnLf+/fvodFoMn9yAPz9/TFs2DAcPXrU7Ouv0+mwcOFCxMbG2uUa5/Xp0wfPnj3D2rVrcePGDcP7ybPF33FQUBDUajXWrFljdH78FBE8X19ftGjRAqtWrcLz589NjmNuWHxG2rVrB1dXV8yZMwepqakm8QNA7dq1ERgYiAULFhg1x2bneckHVBEiOUosFuOnn37C4MGD0bx5c/Tr188wfL5UqVIYN26cYd/atWsD0Hd+DQoKglAoNOqsbE58fDw2b94MQD+smJ9Z+tGjR+jbt6/RUOrAwED8+OOPmDJlCiIjI9G9e3e4urriyZMn2LNnD4YPH46JEyfa5Lx//PFHHDt2DE2aNMH//d//QSQSYdWqVVAqlZg3b16mj7dx40asWLECPXr0QGBgIBITE7FmzRq4ubkZvmjlcjkqVaqEbdu2oVy5cvDy8kKVKlXS7YtQoUIFBAYGYuLEiXj69Cnc3Nywa9cukyoEkLX3JyO1a9fGp59+iiVLluDt27eG4fMPHjwAkHEVp3v37ibb+ApQhw4d4O3tbfRc27Ztw/jx41G3bl24uLigS5cu2YqfN3/+fHTo0AENGzZESEiIYfi8u7u7VXMh+fj4YOLEiZgzZw46d+6Mjh07IiwsDIcOHTI6B0DfHLh//3507twZwcHBqF27NpKTk3Hr1i3s3LkTkZGRJo+x1sKFC/Ho0SOMGTMGu3fvRufOneHp6Yno6Gjs2LED4eHhhvfc1tc4j5/nauLEiRAKhfj000+N7rfF33H37t1Rr149TJgwAQ8fPkSFChWwf/9+xMXFATC+7pYvX44mTZqgatWqGDZsGAICAvDy5UucP38esbGxuHHjRqbOz83NDYsXL8bQoUNRt25dw3xSN27cgEKhwMaNGyEQCLB27Vp06NABlStXxuDBg1G8eHE8ffoUoaGhcHNzw19//ZWp5yVp5OQQNZJ/8END0xu+zZjp8Hnetm3bWM2aNZlUKmVeXl7ss88+Mwy552k0GvbFF18wHx8fxnFchkPp+aGx/M3FxYWVLVuWDRgwgB09etTi43bt2sWaNGnCnJ2dmbOzM6tQoQIbPXo0u3//vtGxzQ3NL1myJOvUqZPJdgBs9OjRRtuuXbvGgoKCmIuLC3NycmItW7Zk586dM9rH0uv68et47do11q9fP1aiRAkmlUqZr68v69y5M7ty5YrR486dO8dq167NJBKJ0fDzQYMGMWdnZ7Ovx927d1mbNm2Yi4sL8/b2ZsOGDWM3btwwGUac3vuT9rkYszyEnT/fJ0+eGLYlJyez0aNHMy8vL+bi4sK6d+/O7t+/zwCwuXPnmo05PZaeOykpifXv3595eHgwAIah9PxrvWPHDqP9+eHnaV+D9Bw/fpw1btyYyeVy5ubmxrp06cLu3r1rVWyM6Ydvz5w5kxUtWpTJ5XLWokULdvv2bVayZEmj4fOMMZaYmMimTJnCypQpwyQSCfP29maNGjViCxYsYCqVyij++fPnWxU/T6PRsLVr17KmTZsyd3d3JhaLWcmSJdngwYNNhtbb8hpP67PPPjMMRbckO3/HjOmnm+jfvz9zdXVl7u7uLDg4mP37778MANu6davRvo8ePWIDBw5kRYoUYWKxmBUvXpx17tyZ7dy5M8vnuX//ftaoUSPD9VKvXj32559/Gu0TFhbGPvnkE1aoUCEmlUpZyZIlWe/evdk///xj8XUhGeMYs1EvRUIIsZPr16+jZs2a2Lx5s2EGbELsbe/evejRowfOnj1raPok+Q/1ESKE5Cppl/LgLVmyBAKBAM2aNXNARKQg+Pi602q1WLZsGdzc3FCrVi0HRUVyAvURIoTkKvPmzcPVq1fRsmVLiEQiHDp0CIcOHcLw4cOzPLqOkIx88cUXSElJQcOGDaFUKrF7926cO3cOs2fPztaQe5L7UdMYISRXOXbsGGbOnIm7d+8iKSkJJUqUwOeff45vv/0WIhH9diP2sWXLFixcuBAPHz5EamoqypQpg1GjRuF///ufo0MjdkaJECGEEEIKLOojRAghhJACixIhQgghhBRY1OCeAZ1Oh2fPnsHV1dUuSwoQQgghxPYYY0hMTESxYsXSXQycEqEMPHv2jEaqEEIIIXlUTEwM/Pz8LN5PiVAGXF1dAehfyKyulUMIIYSQnJWQkAB/f3/D97gllAhlgG8Oc3Nzo0SIEEIIyWMy6tZCnaUJIYQQUmBRIkQIIYSQAosSIUIIIYQUWJQIEUIIIaTAokSIEEIIIQUWJUKEEEIIKbAoESKEEEJIgUWJECGEEEIKLEqECCGEEFJgUSJECCGEkAKLEiFCCCGEFFiUCBFCCCGkwKJEiBBCCMkhsbGxCA0NRWxsrN23M8aMtu/atQu3bt1CcnIytFqtjc8s76LV5wkhhJBsiI2NRUREBMqWLQs/Pz+L29atW4fhw4dDp9NBIBBg9erVCAkJMdk+ZcoUNGnSBH/++Sd+//13Q0JTrlw5+Pv7Q6FQ4Pz584bnl8vlEAgEUCqV0Gg0hu2urq745JNPEBUVhVOnThklRgAgEAggkUggEolQs2ZNnD592nBfnz59IBAIEBgYiDJlyhj+W6RIkQxXc89rOPbxK5MPKZVK1K9fHzdu3EBYWBhq1Khh9WMTEhLg7u6O+Ph4uLm52S9IQgghuZq1yQ0Ao20//PADAgMD0a9fP5NkxNvbG2/evMnxczEnICAAtWvXhr+/P9asWYPExESTfZycnNCiRQscPHjQsO3y5csoWrSo4TXJLaz+/mYFwJgxY1iHDh0YABYWFpapx8bHxzMALD4+3j7BEUIIyVViYmLYiRMnWExMjGHb2rVrmUAgYACYQCBga9euZTExMYZtOXUrV66c2e0tW7Y0u71z586ZOo65m0wmY25ubozjOAaAtWnTxuj1Kly4MAPAGjZsyNauXcsSExNz+i0zy9rv73yfCP3999+sQoUK7M6dO5QIEUIIMbA24Xn06JHZhMfDw8PqZIJPIj6+Va1a1eQ+gUDADh8+bPKcQqGQXbp0ySbbt23bZjaeQYMGWYyVv3l6erIePXqwBQsWsIiICFaxYkUmFAoN9zs7O7OQkBB27tw5ptPpHPb+UiLEGHvx4gUrXrw4u3z5Mnvy5IlViVBqaiqLj4833GJiYigRIoSQPCy7CU92b0KhkK1cuZKNHTvWcHyhUMjWrl1riIVPJHJqu7lqVnoJUq1atcxuFwgErGPHjmzdunVs1qxZrGzZskb3f/vttzn5Vhsp8ImQTqdj7du3Zz/88ANjjFmdCE2fPt3sm02JECGE5G7WJjzR0dFmE57MJEECgYBVrFiR+fr6GiooAoGATZ48mS1ZssRsQsLHGBoaahSjo7ZbmyAJBAKTbRzHsQYNGhhtK1y4MPvqq6/YH3/8wQYOHMjkcjm7cOGC4fnu3r3LDh06xDQaTSbf2azJt4nQ119/neEFeu/ePbZ06VLWuHFjwwtOFSFCCMkfstuHx9nZOVMVHZFIxPz9/Q0JD5808M0+5pIMSwlJbmNNgjRx4kSzr0toaCh78OABGzVqFPP09DS6r0WLFmzNmjVMoVAYjjt8+HAGgPn7+7MjR47Y/dzybSL06tUrdu/evXRvSqWSdevWjQkEAiYUCg03/k0dOHCg1c9HfYQIIcQxspPwcBzHKlWqlKkKz6VLl1ibNm2Mts2YMYOpVCpDPHkhubGFtOdqqRktJibG6P3gOI5VqVLFqI9RoUKF2I4dOxhjjH377bfMy8vLULCwt3ybCFkrKiqK3bp1y3A7cuQIA8B27tyZqYuYEiFSUJj70klvOyH2ZG3CIxAI2NChQzNV4fn4Zm3zVUGWmX5GFy5cYDNmzGD+/v6G7SEhISwpKYmlpKSww4cP50jMBT4R+pi1TWMfo0SI5AXmkpXMJDbmvnTS257Z46e3nZCPrw1LX7BbtmyxOrnhOI7NmzePTZs2zahC0bp1a7Zt2zYWHR1NCU8mffx6nThxwmKTGWOMqVQqNmXKFMPrX65cOXb16tUci5cSoY9QIkTyA2uTGHPbkpOT2U8//WRUxh4+fDibN2+eyXBZjuNY3bp1zX7IlStXjhUpUsRoW9GiRVnr1q1Z586djY5Vvnx51q5dO1a9enWjY3/++ecsMTGR6XQ6FhMTw5YtW8b++usvdufOHRYVFcXevHnDHj58SIlTAWDuWrX0BSuXyy02a328bfjw4YbnuHv3Ljt8+DBdSzaWXpMZf/+JEyfY1q1bWbFixRgAJhaL2cKFC5lWq7V7fNZ+fxeImaWzg2aWJo5g7Qy2QUFBKFmyJHQ6nYMjtq+JEydi6NChYIxh3rx58PLyQmBgILy9vQ03Hx8f+Pr6QiCgJRRzq4+v69jYWJPrVygU4rvvvsP3339vMgszoF9OIjU11eQ+X19fNG7cGOXKlUOnTp3QpEmTfLcURG60bt06jBgxAlqtFkKhEKtWrTK7bMiiRYsQGhqKffv2AQDatWuHjRs3okiRInaLjWaWthGqCBF7srbCExkZafZXr7u7e7b6RqR342eL/fhWoUIFs9stdUz9eDSJPW+lS5dm79+/N7yWixYtYsuXL2f79+9nYWFh7O3btw6d4K0g+fjazkzlx9yN4zg2cuRIptFoWExMDNu7dy9r0KABmzNnDrt27VqOVBiIeR83mVmqFEVHR7Nff/2VyWQyBoD5+PiwgwcP2i0uahqzEUqEiC1kd7ivRCLJdpLQoEEDkyYwoVDI/v33X/bLL79keyK2zM5su2nTJrNxdunSJdvnWqhQIVa3bl0mEolM7pPL5ax8+fIsJCTE6D169+4dJUk28vG1PW/ePLOdnENCQsy+f0WKFDGbiPft29fRp0askFHfoTt37rCqVasato8dO5alpKTYPA5KhGyEEiGSXZkZ7luxYsVsJQB8sjJmzJhMzWDLWPYnYsvs9swkVJYmuuvfv3+GywFYunl7e7Np06axI0eOsNjYWObl5cU8PT1Zw4YN2eDBg9m8efPY/v372YMHD5harbb3ZZJvWBrZZc17wnGcUdIvEolYp06d2IoVK9ijR48cfWrEShn1HWKMsZSUFDZmzBjD/fv27bN5HJQI2QglQiQzrB39snbtWqu/sAUCAfvnn39Y0aJFjba3atWKrV69mp0/f57du3fP6gndbDVSxl4z25rbbq6ikN5yAOY6gPOPyUrS1L59e6Nzefz4MVWP/vPxNZ+Z5q4aNWqw9u3bs5s3bxqujX79+rHGjRuzFStWsNevXzv47EhWpfejK60DBw6wcePG2SUGSoRshBIhYq3s9oEAwKpVq8akUqlJ0sQYY2fOnGGnT5/OdyOprE2cbLFe0pYtWzK1tIJQKGQCgYDVqFGDzZ49m928eZO9f/+ecRzHvLy8WMeOHdn333/Pjh49atQ3qaDIzorsTk5Ohv+/efOm4Zj85IUk73P09ASUCNkIJULEHGsrP6GhoRabbvjhpOZunp6eLCgoKEfn3MgLsrNeUnoJ0ubNm61uvvH09DS7L8dxrHLlymzdunUOfIXsx9pr/quvvmLe3t5G24sXL85Gjhxp0qG+aNGibPz48SwyMtLBZ0fyI0qEbIQSIfKx7M570r9/f/bkyRPGGGOTJk0ybK9evTr79ttv2fnz53NsUcL8wl4J0ty5c9Ptg8RxHHN3dzcavffrr78aYoiNjWWzZs1i4eHhjnhZbCYr1U5nZ2fWpUsX9tdff7Hdu3cbVYKCg4PZ8ePH6TondkWJkI1QIkTSykwHX/5Wvnx59sMPP7Dt27ezCRMmGC02GB8fz1asWMGio6MdeFb5l70SpF69epn02QLAPDw8WM+ePdnWrVtZfHw8+/nnnw33VapUiU2dOpWFhYXlqf5Fmb3m69aty9avX88eP35sOIZarWbVqlVjEydOZC9fvnTg2ZDcIKdmmadEyEYoESrYrO0IOmbMGJPhvlWrVmVHjhxhb9++ZTNmzGAeHh4MAKtTp06e+iLMj7KbIKWX+PI3qVTKmjZtyurXr28yjD8gICDXJgXWXvM9e/Y0Oa8pU6awkydPsubNm7MSJUowpVJpOC5Vfwhj6S/bY2uUCNkIJUIFV1Y6grq4uLBPP/2UnT59mr17945Nnz7dqNmkUqVK7M8//6TJ33IpaxMkS8nBwoULWVBQkMn2ChUqsAEDBrCOHTsaJpMTi8VGHaxzQ6KQlWu+bt26bPHixWz37t2sVatWhu0SiYSdOXPG0adEchFrhtXbEiVCNkKJUMGU3h9s2i/GtL/wf/nlF5aQkMAYY2z16tXMzc3NcH/lypXZtm3bKAHKo6ydOTejSpFUKmW9e/dmM2bMYHPnzjV6jg4dOrCRI0ey58+fO+IU073mV65cadJXqm7duuzw4cPs/PnzrF27dobtYrGYjRo1Kl+NbCS2kdFEi7ZGiZCNUCJUMFn6g/3nn3/Y8uXLmbOzs+GL4rvvvjOZFXXnzp0MAKtSpQrbvn07JUD5UGYqRWPHjmU1atQw2la2bFm2adMmptPp2JUrVwzbnZyc2NSpUw1JdU6xFPv27dtZw4YNDf+uWLEi27dvH9PpdOzUqVOG7SKRiA0fPpxGgBGLqCKUR1EiVDBYMzRYIBCw2rVrG/5dv359w/wnWq3WqMOzVqtlBw4coAQon7O2UhQTE8PWrFljtlrUvHlzdufOHXbq1ClWv359w3YfHx+2bNkyo342to7dmmu+UKFCDNCva/fHH38YXdM6nY4NGzaMDRgwgGZ+JlaxdqJFW6BEyEYoEcr/LHXeS/sHy3Gc4f9dXFzYzz//bOjTodFoWEhICPPx8WH37t1z5KmQXMDaTtccxxn6C4lEIvb111+zxMREtnPnTla2bFnDfoGBgSwiIsLmMVpzzfPNYdWqVTPEoNFojCqglOyTzMqpiRYpEbIRSoTyt4xKtQcPHmQlSpQw3Ne5c2ejyo9KpWL9+vUzfKFs3brVUadCcpGPP+gtNTv9+eefrGvXroZ/lyxZku3bt4+pVCq2YsUKVrhwYVa2bFmbzrac0TV/584d1qRJE8N9AwcOZMnJyYwxxpRKJevduzfr0qULzQBNcj1KhGyEEqH8Lb3Oe9euXTN0ePb19WXbtm0zGvaemprKunXrZvhFv2PHDgeeCcnNMko+1q1bZzT9QpcuXdiTJ09YYmIiu337tuE4SqWShYSEsBs3bmQ5lvSu+Vu3bhmqURKJhK1cudJwzSclJRlGxInFYnb+/PnsvSiE2BklQjZCiVD+k7ZvhKUvqFOnTjEfHx8GgDVt2pS9ffvW6BjJycmGkTJSqZQdOHDAQWdD8or0Fpjlr0GO4wz/L5fL2ezZs5larTYcg5+gUSAQsIMHD2YpDkvX/NKlSw3rf/n7+7OLFy8aHhMXF2foMO3k5MSOHj2avReDkBxAiZCNUCKUv5jrG2FupXO+OaxmzZomi2nGx8ezpk2bMkC/jMA///zjoLMheY01nasFAgFr0KCB4d/Dhw83PD4iIoJ16tSJAfp1uuLi4rIUx8fXfNu2bQ3P17ZtW6NV3589e8aqVq3KAP06a+fOncvei0BIDqFEyEYoEcoZWq2WqVQqlpKSwpKSkuzSATO95gn+C+r69eusfPnyDAArV66c2Zl/ExISWIMGDZibmxv7999/bR4nKTgsNVOdOHGCrV271tBZOe1CrgqFwnCNDho0KMPnsLScAX/Nr1+/3vC83333ndHEjo8ePWIBAQGGxOvWrVs2O3dC7I0SIRuhRCjzFAoF+/fff9ny5cvZ0KFDWZ06dVjRokWZj48P8/LyYosWLTLse/XqVbOLWhYqVIj179/fpiX4jCbzio+PZ7Vq1TI0DURFRVk81rt379j169dtFhspmDLqO/T9998bml8vX75seNy5c+cMj0uvWTaj5QxUKpUhqfrqq69MHn/x4kXm7OzMAgICaHg8yXMoEbIRSoTS9/z5c3bo0CHDfDqMMaNJ1szdfvjhB8O+YWFh6e47a9Ysw77v379n586dy/JSBOl96SgUCta8eXMGgHl7e5usFh4TE2O0qjghtpLevCparZZ16dKFAWAlSpQwarKaOHEiA8CKFStmtonMmsnr+D5Hvr6+Fj/jTp8+zZ49e2bDMyYkZ1j7/S0CIZnw9OlTbN++HUePHkVYWBhevnwJAOjevTsaNGiA6OhoJCYmQi6Xw93dHW5ubnBycoJIJIJWq4VarcaRI0cQFxeHGjVqoEqVKoiJiYFUKoVQKIRIJALHcQgLC8OhQ4fQrVs3w3MfOnQI/fr1g5eXF9q1a4f27dujffv2KFy4sFWx+/n5YfXq1RgxYgS0Wi2EQiFWrVqFwoUL45NPPsGpU6fg5uaGI0eOoHz58obHPX78GK1bt0ZkZCSEQiGGDRtm2xeVFGghISEICgrCw4cPUaZMGfj5+RnuEwgE2LRpE+rWrYuHDx+iX79+OHz4MIRCIb7//nvs378f8fHxePjwIerWrWt03IiICOh0OqNtWq0WDx8+hJ+fH969e4cZM2YAAH744Qe4ubkBAI4dOwYvLy/Url0bANC0aVM7nj0huUAOJWZ5FlWE9F6/fm00zb4tb1KplNWpU4cNGzaMrVixgp0/f56lpqaaxLB8+XLDCu5pb7Vq1WKHDh2y+lzSdljVarXss88+YwCYTCZjp06dMtr33r17rFixYgwAK1OmTLrNZYTYUtq+Pbdu3TKM6Jo8ebJhn9u3b7M3b95YfHx6FaFx48YxQL8MDD8y7eDBg0wikTBvb29qCiN5HjWN2UhBTYTevHnDduzYwdasWcMmTJhgtLTExzeRSGRIFgD9EOBevXqxZcuWsYEDBxr6AHEcx8aMGcPWr1/PxowZw5o2bcpcXV3NHrN48eJs//79JnGp1Wp29uxZ9t133xnF5OzsnOk1jnQ6HRs9erThHD7ua5GUlGSYU6Vy5crUPEByjLm+PX/++afhet+9e7fVxzHX7PbgwQMmFosZAHbkyBHGGGOHDh1iEomEAWCffvqp3Zb2ICSnUCJkIwUlEYqJiWH79+9n8+bNY/Xq1bO6miMQCNiZM2esXok77S/SmJgYdvz4cXb69Gm2fft2NmXKFNahQwfm7e1t2L9Pnz5mR27xXrx4wdq0aWMy3wp/fHOjZXirV682JGhbtmwxuZ9Pkvz8/Iz6ZhBiT+lVcr788ksGgLm6uhot56LT6diWLVvYhAkTzB7v4+UMunfvzgCwjh07MsYYO3z4MJNKpYYkiGaNJvkBJUI2UhASoUmTJlmV8KxatcrsfYsWLcrU9tDQ0HRHsygUCvbVV18Zfsl6eXmxjRs3Gs3qnJa5ofYZjZZhjBnmY5k6darJfUeOHDHEe+zYsWy8uoRkTnqjG1UqFWvWrBkD9KvA8yvU375921B5zaiZmD++UChkd+/eZUeOHDEkQT169KAkiOQblAjZSH5OhOLi4ljLli2trv5s3749U5WfzG7/uHJz9epVVqNGDcM+QUFB7MmTJ+meU2pqKrty5YpVx/fz82MA2NmzZ02Os3DhQiYQCNgXX3yR7deZkMzIqG/P8+fPDU3RPXv2NPxA4KtFfn5+JpOA8jQajeFvavTo0ezcuXOGhV+7d+9OzWEkX6FEyEbyWyIUExPD/vjjDzZkyBDm7OycbgXI3AdxessEWLs9o/l80lKpVGzOnDmGX6zOzs5syZIlZofQ37lzh1WuXJlVrFgxw+O/fv3asJ3/Vf2xixcvGhabJCQnpTeknjH9PEJ8H5/58+czxvTLvgQGBjIAbOjQoWaP+9tvvzEAzN3dnb1+/ZolJiaypk2bsq5du1ISRPIdSoRsJD8lQvwvxrS3ChUqmExoKBQK2fz58y1+EJvrc5CZ7Rn94jXXtyc8PNywrAUA1qBBAxYREWH0PNHR0YZFUs2dU9rjHT9+nAFggYGBtnlxCbExS39PvOXLlxt+tHz11VcsPj6enT592nDt852geYmJiaxo0aIMAFuwYIFhe1JSEiVBJF+iRMhG8noipNPp2Pbt21nx4sXNVn2io6Mt/vrM6IM4O6xZgPLjvj1arZb9+uuvhpFmjRs3Njkuv1yAUCg0HMfcL+qFCxcyQN8xlPfs2TPWrl07dvfuXZufLyG2ptPp2IgRIwx/z97e3mzixIksODiYAfrZ0dN+bk2dOpUB+qUyfvzxRwdGTkjOoETIRvJqIhQTE8M2b97MypUrl26/H765yJ5JT3oxZqZSxLtz5w4D9EPeU1JSjO7T6XSGmXgrV67Mjh49avacPv/8cwaAff/994bHdejQwZBgWeqYTUhuotPp2F9//cV8fX2N/m7c3d0ZAPbHH38wxvTVUrlczgAYhsibGylJSH5CM0sXYGvXrsXw4cPBGDNs8/DwQHx8vNE2oVCIMmXKANDPupx2Rtuc8PFzZjQTLq9ixYrw9fXFq1evcP36dTRo0MBwH8dxWL16NSpXrow7d+7g7NmzaNu2rclzX79+HQBQo0YNAMCaNWtw6NAhSKVSrF69GhzH2fBMCbEPjuNQo0YNvHnzxmh7fHw8GjRogCZNmgAAvvnmG6SkpEAgEEClUqFDhw7o0aOHI0ImJNcRODoAYlvXrl3DsGHDjBIejuNw/fp1rFmzBkKhEAAMy0vkdPKTnrJly0IgML4k0yZrPI7jDMsJXL582eQ4RYoUwYoVKwAAs2bNwtWrV43uVyqVuHfvHgB9IvTo0SOMHz8eADBnzhxUqlTJNidESA4w9wMCAC5cuIDy5ctj5MiR2Lx5MwBAp9Ohffv22L17N2QyWU6HSkiuRIlQHhcbG4vQ0FDExMRg48aNaN68uck+jDE8efIEISEhiIyMRGhoKCIjIxESEuKAiC3j1wKzJlnjE6ErV66YPVafPn3Qu3dvNGvWDIUKFTK67+7du9BoNPD09ETRokUxaNAgJCcno3nz5hg7dqyNz4oQ+zL3A0IgEKBBgwZITU3FqlWrDNvbtWuHPXv2UBJESFo50EyXp+XmPkJpOxandzPXzyY3s6a/0oEDBxign1TOkqSkJLOTLfJDiFu2bMnmzp3LAP1MvRnNUURIbmVu8IFOp2Nr1641jCJr2rQpUygUjg6VkBxDfYTyudjYWAwfPtykJO7s7Iw5c+Zg3LhxRius56YmsIxY01+JrwiFh4cjMTERrq6uJvs4Ozsb/Ts5ORnOzs64ceMGAKB69er4+++/AQBLly5FqVKlbBA9ITnP0gr2ISEhEAgEuHv3Lr7//nvI5XIHR0pI7kOJUB5lqV/AqlWr8Nlnn6FHjx4mH4r5ia+vL0qUKIHo6GhcvXoVLVq0sLhvcnIyvvrqKxw/fhxhYWFGHaXnz5+PXbt2oXfv3jkTOCF2YukHxODBgx0QDSF5B/URyqMePXpksk0oFBr6CPn5+aFFixb5MgnipddhOi21Wo19+/bhwYMHmDx5slEiJBKJ0KdPHxolRgghBRQlQnnQmTNnMGrUKKNtebEJLLv4ROjUqVMIDQ1FbGys2f08PDywbt06AMCyZcsQHx8PgUCA0qVL51ishBBCcidKhPIIfnTY8ePH0bVrV2g0GgBAv379cPz48Vw5CszW+NeAT3j4ROjgwYNo1aoVSpYsaUh4PhYUFIThw4cb/q3T6TBv3jz7B00IISRX4xhLM+EMMZGQkAB3d3fEx8fDzc3NITGsW7fOpGN048aNMXbsWPTs2bNANOukfQ0EAgFWr16NBg0aoEqVKkb7CYVCREZGmq2MJSYmokSJEnj//j2cnJzw7NkzuLu759QpEGITsbGxiIiIQNmyZQtUBZiQzLL2+5sqQrmcpdFhv/76K3r16lUgkqCPXwOdTocRI0bgzp07JvvyM1Gb4+rqaphJevjw4ZQEkTxn3bp1KFmyZIYVUEKI9SgRyuUsjQ57+/atA6JxDEtLb5hLAs3NRJ1WVFQUAKBr1662DZIQO7P0g8BS3zhCiHUoEcrlzM1tk9GXfX5jaemNhg0bom/fvkbb0usw/v79ezx58gSAfg4hQvKS9NbiI4RkXb5PhA4ePIj69etDLpfD09MT3bt3d3RIVmOM4fvvvzfaVhBHh6W39MYXX3wBAPD09DQsI2LJzZs3AQAlSpSAl5eX/QMnxIasXYuPEJI5+XpCxV27dmHYsGGYPXs2WrVqBY1Gg9u3bzs6LKtNnz4dGzZsgFAohFarhYuLC+7du1egkiCepZlza9SoAaFQiHfv3pl8SXws7YzShOQ1/A+CESNG5NlZ4wnJjfJtIqTRaDB27FjMnz/fqEqQV1YW1+l0+OmnnwAAU6ZMwY8//ggPD48C/aFnbuZcJycnVK5cGTdv3sTly5dRvHhxi49PO5EiIXmRpR8EhJCsy7dNY9euXcPTp08hEAhQs2ZNFC1aFB06dMiwIqRUKpGQkGB0c4TU1FSoVCoAQJMmTQCYrp1F9KydYZoSIZIfFIRZ4wnJSfk2EXr8+DEAYMaMGfjuu+9w4MABeHp6okWLFoiLi7P4uDlz5sDd3d1w8/f3z6mQAXyYNDDtEhpqtRqAvvpBTFmTCKnVasNwe2oaI4QQwstzidDkyZPBcVy6t/DwcMPoim+//RaffvopateujfXr14PjOOzYscPi8adMmYL4+HjDLSYmJqdOzWiOEL5qIRaL0bZtW8TExGDv3r05FktewidCV65cgaX5Qe/fvw+lUglXV1daWoMQQohBnusjNGHCBAQHB6e7T0BAAJ4/fw7AuE+QVCpFQEAAoqOjLT5WKpVCKpXaJNbMMDdHCABIJBJIpVIqg6ejSpUqkEgkePfuHR4/fozAwECTffhmserVq2fYqZoQQkjBkecSIR8fH/j4+GS4X+3atSGVSnH//n1DHxu1Wo3IyEiULFnS3mFmmqWJE8VisQOiyVskEglq1KiBS5cu4fLly2YTIRoxRgghxJx8+9PYzc0NI0eOxPTp03H06FHcv3/fsGJ7r169HBydKXNzhAD6DtLHjh3D+PHjsWvXLgdEljdk1E+IOkoTQggxJ98mQgAwf/589O3bF59//jnq1q2LqKgonDhxAp6eno4OzcTHkwbySZGLiwvOnz+PxYsX49ixY44MMVdLLxFijFEiRAghxKx8nQiJxWIsWLAAL1++REJCAo4dO4bKlSs7OiyLQkJCEBkZidDQUGzatAkAIJfLkZycDICGz6eHT4SuXbsGrVZrdN/z58/x5s0bCASCXP3+E0IIyXn5OhHKi/g5QlxcXAAAMpnMkAjR8HnLypcvD2dnZyQnJ+PevXtG9/HVoAoVKkAulzsgOkIIIbkVJUK5VGpqKgB9RUihUACgilB6hEIhateuDUA/jD4tahYjhBBiCSVCuVRKSgoA44oQJULps9RPiEaMEUIIsYQSoVwqbUWIEiHrWEqEqCJECCHEEkqEcqm0FSFqGrMOnwjduHHDsE5bcnIyIiIiAFBFiBBCiKk8N6FiQZG2IvTHH3/g3bt3KFasmIOjyt1Kly4NLy8vxMXF4ebNm6hTpw5u3boFxhiKFCmCwoULOzpEQgghuQxVhHIpviIkl8tRtGhRVKpUCR4eHo4NKpfjOA516tQB8KF5jJrFCCGEpIcSoVwqbdMYsV7aBVgBSoQIIYSkj5rGcqm0TWMzZswAYwz/+9//rFpnrSD7uMM0jRgjhBCSHkqEcqm0FaHZs2cjOTkZgwYNokQoA3widOfOHSQmJuLmzZsAqCJECCHEPGoay6X4ilDaUWM0s3TGihUrhmLFikGn02HHjh1QKBSQy+UoW7aso0MjhBCSC1EilEvxFSGhUAjGGAAaPm8tvsP0mjVrAABVq1Y1LGZLCCGEpEWJUC7FV4T4VegBqghZi28eu3DhAgBqFiOEEGIZJUK5FF8R4slkMqpqWIlPhHiUCBFCCLGEEqFciq8IcRwHgKpBmcE3jfFoxBghhBBLKBHKpfiKkE6nA0D9gzKjUKFCCAgIAKBPJKtWrergiAghhORWNHw+l+IToYoVK+LWrVvQarUOjihvqVOnDh4/fowyZcrA1dXV0eEQQgjJpagilEvxTWMeHh6oUqUKNe9kUpMmTQAA9evXd3AkhBBCcjOqCOVStMRG9owcORIymQydO3d2dCiEEEJyMUqEcim+IvTo0SPs3r0bVapUwaeffurgqPIOsViMYcOGOToMQgghuRw1jeVSfEXozp07mDFjBrZs2eLgiAghhJD8hxKhXEir1UKtVgOA4b80aowQQgixPUqEciG+WQwAVCoVAEqECCGEEHugRCgXokSIEEIIyRmUCOVCfP8gsVhs+H+aWZoQQgixPUqEciE++ZHL5UhOTgZAFSFCCCHEHigRyoX4pjFKhAghhBD7okQoF0o7meJPP/2Ec+fO4ZNPPnFwVIQQQkj+QxMq5kJpK0IBAQGGBUQJIYQQYltUEcqFaHkNQgghJGdQRSgXSlsRWrlyJRQKBfr06YPixYs7ODJCCCEkf6FEKBdKWxGaO3cuoqKi0LhxY0qECCGEEBujprFcKG1FSKFQAKBRY4QQQog9ZKsipFar8eLFCygUCvj4+MDLy8tWcRVoaStC/PB5mlCREEIIsb1MV4QSExPx66+/onnz5nBzc0OpUqVQsWJF+Pj4oGTJkhg2bBguX75sj1gLjLSJEFWECCGEEPvJVCK0aNEilCpVCuvXr0ebNm2wd+9eXL9+HQ8ePMD58+cxffp0aDQatGvXDu3bt0dERIS94s7X+KYxiURi2EaJECGEEGJ7mWoau3z5Mk6fPo3KlSubvb9evXoYMmQIVq5cifXr1+PMmTMoW7asTQItSPiKkFAoNGyjpjFCCCHE9jKVCP35558A9H2D2rdvj5UrV5pNdKRSKUaOHGmbCAsgviLEJ0JyuRwCAfVrJ4QQQmwtS52lxWIxbt68aetYyH/4ilChQoVw5swZqFQqB0dECCGE5E9ZLjMMGDAA69ats2Us5D98RcjV1RVNmjRBq1atHBwRIYQQkj9lefi8RqPBb7/9huPHj6N27domnXkXLVqU7eAKKlpigxBCCMkZWU6Ebt++jVq1agEAHjx4YHQfx3HZi6qA4ytCiYmJWLp0KUqVKoVu3bo5OCpCCCEk/8lyIhQaGmrLOEgafEXo5cuXmD59Opo1a0aJECGEEGIHNBQpF+ITIZ1OB4DmECKEEELsJVuJ0Pv377Fw4UIMHToUQ4cOxaJFixAfH2+r2LLtwYMH6NatG7y9veHm5oYmTZrkiUoW3zRGiRAhhBBiX1lOhK5cuYLAwEAsXrwYcXFxiIuLw+LFixEYGIhr167ZMsYs69y5MzQaDU6cOIGrV6+ievXq6Ny5M168eOHo0NLFV4S0Wi0AmkyREEIIsZcsJ0Ljxo1D165dERkZid27d2P37t148uQJOnfujC+//NKGIWbNmzdvEBERgcmTJ6NatWooW7Ys5s6dC4VCgdu3bzs6vHTxFSE+EaKKECGEEGIf2aoIff311xCJPvS3FolE+Oqrr3DlyhWbBJcdhQoVQvny5bFp0yYkJydDo9Fg1apV8PX1Re3atR0dXrr4ipBarQZAiRAhhBBiL1keNebm5obo6GhUqFDBaHtMTAxcXV2zHVh2cRyH48ePo3v37nB1dYVAIICvry8OHz4MT09Pi49TKpVQKpWGfyckJOREuEb4ihCfCFHTGCGEEGIfWa4I9enTByEhIdi2bRtiYmIQExODrVu3YujQoejXr58tYzQyefJkcByX7i08PByMMYwePRq+vr44c+YMLl26hO7du6NLly54/vy5xePPmTMH7u7uhpu/v7/dzsUSviI0ZMgQ/P333+jfv3+Ox0AIIYQUBBxjjGXlgSqVCpMmTcLKlSuh0WgA6NcgGzVqFObOnQupVGrTQHmvX7/G27dv090nICAAZ86cQbt27fDu3Tu4ubkZ7itbtixCQkIwefJks481VxHy9/dHfHy80XHsSSKRQK1WIzY2FsWLF8+R5ySEEELyk4SEBLi7u2f4/Z3lpjGJRIKlS5dizpw5ePToEQAgMDDQ7s04Pj4+8PHxyXA/hUIBACartgsEAsOwdHOkUqndkjhraLVaQ5MYLbFBCCGE2FeWmsbUajVat26NiIgIODk5oWrVqqhatWqu6svSsGFDeHp6YtCgQbhx4wYePHiASZMm4cmTJ+jUqZOjw7OI7x8EAEePHsVvv/2GmJgYB0ZECCGE5F9ZSoTEYjFu3rxp61hsytvbG4cPH0ZSUhJatWqFOnXq4OzZs9i3bx+qV6/u6PAs4vsHAcCCBQsQEhKS619rQgghJK/KctPYgAEDsG7dOsydO9eW8dhUnTp1cOTIEUeHkSl8IiSRSAz/T8PnCSGEEPvIciKk0Wjw22+/4fjx46hdu7bJl/WiRYuyHVxBxDeNyeVyJCcnA6Dh84QQQoi9ZDkRun37NmrVqgVAv6ZXWhzHZS+qAoyvAslkMkMiRBUhQgghxD6ylAip1WoIBAKsXLkSZcuWtXVMBVraitDLly8BUCJECCGE2Eu+7SydV6WtCPH/T01jhBBCiH1keWZpvrM0sS2+IpR2LiOqCBFCCCH2QZ2lcxm+CiSXy3HgwAEoFArI5XIHR0UIIYTkT9RZOpfhK0JOTk65euJHQgghJD/IciIUGhpqyzjIf9L2ESKEEEKIfWW5jxCxDz4R4jgOGzduxN9//+3giAghhJD8K1uJ0JkzZzBgwAA0bNgQT58+BQD8/vvvOHv2rE2CK4j4pjGlUong4GBMmjTJwRERQggh+VeWE6Fdu3YhKCgIcrkcYWFhUCqVAID4+HjMnj3bZgEWNHxFSCDQvzU0dJ4QQgixnywnQj/++CNWrlyJNWvWQCwWG7Y3btwY165ds0lwBRFfEeITIRo6TwghhNhPlhOh+/fvo1mzZibb3d3d8f79++zEVKClXX0eoESIEEIIsacsJ0JFihTBw4cPTbafPXsWAQEB2QqqIOMrQjxqGiOEEELsJ8uJ0LBhwzB27FhcvHgRHMfh2bNn+OOPPzBx4kSMGjXKljEWKHxFiDEGgCpChBBCiD1leR6hyZMnQ6fToXXr1lAoFGjWrBmkUikmTpyIL774wpYxFih8RYgSIUIIIcT+spwIcRyHb7/9FpMmTcLDhw+RlJSESpUqwcXFxZbxFTh8Rah+/foYPHgwNTMSQgghdpTlRIgnkUhQqVIlW8RC8CERCgwMRN++fR0cDSGEEJK/0czSuQzfNEZLbBBCCCH2R4lQLsNXhB4/fozdu3cjOjrawRERQggh+RclQrkMXxHavXs3Pv30U1rclhBCCLEjSoRyGb4ipFKpANA8QoQQQog90aKruQxfEeITIRo+TwghhNgPLbqay/AVIf71pESIEEIIsR9adDWX4StC/H8pESKEEELshxZdzWX4ihD/X+ojRAghhNgPLbqai2g0Gmg0GgAfEiGqCBFCCCH2k+WZpflFV3/77TfDoqvnz5/HxIkTMXXqVFvGWGCkXXl+xYoV0Gg08PHxcWBEhBBCSP5Gi67mInwVCACCg4MhENDsBoQQQog9cYxf5jyLVCpVvl50NSEhAe7u7oiPj4ebm5tdnysmJgYlSpSAVCo1qg4RQgghJHOs/f7Ocslh6NChOHnypGHR1Xr16uW7JCin8RUhmUyGvXv30qzShBBCiJ1lORF6/fo12rdvD39/f0yaNAk3btywZVwFEl8FEovF6NGjB/r06ePgiAghhJD8LcuJ0L59+/D8+XNMnToVly9fRq1atVC5cmXMnj0bkZGRNgyx4OArQvy8TDRijBBCCLGvbPXG9fT0xPDhw3Hy5ElERUUhODgYv//+O8qUKWOr+AqUtBUhgBIhQgghxN5sMixJrVbjypUruHjxIiIjI1G4cGFbHLbA4StCIpF+MB8lQoQQQoh9ZSsRCg0NxbBhw1C4cGEEBwfDzc0NBw4cQGxsrK3iK1A+ToRoVmlCCCHEvrI8j1Dx4sURFxeH9u3bY/Xq1ejSpQukUqktYytw+KYxfv4gqggRQggh9pXlRGjGjBno1asXPDw8bBhOwcZXhCgRIoQQQnJGtpbYILbFV4SKFy+OL774AqVLl3ZwRIQQQkj+lqlEaPz48fjhhx/g7OyM8ePHp7vvokWLshVYQcRXhAoXLoyRI0c6OBpCCCEk/8tUIhQWFga1Wm34f0s4jsteVAUUXxGSy+UOjoQQQggpGDKVCKVd8mHjxo3w8/MzWRiUMYaYmBjbRFfA8BWhpKQk/PPPPyhfvjz8/PwcHBUhhBCSf2V5+Hzp0qXx5s0bk+1xcXHUtyWL+IrQjRs30KZNG2zdutXBERFCCCH5W5YTIUuL1iclJUEmk2U5oIKMrwjpdDoANGqMEEIIsbdMjxrjO0lzHIdp06YZTfqn1Wpx8eJF1KhRw2YBFiR8IqTVagHQhIqEEEKIvWW6IhQWFoawsDAwxnDr1i3Dv8PCwhAeHo7q1atjw4YNdgjV2KxZs9CoUSM4OTlZnMsoOjoanTp1gpOTE3x9fTFp0iRoNBq7x5ZVfNMYHyNVhAghhBD7ynRFiO8wPXjwYCxduhRubm42D8oaKpUKvXr1QsOGDbFu3TqT+7VaLTp16oQiRYrg3LlzeP78OQYOHAixWIzZs2c7IOKM8RUhSoQIIYSQnJHlCRXXr19vyzgybebMmQBgsfp09OhR3L17F8ePH0fhwoVRo0YN/PDDD/j6668xY8YMSCSSHIzWOnxFiJ+igBIhQgghxL6ynAjx7t69i+joaKhUKqPtXbt2ze6hs+X8+fOoWrUqChcubNgWFBSEUaNG4c6dO6hZs6bZxymVSiiVSsO/ExIS7B4rj68I8a8l9REihBBC7CvLidDjx4/Ro0cP3Lp1CxzHGUaR8ZMp8h1+HeXFixdGSRAAw79fvHhh8XFz5swxVJtyGl8R+uyzz+Dv74+SJUs6JA5CCCGkoMjy8PmxY8eidOnSePXqFZycnHDnzh2cPn0aderUwcmTJ7N0zMmTJ4PjuHRv4eHhWQ3ZKlOmTEF8fLzhlpOTQ/IVoS5dumDSpEnw8fHJsecmhBBCCqIsV4TOnz+PEydOwNvbGwKBAAKBAE2aNMGcOXMwZsyYdJfgsGTChAkIDg5Od5+AgACrjlWkSBFcunTJaNvLly8N91kilUohlUqteg5boyU2CCGEkJyV5URIq9XC1dUVAODt7Y1nz56hfPnyKFmyJO7fv5+lY/r4+NisCtKwYUPMmjULr169gq+vLwDg2LFjcHNzQ6VKlWzyHLbGV4Tu3LkDqVSKOnXq0LpthBBCiB1lORGqUqUKbty4gdKlS6N+/fqYN28eJBIJVq9ebXXVJjuio6MRFxeH6OhoaLVaXL9+HQBQpkwZuLi4oF27dqhUqRI+//xzzJs3Dy9evMB3332H0aNHO6zikxE+ERo5ciQ4jnN4PytCCCEkv8tyIvTdd98hOTkZAPD999+jc+fOaNq0KQoVKoRt27bZLEBLpk2bho0bNxr+zY8CCw0NRYsWLSAUCnHgwAGMGjUKDRs2hLOzMwYNGoTvv//e7rFlFd80BuhHjFE1iBBCCLEvjllaNCwL4uLi4Onpma++wBMSEuDu7o74+Hi7Tx4pEokMVSBfX19DnyZCCCGEZI6139/ZnkcoLS8vL1serkDRaDRGTWE0mSIhhBBif5lKhPgFV62xaNGiTAdTkPH9g3iUCBFCCCH2l6lEyNoh8fmpaSynpO0fBNCs0oQQQkhOyFQixC+4SmyPrwiJxWKo1WqqCBFCCCE5IMszSwPAmTNnMGDAADRq1AhPnz4FAPz+++84e/asTYIrSNJOpjh37lyEhIQ4OCJCCCEk/8tyIrRr1y4EBQVBLpfj2rVrhoVK4+PjMXv2bJsFWFDwFSEnJyd8/fXX+OyzzxwcESGEEJL/ZTkR+vHHH7Fy5UqsWbMGYrHYsL1x48a4du2aTYIrSPhEiJbXIIQQQnJOlhOh+/fvo1mzZibb3d3d8f79++zEVCDxTWNCoRBXr17FixcvHBwRIYQQkv9lOREqUqQIHj58aLL97NmzObLERn7DV4QSEhJQp04dLFmyxLEBEUIIIQVAlhOhYcOGYezYsbh48SI4jsOzZ8/wxx9/YOLEiRg1apQtYywQ+IoQP/UAjRojhBBC7C/LM0tPnjwZOp0OrVu3hkKhQLNmzSCVSjFx4kR88cUXtoyxQOArQnwiRPMIEUIIIfaX5USI4zh8++23mDRpEh4+fIikpCRUqlQJLi4utoyvwPh4QkWqCBFCCCH2l+21xiQSCSpVqmS0LSUlhUY/ZRJfEeLXwKVEiBBCCLG/bE2o+DGlUomFCxeidOnStjxsgcBXhHQ6HQBqGiOEEEJyQqYTIaVSiSlTpqBOnTpo1KgR9u7dCwBYv349SpcujSVLlmDcuHG2jjPf4ytCfCJEFSFCCCHE/jLdNDZt2jSsWrUKbdq0wblz59CrVy8MHjwYFy5cwKJFi9CrVy8IhUJ7xJqv8YlQtWrV0LBhQ5QtW9bBERFCCCH5X6YToR07dmDTpk3o2rUrbt++jWrVqkGj0eDGjRu06nw28E1jdevWxaxZsxwcDSGEEFIwZLppLDY2FrVr1wYAVKlSBVKpFOPGjaMkKJtoiQ1CCCEk52U6EdJqtZBIJIZ/i0QiGjJvA3xFKC4uDg8ePIBGo3FwRIQQQkj+l+mmMcYYgoODIZVKAei/wEeOHGnSuXf37t22ibCA4CtCy5Ytw7Jly/D27Vt4eXk5OCpCCCEkf8t0IjRo0CCjfw8YMMBmwRRkH0+oSMPnCSGEEPvLdCK0fv16e8RR4PEVIQAQCASGihshhBBC7CdTfYSio6MzdfCnT59mav+CLG1FyNnZmTqfE0IIITkgU4lQ3bp1MWLECFy+fNniPvHx8VizZg2qVKmCXbt2ZTvAgiJtRYiaxQghhJCckammsbt372LWrFlo27YtZDIZateujWLFikEmk+Hdu3e4e/cu7ty5g1q1amHevHno2LGjveLOd9ImQjSrNCGEEJIzOMav8pkJKSkpOHjwIM6ePYuoqCikpKTA29sbNWvWRFBQEKpUqWKPWB0iISEB7u7uiI+Ph5ubm92ep1y5coiIiAAAVK1aFTdv3rTbcxFCCCH5nbXf31lafV4ul6Nnz57o2bNnlgMkxviKUHBwcL5KJAkhhJDcLEuJELE9vrP0xIkTUblyZQdHQwghhBQMmZ5ZmtgHLbFBCCGE5DxKhHIJviL0/PlzvHv3zsHREEIIIQWDTZrG3r17h6NHjxrmDSpWrBiCgoLg6elpi8Pne2q1GlqtFgDQpEkTDBkyBOvWrXNwVIQQQkj+l+2K0Lp169CwYUNcvHgROp0OOp0OFy9eRKNGjejL3EofL69Bw+cJIYSQnJHtitC8efNw7do1ky/vH374AbVq1UJISEh2nyLfSzuHEECJECGEEJJTsl0R4jgOiYmJJtsTExNpmQgr8YmQUCgEQDNLE0IIITkl2xWhBQsWoHnz5qhSpQqKFy8OAIiNjcWdO3ewcOHCbAdYEPBNYwKBAFqtlipChBBCSA7JdiLUuXNndOjQAZcuXcKzZ88A6DtL16tXz1DhIOlLWxFSq9WUCBFCCCE5xCajxoRCIRo2bGiy/eLFi6hfv74tniJf4ytCfFMiNY0RQgghOcOu8wj16tXLnofPN/iKkLu7O0aPHo1KlSo5OCJCCCGkYMh2Rah3795mtzPGEBcXl93DFwh8RahYsWL45ZdfHBwNIYQQUnBkOxE6fvw4fv/9d7i4uBhtZ4zh9OnT2T18gUDLaxBCCCGOke1EqEWLFnB1dUWzZs1M7qtWrVp2D18gpO0j9OrVKxQqVIg6mhNCCCE5INuJ0O7duy3ed+zYsewevkDgK0IXLlxA4cKFER4ejvLlyzs4KkIIIST/o0VXcwG+IsSvN0bD5wkhhJCckeWK0Pjx481u5zgOMpkMZcqUQbdu3eDl5ZXl4AoKviLEGANAw+cJIYSQnJLlRCgsLAzXrl2DVqs1NOM8ePAAQqEQFSpUwIoVKzBhwgScPXvWLsPBZ82ahYMHD+L69euQSCR4//690f03btzA3LlzcfbsWbx58walSpXCyJEjMXbsWJvHkl201hghhBDiGFluGuvWrRvatGmDZ8+e4erVq7h69SpiY2PRtm1b9OvXD0+fPkWzZs0wbtw4W8ZroFKp0KtXL4waNcrs/VevXoWvry82b96MO3fu4Ntvv8WUKVNy5fD0tKvPC4VCSCQSB0ZDCCGEFBwc49tjMql48eI4duyYSbXnzp07aNeuHZ4+fYpr166hXbt2ePPmjU2CNWfDhg348ssvTSpC5owePRr37t3DiRMnrD5+QkIC3N3dER8fDzc3t2xEatn48eOxePFiAICrqysSEhLs8jyEEEJIQWHt93eWK0Lx8fF49eqVyfbXr18bvsg9PDygUqmy+hQ2Fx8fn2GfJaVSiYSEBKObvaWtCFGzGCGEEJJzstU0NmTIEOzZswexsbGIjY3Fnj17EBISgu7duwMALl26hHLlytkq1mw5d+4ctm3bhuHDh6e735w5c+Du7m64+fv72z02vo9QgwYNMGDAALs/HyGEEEL0spwIrVq1Cq1bt0bfvn1RsmRJlCxZEn379kXr1q2xcuVKAECFChWwdu1aq485efJkcByX7i08PDzTsd6+fRvdunXD9OnT0a5du3T3nTJlCuLj4w23mJiYTD9fZvEVob59+2L+/Pl2fz5CCCGE6GV51JiLiwvWrFmDxYsX4/HjxwCAgIAAo6U2atSokaljTpgwAcHBwenuExAQkKlj3r17F61bt8bw4cPx3XffZbi/VCqFVCrN1HNkFy2xQQghhDhGtmeWdnFxsdlSGj4+PvDx8bHJsQB9x+1WrVph0KBBmDVrls2Oa2t8RUij0SA1NRUymczBERFCCCEFQ7YSoffv32PdunW4d+8eAKBSpUoICQmBu7u7TYJLT3R0NOLi4hAdHQ2tVovr168DAMqUKQMXFxfcvn0brVq1QlBQEMaPH48XL14A0A9Pt2WyZQt8RWj06NE4cuQI9u3b5+CICCGEkIIhy32Erly5gsDAQCxevBhxcXGIi4vD4sWLERgYiGvXrtkyRrOmTZuGmjVrYvr06UhKSkLNmjVRs2ZNXLlyBQCwc+dOvH79Gps3b0bRokUNt7p169o9tsxKO6EijRojhBBCck6W5xFq2rQpypQpgzVr1kAk0heWNBoNhg4disePH+P06dM2DdRRcmIeoWrVquHWrVsAgJCQkEx1MCeEEEKIKWu/v7PcNHblyhWjJAgARCIRvvrqK9SpUyerhy2QqCJECCGEOEaWm8bc3NwQHR1tsj0mJgaurq7ZCqqgoQkVCSGEEMfIciLUp08fhISEYNu2bYiJiUFMTAy2bt2KoUOHol+/fraMMd+jihAhhBDiGFluGluwYAE4jsPAgQOh0WjAGINEIsGoUaMwd+5cW8aY76WtCDk5OTkwEkIIIaRgyXIiJJFIsHTpUsyZMwePHj0CAAQGBtIXeSYxxgwVoU8++QRVqlRxcESEEEJIwZGpRGj8+PFW77to0aJMB1MQaTQa6HQ6AMC6devg4eHh2IAIIYSQAiRTiVBYWJhV+3Ecl6VgCqK0/YNoRmlCCCEkZ2UqEQoNDbVXHAVW2kQI0DeVUSJJCCGE5IwsjxojtsF3lOY4DnK5HP/++6+DIyKEEEIKDkqEHOzjihANnyeEEEJyDiVCDpZ26DxAiRAhhBCSkygRcjC+IsQv+UaJECGEEJJzKBFyMKoIEUIIIY5DiZCDfdxHiCakJIQQQnIOJUIOlrYiJBKJIJFIHBgNIYQQUrBkeYkNYht8Rahw4cJo2bKlg6MhhBBCChaqCDkYnwjVq1cPf/75p4OjIYQQQgoWqgg5GN80Rstr5G0qlQoajcbRYRBCSL5lr+4jlAg5GF8RkkqlDo6EZEVcXBxePI1Fikrt6FAIISTfk0vEKFLcD15eXjY7JiVCDsZXhDZv3oyoqCicPn3awRERa8XFxeHJk8dwe3kJRWMOQaJ4CQ46R4dFCCH5DoMAKqfCeOPfAU+U9QDAZskQJUIOlnb4PFWF8pYXT2Ph9vISylz6BhyYo8MhhJB8zTn+Pjyen8HDerMQkZoEf7/iKFaqbLaPS52lHSzt8HmaTDHvUKlUSFGp4R1ziJIgQgjJIRwYvGOOQCR3xeHt6/As6mG2j0mJkIOlrQjRZIp5B98xWqJ46eBICCGkYJEoXgAAkhMTcOviqWwfjxIhB6OKUN5GfYIIISRn8Z+7Lm6eeBJ+A2qVKlvHo0TIwdJWhCgRIoQQQqwjkkihVimhUqZkvHM6KBFyMEqECCGEkCzgOIAxMJa9fpqUCDkY3zRWvXp1VKpUycHREEIIIQULDZ93ML4i9N1336Fnz54OjoaQzONmJgAASrhzuP8/F8hEnMk+pZYkIiqeQT3VFSKB6f15yRuFDlV/TcaLJIbG/kKcHWK+knvwgRpLL6pw97UOb1MYirpwqF1MiPENJGjoTx+9hOQWVBFyMFpig+QX0fEMSy5kr9NiXjDiQCqSVOmX4r8+lorOf6bg2nMd2pcRYWx9CWoVFWJfuAaNf1Ng8838/zoRkldQIuRgfEVILpc7OBJCss5TBnjJOcw9q8QbRf4dSbfphgq772kwr43lHy4vknRYcF6Fws4c7o52xtqucsxtI8PO3k44MsAJDMC0UGXOBU0ISRclQg7GV4Q6dOiAgwcPOjgaQrLGScxhajMJ4pXAzJPWfcl336oANzMBP1803X/qiVRwMxMQsi97o0FsKTpehzGHUhFSU4wOZS03bUW910HHgPp+Qvg6G3/EtiwtgqsEeK2gSTgJyS0oEXIwviKkVqvtsqouITlldF0JAj05rLqqRsRbbYb7/9ZNjhLuHL46pkTY8w/7//NYg9lnVajkI8CyjrmjyZgxhuC9KXCXcVgUlH5MZQsJIRECl55qTapjp6M0SFQBbQKojxAhuQX9NToYTaiYfyWn049EKIBRp+L09hVwgFyctX0VagZLI0udJbbttCwWcpjbRoZeO1Lw9XEldvdJf6Z0LzmHPz+Vo/kGBfrsTMG1Ec5IVjEM2JMCqRDY3lMOJ7F1Me4NV+P6i4yTL56HjMOXDaxf22/JBRVORmpx9HMnuEk5xKVYfg+85Bx+aiPF+CNKVFqejO4VRCgk5/DonQ7772vQNkCIVZ1zR4JHCKFEyOFoHqH8y2VOosX7OpYV4WD/D4mC74JEKNTm921eUoiTwR+ujVJLk/DGQtNKnWICXB7mYvh3peVJiIo3vy+b7pZe+FnSs5IYDf1U2BOuwdloDZqUSP8jppG/CD+0lGLKP0qMOJCC18kML5IY1nSRobKv0Orn3RuuwcYbFl5AM0q6W58I3X2txTcnlBhZR2x1JefLBlKU8hBgyL4UrLn2Ia4yXgIE15CYNJkRQhyHEiEHYozRWmMk31nYTopGvykw8agSF4Zm/BHzdWMJQiM12HJLv35bvyoiDK2VuWbiDd3l2NDd9gMO1FqGz/ekoKgLh3ltra/izPtXiW/+UWJMfQn+V0+CIi4cwt/oMOWfVHy2OwXXX2gzdTxCiP1QIuRAarXaaEZMqgjlL0lTXC3eJ/yoIPBqouV9P552J3Ksi/kdzex7d7SLxaYxe2noL0LPSiLsvKvBtttq9KkiTnd/juPwSQUxjj7SN21lpsnK3uacVSHsuQ6hg5zgYmVT4slIDb4+rkSPCiKj/kS1igqxp48Tyi1LwsLzKoysI0GAJ1WGCHE0SoQcKG01CKBEKL/JTB8ce+1rbR8bW5vTWoZ94UmY8k8qelRM/2Mm4q0WE4+lwlMGxCuBoftTcGmYs9mJGS2xVx+ha8+1YABabFSYvf/fGC24mQlwlwLvJ+ubGg880Fe2WpYyPW8nMYd6xYXYE65B2HMtJUKE5AKUCDlQ2o7SzZo1o6Yxkm+U8RLg/+pKsPSiCssuWp48UKlh6LMzBckq4MgAJ5yK0mDWGRW+PJyKlZ2tb+qyVx+htgEieDuZJmRJKoZtdzQo7MyhczmRUcKp1OhLcK8tzKfED52XWN8FihBiR5QIOVDayRRPnTrl4GgIsa1pzaXYeEOFWWeUEHDmqzsTj6Yi7IUOXzeWoG2gCK1KC3EyUotVV9VoXVqEXpXTb1bj2auP0Oh65vsqRb7XYdudJJTxEmBtV+PnbVpShF8uq7H6qhojaktQ3O1D1edQhBr/RmshEwGN/CkTIiQ3oLqsA9HyGiQ/85Jz+KaJFO9SgbdmhpvvuafGL5fVqF9ciB9b6aszQoF+SL2XnMPQv1Lw+F3em6W6ZyUR2gQI8TKZoeLyJAzam4Kvj6Wi658KdNqSAgZgbmspCjnRxy8huQH9JToQLa9B8rsx9SUo5WFaDYqO1yFkfwrcpcDWnnKjhVj93QX4rasMCUqg704FVNq8NQuzgOPwd38nLA6SopKPEHvuqbHwvAoXYrXoWFaEIwOcMDYXdQgnpKCjpjEH4itCz58/R8OGDXH+/HkHR0RI5qU3H5FUxOHJWNMRcSXcBYj72vLjulUQg023rlnMEUp5CNI9b7FQ3wcpN42AI4SYRxUhB+IrQowxo47ThBBCCMkZlAg5EM0qTQghhDhWnk2EZs2ahUaNGsHJyQkeHh7p7vv27Vv4+fmB4zi8f/8+R+KzRtoqEA2dJ4QQQnJenk2EVCoVevXqhVGjRmW4b0hICKpVq5YDUWUOVYQIIYQQx8qzidDMmTMxbtw4VK1aNd39fv31V7x//x4TJ07MocisRyvPE0IIIY6Vr0eN3b17F99//z0uXryIx48fOzocE7TgKiGEEOJY+TYRUiqV6NevH+bPn48SJUpYnQgplUoolUrDvxMSEuwVoqEi5Ovri/Lly9vteQghhBBiXq5qGps8eTI4jkv3Fh4ebtWxpkyZgooVK2LAgAGZimHOnDlwd3c33Pz9/bNyKlbhK0I9evTAhAkT7PY8hBBCCDEvV1WEJkyYgODg4HT3CQgIsOpYJ06cwK1bt7Bz504A+rl6AMDb2xvffvstZs6cafZxU6ZMwfjx4w3/TkhIsFsyREtsEEIIIY6VqxIhHx8f+Pj42ORYu3btMuqDc/nyZQwZMgRnzpxBYGCgxcdJpVJIpTkzGywtsUEIIYQ4Vq5KhDIjOjoacXFxiI6OhlarxfXr1wEAZcqUgYuLi0my8+bNGwBAxYoVM5x3KKfwidAvv/yCGjVqoE+fPg6OiBBCCClY8mwiNG3aNGzcuNHw75o1awIAQkND0aJFCwdFlTl801hSUhJ0ury3yjYhhBCS1+WqztKZsWHDBjDGTG6WkqAWLVqAMZZrqkEADZ8nhBBCHC3PJkL5AU2oSAjJq05GasDNTMCMk7ZZMHrGyVRwMxNwMlJjk+NlxVuFDl4/JeD/DqZkvHM+c/WZFtzMBKy9prL6MUkqhh9PKzHjZCruvtbaMTr7yrNNY/kBLbFB8gNupn6urRLuHO7/zwUyEWeyT6kliYiKZ1BPdYVIYHp/XvAiSYefzqrw90MNouN1kImA0h4CtAsUYW6b9Ed+no7SoOVGBXQM+LapBD+2yr8jRVtsSMapKC3YdDdHh5Jp008qkaIBvmtmfsBMfCrDwvNK7A3X4PE7HQQcUMJdgMb+QvzSUQax0PK1/eCtFjVXJUOhBj6rKsbmTywPkrn8VIt6a5OxoK0UExqZxsIYQ7vNChx/rE8+LP1d3Xqpxdx/lbgYq8XTRAYvOYdyhQQYWVuCXpVFEHAfHlO7mBDdK4gwNVSJvlXEcJGk/3eaombovEWBU1H6GFZdVeN0sBPKFhKm+7i0lBqG2quTcee1DsVdOcSOd7X6sbZEFSEHokVXSX4SHc+w5IL1vybzkn+jNai4PAkrrqhQyUeAsfUlGFhNgqKuAmy9rU73sYlKhkF7U+AkzqFgSZZEx+uw6qoan1cTo5ir6Vdj+Bstqv6ahFlnVCjqymF0XQmG1pKgXCEBtt9VQ5lOQUSjY/h8Twqs/Q2wJ1x/TX1S0fxF88slNUKfaCFLp5Tx1301aq1Oxs67GtQsKsTY+hJ0KCPCzZc69N2VghF/mVbypjSR4kUSw88X0/87VmoYum/TJ0F9q4iwoK0UL5MYWm1S4Mk76/u7fvOPElHxju8fSxUhB6KKEMkvPGUAx3GYe1aJobXE8HbKP7+xXiTp0G1rCtylHC4OdUK5j37xqrUs3cePPZyK+FSGKU2k+PaEMt19ieOsuqKCRgcE1zBNPhRqhq5/piBRxfDvECc08DP+6tToGNIpBmH2GRWuv9Bhflspxh7O+BrYE65B9cIClPY0/Tu6/0aLr4+nYmIjCbbeViMq3vz1N/kfJTQ64OQgJzQv9SHeH1vpUH1lMtaGqTG1uRQl3D88R73iQlTwFmDVVRUmN5EYVYx4ai1Drx0pOPpIiy/qSbC0vRQcx6GoqwDBe1PQelMyTg92hp9b+p8BJyM1WHxBhRWdZBh10DbNq1mVfz6t8iC+IhQQEAA3t7xXRiaE5yTmMLWZBPFKYOZJ677su29VgJuZgJ8vmu4/9YS+v0jIPsf31Zh9RoW3KQwrO8tNkiAA6TaH7AtXY/11NX7uIEMxV9s2CQbvTQE3MwFP3unwyyUVKi1PguzHBJRakojZZ5SGSWR33FGj3pokOM9OgO/8RPzv7xSkqI2/PCPf68DNTEDwXvOvd4sNyYYmUEv4Y/BNJdzMBMOtxYbkTJ3bxusq1FyVBPksfcxD9qXgRZJx5aDhumQIZiYg8r35isLCc0pwMxOw4FzG1yNjDOuvq+HvxqGRv2l9YOUVFSLidJjTWmaSBAGASKBf+cCcK8+0+OG0ElObSVGtcMbNRuFvtAh/o0OPCqYJGV9ZCvAUYGaL9Oe7e/xOBzcpjJIgACjiIkD94vo4XiebJlF9K4sRHc9w7JFpiUurY/hsdwr+eqDBjy2l+LmDzHDe/auK8Vc/J7xKZmi9SYGXSZYrPQlKpk+aAoQYWUeS7nnkBEqEHIivCP3+++/w9fV1cDSEZM/ouhIEenJYdVWNiLcZd5z8rZscJdw5fHVMibDnH/b/57EGs8/qm6CWdXR8X5o/b6vhKQOCAoW4+1qLZRdV+OmsEjvvqpGkslwNepWsw7C/UtG9gggDqtnvw37isVRMP6lEveL6LxUBB3x7QomZp5T4+aISg/amoIyXvl9IERcOyy+rMf6I7X+Be8g4TG8uQUl3/Rfj9OYSw81clcWSxRdUGHkwFdULC/FlfQnKewuw/roajdYl43Xyhy/XUXXEYADWXDXfjLP6mhpSofkKz8fuvNbheRJD4xLmE5Utt9TgAPStIkbkex1+vazCnDNK/HFTjbcKy1/4KWp94lKjiACTm1h3Dey+p+8s3qOiacL142klwl7osKG7HFIzffHSquwjQIISOBtt3Pn8VbIOl55qUdSFQyUf0xSAfw2OPTZ+nI4xDN6Xit33NFjbRYZvzfSjCiojQuggZ8Sl6JOhNxZemzGHUvEulWFd19wxmTA1jTkQnwjREhv5C2MMivS7jeQaTmJY/CWbWWIhh7ltZOi1IwVfH1did5/0+715yTn8+akczTco0GdnCq6NcEayimHAnhRIhcD2nnI4ia2LbW+4GtdfWD9qxUPG4csGGc8g/+SdDm8UDHWLCTDuiBJLP+o7UUjOYVMPGTqWNf2yHfZXKnQMWNnJvn/fV59pcXOkM4r/1xQxo4UUZX5OwvxzKjiJOVwd7oyKPvovN6WGoeaqZPx2XY2ZLaXwdbbdb2EPGYcZLWQ4GalFVLwWM1pk7bwPRWhwcagzahb9kJSMO5yKJRdVmHxciXXd9F+evSuLMf6I0nAuaTsLn4zU4MFbHfpXFVnVTHs2Wn/t1ClqmgiptQw3Xurg48xhzVUVvjmhb3LiOYuBnzvIMKSmaaIz+bgST97pcG2Es9WDBPaEqxHoyZlUjy4/1WLWGRUmN5agTrGMK0uLg2To/KcCbTYp0K2CCAEeArxRMOy9r4GHjMOWT+SQm/n7qvvfsU9HGSdCAo7Dph5ybOqRfvJSt7gQrydZ7vS8554aG2+osbaLzKhZzpEoEXIgvmmMltjIXxRqwGVOoqPDsErSFFc427BY0bOSGA39VNgTrsHZaA2alEj/I6aRvwg/tJRiyj9KjDiQgtfJDC+SGNZ0kaGyr/WjT/aGa7DxhvXZZ0l36xKhV/9VIK491+H2KxV+6SBD78oiaHTA5ptqfHNCiU+3p+DacIEh2QCA38JU2H9fg2095SjsYt8P+6nNpIYkCNAnJF3Li7D+uhoTGoqN4pKKOPSpLMaMU0rce62zaSJkK59XExslQYA+uVt/XYUtt9VY0UkGqYiDTMRhcA0xFpxXYV+4Bp9W+pCMrvqvSjSitnUXd/R/HXaLmukkHZfCoNEBbxUMU/5RYlpzKYbUFEMu0l93Xx5JxdD9qSjlIUCr0h+u938ea7Dskgpz20hRyce6azkmXocrz3SY2NA4br6yVNlHgGnNrVsCqmlJEc6HOKP3jhRsv/MhqXGVAINrSFC1sPn33l3GQSbSD36wtZdJOgw/kIoOZUQIqeX4JjEeJUIOpFAoAAB9+vTBzZs3HRwNIbaxsJ0UjX5TYOJRJS4Mzfgj5uvGEoRGarDllv7Dul8VEYZm8kNyQ3c5NnS3/Q8K3X/fBVoGTGsuxeh6H+Ka1Fg/wmbRBRWWXFBhVRf980e+1+HLw6noVUmE3pXtP1TMXHWA749U20yFo7ib/r7YBNt/0dlC81KmMbvLONQoIsSpKC3uvdGhRhH9PqPqSrDwvAqrrqoMidAbhQ577mlQ0VuAZiWt+4p7q9C/Fp5milhpr4ERtcVGiUhILQkUamDM4VT89K/SkAi9T2UI3peC+n5CTGho/bW8N9x8s9hXx1Lx+J0Ol4c5p9snLa1jjzTouysFdYoJsKmHMyp4C/AiieGXSyp8e0KJgxEanAp2Mlup8pJzeJlk++tj2F+p0OgY1nbNXa0glAg5CGMMSqW+E9/r168dHA2xJSexvtKSF9hjSHdDfxF6VhJh510Ntt1Wo0+V9J+E4zh8UkGMo/91zrSmUpNTPGQfviR6VDD9uOxRUYRFF1S49OxDs9yQfSmQizmssHOTGM9dZvpFJhJkfJ9alzsTocIWqlRFXPTnEp/6Ie4ATwGCyghx5KEWj+J0CPQSYON1/VD2EbWtv7j5JqJUM3M5uhtdA6bH7FFRhDGHgUtPP1wD44+k4q2C4fjnMggzMW/W7nA1irhwaOD3IRk8FanB8stqzGghRfUi1lWW4lIY+uxUwEnMYU8fJ0MTc4Anh0VBMjx5r8PecA0231QjuIZpopaiZpDb+LNh0w0V/nqgwcbuMrPTEzgSJUIOolKpDKM6aOh8/sJxnE2bm/KiOa1l2BeehCn/pJrt9JlWxFstJh5LhacMiFcCQ/en4NIwZ7MTM1pirz5CgV4CiASARmecFPE8/9uWkqZV7tpzLeKVgM/8JLPHnHVGhVlnVOhWXoS9fXPH/GH8d7XGQr/f96k5lzS9TDYfxIv/KhQfJ3ej6khw+GEK1lxTYW4bGVZfU0MmAgZWt/6P0NdZf8y3Kabn6STm4O/GISaBwcNMbmvpGkjRABWWmx8t98ctNf64pUb1wgJcH+kCQF/JOhOlxbBaYqNh62EvtGDQT/Y43cKITPEP+qb4sBHOqFFEiHMxGrxLBVqWFprtZ9eylBB7wzW4+kyL4BrG9+kYw/tUoLSnbUc5Xnuuf18H7U3FoL2mnfWfJjLDyMR3X7ua/XuzF0qEHCTtZIouLi4OjIQQ2yvjJcD/1ZVg6UUVlqUzOZtSw9BnZwqSVcCRAU44FaXBrDMqfHk4FSs7W9/UZa8+QhIhh6YlhAiN1OL2K51Jf5/br/Qf7mm/NAZWl0ChNv1CjYjT4XSUFjWKCFC7qBA1rfx1nxP4L/OYBNMkJEHJ8OCt9ZPeCf97ibQ6lqlqCO9UpBYDqxtvi09luP5CP4FgRW/j96BzORFKuHNYf12NVqVFePBWh4HVxfCUW//c1f7rLxP+xvx5tgnQ97m6/UqH+n7G991+pU/A087580lFMeoUMz3W8ySGvyM0CPTk0KKUyKiz8F/3NdAyoMdHkyhW8RUipKb58sy2O2okqYAhNcTgOH3nfQBQ/lfZMjc8Pu12iZlmtvtvdGCAofnRVhr6CZGkMn8e68LUcBID/f6rHktz+E+DEiEHockUSX43rbkUG2+oMOuM0uzEbAAw8Wgqwl7o8HVjCdoGitCqtBAnI7VYdVWN1qVF6GVlHxt79RECgC/qSRAamYJpJ5Vo4CeE839LD7xPZfjhtP4Xer80zX8/dzDfJLbhugqno7ToVFZkdokNfhmSJ2NdUMojZ5sOXKUcKngL8G+0Fndfaw2de7U6hvFHUpGSieW/+C/j6HiWparC7zfV+F89iVGH6RknlYhXAoNriE2GjQs4DsNrSfBdqBJD/pt3amQmmsUAoGkJEYQccCHWfFVxdF0JNt5QY+6/KnQtL4LPf813qRpmmCQz7TVgqUPzyUgN/o7QoIGfCGs/Gjq+J1wDD5m+WpNWmwAR2gSY/6o+/liDJBXDqi4yo74+Df2FEAmAf2O0OPpIg3aBHx4f898M2gDQOsA04+Bfg5albJse9KkitthMvi5MDU8ZZ/Ka5BRKhBzkyZMnhv+nihDJj7zkHL5pIsVXx5UATH+Z7rmnxi+X1ahfXIgfW+m/OIQC/ZD6GquSMfSvFNQuJkSAmdl1c1KPimIMrqHB+utqVP01CR3KiKBlwIEHGjxNZPi0oggDqmW/QwXfZUfkoNOd1EiCkP2paPxbMnpVEkMm4hAaqYFaC1QvLMCNl9ZVhVqXFmHHXQ0+2a5AxzIiyMUcSrpz+NzKpqoOZUVo/FsyelcWo6gLh7MxWpyN1qKUB4e5bcwnGENrifH9aSWeJjJU9RWgoZlJEdPjLuPQOkCfhL9LYSbVpNrFhJjeXIrpJ5Wo8msyupYTQSYCjjzSIiJOh0b+QnzVOOvt4UkqhmOPNfi0otjqztDpKeYqwNRm+ng7/KFA53IiVCgkwItkht339FWkHhVEZqd9OPpYAyEHdCtfcNKD3NVjqYBYt24dmjRpYvg3dZYm+dWY+hKU8jD9YI+O1yFkfwrcpcDWnnKjX7P+7gL81lWGBCXQd6cCqgyWsMgJ67rKsKaLDIWcOGy4ocamG/pOrb90kGF7L7nFipe13qUwPE1kaOwvzHBpAnsZUlOCtV30HVk33lBj+x01GvkJ8e8Qp0z11xhaS4wpTSSIT2WYd06FqaFKrAuzvtlyXAMJVnSS4foLLZZcVCH8jQ7BNcQ4N8TZ4nD/wi4CdCyr/+K2dsj8x/6vjgQqLSyuHTetuRS7estRvpAA2+6oseaaGmIh8GNLKU4MdMpUn7aPHYrQIFUDfJJBf7rMmNZcir195GgXKMS5GC0Wnldhzz01qvoKsaKjDDt6mVZf4lMZ9oZr0LmcCP65ZI6fnMAxvscuMSshIQHu7u6Ij4+3yTIYsbGxKFmyJHS6D7+uOI5DdHQ0/Pz80nkkyU0UCgXu3buHiqdHwCk+wtHhkDxu/301um1NwcH+crO/0kn6dIyhzM9JeJnM8HyCK9ykmU9KtDqGqr8mQyLUdzq21USj1ui/S4G94Rq8+crV6klE7WHZRRXGHE7FmcFOGc4B5kgK97K412wV7v77N96/jMHgr+fBxc3DZD9rv78LTsqXS0RERBglQYB+KP3Dhw8dFBEhxNFOReo7UVMSlDU772rw5D3DwGriLCVBgL5ZdkE7KW681BmWucgJKi3DwQh9Px5HJkEpaoY5Z5X4tKIoVydB9lCwzjYXKFu2LAQCgVEyJBQKUaZMGQdGRQhxpIVBuWuCubxi7lkl4lIYVl9VwVkMTGmavTmoOpYVY2l7ndn5hOxFIuQQP9nxi25HvtdheG2x2XmF8jtKhHKYn58fVq9ejeHDhxuSoVWrVlGzGCGEZNKUf5QQC4BKPgLMb2ubtavG1M89E3rmpIo+QsxokXumdMhJlAg5QEhICDQaDUaOHAmZTJajbdGEEJJfsOmOr6SQvI/6CDkIP3dQamoqkpLMz0BLCCGEEPuiRMhBaEJFQgghxPEoEXKQtEtsUCJECCGEOAYlQg6StiLk5JQ7Fl4khBBCChpKhByEKkKEEEKI41Ei5CDUR4gQQghxPEqEHIRPhNzc3ODu7u7gaAghhJCCiRIhB+GbxiZMmICKFSs6OBpCCCGkYKJEyEH4ipBcbroCMCGEEEJyBiVCDsJXhGQyWmOIEJL3nIzUgJuZgBknUzPe2QozTqaCm5mAk5E5uNDXR94qdPD6KQH/dzAl453zmavPtOBmJmDtNZXVj0lSMfx4WokZJ1Nx97XWjtHZFy2x4SAJCQkAgIULF2L06NEQCCgnJXkTN1N/LZdw53D/fy6QiUyXjCm1JBFR8Qzqqa4QCfLekjLXnmsx56wSp6O0eJfCUMSFQ+dyIkxvLkVhF8t/u/881uCXyyqcj9HiXSpDITmHqoUFGFtfkm9Xmm+xIRmnorR5cvmL6SeVSNEA3zUzXW/sYZwOs84oceyRBq+SGQo5cWgTIMKM5lIEelm+Bq4912LBOf2181rB4CHjUMFbgJCaYgysbn6B08tPtai3NhkL2koxoZFpLIwxtNuswPHH+uTD0t/VrZdazP1XiYuxWjxNZPCScyhXSICRtSXoVVkEQZrlnWoXE6J7BRGmhirRt4oYLpL0/05T1AydtyhwKkofw6qrapwOdkLZQtavV6bUMNRenYw7r3Uo7sohdryr1Y+1Jfr2dZDk5GQAwNOnTykJIvlCdDzDkgvW/5rMKw48UKPB2mTsuqtBAz8hxtaXoHoRIVZeUaPOmmREx+vMPu6rY6lo87sCV55p0bW8CBMaStCprAivkxlORubdX8/5VXS8DquuqvF5NTGKuRp/Jl95pkWtVUnYcF2NCt76RLZpCSG23laj1uokhD03/37+ckmFumuScfSRFq0D9NdAjwoiaHXA3xGWK197wtUAgE8qmk+Wf7mkRugTLWTplDL+uq9GrdXJ2HlXg5pF9ddthzIi3HypQ99dKRjxl2klb0oTKV4kMfx8Mf2/Y6WGofs2fRLUt4oIC9pK8TKJodUmBZ68M//3YM43/ygRZeHvJydRRchB+ESImsZIfuApAziOw9yzSgytJYa3U/5I7lM1DEP3p0KtA3b1lht9Mf15S43+u1Pwv79Tsb+f8aSoa66qMP+cCoOqi7G6iwwSofGva7WW5Uj8xHqrrqig0QHBNUyTj5D9KUhUAYvaSTGu4YcKzdloDVpsUGDwvhSEjXA2WkD76CMNxhxKRdtAIXb2coKr1PprYE+4BtULC1Da0/Tv6P4bLb4+noqJjSTYeluNqHjzx5n8jxIaHXBykBOal/rwVf9jKx2qr0zG2jA1pjaXooT7h+eoV1yICt4CrLqqwuQmEqOKUdq4e+1IwdFHWnxRT4Kl7aXgOA5FXQUI3puC1puScXqwM/zc0v8MOBmpweILKqzoJMOog7ZpXs2q/PFplQfxnaWlUtOyJyF5jZOYw9RmEsQrgZknlVY9pvtWBbiZCfj5oun+U0/o+4uE7HNsX41zMVq8TGaoU0xg8uu8X1UxqhcW4MADDaLef/hVq9QwfHtCiRLunNkkCADEZrZlVvDeFHAzE/DknQ6/XFKh0vIkyH5MQKkliZh9RgnG9F+QO+6oUW9NEpxnJ8B3fiL+93cKUtTGX56R73XgZiYgeK/517vFhmRDE6gl/DH4phJuZoLh1mJDcqbObeN1FWquSoJ8lj7mIftS8CLJuHLQcF0yBDMTEPnefEVh4TkluJkJWHAu4+uRMYb119Xwd+PQyN+4PvD4nQ43X+rg68xhbAPjpqwmJUToXE6EGy91OBNtXBWadCwVcjGw5RO5SRIEWL4Gwt9oEf5Ghx4VTBMyjY7h8z0pCPAUYGaL9L87Hr/TwU0KoyQIAIq4CFC/uL756nWyaRLVt7IY0fEMxx6ZVrm0OobPdqfgrwca/NhSip87yAzJX/+qYvzVzwmvkhlab1LgZZLlSk+CkumTpgAhRtYx3zyYkygRchAaNUYyIzZBh9AnGsQmOL6MbMnouhIEenJYdVWNiLcZN/381k2OEu4cvjqmNGpa+OexBrPPqlDJR4BlHR1bMX2RpP+iCDDzy5zfzgCcePKhmePYYw1eKxg+qSCGgAMOPlDjp7NKLL2gxPkY23cEnngsFdNPKlGvuP5LRcAB355QYuYpJX6+qMSgvSko46XvF1LEhcPyy2qMP2L7X+AeMg7Tm0tQ0l3/xTi9ucRwM1dlsWTxBRVGHkxF9cJCfFlfgvLeAqy/rkajdcl4nfzh+h9VRwwGffXNnNXX1JAKzVd4PnbntQ7PkxgalzDt38InYKU8OLMVEv7a+Ofxh/f29istbr7UoV2gCF5yDqFPNFhwTomF55T457EGOma5GrT7nv44PSqaNtj8eFqJsBc6bOguh9RMX7y0KvsIkKDUV63SepWsw6WnWhR14VDJx/S65l+DY4+NH6djDIP3pWL3PQ3WdpHhWzP9qILKiBA6yBlxKfpk6I3C/OfVmEOpeJfKsK5r7vj+o6YxB+FHjdE6YyQj666pMPxAKnQMEHDA6s4yhNRy/K+oj4mFHOa2kaHXjhR8fVyJ3X3Sv7a95Bz+/FSO5hsU6LMzBddGOCNZxTBgTwqkQmB7TzmcxNZVTvaGq3H9hfX9bjxkHL5skHE11ttJ//yW+j08/m/7/bcf7r/8VB+HTATUXJWM26+MH9uspBA7e8nh42yb36FXn2lxc6Qziv/XFDGjhRRlfk7C/HMqOIk5XB3ujIo++i83pYah5qpk/HZdjZktpfC1UQyA/jWd0UKGk5FaRMVrMaNF1pLYQxEaXBzqjJpFPyQl4w6nYslFFSYfV2JdN/2XZ+/KYow/ojScS9rOwicjNXjwVof+VUVWNdOe/a+aU6eoaSLEXwNR7xkYY0bNX0D614CvE4cWGxU4HWV8bVb1FWB3HyeUMdPJek+4GoGeHKoVNo7l8lMtZp1RYXJjCeoUy7hD8uIgGTr/qUCbTQp0qyBCgIcAbxQMe+9r4CHjsOUTOeRm/r7q/nfs01HGiZCA47CphxybeqSfvNQtLsTrSZY7Pe+5p8bGG2qs7SIzapZzJEqEHGTGjBkYOXIkvL29HR0KycViE3SGJAgAdAwYcSAVQWVEGbbBO0LPSmI09FNhT7gGZ6M1aFIi/Y+YRv4i/NBSiin/KDHiQApeJzO8SGJY00WGyr7Wjz7ZG67Bxhtqq/cv6W5dItTYXwgPGXD5mQ77wtXolqa5YvsdNW681H/5vUv58Av/1X/NDfPP6ataZwY7oUYRIZ6802HisVQcfaRFrx0pOBlsm6V1pjaTGpIgQJ+QdC0vwvrrakxoKDYkQQAgFXHoU1mMGaeUuPdaZ9NEyFY+ryY2SoIAfXK3/roKW26rsaKTDFIRB5mIw+AaYiw4r8K+cA0+rfThvVn1X5VoRG3rfjDwHd6Lupq+HuUKCVHWS4CIOB1+vqjC2DTXzbkYDQ480CcM71JNr4F1YWoUd+NwsL8cTUqI8DJJh+9Pq7D5phqdtihwa5SzUdNpTLwOV57pMLGhcdwpan2TWGUfAaY1t647RdOSIpwPcUbvHSnYfudDUuMqAQbXkKBqYfPvvbuMg0ykH/xgay+T9J9nHcqIctWPOUqEHMTX1xf+/v7w8/NzdCgkF4t4qzMkQTwt0w/lzY2JEAAsbCdFo98UmHhUiQtDM/6I+bqxBKGRGmy5pf+w7ldFhKGZ/JDc0F2ODd1tX2Z3lnBY2l6G4L2p+GR7CrqWV6OslwDhb3Q48ECDGkUEuP5Ch7Qjl/n3SyQA9vdzQikP/ftUtbAQe/o4ofwvSTgVpcX5GA0a+mf/I9hcdaCYqz6g2mYqHMXd9PfFJuTODtvNS5nG7C7jUKOIEKeitLj3RocaRfT7jKorwcLzKqy6qjIkQm8UOuy5p0FFbwGalbTu9X2r0L8WnhaKWCs7y9DhDwW+PKLEgQgNahQWIiZBh933NKha2PI1oGXA1k/lhvfZTSrEpu4yhL/R4sozHXbd1aBf1Q8J3N5w881iXx1LxeN3Olwe5mx1/7JjjzTouysFdYoJsKmHMyp4C/AiieGXSyp8e0KJgxEanAp2Mjvs3kvO4WWS7a+PYX+lQqNjWNs1dw0Syp2fpAVAjx49EB0djW3btjk6FJKLlS0kwMefU0IOZkvquUVDfxF6VhLh4lMttt3OuErDcRw+SVNpsaZSk5MGVpfgn4FOaBsgxMlIDZZeVCE6XocN3WX47L8vMV/nD2+Sh0z//zWLCg1JEM9JzCEoUP8ld+mpbYbQu8tMv8hEgozvU3+cYecShS1UqYq46M8lPk3lJcBTgKAyQhx/rMWjOH1VZ+N1NZRaYERt6/sl8U1EqRa6cLUqLcKFEGd8UlGE6y90WHpRhesvdPipjRRTmuiv17TVNf4aKOLCmSS7HMehW3l9bB9fA7vD1SjiwqGB34dk8FSkBssvq/FdMymqF7GuShqXwtBnpwJyEbCnjxNqFRXCScwhwFOARUEydK8gwrkYLTbfNP/3maJmkNt4mqtNN1T464EGS9vLTKYncDSqCBGSi/m5CbC6swwjDqRCy/RJ0KrOslxbDeLNaS3DvvAkTPkn1Wynz7Qi3mox8VgqPGVAvBIYuj8Fl4Y5m52Y0RJ79RHitSwtQsvSpucxcI9+0EPd4h++oMp7C/57DvPH8vzvSzLFcRMom+CTbY2FvvjvU3MuaXqZbD4IvuP6x8ndqDoSHH6YgjXXVJjbRobV19SQiWBxskJz+ET2bYrl86xZVIhdvU37vU0L1ff3rFvM3DVg/hr+cA18eL43Ch3ORGkxrJbYqFN22AstGPSTPU63MCJT/EOift8RzqhRRIhzMRq8SwValhaa7WfXspQQe8M1uPpMi+AaxvfpGMP7VKC0p20nPr32XP++DtqbikF7TTvrP01khpGJ7752tfja2QMlQoTkciG1JAgqI8LDOB3KeAlyfRIE6CtW/1dXgqUXVViWzuRsSg1Dn50pSFYBRwY44VSUBrPOqPDl4VSs7Gx9U5e9+gil530qw18P1PBx4tA24MNHaevSInAA7r7WQceYyUij2/8tRVDaI/e8j/wXc4yZUYkJSoYHb60frSj877S0OgZhFmYRPxWpxcDqxtviUxmuv9BPIFjR2/h161xOhBLuHNZfV6NVaREevNVhYHUxPOXWP3e1//rLhL/J3KhMtZbhz9saiAVAz0ofroEGfkI4i/VTCiSrGJw/mqX59ivTa+Cv+xpoGdDjo2kaqvgKEVLTfHlm2x01klTAkBpicBxQ6L9zVv6XZJsbHp92u7mpHe6/0YEBhuZHW2noJ0SSyvx5rAtTw0kM9Kuiv19q26fOECVCDrJ48WJs3boVwcHBGDVqlKPDIbmcn1veSIDSmtZcio03VJh1Rml22DEATDyairAXOnzdWIK2gSK0Ki3EyUgtVl1Vo3VpEXpVtq4+b68+QgCQqGQm88Ao1AyD9qbgfSrwayep0VDmkh4CdCkvwv77Giy9oDKagO/oIw2OPNTCQwa0L2P88csvQ/JkrItJk5q9uUr1yz78G63F3ddaVPqvg7VWxzD+SGqmqlf8l3F0PMtSVeH3m2r8r57EqMP0jJNKxCuBwTXEJsPGBRyH4bUk+C5UiSH/zTs1MhPNYgDQtIQIQg64EGu+qpisYpCJYJTYaXQMYw6l4mGc/votkmapFScxh5CaEvx8SYXvTiixKEhqGG1266UWG26oIRLoBxfw9oRr4CHTV2vSahMgQpsA81/Vxx9rkKRiWNVFZtTXp6G/ECIB8G+MFkcfadAu8MPjY/6bQRsAWgeYZhz8a9CylG3Tgz5VxOhTxXIi5CnjsNZBw+kpEXKQR48e4dKlS2jXrp2jQyHELrzkHL5pIsVXx5UATH+Z7rmnxi+X1ahfXIgfW+mTBaFAP6S+xqpkDP0rBbWLCS3O4ZNTNt5QY+F5JVqUEqGoC4e3Coa/HmjwPIlhbH2J2QnhlneUIex5MsYf1XdKrVlEiCfvddgbroFQAKztIjdp4knbydoRJjWSIGR/Khr/loxelcSQiTiERmqg1gLVCwsMI+Qy0rq0CDvuavDJdgU6lhFBLuZQ0p3D51Y2VXUoK0Lj35LRu7IYRV04nI3R4my0FqU8OMxtY76KN7SWGN+fVuJpIkNVX0GmO6G7yzi0DtAn4e9SmEk1KTRSg6H7U9EmQAQ/Nw5JKobDDzV49I6hZyX9yMeP/dBKitPRGiy5qML5WC0a+wvxMplh9z01UjXAkqAPa5QlqRiOPdbg04pim0y2WcxVgKnNpJh+UokOfyjQuZwIFQoJ8OK/509SAT0qiMyud3f0sQZCDuhWvuCkBwXnTHMZfokNZ2fbDKElJDcaU1+CFVdUiHxvnAhFx+sQsj8F7lJga0+50a9Zf3cBfusqQ/dtKei7U4GzQ5zNlvBzSp1iAlT0FuLwQw3eKhjcpBzqFRdiXT0xOlhYONXPTYCrw53x/Skl9j/Q4HSUFm5SDl3KizCliRT1ihv/En+XwvA0kaGxv9Bhlb8hNSVgDFh0QYWNN/S/0LuVF2F2ayk+3W79DN9Da4kRFa/D1ttqzDunX7aieUmh1YnQuAb69biWXFBh21sdXCQcgmuIMbuV5XmPCrsI0LGsCHvDNVYPmf/Y/9WR4OijFGy9rcaousbHKFdIgMYlhDgVpV9w1Umsbzqa2UKC/lVFJnMLAYCblMOZwc6Yc0aJHXf1i+/KRUCTEkJMbCQ1qtIcitAgVQN8kkF/usyY1lyK6oUFWHlVhXMxWhx8oIGTGKjqK8Tn1cQYbqZqFp/KsDdcg87lRPDPJXP85ASOsXSmuCRISEiAu7s74uPj4eZmu9WUe/fujR07duDnn3/GF198YbPjkpyhUChw7949VDw9Ak7xEY4Oh+Rx+++r0W1rCg72l+fbVentSccYyvychJfJDM8nuMLNzJIWGdHqGKr+mgyJECbrhtlb/10K7A3X4M1XrlZPImoPyy6qMOZwKs4MdspwDjBHUriXxb1mq3D337/x/mUMBn89Dy5uHib7Wfv9nWdTvlmzZqFRo0ZwcnKCh4eHxf02bNiAatWqQSaTwdfXF6NHj865INOhUCgAUEWIEKLvIFyjiICSoCzaeVeDJ+8ZBlYTZykJAvTNsgvaSXHjpc6wzEVOUGkZDkbo+/E4MglKUTPMOavEpxVFuToJsoc8e7YqlQq9evVCw4YNsW7dOrP7LFq0CAsXLsT8+fNRv359JCcnIzIyMmcDtYCaxgghvIVBuWuCubxi7lkl4lIYVl9VwVkMTGmavZGAHcuKsbS9zuJ8QvYgEXKIn2y71oasinyvw/DaYgTXyD0zPueUPJsIzZw5E4C+4mPOu3fv8N133+Gvv/5C69atDdurVauWE+FliE+EaK0xQgjJmin/KCEWAJV8BJjf1jZrV42pn7sm9MwpFX2EmNEih8et5xJ5NhHKyLFjx6DT6fD06VNUrFgRiYmJaNSoERYuXAh/f3+Lj1MqlVAqP0xalZCQYJf4XFxc4OHhAVdXy4vTEUIIsYxNd3wlheR9ebaPUEYeP34MnU6H2bNnY8mSJdi5cyfi4uLQtm1bqFSWJ3ibM2cO3N3dDbf0kqbsOHHiBN69e4cWLVrY5fiEEEIIyViuSoQmT54MjuPSvYWHh1t1LJ1OB7VajZ9//hlBQUFo0KAB/vzzT0RERCA0NNTi46ZMmYL4+HjDLSYmxlanRwghhJBcJlc1jU2YMAHBwcHp7hMQEGDVsYoWLQoAqFSpkmGbj48PvL29ER0dbfFxUqkUUmnBbCMmhBBCCppclQj5+PjAx8fHJsdq3LgxAOD+/fvw8/MDAMTFxeHNmzcoWbKkTZ6DEJa7iqqEEJLvGT53bTQNYp79FI+Ojsb169cRHR0NrVaL69ev4/r160hKSgIAlCtXDt26dcPYsWNx7tw53L59G4MGDUKFChXQsmVLB0dP8jqRSP8bQuVU2MGREEJIwaJyKgIA0KhMV7HPijybCE2bNg01a9bE9OnTkZSUhJo1a6JmzZq4cuWKYZ9Nmzahfv366NSpE5o3bw6xWIzDhw9DLKZJy0j2SCQSyCVivPHvAAbHTYJGCCEFCQOHN/5BUCW9gypVAfzXfzg7aImNDNhriQ2S98XFxeHJ48dwe3UR3jFHIFG8AAfrFqYkhBBiPQYBVE5F8MY/CPG+9fH8QRii7l6FSCJByOQFEEtMJ4K09vs7V/URIiQv8fLyglajwYPEOCTUaeDocAghJN9TJb3DmwdhSHjzDAnv3qB28w5mk6DMoESIkGzw8fXFs0cinN67FiKJFC5unhAKxaDWMkIIsSHGoFEroUpVQJWagvdvX6FQkeKoWr95tg9NiRAh2VStQUuIxRLcunQar55GQqvJwYWKCCGkIOE4iCVSVG/QEjWbtEWxkmWyfUhKhAjJJo7jUKlOE1Ss3RhJ8e+gVikzfhAhhJBM4zgOUrkznFxstzwVJUKE2AjHcXD18HJ0GIQQQjIhzw6fJ4QQQgjJLkqECCGEEFJgUSJECCGEkAKLEiFCCCGEFFjUWToD/MTbCQkJDo6EEEIIIdbiv7czWkCDEqEMJCYmAgD8/f0dHAkhhBBCMisxMRHu7u4W76e1xjKg0+nw7NkzuLq6Zntht7QSEhLg7++PmJiYfLuGWX4/x/x+fgCdY36Q388PyP/nmN/PD7DPOTLGkJiYiGLFikEgsNwTiCpCGRAIBPDz87Pb8d3c3PLthc3L7+eY388PoHPMD/L7+QH5/xzz+/kBtj/H9CpBPOosTQghhJACixIhQgghhBRYlAg5iFQqxfTp0yGVSh0dit3k93PM7+cH0DnmB/n9/ID8f475/fwAx54jdZYmhBBCSIFFFSFCCCGEFFiUCBFCCCGkwKJEiBBCCCEFFiVChBBCCCmwKBFykOXLl6NUqVKQyWSoX78+Ll265OiQbOb06dPo0qULihUrBo7jsHfvXkeHZFNz5sxB3bp14erqCl9fX3Tv3h337993dFg29euvv6JatWqGyc0aNmyIQ4cOOTosu5k7dy44jsOXX37p6FBsZsaMGeA4zuhWoUIFR4dlU0+fPsWAAQNQqFAhyOVyVK1aFVeuXHF0WDZTqlQpk/eQ4ziMHj3a0aHZhFarxdSpU/H/7d17SJN9H8fxz9rdKlPSTuaUaZlaapklSlp0pBKJzlkMWtmfkywzOv0hFWkGhXbAsjtmVGJR2flklkbQYVkL7WRng8oRlXmgGdvv+ePhkXt4wwN1bb90nxfsj/0YF++LIX537dp1DR48GL169UJoaCi2bNnyf+8NpjQOQhIcO3YMmZmZyM7OxoMHDxATE4Pp06fDarXKTlNES0sLYmJisHfvXtkpLlFVVQWj0Yg7d+6gvLwcP3/+xLRp09DS0iI7TTFBQUHYtm0bqqurcf/+fUyePBmzZs3C48ePZacpzmw2Y//+/Rg5cqTsFMVFRUXh48eP7Y9bt27JTlLM169fkZSUhO7du+PSpUt48uQJduzYAT8/P9lpijGbzU7vX3l5OQBgwYIFksuUkZeXh8LCQuzZswdPnz5FXl4etm/fjt27d7s3RJDbxcfHC6PR2P7cbrcLrVYrcnNzJVa5BgBRVlYmO8OlrFarACCqqqpkp7iUn5+f+Pvvv2VnKKqpqUmEhYWJ8vJyMWHCBJGRkSE7STHZ2dkiJiZGdobLrF27VowbN052hltlZGSI0NBQ4XA4ZKcoIiUlRaSlpTmtzZ07V+j1erd28IiQm7W1taG6uhpTp05tX+vWrRumTp2K27dvSyyjX9XY2AgA6Nu3r+QS17Db7SgtLUVLSwvGjh0rO0dRRqMRKSkpTn+PXcmLFy+g1WoxZMgQ6PV61NfXy05SzNmzZxEXF4cFCxZg4MCBiI2NxYEDB2RnuUxbWxuOHDmCtLQ0RW8ALlNiYiIqKipQV1cHAHj06BFu3bqF5ORkt3bwpqtu9vnzZ9jtdvj7+zut+/v749mzZ5Kq6Fc5HA6sXLkSSUlJiI6Olp2jqJqaGowdOxY/fvyAt7c3ysrKEBkZKTtLMaWlpXjw4AHMZrPsFJdISEhAcXExIiIi8PHjR2zatAnjx49HbW0tfHx8ZOf9ttevX6OwsBCZmZnYsGEDzGYzVqxYAY1GA4PBIDtPcadPn8a3b9+wdOlS2SmKWbduHb5//45hw4ZBrVbDbrdj69at0Ov1bu3gIET0G4xGI2pra7vUuRf/ExERAYvFgsbGRpw4cQIGgwFVVVVdYhh6//49MjIyUF5ejp49e8rOcYl/fqoeOXIkEhISEBwcjOPHj2P58uUSy5ThcDgQFxeHnJwcAEBsbCxqa2uxb9++LjkIHTx4EMnJydBqtbJTFHP8+HEcPXoUJSUliIqKgsViwcqVK6HVat36HnIQcrP+/ftDrVajoaHBab2hoQGDBg2SVEW/Ij09HefPn8fNmzcRFBQkO0dxGo0GQ4cOBQCMGTMGZrMZBQUF2L9/v+Sy31ddXQ2r1YrRo0e3r9ntdty8eRN79uyBzWaDWq2WWKg8X19fhIeH4+XLl7JTFBEQENBhKB8+fDhOnjwpqch13r17h2vXruHUqVOyUxS1Zs0arFu3DosWLQIAjBgxAu/evUNubq5bByGeI+RmGo0GY8aMQUVFRfuaw+FARUVFlzv/oqsSQiA9PR1lZWW4fv06Bg8eLDvJLRwOB2w2m+wMRUyZMgU1NTWwWCztj7i4OOj1elgsli43BAFAc3MzXr16hYCAANkpikhKSupw2Yq6ujoEBwdLKnIdk8mEgQMHIiUlRXaKolpbW9Gtm/MYolar4XA43NrBI0ISZGZmwmAwIC4uDvHx8cjPz0dLSwuWLVsmO00Rzc3NTp8637x5A4vFgr59+0Kn00ksU4bRaERJSQnOnDkDHx8ffPr0CQDQp08f9OrVS3KdMtavX4/k5GTodDo0NTWhpKQElZWVuHLliuw0Rfj4+HQ4p6t3797o169flznXKysrCzNnzkRwcDA+fPiA7OxsqNVqLF68WHaaIlatWoXExETk5ORg4cKFuHfvHoqKilBUVCQ7TVEOhwMmkwkGgwF//dW1/mXPnDkTW7duhU6nQ1RUFB4+fIidO3ciLS3NvSFu/Y0atdu9e7fQ6XRCo9GI+Ph4cefOHdlJirlx44YA0OFhMBhkpyni3/YNgDCZTLLTFJOWliaCg4OFRqMRAwYMEFOmTBFXr16VneVSXe3n86mpqSIgIEBoNBoRGBgoUlNTxcuXL2VnKercuXMiOjpa9OjRQwwbNkwUFRXJTlLclStXBADx/Plz2SmK+/79u8jIyBA6nU707NlTDBkyRGzcuFHYbDa3dqiEcPMlHImIiIj+EDxHiIiIiDwWByEiIiLyWByEiIiIyGNxECIiIiKPxUGIiIiIPBYHISIiIvJYHISIiIjIY3EQIqIup7KyEiqVCt++fZOdQkR/OF5QkYg6vYkTJ2LUqFHIz88HALS1teHLly/w9/eHSqWSG0dEf7SudeMSIiL89+bGgwYNkp1BRJ0Avxojok5t6dKlqKqqQkFBAVQqFVQqFYqLi52+GisuLoavry/Onz+PiIgIeHl5Yf78+WhtbcWhQ4cQEhICPz8/rFixAna7vX3bNpsNWVlZCAwMRO/evZGQkIDKyko5O0pELsEjQkTUqRUUFKCurg7R0dHYvHkzAODx48cdXtfa2opdu3ahtLQUTU1NmDt3LubMmQNfX19cvHgRr1+/xrx585CUlITU1FQAQHp6Op48eYLS0lJotVqUlZVhxowZqKmpQVhYmFv3k4hcg4MQEXVqffr0gUajgZeXV/vXYc+ePevwup8/f6KwsBChoaEAgPnz5+Pw4cNoaGiAt7c3IiMjMWnSJNy4cQOpqamor6+HyWRCfX09tFotACArKwuXL1+GyWRCTk6O+3aSiFyGgxAReQQvL6/2IQgA/P39ERISAm9vb6c1q9UKAKipqYHdbkd4eLjTdmw2G/r16+eeaCJyOQ5CROQRunfv7vRcpVL965rD4QAANDc3Q61Wo7q6Gmq12ul1/xyeiKhz4yBERJ2eRqNxOslZCbGxsbDb7bBarRg/fryi2yaiPwd/NUZEnV5ISAju3r2Lt2/f4vPnz+1HdX5HeHg49Ho9lixZglOnTuHNmze4d+8ecnNzceHCBQWqiehPwEGIiDq9rKwsqNVqREZGYsCAAaivr1dkuyaTCUuWLMHq1asRERGB2bNnw2w2Q6fTKbJ9IpKPV5YmIiIij8UjQkREROSxOAgRERGRx+IgRERERB6LgxARERF5LA5CRERE5LE4CBEREZHH4iBEREREHouDEBEREXksDkJERETksTgIERERkcfiIEREREQei4MQEREReaz/AAmaO8jCPqHjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import mpmath as mp\n",
    "import csv\n",
    "\n",
    "def file_reader(filename):\n",
    "    with open(filename) as file:\n",
    "        reader = csv.reader(file, delimiter=\" \")\n",
    "        data  = list(zip(*reader))\n",
    "        # data is a tuple of strings. Tuples are immutable, and we need to perform math on\n",
    "        #   the data, so here we convert tuple to lists of floats:\n",
    "        data0 = []\n",
    "        data1 = []\n",
    "        for i in range(len(data[0])):\n",
    "            data0.append(float(data[0][i]))\n",
    "            data1.append(float(data[1][i]))\n",
    "        return data0,data1\n",
    "\n",
    "first_col48,second_col48 = file_reader(os.path.join(outdir,\"out48.txt\"))\n",
    "first_col64,second_col64 = file_reader(os.path.join(outdir,\"out64.txt\"))\n",
    "first_col96,second_col96 = file_reader(os.path.join(outdir,\"out96.txt\"))\n",
    "\n",
    "for i in range(len(second_col64)):\n",
    "    # data64 = data48*(64/48)**4\n",
    "    # -> log10(data64) = log10(data48) + 4*log(64/48)\n",
    "    second_col64[i] += 4*mp.log10(64./48.)\n",
    "for i in range(len(second_col96)):\n",
    "    # data96 = data48*(96/48)**4\n",
    "    # -> log10(data96) = log10(data48) + 4*log(96/48)\n",
    "    second_col96[i] += 4*mp.log10(96./48.)\n",
    "\n",
    "# https://matplotlib.org/gallery/text_labels_and_annotations/legend.html#sphx-glr-gallery-text-labels-and-annotations-legend-py\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "plt.title(\"Plot Demonstrating 4th-order Convergence\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(r\"$\\log_{10}$(Relative error)\")\n",
    "\n",
    "ax.plot(first_col48, second_col48, 'k--', label='Nx = 48')\n",
    "ax.plot(first_col64, second_col64, 'k-', label='Nx = 64, mult by (64/48)^4')\n",
    "ax.plot(first_col96, second_col96, 'k.', label='Nx = 96, mult by (96/48)^4')\n",
    "\n",
    "legend = ax.legend(loc='lower right', shadow=True, fontsize='x-large')\n",
    "legend.get_frame().set_facecolor('C1')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 4: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-Start_to_Finish-ScalarWave.pdf](Tutorial-Start_to_Finish-ScalarWave.pdf). (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-10-04T19:30:09.819476Z",
     "iopub.status.busy": "2021-10-04T19:30:09.819219Z",
     "iopub.status.idle": "2021-10-04T19:30:12.061343Z",
     "shell.execute_reply": "2021-10-04T19:30:12.060914Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish-ScalarWave.tex, and compiled LaTeX file to\n",
      "    PDF file Tutorial-Start_to_Finish-ScalarWave.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Start_to_Finish-ScalarWave\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}