Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide-Psi4_new_way.ipynb

{
 "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: Head-On Black Hole Collision Gravitational Wave Ringdown\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module implements a basic numerical relativity code to merge two black holes in *spherical-like coordinates*, as well as the gravitational wave analysis provided by the $\\psi_4$ NRPy+ tutorial notebooks ([$\\psi_4$](Tutorial-Psi4.ipynb) & [$\\psi_4$ tetrad](Tutorial-Psi4_tetrads.ipynb)).\n",
    "\n",
    "### Here we place the black holes initially on the $z$-axis, so the entire simulation is axisymmetric about the $\\phi$-axis. Minimal sampling in the $\\phi$ direction greatly speeds up the simulation.\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 *after a short numerical evolution of the initial data* (see [plots at bottom](#convergence)), and all quantities have been validated against the [original SENR code](https://bitbucket.org/zach_etienne/nrpy). Finally, excellent agreement is seen in the gravitational wave signal from the ringing remnant black hole for multiple decades in amplitude when compared to black hole perturbation theory predictions of the exponential dropoff and frequency.\n",
    "\n",
    "### NRPy+ Source Code & Tutorials for this module: \n",
    "* [BSSN/BrillLindquist.py](../edit/BSSN/BrillLindquist.py); [\\[**tutorial**\\]](Tutorial-ADM_Initial_Data-Brill-Lindquist.ipynb): Brill-Lindquist initial data; sets all ADM variables in Cartesian basis: \n",
    "* [BSSN/ADM_Exact_Spherical_or_Cartesian_to_BSSNCurvilinear.py](../edit/BSSN/ADM_Exact_Spherical_or_Cartesian_to_BSSNCurvilinear.py); [\\[**tutorial**\\]](Tutorial-ADM_Initial_Data-Converting_Exact_ADM_Spherical_or_Cartesian_to_BSSNCurvilinear.ipynb): Spherical/Cartesian ADM$\\to$Curvilinear BSSN converter function, for which exact expressions are given for ADM quantities.\n",
    "* [BSSN/BSSN_ID_function_string.py](../edit/BSSN/BSSN_ID_function_string.py): Sets up the C code string enabling initial data be set up in a point-by-point fashion\n",
    "* [BSSN/BSSN_Ccodegen_library.py](../edit/BSSN/BSSN_Ccodegen_library.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-C_codegen_library.ipynb): Implements a number of helper functions for generating C codes from symbolic expressions generated in the following modules/tutorials:\n",
    "    * [BSSN/BSSN_constraints.py](../edit/BSSN/BSSN_constraints.py); [\\[**tutorial**\\]](Tutorial-BSSN_constraints.ipynb): Hamiltonian constraint in BSSN curvilinear basis/coordinates\n",
    "    * [BSSN/BSSN_RHSs.py](../edit/BSSN/BSSN_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb): Generates the right-hand sides for the BSSN evolution equations in singular, curvilinear coordinates\n",
    "    * [BSSN/BSSN_gauge_RHSs.py](../edit/BSSN/BSSN_gauge_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb): Generates the right-hand sides for the BSSN gauge evolution equations in singular, curvilinear coordinates\n",
    "    * [BSSN/Psi4.py](../edit/BSSN/Psi4.py); [\\[**tutorial**\\]](Tutorial-Psi4.ipynb): Generates expressions for $\\psi_4$, the outgoing Weyl scalar\n",
    "        + [BSSN/Psi4_tetrads.py](../edit/BSSN/Psi4_tetrads.py); [\\[**tutorial**\\]](Tutorial-Psi4_tetrads.ipynb): Generates quasi-Kinnersley tetrad needed for $\\psi_4$-based gravitational wave extraction\n",
    "\n",
    "\n",
    "## Introduction:\n",
    "Here we use NRPy+ to generate the C source code necessary to set up initial data for two black holes (Brill-Lindquist, [Brill & Lindquist, Phys. Rev. 131, 471, 1963](https://journals.aps.org/pr/abstract/10.1103/PhysRev.131.471); see also Eq. 1 of [Brandt & Brügmann, arXiv:gr-qc/9711015v1](https://arxiv.org/pdf/gr-qc/9711015v1.pdf)). Then we use it to generate the RHS expressions for [Method of Lines](https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html) time integration based on an [explicit Runge-Kutta fourth-order scheme](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (RK4 is chosen below, but multiple options exist). \n",
    "\n",
    "The entire algorithm is outlined as follows, with links to the relevant NRPy+ tutorial notebooks listed at each step:\n",
    "\n",
    "1. Allocate memory for gridfunctions, including temporary storage for the Method of Lines time integration\n",
    "    * [**NRPy+ tutorial on Method of Lines algorithm**](Tutorial-Method_of_Lines-C_Code_Generation.ipynb).\n",
    "1. Set gridfunction values to initial data \n",
    "    * [**NRPy+ tutorial on Brill-Lindquist initial data**](Tutorial-ADM_Initial_Data-Brill-Lindquist.ipynb)\n",
    "    * [**NRPy+ tutorial on validating Brill-Lindquist initial data**](Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_Exact_Initial_Data.ipynb).\n",
    "1. Next, integrate the initial data forward in time using the Method of Lines coupled to a Runge-Kutta explicit timestepping algorithm:\n",
    "    1. At the start of each iteration in time, output the Hamiltonian constraint violation \n",
    "        * [**NRPy+ tutorial on BSSN constraints**](Tutorial-BSSN_constraints.ipynb).\n",
    "    1. At each RK time substep, do the following:\n",
    "        1. Evaluate BSSN RHS expressions \n",
    "            * [**NRPy+ tutorial on BSSN right-hand sides**](Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb) ([**BSSN Introduction Notebook**](Tutorial-BSSN_formulation.ipynb))\n",
    "            * [**NRPy+ tutorial on BSSN gauge condition right-hand sides**](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb) \n",
    "        1. Apply singular, curvilinear coordinate boundary conditions [*a la* the SENR/NRPy+ paper](https://arxiv.org/abs/1712.07658)\n",
    "            * [**NRPy+ tutorial on setting up singular, curvilinear boundary conditions**](Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb)\n",
    "        1. Enforce constraint on conformal 3-metric: $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ \n",
    "            * [**NRPy+ tutorial on enforcing $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint**](Tutorial-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.ipynb)\n",
    "1. Repeat above steps at two numerical resolutions to confirm convergence to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#initializenrpy): Set core NRPy+ parameters for numerical grids and reference metric\n",
    "1. [Step 2](#ccodegen): Generate C code kernels for BSSN expressions, in parallel if possible\n",
    "    1. [Step 2.a](#rfm_ccodegen): Generate C code kernels for reference metric\n",
    "1. [Step 3](#cparams_rfm_and_domainsize): Set `free_parameters.h`; also output C codes needed for declaring and setting Cparameters\n",
    "1. [Step 4](#bc_functs): Set up boundary condition functions for chosen singular, curvilinear coordinate system\n",
    "1. [Step 5](#mainc): `BrillLindquist_Playground`: The C code `main()` function\n",
    "1. [Step 6](#compileexec): Compile generated C codes & perform the black hole collision calculation\n",
    "1. [Step 7](#visualize): Visualize the output!\n",
    "    1. [Step 7.a](#installdownload): Install `scipy` and download `ffmpeg` if they are not yet installed/downloaded\n",
    "    1. [Step 7.b](#genimages): Generate images for visualization animation\n",
    "    1. [Step 7.c](#genvideo): Generate visualization animation\n",
    "1. [Step 8](#numerror): Plot the numerical error at the end of the simulation\n",
    "1. [Step 9](#compare): Compare ringdown gravitational wave signal with black hole perturbation theory\n",
    "1. [Step 10](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Set core NRPy+ parameters for numerical grids and reference metric \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:28:51.090794Z",
     "iopub.status.busy": "2021-09-23T04:28:51.084185Z",
     "iopub.status.idle": "2021-09-23T04:28:51.411410Z",
     "shell.execute_reply": "2021-09-23T04:28:51.411799Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step P1: Import needed NRPy+ core modules:\n",
    "from outputC import lhrh, outCfunction, outC_function_dict, add_to_Cfunction_dict # NRPy+: Core C code output module\n",
    "import finite_difference as fin  # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface\n",
    "from pickling import pickle_NRPy_env, unpickle_NRPy_env  # NRPy+: Pickle/unpickle NRPy+ environment, for parallel codegen\n",
    "import shutil, os, sys, time     # Standard Python modules for multiplatform OS-level functions, benchmarking\n",
    "\n",
    "EvolOption = \"low resolution\"\n",
    "\n",
    "if EvolOption == \"low resolution\":\n",
    "    FD_order  = 8  # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
    "    # domain_size sets the default value for:\n",
    "    #   * Spherical's params.RMAX\n",
    "    #   * SinhSpherical*'s params.AMAX\n",
    "    #   * Cartesians*'s -params.{x,y,z}min & .{x,y,z}max\n",
    "    #   * Cylindrical's -params.ZMIN & .{Z,RHO}MAX\n",
    "    #   * SinhCylindrical's params.AMPL{RHO,Z}\n",
    "    #   * *SymTP's params.AMAX\n",
    "    domain_size = 150\n",
    "    Nr          = 270\n",
    "    Nth         = 8\n",
    "    Nph         = 2\n",
    "    CFL_FACTOR  = 1.0  # In pure axisymmetry (symmetry_axes = 2 below) 1.0 works fine. Otherwise 0.5 or lower.\n",
    "elif EvolOption == \"high resolution\":\n",
    "    FD_order  = 10  # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
    "    # domain_size sets the default value for:\n",
    "    #   * Spherical's params.RMAX\n",
    "    #   * SinhSpherical*'s params.AMAX\n",
    "    #   * Cartesians*'s -params.{x,y,z}min & .{x,y,z}max\n",
    "    #   * Cylindrical's -params.ZMIN & .{Z,RHO}MAX\n",
    "    #   * SinhCylindrical's params.AMPL{RHO,Z}\n",
    "    #   * *SymTP's params.AMAX\n",
    "    domain_size = 300\n",
    "    Nr          = 800\n",
    "    Nth         = 16\n",
    "    Nph         = 2\n",
    "    CFL_FACTOR  = 1.0  # In pure axisymmetry (symmetry_axes = 2 below) 1.0 works fine. Otherwise 0.5 or lower.\n",
    "else:\n",
    "    print(\"ERROR: EvolOption = \" + EvolOption + \"unknown.\")\n",
    "    sys.exit(1)\n",
    "\n",
    "# Step P2: Create C code output directory:\n",
    "Ccodesrootdir = os.path.join(\"BSSN_Two_BHs_Collide_Ccodes_psi4_new_way\")\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(Ccodesrootdir, ignore_errors=True)\n",
    "# Then create a fresh directory\n",
    "cmd.mkdir(Ccodesrootdir)\n",
    "\n",
    "# Step P3: Create executable output directory:\n",
    "outdir = os.path.join(Ccodesrootdir, \"output\")\n",
    "cmd.mkdir(outdir)\n",
    "\n",
    "# Step 1.a: Enable SIMD-optimized code?\n",
    "#           I.e., generate BSSN and Ricci C code kernels using SIMD-vectorized\n",
    "#           compiler intrinsics, which *greatly improve the code's performance*,\n",
    "#           though at the expense of making the C-code kernels less\n",
    "#           human-readable.\n",
    "#           * Important note in case you wish to modify the BSSN/Ricci kernels\n",
    "#             here by adding expressions containing transcendental functions\n",
    "#             (e.g., certain scalar fields):\n",
    "#           Note that SIMD-based transcendental function intrinsics are not\n",
    "#           supported by the default installation of gcc or clang (you will\n",
    "#           need to use e.g., the SLEEF library from sleef.org, for this\n",
    "#           purpose). The Intel compiler suite does support these intrinsics\n",
    "#           however without the need for external libraries.\n",
    "enable_SIMD = True\n",
    "\n",
    "# Step 1.b: Enable reference metric precomputation.\n",
    "enable_rfm_precompute = True\n",
    "\n",
    "if enable_SIMD and not enable_rfm_precompute:\n",
    "    print(\"ERROR: SIMD does not currently handle transcendental functions,\\n\")\n",
    "    print(\"       like those found in rfmstruct (rfm_precompute).\\n\")\n",
    "    print(\"       Therefore, enable_SIMD==True and enable_rfm_precompute==False\\n\")\n",
    "    print(\"       is not supported.\\n\")\n",
    "    sys.exit(1)\n",
    "\n",
    "# Step 1.c: Enable \"FD functions\". In other words, all finite-difference stencils\n",
    "#         will be output as inlined static functions. This is essential for\n",
    "#         compiling highly complex FD kernels with using certain versions of GCC;\n",
    "#         GCC 10-ish will choke on BSSN FD kernels at high FD order, sometimes\n",
    "#         taking *hours* to compile. Unaffected GCC versions compile these kernels\n",
    "#         in seconds. FD functions do not slow the code performance, but do add\n",
    "#         another header file to the C source tree.\n",
    "# With gcc 7.5.0, enable_FD_functions=True decreases performance by 10%\n",
    "enable_FD_functions = False\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",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\", CoordSystem)\n",
    "rfm.reference_metric()\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",
    "\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",
    "# sinh_width      = 0.4 # If Sinh* coordinates chosen\n",
    "\n",
    "# sinhv2_const_dr sets the default value for:\n",
    "#   * SinhSphericalv2's params.const_dr\n",
    "#   * SinhCylindricalv2's params.const_d{rho,z}\n",
    "sinhv2_const_dr = 0.05 # If Sinh*v2 coordinates chosen\n",
    "\n",
    "# SymTP_bScale sets the default value for:\n",
    "#   * SinhSymTP's params.bScale\n",
    "SymTP_bScale    = 0.5 # If SymTP chosen\n",
    "\n",
    "# Step 2.b: Set the order of spatial and temporal derivatives;\n",
    "#           the core data type, and the CFL factor.\n",
    "# RK_method choices include: Euler, \"RK2 Heun\", \"RK2 MP\", \"RK2 Ralston\", RK3, \"RK3 Heun\", \"RK3 Ralston\",\n",
    "#              SSPRK3, RK4, DP5, DP5alt, CK5, DP6, L6, DP8\n",
    "RK_method = \"RK4\"\n",
    "maximum_l = 2\n",
    "REAL      = \"double\"  # Best to use double here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:28:51.417482Z",
     "iopub.status.busy": "2021-09-23T04:28:51.416973Z",
     "iopub.status.idle": "2021-09-23T04:28:51.418903Z",
     "shell.execute_reply": "2021-09-23T04:28:51.418495Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 5: Set the finite differencing order to FD_order (set above).\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", FD_order)\n",
    "\n",
    "# Directory for reference_metric precomputation header files:\n",
    "rfm_precompute_Ccode_outdir = os.path.join(Ccodesrootdir, \"rfm_files/\")\n",
    "if enable_rfm_precompute:\n",
    "    cmd.mkdir(os.path.join(Ccodesrootdir, \"rfm_files/\"))\n",
    "    par.set_parval_from_str(\"reference_metric::rfm_precompute_Ccode_outdir\", rfm_precompute_Ccode_outdir)\n",
    "\n",
    "# Step 6: Copy SIMD/SIMD_intrinsics.h to $Ccodesrootdir/SIMD/SIMD_intrinsics.h\n",
    "if enable_SIMD:\n",
    "    cmd.mkdir(os.path.join(Ccodesrootdir,\"SIMD\"))\n",
    "    shutil.copy(os.path.join(\"SIMD/\")+\"SIMD_intrinsics.h\",os.path.join(Ccodesrootdir,\"SIMD/\"))\n",
    "\n",
    "# Step 7: Set finite_difference::enable_FD_functions appropriately. Defaults to False\n",
    "if enable_FD_functions:\n",
    "    par.set_parval_from_str(\"finite_difference::enable_FD_functions\", enable_FD_functions)\n",
    "\n",
    "# Step 8: If enable_SIMD, then copy SIMD/SIMD_intrinsics.h to $Ccodesrootdir/SIMD/SIMD_intrinsics.h\n",
    "cmd.mkdir(os.path.join(Ccodesrootdir,\"SIMD\"))\n",
    "if enable_SIMD:\n",
    "    shutil.copy(os.path.join(\"SIMD\", \"SIMD_intrinsics.h\"), os.path.join(Ccodesrootdir, \"SIMD\"))\n",
    "\n",
    "# Step 9: Set the direction=2 (phi) axis to be the symmetry axis; i.e.,\n",
    "#         axis \"2\", corresponding to the i2 direction.\n",
    "#      This sets all spatial derivatives in the phi direction to zero.\n",
    "par.set_parval_from_str(\"indexedexp::symmetry_axes\",\"2\")\n",
    "OMP_pragma_on = \"i1\"  # structure OpenMP loops to parallelize, not over i2 (phi direction),\n",
    "#                       but i1 (theta direction). Results in *very* significant speed-up!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:28:51.422047Z",
     "iopub.status.busy": "2021-09-23T04:28:51.421615Z",
     "iopub.status.idle": "2021-09-23T04:28:51.432342Z",
     "shell.execute_reply": "2021-09-23T04:28:51.432729Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 10: Generate Runge-Kutta-based (RK-based) timestepping code.\n",
    "#       As described above the Table of Contents, this is a 3-step process:\n",
    "#       3.A: Evaluate RHSs (RHS_string)\n",
    "#       3.B: Apply boundary conditions (post_RHS_string, pt 1)\n",
    "#       3.C: Enforce det(gammabar) = det(gammahat) constraint (post_RHS_string, pt 2)\n",
    "import MoLtimestepping.MoL_new_way as MoL\n",
    "# from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict\n",
    "# RK_order  = Butcher_dict[RK_method][1]\n",
    "RHS_string      = \"\"\"\n",
    "Ricci_eval(params, rfmstruct, RK_INPUT_GFS, auxevol_gfs);\n",
    "rhs_eval(params, rfmstruct, auxevol_gfs, RK_INPUT_GFS, RK_OUTPUT_GFS);\"\"\"\n",
    "post_RHS_string = \"\"\"\n",
    "apply_bcs_curvilinear(params, bcstruct, NUM_EVOL_GFS, evol_gf_parity, RK_OUTPUT_GFS);\n",
    "enforce_detgammahat_constraint(params, rfmstruct,                     RK_OUTPUT_GFS);\\n\"\"\"\n",
    "\n",
    "if not enable_rfm_precompute:\n",
    "    RHS_string = RHS_string.replace(\"rfmstruct\", \"xx\")\n",
    "    post_RHS_string = post_RHS_string.replace(\"rfmstruct\", \"xx\")\n",
    "\n",
    "MoL.register_C_functions_and_NRPy_basic_defines(RK_method,\n",
    "        RHS_string=RHS_string, post_RHS_string=post_RHS_string,\n",
    "        enable_rfm=enable_rfm_precompute, enable_curviBCs=True, enable_SIMD=enable_SIMD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ccodegen'></a>\n",
    "\n",
    "# Step 2: Generate C code kernels for BSSN expressions, in parallel if possible \\[Back to [top](#toc)\\]\n",
    "$$\\label{ccodegen}$$\n",
    "\n",
    "In the following code cell, we create a list of Python functions, which each registers a single C code function in `outputC`'s `outC_function_dict` dictionary. These Python functions are defined in \n",
    "\n",
    "1. the [`BSSN.BrillLindquist`](../edit/BSSN/BrillLindquist.py) NRPy+ module, which does the following:\n",
    "    1. Sets up Brill-Lindquist initial data [ADM](https://en.wikipedia.org/wiki/ADM_formalism) quantities in the **Cartesian basis**, as [documented here](Tutorial-ADM_Initial_Data-Brill-Lindquist.ipynb). \n",
    "    1. Converts the ADM **Cartesian quantities** to **BSSN quantities in the desired Curvilinear basis** (set by reference_metric::CoordSystem), as [documented here](Tutorial-ADM_Initial_Data-Converting_ADMCartesian_to_BSSNCurvilinear.ipynb).\n",
    "    1. Sets up the standardized C function for setting all BSSN Curvilinear gridfunctions in a pointwise fashion, as [written here](../edit/BSSN/BSSN_ID_function_string.py), and returns the C function as a Python string.\n",
    "1. the [`BSSN.BSSN_Ccodegen_library`](../edit/BSSN/BSSN_Ccodegen_library.py) NRPy+ module [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-C_codegen_library.ipynb), which contains Python functions for generating C code from symbolic expressions constructed within the following NRPy+ modules/tutorials:\n",
    "    1. [BSSN/BSSN_constraints.py](../edit/BSSN/BSSN_constraints.py); [\\[**tutorial**\\]](Tutorial-BSSN_constraints.ipynb): Hamiltonian constraint in BSSN curvilinear basis/coordinates\n",
    "    1. [BSSN/BSSN_RHSs.py](../edit/BSSN/BSSN_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb): Generates the right-hand sides for the BSSN evolution equations in singular, curvilinear coordinates\n",
    "    1. [BSSN/BSSN_gauge_RHSs.py](../edit/BSSN/BSSN_gauge_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb): Generates the right-hand sides for the BSSN gauge evolution equations in singular, curvilinear coordinates\n",
    "    1. [BSSN/Enforce_Detgammahat_Constraint.py](../edit/BSSN/Enforce_Detgammahat_Constraint.py); [**tutorial**](Tutorial-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.ipynb): Generates symbolic expressions for enforcing the $\\det{\\bar{\\gamma}}=\\det{\\hat{\\gamma}}$ constraint\n",
    "\n",
    "\n",
    "Next, from within a `multiprocessing` environment, we then call all the Python C-code generation functions in this list in parallel (if `multiprocessing` is supported). This is quite useful, as these functions take several seconds to complete.\n",
    "\n",
    "Within each `multiprocessing` process, the current NRPy+ environment is cloned, and a new function is registered to the `outC_function_dict` dictionary. Thus when each process completes, it contains a unique NRPy+ environment, with only its function registered. We address this by saving each process' NRPy+ environment and sending it back in a common binary format known as a `pickle`, using NRPy+'s [`pickling`](../edit/pickling.py) module. The environments are combined in an unpickling such that all functions exist within the same `outC_function_dict` dictionary.\n",
    "\n",
    "To make the current environment fully consistent, we call `reference_metric.py` to register all its associated C functions (stored in globals) and contributions to `NRPy_basic_defines.h`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:28:51.438717Z",
     "iopub.status.busy": "2021-09-23T04:28:51.438212Z",
     "iopub.status.idle": "2021-09-23T04:29:02.701767Z",
     "shell.execute_reply": "2021-09-23T04:29:02.701515Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating C code for psi4, part 1 (SinhSpherical coords)...Generating C code for Spin-weight s=-2 Spherical Harmonics (SinhSpherical coords)...Generating C code for psi4 tetrads (SinhSpherical coords)...Generating C code for psi4, part 0 (SinhSpherical coords)...Generating C code for psi4, part 2 (SinhSpherical coords)...\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "Finished generating C code for Spin-weight s=-2 Spherical Harmonics (SinhSpherical coords) in 0.3 seconds.\n",
      "Generating symbolic expressions for BSSN_RHSs (SinhSpherical coords)...\n",
      "Generating symbolic expressions for BSSN constraints (SinhSpherical coords)...\n",
      "Generating symbolic expressions for 3-Ricci tensor (SinhSpherical coords)...\n",
      "Generating C code for Enforcing det(gammabar)=det(gammahat) constraint (SinhSpherical coords)...\n",
      "Finished generating C code for Enforcing det(gammabar)=det(gammahat) constraint (SinhSpherical coords) in 0.1 seconds.\n",
      "Finished generating symbolic expressions for 3-Ricci tensor (SinhSpherical coords) in 0.3 seconds. Next up: C codegen...\n",
      "Generating C code for 3-Ricci tensor (SinhSpherical coords)...\n",
      "Finished generating symbolic expressions for BSSN constraints (SinhSpherical coords) in 0.5 seconds. Next up: C codegen...\n",
      "Generating C code for BSSN constraints (SinhSpherical coords)...\n",
      "Finished generating symbolic expressions for BSSN_RHSs (SinhSpherical coords) in 1.0 seconds. Next up: C codegen...\n",
      "Generating C code for BSSN_RHSs (SinhSpherical coords)...\n",
      "Finished generating C code for BSSN constraints (SinhSpherical coords) in 2.1 seconds.\n",
      "Finished generating C code for psi4, part 1 (SinhSpherical coords) in 3.5 seconds.\n",
      "Finished generating C code for psi4 tetrads (SinhSpherical coords) in 11.2 seconds.\n",
      "Finished generating C code for psi4, part 0 (SinhSpherical coords) in 15.3 seconds.\n",
      "Finished generating C code for BSSN_RHSs (SinhSpherical coords) in 15.4 seconds.\n",
      "Finished generating C code for psi4, part 2 (SinhSpherical coords) in 16.9 seconds.\n",
      "Finished generating C code for 3-Ricci tensor (SinhSpherical coords) in 16.8 seconds.\n"
     ]
    }
   ],
   "source": [
    "# Step 2: Generate C code kernels for BSSN expressions, in parallel if possible;\n",
    "import BSSN.BSSN_Ccodegen_library as BCL\n",
    "import BSSN.BrillLindquist as bl\n",
    "# Step 2.a: Create a list of functions we wish to evaluate in parallel (if possible)\n",
    "# Create lists for all BSSN functions\n",
    "BSSN_funcs = [bl.BrillLindquist]\n",
    "BSSN_funcs.append(BCL.add_rhs_eval_to_Cfunction_dict)\n",
    "BSSN_funcs.append(BCL.add_Ricci_eval_to_Cfunction_dict)\n",
    "BSSN_funcs.append(BCL.add_BSSN_constraints_to_Cfunction_dict)\n",
    "BSSN_funcs.append(BCL.add_enforce_detgammahat_constraint_to_Cfunction_dict)\n",
    "BSSN_funcs.append(BCL.add_psi4_part_to_Cfunction_dict)  # part 0 of 2\n",
    "BSSN_funcs.append(BCL.add_psi4_part_to_Cfunction_dict)  # part 1 of 2\n",
    "BSSN_funcs.append(BCL.add_psi4_part_to_Cfunction_dict)  # part 2 of 2\n",
    "BSSN_funcs.append(BCL.add_psi4_tetrad_to_Cfunction_dict)\n",
    "BSSN_funcs.append(BCL.add_SpinWeight_minus2_SphHarmonics_to_Cfunction_dict)\n",
    "\n",
    "# Step 2.b: Define master functions for parallelization.\n",
    "#           Note that lambdifying this doesn't work in Python 3\n",
    "def master_func(arg):\n",
    "    if BSSN_funcs[arg] == bl.BrillLindquist:\n",
    "        ret = BSSN_funcs[arg](include_NRPy_basic_defines_and_pickle=True)\n",
    "    elif BSSN_funcs[arg].__name__ == \"add_psi4_part_to_Cfunction_dict\":\n",
    "        which_arg_is_part1 = -1\n",
    "        for i, func in enumerate(BSSN_funcs):\n",
    "            if func.__name__ == \"add_psi4_part_to_Cfunction_dict\":\n",
    "                which_arg_is_part1 = i\n",
    "                break\n",
    "        ret = BSSN_funcs[arg]( includes=[\"NRPy_basic_defines.h\"],\n",
    "                               rel_path_to_Cparams=os.path.join(\".\"), whichpart=(arg-which_arg_is_part1),\n",
    "                               OMP_pragma_on=OMP_pragma_on)\n",
    "    elif BSSN_funcs[arg].__name__ == \"add_psi4_tetrad_to_Cfunction_dict\":\n",
    "        ret = BSSN_funcs[arg](includes=[\"NRPy_basic_defines.h\"],\n",
    "                              rel_path_to_Cparams=os.path.join(\".\"))\n",
    "    elif BSSN_funcs[arg].__name__ == \"add_SpinWeight_minus2_SphHarmonics_to_Cfunction_dict\":\n",
    "        ret = BSSN_funcs[arg](includes=[\"NRPy_basic_defines.h\"],\n",
    "                              rel_path_to_Cparams=os.path.join(\".\"), maximum_l=maximum_l)\n",
    "    else:\n",
    "        if enable_rfm_precompute:\n",
    "            # We use rfm_precompute for all BSSN functions:\n",
    "            par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"True\")\n",
    "            rfm.reference_metric()\n",
    "\n",
    "        if BSSN_funcs[arg].__name__ == \"add_BSSN_constraints_to_Cfunction_dict\":\n",
    "            ret = BSSN_funcs[arg](includes=[\"NRPy_basic_defines.h\"],\n",
    "                                  rel_path_to_Cparams=os.path.join(\".\"), output_H_only=True,\n",
    "                                  enable_rfm_precompute=enable_rfm_precompute, enable_SIMD=enable_SIMD,\n",
    "                                  OMP_pragma_on=OMP_pragma_on)\n",
    "        elif BSSN_funcs[arg].__name__ == \"add_rhs_eval_to_Cfunction_dict\" or \\\n",
    "            BSSN_funcs[arg].__name__ == \"add_Ricci_eval_to_Cfunction_dict\":\n",
    "            ret = BSSN_funcs[arg](includes=[\"NRPy_basic_defines.h\"],\n",
    "                                  rel_path_to_Cparams=os.path.join(\".\"),\n",
    "                                  enable_rfm_precompute=enable_rfm_precompute, enable_SIMD=enable_SIMD,\n",
    "                                  OMP_pragma_on=OMP_pragma_on)\n",
    "        elif BSSN_funcs[arg].__name__ == \"add_enforce_detgammahat_constraint_to_Cfunction_dict\":\n",
    "            ret = BSSN_funcs[arg](includes=[\"NRPy_basic_defines.h\"],\n",
    "                                  rel_path_to_Cparams=os.path.join(\".\"),\n",
    "                                  enable_rfm_precompute=enable_rfm_precompute, OMP_pragma_on=OMP_pragma_on)\n",
    "        else:\n",
    "            print(\"ERROR: DID NOT RECOGNIZE FUNCTION \" + BSSN_funcs[arg].__name__ + \"\\n\")\n",
    "            sys.exit(1)\n",
    "        if enable_rfm_precompute:\n",
    "            par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"False\")\n",
    "            rfm.ref_metric__hatted_quantities()\n",
    "    return ret\n",
    "\n",
    "\n",
    "NRPyEnvVars = []\n",
    "raised_exception = False\n",
    "try:\n",
    "    if os.name == 'nt':\n",
    "        # It's a mess to get working in Windows, so we don't bother. :/\n",
    "        #  https://medium.com/@grvsinghal/speed-up-your-python-code-using-multiprocessing-on-windows-and-jupyter-or-ipython-2714b49d6fac\n",
    "        raise Exception(\"Parallel codegen currently not available in certain environments, e.g., Windows\")\n",
    "\n",
    "    # Step 2.d: Import the multiprocessing module.\n",
    "    import multiprocessing\n",
    "\n",
    "    # Step 2.e: Evaluate list of functions in parallel if possible;\n",
    "    #           otherwise fallback to serial evaluation:\n",
    "    pool = multiprocessing.Pool()\n",
    "    NRPyEnvVars.append(pool.map(master_func, range(len(BSSN_funcs))))\n",
    "    pool.terminate()\n",
    "    pool.join()\n",
    "except:\n",
    "    print(\"FAILED PARALLEL CODEGEN!\")\n",
    "    NRPyEnvVars = [] # Reset, as pickling/unpickling unnecessary for serial codegen (see next line)\n",
    "\n",
    "    # Steps 2.d-e, alternate: As fallback, evaluate functions in serial.\n",
    "    #       This will happen on Android and Windows systems\n",
    "    for i, func in enumerate(BSSN_funcs):\n",
    "        master_func(i)\n",
    "    raised_exception = True\n",
    "\n",
    "if not raised_exception:\n",
    "    unpickle_NRPy_env(NRPyEnvVars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rfm_ccodegen'></a>\n",
    "\n",
    "## Step 2.a: Generate C code kernels for reference metric \\[Back to [top](#toc)\\]\n",
    "$$\\label{rfm_ccodegen}$$\n",
    "\n",
    "In the [reference_metric](../edit/reference_metric.py) NRPy+ module, `register_C_functions_and_NRPy_basic_defines()` registers the following C functions to `outC_Cfunction_dict`:\n",
    "\n",
    "1. `find_timestep()`: Finds the minimum spacing between adjacent gridpoints on our numerical grid $\\min(ds_i)$, and sets the timestep according to the [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673) condition: $\\Delta t \\le \\frac{\\min(ds_i)}{c}$, where $c$ is the wavespeed, and $ds_i = h_i \\Delta x^i$ is the proper distance between neighboring gridpoints in the $i$th direction (in 3D, there are 3 directions), $h_i$ is the $i$th reference metric scale factor, and $\\Delta x^i$ is the uniform grid spacing in the $i$th direction.\n",
    "1. `xx_to_Cart()`: Input = uniformly sampled coordinate xx0,xx1,xx2 (e.g., r,theta,phi in Spherical coordinates). Output = Cartesian coordinate (x,y,z).\n",
    "1. `set_Nxx_dxx_invdx_params__and__xx()`: Sets `Nxx{0,1,2}`, `Nxx_plus_2NGHOSTS{0,1,2}`, `dxx{0,1,2}`, and `invdx{0,1,2}`; and defines `xx[3][]`.\n",
    "1. `Cart_to_xx_and_nearest_i0i1i2()`: Input = Cartesian coordinate (x,y,z). Output = uniformly sampled coordinate xx0,xx1,xx2 (e.g., r,theta,phi in Spherical coordinates), as well as corresponding grid index `i0,i1,i2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:02.713902Z",
     "iopub.status.busy": "2021-09-23T04:29:02.708519Z",
     "iopub.status.idle": "2021-09-23T04:29:03.025608Z",
     "shell.execute_reply": "2021-09-23T04:29:03.025204Z"
    }
   },
   "outputs": [],
   "source": [
    "# Generate & register C function set_Nxx_dxx_invdx_params__and__xx()\n",
    "# Generate & register C function xx_to_Cart() for\n",
    "#               (the mapping from xx->Cartesian) for the chosen\n",
    "#               CoordSystem:\n",
    "# Generate & register the find_timestep() function\n",
    "\n",
    "# Sets reference_metric globals: NRPy_basic_defines_str, rfm_struct__malloc, rfm_struct__define, rfm_struct__freemem\n",
    "if enable_rfm_precompute:\n",
    "    par.set_parval_from_str(\"reference_metric::rfm_precompute_Ccode_outdir\", rfm_precompute_Ccode_outdir)\n",
    "    par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"True\")\n",
    "    par.set_parval_from_str(\"reference_metric::rfm_precompute_to_Cfunctions_and_NRPy_basic_defines\", \"True\")\n",
    "    rfm.reference_metric()\n",
    "\n",
    "rfm.register_C_functions_and_NRPy_basic_defines(enable_rfm_precompute=enable_rfm_precompute,\n",
    "                                                use_unit_wavespeed_for_find_timestep=True)\n",
    "\n",
    "if enable_rfm_precompute:\n",
    "    par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"False\")\n",
    "    rfm.ref_metric__hatted_quantities()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cparams_rfm_and_domainsize'></a>\n",
    "\n",
    "# Step 3: Set `free_parameters.h`; also output C codes needed for declaring and setting Cparameters  \\[Back to [top](#toc)\\]\n",
    "$$\\label{cparams_rfm_and_domainsize}$$\n",
    "\n",
    "First 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": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:03.029927Z",
     "iopub.status.busy": "2021-09-23T04:29:03.029424Z",
     "iopub.status.idle": "2021-09-23T04:29:03.030447Z",
     "shell.execute_reply": "2021-09-23T04:29:03.030838Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.e.i: Set free_parameters.h\n",
    "with open(os.path.join(Ccodesrootdir,\"free_parameters.h\"),\"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "// Set free-parameter values.\n",
    "\n",
    "// Set the default CFL Factor. Can be overwritten at command line.\n",
    "REAL CFL_FACTOR = \"\"\"+str(CFL_FACTOR)+r\"\"\";\n",
    "\n",
    "// Set free-parameter values for BSSN evolution:\n",
    "params.eta = 2.0;\n",
    "//params.eta = 1.0;\n",
    "\n",
    "// Set free parameters for the (Brill-Lindquist) initial data\n",
    "params.BH1_posn_x = 0.0; params.BH1_posn_y = 0.0; params.BH1_posn_z =+0.25;\n",
    "params.BH2_posn_x = 0.0; params.BH2_posn_y = 0.0; params.BH2_posn_z =-0.25;\n",
    "//params.BH1_posn_x = 0.0; params.BH1_posn_y = 0.0; params.BH1_posn_z =+0.5;\n",
    "//params.BH2_posn_x = 0.0; params.BH2_posn_y = 0.0; params.BH2_posn_z =-0.5;\n",
    "params.BH1_mass = 0.5; params.BH2_mass = 0.5;\n",
    "\"\"\")\n",
    "\n",
    "# Append to $Ccodesrootdir/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(Ccodesrootdir,\"free_parameters.h\"),\n",
    "                                        domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bc_functs'></a>\n",
    "\n",
    "# Step 4: 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_new_way.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:03.036844Z",
     "iopub.status.busy": "2021-09-23T04:29:03.036620Z",
     "iopub.status.idle": "2021-09-23T04:29:03.201578Z",
     "shell.execute_reply": "2021-09-23T04:29:03.201371Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evolved gridfunction \"aDD00\" has parity type 4.\n",
      "Evolved gridfunction \"aDD01\" has parity type 5.\n",
      "Evolved gridfunction \"aDD02\" has parity type 6.\n",
      "Evolved gridfunction \"aDD11\" has parity type 7.\n",
      "Evolved gridfunction \"aDD12\" has parity type 8.\n",
      "Evolved gridfunction \"aDD22\" has parity type 9.\n",
      "Evolved gridfunction \"alpha\" has parity type 0.\n",
      "Evolved gridfunction \"betU0\" has parity type 1.\n",
      "Evolved gridfunction \"betU1\" has parity type 2.\n",
      "Evolved gridfunction \"betU2\" has parity type 3.\n",
      "Evolved gridfunction \"cf\" has parity type 0.\n",
      "Evolved gridfunction \"hDD00\" has parity type 4.\n",
      "Evolved gridfunction \"hDD01\" has parity type 5.\n",
      "Evolved gridfunction \"hDD02\" has parity type 6.\n",
      "Evolved gridfunction \"hDD11\" has parity type 7.\n",
      "Evolved gridfunction \"hDD12\" has parity type 8.\n",
      "Evolved gridfunction \"hDD22\" has parity type 9.\n",
      "Evolved gridfunction \"lambdaU0\" has parity type 1.\n",
      "Evolved gridfunction \"lambdaU1\" has parity type 2.\n",
      "Evolved gridfunction \"lambdaU2\" has parity type 3.\n",
      "Evolved gridfunction \"trK\" has parity type 0.\n",
      "Evolved gridfunction \"vetU0\" has parity type 1.\n",
      "Evolved gridfunction \"vetU1\" has parity type 2.\n",
      "Evolved gridfunction \"vetU2\" has parity type 3.\n",
      "Auxiliary gridfunction \"H\" has parity type 0.\n",
      "Auxiliary gridfunction \"psi4_part0im\" has parity type 0.\n",
      "Auxiliary gridfunction \"psi4_part0re\" has parity type 0.\n",
      "Auxiliary gridfunction \"psi4_part1im\" has parity type 0.\n",
      "Auxiliary gridfunction \"psi4_part1re\" has parity type 0.\n",
      "Auxiliary gridfunction \"psi4_part2im\" has parity type 0.\n",
      "Auxiliary gridfunction \"psi4_part2re\" has parity type 0.\n",
      "AuxEvol gridfunction \"RbarDD00\" has parity type 4.\n",
      "AuxEvol gridfunction \"RbarDD01\" has parity type 5.\n",
      "AuxEvol gridfunction \"RbarDD02\" has parity type 6.\n",
      "AuxEvol gridfunction \"RbarDD11\" has parity type 7.\n",
      "AuxEvol gridfunction \"RbarDD12\" has parity type 8.\n",
      "AuxEvol gridfunction \"RbarDD22\" has parity type 9.\n"
     ]
    }
   ],
   "source": [
    "import CurviBoundaryConditions.CurviBoundaryConditions_new_way as CBC\n",
    "CBC.CurviBoundaryConditions_register_C_functions_and_NRPy_basic_defines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mainc'></a>\n",
    "\n",
    "# Step 5: The C code `main()` function for `BrillLindquist_Playground` \\[Back to [top](#toc)\\]\n",
    "$$\\label{mainc}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:03.207038Z",
     "iopub.status.busy": "2021-09-23T04:29:03.206710Z",
     "iopub.status.idle": "2021-09-23T04:29:03.207616Z",
     "shell.execute_reply": "2021-09-23T04:29:03.207793Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_main__BrillLindquist_Playground():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\", \"time.h\"]\n",
    "    desc = \"\"\"// 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",
    "\"\"\"\n",
    "    c_type = \"int\"\n",
    "    name = \"main\"\n",
    "    params = \"int argc, const char *argv[]\"\n",
    "    body = r\"\"\"\n",
    "  paramstruct params;\n",
    "  set_Cparameters_to_default(&params);\n",
    "\n",
    "  // Step 0a: Set free parameters, overwriting Cparameters defaults\n",
    "  //          by hand or with command-line input, as desired.\n",
    "#include \"free_parameters.h\"\n",
    "\n",
    "  // Step 0b: Read command-line input, error out if nonconformant\n",
    "  if((argc != 4 && argc != 5) || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < NGHOSTS || atoi(argv[3]) < 2 /* FIXME; allow for axisymmetric sims */) {\n",
    "    fprintf(stderr,\"Error: Expected three command-line arguments: ./BrillLindquist_Playground Nx0 Nx1 Nx2,\\n\");\n",
    "    fprintf(stderr,\"where Nx[0,1,2] is the number of grid points in the 0, 1, and 2 directions.\\n\");\n",
    "    fprintf(stderr,\"Nx[] MUST BE larger than NGHOSTS (= %d)\\n\",NGHOSTS);\n",
    "    exit(1);\n",
    "  }\n",
    "  if(argc == 5) {\n",
    "    CFL_FACTOR = strtod(argv[4],NULL);\n",
    "  }\n",
    "  if(CFL_FACTOR > 0.5 && atoi(argv[3])!=2) {\n",
    "    fprintf(stderr,\"WARNING: CFL_FACTOR was set to %e, which is > 0.5.\\n\",CFL_FACTOR);\n",
    "    fprintf(stderr,\"         This will generally only be stable if the simulation is purely axisymmetric\\n\");\n",
    "    fprintf(stderr,\"         However, Nx2 was set to %d>2, which implies a non-axisymmetric simulation\\n\",atoi(argv[3]));\n",
    "  }\n",
    "  // Step 0c: 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 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 0e: Find ghostzone mappings; set up bcstruct\n",
    "    driver_bcstruct(&params, &bcstruct, xx);\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",
    "\n",
    "  // Step 0h: Time coordinate parameters\n",
    "  const REAL t_final = 1.5*domain_size; /* Final time is set so that at t=t_final,\n",
    "                                         * data at the origin have not been corrupted\n",
    "                                         * by the approximate outer boundary condition */\n",
    "\n",
    "  // Step 0i: Set timestep based on smallest proper distance between gridpoints and CFL factor\n",
    "  REAL dt = find_timestep(&params, xx, CFL_FACTOR);\n",
    "  //fprintf(stderr,\"# Timestep set to = %e\\n\",(double)dt);\n",
    "  // Add 0.5 to account for C rounding down typecasts to integers:\n",
    "  int N_final = (int)(t_final / dt + 0.5); // The number of sample points in time.\n",
    "  const REAL out_psi4_approx_every_t = 0.5;\n",
    "  int output_psi4_every_N = (int)(out_psi4_approx_every_t*((REAL)N_final)/t_final);\n",
    "  if(output_psi4_every_N == 0) output_psi4_every_N = 1;\n",
    "\n",
    "  const REAL out_2D_approx_every_t = 5.0;\n",
    "  int output_2D_every_N = (int)(out_2D_approx_every_t*((REAL)N_final)/t_final);\n",
    "  if(output_2D_every_N == 0) output_2D_every_N = 1;\n",
    "\n",
    "  // Step 0j: Error out if the number of auxiliary gridfunctions outnumber evolved gridfunctions.\n",
    "  //              This is a limitation of the RK method. You are always welcome to declare & allocate\n",
    "  //              additional gridfunctions by hand.\n",
    "  if(NUM_AUX_GFS > NUM_EVOL_GFS) {\n",
    "    fprintf(stderr,\"Error: NUM_AUX_GFS > NUM_EVOL_GFS. Either reduce the number of auxiliary gridfunctions,\\n\");\n",
    "    fprintf(stderr,\"       or allocate (malloc) by hand storage for *diagnostic_output_gfs. \\n\");\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  // Step 0k: Declare struct for gridfunctions and allocate memory for y_n_gfs gridfunctions\n",
    "  MoL_gridfunctions_struct gridfuncs;\n",
    "  MoL_malloc_y_n_gfs(&params, &gridfuncs);\n",
    "\"\"\"\n",
    "    if enable_rfm_precompute:\n",
    "        body += \"\"\"\n",
    "  // Step 0l: Set up precomputed reference metric arrays\n",
    "  // Step 0l.i: Allocate space for precomputed reference metric arrays.\n",
    "  rfm_struct rfmstruct;\n",
    "  rfm_precompute_rfmstruct_malloc(&params, &rfmstruct);\n",
    "\n",
    "  // Step 0l.ii: Define precomputed reference metric arrays.\n",
    "  rfm_precompute_rfmstruct_define(&params, xx, &rfmstruct);\\n\"\"\"\n",
    "    body += r\"\"\"\n",
    "  // Step 1: Set up initial data to an exact solution\n",
    "  initial_data(&params, xx, gridfuncs.y_n_gfs);\n",
    "\n",
    "  // Step 1a: Allocate memory for non-y_n_gfs gridfunctions\n",
    "  MoL_malloc_non_y_n_gfs(&params, &gridfuncs);\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, gridfuncs.y_n_gfs);\n",
    "  enforce_detgammahat_constraint(&params, &rfmstruct, gridfuncs.y_n_gfs);\n",
    "\n",
    "  // Step 2: Start the timer, for keeping track of how fast the simulation is progressing.\n",
    "#ifdef __linux__ // Use high-precision timer in Linux.\n",
    "  struct timespec start, end;\n",
    "  clock_gettime(CLOCK_REALTIME, &start);\n",
    "#else     // Resort to low-resolution, standards-compliant timer in non-Linux OSs\n",
    "  // http://www.cplusplus.com/reference/ctime/time/\n",
    "  time_t start_timer,end_timer;\n",
    "  time(&start_timer); // Resolution of one second...\n",
    "#endif\n",
    "\n",
    "  int num_R_ext_idx = 0;\n",
    "  for(int R_ext_idx = (Nxx_plus_2NGHOSTS0-NGHOSTS)/4; R_ext_idx<(Nxx_plus_2NGHOSTS0-NGHOSTS)*0.9; R_ext_idx+=10) num_R_ext_idx++;\n",
    "  int *restrict list_of_R_ext_idxs = (int *restrict)malloc(sizeof(int)*num_R_ext_idx);\n",
    "  int dummy = 0;\n",
    "  for(int R_ext_idx = (Nxx_plus_2NGHOSTS0-NGHOSTS)/4; R_ext_idx<(Nxx_plus_2NGHOSTS0-NGHOSTS)*0.9; R_ext_idx+=10) { list_of_R_ext_idxs[dummy] = R_ext_idx; dummy++; }\n",
    "\n",
    "  // Step 3: Integrate the initial data forward in time using the chosen RK-like Method of\n",
    "  //         Lines timestepping algorithm, and output periodic simulation diagnostics\n",
    "  for(int n=0;n<=N_final;n++) { // Main loop to progress forward in time.\n",
    "    const REAL time = n*dt;\n",
    "\n",
    "    // Step 3.a: Check for NaNs -- has the simulation gone bad?\n",
    "    if(n % 10 == 0) {\n",
    "      int num_nans = 0;\n",
    "#pragma omp parallel for reduction(+:num_nans)\n",
    "      LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,  NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,  NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS)\n",
    "        if(isnan(gridfuncs.y_n_gfs[IDX4S(CFGF, i0, i1, i2)])) num_nans++;\n",
    "      if(num_nans > 0) {\n",
    "        char filename[100];  sprintf(filename, \"nan_checker_results.txt\"); FILE *nancheck_file = fopen(filename, \"w\");\n",
    "        LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,  NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,  NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS) {\n",
    "          if(isnan(gridfuncs.y_n_gfs[IDX4S(CFGF, i0, i1, i2)])) {\n",
    "            REAL xCart[3];  xx_to_Cart(&params, xx, i0,i1,i2, xCart);\n",
    "            fprintf(nancheck_file, \"%e %e %e | ERROR, found NAN at time %e ; %d %d %d\\n\", xCart[0], xCart[1], xCart[2], time, i0,i1,i2);\n",
    "          }\n",
    "        }\n",
    "        fclose(nancheck_file);\n",
    "        fprintf(stderr, \"FOUND NANs, exiting...\\n\"); exit(1);\n",
    "      }\n",
    "    }\n",
    "\n",
    "    // Step 3.b: Perform diagnostics.\n",
    "    if(n % output_2D_every_N == 0) {\n",
    "      // Step 3.b.i: Output 2D data file periodically, for visualization\n",
    "      // Evaluate BSSN constraints (currently only Hamiltonian constraint violation computed).\n",
    "      Ricci_eval(&params, &rfmstruct, gridfuncs.y_n_gfs, gridfuncs.auxevol_gfs);  // <- depends on Ricci.\n",
    "      BSSN_constraints(&params, &rfmstruct, gridfuncs.y_n_gfs, gridfuncs.auxevol_gfs, gridfuncs.diagnostic_output_gfs);\n",
    "\n",
    "      char filename[100];\n",
    "      sprintf(filename,\"out%d-%08d.txt\",Nxx[0],n);\n",
    "      FILE *out2D = fopen(filename, \"w\");\n",
    "      LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,\n",
    "                  NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,\n",
    "                  NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS) {\n",
    "        const int idx = IDX3S(i0,i1,i2);\n",
    "        REAL xCart[3];  xx_to_Cart(&params,xx,i0,i1,i2, xCart);\n",
    "        fprintf(out2D,\"%e %e %e %e\\n\",\n",
    "                xCart[1],xCart[2],\n",
    "                gridfuncs.y_n_gfs[IDX4ptS(CFGF,idx)],log10(fabs(gridfuncs.diagnostic_output_gfs[IDX4ptS(HGF,idx)])));\n",
    "      }\n",
    "      fclose(out2D);\n",
    "      }\n",
    "    if(n % output_psi4_every_N == 0) {\n",
    "      // Step 3.b.ii: Decompose psi4 Weyl scalar into spin-weight -2 spherical harmonics.\n",
    "      //              This scalar encodes outgoing gravitational wave content.\n",
    "      psi4_part0(&params, gridfuncs.y_n_gfs, xx, gridfuncs.diagnostic_output_gfs);\n",
    "      psi4_part1(&params, gridfuncs.y_n_gfs, xx, gridfuncs.diagnostic_output_gfs);\n",
    "      psi4_part2(&params, gridfuncs.y_n_gfs, xx, gridfuncs.diagnostic_output_gfs);\n",
    "      driver__spherlikegrids__psi4_spinweightm2_decomposition(&params, gridfuncs.diagnostic_output_gfs,\n",
    "                                                              list_of_R_ext_idxs, num_R_ext_idx, time, xx, xx_to_Cart);\n",
    "    }\n",
    "\n",
    "    // Step 3.c: Step forward one timestep (t -> t+dt) in time using\n",
    "    //           chosen RK-like MoL timestepping algorithm\n",
    "    MoL_step_forward_in_time(&params, &rfmstruct, &bcstruct, &gridfuncs, dt);\n",
    "\n",
    "    // Step 3.d: If t=t_final, output conformal factor & Hamiltonian\n",
    "    //           constraint violation to 2D data file\n",
    "    if(n==N_final-1) {\n",
    "      // Evaluate BSSN constraints (currently only Hamiltonian constraint violation computed).\n",
    "      Ricci_eval(&params, &rfmstruct, gridfuncs.y_n_gfs, gridfuncs.auxevol_gfs);  // <- depends on Ricci.\n",
    "      BSSN_constraints(&params, &rfmstruct, gridfuncs.y_n_gfs, gridfuncs.auxevol_gfs, gridfuncs.diagnostic_output_gfs);\n",
    "      char filename[100];\n",
    "      sprintf(filename,\"out%d.txt\",Nxx[0]);\n",
    "      FILE *out2D = fopen(filename, \"w\");\n",
    "      LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,\n",
    "                  NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,\n",
    "                  NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS) {\n",
    "        int idx = IDX3S(i0,i1,i2);\n",
    "        REAL xCart[3];  xx_to_Cart(&params,xx,i0,i1,i2, xCart);\n",
    "        fprintf(out2D,\"%e %e %e %e\\n\",xCart[1],xCart[2],\n",
    "                gridfuncs.y_n_gfs[IDX4ptS(CFGF,idx)],log10(fabs(gridfuncs.diagnostic_output_gfs[IDX4ptS(HGF,idx)])));\n",
    "      }\n",
    "      fclose(out2D);\n",
    "    }\n",
    "\n",
    "    // Step 3.e: Progress indicator printing to stderr.\n",
    "    //      Note that we're at the end of an iteration, so n=0\n",
    "    //      indicates we are actually at the start of iteration 1.\n",
    "    if((n+1) % 50 == 0 || n == N_final) {\n",
    "      // Step 3.e.i: Measure average time per iteration\n",
    "#ifdef __linux__ // Use high-precision timer in Linux.\n",
    "      clock_gettime(CLOCK_REALTIME, &end);\n",
    "      const long long unsigned int time_in_ns = 1000000000L * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;\n",
    "#else     // Resort to low-resolution, standards-compliant timer in non-Linux OSs\n",
    "      time(&end_timer); // Resolution of one second...\n",
    "      REAL time_in_ns = difftime(end_timer,start_timer)*1.0e9+0.5; // Round up to avoid divide-by-zero.\n",
    "#endif\n",
    "      const REAL s_per_iteration_avg = ((REAL)time_in_ns / (REAL)(n+1)) / 1.0e9;\n",
    "      //  We are at the n+1st iteration now  ------v        -------^\n",
    "      const int iterations_remaining = N_final - (n+1);\n",
    "      const REAL seconds_remaining = (int)(s_per_iteration_avg * (REAL)iterations_remaining);\n",
    "      const int time_remaining__hrs  = (int)(seconds_remaining / 3600.0);\n",
    "      const int time_remaining__mins = (int)(seconds_remaining / 60.0) % 60;\n",
    "      const int time_remaining__secs = (int)(seconds_remaining) - time_remaining__mins*60 - time_remaining__hrs*3600;\n",
    "      //const REAL num_RHS_pt_evals = (REAL)(Nxx[0]*Nxx[1]*Nxx[2]) * 4.0 * (REAL)(n+1); // 4 RHS evals per gridpoint for RK4\n",
    "      //const REAL RHS_pt_evals_per_sec = num_RHS_pt_evals / ((REAL)time_in_ns / 1.0e9);\n",
    "\n",
    "      // Step 3.e.ii: Output simulation progress to stderr\n",
    "      fprintf(stderr,\"%c[2K\", 27); // Clear the line\n",
    "      fprintf(stderr,\"It: %d t=%.2f dt=%.2e | log10H: %.1f | %.1f%%; ETA %dh%dm%ds | t/h %.2f\\r\",\n",
    "              (n+1), (double)((n+1)*dt), (double)dt,\n",
    "              (double)log10(fabs(gridfuncs.diagnostic_output_gfs[IDX4ptS(HGF,(int)(Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2*0.5))])),\n",
    "              (double)(100.0*(REAL)(n+1)/(REAL)N_final), time_remaining__hrs, time_remaining__mins, time_remaining__secs,\n",
    "              (double)(dt * 3600.0 / s_per_iteration_avg));\n",
    "      fflush(stderr); // Flush the stderr buffer\n",
    "    } // End progress indicator if(n % 50 == 0)\n",
    "  } // End main loop to progress forward in time.\n",
    "  fprintf(stderr,\"\\n\"); // Clear the final line of output from progress indicator.\n",
    "\n",
    "  // Step 4: Free all allocated memory\n",
    "  free(list_of_R_ext_idxs);\n",
    "\"\"\"\n",
    "    if enable_rfm_precompute:\n",
    "        body += \"\"\"  rfm_precompute_rfmstruct_freemem(&params, &rfmstruct);\"\"\"\n",
    "    body += r\"\"\"\n",
    "  freemem_bcstruct(&params, &bcstruct);\n",
    "  MoL_free_memory_y_n_gfs(&params, &gridfuncs);\n",
    "  MoL_free_memory_non_y_n_gfs(&params, &gridfuncs);\n",
    "  for(int i=0;i<3;i++) free(xx[i]);\n",
    "\n",
    "  return 0;\n",
    "\"\"\"\n",
    "    # As rfmstruct stores functions of xx, when rfm_precompute is disabled,\n",
    "    #   we always pass xx to a function instead of &rfmstruct.\n",
    "    if not enable_rfm_precompute:\n",
    "        body = body.replace(\"&rfmstruct\", \"xx\")\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"), enableCparameters=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='compileexec'></a>\n",
    "\n",
    "# Step 6: Compile generated C codes & perform the black hole collision calculation \\[Back to [top](#toc)\\]\n",
    "$$\\label{compileexec}$$\n",
    "\n",
    "First we register remaining C functions and contributions to `NRPy_basic_defines.h`, then we output `NRPy_basic_defines.h` and `NRPy_function_prototypes.h`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:03.210499Z",
     "iopub.status.busy": "2021-09-23T04:29:03.210028Z",
     "iopub.status.idle": "2021-09-23T04:29:03.212175Z",
     "shell.execute_reply": "2021-09-23T04:29:03.211754Z"
    }
   },
   "outputs": [],
   "source": [
    "add_to_Cfunction_dict_main__BrillLindquist_Playground()\n",
    "\n",
    "import outputC as outC\n",
    "outC.outputC_register_C_functions_and_NRPy_basic_defines()  # #define M_PI,  etc.\n",
    "# Declare paramstruct, register set_Cparameters_to_default(),\n",
    "#   and output declare_Cparameters_struct.h and set_Cparameters[].h:\n",
    "outC.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(os.path.join(Ccodesrootdir))\n",
    "\n",
    "gri.register_C_functions_and_NRPy_basic_defines()  # #define IDX3S(),  etc.\n",
    "fin.register_C_functions_and_NRPy_basic_defines(NGHOSTS_account_for_onezone_upwind=True,\n",
    "                                                enable_SIMD=enable_SIMD)  # #define NGHOSTS, and UPWIND() macro if SIMD disabled\n",
    "\n",
    "# Output functions for computing all finite-difference stencils.\n",
    "#   Must be called after defining all functions depending on FD stencils.\n",
    "if enable_FD_functions:\n",
    "    fin.output_finite_difference_functions_h(path=Ccodesrootdir)\n",
    "\n",
    "# Call this last: Set up NRPy_basic_defines.h and NRPy_function_prototypes.h.\n",
    "outC.construct_NRPy_basic_defines_h(Ccodesrootdir, enable_SIMD=enable_SIMD)\n",
    "outC.construct_NRPy_function_prototypes_h(Ccodesrootdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we output all the C codes in `outC_function_dict` to files in the directory `Ccodesrootdir`, generate a `Makefile`, and compile the project using a parallel `make` command. If the `make` command fails, a backup serial compilation script is run.\n",
    "\n",
    "To aid in the cross-platform-compatible (with Windows, MacOS, & Linux) compilation and execution, we make use of `cmdline_helper` [(**Tutorial**)](Tutorial-cmdline_helper.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:03.215586Z",
     "iopub.status.busy": "2021-09-23T04:29:03.215113Z",
     "iopub.status.idle": "2021-09-23T04:29:09.658960Z",
     "shell.execute_reply": "2021-09-23T04:29:09.658568Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(EXEC): Executing `make -j34`...\n",
      "(BENCH): Finished executing in 5.210599660873413 seconds.\n",
      "Finished compilation.\n",
      "(EXEC): Executing `taskset -c 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 ./BrillLindquist_Playground 270 8 2 1.0`...\n",
      "\u001b[2KIt: 30612 t=225.01 dt=7.35e-03 | log10H: -2.3 | 100.0%; ETA 0h0m0s | t/h 5427.03\n",
      "(BENCH): Finished executing in 149.44121980667114 seconds.\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd\n",
    "cmd.new_C_compile(Ccodesrootdir, os.path.join(\"output\", \"BrillLindquist_Playground\"),\n",
    "                  uses_free_parameters_h=True, compiler_opt_option=\"fast\") # fastdebug or debug also supported\n",
    "\n",
    "# Change to output directory\n",
    "os.chdir(os.path.join(Ccodesrootdir, \"output\"))\n",
    "# Clean up existing output files\n",
    "cmd.delete_existing_files(\"out*.txt\")\n",
    "cmd.delete_existing_files(\"out*.png\")\n",
    "# Run executable with CFL_FACTOR = 1.0, which is allowed since\n",
    "#   simulation is axisymmetric and all phi derivs are set to zero.\n",
    "# cmd.Execute(\"BrillLindquist_Playground\", \"72 12 2 \"+str(CFL_FACTOR))\n",
    "# cmd.Execute(\"BrillLindquist_Playground\", \"800 16 2 \"+str(CFL_FACTOR))\n",
    "cmd.Execute(\"BrillLindquist_Playground\", str(Nr) + \" \" + str(Nth) + \" \" + str(Nph) + \" \" + str(CFL_FACTOR))\n",
    "os.chdir(os.path.join(\"..\", \"..\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='visualize'></a>\n",
    "\n",
    "# Step 7: Visualize the output! \\[Back to [top](#toc)\\]\n",
    "$$\\label{visualize}$$ \n",
    "\n",
    "In this section we will generate a movie, plotting the conformal factor of these initial data on a 2D grid, such that darker colors imply stronger gravitational fields. Hence, we see the two black holes each of mass $0.5$ initially centered at $z/M=\\pm 0.25$, where $M$ is a mass scale (conventionally the [ADM mass](https://en.wikipedia.org/w/index.php?title=ADM_formalism&oldid=846335453) is chosen), and our formulation of Einstein's equations adopt $G=c=1$ [geometrized units](https://en.wikipedia.org/w/index.php?title=Geometrized_unit_system&oldid=861682626)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='installdownload'></a>\n",
    "\n",
    "## Step 7.a: Install `scipy` and download `ffmpeg` if they are not yet installed/downloaded \\[Back to [top](#toc)\\]\n",
    "$$\\label{installdownload}$$ \n",
    "\n",
    "Note that if you are not running this within `mybinder`, but on a Windows system, `ffmpeg` must be installed using a separate package (on [this site](http://ffmpeg.org/)), or if running Jupyter within Anaconda, use the command: `conda install -c conda-forge ffmpeg`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:09.666760Z",
     "iopub.status.busy": "2021-09-23T04:29:09.666506Z",
     "iopub.status.idle": "2021-09-23T04:29:10.197073Z",
     "shell.execute_reply": "2021-09-23T04:29:10.197468Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!pip install scipy > /dev/null\n",
    "\n",
    "check_for_ffmpeg = !which ffmpeg >/dev/null && echo $?\n",
    "if check_for_ffmpeg != ['0']:\n",
    "    print(\"Couldn't find ffmpeg, so I'll download it.\")\n",
    "    # Courtesy https://johnvansickle.com/ffmpeg/\n",
    "    !wget http://astro.phys.wvu.edu/zetienne/ffmpeg-static-amd64-johnvansickle.tar.xz\n",
    "    !tar Jxf ffmpeg-static-amd64-johnvansickle.tar.xz\n",
    "    print(\"Copying ffmpeg to ~/.local/bin/. Assumes ~/.local/bin is in the PATH.\")\n",
    "    !mkdir ~/.local/bin/\n",
    "    !cp ffmpeg-static-amd64-johnvansickle/ffmpeg ~/.local/bin/\n",
    "    print(\"If this doesn't work, then install ffmpeg yourself. It should work fine on mybinder.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='genimages'></a>\n",
    "\n",
    "## Step 7.b: Generate images for visualization animation \\[Back to [top](#toc)\\]\n",
    "$$\\label{genimages}$$ \n",
    "\n",
    "Here we loop through the data files output by the executable compiled and run in [the previous step](#mainc), generating a [png](https://en.wikipedia.org/wiki/Portable_Network_Graphics) image for each data file.\n",
    "\n",
    "**Special thanks to Terrence Pierre Jacques. His work with the first versions of these scripts greatly contributed to the scripts as they exist below.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:10.205155Z",
     "iopub.status.busy": "2021-09-23T04:29:10.204679Z",
     "iopub.status.idle": "2021-09-23T04:29:12.375868Z",
     "shell.execute_reply": "2021-09-23T04:29:12.376113Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[2KProcessing file BSSN_Two_BHs_Collide_Ccodes_psi4_new_way/output/out270-00030600.txt\r"
     ]
    }
   ],
   "source": [
    "## VISUALIZATION ANIMATION, PART 1: Generate PNGs, one per frame of movie ##\n",
    "\n",
    "import numpy as np\n",
    "from scipy.interpolate import griddata\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.pyplot import savefig\n",
    "from IPython.display import HTML\n",
    "import matplotlib.image as mgimg\n",
    "\n",
    "import glob\n",
    "import sys\n",
    "from matplotlib import animation\n",
    "\n",
    "globby = glob.glob(os.path.join(outdir,'out*-00*.txt'))\n",
    "file_list = []\n",
    "for x in sorted(globby):\n",
    "    file_list.append(x)\n",
    "\n",
    "bound=1.4\n",
    "pl_xmin = -bound\n",
    "pl_xmax = +bound\n",
    "pl_ymin = -bound\n",
    "pl_ymax = +bound\n",
    "\n",
    "for filename in file_list:\n",
    "    fig = plt.figure()\n",
    "    x,y,cf,Ham = np.loadtxt(filename).T #Transposed for easier unpacking\n",
    "\n",
    "    plotquantity = cf\n",
    "    plotdescription = \"Numerical Soln.\"\n",
    "    plt.title(\"Black Hole Head-on Collision (conf factor)\")\n",
    "    plt.xlabel(\"y/M\")\n",
    "    plt.ylabel(\"z/M\")\n",
    "\n",
    "    grid_x, grid_y = np.mgrid[pl_xmin:pl_xmax:300j, pl_ymin:pl_ymax:300j]\n",
    "    points = np.zeros((len(x), 2))\n",
    "    for i in range(len(x)):\n",
    "        # Zach says: No idea why x and y get flipped...\n",
    "        points[i][0] = y[i]\n",
    "        points[i][1] = x[i]\n",
    "\n",
    "    grid = griddata(points, plotquantity, (grid_x, grid_y), method='nearest')\n",
    "    gridcub = griddata(points, plotquantity, (grid_x, grid_y), method='cubic')\n",
    "    im = plt.imshow(gridcub, extent=(pl_xmin,pl_xmax, pl_ymin,pl_ymax))\n",
    "    ax = plt.colorbar()\n",
    "    ax.set_label(plotdescription)\n",
    "    savefig(os.path.join(filename+\".png\"),dpi=150)\n",
    "    plt.close(fig)\n",
    "    sys.stdout.write(\"%c[2K\" % 27)\n",
    "    sys.stdout.write(\"Processing file \"+filename+\"\\r\")\n",
    "    sys.stdout.flush()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='genvideo'></a>\n",
    "\n",
    "## Step 7.c: Generate visualization animation \\[Back to [top](#toc)\\]\n",
    "$$\\label{genvideo}$$ \n",
    "\n",
    "In the following step, [ffmpeg](http://ffmpeg.org) is used to generate an [mp4](https://en.wikipedia.org/wiki/MPEG-4) video file, which can be played directly from this Jupyter notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:12.872428Z",
     "iopub.status.busy": "2021-09-23T04:29:12.869612Z",
     "iopub.status.idle": "2021-09-23T04:29:12.874425Z",
     "shell.execute_reply": "2021-09-23T04:29:12.874102Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"480\" height=\"360\" controls>\n",
       "  <source src=\"BSSN_Two_BHs_Collide_Ccodes_psi4_new_way/output/BH_Head-on_Collision.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## VISUALIZATION ANIMATION, PART 2: Combine PNGs to generate movie ##\n",
    "\n",
    "# https://stackoverflow.com/questions/14908576/how-to-remove-frame-from-matplotlib-pyplot-figure-vs-matplotlib-figure-frame\n",
    "# https://stackoverflow.com/questions/23176161/animating-pngs-in-matplotlib-using-artistanimation\n",
    "\n",
    "fig = plt.figure(frameon=False)\n",
    "ax = fig.add_axes([0, 0, 1, 1])\n",
    "ax.axis('off')\n",
    "\n",
    "myimages = []\n",
    "\n",
    "for i in range(len(file_list)):\n",
    "    img = mgimg.imread(file_list[i]+\".png\")\n",
    "    imgplot = plt.imshow(img)\n",
    "    myimages.append([imgplot])\n",
    "\n",
    "ani = animation.ArtistAnimation(fig, myimages, interval=100,  repeat_delay=1000)\n",
    "plt.close()\n",
    "ani.save(os.path.join(outdir,'BH_Head-on_Collision.mp4'), fps=5,dpi=150)\n",
    "\n",
    "## VISUALIZATION ANIMATION, PART 3: Display movie as embedded HTML5 (see next cell) ##\n",
    "\n",
    "# https://stackoverflow.com/questions/18019477/how-can-i-play-a-local-video-in-my-ipython-notebook\n",
    "\n",
    "# Embed video based on suggestion:\n",
    "#  https://stackoverflow.com/questions/39900173/jupyter-notebook-html-cell-magic-with-python-variable\n",
    "HTML(\"\"\"\n",
    "<video width=\"480\" height=\"360\" controls>\n",
    "  <source src=\\\"\"\"\"+os.path.join(outdir,\"BH_Head-on_Collision.mp4\")+\"\"\"\\\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='numerror'></a>\n",
    "\n",
    "# Step 8: Plot the numerical error at the end of the simulation \\[Back to [top](#toc)\\]\n",
    "$$\\label{numerror}$$\n",
    "\n",
    "Here we plot the log10-absolute value Hamiltonian constraint violation on the $x$-$z$ plane, near the black hole (at $x=y=z=0$). Notice the constraint violation is largest inside the horizon, as expected. Notice also the black hole is a little bit less than $1M$ in radius."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:12.886917Z",
     "iopub.status.busy": "2021-09-23T04:29:12.877392Z",
     "iopub.status.idle": "2021-09-23T04:29:13.021121Z",
     "shell.execute_reply": "2021-09-23T04:29:13.020874Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAEWCAYAAAD7MitWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABUEklEQVR4nO29ebR9WVXf+5lrn3Pv/TXVQYFlaAQFecHeV2KMz2dXMejQ8OTZgIkGNVZssBmagWKFJxH12URjEk2wGCI6RJoYSzSKNCaIJuqjNARBUFFBClEommp+v9+99+y95vtjrb3PPvvu/uzunLu/Y5xx7j1n7bXmbtb3zG7NJarKjBkzZpwHmLEFmDFjxoyhMBPejBkzzg1mwpsxY8a5wUx4M2bMODeYCW/GjBnnBjPhzZgx49xgJrwZvUJEVEQeN8A4t4jI74jIfSLy5R31+SwR+aCIvEpELuZ8/zoR+awuxmoDEXmuiDx3rPF3ETPh5UBEvkxE/oeIXBWR1+V8H4jI94nIX4vIAyLyP0XkxoK+XuQn/ZNSnz1OREZJgPST9FhEHky9fnUMWTrGlwJL4GGq+rL4QxF5noj8kYiEeeQgIl8hIu8UkSsi8ssi8pD4O1X9YeCRwP8GfF7Z4CLyGH+fF5nPXyQi37fdqc3oCjPh5eMDwI8DP1jw/b8C/j7wacD1wFcCxxX9Temhf6aqXk69viivUXbyFn02ETwE+DNVPc18/nbgWcCvZQ8QkY8Bfgp3/z4MuAr8h3QbVb0C/CXw0B5knjEwziXhichHicgHROST/f9/R0TeF5snqvpaVX058Nc5x94EfBvwdar6TnV4s6qWEd7PAh8vIp9ZIM87ROS21P/PFZGf93/HmsNXi8i7vIn19SLyKSLyJhH5kIj8RLsrcUaOzxKRe0TkO0Xkb4Cf8bL8ooj8vIjcDzxji/5vEJGf89f6nSLyL0XE+O8CEflREblXRP5SRJ6ZpzGVYAHY7Ieq+rOq+krggZxj/jHwq6r6elV9EHgO8FQRuS7Tzvr+t4aI/CcR+Rtver/ek2783YtE5D+IyCu95v3fvan+4/6+v01EPqkLOc4rziXhqeqfA98J/Lz3zfwM8LOq+roah38cEAJf4h/cPxWRb6o45irwA8D3byH2pwKPB74cp33eAdwGfAzwZUVk2gK34LSljwBu9589BfhF4EbgxenGfoJuaEUl+PfADcBHAp8JfBXw1f67rwM+H/hE4JOB/6uuwP4efjrwV3WP8fgY4H/F//jn4hT46Ey7dwGfJSIHDfvPwytx9/HhwB+SuZ7AlwH/ErgZOAF+17e7GXcPfqwDGc4tziXhAajqC3Dmzu8DH44jkDp4JG7SfjTwWOBLgOeKyD+oOO6ngEeLyOe3k5jnqeqxqr4auAK8RFXfq6rvBn4baPLL/++8Zhi/npf6zgLfo6onqnrNf/a7qvrLqmpTnwGgqt+oqt9YNaCIBMDTgGer6gOq+g7gR3HmJLiJ/m9V9R5V/SDF7oRsv/8Ip7092vfXBJeB+zKf3QdkNbzvBT4buCIit1b0eW/62gJfkf5SVV/oz/8EeC7wCSJyQ6rJXar6B95iuAs4VtWfU9UIeBnN7vOMDM4t4Xm8APhY4N/7B7AO4gn/vap6TVXfBLwU+IKyg3z/z/OvNvjbjAzZ/y836OtbVPXG1Os5qe/el2Oev6uhrHm4GRdUeGfqs3cCj/B//53MOLXGVNVfwWmk7we+tqFMD+J8sGlcz1nz91uBNwDXq+rdFX3enL62wC/EX3iz/QdF5M+9e+Ad8TGp47u8zzMyOLeEJyKXcabhT+M0tIeUH5HgTf49HWWtG3H9GZxZ+NTM51eAdNrDLTX76wN559JFRPleYIUzlWM8Gni3//s9OO05xqPqdqyq9wH/FXhiQ5neAnxC/I+IfCRwCPxppt3fBX4jq922wFfg3AO34ayEx8RDb9nvjJo4t4QH/FvgblX9Z7gI3vPjL/wv8RHOUW1E5EhElpD4eX4buENEDkXk7+JMtf9SNaCqhsD34PyHabwReJqILL3J9CVbn93E4E2ylwPfLyLXichHAN8O/Lxv8nLgW0XkEeJSfLLXqAonwBkfm7+mR7hnfeHvZeC/fjHwRSLyGSJyCWe6/pKqZjW8pe9/W1zn+3k/7gfuBzroc0YDnEvCE5GnAE8GvsF/9O3AJ4vIP/b/fyXOfPiPwGf4v1+Q6uLpOE3l/TiyfI6q/mbN4V+C02bSeA7wUcAHcSkvv5A9qC785H2wotlPyGYe3h806P+7ReSVqf+fLyLPLzsmhW/GabN/AfwO7jxf6L97AfBqnAb9P4FfxwWHopp9W/Kf5xfg7t/TcX7aa3i/oaq+Bfh6HPG9F0dIef7IgJwIcAv8HM6Mfzfwx8DvddDnjAaQuQDojCnCB3eer6ofUdnYtb8d+Cbg7/vcua7kuBFHTv9MVX89893rgOfWjO53DvGJ1Kr63DHG30WcSw1vxvQgIhdE5AtEZCEij8CZ/nc16OIXcYGOvxSRL+tIpn8BvA34LaCuBj9jwpg1vBmTgM+l+y3cMq5rOFfBt6rq/SUm+uer6m8PJWMWIvIM4HU+xWaM8T8LYCwNcxcxE96MGTPODWaTdsaMGecGU10Inovl4SU9vFQ3XW5gdJBJtTO6dl9ZYztzATq4BBM915MrH2B1cmWr0/uHn31J3/+BesH1P3jTyatU9cnbjNcEO0V4h5cewsd93rdt3Y/2MWG7ILymfVS0r9VfwzFV+s+RlSZulhpNaxfiajDs1sW9eiK8beX6o1f/+NYy3PuBiN9/1SOrGwLLD//zm6tbdYedIrxtsTdEV2Pcyj4bjDkEyRWNV0l+cdOSZio1iUDK+8n2CVsQTA2522BrubqRgki7SFvsHrMPb2SMQnY1oSKDk12eDN3000k3M2pAAYvWetWBiDxZRP5ERN4uIt+1jWznRsPr/IHfU81ubILLQy2Nr6am5/ooGayh5lVbe+xovLoYW9OznSxMSars/CTwD4B7gDeIyK+o6h+36W/vCW+qv+xd++sq++yT6Noc2nIixnKWEl9F31uTVE5/Tqbu+txlKMqqO5P2ScDbVfUvAETkpbgCDDPhDYIJBie2Pb422XX541HVVyVpSb+k15PmNfR4XZN7rTGBqP6J3Cwi6ZJbd6rqnan/H8FmqbB7cMVwW2GvCW+K2l1fZFfY7zZk1wO515586eMKjumd9BpiqpreGKRX1z8H3KuqVUVVO8PeEt5UI7J9oA3ZVWp1Lc+16rq3IoUS8io1cRtEXduM3Qt6Gm9I0lMg6m4F17vZrI34SNY1FBtjL6O0UyW7RnJJ9ZgqzcmuMvJaY9zs+OlXXTRtn8jVRkMuJf6afQyJnsYb0uKxNV818Abg8SLyWL+nyNOAX2kr195peFMlu0HRVt5tAyMtkO6vkbmbaVtq3hYcM2nssKanaBMfXnlfqqGIPBN4Fa4u4Qt9HcNW2CvCOy9k1/Y8CzW7gfL6qtDI31dAeu64ZpOtFgk0TEruhFR2lPRUYdWlb9TVIfz1yoY1sFeE1zk6nOiNzbe2/ZSYsk3HqjVeFbaYtLUiqR1NrK5JrzPsmmYKgBBNUVNgTwhv6prdYL6TpmRXgtoyd+UDKwxK+C6aLK8t0vQGJI9Otage5O4zoqyAnShJ7wXhdY4xf5x69L/VPaZTsquLikndhkBy/XqlEd+JanljjtsSs4bXE6a4ZCyNrkzZxn0lx9Q3ZbssONAKcf9NCKnimDExRv5bU/Qho0s8ngmvc0yd7LpE32bxNsvShkQnE3THtKUEPZF75wnYwEqnmfG2s4Q3xVUUWYztu6vbdnJklx4zZyLm+p+aklhB+z5WX0xdy4Nu5VSEaKIpvjtJeOdJs9sGdYIVXZFd23uyjc/szCRtmp+3q6Q3YTM+hp2oRjIaDYvIo0Tkv4nIH4vIW0TkW8eSpQ8MXiDgzPg5HXbsT2yzwqKon21kad12woQxBjqro4jz4dV5DY0xNbwQ+A5V/UMRuQ74AxF5TVWdq05/OHq63mP/uJ0hu47M2D7PqzJNook2lmlbuQqjQJ6msrTury0mq+kJ0ezD24Sqvgd4j//7ARF5K64UTKs6V1NBH0U9+4BoM1mHXG3RJdGs+21Ger3kp+2IP29buIrHM+EVQkQeA3wS8Ps5390O3A5wcPGmjgbsppssRtHsUmM2NWNzSW/IdJUOk43PoA45FrTpev+LNM6DpqcqnGowthi5GJ2GReQy8J+Bb1PV+7Pfq+qdqnqrqt66PLw0vIC7ji0nQa+5eW3yDtsQckOM7ZLYB1ik1mtojKrhicgSR3YvVtVf6n/A/roefJK0cc5njtmQeax0lbbJxq20q2YrL/pCb6btRDQ9F7QYXZfKxWiEJyIC/DTwVlX9sbHk6AKDJwV3MF7dQp256ON8m6SfDIDeK4rstT9vukGLMaX6dOArgc8RkTf61xf0MpKws367Jv2LavWuXnmfl2l+FW07RUnfuaQvm9+3vhcjRKV77b/v+1SBOGhR5zU0xozS/g5D3JZ9MmM7wqTlHsHEHBO9anojXstoog/ZJKK0vWFfya7AJ1eIjCaU93mTPupgEhva7Ah57hvpKcJKp0kt05SqC+yoCdtoLKX6PKWgr57ILjtOoyrG6fG6Khq6Q6QHPQYzBrwGUw5aTFOqbTFNbboxaq0zzWAjF6/sOiiFKxeaoOkSs2z7Un/hlrLlj5/X8fb9Th5D/lAjRFrvtQ1E5Ev9slQrIrW2etw/wutRsxujsrJo/V/92mS3cVDZ4Pn9dLGGNq+/OuNXjZl7rfoMDnSM3p4zGJT0BgpavBl4KvD6ugfsl0m7J5pdJ8iYuzERlJm2Y/uZe/NlNamgksKY/sddztVTZZC0FFV9K4A02MJgPwhvV4MT2/S9QVRFeRVnx0iWk1VpTtuYlOl2PSyy32jTxQTO8XFts2VkV9jVYIYLWtReWnaziNyd+v9OVb2zB7GAfSC8XSW7DiGqtUmv9jnVCYiUHTsV5K6VnZKA5dhV0msQtLhXVQv9byLyWuCWnK/uUNVXNJVrtwlvDyKxVai7njQ9ifPIr47vq4s2Q6CMAMqinU026m4UNe05CrprEVxFOisAqqq3ddKRx+4S3gQm3k6gw+u0VVHTqStVU5evT/RAenNaSheQ1KsH9Bohy2KbcTT1yn7VJC2lDHWio1LxqtO2DhpeqzN5eB6FKSldTPYBnpveI7gdzS0FrJpar61EFvliEbkH+DTg10TkVVXH7K6G1zEGNdW6HGtbX1vbY7s6hxztYr8X1m+P6V+fYcq3q+pdwF1NjjnXhDd2wc5azav8ZnWitRSkpdQMZvRRreVMfw1I78x3qb+rI7z1Kx/XJpV0xHggIkrfk6mRn0KTKO2gOLeEtwtklxxWRHpV/W2jwW2L7Lgtqw9Xfr/tZC/otzWJDEh6Maam8anK1uZqXzh3hDeFKOPOow6RNp2AfRHMEESQzQUcILk3i0kUa0hhqvXwzg3hjUp0A42dm5pSRk6p7za0yDrydqk99qSl5TYdMgdvROKD8cjP1cObpmax14Q3CW1uBBnqbMBd3gHDyj2QGdg72RVdsxHMXBhT65tuxeO9I7zRSC6bczZE/lvTunhtxqDAh9jzdT4zSfuatHUDI9uiqU+zQwxNfC4tZQraxlnsHeGNgh7JYKjnpiwa3Ijs6prDecfVNktr9FcHTdbsdo0RkrKHCm40XEs7KPaC8HY2h25iqLV0rE20uKp9TlpH0eSsIrupRSxrYUDyG8rHN2/EvevYY6IbFR1Pup0juyxGCHR0DVceapoTZqcJbx8KJRah03PLRGOL+s+L0vb63G47qasSjjsgjdEc/3nXvetF/j2e2+zD6xD7THQ7iTZBmj4JZIe1o1L0pP11TXyuWsps0m4NpWOymyDBFZ1fnYrF67bq29Y/wbJxNetjqqN9tIlKNp1wNQMdTdJRspM+PvdJJfb2FPFV6aYrt7RsJrzpYIJEBxPJG8xBbdIj06YJxiCSFvlxkwyKTM7vN2t442OiZALVOXClbVrk4uX1tTFOTp+1i392QIRd5uAVancl+XdpGfKu/6S0vTRGSHUpwrzSYiyMeN3r7s/QzWD+vcYysjJsEFvRMdukphRpI0XaVtcTt8P+JqntxRhR65ujtENjZJLrok2lZtcURdpaDW2uVW5eHXnqrHJItSsqCVWGJuWgsjKcGbOkXVXb0chxJK1vNmmHwkhkV3uVwkRwhlDztLmm/zcWgs1JmJIlSxBlZDEEkdT5ASoiyWzgI6/NIGjhs2yDLve06Br7QXg9X9sx7l2rMcuOaeDrK/UX5mkMba9PA19d25UXZw+o0SYzbuFYNY7Zps2ubNqThQLhABqeiPwI8EXAKfDnwFer6ofKjpmm3lkXQm9kF+8fMDTZ9f7Lr6lXery0FlJ2XYv8bE3kbto+R0OClOyZc9oX9PYM9jhvYgyxpwXwGuBjVfXjgT8Fnl11wO4SXoc3LPtQjamNtx67IB8v/cprnztetqnW9J9tk2uXM2bVMWX1+wrPO27f4jpP5bnojASrftzaQp1JW+e11TCqr1bV0P/7e8Ajq47ZTZO2gxtUurRq11BT7uyG3UXpKRupKan+a0Vwi8zcMjIs8tW10NgaBSgG8ml1hez9KtR426DDa6E0Sku5WUTuTv1/p6re2WLYrwFeVtVoNwmvJXaW0MrQwzm1SjSukqdpmkSZxpduVif6WSb/jpFeGaaUH9hAe7tXVW8t+lJEXgvckvPVHar6Ct/mDiAEXlw12O4RXsMJvpckB/U0qA0zSM58lpue4ttUlnwv+44SZ39R0KNgklZpfGd8kP6PM5pekQwV4585ZgJkUoXW0eCO5orSXfEAVb2t7HsReQbwhcDnqlar96MSnoi8ECfse1X1Y8eUpVd0nQvVxFys6iqTApKYsy2f1yJXwdZaR8nKh1Z9tSG9uF1KntZt9hiKENpBorRPBp4FfKaqXq1zzNhBixcBT+6yw3Rkq9eqKk1e2WPr9FU2dhpFml2edhe3L4jUVl2vwgBHxgeX98odv2yMsqBKapy883HHpdUc6l+3uii6t5Lzf5vnpCOMkWkAzodX57UlfgK4DniNiLxRRJ5fdcCoGp6qvl5EHtNZf13e2D4fkjZaVEfyZE09RapJOQ9Nc/JqBi3OyFER1c3dqa0uijS9Mpmy6NIHWBYQ2nKMQVd66DD18FT1cU2P2T0fXg76WH7VO7YNAsDZSVCoDeU/6bUJwk/AQr9enW6KtMP4u2xeYF6eYNr8Tv5fLx2rjNDmkVnRD08T4uv6uSmTB1qT31BBjS59eF1j8oQnIrcDtwMcXLzpzPc7SXZ1UCVPTbIrPLyI7DIO72waSuP9a8uQIZVCsku31w60laxWVhXFnRq21CqHWOI2E15L+JycOwEuPfRRG7dn6yTdLZzzo6IN2VVpMngyq9Fn7Dtrpe0VRF1rpZZk2+VoebXRpSk6BvKe4Rbn04epqwjRAEGLNpg84eVhqx+PpqbYyIg1scIVA7U6Sb1XHJMEMNL9F/iT0uZm7YTk9DhFxB0HLFLtstpfbV9hnQBQLEtOoKf3zbu3RdHzPLLYU62HNyoNi8hLgN8FniAi94jI11Yds5VWN8A9UJHktU2bSjQlu6JjciZGmrxEQWxOP/6lxr3wLw3Wn228gnWb9HFqOHtvsiRHzn0vivS2uaQlx1Tdo8Hud1M0GKrrSK76oEXfS8vaYOwo7dN767zja9nmYa1zTNYcqzVOjiYCOdpIDtnlai7eLEw0NW/ubWhhnqRgk6T0zLsk7bMyiK5t5rTmFosSE6vYVBtSfcfd6tlzjeUXq5vmao7GmXvNapJe2b2qe7/roBPNsqHG1y3pTVPD20mTdmj0/csck17TcWq3b6ANFrm2REGiTb+dGt9+ww/oz+OMCXuWjNN5dJIiukLZirS6VJvkB6Smj66p/6/tvRodg/os53p4w6HFdZ7Cw1smQ2MNMOPnSmt1iQYUfx4rZJGuV134Y9XIWptSML4uRWy22oU4bU/ZcI6IpOXNnguJFhe/TJh2HIIGsjazJG6na41QXJv4mHRVFBVZy10ww88WUTirwZWR2hDPS50xphyomTW8GcOh6sH2Zp6ooqQnfiY4oLrxv2zy0uaQmejrRs5Xxrw9WxFYPGn5D6xCloQy5m1e4nGXmMKP4K5CFSI7zeu3H4TX8toO/lDn+LVy2+QQQmXXOfXucv11sX8u9pP5waKl09aClWJCWBxbgmvWaUvGaXvRoWADITwUF4RIBDxrjkrqu7wvNA5cBKAHIJGwOFZMBItjRaxFInd8dChEB4JdCNEBmAiC06wPjzMkXhZpLdLgSq911Q9JnfvbIVpFkgeScapR2t0nvGle103UlTHlH2tKimnkTtq8dI68PmJzM4Tg1BFedGA2ZNPAvc741TKR1TNjpF17cSDCgA0E4x2IcfECiUBCZ4LrwmBt7GuT3H5z4fuqRQx1yayOWTiov6wlepTRGQfTnJi7TXhT0Ozyusojl6JjcqKItcdJd1N1ThnNTgWnuYkgVrnw/ojgOEoIK7wQcO1hB0RLCC+sfX1norJ4913kh4k1R5vxO8Zk50kOA9ZriaKKBnB62Q1ycqNrvLjmNLnlNcvRB8Okv+jAcHKDQY0k4yTJ0OlEq6aact2oZlpLqnOPy9AB6bTOGeyN9OagRfeoeT07I7fsw1HWbZMh+3guUhMy7cxPfGEJaTkzVhSWD4Ys7juBhUEDIbwQcHpJiA4hvOQnVOijqeGmSag+EGLsesw0sYqqIydSyqCwce4qgj10n0UH/rMA7Aksr8LiSohEFjmJCG845OTGg+Q8AIhN6jwS6suM6+Le5cnWkojynvVa64t7IL2p5mvvJuEN+eMhBX8POW6bYxInv6yjs7D2vSkEoXJwf4RESngxILx0kejQYBfeT7dwfYqNI6MkWqJawYTqTM9Yc0yNiTnrUwSviVn3rpEntYVsRHnNWplDDRzfGBBeOMKslMWxBQsX/2aFLgwnNwTYxeZ5xZHohACyGnedyTi0gjIUSZeN3eFYs0nbFRpcx9ba3TTv1SbSv8xpc5NUAGFD01v7y0RdKohZKUfvO4bQ8uBHXubkekN4QbAHzodmQqeZYUG8Cap4/xqOmIJT3UhvyabCnBHbbl5eq2CXmxpfrEnGJurp9cLJQjAnsLxqOLzPcvGv7keXAeHF69DABVNUwESOhDXIXKOiazd1bElGjczdjp57F6Wd19IOilZkNyLRtVrEnZE3W+48T+MykWJWcPBAhFg4ecghGgirSwZ7sCYKDSAKZO1zS/WvxskbXhTn38uRpe4GM3FfiV8QICZW4zuKPzfO3F1dFE5uuYxYZflgRHBNOL0+wC59nwGbyc9593UXftTS2JKkWxVY2AKzSTsgOl+B0FqOzHBnnOj1+ygjw2yb9P+x/8wacWQXwvKq5cJfX8EeLfjgEy4RXgTr01Lia2IPwC7YiMSaldcUA094l3yblHa3XomxKbCk87J03daEEJz4z2L5g82+TOheduE0VLsUwqMDDq4o1//ZA8gqQoPrWV02Lr0mIElWdmtYYxl8/xkCb3rdm97XzqqRFD0vNfsfkvRmk7ZnNNboBia7+LMmD39nz4x6E3QFyysWEyrhDYeERwHRoSc7T2LpCGqydtaTXpS2UsSTYqAbwYNsMGItQyqamoIVf5CeJeykTUyq1uXg2QXIAYSRsLrhCHMaYUJlecXCJUNofDBlI3m5QK48UWXz720Jq/dqww20vyFIT5GZ8CaFLe5F3n08s61hw+PT3xX1lZ00eUUcz2gl3nEvVsE6M/bCux8kvP6I+z7qiPDIRWDTScTRkiRYoUbX/Yr7PE4tWfcPGF0TXcafmJUxvXY2maQXNttJSLLsDJx2aP3nrNxxdilERxAdHhCcwA1/eUzwt6fwiEvYwKWsnDHFCzThDflykD6u6N5UoU6bofaVHYb0+oeIPA94Ci5G/17gGar612XH7AXhNdLuevjh6bTKRA99xSsTxIIeLIgODeER2ANZl3OKiSoADTSVp0dCZHahbJqt/hiDZyo2yHBDlgzhqUqKBHWjnVFxY8fRXPHLzlSw/jNxwVqiAzfTw6MAOV0gVlmcKNGB90E2nHlVP0h9ovNKzqVj9Uh6CjrM0rIfUdXnAIjItwD/D/D1ZQfsPOEN6a8bW0s/M372AY+/t953Fzgt5/C+iKO/vcrpTUd86AmXCI+E1XVrzU7FR0oD0IWuSTDQVM06dcEE0XUKSZrk4plaQHjp1RjxoBpreQnbupcuHKsmBQYi/3eg6AIkFIJTN0Z0wS09u++xS4LTJZfeE3LpnmscP+yIqzcHTow4am0yMqSu4dj3NkYnpAe1iK9P0hvCpFXV+1P/XqLGWe804Q1TtaL3IfLR4MHNO1aRpKSTRIqE1kVjL/o1sZk7n2h2KbJzvjNdF+0MvDCBniG75FZkiS8ZIH73/jXVxHdH/Ld13+tCkzw9RNbfp3Lt4uBtrIFGF8STtiChRSK3LvfM9WxxLYdOYalranczVj+k16DLm0Xk7tT/d/ptHWpBRL4f+CrgPuCzq9rvNOHVQosHfRCS62IipXxRMbkhzg8HcPHeiOUDEdEFwwOPv57TS4bwYspnF2t2BuyB0+bs0pEegTqCE9YEFzhik8Al04mAiLpyUJ704tJQ2d+ieAJootk5x7bizR8FjYz7LvJO70hQq0jgyNCEICKIcXKJdbl5gguyGODqwwJWly6zvGK5fM8pq+sCrj3UnbAJ/TVLm+XpYEbR/RiB9GK0Jr9tfjC3hPsNqz2J7lXVW4u+FJHXArfkfHWHqr5CVe8A7hCRZwPPBL6nbLCdJby+tLvByC5+zz6QdcYvm4BK4nsLTpTl/SesrrvI8Y3GVxxhY1VDVqvTAFioC1jEZBcoYnRNdMbVIzGB9YQHxtgM4W0KqIlm5whPVbDWk564dyvW++wMmqzcEBT1ip4k5qldgoSKifPtxJFeeNGtEFkcw/L+E+zBhWSFSLzsrXTHtjLSo+D7gcilr2hv51pe2k2xbVeqt9Vs+mLg19lHwmtaBr1en+1kqTW2Fnwef1b2fVXf/rO0omJO1xFOe7ggPEppduLJzS/Hilc5JJrdwmt2aaJbWMQoxliM0dTLFQEKjCUwTtMzFYRnvWYXebKLVLDWaXZR5Kqi2EBRK6gxTvsLxRFg7OqzYE4FFqxXheDeXfK0Lye1DECV4Dh1bJyXV0drq7pf2TZZ4it6DrZAY9KrScZdk94Q6X4i8nhV/TP/71OAt1Uds5OE1zU6IbuyPqr6r9IsirSKzMRTvOW50qQcu106zS46JDXpPcFlN9mJzVjjyc4TXRC498UiQgSWQYSIEhj1ZGcJPNnFRGcyszKunhFZ40gv9XcYGVSFlXfXhaGbgBGAFRQDoT9db/LGUVy7UMS4QIbgc/aCOLfQjRmcalKheSOSXHS966BNHx1pgnWS0ceFDBWl/UEReQJO738nFRFa2DXCk+61u63Jrqv7us0kjP1SBrDOlDWhEl40rK5zycXJSoU4GpsQnSODhOyCNdmZhcWIslg6gjtYbL4vPdktjCUQ19aktLwsrAqhdeRmcX9H1rAy7jMRXWuCfi2mNeJiGRiIHCHGxBYXMQCnsVpxydUmgvBIuPZhh5jIJSTbpXB6ndfuijS0uujqmemA+GqTXg2NdmNN9LYYgIxV9f9uesxuEV4dDEV2g/yAVcA/xKLebRK4vxfXLMGJ5cFHHHBygxNUQpwJuGQj3cSZtoouvGbniU4Cp9EZr9kFohwuQwJjOQgijCiHQegIUSwLEyV/O9LbLIFs1Wtz8bsNCNUQWsPKBlgVAmOTReeRWkQCIr+8Q0Wx4oRWwf2mW0F9ZRW7dMENs3LnGl4QwqOAow9ZLr/7hPAo4PQ6gxpHiGdq542BDkzerkmvE2ijoMWg2D/Cq4nW92NC9zFZVx9ngijrMk3qorXRoZ/glk0TNt5H1vi0E+MIT4wzEWN/XRBYlkFEYDQhugMTOeIzEQsTsTCWhViWJsKgpYRnEVY2IDSW0BpCCTCiRF7zc9pbhLGSCnSA1QAxPpiiQOCCGRgf6fWbCWX3u7ULScpRpZevFZm2vS8DK0NLQhpV5iJMTR6P/SG8vjW7kYmuzmJ2E+Jy7qzbxSs6FMJLLogRnG5GYpME44UnvIWLwpqlxRjLYmFZLiIWQcTRMmRpLEeLFUaUi4tTFsZyYEKWXrtbv7vktyIf3kqdNndiF4Q2YKWG42jpNLzwgMgajCyJfPtQXEpJCCDGp+oZVBVBnGYXeY3CF0rQhTN51fgUHVUX0Y3U7bWQupd5ZDE66UFjwqgt81Ba3tgTpgD7Q3h9Ypr3DtjU8sT65WMLITQ+MS3Of4v3ofCaTVoL8mVVQBxRxVHYOPqa9tU589W9H5qIpYlYSMShceZtTHjxe4yVT/4zqgn5BaIYG2DVEKrhwESEooRqIAoIjKJqiYzBGAVVV2zAOHmdlpoKUfvz2qjgYsAeBtil8dFcddtLbnFPByHEwYipJ5TtMTwizh3hNX7QByK7NpG3dVFPl5IRXHUrFE5uCFxBTOM0OxtAuPRaXeK7U7dywQcpYr9dsLAsFhHLIOIwpdktjHWanVguLU5ZmogLxml6R2bFUqLkZcQSZGZrhGDVsNKASA0nuuDYLgmtYWkiVjbA4MkOWPplEivj/hcJCMWntYjL4SNy8iuewAKwh25JnVm5AEa0FK4+fJmsJ3ZanyRBj6bENahrqiHpTUbL6zAPr2ucO8JrhI7vWRGpdRUpjp3x8XaGKt6nF7BRCGD97gMVgvPb+bSSJOVEdFOzi83WlFa3NBGHErI0YUJ4ATbXpI3EYNQ6bc86mVYERD4TemEisLAQC8bl97lghhLZdZ5fFEcTxfkgXeknWZ9bAPjafRi31lZXyuKabqwycRd//XeT+zCYKbyjmt5cALRP1HxQGxFLh2TXdCJtiFFAjvHCehWnwS1WsLw/BCNcucWwuizJdossXb6aBr7OXZJvh0tBCXxScbDW7g4WIYcLR2hHwYoDE3Hd8pilWC4tTlhKxEVzypFZcWhWHMmKQCxHsirV8I51SaSGY7PkxC45liVGlJXfyGLltzI7tQGhNRhRFFcKXsFpdjizXUUgClDU5eNF4jYSQkHc+dslHD9EWD4Il97jatMf3xiggS8sGmmiDbe9Z7n3p8sJ34D0JhPAmIIMOdgPwusaPWvj8UNZZ1LltkkmlSKpSiNuDwqn4dmlC1aY1P4WGvvrUr478SskJPbbpbS7xF9nzgYmlhJlTNmQA/+ZweXlpbHShbO7LZx6rc6KcVqeV7sWfsmEEU3GjSO3gVHCCIxRrI3X74rTUkXcOVkXtU37j6xfTbK4BiaKSx6zfk9HbVuat3n3snPi6VrT692snU3aUTHW9S8at6k8Z4pWpjS94NRpKtGhSTSbeCObKA5W+GPivLtkNYWJNTyXgrLwaSjx6yhYcSFY+SDFikMTJprdRXPCkdfuLpoTR3iJhrfW8iLEEZwajmXJqQZctRHHEnLsN6I4tkts4KK3q2CdowdwqKFfmSHJ0jRrFbCocZEZNepSVULHfxqst3uMq6aER4Hj+0jRlayrIKfuRV5x0CqU3eMxtK0paHljj1+Ec0N4+4QkMgvOeou8eWskicyKTUUrk/p1JIU9k4nu+3KL/73/LtbuUhHZpVn77GLNLia7+D1NeAeZMN0pBusFMVisMWs/ni7AwIm6x3EZ+/KMZWGdPzDW8kTiNbyBPxevXhnvN4oVXr9cTiyJxqcLtyY3/kFwZaZkfU17mKRTIJ/BoeJr908P54LwhtLusg93X+MmWohPyTDhOh0lWS+aTkcxuq5qbHQdrEgSjW3KnLXJa2GihOwOUgGKw5jwPNldlBMuGefXuyQhBuUgY9KeqsGKYNSyTAUtsHBqFhi1HNolGDjw/ryFWA6C0PnzAkOkQiBOmxNjEYwnPB+88KSf5Bv6slPiiS06MEnqTrIFZR8+twzGWPs6OtFOlORLCU9EPrnse1X9w27F6R5Dkl38Xtc/t/V4Asa6QpfOQS94n39irp3NtyPR+kTiyCyJ5hSk/GcLkyI72YzGHslpotldMiccSciRRFyUiEDgKLPm+VgtkZ8ExyiRMW69rIFTDYADliYE63P4jNP0rEoiSxBrdyokgdpU1Jmk1LxsrCZxF8RpvMZXTza4TYlUhn9GhiKiUUlvFwkPuBt4M3Cv/z/9aCjwOX0INeMsCh/elIYnVpFEw0kRHcSl5RySd9eheJPReAI0uNdCLAYlwLp0E3HvBz5AsZSQI4m1vfgFgQiHkl2oaolQrKuBQsQKK4YI4UAiLI5IrRj3rrJepublcUvWYqJbFx4lveg99a6e1InJz18ItxrFr8YYoGp2FqNrX0NgoudXRXjfDnwJcA14KXCXqj7Yu1QdYWqBiu4G8O9+4gQn3qyL96hIbbOosia+RBOCRNuT1KoKkfV7km/n18ou09FZWaXew4TsDgUumoAA4VCWGD+YRQlYEYkS2QiIiBBWEnKEcCwHRJg14ZmICMPCRFgEI5uyBd4EV12fQ5rQk7+NS0K2ARjjT1kV45fZcZhacTGAK+LcYMKJx6X1IlT1x1X1/wC+GXgU8Jsi8nIR+cQuBheRJ4vIn4jI20Xku7roc2tM8z7lIi515PxSmgQvNpBRJTTR6jhTpDNOFk4nDQepv+Poq/vcRWHj11LgQIQAwWQeK+M/CxAORFgKG8cGSVTXbhQdCFLy5MnmzmP9v4pukvqZCxYHeFz5rL3XshiPdxLlu+I1NGoVyFHVvwBeAbwaeBLw0dsOLCIB8JPA5wNPBJ4uIk/ctt8xMOaPmSiYE4s5dRtsS+woK5v4qactNg3jV7qmXWxKBsSmrSbrZZcS+pfLz1sCS6/ZLWW90W0ghkAMh7Jw3yGurTeP0/0Y8WOkxlybspvVlNcmLfmzJxWFjs1ZE/oiAicWs/KZzLtkBezQj3GS41j1GhhVQYuPBJ6GK5/8LpxZ+wOqeq2DsZ8EvN2TKSLyUj/OH3fQ9/lBPNcVsG61wZlIcSHxNRsqW/IpRpxvZ7w/zPQ4M7Naaa1jYo0i9nNayFv7dB40vqEw1WtZ5cN7O/AmnHZ3P/Bo4Bskzl1S/bEtxn4EjkRj3AN8araRiNwO3A5wcOmmLYbbbRQ9QKJgwnhXatbrZ0d+4AJZF/JsgiJSLUTVefofA4nAnLrCqBI5317RUENE2fceE72AVSbtvwLuwk2ny8B1mVfvUNU7VfVWVb11cXRpiCF3E6qu0GWkrcmuqkqtLSgRHMXBiZIV47amUEVjxMiVsXCpg3uLNeA4aIFqbxtQz6C+OdvRLRCR7xARFZGbq9pWaXh/CrxaVd/fjWgbeDcuEBLjkf6zzrC3v9RxsEJJTDNRwFqfZ6aIlaT8kWYfrjj5TNVvm7iuLmxV1i/cK8L4d/f5SgMCXbDSBQHKSiNWEhGgnOgKg+Ew82itNMJiWaGsgJUvFbXy/cSFQaPMmBbZkAk4I3NyTmlsmLFuOZl4H95GPlzG/J+qKdYW+56HJyKPAj4P+Ks67as0vEcD/0lEfltEnisinyrSWeLSG4DHi8hjReQA5yv8lY76Pp9I/XLmPehOwym+fTZFejGi1N9WDRH+pcavkXWvlcKpqs+1sxsanfWfRSinqqyUjWOjVL9pDS9Sv41jiuxsRv5kY2/cuWVJbOOa2Fi7K7uIE8YOyR3/0FS9OsC/AZ5FzatTquGp6g8BPyQi1wG3AV8DPF9E3gr8BvAqVf3bNlKqaigizwRehatg9kJVfUubvjpFi8hddgL1plVmZVNABBsYpw+dprS9ZJI7jQ/jnwgvq1pxKxKsYEWIrPF15wxq3CY7SVFOa1gRsDJBUs8uwHKsSwKUY4184CIisBGBCFHm+TtRy6kqV1VcWXcNONYlx7r0ml7q5UtDhb4aslXj97Jdv1srydKx9TIXf87EE0rWhJe+bAt3zfDtkgTtuO0AVsHOkm5d1D+/m0Xk7tT/d6rqnXUOFJGnAO9W1f9VVw+ru5b2l4EfU9V/7gf6GOCNwH8F/mHNPs5AVX8dt1v4XqFvUzo2vZI1tQtJyihtNuSstpP+PG0i4qorpU3Z+O+IdaXilQac6gKDZaULjlGWGnIcR2qJyW/z5/tYlUjZILvYnD31fabHWGngtnLUrGnrRY9NWs0ntTxNN/bbuX0vNsu8D2XO7j3R0TjH7l5VvbWwL5HXArfkfHUH8N04c7Y26hLeY4FnicitqvqvVPUtIvJmVW1NdjNaIuN/UgPhRUNwLCw+dIJEFjh0y6o877j9W3Gz2mt7eN+eVcF4X561boG+064Mp9btKHZqU9VMLKyMW/t6YH39OusWp8Ym6lItq5ziASsc2a004Ng67c69H3jNbsGJJ75wQ9Nzr8jLGJuxiQszfiXn6DVbXxYq3r4yuBpCaFldf0B05JbgiXXXtG8eOg9Et4GOfvFV9ba8z0Xk43C8FGt3jwT+UESepKp/U9RfXcL7EPC5wL8TkV8F/gk75VHYT4gqNoDVBQNqMcencGoSIkxWFgTiioV6Qkg2sBb3rjgTMfKmbfwK/XaKsWm7sm5joBO7xIrhWMKNqicHEmHj8lCZdJTTuB6e+np4epiQ3ol1Zu2JLljZgFMbvweJaZ3IpeLNWePOwb+c+boO1CQ+ovhaAObqKVhLdHSJ1QXjfxTcnhjugrJbyb0V2OdqKar6R8DD4/9F5B3Arap6b+FB1Cc8UdUQ+EYReQbwO8BOJMWNFantZdxsf+rMs+gAglNBjk9J9nbw2ov4PWnj9nFZ+LQz3wV4jas84gMFod9JLLTuHbvgxLrH5VhcGadjX84pJj3LigghUGUl4UYB0JUuiJB1ifeUhhdrfCsbuAKgqU264427IxVXANS6KiuJPzL202Wd4J7sk0IKAnLtBNTtWhYeidvWcoD0lHOn3THdc65LeM+P/1DVF4nIHwHf1I9I3aNxjboOHdddBDTOrJpK+azsAlaXheBU0PsfcGOYD8cewOIqBCdxVRCvwITebxfvDRE5xoqMICKEkUEkSIqAAsketNesW+RvRJOy7CtdcGqcH+5YlhyZFQBLCTdkXvninsfWaXjHeuD2tNAFV6NDVhpw1R5wahdci5aEapymFwWchO49jAKiyBBFgkY+aGFxGl4kqaKePtF4BfbA700roPfdD+LIbnVZWD6oBKfkRrW7+MHqdNLX7GsSRKOdRWDrD6n6mDrtahGeqv5U5v8/wEVsp4GBImtNUZQakjeRyj7fQPpcUyaYGkGuuwy4zWmC47V/Ctaa3UZeXpzS4Z3/1jpNyti1lifWRWkNQbLBTrzhTqzhGbXAARZn4hqxRJkE4vQmPnGAIvbbxX/HpqxN+RFDddpcZCXls/N5eHYzaCFWNh3mftVJcOxN+8uXwLgodPb6Jde2y3t1njHRa3EuKh63Qo8kmp0w8UQp00Q32nsis4Ek2oxdwrXHPQyxsLyiLrhwWbCH/pgotmPdWleNHFFoaFBVbOC+C0O/YbbXCG3gCnCGnuwWfoOIAxNiVTjRBYd2WW+bRtYR2JOUKes0uoAHwwNCG3A1PCBUw0nkNLvTWMMLA2zkTdrQl4i37txM6DQ8CX06irhrEhzD4QPK8ppy/LiHo4EQLb25669hTIzpYqC5WnWNe9gL+ui7b0KaCW8H0TPpbX1cMjnd6oHoQoBESnDiSkVFF4TIK1pJNNIHM7Cxuuf+VwtW1isYokTb85vpGLwvb72VYoB1e8p6P54VtxOZEUuQ8Y2lN+JepSK1TouLI7LBhu8wDlJEdi2XVfERZ3cOkk5LibW8ODrrfwwOrihmpURHxpXAJ+W7y5BWEzN2EI1uosRRhalqu/tDeDXJqbFvZqLmMqxNVOM3nT69bAhOlct/dRU5WXF63U0cP1QQb+Lapbgdy9SZvaqCBo4wrDGIKpGoS/nwFykKnLYW+P0eF2KxCAuxrIKApYm4JkuW4guFerYJfLVkwOXyeRN3pS7VJJ16EvvsjkP3fi1cElnDcbggjAyn4YIwNE67C43TTkMfnV15zU7X18KsnN/SHsLifcoNb70Pe7jgyqMvEh3IxjXLy7/bzN0bJ+gFNCa7qZLMlLA/hNcnRiC9Rg+v127sQtz+FtdWmCvXEHsT6ve5cekpbGhDcb6ai9q6tbVWndZkrUtREVFWfrvE0LqNeBY2chvyiHWmqhFfrt19DmcLdcZLwhzhBazUJO+n3m936tNPVlGQaHbJqop4eVkSmfWpKJaMhrc+V2ucpmseuAZcwAZ+k6OGTvVRiGTXyWui8u8X4fWl5cV9U6//PlA46XS92bZdQKQQ3nSB4MIS1G1ALaHbm1Xt2pcnIoiBaKEIgq5cbpuFdTCDFFEZR4ZG1C9Ds5xGi6QE/EJcteKlWF+sc5NRrE8vWfmlYknaic+3i6zz2VmVRLNb+ais0+4CbCjoykAkLtpscb47u/bdSeS0V4iDGEp483VEF5dEB67ce5K0PcXAQxeujh76b4QRorR1sV+E1wCtTZWRiS+LJEXF++fsQggvBklZCHPKxoqLjdUXeA0PhcARiFpBrTqNzfvCVhJgrSamrUiA8QS2UMPCui0djShhBeHFeXWxv86qJES38n67Veg1vMgQxUGKOFk6jszGa2XjwEWs3cXVUUJH8IAju6Mg2dhIQn8dTMtnoA+MTbRdY6Lns3+EN5T5OQHfnmTIN/4/vBRgDwwHV5TgVFldFlaXvKl5ChoIalxjWYkr3eATdBHjzVrHFurNycirkSJKtHCm7tIEfv9a569zm/+4dkUmbeiXhiUJxdaw8p+dhkHybq0hCl0bG/qobCROu/ManvgtF13OnfNVgtN0l1fU5dmtYHXdguggDla465VcurEnZgfjT0q7g6Tq/hSxW4SncUSygmkamLawxa98epyc6GnvyJ6nl2F1QTAHwvKK5eA+i10sOb2OJFhhFWThTFsTuQBFFLiu4lUYijNvY2K0PnctXo0hohsBjUDU73HrU18KCC/2xUWaXsbmSHUVBai61Bi1TsNTK2fIjjgNJdHoJDHbwRF3cAwX3h8RHQqnl83GNo3uRBj2nmUJoIPxGvt5K/vT7khxJrzpYiviy7uxZWZvlUlcdxJq6j2t4fnghRpY4jauMaG6lRiB275RZB2ljTyhxf8D6AIIHelhlAj/Q6P47RFdH6F1ROe2TlxvAgSc2XsiXWA0SXvxS8RcQQC8+SpYT3wa+vWycUQ20ey8jyj24aX/90vpgpUioYXDwEWnPcFJ6no1Qp17SkGbsmekJaaqQQHr6zxB7CTh1dbyoNHD3XkKQlYGLfmu6Pg6k8d/Hstvl+uPJXJmbXCsbuL7RGSzAo0cAQLgTUIbO/R9RFeNW9Egfo/XSCDy5CYxycX/y5roiggvqXQSFwHwKybi/1FQr9kReR9dJF6LE8yp999F3pSNl4bFvrtTMCt3ziZ0MkSHrk1wWnIJi+5H1cQt+76HSd+YSGq273xN8Ry06Ba1SK9Vvz04sutMim3GzGiocbwgOhTEBog6n1Z4weXiAcnyMhOJW2mx8A79CKe+xd3GfVtQdQUGbGBBfAUWwPokuE3Cy4gY95NazhbX4dM00SmO6HRNeBL77bz25qKyklRDEV0TYHCsLI6daRZectrdRiLyECZbT9pNX1pTHwUUZg1vLLTU9GDgCF6J5gaUyp+QtG/jkmqV00uGk+sMR/dZLv3NiuOHLIiOAmfWes3OnPjOj0AXrlSShs6sTYpkBuoc/YE6ggp8sb14TE92iTB58iYmeJpw03l1/h2fOuOJT9RHVdMBilC83OvobOA1uwsftBx8KOTkpgVXHr7w37kfx6Rqit0Qpfy61/28R7Qij7EJZ+zxC7DThNdIy2vhuxk1yz6LNr4nT0hxuopbSK/I0pFZ0nVcCNNrQhicdoX60vBONRIETZxgkiK9+HPWql0h4cXvkipbFScDx8RHknLitDdZp51E6+OS2ncRmNN15RP1PszW9e0mNFn7JrteymNpMxmGxE4THrQgPWil7SVdjEmAOQ9RuvDAGd+5d+aHR0J0sHDLzt4dEl40XPmwALt05ZPA+8IQ7NJrejGRGTx5iP88dQFiQo3zENIEU0Z4XrtLUvXSTm715iqsa9x5MjQhSQ5dHKQITpxmd/Fey+KqZXXZcO3mhVthsvLymfyblvuDNvJE3coUbHhsn7UAZ5N2T1BW0aTP8VqvCvCBBxs4YjAriwmFYOU0NV+mbr2u1PvwJCYyv8uXGiCUDd9cTHSurWwQXlZeyRCeeBM2SZz23yVE5zU3VJJk6SQKG7ex8Tm5dwktKoZoGWuEvppxci7513CMH7HsfR2aIHovfDoTXn9INmfpUdPLHzfTreR/VzShqto0XcS+sepC/G5liQUqRAfCyUOWSKRc/uuQ6MBw9eGG6EDWqxAiYOWity53LWdQv0JBF3pWeyogvI3zSZGciXzCcKy1wZqgYjKI1ukmJlq3DU6Ui++1vhKKEF5YutJPmQhhUu6+5iTs4t7ltckbvxOimyC5zEvLBkBfkdv64xd/XoesmrTJHaviwVcjhIewOFaCa44ozKlBRbF++YGoi8gadSsWNkzNVIkpV2tAkvSVs8JWyBeb4n7ta6LR4ceV9bmapM7dWo64ovHimsWsLOHFJeHhetBCDSbHF1p2XfPuSxVJ9UZsW2KIcvbAYD48EXku8HXA+/xH3+13QizEXhEeNNT2im5KD5xZ54Hvqs06TcWZc24DH03Oyy6Fk5sWiMKFD0SoEY5vNERHLrhBQJLTtuGn8zCxDJZNrW1D0HyZsucS7zkBJKshTKzxZSaOIzklOIGj+9xi4NUlA2I8OaccmSnfXdpXWISy69qG4HpFi/EGIzs2XbkD4N+o6r+u23jvCC/GVtpeRybvYMglE2duxn65JKjhNbTVRUOwUi68N3Rl3I8O/Bpb973TqnwffnVGTEhxMCQ41UQrazQJU/40Xbh8wSQ4kurfvbuEZxustcHFiXL4wRU2EK49bOlyC9UTuyfJJAUlT9vYtfsLrTWmIYluAxPQaPOwe4SXY5IUYWsTd0oTI/0A1ZUnNlFF1mpZbEqiqIHTGxeI9asxVkq0dDXj7MJFcJW1BpYk+VpPRIIrPJDyv5VqUylNUQ3OF+irJYs6UiX1fWzGBsdONhO6pXIonNy4TFJPEl8lvn/jfZhNTas217hvDEkcHY7VQOu9WUTuTv1/p6re2WCoZ4rIVwF3A9+hqh8sa7x7hAeNiKhxQGPL8bYaowGRQXl79SshwEdpTbxna0oLEji53u1je/HekOCadcvPjHB6fcCpN2/jTb2TysLROnKrpHxruvkObMiQ1rw0IPH/mciRp0HWJm5sikdw8KDl4P7IjWshumC4evMCxEdn7Xo8F3yJ8/o2Z92ZlJrs9ezo2m+NjoinkXbXNbHW7+9eVb216EsReS1wS85XdwD/EXieH+15wI9SsbnYbhJejCG1vYbjnTmOGsc2nUyFQRI90y7pMpOi4UhMCeP9HmRNDMGxT/sIN2VS2Yw6Zv14KuvmcV9nIthKkkAc9yUr3dAMzCnr9bBH613QXGXnnJzAiuuSvSaNUfeYsvtY9F2n2lXDznogu66itKp6W512IvIC4L9UtdttwoPGpAcdaXulAxW0r/uwV5hWTck7duarCDZY56ihrqqIitsPQ2M/WeRWLRw+ELklWQY0EMJD58+LDlmnesQFOFMpJckpZETciHZa1uWoJDZd/dgn6ot4OhnDC77EU+zLs86PhzfLwWub8freVICm1rVJ5Ms5qOh+Z5+7vHbpNmX3uCOMTnRD9e0hIh+uqu/x/34x8OaqY3af8KCx5lX5kHchT93v62p9GbQl73XCqyR+vPRYyfdmvYVhHPxI2kSpthkST8TNilWg7Ujcj02bv+7reIexJNs51T4hug3NsdksK23fVJvbtp+WaB2U6F2ufvv3+GER+UTc2bwD+OdVB+wH4UFrc7N38qvClg9GE/nTaRuOxCSVCrJpTtoFREuzEWiISzAFfuF+nKyscQS3aPgC7Scp4OnN6nhv2NCbr8m5xUQXqSvbLikyTPsK42cgT1FLXZvOIpfDTOoNbC37UDIPMI6qfmXTY/aH8KC9j81j7MTlQRBfo5Jrpd6Pl2h/knrPe5BzcvWyX+fKEX+faKuZvjTOo9Pq+ZP2I54xpffjno6WYtICU0i0zsN+ER50QnqFXQ84cbJypMdOR2DdByRRzXTb3HOJtbW0ticuMhv73ZIggtem4ojupr+MTXKKFa4KN1gyEeJ+vemcjJcNhLA2rTfMaM2cX9aflpHjzPXKHJO9t2XXf0jspDbqXRRTxP4RHmxNekWYpAaYo9k0rQgdJynDmryUFLlkyCxtvlYuvC9y2vtjN75OBTA2+pV128RqjVNOigJBRWg58Ye6950S3EiPqjBreJ2hdnWLJoGBRuPn//LnPah5E6Toga7yMRX66tJmZobESlEQKd7oLl7ylSGfjZQD2Xyv5ISMrJJzn9KBlCT1JV7Glg66qObf2yLtruqSNCCbOvexql1vtejy/h4aEyU8U92ke4jIl4rIW0TEikhh0mHh8U0vZjqa2DFcIm9+59nvyh7wsn4aI229iSSv3Kb+u2T82GyUTDAiHQxImZ8bRFT1suu/N0zYDSIjVYOv/Met6pw2rkf20rb8Eay6T03u99TQpVYWX6eq19AYS8N7M/BU4KfadtCqjllPWl8VutAcStvkmY05fr0ykyz5Lk1q/j3towOSFRYb7etcz1TfmiczmUmX9tlVyZ0dquw65hzS6rpv2a5TbDFk5+ZnjwrGthiF8FT1rQAypj9sRB9HazQl7KyfrmxlQl5/WTM3HaCo68PLQ1qMtImbZ+5uaKsZMQsmVakJWybvSD+IW2OC5DL78HrA1uXXi27KVB/2rL8uK2dJgABKHO8b0cqz35WVQd/6UqXOJ3eSpM/Zy7cR6c3001SzOyNHVbupoANC6ZOUzl0B0LJFv6r6igb93A7cDnBw8abytik/0FYY6+HPm8ht2meJr0SbzSO4MiSuhJxASWNUTLiE3DLnU1RZurDvukTXpO0Qz0UPhDSY5nXeNLy6i35r9HMncCfApYc+qtZl7HSfgqoR60Yl67RNH1PUNqPxlCKdxxF/FGt5iY+v+FrVvoZdT/7UOZ7R6NLa5RkTuKHq2WRSNp3ATVwmeW07JoxBTcwibX0C2GmTtgydaXtVqGGGFbbNMcly/8/TwrIBirx2RZ+lj0uTXp5fLnN82Xe1x8y2KyI3/91G+XV/7qVkV4Umroy63Vbd7wp3Q39ZBP30W4mJEt5YaSlfLCL3AJ8G/JqIvGoMOXpDk0lSJ6LV4cNTmjxbw8SsRFWbiu8rN8QpkHHwIhB5wZaqvgYmgbHITvzYdV5DY6wo7V3AXUOMVWd3qZ1AVvPJ0yJzNL1SIqggF93op+S4ptc1q+2UaXrUnxiVid5l/ZRp24UD1mw3AKTg2o2FbPHVqWD3TNomvpEMBjNz+0T6/Jv489KoYeZ20ledY+toUUX/1+mjKaY5T0vR5sehV4ygzdbF7hEebBcZpPih2BkiTJ9/kaaX+L1SqSg551cZpa0R3KjsM+665iQobJcJZiRtCxz+tbS7uM+JTtAi9JIs3CHOXVrKINiS+LLoNLo7JPIc4hWaWl/nWb4UrPuJWrpPb9Ko4PNdvNdMn+x667MD7DbhxdjCzM2i6mHqiig6Ma+z510nRaaG/61SpnQ/mR+d2pOxjrmakeeMNlcSGa2dorLFxOzDPzx0+khfmIRpnYP9IDzoXNsrQhVR1Vn9kedzaT1hapB93dJGuTIUHVZEPm3IpYYJCznaXNf+uwbI3uci66DuaqB9Ibqk/wkWR4CR0lL2AXVD63nt6kyMLpEfvYy/LJAhk5d3JkrbNOevwXFlY2WLC5z5LOmjxi9Iy2ueR2LZe1r0bIyenjHQWEn5/orX1uOIfLOIvM1XX/rhqvb7o+HFGMFfU1dLS2sBef6s1tpejpZXOOE3AhrVWl3tpWdNI7vZ+5STmhLLmHtMHtmdIaK1UHlbV7ZBHXKalDk3gixxHl7v44h8NvAU4BNU9UREHl51zP4R3oioOxmaEGPnKFpNkDVPq0ixqs86ctTQjjdkqtNn3bYNMRNdk7F1KJP2G4AfVNUTN6y+t+qA82PSauo1MtJmTGc+nQ2TOWtfpV5Jmwoiyxu/jRlbhrr9lZquBVpoh4GKsjSmZP+PCTxXU3zGa6y0uFlE7k69bm8wzEcDnyEivy8ivyUin1J1wPnU8DqM6m6LsonSVtPrbf+FvOu2zbUsibI27qqHlBeo56MdHVOQIYv6Mt2rqoVVz8uqLuH46yHA3wM+BXi5iHykarF6eT4JD+prMCOja9JL+mrQ54YMedHYbLCjDGUToSiAEn+dJrUc4i0ivbYViCdLdj2O29UU6OralFVdEpFvAH7JE9z/JyIWuBl4X9ExO2XS9u4MnZBJkEatc54KWdcktEbfpZt1cZ5VJi/l13y05PSen83O5pYCkdZ7bYdfBj4bQEQ+GjgA7i07YCc1vLq5TVuhKIo4JdSVqcIUrYrENi0S2go5OXd5paLS8pRWUOkxyNG7ZjfgD25f5zKQ9vtC4IUi8mbgFPinZeYs7CjhjYI6vqoJ+QZHw9DnX2EGt0FZkGLrMYpM8uz3A6Bfa6n/E1HVU+CfNDlmLwhvsBJQZQ9rVtsYYuIXEGxuwMK3PeOPy352pq/MB9vk4TXUqMoIpnKypjTGrvLxtiKIDsz5LjCU33ESwZwc7JQPrw4GudB1xthSjsq0kcLjclJSMhj9Ydwyebm2/HV/DErQ9j5sdlLw2cD+4sHuuzZ4DYy90PCyGETjy2pyeYm0W2h7tVddZLS83Ohs2cqJ1HhFS942Pi8y28ui3nUf7Ey7ponHuRph5r7UidjmLRNrhbF/WBjnx00A2T4g0Qv2kvDS6L3kU5oAilYPTMC3V3gNUrIl16pK3rrn0+SZr2nyVRFauRncfBKOrg1vgTFlH2Uz8hrYe8KDBtpSW9Q1cfuKbtZplkf8VdHYMnR5Po38cdWkVziMyNmJWPPYRhh5ro9O0iOZq3VwLggvRu/EV4UWJm6phlpXy8pqcJljc/svOi7bpq4cDVGVB9d2UueSXgsZigdocUxHGJ3oEuggUdo2OFeEF2N04usRdYIWWdTR/ooPrtnuPGCkOT4doltjijLBOSW8GKMRX0PNqK4fsm6EVtNmXF6/dU3dLq5b3Zy3PkzPEmxTvGEITJVQEswaXgfo6Rr2HtiYsd/YJ7Lrom+do7TdoSe/0WDJy2n0Hb3N9F92Xhva7plUlwofHrQ7j5I5cWZSl0Vqt8U22uOUtc4m2LPATRF2j/BiTCDVY+fRRerJjEGwU2THnJbSD3qalIP69mqcQ5c7YrUpD3Xm2G1QYx6UaXdt+2x6/FbbP3aIXSO6dd8z4fWDHjWRKfj2thq/iyVRdfuo07ZiDkzeET8wdpfsgHkj7h7RM+nB+MTXOXJ8oXX22ihtt8UkKpzcLfrsq/rxkNhZsgMEnU3a3tGzz6kXba/zwEt5afdJbSKe6WsIFE7CtkGLpmuF63TZ9/UY6nrbaap4+0N4sJuk1yfqRltrfr9t+y7QePewjeh7/VUWY2CXtbozY02T7/aM8KDX5U4wIRO3sOpJRsvb9keg4vhtSK9ygreYpFUVVqZKertugmcxxWsM+0h4u4QRUj+61vLaHNPp5O4zP6/h+G2xN5rdxrj9DywiLwOe4P+9EfiQqn5i2TH7S3jnRdNLodB/V2ba1o2ubnmejSZ1kautaTJyjm9uaprHXpLdQMUDVPXL479F5EeB+6qO2V/CizH1BNoh5GtIaLkaWw3Ttk80zs2r22YbbNn/fpKdH3vApWUiIsCXAZ9T1Xb/CQ92Om1l68DAFv61XNLroO8z6CA/b6tUlBHIYW/JzmNgTfozgL9V1T+rajgK4YnIjwBfhNta7c+Br1bVD/U66NQ1vabIOZ+qtJReZOCsHI2PH+q4iWDfyQ5oYtLeLCJ3p/6/U1XvjP8RkdcCt+Qcd4eqvsL//XTgJXUGG0vDew3wbFUNReSHgGcD39n7qFPV9NpGQrv0v+USaM3z6dFfWpcc6qWr5IVtm8mz7XHng+wAW1uYe1X11sKuVG8rO1hEFsBTgf+9zmCj7Fqmqq9W1dD/+3vAI4cbfLCRRsFWpkTOob1N0Kb5dDPKMalr5YMWdV7b4zbgbap6T53GU/DhfQ3wsqIvReR24HaAg4s3dTNiz5pe1/680v5aanm9JA2nn98eViEU9ZdLjH2mi0xNu5sU2XkM58N7GjXNWeiR8OrY3iJyBxACLy7qx9vzdwJcfsijuruK++bT6xGtyHGLO9WEaOqbvEW5Ltv1OzqmKKcC0TBLLVT1GU3a90Z4NWzvZwBfCHyu6kjJUT2RXmOCqLFov1lir7ucVSsuSiOxOW2TJj37QXPRQR7f1m1jtFh72wuBTpHsAGfSTnNt2VhR2icDzwI+U1WvjiFDgl3X9JrI3sS0LbkueZN3GxLsdIlZn6bsRqP6Y+0K2XW7AmaabDyWD+8ngEPgNS5nkN9T1a+vc2BvvqceSsbXlnOb8UuJqV6aShfbMA49qZv47QrN2W1kHnM+T5NL1mgWpR0UoxCeqj5um+N3hfQmiTYpJD2mnRSOVYJOghSF5Nj+2EHQ09id/2DNGl632AXSm4KWVxe18gj7/lHoIVWlae5d18VD+yqU0CX60c5nwuscu0B6jcem5fhN5S4JTlSSXtK4wXgdYOuJ2cUcbNDHuSa7KOqh4+2x04TXGzokva5IuW0/gy83a4se/Gmdr+ccS2mZprJUjlnD6we9Ld4fW9PLGXso0qt9Tbc1wxvKU7ef0py7krHH3qFsyHH7LyU/E965xRRKwxeSXokZXUvu7HNd1L7F898mXaVSqxtoHnZCKNPkjBrQOUrbN3rR9DqMTjYmvRItL2lSM2l4fWyJpldjvMJx82RoiG4KhLYju661nSmT3SArSBR0TjzeYZynlJUtc+7qbvXYGh3763ZmCdmuYaClZU2xd4TXm/nYUepHF1peuj8o0PRKtDx3TImmR/m4ZZjkAvmBNLtO+uyRgAcjd9V5m8YhMWXS6wNtSrJXYshk4yr0QHZ99DFlbXFw2eagxbCYQqAgD33J1aQQgGtfoell+yjopxc0iuAWNO5qvg05b3fZb5eBzhre8Oi15tuWi+UbrcCoOV4vJm6eLMnB1Yc06q8BSn12NfrtcgnZuQ9SnIHOGt5YOG+aXiFqRHChJvFl+xwQXaSdTG697F6RHczFA0bGedP0tuknTShTWqExWH5dw36mGqQY05+ogM5Ly2bkoQ8y7m45Wwutr2PUSjcZKagw1SDF6HKpMhcAHRm9LkGDyWl6XRYByCOdPkiwcS5dl0Q3a3adQmeTdsY+oeuiBANv3Dyjb0xUw5OxtpNoAxF5H/DOHrq+Gbi3h377wC7JCrsl7y7JCv3I+xGq+rBtOhCR38DJVgf3quqTtxmvCXaK8PqCiNxdthnwlLBLssJuybtLssLuyTsFjLIR94wZM2aMgZnwZsyYcW4wE57DnWML0AC7JCvslry7JCvsnryjY/bhzZgx49xg1vBmzJhxbjAT3owZM84NZsIDRORHRORtIvImEblLRG4cW6YyiMiXishbRMSKyCTTEkTkySLyJyLydhH5rrHlKYOIvFBE3isibx5bliqIyKNE5L+JyB/7Z+Bbx5ZplzATnsNrgI9V1Y8H/hR49sjyVOHNwFOB148tSB5EJAB+Evh84InA00XkieNKVYoXAYMlv26JEPgOVX0i8PeAb5r4tZ0UZsIDVPXVqhr6f38PeOSY8lRBVd+qqn8ythwleBLwdlX9C1U9BV4KPGVkmQqhqq8HPjC2HHWgqu9R1T/0fz8AvBV4xLhS7Q5mwjuLrwFeObYQO45HAO9K/X8P86TsHCLyGOCTgN8fWZSdwbkpHiAirwVuyfnqDlV9hW9zB85kePGQsuWhjrwzzi9E5DLwn4FvU9X7x5ZnV3BuCE9Vbyv7XkSeAXwh8Lk6geTEKnknjncDj0r9/0j/2YwOICJLHNm9WFV/aWx5dgmzSYuLKALPAv6Rql4dW549wBuAx4vIY0XkAHga8Csjy7QXEBEBfhp4q6r+2Njy7BpmwnP4CeA64DUi8kYRef7YApVBRL5YRO4BPg34NRF51dgypeEDQM8EXoVzqr9cVd8yrlTFEJGXAL8LPEFE7hGRrx1bphJ8OvCVwOf4Z/WNIvIFYwu1K5iXls2YMePcYNbwZsyYcW4wE96MGTPODWbCmzFjxrnBTHgzZsw4N5gJb8aMGecGM+HN6AUi8koReaSIvE5E/srnj8Xf/bKIPDimfDPOJ2bCm9E5ROQC8FBVvcd/9CFc/hi+9NaHjyPZjPOOmfBmtIaIfIqvIXgkIpd8fbaPBT4LeF2q6Utxqy3AlbWal0PNGAUz4c1oDVV9A27J2PcBPwz8vKq+GVcH7zdSTX8T+D99nbynAS8bWtYZM+AcFQ+Y0Ru+F7d29hj4Fv/ZpwP/ItUmAn4HR3YXVPUdKZfejBmDYSa8GdviocBlYAkciciHAe/yhT/TeClwF/DcYcWbMWONmfBmbIufAp4DPBb4IVyxgN/IaffbwP8LvGQ40WbM2MRMeDNaQ0S+Clip6i94/9z/AP4p8HHZtr7G4L8eWMQZMzYwV0uZ0RlE5BD476o6yZ3UZsyYCW/GjBnnBnNayowZM84NZsKbMWPGucFMeDNmzDg3mAlvxowZ5wYz4c2YMePcYCa8GTNmnBv8/2xLTBPEcme5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "outfile = os.path.join(outdir,'out' + str(Nr) + '.txt')\n",
    "x,y,valuesCF,valuesHam = np.loadtxt(outfile).T #Transposed for easier unpacking\n",
    "\n",
    "pl_xmin = -2.5\n",
    "pl_xmax = +2.5\n",
    "pl_ymin = -2.5\n",
    "pl_ymax = +2.5\n",
    "\n",
    "grid_x, grid_y = np.mgrid[pl_xmin:pl_xmax:100j, pl_ymin:pl_ymax:100j]\n",
    "points = np.zeros((len(x), 2))\n",
    "for i in range(len(x)):\n",
    "    points[i][0] = x[i]\n",
    "    points[i][1] = y[i]\n",
    "\n",
    "grid = griddata(points, valuesCF, (grid_x, grid_y), method='nearest')\n",
    "gridcub = griddata(points, valuesCF, (grid_x, grid_y), method='cubic')\n",
    "\n",
    "grid = griddata(points, valuesHam, (grid_x, grid_y), method='nearest')\n",
    "gridcub = griddata(points, valuesHam, (grid_x, grid_y), method='cubic')\n",
    "\n",
    "# fig, ax = plt.subplots()\n",
    "\n",
    "plt.clf()\n",
    "plt.title(\"x16 Num. Err.: log_{10}|Ham|\")\n",
    "plt.xlabel(\"x/M\")\n",
    "plt.ylabel(\"z/M\")\n",
    "\n",
    "figcub = plt.imshow(gridcub.T, extent=(pl_xmin,pl_xmax, pl_ymin,pl_ymax))\n",
    "cb = plt.colorbar(figcub)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='compare'></a>\n",
    "\n",
    "# Step 9: Compare ringdown gravitational wave signal with black hole perturbation theory  \\[Back to [top](#toc)\\]\n",
    "$$\\label{compare}$$\n",
    "\n",
    "According to black hole perturbation theory ([Berti et al](https://arxiv.org/abs/0905.2975)), the resultant black hole should ring down with dominant, spin-weight $s=-2$ spherical harmonic mode $(l=2,m=0)$ according to\n",
    "\n",
    "$$\n",
    "{}_{s=-2}\\text{Re}(\\psi_4)_{l=2,m=0} = A e^{−0.0890 t/M} \\cos(0.3737 t/M+ \\phi),\n",
    "$$\n",
    "\n",
    "where $M=1$ for these data, and $A$ and $\\phi$ are an arbitrary amplitude and phase, respectively. Here we will plot the resulting waveform at $r/M=33.13$, comparing to the expected frequency and amplitude falloff predicted by black hole perturbation theory.\n",
    "\n",
    "Notice that we find about 4.2 orders of magnitude agreement! If you are willing to invest more resources and wait much longer, you will find approximately 8.5 orders of magnitude agreement (*better* than Fig 6 of [Ruchlin et al](https://arxiv.org/pdf/1712.07658.pdf)) if you simply set `EvolOption=\"high resolution\"`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEhCAYAAADS7c8nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABg3ElEQVR4nO2dd3wVZfb/3096oST0Ii0BRLpJxIaiElx3BVEM2LsEy3dX/apEl11df+uuxr66oomu/StLURQsaIIUEVRI6FUSQOktoaVAkuf3x50Z7r25ZW7aTeC8X6/7ys3MPDNn5s6dzz3Pc55zlNYaQRAEQQgmIcE2QBAEQRBEjARBEISgI2IkCIIgBB0RI0EQBCHoiBgJgiAIQadJiZFSaqJSqkAppZVSRUqpHKVUarDtckYpNV0pVeRlXZFSqsDLugKl1MT6tc4/Sqk04/pmBtsWwbqfptewbbrxWTq/CpRSWUqpOC/bxnnYT5KxrlF912pLba5tY6cpnluTESPjwj4OZAHJwHig2Pi/MZEDxCmlEpwXGv/HAQkeHgQJQAKQ20A2+mICjuuaFmQ7TmmUUnlKqZyGOJbWWmmtFRCP4/NNAbZ4Eh6hbmjIz/dUoUmIkVIqC8fDcbjW+jmtdb7WeobWeiwwIsjmuWMKivuvyFRjXTEwzsM6tNb59WqZPVKBDByimRRsY4KJUipOKVVfopwFWL9c6/lYAGiti7XWuVrrZOAg0GS834a4PjXBh10un28Qjt/kaPRiZHgN6UCGp4e11rqw4a3yjmFPIdVFcgSQDyzzsi7oXpFxUxdrrbNxiOZ1wbUo6KRQTw8UrXW2cZ3r/VheyMfhjTcVGvr62MWjXR4+3wY9flOk0YsRjm4FtNbPBduQAMjFs2eUY7y8rQs2E4BpxvtpSFfdKYnRPddY7jlBAJqGGCXh+BXnF2NgNlMplWoEN0w0lqcafbhm4EOmW5sCt/3E1XLA1mXcyOjuitNa5+IQKud15ljSjLqy1bgGRcYrkKCIVE6OwU3HS1edYUeRMRie6Ry04e0z8GeXn3VZxssMYCky3seZxzaWJ9V2n07rzR8OOA3+e70fDDvynP43PxfnfaaaQQLOA8x2jmWcSzU77WAcL04plWB4v3lArpcfeK2cto8zhKtVAMfK9HBdq43revtsPN0/vq6P8jBQb6xP87VPp83j3K5tptu+fH0fA7XL+ZpMV05jdv7uRy/Xutb3jZ/viD973T+n2j1LtdaN+gUUANNtbjvd2L4ImIhDAMDRzZeKo1siFdBAurHO/D/BaT/pQFEtbI5zO0Y6kOe03n1dkduxa2yr0zVIMl4FZns/Nqe5n7NxHTPdlmVhiK3RRhvHiXM7vvtn4NUufzYb67WxvwQcYx0ax0M11WiTBxTYvQ4+9pngdk20zc/cvBZxTv8XADlO22Sa94Fx/On+jmXHTh82pRvbur+yTDttbOv8SrX5HTTtTXKyd6Kdz8bH/ePr+kx3W6aBNBvPBdPWdOPaTsTp++bv+xiIXcY1L3DaV5Zhj7stAX3Otblv/HwOdux1uabGPjSQ5O355PNc7GwUzBeOh0yOh+Vx5svDB5DkZ5/TcX1IuDx0jWNm1dJuS0SN4znvP8dtnVexDcRW46ZxvxnScBJCH8fJcT9n84Z0W1aE6xfd/f9qn4Evu+zYbOyzyOn/OKNNltOydPNLWct9prq10TY/b7N9mtO1y3Rub9wTmZ4+d2/HsmOnD5vS3feJ64Ml1X1b3ETKWJcUwDFd7lenZQV2PhtP94+N62NHjDzt05Ot1e55X23s2OV0/RLctnH+7tboc67pfePrcwjAXk/X1OXZ6f6/r1dT6KZbhttAq+FOFpkv5Rqimq89BDoox/wZsyvFfSwk21xm7CuJ2oeMz+Dk2JB7/3yOcQyMvy5997Ww1Tye2a2gcdw0PqPi1MkxBJd5KRi/GJX/qDr3bhz3z8CXXXZtXma+0VoXm22c1h+0eTx/+4yrdnY2MNrnczI4ZRyOz6XY6MqIw3EfT63B7uvSzkKttTk22FAD3zmc/A7b+Ww8fodrid19OtsK+P0+2iEFR2CQe7BVrrHOpM4+Zxv78/U52LXX0zXNwogWdno+2brPmoIYmWMXVp+jdoR3KzyHdVeLrjNuInOO0nCM8Rknsjj50E3H8UHU9stgjhs5jxeZ5Dodz2V+UR3YWqyNeSXOLz+2pvpp5xxVlwtMUCdDSuM4GfRg4inC0df+7dhc7GGfBz0ss3M8X/usDVOBVGWMBxpf5lwc9+k4an5fFdeZhSfJ4uT92dD4+2zqI0LW7j5dfljZ+D7aIc7muuIa7NsX/vbn7XOI89HGeZ2nazoNx32Vxsl73lakcKMXI31y0L/arHE7GF+2JGCsdsyxKPZwDDMcOxXHg7fWIZlOH8DjuAVgOD2QHgcKzV8gdWCrGRwR6APGOYrOnRk4RM+dIhzdUCM82emGL7tqanNNj1efzMDx4yKNkz8wpuL4rBpF+L4TDTm3bSwnvwP18dnEmW9q8oxww5yCYev7aBNv52zOPQwGNflO+rXXuEYzcDybRhDAs7TRi5HBWBwqn6ccaUuSDE9prI22pnpnqJPRRN4mqV2Hhy46o52nB7I/8nF9MDmT62FdrWw1hCobmG50DSWYXQzeDHTqovO2zVRcb0yzW3EEDhHz5Z34tasmNtfmeAHuyvyRkGbsx2dEkHHcYhzXxeyOy+XkALGvLrqAjlVTjGuRieOHREZ9HAOHd5jmdKxU4Bmo1Wfj7focxMhqYtyjgX7GzrZOxPF9e8b5mPj+Pvr93AzBn+F0zknG+SZQ+8+gRveNn+9kbe3NwvGZpxJIt7SdgaXG8sLxBSrAMXBWhOOh6D5QWS0YgJORJEWcnBntPugZZ2xTbfASDwPBAdjrcRCSk5E7aXVpq9O+zeuU536MQM/N2I858J5D9SirAozBTm+fgT+7/KyzM0hdbSC3tvs0luU5tbczgJ+FWyCA0z7i/By/2rHs2unrs3V7md+bVC/bxnnYT6ABDHmcDJLwGM3p7bPxc/94uj5JxnE0J8diC5xt9bZPTn6/vNqKve+j3c/NfH6Zn4F7ZFtNP+ca3zfePoea2Ou23wJ8BIJ4eimjoSD4xfhlOB1I1Ce7FhOMZYXakZ5JOI0xfkHHaa0bW5ouoQFRjjlQ+Vpr255fU+mmExoHCThuMGvg0nifhZ+IPUEQTg+cuv4DikgWMRICYQaQpBwzr81++jQc7nxNoowEQThFMJ4HCTgCs3J1gHlDw+rHLOFURGtdqJRKxsgkgGPsKh9HEtuGSAopCELjJQXH2JKnRNF+kTEjQRAEIehIN50gCIIQdE6rbro2bdro7t27B9sMQRCEJkNeXt5+rXXb+j7OaSVG3bt3Z9myZf43FARBEABQSm1riONIN50gCIIQdESMBEEQhKBzWonRunXreOCBB9i1a1ewTREEQRCcaJKh3UZCRDMh4DnAeG0jo27Lli11SUkJrVu35quvviIpSZIGCIIg+EIplae1TvG/Ze1ocp6RkWoiRTtqGj2HIyvsXDtte/XqxfLly4mMjOTuu++mKQqxIAjCqUhTjKZLwZHG3Jzxn4sj1XmcHe+of//+zJkzh9jYWJTyV3NOaOocPHCAHVt/4bgOA/m8BeEkWhOhKujcvRetWrcOtjVNT4y01rlKKefs0GZVzWK7+zjrrLPMfVFcXEx8fHyd2ig0Dg4eOMBvBetIXJJBTPFGQnRFsE0ShEZDlQqjJO5MNlc8S3lZDzp27hJUe5pcNx1Uq055HfCct22NYnzLlFLL9u3b57Lu97//PTfeeGM9WSkEmx1bfyFxSQbNitaKEAmCGyG6gmZFa+n542P8tnk96/IXB9eeoB69lhjjR0m+amZorbO11ila65S2bV0nEQ8dOpQ5c+ZQWOg5uezChQsZOnQozZo1o3///rz55ptUVVXV5SkI9chxHUZM8cZgmyEIjZqY4o2ExsbzzX/fonD9yqDZ0WjEyPBgMn28PJXTzcRe6XGP3HzzzQBMn169UnFpaSk33HAD27dv56677qJFixbce++9vPXWWzU9nNDQKCUekSD4IURXgAohMjqGgjX5/hvUE41mzCjQEgRGvfoMrXWx3eAFd7p3786QIUOYPn06GRmuzlV0dDTffPMNZ5xxBnFxcWitmTp1Ktdee22ghxEEQWj0hEVEcPRwUdCO32g8o0AwCrrNcBIgT16TLa688kry8/M5dOiQtWzLli1UVVXRv39/4uLizGNy/fXXEx4ezv79+/nll19qfgKC0ITJzjtOcdmpPy0if1dlsE1woSHsCeZ0l0bjGdnFqCQ43XhvLi6khpVGx40bR9++fQkPDwegoqKCCy64gCuvvJK333672vZaa4YNG0arVq34/vvva3JIQfBLcZlm/OxScgsd3YwpnULJTI0iqWMoACM+PMayndUfTo8PjWTihZGM+PAYuYWV5NwSQ2rCya/5hNmlJLYKYeKFkS77SYgPsQQmtUcYmSOiiIvyHAqft7OS9OQIq31dHcfTNmlnhZM5IqqGV7F2JGcfIy891rrmdUF85mHy0puREB+4H1Af9jQmmpwYGaVs62zCSJ8+fejTp4/1/9y5c9m9ezejRo3yuL1SigkTJvDAAw/w888/M2TIkLoyRRAskrOPMiE5grdGRRMXpZix7gQHS11/tZrC442kjiFk5JaRl97M57Hc9/PcD+UM/+CYx3aFRVUkd3J9GNblcdy3GTu9hAmzS8kaFe1z3+4Ul2mvYmoX/WSLWrWvDRk5ZSS2CrFEP9j2NARNspuurlm1ahWfffYZAN9//z2hoaGkpnrv+bvjjjto0aIFr7/+usf1e/fu5dprr2Xr1q2Ao9tvxYoVdWy1cKpSWFTFwVLNxAsjrQdqWt9wF8/DDhOSIygsqrK8K7tMvDCS/F1VFBZVjxydse4E4/qF1/txTN4aFU12/omA9gvwzPflAbcRgouIEfD0008zceJEAJYsWcLAgQOJjY31un3z5s255pprmD17NidOuH5RDh06xIUXXsjXX3/Nxo2OsOInnniCCy64gNzc3Po7CeGUwdFN5Xjw15bM1CgyfwjswfzcD+UkxCuPXUkFB6s8ehx1fRwTZ29wwuxSMnLKrP+LyzTxmYcDOqYniss0Iz48RuKrR4jPPGwdIz7zsCWUia8eITvvOCM+PEZ85mFmrDtBdt5xkrOPWv+bOLcDUE95tjEjp8w65oTZpdbysdNLyM4/TkauY725b+f9FhZVWTaP+PCYy/FMW03bsvOO1/oaNQRNrpuuPujSpQuzZ8+moqKCn3/+mdtuu81vm6uvvpr333+fZcuWcf7551vLH3roIbZs2cK8efO46KKLAHjxxRdZuXIlY8eOJT8/nx49etTbuQi+eXBOGSt2N+zA9OAOobxyRWDjHtPHRjN2eilxUaWM6xvOhJSIamMFzywqJ8vtQZOX3sxFLNKTI8jILSN/V6XXsQbn/RQWadKTwj12ueXvqmREoudHRl0ex/l442eXMvECR1fVhJQIhn9wzBpDmrb2RMDeoiemrT1BQlwIObfEGrZV99QOlmrydlaSc0ssuYUVjPiwhIkXRJCX3ozsvOM8s6ictL7h1dr5YkRimHUupuik9Q1n+tgYJswuJblTqEs3nUvbD48xfWwMSR1DDXuOUfCn5i625qU3I7ewgrHTS7zupzEhnhHQtWtXysrK2L9/P9OnT2f8+PF+2/zud7+jsLDQRYi2bdvG+++/zwMPPGAJ0cKFC3n++ef5wx/+QGVlJffee2+9nYdw6pDWNxz9ZAveMsZKkrOPVfOUHh8aScGfmru8PHktjw+NJCO3rNpy9/04hAwyhkZ63M/UNSd8PnDr4jjPLCon8dUjJL56hIzcMh4fGmk9sJM6htIqWlndgVl5x5lQRw/Z3C0V1n69eWpjje5JUwCv63/y//xdgU+GdxbStLPCWbrD3o+k7LzjJHUMtUQ/NSGMpI6hLveHs63F3j+SRoV4Rjg8I4CdO3dyxRVX+N3+559/5qmnnqKiooJJkyZx8cUXA5CdnY1SigceeACAf/7zn0yaNInIyEjKy8tp3rw533zzDbm5uT7HpIT6I1APJdik9Q0nrW84ia0cQQKB/voGx9jMM4vK/YYGx0UpHh8ayYQvSi0voaGPk54U4TN6bkJyBNPXniClUyiFRVXGw1ZXGyPK3VLh0qXXOkZ5DfZIT46g4GAVY6eXUFzm8Eo9XWd3kapJRJwz+bsqeWZROcVl2nEuPew9jgsOVpEQ53rsVlGKpTsqLbtra1swaHoW1wOmGH3//fd88cUXlJSUeN129erVDB8+nOXLl7N8+XIuvfRSlixZAsCTTz7JypUr6dq1K7m5uUyaNIkbbriBoqIiVqxYQcuWLWnevLk1d8kbhw8fZvXq1Rw+XPv+cOHUIDUhrFo0XSA8PtQhFP6YeGEky3ZWVhOU3MIKr110dXkcf6QnRzBt3QmmrT3BOOPBGxelyBwR5fIyw8bNl6+oQ4DMEVEUZbQga2QU453Gb3xhJ1rP23yswqIqR5djahQ5t8SSdlZ18fPWNrFVCIXFrp5YYXEVia2a9uO8aVtfR/Tr14+ff/6Z8vJyRo0a5VUEtNbcd999xMbGsnTpUv76179SVVXFfffdh9aaiIgI+vXrh9aaxx57jG7duvHuu+8SHR3NoEGDmDlzJiUlJXz44Yc+7Tn//PMZOHAgbdq04U9/+pNPcRROPcx+fvNBXVymycgtIz2p5l1S6ckR5BZWsMzGwz8ztfoDOaegwtb4TG2P44+4KEVKp1BHF11KHXXRFVZY40R1MQbVKlpZ+/MWPFBYVEWr6JPBG/lu45hxUYqCgyeDFZwZ1y+c3MKT3Yq5hRUsc5r71VQRMQJiYmI455xzKC11fDFae6ntMX/+fBYtWsRf//pXOnfuzLBhwwBYsWIFL730Eg8//DCHDx/mxx9/JC8vj7/85S9ERp78RZaSksJtt93G5MmTPebD01rz1FNPsWHDBgBOnDjBa6+9xvDhw0WQTiNSE8IYkRDG+NmlqKcOk5x9lKQOodW6rzJyy1FPHXZ5jfjwmMd9xkUp0pMibI1tpBuh2jWJ5muI40xIjqC4TNfZ5M/iMs3Y6SVWZNpbAc5p8mTf2OkljPjwGHFRioT46h5UakIYCfEhxGcetrZz9myu6x9Odv4JkrOPVvMe46IUeenNyMgtc0T/2Zjj1RRokmXHa0pKSopetmyZx3WzZ8/mjjvuoKqqioMHD3rc5o477uCzzz5j165dREVFcfjwYVq2bEloaCgDBw5kw4YNHD58mMcff5x//etf7Nu3j5YtW7rsY/PmzfTq1YuePXu6pBQyRW7+/PnceOONPPDAA+Tl5fHwww9TVlbG1KlTGTu2xjlhT0vy8vJInn1ZsM1o8sxYd4KE+JBGM/N/xroTFBZV+ex6e+6Hcr9dc8JJ8kZ9R96cj2jT/gzG3P2wyzopO97AZGVlceDAAdzLTJhUVFQwe/ZsRo4cSVSU4xdq8+bNiY6OpmvXrmzcuJFu3boRFhbGzJkzGT58eDUhAujZsycdO3akoKCAo0ePWsufeOIJ5s+fzxVXXMFHH33EkCFDuPfee1mwYAERERG89957UiZdCAo5BRWNRogAW2HUIkRNDxEjAzMlkLcuusWLF3PgwAGuueYaa5lSil69etGjRw9KSkpo3bo1+/bto6CgwGe03JVXXonWmq+++gqArVu3Mm/ePJRSvPHGGy7l0M855xyeffZZvvrqK55//nkRJKFBKS7TjWZg/LkfHGHfqT3CmmS0mOAb+UQNTDF69NFHPa7/+eefAawwbpOVK1fyt7/9DXCUnVi7di0AAwcO9HqstLQ0AEuMvvzySwAuueQSunfvXm37+++/n3bt2pGRkYG3bkaA7du38+KLL/L+++9TXi7pUITaExflPSS6oZl4oWOeUrASpwr1i4iRgSlG0dGeBy/z8/Pp2rUrbdq0qbbuzDPPBBye0po1awDo37+/12NddpljHMMUFjOYYcKECR63Dw8Pt+otPfvssx63GT9+PF26dOGRRx7h9ttvp3nz5rz77rtebRAEQWhMiBgZmIKyYsUKtNbVyovn5+eTlJRUrd17771Heno6bdq0oVu3bqxZs4ZWrVrRoUMHr8cKDw/nvPPOs8aezFByMzrPE3fccQchISF89dVXVFa6Rte8/vrrvP3227Ro0YKnn36au+++m6qqKh566CE2b95s7wIIgiAEEREjg7Zt2/Lmm2/Stm1bzjjjDNq1a8fcuXMBOHLkCJs2bfIoRr/++iuff/45vXv3ZsOGDaxZs4b+/fu7jPt44qKLLmLVqlVUVVXRvXt3EhISfApYfHw8SUlJlJWVkZeXZy3fuXMnjz32GMOHD2fz5s1MmjSJt956izVr1hAWFsbVV18tYeGCIDR6RIycSEpK4p577qFt27Z06NCBa6+9ltLSUgoLC9Fa07dv32ptTAHZtGkT69atY/fu3Zxxxhl+j1VRUcGJEyfYtm0bS5Ys4cILL/TbxhxrmjZtmrVs0qRJlJWVkZWV5RIJ2KdPH95//33Wrl3LwIEDJfBBEIRGjYiRwYkTJ7jjjjvo0KED8+bN45VXXuHQoUN8+eWX7N69G8Cj59KpUycA9u/fz5EjRygtLfU67uRMWZkjb9ZLL73E7t27adeund82V155JQC9e/cGHII2ZcoUmjVrRmJiosft+/XrR0FBgRUkIQiC0BgRMTL44IMPWLt2La+99hrx8fFceumldOjQgSlTprBnzx4A2rdvX63diBEjeOWVV7jvvvs4ceIER48etSVGXbt2BRyTXQESEhL8tunbty+tWrXip59+AmDOnDmUl5e7ZA5357XXXgOwAiA8MW/ePG6//XbuvvtuPvnkE/GiBEFocESMgMrKSp5++mmGDBnC6NGjAQgNDeXiiy9m1apVlhh58owiIyN54IEHrBpFhw8fJiYmxu8xTU/GzMLgKaTbnZCQEPr06cOcOXM4dOgQH3/8MYDPshSXXnop7du3t7oQ3XnkkUe47LLLmDJlCjNnziQtLY127dpRWFjo1x5BaEwEmnRVaFyIGOGY57N161YmTpzoEnjQrVs3fvvtN3bt2kVMTAzNmnnP/+TsDdnxjHr27AnAsWOOXGK+ghecadGiBTt37mTdunXk5+cDDsHxxQ033AA4skw4M23aNF588UXCwsKYMmUK+/bt48EHH2T//v0MHDjQa1okoX4ZO72E+ExHTjqrtk+Oa1GaER8e4zkPlVW9LW+suFdFrQ3J2ccaVJDq0nZBxAiAf//735xxxhmWV2TStWtXysvLWb16tccuOmecvSE7YuQ+RmRXjAYMGAA4Qs1//fVXvyIJJz0nc5wKHGNk6enpAEyZMoUxY8YQEhLCyy+/zMMPP8yxY8es+VBCw/P40Ejy0ptZxejM0tiNiYycsqCVtPZ0bP1ki3pLWxTMcz1dOO3FaNOmTeTk5DBhwgTCwlzTx3fr1g1wZF/wJxaBekadOnUiJeVk7kE7AQwAycnJgGOcp7S0lBtvvNFvm969e9O5c2d+++03a9mnn37KoUOHGDx4sBWlZ/LCCy9w5plnsnLlSpfIPSE4mPV6cgulG0o4dTntxejNN98kPDzcY6lxM8jg8OHDfj2jQMVIKWWJXdu2basJoTfM7r2NGzcCVBMSbyQnJzNv3jxLkF5//XVCQ0P5y1/+4nH7GTNmAPDUU0/Z2r9QfxSXacbPLiU9KfAqr55Izj7Kcz+Uk5x9lPjMwy7deoVFVVb34ASnWkOJrx4hf1clY6eX8NwP5YydXkJ2/nEycstIfPWIVQbCvetKPXXY6z5MspYd92hLRo5j3/GZh11ssXPswqIqRnx4zCoL4WxT4qtHyM47eUx/Ho+344Ejg7in/Xi7jv7scr4+E2aXunTPFpdp4jNP3YKbp3XZ8YqKCj7++GNGjhzpUWxMMQLPkXTOBNpNByfHi/bt22dre8AKlNi5cyeA38m1JmeeeSazZs1i5syZ3HvvveTl5XHPPfdw9dVXe9y+f//+XHbZZaxatYry8nKXukwm+fn5PPLII6xbt46ePXvyyiuvuHh7jZVL3qve3TWuXzj3nRNByQnNH/6v+iTh2weHc/vgCPaXVJE2rXpBuHtTIriufzi/Harilpmu6+ffHngJ72cWlZOVd5yDpZriMshMjayWI87cxpnCIs0IPwXizAdgXnoz8ndVkpx9jNSEMJI6hpKcfZTpY2NITQhjwuxSsvOOk54cwcFSR82fzNQoK2P2hNmlJHcKtV3UzdM+wFESPC+9GYVFVSS+etSyZURimJWHzhSBtL7hTB8b4/fYIz48xvSxMSR1DHVUqf3wGAV/am7Zkbezkrz0ZlYhQ1/n4Ot4BQerPO7H23X0Z5fz9cnfVemoBmtcg2lrT9RJ8b/GymntGX333Xfs2bOHm2++2eN65/Lgpgh4I1DPCByh2YESHx9P8+bNOXHC8ets27ZtttqZWcQXL17M8uXLKSkp4eKLLyY01Hsf+yOPPML+/fv59ttvq62bPXs255xzDvPmzWPPnj3k5eVx7rnn8uabbwZ8TkJ1Hh/qSApalNGCvPRYcoyHnadtnF+pCfbGTEwxSOoYSlrfMKauOcGMdSdI6RRqPfAyhkaSY1QTBci4MNJv6QZ/eNqH+X9CfAjpSeFMXeO4t50fvGlnhbN0h71uyuy84yR1DLXGj0xxc/ZoxvYLt9YVl3ncjS087cfbdbRjl/P1SeoYSqtoZVV0zco7zoQmXs3VF01eZpVSWVprzxlG/fDRRx8RFxfHH/7wB2/75t5776W8vJz777/f576cPSM7od3OLF++3Pa2Sinat29v5ZwzJ93646yzzgIcY2RmmiN/Htnw4cOJjIzk7rvvZteuXYSEOH677Nixg3HjxlFVVcWECRN44YUXKC8v5/bbb+fee+9lxYoVjVqUfHkqMeHK5/o2MSE+13dp6Xt9TUjqGErOLbEkvnqE3EJ75b8D4ZxOoRQcrKKwqIplOytJzj5ZZyvFKSAgpVPtgwP87SOxVYhVbjt/VyXPLCqnuExTWFRFag97511wsIqEONff2a2iFEt3VLoIX13gaT/erqMdu9yvz4TkCKavdYhbYVHVKe0ZNekzU0olAelAwGJ07NgxZs6cyfXXX28lLPXE5MmTbe2vJp5R3759WbduHYMHD7a1vYlzEle7YmQGYOzYsYPvv//eVtuIiAj69+9PXl4eq1evZtCgQYCjEGBZWRkjRoyw6i81a9aM6dOn07FjR7Kyshg3bpxE49UxCfEhdRZKXFx2cmLz0p2VnNMplLgoxbi+4WR5KbsdF1W9S9h5P3bW+dtHwcEqEls5znP4B8fIS29GQnwIGTll1fbn7diJrUJcPDqAwuIqxtbSq/N1rs54u47Zecf92uV+fdKTI+jxryMkdwplXC3tb+w09W66BKC4Jg1nzZrF0aNHvXbRBUpNxoyWLVvGoUOHAj5WaenJMQlPJS08ER4eTnR0NAcOHGDLli0Alrj4YsyYMQB89tlnABQXFzNlyhSuuOIK3nnnHZcxq6ioKKtG00033SSZHOqQ7LzjLNtZybh+dfNAMrt+8ndVMmNdBWl9wxnXL5xp605Yc3WKy7S1nSfiopTlxZgi2SpaWe/thkI725Kd7xgXKiyqolW0sjyP/N2uXXSejm0yrl84uYUV1n5zCytYtrPS9tiWJ3wdzx1v17EmdsVFKVI6hTq66FJO3S46aMJipJRK01rPqGn7//u//6NLly5cdNFFdWJPTTyj6OhoWrRoEfCxYmNPdgMF0r579+6cf/757N27l9DQUFtZH8yxJrNrLycnh9LSUv7yl794TAh7/vnnc8EFF7B7924+//xz27YJrjyzyBHtZkZk5RRWkJfezKNnURMKDjqivYZ/cIzpY6NJiA8hLkox99ZYxs8uJfHVIwz/wPe8puv6h5Od74gmMx+8E5IjGDu9hBEfHiMuSpEQ79vehPgQDpRokrOPMnZ6iWVLaoKjmmt85mFrX84VZz0d2yQuyhEQkZFbRnzmYTJyy8hL9z0Xzx++jueOt+tYU7smJEdQXKYbVen3+kA1xV+vSqkEIE5rna+UKtJax/vYNh1HVx5du3ZN3rZtG/v27aNjx4488sgjXovVBcrRo0dp3twRFbNp0yZ69epVJ/v1RFJSEsuXL6dr165s2bLFGsvxx5VXXsnu3btZs2YNUVFRtryy0tJSYmJiaNmyJcXFxYwdO5ZZs2axdu1aK8zcnZ9//plzzz2Xvn37WpVvnTl48CD33nsvX3zxBQBDhgzh448/pmPHjrbOwy55eXkkz5auQnfiMw9b3V9C42fGuhMUFlXVa8XdvFHfkTfnI9q0P4Mxdz/ssk4plae1rvcw2aZ6NyZprfPtbKi1ztZap2itU8wSCzNmzKCystLWhFG71MQzqiktW7YEoGPHjraFCBzRgevXr+f48eM+K9E6Ex0dTadOnWjWrBlaa+bOncvx48d9HnfIkCEkJCRYoevOFBcXk5SUxLRp0ygrK6OyspL58+fTtWtXq7S7IAgneWZRea2jGJsCjUaMlFLpSqlMH69UY7tUILc2x/rvf/9L3759rdQ6dUFoaCjh4Y4bpr7FyAw5dx47skNkZKTV5q677rLdbtiwYYSGhlJcXExRURExMTF+Q93Hjx/Ptm3b2LVrl8vy+++/n+3btzNmzBh+++03ysrKePXVVwkJCeG2224L+JwE4VTluR/KSXz1CKk9wk4LL7bRnKHhwWT4eDkL0DhDvNKBOOO9/xoMwPbt2/n++++54YYbbE8YtYsZxBBoaHegXH/99QBs3bo1oHadO3e23tuNwgNHWqSdO3dSUFAAONIL+bt2v/vd7wB4+eWXrWVr1qxhypQpTJw4kU8++cSy4Y9//CNffvklGzZsYNKkSbbtEmpGUUaL0+Lh1tSZeKFjHpk56fVUp8ndkVrrXEO4srXW2caybK21rZoH06ZNQ2vNddddV+e2mR6Rr1DxusCMoAv0OOZcI3AVCX+Eh4dTUVFh1UY6++yz/bYZNGgQkZGRvP7669ayzMxMAC655JJq26emppKWlsbLL7/MBx98YNs2QRBODZqcGJkopeKUUhON95l2PaMpU6aQnJxcLwEGMTExREVF1bnH5c7+/fsBrCwMdrn44out957S+3jDzKFnlqzwVH7dnZCQEBITEykpKWHHjh1orfn888/RWluVat155plnAHj88cd97ltrTUlJ9ZQ9HjakSjXpqXSCUO9UqTDQwS+F0WTFSGtdrLV+TmutjG48v55ReXk5y5Yts+r71DXR0dH1Pl4Ejpx6ELgYOYeEx8d7DUCshplRfNeuXcTFxfHII4/Yamd6UGYKoiNHjtC9e3evIeU9e/Zk0KBB7Ny502MwQ1VVFS+88ALx8fHExsbSvHlzHnvsMa/Hj1AVlMSdactWQThdKYk7k6qy4CdgbbJiVBPMYnHjxo2rl/3HxMQ0iBiZnkqg3XQtW7ZkzZo1hIWFBTSu1bp1awAOHDhga26SidkdN3fuXObNmwecnLfkjT//+c8A/POf/6y27tFHH+XRRx+lf//+dOrUibKyMjIzM0lJSXHJSmHSuXsvfhnyD47G9xMPSRDcqFJhHI3vx6Zznmb3r4VUVVYRHkCPSV1zWn1DDx48yEUXXUSXLl3qZf8N5RmZod3mg9suISEh9OvXj8jIyIDEyDnLw+7du9m8ebPXOUbOmJ7R+vXr2bt3LwCXX365zzZjxowhNDSU+fPnuyyfN28eL730Evfeey+vv/46SilKSko499xzycvL4+qrr2bWrFkubVq1bs2voeGsHfw3wmLjQZ1Wv70EwTe6iqqyw+zeWsDh/TspOXKIDl1sjXbUC6eVGJWVlVmRaPVBQ4mRObk2kK42k6ysLDp37swFF1xgu43pGYFDjCor7WVPNgMmzj33XHJzc2nVqhVDhw712SYsLIyUlBS2bt1KVVWVNZ/pvvvuAxx1mcwxuZiYGJYvX06PHj3Izc1l586d1aIEBySdR86Md1i77CMio2MJ8ZGlXBBOR3RVJWUlJfTsn8TgC4YHzY7TSozAfjG6mjB+/HirK7A+McuMB1IHyeRvf/sbV111Fddee63tNs6lNMC+CJqZG8rLy9m7dy8jR460lWXhhhtu4MEHH2T79u107dqV9evXs2HDBpo3b14tl2BYWBjz5s3jrLPO4tlnn+XVV191WR8aFsaItDvp0rMv+3b9RkV5OYIgnCQ0PJw2Hc6gz9nnER4h3XQNQosWLWyX964JgTzga0ObNm2YMWMGl156acBto6OjOXLkSEBt3GseBeKRRUVFMX/+fHbs2MGZZ9oLJjC9pw8//JBJkybx4YcfAnDHHXd4jALs2bMn1113HW+88QYjRoxg1KhRrvaHhdEvxbHPoqIioqOj6z38XhCEANFanzav7t2769OduLg4Dei33347oHaA9QqE9u3bW+369Oljq82RI0c0oLt166a11rpnz54a0AsXLvTa5ueff9aAPvvssz2uX7Jkie7Tp49lS8+ePXV+fn5A5yIIpyPAMt0Az+fTakTXvbvpdMQcb2koz8Ac3wL7179Zs2bExsayY8cOjh8/ztatWwkPD+fcc8/12iYlJYXWrVuzYsUKiouLXdYtXbqU1NRUjh07RkpKCj169GDz5s2kpKQwc+bMmpyWIAh1zGklRr5KbJ9uBJqy6MEHHwSoNibjD+cuvUCycnfp0oWKigq+//57KioqePHFF4mI8F7PRSnF1VdfjdaaTz/91Fp+4sQJxo0bR4sWLfjpp59YunQphYWFVn2msWPHVsufJwhCw3NaiZEA77zzDhC4GL3wwguEhoaye/fugNo5R+IFElKfkOAIMTUnv55//vl+25jBDVOmTLGW/d///R9bt26lXbt2LmI4evRosrOzqaysJDs727ZdgiDUDyJGpxlmJJ5zNgY7hIaGUllZyUsvvRRQO+dQa3Oyrh3MYAdzvpE5T8kX559/PiEhIWzfvt1aZtr7xBNPVNv+rrvuYty4cWRmZtravyAI9YeI0WnGihUrAOjatWvAbadNm2ZlUrCLOUE3IiLClndjYua/++233wBsFQKMjIzkwgsvtKrf7ty5k9WrVxMVFcVVV13lsc0TTzxBaWmpX9uKioooKCjwmOlBEITaI2J0mlFYWEibNm1qJEZjx47lvPPOC6iNKQwdO3YMSIzMdocPO3Jm2fWqzj77bNauXYvWmu+//x5wJIgNC/M8i6Ffv3707t2bwsJC1q1bV2390aNHGT16NG3atKFnz560atWKG2+80fbEX0EQ7CFidJoREhLC/v370Q1Ubt6MpmvVqlVA7cxgBXMSsV3xbN++PceOHePVV1/l22+/BfznIjSzO5glLkwqKyu54IILmDVrljXGFBUVxZQpU+jbt2/AiWoFQfCOiNFpxscffwzULHtDTTA9nEBTF5liVFpailLKdiRev379AFi0aBHbtm1j8ODBfsvL33rrrSilmDNnjsvy5557jtWrV9OpUyf+85//MH78eLZv3865557Lpk2bqk2uFQSh5ogYnWZcc801QOABDDWlpmJklnAHR9YIu2H5ZhTe+vXrycvL49xzz/WbLzA+Pp6uXbuyd+9eK5DhxIkT/OMf/wDg008/tewPCwtj4cKFtG/fnm+++YbZs2cHdF6CIHhGxOg048033+S3335rcDGqaTcdENA4VY8ePQBHAcLi4mJ27Nhhq92wYcMIDQ21ivZ99tlnHDt2jOTk5GqTbSMiIpg+fTpAwNGFgiB4RsToNCM8PJwzzjijwY5X2246CKwqbbNmzYiMjLSi73755Rdb7UaNGkVlZaXVffn++++jlOJ//ud/PG5/0UUXccUVV7Bu3TrKfSRfLS4uZu7cuaxfv972OQjC6YiIkVCv1IUYeYuE80b79u2tAA27E21N7+fLL78EYOXKlYwZM8ZnyZEHH3yQvXv38tFHH1Vbp7XmySefpHXr1qSmptK3b19+//vfu8yBEgThJCJGQr1iZmDo0KFDQO2cxSjQPHoDBw60vBWz284fZ5xxBqGhofzrX//it99+Y/v27QwbNsznsVNTU4mNjeWee+7h6NGjLuueffZZ/t//+3+Ao2Dgn//8ZxYtWsSgQYN49913AzofQTgdEDES6pVOnToxd+7cgIsaOgcw+MpJ5wnnFER2Q8KVUsTFxVFcXGxlffAXcRgaGspFF11ERUUF33zzjbV8y5Yt/OUvfwHgm2++4ZNPPuEf//gHc+fOpbi4mLvvvpvVq1cHdE6CcKojYiTUO5dddlnA3o2zADkLkx2cRSSQyb2m97Zw4ULA3ljV7bffDsCMGTOsZU8//TRVVVXcfPPNpKamWsuHDBnCK6+8QlVVFWPGjJFsDoLghIiR0ChxFqNAPSMzBRGcDGW3Q/fu3QHIy8sDHN19/rjooosAWLx4MQDl5eVMmzaNUaNG8cILL1Tb/n/+539ITExk8+bNVuZwQRBEjIRGSm3EqH379oCj681MDGsH04v69ddfgZMTaH3RsWNHYmJi+O233zh69ChLlizh6NGjjB8/3rLDGaUU//73vwH485//7HW/x44dY86cOXzxxRe28vIJQlPntCo7LjQdajNmZHa3mRUkzYKC/ujVqxcAx48fB+zVX1JKMXjwYEvAZs2aZY0/eePyyy+nQ4cOhIeHe7Rv7ty5XHPNNVZ5eLOw4Jw5cxpsfpggNDTiGQmNktp4Rs7zqEJC7N/iPXv2BODIkSPExcXZnt80bNgwdu/eTXh4OHPmzEFrTUVFhdftQ0JCuOeee1izZk21IIk1a9bw+9//niNHjpCZmcm8efPo1q0bixYtok+fPhw7dsz2+QhCU6LJipFSKk0pla6USlVKpfpvITQlahPAkJiYWKNjOtdeuuWWW2wLWa9evaioqGDx4sVs3boVgMGDB/ts84c//AGAxx57zGX57bffzokTJ7jxxht59NFHueSSS9i4cSNDhw5l+/btXH311bbPRxCaEk1SjJRSaUCC1jobKAQy/TQRmhjOE10D9YxMMXIOZLCDWdAP8Dje4422bdsCjtRApaWltGrVyu8k36SkJKKionjvvfes7riFCxeSl5dHbGwsb7zxhtV9FxISwpw5c4iPjyc3N5evv/46oPMShKZAkxQjIFNr/RyA1rpQa50cbIOEukUpZXlEgYqRKQSBllZ3DnYws5vbwRSjDRs2ACeTtfoiNDSU5ORktNYsXboUgP/85z8ApKenW5krTGJjY3nzzTcB+Pzzz23bJghNhSYnRkqpJKDY6b1wimKKUKBiFGj6IGfMrrlASmyYYrRr1y7AnhjBybDwhQsXUlVVxVdffUWHDh1IT0/3uH1aWhoDBgzgp59+8rrPsrIysrOzufHGG5k0aRLFxcW2z0MQgonPb61S6lEgsHTLJzmgta4+0aL2JAAHja66XKXURCBfa51bD8cSgkhERATHjh0LWIwAZs6caUXHBUJ0dDTHjh0LyKsyxcgMLnj55ZdttUtOdjj0CxYsYNWqVezfv5/33nuPPn36eNw+JCSEu+++mwceeIBPPvmEa6+91mX9b7/9xvnnn++SqTwzM5PMzEwefvhh2+cjCEHBDH/19AJa1ubla981fQHpDrOt/+OAIj/bLwOWde3aVQtNh/bt22tAv/XWWw12zDPOOEMDevDgwbbbVFVV6bCwMA3o8PBwXVlZaavd2rVrNaC7du2q33//fQ3oDRs2+GyzZ88erZTSkZGRury83Fp+/Phx3a9fPx0eHq67d++uZ82apZ966ikdHh6uY2NjdWFhoe3zEQRngGW6Hp7l7i+f3XRa60O1eQUiikZkXKaPlxkxdxDId7KxGIhTSnnsG9FaZ2utU7TWKeYvWKFpYI4ZBRpNVxtef/11WrZs6bcgnzNKKdq0aQM4xoK++uorW+169uxJaGgo48aNIz/fcUvPnTvXZ5t27drRp08fysvLrRx6AJMnT2bt2rXMmDGDDRs2MGrUKJ544gny8/MJDw/npptukvRDQqOm0Ux61Y7IODvke1hWjEOkhFOImo4Z1YarrrqKO++8k8rKyoDaderUid27d1NWVsaePXtstYmIiKBnz54UFhZaFWZ79+7tt92oUaNYv3493377LZdffjlVVVU888wznHPOOYwaNcplEm3//v159dVXufXWWznvvPP46aefbE8CFoSGpM4CGJRSLfxvVXu01oU4xozijOMmAIWGhyScQgRDjPLy8rjpppv417/+FVA758mo7pFwvoiJiWHZsmVs2rQJsFdqY8SIEQB89913ACxatIg9e/awadMmj5Ntb775Zrp3787SpUuZOXOmbdsEoSGpkRgppVq4v4BxxvtH6thGT4wFHldKpQMTgOENcEyhgQmGGD344IM8+uijAbfbuHGj9b558+a224WGhvLrr79a0Xt2upLNMuwbNmxAa82UKVMAGDt2rMcuTaUUzz77LABPPfWUx30eP36crKwsLrnkEoYOHcqsWbNsn4Mg1AU19Yz+DOQBbwFvG697gBk4xKFe0VoXa60zjPGgDPGKTk1qOs+oNkRFRTFv3ryABWncuHHW+0A8I3OsSRuVaZ1rMXmjWbNmtG3blksvvRStNQsWLABOekyeGDNmDNHR0axatYrNmze7rCspKWH48OHcc889LFiwgMWLFzN69GjuvPNOTpw4YftcBKE21EiMtNaPAZcD2VrrcVrrccAzWuvLcYiSINQaU4QaMoDBrLt08GBgQ5DPP/+89b5VK/uzIZw9oWuvvdb2HKmzzjqLw4cPU1VVxS+//AI4cuR5Izw8nFGjRgHw6aefuqzLyMhg0aJFREZG8u6777J9+3YmTZrEu+++S+/evdm/f7/t8xGEmlLjMSOt9Rat9Vyl1GVKqcGANpb7DgcSBJsEo5vOFKNAoungZNaH559/3us8IU84pywaP3687XY9evRgzZo1fP7551RUVNCuXTu/KYxuueUWAJdxpQ0bNlglLd5//31uv/12OnXqxNNPP81NN93E1q1bGTlypOW5CUJ9UesABq31d8AWoGbZKQXBC8EUo0Ar0zZr1oyQkJCAMx44pyAKJHNEYmIixcXFvPvuuwC88847fttccMEFAC7RdG+//TYhISHccccdXHfddS7bv/POO7Rp04affvpJxpCEeidgMVJKfauUWqqUekQpNQas+UjP+2srCIEQjDEjM1NBoGKklOK6667jo48+CsiLcA52eOaZZ2y3M5O6rl27FoDf/e53ftu0atWK3r17M3/+fIqLi6moqOCDDz7gmmuu8ShmERERltf0wAMPiHck1Cs18YyygFTtSPUzVyl1rVKqe92aJQjB8YwGDBjANddcQ//+/QNu27FjRw4cOBDQPB5nz8hOWLeJWSJ93759REREWKXS/dG/f3++/fZbpk2bxtq1a9m3bx9Dhw71un1aWhqtWrVi27ZtLFmyxOM2U6dOpX///jRv3pxhw4b5zJ0nCN7wK0aG2Aw2/9daf2JmVzA8ok+A5IaaZyScPgRDjDZt2sStt97K9ddfH3DbTz75hKNHjwbUxtkzCiRowqwkW1JSwvHjx20LgCk8P/30kyUuU6dO9bp9aGgo9913H+A5SvDPf/4z119/PWvXruXo0aOsXbuWCy64gHfffVc8KSEg/IqRITaHDFEaY7y6e9hGCtwJdUowouk+/vhjrr322ho9SLdt2xZwG2fPKBBMETPt7Ny5s612pkeVn59vhYSbhf68YY4l/fDDDy7Lv/jiC6tr8d577+XIkSNs2bKFCy+8kDvvvNMSMUGwg61uOiNy7hOt9ada608BZYjTtU7iVFyfhgqnH8EYM9q3bx9VVVV89NFHAbdNT0/nyiuvDKiNKSpRUVFkZ9vNiFV9Yq1zqXVfdOnSBYAtW7bw448/Ao4wcV/069eP7t27M3nyZCuMXGtNRkYGISEhXHHFFbz++us0a9aM5s2b89lnn1n1l5zz5wmCL2o6z8gUp08McYoHxjp5TtJlJ9SaYHTTmclEndP72CUrK4svvvgioDamqHTo0MG2oICj2J7z2JRdz8gUoyNHjlj1l3r27OmzjVKKCy+8kFWrVlldevPmzWPdunUkJyfz/PPPu9jSqlUrPvzwQ8CRiki66wQ71CSa7lql1FSl1F3mMq31chw1hUzPScK8hVoTDDFqaMxuukCyNoBDIEwhCw8Ptx380LZtW8LCwrjssssoLy8HTpZp94U51mTmw/vss8+Ijo5m4cKFHoM9rrnmGvr27cuOHTv45JNPfO778OHDHnPqCacXNfGMWgHZwO+UUgeMMO8DQIG5gSFOglArgiFG5q/4hspsbQqK8+TXQNuOHDnS9hylkJAQunbtSklJCeAIZbeTSy8lJQU4GUqek5PD4MGDiYyM9NrGHE964QXPNTbXrl3L4MGDadmyJTExMVx++eUcOXLE1nkIpx41EaODQIGRBqg1juJ1CcbkV0GoM4IRwDB8uCPnbrdu3RrkeKZnVBMxMr2pQNIPgSMEfdmyZQBcfPHFttoMGDCAkJAQ9u7dy5YtW9iwYQNLlixhy5YtXtuMHDmS9u3b065du2rrfvnlF4YMGcLKlSsJDQ0lLCyMnJwcrrzyShGk05SAxciInIs3I+q01ssDLaQnCHYYOnQo11xzDSEhdVbpxC/9+/fn5ptvJiHBY63GOqcuPKNAsoSDQwCPHz8O2A8nj4yMpGPHjgBWNoZmzZrRo0cPr21CQkK49tprmTt3LqWlpdbyqqoq0tLSKCkpoU+fPmzdupWSkhKysrJYvHgx48aNs8azhNOHmgYwLNdab61jWwTBhT/84Q/VknrWN/379+fDDz+0VeSuLoiKiiIkJCTgMSM4KUKBhoeb2RsAq8KsHRITExkyZIjlqfbq1ctvd+aoUaMoKSnhpptuspbNmjWLVatWcfPNN/P1119bgRvp6em89tprzJkzh0suuUQCH04zfHY0K6UexTFGVBMOGFkaBEHwglKKW265xVY6H3dMAQvUMxo8eLD1vl+/frbbderUiWXLlvHbb78B2MpSYXYDfv3111RVVRESEsKLL75IQkIC7777brWxrnvuuYfXXnuN9evXW4lbhdMDf6Oe9ic+CIJQI957770atatpN50ZjACBiVGHDh3Ytm2bNVnWTnbymJgYOnXqxM6dO1m1ahWdO3dm0aJFXHHFFR69KqUUH3/8MWeffTaTJk3itttu87jd9u3b+fTTT4mLi+Paa68lNjbW9nkIjROf3XRGup8avxrqJAThdMT0jGrTTedvjpEz7du358SJE2zatIl+/fq5FBT0xaBBgwBH9Nw333wDONIuhYaGetx+8ODB9OzZk507d7J48eJq61esWMGgQYN44IEHuO2222jfvj0vvPCCdOs1cRpuZFgQhDqlpp6Rc6i8vxpIzpjbHjhwgAEDBtgWMrN0xapVq5gzZw5KKX7/+9/7bPPQQw8B8MYbb7gsLy8v57rrruPw4cMA9O7dm7KyMh599FHJLN7EETEShCZKTQMYAN566y1eeOGFgOZTOQtXaWmpS4ScL5KTkwE4evQoK1euRGvtMm7lidGjRwMwcOBAl+XZ2dls2rSJiy66iDlz5rBx40b27NnDBRdcwGuvvRZQGQ6hcSFiJAhNlJoGMADcfffdVu0muziL0eeff247ZZKZHWLEiBHWvCSz684bnTt3pnfv3nz//ffWsoqKCl566SUuvPBCvvvuOyvoo3Xr1nz//ffcfPPNTJo0yarBJDQtaiVGSqkeSqlpSqnNSqlfnF5TpcaRINQvZhmJmoSF1wT3lEN250a1bt0agM2bN3Ps2DHCw8NtBU5cfPHFfPvtt3z++ecAfPvtt2zdutVjF19ISAivv/46YWFhPPzww5JeqAlSW89ouJGJoafWupfT6zqkpIQg1CujR4/m7bffthXVVheYmRRiYmKIjo62nRnDzBDx5JNPAo4yHTExMX7bnX322Rw/fpwPPvjAagewY8cOj9u3aNGC8ePHc/z4ca/ddVprXnzxRa666ipmzJhhy36hgdBa1/gFXOZj3fDa7Ls+XsnJyVoQhJqzZMkSfeutt+oOHTrYblNVVaVDQ0M1oAG9YsUKW+1ycnI0oM844wx94sQJ3bx5cw3ohQsXem1z5MgRHRISotu1a+dx/eTJkzWgw8LCNKCHDBmiCwsLbZ/L6QiwTDfA87m2nlGiUuoZpdTdTuUj7lZKPQMEnt9EEIRGza5du/jggw8CSl/knGEc4L///a+tdr169bKOuW7dOo4cOUJISIgVEOGJZs2a0b9/f/bu3cu+fftc1h07dozHHnuM4cOHM3z4cDp37szPP/9M3759pe5SI6BWYqS1fgvHxNh4YIjxigeytaOUhCAIpxArVqwACLj4oLMYmfvwR5cuXQgLC6OystISi169evnt4hszZgyANafJJDMzk8OHD/O3v/2NOXPmUFBQwNChQykrK+Pqq68OuGS8ULfUOppOOwrtPa+1fsx4Pa+19p7Ktw5QSiUopdKVUmlKqYlKqbj6PJ4gCA5MUXGeOGsHM9gCIDo62labkJAQOnbsSHh4OMuWLbOK/PnjlltuAbDmIoFjOOL5558nNjbW2kdkZCSzZ8+mS5cuhIeHW2U1hOBQb6HdSqnB9bVvIE1rna21nqG1fg54vB6PJQiCgTmnac6cOQG1a9OmjfXeTvCCyaBBg+jTpw+7d+9m8ODBvP76637b9OjRg/bt27N06VJr2X/+8x/KysoYPXq0y9yquLg4Zs2axdGjR7nxxhspKyvzul+tNZMnT+bOO+/kwIEDts+hqVJeXh6wB1wbahvaPdjbC7iubkz0yAi3/+Pq8ViCIBiYnlGgc3nMUhPNmjWz7RmBIxLvyJEjrF69mkGDBhEVFeW3jVKK+Ph4/vvf/7Jnzx4AXnzxRcDRVeeOKXJz586le/fuVFZWetzve++9x/3338+7777L2WefzfTp0zlx4oTtc2lqPPLII5aX2RDU1jO6HodXco+HV72GdiulcpRScUqpVGB6fR5LEAQHZuCCrwqvnjDDu1u3bk2nTp1st4uKirKCEX788UfWrVtnq12nTp0oKyvjp59+ori4mI0bN9KlSxerXIU7d955J5dffjl79uzhb3/7m8dtsrKy6Nu3L0uXLqWkpIRx48YxZswYr+LV2Dh69CiTJ092mYO1Z88ezjnnHI/XderUqQ1pXq3F6BkcwQr3uL+Ax+rAPm+MxVHaYguQpLXOrcdjCYJgYGY9GDp0aEDtzHx4N998M0899ZTtdpGRkZSXl1NZWcmGDRtsBxmYGR5WrlzJv//9b7TWpKWl+Wwzbdo0lFK89dZb1XLczZgxg59++okhQ4aQkpJiBUl88cUXvP/++7bPJxg88MADfPTRR0yePJn777+ft956y1o3efJkli1b5rIM4MiRI9WiEeub2kbTHdJaz/WyzuPyOiIFyDBemUqpid42NAIdlimlljX0xRWEUw0zKCDQyrSJiYlAYONF4PCMnD0PuyXWzSSumzdvJjc3l969e/P4476Hllu2bEmfPn3Ys2cP69evZ9++ffz973/nyy+/ZPz48QBMmDABcOTIMx/g//jHPwI6p4bm1Vdf5ZZbbrEE9ssvv7TWmeNqsbGxzJgxw0oR5aucfL0R6MQk4G4f6wYDg2sy4QlIBzJ9vFKN7RKAiU7tEoAiIM7fMWTSqyDUjrVr12pAP/HEEwG1O3r0qDXpdcaMGbbb/fWvf7XaAXr//v222s2ePVsD+txzz9WhoaH6z3/+s612EydO1IB+44039DXXXKMBHRcXZ02+dadfv34a0LNnz7Z9Tg1JSUmJde0GDBigAd2iRQs9adIk/dhjj+no6GgN6Ntvv93arqqqSi9ZssT5uje+Sa9KqZY4usY8rRsPvA1MUEpdFsh+AbQjOi7Dx8vsiksFcp3aFeKY61TTirSCINjkrLPO4qOPPmLiRK+dER5x9oh2795tu5372JRziLgvunXrBkBJSQmVlZW2M4ybXXkHDx5k1apVhIeHU1xc7LLOmT/+8Y8ALFmyxGpXUlJCSUkJzz33HIcO1b6sm3uQRFFREUopW2M6+/fvt96vXr0acHi3b7/9Nm+//bZ1XVauXGltZ55DQxOQGGlHwbwEo9vrF6XUGKfV6TjSA90LJNWlkW7k4iE4whAlQRDqEaUUN910U8CVVZ3DqQOJpnOOnmvevLnXgnzudOnSBcDKLO4+AdYbycnJREREUFRUxK+//sqtt94KOMTNUwn08847D3BE5JWVldG6dWvS0tL45ptvyMjIYMiQIdXGnwJhxYoVxMTEuAjP119/DcArr7xiLSssLOTXX39ly5YthIaGsnz5cqqqqti7d6/H/e7Zs8cSqrCwMJYvX26t27lzpyVGgZQYqS01GTNKxDFWczkwRCllpgxO0Fqbs8zqrcqrITqFxmTXNKVUOpBVX8cTBKFuMIvs1VSM3Avt+aJly5ZERkZSWOj4jdq5c2db7UJCQujQoQPffvstJ06c4JxzzuHAgQOsWbPGY9kLc2xq48aNVmDG119/zdatWwFHRdsff/yRm266iZkzZ7Jx40a++eYb3njjDR599FG/Htt7771HRUUF119/vTVeN3v2bMDVS0xMTKRbt25MnTqVqqoqrrvuOmJiYqwS86aoeBqzc0+vtGPHDkuMTA+zIQirQZuf9cnghMeMLrnv3Lap13KLWmtJtysITYzXXnuN5OTkgIIYnLvp7MwxMlFK0axZM8rLy4Hq5S98UVFRwapVqwDH/ChfQROxsbF07tyZmTNnkp+fb9lZUFBgbXPjjTeyc+dOK+v4ZZddxnffOR6Z27dvZ8qUKR73XVVV5ZJZfPPmzSQlJTFv3jwAl2OYbN++HYBffvnFZXlSUhJ5eXkMGDCAn376yWVd//79+emnn0hJSWHZsmXs2LHD8kC7d+9uCWt9UxPPaIhSqptSqoVS6hEcXkpLwNmfi6sT6wRBOGUICwtj0KBBLtkY/OEsQBkZGQEdLz4+3nofiBg5e1HZ2dlMmzbN5/a9e/e2PJzQ0FDKyspYuHChFUG4bds27rnnHmt7U4jAUaPJG8uXL2fHjh1WW7PLzRyHKiws5Pjx4y5tnCPlnDG9H091pHr37g1Aeno6EDzPqCZilAXMxRHI0AZHNoRMIEMp9YgxjlRv3XSCIDRN7rvvPnr27Gkrv5yJsxh58gR80bFjR+t9IGJkjjcBfPLJJ2zatMnn9r1792b9+vXAyYf36tWrSUxMRClF9+7dKSgooHv37lbmisGDBxMbG8vBgwe9zp3KyckBHHOzwDHOo7WmrKzMyhSxZcsWy/sDvHoxZi7B6OhoF880Ojqahx56iM8//5y7776bNm3auIwZde/e3ee51yUBi5F2JEbtqbVurR2JUd/Sjomub+EQqdbGe0EQBIvKykqX5KV2CDTTgzNmCqIOHTpw33332W53ySWXACcfxM7i5AnnpLHOY0DHjh1Da01kZCQ//PADl19+Offccw+zZs1iyZIl3HbbbQC8/fbbHvebk5PDgAEDGDhwIODwjEzhSUpyxIg9/PDDXr2XIUOGWO9NASorK3MpH9++fXvCw8O56qqrUErRuXNnF89o7NixPs+9LqnxpFejm26wUwADWuvlIkSCIHjixx9/JCcnh23bttlu4+wZ2UmS6oz50D148KDtqrQAqampLn+7du3qc/ubbrrJer97925LnFJTU7nqqqvYsGEDxcXFJCUlERoayqhRo4iKiuKGG24APHsfr732Gt999x1XXnklzZo1Iyoqii1btlhid/755xMbG8uXX35p5d976aWXGDBggLWPs846y3pveqOjRo1y8RKdhQlgwIABLFq0iH379hEZGWmrPHxdUSMxUkq9CRTjCFwoUko1bBIjQRCaLIGEOjuL0bBhwwI6jlkmPTo6OqBaRd26dWPmzJmWZ+VPjNq1a8fUqVMZN24cs2bNYsGCBTz66KPceeedLh7Z2Wef7dLOHFPauXOny/JNmzbx8MMPc+WVV3LHHXeglKJdu3a8+eab1rhOixYtuPjii13aJScns2rVKq655hrAVWgGDBhAaWkpo0ePpn379kRGRhIREVFNjCZMmEBxcTHZ2dkBZ8uoLQGLkVLqUSBHax2itW6ltQ4FphnBDIIgCD5piGg6OClG999/f0Dpi2JiYrj66qtp3bo1bdq08Zpc1Zlx48YxdepURo4cSfv27bn99tv58ssvXcLBnb0WwBKFzz77jMmTJ1vLn3vuOSIjIznrrLMYOXIk5eXltG3blq5duzJz5kzAIejDhw932V/btm2Bk+U62rVrx44dOyyxM6/fkCFDOO+88xg4cGA1my688EK6du3K8ePHG78YAYVa60+cFxj/S9CCIAhe6dOnD1DzeUY1FaOEhISA2gHMnz+fgQMHWt1VgbRbsGABOTk53HfffYSGhlreh/t5h4SE0K1bN1asWOEyblRQUECvXr146aWXOO+884iIiKBTp07Ex8dbXZXvvPMOd999t0spD1OEnMWoU6dOLoEcAI8//jjz589n8eLF1ZLWmuNGEHgewdpSEzHy5mOf+tWmBEGoMTfeeCPQcGJkRrPVRIz+9Kc/eax95I8//vGPvPzyy2zevJkWLVrQpk0bNm7c6DUD9hlnnEFUVBTLly+3MiLs3buXvXv3Ehsby0svvYRSivbt27Nnzx4rRPvnn39m/vz53H///da+zFB2ZzHyRXh4uMeMFmb7piBGic5BCwBKqe7AEM+bC4IgOGoMDRs2jLAw+3Pta9NN16tXL7Zu3WpFxwWC1prPP//cCq+2S3x8PIcOHeKXX36hV69eKKVo2bKl17lVnTt3tuYKmZNZd+3axa5du6xQa3CI1p49e6xuuscff5wrrrjCZV/mdbUrRt4wu/uaghhlA98ppZYqpb5RSi0FcoB/1q1pgiCcStx1113Mnz8/oDbOAlSTMO+uXbvWKL/amjVrgMDmJ4EjRU9xcTGbN2+2UgX5onPnzuzfv5/Q0FBWrlzJiRMnKCoqoqqqyiUx6/3338+ZZ57JP//peMyOHDnSuh4PP/ywy9yt0aNH8/zzz7uMVwVSkbbJiJF21DBKwTHRNRd4VmvdyykvnSAIQp1gilFYWFhAHlVtefvttxk2bBj9+/cPqF1cXBz79u1j69attsWosrKSvn37Ul5ebnXnPfroo1YSVnB4O87RhNHR0bzxxhs89thjvPDCCyxatMha17JlS5KTk/n8888BR22iM844gw8++MDWOZhi5FwRtiGo8TwjrfUMrfXzZjCDUuruujNLEAThpDcUaBddbTG9uEC9qri4OI4dO8aBAwf43//9X7/bm8ECq1ev5tJLL7XmDJ133nmEhLg+np0jAqOioli3bh2vvPJKtYnEK1euZOTIkVaEXkVFBX369OG2227jtdde82uTKUZmxvOGwq8YGV1x/l7f4vCUBEEQ6oxgiVFNeeihh/juu+9o2bKlraq0znnwsrKy2LFjBwAHDlSPB3MWo+joaMaMGUN5eTkLFy60lmutuf/++4mNjeWjjz4CHGNnc+fOZfTo0Tz44IMsXrzYp02mGDV0TSM7fq/CntBMqKUtgiAILiiliIyMbDJi1KNHD/Ly8nj88cf53//9X+vB7g1nMVq6dKkV+eepKF+LFifjxqKjo62w7wULFjBy5EjAEQTxww8/kJWV5TKhNSwsjA8//JB+/foxYcIEVqxY4bU2lBkA0dCekR0xytBaL/e3kVLqYB3YIwiC4EJTEqMtW7Zw5513smrVKv70pz/53d5ZMHbt2mUVzDv33HOrbevuGUVHRzNkyBAWLFhgLf/oo49o3rw5t9xyS7X2zZs359VXX2XBggWUlpbSrFkzjzY12m46O0IUyHaCIAiBEBUV1WTEaMWKFVYtJH9eETg8luzs7GrBGe6pg6C6ZwRYmSKqqqoAx1hTRkaG17lcV199NS+//LJXIQKs7sWIiAiXchf1japNSdymRkpKil62bFmwzRAEIQC6detGu3btWLp0abBN8cu8efO47LLLgMBy8F100UUsWrQIpRQxMTEec+ktWLDAmjNVVVVV45LgFRUVLFiwgN69e3vMSH748GE6derk7BnlGRHU9UqNo+kEQRAagqbUTedcCjwQzHEaU4w84ewZeRKiXbt2cfCg/9GSPXv2MGLECN577z2vx/nxxx/ZsWMH77//vg3r6wYRI0EQGjVNqZuupmJkjh2dffbZVni3O56SvR46dIjevXvz5ptv8o9//IPu3btbXXbe6Ny5MxdccAGzZs1yWf6f//yHBx98kOPHj9O/f386derErbfeWqPzqQkiRoIgNGr69+9P3759g22GLUwxCjSvnVkgr1+/fl6735w9I+dlO3fuZMOGDSxfvpxBgwZVm5/kiWHDhrFixQorfFtrzQsvvMDixYsDqv1Ul4gYCYLQqPn444/517/+FWwzbNGyZUvWr1/PPffcE1A7M6Dgl19+4emnn/a6b3eUUiQmJlJQUMDGjRttF8M7//zzqaiowBxDX7x4MRs2bOCee+6p8VhUbRExEgRBqCNCQkLo06ePRy/GF2VlZQAsWbKEr776yuM23jyWxMRECgsLOXTokK2JtoCVamjJkiUAzJ49m/DwcMaNGxeQ3XVJwyV7EgRBEDxiljiHwEpsgCM799y5c6moqLAtgm3atGHlypVWafIff/yRs88+22fId30jYiQIghBkzj77bFq1asXBgwcDFqPzzz+fkpIShg0bZtU7ssPAgQNd3tupaFufiBgJgiA0AkwR8le6wX1M57bbbuO2224L+Hi5ubn8+OOP/OUvf+HVV18NuH1dI2NGgiAIjQAzfN2s2OqJvXv3WhVhnTly5Ag///xztQzevliwYAFPPvlkQLWO6pNGLUZKqSSlVJ6H5QlKqYlKqVTjb1wQzBMEQagzLr74Yv7+97+TlZXldZu2bdtWC1KYPn06LVq04Nxzzw0oS0XXrl2pqqpi0qRJxMbGeswU3pA02m46pVQqcBBI8rA6S2s9wtiuEEdWcckaLghCk+Wdd96pUTvnCriBRPF17doVgLy8PEpLSz2GjjckjdYz0lrnaq3z3ZcrpRLctisEghePKAiCUEekpaUxbdq0gNo4i0jz5s1ttzPFKD8/n/j4+AatpOuJRitGPkgCit0XuouUIAhCU+LWW2/lk08+4ZdffgmonbMY1cQzKi4utpVhvL5pimLkaVbXQSDO08ZKqXSl1DKl1DKzvrwgCEJjY/v27UDg84ycBSgQzyg2NpZDhw5x6aWXWolag0mjHTOqK7TW2UA2OEpIBNkcQRAEnwQqRm3atOGSSy6hX79+xMbGBtS2RYsWjB492mvV14akQcVIKZUOJPrYJEdrnetnN568oFZ46LoTBEFoKpj1j2riGc2bN69Gx5w2bRqbNm3i9ddfr1H7uqRBxcjwUmpLPh666oxABkEQhCZJREQE4HuekTeWL1/O4cOHGTZsWEDtVq9ezeTJk3n55Zet4weLJtdNp7UudJ6BbAQuBBZ+IgiC0MgYN24c7dq1Y/To0QG3TUpyzIAJtHK3OVb0wAMP8MYbbwR83Lqk0YqRMc8oyXifiWsX3lil1EQcXlKS1lrmGAmC0KS56667uOuuuxr0mK1btwYCF7H6oNFG0xnzjJ7TWiutdYbzWJLWutBYl6u1fi6YdgqCINQFS5cu5fLLL2f9+vUNdsyxY8fy0EMP8fe//73BjumNRusZCYIgnE48+eST5OTksHfvXqu0g102b95co4i4yMhIXnrppYDb1QciRoIgCI2Ao0ePAtWzctshMdFXkHLToNF20wmCIJxOjBo1CjiZGeF0Q8RIEAShEfDII4+wd+9eunfvHmxTgoKIkSAIQiNAKdUocsQFCxEjQRAEIeiIGAmCIAhBR8RIEARBCDoiRoIgCELQETESBEEQgo6IkSAIghB0RIwEQRCEoCNiJAiCIAQdESNBEAQh6IgYCYIgCEFHxEgQBEEIOiJGgiAIQtARMRIEQRCCjoiRIAiCEHREjARBEISgI2IkCIIgBB0RI0EQBCHoiBgJgiAIQUfESBAEQQg6IkaCIAhC0AkLtgG+UEolAW9prZM9LE81/j0HGK+1Lm5g8wRBEIQ6otGKkVIqFTgIJLktjwNStNbPGf+nAXOBZPd9CIIgCE2DRttNp7XO1Vrne1iVAmQ4/Z8LJBkiJQiCIDRBGq0YeUNrnQuMdVqUYCwvDopBgiAIQq1pcmIE4OYxXQc8521bpVS6UmqZUmrZvn376t84QRAEIWCapBiZGF1zSVrrDG/baK2ztdYpWuuUtm3bNpxxgiAIgm0aNIBBKZUOJPrYJMfohrNLJq5ddoIgCEITpEHFSGudXVf7UkpNBDK01sVKqTgZMxIEQWi6NMluOiOce4aTAKX62FwQBEFo5DT2eUZJxvtMjC48pVQCMN1Ybm5eCMwIhp2CIAhC7Wm0YmSMHeXiFimntS4ElMdGgiAIQpOkSXbTCYIgCKcWIkaCIAhC0BExEgRBEIKOiJEgCIIQdESMBEEQhKAjYiQIgiAEHREjQRAEIeiIGAmCIAhBR8RIEARBCDoiRoIgCELQETESBEEQgo6IkSAIghB0RIwEQRCEoCNiJAiCIAQdESNBEAQh6IgYCYIgCEFHxEgQBEEIOiJGgiAIQtARMRIEQRCCjoiRIAiCEHSU1jrYNjQYSqkjwMZg22GTNsD+YBthE7G1fmgqtjYVO0FsrQndtNZt6/sgYfV9gEbGRq11SrCNsINSapnYWveIrXVPU7ETxNbGjHTTCYIgCEFHxEgQBEEIOqebGGUH24AAEFvrB7G17mkqdoLY2mg5rQIYBEEQhMbJ6eYZCYIgCI0QESNBEAQh6JzyYqSUSlBKTVRKpRp/44JtkzNKqSTDrolKqenO9imlMpVS6UqpOONvWhBN9WlPY7vO7tfSbV1Qr6vxmed5WO71Ggbr+vqwtdHdtz5sbXT3rQ9bG+19W+9orU/pF5Dj9D4ByAq2TU72xAHpTv+nAXlO/2cCGigCJjYCe73a09ius2Gn+2tisK8rkAokOb561dZ5vYbBuL7ebG2M962f69qo7ls/tjbK+7ZBPsNgG1DPH3qC881mLCsKtl1OtqQCBU7/xxk3W5zxf1qwbXSz16M9je06G/akui1zeXg2gmup7V7DYF9fD7Y22vvWywO+Ud63Xu6BRn3f1ufrVO+mSwKK3RcqpRIa3pTqaK1zgbFOixKM5cXO2ymlkhrQLL94sKexXeeDxrU17UgDprlv1Miuq69r2Kiur9y39UZTvG/rjFNdjFp5WHYQxy+5RoHWOt/p3+uA55z+T1BKpQKFRn9xsG9Cb/Y0quvs/FA0+t9buT0oG9t1Bd/XsFFdX5D7tj5oovdtnXG65aZrtBg3X5LWeoS5TGv9nNP6LCAHSGx46xqnPTZ5HHjGeUETPY9Gidy39cZpd9+e6p6Rp185rfDgmjcCMnHt+jC/6ABorQsxukOChQ97GvN1TvXQfRRnvm8M19XA1zVszNdX7tv6oanct3XGqS5G+XhwxY0PstGglJoIZGiti80bznDH5wbVMCf82NMor7Nh80EPyxrNdXXC1zVsrNdX7tt6oIndt3XGKS1G7jeVMTBZbUAwmBiDlDOcfgWlGn+X4eSmm9s1rHUueLWnEV9nTwPUje26Ar6vYWO8vnLf1itN5r6tS06HMaOxxi+4fBx92xOCbZCJcfNPN96biwsxvuRKqULD9mIgUWs91uOOGgAb9jTG61wMLHVeEOzravzCTTLeZ+IILTYjqHxdwwa/vt5sbYz3rTdbG+N96+cegEZ43zYEkihVEARBCDqndDedIAiC0DQQMRIEQRCCjoiRIAiCEHREjARBEISgI2IkCIIgBB0RI+GUxwiHPWUJ5PyMUGJBaHSIGAmnNEqpdII0OdAo2DZdKaWNxJYTjb/T67gw2gw7guQpG7VR5M20MdVTO2O7HKVUwaku7ELwOB0mvQqnKcbDN1Frne22PA6Yq7VOrs/jGxNEk3DU+clws6FIKYXWutZCqbUuVEq1Vkol+EllkwZkubXNV0pNNf71mOvMEKlWOFL/nFKz/oXGg3hGwqnMBNwevgbVcn/VIyNwZFd25yBwTh0e5xkgw88253gRq1Y4MihUywBtCPdBHBkDct3XC0JdIWIknMokeXn4ehOI+iAVt4e4k6fhSShrhJEjzmsWZ8NLXOptPVCAkaLGDbPrrtg9i7Qg1CUiRsIpifHwdU+EmWqMeaQDret7/MPoois2C9EppRKMMawJQHI9ZIcu9FFwbQKQ7b7QuE7LcFyrBLd1ppBWE1RBqGtkzEg4VUnALfOxMYazDHjcfQynnjCrcprBCufgGD+qVYJLpVScFy+lAMd553tY561NktZ6hiFKlhiZJSGMBJ0jMBKjCkJ9IZ6RcKqSgOPh7E5D/sofAWRprWcYrwygVR14ZN6EoRgPXXWGt5TnpU0rOFlOwamAW6pTJmnxjIR6Rzwj4VQmzsMyv+NFhidznZ99H7RRbiAVR/eYu001LhVtdJ35sr/Yw7IJNksjFAIJRlkIs2vR7GpsVAUphVMPESPhVOUgnh/6qThKZeMtFNoIX65VCLMhGp4e4im4iYnxwE/BGOMyuhNTcXg5Znszmi0DyPFiexxu42R+bDTHi0zMcSNrnAvxioQGQrrphFOVQjx7Rq2MeTlxeI4eqytG4PqgN4nDEAylVJphR6YxF+qgsTwOwFiWaeznHFMgtNbPefFUEvEQtIH3br00J9HBaDvBrdBbQ0YeCqcx4hkJpyTGZE5Poc7ZZqaB+pjAaQhJOoboKKXS3SbdZgDnOHWFjQPyDe+o0EkcTO9oqhF4YCfoIcFNXADGunfRmQIIjFNKJTqtz8MINzfsSeVkEEaSh30LQp0hlV6FUxalVBaOrAHFwbbFG2aknSmMhoCak3UzOTmRNc74e53WOsNdHJw8LHfhyWokJeAFwSfSTSecymQCjwfbCF8YIpRg5IgzM0OYIdrP4PBMTOEpBA4opVI9eCnpGGNhJobQ1dnEWkGoT8QzEk5pjDDqGadyNJjhTaV6yMEnXpHQZBDPSDil0Vo/hyNB6KlMmrsQGXiaZyUIjRLxjARBEISgI56RIAiCEHREjARBEISgI2IkCIIgBB0RI0EQBCHoiBgJgiAIQUfESBAEQQg6/x8Ule3eI1L/KAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.pyplot import savefig\n",
    "\n",
    "from matplotlib import rc\n",
    "rc('text', usetex=True)\n",
    "\n",
    "\n",
    "if EvolOption == \"low resolution\":\n",
    "    extraction_radius = \"30.20\"\n",
    "    Amplitude = 1.8e-2\n",
    "    Phase     = 2.8\n",
    "elif EvolOption == \"high resolution\":\n",
    "    extraction_radius = \"33.64\"\n",
    "    Amplitude = 1.8e-2\n",
    "    Phase     = 2.8\n",
    "else:\n",
    "    print(\"Error: output is not tuned for Nr = \"+str(Nr)+\" . Plotting disabled.\")\n",
    "    exit(1)\n",
    "\n",
    "# Transposed for easier unpacking:\n",
    "t,psi4r,psi4i = np.loadtxt(os.path.join(outdir, \"outpsi4_l2_m0-r\"+extraction_radius+\".txt\")).T\n",
    "\n",
    "t_retarded    = []\n",
    "log10abspsi4r = []\n",
    "bh_pert_thry  = []\n",
    "for i in range(len(psi4r)):\n",
    "    retarded_time = t[i]-float(extraction_radius)\n",
    "    t_retarded.append(retarded_time)\n",
    "    log10abspsi4r.append(np.log(float(extraction_radius)*np.abs(psi4r[i]))/np.log(10))\n",
    "    bh_pert_thry.append(np.log(Amplitude*np.exp(-0.0890*retarded_time)*np.abs(np.cos(0.3737*retarded_time+Phase)))/np.log(10))\n",
    "\n",
    "# print(bh_pert_thry)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "plt.title(\"Grav. Wave Agreement with BH perturbation theory\", fontsize=18)\n",
    "plt.xlabel(\"$(t - R_{ext})/M$\", fontsize=16)\n",
    "plt.ylabel('$\\log_{10}|\\psi_4|$', fontsize=16)\n",
    "\n",
    "ax.plot(t_retarded, log10abspsi4r, 'k-', label='SENR/NRPy+ simulation')\n",
    "ax.plot(t_retarded, bh_pert_thry, 'k--', label='BH perturbation theory')\n",
    "#ax.set_xlim([0,t_retarded[len(psi4r1)-1]])\n",
    "ax.set_xlim([0, 1.5*domain_size - float(extraction_radius)])\n",
    "ax.set_ylim([-13.5, -1.5])\n",
    "\n",
    "plt.xticks(size = 14)\n",
    "plt.yticks(size = 14)\n",
    "\n",
    "legend = ax.legend(loc='upper right', shadow=True, fontsize='x-large')\n",
    "legend.get_frame().set_facecolor('C1')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "# Note that you'll need `dvipng` installed to generate the following file:\n",
    "savefig(os.path.join(outdir, \"BHperttheorycompare.png\"), dpi=150)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 10: 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-Two_BHs_Collide-Psi4_new_way.pdf](Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide-Psi4_new_way.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": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-23T04:29:13.123239Z",
     "iopub.status.busy": "2021-09-23T04:29:13.122959Z",
     "iopub.status.idle": "2021-09-23T04:29:16.138662Z",
     "shell.execute_reply": "2021-09-23T04:29:16.138340Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide-\n",
      "    Psi4_new_way.tex, and compiled LaTeX file to PDF file Tutorial-\n",
      "    Start_to_Finish-BSSNCurvilinear-Two_BHs_Collide-Psi4_new_way.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-Two_BHs_Collide-Psi4_new_way\")"
   ]
  }
 ],
 "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}