{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8f402aaa",
   "metadata": {},
   "source": [
    "# Tutorial: `NRPyCritCol` - ADM initial data for scalar field collapse\n",
    "\n",
    "## Author: Leo Werneck"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44597ba6",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "0. [Step 0](#introduction): Introduction\n",
    "    1. [Step 0.a](#problem_statement): Problem statement\n",
    "    1. [Step 0.b](#boundary_value_problem): The boundary value problem\n",
    "1. [Step 1](#initialize_nrpy): Initialize core Python/NRPy+ modules\n",
    "1. [Step 2](#thomas_algorithm): Thomas' algorithm for tridiagonal systems\n",
    "1. [Step 3](#initial_scalar_field_profile): Initial scalar field profile\n",
    "1. [Step 4](#solving_the_hamiltonian_constraint): Solving the Hamiltonian constraint\n",
    "1. [Step 5](#code_validation): Code validation\n",
    "    1. [Step 5.a](#standalone_main_funcion): Generate the `main` function of the standalone C code\n",
    "    1. [Step 5.b](#c_code_generation): Generate all C codes and compile the project \n",
    "    1. [Step 5.c](#running_c_code): Running the C code\n",
    "    1. [Step 5.d](#running_python_code): Running the (trusted) Python code\n",
    "    1. [Step 5.e](#relative_error_plot): Plotting the relative error between the results of the two codes\n",
    "1. [Step 6](#output_to_pdf): Output this module as $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15a9ae39",
   "metadata": {},
   "source": [
    "<a id='introduction'></a>\n",
    "\n",
    "# Step 0: Introduction\n",
    "$$\\label{introduction}$$\n",
    "\n",
    "<a id='problem_statement'></a>\n",
    "\n",
    "## Step 0.a: Problem statement \\[Back to [top](#toc)\\]\n",
    "$$\\label{problem_statement}$$\n",
    "\n",
    "In order to obtain initial data for our simulations, we must solve the Hamiltonian and momentum constraint equations, which are given by, respectively,\n",
    "\n",
    "$$\n",
    "\\begin{alignat}{3}\n",
    "& \\mathcal{H} && = R + K^{2} - K_{ij}K^{ij} - 16\\pi\\rho && = 0\\;,\\\\\n",
    "& \\mathcal{M}^{i} && = D_{j}\\bigl(K^{ij} - \\gamma^{ij}K\\bigr) - 8\\pi S^{i} && = 0\\;,\\\\\n",
    "\\end{alignat}\n",
    "$$\n",
    "\n",
    "where $\\gamma_{ij}$ is the physical spatial metric, $K_{ij}$ is the extrinsic curvature, $K \\equiv \\gamma^{ij}K_{ij}$ is the mean curvature, $R$ is the three-dimensional Ricci scalar, and $D_{i}$ is the spatial covariant derivative compatible with $\\gamma_{ij}$.\n",
    "\n",
    "We now consider a spherically symmetric massless scalar field $\\varphi$ minimally coupled to gravity. Assuming a moment of time symmetry, for which\n",
    "\n",
    "$$\n",
    "\\partial_{t}\\varphi = K_{ij} = K = 0\\;,\n",
    "$$\n",
    "\n",
    "and the momentum constraints are automatically satisfied. Further assuming that the metric is initially conformally flat, i.e.,\n",
    "\n",
    "$$\n",
    "\\gamma_{ij} = \\psi^{4} \\bar\\gamma_{ij} = \\psi^{4} \\hat\\gamma_{ij}\\;,\n",
    "$$\n",
    "\n",
    "where $\\psi$ is the conformal factor, $\\bar\\gamma_{ij}$ is the conformal metric, and $\\hat\\gamma_{ij}$ is the flat spatial metric, the Hamiltonian constraint can be written as (see Eq. 60  of [Werneck et al. 2021](https://arxiv.org/pdf/2106.06553.pdf))\n",
    "\n",
    "$$\n",
    "\\hat\\gamma^{ij}\\hat D_{i}\\hat D_{j}\\psi = -2\\pi\\psi^{5}\\rho = \\pi\\bigl(\\hat\\gamma^{ij}\\partial_{i}\\partial_{j}\\varphi\\bigr)\\psi\\;.\n",
    "$$\n",
    "\n",
    "For a given initial profile of the scalar field $\\phi(t=0,r)$, our task is to solve this last equation for the conformal factor."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdfc03cb",
   "metadata": {},
   "source": [
    "<a id='boundary_value_problem'></a>\n",
    "\n",
    "## Step 0.b: The boundary value problem \\[Back to [top](#toc)\\]\n",
    "<a id='boundary_value_problem'></a>\n",
    "\n",
    "In spherical symmetry we have\n",
    "\n",
    "$$\n",
    "\\partial^{2}_{r}\\psi + \\frac{2}{r}\\partial_{r}\\psi + \\pi\\Phi^{2}\\psi = 0\\;,\n",
    "$$\n",
    "\n",
    "where $\\Phi \\equiv \\partial_{r}\\varphi$. Discretizing this last equation using [second-order accurate, centered finite differences](https://en.wikipedia.org/wiki/Finite_difference_coefficient) and multiplying it by $\\Delta r^{2}$ results in the following [tridiagonal system](https://en.wikipedia.org/wiki/Tridiagonal_matrix):\n",
    "\n",
    "$$\n",
    "\\boxed{\\left(1-\\frac{\\Delta r}{r_{i}}\\right)\\psi_{i-1}+\\left(\\pi\\Delta r^{2}\\Phi_{i}^{2}-2\\right)\\psi_{i} + \\left(1+\\frac{\\Delta r}{r_{i}}\\right)\\psi_{i+1} = 0}\\;.\n",
    "$$\n",
    "\n",
    "We solve this system using a cell-centered grid, i.e.,\n",
    "\n",
    "$$\n",
    "r_{i} = \\left(i-\\frac{1}{2}\\right)\\Delta r\\;,\n",
    "$$\n",
    "\n",
    "with $r_{0} = -\\Delta r/2$. We solve the [boundary value problem](https://en.wikipedia.org/wiki/Boundary_value_problem) by imposing the conditions\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\left.\\partial_{r}\\psi\\right|_{r=0} &= 0\\;,\\\\\n",
    "\\lim_{r\\to\\infty}\\psi &= 1\\;.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The first condition (regularity at the origin) states\n",
    "\n",
    "$$\n",
    "\\frac{\\psi_{1}-\\psi_{0}}{\\Delta r} = 0 \\implies \\psi_{1} = \\psi_{0}\\;.\n",
    "$$\n",
    "\n",
    "The second condition (asymptotic flatness) states\n",
    "\n",
    "$$\n",
    "\\psi_{N} = 1 + \\frac{C}{r_{N}}\\;,\n",
    "$$\n",
    "\n",
    "where $N=N_{r}-1$ and $N_{r}$ is the number of discretization points. Differentiating with respect to $r$ we find\n",
    "\n",
    "$$\n",
    "\\partial_{r}\\psi_{N} = -\\frac{C}{r_{N}^{2}} = -\\frac{1}{r_{N}}\\left(\\frac{C}{r_{N}}\\right) = -\\frac{1}{r_{N}}\\left(\\psi_{N} - 1\\right) = \\frac{1-\\psi_{N}}{r_{N}}\\;,\n",
    "$$\n",
    "\n",
    "which can be written as\n",
    "\n",
    "$$\n",
    "\\frac{\\psi_{N+1}-\\psi_{N-1}}{2\\Delta r} = \\frac{1-\\psi_{N}}{r_{N}}\\Rightarrow \\psi_{N+1} = \\psi_{N-1} - \\frac{2\\Delta r}{r_{N}}\\psi_{N} + \\frac{2\\Delta r}{r_{N}}\\;.\n",
    "$$\n",
    "\n",
    "Using the boundary conditions on the boxed equation above we find\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\left(\\pi\\Delta r^{2}\\Phi^{2}_{1} - 1 - \\frac{\\Delta r}{r_{1}}\\right)\\psi_{1} + \\left(1+\\frac{\\Delta r}{r_{1}}\\right)\\psi_{2} = 0\\quad &(i=1)\\;,\\\\\n",
    "\\left(1-\\frac{\\Delta r}{r_{i}}\\right)\\psi_{i-1}+\\left(\\pi\\Delta r^{2}\\Phi_{i}^{2}-2\\right)\\psi_{i} + \\left(1+\\frac{\\Delta r}{r_{i}}\\right)\\psi_{i+1} = 0\\quad &(1<i<N)\\;,\\\\\n",
    "2\\psi_{N-1} + \\left[\\pi\\Delta r^{2}\\Phi^{2}_{N} - 2 - \\frac{2\\Delta r}{r_{N}}\\left(1+\\frac{\\Delta r}{r_{N}}\\right)\\right]\\psi_{N} = - \\frac{2\\Delta r}{r_{N}}\\left(1+\\frac{\\Delta r}{r_{N}}\\right)\\quad &(i=N)\\;,\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "which can be written in matrix form as (using $N \\to N_{r} - 1$)\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "b_{0} & c_{0} \\\\\n",
    "a_{0} & b_{1} & c_{1} \\\\\n",
    "& a_{1} & \\ddots & \\ddots \\\\\n",
    "& & \\ddots & \\ddots & c_{N_{r}-2}\\\\\n",
    "& & & a_{N_{r}-2} & b_{N_{r}-1} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "x_{0}\\\\\n",
    "x_{1}\\\\\n",
    "\\vdots\\\\\n",
    "\\vdots\\\\\n",
    "x_{N_{r}-1}\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "d_{0}\\\\\n",
    "d_{1}\\\\\n",
    "\\vdots\\\\\n",
    "\\vdots\\\\\n",
    "d_{N_{r}-1}\n",
    "\\end{bmatrix}\\;,\n",
    "$$\n",
    "\n",
    "where $\\vec{a}$, $\\vec{b}$, and $\\vec{c}$ are the lower, main, and upper lower diagonals of the tridiagonal system, $\\vec{d}$ is the vector of source terms, and $\\vec{x}=\\vec{\\psi}$. Explicitly, these are\n",
    "\n",
    "$$\n",
    "\\vec{a} =\n",
    "\\begin{bmatrix}\n",
    "1-\\frac{\\Delta r}{r_{1}}\\\\\n",
    "\\vdots\\\\\n",
    "1-\\frac{\\Delta r}{r_{N_{r}-2}}\\\\\n",
    "2\n",
    "\\end{bmatrix}\n",
    "\\;;\\;\n",
    "\\vec{b} = \n",
    "\\begin{bmatrix}\n",
    "\\bigl(\\pi\\Delta r^{2}\\Phi^{2}_{0} - 2\\bigr) & + \\bigl(1 - \\Delta r/r_{0}\\bigr)\\\\\n",
    "\\bigl(\\pi\\Delta r^{2}\\Phi^{2}_{1} - 2\\bigr)\\\\\n",
    "\\vdots\\\\\n",
    "\\bigl(\\pi\\Delta r^{2}\\Phi^{2}_{N_{r}-1} - 2\\bigr) & - \\frac{2\\Delta r}{r_{N_{r}-1}}\\left(1+\\frac{\\Delta r}{r_{N_{r}-1}}\\right)\\\\\n",
    "\\end{bmatrix}\n",
    "\\;;\\;\n",
    "\\vec{c} =\n",
    "\\begin{bmatrix}\n",
    "1+\\frac{\\Delta r}{r_{0}}\\\\\n",
    "1+\\frac{\\Delta r}{r_{1}}\\\\\n",
    "\\vdots\\\\\n",
    "1+\\frac{\\Delta r}{r_{N_{r}-2}}\n",
    "\\end{bmatrix}\n",
    "\\;;\\;\n",
    "\\vec{s} = \n",
    "\\begin{bmatrix}\n",
    "0\\\\\n",
    "\\vdots\\\\\n",
    "0\\\\\n",
    "-\\frac{2\\Delta r}{r_{N_{r}-1}}\\left(1+\\frac{\\Delta r}{r_{N_{r}-1}}\\right)\n",
    "\\end{bmatrix}\n",
    "\\;.\n",
    "$$\n",
    "\n",
    "This tridiagonal system can be solved using e.g., Thomas' algorith, which is explained in detail [here](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm). Our implementation of Thomas' algorithm is described in [Step 2](#thomas_algorithm) of this tutorial notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5edd0f0b",
   "metadata": {},
   "source": [
    "<a id='initialize_nrpy'></a>\n",
    "\n",
    "# Step 1: Initialize core Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initialize_nrpy}$$\n",
    "\n",
    "Let's start by importing all needed modules from Python and defining the quadratic function for error estimates. We also check if the system has GPU support."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0a4632c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Initialize core Python/NRPy+ modules\n",
    "import sympy as sp           # SymPy: a computer algebra system written entirely in Python\n",
    "import os, sys, shutil       # Standard Python modules for multiplatform OS-level functions\n",
    "sys.path.append(os.path.join(\"..\"))\n",
    "import outputC as outC       # NRPy+: C code output functionality\n",
    "import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
    "\n",
    "# Step 1.a: Create all output directories needed by this tutorial\n",
    "#           notebook, removing previous instances if they exist\n",
    "Ccodesdir = os.path.join(\"ScalarFieldID_Ccodes\")\n",
    "if os.path.exists(Ccodesdir):\n",
    "    shutil.rmtree(Ccodesdir)\n",
    "cmd.mkdir(Ccodesdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74c13e5e",
   "metadata": {},
   "source": [
    "<a id='thomas_algorithm'></a>\n",
    "\n",
    "# Step 2: Thomas' algorithm for tridiagonal systems \\[Back to [top](#toc)\\]\n",
    "$$\\label{thomas_algorithm}$$\n",
    "\n",
    "Consider the system\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "b_{0} & c_{0} \\\\\n",
    "a_{0} & b_{1} & c_{1} \\\\\n",
    "& a_{1} & \\ddots & \\ddots \\\\\n",
    "& & \\ddots & \\ddots & c_{N-1}\\\\\n",
    "& & & a_{N-1} & b_{N} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "x_{0}\\\\\n",
    "x_{1}\\\\\n",
    "\\vdots\\\\\n",
    "\\vdots\\\\\n",
    "x_{N}\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "d_{0}\\\\\n",
    "d_{1}\\\\\n",
    "\\vdots\\\\\n",
    "\\vdots\\\\\n",
    "d_{N}\n",
    "\\end{bmatrix}\\;.\n",
    "$$\n",
    "\n",
    "This tridiagonal system can be solved using Thomas' algorithm, described in detail [here](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm). The algorithm goes as follows: for $i=1,\\ldots,n-1$, do\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "w &= \\frac{a_{i-1}}{b_{i-1}}\\;,\\\\\n",
    "b_{i} &= b_{i} - w c_{i-1}\\;,\\\\\n",
    "d_{i} &= d_{i} - w d_{i-1}\\;,\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "followed by the back-substitution ($N=n-1$)\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "x_{N} &= \\frac{d_{N}}{b_{N}}\\;,\\\\\n",
    "x_{i} &= \\frac{d_{i}-c_{i}x_{i+1}}{b_{i}},\\ {\\rm for}\\ i=N-1,N-2,\\ldots,0\\;.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Note that this algorihm modifies both $\\vec{b}$ and $\\vec{d}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bc7ec560",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2: Thomas' algorithm for tridiagonal systems\n",
    "def NRPyCritCol_add_thomas_algorithm_to_C_func_dict():\n",
    "    desc     = \"\"\"(c) 2023 Leo Werneck\n",
    "\n",
    "Tridiagonal system solver using Thomas algorithm\n",
    "Reference: https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm\"\"\"\n",
    "    includes = [\"NRPy_basic_defines.h\"]\n",
    "    name     = \"tridiagonal_solver_thomas_algorithm\"\n",
    "    params   = r\"\"\"\n",
    "      const int n,              // Number of points\n",
    "      const double *restrict a, // Lower diagonal\n",
    "      double *restrict b,       // Main  diagonal\n",
    "      const double *restrict c, // Upper diagonal\n",
    "      double *restrict d,       // Source   vector\n",
    "      double *restrict x        // Solution vector\n",
    "\"\"\"\n",
    "    body     = r\"\"\"  // Step 1: First pass prepares b and d\n",
    "  if( b[0] == 0.0 ) {\n",
    "    fprintf(stderr, \"(%s) Error: first diagonal element is zero\\n\", __func__);\n",
    "    exit(1);\n",
    "  }\n",
    "  for(int i=1;i<n;i++) {\n",
    "    const double w = a[i-1]/b[i-1];\n",
    "    b[i] -= w * c[i-1];\n",
    "    d[i] -= w * d[i-1];\n",
    "  }\n",
    "\n",
    "  // Step 2: Second pass performs the back-substitution\n",
    "  const int N = n-1;\n",
    "  x[N] = d[N] / b[N];\n",
    "  for(int i=N;i>0;i--)\n",
    "    x[i-1] = (d[i-1]-c[i-1]*x[i])/b[i-1];\n",
    "\"\"\"\n",
    "    outC.add_to_Cfunction_dict(includes=includes, desc=desc, name=name,\n",
    "                               params=params, body=body, enableCparameters=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90b02e8f",
   "metadata": {},
   "source": [
    "<a id='initial_scalar_field_profile'></a>\n",
    "\n",
    "# Step 3: Initial scalar field profile \\[Back to [top](#toc)\\]\n",
    "$$\\label{initial_scalar_field_profile}$$\n",
    "\n",
    "We implement the following scalar field profiles:\n",
    "\n",
    "| Profile ID        | Profile                                                           |\n",
    "|:-----------------:|:-----------------------------------------------------------------:|\n",
    "| Gaussian Pulse    | $$\\eta\\exp\\bigl[-(r-r_{0})^{2}/\\sigma^{2}\\bigr]$$                 |\n",
    "| Gaussian Pulse v2 | $$\\eta r^{3}\\exp\\bigl[-(r-r_{0})^{2}/\\sigma^{2}\\bigr]$$           |\n",
    "| Tanh Pulse        | $$\\eta\\Bigl\\{1-\\tanh\\bigl[(r-r_{0})^{2}/\\sigma^{2}\\bigr]\\Bigr\\}$$ |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e68a74b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3: Initial scalar field profile\n",
    "def NRPyCritCol_add_initial_profiles_to_C_func_dict():\n",
    "    # Step 3.a: Implement symbolic expressions for phi(0,r)\n",
    "    eta, r, r0, sigma = sp.symbols(\"eta r r0 sigma\", real=True)\n",
    "    phi_profiles = {}\n",
    "    phi_profiles[\"gaussian_pulse\"   ] = eta * sp.exp( -(r-r0)**2 / sigma**2 )\n",
    "    phi_profiles[\"gaussian_pulse_v2\"] = eta * r**3 * sp.exp( -(r-r0)**2 / sigma**2 )\n",
    "    phi_profiles[\"tanh_pulse\"       ] = eta * ( 1 - sp.tanh((r-r0)**2/sigma**2) )\n",
    "\n",
    "    # Step 3.b: Get maximum length of keys (Optional)\n",
    "    max_length = 0\n",
    "    for key in phi_profiles.keys():\n",
    "        if len(key) > max_length:\n",
    "            max_length = len(key)\n",
    "\n",
    "    # Step 3.c: Python function that generates C functions for the initial \n",
    "    #           scalar field profile and its radial derivative\n",
    "    def generate_phi_or_Phi_function(var, expr, profile):\n",
    "        desc     = \"\"\"(c) 2023 Leo Werneck\n",
    "\n",
    "NRPyCritCol: Initial \"\"\"+var+\" profile - \"+profile\n",
    "        outCparams = \"includebraces=False,outCverbose=False,preindent=1\"\n",
    "        includes   = [\"NRPy_basic_defines.h\"]\n",
    "        c_type     = \"double\"\n",
    "        name       = \"NRPyCritCol_initial_\"+var+\"_\"+profile\n",
    "        params     = \"\"\"\n",
    "      const double r,    /* Radial coordinate    */\n",
    "      const double eta,  /* Pulse amplitude      */\n",
    "      const double r0,   /* Pulse center         */\n",
    "      const double sigma /* Pulse width          */\"\"\"\n",
    "        body = outC.outputC(expr, \"return\", \"returnstring\",\n",
    "                            params=outCparams).replace(\"_ =\",\"\")\n",
    "        outC.add_to_Cfunction_dict(includes=includes, prefunc=prefunc, desc=desc,\n",
    "                                   name=name, params=params, body=body, c_type=c_type,\n",
    "                                   enableCparameters=False)\n",
    "\n",
    "    # Step 3.d: Finally, generate all C code for computing the initial\n",
    "    #           scalar field profile and its radial derivative, including\n",
    "    #           \"driver\" functions that select which profile to use\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    prefunc  = \"\"\n",
    "    c_type   = \"double\"\n",
    "    name     = \"NRPyCritCol_initial_var\"\n",
    "    params   = \"\"\"\n",
    "      const int id_type, /* Type of initial data */\n",
    "      const double r,    /* Radial coordinate    */\n",
    "      const double eta,  /* Pulse amplitude      */\n",
    "      const double r0,   /* Pulse center         */\n",
    "      const double sigma /* Pulse width          */\"\"\"\n",
    "    body     = \"\"\"\n",
    "  switch (id_type) {\n",
    "\"\"\"\n",
    "    counter = 0\n",
    "    defines = \"\"\n",
    "    for profile, phi in phi_profiles.items():\n",
    "        Phi = phi.diff(r)\n",
    "        generate_phi_or_Phi_function(\"phi\", phi, profile)\n",
    "        generate_phi_or_Phi_function(\"phi_radial_deriv\", Phi, profile)\n",
    "        defines += \"#define \"+profile.upper()\n",
    "        for i in range(max_length - len(profile) + 1):\n",
    "            defines += \" \"\n",
    "        defines += str(counter)+\"\\n\"\n",
    "        body    += \"  case \"+profile.upper()+\"\"\":\n",
    "    return NRPyCritCol_initial_var_\"\"\"+profile+\"(r, eta, r0, sigma);\\n    break;\\n\"\n",
    "        counter += 1\n",
    "    body += r\"\"\"  default:\n",
    "    fprintf(stderr, \"Unknown initial data key (%d)\\n\", id_type);\n",
    "    exit(1);\n",
    "  }\n",
    "\"\"\"\n",
    "    outC.outC_NRPy_basic_defines_h_dict[\"ScalarFieldID\"] = defines\n",
    "    for var in [\"phi\", \"phi_radial_deriv\"]:\n",
    "        newname = name.replace(\"var\", var)\n",
    "        newbody = body.replace(\"var\", var)\n",
    "        outC.add_to_Cfunction_dict(includes=includes, name=newname, params=params,\n",
    "                                   c_type=c_type, body=newbody, enableCparameters=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9823903",
   "metadata": {},
   "source": [
    "<a id='solving_the_hamiltonian_constraint'></a>\n",
    "\n",
    "# Step 4: Solving the Hamiltonian constraint \\[Back to [top](#toc)\\]\n",
    "$$\\label{solving_the_hamiltonian_constraint}$$\n",
    "\n",
    "We now solve the Hamiltonian constraint described above using Thomas' algorithm. Recall the the tridiagonal system is characterized by the following:\n",
    "\n",
    "$$\n",
    "\\vec{a} =\n",
    "\\begin{bmatrix}\n",
    "1-\\frac{\\Delta r}{r_{1}}\\\\\n",
    "\\vdots\\\\\n",
    "1-\\frac{\\Delta r}{r_{N_{r}-2}}\\\\\n",
    "2\n",
    "\\end{bmatrix}\n",
    "\\;;\\;\n",
    "\\vec{b} = \n",
    "\\begin{bmatrix}\n",
    "\\bigl(\\pi\\Delta r^{2}\\Phi^{2}_{0} - 2\\bigr) & + \\bigl(1 - \\Delta r/r_{0}\\bigr)\\\\\n",
    "\\bigl(\\pi\\Delta r^{2}\\Phi^{2}_{1} - 2\\bigr)\\\\\n",
    "\\vdots\\\\\n",
    "\\bigl(\\pi\\Delta r^{2}\\Phi^{2}_{N_{r}-1} - 2\\bigr) & - \\frac{2\\Delta r}{r_{N_{r}-1}}\\left(1+\\frac{\\Delta r}{r_{N_{r}-1}}\\right)\\\\\n",
    "\\end{bmatrix}\n",
    "\\;;\\;\n",
    "\\vec{c} =\n",
    "\\begin{bmatrix}\n",
    "1+\\frac{\\Delta r}{r_{0}}\\\\\n",
    "1+\\frac{\\Delta r}{r_{1}}\\\\\n",
    "\\vdots\\\\\n",
    "1+\\frac{\\Delta r}{r_{N_{r}-2}}\n",
    "\\end{bmatrix}\n",
    "\\;;\\;\n",
    "\\vec{d} = \n",
    "\\begin{bmatrix}\n",
    "0\\\\\n",
    "\\vdots\\\\\n",
    "0\\\\\n",
    "-\\frac{2\\Delta r}{r_{N_{r}-1}}\\left(1+\\frac{\\Delta r}{r_{N_{r}-1}}\\right)\n",
    "\\end{bmatrix}\n",
    "\\;.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2a87e0f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 4: Solving the Hamiltonian constraint\n",
    "def NRPyCritCol_add_Hamiltonian_constraint_solver_to_C_func_dict():\n",
    "    desc     = \"\"\"(c) 2023 Leo Werneck\n",
    "\n",
    "NRPyCritCol: This function solves the Hamiltonian\n",
    "             constraint for the conformal factor\"\"\"\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    name     = \"NRPyCritCol_initial_data\"\n",
    "    params   = r\"\"\"\n",
    "      const int Nr,         /* Number of points */\n",
    "      const double eta,     /* Pulse amplitude  */\n",
    "      const double r0,      /* Pulse center     */\n",
    "      const double sigma,   /* Pulse width      */\n",
    "      const double rmax,    /* Maximum radius   */\n",
    "      double *restrict rr,  /* Radial points    */\n",
    "      double *restrict phi, /* Scalar field     */\n",
    "      double *restrict psi  /* Conformal factor */\"\"\"\n",
    "    body     = r\"\"\"\n",
    "  // Step 1: Allocate memory for auxiliary variables\n",
    "  double *a = (double *)malloc(sizeof(double)*(Nr-1));\n",
    "  double *b = (double *)malloc(sizeof(double)*Nr);\n",
    "  double *c = (double *)malloc(sizeof(double)*(Nr-1));\n",
    "  double *d = (double *)malloc(sizeof(double)*Nr);\n",
    "  \n",
    "  // Step 2: Set ID type\n",
    "  const int id_type = GAUSSIAN_PULSE;\n",
    "\n",
    "  // Step 3: Compute step sizes\n",
    "  const double dr  = rmax/Nr;\n",
    "  const double dr2 = dr*dr;\n",
    "\n",
    "  // Step 4: Build tridiagonal system\n",
    "  for(int i=0;i<Nr-1;i++) {\n",
    "    // Step 4.a: Set helper variables\n",
    "    const double r    = (i-0.5)*dr;\n",
    "    const double Phi  = NRPyCritCol_initial_phi_radial_deriv(id_type, r, eta, r0, sigma);\n",
    "    const double tmp  = dr/(r+dr);\n",
    "  \n",
    "    // Step 4.b: Set arrays\n",
    "    rr[i]  = r;\n",
    "    phi[i] = NRPyCritCol_initial_phi(id_type, r, eta, r0, sigma);\n",
    "    a[i]   = 1 - tmp + (i==Nr-2)*(1 + tmp);\n",
    "    b[i]   = M_PI * dr2 * Phi * Phi - 2;\n",
    "    c[i]   = 1 + dr/r;\n",
    "    d[i]   = 0;\n",
    "  }\n",
    "  // Step 5: Adjust remaining entries\n",
    "  const int N       = Nr-1;\n",
    "  const double r_N  = (N-0.5)*dr;\n",
    "  const double tmp0 = dr/r_N;\n",
    "  const double tmp1 = -2 * tmp0 * (1 + tmp0);\n",
    "  const double Phi  = NRPyCritCol_initial_phi_radial_deriv(id_type, r_N, eta, r0, sigma);\n",
    "  rr[N]             = r_N;\n",
    "  phi[N]            = NRPyCritCol_initial_phi(id_type, r_N, eta, r0, sigma);\n",
    "  b[0]             += 1 - dr/rr[0];\n",
    "  b[N]              = M_PI * dr2 * Phi * Phi - 2 + tmp1;\n",
    "  d[N]              = tmp1;\n",
    "  \n",
    "  // Step 6: Solve the system\n",
    "  tridiagonal_solver_thomas_algorithm(Nr, a, b, c, d, psi);\n",
    "\n",
    "  // Step 7: Free memory of auxiliary variables\n",
    "  free(a);\n",
    "  free(b);\n",
    "  free(c);\n",
    "  free(d);\n",
    "\"\"\"\n",
    "    outC.add_to_Cfunction_dict(includes=includes, desc=desc, name=name,\n",
    "                               params=params, body=body, enableCparameters=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f3d820ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "def register_C_functions():\n",
    "    NRPyCritCol_add_thomas_algorithm_to_C_func_dict()\n",
    "    NRPyCritCol_add_initial_profiles_to_C_func_dict()\n",
    "    NRPyCritCol_add_Hamiltonian_constraint_solver_to_C_func_dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bef9d597",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 5: Code validation \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "We now validate our code against the trusted `NRPy+` solver which is implemented in `Python` and described in [this tutorial notebook](Tutorial-ADM_Initial_Data-ScalarField.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b75a47f6",
   "metadata": {},
   "source": [
    "<a id='standalone_main_funcion'></a>\n",
    "\n",
    "## Step 5.a: Generate the `main` function of the standalone C code \\[Back to [top](#toc)\\]\n",
    "$$\\label{standalone_main_funcion}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "75b22edf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 5: Code validation\n",
    "# Step 5.a: Generate the main function of the standalone C code\n",
    "def NRPyCritCol_add_initial_data_standalone_to_C_func_dict():\n",
    "    c_type   = \"int\"\n",
    "    desc     = \"\"\"(c) 2023 Leo Werneck\n",
    "\n",
    "NRPyCritCol: Initial data - standalone test\"\"\"\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    name     = \"main\"\n",
    "    params   = \"int argc, char **argv\"\n",
    "    body     = r\"\"\"\n",
    "  if( argc != 6 ) {\n",
    "    fprintf(stderr, \"Correct usage is: %s <Nr> <eta> <r0> <sigma> <rmax> \\n\", argv[0]);\n",
    "    exit(1);\n",
    "  }\n",
    "  const int Nr       = atoi(argv[1]);\n",
    "  const double eta   = strtod(argv[2], NULL);\n",
    "  const double r0    = strtod(argv[3], NULL);\n",
    "  const double sigma = strtod(argv[4], NULL);\n",
    "  const double rmax  = strtod(argv[5], NULL);\n",
    "\n",
    "  double *rr  = (double *)malloc(sizeof(double)*Nr);\n",
    "  double *phi = (double *)malloc(sizeof(double)*Nr);\n",
    "  double *psi = (double *)malloc(sizeof(double)*Nr);\n",
    "  \n",
    "  NRPyCritCol_initial_data(Nr, eta, r0, sigma, rmax, rr, phi, psi);\n",
    "  \n",
    "  FILE *fp = fopen(\"output.txt\", \"w\");\n",
    "  \n",
    "  for(int i=0;i<Nr;i++)\n",
    "    fprintf(fp, \"%.15e %.15e %.15e\\n\", rr[i], phi[i], psi[i]);\n",
    "  \n",
    "  fclose(fp);\n",
    "\n",
    "  free(rr);\n",
    "  free(phi);\n",
    "  free(psi);\n",
    "  return 0;\n",
    "\"\"\"\n",
    "    outC.add_to_Cfunction_dict(includes=includes, c_type=c_type, desc=desc,\n",
    "                               name=name, params=params, body=body, enableCparameters=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e766f58",
   "metadata": {},
   "source": [
    "<a id='c_code_generation'></a>\n",
    "\n",
    "## Step 5.b: Generate all C codes and compile the project \\[Back to [top](#toc)\\]\n",
    "$$\\label{c_code_generation}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "928a9ec7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(EXEC): Executing `make -j34`...\n",
      "(BENCH): Finished executing in 0.20 seconds.\n",
      "Finished compilation.\n"
     ]
    }
   ],
   "source": [
    "# Step 5.b: Generate all C codes and compile the project\n",
    "# Step 5.b.i: Register all C functions from this tutorial\n",
    "register_C_functions()\n",
    "NRPyCritCol_add_initial_data_standalone_to_C_func_dict()\n",
    "\n",
    "# Step 5.b.ii: Register all C functions and basic defines from outputC\n",
    "outC.outputC_register_C_functions_and_NRPy_basic_defines()\n",
    "\n",
    "# Step 5.b.iii: Generate core NRPy+ header files\n",
    "outC.construct_NRPy_basic_defines_h(Ccodesdir)\n",
    "outC.construct_NRPy_function_prototypes_h(Ccodesdir)\n",
    "\n",
    "# Step 5.b.iv: Generate all C code and compile the project\n",
    "cmd.new_C_compile(Ccodesdir, \"NRPyCritCol\", mkdir_Ccodesrootdir=False,\n",
    "                  uses_free_parameters_h=False, compiler_opt_option=\"fast\") # fastdebug or debug also supported"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e845e3a8",
   "metadata": {},
   "source": [
    "<a id='running_c_code'></a>\n",
    "\n",
    "## Step 5.c: Running the C code \\[Back to [top](#toc)\\]\n",
    "$$\\label{running_c_code}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4af5ee0a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(EXEC): Executing `taskset -c 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 ./NRPyCritCol 320 0.1 0 1 5`...\n",
      "(BENCH): Finished executing in 0.20 seconds.\n"
     ]
    }
   ],
   "source": [
    "# Step 5: Running the C code\n",
    "# Step 5.a: Change to output directory\n",
    "cwd = os.getcwd()\n",
    "os.chdir(Ccodesdir)\n",
    "\n",
    "try:\n",
    "    # Step 5.b: Try running the code\n",
    "    cmd.delete_existing_files(\"output.txt\")\n",
    "    cmd.Execute(\"NRPyCritCol\", \"320 0.1 0 1 5\")\n",
    "    os.chdir(cwd)\n",
    "except:\n",
    "    # Step 5.c: If an error occured, return to the base directory\n",
    "    os.chdir(cwd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca3c2a01",
   "metadata": {},
   "source": [
    "<a id='running_python_code'></a>\n",
    "\n",
    "## Step 5.d: Running the (trusted) Python code \\[Back to [top](#toc)\\]\n",
    "$$\\label{running_python_code}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "068700b1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generated the ADM initial data for the gravitational collapse \n",
      "of a massless scalar field in spherical coordinates.\n",
      "\n",
      "Type of initial condition: Scalar field: \"Gaussian\" Shell\n",
      "                         ADM quantities: Time-symmetric\n",
      "\n",
      "Parameters: amplitude         = 0.1,\n",
      "            center            = 0,\n",
      "            width             = 1,\n",
      "            domain size       = 5,\n",
      "            number of points  = 320,\n",
      "            Initial data file = ScalarFieldID_Ccodes/output_trusted.txt.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Step 5.d: Generating trusted initial data\n",
    "# Step 5.d.i: Import the ScalarField/ScalarField_InitialData.py\n",
    "import ScalarField.ScalarField_InitialData as SFID\n",
    "\n",
    "# Step 5.d.ii: Generate the initial data\n",
    "outfile   = os.path.join(Ccodesdir, \"output_trusted.txt\")\n",
    "ID_family = \"Gaussian_pulse\"\n",
    "Nr        = 320\n",
    "rmax      = 5\n",
    "eta       = 0.1\n",
    "r0        = 0\n",
    "sigma     = 1\n",
    "SFID.ScalarField_InitialData(outfile, ID_family, Nr, rmax, eta, r0, sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30de9467",
   "metadata": {},
   "source": [
    "<a id='relative_error_plot'></a>\n",
    "\n",
    "## Step 5.e: Plotting the relative error between the results of the two codes \\[Back to [top](#toc)\\]\n",
    "$$\\label{relative_error_plot}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "5c2b34c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAF3CAYAAAB5dDWiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAABcSAAAXEgFnn9JSAABpyklEQVR4nO3dd3QU1dsH8O+kF9IrLYXQQ+8dBAwgKlWKQABBimKh2ZBmV4qivgoKCgEEFETp/EQIvYWWEEoCaRBSSEghvd33j7hDlmySTbLJbOD7OWfPyc7MvfPMZMuzd+7cKwkhBIiIiIhIZwyUDoCIiIjoScMEi4iIiEjHmGARERER6RgTLCIiIiIdY4JFREREpGNMsIiIiIh0jAkWERERkY4xwSIiIiLSMSZYRERERDrGBIuIiIhIx5hgEREREekYEywiIiIiHWOCRURERKRjTLCIiIiIdIwJFhEREZGOMcEiIiIi0jEmWEREREQ6xgSLiIiISMeYYBERERHpGBMsIiIiIh1jgkVERESkY0ywiIiIiHSMCRYRERGRjjHBIqqEoKAgjBo1CrVr14aRkREkSUKbNm2UDktvLFmyBJIkoU+fPhWu49dff0XXrl1hbW0NSZIgSRK++eYbncVIVBNFRETI74eIiAilw9EJf39/+ZieBEZKB0CUn5+PHTt2YM+ePThz5gzi4+ORkZEBW1tbNG7cGD179sS4cePQokULpUNVEx4eju7du+Phw4cAAHt7exgbG8PR0VHhyJ4cK1aswLx58wAARkZGcHZ2hiRJsLS0VCym9evXIyIiAn369KlU4vg08vf3h7+/Pzw8PDBp0iSlw6mQy5cv46+//oKtrS3efvttpcMhfSaIFHT69GnRuHFjAUB+GBsbC3t7e2FgYKC2fPjw4SI7O1vpkGXvvvuuACAaNmwo7t69q3Q4emnx4sUCgOjdu3eFyru4uAgA4s033xQ5OTm6Da6CevfuLQCIxYsXKx1KjVPZ14M++PXXXwUA4e7urmgc4eHh8mdjeHi4orHoypEjR+RjehLwEiEpZvfu3ejTpw9CQkLg4OCAzz//HCEhIcjJyUFiYiJycnJw/vx5vPfee7C2tsaff/6JjIwMpcOWBQUFAQCGDBmCunXrKhzNk+f+/fuIi4sDALz66qswNjZWOCIiIu3xEiEpIjQ0FOPHj0d2djaaN2+OgwcPol69emrbGBoaokOHDujQoQPmz5+PV155RaFoNVMle7Vq1VI4kidT0WSa55iIahylm9Do6TRq1CgBQJiZmYmbN29qXa6goKDYspiYGDFv3jzRvHlzYWFhISwsLETz5s3F/PnzRWxsrMZ6Hm9ej42NFW+++abw8PAQpqamwtnZWYwePVpcv369WFl3d3e1S5ePP44cOaLz+G7duiVeffVV4eHhIUxMTOTLE483qV+5ckWMGTNG1K5dW5iZmYmmTZuKZcuWidzcXLnuEydOiCFDhghXV1dhamoqvL29xffff6/x3Kri//bbb8WLL74omjZtKqytrYWZmZnw8vISU6ZMEVevXtVYToiKXRIqekyaHkUvzTx48ECsXbtWvPTSS6JFixbCzs5OmJqaCjc3NzF27Fhx+vTpMveXlpYmVqxYIXr16iUcHByEsbGxqFu3rujVq5dYvny5/D9SXRoq7fH4pZq8vDyxbt068cwzzwgHBwdhYmIi6tSpI0aOHFnsdVJU0cuQOTk5Yvny5aJ9+/bCxsZG42tMkxdffFEAEMOGDSt1u1u3bsnxHzt2TG3dgQMHxLBhw0TdunWFsbGxsLKyEp6enuLZZ58Vy5YtE4mJiWXGIYT667mkx6+//lru49fm9VXWZaczZ86Il19+WX7vW1hYCDc3N9GrVy/x0UcfiTt37sjblnUMmi4bx8fHiwULFog2bdoIa2trYWpqKjw9PcUrr7xS6ntHCCHu3r0rpk2bJurVqydMTExE3bp1xaRJk0RoaKhOLhFmZ2eLn3/+WQwYMEA4OzsLExMT4erqKrp06SKWLl0qwsLCNJa7ePGimDBhgnBzcxOmpqbC1tZWdO3aVXz99dciKyur1H1ev35dvPzyy8LFxUU+F7NmzRKxsbFaXSJMTU0Vn3/+uejSpYuws7MTJiYmol69emL06NHi1KlTFToPVYUJFlW72NhYuX/VlClTKlWXv7+/sLW1ld+UlpaWwtLSUn5uZ2cnjh8/Xqxc0Q+nPXv2CGdnZwFAWFhYCFNTU3mdtbW1uHz5slrZDh06CBcXF2FsbCzv08XFRX6cPHlSp/Ft3rxZ1KpVS47P0tJSY4K1b98+YWZmJgAIGxsbIUmSvG7MmDFCCCF+/vlnYWhoKCRJkr+sVI93331X4zmeOHGivI2RkZGwt7cXRkZG8jJTU1Oxfft2jWUrkmCdPHlSuLi4CEdHR3kfjo6O8vnt0KFDsfoBCENDQznBUi2TJEmsWrWqxH1duHBB1K9fX97ewMBA2Nvbq9Xx9ddfCyGE2Lp1a6n/dxcXFxEVFSXXnZycLPr06aMWn62trdr/Zd68eRrjUiUY7777rujWrZt87u3s7IQkSVolWH/88YcAIExMTEpNhJYsWSIACE9PT7Uke+nSpWqvDwsLC/l1WNKPiZJERUUJFxcX+bVvbGxc7Nxt3bq13Mdf2QRr/fr1av8PU1NTYW1tXWLi5+LiIq83MDAodgzLli1Tq/+ff/5Re/8bGxurvf9NTEzEhg0bNMZ94cIFYWdnJ29rbm4un39ra2uxbdu2SiVYYWFhokWLFmrvFTs7O2FhYSEve+utt4qVW7lypdo5s7Gxkd8TAESrVq3EvXv3NO5z//79au+tWrVqyZ9ZtWvXFr/88kupCdalS5dEvXr11N5TVlZWasfw2WeflftcVBUmWFTttmzZopbcVFRUVJT84dW8eXNx4sQJed2xY8dEkyZNBABhb29frBN60QTGzs5OdO/eXZw/f14IIURubq74559/RO3atQUA0bNnT437L6uzs67iq1WrlujcubMcnxBCbvUr+uVha2srRo8eLSIjI4UQhb/03n//fXn9559/LoyNjcUbb7wh4uLihBCFLUCTJk2SvzA0tSZ+/PHHYtmyZSIoKEhuCcvPzxdXr14V48aNk5ON6OjoYmUr06lZm1/oa9asEYsXLxYBAQHyDRAFBQUiLCxMvPXWW0KSJGFoaCguXrxYrGxUVJScxNWvX19s3bpVpKeny3UEBweLJUuWiE2bNqmV07aT+4gRI+Qv0W+//VauOyYmRrzyyivysf3444/Fyqr2UatWLVGrVi3x66+/ioyMDCGEEAkJCVq1HGVlZclf0Jr2odKwYUMBQCxatEheFhERIf8ImjNnjtr/Njk5WRw/fly89tprIiAgoMw4itL29aDt8VcmwUpPT5e/nMePHy9u3bolr0tLSxMBAQFi/vz5Yu/evWrltO3kHhgYKMzNzQUA8eqrr4pr166JvLw8IYQQkZGR4rXXXpMTx6LvbSEK37tubm4CgHBzcxP/+9//5OT31KlTwtvbWy1xK2+ClZKSIho1aiR//v30008iOTlZXn/79m2xYsUKsXLlSrVyu3fvlvc5ZMgQuYUrOztb+Pn5yeezW7du8rGq3LlzR05OW7VqJc6ePSuEKPws2b9/v6hXr57aMT3u3r178g/h4cOHi4CAAPnGl7i4OLFw4UL5h9/OnTvLdT6qChMsqnYffvih/CbS9KWsrRkzZsgfEDExMcXWF31Dv/7662rrin55N23aVP7wLmrXrl3yNkUvE6iU9UWrq/jc3d3Fw4cPNe6j6JfHs88+q/EyX8+ePeVtpk6dWmx9Xl6e8PT0FADExx9/rHE/pRk8eHCJZas6wSrL66+/LgDNLaXjx48XAISDg4Nay1NZtEmwzpw5I8e+Zs0ajduoEjBHR0eRmZmpcR8AxK5du7SO7XHTp08XAETXrl01rj916pS8n9DQUHm5qnWkcePGFd63JuVNsMo6/sokWGfPnpV/HBS9hF4WbROsvn37CgDi/fffL3GbN998U05Wivryyy/l5PzatWvFysXExKi1bpX3/aH6DDY1NdX446MkzZo1k390Pp5ACaH+mfnHH3+orZs5c6b8flP9wCsqKChIrSXscaofJS+//HKJ8a1cuVIAEK1bt9b6mKoS7yKkapeYmCj/bW9vX6E6hBD4/fffAQAzZsyAq6trsW3q1auHGTNmAAC2bt1aYl1z586Fubl5seWDBg2CiYkJgEd3DCoR36xZs7Tq5P3uu+9qHKBvwIAB8t/vv/9+sfWGhobo168fACAwMLDM/Txu8ODBAIATJ06Uu2xVKym29PR0bNu2DQDw3nvvoX79+jrdr6ruevXqYerUqRq3+fjjjwEACQkJ+OeffzRu4+3tjRdeeKHCcUyYMAEAcPr0ady6davY+o0bNwIAunbtioYNG8rLbW1tAQAPHz5Eenp6hfdfWZU9/tKojlF117IuRURE4PDhwzAyMpLHcdPE19cXAHDo0CHk5+fLy1WfBy+99BKaNWtWrJyrq6v82VERv/zyCwBg6tSpaNu2rVZlAgMDcf36dQDAhx9+CENDw2LbvPDCC+jUqRMAYMuWLfJyIYT8npgxYwacnZ2LlW3RogVGjhypcd9ZWVn47bffABR+zpVEdT6vXLki34GsJCZYVCOFh4fjwYMHAID+/fuXuN2zzz4LoDCpCw8P17hN586dNS43MjKCk5MTAMj7UiK+7t27a7VP1Qfb41xcXAAUJrMNGjQodZukpCSN669cuYLXXnsNrVq1grW1NQwMDOQRl1977TUAwN27d7WKU9fCwsIwb948tG/fHra2tjA0NJRje+655zTGFhAQgNzcXACoki/wgIAAAMAzzzwDAwPNH7PNmjWTh/dQbf84bf/3JenevTu8vLwAAJs2bVJbl5OTI3/pqb6YVDp16gRHR0fExMSgc+fO+P7773Hjxg0IISoVT3lV9vhL4+XlhaZNmyI3NxedO3fGl19+icuXL6slOhV18uRJAEBBQQGaN28OV1dXjY+BAwcCKEz4VUleTk6O/IOub9++Je6jtHWliYyMxL179wCU77Wveo0aGRmhd+/eJW6n+kwr+pou+nlYkWO6cOECsrKyAAA+Pj4lnk9vb2+5TGRkpJZHVnWYYFG1c3BwkP8ub+KiEh8fL/9d2hhURYd+KFqmKCsrqxLLGxkVjmSi+jJWIj5Nv/Y0Kek4VMdQ0eP8/vvv0a5dO/z4448ICgpCWloabGxs4OLiAhcXF1hbWwOAIi0dO3fuRPPmzbFixQpcvHgRKSkpqFWrFpydneHi4gI7OzuNscXGxsp/u7u76zwu1f+yrPHRVP//yv7vS6NqxXo8wdq3bx8ePHgAExMTjB49Wm2dra0ttmzZAicnJwQHB+ONN95As2bNYGdnhxdffBGbNm0q93uiInRx/CUxNDTE1q1b4enpicjISLz33nto27YtrK2t8eyzz+LHH3+s8Lh7qgSmoKAAcXFxJT4SEhLkMqp9PXjwAHl5eQC0/+woj4q+9lWvUUdHR5iampYZV9HXdEU+D4tSnU8ApZ7Poq1W+jBmIhMsqnZFf2VcunRJwUhqBk1N8dXl+vXrePvtt1FQUICXXnoJ586dQ1ZWFpKSkhAbG4vY2FisXLkSAKq9dSMxMRGTJk1CdnY2+vbtC39/f2RkZCAlJQVxcXGIjY3FH3/8obFsTZnrTBf/e1WCdfv2bbllBXh0efD555+XE9Gi+vfvj/DwcPj5+WHixIlo1KgRUlJSsHv3bkyYMAFt27ZFdHR0peMrTVW/9lu3bo0bN25gx44dmDZtGlq0aIHMzEwcOnQIr732Gpo2bVru7gEA5FYwFxcXiMK+zmU+PDw8dHx0mtWU135RRVsVMzMztTqf+jCNFRMsqnZFL5vs3LmzQnUU/WVb2qWpouuq8tfw4/Q9Pm1t374d+fn5aNasGbZu3YqOHTvK/dJUiv4irk779u1Damoq7OzssHv3bvTu3btYX7qSYivaJ64qLiWo/pdlXTZVra/K/32DBg3kS22qpCopKQl79+4FUPzyYFGWlpaYMGEC1q9fj5CQENy9exdffvklzMzM5JYtJalaXlWXjzRJSUkptQ4TExMMHz4ca9asQVBQEO7fv4/Vq1fD3t4ed+7cwcSJE8sdl+r1lZCQUO6WXXt7ezmxLC2BrWhyW9HXvuo1mpCQgOzs7BK30/SaLvp3RY6pqt+vVYUJFlU7FxcXjBgxAgDw22+/ISQkROuyqlYST09PuYP8v//+W+L2hw4dAlB4WdLT07OiIZebvsenrTt37gAo/KVfUl8i1TFUN1VsTZo0gYWFhcZtSoqtQ4cOcqK4e/fucu1XdR5Ka7Hr0KEDAODIkSMoKCjQuM2NGzfkL5SOHTuWK4byUiVRv//+O3JycvD7778jOzsbjo6Ocj81bdStWxfvvPMO5s6dCwAlds4viTbnrjxULW+q14ImZ8+eLVedDg4OmD59Or788ksAha3sRTvBa3MMqoQ2Pz8f+/fvL9f+TUxM0KpVKwCFr5+SHD58uFz1qri5ucmX6crz2le9pvPy8nD06NESt1O954q+pot+HlbkmIr+sCvv+1VJTLBIEZ988glq1aqFzMxMDB8+vMxfY0lJSRgxYoT8a1SSJLnfyJo1azS2VNy7dw9r1qwBAIwdO1bHR1A6fY9PWzY2NgAK76LU9IWyf/9++Pv7V3NUhVSxhYSEaGzBuHz5snzn0eMsLCwwZswYAMAXX3xR6hf041R9zpKTk0vcRlV3dHQ01q5dq3GbRYsWASjs01LajRC6MGrUKJiamiIpKQm7d++WW7LGjBmjcY7H0looAMgthSUl3SXR5tyVR+vWrQEUvpc0JVLx8fH4+eefNZbV9hgB9ePU5hgaNWokX6JasGBBma1oj/dFVX12/PHHH7h582ax7ePj47F69epS6yzNlClTAABr167VuptGq1at0Lx5cwCFn9+abgbYt2+f/H8o+pkmSRJGjRoFAFi9erVa3zOVa9euYfv27Rr3bWlpiZdffhkA8OWXXyIqKqrUWCvat1fnqmUwCCINdu7cKUxMTOSxgL744gu1sXjy8vLExYsXxcKFC+UB6JKSkuT1d+7ckZd7e3urjaB+4sQJecyWsgbyLG0MGdW0OEVHc1Ypazyk6ohPm6kltBm3p6TxhA4dOiTXP3PmTHmAx7S0NLF69WphYWEhHBwcSqy/KsfBCgkJkQfDHD58uHwOs7OzxbZt24STk5Mcm6bzc+fOHbWBRrdt2yaPh1ZQUCCCgoLEvHnzhJ+fn1q5BQsWCACiYcOGxf5vRRUdaPS7775TG2h06tSpclylDTRa1mCm5TFy5EgBQLRv317et2qwx8ctXbpUDBw4UPj5+amNAZeVlSW2bdsmzwIwduzYcsXwzz//CKBwBO6i74fHaXv8+fn58nu0SZMm4vz586KgoEDk5+eLI0eOiGbNmgl7e3uNr4H169eLbt26idWrV4vbt2/Ly/Py8sSBAwfkEcMfH0MsNDRUrm/btm0lxhYUFCSPvN60aVPx119/qY13dvfuXeHn5yf69u1bbHy6lJQUef8eHh7i0KFD8hh3Z86cES1btqzUQKOpqanFBhpNSUmR19+6dUssXbq02Mj0RQcaHTp0qDzQaE5Ojti0aZM8rp+mgUYjIyPlgUjbtGkjD65aUFAgDh48KNzc3MocaLROnToCgKhTp47w8/MTqamp8vr4+Hixfft2MXToUOHj41Ou81FVmGCRok6cOCGPJK16mJiYCHt7e/nLEyicAmHs2LHyyL0q/v7+alO+PD4Vja2tbbH51YSongSrOuKr6gRLCCHGjBmj9v+xtbUVhoaG8pf1d999p0iCJYQQ7777rlpsRaft8PT0FJs3by71/Fy4cEHUrVtX3sbQ0FA4ODjI03cAj6bKUQkJCZHXq6ZLcXd3F+7u7mrJSHJystqAmUWnelEtK2uqHF0mWEUHgVR96Zek6BREQOE0Lfb29mqxN2vWTOMAuqXJzc2VZzBQfbmrzl3RgSnLc/wHDhxQG6DSwsJC/v80atRIbeaIoh6fW9LU1FQ4ODiofe7UqVNH43yk/fr1k7exsrKSj+Hx18qJEyeEq6trsdeXaoR31UPTAMDnz59XSziKTlVkZWVV6alybt++LZo3by7XoZomqrxT5dja2so/lAGIli1bljiA9J49e9SmyrGyspLPhTZT5Vy7dk00bty4WMxFP1MBiP79+5f7fFQFJlikuLy8PLFlyxYxbtw40bBhQ2FtbS2MjY2Fo6Oj6NGjh1iwYIG4ceNGieXv3bsn5s6dK5o1aybMzc2FhYWFaNasmZg3b16JXwDVlWBVdXzVkWDl5+eLb775RrRq1UqYmpoKKysr0aZNG/H555+LrKysUuuvjpHc/fz8RKdOndTO7QcffCCSk5O1njz2iy++EF26dJG/LOrXry/69OkjVq5cqXHU6dOnT4sXX3xRuLi4qM3LWNJkz3369BF2dnbC2NhY1K5dW4wYMULryZ51JTc3Vzg5OcmxfvrppyVuGx0dLX766ScxduxY0aJFC+Hg4CDPQ9mzZ0/xzTffFBt9Xlt3794VU6dOFZ6enmpfzCVN9qyNM2fOiOeff16ei7JRo0bivffeE6mpqSW+BhITE4Wfn5+YPHmyaN26tXB2dhZGRkbCxsZGdOrUSXz88cdqLeZFJSUlidmzZ4vGjRurJeOa4k1NTRXLly+XJxM3NDQUtWrVEs2aNRPjx48XmzdvFmlpaRr3ExUVJaZOnSrq1q0rT/Y8ceJEnU72/MMPP4g+ffrIE53XqVNHdO3aVXz88cciIiJCY7kLFy6I8ePHi/r16wsTExNhY2MjunTpotVkz8HBwWLMmDHC2dlZmJqaCg8Pj3JN9pyVlSXWrFkjfHx85P+ZhYWFaNiwoXjppZfETz/9JB48eFCh86FrkhDVfG81ERER0ROOndyJiIiIdIwJFhEREZGOMcEiIiIi0jEmWEREREQ6xgSLiIiISMeYYBERERHpGBMsIiIiIh1jgkVERESkY0ywiIiIiHTMSImdljUTtq7Z2trKM6ATERERVTVFpsoxNDSs1v0tXrwYixYtqtZ9EhER0dNLkRYsTn/45HF1dUV6ejrc3NyUDoWIiEgrUVFRsLS0RGxsrM7rViTBOnLkSLXuz8PDo1r39zRKT09Hbm6uTusDAEtLS53V+bTguasYnreK47mrGJ63itPVucvNzZXr0jVFLhHSk8fb2xsAEBwcrJP6Dh8+DADo27evTup7mvDcVQzPW8Xx3FUMz1vF6erc6fq7qyjeRUhERESkY0ywiIiIiHRMrxIsOzs72NvbY8WKFUqHQkRERFRhinRyL0l6ejry8/PRuXNnpUMhIiIiqjC9asFycXEBAFhYWCgcCREREVHF6VWC1bp1awBASEiIwpEQERERVZxeJVi+vr4QQmDdunVKh0JERERUYXqVYI0aNQovvPACDh8+jDlz5iA/P1/pkIiIiIjKTa86uR87dgyzZs3CnTt3sGrVKuzatQtjx45FmzZtYG9vX+Ychr169aqmSImIiIhKplcJVp8+fSBJkvw8PDwcn332mVZlJUlCXl5eVYVGREREpDW9SrAATgRNRERENZ9eJVi//vqr0iEQERERVZpeJVgTJ05UOgQiIiKiStOruwiJiIiIngRMsIiIiIh0TK8uEZYkPDwciYmJAAAHBwd4enoqHBERERFRyfS2BevUqVMYOXIk7Ozs0LBhQ3Tu3BmdO3dGw4YNYWdnh9GjR+P06dNKh0lERERUjN4lWEIIzJ49Gz179sTOnTuRkpICIYTaIyUlBdu3b0ePHj0wZ84cpUMmIiIiUqN3CdbcuXOxatUqOZny9PTEyy+/jLlz52Lu3LkYN24cPD095fWrVq3C3LlzlQ5bo/T0dGzcuBFvvPEGOnfuDFNTU0iShCVLlpRYJi4uDuvWrcOwYcNQr149mJiYwNbWFr1798aGDRvKPU7YpEmTIElSmY+oqKhKHi0RERGp6FUfrIsXL2LVqlWQJAlubm748ccfMXDgQI3bHjx4EDNnzkRERARWrVqF8ePHo23bttUccelCQ0Ph6+tbrjJz587F5s2bYWRkhA4dOqBHjx6Ijo7GiRMncOzYMezZswdbt24tc9oglR49epS47ubNmzhz5gzc3d1Rv379csVJREREJdOrBGvNmjUQQsDJyQknT55EnTp1Stx2wIABOHHiBNq1a4f79+9j9erVWLNmTTVGWzYrKytMmTIFHTt2RMeOHbF3714sWrSo1DIODg749NNP8eqrr8LJyUlefv78efTv3x/bt2/HunXrMG3aNK1imDp1KqZOnapx3ejRo3HmzBmMHz9ebYoiIiIiqhy9ukR49OhRSJKE+fPnl5pcqdSpUwfz5s2DEAJHjx6thgjLx8vLC2vXrsX06dPRrl07GBsbl1lm1apV+OCDD9SSKwDo2LEj3nvvPQDAli1bKh1bamoqdu/eDQCYMGFCpesjIiKiR/Qqwbp37x6A0i9rPU61rarsk6x169YAdHOsO3bsQGZmJjp27IgmTZpUuj4iIiJ6RK8SrPz8fADQun9R0W1VZZ9kYWFhAABXV9dK17Vp0yYAwPjx4ytdFxEREanTqwTLxcUFAHDhwgWtywQEBADQTdKhz3Jzc/HDDz8AAIYMGVKpuqKjo+Hv7w8jIyOMGTNGF+ERERFREXrVyb1Hjx6IiIjAsmXLMH78eNSqVavU7dPS0rB8+XJIkoTu3btXU5TKWLhwIa5fvw5PT0/MmDGjUnVt3rwZBQUFGDRoEJydnctV1tvbW+Py27dvw9XVFYcPH65UbCrp6ekAoLP6niY8dxXD81ZxPHcVw/NWcbo6d+np6bC0tNRFSMXoVYI1adIkbNq0CREREejXrx82bdqERo0aadw2JCQEEyZMQHh4OCRJwqRJk3Qez7Bhw3D9+vVylfHz80OnTp10GsfWrVvx1VdfwczMDL/99hssLCwqVZ/q8iA7txMREVUNvUqw+vbtixEjRmDHjh0ICAhA8+bN0adPH3Tr1k2+BBgbG4tTp07B398fBQUFAICRI0eib9++Oo8nPDwcN2/eLFeZjIwMncZw+PBhTJo0CQYGBtiyZQu6dOlSqfoCAwMRFBQEa2trvPjii+UuHxwcrHG5qmVLV/8H1a+Sqvi/Pul47iqG563ieO4qhuet4nR17qqq9QrQswQLADZu3IjMzEzs27cP+fn5OHz4sMYmQNWI5oMHD4afn1+VxHL58uUqqVdb58+fx5AhQ5CTk4N169Zh6NChla5T1Xo1YsQImJubV7o+IiIiKk6vOrkDgJmZGfbs2YP169ejdevWxeYhVD3atm0LPz8/7N69G6ampkqHrXPXrl3DoEGDkJaWhpUrV2Ly5MmVrrOgoEAeQ4uXB4mIiKqO3rVgqfj6+sLX1xf3799HYGAgEhMTARSOdN6qVatiA3E+SSIiIuDj44PExEQsWbIEb7/9tk7q9ff3x927d1G/fn306dNHJ3USERFRcXqbYKk4OTmhX79+SodRbeLj4+Hj44Po6GjMnTsXixcv1qqcr68vzp07h88//xzDhg3TuI3q8uC4ceM4NQ4REVEV0qsEy87ODpIkYcGCBZg7d67S4ejEsGHDEBMTA+DRCOxr167FgQMHAAC1a9fGzp075e2nT5+O0NBQWFhYICEhQePdkY6Ojli+fLnasqioKNy8eRMpKSka48jKysKOHTsA8PIgERFRVdOrBCs9PR35+fno3Lmz0qHozKVLlxAZGam2LDo6GtHR0QAAd3d3tXVJSUkACu9G3LBhg8Y63d3diyVYZdm1axdSU1PRtm1bNG/evFxliYiIqHz0KsFycXHBvXv3Kj3Okz6JiIgo1/b+/v4V2k9Z5UaNGoVRo0ZVqG4iIiIqH726i1A1mXFISIjCkRARERFVnF4lWL6+vhBCYN26dUqHQkRERFRhepVgjRo1Ci+88AIOHz6MOXPmID8/X+mQiIiIiMpNr/pgHTt2DLNmzcKdO3ewatUq7Nq1C2PHjkWbNm1gb28PQ0PDUsv36tWrmiIlIiIiKpleJVh9+vRRG58pPDwcn332mVZlJUlCXl5eVYVGREREpDW9SrCAR3MMEhEREdVUepVg/frrr0qHQERERFRpepVgTZw4UekQiHRGCIHjoQk4eSsBOfkFsDE3xvWYVFyNToWVWeFbLy07D63r2aJLA3tsPBOJuNRsOFmZ4quRrdDOzU7hIyAioorSqwSLqDKycvMRnpCutqxACNyKT0NuvsCLrevAxKj4jbN3HmQgLbvs/nt2FiZwtTErcX18ahZO3U5EdHImJAnYeu4Ooh5klFnv3aRM7A2KkZ+nZOZiyvrzmNqzAfLyBW7EpqJACDSrbQ1jQwNk5xXgekwqnKxM8e7AprAxNy5zH0REVL30KsF6EucipIoTQsh98sqanPrA1Ri8uyMIKZm5JW6z+8o9rJvYAUaGhUlWWnYe3vjtIo7cvK91TMPa1sUXI1rCxFA9Udty7g6W7ApGTn6B1nWVJikjF8sO3lRbdjA4rth2l6OS8dnwlvjnWiy6ezkiND4NR67lAAAOpQRBAtDVywHPt6qjk7iIiEg7epVgPYlzEVLFnIvNg9+1XGT8sw9u9hYY3bE+rkanIDkjF9djU+FYyxTmxoYIT0iHnaUx7jzILLPOoyH30XDBfgxq4YovhrfCpPXncCkquVxx7bwUjZ2XostVppFzLTStbY3kjBw0r2MN7zo2SMvKgyQB2bn5WPG/EDz8rwWtnp057iaVfSwq12JSMfT/TgIA/u/IbfWVd6MAAJvPRmH3lXsY2qYu+jd3gbGhXg1/R0T0RNKrBOtJnIuQyu9c+AP8HJSL/P9uKI16kFGsNSc541FLVdHLe0YGEsyM1cdLe/zy3/6rsdh/NVZtmYWJIQxKaSXLKyhAVm7ZrVNmxgbo29QZD7PyEBqXhl6NHfHRkBbFYiqqe0NHbDwTiU6e9ni+VR3sCbyHQ9fi5OOvY2MGSEBMchZU99gevh6H9BztB+I9GByHg8Fx6NLAHt+NbYekjByk/tfaZ2JkgMYuVqXGSERE5aNXCVbr1q1x7949hISEoF27dkqHQwq48yADMzZdkJOL8qhra461EzugWW1rteUZOXmY9Mt5nIt4oLHc2E5u+HRoCxgYlJxg5eQVYMHOIPx5KRr5BZqDa1nXBj+Ma4f69uX7gdDIxQofDWkhP3++VZ0yL+kJIfDBziBsOXen2Lp6tSS0dDSEu7s7dly8i/sPs+V1Z8IeoOOnh4qVMTEywIQu7nh/UFP5EioREVWcXiVYvr6+2LdvH9atW4cxY8YoHQ4p4FpMKtKyClucDCSgk6c9zoQVJkY25saY0sMT/jfjcfG/S3uWJoZY+HxzeDpaoo2bLUyNirfCWJgYYdv0LriblImpGwJwM+6hvK5XYyd8PMS71OQKKExAlr3UGgsGN9PYId7IwAAu1qZl9hXTFUmSsPTFFsjMycfhG/FoUdcGYffTUdvWDOM9smBnJqFv36YY3LI23tx6qVjn/8fl5BVg3YlwbD0XhbGd3DCzjxccaplWy7EQET2J9CrBGjVqFDZv3ow9e/Zgzpw5WLZsWZnT49CTZYC3K7ZM64zJ605jdGNjfDCuC/YGxSA0Lg3ju7jDycoUs55piO3/tcxM7OaBWqZlv4wlSUJ9ewv8MbMr1h4Lw/20HNS2McPUnp7larGxtTCBrYVJZQ5RZ0yMDPDNmLbFlh8+fFj+u2U9GxyZ1wcFBQIr/rmJDacikZadB0MDCQ6WJpAk4P7DbKga5dJz8rH2RDgO34xHv6bOyM4rwLV7qTA2NEDT2laoZWqEUR3ql7uVjojoaaNXCRbnIiQAaO9uj8+6m8HCWIIkScUulxkYSBjVoX6F6rY2M8Ycnya6CLNGMTCQMH9AU8zu3xh3kzLhaGUqJ6YPs3Lx2uaLOB6aIG8fdj8dYffD1eo4HZYIAPj1ZARm9vHCAG9X1DI1wq34NHT1coBhGa2ARERPE71KsDgXIalYGPPLuioYGRrAw9FSbZmVmTE2TO6EM2GJ+Pl4WJnDVqRl52HZwZtqNx70b+aMz4e3gq2FMe9SJCKCniVYAOciJFKCgYGEbg0d0a2hI46H3se/1+NRIAQkAHXtzJGbL3AmLFGtlauoQ9fjcejTQ6hra44vR7RC94YO1dYfjYhIH+lVgsW5CImU17ORE3o2ciq2/PVnGuLQtTisPnobAZFJGstGJ2di/LqzqG9vjk+HtkSvxsXrISJ6GuhVgsW5CIn0W//mLujf3AUhcQ/x1tbLuP8wGwlp2cW2u/MgE1M2nMfq8e1hZmwIOwsTNK9jraFGIqInk14lWERUMzR2scL+t3oCKJwD8mp0ChLTc/B/R24h8G4KACA3X2DKhgAAgKGBhI2vdEK3ho6KxUxEVJ3YG5WIKsXM2BAdPOwxwNsVu2b1wLZpXWD12NAZ+QUCXx64wT6WRPTUYIJFRDrVuYEDvnu5LR7v437lbgo839+HnZfuKhMYEVE1UizB2rVrF3bt2oX09NJHmC5LXFwc5syZg7lz5+ooMiKqrD5NnPHl8FaobWNWbN3sbVfQZ9kRrD0ehoISph0iIqrpFEuwhg4diuHDhyMyMlLj+tDQUDRo0ABeXl6l1pOQkIBvvvkG33zzTRVESUQVNapjfZx+vx8Ovt2r2Gj7EYkZ+GTvdXywM4hJFhE9kRS9RFhaf4ycnBxEREQgIiKi+gIiIp1r4mqFMx/0wydDWxRbt/X8HXx9KESBqIiIqhb7YBFRlatlaoTxXdyxcUon9Gykfifhd4dv4WBwrEKRERFVDSZYRFRtejZywsYpnRG0xAcNnB5N2fPGlks4Hlr6FD1ERDUJEywiqnZWZsZYM7693DcrJ68AMzddRNj9NIUjIyLSDSZYRKSIRi5WWD+5IyxMDAEUTiI9Y9MFpGTmKhwZEVHlMcEiIsV08LDHylFt5OchcWmY9Os5pGfnKRcUEZEOMMEiIkUNbOGK1/o8Go7lUlQyPtgZxFHfiahGY4JFRIqbP6AJxnV2k5//ffkeNp+NUjAiIqLKYYJFRIqTJAlLX/RGZ097edlHu68h6L+Jo4mIahqjsjepWpMnT4alpWWx5UWn0Onbt2+J5Ss71Q4R6QcjQwN8N7Ytnvv2BBLSspGTX4BX/QLwx4yuqG9voXR4RETloniCFRAQUOI66b/ZYo8ePVpd4RCRgpytzfDt2DYYv/YsCgQQm5qFCevO4u9ZPWBjbqx0eEREWlN8qhxdPIjoydHNyxGfDG0pP49IzEDrpf/D6qO3kZNXoGBkRETaU6wFKzw8XKldE5Gee7mzG+4/zFabp/CL/TcQl5qFxS94KxgZEZF2FEuw3N3dldo1EdUAb/RtiIDIBzgemiAv23QmEq9092SfLCLSe7yLsAqlp6dj48aNeOONN9C5c2eYmppCkiQsWbKkxDJxcXFYt24dhg0bhnr16sHExAS2trbo3bs3NmzYUKFLogUFBVizZg26du0Ka2trmJiYoF69enj55Zdx+fLlih8gURUyMJCwZkJ7zPNpLC/LzRf4Yv8Ndg0gIr2neCf3J1loaCh8fX3LVWbu3LnYvHkzjIyM0KFDB/To0QPR0dE4ceIEjh07hj179mDr1q0wNDTUqj4hBEaOHImdO3fC3NwcPXv2hI2NDa5evYotW7Zg+/bt+Ouvv/Dcc89V5BCJqpSFiRFm9W0EF2szzN8eCADYGxSD3hecMKpDfYWjIyIqGVuwqpCVlRWmTJmC1atX48KFC/joo4/KLOPg4IBPP/0U9+7dw+nTp7F161YcP34cZ86cgbW1NbZv345169ZpHcPu3buxc+dOeHh44NatWzh48CB+//13XLt2DV9++SVyc3Px+uuvV+Ywiarc8Hb10KnIGFmL/w7GrfiHCkZERFQ6JlhVyMvLC2vXrsX06dPRrl07GBuXfZv5qlWr8MEHH8DJyUlteceOHfHee+8BALZs2aJ1DMeOHQMATJ8+HXXq1FFbN3/+fNjY2CAiIgLx8fFa10lU3QwNJKwa0wZ2FoXvoczcfLy++RIS07IVjoyISDMmWDVI69atAQD37t3TuoypqWmJ6yRJgiRJMDQ0hI2NTaXjI6pKtW3Msfyl1vLzm3EP0X/lUXy27zqSM3IUjIyIqDgmWDVIWFgYAMDV1VXrMj4+PgCANWvWFEvMvvrqKyQnJ2P8+PGlJmJE+qJfMxfM6P1oYuikjFz8dCwMY346g9SsXAUjIyJSx07uNURubi5++OEHAMCQIUO0Lte7d2/Mnz8fy5YtQ8OGDdGrVy9YW1vj6tWruHXrFiZNmiTXS1QTvDuwCRwsTfDlgRvIKyi8m/BG7EN0/+IwGjhawsnKDF+NbAV7SxOFIyWipxkTrBpi4cKFuH79Ojw9PTFjxoxylf3qq69Qt25dzJs3DwcPHpSXN2zYEM8++yzMzc21rsvbW/Mgj7dv34arqysOHz5crthKoppjUlf1PU2ehnPnBeCLHqbYEZqL0zH5AICHWXm4cjcFQAo6f/oPBngYoW99I9iaSlrV+TSct6rCc1cxPG8Vp6tzl56ernE+ZF1gglWKYcOG4fr16+Uq4+fnh06dOuk0jq1bt+Krr76CmZkZfvvtN1hYaD/IYnZ2Nnx9fbFjxw4sWLAAkydPhoODAwICAvDmm29i3LhxiI6Oxvz583UaM1FVszeTMLWFMXILBALi1KfQyS0A9oTl4WxMPj7oZAobLZMsIiJdYYJVivDwcNy8ebNcZTIyMnQaw+HDhzFp0iQYGBhgy5Yt6NKlS7nKf/755/j999/x1ltvYenSpfLyZ555Bnv37kXz5s2xZMkSTJ48GY6OjmXWFxwcrHG5qmWrb9++5YqvJKpfJbqq72nytJ27Dl1z8dy3xxGdnFls3f1MgfW3TfHHjG4wMSq9y+nTdt50ieeuYnjeKk5X566qWq8AJlilUnqU8/Pnz2PIkCHIycnBunXrMHTo0HLXsXHjRgDAyJEji61zc3ND586dcfjwYVy4cAEDBgyobMhE1c7Gwhh/zOiKLeei0NnTARGJ6fA7HYGQuDQAwJW7KXhmuT86e9ojPScPc32aoLGLlcJRE9GTjgmWnrp27RoGDRqEtLQ0fP3115g8eXKF6rl79y4AlDgMg2p5UlJSxQIl0gN1bM0x16cJAKBHI0eM6+yGD3ZexZZzUQCA6ORM/HkpGgBwMDgOS1/0llu0rM2M0aWBveaK9ciVO8lIz85Di3o22HbuDs6GJyI9Ox+eTpZoXtsaB4NjIQTwbHMXTOzmoVY2PTsPK/4Xgtz8Aswb0AQ25prH5LublIEvD9xEwsNs1LE1x/vPNYVjLd5hTFQRTLD0UEREBHx8fJCYmIglS5bg7bffrnBdrq6uiIqKQkBAAFq2bKm2Lj8/H5cuXQIAeHh4VCJiIv0iSRIWPd8c58ITcft+erH1i3epX+o2NJDwgqchXvQqezBgJWw8HYGFf2u+PH86LFHt+YlbCbAxN8bQtnUBFE6X9c6OQOwNjAEARD3IwPrJHSFJ6v3SsvPyMXVDAG7EPhoh/25SBn57tQsMDdiHjai8mGDpmfj4ePj4+CA6Ohpz587F4sWLtSrn6+uLc+fO4fPPP8ewYcPk5UOHDsW3336LRYsWoXv37mjcuHDi3Pz8fHzwwQeIiIiAu7s7OnToUCXHQ6QUcxNDrJ/cCd/+G4rkzFz8cy2uxG3zCwT+up2H49H5MD37L8xMDPHtmLZoUbfqB+AtKBD4dN917Lh4Fy3r2iDsfrrG/mTl8fa2y3h722WN646G3Ifn+/u0quds+AN4ffBo2waOlvjJtwMaOteqVHxETwMmWFVs2LBhiIkp/OWoGuhz7dq1OHDgAACgdu3a2Llzp7z99OnTERoaCgsLCyQkJGDSpEnF6nR0dMTy5cvVlkVFReHmzZtISUlRW75o0SIcPHgQN2/eRKtWrdCtWzfY29vj0qVLCAsLg7m5OX755RcYGfGlQE+e+vYWWPbf6O/5BQKzt13Griv30MDREnVszSEgEBCRhOy8wrsQE7MEkJUFAHht80V8OaIVVI037g6WsLUwRvC9VOTlP7pr0dXGDO4OJXeUFUIgND4NSek5cLQyhZfTo+REdUlu95XCz4bjoQlaH9sLrevgwNUY5OYLrctUVlhCOmb9dhF/z+oOUyPtJpwnelrxW7WKXbp0CZGRkWrLoqOjER1d2B/E3d1dbZ2qL1RGRgY2bNigsU53d/diCVZJHBwccP78eaxYsQI7d+7EuXPnkJOTg9q1a2PixIl499130axZs/IeFlGNY2gg4duxbfHN6DYwKHLJKz41C8N+OFWs1SjqQQbG/nxGq7pb1rVBTl4BXmxTB6/18ZIvvx0MjsXiv4MRm5olb/vOwCZ4rU9D7A+KwVvbLiMnr6CkatUYSICBJCFfCLzepyHmDWiCmJSmGL/2LCITM/Dx0BawszDGm1suIye/eJ1mxoV9zrJyS99fe3c7/DCuHXzXncPNuOITat+IfYiZmy7irX6NYGxoADcH7YeNIXqaSEKI6vv5Q08s1TANJQ3jUF68fbnieO7KLz41Cyt3HENOvsClZFOEJxTvt6Wt0R3qw83BApk5+Vh99LY82nxRjZxrITIxQ2MiBACv9fHC2E5uastcrM1gaCAhMzcftUwf/TbOyy9AXoGAmXFhi1J6dh4epKvPzShJhXM5AkBsahYKNMQEAEaGElytzSBJEgoKBGKKbLv66G1sPhtVrIyVmRHmtDGEu7UBX3PlxPdqxenq3On6u6sotmAR0VPP2doMPu6FH4df9u6FRX8H41jIfah+f8amZqFoTuJYyxQmhhIS03Pky4sq2wLulLm/0Pg0+W8jAwnTejVAq3q2+NH/FtrUt8XsZxvD2FDzuF1FkysAMDI0QNGrdZamRrA0Lfmjva6tdjM3GBhIatt+OLg5bsWn4Wz4A7XtHmblYemZPDzrZogtdwMwuGVtuYM90dOsRidYUVHqv6bc3NxK2JKISDvGhgb4fLj6Hbdrj4fhk72Fszq0rm+LHTO6wsjQAP9ci8OrfgGl1mdubAi/KZ3gbGWKkatP4/7DbHmdmbEBNkzuhM4NHAAAA1toP5F7dTM3KTyOrw7cxN+X7yE7Nx8Ps/Pk9f9E5QOIwz/X4uBkZYruDcseuJjoSVajEywPDw+5r4MkScjLyyujBBFR+U3p4QkLEyNEPkjHjF5eMPqvdenZ5i74YVw7nAt/AB9vF+wLikFK5qPPIVMjA4zt5Ib27nYAgP+93QuX7yajoEBAkoCmrtaoo2WLkj4wNTLEwuebY+HzzQEAn+y5hrUnwottN/+PKzgwuxeszfRz2Aui6lCjEywVdiMjoqokSRJe7qy5hfy5lrXxXMvaAIBuXqW32thZmuCZJs46j08p7wxsChMjAxwLCkdMuii8CxPAvZQsfLT7Gpb/dwcn0dOoRidYbm5uxQbLIyKi6mFiZIB3BjZFB5PCYSaCRX2s+CcEALD9wl0808QZg1vVVjJEIsXU6AQrIiJC6RCIiOg/M/t44dCNeFy5kwwAmPvHZbjZW6BlvaofsJVI35Q+vTwREZGWjAwN8M3oNrA2K/ztnpVbgEW7riocFZEymGAREZHOeDpa4vuX28nPL0Ul41Z88QFLiZ50TLCIiEinejV2gncda/n50P87hRANo8ITPcmYYBERkc6NbF9P/jstOw/Pf3sC5yMelFKC6MmiSCd3Pz+/KqnX19e3SuolIqLyGdKmLr46cBOZufkAgJz8AkzfeAG73+ih9WjyRDWZIgnWpEmTdD68giRJTLCIiPSEvaUJfvJtj+/+vYVz/7VcPUjPwXs7AuH3SicOsUNPPMUuEQohdP4gIiL90bORE36f0RVfjng09dDx0AT8dq74pNFETxpFWrCOHDmixG6JiEgBozrUx8HgOBy+EQ8AWLrrGpq6WqG9u73CkRFVHUUSrN69eyuxWyIiUoAkSfhsWEsM/vY4EtNzkJNfgDd+u4SDs3vBivMV0hOKdxESEVGVc7Uxw+oJ7WFsWNj36l5KFr48cEPhqIiqDhMsIiKqFh097DHrmUby801norAvKEbBiIiqjl4nWDk5OTh16hR27NiBjRs3IjU1VemQiIioEmb28UJTVyv5+dzfr+DOgwwFIyKqGnqZYN2/fx/Tp0+Hra0tevbsiVGjRmHSpEm4e/eu2na//PILevXqhRdffFGhSImIqDxMjAzw4/j2sPpvvsLM3Hz8HnBH4aiIdE/vEqybN2+iffv2WLt2LbKyskodgsHHxwenT5/G3r17cfTo0WqOlIiIKsLT0RJv9G0oPz9wNVbBaIiqhl4lWDk5OXjxxRdx9+5dmJqaYs6cOdi9e3eJ29erV0++I3H//v3VFSYREVXSQO/a8t+h8WkY/sNJhN1PUzAiIt3SqwRr3bp1CA0NhYmJCQ4dOoTly5dj8ODBpZbx8fGBEAJnzpyppiiJiKiy3Bws0Kz2owmhL0YlY8K6c0jPzlMwKiLd0asEa+fOnZAkCdOnT0e3bt20KtO6dWsAwK1bt6oyNCIi0rHBLV3VnkcnZ8J78UGs/CcEBQWcnYNqNr1KsAIDAwGgzFarohwdHQEADx5wlnYioppkcndPvNS+XrHl3/4bis2cTodqOL1KsJKSkgAAzs7OCkdCRERVzdLUCMteao0bHw9UG7oBAL7afwNxqVkKRUZUeXqVYFlbF16Pj4nRfuC5qKjCXzn29pzTioioJjIzNsTvM7ri02Et5GUPs/Pw07EwBaMiqhy9SrAaNGgAALh+/brWZQ4cOAAAaNGiRRlbEhGRvrI2M8a4zu5Y+HxzednOS9HIzstXMCqiitOrBKt///4QQuDHH39EQUFBmdvfvHkTGzduhCRJGDBgQDVESEREVWlk+3owNSr8anqQnoOXVp/G4r+v4tC1OIUjIyofvUqwXnvtNZiamiIsLAzz5s0rcYBRALhw4QIGDRqErKwsWFlZYcqUKdUYKRERVQUbc2M81/LRGFmBd1Ow4XQkpvoF4GxYooKREZWPXiVYdevWxRdffAEhBFatWoW2bdvi008/lddv3boVS5cuRb9+/dC5c2dERERAkiSsWrVK7r9FREQ1m29Xd43LF/59Fbn5ZV/dINIHRkoH8Li33noLqampWLp0KQIDAxEUFARJkgBALdkSQkCSJHz22WeYOHGiUuESEZGOtXWzw2+vdsa283fw9+V78vKQuDQs/99NvD+omYLREWlHr1qwVBYuXIhjx45h0KBBMDQ0lOcjVD0kSULfvn1x9OhRvPvuu0qHS0REOtbNyxGrxrRFxBeDMa6zm7x8zdEwHLkRr2BkRNrRuxYslW7dumHv3r1IT0/HxYsXER8fj7y8PDg5OaFt27aws7NTOkQiIqoGCwY3w7nwBwiNL5yr8P0/g3Bobm/UMtXbrzAi/U2wVCwtLdGzZ0+lwyAiIoVYmBjh/8a1w+BvjyM3XyA2NQstFh/EyPb18NWIVjAwkJQOkagYvbxEWBGXLl1SOgQiIqoijV2s8GrPBmrLtl+4ixGrT+FcOKdKI/1T4xOswMBADB06FB07dlQ6FCIiqkKz+jZEWzdbtWWXopIxas1pbOHchaRn9P4SYUmCg4OxZMkS7Ny5s9TxsoiI6MlgYWKEP2d2w+37aXhp9WkkZeTK6z7bdx1WZkYw0nC5MCEtBydvJSA169H2JoYGGN3RDQNbuFZL7FXhzoMMLDt4E4np2WrLDSQJvRo5YWpPT/kufKp+epFgXb16FQcPHkRkZCSMjIzg6emJESNGoE6dOsW2DQsLw4IFC/DHH3/IdxUCQMOGDas7bCIiqmaSJKGhsxW2z+yGRX9fxclbhYOPPszKw6zfytdV5FhoArZN64IOHjVvLtv8AoGZmy/ganSqxvXHQxNgYmSAid08qjcwkimaYGVlZWHy5Mn4/fffi61755138OGHH2LBggUAgIKCAixevBjLli1Dbm6unFi5u7vjww8/xKRJk6ozdCIiUpCXUy1sntoF60+GY8nuaxWqI79AYOTq0xjg7QIXazPMeqYhnK3NdByp7u0LisG8P64gI6f0eRoX7wrGVwduICuvAG/2bYS3+jeqpggJUDjBGjt2LHbt2qXxEl92djYWLVoEe3t7vPLKKxg0aBCOHj0qb1u3bl0sWLAAU6ZMgbGxcXWHrpX09HT8+eefOHfuHM6dO4fLly8jJycHixcvxpIlSzSWiYuLw549e7Bnzx6cP38e8fHxsLCwQOvWrfHKK6/A19e3Qk2+mzdvxg8//IDAwEAIIdCsWTNMmzYNU6dOZRMyEdVYE7p6IDY1GydvJaCghO4ihgYSWta1QfM61jCUJKRl5+GL/TeQV1C4/cHgwnkOL0Qm4c/XusHUyLDa4i+vwzfi8Nrmi2rL+jRxwkDvwkudAsDKf0Jw/2HhZcP0/5Kwrw+FoEcjB7R3r3mtdTWVYgnWkSNH8Pfff0OSJJiYmGD06NHo1KkTjIyMcPXqVWzatAkpKSn49NNPERQUBH9/fwCAg4MDPvzwQ8yYMQOmpqZKha+V0NBQ+Pr6lqvM3LlzsXnzZhgZGaFDhw7o0aMHoqOjceLECRw7dgx79uzB1q1bYWio/QfAzJkzsXr1apiYmKBr166wtLTEqVOnMG3aNJw8eRLr168v55EREekHQwMJ7w1qWu5y1mbG+PCvq8gpMvVO8L1UNFt4AJIkoVU9G/zs2wGOtZT/nvk94A4+23cdD7PykF+gnkTamBvjyxGt4FKk5a2Rcy1M33gBiek5atuO+PE02rrZ4qcJHeBkpfxxPekUS7A2b94MADA3N4e/vz86dOigtv6dd95Bz549cefOHaxZswYAMGDAAGzcuBGOjo7VHm9FqCah7tixIzp27Ii9e/di0aJFpZZxcHDAp59+ildffRVOTk7y8vPnz6N///7Yvn071q1bh2nTpmkVw44dO7B69WrY2dnhn3/+Qfv27QEAMTEx8PHxwYYNGzBgwACMHTu24gdKRFTDjOpYH/2bu+DwjXh8sf86EtIKk5ECAUAIXIpKxqd7r+Pr0W2qPTYhBK5GpyIuNQvJmbl4b0cgCjQ0zo3v4oZXunuqJVcA0MHDHife7Yvjofdx5OZ9tTssL0UlY+A3x7BqTFu0rGeDiIR01LE1Z8JVBRRLsM6fPw9JkjBr1qxiyRUA1K9fHx9//DEmTpwISZLQsmVL7N69G0ZGetEvXyteXl5Yu3at/Px///tfmWVWrVqlcXnHjh3x3nvv4YMPPsCWLVu0TrB+/PFHAMC8efPk5AoAateujZUrV8LHxwdfffUVEywieurYW5pgZPt6eL5VbfiuO4dzEerjae28FI1D1+KwcnQbPNvcpcrj2XQmEr+cCEdieg5SMnNL3M7QQMLKUa0xpE3dErcxNzGEj7crfLxdYWFiiHUnwuV1iek5GL/urPzc3tIEO2Z2g6ejpW4OhAAomGDduXMHAPDss8+WuI2Pj4/89+uvv16jkquq0Lp1awDAvXv3ytjykQsXLgAA+vTpU2xd7969YWBggMuXLyMqKgpubm7FtiEietKZGRti2/QuiE7ORG6+wOubL+JaTOHdeQ+z8/CqXwD+7+V26NfMGWbGuu2flZ4rkJGTh4uRyfjwr6ulbrv0RW/0buwEO0sT2Jhr3/d44fPNMb13A0zfeAGXopKLrX+QnoMlu4LlS62GBhI8HCxhYlTjh8pUlGIZy8OHDwEAzs7OJW5TdF3z5s2rPCZ9FxYWBgBwddV+3Jb09HQA0Dh3o4mJCWrVqoXU1FRcuXKFCRYRPbUkSUI9OwsAwCfDWmDsT2eQnfeof9brv12Ei7Upts/ohvr2FjrZ55mYPKy9mgvDY/9Dbr76NUBJAtztLXAvJQs5eQXo38wZE7q4V3haIGcrM2x4pRMW/XUV12JScedBJjJzH92FeDTkPo6G3JefezlZYvuMbrCzNKnYwZFyCVZ+fj4kSSq1s3bRu9vs7Z/uOx9yc3Pxww8/AACGDBmidTknJyfcu3cPkZGRaNasmdq6Bw8eIDW18FdaZGSk7oIlIqrB2rnZ4cS7fbH2RBjWHA2Tl8elZuNVvwD8OrkjbM1NkJyZgzNhiTA1MkSXBg4wNzZEQlo2jobcVxvU1NjAAF7OlkjNzIORoYRejZ2Qly+w6XouCgRQ8FhytfiF5ni+VR04WZkiPTsPsalZcLe3qPSci9ZmxvhmTFsAQFZuPu48yMDcP64g8G5KsW1v309Hp88O4YvhrTCkTR0YGbI1q7ye7mtuNcjChQtx/fp1eHp6YsaMGVqX69WrF7Zu3Yr169dj4MCBaut++eUX+W9Vi2JZvL29NS6/ffs2XF1dcfjwYa1jK42q5U1X9T1NeO4qhuet4p7Uc9fZFHjQyAh/hObJy27EPkTXzyt3nKaGgJkRkJFXfN2YJsZwzw5H0PlwteVV9RP4ZY8CpKRIiMsoTPIEgOz/GrZy8wXm/nEF+89exctN9aslS1evufT0dFhaVk3fMyZYpRg2bBiuX79erjJ+fn7o1KmTTuPYunUrvvrqK5iZmeG3336DhYX2zdPz5s3D9u3bsW3bNri5uWHWrFmwsLDAjh07sGjRIhgZGSEvLw8GBvx1QkT0uEGexhjkaYy1QTk4FVP6wJ7ays5/lMQU1cLBAP3dqncMLmcLAyzu+uguxAIhsPJCDq49eHR59FBUPuxMc9G7nhEsjDluorYUT7AmT56sVfZY1naSJOHff//VZWgIDw/HzZs3y1UmIyNDpzEcPnwYkyZNgoGBAbZs2YIuXbqUq3z79u3x66+/4tVXX8WyZcuwbNkyed3gwYNhbGyMv/76S2MfLU2Cg4M1Lle1bPXt27dc8ZVE9atEV/U9TXjuKobnreKehnPXvVc+Vh0Kxe8Bd+QhHQDA0sQQ+UIgK/dRQuJmb4EWda1hIEkQAGJTsnApKgmmRoZq/Z4AoLalhAk9GyMzJx/Te3uhlqniX8vo2C0XK/93ExtOP2o3+yM0D+cfmGDvmz1hqQcx6uo1V1WtV4AeJFgBAQGlrlf1wyptOyFElYxGfvnyZZ3XWR7nz5/HkCFDkJOTg3Xr1mHo0KEVqmf8+PF45pln8PvvvyMkJARmZmbo168fBg8ejJ49ewIo+dIfEREBpkaGeGdgU8zzaYKkjBwIFE6qbGdhjAIBJGUUJl2GkgRbC+MSv5P2Bsbg9d8KR2I3MgBmtDLBxD76NZeujbkxlg5pgR6NnPCq36Pv3ojEDHT69BDef64ZXu7kVuk+YU86RRMsTVPkUKFr165h0KBBSEtLw9dff43JkydXqr66deti9uzZassyMzNx+fJlWFlZoV27dpWqn4joaWBgIMHhsdHdDSVoPeL74Fa18SDdGweD49DBKhX1rfS3e0b/Zs54d2BTfHnghrwsPScfH/51FXeSMvD+oGallCbFEqzw8PCyN3pKRUREwMfHB4mJiViyZAnefvvtKtnPL7/8gvT0dLz22mswNzevkn0QEZG6CV09MKGrh97fFCBJEmb28YJvV3f0We4vz28IAL+cCMfLndzg7sDBSUuiWILl7u6u1K71Wnx8PHx8fBAdHY25c+di8eLFWpXz9fXFuXPn8Pnnn2PYsGFq6wICAoqNlv/333/jnXfegaOjI5YuXaqz+ImI6MliaWqEb8e0xdvbLiEutTDJys0X+Gj3Nayd2KFKuug8CRTvg/WkGzZsGGJiYgA8GoF97dq1OHDgAIDCKWt27twpbz99+nSEhobCwsICCQkJmDRpUrE6HR0dsXz5crVlUVFRuHnzJlJSio9n0rFjR3h5eaFZs2awtLTE1atXERwcDAcHB+zfv7/GzO1IRETK6OrlgLMf9Mdfl6Lx9rbLAIB/b8Tjj4C7GNWxvrLB6SkmWFXs0qVLxQbxjI6ORnR0NIDiLXlJSUkACu9G3LBhg8Y63d3diyVYpZk9ezb8/f1x4sQJZGZmws3NDXPmzMG7775b6kj6RERERQ1pUwc7Lt7F8dAEAMBHe66hbzNnrfugPU2YYFWxiIiIcm3v7+9fof2UVm7lypUVqpOIiKgoSZKwbGRr+Hx9FKlZeUjLzsO3/4bioyEtlA5N7+jv7QtERESkd1xtzPBmv0by89/ORuFukm7HgHwSKNKCFRUVVa37s7W1hbW1dbXuk4iI6Ek1oas7fj0ZgejkTOQVCBy4GoupPRsoHZZeUSTB8vT0rNb9LV68GIsWLarWfRIRET2pTI0MMbRtHfzfkdsAgEPX45hgPUaRS4RCiGp9EBERkW71a+Yi/30+IgkpGbkKRqN/FGnBOnLkSLXuz8PDo1r3R0RE9KRrU88WjrVMkJCWg/wCgUPX4zCifT2lw9IbiiRYvXv3VmK3REREpCMGBhL6NXXBtoA7AID/87+FIW3qwMiQ988BvIuQiIiIKmhKT0+o5nwOu58uJ1vEBIuIiIgqqLGLFUa0e3RZcMHOqzhyIx45eQUKRqUfmGARERFRhc3xaYxapo96HE1efx4vfn8Cadl5CkalPCZYREREVGG1bcyxYHAztWU3Yh/iy/03FIpIP9S4qXKioqJw+PBh3Lx5E8nJyQAKBxJt0qQJnnnmmWJz+xEREVHVGtOxPm7Fp2HdiXB52cYzkRjatg7au9srGJlyakyCFRUVhddffx379u0DgGLjW0lSYS+75557Dv/3f/8HNze3ao+RiIjoaSRJEhY+3xzvDmyKF78/gRuxDwEAXx64iW3Tusjf0U+TGpFgxcbGolu3brh37x7c3Nzg4+ODxo0bw9bWFgCQnJyMkJAQ/O9//8PevXtx6dIlBAQEwNXVVdnAiYiIniImRgZY/II3xv58BgBwLvwBjocmoFdjJ4Ujq341IsFavHgx7t27h+XLl+Ptt9+GgYHmrmNCCKxYsQLvvPMOli5dih9//LGaIyUiInq6dfVyQM9GjjgemgAAWP6/m+jZyPGpa8WqEZ3c9+/fj+effx5z5swpMbkCCpso582bh+eeew579uypxgiJiIhIZa5PE/nvwLspOBgcp2A0yqgRCVZ8fDxatmyp9fatWrXC/fv3qzAiIiIiKkmb+rZ4tvmjuQq/2H8dWbn5CkZU/WpEguXs7IygoCCttw8KCoKT09N3vZeIiEhfzPVpLI/yHpGYgVX/hiobUDWrEQnWoEGDsHfvXnzzzTdlbrty5Urs27cPgwcPrvrAiIiISKOmrtaY0sNTfv7zsTBEJWYoGFH1qhGd3JcuXYrdu3dj7ty5+Pbbb0u9izAyMhK1a9fGkiVLFI2ZiIjoaTf72cbYFxSL6ORM5BUIfHc4FMteaq10WNWiRiRYrq6uOH36NGbOnIkDBw7gp59+KnY3gmpcrEGDBuGHH37gEA1EREQKszAxwlv9G+Gd7YEAgD8vRWNarwZo5GKlcGRVr0YkWADg7u6Offv2ISIiQh7JPSUlBQBgY2ODJk2aoG/fvvDw8FA2UCIiIpINb1sXPxy5hYjEDOQXCCzYeRVbp3WBgcGTPWxDjUmwVDw8PPDKK68oHQYRERFpwcjQAAsGN8erfgEAgHMRD7A78B6GtKmrcGRVq0Z0ciciIqKa69nmLhjg/WjYht1XYhSMpnowwSIiIqIqN7bTozmCT99OQE5egYLRVL0nMsGaPHkyjIxq3NVPIiKiJ1ZnTweYGBWmHek5+bgUlaRwRFXriUywgEd3FRIREZHyzE0M0dnTXn4++qczWHciXMGIqtYTm2ARERGRfunVSH2WlU/2XsOVO8nKBFPFasR1tL59+5Zr++vXr1dRJERERFRRI9rXwy8nwxGTkgUAEAJY8FcQ/n69BwyfsGEbakSC5e/vD0mSynXZ7/GBSImIiEhZ9pYmOPbOM/jtbBQW7woGAFyNTsW58Afo6uWgcHS6VSMSLBsbG9SvXx9btmzRavsPP/wQu3btquKoiIiIqLyMDQ0wsZsH/rkWhxO3EgAAJ27dZ4KlhFatWuHSpUvw9vbWanvVHIVERESkn3o1dnyUYIUmYP4AhQPSsRrRyb1NmzZIT0/HrVu3lA6FiIiIdKBHw0cd3gOjU5CckaNgNLpXIxKs3r17o1WrVoiOjtZq+6FDh2LRokVVHBURERFVVFNXKzjWMgFQ2Nn98I14hSPSrRpxiXD48OEYPny41tsPGTIEQ4YMqcKIiIiIqDIMDCT0auyEPy8WNp7M+f0KjAwN0LepM2qZ1oj0pFQ1ogWLiIiInjyv9fGCUZHhGd7ccgldPvsXG05FKBeUjjDBIiIiIkU0dLbClB6easvSsvOweFcwTt1OUCgq3WCCRURERIqZ/WxjTOzqjtb1bdWWL911Dbn5NXdCaCZYREREpBgzY0MsHdICf7/eHTtmdoNqnPCbcQ+x86J2N7fpI73qRdagQYNyl5EkCWZmZrCxsUGTJk3QvXt3jB49GlZWVlUQIREREVWV9u52GNGuHrZfuAsAeGdHIM5HPMDzreugZ0NHGNSg6XT0qgUrIiKi2CMyMrLUZeHh4bh+/TrOnj0LPz8/TJ8+HXXr1sWPP/6o9OEgPT0dGzduxBtvvIHOnTvD1NQUkiRhyZIlJZZJSkrC+++/j/79+8Pd3R0WFhawsLCAt7c33nnnHSQkVOyadHBwMF566SU4OTnB3NwcLVu2xDfffIOCgprb/EpERE+e159piKKz3f1x4S4m/nIOs3+/DCEECgq0nzZPSXrVgtWrVy9IkoSYmBiEhITIyz09PeHi4gIAiIuLQ0REBIDC1qvGjRvD2dkZycnJCAkJQXZ2NtLS0jBr1iwkJibiww8/VOJQAAChoaHw9fUtV5no6Gh88cUXsLe3h7e3N7p27YqHDx8iICAAy5Ytw+bNm3HixAl4enqWXdl/Tp8+jX79+iEzMxOdOnWCh4cHjh07htmzZ+PUqVPYtm0b524kIiK94Oloieda1MbeoBi15X9fvofAuykoEALvthYwM9Lv7y29asHy9/fHokWLkJCQACsrK3z55ZeIi4vD7du3cerUKZw6dQq3b99GfHw8vvrqK1hZWSEhIQFLly7FlStXkJKSgg0bNsDJyQlCCCxduhQ3b95U7HisrKwwZcoUrF69GhcuXMBHH31UZpn69esjICAA9+/fx7Fjx7B161bs3bsXkZGRmDBhAu7du4f58+drHUNubi7GjRuHzMxMrFy5EmfPnsW2bdsQGhqKrl274o8//sCGDRsqc5hEREQ6tfD55mjvbgcPBwu15eEJ6YhMzMDJe/kKRaY9vUqwoqKiMGLECOTk5OD48eOYP38+HB0di23n4OCAefPm4fjx48jOzsbIkSNx9+5dmJiYYMKECfD394elpSUKCgqwevVqBY6kkJeXF9auXYvp06ejXbt2MDY2LrOMjY0N2rdvDwMD9X+NmZkZPvvsMwDA4cOHtY5h586dCA8PR+vWrTF79mx5ea1atfD9998DAFasWKF1fURERFXN1cYMO2Z2g//8Z3B4bm8YG6q3Vp2IzlMoMu3pVYL1zTffIDk5GbNnz0arVq3K3L5ly5aYPXs2Hjx4gK+//lpe3qxZM0yePBlCCPj7+1dhxNVLlaCZmJhoXWbv3r0AgJEjRxZb165dOzRo0ABXr16VL7sSERHpkwZOtbDkRW8YG0qob2+OD55rinkdTJUOq0x6lWDt27cPkiTh2Wef1bqMaltVIqHi4+MDoLBV7EmQm5srd44fPHiw1uWuXLkCoDCZ0kS1PDAwsHIBEhERVZFxnd0R8skgHH+nL6b18oKlsX73vwL0rJP73buFt2Wam5trXUa17eMTQdetWxcAkJaWpqPoqt+UKVOQn5+PpKQkXLhwAdHR0ejevTu++uorretQJZj16tXTuF61PDIysvIBExERVZGadjOWXiVYJiYmyMzMRHBwMNq3b69VmeDgYAAo1r8pL6/w+qytra1OY6xOGzZsQH7+o458ffr0wa+//goHBwet61AlmBYWFhrXW1paAgAePnyoVX3e3t4al9++fRuurq7l6h9WmvT0dADl629GhXjuKobnreJ47iqG563idHXu0tPT5e9BXdOrBKtZs2Y4ffo0vvnmG7z88sswMio9vLy8PHz99deQJAnNmzdXWxceHg4AGjvJa2vYsGG4fv16ucr4+fmhU6dOFd5nUaokMSYmBidPnsT777+Pli1bYvv27RgwYIBO9kFERES6p1cJ1rhx43D69GlcuXIFzz//PH799VfUrl1b47YxMTGYMmUKrly5AkmSMG7cOLX1x44dA1Byi4s2wsPDyz3MQ0ZGRoX3V5LatWtj5MiR6NixI1q2bIlJkybh1q1bWmXdtWrVQlJSUolxqX4FaDvyvarF8HGq89y3b1+t6imL6leJrup7mvDcVQzPW8Xx3FUMz1vF6ercVVXrFaBnCdb06dOxfv16BAQE4J9//oGnpyf69++PTp06wdnZGQAQHx+P8+fP49ChQ8jJyQEAdOzYEdOnT5frycrKwu+//17uDvOPu3z5cqWOR9fc3d3Rs2dP7Nu3D2fPntXqheXm5oakpCTcvXtX452Zqn5v7u7uOo+XiIjoaaVXCZahoSEOHDiAF154AadPn0ZOTg7279+P/fv3F9tWiMKh8rt164Zdu3apjRuVkJCAjz/+GEDhZb4nieqS5/3797XavnXr1rhy5QouXryI5557rtj6ixcvAoBWw2IQERGRdvRqmAYAsLe3x4kTJ/Djjz+iZcuWEEJofLRs2RKrV6/G8ePHYW9vr1ZHvXr1MH36dEyfPr1SfbD0TX5+Pk6cOAGgcBBTbaiGdNi+fXuxdZcuXUJYWBhatGgBDw8PncVJRET0tNO7BAsovBVz+vTpuHLlCmJiYnDw4EFs2bIFW7ZswcGDBxETE4MrV65g2rRpNe62zbJs3boVQUFBxZY/ePAA06ZNQ1hYGFq2bFnsLktfX180bdoUO3fuVFs+bNgweHp64sqVK2qDsaanp+P1118HAMydO7cKjoSIiOjppVeXCDVxcXGpVD8qpQ0bNgwxMYUTVt67dw8AsHbtWhw4cABAYQf2oknRgQMHMHbsWDRo0AAtW7aEhYUFoqOjcfHiRaSlpaFu3boaJ2eOiorCzZs3kZKSorbc2NgYmzZtQv/+/TFnzhxs27YN7u7uOH78OGJiYjBy5EhMnDixKk8BERHRU0fvEyyg8G6+xMREAIXzEHp6eiockfYuXbpUbBDP6OhoeWDUxzuXT506FZaWljh58iROnjyJ5ORk1KpVCy1atMALL7yA119/HTY2NuWKoVu3bjh//jwWL14Mf39/XLlyBV5eXpg/fz7eeuutJ64VkIiISGl6m2CdOnUKK1euxL///ovU1FS1ddbW1vDx8cHbb7+Nrl27KhShdso7x1+PHj3Qo0ePcu+nrDkXvb29NfbDIiIiIt3Tuz5YQgjMnj0bPXv2xM6dO5GSklKsg3tKSgq2b9+OHj16YM6cOUqHTERERKRG71qw5s6di1WrVsnPPT090bVrV3nA0djYWJw6dUoeqX3VqlWQJAkrVqxQJF4iIiKix+lVgnXx4kU5YXJzc8OPP/6IgQMHatz24MGDmDlzJiIiIrBq1SqMHz8ebdu2reaIiYiIiIrTq0uEa9asgRACjo6OOHnyZInJFQAMGDAAJ06cgLOzM4QQWL16dTVGSkRERFQyvUqwjh49CkmSMH/+fNSpU6fM7evUqYN58+ZBCIGjR49WQ4REREREZdOrBEs1TlR57qJTbasqS0RERKQ0vUqw8vPzARTOSagt1baqskRERERK06sEy8XFBQBw4cIFrcsEBAQAAFxdXaskJiIiIqLy0qsEq0ePHhBCYNmyZUhLSytz+7S0NCxfvhySJKF79+7VECERERFR2fQqwZo0aRKAwtHP+/Xrh9DQ0BK3DQkJQb9+/eTxsFRliYiIiJSmV+Ng9e3bFyNGjMCOHTsQEBCA5s2bo0+fPujWrZt8CVA10Ki/vz8KCgoAACNHjkTfvn2VDJ2IiIhIplcJFgBs3LgRmZmZ2LdvH/Lz83H48GEcPny42HZCCADA4MGD4efnV91hEhEREZVIry4RAoCZmRn27NmD9evXo3Xr1sXmIVQ92rZtCz8/P+zevRumpqZKh01EREQk07sWLBVfX1/4+vri/v37CAwMRGJiIgDAwcEBrVq1gpOTk8IREhEREWmmtwmWipOTE/r166d0GERERERa07tLhEREREQ1HRMsIiIiIh1T5BJhVd315+vrWyX1EhEREZWHJFTjHVQjAwMDSJKk0zolSUJeXp5O6yTtWVlZITc3F15eXjqpLz09HQBgaWmpk/qeJjx3FcPzVnE8dxXD81Zxujp3t2/fhrGxMR4+fKiLsNQodomwpOEXKvMg5VhaWsLY2Fhn9cXGxiI2NlZn9T1NeO4qhuet4njuKobnreJ0de6MjY2rLMFV5BLhkSNHlNgtVSFdf0h4e3sDAIKDg3Va79OA565ieN4qjueuYnjeKq4mnDtFEqzevXsrsVsiIiKiasG7CImIiIh0jAkWERERkY4xwSIiIiLSMSZYRERERDqmyDhYRERERE8ytmARERER6RgTLCIiIiIdY4JFREREpGNMsIiIiIh0jAkWERERkY4xwSIiIiLSMSZYpFcyMzOxaNEiNG7cGGZmZqhTpw5eeeUVREdHKx2a3rpw4QK++OILDB8+HPXq1YMkSZAkSemw9F5GRgb++usvTJkyBU2aNIGZmRksLS3RunVrfPTRR0hLS1M6RL22cuVKDB8+HI0aNYKNjQ1MTU3h7u4OX19fBAUFKR1ejZGYmAhnZ2dIkoSGDRsqHY5e69Onj/z5pulx4MABpUNUw3GwSG9kZWXhmWeewZkzZ1C7dm307NkTEREROHfuHJycnHDmzBk0aNBA6TD1ztChQ/H3338XW863dunWrl2LV199FQDQrFkztGjRAqmpqTh16hQePnyIpk2b4ujRo3B2dlY4Uv3k6OiI9PR0tGrVCnXr1gUABAcHIyQkBMbGxvjzzz/x/PPPKxyl/ps0aRL8/PwghICXlxdu3bqldEh6q0+fPjh69ChGjBiBWrVqFVs/d+5ctGzZUoHISiCI9MSCBQsEANG1a1fx8OFDefmKFSsEANG7d2/lgtNjX3zxhVi4cKHYtWuXiImJEaampoJv7bKtX79eTJs2TVy7dk1t+b1790Tbtm0FADF27FiFotN/J06cEJmZmcWW/9///Z8AIFxcXERubq4CkdUchw4dEgDEtGnTBADh5eWldEh6rXfv3gKACA8PVzoUrbAFi/RCTk4OnJ2dkZKSgosXL6Jt27Zq61u3bo3AwEAEBASgffv2CkVZM5iZmSE7O5stWJVw+vRpdOvWDaampkhNTYWJiYnSIdUoDRs2xO3bt3HlyhW0atVK6XD0UmZmJlq2bAlTU1P89ddfaNy4MVuwyqBqwQoPD4eHh4fS4ZSJfbBIL5w8eRIpKSnw8vIqllwBwMiRIwEAu3fvru7Q6CnUunVrAEB2djYSExMVjqbmMTY2BgAmpqVYunQpwsLCsHr1avl80ZPFSOkAiADgypUrAIB27dppXK9aHhgYWG0x0dMrLCwMQGGiYG9vr3A0NcvGjRtx8+ZNNGrUCI0aNVI6HL0UGBiIFStWYPLkyXJfU9LeunXrkJiYCAMDAzRu3BhDhw6Fm5ub0mEVwwSL9EJUVBQAoF69ehrXq5ZHRkZWW0z09Fq1ahUAYODAgTA1NVU4Gv22bNkyBAcHIz09HdevX0dwcDDq1KmDLVu2wNDQUOnw9E5BQQGmTp0KW1tbfPXVV0qHUyN98sknas/nzZuHhQsXYuHChQpFpBkTLNILqlviLSwsNK63tLQEADx8+LDaYqKn0759+7Bu3ToYGxvj448/VjocvXfw4EH8+++/8nN3d3f4+fmxr2QJvvvuO5w/fx6//vorHBwclA6nRunVqxemTp2Kbt26oXbt2rhz5w62b9+OTz75BIsWLYK1tTXeeustpcOUsQ8WEdF/bty4gfHjx0MIgWXLlsl9sahkhw4dghACSUlJOHbsGBo1aoTevXvj008/VTo0vRMVFYUPP/wQvXv3xqRJk5QOp8b56KOPMH78eDRo0ADm5uZo3LgxPvjgA/z1118AgCVLliAzM1PZIItggkV6QTWmSUZGhsb16enpAAArK6tqi4meLtHR0Rg4cCCSkpIwZ84cvfolXBPY2tqiZ8+e2LdvH9q3b4+FCxfi/PnzSoelV15//XXk5ORg9erVSofyRPHx8UGHDh2QnJyMs2fPKh2OjJcISS+oOijevXtX43rVcnd392qLiZ4eDx48gI+PDyIjIzF58mQsX75c6ZBqLGNjY4wePRoXLlzA7t270bFjR6VD0ht79uyBra0tZsyYobY8KysLQGGS36dPHwDA1q1b4erqWt0h1liNGjVCQEAAYmJilA5FxgSL9ILqUszFixc1rlct55g6pGtpaWkYNGgQrl27huHDh+Pnn3/mVEOV5OjoCAC4f/++wpHon+TkZBw9elTjuqysLHmdKuki7SQlJQF41F9XH/ASIemF7t27w8bGBrdv38bly5eLrd++fTsA4IUXXqjmyOhJlp2djSFDhuDcuXMYMGAA73zTEVWS4OXlpXAk+kUIofERHh4OoPB8qZbVhIE09cX9+/dx/PhxACUP9aMEJlikF0xMTDBr1iwAhf0UVH2ugMJJZQMDA9G7d2/emUQ6k5+fj7Fjx+Lw4cPo2bMn/vzzTw6MqaWTJ0/iwIEDKCgoUFuem5uL7777Dhs3boS5uTlGjx6tUIT0pDl16hT++usv5Ofnqy2PiIjAsGHDkJ6ejhdffLHEoX6UwEuEpDc+/PBDHDp0CKdOnUKjRo3Qs2dPREZG4uzZs3BycsIvv/yidIh6ae/evWrDCeTk5AAAunTpIi9buHAhBg8eXO2x6bPvv/8eO3fuBFB4Seu1117TuN3y5cvlS15UKDQ0FJMnT4ajoyPat28PBwcHJCQkICgoCDExMTAzM8P69etRv359pUOlJ0RISAgmT54MV1dXtGvXDra2toiMjMSFCxeQlZUFb29v/Pzzz0qHqYYJFukNMzMzHDlyBJ9//jl+++03/PXXX7C3t8ekSZPw8ccf69UvE31y//59jXfOFF3GvjDFqfpsAJATLU2WLFnCBOsxvXv3xgcffICjR48iMDAQCQkJMDExgYeHB0aOHIk333wTDRs2VDpMeoJ07twZM2fOxNmzZ3H+/HkkJSXB0tISbdq0wUsvvYSZM2fC3Nxc6TDVcLJnIiIiIh1jHywiIiIiHWOCRURERKRjTLCIiIiIdIwJFhEREZGOMcEiIiIi0jEmWEREREQ6xgSLiIiISMeYYBERERHpGBMsIiIiIh1jgkVERESkY0ywiIiIiHSMCRYRERGRjjHBIiIiItIxJlhEREREOsYEi4iIiEjHmGARERER6RgTLCIiIiIdY4JFRKRnhBCwtraGJEkYM2YMAODq1at444030KRJE1hYWECSJKxbt07hSImoJEZKB0BEROpCQkLw8OFDAECrVq2wePFifPbZZ8jLy1Pbrl27dkqER0RaYIJFRKRnLl68KP+9bds2BAYGomnTppg1axbat28PIQQCAgLg7e2tYJREVBomWEREeubSpUvy34GBgXjllVewevVqGBsby8u7du2qRGhEpCUmWEREeqZoC9bgwYOxdu1aSJKkYEREVF6SEEIoHQQRET3i4OCABw8ewMTEBLdu3UL9+vWVDomIyol3ERIR6ZGIiAg8ePAAADBs2DAmV0Q1FBMsIiI9UrT/1eDBgxWMhIgqgwkWEZEeKdr/qlu3bgpGQkSVwQSLiEiPqBIsW1tbeHl5KRwNEVUUEywiIj2iSrA4iChRzcYEi4hIT8TGxiI2NhYA0L59e4WjIaLKYIJFRKQniva/6tChg4KREFFlMcEiItITRRMstmAR1WxMsIiI9AQ7uBM9OZhgERHpCdUYWOzgTlTzcaocIiIiIh1jCxYRERGRjjHBIiIiItIxJlhEREREOsYEi4iIiEjHmGARERER6RgTLCIiIiIdY4JFREREpGNMsIiIiIh0jAkWERERkY4xwSIiIiLSMSZYRERERDrGBIuIiIhIx5hgEREREekYEywiIiIiHWOCRURERKRjTLCIiIiIdIwJFhEREZGOMcEiIiIi0jEmWEREREQ6xgSLiIiISMeYYBERERHpGBMsIiIiIh1jgkVERESkY0ywiIiIiHSMCRYRERGRjjHBIiIiItIxJlhEREREOsYEi4iIiEjHmGARERER6RgTLCIiIiId+3+bkPlfmDUv0wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Step 5.e: Plotting the relative error between the two initial data\n",
    "# Step 5.e.i: Import required Python modules\n",
    "import numpy as np                # NumPy: mathematical library for the Python programming language\n",
    "import matplotlib.pyplot as plt   # Matplotlib: plotting library for the Python programming language\n",
    "from IPython.display import Image # Image: Display images on Jupyter notebooks\n",
    "\n",
    "# Step 5.e.ii: Read in the initial data files\n",
    "rC ,sfC ,psiC  = np.loadtxt(os.path.join(Ccodesdir, \"output.txt\")).T\n",
    "rPy,sfPy,psiPy = np.loadtxt(os.path.join(Ccodesdir, \"output_trusted.txt\")).T\n",
    "\n",
    "# Step 5.e.iii: Generate the plot\n",
    "fig = plt.figure(figsize=(4,2.5))\n",
    "\n",
    "plt.grid()\n",
    "plt.title(\"Conformal factor vs trusted code\")\n",
    "plt.ylim(-13.05,-12.65)\n",
    "plt.xlabel(r\"$r$\",fontsize=14)\n",
    "plt.ylabel(r\"$\\log_{10}{|\\rm Rel.\\ Error|}$\",fontsize=14)\n",
    "plt.plot(rC,np.log10(np.abs(1-psiC/psiPy)))\n",
    "plt.tight_layout()\n",
    "\n",
    "outfig = os.path.join(Ccodesdir,\"validation.png\")\n",
    "plt.savefig(outfig,dpi=150,facecolor='white')\n",
    "plt.close(fig)\n",
    "Image(outfig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e374956",
   "metadata": {},
   "source": [
    "<a id='output_to_pdf'></a>\n",
    "\n",
    "# Step 6: Output this module as $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{output_to_pdf}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-ADM_Initial_Data-ScalarField_Ccode.pdf](Tutorial-ADM_Initial_Data-ScalarField_Ccode.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": 11,
   "id": "a97581f1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ADM_Initial_Data-ScalarField_Ccode.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-ADM_Initial_Data-ScalarField_Ccode.pdf\n"
     ]
    }
   ],
   "source": [
    "# Step 6: Output this module as LaTeX formatted PDF file\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ADM_Initial_Data-ScalarField_Ccode\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}