{
 "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: Scalar Field Collapse\n",
    "\n",
    "## Authors: Leonardo Werneck & Zachariah B. Etienne\n",
    "\n",
    "## This module sets up spherically symmetric, time-symmetric initial data for a scalar field collapse in Spherical coordinates, as [documented in this NRPy+ module](Tutorial-ADM_Initial_Data-ScalarField.ipynb) (the initial data is shown to satisfy the Hamiltonian constraint [in this tutorial module](Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_ScalarField_initial_data.ipynb)), which is then evolved forward in time. The aim is to reproduce the results from [Akbarian & Choptuik (2014)]( https://arxiv.org/pdf/1508.01614.pdf)  and [Baumgarte (2018)](https://arxiv.org/abs/1807.10342) (which used a similar approach), demonstrating that the Hamiltonian constraint violation during the simulation also converges to zero with increasing numerical resolution.\n",
    "\n",
    "### **Results from this tutorial notebook have been used in the paper [Werneck *et al.* (2021)](https://arxiv.org/pdf/2106.06553.pdf)**\n",
    "\n",
    "The entire algorithm is outlined below, with NRPy+-based components highlighted in <font color='green'>green</font>.\n",
    "\n",
    "1. Allocate memory for gridfunctions, including temporary storage for the RK4 time integration.\n",
    "1. <font color='green'>Set gridfunction values to initial data.</font>\n",
    "1. Evolve the system forward in time using RK4 time integration. At each RK4 substep, do the following:\n",
    "    1. <font color='green'>Evaluate BSSN RHS expressions.</font>\n",
    "    1. Apply singular, curvilinear coordinate boundary conditions [*a la* the SENR/NRPy+ paper](https://arxiv.org/abs/1712.07658)\n",
    "    1. <font color='green'>Apply constraints on conformal 3-metric: $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$</font>\n",
    "1. At the end of each iteration in time, output the <font color='green'>Hamiltonian constraint violation</font>.\n",
    "1. Repeat above steps at two numerical resolutions to confirm convergence to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "* [Akbarian & Choptuik (2015)](https://arxiv.org/pdf/1508.01614.pdf) (Useful to understand the theoretical framework)\n",
    "* [Baumgarte (2018)](https://arxiv.org/pdf/1807.10342.pdf) (Useful to understand the theoretical framework)\n",
    "* [Baumgarte & Shapiro's Numerical Relativity](https://books.google.com.br/books/about/Numerical_Relativity.html?id=dxU1OEinvRUC&redir_esc=y): Section 6.2.2 (Useful to understand how to solve the Hamiltonian constraint)\n",
    "\n",
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "1. [Step 1](#nrpy_core) Set core NRPy+ parameters for numerical grids and reference metric\n",
    "    1. [Step 1.a](#cfl) Output needed C code for finding the minimum proper distance between grid points, needed for [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673)-limited timestep\n",
    "1. [Step 2](#initial_data) Set up ADM initial data for the Scalar Field \n",
    "1. [Step 3](#adm_id_spacetime) Convert ADM initial data to BSSN-in-curvilinear coordinates\n",
    "1. [Step 4](#bssn) Output C code for BSSN spacetime evolution\n",
    "    1. [Step 4.a](#bssnrhs) Set up the BSSN and ScalarField right-hand-side (RHS) expressions, and add the *rescaled* $T^{\\mu\\nu}$ source terms\n",
    "    1. [Step 4.b](#hamconstraint) Output the Hamiltonian constraint\n",
    "    1. [Step 4.c](#enforce3metric) Enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint\n",
    "    1. [Step 4.d](#ccodegen) Generate C code kernels for BSSN expressions, in parallel if possible\n",
    "    1. [Step 4.e](#cparams_rfm_and_domainsize) Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h`\n",
    "1. [Step 5](#bc_functs) Set up boundary condition functions for chosen singular, curvilinear coordinate system\n",
    "1. [Step 6](#main_ccode) The main C code: `ScalarFieldCollapse_Playground.c`\n",
    "1. [Step 7](#visualization) Visualization\n",
    "    1. [Step 7.a](#install_download) Install `scipy` and download `ffmpeg` if they are not yet installed/downloaded\n",
    "    1. [Step 7.b](#movie_dynamics) Dynamics of the solution\n",
    "        1. [Step 7.b.i](#genimages) Generate images for visualization animation\n",
    "        1. [Step 7.b.ii](#gemnvideo) Generate visualization animation\n",
    "    1. [Step 7.c](#convergence) Convergence of constraint violation\n",
    "1. [Step 8](#output_to_pdf) Output this module as $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy_core'></a>\n",
    "\n",
    "# Step 1: Set core NRPy+ parameters for numerical grids and reference metric  \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy_core}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:12:41.106640Z",
     "iopub.status.busy": "2021-06-15T10:12:41.105929Z",
     "iopub.status.idle": "2021-06-15T10:12:41.531078Z",
     "shell.execute_reply": "2021-06-15T10:12:41.530547Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step P1: Import needed NRPy+ core modules:\n",
    "from outputC import lhrh,outputC,outCfunction  # NRPy+: Core C code output module\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import finite_difference as fin  # NRPy+: Finite difference C code generation module\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface\n",
    "import shutil, os, sys           # Standard Python modules for multiplatform OS-level functions\n",
    "\n",
    "# Step P2: Create C code output directory:\n",
    "Ccodesdir = os.path.join(\"BSSN_ScalarFieldCollapse_Ccodes\")\n",
    "# First remove C code output directory if it exists\n",
    "# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
    "# !rm -r ScalarWaveCurvilinear_Playground_Ccodes\n",
    "shutil.rmtree(Ccodesdir, ignore_errors=True)\n",
    "# Then create a fresh directory\n",
    "cmd.mkdir(Ccodesdir)\n",
    "\n",
    "# Step P3: Create executable output directory:\n",
    "outdir = os.path.join(Ccodesdir,\"output\")\n",
    "cmd.mkdir(outdir)\n",
    "\n",
    "# Step 1: Set the spatial dimension parameter\n",
    "#         to three this time, and then read\n",
    "#         the parameter as DIM.\n",
    "par.set_parval_from_str(\"grid::DIM\",3)\n",
    "DIM = par.parval_from_str(\"grid::DIM\")\n",
    "\n",
    "# Step 2: Set some core parameters, including CoordSystem MoL timestepping algorithm,\n",
    "#                                 FD order, floating point precision, and CFL factor:\n",
    "# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,\n",
    "#              SymTP, SinhSymTP\n",
    "CoordSystem     = \"Spherical\"\n",
    "\n",
    "# Step 2.a: Set defaults for Coordinate system parameters.\n",
    "#           These are perhaps the most commonly adjusted parameters,\n",
    "#           so we enable modifications at this high level.\n",
    "domain_size     = 32\n",
    "\n",
    "# sinh_width sets the default value for:\n",
    "#   * SinhSpherical's params.SINHW\n",
    "#   * SinhCylindrical's params.SINHW{RHO,Z}\n",
    "#   * SinhSymTP's params.SINHWAA\n",
    "sinh_width      = 0.2 # If Sinh* coordinates chosen\n",
    "\n",
    "# sinhv2_const_dr sets the default value for:\n",
    "#   * SinhSphericalv2's params.const_dr\n",
    "#   * SinhCylindricalv2's params.const_d{rho,z}\n",
    "sinhv2_const_dr = 0.05# If Sinh*v2 coordinates chosen\n",
    "\n",
    "# SymTP_bScale sets the default value for:\n",
    "#   * SinhSymTP's params.bScale\n",
    "SymTP_bScale    = 0.5 # If SymTP chosen\n",
    "\n",
    "# Step 2.b: Set the order of spatial and temporal derivatives;\n",
    "#           the core data type, and the CFL factor.\n",
    "# RK_method choices include: Euler, \"RK2 Heun\", \"RK2 MP\", \"RK2 Ralston\", RK3, \"RK3 Heun\", \"RK3 Ralston\",\n",
    "#              SSPRK3, RK4, DP5, DP5alt, CK5, DP6, L6, DP8\n",
    "RK_method = \"RK4\"\n",
    "FD_order  = 4        # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
    "REAL      = \"double\" # Best to use double here.\n",
    "CFL_FACTOR= 0.5\n",
    "\n",
    "# Set the lapse & shift conditions\n",
    "LapseCondition  = \"OnePlusLog\"\n",
    "ShiftCondition  = \"GammaDriving2ndOrder_Covariant\"\n",
    "\n",
    "# Step 3: Generate Runge-Kutta-based (RK-based) timestepping code.\n",
    "#       As described above the Table of Contents, this is a 3-step process:\n",
    "#       3.A: Evaluate RHSs (RHS_string)\n",
    "#       3.B: Apply boundary conditions (post_RHS_string, pt 1)\n",
    "#       3.C: Enforce det(gammabar) = det(gammahat) constraint (post_RHS_string, pt 2)\n",
    "import MoLtimestepping.C_Code_Generation 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(Ccodesdir,\"MoLtimestepping/\"))\n",
    "MoL.MoL_C_Code_Generation(RK_method,\n",
    "    RHS_string      = \"\"\"\n",
    "Ricci_eval(&rfmstruct, &params, RK_INPUT_GFS, auxevol_gfs);\n",
    "rhs_eval(&rfmstruct, &params, auxevol_gfs, RK_INPUT_GFS, RK_OUTPUT_GFS);\"\"\",\n",
    "    post_RHS_string = \"\"\"\n",
    "apply_bcs_curvilinear(&params, &bcstruct, NUM_EVOL_GFS, evol_gf_parity, RK_OUTPUT_GFS);\n",
    "enforce_detgammahat_constraint(&rfmstruct, &params,                     RK_OUTPUT_GFS);\\n\"\"\",\n",
    "    outdir = os.path.join(Ccodesdir,\"MoLtimestepping/\"))\n",
    "\n",
    "# Step 4: Set the coordinate system for the numerical grid\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",CoordSystem)\n",
    "rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.\n",
    "\n",
    "# Step 5: Set the finite differencing order to FD_order (set above).\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", FD_order)\n",
    "\n",
    "# Step 6: Copy SIMD/SIMD_intrinsics.h to $Ccodesdir/SIMD/SIMD_intrinsics.h\n",
    "cmd.mkdir(os.path.join(Ccodesdir,\"SIMD\"))\n",
    "shutil.copy(os.path.join(\"SIMD/\")+\"SIMD_intrinsics.h\",os.path.join(Ccodesdir,\"SIMD/\"))\n",
    "\n",
    "# Step 7: Impose spherical symmetry by demanding that all\n",
    "#         derivatives in the angular directions vanish\n",
    "par.set_parval_from_str(\"indexedexp::symmetry_axes\",\"12\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cfl'></a>\n",
    "\n",
    "## Step 1.a: Output needed C code for finding the minimum proper distance between grid points, needed for [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673)-limited timestep \\[Back to [top](#toc)\\]\n",
    "$$\\label{cfl}$$\n",
    "\n",
    "In order for our explicit-timestepping numerical solution to the scalar wave equation to be stable, it must satisfy the [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673) condition:\n",
    "$$\n",
    "\\Delta t \\le \\frac{\\min(ds_i)}{c},\n",
    "$$\n",
    "where $c$ is the wavespeed, and\n",
    "$$ds_i = h_i \\Delta x^i$$ \n",
    "is the proper distance between neighboring gridpoints in the $i$th direction (in 3D, there are 3 directions), $h_i$ is the $i$th reference metric scale factor, and $\\Delta x^i$ is the uniform grid spacing in the $i$th direction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:12:41.537996Z",
     "iopub.status.busy": "2021-06-15T10:12:41.536230Z",
     "iopub.status.idle": "2021-06-15T10:12:41.538945Z",
     "shell.execute_reply": "2021-06-15T10:12:41.538484Z"
    }
   },
   "outputs": [],
   "source": [
    "# Output the find_timestep() function to a C file.\n",
    "rfm.out_timestep_func_to_file(os.path.join(Ccodesdir,\"find_timestep.h\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initial_data'></a>\n",
    "\n",
    "# Step 2: Set up ADM initial data for the Scalar Field \\[Back to [top](#toc)\\]\n",
    "$$\\label{initial_data}$$\n",
    "\n",
    "As documented [in the scalar field Gaussian pulse initial data NRPy+ tutorial notebook](TTutorial-ADM_Initial_Data-ScalarField.ipynb), we will now set up the scalar field initial data, storing the densely-sampled result to file.\n",
    "\n",
    "The initial data function `ScalarField_InitialData` requires `SciPy`, so let's make sure it's installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:12:41.542000Z",
     "iopub.status.busy": "2021-06-15T10:12:41.541591Z",
     "iopub.status.idle": "2021-06-15T10:12:42.134746Z",
     "shell.execute_reply": "2021-06-15T10:12:42.134281Z"
    }
   },
   "outputs": [],
   "source": [
    "!pip install scipy numpy > /dev/null"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next call the `ScalarField_InitialData()` function from the [ScalarField/ScalarField_InitialData.py](../edit/ScalarField/ScalarField_InitialData.py) NRPy+ module (see the [tutorial notebook](Tutorial-ADM_Initial_Data-ScalarField.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:12:42.139833Z",
     "iopub.status.busy": "2021-06-15T10:12:42.139297Z",
     "iopub.status.idle": "2021-06-15T10:12:42.470689Z",
     "shell.execute_reply": "2021-06-15T10:12:42.470439Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generated the ADM initial data for the gravitational collapse \n",
      "of a massless scalar field in Spherical coordinates.\n",
      "\n",
      "Type of initial condition: Scalar field: \"Gaussian\" Shell\n",
      "                         ADM quantities: Time-symmetric\n",
      "                        Lapse condition: Pre-collapsed\n",
      "Parameters: amplitude         = 0.4,\n",
      "            center            = 0,\n",
      "            width             = 1,\n",
      "            domain size       = 35.2,\n",
      "            number of points  = 30000,\n",
      "            Initial data file = BSSN_ScalarFieldCollapse_Ccodes/output/SFID.txt.\n",
      "\n",
      "Output C function ID_scalarfield_ADM_quantities() to file BSSN_ScalarFieldCollapse_Ccodes/ID_scalarfield_ADM_quantities.h\n",
      "Output C function ID_scalarfield_spherical() to file BSSN_ScalarFieldCollapse_Ccodes/ID_scalarfield_spherical.h\n",
      "Output C function ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2() to file BSSN_ScalarFieldCollapse_Ccodes/ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2.h\n",
      "Output C function ID_scalarfield() to file BSSN_ScalarFieldCollapse_Ccodes/ID_scalarfield.h\n"
     ]
    }
   ],
   "source": [
    "# Step 2.a: Import necessary Python and NRPy+ modules\n",
    "import ScalarField.ScalarField_InitialData as sfid\n",
    "\n",
    "# Step 2.b: Set the initial data parameters\n",
    "outputfilename  = os.path.join(outdir,\"SFID.txt\")\n",
    "ID_Family       = \"Gaussian_pulse\"\n",
    "pulse_amplitude = 0.4\n",
    "pulse_center    = 0\n",
    "pulse_width     = 1\n",
    "Nr              = 30000\n",
    "rmax            = domain_size*1.1\n",
    "\n",
    "# Step 2.c: Generate the initial data\n",
    "sfid.ScalarField_InitialData(outputfilename,ID_Family,\n",
    "                             pulse_amplitude,pulse_center,pulse_width,Nr,rmax)\n",
    "\n",
    "# Step 2.d: Generate the needed C code\n",
    "sfid.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(Ccodesdir=Ccodesdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='adm_id_spacetime'></a>\n",
    "\n",
    "# Step 3: Convert ADM initial data to BSSN-in-curvilinear coordinates \\[Back to [top](#toc)\\]\n",
    "$$\\label{adm_id_spacetime}$$\n",
    "\n",
    "This is an automated process, taken care of by [`BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear`](../edit/BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear.py), and documented [in this tutorial notebook](Tutorial-ADM_Initial_Data-Converting_Numerical_ADM_Spherical_or_Cartesian_to_BSSNCurvilinear.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:12:42.472676Z",
     "iopub.status.busy": "2021-06-15T10:12:42.472417Z",
     "iopub.status.idle": "2021-06-15T10:12:42.905412Z",
     "shell.execute_reply": "2021-06-15T10:12:42.904960Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function ID_BSSN_lambdas() to file BSSN_ScalarFieldCollapse_Ccodes/ID_BSSN_lambdas.h\n",
      "Output C function ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs() to file BSSN_ScalarFieldCollapse_Ccodes/ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs.h\n",
      "Output C function ID_BSSN__ALL_BUT_LAMBDAs() to file BSSN_ScalarFieldCollapse_Ccodes/ID_BSSN__ALL_BUT_LAMBDAs.h\n"
     ]
    }
   ],
   "source": [
    "import BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear as AtoBnum\n",
    "AtoBnum.Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear(\"Spherical\",\"ID_scalarfield_ADM_quantities\",\n",
    "                                                               Ccodesdir=Ccodesdir,loopopts=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssn'></a>\n",
    "\n",
    "# Step 4: Output C code for BSSN spacetime evolution \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssn}$$\n",
    "\n",
    "<a id='bssnrhs'></a>\n",
    "\n",
    "## Step 4.a: Set up the BSSN and ScalarField right-hand-side (RHS) expressions, and add the *rescaled* $T^{\\mu\\nu}$ source terms \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssnrhs}$$\n",
    "\n",
    "`BSSN.BSSN_RHSs()` sets up the RHSs assuming a spacetime vacuum: $T^{\\mu\\nu}=0$. (This might seem weird, but remember that, for example, *spacetimes containing only single or binary black holes are vacuum spacetimes*.) Here, using the [`BSSN.BSSN_stress_energy_source_terms`](../edit/BSSN/BSSN_stress_energy_source_terms.py) ([**tutorial**](Tutorial-BSSN_stress_energy_source_terms.ipynb)) NRPy+ module, we add the $T^{\\mu\\nu}$ source terms to these equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:12:42.924208Z",
     "iopub.status.busy": "2021-06-15T10:12:42.922875Z",
     "iopub.status.idle": "2021-06-15T10:12:48.090437Z",
     "shell.execute_reply": "2021-06-15T10:12:48.089955Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating symbolic expressions for BSSN RHSs...\n",
      "(BENCH) Finished BSSN symbolic expressions in 7.168536186218262 seconds.\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "import BSSN.BSSN_RHSs as rhs\n",
    "import BSSN.BSSN_gauge_RHSs as gaugerhs\n",
    "par.set_parval_from_str(\"BSSN.BSSN_gauge_RHSs::LapseEvolutionOption\", LapseCondition)\n",
    "par.set_parval_from_str(\"BSSN.BSSN_gauge_RHSs::ShiftEvolutionOption\", ShiftCondition)\n",
    "\n",
    "print(\"Generating symbolic expressions for BSSN RHSs...\")\n",
    "start = time.time()\n",
    "# Enable rfm_precompute infrastructure, which results in\n",
    "#   BSSN RHSs that are free of transcendental functions,\n",
    "#   even in curvilinear coordinates, so long as\n",
    "#   ConformalFactor is set to \"W\" (default).\n",
    "cmd.mkdir(os.path.join(Ccodesdir,\"rfm_files/\"))\n",
    "par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\",\"True\")\n",
    "par.set_parval_from_str(\"reference_metric::rfm_precompute_Ccode_outdir\",os.path.join(Ccodesdir,\"rfm_files/\"))\n",
    "\n",
    "# Evaluate BSSN + BSSN gauge RHSs with rfm_precompute enabled:\n",
    "import BSSN.BSSN_quantities as Bq\n",
    "par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\",\"True\")\n",
    "\n",
    "rhs.BSSN_RHSs()\n",
    "\n",
    "# Evaluate the Scalar Field RHSs\n",
    "import ScalarField.ScalarField_RHSs as sfrhs\n",
    "sfrhs.ScalarField_RHSs()\n",
    "\n",
    "# Compute ScalarField T^{\\mu\\nu}\n",
    "# Compute the scalar field energy-momentum tensor\n",
    "import ScalarField.ScalarField_Tmunu as sfTmunu\n",
    "sfTmunu.ScalarField_Tmunu()\n",
    "T4UU = sfTmunu.T4UU\n",
    "\n",
    "import BSSN.BSSN_stress_energy_source_terms as Bsest\n",
    "Bsest.BSSN_source_terms_for_BSSN_RHSs(T4UU)\n",
    "rhs.trK_rhs += Bsest.sourceterm_trK_rhs\n",
    "for i in range(DIM):\n",
    "    # Needed for Gamma-driving shift RHSs:\n",
    "    rhs.Lambdabar_rhsU[i] += Bsest.sourceterm_Lambdabar_rhsU[i]\n",
    "    # Needed for BSSN RHSs:\n",
    "    rhs.lambda_rhsU[i]    += Bsest.sourceterm_lambda_rhsU[i]\n",
    "    for j in range(DIM):\n",
    "        rhs.a_rhsDD[i][j] += Bsest.sourceterm_a_rhsDD[i][j]\n",
    "\n",
    "gaugerhs.BSSN_gauge_RHSs()\n",
    "\n",
    "# We use betaU as our upwinding control vector:\n",
    "Bq.BSSN_basic_tensors()\n",
    "betaU = Bq.betaU\n",
    "\n",
    "import BSSN.Enforce_Detgammahat_Constraint as EGC\n",
    "enforce_detg_constraint_symb_expressions = EGC.Enforce_Detgammahat_Constraint_symb_expressions()\n",
    "\n",
    "# Next compute Ricci tensor\n",
    "par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\",\"False\")\n",
    "Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()\n",
    "\n",
    "# Now register the Hamiltonian as a gridfunction.\n",
    "H = gri.register_gridfunctions(\"AUX\",\"H\")\n",
    "\n",
    "# Then define the Hamiltonian constraint and output the optimized C code.\n",
    "import BSSN.BSSN_constraints as bssncon\n",
    "bssncon.BSSN_constraints(add_T4UUmunu_source_terms=False)\n",
    "Bsest.BSSN_source_terms_for_BSSN_constraints(T4UU)\n",
    "bssncon.H += Bsest.sourceterm_H\n",
    "\n",
    "# Add Kreiss-Oliger dissipation\n",
    "diss_strength = par.Cparameters(\"REAL\",\"ScalarFieldCollapse\",[\"diss_strength\"],0.1)\n",
    "\n",
    "alpha_dKOD   = ixp.declarerank1(\"alpha_dKOD\")\n",
    "cf_dKOD      = ixp.declarerank1(\"cf_dKOD\")\n",
    "trK_dKOD     = ixp.declarerank1(\"trK_dKOD\")\n",
    "sf_dKOD      = ixp.declarerank1(\"sf_dKOD\")\n",
    "sfM_dKOD     = ixp.declarerank1(\"sfM_dKOD\")\n",
    "betU_dKOD    = ixp.declarerank2(\"betU_dKOD\",\"nosym\")\n",
    "vetU_dKOD    = ixp.declarerank2(\"vetU_dKOD\",\"nosym\")\n",
    "lambdaU_dKOD = ixp.declarerank2(\"lambdaU_dKOD\",\"nosym\")\n",
    "aDD_dKOD     = ixp.declarerank3(\"aDD_dKOD\",\"sym01\")\n",
    "hDD_dKOD     = ixp.declarerank3(\"hDD_dKOD\",\"sym01\")\n",
    "\n",
    "for k in range(3):\n",
    "    gaugerhs.alpha_rhs += diss_strength*alpha_dKOD[k]*rfm.ReU[k]\n",
    "    rhs.cf_rhs         += diss_strength*   cf_dKOD[k]*rfm.ReU[k]\n",
    "    rhs.trK_rhs        += diss_strength*  trK_dKOD[k]*rfm.ReU[k]\n",
    "    sfrhs.sf_rhs       += diss_strength*   sf_dKOD[k]*rfm.ReU[k]\n",
    "    sfrhs.sfM_rhs      += diss_strength*  sfM_dKOD[k]*rfm.ReU[k]\n",
    "    for i in range(3):\n",
    "        if \"2ndOrder\" in ShiftCondition:\n",
    "            gaugerhs.bet_rhsU[i] += diss_strength*   betU_dKOD[i][k]*rfm.ReU[k]\n",
    "        gaugerhs.vet_rhsU[i]     += diss_strength*   vetU_dKOD[i][k]*rfm.ReU[k]\n",
    "        rhs.lambda_rhsU[i]       += diss_strength*lambdaU_dKOD[i][k]*rfm.ReU[k]\n",
    "        for j in range(3):\n",
    "            rhs.a_rhsDD[i][j] += diss_strength*aDD_dKOD[i][j][k]*rfm.ReU[k]\n",
    "            rhs.h_rhsDD[i][j] += diss_strength*hDD_dKOD[i][j][k]*rfm.ReU[k]\n",
    "\n",
    "# Now that we are finished with all the rfm hatted\n",
    "#           quantities in generic precomputed functional\n",
    "#           form, let's restore them to their closed-\n",
    "#           form expressions.\n",
    "par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\",\"False\") # Reset to False to disable rfm_precompute.\n",
    "rfm.ref_metric__hatted_quantities()\n",
    "end = time.time()\n",
    "print(\"(BENCH) Finished BSSN symbolic expressions in \"+str(end-start)+\" seconds.\")\n",
    "\n",
    "def BSSN_plus_ScalarField_RHSs():\n",
    "    print(\"Generating C code for BSSN RHSs in \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" coordinates.\")\n",
    "    start = time.time()\n",
    "\n",
    "    # Construct the left-hand sides and right-hand-side expressions for all BSSN RHSs\n",
    "    lhs_names = [        \"alpha\",       \"cf\",       \"trK\",         \"sf\",         \"sfM\"   ]\n",
    "    rhs_exprs = [gaugerhs.alpha_rhs, rhs.cf_rhs, rhs.trK_rhs, sfrhs.sf_rhs, sfrhs.sfM_rhs]\n",
    "    for i in range(3):\n",
    "        lhs_names.append(        \"betU\"+str(i))\n",
    "        rhs_exprs.append(gaugerhs.bet_rhsU[i])\n",
    "        lhs_names.append(   \"lambdaU\"+str(i))\n",
    "        rhs_exprs.append(rhs.lambda_rhsU[i])\n",
    "        lhs_names.append(        \"vetU\"+str(i))\n",
    "        rhs_exprs.append(gaugerhs.vet_rhsU[i])\n",
    "        for j in range(i,3):\n",
    "            lhs_names.append(   \"aDD\"+str(i)+str(j))\n",
    "            rhs_exprs.append(rhs.a_rhsDD[i][j])\n",
    "            lhs_names.append(   \"hDD\"+str(i)+str(j))\n",
    "            rhs_exprs.append(rhs.h_rhsDD[i][j])\n",
    "\n",
    "    # Sort the lhss list alphabetically, and rhss to match.\n",
    "    #   This ensures the RHSs are evaluated in the same order\n",
    "    #   they're allocated in memory:\n",
    "    lhs_names,rhs_exprs = [list(x) for x in zip(*sorted(zip(lhs_names,rhs_exprs), key=lambda pair: pair[0]))]\n",
    "\n",
    "    # Declare the list of lhrh's\n",
    "    BSSN_evol_rhss = []\n",
    "    for var in range(len(lhs_names)):\n",
    "        BSSN_evol_rhss.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\",lhs_names[var]),rhs=rhs_exprs[var]))\n",
    "\n",
    "    # Set up the C function for the BSSN RHSs\n",
    "    desc=\"Evaluate the BSSN RHSs\"\n",
    "    name=\"rhs_eval\"\n",
    "    outCfunction(\n",
    "        outfile  = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
    "        params   = \"\"\"rfm_struct *restrict rfmstruct,const paramstruct *restrict params,\n",
    "                      const REAL *restrict auxevol_gfs,const REAL *restrict in_gfs,REAL *restrict rhs_gfs\"\"\",\n",
    "        body     = fin.FD_outputC(\"returnstring\",BSSN_evol_rhss, params=\"outCverbose=False,enable_SIMD=True\",\n",
    "                                  upwindcontrolvec=betaU),\n",
    "        loopopts = \"InteriorPoints,enable_SIMD,enable_rfm_precompute\")\n",
    "    end = time.time()\n",
    "    print(\"(BENCH) Finished BSSN_RHS C codegen in \" + str(end - start) + \" seconds.\")\n",
    "\n",
    "def Ricci():\n",
    "    print(\"Generating C code for Ricci tensor in \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" coordinates.\")\n",
    "    start = time.time()\n",
    "    desc=\"Evaluate the Ricci tensor\"\n",
    "    name=\"Ricci_eval\"\n",
    "    outCfunction(\n",
    "        outfile  = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
    "        params   = \"\"\"rfm_struct *restrict rfmstruct,const paramstruct *restrict params,\n",
    "                      const REAL *restrict in_gfs,REAL *restrict auxevol_gfs\"\"\",\n",
    "        body     = fin.FD_outputC(\"returnstring\",\n",
    "                                  [lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD00\"),rhs=Bq.RbarDD[0][0]),\n",
    "                                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD01\"),rhs=Bq.RbarDD[0][1]),\n",
    "                                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD02\"),rhs=Bq.RbarDD[0][2]),\n",
    "                                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD11\"),rhs=Bq.RbarDD[1][1]),\n",
    "                                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD12\"),rhs=Bq.RbarDD[1][2]),\n",
    "                                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD22\"),rhs=Bq.RbarDD[2][2])],\n",
    "                                   params=\"outCverbose=False,enable_SIMD=True\"),\n",
    "        loopopts = \"InteriorPoints,enable_SIMD,enable_rfm_precompute\")\n",
    "    end = time.time()\n",
    "    print(\"(BENCH) Finished Ricci C codegen in \" + str(end - start) + \" seconds.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='hamconstraint'></a>\n",
    "\n",
    "## Step 4.b: Output the Hamiltonian constraint \\[Back to [top](#toc)\\]\n",
    "$$\\label{hamconstraint}$$\n",
    "\n",
    "Next output the C code for evaluating the Hamiltonian constraint [(**Tutorial**)](Tutorial-BSSN_constraints.ipynb). In the absence of numerical error, this constraint should evaluate to zero. However it does not due to numerical (typically truncation and roundoff) error. We will therefore measure the Hamiltonian constraint violation to gauge the accuracy of our simulation, and, ultimately determine whether errors are dominated by numerical finite differencing (truncation) error as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:12:48.094691Z",
     "iopub.status.busy": "2021-06-15T10:12:48.094159Z",
     "iopub.status.idle": "2021-06-15T10:12:48.096356Z",
     "shell.execute_reply": "2021-06-15T10:12:48.095814Z"
    }
   },
   "outputs": [],
   "source": [
    "def Hamiltonian():\n",
    "    start = time.time()\n",
    "    print(\"Generating optimized C code for Hamiltonian constraint. May take a while, depending on CoordSystem.\")\n",
    "    # Set up the C function for the Hamiltonian RHS\n",
    "    desc=\"Evaluate the Hamiltonian constraint\"\n",
    "    name=\"Hamiltonian_constraint\"\n",
    "    outCfunction(\n",
    "        outfile  = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
    "        params   = \"\"\"rfm_struct *restrict rfmstruct,const paramstruct *restrict params,\n",
    "                      REAL *restrict in_gfs, REAL *restrict auxevol_gfs, REAL *restrict aux_gfs\"\"\",\n",
    "        body     = fin.FD_outputC(\"returnstring\",lhrh(lhs=gri.gfaccess(\"aux_gfs\", \"H\"), rhs=bssncon.H),\n",
    "                                  params=\"outCverbose=False\"),\n",
    "        loopopts = \"InteriorPoints,enable_rfm_precompute\")\n",
    "\n",
    "    end = time.time()\n",
    "    print(\"(BENCH) Finished Hamiltonian C codegen in \" + str(end - start) + \" seconds.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='enforce3metric'></a>\n",
    "\n",
    "## Step 4.c: Enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint \\[Back to [top](#toc)\\]\n",
    "$$\\label{enforce3metric}$$\n",
    "\n",
    "Then enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint (Eq. 53 of [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658)), as [documented in the corresponding NRPy+ tutorial notebook](Tutorial-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.ipynb)\n",
    "\n",
    "Applying curvilinear boundary conditions should affect the initial data at the outer boundary, and will in general cause the $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint to be violated there. Thus after we apply these boundary conditions, we must always call the routine for enforcing the $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:12:48.099864Z",
     "iopub.status.busy": "2021-06-15T10:12:48.099331Z",
     "iopub.status.idle": "2021-06-15T10:12:48.101178Z",
     "shell.execute_reply": "2021-06-15T10:12:48.100767Z"
    }
   },
   "outputs": [],
   "source": [
    "def gammadet():\n",
    "    start = time.time()\n",
    "    print(\"Generating optimized C code for gamma constraint. May take a while, depending on CoordSystem.\")\n",
    "\n",
    "    # Set up the C function for the det(gammahat) = det(gammabar)\n",
    "    EGC.output_Enforce_Detgammahat_Constraint_Ccode(Ccodesdir,exprs=enforce_detg_constraint_symb_expressions)\n",
    "    end = time.time()\n",
    "    print(\"(BENCH) Finished gamma constraint C codegen in \" + str(end - start) + \" seconds.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ccodegen'></a>\n",
    "\n",
    "## Step 4.d: Generate C code kernels for BSSN expressions, in parallel if possible \\[Back to [top](#toc)\\]\n",
    "$$\\label{ccodegen}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:12:48.104958Z",
     "iopub.status.busy": "2021-06-15T10:12:48.104565Z",
     "iopub.status.idle": "2021-06-15T10:13:12.087955Z",
     "shell.execute_reply": "2021-06-15T10:13:12.088357Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating C code for Ricci tensor in Spherical coordinates.Generating C code for BSSN RHSs in Spherical coordinates.Generating optimized C code for Hamiltonian constraint. May take a while, depending on CoordSystem.Generating optimized C code for gamma constraint. May take a while, depending on CoordSystem.\n",
      "\n",
      "\n",
      "\n",
      "Output C function enforce_detgammahat_constraint() to file BSSN_ScalarFieldCollapse_Ccodes/enforce_detgammahat_constraint.h\n",
      "(BENCH) Finished gamma constraint C codegen in 0.1816692352294922 seconds.\n",
      "Output C function Ricci_eval() to file BSSN_ScalarFieldCollapse_Ccodes/Ricci_eval.h\n",
      "(BENCH) Finished Ricci C codegen in 21.17385697364807 seconds.\n",
      "Output C function rhs_eval() to file BSSN_ScalarFieldCollapse_Ccodes/rhs_eval.h\n",
      "(BENCH) Finished BSSN_RHS C codegen in 22.861162185668945 seconds.\n",
      "Output C function Hamiltonian_constraint() to file BSSN_ScalarFieldCollapse_Ccodes/Hamiltonian_constraint.h\n",
      "(BENCH) Finished Hamiltonian C codegen in 39.01502728462219 seconds.\n"
     ]
    }
   ],
   "source": [
    "# Step 4.d: C code kernel generation\n",
    "# Step 4.d.i: Create a list of functions we wish to evaluate in parallel\n",
    "funcs = [BSSN_plus_ScalarField_RHSs,Ricci,Hamiltonian,gammadet]\n",
    "\n",
    "try:\n",
    "    if os.name == 'nt':\n",
    "        # It's a mess to get working in Windows, so we don't bother. :/\n",
    "        #  https://medium.com/@grvsinghal/speed-up-your-python-code-using-multiprocessing-on-windows-and-jupyter-or-ipython-2714b49d6fac\n",
    "        raise Exception(\"Parallel codegen currently not available in Windows\")\n",
    "    # Step 4.d.ii: Import the multiprocessing module.\n",
    "    import multiprocess as multiprocessing\n",
    "\n",
    "    # Step 4.d.iii: Define master function for parallelization.\n",
    "    #           Note that lambdifying this doesn't work in Python 3\n",
    "    def master_func(arg):\n",
    "        funcs[arg]()\n",
    "\n",
    "    # Step 4.d.iv: Evaluate list of functions in parallel if possible;\n",
    "    #           otherwise fallback to serial evaluation:\n",
    "    pool = multiprocessing.Pool()\n",
    "    pool.map(master_func,range(len(funcs)))\n",
    "except:\n",
    "    # Steps 4.d.iii-4.d.v, alternate: As fallback, evaluate functions in serial.\n",
    "    for func in funcs:\n",
    "        func()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cparams_rfm_and_domainsize'></a>\n",
    "\n",
    "## Step 4.e: Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h` \\[Back to [top](#toc)\\]\n",
    "$$\\label{cparams_rfm_and_domainsize}$$\n",
    "\n",
    "Based on declared NRPy+ Cparameters, first we generate `declare_Cparameters_struct.h`, `set_Cparameters_default.h`, and `set_Cparameters[-SIMD].h`.\n",
    "\n",
    "Then 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": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:13:12.093427Z",
     "iopub.status.busy": "2021-06-15T10:13:12.093016Z",
     "iopub.status.idle": "2021-06-15T10:13:12.178095Z",
     "shell.execute_reply": "2021-06-15T10:13:12.178546Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 4.e.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n",
    "par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))\n",
    "\n",
    "# Step 4.e.ii: Set free_parameters.h\n",
    "# Output to $Ccodesdir/free_parameters.h reference metric parameters based on generic\n",
    "#    domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale,\n",
    "#    parameters set above.\n",
    "rfm.out_default_free_parameters_for_rfm(os.path.join(Ccodesdir,\"free_parameters.h\"),\n",
    "                                        domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale)\n",
    "\n",
    "# Step 4.e.iii: Generate set_Nxx_dxx_invdx_params__and__xx.h:\n",
    "rfm.set_Nxx_dxx_invdx_params__and__xx_h(Ccodesdir)\n",
    "\n",
    "# Step 4.e.iv: Generate xx_to_Cart.h, which contains xx_to_Cart() for\n",
    "#               (the mapping from xx->Cartesian) for the chosen\n",
    "#               CoordSystem:\n",
    "rfm.xx_to_Cart_h(\"xx_to_Cart\",\"./set_Cparameters.h\",os.path.join(Ccodesdir,\"xx_to_Cart.h\"))\n",
    "\n",
    "# Step 4.e.v: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n",
    "par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bc_functs'></a>\n",
    "\n",
    "# Step 5: Set up boundary condition functions for chosen singular, curvilinear coordinate system \\[Back to [top](#toc)\\]\n",
    "$$\\label{bc_functs}$$\n",
    "\n",
    "Next apply singular, curvilinear coordinate boundary conditions [as documented in the corresponding NRPy+ tutorial notebook](Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:13:12.181581Z",
     "iopub.status.busy": "2021-06-15T10:13:12.181073Z",
     "iopub.status.idle": "2021-06-15T10:13:12.317547Z",
     "shell.execute_reply": "2021-06-15T10:13:12.317011Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"BSSN_ScalarFieldCollapse_Ccodes/boundary_conditions/parity_conditions_symbolic_dot_products.h\"\n",
      "Evolved parity: ( aDD00:4, aDD01:5, aDD02:6, aDD11:7, aDD12:8, aDD22:9,\n",
      "    alpha:0, betU0:1, betU1:2, betU2:3, cf:0, hDD00:4, hDD01:5, hDD02:6,\n",
      "    hDD11:7, hDD12:8, hDD22:9, lambdaU0:1, lambdaU1:2, lambdaU2:3, sf:0,\n",
      "    sfM:0, trK:0, vetU0:1, vetU1:2, vetU2:3 )\n",
      "Auxiliary parity: ( H:0 )\n",
      "AuxEvol parity: ( RbarDD00:4, RbarDD01:5, RbarDD02:6, RbarDD11:7,\n",
      "    RbarDD12:8, RbarDD22:9 )\n",
      "Wrote to file \"BSSN_ScalarFieldCollapse_Ccodes/boundary_conditions/EigenCoord_Cart_to_xx.h\"\n"
     ]
    }
   ],
   "source": [
    "import CurviBoundaryConditions.CurviBoundaryConditions as cbcs\n",
    "cbcs.Set_up_CurviBoundaryConditions(os.path.join(Ccodesdir,\"boundary_conditions/\"),Cparamspath=os.path.join(\"../\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='main_ccode'></a>\n",
    "\n",
    "# Step 6: The main C code: `ScalarFieldCollapse_Playground.c` \\[Back to [top](#toc)\\]\n",
    "$$\\label{main_ccode}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:13:12.320911Z",
     "iopub.status.busy": "2021-06-15T10:13:12.320506Z",
     "iopub.status.idle": "2021-06-15T10:13:12.322151Z",
     "shell.execute_reply": "2021-06-15T10:13:12.322504Z"
    }
   },
   "outputs": [],
   "source": [
    "# Part P0: Define REAL, set the number of ghost cells NGHOSTS (from NRPy+'s FD_CENTDERIVS_ORDER),\n",
    "#          and set the CFL_FACTOR (which can be overwritten at the command line)\n",
    "\n",
    "with open(os.path.join(Ccodesdir,\"ScalarFieldCollapse_Playground_REAL__NGHOSTS__CFL_FACTOR.h\"), \"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "// Part P0.a: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER\n",
    "#define NGHOSTS \"\"\"+str(int(FD_order/2)+1)+\"\"\"\n",
    "// Part P0.b: Set the numerical precision (REAL) to double, ensuring all floating point\n",
    "//            numbers are stored to at least ~16 significant digits\n",
    "#define REAL \"\"\"+REAL+\"\"\"\n",
    "// Part P0.c: Set the number of ghost cells, from NRPy+'s FD_CENTDERIVS_ORDER\n",
    "REAL CFL_FACTOR = \"\"\"+str(CFL_FACTOR)+\"\"\"; // Set the CFL Factor. Can be overwritten at command line.\\n\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:13:12.328132Z",
     "iopub.status.busy": "2021-06-15T10:13:12.327658Z",
     "iopub.status.idle": "2021-06-15T10:13:12.329926Z",
     "shell.execute_reply": "2021-06-15T10:13:12.329574Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing BSSN_ScalarFieldCollapse_Ccodes/ScalarFieldCollapse_Playground.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/ScalarFieldCollapse_Playground.c\n",
    "\n",
    "// Step P0: Define REAL and NGHOSTS; and declare CFL_FACTOR. This header is generated in NRPy+.\n",
    "#include \"ScalarFieldCollapse_Playground_REAL__NGHOSTS__CFL_FACTOR.h\"\n",
    "\n",
    "#include \"rfm_files/rfm_struct__declare.h\"\n",
    "\n",
    "#include \"declare_Cparameters_struct.h\"\n",
    "\n",
    "// All SIMD intrinsics used in SIMD-enabled C code loops are defined here:\n",
    "#include \"SIMD/SIMD_intrinsics.h\"\n",
    "\n",
    "// Step P1: Import needed header files\n",
    "#include \"stdio.h\"\n",
    "#include \"stdlib.h\"\n",
    "#include \"math.h\"\n",
    "#include \"time.h\"\n",
    "#include \"stdint.h\" // Needed for Windows GCC 6.x compatibility\n",
    "#ifndef M_PI\n",
    "#define M_PI 3.141592653589793238462643383279502884L\n",
    "#endif\n",
    "#ifndef M_SQRT1_2\n",
    "#define M_SQRT1_2 0.707106781186547524400844362104849039L\n",
    "#endif\n",
    "#define wavespeed 1.0 // Set CFL-based \"wavespeed\" to 1.0.\n",
    "#define alpha_threshold (2e-3) // Value below which we rule gravitational collapse has happened\n",
    "\n",
    "// Step P2: Declare the IDX4S(gf,i,j,k) macro, which enables us to store 4-dimensions of\n",
    "//           data in a 1D array. In this case, consecutive values of \"i\"\n",
    "//           (all other indices held to a fixed value) are consecutive in memory, where\n",
    "//           consecutive values of \"j\" (fixing all other indices) are separated by\n",
    "//           Nxx_plus_2NGHOSTS0 elements in memory. Similarly, consecutive values of\n",
    "//           \"k\" are separated by Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1 in memory, etc.\n",
    "#define IDX4S(g,i,j,k) \\\n",
    "( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )\n",
    "#define IDX4ptS(g,idx) ( (idx) + (Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2) * (g) )\n",
    "#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) ) ) )\n",
    "#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \\\n",
    "  for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++)\n",
    "#define LOOP_ALL_GFS_GPS(ii) _Pragma(\"omp parallel for\") \\\n",
    "  for(int (ii)=0;(ii)<Nxx_plus_2NGHOSTS_tot*NUM_EVOL_GFS;(ii)++)\n",
    "\n",
    "// Step P3: Set UUGF and VVGF macros, as well as xx_to_Cart()\n",
    "#include \"boundary_conditions/gridfunction_defines.h\"\n",
    "\n",
    "// Step P4: Set xx_to_Cart(const paramstruct *restrict params,\n",
    "//                     REAL *restrict xx[3],\n",
    "//                     const int i0,const int i1,const int i2,\n",
    "//                     REAL xCart[3]),\n",
    "//           which maps xx->Cartesian via\n",
    "//    {xx[0][i0],xx[1][i1],xx[2][i2]}->{xCart[0],xCart[1],xCart[2]}\n",
    "#include \"xx_to_Cart.h\"\n",
    "\n",
    "// Step P5: Defines set_Nxx_dxx_invdx_params__and__xx(const int EigenCoord, const int Nxx[3],\n",
    "//                                       paramstruct *restrict params, REAL *restrict xx[3]),\n",
    "//          which sets params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for\n",
    "//          the chosen Eigen-CoordSystem if EigenCoord==1, or\n",
    "//          CoordSystem if EigenCoord==0.\n",
    "#include \"set_Nxx_dxx_invdx_params__and__xx.h\"\n",
    "\n",
    "// Step P6: Include basic functions needed to impose curvilinear\n",
    "//          parity and boundary conditions.\n",
    "#include \"boundary_conditions/CurviBC_include_Cfunctions.h\"\n",
    "\n",
    "// Step P7: Implement the algorithm for upwinding.\n",
    "//          *NOTE*: This upwinding is backwards from\n",
    "//          usual upwinding algorithms, because the\n",
    "//          upwinding control vector in BSSN (the shift)\n",
    "//          acts like a *negative* velocity.\n",
    "//#define UPWIND_ALG(UpwindVecU) UpwindVecU > 0.0 ? 1.0 : 0.0\n",
    "\n",
    "// Step P8: Include function for enforcing detgammabar constraint.\n",
    "#include \"enforce_detgammahat_constraint.h\"\n",
    "\n",
    "// Step P9: Find the CFL-constrained timestep\n",
    "#include \"find_timestep.h\"\n",
    "\n",
    "// Step P10: Declare initial data input struct:\n",
    "//           stores data from initial data solver,\n",
    "//           so they can be put on the numerical grid.\n",
    "typedef struct __ID_inputs {\n",
    "    int interp_stencil_size;\n",
    "    int numlines_in_file;\n",
    "    REAL *r_arr,*sf_arr,*psi4_arr,*alpha_arr;\n",
    "} ID_inputs;\n",
    "\n",
    "// Part P11: Declare all functions for setting up ScalarField initial data.\n",
    "/* Routines to interpolate the ScalarField solution and convert to ADM & T^{munu}: */\n",
    "#include \"../ScalarField/ScalarField_interp.h\"\n",
    "#include \"ID_scalarfield_ADM_quantities.h\"\n",
    "#include \"ID_scalarfield_spherical.h\"\n",
    "#include \"ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2.h\"\n",
    "#include \"ID_scalarfield.h\"\n",
    "\n",
    "/* Next perform the basis conversion and compute all needed BSSN quantities */\n",
    "#include \"ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs.h\"\n",
    "#include \"ID_BSSN__ALL_BUT_LAMBDAs.h\"\n",
    "#include \"ID_BSSN_lambdas.h\"\n",
    "\n",
    "// Step P12: Set the generic driver function for setting up BSSN initial data\n",
    "void initial_data(const paramstruct *restrict params,const bc_struct *restrict bcstruct,\n",
    "                  const rfm_struct *restrict rfmstruct,\n",
    "                  REAL *restrict xx[3], REAL *restrict auxevol_gfs, REAL *restrict in_gfs) {\n",
    "#include \"set_Cparameters.h\"\n",
    "\n",
    "    // Step 1: Set up ScalarField initial data\n",
    "    // Step 1.a: Read ScalarField initial data from data file\n",
    "    // Open the data file:\n",
    "    char filename[100];\n",
    "    sprintf(filename,\"./SFID.txt\");\n",
    "    FILE *fp = fopen(filename, \"r\");\n",
    "    if (fp == NULL) {\n",
    "        fprintf(stderr,\"ERROR: could not open file %s\\n\",filename);\n",
    "        exit(1);\n",
    "    }\n",
    "    // Count the number of lines in the data file:\n",
    "    int numlines_in_file = count_num_lines_in_file(fp);\n",
    "    // Allocate space for all data arrays:\n",
    "    REAL *r_arr     = (REAL *)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    REAL *sf_arr    = (REAL *)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    REAL *psi4_arr  = (REAL *)malloc(sizeof(REAL)*numlines_in_file);\n",
    "    REAL *alpha_arr = (REAL *)malloc(sizeof(REAL)*numlines_in_file);\n",
    "\n",
    "    // Read from the data file, filling in arrays\n",
    "    // read_datafile__set_arrays() may be found in ScalarField/ScalarField_interp.h\n",
    "    if(read_datafile__set_arrays(fp,r_arr,sf_arr,psi4_arr,alpha_arr) == 1) {\n",
    "        fprintf(stderr,\"ERROR WHEN READING FILE %s!\\n\",filename);\n",
    "        exit(1);\n",
    "    }\n",
    "    fclose(fp);\n",
    "\n",
    "    const int interp_stencil_size = 12;\n",
    "    ID_inputs SF_in;\n",
    "    SF_in.interp_stencil_size = interp_stencil_size;\n",
    "    SF_in.numlines_in_file    = numlines_in_file;\n",
    "    SF_in.r_arr               = r_arr;\n",
    "    SF_in.sf_arr              = sf_arr;\n",
    "    SF_in.psi4_arr            = psi4_arr;\n",
    "    SF_in.alpha_arr           = alpha_arr;\n",
    "\n",
    "    // Step 1.b: Interpolate data from data file to set BSSN gridfunctions\n",
    "    ID_scalarfield(params,xx,SF_in, in_gfs);\n",
    "    ID_BSSN__ALL_BUT_LAMBDAs(params,xx,SF_in, in_gfs);\n",
    "    apply_bcs_curvilinear(params, bcstruct, NUM_EVOL_GFS, evol_gf_parity, in_gfs);\n",
    "    enforce_detgammahat_constraint(rfmstruct, params,                   in_gfs);\n",
    "    ID_BSSN_lambdas(params, xx, in_gfs);\n",
    "    apply_bcs_curvilinear(params, bcstruct, NUM_EVOL_GFS, evol_gf_parity, in_gfs);\n",
    "    enforce_detgammahat_constraint(rfmstruct, params,                   in_gfs);\n",
    "\n",
    "    free(r_arr);\n",
    "    free(sf_arr);\n",
    "    free(psi4_arr);\n",
    "    free(alpha_arr);\n",
    "}\n",
    "\n",
    "// Step P11: Declare function for evaluating Hamiltonian constraint (diagnostic)\n",
    "#include \"Hamiltonian_constraint.h\"\n",
    "\n",
    "// Step P12: Declare rhs_eval function, which evaluates BSSN RHSs\n",
    "#include \"rhs_eval.h\"\n",
    "\n",
    "// Step P13: Declare Ricci_eval function, which evaluates Ricci tensor\n",
    "#include \"Ricci_eval.h\"\n",
    "\n",
    "//#include \"NRPyCritCol_regridding.h\"\n",
    "\n",
    "REAL rho_max = 0.0;\n",
    "\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 initial data to an exact solution\n",
    "// Step 2: Start the timer, for keeping track of how fast the simulation is progressing.\n",
    "// Step 3: Integrate the initial data forward in time using the chosen RK-like Method of\n",
    "//         Lines timestepping algorithm, and output periodic simulation diagnostics\n",
    "// Step 3.a: Output 2D data file periodically, for visualization\n",
    "// Step 3.b: Step forward one timestep (t -> t+dt) in time using\n",
    "//           chosen RK-like MoL timestepping algorithm\n",
    "// Step 3.c: If t=t_final, output conformal factor & Hamiltonian\n",
    "//           constraint violation to 1D data file\n",
    "// Step 3.d: Progress indicator printing to stderr\n",
    "// Step 4: Free all allocated memory\n",
    "int main(int argc, const char *argv[]) {\n",
    "    paramstruct params;\n",
    "#include \"set_Cparameters_default.h\"\n",
    "\n",
    "    // Step 0a: Read command-line input, error out if nonconformant\n",
    "    if((argc != 4 && argc != 5) || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < 2 || atoi(argv[3]) < 2 /* FIXME; allow for axisymmetric sims */) {\n",
    "        fprintf(stderr,\"Error: Expected three command-line arguments: ./ScalarFieldCollapse_Playground Nx0 Nx1 Nx2,\\n\");\n",
    "        fprintf(stderr,\"where Nx[0,1,2] is the number of grid points in the 0, 1, and 2 directions.\\n\");\n",
    "        fprintf(stderr,\"Nx[] MUST BE larger than NGHOSTS (= %d)\\n\",NGHOSTS);\n",
    "        exit(1);\n",
    "    }\n",
    "    if(argc == 5) {\n",
    "        CFL_FACTOR = strtod(argv[4],NULL);\n",
    "        if(CFL_FACTOR > 0.5 && atoi(argv[3])!=2) {\n",
    "            fprintf(stderr,\"WARNING: CFL_FACTOR was set to %e, which is > 0.5.\\n\",CFL_FACTOR);\n",
    "            fprintf(stderr,\"         This will generally only be stable if the simulation is purely axisymmetric\\n\");\n",
    "            fprintf(stderr,\"         However, Nx2 was set to %d>2, which implies a non-axisymmetric simulation\\n\",atoi(argv[3]));\n",
    "        }\n",
    "    }\n",
    "    // Step 0b: Set up numerical grid structure, first in space...\n",
    "    const int Nxx[3] = { atoi(argv[1]), atoi(argv[2]), atoi(argv[3]) };\n",
    "    if(Nxx[0]%2 != 0 || Nxx[1]%2 != 0 || Nxx[2]%2 != 0) {\n",
    "        fprintf(stderr,\"Error: Cannot guarantee a proper cell-centered grid if number of grid cells not set to even number.\\n\");\n",
    "        fprintf(stderr,\"       For example, in case of angular directions, proper symmetry zones will not exist.\\n\");\n",
    "        exit(1);\n",
    "    }\n",
    "\n",
    "    // Step 0c: Set free parameters, overwriting Cparameters defaults\n",
    "    //          by hand or with command-line input, as desired.\n",
    "#include \"free_parameters.h\"\n",
    "\n",
    "   // Step 0d: Uniform coordinate grids are stored to *xx[3]\n",
    "    REAL *xx[3];\n",
    "    // Step 0d.i: Set bcstruct\n",
    "    bc_struct bcstruct;\n",
    "    {\n",
    "        int EigenCoord = 1;\n",
    "        // Step 0d.ii: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
    "        //             params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
    "        //             chosen Eigen-CoordSystem.\n",
    "        set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &params, xx);\n",
    "        // Step 0d.iii: Set Nxx_plus_2NGHOSTS_tot\n",
    "#include \"set_Cparameters-nopointer.h\"\n",
    "        const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
    "        // Step 0e: Find ghostzone mappings; set up bcstruct\n",
    "#include \"boundary_conditions/driver_bcstruct.h\"\n",
    "        // Step 0e.i: Free allocated space for xx[][] array\n",
    "        for(int i=0;i<3;i++) free(xx[i]);\n",
    "    }\n",
    "\n",
    "    // Step 0f: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
    "    //          params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
    "    //          chosen (non-Eigen) CoordSystem.\n",
    "    int EigenCoord = 0;\n",
    "    set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &params, xx);\n",
    "\n",
    "    // Step 0g: Set all C parameters \"blah\" for params.blah, including\n",
    "    //          Nxx_plus_2NGHOSTS0 = params.Nxx_plus_2NGHOSTS0, etc.\n",
    "#include \"set_Cparameters-nopointer.h\"\n",
    "    const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
    "\n",
    "    // Step 0h: Time coordinate parameters\n",
    "    REAL t_final = 16.0; /* 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",
    "\n",
    "    // Step 0i: Set timestep based on smallest proper distance between gridpoints and CFL factor\n",
    "    REAL dt = find_timestep(&params, xx);\n",
    "    //fprintf(stderr,\"# Timestep set to = %e\\n\",(double)dt);\n",
    "    int N_final = (int)(t_final / dt + 0.5); // The number of points in time.\n",
    "                                             // Add 0.5 to account for C rounding down\n",
    "                                             // typecasts to integers.\n",
    "    int output_every_N = 20;//(int)((REAL)N_final/800.0);\n",
    "    if(output_every_N == 0) output_every_N = 1;\n",
    "\n",
    "    // Step 0j: Error out if the number of auxiliary gridfunctions outnumber evolved gridfunctions.\n",
    "    //              This is a limitation of the RK method. You are always welcome to declare & allocate\n",
    "    //              additional gridfunctions by hand.\n",
    "    if(NUM_AUX_GFS > NUM_EVOL_GFS) {\n",
    "        fprintf(stderr,\"Error: NUM_AUX_GFS > NUM_EVOL_GFS. Either reduce the number of auxiliary gridfunctions,\\n\");\n",
    "        fprintf(stderr,\"       or allocate (malloc) by hand storage for *diagnostic_output_gfs. \\n\");\n",
    "        exit(1);\n",
    "    }\n",
    "\n",
    "    // Step 0k: Allocate memory for gridfunctions\n",
    "#include \"MoLtimestepping/RK_Allocate_Memory.h\"\n",
    "    REAL *restrict auxevol_gfs = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS_tot);\n",
    "\n",
    "    // Step 0l: Set up precomputed reference metric arrays\n",
    "    // Step 0l.i: Allocate space for precomputed reference metric arrays.\n",
    "#include \"rfm_files/rfm_struct__malloc.h\"\n",
    "\n",
    "    // Step 0l.ii: Define precomputed reference metric arrays.\n",
    "    {\n",
    "    #include \"set_Cparameters-nopointer.h\"\n",
    "    #include \"rfm_files/rfm_struct__define.h\"\n",
    "    }\n",
    "\n",
    "    // Step 1: Set up initial data to an exact solution\n",
    "    initial_data(&params,&bcstruct, &rfmstruct, xx, auxevol_gfs, y_n_gfs);\n",
    "\n",
    "    // Step 1b: Apply boundary conditions, as initial data\n",
    "    //          are sometimes ill-defined in ghost zones.\n",
    "    //          E.g., spherical initial data might not be\n",
    "    //          properly defined at points where r=-1.\n",
    "    apply_bcs_curvilinear(&params, &bcstruct, NUM_EVOL_GFS,evol_gf_parity, y_n_gfs);\n",
    "    enforce_detgammahat_constraint(&rfmstruct, &params, y_n_gfs);\n",
    "\n",
    "    // Step 2: Start the timer, for keeping track of how fast the simulation is progressing.\n",
    "#ifdef __linux__ // Use high-precision timer in Linux.\n",
    "    struct timespec start, end;\n",
    "    clock_gettime(CLOCK_REALTIME, &start);\n",
    "#else     // Resort to low-resolution, standards-compliant timer in non-Linux OSs\n",
    "    // http://www.cplusplus.com/reference/ctime/time/\n",
    "    time_t start_timer,end_timer;\n",
    "    time(&start_timer); // Resolution of one second...\n",
    "#endif\n",
    "\n",
    "    // Step 3: Integrate the initial data forward in time using the chosen RK-like Method of\n",
    "    //         Lines timestepping algorithm, and output periodic simulation diagnostics\n",
    "    for(int n=0;n<=N_final;n++) { // Main loop to progress forward in time.\n",
    "\n",
    "        // Step 3.a: Output 2D data file periodically, for visualization\n",
    "        if(n%output_every_N == 0) {\n",
    "            // Evaluate Hamiltonian constraint violation\n",
    "            Hamiltonian_constraint(&rfmstruct, &params, y_n_gfs,auxevol_gfs, diagnostic_output_gfs);\n",
    "\n",
    "            char filename[100];\n",
    "            sprintf(filename,\"out%d-%08d.txt\",Nxx[0],n);\n",
    "            const int i1mid=Nxx_plus_2NGHOSTS1/2;\n",
    "            const int i2mid=Nxx_plus_2NGHOSTS2/2;\n",
    "            FILE *fp = fopen(filename, \"w\");\n",
    "            for( int i0=NGHOSTS;i0<Nxx_plus_2NGHOSTS0-NGHOSTS;i0++) {\n",
    "                const int idx  = IDX3S(i0,i1mid,i2mid);\n",
    "                const REAL xx0 = xx[0][i0];\n",
    "                REAL xCart[3];\n",
    "                xx_to_Cart(&params,xx,i0,i1mid,i2mid,xCart);\n",
    "                const REAL rr = sqrt( xCart[0]*xCart[0] + xCart[1]*xCart[1] + xCart[2]*xCart[2] );\n",
    "                fprintf(fp,\"%e %e %e %e %e %e %e\\n\",xx0,rr,\n",
    "                        y_n_gfs[IDX4ptS(SFGF,idx)],y_n_gfs[IDX4ptS(SFMGF,idx)],\n",
    "                        y_n_gfs[IDX4ptS(ALPHAGF,idx)],y_n_gfs[IDX4ptS(CFGF,idx)],\n",
    "                        log10(fabs(diagnostic_output_gfs[IDX4ptS(HGF,idx)])));\n",
    "            }\n",
    "            fclose(fp);\n",
    "        }\n",
    "\n",
    "        // Step 3.b: Step forward one timestep (t -> t+dt) in time using\n",
    "        //           chosen RK-like MoL timestepping algorithm\n",
    "#include \"MoLtimestepping/RK_MoL.h\"\n",
    "\n",
    "        // Step 3.c: If t=t_final, output conformal factor & Hamiltonian\n",
    "        //           constraint violation to 2D data file\n",
    "        if(n==N_final-1) {\n",
    "            // Evaluate Hamiltonian constraint violation\n",
    "            Hamiltonian_constraint(&rfmstruct, &params, y_n_gfs,auxevol_gfs, diagnostic_output_gfs);\n",
    "            char filename[100];\n",
    "            sprintf(filename,\"out%d.txt\",Nxx[0]);\n",
    "            FILE *out1D = fopen(filename, \"w\");\n",
    "            const int i1mid=Nxx_plus_2NGHOSTS1/2;\n",
    "            const int i2mid=Nxx_plus_2NGHOSTS2/2;\n",
    "            for(int i0=NGHOSTS;i0<Nxx_plus_2NGHOSTS0-NGHOSTS;i0++) {\n",
    "                REAL xCart[3];\n",
    "                xx_to_Cart(&params,xx,i0,i1mid,i2mid,xCart);\n",
    "                const REAL rr = sqrt( xCart[0]*xCart[0] + xCart[1]*xCart[1] + xCart[2]*xCart[2] );\n",
    "                int idx = IDX3S(i0,i1mid,i2mid);\n",
    "                fprintf(out1D,\"%e %e\\n\",rr,log10(fabs(diagnostic_output_gfs[IDX4ptS(HGF,idx)])));\n",
    "            }\n",
    "            fclose(out1D);\n",
    "        }\n",
    "\n",
    "        // Step 3.d: Progress indicator printing to stderr\n",
    "\n",
    "        // Step 3.d.i: Measure average time per iteration\n",
    "#ifdef __linux__ // Use high-precision timer in Linux.\n",
    "        clock_gettime(CLOCK_REALTIME, &end);\n",
    "        const long long unsigned int time_in_ns = 1000000000L * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;\n",
    "#else     // Resort to low-resolution, standards-compliant timer in non-Linux OSs\n",
    "        time(&end_timer); // Resolution of one second...\n",
    "        REAL time_in_ns = difftime(end_timer,start_timer)*1.0e9+0.5; // Round up to avoid divide-by-zero.\n",
    "#endif\n",
    "        const REAL s_per_iteration_avg = ((REAL)time_in_ns / (REAL)n) / 1.0e9;\n",
    "\n",
    "        const int iterations_remaining = N_final - n;\n",
    "        const REAL time_remaining_in_mins = s_per_iteration_avg * (REAL)iterations_remaining / 60.0;\n",
    "\n",
    "        const REAL num_RHS_pt_evals = (REAL)(Nxx[0]*Nxx[1]*Nxx[2]) * 4.0 * (REAL)n; // 4 RHS evals per gridpoint for RK4\n",
    "        const REAL RHS_pt_evals_per_sec = num_RHS_pt_evals / ((REAL)time_in_ns / 1.0e9);\n",
    "\n",
    "\n",
    "        // Step 3.d.ii: Output simulation progress to stderr\n",
    "        if(n%10 == 0) {\n",
    "            fprintf(stderr,\"%c[2K\", 27); // Clear the line\n",
    "            fprintf(stderr,\"It: %d t=%.2f dt=%.2e | %.1f%%; ETA %.0f s | t/h %.2f | gp/s %.2e\\r\",  // \\r is carriage return, move cursor to the beginning of the line\n",
    "                   n, n*dt, (double)dt, (double)(100.0 * (REAL)n / (REAL)N_final),\n",
    "                   (double)time_remaining_in_mins*60, (double)(dt * 3600.0 / s_per_iteration_avg), (double)RHS_pt_evals_per_sec);\n",
    "            fflush(stderr); // Flush the stderr buffer\n",
    "        } // End progress indicator if(n % 10 == 0)\n",
    "    } // End main loop to progress forward in time.\n",
    "    fprintf(stderr,\"\\n\"); // Clear the final line of output from progress indicator.\n",
    "\n",
    "    // Step 4: Free all allocated memory\n",
    "#include \"rfm_files/rfm_struct__freemem.h\"\n",
    "#include \"boundary_conditions/bcstruct_freemem.h\"\n",
    "#include \"MoLtimestepping/RK_Free_Memory.h\"\n",
    "    free(auxevol_gfs);\n",
    "    for(int i=0;i<3;i++) free(xx[i]);\n",
    "\n",
    "    return 0;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:13:12.335088Z",
     "iopub.status.busy": "2021-06-15T10:13:12.332699Z",
     "iopub.status.idle": "2021-06-15T10:13:18.809062Z",
     "shell.execute_reply": "2021-06-15T10:13:18.809526Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Now compiling, should take ~10 seconds...\n",
      "\n",
      "Compiling executable...\n",
      "(EXEC): Executing `gcc -std=gnu99 -Ofast -fopenmp -march=native -funroll-loops BSSN_ScalarFieldCollapse_Ccodes/ScalarFieldCollapse_Playground.c -o BSSN_ScalarFieldCollapse_Ccodes/output/ScalarFieldCollapse_Playground -lm`...\n",
      "(BENCH): Finished executing in 6.512221097946167 seconds.\n",
      "Finished compilation.\n",
      "(BENCH) Finished in 6.542217969894409 seconds.\n",
      "\n",
      "/Users/werneck/Codes/nrpytutorial/BSSN_ScalarFieldCollapse_Ccodes/output\n",
      "Now running, should take ~20 seconds...\n",
      "\n",
      "(EXEC): Executing `./ScalarFieldCollapse_Playground 640 2 2 0.5`...\n",
      "\u001b[2KIt: 810 t=15.90 dt=1.96e-02 | 99.4%; ETA 0 s | t/h 3817.04 | gp/s 5.53e+0505e+14\n",
      "(BENCH): Finished executing in 16.081330060958862 seconds.\n",
      "(BENCH) Finished in 16.09248399734497 seconds.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import time\n",
    "import cmdline_helper as cmd\n",
    "\n",
    "print(\"Now compiling, should take ~10 seconds...\\n\")\n",
    "start = time.time()\n",
    "cmd.C_compile(os.path.join(Ccodesdir,\"ScalarFieldCollapse_Playground.c\"),\n",
    "              os.path.join(outdir,\"ScalarFieldCollapse_Playground\"),compile_mode=\"optimized\")\n",
    "end = time.time()\n",
    "print(\"(BENCH) Finished in \"+str(end-start)+\" seconds.\\n\")\n",
    "\n",
    "# Change to output directory\n",
    "os.chdir(outdir)\n",
    "# Clean up existing output files\n",
    "cmd.delete_existing_files(\"out*.txt\")\n",
    "cmd.delete_existing_files(\"out*.png\")\n",
    "# Run executable\n",
    "\n",
    "print(os.getcwd())\n",
    "print(\"Now running, should take ~20 seconds...\\n\")\n",
    "start = time.time()\n",
    "cmd.Execute(\"ScalarFieldCollapse_Playground\", \"640 2 2 \"+str(CFL_FACTOR),\"out640.txt\")\n",
    "end = time.time()\n",
    "print(\"(BENCH) Finished in \"+str(end-start)+\" seconds.\\n\")\n",
    "\n",
    "# Return to root directory\n",
    "os.chdir(os.path.join(\"../../\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='visualization'></a>\n",
    "\n",
    "# Step 7: Visualization \\[Back to [top](#toc)\\]\n",
    "$$\\label{visualization}$$\n",
    "\n",
    "<a id='install_download'></a>\n",
    "\n",
    "## Step 7.a: Install `scipy` and download `ffmpeg` if they are not yet installed/downloaded \\[Back to [top](#toc)\\]\n",
    "$$\\label{install_download}$$ \n",
    "\n",
    "Note that if you are not running this within `mybinder`, but on a Windows system, `ffmpeg` must be installed using a separate package (on [this site](http://ffmpeg.org/)), or if running Jupyter within Anaconda, use the command: `conda install -c conda-forge ffmpeg`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install scipy > /dev/null\n",
    "\n",
    "check_for_ffmpeg = !which ffmpeg >/dev/null && echo $?\n",
    "if check_for_ffmpeg != ['0']:\n",
    "    print(\"Couldn't find ffmpeg, so I'll download it.\")\n",
    "    # Courtesy https://johnvansickle.com/ffmpeg/\n",
    "    !wget http://astro.phys.wvu.edu/zetienne/ffmpeg-static-amd64-johnvansickle.tar.xz\n",
    "    !tar Jxf ffmpeg-static-amd64-johnvansickle.tar.xz\n",
    "    print(\"Copying ffmpeg to ~/.local/bin/. Assumes ~/.local/bin is in the PATH.\")\n",
    "    !mkdir ~/.local/bin/\n",
    "    !cp ffmpeg-static-amd64-johnvansickle/ffmpeg ~/.local/bin/\n",
    "    print(\"If this doesn't work, then install ffmpeg yourself. It should work fine on mybinder.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='movie_dynamics'></a>\n",
    "\n",
    "## Step 7.b: Dynamics of the solution \\[Back to [top](#toc)\\]\n",
    "$$\\label{movie_dynamics}$$\n",
    "\n",
    "<a id='genimages'></a>\n",
    "\n",
    "### Step 7.b.i: Generate images for visualization animation \\[Back to [top](#toc)\\]\n",
    "$$\\label{genimages}$$ \n",
    "\n",
    "Here we loop through the data files output by the executable compiled and run in [the previous step](#mainc), generating a [png](https://en.wikipedia.org/wiki/Portable_Network_Graphics) image for each data file.\n",
    "\n",
    "**Special thanks to Terrence Pierre Jacques. His work with the first versions of these scripts greatly contributed to the scripts as they exist below.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:13:18.821583Z",
     "iopub.status.busy": "2021-06-15T10:13:18.821066Z",
     "iopub.status.idle": "2021-06-15T10:13:31.949721Z",
     "shell.execute_reply": "2021-06-15T10:13:31.949223Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[2KProcessing file BSSN_ScalarFieldCollapse_Ccodes/output/out640-00000800.txt\r"
     ]
    }
   ],
   "source": [
    "## VISUALIZATION ANIMATION, PART 1: Generate PNGs, one per frame of movie ##\n",
    "import numpy as np\n",
    "from scipy.interpolate import griddata\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.pyplot import savefig\n",
    "\n",
    "import glob\n",
    "import sys\n",
    "from matplotlib import animation\n",
    "\n",
    "globby = glob.glob(os.path.join(outdir,'out640-00*.txt'))\n",
    "file_list = []\n",
    "for x in sorted(globby):\n",
    "    file_list.append(x)\n",
    "\n",
    "for filename in file_list:\n",
    "    fig = plt.figure(figsize=(8,6))\n",
    "    x,r,sf,sfM,alpha,cf,logH = np.loadtxt(filename).T #Transposed for easier unpacking\n",
    "\n",
    "    ax  = fig.add_subplot(321)\n",
    "    ax2 = fig.add_subplot(322)\n",
    "    ax3 = fig.add_subplot(323)\n",
    "    ax4 = fig.add_subplot(324)\n",
    "    ax5 = fig.add_subplot(325)\n",
    "\n",
    "    ax.set_title(\"Scalar field\")\n",
    "    ax.set_ylabel(r\"$\\varphi(t,r)$\")\n",
    "    ax.set_xlim(0,20)\n",
    "    ax.set_ylim(-0.6,0.6)\n",
    "    ax.plot(r,sf,'k-')\n",
    "    ax.grid()\n",
    "\n",
    "    ax2.set_title(\"Scalar field conjugate momentum\")\n",
    "    ax2.set_ylabel(r\"$\\Pi(t,r)$\")\n",
    "    ax2.set_xlim(0,20)\n",
    "    ax2.set_ylim(-1,1)\n",
    "    ax2.plot(r,sfM,'b-')\n",
    "    ax2.grid()\n",
    "\n",
    "    ax3.set_title(\"Lapse function\")\n",
    "    ax3.set_ylabel(r\"$\\alpha(t,r)$\")\n",
    "    ax3.set_xlim(0,20)\n",
    "    ax3.set_ylim(0,1.02)\n",
    "    ax3.plot(r,alpha,'r-')\n",
    "    ax3.grid()\n",
    "\n",
    "    ax4.set_title(\"Conformal factor\")\n",
    "    ax4.set_xlabel(r\"$r$\")\n",
    "    ax4.set_ylabel(r\"$W(t,r)$\")\n",
    "    ax4.set_xlim(0,20)\n",
    "    ax4.set_ylim(0,1.02)\n",
    "    ax4.plot(r,cf,'g-',label=(\"$p = 0.043149493$\"))\n",
    "    ax4.grid()\n",
    "\n",
    "    ax5.set_title(\"Hamiltonian constraint violation\")\n",
    "    ax5.set_xlabel(r\"$r$\")\n",
    "    ax5.set_ylabel(r\"$\\mathcal{H}(t,r)$\")\n",
    "    ax5.set_xlim(0,20)\n",
    "    ax5.set_ylim(-16,0)\n",
    "    ax5.plot(r,logH,'m-')\n",
    "    ax5.grid()\n",
    "\n",
    "    plt.tight_layout()\n",
    "    savefig(filename+\".png\",dpi=150)\n",
    "    plt.close(fig)\n",
    "    sys.stdout.write(\"%c[2K\" % 27)\n",
    "    sys.stdout.write(\"Processing file \"+filename+\"\\r\")\n",
    "    sys.stdout.flush()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='genvideo'></a>\n",
    "\n",
    "### Step 7.b.ii: Generate visualization animation \\[Back to [top](#toc)\\]\n",
    "$$\\label{genvideo}$$ \n",
    "\n",
    "In the following step, [ffmpeg](http://ffmpeg.org) is used to generate an [mp4](https://en.wikipedia.org/wiki/MPEG-4) video file, which can be played directly from this Jupyter notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:13:31.963581Z",
     "iopub.status.busy": "2021-06-15T10:13:31.963042Z",
     "iopub.status.idle": "2021-06-15T10:13:36.778710Z",
     "shell.execute_reply": "2021-06-15T10:13:36.778133Z"
    }
   },
   "outputs": [],
   "source": [
    "## VISUALIZATION ANIMATION, PART 2: Combine PNGs to generate movie ##\n",
    "\n",
    "# https://stackoverflow.com/questions/14908576/how-to-remove-frame-from-matplotlib-pyplot-figure-vs-matplotlib-figure-frame\n",
    "# https://stackoverflow.com/questions/23176161/animating-pngs-in-matplotlib-using-artistanimation\n",
    "from IPython.display import HTML\n",
    "import matplotlib.image as mgimg\n",
    "\n",
    "fig = plt.figure(frameon=False)\n",
    "ax = fig.add_axes([0, 0, 1, 1])\n",
    "ax.axis('off')\n",
    "\n",
    "myimages = []\n",
    "\n",
    "for i in range(len(file_list)):\n",
    "    img = mgimg.imread(file_list[i]+\".png\")\n",
    "    imgplot = plt.imshow(img)\n",
    "    myimages.append([imgplot])\n",
    "\n",
    "ani = animation.ArtistAnimation(fig, myimages, interval=100,  repeat_delay=1000)\n",
    "plt.close()\n",
    "ani.save(os.path.join(outdir,'ScalarField_Collapse.mp4'), fps=5, dpi=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "## VISUALIZATION ANIMATION, PART 3: Display movie as embedded HTML5 (see next cell) ##\n",
    "\n",
    "# https://stackoverflow.com/questions/18019477/how-can-i-play-a-local-video-in-my-ipython-notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"800\" height=\"600\" controls>\n",
       "  <source src=\"BSSN_ScalarFieldCollapse_Ccodes/output/ScalarField_Collapse.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Embed video based on suggestion:\n",
    "#  https://stackoverflow.com/questions/39900173/jupyter-notebook-html-cell-magic-with-python-variable\n",
    "HTML(\"\"\"\n",
    "<video width=\"800\" height=\"600\" controls>\n",
    "  <source src=\\\"\"\"\"+os.path.join(outdir,\"ScalarField_Collapse.mp4\")+\"\"\"\\\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convergence'></a>\n",
    "\n",
    "## Step 7.c: Convergence of constraint violation \\[Back to [top](#toc)\\]\n",
    "$$\\label{convergence}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:13:36.796179Z",
     "iopub.status.busy": "2021-06-15T10:13:36.795771Z",
     "iopub.status.idle": "2021-06-15T10:13:37.838994Z",
     "shell.execute_reply": "2021-06-15T10:13:37.838747Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(EXEC): Executing `./ScalarFieldCollapse_Playground 320 2 2 0.5`...\n",
      "\u001b[2KIt: 400 t=15.71 dt=3.93e-02 | 98.3%; ETA 0 s | t/h 14137.17 | gp/s 5.12e+05.22e+14\n",
      "(BENCH): Finished executing in 4.4767279624938965 seconds.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAJYCAYAAAA6xSjbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAABcSAAAXEgFnn9JSAADdYElEQVR4nOzdd3hUVf7H8fdJJQkQCBB67yBFQEUUQVTEXtG1rGJd69p1f2tdyxb7WlbXir2ga8GGFVRAEZDeey8hkISE9PP7497JTCYzqZOZlM/ree4zuffc8r0zdybznXPuOcZai4iIiIiIiDQ+UZEOQERERERERCJDCaGIiIiIiEgjpYRQRERERESkkVJCKCIiIiIi0kgpIRQREREREWmklBCKiIiIiIg0UkoIRUREREREGiklhCIiIiIiIo2UEkIREREREZFGSgmhiIiIiIhII6WEUEREREREpJFSQigiIiIiItJIKSEUERERERFppJQQioiIiIiINFJKCEVERERERBopJYQiIiIiIiKNlBJCERERERGRRkoJoYiIiIiISCOlhFBERERERKSRUkIoIiIiIiLSSCkhFBERERERaaSUEIqIiIiIiDRSSghFREREREQaKSWEIiIiIiIijZQSQhERERERkUZKCaGIiIiIiEgjpYRQRERERESkkVJCKCISgDHmPmOMNcZMj3Qs1dUQzkFqj66P8hljurnPjzXGdIt0PLXJ5zzHRjoWEQk/JYQi0mD5fOH1n3KNMVuMMZ8aY84xxpgwxjTUjevGcB1Tap8xZpL7uo6tA7E02mvMGPO8z/t8Qznr3eg+R0PDF52ISN0UE+kARETCZKfP38lAR3c6BZhkjDnDWpsXhjiGAvcCG4Enw3A8CY9JwBj37+mRCwOo/DWWBqwENtV+SLXPGHM0cGUlV78R6ApsABbUTkQiIvWDaghFpFGw1rbzTEAScBDwjVt8AvBgxIITiQBr7TPW2n7W2osiHUtNGWMSgReBQmBuhMMREalXlBCKSKNjrS221i4FTgXWuIv/ZIxRqwmR+ukhoCfwMLA0wrGIiNQrSghFpNGy1uYCU9zZZkC/qmxvjDnYGPO6MWaje1/iXmPMLPf+pPgA61vgVXe2a4B7G++r6jkYY04wxnxjjNlnjNlvjFlojLndGBNbye27GWOeNMYsdbfPMcasMMb82xjTJcg2k3zv0TLGjDbGTDXG7DLGZBtjfjfGXOa3zUlunLvdY/xmjDm3gtiijTGXGmO+N8akGWPyjDFbjTFTyrtXzxgz3fN8GscVxphfjTGZxpgsY8xsY8yF5WwfY4y50t1PmjGmwBizxxiz0hjznu+5eZ4LvM1F7w3wunbzWb+k8w5jTKox5nFjzCr3ObE+6yUaY85zr68F7vOWZ4zZZoz52BhzQpDYK32NmXI6lTHGTHbLJrvzZ7vPR7ob6wJjzA3GmKDfI9zn/hL3+c4yxmS4r8OVblmpY1SXMWYk8GdgFRXU9HvOGae5KMCr/s9ROdu2dd8X643zft9pjHnXGFOlz40A+002xtxjjJnvXqMHjDGrjTHPGWN6lLNdpa4ld92WxphHjDFr3di3u++j4ZWMMcoYc4Ex5gv3vPPda/Jr9zoNeB+2MWaDG+MkY0xTY8z9xpjF7vXQ4DvrEak3rLWaNGnS1CAn4D7AOh91Qde5xrMOMCrAttODbHcTUOyz7T4g32d+IdDeb5sdQIZbXuTO+063Vvf83GkvUOD+PQP4ewXncAGQ67N9LpDjM58JjA+w3SS3fANwuXsuxe5z4BvPP9z1/+Zzzv7rXBUktmTgB5/1Ct3z833OHwmy7XS3/AHgY/fvAp/n3jP9LcC20cDXfuvt83uerM/657qvnee13x/gde3ss75nH5e7ZRY44D7XNsBzbH2e22y/uB4NEH+lrzHKucaByW7ZZOAZn/3t9YvhtSCvQTTwrt85pLv7sMDbvseowXs8Hljm7n+MX+wbAqx/q/s8eOLI8H+OfNbt5hP/STj3IVv3dfC9HjKAIdWMfyCw2WdfJdcC3vfkWUG2rey11A3nvepZP8/nGsnDaSnhKRsb4DgpOJ8n/u8J3/lPgLgA23qOewvO/aqeY3quo27Vfe01adIUuiniAWjSpElTbU1ULiF82OdLTb8A204PsM3JPtt8DHR3l8cBf/T5QjcTiPbbdlKwL6tVPDffL3Hv4yYdQAJOkuv7pSvQORyH86W4APiX+6XRuFNfd5+eL7tdgpxDtnucfwNt3LIUvF/Ii4DbcZK5O4Fkd532wJd4E6jkAPF94PPl8Xog0V3eDnjZ59zLJJR4E8J094vrxUCCW9YJ+NQnvt5+216I94v1ZUBTd7kBUoEzgCnlHPO+Cl43T9xZwApgHBDllvXxWe804BHgCM+5+zx39+BNQE8NcIxKXWNULiFMd1+Dm4DmblkrnPv1POcyLsD2f/Epfwxo5S5vDvwf3gSxpgnhQ+4+XgwQe9Dzx5uoTCpnnW4+55AO/AyMcMtigGOBbW75j9WIvRmwzt1+C3Ciz7UwBJiNNykcUp1rCScx/83nHCYCMW7ZAOBHSif5Y/2OEe1zbf+O89nneS8mARfhTZSfKOd5zgK2A6cDsT7vxcSqPm+aNGkK/RTxADRp0qSptiYqSAjdL6db3XX2eL5M+W07PcB2yzxfAvFL+NzyU3y+YJ3tVzapoi+rlTy3pZ74fOP2Kf+TTwzT/cqicJrXWeDKco7xibvOk0HOodQXcZ/yaJ8vuha4M8hzv98tv9Cv7DCfbQPGhzdh3A008Sub7rP90QG2jfd53e/0K/uPu/y/VXw9PMe8r4L1PHFlAJ1q8Prf6u7n2wBllbrGKrjGJ/vEOinI9nMDXQM4iYKnBuqlCo5d7YQQOBjnB40dQIsAsQc9f6qeEC7H/VHBbx3f93qVXk/gDne7fOCgAOXNgPXuOp9V51oCzvFZ75gA5Yk491EHSwj/6HP+yUGOMRwnwc8DUoM8z4XAwdW93jVp0lS7k+4hFJFGxxjTwhhzDPA90MFd/G9rbXElth0M9HdnH7TWFvmvY62dCsxxZ88LQciBYhjgE0OguF/ESXoCOQrojTPswEvlHOp19/H4ctb5p/8C9zn5zp3NJcDQB9baTJwaEIDBfsWeewu3lBPf3e5ja5zazkBmWmt/CHDsPGBakGPvcx/bBdlnqLxhrd1Sg+0/dx8PN8ZEhyKgIDYDrwUp+9R99H8Ox+Mk/ODU4AXyGE7z5GoxTgdQr+DU1P3ZWruvuvuqpMestQcCLP8SJ6EDGFTFfXqu8w+stUv8C621WTgtGABOMMYkB9lPedfSH9zHmdba7/wLrbU5PscIxHO/7HPW2oxAK1hr5+H8QBUHHB1kP19Za38v5zgiEkFKCEWkUfDrNGIv8C3OL9sAbxL8i6u/Ee5jIc59NcF4hrQYUc461eUbw0+BVnCTxOlBtj/CfUwGthljdgSacJJK8HbA4S/dWrs2SJln3Mdl1trsCtZp6bfcc34/BEvSrbXL8Sa8wZ7jX4MsB6epHzhNXH19gdsU0xjzpdthRgdCb2ZFK7idmPzN7ZRljzGm0OcaXuaulkjZ5y+UfrPW2iBlwZ7DYe7jJmvt+kAbusnOvBrE9Rec8RY/s9a+X4P9VFbAa8laW4hTSw1ln4egjDFxeBPpb8tZ1fM5EoX3efVX3rXkeW98X846AcvcHxpGurP3BfuccD8r+rrrBfusqPB6F5HIURfrItJY+A5Mn4dTO/Y78FagWqRypLqPabb8gew9v9inlrNOdVU1Bn+eBCcWaFuJ4yUEWZ5VzjaFVVjHv0dUz/kFq+H02AJ0JPhzXOVjW2t/NsbcgdNb5QR3whizBeeL++tVvF6C2VVeoTHmcJzktIXP4v14O/2JxqkdBaeJZloIYgqkOq9fG/dxG+Wr6PUNyBgzAKeGeD/O/bLhUJ3noTwpOK8hlP88+L6Hg13n5V1LlXkvBfucSMFpXg2V/9EhMcjycq93EYksJYQi0ihYZ0B6cXi+iP5qrR1Z7pqNkLX2EWPMWzj3X40BRuF0gDEJmGSM+QA431pbUIPDlGlq7OE2h3wHJxlcAPwV+NmtVfOs0xPvGJoBu/yvA4LVLNbUszjNE+8F9hpjmvqVe77bGJ+yvBq+XnVZ0GuphnybIp9grf2qBvuqrRhFJATUZFREpGo8v3S3NgHGGvTRyW/92oohrpz1OgZZvsN9DNa8K9I859ep3LVq8Tm21m6z1j5prT3DWtsWp3mf537Gs4GrQ31MH4fjvDZFwMnW2i99k0FXXf6Bw9OEsqKmtsGuz4p0dx//gVNz5z9d4JZ38Vl2bTWPVVs8Q3BA+de5b1l1rnPPNuU918HK9uCt/ayrnxUiEgJKCEVEqmau+xiDdzDyQI51H3/zW+65J64mtTq+MYwOtII7YPjYINt77udpZ4ypjXsca8pzfkebIAOfu4OBe77I+j/HIWetXWytvQLvc+ffkU0oXlePzu7jbmttsKZ+xwZZHupYqmO++9g12MDjbs1dpQZFryURfY6stfnAInf2mHJW9bzOxXif16ooeS+Vs864QAvdGlVP51inVOPYIlJPKCEUEakCa+0ivB163BWoh0djzIk4QyeA0/TPV6b72KKGMSx3Z+8MkjRdSvCahx/wNjd8ooJaRowxle4sI0TedR874gy6Hcj97mMa5XfKUSUV1PqCMz4heBMKjxq/rj48vTm2NcaUucfTGNMJ+HM524cylur42ieGvwZZ5yaC329WLmttN2utCTbh7RV1o8/yJ/12E+nnCLzX+dnGmIP8C92k+XZ39otgvXxW4D338UhjzNgAx0gAbitn+xfcxxPdz7WgIvA5ISIhooRQRKTq7nAfRwMfGGO6AxhjYo0xF+BNAmfhDFzvy9O9fHNjzDk1iOFO9/Fo4G03ScAY08QYcxXwDN4hFEpxe0a8Cqc52JHAj8aYY4wxJZ1iGGN6GGOuMsb8Rvg67vDENwf40J192hhznTEm0Y2rnTHmRZwBtgHuttbmhvDwHxtjXjHGnGCMaeFZaIxJMcbchbc253O/7Tyv64nGmOo2hfT4GcjGqb163xjTx40h2hhzPN4xD4MJ1TVWLW6vsv9yZ68wxjzsSRaMMc3cTnvuw+ntN1I8z9HZxpja7KW1PM/hjDMYC3zpXnNRAMaYQThDo3TH6QTrrmoe40O8NYsfGmPO8vyIZYzpjzNsRptgG+P0wPwtzrX4kTHmLt9ed40xScaYo40xz+KMPSoi9ZASQhGRKrLWfgbcjPOl/HRgnTFmL06vh2/ijMG2GJjoP06htXYN3jH63jPGZBpjNrjTjVWI4SO8Q2WcC2w2xqTj3C/1HE5Tr+fK2f47nKQqC6c281sg2xiTZozJBda624+g9joHKc9lOMN6xAFPAxnu+W3DW2v4qLX2+RAfNwG4BKeHz73GmAxjTAbO/VQP4Hwx/oCy4yO+hjPmYi9gk9sdv+d1reheyFLcmqBb3dmjgJXGmCyc6+srnOFCLiln+5BcYzX0MM7zBE4N1G739duLM3blW8BUtzyUCX1lvYBzXY9yY9vmeY7CFYB7X+ipOD2AdsK55rLd622RG1secKG1dmE1j1GI8z7fjNNr6AfuMfbhtHQ4HLi4nO2LgLOAz3Deiw8AW933xV6cz4/vcX40SqpOjCISeUoIRUSqwVr7BE6y9CbOl61EnOaEv+A0hzvEWhus2/2zgSeAVTi1A13dqUUVY7gLOBnnC1kmThfxy3HGaDsG74DZwbb/GCeB+RtOArnfjSEPWIiT9JwBPFKVuELBTYqOwUkMp+N88WyK0yHOh8DR1trymrpV1/U4NcBfAKtxEsAEnET0U+Asa+1E//ERrbWrcWprP8XpVKUV3te1yj16u4nuSTjnvt/dx1ac5HgIzg8O5QnJNVZdbiJyDk7yPgfnvRGDc0/b5dbai3xi2ReOmPzi+xHn+f3WPX5bvM9ROONYAgzEqTFdgFNrH4/zg8zzwEBr7QfBtq/kMdbhjNn4OE6NpMFJwj8ARllrP61g+0xr7SnAiThNUDe5MSbiXJNfA/+HdyxCEalnTPDxZkVERERCzxhjcBKLTsBF1to3IhySiEijpRpCERERCbc/4iSDhYSwUyAREak6JYQiIiIScsaYd4wxZxtjWvssa2uM+QvworvodWvt9shEKCIioCajIiIiUgvcjkuS3dkcoMBnHuAn4GRrbSYiIhIxSghFREQk5IwxFwEnAAcDqTidAu3D6TzlXeANd/BzERGJICWEIiIiIiIijZTuIRQREREREWmklBCKiIiIiIg0UkoIRUREREREGiklhCIiIoAxZrQx5nNjzG5jTJExxhpjPo50XBKc+xpZY8zYSMcSasaYye65TQ7xfie5+90Qyv3Wl+OLSFkxkQ5AREQk0owxI4Hvcf4vWmAPUATsjWRcdY2beI0FNlhrJ0cyltpkjOkGTAKw1t4XyVjqCz1nIvWXEkIRERG4Eed/4kzgVGttemTDqbPGAvcCM4DJEY3EsdJ9zAnxfrvhnCfAfSHed2Vtxzm/7RE6flV1o3LPWQbOeW2t5XhEpJKUEIqIiMAg9/FdJYP1h7W2X6RjqC3W2v8D/i/ScYSatfYj4KNIxyEiXrqHUEREBBLdx/0RjUJERCTMlBCKSKUZYzobYx42xiwwxmQYYw4YY9YaYz4xxlxkjGkSYJtoY8ylxpjvjTFpxpg8Y8xWY8yU8jqCMMZMdzseuM84rjDG/GqMyTTGZBljZhtjLgyw3cE+HU0MruB8XnfX+y5I+UnGmA/dePOMMXuNMT8aY642xsRVIu5YY8wtxpi5xph9/p1fuOU3u89ntjEm3d3+bP99lXMORxhj3jTGbDTG5LqvyxxjzB3GmKZBtinVWYUx5mz3WOnGmBw3nhuMMeX+j6jO9VDdmCvLff1f99n3XmPMLGPMjcaY+ADrW2OMxWnuBvCqz/Vj3fuiqnL8JPc1neFe7/nGmC3u/C3GmLZBthvrvic811qaMeY7Y8wlxpjoINvc58Y43Z0/xng7xck1xiw3xtwb7HVwtzneGPM/N8Z847y/1hljvjbG3GqMSXHX6+Y+T54mgWP8nidrjJnks99KvQ+MMVFu3E8ZY37xiWOP+5xdZYyJLSf+gJ3KeOL1vIbGmLbGmH8bY9a7z81OY8y7xpgyNYzG6ezkhwDH8EyTg8Xjt5+F7vqPV7DeOHe9YmNMF5/lFXYqU53rpoJYYo0xpxpjXnBfr+3u67HLGDPNGHOeMcYE2G4DlXzOTCU6lTHG9DTGPGeMWW2cz5VMY8x8Y8w9xpjm5TwXnvczxphexphXjDGb3edmizHmRWNMx3KO288991XG+SzMdbf/xRjz90DXi0iDYK3VpEmTpgon4I/AAZwONyyQB6QBBT7Lhvptk4zzJcFTXojTSUexz7JHghxvulv+APCx+3cBzv0n1mf6W4Btl5S3b3edJJzaIAtc7FeWAEzxO06GX9yzgZblxP1PnPvRPHGnu9uP9Tn+DL/nJt3nGH/32dd9AY4TBfzbL8Ysdz+e+RVA1wDbTnbLJwPPuH97OlDx3d9rIb4eqh1zJa/Rm/xeo31Avs/8QqC93zY73KnI53Xe4TN1rsLxhwGbfI5XhNM5Ta7PshsDbPe4T3mx+zr4PiffAc0CbHefWz4duM3d1rO97/PwPRAdYPt7/F6LbPf18F3muV47u8+H5z2T7/c87QDOrcb7oFuA62Gf37IfgYQgz3mpOH2W++73JGCnzzn6vh4ZwBC/bX9z4/Ss43+e/67k9XCru/32QM+/z3qvuuv9EOx9GmS76l43k9zyDQHKxvo99xlApt+y94Go6j5n5R3fLT/H7zXK9JvfBPSvIPaj8V7LmZT+XNoKdAyw/XF+x8mn7Gdimc9iTZoawhTxADRp0lT3J/cLlecL5s/AkZ4vBECcO/8CMMBvuw/wJgvXA4nu8nbAyz7/ZK8KcMzpblk6zhfEi3G/FAKdgE/xfunu7bft7T7/+KOCnNOF7jr7gaZ+ZW+4ZWuB84Hm7vImwKnucgt8VE7cWe40ySfuVkCK+/fzPvHf7okBaI03afJ8GbkvwHEecMt2Atf47DfW/WI03y2f5/8c4P2ime6+Njf5nGMr4EWf12ZcCK+HasdciWv0ZJ+YPwa6+8TzR7xfamcSODna4JZPquZ7pDOwG+8X1nPxXu8GGIBTu3aB33bX+cT9X6CduzwJp6MbzxfZdwMc8z6f66QI50eE1m5Zc+BvPvu+1G/brniT4MeADj5lye5r+CwwPMgxp1fwfEyncu+DTsCbwCmeZe7ypu42W939PB7kOJ7zG+u3vJtPWTrOdTrCLYsBjgW2ueU/BtjvWM/21bke3H20x5ugTQiyToLPtXlJkPfp5ADb1eS6mUTwhPBQnM+mY3E/E9zlKcCf8f4g9+fqPmcVHH8Y3h9xfgYGucuj3GvE85qtoezndsnx3df8E6Cfz+fAOT7P9esBjr3GLZsGHOSzvAkwEOcHlGp9PmjSVNeniAegSZOmuj25X57Wuf8ofwLiKrndYT7/nK8Mso4nYdwNNPErm+6z/dEBto3H+2XxTr+yjni/7I4PcuxpbvkbfstH401aAtYO4XyJ9dSUDC0n7lOCbN/FJ767gqwz2Wc/9/mVdcP5opmDX+2GzzrNgM3u9qeXs+9JQbaf65a/GKLroUYxV2L/y/DWJgVK+E7xOeezA5RvKO/5qMTxPT8ipAW7bgJsk4BTg2iBt4Osc71P3MGSszLXiM86H7rl3/gtP8ddvrKK5+k55vQK1qvwfVDJ443A+8NNkwDllUkIlxOghtHvmujkVzbWU1bd2N39fFXB63ueW56DX20eQRLCEFw3kwiSkFXifM52t10ToKxSz1l5xwe+dMtW4/6g4ld+MN5k99Zgx8epFS/zo5LP85IDxPgsT/XZtn158WvS1BAn3UMoIhU5Guju/n2TtTa/ktud6z5uAV4Kss7d7mNrnOY6gcy01v7gv9Bam4eT1AEM9ivbivOFAJzaoVKMMe2BY9zZN/yKL3Mf37LWbg4UkLV2C977ZY4PEvdSa+3UIGVn4fzinQM8EWSdB4IsB+cLVTTwlbV2YZAYs3BqysqLcTPwWpCyT91H//swq3s9TCI0MZdhnHtF+7uzD1priwLseyowx509r7L7ruTxk/Be7/8Mdt0EcBxOzQsE76b/P3iHHTg/yDp5wKNByj5xH/1fx33uYzM3/tpS3vugQtbaucAunJqvodXczWPW2gMBln+JUxsF3l5mQ83z+XK6MaZZgHLP59PH7vVfGaG6bqrjc/expzGmXQj3izGmBd73/SPW2jJDiVhrfwf+586W9z7+u7W2OMByz/shAejtszwLp9UDODW7Io2KEkIRqcgo93GH++Wsska4jz8E+ceMtXY53rGoRgRaB/i1nGNscx9TApS97j6eEeAL7/k4yck24Fu/siPcx8uMMTuCTThNqsBpehfIzHLiHuY+zrXWZgdawVq7FidhC8QT4/gKYrykghh/s9baIGXBntvqXg+hijkQz7VTiHNfZjDf+K0fKiNwmr0CVCX58cSx2Vq7KtAKbnL7vd/6/pZaa4P1jhrsdZyDU5vZHvjVGHOd26FGmQ5Daqi89wEAxpg4t/OYr40x29wOQEo6JMGpvQGnZr46An6GWGsLcVonQODPkFD4CCfZSMD5IaiEcToYGu/Ovk7lheq6CcgY08wYc5txOvXZ5XYq43ktfJO06r4ewQzDaV4NZT+XfXnex4NN8A6Hgv3f2Obzd8lr7v5g4Olc7CtjzP3GmMNMkM7DRBoajUMoIhXx/Aq8sYrbeb7EVTT48BacJp6pQcrL+9W80H0M9KXgfzi/kjcFzqR0TaDnV/m3AiSrHdzH5u5UkcQgy3eVs00b93FbOeuA89x1DrDcE2OSO1UkWIzVeW6rez2EKuZAPNdOmltzHMwWv/VDxbempCrPS1XeI77r+6vM61jq/721dp8x5jzgbZz7o552izKMMT/idBzynrW2oILYKlLe+wBjTCrOl3/fGrpcnGTVU9PbBucH7OrWZFb3M6TGrLU5xpgPcWrI/4jTDNTjPJwfpnbgTXIqI1TXTRnGmD44iZFvspeDU6Ps+az09JQb6ppl3zjLOzfPecXgJHU7/VcIVttqrS30+c3D/zW/HKdlxBCc1it3A/nGmN9wahZfthqjVBoo1RCKSEWC1SDVaW7Nm6dp0UWe5caYQTj/8KFsc1FwvqABXG2tNZWYJgUJoUyzxUBhVmKdQDwx/quSMY6t5nECqY8x17b6+h75Fqf570U4TYdX43QocwrOe+P38rror6SK3gdP4CSDe4BLce7fSrDWtrHWtrPWtsP7w0moay/DxVP7N9YY4/sDj+eHqbcDNXOOkFdxksENwESglbU2yVqb6r4WvtdDfX09ArLWbsKppZwAPIXbuRVO64aHgTXGmHGRi1Ck9ighFJGK7HAfq9KED7w1AxU1K/KUl1uTUE2ehG+czxdbz5ewBdbaxQG2qe75VoWnmVqHctcq/eXLVzhiDKa6x67NmD3XTmsTYKxBH7V1re3w+bsq5xfx94i1Ntta+4a1dpK1to97rDtwaul8aw5Dzm3ud6Y7e5219lVr7Q6/daJx7jGuz6bjNP+OAi4AMMYMwNt0vCrNRaGWrhs3WfU0CT/PWvtBgBqxkN436Mc3zvLOzVPmGaonZKy1xdbaadbaG6y1I3BqIC/A6Tm4JfC2mpFKQ6SEUEQqMst9bGeMqcq9KJ77y442QQY4dwf59SQ9v1UzvvJ8j9O8KAo4343D08FCsC9hnnueTq6FeDzmu48jgnXoYYzpQeDmouCN8VhTzqDjtaS610Ntxuy51mKAMeWs57nvM9TX2ly8nZOcUsXtADq5TfXKcBOio93Z2niPlGKt3WqtfRhnKAoo29mTp9lgKGqH2uB06Q/we5B1jvRZJ5xKmpLX9L5K9z7dN93ZP/o9LgrWyVI5auu68f28CfZ6HBtkOdT8OZvvs49jylnPE8PCEDRpLpe1Nsta+zbezsbaUnsdEIlEjBJCEanIDzjDDAA8UYVfR991Hzvi3JsRyP3uYxrldyJQLe79gW+5s38ExuEdkuLtIJu94D4eZIy5urz9G2OSqvlr8f9wvvgkATcEWefOcrZ/BefX8dY4Y82VF2OcMaZpNWIMprrXQ63FbK1dhDPsBMBd7pdh/32eiDMUCsA7ld13JY+fg/d6/4tfs8DyfIPTVBKC9xb5J7w1ySGLu4KaVABPr5z+99hmuo8tQhCGZ0w48DbjLmGMiQEeCsFxqiPT5+8WIdif5weoAe4PKRf4La+K2rpuMnz+DvR6NAPuKmf7Gj1n1tp9eHuOvs0YU+Y+YmPMELyd84Ty/VDR55hvL7UBO0kTqc+UEIpIudx7WzyDIB8JfGeMOdJT6+d+eR9rjHnTbQbl2W4OzhhoAE+7vRgmutu0M8a8iHOPCsDd1trcWjoFT7PRQcA/3L+/ttaW6YjAjXsGzn00AM8aY55wa+sA54u0MWakMeZhnA5EqtxBibV2I/CyO3u/MeZWTwJkjGlljHkc536qfUG2X4t3WIrbjTGvG2MO8okxxhgz1BhzD85gy0OrGmM5sVf3eqjtmO9wH0cDHxhjurv7jTXGXID3y+MsvENbhNKdOD9stAJmGmPOMcYkuDEYY8xBxphHjDElw6C4PRve586eZ4x53u15EmNMojHmz8CTbvl71tp5IYz3DmPMl8aYPxpjSprnudf3OcBt7qLP/bZb4j4ONMaMogbcnlE9NcePG2PG+VxHBwFf4PSQGbAn3lq2Cm+t7+UhqCVcgbdm7zmc2rjyfpgqb1+1dd0sx2kaCfCKMWa4p8AYczhO09eW5WwfiufsLpxxBnsB09x7vjHGRLk/6nyB0xJgLfDfauw/mFHGmEXGmJuMMf19rkPjXufPuettARaF8LgidYOtA4MhatKkqe5POB1P5OIdvNfTE2CBz7KhftskU3qA6gKcez6KfZY9EuR4nu3uKyem+6jcINnzfI5ngT9UsH4c8KLfNllu7EV+yztWNW53vaY4A7t79uO5H8bz3DyAM4SCBf4SYHuDU8Pq+1zmuK9JoV+MR/htO5kAA177rTOJcgavrub1UO2YK3mN3uS37704Y/R55hcBHYJsu8FdZ1IN3iPDcL4w+r6maTi1C55lNwbY7nGf8mL3OvB9Hr/Hb9Dyyl7/BBksnNKD2ntehz1+z98yoJ3fdjHACp910t3nbgNwdjXeB8NxBp33vY48NYcFODX7QV8bn+3G+i3v5lPWrZzjl7fvl3z2kY3zA9AG4NFqXh++g8VbnDE5y1t/MuW8T2tw3UwiyHsbp6m87z6y3cm6r9MxwZ7zyj5n5R3fLT+X0u/bDEq/hzYB/St7rVfmmqH0oPYWJ7H1/zzLAEZX9/NBk6a6PKmGUEQqxVr7OtAP55fnZThfdhNw/uF/jPPFbbnfNhk4XyAuw/mCmIWTCO3AqT082lp7G7XPt1lWJt7BiQOy1uZba6/A6WBhMs6v0dE4se/COZf7gcHW2oq6fg92DM+Xq9twEpV8nIRpBnCmtfZuvM2u9gXY3lpr78EZcPw/OM99EU4SvhenJuwRYJS1tsKx4KoRf3Wuh1qN2Vr7BE6N0ps4nXgk4nyR/AUnWTzEWlvRUB/VZq2dD/QH/uIeMwtohtOJ0HTgZgLUCFlrb8ZpzvwhThf6Td1tf8CpKT7OVn7Q8sp6AbgSp+Z0CU5C2BzndfgJuBEYZv06ebHO2H3H4HzxX4/T7LmrO1W5abJ1aq8OxRnmIg2n5VKWOz/KWhuoJ+BwuRYncfZ0PtUF5zyr28nNOzgJhkd1mouWqI3rxlr7GXAUTs3wPpwfANJwWk0Mt9Z+F3xrIATPmbX2PZwOjf6L89kbj/P5sgC4FzjIOmPYhtJvwDk4NYHzcM65Oc4PFAtwehntb639KcTHFakTjLU20jGIiIgftwnpHpzayqP0RURERERqg2oIRUTqpptxksF0wtC7pIiIiDROSghFRCLAGNPMGPOuMWaCMaaFz/KuxphH8HYa8aStvQ53REREpJFTk1ERkQhwk8C9Pos89/o081n2IU4HOIXhiktEREQaFyWEIiIR4I6x9iecgb8PwhmkOwGnM4O5OB1OfGj1IS0iIiK1SAmhiIiIiIhII6V7CEVERERERBopJYQiIiIiIiKNlBJCERERERGRRkoJoYiIiIiISCOlhFBERERERKSRiol0ABJ6xpgdQCKwOdKxiIiIiIhIic5AjrW2XaQD8dCwEw2QMSYzPj6+Wc+ePSMdioiIiIiIuNauXUteXl6WtbZ5pGPxUA1hw7S5Z8+eA5YuXRrpOERERERExDVw4ECWLVtWp1rx6R5CERERERGRRkoJoYiIiIiISCOlhFBERERERKSRUkIoIiIiIiLSSCkhFBERERERaaSUEIqIiIiIiDRSGnZCREQkQqy1aDxgEZH6yRiDMSbSYdSYEkIREZEwKi4uJiMjg71795KXlxfpcEREpAbi4+Np2bIlycnJREXVz8aXSghFRETCxFrLjh07yMjIiHQoIiISAnl5eezYsYPc3FzatWtXL2sMlRCKiIiESVZWVkkymJqaSvPmzYmOjo5wVCIiUh1FRUVkZmaya9cu9u3bR1JSEs2bN490WFWmhFBERCRMMjMzAUhJSaFVq1YRjkZERGoiKiqKVq1aUVhYSHp6OllZWfUyIayfDV1FRETqoZycHACaNWsW4UhERCRUPJ/p2dnZEY6kepQQioiIhIG1lqKiIsDphEBERBoGz2d6UVFRvew5WgmhiIhIGPh+SaiPnQ6IiEhgvp/pSghFRERERESk3lBC2EBZCzNW7Wbqwm2RDkVEREREROoo9TLaQK3dvZ+LX5lDu+ZNOHlwezVPEhERERGRMlRD2EB5mi/vyMxl7e762eORiIiIiIjULiWEDVR8rPel/Xn17ghGIiIiUrEbb7wRYwzGGB577LGA62RlZREVFUVSUhLFxcVhjrB8u3fv5tZbb6Vv374kJCSQkpLCsGHDuO222yq1/Z49e0hNTcUYQ69evYKud+DAAe655x769OlDkyZN6NChA5deeilbt24N1amISCOjhLCBSoyLLvn75zVpEYxERESkYgsXLiz5+9NPPw26jrWWQYMGERVVd77CzJs3j/79+/PYY48RGxvLaaedxsiRI0lPT+eJJ56o1D5uueUW0tLK/3+dm5vLuHHjeOCBB9i/fz+nnXYanTt35tVXX+Xggw9m3bp1oTgdEWlk6s6nqYRUYpz39tBf1qVTUFS3fkkVERHx5UkI+/bty8yZM9mzZ0/QdYYMGRLW2Mqze/duJkyYwIEDB/jkk09YsmQJ7777Ll988QUbNmxg1qxZFe7ju+++47XXXuOKK64od70HH3yQX375hcMPP5xVq1bx3nvv8euvv/LYY4+xe/duLr300lCdlog0IkoIG6iE2Gjiop2Xd39eIQs274tsQCIiIkFs2rSJvXv30rlzZy688EKKior4/PPPy6y3YMECAIYOHRreAMtx7733kpaWxiOPPMKpp55apvzQQw8td/sDBw7wpz/9iQEDBnDrrbcGXS8/P59nnnkGgGeffZamTZuWlN18880MHjyYGTNmMG/evGqeiYg0VkoIGyhjYES3liXzP61Ws1EREambfBM9T1IVqNloXashPHDgAG+++SZJSUlccskl1drH3/72N9atW8fzzz9PbGxs0PVmzpxJRkYGPXv25OCDDy5TfvbZZwMwderUasUhIo2XEsIG7MjerUv+nqn7CEVEpI7yTfQGDx5Mt27dmDZtGnl5eSXrFBUVsWTJEowxDB48OFKhljJ37lyysrI4+OCDSUhI4Msvv+Tmm2/mmmuu4cknn2TbtvLHAl60aBGPPfYYl1xyCaNHjy53Xc9zNGzYsIDlnuWLFi2qxpmISGOmhLABG92rTcnfCzbvIzO3IILRiIiIBObfFPSUU05h//79fP/99yXrrFq1igMHDtCzZ89SzSWrYuzYsSU9mVZ2mjx5ctD9LVu2DIDU1FROP/10TjzxRJ544gmee+45brrpJnr16sU777wTcNvi4mIuv/xyWrRowcMPP1xh7Js2bQKgU6dOAcs9yzdu3FjhvkREfGlg+gZsQIfmtEiMZV9OAUXFll/W7mH8wHaRDktERKQUT+2Xb0L49NNP88knn3DCCScAobl/cMKECXTr1q1K25Q3BMTevXsBp3lrdHQ0zz77LBMnTiQnJ4dnnnmGRx99lIsvvpj+/fuXifvpp5/mt99+49VXX6VVq1YVxrF//34AEhMTA5YnJSUBztAcIiJVoYSwAYuOMozq2YovFu8AnN5GlRCKiNRN1loycwsjHUaVNW8SgzGm2ttnZWWxbt06mjVrRo8ePQCnJq958+Z89tlnWGsxxoTk/sG//OUv1d42EM9YiIWFhTz00ENcc801JWWPPPIIGzduZMqUKTzyyCO89dZbJWWbNm3irrvuYsyYMUyaNCmkMYmIVJUSwgbu8B7ehHD2urJdeIuISN2QmVvIkL99HekwqmzhveNJTgjeGUpFFi1ahLWWwYMHlySWsbGxTJgwgffff5958+YxYsSIStcQ5ufnExcXV+14qsK36WqgTmUuueQSpkyZwowZM0otv/baa8nPz+f555+v8rFycnIClmdnZwPQrFmzSu9TRASUEDZ4I3t4m6Gs2JHJvpx8WiSG5x+liIhIRTyJnn/N3ymnnML777/PJ598wogRIwLWEObn55OUlMQLL7zAjBkz+Oijjxg/fjxTpkwJeKx//vOfrFixokrxXX755Rx55JEBy7p27Qo4zTjbtGlTptzTPHXXrl2lln/22We0aNGCq666qtTy3NxcALZu3crYsWMBePfdd2nXrh1dunQBYMuWLQFj8Sz3xCQiUllKCBu4XqlNaZUUx57sfKyFX9enc7yajYqISB3hf/+gx4knnkhMTAyffvop119/PTt27CAlJYXOnTuXrLN8+XIKCwu5++67ueiii/jf//5HixYtgh7rq6++KlNbV5GxY8cGTQg9wz8cOHCAvLw84uPjS5Wnp6cDBOwEZ9++fUFjyc3NLSnzJImeRHj+/PkBt/Esrys9sIpI/aGEsIEzxjCyRys+X7wdgF/W7VFCKCJSBzVvEsPCe8dHOowqa96kZl8lgtUQpqSkcMQRRzBjxgw+/vjjgOssWbIEgDvuuIPrr7++wmNNnz69RrH669KlC0OGDGHhwoXMmDGD8eNLv36epM5/3EBrbcD9bdiwge7du9OzZ0/WrFlTquyII44gOTmZtWvXsmDBgjIJ9AcffAA4NasiIlWhYScagZE9vc1Gf1mXHsFIREQkGGMMyQmx9W6qSYcynrEFo6OjGTRoUJlyT3Lzj3/8Ayhbi7h48WJSU1O59tprqx1DTd1+++0A3HrrrWzfvr1k+YIFC3jssccAyjQNrY64uDiuu+46wLkH0XPPIMDjjz/OokWLGDNmDMOHD6/xsUSkcVFCGCHGmARjzP3GmFXGmFxjzDZjzCvGmI6hPtbhPVJK/vbcRygiIhJpnrEFe/fuTUJCQpnyU089FXBqziBwDeGxxx5LVFTkvs6cf/75XHzxxSxevJgBAwZw0kknMW7cOEaOHEl6ejpXXHEFEydODMmx7rrrLg477DBmzZpF7969Offccxk5ciS33HILbdq04ZVXXgnJcUSkcVFCGAHGmCbA98DdQFPgE2AzcAnwuzGmRyiP17NNU1o3dTqS8dxHKCIiEmnB7h/06N27N/369SuZD1RD6N8cMxJeffVVXnjhBXr27Mn06dOZM2cOw4YNY/LkybzwwgshO06TJk344YcfuPvuu0lMTOTjjz9m48aNTJo0ifnz55cM2yEiUhUmWDt2qT3GmAeBO4HZwHhr7X53+c3AY8AMa+3YGux/6YABAwYsXbq0ZNm1b8/n80VOU5ZJo7px36kDq38CIiJSZcXFxaxcuRKAvn37RrRWqyHIzMwkOTmZadOmlbl3T0QknKry+T5w4ECWLVu2zFpbZ76M679RmBlj4oDr3NlrPckggLX2cWARMMYYE9KbAHyHn/htg2oIRUSkfvN0KKNeNUVEakYJYfgdASQDa621vwco/8B9DGk3YYd0a1ny9/LtmezPKwzl7kVERMJqyZIltGnThnbt1HO2iEhNKCEMP88d8YEHEvIuD+lPnn1Sm9HM7Rq82MLvm/aGcvciIiJhdeWVV5YZ8F1ERKpOCWH4dXEftwQp9yzvGsqDRkUZhnf11hLO3aCEUERERESksdPA9OHX1H3MCVLuGVioWUU7MsYsDVLUM9DCQ7qlMH3lbgDmbtR9hCIiIiIijZ1qCBsR3xrC3zfto7CoOILRiIiIiIhIpKmGMPw8vYomBilPch+zKtpRsO5q3ZrDAf7Lh3RqQWy0oaDIkpNfxPLtWQzqlFyZmEVEREREpAFSDWH4bXIfOwUp9yzfGOoDJ8RFc1BHbwKoZqMiIiIiIo2bEsLwW+g+DgtS7lm+qDYOPsK3Y5mN6lhGRERERKQxU0IYfjOBDKCnMWZogPKz3ceptXHwEd1SSv6euyEda21tHEZEREREROoBJYRhZq3NB55xZ581xnjuGcQYczPO+IMzrLXzauP4B3duUfL3zsw8sjRAvYiIiIhIo6VOZSLjQeBYYBSw2hjzE864g4cBu4FLa+vALRLjSs3n5BXRvElsbR1ORERERETqMNUQRoC1Nhc4GngAZzzC03ESwsnAMGvtuto6dlxMFLHRpmQ+O181hCIiIiIijZVqCCPEWnsAuMedwioxLoaMAwWAU0MoIiIiIiKNk2oIG6GkuOiSv3NUQygiIiIi0mgpIWyEEuO9FcM5+aohFBERERFprJQQNkKJPjWEuodQRERERKTxUkLYCPkmhLqHUERERETEsXjxYmJiYujUqVOkQwkbJYSNUFKct8moaghFRERERBw33ngjrVq1inQYYaWEsBHSPYQiIiIiIqV9/PHHrFu3jksvrbUhweskJYSNkHoZFRERERHxys/P59Zbb+Wf//wn8fHxkQ4nrJQQNnSrv4V100stSvDtVEb3EIqISB1w4403YozBGMNjjz0WcJ2srCyioqJISkqiuLg4zBEG9vjjj3PmmWfSu3dvkpOTiY+Pp2vXrlx00UUsXrw44DY5OTl8/PHHXHbZZfTt25cmTZqQlJTEkCFDuP/++9m/f3/Q4x04cIB77rmHPn360KRJEzp06MCll17K1q1ba+sUI8YYQ7du3UK2vw0bNmCMYezYsSHbZ3Xdf//9REVFlblGqno91eRagppdT3369An4+uzevZtbb72Vvn37kpCQQEpKCsOGDeO2226rcJ979uwhNTUVYwy9evUKWdzbt28nISGBa665Jug+n3zySdq0acO5555bYZwNjrVWUwObgKUD+vWx9r2LrL23ubVPHGRtXrb1eOSrFbbrHZ/Zrnd8Zm+bssCKiEjtKyoqssuWLbPLli2zRUVFkQ6nzhk7dqwFLGCPOuqogOv89NNPFrCHHXZYmKMLrlWrVrZJkyb20EMPtWeccYY944wzbJ8+fSxgY2Nj7dSpU8ts8+KLL5aca//+/e3EiRPt8ccfb5s1a2YB269fP7tz584y2x04cMCOHDnSArZ9+/b2nHPOsYceeqgFbJs2bezatWvDccphA9iuXbuWWrZ+/XoL2DFjxlR5fzXZNpR27NhhmzZtaidOnFimrKrXU3WvJWtrdj0tX77cAva6664rtXzu3Lm2VatWFrADBw605557rj3hhBNs165dbXR0dIXPzcUXX2yNMRawPXv2DGncf/7zn21MTIxduXJlmbIdO3bY5s2b21mzZllrrb333nttx44dK4zXoyqf7wMGDLDAUlsHcgbPFPEANNXCiwpLB/Ttae3fUpyE8N7m1n5zn/V49ofVJQnhNW/NsyIiUvuUEJavZcuWFrB9+/a10dHRNi0trcw6zzzzjAXslVdeGYEIA/v555/tgQMHyix/9tlnLWDbtm1rCwoKSpVNnjzZXnnllXbZsmWllm/bts0efPDBFrDnnXdemX3eeeedFrCHH364zcrKKln+2GOP1YlEJ9QaakL45z//2QJ2/vz5Zcqqej1V91qytmbX07/+9S8L2GnTppUs27Vrl23durVNTEy0n3zySZltfv3116D7s9bab7/9tuT9XV5CWN24t2zZYqOiouw555xTpuzSSy8ttVwJoaZ6PwFLBwwYYO3X93gTwr+1snbncmuttZNnri9JCCe9Uv6bU0REQkMJYXAbN260gO3cubN94IEHLGBfe+21MutdfvnlFrD/+c9/IhBl1fXs2dMCduHChZXeZtasWRaw8fHxNi8vr2R5Xl6eTU5ODppIDB482AJ27ty5IYm9LmiICWF2drZNTk62Bx10UJW3rer1FOxasrbm19ORRx5pmzVrVmq/V199tQXss88+W4WzcuTk5NiePXvaAQMG2FWrVgVNCGsa9zHHHGNjY2Ptjh07SpYtXrzYxsXF2fnz59u9e/favXv32jvuuMN26NDB7t27t8xzF0h9Twh1D2FDNuZ2SO7i/F1cAJ/fAtaWHodQvYyKiEiELViwAIChQ4dy6qmnAvDpp5+WWW/hwoUADBkyJGyx1URsbCwAcXFxld7Gc255eXns2bOnZPnMmTPJyMigZ8+eHHzwwWW2O/vsswGYOnVqTUIudZ9ddnY2N998M507dyYhIYFhw4aV2v+UKVM47LDDSEpKom3btvz5z3/mwIEDZfY5ffp0jDFMmjQp4DEnTZqEMYbp06eXG9t9991H9+7dAZgxY0bJPafl7TuYzMxMbrjhBjp37kyTJk3o378/TzzxRKl7U+fOnYsxhlGjRgXdz9///neMMdx7772VOu6UKVPIyMjgvPPOq1K8UPXrKdi1BDW7ntLS0pg9ezbHH398SSwHDhzgzTffJCkpiUsuuaTyJ+X629/+xrp163j++edLzjOQmr4Pzj//fAoKCpg8eXLJsjVr1pCfn8+wYcNo2bIlLVu25F//+hfbtm2jZcuWvPLKK1U+n/ompuJVpN6KS4ITH4Z3/uDMb/wZFr5LUvzYklWUEIqISKT5JnqDBw+mW7duTJs2jby8vJLe/oqKiliyZAnGGAYPHhzJcCvljTfeYOXKlfTu3ZvevXtXert169YBzpf/lJSUkuWe52jYsGEBt/MsX7RoUXVDLiU/P59jjjmG9evXc9RRR5GWlsaPP/7IGWecwVdffcXixYu5/fbbGTNmDMcffzw//vgjTz/9NHv27OGtt94KSQz+hg4dyllnncWHH35I27ZtmTBhQknZkUceWen95OXlMW7cONauXcu4cePIz8/nu+++4+abb2bhwoUlycKIESMYNmwYs2fPZunSpQwcOLDUfqy1vPzyy0RFRXHZZZdV6tifffYZQJU7tqnO9RTsWoKaXU9ffPEFRUVFJT/egJM8Z2VlceSRR5KQkMCXX37JN998Q25uLn369OGcc86hQ4cOAY+1aNEiHnvsMS655BJGjx7Nhg0bgp5TTd8Hnuf9888/54477gCca+eHH34otd7kyZP5/PPPmTJlCn369AkaT0OhhLCh63sC9D0JVn7uzH99J80mfFVSrIHpRUQk0nxrCAFOOeUUnn76ab7//ntOOOEEAFatWsWBAwfo1asXTZs2rdZxxo4dy4wZM6q0zauvvlqp2qdHHnmEpUuXkp2dzfLly1m6dCkdOnTgnXfeITo6usLtPf79738DMGHChFJd32/atAmATp06BdzOs3zjxo2VPlZ5Zs+ezbhx41i3bh1JSUmA8yX5kksu4eqrr2bPnj3Mnj2bESNGALBt2zYOPvhg3n77bR544AF69OgRkjh8nX766QwdOpQPP/yQfv36larlqYpffvmFwYMHs3r1alq3bg3A2rVrOeqoo3jttdc4/fTTOf300wG46qqruPLKK3nxxRd58sknS+3nu+++Y926dZxwwgl06dKlUsf+6aefiImJCVi75SsU11Owawlqdj19+umnREdHc+KJJ5YsW7ZsGQCpqamcfvrpfPLJJ6W2+etf/8rLL79cpma0uLiYyy+/nBYtWvDwww9XeE41fR/06NGD1q1bM2fOHHJzc2nSpAmtW7cuk6BPnz6d+Pj4OtEjbTgoIWwMTvgXrPsBCnIgZw99lzwOnAJAjoadEBGRCPP86u+fEH7yySclCaF/0lgdEyZMqPIwBuV1fe9r2rRpfPfddyXzXbt25fXXX2f48OGVPtYXX3zByy+/TGxsLA888ECpMs/wAYmJiQG39SRtWVlZlT5eeaKionjuuedK9gtw0UUXcdttt7FmzRruuuuukmQQoEOHDlxwwQU88cQT/Pjjj7WSEIbSo48+WpIMAvTs2ZO7776bq6++mmeeeaYkITz//PO59dZbeeONN/jXv/5VKrF66aWXALjiiisqdcxdu3axc+dOunfvTkJCQrnr1vR6Ku9agupfT/n5+Xz99deMGjWKVq1alSzfu3cv4E0Wn332WSZOnEhOTg7PPPMMjz76KBdffDH9+/cv9R5++umn+e2333j11VdL7S+YULwP+vbty8yZM1m+fHmFiXljoYSwMWjRGcb+Bb65B4DUVe8wzPRlvu2jGkIRkbrCWsjNiHQUVdckGYyp9uZZWVmsW7eOZs2alSQRY8eOpXnz5nz22WdYazHGhOT+wb/85S/V3rYi3377LQD79u1j8eLF3H///YwZM4YHH3yQO++8s8LtV6xYwYUXXoi1lkceeSTi90l269atTFO5qKgounbtSlpaGuPHjy+zjef12759e1hirK6UlBSOO+64MsvPO+88rr76ambNmkVxcXHJmJcXXngh//nPf/jwww85//zzAec+uo8++oh27dpxyimnVOq4u3btAqBly5YVrluT66k2r6UffviBrKysMufsufeysLCQhx56qNR4f4888ggbN25kypQpPPLIIyVNijdt2sRdd93FmDFjqnwPaE14ms/u3r076Dr33Xcf9913X5giijwlhI3FyGtg4buwy6nSvyP2Xc7Nv4cD+UUl/2xFRCSCcjPgX10jHUXV3bERElpUe/NFixZhrWXw4MEl/4tiY2OZMGEC77//PvPmzWPEiBEhqSEMhxYtWjB69Gi++OILDj/8cO6++27Gjx/PIYccEnSbrVu3MmHCBPbu3cvNN9/MDTfcUGYdTzPZnJycgPvIzs4GoFmzZiE4C+jYsWPA5Z44ApV7yvLy8kISQ23p2jXw+yw5OZkWLVqwb98+9u7dW1JjddVVV/Gf//yHF198sSQhfP3118nPz+eSSy4hJqZyX6czMpwffKryGlX1eqrMtQTVv548nbX43j/ouz8gYKcyl1xyCVOmTCnVZPvaa68lPz+f559/PmAMoYzbV/PmzQEn2RaHEsLGIjrWaTr6mvOLzgizkqbksL84kfyiYuJjKn9/g4iISKh4Ej3/WoxTTjmF999/n08++YQRI0ZUuoYwPz8/aC+M//znP1mxYkWV4rv88sur1GGJR2xsLOeeey7z5s1j6tSpQb/Ap6enM378eDZu3Mgll1zCo48+GnA9zz1qW7ZsCVjuWR4s2amqqKjyO6KvqLwqfHv2rIsGDRrEqFGjmD59OqtXr6Z37968/PLLGGO4/PLLK72f5ORkoHrNeitzPVX2WoLqX09Tp06ld+/e9O3bt9Ryz3qJiYm0adOmzP48TbU9taTgdLDTokULrrrqqlLr5ubmAk5y67mH791336Vdu3YheR94EvMWLVoEXaexUULYmHQ9Epq0gNx9RBvL8KjVzCgeQk5ekRJCERGJCP/7Bz1OPPFEYmJi+PTTT7n++uvZsWMHKSkpdO7cuWSd/Px8kpKSeOGFF5gxYwYfffQR48ePZ8qUKQGP9dVXX1W5U5mxY8dWKyEESu5RC9Y0bf/+/ZxwwgksW7aMM888kxdffDFoix1PIjx//vyA5Z7ldbEHVk+C7rn/y9/mzZvDGU5JxyT+MjMz2bdvHwkJCWWShauuuopZs2bx0ksvceqpp7Js2TKOPfbYKt0rmZqaCjiJW3WUdz1V5VqC6l1PCxcuZNOmTdxyyy1l1vfci3fgwIFSvQN7eM7Zv0Ooffv2BX1P5ubmlpR5ksRQvA889zsGSlwbKyWEjUlUFHQ5HFZ9CcAhUSuYUTyE7PxCWiZVfowkERGpBU2SneaX9U2T5BptHqyGMCUlhSOOOIIZM2bw8ccfB1xn+fLlFBYWcvfdd3PRRRfxv//9r9xf/Ssa5y7UPF9me/bsWaYsLy+P0047jTlz5nD88cdX2HvkEUccQXJyMmvXrmXBggVlEugPPvgAoNL3s4VT+/btAaenWH/p6elBv9wH4kkuCwur3wfCnj17+O677zjmmGNKLX/33XcBOPzww8u8FhMnTuSmm25i8uTJJbVQle1MxiM1NZV27dqxefNmcnJygnaMEkyw66mq1xJU73ryjA0a6Brr0qULQ4YMYeHChcyYMaPMPaae2H07cbHWBoxtw4YNdO/enZ49e7JmzZoax+1vxYoVxMfH079//6DrNDq1Oeq9pshMwNIBAwbYgH7+t7X3Nrf23ub217sPsV3v+Myu3JEZeF0REQmZoqIiu2zZMrts2TJbVFQU6XDqhMLCQpuQkGCjo6NtTk5OmfJHH33UArZbt24WsDfddFOp8jfffNMC9qmnngpXyKX8/PPP9ssvvyzzeubn59unnnrKRkVF2YSEBLtp06ZS5YWFhfaMM86wgB09erTNzs6u1PHuvPNOC9hRo0bZ/fv3lyx/7LHHLGDHjBlTZpt7773XAvbiiy+u1DHWr18fdF/WWjtmzBgL2PXr15cpe/XVVy1g77333jJlXbp0sYD9+OOPS5bt37/fnnXWWRawgP3hhx9KbQPYrl27llqWl5dnY2Njbbt27WxhYWGlzsn/3AA7dOhQm5aWVlK2bt0627FjRwvYDz/8MOD2N998c8n2bdq0sXl5eVU6vrXWTpw40QL2559/LlNWneuputeStVW/ng455BCbkpIS9Hl/6623LGAHDRpkt23bVrL8999/tykpKRaw77//foVxeV6nnj17hiRuX2vWrCl5rkKpKp/vAwYMsMBSWwdyBs+kGsLGpuuokj+HmLXEk092nnoaFRGR8POMLdivX7+A3fCfeuqp3HrrrSUDVfvXEC5evJjU1FSuvfbacIRbxurVq7nkkkto3bo1w4cPp1WrVqSlpbF48WK2b99OkyZNmDx5cqlmrgDPPPMMH330EeA0A/TtkdGX/9AId911F99++y2zZs2id+/ejB49mo0bN/Lrr7/Spk0bXnnllTL78NyfFxsbG6rTrpZ7772Xyy67jLPOOoujjjqKpk2bMmfOHJo3b85pp51WZty6YOLi4pgwYQJTp05lyJAhDBs2jLi4OI444oiAnZkEMnLkSPLz8+nVqxfjxo2joKCA7777jpycHC688ELOPPPMgNv96U9/4oknnsBay8UXXxz0XtXynHTSSUyZMoXp06dzxBFHlCqrzvVU3WsJqnY9bd++nblz53L++ecHrX08//zz+frrr3nttdcYMGAAo0aN4sCBA8yaNYu8vDyuuOIKJk6cWOXnzF913gcenlYCJ510Uo3jaFAinZFqCv1EeTWEhfnWPtiupJZw4l8esTNX7w68roiIhIxqCMt65513LGD/8Ic/BF2nX79+JbUyCxYsKFV20kkn2fPPP7+2wwxq3bp19q9//as94ogjbPv27W1sbKxNSkqyAwcOtNdff71dvXp1wO08tXYVTYFq4XJycuzdd99te/bsaePi4my7du3spEmT7ObNmwMe69RTT7WA/fHHHyt1TrVVQ+gpP+igg2xcXJxt27atvfzyy21aWpq9+OKLK11DaK21O3futH/84x9tu3btbHR0dKVrQH3Pbd++ffaaa66xHTp0sHFxcbZv37720UcfrbDWsXPnzhawK1asqPB4geTk5Njk5GQb6Htada6nmlxLnngqcz298MILFrDvvfdeuedXXFxsX3jhBTt8+HCbmJhok5KS7OGHH24nT55c6eeoohrCqsTtb9y4cTY2Ntbu2LGj0vFURn2vITTWBm6/K/WXMWbpgAEDBixdujTwCq+dCuudttyPFJzD0Ase5LgBbcMYoYhI41NcXMzKlSsBZ2DkUPbS2Fh17dqV66+/nltvvTXSodRJxcXFtGrVihEjRvDNN99EOpx6b/bs2YwaNYoxY8bU6H7Um266iSeffJK5c+dWeqD5SDvllFOYNm0aaWlpJcM21Ddbtmyha9eunH322bz33nsh3XdVPt8HDhzIsmXLlllrB4Y0iBrQf6PGqKu3icJhUcvJ0eD0IiJSz2RmZrJp06Y62atmXfH777+zb98+7r///kiH0iA89NBDAFx33XU12s///d//0bRpU/7xj3+EIqywGD16NE8//XS9TQYBHnnkEaKiovR+CED3EDZGXQ8v+XNY1GqmHqjbA8iKiIj4W7JkCVA3h1moK4YPH45agtXMrFmzePnll1myZAlz5sxh2LBhQe8xrKzU1FRuu+027rvvPhYvXsygQYNCFG3tuf322yMdQo1s376dF154gSuuuKLMGIqihLBx6jiCQmKIoZCmJpeE9KVA5cfRERERibQlS5bQpk0b2rVrF+lQpAFbtWoVr7zyCs2aNeOkk07i2WefDUlz73vuuYd77rknBBFKZbRv354DBw5EOow6S01GG6O4RDYneH8dSdkzL4LBiIiIVN2VV17Jrl27Ih2GNHCTJk3CWktmZiafffYZXbt2jXRIIiGnhLCR2tzU23V3272VHxBWREREREQaDiWEjdT2FsNK/u6UuQB0j4GIiIiISKOjhLCR2pMyjGJrAEgqymD+vF94f+5mcguKIhyZiIiIiIiEixLCRiomqQUrbJeS+Q8+msLtHyzi6e9XRzAqEREREREJJyWEjVRiXAxzir0dyxwatRyAH1bsjlRIIiIiIiISZkoIG6mk+GjmFPcrmT80agVgWbc7k6Itv8PcV2H3ysgFKCIiIiIitU7jEDZSCbEx/OaTEHYw6TwZ+yyHRy0j+qV9zsLYJLh5KSS0jEyQIiIiIiJSq1RD2EglxUezmxasK/YO6Ht69Czamn3elQqyYfvC8AcnItIAGWNK/i4uLo5gJCIiEkq+n+m+n/X1hRLCRioxzqkc9m02GlDm9jBEIyLS8BljiIuLAyA7OzvC0YiISKh4PtPj4uLqZUKoJqONVFJ8NAAvFZ3IkdFLSCCPn4sH8X3RUK5uOZd+2XOcFTO3RjBKEZGGpVmzZuzZs4edO3cCkJSURFSUfpsVEamPiouLyc7OLvlMb9asWYQjqh4lhI1UkltDuMZ24si8p0qVHVO0hZJ6wyzVEIqIhEqrVq3Izs4mNzeXbdu2RTocEREJkSZNmtCqVatIh1EtSggbqYS46KBlK3OaeRsTq8moiEjIREdH06VLF/bs2UNWVhb5+fmRDklERGogLi6OZs2a0apVK6Kjg3+/rsuUEDZSnhrCQDYXtoA4dyZLv2CLiIRSdHQ0qamppKamYq3FWhvpkEREpBqMMfXynkF/SggbqSaxURgDgb6H7LQp3hnVEIqI1JqG8mVCRETqL93J3kgZY8okgwM7NAdgBz7jDmbvgqLCMEYmIiIiIiLhooRQSozs4dwIu8O3htAWw/6dEYpIRERERERqkxJCKdGnbVMA8ogjyzT1FmTqPkIRERERkYZICaEA0Cw+hl6p3rFTthf7NBtVxzIiIiIiIg2SEkIBYECH5nRrlVgyv61YHcuIiIiIiDR0SggbsRuO6Q1AdJThwdMPolmT2JKyHVY1hCIiIiIiDZ2GnQgjY0w/4DRgAjAISAb2ALOAJ6y1P4UznuvH9eKgjsl0bJFA77ZOc9HYaENBkWUnqiEUEREREWnolBCG17dAR2A/8AuQDgwAzgBON8bcbK19MlzBxERHcdyAtqWWJcXHsC+nwK+GUAmhiIiIiEhDpCaj4bUCuAhoY609zlp7rrV2EHAVYIBHjTEDIhlgUpzzG0GpoScyt0YoGhERERERqU1KCMPIWnustfYNa22u3/L/Al8D0cDEiATnSoyLBmCnbw1h5nbKjGIvIiIiIiL1nhLCumOh+9ghkkEkxgeoISw8ALn7IhOQiIiIiIjUGiWEdUcP93FHJINoGu/UEKbTjKIob6+j6lhGRERERKThUUJYBxhjegInu7OfRjKWxDhPP0OG7LhUb4GGnhARERERaXDUy2iEGWNigMlAPPCetXZeFbZdGqSoZ3XjSXLvIQTIimtD81y3QxnVEIqIiIiINDhKCKvAGPMR0L+Km11krZ1TTvlTwJHAOuCa6sYWKknx3ksiI6Y1HT0zGnpCRERERKTBUUJYNd2BvlXcJjFYgTHmTuBqYCdwvLU2vSo7ttYODLLfpTjjG1aZb0KYHt3KW6ChJ0REREREGhwlhFVgrR0aqn0ZY64CHgQygAnW2jWh2ndNJPo0Gd1jfBNC1RCKiIiIiDQ06lQmAowxfwCeBXKAk6y1CyIbkVdSnPc3gh34JITqVEZEREREpMFRQhhmxpgTgdeBQuAMa+3MCIdUim+T0VJjEaqGUERERESkwVFCGEbGmCOADwADnGut/TrCIZWRFO9tMrq1uIW3ICcNCvPCH5CIiIiIiNQa3UMYXp8BCcB64HRjzOkB1vnZWvtSWKPykejTZHRrQXLpwqwd0LJrmCMSEREREZHaooQwvFq4j93dKZiIJYS+4xDuyzeQ2NqpHQRn6AklhCIiIiIiDYYSwjCy1ppIx1AR33sIc/ILIbW9NyHU0BMiIiIiIg2K7iGUUnzvIczOL4JmHbyF6lhGRERERKRBUUIopfjeQ5hfWExRs/bewiwlhCIiIiIiDYkSQinFt8koQEFiO+9MpsYiFBERERFpSJQQSimJPp3KABxISPXOqIZQRERERKRBUUIopcRGRxEX470sDsT7JISqIRQRERERaVCUEEoZvkNPZMb61RBaG3zD7DSY8yLsWFKL0YmIiIiISKho2AkpIyk+hr05BQBkxLT2FhTlQ84eSGodeMMpk2DDTxDXFK6eCS271XqsIiIiIiJSfaohlDKSfHoazSQJYhK8hcGajeakO8kgQP5+mPFILUYoIiIiIiKhoIRQykj0H4uwuc9YhME6ltkyt/T8wrchbU0tRCciIiIiIqGihFDK8K0hzM4vLJ0QBqsh3PJb6XlbDNP/XgvRiYiIiIhIqCghlDKSfGoI9+UUMGdPvLcwaA3hnLLLlnyoDmZEREREROowJYRShm8N4RPfrGL+3gruISwugi3zvPNxzbx/T/9HLUQoIiIiIiKhoIRQyvC9h7Cw2LLDtvQWBkoId6+A/CznbxMFJz3qLVvxGWydX0uRioiIiIhITSghlDKS4kuPRrLDpnhnMreW3WCzT3PRtgNh8LnQfoh32fcPhjhCEREREREJBSWEUoZvk1GAdba9d2b3CsjwSwp9O5TpdCgYA+Pu9i5b+x1snFULkYqIiIiISE0oIZQyEuOiS82vsp3YWJzqXbDsk9Ib+CaEnQ91HnsdC51Hepd/eTvk54Q4UhERERERqQklhFKGf5NRMHxe7JPcLfvY+3dOOqSt8s53OsTdxMC4u7zLdyyGT64Fa0MdroiIiIiIVJMSQimjbEIInxcd5p3Z/Ku32ehWn95FE1tBSg/vfPfRcNhV3vml/4MffTqcERERERGRiFJCKGUk+TUZBVhqu7HR+jQbXf6p8+jboUynQ5yaQV/jH4IeR3vnf3gQlk8NYbQiIiIiIlJdSgiljMS4sjWEYPjCt5Zw6UfO4xa/hNBfdAxMfBVSenqX/e9PGrBeRERERKQOUEIoZSTFl60hBPisyOc+ws2/wr7NpQek93Qo4y+hJZz3LsQnO/MF2fDOebB9YYgiFhERERGR6lBCKGUEuocQnGajhcldvQt+fKT0gPQdhgXfaZs+cPYrznoAGZvgv2Pg42sgc3uIIhcRERERkapQQihl+I9D6GXI6nGyd/b3N7x/tx0I8U3L33HvY+H4v/sssLDgLXh6GEz/J+RnVzdkERERERGpBiWEUkZikCajAGldT/DO2GLv34HuHwxk5NVwwYfQpp93WUEOTP8HPHcE7F4VfFsREREREQkpJYRSRvAaQtjbfAC07Fa2oFOQ+wcD6X0sXDUTTnocElv77Hw9vHwsrJtR+X2JiIiIiEi1KSGUMqKjDE1iA18aOQVFMOD0sgXBOpQJepAYOOQy+PPvMOrP3uW5GfDmmTD/9artT0REREREqkwJoQQUrJYwt6AIBp5eeqH/gPRV0aQ5jH8AznkdYhKcZcWF8On18M09UFxc/vYiIiIiIlJtSggloGD3EebkF0H7oaWbjQYakL6qBpwGl3wOTdt6l838N0y5CPJzarZvEREREREJSAmhBBSshjAnv8hJ/oac713Y69jQHLTjcLj8O0gd6F22fCpMPhGydoTmGCIiIiIiUiJ47yHSqAUbizC3oMj548ibIDoWomJg+CWhO3CLznDZNPjgUlj9tbNs2+/w4jg4/z1oNyh0xxIRERERaeRUQygBHdnL6f0zOsowpHOLkuU5+W5CGBMHo2+GI/7sdBATSvHN4A/vwGFXeZdlboVXJsCqaaE9loiIiIhII6aEUAL68zG9eWXSCL68YTQD2jcrWX7AU0NY26Jj4IR/wQmPgHEv0/z98M4fYMYj6mxGRERERCQElBBKQNFRhnH92tKnbTMSYr01gAfyw5QQehx2JZz3HsQ1deZtMfzwILx1NmSnhTcWEREREZEGRgmhVCghznuZhD0hBOgzHi77unTPpmu/g+dHw6Zfwh+PiIiIiEgDoYRQKpTo0+NoTriajPprOxCunAH9TvYuy9oGr54Is54GayMTl4iIiIhIPaaEUCqUEOsdk/BAfhHWWhZu3sesNWks2ZpBxoGCMAXSAs59E47/h9O7KYAtgq/vgncvgAN7wxOHiIiIiEgDoWEnpEIJcT4JYUEhz89Yx7++WlGyLDba8NQfDuaEQe1rPxhj4PBroNMhMGUSZG5xlq/8HP67GCa+Bh2H1X4cIiIiIiINgGoIpUKJPglhTn4RXy8rPUh8QZHl4WkrseFsttn5ELjqJ+h1nHfZvk3wyvEw50U1IRURERERqQQlhFKhJn5NRvdm55dZZ31aNrPX7glnWJCYAue/D8fc4x2aoigfvrgVPrgEcjPDG4+IiIiISD2jhFAqlFiqyWgR6T4JYUpSXMnfb83ZFNa4AIiKgtG3wMVToWlb7/KlH8ELY2HHkvDHJCIiIiJSTyghlAr5diqzP7eQzNzCkvlrxvYs+fvrpTtI258X1thKdDsSrvoZuh/lXZa+Fl46Bua/oSakIiIiIiIBKCGUCvl2KrPHr7no2cM70bZ5PODcSzhl7pawxlZK01T448dw1O2AcZYV5sKn18HHV0Pe/sjFJiIiIiJSBykhlAr5jkPoKybKkJwQy7mHdClZ9uYvG9mfVxhw/bCIioZxd8KFH0BiK+/yhe/Af4+Cbb9HLjYRERERkTpGCaFUyLfJqK+WSXEYY/jDIZ2JjnJq5LbuO8AdHy4Kb4+jgfQ6Fv70E3Qe6V2WvhZeOg5m/huKiyMXm4iIiIhIHaGEUCrk22TUV8vEWAA6tEgodS/h54u28/rsjWGJrVzJHWHSZ06nM54mpMUF8M098MbpkLk9ktGJiIiIiEScEkKpUNAawkRvD6M3HtuHUT29TTQf/HwZv2/aW+uxVSg61hmW4uKp0KyDd/n6GfDcKFjxReRiExERERGJMCWEUqG4mChi3CahvnwTwugow1PnHVyqg5lr35ofcMzCiOg+Gq6eCf1P8S47kA7vngef3Qz5OZGLTUREREQkQpQQSqUEqiVsmRRbar5103ieOX9Yyf2E2zJyuen9BRQX15EhHxJT4Jw34JSnIDbRu3zuy/Di0RqzUEREREQaHSWEEWaMudsYY93pwkjHE0yg+wh9awg9DumWwh0T+pbMT1+5m2d+WFOrsVWJMTD8YrhyBrQb7F2+ewW8OA5mP6sOZ0RERESk0VBCGEHGmL7AnUAdqUILLrGSCSHAFaN7MH5A25L5J75dxYxVu2sttmpp0wcu/xYOv867rCgPpv0VXj8V9m2OXGwiIiIiImGihDBCjDEGeAHYB3wa2Wgq1iRgk9HACaExhkcmDqFrK6dZprVww7u/szm9jt2nFxMPxz8Ef/wImnoTWDb85HQ4s+BtJ3gRERERkQZKCWHkXA4cBdyCkxTWaYFrCGMDrOlITojl+QuH0yTWucT25RRw5RvzyMwtqLUYq63nOLh6Ngw4zbssLxM+vhreuxCy0yIXm4iIiIhILVJCGAHGmHbAw8B31tq3Ih1PZQS6h7BFkCajHv3bN+cfZw4qmV++PZPLX5tLbkFRyOOrsaRWMPE1OOMFiE/2Ll/xGfxnJKz8KnKxiYiIiIjUEiWEkfEUkABcHelAKishNqbMspQgTUZ9nXFwJ648qkfJ/Jz16Vz39u/kF9bBjluMgSHnwjWzoPsY7/Ls3fDOufDp9ZCXFbn4RERERERCTAlhmBljTgYmAn+31q6OdDyVFbiX0eBNRn393wn9OHt4p5L5b5fv5IrX55KTXxiy+EIquRP88WOY8C+IaeJdPv91eO4I2DgrYqGJiIiIiISSEsIwMsY0Bf4DrAL+FYL9LQ00AT1rum9/iX6dykQZaN6kcgmhMYZ/njmI43x6Hp2xajcXvvRr3Rm43l9UFIy8Cv70I7Qf6l2+byO8eiJMuxMKDkQsPBERERGRUFBCWAXGmI+MMSuqOB3qs4u/A52Bq621eRE6jWrxryFskRhHlDsAfWXEREfxzPkHc9Kg9iXL5m/ax8lP/8zCzftCFWbotenrDE8x5i9gPM+BhdnPwH+Pgi3zIhqeiIiIiEhNlL0xrBYZYy6qhd0usNYuqoX9BtId6FvhWqUlAriJ4bXAG9ba70MRjLV2YKDlbi3hgFAcw6NsQli52kFf8THRPHXewSQnxvL2r5sA2LrvABOfn83tE/pyyRHdia5Ckhk20bFw9P9Bn/Hw0VWQtspZnrYKXj4OjrwJxtwBMRXfUykiIiIiUpeENSEEJhO6QdiNu6+/AWFJCK21Q2uw+Yk4NbKDjDHT/cr6uY93GmMuB76y1v6zBscKOf8moykV9DAaTHSU4aHTD6JLSiIPf7WCYgv5RcU8+PlyPlmwjYfOOIjBnVqEIOJa0HG404T0+wdh9rOABVsEPz0Kq6bBGc9Bu0EV7kZEREREpK4Id0J4SS3sc0Et7LM2DS2nrJ87bQhLJFUQqMlodRljuGpMT4Z2bsF1b/9O2n6n9ezirRmc+sxMxg9oy5+P6c1BHZMr2FMExCY4g9n3O8kZp3DvBmf5zsXwwtEw9g444iaIDvdbS0RERESk6sL6rdVa+1o4j1eXWGvvA+4LVGaMmQxcDPzRWvtm+KKqPP+EMCWp6k1G/Y3s0YppN47moc+X87/ft5Ys/3rZTr5etpMRXVty3qFdmHBQO5Li61iC1XUUXDUTvr0XfnvJWVZc4NQervgCznjeuf9QRERERKQOi0inMsaYEyNxXKm+BL8moy1rUEPoq1XTeB4/dyhvX34Y/ds3L1U2d+NebpmykGEPfMMVr8/l3Tmb2LgnG2tD1eq4huKbwkmPOUNUNPcOq8G2+fD8aJj1DBQXRSw8EREREZGKRKra5SFjTBPgC2ttbkUrG2OirLV1cCTzxiMxhE1GAxnVqzWfX38k3yzfyTPfr2Hx1oySsrzCYr5ZtpNvlu0EoH1yE0b2aMXIHikM7dySnm2SiImOYIe5PY92BrP/6q+wwK3gLcqDr++EFZ/B6f+BlB6Ri09EREREJIhIJYQHgClArtvByhc4yeH6IOsfb4xJstZ+EK4ApbSEuNKXSiiajPqLijIcP7Ad4we0ZfHWDN7+dROfL95OVm7pAey3Z+Ty0e9b+chtZtokNooB7ZszuFMLDuqYzKCOyeFPEpskw+nPQv9TYOqfYb+TvLJptjOY/XH3w4jLnPENRURERETqCBOJ5nfGmMOB14Be7iJPEKtwk0NghrW20F0/BvjRWjsq3LHWR8aYpQMGDBiwdOnSkO1zzvp0zvnv7JL5//5xOMcPbBey/QeTX1jM7HV7+GbZDn5Zl86aXfsrtZ1vkjioYzJDOregR+ukKo2dWG056fDFrbDkw9LLe4yFU5+BFp1rPwYRERERqXMGDhzIsmXLlgUbPi4SIlJDaK2dDfQxxgwGTnCnUThj/PUBbgSyjTHfAbOAFpTfO6fUMv8moylJ4RlzLy4mijF92jCmTxsAdmXlMmd9Or+s28PcDXtZvWs/RcVlf9TILShm/qZ9zN+0r2RZs/gYBndOZnCnFgzp1IKhnVvQLrlJ6INOTIGzX3FqCz+7GQ6kO8vXTYfnRsGEf8DQC8DUwTEXRURERKRRiUgNYSDGmGTgR5xhJE4AWrtFngA3WGt7RiC0eqc2agjX7NrPsY/PKJn/9uaj6JXaLGT7r67cgiKWbc9kydYMFm3JYMnWjKBJYiBtm8cz2E0OD+uewqBOycTHRFe8YWXt3wVTb4SVn5de3mcCnPJvaFb7tawiIiIiUjeohrAc1toMY0ystfZiY4wBRgIn4ySH7YG/RjTARq62O5Wpriax0Qzr0pJhXVqWLPNPEhdt2cfqXfsJ9NvHzsy8Uh3WxMdElSSHh3RPYViXljUb8qJpKvzhLVj0HnxxO+S5neWs+gr+MxJOfBQOOku1hSIiIiISEXUmIfRlnWrL2e50Z4TDEaBNs3jaNo9nZ2Ye3VolklJHEsJAAiWJ+/MKWbwlg4Vb9rFw8z4Wbclg674DZbbNKyzm1/Xp/LreaeYZHWU4qGMyh3ZryaHdW3FYjxSaN6lihzrGwJA/QLfR8Ol1sPZ7Z/mBvfDhZbB8Kpz0OCS1qvY5i4iIiIhUR51pMgpgjFlmrR0Q6Tjqu9poMgqwbFsmXy3ZzslDOtCnbeSbi9bUrqxcFm3OYMHmfczbuJf5m/aSV1j+6CZRBoZ0bsHoXq05sncbDu7Sgtiq9GZqLcybDNPuhIJs7/KkNnDyk9D/5Gqdi4iIiIjUfXWxyWhdSwiXW2v7RzqO+q62EsKGLr+wmMVb9zFn/V7mrHc6rcnKKyx3m6S4aEb2aMWRvVszundrerZpiqlM88/09fDJdbDx59LLB/8BTvgnJLQMvJ2IiIiI1FtKCD0HNeY54HxgK05HMl8B3wE/WGtHhD2gBkYJYWgUFVtW7Mhkzvp0ftuQzuy1e9ibU1DuNu2aN+GoPq0Z1y+VI3u3oWl59x8WF8Oc/8K390Fhrnd5sw5w2tPQ69jQnIiIiIiI1AlKCD0HNSYH8PT37wmgEMgAHgA+t9auC3tgDYQSwtpRXGxZtj2Tn1an8fOa3fy2YS/55TQxjY02HNo9haP7pjKuXyo92jQNvGLaavj4atjyW+nlwyfB+Achvv43zxURERERJYTegxqzBfgfsBs4AjgSSHSLPQGtBD4HPrXW/hT2IOsxJYThcSC/iN82pPPzmjR+Wp3G8u2Z5a7frVUiY93k8LAeKaWHtygqhFlPwQ9/h2KfWsgWXeD056DbkbV0FiIiIiISLkoIPQc1Zhpwh7V2gTsfD4zBO0h9H3dVT3AzgOOtteW31xNACWGk7M7K4+c1u/lhxW5mrNpNxoHgl2tSXDRj+6YyfmBbju6X6u25dOdS+Ogq2LHIu7KJgrF/hdG3QFQVOrARERERkTpFCaHnoMZc6h775SDl3fEmh+NwmpdeaK19J3xR1l9KCCOvsKiY3zfv4/sVu/hhxS5W7MgKum5stGFkj1YcP7Adxw1oS9vEKPjpUfjxUbBF3hV7joMzX4Sk1mE4AxEREREJNSWEnoMaEwe8B5xpKwjAGNMUuAvItdbeF4bw6j0lhHXP1n0H+GHFLqav3MXPa9LILQh+7+HBXVowfkA7Tm29nY7fXgP7NnoLm7WHs1+BrqPCELWIiIiIhJISQt8DG3MoMMpa+2Ql159orZ1Su1E1DEoI67bcgiJ+Xp3G18t28O3yXaRn5wddd0gbeCTmBfrsne5daKLhuPvh8GudQe9FREREpF5QQuh/cGOirLXljwQuVaaEsP4oLCpm3sa9TFu6k6+X7WDL3gMB1rJcEv0Vf419h1h8xkUcegGc/ATExIctXhERERGpPiWEEhZKCOsnay3Lt2cxbekOvl62s0yvpUPMGp6Ne4pOJq1kWXbb4SRe+A6mWdtwhysiIiIiVaSEUMJCCWHDsGlPDl8v28G0pTuYu3Ev1kIKmTwX9ySHRa0oWW+nac0XBz3OISPHMrBDc4yakYqIiIjUSY0+ITTGXFQLu11grV1U8WqNhxLChmdHRi5fLtnO54u2s3Djbu6PeZXzYn4oKc+x8VxdcCMbWx7OiYPac+rQDvRr1zyCEYuIiIiIPyWExhTjHVuwxrtz9/U3a+39Idpng6CEsGHbnnGALxdthzn/5eLM/xJtnLdUvo3mhoLr+LL4MAD6tWvGqUM7cOqQDnRqmRjJkEVEREQEJYQYYy6uhd0usNYurIX91ltKCBuPPYu+pPknlxJblANAkTX8pfAKphSNLbXeId1acurQjpw0qD0pSXERiFREREREGn1CKOGhhLCR2TIP3jwTcveVLPpbwR95teiEMqvGRBmO6tOG04Z24LgBbUmMiwljoLXAWsjZA3s3QOY2yEmD7D3OY84eyM+GggNQmOs+5oGJcicDUdEQHQdxTSG+GcQ3hfjmkJgCTdtB07bQrK0z/mNSGw3zISIiIjVSFxPCev5tUEToNBwu+RLeOB327wTg3tg3OK5nEneln8C6tOySVQuLLd+v2MX3K3aREBvN+IFtOX1oR0b3bk1MdFSETqASrHWSvh2LYMdi2LnMmd+3EfL3hyeG2ERI6eGdWveGtgdBan8N/SEiIiL1lmoIGyDVEDZS6evg9dNg36aSRfbkJ1nS7kw+XrCVqQu3sSsrL+CmrZvGc/rQDpw1vBP929eBzmgK82HrXFj/I2z4GbYvhLzMireLhKgYaNMP2g2GjsOg6yho0x+i6nCCLSIiIhFRF2sI62RCaIzpBvQFWrqL9gIrrbUbIhVTfaKEsBHL2OrUFKatcuZNFJz3LvQ5nqJiyy/r9vDJgq18uXgHWXmFAXcxoH1zzhreidOGdqB10zDWfGVshWUfw5pvYdMvUJBT+W3jm0NyZ0hq7UyJ7mN8M6f2LiYBYptAtHs+tghssTMV5kFellPTmJflTNm7IWsn7N8B+3dVvRaySQvoMtJJDnuOc2oS1dxURESk0VNCWAFjzAXA3UBvnF5EfVlgNXC/tfbtcMdWnyghbOQytsBLx0LWdmc+NhEmfe7UXrlyC4qYvnIXH/2+le9X7KKgqOznQEyUYWzfNpw1rBPj+qcSHxMd+liz05wkcPGHsGlWxeu37A7tBjm1ca16QstuzpTQsnYTrtwMpwY2fR3sWQfpa2HXcti1DIryK96+eUfofRz0Ph56jIG4pNqLVUREROosJYTlMMZcCzwJvABMAVbg1AwCtAD6AROBK4CbrLX/CX+U9YMSQmHHYnjlBMjPcuaT2sDl3zrJk5+92flMXbSND+dtYeGWjIC7S06I5dQhTpPSIZ2SMTVNvrbOg5+fhBWfO7V1gcQ0cWrZuh8FXQ53atma1IHmrL6KCmD3Suf53r4ANs12/rbFwbeJaQJ9JsCgs6HXcU7NpYiIiDQKSgjLYYxZBzxrrX2sgvVuBa611nYPT2T1jxJCAWDNd/D2OVDsNg1t1Rsu+9rpQTOI1Tuz+HD+Vj76fQs7MwPfb9izTRJnDuvEmcM60j45ofLxWAtrv4eZTzr3BgbSsjsMPAN6HQOdDqmfnbXkZsKWObBxFqybDlvnE3T41fhk6H8KDDkXuo1Ws1IREZEGTglhOYwxOcBJ1tofKlhvHPCZtVYjbQehhFBK/P4WfHKNd7738XD+exUmHkXFlplr0vhw/hamLd1BbkHZGi9j4IierTlreEeOH9iu/CEsNvwMX/2f00uov2Yd4KAznanDsIaXFO3fDWu+gVXTnIQ4WOc4rfvAiMtgyB8goUVYQxQREZHwUEJYDmPMTzhNRCdaawNWTRhj4nGak7aw1h4VzvjqEyWEUsr0f8L0f3jnz3wRBp9T6c2zcgv4cvEOPpi/hTnr0wOukxQXzQmD2nPWsE4c1j2FqCg3qcvNgG/uhXmvlt2o0yFwxI3Q98TG0yNnYb6TFC6eAiu/CNxxTmwiHHQWjLwG2g4If4wiIiJSa5QQlsMYMxT4GigGvsC5h3CfW5yMcw/hiUAUcJy1NkBVg4ASQvFjLbxxBqxzK98TUuDaOdC0TZV3tWlPDv/7fQv/m7+VTemBewHt2CKBM4d15IKWy2j3418ha1vpFXodB0fe5PTA2dBqA6siPxtWfgkL33Ga9wZqVtr3JBh9izPWpIiIiNR7SggrYIxpDVwNjAf6UHrYiVXANOB5a21aZCKsH5QQShl7N8J/DocCd5D6gWfAxMnV3p21lt827OXDeVv4fPF29vsMYdGEPP4R+xJnRM8svVHqQDj1aSU3gaSvd2pR578BBwLUwvYY6ySG3dUwQkREpD5TQihhoYRQAvr1v/Dl7d75c9+C/ifXeLcH8ov4etkO/jd/K/NXb+LF2EcZGbW8pDzPxvCf4jNZ1/cyTh/enaP6tCE2upE0Ea2qglxnGI5ZT8POJWXLe46D4+53ht4QERGRekcJoYSFEkIJqLgYXj0BNv/izDdtC9f+6ozhFwoH9lHw+lnEbp9bsmhucR/uKLiCtbZjybJWSXGcMKgdpwzuwCHdfO43FC9rnU5ofnoUtvzmV2icjmfG3QXJnSISnoiIiFRPXUwI693P9MaYCe4QFSJSFVFRcNozEO0O5bB/J0y7KzT7zkmH108tlQzu7TORL0e8wr7E0iPE7MnO581fNnHuC78w6p/f88Bny1iweR/6ccqHMdB3Alz2DVw8Fboe6VNonfsOnxoG3/4N8gPfyykiIiJSGfWuhtAYcxbwvrU2OtKx1FWqIZRy/fwEfHufd/6Sr6Dr4dXf3/5d8PrpsMvnehtxGZz4KERFUVBUzI+rdvPh/C18u2wX+UWBB23vnJLASYM6cNyAthzcuYVqDn1ZC6u/hm/ugd0rSpe16AInPQG9j41MbCIiIlJpdbGGsM4khMaY+yu5an/gTCWEwSkhlHIVFcJLx8D2Bc78QWfD2S9Xb1+F+fDiONi52Lts5LVw/EMBexDNOFDA10t3MHXRdmauSaOoOPDnT+um8RzbP5Vj+7flyN6taRKrtzvgvHYL3oQf/u7U8Po66Cw4/h/QrG1kYhMREZEKKSEshzGmGGeYiYwKVk0EWishDE4JoVRo+VR470Ln7+h4uHVl9e4lnPU0fO3T7HT0LTDu7koNJ7Fnfx5fLtnB1IXbmLMhnWAfRU1ioxjduw3H9k/lyN5t6NgioepxNjR5++HHh2HWM2CLvMvjk2HC32HoBY17SA8REZE6SglhOYwxq4GfrLWXVrDe2cB7SgiDU0IoFSrMh8f7Q447gsuJj8KhV1RtH1k74OnhkL/fmR92kTOsRDXszMzl80Xb+WbZTuZsSA9acwjQvXUSR/RqxZG9WnNY91a0TIqr1jEbhB1LYOoNsHVu6eX9T4VT/g2JKZGJS0RERAJSQlgOY8zrwKHW2n4VrHcWMMVaW+86xAkXJYRSKdPuhNnPOH+3HwJ/+rFq2394BSx+3/k7oSVcPz8kCci+nHymr9zNN8t2MmPV7lJjHAbSs00SI7qmMLxbS4Z3bUn3VkmN6/7D4iJnDMNv/wZ5md7lzTrAGc9DjzGRi01ERERKUUJYDmPMcOAIa+1TFazXGhhorZ0RnsjqHyWEUim7lsN/Rnrn//QTtB9cuW03zITJJ3rnT34CRpRbuV8teYVF/LIunW+X7eTnNWmsT8uucJum8TEM6NCcQR2TOahjc/q3b063VkkN/z7EjK3w0Z9gw0+ll4+6HsbdAzGNuCZVRESkjlBCKGGhhFAq7cVjvM0ND/0TnPhwxdsUFcILY7wDp7cfCld8D1G1n3Bt2ZvDrDV7+HlNGrPW7iFtf16ltosy0DklkV5tmtIztSk92yTRK7Up3Vs3pWViLKah3G9XXAyzn4bvHoDiAu/yziPhnNfV4YyIiEiEKSGUsFBCKJU2b7JzDxpAkxZwy0qIbVL+Nr++AF/e5p2/7FvofEhtRRiUtZZN6TnM3bCXuRvTmbthL6t37a/yfpLiounUMpFOLRPo1DKBzimevxPp3DKR5MTYWoi+lm1f6DTpTVvpXdasPZzzRkReKxEREXEoIZSwUEIolZabCY/1hQJ3cPOzX3GGLwhm/254Zjjkup0BD70QTn+29uOspKzcApZvz2LJ1gyWbMtg6dZM1qXtp6Co+p9zzZrElCSMXVISSyZP4lhnm6Lm58AXtznDVHhExcJJj8LwSRELS0REpDGriwlhTCQPboy5qJKr5gN7gIXW2l21GJJI49KkOQw4HRa+7cz//mb5CeHPj3uTwfhkOPa+2o6wSpo1ieXQ7ikc2t3buU1BUTGb03NYuzubNbv2s3b3fudx136yKuiwBiArt5Dl2zNZvj0zYHm75k1KEsQuKYl0a51Ir9Sm9GzTNLLJYlwinPYMdDwYvrwDigudZqRTb4Btvzs9y0bXw9pPERERCamIJoTAZKAqP91bY8y3wPXW2tW1E5JII3Pwhd6EcO0PsG8TtOgSeN31Pj2RHnUrNG1T+/HVUGx0FD3aNKVHm6YcN8B7D521ln05BWzZe4Ate3PYvDfH/fsAm9Odvw8UFJWzZ8eOzFx2ZOYyZ0N6qeXGQOeWifRObUovd+rbrhl92jYLX6JoDBxyOaQOhPcvgmz397R5kyFjC0x8DeKbhicWERERqZMinRDeD3QDLgL2A18Dm9yyzsB4oBnwBpAHjHKX/WSMGW6t3RrugEUanK6jIKUHpK8DLCx4B8beUXa9wnzY7XNPWs9xYQuxNhhjaJkUR8ukOAZ1Si5Tbq1lT3Z+ScK4KT2HzelOsrgpPYet+w6UO16itbDJXfe7Fd6GDTFRhj5tmzGwQ3MO8ukJNTGuFj+Oux4Of5rhJIVbfnOWrfkWJp8EF0yBpqm1d2wRERGp0yJ6D6ExpicwB/gIuMVam+FX3hx4HDgDOAxYBzwC3AQ8a629PrwR1w+6h1Cq7KfH4Lv7nb9bdIEbFjm1S752LoXnRjl/R8XAX7c36qEMCouK2Z6RW5L0bUrPYdOeHNbu3s+63dnkFxVXel9RBnqlNmV415YM6+KOp9g6KfS9nxbkwkdXwrJPvMtadoML/weteob2WCIiIlJGXbyHMNIJ4fvAMKCPtTbgtydjTBSwCphvrT3HGBMHrAdyrLW9wxdt/aGEUKoscxs8PoCSFtzXzYPWvUqvs2gK/O9y5+/UgXDNrLCGWJ8UFhWzee8BVu/MYs3u/azZuZ/Vu/azcmcW+YWVSxRbJsY6CWLXlozomsKQzsnEx4SgqWlxMUz7P/j1ee+yxFZw/hToNLzm+xcREZGg6mJCGOkmo0cDXwdLBgGstcXGmDk4TUWx1uYbYxYCY8MTokgj0LwDtO4Naauc+W2/l00IPeMOArStM59hdVJMdBTdWyfRvXWS88HlKigqZs2u/SzZmsHSbZks2ZrBsu2Z5OSXvVdxb04B3y7fxbfLneamTWKjGN61JSO7t2Jkz1YM6dSCuJioqgcXFQUT/gnNO8I3dzvLcvbA66fBhR9Al5HVOGMRERGpryKdECYC7SqxXlvAd3C0TKDi7gHrMGPM6cBVwHCc+yR3A3OBx6y1P0cwNGmsOgwrnRAOnli6fKdPjbMSwmqJjY6if3vnnkHPs1tUbFmftp/5m/Yxf+Ne5m0MPJ5ibkExM9fsYeaaPfCNkyCO6JrCyB4pHN6zFYM6ViFBNAaO+LMzNuHHVzu9j+ZnwRtnOvcUdjsidCctIiIidVqkE8LFwFHGmKOstT8GWsEYMxoYA/zms7gzTgJV77hNYF8ELgWygZ+BfUAX4ERgnrtMJLw6HAyL3nX+3vZ72XIlhLUiOsrQK7UZvVKbcc6IzgBk5BTw++a9zN+4l7lukpjn19Q0t6CYn9ek8fOaNAASYqM5rEcKR/ZqzejebejTtmnF9yAOnggJLeHd86EoDwqy4a2z4fz3ofvoWjlfERERqVsinRA+DHwATDPGvO7+vdkt6wychdMDqXHXxRiTjFOr9mHYow2Ne3CSwanAJGttSV/1xpiWQOtIBSaNXIeDvX9vXwjFRRDl3rOWkw5Z27zlSghrVXJiLGP7pjK2r9P7Z15hEYu2ZDB77R5+WbcnYIJ4oKCI6St3M33lbmA5bZrFc2Sv1s7UuzVtmzcJcCSg97Fw3jtOUliYCwU58NZEOP9d6DG2dk9UREREIi6incoAGGNuAP4FxFF2TEKDMyj97dbap9z1ewATge+stXPDGWtNGWM6AWuBHUA/a+2BWjqOOpWRqsvPgX90Auvez3bNL5Da3/l7/U/w2snO3wkt4fb1ZXshlbDJKyxiwaZ9/LIundnr0pi/aV+FndX0aduUI3q1ZnTv1hzWvRVJ8X6/B66bDm//AQrdj6WYJk5NYY8xtXMSIiIijZA6lQnAWvtvY8ynwGU44wy2d4u2AzOBV62163zWX4eTQNZHF+Mkvi/VVjIoUm1xiU4C6Ok8Ztvv3oRw1zLvem0PUjIYYfEx0RzWoxWH9WjFDfQmt6CI+Rv38tOaNH5encaSbRn4/9a3aud+Vu3cz6szNxAbbTi4S0tG92rNEb1bM7hjMjE9xsIF78Pb5zq1hIW5Tq3hxZ9CR/U+KiIi0lBFPCEEsNauB+6KdBxh4BnJe5Yxpj1wAdALyAB+AKbZSFfZSuPWYWjphHDo+c7f6mG0TmsSG82oXq0Z1as1d0yA9Ox8Zq1NY+aaNH5ancaWvaV/fyoossxZn86c9ek89s0qmjWJYVTPVhzVpyvjT3uLNp9c4CSF+fvhzbPh0q+gTd8InZ2IiIjUpjqREDYiA3wePwSSfcpuB6YbY86w1u4Ld2AigHMf4e9vOn9vne9d7tuhTOoApG5LSYrj5MEdOHlwB6y1bNyTw09r0pi5Oo1Za9PIzC3dSXNWbiHTlu5k2tKd3Amcl3IbDxY+RLQthAPp8PrpcNk0aNElIucjIiIitafOJITGmMOB0UBHd9FW4Cdr7ezIRRVyLd3Hx4HZwA3AGuBQnJ5Hx7qPEwNt7M8YE+wmwZ41ilIaL9+OZXYshqICMFGwa7l3eduDwh+XVJsxhm6tk+jWOok/juxKYVExi7dm8PPqNH5ak8bvm/ZSUFS6YcI76X3Jirqap2KfIcpYyNpGxgsns//8z+jYSUmhiIhIQxLxhNAY0wd4AxjhWeQ+Wrd8LnChtXZ1BMIrxRjzEdC/iptdZK2d4/7tGSRsL3CCtTbbnf/OGHMqsAg42xjTx1q7quYRi1RR24MgKtYZl64oz0kE45Kc5oMAGEjtF9EQpWZioqM4uEtLDu7SkuuP6U12XiFz1qczY9Vufly1m3VpzsfSZ8WHk1yYzUOxrwCQnLORTS+cwpXN/8GhfbtydN9UDuuRQnxMdCRPR0RERGooogmhex/dDJyB57cBU4ANOMlgN5yaskNwmlKOsNZuj0ykJboDVb2RJtHn7/04tYRTfJJBAKy1S4wxv+HUFh4FVJgQBuudyK05VLs+qbqYeOcewe0LnPltv0NCC295Sg8nQZQGIyk+hqP7pXJ0P2eIi017cpi+ahfTV+7mw7XjaVGwn9ti3wdgUNQGbsn4F1fMvIVXZ26gaXwMY/q2YfyAtoztm0pyQmwkT0VERESqIdI1hHfhJINPAP9nrc33LTTG3AH8A7gZ+Ctwfdgj9GGtHVrDXWzESQg3BCnfgJMQptbwOCLV1+Fgn4RwPjRt6y1ThzINXpdWiVx0eDcuOrwbuQVF/LZ+OL98axi56z0AxkUv4F77OvcUTmJ/XiGfL9rO54u2ExNlGNmjFccNaMtxA9rSoUVChM9EREREKiOq4lVq1YnASmvtLf7JIIC1tgC4DVgJnBzu4GrB7+5jyyDlKe7j/jDEIhKY732E234v3aGMEsJGpUlsNKP7pDLyquehn/cj+KKYb/hT3Fel1i0stvy8Jo17P13KqH9+z0lP/cS/v13Nqp1Z4Q5bREREqiDSCWF7YH55K7jDMMzHOz5hffap+1hmpGdjTFNgmDv7u3+5SNj4JoQ7lzpJoYcSwsYpKgrOfBE6DCtZ9JeoN/n4mL2cM6ITrZLiymyydFsmT3y7ivFP/Mj4J2bw1HerWbdbv3WJiIjUNZFOCDOBzpVYr7O7bn03FVgOjDLGXONZaIyJxul5NAVYAvwcmfBEcAajj2ni/F1cCJlbvWVKCBuvuEQ4/z1IdnoZNViG/norDx9eyJw7j2XKVYdz5VE96N667D2mq3bu5/FvVjHusRmc+O+feG76Wjan55RZT0RERMLPRHIcdGPMJzhNQU+11n4eZJ0TcRKpqdba08MYXq0wxgzF6UinObAQZ9iJg4EewB7gaGvt4hoeY+mAAQMGLF0abFQKkQq8dCxs+a30stgk+L8tTm2RNF67VsDL4yEvw5lv2hau+AGSnRGDrLWs3b2faUt38tWSHSzemhF0V0M6t+CUwe05aXB72ifrnkMREWn4Bg4cyLJly5YF6xwyEiKdEB4O/IjTq+h7wNt4O1zpCpwH/AGnJnO0tfaXCIQZcsaY7sB9wHigFbAT+Ap40Fq7MQT7V0IoNfPFbTDnhdLLOh0Cl38bmXikblk3Hd48y6lBBmg3GC79KmAPtBvSsvl88XamLtzGih2B7yc0Bo7o2Zozh3Xk+IHtSIqPdH9nIiIitUMJYaAAjLkQ+C+QgDv2oG8xcAD4k7X2zXDHVl8pIZQaW/A2fHx16WXDLoZTn4pMPFL3zHsNpv7ZO9//FJj4erk1yGt27eezRduYunAba3dnB1wnMS6aEw5qz1nDOjKyRyuiokzA9UREROqjupgQRvxnWGvtm8aY6cAVwJFAB7doG/AT8LK1dnOEwhNpnHw6DynR9qDwxyF11/CLYfdK+OVZZ375VPjhITjm7qCb9Eptyo3H9uGGY3qzYkeWmxxuZ5PP/YQ5+UV8OH8LH87fQofkJpwxrCNnDutEzzZNa/uMREREGqWI1xBK6KmGUGqsuAj+0RkKfGpxJn0B3Y6IXExS9xQXwTt/gNVfe5ed+SIMPqfSu7DWMm/jXj6cv5XPFm0jK7cw4HpDOrfg3BGdOXVoB5qqSamIiNRTdbGGUAlhA6SEUELilRNg0yzv/B0bICHYEJrSaOVmwsvHwe4Vznx0PEz6DDofWvVdFRTx7fKd/G/+Vmas2k1Rcdn/T4lx0ZwyuAN/OLQzQzu3wBg1KRURkfqjLiaEYf2Z1RizrgabW2ttz5AFIyLl63CwNyFs3lHJoATWpDmc9y68OA4OpENRHrx7vtPzaIvKjCrks6vYaE4e3IGTB3dgd1YenyzYyv/mb2XZdu+oQzn5Rbw3dzPvzd1Mv3bN+MMhnTnj4E4kJ8aG+sxEREQahbDWEBpjimuyvbVW/d1XgmoIJSRWTYO33aZ/B50NZ78c2XikbtswE14/DYoLnPm2g5yeR+Nrfu/fsm2ZvD93M/+bv4XMAE1K42OiOHFQe/5wSGcO7Z6iWkMREamz6mINoZqMNkBKCCUkrIWZ/4b0dTDmdkjuFOmIpK6b/wZ8ep13vu9JcO6bIRu7MregiC+XbOedOZuZsz494Dq9Upty0eFdOXNYJ91rKCIidY4SQgkLJYQiEjFf3wWznvbOH3kTHHtfyA+zdvd+3vttMx/O28Ke7Pwy5U3jYzhzWEcuOrwrvVKbhfz4IiIi1aGEUMJCCaGIRExxkXMP4aqvvMtOfx6Gnlcrh8svLOabZTt597dN/LQ6LeA6o3q24qLDu3Fs/1RionXngYiIRI4SQgkLJYQiElF5WfDyeNi1zJmPioXz3oHex9XqYdenZfPmLxuZMndzwHsNOyQ34YKRXTn3kM60bhpfq7GIiIgEooRQwkIJoYhE3N6NTs+jOW6tXXQ8nP8e9Dy61g+dk1/IJwu28frsjSz36aHUIy4mijOGduSy0d3p01bNSUVEJHyUEEpYKCEUkTphyzyn59H8LGc+JgEu/AC6HRmWw1trmbtxL6/P3siXi7dTGGBcw6P6tOHyI7szundr9U4qIiK1TgmhhIUSQhGpMzbOhjfPgoJsZz42Cf74EXQ5LKxh7MrM5Z05m3nz143szsorU96nbVMuP7IHpw7tQJPY6LDGJiIijYcSQgkLJYQiUqes/xHemgiFuc58fHO48EPofGjYQ8kvLGbqwm289PP6gM1JWzeN448ju3HhyC600n2GIiISYkoIJSyUEIpInbPmO3jnD1DkDhERFQvH/Q1GXgMRaKpprWX22j289PN6vl+xq0x5fEwUE0d04k9H9aRzSmLY4xMRkYZJCaGEhRJCEamTVn4F710IxQXeZb3Hw+nPQVLriIW1Ztd+Xpm5ng/nbSGvsLhUWXSU4ZTB7bl6bC/6tlMHNCIiUjNKCCUslBCKSJ21/if43xWQtd27rGk7OPMF6DEmcnEB6dn5vPXLRl6bvZG0/WXvMzymXyrXHN2T4V1TIhCdiIg0BEoIJSyUEIpInZa9Bz65pvTg9QCdDoGh58PAMyGhRURCA8gtKOLD+Vv474x1bErPKVN+aPcUrhnbkzF92qhnUhERqRIlhBIWSghFpM6zFn59Hr65x3tfoUd0PPQ7CbqPhmYdoHkHaN4RElPCer9hYVExny/eznPT17JiR1aZ8gHtm3PduF5MGNiOqCglhiIiUjElhBIWSghFpN7YvhCm3gDbfq94XRMNUTEQFe3+HQUmyv072luemALJnbxTiy7QcQQkd6xWiNZapq/czX+mr+G3DXvLlPdr14wbjunN8UoMRUSkAkoIJSyUEIpIvbN9ESx8Bxa9Bzl7aucYLbpAl8OdqccYSOlR5V38tiGd56avDdgzqRJDERGpiBJCCQslhCJSbxXmw5pvYMUXsG+j0/lM5jYoKHsvX411HA5DznPuWUxqVaVNl23L5JkfVvPF4h1lyvq1a8aNx/Zm/AAlhiIiUpoSQgkLJYQi0qBYC7kZkJ0GtgiKi8AW+/1d7PxdlA/7d0HGZsjcChlbYdcyJ7kMJirGGf5i+CTnsQr3Ka7Ykcm/v13Nl0vKJob92zfnhmN6M35AWyWGIiICKCGUMFFCKCLiJ3MbbJwFm2Y7Q1+krQy8XttBMPpmGHCac19iJS3fnslT3wVODAe0b85tx/dlbF/1Sioi0tgpIZSwUEIoIlKBHUtg0buwaArsL5vE0ao3HHkTDD4HomMrvdtl25zE8KulZfd5aLcUbp/QlxHdNI6hiEhjpYRQwkIJoYhIJRUXwbrpMOdFWPVl2fLWfeH4v0PvY6u02/ISw3H9Urnt+L70b9+8mkGLiEh9pYRQwkIJoYhINexYAj89Bks/Avz+N/Y6FsY/BKn9qrTLxVsyeHjaCn5anVZquTFw6pAO3HxcH7q2Sqph4CIiUl8oIZSwUEIoIlIDaWucxHDhO5RKDE00HHIZHH0nJLSo0i5nrU3j4a9WsmDzvlLLY6IM/9/efYe3WZ19HP8eSd4rie3svfdOCAkZ7L0Je5VRVtmlLW3fQoG2FChllw0ttEDYm7AChCSQvffeyxneU+f945E1bDnTlmzr97muXNIzfazI0nM/55z7vnB4O245phvN0xMPu+kiIlK/KSCUiFBAKCJSC7bMh4m/h7WTQ9enNIeT/gZ9zz2ojKTWWr5avI2HJy5jxfb8kG2JcS5+OaYL143pTEqCpzZaLyIi9VB9DAhd0W6AiIhIvdSqP1zxMVzwX2jaKbC+YDu8ezW8djbkrDrg0xljOKFPS764bQz/GD+Atk2T/NuKy7w88c0Kjn7kOybM3ECFVzdrRUQkMhQQioiI1MQY6HUa3PQzHHsPeAJBHKsnwTNHwvcPQ3nJAZ/S7TKcO6Qt39w5lntP701Warx/2/a8En7zznxOf/JHpq7cuY+ziIiI1A4FhCIiIvvjSXDqE970E3Q9PrC+ogQmPQDPHgVrpxzUKRM8bq4c1YlJvx7HTUd3Id4T+EpevCWXi1/8mWv+PZNVO/L3cRYREZHDo4BQRETkQDXtCJe8DeNfhdQWgfU7l8Orp8AHN0FBzkGdMi0xjrtO7Mm3d47lzIGtQ7Z9vWQbJ/7zB+79aBG7C0oPv/0iIiJVKCAUERE5GMZAn7PhVzNg2LVAUGKZua/DU0Nh1qtOjcOD0LZpMo9fOIj3bhzJ4PZN/OvLvZZXp65l7MOTePnHNZRXeGvl1xAREQEFhCIiIocmMQNOfQSu+QZa9gusL9oFH98Kz4+DddMO+rSD2zfl3RtG8tTFg0ISz+QWl3PfJ4s57ckf+Xn1wfVCioiI1EQBoYiIyOFoOwSu/Q5O/CvEBRWZ3zofXjkJ3rkK9m48qFMaYzitf2u+vmMsd5/ck7SgUhRLt+ZxwfM/cdubc9iWW1xLv4SIiMQqBYQiIiKHy+2BI29yhpH2Gx+6beG78OQQ+OL3kL/9oE6bGOfmurFdmHTXOM4f2jZk2wdzN3PMI9/xwg+rKdMwUhEROUQKCEVERGpLRhs490W4aiK0GhhYX14MPz0Nj/WHiX+A/B0Hddqs1AQeOm8A7904kr5t0v3rC0or+MtnSzjl8clMW6VhpCIicvAUEIqIiNS29iPg2klwxlOh2UjLi2DaU/B4f/ji7oMqbA/O/MIPbzqKv5zdlybJcf71K7bnc9ELP3HHhLnk5B94TUQRERFjrY12G6SWGWMW9e7du/eiRYui3RQRESkthJkvw5THoKBqz6CBbsfD8F9Cl2PBdeD3aXcXlPLQxGW8OWM9wV/lGUlx3H1yT84f2g6Xy9R8AhERibg+ffqwePHixdbaPtFuSyUFhI2QAkIRkXqotABmvARTHofCndW3N+sCAy9y5iA27XjAp523YQ9/+GABCzflhqwf2qEpfzm7Hz1aph1mw0VEpLYoIJSIUEAoIlKPleTDvDdg+vNOQftw2o2A/uOhzzmQ3Gy/p6zwWv4zbS3/+HI5+SXl/vUel+Hq0Z249dhuJMd79nEGERGJBAWEEhEKCEVEGgBrYfV3MP0FWP452DCZQl0e6Hq8Exx2Pxnik/d5yq17i7n/k8V8umBLyPo2TZK478w+HNurRQ1HiohIJCgglIhQQCgi0sDsXgfzJ8CCCTX3GsanQq/Tof/50GksuNw1nm7S0u386aOFbNhVFLL+pD4tueeM3rTKSKrhSBERqUsKCCUiFBCKiDRQ1sKWeU5wuPAdyN8Wfr/UltDvPOh/AbTsB6Z68pii0gqe/HYFz/+wmnJv4Ls+Jd7N7cd35xejOuFW0hkRkYhSQCgRoYBQRKQR8FbAmh+c4HDJx1CaF36/7F4w+DIYcFHY+YbLt+Xxx/cXMn3trpD1A9o14eHz+tO9hZLOiIhEigJCiQgFhCIijUxZESz73AkOV34F3vLq+7jjofeZMORK6DAqpNfQWsvbszbyt8+WsLuwzL8+zm245ZhuXD+uC3FulSYWEalrCgglIhQQiog0YgU5sOg9JzjcOD38Plk9YMQNMOBCiAvMF9xVUMoDny7mvdmbQnbv3Sqdh87rT982GXXZchGRmKeAUCJCAaGISIzIWQVzXnf+FWyvvj05E4ZdA8OuhdRs/+pJy7bz+/cWsGVvsX+dx2W4YVwXfnVMVxI8NSesERGRQ1cfA0KND4kwY0yCMea3xpjZxph8Y0yJMWaNMeYFY0znaLdPREQakMwucNw9cMdiOP8/0OWY0O2FOfD93+GffeCT22HPBgCO7tGcL28fw8VHtPfvWu61PPntSk574kfmrN8dyd9CRESiSD2EEWSMSQQmASOAPcBUoBgYDHQE8oCjrbWzDvPnqIdQRCRW5ayCn56BOf+F8tCyE7jiYPDlMPoOyGgLwJSVO/ntu/PZuDuwr8vANaM7c8fx3UmMU2+hiEhtUQ+h/BInGJwBdLTWnmqtPRfoCjwFpAGPRrF9IiLS0GV2gVP/4fQaHvNHSA0qRu8tg5kvweMD4ZM7IHczo7pmMfG2MVw5sqM/D43XwvM/rObkxyczo0p2UhERaVwUEEbWGN/jo9bavZUrrbUVwJ98i8Mi3ioREWl8kpvBmLvgtgVw6qOQ3iawrTIwfGIwfHM/KbaQe8/ow4TrjqRTVop/tzU7Czj/uWnc+9EiCkrCZDYVEZEGTwFhZJUcwD45dd4KERGJHZ4EGHY13DLH6TkMDgzLi2DyI/DEQPj5eYa1TeXzW0dz3ZjOVNastxZenbqWEx/7gSkrd0blVxARkbqjgDCyvvQ93mGM8ef2Nsa4gft8iy9FvFUiItL4eRKcjKO3zIFTHoGUQNZRCnPg87vgmREkrvmau0/pxXs3jqJb81T/Lht3F3HJiz/zwCeLKS33RuEXEBGRuqCAMLJeB97EGRa61hjziTHmHWAlcDXwMHB/FNsnIiKNnScBhl/rBIZjfwtxyYFtu1bB/86H/57PwOQcPrnlKG4+pituV6DI/Ys/rmH8c9PYsKswCo0XEZHapiyjEebrDfwr8Jsqm2YDv7fWTjyIc9WURrRL7969E5RlVERE9itvK3z3N5j9H7BBPX/ueDjyVzDm1yzcUc6dE+axbFuef3N6oodHxg/ghD4to9BoEZGGSVlGGzhjzPvGmKUH+W940PFNgW+AXwG3Am2BZsBZQDbwmTHmgij8aiIiEqvSWsLpj8P1P0KHowLrK0rhx0fhyaH03fUVH9w4kguHtfNvzi0u55evzeK+jzWEVESkIVMP4UEwxswFBhzkYUdba7/zHf8qcAVwu7X2sSrnHgpMBzYDnay1ZYfRTtUhFBGRg2ctLHoPvvw/yN0Uuq3DUXDy33l/SxP+8P5CCksr/JuGdGjK0xcPpmVGYoQbLCLSsKiHsIGz1g601pqD/Pcd+IeKXuQ71Tthzj0TWAO0ATpH6ncSERHxMwb6ngu/mgGj73SGjVZa9yM8N5qztzzOJ9f0pUeLNP+mWet2c9qTk5m6SllIRUQaGgWEkdMcqPxm3VvDPpXrm9Z9c0RERGoQnwLH/glu/Am6nxRYb70w/Xk6vzmGT0atZPzgVv5NO/NLufTFn/nXd6vQ6CMRkYZDAWHk7AJKfc+HVt1ojEkHevgW10WqUSIiIjXK7AIXvwUXvw3NggavFOYQ99ntPLT7dl44uoJ4t3M54bXw9y+W8svXZpFbfMgzH0REJIIUEEaItbYE+MK3+Kgxxn9b1RiTCDwDJANTrLVbotBEERGR8Lqf4PQWHnsPxKX4V5stczl+2mX83GsC/TOK/Ou/WryNM578kSVbcqPRWhEROQgKCCPrDmAbMBBYZoyZaIz5AFgFXILTi3h91FonIiJSE08CjL4Dbp4J/caHbGq68j0+9N7Cgy0nEUc5AGtzCjn7mSl8sVD3OEVE6jNlGY0wY0wL4LfAyUBHwAAbgInAg9bajbXwM5RlVERE6ta6afD5XbB1Qcjq3UntuW3vRXzvDSTlvuvEHtw4rgvGmKpnEWmwvF7L0q15LN6Sy/pdhezIKyYxzk16YhztmyUzoF0TOmel4HLpfS8B9THLqALCRkgBoYiIRIS3Ama9At8+AEW7QzZ9z2D+UnIBy61Tu/DiI9pz/5l9ceviWBqwsgovP67YycfzNzN/2SqSCzeRaXLJNLk0I5cmpoAcm8Zq25rVthW5Ca0Y2a0lY3tkc3yvFjRNid//D5FGTQGhRIQCQhERiajCXU5QOOsVJxOpTwUu3ikfwz/Lz2UrmZw5sDWPjB9AnFszVqRh2bY7jy++/pJtiyfTo3wZg8wK2rt27Pe4UutmtW3N+xVH8SYnMLJXBy4a3p7R3bLUYx6jFBBKRCggFBGRqNgyHz7/LayfGrK62MbxSsVJPFt+OiP6dOHpiwfjUVAo9Z23gl2Lv2XFt/+he84kmpq8wzrdLpvKi+Wn8J+KE+jXuR1/OLUXfdtk1FJjD155hReXMRrSGmEKCCUiFBCKiEjUWAtLP4Gv/ww5K0I25dkkXq04kR19rubeC0brQlTqp5xVlE39F6Xz3yOlLGefu9qEdExqc0jOgpQsSGwCeVsgZwV2zwYM1a+z99pkXig/lee9p3P90T257bjuEftbmL5mFxNnLqFs3QyW7faymWxIbUmH7HRO7teS0/q3JiMpLiJtiVUKCCUiFBCKiEjUVZTDnP/Adw9C/raQTfk2kXmtzmfUZfdCSmZ02idS1e51VHz3EGb+G7hsRbXNpXjYldGHpt1HktDxCGg7DNLbQE1DP8uKYdcqmD8BO/0FTFlByObvKgZwXdntjOvTjkfPH0hKgqcufiuwlpWLpjP7q7fouPtHhpjluE3g+r/MutlsM1ljWzGB42g66CxuO74H2WkJddOeGKeAUCJCAaGIiNQbpQXw0zPYqU9hiveEbCp3J+IZcjmMuCG08L1IJOVuxv7wD+ysf+OyZSGbyqybn80ATL9zGHbipcSnNj20n1GQ4/wd/PwcpjQw9HRaRW+uKbuTLm1b8sa1I2o9KLTrf2bPm7+kaeHaAz5makVvHjFXcPwxx/PLMZ2VCKqWKSCUiFBAKCIi9U5xLsVT/0Xp5CdJt1XnYhnoeSqMvBnaHVFzj4tIbaooh5+exjvpr7jKi0M2rfW24BVOJ2v4+Vx53GDSEmtpGGXRbvjoFljykX/VbG9Xriz9DUN7dub5y4bU2vzaohmv4/n0NuIoq7atOL4pHpcLT3H4IbFea3izYhw/tLmO+y49muZpibXSJlFAKBGigFBEROqr7Tt28M6//sSFFR/RzORX36HVABh6FfQ9DxJSI99AiQ1b5uP96GZcW+aGrN5os3iq4hzcgy7mluN70SK9DgKhinL46Fcw7w3/qkXeDlxWejenjOjL/Wf2PbwMpN4KSibeQ8LPT4asXmy6ktjnZDodeQ6m1UBwuZwe/L0bYfdavD89h2v1NyHH5NpkbvP8gVuvvJQB7ZocepvETwGhRIQCQhERqc/mrN/NFc9/z+n2e65yf04X15bqO8WnQf/zYegvoGW/yDdSGqeyYvj+79gpj2OC5gnusOk8UX4Ou3pcwO0n9adr8zq+GeH1wmd3wsyX/auWettxVul9PHjBEZw1qM2hnbckj7IJVxG36kv/qhybxrMt/syNV1y2/zqIK76i4vO7ce8KJITKsWlcZv7GP647k16t0g+tXeKngFAiQgGhiIjUd+/O2sidb8/D4OVo11z+lPktHfNmh9+57TAY8gvoczbEJ0e2odJ4rPwGPrvLSfQS5O3yMbyYfDW/P28UY7tnR6491sKXf4RpT/lXPVR2AROSxvPNHePISD7IYaplRXhfOgHX1vn+VUu87fii7z+59bzjDjyTaUUZdsaLeCf+H27fnMqV3tZcE/c3Xr7+ODpn133P/d6iMrbt2EFS2R4SmrYhs0l6o5nLqIBQIkIBoYiINAS/e3c+b87Y4F/+35npjNz9Ecx7C0r2Vj8gMQMGXOQEh817RrCldW/O+t28OX0Do7tncVr/1tFuTuOydyN8cXfIvD2ADd5s7i6/hsQex/LPCwbW3jzBg2EtfPZrmPEi4AzRHFPyT045og9/Pfsge8anPAFf/Z9/8auKIcwc8iC/O3PYoQ1BnT8B3rvWvzi5oi9/Sr2XT247utaT35Tu3szGif+kcP1cEou20ty7k3RTCECJjWMhXdiSMRBXhyMZNOokWrVsWas/P5IUEEpEKCAUEZGGoKi0grOensKybU6SmabJcUy8fQzNEypg4Xsw6xXYNCv8we2PdOYa9joD4hp2wovXf1rHvR8totxrcRn47NbR9GypoXmHrbwUfnoGvn8Igko+VFjDqxUn8Uj5eMb06cCTFw0m3lM7iVwOSeEueGIgFDs3QZ4vP5W/ll/CuzeMZEiHA8xqWpxL2aP9iCvdAzi9njP638ffxw88vPmI3/4FfnjIv/jf8mNZOuTP3H+wwWoNyksKmfv2X+m18gVSKN7/AUC5dTE59US84/7A2MF9ai0JT6TUx4CwYb2CIiIi0mgkxbt5+pLBJMW5AdhdWMbd7y7AxiXD4Mvg2m/huh+cHsH4KsPU1k9zei8e7QUT/wC710b+F6gFD3yymD9+sJByr3OD3mvhzekb9nOU7JPXC/PfhqeHwdf3hASDs7zdOL30L9xffhmjenXgqYujHAwCJDeDo273L17h/pLW7OThiUsP+BTFPzzhDwbzbSL/y7iGe8/qd3jBIMDRv4e+5/oXL/F8Q/zMfzF11c7DO6+1rPr+DXY8OJChK5884GAQwGO8HF3wOcM/OZ7n/3YLL09aQmFp+eG1J8aph7ARUg+hiIg0JK/9tI7/+2Chf/mh8/pz/tB2oTuV5MGCd5wkHEFzpPxcHmc46ZhfQ9OOddvgWrJo815OfeLHauvTEz1M/8NxJPoCZTlA1sLKr+HrP8O2BSGbShOa8seC83m7fDQWF31ap/P29UeSHF9HxeAPVmkhPDkY8pwES2+Xj+Gu8uv58vYxdG+Rtu9jC3Io/Udf4r3OEMsnK87h6Osfo2+bjNppW1kx3ldPw7VpBuD00F2S+BQv33HBoQ0dLStm5ysXkrV5Usjq3TaNKdnjadd9EG06dCGzVWdIakLRhnlsX/gdZWunkr1rDk3IDTluvTebJzxX0H3cxVwxsiMJnvr9d6MeQhEREZEqLj2iPaO7ZfmX7/t4MRt3F4bulJDmZBy97ge45lsYdCl4kgLbveUw5zV4cgh8+KsG0WM4a91u//P2zZJJjnfRnN1QvJeJi7Ye0DnenrmB4X/5mt+8M4/yCm9dNbX+270O/n06/Pe80GDQ5WF7z8sYXfgwE8rHYnHRIj2Bl64YVn+CQXCSJY272794jnsy3cxGXv9p3X4PLfvhH/5gcLdNxTXy5toLBgHiEnFd9D9KU5y5rR7j5ZiCT/nPtP23LZz17/4+JBgsty5+zByP9+ZZnParfzLghMvJ6jEKk94KE5dEcucRdDzjd3S75SMy/riKtcPuocAVCJLbu3bwiPcRir+8jzOf/JElW3LD/VjZBwWEIiIiElXGGB46rz9pic4Fen5JOXe9PR+vN8woJmOg7RA482m4cymc/DA06RDYHhwYfvsAVFQvyl1fLNqUSzoFnOiazj+SXuHHhNuZnngT8xOvZfjHx8N7v4Sfn3PmUYYZ0VVcVsE9Hy1ie14JE2Zu5OlJq8L8lFpWUQ65W2DbIqeEQ32w9DN4bjSsnRy6vs85rLlwEsctPZ1tZU522tQEDy9dMYyWGfVw3unASyCrOwBuY/mN5y3em72J/JJ9DIfcuwnXjBf8iy/ZM7lkbN/ab1tqc+LH3OZfHO/+ntcmL6W4rKLmY8JYPf1z2i4JlNr4yfZl4ZlfcNTNL5KZ1WK/xxtPPB1PvYOUX8+nYNC1VJhAb+Atng8Yn/MvznhqMq9NW3tQ7Yp1CghFREQk6lplJHHfmYERVNNW5/Dv/V3UJTWBI34JN8+CM56qHhj+8DC8eBzsWF4nbT5k1sKqb7lw2c3MTriO5+IfY1jOhzQrC9RjbFWxCea/BZ//Bl44Bt64yEmSEmTaqhwKSwMX5E98u4K5G/bUXjvLimDJx/DuNfCvUfBwV7g/Cx7tCf8aCY90cwJW78EFBbWmosyZP/rmRf6ELAB0OQZ++R17Tn2OKz7IIbfYCajiPS5euHxo7fae1Sa3B44JZAk93j2LnqULeX/OphoP8f7wMG6v877YZptQMvhqmiTvp9bgoRpwIdbjBNbNTD5HFE3mrRkHPt/VW7ibtC9uxmWcmxurbBviL3+bgYOPOPi2JDcj5cxHcN/4E+VZvfyrr/Z8zn3mBe75cAGvHUDvqjgUEIqIiEi9cNbANpzcN5BO/sHPl7JqR/7+D3THOUloKgPDjPaBbVvmOr1HPz8ftpctorxeJ8B64Wh47WwGlc3FYw5wmOfyz+GT20N+h6+WbAvZpcJruf2tuYeXYKO0EJZ+6gSBD3eFty6FBW/DtoVQsAMIeg1Lcn0B69E1Z4OtK3s3wiunhNTwIy4FznkRLnufipYDufXNuazf5QylNAaeuHAQR3bJjGw7D1av06HNUP/irZ73eG3aWsLm/MhZBbNf8y8+VX42l4/pVX2/2pKYgekXmmDmue9XUVp+YO/hjf+9iWzvDgDKrJu9Jz/N4C6HWWIluzueqz6DVgP9qy7yTOKfcc/w5w/m8s6sjYd3/hihgFBERETqBWMMD5zVl6xUp4ejpNzLH99fGP5iOJzKwPDGaTD4isD68mL4/C5nfllJXh20vLqd+SVc9eoMxj87lWcmrWDHT2/AMyOcAGvznJB9d9k0KvqcC2c+w5MDP2J0yT+5ufRXTM4cDy2Chv/NfR2mPgmA12v5pkpACLBmZ8HBXQSXl8DaKTDpb/DyyfBge3jzYicILN1HMG6CLiG3zIMXjoVP7gjtqasreVudnt+N0wPrmveGX34H/ccD8NjXy/l++Q7/5l+f0IOT+jaA2nXGwNGBuYSDXStYvi3PH9iG+OFhXNYJ/td7s9nb+yLaNUuu2/YNu9r/dKhrOem5y/lgHz2YlYpmv0X7TZ/6lz/L+gWDRxxdO21KbgZXfOSUovE50z2VJ+Oe5LfvzGHWul2183MaMQWEIiIiUm9kpibwl6AaZ9NW5/Dx/C37OCKMhFQ44wm48A1IDiSrYeXX8OYlThBUx/704UK+XbqdwnWzGTrpErK/uB52LgvZ57uKAVxY+kcuafI67vEvw6BLSMpsxwbbgo+9I/l3+nVO6Y2gC12++hMs/YwFm/ayLdf5PVwGTugdmH+1dOsBBL35253hlg91hldPge8fhPVTwVtlzmV8KvQb7/S8XTURbpkLv98Cdy53srr6WZj5Erx+bt0OIa0ogwlX+LNxAjDwUrjmG8h25t9t2lPEv74LzKc8uW9LbhzXpe7aVNta9vc/TTElpFLE9DXVgxq74iv/88fLz+X84RH4HVsPgtaD/YuXur/ijRnr933M3o2YT+/0L860PRhy8b21267EDLj0XegcCDJPds/gItc33PX2/IOe6xhrFBCKiIhIvXJin5Yc27O5f/kvny7ed2KNmvQ8BW78CbqfHFi35nt4/7o6DVqWbMll+oJl/M3zAh/H/5HhruBA0EDvM3mh16tcWfZbfvL2pnebQPHxrNQE//Md+aXgSYALXocmlcNgLbx7DXNnBspVDO3YjKEdA+fYnruPgDd/B3z5R3isvzPcMlwvYEI69D0PLvgv3LUSzn3R6XlrPwKadXIyYqZmw9nPwpWfQlaPwLEbZ8Dc/x7oS3XwvvoTbPgpsHziX+Gsp502+bzww2p/Xcc2TZJ4ePyAw6/HF0nJWU4ZFZ8WZnf1gLCsGFMYqAU4w/Q58CL2h2voVf6nZ7mnsHz9FrbsLapxd/vZXSRWODcp8mwS84b+nbaZ+ymlcSjiU+CiN505pD53e/5Hac5aHv2qns0jrmcUEIqIiEi9c8/pffwFw7fllvDENysO7USp2XDh/5xerkqL3ofPf1tncwp//OAFvk24g4s8k/wJNAC+qhjMnl9MhvP/wzd7A8MXe7dO9z8PDgh35vkCu5QsuHgCxPsuossKOGn+bWThDM88vlcLmqcFsmbuyAuT/dNa+O5BeLy/M+y0POgC3pPk9Kwce49T0uM3a+C8l6DXaRCXVP1cwToeBdf/6Mx9q/TtA1ByAHM/D9bC9+CnZwLLAy6GETeG7JKTX8KbQT1W14/tTOqh1MqLJpcLUgPvj5ZmFzPWVgkI8zb7n3qtoVWbjiTFR6j+Xt9zsYlOYp5UU8xZ7ilMXFhDmZSyYlg+0b/4QMXlnHvsqLprW1winPUsJDnBcYop4e+e53lx8sqQMi8SSgGhiIiI1DvtM5O5YWxgCNyrU9eyLfcQyxy4XHDmM9Dl2MC6GS84WUhrk7Vs++R+rt12H+kmEHCtpg2Xlf6Oa8t+zZLy1lhrWbw5UCutT3BAmBbIELkzvyQwf7J5Lxj/in/uXgu7g9s87wBwXO8WNE8LBJLb88L0EE5/Hr77G5QF5qKVJWZScsx98JtVcPkHMPoOp6SH+yADKE88nPg38PiC0vxtMPWJgzvH/mxf6tSXrNSiH5z6D2fOXZBXp66luMxJcpKVGs/4oe1qtx2Rkt7K/7Qlu1mbU8j24Pd/bmDI7E4yGNK5ORETn4wZcLF/8VL313y+oIZh3dsXYazTG19s49ja8ay6y4JaKa0FnPKIf3GUexEXub7lr58tOfD5yDFGAaGIiIjUSzeM60KLdCfQKS338tz3qw/9ZJ54OP8/0GZIYN2kv8DMVw6zlT7lJfDBDbSYGbgQzTcpeE/6O39o+RyTvc68sBXb89i4u8hfCgFq7iEsKfeGDpXtdjwc80f/4lGuhbRpkkSnrBSapwf3EJaE1nAsyYfvHwosp2TzXtYN9N/zMMf/1J9iUws1+Zq0C+2tm/IE5G6uef+DUZIHEy6DsgJnOSEDLvhPyDBRcOpX/nvqWv/yL0Z1IjEuQr1mtS0tEBC2ME7v4PSgXkLv3kAily22GUd0jnD21KG/8D/t5VpP+fqf2RHmRoTdPM//fIntwAn92kSkefQ9F3qe5l+82/M/tq1fzuQVO/dxUOxSQCgiIiL1UmKcm+vGBHoJ/zd9HTvzDyMhTEIqXPw2ZHYLrPv8N7B3/1kS96lwF7x2Nsx7w79qtbcls094D9eI6+nSsol//fJteSwK6h1s1yyJ9MQ4/3LT5HhcQZ1eO/NDaw/S73z/046ubXRJcXoim6cHAslyr2V3YdBxPz8LlfPN4pLZcem33LFxNEUksn5XITPX1tJQuqNuh5RsXyOKnKGjteGz38DOoDlg5zwPzTpX323+Fn+gnZbg4bIjO1Tbp8FID5RjaGGc/58ZQfMIt29e43++jWaRmz9YKbsHtuNR/sWL3V/z5eLqw0ZzVs30P19kO3BC7whlejUGTn3UP3Q01RTzoOd5HvtqmXoJw1BAKCIiIvXWRcPb+8tQFJd5eXHymv0csR8pmXDZe5DiG2JXURpayw6w1vLurI1c8Nw0znx6Cpe99DOPfrmM3OKy6ucryIGXjod1U/yrfvb25Nq4vzHqCKfgdvcWgQQay7fls2hzoDRDn1ahRdLdLkOzlMCQupyqAXBGWwoTsv2Lg40ztzItwUNiXOCyrjIDKUV7QodvHnEdX64PvSBesKmWSkUkpsPRvw8sz/0fbJl/eOfcPBfm/S+wPPrX0OOksLtOXBQISM4c1Dok0G5w0oLnEDoB4fSgwH1XUEBYltwyKvMkTVBymeNcs5m4oHqPcOmGQImVgmZ9yQ4a2lzn0lrAyYFh4Ue5F9Fj83v8oF7CahQQioiISL2VFO/mmtGB3qDXpq1lT2HpPo44AE3aw+hAGnxmvgIFzkXixt2FXPHKDO58ex4/r9nFvA17mLxiJ098u5LLXvyZvUVVgsJv/gw5K/2L71SM4dLS3zO8T1fcvq6+bs2DA8I8pqwMXJD2bxcaEEKVxDJVA0Jj2JwaKMvRx7vUt9qEJJbZXplYZuqTgdqACekw8ha+Whxav3BhbQWEAIMuh+yevgULX/7h0JP3WOtkFa3UvHdowBmkoKScyUGv60l9WoXdr8FIC+4hdHoGl27N9b//incFerWTsqI0T7Lr8VjfnNYMU8j2NaHlHWx5KZkFgWRQrXoeEfEm0u+8kKGjN3k+5Imvltb5jy0sLWfCzA3c8+FCnvxmBR/O3RQ6B7SeUUAoIiIi9dqlIzrQJNnp7SkoreD9AyiEvV+DLw/UKCwvgp/+xdqdBZzx1BR+CCpoHmzexr1c+uLPgYB081yY/R//9ue9Z/Drsusow8OJfQI9PN1bpPqf7yksY/b6Pf7lMd0CvX2VqpWeqGJNUh//864lS/zPqyWWyd8BP/0rcODIm8l3pzN1ZU7I+cL1EOYWlx3w0LrC0nLW7PTN73N74ISgoaJrfoAVXx7QeapZ9Y1TJqTS8feBK/ycwO+X76C03Ekmk57o4YjOzQ7tZ9YXQUllWrv2AE58PGvdLqy1xBUEkrg0b90p0q1zJKZD88B7cTDLWL2jwL+8fe1CEnAC2DLrZtjwOswuWhNj4JSH/YFrW7MTz8afa/cmSJD8knLu+3gxR/zla37/zmz+PW0d//hqObe+OZfRD03iqW9XUB8HrCogFBERkXotNcHDBcMCvSBvTt9w+POA4pPhyEASFDv9eW5+9Tt2FQQCsHE9snn8woGMH9LWv27Bpr3c9L/ZWK8Xvvgd+C7v8tK78lDpeYAhLcHDyC5Z/mMyUxPITKmeWTErNZ7erdLDrq+0M0yijuXxPf3P2xYtcYq1EzqPcEdeCfz4z0AiluRMGHEDk5fvoLTCG3K+9bsK2VsY6Pn82+dL6H/vl1z+8vT9vs45+SUc88j3HP3Id/yzstZb1+NCCoQHB80HzFsBX90bWO40xjlvDYKHix7bqwVx7gZ+iRvUQ5jFHtw4PW8bdxeRW1ROlg0E9e06do148yqZDkf6nw91LWPF9jz/cu6qGf7na0w7WmU2iWTTAtJbY4JqE57tnswb09fv44CDtGEGTPwDea9dwtqHRvGLGWcwk0tYmXg5sxKu4/P43/Fq3N+5j2cp/OYhduzM2f85I6yB/7WIiIhILLggqHzAsm15zNtYC3f4h13jZKwETEkuo3d/6N/0wFl9eeXKYZw5sA0Pndefq0YFemGmrMxh3sRXYP00/7rXM66nHGce1zG9mvtrKFbqFtRLWGlMt2xcruoF0/c5ZBRYbDtTYp2fFectga0LAEKGjBbvXA8zXgwcdNTtkJBWbbhopcp5jZv3FPmzuU5esZMlW/LC7l/pzRkb2OobCvfq1LVUeK3TKzP8l4GdNs44+GGj8yfAtgWB5ePvq1ZiolJpuZdvl273L5/Yp8XB/az6KGgOoQsvmTiJiIrLKigqKSGbPf7tSZlRLK3RfoT/6TCzjJXbA/UnvVsCGUY3JnYjqvpf6H96ivtnvpi7loLgDL4Hy1pY8TW8ciq8dBxMe4q0VZ/Qt2Ip7Vw7SDDOuTNNHr1c6xnnnscFnu/4TdwEmtvwIxCiSQGhiIiI1Huds1MZ3jEwDPCtGRsO/6SJGTD8Wv/iVZ7PSaSEW4/txqUjOmB8AYgxhv87rRdjuzvDOxMpodX0v/qP83Y/mWc2tPcvn9SneibF4MQylcb2qD5cFCArbd8BYU6JYZHtGFix0emJCe4hHLbhZahwjrWpLWHYNZRXePl2WSBwSggKWiuHjb4za2PIz9q4u5CaVCbfqbS3qCww/LTtsMCO+dtg70bKK7y8NWM9pz05mUte/CmkVzJEWXFohtK+50HrQTW2Y9rqHPJ82UUTPC7GdA//ujYoCanOnE+flr55hMVlXkr3bsVjAr28cU1aVzs8YtoFAsJ2rh3s2BRIdpOcs8j/fG+T3hFtVjU9T8XGOzdl0k0RR5b9zCfzD6EsirWw6AN4bjT891xY92PttjNKFBCKiIhIgxA8bPTjeZspLD2MO/w+E+JOo9A6gVSWyeX+drO49djqvRmVQaHbZbjO/QktrC+BiSuOT1r+KiQgCRfodasSEBoDR3XNqrYfEDK8NCfMHMLconJme4PauGE6EOghzGY3I3M/829+2X0u5a4EZq/fwx5fEBbnNpwf1Ou6YNNevF7LhJmhgfa2fSTCmL1+N6t3FoSs88+/TM2Gph3961fOmcTx//yB3767gIWbcpmyMod/T1sb/sTTn4NcX6DpioNj/6/GNgB8sTAwn25M92yS4yOfcbNOBNUiDASEFVQElUnZa1MwCdV7nyMmow2FyYGANHmbr8yE10t2/jL/+vLm/aoeGVnxyZheZ/gXz3b/yP+mH8JNpW/+DG9f4e+VrzTD252Hy87nj/YGZo5+CW78CW6dB1d/5dQ/PfkhGHUre9qMpaIehl/1r0UiIiIiYZzSrxVpvvT6+SXlfDp/y36O2LcpK3dy9xdbeKMiML/o3OL3cHnD91x1bZ7GTYPiud7zsX/d9r5Xc9/UQNB0xoDWYQOS7s1DL9r7tckgMzV8Cv6qPYTfLdvOOc9M4elJTjbTvUVlNQSEznHnuSfj8c0522yb8eC2I5i2Oicku+nwTs0YEVTMfOGmvUxdlcPG3UUhbdmyt+aAsGpvIsDkFUHD4YJ6CWdM/jKQeCbcvpUKd8HkfwSWh18bElhWVVRawSfzAu+DcL2zDVZ6cHF6p+REcZkXuzfQs7XdRLggfRjlbQLZQzsUzHeS++xaTaJ13ktea0hqPyBazQsYcIH/6VjXPDZuWMfSrbn7OKCKeW8583KDfFMxiPNK/sT40nt5N+VCLr3h9ww99jxo3st537YbDr3PhCOug+Pvo8m1H2Eyozx8NgwFhCIiItIgJMW7OWNgoDfio3mHMOTLZ/WOfG7872wqvJYXyk+hzDf/z5W3Gea/WeNxv6p4nSTj9NrtsBmcMHOYf1hnYpyLO07oHva4qkNGx+5jWGN2amhymN+9u4DZ6/fw8MRlrN1ZQG7VgHDvesjb6hsyarnAPcm/aULFOMrw8MPyHUxbFUhmMbJLFv3aBEperM0p5IXJq6u1paaAsGogVmn2+j3kVdZrDAoIu5dXT/U/Z/2e6vO4Zr4cWiZj9K/D/vxKExdtJc93jpR4Nyf3a0QBYVr14vTF5RWYvMDrvrMeBITJXQMF6oeYZazNKcAGzR9cY1vSpnk9GMbbcbT/NfUYL6e7p/HlovBzaqvZNBs+utm/uNTbjpNKHuTqsruYaXvSvUUq7904kp4tqyeJqircvOFoU0AoIiIiDcZZg9r4n09dlcPugoOvSbi3sIxr/j3TX9NtT1w2uT3GB3ZY9H74A0vyiF8WSDzz9/IL2eNN8i9fc1RnWmUkhTuSpinxtGkS2DauR/Ma2xecVKagtMKftAVg9c588krK2Uomm2xQMLBhOs3TEhnhWkJHl3OR67WGt8vHAvDV4m3M2RAobD6qaxbtmiWRkRQo3v59mHIbW/YWVVsH8OXiQCCWmuDxlwWp8FqmVgaebYf69+9r1hJPGU2T40hPdILvcq9l+ppdgZNaC/PeCCyPuBFS9h3wvD0rMOzv1P6tGs9wUQjpIawsTl9S5sWVH7gRsssd/YDQ03Gk/3kvs441G7dSHFSQfqHtRIfMlGg0LZTLDf0Df+dnu38MSUZUo7xt8OYl/jm5u2wq15TdyVLbHpeBX4zqyHs3jqJ1k/B/+w2BAkIRERFpMIa0b+ofGlnhtXy15ADv8PuUlnu56X+zQ+a+/WP8QDIHnxnYaW/1oZAA5G4G6yTzqDAe3qsY7d+UlRrP9eO67PNn339WH/q2Seemo7swpEPTGvfLTK1eoqJScJ23OSHDRn+maXIcF3kCvYM/ePuzCadnZm1OIWUVTqbPtAQPfVunY4yhb5t992hsraGHcFLQhfQp/VqGzIf0DwVt0Q/czv9Vgimjl1lHtxZpHNklEMQED2Nl0yzIWelfXNf+7H22bcOuQqYE1VQMnhPZKATNIWyBbw5heQVxBYESG7s99aDnLbsnhS5nSLTbWApWT6NsYyAgXOHqTNPkuJqOjqygbKMDXKvJ37SInDCJm/zKS+CtSyHPCcLLrYubym5lo21OvzYZfHjTUdxzeh9SExr2jQgFhCIiItJguFwmpOj75wsOfB6htZbfvTufH4OCkDuO786p/VuFpPknt4ZzBg3Vc6e34pbjeviX7z65134vCo/p2YJPbh7NXSf23Od+cW6Xv8etquC0/rNtUEC4cQameA8nuab7V71ZcTThHNG5GR5fnb6qQ1czkuK4JSipzpa9xWFrEa7cEWjHiM6ZjOkWOM/kFb7X1xMPrQf61w9yraRLdmpI8Bj8fxHcOzjT9GHsC6t48psVYX8HgHdnBwL3zlkp+wyyG6S0cD2EFcQHBYS5cfUgIHS52JoRmCOYtGU6iTsCSVd2pff0Z+yNuha9sS0DCW7Odk3mu2X7KAPx6Z2wMfA3dV/5Zczz9Oee03vzwU2j6Nc2o+ZjG5CGHc6KiIhIzDm5X0te+2kd4AQUe4vKQoY+1uShict4b04gQ+MZA1pz8zG+ot5BF9+U5kFJHiRUKRWRF7gQJ60ltx3XnbHds7HA4Pa1G4xkpSb4M4IGWxUUiC3zBAWWm+fCnP+SgHPMTpvON97BYc99ZJdAQHbVqE40S0mgqKyCQe2a0LNlGgUlFTzhC8RKyr3sLixjT2Epy7flMa5HcxI8rpCeyi7ZqSElL9blFLJ5T5EzhK7tMNjwM+AEhDuyUxgZFBAu3ZrHzvwSshKBhe/6179VOgqAl6as4YZxXfwBbKW9RWUhxcXPG9q2/gQdtSUkqUyg7ERCUXBAWPPQ40gqbnUE7J4CQP89XxPvDdQJLc3uG61mhWUGXOTPEnqmeyoPLt3KuUPaVt9xwwyY85p/8c3ycbztOpnXrx7OkA7Nqu/fgKmHUERERBqU4R2b+UszlFVYvjmAYaPvztrIv75b5V8e0bkZD4/vHwgiUrLBuAMH5IU5Z25QEhtfj+Kg9k1rPRiE0NITwVYFBWKbk7r5h2RSUQLfP+Tf9m7FaMrwhK3lPqprYMimx+3ivCFtuWxEB/q2ycDjdpGe5CE5PvBazFm/m9Oe/JHrX5/Nnz9exNbcYgpLK/zbO2en0CojKaTN/mQ0bYb41w0yK+jSPJXOWSm0TE/0r5+6KgdWfAlFvsQpxPN5xXAA9hSWhc4zxOnp/e0789mW6wz187gM5wwKc0Hf0AUllUk3RSRTTHFpOUnFgeG6BQn1oIcQSO46yv+8jTfQk77Bm01Wdqtwh0RP3/OwxgmB2pqdFKz4gbIKb/X91v7gf7rE25777VU8e9nQRhcMggJCERERaWA8bhcnBA0b/Ww/w0Y37Crkno8CRbJ7tkzjucuGkuAJCgBdbkhtEVjOC5PBNKSHsG6LgQeXngi2KyiJTnJScmjB9pJAr8yEinEAnNK3Fe6grIaZKfF0b16l57MKYwwtMwIB239/Xu8PAD+ZvyVk2GrztATSEp3e2eC5j5XzsopaBHop27t20D2lGGMMo4J6Caeu3AnzApldv6gYSj7J/uWJi4Jed+C1n9bxRdC6247rFtLeRiO1OZjApXoLs5u4sj14vIE5b4WJLcIdGXEtex1Jia0+8HCR7UiHzOQwR0RRWgu8nQLDqY8on82sdbur7WY3zfY//9o7mItHdttnduCGTAGhiIiINDinBJUX+GH5TvYUhs82WuG13DFhLvm+jJjpiR5eunJY+CGmwfMI87ZW3x40hzBk3zqQXUONwmAZSXHQbli19dO9PVhlnWysY3tkM6hdE/+2EV0yDyjtfeugbKk/rgjM88srLg/JzNg5O5A9MjMl0OYcX+C6qqQJ22zg57fMWwiE9lIuXb0Wlk/0Lwcn6wH4cvE2/zzGzXuKeOCTJf5to7tlceO4rvv9fRqkKjcpWppdpJYG5rsV2zi8CU2i0LDqEpNSWOGu/v+w0NuRDs3qWUAIuLsGao8Ocq1g0rLq2UZL18/0P19gO3PZiI6RaFpUKCAUERGRBufIzpn+IYqlFV4+WxAmgANemLyaGWsDd//vP6tvSPmHEOlBvX55YXodQwLCuh0Gl7WPTKOVMpLioO3wauvfKg/0fgzt0DQk++a5g9tU2z+c4B630irD6T4Oqv/YJTvV/zxcD+GqnQUh2VDdm52L7IFBQWr/Pd+A15n7uM02YYo3dM7Zlr3FzN/o9H5OXrHD356s1HgePX9gvazrVmtCMo3uJqMsEBBusc1IDBraG205mdXnrC60HWlf33oIwSkY79PfrGbZptBhyeRtI6Ew8JmS1ml4/fw9aokCQhEREWlwPG4Xpw8IBHAfzN1UbZ8Nuwr551fL/ctnDGjNmQP3ERDtL9NolaQydSnrgHsIQwPC8rhUPvM665qlxNMpK4XxQ9vy+tVH8M71R3JMzwMbYthqH0Mwd+YHemNDAsKgOYT+HsIdBczxBvUcbZwBQMfMFFJ8wcw57sn+zR9UjKICN+2aJTE0KGto5bDR4GQ2Y7s3J7uGobWNRtBNipZmF03KA72122gWOuw5ygaMPKnauuWmU421OaOqZX+8LmeUQJIphW2LQjbvWDbN/3yrbcoZY4bSmCkgFBERkQbp7KAi9dPX7GLj7sKQ7X/5dAkl5ZW9SQncf+Z+sh2GDBmtEhB6vaEBYXodzyEMCgg9NfSApSfFOW1u0t6/zvQbT4eW2RgDN47rgjEGYwxHdctiaMcDT4ZxoBfxIUNGg9qck18ZEOaHBoSbZoO3ApfL0Lt1Ol3MJga6Asl+KoeLnj+kHSf1Dfx/VAaEwUl1gn92oxXcQ2h2k1kRCAi32GYkxtWfgLBJj9ChvtttEzZXNAmZw1pvxCVSlh0oP9GucBGFpeX+5XULf/Q/X+XpxuigOa+NkQJCERERaZD6t82gc1YgKPhwbmAo4+QVO0ISj9x9ck8y9lccO3gYaNU5hEW7/MManX3rtodwYPsmJHicy7RjezUPO+fRv66y2HZcMu4jb+SzW0Yz5/+O55rRnQ/55++rhzBYjUNGC3xDRrfns8B2otz6LjlL82HHUgD6tM7gbHfgwnuxtwPLrBPcjuvRPKTe5KodBWzYVcjqnYGENl1iIiAMvAYtzG4ybY5/eattRmJcPbqUT26GzQ6UQlno7cgxPetHWYxw4joGetcHuVayZmfgZkPKzvn+5+UtG/mwZBQQioiISANljAkZAvr+nE1Ya8kvKQ/JKjqkQ9OQ3sQahQSEVXoIg0tOxCVDQvqhNvuAZKUm8N6NI/nL2X35+7n9w84pTK8MCMfdDVd+BjdOg+zuuFyGJsn7n4O4L62a7D8gTPC4QuZjhiSVyS+lwmtZvbOAIhJZagO9mJXDRoc2LeA8dyC1/7u+3sG0BA+9W6fTrllySIbKORv2sD4n0AvcOSgYbbSqDBltUS0grD89hACm0xj/8wWmO1eO7Bi9xuyHq20gIdNgsyIwHNlaWhcEEhcFl05prFSYXkRERBqsswa15p9fO/MEV27P5+9fLGPBpj3+iztj4M9n9DmwO/xVewitxV/IL2T+YCvCFvirZX1aZ9CndQYA2WkJIcMlwcmYCoDLBR1HVT38sLRK3/+Q0U5ZKSGva3DQujO/lE27iyj1Ddmd4+1KX9daZ+PGGZDUlJMn34zbOMliyq2LjypGAjC8UzP/MMO+rTNY5wsCJy7cSrnXyTZqDLSvh9kra12VIaMpFPuXt9pmdPLUs76d0b+G3Wux7gSuO/UhEtNqv0ZnrQmaf9vRtY0vN2+EAa1hz3oybK5/W1rn6pl8G5t69i4SEREROXAdMlNCShg8+/0qpqwM9KL8cnRn+rbJOLCTBQ8DrSjxF0oHIpphNJzstOo9dmFLZ9SSqsXpAYZ1DL2471Klh65ZUFKZ3YWlrNyR519em9Q7sOP8t2HC5bhLA3UTn6k4gx00AeCIzoG5jn3aBHpiv16yzf+8bdOketc7VieCegibs4dWJriHsGn9ew3SWsAlb2MufL1+B4MAGe3Ijwt8dlT2XOevme5ftc7bnPZt20a6ZRGngFBEREQatH+MH0iPFtWLrZ85sDW/PalnmCNqkNQUPEGBV3AQGMEahOGEGzJalwFh1eL0GUlxnNY/NJFO1Tl8wUllKrzWXyoCYG+zgYEdK4IKq5PE7aU38Gj5+f51IzoHLtL7tg4E85UJggA6Z8XAcFEIufkQZypIN0X+5S02s/4FhA2JMezNHORfzNg1F4D8NTP865a4uoRkz22sFBAeImNMijHmMmPMk8aYn40xJcYYa4y59wCObWuMecUYs9kYU2yMWW6M+bMx5sBmcIuIiIhfy4xEJlx/JCOCepbG9cjmkfEDDi4ZhDE1l56IckAYrrxCXQaEEFqcvn/bDPq3De1prTqHLz3RQ5w78HrPWb/H/zyueTdIbFLlBwzi6e4v8743kJ0yLcFD71aBXsE+rcPP1YyJDKMACWkQV/13LbcudpJRv5LKNERtA+UkOhYuwlqLa/Mc/7otyb0wERgeHm16Fx26bsB/gF8Bw4EDun1gjOkKzAGuBHKADwE38Cfga2NMIy+oIyIiUvsykuL491XD+eOpvbj75J48e+kQ4tyHcJlTU2KZCJacCCc7TF3C9DoOCDtmBeboDWrflF6t0kNKYFQdMmqMCUksM2d9YMht+8wU6HN2YOeRt8BVX9KyU+/gUzC0Y1M8Qf9vmakJtA6T8TQmEsqAc5MivfoQ5e00wYuLBPUQHpaM7oG5t31Zybbd+WTsCSSkys/qH41mRZySyhy6POAlYIbv36nAfQdw3KtAFvCEtfZWAGOMB5gAnA3cDdxb+80VERFp3BI87sMqtQBUqUUYFARGe8hoFHoIrxrViTnr95Ac7+YXIzuSGOfm6qM68dwPqxnVNZO+bar33mWmxrM110l8klscqOvWITMZRv0NOo2BrO7Q0qkJ2afK/M7g4aKV+rTJYPPe4pB1XbJipIcQnJsUOStDVm2zTm94Yj0qTN8QpXYcSjluPFSQYkrYOOc9WlY4yZu81uBuPWg/Z2gcFBAeImvtKuCaymVjzAn7O8YYMxwYBWwHfhN0rnJjzA3AacAtxpgHrLXlNZxGRERE6kpaUO9fcBCYG+WkMlV6CJPj3YfWA3oQOmen8uktocXG7z6lF9eP7UKT5LiwQ+kyw/Rkgi8jaFwS9D0nZH2vlum4XYYKX/bQI8IEhH1bZ/DV4m0h6zrFypBRCNsjvaUyINSQ0cMTn8y6uM50KVsBQJOFr/g3rbStaduq/tZRrE16F0XWqb7Hj621JcEbrLXbgMlAU+CoSDdMREREqNJD6AsCK8qgYEfQPpEPCJtX6SGs697BfWmaEl/jvKqaEnC0q6FERFK8m2uO6gTAcb1aMKBt9YywVXsik+PdtEyPobQLYd5vW/0BoXoID9f29H7+5y12B+YPLrCd6RwjPdEKCCNrgO9xdg3bK9fHxoBlERGR+ibcHML87YAN2ifyQ0abpcSHlD5MT4xeQLgv4QLCjKS4fQawd5/SiwX3nsCLVwwNG2hWLRvSKSslJhJ9+IUNCJ2SDgoID19Ri/CF5+d5O9NRAaHUgfa+x401bK9c3yECbREREZGq0qsUpw9+BCdTZtz+i7bXNo/bRbPkQLAVzR7CfQk3ZLRD5v4LyKftI8BtnpZAVtB5YyahTKUwSWW2WmdorYaMHr7ETkeEXb8hsSepCbExuy42fsv6o/ITrLCG7QW+x+rFlMIwxiyqYVOXg2mUiIiI+AT3xuRvg4pyyNscfnuEZaclkFNQCtR9htFDlRmmXmJNw0UPlDGGvm3S+W6ZM2y3av3DRi+t+hxCfw+hksoctvZderPTppNlcv3ryqybsuw+UWxVZMVsQGiMeR/odZCHXW6tnV4X7REREZF6IHg4qPU6cwdDSk5ENyBcujUPqL89hFlhAsL2hxkQAlxzVGemrsohJd7NeUPaHvb5GpQw77ktaA5hbWnbLIUFaX3Jyp/qX7fUtqNd86ZRbFVkxWxACHQCehzkMYf7iZa/n/NU3vLKO5CTWWvD3rrw9Rz2DrdNRERE9iE+BRIyoGSvs5y3pUrJiegFhMHDJtOT6uclXHAdwkq1ERAe1S2LWX88jniPi4RY6xVLbYEXgytoHut2Xw9hgkdDRmtD9yHHwPeBgHCBtzOdYmT+IMTwHEJr7UBrrTnIf98d5o9d73us6dZW5fp1h/lzRERE5FBVzTSaG90ahJWO7hlIgT+2e3bU2rEvzcIklelQCwEhOPMMYy4YBHDHke8O9Fbl2DRKiCfe48LliqHkOnUooWPoPMJ5tgvdWhzQDK5GoX7eXmq85gFnAoNr2F65fn5kmiMiIiLVpLWEncuc5/Woh/D0/q1onpZAnNvF4PZNotaOfamLOYQCefFZpBftAoKL0sdsv07taz0YjBtsBQDpXYYzumtWlBsVOXonRdanvsfTjTEhYyqMMS2A0cBuYEqkGyYiIiI+aVUyjQbPIYxiQGiMYUTnTIZ0aFpvyy4kx3tIjg/04nlchlYZMVQzsI4UxAd6h7eoBmHtS0iFo24H44I+Z/OHX4zH446dMCl2ftN6wJeQZgrQHPh75XpjjAd4BogDnrDWlkWnhSIiIhJaeqL+9BA2FMG9hG2bJsXUhXVdKUhs4X+uovR15Nj/g99vhvGvQj294VJXNGT0MPgylVZ+M1TmBL7GGHOS7/kWa+3ZVQ77BTANuNUYcwywGBgGdAamAn+r21aLiIjIPgUHfbvWQPGeoG3Rm0PYUGSmJLBhVxGg4aK1ZWXLU+mz9QMM8G7FaEA1COtEFGqM1gcKCA/PIKoXkW/j+wdhksNYa1cYYwYB9wEnAWfjJJu5H/irtbak7porIiIi+xUc9G2ZF7TBQGqLartLqMygxDK1kWFUYE/mIEaUPI3BkkMGoB5CqT0KCA+DtbbjIR63AaenUEREROqb4ELgpfmB56nNwa1Lp/3p0zqdb5ZuB2BQ+9ip5VaXEuNc7CK9yjoFhFI79KkmIiIiEqymYaEaLnpArj6qM8XlXlITPJw1sPX+D5D9SggT/CkglNqigFBEREQkWE3DQpVQ5oBkJMfx+1N6RbsZjUq4AvQqOyG1Re8kERERkWCeeEgOU4NMAaFESbjeQPUQSm1RQCgiIiJSVXqY4E8BoURJ+IBQl/FSO/ROEhEREakqXPCnOYQSJeGGh6qHUGqLAkIRERGRqsIFf+ohlCjRkFGpSwoIRURERKpKC5MdM9wwUpEICBsQKqmM1BK9k0RERESqUg+h1CPh5guGK0UhcigUEIqIiIhUVTX4c8VBUrPotEVinoaMSl1SQCgiIiJSVdXhoWmtwKXLJomORI+yjErd0TtJREREpKqqPYTKMCpRlBAm+AsXJIocCgWEIiIiIlUlZ4EJuuBWQChRlOBxYUzoOg0ZldqigFBERESkKpcrNAhUQhmJImMMCVWyimrIqNQWvZNEREREwgkOAlVyQqKsao+gegiltiggFBEREQmnw5GB5+1HRq8dIlSfM6geQqktnmg3QERERKReGvtbSG8DGe2g/RHRbo3EuKoBYIKSykgtUUAoIiIiEk5CGoy4IdqtEAE0ZFTqjvqaRURERETquYRqAaEu46V26J0kIiIiIlLPJVbLMqoeQqkdCghFREREROo5DRmVuqKAUERERESknqs6RLRqj6HIodI7SURERESkngvuEfS4DB63LuOlduidJCIiIiJSzwXXIdRwUalNCghFREREROq54CGjyjAqtUnvJhERERGRei64V1BF6aU2KSAUEREREanngusQqodQapPeTSIiIiIi9VzokFH1EErtUUAoIiIiIlLPKamM1BUFhCIiIiIi9VxGUpz/eXqiJ4otkcZGAaGIiIiISD13fJ8WdMpKISXezaUjOkS7OdKI6PaCiIiIiEg9l54Yxzd3jKWk3EtSvIaMSu1RD6GIiIiISAPgchkFg1LrFBCKiIiIiIjEKAWEIiIiIiIiMUoBoYiIiIiISIxSQCgiIiIiIhKjFBCKiIiIiIjEKAWEIiIiIiIiMUoBoYiIiIiISIxSQCgiIiIiIhKjFBCKiIiIiIjEKAWEIiIiIiIiMUoBoYiIiIiISIxSQCgiIiIiIhKjFBCKiIiIiIjEKAWEIiIiIiIiMUoBoYiIiIiISIwy1tpot0FqmTEmNyEhIa1Lly7RboqIiIiIiPisWrWKkpKSPGtterTbUkkBYSNkjCnD6f1dGu22xJjKCHxVVFsRW/SaR55e88jTax4det0jT6955Ok1j7yegNdaGxfthlTyRLsBUieWA1hr+0S7IbHEGLMI9LpHkl7zyNNrHnl6zaNDr3vk6TWPPL3mkVf5mtcnmkMoIiIiIiISoxQQioiIiIiIxCgFhCIiIiIiIjFKAaGIiIiIiEiMUkAoIiIiIiISo1R2QkREREREJEaph1BERERERCRGKSAUERERERGJUQoIRUREREREYpQCQhERERERkRilgFBERERERCRGKSAUERERERGJUQoIGxFjTJIx5j5jzHJjTLExZrMx5mVjTJtot62xMcYkG2POMsa8ZIxZ5nu9C4wx84wxfzLGpEa7jbHAGJNpjNlujLHGmJXRbk9jZ4zJNsY84nvPFxljdhljZhtjHo522xojY8wwY8wE32d5mTFmjzFmsjHmF8YYE+32NVTGmCHGmN8ZY94zxmz0fX7stwaXMeZKY8x0Y0y+773/mTFmZCTa3NAdzGtujHEZY0YbYx4yxswyxuQZY0qMMauMMc8aYzpFuv0N0aG+z6uc4+vK44wxbeuqrY3FYXy2xBljbvN9vuT6PmOWR/IaXnUIGwljTCIwCRgBbAEmAx2B4cAOYIS1dnXUGtjIGGOuAV7wLS4BFgLpwEggDVgKjLXWbo9OC2ODMeZV4HLAAKustV2j26LGyxgzBJgIZAKLCLznewNtrbWeKDav0THGnAu8BbiB2cBKIBsYDXiA/1lrL4leCxsuY8wHwJlV11trawyyjTGPAbcCRcCXQCJwLM5nz3nW2g/qoKmNxsG85saYrsAK3+JWYDpQgXM90wbIA06x1v5YV+1tDA7lfV7l+CuBVwCL8z5vZ63dWItNbHQO8bOlGc5nyhCc6/effJu6Av2A0ZF4rysgbCSMMQ8AfwCmASdYa/N96+8A/gF8b60dF70WNi7GmCtwgr/HrLVLgta3Aj4FBgFvWGsvjlITGz1jzLHA18DzwC9RQFhnjDHZwGIgGbjIWvtRle3DrbXTo9K4RsgY4wE2Ac2BS6y1/wva1gv4EWgGHGOtnRSdVjZcxpjfAinADN+/tUBCTRdtxpjjgK+AHOBIa+0K3/ojge+AQqCTtXZPXbe9oTqY19wY0wX4F/AgMMn6LlSNMQnAs8CVwHqgq7W2LBLtb4gO9n1e5dhsnBvbM4EeQAcUEO7XIXy2GOAb4Gjgz8AD1tryoO2dgVxr7c46broCwsbAGBMPbAcygMHW2jlVts8D+gNDrbWzotDEmOK7SJgKlADp1trSKDep0THGJAELcF7js4DlKCCsM8aYZ4AbgJustc9Euz2NnTGmL877e5m1tmeY7Y8DtwC/tdY+FOn2NTbGmGL2fdH2GXAycLu19rEq2yr/L35trf1HXbe1sdjfa76P45JwelEygHHW2u/ron2N0cG85saY/wLnAH1xAhYFhIfgAD5bzscZCfK2tfb8iDauCs0hbBxG4Xw4rqoaDPq843s8PXJNimnzfI8JOMPrpPbdA3QGrgd0h7gO+S7ALgUKcIYPSd0rOcD9cuq0FVL5/j/Gt/hOmF30/RpB1toinBuAAK2j2ZbGyhhzEnAx8Bdr7apot6eRu9b3+GRUW4EzD0EavgG+x9k1bK9c3z8CbREnUAEnUNkVzYY0RsaY/sCdwCvW2snGmI5RblJjNxRnXuyP1toiY8zJwPE4c6iWAxOstZuj2cBGaDWwCuhhjLk4zJDRS4HdwPtRal8s6YFzc29HDb0j+n6NIGOMC6e3Cpz5hVKLjDEpOMN1lwIafVCHjDFxwFFAOTDdd20zHmeqwCbgQ2vtvH2colYpIGwc2vsea+rKr1zfoYbtUrtu9T1+Ya090Dv9cgB8FwMvAnuA30S3NTGjt+9xew0T5v9qjLnaWvtGZJvVeFlrK3zzlD8B/muMuRMnyUZznKQyi4ErrbW64VT39vn9aq0tMMbsAZoaY9KstXkRa1lsugjn72AHztQMqV334SQkHKfpLnWuM86N1W3A7cBfCB25ea8x5nFr7e2RaIyGjDYOlSUOCmvYXuB7TItAW2KaMeYU4Gqc3sH/i3JzGqObgWHAXdZaDZeLjKa+xzOAk4CbcC7IOgKPAEnAv40xA6PRuMbKWjsFGIvTWzgYuAAn8YAXJ8GJskZHxv6+X0HfsRFhjGkHPOZb/JNuuNYuY8xgnBva/9bczIio/G7NBP6GkzCpC5CFcx1ZBNxmjLkpEo1RQChSS4wxPYHXcdIz3xXJrv5YYIxpDzyAkzH31Sg3J5ZUfk94cC7CnrHW7rDWrrPW3gW8DcQBd0WthY2QMeYinHT7G4AjcAKT7sCrOEOmv/VlXRRp9HxDGd/DuVj+wFr7bJSb1KgYY9wERt/8OrqtiRnB362fW2tvstauttbmWGtfJvCdenckGyMNW77vMbmG7Sm+Rw1lqSO+wqFf4NzxedRa+3iUm9QYPQ3E4ySSkcjJD3oeLqlM5bqxEWhLTDDGdAP+DewETrPWTrfWFlhrV1hrr8MZSjoYuCqa7YwR+/t+BX3H1infXKu3ceYz/4iT8ERq12045bJ+E4kSBwLs/7v1Vd9jG19tzjqlOYSNw3rfY9satleuXxeBtsScoKKiHXD+qHV3rW6chnP38lmndI9fou+xjTHmO9/zC621SjhQOyo/NwqttTvCbF/re2wemebEhAtxel2/qKwpW8UEnL+HMTgJIKTu7PP71ddz1QTYrfmDtc83b/zfOGU/5gKn+zKNSu06HacA/RXGmMurbGvpe3zbGFMCPGit/SKirWucgq/J11bdaK0tNMZsx/lubQ6srMvGKCBsHCqHJg6uYXvl+vkRaEtMMcakAp/jJN54D7i2soiu1Ikm1NwTlRi0LbGGfeTgVZaySTLGJISZt9PM9xgucJFDUxl87K1he+X6pjVsl9qzDKcMSLYxpo21dlOV7fp+rVtP4iSSWQ6caK3dE93mNGoG5yZTTUb4Hl+t+6Y0ftbavcaYNUAnwnyW+26GNPEt1vn3q4aMNg5TcC4QutSQ2OE83+PHEWtRDPDN3/kQGA5MBC6y1lZEt1WNl7XWhPuH82EKTh3OyvVro9jURsVaux7nppMhfDBeuS5cDVQ5NJW920Nr2D7M97i27psS23y9Ud/6FseH2UXfr3XEGPMAcCNOL+3x1trtUW5So2WtHbeP79jKnqx2vnWvRrGpjc1HvsdxYbaNwJkmU4RzY6pOKSBsBHypgZ/yLT7tG8ICgDHmDpz6SN9ba2dFo32NkW8C9hs4BYsnA+coRbM0YpX1qB4xxrSqXOm7AXWnb1FJHmrPh77HMcaYG4I3GGNG4KQoh/CF0qX2Pep7/KNvficAxpgjgetwhrK/FIV2NVrGmNuBP+DcHDnOd2NKpLF5DCgFfuX7bAfAGJNFIKPuK5HIqGs0uq1xMMYkAt/hZKPbghOkdPAt7wBGWGuVpryWGGNuJfDH+j6QW8Ouv9YE7brlK0y/BqeHsM4nXscqY8yrwBU4F79TccpNjMQp2v2CtfaXUWtcI2SMeZjAfORFOLUHWwNH4tzMfd6XYEYOkjHmVELLAg3H6QH/OWjd/dbaT4OOeQwnJX8hTtmPeOB433HnWWs/qNtWN2wH85r7bjTN9m2fhjNcNJwXrbU/1kFzG4VDeZ/XcJ61ONeT7ay1NdW7Fg75s+UqnAyv5Tjv9704362ZOH8HY2uYS16rNIewkbDWFhtjjsZJT3sxcBawC2es9//pj7jWBY/3Pnsf+92LkylQpKH7Bc7w9OtwhrdYnC+r56y1/45iuxola+1dxpipOFl1hwA9cLJYfo8TgL8RzfY1cNk4N0urOqLKPn7W2tuMMXOBX+EEgqXA1zgXdyqQvn8H85o3wbmIBucGyJE1nPM7nKyjEt5Bv8/lsB3KZ8vLxpjVwO98+yXh1Jl9EnjEWltABKiHUEREREREJEZpDqGIiIiIiEiMUkAoIiIiIiISoxQQioiIiIiIxCgFhCIiIiIiIjFKAaGIiIiIiEiMUkAoIiIiIiISoxQQioiIiIiIxCgFhCIiIiIiIjFKAaGIiIiIiEiMUkAoIiIiIiISoxQQioiIiIiIxCgFhCIiIiIiIjFKAaGIiIiIiEiMUkAoIiIiIiISoxQQioiIiIiIxCgFhCIiIiIiIjFKAaGIiEg9ZIzpaIyxxpjvjTHNjDGPGmPWGmPKjTGPR7t9IiLSOHii3QAREREJa6DvsQCYj/Od/QOwBJgRpTaJiEgjo4BQRESkfhrgezwZeAG4xVpbHMX2iIhII6SAUEREpH4a6HucClxvrfVGsS0iItJIaQ6hiIhI/TTQ93iXgkEREakrxlob7TaIiIhIEGNMBrAH2Aq0tvqyFhGROqIeQhERkfqncv7gDwoGRUSkLikgFBERqX8G+h7nRLMRIiLS+CkgFBERqX8qewgVEIqISJ1SQCgiIlL/DPQ9KiAUEZE6paQyIiIi9YgxxgPkAznW2jbRbo+IiDRu6iEUERGpX3oCCah3UEREIkA9hCIiIiIiIjFKPYQiIiIiIiIxSgGhiIiIiIhIjFJAKCIiIiIiEqMUEIqIiIiIiMQoBYQiIiIiIiIxSgGhiIiIiIhIjFJAKCIiIiIiEqMUEIqIiIiIiMQoBYQiIiIiIiIxSgGhiIiIiIhIjFJAKCIiIiIiEqMUEIqIiIiIiMQoBYQiIiIiIiIxSgGhiIiIiIhIjFJAKCIiIiIiEqMUEIqIiIiIiMQoBYQiIiIiIiIxSgGhiIiIiIhIjFJAKCIiIiIiEqMUEIqIiIiIiMQoBYQiIiIiIiIxSgGhiIiIiIhIjFJAKCIiIiIiEqMUEIqIiIiIiMQoBYQiIiIiIiIxSgGhiIiIiIhIjFJAKCIiIiIiEqMUEIqIiIiIiMQoBYQiIiIiIiIx6v8B4L9VIA0L9i0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "\n",
    "os.chdir(outdir)\n",
    "\n",
    "cmd.delete_existing_files(\"out320*.txt\")\n",
    "cmd.Execute(\"ScalarFieldCollapse_Playground\", \"320 2 2 \"+str(CFL_FACTOR),\"out320.txt\")\n",
    "\n",
    "os.chdir(os.path.join(\"..\",\"..\"))\n",
    "\n",
    "outfig = os.path.join(outdir,\"ScalarFieldCollapse_H_convergence.png\")\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "r_640,H_640 = np.loadtxt(os.path.join(outdir,\"out640.txt\")).T\n",
    "r_320,H_320 = np.loadtxt(os.path.join(outdir,\"out320.txt\")).T\n",
    "\n",
    "plt.title(\"Plot demonstrating 4th order\\nconvergence of constraint violations\")\n",
    "plt.xlabel(r\"$r$\")\n",
    "plt.ylabel(r\"$\\log_{10}|\\mathcal{H}|$\")\n",
    "plt.xlim(0,16)\n",
    "plt.plot(r_640,H_640,label=r\"$N_{r} = 640$\")\n",
    "plt.plot(r_320,H_320+4*np.log10(320.0/640.0),label=r\"$N_{r} = 320$, mult by $(320/640)^{4}$\")\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(outfig,dpi=150)\n",
    "plt.close(fig)\n",
    "Image(outfig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='output_to_pdf'></a>\n",
    "\n",
    "# Step 8: Output this module as $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{output_to_pdf}$$\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-BSSNCurvilinear-ScalarField_Collapse.pdf](Tutorial-Start_to_Finish-BSSNCurvilinear-ScalarField_Collapse.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": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:13:37.841111Z",
     "iopub.status.busy": "2021-06-15T10:13:37.840831Z",
     "iopub.status.idle": "2021-06-15T10:13:40.961139Z",
     "shell.execute_reply": "2021-06-15T10:13:40.960676Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish-BSSNCurvilinear-ScalarField_Collapse.tex,\n",
      "    and compiled LaTeX file to PDF file Tutorial-Start_to_Finish-\n",
      "    BSSNCurvilinear-ScalarField_Collapse.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-BSSNCurvilinear-ScalarField_Collapse\")"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}