{
 "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: Gaussian pulse initial data for a massless scalar field, in Curvilinear coordinates\n",
    "\n",
    "## Authors: Leonardo Werneck & Zach Etienne\n",
    "\n",
    "# This tutorial notebook sets up initial data for a massless scalar field in *spherical-like coordinates*, using the *Numerical* ADM Spherical to BSSN Curvilinear initial data module (numerical = BSSN $\\lambda^i$'s are computed using finite-difference derivatives instead of exact expressions).\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This module has been validated to exhibit convergence to zero of the Hamiltonian constraint violation at the expected order to the exact solution (see [plots](#convergence) at bottom).</font>\n",
    "\n",
    "## NRPy+ module to obtain the initial data used in this tutorial: **[ScalarFieldCollapse/ScalarField_InitialData.py](../edit/ScalarFieldCollapse/ScalarField_InitialData.py)**\n",
    "\n",
    "## References\n",
    "\n",
    "* [Tutorial-ADM_Initial_Data-ScalarField_GaussianPulse.ipynb](Tutorial-ADM_Initial_Data-ScalarField.ipynb)\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",
    "* [Step 1](#set_nrpy_core): Set core NRPy+ parameters for numerical grids and reference metric\n",
    "* [Step 2](#adm_id_scalar_field): Set up ADM initial data for the scalar field\n",
    "* [Step 3](#adm_id_spacetime): Convert ADM initial data to BSSN-in-curvilinear coordinates\n",
    "* [Step 4](#validate): Validating that the scalar field initial data satisfies the Hamiltonian constraint\n",
    "    * [Step 4.a](#ham_const_output): Output the Hamiltonian constraint \n",
    "    * [Step 4.b](#bc_functs): Set up boundary condition functions for chosen singular, curvilinear coordinate system\n",
    "    * [Step 4.c](#enforce3metric): Enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint\n",
    "    * [Step 4.d](#cparams_rfm_and_domainsize): Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h`\n",
    "* [Step 5](#mainc): `ScalarField_Playground.c`: The Main C Code\n",
    "* [Step 6](#convergence): Validation: Convergence of numerical errors of the Scalar Field Collapse initial data (Hamiltonian constraint violation) to zero\n",
    "* [Step 7](#latex_pdf_output): Output this module as $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='set_nrpy_core'></a>\n",
    "\n",
    "# Step 1: Set core NRPy+ parameters for numerical grids and reference metric \\[Back to [top](#toc)\\]\n",
    "$$\\label{set_nrpy_core}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:06:47.206313Z",
     "iopub.status.busy": "2021-06-15T10:06:47.202597Z",
     "iopub.status.idle": "2021-06-15T10:06:48.422463Z",
     "shell.execute_reply": "2021-06-15T10:06:48.422872Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step P1: Import needed NRPy+ core modules:\n",
    "from outputC import lhrh,outCfunction,outputC # 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(\"ScalarFieldID_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",
    "thismodule = \"ScalarFieldCollapse_ID_setup\"\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",
    "# Then we set the coordinate system for the numerical grid\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"SinhSpherical\")\n",
    "rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.\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     = \"SinhSpherical\"\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 = 15\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 finite difference derivatives;\n",
    "#           and the core data type.\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",
    "\n",
    "# Step 3: 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 4: 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 5: Set the directions=1,2 (theta,phi) axis to be the symmetry axis;\n",
    "#         This sets all spatial derivatives in the thetha and phi directions to zero.\n",
    "par.set_parval_from_str(\"indexedexp::symmetry_axes\",\"12\")\n",
    "\n",
    "# Step 6: The MoLtimestepping interface is only used for memory allocation/deallocation\n",
    "import MoLtimestepping.C_Code_Generation as MoL\n",
    "from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict\n",
    "RK_method = \"Euler\" # DOES NOT MATTER; Again MoL interface is only used for memory alloc/dealloc.\n",
    "RK_order  = Butcher_dict[RK_method][1]\n",
    "cmd.mkdir(os.path.join(Ccodesdir,\"MoLtimestepping/\"))\n",
    "MoL.MoL_C_Code_Generation(RK_method, RHS_string      = \"\", post_RHS_string = \"\",\n",
    "    outdir = os.path.join(Ccodesdir,\"MoLtimestepping/\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='adm_id_scalar_field'></a>\n",
    "\n",
    "# Step 2: Set up ADM initial data for the scalar field \\[Back to [top](#toc)\\]\n",
    "$$\\label{adm_id_scalar_field}$$\n",
    "\n",
    "As documented [in the time-symmetric Gaussian pulse scalar field initial data NRPy+ tutorial notebook](Tutorial-ScalarField-Gaussian_pulse_initial_data-time_symmetric.ipynb), we generate the scalar field initial data and output it to a file.\n",
    "\n",
    "Our initial data solver requires space matrices methods which are implemented in SciPy, so we make sure that it is installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:06:48.428907Z",
     "iopub.status.busy": "2021-06-15T10:06:48.426149Z",
     "iopub.status.idle": "2021-06-15T10:06:49.005468Z",
     "shell.execute_reply": "2021-06-15T10:06:49.004927Z"
    }
   },
   "outputs": [],
   "source": [
    "!pip install scipy > /dev/null"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we call the [`ScalarField.ScalarField_time_symmetric_initial_data.py()` function](../edit/ScalarField/ScalarField_time_symmetric_initial_data.py) ([NRPy+ tutorial notebook](Tutorial-ScalarField-Gaussian_pulse_initial_data-time_symmetric.ipynb)) to set up the initial data, using the default parameters for initial data. This function outputs the solution to a file named \"outputSFID.txt\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:06:49.010891Z",
     "iopub.status.busy": "2021-06-15T10:06:49.010272Z",
     "iopub.status.idle": "2021-06-15T10:06:49.987348Z",
     "shell.execute_reply": "2021-06-15T10:06:49.987090Z"
    }
   },
   "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.1,\n",
      "            center            = 0,\n",
      "            width             = 1,\n",
      "            domain size       = 22.5,\n",
      "            number of points  = 100000,\n",
      "            Initial data file = ScalarFieldID_Ccodes/output/SFID.txt.\n",
      "\n",
      "Output C function ID_scalarfield_ADM_quantities() to file ScalarFieldID_Ccodes/ID_scalarfield_ADM_quantities.h\n",
      "Output C function ID_scalarfield_spherical() to file ScalarFieldID_Ccodes/ID_scalarfield_spherical.h\n",
      "Output C function ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2() to file ScalarFieldID_Ccodes/ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2.h\n",
      "Output C function ID_scalarfield() to file ScalarFieldID_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",
    "initial_data_family = \"Gaussian_pulse\"\n",
    "pulse_amplitude     = 0.1\n",
    "pulse_center        = 0\n",
    "pulse_width         = 1\n",
    "Nr                  = 100000\n",
    "rmax                = domain_size*1.5\n",
    "\n",
    "# Step 2.c: Generate the initial data\n",
    "sfid.ScalarField_InitialData(outputfilename,initial_data_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": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:06:49.989688Z",
     "iopub.status.busy": "2021-06-15T10:06:49.989268Z",
     "iopub.status.idle": "2021-06-15T10:06:50.674812Z",
     "shell.execute_reply": "2021-06-15T10:06:50.675051Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function ID_BSSN_lambdas() to file ScalarFieldID_Ccodes/ID_BSSN_lambdas.h\n",
      "Output C function ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs() to file ScalarFieldID_Ccodes/ID_ADM_xx0xx1xx2_to_BSSN_xx0xx1xx2__ALL_BUT_LAMBDAs.h\n",
      "Output C function ID_BSSN__ALL_BUT_LAMBDAs() to file ScalarFieldID_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='validate'></a>\n",
    "\n",
    "# Step 4: Validating that the scalar field initial data satisfies the Hamiltonian constraint \\[Back to [top](#toc)\\]\n",
    "$$\\label{validate}$$\n",
    "\n",
    "We will validate that the scalar field initial data satisfies the Hamiltonian constraint, modulo numerical finite differencing error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ham_const_output'></a>\n",
    "\n",
    "## Step 4.a: Output the Hamiltonian constraint \\[Back to [top](#toc)\\]\n",
    "$$\\label{ham_const_output}$$\n",
    "\n",
    "First output the Hamiltonian constraint [as documented in the corresponding NRPy+ tutorial notebook](Tutorial-BSSN_constraints.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:06:50.679746Z",
     "iopub.status.busy": "2021-06-15T10:06:50.679310Z",
     "iopub.status.idle": "2021-06-15T10:07:21.788610Z",
     "shell.execute_reply": "2021-06-15T10:07:21.788113Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function Hamiltonian_constraint() to file ScalarFieldID_Ccodes/Hamiltonian_constraint.h\n"
     ]
    }
   ],
   "source": [
    "# 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",
    "import BSSN.Enforce_Detgammahat_Constraint as EGC\n",
    "enforce_detg_constraint_symb_expressions = EGC.Enforce_Detgammahat_Constraint_symb_expressions()\n",
    "\n",
    "# Now register the Hamiltonian as a gridfunction.\n",
    "H = gri.register_gridfunctions(\"AUX\",\"H\")\n",
    "\n",
    "# Declare scalar field gridfunctions\n",
    "import ScalarField.ScalarField_declare_gridfunctions as sfgfs\n",
    "sf,sfM = sfgfs.declare_scalar_field_gridfunctions_if_not_declared_already()\n",
    "\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",
    "# Then define the Hamiltonian constraint and output the optimized C code.\n",
    "import BSSN.BSSN_constraints as bssncon\n",
    "import BSSN.BSSN_stress_energy_source_terms as Bsest\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",
    "# 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",
    "\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\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bc_functs'></a>\n",
    "\n",
    "## Step 4.b: 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": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:07:21.791677Z",
     "iopub.status.busy": "2021-06-15T10:07:21.791148Z",
     "iopub.status.idle": "2021-06-15T10:07:22.509485Z",
     "shell.execute_reply": "2021-06-15T10:07:22.508996Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"ScalarFieldID_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",
      "\n",
      "Wrote to file \"ScalarFieldID_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='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": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:07:22.587824Z",
     "iopub.status.busy": "2021-06-15T10:07:22.552088Z",
     "iopub.status.idle": "2021-06-15T10:07:22.593846Z",
     "shell.execute_reply": "2021-06-15T10:07:22.593380Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function enforce_detgammahat_constraint() to file ScalarFieldID_Ccodes/enforce_detgammahat_constraint.h\n"
     ]
    }
   ],
   "source": [
    "# Set up the C function for the det(gammahat) = det(gammabar)\n",
    "EGC.output_Enforce_Detgammahat_Constraint_Ccode(Ccodesdir,\n",
    "                                                exprs=enforce_detg_constraint_symb_expressions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cparams_rfm_and_domainsize'></a>\n",
    "\n",
    "## Step 4.d: 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": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:07:22.597369Z",
     "iopub.status.busy": "2021-06-15T10:07:22.596729Z",
     "iopub.status.idle": "2021-06-15T10:07:22.943834Z",
     "shell.execute_reply": "2021-06-15T10:07:22.943345Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.d.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 3.d.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 3.d.iii: Generate set_Nxx_dxx_invdx_params__and__xx.h:\n",
    "rfm.set_Nxx_dxx_invdx_params__and__xx_h(Ccodesdir)\n",
    "\n",
    "# Step 3.d.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 3.d.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='mainc'></a>\n",
    "\n",
    "# Step 5: `ScalarField_Playground.c`: The Main C Code \\[Back to [top](#toc)\\]\n",
    "$$\\label{mainc}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:07:22.947109Z",
     "iopub.status.busy": "2021-06-15T10:07:22.946444Z",
     "iopub.status.idle": "2021-06-15T10:07:22.948397Z",
     "shell.execute_reply": "2021-06-15T10:07:22.947771Z"
    }
   },
   "outputs": [],
   "source": [
    "# Part P0: Define REAL, set the number of ghost cells NGHOSTS (from NRPy+'s FD_CENTDERIVS_ORDER)\n",
    "\n",
    "with open(os.path.join(Ccodesdir,\"ScalarField_Playground_REAL__NGHOSTS.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\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:07:22.953476Z",
     "iopub.status.busy": "2021-06-15T10:07:22.952903Z",
     "iopub.status.idle": "2021-06-15T10:07:22.954936Z",
     "shell.execute_reply": "2021-06-15T10:07:22.955306Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing ScalarFieldID_Ccodes//ScalarField_Playground.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/ScalarField_Playground.c\n",
    "\n",
    "\n",
    "// Step P0: Define REAL and NGHOSTS. This header is generated by NRPy+.\n",
    "#include \"ScalarField_Playground_REAL__NGHOSTS.h\"\n",
    "\n",
    "#include \"rfm_files/rfm_struct__declare.h\"\n",
    "\n",
    "#include \"declare_Cparameters_struct.h\"\n",
    "\n",
    "// Step P1: Import needed header files\n",
    "#include \"stdio.h\"\n",
    "#include \"stdlib.h\"\n",
    "#include \"math.h\"\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",
    "\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 gridfunctions macros\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: Include function for enforcing detgammabar constraint.\n",
    "#include \"enforce_detgammahat_constraint.h\"\n",
    "\n",
    "// Step P8: 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 P9: 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 P10: 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",
    "// 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 2D 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) || 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: ./ScalarField_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",
    "    // 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 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\n",
    "    initial_data(&params,&bcstruct, &rfmstruct, xx, auxevol_gfs, y_n_gfs);\n",
    "    // Step 1a: Compute the energy-momentum tensor\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",
    "    // 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",
    "        int idx = IDX3S(i0,i1mid,i2mid);\n",
    "        REAL rr = sqrt( xCart[0]*xCart[0] + xCart[1]*xCart[1] + xCart[2]*xCart[2] );\n",
    "        fprintf(out1D,\"%e %e\\n\",rr,log10(fabs(diagnostic_output_gfs[IDX4ptS(HGF,idx)])));\n",
    "    }\n",
    "    fclose(out1D);\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": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:07:22.962407Z",
     "iopub.status.busy": "2021-06-15T10:07:22.958228Z",
     "iopub.status.idle": "2021-06-15T10:07:25.086736Z",
     "shell.execute_reply": "2021-06-15T10:07:25.087235Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling executable...\n",
      "(EXEC): Executing `gcc -std=gnu99 -Ofast -fopenmp -march=native -funroll-loops ScalarFieldID_Ccodes/ScalarField_Playground.c -o ScalarFieldID_Ccodes/output/ScalarField_Playground -lm`...\n",
      "(BENCH): Finished executing in 2.8468871116638184 seconds.\n",
      "Finished compilation.\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd\n",
    "cmd.C_compile(os.path.join(Ccodesdir,\"ScalarField_Playground.c\"),\n",
    "              os.path.join(outdir,\"ScalarField_Playground\"),compile_mode=\"optimized\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(EXEC): Executing `./ScalarField_Playground 96 2 2`...\n",
      "(BENCH): Finished executing in 0.4174671173095703 seconds.\n",
      "(EXEC): Executing `./ScalarField_Playground 72 2 2`...\n",
      "(BENCH): Finished executing in 0.2159571647644043 seconds.\n",
      "(EXEC): Executing `./ScalarField_Playground 48 2 2`...\n",
      "(BENCH): Finished executing in 0.21025490760803223 seconds.\n"
     ]
    }
   ],
   "source": [
    "os.chdir(outdir)\n",
    "cmd.delete_existing_files(\"out96.txt\")\n",
    "cmd.Execute(\"ScalarField_Playground\", \"96 2 2\", \"out96.txt\")\n",
    "cmd.delete_existing_files(\"out72.txt\")\n",
    "cmd.Execute(\"ScalarField_Playground\", \"72 2 2\", \"out72.txt\")\n",
    "cmd.delete_existing_files(\"out48.txt\")\n",
    "cmd.Execute(\"ScalarField_Playground\", \"48 2 2\", \"out48.txt\")\n",
    "os.chdir(os.path.join(\"..\",\"..\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convergence'></a>\n",
    "\n",
    "# Step 6: Validation: Convergence of numerical errors (Hamiltonian constraint violation) to zero \\[Back to [top](#toc)\\]\n",
    "$$\\label{convergence}$$\n",
    "\n",
    "The equations behind these initial data solve Einstein's equations exactly, at a single instant in time. One reflection of this solution is that the Hamiltonian constraint violation should be exactly zero in the initial data. \n",
    "\n",
    "However, when evaluated on numerical grids, the Hamiltonian constraint violation will *not* generally evaluate to zero due to the associated numerical derivatives not being exact. However, these numerical derivatives (finite difference derivatives in this case) should *converge* to the exact derivatives as the density of numerical sampling points approaches infinity.\n",
    "\n",
    "In this case, all of our finite difference derivatives agree with the exact solution, with an error term that drops with the uniform gridspacing to the fourth power: $\\left(\\Delta x^i\\right)^4$. \n",
    "\n",
    "Here, as in the [Start-to-Finish Scalar Wave (Cartesian grids) NRPy+ tutorial](Tutorial-Start_to_Finish-ScalarWave.ipynb) and the [Start-to-Finish Scalar Wave (curvilinear grids) NRPy+ tutorial](Tutorial-Start_to_Finish-ScalarWaveCurvilinear.ipynb) we confirm this convergence.\n",
    "\n",
    "We demonstrate convergence along the radial direction in our 1D (spherically symmetric) problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:07:25.095069Z",
     "iopub.status.busy": "2021-06-15T10:07:25.094469Z",
     "iopub.status.idle": "2021-06-15T10:07:25.514885Z",
     "shell.execute_reply": "2021-06-15T10:07:25.515197Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAHgCAYAAACMxVqsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAABJ0AAASdAHeZh94AACvxUlEQVR4nOzdd1zU9R/A8ddxbBABZbkHIOJW3Av33iPT3LnNkaU2NbOy8assNS0TLbc5UnPknmluE/eeCCpT9t3n98fJCYIKChzg+/l43APuO9+33/e5z+f90SilFEIIIYQQQoh0MTN1AEIIIYQQQuQmkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkCLV8rVq1fRaDT069fP1KGIPKhEiRKUKFHC1GFkOY1Gg7+/v6nDeCHyHiCEyAySQItcT6PRpLhotVoKFixI48aNWbx4cZacc/78+Wg0GubPn5/hfZ+M18rKChcXF6pWrcqbb77Jxo0b0el0mR90HvMyj8GL8vf3R6PRZNv5soJSimbNmhmff4mJiam2eVW+CJiaXq/njz/+oEuXLhQtWhRra2vs7OwoW7YsgwcPZt++faYOUQjxFOamDkCIzDJp0iQAEhISOHv2LH/++Sc7duzg8OHDfPvttyaOLrWkeHU6HWFhYQQGBvL777/z66+/4ufnx6JFi/D29jZxlCIjtm3bZuoQnmvGjBns2LEDa2trYmNjTR3OKysoKIiuXbuyb98+8uXLR7NmzShdujRKKS5cuMCSJUv45Zdf+PHHHxk5cqSpwxVCPEESaJFnTJ48OcX1bdu20axZM77//ntGjRqV41rUnowX4O7du7z11lusWLGCpk2bcvjwYVxdXbM/OPFCSpcubeoQnuncuXNMmDCBd955h6VLl3Lt2jVTh/RKio6OpmXLlpw4cYIePXowa9YsnJycUmwTERHBN998Q3h4uImiFEI8kxIilwPU057KZcuWVYBavny5UkqpK1euKED17ds31ba3b99Ww4cPV8WLF1cWFhaqYMGCqlOnTurw4cMptmvYsKHxnE9erly58lLxKqWUTqdT/v7+ClCjR49Otf7+/ftq4sSJysfHR1lbWysHBwfVuHFjtXnz5lTbBgQEKEAFBASov//+W9WrV0/Z2dmpggULqn79+qnQ0FCllFJHjx5Vbdq0UY6OjsrOzk61a9fuqbfl/Pnzqnfv3qpQoULKwsJCeXh4qN69e6vz58+n2nbSpEkKUDt27FArVqxQ1atXVzY2NsrJyUm99tpr6ubNm6n2uXTpkho0aJAqXbq0sra2Vk5OTqp8+fJqyJAh6t69e0qp9D0Gyc+9aNEiVaNGDWVnZ6eKFy+e4v7p3LmzKlmypLK2tlb58uVTderUUb///nuKmJKeN2ldGjZsaNyuePHiKY7/5GOwfft21bBhQ2Vvb6/y5cunWrdurU6fPp3m/Xzu3DnVuXNn5ejoqGxtbVXt2rXV+vXrUxwvIxISElT16tVVuXLlVGxsrCpevLgCVEJCgnGbHTt2PPV2Jn/NJN3ukJAQNWjQIOXu7q4sLS2Vr6+vmjdvXobiSpLe159SKe/TjRs3qoYNGyoHB4cUr6uIiAg1duxYVbhwYWVlZaXKlCmj/ve//6lLly499T3g4cOH6vPPP1eVKlVStra2ys7OTtWqVUstXrw41bZJ99WkSZPUwYMHVevWrZWTk1O63gemTp2qAFW3bl2l0+meuW1sbGyK62FhYWrixInK29tbWVlZKUdHR9W8eXO1ZcuWZ8Z47Ngx1bp1a5U/f35lY2OjGjRooPbt25di+yFDhihArVmzJs1YDhw4oADVpUuXFMsz+36LjY1VkyZNUiVLllSWlpaqRIkS6oMPPlCxsbGpXnNJEhIS1MyZM1XNmjVVvnz5lI2NjapcubL68ccfU93HyT8Hrly5ol577TVVoEABZWVlpapVq6bWrVuX5u1XSqmlS5eqxo0bKycnJ2VlZaWKFy+uevTooQ4dOpRq28WLFyt/f3+VP39+ZWVlpXx8fNSnn36a6jEVuZMk0CLXe1ZC6uPjowC1YsUKpdTTE+jLly+rQoUKKUA1btxYTZw4UfXq1UtZWloqS0vLFG+oAQEBqkOHDgpQHTp0UJMmTTJekhLSF403ydatWxWgXF1dlV6vNy6/evWqKlGihAJU/fr11ZgxY9SgQYOUh4eH0mg06ueff05xnKREo1OnTsrCwkJ16tRJjRs3TtWuXVsByt/fX/3zzz/K1tZWtWjRQo0bN041b95cAapcuXKpPnj+/fdf5eDgoDQajerQoYN67733VKdOnZRGo1EODg7q33//TbF9UhLbrVs3ZWVlpbp166beeecdVb9+fQUoHx+fFB8mt2/fVs7Ozsrc3Fy1b99ejR8/Xo0aNUq1a9dO2draqv/++y/dj0HSudu2bausrKxU165d1YQJE9TQoUON57O2tlbVqlVTffv2VRMnTlSDBg1ShQsXVoD68MMPjduFhoaqSZMmGZPO5OdLnsg+K4Hu0qWLMjc3V+3atVPvvPOOat26tQKUi4uLCgkJSbHPmTNnjElFmzZt1Hvvvadee+01ZWFhYbzdGU2gJ02apMzNzY0f9Gkl0FeuXFGTJk1S+fPnV/nz509xO1evXm3cDlCVKlVS3t7eqnz58mrkyJFq0KBBytHRUQFq/vz5GYotI6+/5PdpmzZtlFarVW3btlXjx49Xr732mlLKkIBVr17dGOf48ePV4MGDlaOjo2rfvn2a7wGhoaGqSpUqClBVq1ZVI0eOVMOHD1elS5dWgPrggw9SbJ+UCDZv3lxZWlqqxo0bq3Hjxqm+ffuqW7duPfP2FitWTAFq06ZNGbqfQkNDla+vrwJU9erV1YQJE9TAgQNVvnz5lEajUbNnz04zxjZt2igbGxtjjN26dVNmZmbK2tpanT171rj9/v37FaA6d+6c5vmHDRumgBSPR2bfb3q9XrVp00YBysvLS7399tvqrbfeUoULF1YdO3ZMM4GOj49XLVq0UIAqU6aMGjJkiBo9erSqWLGiAtQbb7yRYvukzwF/f3/l4uKiatasqcaMGaP69OmjrKyslJmZmdq+fXuKffR6verbt68CVMGCBdXAgQPVxIkT1RtvvKEKFy6sJk2alGL7/v37K0AVKVJEDRgwQL399tuqTp06xvMmf92J3EkSaJHrPS0h3bJli9JoNEqj0airV68qpZ6eQCcljVOnTk2xfN++fUqr1SpnZ2cVGRlpXP6irYDPije52NhYZW5urgB1+fJl4/KGDRsqjUajlixZkmL70NBQValSJWVtba2CgoJSxanVatXOnTuNy3U6nWratKkClJOTk1q4cGGK4w0YMCBVS5Rerzd+IXly+6VLlxo/vJIn3UlJbL58+dTJkydT7PP6668rQC1btsy47IcfflCA+v7771PdJ1FRUSo6OjrVbXvaY5B0bltbW3X06NE0t7l48WKqZXFxcapx48bK3Nw8VQt5Usv30zwrgdZqtWrr1q0p1k2cOFEB6ssvv0yxvHHjxgpQs2bNSrF8w4YNxudPRp57//77rzI3N0/xpSCtBPpZtyO5pBgGDhyoEhMTjcsDAwOVVqtVZcuWTXdsSr3460+j0aiNGzemOt5nn31mTASTPx8vX75s/GLy5HtAUnL05GMRExOjWrRooTQajTp27JhxefLW+icT12e5fv26ApS5ubmKiYlJ935KKTV48GAFqMGDB6f4Yn3+/Hnl4OCgLC0tU7R+J4/xyefL7NmzFaCGDRuWYrm3t7eytLRU9+/fT7E8NjZWOTk5KVdX1xTPmcy+33777TdjA0FcXJxxeWhoqCpTpkyaCXTSa33kyJEpno+JiYlpvpcl/0Vp8uTJKY61adMmBahWrVqlWD5nzhzjF5ewsLAU6xITE9Xt27eN15M3XCR/z0oea1rvcSJ3kQRa5HpJb4RJLWXvv/++6tKli9JqtQpQY8eONW6bVgJ948YNBahixYqp+Pj4VMd/4403FKAWLFhgXJbVCbRSSrm5uSlAHTx4UCml1PHjxxWgunbtmub2a9asUYCaOXNmqjifbIFRSqkFCxYYP6ietHPnzlQfLnv37lWAql27dprnr1evngLUrl27jMuSPiyebIVSSqnt27crQI0bN864LCmBnjNnTprnSC69CfSYMWOee6wnrVy5MtVjrtTLJdC9evVKtf3ly5eNrdNJkhIsT0/PNH/eT/rik97nXnR0tCpTpoyqVKlSiuf3yybQtra2Kjw8PNW6Bg0aKCBFwvssL/P669ixY5rH9PT0VGZmZml+QUp6XiR/D7h3757SarXKz88vzeMlvfbeffdd47KkRLBy5crpup1JDh48qADl5uaWof3i4uKUra2tsre3T5XcKqXUhx9+qAD1ySefpIqxbt26qbaPj49X5ubmqlq1aimWJ335mDFjRorlK1asSPV+mhX3W5MmTVK9jyRZuHBhqgRap9MpZ2dn5e7unuZzOTQ0VGk0GtWtWzfjsqTPgeLFi6dIuJMUK1ZMFShQIMWy8uXLK+CpX8aTq1y5sjI3N0/zF8nExERVoEABVb169eceR+RsMohQ5BmffPIJYCgT5+joSP369Rk4cCBvvPHGM/c7duwYAPXr18fCwiLV+saNG7Nw4UKOHTtGnz59Mj/wp1BKARjLpv3zzz8AhIeHpzkAMSQkBIAzZ86kWufn55dqWaFChQCoVq1aqnWFCxcG4ObNm8ZlR48eBQz3R1oaN27M3r17OXbsGA0aNHju+YsWLQpAaGiocVn79u15//33GTFiBJs3b6ZFixbUrVsXX1/fFy4fV6NGjaeuu379Ol9++SXbtm3j+vXrxMTEpFh/69atFzpnWtJ7Hxw/fhyA2rVrY2aWutJovXr12Lp1a7rPO378eC5fvsyhQ4fSfH6/KC8vLxwcHFItT36b7O3tCQsL4/vvv0+13ZgxY3B0dHyp119aj21kZCQXL16kaNGiaQ7q9Pf3N75XJDl06BA6nQ6NRpPmayshIQFI+7X1rOdXZjp37hzR0dHUrVsXZ2fnVOsbN27M1KlTjfdncmk99ywsLHBzc0vx3APo06cPH330EQsWLGDEiBHG5QsWLABIUT87K+63Y8eOYWZmRp06dVKtq1evXqpl58+f58GDB3h5eTF16tQ0j2ljY5NmDJUrV0ar1aZaXrRoUeP7LcDDhw85deoUbm5uVKlSJc1zJImOjubEiRMULFgwzec9gJWVVZrxiNxFEmiRZyQlnBmVNMrdw8MjzfVJy8PCwl7o+C8iNjaWBw8eAODi4gLA/fv3AdiyZQtbtmx56r5RUVGpluXPnz/VMnNz8+euS/oAhJe7nxwdHZ96juQ1r4sXL86///7L5MmT2bRpE6tWrQIMH2jvvPMOo0aNSvPcz+Lu7p7m8suXL1OjRg1CQ0OpX78+zZs3J3/+/Gi1Wq5evcqCBQuIi4vL8PmeJr33QdL97ObmluZxnrY8Lbt27WLmzJlMnjyZSpUqZSDa50vr9kDq2xQWFpYqYQVDIubo6PhSz6u0Htvn3X9p7ZP02jp06BCHDh1Kcz9I+7X1tOfX0yTdnvv37xMbG4u1tXW69svs1x8YHqsna84XKVKEJk2asGXLFs6cOUPZsmUJDg5m06ZNVK5cmYoVKxq3zYr7LTw8HGdnZ+PzKLm0HtOkGC5cuJDm8+xZMTzrftHr9cbrSfdpUsPCs4SGhqKUIiQk5JnxiNxPJlIRr7ykBDIoKCjN9Xfu3EmxXXbYu3cviYmJuLm5GcvvJZ1/+vTpKEP3qzQvAQEBWRJTdt1PZcuWZdmyZdy/f5/Dhw8zbdo09Ho9o0eP5tdff83w8Z7Wcv3tt99y//59fv31V3bu3MkPP/zAp59+yuTJk2nRosVL3YaXkdSqe/fu3TTXP215Wo4dO4ZSikmTJqWawCephJ2FhQUajcbY8p3ZSpQokebz9Mnn9Ys8r9J6bJO2e9r9lNZ5kvYZO3bsM19bO3bsSFcMz1K0aFGKFStGYmIiu3fvTvd+2fk+1bdvX+Bxq/OiRYtITEw0Ln8ypsy83xwcHHjw4EGaE/yk9ZgmxdCpU6dnxnDlypUM3AMpJSXa6flFKimeKlWqPDOeF23wETmHJNDilZf0k1xS0vqkpDf/qlWrGpcl/eyXFTMG6vV6PvvsMwB69uxpXF6rVi0A9uzZk+nnTI+k+2nnzp1prk/rfnoZ5ubmVKtWjQkTJrBkyRIA1qxZY1z/so/BxYsXAejSpUuqdbt27Upzn6x83JNUrlwZMHTZSd4KlmTv3r3pPlb58uUZOHBgmhd7e3sABgwYwMCBAylQoIBxP61Wm22zYb7I6+9Z8uXLh6enJ7du3eLSpUup1qf1/K1RowZmZmbZ9toaPHgwAFOnTk3zMU4u6VeQMmXKYGtry4kTJ9JsZc7M11/nzp1xcHBg4cKF6PV6FixYgLm5eYr3I8ia+61KlSro9Xr279+fal1az30fHx8cHR05cOBAil/MMpOdnR3ly5fn7t27aXaRSc7e3p5y5coRGBho/BVR5E2SQItXXpEiRWjWrBlXr15N1Wft4MGDLF68GCcnJzp16mRcnpRsXL9+PVNjCQ4OpkePHuzcuZNixYrx/vvvG9f5+flRv359Vq1axbx589Lc/7///iM4ODhTY0pSt25dypQpw969e/njjz9SrPvjjz/Ys2cP3t7eafZTTK8jR46kOXFEUsuTra2tcdnLPgZJLaBPJlSbN29m7ty5ae6TVY97csWKFcPf35+LFy8yZ86cFOs2bdqUof7PTZs2Ze7cuWlekm7LnDlzmDt3rrHvMhhuZ0hISKo+4VnhRV5/z9O/f3/0ej0TJkxIkaBeuXKFH374IdX2rq6u9OrVi8OHD/Ppp5+m+eXh0qVLL9WKmdzYsWOpVKkSe/bsoU+fPmkmxFFRUXzyySd88803AFhaWtKrVy8iIyP56KOPUsX2ww8/YGFhQe/evV86PhsbG7p3786tW7f47rvvOHHiBK1bt041qVNW3G9J/dw//PBD4uPjjcvDw8P59NNPU21vbm7OW2+9xZ07dxg1alSaz9k7d+5w+vTpdMeQlqTuY0OGDEn1HqXX642/AAC8/fbbxMfHM2DAgDQf29DQUOOYEpF7SR9oIYDZs2dTt25d3n33Xf7++2/8/Py4ceMGK1aswMzMjICAAPLly2fcvnbt2tja2vL9999z//59Y3++t956K90/oSYNutHr9capvPfu3Ut8fDw1atRg0aJFFCxYMMU+ixcvpnHjxgwcOJAffviBmjVr4ujoyM2bNzl58iSnTp3in3/+yZLZCzUaDQsWLKBZs2a89tprdOjQAR8fH86dO8eaNWvIly8fv/32W5oD39Lr999/Z86cOdSrV4/SpUvj5OTEpUuXWLduHVZWVowZM8a47cs+BsOHDycgIIBu3brRtWtXChUqxKlTp9i0aRPdu3dn2bJlqfZp0qQJK1asoHPnzrRu3RobGxuKFy+eKUlLcjNnzqRu3boMHz6cDRs2ULFiRS5fvszKlSvp0KEDf/7550vdz8/TpEkTDh06RMuWLWnQoAFWVlZUqlSJdu3aZcn5Mvr6e55x48axZs0aVq5cSdWqVWnRogVhYWEsX76cBg0asHbt2lT7zJgxgwsXLvDxxx/z+++/U69ePdzc3Lh9+zZnzpzh0KFDLFmyhJIlS7707bW1tWXTpk107dqVRYsWsW7duhRTeV+8eJFt27YRERHBjBkzjPtNmzaNPXv2MGPGDA4dOkSjRo24d+8ey5cvJzIykhkzZmRKfGDoxjF37lzee+894/W0ZPb91qdPH5YuXcqmTZsoX7487du3JyEhgZUrV1K9enXOnTuX6rn/0UcfceLECWbPns26deto3LgxhQsXJjg4mAsXLrBv3z4+++wzfH19X/j+ePPNN9mzZw+///47Xl5edOjQARcXF27fvs327dsZMGCA8T19wIABHDlyhFmzZlG6dGlatGhBsWLFePDgAVeuXGH37t3079+f2bNnv3A8IgfIitIeQmQn0lkWTqlnz0R48+ZNNXToUFWsWDFlYWGhChQooDp06JBqcpAkGzduVLVq1VJ2dnbGGDIyE2HSxdLSUhUoUEBVrVpVvfnmm2rjxo3PnJ0sIiJCffbZZ6pq1arKzs5OWVtbqxIlSqjWrVurOXPmqKioKOO2zyr1lnxGsCc96346e/aseuONN5S7u7syNzdX7u7uqlevXikmZEiSfDbA9JzjwIEDaujQoapixYrKyclJWVtbq9KlS6t+/foZJ1FJ7lmPwbPOnWTfvn2qUaNGytHRUdnb26u6deuq1atXP/W+SUxMVO+9954qWbKksU53RmYiTMuTx0hy5swZ1alTJ5U/f35la2uratWqpdavX6++/vprBaSY2ORFPKuMXVRUlBo6dKgqXLiwsRxkWjMRpiWpLnB6XgvJZeT1l54ykuHh4Wrs2LGqUKFCxpkIv/nmm2fORBgXF6d+/PFHVbt2bWNd5aJFi6rGjRur7777zjgTplLPfv2kl06nU8uXL1edOnUyzphoY2OjypQpowYOHJhqpkClDGXZxo8frzw9PZWlpaXKnz+/atq0aZozkT4vxueVK/T09FSAcnZ2TlGT+UmZfb/FxMSojz76SJUoUUJZWlqq4sWLq/fff1/dvHlTgWHypCfp9Xr122+/GWcJtLCwUIUKFVJ169ZVn332mbp+/bpx22e9vyn17HKVCxcuVA0aNFAODg7KyspKlShRQvXs2VMdOXIk1bbr1q1Tbdq0US4uLsrCwkK5ubmp6tWrqw8++ECdOXPmqbdf5A4apaQnuxBC5Ba9evVi8eLFnD17ljJlypg6HCGyzZYtW2jevDkTJ07kiy++MHU44hUnfaCFECKH0ev1aVZb2LZtG8uWLcPX11eSZ5Fn3b59O9Wy+/fvM3HiRIAM9YcXIqtIH2ghhMhh4uPjKVq0KI0aNcLHxwdzc3MCAwPZsmULlpaWzJw509QhCpFl3n77bU6cOEGdOnVwcXHh5s2bbNy4kQcPHjBkyJBsm7xGiGeRBFoIIXIYCwsLhg4dyvbt2zl48CDR0dEULFiQbt26MXHixOfOhiZEbta5c2fu3r3LunXrCAsLw9ramnLlyhlLMAqRE0gf6HSYP38+/fv3T3PdnTt3MjwTlRBCCCGEyL2kBToDpkyZkqoUz9OmAhVCCCGEEHmTJNAZ0KpVK/z8/EwdhhBCCCGEMCGpwpFBkZGR2TbFrRBCCCGEyHkkgc6ARo0a4eDggK2tLe3bt+fChQumDkkIIYQQQmQz6cKRDra2tvTr18+YQB85coRvv/2WOnXqcPToUYoWLfrM/YODgwkJCUmxLCIigvPnz1OhQgWsrKyyMnwhhBBCiDwlLi6OGzdu0LBhQ5OMR3vlqnDo9Xri4+PTta2VlRUajSbNdXv37qVBgwYMHjz4ufPZT548mU8++STDsQohhBBCiKdbs2YNHTp0yPbzvnIJ9M6dO2nUqFG6tj1z5gw+Pj5PXV+7dm1CQkK4ePHiM4+TVgv06dOn6d69O2vWrMHT0zNd8QghhBBCCLh48SIdO3bkyJEjVK1aNdvP/8p14fDx8SEgICBd23p4eDxzfdGiRTl37txzj+Pq6oqrq2ua6zw9PSlXrly64hFCCCGEEI+ZqhvsK5dAu7u7069fv0w51uXLl3FxccmUYwkhhBBCiNxBqnCkw5PdLwA2bNjAkSNHaNmypQkiEkIIIYQQpvLKtUC/iDp16lClShX8/PzInz8/R48eZd68eRQtWpT333/f1OEJIYQQQohsJAl0Orz22mv89ddf/P3330RHR+Ph4cGgQYOYNGkSbm5upg5PCCGEEEJkI0mg02Hq1KlMnTrV1GEIIYQQQogcQPpACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAeamDkAIIcSLCwoK4uTJk9y/fz/Fco1Gk+r/pL/Rumge6h9S0KIgZhozXF1dqVKlCg4ODtkXuBBC5GKSQAshRC5w9f5VDpw6wMlzJzlz6QwXb1zk2p1rRMZHgiVgRcq/luCqg1Ix4HwOjv0Hd5IOVgtoCcQCCzCuKFOmDNWqVcPPz49q1apRpUoV8uXLl823VAghcj5JoIUQIocJDg7m5MmTnDhxghMnTrA5fjPBZYMfb2ALlAHrUuATBmcLAo8anM30cGw2lA4Fu4SUx70A7AZ234ZNURBsDzx4vP7cuXOcq32Oxf8shj+AW+Bt702NKjWMiXXlypWxt7fPwlsvhBA5nyTQQghhYmFhYfz555/M/2s+R8KPEBkSCccM69wBz2LQPB5Khaa8FI40bFNgPMQ62GJnboethS2ucbewS0hMdR6vR5eB12GJtgs7fJwZdXgUV65c4fTu3Ww7f5wtJbYaNi5n+HNef57zd8+zcNtC+A24BT4Ffaherboxqa5evTqWlpZZfC8JIUTOoVFKKVMH8SoKDAykfPnynDp1inLlypk6HCFENouIiGDdunWsXrSIi1u2UDQxEecG8FtjIAiYDcWBq+k4lu7AP2hr1nq8YOJEiImBUqUMF3t7OHAAdu+GvXshKgru3gVX18f7VKuG7vo1/i1lxXr3SP7yiOQ/N9CnNdQ8DrgN3ARugWuCK2MGjmHIkCE4Ozu/8H0ihBDpZeo8ShJoEzH1Ay+EyH4RV6+w//0RxB/eR4GICErdB49kDcWxZmD7IThpnZnnN4+aFavhXrIkJCbbqFChx4lx0qVlS3BxSV8QiYlw+jRUrPh4WXg4ODuDXp9i0/h8tlzydWdXMcVS59vsc40jUfuU44aD9o6WAQUG8O7od/Hy8kpfPEII8QJMnUdJAm0ipn7ghRDZ41zQOb5b9x2RG//g17X3sdY9e/vLJ3dRsnz9x1U0AgIMLcWlSkGJEmBjk/lBRkTAwoWwa5ehlTooKM3NLq1dwI7C8Ry8eZCDtw5y5u4pEjXJPkLOA4sN1T7at2/P22+/Tf369VNUBBFCiMxg6jxKEmgTMfUDL4TIOnFxcXy1eCqzzswlyM6QjDrGwK3/gW0i3LGHy07m6IuWoE6zXmi9vB63Jru6gikTTqXg4kVDIp10uXoVrK0hLAysrIzb6UuVIsLFgdNlC/CHfTCnbuRjy7oDj4/lBhVLV2TCiAl069YNCwsLU9wiIUQeZOo8ShJoEzH1Ay+EyHwPwh8wbXIvSm3bTNeLinIjHlW6AIiG9486U6iKP23Gf0OJgiVNGWrGXL8OZ89C8+aPl128CE9209BqiSlblj0aDXPOnmHHG4mEegC7oMi1IowaNYpBgwbh6OiYndELIfIgU+dRkkCbiKkf+OyQkJBAUFAQ9+7dIyIigpDwEMIiw3C0dsSloAvOzs44Ozvj7u6OVvu0jpVC5HyXL51m6bvdafBvIPVuPV7+nj+srOjFm9XfZETnEdjZ2pksxkx35Qp8952hhfrkSUPLdRqKjYEbh4G9hut2dnYMHDiQ0aNHU6pUqWwLVwiRt5g6j5IEOgO2bt3K559/zpEjR9Dr9Xh7ezN+/Hhee+21DB/L1A98Zrp58ya7D+zk4q61hF08AaF3cIiPxiVex/wKcLg4hskdzGDwYWh1AYLDDEUAgoFQCwssKxQhsloJqlduSbPazahUqRJmZjLTvMjZ7uzdy7Z3utPi5B1cYh4v1wPHfdwo8tl3uHZ+3WTxZZvQUNi373Ef6iNHQKcjzq0AXb6pSeuYNsz8fianT582bO8CmvsaOnXsxNtvv02dOnWkn7QQIkNMnUdJAp1OAQEBDBw4kGbNmtG+fXu0Wi3nzp2jcOHCvPPOOxk+nqkf+Jfx4MEDlq1ZxqHN8ylw4SS1w2LxvwPOsam3HdIWfvZ7fH3unzDwWNrH1WlgT1FYcAM2urrRsmVLWrZsSbNmzShQoEDW3BghXsC5c+c426MHHY4fT7E82AaON6xK/W9/w6Zs7npdZ6qoKPjnH0Ni3b07AEopdi5ZQsCc//F7/aNwF9gCXIOaNWvy888/UzF5ZRAhhHgGU+dRkkCnw9WrV/H19WXQoEFMnz49U45p6gc+oxISEti4cSO//fYb69atI75xPGcPQZn7T98n3NqMuS2Ksb2pL442jjjbOdNx5RHKHDyHc0QC1hFRmKXx9PsFGJx0xRKsWoNf/jp80ncKjf0bS0uVMAldfBzfb5nB4jWLOfbrMVorxfpH6/YVMude1y60mxaAWVZUycgL4uOhUSNijx2iZ7sEVvs+Wn4O2ArWkdb8+OOPDBw4UF7jQojnMnUeJTMRpsPs2bPR6XRMmTIFgKioKOzs7PL8m7xSiqPHjrLou3GU+3s/3wUnEJi08grsKGFIoIPttVwsXwyzJo3xqNuUQt7VsChclPzW1owDxiU/aI9k/+t0cP8+3L1L9JXzhKxehPOm3RTu/Bodbt1i27ZtRHlG8dF96LlzP2v+a8q4Qvlp1/Qtxg8YT758+bLpnhCvLL0etW0bIZ9+yrXzB3hnWAIUACxgYzz84e2N9+TJ1H39Feim8bK2bYP9+7EGVi2HLxua817DRFQZwAtij8cy6O1B7Ny5k9mzZ8t04UKIHE1aoNPBz8+PhIQEJk6cyLvvvsutW7dwcnJixIgRfPLJJy/UV9fU35yeJS4ujlXz53P1889pHnWdag8My38A3rezo2vXrnR5vQu23KaWhx92FapmXtktpQwXMzPi4+MZuWQko9/9hXIhhtUx5jC1PszMb83n7f/HsH7D8vwXGWEC9++jAgKImT4d25s3jYurDYajztDyQUu+HfMtZcuWNWGQudD69dCrl6HuNHDSryj+TW8Sav3oYygBOACedz1ZtXgVFSpUMF2sQogczdR5lCTQ6ZA/f360Wi3R0dGMHz+eSpUqsWrVKhYvXszEiRP54osvnrl/cHAwISEhKZZdvHiRjh075qgEOjY2lp+/+oqEL79kcHQ0ydt39cC5WlUptnU3dnbZWElApyP2u6+5v3AuhU9cMi4+VwBGNoHLmlKsmLiCqpWrZl9MIm9SCg4eJPbH7zFbsQLLhMez8sUAK7RatvWrzoT3fsW3tO/TjyOe7dw56NjRUBYPiCtVnLHDSvLTw52Pt4kG8/3m/Nj3R4YMHCJfkoUQqUgCnc30ej3x8fHp2tbKygqNRoNWq0Wv1zNt2jQmTJhgXN+qVSt27drF3bt3n9mdYPLkyXzyySdprssJCbRSit/nzuLad+MYeSEOp2SzBoc55ONYu6pUmzQTBy8TJ/p79hAzeAA2Zy8aFy33hberQk2XzgRMDsDBwcGEAYpcSykSGjXAYtfeFIvPmsOvGgtshw5l5Ecf4ZLe6bLFs0VEQO/esHat4Xq+fAT+byIDNH/y761/H28XCvXi6/HXl3/hkE9e20KIx0ydQL9ydcJ2796NjY1Nui7nzp0DwObRoKDXn+jn+PrrrxMTE8OxY08pK/HI8OHDOXXqVIrLmjVrsuT2ZdTBwwcp2bMk87aP5KMzj5Pne6VKkbBqFY4PQmm0cKfpk2eA+vWxOXka9c03JNhaA9D9NOzYAGvyraLw6MLsOrTLxEGK3CYsNoxhK4YzJ3ofAAlmhi9mTVrDpJHdGHf9Op/88IMkz5nJwQFWr4bJkw3XIyMp9963HOi6mRXdVlDYprBhuRPsdduL24du/LbzN5OFK4QQT3rlBhH6+PgQEBCQrm09PDwAKFSoEBcuXMDNzS3FeldXVwBCQ0OfeRxXV1fjtjlFSEgIfT/ty0a1EXzgGrCtJHjEWGD1xVeU7jvatNMJP42FBZpx47Do0YPEsWMwX/EHUxpo0JspokpE0WhRIz499ikfDP7A1JGKnOrOHfjxR2Ic7ZngGcTs47NJ0CZQsjY8sIGfikCVAq0IeGc2xYoVM3W0eZeZGUyaBJUrQ79+sGQJGkdHujp2pUOZDvxv5//4ePvHJFgmEOscS9/tfTl95jTThk0zdeRCCPHqJdDu7u7069cvQ/tUq1aNCxcucOvWrRQzZ92+fRsgV7VMJSYmMmvS23jNn0FkCwXFDcud9c6Ezv2cxg3fRJMbZgUsXBjz5Svg4EE+9XLlyIw2nFFnUE6KYxs/pPXRnaz8bq3x1wMhCAyE//0PtWgRmvh4Im3hl7GQYGFYfSUGdvs2Ye87cyhdurRpY32VdOgAV69C/vzGRRZaCybWHMnQ2kPp+HVHdul2wU348usvCToYxMyZM7N3LIYQQjzhlevC8SKSZhr89ddfjcv0ej0BAQE4OztTrVo1U4WWITs2b+T7io4M/vJHWt1WzPoLrOLMGF95PHc+vkPXxkNyR/KcXM2alHAuSeBHgfQv3p/Kd2DZn/DG7q3UqledK1eumDpCYUpKwY4d0KYNlC8PAQFoHo2BuFAA3KOA21D/Rn1Ojz3N1p+2SvJsCsmSZ8Awm2GJEjhu28POSTv5tvK3uGxygVhYsGABNWrUeDyroRBCmIISz6XX61WTJk2URqNRgwcPVjNnzlTNmjVTgJozZ84LHfPUqVMKUKdOncrkaFMLCgpSb7fyUxccjEXilAK1vraHun3zXJafPzudbV3HePt+BlW4UCF17lzeuo0iHRITlVq8WKmqVVM853Wg/iiLqj0AxRBUjT411PHjx00drUguMlIpN7fHj9snnyil06nbt2+rRo0aKUABysbBRi1YsMDU0QohTCQ786i0SAKdTpGRkWr06NHK3d1dWVpaqgoVKqiFCxe+8PGy64H/c8US9WMpbYok4qS7uTq49PssPa/J3LunEnx9jLf1O1Dubm7qv1P/mToykZ10OqX39jY+D6LNUTP9UJ5voRiK8uvtpw4cPGDqKMXTrFqllL394/etjh2VCg9XiYmJ6qOPP1I0NTyOWKH69++vHj58aOqIhRDZzNQJtHThSCd7e3u+//577ty5Q1xcHCdPnqRXr16mDuupwsLC+LBrI3zefJ2Rl3UARJvDnwMa4n0ljBqvjTZxhFmkQAHMt+0ALy8AxgDDY+5SeU5ldh7aacrIRFa6fRuuXEEpxZUHV1ixciWfxcYSDHxsB8XGwIgaUDixISdHnOTQb4eoWaOmqaMWT9OpExw8aHwds2YN1KqF9tIlPDt7Qj3AHWgHAQEB1KhRgzNnzpgwYCHEq0YS6DxGKcWCBQsoVb8UHtd24h1uWH64uBV39mygw687sbLO44Nv3N0N0wYXN4yQ/CgC3jmto9ncZhw6esjEwYlMdeoU9O8PJUpwdfBreHzmgddXXnTv0Z1Pr1+nOPA/ZUtj1Z0r46+wc9ZOKpSX2e1yBV9f+PdfaN3acP3MGahenV7XHeno05ESDiWoEVkDMNSD9fPz448//jBhwEKIV4kk0HlIcHAwbdq0oV+/foQGhjLBEwJdYF2/BlS+EE7pWq1MHWL2KVoUtm5FubsDMG0bDPk7keZNm3P48GETBydeilKGL0itWkGFCjB/PiQkUHTrIcwf3EVnowMvyFegAO998gnXr19n2TfLKFG8hKkjFxnl6GiYbOWDR2UpIyLQdujIUsteHBx8kH3r9/Hhhx+i0WiIjo6mV69ez63LL4QQmUES6DxAKcXqVSv5X9nS7N24EQAPdw++bD2D+MMHaRewC3MLKxNHaQKenmi2bkUVKADAoKvwMDSMdu3acffuXdPGJjJGr4e//4YvviCucgVo2hQ2bTKsAv4A6mrhVgQ4HnTk82Gfc+3aNT7++GMKPHr8RS6l1cLUqfDHH2BnB9WqYdWyDa52rpibm/Ppp5+yfv16zFzMiI+Pp1u3bkRERJg6aiFEHicJdC53+PBhOjaqiuOIrnz5IIrpQI8ePThz5gwjeo+gSrEapg7RtMqVQ/P339CsGQe/+IIEICgoiN69e6PX600dnUjL/fuGZPlRuTmAWF0cMV07wfvvY3UyEDD06Z9hCV5ANyDatwK/NfyN4D+DeW/0e1InOK/p0gUOHIBVq+CJ+u5BHkEwAvCDS5cuMWjQIJRSpolTCPFKeOUmUskrYmNj+fCDD7j/7bf8poX8hnGCtHLJR79Zs9A8WVf1VVa1Kvz9N4OUYk9gIAvXLWRLoS2888U7fPvBt6aO7tV27x4cOZLycu0aACfXz2OjRTB//vsn/8b8ywYXHc0j4YojzKsCswqD/YVi9GnWh27dulGhQgU0OXH2TJF5ypdPeV0p4gYPYLNajb6oHtoAD2H58uX4+/szbNgwk4QphMj7JIHOhQ4dOsS4Xr1458IF2gM8Sp53NiuD39KdaJycTBlejqXRaPjs289YUngJI4/rOBH2HRu2N6V149amDu3VotdD9+5w6BBcv/7UzWbMHsAvfo+uWMKINhCmQHfPiW6e3djRbYQkza+6b77Bau58lmo0lGlmwae1E6Ar8DuMGTOGWrVqUaVKFVNHKYTIgySBzkXi4uL4ZMonXFr0BauvQVLPzsSCBUn8+Sf8O3U1aXy5QTGXYvxxuTIdNx8hygKa0pGq5a/j7upu6tBM4v79+1y8eJGIiAiioqKIjIwkMjKSW2HXeRgRSkxkGLFREaj4WKzNzbDWmhHs4YTW2gILSwssLS2pfPku1gqs0WCWmIhZQiJFcMEiESyionC6do05XaoSbAtOygmXRBf6btuGe1hYiliu5ocjheCIh+Hvv4UfrdCBbYgtFV1b8n6n96laqaokzcKgRAmwtUUTHc2UvxOocEtD//aKhz0gPiCe7t27c+TIERwcHEwdqRAij5EEOpfYtm0bgyYMom/CFZZde7xc36MH5jNmYC4DpdKtfedx6P7oiX0CTDueQO2Pa3Np1iXMzPL2kIDY2FiOHz/OwYMHOXjwILsv7eZWgVtQCAiAmnr4EvADntV72P49eKgBEgyXqLlgl/Dscx8rdpzNVYHDwHrQAiWBw+5wpCkc9YD7FsA9MHtghnO0M6VvFadmqZoM7DKQKhWrSNIsUuvWDcqUgY4d4coVugUqfEKgYw+4/AZc/PUigwcPZsmSJfL8EUJkKkmgc7Do6GjWrFnDr7//yvbo7dASVt6HD/6DcBsN2tlzcOw9yNRh5jpmPV4n/sB+tNNn4H8NisdeZeBXAwmYGGDq0DJVfHw8q1atYueenew6v4sLly+gu6yjFFAX6FwIfkzqveIOsbehYTqOa6mDh8mux2mfnUBfAeyT1hc0/JmaPz+enp4U8i6EczFbxherip+3H16eXhQuXDjPf5kRmahiRTh8GHr0gC1bqBAMh36G17vC371h2bxl+M/xZ+jQoaaOVAiRh2iUDFU2icDAQMqXL8+pU6coV65cinVKKZYsWcK4cWOo5BDC3+1A5Xu8/peIRnQfNQeH4l7ZHHUecvcu+hIlMIuNZWdxaNQWjrx5hKqVq5o6spem0+n4OuBrvtj0BdH2EVTRQN07UPcc1L0KHo+2i7Iwx/k9hWc+b0aUGUHdQrUoM2gQqkIFLtnHkGBhhpmVFRorazRWVphZWaOxtiaiRSM01jboEnTExsVifSyQ+IRE4jSgt7BAY2VJiYLeWNvnx8zWFq29PVcjr6JDh4eDBx4OHtjb20uLoMhciYnw/vvw9dcA6DQwvhl8Wxwsl1pyYNcB6Q8tRB7yrDwqO0gLdA6ilGLHjh188cmHeF/7hz0x4BkEfS/Bb5WhrHNZZrefTYPiDUwdau7n5obZiBHwv//hfw0aRkHvr3tzauGpXJvYKaVYvmo5I5ePxMX5HmuOQM2bYJuY9vZ2FpZEvnEEKy+fxwuPHgUgQ3P1lW3y3E2KUCQjRxQi48zN4auvoGpV1IABaGNiiLACCkN8x3i69ejG0UNHpT+0ECJTSAu0iSR9c+revTthYWHE37+Hx8MrVIkLZcAdKBD7eNsLbhYc3hzAaxVfx0wjP21nmrt3USVLoomJMbRCd4Y1TdfQoXUHU0eWYdu2bWPk1yM563UWCoJbJAT9L+U2ie6uaOs3RFO3LtStC5UqgYWFaQIWIisdP45u1Uo6+Bzjrwt/GZb9B93Nu7N0ydJc+yVZCPGYtEC/4oIuLOfbG1D5Pmif+Cpz3dmcW292x++DmXg5OJokvjzNzQ3NiBHwzTeGVugHMPTXobRt0RatVmvq6NLl4MGDjP/oXdzv7+FCGwyj8wCtQyEeNC6Js3cFQ7Jcty7mJUqAJA7iVVC5MtrKlVmeEE2TBU04cOsAXh5w/8By5sxpJP2hhRAvTZozTUzvDtXuPU6e9cB/pR3Y/c0oCgVFUfvLRVhI8px13n0XZWPDQwsoHQpBXkHMmjfL1FE9l1KKKVOmULtfLTrf28Oyo/DzOrDQa/m0/qdc+/Aaztv2wk8/wRtvQMmSkjyLV46thS3re66nik1p1iyFzUcg8H/DOX78uKlDE0LkctICbWK2NesQePE0tyqVxKxBQ6p2eYsKhUqZOqxXh6srmj/+YIvtLebtGgzAe3+9x4CeA3LsVNBKKUa+N5I5V2Yxzwr6HTMsb33diguv76K4T03TBihEDlLAtgCbfD8nX9hraBVMv6jo16UlM49fIF++fM8/gBBCpEFaoE3s264/U+58KM1XHKXpW9/hLMlz9mvdmg4N38TTyhOAh+UfMumbSSYOKm1KKZpOaMpc/SyWn4R+JwzLw0t44HbsvCTPQqTBtW13bv+1lIcWhg+9/12+y8TevZEhQEKIFyUJtBAYpvn+5bVfDFcsYPqJ6QQFBZk2qCfo9XqGDRvGwcPbWb8SOp81LE+oXJH8/55AU6yYaQMUIgcr3fw1LOb8CoAL0OPPP5k7e7ZpgxJC5FqSQAvxiL9TZX45UJg1SyDRN5EJUyaYOiQjnU7HgAEDWD5nDlt2QrPLj1bUr4/Fzt3g4mLK8ITIFSz7DyCmVy8A6gMR745k5Z6Vpg1KCJErSQItRJIvv+TNTbfocA78f4CFcxZy+vRpU0dFQkICvXr1YsmCBewAaif96tyqFWzaBPnzmzI8IXIVm19+IcrT0F1r3EM9C7/rzpnbZ0wclRAit5EEWogkY8eCjQ0AkyMNXSYmTDBtK3RcXByterdiWegy4s1hU6FChhWvvQZr1oCtrUnjEyLXsbHBfsMGoq0MNR9nbNUzcPIb0h9aCJEhkkALkcTVFUaMAKAh4A+sX7+ey5cvP2uvLJOQkECHzh3Y5rgN6oDNKBv6HD8M69bBokVgaWmSuITI9by8sPptMZfym9OyKPzzy1Hmzp1r6qiEELmIJNBCJPfuu49boQGqwdTFU00SysyZMwnZvhmPcMP1oY2H4uHiAW3bQi6Z6EWInErbvTv5Tl3hfqgHAG+99RYnTpwwcVRCiNxCEmghknuyFbo8LLm/JNt/3r137x5bP/yQnbGwc5UlP/lNZVrLadkagxB5nWuRIixZsgQzMzPi4uLo3bUrkZGRpg5LCJELSAItxJOSt0LvhFjbWHYc3pGtIazo25c/Hj4kH+CVkMDQaG8stdJlQ4jM1rBhQ9775D3eKAqbb17k49d7SH9oIcRzSQItxJOSt0JfA/9psG3Ntmw7/Y0vvmDQhg1YAzqNBn7/Hbp1y7bzC/Gqqeplzu83wCMW+v+9gYBZs0wdkhAih5MEWoi0PGqFvmJjg5kelizJnm4cau5cir7/PuZAjBYCv5mK5lHdWiFE1vBvP4qv6mkAqJgAmtGjuX37tomjEkLkZJJAC5EWV1f45x/Wf/4524ErV65w8ODBrD1nQgKJb70FQIQltO9lQYE3+2btOYUQONs4s3tQc3Y/msyzv07HjammGTwshMgdJIEW4mkqVaJbjx6YmZmBE/yy9JcsPV380aNYxMYCMLYl9Bj6E4UdCmfpOYUQBl0qvMbrXSH4UWn1Kr/8AqdOmTYoIUSOJQm0EM+Qv0B+rN61gtGw5vxSdDpdlp1r3pfDjP+HlyrHwNoDs+xcQoiUOvp0JMTRgp5dQA9YJiYaxh5ERZk6tCwXExPDpUuXuH//fpa+xwmRl5ibOgAhcjIbCxsaRDkw5M8YHMOj2blzJ02aNMn08/xz4R8W5D9GfA2oEKThp0nZN2hRCAFONk40K92MDfoNfFIHPtkPnD0LQ4bAwoWg0Zg6xEwXFRXF9OnT+frrrwkPDzcuz58/P05OTjg7O6f4m9ay5H/t7e3R5MH7SYi0SAItxHNMuOFBo7N3ARiyYE6mJ9Dxung6LOhASAk4UBw+L/MZjRzdMvUcQojn6+bbjQ0XNjC1KdQ/Bk1jgMhIiIsDa2tTh5dp4uLi+Pnnn5k6dSrBwcGp1oeHhxMeHs7Vq1effSBLwB7QAPfB3NzcmGg7FHTAJb8Lzs7O1K1bl9deew1nZ+csuDVCmIYk0EI8h2P7rrD2OADXD64lLi4OKyurTDv+t5u/JcQiBIDit4ozcdLETDu2ECL9OpTpgBYtOjMdPevBCb/38PjsszzT+qzT6Vi4cCGTJk3i2rVrxuUVKlRg6NChJCQkcO/BPW6H3SYoKojgmGBCE0IJ14XzUPOQWPNYlJ0yJM35MCTQAJeB3yAxMZGQkBBC4kKgO3AB2AaLFi1i9OjRtG3blj59+tC6dWssLaWuvcjdJIEW4jnKNu8JfAhAOW0cmzZtokOHDplybKUUP677EdyBOFg1cpX8BCqEiTjZOFHLpRb7QvYRUhl2Va5IjzzwelRKsXr1aiZ8PoGL9y4aWoyB0qVLM2XKFMo3Lk+fNX0ISggi2CwY5azgKY3F/lfgrS2wpTSs9oG7+cDD24O3v36bBw8eEBoayn7Nfk6an4Sy4HzKmQdBD0hISGD16tWsXr0aZ2dnevToQZ8+fahRo4a854lcSQYRCvEc1kVLcj+/obWkqh4WL1mcacfesGEDt7W36XgGft5mT9XT5yAxMdOOL4TImN7Vehv+sYMNZzY8XhEfDwEBkEtmKQw4FsDe63vZunUrNWvWpEuXLlwsfRH6g7aDlp9++okzZ87Qs2dPbCxsOHH3BHcf3kWR9u1zsHKgTIEyWNSpR7WY/Pz0F9z5VkPQn97siunIO9268fnnn/PTTz/xywe/MKDyAGoUrkHwiWB27NjBgAEDsHewh8HwwO8Bs1bPolatWvj4+PDZZ5+laBEXIjfQKJmz1CQCAwMpX748p06doly5cqYORzzHpbq+lN5/hrMFoLKVFffO3cPe3v6lj9vh9Q6s9VnL/NXQ9wRQoACEhOSZn4yFyG3CYsNw+twJtFAqtBSXvr8Et25Bly5w8CB89x2MGWPqMJ/pRvgNSk4vied5T84tOmdcbtXJirhKcThaORI6MdS4PCo+itdXvo67nTse+Txwt3fH3d4dD3vD/252rtha2hk2Vgq+/RbeeSf1if38DPdT587g7Y1SKkXr8trAtXT4I9mvdyHASeA/IMwwrXqfPn3o2rUrDg4OmXqfiLzH1HmUdOEQIh3sazaA/Wfwvg8WVeNYu3YtPXv2fKljKqXYdWUX+ED1W48W+vlJ8iyECTlaO+IS6UKIYwjXbK+RqE/E3MYGgoIMG7z7LtSqZbjkQIGBgXSa2Qmdm45zFufACezi7Rg7diwt+7YkWhONu717in3sLe1Z9/q61AfT6WDKFLh48XElEo0Gxo6FunVh5UrD5coVw/aHDxsu770H+/ejqV07xeEc7BxoWLwhu67tMixwAZo8ulyHXf/tYtdbuxgxYgQdO3akT58+NGvWDHNzSVVEziNdOIRIB9f6LQDDC6ayo2Fq75d1/vx5wvOFYx8HPvceLfTze+njCiFejp+14XWos9Kx5eIWcHaG5cvBwsLQxap7d7h/38RRpnTlyhX69OlD+SrlueBwAQCNVsOo3qO4dOkSn376KXU969KsdDMquFV4/gFv34YmTQwJ9OLF8EuyiaTMzAxfIL7+Gi5dgqNH4cMPoWxZw/qCBaF69cfbKwWff45/iB07++7g2phrTGsyjQquyeIoBrQBxkFs51iWnlpK6w6tKVKkCG+//TbHjx9HfjAXOYkk0EKkgyZZYlvVDDZt2UTso1kDX9Tu3buhGFS9k+yFmPxDRwhhEm0828CjoQjzDs4z/FOjBvzvf4b/b9yA3r1BrzdNgMkEBQUxcuRIypQpw++//w4VABvDupm9ZzJ9+nTc3DJYFnPzZqhcGXY9ain29YV69dLeVqOBKlXg00/h9GnDJSAAkrcaHz8OH3xguA+LF6fYpG+ZoK/NycHHODn0JBPqTqCoQ1HDtlrAG+gKvAN369zlu7XfUaVaFSpVqsQ333zD7du3M3Z7hMgCkkALkR5FihDYrw19O8Ka8pDolsjly5df6pA79uyAQuCX/LNAWqCFMDm/8n5wCQgG++hkYx1GjjTMTgiwcSN8+aVJ4gMIDQ3l/fffp3Tp0sycOZOEhAQAHJob+g4XdSjKoPqDMnbQxERD94uWLQ1jMQAGDIBDhwxJdHqULQtt26ZcduDA4/9v3IDp06FhQyhcmAqTZjJN35irIy6ws+9OBlUdhKO1o2FbK6AS0Bt4G/678x/vvvsuRYsWpUWLFixatIiHDx9m7DYKkUkkgRYiPTQaCvw4l98qw3VHoDScO3fuOTs9nVKK3Xt3w1ZocPXRBA0eHlC4cGZEK4R4Cb6+vvAHMAtKBZd6vEKjgblzwdPTcP3DDx+30maThw8fMm3aNEqVKsUXX3xBdHQ0AM2bN2fmxplEWEcAMKL6CMzNMtB3+MYN8PeHadMM1+3s4Pff4ddfwdb25YIeNgxu3oQff4RGjQxdQADu3oU5c6BFC8xq1KRhiYb83O5ngsYFsfq11XT17YqV1lBz39beFu+C3gDo9Xr+/vtv3hj5Bu7u7nzwwQeSSItsJwm0EOnkbu9O+YLlDVdKG/owv6hr165x68otOAAN7z8a3S7dN4TIEfLly0dRd0OXgtOnT6dc6eAAf/wBVlaGLhw9e0JMTLbEdfToUTw9PXnvvfcICwsDoFatWuzYsYPNmzezKXQTADbmNgyqloHW55gYQ5/mffsM1ytUMAwGfOONzAu+cGFDC/727YYBmXPnQqtWhn7lYEjeH7Eyt6KjzosVqht3h15iXvt5fNXqK86ePsu///7LyJEjye+dH8ZAVI8oPp//OWXKlGHp0qXST1pkG0mghciAVt6tDP8Ugv8u/vfCx9m9ezcAToBj0mAk6b4hRI7h+6jLQmBgYOqVlSrBjBng4gLz54ONTbbENGPGDIIeVQMpX748f/75J/v378ff359LDy6x/vx6AN6o+AbONhmYNtvGxtBHGWDIEEO5Ph+fzA7/MRcXGDgQNmyA4GBDhY8330y5zS+/wGuvkb+oJ/0/+ZMR5xzQhIVRvXp1fvzxRwb98OgLQhEgGm7dusXrr79Ow4YNOXHiRNbFLsQjUhtGiPSKiWHE0su02g6/VYJ/zA6/8KGSEuh8dnboxo9He/QoNGiQWZEKIV5SuXLl2LxlM2cSz9BndR/6V+5Po5KNHm8wcKCh3rFzBhLVl3T4sOE9p379+uzYsQOtVmtcN+PfGcZJUEbVHPX8gyUkPG79BUM3iwoVoH79TI35uRwdoVev1Mu3bjX8jY2FP/80XMzNDZVBOndmSP1OWFlYcenBJZr9rxkTJ04kJCSEPfv2UPnzyrTzaEfARwEUKFAgW2+OeHVIC7QQ6WVtTbFV22h0FerthnuH7j13l6fZcnQL9ADrboUJHjsI1qwxDKoRQuQIvr6+YAGJXRL5/eTvLDy5MOUGGk3q5DkLq3JER0cbu5PUqVMnRfIcGRfJvOOGaiFNSjahvGv5Zx9szRpDP+7k3dA0muxPnp/l8GFYtw769398PycmGiqEDBmCZ/n6TL1YjCVdlzBgwADOnz/P2LFj0dTWgC+sc1pHoRGF+GrmV+h0OtPeFpEnSQItRHppNGiqVgWgWijcv3mf0NDQ5+yU2u3bt7muuQ4+cL7Eee7H5Kx6skKIRwl0HHARLDWWKWbUS9POnYbSbzdvZkk8J06cMCaCfk9095p/fD4RcYbBg6Nrjn76QeLiDLModuoE16/Da68ZktKcyNraUM1j3jxDn+ktW2DoUEgqyafXp+j25ujoyLfNm/NtnxFodYYvF/Fl45lwfQIlO5VkVzYP9hR5nyTQQmREtWoAlMNQYelFBhLu2bPH8MoLB3tze3xd0lkeSgiRbZL6QLMRJphNYG77uU/f+PRpQ9eC//4zJKWPSsplpiNHjhj/r/bofQhAr/T8+O+PAJRyKkVrr9ZpH+DyZcPsgdOnG67nz2+oIpIbZvmzsICmTeGnnwzTqu/ZAx9/bKg/nSQxEXr3ZsyQGUSv9mb2P4XwDQZs4Ea1G/jP86dd33bczKIvOOLVIwm0EBnxqAXaAigPnD13NsOH2LVrFxyGYrOsuHuiGWbTfzCUkBJC5Bj58+encOHCEA6Xzlx69sa+vo8Hwe3fD++/n+nxJCXQTk5OlChRwrh844WNXHhgmHnwrRpvoTXTpt75jz8MyWZSEl69umH2wC5dMj3OLKfVGiZ1+eQTQ7eTJAcPwj1DtzrLU2cYsvk2gbPg3AwzPtsKVa1hfdH1lOpbiimfTXnpibCEkARaiIx4lEADVK0L0y5Ny/AhkgYQ9vL2xnbFahg7Fi5cyLQQhRCZ45mVOJ40fbqhCwfAN98YBr1loqQE2s/PL0V3kh/+/QEAe0t7+lfun3Kn2FgYMcIw+UuEoYsHY8bA3r1QqhR5Sp06cOqUIbGuWNG42Pu+nvf3wpGf4coMsKuewKTbk/Bs5Mm6deuk7J14YZJAC5ERnp5gb5iZrFosXFQXSdCl/+fae/fuGT+Mmzs5PV6R7CdZIUTOUK5cOQDOnj3L1QdX+e6f7/jnxj9pb2xtDStWGOpEA/TrB1euZEoc0dHRxveN5N03ouKjuBxqmBG1f+X+5LfOn3LHiRNh1izD/46OhsGD330HlpaZEleOotFAuXKGrh0nThgaJb780jB9+COxNhaE2QCucKvFLdrPaU+zds1ealIs8eqSBFqIjDAzM/a7q3odnC47EZ0Qne7d9+zZA6WA4uAb+2jyhTJlDP0RhRA5SlILdJwujrKzyvL232/z89Gfn76Dp6dh0BtAWJih5Tcu7qXjOHHiBPpHFT6SJ9D2lvacHXGWda+vY0ytMal3/OADKFTIMEnK8ePQocNLx5JreHrC+PGGrh3Xr8P06Xh+/hOf+n+KOeaggcLecMN1G75dfXnn3XeISGqlFyIdJIEWIqMefYBVDIG4VVHks8yX7l13794NTYB+YH7mmGGhTKAiRI5kHEiYCJXsKwGw5uwa4nXxT9+pSxcY/agSxpEjMG7cS8fxtAGEAFozLW2921LKqZRhRsHkAxhdXAzVQXbvhuLFXzqOXKtoURg1CvP+A/mw4YecHnmaBh71mL8Gjv0Ow4rq+fbi/yhTpgy//fab8cuKEM8iCbQQGdWxIwfbtqUjEB0Tw61bt9K96459O8ADioWDc+Sj8lEyhbcQOZIxgQZKRJUAICw2jK2Xtz57x6++gpo1Df8fOADR6f+VKi1JE6g4OzunGECYwtmzhu4KH3+ccrmXV8oJUwReBbzY6fMZTa5qsE2EGRth03YwCwqib9++1K1b13ifC/E0kkALkVENGxI5ejSbgEQMpezSMxAlPDyck/dPghn43U62QlqghciRnJyc8PDwMFw5D1ZaKwCWBy5/9o6WlrB8ObzzDuzbB7a2LxVHUgt0tWrV0Gg0JOoTGbVxFP/e+tewwW+/GX4ZO3UKpk0zDBIUz6Rp0ADNjh3oihUFoHkEnNJo6A4cOHCA6g2r8+abbxIcHGzaQEWOJQm0EC/A29vb8I8HjD80nk93f/rcffbt24cqYki0qyc1WifrUy2EyHmSWqEvBF4w1lh+bjcOgGLF4Ouvwcrqpc6ffAbCpO4ba8+t5cd/f6TRrJpc6eQPffsaWrk1GkMJvVq1Xuqcr4yGDdH+d8ow4BNwUoplwGJrDY5D4Nfrv+Ll68X06dNJyILa3iJ3kwRaiBdQpEgRbGxsoA0cjTvK9IPTeRj/8Jn77N69G4oZ/q8fYmP4p1y5l26dEkJknaRKHGfOnKFr2a4AhMeFs+XSlowdSCn4++8Mnz/5AMKkGQjP3jtLxXtaDs81o+SaRzPsubjApk3w2We5Y3KUnMLBAQICYNUqKFgQgNdjFf/NhRrFIMI6gjFjxtCgQQOuX79u4mBFTiIJtBAvwGznTlZYWXHkd9Dq4EHMA+YefcZMZcDOXTvB8GshgW2qw+DB0L171gcrhHhhSS3QMTExVLCqgLW5NQDLTz+nG0dy4eHQtSu0aAHz52fo/Mn74larVg2U4v0L7hz71YKywY8Gu/n7G6psNG+eoWOLZDp1Mswk2bYtAI6JWnxLNMbTxhMwdOuoUqUKGzZsMGWUIgeRBFqIF3HtGm3CwqgaBxWuGVp7vj3w7VNrQj98+JDD1w8b5v8GtG/0gTlzDFPpCiFyrOQDCa9duJaiG0dcYjpL1D18+Lhf8vDhhr7K6ZTU/9nZ2ZnixYsbahwPHIhZTKyhy8akSbB1q6FcnXg57u6wdi38/DO2c39j1oT1/PfffwwfPhyABwUf0OajNkx8byKJiYkmDlaYmiTQQryIZDMSVtxteCO9Hn6dpaeWprn5gQMH0BXWGa/XK1Yva+MTQmSK5Al0YGAg3X0NvxpFxEWw5XI6u3EUKgSLFxsS3pgYQ2t0ZGS6dn1yACGVK8O77xqSva1bYfJkw/TWInNoNDBoEGY9e2JjYYO1tTUzZ87k4Ovd+LggaNvAl0FfUrtdbW7fvv3844k8SxJoIV6Ery+6R6WhqlyDgtaGvnNf7f8qzYocu3btMvZ/LmBdAO8C3tkWqhDixRUoUAA3NzcATp8+TRvvNo+7cTyvGkdyTZoYkl2Ac+cMXbieU73n4cOHnA4MxB1D/+e5R+fyw8EfiPx4oqElunHjF7hFIsOOH8fvj9V8sh/2zQMvezjsdxivQV5s+Fu6dLyqJIFOB39/fzQaTZoXC6mv+WqysCDG09A3rpqC1k6Gn3VPBZ9iw4XUb6i7d+829n/+dbcjmvffhx07si1cIcSLS2qFPn36NPaW9rTxagPAn+f+TH83DjB02Urqp7x0qaEb1zOc2r+fRUpxCKjuWZJJOycxetNo2qzoCK6uL3BLxAvRajHzLgNAzVtwfDYMPwrR1aNps7YNAz4ZgE6ne85BRF4jCXQ6fPDBB/z+++8pLrNnzwaguQzaeGWZP5oooQpQ+kFJ44yEX+77MsV2cXFx/BP4DziCmR5abb9uqNW6YkU2RyyEeBHJE2i9Xk/3co+7cfx9KQOVNczMYOFCKFzYcH30aDh6NO1tjx6lTK9e9ACKAGWWz+V2pKHLQP/K/V/wlogXUqECHD4Mb78NGsPkKzM3wMaF4GEFAQRQbGQxzt84b+pIRTaSBDodmjVrxhtvvJHiYmdnB0CvXr1MHJ0wFes6dQCwB+IPn2Ko31AA9lzfwz83/jFu9++//xLvbqgZWzYELGMfDTSUCVSEyBWSStlFR0dz/fp12ni1wcbcUIoyQ9U4wFBubulSQ7/l+HhDf+iHyUpgKgUzZkDt2jiGhACw1cKCt1sYxlq42LrweoXXX/5GiYyxtob//Q+2bzfU+AZaXoJTs6D7KbjtfpuyM8oyacWkdE2sJXI/SaBf0OLFi7Gzs6NDhw6mDkWYSrKBhJanTjGm1hgstZZAylbo7du3G7tv1AlKVp9VpvAWIldIPpDw9OnT2Fna0cb7UTeOs38SmxibsQPWq2f4Fcra2tCt41GDDGFhhoT6rbcgPp5EYDzwYdvKbI4wtFQPqTbE2AdbmIC/P5w8CX36AOAcC8v+gJo3QG+rZ8rpKfhM9eHyg8umjVNkOUmgX0BISAhbtmyhY8eOxpZo8QoqX55EM8NLyOXGDQrlK0Tvir0BQ9/IMyFn0Ov1/Pbbb8YBhK3CXAz/2NhA2bKmiFoIkUFPVuIAjNU4KrlXIigqKOMHHTfOUM5uwADD9X//NcxKumoVAPoiRfDXaPgaiCxvaKE2NzNnWPVhL35DRObInx8WLIA//oACBVC9e1O+wkB49EPCef15vL/35rvd35k2TpGlZLqiF7Bs2TISExPT3X0jODiYkEc/xSW5ePFiVoQmspOVFfsaN2bh1q3siY6mZ3g479Z5l3nH5qFQfL3/a3rZ9+LyrctgGMSPX1LVo6pVZbYwIXIJFxcXXFxcCAkJMU6r3da7LTfH3qSwQ+EXO6hGA6VLG/5XCkaOhKtXDdfbtePwsGHsa90a8sE57TkAuvl2o1A+qfecY3TpAnXrorGxYW7+/PT+tzftfmyHq1Mktxx0vL3jbQ5ePsiiPovQmkmpwbzmlWuB1uv1xMbGpuvytH5MixcvxsXFhWbNmqXrnLNmzaJ8+fIpLh07dszEWyVM5f6wYcwFzgHnz5+nTMEydPTpCMDCkwuZHjAdEsFiiQUfVhtP4Sv3DDtK/2chcpXkAwkBbCxsXjx5fpJGA7//Dk5O8O238OefHLhwwbDOD3QYKjyMrjk6c84nMo+7u6FFGmhYoyE3vjzLzvnWHP0Jqt2CZYHLmPrdVOkXnQe9cgn07t27sbGxSdfl3Llzqfa/fPky//zzD6+99hrm6WxBHD58OKdOnUpxWbNmTSbfMmEK3t6P6zmfP28YgT2h7gQAEvQJ/HX/L9BB9xrd+dSjO2bxjwYQSv9nIXKV5Al0liRDZcrAlSswdixoNIYJVMxBU10DQM3CNalZpGbmn1dkqvyLFlEkMpayD+CfufDhHPj0ncl06dKFsLAwU4cnMtEr9xuyj48PAQEB6drWw8Mj1bLFixcDGau+4erqiqvU7MyTSpcujUajQSll/MJVs0hNGhZvyK5ru9Dn1wMwaNAgQxmkJNICLUSuklSJIyoqihs3blDsUSWGiLgI1p1bR/DDYMbWHvtyJ3nUkgmPZiAsD8rWkKxL63MuMXYsJCTApElYJCbyaSS0AfqsXk3V48eZs3AOzeqk79drkbNplPyukCG+vr7Ex8e/dB/mwMBAypcvz6lTp4xvzCIXUoofnZwoEx7O/QoVeP3kSQAO3zpMt+7duLr/Kl5eXpw7dw7NiROGAUL//QcrVxpqwgohcoUdO3bQ+NHMfxs3bqRly5YA9FzZkyWnlmBvaU/wO8HYWNi89LkePnxIPod8qEEKPKBQvkJcHX0VC61M3JVrHD0KvXvDoy4/D4F3XGH2QOjv0p9fx/xqmJpdvDBT51HyCZ4Bx44d48yZM/Ts2dPUoYicQqOhrU5Hc6BE0gAgIP5aPFf3X4VCcO+Ne/Re3ZvzRW1hyhRYvVqSZyFymSdL2SXpUrYLABZmFpy9dzZTznX8+HFUUUPyDDDMb5gkz7lN1apw5IihRRqwA34Khg0rYN21ADr37kxkZKRpYxQvRT7FM2DRokWATJ4iUgp+NKuYd2QkSm/osvHLL78AYFbcjFAVyqL/FmFu9sr1mBIiz3B1daVAgQLA41J2AK28WrGp1ybuvnOXKh5VMuVcR44cgUfdnS3NLBlSbUimHFdkM2trw6DQbdugqGEyAK97ELMR1ixag5+fHycf/Wopch9JoNNJr9ezdOlSqlatSpkyZUwdjshB4sqXB6AAcPfwYcLDw1m2bBkA9UvX580qb+JT0IeSjiVNGKUQ4mVoNJpUlTgAbC1saeHZIlNbiHed2AU+hv97VuyJi51Lph1bmEDjxobJV/r2xXr+cuqVbAEYBp7XrFmTub/OlSoduVCmN4lNmTIlU47z8ccfZ8pxMouZmRk3b940dRgiB7KuU8fQpxkI2bSJvUeOEBMTA0CH9h0IuBXA7Eh/NAMHGqpvDBokNaCFyIV8fX3Zs2ePsRJHVvVh3fVwl7F5a0zNMVlyDpHNHB1h/nyKABvadeGLL77g448/5mNtLEvWDmLtgbUs+X6JTM6Wi2T6p3h6K1w8i0ajyXEJtBBP49aiBfpx4zAD4g8cYO7duwAUK1aMUZ1HMVY71jCYZOFCWLcOhg41bcBCiBeSNFApIiKCW7duUaRIEeM6pRSBIYGcDjlN93LdX/gcDx8+5MGaB3AZKrWsRCX3Si8dt8hZzMzM+OCDD2gWH0WNKdN4by1877eOKu18+HPm35SVWWpzhUxPoK9cuZLZhxQiRyvi48N5jQYfpYg/eJCjDx4AMGDAALTaR7NPJZWw8/MzTJoghMh1nhxImDyBfufvd/j2wLfYWdjRzrvdC1fjOH78OCpOwSH4eKI0JOVl1Uv5EG9lgWVcAmMOQ/OCN+n3RkW++HqzseKLyLmkD7QQL0mr1XLJwQGAko+SZzMzMwYMGGDYICICkiblkQlUhMi1nlaJA6BesXoAPEx4yMaLG1/4HEeOHDH+7yf14vM0Td++WP4XyP1Khgm5fO/B3uOJbB/RlFUrl5s4OvE8WZJA58uXj9q1azN06FB++ukn9u3bJ+VaRJ4W8mhSBXcMladatmxJ0Uejrjl6FJIGiMgHohC5lru7O46OjkDKShwALT1bYm9pD8DywBdLfs7eO8uuY7sAKFiw4OP3EJF3eXlR4HAgt8cPJ8EMLPQw9azCY9hrLJr2iamjE8+QJSOZGjVqxNGjRzl48CCAcaBFiRIlqFSpEhUrVqRSpUpUqlSJUqVKZUUIQmSr8OrVmfLffxwBwoE333zz8cpDhx7/Lwm0ELmWRqOhXLly7Nu3L1ULtI2FDe2827Hk1BLWnV9HdEI0tha2GTr+uL/HsaHYBugA1WKryUQbrwpzcwp9OZOQNu0I79YBz+B4aoeAz+TJfKqL5KMPvjF1hCINWdICvXbtWm7evMmdO3dYv349Y8aMQSlFREQEa9euZcqUKXTt2hUvLy/y589PvXr1siIMIbJN/vr1mQSsBfK5udG2bdvHK5P6PxcqZLgIIXKt5KXsniw9ljR4MDohmo0XMtaN43r4dTZc2GD4VE6EatWqZUq8IvdwadCSQhdus6KRoW/913Xh48j/0W9yPylzlwNlaR9oNzc3WrduzfDhwwGYM2cOUVFRHDlyhAULFjB+/HgaNGjAnTt3sjIMIbKct7e38f9+/fphYZGsJmxSC7T0fxYi10tKoMPCwlJ9dqXoxnE6Y904iuUvxtzqc+E48K/0f35V2ToUoMu2a3w+uR1f1QVsYIF+AU3HN0V/44apwxPJZHsxWmtra6pUqUKVKpkzY5MQOUH16tWpU6cOISEhjBo16vGK+/chqTKNfCAKkesllbIDQyt0oWS/KlmbW9O+THsW/7eY9efXZ7gbR9SFKFhj+F9aoF9dZhoz3p+0FvPNXzFh3wTQwnnddqI8S2DXtQfaWbMgf35Th/nKkyocQmQCCwsL9v7wA+fr16dQmzbwqBY05uYwYwb07QuNGpk2SCHES3tWJQ6A7r6Pu3FsuLAhQ8dOqsAhAwgFwPgW41nVYRXaWC3fbgaHeD3axYvRlS8H27ebOrxXniTQQmQSTWgozJsHx4/DsWOGhfnzw4gRMH8+1K1ryvCEEJmgUKFCODwqW/lkJQ6AFp4tyGeZD0hfNQ690rPhwgZ0ep0xgfbz85MBhAKATtU6ceDNA4yvb8VGT8My7c1b0KQJvP02xMaaNsBXWJYk0EuXLuXixYtZcWghcq7k3ZKS1XIVQuQdSZU4IO0W6KRuHADrz6/nYfzDZx7vr/N/0WZxG0pPL83pUMPxpPuGSM6vtB+7P71Ar1aODPUC4zPqu++gWrXHDTYiW2VJAt2zZ0/KlCmDk5MTjRs3ZurUqWg0Gm7fvk1CQkJWnFII0ytQAIoXN/x/9KhpYxFCZJmkbhyBgYFpVkdIqsYRkxjz3G4cP/z7AwD3Ht6D+4ZlkkCLJxV1Kcr1aTe4XbYdlYEDSStOn4aaNWH6dNMF94rKkgR606ZNTJs2jbZt23Lv3j0WL16MUopRo0ZhZ2dH2bJl6dq1K5MmTWL58uVp/gwmRK6U9MF39CgEBUGbNjB5MsgvMkLkGUkJdGhoKMHBwanWNy/dHAcrQzePZ1XjCAwOZOvlrQD4mftBvGG5JNAiLfa29qxauYo6ffpQD/jIAhLMgIQE1Lx5EBNj6hBfKVlShaN58+Y0b97ceD0hIYHTp09z8uRJTpw4wcmTJ9m7dy+rVq0CDD+J6XS6rAhFiOxVtSqsWgVXr8KmTbBhg+Hi7w+enqaOTgiRCZJX4ggMDMTNzS3F+qRuHAtPLuSv83/xMP4hdpZ2qY7zw0FD67MGDY7nHAFwcXGRAYTiqczNzQkICMDR0ZGp939gky188zecGd+ZoTY2pg7vlZItZewsLCyMMw/27t3buPzu3bscP36c//77LzvCECLrVa36+P+ffzb81WhSLhdC5GpPVuJo3Lhxqm26+3Zn4cmFxCTG8NeFv4zdOpI8iHnA7yd/B6Ctd1vOLzsPGFqfZQCheBYzMzO+//57dFN0zAydiX9lyDf8f/gWbUSDBg0MG+n1YCZ1IrKSSe9dNzc3WrRowTvvvGPKMITIPMkT5X/+MfwtUwYejdoXQuR+RYoUIV8+Q6WNtAYSwuNuHJZaSy6HXk61/pcjvxCTaPjJfXDFwZw9exaQ7hsifTQaDTMmzeCzEp/BeoiMiKRFixb89ddfqK+/hi5dIC7O1GHmadk+kYoQeZqbGxQuDLduPV4mE6gIkadoNBp8fX05ePDgU8fwWJlb8VfPv6jgWoH81iknvUjUJzLz0EwAyruWJ/+D/MbBiDIDociI98e8T9ECRenfvz+xsbH81qk9bRL0hpVt28KaNWCXuvuQeHmZnkBPmTIlU47z8ccfZ8pxhMh2774LN2/CN98YrssU3kLkOUkJ9NNaoAHqFauX5vI1Z9dwI8IwLfOoGqM4sudx2UtpgRYZ1bt3bxwdHenWrRubWsax5xjUvw5s3QrNmsFff4GTk6nDzHMyPYEOCAh46WNoNBpJoEXuNXo0/PHH4+uSQAuR5yT1g7537x4hISG4uLike9/pBw0lx5xtnOlVsRdDvh8CGAYQFilSJPODFXleu3bt2LRpE63GtaJFj1hWrYSWl4B//kE18kez+W/DL6Qi02R6An3lypXMPqQQuc/hw4a/Wi1UqmTaWIQQmS75QMLAwED8/f2fuu2dyDusPLOSvpX6cuHBBfZe3wvAoKqDsLWwNc5AKAMIxcvw9/dn7897adKnCR3ahrPob+h6BjQnTqKrXw/t1m1QrJipw8wzZIimEFnh0CHD3/LlwdbWtLEIITJd8lJ2z+rGse/6Pgp/W5i3Nr7FXxf+MrY+azVaRlQfQWRkpHEAofR/Fi+rWrVqHFh5AOeN7vSoB/MqG5ZrL1wksU4tOH/epPHlJZJAC5EVKlSAEiWgRg1TRyKEyAJFixbF7tHgrGcl0NULVzdOqrL4v8UsPbUUgM5lO1M0f1GOHz9uHEAo/Z9FZvDx8eHA3wcosbc0b5aB72salpvfukPID9NMG1weYpIqHGFhYfz5558EBgYSFRWFvb095cqVo0OHDjg6OpoiJCEyV8eOhjqcEyaYOhIhRBYwMzOjbNmyHD58+Jmz6VpqLZnddjYlHUuy6eIm1p1fB8DomqMBjN03QBJokXmKFy/O/p37ad6iOWMLnCCsIfjcgzFF17Ej5AxlXcqaOsRcL9tboLdv346npycLFiwgNjYWZ2dnYmNjWbBgAV5eXuzYsSO7QxIi8/n7ww8/GEraCSHypKRuHM9qgQboUb4HVTyqMPvIbACqelSlTtE6wOME2tXVVQYQikzl6urKrp27qB9fn0/ioWcXuBt7D//5/py/fx4e/fIhXky2t0CPGDGCgIAA2rVrl2rd+vXrGTZsmLE/mBBCCJFTJQ0kDA4O5t69exQsWPCp28YlxtGvUj9+Pvozo2uONg4WPPxowLEMIBRZIX/+/GzatIlu3bqxYfsGaArB0cE0/bkBp/dVwn7gUOjUydRh5krZ3gJ948YNmjdvnua6pk2bcuPGjWyOSAghhMi4J6f0fpZ8Vvn4oukX3Bh7gx7lewAQGRnJuXPnAOm+IbKOra0ta9asoWexnvDoR/6pS+9i/9ffqG7d4LffTBtgLpXtCbS/vz8jR47kzp07KZbfuXOHUaNGPbMUkBBCCJFTpLcSR3K2FrZYai0BZAChyDYWFhb8/vvvDCs3DHbDkgoQbQ4anQ769oWZM00dYq6T7Qn0/PnzCQ8Pp0SJEhQsWJBSpUpRsGBBSpYsSXh4OAsWLMjukIQQQogMK168ODY2NkD6E+jkkrpvgJSwE1nPzMyMmTNm8kGdD9h0F1r0hnDLRytHjoTPP5d+0RmQ7X2gCxYsyPLly3n48CEXLlwwVuHw8vIylgQSQgghcrqkShxHjx59oQQ6+QDCwjLgWGQDjUbD1E+n4vQ/J95Z9g6NFGzRaCigFHzwAYSFwZdfgvTHfy6T1YG2s7OjcuXK1KtXj8qVK2NnZ4der+c36YsjhBAil0jqxvGsUnZPIzMQClMZN24cX3b9kmMJUE8p7mi1hhVffw3DhoFOZ9oAc4EcNZFKQkIC/fv3N3UYQgghRLokDSQMCgriwYMH6d4v+QBC6b4hTGH8+PF8+OGHnAVq63RctjOkhOrnn+HAAdMGlwtkexeOKVOmPHVdQkJCNkYihBBCvJzklTjOnDlD3bp107XfsWPHZAChMLkpU6YQGRnJ9L+nU7e1ns2/w7bWPoxN5/P4VZbtLdBTp07l7Nmz3LlzJ9Xl7t272R2OEEII8cKSV+LISDcOmYFQ5AQajYbvvvuO/nX6E3QZ/HrA/KsWREZGmjq0HC/bW6ArVKhAz549adu2bap1sbGxzJ07N7tDEkIIIV5IiRIlsLa2JjY2NkMDCWUAocgpNBoNv8z5haheUawIWMHJyJO0b9+eDRs2YBMcDGPHwi+/QIECpg41R8n2FuhBgwah1+vTXGdhYcGkSZOyOSIhhBDixWi1Wnx8fICMlbJLKmHn5+cnAwiFyWm1Whb9voi2DQ2Nmzt37qR/u3Y8bFgHVq+Ghg3hifk7XnXZnkAPHTqU9u3bp7lOq9VKAi2EECJXyWgljsjISM6fPw9I9w2Rc1hYWLBixQoaN24MwMqgbfyZ77ZhZWAg1KsHV66YMMKcJUdV4RBCCCFym6SBhLdv3yYsLOy528sAQpFTWVtb8+eff1Krdi0SK0DvzjA76Sl6+bIhiX6Bmud5kSTQQgghxEt4shLH8yQfQCgl7EROY29vz8YNG6l4viL6yzCsLXyZVJTj9m1o0ACSPYdfVVmeQJuZmaHVap95sbOzo0yZMgwdOpRLly5ldUhCCCFEpkmeQKenG0dS/2c3NzcKFSqUZXEJ8aIcHR3Zumkr3se84QpMbAYTmzxaef8+NG4Me/aYNEZTy/IE+uOPP6ZixYpotVratm3LmDFjGDNmDG3atEGr1VK5cmWGDx+Or68vAQEBVK1alRMnTmR1WEIIIUSmKFWqFFZWVkD6BhLKDIQiN3BxcWH75u0U/6c4XIUv68Pw1o9WRkRAnz4QH2/KEE0qy8vYFSpUiHv37nH27FlKlSqVYt3Fixfx9/fH19eXr7/+mgsXLlC7dm3ef/99/vrrr6wOTQghhHhp5ubmlClThpMnTz43gY6IiJABhCLXKFy4MDs276COfx2CmgbxUw2IsIKZWy2xW/UH5paWpg7RZLK8Bfrrr79mxIgRqZJnAE9PT0aMGMEXX3wBgJeXF0OHDmX//v1ZHZYQQgiRaZK6cTyvC0fyAYTS/1nkBiVLlmT7pu04b3SGm7CoEhQdEU/vy9+QqE80dXgmk+UJ9M2bNzE3f3pDt7m5OTdu3DBeL1GiBHFxcVkdlhBCCJFpkkrZ3bx5k4iIiKduJzMQityobNmybP1rKw5rHeAWRFrD0lNLGfDnAHR6HaxbB999Z+ows1WWJ9DlypXjp59+SnOa7qCgIH766acUU6FevnwZd3f3rA5LCCGEyDTprcSRlEDLAEKR21SpUoWNqzdiu8oWHs2p8vvJ3/nm87aobt3g7bdh0iR49AtLXpflfaC/+eYbWrVqhaenJx07dsTT0xMw9H9es2YNCQkJzJs3DzBM5T1//nxatWqV1WEJIYQQmebJShw1a9ZMc7ukBFpmIBS5UZ06dVi7bC2tu7Ym/vV4cIPFtzYxzNoGhzhgyhQID4dvvwWzvF0pOcsTaH9/f/bv38+kSZNYtWoVMTExgKFYd9OmTZk8eTJVq1Y1Lrt9+3ZWhySEEEJkKk9PTywsLEhISHjqQMKIiAjOnTsHSPcNkXs1adKEP377g05vdELXW8dJd/DrHcPBZfY4hUTB9OmGKh0//wzP6MKb22XL14MqVaqwdu1aIiMjuX37Nrdv3yYqKoq1a9cak2chhBAit0qqxAFPL2V37Ngx4/+SQIvcrF27diycsxB+A+7BBWdY9uVb4O1t2CAgAHr0gDw8pi1b29ejo6MJDQ0lNDSU6Ojo7Dy1EEIIkaWSunE8LYGWAYQiL+nRowdzv5sLC4Al8PGEuVz49VeoVMmwwcqV0KEDPHxo0jizSrYk0IcOHaJRo0Y4OTlRvnx5ypcvj5OTE40bNzbOyCSEEELkZkkD4q9du0ZUVFSq9Umfd+7u7jKAUOQJAwcO5PtPv4cLEBISQqMePbi2YAEJtWoYNti8GT7+2KQxZpUs75xy8OBB/P39sbS05M0336Rs2bKAYZTykiVLaNCgATt37qRGjRpZHYoQQgiRZZ6sxFG9evUU62UGQpEXjR49msjISD766CNu3bpFvb5tsegWz5bE0pTWOMPkyaYOMUtkeQL9wQcfULhwYfbu3ZuqPN3kyZOpW7cuH3zwAVu2bMnqUIQQQogskzyBPn36dIoEWmYgFHnZBx98QGRkJF999RU3a9yERPBtGcz2rt9TN18+U4eXJbK8C8fBgwcZMmRImrWd3dzcGDx4MAcOHMjqMIQQQogs5eXlZZw47MkZCZMPIJQZCEVeo9FomDZtGsOGDYNVQCQ4X3THt1jdxxspBV98ARcvmizOzJTlCbSZmRmJiU+f6lGn02GWx2sFCiGEyPssLCzwflSF4MmBhMnH+0gLtMiLNBoNM2bMoHfr3jALgpYG0bZtWx4mDSKcOhXefx/q14f//jNtsJkgyzPXOnXqMHPmTK5du5Zq3fXr15k1axZ169ZNY08hhBAid3laJY6k/s8ygFDkZWZmZsybN4/OrToDsH//fjp06MDd8LucuXXCsFFQEDRsCAcPmjDSl5flCfTnn39OeHg4Pj4+9OzZk8mTJzN58mRef/11fHx8CA8P54svvsjqMIQQQogsl1SJ48qVK49b3kg5gFCIvMzc3JzFixfTsmVLALbt24bPZz74eqzkwOiuho1CQ6FJE9ixw4SRvpwsT6CrVKnCgQMHaNmyJWvXrmXKlClMmTKFdevW0bJlSw4cOEClpJqBQgghRC6WfCDh2bNnAQgPDzcOIJT+z+JVYGVlxcqVK2nQoAHkhzCzMADqOK1k30f9QKMx1Idu1QrWrTNprC8qWzoflytXjtWrVxMREcGdO3e4c+cOERERrFq1KsWbjRBCCJGbPVmJA2QGQvFqsrW1Zd26dVQvVh0WAnGgUNTXLmDvF8MN03zHxUGnTrB4sanDzbBML2N3/fr1dG138+bNFNeLFSuW2aEIIYQQ2crLywutVotOpzNW4pAZCMWrysHBgY0bN+Lv78+phafgDVBWioZxP7Fz+rvUHzcdYmPhjTfAywueqJ2ek2V6Al2iRIkXKhCv0+kyOxQhhBAiW1lZWeHl5cXZs2eNLdBJCbSHh4cMIBSvnAIFCrBlyxYaNGjAhcUXoBfoLfU0uvcN2+dMosGIr2DgQMhl3ZsyPYGeN29enpxh6ciRI0yaNInDhw8TFRVFqVKlePPNNxkxYgRardbU4QkhhMghfH19UyTQSSXspPVZvKrc3d3ZunUr1apV496Se9ALdOY6mlybwsZVM2jaZJChX3QukukJdL9+/TL7kCZ35MgR6tSpg5eXFxMmTMDW1paNGzcyevRoLl26xPTp000dohBCiBzC19eXVatWcfnyZYKCgrhw4QIgCbR4tRUrVowZM2bQo0cPWAJmvcxIJJHWB95iTYmitPZqbdgwMhJ+/hnGjoUcPE9Izo0sB5kzZw4Au3fvZuzYsQwZMoQ1a9bQoEED5s+fb9rghBBC5ChJpeyUUixZssS4XBJo8arr3r07Xbp0gUugX6JHi5YEfQKdl3Vmw4UNhv7QHTvCO+9Anz6QkGDqkJ9KEuh0iIiIwNraGkdHxxTLPTw8sLGxMU1QQgghcqTklTh+++034/+SQItXnUajYebMmRQoUAAugO16W8zNzInTxdFpWSc2H1kOd+4YNl60CLp2NSTVOZAk0Ong7+9PREQEQ4YM4cyZM1y7do3Zs2ezatUq3nvvPVOHJ4QQIgfx9vbG7NFPz8ePHwdkAKEQSdzc3Jg5cyYAkYcjqXG1BhZmFsTr4mm3/U3+mjsBkr5srl0LbdpAVJQJI06bJNDpMGjQIEaOHMmCBQvw9fWlRIkSjBw5kh9++IHRo0c/d//g4GACAwNTXC5evJgNkQshhMhu1tbWeHp6plgmrc9CPNa9e3c6d3403fe8/YwrPA5LrSUJ+gQ6bnuT3b9Ogvr1DRtv3w5Nm8KDByaMOLVXLoHW6/XExsam66KUAkCr1VK6dGlatGjBggULWLZsGe3ateOtt95izZo1zz3nrFmzKF++fIpLx44ds/aGCiGEMJknJwmTGQiFeEyj0TBr1ixDVw7g14m/8lur37DSWlG/WH38fJvApk2GmQoBDh4Ef38ICjJd0E945RLo3bt3Y2Njk67LuXPnAJg2bRpffvklS5YsoU+fPnTv3p3Vq1dTr149RowYQWJi4jPPOXz4cE6dOpXikp7EWwghRO70ZAItLdBCpOTm5saMGTMACAkJYdVXq9jRdwfrXl+HrYUt2NrCmjXQvbthh//+g2bNIIfMG5LpZexyOh8fHwICAtK1rYeHB2BoQW7cuDH29vYp1rdv3563336bq1evpvq5LjlXV1dcXV1fPGghhBC5SlIljiSSQAuR2muvvcaKFStYtWoVy5cvp1u3btTuWtu4XllYcPCrUdRycICAAJg6FXLI3BuvXALt7u6e4VrVd+/eTXOmxIRH5VWe1wIthBDi1ZK8BdrDw8PYICOEeCypK8fOnTt58OABw4cPp2HDhri4uKCUYuzmsfxw8Ad+GfozAwcPzlFTfb9yXThehLe3N1u2bOH+/fvGZTqdjuXLl5MvXz5Kly5twuiEEELkNGXKlDFW4pD+z0I83ZNdOUaOHAnA1bCrzDs2D4Xiy/1fEVulQsodH3WzNRVJoNNh4sSJPHjwgJo1a/LVV1/x448/Ur9+fY4cOcKECROwsLAwdYhCCCFyEBsbG5o1awZgrDYghEhbjx496NSpEwDLly/njz/+oKRTSf7u/TdV3Kuwrc82rM2tH+9w6RIMHmyiaA00KqnUhHimzZs388UXXxAYGEhERARlypRhxIgRDBky5IWOFxgYSPny5Tl16lSqvnJCCCFyv4SEBK5fvy6/UgqRDkFBQZQrV44HDx7g4uJCYGAgLi4u6JUeM80T7b1r1xLYpQvlExNNlkdJC3Q6tWjRgp07dxISEkJcXBwnT5584eRZCCFE3mdhYSHJsxDp5O7unqIrx1tvvQWQKnn+9eivfFPwPDza1lQkgRZCCCGEECbXo0cP4zwZy5YtY+XKlSnWLzq5iEHrBvHulneZa3vWBBE+Jgm0EEIIIYQwOY1Gw08//YSzszMAw4YN4969e8b1Fd0qUtC2IADfH/jeFCEaSQIthBBCCCFyBHd3d3788UcgZVUOgApuFdjZbydudm6mCs9IEmghhBBCCJFjvP7660/tyuHr4svOfjuNLdGmIgm0EEIIIYTIMZ7syjF8+PAUXTl8CvqwoOMCU4UHSAIthBBCCCFymORdOYKDg41VOZIUdyxuirCMJIEWQgghhBA5zuuvv06HDh0AWLp0KatWrTJxRI9JAi2EEEIIIXIcjUbD7NmzcXJyAlJX5TAlSaCFEEIIIUSO9GRXjlGjRpk4IgNJoIUQQgghRI7Vs2dP2rdvD8CSJUtYvXq1iSOSBFoIIYQQQuRgT3blGDp0KGFhYSaNSRJoIYQQQgiRo3l4ePDDDz8Ahq4cn3/+uUnjkQRaCCGEEELkeL169TJ25diwYYNJY5EEWgghhBBC5HhPduUwJUmghRBCCCFErpC8K4cpSQIthBBCCCFyjV69euHv72/SGCSBFkIIIYQQuYZGo2HSpEkmjUESaCGEEEIIkau4uLiY9PySQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAZJACyGEEEIIkQGSQAshhBBCCJEBkkALIYQQQgiRAeamDkA8nU6n4+bNm8TExKDT6UwdjhAiE2i1WmxsbChSpAhardbU4QghhHgBkkDnUDqdjosXLhD18CHm8eGYx4WaOiQhRCaIs3Li4cP8xMbE4OnlJUm0EELkQpJA51A3b94k6uFDXC+toMjpn9CYOiAhRKZQwE3fYQSX7sbNmzcpXry4qUMSQgiRQdIHOoeKiYnBPD5ckmch8hgNUOT0T2jjwwl7cJ8HwXdQSpk6LCGEEBkgLdA5lE6nwzwuVJJnIfIgDWARF8r9mET+WRtA4ZJlaNShF/b5nUwdmhBCiHSQFmghhDARjZkZURFhnPhnGzvWLJSWaCGEyCUkgRZCCBPRas0pVMILG7t83Lxyjvt3b5k6JCGEEOkgCbQQQpiYnYMTMQ+jiAi9Z+pQhBBCpIMk0EIIYWIaMw1K6dFLvXchhMgVJIEWQgghhBAiA6QKhzCZzRcTabkoGoDA4Xb4uqScUKLdkmiO3dFx8+182RrXkds6Ptgey/4bOhRQu4iWr5pZU9k99YQXR+/omLwzjr3XE4lNhFJOZgyuZsGomlbZGrMQQgghso8k0MJkTtw1/FxtpoH15xNTJdAngnRUdMveWdqO3tFRL+AhRR3MmNTQCr2CWYfjaTj/If++aUeZgo/j+ftSIu2WRFPFXctHDaywt9RwKVTPzQippCCEEELkZZJAC5M5eVePgxXULWrOuvOJjK/7uNU2NEZxI0LRq0L29jL6aEccNubwz0BbCtgazv1GRQu8Z0Tx/vY4Vna3BSAiTtFndQxtvMz5o7sNZhqp2C2EEEK8KqQPtDCZE3d1VHDV0tbbnH9u6LgfrU+xDsj2Fug91xJpWsrcmDwDeOQzo2Fxc9afTyQq3tC6vPi/BO4+VHzW2AozjYaH8Qq91PAVQgghXgmSQAuTiNcpzt3TU8nNjLbe5ugUbLiQaFx/8lECXck9fU/RBJ3iXrQ+XZdnJbpxOrCxSN2abGsB8To4FWyIa+vlRBys4FakosyMKOy/iMThi0iGrY8hNlESaSGEECIvky4cwiROh+hJ0BtamIvlN6OCqxnrzifSu5IlACeC9FhpoUyB9CXQ+27oaLQgOl3bXhltTwnHtLtclClgxoGbOnR6hdbMsE28TnHwliFxvvWof/OFB3oS9dBhaTQDq1jyRRMtO6/q+PHfeMLiFEu62KYrFiGEEELkPpJA50JjNsVyPChn1Iut7K7l+5bWGd7vpLGLhiFBbuttzsxD8SToFBZaDSfu6ijnamZMYp+nkpuWLb3Tl7S62z/9mMOrWzLsr1gGro1lfF1L9Aqm7o7jTqQhcY551LocFa+IToCh1Sz4oZXh9ncua0G8TjHnSAJT/HV4Fcje7idCCCGEyB6SQOdCx4N07LqWMxLoF3UiSI8GqPCoj3Nbb3O+2BvP7ms6/EtoCQzR06O8RbqP52SjoWmpl386D/Wz5Ea4nq/3x7PgRAIAfoXMGF/Xks/2xGNvaUi+bcw1gOL1Cilj7FnBgjlHEvjnpiTQQgghRF4lCXQulFY9YlN50VhOBuso5aQxJqS1imgpaKth3flECuXTEJsIldye3n0jeRcLMHSzeBCTvr7HLraaZ7Zsf9bEmnfqWBEYoiO/lYYKblre3xYLgPejLiWF8mkIDAE3u5THcX10PTSdsQghhBAi95EEOhd6kS4TOc3Ju3rqFn2cfJtpNLTyNGfd+QRqFTEsT16BY/qBOHZe02FuBhsvJLKimw2tvB63/u7PpD7QSZxsNNQr9vjlsfVyIkUcNPgUNCTQ1Ty0bLmsMwwiLPh4v9uPunq42Mn4XCGEECKvkgRaZLugKD3BDxWVnihR19bbnN9PJrDklKHrRPIW6JN39ey7ruOP7jYs62qDTp9i10zrA52WZacSOHRbzzfNrIz1nruXs2Davnh+PRZP45KPX0Zzj8Zjbgb+JXLOrwRCCCGEyFySQItsdyLIkP1WfKKLRovS5liYwbpzhm4cyWsxnwzWMamhFQ2KG56yZk/kp5nVB3r3tUSm7IqjeWlzCthoOHBTR8DxBFp6ahldy9K4XRUPLQMqWzDveAKJ+mgaFjdn59VEVpxO5L16lhTKJy3QQgghRF4lCbTIdiefMklKfmsN9Ypp2XFVlyK51ivF6RA9ncpm/dO1cD4ztGbw9f54IuMUJZ3MmNrYirdrW2L+RL/p2W2tKZbfjIDj8aw+k0hxRzO+a2HFmFpWTzm6EEIIIfICSaBFtnu3rhXv1k07ydze1y7VsosP9OSz1GRLq25pZzM2v5E6hrRYaDVM8rdikr8kzEIIIcSrRH5nFjneybv6dM9IKIQQQgiR1SQrETnef3d1VHSVQXlCCCGEyBmkC4fI8T5plPvL9gkhhBAi75AW6HTasmUL9erVw9bWFicnJ7p27crVq1dNHZYQQgghhMhmkkCnw/r162nZsiVxcXFMmzaNcePGsWvXLurVq0dISIipwxNCCCGEENlIunCkw4QJEyhVqhT79u3D0tJQC7hdu3ZUrVqVadOm8b///c/EEQohhBBCiOwiLdDP8eDBA06fPk2nTp2MyTNApUqVKFu2LEuXLjVhdEIIIYQQIrtJAv0ccXFxANjY2KRaZ2try+3btwkKCsrusIQQQgghhIlIF47ncHNzw9HRkX379qVYfv/+fU6fPg3ArVu3cHd3f+oxgoODU/WVvnjxYuYHK4QQQgghstwrl0Dr9Xri4+PTta2VlRVmZmYMGTKEL7/8kvfee48BAwYQERHB+PHjjceJiYl55nFmzZrFJ5988tKxCyGEEEII03vlunDs3r0bGxubdF3OnTsHwJQpUxg4cCBfffUV3t7e+Pn5YW5uzsCBAwGwt7d/5jmHDx/OqVOnUlzWrFmT1TdVCCGEEEJkgVeuBdrHx4eAgIB0bevh4QGApaUlc+fO5bPPPuP8+fO4ubnh7e1Nz549MTMzw9PT85nHcXV1xdXV9aVjF0IIIYQQpvfKJdDu7u7069fvhfZ1c3PDzc0NAJ1Ox86dO6lZs+ZzW6CFEEIIIUTe8cp14cgs33zzDXfu3GHcuHGmDkUIkUP8cyMRs08imLo7ztShCCGEyEKvXAv0i1i4cCErV66kQYMG2Nvbs3XrVpYvX86bb75Jly5dTB2eECIH0CvF2M1xVC8s7RJCCJHXyTt9Onh7e/PgwQM+/fRTRo0axblz55g9ezY///yzqUMTr4j5x+PRfBLB1TC9qUPJsK/2xeEzIwq9UqYO5aXMPhxPse8iiUtM+3b8fCSBmoW1lC2ozebIhBBCZDdpgU6HGjVqsGvXLlOHkedsvphIy0XRAAQOt8PXJWXi0W5JNMfu6Lj5dr5si6nfmhgWnEh46vqbY+0p7GDGoVs6FpyIZ8dVHVfD9BSw0VCriJapja3wLvBiCdT+G4n8fSmRMbWscLTWvOhNyFEi4hRf7ovnm2ZWmGkMtym99zGQ4ftZrxRu30Txbh1LmpT8f3t3Hh5Vdf8P/H1nzUxWlhASdkzCEiAQ9qXsCiooKgStKFSUL1patFJotVbhKyqiUvqzKD4WsFZRQK1Fga8g+x6WEPZE9iUbgZB1lnvv+f0xZJIhk2Qmmezv1/P4aO7cO3Pmnpz4vmc+c64Or261YO8VBQLAwNZavHuvH3q2dN8/R1IVvLHdit2XZVhkoGMTDWb01uP3/Y2OdvfU443tViw/bHNuK5JVoOJv+23Y/6w/XtxkqcypIiKieoQBmmrNsXQFAKCRgB+S5VIB+liagh5hNTub9z+99Rjd0fU1hQBm/mhB+xCNM9gt2mPFnisKJnXVo0eYBml5Ah8etCFueT72P+uPbi28b/feKwrm77BhWk9DqQD9VA89Hu+mh7GeTW6uOGqDrAo80V3v3ObpOQa8P88Hrym4USAwuI0WQ1bmo02QBq8PM0IVwLJDNgxblY+Dz/qj012zxD+dkzF+dQF6tdTitaFGBBgknLul4mpO8Wyzn07C1Fg9Pthnw+/6GSBJxX306lYrXhxQut+IiKhhYoCmWpOUriLICAxuo8P6ZBlzBxfP6t0qFLiSI/Bk95qtMhrYRoeBbVy37b4so8AOPFkiBP5hoAFfPqaFQVscmCbH6ND9o3y8s9uGfz9a+tbvVaHVSNDWw4KrlYl2PNRJDz9d8Xny9BwD3p/nDSky2gVLeGu3DSYdsG+6Gc3MjhM3pYce0R/m4ZWtVnwTb3Yek2MVePq7QjwYpcO6eJNzptyd+Bg93t3rmBEf2cHx5/NoqoKE6wr+8YCfdyeHiIjqLQZoqjXH0hV0b6HFuGgdfr/RgqwC1Rl2imana3oG2p0vj9shAfh1iXA3qE3poRPVTIuYFhqcvqF4/RpvbLdg/g7HnS07LM1zbr8wOwDtQzRYlWjDb763OH8u2v/sLH8s2GHD+mQ7DFoJM3vrsWCEEVdzBGZttGDbBRlmvYQ/DjLg5UHGUq97LUfFa9us+DFFRrZFILKpBi8PNOCZXgbnPrlWgde2WfGfM3ak5gkEGyXEttRg0Wg/xIWX3T8XbqlISlfxhwGGMvcp4u4cA96f5x9TZDwYpcPnSXaMjdQ5f58AIDxQg2HtdPghWUaeTSDAIDlfOz1fYOFIR5lJvk3ApIfbIN07QoumJgnfn5GdAXrHJRlnb6ho9YGj325bBXQa4NwtFSsf9u2FFBER1Q0M0FQrbIrA2RsqnovTY1y0Dr/d4Jg9fCrWEbaS7gTo2JaeTbvaFYHbVs++pNbUJJU7y3j38645KWNQGy3ah5TfFiEE0vMEYlp4P1X8aBc9krNUrD4hY8kYI5qbHe0LNZffzsnrCtGluQbvjPLDjyky3txlQ1OThOWH7RjZQYtFo/3wxXE75my2om8rLYa2Kx7y6XkqBvwzHxKAWX0NCPWXsPEXGdP/a0GOVeDFAY7APfPHQqw7JWNWXwO6hmqQVSiw+7KC05lquQF67xVHH5a3D+DdOQbKPs9peSqOpqpYMNyIT4/aYdKXPndmPWBTgBMZCga0dpyLLedlBBmBa7kCE77OQ3KWCn+9o2xmyVg/l9lzx/vRYM8V2fnzjN4GPN6tOPjP3mRBhxAN/jSk9AULERE1DAzQVCtOZaqwq44Z5rbBGnRvocH65OIAfSxNhVELdGrmWRjdc0XBiM8KPNrXMYvrWYD+v3MysgpFqdICd744bse1XIEFIyre9249wrSIC9di9QkZEzrrPQqSANAvQovl4x2znDN669F+aR5e/smKt0cZMe9OgHuiux4R7+dixVG7S4B+dasVigocf97fOVM7s48BT3xTgDe2W/E/vQ0w6SX8mCzjuTg93h9TXKIwd3DFbTtzZ4a4Q5Py34s35xgo+zxvSJHhpwNGdtChUzMN9l9VoKgCWo2jr22KwIFrjjZdK1HbnHJThawCD39VgOm9DHh7lBbbLyr4fwdtyLYKrH7M7PI6HUM0+Pxy8ZcgzXoJ5hJh3aSTEGCQWA9NRNSAMUBTrUhylmg4wtW4aB3+kWCDXRHQayUcS1cQ00LjDD8ViQ3TYvNT5op3BNAywPNg8+VxO/QaID6m/KFy5oaC326wYGBrLabGeh+gK+vZuOLyCK1GQp8ILa7myJgeV9yGED8JnZprcP5W8RJ4Qgh8c9qO+K56CAA3CoofG3OPDl+dkHEkVcHgtjqE+Ek4cE3B9VwVEYGez65nFTpKGYpKJcri6TkGyj/PG1JkjOigg0kv4YW+Bjz/owXT/2vB3MEGqAJ4c6cVqbmO4FxYYim6PJtAgR2Y2VuPv9/vuEh4tIseNkVg+WE7FgxXEFVixY8mJgmFMlBgFy7BuciqCSzbICJq6Big66lVsGEVyl4KrCLToMc0FIevks+3Hf4u+w5HvkfP4Y1jaSokAN3v1DiPi9bh7d027LykYHh7LU5mqi4fi1ekiUnC6I6+/XXOswl8f1bGmLtqae+WlqfiwS8LEGyUsC7e5HHo94W2wa6vFWyU4KcDmt/V3mCjhKzC4tCYWSCQbQE+OWLHJ0fc/x5l5Dv2f/deP0z9TyHaLMlD73ANHojS4elYAzpWMLPsCU/PMVD+ebYrApvPy3h7lCMAz+xjwJXbKhbvtTmXzOsTocHcwQYs3GVzCfUmnQTAdaUQwFGPvfywHfuuugboouWsOb9MRNR4MUDXUxehYofk/ZfVigwXrnWpLs93VylxWa9z93N4IylDQccmkjPIDGitRXOzhPXJMiICJVhkIDas7EBV8qN5wPHx/M1Cz2qgQ82SRyH3P2fsbleGKOm2ReD+LwqQbQF2/cbs1QytL7hbmUNbxlsreR8T9c5/T+mhL3PGvOjTgfgYPX7VVovvzjjWqV6814ZFe2z4Nt6E+6PKPjfNTBJk1fElxECj+0Z5co6Bis/z7ssKcqzAA1HFf9IWjvLDnEFGnMxUEGyU0D1Mi1d+dqzRHF2iNCgiUMLJTCDM37WNLe78fOuu36tbFgGzHm5rrImIqHFggK6n2kODYVUIsO3vugllec9X1va7n8MbSekqBrcpfl6NJOH+SB3WJ9sxoLVje8kVOJbut2L7JQU6DbAxRcbaSa7hbW811EB/cdyOAAPwUCf3w8QiC4xfXYDkLBVbnjKXWsfaWzUZx0LNEgINjgsRT2buwwM1eKGvAS/0NSAjX0Xc8nws3GUrN0B3vrPW8oVstczVVCo6x4Bn5/nHFBldQzWlasebmCQMaVv83FvOy2gdJKFz8+L9eodrsfm8gmu5Ap2aFx97/U65R6i/63NeyFbRpXnNXigREVHdwgBdT02DodLlE94+390lHVWVlqciI18g9q5QNS7asfzY6hOOj9xLzkAnpavYc1nBungTvp5ognLXHa19XQOdma9iy3kFT3TTu61zVVSByesKse+qgu8fN2Ggm+XWvOV/ZzY+21L9t7zWaiQ81lWPL4/bcSJDKXVDksx8FaH+GiiqQJ4NCC7xhbgW/hpEBEqwKuW3c+CdC6RD193fEKeicwx4fp43pMgYF11+H3x9wo6E66rLXREBxwz7O3ts+OdRm3NpOgD49IgNOg0wvL1r24+kqh5/4ZGIiBomBmiqccfSHOm3x10lGmPu0UGvAdafdZRxlKyJTcpQ8Powo3MVCc1deczXNdBfn5Qhq2WXFrz8kxX/PStjfLQONwsF/p1kc3l8Sg/XixFpfg6GtdNi+7SyL0Z631nu7dWtVjweo4NeK2F8tM4ZrH3tnVFGbLsgo/+n+XguzrFE3c1CgSOpCracl3FzXhBybUDrD3IxsasesWEaBBgkbLkgI+G6ivfvK3+Zto5NNOjWQoMt52WXdaWLVHSOAc/O84VbKk7fUPHRg8X9v/OSjAU7rLjvHh2amSTsv6pgZaIdYyO1mH3XutS9wrV4pqceKxLtkNUCDGunw/aLMtaekvHnIQaXcpHD1xXcLBR4uJwZcyIiavj4fwGqcUll3CQl2E/CkLZabLuouIRrVQicylTxSJea+3X94rgdLfylUrecLpKY5ngP65NlrE+WSz1eMkDn2RwzteGB5Qfhvq20+N8RRnx8yIZNv8hQhaPcpLoCdFiABgef88eCHVZ8e9qOZQkCzcwSYkIdN0kBHOsmv9DXgJ/Oyfj2tB2qACKbarDsAT8837fiT0Ce6anHX7dbUWgXpWqGKzrHgGfneUOKjGAjMLht8fO0CtRAqwEW77Uh1yrQoYkGb4404g8DDdC5qX//eJwf2gZrsDLRhu9Oy2gXosGSMUbnWthF1p6yo22whJEdav8GP0REVHskIUT1f15MpZw8eRLdunXDiRMnEBMT4/ZxZJ5BzPZnaqF1dUtyloKhKwuQNiewtptSKRtS7Bj3ZSGOzfR3rjrSWNy2CHT8ex7eHW3E9DjflRyV9MAXBQgwAGsmeVbCU1lWWaD90jz8abABswdU/SYpJ4evwG1tM1xM3Inc2zeRl30TD0+bjchuvX3QWiKihq2iHFXd+E0YqvOS0lWP70hYF227oODxbrpGF54Bx6cKcwcZsHivDWo1XasPb6/FSx7cLryqViY61que2af6X4uIiOo2lnBQnXc8XUGPFvU3fC6+z6/inRqweUOK74pYHeYOrplbZs/sY2B4JiIiAAzQVA/MH9G4AygRERHVLfX3c3EiIiIiolrAAE1ERERE5AUGaCIiIiIiLzBAExERERF5gQGaiIiIiMgLDNBERERERF5ggCYiIiIi8gIDNBERERGRFxigiYiIiIi8wABNREREROQFBmgiIiIiIi8wQBMREREReYEBmojIR/ZdkaGZn4M3d1pruylERFSNGKCJiHxAFQIv/Z8VfVvxzyoRUUPHv/RE9cCqRBuk+Tm4mK3WdlO89u4eKzp/mAdViNpuSpV8fMiGtktyYZXdv49PDtvRv5UWXZpra7hlRERU0xigqdb83y8ypPk5kObn4FSmUurx8asL0PqD3FpoGXAkVcFDqwvQdFEOzAtz0G1ZHv5+wPVjeassMG+zBRHv58K0MAf9P83D5nNypV9z7xUZb2y3INtSv4NmSTlWgUV7bJg32ACNJLndZ+FOK6T5jnPsTkqWgsfXOX4XzAtz0PnDPCzYYUWBvfR5UoVA6OJcvLvHisPXFYz9dz6C3s5B4Ns5uO/zfCSmlf49K1JRn0/rqYdNAZYftpU6NqtAxd/22zB/hLGiU0JERA0AAzTVmmPpjjCjkYAfkksHz2NpCnqE1fxs3k/nZAz8Zz4y8gVeG2rE0rF+GBetw9Uc18A27ftCfLDfhie767F0rB+0koQHvizA7suVC9F7ryiYv8PmNkA/1UOPwlcD0S7YfQitq1YctUFWBZ7ornf7+NUcFW/ttsLf/cO4cltFv0/zsf+qgln9DPjbWD8MbKPF69uteOKbwlL7H7ym4EaBwOA2WgxZmY/ztwReH2bEX4cakXJTxbBV+Th7o3SI9qTP/XQSpsbq8cE+G8Rds+mvbrXixQEGhPjVr/4hIqLK0dV2A6jxSkpXEWQEBrfRYX2yjLmDi2fvbhUKXMkReLJ7zV7j5VgFnv6uEA9G6bAu3lTmrOnBawq+OiFj8b1GzBnkaPfTsXp0W5aHuZut2Dvdt0NLq5GgrYeXuysT7Xiokx5+Ovfncc5PFgxorYWiAjcKSl84fJ5kR7YF2P0bM2JaOC6mZvQGVAH865gdtwoFmpiKn3tDiox2wRLe2m2DSQfsm25GM7PjxE3poUf0h3l4ZasV38Sbncd42ucAEB+jx7t7bdh2UcHIDo4+PpqqIOG6gn884Of9CSIionqpHv4vmRqKY+kKurfQYly0DvuuKMgqUF0eA1DjM9BfHrcjPV9g4UgjNJKEfJtwW7u77pQdWgmY0dvg3OankzC9lwH7riq4ctu7WuU3tlvwx82OcoEOS/OcpS1FNc9310C/sd0CaX4OkrMUTPm2EMHv5CB0cS5e22qBEAJXbqt4+KsCBL2dg5bv5eL9ve5XhbiWo+KZ7wsR9l4ujG/mIGZZHlYcdS1RyLUKvLjJgvZ/c+zTYnEu7v08H0dSyy6HAIALt1QkpasY3cF9H+68JGPdKRl/G1N28MyxOs59WIBrqA0PkKCRAMNdT/1jiowHo3TYdUnG6I46Z3gGgPBADYa10+GHZBl5tuI+9bTPAaB3hBZNTRK+P1P8KcOOSzLO3lDR6oM8tHwvF1+ftGPRHit+833pGXIiImoYOANNtcKmCJy9oeK5OD3GRevw2w2O2cOnYh2BNOlOgI5t6dk1nl0RuG31rHa4qUkqc5Zxy3kZQUbgWq7AhK/zkJylwl/vKKFYMtbPOZN6NE1BdDMNgoyuz9OvlSPRJaYpaBPs+fXpo130SM5SsfqEjCVjjGhudjxvqLn8koDJ6wrRpbkG74zyw48pMt7cZUNTk4Tlh+0Y2UGLRaP98MVxO+ZstqJvKy2Gtise8ul5Kgb8Mx8SgFl9DQj1l7DxFxnT/2tBjlXgxQGOmfWZPxZi3SkZs/oa0DVUg6xCgd2XFZzOVBEXXvYFzt4rjj50t4+iCvxuowXPxunRvZyLpOHtdVi0x4bp/7Vg/nAjmpkk7L2i4KNDNvy+nwH+huLzk5an4miqigXDjfj0qB0mfelzZ9YDNgU4kaFgQGvHufC0z4vEhWuw50pxgJ7R24DHuxXXoMzeZEGHEA3+NIT10EREDRUDNNWKU5kq7KpjhrltsAbdW2iwPrk4QB9LU2HUAp2aeRZC91xRMOKzAo/2vTA7AO1D3AfTlJsqZBV4+KsCTO9lwNujtNh+UcH/O2hDtlVg9WOOj/5TcwXCA0s/R9G267nefRGwR5gWceFarD4hY0JnPdqHePa++0VosXy8CQAwo7ce7Zfm4eWfrHh7lBHz7gS4J7rrEfF+LlYctbsE6Fe3WqGowPHn/Z0ztTP7GPDENwV4Y7sV/9PbAJNewo/JMp6L0+P9EjPFcwdX3LYzd2qNOzQp/V4+PmTHpWwVW54yl3qspLGROvzvCCPe2mXFf88Wh9ZXf2XAmyNdZ643pMjw0wEjO+jQqZkG+68qUFQBrcbRJzZF4MA1R5uulaht9rTPi3QM0eDzy3bnz2a9BHOJsG7SSQgwSKyHJiJqwBigqVYkOUs0HOFqXLQO/0iwwa4I6LUSjqUriGmhcYafisSGabG5gjBWpGVA2c+ZZxMosAMze+vx9/sdAe3RLnrYFIHlh+1YMFxBVDMtCmUBo5ui5KLZysIyljrztWfjiktItBoJfSK0uJojY3pc8YxoiJ+ETs01OH+ruKxECIFvTtsR31UPAeBGifKZMffo8NUJGUdSFQxuq0OIn4QD1xRcz1UREej5rHpWoYBOAwQYXM93VoGKv2634rWhRoT6V/x87UMkDG2nxWNd9GhmdgT6t3bZ0DJAg1n9it//hhQZIzroYNJLeKGvAc//aMH0/1owd7ABqgDe3GlF6p0Lm5L942mfF2liklAoAwV24RKci6yaYPL4HBERUf3EAF1fJdqARHv5+7TUAmNLzNKlKcAmS8XPPc3f9edV+e7366kHehrcP1aBY2kqJMD58f24aB3e3m3DzksKhrfX4mSm6vKxeEWamCSM7lj1X2eTTgJQetWIX3fXY/lhO/ZddYQpk06CVSkdki13gpmpjC/N+Vrbu1blCDZK8NMBzc2aUtuzCovbm1kgkG0BPjlixydH3P8eZeQ79n/3Xj9M/U8h2izJQ+9wDR6I0uHpWAM6uplZ9sRftlrR1CThd/0r/t356oQdM9ZbkPy7ALQOcrzeo130UAHM22LBE90cdc52RWDzeRlvj3L8vs/sY8CV2yoW77Xhs2OO99cnQoO5gw1YuMvmEuo97fMiReXRnF8mImq8GKDrq2wVuFT+l7hKsQjvjwHKPqZ95b/gl5ShoGMTyRlkBrTWorlZwvpkGRGBEiwyEBtWdkAr+dE84Ph4/mahZ7O+oWapzJntiEAJJzOBMH/Xx1vc+fnWndcID5RcygCKFM1wRrgp76gO7lbm0Jbx0iW/F6fe+e8pPfSYGuv+QqXo04H4GD1+1VaL787I+OmcjMV7bVi0x4Zv4024P6rsi5xmJgmy6vgSYuCdWvGULAWfHLHjb2P87pS5OBpikQG7ClzMVhFklND0zsoayxJs6BWudYbnIg9F67Aq0Y6jaSpGd9Rg92UFOVbggajiP2kLR/lhziAjTmYqCDZK6B6mxSs/Oy4go0uUBnna50VuWQTMeritsSYiosaBAbq+CtEA7SoIsC3vetxPqvgYd8o6xsM6XXeS0lUMblP8vBpJwv2ROqxPtmNAa8f2kitwLN1vxfZLCnQaYGOKjLWTXMPbXh/VQPcO12LzeQXXcgU6NS/eXlTTXFRy0DNMi20XbMixCpcvEhbV2Pa8+9x7oCbjWKhZQqDBcSHiycx9eKAGL/Q14IW+BmTkq4hbno+Fu2zlBujOd+7IdyFbdfbltVwBVQC/32TB7zeVPqbD0jzM7u9Y7xkA0vMFmripJbbfqTiR7/z7xxQZXUM1pWrHm5gkDGlb/P62nJfROkhC5+bF+3na50UuZKvo0pwLGBERNWYM0PVVT4P35RMttaXLMzxRmWPKkZanIiNfIPau1RfGRevweZIdq084PnIvOQOdlK5iz2UF6+JN+HqiCcpdq8T5qgY6PkaPd/bY8M+jNuc6vwDw6REbdBpg+J1Z94lddXhvnw2fHLY514G2ygIrEx23c/ZmBY4iRStK1MSdCLUaCY911ePL43acyFDQrYVrX2Tmqwj110BRBfJsQHCJENvCX4OIQPclLCUNvHOBdOh68Q1xurXQ4LvJpWuE/7LVilybwNKxfrinRGlIdDMNfjonIzlLQXSJMorVJ+zQSMWz5BtSZIyLLv/P2dcn7Ei4ruK9e40uq7B42udFjqSqeLKMG8MQEVHjwABNNe5YmiP99rirRGPMPTroNcD6s44yjpJr+CZlKHh9mNG5ioTmrgleX9VA9wrX4pmeeqxItENWCzCsnQ7bL8pYe0rGn4cYnF+i699ah0lddfjzz1Zk5AtENtXgs2N2XMxW8c+HSgd5aX4OhrXTYns5FyO97yz39upWKx6P0UGvlTA+WueyVJsvvTPKiG0XZPT/NB/PxTmWqLtZKHAkVcGW8zJuzgtCrg1o/UEuJnbVIzZMgwCDhC0XZCRcV/H+feUv09axiQbdWmiw5byMZ3o5LvaamzWY0Ln0xcXf9jvWnp7Q2TWY/nGQARtTZPxqZQFm9TWgmVnCD8kyNv4i49leekQEanDhlorTN1R89GBx/++8JGPBDivuu0eHZiYJ+68qWJlox9hILWYPcL3w9LTPAeDwdQU3CwUe7sQ/nUREjRn/L0A1LqmMm6QE+0kY0laLbRcVl3CtCoFTmSoe6VIzv64fj/ND22ANViba8N1pGe1CNFgyxuhcF7nIvx4x4bWtVnye5LgjXo8wDX54wuyyVBwA50073C17V1LfVlr87wgjPj5kw6ZfZKjCUW5SXQE6LECDg8/5Y8EOK749bceyBIFmZgkxoRosGu0ooTDrgRf6GvDTORnfnrZDFUBkUw2WPeCH5/tW/AnIMz31+Ot2KwrtolI1w0Pb6bB3uj/e2G7FskM2ZBUIdGiiwcKRRswd7Hj9DSkygo3A4LbFv0+tAjXQaoDFe23ItTqOeXOkEX8YaIDOTf27p32+9pQdbYMljCzj5jBERNQ4SEKUccstqlYnT55Et27dcOLECcTExLh9HJlnELP9mVpoXd2SnKVg6MoCpM0JrO2mVMqGFDvGfVmIYzP9y71pSEN02yLQ8e95eHe0EdPjKrdiS0Ue+KIAAQZgzSTPSngqyyoLtF+ahz8NNmD2gKrfJOXk8BW4rW2Gi4k7kXv7JvKyb+LhabMR2a23D1pLRNSwVZSjqhu/CUN1XlK66vEdCeuibRcUPN5N1+jCM+D4VGHuIAMW77WVeXvsqhreXouXBlRPOC9pZaIdeo1jiTwiImrcWMJBdd7xdAU9WtTf8Ln4Pr+Kd2rA5g0pvitidZg7uGZumT2zj4HhmYiIADBAUz0wf0TjDqBERERUt9Tfz8WJiIiIiGoBAzQRERERkRcYoImIiIiIvMAATURERETkBQZoIiIiIiIvMEDXUVqtFrKxCXiXG6KGRwCwG5pAle213RQiIqoEBug6ymQyQTYE42rX5xmiiRoQAeBq1+ehGINhLcit7eYQEVElcB3oOqp169awFBYi455JyGpzH/TWW7XdJCLyAbuhCRRjMApuZyHj4unabg4REVUCA3QdpdVqERkVhatXryIrU8GNAju0Oj00Gn5oQFSfqVY7rLcuIePiaQhVcWwUApIkQeL4JiKqFxig6zCtVot27dpBLbyNnevWQKc3oGmL8NpuFhH5WEFeLowmM/wDQ2q7KURE5AEG6HqgVYdoNAtrhfOnEwEAJv9ASJJUu40ioioTQkV+zm3k3rqBLr0Ho0WrdrXdJCIi8kCjDtCpqalYunQpDhw4gEOHDiEvLw/btm3D8OHD3e6/d+9ezJ07F0eOHEFQUBDi4+Px1ltvISAgoFrbaTD6YdSjTwPfAlnp13Aj9QqEUKv1NYmo+kmSBib/AETH9sfICU+xRIuIqJ5o1AH67NmzWLRoEaKiotC9e3fs27evzH0TExMxatQodOnSBR988AGuXr2K9957DykpKdi4cWO1tzWiXSQmv/AKrl1IRn5ONoTg2hxE9Z2k0SCoSTNEtIuCTq+v7eYQEZGHGnWA7t27N7KystC0aVOsW7cOkyZNKnPfV155BU2aNMH27dsRFBQEAGjfvj2ee+45/PTTT7jvvvuqvb0Gox86dO5R7a9DRERERGVr1J8XBgYGomnTphXul5OTg82bN2PKlCnO8AwATz/9NAICArBmzZrqbCYRERER1SGNOkB76vjx45BlGX369HHZbjAY0LNnTxw9erSWWkZERERENa1Rl3B4KjU1FQAQHl56Cbnw8HDs2rWr3OMzMjKQmZnpsu3UqVMAgF9++cVHrSQiIiJqHIryk9VqrZXXbzABWlVV2Gw2j/Y1Go1eLQNXWFjoPO5ufn5+zsfLsmzZMsyfP9/tYxMmTPC4HURERERU7Pjx44iLi6vx120wAXrnzp0YMWKER/uePn0anTt39vi5TSYTAPdXORaLxfl4WV544YVSX1BMTEzElClTsGbNGnTt2tXjtpDv/fLLL5gwYQL+85//IDIysrab06ixL+oO9kXdwb6oO9gXdcepU6cQHx+P6OjoWnn9BhOgO3fujJUrV3q0r7tSDE/2LyrlKCk1NRURERHlHt+iRQu0aNHC7WNdu3ZFTEyMV+2h6hEZGcm+qCPYF3UH+6LuYF/UHeyLuqPk4g41qcEE6JYtW2LatGnV8tzdunWDTqfDoUOHEB8f79xus9mQmJjoso2IiIiIGjauwuGB4OBgjB49Gv/+97+Rm5vr3P75558jLy+v3PWjiYiIiKhhaTAz0JX15ptvAgBOnjwJwBGKd+/eDQD4y1/+4txv4cKFGDRoEIYNG4YZM2bg6tWreP/993Hfffdh7NixNd9wIiIiIqoVjT5Av/baay4/r1ixwvnfJQN0XFwctmzZgnnz5uGll15CYGAgpk+fjrfffrtSrxsaGorXX38doaGhlWs4+Qz7ou5gX9Qd7Iu6g31Rd7Av6o7a7gtJCCFq5ZWJiIiIiOoh1kATEREREXmBAZqIiIiIyAsM0EREREREXmCAJiIiIiLyAgO0j1mtVsybNw8REREwmUzo378/Nm/e7NGx165dQ3x8PEJCQhAUFISHH34Y58+fr+YWN1yV7Ys33ngDkiSV+sfPz68GWt3w5OXl4fXXX8fYsWPRtGlTSJKEVatWeXx8dnY2ZsyYgdDQUPj7+2PEiBE4cuRI9TW4AatKX6xatcrtuJAkCWlpadXb8AYoISEBs2bNQkxMDPz9/dG2bVvEx8cjOTnZo+M5LnynKn3BceFbJ0+exKRJk9CxY0eYzWY0b94cQ4cOxfr16z06vibHRaNfxs7Xpk2bhnXr1uHFF19EVFQUVq1ahQceeADbtm3DkCFDyjwuLy8PI0aMwO3bt/HKK69Ar9djyZIlGDZsGBITE9GsWbMafBcNQ2X7oshHH32EgIAA589arbY6m9tg3bhxAwsWLEDbtm0RGxuL7du3e3ysqqp48MEHcezYMfzxj39E8+bNsWzZMgwfPhyHDx9GVFRU9TW8AapKXxRZsGABOnTo4LItJCTENw1sRBYtWoQ9e/Zg0qRJ6NGjB9LS0vDhhx8iLi4O+/fvR7du3co8luPCt6rSF0U4Lnzj0qVLyM3NxdSpUxEREYGCggJ88803eOihh7B8+XLMmDGjzGNrfFwI8pkDBw4IAGLx4sXObYWFheKee+4RAwcOLPfYRYsWCQDi4MGDzm2nT58WWq1W/PnPf662NjdUVemL119/XQAQmZmZ1d3MRsFisYjU1FQhhBAJCQkCgFi5cqVHx3799dcCgFi7dq1zW0ZGhggJCRFPPPFEdTS3QatKX6xcuVIAEAkJCdXYwsZjz549wmq1umxLTk4WRqNRPPnkk+Uey3HhW1XpC46L6ifLsoiNjRWdOnUqd7+aHhcs4fChdevWQavVulwh+fn5Yfr06di3bx+uXLlS7rF9+/ZF3759nds6d+6MUaNGYc2aNdXa7oaoKn1RRAiBnJwcCC6VXiVGoxEtW7as1LHr1q1DWFgYHn30Uee20NBQxMfH4/vvv4fVavVVMxuFqvRFSbm5uVAUxQctarwGDRoEg8Hgsi0qKgoxMTE4ffp0ucdyXPhWVfqiJI6L6qHVatGmTRtkZ2eXu19NjwsGaB86evQooqOjERQU5LK9X79+AIDExES3x6mqiqSkJPTp06fUY/369cO5c+eQm5vr8/Y2ZJXti5I6duyI4OBgBAYGYsqUKUhPT6+OplI5jh49iri4OGg0rn+q+vXrh4KCAo/rRcl3RowYgaCgIJjNZjz00ENISUmp7SY1GEIIpKeno3nz5uXux3FR/TztiyIcF76Vn5+PGzdu4Ny5c1iyZAk2btyIUaNGlXtMTY8LBmgfSk1NRXh4eKntRduuX7/u9ribN2/CarVW6lhyr7J9AQBNmjTBrFmzsHz5cqxbtw7PPvssvv76a/zqV79CTk5OtbWZSqtKP5Jvmc1mTJs2Df/4xz/w3XffYe7cufj5558xaNAgjz7RoYp98cUXuHbtGiZPnlzufhwX1c/TvuC4qB4vv/wyQkNDERkZiTlz5uCRRx7Bhx9+WO4xNT0u+CVCHyosLITRaCy1vWj1hsLCwjKPA1CpY8m9yvYFAMyePdvl58ceewz9+vXDk08+iWXLluFPf/qTbxtLZapKP5JvxcfHIz4+3vnzhAkTMGbMGAwdOhQLFy7Exx9/XIutq//OnDmD3/72txg4cCCmTp1a7r4cF9XLm77guKgeL774IiZOnIjr169jzZo1UBQFNput3GNqelxwBtqHTCaT2xobi8XifLys4wBU6lhyr7J9UZZf//rXaNmyJbZs2eKT9pFnfN2P5FtDhgxB//79OS6qKC0tDQ8++CCCg4Od398oD8dF9fG2L9zhuKi6zp07Y/To0Xj66afxww8/IC8vD+PHjy/3O0k1PS4YoH0oPDwcqamppbYXbYuIiHB7XNOmTWE0Git1LLlX2b4oT5s2bXDz5s0qt408Vx39SL7FcVE1t2/fxv3334/s7Gxs2rTJo99pjovqUZm+KAvHhW9NnDgRCQkJ5dYx1/S4YID2oZ49eyI5OblUneyBAwecj7uj0WjQvXt3HDp0qNRjBw4cQMeOHREYGOjz9jZkle2LsgghcPHiRYSGhvqqieSBnj174siRI1BV1WX7gQMHYDabER0dXUstoyLnz5/nuKgki8WC8ePHIzk5GT/88AO6du3q0XEcF75X2b4oC8eFbxWVX9y+fbvMfWp6XDBA+9DEiROhKAo++eQT5zar1YqVK1eif//+aNOmDQDg8uXLOHPmTKljExISXEL02bNnsXXrVkyaNKlm3kADUpW+yMzMLPV8H330ETIzMzF27NjqbXgjlpqaijNnzsButzu3TZw4Eenp6fj222+d227cuIG1a9di/PjxbuvdqOrc9YW7cbFhwwYcPnyY46ISFEXB5MmTsW/fPqxduxYDBw50ux/HRfWrSl9wXPhWRkZGqW12ux3/+te/YDKZnBc2dWFcSIKL3PpUfHw8vvvuO7z00kuIjIzEZ599hoMHD+Lnn3/G0KFDAQDDhw/Hjh07XGp5cnNz0atXL+Tm5mLOnDnQ6/X44IMPoCgKEhMTeSVbCZXtC7PZjMmTJ6N79+7w8/PD7t278dVXXyE2NhZ79uyB2WyurbdUb3344YfIzs7G9evX8dFHH+HRRx9Fr169AAC/+93vEBwcjGnTpuGzzz7DhQsX0L59ewCO/7ENGTIEJ06ccLmz1OXLl5GQkIBOnTrV4ruqnyrbF1FRUejVqxf69OmD4OBgHDlyBCtWrEB4eDgSEhIQFhZWi++q/nnxxRexdOlSjB8/3uVLaEWmTJkCABwXNaAqfcFx4VuPPPIIcnJyMHToULRq1QppaWn44osvcObMGbz//vv4wx/+AKCOjAuf35qlkSssLBRz5swRLVu2FEajUfTt21ds2rTJZZ9hw4YJd6f+ypUrYuLEiSIoKEgEBASIcePGiZSUlJpqeoNT2b549tlnRdeuXUVgYKDQ6/UiMjJSzJs3T+Tk5NRk8xuUdu3aCQBu/7lw4YIQQoipU6e6/Fzk5s2bYvr06aJZs2bCbDaLYcOG8a5fVVDZvnj11VdFz549RXBwsNDr9aJt27bi+eefF2lpabXzRuq5or89Zf1ThOOi+lWlLzgufGv16tVi9OjRIiwsTOh0OtGkSRMxevRo8f3337vsVxfGBWegiYiIiIi8wBpoIiIiIiIvMEATEREREXmBAZqIiIiIyAsM0EREREREXmCAJiIiIiLyAgM0EREREZEXGKCJiIiIiLzAAE1ERERE5AUGaCIiIiIiLzBAExERERF5gQGaiIiIiMgLDNBERERERF5ggCYiIiIi8gIDNBERVchisdR2E4iI6gwGaCIicnHvvfdi0KBB2LVrF4YNGwaTyYTZs2fXdrOIiOoMSQgharsRRERUd4SFhSEkJAQ3btzAc889hw4dOiAyMhKjRo2q7aYREdUJutpuABER1R0ZGRnIyMhAQUEBEhIS0Llz59puEhFRncMSDiIickpKSgIAvPLKKwzPRERlYIAmIiKn48ePAwAmT55cyy0hIqq7GKCJiMgpKSkJ4eHh6NixY203hYiozmKAJiIip6SkJMTGxtZ2M4iI6jQGaCIiAgAoioJTp04xQBMRVYABmoiIAAApKSmwWCwM0EREFWCAJiIiAMVfIOzRo0ctt4SIqG7jjVSIiIiIiLzAGWgiIiIiIi8wQBMREREReYEBmoiIiIjICwzQREREREReYIAmIiIiIvICAzQRERERkRcYoImIiIiIvMAATURERETkBQZoIiIiIiIvMEATEREREXmBAZqIiIiIyAsM0EREREREXmCAJiIiIiLyAgM0EREREZEXGKCJiIiIiLzAAE1ERERE5AUGaCIiIiIiLzBAExERERF5gQGaiIiIiMgLDNBERERERF74/5odSBv0KngDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import Image\n",
    "\n",
    "r_96,log10H_96 = np.loadtxt(os.path.join(outdir,\"out96.txt\")).T\n",
    "r_72,log10H_72 = np.loadtxt(os.path.join(outdir,\"out72.txt\")).T\n",
    "r_48,log10H_48 = np.loadtxt(os.path.join(outdir,\"out48.txt\")).T\n",
    "\n",
    "# Set 4 subplots\n",
    "fig = plt.figure()\n",
    "\n",
    "plt.title(\"Plot Demonstrating 4th-order Convergence\")\n",
    "plt.ylabel(r\"$\\log_{10}|H|$\")\n",
    "plt.xlabel(r\"$r$\")\n",
    "plt.xlim(0,3)\n",
    "plt.ylim(-10,-5)\n",
    "plt.plot(r_96,log10H_96,color='black',ls='-',label=r\"$N_{r}=96$\")\n",
    "plt.plot(r_72,log10H_72+np.log10((72.0/96.0)**4),color='green',ls='-.',label=r\"$N_{r}=72$, times (72/96)$^4$\")\n",
    "plt.plot(r_48,log10H_48+np.log10((48.0/96.0)**4),color='red',ls='--',label=r\"$N_{r}=60$, times (48/96)$^4$\")\n",
    "legend = plt.legend(shadow=True,loc=3)\n",
    "legend.get_frame().set_facecolor('C1')\n",
    "plt.tight_layout()\n",
    "\n",
    "outfile = os.path.join(outdir,\"HC_convergence_SFID.png\")\n",
    "plt.savefig(outfile,facecolor='white',dpi=120)\n",
    "plt.close(fig)\n",
    "Image(outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 7: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_ScalarField_initial_data.pdf](Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_ScalarField_initial_data.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": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-15T10:07:25.517284Z",
     "iopub.status.busy": "2021-06-15T10:07:25.517010Z",
     "iopub.status.idle": "2021-06-15T10:07:28.188209Z",
     "shell.execute_reply": "2021-06-15T10:07:28.187743Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish-BSSNCurvilinear-\n",
      "    Setting_up_ScalarField_initial_data.tex, and compiled LaTeX file to PDF\n",
      "    file Tutorial-Start_to_Finish-BSSNCurvilinear-\n",
      "    Setting_up_ScalarField_initial_data.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-Setting_up_ScalarField_initial_data\")"
   ]
  }
 ],
 "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": 4
}