{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Start-to-Finish Example: Numerical Solution of the Scalar Wave Equation, in Curvilinear Coordinates\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module solves the scalar wave equation in *spherical coordinates* (though other coordinates, including Cartesian, may be chosen), using the Runge-Kutta fourth-order scheme. It includes code generation for initial data and scalar wave RHSs, boundary condition drivers, and numerical error plotting. The notebook verifies the solution's convergence towards the exact solution, validating the curvilinear approach.\n",
    "\n",
    "**Notebook Status:** <font color =\"green\"><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This module has been validated to converge at the expected order to the exact solution (see [plot](#convergence) at bottom).\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [ScalarWave/ScalarWaveCurvilinear_RHSs.py](../edit/ScalarWave/ScalarWaveCurvilinear_RHSs.py) [\\[**tutorial**\\]](Tutorial-ScalarWaveCurvilinear.ipynb) Generates the right-hand side for the Scalar Wave Equation in curvilinear coordinates\n",
    "* [ScalarWave/InitialData.py](../edit/ScalarWave/InitialData.py) [\\[**tutorial**\\]](Tutorial-ScalarWave.ipynb) Generating C code for either plane wave or spherical Gaussian initial data for the scalar wave equation \n",
    "\n",
    "## Introduction:\n",
    "As outlined in the [previous NRPy+ tutorial notebook](Tutorial-ScalarWaveCurvilinear.ipynb), we first use NRPy+ to generate initial data for the scalar wave equation, and then we use it to generate the RHS expressions for [Method of Lines](https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html) time integration based on the [explicit Runge-Kutta fourth-order scheme](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (RK4).\n",
    "\n",
    "The entire algorithm is outlined below, with NRPy+-based components highlighted in <font color='green'>green</font>.\n",
    "\n",
    "1. Allocate memory for gridfunctions, including temporary storage for the RK4 time integration.\n",
    "1. <font color='green'>Set gridfunction values to initial data.</font>\n",
    "1. Evolve the system forward in time using RK4 time integration. At each RK4 substep, do the following:\n",
    "    1. <font color='green'>Evaluate scalar wave RHS expressions.</font>\n",
    "    1. Apply boundary conditions.\n",
    "1. At the end of each iteration in time, output the relative error between numerical and exact solutions."
   ]
  },
  {
   "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](#writec): Generate C code to solve the scalar wave equation in curvilinear coordinates\n",
    "    1. [Step 1.a](#id_rhss): C code generation: Initial data and scalar wave right-hand-sides\n",
    "    1. [Step 1.b](#OuterBoundaryConditions): C code generation: Boundary condition driver\n",
    "    1. [Step 1.c](#cparams_rfm_and_domainsize): Generate Cparameters files; set reference metric parameters, including `domain_size`\n",
    "    1. [Step 1.d](#cfl): C code generation: Finding the minimum proper distance between grid points, needed for [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673)-limited timestep\n",
    "1. [Step 2](#mainc): The C code `main()` function for `ScalarWaveCurvilinear_Playground`\n",
    "1. [Step 3](#compileexec): Compile generated C codes & solve the scalar wave equation\n",
    "1. [Step 4](#convergence): Code validation: Plot the numerical error, and confirm that it converges to zero at expected rate with increasing numerical resolution (sampling)\n",
    "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='writec'></a>\n",
    "\n",
    "# Step 1: Using NRPy+ to generate necessary C code to solve the scalar wave equation in curvilinear, singular coordinates \\[Back to [top](#toc)\\]\n",
    "$$\\label{writec}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='id_rhss'></a>\n",
    "\n",
    "## Step 1.a: C code generation: Initial data and scalar wave RHSs \\[Back to [top](#toc)\\]\n",
    "$$\\label{id_rhss}$$\n",
    "\n",
    "\n",
    "We choose simple plane wave initial data, which is documented in the [Cartesian scalar wave module](Tutorial-ScalarWave.ipynb). Specifically, we implement a monochromatic (single-wavelength) wave traveling in the $\\hat{k}$ direction with speed $c$\n",
    "$$u(\\vec{x},t) = f(\\hat{k}\\cdot\\vec{x} - c t),$$\n",
    "where $\\hat{k}$ is a unit vector.\n",
    "\n",
    "The scalar wave RHSs in curvilinear coordinates (documented [in the previous module](Tutorial-ScalarWaveCurvilinear.ipynb)) are simply the right-hand sides of the scalar wave equation written in curvilinear coordinates\n",
    "\\begin{align}\n",
    "\\partial_t u &= v \\\\\n",
    "\\partial_t v &= c^2 \\left(\\hat{g}^{ij} \\partial_{i} \\partial_{j} u - \\hat{\\Gamma}^i \\partial_i u\\right),\n",
    "\\end{align}\n",
    "where $\\hat{g}^{ij}$ is the inverse reference 3-metric (i.e., the metric corresponding to the underlying coordinate system we choose$-$spherical coordinates in our example below), and $\\hat{\\Gamma}^i$ is the contracted Christoffel symbol $\\hat{\\Gamma}^\\tau = \\hat{g}^{\\mu\\nu} \\hat{\\Gamma}^\\tau_{\\mu\\nu}$.\n",
    "\n",
    "Below we generate \n",
    "+ the initial data by calling `InitialData(Type=\"PlaneWave\")` inside the NRPy+ [ScalarWave/InitialData.py](../edit/ScalarWave/InitialData.py) module (documented in [this NRPy+ Jupyter notebook](Tutorial-ScalarWave.ipynb)), and \n",
    "+ the RHS expressions by calling `ScalarWaveCurvilinear_RHSs()` inside the NRPy+ [ScalarWave/ScalarWaveCurvilinear_RHSs.py](../edit/ScalarWave/ScalarWaveCurvilinear_RHSs.py) module (documented in [this NRPy+ Jupyter notebook](Tutorial-ScalarWaveCurvilinear.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:15.367685Z",
     "iopub.status.busy": "2021-09-28T19:19:15.357463Z",
     "iopub.status.idle": "2021-09-28T19:19:16.359409Z",
     "shell.execute_reply": "2021-09-28T19:19:16.359861Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step P1: Import needed NRPy+ core modules:\n",
    "from outputC import lhrh, add_to_Cfunction_dict  # NRPy+: Core C code output module\n",
    "import finite_difference as fin       # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par        # NRPy+: Parameter interface\n",
    "import grid as gri                    # NRPy+: Functions having to do with numerical grids\n",
    "import reference_metric as rfm        # NRPy+: Reference metric support\n",
    "import cmdline_helper as cmd          # NRPy+: Multi-platform Python command-line interface\n",
    "import shutil, os, sys                # Standard Python modules for multiplatform OS-level functions\n",
    "\n",
    "# Step P2: Create C code output directory:\n",
    "Ccodesrootdir = os.path.join(\"ScalarWaveCurvilinear_Playground_Ccodes\")\n",
    "# First remove C code output directory if it exists\n",
    "# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
    "shutil.rmtree(Ccodesrootdir, ignore_errors=True)\n",
    "# Then create a fresh directory\n",
    "cmd.mkdir(Ccodesrootdir)\n",
    "\n",
    "# Step P3: Create executable output directory:\n",
    "outdir = os.path.join(Ccodesrootdir, \"output\")\n",
    "cmd.mkdir(outdir)\n",
    "\n",
    "# Step P4: Enable/disable SIMD. If enabled, code should run ~2x faster on most CPUs.\n",
    "enable_SIMD = True\n",
    "\n",
    "# Step P5: Enable reference metric precomputation.\n",
    "enable_rfm_precompute = True\n",
    "if enable_rfm_precompute:\n",
    "    par.set_parval_from_str(\"reference_metric::rfm_precompute_to_Cfunctions_and_NRPy_basic_defines\", \"True\")\n",
    "else:\n",
    "    par.set_parval_from_str(\"reference_metric::rfm_precompute_to_Cfunctions_and_NRPy_basic_defines\", \"False\")\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 P6: 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 1: Set some core parameters, including CoordSystem, boundary condition,\n",
    "#                                             MoL, timestepping algorithm, FD order,\n",
    "#                                             floating point precision, and CFL factor:\n",
    "\n",
    "# Step 1.a: Set the coordinate system for the numerical grid\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 1.c: 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",
    "# 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     = 10.0 # Needed for all coordinate systems.\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.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    = 1.0 # If SymTP chosen\n",
    "\n",
    "# Step 1.d: Set the order of spatial and temporal derivatives;\n",
    "#           the core data type, and the CFL factor.\n",
    "# RK_method choices include: Euler, \"RK2 Heun\", \"RK2 MP\", \"RK2 Ralston\", RK3, \"RK3 Heun\", \"RK3 Ralston\",\n",
    "#              SSPRK3, RK4, DP5, DP5alt, CK5, DP6, L6, DP8\n",
    "RK_method = \"RK4\"\n",
    "FD_order  = 4        # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
    "REAL      = \"double\" # Best to use double here.\n",
    "CFL_FACTOR= 1.0\n",
    "outer_bc_type = \"RADIATION_OUTER_BCS\"  # can be EXTRAPOLATION_OUTER_BCS or RADIATION_OUTER_BCS\n",
    "radiation_BC_FD_order = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:16.364858Z",
     "iopub.status.busy": "2021-09-28T19:19:16.364326Z",
     "iopub.status.idle": "2021-09-28T19:19:17.352957Z",
     "shell.execute_reply": "2021-09-28T19:19:17.352480Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2: Import the ScalarWave.InitialData module.\n",
    "#         This command only declares ScalarWave initial data\n",
    "#         parameters and the InitialData() function.\n",
    "import ScalarWave.InitialData as swid\n",
    "\n",
    "# Step 3: Import ScalarWave_RHSs module.\n",
    "#         This command only declares ScalarWave RHS parameters\n",
    "#         and the ScalarWave_RHSs function (called later)\n",
    "import ScalarWave.ScalarWaveCurvilinear_RHSs as swrhs\n",
    "\n",
    "# Step 4: Call the InitialData() function to set up initial data.\n",
    "#         Options include:\n",
    "#    \"PlaneWave\": monochromatic (single frequency/wavelength) plane wave\n",
    "#    \"SphericalGaussian\": spherically symmetric Gaussian, with default stdev=3\n",
    "swid.InitialData(CoordSystem=CoordSystem, WaveType=\"SphericalGaussian\")\n",
    "\n",
    "# Step 5: Set the finite differencing order to FD_order (set above).\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", FD_order)\n",
    "\n",
    "# Step 6: Generate SymPy symbolic expressions for\n",
    "#         uu_rhs and vv_rhs; the ScalarWave RHSs.\n",
    "#         This function also declares the uu and vv\n",
    "#         gridfunctions, which need to be declared\n",
    "#         to output even the initial data to C file.\n",
    "# First get into the enable_rfm_precompute environment, if enable_rfm_precompute==True\n",
    "if enable_rfm_precompute:\n",
    "    cmd.mkdir(os.path.join(Ccodesrootdir, \"rfm_files/\"))\n",
    "    par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"True\")\n",
    "    par.set_parval_from_str(\"reference_metric::rfm_precompute_Ccode_outdir\", os.path.join(Ccodesrootdir, \"rfm_files/\"))\n",
    "\n",
    "swrhs.ScalarWaveCurvilinear_RHSs()\n",
    "# Step 6.a: Now that we are finished with all the rfm hatted\n",
    "#           quantities, let's restore them to their closed-\n",
    "#           form expressions.\n",
    "if enable_rfm_precompute:\n",
    "    par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\", \"False\") # Reset to False to disable rfm_precompute.\n",
    "rfm.ref_metric__hatted_quantities()\n",
    "\n",
    "# Step 7: 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 8: Set enable \"FD functions\" parameter. See above for details.\n",
    "par.set_parval_from_str(\"finite_difference::enable_FD_functions\", enable_FD_functions)\n",
    "\n",
    "# Step 9: If enable_SIMD, then copy SIMD/SIMD_intrinsics.h to $Ccodesrootdir/SIMD/SIMD_intrinsics.h\n",
    "if enable_SIMD:\n",
    "    shutil.copy(os.path.join(\"SIMD\", \"SIMD_intrinsics.h\"), os.path.join(Ccodesrootdir, \"SIMD\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mol'></a>\n",
    "\n",
    "## Step 1.b: Generate Method of Lines timestepping code \\[Back to [top](#toc)\\]\n",
    "$$\\label{mol}$$\n",
    "\n",
    "The Method of Lines algorithm is described in detail in the [**NRPy+ tutorial notebook on Method of Lines algorithm**](Tutorial-Method_of_Lines-C_Code_Generation.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:17.357055Z",
     "iopub.status.busy": "2021-09-28T19:19:17.356544Z",
     "iopub.status.idle": "2021-09-28T19:19:17.365910Z",
     "shell.execute_reply": "2021-09-28T19:19:17.365519Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 10: Generate Runge-Kutta-based (RK-based) timestepping code.\n",
    "#       As described above the Table of Contents, this is a 2-step process:\n",
    "#       1.b.A: Evaluate RHSs (RHS_string)\n",
    "#       1.b.B: Apply boundary conditions (post_RHS_string, pt 1)\n",
    "import MoLtimestepping.MoL as MoL\n",
    "# from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict\n",
    "# RK_order  = Butcher_dict[RK_method][1]\n",
    "RHS_string = \"rhs_eval(params, rfmstruct, RK_INPUT_GFS, RK_OUTPUT_GFS);\"\n",
    "if not enable_rfm_precompute:\n",
    "    RHS_string = RHS_string.replace(\"rfmstruct\", \"xx\")\n",
    "RHS_string += \"\"\"\n",
    "if(params->outer_bc_type == RADIATION_OUTER_BCS)\n",
    "  apply_bcs_outerradiation_and_inner(params, bcstruct, griddata->xx,\n",
    "                                     gridfunctions_wavespeed,gridfunctions_f_infinity,\n",
    "                                     RK_INPUT_GFS, RK_OUTPUT_GFS);\"\"\"\n",
    "\n",
    "# Extrapolation BCs are applied to the evolved gridfunctions themselves after the MoL update\n",
    "post_RHS_string = \"\"\"if(params->outer_bc_type == EXTRAPOLATION_OUTER_BCS)\n",
    "  apply_bcs_outerextrap_and_inner(params, bcstruct, RK_OUTPUT_GFS);\"\"\"\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=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:17.370246Z",
     "iopub.status.busy": "2021-09-28T19:19:17.369750Z",
     "iopub.status.idle": "2021-09-28T19:19:17.371819Z",
     "shell.execute_reply": "2021-09-28T19:19:17.371390Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_exact_solution_single_point():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc = \"Exact solution at a single point. params.time==0 corresponds to the initial data.\"\n",
    "    c_type = \"void\"\n",
    "    name = \"exact_solution_single_point\"\n",
    "    params = \"\"\"const paramstruct *restrict params,\n",
    "                const REAL xx0, const REAL xx1, const REAL xx2,\n",
    "                REAL *uu_exact, REAL *vv_exact\"\"\"\n",
    "    body = fin.FD_outputC(\"returnstring\",[lhrh(lhs=\"*uu_exact\",rhs=swid.uu_ID),\n",
    "                                          lhrh(lhs=\"*vv_exact\",rhs=swid.vv_ID)],\n",
    "                          params=\"includebraces=False,preindent=1,outCverbose=False\")\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"))\n",
    "# add_to_Cfunction_dict_exact_solution_single_point()\n",
    "# from outputC import outC_function_dict\n",
    "# print(outC_function_dict[\"exact_solution_single_point\"])\n",
    "# sys.exit(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:17.374939Z",
     "iopub.status.busy": "2021-09-28T19:19:17.374447Z",
     "iopub.status.idle": "2021-09-28T19:19:17.376149Z",
     "shell.execute_reply": "2021-09-28T19:19:17.375812Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_exact_solution_all_points():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc = \"Exact solution at all points. params.time==0 corresponds to the initial data.\"\n",
    "    c_type = \"void\"\n",
    "    name = \"exact_solution_all_points\"\n",
    "    params = \"const paramstruct *restrict params,REAL *restrict xx[3], REAL *restrict in_gfs\"\n",
    "    body = \"\"\"exact_solution_single_point(params, xx0, xx1, xx2,\n",
    "           &in_gfs[IDX4S(UUGF,i0,i1,i2)], &in_gfs[IDX4S(VVGF,i0,i1,i2)]);\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"), loopopts = \"AllPoints,Read_xxs\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:17.381402Z",
     "iopub.status.busy": "2021-09-28T19:19:17.381013Z",
     "iopub.status.idle": "2021-09-28T19:19:17.382918Z",
     "shell.execute_reply": "2021-09-28T19:19:17.382577Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_rhs_eval():\n",
    "    desc=\"Evaluate the scalar wave RHSs\"\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    if enable_FD_functions:\n",
    "        includes += [\"finite_difference_functions.h\"]\n",
    "    if enable_SIMD:\n",
    "        includes += [\"SIMD/SIMD_intrinsics.h\"]\n",
    "    c_type = \"void\"\n",
    "    name = \"rhs_eval\"\n",
    "    params   =\"const paramstruct *restrict params, \"\n",
    "    if enable_rfm_precompute:\n",
    "        params += \"const rfm_struct *restrict rfmstruct, \"\n",
    "    else:\n",
    "        params += \"REAL *xx[3], \"\n",
    "    params += \"const REAL *restrict in_gfs, REAL *restrict rhs_gfs\"\n",
    "    body = fin.FD_outputC(\"returnstring\",[lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"uu\"),rhs=swrhs.uu_rhs),\n",
    "                                          lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"vv\"),rhs=swrhs.vv_rhs)],\n",
    "                          params=\"enable_SIMD=\"+str(enable_SIMD))\n",
    "    loopopts = \"InteriorPoints\"\n",
    "    if enable_SIMD:\n",
    "        loopopts += \",enable_SIMD\"\n",
    "    if enable_rfm_precompute:\n",
    "        loopopts += \",enable_rfm_precompute\"\n",
    "    else:\n",
    "        loopopts += \",Read_xxs\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=os.path.join(\".\"), loopopts = loopopts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='OuterBoundaryConditions'></a>\n",
    "\n",
    "## Step 1.b: Output needed C code for boundary condition driver \\[Back to [top](#toc)\\]\n",
    "$$\\label{OuterBoundaryConditions}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:17.385704Z",
     "iopub.status.busy": "2021-09-28T19:19:17.385312Z",
     "iopub.status.idle": "2021-09-28T19:19:18.441868Z",
     "shell.execute_reply": "2021-09-28T19:19:18.442305Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evolved gridfunction \"uu\" has parity type 0.\n",
      "Evolved gridfunction \"vv\" has parity type 0.\n"
     ]
    }
   ],
   "source": [
    "import CurviBoundaryConditions.CurviBoundaryConditions as CBC\n",
    "CBC.CurviBoundaryConditions_register_NRPy_basic_defines()\n",
    "CBC.CurviBoundaryConditions_register_C_functions(radiation_BC_FD_order=radiation_BC_FD_order)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cparams_rfm_and_domainsize'></a>\n",
    "\n",
    "## Step 1.c: Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h` \\[Back to [top](#toc)\\]\n",
    "$$\\label{cparams_rfm_and_domainsize}$$\n",
    "\n",
    "Here we output `free_parameters.h`, which sets initial data parameters, as well as grid domain & reference metric parameters, applying `domain_size` and `sinh_width`/`SymTP_bScale` (if applicable) as set above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:18.445773Z",
     "iopub.status.busy": "2021-09-28T19:19:18.445259Z",
     "iopub.status.idle": "2021-09-28T19:19:18.447375Z",
     "shell.execute_reply": "2021-09-28T19:19:18.446930Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1.c.i: Set free_parameters.h\n",
    "outstr = r\"\"\"\n",
    "// Free parameters related to physical system:\n",
    "params.time = 0.0; // Initial simulation time corresponds to exact solution at time=0.\n",
    "params.wavespeed = 1.0;\n",
    "\n",
    "// Free parameters related to numerical timestep:\n",
    "REAL CFL_FACTOR = \"\"\"+str(CFL_FACTOR)+\";\\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",
    "outstr += rfm.out_default_free_parameters_for_rfm(\"returnstring\",\n",
    "                                                  domain_size,sinh_width,sinhv2_const_dr,SymTP_bScale)\n",
    "with open(os.path.join(Ccodesrootdir,\"free_parameters.h\"),\"w\") as file:\n",
    "    file.write(outstr.replace(\"params.\", \"griddata.params.\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cfl'></a>\n",
    "\n",
    "## Step 1.d: Output needed C code for finding the minimum proper distance between grid points, needed for [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673)-limited timestep \\[Back to [top](#toc)\\]\n",
    "$$\\label{cfl}$$\n",
    "\n",
    "In order for our explicit-timestepping numerical solution to the scalar wave equation to be stable, it must satisfy the [CFL](https://en.wikipedia.org/w/index.php?title=Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition&oldid=806430673) condition:\n",
    "$$\n",
    "\\Delta t \\le \\frac{\\min(ds_i)}{c},\n",
    "$$\n",
    "where $c$ is the wavespeed, and\n",
    "$$ds_i = h_i \\Delta x^i$$ \n",
    "is the proper distance between neighboring gridpoints in the $i$th direction (in 3D, there are 3 directions), $h_i$ is the $i$th reference metric scale factor, and $\\Delta x^i$ is the uniform grid spacing in the $i$th direction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:18.481660Z",
     "iopub.status.busy": "2021-09-28T19:19:18.468287Z",
     "iopub.status.idle": "2021-09-28T19:19:18.483454Z",
     "shell.execute_reply": "2021-09-28T19:19:18.482960Z"
    }
   },
   "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",
    "rfm.register_C_functions(enable_rfm_precompute=enable_rfm_precompute)\n",
    "rfm.register_NRPy_basic_defines(enable_rfm_precompute=enable_rfm_precompute)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mainc'></a>\n",
    "\n",
    "# Step 2: The C code `main()` function for `ScalarWaveCurvilinear_Playground` \\[Back to [top](#toc)\\]\n",
    "$$\\label{mainc}$$\n",
    "\n",
    "Just as in [the start-to-finish, solving the scalar wave equation in Cartesian coordinates module](Tutorial-Start_to_Finish-ScalarWave.ipynb), we will implement the scalar wave equation via the Method of Lines. As discussed above, the critical differences between this code and the Cartesian version are as follows:\n",
    "1. The CFL-constrained timestep depends on the proper distance between neighboring gridpoints\n",
    "1. The boundary conditions must account for the fact that ghost zone points lying in the domain exterior can map either to the interior of the domain or lie on the outer boundary. In the former case, we simply copy the data from the interior. In the latter case, we apply the usual outer boundary conditions.\n",
    "1. The numerical grids must be staggered to avoid direct evaluation of the equations on coordinate singularities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:18.489438Z",
     "iopub.status.busy": "2021-09-28T19:19:18.488908Z",
     "iopub.status.idle": "2021-09-28T19:19:18.490738Z",
     "shell.execute_reply": "2021-09-28T19:19:18.490254Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_main__ScalarWaveCurvilinear_Playground():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.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: Write test data to gridfunctions\n",
    "// Step 2: Overwrite all data in ghost zones with NaNs\n",
    "// Step 3: Apply curvilinear boundary conditions\n",
    "// Step 4: Print gridfunction data after curvilinear boundary conditions have been applied\n",
    "// Step 5: Free all allocated memory\n",
    "\"\"\"\n",
    "    c_type = \"int\"\n",
    "    name = \"main\"\n",
    "    params = \"int argc, const char *argv[]\"\n",
    "    body = r\"\"\"  griddata_struct griddata;\n",
    "  set_Cparameters_to_default(&griddata.params);\n",
    "  griddata.params.outer_bc_type = \"\"\"+outer_bc_type+r\"\"\";\n",
    "\n",
    "  // Step 0a: Read command-line input, error out if nonconformant\n",
    "  if(argc != 4 || atoi(argv[1]) < NGHOSTS || atoi(argv[2]) < NGHOSTS || atoi(argv[3]) < NGHOSTS) {\n",
    "    printf(\"Error: Expected one command-line argument: ./ScalarWaveCurvilinear_Playground Nx0 Nx1 Nx2,\\n\");\n",
    "    printf(\"where Nx[0,1,2] is the number of grid points in the 0, 1, and 2 directions.\\n\");\n",
    "    printf(\"Nx[] MUST BE larger than NGHOSTS (= %d)\\n\",NGHOSTS);\n",
    "    exit(1);\n",
    "  }\n",
    "  // Step 0b: Set up numerical grid structure, first in space...\n",
    "  const int Nxx[3] = { atoi(argv[1]), atoi(argv[2]), atoi(argv[3]) };\n",
    "  if(Nxx[0]%2 != 0 || Nxx[1]%2 != 0 || Nxx[2]%2 != 0) {\n",
    "    printf(\"Error: Cannot guarantee a proper cell-centered grid if number of grid cells not set to even number.\\n\");\n",
    "    printf(\"       For example, in case of angular directions, proper symmetry zones will not exist.\\n\");\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  // Step 0c: Set free parameters, overwriting Cparameters defaults\n",
    "  //          by hand or with command-line input, as desired.\n",
    "#include \"free_parameters.h\"\n",
    "\n",
    "  // Step 0d: Uniform coordinate grids are stored to *xx[3]\n",
    "  // Step 0d.i: Set bcstruct\n",
    "  {\n",
    "    int EigenCoord;\n",
    "    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, &griddata.params, griddata.xx);\n",
    "    // Step 0e: Find ghostzone mappings; set up bcstruct\n",
    "    bcstruct_set_up(&griddata.params, griddata.xx, &griddata.bcstruct);\n",
    "    // Step 0e.i: Free allocated space for xx[][] array\n",
    "    for(int i=0;i<3;i++) free(griddata.xx[i]);\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",
    "    EigenCoord = 0;\n",
    "    set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &griddata.params, griddata.xx);\n",
    "  }\n",
    "\n",
    "  // Step 0g: Time coordinate parameters\n",
    "  griddata.params.t_final =  0.7*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 0h: Set timestep based on smallest proper distance between gridpoints and CFL factor\n",
    "  griddata.params.dt = find_timestep(&griddata.params, griddata.xx, CFL_FACTOR);\n",
    "  // Ntot = N_outputs * output_every_N -> output_every_N = Ntot/N_outputs; set N_outputs=800\n",
    "  int output_every_N = (int)((griddata.params.t_final/griddata.params.dt + 0.5) / 800.0);\n",
    "  if(output_every_N == 0) output_every_N = 1;\n",
    "\n",
    "  // Step 0i: 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",
    "    printf(\"Error: NUM_AUX_GFS > NUM_EVOL_GFS. Either reduce the number of auxiliary gridfunctions,\\n\");\n",
    "    printf(\"       or allocate (malloc) by hand storage for *diagnostic_output_gfs. \\n\");\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  // Step 0j: Declare struct for gridfunctions and allocate memory for y_n_gfs gridfunctions\n",
    "  MoL_malloc_y_n_gfs(&griddata.params, &griddata.gridfuncs);\n",
    "\"\"\"\n",
    "    if enable_rfm_precompute:\n",
    "        body += \"\"\"\n",
    "  // Step 0k: Set up precomputed reference metric arrays\n",
    "  // Step 0k.i: Allocate space for precomputed reference metric arrays.\n",
    "  rfm_struct rfmstruct;\n",
    "  rfm_precompute_rfmstruct_malloc(&griddata.params, &griddata.rfmstruct);\n",
    "\n",
    "  // Step 0k.ii: Define precomputed reference metric arrays.\n",
    "  rfm_precompute_rfmstruct_define(&griddata.params, griddata.xx, &griddata.rfmstruct);\\n\"\"\"\n",
    "    body += r\"\"\"\n",
    "  // Step 0.l: Set up initial data to be exact solution at time=0:\n",
    "  griddata.params.time = 0.0; griddata.params.n_0 = 0;\n",
    "  exact_solution_all_points(&griddata.params, griddata.xx, griddata.gridfuncs.y_n_gfs);\n",
    "\n",
    "  // Step 0.m: Allocate memory for non y_n_gfs. We do this here to free up\n",
    "  //         memory for setting up initial data (for cases in which initial\n",
    "  //         data setup is memory intensive.)\n",
    "  MoL_malloc_non_y_n_gfs(&griddata.params, &griddata.gridfuncs);\n",
    "\n",
    "  while(griddata.params.time < griddata.params.t_final)\n",
    "    { // Main loop to progress forward in time.\n",
    "\n",
    "    // Step 1: Diagnostics/code validation: Compute log of L2 norm of difference\n",
    "    //         between numerical and exact solutions:\n",
    "    //   log_L2_Norm = log10( sqrt[Integral( [numerical - exact]^2 * dV)] ),\n",
    "    //         where integral is within 30% of the grid outer boundary (domain_size)\n",
    "    if(griddata.params.n%output_every_N == 0) {\n",
    "      REAL integral = 0.0;\n",
    "      REAL numpts   = 0.0;\n",
    "      const int Nxx_plus_2NGHOSTS0 = griddata.params.Nxx_plus_2NGHOSTS0;\n",
    "      const int Nxx_plus_2NGHOSTS1 = griddata.params.Nxx_plus_2NGHOSTS1;\n",
    "      const int Nxx_plus_2NGHOSTS2 = griddata.params.Nxx_plus_2NGHOSTS2;\n",
    "#pragma omp parallel for reduction(+:integral,numpts)\n",
    "      LOOP_REGION(NGHOSTS,Nxx_plus_2NGHOSTS0-NGHOSTS,\n",
    "                  NGHOSTS,Nxx_plus_2NGHOSTS1-NGHOSTS,\n",
    "                  NGHOSTS,Nxx_plus_2NGHOSTS2-NGHOSTS) {\n",
    "        REAL xCart[3]; xx_to_Cart(&griddata.params,griddata.xx,i0,i1,i2, xCart);\n",
    "        if(sqrt(xCart[0]*xCart[0] + xCart[1]*xCart[1] + xCart[2]*xCart[2]) < domain_size*0.3) {\n",
    "          REAL uu_exact,vv_exact; exact_solution_single_point(&griddata.params,\n",
    "                                                              griddata.xx[0][i0],griddata.xx[1][i1],griddata.xx[2][i2],\n",
    "                                                              &uu_exact,&vv_exact);\n",
    "          double num   = (double)griddata.gridfuncs.y_n_gfs[IDX4S(UUGF,i0,i1,i2)];\n",
    "          double exact = (double)uu_exact;\n",
    "          integral += (num - exact)*(num - exact);\n",
    "          numpts   += 1.0;\n",
    "        }\n",
    "      }\n",
    "      // Compute and output the log of the L2 norm.\n",
    "      REAL log_L2_Norm = log10(1e-16 + sqrt(integral/numpts));  // 1e-16 + ... avoids log10(0)\n",
    "      printf(\"%e %e\\n\",(double)griddata.params.time, log_L2_Norm);\n",
    "    }\n",
    "\n",
    "    // Step 2: Evolve scalar wave initial data forward in time using Method of Lines with RK4 algorithm,\n",
    "    //         applying quadratic extrapolation outer boundary conditions.\n",
    "    // Step 2.b: Step forward one timestep (t -> t+dt) in time using\n",
    "    //           chosen RK-like MoL timestepping algorithm\n",
    "    MoL_step_forward_in_time(&griddata);\n",
    "\n",
    "  } // End main loop to progress forward in time.\n",
    "\n",
    "  // Step 3: Free all allocated memory\n",
    "\"\"\"\n",
    "    if enable_rfm_precompute:\n",
    "        body += \"  rfm_precompute_rfmstruct_freemem(&griddata.params, &griddata.rfmstruct);\\n\"\n",
    "    body += r\"\"\"\n",
    "  free(griddata.bcstruct.inner_bc_array);\n",
    "  for(int ng=0;ng<NGHOSTS*3;ng++) free(griddata.bcstruct.pure_outer_bc_array[ng]);\n",
    "  MoL_free_memory_y_n_gfs(&griddata.params, &griddata.gridfuncs);\n",
    "  MoL_free_memory_non_y_n_gfs(&griddata.params, &griddata.gridfuncs);\n",
    "  for(int i=0;i<3;i++) free(griddata.xx[i]);\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",
    "\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": [
    "Next, register all remaining C functions in `outC_function_dict`, and output `finite_difference_functions.h`. Also construct `NRPy_basic_defines.h`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:18.558421Z",
     "iopub.status.busy": "2021-09-28T19:19:18.532882Z",
     "iopub.status.idle": "2021-09-28T19:19:18.631626Z",
     "shell.execute_reply": "2021-09-28T19:19:18.631141Z"
    }
   },
   "outputs": [],
   "source": [
    "# Register all functions needed for scalar wave:\n",
    "add_to_Cfunction_dict_exact_solution_single_point()\n",
    "add_to_Cfunction_dict_exact_solution_all_points()\n",
    "add_to_Cfunction_dict_rhs_eval()\n",
    "add_to_Cfunction_dict_main__ScalarWaveCurvilinear_Playground()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:18.558421Z",
     "iopub.status.busy": "2021-09-28T19:19:18.532882Z",
     "iopub.status.idle": "2021-09-28T19:19:18.631626Z",
     "shell.execute_reply": "2021-09-28T19:19:18.631141Z"
    }
   },
   "outputs": [],
   "source": [
    "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",
    "par.register_NRPy_basic_defines()  # add `paramstruct params` to griddata struct.\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=False,\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": [
    "<a id='compileexec'></a>\n",
    "\n",
    "# Step 3: Compile generated C codes & solve the scalar wave equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{compileexec}$$\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": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:18.635384Z",
     "iopub.status.busy": "2021-09-28T19:19:18.634720Z",
     "iopub.status.idle": "2021-09-28T19:19:20.300605Z",
     "shell.execute_reply": "2021-09-28T19:19:20.300104Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(EXEC): Executing `make -j18`...\n",
      "(BENCH): Finished executing in 1.61 seconds.\n",
      "Finished compilation.\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./ScalarWaveCurvilinear_Playground 16  2 4`...\n",
      "(BENCH): Finished executing in 0.20 seconds.\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./ScalarWaveCurvilinear_Playground 24  2 4`...\n",
      "(BENCH): Finished executing in 0.20 seconds.\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd\n",
    "cmd.new_C_compile(Ccodesrootdir, \"ScalarWaveCurvilinear_Playground\",\n",
    "                  uses_free_parameters_h=True, compiler_opt_option=\"fast\") # fastdebug or debug also supported\n",
    "os.chdir(Ccodesrootdir)\n",
    "if \"Spherical\" in CoordSystem:\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"16  2 4\", os.path.join(\"output\", \"out-lowresolution.txt\"))\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"24  2 4\", os.path.join(\"output\", \"out-medresolution.txt\"))\n",
    "elif \"Cartesian\" in CoordSystem:\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"16 16 16\", os.path.join(\"output\", \"out-lowresolution.txt\"))\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"24 24 24\", os.path.join(\"output\", \"out-medresolution.txt\"))\n",
    "elif \"Cylindrical\" in CoordSystem:\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"16 4 16\", os.path.join(\"output\", \"out-lowresolution.txt\"))\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"24 4 24\", os.path.join(\"output\", \"out-medresolution.txt\"))\n",
    "elif \"SymTP\" in CoordSystem:\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"16 16 2\", os.path.join(\"output\", \"out-lowresolution.txt\"))\n",
    "    cmd.Execute(\"ScalarWaveCurvilinear_Playground\", \"24 24 2\", os.path.join(\"output\", \"out-medresolution.txt\"))\n",
    "else:\n",
    "    print(\"Who ordered that? CoordSystem == \" + CoordSystem + \" not supported!\")\n",
    "    sys.exit(1)\n",
    "\n",
    "# Return to root directory\n",
    "os.chdir(os.path.join(\"..\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convergence'></a>\n",
    "\n",
    "# Step 4: Code validation: Plot the numerical error, and confirm that it converges to zero at expected rate with increasing numerical resolution (sampling) \\[Back to [top](#toc)\\]\n",
    "$$\\label{convergence}$$\n",
    "The numerical solution $u_{\\rm num}(x0,x1,x2,t)$ should converge to the exact solution $u_{\\rm exact}(x0,x1,x2,t)$ at fourth order, which means that\n",
    "$$\n",
    "u_{\\rm num}(x0,x1,x2,t) = u_{\\rm exact}(x0,x1,x2,t) + \\mathcal{O}\\left((\\Delta x0)^4\\right)+ \\mathcal{O}\\left((\\Delta x1)^4\\right)+ \\mathcal{O}\\left((\\Delta x2)^4\\right)+ \\mathcal{O}\\left((\\Delta t)^4\\right).\n",
    "$$\n",
    "\n",
    "Thus the relative error $E_{\\rm rel}$ should satisfy:\n",
    "$$\n",
    "\\left|\\frac{u_{\\rm num}(x0,x1,x2,t) - u_{\\rm exact}(x0,x1,x2,t)}{u_{\\rm exact}(x0,x1,x2,t)}\\right| + \\mathcal{O}\\left((\\Delta x0)^4\\right)+ \\mathcal{O}\\left((\\Delta x1)^4\\right)+ \\mathcal{O}\\left((\\Delta x2)^4\\right)+ \\mathcal{O}\\left((\\Delta t)^4\\right).\n",
    "$$\n",
    "\n",
    "We confirm this convergence behavior by first solving the scalar wave equation at two resolutions: $16\\times 8\\times 16$ (or $16^3$ if `reference_metric::CoordSystem` is set to `Cartesian`), and $24\\times 12\\times 24$ (or $24^3$ if `reference_metric::CoordSystem` is set to `Cartesian`) and evaluating the maximum logarithmic relative error $\\log_{10} E_{\\rm rel,max}$ between numerical and exact solutions within a region $R < 0.1 {\\rm RMAX}$ at all iterations. \n",
    "\n",
    "Since we increase the resolution uniformly over all four coordinates $(x0,x1,x2,t)$, $E_{\\rm rel}$ should drop uniformly as $(\\Delta x0)^4$:\n",
    "$$\n",
    "E_{\\rm rel} \\propto (\\Delta x0)^4.\n",
    "$$\n",
    "\n",
    "So at the two resolutions, we should find that\n",
    "$$\n",
    "\\frac{E_{\\rm rel}(16\\times 8\\times 16)}{E_{\\rm rel}(24\\times 12\\times 24)} = \\frac{E_{\\rm rel}(16^3)}{E_{\\rm rel}(24^3)} \\approx \\left(\\frac{(\\Delta x0)_{16}}{(\\Delta x0)_{24}}\\right)^{4} = \\left(\\frac{24}{16}\\right)^4 \\approx 5.\n",
    "$$\n",
    "\n",
    "Since we're measuring logarithmic relative error, this should be\n",
    "$$\n",
    "\\log_{10}\\left(\\frac{E_{\\rm rel}(16\\times 8\\times 16)}{E_{\\rm rel}(24\\times 12\\times 24)}\\right) = \\log_{10}\\left(\\frac{E_{\\rm rel}(16^3)}{E_{\\rm rel}(24^3)}\\right) \\approx \\log_{10}(5).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:20.358956Z",
     "iopub.status.busy": "2021-09-28T19:19:20.358248Z",
     "iopub.status.idle": "2021-09-28T19:19:20.721967Z",
     "shell.execute_reply": "2021-09-28T19:19:20.721719Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAHHCAYAAACfqw0dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAACQYUlEQVR4nOzdd1hT1xsH8G8WYcoSBBSZLkQEwT1A1LpHnXXjtmqdtdW2/hRr3W3VurVWa7XWra174d57gKICghuQTfb5/UFJjWEkISFB3s/z5FFubu59782997w595xzOYwxBkIIIYQQohWusQMghBBCCCmLKIkihBBCCNEBJVGEEEIIITqgJIoQQgghRAeURBFCCCGE6ICSKEIIIYQQHVASRQghhBCiA0qiCCGEEEJ0QEkUIYQQQogOKIkiKiIiIuDp6WnsMHTm6emJTp06GTuMAs2aNQscDsfYYZBypiTHXf5nk5OT9RyVZjgcDsaNG1cq64qKigKHw0FUVJRBlr9x40ZwOBzEx8cbZPnGUtB+K+vliDZMIonKP7jyX+bm5nBzc0Pbtm2xbNkyZGZmGjvEUrd161YsWbLEIMt+8eIFZs2ahVu3bhlk+foglUrh5+cHDoeDxYsXq7z34MEDzJo166O7GBnCrVu3MGDAALi7u0MoFMLBwQGtW7fGb7/9BrlcbuzwSAlkZWVh5syZ8Pf3h5WVFRwdHREYGIgJEybgxYsXpR6PRCLB0qVLERQUhAoVKsDOzg61a9fGyJEjERMTU+rxlHV07hbtwoULmDVrFtLS0owaB9+oa//A7Nmz4eXlBalUilevXiEqKgoTJ07ETz/9hP379yMgIMDYIZaarVu34t69e5g4caLel/3ixQtERkbC09MTgYGBKu+tW7cOCoVC7+vU1i+//IJnz54V+N6DBw8QGRmJsLCwcvNrRxfr16/H6NGjUalSJQwcOBDVqlVDZmYmTpw4gWHDhuHly5f45ptvjB0m0YFUKkWLFi0QExODwYMH44svvkBWVhbu37+PrVu34tNPP4WbmxsA4LvvvsO0adMMHlOPHj1w6NAh9O3bFyNGjIBUKkVMTAz++ecfNGnSBDVr1jR4DCXRokUL5ObmwszMzNihlPlztzTKkQsXLiAyMhIRERGws7Mz6LqKYlJJVPv27RESEqL8e/r06Th58iQ6deqELl26IDo6GhYWFkaM0DSJRCKYmZmByy15xaJAINBDRCXz5s0bzJ49G19//TX+97//GTscjclkMigUilK7COfk5MDS0rLA9y5duoTRo0ejcePGOHjwIGxsbJTvTZw4EdeuXcO9e/dKJU5Dyc7OhpWVlbHDMIq9e/fi5s2b2LJlC/r166fynkgkgkQiUf7N5/PB5xv2Un/16lX8888/+OGHH9QK9+XLlxu9tqAo718/zc3NjR2OyZ27jDGIRCKtyl5TKEdKi0nczitKeHg4ZsyYgYSEBPzxxx8q78XExKBnz55wcHCAubk5QkJCsH//fpV58m8Vnjt3DuPHj4eTkxPs7OwwatQoSCQSpKWlYdCgQbC3t4e9vT2++uorMMZUlpGdnY0pU6Yoq1Vr1KiBxYsXq82Xf/9+79698Pf3h1AoRO3atXH48GGV+TIzMzFx4kR4enpCKBTC2dkZbdq0wY0bNwAAYWFhOHDgABISEpS3OPNrXPLvP2/btg3fffcdKleuDEtLS2RkZCA1NRVffvkl6tSpA2tra1SoUAHt27fH7du3leuOiopC/fr1AQBDhgxRLn/jxo0A1O9lx8fHK2+prV27Fj4+PhAKhahfvz6uXr2q9n3t2LEDfn5+MDc3h7+/P/bs2aP1/fFp06ahRo0aGDBggNp7GzduRK9evQAALVu2VMb/YTuGc+fOoUGDBjA3N4e3tzd+//13jdf/5s0bDBs2DJUqVYK5uTnq1q2LTZs2qczz/n5ZsmSJcr88ePBAuf769evD3NwcPj4+WLNmTaHr++OPPxAcHAwLCws4ODjgs88+Q2Jioso8YWFh8Pf3x/Xr19GiRQtYWloW+Us0MjISHA4HW7ZsUbkI5wsJCUFERITyb30e4zt37gSHw8Hp06fV1rtmzRpwOByVQkCb8/j06dMYM2YMnJ2dUaVKFeX7K1asgLe3NywsLNCgQQOcPXsWYWFhCAsLU1mOWCzGzJkz4evrC6FQCHd3d3z11VcQi8Vab2e+58+fY9iwYXBzc4NQKISXlxc+//xzlUQmLS0NEydOVO5fX19fLFiwQO3X+suXLxETEwOpVKq2nvc9efIEANC0aVO198zNzVGhQgXl3wW1idJm+/Ljz//Fb2triyFDhiAnJ0ejeHg8HhwdHdXiiYmJQe/evVGhQgU4OjpiwoQJEIlEBa5f0+9h6NChqFSpknK+DRs2qMxT1PWzsDZRly9fRocOHWBvbw8rKysEBARg6dKlyvfv3LmDiIgIeHt7w9zcHC4uLhg6dChSUlIK3JbiGOrclclk+P7775XXKk9PT3zzzTdqx35+u9IjR44gJCQEFhYWyutXUlISunXrBisrKzg7O2PSpElqnwdKVo5osj9nzZqFqVOnAgC8vLyU5cD7TTw0ua7GxsaiR48ecHFxgbm5OapUqYLPPvsM6enpBXwzhWAm4LfffmMA2NWrVwt8PzExkQFgPXv2VE67d+8es7W1ZX5+fmzBggVs+fLlrEWLFozD4bDdu3erLTswMJC1a9eOrVixgg0cOJABYF999RVr1qwZ69evH1u5ciXr1KkTA8A2bdqk/LxCoWDh4eGMw+Gw4cOHs+XLl7POnTszAGzixIkqcQJgdevWZa6uruz7779nS5YsYd7e3szS0pIlJycr5+vXrx8zMzNjkydPZuvXr2cLFixgnTt3Zn/88QdjjLGjR4+ywMBAVrFiRbZ582a2efNmtmfPHsYYY6dOnWIAmJ+fHwsMDGQ//fQTmzdvHsvOzmZXr15lPj4+bNq0aWzNmjVs9uzZrHLlyszW1pY9f/6cMcbYq1ev2OzZsxkANnLkSOXynzx5whhjbPDgwczDw0MZa1xcHAPAgoKCmK+vL1uwYAFbuHAhq1ixIqtSpQqTSCTKef/55x/G4XBYQEAA++mnn9iMGTOYvb098/f3V1lmUS5fvsy4XC67cOGCct2LFi1Svv/kyRM2fvx4BoB98803yvhfvXrFGGPMw8OD1ahRg1WqVIl98803bPny5axevXqMw+Gwe/fuFbv+nJwcVqtWLSYQCNikSZPYsmXLWPPmzRkAtmTJErX94ufnx7y9vdn8+fPZzz//zBISEtidO3eYhYUFq1q1Kps3bx77/vvvWaVKlVhAQAD78JSbM2cO43A4rE+fPmzlypUsMjKSVaxYkXl6erJ3794p5wsNDWUuLi7MycmJffHFF2zNmjVs7969BW5DdnY2EwgELDw8XKN9ru9jPCcnh1lbW7MxY8aoratly5asdu3ayr+1PY/9/PxYaGgo++WXX9j8+fMZY4ytXLmSAWDNmzdny5YtY5MnT2YODg7Mx8eHhYaGKpchl8vZJ598wiwtLdnEiRPZmjVr2Lhx4xifz2ddu3bVejsZY+z58+fMzc1NuczVq1ezGTNmsFq1aim/v+zsbBYQEMAcHR3ZN998w1avXs0GDRrEOBwOmzBhgsp6Bw8ezACwuLi4Ir+zrVu3MgBs9uzZTKFQFDnvzJkz1Y47Tbcv/7NBQUGse/fubOXKlWz48OHK62e+CxcuMABsxIgRTCqVahRPnTp1WOfOndny5cvZgAEDGAA2cOBAneJ89eoVq1KlCnN3d2ezZ89mq1atYl26dGEA2M8//6ycr6jrZ/57p06dUs5/9OhRZmZmxjw8PNjMmTPZqlWr2Pjx41nr1q2V8yxevJg1b96czZ49m61du5ZNmDCBWVhYsAYNGqh8N/nHcFHfrSHP3fxjq2fPnmzFihVs0KBBDADr1q2bynweHh7M19eX2dvbs2nTprHVq1ezU6dOsZycHFa9enVmbm7OvvrqK7ZkyRIWHBysvK69v99KUo5osj9v377N+vbtq/x+88uBrKwsxphm11WxWMy8vLyYm5sbmzNnDlu/fj2LjIxk9evXZ/Hx8Rrtf8YYKxNJFGOM2drasqCgIOXfrVq1YnXq1GEikUg5TaFQsCZNmrBq1aqpLbtt27YqB3Tjxo0Zh8Nho0ePVk6TyWSsSpUqKhfevXv3MgBszpw5KvH07NmTcTgc9vjxY+U0AMzMzExl2u3btxkA9ssvv6hsy9ixY4vcJx07diww8cg/0b29vVlOTo7KeyKRiMnlcpVpcXFxTCgUstmzZyunXb16lQFgv/32m9ryCzv4HR0dWWpqqnL6vn37GAD2999/K6fVqVOHValShWVmZiqnRUVFMQAaJVEKhYI1aNCA9e3bV2Xd7ydRjDG2Y8cOtZM2n4eHBwPAzpw5o5z25s0bJhQK2ZQpU4qNYcmSJQyAMqFljDGJRMIaN27MrK2tWUZGhkpsFSpUYG/evFFZRrdu3Zi5uTlLSEhQTnvw4AHj8XgqhVl8fDzj8Xjshx9+UPn83bt3GZ/PV5keGhrKALDVq1cXuw35x9yHBXRhDHGM9+3blzk7OzOZTKac9vLlS8blclWORW3P42bNmqksUywWM0dHR1a/fn2Vwnvjxo0MgMq5vHnzZsblctnZs2dVtnP16tUMADt//rzW2zlo0CDG5XILvHblX2++//57ZmVlxR49eqTy/rRp0xiPx2PPnj1TTtM0icrJyWE1atRQnlsRERHs119/Za9fv1abt7AkSpPty//s0KFDVT7/6aefMkdHR5VtzT9GK1WqxPr27ctWrFihcg58uMwuXbqoTB8zZgwDwG7fvq11nMOGDWOurq4qiRVjjH322WfM1tZWea0s6vr5YRIlk8mYl5cX8/DwUPlBk7+9+T5cDmOM/fnnn2rXIU2SKEOdu7du3WIA2PDhw1Xm+/LLLxkAdvLkSeW0/Gvo4cOHVebNvzZu375dOS07O5v5+vpqnERpUo5ouj8XLVpU4P7U9Lp68+ZNBoDt2LFDbX3aMPnbefmsra2VvfRSU1Nx8uRJ9O7dG5mZmUhOTkZycjJSUlLQtm1bxMbG4vnz5yqfHzZsmEqVdsOGDcEYw7Bhw5TTeDweQkJC8PTpU+W0gwcPgsfjYfz48SrLmzJlChhjOHTokMr01q1bw8fHR/l3QEAAKlSooLJMOzs7XL58uUQ9aAYPHqx2j1ooFCrbRcnlcqSkpMDa2ho1atRQ3irUVZ8+fWBvb6/8u3nz5gCg3K4XL17g7t27GDRoEKytrZXzhYaGok6dOhqtY+PGjbh79y4WLFhQolj9/PyU8QGAk5MTatSoofIdFObgwYNwcXFB3759ldMEAgHGjx+PrKwstVtUPXr0gJOTk/JvuVyOI0eOoFu3bqhatapyeq1atdC2bVuVz+7evRsKhQK9e/dWHsPJyclwcXFBtWrVcOrUKZX5hUIhhgwZUuw2ZGRkAECBtwIK22Z9H+N9+vTBmzdvVG6N7Ny5EwqFAn369AGg23k8YsQI8Hg85d/Xrl1DSkoKRowYodLup3///irHK5B3q7lWrVqoWbOmyv4ODw8HALX9Xdx2KhQK7N27F507d1Zpy5kv/3qzY8cONG/eHPb29irrbd26NeRyOc6cOaP8zMaNG8EYK/b2t4WFBS5fvqy8pbFx40YMGzYMrq6u+OKLLwq8xfIhTb7HfKNHj1b5u3nz5khJSVEeaxwOB0eOHMGcOXNgb2+PP//8E2PHjoWHhwf69OlTYJuosWPHqvz9xRdfAMg7HrWJkzGGXbt2oXPnzmCMqezjtm3bIj09Xe36V9D180M3b95EXFwcJk6cqNZw+f2y5P3liEQiJCcno1GjRgCg9XXXUOdu/j6dPHmy2nwAcODAAZXpXl5eatergwcPwtXVFT179lROs7S0xMiRIzWKFSi+HAFKvj81va7a2toCAI4cOaJya1pbJtWwvChZWVlwdnYGADx+/BiMMcyYMQMzZswocP43b96gcuXKyr/fL9CA/3agu7u72vR3794p/05ISICbm5vaQV2rVi3l++/7cD0AYG9vr7LMhQsXYvDgwXB3d0dwcDA6dOiAQYMGwdvbu+CNL4CXl5faNIVCgaVLl2LlypWIi4tT6Qb7fpsEXXy4XfknQv525e8HX19ftc/6+voWe/BnZGRg+vTpmDp1qtp3UtJY8+PNj1Uul+Pt27cq7zs4OMDMzAwJCQmoVq2aWiP9wr7vD7+Ht2/fIjc3F9WqVVOLoUaNGioFRGxsLBhjBc4LqDfOrFy5skaN1vPbw2g6NIghjvF27drB1tYWf/31F1q1agUA+OuvvxAYGIjq1asD0O08/nB/F3bc8fl8tUQkNjYW0dHRKknvh+vSZjvfvn2LjIwM+Pv7F7i899d7584djderKVtbWyxcuBALFy5EQkICTpw4gcWLF2P58uWwtbXFnDlzivy8Jt9jYfO+f/7nH29CoRDffvstvv32W7x8+RKnT5/G0qVLsX37dggEArU2rR8e9z4+PuByuWpDl2jyPaSlpWHt2rVYu3Ztgdv64T4u6Pr5ofx2XsV9v6mpqYiMjMS2bdvU1qNV2xoY7txNSEgAl8tVO09cXFxgZ2dX7HUtfxm+vr5q7etq1KihUaxA8eUIUPL9qel11cvLC5MnT8ZPP/2ELVu2oHnz5ujSpQsGDBigzA80USaSqKSkJKSnpysPgPzGmF9++aVatpzvw4Pl/V+vxU1nHzTI00Zh63l/mb1790bz5s2xZ88eHD16FIsWLcKCBQuwe/dutG/fXqP1FPQrau7cuZgxYwaGDh2K77//Hg4ODuByuZg4cWKJu5tqsl0lsXjxYkgkEvTp00d5EU1KSgKQd4LFx8fDzc1NoySiuFgTExPVLhKnTp1Sa4SsiZL0FlUoFOBwODh06FCBMb9fo6fNunx9fcHn83H37l2dYyuKJseCUChEt27dsGfPHqxcuRKvX7/G+fPnMXfuXOU8upzHJd3fderUwU8//VTg+x8m7/o65hUKBdq0aYOvvvqqwPfzk8qS8PDwwNChQ/Hpp5/C29sbW7ZsKTaJ0mb7tN0Xrq6u+Oyzz9CjRw/Url0b27dvx8aNG4vsJVjYgKDFrTv/OBowYAAGDx5c4LwfDo+jz17evXv3xoULFzB16lQEBgbC2toaCoUC7dq10/q6a+hzV9NBVw3VC17T8rEk+1Ob6+qPP/6IiIgI7Nu3D0ePHsX48eMxb948XLp0SaXjSlHKRBK1efNmAFBeaPNrbAQCAVq3bm3QdXt4eOD48ePIzMxUyfbzB4/z8PDQabmurq4YM2YMxowZgzdv3qBevXr44YcflEmULiMM79y5Ey1btsSvv/6qMj0tLQ0VK1ZU/m2IUbPz98Pjx4/V3ito2oeePXuGd+/eoXbt2mrvzZ07F3PnzsXNmzcRGBhY4vhdXFxw7NgxlWl169YFkLcdd+7cgUKhUKmN0vT7dnJygoWFBWJjY9Xee/jwocrfPj4+YIzBy8tLLwVpPktLS4SHh+PkyZNITEwstmbPUMd4nz59sGnTJpw4cQLR0dFgjClv5QH6OY/fP+5atmypnC6TyRAfH69SePr4+OD27dto1aqVXs4BJycnVKhQodju5j4+PsjKyjL4tQrI+2Xv4+NjMsNXCAQCBAQEIDY2VnlLJV9sbKzKj5nHjx9DoVBoPfabk5MTbGxsIJfL9bqP828h3rt3r9Dlvnv3DidOnEBkZKTKcCwFnf+aMNS56+HhAYVCgdjYWGUtFQC8fv0aaWlpGp3jHh4euHfvHhhjKufPh9e1ktBmfxZ2Dmt7Xa1Tpw7q1KmD7777DhcuXEDTpk2xevXqYn+E5DP5NlEnT57E999/Dy8vL/Tv3x8A4OzsjLCwMKxZswYvX75U+8yHt2pKokOHDpDL5Vi+fLnK9J9//hkcDkfjmqN8crlcrUrS2dkZbm5uKu0YrKystK4K5vF4ar8Md+zYodauJH9sHX2O3eLm5gZ/f3/8/vvvyMrKUk4/ffq0Rr+qxo8fjz179qi88rvVRkREYM+ePcoLbknjNzc3R+vWrVVe+dXKHTp0wKtXr/DXX38p55fJZPjll19gbW2N0NDQIpfN4/HQtm1b7N27V2Ww0OjoaBw5ckRl3u7du4PH4yEyMlLte2OM6dxFGgBmzpwJxhgGDhyo8n3ku379unLYBn0f4/lat24NBwcH/PXXX/jrr7/QoEEDlUJTH+dxSEgIHB0dsW7dOshkMuX0LVu2qN2W6t27N54/f45169apLSc3NxfZ2dnabB64XC66deuGv//+G9euXVN7P/877d27Ny5evKj2/QN5x/D7cWs6xMHt27cLfBRLQkICHjx4oNUtFn2IjY0tcHDctLQ0XLx4Efb29mq3M1esWKHy9y+//AIAWh9vPB4PPXr0wK5duwpMHnUtD+rVqwcvLy8sWbJE7VqT/93m13R8eP6W5GkThjh3O3ToUGBc+bWyHTt2LDauDh064MWLF9i5c6dyWk5OTqG3UHWhzf4srBzQ9LqakZGhcu4BeQkVl8vVqE1hPpOqiTp06BBiYmIgk8nw+vVrnDx5EseOHYOHhwf279+vMhDaihUr0KxZM9SpUwcjRoyAt7c3Xr9+jYsXLyIpKUllbKSS6Ny5M1q2bIlvv/0W8fHxqFu3Lo4ePYp9+/Zh4sSJKg0eNZGZmYkqVaqgZ8+eqFu3LqytrXH8+HFcvXoVP/74o3K+4OBg/PXXX5g8eTLq168Pa2trdO7cuchld+rUCbNnz8aQIUPQpEkT3L17F1u2bFFra+Xj4wM7OzusXr0aNjY2sLKyQsOGDTVqJ1CUuXPnomvXrmjatCmGDBmCd+/eYfny5fD39y/wYvC+evXqoV69eirT8m/r1a5dG926dVNODwwMBI/Hw4IFC5Ceng6hUIjw8HBlm7mSGDlyJNasWYOIiAhcv34dnp6e2LlzJ86fP48lS5Zo1OAzMjIShw8fRvPmzTFmzBhlEla7dm3cuXNHOZ+Pjw/mzJmD6dOnIz4+Ht26dYONjQ3i4uKwZ88ejBw5El9++aVO29GkSROsWLECY8aMQc2aNVVGPY6KisL+/fuVv7T0fYznEwgE6N69O7Zt24bs7Gy1x/cAJT+PzczMMGvWLHzxxRcIDw9H7969ER8fj40bN8LHx0fl1+rAgQOxfft2jB49GqdOnULTpk0hl8sRExOD7du3K8fF0cbcuXNx9OhRhIaGYuTIkahVqxZevnyJHTt24Ny5c7Czs8PUqVOxf/9+dOrUCREREQgODkZ2djbu3r2LnTt3Ij4+XllTPH36dGzatAlxcXFF1sgcO3YMM2fORJcuXdCoUSNYW1vj6dOn2LBhA8RiMWbNmqXVdpTU7du30a9fP7Rv3x7NmzeHg4MDnj9/jk2bNuHFixdYsmSJ2q2VuLg4dOnSBe3atcPFixfxxx9/oF+/fspaYW3Mnz8fp06dQsOGDTFixAj4+fkhNTUVN27cwPHjx5Gamqr1MrlcLlatWoXOnTsjMDAQQ4YMgaurK2JiYnD//n0cOXIEFSpUQIsWLbBw4UJIpVJUrlwZR48eRVxcnNbry2eIc7du3boYPHgw1q5di7S0NISGhuLKlSvYtGkTunXrplKLW5gRI0Zg+fLlGDRoEK5fvw5XV1ds3ry50AF/daHN/gwODgYAfPvtt/jss88gEAjQuXNnja+rJ0+exLhx49CrVy9Ur14dMpkMmzdvViblGitR3z49ye/6mf8yMzNjLi4urE2bNmzp0qXKbuUfevLkCRs0aBBzcXFhAoGAVa5cmXXq1Int3LlTbdkfdkHO72b79u1blemDBw9mVlZWKtMyMzPZpEmTmJubGxMIBKxatWps0aJFauOzAChw6AIPDw82ePBgxlhel+ypU6eyunXrMhsbG2ZlZcXq1q3LVq5cqfKZrKws1q9fP2ZnZ6cyREB+N9yCumWKRCI2ZcoU5urqyiwsLFjTpk3ZxYsXWWhoqEpXb8byupb6+fkxPp+vMtxBYV1TPxxmIH97Z86cqTJt27ZtrGbNmkwoFDJ/f3+2f/9+1qNHD1azZk21zxenqHWvW7eOeXt7K4cNyO9e6+HhwTp27Kg2f0H7oDCvX79mQ4YMYRUrVmRmZmasTp06asNBFBUbY4ydPn2aBQcHMzMzM+bt7c1Wr15dYFdzxhjbtWsXa9asGbOysmJWVlasZs2abOzYsezhw4cq8b8/vpKmrl+/zvr166c8du3t7VmrVq3Ypk2bVIbD0Ocx/r5jx44xAIzD4bDExMQCYyzJeZxv2bJlzMPDgwmFQtagQQN2/vx5FhwczNq1a6cyn0QiYQsWLGC1a9dmQqGQ2dvbs+DgYBYZGcnS09N12s6EhAQ2aNAg5uTkxIRCIfP29mZjx45lYrFYOU9mZiabPn068/X1ZWZmZqxixYqsSZMmbPHixSpj5Gg6xMHTp0/Z//73P9aoUSPm7OzM+Hw+c3JyYh07dlTprs5Y4UMcaLJ9hV0nP+yu//r1azZ//nwWGhrKXF1dGZ/PZ/b29iw8PFzle3x/mQ8ePGA9e/ZkNjY2zN7eno0bN47l5ubqFGd+DGPHjmXu7u5MIBAwFxcX1qpVK7Z27VrlPEVdPwsaJ4oxxs6dO8fatGmjvF4HBASoDK+QlJTEPv30U2ZnZ8dsbW1Zr1692IsXL9Suj5oMcfA+fZ+7UqmURUZGMi8vLyYQCJi7uzubPn26yvAi+fu2oGsoY3nHepcuXZilpSWrWLEimzBhAjt8+LDGQxxoUo5ouj8Zyxs+pHLlyozL5art2+Kuq0+fPmVDhw5lPj4+zNzcnDk4OLCWLVuy48ePF7jtheH8uxGEGExgYCCcnJzU2iERYigKhQJOTk7o3r17gbfviPHMmjULkZGRePv2rUpbTULKIpNvE0XKDqlUqnaPOSoqCrdv39ap5xshmhCJRGptH37//XekpqbScUcIMSiTahNFyrbnz5+jdevWGDBgANzc3BATE4PVq1fDxcVFbbA+QvTl0qVLmDRpEnr16gVHR0fcuHEDv/76K/z9/ZXPWSSEEEOgJIrojb29PYKDg7F+/Xq8ffsWVlZW6NixI+bPn1/iwT4JKYynpyfc3d2xbNkypKamwsHBAYMGDcL8+fM1GleMEEJ0RW2iCCGEEEJ0QG2iCCGEEEJ0QEkUIYQQQogOylWbKIVCgRcvXsDGxsYgjz4hhBBCiP4xxpCZmQk3Nze1B8QbU7lKol68eFHss4gIIYQQYpoSExM1fjhwaShXSVT+IzsSExNRoUIFI0dDCCGEEE1kZGTA3d1do0dvlaZylUTl38KrUKECJVGEEEJIGWNqTXFM58YiIYQQQkgZQkkUIYQQQogOKIkihBBCCNEBJVGEEEIIITqgJIoQQgghRAeURBFCCCGE6ICSKEIIIYQQHVASRQghhBCiA0qiCCGEEEJ0QEkUIYQQQogOykQSFR8fj2HDhsHLywsWFhbw8fHBzJkzIZFIjB0aIYQQQsqpMvHsvJiYGCgUCqxZswa+vr64d+8eRowYgezsbCxevNjY4RFCCCGkHOIwxpixg9DFokWLsGrVKjx9+lTjz2RkZMDW1hbp6en0AGJCCCGkjDDV8rtM3M4rSHp6OhwcHIwdBiGEEELKqTJxO+9Djx8/xi+//FLsrTyxWAyxWKz8OyMjw9ChEUIIIaScMGpN1LRp08DhcIp8xcTEqHzm+fPnaNeuHXr16oURI0YUufx58+bB1tZW+XJ3dzfk5hBCCCGkHDFqm6i3b98iJSWlyHm8vb1hZmYGAHjx4gXCwsLQqFEjbNy4EVxu0TlgQTVR7u7uJndPlRBCCCGFM9U2UUa9nefk5AQnJyeN5n3+/DlatmyJ4OBg/Pbbb8UmUAAgFAohFApLGiYhhBBCiJoy0Sbq+fPnCAsLg4eHBxYvXoy3b98q33NxcTFiZIQQQggpr8pEEnXs2DE8fvwYjx8/RpUqVVTeK6MjNBBCCCGkjCsTQxxERESAMVbgixBCCCHEGMpEEkUIIYQQYmooiSKEEEII0QElUYQQQgghOqAkihBCCCFEB5REEUIIIYTogJIoQgghhBAdUBJFCCGEEKIDSqIIIYQQQnRASRQhhBBCiA4oiSKEEEII0QElUYQQQgghOqAkihBCCCFEB5REEUIIIYTogJIoQgghhBAdUBJFCCGEEKIDSqIIIYQQQnRASRQhhBBCiA4oiSKEEEII0QElUYQQQgghOqAkihBCCCFEB5REEUIIIYTogJIoQgghhBAdUBJFCCGEEKIDSqIIIYQQQnRASRQhhBBCiA4oiSKEEEII0QElUYQQQgghOqAkihBCCCFEB5REEUIIIYTogJIoQgghhBAdUBJFCCGEEKIDSqIIIYQQQnRASRQhhBBCiA4oiSKEEEII0QElUYQQQgghOqAkihBCCCFEB5REEUIIIYTogJIoQgghhBAdUBJFCCGEEKIDSqIIIYQQQnRASRQhhBBCiA4oiSKEEEII0QElUYQQQgghOqAkihBCCCFEB5REEUIIIYTogJIoQgghhBAdUBJFCCGEEKKDMpNEdenSBVWrVoW5uTlcXV0xcOBAvHjxwthhEUIIIaScKjNJVMuWLbF9+3Y8fPgQu3btwpMnT9CzZ09jh0UIIYSQcorDGGPGDkIX+/fvR7du3SAWiyEQCDT6TEZGBmxtbZGeno4KFSoYOEJCCCGE6IOplt9lpibqfampqdiyZQuaNGmicQJFCCGEEKJPZSqJ+vrrr2FlZQVHR0c8e/YM+/btK3J+sViMjIwMlRchhBBCiD4YNYmaNm0aOBxOka+YmBjl/FOnTsXNmzdx9OhR8Hg8DBo0CEXdjZw3bx5sbW2VL3d399LYLEIIIYSUA0ZtE/X27VukpKQUOY+3tzfMzMzUpiclJcHd3R0XLlxA48aNC/ysWCyGWCxW/p2RkQF3d3eTu6dKCCGEkMKZapsovjFX7uTkBCcnJ50+q1AoAEAlSfqQUCiEUCjUafmEEEIIIUUxahKlqcuXL+Pq1ato1qwZ7O3t8eTJE8yYMQM+Pj6F1kIRQgghhBhSmWhYbmlpid27d6NVq1aoUaMGhg0bhoCAAJw+fZpqmgghhBBiFGWiJqpOnTo4efKkscMghBBCCFEqEzVRhBBCCCGmhpIoQgghhBAdUBJFCCGEEKIDSqIIIYQQQnRASRQhhBBCiA4oiSKEEEII0QElUYQQQgghOqAkihBCCCFEB5REEUIIIYTogJIoQgghhBAdUBJFCCGEEKIDSqIIIYQQQnRASRQhhBBCiA4oiSKEEEII0QElUYQQQgghOqAkihBCCCFEB5REEUIIIYTogJIoQgghhBAd8LX9QFxcHM6ePYuEhATk5OTAyckJQUFBaNy4MczNzQ0RIyGEEEKIydE4idqyZQuWLl2Ka9euoVKlSnBzc4OFhQVSU1Px5MkTmJubo3///vj666/h4eFhyJgJIYQQQoxOoyQqKCgIZmZmiIiIwK5du+Du7q7yvlgsxsWLF7Ft2zaEhIRg5cqV6NWrl0ECJoQQQggxBRzGGCtupiNHjqBt27YaLTAlJQXx8fEIDg4ucXD6lpGRAVtbW6Snp6NChQrGDocQQgghGjDV8lujhuX5CZRMJsPvv/+O169fFzqvo6OjSSZQhBBCCCH6pFXvPD6fj9GjR0MkEhkqHkIIIYSQMkHrIQ4aNGiAW7duGSAUQgghhJCyQ+shDsaMGYPJkycjMTERwcHBsLKyUnk/ICBAb8ERQgghhJgqjRqWv4/LVa+84nA4YIyBw+FALpfrLTh9M9WGaYQQQggpnKmW3zoNtkkIIYQQUt5pnUTRQJqEEEIIITokUQDw5MkTLFmyBNHR0QAAPz8/TJgwAT4+PnoNjhBCCCHEVGndO+/IkSPw8/PDlStXEBAQgICAAFy+fBm1a9fGsWPHDBEjIYQQQojJ0bpheVBQENq2bYv58+erTJ82bRqOHj2KGzdu6DVAfTLVhmmEEEIIKZyplt9a10RFR0dj2LBhatOHDh2KBw8e6CUoQgghhBBTp3US5eTkVOBgm7du3YKzs7M+YiKEEEIIMXlaNywfMWIERo4ciadPn6JJkyYAgPPnz2PBggWYPHmy3gMkhBBCCDFFWreJYoxhyZIl+PHHH/HixQsAgJubG6ZOnYrx48eDw+EYJFB9MNV7qoQQQggpnKmW31rVRMlkMmzduhX9+vXDpEmTkJmZCQCwsbExSHCEEEIIIaZKqzZRfD4fo0ePhkgkApCXPFECRQghhJDySOuG5Q0aNMDNmzcNEQshhBBCSJmhdcPyMWPGYMqUKUhKSkJwcDCsrKxU3g8ICNBbcIQQQgghpkrrhuVcrnrlFYfDAWMMHA4Hcrlcb8Hpm6k2TCOEEEJI4Uy1/Na6JiouLs4QcRBCCCGElClaJVFSqRTh4eH4559/UKtWLUPFRAghpAwSiUR4+/Yt3r59izdv3ij/n5mZiZycHOUrOzsbOTk5kMvl4HA4ai+hUAhra2tl56X8l52dHVxcXODq6gpXV1dUqFDBpIfVIR8/rZIogUCg7JlHCCGkfBGJRIiNjUVcXJzaKyEhARkZGaUaj4WFBVxdXVG5cmX4+vqievXqypePjw8sLCxKNR5S/mjdJmru3Ll49OgR1q9fDz5f67uBRmWq91QJIcSUKBQKPHnyBPfu3cPdu3dx9+5d3Lt3D7GxscW2e+XxeGCMKV8f6tGjB1q0aAELCwvExMTgp59+KnRZ4eHhqFWrFjIzM/HixQtcvXoVCoUCYrEYEomkyDg4HA6qVq2KwMBABAcHIyQkBMHBwfR4sjLKVMtvrbOgq1ev4sSJEzh69Cjq1Kmj1jtv9+7deguOEEKI4aWmpuLy5cu4ePEiLl26hMuXL2tVqzR58mSMGDECLi4uuHjxIjp06KA2j7m5OSwtLdGpUydEREQAAM6cOYOjR49CLBYrXxKJRPn/Hj16YMyYMQCAI0eOoF27doXGEBoaiipVquDRo0eIiYlBZmYmEhISkJCQgH379innc3d3R3BwMJo2bYrw8HDUrVsXPB5P420l5H1a10QNGTKkyPd/++23EgVkSKaayRJCSGlKTk7GiRMncOzYMZw7dw4PHz5Um0cgECAgIAD+/v4wMzPDunXrVN6vWLEi3NzcULlyZYwbN06ZOGVkZODp06ewt7eHpaUlLC0tYWFhUWDP7uLk9/oGgFevXuHMmTN49eoVXr58iVevXilfSUlJWLFiBXr37g0AOHDgADp16gQgr0Yq/8d+VlaW2jrs7e3RsmVLhIeHIzw8HDVr1qR2VibIVMtvrZOossxUvwRCCDEksViMCxcu4OjRozh27Bhu3LhR4K22902fPh1z584FALx9+xZr166Fv78//Pz84O7uDnNz89IIXWPvJ1zHjh3DnDlzcP/+faSkpKjN26NHD4jFYpw+fVr5+LJ87u7u+PTTT9GjRw80bdqUaqm0cPbsWZw9exbTpk3TKWkuiqmW3zolUTKZDFFRUXjy5An69esHGxsbvHjxAhUqVIC1tbUh4lQSi8Vo2LAhbt++jZs3byIwMFDjz5rql0AIIfqWlZWFgwcPYteuXThw4ACys7NV3q9evToePXqk/NvR0RFNmzZV1j41bNgQnp6epRy1fjHG8ObNG9y/fx/37t3DuXPnEBUVhXPnzqF69eqQyWSYNm0ali1bBhsbG2RkZEAmkyk/7+TkhG7duqF79+4IDw+HmZmZEbfGtL158wZBQUF48eIFfvzxR0yePFmvyzfV8lvrJCohIQHt2rXDs2fPIBaL8ejRI3h7e2PChAkQi8VYvXq1oWIFAEyYMAGxsbE4dOgQJVGEEPKed+/e4Z9//sGuXbtw+PBhiMVilferVq2KOXPmoHXr1nBxccHo0aMRGBiI0NBQ1KxZU++1B6Yov8jLr7WKiIjApk2bVOaxsrKCTCZT2X/29vYYOHAghg8fjjp16pRewGWAXC5Hhw4dcPToUdSqVQtXr15Vay9dUiZbfjMtde3alQ0YMICJxWJmbW3Nnjx5whhj7NSpU8zX11fbxWnl4MGDrGbNmuz+/fsMALt586ZWn09PT2cAWHp6umECJISQUiaVStmBAwdYz549GZ/PZwAKfHl6erJJkyYZO1yTo1Ao2L1799jy5ctZjx49mLW1tXKfcblcNnjwYFapUiWVfdmwYUO2bt06lpGRYezwTcLs2bMZAGZpacnu3btnkHWYavmtdRLl4ODAYmJiGGNMJYmKi4tjFhYW+o3uPa9evWKVK1dmV69eZXFxcRolUSKRiKWnpytfiYmJJvklEEKIth49esSmT5+uVsBbWFgo/1+3bl0WGRnJ7t69yxQKhbFDLhOys7PZ9u3bWffu3VmXLl0YY4zJZDJ28OBB5u3tzbhcrnL/WllZseHDh7P79+8bOWrjOXHiBONwOAwA27hxo8HWY6pJlNZDHCgUigLHCUlKSoKNjY22i9MIYwwREREYPXo0QkJCEB8fr9Hn5s2bh8jISIPERAghpU0mk2HXrl1Yvnw5zp07p5xub2+PwYMHY8iQIbh//z4SExPRvXt3+Pr6GjHassnS0hK9evVCr169lLf+eDwe6tSpg7i4OGUDdisrK2RlZWH9+vVYv349unTpgq+//hpNmjQx8haUnpcvX6Jfv35gjGHo0KEYPHiwsUMqfdpmXb1792YjRoxgjOXVRD19+pRlZmay8PBwFhERodWyvv7660KrnvNf0dHRbOnSpaxp06ZMJpMxxhjVRBFCypW0tDT2448/ssqVK6tdIzkcDtu3b5+xQ/zopaWlsRUrVrDmzZur7H8nJyeVv5s3b87++eefj77mTyqVsrCwMAaA1alTh2VnZxt0faZaE6V1w/KkpCS0bdsWjDHExsYiJCQEsbGxqFixIs6cOaPVaLBv374tsPvp+7y9vdG7d2/8/fffKmN3yOVy8Hg89O/fX61RYGFMtmEaIYQUICEhAcuWLcOaNWvUetc5ODhgxIgRGDVqFLy8vIwUYfl07949LFq0CFu3blX25gsNDcWFCxcglUoBAP7+/vj+++/RtWvXj3LcqRkzZmDOnDmwtrbGtWvXUKNGDYOuz1TLb52HOPjrr79w+/ZtZGVloV69eujfv7/BnlP07NkzldFzX7x4gbZt22Lnzp1o2LAhqlSpotFyTPVLIISQ9z19+hSzZs3C1q1b1ZpPNG7cGGPHjkXPnj0hFAqNFCEB8sqmn3/+GXv37sX9+/eRlpaGJUuWYMWKFcjJyQEANGvWDIsXL0bDhg2NHK3+HDlyBO3btwdjDFu3bkXfvn0Nvk5TLb/L5GCb8fHx8PLyoiEOCCEflefPn2PWrFnYsGEDFAoFAKBVq1aYMmUKEhIS0KRJEwQEBBg5SvIhmUymfJasQqFArVq1kJSUBIlEoqyp6t27N+bNmwdvb29jhlpiz58/R2BgIJKTkzF69GisWrWqVNZrquV32XqCMCGEfISSk5Pxww8/YPny5SqDPf7+++8YOHCgESMjmshPoAAok6f8migHBwekpqZi+/bt2LNnD8aOHYsZM2bAwcHBWOHqTCaToV+/fkhOTkZgYCB+/vlnY4dkdGVyZDVPT08wxrSqhSKEEFMjEokwa9YsVKlSBUuWLFEmUM7Ozli/fn2p3CYh+lW1alXExsZizZo1sLe3R2pqKjgcDqpUqQKpVIolS5agVq1a2L59e7GP3jE1s2fPxpkzZ2BtbY3t27eb3KN/jKFMJlGEEFLWHTp0CNWrV0dkZKRyZGw7OzssWbIECQkJGDZsmEoNByk7+Hw+Ro4ciYcPH2Lw4MFgjCEpKQl2dnbw8fHBmzdv0KdPH3Tr1g3Pnz83drgaOX78OObMmQMAWLt2LapVq2bkiEwDJVGEEFKKEhMT0aNHD3To0AGJiYngcrmwsrLCDz/8gMTEREyYMIF+4X8knJycsHHjRkRFRaFWrVoIDAzEvXv3MHPmTAgEAuzfvx9+fn5YvXq1sg2cKXr16hUGDBgAxhiGDx9ONaTv0alheVpaGnbu3IknT55g6tSpcHBwwI0bN1CpUiVUrlzZEHHqhak2TCOEfPwkEgnmz5+POXPmQCqVgsfjYcKECejYsSPq1asHOzs7Y4dIDEgikeDdu3eoVKkSAODixYvo3bs3kpKSAADNmzfH+vXrUb16dWOGqUYul6Nt27Y4ceIE/P39cfnyZbx8+RI+Pj6lGoeplt9a10TduXMH1atXx4IFC7B48WKkpaUBAHbv3o3p06frOz5CCCnzbt26hWrVqmHmzJmQSqXw9fXFzZs38eOPPyI8PJwSqHLAzMxMmUABwJ9//omkpCQ4OzvDwsICZ8+eRVBQEDZs2GBSbaXmzZuHEydOwNLSEn/99RcuXryIWrVqYcaMGSYVp7FonURNnjwZERERiI2NValy7tChA86cOaPX4AghpCxTKBSYPXs26tWrh2fPngHIu8Wzdu1a1KlTx8jREWMKCwuDm5sb3rx5Ax6Ph4CAAOTk5GDYsGHo168f0tPTjR0izpw5g5kzZwIAVqxYAYVCge7du0MqleLx48eUREGHJOrq1asYNWqU2vTKlSvj1atXegmKEELKuufPnyMoKAgzZ85UFjYjR45EXFwcWrZsaeToiLF1794dN2/eRIsWLZCVlYU7d+6gdevW4HK52LZtG4KCgnDlyhWjxffq1Sv07dsXCoUCgwYNQps2bdC+fXtkZGSgRYsW2LhxI7hcalat9R4QCoUqo4fne/ToEZycnPQSFCGElGW7du2Cr68v7ty5AwCoVKkSzp49izVr1sDKysrI0RFT4ezsjOPHj2Ps2LEA8nrANW/eHO7u7oiLi0PTpk2xaNGiUm90LhKJ8Omnn+LFixeoWbMm5s+fj44dOyIpKQk1atTAnj17aLT8f2mdRHXp0gWzZ89WPh+Iw+Hg2bNn+Prrr9GjRw+9B0gIIWVF/u2Ynj17QiQSgcPhYPjw4YiLi0OzZs2MHR4xQQKBAMuXL8f69ethZmaGBw8e4ODBg+jVqxdkMhm++uordOjQAcnJyaUSD2MMo0aNwqVLl2BnZ4edO3diyJAhuH37NpydnXHo0KEyOVCowWj7xOK0tDTWunVrZmdnx3g8HnN3d2cCgYC1aNGCZWVllfiJyIZkqk+BJoSUfc+ePWN+fn4MAONwOGz69OksLi7O2GGRMuTixYvszJkzjDHGFAoFW7t2LbOwsGAAWNWqVdmVK1cMHsPChQsZAMbj8djRo0dZbm4u69GjB7O0tCyV9RfGVMtvnZ+dd+7cOdy5c0f5AOLWrVvrMbUzDFPtIkkIKdvOnDmDtm3bQiQSwdbWFrt370Z4eLixwyJl3N69exEVFYV//vkHT548gZmZGZYtW4aRI0eCw+HofX3//PMPunTpAsYYfvnlF4wbNw5AXgeJ+/fvG7UzhKmW31onUYmJiXB3dzdUPAZlql8CIaTsWrFiBb744gswxsDhcLBs2TJl4UOIrl68eIEaNWogKysLn376KWQyGf7++28AwODBg7Fq1SpYWFjobX337t1D48aNkZWVhVGjRmHkyJEICgoySLKmC1Mtv7VuE+Xp6YnQ0FCsW7cO7969M0RMhBBi8hQKBUaMGIFx48aBMQYzMzMcOXKEEiiiF25ubli8eDH4fD727NmDtLQ0zJo1C1wuF5s2bULjxo3x5MkTvawrOTkZXbp0QVZWFsLCwtCjRw80atQIAwYMgEQi0cs6PlZaP5jp2rVr2Lp1K2bPno0vvvgC7dq1w4ABA9C5c2dqrU+InjHG8PbtWyQmJiIxMRFJSUlITU1FWlqa2kskEkGhUEChUEAulyv/z+PxYGNjo3xVqFABNjY2sLe3h7u7O6pWrap8OTo6mswvT1OWk5ODVq1a4dKlSwCAihUr4vLly/D29jZyZORjMmrUKNSsWRNdunTB2bNnkZGRgb/++gtjxozB7du3ERwcjJ9++glDhgzR+bzNzc1Fz549ERcXB29vb8yZMwcdO3aEVCqFVCql5zcWQ+c2UYwxREVFYevWrdi1a5dyEK4NGzboO0a9MdXqQFK+KRQKJCUlITo6GtHR0YiJiUF0dDTi4+Px4sULyGSyUovFwsICHh4eCAgIQL169VCvXj0EBQWhYsWKpRaDqcvKykLTpk2VwxcEBATg3LlzsLGxMXJk5GN169YttGvXDq9fv4aPjw9+//13TJ06FRcuXAAAtGnTBmvXroWnp6dWy7127RoGDhyImJgY2NjYYOvWrfj888+RlJSEpk2b4vjx4ybzHEeTLb/10Tr9+vXrLDAwkHG5XH0szmBMtXU/KT9EIhG7evUqW7VqFRs2bBgLCAhg5ubmDIBWr88++4zNmzePrVq1ig0cOLDIeVetWsWOHDnCdu7cyfr166f1ugAwd3d31rVrV/bzzz+z6OhoplAojL0rjSI9PZ01bdpU2XupT58+TCaTGTssUg48fvyYeXl5MQDsyy+/ZFKplC1atEh5/bCysmLLli1jcrm82GVJJBI2a9YsxuPxGADm4uLC1qxZwypWrMgAsFq1arHk5ORS2CrNmWr5rXNNVFJSErZu3YqtW7cqG6T1798fo0eP1mVxpcJkM1nyUWKMIT4+HmfOnMHFixdx4cIFPHjwAHK5vMD5q1atipCQENSsWRNv377Fb7/9BicnJ7i4uMDFxQWVKlVSvnr06IGqVasCAF6/fo0XL16Ay+WCy+WCx+Op/L9KlSrKW+1XrlzB4cOH8ebNG7x9+1b5b1JSEtLT0/H3339DIBDg1q1b2L59O27cuFFgrJ6enmjXrh3atWuH8PDwclEL8+zZM/Tq1QtXrlyBnZ0dDhw4gCZNmhg7LFKOvHz5Ej///DPmzp2rvM326NEjDB8+HGfPngUANGvWDOvXr0eNGjUKXEZMTAwGDhyIa9euAQB69eqFgQMHon///sjMzERwcDAOHz5scrXPJlt+a5t1rV69mrVo0YLxeDxWu3ZtNnfuXBYfH6/37M4QTDWTJR8HhULBYmJi2Jo1a1ifPn2Ys7OzRrU8lSpVYq1atVKOD8MYYzKZrFRrexQKBYuLi2NSqVQ5bdy4cQXGy+FwVP4WCASsU6dObO/evUwikZRazKUpLi6OWVpaMgDM3t6eXb9+3dghEcKkUim7ePEik8vlbMWKFcza2poBYFwul3l4eLDQ0FAWERHBIiMj2aZNm9iCBQuUNVd2dnZs69atTKFQsKNHjzIzMzMWGhpqsuWjqZbfWtdEubu7o2/fvujfvz/q1q1b8iyuFJlsJkvKrLS0NBw7dgz79+/HgQMH1HqscrlcNGzYEM2bN4ezszO2bduGFi1aoE6dOqhVqxZq1qwJW1tbI0VfNMYYXr58iRs3buDGjRu4dOkSTpw4AYlEAh6Ph4iICJw6dQpPnz5VfqZSpUqIiIjAsGHDUK1aNSNGrz8JCQmoXbs2srOzweFwcOjQIbRt29bYYZFyjjGG4cOHY9OmTdiwYQMGDRqEhIQEjB49GocPHy7ys5988gk2bNiAypUrK6edPXsWISEheh02QZ9MtfzWOoli/46FUhaZ6pdAyg7GGG7fvo2DBw/i0KFDuHDhQqHPtXJwcMB3332HSZMmlXKUhpOWlob9+/fj6dOnmDVrFgAgOjoabdq0wevXr1Uawbdo0QKff/45evXqBR6PZ6SISyYxMRF+fn7IysoCh8PBjh076PFWxCTIZDIMGzYMv//+OwDgxx9/xOTJkwHk3fZ7+vQp4uPjERcXp/w3LS0Nw4YNw+eff461a9ciNDQUNWvWNOZmaMxUy2+Nkqg7d+7A398fXC5X2SOlMAEBAXoLTt9M9Usgpk2hUODy5cv466+/sG3bNrx+/brA+VxdXdG5c2eEhoaiSZMm8PDwKLM/OLSRkpICV1dXledpvn9ZqV27NmbPno1PP/20TO2P1NRUeHt7Iz09HRwOB1u2bEHfvn2NHRYhSgqFAlOnTsVPP/0EAKhVqxbatm2Ltm3bIjQ0tMBaJcYYZs+ejVmzZsHd3R23bt0qE8/CM9nyW6N7fhwOe/36tfL/XC6XcTgc5Sv/b+qdRz4WcrmcRUVFsZEjRzJ7e3uVNkBcLpd17tyZrVy5ksXFxbE///yTPXz4sNz2WGOMsbdv37J169axsLAw5X7i8XhMKBQq/65Xrx47cOBAmdhPubm5zNXVVRn7hg0bjB0SIQVSKBRswYIFjM/nq1yn3m9jmZ2dzRQKBZPL5Wz8+PHKeWbPnl0mzkfGTLf81iiJio+PV+7o+Pj4Il+mzFS/BGI67t69y77++mvm6OhYYKPqihUrspEjR5aZC48xnDlzhjVr1ky5zzp37qxs8AqANW7cmJ06dcrYYRZKLpeztm3bKuNdtGiRsUMipFipqalsx44dbPjw4ax27doqnTzGjBnDKleuzJo0aaI8rpctW2bEaLVnquW31r3zTp8+rdKDJ59UKmWnT5/WS1CGYqpfAjGuFy9esEWLFrHAwMACE6cqVaqw7777jt29e5eSJw0pFAp26NAh1rt3byaVStnbt2/Z1KlTVWqmhgwZwt69e2fsUNVMnjxZWeP45ZdfGjscQkrM399fpYZ48+bNxg5Ja6ZafmvdsJzH4+Hly5dwdnZWmZ6SkgJnZ+dCx8AxBSZ7T5WUOrlcjkOHDmHhwoXK8VUAQCAQoEOHDvDz84NIJMLAgQMRGBhYptrymCqxWAxvb2+kpaUhJycHQN7zwdauXYuOHTsaObo88+fPx/Tp0wEAmzdvxoABA4wcESEll5ubi7Nnz+LMmTMIDw9HeHi4sUPSmqmW31o/FIcV0jsvJSUFVlZWegmKEENJTEzE6tWrsWrVKrXhCLp164b169fD0dHRSNF93J4+fQqBQKBMoCwtLfHixQt06tQJgwYNwpIlS2Bvb2+0+D7//HOsXr0aALBgwQJKoMhHw8LCAp988gk++eQTY4fy0dG4Jqp79+4AgH379qFdu3YqDxuWy+W4c+cOatSoUez4FMZkqpksMSzGGI4cOYKff/4ZR48eVXmPy+Xik08+wYQJE9CqVSsIBAIjRVk+iMVirFmzBtOmTUNubi6srKyQk5MDxhhcXV2xevVqdOnSpdTjmj17NmbOnAkAaNCgAS5dukS1j4SYEFMtvzWuicofEJAxBhsbG5Wuk2ZmZmjUqBFGjBih/wgJ0VFOTg42b96MpUuXIjo6WuU9JycnfPHFFxg5ciQqVapkpAjLH6FQiPHjx6NVq1bo3bs3Hjx4AC6XC0dHR7x8+RJdu3bF+PHjsXjx4lJLaH///XdlAuXp6YkLFy5QAkUI0YjWbaIiIyPx5Zdflslbd6aayRL9ev78ORYtWoQ1a9ZAJBIBAGxsbDB06FDY29ujQYMGaNu2LbhcrpEjLd+ys7MxduxYbNq0CZ06dULNmjWxePFiAEBoaCi2b9+u1vZS386cOYOWLVtCoVDAwcEBCQkJsLa2Nug6CSHaM9XyW+cHEJdFpvolEP2Ijo7Gd999hz179qgM9jhhwgTMnj2bvnMTtXXrVrRr1w4ODg7Yt28fBgwYgKysLFSpUgV79uxBSEiIQdb7+PFj+Pn5QSqVwtzcHLGxsahSpYpB1kUIKRlTLb91SqJ27tyJ7du349mzZ5BIJCrvFfbUd1Ngql8CKZnr169j2rRpOH78uMp0b29vzJkzBz179qS2TmUEYwyffvopzp8/j+TkZAiFQqxevRoRERF6XY9EIoGHhwdevXoFHo+HS5cuGSxZI4SUnKmW31rfz1i2bBmGDBmCSpUq4ebNm2jQoAEcHR3x9OlTtG/f3hAxElKgs2fPol27dggJCVFJoOrXr49jx47h8ePH6Nu3LyVQZUhUVBT27duH5ORkODk5QSwWY8iQIfjiiy+Uj5UpKfbvg1vzE6gdO3ZQAkUI0Y22A0vVqFGDbd26lTHGmLW1NXvy5AljjLEZM2awsWPHlmjQKkMz1cG6iHZOnz7NGjVqpPIYFltbW9amTRt25coVY4dHSkChULANGzYwCwsLBoDZ2toqv+cmTZqwpKSkEq/j22+/VQ46ePjwYT1ETQgxNFMtv7WuiXr27BmaNGkCIG/siczMTADAwIED8eeff+onsyOkANevX0ejRo0QGhqKS5cuQSAQYOTIkYiNjcXLly9x9OhR1K9f39hhkhLgcDgYMmQIrly5gurVqyM9PR18Ph8WFha4cOECgoKCcOLECZ2X37NnT/zwww8AgLVr16Jt27b6Cp0QUg5pnUS5uLggNTUVAFC1alVcunQJABAXF6fSmJcQfXnw4AHatGmDkJAQXL58GUDeyPlbtmzBmjVr4O3tXeDTyknZ5e/vjytXrqBr166QyWTIzc2Fs7Mz3r59i08++QRz586FQqHQeHmZmZlo2bIldu3aBQAYPHgwhg4daqjwCSHlhNZJVHh4OPbv3w8AGDJkCCZNmoQ2bdqgT58++PTTT/UeICm/kpKS0LNnT9SuXVvZ5onD4aB3796Ij49Hr169jBwhMSRbW1vs3r0bc+bMgUAgwMaNGzFs2DAoFAp8++236Nq1q9qo8wU5f/48atSogaioKABASEgINmzYYODoCSHlgda98xQKBRQKBfj8vHE6t23bhgsXLqBatWoYNWoUzMzMDBKoPphq636iKjs7G4sWLcKCBQuU4zwBQJs2bbBs2TLUrFnTiNERY3j27BmqVq0KANiwYQPGjBkDsVgMT09PREZGokmTJvDx8VEZJFMqlSIyMhJz585V1pI3b94cp06dAo/HM8p2EEJ0Y6rlN40TRUyGQqHAH3/8gW+++QbPnz8HADg6OsLLywsrV66k9k4EAPDw4UM0aNAAfD5f2bQAACpWrIjGjRujcePG8PPzwzfffIMHDx4o369fvz5Onz5Nt34JKYNMtfzW6LEvd+7c0XiBAQEBOgdDyq/z588jIiICjx8/BpD3+I1Fixahffv2ZXJ0fGI4W7ZsQUZGBgCgVq1aYIzh0aNHSE5Oxt9//42///5b7TONGzfGvn37KIEihOiVRjVRXC4XHA6n2IbjHA4Hcrlcb8Hpm6lmsuXZmzdvMHr0aOzZs0c5rWHDhoiKioK5ubkRIyOmijGGFStWYNKkSZDJZCrv+fj4oF69erhx4wbq1q2LXr16ITw83OCPjyGEGJaplt8a1UTFxcUZOg5SzigUCqxevRpTpkxRtnvicrkYO3Ys5s2bRwkUKRSHw8G4ceMQFBSE8ePHQyQSoVGjRmjUqBGaNGmC2rVrGztEQkg5QW2iSKm7e/cuPvvsM5X2KiEhIdi8eTM1GieEEKLGVMtvnR5jv3nzZjRt2hRubm5ISEgAACxZsgT79u3Ta3Dk45KdnY2vv/4a9erVUyZQFSpUwB9//IErV65QAkUIIaRM0TqJWrVqFSZPnowOHTogLS1N2QbKzs4OS5Ys0Xd85CNx/vx51KlTBwsXLoRMJkOXLl0wYcIEJCYmon///ipd0wkhhJCyQOsk6pdffsG6devw7bffqoy1EhISgrt37+o1OFL2iUQifP7552jWrBni4uLg5uaG/fv3Y9++fViyZIlJVcsSQggh2tCoYfn74uLiEBQUpDZdKBQiOztbL0GRj8O1a9fQrVs35ZhP+Y9qCQsLM25ghBBCiB5oXRPl5eWFW7duqU0/fPgwatWqpY+YSBknkUgwdepUNGjQQJlA+fr64t69e5RAEUII+WhoXRM1efJkjB07FiKRCIwxXLlyBX/++SfmzZuH9evXGyJGUoY8ffoU7du3x6NHjwDkdUefMmUK5s6dC4FAYOToCCGEEP3ROokaPnw4LCws8N133yEnJwf9+vWDm5sbli5dis8++8wQMQLIG8E6vydgvnnz5mHatGkGWyfRzs6dOzF06FBkZmYCAFxdXbFnzx40bNjQyJERQggh+qdVEiWTybB161a0bdsW/fv3R05ODrKyskptNODZs2djxIgRyr9tbGxKZb2kaCKRCFOmTMHKlSsBAIGBgQgMDMTy5cvpkS2EEEI+WlolUXw+H6NHj0Z0dDQAwNLSEpaWlgYJrCA2NjZwcXEptfWR4sXGxqJ169Z49uwZAGDatGmYPXs23bojhBDy0dO6YXmDBg1w8+ZNQ8RSrPnz58PR0RFBQUFYtGiR2nOzPiQWi5GRkaHyIvqzdetW1K5dW5lALVmyBPPmzaMEihBCSLmgdZuoMWPGYMqUKUhKSkJwcLDa7ZqAgAC9Bfe+8ePHo169enBwcMCFCxcwffp0vHz5Ej/99FOhn5k3bx4iIyMNEk95plAoMGXKFJXBVcePH49x48YZLyhCCCGklGn97DwuV73yisPhgDEGDoejHMFcE9OmTcOCBQuKnCc6OrrAx4Fs2LABo0aNQlZWFoRCYYGfFYvFEIvFyr8zMjLg7u5ucs/eKUsyMzPRtWtXnDp1CgAgEAiwefNm9OnTx8iREUII+ViZ6rPzdBpsU1+mTJmCiIiIIufx9vYucHrDhg0hk8kQHx+PGjVqFDiPUCgsNMEi2ouPj0doaKjy9p2TkxOOHz9usNpHQgghxJRpnUR5eHjobeVOTk5wcnLS6bO3bt0Cl8sttZ6B5d25c+fw6aefIjk5GUDeY34OHz4MR0dHI0dGCCGEGIfWSZQxXLx4EZcvX0bLli1hY2ODixcvYtKkSRgwYADs7e2NHd5H77fffsOoUaMglUoRGBiI4cOHY9SoUeDzy8ThQwghhBhEmSgFhUIhtm3bhlmzZkEsFsPLywuTJk3C5MmTjR3aR40xhmnTpmHhwoUAgJ49e2Ljxo009hMhhBACHRqWl2Wm2jDNFDHGMGLECPz6668A8npd3rx5s8COBYQQQoghmWr5XSZqokjpkslk6NOnD3bv3g0AcHR0xP79+ymBIoQQQt6jdamY37W9IGvWrClRMMT4xGIx2rdvr0ygXF1dcfv2bb12KCCEEEI+BlonUe3atcPUqVMhlUqV05KTk9G5c2d6GHAZl5WVhRYtWuD48eMA8h76fOvWLVSuXNnIkRFCCCGmR6eaqD179qB+/fp48OABDhw4AH9/f2RkZODWrVsGCJGUhrS0NLRq1QpXrlwBANSoUQPXr1+nISQIIYSQQmidRDVp0gS3bt2Cv78/6tWrh08//RSTJk1CVFQU3fIpo7KystChQwdcuXIFtra2CAsLw5UrV+Dg4GDs0AghhBCTpVNL4UePHuHatWuoUqUK+Hw+Hj58iJycHH3HRkqBSCRCx44dcfHiRdjb2+PMmTM4deqUSfV+IIQQQkyR1knU/Pnz0bhxY7Rp0wb37t3DlStXcPPmTQQEBODixYuGiJEYiFQqRVhYGM6cOQNzc3McOnSIHuFCCCGEaEjrJGrp0qXYu3cvfvnlF5ibm8Pf3x9XrlxB9+7dERYWZoAQiSEoFAq0bdsWly9fBgA0b94cDRs2NHJUhBBCSNmh9ThRd+/eRcWKFVWmCQQCLFq0CJ06ddJbYMRwGGPo27evcriKFi1a4ODBg0aOihBCCClbtK6J+jCBel9oaGiJgiGGxxjDuHHjsH37dgBArVq1cOzYMXoOHiGEEKIlnUrOa9euYfv27Xj27BkkEonKe/mDNBLTtGjRIqxcuRIA4OTkhLNnz8LMzMzIURFCCCFlj9Y1Udu2bUOTJk0QHR2NPXv2QCqV4v79+zh58iRsbW0NESPRk6NHjyoHRLWwsMCZM2fg6Oho5KgIIYSQsknrJGru3Ln4+eef8ffff8PMzAxLly5FTEwMevfujapVqxoiRqIHjx8/Rp8+fcAYg5+fH/bv34+aNWsaOyxCCCGkzNI6iXry5Ak6duwIADAzM0N2djY4HA4mTZqEtWvX6j1AUnIZGRno0qUL0tLS0KhRI9y4cQOtW7c2dliEEEJImaZ1EmVvb4/MzEwAQOXKlXHv3j0AeY8NoQE3TY9CoUCnTp0QHR0NV1dX7N69G0Kh0NhhEUIIIWWe1g3LW7RogWPHjqFOnTro1asXJkyYgJMnT+LYsWNo1aqVIWIkJTBp0iScPXsWANCxY0e4uroaOSJCCCHk48BhjDFtPpCamgqRSAQ3NzcoFAosXLgQFy5cQLVq1fDdd9/B3t7eULGWWEZGBmxtbZGenl4uHmuyefNmDBo0CADg7u6O+/fvw8bGxshREUIIIdox1fJb6ySqLDPVL8EQbt68ifr160Mul8PS0hIPHjygB0QTQggpk0y1/NbpAcTEtOXm5qJNmzaQy+Xgcrk4evQoJVCEEEKInmncJorH42k0n1wu1zkYoh8jRoxASkoKAGDlypVo2rSpkSMihBBCPj4aJ1GMMXh4eGDw4MEICgoyZEykBE6fPo0tW7YAAJo1a4ZRo0YZOSJCCCHk46RxEnXlyhX8+uuvWLp0Kby8vDB06FD079/fpBuSlzeZmZmIiIgAAPTr1w8///yzcQMihBBCPmIat4kKCQnBqlWr8PLlS0yePBl79uxBlSpV8Nlnn+HYsWOGjJFoaPLkyYiPj4eXlxdWr14NZ2dnY4dECCGEfLS0blhubm6OAQMG4MSJE7h37x7evHmDdu3aITU11RDxEQ39/fffWL9+PQBg48aNNJQBIYQQYmBaD7YJAElJSdi4cSM2btyInJwcTJ061aS6HJY3KSkp6N+/PwBAIBCgdu3aRo6IEEII+fhpnERJJBLs2bMHv/76K86ePYv27dtjyZIlaN++vcY994hhDB48WPkonoULF8LR0dHIERFCCCEfP42TKFdXV9jY2GDw4MFYuXKlsr1Ndna2ynxUI1W6/vzzTxw4cAAAULduXYwfP97IERFCCCHlg8YjlnO5/zWf4nA4au8zxsDhcEx6nChTHfFUVxkZGahSpQoyMzPB4/Fw//591KhRw9hhEUIIIXplquW3xjVRp06dMmQcRAf/+9//lLfxIiMjKYEihBBCShE9O6+MSkxMhLe3N2QyGby9vRETEwOBQGDssAghhBC9M9XyW6MhDj5s96Tv+Yn2vvvuO8hkMvj5+WHnzp2UQBFCCCGlTKMkytfXF/Pnz8fLly8LnYcxhmPHjqF9+/ZYtmyZ3gIk6m7cuIHNmzcDyBsTih7DQwghhJQ+jdpERUVF4ZtvvsGsWbNQt25dhISEwM3NDebm5nj37h0ePHiAixcvgs/nY/r06fS8NgNijKFfv37Kf+vXr2/skAghhJBySas2Uc+ePcOOHTtw9uxZJCQkIDc3FxUrVkRQUBDatm1r8mNGmeo9VW2sX78eI0aMAADcunULdevWNXJEhBBCiGGZavlNDcvLEJlMBicnJ6SlpcHLywtPnz41dkiEEEKIwZlq+a31s/OI8SxZsgRpaWnK/xNCCCHEeLSqiXrw4AGWL1+Oixcv4tWrVwAAFxcXNG7cGOPGjYOfn5/BAtUHU81kNZGRkQFXV1fk5OTA3d0dCQkJBQ56SgghhHxsTLX81niwzUOHDqFbt26oV68eunbtikqVKgEAXr9+jWPHjqFevXrYt28f2rZta7Bgy7M5c+YgJycHADB37lxKoAghhBAj07gmqm7duujatStmz55d4PuzZs3C7t27cefOHb0GqE+mmskWJyUlBa6urpBKpXB2dsbz58/B52uc/xJCCCFlmqmW3xq3iXr06BH69+9f6Pt9+/ZFbGysXoIiqtavXw+pVAoOh4NvvvmGEihCCCHEBGicRHl6euLAgQOFvn/gwAF4eHjoJSjyH5lMhhUrVgAAfv75ZwwfPtzIERFCCCEE0KJN1OzZs9GvXz9ERUWhdevWKm2iTpw4gcOHD2Pr1q0GC7S82rdvHxITE1GxYkWMGjUK5ubmxg6JEEIIIdAiierVqxcqV66MZcuW4ccff1TrnRcVFYXGjRsbLNDyat68eQBACRQhhBBiYrRqXNOkSRM0adLEULGQD9y+fRvXr18HAGoHRQghhJgYGmzThH3//ffK//fr18+IkRBCCCHkQ3pLoqKjo+Ht7a2vxZV7ycnJ2Lt3LwCgRYsWqF69unEDIoQQQogKvSVREokECQkJ+lpcuffTTz9BLpcDyBtckxBCCCGmReOGNpMnTy7y/bdv35Y4mOIcOHAAs2fPxp07d2Bubo7Q0FBlbc3HRCaTYfny5QAAX19fNG3a1MgREUIIIeRDGidRS5cuRWBgYKEjhWZlZektqILs2rULI0aMwNy5cxEeHg6ZTIZ79+4ZdJ3G8ueffyIzMxMA1UIRQgghpkrjx77UqFEDM2bMwIABAwp8/9atWwgODlbegtInmUwGT09PREZGYtiwYTovx1SHjf9QUFAQbt26BQcHB7x9+xZcLrX/J4QQUn6ZavmtcekcEhKi7G5fEA6HAw3zMa3duHEDz58/B5fLRVBQEFxdXdG+fftia6LEYjEyMjJUXqbu5s2buHXrFng8Hvbu3UsJFCGEEGKiNC6hf/zxR0ycOLHQ9+vWrQuFQqGPmNQ8ffoUQN5Djr/77jv8888/sLe3R1hYGFJTUwv93Lx582Bra6t8ubu7GyQ+ffrll18A5A1u2rx5cyNHQwghhJDCaJxEubi46P3ZeNOmTQOHwynyFRMTo0zOvv32W/To0QPBwcH47bffwOFwsGPHjkKXP336dKSnpytfiYmJeo1f35KTk7F582YAwPjx440cDSGEEEKKYtRhsKdMmYKIiIgi5/H29sbLly8BAH5+fsrpQqEQ3t7eePbsWaGfFQqFEAqFeom1NPzyyy+QyWQwNzdH/fr1jR0OIYQQQoqgdRJlb28PDoejNp3D4cDc3By+vr6IiIjAkCFDil2Wk5MTnJycip0vODgYQqEQDx8+RLNmzQAAUqkU8fHxeq8dM6Zdu3YBADw9PekxL4QQQoiJ07qk/t///ocffvgB7du3R4MGDQAAV65cweHDhzF27FjExcXh888/h0wmw4gRI/QSZIUKFTB69GjMnDkT7u7u8PDwwKJFiwDktR36GIhEIkRHRwMAunTpYuRoCCGEEFIcrZOoc+fOYc6cORg9erTK9DVr1uDo0aPYtWsXAgICsGzZMr0lUQCwaNEi8Pl8DBw4ELm5uWjYsCFOnjwJe3t7va3DmA4fPqxs+zV06FAjR0MIIYSQ4mg8TlQ+a2tr3Lp1C76+virTHz9+jMDAQGRlZeHJkycICAhAdna2XoMtKVMdZwIAOnfujH/++QfW1tbIyMgo8JYpIYQQUh6Zavmt9SBEDg4O+Pvvv9Wm//3333BwcAAAZGdnw8bGpuTRlROMMURFRQEAGjduTAkUIYQQUgZofTtvxowZ+Pzzz3Hq1Cllm6irV6/i4MGDWL16NQDg2LFjCA0N1W+kH7Fbt24pH5vTr18/I0dDCCGEEE1onUSNGDECfn5+WL58OXbv3g0g75Ewp0+fRpMmTQDkDV1ANLd//34AQM2aNdGhQwcjR0MIIYQQTejUj75p06Zo2rSpvmMpt/Jvj3711VdwdnY2cjSEEEII0YROSZRcLsfevXuVXfJr166NLl26gMfj6TW48iApKQnXr18Hh8NBx44djR0OIYQQQjSkdRL1+PFjdOjQAc+fP0eNGjUA5D2jzt3dHQcOHICPj4/eg/yY5ddCVatWzaR6HBBCCCGkaFr3zhs/fjx8fHyQmJiIGzdu4MaNG3j27Bm8vLzoeW86+OuvvwAAT548gUwmM3I0hBBCCNGU1jVRp0+fxqVLl5TDGQCAo6Mj5s+fT+2ktJSVlYXz588DAIKCgmBtbW3kiAghhBCiKa1rooRCITIzM9WmZ2VlwczMTC9BlRfHjh1T1j517tzZyNEQQgghRBtaJ1GdOnXCyJEjcfnyZTDGwBjDpUuXMHr0aHrmm5b27dun/H+bNm2MGAkhhBBCtKV1ErVs2TL4+PigcePGMDc3h7m5OZo2bQpfX18sXbrUEDF+lORyuTKJsrS0RP369Y0cESGEEEK0oXWbKDs7O+zbtw+xsbGIiYkBANSqVUvtWXqkaFeuXEFaWhoAIDw8HHy+TqNNEEIIIcRIdC65q1WrhmrVqukzlnIlf5RyAGjbtq0RIyGEEEKILjRKoiZPnqzxAn/66SedgylP8pOoZcuW0fPyCCGEkDJIoyTq5s2bGi2Mw+GUKJjy4smTJ3jw4AH4fD4GDBgAe3t7Y4dECCGEEC1plESdOnXK0HGUKydPngSQ9wxCSqAIIYSQsknr3nmk5O7cuQMASExMVDbOJ4QQQkjZQkmUEdy4cQMA8PTpU9jZ2Rk3GEIIIYTohJKoUsYYw+3btwEAvr6+cHFxMXJEhBBCCNEFJVGl7MWLF8jOzgYAtG/f3sjREEIIIURXlESVsrt37yr/T+NDEUIIIWUXJVGlLL+nI4fDQWhoqJGjIYQQQoiuKIkqZadPnwYAVK1aFdbW1kaOhhBCCCG6oiSqlL1+/RoA0KJFCyNHQgghhJCSoCSqFEmlUjx//hwAMHv2bCNHQwghhJCSoCSqFD169AhSqRQ2Njbw8PAwdjiEEEIIKQFKokrR4cOHAQD+/v70nEFCCCGkjKMkqhQtX74cAGBubm7kSAghhBBSUpRElRKRSIRnz54BAJo0aWLkaAghhBBSUpRElZLz589DoVAAANq0aWPkaAghhBBSUpRElZJLly4p/x8QEGDESAghhBCiD5RElZLExEQAgKWlJezt7Y0cDSGEEEJKipKoUvLixQsAgKOjo5EjIYQQQog+UBJVSl6+fAkAqFixopEjIYQQQog+UBJVSkQiEQAgLCzMuIEQQgghRC8oiSoFjDFlm6iIiAjjBkMIIYQQvaAkqhQkJiYiPT0dfD4fNWvWNHY4hBBCCNEDSqJKwd27dwEAXl5eyMzMNHI0hBBCCNEHSqJKQX4SFRsbi2XLlhk5GkIIIYToAyVRpeDOnTvK//P5fCNGQgghhBB9oSSqFOTXRAGURBFCCCEfC0qiDEwikSAmJkb5NyVRhBBCyMeBkigDi4mJgUwmg0AgAEBJFCGEEPKxoCTKwPJv5dna2gKgJIoQQgj5WFASZWD5jcorVKgAgJIoQggh5GNBJbqB5ddENWvWDF26dEFAQICRIyKEEEKIPlASZWD5SdTIkSPRtGlTI0dDSotMJoNEIjF2GIQQYtLMzMzK9B2aMhF5VFQUWrZsWeB7V65cQf369Us5Is28e/cOSUlJAAB/f38jR0NKA2MMzxISkJySYuxQCCGkTKjo6IiqHh7gcDjGDkVrZSKJatKkCV6+fKkybcaMGThx4gRCQkKMFFXx8muhPDw8kJubi8zMTDg6OsLCwsLIkRFDeZaQgOTkZFSOXgvrlDvgMpmxQyKEEJOk4PCR5RiA57VGQiIWwadadXC5ZaupdplIoszMzODi4qL8WyqVYt++ffjiiy9MOnPNT6Lq1KmDPn364MyZM9i+fTt69epl5MiIIchkMiSnpKBy9Fq4PPnL2OEQQojJs06LBgA89xuFHavno16zNqgWYJp3lwpStlK+f+3fvx8pKSkYMmSIsUMpUnx8PACgWrVqkMnyaiTK8r1fUrT8NlDWKXeKmZMQQki+/GtmytvXOPTXOjx7/MDIEWmuTCZRv/76K9q2bYsqVaoUOZ9YLEZGRobKqzRJpVIAgFAopCSqHKFbeIQQorn8a6ZzFU9kZ6Qh9u51I0ekOaMmUdOmTQOHwyny9f4jUwAgKSkJR44cwbBhw4pd/rx582Bra6t8ubu7G2pTCiSXywEAPB6PkihCCCGkCBwOB+aW1nj17ImxQ9GYUUv0KVOmICIiosh5vL29Vf7+7bff4OjoiC5duhS7/OnTp2Py5MnKvzMyMko1kaIkihBCCNEch8uFXF52avONWqI7OTnByclJ4/kZY/jtt98waNAg5bPoiiIUCiEUCksSYokUlERpEjchpmjjLQmG7BMhboI1PO3KZEsAFZzIDMwMNcOsMHOV6QvPi7HhphQPxlqBa8IdV/QtJUeBqkuysKOXBTpUo+sUIZooU9UiJ0+eRFxcHIYPH27sUDRCNVGEaC4+TQGvpVmFvj+vlRDTmhn2R1GGmGHBeQkWtxGqJFCcyLz2lIvbCDGliWoM+cnl1RFWCHHjKaeniRi+OibCnhgZcqQMDSrz8OMn5qjnykNJJaYrsOGmFAdipYhNVYDH4cDfmYvvWgjR2rvoa8yI/blYf1OKjtX4+KefpXK6oyUXw4PMMOOUmJIoQjRUpkr0X3/9FU2aNEHNmjWNHYpG3k+i+vTpg1evXsHV1dXIURFi2vr689GhmvqlKcil5MlHcTbclECmYOhbp+AkYtEFCT6vbwZLQdE1VArG0HFrDm6/kmNqEyEqWnKw8poEYRuzcX2kFao5lmxb9j2UYcF5MbrV5GNwXTPIFAy/35GizeYcbOhijiFBZgV+7toLOTbelsK8kCv/6BABll2R4GScDOFeZap4IMQoytRZsnXrVmOHoJX3k6jZs2cbORpCyoZ6rjwMCCg4CSgMYwwiGWBRTHJTnN9uSdGlhgDmfPXlBLpwceuVAquvSTC5cdE1YjsfyHAhUY4dvSzQ0y8vIetdm4/qy7MwM0qMrT0si/x8cVp68vBskjUqWv53W3V0iBkC12Tjf1HiApMoxhjGHxJhUIAAJ+IKbnNSy4kHf2cuNt6SUhJFiAbKfsMGE/Z+EkXIx2jlVQlqr8yCcE4G3H7MxNgDuUgTMbX5VlyRwHtpJix+yECDdVk4myBD2MZshG3M1mm9nksy0WlrDo48liFkbRYsfsjEmut543SliRgmHhbB/edMCOdkwHdZJhacE0PB1ON6X9w7Be68VqC1V8Hna1N3HsK9eFh4XoJcadHL2vlAikpWHHSv9V8i4mTFRW8/AfY9lEEsK/rzxantzFNJoABAyOeggy8fSRkMmWL15W++I8W9N3L80KroBLCNNx9/P5KCFbO/CCGURBnU+0lUSkoK3r17p5xGSFk3K0qEsQdFcLPh4MdPzNGjFh9rrkvxyeZsSOX/FcCrrkow7pAIVSpwsbC1OZpX5aPbX7lIylAUuNwcKZCco1B7yRSqhfrDFAX67spBG28+lrYzR6ALDzlShtCN2fjjjhSDAgRY1s4cTavyMf2EGJOPiIvcnguJeedmUW2WZoUK8TqbYdW1oh8uffOVAvVceWoN0xtU5iFHCjxK+W/b3+WyArf3w1dOMYkbALzKVsBSAFh+cDcyU8zw9XExvmkuhIt10Zf9YFce0kTA/bcFfz+EkP9Qfa0BvZ9EVa9eHampqYiOji4zbbqIfjDGkCM1dhT/sRSgxI9LeputwLxzEnziw8Oh/pbKZKFmRR7GHRLhjztSDAkyg0TOMOOUGPXduDg52BJ8bt58AZW4iNgnQpUK6sueGSXGzCj1hOfiMEs0qvLfJetxqgKH+1uire9/0+acEeNJqgI3R/3X7mhUCOBmzcGiCxJMaWwGd9uCk4iY5Lzz1cu+8CSjuQcfLT15eW2jQswKvX34MlOBFlXVkzFXm7z5X2Qy1KmUNy1oTRYS0otPkArqSfi+x6kK7I6WoZefADyualyzT4thwQcmNSr+Nqm3fd5nH7xVwN+ZatEJKQolUQZEvfMIkFezYj0v09hhKGVNt4GVdk2O1Bx/KodEDkxsaKZS2zIiWIBvTopwIFaGIUFmuPZCjpRchnmthMoECgD6Bwgw6YiowGWPrCdAr9rqDbv9nFQLdC87jkoCBQA7HkjR3IMHewsOknP+q0lp7c3H/PMSnEmQo39AwUlSSi4DnwtYmxWdYM4KEyJ0Yw5WX5NgUiFto3JlgLCAUz2/rVXue7fztnS3QK4Gw+J4F5Hc5UgZeu3IgQUfmN9aNaZHKXIsvSzBnz0sICygrdeH7C3y5knOodt5hBSHSnQDoiSKfKwS0vMSlBoVVRMbMx4H3vZc5fsJaXkFsa+DagLA53IKHWuqmiO32G76QME1RrEpCtx5DTgtKniohDfZJb9F1eLf2qiFFyQYHVJwNmrBB8QFJEaif5Mni/eSmaZVS3ZNkCsYPtuZiwdvFTjU3xJuNqr7ZcJhEZq489DDT7NhC/KbQpWfEbII0R2V6AZESRQB8m6fZU23MXYYSh+2lymrLAqoVVEwoI03D181LbiGqLpj4bU5jhYcyBR57YdshEWnEDNDhQjblIM11yWwM1ef19WGi5dZ6jU5LzPzprnZ/PeZt9kKyDWo9LE24xRYSzbibxH+eSTDlu4Waj3qTsbJcPixHLt7WyA+7b8EUqbIqw2LT1PAwYKDCu9t77t/OwZUtKQ0ipDiUIluQO8nUfkPI6YkqvzhcDglvn1majz+bVf0MFmucptJImeIe6dQ1iR52OUVxI9TFWjp9d/nZYq8Ajygkn7b3Pg4cJElgUY1WR+q+W+tWpwGcYV68hHmycOC8xL8r4V6whbowsXZBDkUjKnc7rz8XA5LgWoyV39dts5toqYeFeG3W1IsaSsscGyrZ//WCHbfnqv23vNMObyWZuHntkJMbPTfNsS9y4ullhP1OyKkOFSiG1B+EsXhcJTdhSmJIh+D1t48mPGAZVckaOfLVzZU//WGFOlioOO/g2WGuPHgaMHBuhsSDAkSKNtFbbkjxbuCm0SVSG8/AWadFuPIY5lae6k0EYO1GVTaZr2vsXte4nTthVyj5G7Wv7VRa2+o99TrWUuAnQ9k2B0tU44TlZyjwI4HUnSuzldpm6Rrm6hF58VYfFGCb5qZYUKjgmvewr342NPHQm36yL9F8LDj4NvmQtT5oPH49Zdy2AqB2pREEVIsKtEN6P0kKh8lUeRj4GTFxfRmZog8LUG7LTnoUl2AhykKrLwqQX03LgYE5CUOZjwOZoUJ8cUhEcI35aB3bQHi0xTYeEsKH3sOCuokeOOlHH/cUU9MfOy5aOxe9PkztakZ9j+SotOfOYioK0CwGw/ZEoa7bxTY+UCK+InWhd6m8rbnwt+Zi+NPZRhayIjf7wv15CPUg4fTCerDlvT046PRZR6G7Mtrq1TRkoOVVyWQK4DIMNWER5c2UXuipfjquBjVHLio5cRV219tvPmoZM1FVdu814cmHhahkhUX3Wqq114deypD5xqCEvfgJKQ8oBLdgN6/nTdo0CDI5XKjPhCZEH2aFWYOJ0sull+VYNIRERwsOBgZLMDcVuYQ8P4rgMc1MANjwI8XxfjyqAh1XbjY39cC4w+JCnz8yJ/3ZPjznnrVzOC6gmKTKEsBB6cjrDD3rBg7Hsjw+x0pKgg5qO7IRWSYELbFtHUaGijA/6LEyJUyjUY/nxUmRMtNOWrTeVwODvazxNRjIiy7LEGujKG+Gw8bu1mpNcbXxe3XedeW2FQFBu5Rr9I7NdgSlYoZD6ogMcly3HujwJK2hQ+lQAj5D4eVo2FpMzIyYGtri/T0dFSoUMAANXrWvHlznDt3Djt27EDPnj0Nvj5iXDk5OYiOjkatM6NgmR5r7HBMmoIxOC3KQveafKzron67yVjSRQzey7KwsLUQw+p9ZA3ZNDDxsAhnEmS4PtKKaqJIqcmxrYboFmsQf/ssEmPvwcHJBYMmz1GZp7TLb03RTW8Dose+EJLXrf/D32q/35YiNZchzNO0KsNtzTn4qokZFl2QFPuYmI9NSo4C629IMCdcSAkUIRoyrSvYR+b9NlE5OTng8/kwMyt/v25J+XYpSY5JR0To5SeAowUHN17K8etNKfyduehV2/QuQV83E+LrZuXvtrujJRdZ35jOL3xCygLTu4J9RPKTqIyMDFhZWYHP5yuHOiCkvPC048K9AhfLLkuQmsvgYMHBoLoCzG8thBmPajwIIWUXJVEG9OHDhqlnHimPPO242N/X0thhEEKI3lGbKAOiJIoQQgj5eFESZUD5SRQNtEkIIYR8fCiJMiCqiSKEEEI+XpREGRDVRBFCCCEfL0qiDIhqogghhJCPF5XqBpSfRFWoUAE9evSAvb29kSMihBBCiL5QEmVA+UmUh4cHdu7caeRoCCGEEKJPdDvPgOixL+RjsfGWBJzIDMSnKYwdit5wIjMwK0r94b0Lz4tRc3lWuXzsi9XcDByMLf0BgePTFOBEZmDjLYnB15UlYXBelIktd9S3c8yBXLTZnG3wGErL4ccyWM/NwNvs0j9vV1+ToOrPmRDLPu7ziGqiDCg/ieJyKVclpDjxaQp4Lc0q9P15rYSYZuDHsWSIGRacl2BxGyG47z0/jhOZAQBY3EaIKU1UY9h4S4Ih+0S4OsIKIW7//WBKEzF8dUyEPTEy5EgZGlTm4cdPzFHPVT8/qlZdleBkvAyXk+RIzGAYXFeAjd0Kf5jz8acyzD0rxvWXcigYUN2Ri6+aCNHHXwAg77Evw4PMMOOUGB2qCfQSoylaekkCGyHwmb9q8Rf3ToH1N6Q4MkB1YFh97+f3TTkiwqHHMjwYa42rz+XYdFuCU/FyxKcp4GjBQaMqPMwJF6K6Y+HHjFTOUHd1NqKTFVjURogv3zs+2/ny4evAxbxzEvzU1rzQZTDGMGivCH/ckaJBZS5ODbaCpaD4pwmkiRiq/5KFtzkMO3pZoKfff9sYESjArCgx1lyXYHzDj/cxSlS6G1B+EnXp0iVwOBzUr1/fyBERYvr6+vOx+VNztVfn6ob/zbfhpgQyBUPfOgUnEYsuSJAjLf6XtYIxdNyag613pRhX3wwLW5vjTTZD2MZsxKbIi/28JhacF+NknBy1nbngF3Ml/+2mBJ9szoGAB8wNN8eiNuZoUZWPxAzVGorRIQLceKnAyTiZXmI0NVI5w9LLEgwPMgOPq5okLL0sgZc9Fy29VI8zQ+znfAdiZehYja9cz65oGVp58bG0nTlGBpvhTIIc9dZk496bwo+ZX65I8Cy98JqmUcFmWHNdgkxx4cftNyfE+OOOFB2q8XHthQKf7cyFXFH8cf6/U+JCzwdzPgeD6wrw00WJ2gPIPyZUE2VAHw5xQAgpXj1XHgYEaPegbsYYRDLAQoNfz0X57ZYUXWoIYM5XX06gCxe3Ximw+poEkxsX/ct65wMZLiTKVX6d967NR/XlWZgZJcbWHiV/DM7pCCtUteWAw+HAem5GofPFpykw9qAIXzQww9L2hddGAEAtJx78nbnYeEuKcC/ti4eIvbmIT1MgKsJK68+Whn8eyfA2h6F3bdUkWSpn2HJXitHB6smzIfYzADx9p8DDFAVW//vjYHJjM2ztwVN5nmSf2nzUWZWN+eck+KO7eu3Xm2wFZp8W4+umQvwvSlzgenr48fHFIWDHAymGBqmfV6uvSTD/vARfNzXD/Nbm+P12Xs3quIMirOpUeI3bvTdyrLomwf9aFL7u3rUFWHghr3ZNl+OpLKCaKAOicaLIx27lVQlqr8yCcE4G3H7MxNgDuUgTqf9oWHFFAu+lmbD4IQMN1mXhbIIMYRuzEbZRt/Ynnksy0WlrDo48liFkbRYsfsjEmut57WnSRAwTD4vg/nMmhHMy4LssEwvOiYtt4xT3ToE7rxVo7VXwrZOm7jyEe/Gw8LwEucXURu18IEUlKw661/rvnHey4qK3nwD7Hsr00k7Ew44LDqf4pHH1NQnkDJjdMi/xy5KwIn/YtfHm4+9HUpP48XcyTobmv2XDam4G7OZnoOu2HES/Va+ViYrPOw7M52TAZ1km1lyTYFaUSHkbNt/ehzJ42nHg46Ba9J17JkdyDkNrb/VrtKH284FHMtgKgWZV8463Ju58tQdyV3PkobYzF9HJBddETTsuRo2KXAwIKPz2q7MVFwGVuNj3UL128e+HUow7KML0ZnkJFAAMqmuGjV3NsfaGFPPOFpwcAcCEwyJ8WpOP5h6F32oMduPBwYKDfTEfZ80mQDVRBpWfRCkUeVWtlESRj8msKBEiT0vQ2puHz0PM8TBZjlXXpLj6Qo7zQ60g+LdAWHVVgnGHRGhelYdJjQSIT1Og21+5sDcHqlRQ/x2XIwWSc9RvT9iZc8B/7xbMwxQF+u7KwahgM4yox0WNilzkSBlCN2bjeQbDqGABqtpycSFJjuknxHiZxbCkXeE1BBcS887XotoszQoVosXGHKwqpjbq5isF6rnyVNpVAUCDyjysvSHFoxQF6lTKW8+7XAa5BgmLpYCjUTuVDx1/KkPNilwcjJVh6jERnmcy2JsDY+ubIbKlUC3GYFcefr4E3H+rgL+z8TrFHH8qQ/stOfC252JWqBC5srxbV003ZOPGKGt42uUdOzdfytHujxy42nAQGSbMS2TOiOFkqb6vLiTKC/x+LyTKwQEQVIL2atru54OPpWjjw1c5pj/EGMPrLIbazurnyZXncmy6LcW5IZYoLscLduVh7wdJ1NXncny2KxfTmplhTrjqeTGwrhk4HGDIPhHcbTlqNcM77ktxIVGO6LHWxXY2qefKxflESqKIDqgmirwvW1J4QcnjQuUWUlHzcjmqt620mVdf3mYrMO+cBJ/48HCov6WygKhZkYdxh/IaqA4JMoNEzjDjlBj13bg4OdhSWWAEVOIiYp8IVSqoL3tmlBgzC7g9cHGYJRpV+e8cepyqwOH+lmjr+9+0OWfEeJKqwM1RVqj2b2PcUSGAmzUHiy5IMKWxGdxtC66Aj/n3176XfeEV9M09+GjpycOiCxJ8HmJW6L59malAi6rqBbKrTd78LzIZ6lTKmxa0JgsJ6cUnUTNDzTArrPjbRB+KTVWAxwGG7MvFV03NULcSD7tjpJhzVgKZApjXWnWZ3vZ5MT4wchI19ZgIDhYcXBxmBQeLvJi61eQjaE02ZkaJsenfxt0zo8TgcYHzQ63gZpP33fWuLUCtFaqdFGQKhiepCnStoX4djklRwMGCgwpC3c8VbfZzjpQhKl6OVR2L/j633JXieSbD7JaqNU2MMXxxKBd9avPR2J1fbCLjbc9Fcg7Dm2wFnK3y9lH9yjxkf1PACfivAQFmBd5Wz5UyfHlMhEmNzOBpxy1+3XZcbH5W+j0+SwuV6gb0YU2UQPDx9nghxbOel1noex2q8XGg33/tZJwXZyKnkOtOqAdPpc2J59IsJOcUXAiHuHFxdYS1bgEX4fhTOSRyYGJDM5Vf2COCBfjmpAgHYmUYEmSGay/kSMllmNdKqPKLu3+AAJOOqA8vAAAj6wnQq7b6ueLnpFqge9lxVBIoIK/dR3MPHuwtOCq1Wa29+Zh/XoIzCXL0Dyg4SUrJZeBzAWuzogvSWWFChG7MweprEkwqpDYqVwYIC7i65ifKue/dztvS3QK5GvxQ9y4iuStKlgRQMGB+KyG+/rd3Yw8/AVJzs7H0sgTfNBfC5r3kwf7fhKWwYyqfgjGk5qrOI5YzSBXqNYm2Qo6yZlITLzMVuPVKga+amCkTKAAIqMRDG28eDsbm7TC5guH4Uxk+rcVXJlAA4OvARXtfPv5+9N+OTc1lYADszdXjSMlhyu3WlTb7+WScDGIZ0N638CI4JlmOsQdFaFyFh8F1Vc+HjbekuPtagZ29NGtb9/536lzC5mrzz4khlQPfNNesx529BQe5srzEUZeaVFNHSZQBUU0U+Vgl/NsbqEZF1cTGjMeBtz1X+X5CWt6x7/tBGxQ+l6O8HfOhao7cAtumfKigGqPYFAXuvAacFhU8VMIbPYyX0+Lf2qiFFyQYHVJwA3gLPiAuIDES/Zs8WbxX69i0qmGvCxZ8IFsKtR6Hff0FOPxYjpuv5Gjh8V8M+XcWiyvunqWzQoek+HD/nxpsiTBPzbfzv+NL/TuuVZGHI0/kyJYwZIgZcmWAbwHHwofHXL7CUsOStgHTZj8feCRDiBsXlawLjvFVlgIdt+bAVsjBzt4WKj0JM8QM00+IMbVJ4bWqH9L0Oy1OfJoCiy5IsKKDebE/NvS9blNFpbqB5Nc+AYCLiwvatWuH4OBgI0ZEjC1ruk2h7/E+uBa++bLweT9sQhE/ofCapiKaW5R5FgX0oFMwoI03D181LfhXcnXHwgsdRwsOZAogU8xUamYKMjNUiLBNOVhzXQK7Amo2XG24eJmlXii/zMyb5mbz32feZisg16D8tjbjaFxwvc/NhovYVAUqWal+Nv+2zrsPapPe/dsxoGIBbYre52LNwbGBqjUhiy6I8SqL4cdPVG9T1a1k/AGHHSw44EB9ewHA0ZKDd89Ltnxt9vPBxzIMCSz4zkS6iKH9lhykiYCzQyxVatgAYPEFMSRyhj7+AuWttKR/h1B4l8sQn6aAmw1HpZG6pt9pcf53SozKFbgI8/zvFuKrf4/zt9l5665qy1GpnX4nYrAUGKZZgSmgJMpA3n/4cOvWrdG7d28jRkNMgZUWBaCh5tUXj39/AT9MlqvcZpLIGeLeKZQ1SR52ebE9TlWgpdd/n5cp8i64AXouXH0cuMiSQKOarA/V/LdWLU6DuEI9+Qjz5GHB+bwu3h8KdOHibIIcCsZUCpTLz+WwFKgmc/XXZRu0TVSwW17h/jyTKds7AcCLzLxC0OmDQj/uXV4stZyKruUw53PU9vMfd6QQyxQ67f/3/Xd8qdccxqTIUdGSAyszDsz5gDkfePxOfb7HqarT+Ny8XnlxBbThqenIxZY7DOkiBtsCkmJNaLqf772R41k6Q8cCBjQVyRg6/5mDRykKHB9oqXYLG8irAXwnAmqvVO/ZOvecBHPPSXBzlBUCXf77bNw7BSpacuBkVbIO+c/SFXicqoD3MvUayDEH827Pv/vaBnbvHaZxaQrUKqBG8WNBSZSBvJ9E0WNfyMemtTcPZjxg2RUJ2vnylV3Af70hRboYygEEQ9x4cLTgYN0NCYYECZTtorbckeJdwU2iSqS3nwCzTotx5LFMrb1UmojB2gyF9oZq7J53nl57IdcouZv1b23U2hvqjyrpWUuAnQ9k2B0tU44TlZyjwI4HUnSuzofwvVo0Q7eJ6lNbgG33ZPj1hgQ/tMor3RSM4bdbUjhYcBD8QY+06y/lsBUCtYtJogzJ1YaLQBcuNt2WYnpzobK2794bOY4+kSu79PO4eYnc3hgZXmQqlLU2j1MVOPRYfac2rsJDVHwB0915YMjbdl3HM9J0Px+MlaGSFQchbqr7V65g6LMzFxeT5Nj3mQUauxccx/iGZuhWU/W9N9kMo/4RISJQgK41+PD64Fb59ZdyNK5S8nJoTrhQra3cvTcKzDglxldNzNDYnQerD3LDGy8V6F/I4LUfA0qiDISSKPIxc7LiYnozM0SelqDdlhx0qS7AwxQFVl6VoL7bf+PWmPE4mBUmxBeHRAjflIPetfNuQWy8JYWPPafArtk3Xsrxxx31xMTHnltowZJvalMz7H8kRac/cxBRV4BgNx6yJQx33yiw84EU8ROtC72l4W3Phb8zF8efygoclPBDoZ58hHrwcDpBfQyfnn58NLrMw5B9uXjwNq8WYOVVCeQKIDJMteZK1zZRfz+U4vbrvFoOqQK481qOOWfyejV2qcFXJoJda/DRyouHeeckSM5hqOvCw94YKc49k2NNJ3OVhA4Ajj2VoXMNgUZjIxnSojbmaL8lB41/zcawIAFypXlDHNgK8xLYfLNChTj6RIamG7LxeYgZ5Apg+VUJ/J3zBkd9X9cafGy+I8WjFLnKo1SaVc1L9o8/laklUfrezwdiZWhfja+2f6ccFWP/Qxk6V+cjNZepnQP5PeXqufLUhmnIv7VW24mLbjVVE5Y32Xnjn42tr90AtgVpVsCxameel5TWr8xTW/f1F3Kk5rICe0R+LD7eLTOy95Oo1atXY9asWejXrx/Wrl1rxKgI0Z9ZYeZwsuRi+VUJJh3J644+MliAua3MVXpijWtgBsaAHy+K8eVREeq6cLG/rwXGHxLBvIAr0J/3ZPjznnptweC6gmKTKEsBB6cjrDD3rBg7Hsjw+x0pKgg5qO7IRWSYELbFtHUaGijA/6LEyJUyjdpwzAoTouWmHLXpPC4HB/tZYuoxEZZdliBXxlDfjYeN3azUGuPrale0DJtu/9eF8+YrBW6+yivcq1TgKAt3DoeDvZ9Z4ruTYvx1X4qNt6Wo4cjFH59aoP8HgzTGJMtx740CS4p4zlppae3Nx+H+lpgZJcb/Tokh4AGhHnwsaC1U6VQQ7JY3zMaXR0WYcUoM9woczA4TIjpZgZhk1USkcw0+KlpysP2+DN+1+O97MONx0L+OADseyDC3lWoc+tzP6SKGC4lyjCsgobn1Kq/M+PuRTKVXYT5tR/HPtztaBiEfaqO0l4YdD6SoastBeCED2H4MOMwUhqUtJRkZGbC1tUV6ejoqVCh8fAx9ePfuHRwcHAAAP/zwA7799ltERETgt99+M+h6ifHk5OQgOjoatc6MgmV6rLHDMWkKxuC0KAvda/Kxrkvhj5YobekiBu9lWVjYWohh9Ur+y72smXhYhDMJMlwfaWX0mqiS6rYtB/ffKhD7hWrHi+9Pi/HbLQliv7BW6fX29J0CNZdn4VB/S7QqYZuuwmy/L0X/3blInmqjc9srbQWtyUKYBx8/FzHQrCGIZQyeS7MwrakZJjQqejiEHNtqiG6xBvG3zyIx9h4cnFwwaPIclXlKs/zWxsfb2svI3q+Jyv8/DXFAyiORTP3xF7/fliI1l2nV7b002Jpz8FUTMyy6ICn2MTEfm5QcBdbfkGBOuLDMJVAfPoYnNkWOg7EyhBXwSJJJjc2QJQG2fVDb6W3PxbAgAeafL/xRJyVlZ87BsnbmpZZAHX4sQ2yKAtObl/4Pgt9uSSHgotBhQD4WpnUF+4hQEkVInktJckw6IkIvPwEcLTi48VKOX29K4e/MRa/apndOfN3sv8ESyxNHSy6yihjB2pR5L8tCRF2BcoyyVdekMOMBXzVVL8CtzTh4M7XgIUSKeuCuPnziwwd8DLoKFe18+Ub7TkeHmH30CRRASZTB5CdOXC6XkihSrnnaceFegYtllyVIzWVwsOBgUF0B5rcWqj1wlRBdtPPl4897UrzKYhDy83rhzW1lrnz0DyGGQqW6geQnTjweDzJZXrUxJVGkPPK042J/X80eT0GILn7rajrt6kj5Qm2iDOT9JEoqzevZQUkUIYQQ8vGgJMpA3k+ivLy80Lx5c3h7exs5KkIIIYToC1WNGMj7SdTYsWMxduxYI0dESouCQ6cVIYRoSnnNLIM9YqkmykDeT6JI+WBmltcTJcsxwMiREEJI2ZF/zZSKc40cifboJ7OBUBJV/vD5fFR0dMTzWiMBANYpd8BlGjwUjRBCyiEFh48sxwA8rzUCaa8SoJDnXS/LUp9dSqIM5P0k6osvvsD27dsxc+ZMjBkzxsiREUOq6uEBAHjuN8rIkRBCSNmQ9ioBr5/eAwDIpRIILayMHJHmKIkykPeTqLS0NLx58wYikQEeW09MCofDgYenJ55FX8fNCyfgWtUHXKqNJIQQdYxBKs5V1kDJZFKIcrPhUd3fyIFprswkUY8ePcLUqVNx/vx5SCQSBAQE4Pvvv0fLli2NHVqBaJyo8q12cBM8uX8NsbcuQiA0B4++e0IIKZRcJoNUIoJHNX/UDGpk7HA0Vmau7J06dUK1atVw8uRJWFhYYMmSJejUqROePHkCFxcXY4enhpKo8s3B2RVdIybi8b3rSHr6EFKJCCh7HU8IIcTwOICZ0ALu3jVRrU4IKjhUNHZEGisTpXpycjJiY2Px66+/IiAgrxX//PnzsXLlSty7d4+SKGKS7CtWQv2wDqgf1sHYoRBCCDGAMjHEgaOjI2rUqIHff/8d2dnZkMlkWLNmDZydnREcHFzo58RiMTIyMlRepYWSKEIIIeTjViZKdQ6Hg+PHj6Nbt26wsbEBl8uFs7MzDh8+DHt7+0I/N2/ePERGRpZipP+hJIoQQgj5uBm1JmratGngcDhFvmJiYsAYw9ixY+Hs7IyzZ8/iypUr6NatGzp37oyXL18Wuvzp06cjPT1d+UpMTCy1bXs/ifL19UW9evVQsWLZuc9LCCGEkKJxGDPeOOtv375FSkpKkfN4e3vj7Nmz+OSTT/Du3TtUqFBB+V61atUwbNgwTJs2TaP1paenw87ODomJiSrLMYQTJ06ge/fu8Pf3x/nz5w26LkIIIeRjlpGRAXd3d6SlpcHW1tbY4SgZ9f6Sk5MTnJycip0vJycHAMDlqlaccblcKBQKjdeXmZkJAHB3d9ciypK5d++eSX3hhBBCSFmVmZlpUmVqmWik07hxY9jb22Pw4MH43//+BwsLC6xbtw5xcXHo2LGjxstxc3NDYmIibGxswOHob2D5/Ay5NGq4TFV53wflffsB2gcA7YPyvv0A7QNDbT9jDJmZmXBzc9PbMvWhTCRRFStWxOHDh/Htt98iPDwcUqkUtWvXxr59+1C3bl2Nl8PlclGlShWDxVmhQoVyedK8r7zvg/K+/QDtA4D2QXnffoD2gSG235RqoPKViSQKAEJCQnDkyBFjh0EIIYQQAqCMjBNFCCGEEGJqKInSA6FQiJkzZ0IoFBo7FKMp7/ugvG8/QPsAoH1Q3rcfoH1Q3rbfqEMcEEIIIYSUVVQTRQghhBCiA0qiCCGEEEJ0QEkUIYQQQogOKIkihBBCCNEBJVF6sGLFCnh6esLc3BwNGzbElStXjB1SqTlz5gw6d+4MNzc3cDgc7N2719ghlap58+ahfv36sLGxgbOzM7p164aHDx8aO6xStWrVKgQEBCgH12vcuDEOHTpk7LCMZv78+eBwOJg4caKxQyk1s2bNUnt4fM2aNY0dVql6/vw5BgwYAEdHR1hYWKBOnTq4du2ascMqNZ6enmrHAIfDwdixY40dmkFRElVCf/31FyZPnoyZM2fixo0bqFu3Ltq2bYs3b94YO7RSkZ2djbp162LFihXGDsUoTp8+jbFjx+LSpUs4duwYpFIpPvnkE2RnZxs7tFJTpUoVzJ8/H9evX8e1a9cQHh6Orl274v79+8YOrdRdvXoVa9asQUBAgLFDKXW1a9fGy5cvla9z584ZO6RS8+7dOzRt2hQCgQCHDh3CgwcP8OOPP8Le3t7YoZWaq1evqnz/x44dAwD06tXLyJEZGCMl0qBBAzZ27Fjl33K5nLm5ubF58+YZMSrjAMD27Nlj7DCM6s2bNwwAO336tLFDMSp7e3u2fv16Y4dRqjIzM1m1atXYsWPHWGhoKJswYYKxQyo1M2fOZHXr1jV2GEbz9ddfs2bNmhk7DJMyYcIE5uPjwxQKhbFDMSiqiSoBiUSC69evo3Xr1sppXC4XrVu3xsWLF40YGTGW9PR0AICDg4ORIzEOuVyObdu2ITs7G40bNzZ2OKVq7Nix6Nixo8r1oDyJjY2Fm5sbvL290b9/fzx79szYIZWa/fv3IyQkBL169YKzszOCgoKwbt06Y4dlNBKJBH/88QeGDh0KDodj7HAMipKoEkhOToZcLkelSpVUpleqVAmvXr0yUlTEWBQKBSZOnIimTZvC39/f2OGUqrt378La2hpCoRCjR4/Gnj174OfnZ+ywSs22bdtw48YNzJs3z9ihGEXDhg2xceNGHD58GKtWrUJcXByaN2+OzMxMY4dWKp4+fYpV/2/v/kKaagMwgD9j61RTcWWmm7Etyq1EKnMoJEHhTSOC/kgjBGd/bsJhVLvpKuhi3jq7MAZxlEAkIossWmXb7soQBrMosSyDYCPCcAkrtvNdffINu/lOdt7pnh8c2A6Hl+fcPee87zmnvx+1tbUIh8M4f/48uru7MTg4KDqaEPfu3cPc3Bw6OztFR/nrVswHiIkKXVdXFyYnJ4tqLci/nE4n4vE4vn//jjt37sDr9SIWixVFkfr8+TMuXLiAp0+fYt26daLjCOF2uxd/79q1C83NzbDZbLh9+zbOnj0rMJk2crkcXC4XAoEAAKChoQGTk5O4ceMGvF6v4HTau3nzJtxuNywWi+gofx3vRP2BTZs2Qa/XI5lM5u1PJpOorq4WlIpE8Pl8GB0dRSQSwZYtW0TH0ZwkSdi+fTsaGxvR09OD3bt3IxgMio6liYmJCaRSKezduxcGgwEGgwGxWAx9fX0wGAzIZrOiI2rOZDLB4XBgenpadBRNmM3mJRcMO3fuLKopzX99+vQJz549w7lz50RH0QRL1B+QJAmNjY0YGxtb3JfL5TA2NlZ060GKlaIo8Pl8GBkZwfPnz7F161bRkQpCLpdDJpMRHUMTra2tSCQSiMfji5vL5UJ7ezvi8Tj0er3oiJpLp9N4//49zGaz6CiaaGlpWfJqk6mpKdhsNkGJxJFlGZs3b8bhw4dFR9EEp/P+0KVLl+D1euFyudDU1ITe3l78+PEDp0+fFh1NE+l0Ou9qc2ZmBvF4HBs3boTVahWYTBtdXV0YGhrC/fv3UVZWtrgWrry8HOvXrxecThtXrlyB2+2G1WrF/Pw8hoaGEI1GEQ6HRUfTRFlZ2ZI1cCUlJaioqCiatXF+vx9HjhyBzWbDly9fcPXqVej1epw6dUp0NE1cvHgR+/btQyAQwMmTJzE+Po5QKIRQKCQ6mqZyuRxkWYbX64XBUCT1QvTjgavB9evXFavVqkiSpDQ1NSkvXrwQHUkzkUhEAbBk83q9oqNp4nfnDkCRZVl0NM2cOXNGsdlsiiRJSmVlpdLa2qo8efJEdCyhiu0VBx6PRzGbzYokSUpNTY3i8XiU6elp0bE09eDBA6W+vl5Zu3atsmPHDiUUComOpLlwOKwAUN69eyc6imZ0iqIoYuobERER0crFNVFEREREKrBEEREREanAEkVERESkAksUERERkQosUUREREQqsEQRERERqcASRURERKQCSxQRFbRoNAqdToe5uTnRUYiI8vBlm0RUUA4cOIA9e/agt7cXAPDz5098+/YNVVVV0Ol0YsMREf1HkXzchohWKkmSUF1dLToGEdESnM4jooLR2dmJWCyGYDAInU4HnU6HgYGBvOm8gYEBmEwmjI6Owul0wmg0oq2tDQsLCxgcHITdbseGDRvQ3d2NbDa7OHYmk4Hf70dNTQ1KSkrQ3NyMaDQq5kSJaFXgnSgiKhjBYBBTU1Oor6/HtWvXAACvX79ectzCwgL6+vowPDyM+fl5HD9+HMeOHYPJZMKjR4/w4cMHnDhxAi0tLfB4PAAAn8+HN2/eYHh4GBaLBSMjIzh06BASiQRqa2s1PU8iWh1YooioYJSXl0OSJBiNxsUpvLdv3y457tevX+jv78e2bdsAAG1tbbh16xaSySRKS0tRV1eHgwcPIhKJwOPxYHZ2FrIsY3Z2FhaLBQDg9/vx+PFjyLKMQCCg3UkS0arBEkVEK47RaFwsUABQVVUFu92O0tLSvH2pVAoAkEgkkM1m4XA48sbJZDKoqKjQJjQRrTosUUS04qxZsybvv06n++2+XC4HAEin09Dr9ZiYmIBer8877r/Fi4jo/2CJIqKCIklS3oLw5dDQ0IBsNotUKoX9+/cv69hEVLz4dB4RFRS73Y6XL1/i48eP+Pr16+LdpD/hcDjQ3t6Ojo4O3L17FzMzMxgfH0dPTw8ePny4DKmJqBixRBFRQfH7/dDr9airq0NlZSVmZ2eXZVxZltHR0YHLly/D6XTi6NGjePXqFaxW67KMT0TFh28sJyIiIlKBd6KIiIiIVGCJIiIiIlKBJYqIiIhIBZYoIiIiIhVYooiIiIhUYIkiIiIiUoElioiIiEgFligiIiIiFViiiIiIiFRgiSIiIiJSgSWKiIiISAWWKCIiIiIV/gHyQeQxeU8ELgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import mpmath as mp\n",
    "import csv\n",
    "\n",
    "def file_reader(filename):\n",
    "    with open(filename) as file:\n",
    "        reader = csv.reader(file, delimiter=\" \")\n",
    "        data  = list(zip(*reader))\n",
    "        # data is a tuple of strings. Tuples are immutable, and we need to perform math on\n",
    "        #   the data, so here we convert tuple to lists of floats:\n",
    "        data0 = []\n",
    "        data1 = []\n",
    "        for i in range(len(data[0])):\n",
    "            data0.append(float(data[0][i]))\n",
    "            data1.append(float(data[1][i]))\n",
    "        return data0,data1\n",
    "\n",
    "first_col16,second_col16 = file_reader(os.path.join(outdir, 'out-lowresolution.txt'))\n",
    "first_col24,second_col24 = file_reader(os.path.join(outdir, 'out-medresolution.txt'))\n",
    "\n",
    "second_col16_rescaled4o = []\n",
    "second_col16_rescaled5o = []\n",
    "for i in range(len(second_col16)):\n",
    "    # data16 = data24*(16/24)**4\n",
    "    # -> log10(data24) = log10(data24) + 4*log10(16/24)\n",
    "    second_col16_rescaled4o.append(second_col16[i] + 4*mp.log10(16./24.))\n",
    "    second_col16_rescaled5o.append(second_col16[i] + 5*mp.log10(16./24.))\n",
    "\n",
    "# https://matplotlib.org/gallery/text_labels_and_annotations/legend.html#sphx-glr-gallery-text-labels-and-annotations-legend-py\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "plt.title(\"Demonstrating 4th-order Convergence: \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" Coordinates\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(\"log10(Max relative error)\")\n",
    "\n",
    "ax.plot(first_col24, second_col24, 'k-', label='logErel(N0=24)')\n",
    "ax.plot(first_col16, second_col16_rescaled4o, 'k--', label='logErel(N0=16) + log((16/24)^4)')\n",
    "ax.set_ylim([-8.05, -1.7]) # Manually set the y-axis range case, since the log10\n",
    "                          # relative error at t=0 could be -inf or about -16,\n",
    "                          # resulting in very different-looking plots\n",
    "                          # despite the data being the same to roundoff.\n",
    "if par.parval_from_str(\"reference_metric::CoordSystem\") == \"Cartesian\":\n",
    "    ax.set_ylim([-2.68, -1.62])\n",
    "if par.parval_from_str(\"reference_metric::CoordSystem\") == \"Cylindrical\":\n",
    "    ax.plot(first_col16, second_col16_rescaled5o, 'k.',  label='(Assuming 5th-order convergence)')\n",
    "legend = ax.legend(loc='lower right', shadow=True, fontsize='large')\n",
    "legend.get_frame().set_facecolor('C1')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: 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-ScalarWaveCurvilinear.pdf](Tutorial-Start_to_Finish-ScalarWaveCurvilinear.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": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-28T19:19:20.724049Z",
     "iopub.status.busy": "2021-09-28T19:19:20.723767Z",
     "iopub.status.idle": "2021-09-28T19:19:23.365936Z",
     "shell.execute_reply": "2021-09-28T19:19:23.365481Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish-ScalarWaveCurvilinear.tex, and compiled\n",
      "    LaTeX file to PDF file Tutorial-Start_to_Finish-\n",
      "    ScalarWaveCurvilinear.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-ScalarWaveCurvilinear\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}